Generate and parameterize multi-component systems

The OpenFF Toolkit provides some facilities to prepare topologies from structure files containing multiple molecules, but in other cases external tools are better-suited for the task. In this example, we will use a Python wrapper around the external tool PACKMOL to prepare a system composed of a mixture of organic species.

import nglview
from openff.toolkit import ForceField, Molecule, Topology
from openff.units import unit

from openff.interchange import Interchange
from openff.interchange.components._packmol import UNIT_CUBE, pack_box
from openff.interchange.drivers import get_summary_data

Construct the Topology

The Toolkit provides a class called Topology which stores a collection of molecules. In fact, it can be thought of as a list of Molecule objects. It is similar to an Interchange in that it stores a list of molecules along with positions, box vectors and some other system information, but unlike Interchange a topology is not associated with any force field or parameters.

There are many ways to create a Topology, include:

  • Create one from a (literal) list of molecules via Topology.from_molecules

  • Convert an openmm.app.Topology through Topology.from_openmm

  • Load a multi-component PDB file (modulo some limitations) with Topology.from_pdb

In this example, we’ll use a convenience function provided by openff.interchange that takes a list of Molecule objects and a handful of other options and uses PACKMOL to put everything in a simulation box as a Topology object.

There are many ways to create Molecule objects. Here we’re just using some common organic solvents, so loading them in through SMILES patterns is sufficient.

molecules = [
    Molecule.from_smiles(smi)
    for smi in ["ClC(Cl)(Cl)Cl", "CC1=CC=CC=C1", "CS(=O)C", "CCCCCCO"]
]

The function pack_box takes a number of different arguments (see its docstring for more) covering a number of different use cases. Let’s just pack it with 200 copies of each molecule at density that’s a little lower than what the the mixture might be.

topology = pack_box(
    molecules=molecules,
    number_of_copies=4 * [200],
    mass_density=850 * unit.kilogram / unit.meter**3,
    box_shape=UNIT_CUBE,
)

This Topology object contains copies of each molecule, their positions in the packed box, and the box vectors specifying its periodicity. Let’s have a look!

topology.to_file("system.pdb")
nglview.show_structure_file("system.pdb")

We can get the positions as an array from the PDB file object:

Parametrize with Interchange

Now that we have a topology, we can load a force field and build our Interchange!

force_field = ForceField("openff_unconstrained-2.1.0.offxml")

interchange = Interchange.from_smirnoff(force_field=force_field, topology=topology)

We can visualize it (though, since we can’t see the stored physics parameters, it’ll look the same):

interchange.visualize("nglview")

And we can calculate and compare energies with different MD engines! (The LAMMPS exporter isn’t optimized yet for large systems, so we’re only looking at OpenMM, GROMACS, and Amber.)

get_summary_data(interchange, _engines=["OpenMM", "GROMACS", "Amber"])

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/interchange/components/mdconfig.py:397: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but Amber does not implement a switching function. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(
Bond Angle Torsion Electrostatics vdW RBTorsion
OpenMM 7294.737428 23614.974229 6413.856957 -18759.309105 137408.193219 NaN
Amber 7294.736219 23614.974650 6413.856942 273049.698114 137381.892295 NaN
GROMACS 7294.729492 23614.953125 6413.820932 -18766.899414 137406.373047 0.0