3_committor_analysis_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Committor Analysis Tutorial\n",
    "\n",
    "As with the TPS analysis tutorial, this is based on results that have already been calculated. This was a simulation of 1000 points selected from a TPS simulation, with 10 committor shots performed on each one -- so a total of 10000 MD simulations.\n",
    "\n",
    "The committor is closely related to the reaction coordinate: One way to estimate the reaction coordinate is to follow the increasing isocommittors (surfaces with the same value of the committor). The isocommittor with value of 0.5 is particularly interesting: the corresponds to the transition state, where a trajectory started with random velocities has equal probability of ending in reactants or products.\n",
    "\n",
    "In practice, a full committor analysis is usually very difficult. This one is only partly converged, even with the many, many individual MD simulations included. To reduce the file size, we only saved the first and last snapshot from each committor shot. You can get the file from: https://figshare.com/s/01302bc7a39ec7648ea1\n",
    "\n",
    "Most of this tutorial is focused on how to use [`OpenPathSampling`](http://openpathsampling.org), however, a couple other libraries play an important role:\n",
    "\n",
    "* [`matplotlib`](http://matplotlib.org/), a library for creating plots\n",
    "* [`pandas`](http://pandas.pydata.org/), a library for data analysis\n",
    "* [`nglview`](https://github.com/arose/nglview), which we'll again use to visualize the molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import openpathsampling as paths\n",
    "import matplotlib\n",
    "matplotlib.rcParams.update({'font.size': 18})\n",
    "matplotlib.rcParams.update({'figure.figsize': (8.8, 6.6)})\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "pd.options.display.max_rows = 10\n",
    "\n",
    "import nglview as nv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we open the file, and load the states and CVs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "simulation_storage = paths.AnalysisStorage(\"committor_small.nc\")\n",
    "C_7eq = simulation_storage.volumes['C_7eq']\n",
    "alpha_R = simulation_storage.volumes['alpha_R']\n",
    "phi = simulation_storage.cvs['phi']\n",
    "psi = simulation_storage.cvs['psi']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"File size: {0} for {1} steps, {2} snapshots\".format(\n",
    "    simulation_storage.file_size_str,\n",
    "    len(simulation_storage.steps),\n",
    "    len(simulation_storage.snapshots)\n",
    "))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The initial point for a committor analysis is seen by OPS as a shooting point, so the object to analyze it is `ShootingPointAnalysis`. This is, again, a very time-consuming process (you're analyzing the results of 10000 MD simulations!) It takes ~20 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyzer = paths.ShootingPointAnalysis(steps=simulation_storage.steps,\n",
    "                                       states=[C_7eq, alpha_R])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`pandas` is a convenient and widely-used library for a lot of spreadsheet-like data analysis. In addition, it provides a very nice way to present results in Jupyter notebooks, so OpenPathSampling provides the option to return the committor results as a `pandas.DataFrame` object. The first column shows the index number of the initial shooting snapshot, and the other columns show how many landed in each state (where `NaN`, \"not-a-number\", indicates that 0 shots landed there)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analyzer.to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using a `label_function`, we can convert the index number of the snapshot to a collective variable that is representative of that snapshot. In this particular case, we'll use $\\phi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "phi_hash = lambda x : float(phi(x))\n",
    "analyzer.to_pandas(label_function=phi_hash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's histogram the committor analysis according to the value of `phi`. This requires using the `phi_hash`, and will give a one-dimensional histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist1D, phi_bins = analyzer.committor_histogram(phi_hash, alpha_R, bins=10)\n",
    "bin_widths = [phi_bins[i+1]-phi_bins[i] for i in range(len(phi_bins)-1)]\n",
    "plt.bar(x=phi_bins[:-1], height=hist1D, width=bin_widths)\n",
    "plt.xlabel(\"$\\phi$\")\n",
    "plt.ylabel(\"$p_B$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Is $\\phi$ a good representation of the reaction coordinate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you'll do exactly the same for the other main CV we use in this: $\\psi$. Follow the procedure used above to make a one-dimensional histogram according to a specific CV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR TURN: Do the same for the psi variable.\n",
    "# Remember to start by creating a `psi_hash`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Compare this with the $\\phi$ histogram. Which is a better representation of the reaction coordinate: $\\phi$ or $\\psi$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the one-dimensional histogram enough information? Here we'll make a two-dimensional histogram including both $\\phi$ and $\\psi$. The color will tell us the value of the committor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "ramachandran_hash = lambda x : (float(phi(x)), float(psi(x)))\n",
    "hist2D, bins_phi, bins_psi = analyzer.committor_histogram(ramachandran_hash, alpha_R, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolor(bins_phi, bins_psi, hist2D.T, cmap=\"winter\")\n",
    "plt.clim(0.0, 1.0)\n",
    "plt.xlabel(\"$\\phi$\")\n",
    "plt.ylabel(\"$\\psi$\")\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the most interesting part of the committor is the 0.5 isocommittor, we'll focus on that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# identify the committor 0.5 frames; histogram them by psi\n",
    "# we use a powerful Python trick called a \"list comprehension\" to make `isocommittor_0x5`\n",
    "# list comprehensions are a bit hard to understand at first, but very useful once you're used to them\n",
    "committor = analyzer.committor(alpha_R)\n",
    "isocommittor_0x5 = [s for s in committor if 0.45 < committor[s] < 0.55]\n",
    "\n",
    "plt.hist(psi(isocommittor_0x5))\n",
    "plt.xlabel(\"$\\psi$\")\n",
    "plt.ylabel(\"50% Isocommittor Frequency\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice that this is *not* very sharply peaked in $\\psi$. This means that $\\psi$ alone probably doesn't define the transition state.\n",
    "\n",
    "If $\\psi$ did define the transition state, we'd probably expect the transition state to be around $\\psi = 0.6$ (\\~35 degrees) to $\\psi = 1.0$ (\\~60 degrees). Now, instead of keeping the committor probability fixed and histogramming the values of $\\psi$, let's see what the values of the committor are for a fixed value of $\\psi = \\psi^*$. If $\\psi^*$ defined the transition state, then the committor would always be 0.5 here.\n",
    "\n",
    "Previous, we tried to show that frames with a committor around 0.5 had a certain value of $\\psi$. Now we're trying to see if that value of $\\psi$ always implies that frames have a committor around 0.5. We need to show both directions before we can claim that the transition state is defined by that value of $\\psi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a range around a central psi_star value\n",
    "psi_star_min = 0.6\n",
    "psi_star_max = 0.7\n",
    "\n",
    "psi_star_snapshots = [s for s in committor \n",
    "                      if psi_star_min <= psi(s) <= psi_star_max]\n",
    "plt.hist([committor[s] for s in psi_star_snapshots]);\n",
    "plt.xlabel(\"$p_B$\")\n",
    "plt.ylabel(\"$P(p_B)$ (unnormalized)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try a few different values of `psi_star_min` and `psi_star_max`. \n",
    "\n",
    "* Do you ever get a distribution that is peaked around $p_B = 0.5$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you should select a few frames from the 50% isocommittor surface (i.e., the transition state ensemble) and visualize them using `nglview`. The variable `isocommittor_0x5` contains a list of the appropriate snapshots, so you can pick whichever ones you want. Let's see how many such frames we found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(isocommittor_0x5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Repeat the following several times, changing the number `2` in `isocommittor_0x5[2]` to other choices, and compare the conformations you see. (Don't expect to be able to identify the \"hidden\" variable: The conclusion from [Bolhuis et al.](http://dx.doi.org/10.1073/pnas.100127697) was that the water itself plays an important role in this reaction, and even they couldn't identify exactly what the water-related reaction coordinate should be!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR TURN: modify the snapshot number; compare the conformations\n",
    "md_traj_A = paths.Trajectory([isocommittor_0x5[2]]).to_mdtraj().image_molecules()\n",
    "view_A = nv.show_mdtraj(md_traj_A)\n",
    "view_A.clear()\n",
    "view_A.add_ball_and_stick(\"ACE ALA NME\")\n",
    "view_A.center()\n",
    "view_A"
   ]
  }
 ],
 "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.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}