Experimental

This notebook is experimental. It may use private or unstable APIs, break suddenly, or produce scientifically unsound results.

Prepare a SMIRNOFF ligand with openmmforcefields and Interchange

openmmforcefields is a Python package that provides OpenMM implementations of some small molecule force fields via small molecule template generators.

Validating the implementation of SMIRNOFF force fields

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 Molecule
from openff.units import unit
from openff.units.openmm import ensure_quantity
from openmmforcefields.generators import SMIRNOFFTemplateGenerator

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.

from openff.units import unit
from openff.units.openmm import ensure_quantity

from openff.interchange.drivers.openmm import _get_openmm_energies, _process

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

EnergyReport(energies={'Bond': <Quantity(35.4195157, 'kilojoule / mole')>, 'Angle': <Quantity(42.5424543, 'kilojoule / mole')>, 'Torsion': <Quantity(136.604497, 'kilojoule / mole')>, 'Electrostatics': <Quantity(-2026.02279, 'kilojoule / mole')>, 'vdW': <Quantity(0.0, 'kilojoule / mole')>})

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.

from openff.toolkit import ForceField

from openff.interchange import Interchange
from openff.interchange.drivers.openmm import get_openmm_energies

sage = ForceField("openff_unconstrained-2.1.0.offxml")
interchange = Interchange.from_smirnoff(sage, [molecule])

energy_openff = get_openmm_energies(interchange)

energy_openff

EnergyReport(energies={'Bond': <Quantity(35.4195157, 'kilojoule / mole')>, 'Angle': <Quantity(42.5424543, 'kilojoule / mole')>, 'Torsion': <Quantity(136.604497, 'kilojoule / mole')>, 'Nonbonded': <Quantity(-2026.02279, 'kilojoule / mole')>})

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)

{'Bond': -3.552713678800501e-14 <Unit('kilojoule / mole')>,
 'Angle': 1.4210854715202004e-14 <Unit('kilojoule / mole')>,
 'Torsion': -2.842170943040401e-14 <Unit('kilojoule / mole')>,
 'Nonbonded': 0.0 <Unit('kilojoule / mole')>}