[gmx-users] NVE Temperature Drift

[Johnny Lu]

Wow. Thanks! I never would think of that.
I was thinking that it would be fine as long as the cutoff was long enough,
as some people use infinite cutoff.
The anton comparison of several forcefield ("Systematic Validation of
Protein Force Fields against Experimental Data") used 0.95nm for cut off
for several different force-field.
May be that is because their hardware can only do 0.95nm cutoff, or the way
that the cut off is applied.

I decided to follow the Amber99SB-ILDN paper ("Improved side-chain torsion
potentials for the Amber ff99SB protein force field") and set 10 fs
neighbor list update, pair list cutoff 1.05nm, and both nonbonded cutoff 1
nm, and > 10594 TIP3P water.

; Amber99SB MD
; Run parameters
integrator    = md        ; leap-frog integrator
nsteps        = 32768        ; 2 * 100000 = 200 ps, 1 ns
dt        = 0.001        ; 1 fs
; Output control
nstfout        = 2
nstxtcout    = 2        ; xtc compressed trajectory output every 2 fs
nstlog        = 1000        ; update log file every 2 ps
xtc-grps=protein;

continuation    = yes        ; Restarting after NPT
constraint_algorithm = Lincs    ; holonomic constraints
constraints    = h-bonds    ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter    = 2        ; accuracy of LINCS
lincs_order    = 4        ; also related to accuracy
; Neighborsearching
ns_type        = grid        ; search neighboring grid cells
nstlist        = 10        ; 10 fs
rlist        = 1.05        ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1            ; short-range electrostatic cutoff (in nm)
rvdw        = 1            ; 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
optimize_fft    = yes

tcoupl                   = no
pcoupl                   = no

compressibility = 4.5e-5    ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc        = xyz        ; 3-D PBC
; Dispersion correction
DispCorr    = EnerPres    ; account for cut-off vdW scheme
; Velocity generation
gen_vel        = no        ; Velocity generation is off
verlet-buffer-drift=-1
cutoff-scheme = Verlet;

I also ran a long simulation with the old mdp file, and without the
constraint.
The temperature fluctuation was still about 5 kelvin and stay near 300K.
I guess that is because the water box is still not large enough? Wikipedia
says NVE should approach NPT at the thermodynamic limit.
And mentioned the instaneous temperature v.s. the ensemble average of
temperature. How much of this fluctuation of temperature is real? And how
much of that is caused by molecular dynamics?

The alternating of NVE and NPT step gives near 300K for the NPT step, and
strangely the NVE temperature was around 293K (
http://oi60.tinypic.com/34fjokg.jpg).