{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# System Building: μ-opioid Receptor in Membrane with morphinan antagonist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial, we will showcase how to build a protein system embeded in a membrane with a ligand in one of the compartments for simulating binding. The sample system is a mu-opioid receptor (the protein) in a membrane and a morphinan antagonist (the ligand).\n",
    "\n",
    "Let's start by doing some imports and definitions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. \n",
      "https://dx.doi.org/10.1021/acs.jctc.6b00049\n",
      "Documentation: http://software.acellera.com/\n",
      "To update: conda update htmd -c acellera -c psi4\n",
      "\n",
      "You are on the latest HTMD version (unpackaged : /home/joao/maindisk/software/repos/Acellera/htmd/htmd).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.home import home\n",
    "from os.path import join\n",
    "config(viewer='webgl')\n",
    "datadir = home(dataDir='mor')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Prepare the protein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "View the file as it comes from the OPM database:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b00d446af1b743acacbc6209a6501688",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Molecule(join(datadir, '4dkl.pdb')).view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Retrieve the structure from OPM, do not keep the `DUM` atoms, and print the thickness as calculated by OPM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:54:49,428 - htmd.molecule.molecule - INFO - Removed 2546 atoms. 4836 atoms remaining in the molecule.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "32.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from htmd.util import opm\n",
    "prot, thickness = opm('4dkl', keep=False)\n",
    "thickness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove non-protein atoms, and keep only a monomer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:54:51,462 - htmd.molecule.molecule - INFO - Removed 2574 atoms. 2262 atoms remaining in the molecule.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([   0,    1,    2, ..., 4808, 4809, 4810], dtype=int32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Automatically detecting segments and assigning names to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:54:53,807 - htmd.builder.builder - INFO - Created segment P0 between resid 65 and 263.\n",
      "2018-03-17 01:54:53,808 - htmd.builder.builder - INFO - Created segment P1 between resid 270 and 352.\n"
     ]
    }
   ],
   "source": [
    "prot = autoSegment(prot,'protein')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build the protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:54:57,172 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-03-17 01:54:57,371 - htmd.builder.builder - INFO - One disulfide bond was added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 3005 resid 140 resname CYS chain B segid P0]\n",
      "             B: [serial 3615 resid 217 resname CYS chain B segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:54:57,912 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-03-17 01:54:58,939 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 1 atoms due to bad angles.\n",
      "2018-03-17 01:54:58,941 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 78 atoms.\n",
      "2018-03-17 01:54:58,941 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/tmp-build/log.txt for further information.\n",
      "2018-03-17 01:54:58,943 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "topos  = charmm.defaultTopo() + [join(datadir, 'ff.rtf')]\n",
    "params = charmm.defaultParam() + [join(datadir, 'ff.prm')]\n",
    "\n",
    "prot = charmm.build(prot, topo=topos, param=params, outdir='./tmp-build', ionize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0fbf09d78c23499cb0946b70fc74c412",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "prot.reps.add(sel='segid P0', style='NewCartoon', color=1)\n",
    "prot.reps.add(sel='segid P1', style='NewCartoon', color=2)\n",
    "prot.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Add a sodium atom in the receptor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "28384883aed84f4b99476a2ed5f343ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sod = Molecule(join(datadir, 'sod.pdb'))\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "prot.reps.add(sel='ions', style='VDW', color='green')\n",
    "prot.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Embed the protein into a membrane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "memb = Molecule(join(datadir, 'membrane80by80C36.pdb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center the membrane onto the protein center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcenter = np.mean(prot.get('coords','protein'),axis=0)\n",
    "mcenter = np.mean(memb.get('coords'),axis=0)\n",
    "memb.moveBy(pcenter-mcenter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And now embed. \n",
    "\n",
    "The two are equivalent - `append` with `collisions=True` only\n",
    "adds atoms if they do not clash  sterically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:55:22,309 - htmd.molecule.molecule - INFO - Removed 339 residues from appended Molecule due to collisions.\n"
     ]
    }
   ],
   "source": [
    "mol = prot.copy()\n",
    "mol.append(memb, collisions=True)\n",
    "# mol = embed(prot,memb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize the embedded system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9ff8c68ead914dcba0649b4f505e75e9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')\n",
    "mol.reps.add(sel='ions', style='VDW', color='green')\n",
    "mol.reps.add(sel='lipids', style='Lines')\n",
    "mol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Add a ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "lig = Molecule(join(datadir, 'QM-min.pdb'))\n",
    "lig.set('segid','L');\n",
    "lcenter = np.mean(lig.get('coords'),axis=0)\n",
    "newlcenter=[random.uniform(-10, 10), random.uniform(-10, 10),  43 ]\n",
    "lig.rotateBy(uniformRandomRotation(), lcenter)\n",
    "lig.moveBy(newlcenter-lcenter)\n",
    "mol.append(lig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Put it in a water box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:55:41,667 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/Acellera/htmd/htmd/builder/wat.pdb\n",
      "2018-03-17 01:55:43,163 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n",
      "Solvating: 100%|██████████| 8/8 [00:17<00:00,  2.20s/it]\n",
      "2018-03-17 01:56:03,643 - htmd.builder.solvate - INFO - 9117 water molecules were added to the system.\n"
     ]
    }
   ],
   "source": [
    "coo = mol.get('coords','noh and (lipids or protein)')\n",
    "m = np.min(coo, axis=0) + [0, 0, -5]\n",
    "M = np.max(coo, axis=0) + [0, 0, 20]\n",
    "smol = solvate(mol, minmax=np.vstack((m,M)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "03423611d0e248bf8b536fa43dc97d2b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "smol.reps.add(sel='segid L', style='Licorice')\n",
    "smol.reps.add(sel='water', style='Lines')\n",
    "smol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build with CHARMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:56:35,041 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-03-17 01:56:53,448 - htmd.builder.builder - INFO - One disulfide bond was added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 1212 resid 140 resname CYS chain B segid P0]\n",
      "             B: [serial 2448 resid 217 resname CYS chain B segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 01:56:53,971 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-03-17 01:56:54,837 - htmd.builder.charmm - INFO - Finished building.\n",
      "2018-03-17 01:56:58,567 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 66 ions for the given salt concentration.\n",
      "2018-03-17 01:56:59,001 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2018-03-17 01:56:59,002 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2018-03-17 01:56:59,003 - htmd.builder.ionize - INFO - Placing 80 ions.\n",
      "2018-03-17 01:57:36,506 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-03-17 01:57:54,859 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-03-17 01:57:55,740 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "molbuilt = charmm.build(smol, topo=topos, param=params, outdir='./final-build', saltconc=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize built system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "37951b6d173f4a35b28b225ad80c41b3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "molbuilt.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The `molbuilt` is a `Molecule` object that contains the built system, but the full contents to run a simulation are located in the `outdir` (`./final-build` in this case)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.4"
  },
  "widgets": {
   "state": {
    "1741c23784064c52b59c79251cfaf8e6": {
     "views": [
      {
       "cell_index": 14
      }
     ]
    },
    "223dfc92ccae4887a46893321ac12ed8": {
     "views": [
      {
       "cell_index": 14
      }
     ]
    },
    "416b4c797ae44b82bb5e4406a6e73a49": {
     "views": [
      {
       "cell_index": 32
      }
     ]
    },
    "594aa85570984b3887612502c4597e7c": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    },
    "79e1a424c2554d5fb4bee580cb01833b": {
     "views": [
      {
       "cell_index": 3
      }
     ]
    },
    "9121f11441b2485fab877a116bde2606": {
     "views": [
      {
       "cell_index": 28
      }
     ]
    },
    "927a4600d16a4577b29a137506f3cc17": {
     "views": [
      {
       "cell_index": 3
      }
     ]
    },
    "c1825c2b1a3c4196a7d995f244f98aef": {
     "views": [
      {
       "cell_index": 32
      }
     ]
    },
    "e629564e29704ce6826a175f5c451ac2": {
     "views": [
      {
       "cell_index": 12
      }
     ]
    },
    "f832672917214baa99112a489918050e": {
     "views": [
      {
       "cell_index": 12
      }
     ]
    },
    "fb0fd032ceaf4485aaa81c7d1ad8ee3b": {
     "views": [
      {
       "cell_index": 22
      }
     ]
    },
    "fc1dea2f053a4c3aa924e70e29b91df6": {
     "views": [
      {
       "cell_index": 28
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}