Simulate an Interchange with LAMMPS

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

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
from lammps import lammps
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 openmm.app import PDBFile
from pandas import read_csv

from openff.interchange import Interchange
from openff.interchange.components.mdconfig import MDConfig

# Read a structure from the Toolkit's test suite into a Topology
pdbfile = PDBFile(
    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"]]
off_topology = Topology.from_openmm(pdbfile.topology, unique_molecules=molecules)

# Construct the Interchange with the OpenFF "Sage" force field
interchange = Interchange.from_smirnoff(
    force_field=ForceField("openff-2.0.0.offxml"),
    topology=off_topology,
)
interchange.positions = pdbfile.positions

Tada! A beautiful solvent system:

interchange.visualize("nglview")

Run a simulation

First, we export a .lmp file that can be read by LAMMPS’ read_data command:

interchange.to_lammps("interchange.lmp")

Now we need to write an input file for LAMMPS. Parts of these input files depend on force field parameters, so we should use a sample input file written for our interchange as a starting point. We can generate such a sample file from MDConfig:

mdconfig = MDConfig.from_interchange(interchange)
mdconfig.write_lammps_input(input_file="auto_generated.in", interchange=interchange)
with open("auto_generated.in") as f:
    print(f.read())

units real
atom_style full

dimension 3
boundary p p p

bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333333333 
pair_style lj/cut/coul/long 9.0 9.0
pair_modify mix arithmetic tail yes

read_data out.lmp

thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe

fix 100 all shake 0.0001 20 10 b 2 4
kspace_style pppm 1e-4
run 0

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/interchange/components/mdconfig.py:275: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but LAMMPS may not implement a switching function as specified by SMIRNOFF. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(

That sample file will only perform a single point energy calculation; here’s a more complete file that includes the above parameters but will run an actual MD simulation. Note that out.lmp is changed to interchange.lmp to reflect the filename we used earlier:

⚠️ Don't use example input files in production
Note that this is still just an example! You should carefully write and inspect input files for all production simulations. 
lammps_in = """ # These commands may not be appropriate for all systems
units real
atom_style full

# PBC in 3 dimensions
dimension 3
boundary p p p

# Bonded interactions in Sage force field
bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
improper_style cvff
special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.8333333333

# Non-bonded interactions in Sage force field
pair_style lj/cut/coul/long 9.0 9.0
pair_modify mix arithmetic tail yes

# Load the parameterized system
read_data interchange.lmp

# Thermostat and velocity generation
fix 3 all nvt temp 300.0 300.0 500
velocity all create 300.0 29348 mom yes rot yes

# Output control
dump traj all dcd 10 traj.dcd
thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe

# PME electrostatics in Sage force field
kspace_style pppm 1e-6

# Run for 1000 steps at 2 fs δt
timestep 2
run 1000

"""

with open("lammps.in", "w") as f:
    f.write(lammps_in)

We’ll use LAMMPS’ Python interface to run the simulation. lmp.file("lammps.in") is equivalent to lammps -in lammps.in on the command line. This will run in serial, which is fine for so few steps and such a small system.

lmp = lammps()

lmp.file("lammps.in")

LAMMPS (7 Feb 2024 - Update 1)

And now we visualize! We can construct an MDTraj Topology from our Interchange by using OpenMM as a lingua franca. LAMMPS produces coordinates that are in the central unit cell, so for a simple system like this we just need to make molecules whole to visualize:

traj = md.load(
    "traj.dcd",
    top=md.Topology.from_openmm(interchange.topology.to_openmm()),
)
traj.make_molecules_whole()
nglview.show_mdtraj(traj)