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.
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 0x1680130d0>
forcefield = ForceField("openff-2.1.0.offxml")
forcefield
<openff.toolkit.typing.engines.smirnoff.forcefield.ForceField at 0x168010810>
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 0x16810b300> >
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.
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")
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.141817584586063 <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')>}
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)
if find_executable("gmx"):
pprint(get_gromacs_energies(interchange).energies)
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 | |
---|---|---|---|---|---|
OpenMM | 0.233879 | 7.141818 | 25.580587 | -6.840315 | 9.658657 |
Amber | 0.233886 | 7.141670 | 25.580558 | -6.840840 | 9.727382 |