Download Notebook View in GitHub Open in Google Colab
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.
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)