{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparation of the $\\mu$ opioid receptor with ligand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a complex build system as it has several components, the protein, a sodium ion, the ligand and of course the membrane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Please cite HTMD: Doerr et al.(2016)JCTC,12,1845. https://dx.doi.org/10.1021/acs.jctc.6b00049\n",
      "\n",
      "HTMD Documentation at: https://www.htmd.org/docs/latest/\n",
      "\n",
      "You are on the latest HTMD version (unpackaged : /shared/sdoerr/Work/htmdacellera/htmd).\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from htmd.ui import *\n",
    "from htmd.home import home\n",
    "#get the files\n",
    "shutil.copytree(home()+'/data/mor','/tmp/testmor/pdb')\n",
    "os.chdir('/tmp/testmor')\n",
    "path='./01_prepare/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4dkl.pdb  ff.prm  ff.rtf  membrane80by80C36.pdb  QM-min.pdb  sod.pdb\r\n"
     ]
    }
   ],
   "source": [
    "%ls /tmp/testmor/pdb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-02-20 10:33:51,140 - htmd.molecule.molecule - INFO - Removed 5120 atoms. 2262 atoms remaining in the molecule.\n",
      "2017-02-20 10:33:51,274 - htmd.builder.builder - INFO - Created segment P0 between resid 65 and 263.\n",
      "2017-02-20 10:33:51,275 - htmd.builder.builder - INFO - Created segment P1 between resid 270 and 352.\n",
      "2017-02-20 10:33:51,673 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2017-02-20 10:33:52,150 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2017-02-20 10:33:52,258 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2017-02-20 10:33:52,326 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 1 atoms due to bad angles.\n",
      "2017-02-20 10:33:52,327 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 78 atoms.\n",
      "2017-02-20 10:33:52,327 - htmd.builder.charmm - WARNING - Please check /tmp/testmor/01_prepare/prot/log.txt for further information.\n",
      "2017-02-20 10:33:52,328 - htmd.builder.charmm - INFO - Finished building.\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"
     ]
    }
   ],
   "source": [
    "#Protein 4dkl is taken from opm\n",
    "\n",
    "topos  = charmm.defaultTopo() + ['pdb/ff.rtf']\n",
    "params = charmm.defaultParam() + ['pdb/ff.prm']\n",
    "prot = Molecule('pdb/4dkl.pdb')\n",
    "prot.filter('protein and noh and chain B or water within 5 of (chain B and protein)')\n",
    "pcenter = np.mean(prot.get('coords','protein'), axis=0)\n",
    "prot = autoSegment(prot, sel='protein') \n",
    "\n",
    "prot = charmm.build(prot, topo=topos, param=params, outdir= path+'prot',ionize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prot.view()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2017-02-20 10:41:20,248 - htmd.molecule.molecule - INFO - Removed 311 residues from appended Molecule due to collisions.\n",
      "2017-02-20 10:41:20,896 - htmd.builder.solvate - INFO - Using water pdb file at: /shared/sdoerr/Work/htmdacellera/htmd/builder/wat.pdb\n",
      "2017-02-20 10:41:21,793 - 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-02-20 10:41:43,603 - htmd.builder.solvate - INFO - After removing water molecules colliding with other molecules, 10037 water molecules were added to the system.\n",
      "2017-02-20 10:41:52,504 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2017-02-20 10:43:02,342 - htmd.builder.builder - INFO - One disulfide bond was added\n",
      "2017-02-20 10:43:02,454 - htmd.builder.charmm - INFO - Starting the build.\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-02-20 10:43:02,930 - htmd.builder.charmm - INFO - Finished building.\n",
      "2017-02-20 10:43:05,586 - htmd.builder.ionize - INFO - Adding 14 anions + 0 cations for neutralizing and 70 ions for the given salt concentration.\n",
      "2017-02-20 10:43:05,939 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2017-02-20 10:43:05,940 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2017-02-20 10:43:05,940 - htmd.builder.ionize - INFO - Placing 84 ions.\n",
      "2017-02-20 10:43:43,330 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2017-02-20 10:44:55,163 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2017-02-20 10:44:55,631 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "#Add sodium in the receptor\n",
    "sod = Molecule('pdb/sod.pdb')\n",
    "sod.set('segid','S1')\n",
    "prot.append(sod)\n",
    "\n",
    "#Use a POPC membrane created with vmd and C36\n",
    "memb = Molecule('pdb/membrane80by80C36.pdb')\n",
    "mcenter = np.mean(memb.get('coords'),axis=0)\n",
    "memb.moveBy(pcenter-mcenter)\n",
    "mol = prot.copy()\n",
    "mol.append(memb, collisions=True)  # Append membrane and remove colliding atoms\n",
    "\n",
    "#Add ligand, previously parametrized using gaamp\n",
    "lig = Molecule('pdb/QM-min.pdb') \n",
    "lig.set('segid','L')\n",
    "lcenter = np.mean(lig.get('coords'),axis=0)\n",
    "newlcenter = [np.random.uniform(-10, 10), np.random.uniform(-10, 10),  43]\n",
    "lig.rotateBy(uniformRandomRotation(), lcenter)\n",
    "lig.moveBy(newlcenter - lcenter)\n",
    "mol.append(lig) \n",
    "\n",
    "#Add water\n",
    "coo = mol.get('coords','lipids or protein')\n",
    "m = np.min(coo,axis=0) + [0,0,-5]\n",
    "M = np.max(coo,axis=0) + [0,0,20]\n",
    "mol = solvate(mol, minmax=np.vstack((m,M)))\n",
    "\n",
    "#Build\n",
    "mol = charmm.build(mol, topo=topos, param=params, outdir=path+'/build', saltconc=0.15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.equilibration_v2 import Equilibration\n",
    "md = Equilibration()\n",
    "md.runtime = 10000000\n",
    "md.temperature = 300\n",
    "md.fb_reference = 'protein and resid 293'\n",
    "md.fb_selection = 'segname L and noh'\n",
    "md.fb_box = [-39, 10, -29, 21, 47, 50]\n",
    "md.fb_k = 5\n",
    "md.useconstantratio = True\n",
    "md.write(path+'/build',path+'/equil')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Visualize the flat bottom potential box\n",
    "mol.view('not water')\n",
    "b = VMDBox([-39, 10, -29, 21, 40, 50])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mdx = AcemdLocal()\n",
    "mdx.submit(path+'/equil')\n",
    "mdx.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from htmd.protocols.production_v5 import Production\n",
    "md = Production()\n",
    "md.runtime = 50\n",
    "md.timeunits = 'ns'\n",
    "md.temperature = 300\n",
    "md.fb_reference = 'protein and resid 293'\n",
    "md.fb_selection = 'segname L and noh'\n",
    "md.fb_k = 5\n",
    "md.fb_box = [-25, 25, -25, 25, -10, 45]\n",
    "md.write(path +'/equil','gen/s1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mol.view('not water')\n",
    "b = VMDBox([-25, 25, -25, 25, -10, 45])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}