Running a PyRETIS simulation with OpenMM

In this example, we show the interface between OpenMM and PyRETIS.

Verification status: engine – see Tutorial map.

Tutorial quick start

  • Best starting point: the OpenMM fixture under examples/tests/test-openmm/ and the files included below.

  • Requirements: openmm importable in Python and an OpenMM Simulation object defined in a Python module.

  • Edit first: the openmm_module and openmm_simulation settings in the [engine] section.

  • Run: pyretis run -i retis.toml -p.

  • Expected output: standard RETIS output folders, with propagation delegated to OpenMM.

  • Related check: examples/tests/test-openmm/; see Example test status.

First, we need a Python module that defines an OpenMM Simulation object. The module openmm_sim.py below defines it for this PDB file containing two water molecules: it builds the System with the Amber and TIP3P force fields, and integrates Langevin dynamics with a CustomIntegrator of the BAOAB splitting that holds the thermal energy in its global variable kT. Write a module for a path-sampling run the same way: build the System and the Simulation as for any OpenMM script, and integrate with a CustomIntegrator of velocity Verlet or BAOAB Langevin dynamics with kT, as this one, or with MTSLangevinIntegrator (see the integrator of a path-sampling run below).

import math

from openmm import app
import openmm as mm
from openmm import unit

pdb = app.PDBFile('input.pdb')
forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml')

system = forcefield.createSystem(pdb.topology, nonbondedMethod=app.PME,
                                 nonbondedCutoff=1.0*unit.nanometers,
                                 constraints=app.HBonds, rigidWater=True,
                                 ewaldErrorTolerance=0.0005)
# Langevin dynamics at 300 K with a friction of 1/ps and a 2 fs step, with
# the BAOAB splitting (Leimkuhler and Matthews, Appl. Math. Res. Express
# 2013, 34-56): the velocities it stores are those at the time of the
# positions, which a path-sampling run needs to reverse and draw velocities
# at a frame. The global variable kT (the thermal energy, kJ/mol) gives
# PyRETIS the temperature. The constraint steps keep the velocities and the
# positions on the constraints of a System that has them.
TEMPERATURE = 300.0 * unit.kelvin
FRICTION = 1.0 / unit.picoseconds
TIMESTEP = 2.0 * unit.femtoseconds
thermal_energy = (unit.MOLAR_GAS_CONSTANT_R * TEMPERATURE).value_in_unit(
    unit.kilojoule_per_mole)
damping = math.exp(-FRICTION.value_in_unit(unit.picoseconds**-1)
                   * TIMESTEP.value_in_unit(unit.picoseconds))

integrator = mm.CustomIntegrator(TIMESTEP)
integrator.addGlobalVariable('kT', thermal_energy)
integrator.addGlobalVariable('a', damping)
integrator.addGlobalVariable('b', math.sqrt(1.0 - damping**2))
integrator.addPerDofVariable('x1', 0.0)
integrator.addUpdateContextState()
# B: half a step of the velocities.
integrator.addComputePerDof('v', 'v + 0.5*dt*f/m')
integrator.addConstrainVelocities()
# A: half a step of the positions.
integrator.addComputePerDof('x', 'x + 0.5*dt*v')
integrator.addComputePerDof('x1', 'x')
integrator.addConstrainPositions()
integrator.addComputePerDof('v', 'v + (x - x1)/(0.5*dt)')
integrator.addConstrainVelocities()
# O: the friction and the noise.
integrator.addComputePerDof('v', 'a*v + b*sqrt(kT/m)*gaussian')
integrator.addConstrainVelocities()
# A: half a step of the positions.
integrator.addComputePerDof('x', 'x + 0.5*dt*v')
integrator.addComputePerDof('x1', 'x')
integrator.addConstrainPositions()
integrator.addComputePerDof('v', 'v + (x - x1)/(0.5*dt)')
integrator.addConstrainVelocities()
# B: half a step of the velocities, with the forces of the new positions.
integrator.addComputePerDof('v', 'v + 0.5*dt*f/m')
integrator.addConstrainVelocities()
integrator.setConstraintTolerance(0.00001)

# Check with "python -m openmm.testInstallation" which platforms are available
# Here I use CUDA, but CPU is more commonly available
# platform = mm.Platform.getPlatformByName('CPU')
platform = mm.Platform.getPlatformByName('CUDA')
simulation = app.Simulation(pdb.topology, system, integrator, platform)
simulation.context.setPositions(pdb.positions)

With this simulation object, we can now construct the required .toml options in the engine section:

Engine
------
type = openmm
class = openmm
openmm_simulation = simulation
openmm_module = openmm_sim.py
subcycles = 10

Here we say that:

  • The engine type is openmm, which tells PyRETIS that it has a special system representation.

  • The engine class is also openmm to tell PyRETIS we want to use the OpenMM engine.

  • The openmm_module is the module from which we want to load the simulation. Here it is openmm_sim.py.

  • The openmm_simulation is the name of the object in that file. In this case it is simulation.

Finally, we give the number of subcycles, which indicates how many MD steps OpenMM performs before PyRETIS requests another frame. This value should be relatively high when using GPUs, as it determines how often PyRETIS communicates with the GPU. However, keep in mind that this also lowers the time resolution of your PyRETIS simulation.

After this, the setup can be run as a regular PyRETIS simulation using pyretis run:

pyretis run -i retis.toml -p

Temperature

The engine runs at the temperature the simulation states: the temperatures that its integrator (the global variable kT of a CustomIntegrator, as in openmm_sim.py above, or each integrator its CompoundIntegrator holds) and the forces of its System with a temperature (a thermostat or a barostat) state, which must all be one temperature. A simulation that states no temperature is refused before any engine runs; to state one without changing the dynamics, add to the System an AndersenThermostat at the temperature of the run with collision frequency 0, which a CustomIntegrator applies in its addUpdateContextState() step. A CustomIntegrator that holds a per-degree-of-freedom variable named kT, or a global one in another letter case (KT, kt, Kt), is refused, and so is one without kT that draws random numbers (uniform or gaussian). See the temperature of an OpenMM simulation in the system section.

Stored frames

PyRETIS keeps the frames of each path it stores in .xyz files beside the path’s traj.txt (see The trajectory stores): the frames OpenMM integrates in <n>_openmm.xyz, one file for each propagation (w<k>_<n>_openmm.xyz for a propagation of worker k), and a frame a move takes over from another path in a file of its own. Each frame holds the positions (nm), the velocities (nm/ps) and, on its box line, the periodic box of that frame (nm), all to nine decimals, so a frame read back – by a restart, or by a run that adds a collective variable – gives the configuration, velocities and box of the MD, also when a barostat changes the box. An initial path of a kick initiation is stored the same way, its frames in its traj.xyz and their potential and kinetic energies in its energy.txt, also with [initial-path] kick-parallel = true. A frame staged for a load without a box line, or with a box line of zeros, is read in the periodic box of the OpenMM simulation. The frames of a load folder without energy.txt get their energies from OpenMM: the potential energy at the positions of the frame in the box it is read in, and the kinetic energy of its velocities, in kJ/mol. A frame whose box line has a length of inf or 0 stops the run with a ValueError when it is set in the OpenMM context: OpenMM propagates a periodic cell only.

Random numbers and reproducibility

Each propagation of a move, and each crossing search of a kick, starts from seeds drawn from the random stream of the move or of the kick, which the seeds of the input fix. The parts of the OpenMM simulation that draw random numbers take them: a CustomIntegrator whose expressions use uniform or gaussian (the BAOAB integrator of openmm_sim.py above, MTSLangevinIntegrator), also in a CompoundIntegrator, and an AndersenThermostat share one seed, and each Monte Carlo barostat gets a seed of its own. OpenMM seeds that shared stream with the seed of an integrator, also one that draws no random numbers, so when a part draws, every integrator (a velocity-Verlet CustomIntegrator next to an AndersenThermostat, each integrator of a CompoundIntegrator) takes the shared seed too. OpenMM reads these seeds when it builds a context, so PyRETIS rebuilds the context (Context.reinitialize) for each propagation. A seed the module sets, as integrator.setRandomNumberSeed(42), does not change the paths of a path-sampling run.

The time step of the engine is the step size of the integrator. A VariableVerletIntegrator or a VariableLangevinIntegrator, alone or in a CompoundIntegrator, has no fixed step size, so the times of the paths and the rate are undefined: PyRETIS refuses it before the run.

A path-sampling run needs an integrator whose velocities are at the time of its positions: it reverses the velocities of a frame, or draws new ones, and integrates on from the positions of the frame. The built-in integrators of OpenMM store the velocities half a time step away from the positions (VerletIntegrator, LangevinIntegrator, LangevinMiddleIntegrator, NoseHooverIntegrator, which integrates with a leap-frog scheme too) or have none (BrownianIntegrator), and a path-sampling run refuses them before it starts, alone or in a CompoundIntegrator, together with the other built-in integrators and the accelerated-MD CustomIntegrator classes (AMDIntegrator, AMDForceGroupIntegrator and DualAMDIntegrator), which PyRETIS knows by their global variables and steps, also in a CompoundIntegrator or read from an XML file. Use a CustomIntegrator of velocity Verlet or BAOAB Langevin dynamics with the global variable kT, as the openmm_sim.py above, which integrates Langevin dynamics with the BAOAB splitting, or OpenMM’s MTSLangevinIntegrator, which holds kT. For Brownian dynamics in path sampling, use the internal Langevin engine with high_friction = true, whose path sampling the method validation of PyRETIS covers.

PyRETIS takes the temperature of its velocity draws from the integrator and from the forces of the System with a temperature: from the global variable kT, the thermal energy in kJ/mol, of a CustomIntegrator, divided by the molar gas constant MOLAR_GAS_CONSTANT_R of OpenMM, and also from getTemperature for the integrator of the simulation when its class has it, as those of openmmtools (a temperature Quantity; OpenMM gives a CustomIntegrator a CompoundIntegrator holds as its own CustomIntegrator, without its Python class, whose kT is read), and from the temperature of an AndersenThermostat or a barostat. Compute kT as MOLAR_GAS_CONSTANT_R times the temperature: a kT computed with another value of the gas constant reads as another temperature, which a thermostat at the temperature then differs from. A kT computed from a temperature divides back to it within 1.511 machine epsilons relative to it (3 units in the last place), as a Quantity converted to kJ/mol as in openmm_sim.py above and in MTSLangevinIntegrator (250 K gives 250.00000000000003 K; OpenMM 8.5, the temperatures k/64 K and k/10 K from 1 K to 2000 K). A CompoundIntegrator gives the temperatures of its integrators, that of the integrator it steps with first. These temperatures must be one temperature: equal to four machine epsilons relative to the larger (8.9e-16, for 250 K about 2e-13 K), and the velocities PyRETIS draws take the first of them. When [system] temperature or the temperature of the engine section differs from the temperature of a kT the engine runs at, the message says how PyRETIS reads kT. A CustomIntegrator without kT that draws no random numbers (velocity Verlet, MTSIntegrator) gives no temperature: the simulation then runs at the temperature the forces with a temperature (a thermostat or a barostat, as an AndersenThermostat or a MonteCarloBarostat), or the other integrators of a CompoundIntegrator, give, and a simulation that gives none is refused before the run starts. A CustomIntegrator without kT that draws random numbers is refused before the run starts, also in a CompoundIntegrator: its random numbers act at a temperature PyRETIS does not read; hold its thermal energy in kT. OpenMM runs the path-sampling tasks; an md, md-nve or md-flux task with OpenMM is refused before it starts (run it with an internal engine, GROMACS, CP2K or LAMMPS).

Each propagation then starts from the context parameters and the integrator parameters the simulation holds when the module gives it to PyRETIS: the global and per-degree-of-freedom variables of a CustomIntegrator start each propagation as they were then, and do not carry over from the propagation before. Each propagation also starts at step count 0, so a CMMotionRemover with a frequency above 1 acts at the same steps of each one. Before it starts, a path-sampling run refuses a CustomIntegrator whose step reads a variable before it sets it (a variable it carries from one step to the next, as the friction of a Nose-Hoover thermostat; the message names it) when its steps draw no random numbers: no uniform or gaussian, and no AndersenThermostat or Monte Carlo barostat acting in an addUpdateContextState() step. It warns, once per engine section when it is set up, for such a CustomIntegrator whose steps draw random numbers, for one whose steps set kT (a temperature ramp), for a Monte Carlo barostat with a frequency above 0, whose volume moves are counted from the start of each propagation, and for a CMMotionRemover with a frequency above 1. An AndersenThermostat, a CMMotionRemover and a Monte Carlo barostat act in the addUpdateContextState() step of a CustomIntegrator, and a path-sampling run refuses them, before it starts, in a simulation whose CustomIntegrator has no such step: call integrator.addUpdateContextState() at the start of its steps, as a velocity-Verlet integrator with an AndersenThermostat needs (see the thermostats and barostats of a path-sampling run). The CustomIntegrator of openmm_sim.py above has that step, and its CMMotionRemover, which createSystem adds, has the frequency 1.

A move then gives the same paths whichever worker runs it, with [runner] workers above 1 and with [initial-path] kick-parallel = true. The numbers are reproduced on the same OpenMM build, platform and platform settings only: the CPU platform draws the random numbers of each thread from a stream of its own, so they depend on the number of threads (OPENMM_CPU_THREADS). The single-threaded Reference platform is the one examples/tests/test-openmm/ pins. Rebuilding the context costs time for each propagation, more on a GPU platform, where building a context takes longest.

Tested by

OpenMM coverage is tracked through examples/tests/test-openmm/. Developer runs are optional because they depend on a local OpenMM installation; tagged releases require the recorded 8.5.1 source revision and the strict golden comparison. See Example test status for the current run status before treating this example as verified on a given platform.