Accuracy and stability issues

There are some accuracy and stability issues concerning the molecular dynamics methods described in the Rop and DNA dynamics pages. Here are the relevant parts of a useful discussion I've had on the NAMD mailing list. The gist is that that the proposed integration scheme hovers at the limit of instability. The practical experience, at least with proteins (and with well over 170 nanoseconds worth of total simulation time), is that you may get away with it, but as usually is the case, YMMV ('Your mileage may vary'). Note that if you follow Marc's advice (Message 2 below) the total simulation time will approximately double.

Note (for Message 2) : What we use is more properly called SHAKE-I (which is Verlet-I/r-RESPA with SHAKE (ie. rigidbonds all)). This has been shown in the past to be stable up-to and including a Δt=8fs, see the paper by Izaguirre, Reich & Skeel[1]. For DNA simulations though, it is probably better to follow Marc's advice.

Message 1

From glykos@mbg.duth.gr Tue Apr  5 13:34:00 2005
Date: Tue, 5 Apr 2005 13:33:48 +0300 (EEST)
From: Nicholas M Glykos <glykos@mbg.duth.gr>
To: namd-l@ks.uiuc.edu
Subject: RATTLE and DNA


Dear All,

A student is doing a simulation of the Dickerson DNA dodecamer using the 
scripts that I attach on this message (explicit waters, 22 sodium ions 
placed with the program 'sodium'). The initial minimisation (down to a 
'GRADIENT TOLERANCE: 1.0179') and heating stage appeared to go well. The 
production run went smothly for ~0.7 nanoseconds, and then it died with

______________________________________________________________________

ERROR: Constraint failure in RATTLE algorithm for atom 364!
ERROR: Constraint failure; simulation has become unstable.
ERROR: Exiting prematurely.
==========================================
WallClock: 146567.853583  CPUTime: 142235.640000  Memory: 16432 kB
______________________________________________________________________


He restarted his job, and this time he got to run it for another 1.3 
nanoseconds before the problem re-appeared :

______________________________________________________________________

ERROR: Constraint failure in RATTLE algorithm for atom 375!
ERROR: Constraint failure; simulation has become unstable.
ERROR: Exiting prematurely.
==========================================
WallClock: 255880.024668  CPUTime: 239782.650000  Memory: 14358 kB
______________________________________________________________________


He tried to switch off 'rigidDieOnError', but to no avail. The structure 
(with the exception of the first base pair) appears to be fairly stable 
(and the atoms mentioned above belong to the second base pair, not the 
first). He is currently trying with a lower temperature (300K) [and if 
that fails I suppose that we will have to either reduce the 2fs step, or 
increase 'rigidTolerance'].

What buffles me, is that the simulation protocol he is using is fairly 
standard (and can be found in numerous papers), so I didn't expect to have 
these problems. Any suggestions would be greatly appreciated.

Regards,
Nicholas

Message 2

From qma@oak.njit.edu Tue Apr  5 17:34:33 2005
Date: Tue, 05 Apr 2005 10:01:35 -0400
From: Marc Q. Ma <qma@oak.njit.edu>
To: Nicholas M Glykos <glykos@mbg.duth.gr>
Cc: namd-l@ks.uiuc.edu
Subject: Re: namd-l: RATTLE and DNA

    [ The following text is in the "ISO-8859-1" character set. ]
    [ Your display is set for the "ISO-8859-7" character set.  ]
    [ Some special characters may be displayed incorrectly. ]

Dear Nicholas,

timestep2.0
stepsPerCycle8
nonBondedFreq2
fullElectFrequency4

What it gives is as follows:

1. You are using the Impulse multiple time stepping integrator, also 
called Verlet-I/r-RESPA. This integrator has been shown to be stable 
when outer most time step is 6fs with rigid bonds and rigid waters. 
However, for your long simulation, ~1ns, the numerical resonance and 
nonlinear resonance will make the simulation very unstable even at 6fs. 
For these long simulations, I would suggest to use outermost timestep 
not more than 4fs. (FYI - I have some publications on the stability 
issues on the Impulse integrator, eg.
Q. Ma, J. Izaguirre and R. Skeel, Verlet-I/r-RESPA/Impulse is limited 
by nonlinear instability, SIAM J. Sci. Comput., 24(6):1951-1973, 2003.

The discussion in the paper was made for constant energy simulation. 
However, it gives good indication for other ensembles as well.

2. The above configuration gives you: leapfrog (innermost integrator) 
2fs, Impulse level 2 (middle integrator) for short range nonbonded 
forces (LJ and Coulombic) 2.0*2 = 4fs, and Impulse level 3 (outermost 
integrator) for long range nonbonded forces, 2.0*4 = 8fs!

For Impulse, 8fs timestep is not attainable even with rigid water.

3. Simple fix --> change to timestep 1.0 and keep everything else. This 
should make your outermost timestep to 4, and should be OK for long 
simulations.

4. Or, keep timestep 2.0, make nonBondedFreq 1, and fullElectFreq 2. 
This should make your simulation go faster although asymptotically the 
running time is on the same order.

Hope this helps. Good luck!

Marc

Message 3

From glykos@mbg.duth.gr Tue Apr  5 18:02:34 2005
Date: Tue, 5 Apr 2005 18:02:24 +0300 (EEST)
From: Nicholas M Glykos <glykos@mbg.duth.gr>
To: Marc Q. Ma <qma@oak.njit.edu>
Cc: namd-l@ks.uiuc.edu
Subject: Re: namd-l: RATTLE and DNA



Dear Marc,

Thank you for clarifying this for me. I probably misinterpreted the 
relevant part from NAMD documentation :


* fullElectFrequency < number of timesteps between full electrostatic evaluations >
Acceptable Values: positive integer factor of stepspercycle
Default Value: nonbondedFreq
Description: This parameter specifies the number of timesteps between each 
full electrostatics evaluation. It is recommended that fullElectFrequency 
be chosen so that the product of fullElectFrequency and timestep does not 
exceed  4.0  unless rigidBonds all or molly on is specified, in which case 
the upper limit is perhaps doubled.

This 'perhaps doubled for rigidBonds' made me thing that I could get away 
with it, but obviously this is not the case.

Thanks again,
Nicholas

Message 4

From qma@oak.njit.edu Tue Apr  5 18:25:22 2005
Date: Tue, 05 Apr 2005 11:23:09 -0400
From: Marc Q. Ma <qma@oak.njit.edu>
To: Nicholas M Glykos <glykos@mbg.duth.gr>
Subject: Re: namd-l: RATTLE and DNA

    [ The following text is in the "ISO-8859-1" character set. ]
    [ Your display is set for the "ISO-8859-7" character set.  ]
    [ Some special characters may be displayed incorrectly. ]

Nicholas,

You are welcome!

Stability is always a big issue in MD simulations. To be always on the  
safe side and for production runs, the time steps should be chosen more  
conservative unless you are doing some experiments on the numerical  
schemes.

Impulse integrator is always troubled by the linear and nonlinear  
instabilities. For short (defined relatively) simulations and small  
(again defined relatively) systems, the stability issue is not obvious.  
However, when simulating larger systems and longer time, the stability  
issues become obvious and has to be taken good care of.

You can try hybrid Monte Carlo simulations if you only want to sample  
the system. Or, stop the simulations at 1ns and re-generate an initial  
velocity distribution.

If you use smaller time steps, you can obviously go longer time without  
all these troubles.

Best,

Marc