[gmx-developers] Bug in free energy calculation involving position restraints?

[Berk Hess]

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?
>