Restarting simulations

[Carl Benzer]

Hi,

my first test simulation using prmtop and inpcrd created with AmberTools successfully finished its first 1ns on a GPU. I basically followed the Python tutorial in the OpenMM documentation. Now I am eager to produce longer trajectories ...

With Amber I was used to get a restart file at the end that could be used to continue the simulation. I found a createCheckpoint and loadCheckpoint function in the Python API but I am not sure how to use them. I would really appreciate if someone could explain how to perform simulations in let's say chunks of 500.000 steps, save checkpoint, read checkpoint and continue simulation.

A somewhat similar question is how to do temperature ramps or subsequent simulations with decreasing restraints on atoms in order to slowly start a simulation and adapt the system to the production parameters.

Thanks a lot in advance!

[Carl Benzer]

Mea culpa, I should have done more careful reading last night. Looking through the Wiki and workshop slides, I think I discovered solutions for most of my questions:

Regarding restarting simulations, an example can be found in Lee-Ping's general purpose script for running OpenMM molecular dynamics simulations (https://github.com/leeping/OpenMM-MD).

Regarding temperature ramps for slowly equilibration the system, one can use the following, right?

Code: Select all

simulation.context.setVelocitiesToTemperature(100*kelvin)
simulation.step(100)
simulation.context.setVelocitiesToTemperature(200*kelvin)
simulation.step(100)
...

And finally, positional restraints are mentioned in the 2013 workshop slides on custom forces. Nevertheless I would be grateful if someone could post an example on how to use the original PDB coordinates of the protein to restrain the protein from moving too much during the equilibration of the simulation (first restrain all atoms, then only backbone, then remove restraints).

[Yutong Zhao]

Hi Carl,

Context class's createCheckpoint creates a checkpoint that is only loadable on the exact computer used to create the checkpoint. It basically writes out to a binary stream (or string in case of python) containing various internal states of the system. ie. it tries to maintain the same GPU state (including generated random numbers and what not).

On the other hand, OpenMM 5.1 will have a more portable solution by means of enabling serialization of States that can later used to resume a simulation on any machine. But this resume is not going to be an exact continuation since machine properties can vary.


Regards,
Yutong

[Peter Eastman]

Hi Carl,

What integrator are you using? If it's Verlet, that's a reasonable way of adjusting the temperature. For Langevin, you probably want to adjust the temperature of the integrator's heat bath instead:

Code: Select all

integrator.setTemperature(100*kelvin)

To add restraints to all atoms, the following code should do what you want:

Code: Select all

restraints = CustomExternalForce("k*((x-x0)^2+(y-y0)^2+(z-z0)^2)")
restraints.addGlobalParameter("k", 10*kilojoules/mole/nanometer**2)
restraints.addPerParticleParameter("x0")
restraints.addPerParticleParameter("y0")
restraints.addPerParticleParameter("z0")
for i, pos in enumerate(pdb.positions):
    restraints.addParticle(i, pos)
system.addForce(restraints)

Be sure you insert this code before you create your Simulation. You can select the value of k to set the strength of the restraints. If you want to change it later, for example to slowly turn off the restraints, you can do it like this:

Code: Select all

simulation.context.setParameter("k", 1*kilojoules/mole/nanometer**2)

You can apply restraints to only a subset of atoms by not calling addParticle() for every one. For example, this will omit restraints from water molecules:

Code: Select all

atoms = list(pdb.topology.atoms())
for i, pos in enumerate(pdb.positions):
    if atoms[i].residue.name != 'HOH':
        restraints.addParticle(i, pos)

To apply restraints only to the backbone, the condition should be:

Code: Select all

    if atoms[i].name in ['N', 'CA', 'C']:

As soon as I have time, I'm planning to write a script to automate this, since it's such a common thing to do.

Peter

[Madhuranga Rathnayake]

Hi , a quick question,

I'm using Verlet integrator with Andersen thermostat and MC aniso barostat. In case of temperature ramping I have to reassign temperature in both thermostat and barostat. Also have to add correct velocities on atoms. So I was using following script, but temperature didn't increase. it's still going toward 1K which I initially defined. Any thought on this would be greatly appreciated!

Code: Select all

# Variables
pressure = 1 * u.atmosphere
temperature = 1 * u.kelvin
steps = 500000

# Set input files and force field
pdb = app.PDBFile('ax.pdb')
forcefield = app.ForceField('../../../amoeba-mod.xml')

# Create system for simulation
system = forcefield.createSystem(pdb.topology, polarization='mutual', 
    mutualInducedTargetEpsilon=0.00001, nonbondedMethod=app.PME, 
    nonbondedCutoff=0.4 * u.nanometers, vdwCutoff=0.5*u.nanometer, 
    useDispersionCorrection=True, rigidWater=False, constraints=None,
    ewaldErrorTolerance=0.0005)

# Integrator
integrator = mm.VerletIntegrator(1.0 * u.femtoseconds)
integrator.setConstraintTolerance(0.00001)

# Thermostat
thermostat = mm.AndersenThermostat(temperature, 1.0/u.picoseconds)
system.addForce(thermostat)

# Barostat
barostat = mm.MonteCarloAnisotropicBarostat((pressure,pressure,pressure), temperature)
system.addForce(barostat)

# Set simulation
simulation = app.Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)

simulation.reporters.append(app.StateDataReporter(stdout, 1000, time=True, potentialEnergy=True, temperature=True, separator="\t"))

# Run simulation
for x in range(5):
    tt = (x + 1) *100 * u.kelvin
    thermostat.setDefaultTemperature(tt)
    barostat.setDefaultTemperature(tt)
    simulation.context.setVelocitiesToTemperature(tt)

    print("temperature:", tt)
    simulation.step(10000)

del simulation

Thanks!

Regards,
Madi.

[Peter Eastman]

setDefaultTemperature() just sets the default temperature that will be used for new Contexts. Changing it after you've already created a Context doesn't affect it. This is what you want to do instead:

Code: Select all

simulation.context.setParameter(AndersenThermostat.Temperature(), tt)
simulation.context.setParameter(MonteCarloAnisotropicBarostat.Temperature(), tt)

[Madhuranga Rathnayake]

Thanks Peter, working perfectly now. Appreciate it!!