{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# All distances between two selections\n",
"\n",
"Here we use ``distances.distance_array`` to quantify the distances between each atom of a target set to each atom in a reference set, and show how we can extend that to calculating the distances between the centers-of-mass of residues.\n",
"\n",
"**Last executed:** Feb 06, 2020 with MDAnalysis 0.20.2-dev0\n",
"\n",
"**Last updated:** January 2020\n",
"\n",
"**Minimum version of MDAnalysis:** 0.19.0\n",
"\n",
"**Packages required:**\n",
" \n",
"* MDAnalysis (Michaud-Agrawal *et al.*, 2011, Gowers *et al.*, 2016)\n",
"* MDAnalysisTests\n",
" \n",
"**Optional packages for visualisation:**\n",
"\n",
"* [matplotlib](https://matplotlib.org)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import MDAnalysis as mda\n",
"from MDAnalysis.tests.datafiles import PDB_small\n",
"from MDAnalysis.analysis import distances\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Loading files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The test files we will be working with here feature adenylate kinase (AdK), a phosophotransferase enzyme. (Beckstein *et al.*, 2009) AdK has three domains: \n",
"\n",
" * CORE\n",
" * LID: an ATP-binding domain (residues 122-159)\n",
" * NMP: an AMP-binding domain (residues 30-59)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"u = mda.Universe(PDB_small) # open AdK"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculating atom-to-atom distances between non-matching coordinate arrays"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We select the alpha-carbon atoms of each domain."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LID has 38 residues and NMP has 30 residues\n"
]
}
],
"source": [
"LID_ca = u.select_atoms('name CA and resid 122-159')\n",
"NMP_ca = u.select_atoms('name CA and resid 30-59')\n",
"\n",
"n_LID = len(LID_ca)\n",
"n_NMP = len(NMP_ca)\n",
"print('LID has {} residues and NMP has {} residues'.format(n_LID, n_NMP))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`distances.distance_array`([API docs](https://docs.mdanalysis.org/stable/documentation_pages/analysis/distances.html#MDAnalysis.analysis.distances.distance_array)) will produce an array with shape `(n, m)` distances if there are `n` positions in the reference array and `m` positions in the other configuration. If you want to calculate distances following the minimum image convention, you *must* pass the universe dimensions into the ``box`` keyword."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(38, 30)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dist_arr = distances.distance_array(LID_ca.positions, # reference\n",
" NMP_ca.positions, # configuration\n",
" box=u.dimensions)\n",
"dist_arr.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plotting distance as a heatmap"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Distance (Angstrom)')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAR8AAAEWCAYAAABfWJOFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO29eZhkVZWv/f4ics7KmguoogqQQRwQUVHsj3YAhwacZ7k4YKtc+2uv2iqNQ1/Fq7TatoL9dasXhwaVC4440A6gKDhc0AJKEEGBYqoqqKKGrCnHiFjfH+ckBHnWyYykIjIiMtdbz3kqY5199tlniBV777XXWjIzgiAIZptCsxsQBMH8JJRPEARNIZRPEARNIZRPEARNIZRPEARNIZRPEARNoa2Vj6QvSPqfzW7HI0HSsyVtaHY7moWkuyQ9t95lZxNJp0v6dbPb0a60rPJJX7hhSbslDUr6raS3SXqwzWb2NjP7aI11tdzL+0iZ74ormBu0rPJJeZGZDQAHA58AzgK+3NwmBQFI6mh2G9qdVlc+AJjZTjP7AfAa4I2SjgKQdIGkj6V/L5d0WdpL2i7pV5IKkr4GHAT8UNIeSf+Ylv+WpPsl7ZR0taTHT5wvrfc/JP1X2vO6VtJhVfsfL+mK9DybJX0glRckvU/SHZK2SfqmpKVTXZukD0jamvbOTquSd0v6V0n3pOf4gqReSf3Aj4FV6fXskbQq7SUuT4/9oKSSpIXp549KOm+qeqvO+0JJ66p6m0dX7btL0nsl3Zjet29I6sm5rsMkXZneh62SLpK0OKfs2ZK+nda3W9L1kp44qdgx3nklLUmf+wOSdqR/r57mnr9V0i3puf4k6cmpfOLZTchfVnXM6ZJ+I+lcSduAsx/apX9P23WrpOdUHbNK0g/S9+R2SW+ddM3flPTV9Hw3Szp2qnbPOcysJTfgLuC5jvwe4O/Svy8APpb+/XHgC0Bnuj0DUF5dwN8CA0A3cB6wrmrfBcA24GlAB3ARcEm6bwC4D3gP0JN+Pi7d907gGmB1Wu//Bi7Oub5nAyXgM2nZZwF7gSPT/ecCPwCWpuf4IfDxqmM3TKrvauAV6d+XA3cAJ1fte1kN9T4J2AIcBxSBN6b3rrvqPv4OWJUefwvwtpzrOxx4XnptK9I2nOc9X5Iv8jjwyvTZvRe4E+ic7rzAMuAVQF96Pd8CvjfFe/UqYCPwVEBpOw+u2reK5Ef5NenzWJnuOz19Xv8jfSd6q2T/kLb7NcBOYGnVff8cyXtyDPAAcGLVNY8Ap6T3+uPANc3+3s3qd7zZDZjiJXnw5Zwkvwb4YPr3BTykfP4X8H3g8Frrqtq/GDBgUVW9X6rafwpwa/r3qcANOfXcAjyn6vPK9EvV4ZR9dvri9lfJvgn8z/RLsRc4rGrfXwF3Vh07Wfl8FPi39ItxP4ki/ET64g+nX9Lp6v088NFJ9f4ZeFbVfXxd1b5/Ab5Q4/N8afV9I6t8rqnaVyBR8M+Y6XnTL/mOKdrxU+CdNbZ5HfCS9O/TgXsm7T8d2ET6I5fKfge8HlgDlIGBqn0fBy6ouuafVe17HDDcjO9as7a2GHZN4kBguyP/FHA7cLmk9ZLel1eBpKKkT6Rd7F0kLzfA8qpi91f9PQQsSP9eQ9Kr8DgYuDQdsgySKKMysH9O+R1mtrfq890kv7wrSH7Jr6uq6yepPI+rSJTSk4GbgCtIelNPB243s2011Hsw8J6Jfen+NWmbJsi7Lw9D0v6SLpG0Mb3HX+fh93cy9078YWYVYEMt55XUJ+l/S7o7Pc/VwOL0GT+jamh6c3ps7vOT9IaqIecgcNSkNt/rHLbRUu2RMvEMVwHbzWz3pH0HTnFNPZpHc0ltpXwkPZXk4WXMm2a228zeY2aHAi8G3l01/p7suv/fgJcAzwUWAYdMnKKGZtwLHDrFvpPNbHHV1mNmG3PKL0nncCY4iOSXdCtJb+XxVfUsMrOJL7oXiuC3wJHAy4CrzOxPaX2nkCgmaqj3XuCcSe3vM7OLp7spDv+ctvMJZrYQeB1T3981E38osWiuJrkX0/Eekus+Lj3PMyeqMbNfmdmCdJuY07sXOGxyJZIOBr4IvB1YZmaLgT9OarN33w+UVF1m4hluApZKGpi0L+9dmHe0hfKRtFDSC4FLgK+b2U1OmRdKOjx9EXaS9Dgq6e7NPFxhDACjJPM6fSRflFq5DFgp6V3p5O2ApOPSfV8AzklfZCStkPSSaer7iKQuSc8AXgh8K/3l/yJwrqT90roOlPQ3VdezTNKiiUrMbAi4Dvh7HlI2vwXeNvG5hnq/CLxN0nFK6Jf0gklfoFoZAPYAOyUdCJw5TfmnSHp5+sv/LpLnc02N5xkGBpVM7n94mvJfAt4r6SnpNR6ePq9+EuXyAICkN5H0fKZjP+AdkjolvQp4LPAjM7uX5P5/XFKPkon7N5P0AANaX/n8UNJukl+rD5JMzr4pp+wRwM9IXvj/C3zOzH6R7vs48E9pd/q9wFdJusAbgT9R20sOJD0skonUF5F0m28DTkh3f5ZkMvfytN3XkEze5nE/sIPkV/IikknUW9N9Z5EMI69JhxM/I/mFJy1zMbA+vaaJ4clVJBOfv6v6PEAyFKGGetcCbwX+PW3X7STzGo+Ej5AMAXcC/wV8d5ry3yeZsN1BMmfycjMbr+E855FM/m4lud8/maqwmX0LOAf4P8Bu4HskE8R/Aj5N8u5sBp4A/KaG819L8u5tTet9ZTrEhWR+8BCS53sp8GEz+1kNdc4LJqxBQdA0JJ1NYih4XbPbEswerd7zCYJgjhLKJwiCphDDriAImkL0fIIgaAoNW9Ak6SskpuMtZjbhi/UpEivRGMlCrzeZ2aCk55Gsxu1K951pZldO2/iefuvuz7pOmXNVKvt1WI76rXTllHdWqiin82gd/o5iZ05jHJSzMqa/Y8yV9xR8A5GcJSrd8svuqbjuWoxUOjOyvoLfjkrOkp6dY72u3CveU/TbN1LOtgOgNFZ05YWOSlZW8J9NadSvI/dnuuQ0vOjXPXb3xq1mNtVC0Wn5mxP6bdv22t6f624c/amZnbQv52skDRt2SXomidn7q1XK5/nAlWZWkvRJADM7S9KTgM1mtkmJ0+hPzezA3MpT+petsce/4F0Z+cjS7JvSvSPnhRjwvyR71/jlK45iy/m+M76s5MqXrNzlyuVosc5i9osDcOwKb7EtPLbfX5fX6Wjfw7o2u2V/vedIV37LngMysict9NsxlKO9f7zhca684Fz7o5dsccv+Zcd+rnzzxiWuvG/pUEY20Dvq13FPjh9wl/8cijuyL0R5ga8c7jnjrOvMbJ+cR499Yo/97qcH1VS2uPK2fT5fI2nYsMvMrmaSG4SZXW5mE9/ICQdMzOwGM5v41twM9ErqblTbgqBdMaBS479Wp5l+JH8LfMORvwK43szcnyZJZwBnAHT1+790QTBXMYxxq33Y3so0RflI+iCJR/dFk+SPBz4JPD/vWDM7HzgfkmFXA5sZBC1JO/RqamHWlY+k00kmop9T7Q2sJADUpcAbzCzPazwI5jWGUZ4jy2NmVflIOgn4R5L4MENV8sUk/j/vM7Na/GkAKJSMnm21dUF7t/rlOof8aa9yr2/18CacO7LzmQCMlPzbu6Ow0D/AI8eUtjan+Naxfle+pndHRlbxTHfAH3b6c/137liWkZUq/v0rVfz7t3Vbjo+qc5l57du2w43iQedW/34PlbP3ZGSBPyHeud1vd7nHv87OXVm5PAtYHam4zvXtR8MmnCVdTOKkd6SkDZLeTOKwOABckcZN+UJa/O0kEeU+lMrXTXhdB0HwEAaUsZq2VqdhPR8zO9URu8HfzexjwMca1ZYgmEvMlZ7PvImaFgRzAQPGY84nCILZxtpkSFULba18Kh1iZGl2gnBkSXYqq+h7AVDq9icHx/15W8yZj8xz3SjnuGgUev2Vz1bOtsUqfvu6iv5J+3MudHnnnozsgM6dbtlVvf4K7C1D2cniRZ0jbtk8unv95eCl8eyN7enw71NXty8f78r5UjqT9sWcleOOB0ki7/XLl5zJ5dJAA9fhGJTnhu5pb+UTBPONZIXz3CCUTxC0FaJcU56D1idCagRBG5FMOKumrRbSFEM3SLos/XyBpDurlrwc06hriZ5PELQRyTqfuvZ83kmSX6565euZZvbtep7EI3o+QdBmVEw1bdORujS9gCSd0KzT1j0fVaBzODv1X96blXXt8i0QHR3+QxpZ5t8aL/hYz/Y8K0tOrKA+3wxWGM+Wz1uqv6lrsSsfz3F3GHYCcHnBwQBu3+0nFt0y6Ls1eHTmWONGBv1AZRrNtntz0XfFGNvt37/ewRxXD+e+5g1Lenf68jF8t4visGOhLDTuN73OPZ/zSNydJt/ocyR9CPg5icuTH/xoH4meTxC0EYYoU6hpA5ZLWlu1nTFRT5qEc4uZXTfpFO8HHgM8FVhKkuetIbR1zycI5iO1DKlStk4RyfB44MWSTgF6gIWSvl6VO21U0n8C79231uYTPZ8gaCMMMWbFmrYp6zF7v5mtNrNDgNeShDd+naSVAGna8ZeS5KtvCNHzCYI2Illk2NA+w0WSVpCE9F8HvK1RJwrlEwRtRr0XGZrZL4Ffpn+fWNfKp6CtlU+lA4aWZ38FxhZnH07HiN8NzRs+j+XF+/IypYz6lZRyssRYjp9QpegEpsqxM3R2+v5NfZ2+79Ty7qxv18quQbfsgX2+z9fm3qz1aXHPsFt2Qaff8I0Lcqx0BSctT4/vp5bn1F3q96135e7sAcrxrystyAkmttAvX+nOPjPrb5xvl5ko5+V7ajPaWvkEwXwkLydauxHKJwjaiGTCeW58befGVQTBPGEWJpxnjVA+QdBmlGtf59PShPIJgjZiYoXzXKBhykfSV0jyc22pytX+KeBFwBhwB/AmMxuUdBpwZtXhRwNPNrN1U57D/AiFXiqbrl2+hSlv+Ny513/AXtTC7sEc61WXX8f4rtr9hPIsacMDvo/UfQXfFNRd9K01Hut3Z1PkAOwc7MvIxsv+tXTnWN1K2/12d+7O3qvBiu/bpWH/vg5s9u9VpTMrHx3xs3H33Zfj6zfmvyieP974SGOVQ2WOWLsaeRUXACdNkl0BHGVmRwN/IfEjwcwuMrNjzOwY4PXAndMpniCYjySOpTX7drU0jUydc7WkQybJLq/6eA3wSufQU4FLGtWuIGhnDDE+jetEu9DMOZ+/Bb7hyF8DvCTvoNQz9wyArv4ljWlZELQoZsyZRYZNuQpJHwRKwEWT5McBQ2aW68xmZueb2bFmdmxHT06KiSCYs4hKjVurM+s9H0mnk0xEP8css1D+tcDFtdZlBSg585eurNd/GDmp0N0UOXnk1ZGXUoecdDgFJ3BYISfljxd8C2A8Z2J0tJyV501cFnJvSrZ9lbxryXn3LWdCfCZYh19HZQZvc97zLfvz0JT9eXIqTrqeco77TD0w5k7PZ1aVj6STSCKnPcvMhibtKwCvBp4xm20KgnajHSaTa6GRpvaLgWeTRFPbAHyYxLrVDVyRhAvhGjObcNl/JnCvma1vVJuCoN0xaovP3A400tp1qiP+8hTlfwk8vVHtCYK5QJI6Z26sDZ4bVxEE84a5kzQwlE8QtBHG3Fnh3NbKpzAO/VuyloURxxLUu9Vf7l8Y8y0Towv9SGDecLtn+8yCR40v8F+eohN/q3OPb9nJq2Osw08rs7kn66rQ2+Hfky27c1Lk7Mq+LiM5ljHLmZco5riteK4l447rAkBhLMdKl2MZdEcpMzS6KceAJecWFnPcP+pF9HyCIJh1zBQ9nyAIZp9kwjncK4IgmHUihnMQBE0gmXCOOZ8gCJpArHBuAQpjFfrvzUYO6xjKOuJ0bfDTxOTlYelb6luNvOF21/YRt2ylw7eYde3KsdaMZ9vS4WemoTiSYwnKsbQM781ez+YuP1jXnsGcdjsBv8Y6/FdoZNxvR++gL/cCwJV6/bmNjpxr7x70n6Xnl5Xn69e9I89RL8d65zyfck7KpHoQK5yDIGgaEUA+CIJZxwzGK6F8giCYZZJhVyifIAiaQKxwDoJg1glTe4tgRTG+MGvFGVmStZJ07L/QrUNl37oxvNS3tFRcI5hv3hhZ7NcxstwVuxEO86xdo8t9ZyPt5ziIAcuX7M7IHr1ki1v2lpzmbR1flJH1LXPMVEBft+83tm3Uj7tdcHy7WOlfy8iQ/9pqvNOVl/uyz3hsud++Spdf99gi/353DDnRHZ3ohvUjhl1BEDSJdojPXAtzQ4UGwTwhsXYVa9pqQVJR0g2SLks/P0rStZJul/QNSf6CtzoQyicI2oiJRYa1bDXyTh4+0v4kcK6ZHQ7sAN5c50t4kFA+QdBm1Ct1jqTVwAuAL6WfBZwIfDstciHw0gZdxqznav8oSULACrAFON3MNqW52s8CBOwG/s7M/jDtOSpGcSibg7xjNHtZhVJOrvacZfN52WM8yk4ucICKP/9JJSf1i0cp5yWq5KRn6c2Z6O11cqd3F/wgaMq7eKfd3Z1+DvhiIe9+56S9ce5hR9FvX1k5Q4qcn1L3OeQ8g7wOQ14wsYJzu/PeqXowQ2vXcklrqz6fb2bnV30+jySbzISfzTJg0MwmHuoG4MB9aO6UzHau9k+Z2dFpTvbLgA+l8jtJ0uk8AfgocD5BELhUrFDTBmydSLCZbg9+ryRNdAyua9Z1zHau9l1VH/tJg1ma2W+r5NcAqxvVriBoZ8xEqT6m9uOBF0s6BegBFgKfBRZL6kh7P6uBjfU4mUczMpaeA7wB2Amc4BR5M/DjKY5/MFd7T1d23UkQzHXqscjQzN5PkkcPSc8G3mtmp0n6FvBK4BLgjcD38+qQtB+JElsFDAN/BNaaWU0pW2d9wtnMPmhma0jytL+9ep+kE0iUz1lTHP9grvbOzsjVHswvJuZ86mjtmsxZwLsl3U4yB5TJtSfpBEk/Bf4LOBlYCTwO+CfgJkkfkeSv6q2imYsMLwJ+RJLJFElHk8y6n2xm25rYriBoaertXpEm7Pxl+vd64GnTHHIK8FYzu2fyDkkdJIam5wHfmaqS2c7VfoSZ3ZZ+fAlwayo/CPgu8Hoz+0ut9VW6Cuxdk3VtGNov26Erd/a5deTF4t6z2n/Apd6slWQ0x41ibHFO2ps1vtuAZ2WqjPt173+AHxzt8Uvvd+WH9T2QkT2x7263bH+H374/dGUNH09Ysskt213wrWA/L/nXM17Oyh+7YrNbduMef7i9qbLMlXcuzF7PgYv3+HVrqV9Hn5+XZ7S/OyNTz8xSKc2EVggmZmZnTrGvBHyvlnpmO1f7KZKOJDG13w1M5Gn/EEkX73NpDveSmR3bqLYFQTvTKu4VkhaTzN8eQpUuMbN31HJ8S+RqN7O3AG9pVFuCYK5gBqXWCSb2IxLr9E0kHYoZEY6lQdBmNHvYVUWPmb37kR4cyicI2ohWmPOp4muS3kqyYPjBiTUz217LwaF8gqDNsNZRPmPAp4APki4YTv8/tJaD2175eM/BWwCa52dV7vIfZDlrxEjqcQIMlHxDGiUniBVAd6/vf1UsZofNFSfAGMCqBTtd+UG9/o/O6q7s6oUVxWyAMYAlXh4bYEl3Vn5A1y6nZD59Xf61G1n58u69btkdo/4Nz7My9fdmLVULu/10R/d1+XV0dPhTGuPO41GhkcHEWmfCGXgPcLiZbX0kB7e98gmC+YRZS8353A74v1Q1EMonCNoKUW4da9deYJ2kX/DwOZ/mmtqDIGgMLTTn8z1qXFDoEconCNqIVspeYWYXpmFWH52K/mxm/qSeQyifIGgnLJn3aQVSb/gLgbtIAgGukfRGM7u6luPbWvmUu2HXo7Lj3xEnrczYQn+cXO7O8b86yLeGdHRnfZZGhn1T2sASfy7uKQdscOVdjj9UXpqUZy++1ZU/tSfj6wfAXeOLM7JVRd+H69WL1rryg7uzRo1n9K53y+Ylthvycw+5v+bHD/hufn/oPrjmOgAOX5j1azu01zfQjJb9r4Rn6QO4uyfrC7a01y97lyudOS1k7fo08Hwz+zOApEcDFwNPqeXgtlY+QTDfsNaacO6cUDwAZvYXSTmLWrKE8gmCNqNVhl3AWklfAr6efj4N8LvNDqF8gqDNaCFr198Bfw9MmNZ/BXyu1oND+QRBG2HWGspHUhH4ipmdBnzmkdQxJ5WP616Rkz87z+2ikJNapduZcM6jv9sPQLW003cbWJATxMvjgA4/mNiqoh+sq2xZd4zVHQvcsveV/EBbB3Rk61hZ9CeQx/HdFPbv9N0xOpW9r4d2+K4i27r8dq/s8+s+vC+bk/6Ibj9Q2a0LDnDl/UX/WXquHgNdvrGiXrSCqd3MypIOltRlZv7NmYY5qXyCYC7TQnM+64HfSPoByWpnAMyspp5QKJ8gaCMMUWkda9cd6VbgocSDNavGUD5B0Ga0TseHP5nZt6oFkl5V68Eto0KDIKiBdMK5lm0WeH+NMpfo+QRBu9Hkro+kk0nS5xwo6d+qdi0EarbINDJ7xVdI8vdsMbOjUtlHSVLmVIAtwOlmtknSmSQLlCba9FhgxXThGFWCngecdDMd2Q5djoGJSs4d2NvrRxPbs8CxJo36HcgtOWlv/ti1ypV7VpJSxa+jU3npWW5zpatyAod5XD+23JX/YSjr1rC0+Ee3bMV8K9htw/u5cu961jgB0AA2jPkpch4Y9q1g93dnU+0MFHyL1M6xHlc+3uE/h92j2ffESwNUT1rA1L6JZDHhi4HqXO+7gX+otZJGDrsuAE6aJPuUmR1tZseQxH39EICZfcrMjknl7weuqjUObBDMJ4wkumUtW8PaYPYHM7uQJIrhhenfPwBuN7MdtdbTMOWTerZunySrXojRj9+BPJXEOS0IgskYSezgWrbGc4WkhZKWAtcDX5R0bq0Hz/qEs6RzJN1LMsz60KR9fSS9pdw0q5LOkLRW0trSSM5YKgjmMGa1bbPAorRD8XLgq2Z2HPCcWg+edeVjZh80szUkudrfPmn3i4DfTDXkMrPzzexYMzu2o6e/kU0NgtbEatwaT4eklcCrSaZRZkQzTe0XAa+YJHstMeQKgimozcw+S5PS/wv4Kclcz+8lHUqexcNhVk3tko4ws4nGvQS4tWrfIuBZwOtmVmmNspxkroWcoI95xiQrZysv5Fi7Kp3+7d0znhNQy2n48LjvfLapLxscDOCB7oWuvF9Z95uhiu+SM26+1WiPk09osOynsenJubGjOebFUedV3FvxLY7j5luT8gKB7Sr1ZmQj5t/XUk7dBfl+d53F7IviyepKi6wyTBcYfqvq83qyHYpcGmlqvxh4NrBc0gbgw8Apko4kUQV3A2+rOuRlwOVmFhM5QZCHgdXBkiWpB7ga6CbRA982sw9LuoCkEzDhRXy6ma3LqePfHPFOYK2ZfX+6NjRM+ZjZqY74y1OUv4DEPB8EwZTUZUg1CpxoZnvS6IO/lvTjdN+ZZvbtGuroAR7DQ72fVwB3Ak+UdIKZvWuqg2OFcxC0G3UYdpmZAROxUzrTbaY1Hw0cb2ZlAEmfJwko9tfATdMdHL5dQdBu1G7tWj6xLCXdzqiuRlJR0joSb4MrzOzadNc5km6UdK6knMThACwBqicI+4GlqTKaNjhV9HyCoJ2YWGRYG1vN7NjcqhIlcYykxcClko4i8TC4H+gCzgfOIrFqefwLScbSX5KMBZ8J/LOkfuBn0zWurZVPpROGDsg+iJH9staGcldOJy+v77efr7j7erMWouFu/8dhwULff+gxi7OR9QC6i1mfvD0l3zL2lAV3ufIn9/hpeQ7ryFp8dlV8H8AndN3vysf6s5agY7r9iIqDOdbFx/bf58o9C9ZhndmUNwAjOeEnb+n3oxDu352NcOhFZQRY1u3bO/br9n3j9oxnn/2SrmG3bL2o9wJCMxtMUx6fZGb/mopHJf0n8N4pjvuypB8BT0tFHzCzTenfZ0533hh2BUG7UVFt2xRIWpH2eJDUCzwPuDVdNIgkAS8FfM/hhygADwA7gMMlPbPWy2jrnk8QzEdUn57PSuDCNBB8AfimmV0m6UpJK0iGUet4+HKYh7dD+iTwGuBmHlpJZyQm/GkJ5RME7USdXCfM7EbgSY78xBlU81LgSDOrPfNBFaF8gqCtmDWP9VpYT2Kin4fKpwClBdmfAevLTjhXhmc2vdXZ5U/Gdnc6udoL/qRwIad/3Fv0XQ8Kys7S5uVqL89wum7Usu3+/t5D3LLbS757xZE9mzKybY67CcAPdj/RlS/v8Cdud1eyQbw6nfsB+a4bXq77PPKCseXVsajoTyIvc3K4L+1q8CL9FnGvAIZIrF0/p0oBmdk78g95iGmVj6QO4GSSlYwAtwA/MXPe5iAIGk+OJbEJ/CDdqqlP9gpJBwJXAvcBN5BMQr0Q+HS6fDr7UxgEQeOY2TqfhpJGMHwQSWtIIlPUxHQ9n3OAz5vZeZNO8g7g48Abaz1REAT1oU7WrrqQWsZeRRKBdBVwaa3HTqd8nm5mp08Wmtm/SfrzTBoZBEGdaH72igGS6IX/DXg08F3gUWa2eib1TKd8plqqmZ1pC4JgPrAF+B3wT8CvzcwkvWymlUynfBZJerkjF0mOnqZiXRXKB2ZdGNbsl13yf3+f39xKzvj5qJW+G8B+PXsysjv6/VQzq/r9JfwnL/mDKy86P2n3l7JpXwBO7FvvypfnWN5+PJRt48+2P84tOziWdcUAWLB/9l7fNuq7NPxq2+Gu/NAFW135KYtvzLZvj9++W/audOUvX369f87O7DlvHdvfLfuKpWtd+Xe2+y5SO5x7dfsu/32oFy0w7Ho/ydzO54CLJX3jkVQynfK5iiSuskdNqxiDIKgjxrSuEw1vQjIHfF4aNvW1wPeAVZLOAi41s7/UUs+UysfM3rTPLQ2CoL40v+cDPBg29Z9JPNmPIpl0/hHgd3snMZ2p/d3TnPwzNbYzCII60exhlySlwcgexMz+CHww3dwyk5lumezAFJu/DDYIgsbS/NQ5v5D0PyQdVC2U1CXpREkXUsMynOmGXR/J2ydpyvisQRA0iOYPu04C/pZksvlRwCBJPOcicDlwnpndMF0l++Lb9W7gvLydkr5Cshp6i5kdNWnfe4B/BQM/SksAABjVSURBVFaY2VZJp5FETBNJsvm/MzPfJFRFsVhh0aKsxX/1gqy1a++YH4Aqj0f1b3Plyzuz1q5d41m/JIBDev061nT4Abg8OuV7sZRzXsDtOelwPN+poZJ/T/bmpPbZOj6Qke0s+5axbcN+Sp2+Dt/qeNdY1kJ0+/B+btl79ixx5ev7/fIDheyKkRuHDnJKwqIOfwXJ+t2+BWv3aDaY2NAM37WZIGv+sMvMRkgsXZ9Lg88vB4bNrPYXm30LJjbdlPsFJBry4QclS7CfD9xTJb4TeJaZPQH4KEn4xiAIPOoQTKxemNm4md03U8UD+6Z8ptS/ZnY14KU9Phf4x+rjzey3ZrYj/XgNMKOVkkEwn5jo/Uy3tTrTWbt24ysZAX5/e+r6XgJsNLM/JFEaXd4M/DhvZxqB/wyAzhVNX+cYBLNPGyiWWphuwjk7yH+ESOoDPkAy5MorcwKJ8vnrKdp0PumwrO+IlXPkMQRBjbRYr0bSwcARZvazNBZ0h5n5QZsmMZvBxA4DHgVM9HpWA9dLepqZ3S/paOBLwMlm5s/UTsJMlCrZkeOQk/FhrORfal5e7Ty3i7zMCW7dOQGrOnMCsni52nvkB866bM/jXfmh3ZtrbstJK252y24Z93uUL1mYzZo7mJNP3XsuAEf0+u37m/6sn/KKnMBjK7rWuPLVXf5rM1DIuoX89QJ/Ee62sr+C5OjFG1354Hh2Yn1v2Z+w9+/2I6BFlI+kt5KMQpaSfL9XA18AnlPL8bOmfMzsJuBBc4Sku4BjU2vXQSSesa+vdWl2EMxXcgI8NoO/J0mbcy2Amd0myTc5OjQsdY6ki4H/CxwpaYOkN09R/EPAMhLT3TpJvndfEAStxKiZPbi2I416Wp9IhvuCmZ06zf5Dqv5+C/CWRrUlCOYULTLsAq6S9AGgV9LzgP8X+GGtB0fSwCBoJ2o0s8/SpPT7SBIG3gT8dxKn0n+q9eD2zl4RBPOR1un59AJfMbMvAqQJCHupMdBgWyufSrnAnsGsteFuR+3v2drv1qFOf/bulgE/SFZXIWs12rjbD/hVykl7s6prhyvf61iOhiq+5WSo7FuZugtLXbnnBvHU3jvdsgfkuH8MONe+qujno7+5z8/3fliXb+26diRrwXpctx/QbcR8i+MVO45y5at6stfz6B6/fT/e9gRXfkT/Fle+x3kOR/X5eRW+5UofAa2jfH4OPBeY8DnqJfHt+n9qObitlU8QzDdES1m7eszsQWdHM9uTrueriZjzCYJ2orXmfPZKevLEB0lPYeq47w8jej5B0G60zrDrXcC3JG0i6ZQdALym1oND+QRBu9EiysfMfi/pMcCRqejPZuYvyXcI5RMEbUYr+XYBTwUOIdElT5aEmX21lgPnj/LJmaSzcX/aazgn0FalmPW/Ghr1y+4Y8efetpd8/6Hd5WzAr23jvpXOSzUDcESn79+017KP+tKdT3ZKwq6SHxxt8ZKsBXXTuB/Y69eDfgzxHQv861lUzNZ95d7HuGXXD69w5aOVoivfU8papGbiowewecz3dys55/zj7lUzqnvGtIjykfQ1Ep+udcCEKdSAUD5BMOew+li7JPWQpL/qJtED3zazD6dhUS8hcXe6jsTf0g+PCccCj5suUHweYe0KgnajPgHkR4ETzeyJwDHASZKeDnwSONfMDgd2kIS4yeOPJJPMj4jo+QRBm1GPOZ+0tzKxRqcz3Qw4kSQHO8CFwNnA53OqWQ78SdLvSJTZRN0vrqUNoXyCoN2oXfksnxQh4vw0GB/woDvEdSRJ/v4DuAMYNLOJrAUbgAOnqP/smlviEMonCNqJmeXk2mpmfpJ5wMzKwDGSFgOXAv4Mf/7xV82k/GTaW/mURGFr1moxaFmLSuf2nEvNeZAb+n0rTkdn1r9pZJsfznp9jhXst12HuvKRcrb8jhG/7rwoiXsHbnflnr/WA2N+lNzRin+vbnL8r/La4V0L5PuqrR/OpqY5ss/3AyvkjDvu2rXMlW/vzr4PHQV/1nbjXt9Pr3ehv3xlx1j2+SzoHHVK1gdRf1O7mQ1K+gXwV8BiSR1p72c14IdwBNI5ov8PeCzQRZK3a6+Z1RRcPSacg6DNqId7haQVaY+HNPby84BbgF8Ar0yLvRH4/hTV/DtJfvbbSJxK30IyfKuJUD5B0G7Ux9q1kiTt8Y3A74ErzOwykuSd75Z0O4m5/ctTNsXsdqBoZmUz+0+cXH15tPewKwjmI/Wxdt0IPMmRryeJy1wLQ5K6gHWS/gW4jxl0aKLnEwTtRGt5tb+eRIe8HdgLrAFeXuvBDev5zDBX+5nAaVVtemy6z8t4+hBFozyQnfBcsDjr1b9n1F96r5KfImfZQj8Y24Lu7GTixrJfx6IFfqCtwxZsdeVeupmNnYvdsof2PuDK8wKBXb47GyRr7QN+Cprxsv+btLAje1+357h/3LrVT2IwVvZfueOXZSfK1+32E9feN+RPCj95+b2ufL/ObAqevInv41esd+U37vQtzqNOSqZ7d/nGirrRIu4VwEvN7LPACPARAEnvBD5by8GN7PlcQI252s3sU2Z2jJkdA7wfuGpaxRME8xRVattmgTc6stNrPbiR2SuulnSIs2siV3veLPqpwMUNalYQtD3N9mqXdCrJKuhHSfpB1a6FQM2dhlmdcJ4uV3sagvEkkjFkXh0P5movLvWHJEEwZ5nZIsNG8VuSyeXlwKer5LsBP9yCw6wpn1pytQMvAn4z1ZCrOld79yGrm/8YgmC2afJbb2Z3A3dLei4wbGYVSY8mWSF9U631zKa1qzpX+108lKu92iv2tcSQKwhymVjh3CLWrquBHkkHkmSteD3JXG9NtESu9vTzIuBZwOv2/VyO9SlnAk7jvqVqtORbx4qFrNtAacR3JdhT9N+ALaM5bg2OJWjTHt+ys77HD6iV5+5wx95s+e27ak40AMDGkeww94FhPzDa3j1+QLLtvf457x3JpvzJs2rtHPHr3jnuu6JUnPdhuOxbu/o7fNeIvOBy4+Xse5L37tQLVVqmwy8zG0pToX/OzP5F0rpaD26VXO0ALwMuN7O9jWpTELQ9ta5unh39JEl/RbJM5r9SWc2atyVytaefL2AGXbYgmK8029pVxbtIlsZcamY3SzqUxDesJsK9IgjajRZRPmlIjauqPq8H3lHr8aF8gqDNaHbPR9J5ZvYuST/EUYURyTAI5irN7/l8Lf3/X/elkvZWPgWj0FfKiHu7soGf9nb4Tywv8H6fUwfkBKEq+HV0dmbbBtBb9Ov2rF15gbPy6ugv+NYaz4rT6QRGmwrvnHntc9aQTlm+Ylnbh+c3BVDK8T3bPZ5NkQMw7qS3GcnxMRvOCYI2PF57qp2cS68PdcpesU9NMLsu/f8qSSvSv31nwykIr/YgaCNaZZ2PpLMlbQX+DPxF0gOSPjSTOkL5BEG7YVbb1iAkvRs4HniqmS01syXAccDxkv6h1npC+QRBm9ECPZ/XA6ea2Z0TgtTS9TrgDbVW0t5zPkEw32gNx9LOCc+EaszsAUk1T46F8gmCNqPZE85AXvrk6fY9jLZWPsVihUWLshEH1yzckZHtGfYtIeWcKISPXrLFlXcXshaisudLBqwZ8KMKHr/wNle+vZT1k9rQ50fFe+7Aza58S9n3G/MsVcWi/xaXc6xJ3nUWcn6GLecLMur4QoFvZRp3IjsCjIz5P66LOv3IkSt7dmZkv992sFu2p3ePKx8c8v3GPMZGG/u1agHl80RJuxy5AN/xzqGtlU8QzDuMhk4m19QEs7p4zobyCYI2o9krnOtFKJ8gaDdC+QRBMNs0Il1ys2hr5VMuF9i9JzsRuKGYDXo1sqPmeTAA7t1Te/qTrTv8Sd48V4JbFqxy5Q+MZSecHxjxg3Ud1nOAK39az52ufNni7ETq7nH/now7rg4ApyzNRsi8vX9/t2wpp44nLNnkyvsKWSNJ93LfPeX+fv9+b8wJPrZtNJve55ilG9yya7cd5Mofs6L2vPF5eervcKUzxKyVgontE22tfIJgXjI3dE8onyBoN2LYFQTB7GNADLuCIGgKc0P3hPIJgnYjhl3TIOkrwAuBLWZ2VCo7G3grMBF46ANm9iNJzwM+AXSR+IacaWZXTnsSE6XR7GLLvSPZtCga9a0vlhNkbCQnkJVbR457RR7jOQtEvYBabvAyoEs5lqDyQlc+YlkLzN841iuAoYrvitLnBCo7ecBPUJl3jYuKw678ntFs6pynD9zult3QvcyV/2jT4125dWafz2BOmp1dOWl5lvX4SVX2OAHMDllQc8bgR0Q9rF2S1gBfBfYn6Uudb2afzfuO7vMJHRrZ87kA+HeSC6zmXDObHH5xK/AiM9sk6Sjgp8CBDWxbELQn9fNqLwHvMbPrJQ0A10m6It3nfUfrTiNT51wt6ZAay95Q9fFmoFdSt5n5MUGDYJ6SLDLcd+1jZveR5FvHzHZLuoVZ/sFvRjCxt0u6UdJXJHkr+V4BXJ+neCSdIWmtpLXl3ZFfMJiHVGrcYPnEdyXdzvCqSzsJTwKuTUXTfUfrwmwrn8+T5Gw/hkTrfrp6p6THA58E/nteBWZ2vpkda2bHFgeyK1eDYK4js5o2YOvEdyXdzs/UJS0AvgO8y8x2Mc13tJ7MqvIxs81mVjazCvBF4GkT+yStBi4F3mBmdVmJHgRzjjqmS06jDn4HuMjMvgtTf0frzaya2iWtTMeakORm/2MqX0yS6/l9ZvabmusbE933ZK0NI4uylp2++3w9W8m5A/d31t7bLO7wfXnuG/fP+buiH8hq71jWSjc67jdwSddhrvzoBb7P0hFd92dkx3Vng64B/H50P1d+5a7HZWSP6b3PKQmP6vYzqXxv85Nc+faRvowsL43NxqGs7x7Api2+vLM7axncM5691wBbN/vWwuGcAGajI1n55gHfH68+1Me3S5KALwO3mNlnquTud7QRNNLUfjHwbJJx5wbgw8CzJR1Dopfv4qHh1duBw4EPVaXfeL6Z+eEEg2A+U59gYseTBIK/SdK6VPYB4NSc72jdaaS161RH/OWcsh8DPtaotgTBnKFOSQPN7Nf4+Q0bsqbHI1Y4B0G70eQwqvUilE8QtBtzQ/eE8gmCdkOV5qevqAdtrXxUhk4n04kcH6nuHf7PRbnH98sa3eXfGm+83bXTr2Ok169jyy7fGjI6mrWcWMWvO8/i013wfb46lU3501PIptMBuGHoEFf+l11ZK9hYjrkwL6XO+h1ZHy6A4aGs1TIvdc6O3VnLGEDxft8nbbw3e183DfnWq65Nvnxo2PdVK4xmn8/24QZ+rYyJBYRtT1srnyCYbwiri3tFKxDKJwjajVA+QRA0hVA+QRDMOjHn0xoUStB3f/ZXYNyZz+3fnJ1wBRhb4E9qjg34cnfCebffvnKvX8dQb07e7zGnfE6csi17/UnrjpwVaN5E9GjFn1y9aZef2ufewewk93glLzCa3/DdW3OcgR1XlK34KXJKO33XiL4d/jlLTiC50rh/7T3b/ToK4/51djjp4T3DQT0Ja1cQBE3AYtgVBEETMEL5BEHQJObGqCuUTxC0G7HOJwiC5hDKp/kURyssvCdrbhjvy15W39073TpKi33LkxX8JfyFUvbBd+/yLWmFnEBgQ07QMICiE7Xa8RQB4AH8YGfbcwJZbVy0KCNb2jvklr1j83JXbluyaWVuX+jfp7yhQe9dOdc+lpWNLvYtTAsGfYvUwrv8k44tyJYf7/dv7KK7fPeU0UV++Y7h7PswsqSBAULNoDw3xl1trXyCYF4SPZ8gCJpCKJ8gCGYdA+oQw7kVCOUTBG2Fgc2NOZ+GzYylCce2SPpjlexsSRslrUu3U1L5aVWydZIqaRDrIAiqMZIJ51q2FqclcrWb2UXARQCSngB8z8zWMR0SVshaMgplp1ta8i1S3vFT4Vm78sJa5rg35QYAl2NoMd84lEul7J90rJS1HOWljymP+K9Fx1i27orjNwUgz08N6PANbK61q1L0r6Vr0L/h3YO+pUrl7LUXnSBgAD2bh115YcwPVNYxlH2vrDDDhzZTYs5namaSq30SpwKX1Lc1QTCHmCPKpxVztb8GuDjv4Opc7WNjkas9mG+kjqW1bC1Oq+VqPw4YMrPcLInVudq7uiJXezDPMKBSqW1rcWbV2mVmmyf+lvRF4LJJRV7LFL2eIAhoi15NLbRErvZ0XwF4NfCM2WxTELQX4V4xLTPM1Q7wTOBeM1tf8zlGxui+dVN2R3fW2lC66x63jq69ftS+hcX9XXlhNGvdKAz7KWjK3QtduXJSwngWn1LWnQqA8QW+31N52K9716hj7erxK+/c7Efi63ai/I059QIUHcsYQN9m/4vjWQBV8uvo3ZYTrXFLjimtkvXfK+a1e4c/j9g97p+zOJjN3aSSn9aoLhjYHFnn0xK52tPyvwSe3qj2BMGcIVY4B0HQFGLOJwiCWcesLSxZtdCMdT5BEOwLdVjnI2mNpF9I+pOkmyW9M5UvlXSFpNvS//3AUXWgvXs+xSK2MBs8q7IguxS+Y/QAt4rK8myQLYCRZf4SeTnj7cKov/R+eJmv24dX5KRncTxAyjkr9ceX+ZPcxQFfvrAvG6lsoMeJXgZszJmMrXRkJ6LLS/zzealwAIqjOa+c810ZXeJ/gcq9fvuQn2pneGm2LeWePLcaP5Baqdcv370jO2k/urSRqXMMK/uuQjOkBLzHzK6XNABcJ+kK4HTg52b2CUnvA94HnFWPE04mej5B0E5MhNSoZZuqGrP7zOz69O/dwC3AgcBLgAvTYhcCL23UpbR3zycI5iN1NrWnPphPAq4F9q9ai3c/4K85qQOhfIKgjTDAaje1L5e0turz+WZ2fnUBSQuA7wDvMrNd0kPDSzMzSQ0zrYXyCYJ2wmYUTGyrmR2bt1NSJ4niucjMvpuKN094IkhaCWzZtwbnE3M+QdBmWLlc0zYVSro4XwZuMbPPVO36AfDG9O83At9vyEUAsjZesCTpAeDu9ONyYOssN2G2zzkfrrEZ55yt8x1sZiv2pQJJPyHPJJdlq5mdlFPPXwO/Am7ioURHHyCZ9/kmcBDJd+vVZrZ9X9qcR1srn2okrZ2qizkXzjkfrrEZ52zGNQYx7AqCoEmE8gmCoCnMJeVz/vRF2v6c8+Eam3HOZlzjvGfOzPkEQdBezKWeTxAEbUQonyAImkJbKh9JPZJ+J+kPaTiAj6TyR0m6VtLtkr4hqS7Z26Y43wWS7qzKtFr3LKuSipJukHRZ+rkh1zjF+Rp6jZLuknRTWvfaVNbQsA4553Sz6QaNoy2VDzAKnGhmTyRJw3OSpKcDnyTJiHo4sAN4c4PPB3CmmR2TbtNnWZ057yTxOJ6gUdeYdz5o/DWekNY9sdbmfSRhHY4Afp5+bvQ5IbmvE9f5owacM6iiLZWPJUxE7u5MNwNOBL6dyusWDmCK8zUUSauBFwBfSj+LBl2jd74mMmthHYLm0ZbKBx4cHqwjcXy7ArgDGDSziYTdG0jikzTkfGZ2bbrrnDQD67mS/Khij5zzgH/koeXvy2jgNTrnm6CR12jA5ZKuk3RGKmt0WAfvnDB9Nt2gjrSt8jGzspkdA6wGngY8ZjbPJ+ko4P3peZ8KLKWOEd8kvRDYYmbX1avOR3i+hl1jyl+b2ZOBk4G/l/TM6p2WrAWpdy/TO+eU2XSD+tO2ymcCMxsEfgH8FbBY0kSYkNXAxgae76Q0GpyZ2SjwnyRKsF4cD7xY0l3AJSTDrc/SuGvMnE/S1xt8jZjZxvT/LcClaf2b03AONCKsg3dOM9uc/sBUgC9S5+sMsrSl8pG0QtLi9O9e4Hkkk6S/AF6ZFqtbOICc891a9QURybxEbo75mWJm7zez1WZ2CEka6SvN7DQadI0553tdI69RUn8aPxhJ/cDz0/obFtYh75wT15nysGy6QWNo12BiK4ELJRVJFOg3zewySX8CLpH0MeAGpkhSWKfzXSlpBSBgHfC2Op1vKs6iMdeYx0UNvMb9gUsTvUYH8H/M7CeSfg98U9KbScM6zMI5v6b8bLpBAwj3iiAImkJbDruCIGh/QvkEQdAUQvkEQdAUQvkEQdAUQvkEQdAUQvnMIySZpE9XfX6vpLPTv89O9x9etf9dqezY9POEN/iNki6XdMCsX0QwZwjlM78YBV4uKS/1yk0kCwwneBVw86QyJ5jZ0cBaklQrQfCICOUzvyiRxCv+h5z93yPxKEfSYcBO8vNZXQ0cnrMvCKYllM/84z+A0yQtcvbtAu5NnWZfC3xjinpeSNJTCoJHRCifeYaZ7QK+Crwjp8glJIrnpSROl5P5RRpaZCHw8YY0MpgXtKtvV7BvnAdcT+KlPpnLgE8Ba81sV+oDVc0JZjbb6ZODOUj0fOYhae7tb+KEYDWzIRLn1XNmu13B/CKUz/zl04Br9TKzS8zs+lluTzDPCK/2IAiaQvR8giBoCqF8giBoCqF8giBoCqF8giBoCqF8giBoCqF8giBoCqF8giBoCv8/d4f+kGDv0H4AAAAASUVORK5CYII=\n",
"text/plain": [
"