examples/using_smirnoff_in_amber_or_gromacs/export_with_interchange.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <div class=\"alert alert-block alert-info\">\n",
    "    <b>Note:</b> \n",
    "    Interchange is in the process of replacing ParmEd in many workflows, but it still in an alpha testing phase. Our internal tests indicate it is reliable for many small-molecule systems, but it is not yet reliable for complex, multi-component systems and there are likely still rough edges throughout. Feedback is welcome on the <a href=https://github.com/openforcefield/openff-interchange/issues/>Interchange issue tracker.</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using OpenFF force fields in Amber and GROMACS\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Preparing an OpenFF Topology\n",
    "\n",
    "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.\n",
    "\n",
    "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](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": [],
   "source": [
    "from openff.toolkit import ForceField, Molecule, Topology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<openff.toolkit.topology.topology.Topology at 0x17248b820>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ethanol = Molecule.from_smiles(\"CCO\")\n",
    "cyclohexane = Molecule.from_smiles(\"C1CCCCC1\")\n",
    "\n",
    "# Load the topology from a PDB file and `Molecule` objects\n",
    "topology = Topology.from_pdb(\n",
    "    \"1_cyclohexane_1_ethanol.pdb\",\n",
    "    unique_molecules=[ethanol, cyclohexane],\n",
    ")\n",
    "\n",
    "topology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing an OpenFF ForceField\n",
    "\n",
    "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](https://github.com/openforcefield/openff-forcefields#how-to-cite).\n",
    "\n",
    "Here we will use force field from the \"Sage\" line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<openff.toolkit.typing.engines.smirnoff.forcefield.ForceField at 0x172176fb0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "forcefield = ForceField(\"openff-2.1.0.offxml\")\n",
    "forcefield"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparing an OpenMM System\n",
    "\n",
    "Once a force field and topology have been loaded, an `openmm.System` can be generated natively with the OpenFF Toolkit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4852506f37f8456cb272686e934e3915",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x1721775d0> >"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "omm_system = forcefield.create_openmm_system(topology)\n",
    "omm_system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Preparing an Interchange object\n",
    "\n",
    "To exports to engines other than OpenMM, we will make use of the [Interchange](https://openff-interchange.readthedocs.io/) 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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Interchange with 7 collections, periodic topology with 27 atoms."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from openff.interchange import Interchange\n",
    "\n",
    "interchange = Interchange.from_smirnoff(\n",
    "    force_field=forcefield,\n",
    "    topology=topology,\n",
    ")\n",
    "interchange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exporting to Amber and GROMACS files\n",
    "\n",
    "Once an `Interchange` object has been constructed, its API can be used to export to files understood by GROMACS, Amber, and more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export AMBER files.\n",
    "interchange.to_prmtop(\"system.prmtop\")\n",
    "interchange.to_inpcrd(\"system.inpcrd\")\n",
    "\n",
    "# Export GROMACS files.\n",
    "interchange.to_top(\"system.top\")\n",
    "interchange.to_gro(\"system.gro\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validating the conversion to Amber files\n",
    "\n",
    "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`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openff.interchange.drivers import get_amber_energies, get_openmm_energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Bond': 0.23387858988917706 <Unit('kilojoule / mole')>,\n",
       " 'Angle': 7.141817584586063 <Unit('kilojoule / mole')>,\n",
       " 'Torsion': 25.58058697014608 <Unit('kilojoule / mole')>,\n",
       " 'Nonbonded': -5.556987691195377 <Unit('kilojoule / mole')>}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "openmm_energies = get_openmm_energies(interchange)\n",
    "openmm_energies.energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Bond': 0.0559 <Unit('kilocalories_per_mole')>,\n",
       " 'Angle': 1.7069 <Unit('kilocalories_per_mole')>,\n",
       " 'Torsion': 6.1139 <Unit('kilocalories_per_mole')>,\n",
       " 'vdW': 9.7273816 <Unit('kilojoule_per_mole')>,\n",
       " 'Electrostatics': -15.216371200000001 <Unit('kilojoule_per_mole')>}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "amber_energies = get_amber_energies(interchange)\n",
    "amber_energies.energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Appendix: Validating the conversion to GROMACS and LAMMPS files\n",
    "\n",
    "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:\n",
    "\n",
    "```conda install \"gromacs >=2021=nompi*\" lammps -c conda-forge```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from distutils.spawn import find_executable\n",
    "from pprint import pprint\n",
    "\n",
    "from openff.interchange.drivers import get_gromacs_energies, get_lammps_energies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "if find_executable(\"lmp_serial\"):\n",
    "    pprint(get_lammps_energies(interchange).energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'Angle': <Quantity(7.14183521, 'kilojoule / mole')>,\n",
      " 'Bond': <Quantity(0.233869165, 'kilojoule / mole')>,\n",
      " 'Electrostatics': <Quantity(-15.2353158, 'kilojoule / mole')>,\n",
      " 'RBTorsion': <Quantity(0.0, 'kilojoule / mole')>,\n",
      " 'Torsion': <Quantity(25.5805855, 'kilojoule / mole')>,\n",
      " 'vdW': <Quantity(9.66705048, 'kilojoule / mole')>}\n"
     ]
    }
   ],
   "source": [
    "if find_executable(\"gmx\"):\n",
    "    pprint(get_gromacs_energies(interchange).energies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Bond</th>\n",
       "      <th>Angle</th>\n",
       "      <th>Torsion</th>\n",
       "      <th>Electrostatics</th>\n",
       "      <th>vdW</th>\n",
       "      <th>RBTorsion</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>OpenMM</th>\n",
       "      <td>0.233879</td>\n",
       "      <td>7.141818</td>\n",
       "      <td>25.580587</td>\n",
       "      <td>-15.215644</td>\n",
       "      <td>9.658657</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Amber</th>\n",
       "      <td>0.233886</td>\n",
       "      <td>7.141670</td>\n",
       "      <td>25.580558</td>\n",
       "      <td>-15.216371</td>\n",
       "      <td>9.727382</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>GROMACS</th>\n",
       "      <td>0.233869</td>\n",
       "      <td>7.141835</td>\n",
       "      <td>25.580585</td>\n",
       "      <td>-15.235316</td>\n",
       "      <td>9.667050</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             Bond     Angle    Torsion  Electrostatics       vdW  RBTorsion\n",
       "OpenMM   0.233879  7.141818  25.580587      -15.215644  9.658657        NaN\n",
       "Amber    0.233886  7.141670  25.580558      -15.216371  9.727382        NaN\n",
       "GROMACS  0.233869  7.141835  25.580585      -15.235316  9.667050        0.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from openff.interchange.drivers.all import get_summary_data\n",
    "\n",
    "get_summary_data(interchange)"
   ]
  },
  {
   "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.10.11"
  },
  "category": "parametrization_evaluation"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}