Experimental

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

Preparing a 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 the OpenFF Toolkit and using SMIRNOFF force fields provided by its ForceField class. Before doing novel things, let’s validate that this implementation provides the same result as directly using OpenFF tools.

from openff.toolkit import Molecule
from openff.units import unit
from openff.units.openmm import ensure_quantity
from openmm.app import ForceField as OpenMMForceField
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 = OpenMMForceField()
forcefield.registerTemplateGenerator(smirnoff.generator)

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(34.5182169, 'kilojoule / mole')>, 'Angle': <Quantity(185.281577, 'kilojoule / mole')>, 'Torsion': <Quantity(39.0514253, 'kilojoule / mole')>, 'Electrostatics': <Quantity(-2026.02144, '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.0.0.offxml")
interchange = Interchange.from_smirnoff(sage, [molecule])

energy_openff = get_openmm_energies(interchange)

energy_openff

EnergyReport(energies={'Bond': <Quantity(34.5182169, 'kilojoule / mole')>, 'Angle': <Quantity(185.281577, 'kilojoule / mole')>, 'Torsion': <Quantity(39.0514253, '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': -4.973799150320701e-14 <Unit('kilojoule / mole')>,
 'Angle': 1.1368683772161603e-13 <Unit('kilojoule / mole')>,
 'Torsion': 0.0 <Unit('kilojoule / mole')>,
 'Nonbonded': -0.001349857894865636 <Unit('kilojoule / mole')>}

Comparing to GAFF parametrization

openmmforcefields also provides a (partial) implementation of GAFF, using AM1-BCC in place of RESP to generate partial charges. The API for using GAFFTemplateGenerator is analogous to using SMIRNOFFTemplateGenerator.

from openmm.app import PME
from openmmforcefields.generators import GAFFTemplateGenerator

gaff = GAFFTemplateGenerator(molecules=molecule)
forcefield_gaff = OpenMMForceField()
forcefield_gaff.registerTemplateGenerator(gaff.generator)

system_gaff = forcefield_gaff.createSystem(
    topology=topology.to_openmm(),
    nonbondedMethod=PME,
)

Warning: importing 'simtk.openmm' is deprecated.  Import 'openmm' instead.

From here, one could use from_openmm to import this openmm.System into an Interchange object, which then could be combined with other Interchange objects generated by different paths.

from openff.interchange.interop.openmm import from_openmm

from_openmm(
    molecule.to_topology().to_openmm(),
    system_gaff,
    box_vectors=None,
    positions=ensure_quantity(molecule.conformers[0], "openmm"),
)

Interchange with 5 collections, non-periodic topology with 25 atoms.