{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Average radial distribution functions\n",
"\n",
"Here we calculate the average radial cumulative distribution functions between two groups of atoms.\n",
"\n",
"**Last executed:** Feb 06, 2020 with MDAnalysis 0.20.2-dev0\n",
"\n",
"**Last updated:** February 2020\n",
"\n",
"**Minimum version of MDAnalysis:** 0.17.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"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import MDAnalysis as mda\n",
"from MDAnalysis.tests.datafiles import TPR, XTC\n",
"from MDAnalysis.analysis import rdf\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. [[3]](#References)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"u = mda.Universe(TPR, XTC)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculating the average radial distribution function for two groups of atoms"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A radial distribution function $g_{ab}(r)$ describes the time-averaged density of particles in $b$ from the reference group $a$ at distance $r$. It is normalised so that it becomes 1 for large separations in a homogenous system.\n",
"\n",
"$$ g_{ab}(r) = (N_{a} N_{b})^{-1} \\sum_{i=1}^{N_a} \\sum_{j=1}^{N_b} \\langle \\delta(|\\mathbf{r}_i - \\mathbf{r}_j| - r) \\rangle$$\n",
"\n",
"The radial cumulative distribution function is \n",
"\n",
"$$G_{ab}(r) = \\int_0^r \\!\\!dr' 4\\pi r'^2 g_{ab}(r')$$.\n",
"\n",
"The average number of $b$ particles within radius $r$ at density $\\rho$ is:\n",
"\n",
"$$N_{ab}(r) = \\rho G_{ab}(r)$$\n",
"\n",
"The average number of particles can be used to compute coordination numbers, such as the number of neighbours in the first solvation shell.\n",
"\n",
"Below, I calculate the average RDF between each atom of residue 60 to each atom of water to look at the distribution of water over the trajectory. The RDF is limited to a spherical shell around each atom in residue 60 by `range`. Note that the range is defined around *each atom*, rather than the center-of-mass of the entire group.\n",
"\n",
"If you are after non-averaged radial distribution functions, have a look at the [site-specific RDF class](site_specific_rdf.ipynb). [The API docs for the InterRDF class are here.](https://docs.mdanalysis.org/stable/documentation_pages/analysis/rdf.html#MDAnalysis.analysis.rdf.InterRDF)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"res60 = u.select_atoms('resid 60')\n",
"water = u.select_atoms('resname SOL')\n",
"\n",
"irdf = rdf.InterRDF(res60, water,\n",
" nbins=75, # default\n",
" range=(0.0, 15.0), # distance in angstroms\n",
" )\n",
"irdf.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The distance bins are available at `irdf.bins` and the radial distribution function is at `irdf.rdf`. "
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.1, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, 1.9, 2.1,\n",
" 2.3, 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, 3.9, 4.1, 4.3,\n",
" 4.5, 4.7, 4.9, 5.1, 5.3, 5.5, 5.7, 5.9, 6.1, 6.3, 6.5,\n",
" 6.7, 6.9, 7.1, 7.3, 7.5, 7.7, 7.9, 8.1, 8.3, 8.5, 8.7,\n",
" 8.9, 9.1, 9.3, 9.5, 9.7, 9.9, 10.1, 10.3, 10.5, 10.7, 10.9,\n",
" 11.1, 11.3, 11.5, 11.7, 11.9, 12.1, 12.3, 12.5, 12.7, 12.9, 13.1,\n",
" 13.3, 13.5, 13.7, 13.9, 14.1, 14.3, 14.5, 14.7, 14.9])"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"irdf.bins"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0, 0.5, 'Radial distribution')"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxU9dX48c/JvidAICwhrAEElC2yue9YFaxa666tlS5atba12sXHatufbZ9a7aPV4lKsVXFXVBR3BEFllX3fEiAhELKQPZnz+2NuwiSZTCaYm0nIeb9eeWXuvd+5c4Zlznx3UVWMMcZ0XWGhDsAYY0xoWSIwxpguzhKBMcZ0cZYIjDGmi7NEYIwxXVxEqANordTUVB04cGCowzDGmE5l+fLlB1S1p79rnS4RDBw4kGXLloU6DGOM6VREZFdz16xpyBhjujhLBMYY08VZIjDGmC7OEoExxnRxlgiMMaaLs0RgjDFdnCUCY4zp4iwRGGNMB6eq/PGd9WzKLXHl/pYIjDGmg5u/Lo8nFu5g7Z4iV+5vicAYYzowj0d56MPNDEqNZ8bYvq68hiUCY4zpwN5fn8vG3BJuPWsoEeHufGRbIjDGmHb2/rpczvv7Z8z+fAcV1bXNlvPWBrYwODWe6WP6uRaPJQJjjGlH5VW1/M/cdewuKOPet9Zz8p8/5vEF2zhcWdOk7Hvr6moDmYSHiWsxWSIwxph29NSi7ewrqmD2907kxZmTOa5PEg+8u5GT//wxL3y1G49HAW9t4OEPtzCkZzwXjXGnb6BOp1uG2hhjOqv9JRU89uk2zhuVxqTBPQCYNLgHX2cX8sd5G7j7tTW8tiKH/3fJ8WzMLWFTXgkPXzHW1doAWCIwxph28/cPNlNZ4+Gu849rcH5M/xRenDmZl5fl8Md5Gzj/4YUkx0YytFcCF57gbm0ArGnIGGPaxcbcYl5cms11UwYyKDW+yXUR4fIT+/PRz0/jguP7cOBwFb84d5jrtQFwuUYgItOAh4Fw4ElVfaDR9b8DZziHcUAvVU1xMyZjjGlv3pnBG0iMieTWs4YGLJuaEM1DV4zj3umjSImLapf4XEsEIhIOPAqcA+QAS0Vkrqquryujqj/zKf9TYJxb8RhjTHspKK3iqx0FFJZVcaismpxDZSzccoDfXTgy6A/39koC4G6NYCKwVVW3A4jIHGAGsL6Z8lcC/+NiPMYY47qdB0q56okv2FtUUX8uKiKMk4b24NrJA0IYWfPcTAT9gGyf4xxgkr+CIjIAGAR83Mz1mcBMgIyMjLaN0hhj2sjW/Ye56okvqPEoz944kcE9E+gWF0lsZDgi7rf1H62OMmroCuAVVfU7xU5VZwGzALKysrQ9AzPGmGBsyi3h6ie/BOCFmyYzvHdiiCMKnpujhvYA/X2O051z/lwBvOBiLMYYE5TSyhruf3s9W/KCX/J5/d5irnziC8IE5szsXEkA3E0ES4FMERkkIlF4P+znNi4kIiOAbsASF2MxxpgWqSp3v7aGpxbt4IfPLqfUz7IPja3KLuTKJ74gJiKMl344haG9Etoh0rblWiJQ1RrgFmA+sAF4SVXXich9IjLdp+gVwBxVtSYfY0xIPfvFLuZ+vZeLxvRl58FSfvfm2oDlv9x+kGue/JLk2Ehe/OEUBvqZH9AZuNpHoKrzgHmNzt3T6PheN2MwxphgrNx9iPvfXs9ZI3rx8HfHMjg1noc/2sJJQ1K5dEJ6k/Kfbc5n5rPL6JcSy3M/mEzv5JgQRN02bGaxMabLKyit4ubnVpCWFMODl48lLEy49axMJg3qzu/eXMu2/MP1ZStranltRQ4/eGYZg1ITePGHUzp1EoCOM2rIGGPanaqSc6icX7++hgOHq3j1x1NJjosEIDxMePiKcZz/8Gfc8vxKfnDyID7amMeCTfmUVtUypn8Kz3zvxHad+OUWSwTGmC6ltLKG11fu4csdBSzdUUBusXfi1/+75HiOT09uULZ3cgx/u3wM35+9jJ+//DVpSdFMH9uPs0b04pRhqURHhIfiLbQ5SwTGmC7ji+0H+eUrX5NdUE5aUjQnDuzOpEHdmTy4B5lp/od8njkijed/MInEmEhG90vq0BPDjpYlAmPMMaOiupY/zdvAl9sLOHVYKtNG92Zc/25U1NTy53c38sySXQzsEcecmZOZNKh70B/qU4emuhx5aFkiMMYcE3YdLOXH/13B+n3FTBjQjdmLd/LEwh30TIwmKjyMPYXlfO+kgdx53ghio46NJp22YonAGNPpvbd2H798eTVhYcLTN2Rx5og0iiuq+WTjfuavyyWvuJIHLx9TvyuYacgSgTGmU3ty4Xb+8M4GxvRP4dGrxpHeLQ6ApJhIZoztx4yx/UIcYcdnicAY02kVlVfz0IdbOGN4T/51bRZRETY16mjYn5oxptN6/svdHK6s4RfnDbck8A3Yn5wxplOqrKnl35/v4JTMVEb1TW75CaZZlgiMMZ3Smyv3sr+kkh+eOiTUoXR6lgiMMZ2Ox6P867NtjOyTxElDbSTQN2WJwBjT6Xy0cT/b8kv54WmDj8mZvu3NRg0ZY1xTUFrF++tyuXRCOpHhLX/vrKn1MHvxTsqqaunfPZb0bnH07xZHWlJ0gw/8fy3YRr+UWC44vo+b4XcZlgiMMc2q9ShhwlF9666sqeUHzyxlxe5CdheUcee0EQHL19R6uP3FVby9el+Ta72TYjh9eE/OGNGL2Mhwlu06xP9cNJKIIJKLaZklAmOMX1U1Hi57fDFDeibw9++ObdVzVZXfvL6WFbsLGZ+RwmMLtjFlSA9Oyezpt3x1rYfb5qxk3ppc7j5/BNdPHciewnKyC8rYXVDGF9sP8s7qfcxZmg1ASlwk3z2xv997mdZzNRGIyDTgYSAceFJVH/BT5nLgXkCBr1X1KjdjMqYrW77rEGPSk4P6Jv3Yp9tYnVPE+r3F/PaC4+iREB306zy1aAevLM/htrMy+dFpQ5j+yCJ+9uLXvHvbKfRMbHifqhoPP31hBfPX5fHbC47jB6cMBmBIzwSG9PTu/3vdlIFU13pYvusQCzbnMyY9hbgo+x7bVlyrV4lIOPAocD4wErhSREY2KpMJ3A2cpKqjgNvdiseYrm51TiGXPraYxxdsa7HsptwSHvlkC1kDulHjUd5ctTfo1/l0037+NG8D54/uzW1nZRIbFc4jV42npKKaO15ahcdzZHvy7fmH+fF/lzN/XR73XDiyPgn4ExkexuTBPfjVtBFMG9076HhMy9xsYJsIbFXV7apaBcwBZjQqcxPwqKoeAlDV/S7GY0yX9u7aXACeXLSDw5U1zZarqfVw5ytfkxQTyazrsjghPZmXl+cEvHdVjYe1e4qY89VufvrCSob3TuJvl48hLMzbtzC8dyL3Th/Fwi0H+MfHW3hleQ6XP76EM/+2gE8353PfjFF8/+RBbfdmTau4WbfqB2T7HOcAkxqVGQYgIp/jbT66V1Xfa3wjEZkJzATIyMhwJVhjjmWqyntrc+nfPZbsgnKeXbKLH5/ufyLWU4t28HVOEY9cNY7u8VF8Z0I6v3tzHev2FjWZwfvRhjwe/GAzm/NKqK71ftPvlxLLE9dNaNJ0c8WJ/Vm09QAPfbgFgEGp8dw5bTiXjU+nV1Ln3vO3swt1I1sEkAmcDqQDn4nI8apa6FtIVWcBswCysrK08U2MMYFt2X+YHQdKuf/i0Xy4Po8nFm7n+qkDmnxYb88/zIMfbObckWn1QzMvGtOX+9/ewMvLchg1/UgiyC+p5GcvriI1IZobTx7MqL5JjO6XzIDucfU1AV8iwgOXHM+wXolMHtydia3YGMa4y82moT2Ab7d+unPOVw4wV1WrVXUHsBlvYjDGtKH5a3MRgfNGpnHrWZkUlFbx3Be7G5Qprazh5y9/TXREGH+4eHT9h3RKXBTnjErjzVV7qKrx1Jf/wzvrqaj2MOu6LO46fwQXjenLoNR4v0mgTmJMJLedncmkwT0sCXQgbiaCpUCmiAwSkSjgCmBuozJv4K0NICKpeJuKtrsYkzFd0nvrchmf0Y1eSTFMGNCNk4em8q/PtlNeVQvAodIqrnryS1bnFPHApSc0aaq5bEI6h8qq+XhjHgALNufz5qq9/Pj0IQztldDu78e0LdcSgarWALcA84ENwEuquk5E7hOR6U6x+cBBEVkPfAL8UlUPuhWTMceqsqoaHvl4C1fO+oL8ksoG17ILyli3t5jzRqXVn7v1rEwOHK7kha92k1tUweX/WsKGfcU8fs0EvuVntu6pmT1JS4rm5WU5lFfV8ts31jA4NZ6fnGELvh0LXO0jUNV5wLxG5+7xeazAHc6PMV3a397fRHlVLT8/d3jQe+rW1Hp4ZXkOD36wmf0llYQJ/P6tdTxy1fj6MvPXeUcLnTfqyJDLiYO6M3lwdx5bsI2nFu2gqLyaZ743kSlD/C/gFh4mfHtcOk8s3M49b64lu6CcF26aTHSE7f17LLD52cZ0AGVVNfxrwXaeXLSD6Y8sYsO+4hafs2FfMdMeXshdr60hvVssr/xoCrefPYy3V+/jow159eXmr8tlRO9EBvSIb/D8284aRn5JJeXVtbxw0+Rmk0CdyyakU+tRXl6ew3cmpLdY3nQelgiM6QC+3F5AVa2Hm88YQmF5NTMe/ZzZn+/AW2luyuNR7nxlNYVlVTx+zXhe/fFUsgZ250enDWFYWgK/fWMthytryC+pZNmuQ34nYE0e3J2HvjuW138ylePTW97YZWivBCYM6Eb3+Ch+/a3jvvF7Nh2HJQJjOoDPtuQTHRHGT8/M5N3bTuHkoanc+9Z6fvbiKr/JYO7Xe1mzp4i7zz+OaaP71I/AiYoI44FLTyC3uIK/vreRD9bnodqwWaiOiHDxuH5NagqB/PPq8bzxk5PoFh919G/WdDihnkdgjAE+25zPpME9iIkMJyYynKeuz+LvH27hHx9tYcKAblw7ZWB92YrqWv7y3kZG9U3i2+P6NbnX+IxuXD9lIM8s2cnAHvEM6BHHiN6JbRJnmk38OiZZjcCYENtTWM62/FJOzUytPyci3H5WJqcN68n972xo0Gfw1KId7C2q4DcXHNfsmP1fnDecPkkx7DhQynmjetuYfROQJQJjQmzh5nwATh3WcInmsDDhb5ePITk2kp++sJKyqhoOHK7ksU+3cfZxvZg6JNXf7QBIiI7gT5ccT3REGNPH9HU1ftP5WdOQMS7Zuv8wA3vEtbjk88ItB0hLiibTz8Ss1IRo/n75WK59+kvue2s9EeFCeXUtd53fcmft6cN7se7359nmLaZF9i/EGBcs3naAsx9cwPdmL6W4orrZcrUeZdHWA5yS2bPZ5puTM1P58WlDmLM0m+e+3M3VkzKCns1rScAEw/6VGOOCxz7dRmJMBEu2HeTSfy4mu6DMb7nVOYUUlVc3aRZq7GfnDGN8RgoJ0RHcdpYtx2XaliUCY9rYmpwiFm45wC1nDOU/N04kr7iCb//zc1bsPtSk7MItBxCBk4c2394P3k1Znr9pMh/ecVqrdgozJhiWCIxpY48t2EpiTARXTcpg6pBUXvvJScRHR3DFrC94z9kcps5nm/M5vl8y3YMYlx8TGW7DN40rLBEY04a25x/m3bW5XDdlAIkxkYB3Ru7rPzmJUX2TuPn5Fby2wrvbV3FFNSuzCzklM3BtwBi3tThqSER64t1ScqBveVX9vnthGdM5zfpsO1HhYXzvpIbbLnaPj+K/N07ipv8s446Xvqa0qpaeCdHUepRTMwP3DxjjtmCGj74JLAQ+BGrdDceYziu3qIJXV+Rw5cQMUv2048dHR/D0DSdyy/Mr+N0baxnQI474qHDGZXQLQbTGHBFMIohT1V+5HokxndxTi7bjUbjplMHNlomJDOexayZwx0tf89bXezn7uF5ERVgLrQmtYBLB2yLyLWdvAWOMH1vySnj+y91cdEIf+nePC1g2MjyMh747lnH9U5g82JZyNqEXTCK4Dfi1iFQBdTNjVFWT3AvLmM5hTU4R//x0K++tyyUuMpybzxga1PPCw4Tvnzyo5YLGtIMWE4Gqts2yhcZ0cqpKXnElm/NK2JxXwoLN+SzccoDEmAhuPn0o3ztpoI3xN51SUGsNOXsMn+ocfqqqbwf5vGnAw0A48KSqPtDo+g3AX4E9zqlHVPXJYO5tTHvweJQl2w/y4tJsPt20n+KKmvpraUnR3DltONdMHkCSM1TUmM4omOGjDwAnAs85p24TkZNU9e4WnhcOPAqcA+QAS0Vkrqqub1T0RVW9pfWhm45q6c4C1u4p4lBZNYVlVRSVVzN9TF/OOi6t5Sd3EEVl1Tz7xU5eXJZNdkE5STERnD+6D6P6JZHZK5HMtAS/I4OM6YyCqRF8Cxirqh4AEXkGWAkETATARGCrqm53njcHmAE0TgTmGFJT6+HqJ7+kqsYDQHJsJNW1HtbsKeLMEb06zbr4t7+4kk825TN1SA9+ce5wzhvVm5hI26jdHJuCXYY6BShwHre8ualXPyDb5zgHmOSn3KUiciqwGfiZqmY3LiAiM4GZABkZGUG+vAmFnEPlVNV4uH/GKK6aNIDwMGHOV7u567U1rM4pYkz/lCbPeWlZNhv3lXDPRSNDEHFTa/cU8cmmfH553vCgO3+N6cyCGcD8/4CVIjLbqQ0sB/7YRq//FjBQVU8APgCe8VdIVWepapaqZvXsabMwO7IdB0sBGNEniXBn96zzj+9DVERY/dIKviqqa3ng3Y08/fkOlu4saHI9FB77dBuJ0RFcO2VAqEMxpl20mAhU9QVgMvAa8CowRVVfDOLee4D+PsfpHOkUrrv3QVWtdA6fBCYEE7TpuHYe8CaCgT4boifHRnLOcWm8tXof1bWeBuXnrtpLQWkV0RFhPPzhlnaN1Z9t+YeZt3Yf1021DmDTdTSbCERkhPN7PNAHb9NODtDXOdeSpUCmiAwSkSjgCmBuo9fo43M4HdjQuvBNR7PrYBkJ0RGkJjRcTfOS8f0oKK1iwab8+nOqytOf72BE70R+ce5wFm094EqtoKK6lh//dznvrtnXYtnHPt1GdEQY3z/JxvibriNQjeAO5/ff/Pz8b0s3VtUa4BZgPt4P+JdUdZ2I3OcMRwW4VUTWicjXwK3ADUf1Lkyb+XTTfm76z7L6zt7W2nGglIGpcU06hU8d1pMe8VG8vvJIpXDJ9oNszC3h+ycN4prJA0hNiGrzWoGqcvdra3h3bS5Pf74jYNmcQ2W8sXIPV07MsPkApktptrNYVWc6D89X1QrfayIS1KLozrIU8xqdu8fn8d20PPrItJNNuSXc/NwKSqtq2bK/hFF9gx0XcMTOg6WM7tf0eZHhYVw0pi/Pf7WbovJqkmMj+ffnO+keH8X0sX2JiQznR6cN4Q/vbGDpzgJOHNi9Ld4Ssxfv5PWVexjYI47luw5x4HBls8M+n/hsOyKB1woy5lgUTGfx4iDPmU7s4OFKbnxmaf03+c15Ja2+R3Wth5xD5Qzy6R/w9e1x/aiq8TBvzT52HSzlww15XD0po35Y5tWT/NcKKmtqWbuniL2F5Q36GFSV/SUVfLWjgLdX7+XA4coGz/ti+0H+8M4GzhmZxiNXjcej8PGG/X5jyy+pZM7SbC4Zl07flNhWv3djOrNmawQi0hvvENBYERkH1NX1k4DAq2qZTqWqxsOPn1tBfkklz980iStmfcHG3NYnguyCMmo9ysBU/4nghPRkBveM5/UVe9iSd5iIMOGayUdG5sRGNawV9IiP4oWvdvPK8hwOlXmXuRKB1IRokmMj2VtYTlnVkZXRI8OFbx3fh2snD6BvSiw3P7eCAT3iePDyMSRER9AvJZb31+dx+Yn9m8T29Oc7qK718KPTh7T6fRvT2QWaR3Ae3jb7dOBBn/MlwK9djMm0I1XlnjfX8tWOAh6+YiwTBnRnSM8ENh9FIth10LtB+6BU/98TRIRLxvXjf9/fzJo9RVxwfJ8mWy9ePWkAjy/Yxo2zl1JcUUNEmHDOyDSmje5NaWUtucUV7C+uoLCsmlMyUxnYI54BPeJIio1k7qq9vLo8hzdX7SU2MpzwMGHWtVn1O4WdMzKNF77aTVlVDXFRR/7pH66s4b9f7OL80X0Y1EwSM+ZYFqiP4BngGRG5VFVfbceYTDt6bcUe5izN5pYzhjJjbD8AhvdOZOmO1o/e2eEMHR3QTNMQwMVOIiivrvW7+mZsVDh3ThvBUwt38MPT+vKdrHR6JQa3T+/4jG7cOW04b67ay9xVe5l56mCG9kqov37uyDRmL97Jwi0HOG9U7/rzLy7NpqSihpmnWt+A6ZqCmVk8WkRGNT6pqve5EI9pZ2+s2sPgnvHccc6w+nPDeyfy5qq9FFdUt2os/c6DpSRGR9AjwEbs6d3iOH14T6pqPJyQ3nSWMcDlWf25PKtp800w4qIiuHJiBldObDoD/cRB3UmKieCD9Xn1iaCm1sPTi3YwcVB3v7OejekKgkkEh30exwAXYuP9jwlVNR6W7izgihMzCAs7MtxzeJp35fEteSVMGBD86B3v0NH4FtcTmnVt1tEF/A1Fhodx5ohefLQhj5paDxHhYcxbm8uewnLund7ku44xXUYwM4v/5vPzR+B0wOrQx4CVuw9RUe1hypCGu2QN7+1NBK3tMN55sLTZjmJfURFhIdue8dxRvTlUVs3yXYdQVZ74bDuDU+M5a0SvkMRjTEdwNP8b4/B2IJtO7vNtBwkTmmyX2C8lloToiFZ1GFfVeNhzqJxBPTr2gLJTh/UkKjyM99fn8eWOAtbsKeIHpwxuUCMypqsJZj+CNYA6h+FAT8D6B44Bi7ce4Ph+ySTHNuwHEBGGpSW0qkaQfagMjwbuKO4IEqIjmDq0Bx+sz2PngVJ6xEdxyfh+oQ7LmJAKpo/gQp/HNUCes3yE6cRKK2tYlV3ITc2MlBneO5H31uaiqkHtIVC/2FwnGH557sje/Pr1NewuKOP2szNtnwHT5QXTR7AL6IF3U5lLgOPdDsq476sdBdR4lJOGpPq9PjwtkUNl1eSXVPq93ljd0NHOMA7/7OO8/QHREWFcO9mWmjamxUQgIvfg3SegB5AKzBaR37odmGkb76/LZf663CbnP996gKjwMCYM6Ob3ecOcDuNNQS41sfNgKUkxEXSL6/hLN/dKiuHisX354WlDbHE5YwiuaehqYEzdwnPOHsargD+4GZhpGw+8u5Hc4go+/eXpDSZmLd52kPEDUoiN8t8sUjeEdFNuCadktrwZ0K6DZUENHe0oHrpiXKhDMKbDCGbU0F688wfqRNNogxnTMZVV1bDjYCllVbUNFnIrKK1i/b7iZpuFAHokRJOaEM2mIDuMdxwobbAZjTGm8wi0Mc3/icg/gCJgnbNV5b+BtUBhewVoAntx6W625x/2e21jbgmqMLhnPHOWZrPNKbdk20EApg5tPhEADO+dENQqpJU1tewtLO8UHcXGmKYC1QiW4d2f+HW8i8x9AnwK/AZ40/XITIvyiiv41atreGLhdr/XN+wrBuDBy8cSExHGX97bCMDn2w4QHxXOCemB9xsYnpbE5rzDeDxaf27r/sNc9H+LGuwkll3gHTra3GJzxpiOraVF50wHVvfNflV2kd/rG/YVkxgTwZj0ZH502hD+9sFmlu0sYPHWA0wa3IPI8MAtg8N7J1BeXUv2oTIG9IivX6l0zZ4i7nhpFfNvP5W4qAh2HPCuOmpNQ8Z0ToGahl5yfq8RkdWNf4K5uYhME5FNIrJVRO4KUO5SEVERCc0iNJ3U4m0HANiUW0xZVdOpHev3FnNcnyREhBtPGUSvxGjuem0NOw+WMbXRshL+DO+dBBxZauKt1ftYvO0gl01IJ7ugnL+8twmAXQebblhvjOk8Ao0aus35fWGAMs0SkXDgUeAcvJveLxWRuaq6vlG5ROe1vjya1+nKFm87SHJsJEXl1azJKWKSz1IRHo+yMbekfhXPuKgI7jhnGHe9tgaAk1roHwDIdJZw3pxbwtQhPfjD2+s5vl8yf770BBKiI5i9eCfTRvdmx4FSkmMj6RZg1VFjTMfVbI1AVfc5H+azVXVX458g7j0R2Kqq21W1CpiDd1JaY/cDfwYq/FwzzcguKCPnUDnXTx0IwKrshv33uwrKKKuqZWSfpPpzl01IJ7NXAqkJUfXDQwOJj44go3scG/NKeOjDLeQfruT+i0cTHibcOW04Gd3juPOV1WzYV2wdxcZ0YgEbiVW1FvCISOt3Mfduc5ntc5zjnKsnIuOB/qr6zlHcv0ur6x+48IQ+9O8e2yQR1HUUH+eTCCLCw3jq+hP59w0Tg15kbVhaIku2HWT24p1ccWIGY501++OiIvjLZSewu6CMFbsLO/xic8aY5gW7H8EaEfkAKK07qaq3fpMXFpEwvFtg3hBE2ZnATICMjKYbjnRFi7cdIDUhisxeCYzt341lOxvuKLZhXzHhYUJmWkKD8xmt/MAe0TuRDzfk0S0ukjvPG97g2uTBPbh+ygCeWbKrwy82Z4xpXjCJ4DXnx5f6K9jIHsB3m6l0Gk5ESwRGA586s1F7A3NFZLqqLmvwYqqzgFkAWVlZwbz2MU1VWbztIFOGpCIijO2fwltf7yWvuKJ+D+D1e4sZ0jP+Gy+oNrKvt0Zx1/kj/PYB3DltBCWVNUwb3bvJNWNM5xBMIkhR1Yd9T4jIbc0V9rEUyBSRQXgTwBXAVXUXVbUI79pFdff8FPhF4yRgmtp+oJT9JZVMcTqH65prVu4urP9A3rCvmBMHBb+7WHPOHZnGf2+cxElD/Y8yio+O4MHLx37j1zHGhE4wS0xc7+fcDS09yVmq+hZgPt6tLV9S1XUicp+ITG9VlKaBxXUzg50hoKP6JhEZLvX9BIVlVewtqmjQUXy0IsLDODkztdOsIWSMab1mawQiciXeb/CDRGSuz6UkoMD/sxpS1XnAvEbn7mmm7OnB3NPAkm0H6JscwwCnvT8mMpzj+iSxKvsQAOv9dBQbY0xzAjUNLQb24W2++ZvP+RIgqAllpu15PMqSbQc5Y0SvBt/Sx/ZP4dXlOdR6lPV7LREYY4IXaB7BLlX9FDgbWKiqC/AmhnTA2glCZFNeCYfKqpnaaOXQsf1TKK2qZcv+EjbsK6FnYjQ9E22tfWNMy4LpI8GrHPMAABc1SURBVPgMiBGRfsD7wLXAbDeDMs2r6x+Y0miJiLoO46+zC9mwr9hqA8aYoAWTCERVy/BuU/lPVf0OMMrdsExzlmw7wMAecfRLiW1wflBqPMmxkSzdeYgt+0vapKPYGNM1BJUIRGQK3p3K6mYA227fIVBT6+HL7QVNagMAIsKY/im8u2Yf1bXKcX1aXkLCGGMguERwO3A38Loz/HMw3r0JTDvLLa6gpLKGE9JT/F6v6ycArEZgjAlaixPKnE7iBT7H24FvtLyEOTrF5d6lppvbIH6c008QHRHGIFsEzhgTpEDzCB5S1dtF5C38LCmhqjYprJ2VVFQDkBjjPxHU7Tg2vHciES1sOmOMMXUC1QiedX7/b3sEYlpWUuGtESTG+P9r65EQzZj+KfVLTxhjTDACbVW53Pm9oLkypn0VOzWCpGZqBABv/GRqe4VjjDlGBGoaWkOAVUZV9QRXIjLNaqlGANiaQMaYVgvUNFS3ReXNzu+6pqJrCG4ZatPGissD9xEYY8zRCNQ0tAtARM5R1XE+l34lIiuAZjejN+4oqawhJjKMqAjrCDbGtJ1gJ5Sd5HMwNcjnmTZWUlFttQFjTJsLZmOaG4GnffYtLgS+715IpjnF5TUkBegfMMaYoxHMhLLlwJi6RODsLGZCoNhqBMYYFwT99dISQOiVVNQEHDFkjDFHw9r6O5HiimqSYq1GYIxpW64mAhGZJiKbRGSriDQZZSQiPxKRNSKySkQWichIN+Pp7EoqrI/AGNP2Ak0ouyTQE1X1tUDXRSQceBQ4B8gBlorIXFVd71PseVV93Ck/HXgQmBZk7F1Ocbn1ERhj2l6gr5cXBbimQMBEAEwEtjqrlSIic4AZQH0iUNVin/Lx2ES1ZlXVeKis8ViNwBjT5gJNKPveN7x3PyDb5zgHmNS4kIjcDNwBRAFn+ruRiMwEZgJkZGR8w7A6p5ZWHjXGmKMV1NdLEbkA7/aUMXXnVPW+tghAVR8FHhWRq4DfAtf7KTMLmAWQlZXVJWsNxc46Q0mxViMwxrStFjuLReRx4LvATwEBvgMMCOLee4D+PsfpzrnmzAEuDuK+XVJ9jSDaagTGmLYVzKihqap6HXBIVX8PTAGGBfG8pUCmiAwSkSjgCmCubwERyfQ5vADYElzYXU/d7mQ2j8AY09aC+VQpd36XiUhf4CDQp6UnqWqNiNwCzMe72f3Tzp7H9wHLVHUucIuInA1UA4fw0yxkvOpqBDaPwBjT1oJJBG+LSArwV2AF3pE9TwZzc1WdB8xrdO4en8e3BR9q1xbMXgTGGHM0gllr6H7n4asi8jYQY8tNHB2PR6lVJfIo9hMutlFDxhiXBJpQdqaqfuxvYpmItDihzDT15KLtvPBVNp/84vRWP7e4ogYRSIy2GoExpm0F+lQ5DfgY/xPLgplQZhrZkneYHQdKOVxZQ0IrP9BLKqpJiIogLMy2ojTGtK1AE8r+x/n9TSeWGUehs9XkvsJyMtMSW/Xc4vIa6yg2xrgiUNPQHYGeqKoPtn04x7YiJxHsLapodSLw7k5mzULGmLYX6JOl7pNqOHAiR+YAXAR85WZQx6qisiM1gtYqtkRgjHFJoKah3wOIyGfAeFUtcY7vBd5pl+iOMYXlVYC3RtBaJRU19E6KabmgMca0UjDjGNOAKp/jKuecaaWi8qOvEdjuZMYYtwTzyfIf4CsRed05vhh4xr2Qjk0V1bVUVHsAyC1ufY3AdiczxrglmAllfxSR94CTnVPfU9WV7oZ17KmrDQDsbWWNQFWtRmCMcU1QnyyqulxEsnGWoRaRDFXd7Wpkx5hCp6O4d1IM+4oqUFVEgpsTUF5dS61HbVaxMcYVwSxDPV1EtgA7gAXO73fdDuxYU1cjOK5PImVVtfWriQajrmySJQJjjAuC6Sy+H5gMbFbVQcDZwBeuRnUMKizz9reP6JMEwN6i4JuHjuxOZk1Dxpi2F0wiqFbVg0CYiISp6idAlstxHXPqZhWP6O2dnrGvFYmg2BKBMcZFwXyyFIpIAvAZ8JyI7AdK3Q3r2FPsJIKRdTWCwuBHDh3ZptKahowxbS+YGsEMoAz4GfAesA3/C9GZAArLqgkPEwalxhMRJq2qEdTtRZBkNQJjjAuCGT5a9+3fAzwjImHAlcBzbgZ2rCksryIpJoKI8DDSkmLY15oagVObsM5iY4wbmq0RiEiSiNwtIo+IyLnidQuwHbi8/UI8NhSV15ASFwVAn+SYVnYW1+1OZonAGNP2AjUNPYt3wbk1wA+AT4DvABer6oxgbi4i00Rkk4hsFZG7/Fy/Q0TWi8hqEflIRAYcxXvoFArLqkh22vj7pMSyrxXrDRVXVBMRJsREtn5nM2OMaUmgpqHBqno8gIg8CewDMlQ1qE8wEQkHHgXOAXKApSIyV1XX+xRbCWSpapmI/Bj4C/Ddo3gfHV5ReTXdnBpB3+QY5q8LflJZibO8RLAT0IwxpjUCfcWsXxNBVWuBnGCTgGMisFVVt6tqFTAHb8dzPVX9RFXLnMMvgPRW3L9TKSqvJiXOqREkx1BV4+FgaVULz/Ky5SWMMW4KlAjGiEix81MCnFD3WESKg7h3PyDb5zjHOdecG2lmxrKIzBSRZSKyLD8/P4iX7ngKy6pJ8WkaAoLuMC4ut70IjDHuaTYRqGq4qiY5P4mqGuHzOKktgxCRa/BOUvtrM7HMUtUsVc3q2bNnW750u6j1KMUV1STXNw15E0GwHcYlFTU2YsgY4xo3ex/3AP19jtOdcw2IyNnAb4DpqlrpYjwhU1JRjSo+ncXeDWaC3ZfAmoaMMW5yMxEsBTJFZJCIRAFXcGS7SwBEZBzwL7xJYL+LsYRU3YJzdU1DPeKjiIoIC3rkUHFFtdUIjDGucS0RqGoNcAswH9gAvKSq60TkPhGZ7hT7K5AAvCwiq0RkbjO369TqlqCu6ywWEWcuQXCJwFsjsERgjHGHq+0NqjoPmNfo3D0+j8928/U7iroF55J91grqkxwTVNNQrUc5XGlNQ8YY99gMpXZQ3zQUdyQR9E0OblLZYVtwzhjjMksE7aDI2YsgOTaq/lyflBhyiyuo9WjA59oS1MYYt1kiaAd1fQQNm4ZiqfUo+SWBB0rVJQJbedQY4xZLBO2gqLyauKhwoiKO/HH3dYaQtjSX4MgS1NY0ZIxxhyWCdlBYfmRWcZ0+ycHNLq5bgtpGDRlj3GKJoB0Ulh2ZVVynbnZxSxvU1NcIYq1pyBjjDksE7aCovIrkRh/kSbERxEWFt7hl5ZGN661GYIxxhyWCdlBUXk1KbMMaQd2kstziwDWC4vpNaaxGYIxxhyWCdlBYVt1gDkGdPsmxQdUIYiPDiQy3vypjjDvs06UdFJZXk+w3EcQE1UdgtQFjjJssEbisorqWqhpPgzkEdfqkxLK/pJLKmtpmn1/s7E5mjDFusUTgsvoF5xr1EQCMy0hBFT5Yn9fs861GYIxxmyUClxWWe5eX8NdHcGpmTzK6x/Gfxbuafb53dzKrERhj3GOJwGVFfpaXqBMeJlwzOYOvdhawYZ//3T+9u5NZjcAY4x5LBC7ztwS1r8uz+hMdEcZ/lvivFRTbXgTGGJdZInBZUVnTJah9pcRFcfHYfryxck/9ctW+vLuTWY3AGOMeSwQuO7IXQdPO4jrXThlAeXUtryzPaXC+ssY74shGDRlj3ORqIhCRaSKySUS2ishdfq6fKiIrRKRGRC5zM5ZQKSyvIjxMiI8Kb7bM6H7JTBjQjWeX7MTjsz9Bic0qNsa0A9cSgYiEA48C5wMjgStFZGSjYruBG4Dn3Yoj1ArLvCuPikjActdNGcDOg2V8tiW//lzdyqO2BLUxxk1uftWcCGxV1e0AIjIHmAGsryugqjudax4X4wip5mYVN3b+6D7cn7CBpz/fSU2tsmjrgfqkEMzzjTHmaLmZCPoB2T7HOcCko7mRiMwEZgJkZGR888jaUXF5dbMjhnxFRYRx1cT+/OPjrXy2OZ+YyDAmDurBVRMzmDqkRztEaozpqjpF47OqzgJmAWRlZQXe5LeDKSyrJjWh+Y5iXzedOpjEmEhG9UtifEY3YiKb71cwxpi24mYi2AP09zlOd851KYXlVQztlRBU2cSYSG46dbDLERljTENujhpaCmSKyCARiQKuAOa6+HodUlFZcE1DxhgTKq4lAlWtAW4B5gMbgJdUdZ2I3Cci0wFE5EQRyQG+A/xLRNa5FU8o1HqU4ooaSwTGmA7N1T4CVZ0HzGt07h6fx0vxNhkdk4rLA88qNsaYjsBmFruoyBKBMaYTsETgopYWnDPGmI7AEoGLCsu8exEk+9mUxhhjOgpLBC6ypiFjTGdgicBFRdY0ZIzpBCwRuKgwwO5kxhjTUVgicFFhWTUJ0RFEhtsfszGm47JPKBcVBbngnDHGhJIlAhet2VPIgB5xoQ7DGGMCskTgkm35h9mcd5hzR6aFOhRjjAnIEoFL5q/LBeC80b1DHIkxxgRmicAl763NZWz/FPokx4Y6FGOMCcgSgQv2FJazOqeIaVYbMMZ0ApYIXDB/rbdZaNooSwTGmI7PEoEL3luXy4jeiQxMjQ91KMYY0yJLBG0sv6SSpTsLrFnIGNNpWCJoYx+sz0MVSwTGmE7DEkEbe29dLoNS4xmelhjqUIwxJiiuJgIRmSYim0Rkq4jc5ed6tIi86Fz/UkQGuhmP24rKq1m89QDnjeqNiIQ6HGOMCYpriUBEwoFHgfOBkcCVIjKyUbEbgUOqOhT4O/Bnt+JpDx9tyKPGo9YsZIzpVNzcvH4isFVVtwOIyBxgBrDep8wM4F7n8SvAIyIiqqptHcxLS7N5YuH2tr5tA/mHK+mTHMMJ/ZJdfR1jjGlLbiaCfkC2z3EOMKm5MqpaIyJFQA/ggG8hEZkJzATIyMg4qmBS4iLJTEs4qucGKzMtgfNH9yEszJqFjDGdh5uJoM2o6ixgFkBWVtZR1RbOHdWbc22ClzHGNOFmZ/EeoL/Pcbpzzm8ZEYkAkoGDLsZkjDGmETcTwVIgU0QGiUgUcAUwt1GZucD1zuPLgI/d6B8wxhjTPNeahpw2/1uA+UA48LSqrhOR+4BlqjoXeAp4VkS2AgV4k4Uxxph25GofgarOA+Y1OnePz+MK4DtuxmCMMSYwm1lsjDFdnCUCY4zp4iwRGGNMF2eJwBhjujjpbKM1RSQf2NXKp6XSaLZyB2Qxtp3OEKfF2DYsxuANUNWe/i50ukRwNERkmapmhTqOQCzGttMZ4rQY24bF2DasacgYY7o4SwTGGNPFdZVEMCvUAQTBYmw7nSFOi7FtWIxtoEv0ERhjjGleV6kRGGOMaYYlAmOM6eKO+UQgItNEZJOIbBWRu0IdT2Mi0l9EPhGR9SKyTkRuC3VMzRGRcBFZKSJvhzoWf0QkRUReEZGNIrJBRKaEOqbGRORnzt/zWhF5QURiQh0TgIg8LSL7RWStz7nuIvKBiGxxfnfrgDH+1fn7Xi0ir4tISkeL0efaz0VERSQ1FLEFckwnAhEJBx4FzgdGAleKyMjQRtVEDfBzVR0JTAZu7oAx1rkN2BDqIAJ4GHhPVUcAY+hgsYpIP+BWIEtVR+Ndnr2jLL0+G5jW6NxdwEeqmgl85ByH0myaxvgBMFpVTwA2A3e3d1CNzKZpjIhIf+BcYHd7BxSMYzoRABOBraq6XVWrgDnAjBDH1ICq7lPVFc7jErwfXv1CG1VTIpIOXAA8GepY/BGRZOBUvHtcoKpVqloY2qj8igBinR354oC9IY4HAFX9DO+eIL5mAM84j58BLm7XoBrxF6Oqvq+qNc7hF3h3QgyZZv4cAf4O3Al0yNE5x3oi6Adk+xzn0AE/ZOuIyEBgHPBlaCPx6yG8/5A9oQ6kGYOAfODfTvPVkyISH+qgfKnqHuB/8X4r3AcUqer7oY0qoDRV3ec8zgXSQhlMEL4PvBvqIBoTkRnAHlX9OtSxNOdYTwSdhogkAK8Ct6tqcajj8SUiFwL7VXV5qGMJIAIYDzymquOAUkLflNGA08Y+A2/S6gvEi8g1oY0qOM4Wsh3y2yyAiPwGbzPrc6GOxZeIxAG/Bu5pqWwoHeuJYA/Q3+c43TnXoYhIJN4k8JyqvhbqePw4CZguIjvxNq+dKSL/DW1ITeQAOapaV5t6BW9i6EjOBnaoar6qVgOvAVNDHFMgeSLSB8D5vT/E8fglIjcAFwJXd8A9z4fgTfxfO/9/0oEVItI7pFE1cqwngqVApogMEpEovB1zc0McUwMiInjbtTeo6oOhjscfVb1bVdNVdSDeP8OPVbVDfZNV1VwgW0SGO6fOAtaHMCR/dgOTRSTO+Xs/iw7Wod3IXOB65/H1wJshjMUvEZmGt8lyuqqWhTqexlR1jar2UtWBzv+fHGC88++1wzimE4HTiXQLMB/vf7iXVHVdaKNq4iTgWrzfslc5P98KdVCd1E+B50RkNTAW+FOI42nAqa28AqwA1uD9/9chlh8QkReAJcBwEckRkRuBB4BzRGQL3trMAx0wxkeAROAD5//O4x0wxg7Plpgwxpgu7piuERhjjGmZJQJjjOniLBEYY0wXZ4nAGGO6OEsExhjTxVkiMB2SiNQ6wwHXishbrV1VUkTuFZFfOI/vE5Gz2yCmWBFZ4Cxm6CpnJdWfuHj/C0XkPrfubzoXSwSmoypX1bHOKp0FwM1HeyNVvUdVP2yDmL4PvKaqtW1wr5akAH4TgbNg3Tf1DnCRswSC6eIsEZjOYAnOYoEikiAiH4nIChFZ4yzohXPtNyKyWUQWAcN9zs8Wkcucxzvr1oMXkSwR+dR5fJrPhL6VIpLoJ46rcWbXNheHiAx09kJ4wtl34H0RiXWuneism7/KWUd/rXN+lIh85ZxfLSKZeCdvDfEpe7qILBSRuTgzpkXkDqfGtFZEbvd5/Y3Oe94sIs+JyNki8rl49xWYCPVrB32Kd2kG09Wpqv3YT4f7AQ47v8OBl4FpznEEkOQ8TgW2AgJMwDtbNw5Ics7/wik3G7jMebwTSHUeZwGfOo/fAk5yHicAEY3iiQJyfY6bi2Mg3sXPxjrXXgKucR6vBaY4jx8A1jqP/w/vOjl1rxPr3Getz+udjnchvUHOcd37jXfiXYd35dq61z8e7xe95cDTTmwzgDd87nk18H+h/ru2n9D/WI3AdFSxIrKKI8sff+CcF+BPzjISH+KtKaQBpwCvq2qZeldvbe2aUp8DD4rIrUCKHlnjvk4q4Lu/QXNxgHdhuVXO4+XAQKePI1FVlzjnn/e51xLg1yLyK2CAqpY3E+NXqrrDeXwy3vdbqqqH8S5gd4rP669RVQ/eBPGRqirexDHQ53778a6Caro4SwSmoypX1bHAALwfunV9BFcDPYEJzvU8oDXbPdZw5N99/fNU9QHgB3i/jX8uIiMax9PodQLFUelTrhZv7aFZqvo8MN15jXkicmYzRUsD3ceH7+t7fI49jWKJcV7TdHGWCEyHpt4VJW8Ffu50kibj3RuhWkTOwJsoAD4DLnZG9iQCFzVzy514m1UALq07KSJDnG/Rf8a7am2DRKCqh4BwObLHcHNxNPc+CoESEZnknKrfolJEBgPbVfUfePsgTgBK8C6m1pyFzvuNE+8GPN92zrXGMLzNVaaLs0RgOjxVXQmsBq7Eu/FIloisAa4DNjplVgAvAl/j3aVqaTO3+z3wsIgsw/ttvc7tTqfraqAa/ztdvY+3SYbm4mjBjcATTpNXPFDknL8cWOucHw38R1UP4q2ZrBWRvza+kfN+ZwNf4d3R7knnz6k1zsA7esh0cbb6qDFBEpHxwM9U9dqjfH6C056PiNwF9FHV29oyxlbEkgY8r6pnheL1TcfSFuORjekSVHWFiHwiIuF6dHMJLhCRu/H+v9sF3NCmAbZOBvDzEL6+6UCsRmCMMV2c9REYY0wXZ4nAGGO6OEsExhjTxVkiMMaYLs4SgTHGdHH/HwMxglp6cDIpAAAAAElFTkSuQmCC\n",
"text/plain": [
"