{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# Prepare a SMIRNOFF ligand with `openmmforcefields` and Interchange\n",
    "\n",
    "[`openmmforcefields`](https://github.com/openmm/openmmforcefields) is a Python package that provides OpenMM implementations of some small molecule force fields via small molecule template generators."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1",
   "metadata": {},
   "source": [
    "## Validating the implementation of SMIRNOFF force fields\n",
    "\n",
    "`openmmforcefields` provides SMIRNOFF force fields via its infrastructure, internally calling OpenFF software and converting results to something that is compatible with OpenMM's `ForceField` class. Before doing novel things, let's validate that this implementation provides the same result as directly using OpenFF tools.\n",
    "\n",
    "The process of preparing inputs is similar; we'll prepare a molecule from a SMILES string and use OpenFF 2.1.0 \"Sage\" as the force field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openmm.app\n",
    "from openff.toolkit import ForceField, Molecule\n",
    "from openff.units import unit\n",
    "from openff.units.openmm import ensure_quantity\n",
    "from openmmforcefields.generators import SMIRNOFFTemplateGenerator\n",
    "\n",
    "from openff.interchange import Interchange\n",
    "from openff.interchange.drivers.openmm import (\n",
    "    _get_openmm_energies,\n",
    "    _process,\n",
    "    get_openmm_energies,\n",
    ")\n",
    "\n",
    "molecule = Molecule.from_smiles(\"O=S(=O)(N)c1c(Cl)cc2c(c1)S(=O)(=O)NCN2\")\n",
    "molecule.generate_conformers(n_conformers=1)\n",
    "\n",
    "topology = molecule.to_topology()\n",
    "topology.box_vectors = unit.Quantity([4, 4, 4], unit.nanometer)\n",
    "\n",
    "smirnoff = SMIRNOFFTemplateGenerator(\n",
    "    molecules=molecule,\n",
    "    forcefield=\"openff-2.1.0.offxml\",\n",
    ")\n",
    "forcefield = openmm.app.ForceField()\n",
    "forcefield.registerTemplateGenerator(smirnoff.generator)\n",
    "\n",
    "# Sage uses a 9 Angstrom cutoff and 8 Angstrom switching distance\n",
    "system = forcefield.createSystem(\n",
    "    topology.to_openmm(),\n",
    "    nonbondedCutoff=ensure_quantity(0.9 * unit.nanometer, \"openmm\"),\n",
    "    switchDistance=ensure_quantity(0.8 * unit.nanometer, \"openmm\"),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "`openmmforcefields` has provided us an (OpenMM) force field with SMIRNOFF parameters included as a template generator. The end goal of this setup is to create an `openmm.System`, which we can get a single-point energy of."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_openmmforcefields = _process(\n",
    "    _get_openmm_energies(\n",
    "        system,\n",
    "        box_vectors=None,\n",
    "        positions=ensure_quantity(molecule.conformers[0], \"openmm\"),\n",
    "        round_positions=None,\n",
    "        platform=\"Reference\",\n",
    "    ),\n",
    "    system=system,\n",
    "    combine_nonbonded_forces=False,\n",
    "    detailed=False,\n",
    ")\n",
    "\n",
    "energy_openmmforcefields"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "We can compare this to an analogous series of operations that use OpenFF tools (the toolkit's `ForceField` class and Interchange) to create an `openmm.System` that one would hope has identical contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sage = ForceField(\"openff_unconstrained-2.1.0.offxml\")\n",
    "interchange = Interchange.from_smirnoff(sage, [molecule])\n",
    "\n",
    "energy_openff = get_openmm_energies(interchange)\n",
    "\n",
    "energy_openff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "Manually inspecting the energies shows zero or marginal differences between them. We can also programmically compare these `EnergyReport` objects with `.compare`, which raises an error if there are significant differences, or `.diff`, which reports differences whether or not they are significant.\n",
    "\n",
    "In this case, valence energies are exact to nearly machine precision. Non-bonded energies differ slightly due to approximations in how electrostatics are handled with PME."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_openff.compare(\n",
    "    energy_openmmforcefields,\n",
    "    {\"Nonbonded\": 0.002 * unit.kilojoule_per_mole},\n",
    ")\n",
    "\n",
    "energy_openff.diff(energy_openmmforcefields)"
   ]
  }
 ],
 "metadata": {
  "category": "parametrization_evaluation",
  "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.11.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "86c9b142c8dc60dd36d17e2a57efabbd2ed015b9d3db80dd77f3e0894d5aea85"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}