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 -- 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-requ...@gromacs.org.