[gmx-users] Wrong Nose-Hoover integrator

[Michel Cuendet]

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
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++