Changing Temperature and Pressure

Sometimes you want the temperature and/or pressure of a simulation to vary with time. To see how to do this, let’s begin by loading a PDB file and creating a Simulation for it.

[1]:
from openmm.app import *
from openmm import *
from openmm.unit import *

pdb = PDBFile('villin.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
                                 nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)

When we created the LangevinMiddleIntegrator, we specified the temperature as 300*kelvin. We are free to change that at any time. Let’s run a short simulation, steadily decreasing the temperature.

[2]:
for i in range(10):
    integrator.setTemperature((300-30*i)*kelvin)
    integrator.step(10)

Things are a little more complicated when we have a barostat. The Monte Carlo acceptance criterion for the barostat depends on temperature, so we need to tell both the integrator and the barostat about the new temperature.

Let’s add a barostat to the System. Note the call to reinitialize(). Without that, the existing Simulation will not see any changes we make to the System.

[3]:
system.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
simulation.context.reinitialize(preserveState=True)

Because the barostat is part of the System, it stores all its state information in the Context. We can change that information with setParameter().

[4]:
for i in range(10):
    temperature = (300-30*i)*kelvin
    simulation.context.setParameter(MonteCarloBarostat.Temperature(), temperature)
    integrator.setTemperature(temperature)
    integrator.step(10)

Setting the pressure works exactly the same way, just with a different parameter.

[5]:
simulation.context.setParameter(MonteCarloBarostat.Pressure(), 0.5*bar)