[gmx-users] NVE Vacuum Fluctuation of Energy

[Vitaly V. Chaban]

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
> 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-request at gromacs.org.