export_with_interchange.ipynb
), showcasing how to use OpenFF Interchange. Interchange is well-validated for small molecule interoperability with several formats and is in the process of replacing ParmEd in OpenFF workflows.The Open Forcefield Toolkit can create parametrized System
objects that can be natively simulated with OpenMM. This example shows how you can convert such an OpenMM System
into AMBER prmtop/inpcrd and GROMACS top/gro input files through the ParmEd library.
We start by loading a PDB file containing one copy of ethanol and cyclohexane. Our goal is to create an OFF 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 partial charges and stereochemistry that is not included in an OpenMM topology. 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 Open Force Field Toolkit produces deterministic partial charges that do not depend on the input conformation of parameterized molecules. See the FAQ for more information.
try:
from openmm.app import PDBFile
except ImportError:
from simtk.openmm.app import PDBFile
from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.utils import get_data_file_path
Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing
ethanol = Molecule.from_smiles("CCO")
cyclohexane = Molecule.from_smiles("C1CCCCC1")
# Obtain the OpenMM Topology object from the PDB file.
pdbfile = PDBFile(
get_data_file_path("systems/test_systems/1_cyclohexane_1_ethanol.pdb")
)
omm_topology = pdbfile.topology
# Create the Open Forcefield Topology.
off_topology = Topology.from_openmm(
omm_topology, unique_molecules=[ethanol, cyclohexane]
)
Now we parametrize the OFF Topology
to create an OpenMM System
. Since ParmEd will run with the constraints=HBonds
keyword later, we use the unconstrained version of the Sage force field here.
# Load the "Sage" force field.
forcefield = ForceField("openff_unconstrained-2.0.0.offxml")
omm_system = forcefield.create_openmm_system(off_topology)
First, we convert the OpenMM System
into a ParmEd Structure
. We'll use the atom positions in the PDB to create the coordinate files.
import parmed
# Convert OpenMM System to a ParmEd structure.
parmed_structure = parmed.openmm.load_topology(
omm_topology, omm_system, pdbfile.positions
)
We can then use ParmEd to convert an OpenMM System
to prmtop/inpcrd or top/gro files that can be read by AMBER and GROMACS respectively. ParmEd is capable of converting parametrized files to other formats as well. For further information, see ParmEd's documentation: https://parmed.github.io/ParmEd/html/readwrite.html .
# Export AMBER files.
parmed_structure.save("system.prmtop", overwrite=True)
parmed_structure.save("system.inpcrd", overwrite=True)
# Export GROMACS files.
parmed_structure.save("system.top", overwrite=True)
parmed_structure.save("system.gro", overwrite=True)
ParmEd is generally a reliable and robust library, but we can easily check that everything went as expected during the conversion by loading the exported files into an OpenMM System
and comparing it with the original. Note that you'll have to specify the correct nonbonded method and cutoff settings for the energy comparison to make sense since they are not included in the AMBER prmtop (or GROMACS top/gro) files.
try:
import openmm
except ImportError:
from simtk import openmm
for force in omm_system.getForces():
if isinstance(force, openmm.NonbondedForce):
break
print(force.getCutoffDistance())
print(force.getUseSwitchingFunction())
print(force.getNonbondedMethod() == openmm.NonbondedForce.PME)
0.9 nm False True
try:
from openmm import unit
from openmm.app import PME, HBonds
except ImportError:
from simtk import unit
from simtk.openmm.app import PME, HBonds
from openff.toolkit.tests.utils import (
compare_system_energies,
compare_system_parameters,
)
# Load the prmtop/inpcrd files into a ParmEd Structure.as an OpenMM System object.
amber_structure = parmed.load_file("system.prmtop", "system.inpcrd")
# Convert the Structure to an OpenMM System. Note that by
# default ParmEd will add a CMMotionRemover force to the
# System, and won't constrain the hydrogen bonds.
amber_system = amber_structure.createSystem(
nonbondedMethod=PME,
nonbondedCutoff=9.0 * unit.angstrom,
switchDistance=0.0 * unit.angstrom,
constraints=HBonds,
removeCMMotion=False,
)
# Compare the parameters of the original and converted Systems.
# This raises FailedParameterComparisonError if the comparison fails.
compare_system_parameters(omm_system, amber_system)
# Compare the energies by force.
# This raises FailedEnergyComparisonError if the comparison fails.
amber_energies, omm_energies = compare_system_energies(
amber_system,
omm_system,
amber_structure.positions,
amber_structure.box_vectors,
rtol=1e-3,
)
# Pretty-print the energies by component.
from pprint import pprint
print("System loaded from AMBER files:")
print("-------------------------------")
pprint(amber_energies)
print("\nOriginal OpenMM System:")
print("-----------------------")
pprint(omm_energies)
System loaded from AMBER files: ------------------------------- {'HarmonicAngleForce': Quantity(value=155.70130920410156, unit=kilojoule/mole), 'HarmonicBondForce': Quantity(value=0.4403935372829437, unit=kilojoule/mole), 'NonbondedForce': Quantity(value=2.8259036956507657, unit=kilojoule/mole), 'PeriodicTorsionForce': Quantity(value=19.563232421875, unit=kilojoule/mole)} Original OpenMM System: ----------------------- {'HarmonicAngleForce': Quantity(value=155.70130920410156, unit=kilojoule/mole), 'HarmonicBondForce': Quantity(value=0.44038787484169006, unit=kilojoule/mole), 'NonbondedForce': Quantity(value=2.8251933539722245, unit=kilojoule/mole), 'PeriodicTorsionForce': Quantity(value=19.563230514526367, unit=kilojoule/mole)}