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)

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 21.937131583 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 14038.950049 202.573168 0.688143
1 20 14465.589449 197.740204 0.688143
2 30 15486.442825 175.990017 0.688143
3 40 16281.227935 161.032864 0.688143
4 50 15927.886714 174.384249 0.688143
... ... ... ... ...
95 960 19601.068963 287.616641 0.688143
96 970 19643.432005 288.702015 0.688143
97 980 19423.968322 294.503805 0.688143
98 990 19601.283826 292.745468 0.688143
99 1000 19560.004217 294.456342 0.688143

100 rows × 4 columns