{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Molecular dynamics protocols in HTMD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In this tutorial, we should how to equilibrate and prepare a system for productive simulations in HTMD.\n",
    "\n",
    "First, let's import HTMD:"
   ]
  },
  {
   "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",
    "config(viewer='webgl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build a sample system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a built system to perform simulations. We can easily build one with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 02:36:51,151 - htmd.molecule.readers - INFO - Using local copy for 3PTB: /home/joao/maindisk/software/repos/Acellera/htmd/htmd/data/pdb/3ptb.pdb\n",
      "2018-03-17 02:36:51,421 - htmd.molecule.molecule - WARNING - Residue insertions were detected in the Molecule. It is recommended to renumber the residues using the Molecule.renumberResidues() method.\n",
      "2018-03-17 02:36:51,457 - htmd.molecule.molecule - INFO - Removed 72 atoms. 1629 atoms remaining in the molecule.\n",
      "2018-03-17 02:36:51,555 - propka - INFO - No pdbfile provided\n",
      "2018-03-17 02:37:00,607 - htmd.builder.preparationdata - INFO - The following residues are in a non-standard state: CYS     6  A (CYX), HIS    22  A (HIE), CYS    24  A (CYX), HIS    39  A (HIP), CYS    40  A (CYX), GLU    51  A (GLH), HIS    72  A (HID), CYS   108  A (CYX), CYS   115  A (CYX), CYS   136  A (CYX), CYS   147  A (CYX), CYS   161  A (CYX), CYS   172  A (CYX), CYS   182  A (CYX), CYS   196  A (CYX), CYS   209  A (CYX)\n",
      "2018-03-17 02:37:00,661 - htmd.builder.preparationdata - WARNING - Dubious protonation state: the pKa of 4 residues is within 1.0 units of pH 7.0.\n",
      "2018-03-17 02:37:00,664 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    HIS    39  A (pKa= 7.46)\n",
      "2018-03-17 02:37:00,666 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    GLU    51  A (pKa= 7.06)\n",
      "2018-03-17 02:37:00,667 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    ASP   170  A (pKa= 6.49)\n",
      "2018-03-17 02:37:00,669 - htmd.builder.preparationdata - WARNING - Dubious protonation state:    N+      0T A (pKa= 7.49)\n",
      "2018-03-17 02:37:00,748 - htmd.builder.preparationdata - WARNING - Found N-terminus 80.7% buried (> 50.0% threshold)\n",
      "2018-03-17 02:37:00,749 - htmd.builder.preparationdata - WARNING - Found C-terminus involved in H bonds\n",
      "2018-03-17 02:37:00,810 - htmd.builder.solvate - INFO - Using water pdb file at: /data/joao/maindisk/software/repos/Acellera/htmd/htmd/builder/wat.pdb\n",
      "2018-03-17 02:37:02,115 - htmd.builder.solvate - INFO - Replicating 4 water segments, 2 by 1 by 2\n",
      "Solvating: 100%|██████████| 4/4 [00:03<00:00,  1.10it/s]\n",
      "2018-03-17 02:37:05,791 - htmd.builder.solvate - INFO - 6943 water molecules were added to the system.\n"
     ]
    }
   ],
   "source": [
    "tryp = Molecule(\"3PTB\")\n",
    "tryp.filter(\"protein\")\n",
    "tryp_op = proteinPrepare(tryp)\n",
    "tryp_seg = autoSegment(tryp_op)\n",
    "tryp_solv = solvate(tryp_seg,pad=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "With this solvated system and since HTMD is force-field agnostic, one can either build in CHARMM or Amber (for details on system building, see the corresponding tutorials):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 02:37:32,679 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-03-17 02:37:34,330 - htmd.builder.builder - INFO - 6 disulfide bonds were added\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bond between A: [serial 351 resid 24 resname CYX chain A segid P0]\n",
      "             B: [serial 571 resid 40 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1684 resid 115 resname CYX chain A segid P0]\n",
      "             B: [serial 2612 resid 182 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2147 resid 147 resname CYX chain A segid P0]\n",
      "             B: [serial 2354 resid 161 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 1605 resid 108 resname CYX chain A segid P0]\n",
      "             B: [serial 3005 resid 209 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 2495 resid 172 resname CYX chain A segid P0]\n",
      "             B: [serial 2800 resid 196 resname CYX chain A segid P0]\n",
      "\n",
      "Bond between A: [serial 92 resid 6 resname CYX chain A segid P0]\n",
      "             B: [serial 1989 resid 136 resname CYX chain A segid P0]\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 02:37:34,835 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-03-17 02:37:35,225 - htmd.builder.charmm - WARNING - Failed to set coordinates for 9 atoms.\n",
      "2018-03-17 02:37:35,227 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.\n",
      "2018-03-17 02:37:35,228 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 16 atoms.\n",
      "2018-03-17 02:37:35,229 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/build-charmm/log.txt for further information.\n",
      "2018-03-17 02:37:35,230 - htmd.builder.charmm - INFO - Finished building.\n",
      "2018-03-17 02:37:36,838 - htmd.builder.ionize - INFO - Adding 8 anions + 0 cations for neutralizing and 0 ions for the given salt concentration.\n",
      "2018-03-17 02:37:37,042 - htmd.builder.ionize - INFO - Min distance of ions from molecule: 5A\n",
      "2018-03-17 02:37:37,043 - htmd.builder.ionize - INFO - Min distance between ions: 5A\n",
      "2018-03-17 02:37:37,044 - htmd.builder.ionize - INFO - Placing 8 ions.\n",
      "2018-03-17 02:37:39,577 - htmd.builder.charmm - INFO - Writing out segments.\n",
      "2018-03-17 02:37:41,937 - htmd.builder.charmm - INFO - Starting the build.\n",
      "2018-03-17 02:37:42,333 - htmd.builder.charmm - WARNING - Failed to set coordinates for 9 atoms.\n",
      "2018-03-17 02:37:42,335 - htmd.builder.charmm - WARNING - Failed to guess coordinates for 3 atoms due to bad angles.\n",
      "2018-03-17 02:37:42,337 - htmd.builder.charmm - WARNING - Poorly guessed coordinates for 16 atoms.\n",
      "2018-03-17 02:37:42,338 - htmd.builder.charmm - WARNING - Please check /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/build-charmm/log.txt for further information.\n",
      "2018-03-17 02:37:42,339 - htmd.builder.charmm - INFO - Finished building.\n"
     ]
    }
   ],
   "source": [
    "#tryp_amber = amber.build(tryp_solv, outdir='build-amber')\n",
    "tryp_charmm = charmm.build(tryp_solv, outdir='build-charmm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Equilibration protocol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the `Equilibration` class, which is an MD protocol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.protocols.equilibration_v2 import Equilibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MD protocols are not imported automatically, the user must choose which one to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now let's start an `Equilibration` object, which already has sensible defaults for an equilibration MD simulation and let's define the remaining ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = Equilibration()\n",
    "md.runtime = 1000\n",
    "md.timeunits = 'fs'\n",
    "md.temperature = 300\n",
    "md.useconstantratio = False  # only for membrane sims\n",
    "# this is only needed for setting the flatbottom potential, otherwise remove it\n",
    "# md.fb_reference = 'protein and resid 293'\n",
    "# md.fb_selection = 'segname L and noh'\n",
    "# md.fb_box = [-25, 25, -25, 25, 43, 45]\n",
    "# md.fb_k = 5\n",
    "md.write('./build-charmm/', './equil')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can inspect what the `Equilibration` object has created with the `write` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input  parameters  \u001b[0m\u001b[38;5;34mrun.sh\u001b[0m*  structure.pdb  structure.psf\r\n"
     ]
    }
   ],
   "source": [
    "%ls equil/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Run the equilibration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's use the queue resources of HTMD to run the simulation. One can use the local computer and the `LocalGPUQueue` class to submit a job. The equilibration time is short (see above) for demonstration purposes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 02:39:07,892 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices\n",
      "2018-03-17 02:39:07,953 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3\n",
      "2018-03-17 02:39:08,112 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices\n",
      "2018-03-17 02:39:08,175 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3\n",
      "2018-03-17 02:39:08,179 - htmd.queues.localqueue - INFO - Queueing /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/equil\n",
      "2018-03-17 02:39:08,181 - htmd.queues.localqueue - INFO - Running /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/equil on device 0\n",
      "2018-03-17 02:39:33,879 - htmd.queues.localqueue - INFO - Completed /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/equil\n"
     ]
    }
   ],
   "source": [
    "local = LocalGPUQueue()\n",
    "local.submit('./equil/')\n",
    "local.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Production protocol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's use the `Production` class and do the same as before to perform a short production MD simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from htmd.protocols.production_v4 import Production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "md = Production()\n",
    "md.runtime = 1\n",
    "md.timeunits = 'ns'\n",
    "md.temperature  = 300\n",
    "md.acemd.bincoordinates = 'output.coor'\n",
    "md.acemd.extendedsystem  = 'output.xsc'\n",
    "md.write('equil','prod')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "input  input.coor  input.xsc  parameters  \u001b[0m\u001b[38;5;34mrun.sh\u001b[0m*  structure.pdb  structure.psf\r\n"
     ]
    }
   ],
   "source": [
    "% ls prod/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Run the production"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2018-03-17 02:39:56,043 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices\n",
      "2018-03-17 02:39:56,103 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3\n",
      "2018-03-17 02:39:56,440 - htmd.queues.localqueue - INFO - Trying to determine all GPU devices\n",
      "2018-03-17 02:39:56,501 - htmd.queues.localqueue - INFO - Using GPU devices 0,1,2,3\n",
      "2018-03-17 02:39:56,505 - htmd.queues.localqueue - INFO - Queueing /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/prod\n",
      "2018-03-17 02:39:56,507 - htmd.queues.localqueue - INFO - Running /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/prod on device 0\n",
      "2018-03-17 02:44:05,268 - htmd.queues.localqueue - INFO - Completed /data/joao/maindisk/software/repos/Acellera/htmd/tutorials/prod\n"
     ]
    }
   ],
   "source": [
    "local = LocalGPUQueue()\n",
    "local.submit('./prod/')\n",
    "local.wait()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the above process finished, one can check that a trajectory was produced:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prod/output.xtc\r\n"
     ]
    }
   ],
   "source": [
    "%ls prod/output.xtc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quickly visualize the trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "71be7b468e1b4343b5834edf051f5890",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "traj = Molecule('prod/structure.pdb')\n",
    "traj.read('prod/output.xtc')\n",
    "traj.wrap()\n",
    "traj.align('protein and name CA')\n",
    "traj.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Waters are not for clarity"
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}