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

[David Mobley]

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?

Thanks.


On Tue, Sep 14, 2010 at 4:07 AM, Berk Hess <hess at cbr.su.se> wrote:

> 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?
> >
>
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.
>



-- 
David Mobley
dmobley at gmail.com
504-383-3662
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20100914/44e11498/attachment.html>