The Open Forcefield Toolkit can create parametrized System objects that can be natively simulated with OpenMM. This example shows how you can convert a 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 "Parsley" 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 "Parsley" force field.
In [1]:
from simtk.openmm.app import PDBFile
from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField
In [2]:
ethanol = Molecule.from_smiles("CCO")
cyclohexane = Molecule.from_smiles("C1CCCCC1")
# Obtain the OpenMM Topology object from the PDB file.
pdbfile = PDBFile('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 Parsley force field here.
In [3]:
# Load the "Parsley" force field.
forcefield = ForceField('openff_unconstrained-1.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.
In [4]:
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 .
In [5]:
# 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.
In [6]:
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)
In [7]:
from simtk import unit
from simtk.openmm.app import PME, HBonds
from openforcefield.tests.utils import compare_system_parameters, compare_system_energies
# 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)
In [8]:
# 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)