{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Relaxing an Al (1 1 1) surface under an electric field using SPHInX"},{"metadata":{},"cell_type":"markdown","source":"If you make use of this notebook and/or its recipes, please consider citing the following (for which this notebook has been distributed as supporting material):\n\n> Christoph Freysoldt, Arpit Mishra, Michael Ashton, and Jörg Neugebauer, *Generalized dipole correction for charged surfaces in the repeated-slab approach*, Phys. Rev. B **102**, 045403\n\nThe goal of this notebook is to give a very simple example for using the generalized dipole correction (outlined in detail within the above paper) within the [SPHInX](https://sxrepo.mpie.de) DFT program.\n\nAs one of the simplest use cases, we will add an electric field to both sides of a very thin Al (111) slab and relax the slab. The electric potential in the DFT cell will be plotted at the end of the notebook so you can visualize the generalized dipole correction as it is described in the paper.\n\nThe calculations and analysis are performed using [pyiron](https://pyiron.org). For more details on specific commands/options in pyiron or SPHInX, please consult their individual documentation and tutorials."},{"metadata":{"trusted":true},"cell_type":"code","source":"from pyiron import Project\nfrom pyiron_atomistics.sphinx.base import Group\n\nimport os\n\nimport numpy as np\n\nimport matplotlib.pylab as plt","execution_count":1,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Create the job"},{"metadata":{"tags":[],"trusted":true},"cell_type":"code","source":"pr = Project(\"ChargedRelax\")\njob = pr.create_job(\n job_type=pr.job_type.Sphinx,\n job_name=\"Al_111\"\n)","execution_count":2,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Build an Al (1 1 1) surface and assign it to the job"},{"metadata":{"trusted":true},"cell_type":"code","source":"job.structure = pr.create_surface(\"Al\", \"fcc111\", size=[1,1,4], vacuum=10)\njob.structure.add_tag(selective_dynamics=(True, True, True))","execution_count":3,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"job.structure.plot3d()","execution_count":4,"outputs":[{"output_type":"display_data","data":{"text/plain":"","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"fdb49eda1ce64002b2882c2e8930210f"}},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"NGLWidget()","application/vnd.jupyter.widget-view+json":{"version_major":2,"version_minor":0,"model_id":"9566ca88f97b4654855b12ac030f8e42"}},"metadata":{}}]},{"metadata":{},"cell_type":"markdown","source":"### Freeze the bottom half of the slab"},{"metadata":{"trusted":true},"cell_type":"code","source":"job.structure.selective_dynamics[\n range(len(job.structure)//2)\n] = (False, False, False)","execution_count":5,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Set up the DFT input for a basic geometry optimization"},{"metadata":{"trusted":true},"cell_type":"code","source":"job.calc_minimize()","execution_count":6,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Modify energy cutoff and k-points for a surface calculation\n\nFeel free to use lower settings to get the calculation to run faster on MyBinder, which is a free-to-use server!"},{"metadata":{"trusted":true},"cell_type":"code","source":"job.set_encut(400) # in eV\njob.set_kpoints([9, 9, 1], center_shift=[0.5, 0.5, 0.25])","execution_count":7,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Charge the bottom layer of atoms\nIn this case, we are targeting a specific field on the right and left sides of the slab (`right_field` and `left_field`, respectively). The following cell calculates the charge required to create this field and then distributes this charge evenly across all \"bottom\" (leftmost) atoms of the slab for the initial density."},{"metadata":{"trusted":true},"cell_type":"code","source":"HARTREE_TO_EV = 27.2114\nANGSTROM_TO_BOHR = 1.8897\n\n# atomic units (1 E_h/ea_0 ~= 51.4 V/Å)\nright_field = 0.05\nleft_field = -0.025\n\ncell = job.structure.cell * ANGSTROM_TO_BOHR\narea = np.linalg.norm(np.cross(cell[0], cell[1]))\n\ntotal_charge = (right_field - left_field) * area / (4 * np.pi) # Eqn. 3 from Freysoldt, et al.\n\npositions = [p[2] for p in job.structure.positions]\njob.input.sphinx.initialGuess.rho.charged = Group({})\njob.input.sphinx.initialGuess.rho.charged.charge = total_charge\njob.input.sphinx.initialGuess.rho.charged.z = min(positions)","execution_count":8,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Add the dipole correction to the PAWHamiltonian group\n\nThe `z_field` parameter in the `PAWHamiltonian` input group controls the \"excess\" field that should remain on the other side of the dipole correction layer, i.e. the field that will exist on the bottom (left) side of the slab."},{"metadata":{"trusted":true},"cell_type":"code","source":"job.input.sphinx.PAWHamiltonian.nExcessElectrons = -total_charge\njob.input.sphinx.PAWHamiltonian.dipoleCorrection = True\njob.input.sphinx.PAWHamiltonian.zField = left_field * HARTREE_TO_EV","execution_count":9,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Run the job"},{"metadata":{"trusted":true},"cell_type":"code","source":"job.run(run_again=True)","execution_count":10,"outputs":[{"output_type":"stream","text":":1: DeprecationWarning: pyiron_base.job.generic.run(run_again=True) is deprecated. It is not guaranteed to be in service in vers. 0.4.0\n job.run(run_again=True)\n","name":"stderr"},{"output_type":"stream","text":"The job Al_111 was saved and received the ID: 1\n","name":"stdout"}]},{"metadata":{},"cell_type":"markdown","source":"### Plot the electrostatic potential\nUncomment the commented lines to plot the expected fields on top of the potential."},{"metadata":{"trusted":true},"cell_type":"code","source":"v_electrostatic = job.get_electrostatic_potential()\nplanar_avg = v_electrostatic.get_average_along_axis(2)\n\nx = np.linspace(0, job.structure.cell[2][2], len(planar_avg))\n\n# ind_0 = np.argwhere(x > 9)[0]\n# ind_dipole = np.argmax((2 * planar_avg - np.roll(planar_avg, 1) - np.roll(planar_avg, -1))[np.abs(x - 12) < 3]) + ind_0\n# m_dipole = x[ind_dipole+1]\n\n# E_right = right_field * 51.4 * (x - m_dipole) + planar_avg[ind_dipole]\n# E_left = left_field * 51.4 * (x - m_dipole) + planar_avg[ind_dipole + 1]\n\n# plt.axvline(m_dipole, color=\"k\")\n# plt.plot(x, E_right, 'k--')\n# plt.plot(x, E_left, 'k--')\n\nplt.plot(x, planar_avg)\nplt.xlabel('$z$ [Å]')\nplt.ylabel('Potential [eV]')\nplt.show()","execution_count":11,"outputs":[{"output_type":"display_data","data":{"text/plain":"
","image/png":"iVBORw0KGgoAAAANSUhEUgAAAYMAAAEMCAYAAAAmgtofAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+lklEQVR4nO3dd3icd5Xo8e9RGfXeJRfZjmuakyiOHYcEUiABQkgCIRBIIZCwtA2wJdtYuPfZvctuFpYOgTQgpACBlE0gxklId2wn7iXutqwuWV2jMvO7f8yMrNiSpr1N0vk8jx5JM+/MHI9H73l/59fEGINSSqmZLcXtAJRSSrlPk4FSSilNBkoppTQZKKWUQpOBUkopNBkopZQC0twOIBalpaWmtrbW7TCUUmpK2bhxY5sxpiyWY6dEMqitrWXDhg1uh6GUUlOKiByK9VgtEymllNJkoJRSSpOBUkopNBkopZRCk4FSSik0GSillEKTgVIqTo1dAxxu73c7DGUx25KBiGSKyBsisllEtovIN8O3F4vIGhHZE/5eZFcMSilrvb6/nfd950W+9PBbboeiLGZny2AQuNgYcyawHLhcRFYCdwJrjTELgbXh35VSHveHt47yqXvW0e0f4XB7n9vhKIvZlgxMSG/41/TwlwGuAh4I3/4A8GG7YlBKJc8Yww+f38sdj2zi7DlF3H7hfI71D9M/NOJ2aMpCtvYZiEiqiGwCWoA1xph1QIUxphEg/L18gsfeJiIbRGRDa2urnWEqpSYwHAhy5++28l9/2s2Hl1fzi1tXsLQqH4CGTr/L0Skr2ZoMjDEBY8xyYBawQkROi+Oxdxtj6owxdWVlMa2zpJSyUI9/mE/fv55HNhzhSxefwnc+tpyMtFSqC7OAUEeymj4cWajOGNMpIi8AlwPNIlJljGkUkSpCrQallIc0dg1wy33r2dPSy7euPZ2PnTtn9L6qgkwAGjo1GUwndo4mKhORwvDPWcClwC7gCeCm8GE3AY/bFYNSKn47Grq5+oevUn9sgPtuPvcdiQCgsiATETiqZaJpxc6WQRXwgIikEko6jxpjnhKR14BHReRW4DDwURtjUErF4S9vt/KFB98kLzON33xu1Wj/wFjpqSlU5GVqy2CasS0ZGGO2AGeNc3s7cIldr6uUSszDbxzmn/6wjUUVedx387lUhstB46ku1GQw3UyJzW2UUvYxxnDXs7v54fP7uGhRGT+84WxyMyY/NVQVZrGjoduhCJUTdDkKpWawwZEAdzyyiR8+v4/rz53Nz2+qi5oIAGoKszjaOYAxxoEolRO0ZaDUDNXVP8xtv9zAugMd/O37FvP5dy9ARGJ6bHVBJkMjQdr7hijNzbA5UuUETQZKzUBHOvq5+b43ONIxwHevX85Vy2vienxkrkFD54Amg2lCy0RKzTCbj3Ry9Y9eoa13iF/euiLuRABjk4EOL50utGWg1Azy7PYmvvzwW5TlZfDwzSs4pTw3oecZ2zJQ04MmA6VmiPtfOcA3n9rBGbMK+fmNdZTlJV7eKcpOJzM9RZPBNKLJQKlpLhg0/PvTO/n5ywe4bFkF37v+LLJ8qUk9p4hQXZhFg65PNG1oMlBqGhsYCvCVRzbxx+1N3Hx+Lf/ywWWkpsQ2Yiia0PBS7TOYLjQZKDVNtfcOcusDG9hc38nXP7iMT18wz9LnryrIZHeTLi8/XWgyUGoa2t/ay833rae528+PbziHy0+rtPw1qguzaOkZZHAkQEZacmUn5T4dWqrUNLP+YAfX/PhV+gZHeOi2lbYkAjg+oqi5a9CW51fO0mSg1DTy1JYGbvj5OoqyfTz2+fM5e06Rba9VE04GR3VE0bSgZSKlpgFjDD99cT//8cwu6uYW8bMb6yjK8dn6mpFNbnTHs+lBWwZKTXEjgSD/8vg2/uOZXXzgjCp+9ZnzbE8EcLxMtHZnC+29Wiqa6rRloNQU1j80wpd+/RZrd7Vw+0Xz+fv3LSHFoqGj0WSmp3Jd3Swe3VDPmp3NfHh5NbesnjfuhjjK+2QqLEFbV1dnNmzY4HYYSnlKS4+fW+/fwPaGLv7PVafxyZVzXYljT3MP9796kMfePMrAcIBV80u4ZXUtlyytsGxOg0qMiGw0xtTFdKwmA6Wmnr0tPdx833rae4f4wSfO4pKlFW6HRFf/MA+vP8wDrx6kocvPnOJsbjq/luvqZpGXme52eDOSJgOlprHX97dz2y824EtL5d6b6zhjVqHbIb3DSCDIn7Y3c98rB9hw6Bg5vlQ+Wjebm8+vpbY0x+3wZhRNBkpNU09sbuBvHt3M7OIs7r9lBbOLs90OaVJb6ju5/5WDPLmlgZGg4ZIl5dyyeh7nLyiJeSMdlThNBkpNM2OHjq6YV8zPPlVHQfbUKb20dPv51euHeHDdYdr7hlhckcctq2v58Fk1ZKbr7GW7eCIZiMhs4BdAJRAE7jbGfFdEvgF8FogsavKPxpinJ3suTQZqJgsEDd98cju/eO0QHzyjiv++7swpu/yDfzjAk5sbuPeVg+xs7KYoO51PnDeHT62spTI8b0FZxyvJoAqoMsa8KSJ5wEbgw8B1QK8x5q5Yn0uTgZqpBoYCfPnht1izo5nbLpzPnZc7N3TUTsYY1h3o4N6XD7BmZzOpIlxxehWfXl3LWTbOmp5p4kkGts0zMMY0Ao3hn3tEZCcQ//56Ss1QHX1D3PrAejYd6eRfr1zGLautXXXUTSLCyvklrJxfwuH2fh547SCPrj/Ck5sbOGtOIZ9ePY/LT6skPVXnxTrFkT4DEakFXgROA74K3Ax0AxuArxljjk32eG0ZqJnmcHs/N933Bg2doQ3rLz+tyu2QbNc7OMLvNtZz3ysHONjeT1VBJp9aNZePnzvHkRnV05EnykRjgskF/gL8mzHmMRGpANoAA/xfQqWkT4/zuNuA2wDmzJlzzqFDh2yNUymv2Ha0i5vvW89wIMg9N9VRV1vsdkiOCgYNz+1q4b5XD/DK3nYy0lK45uxZfHp1LQsr8twOb0rxTDIQkXTgKeBPxphvj3N/LfCUMea0yZ5HWwZqpnhlbxu3/3Ij+Zlp/OLWFZxSPrNPfruberjvlQM89tZRhkaCXLiojE+vruXChWXTou/Ebp5IBhIaRPwA0GGMuWPM7VXh/gRE5CvAecaY6yd7Lk0GaiZ4aksDX3lkE/NLc3ng0yt0dM0Y7b2D/HrdYX75+iFaegZZUJbDLavncc3ZNWT7dIm1iXglGVwAvARsJTS0FOAfgY8DywmViQ4Ct0eSw0Q0Gajp7pevHeTrT2ynbm4RP7/x3Ck1h8BJQyNBnt7ayL2vHGBLfRcFWaGhqTeumktVQZbb4XmOJ5KBlTQZqOnKGMP31u7lO39+m0uXlvODT5ytk7BiYIxh46Fj3PPyAf60vYkUET5wRhW3XjDPc8tzuMkTQ0uVUpMzxvDvT+/kZy8d4NqzZ/Gta08nTYdSxkREqKstpq62mCMd/dz/6kEeWX+Exzc1sGJeMZ+5YB6XLq3QfoU4aMtAKRcEgoZ//sNWHnrjCDefX8vXP7hMT1xJ6vEP88j6I9z3ykGOdg4wrzSHWy+Yx0fOmTVjW1taJlLKwwJBw9/+ZjOPvXWUL77nFL723kW6aJuFRgJB/ri9iZ+9uJ/N9V0U5/i4cdVcblpVO+PmK2gyUMqjRgJBvvabzTy+qYGvXbaIL12y0O2Qpi1jDG8c6ODuF/ezdlcLWempXL9iNp991/zRLTunO00GSnnQSCDIVx7dzJObG/i7yxfz+Xef4nZIM8buph5++uI+ntjUgAh85JxZfO6iBcwtmd77K2gyUMpjgkHD3/1uC7/dWM+dVyzhcxctcDukGan+WD93v7ifh9cfIRA0XHt2DV+6eKHn94VIlCYDpTzEGMM3n9zB/a8e5I5LF3LHpYvcDmnGa+n28+O/7OPBdYcJBg3XnTubOy5dSHne9Jrop8lAKQ/5zpq3+e7aPdx6wTz++QNLtbPYQ5q6/Pzohb38et1hfGkp3H7hAj574bxpM6s5nmSgg5qVstFvNhzhu2v38NFzZmki8KDKgkz+z1WnsearF3HRojK+8+e3uezbL7JmR7PboTlOk4FSNnltXzv/+PutrD6lhH+/5nRNBB42rzSHH3/yHB69fRW5GWl89hcb+OwvNtDc7Xc7NMdoMlDKBvtbe/ncrzYytySHH91wjm7SMkWsmFfMU1++gH+4Ygkv7Wnl8v95kWe3N7kdliP0E6qUxfzDAf7qV2+SmiLce9O5FGTponNTSXpqCrdftID//fK7qCnK4rZfbuSffr+VwZGA26HZSpOBUhb75pM72N3cw3c+tpw5JdNzyOJMsKAsl8f+ajW3XzifB9cd5lP3vEFn/5DbYdlGk4FSFnpycwMPvXGYv3r3Ai5aVOZ2OCpJvrQU/uH9S/nu9cvZdLiTa370Kgfb+twOyxaaDJSyyJGOfv7hsa2cPaeQr16mcwmmk6uW1/DgZ8/jWP8QH/nJaxxqn34JQZOBUhYwxvDPf9gGwPc+fpZ2GE9D59YW85vPrSIQDPLJe9ZNu5FG+olVygL/u7WRv7zdytfeu4hZRdpPMF2dUp7H/besoL13iBvveYOu/mG3Q7KMJgOlktQ1MMw3n9zB6TUF3Liq1u1wlM3OnF3I3Z+q40BbH1986E2CQe+v4hALTQZKJemuP+2mvXeQf7/6dFJ1g5oZ4YKFpXz9ymW8tKeNe1854HY4ltBkoFQS9jT38Kt1h7hxVS2nzypwOxzloBvOm8Nlyyr41h93se1ol9vhJE2TgVJJ+N5ze8lOT+XLuknNjCMifOvaMyjK9vHXD79F/9CI2yElRZOBUgna09zDU1sauOn8Wopn2HaKKqQ4x8e3r1vOvtY+fvLCPrfDSYptyUBEZovI8yKyU0S2i8hfh28vFpE1IrIn/L3IrhiUslOkVfCZd813OxTlogsWlvKB06v4+csHaOmZusNN7WwZjABfM8YsBVYCXxCRZcCdwFpjzEJgbfh3paYUbRWosb723kUMjgT5/tq9boeSMNuSgTGm0RjzZvjnHmAnUANcBTwQPuwB4MN2xaCUXX7wvLYK1HHzy3K5/tzZPPTG4Sm7XIUjfQYiUgucBawDKowxjRBKGED5BI+5TUQ2iMiG1tZWJ8JUKiZtvYM8vbWRj507R1sFatRfX7KQ9NQU7np2t9uhJMT2ZCAiucDvgDuMMd2xPs4Yc7cxps4YU1dWpgt+Ke/43cZ6hgOGT5w32+1QlIeU52dy6wXzeGpLI/tae90OJ262JgMRSSeUCB40xjwWvrlZRKrC91cBLXbGoJSVjDE8sv4IdXOLOKU8z+1wlMfcdH4taSnCQ+sOux1K3OwcTSTAPcBOY8y3x9z1BHBT+OebgMftikEpq6070MH+tj6uXzHH7VCUB5XlZfC+0yr57Zv1+Ien1mY4kyYDEemO8tUjIm9P8PDVwKeAi0VkU/jr/cB/AJeJyB7gsvDvSk0JD79xmLzMND5wepXboSiPuuG8OXT2D/P01ka3Q4lLWpT79xljzprsABF5a7zbjTEvAxMt1HJJDLEp5Smd/UM8va2Jj9XNJsuX6nY4yqNWzS9hfmkOD647zDVnz3I7nJhFKxNdG8NzxHKMUlPe45saGBoJcv0K7ThWExMRPnHeHDYeOsbOxpjHzLguWjL4qoisnuwAY8x+C+NRyrP+uK2JRRW5nFqtC9KpyX3knFn40lJ4cN0ht0OJWbRksAe4S0QOisi3RGS5AzEp5TldA8OsP9jBJUsr3A5FTQGF2T6uOK2Sp7Y0Epgi+x1MmgyMMd81xqwCLgI6gPvCaw19XUR0k1c1Y/zl7VZGgoZLl447R1Kpk7x3WSWd/cO8efiY26HEJKahpcaYQ8aYb4U7kz8BXE1oeQmlZoS1O5spzvGxfLauq6hi865FpaSlCH/e2ex2KDGJKRmISLqIXCkiDwLPAG+jHcfKAsYYXny7lY/8+FWu+sHLtPYMuh3SSUYCQV7Y3cp7FpfrTmYqZvmZ6Zw3v5jndk6NebXR5hlcJiL3AvXAbcDTwAJjzMeMMX9wID41jdUf6+faH7/Kjfe+QUPnAG839/Kxn75GY9eA26G9w4ZDx+gaGNYSkYrbxUsq2NPSy+H2frdDiSpay+AfgdeApcaYK40xDxpjpuaSfDOIMYb/+tMuLvqv53l0wxHPdmD91592s6uph3+7+jSe/9t388tbV9DSM8h1P32NIx3e+eN5blcLvtQU3rVI18hS8YlcQKzd5f1SUbQO5PcYY35mjOkQkQtE5BYAESkTkXnOhOgde5p7uOtPu7nwP5/niu++RO+gN7e5+86at/nh8/sYHA7yd7/dwpXff5nNRzrdDusdDrX38eTmBj61ci43nDeXjLRU6mqLefAz59HZP8w3n9zudoij/ryzmfPmF5ObEW2OplLvNLckhwVlOaydAqWiWPsM/hX4e+AfwjelA7+yKygvenzTUS77zov86IW9VBdmsrupm7//7RaM8dZV9w+e28P3ntvL9efO5pU7L+a71y+nvW+QL/z6TU+1EH7yl/2kpaZw6wXvvKY4c3Yhn1w5l+d2tdDU5f6uUQfa+tjf2selOqRUJejSpRWsO9BOj3/Y7VAmFetCdVcDHwL6AIwxDcCMWbJxOBDkrmd3c1pNPuv+8VIevm0Vf3f5Ev53ayP3vHzA7fBGvbqvjbuefZtrzqrh364+ndQU4arlNXzjylOpPzbAmh1NbocIQFOXn99trOe6ulmU52eedP/1584maODRDUdciO6dXtnbBsBFWiJSCbp4STnDAcNLe9rcDmVSsSaDIRO6BDYAIpJjX0je8/imBo50DPCVSxdRlpcBwO0Xzud9p1bw/57ZxfqDHS5HGPLo+iPkZ6bx79ec/o5RL+89tZJZRVnc+/JB94Ib42cv7SdgDLdfuGDc++eW5LD6lBIeWX+EoMutmdf3t1NVkMnckmxX41BT1zlziyjISueF3d4uFcWaDB4VkZ8ChSLyWeDPwM/sC8s7RgJBfvj8Xk6tzufiJcdHk4gId330TCryMvifP0+0cKtzevzD/HF7Ex88s5rM9HcuopaaItx8fi1vHOxga32XSxGGdA0M8+t1h7nqzGpmF098gr3+3Dkc7Rzgpb3uXU0ZY3h9fwfnzSsmtCK7UvFLS01hxbxiXt/vjYvGicQ66ewu4LeENqpZDHzdGPN9OwPziie3NHCgrY8vX7LwpBNCXmY6HzlnFq/ta6el29369jPbmvAPB7l2glUSrzt3Njm+VO552d2lpF7Y3cLAcIAbVs6d9Lj3nlpBUXY6D7/h3iYh+1r7aOsdZOX8EtdiUNPDyvklHO7op6HTW8Omx4p5cxtjzBpjzN8aY/7GGLPGzqC8IhA0fP+5vSypzOOyCToQP7S8mqCBp7a4u3b57zbWM680h7PnFI57f35mOtedO5untjTS7GLiWruzhZIcH8tnF056XEZaKh85ZxZrdjS7NhHt9f3tAJoMVNJWzi8GYN2BdpcjmVi0SWdPRXuCWI6ZqtYf7GB/ax9/9e4FpEww8/SU8jyWVeXz+OYGh6M77khHP+sOdHDt2TWTljNuPr+WgDH8xqWO2dBM3hbesyS2mbzX1c1mJGh4Zps7iXbdgQ4q8jO0v0AlbUllPvmZaby+z7ulomgtgwtE5IlJvp4EljkRqBv+vKMZX2pK1JUqr1pezeYjnRxsc2c+3mNvHkUEro6ykcbckhzOqCnguV3udGRtPHSMbv8IlyyJbSbvKeW5zCnO5sW3W22O7GSh/oJ2Vs4v0f4ClbTUFGHFvBJPtwyizaK5KobnGLIiEK8xxrBmZzOrFpREnWx05ZnV/L9ndvHk5ga+dMlChyIMMcbw2Fv1rJpfQk1hVtTj3724nO89t4eOviGKc3wORHjc2l0tpKcKFywsjel4EeHCRaU89uZRhkaC+NJs27L7JPvb+mjt0f4CZZ2V84v5885mGrsGqCqI/rfqtGgzkP8Sw9drTgXrpH2tvRxq7+fSZdEnG1UXZrGitpg/bDrq+CS0/W19HGrv54NnVMd0/HuWlGMMvLTH+avttTubWTm/hLzM9Jgfc9GicvqHAmw45GzzWvsLlNUin6V1Hh1V5Nyl1hSzZkeolBLr4mQfWl7NvtY+dji8zV1kUtS7YrzaPqOmgJIcH887XCo62NbHvta+dwzPjcWqBSWkpwp/cbhU9Pr+DsrzMqjV/gJlkaVV+eRlpnm2VKTJYAJ/3tnMaTX5MTfnLj+tEsDxk9bLe9qYU5w96Zj9sVJShIsWl/GXt1sdXZ5ibTj5XLIkvmUdcjPSOGduES++7dx8A2MM67S/QFksNUU4z8PzDTQZjKOtd5A3Dx/jsqWVMT+mNDeDRRW5jv5HjwSCvLa/ndWnxNYqiHjP4nKO9Q+zub7TnsDG8dyuZhaW5zIngSvtixaVs7Ox27G5HA1dflp6Bqmr1Y1slLXOm1fCgbY+V4d3TyTa0NKtIrJlnK+tIrIlymPvFZEWEdk25rZviMhREdkU/nq/Vf8QKz23qwVj4NJl8ZU0Vs4vYcPBDoYDQZsie6etR7vo8Y9wQZzJ4MKFZaQIvOBQqWg4EGTDwWMxdxyf6MJFoce96NDaLlvDSfKMWYWOvJ6aOc4LzzfwyhI2Y0VrGXwQuHKcr8jtk7kfuHyc279jjFke/no6vnCd8ecdzVQXZLKsKj+ux62cX0L/UICtR51Z8uGVvW2IhOrq8SjITuecuUU859BaKW839zA4EuSsOYldaS+ryqcsL8OxEtyW+i7SUoQllTNmLUblkMWVeaSlCNsbnO1bjEW00USHJvuK8tgXAe+lvyiCQcNr+9u5aHF53PXiFfNCWT8yEsVuL+9t49Tq/ISGiL57cTnbjjpTetl8JJQclyd4pS0ivGthKS/vcaafY+vRLhZV5J20xpNSycpIS2VhRd7USwYRIrJSRNaLSK+IDIlIQEQS/dd8MVxquldEJrxUFJHbRGSDiGxobXWuU/Zgex89/hHOirJcwnic7DfoHxrhzUOdcfcXRERaE285sOnN5iOdFGWnM7s48bHV5y8o5Vj/MPtbey2M7GTGGLbUd3HGrAJbX0fNXKdW57Ojoctze6HE2oH8A+DjwB4gC/gMkMhCdT8GFgDLgUbgvyc60BhztzGmzhhTV1bm3FryW8Krep4xO7GTgVP9BusPHmMoEIy7vyBiaWU+KQI7HLhC2VzfyRmzCpMamXNqdahkZ/fQ3fpjA3QNDHO6JgNlk1Or82nrHaLFpTW3JhLPQnV7gVRjTMAYcx/wnnhfzBjTHH58kNAS2CvifQ67banvIjM9hVPKchN6vFP9Bq/sbcOXlsK5tcUJPT7Ll8q80hzbT659gyO83dzDmQm0tMZaUJZLeqqws7HHmsAmMHoxUFNo6+uomevU6tCFxvYGd5eTP1GsyaBfRHzAJhH5TxH5ChD3BjciUjXm16uBbRMd65Yt9Z2cWl1AWmpio26d6jdYf7CD5bMLk6prn1pdYHvLYNvRLoIGlifY0orwpaWwsDzP9uS15WgnvtQUFlUmdjGgVDRLq0IDE5xolccj1jPep4BU4IuEtr6cDVw72QNE5CHgNWCxiNSLyK3Af44Zlvoe4CsJR26DkUCQ7Q3dSdWLI/0Gdk45N8awt6WXxRXJjXZZVp3P0c4Buvrt25t1s4XDNJdW5bPT5mSwtb6LxZV5ZKRp57GyR15mOrUl2Z7rRI62UB0QGlUU/nEA+GaMj/n4ODffE2Ncrtjb2svAcCDpzsOV80v47cZ6hgNB0hNsYUymtXeQHv8IC8qS2300MnR2R2N33MNTY7W5votZRVmU5mYk/VzLqvP53Zv1tPYMjm4/aqVg0LD1aBdXnhnbOk9KJWpZdT7bjnorGUSbdPZo+Pu4k8+cCdE5o/XiJK9iz5xVSP9QgCMd/RZEdbJ9LaGlsheUJ1fKWBpOBnbWLjcf6Uy6vyAi0ry2q3VwqKOfHv8IZ9Ro57Gy16nVBRzu6Kfbb1+rPF7RLlv/Ovx9osln08qW+k7yMtKYV5LcFXdkM5RD7TYlg/DwygUJdnJHlOVlUJ6XYVsdvq13kPpjAwnPLzhRpCVjVzKIdPrrSCJlt2WR0XEeKhVFm3QW2WLq8+NMOPu8/eE5a2t9F6fVFEy4q1ms5oaTyaF2eza72dfaS7Yvlcr8zKSfa1l1vm0fyC3h/gKrWgaF2T6qCzJtS15b6zvxpaWwKMm+GKWiiQyV9lK/QawF7cvGue0KKwNx29BIkJ2NPZZMNirN9ZHtS+WQXWWi1j7ml+UknbQgdLW9t6WXwZGABZG906YjXaQInFYT37Iek7GzE3lHYzdLK/Ns6edRaqzyvEzK8jI8Nbw0Wp/BX4nIVkIjgsb2FxwAplWfwe6mHoYCQUtGvYgIc4qzOWxXmailN+kSUcSy6nxGgoY9zdbP7N3X0kttSQ7ZvpjGKcRkWXU++1r78A9bn7zqjw0wJ8kSoVKxOtXGVnkiol0C/ZpQ38ATvLOv4BxjzCdtjs1Rx4dAWlMvnluSbUvLYGAowNHOgYQnxZ1o7Igiqx3tHKA6hq0447G0Kp+ADckrGDQ0dvqpLky+9KZULE6tzmePTa3yRETrM+gyxhwMDxOtB4YBA+SKyBwnAnTK9oZuCrPTmVVkzclrbkkOhzv6CVq8sNr+tnDncZIjiSLmluSQ7Uu15QqlsWvA8pOrXZ3IbX2DDAWCMe0jrZQVaktyCAQNzV3eWJYipva7iHwR+AbQDEQW3THAGfaE5bz6Y/3MKc62bGeruSXZDI0Eaer2W3p1vK81PKzUopZBanipZqtbBkMjQVp6Bi3f+HtOcTY5vlTL423oDK3eWu3BjcrV9BT522jsGkho0yerxVrMvQNYbIzx5uadFmjq8jM/yUlcY80tjowo6rc2GbT0kiLHh69aYWlVPk9ubrDs+QCau/0Yg+VX2ikpwmIbkldD5wCA5WUtpSZSWRBqNTd5ZNezWIdNHAG80+1tg6Yuv6VXsZGT9eEOa4eX7m3tZXZxtqVr7c8uzqbbP0Lf4Ihlzxk5uVbZUIOfV5pr+YS+SLxaJlJOiSSDxi5vJINYWwb7gRdE5H+B0QKXMebbtkTlsB7/MD2DI1QVWHfiqirIJC1FLJ94ZuVIooiK/NDSDk3dfsueu6HLvivtyoIMWnoGCQaNJcNrIdTZneNLJT/LupFPSk0mNyONvMw0mjySDGJtGRwG1gA+IG/M17QQ+c+otDAZpKWmMLs429JkEAgaDrT1Jb0m0YkqwpPXrNyk284afEV+JoGgoa3Puo63hvDIJ6v6jJSKRVVBJo3hCye3xbpQ3TcBRCTHGGPPtFoXNYSTgdVXsXOKszlkYZmooXOAwZGg5S2DSluSwQBF2elk+axf/TOSvFq6BynPsyaBN3Ra29GvVCwq8jOnVstARFaJyA5gZ/j3M0XkR7ZG5qCmcGa2YnmHseaWhFoGVm1vt7fV2mGlEZGTa5OFQ9wabJhjEFE5Gq+1yUuTgXJaqGUwhZIB8D/A+4B2AGPMZuBCm2JyXGOXH5HjJ0WrzCnOpsc/QqdF+wXsa7FmgboT5WSkkZeRZmnLoNHiDvmxRstaPdbE6x8O0N43RI1OOFMOqyzIorV30PZtcmMRz7aXR064yRvT5izQ2OmnNDcDX5q1a9KMLlhn0ciXpi4/WempFGWnW/J8Y1UUZFqaDI52Dth2ci3N9ZEi0GzRFVWjTWVCpaKpKsjEGDyxH3LMQ0tF5HzAiIhPRP6GcMloOmjs9ls6kiiidnQpa2v6DVp7Q5u62NHJWZGfYdl45x7/MD3+EapsOrmmpaZQmptBc7c1f0A6x0C5ZXSugQc6kWNNBp8DvgDUEFqWYjnTaAnrxs4BW5LB7GJr9zWwa4cvCJVeWiw6uTpxpV1ZkGlZ8jqqcwyUS6oKrO+vS1SsyWCxMeYGY0yFMaY8vEjdUjsDc5LVE84iMtNDew5Ymgws2D5yPJX5oTKRFWspRU6u1TYk2IjyPOvKWg2dA7b0GSkVTVX+8SUp3BZrMvh+jLdNOXZMOBtrTkm2ZbOQ23rtbRmMBA3tfUNJP1djpxMtgwxLk0F5nvV9RkpFk5+VRlZ6qieGl046z0BEVgHnA2Ui8tUxd+UD1g8gd4EdE87GqinMYsOhjqSfZ2gkyLH+YVuTAYTmGiT7Gg2dA6QIlNsUK0BFXibH+ofxDweSXpqjodO+kU9KTUZEQsNLPbA+UbRLIR+QSyhpjJ153A18ZLIHisi9ItIiItvG3FYsImtEZE/4e1Fy4ScvUt+262RQmuujrSf5q+328Gxbu5JBJBlacbXd0DVAZX4maTbuGFYRjrfVglEYDZ0D2l+gXFNZ4I2JZ9H2M/hLePbxyvD3bwP/bYz5tjFmT5Tnvh+4/ITb7gTWGmMWAmvDv7sqUquzq0xUkpvBwHAg6UXgIic9u/oMxq5PlCwnJnCNTpRLMl5jTHgTHu0vUO6YEslgjDwReQvYBmwXkY0ictpkDzDGvAicWB+5Cngg/PMDwIfjiNUWdk04iygNn7zbe5NrHYwmA5taBmW5GaGx+xaMKGro9Ns2rDTCqiU0OvqGGBwJ6rBS5Zqq8ByfgMUbYcUr1mRwN/BVY8xcY8xc4Gvh2+JVYYxpBAh/L5/oQBG5TUQ2iMiG1tbWBF4qNnZNOIsoyfUBoTkCybA7GYyO3U/yCiUYNDR12b99pFVLUjQ40Nmt1GQqC7JCgzeSPEckK9YzYI4x5vnIL8aYFwBbdw43xtxtjKkzxtSVlZXZ9jp2TTiLiJR12ixKBpHkYoeK/OTH7ke2j7R7x7D8rDQy0lKSnrmpcwyU26ryvbGvQazJYL+I/IuI1Ia//hk4kMDrNYtIFUD4e0sCz2Gppi57JpxFRE7eSZeJegcpzE4nI82+QVwV+cmP3XdiWCmERmFYUWvV2cfKbV7Z5CbWZPBpoAx4LPxVCtySwOs9AdwU/vkm4PEEnsNSjTYPKyzJsa5lUGpT53GEFWP3R3c4szHBRlTkJd+SaekZxJeaYst6T0rFosojS1JEm2eQSWgpilOArcDXjDExLcEpIg8B7wZKRaQe+FfgP4BHReRWQhvmfDTx0JMXmXBm1xwDAF9aCgVZ6UnXA+2cfRxhxdj9SN9Ieb69sUJoeOnW+s6knqO9d5DiHJ9uaqNcU5zjw5ea4vpcg2ib2zwADAMvAVcQWoLijlie2Bjz8QnuuiTW4OzWNDrHwN6r2JJcH20WlInOnFVoTUATGDt2P7KuUrwi5bDibPv6NiIq8jJY0+3HGJPwyby9b8jWfhilorGq5JmsaMlgmTHmdAARuQd4w/6QnGP3hLOI0twMS0YT2TWSKKJyzNj9hJNB3yBF2em2TjiLqCzIxD8cpNs/QkFWYmWe9t5BSmxucSkVTWVB5mh/m1ui/cWOloSMMcnNmvKgSH3c6h3OTlSa60uqTNQ3OEL/UMD2ZFBhwXDN9t4hinOcudIut2CuQVvvEKUOxavURMryMpLuV0xWtGRwpoh0h796gDMiP4tItxMB2imyKFtpnr0ng9LcjKTKRJEPid19BlZM5AqVXZy50rYm3kEtEynXleb4XE8Gk5aJjDHTYjG6ibT3DpKVnkq2L1q1LDklORl0DQwzNBJMaHKb3RPOIvKz0shMT0nu5No7yOLKPAujmtjoEhoJtmT6h0bwDwe1TKRcV5KbQbd/JOFzhBVm9Jq97b3OdB5GWh4dCS4P7VQyEJHwxLPEr1Da+5wrE0XKWolOPIt0dpdomUi5LHIeOtaf/KKWiZrRyaDNoZJGaZKzkCOdz3YnAwgtO92W4Ml1JBCks394dG6F3TLTUynISk+4JRP5/7B7/oZS0Vg1HykZMzoZtPcOOtJ5WBrO+gkng55BUlOEIgeGaxbn+EaXy45XR/iqptTBGnyocz6xq6nRloH2GSiXWbVSQTJmeDJwpqRxvGWQeJmoJMdHaor9E6NKcjOSPrkWO9QygNAVVaJJNpL0tM9AuS1Sqky0lGyFGZsMjDHhkST2nwhKRpexTrxl4ESJCEKjGo71DyW0nG7kg+zklXZJri/hrTrbtM9AeUSJRQtaJmPGJoNu/wjDAeNISSPHl0pmekpSfQZOJYOS3AyCBjoT6Mg6XoN3OBkk2jLoHSI3Iy3pbTOVSlZ+ZhrpqWLJHuSJmrHJwMmrWBEJlzMSLxPZPccgIlI2S+RD6VaZ6Fj/MCOBYNyP1TkGyitEJNRfpy0D50XedKdGvpQmOMMwGDS09Q5S6ljLIPGOrI6+IVIEChNcGiIRkVZIRwItmfbeIS0RKc8oyUm8v84KMzYZtDk8kqQswcXqugaGGQ4Yx1oGo9t0JjCiqL1vkOKcDFIc6OiOiNRaE+l4a+sddLQVo9Rkkun/ssKMTQaRk51TY8xDWT/+E6yTcwxgTJkogcTV5sKVdkkS8bb3DTnav6HUZEpzMxIe1m2FmZsMwicPJ8buQ2gWcnvfEME4R+lEJoA5lbSKsn2IJDbyqcOF5aATHYURDBpX4lVqIiU5ic+ZscIMTgaDFGSlO7YOSElOBoGgoXMgpr2BRh1vGThz0kpNEYqzE2uuurEcdGmCfRzd/mECQeNYn5FS0RTn+ugfCjAwFHDl9WdsMmjrc7akEekAjveKO3KSc3LJhJIEZ/W60SGbn5lOWorE3bx2us9IqWhKcxLvr7PCjE0GoatYB5NB+CQZ7yY3bb2DpKUI+ZnOjdBJZEmKwZEAPYMjjieDlJTIkLz4kle7rkukPMbtJSlmcDIYcrREcLxlEO9JK1TXdnqETrxxHp+34fzJtSSB/SLaXZgtrdRkSpIYyWeFGZsMnO48jFwxx9vR2dY76HhduzQn/j6D4xPOnD+5lubG35Jxep6JUtEcP0doy8AxgaCho9+5HbkgNEonNUXiTwZ9Q45NOIsoyT2+GU+sRneNc+FKO5FRGG29Q4hAUbZz5TelJqNlIhcc6x/CGGdPXJHadltPnCetHmeW2R4rcnUfz0Ybo1faLpWJ4u6Y7xukKNtHWuqM/BNQHpTtSyMrPZUOl8pE9u73OAEROQj0AAFgxBhT5+TrH9/hytkTV1lufEtSRFZWdbplMHa4ZmQ3sWjcLBOV5ProCw/Jy/LFtuicLkWhvCjRkXxWcCUZhL3HGNPmxgsfv4p19mRQlpcR12iivqEA/uGg46WXRDqy2vuGSE8V8jOd/0iNHZI3y5cd02Oc2vJUqXiU5GbQ5tKSFDOyjRx5s52+MizLyxjdzzgWkdnHTrdgElmSoj3c0S3i3KiniERqrW0O7WWhVDxKXFy51K1kYIBnRWSjiNw23gEicpuIbBCRDa2trZa+uFv17bLwyqWxLkkxun6S02WiBPZjbe9zZte48RxfdjuOeLVMpDyoJMfn2m5nbiWD1caYs4ErgC+IyIUnHmCMudsYU2eMqSsrK7P0xdt7nV9qGUITnIYDhq4Yl6Ro7XGnBZOflUZaisT1oWx3cZ2feLcVHQ4E6RoY1mGlynMic3yMiX+nwWS5kgyMMQ3h7y3A74EVTr6+G0stw/GVR2O94o5c6Tq1YmmEiMTdkdXeO+jabN54y0THdMKZ8qjSXB9DgSA9gyOOv7bjyUBEckQkL/Iz8F5gm5MxtPe6s3RxZE+CWPsNIsNQ3Si/FOfEt5xue697ZaLIkLxYa61tve7NiVBqMsksIZ8sN1oGFcDLIrIZeAP4X2PMH50MwK2SRuQKP9YRRe19gxRmp5Puwlj40jg24+kfGmFgOODqlXY8G4O09PgB51tcSkUzOpLPhU5kx8cBGmP2A2c6/bpjtfcOcsasQsdfdzQZxNoy6B10rZOzJMfHofb+mI6NtGBKXazBl8Qxh6O5O5QMYp1DoZRTRjdrcqETeUYOLXVrjHl+Zhq+tJQ4ksGQi3X42Gf1NkVOrgXunVxL41iSoqkr9O8qz9NkoLxldNvZGVImcpV/OLTUshsnWRGhLDf2uQZtLnbKFuccn9UbTeRKu9LFK+2SOBara+7xU5Ljc2xjI6ViVZQTGuHoRploxv01tI2uVunSMMg4ZiG71dENY5akiOEEe7zs4m6ZqKMvtiF5zV1+yrVEpDwoIy2Vgqx0WuKYnGqVGZcMmrrCV7EulTRibRkMjYTHwrtVJgrX/2OZa9DU5ScjLYUCh+dtjFWS42M4YOj2Rx+S19Ttp9LFxKXUZKoKMmkMn6ecNOOSQUP4Ta4uzHLl9SOzkKPpGF0S2qUyURxj95t7BqksyHRlKYqI0jhGYTR3D7p2MaBUNKFkMOD46864ZNAUfpNdaxnkZdDeN8RIYPK9AtpcWkwvIp4lKZq7/FS43Bk7OvEsSktmOBCkvW9QO4+VZ1UWZI1WMJw045JBQ6ef3Iw0R/cUHqss14cx0csvbS7v0VuaFzq5xlK7bOr2uzqSCI6/Ty3dk8fb0jOIMe5dDCgVTXVBJu19Q/iHow/esNKMSwaNXQNUuXgiiHXiWWTCV5lLySDbl0ZhdnrU5qoxhmYP1OAjZb+jnZPPjfDCyCelJhO5UIl8Vp0y45JBU5ff1avCWCeeubXnwlg1hVkcPTZ5MugaGGZwJOj6BK6CrHTyMtKixtscbn6Xawey8qjIhU1DpyYDWzV0+akucKfzGKAsN3TSjJYM2noHyUpPJSfDvf2Hqguzon4gmzw0m7emKIujnZMngyZtGSiPi1ysNnU724k8o5LB0EiQtl53R5JEavGxlIncXlWzpjCLhign1+Zwjd4LNfiawizqo7UMugdJTxXXFtVTKppIGdvp4aUzKhk0d/sxBqoL3TtxZfvSyM1Ii6ll4FbncURNYRY9gyN0+yfefyFSdnF7NBHE1jJo7vZTnufuMFilJpPtS6MgK51GLRPZp3F0wpl7ZSIIze6Nngzcm30cMdopO8nVdqTs4oUafE1hFj3+yZOX231GSsXCjYlnMywZhE5q1S6fDGKZeObmZjERkRbUZKWi5m4/RdnpZKanOhXWhGqKoiev5h6/9hcoz3Nj4tkMSwahTFvl0uzjiLK8yZekCE2MGnJ9vf2a0VENkycDL3Qew/F4J00GXX5PtGKUmkxVofMTz2ZUMmjq8pOXEarZuyna+kRHOvoJBA21JTkORnWy0twMfKkpHJ2kdtnkpWQQaRlMkLx6/MP0DQW0ZaA8ryrf+YlnMyoZNHQOUOVi53FEWV4G3f6RCf+jD7b3AVBb6m4ySEkRqgozJ+2Ube4e9MzJtTQnA19ayoTxRkY+eSV5KTWRSPXCyYlnMyoZNHb5Xe88huMTzybqN9jfGkoG811OBgDVBRMPLx0OhIbqurl09VgpKTLpRDnd4UxNFZHhpU5OPJtxycDtzmM4vo7ORHsMH2zvoyArnSIPjIWvnmSuQWt4nR+31yUaq6Ywi/oJ4nV7+XKlYlXlwsSzGZMMBkcCtPUOUuWBlkFkxczGCU5aB9r6XC8RRdQUZdHc7Wd4nFVWvbjOz6Qtgx73N+FRKhaR85S2DGzQHN731s1F6iLml+UgAm839457/8G2fk+UiABqCjMJGsYd2eDFsktNURZtvYPj9sc0d/nJy0wj2+fuAAKlosnypVKYne7oiCJXkoGIXC4iu0Vkr4jc6cRrRsbseqEDOScjjbnF2exu7j7pPv9wgIauAddHEkVUTzK8NPJB9VQymCzebp1joKaOynxn5xo4ngxEJBX4IXAFsAz4uIgss/t1R+cYeKBlALC4Mo9djT0n3X6ovR9jYF6ZN5LB6Ml1nA9lc88gaSni2n7S45lseGlz96CnEpdSk6kuzHJ0FrIbLYMVwF5jzH5jzBDwMHCV3S96PBm432cAsKQynwPtfQwMvbOccaAtNJJonsdaBuPV4Zu7/JTnZZCS4p11fiabeOalCXJKRVPp8JIUbiSDGuDImN/rw7fZqrFrgPzMNFeXhB5raVUexsCelne2Do7PMch2I6yTZKanUpLjG3fimRd2ODtRZUEmKXJyy6B3cITmbj+zirxxMaBUNNUFmXQ4OPHMjWQw3mWkOekgkdtEZIOIbGhtbU36RRs6/Z5pFQAsrswHOKlUdKC1j9LcDPJc2pZzPBMNLz1yrN/VvSHGk56aQmV+5kktg/UHOggaWDGv2KXIlIpPZE6UU53IbiSDemD2mN9nAQ0nHmSMudsYU2eMqSsrK0v6Rb0y+zhiTnE2Wemp7Gx6ZyfygfY+5nmkVRAx3r4GRzsHONIxwNlzi1yKamI1RSfPNXhtfzu+1BTOnuO9eJUaT7XD+xq4kQzWAwtFZJ6I+IDrgSfsfMH23kF2NXVzek2BnS8Tl9QUYdE4ncgH2vo8M5IoorowtE+AMccbcK/ubQNg9SklboU1ofHmGry2r53lswvJ8rm/uqpSsZhflsstq2sdW8re8WRgjBkBvgj8CdgJPGqM2W7naz67o5mggctPq7TzZeK2pCKPXU3doyfZ3sERWnsGPTOSKKK6MJP+oQBdA8f3CXhtXzslOT4Wlee5GNn4aoqyaOr2MzQSmijXNTDM9oYuVi7wXuJSaiKVBZn865WnsrDCmb8xV+YZGGOeNsYsMsYsMMb8m92v98y2JuYUZ7OsKt/ul4rLkqo8jvUPj65getBjI4kiIi2VLfVdABhjeGVfGysXlHhqJFHEubXFBIKGZ7Y1AvBGuL9g1XxNBkpNZNrPQO7qH+bVvW1ccXql57Y6XBLuRN7ZFCoVjQ4r9VjL4IKFpRRlp/PQG4cB2N/WR3P3IKsXlLoc2fguXFjG/LIcfv7SAYwxvL6/HV9aCmfNKXQ7NKU8a9ongzU7mxkJGq44rcrtUE6ypDLU/Nsd7kSOtAzmFnsrGWSmp3Jd3Wye3dFMU5d/tL/gfI+WXVJShFsvmMfWo12sP3iM1/a1c86cIk/sxqaUV037ZPDHbY1UF2Ry5izvdB5HFOX4qMjPYFdjD4GgYdORTqoKMj3ZyfmJ8+YQCBoeXn+YV/e1U12QydwSb416Guuas2ZRmJ3Ofz+7m51N3azyaOJSyiumdTLoHRzhxT1tvO8075WIIpZU5rP+UAcf/cmrrN3V4rlO7oi5JTlctKiMh944zGv72zn/lFLPvqcQWujrk+fNZd2BDoxBk4FSUUzrZPDcrhaGRoKeLBFFLKnM40jHAAfa+vifjy3n6x+0fZmmhH1y5Vyauwfp7B/2bIlorBtXzSU9VchMT+HMWYVuh6OUp3ljbQabvLC7hbK8DM7x4MSoiE+cN4f01BRuXl07uumNV128pJzqgkwauvyc79HO47HK8zO5/cIF9A2N4Eub1tc9SiVNxk4k8qq6ujqzYcOGuB83HAhyuKOfBWW5NkQ1M/1uYz2v7W/nro+e6XYoSqkoRGSjMaYulmOndcsgPTVFE4HFrj1nFteeM8vtMJRSFtO2s1JKKU0GSimlNBkopZRCk4FSSik0GSillEKTgVJKKTQZKKWUQpOBUkoppsgMZBFpBQ4l+PBSoM3CcJygMdtvqsULGrMTplq8MHnMc40xMW0iPyWSQTJEZEOs07G9QmO231SLFzRmJ0y1eMG6mLVMpJRSSpOBUkqpmZEM7nY7gARozPabavGCxuyEqRYvWBTztO8zUEopFd1MaBkopZSKQpOBUkqp6ZMMRORyEdktIntF5M5x7hcR+V74/i0icrYbcY6JZ7aIPC8iO0Vku4j89TjHvFtEukRkU/jr627EOiaegyKyNRzLSVvPefA9XjzmvdskIt0icscJx7j+HovIvSLSIiLbxtxWLCJrRGRP+Pu4e7dG+9w7HPN/iciu8P/970WkcILHTvo5cjDeb4jI0TH/9++f4LFeeo8fGRPvQRHZNMFj43+PjTFT/gtIBfYB8wEfsBlYdsIx7weeAQRYCaxzOeYq4Ozwz3nA2+PE/G7gKbff3zHxHARKJ7nfU+/xOJ+RJkKTcDz1HgMXAmcD28bc9p/AneGf7wS+NcG/adLPvcMxvxdIC//8rfFijuVz5GC83wD+JobPjWfe4xPu/2/g61a9x9OlZbAC2GuM2W+MGQIeBq464ZirgF+YkNeBQhGpcjrQCGNMozHmzfDPPcBOoMateCziqff4BJcA+4wxic5kt40x5kWg44SbrwIeCP/8APDhcR4ay+feFuPFbIx51hgzEv71dcAz+6NO8B7HwlPvcYSICHAd8JBVrzddkkENcGTM7/WcfGKN5RhXiEgtcBawbpy7V4nIZhF5RkROdTaykxjgWRHZKCK3jXO/Z99j4Hom/sPx0nscUWGMaYTQhQNQPs4xXn6/P02olTieaJ8jJ30xXNa6d4JSnFff43cBzcaYPRPcH/d7PF2SgYxz24ljZmM5xnEikgv8DrjDGNN9wt1vEiprnAl8H/iDw+GdaLUx5mzgCuALInLhCfd79T32AR8CfjPO3V57j+Ph1ff7n4AR4MEJDon2OXLKj4EFwHKgkVDZ5USefI+BjzN5qyDu93i6JIN6YPaY32cBDQkc4ygRSSeUCB40xjx24v3GmG5jTG/456eBdBEpdTjMsfE0hL+3AL8n1IQey3PvcdgVwJvGmOYT7/DaezxGc6TEFv7eMs4xnnu/ReQm4IPADSZcvD5RDJ8jRxhjmo0xAWNMEPjZBHF48T1OA64BHpnomETe4+mSDNYDC0VkXvgq8HrgiROOeQK4MTziZSXQFWmGuyFc87sH2GmM+fYEx1SGj0NEVhD6/2p3Lsp3xJIjInmRnwl1Fm474TBPvcdjTHgV5aX3+ARPADeFf74JeHycY2L53DtGRC4H/h74kDGmf4JjYvkcOeKE/qyrJ4jDU+9x2KXALmNM/Xh3JvweO9Er7sQXoZEsbxPq+f+n8G2fAz4X/lmAH4bv3wrUuRzvBYSam1uATeGv958Q8xeB7YRGMLwOnO9ivPPDcWwOx+T59zgcUzahk3vBmNs89R4TSlSNwDChK9FbgRJgLbAn/L04fGw18PSYx570uXcx5r2E6uuRz/NPTox5os+RS/H+Mvw53ULoBF/l9fc4fPv9kc/vmGOTfo91OQqllFLTpkyklFIqCZoMlFJKaTJQSimlyUAppRSaDJRSSqHJQCmlFJoMlLKMiHwmvGzwLW7HolS8NBkoZZ1rgYuBj7odiFLx0mSgVJxEpFZEBsbZWGQdoTWE1o05Niu8wciQR9Y8UmpcmgyUSsw+Y8zyE27LBV4CCiI3GGMGwsd5YcE+pSakyUCpcYjIc2O2F/SLyKSlHxFJIbTY2Y3A1SKS6kigSllEk4FS4zDGXBy+ov8poUXMTlpi/AQXA1uMMQcJLRB2sa0BKmUxTQZKTUBEbiS0F8INxphAlMNv4Pgy2Q+Ff1dqykhzOwClvChcFroBuMoYMxzl2CxC++JeIiL/SegiK09EsowxA/ZHq1TytGWg1AlE5IPA54FrjDH+GB7yIeAZY8wcY0ytMWYO8CRwpZ1xKmUlTQZKnewBQtsbvhLuQL41yvE3ENpacKzfA5+0Izil7KBlIqVOYIwpifP4D41z2+85OUEo5VnaMlAqfgGgYJxJZyeJTDoD0oGgzXEplTDd9lIppZS2DJRSSmkyUEophSYDpZRSaDJQSimFJgOllFJoMlBKKYUmA6WUUmgyUEopBfx/WldazNm8Z/MAAAAASUVORK5CYII=\n"},"metadata":{"needs_background":"light"}}]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.8.8","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":4}