[MMTK] HarmonicDihedralRestraint problems

[Qadir Timerghazin]

Oops, should have CC'ed to the main list....

---------- Forwarded message ----------
Hi Konrad,

Thanks for a prompt reply!

> Your example is not very useful since it tries to define a dihedral
> restraint for a system of three atoms. Two of the atoms in the dihedral
> definition are the same, so you can't expect any meaningful constraint.

In fact, there are _two_ water molecules in my example, so there are 6
atoms. The constrain is defined with:
universe[0].H1, universe[0].O, universe[1].O, universe[1].H1 - note
[0] and [1] indexes.

> Did you check if the problem occurs as well with a well-defined dihedral?

The actual system I am working with is a custom molecule with a custom
FF modification file, which is inconvenient to post. Below is a
protein example, which segfaults with the dihedral constraint.

Thanks a lot!

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
# A simple Molecular Dynamics simulation.
#

import sys
from MMTK import *
from MMTK.Proteins import Protein
from MMTK.ForceFields import Amber99ForceField
from MMTK.Dynamics import VelocityVerletIntegrator, Heater, \
                         TranslationRemover, RotationRemover
from MMTK.Visualization import view
from MMTK.Trajectory import Trajectory, TrajectoryOutput, \
                           RestartTrajectoryOutput, StandardLogOutput, \
                           trajectoryInfo

from MMTK.ForceFields.Restraints import HarmonicDistanceRestraint
from MMTK.ForceFields.Restraints import HarmonicDihedralRestraint
from MMTK.ForceFields.Restraints import HarmonicAngleRestraint


# Define system

universe = InfiniteUniverse()

universe.protein = Protein('bala1')

#print universe.protein.backbone()[0]

#mypsi = universe.protein.residues()[1].phiPsi()[]/Units.deg

mybckbn = universe.protein.residues()[1].backbone()

ffrest = HarmonicDihedralRestraint(mybckbn.O,mybckbn.C,mybckbn.C_alpha,mybckbn.N,
90.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)

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


# Initialize velocities
universe.initializeVelocitiesToTemperature(300.*Units.K)

# Create integrator
integrator = VelocityVerletIntegrator(universe, delta_t=1.*Units.fs)

# Heating and equilibration
integrator(steps=1000,
                   # Heat from 50 K to 300 K applying a temperature
                   # change of 0.5 K/fs; scale velocities at every step.
          actions=[Heater(50.*Units.K, 300.*Units.K, 0.5*Units.K/Units.fs,
                          0, None, 1),
                   # Remove global translation every 50 steps.
                   TranslationRemover(0, None, 50),
                   # Remove global rotation every 50 steps.
                   RotationRemover(0, None, 50),
                   # Log output to screen every 100 steps.
                   StandardLogOutput(100)])
#