Tutorial_QMMM_String_eABF.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run This Cell First!\n",
    "# type: Ctrl + Enter\n",
    "\n",
    "from IPython.display import HTML\n",
    "\n",
    "hide_me = ''\n",
    "HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    "  if (code_show) {\n",
    "    $('div.input').each(function(id) {\n",
    "      el = $(this).find('.cm-variable:first');\n",
    "      if (id == 0 || el.text() == 'hide_me') {\n",
    "        $(this).hide();\n",
    "      }\n",
    "    });\n",
    "    $('div.output_prompt').css('opacity', 0);\n",
    "  } else {\n",
    "    $('div.input').each(function(id) {\n",
    "      $(this).show();\n",
    "    });\n",
    "    $(\n",
    "    'div.output_prompt').css('opacity', 1);\n",
    "  }\n",
    "  code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "<form action=\"javascript:code_toggle()\"><input style=\"opacity:1\" type=\"submit\" value=\"Click here to toggle on/off the raw code.\"></form>\n",
    "\n",
    "<p>You can find the video at: https://youtu.be/2eCHD6f_phE</p>\n",
    "\n",
    "<style>\n",
    ".video-responsive{\n",
    "overflow:hidden;\n",
    "padding-bottom:56.25%;\n",
    "position:relative;\n",
    "height:0;\n",
    "}\n",
    ".video-responsive iframe{\n",
    "left:0;\n",
    "top:0;\n",
    "height:100%;\n",
    "width:100%;\n",
    "position:absolute;\n",
    "}\n",
    "</style>\n",
    "\n",
    "<div class='imgright' style='width: 80%'>\n",
    "<div class=\"video-responsive\">\n",
    "<p align=\"center\">\n",
    "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/2eCHD6f_phE\" \n",
    "frameborder=\"0\" allow=\"autoplay; encrypted-media\" allowfullscreen></iframe>\n",
    "</p>\n",
    "</div>\n",
    "</div>\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NAMD-QM/MM Tutorial: Advanced Module\n",
    "====================\n",
    "Free Energy Profile of Reaction Mechanisms\n",
    "=======================\n",
    "\n",
    "<img width=\"800\" src=\"toolkit/3ala_orbital.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will find a workflow that will help you set up, run and analyse a QM/MM SMD simulation of a reaction mechanism, a String Method optimization of the chosen reaction coordinate, and an eABF calculation of free energy along that optimized string.    \n",
    "    \n",
    "We will use a structure and parameters created with VMD following the instructions from the tutorial PDF. If you have not gone through (or can not go through) the PDF and created the initial files, example outputs can be found in the example-output folder, under \"OMPDecarb_initial\".\n",
    "    \n",
    "The notebook is divided in main section for clarity:\n",
    "\n",
    "0. Notebook setup\n",
    "1. SMD preparation and execution\n",
    "2. String Method preparation and execution\n",
    "3. String Method analysis\n",
    "4. eABF preparation and execution\n",
    "5. eABF analysis\n",
    "\n",
    "Throughout the notebook, we will attempt to explain all aspects of the methods being used, from texts and images explaining the overall processes being executed, to comments in python code being ran in the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook Setup\n",
    "====\n",
    "\n",
    "The following code block loads all packages needed for the entire notebook. Packages can be loaded in any cell, as needed, but aggregating all of them in a single cell at the beginning of the notebook is considered good practice, allows us to check if all modules have been installed correctly, and ensures they will be available for the entire notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-02T17:03:01.943395Z",
     "start_time": "2018-03-02T17:03:01.289025Z"
    }
   },
   "outputs": [],
   "source": [
    "# Python packages for file operations, regular expressions, and basic data structures\n",
    "import os, sys, shutil, glob, re, subprocess, fileinput\n",
    "from copy import copy\n",
    "from math import floor\n",
    "# Packages optimized for analysis of large datasets\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "# Packages for interactive buttons and sliders, to control interactive plots.\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import interact\n",
    "from IPython.display import display\n",
    "# Plotting packages\n",
    "import plotnine as p9\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "# Package to load and analyse MD trajectories. \n",
    "import MDAnalysis\n",
    "\n",
    "# This is necessary for the interactive plot we will use along the tutorial.\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set paths for tools used along the tutorial\n",
    "-----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-02T17:46:31.181309Z",
     "start_time": "2018-03-02T17:46:31.176945Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "catdcdPath = \"/usr/local/lib/vmd/plugins/LINUXAMD64/bin/catdcd5.2/catdcd\"\n",
    "vmdPath = \"vmd\"\n",
    "\n",
    "namdPath = \"/home/ubuntu/NAMD_Git-2018-04-03_Linux-x86_64-multicore/namd2\"\n",
    "namdPathMPI = \"/home/ubuntu/NAMD_Git-2018-04-04_Source/Linux-x86_64-g++.mpi/namd2\"\n",
    "namdPathGPU = \"/home/ubuntu/NAMD_2.12_Linux-x86_64-multicore-CUDA/namd2\"\n",
    "psfgenPath = \"/home/ubuntu/NAMD_2.12_Linux-x86_64-multicore-CUDA/psfgen\"\n",
    "\n",
    "qmBinPath = \"/home/ubuntu/orca_4012/orca\"\n",
    "\n",
    "tempDir = \"./tmp\"\n",
    "\n",
    "if True:\n",
    "    if sys.platform == 'linux':\n",
    "        catdcdPath = \"/home/melomcr/Projects_melomcr/vmd-1.9.4a12/plugins/LINUXAMD64/bin/catdcd5.2/catdcd\"\n",
    "        vmdPath = \"vmd\"\n",
    "        namdPath = \"namd2\"\n",
    "        namdPathGPU = \"/Projects/namd2/bin/latest/Linux64-multicore-CUDA/namd2\"\n",
    "        psfgenPath = \"psfgen\"\n",
    "\n",
    "    elif sys.platform == 'darwin':\n",
    "        catdcdPath = \"/Applications/VMD\\ 1.9.4.app/Contents/vmd/plugins/MACOSXX86/bin/catdcd5.2/catdcd\"\n",
    "        vmdPath = \"/Applications/VMD\\ 1.9.4.app/Contents/Resources/VMD.app/Contents/MacOS/VMD\"\n",
    "        namdPath = \"/Users/melomcr/Programs/NAMD_2.12_MacOSX-x86_64-multicore/namd2\"\n",
    "        psfgenPath = \"/Users/melomcr/Programs/NAMD_2.12_MacOSX-x86_64-multicore/psfgen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SMD preparation and execution\n",
    "========\n",
    "\n",
    "We will now prepare TCL files and call them using VMD in order to combine the files prepared in the first section of this tutorial (provided in the tutorial's PDF file). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is changing the default name of the ligand created for CGenFF's parameter file, from \"/scrat\" to \"OMP\". We will do that in the next cell with a simple command called **sed**. We will also change the PDB we created before for the new ligand as to match the new residue name.\n",
    "\n",
    "The next cell is a \"magic\" cell. Within Jupyter Notebooks, one can run individual cells, an even single lines within a cell, using a different interpreter other than Pyhton. The re are magics for bash, R, HTML, and many other languages. \n",
    "\n",
    "In the next cell, we begin by indicating the entire cell will be ran using **bash**, the most widelly used shell in linux. By using two percent signs, we indicate the entire cell will be interpreted by bash. Had we used only one percent sign at the begining of a line, only that line would be interpreted by bash."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "## For linux environments\n",
    "# sed -i 's|/scrat|OMP|g' setup/toppar/ligand.str\n",
    "# sed -i 's|BMP|OMP|g' setup/ligand.pdb\n",
    "\n",
    "## For MacOS environments\n",
    "sed -i '' 's|/scrat|OMP|g' setup/toppar/ligand.str\n",
    "sed -i '' 's|BMP|OMP|g' setup/ligand.pdb\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we must combine the new ligand structure, its parameters, and the protein chain we previously isolated from the PDB (a step carried out along the initial PDF tutorial).\n",
    "\n",
    "NAMD provides a tool called **psfgen** that is capable of performing all these steps. In the next cell, we will run psfgen with all necessary files and options. Since psfgen opens an interactive shell, the syntax below will look like we are creating a file, but we are in fact sending all lines directly to psfgen's interpreter.\n",
    "\n",
    "You will notice that the cell's magic is inicialized with an argument defined in a previous cell: the path to psfgen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash -s \"$psfgenPath\"\n",
    "\n",
    "# We move the execution of this bash shell into the setup directory.\n",
    "cd setup/\n",
    "\n",
    "# Now, we pass the following lines to psfgen\n",
    "$1 << endmol\n",
    "\n",
    "# Indicates to psfgen that the following residues and atoms must have their names changed.\n",
    "pdbalias residue HIS HSD\n",
    "pdbalias atom ILE CD1 CD\n",
    "\n",
    "# Loads necessary topology files. These are used to write the PSF file, with the correct bonds, angles, etc.\n",
    "topology toppar/top_all36_prot.rtf\n",
    "topology toppar/top_all36_cgenff.rtf\n",
    "topology toppar/toppar_water_ions_namd.str\n",
    "topology toppar/ligand.str\n",
    "\n",
    "# Creates a segment with namd  \"ENZY\" containing the protein chain.\n",
    "segment ENZY {\n",
    "        auto angles dihedrals\n",
    "        pdb enzyme.pdb\n",
    "}\n",
    "\n",
    "# Creates a segment named \"OMP\" containing the new ligand.\n",
    "segment OMP {\n",
    "        auto angles dihedrals\n",
    "        pdb ligand.pdb\n",
    "}\n",
    "\n",
    "# Creates all angles and dihedrals for the final PSF file, combining all segments.\n",
    "regenerate angles dihedrals\n",
    "\n",
    "# Indicates which files will be used to determine atom coordinates for each segment.\n",
    "coordpdb enzyme.pdb  ENZY\n",
    "coordpdb ligand.pdb  OMP\n",
    "\n",
    "# Writes the final PSF and PDB files for the full system: protein + new ligand.\n",
    "writepsf system.psf\n",
    "guesscoord\n",
    "writepdb system.pdb \n",
    "\n",
    "# Closes the pipe that sends commands to psfgen\n",
    "endmol\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our new system is complete, we must solvate and ionize it. The following cell creates the necessary TCL files for both processes, and runs VMD in text mode to carry out each action."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$vmdPath\"\n",
    "\n",
    "# Move into setup folder, where the system will be built and main files stored.\n",
    "cd setup\n",
    "\n",
    "# BASH syntax for creating a file called \"solvate_ionize.tcl\". The file will have all \n",
    "#     the following lines, untill a line containing \"EOF\" is found.\n",
    "cat <<EOF >solvate_ionize.tcl\n",
    "\n",
    "### this function translates the system coordinates\n",
    "proc moveby {sel offset} {\n",
    "    foreach coord [\\$sel get {x y z}] {\n",
    "        lappend newcoords [vecadd \\$coord \\$offset]\n",
    "    }\n",
    "    \\$sel set {x y z} \\$newcoords\n",
    "}\n",
    "#################################################\n",
    "\n",
    "## We first center the system before adding a box of water arround it.\n",
    "\n",
    "# Loads system coordinates\n",
    "mol new system.psf\n",
    "mol addfile system.pdb\n",
    "# Selects all atoms in the system\n",
    "set sel [atomselect 0 \"all\"]\n",
    "# Measures the geometric center of all atoms, and multiplies the \n",
    "#  resulting vector by -1. This is used as the offset vector to\n",
    "#  move all atoms in the system.\n",
    "set vec [ vecscale -1 [measure center \\$sel]]\n",
    "# Centers the entire system in the origin of the coordinate axis.\n",
    "moveby \\$sel \\$vec\n",
    "# Over-Writes PDB with centered system.\n",
    "\\$sel writepdb system.pdb\n",
    "\n",
    "## solvating the box\n",
    "# Loads the required VMD package and creates a water box with 80 Angstroms of side.\n",
    "package require solvate\n",
    "solvate system.psf system.pdb  -minmax {{-40.0 -40.0 -40.0 } {40.0 40.0 40.0}} -o system.solv\n",
    "\n",
    "## minimun sovation box\n",
    "# Alternatively, a minumum box can be created with the command below, but it will \n",
    "#   not (necessarily) create a water cube.\n",
    "#solvate tmp.psf tmp.pdb  -t 5 -o tmp.solv\n",
    "\n",
    "## Adding Ions to the system.\n",
    "# Loads VMD package\n",
    "package require autoionize\n",
    "autoionize -psf system.solv.psf -pdb system.solv.pdb -neutralize -o system.solv.ionized\n",
    "\n",
    "# End Of File\n",
    "EOF\n",
    "\n",
    "# Runs VMD in text mode, and executes the script we just created.\n",
    "eval \"$1\" -dispdev text -e solvate_ionize.tcl\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img width=\"450\" src=\"toolkit/systemWaterBox.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare the classical equilibraiton\n",
    "-------\n",
    "\n",
    "Now that we have the complex we wish to study, we must equilibrate the system. This will be done classically, for a short period of time, to keep the tutorial short. Ideally, after a long classical equillibration, a QM/MM equilibration should also be ran.\n",
    "\n",
    "(If you did NOT execute the previous steps, or if you just wish to check your output, previously-ran results can be fount in the *example-output* folder, under *OMPDecarb_initial*. Copy the contents of this folder to the \"setup\" folder in the main tutorial directory, with the exception of the file \"ligand.str\", which should be placed under the \"setup/toppar\" folder.)\n",
    "\n",
    "First, we copy the structure (PSF) file, and create the *extended system* (XSC) file (which allows the use of periodic boudary conditions). The latter indicates the dimensions of the periodic box, and is created based on the limits of the solvation box created in previous cells of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "\n",
    "# Move into the folder where we will run the classical equillibration.\n",
    "cd run.mm/\n",
    "\n",
    "# Copy the PSF file from our final system.\n",
    "cp ../setup/system.solv.ionized.psf decarboxylase.0.psf\n",
    "\n",
    "# Create the extended system file for periodic boundary conditions.\n",
    "#   We are using the same box dimentions we used to create the water box, with\n",
    "#   an extra 0.5 A of space.\n",
    "cat <<EOF >decarboxylase.0.xsc\n",
    "# NAMD extended system configuration restart file\n",
    "#\\$LABELS step a_x a_y a_z b_x b_y b_z c_x c_y c_z o_x o_y o_z s_x s_y s_z s_u s_v s_w\n",
    "0 80.50 0 0 0 80.50 0 0 0 80.50 0 0 0 2.1166e-06 2.1166e-06 2.1166e-06 0 0 0\n",
    "EOF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are using parameters for our ligand which are \"less than perfect\", we will constrain most of the ligand molecule, except the atoms we added or modified. We will also constrain the backbone of the protein chain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$vmdPath\"\n",
    "\n",
    "# Move into the folder where we will run the classical equillibration.\n",
    "cd run.mm/\n",
    "\n",
    "cat <<EOF >prepare.tcl\n",
    "\n",
    "# For this TCL script, we load the system we just prepared\n",
    "mol new ../setup/system.solv.ionized.psf\n",
    "mol addfile ../setup/system.solv.ionized.pdb\n",
    "\n",
    "# Now we use the two fields available in PDB files, \"beta\" and \"occupancy\", \n",
    "#   and set tehir values to zero for all atoms.\n",
    "[atomselect 0 all] set beta 0.0\n",
    "[atomselect 0 all] set occupancy 0.0\n",
    "\n",
    "# The \"beta\" field will be used to indicate which atoms will have their \n",
    "#    movements contrained during the equilibration.\n",
    "[atomselect 0 \"backbone\" ] set beta 1.0\n",
    "[atomselect 0 \"resname OMP and not hydrogen and not name OA OB C7\" ] set beta 1.0\n",
    "\n",
    "# Now we write the constraints file:\n",
    "set sel [atomselect 0 all]\n",
    "\\$sel writepdb decarboxylase.0.pdb\n",
    "\n",
    "# And exit VMD\n",
    "quit\n",
    "\n",
    "EOF\n",
    "\n",
    "eval \"$1\" -dispdev text -e prepare.tcl\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the Classical Equillibration\n",
    "--------\n",
    "\n",
    "Now that the contrained equillibration is set up, we will execute it!\n",
    "\n",
    "(WARNING: this may take several hours (~5.5 hours in amazon's g2.2xlarge node). For the purpose of this tutorial, we have provided pre-ran results in the \"example-output\" folder, under \"run.mm\". Just copy its contents into the *run.mm* folder in the main tutorial folder.)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$namdPathGPU\"\n",
    "\n",
    "cd run.mm/\n",
    "\n",
    "#eval $1 decarboxylase.conf > output.log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare SMD - Step 1: QM/MM\n",
    "-------\n",
    "\n",
    "In order to select a small QM region for this chemical reaction, we will only use parts of two residues from the protein chain, and a part of our ligand molecule. As covered in the basic QM/MM tutorial, while making a \"manual\"selection of QM atoms, it is very important to make sure that the total charge of the QM region is zero. Otherwise, fractional charges would be added to the system, wich would cause serious problems for PME calculations.\n",
    "\n",
    "First, we create the PDB with the selection of atoms for our QM region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$vmdPath\"\n",
    "\n",
    "# Move into the folder where we will run the classical equillibration.\n",
    "cd run.smd/\n",
    "\n",
    "# Copy the PSF file from our final system.\n",
    "cp ../setup/system.solv.ionized.psf decarboxylase.0.psf\n",
    "\n",
    "# Copy the final XSC and COOR files from our classical equillibration.\n",
    "cp ../run.mm/decarboxylase.1r.xsc decarboxylase.0r.xsc\n",
    "cp ../run.mm/decarboxylase.1r.coor decarboxylase.0r.coor\n",
    "\n",
    "# We now create the TCL file that will prepare the QM region.\n",
    "cat <<EOF >prepare.tcl\n",
    "\n",
    "# Load the structure from the results of the classical equillibration.\n",
    "# This isn't strictly necessary since we will only be using this script to create a PDB\n",
    "#    file that will mark the QM region. Atom positions will be taken from \n",
    "#    the decarboxylase.0r.coor file copied above.\n",
    "mol new decarboxylase.0.psf\n",
    "mol addfile ../run.mm/decarboxylase.1r.coor type namdbin\n",
    "\n",
    "# We will once more set all beta and occupanciy fields to zero.\n",
    "# The QM region will be defined by setting the beta value to 1 for all QM atoms.\n",
    "# Since we will create QM-MM bonds (where one atom of the bond is treated classically \n",
    "#   and the other is treated quantum mechanically), we will also use the occupancy column,\n",
    "#   marking BOTH atoms that participate in a QM-MM bond.\n",
    "[atomselect 0 all] set beta 0.0\n",
    "[atomselect 0 all] set occupancy 0.0\n",
    "[atomselect 0 \"resname OMP and name N1 C2 O2 N3 C4 O4 C5 C6 C7 OA OB H H2\" ] set beta 1.0\n",
    "[atomselect 0 \"resname OMP and name N1 C1'\" ] set occupancy 1.0\n",
    "[atomselect 0 \"resname ASP and resid 70 and sidechain\" ] set beta 1.0\n",
    "[atomselect 0 \"resname ASP and resid 70 and name CA CB\" ] set occupancy 1.0\n",
    "[atomselect 0 \"resname LYS and resid 72 and sidechain\" ] set beta 1.0\n",
    "[atomselect 0 \"resname LYS and resid 72 and name CA CB\" ] set occupancy 1.0\n",
    "\n",
    "# We load the topotools package to guess element names based on the atom's mass.\n",
    "# We do this for all atoms so the PDB file contain their elements, which will be used \n",
    "#   by QM packages in their calculations.\n",
    "package require topotools\n",
    "topo guessatom element mass\n",
    "\n",
    "# Write all atoms in a new PDB file for the entire system.\n",
    "set sel [atomselect 0 all]\n",
    "\\$sel writepdb decarboxylase.0.pdb\n",
    "\n",
    "# Exit VMD\n",
    "quit\n",
    "EOF\n",
    "\n",
    "# Run VMD with the TCL\n",
    "eval \"$1\" -dispdev text -e prepare.tcl\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Edit Charges\n",
    "\n",
    "Now, we will need to edit the charges from the OMP molecule (the ligand) in order to be able to select only part of it, without creating fractional charges.\n",
    "We do this by changing the charges of two atoms of our molecule using VMD. Our case is so simple that this modification could easily be done \"by hand\", by editing the PSF file directly, but we will make use of a more flexible strategy here, so it can be adapted for more complex cases in the future."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$vmdPath\"\n",
    "\n",
    "# Move into the folder where we will run the classical equillibration.\n",
    "cd run.smd/\n",
    "\n",
    "# We now create the TCL file that will prepare the QM region.\n",
    "cat <<EOF >edit_charges.tcl\n",
    "\n",
    "mol new decarboxylase.0.psf\n",
    "mol addfile decarboxylase.0.pdb\n",
    "\n",
    "# Selects all the OMP atoms in the QM region.\n",
    "set sel  [atomselect 0 \"resname OMP and beta 1\"]\n",
    "\n",
    "# Select the atoms that will have their charge slightly shifted \n",
    "set selN [atomselect 0 \"resname OMP and name N1\"]\n",
    "set selC [atomselect 0 \"resname OMP and name C1'\"]\n",
    "\n",
    "\\$selC set charge  0.130\n",
    "\\$selN set charge -0.459\n",
    "\n",
    "set all [atomselect 0 all]\n",
    "\\$all writepsf decarboxylase.0.psf\n",
    "\n",
    "# Exit VMD\n",
    "quit\n",
    "EOF\n",
    "\n",
    "# Run VMD with the TCL\n",
    "eval \"$1\" -dispdev text -e edit_charges.tcl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare SMD Step - 2: Constraints\n",
    "-------\n",
    "\n",
    "Using the results of the classical equillibration, we can setup and start the QM/MM SMD simulation. This SMD simulation uses a slight variation of the classical SMD, called \"conditional-SMD\" (cSMD), which was optimized to induce chemical reactions during QM/MM MD simulations. In this case, one atom is directed towards another during the simulations, much like the traditional SMD, but the direction in which it is pulled is constantly updated according to the movement of the targeted atom, and the force is only applied wile both atoms are separated by a cutoff. Once atoms are closer than the cutoff (which shoudl be slightly larger than the expected bond length), the force stops being applied. For reactions involving multiple bond formation this becomes especially relevant, since some bonds may form before others, and the conditional nature of cSMD will make the atoms stop being \"pushed\" once they form the bonds they are being induced to form.\n",
    "\n",
    "The same could be accomplished using colvars, but this way provided a faster and easier interface for common bio-molecular systems. All that is needed is a text file in the following format:\n",
    "\n",
    "Atom1 Atom2 Force(Kcal/Mol/A) Speed(A/step) Cutoff(A)\n",
    "\n",
    "In the following cells, we will create the necessary files.\n",
    "\n",
    "<img width=\"450\" src=\"toolkit/Colvars_OMP-01.png\" />\n",
    "\n",
    "It is important to point out that we will NOT separate the carbons C6 and C7 of OMP, but rather orient the key groups to create the perfect scenario for the reaction to occur. We will 1) steer the oxygen atoms of the carboxylic group on OMP and the aspartic acid (resid 70) towards one another, and 2) steer the hydrogen atom of the lysine (resid 72) towards the carbon C6 of OMP. This mechanism was suggested by Van Vleet, *et al* [1].\n",
    "\n",
    "[1]: https://pubs.acs.org/doi/abs/10.1021/bi701664n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$vmdPath\"\n",
    "\n",
    "# Move into the folder where we will run the classical equillibration.\n",
    "cd run.smd/\n",
    "\n",
    "# Create \"reaction_1.rct\" file, which will prepare the system for the \n",
    "#    reaction, bringing together the relevant chemical groups.\n",
    "cat <<EOF >reaction_1.rct\n",
    "#### OMP CO2 <-> ASP O2\n",
    "3304  980 100  0.0002  2.5\n",
    "3303  979 100  0.0002  2.5\n",
    "#\n",
    "#### NH3 LYS <-> OMP C6\n",
    "1019 3301 100  0.0002  2.5\n",
    "EOF\n",
    "\n",
    "\n",
    "# Create \"reaction_2.rct\" file, which will induce the chemical reaction.\n",
    "cat <<EOF >reaction_2.rct\n",
    "#### OMP CO2 <-> ASP O2\n",
    "3304  980 100  0.001  2.0\n",
    "3303  979 100  0.001  2.0\n",
    "#\n",
    "####   H LYS <-> OMP C6\n",
    "1020 3301 100  0.001  1.2\n",
    "EOF\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We copy the configuration files from the template directory so that we can update the path to the QM package. This strategy can be expanded to automatically change other relevant values for the simulation. For the tutorial, the only variable we set is the path to the QM package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create directory in RAM for faster QM calculation and I/O during SMD simulations.\n",
    "# This path is used in the \"decarboxylase.*.conf\" files.\n",
    "if sys.platform == 'linux':\n",
    "    os.makedirs(\"/dev/shm/NAMD\", exist_ok=True)\n",
    "\n",
    "# Get configuration files for SMD calculations.\n",
    "for filename in glob.glob('./templates_smd/*'):\n",
    "    shutil.copy(filename, \"./run.smd/\")\n",
    "    \n",
    "    with fileinput.FileInput(\"./run.smd/\" + os.path.basename(filename), inplace=True) as file:\n",
    "        for line in file:\n",
    "            print(line.replace(\"QMBINPATH\", qmBinPath), end='')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the SMD steps\n",
    "----------\n",
    "\n",
    "The following cell is a *bash* cell that can be used to run NAMD. Remember that this process can take **several hours**, and for the purposes of the QM/MM Workshop, we have provided pre-ran files in the *example-output/run.smd* folder. You can simply copy the pre-ran files to the *run.smd* folder of the main tutorial folder.\n",
    "\n",
    "If you still wish to proceed, just remove the comment (the *hash* character \"#\") from the lines bellow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$namdPath\"\n",
    "\n",
    "cd run.smd/\n",
    "\n",
    "#eval $1 decarboxylase.1.conf > output.1.log\n",
    "\n",
    "#eval $1 decarboxylase.2.conf > output.2.log\n",
    "\n",
    "#eval $1 decarboxylase.3.conf > output.3.log\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine SMD results into unique DCD and XST files\n",
    "--------\n",
    "\n",
    "We now combine all the trajectory output for the 3 SMD steps we performed so they can be loaded at once. We also read only 1 out of every 100 steps (this is controlled by the *stride* option in catdcd), to remove small fluctuations, allowing us to handle less data.\n",
    "\n",
    "If you chose the copy the pre-ran results instead of executing the *previous* cell, there is no need to execute the *following* cell, since the pre-ran results are already combined and reduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$catdcdPath\"\n",
    "\n",
    "cd run.smd/\n",
    "\n",
    "#eval \"$1\" -o smd.dcd -otype dcd -stype psf -s decarboxylase.0.psf -stride 10 -dcd decarboxylase.1.dcd -dcd decarboxylase.2.dcd -dcd decarboxylase.3.dcd \n",
    "\n",
    "# Now we combine the trajectory XST file, to get the box size along the simulation. This is only really needed if \n",
    "#    the periodic cell was not constant (as in NPT simulations), but does not \"hurt\" in NVT simulations, so we \n",
    "#    provide the code to here keep the workflow flexible. \n",
    "\n",
    "# Gets only the data, without the header.\n",
    "cat decarboxylase.1.xst | grep -v ^# >  smd_full.xst\n",
    "cat decarboxylase.2.xst | grep -v ^# >> smd_full.xst\n",
    "cat decarboxylase.3.xst | grep -v ^# >> smd_full.xst\n",
    "\n",
    "# Add header to combined file:\n",
    "#echo '# NAMD extended system trajectory file' > smd.xst\n",
    "#echo '#$LABELS step a_x a_y a_z b_x b_y b_z c_x c_y c_z o_x o_y o_z' >> smd.xst\n",
    "\n",
    "# Gets only one out of every 100 lines, like the DCD stride above.\n",
    "#sed -n '1p;0~10p' smd_full.xst >> smd.xst\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# String Method preparation and execution\n",
    "\n",
    "Set input files, options, and atoms for colvars\n",
    "--------\n",
    "\n",
    "The main options, files and paths used throughout the notebook are set in the following cell.\n",
    "All of the files will be *created* along the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the PDB and PSF files for future use along the notebook.\n",
    "pdbFile = \"./run.smd/decarboxylase.0.pdb\"\n",
    "psfFile = \"./run.smd/decarboxylase.0.psf\"\n",
    "\n",
    "# Define the location of the combined SMD data.\n",
    "dcdFile = \"./run.smd/smd.dcd\"\n",
    "xstFile = \"./run.smd/smd.xst\"\n",
    "\n",
    "## Calculation options\n",
    "\n",
    "# Location of String Method setup and output folders.\n",
    "basePath = \"./String_eABF\"\n",
    "setupPath = basePath + \"/setup/\"\n",
    "outputPath = basePath + \"/output/\"\n",
    "\n",
    "# Indicates if the string replicas will be executed in parallel or sequential form.\n",
    "parallelExec = False\n",
    "sequentialExec = True\n",
    "\n",
    "# Number of iterations for the String Method\n",
    "stringNumIter = 5\n",
    "# Number of walkers (replicas) per image.\n",
    "stringNumWalkers = 5\n",
    "\n",
    "# Number of eABF walkers per image.\n",
    "eABFWalkersPerImg = 1\n",
    "\n",
    "# Number of eABF steps per walker.\n",
    "eABFnumSteps = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we define the collective variables (colvar) that will be used to describe the reaction coordinate. Each colvar will be described by a distance between a pair of relevant atoms. They are indicated in the image below by the dashed red lines.\n",
    "\n",
    "<img width=\"450\" src=\"toolkit/Colvars_OMP-02.png\" />\n",
    "\n",
    "The pairs of atoms should be defined in a data structure that takes three elements: the index of the first atom, as indicated in VMD, the index of the second atom, and a label for the colvar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Decarboxylase\n",
    "\n",
    "# We will create a list of colvars.\n",
    "atomDistList = []\n",
    "\n",
    "# This list should be built with VMD's 0-based atom index.\n",
    "atomDistList.append( (3301,3302,\"OMP C6 <-> OMP C7\") )\n",
    "atomDistList.append( (3301,1020,\"OMP C6 <-> Lys62 HZ1\") )\n",
    "atomDistList.append( (1019,1020,\"Lys61 NZ <-> Lys62 HZ1\") )\n",
    "atomDistList.append( (3304,980,\"OMP OA <-> Asp60 OD2\") )\n",
    "atomDistList.append( (3303,979,\"OMP OB <-> Asp60 OD1\") )\n",
    "\n",
    "atomDistList"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Calculate colvar values during SMD trajectory\n",
    "----\n",
    "We will now use the MDAnalysis module to calculate the values of collective variables along the SMD trajectory. As we direct atoms' movements to induce the chemical reaction, the relevant distances between atom pairs will oscilate and vibrate during the MD simulation.\n",
    "\n",
    "#### Working with large volumes of data: Pandas\n",
    "\n",
    "In interpreted languages like Python and R, it is important to use packages that provide efficient ways to load and process large amounts of data. In Python, the package *Pandas* provides functions and data structures that are highly optimized. It'scounterpart in R is the package *data.tables*.\n",
    "\n",
    "In this notebook we use Pandas *Data Frames*, an optimized data structure for multi-dimensional datasets, so that these tools become familiar and integrated into the workflow, even though the small amounts of data created in this tutorial do not call for such optimizations.\n",
    "\n",
    "#### Re-arranging data for faster and easier statistics calculation and plotting:  Melt/Cast\n",
    "\n",
    "Another advantage of using Pandas (or similar packages) is the integration and availability of tools that **Melt** and **Cast** datasets. \n",
    "\n",
    "Usual datasets are created in a \"wide\" format, where multiple columns describe values associated with a series of observations (in rows). One can *melt* the data fromat a \"wide\" into a \"long\" form which distinguishes measured and identifying variables. The data in its new format can then be ‘cast’ into a new shape, whether it be a data frame, list, or high dimensional array ([Some literature][1]).\n",
    "\n",
    "These concepts and tools make it very easy to analyse and process data, calculate statistics, and create plots and visualizations to help interrogate the data.\n",
    "\n",
    "[1]: https://www.jstatsoft.org/article/view/v021i12 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-02T17:46:32.120047Z",
     "start_time": "2018-03-02T17:46:32.008501Z"
    }
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "if not os.path.exists(tempDir):\n",
    "    os.mkdir(tempDir)\n",
    "\n",
    "numColvars = len(atomDistList)\n",
    "keyNameDic = { \"dist\"+str(indx) : str(atomDistList[indx][2]) for indx in range(numColvars) }\n",
    "    \n",
    "# Open the PSF and SMD DCD files\n",
    "u = MDAnalysis.Universe(psfFile,dcdFile)\n",
    "\n",
    "# Creates atom selection objects\n",
    "selPairList = []\n",
    "for cvSel in atomDistList:\n",
    "    selPairList.append((u.atoms[cvSel[0]], u.atoms[cvSel[1]] ) )\n",
    "\n",
    "# Initialize the data-frame with standard column names\n",
    "cvSMD = pd.DataFrame(columns=[ \"dist\"+str(i) for i in range(numColvars) ])\n",
    "\n",
    "def calcCV(ts):\n",
    "    \n",
    "    return [ np.linalg.norm( selPairList[cvIndx][0].position - selPairList[cvIndx][1].position ) \\\n",
    "         for cvIndx in range(numColvars) ]\n",
    "\n",
    "\n",
    "# Process trajectory\n",
    "cvSMD = pd.DataFrame([calcCV(ts) for ts in u.trajectory])\n",
    "cvSMD.rename(index=str, columns={ cvIndx:\"dist\" + str(cvIndx) for cvIndx in range(numColvars) }, inplace=True)\n",
    "cvSMD.reset_index(drop=True, inplace=True)\n",
    "\n",
    "numFrames = len(cvSMD)\n",
    "\n",
    "cvSMDMelt = pd.melt(cvSMD, var_name=\"colvar\", value_name=\"Distance\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select initial string images\n",
    "\n",
    "In order to initiate the string optimization, we must choose frames from the trajectory that represent significant steps during the reaction, keeping a good resolution of the chemical process being studied. This process can be tedious and repetitive if done manually, particularly for transformations that can be abrupt and rare. \n",
    "\n",
    "To streamline the selection of frames from the SMD simulation, we created a graphical interface and tools (in the cell below) to help analyze colvar values along the trajectory, and select the best frames to initialize the string optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Attention**\n",
    "For simplicity, the code inside the following cell has been hidden, but you can select it using your keyboard arrow keys, clicking the Down Arrow key *Once*, and then *run* it like any other cell, by typing *Ctrl + Enter*. It may take a few seconds to load."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "hide_me\n",
    "%matplotlib notebook\n",
    "\n",
    "class imageDat:\n",
    "    ''' Holds data on selected images for String Optimization'''\n",
    "    \n",
    "    frame = 0\n",
    "    ftype = 0\n",
    "    selected = False\n",
    "    \n",
    "    def __init__(self, frame=0, ftype=0, selected=False):\n",
    "        self.frame = frame\n",
    "        self.ftype = ftype\n",
    "        self.selected = selected\n",
    "\n",
    "allImages = []\n",
    "numImgs = len(allImages)   # Will initialize as zero\n",
    "defaultNumImgs = 10\n",
    "\n",
    "# Sets widgets to select number of images\n",
    "imagesText = widgets.BoundedIntText(\n",
    "    value=defaultNumImgs, step=1, min=3, max=40, \n",
    "    description=\"Initial Images:\",\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "imagesSlider = widgets.IntSlider(\n",
    "    value=numImgs, step=1, readout=False, min=3, max=40\n",
    ")\n",
    "offsetText = widgets.BoundedIntText(\n",
    "    value=numImgs, step=1, min=0, max=numFrames,\n",
    "    description=\"Offset:\",\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "resetSelB = widgets.Button(\n",
    "    value=False,\n",
    "    description='Reset Selection',\n",
    "    disabled=False,\n",
    "    button_style='danger', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Resets all selections and allows a new set of evenly spaced images to be created.',\n",
    "    icon='eraser'\n",
    ")\n",
    "\n",
    "mylink = widgets.jslink((imagesText, 'value'), (imagesSlider, 'value'))\n",
    "\n",
    "hboxImgs = widgets.HBox([imagesText,imagesSlider,offsetText,resetSelB])\n",
    "\n",
    "# Sets widgets to automatically create images in high variability regions.\n",
    "activateInterTB = widgets.ToggleButton(\n",
    "    value=False,\n",
    "    description='Add Intermediates',\n",
    "    disabled=False,\n",
    "    button_style='info', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Adds intermediate images in regions of high variability. \\\n",
    "This will lock the number of evenly spaced images selected above.',\n",
    "    icon='square' # \"check-square\" vs \"square\"\n",
    ")\n",
    "defaultMaxDelta = 0.3\n",
    "maxDeltaFT = widgets.BoundedFloatText(\n",
    "    value=defaultMaxDelta, min=0, \n",
    "    description=\"Maximum delta:\",\n",
    "    style={'description_width': 'initial'},\n",
    "    disabled=True\n",
    ")\n",
    "addLastTB = widgets.ToggleButton(\n",
    "    value=True,\n",
    "    description='Add last frame',\n",
    "    style={'description_width': 'initial'},\n",
    "    disabled=True,\n",
    "    button_style='info', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Add an image at the last step of the simulation.',\n",
    "    icon='check-square' # \"check-square\" vs \"square\"\n",
    ")\n",
    "calcInterB = widgets.Button(\n",
    "    value=False,\n",
    "    description='calculate',\n",
    "    disabled=True,\n",
    "    button_style='warning', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='(Re)Calculate intermediate images. This may take some time.',\n",
    "    icon='retweet'\n",
    ")\n",
    "\n",
    "hboxInter = widgets.HBox([activateInterTB, maxDeltaFT, addLastTB, calcInterB])\n",
    "\n",
    "# Sets widgets to indicate total number of images\n",
    "totalImagesText = widgets.Text(\n",
    "    value='0',\n",
    "    description='Total Images:',\n",
    "    disabled=True\n",
    ")\n",
    "\n",
    "imgSelect = widgets.Dropdown(\n",
    "    options=[('Zero', 0)],\n",
    "    value=0,\n",
    "    description='Image:',\n",
    ")\n",
    "\n",
    "imgChgFrameText = widgets.BoundedIntText(\n",
    "    value=0, min=0, max=numFrames,\n",
    "    description=\"Frame:\",\n",
    "    style={'description_width': 'initial'}\n",
    ")\n",
    "\n",
    "addImg = widgets.Button(\n",
    "    description='',\n",
    "    disabled=False,\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Add Image',\n",
    "    icon='plus'\n",
    ")\n",
    "\n",
    "rmImg = widgets.Button(\n",
    "    description='',\n",
    "    disabled=False,\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Remove Image',\n",
    "    icon='minus'\n",
    ")\n",
    "\n",
    "hboxFT = widgets.HBox([totalImagesText, addImg, rmImg, imgSelect, imgChgFrameText])\n",
    "\n",
    "def resetImages(imgCount = defaultNumImgs):\n",
    "    '''Creates main images. Erases all prior image data.'''\n",
    "    global allImages, numImgs\n",
    "    \n",
    "    if imgCount > 40:\n",
    "        raise Exception('Too many images!')\n",
    "    \n",
    "    offset = offsetText.value\n",
    "    \n",
    "    # Initializes number of frames between images as one tenth of total frames, rounding down, and \n",
    "    # making sure we have at least 1 (one) frame between images\n",
    "    imgStride = max(floor( (numFrames-offset)/imgCount ),1)\n",
    "    allImages = [ imageDat(frame=(x*imgStride) + offset) for x in range(imgCount)]\n",
    "    numImgs = len(allImages)\n",
    "    \n",
    "resetImages()\n",
    "\n",
    "def toggle_addLast(*args):\n",
    "    if addLastTB.value:\n",
    "        addLastTB.icon = \"check-square\"\n",
    "    else:\n",
    "        addLastTB.icon = \"square\"\n",
    "\n",
    "addLastTB.observe(toggle_addLast,\"value\")\n",
    "\n",
    "\n",
    "# Displays all widgets\n",
    "ui = widgets.VBox([hboxImgs, hboxInter, hboxFT])\n",
    "display(ui)\n",
    "\n",
    "# Build basic plot\n",
    "xrange = np.arange(0, numFrames, 1)\n",
    "\n",
    "fig = plt.figure()\n",
    "gs = gridspec.GridSpec(2, 2, width_ratios=[9, 1], height_ratios=[7,1])\n",
    "ax = plt.subplot(gs[0,0])\n",
    "ax1 = plt.subplot(gs[0,1])\n",
    "lax = plt.subplot(gs[1,:])\n",
    "\n",
    "ax.set(xlabel='Frame', ylabel='Distance (A)')\n",
    "ax1.set(xlabel='', ylabel='Max. Variation (A)')\n",
    "ax1.set_ylim([0,1.5])\n",
    "\n",
    "plt.show()\n",
    "\n",
    "colors = ['b','g','r','c','m']\n",
    "styles = ['-','-.','--',':']\n",
    "styleTupls = []\n",
    "\n",
    "def getStyleTupl(i):\n",
    "    cIndx = i%len(colors)\n",
    "    sIndx = int(i/len(colors))\n",
    "    return (colors[cIndx],styles[sIndx])\n",
    "\n",
    "deltaLines = []\n",
    "# Populate with selected colvar values\n",
    "for key, grp in cvSMDMelt.groupby(['colvar']):\n",
    "    colStyl = getStyleTupl(len(styleTupls))\n",
    "    styleTupls.append( colStyl )\n",
    "    ax = grp.plot(ax=ax, kind='line', x=xrange, y='Distance', label=keyNameDic[key], \n",
    "                  color=colStyl[0], style=colStyl[1])\n",
    "    \n",
    "    # Initializes points used for CV deltas\n",
    "    deltaLines.append(ax1.axhline(y=0, linewidth=3, \n",
    "                                  label=key, color=colStyl[0], linestyle=colStyl[1]) )\n",
    "    \n",
    "    plt.draw()\n",
    "\n",
    "ax.set_xlim([-2,numFrames+2])\n",
    "\n",
    "h,l = ax.get_legend_handles_labels()\n",
    "lax.legend(h,l, borderaxespad=0, ncol=3)\n",
    "lax.axis(\"off\")\n",
    "\n",
    "ax.legend().set_visible(False)\n",
    "ax1.legend().set_visible(False)\n",
    "\n",
    "linesList = []\n",
    "\n",
    "def setMainImages(*args):\n",
    "    resetImages(imgCount= imagesText.value)\n",
    "    update_plot()\n",
    "\n",
    "imagesText.observe(setMainImages, 'value')\n",
    "\n",
    "def update_TotImgs(*args):\n",
    "    global allImages\n",
    "    totalImagesText.value = str(len(allImages))\n",
    "\n",
    "imagesText.observe(update_TotImgs, 'value')\n",
    "\n",
    "# Determines the maximum variation in ColVar values between sequential images.\n",
    "def find_max_delta(allImages):\n",
    "    \n",
    "    numImgs = len(allImages)\n",
    "    \n",
    "    maxDeltaCV = [0]*len(cvSMD.columns)\n",
    "    maxDeltaFrames = [ [0,0] for i in range(len(cvSMD.columns)) ]\n",
    "    \n",
    "    for indx in range(1,numImgs):\n",
    "        \n",
    "        frameFin = allImages[indx].frame\n",
    "        frameIni = allImages[indx -1].frame\n",
    "        \n",
    "        for colIndx in range(len(cvSMD.columns)):\n",
    "            \n",
    "            delta = abs( cvSMD.iloc[frameFin,colIndx] - cvSMD.iloc[frameIni,colIndx] )\n",
    "            \n",
    "            if maxDeltaCV[colIndx] < delta:\n",
    "                maxDeltaCV[colIndx] = delta\n",
    "                maxDeltaFrames[colIndx][0] = frameIni\n",
    "                maxDeltaFrames[colIndx][1] = frameFin\n",
    "    \n",
    "    return(maxDeltaCV,maxDeltaFrames)\n",
    "\n",
    "def toggle_reset(*args):\n",
    "    \n",
    "    # Resets interface\n",
    "    imagesText.disabled = False\n",
    "    imagesSlider.disabled = False\n",
    "    offsetText.disabled = False\n",
    "    \n",
    "    activateInterTB.value = False\n",
    "    toggle_intermediate()\n",
    "    \n",
    "    # Resets user input data\n",
    "    maxDeltaFT.value = defaultMaxDelta\n",
    "    imagesText.value = defaultNumImgs\n",
    "    \n",
    "    # Resets image selection\n",
    "    resetImages()\n",
    "    \n",
    "    # updates plot\n",
    "    update_plot()\n",
    "    update_TotImgs()\n",
    "    \n",
    "resetSelB.on_click(toggle_reset)\n",
    "\n",
    "def update_plot():\n",
    "    \n",
    "    global allImages,numImgs, maxDeltaFrames\n",
    "    \n",
    "    # Initialize Images on the plot\n",
    "    \n",
    "    lineIndx = 0\n",
    "    for imageObj in allImages:\n",
    "        \n",
    "        # We determine the line color\n",
    "        color = ''\n",
    "        # Red for main images (equally spaced)\n",
    "        if imageObj.ftype == 0:\n",
    "            color = 'r'\n",
    "        elif imageObj.ftype == 1:\n",
    "            color = 'g'\n",
    "        # Blue for interspaced and last frame (if requested)\n",
    "        else:\n",
    "            color = 'b'\n",
    "        \n",
    "        if imageObj.selected:\n",
    "            color = 'gold'\n",
    "        \n",
    "        # If we can use previously created vertical line objects for the plot, update their values\n",
    "        if lineIndx < len(linesList):\n",
    "            linesList[lineIndx].set_xdata(imageObj.frame)\n",
    "            linesList[lineIndx].set_color(color)\n",
    "        else:\n",
    "            # If not, create new vertical lines in the plot.\n",
    "            linesList.append( ax.axvline(x=imageObj.frame, color=color, linestyle=\"--\") )\n",
    "        lineIndx += 1\n",
    "    \n",
    "    while lineIndx < len(linesList):\n",
    "        linesList[lineIndx].set_xdata(allImages[0].frame)\n",
    "        linesList[lineIndx].set_color(\"r\")\n",
    "        lineIndx += 1\n",
    "    \n",
    "    maxDeltaCV,maxDeltaFrames = find_max_delta(allImages)\n",
    "    \n",
    "    index = 0\n",
    "    for delta in maxDeltaCV:\n",
    "        deltaLines[index].set_ydata(delta)\n",
    "        index += 1\n",
    "    \n",
    "    ax1.set_ylim([0,max(maxDeltaCV)*1.1])\n",
    "    \n",
    "    fig.canvas.draw()\n",
    "\n",
    "def calcIntermediates():\n",
    "    \n",
    "    global allImages,numImgs\n",
    "    \n",
    "    # Gets current Maximum Delta\n",
    "    maxD = float(maxDeltaFT.value)\n",
    "    \n",
    "    # Clears all intermediate images, and added last frame.\n",
    "    allImages = [img for img in allImages if img.ftype in [0,1]]\n",
    "    numImgs = len(allImages)\n",
    "    \n",
    "    # Add last frame as main image, if requested.\n",
    "    if addLastTB.value:\n",
    "        if (allImages[-1].frame != (len(cvSMD)-1)):\n",
    "            allImages.append( imageDat(frame=len(cvSMD)-1, ftype=5) )\n",
    "            numImgs = len(allImages)\n",
    "    \n",
    "    interImgs = []\n",
    "    \n",
    "    # Loops over all Equally Spaced \"gaps\"\n",
    "    for indx in range(1,numImgs):\n",
    "        \n",
    "        frameFin = allImages[indx].frame\n",
    "        frameIni = allImages[indx -1].frame\n",
    "        \n",
    "        tmpIni = frameIni\n",
    "        tmpFin = frameFin\n",
    "        \n",
    "        cvCheck = np.asarray([False]*numColvars)\n",
    "        \n",
    "        while tmpIni != frameFin:\n",
    "            \n",
    "            # Looks for the first cv index that needs to be checked.\n",
    "            if len(np.where(cvCheck == 0)[0]):\n",
    "                colIndx = np.where(cvCheck == 0)[0][0]\n",
    "            else:\n",
    "                # If all indices are true, update the tmpIni and reset the count.\n",
    "                cvCheck = np.asarray([False]*numColvars)\n",
    "                # In this order!\n",
    "                tmpIni = tmpFin\n",
    "                tmpFin = frameFin\n",
    "                \n",
    "                # If we found a new image, add it to the list\n",
    "                if tmpIni > frameIni and tmpIni != frameFin:\n",
    "                    interImgs.append(tmpIni)\n",
    "\n",
    "            # If the difference in this colvar is less than the maximum,\n",
    "            # skip to the next.\n",
    "            tmpdelta = abs( cvSMD.iloc[tmpFin,colIndx] - cvSMD.iloc[tmpIni,colIndx] )\n",
    "            if tmpdelta <= maxD:\n",
    "                cvCheck[colIndx] = True\n",
    "                continue\n",
    "            \n",
    "            # Ugly trick to get the first frame, between two images, with a ColVar\n",
    "            # variation larger than maxD. \n",
    "            # R + data.table would do this much more \"cleanly\", but I'll stop complaining...\n",
    "            baseVal = cvSMD.iloc[tmpIni,colIndx]\n",
    "            queryStr = \"abs(\" + str(baseVal) + \"-\" + str(cvSMD.columns[colIndx]) + \") >= \" + str(maxD)\n",
    "            indxV = cvSMD.iloc[tmpIni:tmpFin,:].query(queryStr)\n",
    "            \n",
    "            # In case there is at least one frame where this ColVar has a variation greater than maxD,\n",
    "            # we store the frame.\n",
    "            if len(indxV) > 0 :\n",
    "                # We also make sure that the image is *after* the initial frame.\n",
    "                tmpNewFrame = (indxV.index[0] - 1)\n",
    "                if tmpNewFrame > tmpIni:\n",
    "                    tmpFin = tmpNewFrame\n",
    "                    # If we found a new intermeiate image, uncheck all colvars and check the current one.\n",
    "                    cvCheck = np.asarray([False]*numColvars)\n",
    "                    cvCheck[colIndx] = True\n",
    "                    \n",
    "            else:\n",
    "                cvCheck[colIndx] = True\n",
    "    \n",
    "    # Add new images to the global list\n",
    "    for frame in interImgs:\n",
    "        allImages.append( imageDat(frame=frame, ftype=2) )\n",
    "    \n",
    "    # Sort all images\n",
    "    allImages = sorted(allImages, key=lambda x: x.frame)\n",
    "    numImgs = len(allImages)\n",
    "\n",
    "def toggle_calcInter(*args):\n",
    "    calcIntermediates()\n",
    "    update_plot()\n",
    "\n",
    "calcInterB.on_click(toggle_calcInter)\n",
    "calcInterB.on_click(update_TotImgs)\n",
    "\n",
    "def update_offset(*args):\n",
    "    '''Call function to reset all initial images using the new value for offset.'''\n",
    "    resetImages(imagesText.value)\n",
    "    update_plot()\n",
    "\n",
    "offsetText.observe(update_offset,'value')\n",
    "\n",
    "# Connects inputs\n",
    "def toggle_intermediate(*args):\n",
    "    # Activate user input\n",
    "    \n",
    "    # Change icon and Lock Image Selection.\n",
    "    if activateInterTB.value:\n",
    "        activateInterTB.icon = \"check-square\"\n",
    "        maxDeltaFT.disabled = False\n",
    "        addLastTB.disabled = False\n",
    "        calcInterB.disabled = False\n",
    "        \n",
    "        imagesText.disabled = True\n",
    "        imagesSlider.disabled = True\n",
    "        offsetText.disabled = True\n",
    "    else:\n",
    "        activateInterTB.icon = \"square\"\n",
    "        maxDeltaFT.disabled = True\n",
    "        addLastTB.disabled = True\n",
    "        calcInterB.disabled = True\n",
    "        \n",
    "activateInterTB.observe(toggle_intermediate, \"value\")\n",
    "\n",
    "def update_imgSel(*args):\n",
    "    \n",
    "    global allImages\n",
    "    \n",
    "    optTouplList = [ (str(indx) + \" (\" + str(img.frame) + \")\",indx) for indx,img in enumerate(allImages) ]\n",
    "    \n",
    "    imgSelect.options = optTouplList\n",
    "    \n",
    "    imgSelect.value = optTouplList[0][1]\n",
    "\n",
    "totalImagesText.observe(update_imgSel,'value')\n",
    "\n",
    "def selection_update(*args):\n",
    "    \n",
    "    global allImages\n",
    "    \n",
    "    imgChgFrameText.value = allImages[ imgSelect.value ].frame\n",
    "    \n",
    "    for img in allImages:\n",
    "        img.selected = False\n",
    "    allImages[ imgSelect.value ].selected = True\n",
    "    \n",
    "    update_plot()\n",
    "\n",
    "imgSelect.observe(selection_update,'value')\n",
    "\n",
    "def update_image(*args):\n",
    "    \n",
    "    global allImages\n",
    "    \n",
    "    newVal = imgChgFrameText.value\n",
    "    \n",
    "    if newVal >= numFrames -1 or newVal < 0:\n",
    "        return\n",
    "    \n",
    "    switchOccurred = False\n",
    "    \n",
    "    if newVal < allImages[ imgSelect.value ].frame:\n",
    "        # If there is a next image, we make sure we don't change the order of images\n",
    "        if imgSelect.value - 1 >= 0:\n",
    "            if newVal < allImages[ imgSelect.value - 1 ].frame:\n",
    "                switchOccurred = True\n",
    "        \n",
    "    if newVal > allImages[ imgSelect.value ].frame:\n",
    "        # If there is a next image, we make sure we don't change the order of images\n",
    "        if imgSelect.value + 1 <= numImgs - 1:\n",
    "            if newVal > allImages[ imgSelect.value + 1 ].frame:\n",
    "                switchOccurred = True\n",
    "    \n",
    "    allImages[ imgSelect.value ].frame = newVal\n",
    "    if switchOccurred:\n",
    "        allImages = sorted(allImages, key=lambda x: x.frame)\n",
    "        update_imgSel()\n",
    "        \n",
    "        for indx,img in enumerate(allImages):\n",
    "            if img.frame == newVal:\n",
    "                imgSelect.value = indx\n",
    "        \n",
    "        selection_update()\n",
    "        \n",
    "    update_plot()\n",
    "\n",
    "imgChgFrameText.observe(update_image, 'value')\n",
    "\n",
    "def add_img(*args):\n",
    "    \n",
    "    global allImages,numImgs\n",
    "    \n",
    "    selectedImg = imgSelect.value\n",
    "    \n",
    "    iniFrame = allImages[ selectedImg ].frame\n",
    "    if selectedImg < numImgs-1:\n",
    "        finFrame = allImages[ selectedImg + 1 ].frame\n",
    "    else:\n",
    "        if iniFrame < numFrames-1:\n",
    "            finFrame = numFrames-1\n",
    "        else:\n",
    "            return\n",
    "    \n",
    "    newframe = int(np.floor((iniFrame + finFrame)/2))\n",
    "    \n",
    "    if (newframe != iniFrame):\n",
    "        allImages.append( imageDat(frame=newframe, ftype=1) )\n",
    "        \n",
    "        # Sort all images\n",
    "        allImages = sorted(allImages, key=lambda x: x.frame)\n",
    "        numImgs = len(allImages)\n",
    "        \n",
    "        update_plot()\n",
    "        update_TotImgs()\n",
    "        imgSelect.value = selectedImg + 1\n",
    "\n",
    "addImg.on_click(add_img)\n",
    "\n",
    "def rm_img(*args):\n",
    "    \n",
    "    global allImages,numImgs\n",
    "    \n",
    "    selectedImg = imgSelect.value\n",
    "    \n",
    "    # If tehre are only three, don't do anyhting or we would get less images than the minimum.\n",
    "    if numImgs == 3:\n",
    "        return\n",
    "    \n",
    "    # After deleting the image, select the one before it.\n",
    "    newSel = 0\n",
    "    if selectedImg == 0:\n",
    "        newSel = 1\n",
    "    else:\n",
    "        newSel = selectedImg -1\n",
    "    \n",
    "    # deletes the image\n",
    "    del allImages[selectedImg]\n",
    "    \n",
    "    numImgs = len(allImages)\n",
    "    \n",
    "    # Updates widget and plot\n",
    "    update_plot()\n",
    "    update_TotImgs()\n",
    "    imgSelect.value = newSel\n",
    "    \n",
    "    \n",
    "rmImg.on_click(rm_img)\n",
    "    \n",
    "update_plot()\n",
    "\n",
    "update_TotImgs()\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "###  2 TABS = all frames (Value vs Frame) + MaxVar | All Images (Value vs Image Index) + Min/Avr/Max Var\n",
    "### Color backgroud by Max Var ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Steps:\n",
    "- Set the Initial Images value to 5.\n",
    "- Set an Offset of 150 frames.\n",
    "- Check the ``Add intermediates`` box and set the maximum delta to *0.6*.\n",
    "- Click *calculate* (this button should be used every time you change the delta).\n",
    "    - The *Total Images* should indicate 14 images were created. \n",
    "\n",
    "Now we have selected 14 images, and kept the maximum colvar distance to reasonable levels. You can fine tune this initial set of frames using the drop-down box, selecting an individual image and either delete it or change its frame. New images can also be added using this interface. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-02T18:38:44.458974Z",
     "start_time": "2018-03-02T18:38:44.450100Z"
    }
   },
   "outputs": [],
   "source": [
    "frames = [f.frame for f in allImages]\n",
    "\n",
    "print(\"We selected\",numImgs,\"images:\",frames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create _position_ and _periodic box_ files to initiate String calculation\n",
    "\n",
    "Now that we have selected the frames from our SMD simulations that will compose the statirng points of the string optimization, we must extract and create the necessary files for the string walkers. Both the COOR (a binary format used by NAMD to set atoms' positions) and the XSC (a text file that indicates to NAMD the size of the periodic box) need to be created for each image.\n",
    "\n",
    "The following cell has python code that creates the necessary directories for setup and output, and runs the *CATDCD*  software (distributed along with VMD) that specializes in reading and writing MD trajectory files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-02T17:07:16.390992Z",
     "start_time": "2018-03-02T17:07:15.692547Z"
    }
   },
   "outputs": [],
   "source": [
    "# Erases setup files if they exist. We start from scratch to avoid future errors.\n",
    "if os.path.exists(basePath):\n",
    "    shutil.rmtree(basePath)\n",
    "\n",
    "# Creates setup directory (where coor and xsc files will be written), and base output directory\n",
    "os.makedirs(setupPath)\n",
    "os.makedirs(outputPath)\n",
    "\n",
    "# Creates output directories for each image\n",
    "if parallelExec:\n",
    "    numOutFolders = numImgs*stringNumWalkers\n",
    "elif sequentialExec:\n",
    "    numOutFolders = numImgs\n",
    "\n",
    "for i in range(numOutFolders):\n",
    "    os.makedirs(outputPath + \"/\" + str(i))\n",
    "\n",
    "print(\"Total selected frames:\",len(frames))\n",
    "\n",
    "# Main loop that creates files for each image.\n",
    "stringFileList = []\n",
    "stringIndex = 0\n",
    "for frame in frames:\n",
    "    \n",
    "    # Creates a list of DCD files that will be used to write configuration files for the String.\n",
    "    stringFileList.append(\"image_\"+str(stringIndex)+\".coor\")\n",
    "    \n",
    "    # Writes selected frames to DCD files for String input\n",
    "    print(\"Writing frame \" + str(frame) + \" to file \" + stringFileList[-1])\n",
    "    command = catdcdPath + \" -o \" + setupPath + stringFileList[-1] + \" -otype namdbin -stype psf -s \" + psfFile + \\\n",
    "        \" -first \" + str(frame+1) + \" -last \" + str(frame+1) + \" -dcd \" + dcdFile\n",
    "    process = subprocess.Popen(command, shell=True)\n",
    "    process.wait()\n",
    "\n",
    "    # Exact frames from the XST file for String input\n",
    "    # Each XSC file is a text file with the same header, which we create below.\n",
    "    sedFile = setupPath + \"image_\" + str(stringIndex) + \".xsc\"\n",
    "    tmpFile = open(sedFile,\"w\")\n",
    "    tmpFile.write(\"\"\"# NAMD extended system trajectory file\n",
    "#$LABELS step a_x a_y a_z b_x b_y b_z c_x c_y c_z o_x o_y o_z\n",
    "\"\"\")\n",
    "    tmpFile.close()\n",
    "    # Now we get the data from the chosen frame and append it to the file for the current image.\n",
    "    sedLine = \"cat \" + xstFile + \" | grep -v ^# | sed -n \" + str(frame + 1) + \"p >> \" + sedFile \n",
    "    process = subprocess.Popen(sedLine, shell=True)\n",
    "    process.wait()\n",
    "    \n",
    "    stringIndex += 1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create _configuration_ files to initiate String calculation\n",
    "\n",
    "We now create the configuration files needed to initiate all walkers from the string optimization, one file per string image. All files have the same basic syntax, but differ in individual values that describethe collective variables starting points, since each image represents a different point along the chemical reaction being studyed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create the common block of text for defining colective variables, which will be used in all configuration files, and use python's \"format\" capability to dynamically insert data for each colvar, in each image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# In Python's \"format\" module, every keyword enclosed in curly brackets will be substituted by a value \n",
    "#  provided in a dictionary. If the final string needs actual curly brackets, they can be excaped by using \n",
    "#  two consecutive brackets.\n",
    "\n",
    "# Using example values, the block should read:\n",
    "#    # Human Readable Label - just a comment\n",
    "#    colvar {\n",
    "#        name colvar_name\n",
    "#        distance {\n",
    "#            group1 {\n",
    "#                atomNumbers 123\n",
    "#            }\n",
    "#            group2 {\n",
    "#                atomNumbers 456\n",
    "#            }\n",
    "#        }\n",
    "#    }\n",
    "#\n",
    "#    harmonic {\n",
    "#        name bias_name\n",
    "#        colvars colvar_name\n",
    "#        centers 1.234\n",
    "#        forceConstant 0.0\n",
    "#    }\n",
    "\n",
    "# Therefore, using keywords for values that will be substituted later, our string variable will read:\n",
    "\n",
    "colvarBlock = \"\"\"\n",
    "# {CVlabel}\n",
    "colvar {{\n",
    "    name swarms_atmDist{CVcolIndx}\n",
    "    distance {{\n",
    "        group1 {{\n",
    "            atomNumbers {atmIndx1}\n",
    "        }}\n",
    "        group2 {{\n",
    "            atomNumbers {atmIndx2}\n",
    "        }}\n",
    "    }}\n",
    "}}\n",
    "\n",
    "harmonic {{\n",
    "    name bias_swarms_atmDist{CVcolIndx}\n",
    "    colvars swarms_atmDist{CVcolIndx}\n",
    "    centers {CVcenter}\n",
    "    forceConstant 0.0\n",
    "}}\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, we run the loop that uses the string we created in the precious cell to write each configuration file, using the colvar values from the SMD frames previously selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loop over all images.\n",
    "for imgIndx in range(len(frames)):\n",
    "    \n",
    "    # We create one file per image in the string.\n",
    "    stringFile = open(setupPath + \"colvars_\" + str(imgIndx) + \".conf\",'w')\n",
    "    \n",
    "    # Common colvar communication freqeuncy\n",
    "    stringFile.write(\"\\ncolvarsTrajFrequency 5\\n\\n\")\n",
    "    \n",
    "    # Loop over all colvars\n",
    "    for colIndx in range(len(cvSMD.columns)):\n",
    "        \n",
    "        # We aggregate the information for the currently selected colvar and image in a python dictionary.\n",
    "        currInfo = {}\n",
    "        currInfo['CVlabel'] = atomDistList[colIndx][2]\n",
    "        currInfo['CVcolIndx'] = str(colIndx)\n",
    "        # Convert atom indices from VMD 0-base index to ColVar 1-base index.  \n",
    "        currInfo['atmIndx1'] = atomDistList[colIndx][0] + 1\n",
    "        currInfo['atmIndx2'] = atomDistList[colIndx][1] + 1\n",
    "        # Get colvar value from the selected SMD frame\n",
    "        currInfo['CVcenter'] = cvSMD.iloc[frames[imgIndx],colIndx]\n",
    "        \n",
    "        # Finally, we substitute the values in the common string block and write it to the file.\n",
    "        # We write one block of text per colvar.\n",
    "        stringFile.write( colvarBlock.format(**currInfo) )\n",
    "    \n",
    "    \n",
    "    stringFile.close()\n",
    "######## end-of-loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create _configuration_ file to initiate eABF free energy calculations\n",
    "\n",
    "We can also create the configuration file needed to initiate all walkers for the eABF free energy calculation. We have not ran the string simulations yet, but we can decide on file names right now :).\n",
    "For eABF, a single file is needed since no unique (walker-specific) values will be written into them. The only information needed are the atom indices for atoms involved in collective variables, and file names that will hold the end results of the string optimization. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Path Collective Variables use the RMSD of relevant atoms to describe the reaction coordinate. For that, we need \n",
    "# the atom numbers of all atoms involved in colvars we set at the begining of the notebook.\n",
    "# We can use python's flexibility to create a unique set of atom indices in one line.\n",
    "# Warning: Computational tricks ahead! Using \"list comprehension\", we can create a list by calling each colvar\n",
    "#    (or cv) in atomDistList, select only the frist two items (the 0-based atom indices) and create\n",
    "#    a set with them shifted to 1-based index. List comprehensions, automatically combine all pairs in a list \n",
    "#    (by doing this operation within brackets), and the resulting list of sets is unpacked by placing an '*' in \n",
    "#    fornt of it. The unpacked list of sets is used by the 'set.union' operation to create a union of all sets. \n",
    "\n",
    "atmsUnion = set.union(*[ set([cv[0]+1,cv[1]+1]) for cv in atomDistList ])\n",
    "\n",
    "# An equivalent operation that would be easier to read can be found below, but we can take the oportunity to \n",
    "# get used to pyhton's toolset, writing 1 line of code instead of 4.\n",
    "## atmsUnion = set()\n",
    "## for cv in atomDistList:\n",
    "##     for atm in [cv[0]+1,cv[1]+1]:\n",
    "##         atmsUnion.add(atm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know the relevant atoms whose positions we must track, we follow a similar approach as the one used for the string configuration files and use easily edditable text blocks to create the final eABF configuration file. This configuration file creates two collective variables that are *scripted*, meaning they are not coded in the colvar module within NAMD, but are calculated from functions provided in a TCL script called **pathCV.tcl**. You can find more inforamtion in the tutorial PDF file that accompanies this notebook. \n",
    "Both collective variables use the same set of PDB files that result from the string optimization.\n",
    "\n",
    "We will configure the eABF simulation to sample only along the 's' variable. This will not give us the best possible estimate, but will be computationally cheaper since we need only sample one dimension. The eABF algorithm will apply forces on the virtual particle in a region only after acquiring 100 samples in that region. We are using multiple walkers initialized along the string to speed up the sampling process. The 'z' colective variable is only used with a harmonic restraint to keep the eABF virtual particle from wandering away from the main path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Sets common text block that will be used to set all points along the optimized string. \n",
    "# This defines the path that eABF will use to calculate the free energy.\n",
    "# The list of atoms is common to all points that define the path, so we add it to the common text block now.\n",
    "allAtmsList = [atom for atom in atmsUnion]\n",
    "# We need to sort the atom indices! ColVars modules needs ordered indices for RMSD calculation.\n",
    "allAtmsList.sort()\n",
    "# Now we create a list of strings, so we can use the \"join\" method to \n",
    "# join all atom numbers with spaces in between, passing a list of strings to the join method.\n",
    "allAtmsStr = ' '.join([str(atom) for atom in allAtmsList])\n",
    "pathNodeStr = \"\"\"\n",
    "rmsd {{\n",
    " atoms {{\n",
    "   atomNumbers {{ \"\"\" + allAtmsStr + \"\"\" }}\n",
    " }}\n",
    " refpositionsfile ./string/string_{imgIndx}.pdb\n",
    "}}\n",
    "\"\"\"\n",
    "\n",
    "# Only one configuration file is needed for eABF.\n",
    "pathCVFile = open(setupPath + \"pathCV.conf\",'w')\n",
    "\n",
    "# Common colvar communication freqeuncy\n",
    "pathCVFile.write(\"colvarsTrajFrequency 5\\n\\n\")\n",
    "\n",
    "# Now we set the two collective variables needed for eABF, the path collective variables.\n",
    "\n",
    "# Colvar 's' defines the position of the simulations along the path.\n",
    "begCVs = \"\"\"\n",
    "colvar {\n",
    "   name                  s\n",
    "\n",
    "   width                 0.05\n",
    "   lowerboundary         0.0\n",
    "   upperboundary         1.0\n",
    "   lowerwallconstant     10.0\n",
    "   upperwallconstant     10.0\n",
    "\n",
    "   scriptedFunction      pathCVs\n",
    "   extendedLagrangian    on\n",
    "   extendedFluctuation   0.01\n",
    "\n",
    "\"\"\"\n",
    "pathCVFile.write(begCVs)\n",
    "\n",
    "# Writes all the definitions of points along the path.\n",
    "for imgIndx in range(len(frames)):\n",
    "    pathCVFile.write( pathNodeStr.format(**{'imgIndx': imgIndx}) )\n",
    "\n",
    "# Closes the first colvar block\n",
    "pathCVFile.write(\"}\\n\\n\")\n",
    "\n",
    "# Colvar 'z' defines the distance of the simulations 'perpendicular' to the path.\n",
    "begCVz = \"\"\"\n",
    "colvar {\n",
    "   name                  z\n",
    "\n",
    "   width                 0.05\n",
    "   lowerboundary        -0.5\n",
    "   upperboundary         1.5\n",
    "\n",
    "   scriptedFunction      pathCVz\n",
    "   extendedLagrangian    on\n",
    "   extendedFluctuation   0.01\n",
    "\"\"\"\n",
    "pathCVFile.write(begCVz)\n",
    "\n",
    "# Writes all the definitions of points along the path.\n",
    "for imgIndx in range(len(frames)):\n",
    "    pathCVFile.write( pathNodeStr.format(**{'imgIndx': imgIndx}) )\n",
    "\n",
    "# Closes the second colvar block\n",
    "pathCVFile.write(\"}\\n\\n\")\n",
    "\n",
    "# Writes the configuration for the ABF procedure. We call this \"extended\"-ABF because the colvars \n",
    "# used in the sampling have been configured to use the extended-lagrangian option.\n",
    "configPathCV = \"\"\"\n",
    "abf {\n",
    "    #colvars s z\n",
    "    colvars s\n",
    "    fullSamples 100\n",
    "\n",
    "    # Multiple-replica ABF\n",
    "    shared yes\n",
    "    sharedFreq 5\n",
    "    \n",
    "    # CZAR is the default\n",
    "    CZARestimator yes\n",
    "    writeCZARwindowFile yes\n",
    "    \n",
    "}\n",
    "\n",
    "harmonic {\n",
    "   colvars               z\n",
    "   centers               0.0\n",
    "   forceConstant         0.5\n",
    "}\n",
    "\n",
    "\"\"\"\n",
    "pathCVFile.write(configPathCV)\n",
    "\n",
    "pathCVFile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Populate the String/eABF work directory with the necessary files.\n",
    "\n",
    "We now copy essential files from the decarboxylase system (such as PDB, PSF, topology) and configuration files from the 'templates' directory.\n",
    "Some of files need specific data that is only knwon after we setup our system, decide on parameters, chose images, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Get structure data and NAMD parameters\n",
    "shutil.copy(pdbFile, basePath)\n",
    "shutil.copy(psfFile, basePath)\n",
    "\n",
    "try:\n",
    "    shutil.copytree('./setup/toppar', basePath + '/toppar', )\n",
    "except:\n",
    "    pass\n",
    "\n",
    "# Get TCL code and parameters for String and eABF calculation\n",
    "for filename in glob.glob('./templates/*'):\n",
    "    shutil.copy(filename, basePath)\n",
    "\n",
    "# Create directory in RAM for faster QM calculation and I/O.\n",
    "# This path is used in the \"namd.conf\" file.\n",
    "if sys.platform == 'linux':\n",
    "    os.makedirs(\"/dev/shm/NAMD_SM_qmmm\", exist_ok=True)\n",
    "\n",
    "with fileinput.FileInput(basePath + \"/namd.conf\", inplace=True, backup='.bak') as file:\n",
    "    for line in file:\n",
    "        print(line.replace(\"QMBINPATH\", qmBinPath), end='')\n",
    "\n",
    "with fileinput.FileInput(basePath + \"/string_param.conf\", inplace=True, backup='.bak') as file:\n",
    "    for line in file:\n",
    "        print(line.replace(\"STRINGITERATIONS\", str(stringNumIter)), end='')\n",
    "\n",
    "with fileinput.FileInput(basePath + \"/string_param.conf\", inplace=True, backup='.bak') as file:\n",
    "    for line in file:\n",
    "        print(line.replace(\"NUMWALKERS\", str(stringNumWalkers)), end='')\n",
    "\n",
    "if parallelExec:\n",
    "    eABFWalkersPerImg = stringNumWalkers\n",
    "elif sequentialExec:\n",
    "    eABFWalkersPerImg = 1\n",
    "    \n",
    "with fileinput.FileInput(basePath + \"/eABF.conf\", inplace=True, backup='.bak') as file:\n",
    "    for line in file:\n",
    "        print(line.replace(\"EABFSTEPS\", str(eABFnumSteps)), end='')\n",
    "\n",
    "with fileinput.FileInput(basePath + \"/eABF.conf\", inplace=True, backup='.bak') as file:\n",
    "    for line in file:\n",
    "        print(line.replace(\"EABFWPERIMG\", str(eABFWalkersPerImg)), end='')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize any file in the notebook, use the magic (USE quotes arounf the file name):\n",
    "\n",
    "**%load \"file-name\"**\n",
    "\n",
    "To save modifications made to file, *add* the following magic to the top of the cell (DO NOT use quotes around the file name!):\n",
    "\n",
    "**%%writefile file-name**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Run the calculation!\n",
    "\n",
    "For fully parallel calculations in clusters (using multiple nodes), where all walkers from all swarms of all images walk at the same time, use the file *string_parallel.tcl*.\n",
    "For single-node calculations, where wlakers of a same image walk sequentially, but walkers of different images walk simultaneously, use *string_sequential.tcl*. Warning: The node needs to have at least as many processor cores as there are images in the string.\n",
    "\n",
    "An example of a command that can be ran in the terminal window would be (if we had 8 images):\n",
    "\n",
    "__mpirun -n 8 namd2 +replicas 8 string_sequential.tcl +stdout output/%1d/job00.%1d.log__\n",
    "\n",
    "The following cell is a *bash* cell that can be used to run NAMD. Remember that this process can take **several hours**, and for the purposes of the QM/MM Workshop, we have provided pre-ran files in the *example-output* folder. You can simply copy the pre-ran files to the *String_eABF/output* folder.\n",
    "\n",
    "If you still wish to proceed, just remove the comment (the *hash* character \"#\") from the NAMD line bellow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$numImgs\"  \"$namdPathMPI\"\n",
    "\n",
    "echo \"Calculating $1 images...\"\n",
    "\n",
    "cd String_eABF/\n",
    "\n",
    "echo \"mpirun -n $1 $2 +replicas $1 string_sequential.tcl +stdout output/%1d/job00.%1d.log\"\n",
    "#mpirun -n $1 $2 +replicas $1 string_sequential.tcl +stdout output/%1d/job00.%1d.log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# String Method analysis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running the string optimization, we will parse each output file in order to get the evolution of the restraints applied by the string method on the atom distances. As the string iterations progress, we will notice that the harmonic restraints are changing very little, which indicates that the method has converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# We are looking for the following pattern inside the log file:\n",
    "# Updating bias_swarms_atmDist1\n",
    "# colvars:   # centers = { 4.14177 }\n",
    "    \n",
    "# First, we compile a Regular Expression to look for the correct format.\n",
    "# Regular expression are a powerful tool to match patterns in text files that have complex formats.\n",
    "regexVal = re.compile(\"^colvars:   # centers = { ([0-9\\.]*) }$\")\n",
    "regexName = re.compile(\"^Updating bias_swarms_atmDist([0-9]*)$\")\n",
    "\n",
    "# Create a Pandas dataframe with systematic column names (so it's easy to access in the future).\n",
    "cvString = pd.DataFrame(columns=[ \"dist\"+str(i) for i in range(numColvars) ] + [\"Image\",\"Iter\"])\n",
    "\n",
    "# Initialize with colvar values taken from SMD simulation that initiated the String optimization.\n",
    "for i in range(numImgs):\n",
    "    cvString.loc[i] = list(cvSMD.iloc[frames[i],:]) + [i, 0]\n",
    "\n",
    "# Loop over all images\n",
    "for image in range(numImgs):\n",
    "    \n",
    "    if parallelExec:\n",
    "        foldIndx = image*stringNumWalkers\n",
    "    else:\n",
    "        foldIndx = image\n",
    "    \n",
    "    logFiles = []\n",
    "    folderStr = outputPath + str(foldIndx)\n",
    "    for logFN in glob.iglob(folderStr + \"/job*.log\"):\n",
    "        logFiles.append(logFN)\n",
    "    \n",
    "    ## Assuming there is only one log file per image, for now.\n",
    "    fileName = logFiles[0]\n",
    "    \n",
    "    # Set temporary variables\n",
    "    colvarName = \"\"\n",
    "    numIterations = 1\n",
    "    tmpDF = pd.DataFrame()\n",
    "    \n",
    "    with open(fileName,'r') as fileHandle:\n",
    "        print(\"Processing file\",fileName)\n",
    "        \n",
    "        for line in fileHandle:\n",
    "            matchVal = regexVal.match(line.strip())\n",
    "            matchName = regexName.match(line.strip())\n",
    "            \n",
    "            # First we check if the line provides a name for a colvar (the value will follow in the next line)\n",
    "            if matchName:\n",
    "                colvarName = line.strip().split(\" \")[1]\n",
    "                cvIndx = matchName.group(1)\n",
    "            \n",
    "            # We only consider matching value strings after we match at least one name string.\n",
    "            if matchVal and colvarName:\n",
    "                colvarValue = matchVal.group(1)\n",
    "                \n",
    "                tmpDF.loc[0,\"dist\" + str(cvIndx)] = round(float(colvarValue),4)\n",
    "                \n",
    "            # After accumulating values for all colvars, append to data frame.\n",
    "            if len(tmpDF.columns) == numColvars:\n",
    "                tmpDF.loc[0,\"Image\"] = image\n",
    "                tmpDF.loc[0,\"Iter\"] = numIterations\n",
    "                \n",
    "                cvString = cvString.append(tmpDF)\n",
    "                \n",
    "                numIterations += 1\n",
    "                # Resets temporary data frame\n",
    "                tmpDF = pd.DataFrame()\n",
    "\n",
    "cvString.reset_index(drop=True, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the raw data has been read from all output files, we process the data to determine what is the variation of the string harmonic restraint value along the optimization iterations. The value of the restraint for each colvar is compared to the value at the last iteration, so the deviation will approach zero as the optimization progresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sets the \"Image\" column in the data frame as an integer instead of float (for better visualization)\n",
    "cvString[\"Image\"] = cvString[\"Image\"].astype('int', copy=False)\n",
    "# First, we create a \"melt\" version of the dataset for easyer analysis and plotting.\n",
    "cvStringMelt = pd.melt(cvString, id_vars=[\"Image\",\"Iter\"], var_name=\"colvar\", value_name=\"Distance\")\n",
    "\n",
    "# We create new column names\n",
    "cols = [ \"dist\"+str(i) for i in range(numColvars) ]\n",
    "colsTrg = [ \"dist\"+str(i)+\"Dev\" for i in range(numColvars) ]\n",
    "\n",
    "# Now we define a function that will calculate the deviation of a restraint when compared to the \n",
    "# restraint value at the last string iteration. This used groups created by the Pandas DataFrame below. \n",
    "def calculate_deviation(group):\n",
    "    # We use the new column names to hold the deviation values for each colvar.\n",
    "    group[colsTrg] = abs(group.loc[:, cols] - group.loc[group.Iter==max(group.Iter), cols].values[0])\n",
    "    return group\n",
    "\n",
    "# For each Image, apply the function above to calculate the progression of the restraints over all colvars.\n",
    "cvStat = cvString.groupby(['Image']).apply(calculate_deviation)\n",
    "\n",
    "# Melts the dataframe for easyer analysis and plotting, saving only the deviations from the final colvar value.\n",
    "cvStatMelt = pd.melt(cvStat.loc[:,[\"Image\",\"Iter\"] + colsTrg], id_vars=[\"Image\",\"Iter\"], \\\n",
    "                 var_name=\"colvar\", value_name=\"Deviation\")\n",
    "\n",
    "# Calculates the mean deviation per iteration, for each colvar. The mean is calculated accross all walkers of an image.\n",
    "cvStat = cvStatMelt.loc[:,cvStatMelt.columns != \"Image\"].groupby([\"Iter\",\"colvar\"], as_index=False).mean()\n",
    "# Similarly, determines the maximum deviation accross all walker of an image.\n",
    "cvStat['MaxDev'] = cvStatMelt.loc[:,cvStatMelt.columns != \"Image\"].groupby([\"Iter\",\"colvar\"], as_index=False).max()['Deviation']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "pd.set_option('mode.chained_assignment',None)\n",
    "\n",
    "radioStat = widgets.RadioButtons(\n",
    "    options=['Mean Deviation', 'Maximum Deviation'],\n",
    "    value='Mean Deviation',\n",
    "    #description='Pizza topping:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "stringCtl = widgets.HBox([widgets.Label(value=\"Statistic:\"), radioStat])\n",
    "\n",
    "fig = None\n",
    "axs = None\n",
    "\n",
    "def plotString(*args):\n",
    "    \n",
    "    global fig, axs\n",
    "    \n",
    "    if radioStat.value == \"Mean Deviation\":\n",
    "        selStr = 'Deviation'\n",
    "    else:\n",
    "        selStr = 'MaxDev'\n",
    "    \n",
    "    p = p9.ggplot(cvStat, p9.aes(x=\"Iter\", y=selStr, color=\"colvar\")) + \\\n",
    "        p9.geom_point() + p9.geom_path() + p9.theme_linedraw() + \\\n",
    "        p9.labs(x=\"String Iteration\", y=\"Deviation (A)\") +\\\n",
    "        p9.scale_colour_manual(values = [\"blue\", \"green\", \"red\", \"cyan\", \"magenta\"], \n",
    "                               labels=[ keyNameDic['dist'+str(i)] for i in range(len(keyNameDic)) ])\n",
    "        \n",
    "    if fig is None:\n",
    "        fig, plot = p.draw(return_ggplot=True)\n",
    "        axs = plot.axs\n",
    "    else:\n",
    "        for artist in plt.gca().lines +\\\n",
    "                        plt.gca().collections +\\\n",
    "                        plt.gca().artists + plt.gca().patches + plt.gca().texts:\n",
    "            artist.remove()\n",
    "        p._draw_using_figure(fig, axs)\n",
    "    \n",
    "radioStat.observe(plotString, 'value')\n",
    "\n",
    "display(stringCtl)\n",
    "\n",
    "radioStat.value = 'Maximum Deviation'\n",
    "plt.tight_layout()\n",
    "radioStat.value = 'Mean Deviation'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will notice that the **Mean** deviation is less than 0.2 Angstroms from the 3rd iteration onwards, for all colvars, and the **Maximum** deviation is less than 0.25 A from the 7th iteration onwards.\n",
    "For a production run, one could require the mean deviation to be less than 0.2 for the last 10 iteration rounds, and the maximum deviation to be less than 0.2 for the last 5, for example, but for the purposes of this tutorial, we will consider this string to be converged \"enough\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following interactive plot will show you the evolution of the constraint values themselves during the string optimization iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "pd.set_option('mode.chained_assignment',None)\n",
    "\n",
    "playString  = widgets.Play(\n",
    "    interval=1500,\n",
    "    value=1,\n",
    "    min=0,\n",
    "    max=numIterations-1,\n",
    "    step=1,\n",
    "    description=\"Iterate over the string optimization\",\n",
    ")\n",
    "\n",
    "playStringSlider = widgets.IntSlider(value=1, min=0, max=numIterations-1)\n",
    "widgets.jslink((playString,'value'),(playStringSlider,'value'))\n",
    "\n",
    "showInitCB = widgets.Checkbox(\n",
    "    value=False,\n",
    "    description='Show initial',\n",
    "    disabled=False\n",
    ")\n",
    "showInitTB = widgets.ToggleButton(\n",
    "    value=False,\n",
    "    description='Show Initial Values',\n",
    "    disabled=False,\n",
    "    button_style='info', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Shows the initial values of colvars for all images.',\n",
    ")\n",
    "\n",
    "stringCtl = widgets.HBox([widgets.Label(value=\"Iteration:\"), playString, playStringSlider, showInitTB])\n",
    "\n",
    "maxY = max(cvStringMelt['Distance'])\n",
    "\n",
    "fig = None\n",
    "axs = None\n",
    "\n",
    "def plotString(*args):\n",
    "    \n",
    "    global fig, axs\n",
    "    \n",
    "    newIter = playStringSlider.value\n",
    "    #newValue = float(args[0]['new'])\n",
    "    #print(newValue)\n",
    "    \n",
    "    p = p9.ggplot(cvStringMelt.loc[cvStringMelt['Iter'] == newIter,:], \\\n",
    "        p9.aes(x=\"Image\", y='Distance', color=\"colvar\")) + \\\n",
    "        p9.geom_point() + p9.geom_path() + p9.theme_linedraw() + \\\n",
    "        p9.ylim([0,maxY]) +\\\n",
    "        p9.scale_colour_manual(values = [\"blue\", \"green\", \"red\", \"cyan\", \"magenta\"], \n",
    "                               labels=[ keyNameDic['dist'+str(i)] for i in range(len(keyNameDic)) ])\n",
    "    \n",
    "    if showInitTB.value:\n",
    "        p += p9.geom_point(cvStringMelt.loc[cvStringMelt['Iter'] == 0,:], \\\n",
    "            p9.aes(x=\"Image\", y='Distance', color=\"colvar\"))\n",
    "        p += p9.geom_path(cvStringMelt.loc[cvStringMelt['Iter'] == 0,:], \\\n",
    "            p9.aes(x=\"Image\", y='Distance', color=\"colvar\"), linetype=\"dashed\")\n",
    "            \n",
    "    if fig is None:\n",
    "        fig, plot = p.draw(return_ggplot=True)\n",
    "        axs = plot.axs\n",
    "    else:\n",
    "        for artist in plt.gca().lines +\\\n",
    "                        plt.gca().collections +\\\n",
    "                        plt.gca().artists + plt.gca().patches + plt.gca().texts:\n",
    "            artist.remove()\n",
    "        p._draw_using_figure(fig, axs)\n",
    "    \n",
    "playStringSlider.observe(plotString, 'value')\n",
    "showInitTB.observe(plotString, 'value')\n",
    "\n",
    "display(stringCtl)\n",
    "\n",
    "playStringSlider.value = 0\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will notice that during the optimization process, the values of the harmonic restraints change, making the string \"smoother\" over iteraton steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following interactive plot allows a more precise analysis of the progression of the string optimization. Here you can compare the values of colvars plotted on the *X* and *Y* axis.\n",
    "For chemical reactions, such as the transfer of the Lysine 61 hydrogen atom to the OMP carbon, it is useful to compare two distances: the distance between the hydrogen and the Lysine nitrogen, and the distance between the hydrogen and the OMP carbon. As these distances change, their difference will show the chemical reaction taking place. In our example, the hydrogen from Lysine 61 binds to the carbon that looses the carboxylic group, so the distance between OMP's carbons 6 and 7 is also of interest. By plotting these two measures, we can track the progression of the reaction being studied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- From the **X-Final** dropdown, select *OMP C6 <-> Lys61 HZ1* \n",
    "- From the **X-Initial** dropdown, select *Lys61 NZ <-> Lys61 HZ1*\n",
    "- From the **Y-Final** dropdown, select *OMP C6 <-> OMP C7*\n",
    "- From the **Y-Initial** dropdown, select *None*\n",
    "\n",
    "The images are marked with boxes carrying the image index. In this case (due to this selection of colvars, in this particular order) the plot will show the evolution of of the reaction from the bottom right to the top left. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "pd.set_option('mode.chained_assignment',None)\n",
    "\n",
    "playString  = widgets.Play(\n",
    "    interval=1500,\n",
    "    value=1,\n",
    "    min=0,\n",
    "    max=numIterations-1,\n",
    "    step=1,\n",
    "    description=\"Iterate over the string optimization\",\n",
    ")\n",
    "\n",
    "playStringSlider = widgets.IntSlider(value=1, min=0, max=numIterations-1)\n",
    "widgets.jslink((playString,'value'),(playStringSlider,'value'))\n",
    "\n",
    "stringCtl = widgets.HBox([widgets.Label(value=\"Iteration:\"),playString, playStringSlider])\n",
    "\n",
    "# We create an inverse dictionary to populate the drop-down boxes\n",
    "cvNameKeyDict = { val:key for key,val in keyNameDic.items() }\n",
    "# And add the option of not chosing \n",
    "cvNameKeyDictIni = dict(cvNameKeyDict)\n",
    "cvNameKeyDictIni.update({'None':None})\n",
    "\n",
    "selXi = widgets.Dropdown(\n",
    "    options=cvNameKeyDictIni,\n",
    "    value='dist0',\n",
    "    description='X-Initial',\n",
    ")\n",
    "\n",
    "selXf = widgets.Dropdown(\n",
    "    options=cvNameKeyDict,\n",
    "    value='dist1',\n",
    "    description='X-Final',\n",
    ")\n",
    "\n",
    "dropBoxesX = widgets.HBox([widgets.Label(value=\"X-Diff:\",style={'value_width': 'initial'}),\n",
    "                          selXf,selXi])\n",
    "\n",
    "selYi = widgets.Dropdown(\n",
    "    options=cvNameKeyDictIni,\n",
    "    value='dist2',\n",
    "    description='Y-Initial',\n",
    ")\n",
    "\n",
    "selYf = widgets.Dropdown(\n",
    "    options=cvNameKeyDict,\n",
    "    value='dist3',\n",
    "    description='Y-Final',\n",
    ")\n",
    "\n",
    "dropBoxesY = widgets.HBox([widgets.Label(value=\"Y-Diff:\",style={'value_width': 'initial'}),\n",
    "                          selYf,selYi])\n",
    "\n",
    "fig = None\n",
    "axs = None\n",
    "\n",
    "def plotString(*args):\n",
    "    \n",
    "    global fig, axs\n",
    "    \n",
    "    #newValue = float(args[0]['new'])\n",
    "    #print(newValue)\n",
    "    newValue = float(playStringSlider.value)\n",
    "    \n",
    "    tmpDF = cvString.loc[cvString['Iter'] == newValue,:]\n",
    "    \n",
    "    # Verifies we have a selection for difference between colvars. If not, only one colvar value is shown.\n",
    "    if selXi.value:\n",
    "        tmpDF['xvals'] = tmpDF[selXf.value] - tmpDF[selXi.value]\n",
    "    else:\n",
    "        tmpDF['xvals'] = tmpDF[selXf.value]\n",
    "    \n",
    "    if selYi.value:\n",
    "        tmpDF['yvals'] = tmpDF[selYf.value] - tmpDF[selYi.value]\n",
    "    else:\n",
    "        tmpDF['yvals'] = tmpDF[selYf.value]\n",
    "    \n",
    "    p = p9.ggplot(tmpDF, p9.aes(x=\"xvals\", y=\"yvals\", label=\"Image\")) + \\\n",
    "        p9.geom_point() + p9.geom_path() + p9.theme_linedraw() +\\\n",
    "        p9.labs(x=\"Colvar Value (A)\", y=\"Colvar Value (A)\") + p9.geom_label()\n",
    "            \n",
    "    if fig is None:\n",
    "        fig, plot = p.draw(return_ggplot=True)\n",
    "        axs = plot.axs\n",
    "    else:\n",
    "        for artist in plt.gca().lines +\\\n",
    "                        plt.gca().collections +\\\n",
    "                        plt.gca().artists + plt.gca().patches + plt.gca().texts:\n",
    "            artist.remove()\n",
    "        p._draw_using_figure(fig, axs)\n",
    "    \n",
    "playStringSlider.observe(plotString, 'value')\n",
    "\n",
    "selXf.observe(plotString, \"value\")\n",
    "selXi.observe(plotString, \"value\")\n",
    "selYf.observe(plotString, \"value\")\n",
    "selYi.observe(plotString, \"value\")\n",
    "\n",
    "display(stringCtl)\n",
    "display(dropBoxesX)\n",
    "display(dropBoxesY)\n",
    "\n",
    "playStringSlider.value = 0\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot above, the 0 (zero) in the X axis indicates the point where the lysine hydrogen atom is at equal distance from the lysine nitrogen, and the OMP carbon. \n",
    "You will notice that, after the string optimization, the final reaction coordinate **first** separates the carbon atoms C6 and C7 (seen by the increased distance between them in the Y axis from image 4 through 8), **and then** transfer the hydrogen atom (seen by the  change in X value, from image 9 to 11)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# eABF preparation and execution\n",
    "\n",
    "The files for eABF have already been prepared in a previous step. Now, we just create output directories and copy *coordinate* and *periodic box* files from the results of the String Optimization. From the coordinate files, we create PDB files that will define the string path over which multiple ABF walkers will \"walk\". In order to optimize the sampling, walker are initialized *along* the string, using the optimized string images as starting points, instead of the classical approach where all walkers are initialized at the middle of the string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gets the output path, but removes the last \"slash\" so a new output directory can be created for the eABF results.\n",
    "if outputPath[-1] == \"/\":\n",
    "    tmpOutputPath = outputPath[0:len(outputPath)-1]\n",
    "\n",
    "# Creates individual output directories for each eABF walker.\n",
    "for i in range(numImgs*eABFWalkersPerImg):\n",
    "    os.makedirs(tmpOutputPath + \"_eABF/\" + str(i),exist_ok=True)\n",
    "\n",
    "# Creates the directory where we will place the results of the string optimization.\n",
    "os.makedirs(basePath + \"/string\", exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell is a BASH cell that simply runs two bash scripts that collect the string optimization results (copies files to the \"string\" folder), and transforms the final namdbin files into PDB files (runs \"catdcd\" on the some of the files). They could be written entirely in python into the cell, or as a bash script directly into a BASH cell, but we chose to keep this code in separate files to keep the notebook \"small\", since there is nothing special about this code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$numImgs\"  \"$catdcdPath\" \"$psfFile\"\n",
    "\n",
    "cd String_eABF/\n",
    "\n",
    "../toolkit/collect-configuration $1 decarb\n",
    "\n",
    "cd string\n",
    "\n",
    "for (( i = 0; i < ${1}; i+=1))\n",
    "do\n",
    "\n",
    "    eval \"$2\" -o string_${i}.pdb -otype pdb -stype psf -s ../../${3} -namdbin string_${i}.coor \n",
    "\n",
    "done\n",
    "\n",
    "cd ../../"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run eABF simulation\n",
    "\n",
    "Once more, pre-ran results have been provided in the *example-output* folder. You can copy those results into the *String_eABF/output_eABF* folder and proceed with the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$numImgs\"  \"$namdPathMPI\"\n",
    "\n",
    "cd String_eABF/\n",
    "\n",
    "#mpirun -n $1 $2 +replicas $1 eABF.tcl +stdout output_eABF/%1d/job00.%1d.log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running the simulation with multiple walkers, one must combine all the outputs into one, so that the final calculated PFM is the result of all sampling done by all walkers.\n",
    "The easiest way to do that is by running a new colculation that takes all previous results as inputs, and then combines them. The following cell runs NAMD with a script (*eABF_combine.tcl*) that has been designed for that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash -s \"$numImgs\"  \"$namdPathMPI\"\n",
    "\n",
    "cd String_eABF/\n",
    "\n",
    "#mpirun -n 1 $2 +replicas 1 eABF_combine.tcl +stdout output_eABF/%1d/job00_combine.%1d.log"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# eABF analysis\n",
    "\n",
    "The following plot combines the colvar transformations with the eABF results, showing how the PMF changes along the reaction coordinate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read PMF from czar file\n",
    "\n",
    "eABFres = pd.DataFrame(columns=['bin','pmf'])\n",
    "\n",
    "regexBinPMF = re.compile(\"^([0-9\\.]*) {1,}([0-9\\.]*)\")\n",
    "\n",
    "with open(\"./String_eABF/output_eABF/0/decarb_combined0.job0000.pmf\",'r') as infile:\n",
    "    lineIter = 0\n",
    "    for line in infile:\n",
    "        if line.startswith(\"#\") or not line.strip():\n",
    "            continue\n",
    "        \n",
    "        matchBinPMF = regexBinPMF.match(line.strip())\n",
    "        \n",
    "        eABFres.loc[lineIter,\"bin\"] = float(matchBinPMF.group(1))\n",
    "        eABFres.loc[lineIter,\"pmf\"] = round(float(matchBinPMF.group(2)), 2)\n",
    "        \n",
    "        lineIter += 1        \n",
    "\n",
    "eABFres['pmf'] = eABFres.pmf.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot free energy\n",
    "\n",
    "p9.ggplot(eABFres, p9.aes(x='bin', y='pmf', group=1))  + p9.geom_point() + p9.geom_path() +\\\n",
    "p9.theme_light() + p9.labs(x=\"Reaction Coordinate\", y = \"PMF (kcal/mol)\") +\\\n",
    "p9.theme(axis_text_x = p9.element_text(angle = -45, hjust = 0.5)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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": 2
}