{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from htmd.ui import *\n",
    "config(viewer='ngl')\n",
    "os.chdir('/webdata/nc983hu3brda/') # Don't use this command"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using docking to generate starting poses for simulations\n",
    "\n",
    "Download the files for this tutorial from this [link](http://pub.htmd.org/nc983hu3brda/bentryp.tar.gz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dock the protein with the ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prot = Molecule('bentryp/trypsin.pdb')\n",
    "prot.center()\n",
    "lig = Molecule('bentryp/benzamidine.pdb')\n",
    "poses, scores = dock(prot, lig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the docked poses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mol = Molecule()\n",
    "mol.append(prot)\n",
    "for i, p in enumerate(poses):\n",
    "    mol.append(p)\n",
    "mol.view(sel='protein', style='NewCartoon', hold=True)\n",
    "mol.view(sel='resname MOL', style='Licorice', color=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build systems from docked poses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "molbuilt = []\n",
    "for i, p in enumerate(poses):\n",
    "    prot = Molecule('bentryp/trypsin.pdb')\n",
    "    prot.filter('chain A and (protein or water or resname CA)')\n",
    "    prot.set('segid', 'P', sel='protein and noh')\n",
    "    prot.set('segid', 'W', sel='water')\n",
    "    prot.set('segid', 'CA', sel='resname CA')\n",
    "    prot.center()\n",
    "    from moleculekit.util import maxDistance\n",
    "    D = maxDistance(prot, 'all')\n",
    "    \n",
    "    ligand = p\n",
    "    ligand.set('segid','L')\n",
    "    ligand.set('resname','MOL')\n",
    "    \n",
    "    mol = Molecule(name='combo')\n",
    "    mol.append(prot)\n",
    "    mol.append(ligand)\n",
    "    \n",
    "    D = D + 15\n",
    "    smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])\n",
    "    topos  = ['top/top_all22star_prot.rtf', 'top/top_water_ions.rtf', './bentryp/benzamidine.rtf']\n",
    "    params = ['par/par_all22star_prot.prm', 'par/par_water_ions.prm', './bentryp/benzamidine.prm']\n",
    "\n",
    "    molbuilt.append(charmm.build(smol, topo=topos, param=params, outdir='./docked/build/{}/'.format(i+1), saltconc=0.15))\n",
    "    if i==1: # For time purposes lets only build the two first\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from ipywidgets.widgets import Box\n",
    "w = []\n",
    "for i, m in enumerate(molbuilt):\n",
    "    m.view(sel='protein', style='NewCartoon', hold=True)\n",
    "    m.view(sel='water', style='Lines', hold=True)\n",
    "    h = m.view(sel='resname MOL', style='Licorice', color=0)\n",
    "    w.append(h)\n",
    "b = Box(children=(w[0],w[1]))\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate the build systems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.equilibration_v1 import Equilibration\n",
    "md = Equilibration()\n",
    "md.numsteps = 1000\n",
    "md.temperature = 298\n",
    "\n",
    "builds = natsort(glob('docked/build/*/'))\n",
    "for i, b in enumerate(builds):\n",
    "    md.write(b, 'docked/equil/{}/'.format(i+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mdx = AcemdLocal()\n",
    "mdx.submit(glob('./docked/equil/*/'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mdx.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Create the production folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.production_v1 import Production\n",
    "md = Production()\n",
    "md.acemd.run = '50ns'\n",
    "md.temperature = 298\n",
    "\n",
    "equils = natsort(glob('docked/equil/*/'))\n",
    "for i, b in enumerate(equils):\n",
    "    md.write(b, 'docked/generators/{}/'.format(i+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mdx = AcemdLocal()\n",
    "mdx.submit(glob('./docked/generators/*/'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mdx.wait()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}