Date: Fri, 22 Nov 2002 08:33:34 -0800
From: "David A. Case"
Subject: Re: Restraint Ambiguities

On Thu, Nov 21, 2002, Neema Salimi wrote:
>
> First, how is the restraint "bond" actually calculated? It is based on a
> flat well potential between the two centers of mass for the grouped
> atoms, or is it based on individual atom positions?

See the IR6 variable, p. 116.

>
> Second, how does makeDIST_RST calculate the r3 (and therefore r4)
> distances from the upper bound in these cases where the distance is
> between groups? In my case, I defined a restraint with an upper bound of
> 5.5 Angstroms as a contact between any of the heavy side chain atoms
> between a PHE and a TYR. makeDIST_RST created an .rst file with r3 =
> 10.76 Angstroms. Similarly, in a manual example (pg 136) a restraint is
> defined between the amide proton and any one of the three ALA C-beta
> protons. The upper bound is given as 5.5 Angstroms and r3 = 6.22
> Angstroms. I guess what I am asking is does sander need r3 to equal
> 10.76 Angstroms to achieve the desired result, or should I edit the .rst
> file to make the distance 5.5 Angstroms. Basically, this question
> relates to the first question. Thanks for your help.
>

If ir6=1 (the default), the program uses an average based on the inverse
6th power of the individual distances. For your PHE-TYR example, there are
7 side chain atoms on phe, and 8 on tyr. The program sums r**-6 for each
of the 56 pairs of atoms, since the "ambiguous constraint" model is really
that each pair might be contributing to an overlapped NOESY peak, where
peak intensities are proportional to r**-6.

To ensure that the calcuated intensity is equal to or greater than the
observed intensity (for an upper bound calculation), the program requires:

sum-over-pairs( r**-6 ) >= robs**-6

where robs**-6 has the role of an observed intensity. If the observed
intensity is such that for a non-overlapped peak the distance (robs) would be
5.5 Ang., then (in a worst-case scenario) all 56 pairs could be contributing
equally to the observed intensity:

56*(r**-6) >= (5.5)**-6

or

r <= 5.5*(56)**1/6 = 5.5*1.955 = 10.76

When there are many pairs of atoms contributing to a peak (as in this example)
the model breaks down, since the liklikhood that all would contribute equally
is negligible, and hence the upper bound of 10.76 is too generous. There are
various things that could be done, including modifying the
correct_upperbound() routine in makeDIST_RST.c

For your case, the "NMR" model may not be so good; you might want to consider
using the -nocorr option to makeDIST_RST (to not make these r3 "corrections"),
and to set ir6=0 (by manually editing the RST file made by makeDIST_RST).
That will give you a constraint between the centers of mass of the two groups
of particles, which may be more what you have in mind.

[Aside: the example cited above from the manual, for a single proton vs.
3 methyl protons, is actually out-of-date: it came from an earlier version of
the program that used a slightly different way of computing r3.]

..hope this helps...dac