VMD-L Mailing List

From: Pavan Ghatty (pavan.vmd_at_gmail.com)
Date: Wed Aug 23 2006 - 08:55:04 CDT

Hello all,

I started off with a protein in water. Minimized it (1000 steps), performed
temperature rescaling (1000 steps) and then a 6ns NPT simulation. All went
well. I then used the output of the 6ns NPT run to perform another 6ns NPT
run (restarting essentially). But then the run crashed complaining about the
"ERROR: Constraint failure in RATTLE algorithm for atom ". I checked in the
archive for related problems and found the following fixes:

On Tue, 5 Apr 2005, Marc Q. Ma wrote:
*> Dear Nicholas, *
*> *
*> timestep 2.0 *
*> stepsPerCycle 8 *
*> nonBondedFreq 2 *
*> fullElectFrequency 4 *
*> *
*> 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:

Now, I tried all the fixes but it just wouldn't run.... at that point I
found that I was using the same input file, which is attached, without
changing the cellbasisvector values. The values initially were 91, 89 and 89
along x, y and z axes. I then found the box dimensions from the output of
the 6ns run and changed them to the new values which were, 87.9, 100.7 and
87.5 . The simulation is now running fine. Now, the question is, should I
have done something else or is what I did the right thing. Also, is there a
better way. Mighty Thanks for any suggestions/comments.

--------------------------------------------------------------
INPUT FILE
---------------------------------------------------------------
firsttimestep 0

# -------------
# Force field
# -------------
paratypecharmm on
parameters ../par_all27_prot_lipid.prm

# ------------------
# Molecular system
# ------------------
structure ../ionxlor.psf

# -----------------
# Initialization
# -----------------
{
  temperature 310.0
  coordinates coordiates.file
  cellBasisVector1 87.9 0.0 0.0
  cellBasisVector2 0.0 100.7 0.0
  cellBasisVector3 0.0 0.0 87.5
  cellOrigin 0.0 0.0 0.0
}
wrapall on

# -------------------
# MD Approximations
# -------------------
switching on # vdw: switch function, charge: shift function
switchdist 10 # cutoff starting to smooth
cutoff 12 # cutoff ending to smooth
pairlistdist 14 # cutoff for nonbonded neighbor list
stepspercycle 2 # neighbour list updating frequency
exclude scaled1-4 # 1-2,1-3,1-4 interactions are excluded
1-4scaling 1.0 # 1.0 for Charmm, 0.8333 for Amber
margin 1.0 # adjust cell size for NPT
#PME on
#PMEGridSizeX 72
#PMEGridSizeY 72
#PMEGridSizeZ 72

# ---------------
# MD integrator
# ---------------
seed 31415 # for starting a MD
timestep 2.0 # integration step, fs
nonbondedFreq 1 # timesteps between nonbonded evaluation
fullElectFrequency 4 # timesteps for full elec. calculation
rigidBonds all # all X-H bonds are constrained
useSettle on # Settle algorithm for water

langevin on # langevin dynamics
langevinTemp 310.0 # target temperature
langevinDamping 1.0 # damping coefficient of 1/ps
{
useGroupPressure yes # needed for rigid bonds
useFlexibleCell no # no for water box, yes for membrane
useConstantArea no # no for water box, maybe for meebrane
useConstantRatio no # x-y ratio constant while x,y,z change

LangevinPiston on
LangevinPistonTarget 1.0325 # pressure in bar -> 1 atm
LangevinPistonPeriod 200 # ocillation period around 200 fs
LangevinPistonDecay 100 # ocillation decay time of 100 fs
LangevinPistonTemp 310.0 # the same as langevinTemp
}

# --------
# output
# --------
dcdfile dcdfile
dcdunitcell yes

restartname restart
restartsave yes

dcdfreq 20
xstfreq 20
restartfreq 50

binaryrestart no # all restart file no binary
binaryoutput no # all outputs no binary

outputname output
outputEnergies 20
outputPressure 20
outputtiming 20

# --------
# script
# --------
run 100

*