Coarse grained instability: atoms moving too fast

From: Marlon Sidore (marlon.sidore_at_gmail.com)
Date: Mon Oct 05 2015 - 06:47:41 CDT

Hello,

I am experiencing an "ERROR: Atoms moving too fast; simulation has become
unstable (x atoms on patch y pe z)." error which makes simulations crash
after an undetermined amount of steps (from a few thousands to tens of
millions).

The system is a MARTINI coarse-grained membrane system with a membrane
protein. The atoms moving too fast can be any kind of bead, water, lipid or
protein, although it seems to happen more frequently on a protein bead.

In order to reduce the crashes, I had to reduce the integration step from
20 (usual with MARTINI) to 10fs, but it still crashes at least once every
50 million steps (and that's very optimistic).
The crashes are also less frequent when using a smaller integration step.

Even though I can restart the simulation without crashing at the same
moment, I would really like to understand why it is happening and try to
fix it (the crashes seem pretty random), so I can get back to a 20fs
timestep.

Is there something specific to NAMD that could cause these crashes ? (as
GROMACS doesn't seem to crash like that, however I can't use GROMACS for
that)

Here is a classic input file I use (it has restraints and uses colvars, but
the crashes also happen without):
#############################################################
## JOB DESCRIPTION ##
#############################################################

# ha8a

#############################################################
## ADJUSTABLE PARAMETERS ##
#############################################################
set inName #NAME#-#PREVIOUS#
set outName #NAME#-#NEXT#

set outputname $outName
set restart 1 ;#0 for new sim, 1 for continuation

set pvmode "p" ;#p to activate pression control
set temode "t" ;#t to activate temperature control

proc get_first_ts { xscfile } {
  set fd [open $xscfile r]
  gets $fd
  gets $fd
  gets $fd line
  set ts [lindex $line 0]
  close $fd
  return $ts
}

set temp 325.0
cosAngles on

structure DimerAqpZ.psf
coordinates edge_init.pdb
bincoordinates $inName.restart.coor
binaryoutput yes ;# Use flipinpdb to use the .coor
binaryrestart yes

if {$restart == 1} {
  binvelocities $inName.restart.vel
  extendedSystem $inName.restart.xsc
  set currenttimestep [get_first_ts $inName.restart.xsc]
  COMMotion yes
} else {
  temperature $temp
  set currenttimestep 0
  COMmotion no ;#Deffault for JH and out of the loop
}

firsttimestep $currenttimestep

#############################################################
## SIMULATION PARAMETERS ##
#############################################################

# Input
paraTypeCharmm on
parameters
/scratch/cnt0023/ism6198/sidorem/martini-par/martini-protein-bonds.par
parameters
/scratch/cnt0023/ism6198/sidorem/martini-par/martini-protein-angles-cos.par
parameters
/scratch/cnt0023/ism6198/sidorem/martini-par/martini-protein-dihedrals.par
parameters /scratch/cnt0023/ism6198/sidorem/martini-par/martini-all-nonb.par
parameters
/scratch/cnt0023/ism6198/sidorem/martini-par/martini-lipids-bonds-angles-dihedrals.par

# Force-Field Parameters
exclude 1-2
#1-4scaling 1.0 ;#useful when exclude == 1-4
cutoff 12.0
martiniSwitching on
PME off ;# Default off
switching on
switchdist 9.0
pairlistdist 24
dielectric 15.0

# Integrator Parameters
timestep 10.0
nonbondedFreq 1 ;# Default 1
stepspercycle 20

#Constraints and restraints

#Constraints and restraints

fixedAtoms off ;# switch to off after PR simulation
#fixedAtomsFile mono1fixedatomshelicesBAS.cnst ;# Ref atoms and positions
for PR
#fixedAtomsCol O ;# Read occupancy field for constrained atoms A value of 0
in the fcorresponding field indicates that the atom is not fixed
extraBonds yes ;#For elastic network
extraBondsFile ElastDimer.dat #For elastic network
constraints on
consref edge_init.pdb
conskfile harmonichelicesBAS.cnst
conskcol O
consexp 2
constraintScaling 1

# CONSTANT-T
if {$temode == "t"} {
langevin yes ;# do langevin dynamics
langevinDamping 1 ;# damping coefficient(gamma)5/ps (1 pour
production)
langevinTemp $temp
langevinHydrogen off ;# don't couple langevin bath to hydrogens
}

# Constant Temperature Control ONLY DURING EQUILB
#reassignFreq 500; # reassignFreq: use this to
reassign velocity every 500 steps

# Periodic Boundary Conditions

cellBasisVector1 200.0 0.0 0.0
cellBasisVector2 0.0 208.0 0.0
cellBasisVector3 0.0 0.0 93.0
cellOrigin 0.0 0.0 0.0
wrapWater on
wrapAll on

#
#PME yes
#PMEGridSizeX 256
#PMEGridSizeY 256
#PMEGridSizeZ 256

# CONSTANT-P
useGroupPressure no ;# yes needed for rigid bonds and yes for JH; rigid
bonds not in used for CG
useFlexibleCell yes ;# yes for membrane anf yes for JH
#useConstantArea yes ;# commented for JH
useConstantRatio yes ;# yes for JH

if {$pvmode == "p"} {
langevinPiston on
langevinPistonTarget 1.01325 ;# in bar -> 1 atm
#one may need to bump up the pressure constants at first
langevinPistonPeriod 2000 ;#usually 2000
langevinPistonDecay 1000 ;#usually 1000
#langevinPistonPeriod 500 ;# JH setting
#langevinPistonDecay 250 ;# JH setting
langevinPistonTemp $temp
}
#ExcludeFromPressure on ;# permet d'exclure des atomes (grains) de
la pression, pourrait être utile
#ExcludeFromPressureFile Poccupancy.cnst ;# 0 pour not excluded, 1 pour
excluded
#ExcludeFromPressureCol O ;# reads occupancy

# Output
outputname $outName
restartname $outName.restart
restartfreq 10000 ;#sécurité toutes les 5000 steps
#restartsave yes ;#pour obtenir une chronologie dans les fichiers
.restart, utile si le dernier .restart est à une step trop lointaine du
crash

DCDfile $outName.dcd
DCDfreq 10000

outputTiming 10000 ;#CPU & wallclock times, energy and pressure
output to stdout
outputEnergies 10000
outputPressure 10000

XSTfile $outName.xst
XSTFreq 10000

#############################################################
## EXTRA PARAMETERS ##
#############################################################

#Colvars
colvars on
colvarsConfig pl60.in ;#Commencer une simulation
colvarsInput $inName.restart.colvars.state ;#Continuer une
simulation

#############################################################
## EXECUTION SCRIPT ##
#############################################################

run #EXTEND#

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:22:06 CST