openmmforcefields and Interchange¶openmmforcefields is a Python package that provides OpenMM implementations of some small molecule force fields via small molecule template generators.
openmmforcefields provides SMIRNOFF force fields via its infrastructure, internally calling OpenFF software and converting results to something that is compatible with OpenMM's ForceField class. Before doing novel things, let's validate that this implementation provides the same result as directly using OpenFF tools.
The process of preparing inputs is similar; we'll prepare a molecule from a SMILES string and use OpenFF 2.1.0 "Sage" as the force field.
import openmm.app
from openff.toolkit import ForceField, Molecule
from openff.units import unit
from openff.units.openmm import ensure_quantity
from openmmforcefields.generators import SMIRNOFFTemplateGenerator
from openff.interchange import Interchange
from openff.interchange.drivers.openmm import (
_get_openmm_energies,
_process,
get_openmm_energies,
)
molecule = Molecule.from_smiles("O=S(=O)(N)c1c(Cl)cc2c(c1)S(=O)(=O)NCN2")
molecule.generate_conformers(n_conformers=1)
topology = molecule.to_topology()
topology.box_vectors = unit.Quantity([4, 4, 4], unit.nanometer)
smirnoff = SMIRNOFFTemplateGenerator(
molecules=molecule,
forcefield="openff-2.1.0.offxml",
)
forcefield = openmm.app.ForceField()
forcefield.registerTemplateGenerator(smirnoff.generator)
# Sage uses a 9 Angstrom cutoff and 8 Angstrom switching distance
system = forcefield.createSystem(
topology.to_openmm(),
nonbondedCutoff=ensure_quantity(0.9 * unit.nanometer, "openmm"),
switchDistance=ensure_quantity(0.8 * unit.nanometer, "openmm"),
)
openmmforcefields has provided us an (OpenMM) force field with SMIRNOFF parameters included as a template generator. The end goal of this setup is to create an openmm.System, which we can get a single-point energy of.
energy_openmmforcefields = _process(
_get_openmm_energies(
system,
box_vectors=None,
positions=ensure_quantity(molecule.conformers[0], "openmm"),
round_positions=None,
platform="Reference",
),
system=system,
combine_nonbonded_forces=False,
detailed=False,
)
energy_openmmforcefields
We can compare this to an analogous series of operations that use OpenFF tools (the toolkit's ForceField class and Interchange) to create an openmm.System that one would hope has identical contents.
sage = ForceField("openff_unconstrained-2.1.0.offxml")
interchange = Interchange.from_smirnoff(sage, [molecule])
energy_openff = get_openmm_energies(interchange)
energy_openff
Manually inspecting the energies shows zero or marginal differences between them. We can also programmically compare these EnergyReport objects with .compare, which raises an error if there are significant differences, or .diff, which reports differences whether or not they are significant.
In this case, valence energies are exact to nearly machine precision. Non-bonded energies differ slightly due to approximations in how electrostatics are handled with PME.
energy_openff.compare(
energy_openmmforcefields,
{"Nonbonded": 0.002 * unit.kilojoule_per_mole},
)
energy_openff.diff(energy_openmmforcefields)