Can I get some clarification on what exactly this affects? Is it all free energy calculations in 4.x using position restraints, or is it only if position restraints vary as a function of lambda?<br><br>Thanks.<br><br><br><div class="gmail_quote">
On Tue, Sep 14, 2010 at 4:07 AM, Berk Hess <span dir="ltr"><<a href="mailto:hess@cbr.su.se">hess@cbr.su.se</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Hi,<br>
<br>
This is indeed a bug.<br>
I don't understand how this went unnoticed. I implemented and tested<br>
this code, but apparently I did not<br>
test it well enough.<br>
Your fix produces the correct potential and forces for your particular<br>
case, but not the correct virial.<br>
I committed a general solution for 4.5.2 and 4.0.8 (if this gets released).<br>
If you need a proper fix now, replace the loop that calls harmonic by<br>
the code below.<br>
<br>
Berk<br>
<br>
<br>
for (m=0; (m<DIM); m++)<br>
{<br>
real kk;<br>
kk = (1 - lambda)*pr->posres.fcA[m] +<br>
lambda*pr->posres.fcB[m];<br>
fm = -kk*dx[m];<br>
f[ai][m] += fm;<br>
vtot += 0.5*kk*dx[m]*dx[m];<br>
*dvdlambda +=<br>
0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]<br>
-fm*dpdl[m];<br>
<br>
/* Here we correct for the pbc_dx which included rdist */<br>
vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm;<br>
<div><div></div><div class="h5"> }<br>
<br>
On 09/13/2010 10:52 PM, JR Schmidt wrote:<br>
> I believe I found a bug involving free energy calculations with<br>
> position restraints. It appears that the "interpolation" of position<br>
> restraints with lambda are not working correct (or at least in any way<br>
> that makes sense to me). I did a test case of a single molecule, with<br>
> a structure "A" and "B": setting lambda=0 and minimizing yields<br>
> structure "A", as expected, but setting lambda=1 does not yield<br>
> anything close to structure "B".<br>
><br>
> After checking the obvious, I went to the code. In bondfree.c,<br>
> posres() subroutine I see the culprit. In the present simple case of<br>
> no refcoord_scaling, the algorithm is:<br>
> 1) Calculate a "ref", which is in this case 0 (the origin)<br>
> 2) Calculate an expected position or the constrained atom, stored in<br>
> "rdist", which is a linear interpolation between posA and posB.<br>
> 3) Calculate "dx", which is the difference between the "rdist" and<br>
> the current coordinates of the molecule.<br>
><br>
> Makes sense, obviously: when dx = 0, the positions coincide with the<br>
> expected interpolated positions.<br>
><br>
> But here's the rub. Later on, the code calls the "harmonic" function<br>
> to evaluate the energy and force due to the constraint. The<br>
> "harmonic" function takes the force constants and coordinates of<br>
> structures A/B, as well as the current coordinates as parameters,<br>
> along with lambda. But INSTEAD, the posres() subroutine sends all the<br>
> positions relative to strucutre A, EXCEPT for the current position,<br>
> for which it sends "dx". This is incorrect, since "dx" is relative to<br>
> the expected interpolated position, rather than that of structure A!<br>
><br>
> Changing the line :<br>
> rdist[m] = (1 - lambda)*posA + lambda*posB;<br>
> to<br>
> rdist[m] = posA;<br>
> seems to yield the expected results (and uses the coordinates of<br>
> structure A as a consistent reference). After making this change,<br>
> setting lambda = 0 and minimizing yields structure A, setting lambda=1<br>
> yields structure B, and lambda=0.5 yields a linear interpolation<br>
> halfway in between.<br>
><br>
> Am I misunderstanding the way these restraints are supposed to work?<br>
><br>
<br>
--<br>
</div></div><div><div></div><div class="h5">gmx-developers mailing list<br>
<a href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
<a href="http://lists.gromacs.org/mailman/listinfo/gmx-developers" target="_blank">http://lists.gromacs.org/mailman/listinfo/gmx-developers</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-developers-request@gromacs.org">gmx-developers-request@gromacs.org</a>.<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>David Mobley<br><a href="mailto:dmobley@gmail.com" target="_blank">dmobley@gmail.com</a><br>504-383-3662<br><br><br>