[gmx-developers] Drift in Conserved-Energy with Nose-Hoover thermostat

[Bernhard]

Dear Gromacs Users and Developers,

I have a problem regarding energy conservation in my 10ns NVT 
protein+water+ions (22765 atoms) production (minimization and 
equilibration for more than 15ns was carried out in NPT before) 
simulations using a Nose-Hoover thermostat (tau=2.5ps).
On first glance everything looks fine - the potential, kinetic and total 
energy are nearly perfectly constant (with normal fluctuations) - but 
when I checked the "Conserved-Energy" quantity that g_energy outputs I 
had to recognize a significant (nearly perfectly) linear downward drift 
of this "to-be-conserved" quantity of around 1.7 kJ/mol/ps (7.48*10^-5 
kJ/mol/ps per atom).
This appears somehow disturbing to me since I would expect that this 
Conserved-Energy is the conserved energy of the extended Nose-Hoover 
Hamiltonian - which should by definition be conserved.

If it would be a drift caused by normal round-off error due to single 
precision I would expect it to grow with Sqrt(N) and not with N (linear) 
(N=number of steps).
So I would like to know if this is a normal behaviour and also what 
could cause this (buffer size, precision, constraints etc)?
Also I would like to know, if I am correct with my guess that the 
"Conserved-Energy" quantity is in this case the energy of the extended 
Nose-Hoover Hamiltonian?
The .mdp file is atatched (don't be confused about rlist=1 - since Im 
using the Verlet-scheme the verlet-buffer-drift option should be by 
default active and determine the rlist value (Verlet buffer-size) 
automatically).

Best regards,
Bernhard

; Run parameters
integrator    = md        ; leap-frog integrator
nsteps        = 10000000    ; 10000 ps = 10 ns
dt        = 0.001        ; 1 fs
; Output control
nstxout        = 5000        ; save coordinates every ps
nstvout        = 5000        ; save velocities every ps
nstxtcout    = 1000        ; xtc compressed trajectory output every ps
nstenergy    = 1000        ; save energies every ps
nstlog        = 5000        ; update log file every ps
; Bond parameters
continuation    = yes        ; continue from NPT
constraint_algorithm = lincs    ; holonomic constraints
constraints    = h-bonds    ; all bonds (even heavy atom-H bonds) 
constrained
lincs_iter    = 1        ; accuracy of LINCS
lincs_order    = 4        ; also related to accuracy
; Neighborsearching
cutoff-scheme    = Verlet    ; Verlet cutoff-scheme instead of 
group-scheme (no charge-groups used)
ns_type        = grid        ; search neighboring grid cells
nstlist        = 10        ; 10 fs
rlist        = 1.0        ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.0        ; short-range electrostatic cutoff (in nm)
rvdw        = 1.0        ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype    = PME        ; Particle Mesh Ewald for long-range 
electrostatics
pme_order    = 4        ; cubic interpolation
fourierspacing    = 0.12        ; grid spacing for FFT
; Temperature coupling is on
tcoupl        = nose-hoover    ; modified Berendsen thermostat
tc-grps        = Protein Non-Protein    ; two coupling groups - more 
accurate
tau_t        = 2.5    2.5    ; time constant, in ps
ref_t        = 300     300    ; reference temperature, one for each 
group, in K
; Pressure coupling is off
pcoupl        = no         ; no pressure coupling in NVT
; Periodic boundary conditions
pbc        = xyz        ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres    ; account for cut-off vdW scheme
; Velocity generation
gen_vel        = no        ; don’t assign velocities from Maxwell 
distribution