[gmx-users] problems with gromacs angle restraints
dmobley at gmail.com
Sat Nov 12 20:53:32 CET 2005
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
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.
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.cpA =old;
> 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
> > > > code. I'm particularly interested in making sure that angle
> > > > 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
> > > > to restrain the orientation of the small molecule relative to the
> > > > protein. I'm computing the free energy for turning these on, but
> > > > 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
> > > > 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
> > > > either direct me to it or check?
> > > >
> > > > FYI, here, I'm using angle restraints where (see p. 64) I'm
> > > > 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
> > > > 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 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
> 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...
More information about the gromacs.org_gmx-users