docs/tutorials/02-tutorial-cb6-but-pbc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pAPRika tutorial 2 - Explicit Solvent\n",
    "\n",
    "In this example, we will setup and simulate butane (BUT) as a guest molecule for the host [cucurbit[6]uril](https://en.wikipedia.org/wiki/Cucurbituril) (CB6). \n",
    "\n",
    "Unlike the previous tutorial, in this one, we will solvate the host-guest complex in water with additional Na<sup>+</sup> and Cl<sup>-</sup> ions. This tutorial builds off tutorial 1. I've marked new pieces of information with a blue circle, like so:\n",
    "\n",
    "🔵 Since we have a prepared host-guest-dummy setup from the first tutorial, I'm going to skip the initial `tleap` steps and go right into initializing the restraints. \n",
    "\n",
    "🔵 The default mode of *pAPRika* is to detect whether a simulation has completed successfully, and if so, skip that window. This helps when launching mulitple jobs on a cluster or queuing system that might kill jobs after a pre-determined time. That means that running this notebook with the same set of windows as the first tutorial will not re-run the simulations and the results won't change. The easiest thing to do is to rename or delete the existing `windows` directory, if you would like to compare the results between the two methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "data = \"../../../paprika/data/cb6-but\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine the number of windows\n",
    "\n",
    "Before we add the restraints, it is helpful to set the $\\lambda$ fractions that control the strength of the force constants during attach and release, and to define the distances for the pulling phase.\n",
    "\n",
    "The attach fractions go from 0 to 1 and we place more points at the bottom of the range to sample the curvature of $dU/d \\lambda$. Next, we generally apply a distance restraint until the guest is ~18 Angstroms away from the host, in increments of 0.4 Angstroms. This distance should be at least twice the Lennard-Jones cutoff in the system. These values have worked well for us, but this is one aspect that should be carefully checked for new systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:14.528611Z",
     "start_time": "2018-10-18T01:18:14.522595Z"
    }
   },
   "outputs": [],
   "source": [
    "attach_string = \"0.00 0.40 0.80 1.60 2.40 4.00 5.50 8.65 11.80 18.10 24.40 37.00 49.60 74.80 100.00\"\n",
    "attach_fractions = [float(i) / 100 for i in attach_string.split()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:24.616399Z",
     "start_time": "2018-10-18T01:18:24.610179Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "initial_distance = 6.0\n",
    "pull_distances = np.arange(0.0 + initial_distance, 18.0 + initial_distance, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These values will be used to measure distance relative to the first dummy atom, hence the addition of `6.00`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:24.924746Z",
     "start_time": "2018-10-18T01:18:24.922190Z"
    }
   },
   "outputs": [],
   "source": [
    "release_fractions = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:25.245830Z",
     "start_time": "2018-10-18T01:18:25.242238Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are [15, 18, 0] windows in this attach-pull-release calculation.\n"
     ]
    }
   ],
   "source": [
    "windows = [len(attach_fractions), len(pull_distances), len(release_fractions)]\n",
    "print(f\"There are {windows} windows in this attach-pull-release calculation.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could specify the number of windows for each phase and the force constants and targets will be linearly interpolated. Other ways of specifying these values are documented in the code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add restraints using pAPRika\n",
    "\n",
    "`pAPRika` supports four different types of restraints:\n",
    "\n",
    "\n",
    "- Static restraints: these six restraints keep the host and in the proper orientation during the simulation (necessary),\n",
    "- Guest restraints: these restraints pull the guest away from the host along the $z$-axis (necessary),\n",
    "- Conformational restraints: these restraints alter the conformational sampling of the host molecule (optional), and\n",
    "- Wall restraints: these restraints help define the bound state of the guest (optional).\n",
    "\n",
    "More information on these restraints can be found in:\n",
    "\n",
    "Henriksen, N.M., Fenley, A.T., and Gilson, M.K. (2015). Computational Calorimetry: High-Precision Calculation of Host-Guest Binding Thermodynamics. J. Chem. Theory Comput. 11, 4377–4394. [DOI](http://doi.org/10.1021/acs.jctc.5b00405)\n",
    "\n",
    "In this example, I will show how to setup the static restraints and the guest restraints.\n",
    "\n",
    "We have already added the dummy atoms and we have already defined the guest anchor atoms. Now we need to define the host anchor atoms (H1, H2, and H3) in the above diagram. The host anchors should be heavy atoms distributed around the cavity (and around the pulling axis). One caveat is that the host anchors should be rigid relative to each other, so conformational restraints do not shift the alignment of the pulling axis relative to the solvation box. For CB6, I have chosen carbons around the central ridge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:25.912392Z",
     "start_time": "2018-10-18T01:18:25.906786Z"
    }
   },
   "outputs": [],
   "source": [
    "G1 = \":BUT@C\"\n",
    "G2 = \":BUT@C3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:26.269326Z",
     "start_time": "2018-10-18T01:18:26.266155Z"
    }
   },
   "outputs": [],
   "source": [
    "H1 = \":CB6@C\"\n",
    "H2 = \":CB6@C31\"\n",
    "H3 = \":CB6@C18\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll also make a shorthand for the dummy atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:30.137541Z",
     "start_time": "2018-10-18T01:18:30.131221Z"
    }
   },
   "outputs": [],
   "source": [
    "D1 = \":DM1\"\n",
    "D2 = \":DM2\"\n",
    "D3 = \":DM3\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Static restraints\n",
    "\n",
    "These harmonic restraints are constant throughout the entire simulation. These restraints are used to control the distances and angles between the host and guest relative to the dummy atom. We have created a special class for these restrains, `static_DAT_restraint`, that uses the initial value as the restraint target (this is why the starting structure should be a reasonable facsimile of the bound state).\n",
    "\n",
    "Note that these restraints are not \"attached\" and they don't need to be \"released\" -- their force constants do not change in magnitude.\n",
    "\n",
    "The first three static restraints affect the translational distance, angle, and torsion angle between the host and the dummy atoms. These control the position of the host, via the first anchor atom, from moving relative to the dummy atoms.\n",
    "\n",
    "There is no *correct* value for the force constants. From experience, we know that a distance force constant of 5.0 kcal/mol/$\\overset{\\circ}{\\text{A}}$$^2$ won't nail down the host and yet it also won't wander away. Likewise, we have had good results using 100.0 kcal/mol/$rad^2$ for the angle force constant.\n",
    "\n",
    "![static1](images/static-restraints-1.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:32.038255Z",
     "start_time": "2018-10-18T01:18:32.033456Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika import restraints\n",
    "static_restraints = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:18:59.052016Z",
     "start_time": "2018-10-18T01:18:59.014498Z"
    }
   },
   "outputs": [],
   "source": [
    "import parmed as pmd\n",
    "structure = pmd.load_file(\"complex/cb6-but-dum.prmtop\", \"complex/cb6-but-dum.rst7\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:00.187585Z",
     "start_time": "2018-10-18T01:19:00.121551Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D1, H1],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 5.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:00.644783Z",
     "start_time": "2018-10-18T01:19:00.602812Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D2, D1, H1],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 100.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:01.059199Z",
     "start_time": "2018-10-18T01:19:01.018767Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D3, D2, D1, H1],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 100.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next three restraints control the orientation of the host relative to the dummy atoms. These angle and torsion restraints prevent the host from rotating relative to the dummy atoms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "![static2](images/static-restraints-2.png)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:02.672360Z",
     "start_time": "2018-10-18T01:19:02.611492Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D1, H1, H2],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 100.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:03.442798Z",
     "start_time": "2018-10-18T01:19:03.378864Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D2, D1, H1, H2],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 100.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:03.865544Z",
     "start_time": "2018-10-18T01:19:03.819890Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.static_DAT_restraint(restraint_mask_list = [D1, H1, H2, H3],\n",
    "                                    num_window_list = windows,\n",
    "                                    ref_structure = structure,\n",
    "                                    force_constant = 100.0,\n",
    "                                    amber_index=True)\n",
    "\n",
    "static_restraints.append(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Guest restraints\n",
    "\n",
    "Next, we add restraints on the guest. These restraints control the position of the guest and are the key to the attach-pull-release method. During the attach phase, the *force constants* for these restraints is increased from zero. During the pull phase, the *target* for the distance restraint is increased (in the orange box, below), translating the guest away from the host cavity. And during the release phase, the *force constants* are reduced from their \"full\" value back down to zero.\n",
    "\n",
    "![guest](images/guest-restraints.png)\n",
    "\n",
    "We use the class `DAT_restraint` to create these three restraints. We will use the same anchor atoms as before, with the same distance and angle force constants. Note that unlike `static_DAT_restraint`, we will first create the restraint, *then* set the attributes, *then* initialize the restraint which does some checks to make sure everything is copacetic.\n",
    "\n",
    "There are two additional convenience options here:\n",
    "\n",
    "- `auto_apr = True` sets the force constant during pull to be the final force constant after attach and sets the initial restraint target during pull to be the final attach target.\n",
    "\n",
    "- `continuous_apr = True` sets the last window of attach to be the same as the first window as pull (and likewise for release)\n",
    "\n",
    "Also note, due to a quirk with AMBER, we specifcy angle and torsion targets in degrees but the force constant using radians!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:04.952877Z",
     "start_time": "2018-10-18T01:19:04.947610Z"
    }
   },
   "outputs": [],
   "source": [
    "guest_restraints = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:05.810565Z",
     "start_time": "2018-10-18T01:19:05.776595Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.DAT_restraint()\n",
    "r.mask1 = D1\n",
    "r.mask2 = G1\n",
    "r.topology = structure\n",
    "r.auto_apr = True\n",
    "r.continuous_apr = True\n",
    "r.amber_index = True\n",
    "\n",
    "r.attach[\"target\"] = 6.0                            # Angstroms\n",
    "r.attach[\"fraction_list\"] = attach_fractions\n",
    "r.attach[\"fc_final\"] = 5.0                          # kcal/mol/Angstroms**2\n",
    "\n",
    "r.pull[\"target_final\"] = 24.0                       # Angstroms\n",
    "r.pull[\"num_windows\"] = windows[1]\n",
    "\n",
    "r.initialize()\n",
    "guest_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:06.987350Z",
     "start_time": "2018-10-18T01:19:06.947188Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.DAT_restraint()\n",
    "r.mask1 = D2\n",
    "r.mask2 = D1\n",
    "r.mask3 = G1\n",
    "r.topology = structure\n",
    "r.auto_apr = True\n",
    "r.continuous_apr = True\n",
    "r.amber_index = True\n",
    "\n",
    "r.attach[\"target\"] = 180.0                          # Degrees\n",
    "r.attach[\"fraction_list\"] = attach_fractions\n",
    "r.attach[\"fc_final\"] = 100.0                        # kcal/mol/radian**2\n",
    "\n",
    "r.pull[\"target_final\"] = 180.0                      # Degrees\n",
    "r.pull[\"num_windows\"] = windows[1]\n",
    "\n",
    "r.initialize()\n",
    "guest_restraints.append(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:07.409139Z",
     "start_time": "2018-10-18T01:19:07.385980Z"
    }
   },
   "outputs": [],
   "source": [
    "r = restraints.DAT_restraint()\n",
    "r.mask1 = D1\n",
    "r.mask2 = G1\n",
    "r.mask3 = G2\n",
    "r.topology = structure\n",
    "r.auto_apr = True\n",
    "r.continuous_apr = True\n",
    "r.amber_index = True\n",
    "\n",
    "r.attach[\"target\"] = 180.0                          # Degrees\n",
    "r.attach[\"fraction_list\"] = attach_fractions\n",
    "r.attach[\"fc_final\"] = 100.0                        # kcal/mol/radian**2\n",
    "\n",
    "r.pull[\"target_final\"] = 180.0                      # Degrees\n",
    "r.pull[\"num_windows\"] = windows[1]\n",
    "\n",
    "r.initialize()\n",
    "guest_restraints.append(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create the directory structure\n",
    "\n",
    "We use the guest restraints to create a list of windows with the appropriate names and then we create the directories."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:08.673439Z",
     "start_time": "2018-10-18T01:19:08.667229Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from paprika.restraints.restraints import create_window_list\n",
    "\n",
    "window_list = create_window_list(guest_restraints)\n",
    "for window in window_list:\n",
    "    os.makedirs(f\"windows/{window}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, we ask pAPRika to write the restraint information in each window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In each window, we create a file named `disang.rest` to hold all of the restraint information that is required to run the simulation with AMBER. We feed the list of restraints to pAPRika, one by one, and it returns the appropriate line.\n",
    "\n",
    "The functional form of the restraints is specified in section 25.1 of the AMBER18 manual. Specifically, these restraints have a square bottom with parabolic sides out to a specific distance and then linear sides beyond that. The square bottom can be eliminated by setting `r2=r3` and the linear extension can be eliminated by setting the the `r4 = 999` and `r1 = 0`, creating a harmonic restraint "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:10.707258Z",
     "start_time": "2018-10-18T01:19:10.674546Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.restraints.amber import amber_restraint_line\n",
    "\n",
    "host_guest_restraints = (static_restraints + guest_restraints)\n",
    "for window in window_list:\n",
    "    with open(f\"windows/{window}/disang.rest\", \"w\") as file:\n",
    "        for restraint in host_guest_restraints:\n",
    "            string = amber_restraint_line(restraint, window)\n",
    "            if string is not None:\n",
    "                file.write(string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is a good idea to open up the `disang.rest` files and see that the force constants and targets make sense (and there are no `NaN` values). Do the force constants for the guest restraints start at zero? Do the targets for the pull slowly increase?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save restraints to a JSON file\n",
    "\n",
    "🔵 The previous code block writes restraints in a format that *AMBER* uses to run simulations in each window, but *pAPRika* uses other information stored in the `DAT_restraint` class in the analysis module. By saving the restraints in a file, it is easy to chunk up setup, simulation, and analysis in different sessions or*python* scripts.\n",
    "\n",
    "Note: the API for saving and loading restraints is not stable and may be moved into another submodule. The output notes that *pAPRika* is not including the entire topology in the restraint JSON, but will contain a reference to the file that was used to setup th restraints, in this case, `cb6-but-dum.prmtop`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:14.391701Z",
     "start_time": "2018-10-18T01:19:14.386325Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.io import save_restraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:18.453177Z",
     "start_time": "2018-10-18T01:19:18.448137Z"
    }
   },
   "outputs": [],
   "source": [
    "save_restraints(host_guest_restraints, filepath=\"windows/restraints.json\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare host-guest system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Translation of the guest (before solvation)\n",
    "\n",
    "For the attach windows, we will use the initial, bound coordinates for the host-guest complex. Only the force constants change during this phase, so a single set of coordinates is sufficient. For the pull windows, we will translate the guest to the target value of the restraint before solvation, and for the release windows, we will use the coordinates from the final pull window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:29.723830Z",
     "start_time": "2018-10-18T01:19:28.670921Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In window p000 we will translate the guest 0.0 Angstroms.\n",
      "In window p001 we will translate the guest 1.1 Angstroms.\n",
      "In window p002 we will translate the guest 2.1 Angstroms.\n",
      "In window p003 we will translate the guest 3.2 Angstroms.\n",
      "In window p004 we will translate the guest 4.2 Angstroms.\n",
      "In window p005 we will translate the guest 5.3 Angstroms.\n",
      "In window p006 we will translate the guest 6.4 Angstroms.\n",
      "In window p007 we will translate the guest 7.4 Angstroms.\n",
      "In window p008 we will translate the guest 8.5 Angstroms.\n",
      "In window p009 we will translate the guest 9.5 Angstroms.\n",
      "In window p010 we will translate the guest 10.6 Angstroms.\n",
      "In window p011 we will translate the guest 11.6 Angstroms.\n",
      "In window p012 we will translate the guest 12.7 Angstroms.\n",
      "In window p013 we will translate the guest 13.8 Angstroms.\n",
      "In window p014 we will translate the guest 14.8 Angstroms.\n",
      "In window p015 we will translate the guest 15.9 Angstroms.\n",
      "In window p016 we will translate the guest 16.9 Angstroms.\n",
      "In window p017 we will translate the guest 18.0 Angstroms.\n"
     ]
    }
   ],
   "source": [
    "import shutil\n",
    "\n",
    "for window in window_list:\n",
    "    if window[0] == \"a\":\n",
    "        shutil.copy(\"complex/cb6-but-dum.prmtop\", f\"windows/{window}/cb6-but-dum.prmtop\")\n",
    "        shutil.copy(\"complex/cb6-but-dum.rst7\", f\"windows/{window}/cb6-but-dum.rst7\")\n",
    "    elif window[0] == \"p\":\n",
    "        structure = pmd.load_file(\"complex/cb6-but-dum.prmtop\", \"complex/cb6-but-dum.rst7\",\n",
    "                          structure = True)\n",
    "        target_difference = guest_restraints[0].phase['pull']['targets'][int(window[1:])] - guest_restraints[0].pull['target_initial']\n",
    "        print(f\"In window {window} we will translate the guest {target_difference.magnitude:0.1f}.\")\n",
    "        for atom in structure.atoms:\n",
    "            if atom.residue.name == \"BUT\":\n",
    "                atom.xz += target_difference.magnitude\n",
    "        structure.save(f\"windows/{window}/cb6-but-dum.prmtop\")\n",
    "        structure.save(f\"windows/{window}/cb6-but-dum.rst7\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before running a simulation, open each window and check that the position of the host and dummy atoms are fixed and that the position of the guest is bound during the attach windows, moves progressively further during pull, and is away from the host during the release windows."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solvate each window\n",
    "\n",
    "🔵 Unlike the previous tutorial we will solvate the system in each window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:19:33.178637Z",
     "start_time": "2018-10-18T01:19:33.172333Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.build.system import TLeap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, this will use cubic periodic boundary conditions; we usually use octahedral periodic boundary conditions to optimize space. Based on experience, 2000 waters is a resaonable amount for CB6, but sometimes it is necessary to run with 2500 or 3000 waters depending on the size of the host and guest molecules. The code below first neutralizes the system (which does not do anything here) and then adds an *additional* 6 Na+ and 6 Cl- ions. The counterion species can be changed via `counter_cation` and `counter_anion`, and the `add_ions` string can take concentrations in the form of `0.150M` or `0.150m` for Molarity or molality, respectively.\n",
    "\n",
    "Note that when running simulations with periodic boundary conditions, be mindful of the box size relative to the Lennard-Jones cutoff."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:21:32.919048Z",
     "start_time": "2018-10-18T01:19:34.659504Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solvating system in window a000.\n",
      "Solvating system in window a001.\n",
      "Solvating system in window a002.\n",
      "Solvating system in window a003.\n",
      "Solvating system in window a004.\n",
      "Solvating system in window a005.\n",
      "Solvating system in window a006.\n",
      "Solvating system in window a007.\n",
      "Solvating system in window a008.\n",
      "Solvating system in window a009.\n",
      "Solvating system in window a010.\n",
      "Solvating system in window a011.\n",
      "Solvating system in window a012.\n",
      "Solvating system in window a013.\n",
      "Solvating system in window p000.\n",
      "Solvating system in window p001.\n",
      "Solvating system in window p002.\n",
      "Solvating system in window p003.\n",
      "Solvating system in window p004.\n",
      "Solvating system in window p005.\n",
      "Solvating system in window p006.\n",
      "Solvating system in window p007.\n",
      "Solvating system in window p008.\n",
      "Solvating system in window p009.\n",
      "Solvating system in window p010.\n",
      "Solvating system in window p011.\n",
      "Solvating system in window p012.\n",
      "Solvating system in window p013.\n",
      "Solvating system in window p014.\n",
      "Solvating system in window p015.\n",
      "Solvating system in window p016.\n",
      "Solvating system in window p017.\n"
     ]
    }
   ],
   "source": [
    "for window in window_list:\n",
    "    print(f\"Solvating system in window {window}.\")\n",
    "    structure = pmd.load_file(f\"windows/{window}/cb6-but-dum.prmtop\",\n",
    "                              f\"windows/{window}/cb6-but-dum.rst7\")\n",
    "    \n",
    "    if not os.path.exists(f\"windows/{window}/cb6-but-dum.pdb\"):\n",
    "        structure.save(f\"windows/{window}/cb6-but-dum.pdb\")\n",
    "    \n",
    "    system = TLeap()\n",
    "    system.output_path = os.path.join(\"windows\", window)\n",
    "    system.output_prefix = \"cb6-but-dum-sol\"\n",
    "    \n",
    "    system.target_waters = 2000\n",
    "    system.neutralize = True\n",
    "    system.add_ions = [\"Na+\", 6, \"Cl-\", 6]\n",
    "    system.template_lines = [\n",
    "        \"source leaprc.gaff\",\n",
    "        \"source leaprc.water.tip3p\",\n",
    "        f\"loadamberparams {data}/cb6.frcmod\",\n",
    "        \"loadamberparams ../../complex/dummy.frcmod\",\n",
    "        f\"CB6 = loadmol2 {data}/cb6.mol2\",\n",
    "        f\"BUT = loadmol2 {data}/but.mol2\",\n",
    "        \"DM1 = loadmol2 ../../complex/dm1.mol2\",\n",
    "        \"DM2 = loadmol2 ../../complex/dm2.mol2\",\n",
    "        \"DM3 = loadmol2 ../../complex/dm3.mol2\",\n",
    "        \"model = loadpdb cb6-but-dum.pdb\",\n",
    "    ]\n",
    "    system.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "![solvated](images/cb6-but-dum-sol-diagram.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation\n",
    "\n",
    "🔵 This is going to look very similar to the the implicit solvent simulation setup, except instead of using `config_gb_min` and `config_gb_sim`, we will use `config_pdb_min` and `config_pbc_sim`. \n",
    "\n",
    "For this part, you need to have the AMBER executables in your path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:21:32.924748Z",
     "start_time": "2018-10-18T01:21:32.921317Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.simulate import AMBER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a quick minimization in every window. Note that we need to specify `simulation.cntrl[\"ntr\"] = 1` to enable the positional restraints on the dummy atoms.\n",
    "\n",
    "I'm using the `logging` module to keep track of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:21:33.016051Z",
     "start_time": "2018-10-18T01:21:32.926515Z"
    }
   },
   "outputs": [],
   "source": [
    "import logging\n",
    "from importlib import reload\n",
    "reload(logging)\n",
    "\n",
    "logger = logging.getLogger()\n",
    "logging.basicConfig(\n",
    "    format='%(asctime)s %(message)s',\n",
    "    datefmt='%Y-%m-%d %I:%M:%S %p',\n",
    "    level=logging.INFO\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:24:19.474204Z",
     "start_time": "2018-10-18T01:21:33.018199Z"
    },
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Energy Minimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-10-05 12:18:34 PM Running minimization in window a000...\n",
      "2020-10-05 12:19:00 PM Running minimization in window a001...\n",
      "2020-10-05 12:19:26 PM Running minimization in window a002...\n",
      "2020-10-05 12:19:52 PM Running minimization in window a003...\n",
      "2020-10-05 12:20:18 PM Running minimization in window a004...\n",
      "2020-10-05 12:20:45 PM Running minimization in window a005...\n",
      "2020-10-05 12:21:11 PM Running minimization in window a006...\n",
      "2020-10-05 12:21:37 PM Running minimization in window a007...\n",
      "2020-10-05 12:22:03 PM Running minimization in window a008...\n",
      "2020-10-05 12:22:29 PM Running minimization in window a009...\n",
      "2020-10-05 12:22:56 PM Running minimization in window a010...\n",
      "2020-10-05 12:23:22 PM Running minimization in window a011...\n",
      "2020-10-05 12:23:48 PM Running minimization in window a012...\n",
      "2020-10-05 12:24:14 PM Running minimization in window a013...\n",
      "2020-10-05 12:24:40 PM Running minimization in window p000...\n",
      "2020-10-05 12:25:07 PM Running minimization in window p001...\n",
      "2020-10-05 12:25:33 PM Running minimization in window p002...\n",
      "2020-10-05 12:25:59 PM Running minimization in window p003...\n",
      "2020-10-05 12:26:25 PM Running minimization in window p004...\n",
      "2020-10-05 12:26:51 PM Running minimization in window p005...\n",
      "2020-10-05 12:27:17 PM Running minimization in window p006...\n",
      "2020-10-05 12:27:44 PM Running minimization in window p007...\n",
      "2020-10-05 12:28:10 PM Running minimization in window p008...\n",
      "2020-10-05 12:28:36 PM Running minimization in window p009...\n",
      "2020-10-05 12:29:02 PM Running minimization in window p010...\n",
      "2020-10-05 12:29:28 PM Running minimization in window p011...\n",
      "2020-10-05 12:29:54 PM Running minimization in window p012...\n",
      "2020-10-05 12:30:20 PM Running minimization in window p013...\n",
      "2020-10-05 12:30:47 PM Running minimization in window p014...\n",
      "2020-10-05 12:31:13 PM Running minimization in window p015...\n",
      "2020-10-05 12:31:39 PM Running minimization in window p016...\n",
      "2020-10-05 12:32:05 PM Running minimization in window p017...\n"
     ]
    }
   ],
   "source": [
    "for window in window_list:\n",
    "    simulation = AMBER()\n",
    "    simulation.executable = \"mpirun -np 4 pmemd.MPI\"\n",
    "\n",
    "    simulation.path = f\"windows/{window}/\"\n",
    "    simulation.prefix = \"minimize\"\n",
    "\n",
    "    simulation.topology = \"cb6-but-dum-sol.prmtop\"\n",
    "    simulation.coordinates = \"cb6-but-dum-sol.rst7\"\n",
    "    simulation.ref = \"cb6-but-dum-sol.rst7\"\n",
    "    simulation.restraint_file = \"disang.rest\"\n",
    "\n",
    "    simulation.config_pbc_min()\n",
    "    simulation.cntrl[\"ntr\"] = 1\n",
    "    simulation.cntrl[\"restraint_wt\"] = 50.0\n",
    "    simulation.cntrl[\"restraintmask\"] = \"'@DUM'\"\n",
    "\n",
    "    logging.info(f\"Running minimization in window {window}...\")\n",
    "    simulation.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For simplicity, I am going to skip equilibration and go straight to production!\n",
    "\n",
    "🔵 I'm also going to run a little bit longer in each window, because the additional atoms increases the convergence time. The total simulation time per window, in nanoseconds, is (nstlim) × (dt in fs / 1000) = (50000) × (0.002 / 1000) = 0.1 ns. This is still very short and just for demonstration. A production simulation might require as much as 15 ns per window for binding free energies and up to 1 us per window for binding enthalpies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Production Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:22.448474Z",
     "start_time": "2018-10-18T01:24:19.476215Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-10-05 12:42:12 PM Running production in window a000...\n",
      "2020-10-05 12:42:25 PM Running production in window a001...\n",
      "2020-10-05 12:42:38 PM Running production in window a002...\n",
      "2020-10-05 12:42:51 PM Running production in window a003...\n",
      "2020-10-05 12:43:04 PM Running production in window a004...\n",
      "2020-10-05 12:43:18 PM Running production in window a005...\n",
      "2020-10-05 12:43:31 PM Running production in window a006...\n",
      "2020-10-05 12:43:44 PM Running production in window a007...\n",
      "2020-10-05 12:43:57 PM Running production in window a008...\n",
      "2020-10-05 12:44:11 PM Running production in window a009...\n",
      "2020-10-05 12:44:24 PM Running production in window a010...\n",
      "2020-10-05 12:44:37 PM Running production in window a011...\n",
      "2020-10-05 12:44:50 PM Running production in window a012...\n",
      "2020-10-05 12:45:04 PM Running production in window a013...\n",
      "2020-10-05 12:45:17 PM Running production in window p000...\n",
      "2020-10-05 12:45:30 PM Running production in window p001...\n",
      "2020-10-05 12:45:43 PM Running production in window p002...\n",
      "2020-10-05 12:45:57 PM Running production in window p003...\n",
      "2020-10-05 12:46:10 PM Running production in window p004...\n",
      "2020-10-05 12:46:23 PM Running production in window p005...\n",
      "2020-10-05 12:46:36 PM Running production in window p006...\n",
      "2020-10-05 12:46:50 PM Running production in window p007...\n",
      "2020-10-05 12:47:03 PM Running production in window p008...\n",
      "2020-10-05 12:47:16 PM Running production in window p009...\n",
      "2020-10-05 12:47:29 PM Running production in window p010...\n",
      "2020-10-05 12:47:43 PM Running production in window p011...\n",
      "2020-10-05 12:47:56 PM Running production in window p012...\n",
      "2020-10-05 12:48:09 PM Running production in window p013...\n",
      "2020-10-05 12:48:22 PM Running production in window p014...\n",
      "2020-10-05 12:48:35 PM Running production in window p015...\n",
      "2020-10-05 12:48:49 PM Running production in window p016...\n",
      "2020-10-05 12:49:02 PM Running production in window p017...\n"
     ]
    }
   ],
   "source": [
    "for window in window_list:\n",
    "    simulation = AMBER()\n",
    "    simulation.executable = \"pmemd.cuda\"\n",
    "    simulation.gpu_devices = 0\n",
    "    \n",
    "    simulation.path = f\"windows/{window}/\"\n",
    "    simulation.prefix = \"production\"\n",
    "\n",
    "    simulation.topology = \"cb6-but-dum-sol.prmtop\"\n",
    "    simulation.coordinates = \"minimize.rst7\"\n",
    "    simulation.ref = \"cb6-but-dum-sol.rst7\"\n",
    "    simulation.restraint_file = \"disang.rest\"\n",
    "\n",
    "    simulation.config_pbc_md()\n",
    "    simulation.cntrl[\"ntr\"] = 1\n",
    "    simulation.cntrl[\"restraint_wt\"] = 50.0\n",
    "    simulation.cntrl[\"restraintmask\"] = \"'@DUM'\"\n",
    "    simulation.cntrl[\"nstlim\"] = 50000\n",
    "    \n",
    "    logging.info(f\"Running production in window {window}...\")\n",
    "    simulation.run(overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis \n",
    "\n",
    "🔵 We will perform the analysis using restraints saved previously to a JSON file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the simulation is completed, we can using the `analysis` module to determine the binding free energy. We supply the location of the parameter information, a string or list for the file names (wildcards supported), the location of the windows, and the restraints on the guest.\n",
    "\n",
    "In this example, we use the method `ti-block` which determines the free energy using **t**hermodynamic **i**integration and then estimates the standard error of the mean at each data point using blocking analysis. Bootstrapping it used to determine the uncertainty of the  full thermodynamic integral for each phase.\n",
    "\n",
    "After running `compute_free_energy()`, a dictionary called `results` will be populated, that contains the free energy and SEM for each phase of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:22.453521Z",
     "start_time": "2018-10-18T01:40:22.450575Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.io import load_restraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:22.729689Z",
     "start_time": "2018-10-18T01:40:22.455341Z"
    }
   },
   "outputs": [],
   "source": [
    "guest_restraints = load_restraints(\"windows/restraints.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:23.032641Z",
     "start_time": "2018-10-18T01:40:22.731725Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika import analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:33.699925Z",
     "start_time": "2018-10-18T01:40:23.034510Z"
    }
   },
   "outputs": [],
   "source": [
    "free_energy = analysis.fe_calc()\n",
    "free_energy.topology = \"cb6-but-dum-sol.prmtop\"\n",
    "free_energy.trajectory = 'production*.nc'\n",
    "free_energy.path = \"windows\"\n",
    "free_energy.restraint_list = guest_restraints\n",
    "free_energy.collect_data()\n",
    "free_energy.methods = ['ti-block']\n",
    "free_energy.ti_matrix = \"full\"\n",
    "free_energy.boot_cycles = 1000\n",
    "free_energy.compute_free_energy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:33.767780Z",
     "start_time": "2018-10-18T01:40:33.701759Z"
    }
   },
   "outputs": [],
   "source": [
    "free_energy.compute_ref_state_work([\n",
    "    guest_restraints[0], guest_restraints[1], None, \n",
    "    None, guest_restraints[2], None\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:33.811304Z",
     "start_time": "2018-10-18T01:40:33.769779Z"
    }
   },
   "outputs": [],
   "source": [
    "binding_affinity = -1 * (\n",
    "free_energy.results[\"attach\"][\"ti-block\"][\"fe\"] + \\\n",
    "free_energy.results[\"pull\"][\"ti-block\"][\"fe\"] + \\\n",
    "free_energy.results[\"ref_state_work\"]\n",
    ")\n",
    "\n",
    "sem = np.sqrt(\n",
    "free_energy.results[\"attach\"][\"ti-block\"][\"sem\"]**2 + \\\n",
    "free_energy.results[\"pull\"][\"ti-block\"][\"sem\"]**2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:40:33.961191Z",
     "start_time": "2018-10-18T01:40:33.813490Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The binding affinity for butane and cucurbit[6]uril = -11.62 +/- 1.47 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "print(f\"The binding affinity for butane and cucurbit[6]uril = {binding_affinity:0.2f} +/- {sem:0.2f} kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the first tutorial, using implicit solvation, the binding affinity estimate was -9.00 +/- 5.53 kcal/mol. Now, the affinity estimate is more favorable, by almost 2 kcal/mol, and we've halved the uncertainty. As before, with more windows and more simulation time per window, this value would continue to be refined.\n",
    "\n",
    "The experimental value is $-RT \\ln (280*10^3 M)$ = -7.44 kcal/mol."
   ]
  }
 ],
 "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.7.8"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}