[gmx-users] Energy conservation issue
Yujie Wu
yujie.wu at hec.utah.edu
Thu Oct 27 08:08:03 CEST 2005
I am having an energy conservation problem with my system, and I am not
sure exactly what causes the problem. Previously, I thought it could be
due to the single-precision. But I just noticed the recent GROMACS paper
published on JCC, saying:
(quote)"By taking care of summation order and sometimes using
intermediate double precision variables it has been possible to use
single precision by default while still conserving energy perfectly in
microcanonical ensemble runs."
This is good (if true), but it worries me because I am afraid that I did
something stupid that caused the energy drifting. To confirm the problem,
I ran an NVE simulation on a 216-SPC-water system (version 3.1.4, single
precision, serial). The energy conservation is bad (-32 kJ/mol/ns). For
the same system and equivalent simulation conditions, DL_POLY (which
uses double precision) can give really almost perfect conservation. I
still cannot see what is wrong with my GROMACS run, but I could be
easily blind to my own mistakes. So I post my mdout.mdp file below, and
I would appreciate it if anyone please points out the problem. Many
thanks in advance.
Sincerely,
Yujie Wu
--------------------------- mdout.mdp file pasted below -------------------------------
; RUN CONTROL PARAMETERS =
integrator = md
; start time and timestep in ps =
tinit = 0
dt = 0.002
nsteps = 3250000
; mode for center of mass motion removal =
comm-mode = Linear
; number of steps for center of mass motion removal =
nstcomm = 0
; group(s) for center of mass motion removal =
comm-grps = SOL
; LANGEVIN DYNAMICS OPTIONS =
; Temperature, friction coefficient (amu/ps) and random seed =
bd-temp = 300
bd-fric = 0
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS =
; Force tolerance and initial step-size =
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells =
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints =
fcstep = 0
; Frequency of steepest descents steps when doing CG =
nstcgsteep = 1000
; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 500
nstvout = 500
nstfout = 500
; Output frequency for energies to log file and energy file =
nstlog = 500
nstenergy = 500
; Output frequency and precision for xtc file =
nstxtcout = 0
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps =
; Selection of energy groups =
energygrps =
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist = 1
; ns algorithm (simple or grid) =
ns-type = simple
; Periodic boundary conditions: xyz or no =
pbc = xyz
; nblist cut-off =
rlist = 0.9
domain-decomposition = no
; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 0.9
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon-r = 1
; Method for doing Van der Waals =
vdw-type = Cut-off
; cut-off lengths =
rvdw-switch = 0
rvdw = 0.9
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid =
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used =
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters =
pme_order = 6
ewald_rtol = 1e-06
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling =
tcoupl = No
; Groups to couple separately =
tc-grps = SOL
; Time constant (ps) and reference temperature (K) =
tau-t = 0.5
ref-t = 298.15
; Pressure coupling =
Pcoupl = No
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau-p = 1
compressibility = 4.5e-5
ref-p = 1.01325
; SIMULATED ANNEALING CONTROL =
annealing = no
; Time at which temperature should be zero (ps) =
zero-temp_time = 0
; GENERATE VELOCITIES FOR STARTUP RUN =
gen-vel = yes
gen-temp = 298.15
gen-seed = 173529
; OPTIONS FOR BONDS =
constraints = all-bonds
; Type of constraint algorithm =
constraint-algorithm = Lincs
; Do not constrain the start configuration =
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations =
Shake-SOR = no
; Relative tolerance of shake =
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix =
lincs-order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials =
morse = no
; ENERGY GROUP EXCLUSIONS =
; Pairs of energy groups for which all non-bonded interactions are excluded =
energygrp_excl =
; NMR refinement stuff =
; Distance restraints type: No, Simple or Ensemble =
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal =
disre-weighting = Conservative
; Use sqrt of the time averaged times the instantaneous violation =
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file =
nstdisreout = 100
; Orientation restraints: No or Yes =
orire = no
; Orientation restraints force constant and tau for time averaging =
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file =
nstorireout = 100
; Free energy control stuff =
free-energy = no
init-lambda = 0
delta-lambda = 0
sc-alpha = 0
sc-sigma = 0.3
; Non-equilibrium MD stuff =
acc-grps =
accelerate =
freezegrps =
freezedim =
cos-acceleration = 0
; Electric fields =
; Format is number of terms (int) and for all terms an amplitude (real) =
; and a phase angle (real) =
E-x =
E-xt =
E-y =
E-yt =
E-z =
E-zt =
; User defined thingies =
user1-grps =
user2-grps =
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
More information about the gromacs.org_gmx-users
mailing list