{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# (Deprecated) Analysing pore dimensions with HOLE\n",
"\n",
"Here we use `HOLE` to analyse pore dimensions in a membrane.\n",
"\n",
"
\n",
" \n",
"**Warning**\n",
"\n",
"The `MDAnalysis.analysis.hole` is now deprecated in favor of `MDAnalysis.analysis.hole2`, and will be removed in 2.0.0. Please see the [updated notebook](hole2.ipynb) for the updated tutorial.\n",
"
\n",
"\n",
"**Last executed:** May 18, 2021 with MDAnalysis 1.1.1\n",
"\n",
"**Last updated:** January 2020\n",
"\n",
"**Minimum version of MDAnalysis:** 0.18.0\n",
"\n",
"**Packages required:**\n",
" \n",
"* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n",
"* MDAnalysisTests\n",
"* [HOLE](http://www.holeprogram.org)\n",
"* matplotlib\n",
"* numpy\n",
"\n",
"
\n",
" \n",
"**Note**\n",
"\n",
"The classes in `MDAnalysis.analysis.hole` are wrappers around the HOLE program. Please cite (Smart *et al.*, 1993, Smart *et al.*, 1996) when using this module in published work.\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:27.248912Z",
"iopub.status.busy": "2021-05-19T06:17:27.248178Z",
"iopub.status.idle": "2021-05-19T06:17:27.972046Z",
"shell.execute_reply": "2021-05-19T06:17:27.972693Z"
}
},
"outputs": [],
"source": [
"import MDAnalysis as mda\n",
"from MDAnalysis.tests.datafiles import PDB_HOLE\n",
"from MDAnalysis.analysis import hole\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Using HOLE with a PDB file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`MDAnalysis.analysis.hole.HOLE` is a wrapper class ([API docs](https://docs.mdanalysis.org/stable/documentation_pages/analysis/hole2.html)) that calls the [HOLE](http://www.holeprogram.org) program. This means you must have installed the program yourself before you can use the class. You then need to pass the path to your ``hole`` executable to the class; in this example, ``hole`` is installed at ``~/hole2/exe/hole``.\n",
"\n",
"HOLE defines a series of points throughout the pore from which a sphere can be generated that does not overlap any atom (as defined by its van der Waals radius). (Please see (Smart *et al.*, 1993, Smart *et al.*, 1996) for a complete explanation). By default, it ignores residues with the following names: \"SOL\", \"WAT\", \"TIP\", \"HOH\", \"K \", \"NA \", \"CL \". You can change these with the ``ignore_residues`` keyword. Note that the residue names must have 3 characters. Wildcards *do not* work.\n",
"\n",
"The PDB file here is the experimental structure of the Gramicidin A channel. Note that we pass `HOLE` a PDB file directly, without creating a `MDAnalysis.Universe`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:27.977976Z",
"iopub.status.busy": "2021-05-19T06:17:27.977246Z",
"iopub.status.idle": "2021-05-19T06:17:28.473011Z",
"shell.execute_reply": "2021-05-19T06:17:28.473797Z"
}
},
"outputs": [],
"source": [
"h = hole.HOLE(PDB_HOLE, executable='~/hole2/exe/hole',\n",
" logfile='hole1.out',\n",
" sphpdb='hole1.sph',\n",
" raseed=31415)\n",
"h.run()\n",
"h.collect()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This will create several outputs in your directory:\n",
"\n",
" - **hole1.out**: the log file for HOLE. \n",
" - **hole1.sph**: a PDB-like file containing the coordinates of the pore centers.\n",
" - **simple2.rad**: file of Van der Waals' radii\n",
" - **run_n/radii_n_m.dat.gz**: the profile for each frame\n",
" - **tmp/pdb_name.pdb**: the short name of a PDB file with your structure. As `hole` is a FORTRAN77 program, it is limited in how long of a filename that it can read. \n",
" \n",
" \n",
" The pore profile itself is in a dictionary at `h.profiles`. There is only one frame in this PDB file, so it is at `h.profiles[0]`."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:28.519061Z",
"iopub.status.busy": "2021-05-19T06:17:28.517795Z",
"iopub.status.idle": "2021-05-19T06:17:28.522083Z",
"shell.execute_reply": "2021-05-19T06:17:28.522542Z"
},
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"425"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"len(h.profiles[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Each profile is a ``numpy.recarray`` with the fields below as an entry for each `rxncoord`: \n",
"\n",
" - **frame**: the integer frame number\n",
" - **rxncoord**: the distance along the pore axis in angstrom\n",
" - **radius**: the pore radius in angstrom"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:28.527327Z",
"iopub.status.busy": "2021-05-19T06:17:28.526704Z",
"iopub.status.idle": "2021-05-19T06:17:28.528966Z",
"shell.execute_reply": "2021-05-19T06:17:28.529476Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"('frame', 'rxncoord', 'radius')"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h.profiles[0].dtype.names"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can then proceed with your own analysis of the profiles. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:28.533622Z",
"iopub.status.busy": "2021-05-19T06:17:28.533019Z",
"iopub.status.idle": "2021-05-19T06:17:28.535180Z",
"shell.execute_reply": "2021-05-19T06:17:28.535547Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The pore is 42.4 angstroms long\n"
]
}
],
"source": [
"rxncoords = h.profiles[0].rxncoord\n",
"pore_length = rxncoords[-1] - rxncoords[0]\n",
"print('The pore is {} angstroms long'.format(pore_length))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Both `HOLE` and `HOLEtraj` (below) have the `min_radius()` function, which will return the minimum radius in angstrom for each frame. The resulting array has the shape (#n_frames, 2)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:28.540087Z",
"iopub.status.busy": "2021-05-19T06:17:28.539246Z",
"iopub.status.idle": "2021-05-19T06:17:28.542614Z",
"shell.execute_reply": "2021-05-19T06:17:28.543158Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"array([[0. , 1.19707]])"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"h.min_radius()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The class has a convenience `plot()` method to plot the coordinates of your pore."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-05-19T06:17:28.560013Z",
"iopub.status.busy": "2021-05-19T06:17:28.556353Z",
"iopub.status.idle": "2021-05-19T06:17:28.915271Z",
"shell.execute_reply": "2021-05-19T06:17:28.916045Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZfbA8e9JMmkQpEUEAgQBqUoRUFYs2EDsrrpgY0V/WHCxd3cV26JiWcWu2AXdFSwrKBasLCX0Lh0CSO/pk/P7Y25wjJMwCZncO5PzeZ55Mvfed+49ucY5vOW+r6gqxhhjTDji3A7AGGNM9LCkYYwxJmyWNIwxxoTNkoYxxpiwWdIwxhgTtgS3A4i0hg0bamZmptthGGNM1Jg5c+ZWVU0PdSzmk0ZmZiZZWVluh2GMMVFDRNaUdcyap4wxxoTNkoYxxpiwWdIwxhgTtpjv0wilsLCQ7Oxs8vLy3A6lTMnJyWRkZODz+dwOxRhj9quRSSM7O5u0tDQyMzMREbfD+QNVZdu2bWRnZ9OyZUu3wzHGmP1qZPNUXl4eDRo08GTCABARGjRo4OmakDGmZqqRSQPwbMIo4fX4jDE1U41NGsYYE6umfDqDDx7/JCLntqThki+++IK2bdvSunVrRowY4XY4xpgYMu2/Mxn3r88jcm5LGi7w+/0MHTqUiRMnsmjRIsaMGcOiRYvcDssYEyP8RcXEJ0Tm692ShgumT59O69atOfzww0lMTGTAgAF88klkqpLGmJrH7/cTnxAfkXPXyCG3wV646Q1WzF1dpeds1TmT65+5sszj69evp1mzZvu3MzIymDZtWpXGYIypufxFkUsaVtNwQah12W20lDGmqkSyearG1zTKqxFESkZGBuvWrdu/nZ2dTZMmTao9DmNMbLKaRozp0aMHy5YtY9WqVRQUFDB27FjOOecct8MyxsSISCaNGl/TcENCQgKjRo2ib9+++P1+Bg8eTMeOHd0OyxgTI6x5Kgb179+f/v37ux2GMSYGWfOUMcaYsBVb0jDGGBMue7gvAkINe/USr8dnjPEuf5GfuFjoCBeR0cBZwGZV7eTs+wBo6xSpC+xU1S4hPrsa2AP4gSJV7V7ZOJKTk9m2bZtnp0cvWU8jOTnZ7VCMMVEolkZPvQmMAt4u2aGqfyl5LyJPArvK+XwfVd16sEFkZGSQnZ3Nli1bDvZUEVOycp8xxlRUzIyeUtUfRCQz1DEJ/JP/YuDkSMfh8/lsRTxjTMyqKaOnjgc2qeqyMo4rMElEZorIkPJOJCJDRCRLRLK8XJswxphIqClJYyAwppzjx6lqN+AMYKiInFBWQVV9RVW7q2r39PT0qo7TGGM8LeZHT4lIAnAB8EFZZVR1g/NzMzAe6Fk90RljTHTxFxbFfE3jVGCJqmaHOigitUQkreQ9cDqwoBrjM8aYqOEvKiY+PgaShoiMAf4HtBWRbBG5yjk0gFJNUyLSREQmOJuNgJ9EZC4wHfhcVb+orriNMSaa+Iv8JPhiYMitqg4sY/9fQ+zbAPR33q8EOkc0OGOMiRGRfLjPK81TxhhjqkjMd4QbY4ypOjVlyK0xxpgqYEnDGGNM2Kx5yhhjTFhU1WoaxhhjwlNcXAxgScMYY8yB+YssaRhjjAmTv8gPYH0a1W3H5l3s27XP7TCMMaZCivcnDatpVKvLMq/j/UfGuR2GMcZUiDVPucSX5KOwoMjtMIwxpkKsecolvsQECvML3Q7DGGMqxG/NU+7wJfkozLeahjEmupQ0T9mEhdUsITGBokJLGsaY6GLNUy7xJSVQYM1TxpgoY81TLvEl+SiyjnBjTJQpKrSk4QrrCDfGRCNrnnKJdYQbY6KRPafhkoTEBHtOwxgTdUpqGpFaI7xak4aIjBaRzSKyIGjfAyKyXkTmOK/+ZXy2n4gsFZHlInJXpGP1JVnzlDEm+sTaNCJvAv1C7H9aVbs4rwmlD4pIPPA8cAbQARgoIh0iGWigecqShjEmusRU85Sq/gBsr8RHewLLVXWlqhYAY4FzqzS4UnyJCTZ6yhgTdWpKR/gNIjLPab6qF+J4U2Bd0Ha2sy8kERkiIlkikrVly5ZKBWQd4caYaFSSNGL5ifAXgVZAF2Aj8GSIMhJin5Z1QlV9RVW7q2r39PT0SgWV4LOOcGNM9Imp5qlQVHWTqvpVtRh4lUBTVGnZQLOg7QxgQyTjso5wY0w02t88FR+jzVMi0jho83xgQYhiM4A2ItJSRBKBAcCnkYzLOsKNMdFo/xPhERpymxCRs5ZBRMYAJwENRSQbuB84SUS6EGhuWg1c45RtArymqv1VtUhEbgC+BOKB0aq6MJKxWke4MSYa+Z2JViP1nEa1Jg1VHRhi9+tllN0A9A/angD8YThupPiSfBQV+ikuLiYuzvUKmTHGhKWkTyPBF5mvd/s2LENCYuCGl1T1jDEmGpQs6RDrQ249x5fkA7B+DWNMVPmtT8NqGtXKlxS44ZY0jDHRxF8YQ3NPRRNfSfOUdYYbY6KILcLkkt+apyxpGGOiR6SH3FrSKENJR7g9FW6MiSZFER5ya0mjDNYRboyJRsVenUZERGo5U5bHpETrCDfGRKGiwiLi4iRiz5eFfVYRiRORS0TkcxHZDCwBNorIQhF5QkTaRCRClyRYR7gxJgoVFfojNtwWKlbTmExgNtq7gcNUtZmqHgocD0wFRojIZRGI0RUlzVMFeVbTMMZED39hUcT6M6Bi04icqqp/+AZV1e3AR8BHIuKrsshclpSSCFjSMMZEF39RccT6M6ACNY1QCaOEiNQ9UJlok+gkjfzcApcjMcaY8BW5XdMQkaOBs4BnAT+BNbo7Br06AalAqBX3otb+moYlDWNMFAn0abjbPPUygenK1wJ7gIUEOsEXE1jXoouqbo5YhC6xmoYxJhr5/f6IzXAL4SWNKcDtwCwCNYpXVfVDABG5PRYTBlhNwxgTnfyF/ojNcAthJA1VHSYiqaqaIyL1gftE5GbgQcpZpzvaJVlNwxgThTwx5FZVc5yf21X1FgLNUpcAjUTkpIhF5yJfkg8RIT833+1QjDEmbJEeclupOoyqrlHVy4HjgLtE5IeqDct9IkJiss+ap4wxUSXSQ27DrsOIiKjq75qjVHUO0E9E+pRVJpolpiRa85QxJqpEeshthZ4IF5G/iUjz4J0ikgjEichbwKDyTiAio0Vks4gsCNr3hIgsEZF5IjK+5JmPEJ9dLSLzRWSOiGRVIO5KS05NspqGMSaqRHrIbUWSRj8Cz2mMEZENIrJIRFYCywj0cTytqm8e4BxvOucJ9hXQSVWPAn4hME1JWfqoahdV7V6BuCstMSWR/DxLGsaY6OEvcn/ILQCqmge8ALzgTBfSEMhV1Z0VOMcPIpJZat+koM2pwIXhni/SklISyc+xpGGMiR7+Qj++Wt6YsHA/VS1U1Y0VSRhhGgxMLOuywCQRmSkiQ8o7iYgMEZEsEcnasmVLpYOxPg1jTLTxxJDb6iAi9wJFwHtlFDlOVbsBZwBDReSEss6lqq+oandV7Z6enl7pmJJSEq1PwxgTVbzUER4xIjKIwPxWl5Y1+kpVNzg/NwPjgZ6RjstqGsaYaFPslVluI0VE+gF3AueUPEQYokwtEUkreQ+cDiwIVbYqWU3DGBNtPFPTEJHLRWSLiGQ7NQNE5FgReVhEZoZ5jjHA/4C2znmuAkYBacBXznDal5yyTURkgvPRRsBPIjIXmA58rqpfhP1bVlKS1TSMMVHGC7PclvgH0B9YBdwgIl8B7YAxwE3hnEBVB4bY/XoZZTc410NVVwKdKxBrlUhMtqRhjIku/iI/CQkeGHIL7FXVGQAiMhzYBBwRgRFUnmHNU8aYaOP6LLdBDnOGui51XtmxnDAAklKtpmGMiS6RHnJbkTPfDxwFXAocCaSJyNfAbGC2qr4fgfhclVwrmcL8QvxF/oiORjDGmKri+nKvJVT1leBtEckgkESOJPDsRMwljdS0FABy9uSSVq+2y9EYY8yBRXrIbaXrMKqaDWQDEw5UNlqlpCUDkGtJwxgTJTwz5LYm+q2mkedyJMYYc2Cq6qlZbmuclKDmKWOM8bpifzFARGe5rXDSEJGLgp7Ovk9ExolIt6oPzX2pQc1TxhjjdUWFRQCem0bk76q6R0R6A32Bt4AXqzYsb0ix5iljTBQpzA8kjaSUxIhdozJJw+/8PBN4UVU/ASIXoYtK+jSspmGMiQYFzqJxvmRfxK5RmaSxXkReBi4GJohIUiXP43klo6esT8MYEw0K8goBSPRY0rgY+BLo5zwRXh+4vUqj8ojfahrWPGWM8b7fkkbkGn8q3MXuTF8+Lmh7I7CxKoPyCl+Sj/iEeKtpGGOiQknzVCRrGhVOGiLyj1D7VfXBgw/HW0SE1LRk69MwxkSF6mieqsxg3n1B75MJrLi3uGrC8Z6UtBRy9lrSMMZ4X8ms3F5rnnoyeFtERgKfVllEHpOalmJ9GsaYqFBS0/Da6KnSUoHDq+A8npSSlmx9GsaYqODVPo35gDqb8UA68FBVBuUlKWkp1qdhjIkKnhw9RaAPoyRpFAGbVNVfTvmolpqWwvYNO9wOwxhjDshTz2mIyE/O2wVBryXADhHZHeY5RovIZhFZELSvvoh8JSLLnJ/1yvhsPxFZKiLLReSucOM+WNY8ZYyJFoXV0DwVdtJQ1d7OzzRVrVP6FeZp3gT6ldp3F/CNqrYBvnG2f0dE4oHnCSz21AEYKCIdwo39YKTWtuYpY0x0qI7mqWqd/kNVfwC2l9p9LoFJD3F+nhfioz2B5aq6UlULgLHO5yIuJS2FnD15qOqBCxtjjIs89ZyGiNxS3nFVfaqSMTRynipHVTeKyKEhyjQF1gVtZwPHlHVCERkCDAFo3rx5JcMKSE1LwV/kpzC/MKLZ2xhjDtb+CQuTPNA8BaQ5r+7AdQS+yJsC1xJoMookCbGvzH/6q+orqtpdVbunp6cf1IVt0kJjTLQoyCskMdmHSKivzKoRdk1DVYcDiMgkoJuq7nG2HwD+fRAxbBKRxk4tozGwOUSZbKBZ0HYGsOEgrhm24EkL66YfUh2XNMaYSinIK4h4i0hl+jSaAwVB2wVA5kHE8CkwyHk/CPgkRJkZQBsRaSkiicAAqukpdFvy1RgTLQqdmkYkVeY5jXeA6SIynkAT0fnA2+F8UETGACcBDUUkG7gfGAF8KCJXAWuBi5yyTYDXVLW/qhaJyA0EpmSPB0ar6sJKxF5htuSrMSZaFOR7MGmo6iMiMhE43tl1parODvOzA8s4dEqIshuA/kHbE4AJFQz3oNmSr8aYaFGQV4Avws1TlalpAKxyPpsMpInICc5w2phjS74aY6JFgRebp0TkauBGAp3Rc4Bjgf8BJ1dtaN5Q0jy1b7clDWOMt1VH0qhMR/iNQA9gjar2AboCW6o0Kg+p26guCb54Nq741e1QjDGmXPk5+SSleG/0VJ6q5gGISJKqLgHaVm1Y3pGY5KPlkc1ZmrXC7VCMMaZcObtzSa2TGtFrVCZpZItIXeBj4CsR+YRqembCLUd0b80vWSsoLi52OxRjjClTzu4cUuukRPQaFUoaEnjMcJiq7lTVB4C/A68Ter6omNGpdzv27cph9D3vux2KMcaUad/u3P2DdyKlQklDA7P2fRy0/b2qfupMIhiz+gw8jv5Xn8IHj3/ClE9muB2OMcb8gap6tnlqqoj0qPJIPCw+Pp4bRl1Fqy6ZPD3kJXZs2ul2SMYY8zsFeQX4i/zU8lLzlKMP8D8RWSEi80RkvojMq+rAvMaX6OOud4axb3cuTw15yaZKN8Z4So7zWIAXaxpnAK0IPJdxNoHlX8+uyqC8KrNjM67+56VM/WwmX4z+1u1wjDFmv337k4bHahqquibUKxLBedF5w86gy8mdeOGmN9hgz24YYzwiZ3cO4MGkUdPFxcVx+xtDiU+I57FBo/D7/W6HZIwx+5unanmwearGO7RZQ/426moWTVnKh49XywztxhhTrhyvNk+ZgJMv6c2JF/fi7Qc+YPnsVW6HY4yp4UrW/EnxynMaInJH0PuLSh17tCqDigYiwrAX/o86Desw4vJn96/Na4wxbvitecojSYPAankl7i51rF8VxBJ16tRP47bR17NmUbY9LW6McVVJR7hnahqAlPE+1HaN0aNvF865vi8fPfM5s76Z73Y4xpgaKi8nHxHx1Cy3Wsb7UNs1yv89fjkZRzRm5JXPs3fnPrfDMcbUQHn78klKTSQwRWDkVCRpdBaR3SKyBzjKeV+yfWSE4osKyalJ3PXOMLZt3MGoYa+7HY4xpgbKz8knOTUp4tcJO2moaryq1lHVNFVNcN7XUdU04I4DnqAcItJWROYEvXaLyE2lypwkIruCyvzjYK5Z1dr2aM1lf7+Qb979ke8/nOJ2OMaYGiY/t4AkLyWNA7j5YD6sqktVtYuqdgGOBnKA8SGK/lhSTlUfPJhrRsIl91xAu56t+dd1r7B1/Ta3wzHG1CB5OYHmqUirqqRRlY1opwAronFqkviEeO58+28U5hcx8qoXbVJDY0y1yc/JJ7lWcsSvU1VJoyq/HQcAY8o41ktE5orIRBHpWNYJRGSIiGSJSNaWLdW7fHnGEU0YMvIKZk6ay6cvfFmt1zbG1Fz5OQXeqmmIyJ6gzu/dpTrCm1RFMCKSCJwD/DvE4VlAC1XtDDxH0GJQpanqK6raXVW7p6enV0VoFXLWNafR44yuvHrHO6xdsr7ar2+MqXny9uV5riM8LajzO/iVpqoJVRTPGcAsVd0U4vq7VXWv834C4BORhlV03SolItz62nUkpSbx2BXPUVRY5HZIxpgYF6hpeChplBCRZBHpJCIdRaSqG9AGUkbTlIgc5qxRjoj0JBC7Z3ubGzSux00vX8MvWSt47+GP3A7HGBPjPNcRLiIJIvI4kA28DbwLrBORx0XEd7CBiEgqcBowLmjftSJyrbN5IbBAROYCzwID1OM9zcdfcAynDTqR9x8dx+Jpy9wOxxgTwwLPaXirI/wJoD7QUlW7qWpXAiv41QVGHmwgqpqjqg1UdVfQvpdU9SXn/ShV7aiqnVX1WFWNiochhj5zJekZDXjsiufI3ZfndjjGmBiVn1NAspdqGgSWdf0/Vd1TskNVdwPXAf2rOrBYUeuQWtz+5lA2LP+VV29/x+1wjDExSFXJ25fnuT4NDdUcpKp+avjcUwfS+cSOXHTr2Xz20iSmT5ztdjjGmBhTWFBEcbF6LmksEpErSu8UkcuAJVUXUmwa9NAAWh7ZnCeveoFdW3e7HY4xJobk5+QDeGvILTAUGCoi34nIkyIyUkS+B4YRaKIy5UhM8nHXO8PYs30vz1z7ij0tboypMvuTRi0PJQ1VXa+qxwAPAquBNcBwVe2pqvYEWxgOP6oFf31oAD+Nm8akt75zOxxjTIzIywmsHFodzVNhP5QnIs/xx76LNiJyHoCqDqvKwGLVn285i2kTZvH8sNF06t2Opq0bux2SMSbKrV+2EYDadVMjfq2KNE9lATOd1zlB70teJgzx8YFJDRN88Yy47Fl7WtwYc1CKi4t5/5GPaNi0Pl1PPSri16tI89RbJS9gR/C2s8+E6dBmDbnp5WtYMn057wwPNc2WMcaE56On/sui//3CXx8aQGLSQT9nfUCVneXWenEP0gkX9qLf4JMZ88/xzP1+odvhGGOi0Iq5qxl97/scd35PTh90UrVcs6qmRjeVcP0zf6VJ68N47PLn2LNjr9vhGGOiSH5uPiMue5Y6Detw88vXRHxt8BKVmhqdUmuEO/tMBaXUTuHu925k+687efqal20YrjEmbK/d+R6rF67jttHXc0jDOtV23cpOjf67NcJVtfoijjFtu7firw8N4Mf/TOXLNya7HY4xJgpM+WQGH4+ayPnD+tOjb5dqvbY1T3nAxbefQ5eTO/H8jaPJdobOGWNMKJvXbmHk4Odp060lVz92WbVf35KGB8TFxXHnWzfgS/Lx6CXPUJBf6HZIxhgPKios4pFL/oW/qJh7x95cLaOlSrOk4RENmzbg1teuY9nMlbx257tuh2OM8aC37v+QRVOWctPL17j2YLAlDQ857ryeXHDjmYx/dgI/jpvmdjjGGA/JmjSXsSPGc8ZVp9BnwHGuxWFJw2OufuxS2vZoxZNXvcDGlX9YKt0YUwNt27iDx654jsyOzbj+X1e6GoslDY/xJfq4d+zNiAgPD3ja+jeMqeH8fj8jLn+W3D253PfBzdUy/Xl5LGl4UOOWjbht9PX8krXCVvszpoYb8+h45ny7gBueu4oWHZq5HY53koaIrBaR+SIyR0SyQhwXEXlWRJaLyDwR6eZGnNXluPN68uebzuTjURP58aOpbodjjHHBvB8W8c7wDzn5kt70vbKP2+EAHkoajj6q2kVVu4c4dgbQxnkNAV6s1shccNWIS2nXszUjrX/DmBpnx6adPHrJMxx2eCNufHFItU0TciBeSxrlORd4WwOmAnVFJKYXoyjp34iLi+Ohvzxl/RvG1BD+Ij+PDHyGPdv38vcPbyE1LcXtkPbzUtJQYJKIzBSRISGONwXWBW1nO/v+QESGiEiWiGRt2bIlAqFWn8MyD+X2N4aybOZKXrntbbfDMcZUg9fvfo+53y3kppeuoXWXlm6H8zteShrHqWo3As1QQ0XkhFLHQ9XNQs7wp6qvqGp3Ve2enp5e1XFWuz+d24M/33wWnzz/BT/8539uh2OMiaDv//0//v3kZ5x9XV9Ou+JEt8P5A88kDVXd4PzcDIwHepYqkg0EDx3IADZUT3Tuu+qfl9DumDY8efWLrF9u81MZE4vWLM5m5ODnaX9sG657epDb4YTkiaQhIrVEJK3kPXA6sKBUsU+BK5xRVMcCu1S1xnx7+hJ93Df2ZuLj43jwwifJ3ZfndkjGmCq0b3cOwy94guRayfzj37fiS6z+eaXC4YmkATQCfhKRucB04HNV/UJErhWRa50yE4CVwHLgVeB6d0J1T6MW6dz93o2smr+Wp4e8ZOtvGBMjVJWRg19g/fJfuW/szTRs2sDtkMqU4HYAAKq6EugcYv9LQe8VGFqdcXlRj35dufLhgYy+932OOLoVF95yttshGWMO0odPfMpP46Yx5Ikr6HxSR7fDKZdXahqmAgbcdR7H//kYXr3jHWZ/O9/tcIwxB2HWN/MZfc97nHhxLy685Sy3wzkgSxpRSES4bfRQmrVrysN/eZpNa6J7WLExNdXmdVt5dODTZLRtwq2vXeeZB/jKY0kjSqWmpfDA+DvwF/l54IInyM/NdzskY0wFFOQX8tBFT1KYX8QD424npbZ3HuArjyWNKJbRpjF3vzuMFXNW8/Q1L1vHuDFR5IUb32DJ9OXc/uZQmrUN+ZyyJ1nSiHLHnHk0VzxwMd+8+yMfPP6J2+EYY8Lw2Ytf8vkrX/GXO86l9/nHuB1OhXhi9JQ5OJfcewFrl2Tz+t3vkZ7RgFMuPd7tkIwxZZg2YRaj/vY6x5zZjSsfHuh2OBVmSSMGxMXFcdvooWzfuJORg5+n3mF16XbKkW6HZYwpZfnsVTz8l6do1SWTe8fcRHxCvNshVZg1T8WIxCQfD4y7nYy2TRj+5ydYOW+N2yEZY4JsXreVe8/6J2n1a/PQZ3dHTcd3aZY0YkjturV45PN7SE1L4Z7+j7B5rQ3FNcYL9u3ax71nPkrevjwe+fweGjSu53ZIlWZJI8Yc2qwhj064h9y9edzT/1H27NjrdkjG1GgFeQU88OeRrFuygfv/cxstOzV3O6SDYkkjBrU8sgXDx9/B+mUbeeCCJ2zxJmNc4i/y8+il/2LOtwu4bfT1dDv1KLdDOmiWNGJUlz6duP3NG5j3/SIeH/QcxcXFbodkTI2iqjxzzcv8PH461z9zJadeVnqJoOhko6di2MkDe7M1exuv3vkudeqn8bfnr46KaQqMiXaqymt3vssXb0zmsr9fyPnD+rsdUpWxpBHjLrrtHHZt3cOHTwQe/Lth1FXExVkF05hIUVXeuG8MH478lHOu78sVD1zsdkhVypJGjBMRrh5xKSLwweOfUFysDHvhakscxkSAqvLK7e/wn6c+o//VpzD02cExV7u3pFEDiAhX/fNSJC6OsSPGU+wv5qaXh1jiMKYKqSov3PgGH4+ayDnX92Xos4Nj8v8xSxo1hIgw+JGBxMUJ7z86Dr/fzy2vXBuVT6Qa4zXFxcU8e92rfP7q1/z55rO4ZuQVMVfDKGFJowYREf760AASfAm8PfxDdm7exX1jb47aJ1ON8YKCvAKeGPwC3439mQF3nc/gRwbGbMIAG3Jb44gIl99/ETe+OISsL+Zw28nD2bF5l9thGROVdm7Zxe2nPsh3Y3/mqn9eGvMJAzySNESkmYhMFpHFIrJQRG4MUeYkEdklInOc1z/ciDVWnHXNadw/7nbWLFzHTcfdy/rlG90OyZiokr1sI8N63cvyWSv5+4e3MODO82I+YYBHkgZQBNyqqu2BY4GhItIhRLkfVbWL83qwekOMPX86pwePf3M/e3fmcNNx97F8ziq3QzImKqxdsp5bT7qf3D25PPndcE64sJfbIVUbTyQNVd2oqrOc93uAxUD0LGUVxTocewT/+vlhfMk+7jhlOEtnLHc7JGM8bc2iddzW5360uJiRkx+gXc82bodUrTyRNIKJSCbQFZgW4nAvEZkrIhNFpGM55xgiIlkikrVli830eiAZRzThqe8fpFbdWtxx2oMsnLLU7ZCM8aRVC9ZyW58HkLg4Rk4eTosOzdwOqdp5KmmISG3gI+AmVd1d6vAsoIWqdgaeAz4u6zyq+oqqdlfV7unp6ZELOIYclnkoT343nHqN6nJX34eY+/1Ct0MyxlNWzV/DHacMJ94Xz5OTH6B5u5rZGOKZpCEiPgIJ4z1VHVf6uKruVtW9zvsJgE9EGlZzmDHt0GYNefK74RzavCH39n+UmV/NdTskYzxhzuQF3HzCP0hITGDk5OFkHNHE7ZBc44mkIYEhB68Di1X1qTLKHOaUQ0R6Eoh9W/VFWTM0aFyPkZOH07RNY/5+9ggmvfWd2yEZ46pv3/+Ru/s9TIMm9Xjmp4fJaNPY7ZBc5YmkARwHXA6cHDSktr+IXCsi1zplLgQWiMhc4FlggKqqWwHHsiZ4iIMAABARSURBVHqHHsIT395Px97teOLK53n+xtEUFRa5HZYx1crv9/PW/R/wz8uepcOf2vLMTw/TqIU1d0usf+92795ds7Ky3A4jKvmL/Lx6xzt89MznHHlCe+5+90bSMxq4HZYxEbd1w3ZGXPYsc79byGmDTuSml64hMcnndljVRkRmqmr3UMe8UtMwHhSfEM+1T/2VO9/+G8tmrmTIUbfy/YdT3A7LmIia98Miru1yG0unL+f2N4Zyxxs31KiEcSCWNMwBnXrZCbw0+wky2jbh4QFP89ig59i7c5/bYRlTpXZt3c1Lt7zJnac9SJ2GdRg1YwSnDzrJ7bA8x5qnTNj8RX7ee/gj3nvkIxKTfZx62QmcM7QfLTs1dzs0Yypt87qtjP/XBP770iQK8go4fdBJ/N/jl1OnQZrbobmmvOYpSxqmwlbMXc24f33Od2N/piCvkE6929Hr7O50O/UoDu/cwpNrCKgqcyYv4MePplGnfm1adcmkc5+O1KkfmS+GgvxCsr6YQ3KtJNr2aEWtQ2pF5DrlUVX2bN/L8jmr2bcrhx79upCcmhSRaxXkF/LLjOXM/3EJ2zZsp9c5gb8HL87FpKqsXrCWKZ9k8fMn01k2cyVxcUKfS3oz8O4LaNE+w+0QXWdJw5JGROzetocv35jMpLe/Y/WCdQAc0jCNdse0ofHhjUjPaECtQ1Jp3bUlrbu2rPa1O3Zt3c2yWauY/8Mipnwyg9UL15GcmkRBXgHFxUpcnNC8QwaHZR5KUmoi/qJiGjSuR4c/teWEC48lwVexlQOKi4tZMm0Z377/E9998DO7tu7Zf6x+43q06JBB99M7c8xZR9O8XdMDfqEWFRaxZNoy/P5iOvQ6Al/igdvVt/+6I/Bl+PE05v+wmPzcgv3HUuukcPLA3px6+Ym069m6wv89igqL+Hn8dBb8tIQt67fhS/KRty+PzWu3sm7xegoLAiPsfIkJFBYU0bx9U447rydHndiRDr2OIDWteqfg9/v9bFyxidUL17FxxSa2ZG/j19WbWT57FVvWBUbrtz+2DX86tycnXtSLxoc3qtb4vMyShiWNiNu2cQezvp7HrK/nsXLuGjau3ETu3rz9xxs0qcefbz6bs645NSLrd/y6ejNZX85l+ayVrF2ynrWLs/d/acfFx9G2RyvOHHIafQYch8QJS2esYNZX81g2eyWb126lMK8QiRO2Zm8nZ08ujVqkc+EtZ9N3cB9SaiWXed3sZRv56aOpzPx6HstmrmTfrhwSk30ce3Z3+g0+GRFh2cyVZC/bwPJZq1g5bw0ATVo14pgzj6ZHvy5kHNGEhMQE8nML2LZhO6vmrWXOdwuY/c18cnbnApBcK4kufTpx9Omdad0lkwZN6uNLSqAwv4hfV29m3veLmPX1PBZPXYaq0qRVI3r270ajFum06NiM+IR4Jr01mR//M5WCvEJS01Lo2LsdHXu15Zgzu9GqS2aZSSx3Xx4TX/uG/zz1GVvWbSOldjKNWqRTWFBEcq0k6jeuR8uOzWjf6wiOPL49KWkpfP/hFD5/5SsWT11Gsb+YuPg4Wh7ZnBYdMsg4ognN2jahefsMMjs1q/KaaX5uPl+/8wNjH/uYX1dt3r+/JO7MTs3o0qcTx57dnQaN61XptWOFJQ1LGtVOVcndm8fubXtYPHUZE177mjnfLiA1LYXjzu9Jn4G96XbKkRX+166/yE9+bgHbN+5g+exVLPrfL2R9OYd1SzcAUKdBGs3bN6VZ26Y0b9+Ulkc2p/2x4f8rV1WZPmEW7/9zPIumLCWtXi3aHdOGuo0OoVadVFAozC9k4+rNrF2Uzdb12wFo3bUl7Xq2puNx7eh1TvdA2RA2r9vKtP/OZOrnM5n9zQIK8wtDlmvUIp2upxxJzzO6Ep8QT9aXc8iaNJeNKzeFLB8XJxzRvRU9+3ej9/k9yezUPGQS2LtzHzMnzWX2N/NZOGUpqxcGaohp9WuT2akZTQ4/jMTkQI1mz4697NyymxVzVrNn+16OPKE9F992Lj37dw37iz53by6Lpy5j3g+LWDJ9OdlLN7B57VZKvnfqN64XqI2c0J7DO2eSnlGfxJRE4uMr+Hfh9zP3u0VMHvMTP42bxt6d+ziieyvOvvZ0Wh7ZnIwjGrvSRBitLGlY0vCExdOWMeGVr/hx3DT27cohrV4tjj69Mz36daVHvy7Ua1R3f9nd2/awfPYq5kxewOxvF7B57Vb27thLQd7vv2QTk310Pqkj3ft2oecZXWnapnGVtaMvnLKUz178knVLN7Bj005yducSFyfE+xJIb9aAFh0yaNezDT37d6Vxy4o3beTuy2PJtGVsWrOV4iI/vmQfdQ89hBYdMji0WegZcjat2cKaRdns2LQTf6GfuIR40jPq07ZHa2rXrfiX4o5NO5k+cTaLpixl9aJsNq3ejL/QT3Gxkla/Noek1+GwzHTOub4fHf/UtsLnDyU/N5/1y35lxZzVTPl0BjMmzv5dMxoE/rumpqWQ3rwhnU/sSJeTO9HxT21/9zvu253DvO8DTY9TP8ti55bdpNROpvcFx9D3yj4cdUIHT/apRANLGpY0PKUgv5AZE2cz5dMZZH0xh+2/7gSgbnodatWtRe6e3P374uLjaHdMG1q0z6B23VRS0lJITk2iTsM0WnXOpHmHDBtDH+UK8gtZuziblXPXsHPzLvJzCsjdm0vOnjyyf9nAoilL9/eXpNWvjS/JR3GRn51bAnOaptZJcWpYx3DsWd1ISolMZ39NYknDkoZnFRcXs3LuGmZ+NY+NK35l355cklMSadauKYd3zqT9sW3KbOoxNUN+bj4Lf17Kslmr+HXVJor9xYgIh7U8lDZHH85RJ3YIa5CACZ8lDUsaxhgTNptGxBhjTJWwpGGMMSZsljSMMcaEzZKGMcaYsFnSMMYYEzZLGsYYY8JmScMYY0zYLGkYY4wJW8w/3CciW4A1bscBNAS2uh2ER9m9KZ/dn/LZ/SlfZe5PC1VND3Ug5pOGV4hIVllPWNZ0dm/KZ/enfHZ/ylfV98eap4wxxoTNkoYxxpiwWdKoPq+4HYCH2b0pn92f8tn9KV+V3h/r0zDGGBM2q2kYY4wJmyUNY4wxYbOkEUEi8oSILBGReSIyXkTqBh27W0SWi8hSEenrZpxuEZGLRGShiBSLSPdSx2r8/QEQkX7OPVguIne5HY/bRGS0iGwWkQVB++qLyFcissz5Wc/NGN0iIs1EZLKILHb+v7rR2V+l98eSRmR9BXRS1aOAX4C7AUSkAzAA6Aj0A14QkXjXonTPAuAC4IfgnXZ/Apzf+XngDKADMNC5NzXZmwT+JoLdBXyjqm2Ab5ztmqgIuFVV2wPHAkOdv5cqvT+WNCJIVSepapGzORXIcN6fC4xV1XxVXQUsB3q6EaObVHWxqi4NccjuT0BPYLmqrlTVAmAsgXtTY6nqD8D2UrvPBd5y3r8FnFetQXmEqm5U1VnO+z3AYqApVXx/LGlUn8HAROd9U2Bd0LFsZ58JsPsTYPchPI1UdSMEvjiBQ12Ox3Uikgl0BaZRxfcn4WCDq+lE5GvgsBCH7lXVT5wy9xKoOr5X8rEQ5WNy7HM49yfUx0Lsi8n7cwB2H0yFiUht4CPgJlXdLRLqz6jyLGkcJFU9tbzjIjIIOAs4RX97KCYbaBZULAPYEJkI3XWg+1OGGnN/DsDuQ3g2iUhjVd0oIo2BzW4H5BYR8RFIGO+p6jhnd5XeH2ueiiAR6QfcCZyjqjlBhz4FBohIkoi0BNoA092I0aPs/gTMANqISEsRSSQwOOBTl2Pyok+BQc77QUBZNdiYJoEqxevAYlV9KuhQld4feyI8gkRkOZAEbHN2TVXVa51j9xLo5ygiUI2cGPossUtEzgeeA9KBncAcVe3rHKvx9wdARPoDzwDxwGhVfcTlkFwlImOAkwhM970JuB/4GPgQaA6sBS5S1dKd5TFPRHoDPwLzgWJn9z0E+jWq7P5Y0jDGGBM2a54yxhgTNksaxhhjwmZJwxhjTNgsaRhjjAmbJQ1jjDFhs6RhjDEmbJY0jDHGhM2ShjEmbCLynIjMEpEebsdi3GFJwxgTFhGpRWCG1GsIzKdmaiBLGsZUMxF5QERuC9qeUsnz1BWR6yvxuYEiMkdE5ouIisjsEGVSROT74MWvVHUf0Bj4DnhWRBJF5AcRsYlPaxBLGibmSYArf+vhXFtV/1TJ09cFKpw0VHUMcBGwh8AEd2eGKDYYGKeq/pIdItIASHU+53cWhvoG+EvFQzfRypKGiQoikumst/6Ws+b6f0Qk1Tl2i4gscF43BZVfLCIvALOAZiJymYhMd/6V/XKoJWRF5Arn/HNF5J2g/X+4RgWvfa+z1vfXQNtS19wb9JlXnfWdJ4lISlCZj0VkpnNsiLN7BNDK+X2ecMqF8zt2JDDz6YOqerWqhppu/VL+OBvqfcBIYCGB5WchMFngpSE+b2KVqtrLXp5/AZkEFiA6ztkeDdwGHE1gVs9aQG0CX2hdnfLFwLFO+fbAZ4DP2X4BuKLUNToCS4GGznZ952dZ1wj32iXlUoE6BJavvS3ounudzxQBXZx9HwKXBZUpiSWFwNrqDZzPLAgqE87v6ANmAseXc68TgV9D3P9pBBaGGgX8n7M/Htji9t+HvarvZW2RJpqsU9WfnffvAsOAQmC8BtrbEZFxwPEE/iW9RlWnOuVPIfDlPcNZySyFPy5GczLwH1XdCqC/TR/du4xrSJjXPt4pl+OUK2tNjFWqOsd5P5PAF3WJYc5U8hBYmKkN8Gupz4fzO/YDFqrqj2XEAIFpx3eW2vcwgZqJishiAgkWVfWLSIGIpGlgXWoT4yxpmGhSeh5/JfSSqCX2Bb0X4C1Vvbuc8hLiGiX7yyofzrUp47yl5Qe99xP40kdETgJOBXqpao6IfAcklxHPgX7HnsDkA8SRG3x+EekCXAD0FpHnnWPzg8onAXkHOKeJEdanYaJJcxHp5bwfCPwE/ACcJyKpzpDQ8wksRFPaN8CFInIogIjUF5EWIcpc7HT4IiL1nf1lXSPca/8AnO+MSEoDzq7g730IsMNJGO2AY539e4C0Cv6Oe4BelENVdwDxIlKSOB4DzlbVTFXNBDrj1DSce7VFVQsr+DuZKGU1DRNNFgODRORlYBnwovNF+ia/LQf7mqrOFpHM4A+q6iIRuQ+Y5IxmKgSGAmuCyiwUkUeA70XED8wG/qqqs0JdAyDMa88SkQ+AOc71ymsaCuUL4FoRmUegz2Wqc95tIvKziCwAJqrq7Qf6HYFXgdFOE9MK4FJV3RXimpMI1CyKgVqq+k3Q77NJRGo5SbUPMKGCv4+JYrZyn4kKzhfxf1W1k8uhxAwR+ZxAP8W0EMe6Areo6uUHOMc44G5VXRqhMI3HWPOUMTWQiFxKYO36GaGOOzWpyaGG7AadIxH42BJGzWI1DWOMMWGzmoYxxpiwWdIwxhgTNksaxhhjwmZJwxhjTNgsaRhjjAmbJQ1jjDFhs6RhjDEmbP8PY6Wxj264wIEAAAAASUVORK5CYII=\n",
"text/plain": [
"