Alanine dipeptide committor stuff

0_README.md

Alanine dipeptide committor stuff

This directory is part of a collaboration with Max Welling and his colleagues Rianne van den Berg and Thomas Kipf. We're playing with ways to use deep learning neural networks in the context of chemistry.

This is mainly for me to create some useful data for them, and for me to explain some of the tools that I think they will find useful. Some of this may get re-purposed into OpenPathSampling examples.

Here I am taking trajectories from a previously performed (flex-length) TPS calculation on the C7eq to alpha_R states of alanine dipeptide. I select a subset of the snapshots, and run a committor analysis on them. In the end, I give the experiment-by-experiment results, where an "experiment" is the mapping of an initial MD snapshot to the state it ended up in.

• select_snapshots.ipynb: Select the initial set from snapshots. Input requires alanine_dipeptide_tps.nc (calculated for the OPS paper). Output is snapshots.nc.
• commttor_simulation.ipynb: Committor analysis for the selected snapshots. Output is committor_simulation.nc.
• committor_analysis.ipynb: Analysis of the committor simulation. Output is in committor_results.nc. This output includes the mapping of point in phase space to final state (binary value) that we discussed.
• analysis_help.ipynb: How to load the output from the former and get relevant information from it.
• mdtraj_tricks.ipynb: Introduction to a few potentially relevant aspects of MDTraj. Examples of how to create descriptors.

(Rianne and Thomas, the parts that are most relevant to you are mdtraj_tricks.ipynb and analysis_help.ipynb. The rest provides a record of how I'm creating these files.)

---

Information on various programs that interact with each other, to give context for Rianne and Thomas:

• OpenMM: Molecular dynamics code. Unless you're actually doing simulations, you won't need to worry about this.
• MDTraj: Trajectory analysis code, written by some of the same people as OpenMM, but useful for many other packages. You should use this for descriptors.
• OpenMMTools: Useful shortcuts for working with OpenMM. Again, you won't need this unless you actually do simulations.
• OpenPathSampling (OPS): What we use for our path sampling simulations. The files I create are readable by OPS, so you'll be working with some OPS objects.

If you can install OpenPathSampling using conda (you'll need to add the omnia channel), you should get all of these automatically. There are a few things I use that still aren't part of an OPS release, so you might want prefer to use a git-based local install of OPS in practice.

1_select_snapshots.ipynb

2_committor_simulation.ipynb

3_committor_analysis.ipynb

4_analysis_help.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis help\n",
    "\n",
    "This covers stuff that you will need to know in order to use the `committor_results.nc` file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import openpathsampling as paths\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "pd.options.display.max_rows = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "storage = paths.Storage(\"committor_results.nc\", \"r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "phi = storage.cvs['phi']\n",
    "psi = storage.cvs['psi']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 49.6 s, sys: 239 ms, total: 49.8 s\n",
      "Wall time: 51.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "C_7eq = storage.volumes['C_7eq']\n",
    "alpha_R = storage.volumes['alpha_R']\n",
    "experiments = storage.tag['experiments']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `experiments` object is a list of tuples `(snapshot, final_state)`. Each `snapshot` is an OPS snapshot object (a point in phase space), and the `final_state` is either the `C_7eq` object or the `alpha_R` object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Directly obtaining a committor analysis\n",
    "\n",
    "As it happens, `experiments` is in precisely the correct format to be used in one of the approaches to constructing a committor analysis.\n",
    "\n",
    "This section requires OpenPathSampling 0.9.1 or later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 44 s, sys: 143 ms, total: 44.2 s\n",
      "Wall time: 49.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "committor_analyzer = paths.ShootingPointAnalysis.from_individual_runs(experiments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before going further, let's talk a little bit about the implementation of the `ShootingPointAnalysis` object. The main thing to understand is that the purpose of that object is to histogram according to configuration. The first snapshot encountered is kept as a representative of that configuration.\n",
    "\n",
    "So whereas there are 10000 snapshots in `experiments` (containing the full data, including velocities), there are only 1000 entries in the `committor_analyzer` (because, in this data set, I ran 1000 snapshots with 10 shots each.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Per-configuration results\n",
    "\n",
    "The `.to_pandas()` function creates a pandas table with configurations as the index, the final states as columns, and the number of times that configuration led to that final state as entries. With no argument, `to_pandas()` using the an integer for each configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>C_7eq</th>\n",
       "      <th>alpha_R</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>10.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>995</th>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>996</th>\n",
       "      <td>4.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>997</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>998</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>999</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     C_7eq  alpha_R\n",
       "0      8.0      2.0\n",
       "1      5.0      5.0\n",
       "2      9.0      1.0\n",
       "3      9.0      1.0\n",
       "4     10.0      NaN\n",
       "..     ...      ...\n",
       "995    3.0      7.0\n",
       "996    4.0      6.0\n",
       "997    9.0      1.0\n",
       "998    8.0      2.0\n",
       "999    8.0      2.0\n",
       "\n",
       "[1000 rows x 2 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "committor_analyzer.to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also pass it a function that takes a snapshot and returns a (hashable) value. That value will be used for the index. These collective variables return numpy arrays, so we need to cast the 1D array to a `float`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>C_7eq</th>\n",
       "      <th>alpha_R</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.824873</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.576733</th>\n",
       "      <td>5.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.521896</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.444523</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.977096</th>\n",
       "      <td>10.0</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.419521</th>\n",
       "      <td>3.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.353020</th>\n",
       "      <td>4.0</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.492791</th>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.054834</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.254955</th>\n",
       "      <td>8.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1000 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "          C_7eq  alpha_R\n",
       "0.824873    8.0      2.0\n",
       "0.576733    5.0      5.0\n",
       "0.521896    9.0      1.0\n",
       "0.444523    9.0      1.0\n",
       "0.977096   10.0      NaN\n",
       "...         ...      ...\n",
       "0.419521    3.0      7.0\n",
       "0.353020    4.0      6.0\n",
       "1.492791    9.0      1.0\n",
       "1.054834    8.0      2.0\n",
       "1.254955    8.0      2.0\n",
       "\n",
       "[1000 rows x 2 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi_hash = lambda x : float(psi(x))\n",
    "committor_analyzer.to_pandas(label_function=psi_hash)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also directly obtain the committor as a dictionary of (representative) snapshot to committor value. The committor here is defines as the probability of ending in a given state, so you must give the state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "committor = committor_analyzer.committor(alpha_R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{<openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f116b90>: 0.0,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f192f50>: 0.4,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f31a990>: 0.0,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f34ef10>: 0.3,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f5c13d0>: 0.3,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10f9e1950>: 0.0,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x10fd170d0>: 0.3,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x1100f4c50>: 0.1,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x1101f7410>: 0.2,\n",
       " <openpathsampling.engines.openmm.snapshot.Snapshot at 0x1103bd390>: 0.0}"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# show the first 10 values\n",
    "{k: committor[k] for k in committor.keys()[:10]}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Committor histogram in 1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "hist1D, bins = committor_analyzer.committor_histogram(psi_hash, alpha_R, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAFnCAYAAADzOqBQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAGSNJREFUeJzt3X+U5WddH/D3ZxNkTSpp7RRQccUU3S79AWTAklIKbQrR\nohxAEAeGKqn0pKQ9urUirVJaTsUDFlLoMZValdDAVHoO9aSUQwBBQMKPkiEgdFg5NnRDXKJDzFoI\na5B9+se9W4bJzM7cZ++de+/s63XOnJP57v0+9/Nk7tx53+/3+VGttQAA9Dgw7QIAgPklSAAA3QQJ\nAKCbIAEAdBMkAIBuggQA0E2QAAC6CRIAQDdBAgDoJkgAAN0ECQCg20wEiar6/qr6dFUdq6p/MO16\nAIDdqWlv2lVVFyT5X0mekOT/JrklyWNba3dPtTAAYEezcEXie5J8srX2+dbal5K8LcmVU64JANiF\nWQgS35rkjg3f/36Sb5tSLQDACM4pSFTV46vqxqq6o6pOV9VTt3jMNVV1W1V9uao+VFWP2fyQLZqe\n7v0WAGBXzvWKxMVJbk1yTbb4419Vz07yqiQvTfKoJB9PclNVLWx42B1JHrLh+29LcuIc6wIA9sDY\nBltW1ekkT2ut3bjh2IeSfLi19uPD7yvJ7Ule21p75fDYmcGWT8xgsOX/TPI3Wmt/tM3z/PkMxlB8\nNsmpsRQPAOeHg0kemuSm1toXxtHgheNoZCtVdb8ki0lefuZYa61V1buSXL7h2Fer6ieT/FYGtzle\nsV2IGLoyyRsnUjQAnB+em+RN42hoYkEiyUKSC5Lcuen4nUkObzzQWntrkrfust3PJskNN9yQI0eO\nnGOJ03f06NFce+210y5jbPRndu2nviT6M8v2U1+S/dWftbW1LC8vJ8O/peMwySCxncq5DaY8lSRH\njhzJZZddNp6KpuiSSy7ZF/04Q39m137qS6I/s2w/9SXZf/0ZGtvQgElO/1xP8tUkD9p0/IG571UK\nAGAOTSxItNa+ksEqlVecOTYcbHlFkpsn9bwAwN45p1sbVXVxkofla2tBXFpVj0hyV2vt9iSvTnJ9\nVd2S5CNJjia5KMnrz+V5AYDZcK5jJB6d5D0ZjHloGawZkSTXJ7mqtfbm4ZoRL8vgFsetSa5srf3h\nOT7vvrG0tDTtEsZKf2bXfupLoj+zbD/1Jdl//Rm3qW/aNaqquizJLbfccst+HPwCABOzurqaxcXF\nJFlsra2Oo81pzNoYi6NHj+aSSy7J0tKStAgAZ7GyspKVlZWcPHly7G27IgEA5wlXJDY4ceJEVlfH\n8v/gPhYWFnLo0KGJtA0A+8ncBolnPOOZuffeyWy1cfDgRTl2bE2YAIAdzG2QGISIG5KMe5nstZw6\ntZz19XVBAgB2MLdBYuBIEuMkAGBaJrlENgCwz835FYmjSa5OYvonAGxnktM/5zxIXBu3NgDg7M6s\nubRh+ufYuLUBAHQTJACAboIEANBNkAAAugkSAEC3OZ+1MTlra2tjb9MeHgDsN4LEfZxIciDLy8tj\nb9keHgDsN3MeJCaxINXdSU5n/Pt42MMDgOmwINW2JrkglX08ANgfJrkg1ZwHiflj7AUA+4kgsWeM\nvQBg/xEk9oyxFwDsP4LEnjP2AoD9w4JUAEA3QQIA6CZIAADdjJHYJyYxrTQxtRSAsxMk5t7kppUm\nppYCcHZzHiQmsUT2vJnUtNLE1FKA/cES2dua5BLZ88a0UgC2Nsklsg22BAC6CRIAQDdBAgDoNudj\nJNgLdiwFYDuCBGdhx1IAzk6Q4CzsWArA2QkS7MJkppa6ZQIw/wQJpsAtE4D9QpBgCtwyAdgvBAmm\nyGqcAPNuzoOEvTYAYCf22tiWvTYAYCf22gAAZpIgAQB0EyQAgG6CBADQTZAAALoJEgBAN0ECAOgm\nSAAA3QQJAKDbnK9sCfc1ie3JE1uUA2xFkGAfmdz25IktygG2Ikiwj0xqe/LEFuUAWxMk2IdsTw6w\nV+Y8SNhGHAB2YhvxbdlGnL01iYGcBnECkzbJbcTnPEjAXpncQE6DOIF5JkjArkxqIKdBnMB8EyRg\nJAZyAmxkZUsAoJsgAQB0EyQAgG6CBADQTZAAALqZtQEzwI6lwLwSJGCq7FgKzDdBAqZq8juWvv/9\n78+RI+Nt25UO4AxBAmbCJBa6sqw3MHmCBOxblvUGJk+QgH3Pst7A5Jj+CQB0EyQAgG6CBADQbc7H\nSBxNcnWSpWkXAgAza2VlJSsrKzl58uTY257zIHFtDCIDgLNbWlrK0tJSVldXs7i4ONa23doAALoJ\nEgBAN0ECAOgmSAAA3QQJAKCbIAEAdBMkAIBuggQA0E2QAAC6CRIAQDdBAgDoNud7bQDTsra2NvY2\nFxYWcujQobG3C0yOIAGM6ESSA1leXh57ywcPXpRjx9aECZgjggQworuTnE5yQ5IjY2x3LadOLWd9\nfV2QgDkiSACdjiS5bNpFAFNmsCUA0E2QAAC6CRIAQDdjJICZMolppYmppTApggQwIyY3rTQxtRQm\nRZAAZsSkppUmppbC5AgSwIwxrRTmicGWAEA3QQIA6CZIAADdBAkAoNucD7Y8muTqJEvTLgQAZtbK\nykpWVlZy8uTJsbc950Hi2hjdDQBnt7S0lKWlpayurmZxcXGsbbu1AQB0EyQAgG6CBADQTZAAALoJ\nEgBAN0ECAOgmSAAA3QQJAKDbnC9IBTBdx48fz/r6+kTaXlhYyKFDhybSNoyLIAHQ6fjx4zl8+EhO\nnbpnIu0fPHhRjh1bEyaYaYIEQKf19fVhiLghyZExt76WU6eWs76+Lkgw0wQJgHN2JPb94XxlsCUA\n0M0VCeC8sba2NtPtwTwSJIDzwIkkB7K8vDztQmDfESSA88DdSU5n/IMi35bkJWNsD+aPIAGcR8Y9\nKNKtDTDYEgDoJkgAAN0ECQCgmyABAHQTJACAboIEANBNkAAAugkSAEA3QQIA6CZIAADdBAkAoJsg\nAQB0EyQAgG6CBADQTZAAALoJEgBAN0ECAOgmSAAA3QQJAKCbIAEAdJuJIFFVb6mqu6rqzdOuBQDY\nvZkIEklek+R50y4CABjNTASJ1tp7k3xx2nUAAKOZiSABAMynkYNEVT2+qm6sqjuq6nRVPXWLx1xT\nVbdV1Zer6kNV9ZjxlAsAzJKeKxIXJ7k1yTVJ2uZ/rKpnJ3lVkpcmeVSSjye5qaoWNjzmhVX1sapa\nrar7d1UOAEzdhaOe0Fp7e5K3J0lV1RYPOZrkda21Nwwfc3WSpyS5Kskrh21cl+S6TefV8AsAmBMj\nB4mzqar7JVlM8vIzx1prrareleTys5z3ziR/LcnFVXU8ybNaax/e+RmPJrlkw/dLwy8AOL+trKxk\nZWXl646dPHly7M8z1iCRZCHJBUnu3HT8ziSHtzuptfakvqe7NsllfacCwD62tLSUpaWv/3C9urqa\nxcXFsT7PXs3aqGwxngIAmG/jDhLrSb6a5EGbjj8w971KAQDMubHe2mitfaWqbklyRZIbk/8/IPOK\nJK8d53MBMFuOHz+e9fX1ibS9sLCQQ4cOTaRtzs3IQaKqLk7ysHxthsWlVfWIJHe11m5P8uok1w8D\nxUcyGBF5UZLXj6ViAGbO8ePHc/jwkZw6dc9E2j948KIcO7YmTMygnisSj07yngzGPLQM1oxIkuuT\nXNVae/NwzYiXZXCL49YkV7bW/nAM9QIwg9bX14ch4oYkR8bc+lpOnVrO+vq6IDGDetaReG92GFux\nzToRAOx7R2I23fll3NM/99jRJFfH2hEAsL0za0rMwzoSe8w6EgCwkzNrSszzOhIAwD4051ckAPa3\ntbW1sbdpKiXjJEgAzKQTSQ5keXl57C2bSsk4CRIAM+nuJKcz/umUplIyXoIEwEwznZLZNudBwvRP\nANiJ6Z/bMv0ToMe4B3FOYlAo4zPJ6Z9zHiQAGM3kBnFyfhIkAM4rkxrE+bYkLxlje8wLQQLgvDTu\nQZxubZyvrGwJAHQTJACAboIEANBNkAAAus35YEsLUgHATixItS0LUgHATia5IJVbGwBAN0ECAOgm\nSAAA3QQJAKCbIAEAdBMkAIBuggQA0E2QAAC6CRIAQLc5X9nSEtkA9Dt+/HjW19fH3u7CwkIOHTo0\n9nZ7WSJ7W5bIBqDP8ePHc/jwkZw6dc/Y2z548KIcO7Y2M2Fikktkz3mQAIA+6+vrwxBxQ5IjY2x5\nLadOLWd9fX1mgsQkCRIAnOeOxNXtfgZbAgDdBAkAoJtbGwDMhbW1tZlu73wlSAAw404kOZDl5eVp\nF8IWBAkAZtzdSU5n/LMr3pbkJWNs7/wkSAAwJ8Y9u8KtjXEw2BIA6CZIAADd5vzWhr02AGAn9trY\nlr02AGAnk9xrw60NAKCbIAEAdBMkAIBuggQA0E2QAAC6CRIAQDdBAgDoJkgAAN0ECQCgmyABAHQT\nJACAboIEANBNkAAAus357p+2EQeAndhGfFu2EQeAndhGHACYSYIEANBNkAAAugkSAEA3QQIA6CZI\nAADdBAkAoJsgAQB0EyQAgG6CBADQTZAAALoJEgBAN0ECAOgmSAAA3QQJAKCbIAEAdBMkAIBuF067\ngHNzNMnVSZamXQgAzKyVlZWsrKzk5MmTY297zoPEtUkum3YRADDTlpaWsrS0lNXV1SwuLo61bbc2\nAIBuggQA0E2QAAC6CRIAQDdBAgDoJkgAAN0ECQCgmyABAHQTJACAboIEANBNkAAAugkSAEA3QQIA\n6CZIAADdBAkAoNuF0y4AAPajtbW1ibS7sLCQQ4cOTaTtHoIEAIzViSQHsry8PJHWDx68KMeOrc1M\nmBAkAGCs7k5yOskNSY6Mue21nDq1nPX1dUECAPa3I0kum3YRE2ewJQDQTZAAALoJEgBAN0ECAOgm\nSAAA3QQJAKCbIAEAdJvzdSSOJrk6ydK0CwGAmbWyspKVlZWcPHly7G3PeZC4NufDYh8AcC6Wlpay\ntLSU1dXVLC4ujrVttzYAgG6CBADQTZAAALoJEgBAN0ECAOgmSAAA3QQJAKCbIAEAdBMkAIBuggQA\n0E2QAAC6CRIAQDdBAgDoJkgAAN0ECQCgmyABAHQTJACAboIEANBNkAAAugkSAEA3QQIA6CZIAADd\nBAkAoJsgAQB0EyQAgG6CBADQTZAAALoJEgBAN0ECAOgmSAAA3QQJAKCbIAEAdBMkAIBuggQA0E2Q\nAAC6CRIAQDdBAgDoNvUgUVUPqar3VNWnqurWqnrmtGsCAHZn6kEiyZ8m+fHW2l9OcmWSf1dV3zjl\nmvbQyrQLGDP9mV37qS+J/syy/dSXZP/1Z7ymHiRaa59vrX1i+N93JllP8s3TrWov7bcXqP7Mrv3U\nl0R/Ztl+6kuy//ozXlMPEhtV1WKSA621O6ZdCwCws5GDRFU9vqpurKo7qup0VT11i8dcU1W3VdWX\nq+pDVfWYXbT7zUmuT/KCUWsCAKaj54rExUluTXJNkrb5H6vq2UleleSlSR6V5ONJbqqqhQ2PeWFV\nfayqVqvq/lX1DUn+W5KXt9Y+3FETADAFF456Qmvt7UneniRVVVs85GiS17XW3jB8zNVJnpLkqiSv\nHLZxXZLrzpxQVStJfrO19qZdlHDwa/+5Nmr5u3DbhNrert2TSVYn0O449LS9m/7s9f/jc3GmP/NU\n83Ztn+trbbt2x2XUdkfpz6zUfDab+zMPNW/X7rhea1u1PU67bXfU/kzy93rQ5tpaX9sbzjt4tseN\nolq7z0WF3Z9cdTrJ01prNw6/v1+Se5L84Jljw+OvT3JJa+3pW7TxuCTvTfKJJJXBVY7ntdY+tc1z\nPifJG7uLBgCeu8sP7zsa+YrEDhaSXJDkzk3H70xyeKsTWmsfGLGOm5I8N8lnk5wavUQAOG8dTPLQ\nDP6WjsW4g8R2zlxpOGettS8kGUuKAoDz0M3jbGzc0z/Xk3w1yYM2HX9g7nuVAgCYc2MNEq21ryS5\nJckVZ44NB2RekTEnIABg+ka+tVFVFyd5WAa3K5Lk0qp6RJK7Wmu3J3l1kuur6pYkH8lgFsdFSV4/\nlooBgJkx8qyNqnpCkvfkvmMerm+tXTV8zAuTvCiDWxy3JvknrbWPnnu5AMAsGfnWRmvtva21A621\nCzZ9XbXhMde11h7aWvvG1trlo4aIUVfGrKpnVdXa8PEfr6rvG7VfkzJKX6rqx6rqfVV11/DrnbtZ\nFXQv9axaOjzvh4crob5l0jWOouO1dklV/WJV/f7wnE9X1ffuVb1n09GXnxjWf09VHa+qV1fV/feq\n3rPZzQq6W5zzxKq6papOVdXvVtWP7EWtOxm1L1X19Kp6R1X9QVWdrKqbq+rJe1XvTnp+NhvOfVxV\nfaWqxrnIxDnpfK19Q1X9XFV9dvh6+99V9aN7UO5OdfX05bk12Hn7S8P3tV8ZrjS9azO110ayu5Ux\nNz3+8gxmcfxykkcm+Y0kv1FVD9+birc3al+SPCGDvjwxyWOT3J7kHVX1LZOvdmcd/Tlz3nck+YUk\n75t4kSPoeK3dL8m7khxK8owMpjS/IMnU94bp6Mtzkvz88PF/KYMF456d5Of2pOCdnXUF3c2q6qFJ\n3prkN5M8IslrkvynqnrS5ErctZH6kuRvJXlHku9LclkGV4D/+/AW8iwYtT9Jkqp6QAbbILxrQnX1\n6unPf03yt5M8P8l3J1lKcmwi1Y1m1N+bx2XwM/nlJA9P8swk35PkP470rK21mfpK8qEkr9nwfSX5\nXJIXbfP4/5Lkxk3HPpjkunnryxbnH8hgSbXlafeltz/DPrw/g1+4X0vylmn34xxea1cn+UySC6Zd\n+xj68u+TvHPTsX+b5H3T7ssWtZ5O8tQdHvOKJJ/YdGwlydumXf+ofdnmvE8m+dlp138u/Rn+PP51\nBuF1ddq19/YnyfcmuSvJn512vWPoy08m+cymY/84yfFRnmumrkgMP/EtZvCpIknSBj17V5LLtznt\n8tw34d50lsfvic6+bHZxkvtl8KKdqnPoz0uT/EFr7dcmW+FoOvvzAxmG1Kr6fFX9TlX986qa6u9R\nZ19uTrJ45vZHVV2a5O8l+R+TrXZiHpsZfB8Yh6qqJN+UGXgf6FVVz09yaQZBYt79QJKPJvnpqvpc\nVR2rql+oqrEtOb2HPpjk22s4HKCqHpTBVYmR3gf2akGq3Rp5ZcwkD97m8Q8eb2kj6+nLZq/I4LL5\nLFwKHLk/w8tmz8/gUvOs6fn5XJrk7yS5IYPLzt+VwZ4xFyT5N5Mpc1d6VpRdGd72+O3hH6oLkvxS\na+0VE610crZ7H3hAVd2/tfYnU6hpXH4qgw8Vb552IT2q6ruSvDzJ32ytna4tt2iaK5cmeXwGKys/\nLYPfv/+Q5M8l+bEp1jWy1trNVbWc5NeHQejCJDdmcFVi12bqisRZjLoy5thW0pyAXdVWVS9O8kMZ\n7GVy78Sr6rdlf6rqzyT5z0le0Fr7oz2vqt/Zfj4HMvjj9A9bax9rrb05gzEF/2ivihvRtn2pqicm\n+RcZ3K55VAZjPr6/qn52z6qbvDN/sWb1vWBHw7EsL0nyrNba+rTrGdXwat0bk7y0tfZ7Zw5PsaRx\nOJDBbYPntNY+2gYbWf7TJD86K4OVd2s4lvA1Sf5VBuNxrkzynUleN0o7s3ZFomdlzM+P+Pi90r3K\nZ1X9swymz17Rttm8bApG7c9fTPIdGQwSO/PGcSBJqureJIdba7dtcd5e6fn5nEhy7/C2wRlrSR5c\nVRe21v50/GXuSk9fXpbkDRtuOX1qGP5el+leXem13fvAH894EN9WVf1wBoPentlae8+06+n0TUke\nneSRVfWLw2MHMrhjc2+SJ7fWfmtaxXU6keSO1toXNxxbyyAgPSTJ72151mx6cZLfbq29evj9J2uw\nfMP7q+pnWmu7+js6U1ckWt/KmB/c+PihJw2PT01nX1JVP5XkZ5Jc2Vr72KTr3K2O/qwl+asZzKR5\nxPDrxiTvHv737RMu+aw6fz4fyGAxto0OJzkxxRDR25eLMvhUtdHp4anz+Ilxq/eBJ2fK7wO9qmop\nya8kWRp+4p1Xf5zkr+Tr3wd+Kcmnh//94emV1u0DSb61qi7acOxwBr8/n5tOSd22ex9oGeXK0bRH\nlm4xivSHknw5yd/PYFra65J8IclfGP77G5K8fMPjL09ybwaXlg5ncInmVJKHz2FfXjSs/ekZfLo6\n83XxtPvS058tzp+1WRuj/nweksEsmtdkMD7iKRl8En7xHPblpUnuzmDK50MzCN+fSfKmafdlWN/F\nGfyheWQGb2w/Mfz+24f//vMZLIJ35vEPTfLFDMYVHU7ywuH7wt+dw74sDWu/etP7wAOm3Zee/mxx\n/kzN2uj4+Vyc5P8k+fUkRzKYrnssgzFG89aXH0nyJ8PX2ncmeVwGK1LfPNLzTrvj2/zPeGEG24R/\nOYNPFI/e8G/vTvKrmx7/gxkk3C8n+UQGn+an3o9R+5LktgwuUW/++pfT7kfvz2bTuTMVJDpfa389\ng0/592Twh/enM1whdtpfI77WDmRw7/13k3xpeN5rMzt/rJ4wfCPc/LvwqxteS+/e4pxbhv3/TJLn\nTbsfPX3JYN2Ird4Htv3dmuX+bHH+rAWJntfad2cwK+iLGYSKVya5/5z25ZokvzPsy+cyWFfiW0Z5\n3pGXyAYAOGOmxkgAAPNFkAAAugkSAEA3QQIA6CZIAADdBAkAoJsgAQB0EyQAgG6CBADQTZAAALoJ\nEgBAN0ECAOj2/wBVzmBZjnnmUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c604e90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bin_widths = [bins[i+1]-bins[i] for i in range(len(bins)-1)]\n",
    "plt.bar(left=bins[:-1], height=hist1D, width=bin_widths, log=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Committor histogram in 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ramachandran_hash = lambda x : (float(phi(x)), float(psi(x)))\n",
    "hist2D, bins_phi, bins_psi = committor_analyzer.committor_histogram(ramachandran_hash, alpha_R, bins=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAFkCAYAAAAe3CMfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XucXVV99/HPlxguCRCoEcIdoxChtpEMUOOVEgEpaiva\n4ijVgjXSprXPWEqV1iJ44cEKqT4lJdQWEtF5FW1F0JYoKFgMCCQErSZabgZICEQgXJJASH7PH3tP\nOHNyZubsNbPP2efM9/167RecPft31lpzZvKbtfbaaykiMDMzs/bbqd0VMDMzs4yTspmZWUU4KZuZ\nmVWEk7KZmVlFOCmbmZlVhJOymZlZRTgpm5mZVYSTspmZWUU4KZuZmVWEk7KZmVlFFE7Kkt4o6VpJ\nD0vaJukdTcS8T9IKSc9KWiPpXyT9WlqVzczMypeY746TtEzSZkm/kPSBImWm9JQnAyuAecCIC2dL\nej2wCPhn4Ejg3cCxwOUJZZuZmbVK0Xx3KPAt4EZgJvAF4EuSTmi2QI1mQwpJ24Dfi4hrh7nmL4Gz\nIuKwmnN/BpwTEQcnF25mZtYiTea7i4CTI+I3a871A1Mi4neaKacV95RvBQ6SdDKApH3JesvfbkHZ\nZmZmrfJa4Ia6c0uA2c2+wUvGtDoNRMRSSacD/yZp17zMa4E/GypG0kuBk4AHgM1l19HMzEqzK3Ao\nsCQifjXWby7pYGBqYvj6iFg9htWZBqyrO7cO2FPSLhHx3EhvUHpSlnQk2bj6J4HvAPsBnwcWAn88\nRNhJwFfKrpuZmbXM+4CvjuUbSjqYSRN/ycYtqW/xnKTDxzgx11P+36buFZeelIGPAbdExCX56/+R\n9KfAf0v6m4io/6sCsh4yV111FUcccUQLqlgdfX19zJ8/v93VaDm3e3xxu8ePlStXcvrpp0P+7/oY\nm8rGLXDVqXBEwc7yyvVw+n/sQtbLHquk/Aiwb925fYCnIuL5Zt6gFUl5ElD/Z8w2sr8atOPlQD5k\nfcQRRzBr1qwSq1Y9U6ZMGXdtBrd7vHG7x6XybkW+aioctX+xmPQ5zsO5FTi57tyJ+fmmpDynPFnS\nTEmvyU9Nz18flH/9QkmLakKuA06VdJakl+ePSH0B+FFEPFK0fDMzs0FCaccIEvLdZcArJF0kaUY+\nKvxu4BKalNJTPhr4PtnfGQFcnJ9fBJxJdqP7oIGLI2KRpN3JnvP6PPAk2TNcH0so28zMrFWK5rsH\nJJ1CloQ/AjwEfDAi6mdkD6lwUo6Imxmmhx0RZzQ4dylwadGyzMzMRjSQMovGjHRJWr67GegpWJvt\nWnFP2Qro7e1tdxXawu0eX9xuG1NNDkfvEFNBTsoVM15/ad3u8cXttjFXzsStlnNSNjOzzuaespmZ\nWUWUdE+5HbyfspmZWUW4p9wmqubIyejs+kJa3OYu/THc75niMWt3H/t6jKXdm1qUaEfP7Fw8JuX7\nB9X/Ho7CKDb1624evjYzM6uILhq+dlI2M7PO5qRsZmZWIRUdji7KSdnMzDpbF91T9uxrMzOzinBP\n2czMOpvvKZuZmVVEFw1fOymbmVlnc0/ZzMysKhJ6ylSzp+yJXmZmZhXhnrKZmXU2D1+bmZlVhCd6\nmZmZVYR7ymZmZhURJPSUS6nJqDkpd5K5y4rHXDcjraxHJxeP2SNxW7/dErZ8TN0mspXb+nXjFoKp\nn3HK1o1PJ8TY+FVSkpU0DzgbmAbcDfx5RNwxxLUvAc4F3g8cAKwCPhYRS5otz7OvzczMGpB0GnAx\ncB5wFFlSXiJp6hAhnwE+BMwDjgAWAt+QNLPZMp2Uzcyssw1M9Cp6jKwPWBgRiyNiFXAWsBE4c4jr\nTwc+ExFLIuKBiLgM+E/gL5ttSuGkLOmNkq6V9LCkbZLe0UTMzpI+I+kBSZsl3Sfpj4qWbWZmtoNI\nPIYhaSLQA9y4vZiIAG4AZg8RtgvwXN25TcAbmm1Kyj3lycAK4F+Bf28y5mvAy4AzgHuB/XAv3czM\nxkI5j0RNBSYA6+rOrwOGmqyzBPiopP8my3VvAU6lQL4rnJQj4nrgegBJI7ZK0luBNwLTI+LJ/PTq\nouWamZk1NFLP9xt3wTV3DT731ObU0jRMaX8BXE42wWsbWWL+V7IOaVNaMfv67cCdwF9L+kPgWeBa\n4BMRkfxdMTMzy4zQU/69WdlR6ycPwVvnD/em64GtwL515/dhx94zABGxHjhV0s7ASyNiraT/C9w/\nQgO2a8UQ8nSynvKvA79H9pfEu4F/bEHZZmZmhUXEFmAZMGfgXD46PAdYOkLs83lCngi8C7im2XJb\n0VPeiawb/96IeAZA0keBr0maFxH1N8XNzMyaV96KXpcAiyQtA24nm409CbgSQNJi4KGIODd/fSzZ\n88krgAPJHqUS8PfNVqsVSXkt8PBAQs6tJKvogWRj7g319fUxZcqUQed6e3vp7e0to55mZjYK/f39\n9Pf3Dzq3YcOG8gsuae3riLg6fyb5ArJh7BXASRHxWH7JgUDtSka7Ap8GXg48A3wbOD0inmq2Wq1I\nyj8E3i1pUkRszM/NIOs9PzRc4Pz585k1a9Zwl5iZWUU06jQtX76cnp6ecgsuce3riFgALBjia8fX\nvf4B2a3aZCnPKU+WNFPSa/JT0/PXB+Vfv1DSopqQrwK/Aq6QdISkNwGfA/7FQ9dmZjZqA2tfFzra\nXenGUiZ6HQ3cRXYDPMiWIFsOnJ9/fRpw0MDFEfEscAKwF3AH8GXgm2QTvszMzEZvDBcOaaeU55Rv\nZphkHhE7PI8VEb8ATipalpmZ2XjiXaLMzKyzlTTRqx2clDvJ5QmTJXZP3Gpva8IPbEoMpG3DmLqt\nX0pZmxN/TfZOWBtnS+LSATN+VTxmzR7FY1q5HeU+z6bFpWwTufPWtLKen5AWZ2OrxIlereakbGZm\nnc09ZTMzs4roop6yd2oyMzOrCPeUzcyswyUMX+PhazMzs7HXRcPXTspmZtbZPNHLzMysItxTNjMz\nq4guSsqefW1mZlYR7imbmVmH8+xrMzOzauii4WsnZTMz62wD+ykXjakgJ2UzM+tsfiTKzMysIjx8\nbR3juAfS4n54cPGY1G3sUv5g3TdxW78J24rHpLbruYRfr723pJW1bL/iMfs9UzwmZTtKSNsy85HE\nbSJbuT2n2RjzI1FmZtbh9OIQdrNHk70BSfMk3S9pk6TbJB0zwvX/R9IqSRslrZZ0iaRdmm2Jk7KZ\nmXW2SDxGIOk04GLgPOAo4G5giaSpQ1z/XuDC/PpXAWcCpwGfabYpTspmZtbZivaSm58Y1gcsjIjF\nEbEKOAvYSJZsG5kN3BIR/xYRqyPiBqAfOLbZpjgpm5lZZyuhpyxpItAD3Li9mIgAbiBLvo0sBXoG\nhrglTQd+B/h2s03x7AYzM+ts5TwSNRWYAKyrO78OmNHwLSP686HtWyQpj78sIi5qtlpOymZm1t2u\nvx2W3DH43NObUt9NDNHPlnQccC7ZMPftwCuBL0paGxGfbubNnZTNzKzzDTccfdKx2VFr1Wo4fdj5\nV+uBrcC+def3Ycfe84ALgMURcUX++qeSdgcWAk0l5cL3lCW9UdK1kh6WtE3SOwrEvl7SFknLi5Zr\nZmbW0MAym4WOEd4yYguwDJgzcC4fkp5Ddu+4kUlA/WII2/LQpsbXUyZ6TQZWAPMosCaKpD2BRWQ3\nyc3MzMZGSY9EAZcAcyW9X9KrgMvIEu+VAJIWS/pszfXXAX8i6TRJh0o6gaz3/M18ktiICg9fR8T1\nwPV5hYrcWV8IfIXsr4bfLVqumZlZQyWtfR0RV+cTty4gG8ZeAZwUEY/llxwI1C4h9ymyHPcp4ADg\nMeBa4G+brVZL7ilLOgOYDrwP+EQryjQzMxutiFgALBjia8fXvR5IyJ9KLa/0pCzpMOCzwBsiYlux\nzrWZmdkIvCFFcyTtRDZkfV5E3Dtwutn4vr4+pkyZMuhcb28vvb29Y1dJMzMbE/39/fT39w86t2HD\nhhaUnDB8nbQTTvnK7invARwNvEbSpfm5nchuRz8PnBgRNw0VPH/+fGbNmlVyFc3MbCw06jQtX76c\nnp6ecgt2T7lpTwGvrjs3D/ht4F3AAyWXb/+zT1rcpoQfjZcmPoy/U8Jvx9M7p5W1y9biMS8krka7\nNmHrwd2fTytr54R2pXwPU78XzyVsf7lxYlpZ058oHrMpsayUz9jGXkkTvdqh8L+8kiaTrVIy0KLp\nkmYCj0fEg5IuBPaPiA/kU8B/Vhf/KLA5IlaOsu5mZmbjvqd8NPB9Xvw2XJyfX0S2c8Y04KAxqZ2Z\nmdk4kvKc8s0Ms+hIRJwxQvz5wPlFyzUzM2toPA9fm5mZVco4H742MzOrlor2fItyUjYzs87mnrKZ\nmVlFdNE95cSHDs3MzGysuadsZmadzcPXZmZmFdFFw9dOymZm1vkq2vMtyknZzMw6m3vKZmZmFdFF\n95Q9+9rMzKwi3FPuJPs9Uzzm1Y+mlfVIwpZ0zydsz5cat3cLt4l8cM+0slIc90Ba3M9eVjxmzR7F\nYw5O3LA+5edp4ra0sh5K+LxSf3Z71haPWbZfWlk2NA9fm5mZVYSHr83MzCokCh5NkjRP0v2SNkm6\nTdIxw1z7fUnbGhzXNVuek7KZmXW24MUh7KaPkd9W0mnAxcB5wFHA3cASSVOHCHknMK3meDWwFbi6\n2aY4KZuZmTXWByyMiMURsQo4C9gInNno4oh4MiIeHTiAE4Fnga83W6CTspmZdbaiQ9dNDGFLmgj0\nADduLyYigBuA2U3W7EygPyKanpnqiV5mZtbZypl9PRWYAKyrO78OmDFSsKRjgV8HzihSLSdlMzPr\nbCMl5ZuXwg+WDj737MbU0kRzU8U+CPxPRCwr8uZOymZm1vmGS5Nvel121Lr3fuj7m+HecT3ZJK19\n687vw46950Ek7QacBvztcNc14nvKZmbW2QrPvB55uDsitgDLgDkD5yQpf710qLjcacDOwFeKNsU9\nZTMzs8YuARZJWgbcTjYbexJwJYCkxcBDEXFuXdwHgWsi4omiBTopm5lZZytpRa+IuDp/JvkCsmHs\nFcBJEfFYfsmBwAu1MZIOA14HnFCwRkDC8LWkN0q6VtLD+Uol7xjh+ndK+o6kRyVtkLRU0okplTUz\nM9tBCcPX2986YkFEHBoRu0XE7Ii4s+Zrx0fEmXXX/29ETIiI76U0JeWe8mSyvxbm0dzfJm8CvgOc\nDMwCvg9cJ2lmQtlmZmaDlfCccrsUHr6OiOuB62H7Te+Rru+rO/U3kn4XeDvZkmVmZmajU9Fdn4pq\n+T3lPJHvATze6rI73h7PFY9Ztn9aWUevSYtL8ULCgM19e6eVtTWhrCMfG/maRp7YrXhM6ufVk/B5\n7by1eMxLErdT3P/p4jEpPxfQ2m1HvQ1jNXiXqFH5K7Ih8KYX6DYzMxsPWtpTlvRe4BPAOyJifSvL\nNjOzLlXOMptt0bKkLOk9wOXAuyPi+83E9PX1MWXKlEHnent76e3tLaGGZmY2Gv39/fT39w86t2HD\nhvIL7qLh65YkZUm9wJeA9+QTxZoyf/58Zs2aVV7FzMxszDTqNC1fvpyenp5yCx7PPWVJk4FXki3K\nDTA9f7zp8Yh4UNKFwP4R8YH8+l5gEfAR4HZJA+uIboqIp0bdAjMzs4r2fItKmeh1NHAX2ZqgAVwM\nLAfOz78+DTio5vq5ZNtfXQqsqTn+Ia3KZmZm3SnlOeWbGSaZR8QZda9/O6FeZmZmzRnPw9dmZmaV\n4oleZmZmFeGespmZWUW4p2xmZlYhFe35FtWOZTbNzMysAfeUzcyss3n42szMrCI80cva4hcvLR5z\n+K/SykrZyu7ZndPK2pJwF+WgxMXgfjll5GvqPbxnWlkHJNRx781pZW3YtTUxKVswQtpWmxMTt4k8\n4b7iMTcdmlbWEwnfQxt77imbmZlVhHvKZmZmFVLRnm9Rnn1tZmY2BEnzJN0vaZOk2yQdM8L1UyRd\nKmlNHrNK0lubLc89ZTMz62wlDV9LOo1s06W5wO1AH7BE0uERsb7B9ROBG4BHgFPJNl86BHiy2Wo5\nKZuZWWcrb6JXH7AwIhYDSDoLOAU4E/hcg+s/COwFvDYitubnVheploevzcyssw30lIsew8h7vT3A\njduLiQiynvDsIcLeDtwKLJD0iKSfSPq4pKZzrXvKZmbW2crpKU8FJgDr6s6vA2YMETMdOB64CjgZ\nOAxYkL/Pp5uplpOymZl1vuF6vrffDHf8YPC5Tc+mliSGTuk7kSXtuXmv+i5JBwBn46RsZmYGHPvm\n7Ki1+h74TN9wUeuBrcC+def3Ycfe84C1wPN5Qh6wEpgm6SUR8cJIVfU9ZTMz62yReAz3lhFbgGXA\nnIFzkpS/XjpE2A+BV9admwGsbSYhg5OymZl1uhKScu4SYK6k90t6FXAZMAm4EkDSYkmfrbn+n4CX\nSvqCpMMknQJ8HPjHZpvi4WszM+tsJT2nHBFXS5oKXEA2jL0COCkiHssvORB4oeb6hySdCMwH7gYe\nzv+/0eNTDTkpm5lZZytxQ4qIWEA2g7rR145vcO5HwOsK1mY7D1+bmZlVhHvK7dKztnjMsv2Kxxz3\nQPEYSNvK7mUb08pau3vxmJP/N62sf3ht8ZgZidtfpmybeeuBaWV99saRr6n36TcVj7nn14rHALzu\nweIxN0xPK+uRhJ+nmY+klZW65aONsYTha7xLlJmZ2djroq0bCw9fS3qjpGslPSxpm6R3NBFznKRl\nkjZL+oWkD6RV18zMrE55s69bLuWe8mSyGWjzaKJZkg4FvkW2fuhM4AvAlySdkFC2mZnZYCWsfd0u\nhYevI+J64HrY/iD1SP4EuC8izslf/1zSG8h23/hu0fLNzMwGKXH2dau1Yvb1a8l21ai1hKF32TAz\nMxuXWjHRaxqNd9nYU9IuEfFcC+pgZmbdKkiY6FVKTUatXbOvB757w35b+vr6mDJlyqBzvb299Pb2\nllUvMzNL1N/fT39//6BzGzZsaE3hFU2yRbUiKT9C4102noqI54cLnD9/PrNmzSqtYmZmNnYadZqW\nL19OT09PuQV30SNRrUjKt5Jt9lzrxPy8mZnZ6IzniV6SJkuaKek1+anp+euD8q9fKGlRTchlwCsk\nXSRphqQ/Bd5NtvuGmZnZ6HTRI1Eps6+PBu4i22cygIuB5cD5+denAQcNXBwRDwCnAG8he765D/hg\nRNTPyDYzMxvXUp5TvplhknlEnDFETMk3FczMbFzqouFrr31tZmYdzhtSmJmZVYN7yjZqa/Zodw2G\n95v16700YV3ClnkAMxPKWjU1rayT7i0e88BeaWX1v7p4zH7PpJX13ncVj0nZPnTfxPo9Orl4zNxl\naWVdN6N4zMStaWVZNfiRKDMzs4roop5yK9a+NjMzsya4p2xmZp2ti9a+dk/ZzMw6XxQ8miRpnqT7\nJW2SdJukY4a59gOStknamv93m6SNRZrhpGxmZp2tpBW9JJ1GtkDWecBRwN3AEknDzTTdQLaI1sBx\nSJGmOCmbmVlnK9pLbr633AcsjIjFEbEKOAvYCJw5XG0i4rGIeDQ/HivSFCdlMzOzOpImkq1EeePA\nuYgI4AZg9jChu0t6QNJqSddIOrJIuU7KZmbW2coZvp4KTADqF1JYRzYs3cjPyXrR7wDeR5Zjl0o6\noNmmePa1mZl1tpGGo39yI/zke4PPbU5cCCdbn7NhaRFxG3Db9gulW4GVwFyy+9IjclI2M7MON0LP\n99VvyY5aa38Bl88d7k3XA1uBfevO78OOveeGIuIFSXcBr2zmevDwtZmZdboSJnpFxBayLYrnDJyT\npPz10maqJWkn4NVA02vauqdsZmadrby1ry8BFklaBtxONht7EnAlgKTFwEMRcW7++hNkw9f3AHsB\n55A9EvWlZqvlpGxmZtZARFydP5N8Adkw9grgpJrHnA4EXqgJ2Ru4nGwi2BNkPe3Z+eNUTXFSNjOz\nzlbihhQRsQBYMMTXjq97/VHgowVrMoiTcrvs/3TxmIuXFI9J2dIPoPd/iscc+3BaWc9OLB7zs5el\nlbXPs8Vj9kv4rLLA4iGHPplWVM+a4jEvK7T6X2Zt4vacu70w8jX1NiX8XEDaNoyrp6SVZdXQRWtf\nOymbmVnnq2iSLcpJ2czMOlt5E71azknZzMw6W4n3lFvNzymbmZlVhHvKZmbW2cZ7T7nIps/59f9H\n0ipJG/OdMy6RtEtalc3MzGqUtJ9yOxTuKdds+jyXF1c4WSLp8IhY3+D69wIXAn8E3AocDiwCtgFn\nJ9fczMwMGHHt66FiKiilp1x00+fZwC0R8W8RsToibgD6gWOTamxmZlarhLWv26VQUk7c9Hkp0DMw\nxC1pOvA7wLdTKmxmZjbIOB6+Hm7T5xmNAiKiP1879JZ8h40JwGURcVHRypqZmXWzsXokashNnyUd\nB5xLNsx9FHAq8DZJfztGZZuZ2XjWRcPXRXvKKZs+XwAsjogr8tc/lbQ7sBD49HCF9fX1MWXK4DVp\ne3t76e3tLVhtMzMrW39/P/39/YPObdiwofyCx+uKXhGxJd9Xcg5wLQza9PmLQ4RNIptpXWtbHqr8\nnnRD8+fPZ9asWUWqaGZmbdKo07R8+XJ6enrKL7yiPd+iUhYPKbTpM3Ad0CdpBfAj4DCy3vM3h0vI\nZmZmTRnPu0QlbPr8KbKe8aeAA4DHyHrZ4/ue8rKEbf1StmGcd0fxGIDHJrUmBuDx3YrHHJw4JPbQ\nnsVjjknYFhFgxq+Kx9x0aFpZ3zq8eMzcZcVjUr/va/YoHvNgwmcFsFNF/7U1a0LSMpsFN30eSMif\nSinLzMxsWF20zKbXvjYzs842Xid6mZmZVY57ymZmZlUxvte+NjMzq44SFw8puitiTdx7JG2T9B9F\nmuKkbGZm1kDNrojnka1IeTfZrohTR4g7BPh74AdFy3RSNjOzzlbehhRFd0VE0k7AVcDfAfcXbYqT\nspmZdbYShq8Td0WErFf9aM3S0oV4opeZmXW2ch6JKrwroqTXA2cAM4tV5kVOymZm1vmG6/necz3c\ns2TwueefSS2p4a6I+UZLXwY+FBFPpL65k7KZmXW2kda+fsXJ2VFr/Ur4xunDvWvRXRFfARwCXJdv\n1AT5LWJJzwMzImLEe8y+p2xmZlYnIrYAA7siAoN2RVzaIGQl8BvAa8iGr2eS7fPwvfz/H2ymXPeU\nzcyss5W3olfTuyJGxPPAz2qDJT1JNj9sZbPVclI2M7POVtLa1wm7Io6ak3K3WzEtLW7y88Vj3tX0\nH4ODrRr2OfzGVk9JKytlu8IXEu/yXNdwgubw3nJvWlkpW4HecnDxmJ61xWMAtkwoHpO6jeXhCVtm\nWmcrce3rIrsiNvj6GQVr5aRsZmYdzrtEmZmZVUhFd30qyrOvzczMKsI9ZTMz62wevjYzM6uIEid6\ntZqTspmZdTb3lM3MzCqii3rKnuhlZmZWEe4pm5lZ56vocHRRTspmZtbZxvvwtaR5ku6XtEnSbZKO\nGeH6KZIulbQmj1kl6a1pVTYzM6sxMNGr6FFBhXvKkk4DLgbm8uKuGUskHR4R6xtcPxG4AXgEOBVY\nQ7bn5JOjqLeZmVmmi3rKKcPXfcDCiFgMIOks4BTgTOBzDa7/ILAX8NqI2JqfW51QrpmZ2Y66KCkX\nGr7Oe709wI0D5yIiyHrCs4cIeztwK7BA0iOSfiLp45I889vMzKxG0Z7yVGACsK7u/DpgqH3qpgPH\nA1cBJwOHkW2DNQH4dMHyx7dP3lQ85s7908ranDCIsv/TaWWlbN241+a0slK3A0zxzoStLN/7rrGv\nx1Ce2K14zG5b0sr6zI0jX1Nv4taRr2lkWeLPfIqF3yoe07NmFAXOHUVsN0u5R9wl95SHIIYeDNiJ\nLGnPzXvVd0k6ADibEZJyX18fU6YM3je3t7eX3t7e0dfYzMzGVH9/P/39/YPObdiwofyCx/GKXuuB\nrcC+def3Ycfe84C1wPN5Qh6wEpgm6SUR8cJQhc2fP59Zs2YVrKKZmbVDo07T8uXL6enpKbfg8XpP\nOSK2AMuAOQPnJCl/vXSIsB8Cr6w7NwNYO1xCNjMza0oXPRKVMtnqEmCupPdLehVwGTAJuBJA0mJJ\nn625/p+Al0r6gqTDJJ0CfBz4x9FV3czMjBd7ykWPCip8TzkirpY0FbiAbBh7BXBSRDyWX3Ig8ELN\n9Q9JOhGYD9wNPJz/f6PHp8zMzMatpMeSImJBRBwaEbtFxOyIuLPma8dHxJl11/8oIl4XEZMi4rCI\nuKjuHrOZmVm6koaui6xgKemdku6Q9ISkZyTdJen0Is3ws8JmZtbZShq+rlnB8jzgKLLR3iX5aHEj\nvyJ7qui1wG8AVwBXSDqh2aY4KZuZWWcrb6LX9hUsI2IVcBawkWwFyx2rEfGDiPhmRPw8Iu6PiC8C\nPwbe0GxTnJTNzKyzldBTTlzBsv495gCHAzc32xRv3WhmZp2tnMVDUlawRNKeZBOadyGb9PynEfG9\nZqvlpGxmZt3t4W9lR60ticsCD7+CJcDTwExgd7I1POZLui8iftDMmzspm5lZ5xsuTe7/tuyoteGn\n8N+nDveOKStYDgxx35e//LGkI8nW5mgqKfuespmZdbYSJnolrmDZyE5kQ9lNcU+5XVJ2fPrkca0p\nB2BTwo/G/zs2raw3/7J4zNSNaWW94vHiMe95d1pZPWuLx7ztF2llHZ2w81BK/Z7euXgMwPQnisek\n7jp29J0jX1Mv5XcLYNl+xWM+/LaRrxmKV3dorLy1ry8BFklaBtxONht70AqWwEMRcW7++mPAncC9\nZIn4FOB0slnbTXFSNjOzzlbSLlFFV7AEJgOX5uc3AauA90XE15utlpOymZnZECJiAbBgiK8dX/f6\nE8AnRlOek7KZmXW2Ltq60UnZzMw6X0W3YizKSdnMzDqbe8pmZmYVUdJEr3ZwUjYzs87WRT1lLx5i\nZmZWEe4pm5lZZ/PwtZmZWYVUdDi6KCdlMzPrbO4pm5mZVUQXTfRyUjYzs87WRT1lz742MzOrCPeU\n2yVli76U7eVStucDuPHlxWO+84q0st5y38jX1JvyXFpZL9lWPOar/55W1r8eVTzmhulpZV1+XfGY\n/f+yeMwurOdKAAAQM0lEQVTnv1M8BuBrv54Wl+LO/VtX1uU9rSvLhtZFw9dJPWVJ8yTdL2mTpNsk\nHdNk3HskbZP0HynlmpmZNTQwhN3sUVGFk7Kk04CLgfOAo4C7gSX5npPDxR0C/D3wg4R6mpmZNRaJ\nRwWl9JT7gIURsTgiVgFnARuBM4cKkLQTcBXwd8D9KRU1MzNrqGgvucK95UJJWdJEoAe4ceBcRARw\nAzB7mNDzgEcj4oqUSpqZmQ2pi3rKRSd6TQUmAOvqzq8DZjQKkPR64AxgZuHamZmZjSNjNftaNPi7\nQ9LuwJeBD0XEE2NUlpmZ2Yu66Dnlokl5PbAV2Lfu/D7s2HsGeAVwCHCdpIHvwE4Akp4HZkTEkPeY\n+/r6mDJlyqBzvb299Pb2Fqy2mZmVrb+/n/7+/kHnNmzY0JrCKzocXVShpBwRWyQtA+YA1wLkyXYO\n8MUGISuB36g79xlgd+AjwIPDlTd//nxmzZpVpIpmZtYmjTpNy5cvp6en5Oe5S3xOWdI84GxgGtnT\nRn8eEXcMce0fA+8HXp2fWgacO9T1jaQMX18CLMqT8+1ks7EnAVfmlVoMPBQR50bE88DP6ir9JNn8\nsJUJZZuZmQ1W0vB1zSPAc3kx3y2RdHhErG8Q8mbgq8BSYDPwMeA7ko6MiKZWciqclCPi6vyZ5AvI\nhrFXACdFxGP5JQcCLxR9XzMzs4rZ/ggwgKSzgFPIHgH+XP3FEfGHta/znvO7yEaTr2qmwKSJXhGx\nAFgwxNeOHyH2jJQyzczMGiqhp1zzCPBnt4dEhKSRHgGuNRmYCDzebLW8IYWZmXW2cp5THu4R4GlN\n1uwi4GGytTya4g0pzMys8w3X8338G/D4NYPPbX0qtaSGjwDvcJH0MeAPgDfn86ua4qRsZmadbaSe\n797vzI5aG38Mq9463LsWfQR4O0lnA+cAcyLip8NdW89JuV2Onls85rr+ka+pt2aP4jEA//Da4jF/\nfUtaWSnbMM5oNPGxCdc1XHhueBc3e/uozsJvFY95w+q0slLaleK4B9LiliVsp9jKbRHj/LQ4nTe2\n9bA0JdxTTngEmPyavwLOBU6MiLuKVcpJ2czMbChNPwKcvz6H7MmkXmC1pIFe9jMR8WwzBTopm5lZ\nZytp8ZCER4D/hGy29dfr3ur8/D1G5KRsZmadrcS1r4s8AhwRLy9WiR05KZuZWWcrcZnNVnNSNjOz\nDpfQU6Y7dokyMzOrli7qKXtFLzMzs4pwT9nMzDpbiRO9Ws1J2czMOlsXDV87KZuZWWcLEnrKpdRk\n1JyUzcys81U0yRbliV5mZmYV4Z6ymZl1Nk/0MjMzqwhP9LK2eHtv8ZjvfjmtrHuH3JlsaDdMTyvr\n1x8tHnPPr6WVdfCGtLgUH35b68pK8fnvFI+56dC0ss4+sXhMSv0APryseIy3YOxs7imbmZlVhHvK\nZmZmVdE9a1979rWZmVlFuKdsZmadzcPXZmZmFdFFE72Shq8lzZN0v6RNkm6TdMww1/6xpB9Iejw/\nvjvc9WZmZoVE4lFBhZOypNOAi4HzgKOAu4ElkqYOEfJm4KvAccBrgQeB70jaL6XCZmZmgwysfV3o\naHelG0vpKfcBCyNicUSsAs4CNgJnNro4Iv4wIi6LiB9HxC+AP87LnZNaaTMzs0G6oJcMBZOypIlA\nD3DjwLmICOAGYHaTbzMZmAg8XqRsMzOzVit4u/ZISV/Pr98m6SNFyyvaU54KTADW1Z1fB0xr8j0u\nAh4mS+RmZmajU3jourmJYQm3aycB9wJ/DaxNacpYzb4WTQwISPoY8AfAmyPi+ZGu7+vrY8qUKYPO\n9fb20tubsNykmZmVqr+/n/7+/kHnNmxowdK25T0Stf12LYCks4BTyG7Xfm6Ht4y4E7gzv/aigjUC\niifl9cBWYN+68/uwY+95EElnA+cAcyLip80UNn/+fGbNmlWwimZm1g6NOk3Lly+np6en3IJLeCSq\n5nbtZ7eHRISkIrdrCys0fB0RW4Bl1EzSkqT89dKh4iT9FfA3wEkRcVdaVc3MzBoo55GosbhdW1jK\n8PUlwCJJy4Dbybr3k4ArASQtBh6KiHPz1+cAFwC9wGpJA73sZyLi2dFV38zMbISe8qavZUeteGoU\nhZU3f7twUo6Iq/Ob3BeQDWOvIOsBP5ZfciDwQk3In5DNtv563Vudn7+HlemEP2xdWWsubl1Z9+2d\nFveRHxWP+Yu3ppV1Xf/I19RL2Z4zVcp2iiltStXKbSJb6aYrRxH8R2NUiXFmt9/PjlpbVsD6Nw4X\nlXy7djSSJnpFxAJgwRBfO77u9ctTyjAzM2tKCRO9ImJLPiI8B7gWBt2uTdhwvjle+9rMzDpbeWtf\nF71dOxE4kmyIe2fgAEkzyW7X3ttMgU7KZmbW2Up6JCrhdu3+wF017352ftwMDBpFHoqTspmZdbaB\nta+LxjRzWbHbtb8kcaOnAaMKNjMzs7HjnrKZmXW+Cm8yUYSTspmZdbbyJnq1nJOymZl1tvLWvm45\nJ2UzM+ts7imbmZlVRBf1lD372szMrCLcUzYzsw6XMHyNh6/NzMzGXhcNXzspm5lZZ/NEL7MGnt4l\nMW7n4jHvXJVW1uU9xWNauSVlnJ8Wp/OKxzx9YfGYPT5ePAbS2pXSJoC5y4rHfPKmtLL2/8viMcf9\nUVpZUNneXdu5p2xmZlYRXdRT9uxrMzOzinBP2czMOl9Fh6OLclI2M7POVuLWja3mpGxmZp3NE73M\nzMwqoosmejkpm5lZZ+uinrJnX5uZmVWEk3Ll9Le7Au1x/e3trkF7fOOudtegPfp/0u4atMk4/f0u\nnV4cwm72qOja10lJWdI8SfdL2iTpNknHjHD970tamV9/t6ST06o7HozTX9old7S7Bu1xjZPy+DJO\nf7/LFolHE1qd7wonZUmnARcD5wFHAXcDSyRNHeL62cBXgX8GXgNcA1wj6ciiZZuZme2gaC+5yYlh\n7ch3KT3lPmBhRCyOiFXAWcBG4Mwhrv8L4L8i4pKI+HlEnAcsB/4soWwzM7PByusptzzfFUrKkiYC\nPcCNA+ciIoAbgNlDhM3Ov15ryTDXm5mZtVW78l3RR6KmAhOAdXXn1wEzhoiZNsT104YpZ1eAlStX\nFqxeN9hA9odVB/rpg2lxGyfCM5tg1ermYzY9mVbW6t2Kx/z4obSymvHU5sHvv3Zj4hsl/MyseLg1\n5QAsXzP49YbNO54bq7LWJ+wglvwZF63j6H6/l3fgPw01/47vWlohsSrhueMRf05ale8Gi4imD2A/\nYBvwW3XnPwcsHSLmOeC0unN/CqwZppz3kj4g4cOHDx8+qne8t0i+aTInHQw8O4o6bQYObme+qz+K\n9pTXA1uBfevO78OOfx0MeKTg9ZB1998HPED2TTMzs860K3Ao2b/rYyoiVks6gqxXm2J9RAw1RNeq\nfDdIoaQcEVskLQPmANcCSFL++otDhN3a4Osn5OeHKudXZDPYzMys8y0t643zpFrg3lfT79uSfFcv\nZZnNS4BFeWVvJ5udNgm4EkDSYuChiDg3v/4LwM2SPgp8G+glu3n+oYSyzczMWqXl+a5wUo6Iq/Nn\ntC4g66avAE6KiMfySw4EXqi5/lZJvcBn8uN/gd+NiJ8VLdvMzKxV2pHvlN+INjMzszbz2tdmZmYV\n4aRsZmZWEZVIypK+KemX+QLeayQtlrTfCDE3SdpWc2yVtKBVdR4Lie3eRdKlktZLelrS1yXt06o6\nj5akQyR9SdJ9kjZK+l9Jn8xXzxkurqM/71G0u6M/bwBJ50r6oaRnJT3eZMwVdZ/3Nkn/WXZdx1JK\nu/O4C/J/DzZK+q6kV5ZZz7EmaW9JX5G0QdIT+c/95BFiOvr3eyxVIikD3wN+HzgcOBV4BfC1EWIC\nuJzs5vs0sge9zymxjmVIafc/AKcA7wLeBOwP/HuJdRxrryLbM+1DwJFksxnPIpsUMZxO/7xT293p\nnzfAROBq4J8Kxv0XL37e08hmsnaSwu2W9Ndk6yR/GDiWbGGMJZJ2LqWG5fgqcATZo0GnkP3cLhwh\nptN/v8fOWK+wMkartLydbEbbhGGu+T5wSbvr2sp2A3uSrRjzzppzM8hWnTm23fUfRbvPBu4Z4Zpu\n/LyHbXe3fd7AB4DHm7z2CuA/2l3nNrR7DdBX9zOwCfiDdrejyfq/Kv/5PKrm3En5v2vThonrut/v\n1KMqPeXtJP0a2WpeP4yIrSNc/j5Jj0n6iaTPSkpY2Lgammx3D9ljbLULpP+c7MH5Tt7gYy+gmeG9\nrvm8cyO1u1s/72YdJ2mdpFWSFuS/I11L0svJeom1n/dTwI/onM97NvBERNRuFH4DWU/4t0aI7bbf\n7yQpi4eUQtL/JRu2mUS2+snbRgj5CvBLsr8sf5NsPdLDgXeXWM0xV7Dd04Dn81/UWsUWPK+Q/H7Z\nnwEfHeHSrvi8BzTZ7q77vAv4L7Jh+vvJbutcCPynpNmRd6260DSy5DW6DQ3aaxrwaO2JiNia31Mf\nrg1d9fs9GqX1lCVd2GCiRv2N/MNrQj5Htin0CWTrjX55uPePiC9FxHcj4qcR0Q+8H3hn/tdm25Td\n7qGKJftlbpuEdiPpALJ/fP8tIv51uPfvos+7ULuHKpYO/LyLiIirI+Jb+ed9Ldkfq8cCx41VG1KU\n3e6hiqXzP+9h21DV3+92KLOn/Hmy+0LDuW/gfyLicbKhvHskrQIelPRbEfGjJsv7EdkH/0qyv67b\npcx2PwLsLGnPut5ToQXPS1Ko3ZL2J5vodktEfDihvI78vAu2u2s+79GKiPslrSf7vL8/Vu+boMx2\nP0L2M70vgz/ffYC7Gka0TrPtfoSsvttJmgDsTbGf2ar8frdcaUk5sk0lfpUYPiH/7y4FYo4i+0ts\nbWKZY6Lkdi8jmzAxB/gGQP7X6cEUWPC8DEXanfcUvwfcAZyZWGTHfd4J7e6Kz3ssSDoQeCkd9Hkn\nvPf9kh4h+7x/DCBpT7J7sZeWUWaBujXVbkm3AntJOqrmvvIcsgTbbAcLKvL73RbtnmkGHAPMA2aS\n/WNzPHAL8HNgYn7N/sBK4Oj89XTgb4FZwCHAO4B7gO+1uz1ltjs/t4DsL8fjyCYC/RD473a3p0C7\n9yNbD/a7efv2HThqrunGz7twu7vh887bcFD+c/53wIb8/2cCk2uuWUW2RjDAZLLbOr+Vf95zgDvz\n783EdrenrHbnr88hS35vB34DuCb/udm53e0p0O7/zD+vY4DX5/+mfbnm6133+z2m37+2VwBeTTbb\n8DFgI3Av8I/AfjXXHEJ2v/VN+esDgZtqYn5ONhFk93a3p8x25+d2Af4f2V6fT5M917xPu9tToN0f\nyNtUe2wDtnb551243d3weedtuKJB2+vbuRV4f/7/uwLXkw2FbiYbFv0n4GXtbkuZ7a4590myCU8b\nyfYgfmW721Kw3XsBV5H9IfIE8M/ApJqvd93v91ge3pDCzMysIir3nLKZmdl45aRsZmZWEU7KZmZm\nFeGkbGZmVhFOymZmZhXhpGxmZlYRTspmZmYV4aRsZmZWEU7KZmZmFeGkbGZmVhFOymZmZhXx/wF9\nmtmuJnbVfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c6049d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# not the best, since it doesn't distinguish NaNs, but that's just a matter of plotting\n",
    "plt.pcolor(bins_phi, bins_psi, hist2D.T, cmap=\"winter\")\n",
    "plt.clim(0.0, 1.0)\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining information from the snapshots\n",
    "\n",
    "The information `committor_results.nc` should be *everything* you could want, including initial velocities for every system. In principle, you'll mainly access that information using collective variables (see documentation on using MDTraj to create OPS collective variables). However, you may decide to access that information directly, so here's how you do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# let's take the first shooting point snapshot\n",
    "# experiments[N][0] gives shooting snapshot for experiment N\n",
    "snapshot = experiments[0][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenMM-based objects come with units. So `snapshot.coordinates` is a unitted value. This can be annoying in analysis, so we have a convenience `snapshot.xyz` to get the version without units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Quantity(value=array([[ 0.2793301 , -0.12214785, -0.202672  ],\n",
       "       [ 0.2449614 , -0.05111128, -0.1274817 ],\n",
       "       [ 0.33534443, -0.01180077, -0.08093405],\n",
       "       ..., \n",
       "       [-1.46748817, -0.68787634,  0.40931037],\n",
       "       [-1.55375791, -0.72757316,  0.42131069],\n",
       "       [-1.45325804, -0.68954408,  0.31466872]], dtype=float32), unit=nanometer)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snapshot.coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.2793301 , -0.12214785, -0.202672  ],\n",
       "       [ 0.2449614 , -0.05111128, -0.1274817 ],\n",
       "       [ 0.33534443, -0.01180077, -0.08093405],\n",
       "       ..., \n",
       "       [-1.46748817, -0.68787634,  0.40931037],\n",
       "       [-1.55375791, -0.72757316,  0.42131069],\n",
       "       [-1.45325804, -0.68954408,  0.31466872]], dtype=float32)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snapshot.xyz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For velocities, we don't have the convenience function, but if you want to remove units from velocities you can do so with `velocity / velocity.unit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Quantity(value=array([[ -3.21658134e-01,  -1.10456896e+00,   7.57985592e-01],\n",
       "       [ -2.63174623e-01,   7.86199495e-02,  -5.87066524e-02],\n",
       "       [ -4.26795304e-01,   3.19573097e-03,  -2.02318802e-01],\n",
       "       ..., \n",
       "       [  2.87866861e-01,   5.24361193e-01,   2.04794154e-01],\n",
       "       [ -3.03164870e-01,   3.03924894e+00,  -1.35492578e-01],\n",
       "       [ -1.56592965e+00,  -4.29962158e+00,  -5.62753141e-01]], dtype=float32), unit=nanometer/picosecond)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snapshot.velocities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -3.21658134e-01,  -1.10456896e+00,   7.57985592e-01],\n",
       "       [ -2.63174623e-01,   7.86199495e-02,  -5.87066524e-02],\n",
       "       [ -4.26795304e-01,   3.19573097e-03,  -2.02318802e-01],\n",
       "       ..., \n",
       "       [  2.87866861e-01,   5.24361193e-01,   2.04794154e-01],\n",
       "       [ -3.03164870e-01,   3.03924894e+00,  -1.35492578e-01],\n",
       "       [ -1.56592965e+00,  -4.29962158e+00,  -5.62753141e-01]], dtype=float32)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "snapshot.velocities / snapshot.velocities.unit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that snapshots include coordinates and velocities. We have several sets of initial velocities for each initial snapshot. Taking the second shooting snapshot and comparing coordinates and velocities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "snapshot2 = experiments[1][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(snapshot.coordinates == snapshot2.coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.any(snapshot.velocities == snapshot2.velocities)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

5_mdtraj_tricks.ipynb