[gmx-developers] Bug in free energy calculation involving position restraints?
Berk Hess
hess at cbr.su.se
Tue Sep 14 11:07:19 CEST 2010
Hi,
This is indeed a bug.
I don't understand how this went unnoticed. I implemented and tested
this code, but apparently I did not
test it well enough.
Your fix produces the correct potential and forces for your particular
case, but not the correct virial.
I committed a general solution for 4.5.2 and 4.0.8 (if this gets released).
If you need a proper fix now, replace the loop that calls harmonic by
the code below.
Berk
for (m=0; (m<DIM); m++)
{
real kk;
kk = (1 - lambda)*pr->posres.fcA[m] +
lambda*pr->posres.fcB[m];
fm = -kk*dx[m];
f[ai][m] += fm;
vtot += 0.5*kk*dx[m]*dx[m];
*dvdlambda +=
0.5*(pr->posres.fcB[m] - pr->posres.fcA[m])*dx[m]*dx[m]
-fm*dpdl[m];
/* Here we correct for the pbc_dx which included rdist */
vir_diag[m] -= 0.5*(dx[m] + rdist[m])*fm;
}
On 09/13/2010 10:52 PM, JR Schmidt wrote:
> I believe I found a bug involving free energy calculations with
> position restraints. It appears that the "interpolation" of position
> restraints with lambda are not working correct (or at least in any way
> that makes sense to me). I did a test case of a single molecule, with
> a structure "A" and "B": setting lambda=0 and minimizing yields
> structure "A", as expected, but setting lambda=1 does not yield
> anything close to structure "B".
>
> After checking the obvious, I went to the code. In bondfree.c,
> posres() subroutine I see the culprit. In the present simple case of
> no refcoord_scaling, the algorithm is:
> 1) Calculate a "ref", which is in this case 0 (the origin)
> 2) Calculate an expected position or the constrained atom, stored in
> "rdist", which is a linear interpolation between posA and posB.
> 3) Calculate "dx", which is the difference between the "rdist" and
> the current coordinates of the molecule.
>
> Makes sense, obviously: when dx = 0, the positions coincide with the
> expected interpolated positions.
>
> But here's the rub. Later on, the code calls the "harmonic" function
> to evaluate the energy and force due to the constraint. The
> "harmonic" function takes the force constants and coordinates of
> structures A/B, as well as the current coordinates as parameters,
> along with lambda. But INSTEAD, the posres() subroutine sends all the
> positions relative to strucutre A, EXCEPT for the current position,
> for which it sends "dx". This is incorrect, since "dx" is relative to
> the expected interpolated position, rather than that of structure A!
>
> Changing the line :
> rdist[m] = (1 - lambda)*posA + lambda*posB;
> to
> rdist[m] = posA;
> seems to yield the expected results (and uses the coordinates of
> structure A as a consistent reference). After making this change,
> setting lambda = 0 and minimizing yields structure A, setting lambda=1
> yields structure B, and lambda=0.5 yields a linear interpolation
> halfway in between.
>
> Am I misunderstanding the way these restraints are supposed to work?
>
More information about the gromacs.org_gmx-developers
mailing list