Re: RATTLE and DNA

From: Nicholas M Glykos (glykos_at_mbg.duth.gr)
Date: Tue Apr 05 2005 - 10:02:24 CDT

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

On Tue, 5 Apr 2005, Marc Q. Ma wrote:

> Dear Nicholas,
>
> timestep2.0
> stepsPerCycle8
> nonBondedFreq2
> fullElectFrequency4
>
> The above configuration will not work at all. 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
> On Apr 5, 2005, at 6:33 AM, Nicholas M Glykos wrote:
>
> >
> > 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
> >
> >
> > --
> >
> >
> > Dr Nicholas M. Glykos, Department of Molecular
> > Biology and Genetics, Democritus University of Thrace,
> > Dimitras 19, 68100 Alexandroupolis, GREECE, Fax ++302551030613
> > Tel ++302551030620 (77620), http://www.mbg.duth.gr/~glykos/
> > <heat.namd><equi.namd>
>

-- 
            Dr Nicholas M. Glykos, Department of Molecular 
        Biology and Genetics, Democritus University of Thrace,
    Dimitras 19, 68100 Alexandroupolis, GREECE, Fax ++302551030613
     Tel ++302551030620 (77620),  http://www.mbg.duth.gr/~glykos/

This archive was generated by hypermail 2.1.6 : Wed Feb 29 2012 - 15:40:38 CST