[MMTK] another HarmonicDihedralRestraint problem

[Qadir Timerghazin]

Hello,

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
problem.

Thanks in advance!

Cheers,
Qadir

--
Dr. Qadir K. Timerghazin,
NSERC Postdoctoral Fellow

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


#!/usr/bin/python

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,
#80.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)

# 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,
170.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)


ff99 = Amber99ForceField()
ffmy = ff99 + ffrest
#ffmy = ff99
universe.setForceField(ffmy)

universe.initializeVelocitiesToTemperature(300.*Units.K)
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),
StandardLogOutput(100)])

for i in range(10):
  integrator(steps=50)
  print "Angle:",
universe.dihedral(mybckbn.O,mybckbn.C,mybckbn.C_alpha,mybckbn.N)/Units.deg