[gmx-users] error in improper dihedral perturbation code
Chris Oostenbrink
c.oostenbrink at few.vu.nl
Wed Oct 5 16:33:56 CEST 2005
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
========================================================================
More information about the gromacs.org_gmx-users
mailing list