diff --git a/notebooks/tutorial/Screening_example.ipynb b/notebooks/tutorial/Screening_example.ipynb new file mode 100644 index 00000000..74e9620c --- /dev/null +++ b/notebooks/tutorial/Screening_example.ipynb @@ -0,0 +1,382 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b76f62f3-5751-4342-a1ed-69d54a15f35a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "# Screening\n", + "\n", + "Screening is enabled by default for different integral and evaluation calculations in `gbasis`. This helps reduce computation time and memory usage without sacrificing accuracy. This notebook demonstrates the screening feature of `gbasis` on a small molecule (formaldehyde). It includes two examples of screening: one-index screening, where only one contraction is involved in the evaluation, and two-index screening, where two Gaussian contractions are involved in the integration." + ] + }, + { + "cell_type": "markdown", + "id": "b9555f86-81b4-4c8f-a67f-38ae48d4225f", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "## one-index screening\n", + "\n", + "When evaluating a property on a set of grid points, not all points yield non-negligible values for all basis functions. In fact, many points correspond to zero or negligible contributions from a given basis function. This can be estimated by introducing an energy tolerance $\\epsilon$ and calculating a cutoff distance $d_{s,p}$ for each primitive $s$ and each grid point $p$. Then, when evaluating a basis function, the calculation is restricted to points that lie within the maximum cutoff distance over all primitives associated with that basis function. The screening is available for the following functions:\n", + "\n", + "- [`evaluate_basis()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.eval.evaluate_basis)\n", + "- [`evaluate_deriv_basis()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.eval_deriv.evaluate_deriv_basis)\n", + "- [`evaluate_density()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_density)\n", + "- [`evaluate_deriv_reduced_density_matrix()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_deriv_reduced_density_matrix)\n", + "- [`evaluate_deriv_density()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_deriv_density)\n", + "- [`evaluate_density_gradient()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_density_gradient)\n", + "- [`evaluate_density_laplacian()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_density_laplacian)\n", + "- [`evaluate_density_hessian()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_density_hessian)\n", + "- [`evaluate_posdef_kinetic_energy_density()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_posdef_kinetic_energy_density)\n", + "- [`evaluate_general_kinetic_energy_density()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_general_kinetic_energy_density)" + ] + }, + { + "cell_type": "markdown", + "id": "5a9267ea-cc58-42ff-98db-de540aa37a2e", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "### Loading data\n", + "To evaluate properties dependent on the system state, we need its wavefunction or an equivalent representation. In `gbasis`, this is done through the one-electron reduced density matrix (1-RDM). As a case example, the following code will load this data from a `.fchk` file corresponding to a UB3LYP/Aug-cc-pVTZ calculation on formaldehyde. We will be using the [`iodata`](https://iodata.readthedocs.io/en/latest/) package to load the data.\n", + "\n", + "### Defining the grid\n", + "We also need to define an auxiliary set of grid points to evaluate the properties. The loaded system is planar and all of its atoms lie on the $xy$ plane ($z=0$) in this specific case. So for simplicity, we will perform all successive analysis on the $xy$ plane. We define a grid of points in the $xy$ plane spanning the area enclosed by: $-5 \\leq x \\leq 5$ and $-5 \\leq y \\leq 5$. This grid will be used to evaluate the properties of interest, in this case electron density" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "7fa4fa20-9d16-450a-9e19-4e25da319b8f", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from iodata import load_one\n", + "from gbasis.wrappers import from_iodata\n", + "\n", + "## load molecule info and basis from fchk file\n", + "mol = load_one(\"ch2o_q_0.fchk\")\n", + "basis = from_iodata(mol) \n", + "rdm = mol.one_rdms[\"scf\"]\n", + "atcoords = mol.atcoords\n", + "atnums = mol.atnums\n", + "\n", + "## make grid\n", + "grid_1d = np.linspace(-5, 5, 50)\n", + "grid_x, grid_y = np.meshgrid(grid_1d, grid_1d)\n", + "grid_z = np.zeros_like(grid_x)\n", + "points = np.vstack([grid_x.ravel(), grid_y.ravel(), grid_z.ravel()]).T" + ] + }, + { + "cell_type": "markdown", + "id": "7afed28d-1cb0-469d-ae57-2560f56dd344", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "### Electron density\n", + "\n", + "The electron density is the most fundamental property that can be obtained from the 1-RDM, and it can be evaluated at an array of points using the [`evaluate_density()`](http://gbasis.qcdevs.org/_autosummary/gbasis.evals.html?highlight=evaluate_density#gbasis.evals.density.evaluate_density) function. In this section, we calculate the electron density as an example of a one-index evaluation and display it on the molecular plane. For this example, screening is disabled with `screen_basis = False`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "ad3deaec-617b-4469-a6e7-862415d6955f", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAJjCAYAAADkjmzrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcHpJREFUeJzt3QeYVPX1//GztKWD9N6VYgGpgsYSEYgmkdiwRIqI0QCKGBVQWcCCRoJYCIg/URL1L2IXFUUUjQGBQEAgggFFFpAmsksTEPb/nO/mDndmp9yZuVPv+/U8k2Vm7ty5OzNmP3v2fM/NKSoqKhIAAADAw0ql+gAAAACAVCMUAwAAwPMIxQAAAPA8QjEAAAA8j1AMAAAAzyMUAwAAwPMIxQAAAPA8QjEAAAA8j1AMAAAAzyMUA4hJTk6OjBs3LtWHgTDOP/98cwEAREYoBuDz/PPPm7Ab6vLFF18k5TgOHjxoAvfChQsl3WzatMnvNSlbtqzUqlVLevToIWPGjJHNmzdLutq2bZt5XVeuXJnqQwGAtFMm1QcAIP1MmDBBmjdvXuL2Vq1aJS0Ujx8/3vw7XSud11xzjVx88cVy/Phx+fHHH2XZsmUyZcoUefzxx+XZZ5+Vq6++OtWHKB9++GGJUKyva7NmzaRDhw4pOy4ASEeEYgAl/OpXv5LOnTtLpjhw4IBUqlQpqc/ZsWNH+f3vf+9323fffSe9evWSAQMGSNu2baV9+/aSSuXKlUvp8wNAJqF9AoBrtm7dKjfccIPUrVtXcnNz5dRTT5WZM2eW2O6nn34yf8Y/5ZRTpHz58lK/fn257LLLZOPGjaY9oXbt2mY7rWpabQpW//LAgQOlcuXKZlut1FapUkWuu+46Xzi+4447pHHjxub5W7duLZMmTZKioiK/59f9DRs2TN5880057bTTfMc6b968uL7/pk2bmhaUI0eOyJ///Ge/+/bu3SsjRozwHZtW3R955BFTaQ5szdBjnjFjhrRs2dJs26VLF1OJttu+fbsMGjRIGjVqZLbR1/DSSy81+wjWU6ytKLofpY+zXlc93ry8PNMGsmvXrhLf00033STVq1c37xkAZDMqxQBKKCgokN27d/vdpgGqZs2aIR+zY8cOOeuss3yBU4Pt+++/L4MHD5bCwkITCNWxY8fk17/+tSxYsMC0GNx2222yb98+mT9/vqxZs0Z69uwp06ZNk1tuuUV+97vfmbCszjjjDN9z/fzzz9K7d28555xzTICsWLGiCb6//e1v5ZNPPjHPqe0BH3zwgdx5550mrD/22GN+x/v555/L66+/Ln/84x9NsH7iiSfk8ssvNz3B4b7PSLp3727CrH4/9naQ8847zxzHH/7wB2nSpIksWrRIRo8eLd9//71pu7B76aWXzGui2+rrqQFbX4dvvvnGhFelx7p27VoZPny4aYfYuXOneU49fr0eSCvX2hYzduxYE3R/8YtfmNu1F1pfR71v9uzZ5r2zaLh/9dVXzXPpLy8AkNWKAOB/nnvuOS2pBr3k5ub6bau35eXl+a4PHjy4qH79+kW7d+/22+7qq68uqlatWtHBgwfN9ZkzZ5rHTp48ucTzHz9+3HzdtWtXif1bBgwYYO4bNWqU3+1vvvmmuf2BBx7wu/2KK64oysnJKdqwYYPfsZcrV87vtlWrVpnbn3zyybCv0bfffmu2e/TRR0Nuc+mll5ptCgoKzPX777+/qFKlSkVff/2133b6PZQuXbpo8+bNfvuuWbNm0Z49e3zbvfXWW+b2d955x1z/8ccfIx6DOu+888zFsmzZMvM4fZ8Dde/evahbt25+t73++utm+08++STs8wBANqB9AkAJU6dONVVH+0WrvqFoznzttdfkN7/5jfm3Vpmti1Z0tfK8YsUKs61up9MatMIZSKuiTmkl2e69996T0qVLy6233up3u7ZT6DEFHr9WpLWia9FKdNWqVU01Nl7a3qG02qvmzJljKrMnnXSS32ujx6CV888++8zv8f369TPbWqyqrnVsFSpUMP3C2hKhi/zc0L9/f1myZIlpS7G8+OKLpt1Dq9wAkO1onwBQQteuXaNaaKe9qNozq32weglG/7yvNHRpr2+ZMrH/348+VntpAxe5NWjQwLRCBLYNWPfbaQtDIA2iboTM/fv3m6/Wsfz3v/+VL7/80tcrHeq1CXVsVkC2jk17iLUfWQO/9m9r24q2pGiwrVevXkzHrEFcW1w0CGuLhf4iM3fuXLn99tuj+mUFADIVoRhA3KzFYjqNQScvBGPvCY6XhsJSpeL7Q5dWlYMJXJQXC+2NrlOnjqk8W6/PRRddJHfddVfQ7XXBYbTHpgFWK/O6WFB7p++77z6ZOHGifPzxx3LmmWdGfcwavDVYW6FYe4kPHz5cYsIGAGQrQjGAuGkFVKui2gqgLQHhaMuC/pn+6NGjvkVjgWKpTOrkh48++si0LNirxevWrfPdnwyLFy821XB7mNTvWavHkV6baOl+tVqsF61G6+LCv/zlL/LCCy/E9LpqpVknWOikCw3HGq51KgcAeAE9xQDippVNnVCg/cJaJQ1kH/Wl22k/7VNPPRWyEqrTJJS2ZDil49k0lAfuV6dOaBjU2cuJpi0aOjJO+3116oXlqquuMmFZK7qB9HvUaRrR0GkWgSPSNCDrLwNa3Q3FmuUc6nXV10j7vbU149NPP6VKDMBTqBQDKEEXpVkVVjsd39WiRYugj3n44YfNOLRu3brJkCFDpF27drJnzx6zwE4ruPpvqxr5t7/9TUaOHClLly41i8h0vrBuo+PRtFKpC8n08ToiTFsLatSoYeYJ6yUUbSW44IIL5J577jGzevXEGXpGt7feesu0GtgX1blBvy+tyGprhIZMra7qLwUawP/+97/7tYtoQH777bdNe4KG5k6dOpnvefXq1aZNQY9Xw6hTX3/9tVx44YUmbOvrpD3Wb7zxhhmLF+5Mevoa6Mzh6dOnmwCtIVnfL+vshVq518frLxb6i46etQ8APCPV4y8AZMZItsBRXsFGpu3YsaNo6NChRY0bNy4qW7ZsUb169YouvPDCohkzZvhtp+PZ7rnnnqLmzZv7ttPRaRs3bvRts2jRoqJOnTqZ0Wn259KRbDreLJh9+/YV3X777UUNGjQw+z355JPN2DJr1Jv92PU4AzVt2tTsPxxrbJp1KVOmTFGNGjXMOLPRo0cXfffddyGPTe9v1aqV+Z5q1apV1KNHj6JJkyYVHTlyxG/fwUat2V8DHXunx9+mTRvzWujIO33+V155JexINmu8W7t27cxxBxvPtnTpUnN7r169wr4OAJBtcvR/Uh3MAQDpYdWqVaY3Wav5119/faoPBwCShp5iAIDPM888Y+YsW2cSBACvoKcYACDvvPOO/Oc//zFzpvVUz9aiPADwCtonAADSrFkzs1BPz0CoCwUDT4ICANmOUAwAAADPo6cYAAAAnkcoBgAAgOd5aqGdDtnftm2b6ZWL5TSyAAAAdtqFqqeXb9CggZQqRa0xk3kqFGsgbty4caoPAwAAZJn8/Hxp1KhRqg8DcfBUKLZWUz8+5jmpUL5iqg8HAABkuEM/HZTbHhrExJYs4KlQbLVMaCAmFAMAALfQlpn5aH4BAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4MjUqVOlWbNmUr58eenWrZssXbo07PZTpkyR1q1bS4UKFaRx48Zy++23y08//RTVPnX7oUOHSs2aNaVy5cpy+eWXy44dO/y22bx5s1xyySVSsWJFqVOnjtx5553y888/SzQIxQAAAIho9uzZMnLkSMnLy5MVK1ZI+/btpXfv3rJz586g27/00ksyatQos/1XX30lzz77rNnHmDFjotqnBul33nlH5syZI59++qls27ZNLrvsMt/9x44dM4H4yJEjsmjRIpk1a5Y8//zzMnbsWIlGTlFRUZF4RGFhoVSrVk1mTJgtFcpXTPXhAACADHfop4Ny09h+UlBQIFWrVpVs1q1bN+nSpYs89dRT5vrx48dN9Xf48OEm/AYaNmyYCcMLFizw3XbHHXfIkiVL5PPPP3e0T31da9eubQL2FVdcYbZZt26dtG3bVhYvXixnnXWWvP/++/LrX//ahOW6deuabaZPny5333237Nq1S8qVK+fo+6NSDAAAkIa0bUALeom6aOAMvO3w4cNBj0WrsMuXL5eePXv6bitVqpS5ruE0mB49epjHWO0Q33zzjbz33nty8cUXO96n3n/06FG/bdq0aSNNmjTxbaNfTz/9dF8gVlpt1u9n7dq1jl/vMo63BAAAQNICcdNmdWTnjn0Jew7tz92/f7/fbXl5eTJu3LgS2+7evdu0KdiDp9LrWrkN5tprrzWPO+ecc0QbE7TH9+abb/a1TzjZ5/bt202lt3r16iW20fusbYLtw7rPKUIxAABAmtEqqgbiL9feI1WqlHd9//v2/SRnnPqg5Ofn+7V95ObmuvYcCxculIceekj++te/mjaJDRs2yG233Sb333+/3HfffZJuCMUAAABpSgNxlaruh2KLBmInvdC1atWS0qVLl5j6oNfr1asX9DEafK+//nq58cYbzXVtcThw4IDcdNNNcs899zjap37VXxD27t3rVy0O3CZwYoW1z1DHFgw9xQAAAAhLWxg6derkt2hOF8Xp9e7duwd9zMGDB02PsJ2GYKXtFE72qfeXLVvWb5v169ebEWzWNvp19erVfhMr5s+fb8J+u3btxCkqxQAAAIhIR6cNGDBAOnfuLF27djUziLXyO2jQIHN///79pWHDhjJx4kRz/Te/+Y1MnjxZzjzzTF/7hFaP9XYrHEfap04NGzx4sNmuRo0aJujqZAoNwjp5QvXq1cuEX61K//nPfzZ9xPfee6+ZbRxNOwihGAAAABH169fPjDjT+b8aPDt06CDz5s3zLWrT6q29MqzBNCcnx3zdunWrGa2mgfjBBx90vE/12GOPmf3qSTt0OoZOltA+ZYsG7Llz58ott9xiwnKlSpVM0J4wYYJEgznFAAAAaTan2Mos326+PyE9xfsKf5LmTe7zxHxlp+gpBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOdlbCh++OGHJScnR0aMGJHqQwEAAECGy8hQvGzZMnn66afljDPOSPWhAAAAIAtkXCjev3+/XHfddfLMM8/ISSedlOrDAQAAQBbIuFA8dOhQueSSS6Rnz54Rtz18+LAUFhb6XQAAAIBAZSSDvPzyy7JixQrTPuHExIkTZfz48Qk/LgAAAGS2jKkU5+fny2233SYvvviilC9f3tFjRo8eLQUFBb6L7gMAAADI2Erx8uXLZefOndKxY0ffbceOHZPPPvtMnnrqKdMqUbp0ab/H5ObmmgsAAACQFaH4wgsvlNWrV/vdNmjQIGnTpo3cfffdJQIxAAAAkHWhuEqVKnLaaaf53VapUiWpWbNmidsBAACArOwpBgAAAMTrleJgFi5cmOpDAAAAQBagUgwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAABHpk6dKs2aNZPy5ctLt27dZOnSpSG3Pf/88yUnJ6fE5ZJLLvFt8/rrr0uvXr2kZs2a5r6VK1eW2M+MGTPMvqpWrWq22bt3b4ltHnzwQenRo4dUrFhRqlevLrEgFAMAACCi2bNny8iRIyUvL09WrFgh7du3l969e8vOnTuDbq+B9/vvv/dd1qxZI6VLl5Yrr7zSt82BAwfknHPOkUceeSTk8x48eFD69OkjY8aMCbnNkSNHzH5vueUW8eRpngEAAJAckydPliFDhsigQYPM9enTp8u7774rM2fOlFGjRpXYvkaNGn7XX375ZVPJtYfi66+/3nzdtGlTyOcdMWKE+bpw4cKQ24wfP958ff755yVWVIoBAAA8qrCw0O9y+PDhkJXY5cuXS8+ePX23lSpVylxfvHixo+d69tln5eqrr5ZKlSpJOqJSDAAAkKY2Hq0llY5UcH2/B44eMl8bN27sd3teXp6MGzeuxPa7d++WY8eOSd26df1u1+vr1q2L+Hzae6ztExqM0xWhGAAAwKPy8/PNAjZLbm6uJIKG4dNPP126du0q6YpQDAAA4FEaiO2hOJRatWqZRXI7duzwu12v16tXT8LRxXTaTzxhwgRJZ/QUAwAAIKxy5cpJp06dZMGCBb7bjh8/bq5379497GPnzJljepV///vfSzqjUgwAAICIdBzbgAEDpHPnzqYNYsqUKaYKbE2j6N+/vzRs2FAmTpxYonWib9++ZhZxoD179sjmzZtl27Zt5vr69evNV60+WxXo7du3m8uGDRvM9dWrV0uVKlWkSZMmvgkXug9rX9r7bM07btWqlVSuXDnyN0coBgAAgBP9+vWTXbt2ydixY01I7dChg8ybN8+3+E4DqU6ksNOQ+/nnn8uHH34YdJ9vv/22L1QrnU4RuOBPR79ZI9fUueeea74+99xzMnDgQPNvPaZZs2b5tjnzzDPN108++cSc+MOJnKKioiLxCB01Uq1aNZkxYbZUKF8x1YcDAAAy3KGfDspNY/tJQUGBo97caDPLRxunSaUqCZg+se+Q9Gx5i+vHncnoKQYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIAjU6dOlWbNmkn58uWlW7dusnTp0rDb7927V4YOHSr169eX3NxcOeWUU+S9997z3T9x4kTp0qWLVKlSRerUqSN9+/aV9evX++7fs2ePDB8+XFq3bi0VKlSQJk2ayK233ioFBQW+bVatWiXXXHONNG7c2GzTtm1befzxxyVaZaJ+BAAAADxn9uzZMnLkSJk+fboJxFOmTJHevXubEKuBNtCRI0fkoosuMve9+uqr0rBhQ/nuu++kevXqvm0+/fRTE5o1GP/8888yZswY6dWrl/znP/+RSpUqybZt28xl0qRJ0q5dO/P4m2++2dym+1TLly83z/HCCy+YYLxo0SK56aabpHTp0jJs2DDH319OUVFRkXhEYWGhVKtWTWZMmC0VyldM9eEAAIAMd+ing3LT2H6mclm1alXXM8tHG6dJpSoVxG0H9h2Sni1vieq4u3XrZsLrU089Za4fP37chFCt5I4aNarE9hqeH330UVm3bp2ULVvW0XPs2rXLBFwNy+eee27QbebMmSO///3v5cCBA1KmTPD6rgbtr776Sj7++GNxivYJAAAAj9Lwbb8cPnw46HZa9dWKbM+ePX23lSpVylxfvHhx0Me8/fbb0r17dxNQ69atK6eddpo89NBDcuzYsZDHY7VF1KhRI+w2GuRDBWJrm3D7CIb2CQAAgDS16WBNqVja/b9uHzx40HzVSq9dXl6ejBs3rsT2u3fvNmFWw62dXtdKcDDffPONqdRed911po94w4YN8sc//lGOHj1qnieQVp5HjBghZ599tgnQwehx3H///aY9IhRtn9BWj3fffVeiQSgGAADwqPz8fL/2idzcXNf2rSFXWyFmzJhh+ns7deokW7duNS0VwUKxVpTXrFkjn3/+edD9aSX7kksuMb3FwYK70sdfeumlZv/amxwNQjEAAIBHaSB20lNcq1YtE2x37Njhd7ter1evXtDH6MQJ7SXWx1l0MsT27dtNO0a5cuV8t+uCuLlz58pnn30mjRo1KrGvffv2SZ8+fcyUijfeeCNoj7IuzrvwwgtNFfnee++VaNFTDAAAgLA0wGqld8GCBX6VYL2ufcPBaBuEtkzodpavv/7ahGUrEOu8Bw3EGnS11aJ58+ZBK8Ra9dXHaJ+yjoMLtHbtWrngggtkwIAB8uCDD0osCMUAAACISMexPfPMMzJr1iwz2eGWW24xEyAGDRpk7u/fv7+MHj3at73er3OGb7vtNhOGtcdXF9ppm4RF/62j1F566SVTBdYqsl4OHTrkF4j1eZ599llz3drGWrCnLRMaiHU7PUbrfp1kEQ3aJwAAABBRv379TNAcO3asCZ0dOnSQefPm+Rbfbd682UyksOgivg8++EBuv/12OeOMM8ycYg3Id999t2+badOmma/nn3++33M999xzMnDgQFmxYoUsWbLE3NaqVSu/bb799ltzIhGdV6zHpeFaL5amTZvKpk2bxCnmFAMAAKTpnOL/W/2KVKySgOkT+w7Kjadf5fpxZzIqxQAiWvnV9zE9rkPb+q4fCwAAiUAoBhAxBJ/R45So9/Hloq9L7IeQDABIV4RiAK4F4XCPDwzJBGQAQDohFAMeZg+p8YbgSAL3v3LR175/E5ABAKlGKAY8HIYTHYTDsT+3FZAJxwCAVCEUAx6SzMpwNKxjoXoMAEgVQjHgEelQHY7EOjZ7/zHhGACQDIRiIMtlQhgORDgGACQboRjIUunaKhENwjEAIFkIxUAWysTqcDThmGAMAHDbiRNUA8h4Ghj1oiEyWwJx0AV5MZ5hDwCAUKgUA1ki26rDTidVUDUGALiBSjGQBbwSiO2oGgMA3EQoBjKcFwOxhWAMAHAL7RNABvNyILbQTgEAcAOhGMhQyQrES7/dGtPjujZvKMmkrwPTKQAAsSIUAxkoGYHYHobbd2wW1WNXrdjke3wywzHBGAAQK0IxkKESEYjjCcJ29scuXbHJ9+9kBGR7OwXBGADgFKEYyDCJWFTmVhgOxtpfsqvHGo7pMwYAOMX0CcDDVWINqVZQ1fDqdiC2s+/f/ryJxHQKAIBThGIgg1hnq3NDssJwqsMxwRgAkFWheOLEidKlSxepUqWK1KlTR/r27Svr169P9WEBSeNmqLMH4lQJDMeJRDAGAGRNKP70009l6NCh8sUXX8j8+fPl6NGj0qtXLzlw4ECqDw1IGjeqxOkQiO0IxgCAdJAxC+3mzZvnd/355583FePly5fLueeeG/Qxhw8fNhdLYWFhwo8TSAS3glysgfiLLVtK3HZWo0auHJN1PNZCvEQuwmNkGwAg40NxoIKCAvO1Ro0aYVsuxo8fn8SjAtK3ShxLILaH4Y5tG/vf91V+0MfEGpZ9FeP/jXBLVDgmGAMAMrp9wu748eMyYsQIOfvss+W0004Lud3o0aNNeLYu+fnBf4gD2V4ljicQaxgODMT22+0X63HWJV3bKWilAABkRSjW3uI1a9bIyy+/HHa73NxcqVq1qt8F8FqVONpAbA+0wcJwOMECciwIxgCAZMu4UDxs2DCZO3eufPLJJ9LIxZ5GIBtHsMUSiMNVh6MRWDmOFsEYAJBMGROKi4qKTCB+44035OOPP5bmzZun+pCAjBBLIHZLvOGaYAwASJZSmdQy8cILL8hLL71kZhVv377dXA4dOpTqQwMSIhUhzQqwDVatkKG9e/gubT58t8S2VbZ/77dNl7//X9j90koBAEhnGROKp02bZhbLnX/++VK/fn3fZfbs2ak+NCBh4m2diKZK7GaFONzzxCIZM5UJxgDgbRnVPhHsMnDgwFQfGpDRnATVDft+kEW7NvsusXBj8R0n+AAAiNdDMeAlyQ5lTqrEXVo2NhelwXj5nq1JD8aKYAwASARCMeDh1gmn4bR5nRpBw3Es0n3hnSIYA4D3ZOwZ7QDEJ5ppE02/WCLl9xafRdLSY98+v+v5B/zvD8csvPsqP6az33FKaABAIhCKAQ9zWrVt/dECc4nE6jfuUbuJ42BOMAYApAPaJ4AsPGFHvNMa1u7dEfVjGtao5tdzHIlbZ71LNFopAMAbCMWAB4UbwRYs0H547xiZsviffpeZr78a9PHJDsaJ7i9WBGMAyH6EYgAltGlYJ67HR7MQz435yARjAEC8CMUASli3dWfSnzOdz3hnIRgDQPYiFAPw43SRnJvirRYnq79YEYwBIDsRigGUcGr1unE9ftnGfEm2ZPUXK4IxAGQfQjHgUSu+yk+7inM8Z7qzEIwBALFgTjGQRePYnNLZwOEC6Lb2HWXqB4vMBIlQi+YK69c3UyjcbKGIN6hb84uTxReMF33NHGMAyHBUigG41gqh26aidSJVbRT2cKy/0FA1BpDtpk6dKs2aNZPy5ctLt27dZOnSpSG3ff755yUnJ8fvoo+zGzdunLRp00YqVaokJ510kvTs2VOWLFnit81vf/tbadKkiXls/fr15frrr5dt27b57l+/fr1ccMEFUrduXbNNixYt5N5775WjR49G9b1RKQYQtgVCq8UadENVjANDcLwL9WI9y12gRJ/tLhBnv0svbv6CwvsJFJs9e7aMHDlSpk+fbgLxlClTpHfv3iaU1qkTfJRn1apVzf0WDcZ2p5xyijz11FMmyB46dEgee+wx6dWrl2zYsEFq165tttHAO2bMGBOIt27dKn/605/kiiuukEWLFpn7y5YtK/3795eOHTtK9erVZdWqVTJkyBA5fvy4PPTQQ+IUoRiAo2CcyCDsZgtFMk8DHYhgnF7B161WJG2PCYX3GV4yefJkEzYHDRpkrms4fvfdd2XmzJkyatSooI/REFyvXr2Q+7z22mtLPMezzz4rX375pVx44YXmtttvv913f9OmTc1z9e3b11SCNRBroNaLfZuFCxfKP/7xj6i+P0Ix4FGmr/irfMfj0BIVhBOFYOyNEJyMHvxwzxEYmHnPkWkKCwv9rufm5ppLoCNHjsjy5ctl9OjRvttKlSpl2h0WL14ccv/79+83IVWrtlrJ1crtqaeeGnRbfY4ZM2ZItWrVpH379kG32bNnj7z44ovSo0cPE4iD0SrzvHnz5LLLLpNoEIoBRJTsAOxWCwXBOLOlKgRHw3481ntu4b2HG344UEkO5FRyfb8/HShuY2jc2L8wkpeXZ/p8A+3evVuOHTtm+nbt9Pq6deuCPkfr1q1NFfmMM86QgoICmTRpkgmza9eulUa2/4+fO3euXH311XLw4EHTIjF//nypVauW377uvvtu02ah25x11lnmMYF03ytWrJDDhw/LTTfdJBMmTIjqNSEUA0grbrVQpEMwtlcSCUixtUKkWwiONSArPgNIR/n5+abv1xKsShyr7t27m4s9tLZt21aefvppuf/++323a8/wypUrTfB+5pln5KqrrjKL7ex9ynfeeacMHjxYvvvuOxk/frzpIdZgbO9R1p7nffv2mZ5i3V5D+F133eX4eAnFgIdF20KRqZI9qi1U1VgRjE4IDI2ZFIAjCfxeqCIjXWkgtofiULRyW7p0admxY4ff7Xo9XM+wnbY7nHnmmaa9wU4nT7Rq1cpctAp88sknm75ie6uGPr9edGGeBmutcH/xxRd+oduqerdr185UtbVafMcdd5jjdoKRbECa0R+e8dBqaLQBMNEn8oiFGyfySPWoNntA4mQfxd+7/WK9LvbXJ1sFfp/21wHIBOXKlZNOnTrJggULfLdpn7BetwfTcDSorl692rRIhKP71RaIcPerSNvoQjxrWyeoFANpRKtHyf4haZ3IY0UaVYzdbqHwC8YrNiW1jcLLVeNM6AlOpzaLbP88IPONHDlSBgwYIJ07d5auXbuakWwHDhzwTaPQloaGDRvKxIkTzXXt6dXKr1aA9+7dK48++qhpf7jxxhvN/frYBx980Mwh1qCs7RM6B1nHrl155ZVmG22jWLZsmZxzzjlmjvHGjRvlvvvuk5YtW/rCuC680yr06aefbto//vWvf5kqc79+/UIuxguGUAwg4hnuskk6BONsDkPZ3BKRCPbXxz7JIps+E8ge/fr1k127dsnYsWNl+/bt0qFDBzPlwVp8t3nzZjORwvLjjz+aEW66rQZarTTrbGFtb1Da1qCL9GbNmmUCcc2aNaVLly5mlJo1oaJixYry+uuvmwWAGqI1PPfp08ecnMPqfy5Tpow88sgj8vXXX0tRUZGZdjFs2DC/UW5O5BTpoz00dkTHfMyYMFsqlK+Y6sMBEnaqZ6tNQANgNDQYp0u12KoUuzGFIphVKQzGwdplMjkEEYTdlQ2fCS859NNBuWlsPzNdwUlvbrSZ5ZEv3pXylRMwfWL/Abn7rEtcP+5MRqUYyEIa9mLtn02XNopEtVDYJXsihZMpFZkShgjCicPkEiA1WGgHwMeqyqbjwju3WVX0VC2+CxRsEVa6CVwc5pVFcqnC4kwguagUA1lMWwSibaFIt/5it07kkU4zjDOlx5RqcOpRNQaSh0oxkKXiCXkaQtOhWpyMNo50qxinunpMNTg9pftfEYBsQKUYSNOFNukQQtKlvzhpwThNFt85rR7bRVNBDBeq0uFzB2eTS6gaA+4iFANZPqs4lhaKdJtfnMgWimDj2lQ6huNwoTVUWI5mH5kumdX+VH4+vDbvGkgWQjGQxeKZQpEu/cXWFIpkBuN06zP2ctANJtxnOpZfAGM6hiBnjUzm54WqMeA+QjHgAbFWi33BOMXV4mSMZ8uGYOylEJys8BtK4PNbnxdLsj439qoxwRiID6EYyPK+4nirxZZUt1Hoc2s4T0a1WBGMUyvwM5vqEBxJ4PEFVpIT+RkiGAPuIBQDHugrdqVanAb9xakIxird+4yzQaaF4FhDcqI+QwRjIH6MZAM8wI0fxBpErVFtqRzXZoLxli1J7XVO57FtmU5fU/tpya1LtrF/hhL1OeJkH0B8CMVAGtPKj9vV4mw4651VrSYYZ1cYznb27zNR4diaKc08YyB6tE8AHmmhcKu3OF3aKZI9lULRThE/+2cwWUE4nl+eEvHZsr7vRPas004BRI9QDHhMPL3Fdl4NxopFeOkbhIMF4Fg/n9ZnK5Bbn7VE/5JFMAaiQygGPMTNarE9HOjiN5WKcJwOwVgRjlMTht0MwU73Y33eLfF+7hIZjgnGgHP0FAMe6yt2q7c4nfqMU9FjHKxHFMnpF7YWWlrvt77/9kuiBT5X4PHEKlH9xizAA5yhUgx4bDSb29XidDn7nb1ibB1PslA1TmxVOPBzlerTjoc6lsB2i1g+g4nqN6ZiDERGKAY8yq3e4kCp7jG2jiHZ4TgZi6e8FIbTOQiHEniM8bRZJKJvnWAMhEcoBjKhWuzS2e2yvVocKhxTNfZmEF60a7O4pUftJq5UkdMlGAMoiVAMeFg2VovToaUisGqcDeE4kWecs4fheD43wUJwl5bxfw6XbcwPGbCdhuVYP4uJWIRn5hgv+ppqMRCAUAx4VLZXi9OlamyxQk0mBeREn3o5njCcqAAcTKj92sOyk3Acz2cxEVVj2igAf4RiwKMtFF6pFqfDQrxMqx4nepSaW2E4USHYKfvzL9p44rgiBWTf1Ir/9RxHUzV2KxjTXwyURCgGPMwr1eJ0qBqnezhO9kzhVIbhpQXfRNyma7UWUe3TOqZoqsex/KLmZjsFwRjwRygG4JlqcbpUjUOFY0uiQ3KwX4QSeaa5VIThSMH3rEbh9/XFlm9iCsvRVo9j/UXNraoxwRg4gVAMZJAvE9BC4bVqcTpVjYOF0WAhWcUTfJIdguMNw24F4UjBN5xgj/1iS36J5wgXkqOpHsfSUuF2MAa8jlAMePhEHl6uFqdT1dgu2HsQKijHs89sCcNuBuFIAvdtD8nRhONwbRXRnrbczWDMRAp4HaEYgGerxelWNU6XUJvuYdgehN0Iwd8cWut3vUWFUx09zv7c9naLUAFZvzenVeNUBGNFGwW8jFAMZJhEtFAkulqcKWKdCuB1bvQMJyIMB4bdULrVbup3fcmu4I8LF5at44lUPXZaNY41GMeD/mJ4HaEYyCCJbKFIVLU4G8KxIiCHD8SJCsPRBOHAEBwYdp0K9bhgYTkwKAcLx8ECsr1q7GYw1qkU9BcDsSEUA4DDtgpFOE6vMOxWEFb7jy42XyuX7e6wovyd7/lDhWN7QA4WjO1TKoKF42iDsaKNAogNoRjIQJnUQmH6itN8sV0khOPk9Q07CcP2IBwpBFtB14nTq50sqwv+G/QxwYKy/bmtKnKwFgv9PsK1VUSqGkcTjN3oL6ZaDK8iFAMZhhaK1En3BXmJlElhODDUath1Kti2ToKydTz2Fgt7QA5sqwgXjFVgOI4lGMeLajG8hlAMAFHy0oK8WNskYmmViCcMxxOE4w3KwcNx8NYKe9U4VDtFqKpxtK0UblSLCcbwEkIxkKnV4gxqochW6TbjONvCsLKCZbAw7DQIV/j5I4nHoTI9S9xmPVcs4ThcMI7UTuE0GNNGAUSPUAzADy0U3u43jicMJ6I6HCkMOwnCjSq2l1htObiqRKi2h+Row3GwYBxtn3G0wTheVIvhFYRiAPB4OI5nznAyWyUiheHA8BpPGA63jy0HPyoRkJ2E42DB2GmfcTzBOJ4xbVSL4SWlUn0AAGKjlZtE/bDSH6BuVJjs9Ae4F2hY8fUcb9mStmf0s47NXhmOtVXCjUCsYdFeHbYHYg2ZVtDU8BkYiDUMW4FYQ6x1SRT7c9ifO/D47MdtfV/W92pnvR6Bp6z2G9tmW7BosX/OIonnrz/mFNAJPMU8kC6oFANIuEw53XOiK8eprh4HvgfxjslLZqtEsCBs5yQE16z2J4nWDwWTwt6vz2tvsQhVObaqxvYpFdEswNOKcTBWxThZp4EGshmhGEDWCVZVc0uos4+FEhg87WfJS0ZIdjsIJysQh2uVsFeFExGEwz0+WEi2jiNUOA4MxifaKZwHYxXu7HeJ7i821eJFX9NbjKxGKAYyGFMoYg9rMe//f2cfiyUgB4bSwCpyPEE5XCXezROnOJ077EYgDlcdDheI4w3C4dj3HRiQA8Oxm8HYSX+xE1SLgdAIxQCyZgpFogOxfd/hTrTgVKiwGlhNjmdfqXp9I/UPx9oukaowHO75goXjwKpxLME4UKQ2Cv3cJHoaBZMokM0IxQCyQqyBONjiJkuwP2Fb7M8Tb/U4ULqdEjuV1eFIgTjZQTgYPQYnVePog7H7bRTxVIuZRIFsRygGMlyiWyiyNRA7Oa2whpJgIp2JzK1wnA7cqg7H0y6RrmHYadXYGuMWLhgr+7g2i9ttFG7NLgayESPZACR1NFuqA7GGDHt4C3c2Nev+wIt9H3Z6DNbFPqosE0Uzas3rgdjJsVnj25R9ZFu4cW2hRrVFej8SPaKN8WzeNnXqVGnWrJmUL19eunXrJkuXLnX0uJdffllycnKkb9++Ibe5+eabzTZTpkwJev/hw4elQ4cOZpuVK1cG3WbDhg1SpUoVqV69ukSLUAwgY8USiJ2E4UjsYSVU+4V9vmymhWP76+pG/7BeAmcPZ2Mgth9jsON0MxjHM7s4nsWzXv6LFERmz54tI0eOlLy8PFmxYoW0b99eevfuLTt37gz7uE2bNsmf/vQn+cUvfhFymzfeeEO++OILadCgQcht7rrrrrD3Hz16VK655pqwzxMOoRjIEl7r9Yu2iqmXwDBsBbZoL8q+r3DBONzJF9KNFeCjCcPBXlenC+oCT8RhnQgj2Mk3QgXNdOZ2MI72pB7hmDPdUS1GlCZPnixDhgyRQYMGSbt27WT69OlSsWJFmTlzZsjHHDt2TK677joZP368tGgRfJ3G1q1bZfjw4fLiiy9K2bJlg27z/vvvy4cffiiTJoWeHX7vvfdKmzZt5Kqrrorhu6OnGMievuIE/pByazSb9jsm64QR4fqGw51KOJLQp+otfq5gC6N8AWZjevYbR7OQzkK7RHy9xuEW34U6JXTgwrtw0ygSueiOBXfZpbCw0O96bm6uuQQ6cuSILF++XEaPHu27rVSpUtKzZ09ZvPjEL3aBJkyYIHXq1JHBgwfLP/7xjxL3Hz9+XK6//nq588475dRT/XvqLTt27DBh/M033zQhPJiPP/5Y5syZY9oqXn/9dYkFoRhAUkazuXFWu2hPGOE0DNsrdZF0q10cXHRagPIPLcXzZcOF43RajBdLGFYE4vgmVNhHtkUKxoGCndiDRXfZ7fj+snK8qKz7+z1QvM/Gjf3/O87Ly5Nx48aV2H737t2m6lu3bl2/2/X6unXrgj7H559/Ls8++2zI/l/1yCOPSJkyZeTWW28Nen9RUZEMHDjQ9Bt37tzZtGIE+uGHH8w2L7zwglStWlViRSgGsohWb7K95y+akWD2dgeLhg0NH/uPbpNmZVaY29pVqun4+XceLn7MqdU7mq9r9564r87/drPzh6ohw3HgpIpUBONYw7AiECcvGKuS1WL/oOtGtRjelp+f7xckc4NUiWOxb98+UwF+5plnpFatWkG30crz448/bvqTdfFcME8++aTZl71CHUiryNdee62ce+65cR0zoRjIEoluoUg1q9fVSWCrU7PQL5zZK8G1iuZKrTIidXI1wdY04SSa8KWhRsNMndziKp4G600/FwdktXZvA/P8akthgSwtSJ+qcTxhOJ5AbPFyIA7GCsYWKxhbklUtXrpiU0wtFJz2OTtoIHZSXa1Vq5aULl3atDLY6fV69eqV2H7jxo2mqvub3/zGr1VCaWV4/fr1pp1CF+k1aXLis6vV6DvuuMNMoNDHa1uEtmcEhnWtGmuv8qxZs8w2b7/9tq/fWKvL+lz6PDNmzJAbbrjB0WtBKAaQ9pwsUrMHYnVq9W3mq1aEVddKVrgoDsKxBq/iBV/F/y6u+p0IyDsP/yDiNwVIV0kXytIfnFWNExmO4z3bXzyBOFh7ilcDcbATfdhPCa3crhYDbihXrpx06tRJFixY4BurpsFTrw8bNqzE9rrgbfXq1SUWwmnVV6vD2rahlWTtSbbTaRZ6uy7mU0888YQ88MADvvu3bdtmttFJGDoSTmlo1jBteeutt0xbxqJFi6RhQ+e/8BGKASR1sV2snIS5wEB8IggXiycMhwvIJ0LOiYBs2Xpwr0jNdr6WinAL8RIRjt049bUbgTjYmeq8FojdrBYHE6odhxYKuGXkyJEyYMAAU6Xt2rWrqeYeOHDAF2D79+9vQujEiRPNHOPTTjvN7/HW7GDr9po1a5qLnU6f0Mpz69atzXV7FVlVrlzZfG3ZsqU0+t/num3btn7b/Otf/zKLAAOfPxJCMZBlEtFX7NZiu0RWie2BWFsaitsj/INXogLXif2eqABq0NFQvlQ0KP4nYjAOFY5jCciBr1kiA7ElmkBsVYm9GogDq8XFwTjaarGzSRROWijiOfUzvKVfv36ya9cuGTt2rGzfvt2cSGPevHm+xXebN282YTRTEYqBLJKtfcXhQl2oQJyMMBzI/3kmxRSMA79fJ9VjN0NwtIHYOjFHoFgW1nlJuDaKeKrFsYhnCoX+/w19xd4ybNiwoO0SauHChWEf+/zzz0fcf7DpEnZ6Nj3tGQ5HJ1HoJVqEYgBpK1KVOFIgTmXlsfi5T7RV2IOxVvzCzTV2Uj0Oto2bnAbicFhYF3sbRaRqsQr85cqcXnxj8loomFeMbJO5NW4AnhAq8KVzILZYx6DHpBXjhhWrS/Vy/zHhJtLZ8AJZZ5kLdnFbuLPURXO2umAIxOG/78Cz3VkCX+doTlMe7wlzAK8gFAPIuBFsmRCIQwXjs6t+E3MwTgYrEIcT7aQJZZ2+OVA6vVepYP/+g70+ga9nYHU+2GcnllOKx3vaZyAbEIqBLJSoP2km6+xX4X6oZ1IgDhaM9Vg1GOuxBwbjVIZj6/ndCMQsrItPPNXiRPzlAPAKQjGQZRK16CXZK9OD/XC3h0Zr7Fq6B+JQwVjDvD0Yp7Jq7HTCRDyB2OsL62KtFuvrGc1pyMO1UMR7mvWQJ/HIwsW98CZCMYCMolViK5RpO0ImBOJQwTgwbKYyGDvtUQ11tjpFH3Hi6esfqYXCLLiLoYUC8DpCMYC0EuqHuf0Hf2DlLJMCVuCxxrOAyg1OWiYiTZpwo5IJ/18iQrWcBEr25wXIZoRiAGknVF+kPQAEnq0uE+n3ECxQ6veZjGpxtM8RbZWYxXXORXpNkvGLB4vt4HWEYgBpL1R408CViQFLjzkwLAarxCYyGDvtI3YyjxiJFaotJV36ioFsQSgGkBE0vGVrOLMqsPbvL5H9xdEE4sBjjGYEG+LDawgkF6EYQFr3EweGQp3YYC2wy8QqcagWinDBOJWB2MkvIiywc1+kaR2BnxO35hXHggkUyBaEYgBp30+cjYuJgrVQBAvGiagWR/t6UiVODqe/MITr7bYwrxiIHqEYyFKJOoFHMgULgzrfN5srjqGmUbgRjJ1OmrBQJU4fbs0rBhAaoRjIQok8gUeyzmpnsYc4PT1ywk8EsWOHyP33i5x3nkjduiLlyolUqiRy6qkigweLvP++SFGRK09lhcTAKRSBs2jdqJTHGqqdVCUzqkq8apXIH/8ocvrpItWrF7+/+j7/8pcikyaJFBRIOnH79WSxHRAaoRhAWgjsfwwMcRoSG1asbk56kbCK41//KtKsmcjYsSKffSayc6fI0aMiBw+K/Oc/IjNnilx8sch37o2DswJ+sEqgk77RRC2s0+cOF4j1WKOpEqfczz+L3HabSIcOItOmiaxZUxyA9f3V9/mTT0TuvFOkVSuRDz9M9dHG9Dqm8jThQDYgFANIG4F9kIEhTlsnEha4/vxnkaFDRX76qfh66dIiv/2tyPjxIhMmiPTvX1xRdFlgwLeCcbQV2khS3Zed8taJ4cNFnnjixPUGDURGjBB54AGRfv2K32+1e7fIb34j8s9/ZtTr5Nb7y6xieFmZVB8AAKScVoHHjDlxvU4dkXnzRM480387rSo+/7xIxYoJOQytvK4u+G9C9u1pixaJTJ9+4nrHjsWV4apVT9z28cciF10kcvy4yJEjIn/4g8iXX4qUonYEeEXG/dc+depUadasmZQvX166desmS5cuTfUhAch0Tz4pcuzYiev65/XAQKzKlhUZMqQ4NCNzzJhR8q8C9kCstKdYK8aWtWtF/vGP5BwfgMwMxQMGDJDPtNcuBWbPni0jR46UvLw8WbFihbRv31569+4tO7UfDABitWDBiX+fdJJI376pPBq4zR5u9f298MLg29lDceDjAGS9qENxQUGB9OzZU04++WR56KGHZOvW5PUeTZ48WYYMGSKDBg2Sdu3ayfTp06VixYoyUxe/AECs7P8/dsop/Mk823xvO7FE0zC92oH32R8HIOtF/f/8b775pgnCt9xyi6ncaivDr371K3n11VflqPbbJciRI0dk+fLlJpBbSpUqZa4vXhx8duPhw4elsLDQ7wIAAAAEiqkcUrt2bdPGsGrVKlmyZIm0atVKrr/+emnQoIHcfvvt8t//ur9QZPfu3XLs2DGpG7D6W69v37496GMmTpwo1apV810aN+YMPwCCaNjwxL+//tq1OcRIE/Vtc7s3hzn1ceCoPfvjAGS9uP5G+P3338v8+fPNpXTp0nLxxRfL6tWrTWvDY489Jqk2evRo0+5hXfLz81N9SADSkb3H9McfRd56K5VHA7f94hcn/r1nT/GkiWBeeSX04wBkvahDsbZIvPbaa/LrX/9amjZtKnPmzJERI0bItm3bZNasWfLRRx/JK6+8IhN0rqeLatWqZYL3Dj3blI1er1evXtDH5ObmStWqVf0uAFDCsGEn5tSqW24pPvNZIG0R+7//Kz7ZAzLHTTf5X7/7bpF9+/xvW7hQV3OfuN6uHaEY8Jio5xTXr19fjh8/Ltdcc40Zh9ZBzw4U4IILLpDqevpMF5UrV046deokCxYskL7/Wxmux6HXh+kPNACIlZ7CWU/tbM0q1paszp1Ffv3r4tFsOTkiGzaIfPBB8WmgbWsbkAF69CieO/z008XX//UvkbZtRa66SisuIqtXi7z66omxfHrqZx3jxoJLwFOiDsXaFnHllVeaOcGhaCD+9ttvxW3ax6wj4Tp37ixdu3aVKVOmyIEDB8w0CgCIy+jRIpUqidx1l67SLT4t8JtvFl+Q+Z56qnjOtH61Jo4Ea/OrWVPkpZdEzj476YcIIMNCsS6oS5V+/frJrl27ZOzYsWZxnVap582bV2LxHYDsteXgKqlZLUE7v/VWkSuvLK4SfvSRyPr1xT3GWjnUcV3nn188yzbcWK8o/VAwyffvbD6bnX6fKT3Vc5kyxSdpGTy4+P399FMRXWdy6FDx7GL9a8HFFxe3WlRL1Acsus9DKqxasSmlzw+kUsad5llbJWiXALJPj9pNZNHGzdKlZfGUmK7VWsgXW76RsxqdmBrzz8IW8rvaBYkNWDpxIC+v+JJEViCuXLa7+bpkl/8khC+2xL5QONhr6YQeQ7faTcMes56aOtgvLo0qtpe0pC1/f/2rZAJ9HZ3Sz4e+z27o2tw2jQXwEBqmgCy08qvEnHRg6bfJO1lPoBYVTpW9R9rJzsM/RBUW0p1+L0sPNA0aiPV7tgfieENPNMHaeu5QrGMNdKgM/dZucvv1XPEVU5iAUAjFQJY6o8cpCdlv+47NJFmKK5z+P8Q3/dwxLf7M7Ab79xAqEFviDcTW46OtOAdWrAOFavkI9otLNrxnqeJma81ZjRq5ti8gmxCKAaS9eFoH0llgcAwWiN383qMNxlSLkyvSLw2hXm8A7iAUA0irvuJlG/PDVki3Htzra6HI5MqjdexW60S4QOxWr2iswThStRjuC9cipO9HpF+c9L8j/e8pGb5c9LV0aMvZ/5D5CMUAMoL+4Lf6iq0WimxhVQCDVWbdDMSJ2mewP+1rtZgWivhEU3FPxOcE8BpCMYC0F+oHfqZWi/WY7YExsPLn9jQBp/3asVSL+ZN+5vQTs8gOCI9QDCDtBLZQWOwhzmo7yGShvodk9lBH81zh2iiiWXCH6Ln1y0e4RXbMKIbXEYoBpJVQfZD2qmlgQMikanHgsa7d28DveiL6iN3oLw636C7aBXeZ9H4li/01SeUvEswohpcRigFkHHvFMhMrkYHH7Pb4tUQuvItl0V0mvkepFuoXimAndQn8vIT6SwuA8AjFANJSqB/sO3+o6ptCYbUfZEpvsb2XWI/dXiXWcJPoPuJ4g7EV3oMFY60Wh1pwF0wmvF/JEu61CPaaRhqVF+wvLonqJ9bJE0C2IBQDWXg2u0ScuEPPZpesE3dEaqHQYKxTKDIpGFuBWMfJ2QOxBpxktkykKhgzicIZfZ0Cf5Fw69TfiTppB+PYkC0IxQAyqlqcicHYHoh1nJwViPV7sKrD6TBSK9pgHApnuXOnl9hJlThY60Qs84lZZAcQigGkqXA/2DMpGIcLxCqaMKyBJ/CSymAcqlocLNBZ1U+CcXhOq8TRcNo6wSI7eB2hGEBGneEummCcyrBlPb8bgdgegPU1sS72+9wMydEsvoslGCP+KnG0C+wS0TpBPzGyDaEYQNqLNRirZIdj+/MFBmINNU4DsT3oBgZhi/12t0Oyk2Acqb84VDCmWhz8+w1VJY5Gslsn6CdGNimT6gMAkBmL7FJFf8gv2rU5bHhbWvBNcdis2U5E/iNLpal0rfRdQPgqDiE1q/0pIcdpD8IWDcT/LGwhe480MMe3UyJPl7AH2WgDTuD2+rpZ++vSsrFEy3ptNRif1ahxyGD8zaG1Jhh3q920RKjbf3Rx0Mfp69SoYvsSr2Gi3p90FqxKbA/EgWc8DPaLSrx/JaB1AqBSDMChZE2eCCXcD/1QFWPrYgUPt9sqrH1ZbRL2cWt6KQ7E7RxVh4O1SMQrWAU5WtYiwHgqxvQXR18ljiTYZynYZ0b7iRM1dQLeNHXqVGnWrJmUL19eunXrJkuXLg257euvvy6dO3eW6tWrS6VKlaRDhw7y97//3W+bgQMHSk5Ojt+lT58+ftvo8wVu8/DDD/ttU1RUJJMmTZJTTjlFcnNzpWHDhvLggw9G9b1RKQaQ9iJVi4NVjPceEVO51EqlFYyblVmhfyw2/95y8HpTqYylMhlYFbZaJOxBsDggFleIreMLxh6EE8Xat1U5jrVq/MWWbyJWjEPRYHx6tZP9QmCFnz/yXMU4MBBHWyUGUmn27NkycuRImT59ugnEU6ZMkd69e8v69eulTp06JbavUaOG3HPPPdKmTRspV66czJ07VwYNGmS21cdZNAQ/99xzvusaagNNmDBBhgwZ4rtepUoVv/tvu+02+fDDD00wPv3002XPnj3mEg1CMYCMESnQ2YNxnZqFvmBqqVy7u2w6sNgXkHce/ljqBAll4Vh9wsoKwlsPVjMtEid8F3H+cDxtEnGH443Fv2BEG46tYKyChePiiRRrQ7ZRhArGXmdViUONsrMLtcAu1s9QrP3ELLLzpsmTJ5tgqsFWaTh+9913ZebMmTJq1KgS259//vklguusWbPk888/9wvFGoLr1asX9rk1BIfa5quvvpJp06bJmjVrpHXr1ua25s2bR/390T4BZIlE9ROnC3sbQDj2VgoNadZFaUjWBW962Z0z1IRarSKv+PFjxxfdXh/37u5f+/alLRLW85je4Qjzh91uk4hWPC0VThbgsfDOeZU43OK6eMawJWMUG4vsskNhYaHf5fDhw0G3O3LkiCxfvlx69jzxmS1VqpS5vnhx8LUDge0NCxYsMFXlc8891+++hQsXmuqxBtpbbrlFfvihuPBgp+0SNWvWlDPPPFMeffRR+fnnn333vfPOO9KiRQtTidYwrO0WN954I5ViANmzyC5UG4XTirEV3LSqaf8T9ImFYcUBZLc4X+VfueyJsBJq8VOkVgnre0mlwJaKaCrH4RbgOVl4F1gxVqHaKFQ2tFI4aZuwC/UZi2YMmwrXT8wJOzJDmcIcKXssx/X9/nygeJ+NG/v/N5yXlyfjxo0rsf3u3bvl2LFjUrduXb/b9fq6detCPk9BQYHp79WwXbp0afnrX/8qF110kV/rxGWXXWbC7MaNG2XMmDHyq1/9ygRt3V7deuut0rFjR9OOsWjRIhk9erR8//33pnKtvvnmG/nuu+9kzpw58re//c0c5+233y5XXHGFfPzxx45fE0IxgKQssvtiyxbp2Db6XtZ4grEKDMfqROW4uAfWHt6cVOcCezydnKY5Gb3DyQzHkYJxqNc3WDA+sfCuuJUi23qMQwXiwLaJwBN1OJ04Ec8Cu1iqxLROZJf8/HypWrV47UOoft54aNvDypUrZf/+/aZSrD3JWtW1Wiuuvvpq37baC3zGGWdIy5YtTfX4wgsvNLfrYyx6v/Yn/+EPf5CJEyea4z1+/LgJ3RqIdaGdevbZZ6VTp06mMm21VERCKAaypHXCS5wG42Dh2B7g7FVN+21OZXIYdmMxXqSRbaGqxpF6jLNl8V2wFpBIgdgSLBBHM3EiknirxLROZA8NxPZQHEqtWrVM5XbHjh1+t+v1cP3A2mLRqlUr82+dPqH9vxpmA/uNLRqY9bk2bNjgC8WBdJGftk9s2rTJBN769etLmTJlfIFYtW3b1nzdvHmz41BMTzGQJbK5nzieHuNgo8XsVTd737HTQGzfR7hWiUwJxHZWn7PTfuNIPcahxrU5OR10YItBqs9SmIxAHPg6hfqcxdM2YWE2MaJRrlw5U3nVaq9FK7R6vXt35y1oVlU3lC1btpieYg26oWjlWcO2NfHi7LPPNiFZ2y8sX39d/BeNpk3927jCoVIMICP6iSNNU4ilH1aFGjEWipMwHHh8mV6NV05aVUJNpghXMVarC4oX6QS2U2Ri1ThUaI8mEAf+chYqEKeiSkzrhLeNHDlSBgwYYGYPd+3a1YxkO3DggG8aRf/+/U3/sFaClX7VbbUdQoPwe++9Z+YU66QIpS0V48ePl8svv9xUmzXU3nXXXaaybE2n0N7iJUuWyAUXXGBaMfS69gv//ve/l5NOOslso4v9tOf4hhtuMMekwXvo0KGmd9lePY6EUAxkuES3TqT6pB2OA1wUY8aC9Rs7le1hOJ6WCjf7jJ20U6h0CsduB+Jwn81wgTjRVWJaJ7yrX79+smvXLhk7dqxs377dtEPMmzfPt/hOWxW0gmvRwPzHP/7RVH8rVKhg5hW/8MILZj9K2zG+/PJLM6Zt79690qBBA+nVq5fcf//9vt5m/fryyy+bxX8arHVBnoZie5+xPqdOoBg+fLiZbKEnCtHFen/5y1+i+v5yinRGhkfoqJFq1arJjAmzpUL5iqk+HCDtR7FppditUOzWQrtwrBN8xHJyilgkIwhHqvwl+jW12E+eEun11WAcrgpvneQjcDqFdUrowMkUyppnHGqmdCrDsdMwHEsgjrZKHCkQa5U41kBsVYkJxf4O/XRQbhrbz0xZcNKbG21mefT1+VKhUiVx26EDB+TOyy5y/bgzGZViAJ5cgJcuYdjJn7tDBR39RSPc490MzNFUjd1YgKecVo1TWTkO1zuc7EAciRsj2AjEyGaEYiCDJbJ1ws0qcSp7YQM5CcuRFjG5HYSd/Lk7mHCPCwzMbgXkaEfixdpnHG07RbCQmsiAHE11ONoe4nBtOuE+eyyuA+JDKAYynNemTjgRLjhYvcfx7COVQdipwP1/4WJAjuZ00W70GatQM41VpIAcbzh2MvEiUnXYaSAOdoKOSL+oJaNKzAI7eAGhGICnJHMRnD2sJDoER2J/fntAjickx3KGwWjbKcJVjZVVOVZOqsfBBIbmaMa+hQrD4WYQhwvEwThpm0hGlZjWCWQ7QjGQoRLdOuGmZCyySyfRBJWUB+Q42yycjm6LNAov8EQq0fQaB4bjcAE5mFhmH8cShqPtIXYrEMezuE5RJYZXEIqBDJbI1olM7CdOtXSqDMcbkKMJx04X4YU7u6A9KAYLx4FV41Dh2Gn1OFqBJxIJfM5wrRJOqsOJDMRuoEoMLyAUA4AHw3Aw1rHHHY4j9Bo7qRorJy0VllAB2d57HEyo0OwkBDsJw4kIxBangTjeKjGBGF5BKAYyUCa1TmSzbAnDiQjH8VaNnSzEszgNyIECWy7swj0u2jBs/14Ce4djCcRO5hH79k/bBOAYoRjIULROpE62hmE3w7FbVeNI/cbhAnKwk4BEG3xDheBgz+1k1JrT04THG4hpmwCiRygGgCxaRJeO4TjSQrxgp92Opt84VEhdXVB8hrxgIgXmYAE4cP9OqsKpDMRUiYHoEIqBDJNpJ+zQIJUNvBiG3QrHgQvx3A7HTgKyXWC7RSihHh9tGHajfzialgm3AjFVYngNoRjIQJl2wo5MHsdGGE7vcGyx+o7DheRowq7TIBzsWJz2DiciEMc7fs1CIIYXEYoBGCyw817PsNvhOOp+Y4eL8eyniy5+3tCL8oJVkSOF5FAC9xHq+eIJw24H4njRNgEvIxQDGdY6wQK7xCEMxxmOE7QYz2n1OFxgDawkOxUu/IY7G12oMOy0OhxrIKZKDMSOUAwAGdIm4aQSmMpfbOIJx07OihcqHBc/d2NXwm20nFSFo60OB/6ClqxATJUYXkcoBpCQBXaZdHrndAzEoQJwpOCzNEJwTkZo1tcx1n7jaCvHwdorio8hcZ+7aKrCsVSHVbIrxCyuAwjFQMZI5NSJRMiEqRPpEoZjDcDRPkZ/+Qn2XIkIym7OOI4UkAODabRV5GgCcLDnc7s6TCAGUoNQDGSQRPQTJ6pKrNK1SpwuvcP2gOpGL2gkwZ4jMCi7/VlwJRw7bK1wEpCj5TT8xhqEY/3ljEAMuI9QDHhYoiZOpHOVOB2qw8kOw+EEPv/SBAXkVITjWENtrGINw6kKxBYCMVCMUAx4vHXCK1XiVFeH0ykIh2MdW6IqyG6HY4vTkJzqIOxGdditz5BWiQnEwAmEYsCjrRNeqhKnsjqcKWE4kP1YExGQ3QjHlmSG5MDnCXY8ifgsul0dZtIEUBKhGPCwbK8Sp7I6nKlhOJqAnKhwHO1nKFgotS/SCyVUcA4WfCM9XyI/i4kKxFSJAX+EYiDNZdLUiXSqEqeqOpxNYThiQHaxemx/n+IJyNEE13DBOZbgm4hfzAjEQPIQigGPtk5ka5U4HarD2RiGnfYfu1k9disgJzr4RjrxRjqEYUUgBsIjFAPIiipxOoRhLwXiZFSPkxmQUx2Ck/XLFYEYCI1QDKR560QmVImtQJyqkJIOC+m8GIaTWT0O9v5ap5O2S+Zn0O1fxBL1yxUVYsAZQjGAuHg1EHu9OpzKcBwyJAdUkt36XAbuM9wxpNtniUAMOEcoBjzE7Sqx1wMxYTi1rRWxhGS39p0JnyUCMRAdQjGQptJ56kSqw7AiEGeeZFSP7VJ51sJUf44IxED0CMVAGnOzn9itKrFXAzFhOHPDcTpJRtsNgRiIDaEYgGMEYgJxprZWpFqyetAJxEDsCMWAB7hRJSYQE4iTXT3O9JAc+H0k+jNEIAbiQygGsnwUm4aMeBGICcTJEvhaB4bkdA/IyQ7CijAMuINQDHhAPEHCi4GYMJzGITkgdKY6KKciBNsRiAH3EIqBLBZvlTgdArGFQBxbWFJunwAmnSvJbgfmcPsOdjzJRCAG3EUoBrJcvKEg1YFYq8QE4tCBNxwrLK10uH0mhudI71WwyrKb+0/1+08gBtxDKAY8cGrnWKvEqQ7EyZTugTgwCEcThpxsq5+7bKwup+v7GSuqw0DiEIqBLBVP64TVNpFqbp2NLJMDcbKqgoH7DqwuZ0tIzlRUh4HEIxQDWSyW1ol06iNWiW6dSMdAnA4ByP682VpFzhRUh4HkIBQDKCFdAnGipVsgTocwHAxV5NQgDAPJRSgG4Mk+4nQJxOkahMOhipx4BGIg+QjFQBaK5Qx26dJHnIypE+kWiDM5+BCQ3ZUNnwkgUxGKAaRdH3EipUMgzsTqsBME5Nhk6+cByDSEYiCNaJBIFQJx4nkp/IQLyMrrITmeEXsAEqNUgvYLIEbxhoVoR7GlYx9xIlon0iUQa/jxWgCyvmf7966vh/3iBYHfb+BrAmSCqVOnSrNmzaR8+fLSrVs3Wbp0achtn3nmGfnFL34hJ510krn07Nkz6PZfffWV/Pa3v5Vq1apJpUqVpEuXLrJ582a/bRYvXiy//OUvzf1Vq1aVc889Vw4dOmTuW7hwoeTk5AS9LFu2zPH3RqUYyEJO+4nTrY84UVIZiL1UHY55mkUWV5KpCCObzJ49W0aOHCnTp083gXjKlCnSu3dvWb9+vdSpU6fE9hpWr7nmGunRo4cJ0Y888oj06tVL1q5dKw0bFv//8caNG+Wcc86RwYMHy/jx403g1ft1e3sg7tOnj4wePVqefPJJKVOmjKxatUpKlSqu7er+v//e/y+t9913nyxYsEA6d+7s+PsjFAMel25V4mwMxAQh74TkYFVv3n9ki8mTJ8uQIUNk0KBB5rqG43fffVdmzpwpo0aNKrH9iy++6Hf9//7v/+S1114zYbV///7mtnvuuUcuvvhi+fOf/+zbrmXLln6Pu/322+XWW2/1e47WrVv7/l2uXDmpV6+e7/rRo0flrbfekuHDh5tqsVOEYgCOLNpV/KesHrWbZMxZ7FIViAnD3grJVIORyQoLC/2u5+bmmkugI0eOyPLly0211qKVWm2J0EquEwcPHjSBtUaNGub68ePHTai+6667TMX53//+tzRv3tw8R9++fc02O3fulCVLlsh1111nKsJaWW7Tpo08+OCDpsIczNtvvy0//PCDL7w7RSgGPCqa1gkNxF1aFleUF2080eeVqIDsRj8xgdhbITlZYTnU8/J+I1HK7S+ScseLXN/vsYPF+2zc2P+vhXl5eTJu3LgS2+/evVuOHTsmdevW9btdr69bt87Rc959993SoEEDE6StwLt//355+OGH5YEHHjDtFfPmzZPLLrtMPvnkEznvvPPkm2++MdvqMU2aNEk6dOggf/vb3+TCCy+UNWvWyMknn1zieZ599lkTshtF+bOEUAx4mJPWCatCbLHCcWBATkYVOVoE4uwU6vUNF5ad0EDt5PG8v8gm+fn5po/XEqxK7AYNvi+//LLpM7b6hbVSrC699FLTIqE09C5atMi0Zmgotrb5wx/+4Kv8nnnmmaYFQ9s2Jk6c6Pc8W7ZskQ8++EBeeeWVqI+RUAwgInsQDnX7so35fgE6loDs1tQJrRITiL0n3tfeOn017yG8RAOxPRSHUqtWLSldurTs2LHD73a9bu/nDUYrvBqKP/roIznjjDP89qmL5tq1a+e3fdu2beXzzz83/65fv/i/x2DbBE6oUM8995zUrFnTTLOIFiPZAJTQYNUKGdq7h7m8+Purpd2775bYpur338uI7mf7LsM/+dCE5FABOlmstolkIhBnB8ajAaHpYrZOnTqZCq1Fq7h6vXv37iEfpwvo7r//ftMWETgJQvep49d0eoXd119/LU2bNjX/1vFv2nIRbhtLUVGRCcW6iK9s2bISLSrFQBrRH8harUqnRUSZKJlVYgIxAK8YOXKkDBgwwITbrl27mpFsBw4c8LU1aBjVUWtWS4P2CI8dO1ZeeuklE263b99ubq9cubK5qDvvvFP69etn5g5fcMEFJjy/8847ps1C6fQI3UZ7ndu3b2/aK2bNmmX6mF999VW/4/v444/l22+/lRtvvDGm749QDHhQMk7Yoa0U0bRQuNE6key2CQIxAC/p16+f7Nq1ywRdDbgaUDXEWovvtJ3Bmh2spk2bZqZWXHHFFSEX8/3ud78z/cMapHXsmo5a07Ft9skSI0aMkJ9++sn0He/Zs8eE4/nz55cY3aYL7HRChU6niAWhGEAJa/fukN/F8XhtodAe42xum9BATBgG4DXDhg0zl2Cs6q5l0yZn/798ww03mEs4OqM42CxkO61Ix4NQDCDl4p1NnMzxa1SHASA7EYoBlGh7+F1D/9N1Nv1iiZTfW+B3W/l9+1x93nhbJwjEAIB4EIoBlLBu606/660/WmAu6VwlTjQCMQBkN0ayAR4VKoy6cQIO7SeOZj+xVomT1TZBIAaA7EelGEhDX8Y5lk3DYvuOzcKG0HCneT61uv9pPD+8d4z855JLSswpvuEy/xXF1uI6p4E4nokTBGIAgJuoFANpJt7wFU1IjLd1wW7rngIThqMJxLEiEAMA3EYoBjzKjdMp2zWuVM3xtlYgjucYCMQAADcRioE0ZYWyRC9Ac1Kx/XbnnhJtEqu+2xbTccUbiJNxgg4CMQB4D6EY8HALhRVMnQRjDcLWRVskOtVomJJAnGgEYgDwJkIx4HFOAmqrKjXN12h6huN5vlT1EROIAcC7mD4BeHgKRWAVt2Pbxubf29p3lKkfLPK7v0fA9vvq1S+xTaT9Z8IZ6wjEAOBNVIqBNJXMKRRuL7rL1IV1BGIA8K6MCMWbNm2SwYMHS/PmzaVChQrSsmVLycvLkyNHjqT60ICsoYHVzRFtmbawjkAMAN6WEe0T69atk+PHj8vTTz8trVq1kjVr1siQIUPkwIEDMmnSpFQfHpC2LRQaJpdG0UJhTurxvxBrtVLEyh6w03VhHT3EAICMCsV9+vQxF0uLFi1k/fr1Mm3aNEIxspqGtZVffR/XPmIKxlu2+PUYp2rKBIvqAADJkhGhOJiCggKpUaNG2G0OHz5sLpbCwsIkHBmQgGAc54K7WIKxclI1DtVyQSAGAGSSjAzFGzZskCeffDJilXjixIkyfvz4pB0XkM6TKKINxoFV43DbuIVADADw5EK7UaNGSU5OTtiL9hPbbd261bRSXHnllaavOJzRo0ebirJ1yc93fxERkAxuBrho+3Q19FrB1/q3/eL2cRGIAQCeqxTfcccdMnDgwLDbaP+wZdu2bXLBBRdIjx49ZMaMGRH3n5ubay5AtnClWvzt1qjmFydjbBuBGADg6VBcu3Ztc3FCK8QaiDt16iTPPfeclCqVEdPkgLRadBdvMHZbsk7KQSAGAESSEclSA/H5558vTZo0MX3Eu3btku3bt5sL4DVWwIuHFUI1lCZ67Fkw9uclEAMA0kFGLLSbP3++WVynl0YBf8ItKipK2XEBmVottodRq2qsEl05tgfwZJyymUAMAMiqUKx9x5F6jwGvcGNEW6Rw7HZATlZV2I5ADADIulAMwP1Fd4HsgdWN6nGyq8J2BGIAQLQIxUAGt1G4HYydtlY47UNOVRhWBGIAQDQIxUCGSnQwLlE9TmHl1wmqwwCArJ8+ASA4KwC6MZEiEg3C1iXdEIgBAPEiFAMZLpnBON3o90wgBgC4gVAMZAEvBmN7GCYQAwDiRSgGsoSXgjHVYQCA2wjFQBbJ9mBMuwQAIFGYPgFkGSss6gk+VKImUyQbYRgAkEiEYsADI9syORwzexgAkAyEYsALVeMEzzNOBCrDAIBkIhQDHpBJLRWEYQBAKhCKAQ9J55YKwjAAIJUIxYDHWypSGY7pFwYApAtCMeBRwcJxogNysFFxhGEAQDogFAMeFxhKrb5jt0JyYBAmBAMA0hGhGEDI0BpYRY53fwAApCtCMYCQCLQAAK/gNM8AAADwPEIxAAAAPI9QDAAAAM8jFAMAAMCRqVOnSrNmzaR8+fLSrVs3Wbp0acht165dK5dffrnZPicnR6ZMmVJiG+u+wMvQoUPN/Zs2bQp6v17mzJnj28+yZcvkwgsvlOrVq8tJJ50kvXv3llWrVkk0CMUAAACIaPbs2TJy5EjJy8uTFStWSPv27U343LlzZ9DtDx48KC1atJCHH35Y6tWrF3QbDbPff/+97zJ//nxz+5VXXmm+Nm7c2O9+vYwfP14qV64sv/rVr8w2+/fvlz59+kiTJk1kyZIl8vnnn0uVKlXMsR09elScIhQDAAAgosmTJ8uQIUNk0KBB0q5dO5k+fbpUrFhRZs6cGXT7Ll26yKOPPipXX3215ObmBt2mdu3aJjBbl7lz50rLli3lvPPOM/eXLl3a7369vPHGG3LVVVeZYKzWrVsne/bskQkTJkjr1q3l1FNPNcF9x44d8t1334lThGIAAACPKiws9LscPnw46HZHjhyR5cuXS8+ePX23lSpVylxfvHixK8eiz/HCCy/IDTfcYNojgtFjWLlypQwePNh3mwbhmjVryrPPPmv2cejQIfPvtm3bmvYMp5hTDAAAkKZy9x6V3MPOWwCcOn7oqK89wS4vL0/GjRtXYvvdu3fLsWPHpG7dun6363Wt1LrhzTfflL1798rAgQNDbmOF3R49evhu01aJhQsXSt++feX+++83t5188snywQcfSJkyzqMuoRgAAMCj8vPzpWrVqr7rodockkEDr/YJN2jQIOj9WgF+6aWX5L777itxu1aOzz77bPl//+//mfA+adIkueSSS0zPcoUKFRw9P6EYAADAozQQ20NxKLVq1TL9vdqna6fXQy2ii4b2/n700Ufy+uuvh9zm1VdfNYv3+vfv73e7BmWdUqFtHNrSYd2mUyjeeust09PsBD3FAAAACKtcuXLSqVMnWbBgge+248ePm+vdu3eXeD333HNSp04dU90NV0n+7W9/axbn2WlQ1jBs70O2rusxOkUoBgAAQEQ6ju2ZZ56RWbNmyVdffSW33HKLHDhwwEyjUFrBHT16tG97XfSmi+L0ov/eunWr+feGDRv89qvBVUPxgAEDQvYA62M+++wzufHGG0vcd9FFF8mPP/5oZhvrcel8ZD0m3dcFF1wgTtE+AQAAgIj69esnu3btkrFjx8r27dulQ4cOMm/ePN/iu82bN/vaF9S2bdvkzDPP9F3XPl+96Lg1XRhn0bYJfaxOnQhFx741atRIevXqVeK+Nm3ayDvvvGPmF2vVWo9Bn1ePrX79+uJUTlFRUZF4hI4aqVatmsyYMFsqlK+Y6sMBAAAZ7tBPB+Wmsf2koKDAUW9utJll+hPvSoUKlcRthw4dkJtvvcT1485ktE8AAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxQDAADA8wjFAAAA8DxCMQAAADyPUAwAAABHpk6dKs2aNZPy5ctLt27dZOnSpWG3nzNnjrRp08Zsf/rpp8t7773nd//+/ftl2LBh0qhRI6lQoYK0a9dOpk+f7rfN9u3b5frrr5d69epJpUqVpGPHjvLaa6/5bbNixQq56KKLpHr16lKzZk256aabzL6jQSgGAABARLNnz5aRI0dKXl6eCaHt27eX3r17y86dO4Nuv2jRIrnmmmtk8ODB8u9//1v69u1rLmvWrPFto/ubN2+evPDCC/LVV1/JiBEjTEh+++23fdv0799f1q9fb25bvXq1XHbZZXLVVVeZfapt27ZJz549pVWrVrJkyRKzv7Vr18rAgQMlGoRiAAAARDR58mQZMmSIDBo0yFfRrVixosycOTPo9o8//rj06dNH7rzzTmnbtq3cf//9psr71FNP+QXnAQMGyPnnn28q0Frh1bBtr0DrNsOHD5euXbtKixYt5N577zUV4eXLl5v7586dK2XLljVV7NatW0uXLl3MsWk1ecOGDeIUoRgAAMCjCgsL/S6HDx8Out2RI0dMCNWKrKVUqVLm+uLFi4M+Rm+3b6+0smzfvkePHqYCvHXrVikqKpJPPvlEvv76a+nVq5ffNlql3rNnjxw/flxefvll+emnn0yQVnrM5cqVM8dj0VYM9fnnnzt+LQjFAAAAaarU3v1S6sd97l/2FvfbNm7cWKpVq+a7TJw4Mehx7N69W44dOyZ169b1u12va89vMHp7pO2ffPJJU3XWnmINtlpZ1orvueee69vmlVdekaNHj5pe4dzcXPnDH/4gb7zxhmmXUL/85S/NPh999FET3n/88UcZNWqUue/77793/FqXcbwlAAAAskp+fr5UrVrVd11DZzJpKP7iiy9Mtbhp06by2WefydChQ6VBgwa+KvN9990ne/fulY8++khq1aolb775pukp/sc//mEW75166qkya9Ys0588evRoKV26tNx6660mgNurx5EQigEAADxKA7E9FIeiYVTD5o4dO/xu1+s6FSIYvT3c9ocOHZIxY8aYqu8ll1xibjvjjDNk5cqVMmnSJBOKN27caHqQdXGehl+lPccaiLWibE2quPbaa81F968TKnJyckwPtPYgO0X7BAAAAMLS1oZOnTrJggULfLdpf69e7969e9DH6O327dX8+fN922tLhF4Cq7kavnXf6uDBg+ZruG3stDpcuXJl04OsY+B0TJtTVIoBAAAQkbYn6KSIzp07m0kQU6ZMkQMHDphpFNbotIYNG/r6km+77TY577zz5C9/+YupBOsCuX/9618yY8YMc79WqPV+nU6hC+O0feLTTz+Vv/3tb6bKq3TGsfYOax+xVo+1r1jbJzRc69QJi1aTdUGeBmK9T/f58MMPmykVThGKAQAAEFG/fv1k165dMnbsWLOwrUOHDmYmsLWYbvPmzX4VXQ2pL730khmhpm0SJ598sgm0p512mm8bDcraB3zdddeZ6RIajB988EG5+eabzf06ak1P+KEL537zm9+YE3JoSNYe4osvvti3Hx3hpvOT9X4N0k8//bQ54Uc0cop0/oVH6KgRXVk5Y8JsqVC+YqoPBwAAZLhDPx2Um8b2k4KCAke9uemSWRJ13JmMnmIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5GReKDx8+bM6gkpOTIytXrkz14QAAACALZFwovuuuu6RBgwapPgwAAABkkYwKxe+//758+OGHMmnSpFQfCgAAALJIGckQO3bskCFDhsibb74pFStWdNxqoRf7ecQBAACAjKwUFxUVycCBA+Xmm2+Wzp07O37cxIkTpVq1ar5L48aNE3qcAAAAyEwpDcWjRo0yC+bCXdatWydPPvmk7Nu3T0aPHh3V/nX7goIC3yU/Pz9h3wsAAAAyV0rbJ+644w5TAQ6nRYsW8vHHH8vixYslNzfX7z6tGl933XUya9asoI/V7QMfAwAAAKRVKK5du7a5RPLEE0/IAw884Lu+bds26d27t8yePVu6deuW4KMEAABAtsuIhXZNmjTxu165cmXztWXLltKoUaMUHRUAAACyRUYstAMAAADE65XiQM2aNTMTKQAAAAA3UCkGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4MjUqVOlWbNmUr58eenWrZssXbo07PZz5syRNm3amO1PP/10ee+99/zuLyoqkrFjx0r9+vWlQoUK0rNnT/nvf//rt82ePXvkuuuuk6pVq0r16tVl8ODBsn//fr9tvvzyS/nFL35hnqdx48by5z//WaJFKAYAAEBEs2fPlpEjR0peXp6sWLFC2rdvL71795adO3cG3X7RokVyzTXXmBD773//W/r27Wsua9as8W2j4fWJJ56Q6dOny5IlS6RSpUpmnz/99JNvGw3Ea9eulfnz58vcuXPls88+k5tuusl3f2FhofTq1UuaNm0qy5cvl0cffVTGjRsnM2bMkGjkFGlE9wh90apVqyYzJsyWCuUrpvpwAABAhjv000G5aWw/KSgoMJXMTMkssRx3t27dpEuXLvLUU0+Z68ePHzdV2eHDh8uoUaNKbN+vXz85cOCACbKWs846Szp06GBCsEbQBg0ayB133CF/+tOfzP16PHXr1pXnn39err76avnqq6+kXbt2smzZMuncubPZZt68eXLxxRfLli1bzOOnTZsm99xzj2zfvl3KlStnttHjefPNN2XdunWOXxMqxQAAAGlKw2uiLlb4tl8OHz4c9DiOHDliqrDa3mApVaqUub548eKgj9Hb7dsrrQJb23/77bcmyNq30V8ENHxb2+hXbZmwArHS7fW5tbJsbXPuuef6ArH1POvXr5cff/zR8WtdxvGWAAAASAoNePXq1ZPbHhqUsOeoXLmyqfTa5eXlmdaDQLt375Zjx46ZKq6dXg9VjdXAG2x7vd2637ot3DZ16tTxu79MmTJSo0YNv22aN29eYh/WfSeddFLY18G3X0dbAQAAIGl0wZhWUrVCmyjavpCTk+N3W25urngVoRgAACBNg7Fe0kGtWrWkdOnSsmPHDr/b9bpWtIPR28Ntb33V23T6hH0b7Tu2tglcyPfzzz+biRT2/QR7HvtzOEFPMQAAACK2c3Tq1EkWLFjgu00X2un17t27B32M3m7fXukECWt7bXnQ0GrfRvuatVfY2ka/7t271/QzWz7++GPz3Np7bG2jEymOHj3q9zytW7d23DqhCMUAAACISMexPfPMMzJr1iwzFeKWW24x0yUGDSrue+7fv7+MHj3at/1tt91mJkX85S9/MX3H2qv8r3/9S4YNG2bu19aNESNGyAMPPCBvv/22rF692uxDJ0ro6DbVtm1b6dOnjwwZMsTMRP7nP/9pHq+TKXQ7de2115rQrqPfdHSbjo57/PHHzfFGg/YJAAAARKQj1nbt2mVOtqEL2LTFQUOvtaht8+bNZiqEpUePHvLSSy/JvffeK2PGjJGTTz7ZjEk77bTTfNvcddddJljr3GGtCJ9zzjlmn/a2kRdffNEE4QsvvNDs//LLLzezje0TKz788EMZOnSoqWZrq4ceo32WsRPMKQYAAEizOcVIPtonAAAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmEYgAAAHgeoRgAAACeRygGAACA5xGKAQAA4HmeOs2zdfI+PfsMAABAvKxM4aETBGctT53mecuWLdK4ceNUHwYAAMgy+fn50qhRo1QfBuLgqVB8/Phx2bZtm1SpUkVycnIkGxQWFpqgr/8xcs719MZ7lTl4rzIL71fmyMb3SmPUvn37pEGDBlKqFF2pmcxT7RP6Yc3W3+L0/1yy5f9gsh3vVebgvcosvF+ZI9veq2rVqqX6EOACfqUBAACA5xGKAQAA4HmE4gyXm5sreXl55ivSG+9V5uC9yiy8X5mD9wrpzFML7QAAAIBgqBQDAADA8wjFAAAA8DxCMQAAADyPUAwAAADPIxRnocOHD0uHDh3MWftWrlyZ6sNBEJs2bZLBgwdL8+bNpUKFCtKyZUuzIvvIkSOpPjSIyNSpU6VZs2ZSvnx56datmyxdujTVh4QAEydOlC5dupgzlNapU0f69u0r69evT/VhwYGHH37Y/HwaMWJEqg8F8EMozkJ33XWXOd0k0te6devMaceffvppWbt2rTz22GMyffp0GTNmTKoPzfNmz54tI0eONL+krFixQtq3by+9e/eWnTt3pvrQYPPpp5/K0KFD5YsvvpD58+fL0aNHpVevXnLgwIFUHxrCWLZsmfn/vTPOOCPVhwKUwEi2LPP++++bH+ivvfaanHrqqfLvf//bVI2R/h599FGZNm2afPPNN6k+FE/TyrBWIJ966ilzXX95ady4sQwfPlxGjRqV6sNDCLt27TIVYw3L5557bqoPB0Hs379fOnbsKH/961/lgQceMD+bpkyZkurDAnyoFGeRHTt2yJAhQ+Tvf/+7VKxYMdWHgygVFBRIjRo1Un0YnqbtK8uXL5eePXv6bitVqpS5vnjx4pQeGyL/96P4byh9aWX/kksu8fvvC0gnZVJ9AHCHFvwHDhwoN998s3Tu3Nn0rCJzbNiwQZ588kmZNGlSqg/F03bv3i3Hjh2TunXr+t2u17XlBelJq/nan3r22WfLaaedlurDQRAvv/yyaUfS9gkgXVEpTnP651pdkBDuoj+sNVDt27dPRo8enepD9jSn75fd1q1bpU+fPnLllVeaSj+A6CuQa9asMcEL6Sc/P19uu+02efHFF83iVSBd0VOcAX1yP/zwQ9htWrRoIVdddZW88847JnRZtOJVunRpue6662TWrFlJOFo4fb/KlStn/r1t2zY5//zz5ayzzpLnn3/e/KkeqW2f0NajV1991UwzsAwYMED27t0rb731VkqPDyUNGzbMvC+fffaZmeaC9PPmm2/K7373O/PzyP7zSX9e6f/n6cQk+31AqhCKs8TmzZulsLDQd13Dlq6Y1x/uunCoUaNGKT0+lKQV4gsuuEA6deokL7zwAj8U0oT+99K1a1fz1xfrT/NNmjQx4YuFdulDf3Tp4sc33nhDFi5cKCeffHKqDwkh6F8xv/vuO7/bBg0aJG3atJG7776blhekDXqKs4T+0LarXLmy+arzbwnE6RmItULctGlT00esFWZLvXr1UnpsXqfTW7QyrL35Go51dbyO+dIf4kivlomXXnrJVIl1VvH27dvN7dWqVTOzv5E+9P0JDL6VKlWSmjVrEoiRVgjFQAroXFVdXKeXwF9a+ONNavXr18/8kjJ27FgTtHRs1Lx580osvkNq6fhCpb9c2j333HNm0TEARIv2CQAAAHgeq3oAAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAYAAIDnEYoBAADgeYRiAAAAeB6hGAAAAJ5HKAaA/9m1a5fUq1dPHnroId9tixYtknLlysmCBQtSemwAgMTKKSoqKkrwcwBAxnjvvfekb9++Jgy3bt1aOnToIJdeeqlMnjw51YcGAEggQjEABBg6dKh89NFH0rlzZ1m9erUsW7ZMcnNzU31YAIAEIhQDQIBDhw7JaaedJvn5+bJ8+XI5/fTTU31IAIAEo6cYAAJs3LhRtm3bJsePH5dNmzal+nAAAElApRgAbI4cOSJdu3Y1vcTaUzxlyhTTQlGnTp1UHxoAIIEIxQBgc+edd8qrr74qq1atksqVK8t5550n1apVk7lz56b60AAACUT7BAD8z8KFC01l+O9//7tUrVpVSpUqZf79j3/8Q6ZNm5bqwwMAJBCVYgAAAHgelWIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgOcRigEAAOB5hGIAAAB4HqEYAAAAnkcoBgAAgHjd/wdRcUCBbriqvAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "from gbasis.evals.density import evaluate_density\n", + "\n", + "## compute electron density without screening\n", + "rho_unscreened = evaluate_density(rdm, basis, points, screen_basis=False)\n", + "\n", + "## plot electron density\n", + "X = points[:, 0].reshape(50, 50)\n", + "Y = points[:, 1].reshape(50, 50)\n", + "rho = rho_unscreened.reshape(50, 50)\n", + "\n", + "fig1, ax1 = plt.subplots(figsize=(7, 6), constrained_layout=True)\n", + "rho_plot = ax1.contourf(X, Y, rho, levels=np.linspace(0, 0.8, 10), cmap=\"viridis\", alpha=0.5)\n", + "ax1.set_title(\"Electron Density\")\n", + "ax1.set_xlabel(\"x\")\n", + "ax1.set_ylabel(\"y\")\n", + "\n", + "## add atomic symbols\n", + "symbols = [\"O\", \"C\", \"H\", \"H\"]\n", + "xy_coords = atcoords[:, :2]\n", + "for coord, symbol in zip(xy_coords, symbols):\n", + " ax1.text(\n", + " coord[0], coord[1], symbol,\n", + " color=\"red\", ha=\"center\", va=\"center\",\n", + " fontsize=14, fontweight=\"bold\"\n", + " )\n", + "\n", + "## colorbar\n", + "fig1.colorbar(rho_plot, ax=ax1, orientation=\"vertical\", fraction=0.046, pad=0.04)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8773b488-a18a-4b7b-813c-328d2bf7c62d", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "Now we calculate the electron density for the same system, but with screening enabled (the default) and using the default tolerance value ($10^{-8}$). The screening tolerance can be modified by setting `tol_screen = X` inside the `evaluate_density()` function. We then compare the screened and unscreened electron densities, calculate their difference, and plot it." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b163800d-f35b-499d-83aa-825346b1b393", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Are the screened and unscreened electron densities identical within 1e-4: True\n", + "\n", + "Maximum difference between screened and unscreened electron densities: 2e-09\n", + "\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAJjCAYAAADkjmzrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdi9JREFUeJzt3Qu8TPX+//HPdtnbnZBbRKVDShQR3U9C6aJOpauSVE6kdqeiTuiq4kgXkf5dnOSX6H5zEkeXQxFd6ESpROSW3HZC7P/j83XWbmbvmb1n9p6Ztb7f9Xo+HpNmZs3MmrW+6/t97+9813dl5efn5wsAAAAQYuX8XgEAAADAb4RiAAAAhB6hGAAAAKFHKAYAAEDoEYoBAAAQeoRiAAAAhB6hGAAAAKFHKAYAAEDoEYoBAAAQeoTikMjKypLhw4f7vRooxoknnmhuQXT55ZdLs2bNoh7btm2bXHnlldKgQQNTvq6//nrz+Nq1a+Xcc8+VOnXqmMfHjBkjYTdv3jzJzs6WH374we9VscLs2bNN2dF/4U6d9t///lcqVKggixcv9nW9gHgIxRZ75plnTMMR7/bRRx9lZD1+/fVXE7iD2IAtX748aptUrFhR6tatK507d5Zbb71VVqxYIUG1evVqs10/++yzlL6vvmfkNqlSpYrsv//+csYZZ8jTTz8tO3bsSOh97r33XlMG+/fvL88++6xceuml5vEbbrhB/vWvf8mQIUPM4927d5ewu+222+TCCy+Upk2b+r0qgG9atWolPXr0kKFDh/q9KkBMFWI/DJvceeedcsABBxR5vHnz5hkLxXfccYf5/6D2dGogOe2002TPnj3yyy+/yPz5800P5kMPPSRPPvmkXHDBBX6vorzzzjtFQrFuV+2hbdu2bco/b9y4cVKtWjUTgletWmWC7BVXXGG2yxtvvCFNmjQpWPaJJ54w2y7SrFmz5Oijj5Zhw4YVefyss86Sv/3tbylfZxvpHzXvvvuuzJkzx+9VAXx3zTXXmLr422+/lYMOOsjv1QGiEIodcOqpp0r79u3FFnl5eVK1atWMfuaRRx4pl1xySdRj+lN2165d5bLLLpNDDjlE2rRpI37Sn9czSYc4aK+5R3tvnnvuOendu7ecd955Ub80aA97YevWrTM9P7Eer1WrVsrW8/fffzeBPNPbJ1W091174vUPiEzQbbVz506pVKlSRj4P7tSTmdClSxfZZ599ZOLEiaZDBwgShk+EmPYOas9g/fr1JScnRw499FB56qmniiz322+/mZ/c//SnP5mGtmHDhnLOOeeYv/R1eMK+++5rltNeTe8neW/8so5F1d5IXVZ7B6pXry4XX3xxQaV/4403mh5J/fwWLVrIqFGjJD8/P+rz9f0GDBggr7zyihx22GEF6zp9+vQyfX/9KVt//tcA8cADD0Q9t2nTJjNG1ls37XW///77o3pLvaEZus4TJkwwvR667FFHHWV6oiOtWbNG+vTpI40bNzbL6DbU3lR9j1jj73Qoir6P0td521XXV3tmNaSuX7++yHe66qqrTCDVfVYaum90nPDHH38sM2bMiDmm2Bvv+f3338ubb74ZtW76r+6/sWPHFjxe2m2qPdbeNtWxiGrJkiUmzNeuXduURf1j8LXXXov6Dt56/Oc//5Hc3FxTPjVcnH322TG32dtvvy0nnHCCKZs1atQw233y5MlRy+j20GEgNWvWNMNNdHl9/0Rouf3zn/8ctS3UJ598It26dTN/mFSuXNn82qPHYyTdNvprRuvWrc331e+i66GvLXx86B80elzo9vKOjUSPcf21QMuV7hNdTvfRzTffXGQoTTLHYqKf/eOPP0rPnj3NPqpXr54ZfpPoEJ7SDBsqzCsvkceilvXTTz9dPvzwQ+nQoYPZ9gceeKD885//jHrtrl27TL138MEHm2V0HP2xxx4bdex45fb88883+0/3tdZ1OqSm8LppOb/oootMaNT38UyaNEnatWtnXqtlX3/ZWrlyZZHvkkg59T5r2bJl5rjW+kKX13pGf/UrLNHP9upAXU632QcffBBzP2jdpfXcq6++GvN5wE/0FDtg8+bNsmHDhqjHtNLTCjoePRlKe668Rk4raw0Hffv2lS1bthScNLV7927TOMycOdNUhoMGDZKtW7eaSl9PltC/+vVneB1XqqFDw7I6/PDDo3r6tPHXSl7DjlbWGpzOPPNM+fe//20+U4cH6M/3N910k2lMH3zwwaj11cbppZdekr/+9a8mvDz88MPyl7/8xYwJLu57lqRTp06mIo9sxLRh0MZE1+Pqq682vXz607eOkf3pp5+KnDimAUq3iS6r21MDtm6H7777rqCHVdf1yy+/lIEDB5oGV3tT9TN1/QufwKa051p7UbT3VoPucccdZx7XsdC6HfW5KVOmmH3n0XA/bdo081ll6SXUscHawOlwjlNOOSXmuulYYQ0vGvL1Dxt1xBFHFIwt1tdpj3Npt6n2rmqw1++ugUobY91+xxxzjOy3334yePBgE6JeeOEFE6hefPFFU/4i6bbWcKFhTwOPfoZuL91ukYFIg5sGNl0XDQiffvqpCXkaTrzhIPprjAYDfa9y5cqZ9dOgqw2/BoB49PvqPtZfKiLp/tdfKfS40++in6vrqGU8kh4buo76+frHih5L+pnaix/565Cuo24L/X4asrVMJXqMa/DWY1GPMd3eun8XLVpkjsGvv/7aBOBkj8VEP3v79u1y8sknm9ded9110qhRI1OG9Pv4TUOj/gGm66y/Jmmg1xCp5UDLixcwR4wYYfaNlgP9bvoHy8KFCwuOnS+++MIcv1oX6PbVfaOdBK+//rrcc889UZ+pv9BowNbx+l7ngC5z++23m1Ctn6N/2D3yyCNy/PHHm7Lq/SqTbDnV99M/xHT9dX3/3//7f+aPEv1D1ZPoZ+sQND2utX7Sfat1n5YpPW4jh2F5dB01FOv20j9EgcDIh7WefvpprTVj3nJycqKW1ceGDRtWcL9v3775DRs2zN+wYUPUchdccEF+zZo183/99Vdz/6mnnjKvHT16dJHP37Nnj/l3/fr1Rd7fc9lll5nnBg8eHPX4K6+8Yh6/++67ox4/99xz87OysvKXLVsWte7Z2dlRj33++efm8UceeaTYbfT999+b5UaOHBl3mbPOOssss3nzZnP/rrvuyq9atWr+119/HbWcfofy5cvnr1ixIuq969Spk79x48aC5V599VXz+Ouvv27u//LLLyWugzrhhBPMzTN//nzzOt3PhXXq1Cm/Y8eOUY+99NJLZvl///vfxX6O7iddTvdbLN76nn322VH7sWnTplHL6f0ePXoUeb2+9tprr416LNltWqNGjfx169ZFLXvyySfnt27dOv+3336LKoOdO3fOP/jgg4scF126dCkoo+qGG24wn7Vp0yZzX/+tXr262Y7bt2+P+izvdfqvvne3bt2i3kuPjwMOOCD/lFNOyS/Ou+++G1UWPC+//LJ5XPdxPLNmzTLLXHfddUWei1wXXaZcuXL5X375ZdQyiR7jzz77rHn9Bx98ELXc+PHjzXv/5z//SfpYTPSzx4wZY177wgsvFCyTl5eX37x584TKcjK8cl+YV1607EWWbX3s/fffL3hMy6PWqzfeeGPBY23atIl5DEQ6/vjjTTn74Ycf4u5Db90uvPDCqGWWL19uyuw999wT9fiiRYvyK1SoUPB4MuXU+6wrrrgi6j31eNe6LNnP3rlzZ369evXy27Ztm79jx46C5SZMmGA+J7JO80yePNk89/HHHxe77YBMY/iEA/Snau11jLxpr0w82rZpz5rONqD/r73M3k17dLXnWXsOlC6nPU/a61ZYrJ8i49Ge5EhvvfWWlC9f3vQORdJeR12nwuuvPdKRJ2VoT7T2MGiPRFnp8A6lvb1q6tSppmdHexkjt42ug/acv//++1Gv79Wrl1nW4/XqeuumPyfqeFgddqAn+aWC9sLqT6Xa4+TRn8+1V0Z7ZFO5PVIh2W2qPY/esBy1ceNG0xOmPVa6Xt7rf/75Z1Nmv/nmG9MrG0l75SLLqH6+fpY3LZoeJ/pe2lNbuGfde52eJKfvrb3G+lne5+rQH+3h1PUufAJiJH2Niiwfyuth0xMa9Sf4WPTY0/UofCJj5Pp5dJ9Hju9O5hjXfaO9wy1btoxaTnsYlf6ak8yxmMxnaz2gQ4m0R9ajvyTpvvObbk/vWFZaHnXYQ2Sdo/tRf8HQMhKL9qxqGdFfI/TXkZLqTz0JLZL2yGv50nIfuR11GkTtUfb2TWnKaeHP0u+qr9Xe22Q+W3vG9ZcPfb/Icf/aq67DMmLxjofCv3ACfmP4hAP0Z7FkTrTTilrHd+pP5HqLRSs5paFLGwKdW7K09LX6M3skDSb6U6n+/BpJG2fv+UiFGxSvYk1FyNT5dpW3Ltq46E+ekaEs1raJt25ehe+tm/78rz9JauDX8ZX6s7IOSdFgqw1MaWgQ158pNQjrEAsNGhqwdEhDMn+sJLI9UiHZbVp4NhX9KVsDlv6Uq7d476FDKxLdL94fFDo2trj1VvrzeTy67QuH3sIKj5PXEKvBX8ej6jAFHWOpw0A01Gh58dZPjxH9CbokhbdXMse4fsevvvqq1OW98LGYzGfrca7jmAuXWa1zSqLDhfSPpUj6HfSP7VRIpM7RYUx6boCeb6HlSMfz6vAhb/iYF6CLK2PF7UfdN1p2NITG4g3PKk05Le740D9yEv1sr64uvJw+r+OwizseylpXAalGKA4hr8dAZ2OIV4lGjgkuK23kdXxbWcRr6AqHjdLQsdE6ls4b26bbR8cD6olGsWgDmOy6aYDVnjMdn6ljpzXY6Vg+7f3UsbjJ0gZMg7UXinUssZ6cVHiGjdLwJtZP5ZR+yW5T7V0v/Hql07xpb2Mshdc3FWXG+9yRI0fGnRbP61mPxRtjW/iPNw0Dus90bLCOLfWmw/vHP/5hHivuPWOJt70SOcZ1WT2Rb/To0TGXKzwmtKTtmqn6Rcekn3TSSVGP6cmfscboFxfA9NeDWBIpPzq2Vv940fGxOgZfx+XqHznjx483Y3CTFWs/6nrrL2ex1scrJ6Upp4nsx0Q+uzS84yFy9hsgCAjFIaS9KdoLqI2B/hRaHP2ZVH+m1594Y03LVdq/9nXmB527VX++juyR1LO0veczYe7cuaZRiwyT+p21t7SkbZMsfV/tLdab9sJo46UhSM/uLs121Z5m7aXSmS40HGu49k4AKgs90UnFC5+lUdZt6vU4aRlM1X7xhgDoHwHx/gDwltE/mErzuTokwQtrseivBnrTE5r0hE2d/eP55583gUo/W8Oy9oYm0ltclmP8888/Nz+zp6LnLpnP1uNct78GscjPXrp0aYmfo1MoFp7lobhfXryeUO3FjpwysKxXGdR9ozM36E3LuAZlPQFP96FXbkt7BTfdN7pttAe58B+OhZcrSzkty2d7dbXWad6QG6Vthpb7WFNd6uPaUVLc+wJ+YExxCOlf/frTrY77i1VZR05bpcvpuK9HH300bo+CjgH0GptE6fRs2mgWfl/tZdHGUc+iTjdtDHXcm46D01kvPDqGTsOyBpLC9DvqDADJ0JkXCk+Rpg2OBofipp7y5iiNt111G2lPiw7NeO+991LSS6zBTHu7dFYODUmpUtZtqj35OsTg8ccfN7NVFBZrqrWS6OwPug+0x77w/vHKtp4lr/tKZ03xhpUk87k6nEN7WiOnUPN6ygr3WHs9fF6Z0GNPl/EujBNr/VJxjOu+0fHYeoGWwnR2CB2XmoxkPlvrAb1IjfaaRx4v8YZdFA65GgAjb8XNuuIFx8jx6/rddL7c0vLGjEf2nuofWN4+1D8QNCTrzBWFr56ZyC8WOouNbk8tA4WX1/ve55e1nJbls3Xonn5P7R3XIS0enTUlXt21YMEC8wd8vDHHgF/oKXaA/rzl9bBG0ulx4o3puu+++8yJEh07dpR+/fqZk0q0R0pPgNEeXG+snvZG6tycOt/rvHnzzMkY2pDoMjolk/ZU6k9++nqd6kr/8teeEx1DV9w4Oh1KoD996lydOhWV9iboz4/6M6QONUj1lY70e2mPrP4kqBW19q56JzJpz2jkz7kakHXuWx2e4E3BpN9Zp6nSxlvXN5mf/XRaKw2YGj50O+kY65dfftlMW1XclfR0G2iPljY2Gt40JOv+8sYdaq+pvl7/sNDGS6/alwz9LtqIa0PmXdFO5zTVfaEnX6VSKrapnlCq09HpT/1aZrVs6zbUsK1z3WpvZzK0V03/CNMePZ2b2JsfVt9Hg5mGJe3N0j8S9A8QbcS1N1CDrm4vPX70PXT4Q3H0GNH9Hdkbqu/92GOPmWnkdD/rLyYaSvX9NCgqPT50fKpOeaa9cDpeVcuvTq+lz0VOx1eWY1w/Q6dz0xOldHmd9k7/YNU6RR/XcpHsxYES/Wx9Tsuv1jMalPSkOz0evT+0U0n/CNJxtDrFmpZHPWY0rGqgK+3l3vV76R9rWp613tM/frQ8R+4b3X9abnVaPj2BUI9fLe86x3dJl3DXsnH33Xeb6QL1NTruXOsC7WnVMqXvp0OKUlFOS/vZWg/pcjolm/YU6/kOuoxOBxer/dEeZP0jXtsPIHAyPt8FMjIlW+GpvGJNmbZ27VozdVaTJk3yK1asmN+gQQMz7ZVOpRNJp/W57bbbzNQ+3nI6ddq3335bsMycOXPy27VrZ6ZrivwsncpLp+KKZevWrWaarEaNGpn31SmFdNqyyCmF4k3x5U2bpO9fHG+KL++mUwnVrl3bTMM1ZMiQItMkRa6bPq9TQ+l3qlu3rpn6a9SoUWYKosj3jjXVWuQ20GmpdP1btmxptoVOSaWfHzkNVawp2bzp3Vq1amXWO9b0bPPmzTOPd+3aNT9R3pRM3q1SpUr5jRs3zj/99NPNFHyRU56lakq2VGxTpWWud+/epgxqmdlvv/3Mek+bNq3IcVF4ujOd3ivWNF+vvfaaWY/KlSubqeA6dOiQ/3//939Ry3z66af555xzjpmySqfl0u9+/vnn58+cOTO/JAsXLjSfGznlmT6m02/tv//+5v10Siv9Hp988knUa3///XezLbTs6Dbbd99980899dT8BQsWlLi9kznGdfvff//9+YceeqhZn3322cccz3fccUfBVIXJHouJfrYeg2eeeWZ+lSpVTJkYNGhQ/vTp01M+JZvS7abHnm5L3fY61WS8Kdlile3Cx6hOKanlpVatWqb86H7Sqcq88uxZvHixmfJMl9PjrUWLFvm33357wtMkvvjii/nHHnusqT/0pp+j23bp0qVJl9N4nxVrOyTz2Y899phpI/Rz27dvb6azi1Wnvf322+Zzvvnmm5jfFfBTlv7H72AOoHS0V1N/dtfefO3xQzDpLwXehSmAMNMeZ/3FRHubgaAhFAMW059p9ad4vYy0NwYZwaMnq+rQIx0GkamTSIGg0an/dPiTDhtJdJo6IJMYUwxYSMcH/ve//zUnJGkwJhAHm46tjTwJCQgjnYc+2ROVgUyipxiwkM7FqieZ6bRp+pN8Ki+0AQBAGBGKAQAAEHrMUwwAAIDQIxQDAAAg9EJ1op1OfK9XT9Lxl6m4nCkAAAg3HYWqF+DRaRf1QiqwV6hCsQZiveQqAABAKq1cuVIaN27s92qgDEIVir0z9LXg6mUv1ZSHpovNeg3q7vcqZGQbBuF7IvVsP/4yjeMACJ4tW7aYDjdmAbJfqEKxN2RCA7HeJo18QypXqiK2uuSm0/1ehYxtQ++PGLjD9uMvjMc7gPgYlmm/cmFukG3mdwOp2y+T29D2/QXYfLwDQBiEMhTb/pOt3w0kARVlRRmy53gHgLAIZSi2md8NJGEGZUUZsud4B4AwIRRbxO8G0u8w4/fno+zYh/Yc7wAQNoRiS/jdQBJmUFaUIXuOdwAIo1DNPmErPxtIggxSgXKUOAKx/+WUfQCEE6E44AjERdeJBssuQSxHQUXZzjzKJwAPodgHNjR8NBRIBcqRnXVAIvst6N8hEZRPAJEIxRlmQ0NCQ4FUcLUc2XAMp+o7uroPASAWTrTLIBsaUxpBpIKL5UiPXxuO4VRy/Tu7/N2AdBk7dqw0a9ZMKlWqJB07dpR58+YVu/yYMWOkRYsWUrlyZXM57BtuuEF+++23pN5Tl7/22mulTp06Uq1aNfnLX/4ia9eujVpmxYoV0qNHD6lSpYrUq1dPbrrpJvn999+T+m6E4gyxofK1JcjYsp5h413l0MX9Y8Pxm04uh2NXvxeQDlOmTJHc3FwZNmyYLFy4UNq0aSPdunWTdevWxVx+8uTJMnjwYLP8V199JU8++aR5j1tvvTWp99Qg/frrr8vUqVPlvffek9WrV8s555xT8Pzu3btNIN65c6fMmTNHJk6cKM8884wMHTo0qe+XlZ+fny8hsWXLFqlZs6ZMuHOKVK5UJWOfa0Ola1uQsWGbus62MlNalDX3T3gtXJZd+m7IXLbYvHmz1KhRQ1zWsWNHOeqoo+TRRx819/fs2WN6fwcOHGjCb2EDBgwwYXjmzJkFj914443y8ccfy4cffpjQe+p23XfffU3APvfcc80yS5YskUMOOUTmzp0rRx99tLz99tty+umnm7Bcv359s8z48ePllltukfXr10t2dnZC348xxWlmQ+UalnCDsgljObHh+PWDa9tFv08YyzeCT4cNaO9numi/aFZWVtRjOTk55laYrseCBQtkyJAhBY+VK1dOunTpYsJpLJ07d5ZJkyaZ4RAdOnSQ7777Tt566y259NJLE35PfX7Xrl3mMU/Lli1l//33LwjF+m/r1q0LArHS3ub+/fvLl19+KUcccURC24NQHPKGw9aGwLWeqqCztZyUBeUrXAjGCGIgbtqsnqxbuzVtn6Hjc7dt2xb12LBhw2T48OFFlt2wYYMZphAZPJXe157bWC666CLzumOPPdYEcB3je8011xQMn0jkPdesWWN6emvVqlVkGX3OWybWe3jPJYpQHOIGlQYAsPPYRXqw7xEk2ouqgfiLL2+T6tUrpfz9t279TQ4/9B5ZuXJl1LCPnBi9xKU1e/Zsuffee+Wxxx4zwySWLVsmgwYNkrvuuktuv/12CRpCcUgrVgIxEhWmsmLDsQsgXDQQV6+R+lDs0UCcyFjounXrSvny5YvM+qD3GzRoEPM1Gnx1qMSVV15p7usQh7y8PLnqqqvktttuS+g99V/9A2HTpk1RvcWFlyk8Y4X3nvHWLRZmnwhho+pKyHHleyAYbDh2AcAv2dnZ0q5du6iT5vSkOL3fqVOnmK/59ddfzRjhSBqClQ6nSOQ99fmKFStGLbN06VIzBZu3jP67aNGiqBkrZsyYYcJ+q1atEv6O9BSHrFElSCIZYSkvNhy7COaxQNlBmOTm5spll10m7du3NyfO6RzE2vPbp08f83zv3r1lv/32kxEjRpj7Z5xxhowePdqc6OYNn9DeY33cC8clvafO7NG3b1+zXO3atU3Q1ZkpNAjrSXaqa9euJvxqr/QDDzxgxhH//e9/N3MbJzMchFBcDNcqu7AEHADBYftJsSXVm7Z/PyAZvXr1MlOc6fy/Gjzbtm0r06dPLzipTXtvI3uGNZjq7Bb676pVq8zUahqI77nnnoTfUz344IPmffWiHTt27DAzS+g4ZY8G7DfeeMPMNqFhuWrVqiZo33nnnUl9P+YpjsO1Ss7FQOzaPgo6F8uQohwhFeWechRe6Zqn2Hvf71fclZYxxVu3/CYH7H97KOZXThRjikNQubkYZlzbR/AH5Qip4mI9C4QNodjxRtLFitq1fQR/UI6Qai7Wt0CYEIodbiRdrKBd20c2cWnbu/RdEKyy4mK9C4QFodjRRtLFitm1fQR/UI6Qbi7Wv0AYEIodbCRdrJBd20fwpwxRjpApLtbDgOtCH4pdayRdrIhd20fIPMoQ/OBifQy4LNTzFLvWULpQAbu2T1yj+8e2ckaZgp+YxxiwR2hDsWuVlG1BxcV9gOChjCEICMaAHUIZinsN6i4uIRADlCkE+xeSZF9PeQYyL5Sh2CUEYoSpvFF2EBb0LgOZRyi2mN8BpTSo5O0SlDJGuUEYEYyBzCIUWyooYQVuCWK5IhQgzAjGQOYQii0UxOCSCCr24AlyWaK8AHsRjIHMCP08xbYJcogpDhV68AS1LHGRDQSR32UyqMcr4BJCsUVsrRT9bkxgB8Iw4GYbANiC4ROWsLUyJOSEszwlut+99aCcAIlhKAWQPvQUW4BADBcDsbcs5QQIR5sABB2hOOBsrfwIOsEUpEAM2CZI5dvWtgEIMkJxgNla6QWp4QAAV9naRgBBRSgOKFsrOwJxcNFLDLjH1rYCCCJCMVKGUBTO/cO4YIRJEMs6wRhIDUIxnG0okD5eEGa/I4yCWO4JxkDZEYrhZAOB1O8rgjDwhyAeBwRjoGwIxXCuYUBq9xlBGLAHwRgoPUIxSo2g5DbCMGAngjEQslB83333SVZWllx//fV+r0ooEZbc3H8MkQASF+TjRIMx4RgIQSieP3++PP7443L44Yf7vSqhFOSGAKXbjwRhwE0EY8DhULxt2za5+OKL5YknnpB99tnH79UJHYKTWwjDgPsIxoCjofjaa6+VHj16SJcuXUpcdseOHbJly5aoG0qP8AQAdtaLBGOgZBXEIs8//7wsXLjQDJ9IxIgRI+SOO+5I+3qFgS0VPwCg7MGYOh9hZE0oXrlypQwaNEhmzJghlSpVSug1Q4YMkdzc3IL72lPcpEkTCbqg/UVP5QgE7/jmuES6yyllDGFjzfCJBQsWyLp16+TII4+UChUqmNt7770nDz/8sPn/3bt3F3lNTk6O1KhRI+oWdARiADbWFWHmaj1JGUPYWBOKTz75ZFm0aJF89tlnBbf27dubk+70/8uXLy+2C1oF5GpFD7giaHVGmLlaXzK1G8LEmlBcvXp1Oeyww6JuVatWlTp16pj/t13QKh1XK3jANUGrO8LM5XqTcoYwsCYUuyxolY3LFTui0QvkBvZhcLg8zSHlDK6zOhTPnj1bxowZIzYLWiXjamWO4hGO7cf+CxZX61LKGVxmdSi2XdAqF1crcdhbJl2UzuOM/RcsrtaplDO4ilDsk6BVKq5W3kgevcZ2Y98Fi6t1K+UMLiIU+yBolYmrlTbcKqdIHPsumOOMXRtvTDmDa6y5eIcrglaJuFRBI33llXJiHy6+EFyx9kvQ2oZEUc7gEkJxBgWt0qMiC7dkyiPh2I19zP4LrpL2TdDaj0gEY7iCUJwhQavQqMCA8CmpHqJeCK6y7pt0t0EEY7iAMcUZQCCGCyg3bm/LIK4TUicT45mD1tYBySIUp1nQKgkaPpSmLFBu3Mb+DQ+CMRAfoThElQMNH0qDcuP2dg3KeiBzCMYoi7Fjx0qzZs2kUqVK0rFjR5k3b17cZU888UTJysoqcuvRo0fBMi+99JJ07dpV6tSpY5777LPPirzPhAkTzHvVqFHDLLNp06Yiy9xzzz3SuXNnqVKlitSqVatU341QHJJKgYYPpUG5cRv7N7wIxiiNKVOmSG5urgwbNkwWLlwobdq0kW7dusm6detiLq+B96effiq4LV68WMqXLy/nnXdewTJ5eXly7LHHyv333x/3c3/99Vfp3r273HrrrXGX2blzp3nf/v37l/r7caJdCCoDGj77+TH7A+XGbexfaBlIZ3vFyXfuGT16tPTr10/69Olj7o8fP17efPNNeeqpp2Tw4MFFlq9du3bU/eeff9705EaG4ksvvdT8u3z58rife/3115t/Z8+eHXeZO+64w/z7zDPPSGnRU5xiBGKks0xl6mpzlJvM8Gs7s3/hoccYW7Zsibrt2LEjbk/sggULpEuXLgWPlStXztyfO3duQp/15JNPygUXXCBVq1aVIKKn2OGDn4YvHGUt1fuZcuM29i8Ko8c42L7dVVeq7qyc8vfN27Xd/NukSZOox4cNGybDhw8vsvyGDRtk9+7dUr9+/ajH9f6SJUtK/Dwde6zDJzQYBxU9xSlCIIZf5SqVvceUG7exfxEPPcbhtXLlStm8eXPBbciQIWn5HA3DrVu3lg4dOkhQEYodPNhp+MKprOGYcuOPTG139i9KQjAOJ53RIfKWk5MTc7m6deuak+TWrl0b9bjeb9CgQbGfoSfT6Xjivn37SpAxfMKxg5yGzx2lLVuleR3lxm3sXwRlKAXslZ2dLe3atZOZM2dKz549zWN79uwx9wcMGFDsa6dOnWrGKl9yySUSZIRih9DwuSOTjRLlxn/sA4QlGPsxkw5SJzc3Vy677DJp3769GQYxZswY0wvszUbRu3dv2W+//WTEiBFFhk5okNa5iAvbuHGjrFixQlavXm3uL1261Pyrvc9eD/SaNWvMbdmyZeb+okWLpHr16rL//vsXzHCh7+G9l4599uY7bt68uVSrVi2h78fwCUdQwbiDQAzAbwylQCy9evWSUaNGydChQ6Vt27YmeE6fPr3g5DsNpDofcSQNuR9++GHcoROvvfaaHHHEEQUX9NDZKfS+Tvfm0f/Xx3Q6OHX88ceb+/paj66TPqYnCm7bts38v94++eQTSVRWfn5+voSETjVSs2ZNM5Bcx824cmATbNxBIAYQpjrJhXooHdki8n3f/XacVK2ehtkntm6XLgf1T/l624yeYsu5UKFgLwIxgKChxxhhQii2GMHGHQRiAEFFMEZYEIotRbBxB4EYQNBp3ZHO+oNgjCAgFFuIYOMOAjEAmxCM4TJCsWUINu7IVAOQ7h4eAOFCMIarmKfYIgSb4AtahU6ZAWDjXMbUXfADPcWWoIIIPgIxgDChxxiuIRRb8vM3kAzKTDDqh8gb4CKCMVxCKC4lAjGCWnlTZvwXpPIApBvBGK4gFAMOVdoEYv8FqTwAmUIwhgsIxaVALzGCWFlTXvwXpPIAZFo6Z7rh2EImEIqTRCBGECtpyov/glQeAD8RjGErQnESCMQIYuVMefFfkMoD4HKvMcca0olQnCACMYJYKVNe/Bek8gAEDcEYNiEUJ4BAjCBWxpQX/wWpPABBRTCGLQjFJSAQI4iVMOXFf0EqD0AYh1NwDCLVCMXFIBAjiJUv5cW+8sA+A/YiGCPICMVxcKAhiGWBcOW/IJUHwEYEYwRVBb9XIIgyeYARcoItSJUtZcUfQSoDgCu0PkvlsaXvRR2JsqKn2EccwMEWpDBEWfFHkMoA4Bp6jBE0hGIg4JUrgRiAq1J9Al6Q6m7Yh1DsE4JOcAWpUqWcAAgDgjGCgFAMBLQyJRADCJNU9hoHqS6HPQjFQAArUQIxgLBKVTgOUp0OOxCKfUDgCZ4gVZ6UDwBITV0YpLodwUcoRugFqdIkEAdHkMoFEFYEY2QSoRihFqTKkkDsJvYrUDYEY2QKoRihFaRKkuAEAPERjJEJXNEOzrC1wiMQA0BmroLHle9QHHqK4QQCMQC4jzoT6UQozjAO6NSzNRAjuChTQHiuggd4CMUZxEGceoQXAAin0oZj2g3EQyjOEAJx6lGxIR0oV4BdaF+RKoTiDOCATT2CCwDAQzuLVCAUpxkHauoRiAEAhdHeoqwIxWnEAZp6BGKkE+ULsBvtLsqCUJwmHJip51pgoYwAQOpRt6K0CMVpwAGZeq4FYgBA+tAOozQIxSnGgZh6BGJkAuUMcAvzGSNZXOY5hTj4Uo+gApvLl6t1guvHZTL7LQjbwtVyls7LQ3O5Z8RCKE4RDi7Af0EIKK4Lwza2LTDZtr5BCMZsL8RCKAZ8QkOWnDCEsaAL0z6w7fi0bX39DMZsJ8TDmGIAgaaNWJjCWFCFcR/Y9p1tW18/EIhRHEIxgMAGYRr5YAjzfrDtu9u2vkCQMHwCsECQGrp09bQE6TviD+yX+EMTgrptGCKAdBo7dqyMHDlS1qxZI23atJFHHnlEOnToEHf5TZs2yW233SYvvfSSbNy4UZo2bSpjxoyR0047zTw/YsQI89ySJUukcuXK0rlzZ7n//vulRYsW5nl9zbBhw+Sdd96RFStWyL777is9e/aUu+66S2rWrGmW+fzzz+W+++6TDz/8UDZs2CDNmjWTa665RgYNGpTUdyMUAz4KaqNanFSdrGLjdw8b9pG9QdO29YUdpkyZIrm5uTJ+/Hjp2LGjCbfdunWTpUuXSr169Yosv3PnTjnllFPMc9OmTZP99ttPfvjhB6lVq1bBMu+9955ce+21ctRRR8nvv/8ut956q3Tt2lX++9//StWqVWX16tXmNmrUKGnVqpV5vQZefUzfUy1YsMB8xqRJk6RJkyYyZ84cueqqq6R8+fIyYMCAhL9fVn5+fr6ExJYtW8xfFZs3b5YaNWqktCGg8kkPGmV7lHQMhG1f2l4nhG1/lWa/2rCNbC+HLmWL0r7vu9+Ok6rVK0uq5W3dLl0O6p/Uenfs2NGE10cffdTc37NnjwmhAwcOlMGDBxdZXsOz9iprL3DFihUT+oz169ebgKth+fjjj4+5zNSpU+WSSy6RvLw8qVAhdv+uBu2vvvpKZs2aJYmipxhAWnqQbQgMiI19504PrG3ri8zT8B0pJyfH3GL1+mqP7JAhQwoeK1eunHTp0kXmzp0b871fe+016dSpkwmor776qhn6cNFFF8ktt9xienFj0ZCuateuHXedvSAfLxB7yxT3HrEQigGkHKHKXuw794KmbeuLaMt/rSNVyldJ+fv++uuv5l/t6Y00bNgwGT58eJHldazu7t27pX79+lGP633tCY7lu+++Mz21F198sbz11luybNky+etf/yq7du0yn1OY9jxff/31cswxx8hhhx0W8z11PXQ8sQ6PiEeHT+hQjzfffFOSQSgGABgEYne3FcEY8axcuTJq+EROjF7i0tKQq0MhJkyYYHqG27VrJ6tWrTJDKmKFYu1RXrx4sTlhLl6vdo8ePczY4ljBXenrzzrrLPP+OjY5GYRiAEhjaPIziNgW3BCc8hD2AB2E4zdTNBAnMqa4bt26JtiuXbs26nG936BBg5ivadiwoRlLHDlU4pBDDjEzV+hwjOzs7ILH9YS4N954Q95//31p3LhxkffaunWrdO/eXapXry4vv/xyzDHKenLeySefbHqR//73v0uyCMUpEIaDxg806HCBX40rx09mJLpfbdsfYe9ZDvN3j0cDrPb0zpw500yJ5vUE6/14MzzoMIjJkyeb5XT8sfr6669NWPYCsc73oCfqadCdPXu2HHDAATF7iHWWC+3F1nHKlSpVKrLMl19+KX/+85/lsssuk3vuuUdKg1BcRhw46WFbAwJkaiq7ZD4HAFIpNzfXhM727dubuYl1SjadAaJPnz7m+d69e5tp13TuYdW/f38zU4XOF6zB95tvvpF7771XrrvuuqghExqc9UQ87QXWXmSlM2/ovMUaiHUYhI6B1inX9L53cqCeuKe90DpkQgOxBmddR+899DldJlGE4jIgEKcHjTpcl66AzLEDIJ169eplpkwbOnSoCZ5t27aV6dOnF5x8pxfX8HqEvZP4/vWvf8kNN9wghx9+uAnMGpB19gnPuHHjzL8nnnhi1Gc9/fTTcvnll8vChQvl448/No81b948apnvv//eXKhD5yvW9dLQrDePXihk+fLlCX8/5ikuZcNCIE4PGnWEVVnrFI4df7g6fELRzgVjnuL/t+gFqVI9DbNPbP1Vrmx9fsrX22Z/xHkkjIoiPWxsNIBUlv/SHgMcO0gHyhXChlCcZOglEKcHlS9QunDMsQMAqUEoTgKBOD1o1IHShWOOHQBIHUJxggjE6UGjDpQuHHPsIBMoZwgTZp9IAIE4PahsgcRxvABAeoUyFE95aLpUrlQlocBLIE4PGngAABAkoQzFiYYzAnF6EIgBIDxXtytc59O2IqgYU4yMIhADgH1SOV0g7QCCilCMjKEiBIDw1OHFLU97gCAiFCMjqAABIDx1eSLL0S4gaAjFSDsqPgBwRyrnz6Z9QJAQipFWVHgA4J54Y4VLU+fTTiAoCMVIGyo6AAhHHV/W+p72AkFAKAYAAL4HWoIx/EYoRlpQuQGA+1Jd19N2wE+EYqQclRoAoLRoQ+AXQjFSisoMAFBWtCXwA6EYKUMlBgBIFdoUZBqhGAAABBLBGJlEKAYAoJQIbenHNkamEIoBAECgEYyRCdaE4hEjRshRRx0l1atXl3r16knPnj1l6dKlfq8WACDkCGyZwXZGulkTit977z259tpr5aOPPpIZM2bIrl27pGvXrpKXl+f3qgEAgAwgGCOdrAnF06dPl8svv1wOPfRQadOmjTzzzDOyYsUKWbBgQdzX7NixQ7Zs2RJ1AwAg1QhrmcO2hoQ9FBe2efNm82/t2rWLHXJRs2bNgluTJk0yuIYAACAdCMZIBytD8Z49e+T666+XY445Rg477LC4yw0ZMsSEZ++2cuXKjK5nmFBBAQAyiXYHqVZBLKRjixcvXiwffvhhscvl5OSYGwAUdslNp8d8nIYWpaVlJ165QnqwzRHqUDxgwAB544035P3335fGjRv7vToALJJI41l4GUIykkF5yTyCMUIXivPz82XgwIHy8ssvy+zZs+WAAw7we5UQokYgkxVuGLanTfvPey37BWE5Nmws6wRjhCoU65CJyZMny6uvvmrmKl6zZo15XE+gq1y5st+rF1o2Vp5BFpTtWf/bRXLK47cW3J9z/iD5rn2XqGWqblwrZ993ZcH9L7pcKF90vSij6wkAHoIxQnOi3bhx48zJcieeeKI0bNiw4DZlyhS/Vy20ghLgXMH2TA8aSSA8qEcRilCswydi3XTuYmQeFU9qsT2Dj3AN2IH6FM6HYgQHFU5qsT0BILWoV1EahGIkhYomtdieAPzkch3k8ndDyE+0g/+oYMK7LRstXSg5edGXSc/evk2CjiEPQLhx8h2SQSiGkyEuiGzehs0+/8DcAMA2BGMkilAMp8Ocn9huABAMBGMkglCMYhHsEufytkpknuKgoQEEEIlgjJIQihHKkJfs5Xxd3hZcsQ1AWBCMURxCMeDY5U5T8v1nVxN5XKyWjsZP3y9MZQLhELagGLbvi8QxJRsQB+HHfunYhzSmgP2o3xELoTgOGr5wo8J0B8EYKFkY67wwfmcUj1AMFEJFKTLj+bni2j5N9X4lGAP2o75HJMYUAxGoIP+3DQ5qLZMeeL3Y5fJq1y9xGdfHEjLGGC4J61jbsH5vFEVPMfA/hJtwbAN6jAGEse5LlbFjx0qzZs2kUqVK0rFjR5k3b17cZZ955hnJysqKuunrIg0fPlxatmwpVatWlX322Ue6dOkiH3/8cdQyZ555puy///7mtQ0bNpRLL71UVq9eXfD80qVL5aSTTpL69eubZQ488ED5+9//Lrt27UrquxGKASrE0G0DgjEQW5jqgUwMs3LNlClTJDc3V4YNGyYLFy6UNm3aSLdu3WTdunVxX1OjRg356aefCm4//PBD1PN/+tOf5NFHH5VFixbJhx9+aAJ3165dZf369QXLaOB94YUXTPh98cUX5dtvv5Vzzz234PmKFStK79695Z133jHLjBkzRp544gmznskgFCP0qATDuQ0IxgBiCWN9mKjRo0dLv379pE+fPtKqVSsZP368VKlSRZ566qm4r9He4QYNGhTctDc30kUXXWR6h7V399BDDzWfsWXLFvniiy8Klrnhhhvk6KOPlqZNm0rnzp1l8ODB8tFHHxX0BOtrdZ00pOsy2rN88cUXywcffJDU9yMUI9So/MK9DQjGQDDpseTn8RSmelEDaORtx44dMZfbuXOnLFiwwARYT7ly5cz9uXPjn5y9bds2E1SbNGkiZ511lnz55Zdxl9XPmDBhgtSsWdME3Fg2btwozz33nAnH2kMcy7Jly2T69OlywgknSDIIxQitMFV6iI9gDARL5DHkZzgOShvxc15VWb+tWspv+r5Kw6qG0Jr/u40YMSLmemzYsEF2795dpKdX769Zsybma1q0aGF6kV999VWZNGmS7Nmzx4TZH3/8MWq5N954Q6pVq2bGAz/44IMyY8YMqVu3btQyt9xyixl3XKdOHVmxYoV5z8L0vfU9Dj74YDnuuOPkzjvvTGpbM/sEAKQYs1KEV7wAR3koG46p9Fm5cqUZ9+vJycmRVOnUqZO5RYbWQw45RB5//HG56667osYMf/bZZyZ461jg888/35xsV69evYJlbrrpJunbt68Zk3zHHXeYMcQapnV4RuSY561bt8rnn39ulh81apTcfPPNCa8voRgxUfnYK5leFfZzerBdw4VfB9I/RRrHVPpoII4MxfFoz2358uVl7dq1UY/rfR0rnAgd7nDEEUeY4Q2RtAe4efPm5qZjh7Wn98knn5QhQ4ZEfb7e9MQ8Ddbaw63jiiNDtz6mdLyz9mpfddVVcuONN5r1TgTDJ1AElY+9aJz932YcP0DqMCNEcGRnZ0u7du1k5syZBY/pcAi9HxlMi6NBVWeZ0GnViqPvG29ss/e8KmkZPRHPWzYR9BQjCpWPvQjEySMQIxW4+EN6cDwFT25urlx22WXSvn176dChg5n6LC8vz8z8oHRIw3777VcwLlnH9GrPr/YAb9q0SUaOHGmGP1x55ZXmeX3tPffcY2aL0KCswyd0HuRVq1bJeeedZ5bRYRTz58+XY4891sxjrNOx3X777XLQQQcVhHE98U57oVu3bm2Gf3zyySeml7lXr15xT8aLhVAcQ1grNyoge4W1zAZpe3H8AKnD8RRMvXr1MvMHDx061Jxc17ZtWzPLg3fynZ4ApzNSeH755RczhZsuq4FWe5rnzJljhjcoHdawZMkSmThxognEehLdUUcdZaZS0+nZlE759tJLL5k5hzVEa3ju3r27uTiHN/65QoUKcv/998vXX38t+fn5ZraLAQMGmKnckkEohkEFZC8CcXLoHQaA0hswYIC5xTJ79uyo+zqThN7i0ZkiNPAWR3t/Z82aVWJY11tZMaYYNOoWIxAnh0CMdKEsAPYjFIccFbm9CMTJIRADAIpDKA4xGnV7EYiTQyAGAJSEUBxSNOr2IhAnh0CMTKFsAHbjRLuQodK2G4E4OQRiAECi6CkGLEHAg6vl2rsBQUM9GS6E4kKomBFEBGKEgQv1L8eWe9in4UEoBgKOQIwwcSEYwz3Um+FAKAYCjECMMJZv24Mxx5mb2K/uIxQDAUUgRpgxzhhBRD3qNkIxEEAEYmAvW4Mxx5y72LfuIhQ7UPnCLQRihEEy5Zy6GUFDveomQnHI0LgEG4EYcKfu4vhzG/vXPYRiICAIxIB7wRhuo551C6EYCAACMeBmMOZYdN+Uh6b7vQpIES7zHELJNCpU6HY18uwvhKHM6+sp6wBSLZShuNeg7lKjRg2/VwMgEKcRve/B5W3Lsuwj77U27Bc/1jHWtrVhWwF+YvhEHFQesAnlNRqB2I55gcOwXYPEtqEnQKYRiotBhQ0byhjlNBqBuGzbjGDsNoIxEB+hGHHRWGVuO5d2W7OPohGIU7PNCMZuIxgDsRGKERONVLC3eVmCtKsIxKndZgRjtxGMgaIIxSUIY0Udxu8cFImEXfZPUQTi9GwzgrHbCMZANEIxotAoBTscs3+KIhCn94Q6gpPb2L/AHwjFgCXh2PWwVhoE4sxsr0wGJ9f3QRARjIG9CMUJCEslHZbvaSP2TVEE4sxfMAPuYv8ChGKEJBDALQRif7ZXpoKT6/sjqAjGCDtCMWiAYBUCsb/bi+Bkh9KWa/YvwiyUl3kGAJRessEpXX94pON99bu5/ocSgNjoKQ45Kn8A6WZD72PkrBx+XPIagP8IxQmgcgSAsgnyNHDxPo+6HwgXQnGI0UsMIJOCFjIT6RGm1xgID0JxSBGIAfghKAEz2fUgHAPuIxQDADIqFeGyLH/Yl+XzbQrGdH4AySEUO1QBJoqKEkAY69ZU9fbSawy4iVAMAPBFrGCZrrBJiAVQEuYpDhF6iAEETSLzAnMhCgCZQCgOSYVKIAZgQzAua13lUr0NILMIxSGoWAnEAIKutFeSc6muTgfdpoluI9oKhB1jih1HJQeEF4ExvVzavrQVAKEYAJzkUmBD+gKvPkcgBvZi+ITDqOiAcCIQoyS0D0BR9BQ7igoPCCcCMWKJPJGR9gGIjVDsICo8IJwIxCgObQNQPEKxY6j0gHAiEAPIhLFjx0qzZs2kUqVK0rFjR5k3b15Cr3v++eclKytLevbsGXeZa665xiwzZsyYmM/v2LFD2rZta5b57LPPYi6zbNkyqV69utSqVUuSRSgGAMsRiP3DtkeYTJkyRXJzc2XYsGGycOFCadOmjXTr1k3WrVtX7OuWL18uf/vb3+S4446Lu8zLL78sH330kTRq1CjuMjfffHOxz+/atUsuvPDCYj+nOIRih9BLDIQPoQxApowePVr69esnffr0kVatWsn48eOlSpUq8tRTT8V9ze7du+Xiiy+WO+64Qw488MCYy6xatUoGDhwozz33nFSsWDHmMm+//ba88847MmrUqLif9fe//11atmwp559/fim+HaHYGQRiIHwIxADKasuWLVG3HTt2xFxu586dsmDBAunSpUvBY+XKlTP3586dG/f977zzTqlXr5707ds35vN79uyRSy+9VG666SY59NBDYy6zdu1aE8afffZZE8JjmTVrlkydOtUM7ygtpmQDAAsRiIFw2LOtouzJr5j6983b+55NmjSJenzYsGEyfPjwIstv2LDB9PrWr18/6nG9v2TJkpif8eGHH8qTTz4Zd/yvuv/++6VChQpy3XXXxXw+Pz9fLr/8cjPeuH379mYoRmE///yzWWbSpElSo0YNKS1CMYBQK+3lhf1EIAaQKitXrowKkjk5OSl5361bt5oe4CeeeELq1q0bcxnteX7ooYfM+GQ9eS6WRx55xLzXkCFD4n6W9iJfdNFFcvzxx5dpnRk+ASD0NGR6N/jH1u1v63oDSgNx5C0nTijWYFu+fHkzlCGS3m/QoEGR5b/99lvTq3vGGWeYnmC9/fOf/5TXXnvN/L8+/8EHH5iT9Pbff/+CZX744Qe58cYbzQwX3rAIHZ6h66XPN2/e3DyuvcaXXXZZwTI61th7Dx2qsXnzZvP/xY13LoyeYgCIEXBs6z12hY0990AYZGdnS7t27WTmzJkF06rpeGC9P2DAgCLL6wlvixYtKnIinPb6au+wDtvQnuTIMcpKZ7PQx/VkPvXwww/L3XffXfD86tWrzTI6E4ZOCac0NOvQDs+rr75qhmXMmTNH9ttvv4S/I6EYAEro/SOkZZaNwdjGdQaSlZuba3pntZe2Q4cOZj7hvLy8ggDbu3dvE0JHjBhh5jE+7LDDol7vzR3sPV6nTh1zi6SzT2jPc4sWLcx97UWOVK1aNfPvQQcdJI0bNzb/f8ghh0Qt88knn5iTAAt/fkkIxQBQAnqPM4+QCQRPr169ZP369TJ06FBZs2aNuZDG9OnTC06+W7FihQmjtiIUAwACybZgbNv6AqUxYMCAmMMl1OzZs4t97TPPPFPi+8eaXSKSjjXWGSmKozNR6C1Z9sZ5AIDzijuJjRPcAKQSoRgAEGg2hV+b1hVANEIxACDwbAqbNq0rgD8QigEAToZNHd/r3QCgJJxoBwCw8mS2IIddTroD7ENPMQDAKgxPAJAOhGIAgHVsCMY2rCOAPxCKAQBWInQCSCVCMQDAWgRjAKlCKAYAWC3IwTjI6wYgGqEYAGA9wieAsiIUAwCcENRgrOsV1HUD8AdCMQDAGUEOn144DvI6AmFGKAaASGvXitx1l8gJJ4jUry+SnS1StarIoYfK0VMflkZLPhHJz/d7LVGMYkPn55+L/PWvIq1bi9SqtXf/6n7+859FRo0S2bw5Y+tIQAaChSvaAYDnscdEbrxR5Lffoh/ftUvkv/+V5vJfaT5/hrw8+P9JXu36fq0lSnNFud9/37tvH3646MLr1u29/fvfIvffL/LccyJdu2Z0XT1cBQ/wD6EYANQDD4jccssf98uXF+nRQ6RdO5GsLJFly2T7S69J5W2b/FxLlDYYDxwoMn78H082aiRy/vkideuKLFokMm2ayO7dIhs2iJxxhsisWSLHHOPLOivCMZB5hGIA+O9/RW699Y/79eqJTJ8ucsQRUYu9dMgrctAnM+X37JzMryNKHzLnzBG5OSIQH3nk3l7hGjX+eExD8CmniOzZI7Jzp8jVV4t88YVIOUYZAmFh3dE+duxYadasmVSqVEk6duwo8+bN83uVANjukUf29hJ6xo0rEohVfvkKsqxjN9lRrVZm1w9lM2FC0V8FIgOx0jHFvXr9cf/LL0U++CAz6wfAzlB82WWXyfvvvy9+mDJliuTm5sqwYcNk4cKF0qZNG+nWrZus07FgAFBaM2f+8f/77CPSs6efa4NUiwy3un9PPjn2cpGhuPDrADgv6VC8efNm6dKlixx88MFy7733yqpVqyRTRo8eLf369ZM+ffpIq1atZPz48VKlShV56qmnMrYOABwUWY/96U/8ZO6an3764/+bNo2/XOHnIl8HwHlJ1/yvvPKKCcL9+/c3Pbc6lOHUU0+VadOmyS49QztNdu7cKQsWLDCB3FOuXDlzf+7cuTFfs2PHDtmyZUvUDQAAACisVN0h++67rxnG8Pnnn8vHH38szZs3l0svvVQaNWokN9xwg3zzzTeSahs2bJDdu3dLfZ1PMoLeX7NmTczXjBgxQmrWrFlwa9KkScrXC4AD9tvvj///+mvmIXZNw4Z//P+KFfGX++GH+K8D4Lwy/Ub4008/yYwZM8ytfPnyctppp8miRYvM0IYHH3xQ/DZkyBAz3MO7rVy50u9VAhBEkWNMf/lF5NVX/VwbpNpxx/3x/xs37p1pIpYXXoj/OgDOSzoU6xCJF198UU4//XRp2rSpTJ06Va6//npZvXq1TJw4Ud5991154YUX5M4770zpitatW9cE77V6takIer9BgwYxX5OTkyM1atSIugFAEQMG7J2X2NO//94rnxWStft3af7xvySHuYrtctVV0fd1PuqtW6Mfmz1bz+b+436rVoRiIGSSnqe4YcOGsmfPHrnwwgvNdGht27YtssxJJ50ktfTymSmUnZ0t7dq1k5kzZ0rP/50Zruuh9wdogwYApXXooXsv7ezNVaxDstq3Fzn99L1Ts/3v4h3n/O/iHT8d3FZ2+L3OSFznznvnHX788b33P/lE5JBDYl+8Q+mln3UaN064BEIl6VCswyLOO+88M09wPBqIv//+e0k1HcesU8K1b99eOnToIGPGjJG8vDwzGwUAlMmQISJVq4rcfLOepbv3ssCvvLL39j+VfV1BlMmjj4pUrLj3X2/GkVjD/OrUEZk82Zer2QGwLBTrCXV+6dWrl6xfv16GDh1qTq7TXurp06cXOfkOAErluutEzjtvby/hu++KLF26d4yx9hw2bSpLazWTHw4/TvL2qef3miJB3uWSzZXt9CItffvu3b/vvSei55ls37537mL9teC00/YOtahZMxDrDCCzrLvMsw6VYLgEgLTRGQeGDdt7K2Q+YcUqMcOlDvl77DEJGoIw4D/rQjEAAC6ETBvWEQgTQjEAwClBDptBXjcg7AjFAABnBDV0BnW9APyB+WYAAE4IavAM6noBiEYoBgBYj+AJoKwIxQAAqwU5EAd53QBEIxQDAKxF6ASQKpxoBwAxELaCL+j7KOjrByAaoRhAqBFc7MR+A5BqDJ8AAFh76WZz+eYAIrQD9qGnGABgjVhhM6jBGIBd6CkGAFjBlt5XW9YTQDRCMQAg8AiaANKNUAwAsDIQBzEoB3GdACSGUAwACCxCJoBMIRQDAALJtkBs2/oCiEYoBgAEjm0B07b1BUpr7Nix0qxZM6lUqZJ07NhR5s2bF3fZl156Sdq3by+1atWSqlWrStu2beXZZ5+NWubyyy+XrKysqFv37t2jltHPK7zMfffdF7VMfn6+jBo1Sv70pz9JTk6O7LfffnLPPfck9d2Ykg0AEBiESyC4pkyZIrm5uTJ+/HgTiMeMGSPdunWTpUuXSr169YosX7t2bbntttukZcuWkp2dLW+88Yb06dPHLKuv82gIfvrppwvua6gt7M4775R+/foV3K9evXrU84MGDZJ33nnHBOPWrVvLxo0bzS0ZhGIAQCDYGohtXW8gWaNHjzbBVIOt0nD85ptvylNPPSWDBw8usvyJJ55YJLhOnDhRPvzww6hQrCG4QYMGxX62huB4y3z11Vcybtw4Wbx4sbRo0cI8dsABByT9/Rg+AQCWIYQBSJUtW7ZE3Xbs2BFzuZ07d8qCBQukS5cuBY+VK1fO3J87d26Jn6PDG2bOnGl6lY8//vio52bPnm16jzXQ9u/fX37++ecir9fhEnXq1JEjjjhCRo4cKb///nvBc6+//roceOCBpidaw7AOt7jyyivpKQaAsARjruQGuK/CliypuDsr5e/7e97e92zSpEnU48OGDZPhw4cXWX7Dhg2ye/duqV+/ftTjen/JkiVxP2fz5s1mfK+G7fLly8tjjz0mp5xyStTQiXPOOceE2W+//VZuvfVWOfXUU03Q1uXVddddJ0ceeaQZjjFnzhwZMmSI/PTTT6bnWn333Xfyww8/yNSpU+Wf//ynWc8bbrhBzj33XJk1a1bC24RQDACWIhgDKKuVK1dKjRo1ih3PWxY67OGzzz6Tbdu2mZ5iHZOsvbre0IoLLrigYFkdC3z44YfLQQcdZHqPTz75ZPO4vsajz+v45KuvvlpGjBhh1nfPnj0mdGsg1hPt1JNPPint2rUzPdPekIqSMHzCETSMQDgxlAJAWWggjrzlxAnFdevWNT23a9eujXpc7xc3HliHWDRv3tzMPHHjjTea3lsNs/FoYNbPWrZsWdxl9CQ/HT6xfPlyc79hw4ZSoUKFgkCsDjnkEPPvihUrivn2hdY14SUReARjIJwIxv5h2yMssrOzTc+r9vZ6tIdW73fq1Cnh9/F6deP58ccfzZhiDbrxaM+zhm1vxotjjjnGhGQdfuH5+uuvzb9NmzZNeN0IxY4hGAPhRDgDkG65ubnyxBNPmBkkdMYHPSkuLy+vYDaK3r17m/G+Hu0RnjFjhhnzq8v/4x//MPMUX3LJJeZ5HVJx0003yUcffWR6fTVgn3XWWaZn2ZudQscW69Rvn3/+uXmf5557zowX1vfYZ599zDJ6sp+OOb7iiivk008/NScE6vAKHbsc2XtcEsYUOxqMaSCB8GGMMYB06tWrl6xfv16GDh0qa9asMUMipk+fXnDynQ5V0B5cjwbmv/71r6b3t3Llyma+4kmTJpn3UToc44svvjAhe9OmTdKoUSPp2rWr3HXXXQXDOPTf559/3pz8pz3MekKehuLIccb6mToDxcCBA83MFnqhED1ZT0N4MrLydY6MkNCpRmrWrGnOhIwcVO4qgjFclOrQ5+JxYmswLs2+8Pu72lJ+IreTLetsi+2//SpXDe2V8mzhZZaRL82QylWrSqptz8uTm845JTSZKBEMnwAAxxB6UFJA9vuPCSCICMUOo9IDwotgjJIQjoFohGLHUeEB4UUwTi9btm9J7QDhGNiLUGx5ZZcIKjsArta5+jqX6ms/EY4Rdsw+ERLMSAEgqCLrptKeEFZ4WcLdXmwHIHGE4hBNb0QwBhA0xdVJserfROuweEEbAOJh+ETI0DgACIpYAbek0Fuan/gZYgEgEYTiECIYA/CbHyE1lZ9JyAbcQyguARUfALhTr5a115heZ8BdhOKQorcYgB9SFSjLWoeVZj1sC8PU80ByCMUAgIwIWqhM5qS9oK07gNQjFCfA1cqQXgQAmVKW6dXSqaTA62r9D6AoQnHIEYwBpJsNwbLwOtI7DIQP8xQDANImncEy1XOvE4KBcCMUg4t6INQX1XG9/Kd6e7m8rVxS2n3O/kWYMXwCBsMoYJNUN9yul/9UbS8Ck9vYvwg7QjEAKxGMM7u9MhmYXN8XQUQgBgjFCQlLBR2W72kb77K27J+iCMaZ2V4EJrexf4G9CMUIVSiwSawgzP4pimCc3u2V6cDk+vYPGgIx8AdCMYqgUfJXSb3C7J+iCMbpu2hFJrm+3YOGQAxEIxSXgEoamZLMEAnKZVEE49RuLwKx2wjEQFGEYsREAxX87c0+KopgnJrtRSB2G4EYiI1QjLhoqNKvrCfQcQJeUQTjsl/JLZNc375BQyAG4iMUF4PKGragrEYjGJduexGI3UYgBooXyivaTXloulSuVKXYCprKGrZx/cpsyeLKd8khELvF5bIKpEsoQ3HYUVkGB8HNLmzf1ChrmbdtH1BuADswfALwGT/124Xt62+ZJ1wiaHoN6u73KiBFCMVAABCM7cL29YeNgZiy4jb2r1sIxUBAEIztwvbNLBsDMdxGHeAeQnHIKt0wfEebEYztwvbNTFm3td6ifLiLfesmQjEQMARju7B908vWQAx3ccy7i1AMBBDB2C5s3/SwORBTJtzEfnUboRgIKIKxXdi+qWVzIIabOMbdRygGAoxgbBe2b2rKt+2BmHLgHvZpOBCKHauM4R6CsV3YvmVDHYyg4ZgOD0IxYAGCsV3YvsmV7cib7dj3bmF/hguhOGQ4wO1FME6eCyELAJAZhOIQCkMYchUhL3lsM2QKdStgN0JxSFF524uQlzy2GQCgJITiECMY24uQlzy2GQCgOITikCMY28uVE5Myie2FdKEuBexHKAaVueUIeslhewEAYiEUxxDGRpNgbLcwltmyYHsBAAojFKMAwdhuBL3ksL2AYLdHtEnINEIxolAJ2Y2glxy2F1KBejN92LbBM3bsWGnWrJlUqlRJOnbsKPPmzYu77BNPPCHHHXec7LPPPubWpUuXmMt/9dVXcuaZZ0rNmjWlatWqctRRR8mKFSuilpk7d678+c9/Ns/XqFFDjj/+eNm+fbt5bvbs2ZKVlRXzNn/+/IS/G6EYRVAJ2Y2glxy2FxDsNog2KTimTJkiubm5MmzYMFm4cKG0adNGunXrJuvWrYu5vIbVCy+8UP7973+bUNukSRPp2rWrrFq1qmCZb7/9Vo499lhp2bKlWf6LL76Q22+/3YRuj762e/fu5rUaqjXoDhgwQMqV2xtjO3fuLD/99FPU7corr5QDDjhA2rdvn/D3y8rPz8+XkNiyZYv5K2TCnVOkcqUqxS7LQUhYQHik43jn+AFSezz6cUwlUjd42WLz5s2mBzNVvPcd+dIMqVy1qqTa9rw8uemcU5Ja744dO5pe3EcffdTc37Nnjwm6AwcOlMGDB5f4+t27d5seY3197969zWMXXHCBVKxYUZ599tm4rzv66KPllFNOkbvuuiuh9dy1a5fst99+Zr00YCeKnmIASDECMZC648ev6SfD0jmm4TvytmPHjpjL7dy5UxYsWGCGQHi0p1bva09uIn799VcTWGvXrl0Qqt98803505/+ZHqc69WrZ4L3K6+8UvAa7YX++OOPzXPaI1y/fn054YQT5MMPP4z7Oa+99pr8/PPP0qdPnyS2BKEYIRaWCq8kYd4O6TiZh0AM2D8Xe5Dqxext+ZK9NQ23bXsHCmhPr/ZI1/zfbcSIETHXY8OGDaanV0NpJL2/Zs2ahL7LLbfcIo0aNSoI1hp4t23bJvfdd58ZHvHOO+/I2WefLeecc4689957ZpnvvvvO/Dt8+HDp16+fTJ8+XY488kg5+eST5Ztvvon5OU8++aQJ2Y0bN05iS4tUSGppwDFa8RFiwrkdGDIBBJefx1KQAnEmrFy5Mmr4RE5OTlo+R4Pv888/b8YNe+OFtadYnXXWWXLDDTeY/2/btq3MmTNHxo8fb3qEvWWuvvrqgp7fI444QmbOnClPPfVUkRD/448/yr/+9S954YUXkl5HeooRemGrAOMJ03YgEAOxhakeiCWM318DceQtJ04orlu3rpQvX17Wrl0b9bjeb9CgQbGfMWrUKBOKtSf48MMPj3rPChUqSKtWraKWP+SQQwpmn2jYsKH5t7hlIj399NNSp04dM5tFsgjFQEgrwuK2Q/1vF8klN59RcDvwk3eLLFt149qoZQ5/Z7LYgEAMIBbageJlZ2dLu3btTA+tR3tx9X6nTp3ivu6BBx4wJ8jpsIfCM0Hoe+qJe0uXLo16/Ouvv5amTZua/9fp33TIRXHLeHTuCA3FehKfnryXLEIx8D9UiO5vBwIxEM5jvyRh/u7JyM3NNXMPT5w40cwt3L9/f8nLyysY1qBhdMiQIQXL33///Wb2Bx3moOFWxx7rTccRe2666SYz1Zu+77Jly8zMFK+//rr89a9/Nc/rXMO6zMMPPyzTpk0zy+h7LlmyRPr27Ru1frNmzZLvv//eTMdWGowpBkI+tjaWUy7oJPK4OIUT6gDEQiBOXK9evWT9+vUydOhQE251/K/2AHsn3+lwBm/uYDVu3Dgza8W5554b9T46z7GeOKf0xDodP6xjg6+77jpp0aKFvPjii2buYs/1118vv/32mxl3vHHjRjM/8owZM+Sggw4qcoKdzlChcx6XBqG4mAaPAyWcCMbuIRADxQtrexfW710WAwYMMLdY9CS6SMuXL0/oPa+44gpzK47Og1zSXMiTJ5dtGB/DJ4AYqCjdQSAGEAv1PAqjpxiIgx7jPzRaulBy8rZEPZa9/Y8xYUFFIAZKFsZwGMbvjJIRihHaUOjydyvrdtDZJ06JuN/s8w/MzSY0egBioW5APIRihDoYp6sSZZsBQPAQiFEcQjFKRDAuW8Xrwrabc/4g+a79H9e79+YpPvu+0k17kwmcLAsgEvUBSkIoRkIIxqmriNmOAJBZBGIkgtknkDAqldRtR7YlAGQG9S0SRShGUqhcUodtmX70ygPhRj2LZBCKkTQqmdRhWwLwk8t1kMvfDelBKEapUNmkDtvSDvQ6A/agXkVpcKIdSo2T79zdlmsPai2THni92GXyatcvcZkgYBYKIFw43lFa9BSjTKh8UodtCQBlQz0K53uKly9fLnfddZfMmjVL1qxZI40aNZJLLrlEbrvtNsnOzvZ79UIvaL2cNmL7Bbu3mP0DV7lUtgnECEUoXrJkiezZs0cef/xxad68uSxevFj69esneXl5MmrUKL9XDyHgUsMRVpH7MJHGk32O0ohXtihP6UUgRmhCcffu3c3Nc+CBB8rSpUtl3LhxhGIAKQnIhBbATgRihCoUx7J582apXbt2scvs2LHD3DxbtmzJwJqFE0MoYCvKLVKFcJZ5bHNI2E+0W7ZsmTzyyCNy9dVXF7vciBEjpGbNmgW3Jk2aZGwdAQBA+hCI4VQoHjx4sGRlZRV70/HEkVatWmWGUpx33nlmXHFxhgwZYnqUvdvKlSvT/I0AAGFEQMsstjecGz5x4403yuWXX17sMjp+2LN69Wo56aSTpHPnzjJhwoQS3z8nJ8fcAACA/QjDcDYU77vvvuaWCO0h1kDcrl07efrpp6VcOStHfgAAHENQywy2M9LNihPtNBCfeOKJ0rRpUzPbxPr16wuea9Cgga/rBgAA0otAjEywIhTPmDHDnFynt8aNG0c9l5+f79t6AQCA9CIQI1OsGIOg4441/Ma6AQDgFwJberF9kUlWhGLYgcoLAJAqtCnINEIxAAAIFAIx/EAoRkpRkQEAyoJ2BH4hFCPlqNAAAKVB+wE/EYqRFlRsABCOuj5V9T3tBvxGKEbaUMEBQDjq+LLW97QXCAJCMdKKig4AwlG3l7a+p51AUBCKkXZUeAAQjjo92fqe9gFBQihGRlDxAYD9EqnLEx1nTLuAoCEUI2OoAAHAXqnsBaY9QBARigEAQLFSOV6YQIygIhQXY9LIN/xeBedQGQKAXVI5swRtAIKsgoRQr0HdpUaNGgkFYv2Xgzi1dHvyBwcAhAftKGxAT3ECCHCpRwUJ+HuhBABANEJxggjGqUfjDpQuDHPsIFMoawgTQnESCMapR4ULlK5nmGMHAFKLUJxk+CUYpx6NO1D6eV45fgAgNQjFpUAwTj0adqD0IZfjBwDKjlBcSgTj1KNhR1iloseX4wepRplC2BCKy4BgnHr8HIywlfVUlneOHQAoPUJxGRGM04PGHS5KRxCO9RkAgOQRilOAYJweNO5wRaZ/AeHYyWz9791cEvYy5OI+TZWxY8dKs2bNpFKlStKxY0eZN29e3GW//PJL+ctf/mKWz8rKkjFjxhRZxnuu8O3aa681zy9fvjzm83qbOnVqwfvMnz9fTj75ZKlVq5bss88+0q1bN/n888+T+m6hvKIdALgeKpL9fAJA2bEN3cGVV2ObMmWK5Obmyvjx400g1pCr4XPp0qVSr169Isv/+uuvcuCBB8p5550nN9xwQ8z31DC7e/fugvuLFy+WU045xbxGNWnSRH766aeo10yYMEFGjhwpp556qrm/bds26d69u5x55pny2GOPye+//y7Dhg0z67Zy5UqpWLGiJIJQDAAgBDi+vfz+I81GbLOiRo8eLf369ZM+ffqY+xqO33zzTXnqqadk8ODBRZY/6qijzE3Fel7tu+++Uffvu+8+Oeigg+SEE04w98uXLy8NGjSIWubll1+W888/X6pVq2buL1myRDZu3Ch33nmnCdFKQ/Hhhx8uP/zwgzRv3lwSwfAJAIBBCHBzO9m2vsisLVu2RN127NgRc7mdO3fKggULpEuXLgWPlStXztyfO3duStZFP2PSpElyxRVXmOERseg6fPbZZ9K3b9+Cx1q0aCF16tSRJ5980rzH9u3bzf8fcsghZnhGougpBpD2xtemHrWws60HNNNsC5i2rS+Kytm0S3J27Er5++7Zvvc9vZ5Vj/awDh8+XArbsGGDGeZQv379qMf1vvbUpsIrr7wimzZtkssvvzzuMl7Y7dy5c8Fj1atXl9mzZ0vPnj3lrrvuMo8dfPDB8q9//UsqVEg86hKKAZRJoldei0ToCjaCsRsB07b1hT90zG2NGjUK7ufk5Pi2Lhp4dZxwo0aNYj6vPcCTJ0+W22+/vcjj2nN8zDHHyP/93/+Z8D5q1Cjp0aOHGbNcuXLlhD6fUAwgaam+0ISLAUy/k82hhGAczbZ9adv6wj8aiCNDcTx169Y143vXrl0b9bjeLzzmtzR07O+7774rL730Utxlpk2bZk7e6927d9TjGpR1lgodxqFDOrzHdBaKV199VS644IKE1oFQDPiIhit6OxDCgoVgbOdxatv6wg7Z2dnSrl07mTlzphmmoPbs2WPuDxgwoMzv//TTT5sZLLR3t7ieZJ1hovDJeRqUNQxHjkP27us6JooT7QCE6uIWSE7Y90Vx3z+I2yaI6wR35ObmyhNPPCETJ06Ur776Svr37y95eXkFs1FoD+6QIUMKlteT3vSkOL3p/69atcr8/7Jly6LeV4OrhuLLLrss7hhgfc37778vV155ZZHndAq3X375xcxtrOul8yPrOul7nXTSSQl/P0IxgEAiIAdHWPeBbd/btvX1A798lE2vXr3MWN2hQ4dK27ZtTcCdPn16wcl3K1asiJpTePXq1XLEEUeYmz6ur9X/LxxsddiEvlZnnYhHp31r3LixdO3atchzLVu2lNdff12++OIL6dSpkxx33HHms3XdGjZsmPD3y8rPz8+XkNCpRmrWrCmbN28udvxMaQ4aKqP0cLkCo8y4X0Zc28dB3tZ+7bugbBPXylo6RO6rVG6vRLNFad93/MNvSuXKVSXVtm/Pk2uu65Hy9bYZY4pTxPaTagBb2RZebBaWMca21eW2ra8fCpdb2mzEQihOIQ6y1AtLI4z0K+nYTEc5c7FOcO37pHN7lFSmUvEHHfujZLQhSBShOMVcbAT9RjBGJlDOkGqpagtoU0ov3jHNNkUsnGiXBjSsqUcFBgBIBm0xkkUoThMOxtRzLRhTRgAgPahfURqE4jTioEw914IxgoXyBdiPthelRShOMw7O1CO4AABitbe0uSgLQnEGcJCmHsEYAOChnUUqEIozhAM29QjGSAfKFWAX2lekCqE4gzhwU48AAwDhxRVokUqEYliPCg4AwoXxw0gHQnGGcRCnh63BmPIQTLaWJyAMqDeRLlzRDqENMkGpWLkKIgBkps6mrkVx6ClGaAWpcgxKQAeAICIQIxMIxQi1IFWSBGMAKIpAjEwhFCP0glRZEozdKxfsU6D0CMTIJEKxD2gkgydIlSblAwAIxMg8QjEQwMqTYAwgzAjE8AOhGAhoJUowBhA2qZp/OEh1OexBKPYJgSe4glSZUk4AhEEqL8YRpDocdiEUAwGvVAnGAFyWyjouSHU37EMoBiyoXAnG/ghSGQBck+pLNXO8oqwIxT4i6ARfkCpZygsAV6S6PgtSXQ17EYrjHFyZOsAIOkgG5cW/+iCT9QLgMgIxgopQXAyCMYJY4VJe/BW08gDYhECMICMUBwRBJ9iCVvFSXuwqD+wvgECM4CMUB+igo+EMtqBVwJQXfwWtPABhOaFOcfwhHQjFCSAYI6gVMeXFX0ErD0DQpKOO4rhDuhCKE0QwRlArZMqLv4JWHoCgIBDDNoTiJBCMEdSKmfLiL2amANI7XEJxjCHdCMVJIhgjqBU05cV/QSsTgCv1EMcWMoFQXAocnAhqWSAY+y9oZQLIFAIxbEcoLiXmMEZQK2zKjP+CViYAG4dLKI4lZBKhuAwIxghqxU2Z8R9XwkMYpLOu4bhBphGKLUHICb6gVeCUGQDpRCCGayr4vQKAS+JV5H4FVP1cGhcAqUYghosIxRYh4NjLz7BMuQGQSgRiuIrhE5bhJ3G3MC4dgE0IxHAZodhCBBy3EIwB2IBADNcRii1FwHELwRhAkBGIEQaEYosRcNxCMAYQRARihAWh2HIEHLcQjAGE4aIcikCMoCEUO4CA4xaCMQC/pbt+IBAjiJiSzRFMu+UW3ZdM1xYe9MYhSAjEKM7YsWNl5MiRsmbNGmnTpo088sgj0qFDh7jLT506VW6//XZZvny5HHzwwXL//ffLaaedVvD8tm3bZPDgwfLKK6/Izz//LAcccIBcd911cs011xQso5910003yYwZM2Tr1q3SokULue222+Qvf/lLwTILFy6UW265RebPny/ly5c3z40ePVqqVasmiaKn2KGDm54/t6SibBW+1HCsSw9TbtzG/kUyCMQozpQpUyQ3N1eGDRtmQqiG4m7dusm6detiLj9nzhy58MILpW/fvvLpp59Kz549zW3x4sUFy+j7TZ8+XSZNmiRfffWVXH/99TJgwAB57bXXCpbp3bu3LF261Dy2aNEiOeecc+T8888376lWr14tXbp0kebNm8vHH39s3u/LL7+Uyy+/XJKRlZ+fny8hsWXLFqlZs6Zs3rxZatSo4WzDQ6XjjmTLVbL7PvL9KTf+yFTdwf5FSQjEwcoW3vuOf/hNqVy5qqTa9u15cs11PZJa744dO8pRRx0ljz76qLm/Z88eadKkiQwcOND09hbWq1cvycvLkzfe+KNsHX300dK2bVsZP368uX/YYYeZ5bQ32dOuXTs59dRT5e677zb3tbd33LhxcumllxYsU6dOHdPrfOWVV8qECRPM63/66ScpV25vf6+G58MPP1y++eYbE5YTQU+xgwd7kAI6MlOuInt+S4ty4zb2L4pDIA4vDd+Rtx07dsRcbufOnbJgwQLTI+vRAKr3586dG/M1+njk8kp7liOX79y5s+kBXrVqlWg/7b///W/5+uuvpWvXrlHLaC/1xo0bTRB//vnn5bfffpMTTzzRPK/rnJ2dXRCIVeXKlc2/H374YcLbglDs6EFPA+i+wsMgUoFy47Z0zyYAOxGIg63cpm1S7petqb9t2mbeX3t6tUe65v9uI0aMiLkeGzZskN27d0v9+vWjHtf7OuY3Fn28pOV1THKrVq2kcePGJth2797djFs+/vjjC5Z54YUXZNeuXaZ3OCcnR66++mp5+eWXC3qA//znP5v31LHOGt5/+eWXgp5r7T1OeFsnvCSsO/hp/NwsU6kOwoVRbjLHr23NPoaHQIyVK1eaIRSb/3cbMmRIRj9fQ/FHH31keou1J/of//iHXHvttfLuu+8WLKNDIzZt2mQe++STT8w4ZB1TrEMk1KGHHioTJ040r61SpYo0aNDAnLCnATyy97gkzD5h8cwBiWB2ATdkeh9SbtzHPgaBGErHEycyprhu3bpmVoe1a9dGPa73NYTGoo8Xt/z27dvl1ltvNb2+PXr0MI/pOODPPvtMRo0aZYZefPvtt2YMs56cp+FX6Ql+H3zwgelR9sYmX3TRReam71+1alXJysoys08ceOCBCW8LeopDUBkEJaDDLpQb97GPw4tAjGRlZ2ebE+BmzpxZ8JiO79X7nTp1ivkafTxyeaXTqnnL65AIvRXuzdXwre+tfv31V/NvcctE0t5hPTFPxyBXqlRJTjnllIS/I6E4JJUCjR9Kg3Lj/rYNynogcwjEKK3c3Fx54oknzFAFnT6tf//+ZnaJPn36FEydFjn8YtCgQWZ6NB3WsGTJEhk+fLgZ/qBTrintoT7hhBPMHMSzZ8+W77//Xp555hn55z//KWeffbZZpmXLlmbssI4jnjdvnuk51vfTcK3Tu3m0N1mnidOT9LQHWT9Dx0fXqlUr4e9HKE6zIFUONH4oDcqN+9jH4UEgRln06tXLDGsYOnSomVZNhzlo6PVOpluxYkXUiW06a8TkyZPNlGk65GHatGnmIh06DZtHZ5LQad4uvvhic8LdfffdJ/fcc0/BxTsqVqwob731luy7775yxhlnmOEVGpo1mEdeBEQDs/YKt27d2nze448/bi4CkgzmKQ5ho0OlhUzMgQx76gQP+9htBGI75ymecOcUqVypiqTa9t9+lauG9vIlEwUVJ9plCCffwXaUm9QJSl0AOwWx/FA3wAWE4gwiGMPm8kh5sQ/7zE5BaSeAsCEUZxjBGLahjNiJ/RZ8QWkLyoqyBlcQin1AMIYNKBf2Yt8FT1Dq/FSinME1zD7hkyBVJi5W1nCnfCI57LtgcfXS2pQzuIhQ7KMgVSouVtpIXrovIY30Hmvsu2BxtV6lnMFVhGKfBalycbUCh31lEclj/wWLq/Up5QwuIxQHQJAqGVcrcsRH77D92H/B4mo9SjmD66wLxTt27DBXUcnKyjJXUnFFkCobVyt0BLvcoXTYh8Hiav1JOUMYWBeKb775ZmnUqJG4KEiVjqsVO+CSINUZcLfepJwhLKwKxW+//ba888475rrbrgpS5eNqBQ+4IEh1Rdi5OsOEopwhTKwJxWvXrpV+/frJs88+K1WqVEl4qIVeOzzyZoMgVUKuVvSAzYJUR8BdlDOEjRUX78jPz5fLL79crrnmGmnfvr0sX748odeNGDFC7rjjjrSvn+u4wAeQOhxL8BtlEAhgT/HgwYPNCXPF3ZYsWSKPPPKIbN26VYYMGZLU++vymzdvLritXLkybd/FdfQYA4C9daM3ywyBGAhoT/GNN95oeoCLc+CBB8qsWbNk7ty5kpOTE/Wc9hpffPHFMnHixJiv1eULvwalR4+xe9ingPs4xgELQvG+++5rbiV5+OGH5e677y64v3r1aunWrZtMmTJFOnbsmOa1RCRClHs9XJE9XexbwC0c04BjY4r333//qPvVqlUz/x500EHSuHFjn9YqvAjG7v7k6z3O/gXsx3EMODr7BILFlnF0KB2Xp5gCUiXIxwiBGAhJKG7WrJmZkUKvbAf/BLlBQGr2G+EYsA+BGAhRKEZwEJjCsb+8cMz+BoKNQAw4PqYYwcYY43DhxDyEXVD/OOR4BMqGnmI43UggvfuI3mOETRDLO/MPA6lBKIbTjQUyg30P+IMwDKQOoTigbK3oCEfhLVPsewCAzQjFAUYwRqoRjAF3yretbQQQVITigLO10gta44E/EIwBuzGGGEgPZp+wgFZ+NgYNZqVApJLKMGUFKBnHCZA+hGJLEIxhU3kqzXtziWmgeBwbQHoxfMIitlaINob5MAhqeaK8IIj8LpdBPV4BlxCKLWNrxeh3g4LixyYGrVwx/zHwh6Adn4CrCMUWsrWCJOTYE5CDUsYIxwi7oByLQBgQii1la0VJwLFHkEIy5QZh5PdxB4QNodhitlaYBBw7+V3eKDcIE7+PNyCMmH3CcrbOSgE7+V3ekv1sggVSIRVlnrIIBF8oe4qnPDRdXBKEn7dLO1aUMaNIJ8oWgsC2+hkIq9D2FLs4f67fvXhlEW+9XdtHtrOxfLl4rMMelD3AHqHsKba5gQ9bBUxPMlKBMgQ/uFYfA64LdSh2tbF0sSJ2cT8hsyhDyCQX62HAdaEPxa42li5WyC7uJ2QWZQiZ4GL9C4QBodjhxtLFitnF/WQLV7a9K98DmcGMJ0B4EIodbyxdrKBd3E/ILMoQ0sHF+hYIE0JxCBpLFytqF/cTMosyhFRysZ4FwoZQHJLG0sUK28X9FEQuzwDi6vdCZrlYvwJhFNp5isM4jszmeYzhjzCUF+YxTi/Xt6/L3w0IG3qKQxYeqMCRjLCUF1uOXxvZXoa8K4bGuwFhM3bsWGnWrJlUqlRJOnbsKPPmzSt2+alTp0rLli3N8q1bt5a33nor6vn8/HwZOnSoNGzYUCpXrixdunSRb775JmqZjRs3ysUXXyw1atSQWrVqSd++fWXbtm1Ry3zxxRdy3HHHmc9p0qSJPPDAA0l/N0JxCBtWVypyV74HgsGW4xcA/DJlyhTJzc2VYcOGycKFC6VNmzbSrVs3WbduXczl58yZIxdeeKEJsZ9++qn07NnT3BYvXlywjIbXhx9+WMaPHy8ff/yxVK1a1bznb7/9VrCMBuIvv/xSZsyYIW+88Ya8//77ctVVVxU8v2XLFunatas0bdpUFixYICNHjpThw4fLhAkTkvp+Wfka0UNCN1rNmjVlwp1TpHKlKhL2sGZ7CLBlO7vA9rKSDMoVgNJki82bN5ueTFsyy/bffpWrhvZKar07duwoRx11lDz66KPm/p49e0yv7MCBA2Xw4MFFlu/Vq5fk5eWZIOs5+uijpW3btiYEawRt1KiR3HjjjfK3v/3NPK/rU79+fXnmmWfkggsukK+++kpatWol8+fPl/bt25tlpk+fLqeddpr8+OOP5vXjxo2T2267TdasWSPZ2dlmGV2fV155RZYsWZLwNqGnOMQBgsYfsPf4Req5fFIp7KXhNV03L3xH3nbs2BFzPXbu3Gl6YXV4g6dcuXLm/ty5c2O+Rh+PXF5pL7C3/Pfff2+CbOQy+oeAhm9vGf1Xh0x4gVjp8vrZ2rPsLXP88ccXBGLvc5YuXSq//PJLwtuaE+1CfnIJJ98hEWErJ7Ycv0idMJVv2EEDXoMGDWTQvX3S9hnVqlUzPb2Rhg0bZoYeFLZhwwbZvXu36cWNpPfj9cZq4I21vD7uPe89Vtwy9erVi3q+QoUKUrt27ahlDjjggCLv4T23zz77FLsdCt43oaXgdMNqY+CxYbu6xtvmtpWV0vK+J2UNgB/0hDHtSdUe2nTR4QtZWVlRj+Xk5EhYEYrTjGAM1xQuz66XG1uO4UyI3NcubRPXyzDsDsZ6C4K6detK+fLlZe3atVGP633t0Y5FHy9uee9ffUxnn4hcRscde8sUPpHv999/NzNSRL5PrM+J/IxEMKY4A2ypcF1q5JA5YZiiypZjOF1cHmvr6vcC0jGco127djJz5syCx/REO73fqVOnmK/RxyOXVzqDhLe8DnnQ0Bq5jI5r1rHC3jL676ZNm8x4Zs+sWbPMZ+vYY28ZnZFi165dUZ/TokWLhIdOKEJxhthS8doQamxYxzBzNRzbcgynksthWLn83YB0yM3NlSeeeEImTpxoZoXo37+/mV2iT5+945579+4tQ4YMKVh+0KBBZqaIf/zjH2bcsY5V/uSTT2TAgAHmeR26cf3118vdd98tr732mixatMi8h84ooVO3qUMOOUS6d+8u/fr1M3Mi/+c//zGv15kpdDl10UUXmdCuU7/p1G06ddxDDz1k1jcZhOIMsqUCdjHQIPNcLEe2HMNl5XoYBlA6vXr1klGjRpmLbejwhs8++8yEXu+kthUrVshPP/1UsHznzp1l8uTJZr5gndN42rRpZpq0ww47rGCZm2++2UzppvMO63RvelEOfc/IYSPPPfecuQDIySefbKZiO/bYY6PmINYZK9555x0zBlt7s3WKN13HyLmME8E8xQEStBAR1EYxaNsJdpajIApS2S5pvwVpXVP9HV34brB/nmJkHqHYAn5W0EELNDRWdgpaOQoyyrg/XD2JEOlHKHYHwycs4GegoHFAKlCOEscfEP5w/WRRACUjFFuCYAzbUY4SRzAGgMwjFFsk7ME4COuAsmEfJo5gDACZRSi2jN/BmFADZA7BGAAyh1BsIb8bSoIxyoLyY9fxDgBhEcpQ3GtQd7Gd3w1lpoMNQcot7E+7jncACINQhmJXGmW/G0oXtiH8Q/mx63gHANeFNhS70ij73VAyzhhlQdmx63gHAJeF6uIdOrF2rVq1ZOXKlVETbE95aLrYLghDQtK1HYPw3ZB+LhyHmcDxAATv4h1NmjSRTZs2mYt4wF6hCsU//vijKbgAAACppB1ujRs39ns1UAahCsV79uyR1atXS/Xq1SUrK0tc+gu1cO83god9ZQ/2lV3YX/ZwcV9pjNq6das0atRIypUL9ahU61WQENHC6upfcVq5uFLBuI59ZQ/2lV3YX/ZwbV8xbMIN/EkDAACA0CMUAwAAIPQIxZbLycmRYcOGmX8RbOwre7Cv7ML+sgf7CkEWqhPtAAAAgFjoKQYAAEDoEYoBAAAQeoRiAAAAhB6hGAAAAKFHKHbQjh07pG3btuaqfZ999pnfq4MYli9fLn379pUDDjhAKleuLAcddJA5I3vnzp1+rxpEZOzYsdKsWTOpVKmSdOzYUebNm+f3KqGQESNGyFFHHWWuUFqvXj3p2bOnLF261O/VQgLuu+8+0z5df/31fq8KEIVQ7KCbb77ZXG4SwbVkyRJz2fHHH39cvvzyS3nwwQdl/Pjxcuutt/q9aqE3ZcoUyc3NNX+kLFy4UNq0aSPdunWTdevW+b1qiPDee+/JtddeKx999JHMmDFDdu3aJV27dpW8vDy/Vw3FmD9/vqn3Dj/8cL9XBSiCKdkc8/bbb5sG/cUXX5RDDz1UPv30U9NrjOAbOXKkjBs3Tr777ju/VyXUtGdYeyAfffRRc1//eGnSpIkMHDhQBg8e7PfqIY7169ebHmMNy8cff7zfq4MYtm3bJkceeaQ89thjcvfdd5u2acyYMX6vFlCAnmKHrF27Vvr16yfPPvusVKlSxe/VQZI2b94stWvX9ns1Qk2HryxYsEC6dOlS8Fi5cuXM/blz5/q6bij5+FEcQ8GlPfs9evSIOr6AIKng9wogNbTD//LLL5drrrlG2rdvb8aswh7Lli2TRx55REaNGuX3qoTahg0bZPfu3VK/fv2ox/W+DnlBMGlvvo5PPeaYY+Swww7ze3UQw/PPP2+GI+nwCSCo6CkOOP25Vk9IKO6mjbUGqq1bt8qQIUP8XuVQS3R/RVq1apV0795dzjvvPNPTDyD5HsjFixeb4IXgWblypQwaNEiee+45c/IqEFSMKbZgnNzPP/9c7DIHHnignH/++fL666+b0OXRHq/y5cvLxRdfLBMnTszA2iLR/ZWdnW3+f/Xq1XLiiSfK0UcfLc8884z5qR7+Dp/QoUfTpk0zsxl4LrvsMtm0aZO8+uqrvq4fihowYIDZL++//76ZzQXB88orr8jZZ59t2qPI9knbK63zdMakyOcAvxCKHbFixQrZsmVLwX0NW3rGvDbueuJQ48aNfV0/FKU9xCeddJK0a9dOJk2aRKMQEHq8dOjQwfz64v00v//++5vwxYl2waFNl578+PLLL8vs2bPl4IMP9nuVEIf+ivnDDz9EPdanTx9p2bKl3HLLLQx5QWAwptgR2mhHqlatmvlX578lEAczEGsPcdOmTc04Yu1h9jRo0MDXdQs7nb1Fe4Z1bL6GYz07Xqf50kYcwRoyMXnyZNNLrHMVr1mzxjxes2ZNM/c3gkP3T+HgW7VqValTpw6BGIFCKAZ8oPOq6sl1eiv8Rws/3virV69e5o+UoUOHmqCl00ZNnz69yMl38JdOX6j0j8tITz/9tDnpGACSxfAJAAAAhB5n9QAAACD0CMUAAAAIPUIxAAAAQo9QDAAAgNAjFAMAACD0CMUAAAAIPUIxAAAAQo9QDAAAgNAjFAMAACD0CMUAAAAIPUIxAAAAQo9QDAD/s379emnQoIHce++9BY/NmTNHsrOzZebMmb6uGwAgvbLy8/Pz0/wZAGCNt956S3r27GnCcIsWLaRt27Zy1llnyejRo/1eNQBAGhGKAaCQa6+9Vt59911p3769LFq0SObPny85OTl+rxYAII0IxQBQyPbt2+Wwww6TlStXyoIFC6R169Z+rxIAIM0YUwwAhXz77beyevVq2bNnjyxfvtzv1QEAZAA9xQAQYefOndKhQwczlljHFI8ZM8YMoahXr57fqwYASCNCMQBEuOmmm2TatGny+eefS7Vq1eSEE06QmjVryhtvvOH3qgEA0ojhEwDwP7NnzzY9w88++6zUqFFDypUrZ/7/gw8+kHHjxvm9egCANKKnGAAAAKFHTzEAAABCj1AMAACA0CMUAwAAIPQIxQAAAAg9QjEAAABCj1AMAACA0CMUAwAAIPQIxQAAAAg9QjEAAABCj1AMAACA0CMUAwAAQMLu/wMrtwA9u+pISwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "## compute electron density with screening\n", + "rho_screened = evaluate_density(rdm, basis, points) ## screen_basis=True, tol_screen=1e-8: defaults\n", + "print(f\"Are the screened and unscreened electron densities identical within 1e-4: {np.allclose(rho_unscreened, rho_screened, atol=1e-4)}\\n\")\n", + "print(f\"Maximum difference between screened and unscreened electron densities: {np.max(np.abs(rho_unscreened - rho_screened)):.0e}\\n\")\n", + "\n", + "## calculate the difference between screened and unscreened electron density for plotting\n", + "delta_rho = (rho_screened - rho_unscreened).reshape(50, 50)\n", + "\n", + "## plot the difference in screened and unscreened electron densities\n", + "fig2, ax2 = plt.subplots(figsize=(7, 6), constrained_layout=True)\n", + "delta_rho_plot = ax2.contourf(X, Y, delta_rho, levels=np.linspace(0, 0.8, 10), cmap=\"viridis\", alpha=0.5)\n", + "ax2.set_title(\"Electron Density Difference (screened - unscreened)\")\n", + "ax2.set_xlabel(\"x\")\n", + "ax2.set_ylabel(\"y\")\n", + "\n", + "## add atomic symbols\n", + "for coord, symbol in zip(xy_coords, symbols):\n", + " ax2.text(\n", + " coord[0], coord[1], symbol,\n", + " color=\"red\", ha=\"center\", va=\"center\",\n", + " fontsize=14, fontweight=\"bold\"\n", + " )\n", + "\n", + "## color bar\n", + "fig2.colorbar(delta_rho_plot, ax=ax2, orientation=\"vertical\", fraction=0.046, pad=0.04)\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "c82ee6c8-1845-4d8d-a6bf-9c5b8a4fcc58", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "We can see that the screened and unscreened electron densities are identical, with the maximum difference being negligible. The plot clearly shows that there is practically no difference in accuracy between the two electron densities." + ] + }, + { + "cell_type": "markdown", + "id": "600d2568-8ba9-435d-881a-25765fcbcffe", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "## Two-index screening\n", + "When calculating one-electron integrals between two basis functions, the overlap is zero or negligible if the functions are far apart. We can estimate a cutoff distance $d_{s,t}$ for each pair of basis functions $s$ and $t$ using an energy tolerance $\\epsilon$, and return zero instead of performing the integral when the distance between the basis functions exceeds this cutoff. The screening is available for the following integrals:\n", + "\n", + "- [`angular_momentum_integral()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.angular_momentum.angular_momentum_integral)\n", + "- [`kinetic_energy_integral()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.kinetic_energy.kinetic_energy_integral)\n", + "- [`moment_integral()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.moment.moment_integral)\n", + "- [`momentum_integral()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.momentum.momentum_integral)\n", + "- [`overlap_integral()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.overlap.overlap_integral)\n", + "- [`overlap_integral_asymmetric()`](http://gbasis.qcdevs.org/_autosummary/gbasis.integrals.html?highlight=evaluate_density#gbasis.integrals.overlap_asymm.overlap_integral_asymmetric)" + ] + }, + { + "cell_type": "markdown", + "id": "0e73f529-e1d9-4845-8593-2e61897f8869", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "### Dipole moment\n", + "To illustrate two-index screening, we calculate the dipole moment integrals with and without screening. We then use the 1-RDM to compute the electronic dipole moment of the system. By comparing both the dipole moment integrals and the resulting electronic dipole moment of the molecule, we can assess the effect of screening on these quantities." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "5262feae-adb0-4260-82ac-f694c74a6ab4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Electric Dipole Moment of the molecule (unscreened): [ 9.39793535e-01 -2.65988737e-08 2.02253056e-07]\n", + "\n", + "Electric Dipole Moment of the molecule (screened): [ 9.39793535e-01 -2.65988737e-08 2.02253056e-07]\n", + "\n", + "Maximum difference between screened and unscreened dipole moment integrals: 0e+00\n", + "\n", + "Are screened and unscreened dipole moment matrices identical within 1e-4? True \n", + "\n" + ] + } + ], + "source": [ + "from gbasis.integrals.moment import moment_integral\n", + "\n", + "## set the orders of the moment integrals\n", + "order = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])\n", + "\n", + "## calculate center of mass\n", + "center_of_mass = mol.atcorenums @ mol.atcoords / np.sum(mol.atcorenums)\n", + "\n", + "## compute dipole moment integrals of a set of molecular orbitals without screening\n", + "electric_moment_integral_unscreened = moment_integral(basis, center_of_mass, order, screen_basis=False)\n", + "\n", + "## calculate the electronic dipole moment using 1RDM and dipole moment integral\n", + "electric_moment_unscreened = np.einsum(\"ij,ijk->k\", rdm, electric_moment_integral_unscreened)\n", + "print(f\"{'Electric Dipole Moment of the molecule (unscreened): ':>40s}{electric_moment_unscreened}\\n\")\n", + "\n", + "## compute moment integrals of a set of molecular orbitals with screening and calculate electronic dipole moment\n", + "electric_moment_integral_screened = moment_integral(basis, center_of_mass, order) # screen_basis=True, tol_screen=1e-8: defaults\n", + "electric_moment_screened = np.einsum(\"ij,ijk->k\", rdm, electric_moment_integral_screened)\n", + "print(f\"{'Electric Dipole Moment of the molecule (screened): ':>40s}{electric_moment_screened}\\n\")\n", + "\n", + "print(f\"Maximum difference between screened and unscreened dipole moment integrals: {np.max(np.abs(electric_moment_integral_unscreened - electric_moment_integral_screened)):.0e}\\n\")\n", + "print(f\"Are screened and unscreened dipole moment matrices identical within 1e-4? {np.allclose(electric_moment_integral_screened, electric_moment_integral_unscreened, atol=1e-4)} \\n\")" + ] + }, + { + "cell_type": "markdown", + "id": "4ccb6c3c-c187-4602-b385-04276381bfef", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "slide" + }, + "tags": [] + }, + "source": [ + "Here, we see that both dipole moments, as well as the dipole moment integrals, are identical, indicating that screening maintains acceptable accuracy." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "gbasis", + "language": "python", + "name": "env" + }, + "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.13.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}