docs/tutorials/01-tutorial-cb6-but.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pAPRika tutorial 1 - An introduction to APR with Implicit 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). CBs are rigid, symmetric, cyclic host molecules with oxygen atoms around the portal edge of the cavity. We will run the simulation in implicit solvent, using the [Generalized-Born](https://en.wikipedia.org/wiki/Implicit_solvation#Generalized_Born) model, for speed and simplicity, using *AMBER*. This tutorial assumes familiarity with basic MD procedures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Setup\n",
    "\n",
    "The very first step in the calculation is creating a coordinate file for the bound host-guest complex (we usually use PDB format for this part because it works well with `tleap`). This does not have to perfectly match the bound state by any means (we will later minimize and equilibrate), but this should be a reasonable illustration of the bound complex. This file can be created by hand (in a program like Chimera, VMD, or PyMOL) or by docking the guest into the host cavity (with MOE, AutoDock, DOCK, ...).\n",
    "\n",
    "In this example, this file is called `cb6-but.pdb`, in the `complex/` directory, and this is what it looks like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![system](images/cb6-but.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the coordinate file, we need separate `mol2` files for the host and guest molecule that contain the partial atomic charges. For cyclic hosts like CB6, you can specify either a single residue for the entire molecule (this is what we do here) or you can provide coordinates and charges for a single monomer and have `tleap` build the structure (this is a little more tricky).\n",
    "\n",
    "The `mol2` files that we use here (`cb6.mol2` and `but.mol2`) were created by running `antechamber -fi pdb -fo mol2 -i <pdb> -o <mol2> -c bcc` (I added `pl 10` for the host, which reduces the number of paths that `antechamber` takes to traverse the host; see `antechamber` help for more information)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "if not os.path.isdir(\"complex\"):\n",
    "    os.makedirs(\"complex\")\n",
    "data = \"../../../paprika/data/cb6-but\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create AMBER coordinate (`.rst7` or `.inpcrd`) and parameter files (`.prmtop` or `.topo`) for the host-guest complex\n",
    "\n",
    "In this example, we will use GAFF parameters for both the host and guest. For the host, `parmchk2` has identified two parameters that are missing from GAFF and added the most similar ones into the supplementary `cb6.frcmod` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:01.109793Z",
     "start_time": "2018-10-18T01:05:00.780271Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.build.system import TLeap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:01.884510Z",
     "start_time": "2018-10-18T01:05:01.814547Z"
    }
   },
   "outputs": [],
   "source": [
    "system = TLeap()\n",
    "system.output_path = \"complex\"\n",
    "system.pbc_type = None\n",
    "system.neutralize = False\n",
    "\n",
    "system.template_lines = [\n",
    "    \"source leaprc.gaff\",\n",
    "    f\"loadamberparams {data}/cb6.frcmod\",\n",
    "    f\"CB6 = loadmol2 {data}/cb6.mol2\",\n",
    "    f\"BUT = loadmol2 {data}/but.mol2\",\n",
    "    f\"model = loadpdb {data}/cb6-but.pdb\",\n",
    "    \"check model\",\n",
    "    \"savepdb model vac.pdb\",\n",
    "    \"saveamberparm model vac.prmtop vac.rst7\"\n",
    "]\n",
    "system.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running `tleap`, it is always a good idea to check `leap.log`. `pAPRika` does some automated checking of the ouptut, but sometimes things slip through. (By default, `tleap` will append to `leap.log`.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare the complex for an APR calculation\n",
    "\n",
    "Now we are ready to prepare the complex for the attach-pull-release calculation. This involves:\n",
    "\n",
    "- Aligning the structure so the guest can be pulled along the $z$ axis, and\n",
    "- Adding dummy atoms that are used to orient the host and guest.\n",
    "\n",
    "To access the host-guest structure in Python, we use the ParmEd `Structure` class. So we start by loading the vacuum model that we just created. Then, we need to define two atoms on the guest that are placed along the $z$ axis. These should be heavy atoms on either end of the guest, the second atom leading the pulling.\n",
    "\n",
    "These same atoms will be used later for the restraints, so I will name them `G1` and `G2`, using AMBER selection syntax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Align the pulling axis\n",
    "<a id='guest'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:05.683159Z",
     "start_time": "2018-10-18T01:05:05.678056Z"
    }
   },
   "outputs": [],
   "source": [
    "import parmed as pmd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:06.344770Z",
     "start_time": "2018-10-18T01:05:06.277675Z"
    }
   },
   "outputs": [],
   "source": [
    "structure = pmd.load_file(\"complex/vac.prmtop\",\n",
    "                          \"complex/vac.rst7\",\n",
    "                          structure=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:06.834369Z",
     "start_time": "2018-10-18T01:05:06.830963Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.build import align"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:07.293389Z",
     "start_time": "2018-10-18T01:05:07.290781Z"
    }
   },
   "outputs": [],
   "source": [
    "G1 = \":BUT@C\"\n",
    "G2 = \":BUT@C3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:07.879912Z",
     "start_time": "2018-10-18T01:05:07.676530Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/peanut/pAPRika/paprika/build/align.py:255: RuntimeWarning: invalid value encountered in true_divide\n",
      "  x = np.cross(vector, ref_vector) / np.linalg.norm(np.cross(vector, ref_vector))\n"
     ]
    }
   ],
   "source": [
    "aligned_structure = align.zalign(structure, G1, G2)\n",
    "aligned_structure.save(\"complex/aligned.prmtop\", overwrite=True)\n",
    "aligned_structure.save(\"complex/aligned.rst7\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the origin is shown as a grey sphere, with the $z$ axis drawn as a blue arrow. The coordinates used for this example were already aligned, so Python warns that the cross product is zero, but this won't be the case in general.\n",
    "\n",
    "![aligned](images/aligned.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we add the dummy atoms. The dummy atoms will be fixed in place during the simulation and are used to orient the host and guest in the lab frame. The dummy atoms are placed along the $z$ axis, behind the host. The dummy atoms are used in distance, angle, and torsion restraints and therefore, the exact positioning of these atoms affects the value of those restraints. For a typical host-guest system, like the one here, we generally place the first dummy atom 6 Angstroms behind the origin, the second dummy atom 9 Angstroms behind the origin, and the third dummy atom 11.2 Angstroms behind the origin and offset about 2.2 Angstroms along the $y$ axis. After we add restraints, the positioning of the dummy atoms should be more clear.\n",
    "\n",
    "Note, these dummy atoms do not interact with the other atoms in the system, and therefore, there is no problem placing them near host atoms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add dummy atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:11.611676Z",
     "start_time": "2018-10-18T01:05:11.474577Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.build import dummy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:11.865268Z",
     "start_time": "2018-10-18T01:05:11.834672Z"
    }
   },
   "outputs": [],
   "source": [
    "structure = pmd.load_file(\"complex/aligned.prmtop\",\n",
    "                          \"complex/aligned.rst7\",\n",
    "                          structure=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:12.113747Z",
     "start_time": "2018-10-18T01:05:12.108580Z"
    }
   },
   "outputs": [],
   "source": [
    "structure = dummy.add_dummy(structure, residue_name=\"DM1\", z=-6.0)\n",
    "structure = dummy.add_dummy(structure, residue_name=\"DM2\", z=-9.0)\n",
    "structure = dummy.add_dummy(structure, residue_name=\"DM3\", z=-11.2, y=2.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:12.385810Z",
     "start_time": "2018-10-18T01:05:12.348801Z"
    }
   },
   "outputs": [],
   "source": [
    "structure.save(\"complex/aligned_with_dummy.prmtop\", overwrite=True)\n",
    "structure.save(\"complex/aligned_with_dummy.rst7\", overwrite=True)\n",
    "structure.save(\"complex/aligned_with_dummy.pdb\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![aligned_dummy](images/aligned_with_dummy.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we solvate the system in `tleap`, we will need `frcmod` files for the dummy atoms (otherwise `tleap` will use GAFF parameters and the dummy atoms will *not* be non-interacting). There is a convenient method in `paprika` to write a `frcmod` file that only contains a `MASS` section. For convenience, I am also going to write `mol2` files for each of the dummy atoms. This makes it easy to build up the system, piece-by-piece, if we have a separate `mol2` file for each component of the system: host, guest, dummy atoms.\n",
    "\n",
    "In principle, an APR calculate can be completed without dummy atoms, by simply lengthening the distance between the host and guest, but the addition of dummy atoms permits an easier way to think about dissociating the guest from the host. (Also, in the absence of dummy atoms, it is challenging to pull the guest straight out of the cavity without adding additional restraints.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:14.578465Z",
     "start_time": "2018-10-18T01:05:14.568671Z"
    }
   },
   "outputs": [],
   "source": [
    "dummy.write_dummy_frcmod(filepath=\"complex/dummy.frcmod\")\n",
    "dummy.write_dummy_mol2(residue_name=\"DM1\", filepath=\"complex/dm1.mol2\")\n",
    "dummy.write_dummy_mol2(residue_name=\"DM2\", filepath=\"complex/dm2.mol2\")\n",
    "dummy.write_dummy_mol2(residue_name=\"DM3\", filepath=\"complex/dm3.mol2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all the pieces are in place to build the system for an APR calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Put it all together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:18.731348Z",
     "start_time": "2018-10-18T01:05:18.644544Z"
    }
   },
   "outputs": [],
   "source": [
    "system = TLeap()\n",
    "system.output_path = \"complex\"\n",
    "system.pbc_type = None\n",
    "system.neutralize = False\n",
    "\n",
    "system.template_lines = [\n",
    "    \"source leaprc.gaff\",\n",
    "    f\"loadamberparams {data}/cb6.frcmod\",\n",
    "    \"loadamberparams dummy.frcmod\",\n",
    "    f\"CB6 = loadmol2 {data}/cb6.mol2\",\n",
    "    f\"BUT = loadmol2 {data}/but.mol2\",\n",
    "    \"DM1 = loadmol2 dm1.mol2\",\n",
    "    \"DM2 = loadmol2 dm2.mol2\",\n",
    "    \"DM3 = loadmol2 dm3.mol2\",\n",
    "    \"model = loadpdb aligned_with_dummy.pdb\",\n",
    "    \"check model\",\n",
    "    \"savepdb model cb6-but-dum.pdb\",\n",
    "    \"saveamberparm model cb6-but-dum.prmtop cb6-but-dum.rst7\"\n",
    "]\n",
    "system.build()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have AMBER coordinates and parameters for the `cb6-but` system with dummy atoms in the appropriate place and with the proper \"dummy\" parametesr."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![anchor](images/anchor-atoms-diagram2.png)"
   ]
  },
  {
   "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": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:22.941964Z",
     "start_time": "2018-10-18T01:05:22.935422Z"
    }
   },
   "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": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:05:23.130265Z",
     "start_time": "2018-10-18T01:05:23.126754Z"
    }
   },
   "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": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:29.480278Z",
     "start_time": "2018-10-18T01:06:29.477124Z"
    }
   },
   "outputs": [],
   "source": [
    "release_fractions = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Later, I will explain why there are no release windows in this calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:36.584915Z",
     "start_time": "2018-10-18T01:06:36.578561Z"
    }
   },
   "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",
    "<a id='host_static'></a>\n",
    "For the host-guest complex, we will apply four different types ot restraints in `pAPRika`:\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": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:38.198698Z",
     "start_time": "2018-10-18T01:06:38.193577Z"
    }
   },
   "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": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:38.868388Z",
     "start_time": "2018-10-18T01:06:38.863487Z"
    }
   },
   "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/Angstrom$^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/radian$^2$ for the angle force constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![static1](images/static-restraints-1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:40.535839Z",
     "start_time": "2018-10-18T01:06:40.529491Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika import restraints\n",
    "static_restraints = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:40.796728Z",
     "start_time": "2018-10-18T01:06:40.764151Z"
    }
   },
   "outputs": [],
   "source": [
    "structure = pmd.load_file(\"complex/cb6-but-dum.prmtop\",\n",
    "                          \"complex/cb6-but-dum.rst7\",\n",
    "                          structure = True\n",
    "                          )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:41.092360Z",
     "start_time": "2018-10-18T01:06:41.051265Z"
    }
   },
   "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": 23,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:41.474957Z",
     "start_time": "2018-10-18T01:06:41.434306Z"
    }
   },
   "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": 24,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:41.764205Z",
     "start_time": "2018-10-18T01:06:41.724550Z"
    }
   },
   "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": [
    "![static2](images/static-restraints-2.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:42.818533Z",
     "start_time": "2018-10-18T01:06:42.767475Z"
    }
   },
   "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": 26,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:43.042473Z",
     "start_time": "2018-10-18T01:06:43.002070Z"
    }
   },
   "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": 27,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:06:43.412015Z",
     "start_time": "2018-10-18T01:06:43.371189Z"
    }
   },
   "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": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Guest restraints\n",
    "<a id='guest'></a>\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 specify angle and torsion targets in degrees but the force constant using radians!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:25.544817Z",
     "start_time": "2018-10-18T01:07:25.540049Z"
    }
   },
   "outputs": [],
   "source": [
    "guest_restraints = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:25.774467Z",
     "start_time": "2018-10-18T01:07:25.739851Z"
    }
   },
   "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": 30,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:25.952457Z",
     "start_time": "2018-10-18T01:07:25.930000Z"
    }
   },
   "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": 31,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:26.096410Z",
     "start_time": "2018-10-18T01:07:26.074002Z"
    }
   },
   "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": 32,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:40.853027Z",
     "start_time": "2018-10-18T01:07:40.844270Z"
    }
   },
   "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",
    "    if not os.path.isdir(f\"windows/{window}\"):\n",
    "        os.makedirs(f\"windows/{window}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, we ask pAPRika to write the restraint information in each window\n",
    "\n",
    "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": 33,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:42.277892Z",
     "start_time": "2018-10-18T01:07:42.243610Z"
    }
   },
   "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": [
    "## Prepare host-guest system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Translation of the guest\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": 34,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:48.115289Z",
     "start_time": "2018-10-18T01:07:47.087840Z"
    }
   },
   "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\", overwrite=True)\n",
    "        structure.save(f\"windows/{window}/cb6-but-dum.rst7\", overwrite=True)"
   ]
  },
  {
   "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": [
    "## Simulation\n",
    "<a id='simulate'></a>\n",
    "Since we are going to run an implicit solvent simulation, we have everything ready to go. `paprika` has an `amber` module that can help setting default parameters for the simulation. There are some high level options that we set directly, like `simulation.path`, and then we call the function `config_gb_min()` to setup reasonable default simulation parameters for a minimization in the Generalized-Born ensemble. After that, we directly modify the simulation `cntrl` section to apply the positional restraints on the dummy atoms. \n",
    "\n",
    "For this part, you need to have the *AMBER* executables in your path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:49.559414Z",
     "start_time": "2018-10-18T01:07:49.554039Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika.simulate import AMBER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'm using the `logging` module to keep track of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:07:51.039563Z",
     "start_time": "2018-10-18T01:07:51.025208Z"
    }
   },
   "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:13:27.032804Z",
     "start_time": "2018-10-18T01:07:53.139754Z"
    },
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Energy Minimization\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "pycharm": {
     "name": "#%% \n"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-10-05 12:11:35 PM Running minimization in window a000...\n",
      "2020-10-05 12:11:39 PM Running minimization in window a001...\n",
      "2020-10-05 12:11:43 PM Running minimization in window a002...\n",
      "2020-10-05 12:11:47 PM Running minimization in window a003...\n",
      "2020-10-05 12:11:51 PM Running minimization in window a004...\n",
      "2020-10-05 12:11:55 PM Running minimization in window a005...\n",
      "2020-10-05 12:11:59 PM Running minimization in window a006...\n",
      "2020-10-05 12:12:04 PM Running minimization in window a007...\n",
      "2020-10-05 12:12:08 PM Running minimization in window a008...\n",
      "2020-10-05 12:12:12 PM Running minimization in window a009...\n",
      "2020-10-05 12:12:16 PM Running minimization in window a010...\n",
      "2020-10-05 12:12:20 PM Running minimization in window a011...\n",
      "2020-10-05 12:12:24 PM Running minimization in window a012...\n",
      "2020-10-05 12:12:28 PM Running minimization in window a013...\n",
      "2020-10-05 12:12:33 PM Running minimization in window p000...\n",
      "2020-10-05 12:12:37 PM Running minimization in window p001...\n",
      "2020-10-05 12:12:41 PM Running minimization in window p002...\n",
      "2020-10-05 12:12:45 PM Running minimization in window p003...\n",
      "2020-10-05 12:12:49 PM Running minimization in window p004...\n",
      "2020-10-05 12:12:53 PM Running minimization in window p005...\n",
      "2020-10-05 12:12:57 PM Running minimization in window p006...\n",
      "2020-10-05 12:13:01 PM Running minimization in window p007...\n",
      "2020-10-05 12:13:05 PM Running minimization in window p008...\n",
      "2020-10-05 12:13:09 PM Running minimization in window p009...\n",
      "2020-10-05 12:13:13 PM Running minimization in window p010...\n",
      "2020-10-05 12:13:17 PM Running minimization in window p011...\n",
      "2020-10-05 12:13:21 PM Running minimization in window p012...\n",
      "2020-10-05 12:13:25 PM Running minimization in window p013...\n",
      "2020-10-05 12:13:29 PM Running minimization in window p014...\n",
      "2020-10-05 12:13:33 PM Running minimization in window p015...\n",
      "2020-10-05 12:13:37 PM Running minimization in window p016...\n",
      "2020-10-05 12:13:41 PM Running minimization in window p017...\n"
     ]
    }
   ],
   "source": [
    "for window in window_list:\n",
    "    simulation = AMBER()\n",
    "    simulation.executable = \"sander\"\n",
    "\n",
    "    simulation.path = f\"windows/{window}/\"\n",
    "    simulation.prefix = \"minimize\"\n",
    "\n",
    "    simulation.topology = \"cb6-but-dum.prmtop\"\n",
    "    simulation.coordinates = \"cb6-but-dum.rst7\"\n",
    "    simulation.ref = \"cb6-but-dum.rst7\"\n",
    "    simulation.restraint_file = \"disang.rest\"\n",
    "\n",
    "    simulation.config_gb_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(overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Production\n",
    "For simplicity, I am going to skip equilibration and go straight to production!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:14:15.079063Z",
     "start_time": "2018-10-18T01:13:27.034729Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-10-05 12:13:45 PM Running production in window a000...\n",
      "2020-10-05 12:13:48 PM Running production in window a001...\n",
      "2020-10-05 12:13:52 PM Running production in window a002...\n",
      "2020-10-05 12:13:55 PM Running production in window a003...\n",
      "2020-10-05 12:13:59 PM Running production in window a004...\n",
      "2020-10-05 12:14:02 PM Running production in window a005...\n",
      "2020-10-05 12:14:06 PM Running production in window a006...\n",
      "2020-10-05 12:14:09 PM Running production in window a007...\n",
      "2020-10-05 12:14:13 PM Running production in window a008...\n",
      "2020-10-05 12:14:16 PM Running production in window a009...\n",
      "2020-10-05 12:14:19 PM Running production in window a010...\n",
      "2020-10-05 12:14:23 PM Running production in window a011...\n",
      "2020-10-05 12:14:26 PM Running production in window a012...\n",
      "2020-10-05 12:14:30 PM Running production in window a013...\n",
      "2020-10-05 12:14:33 PM Running production in window p000...\n",
      "2020-10-05 12:14:37 PM Running production in window p001...\n",
      "2020-10-05 12:14:40 PM Running production in window p002...\n",
      "2020-10-05 12:14:44 PM Running production in window p003...\n",
      "2020-10-05 12:14:47 PM Running production in window p004...\n",
      "2020-10-05 12:14:51 PM Running production in window p005...\n",
      "2020-10-05 12:14:54 PM Running production in window p006...\n",
      "2020-10-05 12:14:57 PM Running production in window p007...\n",
      "2020-10-05 12:15:01 PM Running production in window p008...\n",
      "2020-10-05 12:15:04 PM Running production in window p009...\n",
      "2020-10-05 12:15:07 PM Running production in window p010...\n",
      "2020-10-05 12:15:11 PM Running production in window p011...\n",
      "2020-10-05 12:15:14 PM Running production in window p012...\n",
      "2020-10-05 12:15:17 PM Running production in window p013...\n",
      "2020-10-05 12:15:21 PM Running production in window p014...\n",
      "2020-10-05 12:15:24 PM Running production in window p015...\n",
      "2020-10-05 12:15:27 PM Running production in window p016...\n",
      "2020-10-05 12:15:31 PM Running production in window p017...\n"
     ]
    }
   ],
   "source": [
    "for window in window_list:\n",
    "    simulation = AMBER()\n",
    "    simulation.executable = \"sander\"\n",
    "    \n",
    "    simulation.path = f\"windows/{window}/\"\n",
    "    simulation.prefix = \"production\"\n",
    "\n",
    "    simulation.topology = \"cb6-but-dum.prmtop\"\n",
    "    simulation.coordinates = \"minimize.rst7\"\n",
    "    simulation.ref = \"cb6-but-dum.rst7\"\n",
    "    simulation.restraint_file = \"disang.rest\"\n",
    "\n",
    "    simulation.config_gb_md()\n",
    "    simulation.cntrl[\"ntr\"] = 1\n",
    "    simulation.cntrl[\"restraint_wt\"] = 50.0\n",
    "    simulation.cntrl[\"restraintmask\"] = \"'@DUM'\"\n",
    "    \n",
    "    logging.info(f\"Running production in window {window}...\")\n",
    "    simulation.run(overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "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": 39,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:14:19.273341Z",
     "start_time": "2018-10-18T01:14:19.268789Z"
    }
   },
   "outputs": [],
   "source": [
    "from paprika import analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:14:29.940490Z",
     "start_time": "2018-10-18T01:14:20.084861Z"
    }
   },
   "outputs": [],
   "source": [
    "free_energy = analysis.fe_calc()\n",
    "free_energy.topology = \"cb6-but-dum.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": "markdown",
   "metadata": {},
   "source": [
    "But what about release? The guest's rotational and translational degrees of freedom are still restrained releative to the frame of reference of the host. The work to release these restraints is the difference between the chemical potential of this state, and the chemical potentials of the separate host and guest at standard concentration. This can be calculated analytically without doing additional simulation (see [equation 8](https://pubs.acs.org/doi/10.1021/acs.jctc.5b00405)), using the function `compute_ref_state_work`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:17:22.534871Z",
     "start_time": "2018-10-18T01:17:22.450355Z"
    }
   },
   "outputs": [],
   "source": [
    "free_energy.compute_ref_state_work([\n",
    "    guest_restraints[0], guest_restraints[1], None, \n",
    "    None, guest_restraints[2], None\n",
    "])\n",
    "\n",
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-10-18T01:17:23.400457Z",
     "start_time": "2018-10-18T01:17:23.394759Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The binding affinity for butane and cucurbit[6]uril = -2.98 +/- 5.82 kcal/mol\n"
     ]
    }
   ],
   "source": [
    "print(f\"The binding affinity for butane and cucurbit[6]uril = {binding_affinity.magnitude:0.2f} +/- {sem.magnitude:0.2f} kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='results'></a>\n",
    "There is a large uncertainty associated with this calculation because we only simulated for a very short amount of time in each window and we used a large amount of spacing between each window in the pull phase, but the uncertainty will go down with more time.\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
}