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

[Berk Hess]

Hi,

This is only with position restraints where the reference position (not
the force constant)
changes as a function of lambda. And then only when lambda is non-zero.
What happens is that, in all cases I think, the reference position is
rA+2*lambda*(rB-rA).
Note the erroneous factor of 2 (which should be 1).
I would think that in most cases the error would be pretty obvious,
since atoms don't
go where you expect them to go. That's why I am surprised this went
unnoticed for so long.
But apparently nobody has used it up till now.
I did test it after implementing it, so I am surprised that I myself did
not notice it,
although I can't exclude that the initial code was correct and the
mistake came in later.

The usual conclusion is that we need more test sets.

Berk

On 09/14/2010 04:36 PM, David Mobley wrote:
> 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
> <mailto: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 <mailto: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
>     <mailto:gmx-developers-request at gromacs.org>.
>
>
>
>
> -- 
> David Mobley
> dmobley at gmail.com <mailto: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/d3fab861/attachment.html>