[gmx-users] NVE Vacuum Fluctuation of Energy
Vitaly V. Chaban
vvchaban at gmail.com
Wed Sep 2 20:16:04 CEST 2015
The major source is rounded numbers -- in all functions, not only in
the motion propagation algorithm.
On Tue, Sep 1, 2015 at 11:20 PM, Johnny Lu <johnny.lu128 at gmail.com> wrote:
> Dear Users,
> I was trying to run NVE simulation of alanine peptide in vacuum with double
> precision gromacs 5.1. [plot of energy at
> 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)
> 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
> 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
> 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
> Gromacs Users mailing list
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.
More information about the gromacs.org_gmx-users