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