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 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 openmm.app import PDBFile

from openff.interchange import Interchange

# 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_datafile("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 with the to_lammps_input() method:

interchange.to_lammps_input("interchange_pointenergy.in", data_file="interchange.lmp")
with open("interchange_pointenergy.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 interchange.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.11/site-packages/openff/interchange/components/mdconfig.py:277: 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(

Note that the read_data line must match the data file name - in this case, interchange.lmp. That’s why to_lammps_input needs the data file name. We could alternatively use the to_lammps method to write both files in a consistent way:

interchange.to_lammps("interchange")  # writes interchange.lmp and interchange_pointenergy.in

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.

⚠️ 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 (29 Aug 2024)

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)