![[Home]](/wiki.png)
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