[gmx-developers] constraining velocities

[Berk Hess]

Hi,

After speaking with several people at the computational chemistry Gordon 
conference
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 
velocities
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

settle
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.

Berk.