[MMTK] constrained bonds - nan problem

Qadir Timerghazin qadirt at uwaterloo.ca
Wed May 6 01:59:19 UTC 2009

Hi Konrad,

Recently, I've been trying to perform simulations with some bonds constrained,
but it looks like there is an elusive bug which leads to nan's during
the integration (or segfaults at some systems, 64-bit Linux mostly). The bug is
quite subtle - it shows up and disappears depending on the molecular system,
force field, restraints applied, thermostat and the computer/OS.  When a
trajectory is written during the integration it always leads to nan's, but at
different times. The following usually 'helps' to get nan's:

* Dihedral restraints
* Nose thermostat
* Repeated calls of the integrator

Finally, I was able to come up with a reproducible example to demonstrate the
problem (attached). In this example, if a trajectory file is written, I
consistently get nan's (Mac OS, 32-bit Linux) or segfault (64-bit Linux). If I
remove the trajectory output as action, it works OK on the most of the computers
I tested it. (on Mac OS it takes some time before the nan's show up).

Thanks a lot!

Dr. Qadir K. Timerghazin,
NSERC Postdoctoral Fellow

Department of Chemistry
University of Waterloo
200 University Avenue West
Waterloo, ON, Canada N2L 3G1


from MMTK import *
from MMTK.Proteins import Protein
from MMTK.ForceFields import Amber99ForceField
from MMTK.Dynamics import VelocityVerletIntegrator

from MMTK.Trajectory import Trajectory, TrajectoryOutput, StandardLogOutput
from MMTK.Minimization import SteepestDescentMinimizer
from MMTK.ForceFields.Restraints import HarmonicDihedralRestraint
from MMTK.Environment import NoseThermostat

universe = InfiniteUniverse(Amber99ForceField())
universe.protein = Protein('bala1')

minimizer = SteepestDescentMinimizer(universe)
minimizer(steps = 1000)

# constraining all C-H bonds
for bUnit in universe.protein.bondedUnits():
   for bond in bUnit.bonds:
       bondedAtoms = (bond.a1.symbol, bond.a2.symbol)
       if "C" in bondedAtoms and "H" in bondedAtoms:
           dist = universe.distance(bond.a1, bond.a2)
           print  "Constraining:", bond.a1.symbol, bond.a2.symbol, dist
           universe.protein.addDistanceConstraint(bond.a1, bond.a2, dist)

# setting up a dihedral restraint
mybckbn = universe.protein.residues()[1].backbone()
ffrest = HarmonicDihedralRestraint(mybckbn.O,mybckbn.C,mybckbn.C_alpha,mybckbn.N,
ff99 = Amber99ForceField()
ffmy = ff99 + ffrest

universe.thermostat = NoseThermostat(300.*Units.K)

integrator = VelocityVerletIntegrator(universe, delta_t=1.*Units.fs)
trajectory = Trajectory(universe, "bala.nc", "w")

for i in range(20):
	# nan's:
	integrator(steps=6000, actions=[ TrajectoryOutput(trajectory,
("time", "energy", "thermodynamic", "configuration"),0, None, 1),
	# works OK:
	#integrator(steps=6000, actions=[StandardLogOutput(100)])

More information about the mmtk mailing list