[gmx-developers] Bug in free energy calculation involving position restraints?
JR Schmidt
schmidt at chem.wisc.edu
Mon Sep 13 22:52:35 CEST 2010
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?
--
J.R. Schmidt
Assistant Professor of Chemistry
Room 8305D
Department of Chemistry
University of Wisconsin-Madison
1101 University Ave
Madison, WI 53706
Phone: (608) 262-2996
Fax: (608) 262-9918
E-mail: schmidt at chem.wisc.edu
http://www.chem.wisc.edu/users/schmidt
More information about the gromacs.org_gmx-developers
mailing list