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
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],
    target_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.11/site-packages/openff/interchange/components/mdconfig.py:401: 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 7487.617647 26294.876129 7022.396693 -20126.169323 75022.376355 NaN
Amber 7487.616946 26294.876439 7022.396730 314903.626650 75003.279376 NaN
GROMACS 7487.617188 26294.876953 7022.368787 -20133.988037 75020.966675 0.0