AMBER Archive (2008)

Subject: AMBER: What's making my system lose total energy?

From: David Cerutti (dcerutti_at_mccammon.ucsd.edu)
Date: Mon Mar 17 2008 - 19:20:30 CDT

I thought that this input would order up a very nice, stable simulation
in the NVE ensemble. Now, I find that the system (256 glycerol
molecules) is actually losing potential energy at the rate of roughly 1
kcal/60 ps. The typical problem with NVE simulations is that they gain
energy (significantly) if the time step is set too high. But, the
opposite is true here. The input file:

>>>>>
An NVE MD run file
  &cntrl

   nmropt = 0,
   ntx = 5, irest = 1, ntrx = 1, ntxo = 1,
   ntpr = 100, ntwx = 10000, ntwv = 0, ntwe = 1000,
   iwrap = 0, ioutfm = 0,

   ntf = 2, ntb = 1,
   cut = 9.0,

   ibelly = 0, ntr = 0,

   imin = 0,
   nstlim = 20000000,
   nscm = 1000,
   t = 0.0, dt = 0.001,

   temp0 = 298.15, tempi = 298.15,
   ig = 9376841,
   ntt = 0,
   gamma_ln = 0.0,
   vlimit = 20.0,

   ntp = 0, pres0 = 1.0, comp = 44.6,
   taup = 1.0,

   ntc = 2, tol = 0.000005, watnam = 'WAT ',

  &end
<<<<<

Note that I have set the SHAKE tolerance to be half the default value,
along with a conservative time step and interaction cutoff. Here is some
sample output:

  NSTEP = 100 TIME(PS) = 25182.850 TEMP(K) = 307.65 PRESS = 0.0
  Etot = 4401.8340 EKtot = 2659.7186 EPtot = 1742.1154
  BOND = 417.5512 ANGLE = 1388.9901 DIHED = 601.5975
  1-4 NB = 363.9760 1-4 EEL = 16825.8254 VDWAALS = -2246.5884
  EELEC = -15609.2364 EHBOND = 0.0000 RESTRAINT = 0.0000
  Ewald error estimate: 0.7533E-05

...

  NSTEP = 106400 TIME(PS) = 25289.150 TEMP(K) = 296.16 PRESS = 0.0
  Etot = 4398.9124 EKtot = 2560.4280 EPtot = 1838.4844
  BOND = 407.1341 ANGLE = 1415.2345 DIHED = 607.2740
  1-4 NB = 348.7719 1-4 EEL = 16857.0617 VDWAALS = -2249.7184
  EELEC = -15547.2735 EHBOND = 0.0000 RESTRAINT = 0.0000
  Ewald error estimate: 0.2851E-04

Note the difference in the values of Etot. In an earlier run with this
script, I had conducted 20ns of NVE MD (for a diffusion constant
calculation of the glycerol) and seen the temperature go down about 15K
over the course of the run. Naturally, I assumed that this was a problem
with the system not being settled into a true potential energy equilibrium
at the start of the run, meaning that I would have to run more NPT
dynamics and then some more NVT to make sure I had the right box volume
and a stable amount of potential energy before starting NVE. So, I did
all that, and found that NVT was giving stable levels of potential energy.
That's when I found that it was this issue of losing TOTAL energy all
along. Can anyone help?

Thanks!

Dave
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu