{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8e2d6403",
   "metadata": {},
   "source": [
    "# openCOSMO-RS \n",
    "\n",
    "This tutorial demonstrates how to perform an OpenCOSMO-RS calculation using the ORCA python interface (OPI) and combining it with the [openCOSMO-RS_py](https://github.com/TUHH-TVT/openCOSMO-RS_py) project to calculate sigma profiles and solubilities.\n",
    "\n",
    "In this notebook we will:\n",
    "1.  Import the required Python dependencies\n",
    "2.  Define a working directory\n",
    "3.  Prepare and visualize input structures\n",
    "4.  Run ORCA calculations with COSMO solvation\n",
    "5.  Parse .orcacosmo output files and extract sigma profiles\n",
    "6.  Compute solubility using openCOSMO-RS (non-iterative and iterative methods)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "070ef794",
   "metadata": {},
   "source": [
    "> **Note:** For this notebook to work you will have to install opencosmors_py from GitHub into your OPI `venv` by activating the `.venv` and using:\n",
    "\n",
    "```\n",
    "uv pip install git+https://github.com/TUHH-TVT/openCOSMO-RS_py\n",
    "```\n",
    "\n",
    "For more details visit the [openCOSMO-RS_py](https://github.com/TUHH-TVT/openCOSMO-RS_py) project on GitHub."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9efa50ca",
   "metadata": {},
   "source": [
    "## Step 1: Import Dependencies\n",
    "\n",
    "We start by importing the modules needed for:\n",
    "- Interfacing with ORCA input/output\n",
    "- Plotting results\n",
    "- Numerical calculations and data handling\n",
    "- Handling for directory\n",
    "- Performing COSMO-RS calculations using the openCOSMO-RS library\n",
    "> **Note:** We additionally import modules for visualization/plotting like `py3Dmol`. For this, it might be necessary to install `py3Dmol` into your OPI `venv` (e.g., by activating the `.venv` and using `uv pip install py3Dmol`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "df4e7354",
   "metadata": {},
   "outputs": [],
   "source": [
    "# > Import pathlib for directory handling\n",
    "from pathlib import Path\n",
    "import shutil\n",
    "\n",
    "# > Import necessary libraries for numerical computations\n",
    "import numpy as np\n",
    "from scipy.optimize import fsolve\n",
    "\n",
    "# > Import the openCOSMO-RS library components\n",
    "from opencosmorspy.parameterization import openCOSMORS24a\n",
    "from opencosmorspy.cosmors import COSMORS\n",
    "from opencosmorspy.input_parsers import SigmaProfileParser\n",
    "\n",
    "# > OPI imports for performing ORCA calculations and reading output\n",
    "from opi.core import Calculator\n",
    "from opi.output.core import Output\n",
    "from opi.input.structures.structure import Structure\n",
    "from opi.input.simple_keywords import SolvationModel, Solvent\n",
    "\n",
    "# > Import libraries for visualization\n",
    "from matplotlib import pyplot as plt\n",
    "import py3Dmol"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb47b6f0",
   "metadata": {},
   "source": [
    "## Step 2: Define Working Directory\n",
    "\n",
    "All actual calculations will be performed in a subfolder **opencosmors**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dacfcb12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# > Calculation is performed in `opencosmors`\n",
    "working_dir = Path(\"opencosmors\")\n",
    "# > The `working_dir`is automatically (re-)created\n",
    "shutil.rmtree(working_dir, ignore_errors=True)\n",
    "working_dir.mkdir()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7cd84d3",
   "metadata": {},
   "source": [
    "## Step 3: Prepare and Visualize Input Structures\n",
    "\n",
    "We use **water** as our example molecule. The 3D structure in Cartesian coordinates is defined in XYZ format and visualized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "58f5cf3d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_17514655692888129\"  style=\"position: relative; width: 400px; height: 400px;\">\n        <p id=\"3dmolwarning_17514655692888129\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n        </div>\n<script>\n\nvar loadScriptAsync = function(uri){\n  return new Promise((resolve, reject) => {\n    //this is to ignore the existence of requirejs amd\n    var savedexports, savedmodule;\n    if (typeof exports !== 'undefined') savedexports = exports;\n    else exports = {}\n    if (typeof module !== 'undefined') savedmodule = module;\n    else module = {}\n\n    var tag = document.createElement('script');\n    tag.src = uri;\n    tag.async = true;\n    tag.onload = () => {\n        exports = savedexports;\n        module = savedmodule;\n        resolve();\n    };\n  var firstScriptTag = document.getElementsByTagName('script')[0];\n  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n});\n};\n\nif(typeof $3Dmolpromise === 'undefined') {\n$3Dmolpromise = null;\n  $3Dmolpromise = loadScriptAsync('https://cdn.jsdelivr.net/npm/3dmol@2.5.0/build/3Dmol-min.js');\n}\n\nvar viewer_17514655692888129 = null;\nvar warn = document.getElementById(\"3dmolwarning_17514655692888129\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nviewer_17514655692888129 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17514655692888129\"),{backgroundColor:\"white\"});\nviewer_17514655692888129.zoomTo();\n\tviewer_17514655692888129.addModel(\"3\\n\\n  O  -0.0007948665470900  0.4014278382603300  0.0000000000000000\\n  H  -0.7647999815056600  -0.2022173883142000  0.0000000000000000\\n  H  0.7655948480527500  -0.1992104499461300  0.0000000000000000\\n\\n\",\"xyz\");\n\tviewer_17514655692888129.setStyle({},{\"stick\": {\"radius\": 0.1}, \"sphere\": {\"scale\": 0.3}});\n\tviewer_17514655692888129.zoomTo();\nviewer_17514655692888129.render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_17514655692888129\"  style=\"position: relative; width: 400px; height: 400px;\">\n",
       "        <p id=\"3dmolwarning_17514655692888129\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n",
       "        </div>\n",
       "<script>\n",
       "\n",
       "var loadScriptAsync = function(uri){\n",
       "  return new Promise((resolve, reject) => {\n",
       "    //this is to ignore the existence of requirejs amd\n",
       "    var savedexports, savedmodule;\n",
       "    if (typeof exports !== 'undefined') savedexports = exports;\n",
       "    else exports = {}\n",
       "    if (typeof module !== 'undefined') savedmodule = module;\n",
       "    else module = {}\n",
       "\n",
       "    var tag = document.createElement('script');\n",
       "    tag.src = uri;\n",
       "    tag.async = true;\n",
       "    tag.onload = () => {\n",
       "        exports = savedexports;\n",
       "        module = savedmodule;\n",
       "        resolve();\n",
       "    };\n",
       "  var firstScriptTag = document.getElementsByTagName('script')[0];\n",
       "  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n",
       "});\n",
       "};\n",
       "\n",
       "if(typeof $3Dmolpromise === 'undefined') {\n",
       "$3Dmolpromise = null;\n",
       "  $3Dmolpromise = loadScriptAsync('https://cdn.jsdelivr.net/npm/3dmol@2.5.0/build/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_17514655692888129 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_17514655692888129\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "viewer_17514655692888129 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17514655692888129\"),{backgroundColor:\"white\"});\n",
       "viewer_17514655692888129.zoomTo();\n",
       "\tviewer_17514655692888129.addModel(\"3\\n\\n  O  -0.0007948665470900  0.4014278382603300  0.0000000000000000\\n  H  -0.7647999815056600  -0.2022173883142000  0.0000000000000000\\n  H  0.7655948480527500  -0.1992104499461300  0.0000000000000000\\n\\n\",\"xyz\");\n",
       "\tviewer_17514655692888129.setStyle({},{\"stick\": {\"radius\": 0.1}, \"sphere\": {\"scale\": 0.3}});\n",
       "\tviewer_17514655692888129.zoomTo();\n",
       "viewer_17514655692888129.render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# > define cartesian coordinates in Angstroem as python string \n",
    "xyz_data = \"\"\"\\\n",
    "3\n",
    "\n",
    "  O  -0.0007948665470900  0.4014278382603300  0.0000000000000000\n",
    "  H  -0.7647999815056600  -0.2022173883142000  0.0000000000000000\n",
    "  H  0.7655948480527500  -0.1992104499461300  0.0000000000000000\\n\n",
    "\"\"\"\n",
    "# > Visualize the input structure\n",
    "view = py3Dmol.view(width=400, height=400)\n",
    "view.addModel(xyz_data, 'xyz')\n",
    "view.setStyle({}, {'stick': {'radius': 0.1}, 'sphere': {'scale': 0.3}})\n",
    "view.zoomTo()\n",
    "view.show()\n",
    "\n",
    "# > Write the input structure to a file for reading\n",
    "with open(working_dir / \"struc.xyz\",\"w\") as f:\n",
    "    f.write(xyz_data)\n",
    "\n",
    "# > Load the molecular structure from XYZ file\n",
    "structure = Structure.from_xyz(working_dir / \"struc.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ab68e21",
   "metadata": {},
   "source": [
    "## Step 4: Calculation with openCOSMO-RS Solvation\n",
    "\n",
    "We run a caculation using the **openCOSMO-RS solvation model** with ethanol as the solvent. First, we set up the calculator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "bf0c5f26",
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_calc(basename : str, working_dir: Path, structure: Structure, ncores: int = 4) -> Calculator:\n",
    "    # > Set up a Calculator object\n",
    "    calc = Calculator(basename=basename, working_dir=working_dir)\n",
    "    # > Assign structure to calculator\n",
    "    calc.structure = structure\n",
    "\n",
    "    # > COSMO-RS simple keyword\n",
    "    calc.input.add_arbitrary_string(\"!COSMORS(ethanol)\")\n",
    "\n",
    "    # Define number of CPUs for the calcualtion\n",
    "    calc.input.ncores = ncores # > 4 CPUs for this ORCA run\n",
    "\n",
    "    return calc\n",
    "\n",
    "basename = \"opencosmors\"\n",
    "calc = setup_calc(basename=basename,working_dir=working_dir,structure=structure)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e546858",
   "metadata": {},
   "source": [
    "Next, we run the calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f75d191d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running ORCA calculation ...   Done\n"
     ]
    }
   ],
   "source": [
    "def run_calc(calc: Calculator) -> Output:\n",
    "    # > Write the ORCA input file\n",
    "    calc.write_input()\n",
    "    # > Run the ORCA calculation\n",
    "    print(\"Running ORCA calculation ...\", end=\"\")\n",
    "    calc.run()\n",
    "    print(\"   Done\")\n",
    "\n",
    "    # > Get the output object\n",
    "    output = calc.get_output()\n",
    "    \n",
    "    return output\n",
    "\n",
    "output = run_calc(calc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06069709",
   "metadata": {},
   "source": [
    "Finally, we check that the calculation did terminate normally and converged:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "410ae0b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_termination(output: Output):\n",
    "    # > Check for proper termination of ORCA\n",
    "    status = output.terminated_normally()\n",
    "    if not status:\n",
    "        # > ORCA did not terminate normally\n",
    "        raise RuntimeError(f\"ORCA did not terminate normally, see output file: {output.get_outfile()}\")\n",
    "    \n",
    "check_termination(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e129732b",
   "metadata": {},
   "source": [
    "## Step 5: Analyze COSMO Output and Sigma Profiles\n",
    "\n",
    "After the ORCA OpenCOSMO-RS calculation completes, we parse the resulting .orcacosmo files to:\n",
    "\n",
    "- Visualize the sigma profile (σ-profile)\n",
    "\n",
    "- Extract molecular descriptors such as dipole moment, surface area, volume, and hydrogen-bonding moments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e213f534",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAiCFJREFUeJzt3Xl8VOXZN/DfmTX7CiEBAgShLIKIIAKtLCqgWKutrVat9q3a6qPWKuXxrdpXa32qtvVRqq1SLaJVq7ZFW1upJVVWBRUMimzKDiEheyaZJLOe94+Z+8wkmUxmOWeWM7/v5+OnzSyZM4fJzDXXfV/XJcmyLIOIiIiI0p4h2QdAREREROpgYEdERESkEwzsiIiIiHSCgR0RERGRTjCwIyIiItIJBnZEREREOsHAjoiIiEgnGNgRERER6QQDOyIiIiKdYGBHRLry17/+FZIk4bXXXut33bRp0yBJEv7973/3u+60007DWWedFfHjPPXUU3j++efjOVQiItUxsCMiXVmwYAEkScL69et7Xd7S0oJdu3YhNze333UnTpzAoUOHsHDhwogfh4EdEaUiBnZEpCtDhgzBlClTsGHDhl6Xb9y4ESaTCTfccEO/wE78HE1gpwVZltHd3Z3UYyCi9MbAjoh0Z+HChdi/fz/q6uqUyzZs2ICzzz4bS5cuxY4dO9DR0dHrOqPRiHPPPRcPPPAAzjnnHJSUlKCgoABnnXUWVq1aBVmWlduPGTMGu3fvxsaNGyFJEiRJwpgxY5TrbTYbli9fjqqqKlgsFowYMQJ33HEH7HZ7r+OUJAm33XYbVq5ciUmTJsFqteKFF17Q7sQQke6Zkn0ARERqW7hwIZ544gls2LABV111FQBfVu6rX/0qvvzlL0OSJGzevBlLly5VrjvrrLNQWFiII0eO4KabbsKoUaMAANu2bcMPf/hD1NbW4r777gMAvPHGG/jmN7+JwsJCPPXUUwAAq9UKAOjq6sL8+fNx4sQJ3HPPPTjjjDOwe/du3Hfffdi1axf+85//QJIk5Vj/9re/YfPmzbjvvvtQXl6OsrKyhJ0nItIfBnZEpDvz58+HwWBQArvm5mZ89tln+PWvf428vDycddZZWL9+PZYuXYrjx4/j8OHD+Na3vgUAWL16tfJ7vF4vFixYAFmW8Zvf/Ab/7//9P0iShOnTpyM7OxsFBQWYPXt2r8d+4okn8Omnn+KDDz7AzJkzAQDnn38+RowYgW9+85t4++23cdFFFym37+zsxK5du1BcXJyAM0NEeselWCLSneLiYkybNk3ZZ7dx40YYjUZ8+ctfBuAL/MS+ur776959911ccMEFKCwshNFohNlsxn333Yfm5mY0NDQM+tj//Oc/MWXKFJx55plwu93Kf0uWLIEkSf32/p133nkM6ohINQzsiEiXFi5ciM8//xwnT57E+vXrMWPGDOTl5QHwBXY1NTVob2/H+vXrYTKZ8JWvfAUffvghFi9eDAB49tln8d577+Gjjz7CvffeCwARFTacOnUKn376Kcxmc6//8vPzIcsympqaet2+oqJC5WdORJmMS7FEpEsLFy7EY489hg0bNmDDhg3KfjoA+MpXvgIA2LRpk1JUkZeXh1dffRVmsxn//Oc/kZWVpdz+b3/7W8SPO2TIEGRnZ+O5554b8PpgwfvtiIjixcCOiHRp3rx5MBqN+Otf/4rdu3fjV7/6lXJdYWEhzjzzTLzwwgs4cuQIrr76agC+IMtkMsFoNCq37e7uxosvvtjv91ut1pAZvK9+9at46KGHUFpaiqqqKg2eGRHRwLgUS0S6JFqV/O1vf4PBYFD21wnz58/H66+/DiCwv+7iiy9GZ2cnrr76alRXV+PVV1/Fueeeq1S8Bps6dSo++eQTvPbaa/joo4+wa9cuAMAdd9yBCRMmYN68eXjsscfwn//8B+vWrcMf/vAHXHHFFfjggw80fuZElMmYsSMi3Vq4cCE++ugjTJ8+HQUFBb2umz9/Ph5//HFYLBbMnTsXgK+Q4bnnnsMvf/lLXHLJJRgxYgS+//3vo6ysDDfccEOv+z/wwAOoq6vD97//fXR0dGD06NE4cuQIcnNzsXnzZjzyyCN45plncPjwYWRnZ2PUqFG44IILevW7IyJSmyQHd90kIiIiorTFpVgiIiIinWBgR0RERKQTDOyIiIiIdIKBHREREZFOMLAjIiIi0gkGdkREREQ6kZF97LxeL06ePIn8/HyO8yEiIiLNybKMjo4ODB8+HAaDdnm1jAzsTp48icrKymQfBhEREWWY48ePY+TIkZr9/owM7PLz8wEAhw8fRklJSZKPRr9cLhfWrVuHxYsXw2w2J/twdI3nOjF4nhOH5zoxeJ4Tp6WlBVVVVUoMopWMDOzE8mt+fn6/MUOkHpfLhZycHBQUFPANQ2M814nB85w4PNeJwfOcOC6XCwA03wLG4gkiIiIinWBgR0RERKQTDOyIiIiIdIKBHREREZFOMLAjIiIi0gkGdkREREQ6wcCOiIiISCcY2BERERHpRFIDu6effhpnnHEGCgoKUFBQgDlz5uBf//pX2Pts3LgRM2bMQFZWFsaOHYuVK1cm6GiJiIiIUltSA7uRI0fikUcewfbt27F9+3acd955uPTSS7F79+6Qtz98+DCWLl2Kc889FzU1Nbjnnntw++23Y82aNQk+ciIiIqLUk9SRYpdcckmvn3/xi1/g6aefxrZt23D66af3u/3KlSsxatQorFixAgAwadIkbN++HY8++iguv/zyRBwyERERUcpKmVmxHo8Hf/nLX2C32zFnzpyQt9m6dSsWL17c67IlS5Zg1apVcLlcA865czgccDgcys82mw2Ab26bmN1G6hPnludYezzXicHznDg814nB85w4iTrHSQ/sdu3ahTlz5qCnpwd5eXl44403MHny5JC3ra+vx7Bhw3pdNmzYMLjdbjQ1NaGioiLk/R5++GE88MAD/S5fv349cnJy4n8SFFZ1dXWyDyFj8FwnBs9z4vBcJwbPs/a6uroS8jhJD+wmTJiAnTt3oq2tDWvWrMF3v/tdbNy4ccDgTpKkXj/Lshzy8mB33303li1bpvxss9lQWVmJhQsXorS0VIVnQaG4XC5UV1dj0aJFA2ZTSR0814nB85w4ej3XnQ43rCYDzMbUaEqh1/OcipqbmxPyOEkP7CwWC8aNGwcAmDlzJj766CP85je/we9///t+ty0vL0d9fX2vyxoaGmAymcIGaFarFVartd/lZrOZL+QE4HlOHJ7rxOB5Thw9nWtbjwvzH92EiRUF+PNNobccJYueznOqStT5TY2vDEFkWe61Hy7YnDlz+qWL161bh5kzZ/IFSUREKe1wox22Hjd2nWhP9qGQjiU1sLvnnnuwefNmHDlyBLt27cK9996LDRs24JprrgHgW0K97rrrlNvffPPNOHr0KJYtW4a9e/fiueeew6pVq7B8+fJkPQUiIqKItHf7Ns93uzxwe7xJPhrSq6QuxZ46dQrXXnst6urqUFhYiDPOOANvv/02Fi1aBACoq6vDsWPHlNtXVVVh7dq1uPPOO/G73/0Ow4cPxxNPPMFWJ0RElPLaugNVkXaHB4U5KbdoRjqQ1MBu1apVYa9//vnn+102f/58fPzxxxodERERkTbagwK7DocLhTncQkTq49cFIiKiBLAFBXadDncSj4T0jIEdERFRArR1OZX/39nDwI60wcCOiIgoAXovxTKwI20wsCMiIkqA9l7FEwzsSBsM7IiIiBKgrStojx2XYkkjDOyIiIgSoJ3FE5QADOyIiIgSILgqtoMZO9IIAzsiIqIEaGPGjhKAgR0REZHGXB4vupwe5WfusSOtMLAjIiLSWPD+OoAZO9IOAzsiIiKNBVfEAuxjR9phYEdERKSxfhm7HtcAtySKDwM7IiIijdm4FEsJwsCOiIhIY23dvjmx+VkmACyeIO0wsCMiItJYu3+P3cjiHADcY0faYWBHRESksfZuXyA3oigbgG8pVpblZB4S6RQDOyIiIo2JpdiRxb7ATpbRq68dkVoY2BEREWlMVMUOK8iC0SABAOxcjiUNMLAjIiLSmNhjV5RjRp7VV0DBfXakBQZ2REREGhMZu6LsQGDHyljSAgM7IiIijYnArjA4sGPGjjTAwI6IiEhjbf7AriDbjDx/L7sOZuxIAwzsiIiINKYsxeYwY0faYmBHRESkoR6XB063F4B/KVaZPsF5saQ+BnZEREQaavNXxBoNEvKsJuQzY0caYmBHRESkoeDCCUmS2O6ENMXAjoiISEPBgR2AoKVYBnakPgZ2REREGmrr8o0TKxCBHZdiSUMM7IiIiDQU3JwYAPKZsSMNMbAjIiLSUL+lWKvvf7nHjrTAwI6IiEhD3GNHicTAjoiISEPBzYkBIM9qBADYnQzsSH0M7IiIiDQ00FIsM3akBQZ2REREGhINigv6LMVyjx1pgYEdERGRhvpWxYp2J063Fw63J2nHRfrEwI6IiEhDtn5LsSblOruDgR2pi4EdERGRhtpEYOcvnjAaJORYfAUU3GdHamNgR0REpBFZloOWYi3K5YF5sa6kHBfpFwM7IiIijdidHni8MoDAUizAXnakHQZ2REREGhFzYi1GA7LMgY/cfM6LJY0wsCMiItJIe9D+OkmSlMuVjB0DO1IZAzsiIiKN9G1OLCh77LgUSypjYEdERKSR9q6BAjv/9Alm7EhlDOyIiIg00rc5sZDP4gnSCAM7yhhHmuy4/++f4UBDZ7IPhYgyRNsAS7G5Vn8fO2bsSGUM7ChjvLTtKF7YehRX/H4rPqttT/bhEFEGEBm7Ai7FUoIwsKOMITYpt9iduPrZbag51prkIyIivVOWYnP6BHZciiWNMLCjjNHjH7ZtNRlg63HjO3/4AB8cak7yURGRng1UPME+dqSVpAZ2Dz/8MM4++2zk5+ejrKwMl112Gfbv3x/2Phs2bIAkSf3+27dvX4KOmtJVj8sX2C1fPAFzTyuF3enBd1d/iC1fNCX5yIhIrwZtd8LAjlSW1MBu48aNuPXWW7Ft2zZUV1fD7XZj8eLFsNvtg953//79qKurU/4bP358Ao6Y0lmPywsAKM614Ln/czYWTBiKHpcX17/wEd7ZeyrJR0dEejT4UixnxZK6TMl88LfffrvXz6tXr0ZZWRl27NiBefPmhb1vWVkZioqKNDw60htH0FJsltmI3187Az/8Uw3W7TmFm1/agSe+PR0XTa1I8lESkZ60dftGig2UseNSLKktqYFdX+3tvkrFkpKSQW87ffp09PT0YPLkyfjpT3+KhQsXDnhbh8MBh8Oh/Gyz2QAALpcLLhe/LWlFnNtUOcfdTl9gZ5ZkuFwuGACsuGIq/nuNhLd21eO2V2rwS4cLl05Lv+Au1c61XvE8J45ezrXYY5djkno9F3/CDp097qQ+R72c53SQqHMsybIsJ+SRBiHLMi699FK0trZi8+bNA95u//792LRpE2bMmAGHw4EXX3wRK1euxIYNGwbM8v3sZz/DAw880O/yP/3pT8jJyVHtOVBq++UnRpzskvBfkzyYWBR42Xtl4JWDBnzYaIAEGVeO9WLOsJT4syCiNOaVgWXbjJAh4cEZbhRYAtd1uoB7t/uiu8dnu2GQBvglpBtdXV24+uqr0d7ejoKCAs0eJ2UCu1tvvRVvvfUWtmzZgpEjR0Z130suuQSSJOHNN98MeX2ojF1lZSXq6upQWloa13HTwFwuF6qrq7Fo0SKYzebB76CxRSu24EhzF1658WzMHF3c6zqvV8YDb+3Fnz48AQC4Zf5Y3LZwLMzG9CgcT7VzrVc8z4mjh3Pd1uXC2Q+vBwDsvv8CWEyB9xOH24spD/wHAPDxvQuRn5Wc56iH85wumpubUVFRoXlglxJLsT/84Q/x5ptvYtOmTVEHdQAwe/ZsvPTSSwNeb7VaYbVa+11uNpv5Qk6AVDnPDreveCIvyxryeH7x9TOQYzHjD1sO46mNh7D5QDMev3IaxpXlJ/pQY5Yq51rveJ4TJ53PdZfbt78u22xEbnbvzyCzGbAYDXB6vOjxSChJ8nNM5/OcLhJ1fpOajpBlGbfddhtef/11vPvuu6iqqorp99TU1KCiIv32RVFiiXYnWebQL3tJkvDTr07Gk1dNR2G2Gbtq23HxE1uw+r3D8HpTIrFNRGlkoIpYQamMZQEFqSipGbtbb70Vf/rTn/D3v/8d+fn5qK+vBwAUFhYiOzsbAHD33XejtrYWf/zjHwEAK1aswJgxY3D66afD6XTipZdewpo1a7BmzZqkPQ9KDyJjZzUZw97ukmnDcfaYEvz3Xz/B5i+a8MA/9uCdvQ349bfOQEVhdiIOlYh0oG2A5sRCrtWIFntgKg6RGpKasXv66afR3t6OBQsWoKKiQvnvtddeU25TV1eHY8eOKT87nU4sX74cZ5xxBs4991xs2bIFb731Fr7xjW8k4ylQmpBledCMXbDywiz88fpZePDS05FlNmDLgSYseXwT/r6zVutDJSKdGKg5scB5saSFpGbsIqnbeP7553v9fNddd+Guu+7S6IhIr1weGWI11WoOn7ETJEnCtXPGYO64IVj250/wyfE2/OjVnajecwr/c9kUFOVYBv8lRJSxBgvsxFgxOwM7UlF6lPwRxUnMiQUiy9gFO21oHtbcPAd3XvAlGA0S/vlpHZas2ITPT3WofZhEpCODZuyU6RMM7Eg9DOwoI4hlWEnyVaJFy2Q04EcXjMcbt8zF2KG5OGVz4JlNh9Q+TCLSkUGLJzgvljTAwI4ygsMlCicMkKTYO4GeMbIId17wJQDA0ebBZxoTUeZqH6R4ghk70gIDO8oIYk5sVoT768IZVeKbVnKspSvu30VE+jXQnFghX5kXy3FepB4GdpQRevwZu6xBWp1EQgR2p2wOZYmXiKgvZY/dAIVWeVb2sSP1MbCjjBBNq5PBFOWYlTfkE63dcf8+ItKn9m5fwDbYUiz72JGaGNhRRlAydiosxUqShEp/1u44l2OJaADtXeGXYpmxIy0wsKOMIDJ2kfawG8yoEt8ECq332Z1o7cJ5/7sBL249ounjEJH6lKrYgfbYsXiCNMDAjjJCYJyYOi/5RBVQbNjfiEONdrz5yUlNH4eI1OXyeGF3+r5QcvIEJRIDO8oIgT126mTsErUUe7LNt4ePe3CI0ovI1gFAAffYUQIxsKOMICZPZKmUsatMUMaOgR1RehKBXX6WCUZD6N6ZeVbfF01m7EhNDOwoI6hZPAEElmKPt3RFNPM4VifbegAAHT3sc0WUTtoGaU4M9F6K1fJ9hDILAzvKCGq2OwGAEUXZkCTA7vSgtUu7oKvWn7HjGz9RerENMicWCCzFeryysg+YKF4M7CgjOERVrAoNigFf5m9YfhYA7ZZjPV4Z9TZfxs4rA11ONkMmSpYWuxNuT+TB12BzYgEgx2yEmHDI7RakFgZ2lBHEt2G1MnaA9pWxDR098HgDWTq+8RMlx4nWLpzz0H9w80sfR3yftkF62AGAwSAhz8JedqQuBnaUEdSuigW0r4yt7TPVgvvsiJJjX10HXB4Z7x9sinhLxGBTJ4Q89rIjlTGwo4ygdvEEAFT6mxRrFti19Qns+I2eKCla/Nm3LqcHJ9t7IrqPMic2O/ScWEFMn+hw8IsbqYOBHWUE0e5ErQbFgPZLsaIiVuBSLFFytNqdyv//4lRHRPdp6x58KRZgxo7Ux8COMoIWS7HaB3ZciiVKBSJjBwAHGjojuk8kVbEA58WS+hjYUUYQS7FaZOzq2nvgiqJaLlJ9Azt+oydKjuCMXaSBXSRVsUDQvFgGdqQSBnaUERxu9TN2Q/OtsJoM8Hhl1LVFtu8mGmKPXbH/g4FLsUTJ0WIPZMu/iDCwi6RBMRC0x45/36QSBnaUEbQonpAkSdPRYiJjN6E8HwCXYomSpbWr9x67SCpj2yNeig1MnyBSAwM7yghqT54QtNpn19Hjgs3/DX5ieYHvMr7xEyVF8FKsrceNxk7HoPeJPLDzz4tlxo5UwsCOMkKgQbF6GTsAqCz2tzxpVTewExWxRTlmlBVYAXCphihZRPGEyeAbE3HgVPjl2B6XR3nPKRxkj10e99iRyhjYUUZQMnYqjRQTtFqKFcuwwwuzkZ8l9thxKZYo0dwer5J9mzqyEMDg++zE7Q0SlMkSAxFLsfziRmphYEcZQXx7tmq0FKt2k2JRODG8KBv5bIdAlDTt3S6ILXVnjykBMHhlbPAyrMGf5RtIIGPHL26kDgZ2lBG0ytiNKtUmsBMZuxFFWUo7BH6jJ0q81qCZrxOG+QqZvmgI36Q40opYAMoXN7vDE89hEikY2JHuybKsWfFEZbEvsGvtcsGm4lLpyeCMXRaXaoiSRbQ6Kcm1YPywPADRZewGwz12pDYGdqR7Lo8Mr38pxapy8USu1YTSXN8sSDWzdqJ4YnhRNvtcESVRi78itjjHjNOG+gK7pk5nr0rZvpTALif8nFiAfexIfQzsSPfEnFhA/YwdECigUDOwE3vsRhRnBy3Fcg8OUaKJpdiSXAtyrSaMKPJVwh9oHDhr19YV2ZxYIHikGP++SR0M7Ej3xDKsJAEWo5aBXfcgt4yM2+NFvc2XsRtRlI0C/1Ksw+2F063+6DIiGlggY+fLvo0r82XtvgjT8iQwJzZ8RSwQGCnW4/JqMpqQMg8DO9I9R9CcWEkKX6EWi1Elvm/warU8aehwwOOVYTZKGJpnRa41sHzMfThEiSWWXEv8Wy7Glw2+z06ZE5s9+FJsrjUQ/Nn5900qYGBHuqfFnNhgak+fEIUT5YVZMBgkmIwG5Fh8x87lWKLEEs2Ji0Vg5y+gCFcZ2xZF8YTZaFC2iHCfHamBgR3pnjInVuVWJ4KyFKvS9InaoObEAlueECWHkrHrsxQbScZusKkTAufFkpoY2JHuadXqRBAtT060dMPrHXw4+GBERazYpA2wco4oWVr8PelExm7cUF8vu7r2ngEz6NG0OwECX9wY2JEaGNiR7ikZO42WYisKs2AySHB6vDjV0RP37wvuYSdwrBhRcgT22Pn+BgtzzCjL981vPthoD3mf9igaFANQ9tF28osbqYCBHeme2GNnNWnzcjcZDRhR7C+gaI5/OfZkUKsTgUuxRMnR2qcqFgjaZ3cq9D47pXgi4qVY/983M3akAgZ2pHsiY6d2c+JgYjn2eGv8LU9qQ2TsCpixI0o4p9urBFuiKhYAxg0deJ+dLMtRFU8AQXvs+MWNVMDAjnQvsMdOw8BOxcrY2qA5sUKgiSnf+IkSRTQaNkiBL1cAMM4/MzZUYGd3euDx77WNfo8dv7hR/BjYke6JyRNZGi3FAoGWJ/FOn7D1uJTl1gpWxRIlldLqJMcCgyHQA1P0svsiRGAnlmEtRgOyI/wyqXxx4983qYCBHeme1sUTgHqBXZ2/IrYox9yrcakonrDxjZ8oYZSpE7m9Gw2LlifHW7uUFQFBZPkKss0RN0TPUzJ2nkFuSTQ4Bnake+KNV6viCQCoVGn6xMkQPeyA4Dd+BnZEidJq92XfSnJ6B3aluRYU55ghy8DBPjNj26MYJyZwXiypiYEd6Z7DnbiMXUOHA93O2L91hyqcAIKXYvnGT5QogakTvffKSZKE8WWh99nZlIrYwceJCexjR2piYEe659C4QTHg2yQt3pxPxDGB4mSIwgkAKOAeO6KE6zsnNtg4peVJ78CuLcoedgAbkJO6GNiR7iWiKlaSpMA+uzgCu9oQPewAtkMgSoaWED3shIFankQ7dQJg1Tupi4Ed6V4iiieAQC+7eJoUh5o6AXApligZWrsGztgpTYobejcpjimwy2JVLKmHgR3pXo/GkyeEUaWil13sTYrFnNiBAzu+8RMlSriMndhjd6S5C07/Pl4AUTcnBoB8kZFnxo5UwMCOdE+pitU6YxfnUqzb40W9zRfYjSgaoCrW6YbX3/yUiLQVLmM3rMCKPKsJHq+MI82BmbHM2FGyJTWwe/jhh3H22WcjPz8fZWVluOyyy7B///5B77dx40bMmDEDWVlZGDt2LFauXJmAo6V0pVTFapyxq/Tvi4u1l11DhwMerwyzUcLQPGuv60TXe1kG7E6++RMlgmh30rePHeDbVyv62QXvs7NFOScWAHKtvi+d/OJGakhqYLdx40bceuut2LZtG6qrq+F2u7F48WLY7fYB73P48GEsXboU5557LmpqanDPPffg9ttvx5o1axJ45JROElE8AQRanhxr6YIsR//mLPbXlRdm9epyD/iWkc1G32VcjiVKDLEU27ePnaBMoAiqjI2lKlYsxcoy0OVik2KKT+QdFDXw9ttv9/p59erVKCsrw44dOzBv3ryQ91m5ciVGjRqFFStWAAAmTZqE7du349FHH8Xll1+u9SFTGkpU8cSI4mxIEtDl9KDZ7sSQPlm3wdQO0JwY8GUH8qwmtHa5uA+HKAG6nR50+4Osvn3shHFl/QsoYlmKzTIbYDRI8HhldPa4lSpZolik1Kunvb0dAFBSUjLgbbZu3YrFixf3umzJkiVYtWoVXC4XzOb+f0wOhwMOh0P52WazAQBcLhdcLlYZakWc22Sf4x7/0qVJkjU9FgOA8oIs1LX34HCDDYXWoqjuf9y/T2d4oTXkcYrArrWzB66S3n3uUuVc6x3Pc+Ik+1w3+ve7mo0SrIbQ7x1Vpb4vYQdOdSjXt3f7sny5ZimqY8+zGtHe7UZrZzdKc7T9Ehos2ec5kyTqHKdMYCfLMpYtW4avfOUrmDJlyoC3q6+vx7Bhw3pdNmzYMLjdbjQ1NaGioqLffR5++GE88MAD/S5fv349cnJy4j94Cqu6ujqpj9/cbgQgoWb7B7B9ru1j5Xh9j/WP9Vtxckh0y7FbDxkAGGBvrMXatcf7XS87fb/73c1bUV8c+ncn+1xnCp7nxEnWuT5hBwATsg1e/Otf/wp5m+Ye320ONHTgH2+thQTA1u37O93+/iZ8HvnwCRj97x3r1m/C5/nxHn30+JrWXldXfCMnI5Uygd1tt92GTz/9FFu2bBn0tn0HK4v9TAMNXL777ruxbNky5WebzYbKykosXLgQpaWlcRw1heNyuVBdXY1FixaFzKQmyiN7NgE9PVhw7pcxdUShpo+1seczHKw5idJRE7B0wdio7vvGix8Dp5owb+YULJ05st/1L9d9hBP2Vkw6YzqWTi3vdV2qnGu943lOnGSf6y0HmoFPd6CiJB9Ll84NeRuvV8avPnsHPS4vps6ej+IcC+Rt6wEAX//qhVG1WHrq0PtoOdWJaTPPwZdPS9znUrLPcyZpbm5OyOOkRGD3wx/+EG+++SY2bdqEkSP7f6AFKy8vR319fa/LGhoaYDKZBgzSrFYrrNb++53MZjNfyAmQ7PMsqmLzsq2aH8eYIb49NyfbHVE/Vr3Nt12gsjQv5H0Lsn1f/7tc8oC/O9nnOlPwPCdOss61zeHbX1eSG/5947Shedh90oYjLQ5Yzb6/0WyzEXnZ0e2xzfdXvve4B/771hJf09pL1PlNalWsLMu47bbb8Prrr+Pdd99FVVXVoPeZM2dOv5TxunXrMHPmTL4oKSSlKtak/b6VyhLfnptjMbQ8qR1gTqxQwOkTRAkTbk5ssPFBBRRt/v110RROCHlsQk4qSWpgd+utt+Kll17Cn/70J+Tn56O+vh719fXo7g507r/77rtx3XXXKT/ffPPNOHr0KJYtW4a9e/fiueeew6pVq7B8+fJkPAVKcbIsB7U70f7lHtzyJBq2Hpfyhl4RoioWCGpiyqpYIs21dIkeduGDtPHDfBviDpzqjKkiVuC8WFJLUgO7p59+Gu3t7ViwYAEqKiqU/1577TXlNnV1dTh27Jjyc1VVFdauXYsNGzbgzDPPxIMPPognnniCrU4oJJdHhuj3qfXkCSAwfaKuvRsuj3eQWwfU+UeJFeWYkTtAqwOOFSNKnNZBetgJpw31NyluDArsomhOLORz+gSpJKl77CJp4vr888/3u2z+/Pn4+OOPNTgi0hsxJxbQflYsAAzNsyLLbECPy4uTbd0YXZob0f1q23wZvr6jxIKJPTg2LsUSaa7FP04s1NSJYOOHBaZPtMbQnFhgxo7UwlmxpGsOf3NiSUpMYCdJEiqLo1+OrfVn7IaHDez4jZ4oUSLdYze6JAdmo4Qupwd763w9UmML7Hz36WBgR3FiYEe6JvbXWU2GAdvhqC2WfXYnlcKJgQM78Y2eS7FE2hPjxIoHWYo1GQ0Y66+G33GkFQBQFENgp8yL5d83xYmBHemaw52YObHBKuMI7IYPUBELAAVZ4hs9l2KJtNbaFVnGDgiMFvvcP1osloxdPoujSCUM7EjXlDmxCWh1IojA7kRL9yC3DAgEdoMvxTJjR6QtWZbRahdVsZEHdmLbeCzFE2Iplhk7ihcDO9K1RLY6EWJbih18j10e99gRJYTd6YHTX9U+WFUsECigEOLqY8eMHcWJgR3pmsjYWROYsYs2sHN7vKj3DxyPpCqWGTsibYnCiSyzAdmWwd87RMZOiK8qllstKD4M7EjXAnvsEvdSF9Mn2rtdSl+rcE51OODxyjAbJQzNG3gMkViKdXq8SiaSiNTXEmEPO6FqSC4MQbVZce2x4xc3ihMDO9I1JWOXwOKJHIsJQ/J8HwjHI8jaif11FYXZMBgGrtzNtQTaTnKDNZF2Iu1hJ1hNRowJ6llZFGFAGIx97EgtDOxI1wJ77BIX2AGBAopoArtwFbEAYDRIbHlClACR9rALFrwcG88eO5dHVlYaiGLBwI50TUyeyEpAc+Jg0eyzq42gIlYIVMZyHw6RViLtYRcsOLAryIp+qFOvjDy/uFEcGNiRrintThKdsfNPnzjQ0DnobSNpTiwoyzV84yfSTDQ97ARRGZtvNcFkjP6j1WiQkOsv1OByLMWDgR3pmljSSMQ4sWBnV5UAAP6+8yQONYYP7iJpdSKIjJ2NgR2RZlpED7soMnZTRxQBAEb6s/WxyGOvSlIBAzvStWRl7OaNH4L5XxoKp8eL+9/cDVl0Lg0hkubEQqDlCZdiibQS2GMX+V65cWV5ePnGc/D0NWfF/LgsoCA1MLAjXXMkoUExAEiShAe+djosJgM2f9GEtbvqB7xtrbIUG754AghqUsw3fiLNRFsVK3x53BCMGZI7+A0HkJfF6RMUPwZ2pGvJqooFgDFDcvFf808DADz4zz0hgzFbj0tZdokkY1fApRoizbVG2cdOLXlW7rGj+DGwI11L1lKs8F8LTsOokhzU23rwxDtf9LteLMMW55iRYxm8ko5LsUTaa40xYxcvpZ0RAzuKAwM70rWeJBVPCFlmIx742ukAgFVbDmN/fUev66PZXwdwDw6R1rxeGa1dvi9O0VTFqiHPyqVYih8DO9I1RxImT/S1cGIZFk8eBo9Xxv/7+2e9Cilqo6iIBVgVS6S1jh43PF7f32hRTvSNhuOhjBXjvFiKAwM70rVkNSju675LJiPbbMSHh1vwRk2tcnk0PeyA4KVYBnZEWhCFE3lWE6ymxH4hZJ9KUgMDO9K1ZBZPBBtZnIMfnj8OAPDQ2r1o7/Z9I490nJgQGBTOb/REWlCmTkTR6kQtSh87brWgODCwI11LdvFEsBu/MhanDc1FU6cT/7tuP4Do99jlc1YskaaSVRELMGNH6mBgR7rWk6Q+dqFYTAY8eOkUAMBL247is9p21LZyKZYolcTaw04N+exTSSpI/qcdkYYcbn/xRIL3ygxk7rgh+Nq04fDKwL1v7EK9zVc8EXlgxzd+Ii2lQsbOzr9vikNMgZ3L5cLx48exf/9+tLS0qH1MRKpJ1uSJcH568STkWU345EQ7vDJgNkoYkmeN6L7BkydE5R4RqSeZGTv2sSM1RPxp19nZid///vdYsGABCgsLMWbMGEyePBlDhw7F6NGj8f3vfx8fffSRlsdKFLUed+rssRPKCrKwbNGXlJ8rCrNhMEgR3Vdk7ABm7Yi0EJgTm4TALot77Ch+EQV2jz/+OMaMGYNnn30W5513Hl5//XXs3LkT+/fvx9atW3H//ffD7XZj0aJFuPDCC/HFF/077BMlg7LHLkWWYoXr5ozGpIoCAJFXxAK+JWWLv3ULAzuKxL931+Oi32zG56c6Br8xocXuqzgvTsJSbL5oUMy/bYrD4DOMALz//vtYv349pk6dGvL6WbNm4frrr8fKlSuxatUqbNy4EePHj1f1QImiJctyShVPBDMZDfjV5Wfgjtdq8I2zRkZ133yrCc1up3+sWGR78yhzvbnzJPbW2fD2Z/X40rD8ZB9OyhPjxEqS0O4k1z8rtsvpgccrwxhhJp8oWESB3V/+8peIfpnVasUtt9wS1wERqcXlkSG2oaVK8USwqSML8c6PF0R9v/wsE5rtTlbGUkTsTt/rRBTqUHhiKTYZGbu8PlstCrMTH1xS+osqjWG325X/f+TIEbWPhUhVDv/UCQCwpljGLh6BlidsUkyDExWWp9oZ2EWipSt5e+ysJiMsRm61oPhE/Gl3++23Y+TIkXj66acBAFdffbVmB0WkBtGcWJIAa5JHiqkpj02KKQp2h+8LTh0Du0F5vLIyFSYZVbEACygofhF/2q1btw6nTp1CTU0N1qxZo+UxEalC7K+zmgyQJP3sVRGVsQzsKBJd/qXYU1yKHVR7twuyf/tGUZKWQZXpEw5m5Ck2Ee2xA4CRI0fCYrHg6aefxte//nWcPHlSy+MiiptYik2lVidq4PQJikanP2PXbHeix+XR3d+DmsSc2MJsM0zG5GT5mZGneEX8yp00aRJcLheMRiN+//vfo6ioSMPDIoqfMic2BQsn4hGYPsFv9DQ4kbEDgAabI4lHkvpak7i/TsjjdBmKU8SB3ZNPPgmz2ZcpqKiowM6dO/vdhg2KKZUoS7E6KpwAuBRLkfN6ZXQ5A0VErIwNr0WpiE1eNWq+lXvsKD5xf+I1NDTgsccew5QpUzB79mw1jolIFWJOrF4zdgzsaDBdLk+vn+vau5N0JOkhmVMnBGbsKF4xBXYejwd///vfcdlll6GyshLPPvssLrvsMmzfvl3t4yOKWao2J45XnpV77CgyXX2Cg3pWxoalzIlNQg87IVA8wb9vik3Y4onGxkY89thjKCkpwR133IH9+/dj9erVeOmllwAAV1xxBbxeL9asWYPJkycn5ICJIiX22Fl1tlk8kLHjHjsKr29wwKXY8FIqY8cvbhSjsKmMq6++Gl1dXQCAESNGYPbs2Th58iSee+45nDx5Ek8++WRCDpIoFoGMnV4DO77xU3jB++sAZuwGo8yJTWJgl8+MHcUpbGC3b98+XHPNNbj++uvR0tKCH/zgB/j5z3+Oiy++GEajvj4sSX96RLsTHTUnBoKrYvnGT+ExYxed1q7kF08o7U74900xCvuJ99Of/hRf//rXMX/+fDzyyCM4cuQIpkyZgnPOOQe//e1v0djYmKjjJIqafpdiOVKMIiNanYhh8szYhdeSxDmxQi6rYilOYffY3XTTTbjmmmtgtVqVVieNjY146aWX8Oyzz+LOO++E1+tFdXU1KisrkZ+fn5CDJoqEQ+cZu44eN2RZ1tVUDVKXGCc2ujQHhxrtaOhwwOOVlUCPekuFPnbMyFO8Bv3Ey8vLU4I6ABg6dCjuvPNOfPLJJ9i2bRv+67/+Cw8++CDKysrwta99TdODJYqG0qBYpxk7t1dWWrpQ4njlZB9B5Oz+4GB0SQ6MBgker4ymTjYpHoiSsUtm8YS/6p0ZO4pVXKmMGTNm4Le//S1OnjyJF154AW43X4iUOhw6bXeSYzZCJOlsXI5NqB//ZRf+3w6jktlJdXZ/8UR+lhll+VYAQB2XY0NyebxKQVJJMtudMGNHcVLlE89iseCKK67A2rVr1fh1RKrQa1WswSBxnmQS9Lg8WPtZPTpdEvbVdyT7cCIiMna5VhPKC7MAAPVsUhySCNYNElCQnQLFE/zSRjGKKLC7+eabcfz48Yh+4WuvvYaXX345roMiUoNSPKGzPXYAUJDF5ZpE+/REO9z+ddi2rvT40LX7iydyLUaUF4jAjhm7UFr9rU6KcixJ3YMYvMdOltNo3Z9SRtjiCWHo0KGYMmUK5s6di6997WuYOXMmhg8fjqysLLS2tmLPnj3YsmULXnnlFYwcORLPPPOM1sdNNCileEJnGTsAzNglwY6jrcr/t6XJee/yF08EZ+zq2PIkpFSYEwsEvrR5Zd9SuvhbJ4pURK+YBx98ED/84Q+xatUqrFy5Ep999lmv6/Pz83HBBRdg1apVWLx4sSYHShQtvbY7ATh9Ihk+PhYI7Nq70+O8B5Zijcix+P4OTjFjF1IqVMQCvj3BZqMEl0eGrdvFwI6iFvErpqysDHfffTfuvvtutLW14ejRo+ju7saQIUNw2mmnseUCpRy9NigGggI7brBOCFmW8fHRNAzs/EuxORaT8pph8URoqdDDDgAkSUJBlhnNdidsPS4MR3ZSj4fST0yfeEVFRZg2bRpmz56NcePGxRXUbdq0CZdccgmGDx8OSZLwt7/9LeztN2zYAEmS+v23b9++mI+B9EmvxRMAkKc0KWZglwhHm7vQbA9UwrZ3p8d5F33s8qwmVBT6AgROnwgtFebEChwbSPFIeo7Xbrdj2rRp+N73vofLL7884vvt378fBQUFys9Dhw7V4vAojem1jx3ApdhEC16GBdIxY9e7eIKNrftr6Up+DztBVOXa0uR1Rqkl6YHdRRddhIsuuijq+5WVlaGoqEj9AyLdEBk7PVbF8ht9YonCiaF5FjR2OtPmAze4eKKswNfHzuH2oq3LlRIBTCpRMnZJXooFAgUU7FNJsUh6YBer6dOno6enB5MnT8ZPf/pTLFy4cMDbOhwOOByBbus2mw0A4HK54HLxD0cr4twm6xyLBsUmSdbdv3OOP1i1dTt7vY719jxTxY4jLQCAc8eV4vWddWjzn/dU1+nwHaPVCBjhRUmuGS12F443dyLPktojIBP9mm72T+QoyDIk/d821+L7+27tdGh+LHzvSJxEneO0C+wqKirwzDPPYMaMGXA4HHjxxRdx/vnnY8OGDZg3b17I+zz88MN44IEH+l2+fv165OTkaH3IGa+6ujopj9vaYQQgYceHW9G4JymHoJljdRIAI744fBxr1x5VLk/WudazHjew/5TvtVTUdQKAEfUtHWnRkL2t0/83sO091OYA2bLv53++uwWHi9OjR1qiXtNH6nzn5uCeT7G27pOEPOZA2psMAAzY8elulLZ8Nujt1cD3Du11dXUl5HHSLrCbMGECJkyYoPw8Z84cHD9+HI8++uiAgd3dd9+NZcuWKT/bbDZUVlZi4cKFKC0t1fyYM5XL5UJ1dTUWLVrUa95wovzsk/WA04Xz58/D+GF5CX98LTl3nsSaI58ht3goli6dkfRzrWfvHWyG/NEOjCjKwuUXTMNzn38AF8xYunRJsg9tUP/90X8AeHHRovNQUZiFv7V8jNr9Taj80lQsPXtk0o6rocOBoXmWsPv8Ev2a/tXeTQB6sGjeHJxZWaT544Xz6dv7sa3hKMpHjcXSCycMfoc48L0jcZqbmxPyODEFdn/961/x5z//GceOHYPT2Xtm4scff6zKgUVj9uzZeOmllwa83mq1wmq19rvcbDbzhZwAyTrPDreveCIv26q7f+fCHN/r2e709HpufE2r75MTvvFhM0aXoDTfV4DQ4XDDYDQldULBYFweL5z+v4HCnCyYzWZUFPlWKBrtrqS9Tt7+rA43v/Qx/nvJBNy6cNygt0/Ua7rVP02krDAn6X9DRcrftzdhx8L3Du0l6vxGvav8iSeewPe+9z2UlZWhpqYGs2bNQmlpKQ4dOhRTEYQaampqUFFRkZTHptQky3KgeMKsx+IJ0e6E+2K0JipiZ4wu7jVDNNULKEThBADkWH2V4RUFyZ8Xu6u2HQCw7VBisheR6HF50OX0na9UKCphcRTFI+qM3VNPPYVnnnkGV111FV544QXcddddGDt2LO677z60tLREfQCdnZ04cOCA8vPhw4exc+dOlJSUYNSoUbj77rtRW1uLP/7xjwCAFStWYMyYMTj99NPhdDrx0ksvYc2aNVizZk3Uj0365fbK8I/1RJZJz+1O+MavJa9XVgK7s0YVw2w0wGqQ4fBKaO9O7crSTn+rE4vJALPR9+VGGSuWxCbFIjN2uMmetGPoS0ydMBkk5KfApAel3Qm/uFEMon4FHzt2DHPnzgUAZGdno6PDt0xx7bXXYvbs2fjtb38b1e/bvn17r4pWsRfuu9/9Lp5//nnU1dXh2LFjyvVOpxPLly9HbW0tsrOzcfrpp+Ott97C0qVLo30qpGMiWwfoNWMXGBRO2jnQ2ImOHjeyzUZMrMgHvB7kmACHE2hL+Yydf5yYJfDFRgR2p5LYpLjdH9jVtnWjx+VJiT6TytSJ3PD7/hJFaXeS4q8xSk1RB3bl5eVobm7G6NGjMXr0aGzbtg3Tpk3D4cOHIcvRV1ktWLAg7P2ef/75Xj/fdddduOuuu6J+HMosojmxJOm1j53vjb/L6YHb403y0eiXGCM2rbIQZqMBLn9g1+pM/SbFdmegh51QkRIZO18QJcvAsZYufGlY8tuutNp9/5ap0MMOCM7Y8YsbRS/qT7zzzjsP//jHPwAAN9xwA+68804sWrQIV155Jb7+9a+rfoBEsQhuTpwK38DVFjwYnFk77YjGxGeNKlYuyzb5voi2dTlD3idV2JWMXeC1Msy/x66jx61cn2htXYGA+FBjaizHBqZOpEbxACfLUDyiztg988wz8Hp9GYKbb74ZJSUl2LJlCy655BLcfPPNqh8gUSwcbv3OiQV8+6asJgMcbi86etzIZTWbJnYEFU4IOf53zVRfJhOBmyicAHyZ3jyrCZ0ON+ptPThtaOLbAAUHxKmyzy6V5sQCwSPF3Bz/RlGLOrAzGAwwGAKJviuuuAJXXHGFqgdFFC+xFKvHZVghP8sMR6cDHT1ulOczsFNbq92pZJSmj+of2KX+UqwvsMvrUwxQXpiFAw2dqG9PUmAXdN4ON3Um/PFDUfbYpcpSrD9j5/R44XB7dfsFlbQR06fe5s2b8Z3vfAdz5sxBbW0tAODFF1/Eli1bVD04oljpPWMHBN78uVyjjZrjvmzd2KG5vTI5IrALXlJMRXZ/u5McS++/gXKl5Uni99k53IG2IkAKZey6Uitjl2sxQSTpUj0zTKkn6sBuzZo1WLJkCbKzs1FTU6PMYO3o6MBDDz2k+gESxUJk7PTY6kRgZay2Qu2vA4Ac/x67VM/YdfkzdrkhMnYAUJ+Eytj2PsFwqgR2qZaxMwS1XWEBBUUr6sDuf/7nf7By5Uo8++yzvbooz507NylTJ4hCEcUTWTpsdSLksZedpj4+2gag9/46AMj2f1dI9cCu05+xCy6eAIIrYxPfpFj0sMv2Z9KbOp0pcR5TLWMHsJcdxS7qT739+/eHnMlaUFCAtrY2NY6JKG7KHjsdL8XmWzl9Qitujxc7j7cB6B/Y5Yql2BQISMLpClE8AQQqY+vbHQk/JlE4UVGYhaH5vrFZR1Iga9fib3eSSg2n2cuOYhV1YFdRUdFrUoSwZcsWjB07VpWDIopXIGOn48BOZOy4FKu6ffUd6HZ5kJ9lwrg+BQbZ6VIVK4onBsjY1duSl7EryjGjakgugNRYjlWqYlNkKRbgdBmKXdSB3U033YQf/ehH+OCDDyBJEk6ePImXX34Zy5cvxy233KLFMRJFrccd6GOnV1yK1Y7YXzd9VDEMht6tJnKVPnYpHtiJ4ok+e+ySmbFr7/YFUEU5FoxNkcBOluWU62MHcCmWYhd1u5O77roL7e3tWLhwIXp6ejBv3jxYrVYsX74ct912mxbHSBQ1hyie0HXGjkuxWhHzYWf0KZwAAhm7VNgbFo5SPNGnKlZk7Jo6HXC6vbAk8MtPKmbsupweON2+94uU2mOXFehlRxSNqAI7j8eDLVu24Mc//jHuvfde7NmzB16vF5MnT0ZeXuL7IRENRGTssnScsRPtTjqZsVOdUhE7uqjfdaLdSbfLA4fbA2uKVl6Laum+VbEluRZYjAY4PV6csvWgsiQnYcckspxF2ZaUCexERazVZFCKOlKBWIplxo6iFdWnntFoxJIlS9De3o6cnBzMnDkTs2bNYlBHKacnAzJ2ovEsl2LV1WDrwYnWbkgScGZlUb/rs4xQeoylctauS5kV2/tvQJIkDCv0FS6cSnDLE1E8UZxjxtihgcAuljnjagmuiE2lCQ9iKZYZeYpW1OmMqVOn4tChQ1ocC5FqHBnQ7iSwFMvATk1iGXbCsHzlHAczSIFsaSoXUHSGmBUrVBRkAwDqEtykuC1oKbayJAcGyXecjZ2J3+8HAI0dDjy67nMAwJA8a1KOYSCB1xj/vik6UX/q/eIXv8Dy5cvxz3/+E3V1dbDZbL3+I0oFoio2VZfJ1MCqWG2IZdi+bU6CFfqzKalcQNEl+thZ+wd2wwqTM31CZMeKciywmowYWexbBj7cmPjl2PX7GnDRbzZh0+eNsJgMuO28cQk/hnBYPEGxirp44sILLwQAfO1rX+uVthaDij0ez0B3JUoYh1ssxeo3Y5fHkWKaGGjiRDBfYNed0kuxot1J35FiQHDLk8QGduJ8FeX4gpaqIbk41tKFw012nDO2NCHH0OPy4JF/7cPz7x8BAEwsz8dvvj0dE8rzE/L4kSpg1TvFKOrAbv369VocB5GqMqGPHd/41edwe/BZrW/lIVzGTlQspmpgJ8sy7P5Mbl6IjF2y5sW2dvUe3VU1JBcbP29MWAHFvnobfvTKTuw/1QEA+N6Xx+D/XjgxJd8n2KCYYhV1YDd//vwBr9u5c2c8x0KkmoyYPOF/4+90uJO6+VxPPqu1wenxojTXgtGlA1eLFqX4UqzD7YXX/5Lo28cOSN68WHG+xFK2KKA4pHFgJ8synn//CB7+1z443V4MybPi1986AwsnlGn6uPHgUizFKurArq/29na8/PLL+MMf/oBPPvmES7GUEjKh3YnIxHi8Mrpd/LtTw8dBjYnDVUgW+JvZpWrGrjNo32VOiC835UnYY9ft9ChbJMTorkS0PGnscGD5Xz7Bxs8bAQDnTSzDr755RsoVS/SVz+IJilHMgd27776L5557Dq+//jpGjx6Nyy+/HKtWrVLz2IhilglLsTkWI4wGCR6vzOVYlSiNicMswwKBjF2qBnaicCLHYuw3OQMI7LE7ZeuB1yuHvI3a2vxTJ0wGSWmaLAK7o812eLwyjCofx46jrbjpxe1o6nTCajLg3osn4drZo1OqrclAxFJst8sDl8cLs1G/X1JJXVEFdidOnMDzzz+P5557Dna7HVdccQVcLhfWrFmDyZMna3WMRFFTlmJ1nLGTJAl5VhPau10M7FQgyzK2R1ARCwSWyVI1sAsUToR+ix+aZ4VBAtxeGU12B8ryszQ/pla7KJwI9IsbXpgNi8kAp9uL2tZujAqz/B2t9fsb8F8v7UCPy4uJ5fl44qrp+NKw1CqQCEdk7ADfPtpUmopBqS3iT72lS5di8uTJ2LNnD5588kmcPHkSTz75pJbHRhSzQFWsfjN2QGA5tpMtT+J2orUbjR0OmAwSzhhZGPa2hSm+FBsonAj9+jcZDRia71uKTNRybJsyJzbQG9BgkFBVKvbZdar2WH+rqcX3X9iOHpcX8780FK/fMjetgjrA928kMpssoKBoRBzYrVu3DjfeeCMeeOABXHzxxTAa9f2BSenNkQFLsUByetm5PF54vPor1hDLsKcPLxj0dRPoY+fU/LhiYXeKpdiBF2VEZWyimhSLwoninN5Nn9XeZ/fclsO447WdcHtlXHrmcPzhuzPDnodUxibkFIuIA7vNmzejo6MDM2fOxDnnnIPf/va3aGxs1PLYiGLWkwGTJ4DAPpxEzItt73bhf9ftx5kPrMONL3yk+eMl2sfKfNjwy7BAILBL1YxdlzInduAAtTxon10iBCpiey8pjlEpsJNlGb/+9z78/J97AAD/Z+4YPH7FmWm9N00U6bAylqIR8St+zpw5ePbZZ1FXV4ebbroJr776KkaMGAGv14vq6mp0dHRoeZxEUenJkKVYkbHTcim2y+nG0xsOYt6v1uPJdw/A7vTg/YPNumuxsiPCwgkg9QM7ZZxYiFYnQkVhYseKtQbNiQ02VoXAzuOVcc8bu/C79QcBAP+9ZALuv2RyQopCtMRedhSLqL/K5OTk4Prrr8eWLVuwa9cu/PjHP8YjjzyCsrIyfO1rX9PiGImiFhgplr7f1iORp2GTYqfbiz9uPYL5v96AX769D+3dLowrywPg28Oot1FmR5u7APgmEQwmOLBLxQC3y78UG2pOrDDMvxR7KkGBXd+pE0KV6GUX41ixHpcHt778MV758DgMEvDwN6bi1oXj0qLydTDsZUexiOtTb8KECfjVr36FEydO4JVXXlHrmIjiIstyRrQ7AYL22KkY2Lk9Xvxl+3Gc978bcN/fd6Oxw4FRJTl4/Mpp+Pcd85DvzwI12JIzuF0LXq+sZLn6LhWGIqZ+uDyp2UOwM4KlWNHyJGEZO3tgTmwwscfuZHu38ncbqY4eF763+iO8vbseFqMBT11zFq6aNUqdA04BWvx9k/6psqPUaDTisssuw2WXXabGryOKi9srK133s0x6D+wC0yfUsG53PX759j4c9GdPyvKtuP388bhiZiUs/uzn0AIrOhrdaOxwKBm8dNfhcEMk3oLbTAwkx2KE2SjB5ZHR1uVKuc35XYO0OwESP32ibYCMXWmuBflZJnT0uHG0uSvima22HheuemYbdp+0Ic9qwjPXzcDc04aoftzJxKVYikVqvRsRqSD4W79V58UTot1Jh8MNxNlI/4tTHfjBizsA+D58b1lwGq6bM6Zf1nNonhWHGu1o7NRPxk58cFpNhoiyvJIkoTDbjKZOJ9q7XRhelK31IUbF7m9QHLZ4ImherCzLmi9dtvWZEytIkoSxQ3LxyYl2HG7qjDiwe/XDY9h90obSXAteuH4WpowI36ImHQWKJ5ixo8gxsCPdEc2JJUn/e+wKgpdq4gzsRJZuXFke3rhlrpIN7Ev0P2vs0E9gJ5a6xJ6mSIjALhXnxdojKJ4QGbtulwe2bjcKcyJ/7rEQ56koxDmuUgK7roh/37rdpwAAP7pgvC6DOiCQkeceO4qGvj/1KCMFF07oYQN1OGouxda3dwMAxpflDRjUAfoM7MQHZ0EEy7BCKlfGRlI8kWU2KsuidbZuzY+ptSsweaKvqiG+Jf3DETYpbu50KFXMF0waptIRpp7AUiwzdhQ5BnakOw63COz0vb8OCJo8ocJSTb2/GEJkcgaiy8DOH5xFm7ELvm8qiaTdCdB7OVZLsiyjPcTkCUFUxkba8uSdfQ2QZWDKiIKUWwZXE/vYUSwY2JHuiKVYvTcnBtStmhMZu4rBArs8X2DX0JGYTfeJIPYwFYTJVPYlMk9iVFYqEcUTYiTVQMS/tdaBnd3pgcvjq07pu8cOiL6XXfUe3zKsnrN1AIsnKDb6/+SjjCMydnpvdQKouxQr2l6UF4bPgJT5szzM2KXuUqwonsgZLGOXoMpYUThhMRlCfuES0ydEMUo4PS4PNn/hm3q0aLK+Azu2O6FYMLAj3VEydhmwFKvmrFjx4S6W5wYiMnZNeqqK9S91RdLqRAjMi03BwM6fscsLUxULAOUFviBe64xd8JzYUPte86wmlPmX+I8MkrXb8kUTelxejCjKxuSKAvUPNoWwQTHFgoEd6U6mzIkFAoFIj8sLjzf23yPLspKxG3Qp1v8B3Gx3wh3Pg6YQsTk9mqXYtMjYDdJfr7zQ92+pdZPiQEXswM2fqyJcjg0sw5bpvjiqICgj7/Wm3oQTSk36/+SjjCMydtYMWIrNC1pq64ljAEJrlwtO/3zdsoLwfVNKci0wSIAsAy321NtfFgulKjZbH1Wxot1J3qBLsb6M3SmNl2LFnNhQhRPCWDFaLExg5/HKeGefL7BbNLlcxSNMTeKLmywDnU4ux1JkGNiR7mTKnFgAMBkNyPYHsN1xBHZ1/sKJIXmWQauJjQYJpUoBhT6WY5U9dlEVT6RmYOfxBsac5URYPKF5xm6AqRPBIsnY7TzehqZOJ/KzTDhnbIm6B5mCssxGZeILCygoUvr/5KOM43CLqlj9Z+yAoOXYOAK7eqVwIvwyrCD22emlgCKQsUv/pdjg2bWDtTsZ5t9P2d7tQrdTu5m3bfbQUyeCRdLLTizDLpxQBrMxMz6+2MuOopUZfxmUUQJ77DIssIvjfV+piC2IrCeYWK7VS2CnTJ6IonhCZJ9SrXhCLMMaDdKgWeuCLJOS1dOyMlZk7MJNt1Aydo12yHLo/WTVe+oBABfovBo2GHvZUbQY2JHu9Ih2JxmwFAsAef5v9N2e2DeS10dYOCEoGTudVMbGkrELrlhMpY3tIrDLsRgHLS6QJEnJ0orleC20DjAnNtiokhwYJF/Pu1BfGA41duJgox1mo4QFE4ZqdqypRrQ0YssTilRmfPJRRgk0KM6MjF2BCkuxddEuxeps+kQ8VbGyrE67GbWIithw48SCifY2WhZQtIeZEytYTAZUluQACF1A8Z+9vmXY2WNLo/p3Snfi75t77ChSDOxIdxwZVDwBBJZi49mCIz7UI87Y6Siw83pldMRQFWs1GZXClfYUWo4VPexyB+lhJ5QnoIAiUBU7cMYOCF9AkSnTJvpiLzuKVmZ88lFGybQ9dqKlRXwZO98yXCZm7OxON8RKarSZoFQsoFDGiQ1SOCEkYl5sJFWxwMCBXXOnAzuOtgLIrP11QHDGLnWywpTaGNiR7gSqYjPj5S324PTEuMeud3PiyIon9DQvVsyJtRgNUWd5lQKKFJoX2xnlUmwi5sUGJk+Ez9iJmbGHGnsHdu/ua4BXBk4fXoARRZG9RvWiQNljlzpfHii1ZcYnH2WUTMvYKUuxMWbsOhxudPlbXQw2TkzQU8YuMCfWFPUkg4JUzNg5ol2K9Y8V02iPndcrK7NiB8/YhW55kqnLsACXYil6DOxIdzJp8gQQlLGLcaVGZGqKcszIHqShrVDmDwDtTo9ShZmuYmlOLKTiUmynUhUb3VKsVnvsOhyBpe7CQaqOq/zTJ461dCnj6npcHmz+ogkAsCjDlmEBLsVS9BjYke5kWruT/Dj32AV62EWWrQOAXEugcKApzVueiKXY/ChanQiiyjOVetmJ7GvEe+z8S7FNnQ64NJj9KwpLss3GQbPoFQVZsJoMcHlk1Ppfl+8fakG3y4PhhVk4fXiB6seX6pR2J47UeY1Rakv6J9+mTZtwySWXYPjw4ZAkCX/7298Gvc/GjRsxY8YMZGVlYezYsVi5cqX2B0ppQxkpljEZO7EUG9seu/ooCycAX/8zvSzHBjJ2kVfECiIDlUqtKJSq2Aizr6W5FpiNEmRZmxFxgR52gwfOBoOkFFAc8RdQvLO3AYCvaCLapXI9UBoUM2NHEUp6YGe32zFt2jT89re/jej2hw8fxtKlS3HuueeipqYG99xzD26//XasWbNG4yOldKEUT2RIxq7Cv5m8JcaVtLoomxMLegnsOmJoTiyk4lKs3RFdVazBIKEsX7sCisDUifCFE4JSGdvcBa8MvLu/EUBmLsMCQSPFuMeOIhT9V1SVXXTRRbjooosivv3KlSsxatQorFixAgAwadIkbN++HY8++iguv/xyjY6S0kmmNSiuKvV9ELa7JHQ53Sg0Rxeg1Ec5TkzQy/QJW0/0zYmFVBwr1iWqYiMsngB8QX1tW7c2gV0UGTsgENgdbe5Clx1o6nQi32rCOVWlqh9bOihIwawwpba0S2ls3boVixcv7nXZkiVLsH37drhcfOFToEFxpgR2hTlmZa/X0ebox0LFm7FrsKV5YBdUFRutVKyKFUuxkRZPAMAwDceKiaB3sIpYIdDLrgu7WnwfUfMnDIUlQzLwfYmtFh097gFn6BIFS3rGLlr19fUYNqx3Sn7YsGFwu91oampCRUVFv/s4HA44HIEPH5vNBgBwuVwMBjUkzm2iz3G3P7AzSd6M+fcdVZKNtloXDjbYMKkiP6r71rX5PsyH5JmiOl+lub4P6lO27rQ+zyKjlGc2DPo8+r6m8y2+YKO1y5ky50AsLWeZpIiPqTzft0x6vMWu+vNo9vc6LMiK7PU1qtgXZB5ussPj9O2pO2/CkJQ5v4mW7f9+6vbKsHX1RBWwRyJZ79OZKFHnOO0COwD9NtCKbzEDbax9+OGH8cADD/S7fP369cjJyVH/AKmX6urqhD6ezW4EIOGDre/heIb881ocBgAGvPPhpzDUfhLVfU80+87XFzs/RMfnkd+v/pQEwIjdB49j7dqjUT1mKvn8iO/cHT24H2vt+yK6j3hNH+0AABNOtdiwdu1azY4xGidP+f499+3aCXNtTUT3sdX7/i137DuCtTik6vHsPOw7v00nj2Ht2iOD3r7TBQAmnGzvASDBIMlwHKnB2hORPRe9kWXAACO8kPC3t9ahyKrN4yT6fToTdXV1JeRx0i6wKy8vR319fa/LGhoaYDKZUFoaeg/G3XffjWXLlik/22w2VFZWYuHChQPeh+LncrlQXV2NRYsWwRzlvq943L3jHQAeLDpvAUaXZEZk98V/Psf2jUdgKRmJpUunRny/LqcbXVvfBQB866uLlWWfSGTtb8Srh2pgyCnC0qWzoz7mVPHqqe1ASwvmzDgTS6f1z/gH6/uaPtJsx2OfvQenZMLSpUsSdMTh/fbge0CnHfPmzsKcsZG9v+V/0YS/HP4YTnMBli6dq+rxvPOXXUB9HWZOnYilXx4z6O1lWcavdq9Hu78KdNaYEnzza2erekzp5mefrEdbtwtnf3kexpflqfq7k/U+nYmam5sT8jhpF9jNmTMH//jHP3pdtm7dOsycOXPAF6XVaoXV2v9rjtls5gs5ARJ9nkVVbH62NWP+fceW+ZZfj7V2R/Wcm9t8WxTyrCaU5EdXPFFR5AuamzqdaX2eO/xVpCV5WRE/D/GaHpLvOwd2hwcwGGE2Jn8fWJfT9/ovyIn8+YwZ6nv9HG/thskU/QSOcGwxnN+qIXnYebwNgK8aNp1fX2ooyDajrduFbres2bng56H2EnV+k/4u1NnZiZ07d2Lnzp0AfO1Mdu7ciWPHjgHwZduuu+465fY333wzjh49imXLlmHv3r147rnnsGrVKixfvjwZh08pxuXxwuNvc59lyoziCQAYU+oLMI42R5fqVypioyycAALFE02dDni96bupW/QHi6d4wvd7UmOPUleUfewAYERxNiTJtz+1qVPdubetEc6JDSZmxgLA+ROHqno86Yi97CgaSQ/stm/fjunTp2P69OkAgGXLlmH69Om47777AAB1dXVKkAcAVVVVWLt2LTZs2IAzzzwTDz74IJ544gm2OiEAgebEAGA1J/3lnTBiybmx06mMlIpErBWxAFCa6wvs3F5ZaUKbjkR/sFjanRgNkrJ83ZYigZ3dEd3kCQCwmozK5JHjreruA2qPcE5sMFEZOyJHxoii6DLJesRedhSNpC/FLliwIGwJ9/PPP9/vsvnz5+Pjjz/W8KgoXYkedgBgzaD2CAXZZuSaZNjdEo402TFlRGFE9xOD36MZJyZYTAaU5FrQYneisdOB0jyNdnVrSJZldIg+djE0KAZ8TYo7etwp0fLE6fbC6R8Llhtl9WRlSQ7q2ntwvKULZ40qVu2YAhm7yM/vZdNHYMP+BpyVk5g9SalOfHmwxToQmjJK5nzyUUZQxomZDBk3fmioPzY70myP+D51MYwT6/WYeek9faLL6VGW7mPJ2AGpNX1CLMMCQE4UDYoBoLLYl/U93qJexs7jlZUsU2F25EuxlSU5ePX7szClOH2X+NWkZOxS4DVGqY+BHemKMk4sQ5oTBxua7fsQFDM2IxHPHjsg/ceKiaDDbJSQFePSvVhibE+B6RN2p++LjcVkiLqQY5R/Of+YioGdrdsFsSATzVIs9aZMn+BSLEWAgR3pSo8ydSLzXtpDs/yBXRQFFPHssQN0ENh1B8aJxZrhTamMnSP6wglhVKlvL9vxFvWmT4h9h3lWU0pUDKer4OkTRIPhXxrpisOdWePEgilLsbFk7KKcE6s8ZroHdqJwIsb9dUBgiTEV5sWKwploCicEsRSrZsauNYbCCeqPS7EUDQZ2pCuieCKTWp0IgYxdZIFdj8uDZrvvgzfmjJ1/j11DugZ2/g/KaBoz95VSGTv/Umy0hRNAYCm2rr0bLo93kFtHpj3KObEUWmAplhk7GhwDO9KVzF6K9f1vU6dTmRcaToPNF4xZTYaYP3jLCnSSsYuxcAJIrcBOZOyiLZwAfNlXq8kArwycbFNnOVZk7KLpYUf9FYiq2BR4jVHqy7xPP9I1kbGzZmDGLssElOb6PkCPNA2+nCZanVQUZsW8v0ypiu1M08AujubEglI80Z38Xn6iKjYvhqVYSZJQqXIBhWh1UhjHUjcB+f4vHpF8YSNiYEe6IvbYZVJz4mBiAsXhCJZj4211Auhgj123vjJ2ojlxTgzFEwBQWaxuAUU7M3aqUCZPcCmWIpCZn36kW8oeuwwsngCA0f7ALpICinqlIjb2zv4isGvvdilBdTpRo3iiyH/fVCiesMdRPAGo3/KklXvsVMHiCYoGAzvSlcAeu8wM7ETGLpICiro4e9gBvmyV2ehbxlV7xmgiBNqdxL4UW5BKGbs4iicAKEuxao0VE+1Oipixi4t4jTnc3rT8AkWJxcCOdKVHtDvJoHFiwcbElLGLPbCTJClQGevfsxcJWZbx4D/34LHqz2N+bDV0ONRod5JCgV0cxRNAUGCnUsauTbQ74R67uATvmWQvOxpMZn76kW4pxRMZusdOLKVF0qS4Lo45scGG+u8fzT67g412rNpyGE+88wVORREQqi24QXGsxDKjw+1VMsbJohRPxJixG6V6YOefE5vLwC4eRoOEfCsrYykymfnpR7rlEEuxGVgVCwT22LXYnYNmkOr9xRPx7LEDYquM3VdvU/7/tkPJG/Qe2GMX+1JsntUEo8G3HJ3srJ1SPBHjHjuRsWvtcqkyvkq0O4lmTiyFxl52FCkGdqQrmTwrFvAFGaKgIdxyrMvjVZoKDyu0xvWYsVTG7q0LDuxa4nr8eKhRFStJkrIcm+wCCnscI8UA3+unxN8yR42snWhQXMziibgFxooxY0fhMbAjXcnkBsVCVWkugPAFFI0dDsgyYDJIGJKbjMCuQ/n/HyQ1Y+cLhPLjCOyA1NlnZ3fGVxULqNfyxOXxosMfaLJ4In6Bylhm7Ci8zP30I13K9KpYABgzxN/LLkzGTlTEDivIgsEQW3NiId6M3aEme1SFF2qRZTmQsYtjKdZ3/xQJ7PxLsbkxFk8A6hVQBJ8LNiiOX6CXHTN2FB4DO9KVQPFEJgd2vozd0TAFFGpUxArRzott63IqgaWo4t12OPHLsd0uD9xe33zdeJZigeBedslt+aJk7GIsngDUa3kizkVBVmAPIsWugNMnKEIM7EhXRLsTa4a2OwGAMf6l2HAZOzFOLJ4edkK0GTuxDFtZko3zJw0DkJwCCrGkZTRIMU9qEFJlKbZLydjFHtip1aQ4UBHLZVg15CvzYrkUS+Fl7qcf6ZIjwydPAIHALtweu0BFbPyBXVl+oCpWluVBby+WYSeWF2D22FIASQrsREVslinmWblCYF5sspdi/X3s4ghU1Wp5okyd4DKsKgJVsczYUXgM7EhXMr1BMRDYY9fW5RpwaTAwdSK+VidAIGPndHsjasUgWp1MqijArDElkCTgUGPi99kF9tfFH3ikQsZOlmVlKTYvruIJsRTbDa938EB9IEpzYhZOqIJjxShSmfvpR7qU6bNiASDHYsKwAl+wNdByrJp77LLMRmWZKJLlWLEUO7kiH4U5ZkwqLwCQ+H12ooN/vPvrgNQI7HpcXog4LNY+dgBQUZQFo0GC0+2NeN9kKG2cE6uqQLsTLsVSeAzsSFccrIoFMPhyrBpzYoNFus/O7fFi/ylfYDepwhfQieXYRLc9UaM5sZAKfexEtg4AcuJ4/ZuNBiXgj6eAoq3bl7ErZsZOFVyKpUgxsCNdEe1OMrl4AgCq/JWxR5r6fzB7vbIyxkuNjB0QXBkbfjn1cJMdTrcXuRajsuQ3e2wJgMTvs1OjObGQChk7UTiRYzHG3cJGKaCIYDTdQMQeO7Y6UQf72FGkMvvTj3SnJ8MnTwijw2TsmuwOuL0yDFIgIItXpBm7vfW+bN2E8nwl+JhV5dtnd7DRPmhgqCabikuxYh9ZMgO7TqVwIv4MZGCfXeyBHadOqIt97ChSDOxIVxycPAEAqPIXUIQaKyb215XlZ8FkVOc8leX7Mn+DzYsVFbFiGRbwBUVin90HCRwvJjJ2Yu9SPFIiY6cUTsT/pWZUafwtT1pZPKGqfKWPHTN2FF5mf/qR7jBj5yOaFB9usvdrQaL2/jogioydaHUSFNgBwDlJWI4N7LFTdyk2kpYvWrA7xVKsChk7/1LsiTjGirF4Ql0F/i8gnQ433B5vko+GUhkDO9INl8cLj78sMMuU2YHd6BJfYGfrcSt7nQSRsSsvSF5gN7kiv9flSgFFAitjxV6lAhUydiJ48XhlZUk00UQPu3hanQhiXmw8GTu2O1FX8DzjZL3GKD0wsCPdEIUTAGDN8KXYbItRKYzo2/IkWRm7FrsTp2y+6yeU98nY+ffZHWjojGrmbDzUzNhlmY2w+At2krUcqzQnVmMp1p+xq7f19Pq7ikZbN/fYqcliMihbTLgcS+Fk9qcf6YroYQewKhYIannSJ7BTc+qEIIowwgVl+/zZutGlOf2ySkU5FkwU++wOJ2Y5Vs2qWCB4XmxyA7t45sQKJbkWZXpFbVv0y7EOtwdd/qXhomxm7NQiXqvJnnBCqY2ffqQbwa1O4h0RpQdin93RPpWxas6JFUTGrqXLCdcA+3/2KKPE8kNen+i2J0pVrErtOMQ+u2RNBhB77HJVyNhJkhTXzFhREWuQ1ClOIR/2sqNIMLAj3XCwcKKXMf7KxsN9epEFpk7EP05MKMm1wGiQIMu+JddQxMSJSX0KJ4RzqkSj4sTss+tQsUExkPzKWFEVq0bxBACMLBYFFNEHdsE97OLtqUcBYj8oe9lROAzsSDd62OqklzFD+i/FyrKs7LFTcynWaJBQmutbchtoOTZ4Rmwo51T5MnZfNHSiaZC2KfGSZTmoeEKlpVj/XrK2pO2xUy9jByCujJ0onODUCXUFWp4wY0cD4ycg6YbDzXFiwaqCAjvRgqOty6VkNssK1GlOLIQroHB5vPjiVCcAYPIAgV1xrkVZptU6a+dwe+H0LxmrtRRbkOSMnbLHToWqWAAYVeLL6B6PoeWJkrFj4YSqAkuxzNjRwBjYkW6I4olMb3UijCrJgSQBHQ43mv3LoyJbNyTPAqvK5ylcYHeo0Q6nx4s8qwkjigZeAhZtT7TeZyf2wRkkINeiznkQRQJJC+yc6hVPAIFedjHtseOcWE0ElmKZsaOBMbAj3VCKJ7gUC8CXuRzu30cnlmPrbb7si5qFE0K4ebFiGXZi0CixUBIW2PWIqRNm1QptCpNeFSuWYtXK2PnHirV0Rd10WWTsijgnVlUiY8d2JxQOPwFJN5TiCWbsFGP8o8VELzulh12BeoUTQriM3Z4Qo8RCmZWgfXbtYn+dSoUTAFCYndxsSpeSsVPn9S+KJzoc7qizkIGpE8zYqUlUGLMqlsJhYEe6wYxdf6NLRcsT33JavQaFE0KZCOxCBGSiInZiRehWJ0JJ0D67DzWcQqE0J1apcAIIBDFt3aGrgrXW6c/Y5aiUscu2GJVgPdrl2MDUCWbs1CRer1yKpXD4CUi6oeyxY/GEosof2B1u7pOx02IpNt/3O0Nl7PZGmLEDErMcq3ZzYiB12p3kqVQVCwRGi0VbQCEydpw6oS72saNIMLAj3Qi0O2FgJ/RteaJlxm6gpdimTgcaOxyQpIGbEwdLRKPiQHNiFZdic1KjKlatPnZA7C1PWv0Zu0IuxapKLMVyjx2Fw8COdKPHHZg8QT5V/j12ouVJnX+cWHlB4gK7ff5l2DGluREFHbP8jYo/P9WJZo322WmZsUt28UTfcW3xUAooWqML7No5J1YTylIsM3YUBj8BSTcCS7F8WQuVJTkwSL5xU42dDpyy+QIlbZZifYGd3elRskdAYBk2kmwdkJh9dh0qjxMDAoFdR48bHm90VaTx8nhldPsz1jkqFU8AwMigythoiIwd58SqK1Cgw4wdDYyfgKQbSoNiVsUqrCYjhvv7xn1W245Of8ClRWCXazEi278MHpy1i2Z/nSCmUGi1HKtF8URhUJCY6M3tYn8doF67E6B3y5NoBKpimbFTU0HQ5IloW9BQ5mBgR7rhYPFESGP8BRRbD/qCpMJss6r7sARJkpRpFsGVsXvrw8+IDSVQQKFNxk5ZilVxj53ZaFBajSR6n12X0/elxmiQVN2KIJoUn2jtjjgL2e30KK2HGNipS4wU88q+zDhRKAzsSDc4KzY00ctOBElaFE4IokmxyNg53V4caPC3OolwKRYI9LPbf6oDLXb124coxRMqZuyAQMuTRAd2nUrhhFG1hsuAby+m2SjB7Q3szxyMaPdiMkiq7vcj33ub2ej792XLExoIPwFJN1gVG5rI2O0+2Q5Am2VYoW8BxcHGTrg8MvKzTBhZHHlT5NI8KyYME/vs1F+OFR+KospQLWLPXluiM3YaFE4AvgygGAEXacuTVntgGVbNIJN8WXEWUNBgGNiRbojiCVbF9lblb3kiVtI0zdj1CeyU/XXlBVF/yAfanqi/HKvssVN55JXY3J7ojJ2YE6tm4YRQGeU+O5Gx49QJbbDlCQ2Gn4CkG0q7E2bsehG97AQtxokJfZdi9yn76yJfhhXO0bBRsagqVH0p1l8F2t6V2OkTogpZzcIJIdqWJ22cE6sppUkxl2JpAAzsSDdYPBFaZbGv5YmQiIxdQ4evEbLS6iSKwglBVMbuq1d/n10gY6duIJSs6RNiI32uBkUxlVE2KeacWG1xKZYGkxKB3VNPPYWqqipkZWVhxowZ2Lx584C33bBhAyRJ6vffvn37EnjElIp6lHYnKfGyThkWk0EZ6A5ou8eub1VsLK1OhNI8K740LA8A8IGKWbselwdOf9Wm2kuxRUmaPhHI2Kn/pSbalietnBOrKfFlhEuxNJCkfwK+9tpruOOOO3DvvfeipqYG5557Li666CIcO3Ys7P3279+Puro65b/x48cn6IgpVXFW7MBGlwYCO22rYgPzYhs7HGjqdEKSoBRCRGuOfzl2q4qBnch0SBKQp3KGqyBJ0ye0GCcmVBaLjF1kxROcOqGtfCuXYim8pAd2jz32GG644QbceOONmDRpElasWIHKyko8/fTTYe9XVlaG8vJy5T+jkR/mmc7h4kixgVQF7bNLRFVsU6dTqcKtKs1Fdoyb+uecpv4+O5HpyLeaYDCoW7WZrKVY0cdOyz12TZ2OXo2QB9JqZ/GElkTGzsaMHQ0gqZ+ATqcTO3bswOLFi3tdvnjxYrz//vth7zt9+nRUVFTg/PPPx/r167U8TEoTbHcyMNHyJM9qUpqcaqE0z/dh7vHKSkPkWJZhhXOC5sY2qTQ3NtCcWP3zIJYfE93uRFmK1aAqtjDHrFRinmgdPGsnnjuXYrWh7LFjxo4GkNTukU1NTfB4PBg2bFivy4cNG4b6+vqQ96moqMAzzzyDGTNmwOFw4MUXX8T555+PDRs2YN68eSHv43A44HAEPhRsNt++H5fLBZeLfxxaEec2UedY7LEzSd6M+3cd7FyPKfFl6YYXZml+bopzzGjtcmHj/gYAwJfKcmN+zDyLhInD8rDvVCfe+7wBS6eWx318rZ2+wo58qynq4xrsPOf6m2O3dzkT+hrs6PFlybLNkiaPW1mcjT11HTjUYENVSfiMb6vd916bbzHEdSyJfv9IF7kWdV9jPM+Jk6hznBJtwfv2t5JlecCeVxMmTMCECROUn+fMmYPjx4/j0UcfHTCwe/jhh/HAAw/0u3z9+vXIyckJcQ9SU3V1dUIex95tBCBh23ubcUC71caUNtC59srA4hEGjCtsx9q1azU9hizZ9++w71QnAKCzdj/Wro29uGmYZMA+GPDnjTuB4964j+/jJgmAEa6u2M/FQOf5eCcAmHCqtUPz8xzsi0MGAAYcO/QF1vZ8rvrvNzl8v3/dezvgOBR+tNiJBt+///5dNZCPxT/PNFHvH+niUKPv9XvoRB3Wrq1V7ffyPGuvqyu6mcuxSmpgN2TIEBiNxn7ZuYaGhn5ZvHBmz56Nl156acDr7777bixbtkz52WazobKyEgsXLkRpaWn0B04RcblcqK6uxqJFi2A2a78ss+yDagAyllxwHoYVZFZkF8m5/mqCjuW1hu2oOxhoKvydixdgeFHsvfMsexuw8U87UefJx9KlX477+No/Og58sRdjhg/D0qXTo7rvYOf5WEsXHt21BQ7ZiKVLl8R9rJH65592Ak0NOOuMKVg6q1L13/+pYT8+fe8o8iuqsHTpxLC3fXDXBgBOLF7wZUyOYxk+0e8f6cK6rwEvHdgJa14Rli6dHffv43lOnOZm9XtyhpLUwM5isWDGjBmorq7G17/+deXy6upqXHrppRH/npqaGlRUVAx4vdVqhdVq7Xe52WzmCzkBEnGeXR6vMqQ8PzsrY/9dU+E1HdwAuSDLhFFD8uMaLTV3XBkkCTjUZEdrtwdlcQbtdqfvdVKYY435XA10nofk+1YAul1eeCUDrKbY97x5vTKa7A6U5Q/+fLv9FeHxPKdwxgz1VTXXtjnC/n5ZlpXCkaEFOaocSyq8plNJSZ7v76vT4VH1vPA8ay9R5zfpS7HLli3Dtddei5kzZ2LOnDl45plncOzYMdx8880AfNm22tpa/PGPfwQArFixAmPGjMHpp58Op9OJl156CWvWrMGaNWuS+TQoyUThBABYzayKTSZRGQv4CifinRdamGPG5IoC7D5pw9ZDzbj0zBFx/T6tmhMDvnFPkgTIsq8ytiw/tsDO7fHihhe2Y+PnjXj9lrk4a1Rx2Nt3OrQbKQb49tgBg/eyszs9cHl8gTOLJ7QhClnYoJgGkvTA7sorr0RzczN+/vOfo66uDlOmTMHatWsxevRoAEBdXV2vnnZOpxPLly9HbW0tsrOzcfrpp+Ott97C0qVLk/UUKAWIHnYA250kW9/ATg1zxpZi90kbth1qiT+wE1WxGlQHGwy+Ie3t3S7Yul0RZdtCeeRf+7Dx80YAQPWeU4MGdqINSZ4G7U6A3mPFwu2BbvM3J7aYDMhmdbomAiPF3GH/LShzJT2wA4BbbrkFt9xyS8jrnn/++V4/33XXXbjrrrsScFSUThzuQA87vtElV+/ALrbGxH3NHluKP2w5rEo/O9H/S4t2J4AvU9Xe7cKhRjvGlUX//P++sxZ/2HJY+XnHkdZB72N3+F7/ORoFdiOKsyFJvn55zXYnhuT139oC9J4Ty79DbRT4M3ZOjxcOt5ftnagfpjZIFzh1InUMzVM/YzdrbAkMEnC4yY769p64flcgY6dNEDRjtC+79uM/f4Kdx9uiuu9nte2466+fAgAumTYcAPDJiTZlBNpA7E7t+tgBgNVkRLl/b2O4mbEisCtmc2LN5FpMyuxnLsdSKAzsSBcCzYn5kk42MS/WIAFfinGUWF8FWWZMGVEIIP4pFB3KHjttMnYPXjoFs8aUoMPhxrWrPsAnEQZ3zZ0O3PTiDjjcXiyYMBQrrjwTJbkWONxefOaf4jGQLod2kycEMVos3D47MSe2kPvrNGMwSMqSu62b0yeoP34Kki6IpVhm7JLvtKF5+OaMkbjzgi+p+u8xW8yNPRhfYKcsxWo0gSPXasLq753tC+563PhOBMGd2+PFbX+qQW1bN6qG5OI3354Oo0FS9tZtP9Iy4H2dbi+cHl9GL1eDWbFCpX+f3Z2v7cTsh97B1596D7e+/DH+5597sGrLYazdVYddtb4AlHNitaXss2PGjkJIiT12RPESS7EsnEg+SZLw6Lemqf5754wtxTObDmFrnBm7wEgx7d7+RHD3f1Z/iI+OtOI7qz7AyzeegzNGFoW8/UNr92HroWbkWox45toZyszZmWOK8Z+9p7D9SCt+ELr/eq/5rTlW7b7YLJ1ajrW76tDt8qDe1oN6Ww9q0BbytkXZXIrVku9LSTfHilFIDOxIFzgnVv/OriqB0SDhWEsXatu6MSLGpsdKuxMNZ+YCIribhf/z3IfYfrQV3/nDB3j5xtmYOrKw1+1e//gEnnvPVyzxv1ecifFBy9dnj/Fl7HYcbR2wAtLu9L32LSYDzEbtvticP2kYPntgCZo7HTjZ3oP69m6cbOtBXXs36tp7fP+1daPL5cGiyZE3mKfoiZYnHT1ciqX+GNiRLjj8m8uz4mgIS6ktz2rC1BGF2Hm8DdsONuPyGSOj/h0Ot0fJ7mod2AG+Y37++qDgzp+5E/sFd51ox92v7wIA3H7eOFw4pfcs3CkjCmExGdBsd+JIcxeqhuT2ewy7Q9vCiWBGg4Sygixfk+jKIs0fj0LjUiyFw3Ur0gWRsWNzYn1T9tnFuBwbnOHI06gqti8R3M0YXYz2bheu+cMH+Ky2HU2dDtz04nY43F6cP7EMd1zwpX73tZqMOMMfBA60z04J7DQsnKDUIr6UsHiCQuGnIOkC251khjmnxVdAIfYk5VtNMBoS12ctz2rC8987G2eNKlKCuxte2I6T7T0YOyQXj3/7TBgGOJ4ZQcuxoYgedloWTlBqEftDO5ixoxAY2JEuKBk7Fk/o2szRxTAZJNS2dQ863ioUrZsTh5OfZcYL18/CdH9w98nxNuRZTXjmuhlhl4Vnji4BAGwfKLDzF09oWThBqSU/i0uxNDB+CpIu9LDdSUbItZpwhr/4IJblWCVjl6Bl2L5EcHf2mGJYTAY8fuWZg06nEA2PDzR0KiO7gmk9ToxSj2iuzaVYCoWBHemCQ1mK5Uta78Ry7LYYlmNtGjcnjkRBlhl/vmkOPrr3goiqR0tyLRg71Fc0EWo5tlOME0tA8QSlBhZPUDj8FCRdUDJ2rIrVvTljhwDwTaCQZTmq+3Zo3Jw4UpIkKb3qIjHTn7X7KMTc2C4WT2ScArY7oTAY2JEuOFg8kTFmjC6G2SjhZHtP2LmloSSiObEWZo7x7bPbcbR/ZWyg3Ul6PSeKXaAqlhk76o+BHemC2GfEpVj9y7YYcaa/h1q01bGJak6sNpGx++REuzI+TxANilk8kTm4FEvh8FOQdGFffQeAwDxL0rc5/n5226IsoBCbzZO5xy4WVUNyUZprgdPtxWe1tl7XKcUTzNhlDPHFhEuxFAoDO0p7XU43dp/0fdid7V+yIn2bfVqgUXE0++wCGbv0CoIkScJZo0U/u97LsUrxBPfYZYyiXF9g1+X0oL2LWTvqjYEdpb2dx9rg8coYXpiF4THOD6X0ctaoYliMBpyyOXC4yR7x/ZQ9dmm2FAsElmO39ymgEMUTeVyKzRgFWWaM9Y+X2x5i3yVlNgZ2lPZE49aZzNZljCyzEdNHFQGIrp9doEFx+mW3ZgZNoAjOUnb6A7scLsVmlFlVvve7Dw4zsKPeGNhR2vvIP0PzbP8HH2UGpZ/docg/2NI5YzdlRCEsJgOa7U4caQ5UA3f5iydymbHLKOeMZWBHoTGwo7Tm8cqoOdYGAJgxmhm7TDJ7bGBubKT77FKhQXGsrCYjzhjhm7qx/Ujgw1yMFGO7k8wyq8r3+v+stl1peUMEMLCjNLev3oZOhxv5VhMmlIcfzUT6Mn1UEawmA5o6HTjY2BnRfZSq2DTM2AGB7QbB++zsbFCckUYUZWNEUTY8XhkfHws9R5gyEwM7SmviA2766GIYDVKSj4YSyWoyKnNUt0awHOvyeNHt8i1bpuMeOyCogCJow3wXR4plrHP8++w+5HIsBWFgR2lNFE6cPZr76zKRWI6NZG5scM+vvDTNbolA9mCjHa12J2RZVpZi0/U5UexYQEGhMLCjtCb2Gs1g4URGChRQDL7PThRO5FlNMBnT862vONeC04b62lzsONqKHpcXXv/TZh+7zCMCu53H29Dj8gxya8oU6fnuRgSgtq0bde09MBkkZcQUZZZpI4uQYzGi2e7Ertr2sLdN1+bEfc30FwltP9qqZOsAIIdzkjNO1ZBcDMmzwun24tMT4V//lDkY2FHaEtm604cXsIdXhrKYDFgwYSgA4N+768PeVhRO5Kdp4YQwY0xgAoVd6WFnhIF7TDOOJElB++yiG69H+sXAjtKW6F/HxsSZbcnp5QCAf+8+FfZ2gVYn6f0lQBRQfHKiHa3+cVL8YpO5uM+O+mJgR2lLVMSyMXFmWzixDGajhAMNnTjQMHDbk3RuThysakguSnMtcLq9SpaG48Qyl2hUvONoK1web5KPhlIBAztKS+3dLuw/1QGAjYkzXUGWGXNPGwIg/HJsOjcnDiZJklIdu+nzJgDM2GWyL5XlozDbjC6nB7tP2pJ9OJQCGNhRWqo51gpZBsaU5mBovjXZh0NJduEUsRwbJrBTmhOnfxAk5saK/mVsdZK5DAYJZ4/hPjsKYGBHaUkswzJbRwBwwaRhkCTg0xPtqG3rDnkbvWTsgMDr3ulfesvhUmxGY6NiCsbAjtKS6LzP/XUEAEPzrTjbH+ysGyBrJxoUp/seOwCYMqIAFlPg7ZtzYjPbrKDAzuuNbG4y6RcDO0o7TrcXO4+3AQgsSREtPn0YgIGXY5XiiTSvigV849SmjSxUfs5lxi6j+Vo+GWHrcSt7jylzMbCjtLP7ZDt6XF4U55hx2tC8ZB8OpQjR9uTDwy1o7nT0uz7QoDj9M3ZA720ILJ7IbCajQSmo4XIsMbCjtLPjaGB/nSSxKSv5VJbk4PThBfDKwDt7G/pdrxRP6GCPHRDoZweweIK4z44CGNhR2gk0JuYyLPUmsnZvh1iOFRm7fB1UxQJQMjQAiycImFXlm5v8weGWQecmk74xsKO0IsuykrFj4QT1JdqebPmiCZ0Od6/r9NKgWCjOtWBcmW8rAosn6IyRhbCYDGjqdOBQkz3Zh0NJxMCO0sqR5i40dTphMRkwZUTh4HegjDK+LA9VQ3Lh9HixYX9gOdbt8cLu9ADQz1IsAFx+1kjkWU04axS/5GS6LLMR0yuLAHA5NtMxsKO0st2/DDttZCGsJi4/UW+SJAWWYz8LLMeKVieAfpZiAeC/FpyGT+5fjKkj+SWHuM+OfBjYUVphY2IazBJ/25P1+xrQ4/Jl6cT+uhyLEWajvt72jAYWEJGP2GfHwC6z6esdjnTvIzYmpkFMG1mE8oIs2J0evH/QN0s1ME5MP8uwRH2dNboIJoOE2rZunGjtSvbhUJIwsKO00dzpwKFG36bg4IpAomAGgxRoVvzZKQBAR49+mhMTDSTHYlL2HjNrl7kY2FHaENWw48vyUJRjSfLRUCq70L/PrnrvKbg9Xt01JyYaCPfZEQM7ShsisJs5hvvrKLxZVSUoyjGjxe7E9qOtumtOTDSQWQzsMh4DO0obSmNiLsPSIExGA86f6FuOffuz+qCMHZdiSd9mji6BJAGHmuxo6OhJ9uFQEjCwo7TQ4/JgV207AOBsZuwoAqJZcfWeU2jvFlMnmLEjfSvMMWNieQEAZu0yFQM7SgufnmiHyyNjaL4VlSXZyT4cSgPnjh+CHIsRtW3deP9gMwAWT1Bm4D67zMbAjtKCWIY9e0wxJIl9u2hwWWYjFkwYCiCwP5PFE5QJGNhltpQI7J566ilUVVUhKysLM2bMwObNm8PefuPGjZgxYwaysrIwduxYrFy5MkFHSsmiFE6wMTFFQUyhEFg8QZngbH9gt6++A21dziQfDSVa0gO71157DXfccQfuvfde1NTU4Nxzz8VFF12EY8eOhbz94cOHsXTpUpx77rmoqanBPffcg9tvvx1r1qxJ8JFToni9sjJKbCYbE1MUFk4sg9kYyPAyY0eZYEieFacNzQUAfOSf1kOZI+mB3WOPPYYbbrgBN954IyZNmoQVK1agsrISTz/9dMjbr1y5EqNGjcKKFSswadIk3Hjjjbj++uvx6KOPJvjIKVEONHbC1uNGjsWIyRUFyT4cSiMFWWbMPW1I4GfusaMMERgv1pzkI6FES+q7nNPpxI4dO/CTn/yk1+WLFy/G+++/H/I+W7duxeLFi3tdtmTJEqxatQoulwtmc+TfyN/d14D8Qlf0B04R8Xjc+KRZgmnPKRiNsb/Uth3yvTGdWVkEk87mfJL2LpxSjo2fNwJgxo4yxzlVJXjlw2PYsL8RM0bXD3g7td6nCRhXlotxZfnJPozkBnZNTU3weDwYNmxYr8uHDRuG+vrQL8T6+vqQt3e73WhqakJFRUW/+zgcDjgcDuVnm80GAPjxXz+DwZoT79OgsIx47vNPVPlN0ysL4XIxEA9FnBeen/4WjPf19ZJlINcixXWOeJ4Th+c6PmdV+lY3vmjoxM0v7Rjk1uq9T2ey2xeehh+ed9qA1yfqtZwS4XnfKkdZlsNWPoa6fajLhYcffhgPPPBAv8tH58kwZcnRHi4lQbZRRpntc6xd+3myDyWlVVdXJ/sQUtI3Rktod0n4bOsG7FahqJrnOXF4rmO3ZIQBn9vYRSBRGo99jrVr9w94fVdXV0KOI6mB3ZAhQ2A0Gvtl5xoaGvpl5YTy8vKQtzeZTCgtLQ15n7vvvhvLli1TfrbZbKisrMRfbp034H0ofi6XC9XV1Vi0aFFUS+QUPZ7r8Jaq9Ht4nhOH5zp+kbzueZ4Tp7k5MfsdkxrYWSwWzJgxA9XV1fj617+uXF5dXY1LL7005H3mzJmDf/zjH70uW7duHWbOnDngi9JqtcJqtfa73Gw284WcADzPicNznRg8z4nDc50YPM/aS9T5TfpO9GXLluEPf/gDnnvuOezduxd33nknjh07hptvvhmAL9t23XXXKbe/+eabcfToUSxbtgx79+7Fc889h1WrVmH58uXJegpEREREKSHpe+yuvPJKNDc34+c//znq6uowZcoUrF27FqNHjwYA1NXV9eppV1VVhbVr1+LOO+/E7373OwwfPhxPPPEELr/88mQ9BSIiIqKUkPTADgBuueUW3HLLLSGve/755/tdNn/+fHz88ccaHxURERFRekn6UiwRERERqYOBHREREZFOMLAjIiIi0gkGdkREREQ6wcCOiIiISCcY2BERERHpBAM7IiIiIp1gYEdERESkEwzsiIiIiHSCgR0RERGRTqTESLFEk2UZANDR0QGz2Zzko9Evl8uFrq4u2Gw2nmeN8VwnBs9z4vBcJwbPc+J0dHQACMQgWsnIwK65uRkAUFVVleQjISIiokzS3NyMwsJCzX5/RgZ2JSUlAIBjx45penIznc1mQ2VlJY4fP46CgoJkH46u8VwnBs9z4vBcJwbPc+K0t7dj1KhRSgyilYwM7AwG39bCwsJCvpAToKCggOc5QXiuE4PnOXF4rhOD5zlxRAyi2e/X9LcTERERUcIwsCMiIiLSiYwM7KxWK+6//35YrdZkH4qu8TwnDs91YvA8Jw7PdWLwPCdOos61JGtdd0tERERECZGRGTsiIiIiPWJgR0RERKQTDOyIiIiIdIKBHREREZFO6DKwa21txbXXXovCwkIUFhbi2muvRVtbW9j7yLKMn/3sZxg+fDiys7OxYMEC7N69u9dtbrrpJpx22mnIzs7G0KFDcemll2Lfvn0aPpPUpsV5bmlpwQ9/+ENMmDABOTk5GDVqFG6//Xa0t7dr/GxSm1av6WeeeQYLFixAQUEBJEka9Hfq0VNPPYWqqipkZWVhxowZ2Lx5c9jbb9y4ETNmzEBWVhbGjh2LlStX9rvNmjVrMHnyZFitVkyePBlvvPGGVoefNtQ+z7t378bll1+OMWPGQJIkrFixQsOjTx9qn+dnn30W5557LoqLi1FcXIwLLrgAH374oZZPIW2ofa5ff/11zJw5E0VFRcjNzcWZZ56JF198MfoDk3XowgsvlKdMmSK///778vvvvy9PmTJF/upXvxr2Po888oicn58vr1mzRt61a5d85ZVXyhUVFbLNZlNu8/vf/17euHGjfPjwYXnHjh3yJZdcIldWVsput1vrp5SStDjPu3btkr/xjW/Ib775pnzgwAH5nXfekcePHy9ffvnliXhKKUur1/Tjjz8uP/zww/LDDz8sA5BbW1s1fiap5dVXX5XNZrP87LPPynv27JF/9KMfybm5ufLRo0dD3v7QoUNyTk6O/KMf/Ujes2eP/Oyzz8pms1n+61//qtzm/fffl41Go/zQQw/Je/fulR966CHZZDLJ27ZtS9TTSjlanOcPP/xQXr58ufzKK6/I5eXl8uOPP56gZ5O6tDjPV199tfy73/1Orqmpkffu3St/73vfkwsLC+UTJ04k6mmlJC3O9fr16+XXX39d3rNnj3zgwAF5xYoVstFolN9+++2ojk13gd2ePXtkAL3eRLdu3SoDkPft2xfyPl6vVy4vL5cfeeQR5bKenh65sLBQXrly5YCP9cknn8gA5AMHDqj3BNJEIs/zn//8Z9liscgul0u9J5BGEnGu169fn5GB3axZs+Sbb76512UTJ06Uf/KTn4S8/V133SVPnDix12U33XSTPHv2bOXnK664Qr7wwgt73WbJkiXyt7/9bZWOOv1ocZ6DjR49moGdrP15lmVZdrvdcn5+vvzCCy/Ef8BpLBHnWpZlefr06fJPf/rTqI5Nd0uxW7duRWFhIc455xzlstmzZ6OwsBDvv/9+yPscPnwY9fX1WLx4sXKZ1WrF/PnzB7yP3W7H6tWrUVVVhcrKSnWfRBpI1HkGfIOTCwoKYDJl5GjjhJ7rTOJ0OrFjx45e5wgAFi9ePOA52rp1a7/bL1myBNu3b4fL5Qp7m0w971qdZ+otUee5q6sLLpdL80H2qSwR51qWZbzzzjvYv38/5s2bF9Xx6S6wq6+vR1lZWb/Ly8rKUF9fP+B9AGDYsGG9Lh82bFi/+zz11FPIy8tDXl4e3n77bVRXV8Nisah09OlD6/MsNDc348EHH8RNN90U5xGnr0Sd60zT1NQEj8cT1Tmqr68PeXu3242mpqawt8nU867VeabeEnWef/KTn2DEiBG44IIL1DnwNKTluW5vb0deXh4sFgsuvvhiPPnkk1i0aFFUx5c2gd3PfvYzSJIU9r/t27cDACRJ6nd/WZZDXh6s7/Wh7nPNNdegpqYGGzduxPjx43HFFVegp6cnzmeXOlLlPAOAzWbDxRdfjMmTJ+P++++P41mlplQ615ks2nMU6vZ9L+d570+L80z9aXmef/WrX+GVV17B66+/jqysLBWONr1pca7z8/Oxc+dOfPTRR/jFL36BZcuWYcOGDVEdV9qsbd1222349re/HfY2Y8aMwaeffopTp071u66xsbFftCyUl5cD8EXUFRUVyuUNDQ397iOqEsePH4/Zs2ejuLgYb7zxBq666qpon1JKSpXz3NHRgQsvvBB5eXl44403YDabo30qKS9VznWmGjJkCIxGY79v2OHOUXl5ecjbm0wmlJaWhr1Npp53rc4z9ab1eX700Ufx0EMP4T//+Q/OOOMMdQ8+zWh5rg0GA8aNGwcAOPPMM7F37148/PDDWLBgQcTHlzYZuyFDhmDixIlh/8vKysKcOXPQ3t7eqxz7gw8+QHt7O+bOnRvyd1dVVaG8vBzV1dXKZU6nExs3bhzwPoIsy3A4HOo8yRSQCufZZrNh8eLFsFgsePPNN3X7zTAVznUms1gsmDFjRq9zBADV1dUDnqM5c+b0u/26deswc+ZM5cvHQLfJ1POu1Xmm3rQ8z7/+9a/x4IMP4u2338bMmTPVP/g0k8jXdEwxRlSlFmniwgsvlM844wx569at8tatW+WpU6f2aw0xYcIE+fXXX1d+fuSRR+TCwkL59ddfl3ft2iVfddVVvVpDHDx4UH7ooYfk7du3y0ePHpXff/99+dJLL5VLSkrkU6dOJfT5pQotzrPNZpPPOecceerUqfKBAwfkuro65b9MbSsjy9qca1mW5bq6OrmmpkZ+9tlnZQDypk2b5JqaGrm5uTlhzy2ZRMuCVatWyXv27JHvuOMOOTc3Vz5y5Igsy7L8k5/8RL722muV24uWBXfeeae8Z88eedWqVf1aFrz33nuy0WiUH3nkEXnv3r3yI488wnYnGpxnh8Mh19TUyDU1NXJFRYW8fPlyuaamRv7iiy8S/vxShRbn+Ze//KVssVjkv/71r73ejzs6OhL+/FKJFuf6oYcektetWycfPHhQ3rt3r/y///u/sslkkp999tmojk2XgV1zc7N8zTXXyPn5+XJ+fr58zTXX9GvjAEBevXq18rPX65Xvv/9+uby8XLZarfK8efPkXbt2KdfX1tbKF110kVxWViabzWZ55MiR8tVXXz1gu4lMoMV5Fm03Qv13+PDhxDyxFKTFuZZlWb7//vtDnuvg36N3v/vd7+TRo0fLFotFPuuss+SNGzcq1333u9+V58+f3+v2GzZskKdPny5bLBZ5zJgx8tNPP93vd/7lL3+RJ0yYIJvNZnnixInymjVrtH4aKU/t83z48OGQr92+vyfTqH2eR48eHfI833///Ql4NqlN7XN97733yuPGjZOzsrLk4uJiec6cOfKrr74a9XFJsuzfvUdEREREaS1t9tgRERERUXgM7IiIiIh0goEdERERkU4wsCMiIiLSCQZ2RERERDrBwI6ISEVerxc/+MEPUFFRgR/84Afwer3JPiQiyiAM7IiIVPTvf/8bn3/+OdauXYu9e/fi3//+d7IPiYgyCAM7Isp4TqcT48aNw3vvvRf37yosLERxcTHGjx+P0tJSlJSUDHofh8OBUaNGYceOHXE/PhFlNgZ2RJTxnnnmGYwePRpf/vKXo77vz372M3z7299Wfp47dy6cTicKCwvh8XhwzjnnKNd1dnbioosuwrx58zBhwgQ8//zzAACr1Yrly5fj//7f/xv3cyGizMbAjogy3pNPPokbb7wxpvu++eabuPTSS5WfXS4Xtm/fjrvuugsffvgh3G63cl1OTg7efPNNbNq0CatWrcLTTz+tXHfNNddg8+bN2Lt3b+xPhIgyHgM7ItKN5uZmfPvb30ZxcTEkSer1n8iO9fXxxx/jwIEDuPjii3tdXltbiyuvvBLFxcUoLS3FpZdeiiNHjvS6zfHjx/HZZ5/hoosuUi576623YDab8fOf/xwmkwlvvfWWcp3BYIDZbEZjYyPuv/9+rFixQrmutLQUc+fOxSuvvBL3eSCizMXAjoh044477sD777+P1157DXv27FGycE8++STmzZsX8j6bNm3Cl770JRQUFCiXdXV1YeHChcjLy8OmTZuwZcsW5OXl4cILL4TT6VRu9+abb2LevHkoKipSLlu9ejWuuuoqmM1mXHXVVVi9enWvx9u2bRuuuuoqrFixAnPmzOl13axZs7B58+Z4TwMRZTAGdkSkCzabDS+//DJ+/etfY/HixZg0aRKefvppjBgxAi6XC2PHjg15vyNHjmD48OG9Lnv11VdhMBjwhz/8AVOnTsWkSZOwevVqHDt2DBs2bFBu9/e//73XMmxDQwPWrl2L73znOwCA73znO3jrrbfQ0NCgXH/uueeivr4e3//+93HhhRf2etwRI0b0ywoSEUXDlOwDICJSw8GDByHLMubOnatcZjKZMGvWLHz66acD3q+7uxtZWVm9LtuxYwcOHDiA/Pz8Xpf39PTg4MGDAHyB5MaNG/Hss88q17/44ouYOHEipk2bBgA488wzMXHiRLz00ktYtmwZysrK4HK5BjyW7OxsdHV1Rf6kiYj6YGBHRLpgNpsBAB6Pp9flHo8HRqNxwPsNGTIEu3bt6nWZ1+vFjBkz8PLLL/e7/dChQwEA//rXvzBp0iSMHj1auW716tXYs2cPTCZTr9+1evVqLFu2bNDn0NLSovx+IqJYMLAjIl047bTTkJWVhffeew9jxowBEKhQDRdUTZ8+HU8//TRkWYYkSQCAs846C6+99hrKysp67b0L9ve//x1f+9rXlJ8/+ugj7NmzBxs2bOjVu66trQ3z5s3D9u3bMXPmzLDP4bPPPsP06dMjfcpERP1wjx0R6UJ2djZuu+023HXXXXj77bexZ88efP/730dPTw9uuOGGAe+3cOFC2O127N69W7nsmmuuwZAhQ3DppZdi8+bNOHz4MDZu3Igf/ehHOHHiBNxuN/71r3/12l+3evVqzJo1C/PmzcOUKVOU/77yla9gzpw5/YooQtm8eTMWL14c34kgoozGwI6IdOMXv/gFrrjiClx33XU466yzcODAAfz73//uVbXaV2lpKb7xjW/0WnbNycnBpk2bMGrUKHzjG9/ApEmTcP3116O7uxsFBQXYuHEj8vLyMGPGDAC+vXevvPIKLr/88pCPcfnll+OVV15BT0/PgMexdetWtLe345vf/GZsT56ICIAky7Kc7IMgIkqmXbt24YILLghZMBHK7bffDrfbjaeeekq1Y/jWt76F6dOn45577lHtdxJR5uEeOyLKeFOnTsWvfvUrHDlyBFOnTh309lOmTOnXgy4eDocD06ZNw5133qna7ySizMSMHREREZFOcI8dERERkU4wsCMiIiLSCQZ2RERERDrBwI6IiIhIJxjYEREREekEAzsiIiIinWBgR0RERKQTDOyIiIiIdIKBHREREZFOMLAjIiIi0gkGdkREREQ68f8BpHpd+YJ4qo4AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma_moments: [4.28977468e+01 7.90903554e-05 6.05283055e+01 3.89164821e+00\n",
      " 1.20345287e+02 1.77937287e+01 2.66092559e+02]\n",
      "hb_donor_moment: [5.11873012 2.60564439 0.85940158]\n",
      "hb_acceptor_moment: [5.76486365 3.31610676 1.37737025]\n",
      "energy_dielectric: -32.83919025918012\n",
      "dipole_moment: None\n",
      "area_total: 42.897746777517476\n",
      "volume: 25.36550524969701\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAXvRJREFUeJzt3Xd4VFX+BvD3Tk3vlSQkoYMJXZrSVLqIoiIC6lpWrIjIqqzuD7AA6q6yiwLCKhYEUbFjIatUA1KkJ5RAQk0nfVKm3N8fk5kQU5hJZubembyf5+FZMrmZfOeYHd6cc8/3CKIoiiAiIiIit6eQugAiIiIicgwGOyIiIiIPwWBHRERE5CEY7IiIiIg8BIMdERERkYdgsCMiIiLyEAx2RERERB6CwY6IiIjIQzDYEREREXkIBjsicisffPABBEFo8s/WrVsBAIsWLcLXX3/d5Nfv27fPtYXbKSEhAX/5y1+kLoOI3IxK6gKIiFpizZo16NatW4PHe/ToAcAc7O644w7ceuutLq6MiEg6DHZE5JaSkpLQv39/qcsgIpIVLsUSkccRBAEVFRX48MMPrUu0I0aMqHdNWVkZHn30UYSFhSE0NBSTJ0/GpUuX6l2zYcMGjB49GtHR0fD29kb37t3x/PPPo6Kiot51f/nLX+Dn54eMjAyMHz8efn5+iIuLwzPPPIPq6up6116+fBmPPfYYYmJioNFo0KFDB7zwwgsNriMiagkGOyJyS0ajEQaDod4fo9EIANi1axe8vb0xfvx47Nq1C7t27cLy5cvrff1DDz0EtVqNdevW4fXXX8fWrVsxY8aMetecOnUK48ePx3vvvYeffvoJs2fPxmeffYaJEyc2qEev1+OWW27BjTfeiG+++QYPPPAA3nrrLbz22mvWa6qqqjBy5Eh89NFHmDNnDjZt2oQZM2bg9ddfx+TJk50wSkTU1nAplojc0qBBgxo8plQqYTAYMGjQICgUCoSHhzd6HQCMHTsW//nPf6wfX758Gc8++yxycnIQFRUFAHjxxRetnxdFEddddx26d++O4cOH4/Dhw+jZs6f18zU1NVi4cCHuvPNOAMCNN96Iffv2Yd26dfi///s/AMCHH36Iw4cP47PPPrNeN2rUKPj5+eG5555DSkoKRo0a1cqRIaK2jDN2ROSWPvroI+zdu7fen99//93mr7/lllvqfWwJaWfPnrU+dubMGUybNg1RUVFQKpVQq9UYPnw4ACA9Pb3e1wuC0GAmr2fPnvWe79dff4Wvry/uuOOOetdZdr/+8ssvNtdPRNQYztgRkVvq3r17qzZPhIaG1vtYq9UCACorKwEA5eXlGDp0KLy8vPDKK6+gS5cu8PHxwfnz5zF58mTrdRY+Pj7w8vJq8JxVVVXWjwsLCxEVFQVBEOpdFxERAZVKhcLCwha/HiIigMGOiKhRv/76Ky5duoStW7daZ+kAoLi4uMXPGRoait9//x2iKNYLd3l5eTAYDAgLC2tNyUREXIolIs+k1WobzKrZwxK8LDN5Fu+++26Ln/PGG29EeXl5g8bJH330kfXzREStwRk7InJLR48ehcFgaPB4x44dER4ejuTkZGzduhXfffcdoqOj4e/vj65du9r8/EOGDEFwcDAeeeQRzJ8/H2q1Gp988gkOHTrU4prvvfdevPPOO7jvvvuQlZWF5ORk7Ny5E4sWLcL48eNx0003tfi5iYgABjsiclP3339/o4+vXr0aDz30EP7973/j8ccfx9SpU6HT6TB8+HDrcWO2CA0NxaZNm/DMM89gxowZ8PX1xaRJk7Bhwwb07du3RTV7eXlhy5YteOGFF/DGG28gPz8fMTExmDt3LubPn9+i5yQiupIgiqIodRFERERE1Hq8x46IiIjIQzDYEREREXkIBjsiIiIiD8FgR0REROQhGOyIiIiIPASDHREREZGH8Pg+diaTCZcuXYK/v3+D8xmJiIiIXEEURZSVlaFdu3ZQKJw3r+bxwe7SpUuIi4uTugwiIiIinD9/HrGxsU57fo8Pdv7+/gCAzMxMhISESFyN59Pr9di8eTNGjx4NtVotdTkejWPtOhxr1+FYuw7H2rUuX76MxMREay5xFo8PdpblV39/fwQEBEhcjefT6/Xw8fFBQEAA3yicjGPtOhxr1+FYuw7H2rX0ej0AOP22MG6eICIiIvIQDHZEREREHoLBjoiIiMhDMNgREREReQgGOyIiIiIPwWBHRERE5CEY7IiIiIg8BIMdERERkYdgsCMiIiLyEAx2RERERB6CwY6IiIjIQzDYEREREXkIBjsiIiIiD8FgR0SyU1RRI3UJRERuicGOiGTl6wMX0eflFKzdfVbqUoiI3A6DHRHJytGLJQCA3zIKJK6EiMj9MNgRkazo9EYAwJn8CokrISJyPwx2RCQrlTXmYJdZWAGjSZS4GiIi98JgR0SyYgl2NQYTLhZVSlwNEZF7YbAjIlmxLMUCwOmCcgkrISJyPwx2RCQrlTUG699P5zHYERHZg8GOiGRFV1M3Y3emgBsoiIjswWBHRLJSecVS7Jl8ztgREdmDwY6IZKXyihm702x5QkRkFwY7IpKVK5di88uqUVqll7AaIiL3wmBHRLJimbFTKgQAbFRMRGQPBjsikg2D0YQaowkA0DnCDwDvsyMisgeDHRHJxpUbJ5JiAgFwxo6IyB4MdkQkG5ZlWIUAdIvyBwCc5owdEZHNGOyISDYsGye81Up0tC7FcsaOiMhWDHZEJBvWYKdRoWOYOdhlFlbAaBKlLIuIyG0w2BGRbFjusfPRKBET7A2NSoEagwkXiyolroyIyD0w2BGRbFjusfPRKKFUCEgM9QUAnC7gfXZERLZgsCMi2dDVGAAAXmolAKBDeG2wy2OwIyKyBYMdEcnGlUuxANAxvHYDRQE3UBAR2YLBjohk48qlWKBuxo5NiomIbMNgR0SyceWuWKBuxu40W54QEdmEwY6IZMO6FPune+zyy6pRWqWXrC4iInfBYEdEsmHZPOFduxTr76VGuL8WABsVExHZgsGOiGSjbilWaX2sI++zIyKyGYMdEclG1Z+WYgGgQziPFiMishWDHRHJRuMzdpYNFJyxIyK6GgY7IpKNxoJdXcsTztgREV0Ngx0Rycaf+9gBQMcw84xdZmEFjCZRkrqIiNwFgx0RyYal3Ym3WmV9LCbYGxqVAjUGEy4WVUpVGhGRW2CwIyLZ0DUyY6dUCEgMrT0ztoD32RERNYfBjohko/JPfewsOkbUBrs8BjsiouYw2BGRbFg3T6jrB7sOtffZnSngBgoiouYw2BGRbFiPFPvTjJ1lZyxn7IiImsdgR0SyUbcrVlXvcUsvO87YERE1j8GOiGShxmCCobadyZ/vsbPM2OWXVaO0Su/y2oiI3IXkwW779u2YOHEi2rVrB0EQ8PXXX9f7vCiKWLBgAdq1awdvb2+MGDECx44dk6ZYInIay2wd0PAeO38vNSL8tQDYqJiIqDmSB7uKigr06tULb7/9dqOff/311/Hmm2/i7bffxt69exEVFYVRo0ahrKzMxZUSkTPp9OYdsSqFAI2q4VtT3QkUvM+OiKgpqqtf4lzjxo3DuHHjGv2cKIpYunQpXnjhBUyePBkA8OGHHyIyMhLr1q3DzJkzXVkqETlRZSPHiV2pQ7gfdp+5zDNjiYiaIXmwa05mZiZycnIwevRo62NarRbDhw9Hampqo8Guuroa1dXV1o9LS0sBAHq9Hno9781xNssYc6ydz9PGulRn/v+tj1rZ6GtKCPEGAGTklrn8NXvaWMsZx9p1ONau5apxlnWwy8nJAQBERkbWezwyMhJnz55t9GsWL16MhQsXNnh8y5Yt8PHxcXyR1KiUlBSpS2gzPGWsT5cCgApGfRV++OGHBp8vLBIAKHE4M7fRz7uCp4y1O+BYuw7H2jV0Op1Lvo+sg52FIAj1PhZFscFjFvPmzcOcOXOsH5eWliIuLg4jR45EaGioU+sk828kKSkpGDVqFNRqtdTleDRPG+sdpwqAY38gLCgA48cPbvD55CId3j2+E4V6JcaMHQ2lovH3AGfwtLGWM46163CsXauwsNAl30fWwS4qKgqAeeYuOjra+nheXl6DWTwLrVYLrVbb4HG1Ws0fXBfieLuOp4x1jckc1Hy1qkZfT3xYADQqBWoMJuSVG9A+1PUz8J4y1u6AY+06HGvXcNUYS74rtjmJiYmIioqqN01cU1ODbdu2YciQIRJWRkSOVlm7K/bPp05YKBUCEkNrT6DgBgoiokZJPmNXXl6OjIwM68eZmZk4ePAgQkJC0L59e8yePRuLFi1C586d0blzZyxatAg+Pj6YNm2ahFUTkaM1dU7slTpG+OJEbhlO55djZLcIV5VGROQ2JA92+/btw8iRI60fW+6Pu++++/DBBx/g2WefRWVlJR577DEUFRVh4MCB2Lx5M/z9/aUqmYic4GrtTgCgQxiPFiMiao7kwW7EiBEQRbHJzwuCgAULFmDBggWuK4qIXK7unNjmZ+wA4HQel2KJiBoj63vsiKjt0OktS7FN/77JGTsiouYx2BGRLNgyY2c5Viy/rBqlVWyqSkT0Zwx2RCQLuhrzrtjm7rHz91Ijwt/czuhMPmftiIj+jMGOiGTBll2xQN2s3Rm2PCEiaoDBjohkoUp/9aVYAOgYbr7Pjr3siIgaYrAjIlnQ2dDuBAA61AY7LsUSETXEYEdEsmDvUixn7IiIGmKwIyJZqNsV23x7zU61M3ZZhToYTU33wCQiaosY7IhIFir1ti3FtgvyhkalQI3BhItFla4ojYjIbTDYEZEs6GzoYwcASoWADmFcjiUiagyDHRHJQmVtH7urBTuA99kRETWFwY6IJCeK4hVHitkQ7MIsLU+4M5aI6EoMdkQkuWqDCWLtPoir3WMHAJ0izMEuI6/MmWUREbkdBjsikpxlRyxw9V2xANAt2h8AcDynDKLInbFERBYMdkQkOcsyrEalgFIhXPX6DmF+UCkElFUZcKmkytnlERG5DQY7IpKcZeOELffXAeYAaFmOPZ5d6rS6iIjcDYMdEUnO1lYnV+oaVbccS0REZgx2RCS5ShvPib1St6gAAAx2RERXYrAjIslZ7rGzZ8bOuoGCS7FERFYMdkQkOeuMnY332AFAt9ql2DMFFag2GK9yNRFR28BgR0SS01mXYq/e6sQiKsALgd5qGE0iMvJ4AgUREcBgR0QyUGlZirVjxk4QBOus3fFs3mdHRAQw2BGRDNhzTuyVLMHuRC6DHRERwGBHRDKga8GuWADoFm3eGZvODRRERAAY7IhIBlqyeQKo62V3gi1PiIgAMNgRkQy0pEExAHSNNAe7vLJqFJZXO7wuIiJ3w2BHRJKzbJ6wZ1csAPhqVYgP9QHAWTsiIoDBjohkoLKFM3ZA3awdT6AgImKwIyIZ0NXuirX3HjugbgPF8RxuoCAiYrAjIsm1dFcsUNfyhDN2REQMdkQkA1UtOCvWwhLsTuaWwWgSHVoXEZG7YbAjIsm1ZsYuPtQXXmoFqvQmnC2scHRpRERuhcGOiCRX1+7Evl2xAKBUCOgSyX52REQAgx0RyYC13UkLNk8Adcux6Qx2RNTGMdgRkeRa0+4EALpG1e6M5dFiRNTGMdgRkaRMJvGKBsUtC3bdLUeL5XLGjojaNgY7IpJUlcFo/XvLZ+zMwe5soQ4V1QaH1EVE5I4Y7IhIUpaNEwDgpWpZsAv10yLcXwvA3PaEiKitYrAjIklZ7q/zUiugUAgtfh42KiYiYrAjIolV6lve6uRK1mDHDRRE1IYx2BGRpKzNiVvY6sSim2VnLGfsiKgNY7AjIknpasybHVq6I9aiW3TdUqwo8mgxImqbGOyISFKt7WFn0SnCD0qFgJJKPXJKqxxRGhGR22GwIyJJtfbUCQutSokOYb4AuBxLRG0Xgx0RSUrnoBk7oK6f3fFsBjsiapsY7IhIUnVLsa3bFQsA3aPNGyhO5HBnLBG1TQx2RCQpnbWPXetn7NjLjojaOgY7IpJUXR87xy3FZuSVo8ZgavXzERG5GwY7IpJUZW27E0cEu5ggb/hrVTCYRJwpKG/18xERuRsGOyKSlLVBsQOCnSAIdf3suIGCiNogBjsiklSlg06esOjK++yIqA2TfbAzGAx48cUXkZiYCG9vb3To0AEvvfQSTCbeP0PkCRzZ7gS48mixq++M1RtNWPq/k/j+8CWHfG8iIqm1vr+Ak7322mtYuXIlPvzwQ1xzzTXYt28f7r//fgQGBuKpp56SujwiaiVrg2IHtDsBrtgZa8NS7EvfpeHj3WehVSkwomsE/LSyf0skImqW7N/Fdu3ahUmTJmHChAkAgISEBKxfvx779u2TuDIicgRHHSlm0aU22OWUVqFYV4MgH02j1320Kwsf7z4LAKg2mPC/tFzc2ifGITUQEUlF9sHu+uuvx8qVK3Hy5El06dIFhw4dws6dO7F06dJGr6+urkZ1dbX149JS83KMXq+HXq93RcltmmWMOdbO5yljXVFjrl+jEB3yWryVQGyQFy4UV+HohSIMTAxpcM2OjAIs/C4NANAx3Ben8yvwzcELmJAU0ehzespYuwOOtetwrF3LVeMsiKIouuQ7tZAoivj73/+O1157DUqlEkajEa+++irmzZvX6PULFizAwoULGzy+bt06+Pj4OLtcIrLTqweUyKsS8GQPAzoFOuY5Vx9X4GiRArcnGDEsuv5bXG4l8NYRJSqNAgaEm3BjOxMWH1JBKYh4uZ8RvmrH1EBEdCWdTodp06ahpKQEAQEBTvs+sp+x27BhA9auXYt169bhmmuuwcGDBzF79my0a9cO9913X4Pr582bhzlz5lg/Li0tRVxcHEaOHInQ0FBXlt4m6fV6pKSkYNSoUVCr+S+kM3nKWC9J2w5UVWHksOuQHOOYZHdccwpHt2VCFdYe48dfY328WKfHHe/+jkqjDv3aB+H9+/tDq1Lgq5xUHM8tB2J7Yny/2AbP5ylj7Q441q7DsXatwsJCl3wf2Qe7v/3tb3j++ecxdepUAEBycjLOnj2LxYsXNxrstFottFptg8fVajV/cF2I4+067j7Wls0TAT5ah72OHu2CAAAnciusz6k3mjBrw36cvaxDTJA33r23P/y8ze8VE3vH4PjPJ/DjsTxMG5TY5PO6+1i7E46163CsXcNVYyz7dic6nQ4KRf0ylUol250QeYi6BsWO+z2ze22T4pO5ZTCZRIiiiP/75hh2nSmEr0aJ9/7SH2F+db8A3twzGgDwW0YBCsqrG31OIiJ3IPtgN3HiRLz66qvYtGkTsrKy8NVXX+HNN9/EbbfdJnVpRNRKRpNoPdPVUQ2KASAh1BcalQK6GiPOF+nwQWoW1u85B0EA/nN3H2uvO4v4UF/0ig2ESQR+PJLtsDqIiFxN9sFu2bJluOOOO/DYY4+he/fumDt3LmbOnImXX35Z6tKIqJUsy7CA49qdAIBKqUDnCD8AwKrtZ/Dy9+YdsPPGdcON3SMb/ZqJvdoBAL47xGBHRO5L9sHO398fS5cuxdmzZ1FZWYnTp0/jlVdegUbTeG8qInIfuhoDAEAQAK3KsW9Hllm5T34/B5MI3NkvFn8d2qHJ68cnm5dj9569jOySSofWQkTkKrIPdkTkuazNidVKCILg0Oe2nEABAAMSQvDKbUnNfo92Qd64NiEYoghsOsxZOyJyTwx2RCSZuo0TjluGteifEAwAiAvxxooZfaFVXf17WJdjGeyIyE0x2BGRZJwZ7Pq0D8bGR4fguyeuR6hfwxZIjRmXFA2FABw6X4xzhTqH10RE5GwMdkQkmSq9ZSnWOS01+8UHN3lWbGPC/bUY3NHcyPz7I5ecUhMRkTMx2BGRZJw5Y9dSE3tydywRuS8GOyKSjGVXrCNbnbTW2KQoqBQC0rNLkZFXJnU5RER2YbAjIslYdsU6sjlxawX5aDC0cxgAztoRkfthsCMiyVgaFMtpKRao2x37/eFLEEVR4mqIiGzHYEdEkrHcYyenpVgAGNUjEhqVAqfzK5CezeVYInIfDHZEJBlrg2KNc3bFtpS/lxo3dI0AAHx3mLtjich9MNgRkWQsM3ZeMrrHzuLmXuYjxr47xOVYInIfDHZEJBnLPXZyW4oFgBu6RcBHo8SFokocvlgqdTlERDZhsCMiyVTKsN2JhY9GhZu6RwIANh3JkbgaIiLbMNgRkWTk2KD4SpbdsT8czYGJq7FE5AYY7IhIMtZ2JzK8xw4AhnUJg7+XCrml1cjk5lgicgMMdkQkGbm2O7HQqpQYc00UAOCPAr5dEpH88Z2KiCRjPXlCZu1OrmRZjj1YKMDE9VgikjkGOyKSjJx3xVoM6RgKb7UC5QYBmYU6qcshImoWgx0RSUZXuytWrvfYAYBaqUD36AAAwNGLJRJXQ0TUPAY7IpKM3HfFWiS1qw12l9jPjojkjcGOiCRT5QZLsQCQHMNgR0TugcGOiCShN5qgN5o3I/io5bt5AgCuqZ2xS8sug5EbKIhIxhjsiEgSlmVYQP5LsR3CfKFRiNDVGHEmv1zqcoiImsRgR0SSsLQ6USoEqJWCxNU0T6kQEOtr/vvhC9xAQUTyxWBHRJKwtjpRKyEI8g52ABDnZ16CPcKdsUQkYwx2RCQJa6sTmS/DWrT3ZbAjIvljsCMiSVTK/DixP7PM2KVdKoXBaJK4GiKixjHYEZEkdG5wnNiVwr0AX40SlXojTudXSF0OEVGjGOyISBLWYKd2j7chhQD0qG17cvhCsbTFEBE1wT3eUYnI49Q1J3aPGTsASG7Ho8WISN4Y7IhIEu5ynNiVkmpPoDjMYEdEMsVgR0SSsOyKdZfNE0DdmbHcQEFEcsVgR0SSqLTeY+c+wS4+xAf+WhWqDSacyuMJFEQkPwx2RCQJS4Nid1qKVSgEXFO7HHuEJ1AQkQwx2BGRJHRu1sfOomdsEAA2KiYieWKwIyJJ1DUodp9dsQCQHBMIgBsoiEieWvSOqtfrkZOTA51Oh/DwcISEhDi6LiLycLrapVgvN7rHDqgLdunZpdAbTVAr+fsxEcmHze9I5eXlePfddzFixAgEBgYiISEBPXr0QHh4OOLj4/HXv/4Ve/fudWatRORB3O1IMYv4UB/4e6lQYzDhZG6Z1OUQEdVjU7B76623kJCQgNWrV+OGG27Al19+iYMHD+LEiRPYtWsX5s+fD4PBgFGjRmHs2LE4deqUs+smIjdXqXe/dicAIAiCddaOGyiISG5sWopNTU3Fli1bkJyc3OjnBwwYgAceeAArV67Ee++9h23btqFz584OLZSIPIvODdudWCTHBiL1dCGOXCzBVKmLISK6gk3B7vPPP7fpybRaLR577LFWFUREbYO7bp4AgJ4xQQC4M5aI5Meuu34rKiqsf8/KynJ0LUTUhtQdKeZ+mw8sS7HHs8tQY+AJFEQkHza/o86aNQuxsbFYsWIFAGDatGlOK4qIPJ+1QbHa/Wbs4kK8EeitRo2RGyiISF5sDnabN29Gbm4uDhw4gI0bNzqzJiJqA9x1VyxQfwPFYW6gICIZsTnYxcbGQqPRYMWKFfjwww9x6dIlZ9ZFRB5MFEXoatxzV6xFcmztzljeZ0dEMmJzsOvevTv0ej2USiXeffddBAUFObEsIvJk1QYTTKL5715uGux6WlqeXCyWthAioivYHOyWLVsGtVoNAIiOjsbBgwcbXMMGxURki6ra++sAwMcN250AQFJtsDuRU4Zqg/EqVxMRuUart6Pl5eXhzTffRFJSEgYNGuSImojIw1l2xGqUCqjc9Eiu2GBvBPuooTeKOJHDDRREJA8tekc1Go345ptvcOuttyIuLg6rV6/Grbfein379jm6PiLyQHWtTtxztg4wb6BI4gYKIpKZZvsM5Ofn480330RISAhmz56NEydOYM2aNVi7di0AYMqUKTCZTNi4cSN69OjhkoKJyP25847YK/WMDcSOUwU8WoyIZKPZGbtp06ZBp9MBAGJiYjBo0CBcunQJ77//Pi5duoRly5a5pEgi8iyWHbHueJzYlZJ5AgURyUyzwe748eOYPn06HnjgAVy+fBkPP/wwXnrpJUyYMAFKpevekC9evIgZM2YgNDQUPj4+6N27N/bv3++y709EjmVtTuzmM3aWlicnc8vqbQghIpJKs8HuxRdfxG233Ybhw4djyZIlyMrKQlJSEgYOHIi3334b+fn5Ti+wqKgI1113HdRqNX788UekpaXhX//6F9utELkxT1mKbRfohVBfDQwmEce5gYKIZKDZe+xmzpyJ6dOnQ6vVWlud5OfnY+3atVi9ejWefvppmEwmpKSkIC4uDv7+/g4v8LXXXkNcXBzWrFljfSwhIcHh34eIXKdu84T7HSd2JUEQkBwbiK0n8nHkQjF6xwU1eW21wYgF3x7DgXPFWP/XQQj21biuUCJqM666K9bPz88a6gAgPDwcTz/9NA4dOoTdu3fj0Ucfxcsvv4yIiAjccsstDi/w22+/Rf/+/XHnnXciIiICffr0werVqx3+fYjIdXTWc2Lds9XJlWw5WqykUo/73t+D9XvO43hOGfZkXXZVeUTUxrTq1+V+/fqhX79+ePPNN/H111/jgw8+cFBZdc6cOYMVK1Zgzpw5+Pvf/449e/Zg1qxZ0Gq1uPfeextcX11djerqauvHpaWlAAC9Xg+9Xu/w+qg+yxhzrJ3Pnce6oqoGAOClUrhF/c2NdY8oPwDAkQvFjX4+u6QKD330B07mlVsfu1RU4RavWwru/HPtbjjWruWqcRZEURRd8p1aSKPRoH///khNTbU+NmvWLOzduxe7du1qcP2CBQuwcOHCBo+vW7cOPj4+Tq2ViGzz03kBP15QYkikCXd1MEldTqsUVwPz/1BBARGvDTDiytsGs3XAynQlimsEBKhFRPmIOFmiwE0xJkxs796vm4jso9PpMG3aNJSUlCAgIMBp38emGbtHHnkEL7zwAuLi4q567YYNG2AwGDB9+vRWFweYjy/7c4+87t27Y+PGjY1eP2/ePMyZM8f6cWlpKeLi4jBy5EiEhoY6pCZqml6vR0pKCkaNGlVvCZ8cz53H+ujPJ4ELWejWMRHjx3WVupyram6sRVHEspPbUFBeg/jeQ9Cn9j673zMv4x/rDqK0xoAOYb54796++P5wNv71vwz4R8Ri/PgkCV6J/Lnzz7W74Vi7VmFhoUu+j03BLjw8HElJSRgyZAhuueUW9O/fH+3atYOXlxeKioqQlpaGnTt3Yv369YiNjcWqVascVuB1112HEydO1Hvs5MmTiI+Pb/R6rVYLrVbb4HG1Ws0fXBfieLuOO451tdG8UODr5V61NzXWPWOD8OvxPKTnVGBAh3B8f/gS5mw4hBqjCf3jg/Hf+/ojyEeDdsHmW0Pyy2rc6nVLwR1/rt0Vx9o1XDXGNgW7l19+GU8++STee+89rFy5EkePHq33eX9/f9x000147733MHr0aIcW+PTTT2PIkCFYtGgRpkyZgj179mDVqlUODY9E5FqVHnCk2JWSYgLx6/E8HL5Qgvd2ZuKVTWkQRWDMNZH499Q+8KptxBwZ4AUAyC2tkrJcIvJgNm+eiIiIwLx58zBv3jwUFxfj7NmzqKysRFhYGDp27AhBEJxS4LXXXouvvvoK8+bNw0svvYTExEQsXbrUYUu9ROR6ll2xPm5+8oRFz9qdsd8duoSNf1wAANw7OB7zJ14DpaLuvTEq0LyakMNgR0RO0qJdsUFBQS5tEHzzzTfj5ptvdtn3IyLnqmtQ7N597CwsJ1DUGM0bIp4d2xWPDm/4C69lxq6sygBdjcFjXj8RyYf7N5EiIrdjPSvWQ5ZiIwO80CXSDyqFgDen9MJjIzo1uorhp1VZT9vILa1u8Hkiotbir4tE5HLWe+w8ZCkWAD5/ZAh0NQZEB3o3eY0gCIgK8MKZggrklFQhMczXhRUSUVvAGTsicrlKvWecFXulQG91s6HOIiLAfJ9dXhnvsyMix2OwIyKX03nYrlh7RNXeZ5dTwmBHRI7HYEdELudpmyfsUdfyhPfYEZHjtehd9YsvvsBnn32Gc+fOoaampt7n/vjjD4cURkSeS+eB99jZir3siMiZ7J6x+89//oP7778fEREROHDgAAYMGIDQ0FCcOXMG48aNc0aNRORBRFG03mPXJpdiA2uXYhnsiMgJ7A52y5cvx6pVq/D2229Do9Hg2WefRUpKCmbNmoWSkhJn1EhEHqRKb7L+3ZM2T9gqsnbzBGfsiMgZ7A52586dw5AhQwAA3t7eKCsrAwDcc889WL9+vWOrIyKPY+lhB7Ttpdi80mqIoihxNUTkaewOdlFRUSgsLAQAxMfHY/fu3QCAzMxMvkkR0VVZ7q/TqhRQKJxzFKGcRfibg12N0YQinV7iaojI09gd7G644QZ89913AIAHH3wQTz/9NEaNGoW77roLt912m8MLJCLPUuWBPezsoVEpEOqrAcCWJ0TkeHbvil21ahVMJvM9Mo888ghCQkKwc+dOTJw4EY888ojDCyQiz6Jrw61OLCIDvFBYUYPcsir0QIDU5RCRB7H7nVWhUEChqJvomzJlCqZMmeLQoojIc7Xl5sQWkQFapGUDuZyxIyIHa1GD4h07dmDGjBkYPHgwLl68CAD4+OOPsXPnTocWR0Sep1Jv3jzRVpdiAbY8ISLnsTvYbdy4EWPGjIG3tzcOHDiA6mpz9/SysjIsWrTI4QUSkWexzNh5tcEdsRaWDRQ8fYKIHM3uYPfKK69g5cqVWL16NdRqtfXxIUOG8NQJIrqquuPE2m6ws8zYsZcdETma3cHuxIkTGDZsWIPHAwICUFxc7IiaiMiDVbbxXbEAmxQTkfPYHeyio6ORkZHR4PGdO3eiQ4cODimKiDxX3TmxbXtXLMBgR0SOZ3ewmzlzJp566in8/vvvEAQBly5dwieffIK5c+fisccec0aNRORB6nbFtmjvlkeIqg12BeU10BtNV7maiMh2dv/K/Oyzz6KkpAQjR45EVVUVhg0bBq1Wi7lz5+KJJ55wRo1E5EHqGhS33Rm7YB8N1EoBeqOIvLJqxAR5S10SEXkIu95ZjUYjdu7ciWeeeQYvvPAC0tLSYDKZ0KNHD/j5+TmrRiLyIJazYtviObEWCoWACH8vXCyuRE5JFYMdETmMXcFOqVRizJgxSE9PR0hICPr37++suojIQ+m4KxaAeQPFxeJK5PE+OyJyILtvcklOTsaZM2ecUQsRtQFsd2LGJsVE5Ax2B7tXX30Vc+fOxffff4/s7GyUlpbW+0NE1BxLu5O23KAYuHJnLJsUE5Hj2H338tixYwEAt9xyCwRBsD4uiiIEQYDRaHRcdUTkceqWYtvu5gmALU+IyDnsfmfdsmWLM+ogojaCS7FmlpYnOSUMdkTkOHYHu+HDhzf5uYMHD7amFiJqA6y7Ytt4sIuwnD5RxmBHRI7T6g6hJSUlWL58Ofr27Yt+/fo5oiYi8mCV1pMn2naws8zY5XLGjogcqMXB7tdff8WMGTMQHR2NZcuWYfz48di3b58jayMiD8SzYs0s99hV1BhRXm2QuBoi8hR2LcVeuHABH3zwAd5//31UVFRgypQp0Ov12LhxI3r06OGsGonIg9QdKda2g52vVgV/rQpl1QbklFShUwSbvBNR69k8Yzd+/Hj06NEDaWlpWLZsGS5duoRly5Y5szYi8jBGk4hqg/ls1La+KxYAImt72bFJMRE5is3vrJs3b8asWbPw6KOPonPnzs6siYg8lGUZFuA9doD59ImMvHI2KSYih7F5xm7Hjh0oKytD//79MXDgQLz99tvIz893Zm1E5GEsGycEAfBSt3rvltuz3GfHYEdEjmLzO+vgwYOxevVqZGdnY+bMmfj0008RExMDk8mElJQUlJWVObNOIvIAV+6IvbLBeVtlCXZ5PH2CiBzE7l+ZfXx88MADD2Dnzp04cuQInnnmGSxZsgQRERG45ZZbnFEjEXkInd68+7Ot74i1YJNiInK0Vq2FdO3aFa+//jouXLiA9evXO6omIvJQ3BFbn/VYMTYpJiIHcchNLkqlErfeeiu+/fZbRzwdEXmoKjYnrifScvoEZ+yIyEF49zIRuUzdjB1bnQBAlKXdSVk1TCZR4mqIyBMw2BGRy+gsp05wxg4AEOanhSAABpOIwooaqcshIg/AYEdELlNZw80TV1IrFQjzq12OZcsTInIABjsichnLUqwXg52V9T47BjsicgAGOyJymUouxTYQxSbFRORADHZE5DKWBsVciq1jbXnCJsVE5AAMdkTkMtwV25A12LHlCRE5AIMdEbmMjjN2DUSxSTERORCDHRG5jI67YhuIqN08wWPFiMgRGOyIyGWKdHoAQJCPRuJK5MPSpJi7YonIERjsiMhlimqb8Ib4qiWuRD4sS7FFOj2qDUaJqyEid8dgR0Quc7k22AVzxs4q0FsNjcr8VpzHnbFE1EoMdkTkMkU6y4wdg52FIAh1Gyi4HEtErcRgR0QuUaU3WnfFBjPY1WM5fYJNiomotdwu2C1evBiCIGD27NlSl0JEdrDM1qkUAvy17GN3JTYpJiJHcatgt3fvXqxatQo9e/aUuhQispP1/jpfDQRBkLgaeeFSLBE5itsEu/LyckyfPh2rV69GcHCw1OUQkZ2KKsytTkK4caIBy4wde9kRUWu5TbB7/PHHMWHCBNx0001Sl0JELXBZZ5mxY6uTP4tkLzsichC3uNHl008/xR9//IG9e/de9drq6mpUV9fdp1JaWgoA0Ov10Ov1TquRzCxjzLF2Pncb64LSSgBAkLfabWq2cPZYh/qYT+LIKalyu7FxNHf7uXZnHGvXctU4yz7YnT9/Hk899RQ2b94MLy+vq16/ePFiLFy4sMHjW7ZsgY+PjzNKpEakpKRIXUKb4S5jvee8AoACZQXZ+OGHi1KX0yLOGuuCKgBQIbu4Aps2/QDegug+P9eegGPtGjqdziXfRxBFUXTJd2qhr7/+GrfddhuUyrqzJY1GIwRBgEKhQHV1db3PNTZjFxcXh+zsbISGhrq09rZIr9cjJSUFo0aNglrNJTdncrexXvh9Otb+fh6Pj+iA2Td2krocuzh7rKv0RiS/9AsAYP/fRyLAW/7/PZ3F3X6u3RnH2rUKCwsRHR2NkpISBAQEOO37yH7G7sYbb8SRI0fqPXb//fejW7dueO655+qFOgDQarXQarUNnketVvMH14U43q7jLmNdXGkAAIT5e7lFvY1x1lir1WoEeqtRUqlHYaURoQFcXXCXn2tPwLF2DVeNseyDnb+/P5KSkuo95uvri9DQ0AaPE5F88dSJ5kUFeKGkUo/c0ip0ifSXuhwiclNusyuWiNzb5dp2JzwntnERltMn2PKEiFpB9jN2jdm6davUJRCRnYoqOGPXHEuT4rwynj5BRC3HGTsicjpRFOudPEENsUkxETkCgx0ROV1FjRE1RhMAnjzRFDYpJiJHYLAjIqezLMN6qRXw1iivcnXbxPNiicgRGOyIyOksy7CcrWtapGXzBIMdEbUCgx0ROV3dObEMdk2xzNjll1XDaJJ133gikjEGOyJyOu6IvbpQPy2UCgEmESgo585YImoZBjsicjrrjlguxTZJqRAQ7mdejuV9dkTUUgx2ROR0PHXCNpFsUkxErcRgR0ROx1MnbGPpZZfLJsVE1EIMdkTkdHX32PGg8eZEWXrZccaOiFqIwY6InI67Ym1jPX2C99gRUQsx2BGR03FXrG0i2aSYiFqJwY6InI6bJ2xj2TzBYEdELcVgR0ROZTKJKNKZN0/w5Inm1R0rxs0TRNQyDHZE5FRlVQbrSQpBDHbNiqzdPFFSqUeV3ihxNUTkjhjsiMipLBsn/LUqaFR8y2mOv1YFb7USAJdjiahl+C5LRE51ucK8rMgdsVcnCIK15QmbFBNRSzDYEZFTWZsTM9jZJMK/9vQJztgRUQsw2BGRU1lbnfiwObEtYoN9AAAXiiolroSI3BGDHRE5FZsT2ych1BzssgoqJK6EiNwRgx0ROVXdjB2DnS3iw3wBAGcLdRJXQkTuiMGOiJzqcgVn7OxhnbEr5IwdEdmPwY6InIqnTtgnPsQ8Y5dXVg1djUHiaojI3TDYEZFTWWfsuBRrk0AfNYJrN5pwOZaI7MVgR0ROZT1OjDN2NosPtdxnx+VYIrIPgx0ROZVlxi7El+1ObGW5zy6zgDN2RGQfBjsichqD0YSSytoGxVyKtRln7IiopRjsiMhpimtDnSAAQQx2NksI485YImoZBjsichpLD7sgbzWUCkHiatxH3Ywdl2KJyD4MdkTkNOxh1zIJtcEuu6QKVXqjxNUQkTthsCMip7H2sOMyrF2CfdQI8FIBAM5d5qwdEdmOwY6InKaQM3YtIggCEmqPFuOZsURkDwY7InIanhPbcrzPjohagsGOiJzmckVtqxPO2NmNZ8YSUUsw2BGR09SdE8vmxPbijB0RtQSDHRE5Dc+JbTnO2BFRSzDYEZHT1M3YMdjZyzJjd6m4EtUGtjwhItsw2BGR07CPXcuF+Wngq1HCJALnL1dKXQ4RuQkGOyJyGu6KbTlBEHhmLBHZjcGOiJyiSm9ERY15CZEzdi1Td2YsN1AQkW0Y7IjIKYp15lYnSoVgPUWB7MMZOyKyF4MdETnFlTtiBUGQuBr3VLczljN2RGQbBjsicgr2sGs9ztgRkb0Y7IjIKdjDrvUSa8+LvVBUCb3RJHE1ROQOGOyIyCnYw671Ivy18FIrYDSJuFjElidEdHUMdkTkFOxh13qCICChdjmWJ1AQkS0Y7IjIKSw97EIZ7FolvnYDBc+MJSJbMNgRkVMU8h47h+CMHRHZg8GOiJyC99g5hmVnbFYBgx0RXR2DHRE5xeUKc4Ni3mPXOglciiUiOzDYEZFT8JxYx4ivbXlyvkgHA1ueENFVyD7YLV68GNdeey38/f0RERGBW2+9FSdOnJC6LCJqhiiKuKyz7Iplg+LWiA7wgkalgN4oIrukSupyqIWMJhHvbMnA2t1nYTSJUpdDHkz2wW7btm14/PHHsXv3bqSkpMBgMGD06NGoqOD9JkRypasxosZgnl3iPXato1AIaB9iOVqM73vu6sej2Xjj5xN48eujuGNlKjLyyqUuiTyU7E/m/umnn+p9vGbNGkRERGD//v0YNmyYRFURUXMsPey0KgW81UqJq3F/CaE+yMgrR1ahDkM7S10N2UsURazekWn9+MC5Yoz/zw48M6oLHhraAUoFz1Imx5H9jN2flZSUAABCQkIkroSImnLljlhB4D9arWU9M5Y7Y93SH+eKcOh8MTQqBb574noM6xKOGoMJi388ztk7cjjZz9hdSRRFzJkzB9dffz2SkpIavaa6uhrV1dXWj0tLSwEAer0eer3eJXW2ZZYx5lg7n5zHOr/UfPxVkLdalvXZS+qxjgvSAgAyC8o9YjybI/VYO8OqbacBAJN6RaNbpA/+O6M3vvjjEhb9eMI6ezf7xo54YEiCS2fvPHGs5cxV4yyIoug2d3E+/vjj2LRpE3bu3InY2NhGr1mwYAEWLlzY4PF169bBx8fH2SUSEYC9+QLWZijRJdCEx3twJ2drHS8WsCJdiShvEfN6G6Uuh+xQUAW8ckAJEQKe72VA9BX/DBVVAxvOKJBebF48i/cTMb2TEZHeEhVLTqXT6TBt2jSUlJQgICDAad/HbYLdk08+ia+//hrbt29HYmJik9c1NmMXFxeH7OxshIaGuqLUNk2v1yMlJQWjRo2CWs3dkM4k57Fek3oWi348gQnJUVg6pafU5bSa1GN9vkiHG97cCY1KgSP/uBEKD74nS+qxdrSXNx3HR7vPYWinULx/X78GnxdFERsPXMKrP5xAebUBGpUC/xjfDVOvbXzywpE8bazlrrCwENHR0U4PdrJfihVFEU8++SS++uorbN26tdlQBwBarRZarbbB42q1mj+4LsTxdh05jnVplXlWKcxPK7vaWkOqsW4f6g+1UkCNwYTCSiPaBXn+lI4cf67tVVKpxxd/XAQAPDy8Y5Ov5+6BCRjRLRLPbzyCbSfzsfD7dEzqG4sAL9e8fk8Ya3fgqjGW/eaJxx9/HGvXrsW6devg7++PnJwc5OTkoLKyUurSiKgJdT3s2OrEEVRKBeKC2fLE3Xy65xx0NUZ0jfTH9Z3Cmr02OtAbH9x/LWKCvGEwiTh2sdRFVZKnkX2wW7FiBUpKSjBixAhER0db/2zYsEHq0oioCdZTJxjsHCa+9mixrAIeLeYO9EYTPkjNAgA8ODTRpt3hgiAgOSYQAHD0YokzyyMP5hZLsUTkXix97IJ5nJjDmFue5OMsZ+zcwg9HspFdUoUwPy0m9W5n89clxQTgp2M5OMJgRy0k+xk7InI/V/axI8dICOVSrLsQRRHv7TQ3JL53cDy0KtubdCdZZuwuMdhRyzDYEZHDXeZSrMPFh9U2KS7kUqzc7c0qwuELJdCqFJg+sL1dX2tZis0sqEB5tcEZ5ZGHY7AjIocymUQU6cyNOBnsHCeh9vSJrMIK3qIic6t3nAEA3N4vFqF+Dbs0NCfUT4t2gV4QReAYl2OpBRjsiMihyqoMMJrMwSPIhy0UHCUmyBtKhYAqvQl5ZdVX/wKSRGZBBf6XngsAeOC65ttzNcWyHMv77KglGOyIyKEsrU78tCq77i2i5mlUCsTU9q/L4pmxsrXmt0yIInBDtwh0ivBr0XNwZyy1BoMdETmUdUesL2frHM3S8oT32clTsa4Gn++7AAB46PqWzdYBnLGj1mGwIyKHsvawY6sTh7vyPjuSn3V7zqFSb0T36AAM7tjyIywtwe4MN1BQCzDYEZFD8dQJ5+GMnXzVGEz4sLYh8UPX29aQuCnh/lpEBZg3UKRn8wQKsg+DHRE5FGfsnCcxjDN2crXpyCXkllYjwl+Lib1sb0jcFOty7AUux5J9GOyIyKE4Y+c88aF1vezY8kQ+RFHE6u3mhsT3DUmARtX6f1q5gYJaSvZHihGRe+E5sc4TF+INQQDKqw0oKK9BuL99PdLIsUwmEVtP5mHltjNIyy6Ft1ppd0PipiTHBgDgBgqyH4MdETnU5Qpzc2KeE+t4WpUS7QK9cbG4EmcLKxjsJFJjMOHbQ5ewavtpnMwtBwCoFAKeG9sVQQ76uU9qZ56xO51fDl2NAT4a/nNNtuFPCrV5oiiixmhizzUHqTsnlu1OnCEhzAcXiyuRVahD/4QQqctpU8qq9Fi/5xze35mFnNIqAOZ+jdMGtsf91yUgOtDbYd8rIsALEf5a5JVVI+1SKf9bk80Y7KhN234yH69sSsOl4iqsuf9aXMs3z1azLMVyxs454kN98VtGIc5yA4XL5JZW4f3fMrFu9zmU1bYfifDX4oHrEzFtYHsEeDnnl5jkmED8cjwPRy+WMNiRzRjsqE3KyCvHoh/S8evxPOtjj3/yBzbNGsrlrVa6rOM9ds6UUNvyJMsBLU8uFlfisbX70TM2CM+N6wY/rXP/SThXqEON0dTiExlcyWQS8dvpAny65zw2p+VAbzRvVukU4YeHh3bApD7tnD7Ln1Qb7I5cZMsTsh2DHbUpxboa/PuXU/h411kYTCJUCgH3DI7HjlMFyMgrx6z1B/DxgwOgUnLDeEsYjCaUVNbeY8dg5xR1O2NbP2P38a6zOHShBIculGDryTz8845eGNih5Y11m1Osq8HEt3eissaIrx4fgmtq7yGTm+ySSny+7wI27D2Pi8WV1sf7xwfjkeEdcUO3CCgULe9RZw/ujKWWYLCjNkFvNOGT3Wex9JdTKNaZg8eN3SLw9wnd0THcDxl5Zbjl7d+w60wh3kw5iWfHdpO4YvdUUqmHpQtHkDfvsXMGy+kTmQUVEEWxxY1wRVHEj0ezAQC+GiXOX67E1NW78eB1iZg7piu81I6djfow9aw19D/z2SF888R1srmvVW804dfjediw9zy2nsiDqfZnOMBLhdv6xGDKtXGSBNHkWPP3PJVXhsoaI7w18hgvkjcGO/J4W07k4ZXv03A63zzD0TXSHy/e3B1DO4dbr+kU4Y8lt/fErPUHsHzrafRtH4ybekRKVbLbspwTG+it5qynk7QPMS/FllUZUKzTt3hmND27DGcLddCqFNgydwTeTDmJT/eex393ZmLryXy8OaUXesYGOaTmimoD1qSa+7yplQKO55Rh6f9O4TmJf4EqqdRj5bbT+HzfBRSUV1sfH5gYgqkD4jAuKdrhAdceEf5ahPlpUVBejbTsUvSLD5asFnIffOclj/Z/3xzF/Wv24nR+BUJ8NXjl1iRsmnV9vVBncUuvdrhvcDwAYM5nB3GOxzbZzRLsQrkM6zTeGiWiArwAtO4Eip9qZ+uGdwlHRIAXltzeE+//pT/C/bXIyCvHbctT8VbKSeiNplbXvO73cyjW6ZEY5ov/TO0DAHh322nsP1vU6uduqYpqA+57fw9WbD2NgvJqhPlpMHN4B/z6zHBsmDkYt/WJlTTUAYAgCEiOMfez43Is2YrBjjxWiU6Pj3efBWA+u3HL3BGYMSi+2ZmkFyb0QO+4IJRWGfDYuv2o0htdVa5HKOKpEy7hiDNjfziaAwAYlxxlfeyGbpHYPHsYJvSMhtEk4t+/nMLk5ak4lVvW4u9TpTdi1Y4zAIBHh3fEuORoTO4bA5MIzP38EHQ1rj/kvsZgwiNr9+Pg+WIE+aixYnpf7Jp3I+aN644O4fLa2GG5z46NislWDHbksXZnFkIUgY7hvnjx5h4ItOGeL41KgeXT+yLYR42jF0ux8LtjLqjUc7A5sWtYzow9dKG4RV+fkVeGjLxyqJUCbuhW/5aDYF8N3pnWF/+5uw8CvdU4crEEE5btxIa951r0vb7YfwH5ZdVoF+iFW/vEAADmT7wG0YFeyCyowOs/nWjR87b0SDWTScQznx/CjlMF8FYr8f5frsW45GioZXrrQBI3UJCd5PmTTOQAu04XAgCGdAyz6+vaBXnjP3f3gSAA6/ecxxf7LzijPI/E5sSuMar2/s+vDlxs0azyj0fMs3XXdQpr8heeW3q1w+anh2F4l3DUGEx44aujdocLg9GEldtOAwAeHtbBeoZqoLcar93eEwDwQWoWfssosPk5i3U1ePijfei1cDPW/X7OroAniiIWfncM3x26BJVCwIoZfdG3vbzvW6vbQFHOFQSyCYMdeazU0+Z/LIZ0tL99w9DO4Zh9YxcAwAtfHUF6NvtI2cJyjx2XYp1rRNcItAv0QrFOj59ql1Tt8WPt14xPim72usgAL3xw/7UYnxwFg0nEnM8Ootpge7j49tAlXCiqRKivBnddW/8M1WFdwjFjkPmxZ784jNIq/VWf7+jFEty8bCc2p+WitMqAv391BI+s3W9tin01y37NwIe7zLdn/GtKL4zoGmHza5FKVIAXQn01MJpEvg+RTRjsyCPllVXhZG45BAEY1MK+XE/e0AnDu4Sj2mDCo2v3N/oPj8FowvnLOuw4lY+1u89i9fYzyCuram35bsvyD2wIl2KdSqkQMHWAORR98vtZu772bGEF0rJLoVQI1pm/5giCgJcnJSHMT4OTueV4K+WUTd/HZBKxfKt5tu7BoYmNtuqYN6474kPNR6S9/F1as8+3Ye85TF6RigtFlWgf4oMnRnaCWing52O5GPfvHUi9yqzf2t1n8WbKSQDA/Ik9MKl3jE2vQ2qCIHA5luzCdifkkSzLsD2iA1o8e6RQCFh6V2/cvGwnsgp1mP3pQQztHIazhTpkFVbgXKEO54t01o70Fv/deQbv3tMfveOCWvsy3M5lbp5wmbuujcO/fzmFvVlFOJlbhi6R/jZ9nWW2blCHEJv/O4X6afHqbcmY+fF+rNp+GqN6RKBffPNHXG1Oy0FGXjn8vVSYMSi+0Wt8tSr8885emPLuLny+/wLGXBOF4Z3rP2+V3oj53xzDhn3nAQA3dY/Av6b0RqC3GmOTojDr0wM4k1+B6e/9jpnDOmLOqC7WJV+LH45k4x/fHAVg/oXt/usSbXrdcpEcE4htJ/O5gYJswhk78kh199e1rot+sK8G70zvC7VSwK/H87DwuzR8kJqFrSfycaagAnqjCI1SgY7hvrixWwQ6hPsit7QaU97dhY1t8N48zti5TmSAF27qbl5KXPe77RsbLMFu3FWWYf9szDVR1t2sz3zW/G5WURTxzhbzbN1fhiQ0e5bqtQkheHhoBwDA818esS7nA8D5yzrcsTIVG/adh0IA/jamK1bd0996X2BSTCC+f/J63D2gPUQRWLntNO5YmYrMgro2ML9lFGD2pwchisC0ge0xZ1QXu163HCRZd8ZyKZaujjN25JFSW7hxojG944Lwzzt7Yc1vWYgO9EJ8qC/iQ31q//giKsALytojhsqq9Hh6wyH8Lz0Xz3x+COnZpXh+XLc206yXM3auNX1gPH4+louNf1zAc2O7XfVkgovFlTh0vhiCAIy+xv4G3PMnXoPUjEJkFerw+k8nsOCWaxq9bvupAhy5WAJvtdKm2bGnR3XBlhN5OJlbjgXfpWOMP7D1ZD7mfnEUJZV6hPhq8J+pfXB954b/f/bRqLB4cjKGdwnH818exuELJZjwnx1YcMs16Bblj4c/2ocaownjkqLw8qSkFp/UISXrBorcMlTpjZL31yN5Y7Ajj3P+sg7nLuugVAi4NrH55SJbTeodY9M9Of5eaqy6px+W/u8k/vNrBv67MxMncsuw7O4+CGoDs1hFte1OQhjsXOL6TmFoH+KDc5d1+O7wJUzpH9fs9ZaNFtfGhyDC38vu7xforcbrd/TEve/vwQepWRjdIxJDOjUMW+9syQAA3D2gvU0/C15qJf51Z2/ctvw3/HgsFxdDFDi8+wAA8y9Wy6f3Rbsg72afY2xSFHrFBWLOhkPYdaYQz35xGBqlAjVGE4Z0DMXSqb2tv4C5m3aBXgjx1eByRQ1O5JShVxu8zYNs1zamEahN2XXGPFvXKzYQflrX/+6iUAiYM7orlk/vC2+1EjtOFWDSO7/hZCuavLqDaoMR5dXm5TkuxbqGQiHg7tpNFLYsx1pOmxibFHWVK5s2rEs4pg00f8+/fXEYZX/aVLQ36zL2ZF6GWing4WEdbH7e5NhAPHFDJwDA4cvmf5ruHRyPz2YOvmqos4gO9MbahwbiubHdoFIIqDGakBwTiFX39pfNubQtIQgCrmlnPoGC99nR1TDYkcex7I67rpGZBFcanxyNjY8OQWywN84W6nDbO79h8zH7W1O4i2Kd+R94pUKAvxcXA1zlzv6xUCsFHDxfjGOXmv5HP6+0Cvtqj/BqTbADgL+P7464EG9cLK7Eq5vS633OMlt3R79YRAXaNyv4+MhOGJgYDK1CxD/vSMZLk5IabIS4GqVCwKMjOuLrx6/DnFFd8OEDAyT5Bc/RkrkzlmzEYEceRRRF6/11g1u5ccIRerQLwLdPXI9BHUJQUWPEwx/vx7JfTrW4a74jiaKIjLxyrN19Ft8fvgSTqXU1WXvY+aihcNMlL3cU5qfFmGvMQa25Wbufj+VAFM1Lm7bOgDXFT6vCG3f0giAAn+49jy3H8wCYQ8fWE/lQCMAjwzva/bxqpQIf/qU/Fl1rxKRe9m3u+LOkmEDMurGzx9wWwKPFyFYMduRRTudXIK+sGhqVQjYd5UN8Nfj4wYG4b7C55cO/Uk7ivjV7cam40qV1iKKIc4U6bNh7Dk99egADFv2Cm97chhe/Poon1h3AXat21dtNaK8ia7DzjH9I3YllafTrAxety+F/VrcbtnWzdRaDOoTigdqNEc9tPIxiXQ2WbzXP1k3s1Q7xob4tel6lQoCdk3RtgmVn7MncMruaRFPb4/7z00RX2FV72kT/+GBZ7RxTKxVYOCkJPdoF4B/fHMP2k/kY/dZ2vDihO+66Nq5FO/UuV9QgtxLIyCuHSqWCCMAyEShChCgCJlHEydwypGYUIvV0IS7+KUxqVQr0aR+EIxdKsDerCOP+vR1/G9MN9w9JsHvWjTtipTO4Qyg6hPniTEEFvj14yRr0LArLq/F75mUA9rc5ac7fxnTFlhN5OJNfgUfX/oHdmebZ8sdGdHLY9yCz2GBvBPmoUazT40ROGXrGBkldEskUgx15lFQH9a9zlruubY9+8SF49otD+ONcMZ7/8gg2HcnGktt7IsbG5bGD54uxescZ/HgkGyZRBRxMtfn7qxQC+rQPwuAOoRjcMQx92gfBS63E+cs6PP/lYfyWUYiXv0/DT0ez8fodvayHzdviMnvYSUYQBEwb2B6vbErHJ7+fxd0D6v+ykJKWC6NJxDXtAtA+1Mdh39dLrcSbU3pj8vLfrJuWRvWIRNco25olk+0EQUByTCB2nCrA0YulDHbUJAY78hgmk2j9x2WwA/rXOUunCD98/sgQrPktE2/8fAI7ThVgzFvb8cKE7pjaxOyd0STif+m5+O+OM9ibVWR93FspQqNRQxAECMAV/wsAAgQBiA70wuCOoRjSMQz944Ph28iN5HEhPlj74ECs33Mer25Kw96sIoxduh1/G9MV91+XaFObCJ4TK63b+8bi9Z9P4NilUhy+UFKvJYajl2Gv1DsuCI+N6IS3azdNPD6Ss3XOck07c7DjfXbUHAY78hjpOaUo1unhp1WhV21DT7lSKgQ8NLQDbugWgb99cRj7zxZh3pdH8MORbCyenIzYYPOsiq7GgC/2X8D7OzORVagDAKiVAm7pFYP7BsUh88AOjB8/Bmp10539bWWZ9RnWJQzzvjyCHacK8MqmdPx4NAdv3NETHcL9mv16yz12oQx2kgj21WBCcjS+OnAR634/Zw12JTo9UmtvURjrwGXYK826sTNyS6sQFejVJo/ScxXujCVbMNiRx0jNMM/WDUgMcZuTHjqE++GzmYMbzN79bUxXFJTXYO3vZ61tRAK91Zg+sD3uG5KAyAAv6PV6ZB5wfE2xwT746IEB+HTveby6KR37zxZh3L934IHrExEd6AWtSgGNSgGtSlnv72dqN15wxk460wa2x1cHLuLbQ5fwws3dEeClxv/Sc6E3iugS6YdOEc2H85bSqBR4485eTnluqmMJdidyylBjMNndCobaBgY78hiWWQm53l/XFMvs3Y3dI/G3zw9h39kiLPguzfr5+FAfPHh9Iu7oFwsfjWv+LysI5sa3w7qE4/mNh7HjVAFWbD1t09cG+7R+9pBapn98MLpE+uFkbjm+PnAR9w5OsC7DOmu2jlwnLsQbgd5qlFTqcTK3zLpTluhKDHbkEfRGE/bU7vqTQ/+6lkgM88WGmYPxYWoW3ko5ia5R/nhoaAeM6hEp2VFIMUHe+OiBAfjqwEXsPFWAaoMJ1QZj7f+a/9RYHtObEOKrwdDO4ZLUSrXL6QPaY8F3afhk9zlM7huL7afyATjn/jpyLUEQkBQTgN8yCnHkYgmDHTWKwY48wuELJaioMSLIR43uUQFSl9NiSoWAB65PxP3XJcjmsHJBEDC5bywm942VuhSywW19Y7Hkp+M4kVuGf/58AjUGExLDfNGNO1U9QlJMoDXY3S11MSRLXKAnj2DpXze4Q6hHnHogl1BH7ifQW42JPdsBAD5IzQJgPkKMP1OeIamdeZbuGDdQUBMY7MgjyL1/HZErTR8UX+9jLsN6DssGivTaDRREf8ZgR26vSm+0Hm4+pJN8+9cRuUqv2ED0iDbfkhAT5G0NA+T+4kN9EOqrQY3BhDW/ZUpdDskQgx25vT/OFqHGYEJkgBYd7DgpgchTCYKAR0d0BGBugcJlWM8hCAL+NqYrAOBfm0/iZG6ZxBWR3DDYkdurW4YN4z9gRLUm9mqHP/4xCo/VBjzyHHddG4eRXcNRYzRhzmcHoTdySZbqMNiR27P0r3PXNidEzhLiq+EvOx5IEAQsub0nAr3VOHqxFG//miF1SSQjDHbk1sqrDTh0wbw7jBsniKitiAzwwsu3JgEA3t6SgcMXiqUtiGSDwY7c2t7MyzCaRLQP8bGer0pE1BZM7BmNCcnRMJpEzPnsEKr0RqlLIhlgsCO39luGeRn2uk6crSOitkUQBLx8axLC/LTIyCvHvzafkLokkgEGO3Jrlo0TgzuyzQkRtT0hvhosmZwMAPjvzkzr0YrUdjHYkdsqqqhBWnYpAPOJE0REbdFNPSJxZ79YiCIw9/NDqKg2SF0SSchtgt3y5cuRmJgILy8v9OvXDzt27JC6JJLY7jPm2boukX4I99dKXA0RkXT+b2IPxAR549xlHRb9kC51OSQhtwh2GzZswOzZs/HCCy/gwIEDGDp0KMaNG4dz585JXRpJ6Mr+dUREbZm/lxpv3NETAPDJ7+ew7WS+xBWRVNwi2L355pt48MEH8dBDD6F79+5YunQp4uLisGLFCqlLIxfSG03ILKjAluN5eH9nJlLScgGwfx0REWA+UvEvQxIAAM99cRglOr20BZEkVFIXcDU1NTXYv38/nn/++XqPjx49GqmpqTY/z6/H8+AfyB9yZzMaDThUKECVlgulsjU/XiJyS6uRWVCBrMIKZBVU4EJRJQwmsd5VGpUCgxIZ7IiIAOC5sd2w7WQ+Mgsq8OzGQ7itT2yT1zru/ZoAoHu0P+JDpT/WUvb/JQsKCmA0GhEZGVnv8cjISOTk5DS4vrq6GtXV1daPS0vNN9c/88VRKLTsc+YaSrx/8pBTntlLrUBCiA/iQ32QEOqLoZ1D4aMG9Pq2F9otr7ktvnZX41i7Dse6dVQC8NrkazB19R78fCwXPx/LvcpXOO/9uq35x4RuuHdQ+yY/76qfadkHO4s/H4sjimKjR+UsXrwYCxcubPB4vJ8IlZfY4HGSL1+ViAgvINxbRLgXEO4lIkADKIQaAMWAAShMP4W2fp9wSkqK1CW0GRxr1+FYt84diQL25rvF3VYe48KpY/jh8tEmP6/T6VxSh+yDXVhYGJRKZYPZuby8vAazeAAwb948zJkzx/pxaWkp4uLi8PnjwxAayiU7Z9Pr9UhJScGoUaOgVqulLsejcaxdh2PtOhxrxxhvwzUca9cqLCx0yfeRfbDTaDTo168fUlJScNttt1kfT0lJwaRJkxpcr9VqodU2bH2hVqv5g+tCHG/X4Vi7DsfadTjWrsOxdg1XjbHsgx0AzJkzB/fccw/69++PwYMHY9WqVTh37hweeeQRqUsjIiIikg23CHZ33XUXCgsL8dJLLyE7OxtJSUn44YcfEB8fL3VpRERERLLhFsEOAB577DE89thjUpdBREREJFvcMkNERETkIRjsiIiIiDwEgx0RERGRh2CwIyIiIvIQDHZEREREHoLBjoiIiMhDMNgREREReQgGOyIiIiIPwWBHRERE5CEY7IiIiIg8hNscKdZSoigCAMrKyqBWqyWuxvPp9XrodDqUlpZyvJ2MY+06HGvX4Vi7DsfatcrKygDU5RJn8fhgV1hYCABITEyUuBIiIiJq6woLCxEYGOi05/f4YBcSEgIAOHfunFMHksxKS0sRFxeH8+fPIyAgQOpyPBrH2nU41q7DsXYdjrVrlZSUoH379tZc4iweH+wUCvNthIGBgfzBdaGAgACOt4twrF2HY+06HGvX4Vi7liWXOO35nfrsREREROQyDHZEREREHsLjg51Wq8X8+fOh1WqlLqVN4Hi7DsfadTjWrsOxdh2OtWu5arwF0dn7bomIiIjIJTx+xo6IiIiorWCwIyIiIvIQDHZEREREHoLBjoiIiMhDeESwKyoqwj333IPAwEAEBgbinnvuQXFxcbNfI4oiFixYgHbt2sHb2xsjRozAsWPH6l0zc+ZMdOzYEd7e3ggPD8ekSZNw/PhxJ74S+XPGWF++fBlPPvkkunbtCh8fH7Rv3x6zZs1CSUmJk1+NvDnr53rVqlUYMWIEAgICIAjCVZ/TEy1fvhyJiYnw8vJCv379sGPHjmav37ZtG/r16wcvLy906NABK1eubHDNxo0b0aNHD2i1WvTo0QNfffWVs8p3O44e72PHjuH2229HQkICBEHA0qVLnVi9e3H0WK9evRpDhw5FcHAwgoODcdNNN2HPnj3OfAluw9Fj/eWXX6J///4ICgqCr68vevfujY8//tj+wkQPMHbsWDEpKUlMTU0VU1NTxaSkJPHmm29u9muWLFki+vv7ixs3bhSPHDki3nXXXWJ0dLRYWlpqvebdd98Vt23bJmZmZor79+8XJ06cKMbFxYkGg8HZL0m2nDHWR44cESdPnix+++23YkZGhvjLL7+InTt3Fm+//XZXvCTZctbP9VtvvSUuXrxYXLx4sQhALCoqcvIrkZdPP/1UVKvV4urVq8W0tDTxqaeeEn19fcWzZ882ev2ZM2dEHx8f8amnnhLT0tLE1atXi2q1Wvziiy+s16SmpopKpVJctGiRmJ6eLi5atEhUqVTi7t27XfWyZMsZ471nzx5x7ty54vr168WoqCjxrbfectGrkTdnjPW0adPEd955Rzxw4ICYnp4u3n///WJgYKB44cIFV70sWXLGWG/ZskX88ssvxbS0NDEjI0NcunSpqFQqxZ9++smu2tw+2KWlpYkA6r2B7tq1SwQgHj9+vNGvMZlMYlRUlLhkyRLrY1VVVWJgYKC4cuXKJr/XoUOHRABiRkaG416AG3HlWH/22WeiRqMR9Xq9416AG3HFWG/ZsqVNBrsBAwaIjzzySL3HunXrJj7//PONXv/ss8+K3bp1q/fYzJkzxUGDBlk/njJlijh27Nh614wZM0acOnWqg6p2X84Y7yvFx8cz2NVy9liLoigaDAbR399f/PDDD1tfsBtzxViLoij26dNHfPHFF+2qze2XYnft2oXAwEAMHDjQ+tigQYMQGBiI1NTURr8mMzMTOTk5GD16tPUxrVaL4cOHN/k1FRUVWLNmDRITExEXF+fYF+EmXDXWgPmw5ICAAKhUHn+ccaNcOdZtSU1NDfbv319vjABg9OjRTY7Rrl27Glw/ZswY7Nu3D3q9vtlr2vq4O2u8qSFXjbVOp4Ner3f6QfZy5oqxFkURv/zyC06cOIFhw4bZVZ/bB7ucnBxEREQ0eDwiIgI5OTlNfg0AREZG1ns8MjKywdcsX74cfn5+8PPzw08//YSUlBRoNBoHVe9enD3WFoWFhXj55Zcxc+bMVlbsvlw11m1NQUEBjEajXWOUk5PT6PUGgwEFBQXNXtPWx91Z400NuWqsn3/+ecTExOCmm25yTOFuyJljXVJSAj8/P2g0GkyYMAHLli3DqFGj7KpPtsFuwYIFEASh2T/79u0DAAiC0ODrRVFs9PEr/fnzjX3N9OnTceDAAWzbtg2dO3fGlClTUFVV1cpXJy9yGWsAKC0txYQJE9CjRw/Mnz+/Fa9KnuQ01m2ZvWPU2PV/fpzj3jRnjDc1zplj/frrr2P9+vX48ssv4eXl5YBq3Zszxtrf3x8HDx7E3r178eqrr2LOnDnYunWrXXXJdp3riSeewNSpU5u9JiEhAYcPH0Zubm6Dz+Xn5zdIxxZRUVEAzAk6Ojra+nheXl6Dr7HsSOzcuTMGDRqE4OBgfPXVV7j77rvtfUmyJZexLisrw9ixY+Hn54evvvoKarXa3pcie3IZ67YqLCwMSqWywW/VzY1RVFRUo9erVCqEhoY2e01bH3dnjTc15Oyx/uc//4lFixbhf//7H3r27OnY4t2MM8daoVCgU6dOAIDevXsjPT0dixcvxogRI2yuT7YzdmFhYejWrVuzf7y8vDB48GCUlJTU2379+++/o6SkBEOGDGn0uRMTExEVFYWUlBTrYzU1Ndi2bVuTX2MhiiKqq6sd8yJlQg5jXVpaitGjR0Oj0eDbb7/12N8G5TDWbZlGo0G/fv3qjREApKSkNDlGgwcPbnD95s2b0b9/f+svH01d09bH3VnjTQ05c6zfeOMNvPzyy/jpp5/Qv39/xxfvZlz5c92izGHXVguZGjt2rNizZ09x165d4q5du8Tk5OQGbSG6du0qfvnll9aPlyxZIgYGBopffvmleOTIEfHuu++u1xbi9OnT4qJFi8R9+/aJZ8+eFVNTU8VJkyaJISEhYm5urktfn5w4Y6xLS0vFgQMHisnJyWJGRoaYnZ1t/dPWW8s4eqxFURSzs7PFAwcOiKtXrxYBiNu3bxcPHDggFhYWuuy1ScnSpuC9994T09LSxNmzZ4u+vr5iVlaWKIqi+Pzzz4v33HOP9XpLm4Knn35aTEtLE997770GbQp+++03UalUikuWLBHT09PFJUuWsN1JLWeMd3V1tXjgwAHxwIEDYnR0tDh37lzxwIED4qlTp1z++uTEGWP92muviRqNRvziiy/qvTeXlZW5/PXJiTPGetGiReLmzZvF06dPi+np6eK//vUvUaVSiatXr7arNo8IdoWFheL06dNFf39/0d/fX5w+fXqDFg4AxDVr1lg/NplM4vz588WoqChRq9WKw4YNE48cOWL9/MWLF8Vx48aJERERolqtFmNjY8Vp06Y12WqirXDGWFvabjT2JzMz0zUvTIacMdaiKIrz589vdKyvfB5P984774jx8fGiRqMR+/btK27bts36ufvuu08cPnx4veu3bt0q9unTR9RoNGJCQoK4YsWKBs/5+eefi127dhXVarXYrVs3cePGjc5+GW7D0eOdmZnZ6M/wn5+nLXL0WMfHxzc61vPnz3fBq5E3R4/1Cy+8IHbq1En08vISg4ODxcGDB4uffvqp3XUJolh79x4RERERuTXZ3mNHRERERPZhsCMiIiLyEAx2RERERB6CwY6IiIjIQzDYEREREXkIBjsiIgcwmUx4+OGHER0djYcffhgmk0nqkoioDWKwIyJygJ9//hknT57EDz/8gPT0dPz8889Sl0REbRCDHRG1WTU1NejUqRN+++23Vj9XYGAggoOD0blzZ4SGhiIkJOSqX1NdXY327dtj//79rf7+REQAgx0RtWGrVq1CfHw8rrvuOru/dsGCBZg6dar14yFDhqCmpgaBgYEwGo0YOHCg9XPl5eUYN24chg0bhq5du+KDDz4AAGi1WsydOxfPPfdcq18LERHAYEdEbdiyZcvw0EMPtehrv/32W0yaNMn6sV6vx759+/Dss89iz549MBgM1s/5+Pjg22+/xfbt2/Hee+9hxYoV1s9Nnz4dO3bsQHp6estfCBFRLQY7InJ7hYWFmDp1KoKDgyEIQr0/ltmxP/vjjz+QkZGBCRMm1Hv84sWLuOuuuxAcHIzQ0FBMmjQJWVlZ9a45f/48jh49inHjxlkf27RpE9RqNV566SWoVCps2rTJ+jmFQgG1Wo38/HzMnz8fS5cutX4uNDQUQ4YMwfr161s9DkREDHZE5PZmz56N1NRUbNiwAWlpadZZuGXLlmHYsGGNfs327dvRpUsXBAQEWB/T6XQYOXIk/Pz8sH37duzcuRN+fn4YO3YsampqrNd9++23GDZsGIKCgqyPrVmzBnfffTfUajXuvvturFmzpt732717N+6++24sXboUgwcPrve5AQMGYMeOHa0dBiIiBjsicm+lpaX45JNP8MYbb2D06NHo3r07VqxYgZiYGOj1enTo0KHRr8vKykK7du3qPfbpp59CoVDgv//9L5KTk9G9e3esWbMG586dw9atW63XffPNN/WWYfPy8vDDDz9gxowZAIAZM2Zg06ZNyMvLs35+6NChyMnJwV//+leMHTu23veNiYlpMCtIRNQSKqkLICJqjdOnT0MURQwZMsT6mEqlwoABA3D48OEmv66yshJeXl71Htu/fz8yMjLg7+9f7/GqqiqcPn0agDlIbtu2DatXr7Z+/uOPP0a3bt3Qq1cvAEDv3r3RrVs3rF27FnPmzEFERAT0en2TtXh7e0On09n+oomImsBgR0RuTa1WAwCMRmO9x41GI5RKZZNfFxYWhiNHjtR7zGQyoV+/fvjkk08aXB8eHg4A+PHHH9G9e3fEx8dbP7dmzRqkpaVBpVLVe641a9Zgzpw5V30Nly9ftj4/EVFrMNgRkVvr2LEjvLy88NtvvyEhIQFA3Q7V5kJVnz59sGLFCoiiCEEQAAB9+/bFhg0bEBERUe/euyt98803uOWWW6wf7927F2lpadi6dWu93nXFxcUYNmwY9u3bh/79+zf7Go4ePYo+ffrY+pKJiJrEe+yIyK15e3vjiSeewLPPPouffvoJaWlp+Otf/4qqqio8+OCDTX7dyJEjUVFRgWPHjlkfmz59OsLCwjBp0iTs2LEDmZmZ2LZtG5566ilcuHABBoMBP/74Y73769asWYMBAwZg2LBhSEpKsv65/vrrMXjw4AabKBqzY8cOjB49unUDQUQEBjsi8gCvvvoqpkyZgnvvvRd9+/ZFRkYGfv7553q7Vv8sNDQUkydPrrfs6uPjg+3bt6N9+/aYPHkyunfvjgceeACVlZUICAjAtm3b4Ofnh379+gEw33u3fv163H777Y1+j9tvvx3r169HVVVVk3Xs2rULJSUluOOOO1r24omIriCIoihKXQQRkRSOHDmCm266qdENE42ZNWsWDAYDli9f7rAa7rzzTvTp0wd///vfHfacRNR28R47ImqzkpOT8frrryMrKwvJyclXvT4pKalBD7rWqK6uRq9evfD000877DmJqG3jjB0RERGRh+A9dkREREQegsGOiIiIyEMw2BERERF5CAY7IiIiIg/BYEdERETkIRjsiIiIiDwEgx0RERGRh2CwIyIiIvIQDHZEREREHoLBjoiIiMhDMNgREREReYj/B0tqWM5NFhaGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sigma_moments: [ 8.95401962e+01 -3.52247325e-04  4.47559164e+01  1.76051131e+01\n",
      "  8.17805643e+01  4.88191312e+01  1.85692127e+02]\n",
      "hb_donor_moment: [2.30899094 1.1905562  0.38328986]\n",
      "hb_acceptor_moment: [4.83674399 2.87227988 1.31135965]\n",
      "energy_dielectric: -25.107211968398246\n",
      "dipole_moment: None\n",
      "area_total: 89.54019619110814\n",
      "volume: 68.36643463429641\n"
     ]
    }
   ],
   "source": [
    "def process_molecule(file_path: Path, name: str) -> None:\n",
    "    \"\"\"Process .orcacosmo file to visualize sigma profile and extract descriptors\"\"\"\n",
    "    spp = SigmaProfileParser(file_path)\n",
    "    \n",
    "    # > Generate sigma profile\n",
    "    sigmas, areas = spp.cluster_and_create_sigma_profile()\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.title(name)\n",
    "    plt.plot(sigmas, areas)\n",
    "    plt.xlim([-0.03, 0.03])\n",
    "    plt.xlabel(\"σ (e/Ų)\")\n",
    "    plt.ylabel(\"Area (Ų)\")\n",
    "    plt.grid(True)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "    # > Calculate descriptors\n",
    "    spp.calculate_sigma_moments()\n",
    "    descriptors = {\n",
    "        'sigma_moments': spp['sigma_moments'],\n",
    "        'hb_donor_moment': spp['sigma_hydrogen_bond_donor_moments'][2:5],\n",
    "        'hb_acceptor_moment': spp['sigma_hydrogen_bond_acceptor_moments'][2:5],\n",
    "        'energy_dielectric': spp['energy_dielectric'],\n",
    "        'dipole_moment': spp['dipole_moment'],\n",
    "        'area_total': spp['area'],\n",
    "        'volume': spp['volume']\n",
    "    }\n",
    "    for key, value in descriptors.items():\n",
    "        print(f\"{key}: {value}\")\n",
    "\n",
    "# > Define paths to .orcacosmo files for the solute and solvent\n",
    "cosmo_files = {\n",
    "    'Water': working_dir / (basename + '.solute.orcacosmo'),\n",
    "    'Ethanol': working_dir / (basename + '.solvent.orcacosmo')\n",
    "}\n",
    "\n",
    "# > Loop through and process both solute and solvent\n",
    "for custom_name, file_path in cosmo_files.items():\n",
    "    process_molecule(file_path, name=custom_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17b83733",
   "metadata": {},
   "source": [
    "## Step 6: Predict Solubility of Paracetamol Using openCOSMO-RS\n",
    "\n",
    "We use the openCOSMO-RS to calculate the activity coefficient (ln(γ)) and predict solubility using both:\n",
    "\n",
    "- A non-iterative method\n",
    "\n",
    "- An iterative method that solves the full equilibrium condition\n",
    "\n",
    "The results are printed and compared."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fd6cca0",
   "metadata": {},
   "source": [
    "### Structure of Paracetamol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "69545af3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_17514655796399615\"  style=\"position: relative; width: 400px; height: 400px;\">\n        <p id=\"3dmolwarning_17514655796399615\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n        </div>\n<script>\n\nvar loadScriptAsync = function(uri){\n  return new Promise((resolve, reject) => {\n    //this is to ignore the existence of requirejs amd\n    var savedexports, savedmodule;\n    if (typeof exports !== 'undefined') savedexports = exports;\n    else exports = {}\n    if (typeof module !== 'undefined') savedmodule = module;\n    else module = {}\n\n    var tag = document.createElement('script');\n    tag.src = uri;\n    tag.async = true;\n    tag.onload = () => {\n        exports = savedexports;\n        module = savedmodule;\n        resolve();\n    };\n  var firstScriptTag = document.getElementsByTagName('script')[0];\n  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n});\n};\n\nif(typeof $3Dmolpromise === 'undefined') {\n$3Dmolpromise = null;\n  $3Dmolpromise = loadScriptAsync('https://cdn.jsdelivr.net/npm/3dmol@2.5.0/build/3Dmol-min.js');\n}\n\nvar viewer_17514655796399615 = null;\nvar warn = document.getElementById(\"3dmolwarning_17514655796399615\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nviewer_17514655796399615 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17514655796399615\"),{backgroundColor:\"white\"});\nviewer_17514655796399615.zoomTo();\n\tviewer_17514655796399615.addModel(\"20\\nenergy: -32.783632377928 gnorm: 0.000608672974 xtb: 6.7.1 (edcfbbe)\\nC      2.88225    1.70760    0.08636\\nC      1.87583    0.57255    0.08465\\nH      3.88408    1.29416    0.14629\\nH      2.79182    2.28846   -0.82886\\nH      2.71377    2.36342    0.93735\\nO      2.19537   -0.59390    0.10126\\nN      0.59013    1.02699    0.05895\\nC     -0.60624    0.29815    0.04364\\nH      0.46337    2.02968    0.04757\\nC     -1.80437    1.01047    0.01777\\nC     -1.86108   -1.74887    0.03135\\nH     -1.78549    2.09153    0.01346\\nC     -3.01629    0.35366   -0.00279\\nC     -0.64533   -1.09159    0.05178\\nH      0.27418   -1.65102    0.07346\\nC     -3.05363   -1.03730    0.00247\\nH     -3.94534    0.90006   -0.02339\\nO     -4.26895   -1.65165   -0.02125\\nH     -1.88115   -2.83055    0.03727\\nH     -4.14859   -2.60876   -0.01690\\n\\n\",\"xyz\");\n\tviewer_17514655796399615.setStyle({},{\"stick\": {\"radius\": 0.1}, \"sphere\": {\"scale\": 0.3}});\n\tviewer_17514655796399615.zoomTo();\nviewer_17514655796399615.render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_17514655796399615\"  style=\"position: relative; width: 400px; height: 400px;\">\n",
       "        <p id=\"3dmolwarning_17514655796399615\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n",
       "        </div>\n",
       "<script>\n",
       "\n",
       "var loadScriptAsync = function(uri){\n",
       "  return new Promise((resolve, reject) => {\n",
       "    //this is to ignore the existence of requirejs amd\n",
       "    var savedexports, savedmodule;\n",
       "    if (typeof exports !== 'undefined') savedexports = exports;\n",
       "    else exports = {}\n",
       "    if (typeof module !== 'undefined') savedmodule = module;\n",
       "    else module = {}\n",
       "\n",
       "    var tag = document.createElement('script');\n",
       "    tag.src = uri;\n",
       "    tag.async = true;\n",
       "    tag.onload = () => {\n",
       "        exports = savedexports;\n",
       "        module = savedmodule;\n",
       "        resolve();\n",
       "    };\n",
       "  var firstScriptTag = document.getElementsByTagName('script')[0];\n",
       "  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n",
       "});\n",
       "};\n",
       "\n",
       "if(typeof $3Dmolpromise === 'undefined') {\n",
       "$3Dmolpromise = null;\n",
       "  $3Dmolpromise = loadScriptAsync('https://cdn.jsdelivr.net/npm/3dmol@2.5.0/build/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_17514655796399615 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_17514655796399615\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "viewer_17514655796399615 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17514655796399615\"),{backgroundColor:\"white\"});\n",
       "viewer_17514655796399615.zoomTo();\n",
       "\tviewer_17514655796399615.addModel(\"20\\nenergy: -32.783632377928 gnorm: 0.000608672974 xtb: 6.7.1 (edcfbbe)\\nC      2.88225    1.70760    0.08636\\nC      1.87583    0.57255    0.08465\\nH      3.88408    1.29416    0.14629\\nH      2.79182    2.28846   -0.82886\\nH      2.71377    2.36342    0.93735\\nO      2.19537   -0.59390    0.10126\\nN      0.59013    1.02699    0.05895\\nC     -0.60624    0.29815    0.04364\\nH      0.46337    2.02968    0.04757\\nC     -1.80437    1.01047    0.01777\\nC     -1.86108   -1.74887    0.03135\\nH     -1.78549    2.09153    0.01346\\nC     -3.01629    0.35366   -0.00279\\nC     -0.64533   -1.09159    0.05178\\nH      0.27418   -1.65102    0.07346\\nC     -3.05363   -1.03730    0.00247\\nH     -3.94534    0.90006   -0.02339\\nO     -4.26895   -1.65165   -0.02125\\nH     -1.88115   -2.83055    0.03727\\nH     -4.14859   -2.60876   -0.01690\\n\\n\",\"xyz\");\n",
       "\tviewer_17514655796399615.setStyle({},{\"stick\": {\"radius\": 0.1}, \"sphere\": {\"scale\": 0.3}});\n",
       "\tviewer_17514655796399615.zoomTo();\n",
       "viewer_17514655796399615.render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# > define cartesian coordinates in Angstroem as python string \n",
    "xyz_data = \"\"\"\\\n",
    "20\n",
    "energy: -32.783632377928 gnorm: 0.000608672974 xtb: 6.7.1 (edcfbbe)\n",
    "C      2.88225    1.70760    0.08636\n",
    "C      1.87583    0.57255    0.08465\n",
    "H      3.88408    1.29416    0.14629\n",
    "H      2.79182    2.28846   -0.82886\n",
    "H      2.71377    2.36342    0.93735\n",
    "O      2.19537   -0.59390    0.10126\n",
    "N      0.59013    1.02699    0.05895\n",
    "C     -0.60624    0.29815    0.04364\n",
    "H      0.46337    2.02968    0.04757\n",
    "C     -1.80437    1.01047    0.01777\n",
    "C     -1.86108   -1.74887    0.03135\n",
    "H     -1.78549    2.09153    0.01346\n",
    "C     -3.01629    0.35366   -0.00279\n",
    "C     -0.64533   -1.09159    0.05178\n",
    "H      0.27418   -1.65102    0.07346\n",
    "C     -3.05363   -1.03730    0.00247\n",
    "H     -3.94534    0.90006   -0.02339\n",
    "O     -4.26895   -1.65165   -0.02125\n",
    "H     -1.88115   -2.83055    0.03727\n",
    "H     -4.14859   -2.60876   -0.01690\n",
    "\n",
    "\"\"\"\n",
    "# > Visualize the input structure\n",
    "view = py3Dmol.view(width=400, height=400)\n",
    "view.addModel(xyz_data, 'xyz')\n",
    "view.setStyle({}, {'stick': {'radius': 0.1}, 'sphere': {'scale': 0.3}})\n",
    "view.zoomTo()\n",
    "view.show()\n",
    "\n",
    "# > Write the input structure to a file for reading\n",
    "with open(working_dir / \"paracetamol.xyz\",\"w\") as f:\n",
    "    f.write(xyz_data)\n",
    "\n",
    "basename_paracetamol = \"opencosmors_paracetamol\"\n",
    "structure_paracetamol = Structure.from_xyz(working_dir / \"paracetamol.xyz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a195773",
   "metadata": {},
   "source": [
    "### openCOSMO-RS calculation for Paracetamol:\n",
    "\n",
    "We can re-use the helper functions defined above to run the orca calculation and obtain the `.orcacosmo` file for paracetamol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "32bd3d93",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running ORCA calculation ...   Done\n"
     ]
    }
   ],
   "source": [
    "calc = setup_calc(basename=basename_paracetamol,working_dir=working_dir,structure=structure_paracetamol)\n",
    "output = run_calc(calc)\n",
    "check_termination(output)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1303575e",
   "metadata": {},
   "source": [
    "### Constants and solute specific experimental data\n",
    "\n",
    "Solute specific constants are required and can be obtained for example from [NIST](https://webbook.nist.gov/cgi/cbook.cgi?ID=C103902&Mask=4#Thermo-Phase)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4f4d1e61",
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 8.314 # > Gas constant J/(mol·K)\n",
    "TEMP = 298.15 # > Temperature for the solubility calculation (K)\n",
    "# > Solute specific: \n",
    "DELTA_H_FUSION = 27.1e3 # > Fusion enthalpy of the solute (J/mol)\n",
    "T_FUSION = 443.6 # > Fusion temperature of the solute (K), aka melting point of the solute"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "440be5df",
   "metadata": {},
   "source": [
    "Now the solubilities of paracetamol in water and in ethanol are calculated. The resulting solubility is expressed in mole fractions\n",
    "\n",
    "$x_\\text{solvent}=(\\frac{n_{\\text{solute}}}{n_{\\text{solute}} + n_{\\text{solvent}}})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "65652ba9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solubility results in mole fractions:\n",
      "\n",
      "Solute: Paracetamol\n",
      "\n",
      "Solvent                Non-Iterative       Iterative    Experimental\n",
      "Water                    7.42369e-04     1.94890e-04     1.77260e-03\n",
      "Ethanol                  1.71201e-01     9.11733e-02     7.58600e-02\n"
     ]
    }
   ],
   "source": [
    "def compute_ln_gamma(crs: COSMORS, mole_fractions: list[float], temp: float) -> float:\n",
    "    \"\"\"Compute ln(gamma) for a given composition at specified temperature using openCOSMO-RS\"\"\"\n",
    "    crs.clear_jobs()\n",
    "    crs.add_job(np.array(mole_fractions), temp, refst='pure_component')\n",
    "    return crs.calculate()['tot']['lng'][0][0]\n",
    "\n",
    "def solubility(crs: COSMORS, delta_h_fus: float, t_fus: float, temp: float, iterative: bool) -> float:\n",
    "    \"\"\"Estimate solubility from openCOSMO-RS ln(gamma)\"\"\"\n",
    "    rhs = -delta_h_fus / R * (1 / temp - 1 / t_fus)\n",
    "\n",
    "    if not iterative:\n",
    "        ln_gamma_inf = compute_ln_gamma(crs, [0.0, 1.0], temp)\n",
    "        return np.exp(rhs - ln_gamma_inf)\n",
    "    \n",
    "    def equation(x_guess: float) -> float:\n",
    "        \"\"\"Equilibrium condition to be solved\"\"\"\n",
    "        x_guess = max(1e-15, x_guess)\n",
    "        x_guess = min(1, x_guess)\n",
    "        x = np.array([x_guess, 1 - x_guess])\n",
    "        ln_gamma = compute_ln_gamma(crs, x, temp)\n",
    "        return np.abs(ln_gamma + np.log(x_guess) - rhs)\n",
    "\n",
    "    result = fsolve(equation, 1e-5)\n",
    "    \n",
    "    x_sol = max(1e-15, result[0])\n",
    "    x_sol = min(1, x_sol)\n",
    "    return x_sol\n",
    "\n",
    "def run_solubility_analysis(crs: COSMORS, solute_path: str, solvent_paths: list[str], solute_label: str, solvent_labels: list[str]) -> None:\n",
    "    \"\"\"Run solubility analysis comparing non-iterative and iterative results for given solute/solvents.\"\"\"\n",
    "    if solvent_labels is None:\n",
    "        solvent_labels = [Path(p).stem for p in solvent_paths]\n",
    "\n",
    "    print(\"Solubility results in mole fractions:\\n\")\n",
    "    print(f\"Solute: {solute_label}\\n\")\n",
    "    print(f\"{'Solvent':<20} {'Non-Iterative':>15} {'Iterative':>15} {'Experimental':>15}\")\n",
    "\n",
    "    # > experimental reference data for comparison, taken from https://doi.org/10.1021/je990124v (for T = 298.15 K, converted from g/kg to mole fractions)\n",
    "    exp_ref = {\"Water\" : 0.001772602, \"Ethanol\" : 0.075859954}\n",
    "\n",
    "    for solvent_path, solvent_name in zip(solvent_paths, solvent_labels):\n",
    "        crs.clear_molecules()\n",
    "        crs.add_molecule([solute_path])\n",
    "        crs.add_molecule([solvent_path])\n",
    "\n",
    "        x_non_iter = solubility(crs, DELTA_H_FUSION, T_FUSION, TEMP, iterative=False)\n",
    "        x_iter = solubility(crs, DELTA_H_FUSION, T_FUSION, TEMP, iterative=True)\n",
    "\n",
    "        print(f\"{solvent_name:<20} {x_non_iter:>15.5e} {x_iter:>15.5e} {exp_ref[solvent_name]:>15.5e}\")\n",
    "\n",
    "\n",
    "# > Set up openCOSMO-RS with the 2024a parameterization\n",
    "crs = COSMORS(par=openCOSMORS24a())\n",
    "\n",
    "# > Define file paths and labels\n",
    "solute_file = working_dir / (basename_paracetamol + '.solute.orcacosmo')\n",
    "solute_label = \"Paracetamol\"\n",
    "solvent_files = [working_dir / (basename + '.solute.orcacosmo'), working_dir / (basename + '.solvent.orcacosmo')]\n",
    "solvent_labels = [\"Water\", \"Ethanol\"]\n",
    "\n",
    "# > Run solubility analysis\n",
    "run_solubility_analysis(crs, solute_file, solvent_files, solute_label=solute_label, solvent_labels=solvent_labels)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c13c0abc",
   "metadata": {},
   "source": [
    "For water the non-iterative solver gets closer to the experimental value, while for ethanol the iterative approach agrees better with the experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cf68c7",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook we showed how OPI can be employed together with `openCOSMO-RS_py` to calculate the solubilitiy of paracetamol in water and in ethanol."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "orca-pi",
   "language": "python",
   "name": "python3"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}