[gmx-users] problems with gromacs angle restraints

[David Mobley]

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.

Thanks,
David