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)

Example 5-1

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:

\[\begin{split}E(r) = \begin{cases} 0 & r\le2\\ 100(r-2)^2 & r>2 \end{cases}\end{split}\]

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)
...

Example 5-2

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]))

Example 5-3

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 to Context.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())

Example 5-4

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.