{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Note: \n", " Please take a look at the other notebook in the folder (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.
" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Convert Open Forcefield System to AMBER and GROMACS input files\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create an OpenMM System\n", "\n", "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.\n", "\n", "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](https://open-forcefield-toolkit.readthedocs.io/en/latest/faq.html#the-partial-charges-generated-by-the-toolkit-don-t-seem-to-depend-on-the-molecule-s-conformation-is-this-a-bug) for more information." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "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\n" ] } ], "source": [ "try:\n", " from openmm.app import PDBFile\n", "except ImportError:\n", " from simtk.openmm.app import PDBFile\n", "\n", "from openff.toolkit.topology import Molecule, Topology\n", "from openff.toolkit.typing.engines.smirnoff import ForceField\n", "from openff.toolkit.utils import get_data_file_path" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "ethanol = Molecule.from_smiles(\"CCO\")\n", "cyclohexane = Molecule.from_smiles(\"C1CCCCC1\")\n", "\n", "# Obtain the OpenMM Topology object from the PDB file.\n", "pdbfile = PDBFile(\n", " get_data_file_path(\"systems/test_systems/1_cyclohexane_1_ethanol.pdb\")\n", ")\n", "omm_topology = pdbfile.topology\n", "\n", "# Create the Open Forcefield Topology.\n", "off_topology = Topology.from_openmm(\n", " omm_topology, unique_molecules=[ethanol, cyclohexane]\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Load the \"Sage\" force field.\n", "forcefield = ForceField(\"openff_unconstrained-2.0.0.offxml\")\n", "omm_system = forcefield.create_openmm_system(off_topology)" ] }, { "cell_type": "markdown", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "### Convert OpenMM System to AMBER and GROMACS files\n", "\n", "\n", "First, we convert the OpenMM `System` into a ParmEd `Structure`. We'll use the atom positions in the PDB to create the coordinate files.\n", "\n", "
\n", " Warning: ParmEd's Structure model is inspired by AMBER, and some information in an OpenMM System are not directly translatable into a Structure. In particular, as of today (4/2/2019), long-range interaction treatment method (e.g., PME, CutoffPeriodic) and parameters (e.g., cutoff and cutoff switching distance, PME error tolerance) are known to be lost during the conversion.\n", "
" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import parmed\n", "\n", "# Convert OpenMM System to a ParmEd structure.\n", "parmed_structure = parmed.openmm.load_topology(\n", " omm_topology, omm_system, pdbfile.positions\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 ." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Export AMBER files.\n", "parmed_structure.save(\"system.prmtop\", overwrite=True)\n", "parmed_structure.save(\"system.inpcrd\", overwrite=True)\n", "\n", "# Export GROMACS files.\n", "parmed_structure.save(\"system.top\", overwrite=True)\n", "parmed_structure.save(\"system.gro\", overwrite=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Validate the conversion\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9 nm\n", "False\n", "True\n" ] } ], "source": [ "try:\n", " import openmm\n", "except ImportError:\n", " from simtk import openmm\n", "\n", "for force in omm_system.getForces():\n", " if isinstance(force, openmm.NonbondedForce):\n", " break\n", "print(force.getCutoffDistance())\n", "print(force.getUseSwitchingFunction())\n", "print(force.getNonbondedMethod() == openmm.NonbondedForce.PME)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "try:\n", " from openmm import unit\n", " from openmm.app import PME, HBonds\n", "except ImportError:\n", " from simtk import unit\n", " from simtk.openmm.app import PME, HBonds\n", "\n", "from openff.toolkit.tests.utils import (\n", " compare_system_energies,\n", " compare_system_parameters,\n", ")\n", "\n", "# Load the prmtop/inpcrd files into a ParmEd Structure.as an OpenMM System object.\n", "amber_structure = parmed.load_file(\"system.prmtop\", \"system.inpcrd\")\n", "\n", "# Convert the Structure to an OpenMM System. Note that by\n", "# default ParmEd will add a CMMotionRemover force to the\n", "# System, and won't constrain the hydrogen bonds.\n", "amber_system = amber_structure.createSystem(\n", " nonbondedMethod=PME,\n", " nonbondedCutoff=9.0 * unit.angstrom,\n", " switchDistance=0.0 * unit.angstrom,\n", " constraints=HBonds,\n", " removeCMMotion=False,\n", ")\n", "\n", "# Compare the parameters of the original and converted Systems.\n", "# This raises FailedParameterComparisonError if the comparison fails.\n", "compare_system_parameters(omm_system, amber_system)\n", "\n", "# Compare the energies by force.\n", "# This raises FailedEnergyComparisonError if the comparison fails.\n", "amber_energies, omm_energies = compare_system_energies(\n", " amber_system,\n", " omm_system,\n", " amber_structure.positions,\n", " amber_structure.box_vectors,\n", " rtol=1e-3,\n", ")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "System loaded from AMBER files:\n", "-------------------------------\n", "{'HarmonicAngleForce': Quantity(value=155.70130920410156, unit=kilojoule/mole),\n", " 'HarmonicBondForce': Quantity(value=0.4403935372829437, unit=kilojoule/mole),\n", " 'NonbondedForce': Quantity(value=2.8259036956507657, unit=kilojoule/mole),\n", " 'PeriodicTorsionForce': Quantity(value=19.563232421875, unit=kilojoule/mole)}\n", "\n", "Original OpenMM System:\n", "-----------------------\n", "{'HarmonicAngleForce': Quantity(value=155.70130920410156, unit=kilojoule/mole),\n", " 'HarmonicBondForce': Quantity(value=0.44038787484169006, unit=kilojoule/mole),\n", " 'NonbondedForce': Quantity(value=2.8251933539722245, unit=kilojoule/mole),\n", " 'PeriodicTorsionForce': Quantity(value=19.563230514526367, unit=kilojoule/mole)}\n" ] } ], "source": [ "# Pretty-print the energies by component.\n", "from pprint import pprint\n", "\n", "print(\"System loaded from AMBER files:\")\n", "print(\"-------------------------------\")\n", "pprint(amber_energies)\n", "\n", "print(\"\\nOriginal OpenMM System:\")\n", "print(\"-----------------------\")\n", "pprint(omm_energies)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.12" } }, "nbformat": 4, "nbformat_minor": 4 }