07_ligands/01_ibuprofen_gas_phase.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ibuprofen in gas phase\n",
    "\n",
    "This notebook demonstrates a simple use case of the [Open Force Field (OFF)](https://openforcefield.org/) toolkit. [OFF](https://openforcefield.org/) is a relatively new initiative and aims to bring more transparency into the process of generating force fields (for small molecules) through open science and open source principles.\n",
    "\n",
    "We will create a force field for a small molecule, for which no standard force fields are available, using version 2.0.0 of the [Sage force field](https://openforcefield.org/force-fields/force-fields/). Sage covers a very broad set of organic small drug-like molecules. This parameterization is comparable to older initiatives like GAFF, OPLS or CGENFF. Unlike these older models, Sage makes extensive use of modern cheminformatics tools, which makes the force field easier to extend and improve. More details can be found here: [10.1021/acs.jctc.8b00640](https://pubs.acs.org/doi/10.1021/acs.jctc.8b00640)\n",
    "\n",
    "This notebook uses ibuprofen as a simple test case. In principle any small drug-like molecule could be used instead and your are encouraged to try your favorite compound from [PubChem](https://pubchem.ncbi.nlm.nih.gov/).\n",
    "\n",
    "The imports below will generate some warnings, which you can safely ignore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python built-in modules\n",
    "from sys import stdout\n",
    "\n",
    "# Popular scientific packages for Python\n",
    "import mdtraj\n",
    "\n",
    "# MD related packages\n",
    "import nglview\n",
    "\n",
    "# OpenFF package.\n",
    "# Do not use from openff.xxx import ... to avoid name collisions.\n",
    "import openff.toolkit.topology\n",
    "import openff.toolkit.typing.engines.smirnoff\n",
    "import pandas\n",
    "import requests\n",
    "\n",
    "# Custom functions defined in the current directory\n",
    "from ligands import *\n",
    "from openmm import *\n",
    "from openmm.app import *\n",
    "from openmm.unit import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Download a molecule from PubChem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cid = \"3672\"  # This is the pubchem ID for ibuprofen\n",
    "fn_sdf = f\"CID_{cid}.sdf\"\n",
    "if not os.path.isfile(fn_sdf):\n",
    "    url = f\"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/{cid}/record/SDF/?record_type=3d&response_type=save\"\n",
    "    with open(f\"CID_{cid}.sdf\", \"w\") as f:\n",
    "        f.write(requests.get(url).text)\n",
    "fn_pdb = f\"CID_{cid}.pdb\"\n",
    "convert_sdf_to_pdb(fn_sdf, fn_pdb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Assigning Sage 2.0.0 parameters with the SMIRNOFF engine\n",
    "\n",
    "Only in the following code, the Open Force Field toolkit is used.\n",
    "\n",
    "The first step is to maken an OpenFF Molecule object, here by loading the SDF file.\n",
    "OpenFF Molecule objects can be made in [many different ways](https://open-forcefield-toolkit.readthedocs.io/en/stable/users/molecule_cookbook.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the topology from the SDF file.\n",
    "molecule = openff.toolkit.topology.Molecule(fn_sdf)\n",
    "# Quick and dirty visualization of the molecule.\n",
    "molecule.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will derive an OpenMM system object, using the Sage force field.\n",
    "This looks simple, yet a lot is going on under the hood to derive force field parameters for this specific molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "force_field = openff.toolkit.typing.engines.smirnoff.ForceField(\"openff-2.0.0.offxml\")\n",
    "system = force_field.create_openmm_system(molecule.to_topology())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Perform a short molecular dynamics simulation\n",
    "\n",
    "The code below is essentially a copy from the molecular dynamics simulations on alanine dipeptide in section 2.\n",
    "\n",
    "You should receive warnings about duplicates atoms, which can be ignored. These are caused by a few too simple choices in the PDB file generated by OpenBabel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup the MD\n",
    "pdb = PDBFile(fn_pdb)\n",
    "\n",
    "integrator = LangevinIntegrator(300 * kelvin, 1 / picosecond, 2 * femtoseconds)\n",
    "simulation = Simulation(pdb.topology, system, integrator)\n",
    "simulation.context.setPositions(pdb.positions)\n",
    "simulation.minimizeEnergy()\n",
    "\n",
    "# Write the initial state back to a PDB, could be useful\n",
    "# for debugging.\n",
    "with open(\"init_01.pdb\", \"w\") as f:\n",
    "    PDBFile.writeFile(simulation.topology, pdb.positions, f)\n",
    "\n",
    "# Set the reporters collecting the MD output.\n",
    "simulation.reporters = []\n",
    "simulation.reporters.append(DCDReporter(\"traj_01.dcd\", 100))\n",
    "simulation.reporters.append(\n",
    "    StateDataReporter(stdout, 1000, step=True, temperature=True, elapsedTime=True)\n",
    ")\n",
    "simulation.reporters.append(\n",
    "    StateDataReporter(\n",
    "        \"scalars_01.csv\",\n",
    "        100,\n",
    "        time=True,\n",
    "        potentialEnergy=True,\n",
    "        totalEnergy=True,\n",
    "        temperature=True,\n",
    "    )\n",
    ")\n",
    "simulation.step(10000)\n",
    "\n",
    "# The last line is only needed for Windows users,\n",
    "# to close the DCD file before it can be opened by nglview.\n",
    "del simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the trajectory.\n",
    "view = nglview.show_mdtraj(mdtraj.load(\"traj_01.dcd\", top=\"init_01.pdb\"))\n",
    "view.clear_representations()\n",
    "view.add_licorice()\n",
    "view.add_unitcell()\n",
    "view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot temperature as an initial verification of convergence.\n",
    "df = pandas.read_csv(\"scalars_01.csv\")\n",
    "df.plot(kind=\"line\", x='#\"Time (ps)\"', y=\"Temperature (K)\")"
   ]
  }
 ],
 "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}