{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# RDF Fitting (Optional)\n", "\n", "This tutorial runs through some of the functionality that can be used to determine solvation radii from the RDF. It covers the default functionality and the use of custom functions to determine the radii.\n", "\n", "It is recommended that you understand the basics of `solvation_analysis` before exploring this tutorial." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# imports\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import MDAnalysis as mda\n", "\n", "from solvation_analysis.rdf_parser import plot_scipy_find_peaks_troughs, identify_cutoff_scipy\n", "from solvation_analysis.tests import datafiles\n", "from solvation_analysis import Solute\n", "\n", "from scipy.signal import find_peaks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets make a convenience function to allow us to quickly load a precomputed RDF from a `.npz` file" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def npz_loader(npz):\n", " arr = list(np.load(npz).values())[0]\n", " return arr" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now lets use the functionality in `solvation_analysis.rdf_parser` to quickly plot the data and determine the peaks and troughs of an the RDF.\n", "\n", "We are going to be considering the RDF of the Butyro Nitride molecules in our test system. See the basics tutorial and the README in tests/data/ for a detailed description." ] }, { "cell_type": "markdown", "source": [ "## The Default Minima Finding Method" ], "metadata": { "collapsed": false } }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "solvation radius is 2.7\n" ] }, { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMtUlEQVR4nO3dd3hUVfrA8e+bQkIgCYReE0CKSAlIUVGKKCAqgl3B7mJbXXdtuFhQ17Ksrq67lh923WAXe0fBRVEgdCkiSgIkQCgJCenJ+/vj3oEQUiZlZlLez/PMM3PbuW8myTtnzj33HFFVjDHGNDxBgQ7AGGOMb1iCN8aYBsoSvDHGNFCW4I0xpoGyBG+MMQ2UJXhjjGmgLMGbSolIloh0r+19A0FE/ioizwc6DgAR6eq+X8G1WOYIEdnkljtZRD4TkcuqWdYsEflvbcVWjfO/LCJ/C9T5GwJL8HWAiGwRkRz3n3KH+4fdvMT2l0UkX0Qy3cdaEXlYRKJL7HO5iBS5ZXge/6mN+FS1uar+Vtv7BoKqPqSqVwc6DgBVTXbfr6JaLPZ+4D9uue+r6mmq+kotlm/qEUvwdceZqtociAcGAXeW2j5bVSOBNsAVwHHA9yLSrMQ+i91/bM/jj/4I3NQpscDPgQ7C1A2W4OsYVd0BfIGT6MvanquqS4FJQCucZF8lIrJARP4mIj+4Nf2PRKSViCSIyH4RWSoicSX2VxE5yn39sog8JSKfuN8mfhKRHhXs+7TbTJAlIt+LSHsReUJE9onIBhEZVOLYGSKy2S13nYhMqeBnOOzru4iMFpFtJZbvEJHtblkbRWSsu/5gs4OIxLnxXiYiySKyW0RmliijqYi84sa6XkRuL3mOUvGIiDwuIrtEJENEVotIvxLlPCYiSe62Re46z/lDSvxeHhaRJe5+H4hIjLvtExG5sdQ5V4vI5BLLm4HuwEfu+x3mlnm1u/1y99yPuj/T7yJyWonju4nIQvc9+wpoXd77X+IYz88wXURSRCRVRG4psT2oxO91j4i85fmZ3O1vi/OtNUNEvhORY8o5T6SIfCsiT7rv9UT3byTT/T3fWlmsjZEl+DpGRDoDpwG/VrSfqmYCXwEnVfNUFwKXAJ2AHsBi4CUgBlgP3FvBsRcB9wEt3TgfrGDf84G7cJJFnnue5e7yO8A/S+y7GefniXbL/6+IdKjiz4WI9Ab+CAx1v/WMB7ZUcMiJQG9gLHCPiBztrr8XiMNJmqcC0yooYxwwEugFtAAuAPa42x4FjgVOwHl/bweKyynnUuBKoCNQCDzprn+l5PlFZCDO7+5TzzpV7QEk434bVNW8MsofDmzEef9nAy+IiLjb5gKJ7rYHgKq03Y8BeuK8DzNE5BR3/U3AZGCU+zPtA54qcdxn7nFtcf4uEkoXLCKtgPnA96p6kzrjq7wAXOP+fvsB31Qh1kbDEnzd8b6IZAJbgV1UnGA9UnAShsdxIpJe4nFcBce+pKqbVTUD559ss6p+raqFwNs4zUTleU9Vl7j7JlDOtw3XPFVNVNVcYB6Qq6qvuu3Ob5Y8j6q+raopqlqsqm8Cm4BhFZRdniIgDOgrIqGqukVVN1ew/32qmqOqq4BVwEB3/fnAQ6q6T1W3cSjZlqUAiAT6AKKq61U1VUSCcBL2n1R1u6oWqeoP5SRfgNdUda2qHgDuBs4X5yLsB0BPEenp7ncJ8Kaq5lf6bhwuSVWfc9//V4AOQDsR6QoMBe5W1TxV/Q74qArl3qeqB1R1DU5F4SJ3/TXATFXd5v7Ms4BzPd9aVPVFVc0ssW2glLi2hPOhsBB4W1XvKrG+AOf3G+X+fpZX7W1oHCzB1x2T3drIaJwkUenXY5wa3N4Syz+qaosSjx8rOHZnidc5ZSw3p3w7SrzOrmRfr88jIpeKyErPBxROzcyb9+EwqvorcDNOwtglIm+ISMcKDinv5+mI84HrUfJ16XN+A/wHp3a6U0TmiEiUG384zrcTb5Q8RxIQCrR2E+BbwDT3Q+Mi4DUvyyzp4M+qqtnuy+a4tWv3g6Xk+b1VOm7P+x0LzCvxO12P8wHcTkSCReQRt/lmP4e+ZZX8nZ8ONAWeLXW+c4CJQJLbrHR8FWJtNCzB1zGquhB4GedrfbnE6WVzCvA/P4TlcyISCzyH07TSSlVbAGsBKeeQA0BEieX2JTeq6lxVPREnwSjw92qElQp0LrHcpaKdVfVJVT0WOAanqeY2YDeQi9MM5o2S5+iKU1Pd7S6/AkzFaUrKVtXFXpbpjVSgpRx+0b5rFY4vHXeK+3orcFqpike4qm4HLgbOwvk7jsZpDoPDf+fPAZ8Dn5aMTVWXqupZOE077+N8+JlSLMHXTU8Ap4pIfOkN7oWzY3H+qPfhfB1uCJrhJOI0ABG5AqcGX56VwEQRiRGR9jg1dtxje4vIySIShpNcc3BqjVX1FnCniLQUkU44Hz5lEpGhIjJcREJxPnxygSJVLQZeBP4pIh3dWuvxbmxlmSYifUUkAqfL4zuebpRuQi8GHqN6tfdyqWoSsAy4T0SaiMiJwJlVKOJuEYlwL5JegdP8Bk7N+0H3AxwRaSMiZ7nbInGuy+zB+bB+qJyy/4hz3eBjcS5ONxGRqSISraoFwH6q9/tt8CzB10Gqmga8itMG63G720a/192WCJxQ6it1vaWq63AS12KcZpz+wPcVHPIaTnv5FuBLDiUUcNrfH8Gp+e7AqeX9tRph3Q9sA34Hvsa5KFxe23kUTm1zH04TxR4OfQu7FVgDLMX5/f2d8v/3XsP5BrcDp2nnplLbX8V5bzw9gZ4VkdLNF9V1Mc5F2L0414BercKxC3EuuM8HHlXVL931/wI+BL50/35/dM+BW34SsB1Y5247gntRdTrOt4EPcN6XS4AtbtPOtVR8AbzREpvwwxjviMh1wIWqOspH5S8A/quq5d5pKyKXAtPd5qeAE6c77e9AqHvR3dQhVoM3phwi0kGcW/+D3K6Xt+D0BApUPBHA9cCcQMVg6hdL8MaUrwnwf0AmTj/rD4CnAxGIiIzHuT6xE6e/uj/PPVUOHwLD87A7Zus4nzbRiMgWnH+OIqBQVYf47GTGGGMOE+KHc4xR1d2V72aMMaY2+SPBe61169YaFxcX6DCMMabeSExM3K2qbcra5usErzjdoxT4P1U94uKQiEzH6QJF165dWbZsmY9DMsaYhkNEyr3j2NcXWUeo6mCcwbNuEJGRpXdQ1TmqOkRVh7RpU+aHkDHGmGrwaYJX1RT3eRdO97LqDBxljDGmGnyW4EWkmYhEel7jDCO61lfnM8YYczhftsG3wxlFznOeuar6uQ/PZ4zxs4KCArZt20Zubm6gQ2nwwsPD6dy5M6GhoV4f47MEr868nAMr3dEYU29t27aNyMhI4uLiODRviKltqsqePXvYtm0b3bp18/o4u5PVHxISIC4OgoKc54QjJq0xpl7Kzc2lVatWltx9TERo1apVlb8p1al+8A1SQgJMnw7Z7twKSUnOMsDUqYGLy5haYsndP6rzPlsN3tdmzjyU3D2ys531xhjjQ5bgfS05uWrrjTG1ZtasWTz6qDMs/4YNG4iPj2fQoEFs3uztDIr1myV4X+tazqxn5a03xlSbqlJcXFzmtvfff5+zzjqLFStW0KOHtzMo1m+W4H3twQchIuLwdRERznpjTI1t2bKFo48+muuvv57BgwfzwAMP0Lt3b0455RQ2btwIwKeffsoTTzzB888/z5gxYwIcsf/YRVZfcy+kbpk2k64kExTb1UnudoHVNECjRx+57vzz4frrnUtPEyceuf3yy53H7t1w7rmHb1uwwLvzbty4kZdeeomrrrqKyy+/nBUrVlBYWMjgwYM59thjmThxItdeey3Nmzfn1ltvrdoPVY9ZDd4fpk6lG1sIphi2bLHkbkwti42N5bjjjuN///sfU6ZMISIigqioKCZNmhTo0ALKavB+1K9foCMwxrcqqnFHRFS8vXVr72vspTVr1uzga+u2eYjV4P0kNBTOPDPQURjTsI0cOZJ58+aRk5NDZmYmH330UaBDCiirwfuBKhQUwOLFgY7EmIZt8ODBXHDBBcTHxxMbG8tJJ50U6JACyqdzslbVkCFDtCFO+FFUBCHuR2kderuNqbH169dz9NFHBzqMRqOs91tEEsub79qaaPwgOBhmzICwsEBHYoxpTCzB+0l4OOTlQTn3YBhjTK2zBO8HO3fCrFnO67y8gIZijGlELMH7QcmkbvMiGGP8xRK8H+TnO89Tp0KJ7rrGGONTluD9wJPgJ02CJk0CG4sxpvGwBO8HngT/ySeQnh7QUIwxjYgleD9o2RI6dYJXX4VGMgy1MQ3Oli1bmDt37sHll19+mT/+8Y+1UnbJstLS0hg+fDiDBg3if//7X43KtQTvB7Gx8OKLzmu7yGpM/VQ6wddEUVFRudvmz59Pnz59WLFiRY3vxLWhCvyg5J2sluBNQ3Xz5zezcsfKWi0zvn08T0x4osJ9Dhw4wPnnn8+2bdsoKiri7rvv5o477uDiiy/m22+/paCggDlz5nDnnXfy66+/ctttt3Httdeiqtx+++189tlniAh33XUXF1xwQbnrZ8yYwfr164mPj+eyyy6jZcuWpKSkMGHCBDZv3syUKVOYPXt2uXE2b96cv/zlL3zxxRc89thjbNq0iYcffpgOHTrQq1cvwsLCWLlyJbfffjs5OTnEx8ezePFimjZtWu33zxK8H3z1FZx2mvPaErwxtevzzz+nY8eOfPLJJwBkZGRwxx130KVLFxYvXsyf//xnLr/8cr7//ntyc3M55phjuPbaa3nvvfdYuXIlq1atYvfu3QwdOpSRI0fyww8/lLn+kUce4dFHH+Xjjz8GnGaVlStXsmLFCsLCwujduzc33ngjXbp0KTPOAwcO0K9fP+6//35SU1O5+OKLSUxMJDo6mjFjxjBo0CDi4+O5//77WbZsGf/5z39q/N5YgvcDz0VWgJycwMVhjC9VVtP2lf79+3Prrbdyxx13cMYZZxxs1vCMBd+/f3+ysrKIjIwkMjKS8PBw0tPTWbRoERdddBHBwcG0a9eOUaNGsXTp0nLXR0VFHXHusWPHEh0dDUDfvn1JSkoqN8EHBwdzzjnnAPDTTz8xevRo2rRpA8AFF1zAL7/8UuvvjSV4P/Ak+GeeKXvGG2NM9fXq1YvExEQ+/fRT7rzzTsaNGwdAmDv4U1BQ0MHXnuXCwkLKG2ixKgMwliw3ODiYwsLCcvcNDw8nODj44LI/xq23i6x+4EnwJ5/sTGpgjKk9KSkpREREMG3aNG699VaWL1/u1XEjR47kzTffpKioiLS0NL777juGDRtW7vrIyEgyMzNrJebhw4ezYMEC9uzZQ0FBAW+//XatlFua1eD9wJPg33wTzj4bjjkmsPEY05CsWbOG2267jaCgIEJDQ3nmmWc4t/TkrmWYMmUKixcvZuDAgYgIs2fPpn379uWub9WqFSEhIQwcOJDLL7+cli1bVjvmDh06MGvWLI4//ng6dOjA4MGDK+xZU102HrwfLF0KL73kNNHMng233RboiIypHTYevH/ZePB10NCh8OSTzmvrRWOM8RdrovGD7GynmSY42HrRGNPQDR8+nLxS44K/9tpr9O/f3++xWIL3gyefhDvvdG52shq8MQ3bTz/9FOgQDrImGj/wXGSNirIEb4zxH5/X4EUkGFgGbFfVM3x9vrrI0zzz3XfOwGPGGOMP/mii+ROwHjjyNrBGIj/fGQfeukcaY/zJp000ItIZOB143pfnqes8Cf6tt+DDDwMdjTENR3p6Ok8//bTPz3P55Zfzzjvv+Pw8tc3XbfBPALcDxT4+T502cSLcdRc89pjTF96YxiphTQJxT8QRdF8QcU/EkbAmoUbllZfgfXHTUH3kswQvImcAu1Q1sZL9povIMhFZlpaW5qtwAmrcOLj1VggPt26SpvFKWJPA9I+mk5SRhKIkZSQx/aPpNUryM2bMYPPmzcTHxzN06FDGjBnDxRdfTP/+/cnNzeWKK66gf//+DBo0iG+//RY4cqKOM844gwULFgDwwgsv0KtXL0aPHs0f/vCHw/b77rvvOOGEE+jevfvB2nxqaiojR44kPj6efv361XiCjtrmyzb4EcAkEZkIhANRIvJfVZ1WcidVnQPMAedOVh/GEzC7djljwoeHQ0ZGoKMxJjBmzp9JdkH2YeuyC7KZOX8mU/tPrVaZjzzyCGvXrmXlypUsWLCA008/nbVr19KtWzcee+wxwBnKYMOGDYwbN67CERtTUlJ44IEHWL58OZGRkZx88skMHDjw4PbU1FQWLVrEhg0bmDRpEueeey5z585l/PjxzJw5k6KiIrKzs8stPxB8VoNX1TtVtbOqxgEXAt+UTu6NxU03wZgx0LSpdZM0jVdyRnKV1lfHsGHD6NatGwCLFi3ikksuAaBPnz7ExsZWmOCXLFnCqFGjiImJITQ0lPPOO++w7ZMnTyYoKIi+ffuyc+dOAIYOHcpLL73ErFmzWLNmDZGRkbX2s9QG6wfvB56LrNZEYxqzrtFdq7S+Opo1a3bwdXnjbIWEhFBcfOiyYK5b66psXK6SQwN79h05ciTfffcdnTp14pJLLuHVV1+tduy+4JcEr6oLGmsfeDiU4P/1L3CbAY1pdB4c+yARoRGHrYsIjeDBsQ9Wu8yKhvAdOXIkCQlO+/4vv/xCcnIyvXv3Ji4ujpUrV1JcXMzWrVtZsmQJ4NT+Fy5cyL59+ygsLOTdd9+t9PxJSUm0bduWP/zhD1x11VVeD1XsLzZUgR94Eny7doGOxJjA8bSzz5w/k+SMZLpGd+XBsQ9Wu/0doFWrVowYMYJ+/frRtGlT2pX4J7v++uu59tpr6d+/PyEhIbz88suEhYUxYsQIunXrRv/+/enXrx+DBw8GoFOnTvz1r39l+PDhdOzYkb59+x6crak8CxYs4B//+AehoaE0b968ztXgbbhgP/DM4nTvvbBkCdxxR0DDMabWNLThgrOysmjevDmFhYVMmTKFK6+8kilTpgQ6rINsuOA66MYbnceXX8Lddwc6GmNMeWbNmnWwy2O3bt2YPHlyoEOqEWui8QN3nl1+/hkKCpwukyWmZjTG1BGPPvpooEOoVZbg/WDjRmjWzOlFA05XyRIX+40xxicswfvBWWdBfDyccIKzbAneGOMP1gbvByX7wYPd7GSM8Q9L8H7gSfCXXAJ790KHDoGOyBjTGFiC9wNPgm/a1JnwI8jedWPqnPo6JHBFLNX4gSfB//KL0wc+ufaG3jCmfklIgLg4p5YTF+csG5+xBO8HTz0FF13kJPbZs2Hr1kBHZEwAJCTA9OmQlASqzvP06TVO8lu2bKFPnz5cdtllDBgwgHPPPZfs7GwSExMZNWoUxx57LOPHjyc1NRWA5557jqFDhzJw4EDOOeecMkeAvPvuu7n88sspLi5mxowZ9O3blwEDBnDrrbfWKFZ/swTvB1OnwvDhhy6y2oBjplGaORNKJ9PsbGd9DW3cuJHp06ezevVqoqKieOqpp7jxxht55513SExM5Morr2Sme56zzz6bpUuXsmrVKo4++mheeOGFw8q6/fbb2bVrFy+99BLp6enMmzePn3/+mdWrV3PXXXfVOFZ/sm6SPqYKixc730atF41p1Mprm6yFNssuXbowYsQIAKZNm8ZDDz3E2rVrOfXUUwFnhqcObu+GtWvXctddd5Genk5WVhbjx48/WM4DDzzA8OHDmTNnDgBRUVGEh4dz9dVXc/rpp3PGGfVrzESrwftYXh6MGAGvvupcZAVL8KaR6lrOsMDlra8CETlsOTIykmOOOYaVK1eycuVK1qxZw5dffgk4F1P/85//sGbNGu69996DwwWDM757YmIie/fuBZyhhZcsWcI555zD+++/z4QJE2ocqz9Zgvex/HznuWQ/eGuiMY3Sgw9CxOHDBRMR4ayvoeTkZBYvXgzA66+/znHHHUdaWtrBdQUFBfz8888AZGZm0qFDBwoKCg4OJ+wxYcIEZsyYwemnn05mZiZZWVlkZGQwceJEnnjiCVauXFnjWP3Jmmh8rGSC79bNWQ6xd900RlPdYYFnznSaZbp2dZL71OoPF+xx9NFH88orr3DNNdfQs2dPbrzxRsaPH89NN91ERkYGhYWF3HzzzRxzzDEHm2FiY2Pp37//EePJn3feeWRmZjJp0iTmzp3LWWedRW5uLqrK448/XuNY/anS4YJF5FHgJVX92dfBNMThglNSoFMn+L//czoMGNOQ1IXhgrds2cIZZ5zB2rVrAxqHP/hiuOANwBwR+UlErhWRikfAN4fJy3OemzRxau833ACffRbYmIwxjUOlCV5Vn1fVEcClQBywWkTmisgYXwfXELRpA+++60z6ERQETz8NS5cGOipjGo64uLhGUXuvDq8usopIMNDHfewGVgF/EZE3fBhbg9C8OZx9ttNNMiTEeVgvGtOQ1KVZ4Rqy6rzPlSZ4EfknTjPNROAhVT1WVf+uqmcCg6p8xkYmPR2++AJ273aWmza1BG8ajvDwcPbs2WNJ3sdUlT179hDu6YrnJW/6c6wF7lLVI+/nhWFVOlsjtGEDTJjgtLtPmOB0lbRukqah6Ny5M9u2bSMtLS3QoTR44eHhdO7cuUrHeJPgp6rqiyVXiMh8VR2rqhlVOlsjVLKbJEBUVOBiMaa2hYaG0q1bt0CHYcpRboIXkXAgAmgtIi0Bz61iUUBHP8TWIJRO8L/+GrhYjDGNS0U1+GuAm3GS+fIS6/cDT/kwpgaldII3xhh/KTfBq+q/gH+JyI2q+m8/xtSglE7wDzwAwcHw178GLiZjTONQURPNyar6DbBdRM4uvV1V3/NpZA3ECSc4vWh69HCW5893ni3BG2N8raImmlHAN8CZZWxTwBK8F9q2hXHjDi2Hh8O+fYGLxxjTeFTURHOv+3yF/8JpeDZvhpUrYeJEpw98eLj1gzfG+Ic3Nzr9SUSixPG8iCwXkXGVHWccX34J554LGW6HUkvwxhh/8WaogitVdT8wDmgLXAE84tOoGpCCAufZc5G1TRvrC2+M8Q9vbnTy9H+fiDNs8CopPX2KKVfpXjT/tv5Ixhg/8aYGnygiX+Ik+C9EJBIo9m1YDYf1gzfGBIo3Cf4qYAYw1B2PpglOM02FRCRcRJaIyCoR+VlE7qthrPWSJ8GHhjrPc+fCWWcFLh5jTONRaRONqhaLyE6gr4hUZbK5POBkVc0SkVBgkYh8pqo/VjfY+ujqq51ukp5GrV9/hQ8/hKIi54YnY4zxlUoTtoj8HbgAWAcUuasV+K6i49QZPzTLXQx1H41uTNHOnZ2Hh2e0z9xcaNYsMDEZYxoHb2rkk4HeqppX1cLdiUISgaOAp1T1pzL2mQ5MB+jatWtVT1Hn/fCDM7/whRc6y02bOs+W4I0xvuZNG/xvOLXvKlPVIlWNBzoDw0SkXxn7zFHVIao6pE2bNtU5TZ32yivw5z8fWvbU4G1MeGOMr3lTg88GVorIfJx2dQBU9SZvT6Kq6SKyAJiAM4FIo5Gff3gPmlatoGdPKLZ+SMYYH/MmwX/oPqpERNoABW5ybwqcAvy9quXUd/n5h3rQgDM/69lHDN1mjDG1z5teNK+4Cbqrqm6sQtkdgFfcdvgg4C1V/biacdZbpWvwxhjjL96MRXMmsBL43F2OF5FKa/SqulpVB6nqAFXtp6r31zjaeqh0gl+xAkaPdgYgM8YYX/KmiWYWzuTaCwBUdaWI2CSMXvrPfw4fXOzAAVi4EGyOYmOMr3mT4AtVNaPU8DONrj97dXXpcvhyyW6SxhjjS950k1wrIhcDwSLSU0T+Dfzg47gajDffhPffP7Rs3SSNMf7iTYK/ETgGp4vk6ziTbt/sw5galMcfh2efPbRc8k5WY4zxJW960WQDM92HqaK8vMMvskZGwqBBNia8Mcb3KqzBi8hl7gxOB9zHMhG51F/BNQSle9G0bQvLl8PkyQELyRjTSJRbg3cT+c3AX4DlOBN/DAb+ISKo6qt+ibCes37wxphAqagGfz0wRVW/VdUMVU1X1W+Ac9xtxgv5+RAWdmhZFYYMObxd3hhjfKGiNvgoVd1SeqWqbhERa0H20uLFhw9VIAJr1kBSUuBiMsY0DhUl+Io68lknPy917HjkuvBw60VjjPG9ihL80SKyuoz1AnT3UTwNziOPwNChMHbsoXVNm1qCN8b4XoUJ3m9RNGCzZsHNNx+e4MPD7UYnY4zvlZvgVdVaiWtItexeNCee6IwJb4wxvlSVSbRNFRUVOUm+ZC8agP/+NzDxGGMaF2+GKjDVlJ/vPFs/eGNMIHgzHvwZImIfBNVQXoK/9FK44AL/x2OMaVy8aaK5EPiXiLwLvKSq630cU4MRHQ3p6Ucm+F27YN++gIRkjGlEKq2Zq+o0YBCwGXhJRBaLyHQRifR5dPWciJPkPWPAe1g3SWOMP3jV9KKq+4F3gTdw5lqdAiwXkRt9GFu9t3s33H77kdPzWTdJY4w/eNMGP0lE5gHfAKHAMFU9DRgI3Orj+Oq1tDT4xz9gY6mpyu1OVmOMP3jTBn8u8Liqfldypapmi8iVvgmrYfBcZC05Fg04g40ZY4yveZPgU0sndxH5u6reoarzfRRXg1BeL5obbvB/LMaYxsebNvhTy1h3Wm0H0hCVTPCPLHqEV1fZEPrGGP8pN8GLyHUisgboIyKrSzx+B8oahMyUUjLBP/nTk8z+fjYATzwBbdpAYWHgYjPGNHwVNdHMBT4DHgZmlFifqap7fRpVAzFmjJPE84vySF2YSmpWKrsO7KKgoC27dzsXWps3D3SUxpiGqqImGnUn/LgByCzxQERifB9awxAcDClZ2w4uf5f0HeHhzmvrSWOM8aWKEvxc9zkRWOY+J5ZYNpVITITrroPEzYcG5ly4ZaEleGOMX1Q0XPAZ7nM3/4XTsPzyizP3apdJToI/KuYoFiQtYLh7D7AleGOML5Wb4EVkcEUHqury2g+nYfFcZN2Vl4wgTO0/lfsW3kfrY3YzbVrrI4YwMMaY2lTRRdbHKtimwMm1HEuD40nwO3KSaN+8PeN6jOO+hfeR3eY7Xnvt7MAGZ4xp8Cpqohnjz0AaIk+CT8lOJrZFLEM6DiEiNIIFWxZw9tGNI8GvXAn9+kGITS1jjN9V1A/+ZPf57LIe/guxfgsPh+2ZScRGx9IkuAkndDmBz9YtpGlT+OqrQEfnW2+9BYMGwYoVgY7EmMapol40o9znM8t4nFFZwSLSRUS+FZH1IvKziPypxtHWMzfcAAeyi9mWlUzX6K4AjI4dza9Zq8mVPQ16RMkdO+D66+HYYyE+PtDRGNM4VdREc6/7fEU1yy4EblHV5e7Y8Yki8pWqrqtmefXSrgO7yC/KJzY6FoBRce7nZuz/yM2dHLjAfEgVpk+HrCx47TXIzoagIIi0GQSM8StvhgtuJSJPishyEUkUkX+JSKvKjlPVVE9PG1XNBNYDnWoecv3x1ltw5V+cLpKxLZwEP7TjUMKDm0LcgoZXg09IgLg4CAriyY/imHduAlFR0LYtvPxyoIMzpvHxZrCxN4A04BycoYPTgDerchIRicOZFeqnMrZNF5FlIrIsLS2tKsXWeYmJ8NUSJ8F7mmjCQsIY0u4EiFvQsPrBJyQ41fakJAQljiQmzJtOpwUJ9OwJb78d6ACNaXy8SfAxqvqAqv7uPv4GtPD2BCLSHGc2qJvdmaEOo6pzVHWIqg5p06aN14HXB/n5EByTDHCwiQZgdNwoaLeaTkfV/yF9UlJg7VrIu3Wm0xZTgmRnw8yZnH8+LFrk7GuM8R9vEvy3InKhiAS5j/OBT7wpXERCcZJ7gqq+V5NA66P8fJCWSUSHRRMdHn1w/bheo0GUok7/C1xwteCll6BTJ+jfH0J3JJe5jyYnMeWcQlTh3Xf9HKAxjVxF3SQzRWQ/cA3OuDT57uMN4M+VFSwiArwArFfVf9ZOuPVLfj5odNLB5hmPYZ2GER4Szre/LwxQZDWXne3MN3vccfDi6/vY0aJJmfslRcFrKX+lXz9rpjHG38pN8KoaqapR7nOQqoa4jyBVjfKi7BHAJcDJIrLSfUystcjrgWbNIKhF8sELrB5hIWHk/TaMt39cHKDIai4iwqmRP/DkVh7deyIzTimkMPzwJK8REXxx5Uge//Fxbvn7WubMCVCwxjRS3jTRICItRWSYiIz0PCo7RlUXqaqo6gBVjXcfn9Y85PrjiScgrF3SYe3vHiH7jyKDsps16rqiIue564AtXP6/49m2fxtXzP6KkOdfhNhYEIHYWGTOHM7527tEhUXxws7r6NW7OLCBG9PIVHoDuYhcDfwJ6AysBI4DFmNj0VRqf95+0nPTj2iiAWiS24ls2UFhcSEhQfXnPn5VOP10GDYMkgfPYm/OXhZftZiB7QdCN2Dq1MP2bw3MPmU2V390NTNef4Vmm67g3nsDEbkxjY83Nfg/AUOBJHd8mkE4XSVNJWY8cmQPGo+w/E6oFLMja4e/w6qRefPgiy8gpFUyCWsS+MPgPzjJvQJXDLqCE7qcwFObbmPW7D3Wm8YYP/Emweeqai6AiISp6gagt2/Dahh+XH/4TU4lRRR0BmD7/u1VL7jEDUXExTnLfpCTA7fc4gweltbTGWz0lhNuqfS4IAnimdOfIU/SYewM601jjJ94k+C3iUgL4H3gKxH5ALA6mBeygp0afFlNNBdMdG7q3bZ/2xHbKlTihiJUnefp0/2S5F98EbZsgfsfTeOFFc8xtf/UMn+2sgxoN4Cbj7sZjn2e5z6rvxeXjalPKk3wqjpFVdNVdRZwN07Xx8k+jqtByA5NQoqb0L55+yO2zbjercFnVrEGP/PIG4pwbyjytQULoFs3WB76JLmFudwx4o4qHT9r9Cya0YY1UbNJSqp8f2NMzXh1dc+d3elEnIk+vlfVfJ9G1UDkNEmiaX4XguTIz9Hw4laEBYdVvYkmuZyeN+Wtr0UPPAC/JmdyydL/MLnPZI5uc3SVjm/epDlTj7mCOcWPsWxjCrGxHX0UqTEGvBts7B7gFaAVTqeIl0TkLl8H1hBodDItpOwmjAsuEHR/R7ZlVrGJpms5TSLlra9FffrAhub/R3puOneeeGe1yrh1zNUQVMSG8JdqOTpjTGnetMFfBAxV1XvdIYSPA6ZWcowBIjokMW74kRdYAWJiICirc9Vr8A8+6NxldNiJIpz1PpSYCM+/nMdji//J2G5jGdppaLXK6dmqJyd3O5nnlj9Hxn7rF2+ML3mT4LcA4SWWw4DNPommAckvyiclM6XMLpLgJPiifZ2qfpF16lSYM+ewG4qYM+eI/ue1LSEBrnt6LjuyUplx4owalXVZv+kkZSRx/WMNfEorYwKsorFo/i0iTwJ5wM8i8rKIvASsBbL8FWB9tX3/dhRl/Y9lN53ExEDBns5sz9yOqlat8KlTKdi0hc8/LWZCny3knO37L1Q//giRw9+hR8sejO02tkZlXTBgMqH5rflg2xyq+qMbY7xX0UXWZe5zIjCvxPoFPoumAUnKcLqJZKeWX4MnsxO5hbnszdlLq4hK51A5TN++kJoKBw7Apk0wYEBNIy5fXh4sW32A4nHzuazXdTjjyFVfWEgYJ8dcwRfBj/PlD6mMH9GhliI1xpRU0WBjr3gewOs4iT4RmOuuMxVIznB6tcQElZ3gTzwRpk1y+sJXtaukqjO2+nHHOcvr11c/Tm+sXAkFnb+mSPI4s/eZtVLm/ZOvhuBC/vbJy7VSnjHmSN70ohkNbAKeAp4GfvFmsLHGLindqcHHhHQuc/ugQXD9tOrdzZqZ6XR9HzXKaYbfsKFmsVZm1Sqg90dEhkZxYtcTa6XMYT160Wr/GBbnPUdhkV1sNcYXvLnI+hgwTlVHqepIYDzwuG/Dqv+SMpIgsz0RTcLL3F5QAIX7qnc3a2qq89y9u3ON1dcJ/uo/FNNmxMec1msCTYLLHve9Om4+aTpFUb/z9W9f11qZxphDvEnwoaq60bOgqr8Aob4LqWFIzkgmoiCWzmVX4NmyBUYO6oAgVW6i8QzW1aEDnHQSRHkzOn8NLEtZRlrOTs7sVTvNMx63nTGF1hGteW7F/9VqucYYhzd3siaKyAvAa+7yVJy2eFOBpIwkTj9xINedV/b2mBigOJRIaVflGnyHDs6gX716wauv1jzWiuzYAZc99BFBbYM47ajTarXssJAwJnW5gpfW/ZPkfSl0bWl3thpTm7ypwV8L/AzchDN08Dp3nSmHqrI1Y2uFA3G1aOG0nzcv7lzlGnyfPvDoo5T77aA2LV4MG4o/YkDLEVXu6eONwTodlSLunvdCrZdtTGNXYYIXkSAgUVX/qapnuwOPPa6qeX6Kr17KyMsgpzCHN57rWO48pMHBTpIPy+9U5Yuse/Y43SMB1q1zhu/92kfN2F/8mAztV3F+fO02z3hcOfkogrecyru/z6GwuNAn5zCmsaowwatqMbBKpJwBVUyZUjKdRvLtGzqyd2/5+8XEQEhO1e9mvfFGGOjOsdG6Nfz8s/Pwha+TPwZgSl/fJPimTWFE2HUcCNnGB+sb1YyOxvicN000HXDuZJ0vIh96Hr4OrD7zJHgyO9Kkgk4nDz4II+M7sy93H9kF2eXvWLr8FKcdHqBNG2jZ0jc9aQoK4PcmH9Gi6Ch6t/LdHC9/Pv0M2N+Rh7981mfnMKYx8uYi630+j6KBSc10+zFWkuAvuADyVnXihS1OX/ierXp6V34qxMc7r0XcUR59kOB/25aFxn7DiW1vqPHdqxWZOD6U6BevZnnUA/y+73e6tezms3MZ05hUNBZNuIjcDJwH9MEZB36h5+GvAOujQzX4DhUm+K1bISul6hN/pKYeqsGD7xL8utyv0OB8/nK6b5pnPJo0gVUvXo2IMCdxjk/PZUxjUlETzSvAEGANcBrODU/GCymZKTQPjWLUCc1o1678/f7+d7jzhqrd7JSV5dzJWjLBn3wynHoqFNbyNcqPN31CdFh0rd29WpHYll04s9eZPL/8BfKLbD4ZY2pDRQm+r6pOU9X/A84FTvJTTPVeSlYKnaM7smABjKxgUIdWrWD/dnc8Gi970ojA44/DKaccWjdtmtMfPsSr+bm8o6q8+v3ndMw7hdBg/9zX9usb17I7J4156+dVvrMxplIVJfgCzwtVtf5rVZCamUrHyMpv2omJAfKbE9Uk2usmmmbN4Oab4dhjD1+vCvm1WPH9Zu06CiO2M6Dp+NortBIndhiHpHdj9qJHKVYbn8aYmqoowQ8Ukf3uIxMY4HktIvv9FWB9lJKZQmhuR3r1cmZCKk9MjPPctqn3XSXT0mDjxsObY4qLoX17uOeeGgRdyiuLvgBg2vH+S/AXXhCEfjuL5TuX8fqa1/12XmMaqoqGCw5W1Sj3EamqISVe+3j0k/pLVZ02eO3Apk0Vt4t7EnyrEO/vZk1IcC6q7i/xERsU5PSHr80LrQu3fYHs6cP44/x3C8RJJ0H7XdOIOjCYGfNnVKnrqDHmSN70gzdVsC93H3lFeUQHOU00FfWiGTIE3nkHerT1vgafkgJhYU7f95JqsydNTkEOW4O/o1POeEL9OKxccDDc+Mcg9r/1T7bt38bji23QUmNqwhJ8LfN0kYySyhN8u3ZwzjnQvXUndmTt8OpWfU8XydLd0vv0gc2bnZuTamrhlu/Q4FwuGOK/5hmPm26C1/42ism9p/DwoofZkbXD7zEY01BYgq9lnpucIqk8wRcWwldfQZPczhRrMTuzdlZavucu1qz8LB5Z9Ah7svcAToIvLHSSfE19+dsXhAWHcf+Vo2peWBU1b+70Cpp96t/JL8rn7m/u9nsMxjQUluBrmacG37N9B04/vfKx2seNg/VLve8Ln5oKHTvCc4nPcef8O5ny5hTyCvMYPhzuvNPpZVNTH6//guM7nUREaETNC6um797vSdedf+SFFS+waseqgMVhTH1mCb6WeRL8OeM78PHHVHijU0gIREcDGd7fzTp7tjPYWMKaBNo2a8v/kv/HlR9eSc+eykMPQZcuNYt/a8ZWNmWs45dP/d88U1JeHmx+4W6ah7Tkmo+vITMvM6DxGFMf+SzBi8iLIrJLRNb66hx1UUpmCi3CW3hd+42JgYI93tfgzzgD2vZdT2JqIneeeCcPj32YuWvmcs+397B/PyQn1yh8Ptv0JQAndQhsgr/qKoht15K2Pz3LspRlnPraqezL2RfQmOqNhASIi3O6V8XFOcvGN+r4e+3LGvzLwAQfll8npWSl0DGyI08/7TSlZGVVvH9MDGSltaZJcJNK72bdvx/mz4cXliYQJEFc2O9C7hhxB1cPupq//e9vHHv1S1xySc3if2fFF7C/I5OO61ezgmooLMzp17/5o/O4Pe4dVuxYwehXRnt1naJRS0iA6dMhKcm5+y0pyVmuY4mnQagH73Ut3tx+OFX9TkTifFV+XZWamUqH5h3I2O60l1fWzTAmBvbtFTpFdmJbZsU1+DVr4JRTi2n3UAKndj+V9s3bA/D06U+TlJHE1zqd9I0DgGMrLKc8RcVFfJ/6NWw+ixP/5LvRI7116aXw8MMw7+HJfPTRx0x5azIjXx7JV5d8VeFsWQ2VKuzcqXy/ahepeZto3X0baQfSePHNNHbn7CJHM1ie8AFds3MOPzA7m+Trrmbo0rcI1WY0kQi6d2nG8PjmxDSNYcl3LWkbGUOnmBiOat+evl070L1Tc8LDffc3kLAmgZnzZ5KckUzX6K48OPZBpvaf6rPzFRYq23ZlsTl1L9v37COu934ycjNYvHw/G37fT3ZhFtlFWeQUZVEUdIB+g7LJKcxh+Zpsdu7JQYPyUSlAg/IJbpJP565FzL9nM52yS/V8y84m+foruKRgDp0iO7H2h86kb+1E6/QJxBT3Jjwc+vZ1mloBHnjA6Thxzz2Hjy9VW3yW4BurlMwURsaOPDhsQGUJ/pFHnC6PN62sfGan1FSgyw/szNvCtAEPHFwfGhzKW+e9Rdu/d2R355c4cODYal1sXZqylGzdR8y+8XStA/kzJATefdfpGTTuqFP59KIvmZgwkUH/N4gr4q/gD4P/QO/WvhunPtDSc9NZkbqC2a8tZ3nqcvYEbaAoehOEudcjVrg7RgUREt6KJsUt6JyVU2ZZnTNzyQxKoigom+LgA6RkHmDB91kUaZGzQx6wG/jFWQzVCLrGdKB9s45sXdeJjpGd6N66Ez3adqB7u3b0i2tHtzbtaNm0JUFStYaAhDUJTP9o+sEb2ZIykpj+0XSASpN8XkEBKXszyA9KJ7MgnY1J+1i9aR+p6XvZlbmXPTl7Sc/bQ6+Be0jP38MvW/ewJ3svGrYPgksk48VlFB4EFDYjKK8Z+5IjaBrSlLzQCMIiwgkmjCCNJEibECahDGwfTIf0jWXG2Hm/01d5yfYlbIl8j6Jj8theHEbsxn/SYtN1REUd+uBcuBBWr4a//MX7968qAp7gRWQ6MB2ga13IKjXguYu1Y2RH8vKcBBVUyd/+4MHOc+ffOrNk+5IK901NBQb8l4iQCCb3mXzYthbhLYhvPp6lfd7nl01PMii+6q1vX/z6BYLw7G2nVvlYXxkwwHkAbP1+BNn/XkToWbN4PPtfPLb4MQa1HMm1x1/GqG4j6NmqZ5WTTV2xfTt8umg7H61ZwPI9C9gTuZDcZpsObm8S3ZkOcgzdok6gf8eenNC7J/HdutK2WVtimsYQHBTs7PhGnNNUUEpQbCzZ/1x52DpVJTMvi21797I5ZQ9bdu3h97QdbN27g5CWqRCZym9pKaSwlOTC9/lxTy7sAdYfKkMQyIsmpKAFIUUtCdNowqU5/XpF0q1TcwpzIvhtUxjhIWGEh4bRJDiUTzIeJlsPv0s5uyCbq+Zdx7MfLyY9O5NOPfZTGJTJlpRMknZmUBSagTbJgNCyP8AOaUowrQhLb0W7qFb0ih5Aek4MMUGtaNW0Je2iYmgb1YKThkUTExFNcGE0kU0iaRXZnKimEQRX9g9bUtefyn2vF16+8OB7vHX/Vq775Do+DbqBsyZ/yb8nvQA48xv7aqpNj4AneFWdA8wBGDJkiAY4nBrZk7OHguICOkZ2JDm/4j7wHhs2wJIlMKDbQN5Y+wY7snYcbHopbWtKPhzzFpP7TKF5k+ZHbD+z59kszXqfz1cvY1D8sCrH/8HGDxjeeTjnnVH7k2vXhgsugMzM/rzxxrss/2gHWUe9zIrBz3PNvqsACCMSUo4l+sAQhnYZxAUjB3P+2J6ENQkOcORlKyou4oetP3DLcx+wdP+H0GoThEJwTAs6F5/E9DFXcGzHwcS3H0S75m29K/TBB5124OwSCTQiwllfiogQFR5J346R9O0YW2GxubnKig37WL91B7+n7SSmy06CInfya8oevv5+H9nF6eToPvLYT7puY11GFmsPZLE/5wDZefkQkge5FYeeV5zJovTXIT+S7N2RtGsZSUxEC4qbxNI0KJpmRNNMo2kR1pJxI1vQpU0L8vZHU5gZQ7d2MfTo1JKYqPAjbgL0GS/eaxGha3RXPrroI5786Ulu/+p2Bj47kDfPfZMRXUf4PMSAJ/iGxNNFskPzDsTEw4UXVn7MZ585X88WbhwP3MmXm7/k0oGXlrnvsozPoN0+pg0o+2vsFSPOYNaKEJIi3gOqluB/3vUzK3as4JLWT5Dv5YeTv4WGwnXXOQ/V9mzZMoOVq26nx3E/k5i6jDcXLSMxbylp7Z/kk6B8PlkCl30fwfC4gQzrNJSorOH0ix7Oicd0p2NH8V8iKCEru5BH3vyGVxPfJLvLR+zJTSOkSSg9W5/MxN7XceFxoxnaZcChGnlVTXX/NmbOdLpUde3qJJypNWvfDg8Xjo+P4fj4GKDv4RvPLv+44mLIzYXsbGX/gQJy8goY9+HRpBzYesS+nSNjWfHHLcTEVP7Nt06ownsdJEHcfNzNjIwdyXlvn8f575zPlj9t8flQ3KLqm0qziLwOjAZaAzuBe1X1hYqOGTJkiC5btswn8fjDF79+wYSECSy6YpHXn86vvAKXXw6bfi1mxLwOjO02lrnnzC1z3/EvnMeSnQtJm5FCSFDZn82nvnYqyRnJbLhhQ5Wm2fvr/L/yyP9m0/SZ7WRsb1erY8v7W0FRAUu2rOP1b1ewOm0Fxe0SWZ66nJxC9+t9diuapp3IwOjRXHzCKK4/pwYJ1QtFxcU89/kP/HvB66yXt9GINCQ/ijGdTueaUZOZcNQEosIaz/h9pdvgASJCI5hz5hyfXmitKz7c+CFnvXEW753/HlOOnlLj8kQkUVWHlLXNl71oLvJV2XWVpwbvzVjwHp4RJdP3BTG+x3g++/UzirX4iLbk9Nx0FqZ+xDXHXlNucgcY2eZs7vntetbvXk/fNn3L3a+kYi1m7pq5hG0/lTFD63dyB+ei84geAxnRYyBwOQCFxYV8umwtX6//iZ9SfuTn4O/4MewDflwH9/zWglGxo2ieNparTx7LqL5H13gO2qLiIr7d/D0f/voOb//8LjsOpECTcLrmnsnlvS/i9imn0SwsvOY/bD3kSeL+7EVTl0zsOZHOUZ15ZtkztZLgK1LP/5XrloNNNJEduOoqWLYMVlVyl70nwe/dC+N7jOe11a+xPHU5Qzoe/oH81s9vkVeUxwmR0yosb917Z0Gn63lv/XteJ/jvk78nKSMJfnqQ827x6pB6JyQohEnD4pk0LB64BoDf927l8/ULSdyzgC9++YZtBz4g4R0IzWtP7/ATObHHUKYMHcbxsccSGRYJlN+9r1iLWZe2jk/XLuKdJYtYuX8+BWE7CA8JZ8JRE+hVeC43jZ9Ep9aRgXsT6pCp/ac2moReWkhQCNMHT+eeBffw695fOSrmKN+dy2clN0IpmSm0DG9JeEg4WVnezbBUMsGfOsLpvfL5r58fluBVlWeWPgs7BrDxmyFQwRSp8T068sam43jn53ncNfIur+L+7+r/EqrN4LfJnHWWV4c0CN1iunDdiGmA86H55dLfeerTb1i4cz5rI39krb7Ds786PUXaNe1MXk4I6SSjOF0LkzKSuHTepcycP5PdWRkcKEp3Cs5sT+TeUZzTfQr/d8tEosItqZvDXTX4Ku5beB9zEucw+9TZPjuPJfhalJp1aKo+by9Udu/uzPrUowdEN2vLsR2O5YvNXxyWnJemLGXlzhWw7Gk6XlNx00HPnsCHZ7Oqy+0kpScR26Li3hF5hXm8te4tWqVNZujJzWjRovKYG6pxQ7sxbuhVwFWkpcFnC3fTtMcy1u9fwmc//saPBW9CaNFhxxRrMbsO7GJQ8KX8tugEzh12In+4tBv9+wfmIq6pHzpGdmRyn8m8uOJF7h9zP+Ehvmmuqw/XqusNTx948D7Bh4U5feGjo53l8T3Gs3jrYjJyMw7u8+yyZ4kIbg6rp9Gxkub9nj2BDU673rwNlU9e/emmT0nPTefFm6cxt+xru41SmzZw6bmtOW/QBO4ZdQ+L73gZCc0rc9/cwlwW3vYsKZ9dyr/v786AAZbcTeWuHXIte3L28O66d312Dkvwtag6CR6cnjSeGx7GHzWeIi1i/u/zAdiXs4831r7BidFTIT+y0tuZe/QA9h5FW/p5leD/u+a/tG3WllN7nELzI7vWmxLKGx6ha3RXQkKOnITFmIqc3O1kesb05Jllz/jsHJbga0mxFh/WRDNxIkye7N2x994Lr73mvD6+8/FENonki1+dSa9fXfUqOYU5HJN9HQCxFbe4EBEBn38OFw86m0XJi9h1YFe5++7L2cfHv3xM4YqLeOYpa62rzINjHzxilNCI0AgeHHvkTUTGVCZIgrjm2Gv4fuv3rNm5xjfn8EmpjdDu7N0UFhfSoblTxf7zn+G227w7NibGucgKThe/sd3H8vnmz1FVnk18luM6H8d91w/kxx8PXZStyPjxcNmwKRRrMQmryx/Z7p1175BflM/eBdNo3dq7WBuzqf2nMufMOcRGxyIIsdGxjabvtvGNy+MvJyw4jGeXPeuT8i3B1xLPVH2eGnxhoTP6nzdKJnhw2uGTM5J5dtmzbNi9geuGXEdkJAwf7l15P/8MP74/kFO6n8LMb2ayLm3dEfvkFOTw9LKnaVHUi7C9x3LGGd6V3dhN7T+VLTdvofjeYrbcvMWSu6mRVhGtOP+Y85m7di55hWVf46kJS/C1pPRNTvHxcO653h1bVoIHuPWrW2kZ3pIpvc7j7rsr71Pv8fnncN11wpNjXqV5k+Zc+M6F5BQcGqSpqLiIS+ZdwsodK2H+w5w+UYi0nnzGBMR9o+9j5TUrCQsJq/WyLcHXktIJvioXWUsn+G4tu9GrVS+yC7K5Iv4Ktv7elL/9zRlW1Bs9ezrPmSkdeHXKq6zZtYZbvjx0B9OtX97Ku+vf5Yaj/kn6D2dz/vnelWuMqX3dWnartDtzdVmCryWeBO8ZCTI/3+kC6Y0HHjiydn7aUacBcM2Qa1jjXn/xDJtbGU+C/+UXmHDUBG45/haeWfYM761/j8cXP84TPz3BzcNv5ubhf+Yvf4HTT/euXGNM/WJdJ2pJSmYKrZq2Ovg1q6DA+xp8mzZHrrtr5F2c0esMerXqxSurITgY+vTxrrzu3Z3R+Da5w4k/NPYhFiYt5NJ5l3Kg4ADn9j2Xx8Y/RpDAY495V6Yxpv6xGnwtKdlFEqrWRLNuHcyaBbt3H1rXOqI1p3Q/BXCm6uvTx/tvBGFhzsilngTfJLgJb5zzBsFBwYzoMoLXprzGP2YH8dNP3pVnjKmfrAZfS0re5AROF8nOnb07dtMmuO8+mDSJMrsr/v479O9ftXi+/hratTu03COmB7/e+CstwluwIjGUGTPgjju875ljjKl/LMHXkpTMFI5pe8zB5dtv9/7YkgOOlWX1ajhwoGrx9Ohx5Lo2zdpQXAw33uhM8DtzZtXKNMbUL9ZEUwuKiovYkbWDjs2dGnxaGuzZ4/3xrdwZ8so7RoQqDyOwdq1TQ/d8aCSsSSDuiTiCHwhiyQlxTJqZYF0jjWngLMHXgl/2/EKRFtEjxqk2/+tfTvNInpf3LVRUg3/vPbjiisOnffRGUhLMnu3M+eqZQScpIwlQaJHEaxnTSVhT/l2uxpj6zxJ8LViwZQEAo2JHAU6Xx6pcFG3Z0nkuK8F//TXMmwdNm1YtJk9XyX/8A254b+Zh06OBM4v9zPnWRmNMQ2Zt8LVgQdICOkd1pnvL7oDTZj6iChOmh4XBvn0QVca0nKtXOxdYqzpSYbduMHAgfPop5A9MhjKOT85Irlqhxph6xWrwNaSqLNiygNFxoxER9u1zJlgfOLBq5bRoceRM8qpOF0lvb3AqKTQUVq50mokqGubWGNNwWYKvofW717PrwC7GxI0BDg0nUNUE/+9/w1NPHb4uORn27696F8nSHjrFhrk1pjGyBF9Dnvb30XGjAejdG55/HoYNq1o5778Pb7xx+Lo9e5y2/Kp+WJRmw9wa0zhZG3wNLdiygC5RXejWohsA7dvDVVdVvZyYGOeO1pIGD4b162shSBr3LPbGNFZWg6+B0u3vAJ98Aps3V70sz9ACP/xQy0EaYxotS/A1sH73etKy0w42zxQWwjnnwNNPV72sv/7VmY5v8mSn7R1g7FhnjBpjjKkOS/A18O3v3wKH2t9/+cXptVKdNvNWreDjj+G88+Db3QnEPh7HNycF8ThxdkOSMaZarA2+BhYkHd7+Xt0eNB69e8MJ1zp3nWYXZIPAfpKY/tF0AGtDN8ZUidXgq8nT/j6m25iD7e+rVjn9z48+uvrlzpxvd50aY2qHJfhqWpe2jt3ZuxkdO/rgulWrnOTu7TjwZSnv7lK769QYU1XWRFNNpfu/A7zwgjOSZE10je7qDgp25HpjjKkKq8FX04KkBXSN7kpci7iD6zp0qN6wAiU9ONbuOjXG1A5L8NVQVv/3tWud4XmrMg58WeyuU2NMbbEmmirKL8pnxtcz2J29m7Hdxh5c/9VXzgQbV1xR83PYXafGmNrg0xq8iEwQkY0i8quIzPDlufzht32/ceKLJ/L4j4/zx6F/5KJ+Fx3ctmqVM0xBmzYBDNAYY0rwWQ1eRIKBp4BTgW3AUhH5UFXXVXxk3VKsxezP289nmz7j2k+uJUiCePf8dzn76LMP22/lypoPCmaMMbXJl000w4BfVfU3ABF5AzgLqPUEf8y/hrA5KeeI9e3aQfNIyMmBlO1HHtehA0Q0cya0Tk1VEAUUAA0qIKJlOpmF6RRrMQBNdx9Hpx9f5+5X47gbeOcdp1vkm28647aPH1/bP5kxxlSfLxN8J2BrieVtwPDSO4nIdGA6QNeu1esKeFSLPmQkHTkBavcoaN0aMjIgv4wE3yPamS5vbzAUpoKo4Jn6SApDOLFrS47qFMO+lJYs/aYjnfafTVBc6MHjw8Od55Yt4fzz4bLLqhW+Mcb4hKiqbwoWOQ8Yr6pXu8uXAMNU9cbyjhkyZIguW7bMJ/EYY0xDJCKJqjqkrG2+vMi6DehSYrkzkOLD8xljjCnBlwl+KdBTRLqJSBPgQuBDH57PGGNMCT5rg1fVQhH5I/AFEAy8qKo/++p8xhhjDufTG51U9VPgU1+ewxhjTNlsqAJjjGmgLMEbY0wDZQneGGMaKEvwxhjTQPnsRqfqEJE04MjZLirWGtjtg3Bqk8VYOyzG2lEfYoT6EWddiDFWVcsc5rBOJfjqEJFl5d3FVVdYjLXDYqwd9SFGqB9x1vUYrYnGGGMaKEvwxhjTQDWEBD8n0AF4wWKsHRZj7agPMUL9iLNOx1jv2+CNMcaUrSHU4I0xxpTBErwxxjRQ9TbB1/UJvUWki4h8KyLrReRnEflToGMqj4gEi8gKEfk40LGUR0RaiMg7IrLBfU+PD3RMpYnIn93f9VoReV1EwutATC+KyC4RWVtiXYyIfCUim9znlnUwxn+4v+vVIjJPRFoEMMQyYyyx7VYRURFpHYjYKlIvE3yJCb1PA/oCF4lI38BGdYRC4BZVPRo4DrihDsbo8SdgfaCDqMS/gM9VtQ8wkDoWr4h0Am4ChqhqP5whsi8MbFQAvAxMKLVuBjBfVXsC893lQHqZI2P8CuinqgOAX4A7/R1UKS9zZIyISBfgVCDZ3wF5o14meEpM6K2q+YBnQu86Q1VTVXW5+zoTJyF1CmxURxKRzsDpwPOBjqU8IhIFjAReAFDVfFVND2hQZQsBmopICBBBHZjBTFW/A/aWWn0W8Ir7+hVgsj9jKq2sGFX1S1UtdBd/xJkRLmDKeR8BHgduB+pkb5X6muDLmtC7ziVPDxGJAwYBPwU4lLI8gfMHWhzgOCrSHUgDXnKbkp4XkWaBDqokVd0OPIpTk0sFMlT1y8BGVa52qpoKTkUEaBvgeCpzJfBZoIMoTUQmAdtVdVWgYylPfU3wUsa6OvkJKiLNgXeBm1V1f6DjKUlEzgB2qWpioGOpRAgwGHhGVQcBBwh8s8Jh3Hbss4BuQEegmYhMC2xU9Z+IzMRp7kwIdCwliUgEMBO4J9CxVKS+Jvh6MaG3iITiJPcEVX0v0PGUYQQwSUS24DRznSwi/w1sSGXaBmxTVc83oHdwEn5dcgrwu6qmqWoB8B5wQoBjKs9OEekA4D7vCnA8ZRKRy4AzgKla927Y6YHzYb7K/f/pDCwXkfYBjaqU+prg6/yE3iIiOG3G61X1n4GOpyyqeqeqdlbVOJz38BtVrXO1TlXdAWwVkd7uqrHAugCGVJZk4DgRiXB/92OpYxeCS/gQuMx9fRnwQQBjKZOITADuACapanag4ylNVdeoaltVjXP/f7YBg92/1TqjXiZ49+KLZ0Lv9cBbdXBC7xHAJTi14pXuY2Kgg6rHbgQSRGQ1EA88FNhwDud+u3gHWA6swfnfCvht7CLyOrAY6C0i20TkKuAR4FQR2YTTA+SROhjjf4BI4Cv3f+fZOhhjnWdDFRhjTANVL2vwxhhjKmcJ3hhjGihL8MYY00BZgjfGmAbKErwxxjRQluCN34hIkdvlba2IfFTVEQJFZIGIDHFff1rZ8SKypawR/tz1a9zHOhH5m4iEuds6isg7FZTZQkSur0rctUkc37jj83jWTXFHM+xTYl0bEfk8MFGausISvPGnHFWNd0db3AvcUN2CVHViDQccG6Oq/XEGruuO22ddVVNU9dwKjmsBBCzBAxOBVaWGvbgIWESJ0StVNQ1IFZERfo7P1CGW4E2gLMYdIE5EhonID+5AYj947lgVkaYi8oY7JvibQFPPwSVr5yLyvogkumOxT69KEKqaBVwLTHbHSY/zjPktIseIyBL3W8dqEemJc1NQD3fdP0SkuYjMF5Hl7jeCs9xj48QZt/45N64vRaSpu+0oEflaRFa5x/Vw198mIkvdc91XTshTKXHnqTvW0QjgKo4cnvh9d3/TWKmqPezhlweQ5T4HA28DE9zlKCDEfX0K8K77+i/Ai+7rATiDTg1xl7cArd3XMe5zU2At0Kr0PqXiOGI9sBIYDsQBa911/8YZBwWgiVv+we3u+hAgyn3dGvgVZzC8ODfeeHfbW8A09/VPwBT3dTjO0MLjcL5FCE7F62NgZBmxJwGRJZanAS+4r3/AuV3es60TsCbQv3d7BO4R4uXngDG1oamIrMRJfok4kzoARAOvuDVkBULd9SOBJwFUdbU7TEFZbhKRKe7rLkBPYE8VYytrhNLFwExxxsx/T1U3OcPMHHHcQyIyEmfI5U5AO3fb76q60n2dCMSJSCTQSVXnuT9XLoCIjMNJ8ivc/Zu7P8d3pc4Xo878Ah4X4Qz5DM6AcRfhDJcAziBiHSv+sU1DZk00xp9yVDUeiMWpEXva4B8AvlWnbf5MnFqtR4VjaYjIaJxa//GqOhAnQVZpqjw36cbhzBx06MSqc4FJQA7whYicXMbhU4E2wLHuz7azxPnzSuxXhFPbL+uDBHf9w+pco4hX1aNU9YUy9isUkSA37lbAycDz7oiGtwEXyKFPoXA3dtNIWYI3fqeqGTjT290qzpDK0cB2d/PlJXb9DrcNWUT64TTTlBYN7FPVbLcXyXFVicVtw34aeF9V95Xa1h34TVWfxBmBcQCQiTMIVsnz71LVAhEZg/PhVS51Lo5uE5HJ7jnCxBlb/AvgSjceRKSTiJQ1EcdGnIvCAOcCr6pqrDqjGnYBfgdOdLf3wmmyMo2UJXgTEKq6AliFc2FwNvCwiHyP0z7v8QzQ3G2auR1YUkZRnwMh7j4P4Ezv5o1v3YupS3CG+r2mjH0uANa6zUp9cJLpHuB7t6vnP3AmohgiIstwPow2eHHuS3CalVbjtJu3V2f2p7nAYhFZgzMyZWQZx34CjHZfXwTMK7X9XeBi9/UYd3/TSNloksbUI+JM0PGqqp7qxb7fAWeV/mZiGg+rwRtTj6gzh+pzJW90KouItAH+acm9cbMavDHGNFBWgzfGmAbKErwxxjRQluCNMaaBsgRvjDENlCV4Y4xpoP4f0uyk7qE2p+IAAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# example rdf plot\n", "\n", "easy_rdf_bins = npz_loader(datafiles.easy_rdf_bins[\"bn_all\"])\n", "easy_rdf_data = npz_loader(datafiles.easy_rdf_data[\"bn_all\"])\n", "\n", "#find the cutoff point of the RDF\n", "cutoff = identify_cutoff_scipy(easy_rdf_bins, easy_rdf_data)\n", "print(f\"solvation radius is {cutoff}\")\n", "\n", "# plot the result to check\n", "fig, ax = plot_scipy_find_peaks_troughs(easy_rdf_bins, easy_rdf_data)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This approach works well for RDFs with well defined minima and is the default minima-identifying method used in `solvation_analysis`. Lets now confirm this by going though a typical calculation of solution properties." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning: importing 'simtk.openmm' is deprecated. Import 'openmm' instead.\n" ] }, { "data": { "text/plain": "{'PF6': 2.8500000000000005, 'BN': 2.6500000000000004, 'FEC': 2.75}" }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "# instantiate Universe\n", "u = mda.Universe(datafiles.bn_fec_data, datafiles.bn_fec_dcd_wrap)\n", "\n", "# define solute AtomGroup\n", "li_atoms = u.atoms.select_atoms(\"type 22\")\n", "\n", "# define solvent AtomGroups\n", "PF6 = u.atoms.select_atoms(\"byres type 21\")\n", "BN = u.atoms.select_atoms(\"byres type 5\")\n", "FEC = u.atoms.select_atoms(\"byres type 19\")\n", "\n", "solute = Solute.from_atoms(li_atoms, {'PF6': PF6, 'BN': BN, 'FEC': FEC}, solute_name=\"Li\")\n", "# run analysis\n", "solute.run()\n", "solute.radii" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here the slight difference between the precomputed **BN** radii (2.7) and those we have generated using `Solute` (2.65) is due to the small number of frames analyzed and can be ignored." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using a Custom Function to Identify Minima\n", "\n", "Due to the roughness of some RDF data, the univariate spline interpolation and critical point finding algorithms may sometimes fail to find the correct minima (if there even is one). One can either set the radii manually as is done for **PF6** in the basics tutorial, or you can supply a custom function of arbitrary complexity. \n", "\n", "This custom function must take the rdf bins and rdf data return a single scalar value (the minima), e.g. `foo(numpy.ndarray, numpy.ndarray) -> float`. Lets try constructing a custom cutoff finder based on Scipy's signal processing functionality. First, we'll see how the function does as a standalone operation, then, we'll integrate it into our `Solute`.\n", "\n", "We will use `scipy.signal.find_peaks` on the **negated data** to find the minima. We will also require that the peaks have a distance of **5 samples** between them with the `distance=5` kwarg." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# find peaks using scipy\n", "peaks, _ = find_peaks(-easy_rdf_data, distance=5)\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "solvation radius is 2.7\n" ] }, { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7l0lEQVR4nO3dd5xTZdbA8d9hKMNQBUZFygwoivQmYgXEigrY9cW6Krq6uq5rYeW1r10XdVd3X9YCCiqI4tobiIigFKUjKEsRQSkiRdoMnPeP58YJM0nmJuQmk8z5fj75JLn35t4zmZmTJ8997nlEVTHGGJN9qqQ7AGOMMcGwBG+MMVnKErwxxmQpS/DGGJOlLMEbY0yWsgRvjDFZyhK8KZeIbBGRlsneNh1E5DYReSbdcQCISHPv/cpJ4j6PEpFvvf0OEJH3ROSSBPd1l4iMTFZsCRx/uIj8NV3HzwaW4CsAEVkmItu8f8ofvT/s2mHrh4vIThHZ7N3micgDIlIvbJtLRWSXt4/Q7R/JiE9Va6vqf5O9bTqo6v2qekW64wBQ1RXe+7Uribu9B/iHt983VPUUVR2RxP2bDGIJvuI4XVVrA52AzsBfSq1/WFXrAPnAZUAP4HMRqRW2zVTvHzt0+0MqAjcVSgEwP91BmIrBEnwFo6o/Ah/gEn2k9dtVdTrQD2iIS/ZxEZGJIvJXEZnitfTfEpGGIjJKRDaJyHQRKQzbXkXkIO/xcBF5SkTe8b5NfCkiB8bY9mmvm2CLiHwuIvuLyOMiskFEvhGRzmGvHSwiS7z9LhCRM2L8DHt8fReRXiKyMuz5rSLyg7evRSLSx1v+W7eDiBR68V4iIitEZJ2IDAnbR00RGeHFulBEbgk/Rql4RESGisgaEdkoInNEpF3Yfh4TkeXeusnestDxq4b9Xh4QkWnedv8RkQbeundE5LpSx5wjIgPCni8BWgJvee93DW+fV3jrL/WO/aj3My0VkVPCXt9CRD713rOPgEbR3v+w14R+hkEiskpEVovIn8PWVwn7va4XkTGhn8lb/6q4b60bRWSSiLSNcpw6IvKJiDzpvdd9vb+Rzd7v+abyYq2MLMFXMCLSFDgF+C7Wdqq6GfgIOCbBQ50PXAQ0AQ4EpgLPAw2AhcCdMV57AXA3sI8X530xtj0X+F9cstjhHecr7/lY4G9h2y7B/Tz1vP2PFJHGcf5ciMghwB+Aw7xvPScBy2K85GjgEKAPcIeIHOotvxMoxCXNE4ALY+zjROBY4GCgPnAesN5b9yjQFTgS9/7eAuyOsp+Lgd8BBwDFwJPe8hHhxxeRjrjf3buhZap6ILAC79ugqu6IsP/DgUW49/9h4FkREW/dS8BMb929QDx9972BVrj3YbCIHO8tvx4YAPT0fqYNwFNhr3vPe92+uL+LUaV3LCINgfHA56p6vbr6Ks8CV3m/33bAhDhirTQswVccb4jIZuB7YA2xE2zIKlzCCOkhIr+E3XrEeO3zqrpEVTfi/smWqOrHqloMvIrrJormdVWd5m07iijfNjzjVHWmqm4HxgHbVfUFr995dPhxVPVVVV2lqrtVdTTwLdA9xr6j2QXUANqISDVVXaaqS2Jsf7eqblPV2cBsoKO3/FzgflXdoKorKUm2kRQBdYDWgKjqQlVdLSJVcAn7j6r6g6ruUtUpUZIvwIuqOk9VfwVuB84VdxL2P0ArEWnlbXcRMFpVd5b7buxpuar+23v/RwCNgf1EpDlwGHC7qu5Q1UnAW3Hs925V/VVV5+IaChd4y68ChqjqSu9nvgs4O/StRVWfU9XNYes6Sti5JdyHwqfAq6r6v2HLi3C/37re7+er+N6GysESfMUxwGuN9MIliXK/HuNacD+HPf9CVeuH3b6I8dqfwh5vi/C8NtH9GPZ4aznb+j6OiFwsIrNCH1C4lpmf92EPqvodcAMuYawRkVdE5IAYL4n28xyA+8ANCX9c+pgTgH/gWqc/icgwEanrxZ+L+3biR/gxlgPVgEZeAhwDXOh9aFwAvOhzn+F++1lVdav3sDZe69r7YAk/vl+l4w693wXAuLDf6ULcB/B+IpIjIg963TebKPmWFf47PxWoCfyr1PHOAvoCy71upSPiiLXSsARfwajqp8Bw3Nf6qMSNsjke+CwFYQVORAqAf+O6Vhqqan1gHiBRXvIrkBf2fP/wlar6kqoejUswCjyUQFirgaZhz5vF2lhVn1TVrkBbXFfNzcA6YDuuG8yP8GM0x7VU13nPRwADcV1JW1V1qs99+rEa2Ef2PGnfPI7Xl457lff4e+CUUg2PXFX9AfgfoD/u77gerjsM9vyd/xt4H3g3PDZVna6q/XFdO2/gPvxMKZbgK6bHgRNEpFPpFd6Js664P+oNuK/D2aAWLhGvBRCRy3At+GhmAX1FpIGI7I9rseO99hAROU5EauCS6zZcqzFeY4C/iMg+ItIE9+ETkYgcJiKHi0g13IfPdmCXqu4GngP+JiIHeK3WI7zYIrlQRNqISB5uyOPY0DBKL6HvBh4jsdZ7VKq6HJgB3C0i1UXkaOD0OHZxu4jkeSdJL8N1v4Fred/nfYAjIvki0t9bVwd3XmY97sP6/ij7/gPuvMHb4k5OVxeRgSJST1WLgE0k9vvNepbgKyBVXQu8gOuDDbnF66P/2Vs3Eziy1FfqjKWqC3CJayquG6c98HmMl7yI6y9fBnxISUIB1//+IK7l+yOulXdbAmHdA6wElgIf404KR+s7r4trbW7AdVGsp+Rb2E3AXGA67vf3ENH/917EfYP7Ede1c32p9S/g3pvQSKB/iUjp7otE/Q/uJOzPuHNAL8Tx2k9xJ9zHA4+q6ofe8ieAN4EPvb/fL7xj4O1/OfADsMBbV4Z3UnUQ7tvAf3Dvy0XAMq9r52pinwCvtMQm/DDGHxH5PXC+qvYMaP8TgZGqGvVKWxG5GBjkdT+lnbjhtEuBat5Jd1OBWAvemChEpLG4S/+reEMv/4wbCZSuePKAa4Bh6YrBZBZL8MZEVx34P2Azbpz1f4Cn0xGIiJyEOz/xE268eiqPPVD2LIERutkVsxVcoF00IrIM98+xCyhW1W6BHcwYY8weqqbgGL1VdV35mxljjEmmVCR43xo1aqSFhYXpDsMYYzLGzJkz16lqfqR1QSd4xQ2PUuD/VLXMySERGYQbAkXz5s2ZMWNGwCEZY0z2EJGoVxwHfZL1KFXtgiueda2IHFt6A1UdpqrdVLVbfn7EDyFjjDEJCDTBq+oq734NbnhZIoWjjDHGJCCwBC8itUSkTugxrozovKCOZ4wxZk9B9sHvh6siFzrOS6r6foDHM8YkQVFREStXrmT79u3pDsWEyc3NpWnTplSrVs33awJL8Orm5exY7obGmApl5cqV1KlTh8LCQkrmAjHppKqsX7+elStX0qJFC9+vsytZjTF72L59Ow0bNrTkXoGICA0bNoz7W5UleGNMGZbcK55EfieW4FNk8mSYZ6eYjTEpZAk+RY45Bg47LN1RGJP57rrrLh591JXa/+abb+jUqROdO3dmyRK/syJWHpbgU+Sww6BXr3RHYUzmUVV2794dcd0bb7xB//79+frrrznwQL+zIlYeFaoWTTbLzQUbdWaMP8uWLeOUU06hd+/eTJ06lQEDBvDSSy/RrFkz8vPz6dq1K++++y6PP/44OTk5TJo0iU8++STdYVc4luBT5LOsmBrbVEaRvnmeey5ccw1s3Qp9+5Zdf+ml7rZuHZx99p7rJk70d9xFixbx/PPPc/nll3PppZfy9ddfU1xcTJcuXejatSt9+/bl6quvpnbt2tx0003x/VCVhHXRGGMqpIKCAnr06MFnn33GGWecQV5eHnXr1qVfv37pDi1jWAs+hdq1S3cExsQvVos7Ly/2+kaN/LfYS6tVq9Zvj23YZmKsBZ8i1arB6aenOwpjMs+xxx7LuHHj2LZtG5s3b+att95Kd0gZw1rwKaAKRUUwdWq6IzEm83Tp0oXzzjuPTp06UVBQwDHHHJPukDJGoHOyxqtbt26ajRN+7NoFVb2P0gr0dhsT0cKFCzn00EPTHYaJINLvRkRmRpvv2rpoUiAnBwYPhho10h2JMaYysQSfIrm5sGMHRLlewxhjks4SfAr89BPcdZd7vGNHWkMxxlQiluBTIDyp29WsxphUsQSfAjt3uvuBAyFsaK8xxgTKEnwKhBJ8v35QvXp6YzHGVB6W4FMglODfeQd++SWtoRiTNd58800efPBBwFWVXLBgwW/r7rjjDj7++OO497ls2TLaJemS8759+/JLOf/ww4cPZ9WqVUk5XiSW4FNgn32gSRN44QWwktUm24yaO4rCxwupcncVCh8vZNTcUYEfs7i4mH79+jF48GCgbIK/5557OP744wOPI5Z3332X+vXrx9zGEnwWKCiA555zj+0kq8kmo+aOYtBbg1i+cTmKsnzjcga9NWivk/wLL7xAhw4d6NixIxdddBEAl156KTfeeCO9e/fm1ltvZfjw4fzhD39gypQpvPnmm9x888106tSJJUuWcOmllzJ27FgApk+fzpFHHknHjh3p3r07mzdvZtmyZRxzzDF06dKFLl26MGXKlJjxTJw4kWOPPZYzzjiDNm3acPXVV/9Wo/7ll1+mffv2tGvXjltvvfW31xQWFrJu3TqWLVvGoYceypVXXknbtm058cQT2bZtG2PHjmXGjBkMHDiQTp06sW3bNgYPHkybNm3o0KFDcipkqmqFuXXt2lWzUXGx6vjxqqD68cfpjsaY2BYsWOB724KhBcpdlLkVDC1I+Pjz5s3Tgw8+WNeuXauqquvXr1dV1UsuuURPPfVULS4uVlXV559/Xq+99trf1r366qu/7SP0fMeOHdqiRQudNm2aqqpu3LhRi4qK9Ndff9Vt27apqurixYs1lHuWLl2qbdu2LRPTJ598ojVq1NAlS5ZocXGxHn/88frqq6/qDz/8oM2aNdM1a9ZoUVGR9u7dW8eNG+fem4ICXbt2rS5dulRzcnL066+/VlXVc845R1988UVVVe3Zs6dOnz79t5/z4IMP1t27d6uq6oYNG8rEEel3A8zQKDnVWvAp8NFH0KePe2wteJNNVmxcEddyPyZMmMDZZ59No0aNAGjQoMFv68455xxycnJ872vRokU0btyYw7z5MuvWrUvVqlUpKiriyiuvpH379pxzzjl7dO9E0717d1q2bElOTg4XXHABkydPZvr06fTq1Yv8/HyqVq3KwIEDmTRpUpnXtmjRgk6dOgHQtWtXli1bVmabunXrkpubyxVXXMHrr79OXl6e758zGkvwKRA6yQqwbVv64jAm2ZrXax7Xcj9UNWp54FpxjjOOtq+hQ4ey3377MXv2bGbMmMHO8H/SKErvR0RQn8WlaoTVKcnJyaG4uLjMNlWrVmXatGmcddZZvPHGG5x88sm+9h2LJfgUCP3t/POfNi+ryS739bmPvGp7tjTzquVxX5/7Et5nnz59GDNmDOvXrwfg559/Lvc1derUYfPmzWWWt27dmlWrVjF9+nQANm/eTHFxMRs3bqRx48ZUqVKFF198kV27dpV7jGnTprF06VJ2797N6NGjOfroozn88MP59NNPWbduHbt27eLll1+mZ8+evn/W8Li3bNnCxo0b6du3L48//jizZs3yvZ9oLMGnQCjBH3ecmwDBmGwxsP1Ahp0+jIJ6BQhCQb0Chp0+jIHtBya8z7Zt2zJkyBB69uxJx44dufHGG8t9zfnnn88jjzxC586dWRI2VK169eqMHj2a6667jo4dO3LCCSewfft2rrnmGkaMGEGPHj1YvHixr28GRxxxBIMHD6Zdu3a0aNGCM844g8aNG/PAAw/Qu3dvOnbsSJcuXejfv7/vn/XSSy/l6quvplOnTmzevJnTTjuNDh060LNnT4YOHep7P9FYueAUGD4cLrsM7rkHzjwT2rZNd0TGRGflgsuaOHEijz76KG+//XZa44i3XLBN+JECbdvC738Pd9zhqkpagjfGpIJ10aTAYYfBk0+6xzaKxpjM06tXr7S33hNhLfgU2LrV9cPn5NgoGpMZYo1kMemRSHe6teBT4MknXbkCEWvBm4ovNzeX9evXJ5RQTDBUlfXr15ObmxvX66wFnwKhUTR161qCNxVf06ZNWblyJWvXrk13KCZMbm4uTZs2jes1gSd4EckBZgA/qOppQR+vIgp1z0ya5FryxlRk1apVo0WLFukOwyRBKlrwfwQWAnVTcKwKaedOVwfeRs8YY1Ip0D54EWkKnAo8E+RxKrpQgh8zBt58M93RGGMqi6Bb8I8DtwB1Aj5Ohda3LzRtCo89Bg0auJmdjDEmaIG14EXkNGCNqs4sZ7tBIjJDRGZk60mdE0+Em25yFznZMEljTKoE2UVzFNBPRJYBrwDHicjI0hup6jBV7aaq3fLz8wMMJ33WrIHVq12Ct1E0xphUCSzBq+pfVLWpqhYC5wMTVPXCoI5XkV1/PfTuDTVrWoI3xqSOXeiUAqGTrNZFY4xJpZRc6KSqE4GJqThWRRRK8E88AT7KThtjTFLYlawpEErw++2X7kiMMZWJddGkQCjBf/IJPPRQuqMxxlQWluBT4Lrr3O3DD+H229MdjTGmsrAumhQ46yx3P38+FBW5fvg4JoY3xpiEWIJPgUWLoFYtN4oG3FDJOCeHN8aYuFmCT4H+/aFTJzjySPfcErwxJhWsDz4FwsfBg13sZIxJDUvwKRBK8BddBD//DI0bpzsiY0xlYF00KRBK8DVrupsxxqSCteBTIJTgFy+GW2+FFSvSHZExpjKwBJ8CTz0FF1zgEvvDD8P336c7ImNMZWBdNCkwcKC7nzzZ3VvBMWNMKlgLPmCqMGUKrFplo2iMMallCT5gO3bAUUfBCy+UnGC1BG+MSQVL8AHbudPdh4+Dty4aY0wqWB98wMITfIsW7nlVe9eNMSlQbgteRB4VkbapCCYbhSf4KlWgWjUQSW9MxpjKwU8XzTfAMBH5UkSuFpF6QQeVTXbscPfVq7tkf+218N576Y3JGFM5lJvgVfUZVT0KuBgoBOaIyEsi0jvo4LJBfj689hr06uVa8E8/DdOnpzsqY0xl4Oskq4jkAK292zpgNnCjiLwSYGxZoXZtOPNMKCx0fe9Vq9ooGmNMavjpg/8brpumL3C/qnZV1YdU9XSgc9ABZrpffoEPPoB169zzmjUtwRtjUsNPC34e0FFVr1LVaaXWdQ8gpqzyzTdw8skwY4Z7nptrwySNManhJ8EPVNWt4QtEZDyAqm4MJKosEj6KBqBu3fTFYoypXKKOyBaRXCAPaCQi+wChwX11gQNSEFtWKJ3gv/sufbEYYyqXWJfcXAXcgEvmX4Ut3wQ8FWBMWaV0gjfGmFSJmuBV9QngCRG5TlX/nsKYskrpBH/vvZCTA7fdlr6YjDGVQ6wumuNUdQLwg4icWXq9qr4eaGRZ4sgj3SiaAw90z8ePd/eW4I0xQYvVRdMTmACcHmGdApbgfdh3XzjxxJLnubmwYUP64jHGVB6xumju9O4vS1042WfJEpg1C/r2dWPgc3NtHLwxJjX8XOj0RxGpK84zIvKViJxY3uuM8+GHcPbZsNEbUGoJ3hiTKn7Gwf9OVTcBJwL7ApcBDwYaVRYpKnL3oZOs+fk2Ft4Ykxp+KpOHxr/3BZ5X1dkiVvDWr9KjaP5u45GMMSnipwU/U0Q+xCX4D0SkDrA72LCyh42DN8aki58EfzkwGDjMK1lQHddNE5OI5IrINBGZLSLzReTuvYw1I4USfLVq7v6ll6B///TFY4ypPMrtolHV3SLyE9BGROKZbG4HcJyqbhGRasBkEXlPVb9INNhMdMUVbphkqFPru+/gzTdh1y53wZMxxgSl3IQtIg8B5wELgF3eYgUmxXqdqiqwxXtazbtpwpFmqKZN3S0kNPH29u1Qq1Z6YjLGVA5+WuQDgENUdUe8O/cmCpkJHAQ8papfRthmEDAIoHnz5vEeosKbMgVWrIDzz3fPa9Z095bgjTFB89MH/19c6ztuqrpLVTsBTYHuItIuwjbDVLWbqnbLz89P5DAV2ogR8Kc/lTwPteCtJrwxJmh+WvBbgVleDfjfWvGqer3fg6jqLyIyETgZN4FIpbFz554jaBo2hFatYLeNQzLGBMxPgn/Tu8VFRPKBIi+51wSOBx6Kdz+ZbufOkhE04OZnPbNM6TZjjEk+P6NoRngJurmqLopj342BEV4/fBVgjKq+nWCcGat0C94YY1LFTy2a04FZwPve804iUm6LXlXnqGpnVe2gqu1U9Z69jjYDlU7wX38NvXq5AmTGGBMkP100d+Em154IoKqzRKRFgDFllX/8Y8/iYr/+Cp9+CmvXpi8mY0zl4CfBF6vqxlLlZyrdePZENWu25/PwYZLGGBMkP8Mk54nI/wA5ItJKRP4OTAk4rqwxejS88UbJcxsmaYxJFT8J/jqgLW6I5Mu4SbdvCDCmrDJ0KPzrXyXPw69kNcaYIPkZRbMVGOLdTJx27NjzJGudOtC5s9WEN8YEL2YLXkQu8WZw+tW7zRCRi1MVXDYoPYpm333hq69gwIC0hWSMqSSituC9RH4DcCPwFW7ijy7AIyKCqr6QkggznI2DN8akS6wW/DXAGar6iapuVNVfVHUCcJa3zviwcyfUqFHyXBW6dduzX94YY4IQqw++rqouK71QVZeJiPUg+zR16p6lCkRg7lxYvjx9MRljKodYCT7WQD4b5OfTAQeUXZaba6NojDHBi5XgDxWRORGWC9AyoHiyzoMPwmGHQZ8+Jctq1rQEb4wJXswEn7Iosthdd8ENN+yZ4HNz7UInY0zwoiZ4VbVe4r2kGnkUzdFHu5rwxhgTpHgm0TZx2rXLJfnwUTQAI0emJx5jTOXip1SBSdDOne7exsEbY9LBTz3400TEPggSEC3BX3wxnHde6uMxxlQufrpozgeeEJHXgOdVdWHAMWWNevXgl1/KJvg1a2DDhrSEZIypRMptmavqhUBnYAnwvIhMFZFBIlIn8OgynIhL8qEa8CE2TNIYkwq+ul5UdRPwGvAKbq7VM4CvROS6AGPLeOvWwS23lJ2ez4ZJGmNSwU8ffD8RGQdMAKoB3VX1FKAjcFPA8WW0tWvhkUdgUampyu1KVmNMKvjpgz8bGKqqk8IXqupWEfldMGFlh9BJ1vBaNOCKjRljTND8JPjVpZO7iDykqreq6viA4soK0UbRXHtt6mMxxlQ+fvrgT4iw7JRkB5KNbBy8MSadoiZ4Efm9iMwFWovInLDbUiBSETJTSrQE//jjkJ8PxcUpD8kYU4nE6qJ5CXgPeAAYHLZ8s6r+HGhUWaJ3b5fERfZcXlTkRths3w61a6cnNmNM9ovVRaPehB/XApvDbohIg+BDyw45OVCl1Lucm+vubSSNMSZIsRL8S979TGCGdz8z7Lkpx8yZ8Pvfw6pVey63BG+MSYWoCV5VT/PuW6hqS+8+dLMJP3xYvNjNvbp5857LQ1e2WoI3xgQpah+8iHSJ9UJV/Sr54WSXaCdZDzoILrywbAkDY4xJplgnWR+LsU6B45IcS9aJluB79HA3Y4wJUqwZnXqnMpBsZOPgXR2edu2gqk0tY0zKxRoHf5x3f2akW+pCzGy5uWVLFXzxheue+eij9MSUKmPGQOfO8PXX6Y7EmMopVruqJ67A2OkR1inweqwdi0gz4AVgf2A3MExVn0gwzox07bWRyxJUq+ZOsGZzRckff4RrroGuXaFTp3RHY0zlFKuL5k7v/rIE910M/FlVv/Jqx88UkY9UdUGC+8sa2T6KRhUGDYItW+DFF2HrVnctQB2bQcCYlPJTLrihiDwpIl+JyEwReUJEGpb3OlVdHRppo6qbgYVAk70POXOMGQOXXFJ2eWgcfLa24IcPh7feggcfhLp1Yd993TJjTGr5KTb2CrAWOAtXOngtMDqeg4hIIW5WqC8jrBskIjNEZMbatWvj2W2FN3MmjI7wTmX7hU7vvQe9esH110OTJtCqFbz6arqjMqby8TO2oYGq3hv2/K8iMsDvAUSkNm42qBu8maH2oKrDgGEA3bp1U7/7zQQ7d0YeQVO3Llx9NRx6aOpjSrZVq+Dnn+HXX93tmGPch9qmTSUlGs49F+66y217wAFpDdeYSsVPC/4TETlfRKp4t3OBd/zsXESq4ZL7KFWNeVI2G+3cCTVqlF1euzb8859w7LGpjymZnn/etdDbt3fj+vv0cUXUQnPRhpxzjuuXf+219MVqTGUU60rWzbjRMgLcCIz0VlUBtgB3xtqxiAjwLLBQVf+WlGgzTLQWPMCuXbB7d9khlJli61Y332yPHnDjjVCrlrs1iFCG7tBD3Vj4V1+F62wWX2NSJtYomr0d83AUcBEwV0RmectuU9V393K/GaNWLdh//+jrbrjBnYjMRHl5rkXesCG0bVv+9k895U62GmNSx9f1hSKyD9AKyA0tKz2NX2mqOhnX+q+0Hn88+rpMnnh71y5XBjmeLqZM744yJhP5GSZ5BTAJ+AC427u/K9iwsl9ubmYOk1SFU0+FO+6I/7UTJsDddyc/JmNMZH5Osv4ROAxY7tWn6YwbKmnKcffd8Oc/R16XqS34cePggw8S626ZPNm9J6Xr4xtjguEnwW9X1e0AIlJDVb8BDgk2rOwwZYq7RVKzZuYl+G3b3AdWu3ZumGe8bDSNManlJ8GvFJH6wBvARyLyH8DaYD7EGkUzaBD067d3+091S/i552DZMnjyycSqQx56qBtS+corSQ/NGBNBuQleVc9Q1V9U9S7gdtzQxwEBx5UVYiX4P/0JBg5MfN/vvefGoL//fuL7iNfEidCihZtMPFEXXOC+1SxfnrSwjDFR+GnBIyJdROR6oAOwUlV3BhtWdoiV4H/9FTZsSHzfM2e6+/HjE99HvO69F555Zu/2cf757kNi2bKkhGSMiaHcL9oicgdwDiXlgZ8XkVdV9a+BRpYFDjgAGjeOvO7cc+Gnn2BGgtOXh7pIUjmRRuvW7rY3WrSAJUvc1a7GmGD5SQ8XAJ3DTrQ+CHwFWIIvx3/+E31dgwbwzTeJ77t1a9fd8cADie8jHjNnwrx5cN55JcXSEiUCRUXuJLOVEDYmOH66aJYRdoETUANYEkg0lUiDBrB+feKvHzAAXnopaeGUa9QoN3Kmiq9Ovdi2boVmzeDRR/d+X8aY6GJN2fd3EXkS2AHMF5HhIvI8MA9Xi8aU4+yz4bEoU5c3aAAbN0JxcWL7LipyJ1hPPjk1F0x98QV065ac+WXz8lx5g5dfdsMmjTHBiNVFE+odngmMC1s+MbBossxnn0F+fuR1oaJcv/wCjRrFv+82bWD1aney9ttvoUOHhMMs144drovm+uuTt88LLoArr4SvvnLT+hljki9qC15VR4RuwMu4RD8TeMlbZsoRaxTN0UfDww8n1iJWdWPge/RwzxcuTDxGP2bNcj/LEUckb59nnukqab78cvL2aYzZk59aNL2Ab4GngKeBxSJipaN8iJXgO3eGm292k3/Ea/Nm14/ds6c7Ybk3J2v9mD3b3Yc+UJKhQQM46SQ3Ocju3cnbrzGmhJ9RNI8BJ6rqIgARORjXorcv1uWIleCLimDFCtc9Ez45hh+rV7v7li2hoCD4BB+66jZa6eNE3X57ZhZcMyZT+BkTUS2U3AFUdTGQodNUpI6q6xdv2jTy+mXL4KCD3OTU8QqVKGjc2E2Rl8i3gHglO7kDdO/uvoUkY2SOMaYsP/9aM0XkWRHp5d3+jeuLNzGIuBOTv/995PWhk6w//xz/vhs3dkW/Dj4YXngB/u//Eo+zPD/+6EYDzQzoN75wIQwenPhoImNMdH4S/NXAfOB6XOngBd4ysxfq13cfAomMhW/d2o0hj/btIJmmTnXVH4uKgtn/woXw0EPwySfB7N+YyixmgheRKsBMVf2bqp7pFR4bqqo7UhRfxvrlFzf879VXI6/PyXFJPpEW/Pr1bngkwIIFrnzvxx8nGmlsU6e68widOwez/1NOcVezjh4dzP6NqcxiJnhV3Q3MFpHmKYona2zb5sZ4x0rgDRokluCvuw46dnSPGzWC+fPdLQhffOGSe40awey/Zk3o3x9ef92dlDbGJI+fLprGuCtZx4vIm6Fb0IFlulCyijXO/b774PLL49/3qlUlRczy82GffYIZSVNU5IqhJXP8eyTnnecqa370UbDHMaay8TNM0mbRTICfBH/eeYnte/Vq6NTJPRZxffJBJPiffnIlBY4+Ovn7DnfiiVBYaFP5GZNsURO8iOTiTqYeBMwFnlVVG+vgk58E//33bpTKYYfFt+/Vq13fdUjr1m4CkGRr2hSmT0/+fkurXt2VELbhksYkV6wW/AigCPgMOAVogxtFY3yoUcON8d5vv+jbPPSQqwgZTz/8li3uStbwOvPHHeeGGRYXJ7c+vGrq6rZXqeKOt3Ur1KqVmmMak+1ipYM2qtoeQESeBaalJqTscNBBboq7WBo2dKNtdu1yo2r8EIGhQ90FTiEXXuhuydaunZuY5M47k7/vSHr2dJOk2JytxiRHrC/Fv418tq6ZYDRo4FqtGzf6f02tWnDDDWUrMKomdxTKypVuCGb9+snbZ3natnVX9oaGgBpj9k6sBN9RRDZ5t81Ah9BjEdmUqgAz1ZQp7krTWFeAJnI169q1sGjRnld+7t7tSgnccUdisUby+efu/qijkrfP8px3nuuiefvt1B3TmGwWq1xwjqrW9W51VLVq2OMUVD/JbBs3ujrtsS7BDyX4eK5mHTXKnVTdFPYRW6WKGw+fzJE0n3/uJuYIjbdPhWOOcecWRlgxamOSwsYtBMTPKJpu3WDsWNdf79eqVe4E7j777Lk82UMlJ0+Gww93NdtTJScH/vAHNyIoqNo3xlQmSRxzYcL5SfD77QdnnRXfflevdq3c0qNbWreGN990FyftbVJWdQXGmjXbu/0k4vrroXnzknH+xpjEWYIPiJ8EX1zsimwVFkKrVv72G34Va7jWrd3+lixxj/eGCNx2297tI1G1awczIsiYysi6aAKy//5w6qnl12o/8UQ3Ft6v1avdUMLSDj8c/vKX5Iwh/+9/4xvZE4Rnn3XT+tmk3MYkzlrwAenTx91iqVrVzeYUzyiahx921RdLO/hguP/++GKM5uqr3Widr79Ozv4SsWMHjBvn6tOceGL64jAmkwXWgheR50RkjYjMC+oY2SDeipKnneYuCIpk0yY3DeDeKC52JYKPPHLv9rO3Lr/cTUf4v/9rrXhjEhVkF81w4OQA91+hPf2060rZsiX2dg0a+B8muWkTjB/vKi9G0rcvXHRRfHGWNneuiznoAmPlqVHDjeufPj2xaQ2NMQF20ajqJBEpDGr/Fd3Gja6/vLwRLfG04OfOheOPd8MIT47w0dmqFXz4Yfyxhps82d2n8gKnaC6+GB54AG69FU4/PXV1cSoqVVizxl1hnJMDxx7rll9xhbvyeNMmNw/Btm3u/M9jj7n1oW98tWq5k/7Vq8NJJ5WUqn7sMfd3mJ/vRnaFbsmcA2DU3FEMGT+EFRtX0Lxec+7rcx8D2w9M3gHKsWuXaxitW+fKg/To4ZaPH+/mbdixo+RWtWpJd+fgwfDOO3t+i9x//5IJdgYOLPmfy8lxt1atSsqUXH45zJnj3svcXHdr08Z1tQLce68bOHHHHZEHT+wt64MPSGgUTXkJ/sEH/Seu1avdfaSTrOD+sIYPd5f6J3qy9fPPXRXJ5hVgipeqVd10gUuWuPeouBgmTIATTqhcyf6vf3UJZe7cksZAnz4lSWbJEncFcN267oK3mjWhZcuS19ev7xoca9e6v8uiIpdkwL3uppvKHvOWW1wxvK1bXZI6+GD399W0qUv+zZq5EU9+jJo7ikFvDWJr0VYAlm9czqC3BgHEneR373bfMGvVcsl0xQqYPdv9bOvWldyeesq9D/ff7z7ANmzYM0kXFbm/r7Fj4V//cstEXCKuV68kwTds6K5TCa8VFX4NyhFHuO3BfYiErioPyc93tx073Afvhg17Drz49FP3AXDjjXG9Db6JBtjB6bXg31bVdjG2GQQMAmjevHnX5cuXBxZPKt12GzzySHLnMv3739048TVr3B9NaWPGuMv9Z81K/ArU+fNda/Ckk/Yq1ECMHOm6oDp2dP9YnTu78fKdOsUejpoJfvgBvvzSTbAyY4b7kA6VizjzTPfh3rGjS8xt2ri6Pclo8am6Y61d625r1rgS1u3bu5FZS5e6D5Ply13yCvnHP+Daa92HzqmnuoRbq5ZLXvXru4R19NHudZ1GFPKLlv2/LqhXwPjTl/HOOy65VqniPsRXr3azljVp4kaY3XKL+6DZtg22b3ev/fZbl3gffRRuvrlkn7m57kNu+nSXaMeNcx+EjRqV3PbZxzUScnJcZdbQ66pWzcyGg4jMVNVukdalvQWvqsOAYQDdunXLmtNpO3f6SzrffAPTprmx3+XVQ1+1yv0RNmwYeX1oLP233yae4Nu2dbeK6Lzz3D/kK6+4f/xQy+v7713L8oEH4J//dP/Effq4bp0jj0xuCeUgDB7sWsvgYm3fHrp3dwm1ShX3LSaoxCPiWuK1a0OLFmXXt2jhhs3u2OHuV692HwKhYne1arn3+tdfS0pZL15cUjBu0SL4ZfcKiBD/io0rmD0b/liqCHm1au5316SJ+72efLJrjeflufs6dUqK4J1/PvTq5X7n+flum/D36owz3C2aSCPSsknaW/DhunXrpjNmzAgsnlQaOdJdxPTss7G3GzrUtXY2bCi/cuOll7ouimgjZX791dVxOfnkPb+i+zV9uvtwOPvsit8iVoVly9zX23793D/1a6+5QmUrVsBnn7lvTwcc4FqRVavCpEkuYbZsGflq4FTYvt2dNH7hBfeNrLDQ9eHOmOESZceOrjWZLXbvhsInCvl+U+QW/OJrlrFli9suNP9AgwY2+Us8YrXgUdVAbsDLwGpc2eGVwOXlvaZr165a2QwfrgqqS5aUv+2CBaoTJgQXyzXXqNaurVpUFNwxUmXjRtWxY1X/9reSZd26ufcaVAsKVK+8UvWdd4KPZfdu1SlTVK+6SrV+fXf8Jk1Ux48P/tgVwcg5IzXvvjzlLn675d2XpyPnjEx3aFkBmKFRcmqQo2guCGrf2SS8ZHB5re5DD3W3WJYudX3o4ROC+KHqWr89e1b8Lg0/6tYtW+fn5Zfhu+9cF8LEiTB6tOvb7dvXrR8xwvXNRjuJHa8dO9xJu3Xr3O+jenXXn37JJW4WLr+TvGS60InUdI6iqbSiZf503LKpBf+736l26FD+dpMnuxbdBx+Uv+3Ikarz58fe5vLLVffbz1+M4b74wsUxfHj8r81URUWqa9e6x4sWuZ9fRPXoo1Xvv1/1s89Ut2/f8zUj54zUgqEFKneJFgwtKNMK/fFH1aefVj3mGNVevUqWv/+++1ZhTLKRjhZ8Zbdli78ZlvxO+rFjhzsRe889JUPcImnVCn76yY2JLq8OTrgxY1wLs39//6/JdFWrupNz4IYBzp/vhs2NHVtSbO2111yre+5cuP+tUbxePIidGjbc70033K/O0oE8/rgb9rZ7tyv4dvHFJf3KFXFUksl+diojIH5H0bRs6Wqfn3JK7O1+/NHdlzc0LnwkTTzmznVJKJVT9FU0bdq4C07mzHEjRcaNcyM0wM3Q9craIb8l95CtxVsZMn4IS5e6UU5DhrjXL1jgir9l4rA7kz2sBR8Qvwm+Rg3o0qX87cq7yCkkPMGXnrc1lg8+sLlQw+Xnw4ABJc+vugp+f/cKIo05W7FxBdcOcdcoWEI3FYm14APiN8GDO7kXuioxmlCCL68Ff+CB7j6eFnyoG8HvlYmVVfN6kS/vbV6vecZeJGOymyX4gPTtu2cLMJY774QXX4y9zZIl7r6gIPZ2eXnw/vtw2WX+jr17N7Rr58Zkm9ju63MfedXy9liWVy2P+/rcl6aIjInNumgC8qc/+d/WT8Gxq65yQ+1CJ2VjieeE3pQprr84dLLRRGfD/UymsQQfkOJiN87Zz9d2Pwm+Th1XG8SP+fNdVcirrip/2zFj3JWTp53mb9+V3cD2Ay2hm4xhXTQB6dTJXfLvR3kJvrgYbr/dVc3z4/333axM0erGh+za5YYE9u2b/TU5jKmMLMEHJJ6TrOUl+MWLXcnYOXP87c/vUMnPP3cnb889199+jTGZxRJ8QHbu9D9hwr33xm6dz53r7jt08Le/UIJfvDj2dgcc4AqdnXqqv/0aYzKL9cEHpKjIfws+Um33cHPmuP781q397a9lS1eNr7wW/EEHlcz6Y4zJPtaCD0g8XTQLFsBdd7miVJHMneuSu99vBDVquBmZYiX4hx5yE0wYY7KXJfiA3Hyzm2zCj2+/hbvvjl7nfelSNwlEPD7+GIYNi7xu2jQ3ycS4cfHt0xiTWayLJiC33OJ/2/IKjs2ZE38ZgdAVraXt3u2mQ2vc2NVNMcZkL2vBB2DtWli/3v/2oSn4or0mkTIC8+bBrbeWfGiMmjuKwscLybm3CtOOLKTfkFE2NNKYLGcJPgBPPOFmnt+xw9/2sVrwr7/uyg5s3Vp2XSzLl8PDD7s5X0Oz2i/fuBxQqL+cFzcOYtTcUfHt1BiTUSzBB2D27PhOiu6zj7uPlOA//tj1ldesGV8MoaGSjzwC174+hK1FpcrcFrkyt8aY7GV98AGYMweOOsr/9jVquKtOI03QMWeOO8Eab6XCFi3cBM7vvgs7O0af1d4Yk72sBZ9kGza40TAdO8b3uvr1y84kr+qGSPq9wClctWowa5brJopV5tYYk70swSdZqJxAvAn+73+Hp57ac9mKFW7qvXiHSJZ2//FW5taYysgSfJIdcgg88wx07x7f6954A155Zc9l69e7vvx4PyxKG9h+IMNOH0ZBvQIEoaBeAcNOH2ZVEY3JcuIm5a4YunXrpjNmzEh3GGlxzjnuitb589MdiTEmk4jITFXtFmmdteCT7J13SmZfikeotMCUKcmPyRhTOVmCT6LiYjjrLHj66fhfe9ttbjq+AQNKShb06eNq1BhjTCIswSfR4sVu1EoifeYNG8Lbb7uumk/WjaJgaCETjqnCUArtgiRjTEJsHHwSJTqCJuSQQ+DIq91Vp1uLtoLAJpYz6K1BAHZS1BgTF2vBJ9Hs2W78+aGHJr6PIePtqlNjTHJYgk+i2bNdcvdbBz6SaFeX2lWnxph4WRdNEj37rKskuTea12vuFQUru9wYY+JhLfgkatw4sbIC4e7rY1edGmOSwxJ8ksyb58rzxlMHPhK76tQYkyx2JWuSDB0KN94Ia9aUP4m2McYkS9quZBWRk0VkkYh8JyKDgzxWus2eDfvvb8ndGFNxBJbgRSQHeAo4BWgDXCAibYI6XrrNmrX3RcGMMSaZghxF0x34TlX/CyAirwD9gQXJPtDcuXD++WWXP/II9O0LX34Jv/td2fVPPw09e8KECW4i6tJGjIBu3dwVprfeWnb92LFuWOTo0S6Gk07a+5/FGGOSJcgE3wT4Puz5SuDw0huJyCBgEEDz5okNBaxZE9pE+G5Qr567z8uLvD40kXWdOpHX5+WV7CfS+txcd7/PPnDuuXDJJfHHbowxQQnsJKuInAOcpKpXeM8vArqraoS2spPJJ1mNMSYd0nWSdSXQLOx5U2BVgMczxhgTJsgEPx1oJSItRKQ6cD7wZoDHM8YYEyawPnhVLRaRPwAfADnAc6pq8xUZY0yKBFqLRlXfBd4N8hjGGGMis1IFxhiTpSzBG2NMlrIEb4wxWcoSvDHGZKkKVU1SRNYCZWe7iK0RsC6AcJLJYkwOizE5MiFGyIw4K0KMBaoascxhhUrwiRCRGdGu4qooLMbksBiTIxNihMyIs6LHaF00xhiTpSzBG2NMlsqGBD8s3QH4YDEmh8WYHJkQI2RGnBU6xozvgzfGGBNZNrTgjTHGRGAJ3hhjslTGJviKPqG3iDQTkU9EZKGIzBeRP6Y7pmhEJEdEvhaRt9MdSzQiUl9ExorIN957ekS6YypNRP7k/a7nicjLIpJbAWJ6TkTWiMi8sGUNROQjEfnWu9+nAsb4iPe7niMi40SkfhpDjBhj2LqbRERFpFE6YoslIxN8hkzoXQz8WVUPBXoA11bAGEP+CCxMdxDleAJ4X1VbAx2pYPGKSBPgeqCbqrbDlciOMFNwyg0HTi61bDAwXlVbAeO95+k0nLIxfgS0U9UOwGLgL6kOqpThlI0REWkGnACsSHVAfmRkgidsQm9V3QmEJvSuMFR1tap+5T3ejEtITdIbVVki0hQ4FXgm3bFEIyJ1gWOBZwFUdaeq/pLWoCKrCtQUkapAHhVgBjNVnQT8XGpxf2CE93gEMCCVMZUWKUZV/VBVi72nX+BmhEubKO8jwFDgFqBCjlbJ1AQfaULvCpc8Q0SkEOgMfJnmUCJ5HPcHujvNccTSElgLPO91JT0jIrXSHVQ4Vf0BeBTXklsNbFTVD9MbVVT7qepqcA0RYN80x1Oe3wHvpTuI0kSkH/CDqs5OdyzRZGqClwjLKuQnqIjUBl4DblDVTemOJ5yInAasUdWZ6Y6lHFWBLsA/VbUz8Cvp71bYg9eP3R9oARwA1BKRC9MbVeYTkSG47s5R6Y4lnIjkAUOAO9IdSyyZmuAzYkJvEamGS+6jVPX1dMcTwVFAPxFZhuvmOk5ERqY3pIhWAitVNfQNaCwu4VckxwNLVXWtqhYBrwNHpjmmaH4SkcYA3v2aNMcTkYhcApwGDNSKd8HOgbgP89ne/09T4CsR2T+tUZWSqQm+wk/oLSKC6zNeqKp/S3c8kajqX1S1qaoW4t7DCapa4Vqdqvoj8L2IHOIt6gMsSGNIkawAeohInve770MFOxEc5k3gEu/xJcB/0hhLRCJyMnAr0E9Vt6Y7ntJUda6q7quqhd7/z0qgi/e3WmFkZIL3Tr6EJvReCIypgBN6HwVchGsVz/JufdMdVAa7DhglInOATsD96Q1nT963i7HAV8Bc3P9W2i9jF5GXganAISKyUkQuBx4EThCRb3EjQB6sgDH+A6gDfOT97/yrAsZY4VmpAmOMyVIZ2YI3xhhTPkvwxhiTpSzBG2NMlrIEb4wxWcoSvDHGZClL8CZlRGSXN+Rtnoi8FW+FQBGZKCLdvMfvlvd6EVkWqcKft3yud1sgIn8VkRreugNEZGyMfdYXkWviiTuZxJng1ecJLTvDq2bYOmxZvoi8n54oTUVhCd6k0jZV7eRVW/wZuDbRHalq370sONZbVdvjCte1xBuzrqqrVPXsGK+rD6QtwQN9gdmlyl5cAEwmrHqlqq4FVovIUSmOz1QgluBNukzFKxAnIt1FZIpXSGxK6IpVEakpIq94NcFHAzVDLw5vnYvIGyIy06vFPiieIFR1C3A1MMCrk14YqvktIm1FZJr3rWOOiLTCXRR0oLfsERGpLSLjReQr7xtBf++1heLq1v/bi+tDEanprTtIRD4Wkdne6w70lt8sItO9Y90dJeSBhF156tU6Ogq4nLLlid/wtjeVlarazW4puQFbvPsc4FXgZO95XaCq9/h44DXv8Y3Ac97jDriiU92858uARt7jBt59TWAe0LD0NqXiKLMcmAUcDhQC87xlf8fVQQGo7u3/t/Xe8qpAXe9xI+A7XDG8Qi/eTt66McCF3uMvgTO8x7m40sIn4r5FCK7h9TZwbITYlwN1wp5fCDzrPZ6Cu1w+tK4JMDfdv3e7pe9W1efngDHJUFNEZuGS30zcpA4A9YARXgtZgWre8mOBJwFUdY5XpiCS60XkDO9xM6AVsD7O2CJVKJ0KDBFXM/91Vf3WlZkp87r7ReRYXMnlJsB+3rqlqjrLezwTKBSROkATVR3n/VzbAUTkRFyS/9rbvrb3c0wqdbwG6uYXCLkAV/IZXMG4C3DlEsAVETsg9o9tspl10ZhU2qaqnYACXIs41Ad/L/CJur7503Gt2pCYtTREpBeu1X+EqnbEJci4psrzkm4hbuagkgOrvgT0A7YBH4jIcRFePhDIB7p6P9tPYcffEbbdLlxrP9IHCd7yB9Sdo+ikqgep6rMRtisWkSpe3A2B44BnvIqGNwPnScmnUK4Xu6mkLMGblFPVjbjp7W4SV1K5HvCDt/rSsE0n4fUhi0g7XDdNafWADaq61RtF0iOeWLw+7KeBN1R1Q6l1LYH/quqTuAqMHYDNuCJY4cdfo6pFItIb9+EVlbqToytFZIB3jBriaot/APzOiwcRaSIikSbiWIQ7KQxwNvCCqhaoq2rYDFgKHO2tPxjXZWUqKUvwJi1U9WtgNu7E4MPAAyLyOa5/PuSfQG2va+YWYFqEXb0PVPW2uRc3vZsfn3gnU6fhSv1eFWGb84B5XrdSa1wyXQ987g31fAQ3EUU3EZmB+zD6xsexL8J1K83B9Zvvr272p5eAqSIyF1eZsk6E174D9PIeXwCMK7X+NeB/vMe9ve1NJWXVJI3JIOIm6HhBVU/wse0koH/pbyam8rAWvDEZRN0cqv8Ov9ApEhHJB/5myb1ysxa8McZkKWvBG2NMlrIEb4wxWcoSvDHGZClL8MYYk6UswRtjTJb6fz6dzMPKLQLsAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(easy_rdf_bins, easy_rdf_data, \"b--\", label=\"rdf\")\n", "plt.plot(easy_rdf_bins[peaks], easy_rdf_data[peaks], \"go\", label=\"critical points\")\n", "plt.xlabel(\"Radial Distance (A)\")\n", "plt.ylabel(\"Probability Density\")\n", "plt.title(\"RDF minima using scipy.find_peaks\")\n", "plt.legend()\n", "\n", "print(f\"solvation radius is {easy_rdf_bins[peaks[0]]}\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks great! Our result is very close to our previous spline interpolation. Now lets try integrating it into our `Solute`.\n", "\n", "First we need to **wrap our function** so that it's signature is `foo(bins, data) -> cutoff radii`. Remember that our function must return the radii at the **first minima** as a float and that find_peaks requires the **negated** data to find minima rather than maxima." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def find_peaks_wrapper(bins, data, **kwargs):\n", " peaks, _ = find_peaks(-data, **kwargs)\n", " radii = bins[peaks[0]]\n", " return radii" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice we passed our `**kwargs` to the the find peaks function. This lets us use the `kernel_kwargs` argument to `Solute` to pass arbitrary kwargs to our method." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": "" }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solute_scipy_peak = Solute.from_atoms(\n", " li_atoms,\n", " {'PF6': PF6, 'BN': BN, 'FEC': FEC},\n", " rdf_kernel=find_peaks_wrapper,\n", " kernel_kwargs={\"distance\": 5}\n", ")\n", "solute_scipy_peak.run()" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'PF6': 2.75, 'BN': 2.55, 'FEC': 2.8500000000000005}\n" ] }, { "data": { "text/plain": "(
,\n )" }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": "
", "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGi0lEQVR4nO29eZwU1dX//z4Mu+wyiiwCKogo++oGqNGIG65BghqNBnkSTfzl0WhivonZjCHGGKPRGOPyJCIaIy4RFRUQlR1ENgEBMQyowAADyDpwfn/cKqemp7qneuh1+rxfr35Vd9261ad6oD51zrn3XFFVDMMwDCOWOtk2wDAMw8hNTCAMwzCMUEwgDMMwjFBMIAzDMIxQTCAMwzCMUEwgDMMwjFBMIAzDMIxQTCAKGBEZJiIl2bYjHYiIishxaTjv6SKyItXnjfNda0Xka977n4jIY5n43nQgIpeIyDoR2SkifbJtjxENEwgjEqkUExHpLSLzRWSXt+2divOmg1ihUdV3VfX4TNuhqner6g3VHSci00Sk2uOywL3ATaraRFU/ONSTede5xxOczSLygogc5bU9KSL7vDb/NTLQ90oR+UhEvhSR1SJy+qHaU1sxgTAyiojUB14C/gm0BJ4CXvL2G7WXjsDSmnQUkaI4TTepahOgK9AC+GOgbZwnRv7rWe9cZwO/A64DmgJDgDU1sasQMIGoJYjI7SKyXkR2iMgKETnL299ARO4XkQ3e634RaRDnHJWelr0nsV+LyGHAa0DbwBNZWxGpIyJ3eE9hpSLynIi0qsbUYUBd4H5V3auqDwACnBnHpvNEZJl3XetF5NZA23dEZJWIbBGRl0WkbUj/wSLyefAm44U7FnnvB4rITBHZJiKficiDvliJyHSvy4f+U2isJyUiJ3hPs9tEZKmIXBTz+z0kIq969s8WkWPj/TAicrWIfOr9lnfGtN0lIv/03jcUkX96x20TkbkicqSI/AY4HXjQs/dB7/g/eeGd7Z7HdnrMeZ8Tkf/zbFwqIv0D7R28p/NN3vc9GGj7tvckvlVE3hCRjiHX1EBEdgJF3u+4OuLv9rCITBKRL4Ez4v1mAKq6Bfg3cFKi4zx+AfxSVWep6kFVXa+q6yP0K0hMIGoBInI8cBMwQFWbAl8H1nrNdwKDgd5AL2Ag8NNkzq+qXwLDgQ2BJ7INwPeBi4GhQFtgK/BQNac7EViklYuALfL2h/F34Ebvuk4CpgCIyJnAb4FvAEcBnwITQmyfBXxJZQH6JjDee38A+P+A1sDJwFnAd72+Q7xjegWfQn1EpB7wCjAZOAK4GXja+3v4jMLdlFoCq4DfhF2kiHQHHgauxv2WhwPt4/wm3wKaAx2848YCu1X1TuBdKkI5N3nHz8X9/Vt51/0vEWkYON9FuN+uBfAy4AtLEfAf3G/bCWjnHYeIXAz8BLgUKPa+95lYQ72HgCbex16qemzE3+2b3m/VFHgvzu+AZ0tr4DIgYejKu57+QLH3YFHiPRA0StSvkDGBqB0cABoA3UWknqquVdXVXtto3BPTRlXdhLtZXZ2i770RuFNVS1R1L3AXcLmI1E3QpwlQFrOvDHcjCGM/7rqaqepWVV3g7R8NPK6qC7zv/jFwsoh0CjnHM7gbNSLSFDjP24eqzveeJstVdS3wV5zgRWGwdz33qOo+VZ2Cu6GOChzzgqrOUdVy4GncjTqMy4H/qOp073r+H3AwzrH7ccJwnKoe8K5hezwjVfWfqlrqXeMfcP9Wgjfj91R1kqoeAP6Be5AA9zDRFrhNVb9U1T2q6t+sbwR+q6ofedd2N9A7zIsIIcrv9pKqvu895e+Jc54HRGQb8CHwGfDDQNutnneyTUQ2e/uOBOrhfuvTcX+LPiT5wFRImEDUAlR1FXAL7ga9UUQmBMItbXFPgD6fevtSQUdgov8fEfgIJ1ZHJuizE2gWs68ZsCPO8Zfhbuifisg7InKyt7/SdanqTqAU95Qby3jgUnGhtUuBBar6KYCIdBWR/3hhqO24G13rBPYHaQusU9XgjfzTGBs+D7zfhbsxxj1X4Hq+9K4njH8AbwATxIUNx3lP5aGIyP96oaAy7+/UnMrXGGtjQ0/kOwCfegIQS0fgT4G//RZcqDDs948lyu+2jur5vqq2UNV2qjraewDyuddra6Gq/rXu9rZ/VtXPVHUzcB/u35cRgglELUFVx6vqabj/uIpLxAFs8Pb5HO3tC2MX0DjwuU3wK0KOXwcMD/xHbKGqDauJ6S4FeoqIBPb1JE4CU1XnquoIXCjiReA5r6nSdYnLkxwOVPluVV2GuwENp3J4CVxYZznQRVWb4cImEnuOOGwAOohI8P/R0WE2ROAz3A0ZABFpjLueKqjqflX9hap2B04BLgCu8ZuDx3r5httxobiWqtoC57FFucZ1wNFxPMJ1uNBf8G/fSFVnRDhvlN8t5esQqOpWoCQd566tmEDUAkTkeBE503tC3oN7UjrgNT8D/FREir1Y7c9wI4jCWAh8U0SKRORcKodavgAOF5HmgX2PAL/xwwred4yoxtxpnm3f9xKYfpx8Ssh11ReR0SLSXFX3A9sD1zUeuE7ckNkGuCf/2V6YKIzxuJzJEOBfgf1NvfPuFJFuwP/E9PsCOCbOOWfj8hs/EpF6IjIMuJCQXEgEngcuEJHTxCXJf0mc/58icoaI9PBi6ttxISf/d4m1tylQDmwC6orIz6jqwcVjDk647hGRw8Qlx0/12h4BfiwiJ3o2NReRKyKeN5W/W7I8AdwsIkeISEuc5/2fDHxvXmICUTtoANwDbMaFC47APQkD/BqYh0sELwYWePvC+AHuP+o2XIz/Rb9BVZfjxGaNF1ZoC/wJl9ScLCI7gFnAoESGquo+XGL7Gu97vg1c7O0P42pgrRf+GQtc5Z3nbVyc/t+4m9ixwJUJvvoZ3AiqKV5owedWnFexA/gb8GxMv7uAp7xr/kbItVyE80w2A38BrvF+q6RQ1aXA93BC9hku4R9v3kkbnKBsx4X13qFC9P+EywNtFZEHcKGo14CVOC9qD9HCN3g5iQuB44D/evaM9Nom4rzUCd7fZgnud4hy3pT9bjXgV7ik/Urcb/cBcQYOGCBqK8oZhmEYIZgHYRiGYYRiAmGkHC9vsDPkVaOZtIZhZAcLMRmGYRihJJrQlHe0bt1aO3XqlG0zjEyywiusenzG6+cZRq1g/vz5m1W1OKytVglEp06dmDdvXrbNMDLJsGFuO21aNq0wjLxFRD6N12Y5CMMwDCMUEwjDMAwjFBMIwzAMI5RalYMIY//+/ZSUlLBnT7yCkEY+0LBhQ9q3b0+9enFr0hmGkWJqvUCUlJTQtGlTOnXqROX6cEa+oKqUlpZSUlJC586ds22OYRQMtT7EtGfPHg4//HAThzxGRDj88MPNCzSMDFPrBQIwcagF2N/QMDJPQQiEYRhGbUQVXn8dfve76o+tCSYQWeSuu+7i3nvvBWD58uX07t2bPn36sHr16mp6GoZRyBw4AM8+C337wvDh8MgjsHt39f2SxQQig6gqBw+GLzP84osvMmLECD744AOOPfbYDFtmGEa+8OabrrLMlVc6UXj8cVdxplGj1H9XrR/FlG3Wrl3L8OHDOeOMM5g5cyYXX3wx48ePp0OHDhQXF9OvXz8mTZrE/fffT1FREdOnT2fq1KnZNtswjBxEFb77XTh4EJ5/Hi6+GIqK0vd9BSUQt9wCCxem9py9e8P99yc+ZsWKFTzxxBNcf/31XHvttXzwwQeUl5fTt29f+vXrx3nnncfYsWNp0qQJt956a2oNNAyj1rB4MaxaBX/9K1x2Wfq/z0JMGaBjx44MHjyYd999l0suuYTGjRvTrFkzLrroomybZhhGHvGvf0GdOs5zyAQF5UFU96SfLg477LCv3ttwTcMwfFThL3+BwYOhX7/qj/3Xv2DoUDjiiMzYZx5EBhkyZAgTJ05k9+7d7Nixg1deeSXbJhmGkUXWr4ebboJTT4Xx4xMfu2yZS0ZffnlmbIMC8yCyTd++fRk5ciS9e/emY8eOnH766dk2yTCMLLJypdsedRSMHu1yDL/5jQsjxfL88yACl16aQQNVNW0v4FxgBbAKuCOkfQSwCFgIzANOi9o37NWvXz+NZdmyZVX2GflJ6N9y6FD3Mow85C9/UQXVNWtUb7zRvb/wQtWysqrHnnSS6pAhqbcBmKdx7qlpCzGJSBHwEDAc6A6MEpHuMYe9DfRS1d7At4HHkuhrGIaR16xYAY0bQ6dObrLbQw/BpEkuCV1eXnHc8uWwZElmw0uQ3hzEQGCVqq5R1X3ABJzH8BWqutNTMIDDAI3a1zAMI99ZuRK6dnWhI3BzHB5/HKZOhdtvrzju+efdNqPhJdIrEO2AdYHPJd6+SojIJSKyHHgV50VE7uv1HyMi80Rk3qZNm1JiuGEYRk356pE3AitWuFnRQa65Bm6+Ge67ryJx/fzzcMop0C70Lpg+0ikQYeM5q/x0qjpRVbsBFwO/Sqav1/9RVe2vqv2Li4traqthGMYh88Yb0Lo1bNxY/bF798LatVUFAuAPf4DTT4cbbnDi8OGHcMUVKTe3WtIpECVAh8Dn9sCGeAer6nTgWBFpnWxfwzCMXOC992DLFnj77eqPXb3alczo2rVqW716bs5Dq1bwjW+4fZkOL0F6BWIu0EVEOotIfeBK4OXgASJynHgzx0SkL1AfKI3S1zAMI9dYtcptp0yp/lh/iGuYBwFw5JHwwgtOLAYNgqOPTo2NyZA2gVDVcuAm4A3gI+A5VV0qImNFZKx32GXAEhFZiBu1NNIbeRXaN1225hLDhg1j3rx5Nep79913V/p8yimnpMKkGtOpUyc2b96cE7YYRiZIRiBWrHDbMA/CZ+BAmDmz+kl06SKtE+VUdRIwKWbfI4H3vwNCl7oI62sk5u677+YnP/nJV59nzJiR8u8oLy+nbt3k/9mkwxbDyCVU4eOP3bDVNWtcfqFTp/jHr1wJbdpAs2aJz9u3byqtTA4rtZFmvvzyS84//3x69erFSSedxLPPPgvA22+/TZ8+fejRowff/va32bt3b6V+Dz/8MD/60Y+++vzkk09y8803A3DxxRfTr18/TjzxRB599FEA7rjjDnbv3k3v3r0ZPXo0AE2aNAHcZMjbbruNk046iR49enxlw7Rp0xg2bBiXX3453bp1Y/To0WjIEIxhw4bxk5/8hKFDh/KnP/2JV155hUGDBtGnTx++9rWv8cUXXwBQWlrKOeecQ58+fbjxxhsrncu3Zdq0aVxwwQVf7b/pppt48sknv7qG7t2707NnT6tqa+QdW7ZAWRmMHOk+V1e1f8WKxN5DLlBYpTayUO/79ddfp23btrz66qsAlJWVsWfPHq699lrefvttunbtyjXXXMPDDz/MLbfc8lW/yy+/nJNPPplx48YB8Oyzz3LnnXcC8Pjjj9OqVSt2797NgAEDuOyyy7jnnnt48MEHWRhyfS+88AILFy7kww8/ZPPmzQwYMIAhQ4YA8MEHH7B06VLatm3Lqaeeyvvvv89pp51W5Rzbtm3jnXfeAWDr1q3MmjULEeGxxx5j3Lhx/OEPf+AXv/gFp512Gj/72c949dVXvxKvKGzZsoWJEyeyfPlyRIRt27ZF7msYuYAfXhoxAl591YWZrrsu/vErVsAll2TGtppiHkSa6dGjB2+99Ra333477777Ls2bN2fFihV07tyZrt7jw7e+9S2mT59eqV9xcTHHHHMMs2bNorS0lBUrVnDqqacC8MADD9CrVy8GDx7MunXr+PjjjxPa8N577zFq1CiKioo48sgjGTp0KHPnzgVg4MCBtG/fnjp16tC7d2/Wrl0beo6R/mMRUFJSwte//nV69OjB73//e5Yudemh6dOnc9VVVwFw/vnn07Jly8i/U7NmzWjYsCE33HADL7zwAo0bN47c1zByAV8gunaFM890AhFvTsSWLbB5s3kQuUUW6n137dqV+fPnM2nSJH784x9zzjnnRF4HYuTIkTz33HN069aNSy65BBFh2rRpvPXWW8ycOZPGjRszbNgw9uzZk/A8YWEjnwYNGnz1vqioiPLg/P4AwZLlN998Mz/84Q+56KKLmDZtGnfddddXbdWVM69bt26lZVd92+vWrcucOXN4++23mTBhAg8++CBTomT6DCNHWLXKzYju3NkJxIQJLs8QNkqpuhFMuYJ5EGlmw4YNNG7cmKuuuopbb72VBQsW0K1bN9auXcsq75HjH//4B0OHDq3S99JLL+XFF1/kmWee+eoJvqysjJYtW9K4cWOWL1/OrFmzvjq+Xr167N+/v8p5hgwZwrPPPsuBAwfYtGkT06dPZ+DAgTW+prKyMtp5UzqfeuqpSt/z9NNPA/Daa6+xdevWKn07duzIsmXL2Lt3L2VlZbztDRjfuXMnZWVlnHfeedx///2hoTLDyGVWrYIOHaBhQycQEH80ky8Q5kEUOIsXL+a2226jTp061KtXj4cffpiGDRvyxBNPcMUVV1BeXs6AAQMYO3Zslb4tW7ake/fuLFu27Ksb+rnnnssjjzxCz549Of744xk8ePBXx48ZM4aePXvSt2/fr27UAJdccgkzZ86kV69eiAjjxo2jTZs2LF++vEbXdNddd3HFFVfQrl07Bg8ezCeffALAz3/+c0aNGkXfvn0ZOnQoR4cM3O7QoQPf+MY36NmzJ126dKFPnz4A7NixgxEjRrBnzx5UlT/+8Y81ss0wssWqVXDcce79Mce4eQtTpsD//E/VY1escGtJH3NMZm1MFkkUfsg3+vfvr7FzCD766CNOOOGELFlkpJLQv+WwYW47bVqmzTGMShQXu9nOf/2r+3zddfDKK67sRuz6Dldc4cpn+J5ENhGR+araP6zNQkyGYRiHyLZtLunsexDgwkylpW4RoFjCivTlIiYQhmEYh8jq1W4bFIgzznDb2DzEwYNuQp0JRI5Qm8JohYr9DY1cxh/iGhSI9u1dEjpWINatgz17cj9BDQUgEA0bNqS0tNRuMHmMqlJaWkrDhg2zbYphhOJPRYpNOp95JrzzTuXV4fwaTPngQdT6UUzt27enpKQEW0wov2nYsCHt27fPthmGEcqqVdC2LQSmCwFw1lluKdGnn4Zvfcvty5chrlAAAlGvXj06d+6cbTMMw6jFBIe4BrnoIhg6FMaMce2nnuo8iKZNXaG+XKfWh5iM3OHzz8GrE2gYtYp4AlG/Pvz739CxI1x8savyGrsOdS5jAmFkjN/+Fq68EnbuzLYlhpE6duyAL74IFwiAww93xfsOHIDzz3fDXvMh/wAmEEYG8eeybdmSVTMMI6WEDXGNpUsXmDjRHfvZZ/mRfwATCCNDlJbCokXuvQmEUZvwh7h26ZL4uKFDwa+A37t3Wk1KGbU+SW3kBt5SEoAJhFG78AXi2GOrP/baa90Eug4d0mpSyjCBMDJCsFSSCYRRm1i1Co480o1MikLHjum1J5VYiMnICFOnQvfu7r0JhFGbiDeCqTZgAmGknU2bYMkSV+kSTCCM2oUJhGEcAv5qqued5xZTMYEwagu7dsH69SYQNUJEzhWRFSKySkTuCGkfLSKLvNcMEekVaFsrIotFZKGIzIvta+QPU6dC48bQvz+0bAkhC80ZRl6yZo3b1laBSFuSWkSKgIeAs4ESYK6IvKyqywKHfQIMVdWtIjIceBQYFGg/Q1U3p8tGIzNMmwannQb16kGrVuZBGPnDnj2uPHfjxuHtYVVcaxPp9CAGAqtUdY2q7gMmACOCB6jqDFX1nydnAVaNrZaxcSMsXVpRG98EwsgnRo2CgQNh9+7w9g8+cNsoQ1zzkXQKRDtgXeBzibcvHtcDrwU+KzBZROaLyJh4nURkjIjME5F5VrE19/DnP/grg5pAGPnCvn3wxhvuAef226u2L14M48bB8OEudFobSadAhJWiCl2UQUTOwAlE8M9wqqr2BYYD3xORIWF9VfVRVe2vqv2Li4sP1WYjxUyb5kog9+vnPptAGPnC7NnOc+jdG/78Z5g8uaJt1y4YORKaN4cnnsiaiWknnQJRAgTnC7YHNsQeJCI9gceAEapa6u9X1Q3ediMwEReyMvKMqVPh9NNd/gFMIIz8YcoUV3F10iQ3h+faa13JGIBbboGPPoJ//MNNkqutpFMg5gJdRKSziNQHrgReDh4gIkcDLwBXq+rKwP7DRKSp/x44B1iSRluNNPDFF+4/kR9eAicQu3a55J9h5DJTpkDfvnDUUfDPf8LmzTB2rCtZ/7e/wR13wNlnZ9vK9JI2gVDVcuAm4A3gI+A5VV0qImNFZKx32M+Aw4G/xAxnPRJ4T0Q+BOYAr6rq6+my1UgPfnkNP0ENTiDAhroauc2uXTBzplsyFKBPH/jlL+H55+Gaa2DwYPe5tpPWWkyqOgmYFLPvkcD7G4AbQvqtAXrF7jfyi5kz3fDAvn0r9vkCsWWLezIzjFzk/fdh//4KgQC47TZ47TX48EMYP74ibFqbsWJ9Rtr4/HNo1w7qBv6VBQXCMHKVKVPcv9vTTqvYV1TkEtVlZXDEEdmzLZNYqQ0jbZSWutW0gvjDAS3EZOQyU6bAoEHQpEnl/Q0aFI44gAmEkUbCBMI8CCPXKSuDefMqh5cKFRMII22YQBj5yPTprryGCYQJhJFGSksrBMGnWTMXyzWBMHKVKVNc1eHBg7NtSfYxgTDSwt698OWXVT0IEZeHMIEwcpUpU+DUU51IFDrVCoSI3CsiJ2bCGKP24AtArECAzaY2cpfNm2HRIgsv+UTxIJYDj4rIbG+SW/N0G2XkP35JAhMII5/wJ3eaQDiqFQhVfUxVTwWuAToBi0RkvFdgzzBCMYEw8pEpU6BpU7e4lRExB+Et/tPNe20GPgR+KCIT0mibkceYQBj5yLRprrhkXZtCDETLQdyHCzOdB9ytqv1U9XeqeiHQJ90GGvmJCYSRa+zYAb16ueqsYezcCcuXuwlyhiOKB7EE6KWqN6rqnJg2K8FthOILROwwV3CjmMrK4MCBzNpkFDYvvugS0C++GN6+cCGoVqxdYkQTiNGquiu4Q0TeBlDVsrRYZeQ9paWuLEHYWr6+aGzbllGTjAJn/Hi3nRP7mOuxYIHbBotLFjpxI20i0hBoDLQWkZZUrBDXDGibAduMPGbLFhdekpB1BYOzqcNCUIaRajZuhDffdLWVlixx5bxjH17mz4c2bazKcJBEHsSNwHxcYnqB934+8BLwUPpNM/KZsDIbPlZuw8g0zz3nQpo/+YnbfvBB1WMWLLDwUixxBUJV/6SqnYFbVbVz4NVLVR/MoI1GHmICYeQSTz8NPXu6ZUMB5s6t3L5rFyxbZuGlWBKFmM5U1SnAehG5NLZdVV9Iq2VGXlNa6tbxDcMEwsgka9bArFnwu9+58FH79lXzEIsWuQJ95kFUJtFo36HAFODCkDbFrSVtGKGYB2HkCn5y+sor3XbAgKoehCWow4krEKr6c297XebMMWoDqokT0C1auK0JhJFuVF14acgQOPpot2/gQJg40f378x9W5s+H1q2dd2FUEGWi3A9EpJk4HhORBSJyTiaMM/KT7duhvDx8DgS4WarNm5tAGOln4UI3+e2b36zYN2CA286bV7HPT1CHjborZKLMg/i2qm4HzgGOAK4D7kmrVUZek6iSq4/NpjYywfjx7oHk8ssr9vl1lvw8xJ49buirhZeqEkUgfE09D3hCVT8M7DOMKiQqs+Fja0IY6ebAAXjmGRg+vPK/xebNoVu3ijzEkiXO47UEdVWiCMR8EZmME4g3RKQpcDDKyUXkXBFZISKrROSOkPbRIrLIe80QkV5R+xq5SxSBaNUKtm7NjD1GYfLBB7B+PYwcWbVtwADnQai6/AOYBxFGFIG4HrgDGOCV3KiPCzMlxKsA+xAwHOgOjBKR2IGPnwBDVbUn8Cvg0ST6GjlKVIEwD8JIJ/6N/5RTqrYNHAiff+4EZMEC59F26pRR8/KCaovaqupBEfkC6C4iyRTBHQisUtU1AF5p8BHAssC5ZwSOnwW0j9rXyF1MIIxcINGN309Uz5njhKRvX0tQhxFlFNPvgPeBnwK3ea9bI5y7HbAu8LnE2xeP64HXku0rImNEZJ6IzNu0aVMEs4x04wtEy5bxj/EFQjUzNuUy48fDpVWmolawbh2MGuXKURvRSXTj79UL6tWD99+HxYstvBSPKCGmi4HjVfU8Vb3Qe10UoV+YHofeDrzV6a4Hbk+2r6o+qqr9VbV/cXFxBLOMdFNa6hKBiRZdadXKJRF37MicXbnK88+7cfl794a3T54MEybAjBnh7UZV9u1LfONv2NCV3nj6aXesJajDiSIQa4B6NTh3CdAh8Lk9sCH2IBHpCTwGjFDV0mT6GrlJlCqtNpu6gsWL3XZDnH/h69e77TILsEZm2bLqb/wDB8IXX7j35kGEE0UgdgELReSvIvKA/4rQby7QRUQ6i0h94Erg5eABInI0rmTH1aq6Mpm+Ru6SqMyGjwmEY9cuWL3avTeBSB1RRib5eYhmzeDYY9NvUz4SJen8MjW4OatquYjcBLwBFAGPq+pSERnrtT8C/Aw4HPiLuEBhuRcuCu2brA1GdigtdWULEmEC4Vi2rCIP4wtBLCYQybNgATRtmvjGP9BbD7NPH6gT5VG5AIkyiukpEWkEHK2qK5I5uapOAibF7Hsk8P4G4IaofY38oLQUjj8+8TEmEA4/vATRBELVRttEwU9QJ7rxd+sGRx4Jp5+eObvyjSijmC4EFgKve597i4iFe4y4RAkx+SOcCl0glixxCdOGDeMLxIYNLuG/dWtFzNyIT3k5fPhh9XmFoiJYuhT+3//LjF35SBTH6i7cvIRtAKq6EOicNouMvGb/flesL6pAFMJs6nXr4ucXFi9262a0axcuEPv2ueUyTz7ZfbYwU/UsX+7qK0VJPB9+ONSvn36b8pUoAlGuqmUx+2z0uhGK7xHEq+Tq06iRexWCB3H11W4eQxhLlkCPHk4gwkTks8/c9uyz3XapZeKqxU9Q29DVQydKknqJiHwTKBKRLsD3ARuRbYQSpZKrT6HMpl63Dv77XzfRrUmTiv2lpU4ATjrJeQqxq5xBhVfRv79bR8M8iOpZsAAaN4auXbNtSf4TxYO4GTgR2As8A2wHbkmjTUYeE6XMhk9tEIiDB+G44+DJJ+Mfs3mzi4u/917l/UuWuG2PHtC2rROD2JnlvkC0a+dCUSYQ1bNgAfTu7XIMxqFRrUCo6i5VvVNVB3hDUO9U1T2ZMM7IPwpNILZudfMYPvggvH3fPpeTAZgypXKbP4LppJOcAOzZUzUnYwKRHAcOuL+FhZdSQ0KBEJFveSvIfem95onINZkyzsg/Ck0gNm9223hlwILXFysQS5a4ZH3btk4AoGqiesMGaNDA/Vbdu7vvs5Jj8fn4Y/jyS5sZnSriCoQnBLcA/wu0xRXL+xHwAxMJIx6FKhD+Nl579+4u9BH0EBYvduElkQqBiE1Ur1/v2kTcOcC8iEQsWOC2JhCpIZEH8V3gElWdqqplqrpNVacAl3lthlGF0lJXJTOYjI1HbRAI/2k+3lO9LxBXXOHyC9Onu8+qzoM46ST3uW1bt431INavr2g78US3NYGIz/z5bk5Jd1s9JiUkEohmqro2dqe3r1m6DDLym9JSd+OPMtu3VSsXd9+9O/12pYvqQkx++/nnu2G9fphp3TqXm+jRw31OJBC+d9GunSsfYQIRnwULXJXWRJWEjegkEohE/23z+L+0kU6iVHL1qQ2zqYMCEba2hd/erh2cdlqFQAQT1ODyDK1bVxYI1coC4YeZTCDCOXjQCYSFl1JHIoE4IbBedPC1GOiWKQON/CJKmQ0f/7h48ft8wPcc9u0LX9AnmJM580wXVtq4sWKIqy8QUHWy3LZtzrtqF1gqq3t3mywXjzVrnFdmI5hSRyJH7ISMWWHUGkpL3byAKPjrO+XzqJyguG3a5EJAse1NmzoP4Ywz3L5p05wH0aGDm/zm48+F8PHFIlYgnngiOSEuFJ5/3m0HD86uHbWJuB6Eqn6a6JVJI438IZkbV20QiKDtYdexeXPF79GvnxOLKVMqJ6h9YusxBedA+PjJ148+OnTbaxPbtsG4cXDeeVV/V6PmWBV0I2WoFp5AbN5ckUsJC5Vt3lyxNkbdujB0KLz5prvB+wlqn3btXPhp/373OUwgbCRTOH/8oxtC/KtfZduS2oUJhJEyvvzSxeKjCkSrVq5ef74LxAleMDaeBxFcPOnMM12sfN++cA9CtaJAny8QRx1VcUyHDnDYYSYQQTZvhvvug8suswR1qomyHsQFImJCYlSLn5CtrpKrT1GRE5N8FohNmxILROzqen4eAqp6EP5QVz/3sH69+30aNqw4pk4d930mEBX87nfu4eSXv8y2JbWPKDf+K4GPRWSciFji2ohLMrOofYqL81cg9uxxI5c6d3ZrClQXYgI3Rr9VKyeO3WLGAsaW2wgOcQ1iQ10r+OwzePBBuOoqmxyXDqIU67sK6AOsBp4QkZkiMkZEmlbT1Sgwkin17ZPPAuELYnFx+HXs3Qs7dlT+PerUcZPm+vat7BlAcgKxfj2Uxa7SUoD85jeuUu5dd2XbktpJpNCRqm4H/g1MAI4CLgEWiMjNabTNyDMKzYPw7W7dOvw6/N8j6EEA/PWv8NZbVc93+OGuTIkvEBs2xBcIqJhLUah8+ik8+ihcfz0cc0y2ramdRMlBXCQiE4EpQD1goKoOB3oBt6bZPiOPKDSB8ENKrVu7V+x1BNuDNGoEzUKK1dSpUzEXYv9+t/50mECcfLILaU2YcOjXkM8884z7ne68M9uW1F6ieBCXA39U1Z6q+ntV3QhunQjg22m1zsgrkk1SgxOILVtcHf98wxcEP8QUm4OI50Ekwp9N/fnnbkRTmEC0bg0jR7pFivy1JgqR2bPdqnEdOmTbktpLFIH4TFWnB3eIyO8AVPXtRB1F5FwRWSEiq0TkjpD2bl5OY6+I3BrTtlZEFovIQhGZF8FOI8uUlrqJYMksAl9cXDF/It8IeghhnlA8DyIR/mQ5P8zkj2yK5aabXIL8H/9IzuZMsmhR+gRMFWbNgkGD0nN+wxFFIM4O2Te8uk4iUgQ85B3bHRglIrHjDLbg1ri+N85pzlDV3qraP4KdRpapSfmHfJ4st3mzK6DXqpUTge3bXWI62A7JCYQfYgqbJBdk4EAYMMCN4AkrEphtZs+GXr2gfXu45Ra36l4qWbfOeVkmEOkl0YJB/+MX5osp1vcJsCjCuQcCq1R1jaruwyW4RwQPUNWNqjoX2H8I12DkCOvWwZFHJtcnnwVi06aKIav+dQQ9IV8gkgm5tWvnxvT7pTTiCQQ4L2L58qor1WUbVbj1VvdvYcQIeOgh6NIFLr44dYUGZ892WxOI9JLIgxgPXAi85G39Vz9v6Gt1tAPWBT6XePuiosBkEZkvImPiHeQNuZ0nIvM25eNdppawZ49z+U85Jbl++SwQwTkOYdexebNLRicTcvMFYe5c1y+R9/GNb7j2Bx9Mzu5089JL8N578ItfuBDYp5+6RPK777pCei+/fOjfMXu2K4DYs+ehn8uITyKBUG9xoO8BOwIvRCTKM1HYkjHJOMOnqmpfXIjqeyIyJI6Rj6pqf1XtX+z/LzUyzsyZLrwSnCkchXwWiE2bKuz3b+TB64idRR0FXyDmzHHhpkQLLzVsCN/5jrvhfpoj5TP374fbb3ezva+/3u1r29bVSFq0yE0OvPhiuOeeQwuNzZrlih8mI75G8lTnQQDMB+Z52/mBz9VRAgTHF7QHNsQ5tgqqusHbbgQm4kJWRo4ydaobpjkkVMbj4+csNm5MvU3pJooHkaxA+Enpzz9PHF7yGTvWbR9+OLnvSRePPQYrV7ryF7GrurVr55ZcHTkSfvxjuPpq53kmy/79bmlRCy+ln0Tlvi/wtp1V9Rhv67+iTEuZC3QRkc4iUh9XsiOScykih/kztUXkMOAcoMCnBeU2U6e62cHNmyfXr149F6PPRw8iTCCCQ11rIhBBUYgiEEcf7eL8jz2W/aVbd+xwM5qHDIELLgg/plEjGD8efv1rePppuO225L9n8WInLCYQ6SfugkEikrAuoqouqKa9XERuAt4AioDHVXWpiIz12h8RkTY4b6QZcFBEbsGNeGoNTBTnX9cFxqvq65Gvysgou3a5mPAtt9Ssfz5OllN1AuALg78Od6wH4Zfnjkrjxm4RoW3b4g9xjeWmm2DiRHdzvueeaOuBp4Nx45wn+J//JLZBxOUkZs8On1FeHZagzhyJVpT7Q4I2Bc6s7uSqOgmYFLPvkcD7z3Ghp1i242ZqG3nA++87tz/Z/INPPgpEWZmrAeR7CEVFVT2h4GJBydCunROIKB4EuN/9hhvcDfrgQbfNtEh88QX84Q9w5ZVu+G0UBg2CV15xv2UynuesWXDEEdCxY81sNaITVyBUtYb/3Y1CY+pUd4M87bSa9S8udnHrfMIPJQXHRQRnU+/Z44arJhtiAuc5LF0aXSBEXH2nBg3g3nudR/fnP7ucUKbwQ1y/+EX0PgO9rOLcufC1r0XvN3u2Gw2VLU+pkEgUYjpTVaeIyKVh7ar6QvrMMvKJqVPdU2PsesxRKS52Xkg+ESzU5xP0hGpSZsPHF4aoAgFODP78Zxfjv/deJ1CPPuqEO90cOOC+62tfc6UvouJ7GnPmRBeIrVthxQq45prk7TSSJ1GIaSiuQN+FIW0KmEAY7NjhngBvv73m5ygudjfUgwcz+9R7KITNkm7d2t284rVHpSYCAe6Jetw4l8f45S/dgkQ1zQslw+uvw3//61Z1S4YWLeD44ytyClGYM8dtLf+QGRKFmH7uba/LnDlGvvHee+4Jsqb5B3ACceCAezqsScw+GwQL9fkEPaFDEYizznKTympShE7EhXnGj3d/m0wIxMMPQ5s2cNFFyfcdNAjeeMMl/aOEjGbPdsdFzXMYh0aUct+Hi8gDIrLAm9X8JxHJk//GRrqZOtUNVU12BnWQfJwsFyYAQU/Ib6+J4J1xBrzzzqFNAuvVy01MSzeffgqTJrkkeb16yfcfONAluNetq/5YcAJxwgnh5dKN1BPFoZ8AbAIuw5X+3gQ8m06jjPxh6lSXMGzcuObnyFeBaNAADjusYl/r1hWe0KHkIFJBz56wapVLlKeTv/3NPdF/5zs16++HivzQUSJUKxLURmaIIhCtVPVXqvqJ9/o10CLNdhl5QFkZLFhwaOElyE+B8MtsBMMiweuoSaG+VNKzp7uhpqo4Xhj797vRS+ed5ybs1YSePZ2nFCUPsXq1E17LP2SORElqn6kiciXwnPf5cuDV9Jlk5AvTp7twSiEKRNgs6eBs6s2bXRK2JmGXVOAXsVu0qGI4aap56SUXHvLLfdSE+vWhT59wD+LgQVi2zC0otW2b81bBBCKTJBrmugM3WkmAHwL/9JrqADuBn6fdOiOnmTrVhVkO1eUPK3SX6wRnUfsEr6MmZTZSSadO0KRJevMQDz/sPIdzzz208wwa5DyR8vLK9Zuuuw7+7/8qH9umTfKz042ak6gWU1NVbeZt66hqXe9VR1UtRWTw/vtOHBo2PLTzNGjgko75JBCbNsX3IDZtqlkl11RSp44b5ppqgfj8c1c99sc/dutQjBlz6HMtBg50k/uWLavYt2yZKxV+3XXw5pswb57LqaxeXbUIoJE+Iv3UItIS6AJ8dSuIXYbUKDxSGQ8+4oj8EogwD8H/7IeYjjoq83YF6dkTnnsu+hDSeKjC44+7kt1+WfGiIjj9dLjxxkO30/83NHt2RWjsN79xAx/Gjcuu0BY6UYa53gBMxxXd+4W3vSu9Zhn5wPbtNZ89HUtxcf6U/N63zyXoY0NMjRq5UU25EGICd7PdurVi+dKa8OWXcO21bhhr+/au3tJ777m//fTpqbnGY491yXw/D/HxxzBhAnz3u9n/DQudKB7ED4ABwCxVPUNEuuGEwihwduxI3Xj04mJYuzY150o3iYaw+uU2ckUgwIWZ2oeVxKyGjz6CK65w4Z677oKf/jQ9pTtEXJjJH8l0990u7Pi//5v67zKSI8ow1z2qugdARBqo6nLg+PSaZeQ6+/e7ej+p9CDyJcQUVqjPp7jYlZ3YtSv7AtGjh9vWJA8xbZqbrbxxo5vp/POfp7eu08CBbkju4sUu93Djjcmvb26knigeRImItABeBN4Uka0ksTKcUTvZscNtU+lBbN586PHyTBBWqM+nuNglVOO1Z5LmzV1J7JoIxH33uf5z5iRfE6omDBzohrWOHu2S0DVZSMhIPdUKhKpe4r29S0SmAs0BW7ynwNm+3W1T6UHs3+9i+y1apOac6SJRnaXWrStyKblQV6pnz+QFYu9eN0LpmmsyIw5QMVdj8WKXe4i6WJKRXiLVzhSRviLyfaAnUKKq+9JrlpHr+B5EKgUC8iPMFFaozye4L9seBDiBWL7c3fSj8t57Ljk9fHj67IqluBg6d3YTCw+lMrCRWqKMYvoZ8BRwOG4p0CdE5KfpNszIbXwPIpUhJsgPgUhURiMXBeLAAZdwjsprr7kZzmdWu2ZkarnjDreWRU3LdhipJ0oOYhTQJ5CovgdYAPw6nYYZuU0hexCJymjErg+RbYIjmXr3jtbntddgyJDKhQgzwZgxmf0+o3qihJjWEpggBzQAVqfFGiNvSEeSGvJDIPxCfWH4+0WgZcvM2RSP445zM92j5iH++183rPVQy2cYtYNEtZj+jKvFtBdYKiJvep/PBt7LjHlGrpKOJDXkh0AkmuPgX0fLlrlREqJuXVe7KFYgFi92bSecUHn/697wk0zmH4zcJdE/YW+wHvOBiYH909JmjZE3pNqDCM5CznU2b3bDR8PwhSMXRjD59OzpFvXxmTIFzj/fhclWrqws8q+95nIAscJhFCaJivU95b+AZ3BCMR8Y7+2rFhE5V0RWiMgqEbkjpL2biMwUkb0icmsyfY3skmoPAvJnslxYoT4f34PIhfyDT8+eriz3F1+4EUoXXuiGkX7+uZu17LNvH7z1lvMecn0uipEZooxiGgZ8DDwE/AVYKSJDIvQr8voMB7oDo0Ske8xhW4DvA/fWoK+RRXbscE/9qQyj5EPBPtXEIabmzd1vkmsCARWL+3ToADNmwNVXuwlxn3zi2t9/H3butPCSUUGUJPUfgHNUdaiqDgG+DvwxQr+BwCpVXePNm5gAjAgeoKobVXUusD/ZvkZ2SWWhPp988CB27HBP2vGS1CIuRNOhQ2btSoRfcuOnP3Ui/PbbrozFb39bedby66+7kVmZHt5q5C5Rnv/qqeoK/4OqrhSRKOtktQOCS5GXAFGLQx9KXyMDpLJQn09xcXoXuEkFfrnrRKW833wzt2aDFxc70RJx+Qd/dnS7dm7uwc9+Bu+84/IPp52WeuE38pcoHsR8Efm7iAzzXn/D5SKqIyyKqRHtitxXRMaIyDwRmbcp1x8/axE7dqTPg9Co/0qywKxZbptoGc9jjsneWtTxeOMNV1cpdhLarbc6b+c733Ejmyy8ZASJIhBjgaW4XMEPgGXevuooAYKOdnuiF/mL3FdVH1XV/qravzie32+knHSFmPbscXHwXGXmTJdfOO64bFuSHN26ufBSLI0auUV5Pv7YfTaBMIIkDDGJSB1gvqqeBNyX5LnnAl1EpDOwHrgS+GYG+hoZYMeO1BdyC86FyNUwx4wZbpnV2jTKZ+RIeOght7CQrfdsBEnoQajqQeBDEUm6OoqqlgM34Vag+wh4TlWXishYERkLICJtRKQE+CHwUxEpEZFm8foma4ORPtLlQUDuJqpLS2HFCjjllGxbklpE4D//gXffrV3CZxw6UZLUR+FmUs8BvvR3qupF1XVU1UnApJh9jwTef44LH0Xqa+QO6UhS+6ue/fe/8de6nj3bhUuaN0/td0fBzz+cfHLmvzvdNG+end/UyG2iCIQtL2pUIR0ehB/XX7EivH37djfK5q674M47U/vdUZg5062qNmBA5r/bMLJBolpMDXHJ6OOAxcDfvdCPUeD4y42m2oM47DA3omblyvD2FSugvBzWrEnt90Zlxgzo1SvzVU4NI1skykE8BfTHicNw3IQ5w0h5qe8gXbvG9yD8NQ1KSlL/vdVRXu6GidbG8JJhxCNRiKm7qvYAEJG/A3MyY5KR66S6UF+Q44+H8ePD16Zevtxt162r2i/dLFniVlmrbQlqw0hEIg/iq/IXFloygqTbg9i2LXwkUzY9iBkz3NY8CKOQSORB9BIRr2YnAjTyPgugqpqG50cjH0hHJVef449325Urq07s8j2IHTugrCyzo25mzoQ2baBTp8x9p2Fkm0TlvotUtZn3aqqqdQPvTRwKmHSHmKBqHmL/fli1Co491n3OtBcxY4bzHmyegFFIRCm1YRiVSKcHcfTR0KBB1ZFMq1e7RPHZZ7vPmcxDfPGFGzll+Qej0DCBMJImnR5EUZGbDxHrQfj5h2wIxMyZbmv5B6PQMIEwkiadHgS4RHWsB+HnH844A+rUyWyIaeZMt05Cv36Z+07DyAVMIIykSecoJnB5iFWrXEjJ56OPXHHAli1dsjiqB7FsWUWJjJoycyb07QsNGx7aeQwj3zCBMJJmxw53s6wXZdmoGtC1q0tK+4vzgPMgTjjBve/QIboHMWYMfOMbNbdl3z6YO9fCS0ZhYgJhJM327enJP/jEjmRSdQLRrZv73L59NA+irMx5D+vWuVLWNWHRIldWxATCKERMIIykScdqckG6dnVbPw+xYYP7zqAHsW5d9SvPTZ0KBw6497Nn18yWefPcNtEKcoZRWzGBMJImHZVcg7Ru7Zbs9D0IfwRT0IP48kvnISRi8mRXWK9+/ZoLxPz5Lu/RsWPN+htGPhOl3LdhVCIda0HEEhzJ5I9gCnoQ4LyIFgnOMXkynHmmm8dQ00T1/Plu9JJNkDMKEfMgjKRJtwcBLg8R9CCaNXOjl6BiYaFEierVq93rnHPcEqHz5lUeFeUzdSpceKFLiseyd68r0mfDW41CxQTCSJpMeRDr18POnRUjmPyn+KAHEY8333Tbs892q9Pt2uVu9rE8+KBbbnP+/KptS5Y44TCBMAoVEwgjaTLlQQB8/LHzIPz8A8BRR1U/WW7yZFe2o2tX50FA1TzEnj3wxhvu/bRpVc/hi4YJhFGomEAYSZMpDwLcHITPPqvIPwDUretEIp4HUV4OU6a48JIIdO7sEt+xeYi333bJ7nr14J13qp5n/nxo0cL1N4xCxATCSIrycti9O/0exHHHuZv7yy+7z0EPAiqGuoYxd64b4XTOOe6ziPMiYj2Il15y13HNNfDee1XzEPPnuxnUlqA2ChUTCCMp0lmoL0ijRi5E9NZb7nPQgwCXqI4XYpo82d3UzzqrYt+gQS5UtW2b+3zwILzyCgwfDl//ust1LFhQcfy+fbB4sYWXjMLGBMJIinTXYQpy/PFuJFG9enDMMZXbvposF9Jv8mQYMMDNpfDx8xBz57rtnDnw+ecwYgQMHer2BfMQS5c6kTCBMAqZtAqEiJwrIitEZJWI3BHSLiLygNe+SET6BtrWishiEVkoIvPSaacRnXRXcg3i5yG6dHF5hyDt27uRSbFDV7dtc6EkP7zkM2CA8yr8PMRLL7lzDh/uVq474YTKeQhLUBtGGgVCRIqAh4DhQHdglIh0jzlsONDFe40BHo5pP0NVe6tq/3TZaSRHpkJMUDGSKTb/ABVDXffurbzfL68RKxDNmzsR8PMQL70EQ4a4WdIAw4bBu+9WCM78+a6Pv4KdYRQi6fQgBgKrVHWNqu4DJgAjYo4ZAfyfOmYBLUTkqDTaZBwi2fAgYvMPUDFZLlYgJk+GJk0qQkpBBg1yHsTKlS4fMSLwr3HYsMp5iAULoE8fS1AbhU06BaIdEBxnUuLti3qMApNFZL6IjIn3JSIyRkTmici8TZs2pcBsIxGZ9CB69XLJ6rClPr/yIPZU7DtwwCWezzorvBT54MFQWgp//KP7HBSIYB5i/3748EMLLxlGOmsxhT17xeYUEx1zqqpuEJEjgDdFZLmqTq9ysOqjwKMA/fv3r6a+p3GoZNKDOPJId0Nv1KhqW5s2brJc0IN48003+/r++8PP53sVjz3mxCdYgO/IIyvyEF//ujuvCYRR6KTTgygBOgQ+twc2RD1GVf3tRmAiLmRlZJlMehAQLg7gEsxt21YWiCefdCOXLrwwvM+JJ7rqruXllb0Hn6FDXR5izhz32QTCKHTSKRBzgS4i0llE6gNXAi/HHPMycI03mmkwUKaqn4nIYSLSFEBEDgPOAUIq6RiZJpPDXKujfXvY4wnE1q3w4oswejQ0aBB+fFGRG80E4QIxbJi7vr/9zV3fccelw2rDyB/SFmJS1XIRuQl4AygCHlfVpSIy1mt/BJgEnAesAnYB13ndjwQmissQ1gXGq+rr6bLViM727eldbjQZOnSAvR+6988847yJ665L3Ofyy50H0adP1TY/DzF3rhvhVMdmCRkFTlrXg1DVSTgRCO57JPBege+F9FsD9EqnbUbNSPdqcsnQvr0TBQWeeAJ69oTevRP3+d733CuMNm3ckNrlyy28ZBhgM6mNJMlEJdeodOjgSmaUlbn1Hq677tCHpfpehAmEYZhAGEmSiUquUfGHun661iWtR48+9HNecIE7V9jQWsMoNGzJUSMpcsmDaN8e9gBbt8GFl0Bx8aGf84IL3BKlwTpOhlGomAdhJEUuehBQfXI6GUwcDMNhAmEkRS4lqdu0cTMt69dzRfcMw0gtFmIykmL79tzxIIqK3NN+8xZVq70ahnHo2H8rIylyyYMA6NEj2xYYRu3FQkxGZMrL3RoMueJBGIaRXkwgjMjs3Om2ueRBGIaRPkwgjMhkspKrYRjZxwTCiEymK7kahpFdTCCMyJgHYRiFhQmEERnzIAyjsDCBMCKTS2tBGIaRfkwgjMj4ISbzIAyjMDCBMCJjHoRhFBYmEEZkLEltGIWFCYQRmR073HrP9etn2xLDMDKBCYQRmVxaC8IwjPRjAmFEJpfWgjAMI/2YQBiRybVKroZhpBcTCCMyubQWhGEY6SetAiEi54rIChFZJSJ3hLSLiDzgtS8Skb5R+xqZxzwIwygs0iYQIlIEPAQMB7oDo0Ske8xhw4Eu3msM8HASfY0MYx6EYRQW6VxRbiCwSlXXAIjIBGAEsCxwzAjg/1RVgVki0kJEjgI6ReibMvr3h92703Hm2sUnn8CwYdm2wjCMTJFOgWgHrAt8LgEGRTimXcS+AIjIGJz3wdFHH10jQ7t1g717a9S1oDjxRLjmmmxbYRhGpkinQEjIPo14TJS+bqfqo8CjAP379w89pjr++c+a9DIMw6jdpFMgSoAOgc/tgQ0Rj6kfoa9hGIaRRtI5imku0EVEOotIfeBK4OWYY14GrvFGMw0GylT1s4h9DcMwjDSSNg9CVctF5CbgDaAIeFxVl4rIWK/9EWAScB6wCtgFXJeob7psNQzDMKqSzhATqjoJJwLBfY8E3ivwvah9DcMwjMxhM6kNwzCMUEwgDMMwjFBMIAzDMIxQTCAMwzCMUMTliWsHIrIJ+DTJbq2BzWkwJ5WYjanBbEwN+WAj5IeduWBjR1UtDmuoVQJRE0Rknqr2z7YdiTAbU4PZmBrywUbIDztz3UYLMRmGYRihmEAYhmEYoZhAeIX+chyzMTWYjakhH2yE/LAzp20s+ByEYRiGEY55EIZhGEYoJhCGYRhGKAUrECJyroisEJFVInJHtu0JQ0QeF5GNIrIk27bEQ0Q6iMhUEflIRJaKyA+ybVMsItJQROaIyIeejb/Itk3xEJEiEflARP6TbVvCEJG1IrJYRBaKyLxs2xOGt3Tx8yKy3Pt3eXK2bQoiIsd7v5//2i4it2TbrjAKMgchIkXASuBs3KJFc4FRqpqWNa9riogMAXbi1u0+Kdv2hOGtIX6Uqi4QkabAfODiXPotRUSAw1R1p4jUA94DfqCqs7JsWhVE5IdAf6CZql6QbXtiEZG1QH9VzfbkrriIyFPAu6r6mLeeTGNV3ZZls0Lx7kXrgUGqmuwk37RTqB7EQGCVqq5R1X3ABGBElm2qgqpOB7Zk245EqOpnqrrAe78D+Ai3pnjOoI6d3sd63ivnnoxEpD1wPvBYtm3JV0SkGTAE+DuAqu7LVXHwOAtYnYviAIUrEO2AdYHPJeTYTS0fEZFOQB9gdpZNqYIXulkIbATeVNWcsxG4H/gRcDDLdiRCgckiMl9ExmTbmBCOATYBT3ihusdE5LBsG5WAK4Fnsm1EPApVICRkX849UeYTItIE+Ddwi6puz7Y9sajqAVXtjVvffKCI5FTITkQuADaq6vxs21INp6pqX2A48D0vDJpL1AX6Ag+rah/gSyBXc4z1gYuAf2XblngUqkCUAB0Cn9sDG7JkS97jxfX/DTytqi9k255EeOGGacC52bWkCqcCF3kx/gnAmSLyz+yaVBVV3eBtNwITceHaXKIEKAl4iM/jBCMXGQ4sUNUvsm1IPApVIOYCXUSks6fiVwIvZ9mmvMRLAP8d+EhV78u2PWGISLGItPDeNwK+BizPqlExqOqPVbW9qnbC/XucoqpXZdmsSojIYd5ABLywzTlATo2wU9XPgXUicry36ywgZwZMxDCKHA4vQZrXpM5VVLVcRG4C3gCKgMdVdWmWzaqCiDwDDANai0gJ8HNV/Xt2rarCqcDVwGIvxg/wE29N8VzhKOApb8RIHeA5Vc3JYaQ5zpHARPdMQF1gvKq+nl2TQrkZeNp7+FsDXJdle6ogIo1xoyhvzLYtiSjIYa6GYRhG9RRqiMkwDMOoBhMIwzAMIxQTCMMwDCMUEwjDMAwjFBMIwzAMIxQTCCNvEJEDXvXLJSLyij+3IYn+00Skv/d+UnX9vcqlrePsX+y9lonIr0WkgdfWVkSeT3DOFiLy3WTsTiXimOLVLPL3XSIiKiLdAvuKRSQXh7AaGcQEwsgndqtqb6+y7RbgezU9kaqed4hF3M5Q1R64mcTH4C0dqaobVPXyBP1aAFkTCOA84MOYciijcBVur/R3qOom4DMROTXD9hk5hAmEka/MxCuwKCIDRWSGV5xthj+LVkQaicgEEVkkIs8CjfzOQe9ARF70is8tTbYAnVcldixwsYi0EpFO4q3fISIneutQLPRs6ALcAxzr7fu9iDQRkbdFZIHnkYzw+nby1jL4m2fXZG8WOCJynIi8JW59iwUicqy3/zYRmet9V7w1L0YDLwV+hya4yY7XExAIjxe9441CRVXtZa+8eAE7vW0RrsDZud7nZkBd7/3XgH9773+ImyUP0BMox61lALAWaO29b+VtG+FKRxwee0yMHVX2AwuBQUAnYIm378/AaO99fe/8X7V7++vi1n4AaA2swhWT7OTZ29trew64yns/G7jEe98QaIwre/Go17cO8B9gSIjtnwJNA5+vAv7uvZ8B9A20tQMWZ/vvbq/svQqy1IaRtzTyynl0wi1M9Ka3vzmulEYXXFXeet7+IcADAKq6SEQWxTnv90XkEu99B6ALUJqkbWEVgmcCd3rrPLygqh97ZSpi+93tVUU9iLspH+m1faKqC73384FOXi2kdqo60buuPQAicg5OJD7wjm/iXcf0mO9rpW7dDp9RuDLj4IoEjgIWeJ83Am0TX7ZRm7EQk5FP7FZXsrsj7oncz0H8CpiqLjdxIe6p2idhLRkRGYbzOk5W1V64G2zDRH1CztEUJ1org/tVdTyunPNu4A0ROTOk+2igGOjnXdsXge/fGzjuAM7bCBMivP2/VZej6a2qx2l43a5yEanj2X04cCbwmFdF9jZgpFSoWEPPdqNAMYEw8g5VLQO+D9zqlRpvjlu2EeDawKHT8WLo3voPPUNO1xzYqqq7vFE8g5OxxYvh/wV4UVW3xrQdA6xR1Qdw1YJ7AjuApjHfv1FV94vIGTjxi4u65HKJiFzsfUcDr/DbG8C3PXsQkXYickTIKVbgkuoAl+OWs+2oqp1UtQPwCXCa196VHKvWamQWEwgjL1HVD4APcYnVccBvReR9XH7C52GgiRda+hEwJ+RUrwN1vWN+BURdp3qql4yeA/yX8KqcI4ElXlisG+5mXAq87w3V/T3wNNBfRObhxCxKGfKrcWGxRbi8QRtVnQyMB2aKyGLcOghNQ/q+iqsQDC6cNDGm/d/AN733Z3jHGwWKVXM1jAJCRI7CCdXZEY6dDoyI9YyMwsE8CMMoIFT1M+BvwYlyYYhIMXCfiUNhYx6EYRiGEYp5EIZhGEYoJhCGYRhGKCYQhmEYRigmEIZhGEYoJhCGYRhGKP8/M9FxwKzsLI4AAAAASUVORK5CYII=\n" }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "print(solute_scipy_peak.radii)\n", "solute_scipy_peak.plot_solvation_radius('Li', 'PF6')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fantastic! We have implemented our own function that can find the tricky **PF6** minima.\n", "\n", "Now you know how to write your own custom functions and wrappers to calculate solvation radii however you like. The possibilities are endless." ] } ], "metadata": { "kernelspec": { "name": "solvation_analysis", "language": "python", "display_name": "solvation_analysis" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }