[MMTK] Molecular Dynamics Question

Ken Warnock kwarnock at mit.edu
Tue Aug 4 15:22:08 UTC 2009


Hi all,

I am fairly new to MD simulations and had a quick question about an NPT
simulation I'm trying to run.  I am currently trying to take a glycan (I've
successfully made the database entry), solvate it, run EM, equilibrate the
system, and then run NPT.  My problem is that during the production run the
Nose Energy term increases rapidly and starts to blow up toward infinity
(It's around 100,000 after 1ns and keeps going up from there).  After
checking the trajectory, it seems like everything in the system gets
"locked" in place and doesn't actually do anything during the simulation.

In checking the examples for NPT simulations, I've noticed that
VelocityScaler and BarostatReset are only used during equilibration, and not
during the production run.  Should I be calling these during the production
run, or is something else going on?  The relevant part of my code is
produced below.  Thank you very much!

Ken

    #Create molecule and atom lists and run EM with all non-H atoms fixed in
place
    atoms = universe.atomList()
    molecules = universe.objectList()
    for atom in atoms:
        if not atom.name.startswith('H'):
            atom.fixed = 1
    minimizer = SteepestDescentMinimizer(universe, threads = 2, actions =
[StandardLogOutput(1000)])
    minimizer(steps = 1000)

    #Release all water molecules and run EM
    if solvated:
        for mol in molecules:
            if mol.name == 'water':
                mol.O.fixed = 0
    minimizer(steps = 10000)

    #Release everything and run EM.
    for atom in atoms:
        atom.fixed = 0
    minimizer(steps = 30000)

    #Gradually heat the system to 310.
    universe.initializeVelocitiesToTemperature(20)
    integrator = VelocityVerletIntegrator(universe, threads = 2, actions =
[VelocityScaler(temperature, 0.1*temperature,0, None, 100),

TranslationRemover(skip = 100),

Heater(20, temperature, .1),

StandardLogOutput(1000)])

    #Add the environment objects and equilibrate the system
    integrator(steps = 10000)
    universe.addObject(barostat)
    universe.addObject(thermostat)
    integrator = VelocityVerletIntegrator(universe, threads = 2, actions =
[VelocityScaler(temperature, 0.1*temperature,0, None, 100),

BarostatReset(skip = 100),

TranslationRemover(skip = 100),

StandardLogOutput(5000)])

    #Run the simulation and save it to a trajectory file.
    integrator(steps = 40000)
    trajectory = Trajectory(universe, outputName + '.nc', 'w')
    integrator = VelocityVerletIntegrator(universe,threads = 2, actions =
[StandardLogOutput(20000),

TrajectoryOutput(trajectory, skip = 1000)])
    integrator(steps = 40000000)
    trajectory.close()
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://starship.python.net/pipermail/mmtk/attachments/20090804/022258e1/attachment.htm>


More information about the mmtk mailing list