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