[MMTK] another HarmonicDihedralRestraint problem

Qadir Timerghazin qadirt at uwaterloo.ca
Mon Apr 27 22:47:10 UTC 2009


I seem to have another problem with dihedral harmonic restraints -
everything works perfectly if the angle is constrained at a value
sufficiently far from 180. However, if the constraint is around
170-180 degrees then the integration routine goes crazy and eventually
crashes. By the looks of it, this is somehow related to the
constrained angle measurement becoming ambiguous: it switches from
around -180 to +180 and back. This behaviour is fairly consistent for
different systems. Below is an artificial example to demonstrate the

Thanks in advance!


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.ForceFields.Restraints import HarmonicDistanceRestraint
from MMTK.ForceFields.Restraints import HarmonicDihedralRestraint
from MMTK.ForceFields.Restraints import HarmonicAngleRestraint

universe = InfiniteUniverse()
universe.protein = Protein('bala1')
mybckbn = universe.protein.residues()[1].backbone()

# The following works OK - the angle oscillates around 80 deg,
energies are normal
#ffrest = HarmonicDihedralRestraint(mybckbn.O,mybckbn.C,mybckbn.C_alpha,mybckbn.N,

# The following leads to exponentially increasing energies.
# The angle value wildly jumps between ~-180 and ~+180
ffrest = HarmonicDihedralRestraint(mybckbn.O,mybckbn.C,mybckbn.C_alpha,mybckbn.N,

ff99 = Amber99ForceField()
ffmy = ff99 + ffrest
#ffmy = ff99

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

integrator(steps=2000, actions=[ TrajectoryOutput(trajectory, ("time",
"energy", "thermodynamic", "configuration"),0, None, 1),

for i in range(10):
  print "Angle:",

More information about the mmtk mailing list