<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#ffffff" text="#000000">
Hi,<br>
<br>
This is only with position restraints where the reference position (not
the force constant)<br>
changes as a function of lambda. And then only when lambda is non-zero.<br>
What happens is that, in all cases I think, the reference position is
rA+2*lambda*(rB-rA).<br>
Note the erroneous factor of 2 (which should be 1).<br>
I would think that in most cases the error would be pretty obvious,
since atoms don't<br>
go where you expect them to go. That's why I am surprised this went
unnoticed for so long.<br>
But apparently nobody has used it up till now.<br>
I did test it after implementing it, so I am surprised that I myself
did not notice it,<br>
although I can't exclude that the initial code was correct and the
mistake came in later.<br>
<br>
The usual conclusion is that we need more test sets.<br>
<br>
Berk<br>
<br>
On 09/14/2010 04:36 PM, David Mobley wrote:
<blockquote
cite="mid:AANLkTimeuhMPdj9dXLaGQW6bVuZcb4S2cmEymo9ju_wC@mail.gmail.com"
type="cite">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 moz-do-not-send="true" href="mailto:hess@cbr.su.se">hess@cbr.su.se</a>></span>
wrote:<br>
<blockquote class="gmail_quote"
style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; 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 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 class="h5">gmx-developers mailing list<br>
<a moz-do-not-send="true" href="mailto:gmx-developers@gromacs.org">gmx-developers@gromacs.org</a><br>
<a moz-do-not-send="true"
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 moz-do-not-send="true"
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 moz-do-not-send="true" href="mailto:dmobley@gmail.com"
target="_blank">dmobley@gmail.com</a><br>
504-383-3662<br>
<br>
<br>
</blockquote>
<br>
</body>
</html>