[gmx-users] problems with gromacs angle restraints

[David Mobley]

David,

Thanks for your response. Sorry about my confusion about the case statments;
I don't use much C.

Anyway, I've now identified the exact problem. It seems that angres computes
angles using essentially the cosine of the angle cos_angle(r_ij,r_kl), while
g_angle (which is what I use to compute the angles to restrain to) computes
cos_angle(r_ij, r_kj). For the angres case, when j=k, this is
cos_angle(r_ij,r_jl). Renumber l as k and this is cos_angle(r_ij,r_jk).

So g_angle computes the cosine of the angle between r_ij and r_kj while the
angle restraints (for the case where the two middle atoms are the same) uses
the cosine of the angle between r_ij and r_jk, which turns out in this case
to give different answers; is the cos_angle routine using a dot product
without taking the magnitude? That is, one can write cos(theta)=|r_ij dot
r_jk|/|r_ij||r_jk|, and if cos_angle were using this to compute the angle,
there would be no problem. But I think (based also on the documentation of
the angle restraints) it's using cos(theta)=r_ij dot r_jk/|r_ij||r_jk|, so
the result is dependent on the sign of r_ij dot r_jk. And r_ij dot r_jk has
a different sign from r_ij dot r_kj.

It's not clear to me that this is exactly a bug, but it's possible it is one
in g_angle. At the very least, g_angle is using a different convention to
compute angles than the angle restraints are, and this should at least be
documented.

If you want to fix it, you could interchange the atom ordering used by
g_angle in the angle restraint case.

In the meantime, a temporary fix for my purposes is to impose an angle
restraint on atoms AB-CB rather than AB-BC as I was doing before. (I've
tried this and it works).

Again, though, if you decide this isn't a bug, you should at least document
it somewhere, as it's not at all obvious from the documentation that g_angle
is not computing the same angles as the angle restraints use.

Let me know if you want me to submit a bugzilla on this.

Thanks,
David


On 11/12/05, David <spoel at xray.bmc.uu.se> wrote:
>
> On Fri, 2005-11-11 at 15:17 -0800, David Mobley wrote:
> > David,
> >
> > I'm not using FEP in these calculations.
> >
> > Looking at these routines they look fine, but I'm wondering where
> > pdihs.phiA and pdihs.phiB are set. My suspicion is that either (a)
> > they are converted from degrees into radians twice (as this conversion
> > happens once in the dopdihs_min routine, and possibly also when
> > they're first set), or (b) they are never set and thus assumed zero.
> Have you checked the tpr file with gmxdump -s topol.tpr | less
> Search for ANGR
> >
> > I haven't been able to find where these are set for the F_ANGRES case.
> > I have to admit I'm a bit of a novice to the code, but it looks like
> > they're set for F_ANGRESZ in convparm.c, but the case statement there
> > for F_ANGRES is empty. In fact, I can't find a case statement anywhere
> > in the code for F_ANGRES which does anything. There are just a bunch
> > of empty case statements. (i.e. "grep -r F_ANGRES *" in the src
> > directory only turns up empty case statements. Is this the problem?
> These are not empty case statements, but rather the control flow of the
> program "fall through", i.e. all three are treated like dihedrals.
> case F_PDIHS:
> case F_ANGRES:
> case F_ANGRESZ:
> new->pdihs.phiA=old[0];
> new->pdihs.cpA =old[1];
> new->pdihs.mult=old[2];
> etc.
>
> There is no conversion here, that means that if you define the angle
> restraints in degrees in the top file it *should* be fine.
>
> >
> > If you can't help me out here, I'll know a little more later, because
> > I'm trying a couple test runs where I vary the preferred angle
> > drastically to see if it makes any difference in the restraining
> > energies I'm getting.
> >
> > Thanks,
> > David
> >
> >
> >
> >
> > On 11/11/05, David <spoel at xray.bmc.uu.se> wrote:
> > > On Fri, 2005-11-11 at 13:23 -0800, David Mobley wrote:
> > > > Dear all,
> > > >
> > > > I'm hoping someone can take a look at the source code for computing
> > > > angle restraint energies, or direct me to where to look in the
> source
> > > > code. I'm particularly interested in making sure that angle
> restraints
> > > > are functioning as advertized in the manual.
> > > >
> > > > Here's the situation:
> > > > I have a system with a small molecule in a binding site in the
> > > > protein. I'm trying to impose distance, dihedral, and angle
> restraints
> > > > to restrain the orientation of the small molecule relative to the
> > > > protein. I'm computing the free energy for turning these on, but
> never
> > > > mind that. I've previously done the calculations in AMBER and am now
> > > > doing them in GROMACS with the AMBER ff. Distance and dihedral
> > > > restraints are working as I expect, but angle restraints are not.
> > > > Particularly, the if I begin a calculation with the small molecule
> > > > having a particular set of angles according to g_angle, then impose
> > > > angle restraints to keep it there, at the end of the restrained
> > > > calculation the angle restraints have forced the ligand AWAY from
> its
> > > > initial orientation, at a considerable energetic cost.
> > > >
> > > > My guess is that either:
> > > > (a) The preferred angle read from the topology file is being
> > > > interpreted as being in radians even though it should be in degrees
> > > > (to be consistent with the dihedral restraint format), or
> > > > (b) the functional form being used for angle restraints is NOT that
> > > > given around p. 64 of the manual.
> > > >
> > > > I haven't been able to find the relevant section of code. Can
> someone
> > > > either direct me to it or check?
> > > >
> > > > FYI, here, I'm using angle restraints where (see p. 64) I'm
> computing
> > > > the angle ABBC (so atom k is identical to atom l) because I want to
> > > > restrain the angle ABC. According to the functional form given in
> the
> > > > manual this should work, but if in fact GROMACS is doing something
> > > > different than what's in the manual it might not work.
> > >
> > > Check routines do_pdihsmin and angres in src/gmxlib/bondfree.c
> > > It might well be that there recently was a bug fix in this bit of the
> > > code in case you are using free energy calculations. Without FEP it
> > > should be fine.
> > >
> > > >
> > > > Thanks,
> > > > David
> > > > _______________________________________________
> > > > gmx-users mailing list
> > > > gmx-users at gromacs.org
> > > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > > Please don't post (un)subscribe requests to the list. Use the
> > > > www interface or send it to gmx-users-request at gromacs.org.
> > > --
> > > David.
> > >
> ________________________________________________________________________
> > > David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> > > Dept. of Cell and Molecular Biology, Uppsala University.
> > > Husargatan 3, Box 596, 75124 Uppsala, Sweden
> > > phone: 46 18 471 4205 fax: 46 18 511 755
> > > spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
> > >
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> > >
> > > _______________________________________________
> > > gmx-users mailing list
> > > gmx-users at gromacs.org
> > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > Please don't post (un)subscribe requests to the list. Use the
> > > www interface or send it to gmx-users-request at gromacs.org.
> > >
> > _______________________________________________
> > gmx-users mailing list
> > gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> --
> David.
> ________________________________________________________________________
> David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596, 75124 Uppsala, Sweden
> phone: 46 18 471 4205 fax: 46 18 511 755
> spoel at xray.bmc.uu.se spoel at gromacs.org http://xray.bmc.uu.se/~spoel
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
>
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20051112/b98f4998/attachment.html>