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

[JR Schmidt]

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