[gmx-developers] constraining velocities

Berk Hess hessb at mpip-mainz.mpg.de
Mon Oct 16 15:11:04 CEST 2006


After speaking with several people at the computational chemistry Gordon 
I realized that Gromacs loses a lot of accuracy when using constraints 
by doing:
v_(n+1/2) = (x_(n+1) - x_(n))/dt
I have now fixed this (in the development branch) by correcting the 
using the Lagrange multipliers of the constraints.
The error still increases with decreasing timestep, but it has decreased 
a lot.

I have done some NVE tests with
single precision v_(n+1/2) = (x_(n+1) - x_(n))/dt (first column)
single precision with the new velocity correction (second column)
double precision (last column)

The results are the drift in the total energy in kJ/mol /ns /dof

tetrapeptide in vacuum, flexible, nrdf: 123

dt=0.001  0.04 <0.01 <0.01

tetrapeptide in vacuum, bond-constraints, nrdf: 82

shake tol=1e-5
dt=0.002  4.64 -0.18 -0.29

lincs iter=2, order=8
dt=0.001  2.73  0.27 <0.01
dt=0.002  0.47  0.13 <0.01

216 SPC, nrdf: 1293

dt=0.001 -0.14 -0.04
dt=0.002 -0.06 -0.02
dt=0.004  0.01  0.01

Especially shake seems to profit from this change.
Velocity verlet seems to be more accurate than the old Gromacs code,
but less accurate than my new code.


More information about the gromacs.org_gmx-developers mailing list