[gmx-users] "ideal" mdp file for NVE simulation as textbook example

[Anton Feenstra]

Marc Baaden wrote:

> Hi David,
> 
> 
>>>>David said:
> 
>  >> We have done quite a few simulations of water clusters in vacuum with
>  >> settle that give perfect energy conservation (of course there are small
>  >> fluctuations, but no trend like you have here). For these systems we are
>  >> lacking all the bonded interactions of course. Maybe you want to try to
>  >> reduce your time step to 0.2 fs (unless you have dummies). You should
>  >> then do the same simulation length (10 ps) and compare.
> 
> thanks for the suggestion. I tried it and indeed this seems a clear source
> of energy drift. I observe the following tendency:
> 
> timestep  energy drift (over 10 ps)
> --------  ------------
> 0.1 fs    - 0.04 kJ/mol
> 0.2 fs    - 0.15  "
> 1   fs    - 4.0   "
> 2   fs    - 11.7  "
> 
> I'd be tempted to conclude that constant energy simulations would thus be
> rather utopic for any protein/peptide system bigger than a couple of
> residues and aiming at roughly nanosecond timescale. Or are there other
> tricks that would work (eg dummies .. but I think they'd only buy an
> order of magnitude) ?

Do also look at JComputChem1999,20,786 (if you haven't already), 
especially Figure 2 which shows drift versus dt for water (spc216) and a 
protein in vacuum (805 atoms). There are several factors to consider.

First, in order to actually observe the drift vs dt at low dt (<2fs), 
several 'tricks' are needed to reduce the noise from other sources, but 
I now see that you got all these already (double precision, shift 
function (no cutoff, in your case), neighborlist update each step, and 
no P or T coupling). I think, but cannot find it in my paper, that I 
also had to use a LINCS order of 8.

(actually, there was another interesting observation that, even at 
double precision, at very low dt (somewhere below 0.1fs), you will see 
the drift *increasing* again due to numerical errors of adding very 
small changes to atomic coordinates.)

Second, energy drift is diffusive, so you'll need several runs and 
average. (For this, also check two papers by A.K.Mazur, 
JComputPhys1997,136,354 and JPhysChem1998,102,473.) I'd say the average 
over 10x1ps would be better than the trend over 1x10ps.

Finally, the interesting parameter is actually the order of the drift 
vs. dt (as in the plot). You'd want that at 2nd order (which is, roughly 
speaking, the order of acurracy of Leap-Frog in energy), and not higher. 
You can see in the plot, that the dummy atoms increase the timestep at 
which a 2nd to higher order transition occurs.

I don't know whether comm removal (-1, angular) would do anything, I 
think it shouldn't but you may check...


-- 
Groetjes,

Anton

* NOTE: New Phone & Fax numbers (below) *

  _____________ _______________________________________________________
|             |                                                       |
|  _   _  ___,| K. Anton Feenstra                                     |
| / \ / \'| | | Dept. of Pharmacochem. - Vrije Universiteit Amsterdam |
|(   |   )| | | De Boelelaan 1083 - 1081 HV Amsterdam - Netherlands   |
| \_/ \_/ | | | Tel: +31 20 59 87608 - Fax: +31 20 59 87610           |
|             | Feenstra at few.vu.nl - www.few.vu.nl/~feenstra/         |
|             | "If You See Me Getting High, Knock Me Down" (RHCP)    |
|_____________|_______________________________________________________|