[gmx-users] Re: NVE problems in gmx-users digest, Vol 1 #1174

[istvan at kolossvary.hu]

Hi, 

With Michael Shirts' help I made some progress, hope this is useful for all 
NVE "enthusiasts". BTW, I am trying to use NVE for spectral analysis (see my 
other thread on spectral density in Vol 1 #1185, topic 5). Here is a set of 
parameters that I believe are quite conservative and I hoped I could run 
long NVE simulations on a small, totally unconstrained (except settle on the 
waters) system: alanine dipeptide in a 3x3x3 nm box.
title               =  Yo
cpp                 =  /lib/cpp
constraints         =  none
integrator          =  md
dt                  =  0.00025  ; ps !
nsteps              =  4000000  ; total 1000 ps.
nstcomm             =  1
comm_grps           =  Protein  SOL
nstxout             =  1000
nstvout             =  1000
nstfout             =  0
nstlog              =  1000
nstenergy           =  100
nstlist             =  10
ns_type             =  grid
pbc                 =  xyz
rlist               =  1.1
; Long-range
coulombtype         =  pme
pme_order           =  6
fourierspacing      =  0.1                      ; 1 angstrom grid spacing
rcoulomb            =  1.0                      ; 10 A cut-off for coulomb
ewald_rtol          =  2.20905e-5
vdw-type            =  Switch
rvdw-switch         =  0.9
rvdw                =  1.0
; NVE - no temperature control
Tcoupl              = no
tc-grps                 =  Protein      SOL
tau_t               =  0.1      0.1
ref_t               =  300      300
; Energy monitoring
energygrps          =  Protein  SOL
; NVE - no pressure control
Pcoupl              = no
Pcoupltype          =  isotropic
tau_p               =  0.5
compressibility     =  4.5e-5
ref_p               =  1.0
; Generate velocites is off at 300 K.
gen_vel             =  no
gen_temp            =  300.0
gen_seed            =  173529
; CONSTRAINTS
constraints              = none
; constraint-algorithm     = Lincs
; unconstrained-start      = no
; lincs_order              = 4
; lincs_warnangle          = 30
morse                    = yes
; Shake-SOR                = no
; shake-tol                = 1e-04 

ewald_rtol=2.20905e-5 comes from the fact that in Gromacs 
erfc(sigma*rcutoff)=ewald_rtol and with a 1nm rcutoff I used a rather 
conservative heuristic of rcutoff = 3/sigma. Also note that with 
rlist>(rcoulomb, rvdw) grompp gives an error message that can be commented 
out in readir.c (info from Michael Shirts). 

A good way of testing energy conservation is using delta_E = log(1/N * 
Sum|(Ei-E0)/E0|), see JCP 103, 1995, 9444-9459. Ei is the current total 
energy, E0 is the initial total energy and N is the number of sample points 
on the stretch of the trajectory considered. delta_E < -2.5 is indicative of 
a stable constant energy simulation. With the above settings I was able to 
keep delta_E below -2.5 for 33 ps. I guess this can be good or this can be 
bad depending on your perspective. As far as I am concerned I am rather 
disappointed that after only 33 ps the simulation should be considered 
unstable. Can anybody comment on
1. Whether I missed something and even more conservative settings should be 
used for NVE,
2. It concerns me a lot that we typically run hundreds or thousands of 
picoseconds with T/P coupling whereas the underlying integrator cannot 
conserve energy beyond a few tens of picoseconds. 

Many thanks, 

  Istvan Kolossvary