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

[Anton Feenstra]

istvan at kolossvary.hu wrote:
> 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.
[...]
> 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,

Are you running in double precision? I have found that in single precision,
setting the timestep very short actually increases the energy drift, possibly
due to increased accumulation of numerical addition errors in the integration.
I've seen it begin at about 0.5-0.1fs timestep! This effect obviously will
also exist in double precision, but then at <1e-6 fs... ;-)

While PME is great for accurate long-range forces, it still has a transition
aread from 'real-space' to 'reciprocal-space' interactions around the cut-off,
where still discontinuities may be in the potential, or its derivatives. Not
sure how (if?) that may be optimized further...

I've used shifted potentials instead of a switch function, although I was
only interested in measuring energy conservation and not in the results of
that particular simulation... Shifted potentials, I'm told, give 'artefacts',
as in J.Phys.Chem. 98:8785-8792 1994 and J.Chem.Phys. 95:8430-8441 1991.

I've user rms-drift over several simulations to measure energy conservation
accurately. Drift is as reported by g_energy from a linear fit to the energy,
then rms averaged over 5-10 runs. I'm not sure how drift would compare to your
measure, though.

> 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.

For spc216 water box at 0.5fs timestep, I've found drifts of around 0.01 kJ/mol/ps,
and still only about 0.1 kJ/mol/ps at 2 fs. Values are essentially the same for
a 85 residue protein in a water box.

> Many thanks,
>  Istvan Kolossvary
> _______________________________________________
> gmx-users mailing list
> gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-users-request at gromacs.org.


-- 
Groetjes,

Anton
  _____________ _______________________________________________________
|             |                                                       |
|  _   _  ___,| K. Anton Feenstra                                     |
| / \ / \'| | | Dept. of Pharmacochem. - Vrije Universiteit Amsterdam |
|(   |   )| | | De Boelelaan 1083 - 1081 HV Amsterdam - Netherlands   |
| \_/ \_/ | | | Tel: +31 20 44 47608 - Fax: +31 20 44 47610           |
|             | Feenstra at chem.vu.nl - www.chem.vu.nl/~feenstra/       |
|             | "If You See Me Getting High, Knock Me Down"           |
|             | (Red Hot Chili Peppers)                               |
|_____________|_______________________________________________________|