[gmx-users] Wrong Nose-Hoover integrator
Michel Cuendet
michel.cuendet at epfl.ch
Thu Jan 27 11:32:56 CET 2005
Hi everyone,
The integrator for the Nose-Hoover dynamics (NVT) in gromacs is (as far
as I could deduce from the code) :
xi(t) = xi(t-dt) + Qinv*( T(t-dt) - To)*dt
v(t+dt/2) = F(t)/m*dt + v(t-dt/2)*[1 - xi(t)*dt]
x(t+dt) = x(t) + v(t-dt/2)*dt
This is not a Leap-Frog. It is not time-reversible.
A correct Leap-Frog scheme for the NH system would be (Holian et al.
Phys Rev E 52(3) : 2338, 1995):
xi(t) = xi(t-dt) + Qinv*( T(t-dt/2) - To)*dt
v(t+dt/2) = { F(t)/m*dt + v(t-dt/2)*[1 - xi(t)*dt/2] } / [1 + xi(t)*dt/2]
x(t+dt) = x(t) + v(t-dt/2)*dt
And this seems to have an influence : read the following to see how I
came across this problem.
I added a few lines in the gromacs 3.2.1 code so that it integrates the
work of the NH thermostats on the system, to monitor the heat exchanges
between the thermostats and the different parts of the system:
Wthermo = - Ndf * kB * integral( xi_t * T_t * dt)
To test this I computed the work of one thermostat on a box of water,
for different time constants tau_t. When run without thermostat, this
system shows a total energy drift of about 285 KJ/mol on 140 ps. The
thermostat should compensate for that, independently of tau_t. This is
not what I observed (see table below).
Work Work
tau_t 3.2.1 Right Leap-Frog
-------------------------------
1.5 288 280
0.5 291 275
0.2 107 287
0.15 -38 284
0.10 -539 261
0.05 -3298 299
I think that if gromacs wants to earn the respect its maturity and
reputation implies, it should have correct integrators for all
ensembles. The problem with the Leap-Frog proposed above is that one
needs the temperature at t+dt/2, which is not directly available in
gromacs, and would require extra communication in mpi jobs. I would
suggest computing exactly the temperature at half timesteps for the
integrator, and average to get T(t) for output, etc...
I would even recommend implementing correct integrators for all extended
ensembles: NH, NPT, including constraints. A paper (Martyna et al.
Mol. Phys. 87(5) : 1117, 1996) explicitly provides fortran code models
for each of these ensembles. Their approach, based on the Trotter
expansion of the Liouvillean, additionally allows easy implementation of
multiple timesteps, which would be a nice feature in gromacs !!
Even better for NVT : the Nose-Poincare thermostat (Bond et al. J. Comp.
Phys. 151 : 114, 1999) is a hamiltonian formulation of NH, which allows
true symplectic integration (pseudo-code given). Symplectic is what we
need : time-reversibility, built-in energy conservation, and
conservation of phase space volume and geometry.
Cheers,
Michel
--
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Michel Cuendet, Ph.D. student
Laboratory of Computational Biochemistry and Chemistry
Swiss Federal Institute of Technology in Lausanne (EPFL)
CH-1015 Lausanne
Switzerland Phone : +41 1 693 0324
lcbcpc21.epfl.ch/~mitch +41 1 633 4194
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
More information about the gromacs.org_gmx-users
mailing list