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. For DNA simulations though, it is probably better to follow Marc's advice.
From email@example.com Tue Apr 5 13:34:00 2005 Date: Tue, 5 Apr 2005 13:33:48 +0300 (EEST) From: Nicholas M Glykos <firstname.lastname@example.org> To: email@example.com 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
From firstname.lastname@example.org Tue Apr 5 17:34:33 2005 Date: Tue, 05 Apr 2005 10:01:35 -0400 From: Marc Q. Ma <email@example.com> To: Nicholas M Glykos <firstname.lastname@example.org> Cc: email@example.com 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
From firstname.lastname@example.org Tue Apr 5 18:02:34 2005 Date: Tue, 5 Apr 2005 18:02:24 +0300 (EEST) From: Nicholas M Glykos <email@example.com> To: Marc Q. Ma <firstname.lastname@example.org> Cc: email@example.com 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
From firstname.lastname@example.org Tue Apr 5 18:25:22 2005 Date: Tue, 05 Apr 2005 11:23:09 -0400 From: Marc Q. Ma <email@example.com> To: Nicholas M Glykos <firstname.lastname@example.org> 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