[gmx-users] error in improper dihedral perturbation code

[Chris Oostenbrink]

Hello,

We have found an error in the calculation of the derivative of the 
potential energy with respect to the coupling parameters lambda in cases 
where the zero-energy reference angle of an improper dihedral was 
modified. This could be traced back to the (confusing) fact that the 
zero-energy angle is given in degrees, while the improper dihedral force 
constant is given in kJ/mol/rad2.

In file gmxlib/bondfree.c, subroutine idihs, the contribution of a 
perturbation in the improper dihedral angle to the derivative of the 
hamiltonian with respect to lambda is calculated as

dVdl += 0.5*(kB-kA)*dp2 + (pA-pB)*kk*dp;

where in the second term pA and pB are the zero-energy reference values 
of an improper dihedral angle in states A and B. The value has the unit 
degree. kk is the lambda-weighted force constant, with unit kJ/mol/rad2. 
dp is the difference between the actual value at time t of the improper 
dihedral angle and the lamdba-weighted zero-energy reference value. dp 
has unit rad.

Effect: free energy calculations involving changes in improper dihedral 
reference angles will be incorrect.

Fix: Add a multiplication by DEG2RAD in the second term of the indicated 
line of code. A cleaner solution would probably be to convert the 
reference values to radians immediately after reading them in from the 
topology file and removing the conversion factor in the rest of the code.

Cheers,

Chris de Graaf and Chris Oostenbrink

-- 
========================================================================
Chris Oostenbrink

Computational Medicinal Chemistry and Toxicology
Department of Pharmacochemistry
Vrije Universiteit                                  Tel: + 31 20 5987606
De Boelelaan 1083 P 262                             Fax: + 31 20 5987610
NL-1081 HV Amsterdam
The Netherlands                          e-mail: c.oostenbrink at few.vu.nl
========================================================================