01-convert-benchmarkset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert a host-guest pair from `https://github.com/MobleyLab/benchmarksets/`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**CB7-memantine**\n",
    "\n",
    "The initial bonded and Lennard-Jones parameters are GAFF v1.7 with partial charges generated using RESP, and conformation given by docking the host and guest with MOE. For more information, see [here](https://github.com/MobleyLab/benchmarksets/tree/d9bd05719fe42a390442d3984eccec591ec32950/input_files)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conversion strategy is outlined in the [README](README.md). Briefly, we will download the files from the `benchmarksets` repository, create a PDB, extract the host and guest as topologies, re-parameterize those molecules with SMIRNOFF99Frosst bond, angle, torsion, and Lennard-Jones parameters, then parameterize the water and ions with TIP3P and Joung-Cheatham, merge the structures, and write out a new set of combined parameters and coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:16.323922Z",
     "start_time": "2018-03-06T22:37:16.239331Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import os as os\n",
    "import urllib.request\n",
    "\n",
    "import parmed as pmd\n",
    "\n",
    "from openforcefield.typing.engines.smirnoff import ForceField, unit\n",
    "from openforcefield.utils import mergeStructure\n",
    "\n",
    "from smirnovert.utils import (create_pdb_with_conect, prune_conect, split_topology, create_host_guest_topology,\n",
    "                    create_host_mol2, convert_mol2_to_sybyl_antechamber,\n",
    "                    load_mol2, check_unique_atom_names,\n",
    "                    check_bond_lengths,\n",
    "                    extract_water_and_ions, create_water_and_ions_parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we begin, let's specify the starting file names and prefix for the intermediary files. We are also going to write a bunch of temporary files that can be cleaned up later, but for debugging, I leave them. Also, the `utils.py` functions using the `logging` module, so we can specify how much information we want. Here I'll set the logging level to `INFO`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:16.412621Z",
     "start_time": "2018-03-06T22:37:16.325376Z"
    }
   },
   "outputs": [],
   "source": [
    "import logging\n",
    "from importlib import reload\n",
    "# `logging` needs to be reloaded, because `jupyter notebook` itself \n",
    "# uses the logging module to print messages to standard output...\n",
    "reload(logging)\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)\n",
    "logging.basicConfig(\n",
    "    format='%(asctime)s %(message)s', datefmt='%Y-%m-%d %I:%M:%S %p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:16.507618Z",
     "start_time": "2018-03-06T22:37:16.414078Z"
    }
   },
   "outputs": [],
   "source": [
    "test_case = 'cb7-1/'\n",
    "destination = './tests/' + test_case\n",
    "reference_prmtop = 'cb7-1.prmtop'\n",
    "reference_inpcrd = 'cb7-1.rst7'\n",
    "\n",
    "\n",
    "prefix = 'cb7-1'\n",
    "host_resname = 'CB7'\n",
    "guest_resname = 'MOL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's download a fresh host-guest example from David Mobley's `benchmarksets` GitHub repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:17.176960Z",
     "start_time": "2018-03-06T22:37:16.509116Z"
    }
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(destination):\n",
    "    os.makedirs(destination)\n",
    "    \n",
    "request = urllib.request.urlretrieve(\"https://rawgit.com/MobleyLab/benchmarksets/master/input_files/cb7-set1/prmtop-rst7/cb7-1.prmtop\", destination + reference_prmtop)\n",
    "request = urllib.request.urlretrieve(\"https://rawgit.com/MobleyLab/benchmarksets/master/input_files/cb7-set1/prmtop-rst7/cb7-1.rst7\", destination + reference_inpcrd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, it is useful to visualize the structure with e.g., `nglview`.\n",
    "\n",
    "Before we start the conversion, let's grab the box vectors from the files -- this will come in handy later, because we lose this information when we merge the structures down below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:17.487296Z",
     "start_time": "2018-03-06T22:37:17.178518Z"
    }
   },
   "outputs": [],
   "source": [
    "reference = pmd.load_file(destination + reference_prmtop, xyz=destination + reference_inpcrd)\n",
    "box = reference.box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll create a proper PDB from these AMBER files, and delete the `CONECT` records that are solvent-solvent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:17.691898Z",
     "start_time": "2018-03-06T22:37:17.489040Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:17 PM Creating ./tests/cb7-1/cb7-1.pdb with CONECT records...\n"
     ]
    }
   ],
   "source": [
    "create_pdb_with_conect(solvated_pdb=destination + reference_inpcrd,\n",
    "                      amber_prmtop=destination + reference_prmtop,\n",
    "                      output_pdb=destination + prefix + '.pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to prune the `CONECT` records to deal with https://github.com/openforcefield/openforcefield/issues/68."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:17.778923Z",
     "start_time": "2018-03-06T22:37:17.693626Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:17 PM Pruning water-water CONECT records...\n"
     ]
    }
   ],
   "source": [
    "prune_conect(input_pdb=prefix + '.pdb',\n",
    "            output_pdb=prefix + '.pruned.pdb',\n",
    "            path=destination)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll split the PDB into separate topology objects, and extract the host-guest topology for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:18.971509Z",
     "start_time": "2018-03-06T22:37:17.780560Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:17 PM Splitting topology into components...\n",
      "2018-03-06 02:37:18 PM Creating a combined topology for the host and guest molecules...\n"
     ]
    }
   ],
   "source": [
    "components = split_topology(file_name=destination + prefix + '.pruned.pdb')\n",
    "hg_topology = create_host_guest_topology(components, \n",
    "                                         host_resname=host_resname, \n",
    "                                         guest_resname=guest_resname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also need to create `mol2` files **with SYBYL atom types** for the host and guest, that we can use later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:19.211455Z",
     "start_time": "2018-03-06T22:37:18.973783Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:19 PM Writing a `mol2` for the host molecule...\n",
      "2018-03-06 02:37:19 PM Writing a `mol2` for the host molecule...\n",
      "2018-03-06 02:37:19 PM Converting ./tests/cb7-1/CB7.mol2 to SYBYL atom types via Antechamber...\n",
      "2018-03-06 02:37:19 PM Converting ./tests/cb7-1/MOL.mol2 to SYBYL atom types via Antechamber...\n"
     ]
    }
   ],
   "source": [
    "create_host_mol2(\n",
    "    solvated_pdb=destination + prefix + '.pruned.pdb',\n",
    "    amber_prmtop=destination + reference_prmtop,\n",
    "    mask=host_resname,\n",
    "    output_mol2=destination + host_resname + '.mol2')\n",
    "\n",
    "create_host_mol2(\n",
    "    solvated_pdb=destination + prefix + '.pdb',\n",
    "    amber_prmtop=destination + reference_prmtop,\n",
    "    mask=guest_resname,\n",
    "    output_mol2=destination + guest_resname + '.mol2')\n",
    "\n",
    "convert_mol2_to_sybyl_antechamber(\n",
    "    input_mol2=destination + host_resname + '.mol2',\n",
    "    output_mol2=destination + host_resname + '-sybyl.mol2',\n",
    "    ac_doctor=False)\n",
    "\n",
    "convert_mol2_to_sybyl_antechamber(\n",
    "    input_mol2=destination + guest_resname + '.mol2',\n",
    "    output_mol2=destination + guest_resname + '-sybyl.mol2',\n",
    "    ac_doctor=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to separate out waters and ions from the starting AMBER files.\n",
    "\n",
    "It is a little confusing that the following function runs with `dummy_atoms=True`. That setting means don't specifically strip dummy atoms from the reference coordinates, and is correct when dummy atoms are not present. Yet the second function call needs `dummy_atoms=False` because dummy atoms do *not* need parameters, if they are absent. This could be made more clear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:20.606655Z",
     "start_time": "2018-03-06T22:37:19.213407Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:19 PM Extracting water and ions from cb7-1.prmtop...\n",
      "2018-03-06 02:37:19 PM Creating parameters for the waters and ions...\n"
     ]
    }
   ],
   "source": [
    "extract_water_and_ions(\n",
    "    amber_prmtop=reference_prmtop,\n",
    "    amber_inpcrd=reference_inpcrd,\n",
    "    host_residue=':' + host_resname,\n",
    "    guest_residue=':' + guest_resname,\n",
    "    dummy_atoms=True,\n",
    "    output_pdb='water_ions.pdb',\n",
    "    path=destination)\n",
    "\n",
    "create_water_and_ions_parameters(\n",
    "    input_pdb='water_ions.pdb',\n",
    "    output_prmtop='water_ions.prmtop',\n",
    "    output_inpcrd='water_ions.inpcrd',\n",
    "    dummy_atoms=False,\n",
    "    path=destination)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll create `mol2` files with **SYBYL atom types** for the water and ions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:20.640396Z",
     "start_time": "2018-03-06T22:37:20.608512Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:20 PM Loading ./tests/cb7-1/CB7-sybyl.mol2...\n",
      "2018-03-06 02:37:20 PM Loading ./tests/cb7-1/MOL-sybyl.mol2...\n",
      "2018-03-06 02:37:20 PM Checking all atoms have unique names...\n",
      "2018-03-06 02:37:20 PM Checking all atoms have unique names...\n"
     ]
    }
   ],
   "source": [
    "host = load_mol2(\n",
    "    filename=destination + host_resname + '-sybyl.mol2',\n",
    "    name=host_resname,\n",
    "    add_tripos=True)\n",
    "\n",
    "guest = load_mol2(\n",
    "    filename=destination + guest_resname + '-sybyl.mol2',\n",
    "    name=guest_resname,\n",
    "    add_tripos=False)\n",
    "\n",
    "check_unique_atom_names(host)\n",
    "check_unique_atom_names(guest)\n",
    "molecules = [host, guest]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can create the OpenMM system for the host and guest with SMIRNOFF99Frosst parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:21.485767Z",
     "start_time": "2018-03-06T22:37:20.641865Z"
    }
   },
   "outputs": [],
   "source": [
    "ff = ForceField('forcefield/smirnoff99Frosst.ffxml')\n",
    "system = ff.createSystem(\n",
    "    hg_topology.topology,\n",
    "    molecules,\n",
    "    nonbondedCutoff=1.1 * unit.nanometer,\n",
    "    ewaldErrorTolerance=1e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll convert the OpenMM system to a ParmEd structure and check for bad bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:21.597315Z",
     "start_time": "2018-03-06T22:37:21.487221Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-06 02:37:21 PM Checking structure for bonds >4 A...\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check if the host-guest parameter file already exists...\n",
      "Check if the host-guest coordinate file already exists...\n"
     ]
    }
   ],
   "source": [
    "hg_structure = pmd.openmm.topsystem.load_topology(\n",
    "    hg_topology.topology, system, hg_topology.positions)\n",
    "\n",
    "check_bond_lengths(hg_structure, threshold=4)\n",
    "\n",
    "try:\n",
    "    hg_structure.save(destination + 'hg.prmtop')\n",
    "except OSError:\n",
    "    print(\n",
    "        'Check if the host-guest parameter file already exists...')\n",
    "\n",
    "try:\n",
    "    hg_structure.save(destination + 'hg.inpcrd')\n",
    "except OSError:\n",
    "    print(\n",
    "        'Check if the host-guest coordinate file already exists...')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll load the water and ions into a ParmEd structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:21.910313Z",
     "start_time": "2018-03-06T22:37:21.599429Z"
    }
   },
   "outputs": [],
   "source": [
    "water_and_ions = pmd.amber.AmberParm(\n",
    "    destination + 'water_ions.prmtop',\n",
    "    xyz=destination + 'water_ions.inpcrd')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's merge the host-guest structure with the dummy atoms, the waters, and the ions, and set the box coordinates of the merged structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:22.322554Z",
     "start_time": "2018-03-06T22:37:21.911933Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Check if solvated parameter file already exists...\n",
      "Check if solvated coordinate file already exists...\n"
     ]
    }
   ],
   "source": [
    "merged = mergeStructure(hg_structure, water_and_ions)\n",
    "merged.box = reference.box\n",
    "try:\n",
    "    merged.save(destination + 'smirnoff.prmtop')\n",
    "except:\n",
    "    print('Check if solvated parameter file already exists...')\n",
    "try:\n",
    "    merged.save(destination + 'smirnoff.inpcrd')\n",
    "except:\n",
    "    print('Check if solvated coordinate file already exists...')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the parmaeters in a ParmEd structure, we can do quick comparisons between the \"reference\" (i.e., GAFF v1.7 here) parameters and SMIRNOFF99Frosst parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-06T22:37:22.362494Z",
     "start_time": "2018-03-06T22:37:22.324061Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      2 N2   (2   )      30 C2   (30  )     1.3350   490.0000\n",
      "      2 N2   (n   )      30 C2   (c   )     1.3789   427.6000\n"
     ]
    }
   ],
   "source": [
    "for smirnoff_bond, reference_bond in zip(merged.bonds, reference.bonds):\n",
    "    smirnoff_atom1, smirnoff_atom2 = smirnoff_bond.atom1, smirnoff_bond.atom2\n",
    "    reference_atom1, reference_atom2 = reference_bond.atom1, reference_bond.atom2\n",
    "    \n",
    "    if smirnoff_bond.type is not None :\n",
    "        if smirnoff_atom1.name == 'N2' and smirnoff_atom2.name == 'C2':\n",
    "            atom1 = smirnoff_atom1\n",
    "            atom2 = smirnoff_atom2\n",
    "            bond = smirnoff_bond\n",
    "            print(\n",
    "                f'{atom1.idx + 1:7d} {atom1.name:4} ({atom1.type:4}) {atom2.idx + 1:7d} '\n",
    "                f'{atom2.name:4} ({atom2.type:4}) {bond.type.req:10.4f} {bond.type.k:10.4f}'\n",
    "            )\n",
    "        if reference_atom1.name == 'N2' and reference_atom2.name == 'C2':\n",
    "            atom1 = reference_atom1\n",
    "            atom2 = reference_atom2\n",
    "            bond = reference_bond\n",
    "            print(\n",
    "                f'{atom1.idx + 1:7d} {atom1.name:4} ({atom1.type:4}) {atom2.idx + 1:7d} '\n",
    "                f'{atom2.name:4} ({atom2.type:4}) {bond.type.req:10.4f} {bond.type.k:10.4f}'\n",
    "            )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nice, these parameters look identical, but note the difference in atom type!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:smirnoff-host-guest]",
   "language": "python",
   "name": "conda-env-smirnoff-host-guest-py"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}