{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Calculating the RDF atom-to-atom\n",
"\n",
"We calculate the site-specific radial distribution functions of solvent around certain 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.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"
]
},
{
"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",
"import numpy as np\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)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"u = mda.Universe(TPR, XTC)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculating the site-specific radial distribution function"
]
},
{
"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. See [the tutorial on averaged RDFs](average_rdf.ipynb) for more information. The `InterRDF_s` class allows you to compute RDFs on an atom-to-atom basis, rather than simply giving the averaged RDF as in `InterRDF`.\n",
"\n",
"Below, I calculate the RDF between selected alpha-carbons and the water atoms within 15 angstroms of CA60, *in the first frame of the trajectory*. The water group does not update over the trajectory as the water moves towards and away from the alpha-carbon. \n",
"\n",
"The RDF is limited to a spherical shell around each atom by `range`. Note that the range is defined around *each atom*, rather than the center-of-mass of the entire group.\n",
"\n",
"If `density=True`, the final RDF is over the average density of the selected atoms in the trajectory box, making it comparable to the output of `rdf.InterRDF`. If `density=False`, the density is not taken into account. This can make it difficult to compare RDFs between AtomGroups that contain different numbers of atoms."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"ca60 = u.select_atoms('resid 60 and name CA')\n",
"ca61 = u.select_atoms('resid 61 and name CA')\n",
"ca62 = u.select_atoms('resid 62 and name CA')\n",
"water = u.select_atoms('resname SOL and sphzone 15 group sel_a', sel_a=ca60)\n",
"\n",
"ags = [[ca60+ca61, water], [ca62, water]]\n",
"\n",
"ss_rdf = rdf.InterRDF_s(u, ags,\n",
" nbins=75, # default\n",
" range=(0.0, 15.0), # distance\n",
" density=True,\n",
" )\n",
"ss_rdf.run();"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Like `rdf.InterRDF`, the distance bins are available at `ss_rdf.bins`."
]
},
{
"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": [
"ss_rdf.bins"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`ss_rdf.rdf` contains the atom-pairwise RDF for each of your pairs of AtomGroups. It is a list with the same length as your list of pairs `ags`. A result array has the shape `(len(ag1), len(ag2), nbins)` for the AtomGroup pair `(ag1, ag2)`. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"There are 1041 water atoms\n",
"The first result array has shape: (2, 1041, 75)\n",
"The second result array has shape: (1, 1041, 75)\n"
]
}
],
"source": [
"print('There are {} water atoms'.format(len(water)))\n",
"print('The first result array has shape: {}'.format(ss_rdf.rdf[0].shape))\n",
"print('The second result array has shape: {}'.format(ss_rdf.rdf[1].shape))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Index the results array to get the RDF for a particular pair of atoms. `ss_rdf.rdf[i][j][k]` will return the RDF between atoms $j$ and $k$ in the $i$th pair of atom groups. For example, below we get the RDF between the alpha-carbon in residue 61 (i.e. the second atom of the first atom group) and the 571st atom of water."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"array([0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0.41218412, 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0.19906709,\n",
" 0.18622866, 0. , 0.1640152 , 0. , 0. ,\n",
" 0. , 0.13003857, 0. , 0. , 0. ,\n",
" 0. , 0. , 0.09591515, 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0.05427213, 0. , 0. , 0. ,\n",
" 0. , 0. , 0.04435226, 0.04296637, 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ,\n",
" 0. , 0. , 0. , 0. , 0. ])"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ca61_h2o_571 = ss_rdf.rdf[0][1][570]\n",
"ca61_h2o_571"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO29eZxkdXnv//50VVevswA9IMyAg4ALiwuOuKC4KxgF43bFDaO5/O7rSsSoUdRcoyQxaowxUcyN2wWNxihqJIpxR9xhAIUZZZlBkBnAmWGZmV6ma3t+f5zvqT5dXcuppquruup5v1796qpzTp16qpfzOc/yfR6ZGY7jOE7/MtBpAxzHcZzO4kLgOI7T57gQOI7j9DkuBI7jOH2OC4HjOE6f40LgOI7T57gQOC0j6QpJf9ppO5zW8d+dUwsXgh5A0m2SZiRNSrpb0sWSxhP7L5aUl7Q/fG2R9HeS1iSOea2kUjhH/PWxNti6Ii5Ekl4haXP4Odwl6VuSnlx1zGslmaT/UeP1o5I+LmmPpL2Srkzse7qkH4btty3Dx0mFpPeEz3N+1fbzw/b3hOc3JT+zpFOrfw5h235J2fB8naQvhM98n6TPJ479oKQ7JO2TdLukd1a9f0bS30i6M5zzOklra9j//WBHdsl+KH2CC0Hv8AIzGwceDTwGeEfV/g+a2SpgHfAnwBOAn0oaSxzzczMbT3ydtyyWdxmS3gx8BHgfcBhwFPBx4KyqQ88B7gVeU+M0nwAOBh4Rvv95Yt8U8BngL5bU8KXhZhZ+nnPC9pgrgdMSz08Dbqyx7edmVgzPvwrcTfSzPBT4UOLYTwMPN7PVwJOAV0p6UWL/e8P2JwKrgVcDB5IGSnolMJjuIzrVuBD0GGZ2N/BtIkGotf+AmV0NnAkcQiQKi+EYSVeFu7ivSzo43iHpCZJ+Jul+Sb+W9LSw/W+BpwAfiz0OSe+V9NGwf1DSlKS/D89HJB2Iz13vvGHfGkmfDnfvO8MdZCbse62kn0j6ULgb/Z2kM2p9qOAlXQi8wcy+amZTZlYws/8ys79IHPdg4KnAucBzJT0ose/h4ed7rpntNrOSmV2T+B1cZWafA25N84OW9OXg6e2VdKWkExL7LpZ0kaRvhrvlX0o6JrH/2ZJuDK/9GKAmb3c1MBq/R/g+HLbHVAvBU4AP1Nh2ZTjHc4Ajgb8ws73h53ld4udxk5lNJV5bBo4Nrz0IeBPwP83sdovYYmYVIQi/s78C3tbkszl1cCHoMSRtAM4AtjU6zsz2A98l+oddDK8BXgccDhSBfw7vvx74JvA3RHfCbwW+Immdmb0L+DFwXsLj+BHwtHDOxxHdNcYXlCcCN5nZvY3OG469ONhxLJFH9BwgGYJ6PHATMAF8EPi0pFoXxScSXfi+luLzbzazrwC/BV6Z2HcKcDvw3hAaukHSi5ucrxHfAo4jupO+Fvh81f6XE901H0T0e/9bAEkTRHfif0n0ubcDp6Z4v88x5xWcE54nuRI4QdLBkgaATcB/AGsT204Nx0Hkfd4EXCLpHklXS3pq8oSSLpA0CewAxoAvhF0nEf1eXxLE8GZJb6iy533AvxD97TiLwIWgd/hPSfuBO4BdRHdIzbiT6KIa84Rwtx1/PaHBaz8X7symgP8DvCzcgb8KuNzMLjezspl9F9gMPK/OeX4OHCfpECIB+DSwXlGO46lEQkGj80o6LJz/TeEOfhfwj0QXyJjbzeyTZlYCLiESsMNq2HMIsCcR0qjHa5i7WH2B+eGUDcCJwF7gCOA8oovgI5qcsyZm9hkz229ms8B7gEcpkd8Bvha8jCKRSMTe4POArWZ2qZkViMJdaS6W/wacLWmQ6Gf4b1X23A78nugm4lHALWY2A/w0sS0H/DK8ZAORMP8QeBDwD8DXg1DF53w/sAo4mUh49iZeuwZ4KHA08BLgPZKeDSBpE5HofDTF53Lq4ELQO7ww5ACeBjyc6A6wGeuJYtwxvzCztYmvXzR47R2Jx7cTxWcngAcDL00KCvBkogvvAsIFZDPRRf80ogv/z4j+uZNC0Oi8Dw7vf1di378S3UHHVC6AZjYdHo6zkHuAiUYJR0mnEl2Uvhg2fQE4SVJ8AZ4BCsDfmFnezH5EdBF8Tr1zNnivjKT3S9ouaR9wW9iV/P0mL+7TzH2uI0j8nizqMJn8vdXEzH5P5Fm8j+giX+s1cXjoNCIvD+AniW1XBeGC6Odxm5l9OoSFvhjsmOedhLDPdeH49yZeC3Chmc2Y2fVEP/fnBc/j48D5KYTbaYALQY8RLjoXMz8Zt4Bwx/0s5v6JW+XIxOOjiC58e4j+wT9XJShj4Y4PoFa72x8BzyAK6Vwdnj+XKMQShxcanfcOYBaYSOxbbWYnLHyrpvw8nOuFDY45hyjW/itJdzN353tO+H59jdcsts3vK4iS1M8iujPeGLY3i/UD3EXi9xRCYUfWP3wenwXeEr7XIhaCpzD3N/TjxLYrE8dez8LP3+jnkQXiPEf8s0weHz9eTQhLhd9DnMfYIWmxIc++xIWgN/kI8GxJj6reIWlI0mOB/wTuA/7fIt/jVZKOlzRKlFy9NIRd/g14gaTnhrvZYUlPC7kLgD8AD6k614+IQiu/MbM8cAVRfP93ZrY7HFP3vGZ2F/Ad4B8krZY0IOmY6jh0GsxsL/Bu4CJJL1RUBjoo6QxFZY7DwMuIksSPTnz9GfCK4ElcSRQ6eYekbPAgnk6UxCfYN0zkxSh8llwdk1YRCdM9wCjRXXpavkkUy39RsOuNRKGZNPwHkQfzpTr7ryQS7tOIQkIANxB5Sk9nvhB8DThI0jnhd/cSopDPT8PP4v+TdJAiTgHeAHwfwMy2EwnMu8Lf7iOIwlXfYC70Fv8O4vDjY5kTZycFLgQ9SLh4fpboghbztpBDuCfsuwZ4UlW1Rit8jsjzuJsoufrG8N53EN3BvhPYTXS3/hfM/a39E1Hi7z5J/xy2/QwYYe7i8Rui8sDKxSTFeV9DFJf+DZHAXUqdcFQzzOwfgDcTJVnj9zqPSDxfSBSu+KyZ3R1/EZWDZoHTQzz+LKIL017gk8BrzOzG8BanhXNcTuRNzRAJWS0+SxR62xk+W6NwXfXn2AO8FHg/0e/9OOYu2s1eO2Nm3wuhu1r7byb62dxtZveHbWXgKqI79Z8ljr2XqIrqrUQ/jwuAs4J9AH9MlMjeTyT4H2V+zP9sovDfPUTi9n/M7PshlJT8HcQ3DX8INxROSmQ+mMZxHKevcY/AcRynz3EhcBzH6XNcCBzHcfocFwLHcZw+Z8V16ZuYmLCNGzd22gzHcZwVxTXXXLPHzNbV2rfihGDjxo1s3ry502Y4juOsKCTdXm+fh4Ycx3H6HBcCx3GcPseFwHEcp89xIXAcx+lzXAgcx3H6HBcCx3GcPseFwHEcp89xIXBawsz40uY7mC2WOm2K4zhLhAuB0xJb79zH2y69np/csqf5wY7jrAhcCJyWiD2BA4Vyhy1xHGepcCFwWqJQigYZFcsuBI7TK7RVCCSdLukmSdskXdDguBdLMkmb2mmP88ApBiGIBcFxnJVP24RAUga4CDgDOB44W9LxNY5bBZyPD5teERRK5XnfHcdZ+bTTIzgF2GZmt4ZB0l8kGuhdzV8DHyAaVu50ObEAFF0IHKdnaKcQrAfuSDzfEbZVkHQycKSZfbPRiSSdK2mzpM27d+9eekud1BTLHhpynF6jY8liSQPAh4G3NDvWzD5hZpvMbNO6dTXnKjjLhIeGHKf3aKcQ7ASOTDzfELbFrAJOBK6QdBvwBOAyTxh3N3NVQ+4ROE6v0E4huBo4TtLRknLAy4HL4p1mttfMJsxso5ltBH4BnGlmPn6siym6R+A4PUfbhMDMisB5wLeB3wJfMrOtki6UdGa73tdpLx4acpzeo60zi83scuDyqm3vrnPs09ppi7M0VEJDnix2nJ7BVxY7LRGvKPaqIcfpHVwInJYoVFYWe2jIcXoFFwKnJSoLyrzXkOP0DC4ETkt4ryHH6T1cCJyW8Kohx+k9XAiclvCqIcfpPVwInJaYqxpyj8BxegUXAqclvMWE4/QeLgROS3iOwHF6DxcCpyW815Dj9B4uBE5LFMqeLHacXsOFwGmJQtE9AsfpNVwInJbwCWWO03u4EDgt4S0mHKf3cCFwWmKuasg9AsfpFVwInJYoevdRx+k5XAiclvCqIcfpPVwInJbwqiHH6T1cCJyW8F5DjtN7uBA4LVH0XkOO03O4EDgtkfcWE47Tc7gQOC2RnFBm5l6B4/QCLgROSyQXkpU8POQ4PYELgdMS+eKcEPiiMsfpDVwInJYolg0pelzwNhOO0xO4EDgtUSwZo4OZymPHcVY+LgROasyMfKnMSC4LeOWQ4/QKLgROauLk8Egu+rNxIXCc3sCFwElNvIhsdDDyCDw05Di9gQuBk5p4MdlILsoRuEfgOL2BC4GTmtgDGBmMhcA9AsfpBVwInNQUgwcwGjwCn1LmOL2BC4GTGg8NOU5v4kLgpMZDQ47Tm7gQOKmJQ0GV0JALgeP0BC4ETmryxXgdgS8oc5xewoXASU3sEcyFhlwIHKcXcCFwUhPnBOKVxT6lzHF6AxcCJzWFStWQh4Ycp5doqxBIOl3STZK2Sbqgxv7/JekGSb+S9BNJx7fTHueBESeHR71qyHF6irYJgaQMcBFwBnA8cHaNC/0XzOwkM3s08EHgw+2yx3ngxPMHRipVQ+4ROE4v0E6P4BRgm5ndamZ54IvAWckDzGxf4ukY4LeYXUyh6AvKHKcXybbx3OuBOxLPdwCPrz5I0huANwM54Bm1TiTpXOBcgKOOOmrJDXXSMdd91ENDjtNLdDxZbGYXmdkxwNuBv6xzzCfMbJOZbVq3bt3yGuhUKFS1mPBeQ47TGzT1CCStA/4nsDF5vJm9rslLdwJHJp5vCNvq8UXgX5rZ43SO2AMYzblH4Di9RJrQ0NeBHwPfA0otnPtq4DhJRxMJwMuBVyQPkHScmd0Snv4RcAtO1xInh4d9QZnj9BRphGDUzN7e6onNrCjpPODbQAb4jJltlXQhsNnMLgPOk/QsoADcB5zT6vs4y0ch5Ahy2QEG5L2GHKdXSCME35D0PDO7vNWTh9dcXrXt3YnH57d6TqdzxFVDucwAg5kB9wgcp0dIkyw+n0gMDkjaH772NX2V03PEyeFsRQjcI3CcXqCpR2Bmq5bDEKf7iS/82QGRzcirhhynR0i1jkDSmcBp4ekVZvaN9pnkdCtxKGjQQ0OO01M0DQ1Jej9ReOg34et8SX/XbsOc7qNYMgYEmQExOCAPDTlOj5DGI3ge8GgzKwNIugS4DnhHOw1zuo9CuUw2E907ZDMD3mvIcXqEtCuL1yYer2mHIU73UygauSAEgxn3CBynV0jjEfwdcJ2kHwIiyhUsaCnt9D7FcplsRgCeI3CcHiJN1dC/S7oCeFzY9HYzu7utVjldSaFkZAfi0JB8Qpnj9Ah1Q0OSHh6+nwwcTtQ9dAdwRNjm9BnFUpmcewSO03M08gjeTNT6+R9q7DPqtIx2epdCaS5ZPDjgQuA4vUJdITCzc8PDM8zsQHKfpOG2WuV0JYWyVXIE2YzIF10IHKcXSFM19LOU25weJwoNxVVDA5UmdI7jrGzqegSSHkQ0ZWxE0mOIKoYAVgOjy2Cb02UUSpaoGlKlCZ3jOCubRjmC5wKvJRookxwqvx94ZxttcrqUQqk8VzU0MOC9hhynR2iUI7gEuETSi83sK8tok9OlFEuJBWXZAZ9H4Dg9QpoFZSdKOqF6o5ld2AZ7nC6mUCqTy8ZVQyLvVUOO0xOkSRZPAlPhqwScQTS/2OkzoqqhxIIy9wgcpydIs7J43joCSR8iGj/p9BnVC8o8R+A4vUHapnNJRokSyE6fkUwWD2YGfB2B4/QITT0CSTcQrSSGaAj9OsDzA31IMVE+mh3wXkOO0yukSRY/P/G4CPzBzIptssfpYgrlslcNOU4P0jQ0ZGa3A4cAZwEvAk5qt1FOd1IoJhaUhaohMxcDx1nppBlV+W7gEiIxmAAulvSX7TbM6T6KVRPKAEoeHnKcFU+a0NArgUfFjefCDONfAX/TTsOc7qOQXFAWvhfLRjbTSascx3mgpKkauhNIdhsdAna2xxynm4mqhuZ6DQG+qMxxeoBGTec+SlQttBfYKum74fmzgauWxzynm4iqhuJeQ6pscxxnZdMoNLQ5fL8G+Fpi+xVts8bpWswsVA0FjyC0mii6R+A4K55mTeccB4iSwmbMm1AGHhpynF6gUWjoS2b2sqoFZRXM7JFttczpKuLFY8kJZeChIcfpBRqFhs4P35/f4BinT4jnEy+sGnKPwHFWOo1CQ3dJygAXm9nTl9EmpwsphDv/BVVDRfcIHGel07B81MxKQFnSmmWyx+lS4qTwXNWQewSO0yukWVA2CdwQyken4o1m9sa2WeV0HfGg+mSvIZjzFBzHWbmkEYKvhq8k/t/fZ8SD6pO9hmAud+A4zsoljRCsNbN/Sm6QdH69g53eJA4BVfca8qohx1n5pGkxcU6Nba9dYjucLicOAc1NKAsegecIHGfF02gdwdnAK4CjJV2W2LUauLfdhjndRRwCSk4og7mQkeM4K5dGoaGfAXcRtZ5Ozi3eD1zfTqOc7qNSPlq9oMzbUDvOiqfROoLbgdslPQuYMbOypIcCDwduWC4Dne6gWGdBmSeLHWflkyZHcCUwLGk98B3g1cDFaU4u6XRJN0naJumCGvvfLOk3kq6X9H1JD27FeGf5mPMI5vca8vJRx1n5pBECmdk00ZjKj5vZS4ETmr4oWpV8EXAGcDxwtqTjqw67DtgU+hZdCnywFeOd5aNQnl8+OtdryD0Cx1nppBICSU8kmlT2zbAtzUyqU4BtZnarmeWBLxLNPa5gZj8MIgPwC2BDOrOd5SYuE409gWylaqg7PIItO/dS7hJbHGelkUYI3gS8A/iamW2V9BDghyletx64I/F8R9hWj9cD36q1Q9K5kjZL2rx79+4Ub+0sNXEuYDAbCUCui6qGbt09yfM/+hM+8r2bO22K46xImgqBmf3IzM40sw+E57cudXsJSa8CNgF/X8eGT5jZJjPbtG7duqV8aycl1eWj2S7qPnrvVB6Ai67YzpadeztsjeOsPOoKgaSPhO//Jemy6q8U594JHJl4voEas45DVdK7gDPNbLY1853lohIainMElRYTnQ/HTOdLlcdv/fKvyXeBl+I4K4lG6wg+F75/aJHnvho4TtLRRALwcqIFahUkPQb4V+B0M9u1yPdxloFKaKgLy0djITj/mcfx4e/ezEd/cAtvec7DOmyV46wcGq0juCZ8/9FiTmxmRUnnAd8mSi5/JuQYLgQ2m9llRKGgceDLkgB+b2ZnLub9nPZSqJpQlhkQA+qOXkMzhSIAz3/k4dx2zxQfv2I7zzn+QZy0wbunO04aGrWYqDmiMibNqEozuxy4vGrbuxOPn5XOTKfTxGWicdUQRHmCbug1FHsEo7ksf/X8E/jptj289cu/5rI/O5WhbJoCN8fpbxoli58PvAD47/D1yvD1Laou7k7vM1c1NPcnk8sMUOiCCWUzQQhGchnWjA7yvj8+iZv+sJ+v/+rODlvmOCuDZi0mkPRsM3tMYtfbJV0LLFgp7PQu1aMqIQoTdUPV0JxHEN39P/m4CQB27/faA8dJQ9oFZacmnjwp5eucHmKuaigRGhoY6JqqocGMKrblMgNkB8R0vthhyxxnZZBmMM3rgc8k5hbfD7yufSY53UihVGZAUZI4JpdRV1QNzeSLjAzO5QIkMZrLMDVbavAqx3FimgpBqB56VCwEZuYrdvqQQrlcWUQWk80MdEWvoel8idHc/D/l0VzWPQLHSUkajwBwAeh3iiWrzCmOyWbUFb2GpgulSn4gZnQow1TePQLHSYPH+p1UFErleRVDEFcNdd4jmMmXGKkSgrFctlJN5DhOY1wInFQUSlbpMxQTVQ11gUeQLy70CHIZpmY9NOQ4aWi0oOxFjV5oZl9denOcbqVYKlf6DMVEVUPd4RGsGc3N2zY2lPXyUcdJSaMcwQsa7DPAhaCPKJZtXukoRKGhbmgxMZ0vcfia+R7BSC7DlCeLHScVjRaU/clyGuJ0N/lSudJnKCabUVd0+oyqhqpzBBmmG5SPbr7tXtaODnLsoavabZ7jdD2pqoYk/RHReMrheJuZXdguo5zuo1gqz+szBFH5aDdU5swUFiaLm5WPvu0r1/OIw1dz0StObrd5jtP1NE0WS/q/wP8A/gwQ8FLAh8z3GcWSVaaTxeQy6pJ1BEXGhubf04wNZZjOlzCrHbraN1Ng30xhOcxznK4nTdXQk8zsNcB9ZvZe4InAQ9trltNt5EvlhVVDXZAsLpeNA4XyvJXFEHkExbKRr2Pf5GyRSa8qchwgnRDMhO/Tko4ACsDh7TPJ6UaKJVtYNZRRx5PFM4X5Dedi4ue18gTFUpkDhTKTB1wIHAfSCcE3JK0lGiJzLXAb8O/tNMrpPorlcs2qoU7PI6juPBozFlpO1KocinsQ+ToDx4lI02vor8PDr0j6BjDs7Sb6j3zJGMktXFDW6XkEcUJ4pLrX0FAkDLVWF++fjXIDHhpynIhGC8qeYWY/qLWwTJIvKOszoqqh6tDQQMfnETT3CBYKQewRTM4WMTPCmFTH6VsaeQRPBX5A7YVlvqCsz4hyBDVCQx3OEUwnppMlGankCBbe9ceeQNmIEs05H2fp9DeNFpT9VfjuC8scCrUWlA10fh5BHPoZHUzvESRDQvtnCy4ETt/TKDT05kYvNLMPL705TrdSqJEsznZBi4k4R7BgHkHIEdRaVJZMEk/NlsAXFzt9TqPQUPzv8TDgccBl4fkLgKvaaZTTfdQqH81lRKFc7micPS4frdWGGuZCR0km5wmBJ4wdp1Fo6L0Akq4ETjaz/eH5e4BvLot1TtcQhYYWegRmUCrbgrDRclEvWRx7BLUu9Mn1A/t9LYHjpFpHcBiQTzzPh21OH1GoM6EM6OhMgrpCMJiZtz/JlHsEjjOPNE3nPgtcJelr4fkLgUvaZ5LTjUTzCBZWDUHkLQwPdibhOlNZRzD//bOZAXLZgZoLypKhIV9L4DjpFpT9raT/Bp4cNv2JmV3XXrOcbqNQsoWhoeAhdLKEdDpfIjOgiiglGctlai4om5wtIoGZC4HjQMo21GZ2jaQ7CG2oJR1lZr9vq2VOVxFVDS1cUAZ0tAPpdL7E6GCmZrJ6NJetLB5LMjVb5JCxIfZMzi46NGRm3DuV55DxoUW93nG6iTRtqM+UdAvwO+BH4fu32m2Y0z2UyoYZ9UNDHcwR1BpcHzOay9QsH52cLTIxnkNavEfw7a1/4Inv/wH3TeWbH+w4XU6aZPFfA08Abjazo4FnAb9oq1VOVxEvGqs1oQyg0MEpZdOFhdPJYkaHsnUXlK0eHmQ8l120ENx+zxT5Ypk/7D+wqNc7TjeRRggKZnYPMCBpwMx+CGxqs11OFxELQa0JZUBH+w3N5IsLGs7FROMqa3sEY0MZxoayiw4N7Q1DbfZO+3AbZ+WTJkdwv6Rx4Erg85J2AVPtNcvpJuLVw7UWlEHnk8V1PYJclvunZxZsn5otsfGQLGNDmUV7BBUh8ClnTg+QxiM4C5gG/hz4b2A7tRvROT3KXGho4YSy5P5O0EgIonGVtT2CVcNZxocHmWww4L4R+8JCtH2+IM3pAZoKgZlNmVnZzIpmdgnwMeD09pvmdAtxMrjWhDLorEcwky8tGFMZM5rL1M4RHCgylssyPpRh8sDi7ujdI3B6ibpCIGm1pHdI+pik5yjiPOBW4GXLZ6LTaeLy0HpVQx0tHy0UG4aGqnMEpbIxUygxPpxlfKh2eWkaXAicXqJRjuBzwH3Az4E/Bd4JCHihmf1qGWxzuoS6oaHKyuJOl482SBYXSvOa4sUrjceHsowNLb5qaF8QgH0uBE4P0EgIHmJmJwFI+hRwF3CUmXm9XJ8RX+jr9Rrq5NzihsnioSxWNXwmbjg3NhR5BJ4sdpzGOYLKX7iZlYAdLgL9yVzVUL3QUGc8ArMozFM3WRy2J/sNxeWi40NxaCgaV9nq+7oQOL1EI4/gUZL2hccCRsJzAWZmq9tundMV5JstKOtQjuBAoYzZwoZzMXHIaHq2BOPRtv2z80NDxbIxW2ytad5UvkQpJNBdCJxeoNE8Ap/f5wD1k8WdLh+tTCercxFv6BEMZ1k1HP35T84WWxKC5MXfhcDpBdKsI1g0kk6XdJOkbZIuqLH/NEnXSipKekk7bXEWT7HcnaGhuVkEte9nRocWTimLhWAsl52ba9xiniBOEK8ZGfRksdMTtE0IJGWAi4AzgOOBsyUdX3XY74HXAl9olx3OA6dbQ0P1xlTGxB5BclFZPJEsDg0lt6Ul9gKOOnjUPQKnJ2inR3AKsM3MbjWzPPBFolXKFczsNjO7Huhc2YnTlEqyeEGvobhqqNMeQb0cQTyucqFHkAwNteoRJIVgtljmQGFxaxEcp1topxCsB+5IPN8RtrWMpHMlbZa0effu3UtinJOeYh2PIBaGTi0oq+QI6q4jyM47DubaTsdN54CaU8waEQvBkQePAr6WwFn5tDVHsFSY2SfMbJOZbVq3bl2nzek78nWSxYPZziaLZ5p4BJUB9okcweRsicGMGMpmGA/7Ww0N7Ut4BNC+hPHnfnE7F/1wW1vO7ThJ2ikEO4EjE883hG3OCqNe99FOj6psFhqKPYKZqqqh8eAJjA8Nhm2thXb2zhSQYP1BI5Xn7eA/r9vJV67Z0ZZzO06SdgrB1cBxko6WlANeDlzWxvdz2kQ8b6C6xcRgh6uGYo+g7jqCwYU5gmgWQSQEY8EjmJxt7UK+d6bA6uFB1o4MVp63gz2Ts+yenG3LuR0nSduEwMyKwHnAt4HfAl8ys62SLpR0JoCkx0naAbwU+FdJW9tlj7N48nU8gsyAGFAXrCOokyMYGBAjg5kFOYLYI4g9hlZbUe+dKbBmZJDV7RaC/bPsP1D0ZLTTdlINr18sZnY5cHnVtncnHl9NFDJyuphinQllEHkJneo1NF1oHBqC6K5/qmodQSwEAwNiLJdZVNXQmpFB1rRRCGbypYrd907lOWLtyJK/h+PErKH4wegAABbASURBVIhksdNZ4tBPddUQRI3oOhkakmAoW//PeDSXrYSQYH5oCKLmc5OLSBavHsmyOpSf7ptZ+uE0exIhoT0eHnLajAuB05R6VUMQVQ51LjRUYnQwU2kxXYvRqjv+ydki48NzQjA+nGVyEeWja0YGyWYGGB/KtsUjcCFwlhMXAqcp9bqPQtRvqHNVQ/UH18eMDWXntZiYPFBkPPGa8UUMsN87U6yEhdaMDLZJCPJzj/fnGxzpOA8cFwKnKcVyGSlKDlczmFEHF5TVb0EdE42rrCofHZ4vBK2EhswshIYiIVjdNiGY8wK8cshpNy4ETlPypXJNbwAiL6GjoaEUQjAdqoLKZWMqX1qYI2jBIzhQKJMvlRMeQbYtK4v37I8u/kPZAQ8NOW3HhcBpSrFkC6aTxWQz6livoWhMZWMhGMtlmS5EF/q5MZVzr2l1StneROfR+Hu7PILVw1ketGZ4XpjIcdpBW8tHnd6gWCovWEwWMzgw0NFeQ009gqE5jyBeWBavKI4et5YjWD4hyDMxPsRBY7mKd+A47cI9Aqcp+ZLVDw1l1dEWEyODje9lRnPZiicQryAeS3gEY0PZllpMVAvB6uH2eQQT40McMpbjnikXAqe9uBA4TSmWygtWFcdEVUOdm0eQJkdwoFCmVLbKCuLxoWSyOEO+VGa2mE4M9tXwCGYKJfLFpf0Z7JmcZWJVjolVQx4actqOC4HTlGLZai4mg7hqqHMeQTMhSLaiTg6uj4kfp/UKFoSGRqPv+w4srVcQh4Ymxoe4bzrfsfCb0x+4EDhNKTSpGip2qMVEmmRx3Ip6Jl+qtJuurhqC9MNpYiFYPTznESS3LwX5Ypm9MwUmxodYN57DLGoz4TjtwoXAaUqhVK7ZZwiiXkP5DngEZpYqWVyZS5wvVS72qxLrCOLHaWcSVIQgsY4guX0piHMCsUcAvpbAaS9eNeQ0pVhqEBoa6MyCstlimbLV7zwaMzeuslhJGtf0CFK2mdg7U2DVULayuK4dHkG8knhiPMdBY7lom+cJnDbiQuA0pVBuUDWUGehIjqAyi2AwbY6gNG9wfWV/eJx2dXFyVTHMCcFSLiqLF5BNrBrioNEgBF5C6rQRFwKnKYVig6qhjDpSNZSmBTXM5QjiZHF2QPO6la6KhaCFHMGaGkKwlB5BHAaaGBvioLHo/L662GknniNwmlIsl8nWyREMdmgeQTx+Ms3KYog8grgFdbJb6WKSxUkhiJPGe6fb4RHkGB/K9nSbiel8kVd/+pf89q59nTalr3EhcJpSKFllUH01nSofnZtX3GxB2VyOIDmdLCZuQLdYjyCXHWBkMLO0yeLJPKO5DKO5SLQmxoe4p0dzBFvv3MePb9nDD27c1WlT+hoXAqcpUdVQvdBQZxaUNRtcHxPvnw5VQ9VCMDeuMmWO4MB8IYAoPLSU6wjiVcUxE+O5nq0a2rZrct53pzO4EDhNaVY11IkWE80G18ckq4Ki0ND84zNhrnHaZPHemUJlEVnMUvcbioQgV3k+Md67q4u3BwHYvtuFoJO4EDhNKZSbLCjrYo9gKDvAgCLhmJwtMT48uOCY8eFsqvLR2WKJA4VyZURlzJILwf58lUcw1LM5gm1BALbvmsSsfTcUrQ4f6je8ashpSqOVxVFoqBM5gugfe7RJ0zlJjOWixnKTBwqsXzu84JioFXXzFhPV7SViVo8MsvP+mbSmN2XP5CyP3XhQ5fnEqhz3TuUpl42BOiG6lcr23ZMMKFrwd/e+Axy+ZmTJ3+NLm+/gbZdez0Gjgxx76DjHHjrOSevXcvYpRzYcc9pPuEfgNKVYMrJ1LkCDGXWmaqiQLjQUHxOVj5YqOYEkY0MZJlPE+PdVrSqOWTMyuGTrCIqlMvdOL/QISmXjvuneCg8dKJTYcd8Mpxx9MADbd0215X1+sf0e1owMcvqJDwLgW1vu5p1fu4Ff79jblvdbibgQOE1pXDU0gBmUlnk4TdrQEIRW03GyeHihEIynbEVdzyNYytDQvdN5zFiQI4DeW1186+4pzOC5J0QX6G279rflfbbeuY/HPvgg/u5Fj+TL/+tJfPONTwHghh33t+X9ViIuBE5TGlcNqXLMcjKdcmUxhLnFs0Um8wurhiD9lLL6oaHo9UuRK5lrLzHfI4DeW1QW5wee8JBDWDWcZfvupfcIDhRKbNs9yQlHrK5sO2LNMAeP5bhhp3sEMS4ETlOaTSiD5ReCmXyRkcFMqpj5WC7LPZOzmNEWIai0mUhZedSIymKyhBCsW5Wbt69X2L4ryg8cPTHGMevG21JCeuPd+ymVbZ4QSOLE9Wu4YacvYotxIXCa0rjXUHQhXu5FZWlmEcSMDmXYFXr1jNUQgrGU4yr3zUTH1BWCJQgPzXUe7f3Q0Pbdkxx58CjDgxmOPXS8LSWkW8Jd/wlHrJm3/ZHr13DLH/ZzoJB+Ol0v40LgNKXQaEJZplMeQfNZBDGjuQy7gxDU8wj2t+AR1EoWJ/c/ECqhoVVzHsHq4UGyA+o5j2DbrkmOWTcOwDHrxtm1f3bJB/xsvXMfa0YG2XDQ/GqkE9evoVg2b20RcCFwGlIqG2Y06DUUcgQdSBan9ghyWYrBvnpCkC+Wm4rZ3pkCY7nMAu9oSYVgcpZcdqDSDA9gYEAcMt5bQ+xLZeN3e6Y4Zt0YAMceGgnC9iUOD229cy8nHLF6QZnoSRsiD2GL5wkAFwKnCfHFcTBbr3w0+hNa7kVl04USI036DMWM5eYPq1+wP2Xjueo+QzFLKQS7J2dZNz604MLVa4vKdt43w2yxXBGAWBCWMmFcKJW58e798/IDMZ4wno8LgdOQihA0mFCWPG65mMkXGU1RMQQwWmNGcZJ4W7MpZXurZhHELK1HkJ+XH4jptTYTcT4gDg0ddfAogxktacJ4265J8sUyJ65fs2CfJ4zn40LgNCROAjfqNQQs++rilkJDCcGouY5gON2UsnpCsJTjKvfsn+WQRMVQTK95BPEFPxaCbGaAjYeMLWnCeOud0UW+lkcAcNL61Z4wDrgQOA2JVw036jUEy1811FKyeN5EsoWvSRsa2lcnNDQ8mCGXHViSqqHqhnMxE6ty3DOZb2s/nuVk++5JDhmbG8UJUZ5gKXMEW3buZWQww9ET4zX3nxQSxjfe3Z6FbCsJFwKnIfGdfqMJZQD5DiwoS+sRJHMEq4ZqNJ1rITRUSwhgaVYXl8vGPVPz20vErBsfIl8qV0pYVzrbdk1yzKHzL9DHrBvn9nunyReX5m/pN3fu4xGHr6rMl64mDhl5nsCFwGlCnARuNKEsedxyMZ0vNh1KExN7BAOC4cGFn2O84hE0DhHU8whgaWYS7J0pUCpbTSGorCWY6o3w0Pbdc6WjMcceOk6pbNx+zwNPGJfLxtY799bMD8SsXzvCQaODbPGeQy4ETmMqHkGDXkNApTxzuZgppA8NxR5B9ZjKyv4QLpqcrX8hL5TKTOVLbfUIkkPrq6kIQQ+UkN4zOct904VKxVBMLAxLkSe4/d5ppvKluvkBiBLGJ21Y6x4BLgROE+aqhronNBTV/FvqqqFYMFbVqBiKtkcX90atqPfVaS8RsxRCUBlaXydHAL2xujguEY1LRmMesoQlpFvvrL2iuJqT1q/mZk8YuxA4jZmrGmrca2g5k8Vpp5PFxK2na60hiLbPzTWuR70+QzFL4xFEF/l1NUJDh4z1TuO56oqhmLGhLEesGV6SEtItO/cxmBEPPWxVw+PihPFNfZ4wdiFwGjJXNVSnfDQb9xpaPo9guhCG0qRdUBYu9LVKRyESueHBgYaN51IJwfQDFIIQ9qlVPnrwWI4B9YYQbN89yfDgAOvXLhxCc8wS9RzaeudeHnrYKnJ1QpoxcQ7h+j4PD7kQOA0pFBuXj8ZJ5OUMDbUyiyA6LhKAWovJYpp1IK3XZyhm9cgg+2eLlB9ArmTP5CyZAbG2xntkBsTBY7meEIJtuyZ5yMR4zc6xx6wbf8BjK82MrXfua5gfiPGEcURbhUDS6ZJukrRN0gU19g9J+o+w/5eSNrbTHqd14iRwowllsEJCQw08iGYdSOc8gtrnWD2cxax5CWoj9kzOcshYrm5r7YnxIXbv74UcweSCRHHMsYeOV8ZWLpa79x3g3ql8w4qhmLkVxi4EbUFSBrgIOAM4Hjhb0vFVh70euM/MjgX+EfhAu+xxFsdcr6FmVUPd6xHEglEvNATBI2hwEa83pjJmKdpMRO0lFoaFYnphdfFMvsTO+2cW5Adi4u0PJE+wZWfjFcXVnLR+Td8njNs5vP4UYJuZ3Qog6YvAWcBvEsecBbwnPL4U+JgkWRuWT37p6jv45I9vXerT9jzxXXL9XkPR3euHvnMzn/rx75bFplaFIJcdYDCjhqGhsaEsP9t+D8/+8I9q7o/nBTfKEQC86tO/ZKhJXLoed9w3zSlHH1J3/8R4jqt+d29dG1cChVIZMxp6BABv/fKvWT1c+2fdjPtnCkjwiMPTC0GxbJz+kSvrhkC7hTc+8zhe8Kgjlvy87RSC9cAdiec7gMfXO8bMipL2AocAe5IHSToXOBfgqKOOWpQxa0cHOe6w2n98TmNOG67/s1s3PsTrTj2au/fNLKtNTz52guMPb+76x7zreY9g08aD6+5/7ZM28o3r72x4jmPXjTOUrS0+pxx9MC8+eQMzhcWHho47bJwXn7yh7v6Xn3IUhZJhrOw2E4998ME8+diJmvsmxnP876cdw20PcFHZCUesSV1M8JSHruOlj93QtNdUN1DvRuSBonb1LpH0EuB0M/vT8PzVwOPN7LzEMVvCMTvC8+3hmD21zgmwadMm27x5c1tsdhzH6VUkXWNmm2rta6cftBM4MvF8Q9hW8xhJWWANcE8bbXIcx3GqaKcQXA0cJ+loSTng5cBlVcdcBpwTHr8E+EE78gOO4zhOfdqWIwgx//OAbwMZ4DNmtlXShcBmM7sM+DTwOUnbgHuJxMJxHMdZRtqZLMbMLgcur9r27sTjA8BL22mD4ziO05jurpVyHMdx2o4LgeM4Tp/jQuA4jtPnuBA4juP0OW1bUNYuJO0Gbm/xZRNUrVbuQtzGpWMl2Ok2Lg1uY3oebGbrau1YcUKwGCRtrreirltwG5eOlWCn27g0uI1Lg4eGHMdx+hwXAsdxnD6nX4TgE502IAVu49KxEux0G5cGt3EJ6IscgeM4jlOffvEIHMdxnDq4EDiO4/Q5PS8Ekk6XdJOkbZIu6LQ91Ug6UtIPJf1G0lZJ53fapnpIyki6TtI3Om1LLSStlXSppBsl/VbSEzttUzWS/jz8nrdI+ndJw522CUDSZyTtCsOi4m0HS/qupFvC94O60Ma/D7/v6yV9TdLabrMxse8tkkxS7fFsHaSnhUBSBrgIOAM4Hjhb0vGdtWoBReAtZnY88ATgDV1oY8z5wG87bUQD/gn4bzN7OPAousxWSeuBNwKbzOxEovbs3dJ6/WLg9KptFwDfN7PjgO+H553kYhba+F3gRDN7JHAz8I7lNqqKi1loI5KOBJ4D/H65DUpDTwsBcAqwzcxuNbM88EXgrA7bNA8zu8vMrg2P9xNdvNZ31qqFSNoA/BHwqU7bUgtJa4DTiGZcYGZ5M7u/s1bVJAuMhIl8o0DjQcnLhJldSTQTJMlZwCXh8SXAC5fVqCpq2Whm3zGzeNjwL4gmIXaMOj9HgH8E3gbdOXC614VgPXBH4vkOuvAiGyNpI/AY4JedtaQmHyH6Qy532pA6HA3sBv5fCF99StJYp41KYmY7gQ8R3RXeBew1s+901qqGHGZmd4XHdwOHddKYFLwO+FanjahG0lnATjP7dadtqUevC8GKQdI48BXgTWa2r9P2JJH0fGCXmV3TaVsakAVOBv7FzB4DTNH5UMY8Qoz9LCLROgIYk/SqzlqVjjBCtivvZgEkvYsozPr5TtuSRNIo8E7g3c2O7SS9LgQ7gSMTzzeEbV2FpEEiEfi8mX210/bU4FTgTEm3EYXXniHp3zpr0gJ2ADvMLPamLiUShm7iWcDvzGy3mRWArwJP6rBNjfiDpMMBwvddHbanJpJeCzwfeGUXzjw/hkj4fx3+fzYA10p6UEetqqLXheBq4DhJR0vKESXmLuuwTfOQJKK49m/N7MOdtqcWZvYOM9tgZhuJfoY/MLOuupM1s7uBOyQ9LGx6JvCbDppUi98DT5A0Gn7vz6TLEtpVXAacEx6fA3y9g7bURNLpRCHLM81sutP2VGNmN5jZoWa2Mfz/7ABODn+vXUNPC0FIIp0HfJvoH+5LZra1s1Yt4FTg1UR32b8KX8/rtFErlD8DPi/peuDRwPs6bM88grdyKXAtcAPR/19XtB+Q9O/Az4GHSdoh6fXA+4FnS7qFyJt5fxfa+DFgFfDd8L/zf7vQxq7HW0w4juP0OT3tETiO4zjNcSFwHMfpc1wIHMdx+hwXAsdxnD7HhcBxHKfPcSFwuhJJpVAOuEXSf7XaVVLSeyS9NTy+UNKzlsCmEUk/Cs0M20ropPq/23j+50u6sF3nd1YWLgROtzJjZo8OXTrvBd6w2BOZ2bvN7HtLYNPrgK+aWWkJztWMtUBNIQgN6x4o3wReEFogOH2OC4GzEvg5oVmgpHFJ35d0raQbQkMvwr53SbpZ0k+AhyW2XyzpJeHxbXE/eEmbJF0RHj81saDvOkmratjxSsLq2np2SNoYZiF8Mswd+I6kkbDvcaFv/q9CH/0tYfsJkq4K26+XdBzR4q1jEsc+TdKPJV1GWDEt6c3BY9oi6U2J978xfOabJX1e0rMk/VTRXIFToNI76Aqi1gxOv2Nm/uVfXfcFTIbvGeDLwOnheRZYHR5PANsAAY8lWq07CqwO298ajrsYeEl4fBswER5vAq4Ij/8LODU8HgeyVfbkgLsTz+vZsZGo+dmjw74vAa8Kj7cATwyP3w9sCY8/StQnJ36fkXCeLYn3expRI72jw/P4844Fe7cSda6N3/8kohu9a4DPBNvOAv4zcc5XAh/t9O/avzr/5R6B062MSPoVc+2Pvxu2C3hfaCPxPSJP4TDgKcDXzGzaou6trfaU+inwYUlvBNbaXI/7mAkgOd+gnh0QNZb7VXh8DbAx5DhWmdnPw/YvJM71c+Cdkt4OPNjMZurYeJWZ/S48fjLR550ys0miBnZPSbz/DWZWJhKI75uZEQnHxsT5dhF1QXX6HBcCp1uZMbNHAw8muujGOYJXAuuAx4b9fwBaGfdYZO7vvvI6M3s/8KdEd+M/lfTwanuq3qeRHbOJ40pE3kNdzOwLwJnhPS6X9Iw6h041Ok+C5PuXE8/LVbYMh/d0+hwXAqersaij5BuBt4Qk6Rqi2QgFSU8nEgqAK4EXhsqeVcAL6pzyNqKwCsCL442Sjgl30R8g6lo7TwjM7D4go7kZw/XsqPc57gf2S3p82FQZUSnpIcCtZvbPRDmIRwL7iZqp1ePH4fOOKhrA88dhWys8lChc5fQ5LgRO12Nm1wHXA2cTDR7ZJOkG4DXAjeGYa4H/AH5NNKXq6jqney/wT5I2E92tx7wpJF2vBwrUnnT1HaKQDPXsaMLrgU+GkNcYsDdsfxmwJWw/Efismd1D5JlskfT31ScKn/di4CqiiXafCj+nVng6UfWQ0+d491HHSYmkk4E/N7NXL/L14yGej6QLgMPN7PyltLEFWw4DvmBmz+zE+zvdxVLUIztOX2Bm10r6oaSMLW4twR9JegfR/93twGuX1MDWOAp4Swff3+ki3CNwHMfpczxH4DiO0+e4EDiO4/Q5LgSO4zh9jguB4zhOn+NC4DiO0+f8/82Py3iCipyzAAAAAElFTkSuQmCC\n",
"text/plain": [
"