Note: Interchange is in the process of replacing ParmEd in many workflows, but it still in an alpha testing phase. Our internal tests indicate it is reliable for many small-molecule systems, but it is not yet reliable for complex, multi-component systems and there are likely still rough edges throughout. Feedback is welcome on the Interchange issue tracker.

Using OpenFF force fields in Amber and GROMACS

The Open Forcefield Toolkit can create parametrized openmm.System objects that can be natively simulated with OpenMM. This example shows the Interchange project can enable parallel workflows using Amber and GROMACS.

Preparing an OpenFF Topology

We start by loading a PDB file containing one copy of ethanol and cyclohexane. Our goal is to create an OpenFF Topology object describing this system that we can parametrize with the SMIRNOFF-format “Sage” force field.

The two Molecule objects created from the SMILES strings can contain information such as formal charges and stereochemistry that is not included in a PDB file. These objects are passed to the PDB loading function via the unique_molecules argument. In this example, partial charges are not explicitly given, and ForceField will assign AM1/BCC charges as specified by the “Sage” force field. Note that the OpenFF Toolkit produces partial charges that do not depend on the input conformation of parameterized molecules. See the FAQ for more information.

from openff.toolkit import ForceField, Molecule, Topology

ethanol = Molecule.from_smiles("CCO")
cyclohexane = Molecule.from_smiles("C1CCCCC1")

# Load the topology from a PDB file and `Molecule` objects
topology = Topology.from_pdb(
    "1_cyclohexane_1_ethanol.pdb",
    unique_molecules=[ethanol, cyclohexane],
)

topology

<openff.toolkit.topology.topology.Topology at 0x7fb3eb8c7a30>

Preparing an OpenFF ForceField

Once the ForceField class is imported, the only decision to make is which force field to use. An exhaustive list of force fields released by the Open Force Field Initiative can be found here.

Here we will use force field from the “Sage” line.

forcefield = ForceField("openff-2.1.0.offxml")
forcefield

<openff.toolkit.typing.engines.smirnoff.forcefield.ForceField at 0x7fb3ebbab3d0>

Preparing an OpenMM System

Once a force field and topology have been loaded, an openmm.System can be generated natively with the OpenFF Toolkit.

omm_system = forcefield.create_openmm_system(topology)
omm_system

<openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x7fb3dfffae50> >

Preparing an Interchange object

To exports to engines other than OpenMM, we will make use of the Interchange project. There is a high-level Interchange.from_smirnoff function that consumes OpenFF Toolkit and ForceField objects and produces an Interchange object which can then be exported to formats understood by other molecular simulation engines. This extra step is needed to provide a clean interface between applied parameters and engines. Note also that this step does not require an OpenMM System to be generated; ForceField.create_openmm_system does not need to be called to use Amber and GROMACS.

from openff.interchange import Interchange

interchange = Interchange.from_smirnoff(
    force_field=forcefield,
    topology=topology,
)
interchange

Interchange with 7 collections, periodic topology with 27 atoms.

Exporting to Amber and GROMACS files

Once an Interchange object has been constructed, its API can be used to export to files understood by GROMACS, Amber, and more.

# Export AMBER files.
interchange.to_prmtop("system.prmtop")
interchange.to_inpcrd("system.inpcrd")

# Export GROMACS files.
interchange.to_top("system.top")
interchange.to_gro("system.gro")

Validating the conversion to Amber files

The Interchange project includes functions that take in an Interchange object and call out to simulation engines to run single-point energy calculations (with no minimization or dynamics) for the purpose of validating the export layer with each engine. Under the hood, each of these functions calls API points like those used above while converting to files understood by each engine. These rely on having each engine installed and accessible in the current $PATH.

from openff.interchange.drivers import get_amber_energies, get_openmm_energies

openmm_energies = get_openmm_energies(interchange)
openmm_energies.energies

{'Bond': 0.23387858988917706 <Unit('kilojoule / mole')>,
 'Angle': 7.141817584586029 <Unit('kilojoule / mole')>,
 'Torsion': 25.58058697014608 <Unit('kilojoule / mole')>,
 'Nonbonded': 2.818341523387724 <Unit('kilojoule / mole')>}

amber_energies = get_amber_energies(interchange)
amber_energies.energies

{'Bond': 0.0559 <Unit('kilocalories_per_mole')>,
 'Angle': 1.7069 <Unit('kilocalories_per_mole')>,
 'Torsion': 6.1139 <Unit('kilocalories_per_mole')>,
 'vdW': 9.7273816 <Unit('kilojoule_per_mole')>,
 'Electrostatics': -6.84084 <Unit('kilojoule_per_mole')>}

Appendix: Validating the conversion to GROMACS and LAMMPS files

If GROMACS and/or LAMMPS are installed on your machine, the same comparisons can also take place with those engines. They are available via conda by running a command like:

conda install "gromacs >=2021=nompi*" lammps -c conda-forge

from distutils.spawn import find_executable
from pprint import pprint

from openff.interchange.drivers import get_gromacs_energies, get_lammps_energies

if find_executable("lmp_serial"):
    pprint(get_lammps_energies(interchange).energies)

{'Angle': <Quantity(1.7069354, 'kilocalorie_per_mole')>,
 'Bond': <Quantity(0.33085445, 'kilocalorie_per_mole')>,
 'Electrostatics': <Quantity(-1.659727, 'kilocalorie_per_mole')>,
 'Torsion': <Quantity(6.113907, 'kilocalorie_per_mole')>,
 'vdW': <Quantity(2.30186395, 'kilocalorie_per_mole')>}

if find_executable("gmx"):
    pprint(get_gromacs_energies(interchange).energies)

{'Angle': <Quantity(7.14183187, 'kilojoule / mole')>,
 'Bond': <Quantity(0.233864918, 'kilojoule / mole')>,
 'Electrostatics': <Quantity(-6.85864258, 'kilojoule / mole')>,
 'RBTorsion': <Quantity(0.0, 'kilojoule / mole')>,
 'Torsion': <Quantity(25.5805855, 'kilojoule / mole')>,
 'vdW': <Quantity(9.66705143, 'kilojoule / mole')>}

Finally, there is a helper function get_summary_data that will attempt to run drivers of each engine. A summary reported is prepared as a Pandas DataFrame.

from openff.interchange.drivers.all import get_summary_data

get_summary_data(interchange)
Bond Angle Torsion Electrostatics vdW RBTorsion
OpenMM 0.233879 7.141818 25.580587 -6.840315 9.658657 NaN
Amber 0.233886 7.141670 25.580558 -6.840840 9.727382 NaN
GROMACS 0.233865 7.141832 25.580585 -6.858642 9.667051 0.0
LAMMPS 1.384295 7.141818 25.580587 -6.944298 9.630999 NaN