Using OpenFF force fields in Amber and GROMACS

The OpenFF Toolkit can create create parametrized openmm.System objects and also, with an extra API call to Interchange, input files for simulations in Amber and GROMACS. This example prepares a system, parametrizes it with Sage, shows how simple it is to write Amber and GROMACS files, and also validates the results.

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 pprint import pprint
from shutil import which

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.visualize()

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.2.0.offxml")
forcefield

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/smirnoff99frosst/smirnoff99frosst.py:11: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  from pkg_resources import resource_filename

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

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 0x7f5e4765f720> >

Preparing an Interchange object

To export to engines other than OpenMM, we will make use of the Interchange project. There is a high-level ForceField.create_interchange method which acts much like ForceField.create_openmm_system. It uses OpenFF Topology and ForceField objects to produce 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 create an openmm.system under the hood; ForceField.create_interchange is an entirely in-memory conversion without calls to external engines or conversion tools.

interchange = forcefield.create_interchange(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_amber(prefix="system")

# Export GROMACS files.
interchange.to_gromacs(prefix="system")

# List the produced input files
!ls system*

system.gro     system.prmtop  system_pointenergy.in
system.inpcrd  system.top     system_pointenergy.mdp

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:401: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but Amber does not implement a switching function. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(

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)
pprint(openmm_energies.energies)

{'Angle': <Quantity(8.66086395, 'kilojoule / mole')>,
 'Bond': <Quantity(0.728703487, 'kilojoule / mole')>,
 'Nonbonded': <Quantity(2.81834152, 'kilojoule / mole')>,
 'Torsion': <Quantity(24.8254772, 'kilojoule / mole')>}

amber_energies = get_amber_energies(interchange)
pprint(amber_energies.energies)

{'Angle': <Quantity(8.66088, 'kilojoule / mole')>,
 'Bond': <Quantity(0.7288528, 'kilojoule / mole')>,
 'Electrostatics': <Quantity(-6.84084, 'kilojoule / mole')>,
 'Torsion': <Quantity(24.8253456, 'kilojoule / mole')>,
 'vdW': <Quantity(9.6667136, 'kilojoule / mole')>}

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:401: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but Amber does not implement a switching function. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(

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:

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

from openff.interchange.drivers import get_gromacs_energies, get_lammps_energies

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

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

{'Angle': <Quantity(8.66089058, 'kilojoule / mole')>,
 'Bond': <Quantity(0.728675842, 'kilojoule / mole')>,
 'Electrostatics': <Quantity(-6.8472681, 'kilojoule / mole')>,
 'RBTorsion': <Quantity(0.0, 'kilojoule / mole')>,
 'Torsion': <Quantity(24.8254776, 'kilojoule / mole')>,
 'vdW': <Quantity(9.66704929, '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)

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:401: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but Amber does not implement a switching function. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:277: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but LAMMPS may not implement a switching function as specified by SMIRNOFF. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(
Bond Angle Torsion Electrostatics vdW RBTorsion
OpenMM 0.728703 8.660864 24.825477 -6.840315 9.658657 NaN
Amber 0.728853 8.660880 24.825346 -6.840840 9.666714 NaN
GROMACS 0.728676 8.660891 24.825478 -6.847268 9.667049 0.0
LAMMPS 0.728703 8.660864 24.825477 -6.944297 9.630999 NaN