{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Please cite -- HTMD: High-Throughput Molecular Dynamics for Molecular Discovery\n",
      "J. Chem. Theory Comput., 2016, 12 (4), pp 1845-1852. \n",
      "http://pubs.acs.org/doi/abs/10.1021/acs.jctc.6b00049\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",
    "config(viewer='ngl')\n",
    "#os.chdir('/webdata/73hboiwia98hdj209jq0/')  # Skip this command."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# System Building μ-opioid Receptor in Membrane"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Prepare the protein\n",
    "\n",
    "(Download the required files for the tutorial from this [link](http://pub.htmd.org/73hboiwia98hdj209jq0/mor.tar.gz).)\n",
    "\n",
    "View the file as it comes from the OPM database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Molecule('mor/4dkl.pdb').view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Retrieve and clean up `DUM` marker atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:20:06,384 - 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')\n",
    "thickness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove non-protein atoms, and keep only a monomer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:20:10,482 - htmd.molecule.molecule - INFO - Removed 2574 atoms. 2262 atoms remaining in the molecule.\n"
     ]
    }
   ],
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:20:14,499 - htmd.builder.builder - INFO - Created segment P0 between resid 65 and 263.\n",
      "2017-01-10 23:20:14,501 - 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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:20:19,379 - htmd.builder.charmm - INFO - Writing out segments.\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": [
      "2017-01-10 23:20:20,584 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2017-01-10 23:20:20,753 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2017-01-10 23:20:20,865 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 1 atoms due to bad angles.\n",
      "2017-01-10 23:20:20,867 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 78 atoms.\n",
      "2017-01-10 23:20:20,868 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/SANDBOX/03_PreparationBuilding_TGiorgino/03c_out_morbuild/prot/log.txt for further information.\n",
      "2017-01-10 23:20:20,869 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "topos  = charmm.defaultTopo() + ['mor/ff.rtf']\n",
    "params = charmm.defaultParam() + ['mor/ff.prm']\n",
    "\n",
    "prot = charmm.build(prot, topo=topos, param=params, \n",
    "                    outdir='./03c_out_morbuild/prot/', ionize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sod = Molecule('mor/sod.pdb')\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "prot.reps.add(sel='ions', style='VDW')\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": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "memb = Molecule('mor/membrane80by80C36.pdb')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center the membrane onto the protein center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:20:51,502 - htmd.molecule.molecule - INFO - Removed 305 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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')\n",
    "mol.reps.add(sel='ions', style='VDW')\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": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lig = Molecule('mor/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": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:35:27,474 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/Acellera/htmd/htmd/builder/wat.pdb\n",
      "2017-01-10 23:35:28,940 - htmd.builder.solvate - INFO - Replicating 8 water segments, 2 by 2 by 2\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solvating: 100% (8/8) [############################################] eta 00:00 /\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:36:06,230 - htmd.builder.solvate - INFO - After removing water molecules colliding with other molecules, 9116 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": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "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": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:36:50,359 - htmd.builder.charmm - INFO - Writing out segments.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 1212 resid 140 resname CYS chain  segid P0]\n",
      "             B: [serial 2448 resid 217 resname CYS chain  segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-01-10 23:38:48,651 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2017-01-10 23:38:48,830 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2017-01-10 23:38:49,628 - htmd.builder.charmm - INFO - Finished building.\n",
      "2017-01-10 23:38:53,686 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 66 ions for the given salt concentration.\n",
      "2017-01-10 23:38:54,202 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2017-01-10 23:38:54,203 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2017-01-10 23:38:54,204 - htmd.builder.ionize - INFO - Placing 80 ions.\n",
      "2017-01-10 23:39:48,897 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2017-01-10 23:41:40,310 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2017-01-10 23:41:41,095 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "molbuilt = charmm.build(smol, topo=topos, param=params, \n",
    "                        outdir='./build/',\n",
    "                        saltconc=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Visualize built system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "molbuilt.view()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "widgets": {
   "state": {
    "1741c23784064c52b59c79251cfaf8e6": {
     "views": [
      {
       "cell_index": 14.0
      }
     ]
    },
    "223dfc92ccae4887a46893321ac12ed8": {
     "views": [
      {
       "cell_index": 14.0
      }
     ]
    },
    "416b4c797ae44b82bb5e4406a6e73a49": {
     "views": [
      {
       "cell_index": 32.0
      }
     ]
    },
    "594aa85570984b3887612502c4597e7c": {
     "views": [
      {
       "cell_index": 22.0
      }
     ]
    },
    "79e1a424c2554d5fb4bee580cb01833b": {
     "views": [
      {
       "cell_index": 3.0
      }
     ]
    },
    "9121f11441b2485fab877a116bde2606": {
     "views": [
      {
       "cell_index": 28.0
      }
     ]
    },
    "927a4600d16a4577b29a137506f3cc17": {
     "views": [
      {
       "cell_index": 3.0
      }
     ]
    },
    "c1825c2b1a3c4196a7d995f244f98aef": {
     "views": [
      {
       "cell_index": 32.0
      }
     ]
    },
    "e629564e29704ce6826a175f5c451ac2": {
     "views": [
      {
       "cell_index": 12.0
      }
     ]
    },
    "f832672917214baa99112a489918050e": {
     "views": [
      {
       "cell_index": 12.0
      }
     ]
    },
    "fb0fd032ceaf4485aaa81c7d1ad8ee3b": {
     "views": [
      {
       "cell_index": 22.0
      }
     ]
    },
    "fc1dea2f053a4c3aa924e70e29b91df6": {
     "views": [
      {
       "cell_index": 28.0
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}