Simulate an Interchange with OpenMM

In this example, we’ll quickly construct an Interchange and then run a simulation in OpenMM.

We need an Interchange to get started, so let’s put that together quickly. For more explanation on this process, take a look at the packed_box and protein_ligand examples.

import time

import mdtraj as md
import nglview
import openmm
from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.utils import get_data_file_path
from openff.units import unit
from pandas import read_csv

from openff.interchange import Interchange

# This prepared PDB file from the toolkit's test suite is a box of solvents
pdb_path = get_data_file_path(
    "systems/packmol_boxes/propane_methane_butanol_0.2_0.3_0.5.pdb"
)
molecules = [Molecule.from_smiles(smi) for smi in ["CCC", "C", "CCCCO"]]

# The OpenFF Toolkit can directly read PDB files!
topology = Topology.from_pdb(pdb_path, unique_molecules=molecules)

# Construct the Interchange with the OpenFF "Sage" force field
interchange = ForceField("openff-2.0.0.offxml").create_interchange(topology)
/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/pandas/core/computation/expressions.py:21: UserWarning: Pandas requires version '2.8.4' or newer of 'numexpr' (version '2.7.3' currently installed).
  from pandas.core.computation.check import NUMEXPR_INSTALLED

Tada! A beautiful solvent system:

interchange.visualize()

Run a simulation

We need OpenMM System and Topology objects to run our simulation, as well as positions, so lets prepare them first. We could just reuse the positions from the PDBFile and not have to worry about the units here, but in case you got your positions from somewhere else here’s how to do it in the general case:

Let’s choose parameters for the simulation and use them to prepare an Integrator:

# Length of the simulation.
num_steps = 1000  # number of integration steps to run

# Logging options.
trj_freq = 10  # number of steps per written trajectory frame
data_freq = 10  # number of steps per written simulation statistics

# Integration options
time_step = 2 * openmm.unit.femtoseconds  # simulation timestep
temperature = 300 * openmm.unit.kelvin  # simulation temperature
friction = 1 / openmm.unit.picosecond  # friction constant

integrator = openmm.LangevinIntegrator(temperature, friction, time_step)

Put the parts together, specify initial conditions, and configure how the simulation is recorded:

Configure how the simulation is recorded:

simulation = interchange.to_openmm_simulation(integrator=integrator)

simulation.context.setVelocitiesToTemperature(temperature)

pdb_reporter = openmm.app.PDBReporter("trajectory.pdb", trj_freq)
state_data_reporter = openmm.app.StateDataReporter(
    "data.csv",
    data_freq,
    step=True,
    potentialEnergy=True,
    temperature=True,
    density=True,
)
simulation.reporters.append(pdb_reporter)
simulation.reporters.append(state_data_reporter)

We’re using a PDB reporter for simplicity but you should use something more space-efficient in production. Finally, run it!

print("Starting simulation")
start = time.process_time()

# Run the simulation
simulation.step(num_steps)

end = time.process_time()
print(f"Elapsed time {end - start} seconds")
print("Done!")
Starting simulation
Elapsed time 22.07283057 seconds
Done!

We can take visualize the trajectory with NGLView:

nglview.show_mdtraj(md.load("trajectory.pdb"))

And read the produced data with Pandas:

read_csv("data.csv")
#"Step" Potential Energy (kJ/mole) Temperature (K) Density (g/mL)
0 10 13867.741483 206.746644 0.688143
1 20 14502.333300 196.010638 0.688143
2 30 15613.744083 172.568685 0.688143
3 40 16208.128980 164.406320 0.688143
4 50 15916.101393 176.234882 0.688143
... ... ... ... ...
95 960 19585.437420 278.951918 0.688143
96 970 19425.163372 283.292249 0.688143
97 980 19546.784662 280.827633 0.688143
98 990 19590.100682 278.845496 0.688143
99 1000 19558.247286 279.062159 0.688143

100 rows × 4 columns