[MMTK] constrained bonds - nan problem

Guido Wagner wagner at chemie.uni-frankfurt.de
Wed May 6 08:37:30 UTC 2009


Hi,

I have the same problem (Suse Linux 10.3 (32-bit), MMTK 2.5.24). A 
debugger (gdb) shows that somehow atom-inices get corrupted in the 
C-part of the forcefield, what causes a segfault, when acessing an array 
(eg. the one with the atom koordinates). I already checked the bondlist 
of the universe for wiered indices, but everything seems to be fine. 
 From time to time I also get a double-free-error in the array-module 
(tried Numeric and N from Scientific Python compiled with numpy-support) 
when the temporary gradient-array in the integrator or the minimizer is 
deallocated.

Guido

Qadir Timerghazin schrieb:
> Hi Konrad,
>
> Recently, I've been trying to perform simulations with some bonds constrained,
> but it looks like there is an elusive bug which leads to nan's during
> the integration (or segfaults at some systems, 64-bit Linux mostly). The bug is
> quite subtle - it shows up and disappears depending on the molecular system,
> force field, restraints applied, thermostat and the computer/OS.  When a
> trajectory is written during the integration it always leads to nan's, but at
> different times. The following usually 'helps' to get nan's:
>
> * Dihedral restraints
> * Nose thermostat
> * Repeated calls of the integrator
>
> Finally, I was able to come up with a reproducible example to demonstrate the
> problem (attached). In this example, if a trajectory file is written, I
> consistently get nan's (Mac OS, 32-bit Linux) or segfault (64-bit Linux). If I
> remove the trajectory output as action, it works OK on the most of the computers
> I tested it. (on Mac OS it takes some time before the nan's show up).
>
> 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
>
> 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.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,
> 170.0*Units.deg,30.0*Units.kcal/Units.mol/Units.rad/Units.rad)
> ff99 = Amber99ForceField()
> ffmy = ff99 + ffrest
> universe.setForceField(ffmy)
>
>
> universe.initializeVelocitiesToTemperature(300.*Units.K)
> universe.adjustVelocitiesToConstraints()
> 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),
> StandardLogOutput(100)])
> 	# works OK:
> 	#integrator(steps=6000, actions=[StandardLogOutput(100)])
>
> _______________________________________________
> mmtk maillist  -  mmtk at starship.python.net
> http://starship.python.net/mailman/listinfo/mmtk
>
>   




More information about the mmtk mailing list