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)