Date: Fri, 9 May 2003 10:38:53 -0400 (EDT)
From: darden
Subject: Re: Etot is not constant in NVE ensemble

Dear Hungie
For comparison I have run a water box of 216 molecules for 5ns NVE with
either a 2fs or 1fs time step with standard default. I got a drift of
~2kcals per mole out of ~1750 kcals with a 2fs time step and ~1/2 kcal per
mole with a 1fs time step. This latter number is comparable to yours
relative to total energy. Also note that the energy conservation behaved
quadratically with respect to timestep. This does not happen when
you have a serious energy drift problem. However if you're unhappy with
your results....
Tom Cheatham (he can agree/disagree) had a similar problem and found that
by increasing the default ewald coefficient a bit he got rid of the
drift. It might be that for highly charged systems the discontinuity in
direct sum energy/forces at the cutoff [i.e. erfc(beta*r_cutoff)/r_cutoff
is too large] could cause heating. By increasing beta (ewald coefficient)
this problem is alleviated. I don't recall how much Tom C increased beta.
To be correct [i.e. keep a low error vs exact ewald] you then need to
increase reciprocal sum parameters slightly. However since the
reciprocal pme sum is smooth (no discontinuities) it doesn't seem to me
that the reciprocal sum can cause heating. That is your problem is
probably due to direct sum discontinuities and can be helped by increasing
beta.

Hope this help
Tom Darden

On Fri, 9 May 2003, A. Hungie wrote:

> Dear All,
>
> I have run MD simulation of 14 basepairs DNA (in water and nuetralized by
> Na+ ions) using protocol as in Amber tutorial. After doing 6 rounds of
> minimizations (reducing restrain energy from 25 to 0 kcal/mol.A^2), I ran MD
> in NPT ensemble for 100 ps which I obtained canstant box size. Then I
> switched to NVE ensemble for 5.6 ns (800 ps for each run). I found that
> total energy is not constant, i.e. it increases about 6 kcal/mol (from
> -37828 kacl/mol to -37822 kcal/mol) within 5.6 ns. I used Amber 6 on Linux
> cluster and input file as below.
>
> Is the obtained total energy normal and acceptable? Is it equilibrium?
>
> #my input file.
> ------------------------------------------
> &cntrl
> ntx = 7, irest = 1, nmropt = 0,
> imin = 0,
> maxcyc = 1000, ncyc = 5000,
> ntmin = 1, dx0 = 0.1, dxm = 0.5, drms = 0.0001,
>
> nstlim = 800000,
> dt = 0.001, ndfmin = 0, t = 0.0,
> timlim = 999999., ntcm = 1, nscm = 2500,
>
> ntpr = 1000, ntwr = 1000,
> ntwx = 1000, ntwv = 0, ntwe = 0,
> ntwprt = 0, ioutfm = 0, ntrx = 1, ntxo = 1,
>
> cut = 8.0, dielc = 1.0, nsnb = 10,
> scnb = 2.0, scee = 1.2, iwrap = 0,
>
> ntb = 1, ntp = 0, npscal = 1,
> pres0 = 1.0, comp = 44.6, taup = 2.0,
>
> ntc = 2, tol = 0.000001,
> ntf = 2,
>
> ibelly = 0, ntr = 0,
>
> temp0 = 300.0, tempi = 300.0,
> ig = 71277, heat = 0.0, tautp = 4.0,
> ntt = 0, vlimit = 20.0, dtemp = 1.0,
> &end
>
> &ewald
> a = 42.0151161, b = 43.0746222, c = 69.8275997,
> alpha = 90.000, beta = 90.000, gamma = 90.000,
> nfft1 = 45, nfft2 = 45, nfft3 = 72,
> order = 4, dsum_tol = 0.00001, eedmeth = 1,
> opt_infl = 1, vdwmeth = 1, use_pme = 1,
> frc_int = 0, nbflag = 1, skinnb = 1.0
> &end
> -----------------------------------
>
> Thank you very much in advance.
>
> Hungie