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