[gmx-users] NVE Vacuum Fluctuation of Energy

[Johnny Lu]

Dear Users,

I was trying to run NVE simulation of alanine peptide in vacuum with double
precision gromacs 5.1. [plot of energy at
http://postimg.org/image/6b6kksax7/]
The total energy of the simulation still has fluctuation about 1kcal/mol.
There should be no fluctuation of total energy in an ideal world.

What are the possible source of this fluctuation? I'm testing the error
contributed by the integrator by halving the time-step.

How big is the error caused by integrator, error caused by h-bond
constraint, and fluctuation caused by thermostat relative to each other?

Previously I used the same mdp file, except with 1fs time step and with
H-bond constraint, to simulate alanine peptide in water. The total energy
of the whole system seemed conserved. I now wonder if that conservation is
actually a cancellation of different numerical errors.

When I run alanine peptide in vacuum, H-bond constraint made the energy
drop quickly. I set constraint to none, and the energy drop was lower.

The question I'm interested in is verification of a theory. It is affected
by the shape and pattern of these errors. Without knowing the pattern of
these errors, I hope the errors are small enough.

Thank you again.

The mdp file is below. the box size is large enough such that the distance
between mirror image of the molecules is larger than the cut off for vdw
and coulomb forces.

title        = Alanine
; Run parameters
integrator    = md        ; leap-frog integrator
nsteps        = 300000000    ; 300 ns
dt        = 0.0005        ; 0.5 fs

; Output control
nstxout        = 10            ; save coordinates every 1.0 ps
nstvout        = 10            ; save velocities every 1.0 ps
nstfout        = 10

nstenergy    = 500            ; save energies every 1.0 ps
nstlog        = 500        ; update log file every 1.0 ps

; Bond parameters
continuation            = no        ; Restarting after NPT equilibration 1
constraint_algorithm    = lincs        ; holonomic constraints
constraints                = none    ; all bonds (even heavy atom-H bonds)
constrained
lincs_iter                = 2            ; accuracy of LINCS
lincs_order                = 4            ; also related to accuracy

; Neighborsearching
cutoff-scheme   = Verlet
ns_type            = grid        ; search neighboring grid cells
nstlist            = 20        ; 20 fs, largely irrelevant with Verlet
scheme
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.1        ; grid spacing for FFT

tcoupl        = no                ; no thermostat
pcoupl                = no        ; no barostat


; Periodic boundary conditions
pbc        = xyz        ; 3-D PBC

; Dispersion correction
DispCorr    = EnerPres    ; account for cut-off vdW scheme

; Velocity generation
gen_vel        = yes        ; Velocity generation is off

verlet-buffer-drift = 0.00000005