[MMTK] constrained bonds - nan problem

Qadir Timerghazin qadirt at uwaterloo.ca
Wed May 6 19:18:31 UTC 2009

> What you should do is heat up your system slowly. You can initialize the
> velocities to a low temperature, say 20 K, and then run NVE MD with a
> heater. There is an example script coming with MMTK that shows how to do it.
> Once your NVE MD is stable at the desired temperature of 300 K (absence of a
> notable energy drift), you can add the thermostat to do NVT simulations. You
> should never add a Nose thermostat to a system that is at a very different
> temperature - this is almost a guarantee for disaster.

I actually do that for the 'real' system I am working on. However, a
proper heating of the system in the example I've given does not seem
to help the nan problem... For now, I am trying ether not to write a
trajectory during the integration (writing snapshots is OK) or do not
cal the integrator repeatedly to avoid the problem.

Thanks a lot!



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

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


hintegrator = VelocityVerletIntegrator(universe, delta_t=1.*Units.fs)
hintegrator(steps=10000, actions=[Heater(20.*Units.K, 300.*Units.K,
                          0, None, 1), StandardLogOutput(100)])

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