[MMTK] HarmonicDihedralRestraint problems

[Qadir Timerghazin]

Hello All,

I am trying to use HarmonicDihedralRestraint, but while the distance
and angular restraint work, the dihedral restraints generate
segmentation faults or nan's.
Below is the simplest example demonstrating the problem. There is no
meaningful error generated by python interpreter, so I am not sure how
to debug this...

Thanks in advance!

Qadir Timerghazin


--
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.ForceFields import Amber94ForceField
from MMTK.Trajectory import Trajectory, TrajectoryOutput, \
                            StandardLogOutput
from MMTK.Dynamics import VelocityVerletIntegrator
from MMTK.ForceFields.Restraints import HarmonicDistanceRestraint
from MMTK.ForceFields.Restraints import HarmonicDihedralRestraint
from MMTK.ForceFields.Restraints import HarmonicAngleRestraint


universe = InfiniteUniverse()

universe.molecule = Molecule('water')
universe.molecule = Molecule('water',position=Vector(3.0*Units.Ang,
1.0*Units.Ang, 0.5*Units.Ang))


ffmy = Amber94ForceField()

# Works:
#ffrest = HarmonicDistanceRestraint(universe[0].O, universe[1].O, 5.0*Units.Ang,
#1.0*Units.kcal/Units.mol/Units.Ang**2)

# Works:
#ffrest = HarmonicAngleRestraint(universe[0].H2, universe[0].H1, universe[1].O,
#90.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)

# Potential and kinetic energy are 'nan' after the first step:
ffrest = HarmonicDihedralRestraint(universe[0].H1, universe[0].O, universe[1].O,
universe[1].H1, 90.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)

# Different angle - Segmentation fault:
#ffrest = HarmonicDihedralRestraint(universe[0].H1, universe[0].H2,
universe[1].O,
#universe[1].H1, 0.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)

fftotal = ffmy + ffrest

universe.setForceField(fftotal)

universe.initializeVelocitiesToTemperature(300*Units.K)
integrator = VelocityVerletIntegrator(universe)
trajectory = Trajectory(universe, "water.nc", "w")
integrator(delta_t = 1.*Units.fs, steps = 50,
           actions = [StandardLogOutput(5),
                      TrajectoryOutput(trajectory, ("time", "energy",
                                                    "thermodynamic",
                                                    "configuration"),
                                       0, None, 1)])
trajectory.close()

----------------------------------------------------------------------------------------