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 pandas import read_csv

# 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 22.024500716 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 13992.528551 207.713905 0.688143
1 20 14574.410217 196.856242 0.688143
2 30 15594.053169 176.121275 0.688143
3 40 16331.516629 161.445346 0.688143
4 50 16026.769608 175.514175 0.688143
... ... ... ... ...
95 960 19839.297982 284.935489 0.688143
96 970 20020.902433 280.410168 0.688143
97 980 19864.318410 282.813750 0.688143
98 990 19728.780797 288.115613 0.688143
99 1000 19821.893229 285.083195 0.688143

100 rows × 4 columns