Experimental

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

Prepare a GAFF ligand with openmmforcefields and Interchange

openmmforcefields is a Python package that provides OpenMM implementations of some small molecule and biopolymer force fields via template generators. It provides classes that port SMIRNOFF, GAFF, and Espaloma force fields into a format that be used alongside the many force fields that ship with OpenMM by default.

Comparing to GAFF parametrization

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

import openmm.app
from openff.toolkit import Molecule, Quantity, unit
from openmmforcefields.generators import GAFFTemplateGenerator
from rich.pretty import pprint

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 = Quantity([4, 4, 4], unit.nanometer)

gaff = GAFFTemplateGenerator(molecules=molecule)

# If using this alongside another force field, we'd load that in here. But
# in this case, we're only using GAFF for this one molecule, so a "blank"
# force field is a sufficient starting point
forcefield_gaff = openmm.app.ForceField()
forcefield_gaff.registerTemplateGenerator(gaff.generator)

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

%env INTERCHANGE_EXPERIMENTAL=1

env: INTERCHANGE_EXPERIMENTAL=1

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.

imported = Interchange.from_openmm(
    topology=topology.to_openmm(),
    system=system_gaff,
    positions=molecule.conformers[0].to_openmm(),
)

This Interchange contains an internal representation of approximately everything contained in the openmm.System. Topological information (atoms, their elements, bonds, the graph they compose, etc.), forces (including non-bonded settings), and box vectors.

imported.visualize()

We can use Interchange’s machinery to evaluate each of these objects - the openmm.System created by openmmforcefields and the Interchange object created from it.

energy_openmmforcefields = _process(
    _get_openmm_energies(
        system=system_gaff,
        box_vectors=topology.to_openmm().getPeriodicBoxVectors(),
        positions=molecule.conformers[0].to_openmm(),
        round_positions=None,
        platform="Reference",
    ),
    system=system_gaff,
    combine_nonbonded_forces=True,
    detailed=False,
)

pprint(energy_openmmforcefields)

EnergyReport(
│   energies={
│   │   'Bond': <Quantity(47.8549177, 'kilojoule / mole')>,
│   │   'Angle': <Quantity(86.2393125, 'kilojoule / mole')>,
│   │   'Torsion': <Quantity(226.72695, 'kilojoule / mole')>,
│   │   'Nonbonded': <Quantity(-2030.29009, 'kilojoule / mole')>
│   }
)

energy_imported = get_openmm_energies(imported, detailed=False)

pprint(energy_imported)

EnergyReport(
│   energies={
│   │   'Bond': <Quantity(47.8549177, 'kilojoule / mole')>,
│   │   'Angle': <Quantity(86.2393125, 'kilojoule / mole')>,
│   │   'Torsion': <Quantity(226.72695, 'kilojoule / mole')>,
│   │   'Nonbonded': <Quantity(-2030.29252, 'kilojoule / mole')>
│   }
)

Except for some small differences in non-bonded energies due to PME errors, the energies evaluated from each of these objects should match quite closely.

energy_imported.compare(
    energy_openmmforcefields,
    {"Nonbonded": abs(energy_openmmforcefields["Nonbonded"] * 1e-5)},
)

energy_imported.diff(energy_openmmforcefields)

{'Bond': <Quantity(-1.70530257e-13, 'kilojoule / mole')>,
 'Angle': <Quantity(9.9475983e-14, 'kilojoule / mole')>,
 'Torsion': <Quantity(0.0, 'kilojoule / mole')>,
 'Nonbonded': <Quantity(-0.00242279978, 'kilojoule / mole')>}