5. Advanced Simulation Examples¶
In the previous chapter, we looked at some basic scripts for running simulations and saw lots of ways to customize them. If that is all you want to do—run straightforward molecular simulations—you already know everything you need to know. Just use the example scripts and customize them in the ways described in Section 3.6.
OpenMM can do far more than that. Your script has the full OpenMM API at its disposal, along with all the power of the Python language and libraries. In this chapter, we will consider some examples that illustrate more advanced techniques. Remember that these are still only examples; it would be impossible to give an exhaustive list of everything OpenMM can do. Hopefully they will give you a sense of what is possible, and inspire you to experiment further on your own.
Starting in this section, we will assume some knowledge of programming, as well
as familiarity with the OpenMM API. Consult this User’s Guide and the OpenMM
API documentation if you are uncertain about how something works. You can also
usethe Python help
command. For example,
help(Simulation)
will print detailed documentation on the Simulation
class.
5.1. Simulated Annealing¶
Here is a very simple example of how to do simulated annealing. The following lines linearly reduce the temperature from 300 K to 0 K in 100 increments, executing 1000 time steps at each temperature:
...
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
for i in range(100):
integrator.setTemperature(3*(100-i)*kelvin)
simulation.step(1000)
This code needs very little explanation. The loop is executed 100 times. Each
time through, it adjusts the temperature of the LangevinMiddleIntegrator
and then
calls step(1000)
to take 1000 time steps.
5.2. Applying an External Force to Particles: a Spherical Container¶
In this example, we will simulate a non-periodic system contained inside a spherical container with radius 2 nm. We implement the container by applying a harmonic potential to every particle:
where r is the distance of the particle from the origin, measured in nm.
We can easily do this using OpenMM’s CustomExternalForce
class. This class
applies a force to some or all of the particles in the system, where the energy
is an arbitrary function of each particle’s (x, y, z)
coordinates. Here is the code to do it:
...
system = forcefield.createSystem(pdb.topology, nonbondedMethod=CutoffNonPeriodic,
nonbondedCutoff=1*nanometer, constraints=None)
force = CustomExternalForce('100*max(0, r-2)^2; r=sqrt(x*x+y*y+z*z)')
system.addForce(force)
for i in range(system.getNumParticles()):
force.addParticle(i, [])
integrator = LangevinMiddleIntegrator(300*kelvin, 91/picosecond, 0.004*picoseconds)
...
The first thing it does is create a CustomExternalForce
object and add it to the
System
. The argument to CustomExternalForce
is a mathematical expression
specifying the potential energy of each particle. This can be any function of x,
y, and z you want. It also can depend on global or per-particle
parameters. A wide variety of restraints, steering forces, shearing forces,
etc. can be implemented with this method.
Next it must specify which particles to apply the force to. In this case, we
want it to affect every particle in the system, so we loop over them and call
addParticle()
once for each one. The two arguments are the index of
the particle to affect, and the list of per-particle parameter values (an empty
list in this case). If we had per-particle parameters, such as to make the
force stronger for some particles than for others, this is where we would
specify them.
Notice that we do all of this immediately after creating the System
. That is
not an arbitrary choice.
Warning
If you add new forces to a System
, you must do so before creating the Simulation
.
Once you create a Simulation
, modifying the System
will have no effect on that Simulation
.
5.3. Extracting and Reporting Forces (and other data)¶
OpenMM provides reporters for three output formats: PDB, PDBx/mmCIF and DCD. All of those formats store only positions, not velocities, forces, or other data. In this section, we create a new reporter that outputs forces. This illustrates two important things: how to write a reporter, and how to query the simulation for forces or other data.
Here is the definition of the ForceReporter
class:
class ForceReporter(object):
def __init__(self, file, reportInterval):
self._out = open(file, 'w')
self._reportInterval = reportInterval
def __del__(self):
self._out.close()
def describeNextReport(self, simulation):
steps = self._reportInterval - simulation.currentStep%self._reportInterval
return = {'steps': steps, 'periodic': None, 'include':['forces']}
def report(self, simulation, state):
forces = state.getForces().value_in_unit(kilojoules/mole/nanometer)
for f in forces:
self._out.write('%g %g %g\n' % (f[0], f[1], f[2]))
The constructor and destructor are straightforward. The arguments to the constructor are the output filename and the interval (in time steps) at which it should generate reports. It opens the output file for writing and records the reporting interval. The destructor closes the file.
We then have two methods that every reporter must implement:
describeNextReport()
and report()
. A Simulation object
periodically calls describeNextReport()
on each of its reporters to
find out when that reporter will next generate a report, and what information
will be needed to generate it. The return value should be a dictionary with the following content:
steps
(int): The number of time steps until the next report. We calculate this as (report interval)-(current step)%(report interval). For example, if we want a report every 100 steps and the simulation is currently on step 530, we will return 100-(530%100) = 70.include
(list of strings): The types of information that need to be included in the next report. The values in the list correspond to arguments toContext.getState()
. Allowed values include'positions'
,'velocities'
,'forces'
,'energy'
,'parameters'
,'parameterDerivatives'
, and'integratorParameters'
.periodic
(bool, optional): Whether the positions should be wrapped to the periodic box. If None or not set, it will automatically decide whether to wrap positions based on whether the System uses periodic boundary conditions.
When the time comes for the next scheduled report, the Simulation
calls
report()
to generate the report. The arguments are the Simulation
object, and a State
that is guaranteed to contain all the information that was
requested by describeNextReport()
. A State object contains a
snapshot of information about the simulation, such as forces or particle
positions. We call getForces()
to retrieve the forces and convert
them to the units we want to output (kJ/mole/nm). Then we loop over each value
and write it to the file. To keep the example simple, we just print the values
in text format, one line per particle. In a real program, you might choose a
different output format.
Now that we have defined this class, we can use it exactly like any other reporter. For example,
simulation.reporters.append(ForceReporter('forces.txt', 100))
will output forces to a file called “forces.txt” every 100 time steps.
5.4. Computing Energies¶
This example illustrates a different sort of analysis. Instead of running a
simulation, assume we have already identified a set of structures we are
interested in. These structures are saved in a set of PDB files. We want to
loop over all the files in a directory, load them in one at a time, and compute
the potential energy of each one. Assume we have already created our System
and
Simulation
. The following lines perform the analysis:
import os
for file in os.listdir('structures'):
pdb = PDBFile(os.path.join('structures', file))
simulation.context.setPositions(pdb.positions)
state = simulation.context.getState(energy=True)
print(file, state.getPotentialEnergy())
We use Python’s listdir()
function to list all the files in the
directory. We create a PDBFile
object for each one and call
setPositions()
on the Context to specify the particle positions loaded
from the PDB file. We then compute the energy by calling getState()
with the option energy=True
, and print it to the console along
with the name of the file.