{ "cells": [ { "cell_type": "markdown", "id": "2c4c131d-c678-454c-b878-e25905d61281", "metadata": {}, "source": [ "## Initial imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "2b77f49c-4a9f-4803-bdf7-b597f924bc41", "metadata": {}, "outputs": [], "source": [ "import scm.plams as plams\n", "import os\n", "import numpy as np" ] }, { "cell_type": "markdown", "id": "a96ee558-b1d3-4ce2-9a4e-e72e5f18d2f2", "metadata": {}, "source": [ "## Create initial benzene system" ] }, { "cell_type": "code", "execution_count": 2, "id": "b8c1fb9d-71dd-4efb-b3a8-f8c0805fc2ad", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHsAAAB7CAYAAABUx/9/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABhbElEQVR4nO29d3RVxfoG/Mze5+zT03snIQ1CCSGhhhA6GMDQQZCqFC+gwlVARKRcURQVpApSpAtKJ7QQAiEBAimUhJCQ3ntyepvvDzjnEpIo4r2/7/vW9VmLtcjZs2fPnnfPzNtfQinF3/jfAPP/9gD+xv8d/ib2/xD+Jvb/EP4m9v8Q/ib2/xD+Jvb/EP4m9v8Q/ib2/xD+Jvb/EHiv2pAQUgcglVLalxAS10KTzgBSW7n9P3rtd8bwutf+vz6+Pxo7KKVWrYzRDPKq6lJCSN2rdPg3/u/xqrR55ZWN1r/u/1dBCLEHoKGUNvxOGyGlVP1fer4EwHAATgBsOY5zZ1lWpFarEyilRwBYANAAKKT/PUNE6iu1opS+0j8AcS/839rS0vK8RCKpFAgEdwH8E4DDK/QhtLCw2G5jY/PYwsLiRwCiV33+C30QsVi8WiaTFVtYWFSKRCINx3EakUi0sIW2XWUymZJhGCqVShuFQmGtQCBQiMXiT1/hOb4A5gIY+PxvFsBEkUh0l8/nNwgEgnI+n3+VZVmlo6OjUiQSUXt7ezpy5Eg6fvx42qFDBx2Px6NCoVDPMAwVi8V6AMtM4wLwFgDPP/v+f0Sb3233Oh1aWloemTp1qiY8PJxGRETQCRMm6MVicSWATgB4rfUhlUo3Dho0SHnz5k06dOhQpUwm2/IaLzbJz89P/sMPP1APDw/6z3/+k86cOZPyeDwtANfnbaQA3pbJZKqlS5dSo9FIFy9eTIcMGUJzc3Opo6OjHEDv33lGZ5FIJJ8yZYrS1dVVLhQKPxMIBJednJx006dPp6tWraIfffQR7dWrFxWLxVQgENChQ4fSb7/9ln744Ye0W7du1MPDg3p4eFA+n0/Pnj1Lf/nlFyqVSo18Pv+4tbW1IioqqkEkEjUC6PR/RezX4sZZlu04dOhQ7smTJ7hy5QoOHTrEhoaG2rEsm8pxnFYoFD4ghHQ3tSeEMIQQKcdxPRctWiTq0aMHFi1aJOI4rsefeS4hhPB4vA5vvvmmZOvWraivr0daWhoqKirg4eHB5/F4mSzL/sDj8cr9/f03W1hYCHv16gVCCHr37g2NRgNKKTp16kQBeLb2HLFYPHvZsmXiffv2ia5cuSJhGOYzZ2fn/osXL+Z16tQJFhYWcHFxwdixY9GzZ09YWVlh0KBBOH/+PH766Sc4ODhg1KhRGDJkCAIDAzFu3Di4urpCp9MRoVA46sKFC+LTp0/LVqxYIZHJZB+8Dg1eB3/mzDZDrVb/snLlykX19fXirKwsODk5IS8vDxMnToSXlxdJS0trf/ny5QQej/cZy7JtGIaZRAhhGxoayNSpU/Vr167l7dixQ6lWqy+82C8hxB9AKIA0Sun957+xAOaLRKK1AMR6vd64adMmhIaGYubMmWCYZ99r//79cebMGent27ffmz9/Puzt7ZGcnIxZs2Zh3Lhx+Pnnn9HQ0IDg4GAoFAopx3EjCCFXKKVlhBARgHFisfhdSqkjpdT10qVLZOnSpXj48CEYhiGWlpbN5kGpVKKkpASBgYFITk5GWloaFi9eDKlUam7Trl07JCUloX///nBxcUFDQwPWrl2LgwcPoqCgQG8wGJrwGoQQP5lMdk6pVHpJpdLY+vr6Nymlyteh08t4ZW78pQExLMuuY1l2EaWUAYCIiAi88cYbIIQAAKqqqrB+/XowDAOtVouQkBAMHDgQ8fHxhlu3bhkJIVs0Gs1iSqn+eZ/9xGLx6YEDBxpiY2PZxsbGtwGcEAgEJ+zs7IZGRUWxfn5+ePjwIU6dOoVly5aZCQ0ABoMBq1evxowZM+Dh4WH+/fLly7h48SLCw8MxePBgcByH+vp6xMbG6hMTE6u0Wu0QjuNOuLm52fXp00dqaWmJgoICxMbGQqPRgGVZ+Pj4oKqqClZWVnjnnXegVqtx8uRJpKWlwcbGBpRSVFVVoUOHDnjrrbeajMuE1atXw8/PD6tXr8aXX36J27dv6wwGQ5FcLu9BKS03tbOysrq7YsWK4HfffZe89dZb6gsXLnytVqs//dNEagF/Rs6Oo5T2BQBKqVEkErk4OjqSwMBADBgwADxe067s7OwQEREBlmXx66+/YsCAAcjKysLo0aNZS0tLXL58OdhEaACwsbH5ZMOGDeKpU6fi9OnTmDx58jGtVlsmFout33//fZbP5wMAHj9+jF69ejWb0Ly8PEil0iaE1uv1iI+Px5w5c+Dt7W3+3dLSEtHR0TyWZe0SEhKuRkREyIYNG2Z+gTZt2qB79+7YuHEj+Hw+PD09wXEc0tLScPToURQVFcHNzQ2ffPKJeRVXVlbi6NGjOHToECZNmmT+6E2wsLDAkiVL0L9/f4hEIrzxxhtlcrnc78U5AACj0eg+ZMgQIpVKMWjQIOG1a9d8/gxtfg+vrUFjGKaDUqkkQUFBzQhtQseOHZGdnQ2hUIiwsDCkpKQAAPr06cNSSsMIIW1MbXU6nVNeXh4AoLCwEJ6enpgyZYqTUCgU/Prrr8jOzsajR4/Q0NAAiUTS7FkKhQIvbrVGoxG3b9+Gvb19E0K/iH79+vEMBoN1REREsxcQCASYMGEC6uvrERkZicmTJ+OTTz5BWloapFIpxo0b12S7tre3xzvvvIOnT5/C9B4vQiqVYu3atbh48SI+//xz6PV6ewD8F9sQQhi9Xn9t1KhRuhUrVmDZsmWa+vr6PS0O/jXwWmf2c9RTSmEwGFptUF5ejsLCQgwfPhzx8fEghODJkyfw9fWFu7u7Jjs7OwhALiGkPcdx3ps3b8b27duh1Woxc+ZMODs7o6ysDBcuXEB+fj5kMhkKCgpQWloKb29v2NnZmZ9lZ2eH4uJi6HQ6XLt2DQkJCVAqlejVq1er45NKpbCyskJdXV2LH5C7uzsYhkFlZSUcHR1hZWUFqVSK/v37N1m5JSUlUKvV8PT0RK9evZCUlIQ2bdo06UutVqOiogITJ07EjBkzkJKSopfL5e4AsgCAEBLIcdwFqVRqLRQK+QcPHgQhxCgQCLYRQgYBKAXgDyCfUlr9ShR6Ca9NbKVSuUsgEHRPTU3le3o2ZWyLiopw/PhxVFZWwtXVFbdu3YLRaERAQACuX78OX19fKJVKAkAJABzHLezbty9v0KBBqKurg6WlJfh8Ps6fP4/79+/j/fffh4uLCwBAq9Xi+vXr2LhxIz744ANYW1sDAFxcXCCTybBp0yaIRCLMmDEDT58+RUlJSavvYDQaoVQqIRQKW23DMAxe5GsaGxvh7Oxs/vvcuXNIT0837yrh4eHIzMxs0seDBw+Ql5cHW1tbeHl5wdXVlcrlci2AQgAghNjy+fzro0ePtgkLCyOmD4lSKkpISPA8depUEsuy1MnJiSsrK2MJIVGU0rjfJVAL+Csr+6hcLv9XYmKiS+fOnWEieGlpKbZt24aoqCiEhoaCZVkAz7bmvXv3QqFQoLKyEpWVlQRAAgDweLxuvr6+PB6PB7FYjBMnTqCsrAxFRUVYvnw5ZDKZ+aEcx6F///5QKpU4cuQI3nnnHbAsi7KyMqjVavB4PCxYsAByuRyEEKSkpCA8PBxubm7NXiAjIwOUUvMH8zJMfb64g0gkEtTU1EAikUAulyMhIQGFhYWwtLREeHg4Hj58CJ1Oh6qqKiiVSiQlJSElJQXe3t7Izs5GeXk5Hj16pHzOdGkBgGXZdzt27Cju1q1bk4P+ucjIPHr0yFahUNBLly6R9PR0TJ8+fSuAwD9LsNc+symlKp1O102v1xf88MMPOHLkCB4/foxffvkFAwcORPfu3UEIwa1bt3DmzBmo1WrMnTsXOp0Ou3bt0hBCvqb/VmGqNBoNAOD48ePw9fVFZGQkOnfu3ITQLyIiIgLZ2dlYunQpli5diu3bt0OhUGDAgAE4dOgQvvrqK+Tk5MDf3x8//PADtm3bhsbGRvP95eXlOHjwoFGn05XfvHnT+HL/BoMBv/32GyQSiflIAIDg4GBcu3bt5bkA8GynePjwIUpKSrBhwwYcPHgQUqkU//znP+Hi4gI+n4/AwEB06dJF4uDg8C3HccWEkEiO46b37NlT1Npch4eHQy6Xk65du+LBgwcghPBba/t7eOWV3RK3RyktIoR4AXgzKSnp16ysLDQ0NGD27NkAgAsXLqCsrAxjxozBxo0bER0djbZt2yI7O1uu1+vXmPpRqVT7b968GRQUFCQpKCjAL7/8gi+//NK8dbcECwsLCAQCzJ8/H6dOncLs2bPx0Ucf4fr163BycsKKFSvM27NOp0NMTAy++uor9OrVC5WVlUhLS4Onp6exoKCg8tSpU/a5ubno3bs3LC0tUVhYiIsXL4JhGPTr1w+ZmZm4ePEipk6dCoZhkJ6ejtjYWERERKB3794IDAyEhYUF6uvroVKpQAgxrFy5kuU4DsCzD7isrAzLly9/kanjMjIyHH/66aczer2ePnjwANXV1ejYsWOzY0UqlZoZxtWrV0Ov1//rj2jTEl5L9HrpQZQQ0gAA06dPx/79+yEQCAAAaWlpuHLlCtq1awdbW1scPHgQ7u7uyM7ONlJKjYQQP4FA8AkhZOyTJ09Et2/fRps2bbB06VLweDxUVFRAr9fj1q1bSExMRFVVFUQiEbp06YKuXbtCo9HgwIEDKC0txeLFi2FiGMeNG9dENOPz+Rg+fDhKS0uhUCjQv39/ZGZmYu7cubzly5f79uzZU9ulSxfh5s2bYTQaYW9vj379+qFz585gWRZhYWHIycnBjz/+CKPRCKPRiGvXruHatWsIDAyEnZ0dcnJyEBgYCHt7e1hYWDBGoxHXr1/H7du3UVlZCXd3dzx58gQdO3YEy7IoLCzEr7/+CkKI2MrKCgaDAffv38eJEycwcOBA9O3b18wE5uXlwdHRER4eHggJCdEkJyc34f5eVfT6K2c2AIBl2RkikWiTTCZDSUkJGhoaoNVqwXEcbGxscObMGbi4uOD8+fNmbhqAHcMw+3k83sjIyEhRjx49WLlcjh07dsDCwgKNjY2Qy+UoLy9HaWmpmVju7u6or69HYmIivvvuOxBCwHEcampqUFxcjM6dO6N79+4tKjUAIDIyErt27UJqaireeust8Hg8+Pr6clqt1vDll1+ioaEBt2/fxrRp05rd6+Pjg7CwMDx48ABSqRQqlQoNDQ3Iz8+Hi4sLRowYgdDQUFy6dAlqtZps2LABTk5OeOONN2BtbY3i4mLExcXh1q1bmDVrFvbs2YNNmzahf//+eOONN+Dq6opRo0ahuroaO3bsAI/HQ3h4ONRqNa5du4bx48cDADp37ixIT08fBuBPK1r+kqcKIcSFZdnNH3zwgXjw4MGIi4uDu7s77t69CwAYNWoUtm/fDhcXF5SWlqJLly7IycnB9OnTia+v71uWlpbSvn37spaWlnB1dcWnn36Kvn37orKy0lBZWZlFKS0Ri8WYPXs2/P39IRaL4ezsjO7du4PH42HYsGHw9/eHTCaDj48POI5rIvu+DDs7O3Ach1WrVqFdu3YAgLZt25LU1NR6d3d3+f79+2mPHq2r68PCwqBSqbBgwQKcOnUKtra20Ol04PF4SEtLw+effw65XI7U1FR06dIFM2bMQEBAABwdHdGlSxcsWLAAHMfhzJkzqK2txYgRI+Dg4IDIyEjU1tYCAGxtbTFjxgxcuHABWVlZ2LJlC/z8/NC2bVsAZv7gtUylf2ll83i8d0NCQuDg4AB7e3vk5+fj/v37OHXqFMRiMTp27IgFCxYAeKZh2r17N/r164f27dujXbt2OHDgAE6dOoVx48aZ+kNISAjc3NzY9evXOwMQRkdHN1up169fR0REBMLCwvD9999jwoQJKCoqgkgkgkKhaHW8ZWVlsLa2biIjZ2dnK+Vy+Rq5XB4nFouPcxzXsgYGzxQtRqMRH3zwARiGwbhx45CdnY0BAwaY+9+6dSsMBgMGDhzY7P7S0lJ4eXkhJiYGlpaWiIqKQnh4OH788UfMmjXL3M7R0RGWlpbYv38/Bg0aBJMxBwBSUlLUGo3m9B/RpiX8pZXNcVwvDw8P4f79+7F582ZYWVlh7NixcHZ2xqFDh7By5Urs3LkTP/zwA7777juEhoaaJ4EQghEjRuDevXtQqVRN+nV0dASfzxdJpVK+ra1ts+fm5OSgQ4cOkEgkmD9/PvR6PTw9PTFixAjcuHGjRUUPpRRXr16Fj4+P+XppaSkyMjIIgD2U0lSDwXAtMzNT3+zm58jMzIRYLMaMGTPw008/4eeffzavOABwcnJCaGgo2rdvbxY5TUhLS8Pu3bshFArBMAzs7Oxga2uLxMREzJw5s5lo6ODggKioKPTu3dtM6KdPnyIlJcVgMBi2/x5dWsNfWtlGo5FNTEzE+PHjMXLkSLzzzjuwtLTEggULcPv2bVy6dAlarRaRkZHw8/ODSb9tgoWFBWxtbVFWVtZE42Q0GqHX63mtqWGBf4s7MpkM4eHh5vsSExOxe/dujBkzBlZWVgCeqVJPnDiB4uJi1NfXIzk5GU5OTigrK6MSiaSivr5+sVgsjjIYDP5Xr17l1dXVYeDAgU3ka4VCgdjYWPTu3RunTp1CZmYmxowZAy8vrybjsrS0RF1dXbPxJiUlYc+ePYiKisK6deuwZs0azJkzp9lHYXq33NxcsCwLJycn6PV6JCcnG+7cuaPR6XSjKKVlv0uYVvDaohchJIhl2R4VFRVYtmwZrK2tMXXqVFy6dAkhISGQSCQwGAywsrJC+/btW+1Xp9M1e+HHjx+DYRhoNBrk5eU1m1A/Pz+kpaWZRTO1Wg2FQgFra2tMnz4dGzduxL/+9S/Y29uDYRiUl5ejY8eOWL58OUQiEUpKSrB7924MGDCAFBcXe+bn538ybNgweHt7Q6FQIDExEd988w2mTp2KNm3a4MGDBzh9+jT8/f1ha2sLgUCA0tJS8Hg8eHt7m6UPAPD19cX58+dhMBiavBfHcXj69Kl5vEajESkpKejatWuzOXn06BEMBgMYhsGRI0cAACzLEgAKAI/+iDat4bVFL5FItHnw4MGCuLg4/OMf/8DIkSOxefNmDBkyBAAgFovR0NCABw8emLnzl1FUVASVSgVXV1fzb3V1dThy5AgcHR3B4/GwZ88eLFq0qIlypVevXti4cSM6d+6MmpoaHDp0CEKhEDY2Nhg8eDCqq6vh7OwMjUYDW1tbzJ49uwnj5uLiggULFuCLL76Ai4sLFi9ebDbmSKVSvPnmm/D398euXbug1+shFothNBqhVqtx+/ZtdO7cGV988QXWrl2Lc+fOITo62ty3ra0tjEaj4cSJE3TUqFE80xY8ZMgQrFq1CitXroRGo1GrVKo5R44c2aLVasWhoaHg8/nQ6/VITU3Fr7/+iunTp8PX1/fF6WLOnTtnHR8fvx3AG79Hm1Zp+Kr27Bc7JIQ4cxz3dM2aNcKqqir8+OOPZlNg27Zt8fjxYzx48ACenp4oLS1FUFBQM9n3ypUrOH/+PAAgKCgInTp1Qn5+Pm7dugWJRIKCggJwHIeAgAAUFRWhX79+TUQvuVwOtVoNhmFw/vx59OjRA9HR0Th//jyGDx+OmJgYGI1GLFmyxLydv4xDhw6BZVkzg/gyvv32W+j1ekOfPn1Yd3d3PH36FJcuXUKfPn1w8uRJxMXFYd68eZgzZw4AQC6XY9OmTWBZltbW1sLS0lIdFRUlsrKyQklJCS5evCivr69P0el0QwB0YFn2kkwmE2s0GtbKysrMkQ8bNgx9+vRpNh61Wo1PP/1Uo9PpvF7cyv/bcrabpaWlhuM4oYuLCz799FM8evQIDx8+RHZ2Njw9PSEWi2FpaYnKykrNvXv3uKdPn5KIiAhYW1sjMzMTqampyM3NBSEEQUFBUKlUcHd3x/Tp03HkyBHo9XrweDzo9XqzrjknJwcikQj9+/dH+/btUVNTg6+//hoeHh5gGAZeXl5gWRanT59GmzZtUFRU1CqhAcDNzQ2lpaWtXg8ODkZ1dTXbrVs3AM92hJCQEHz77beYN28ebty4AaPRiNOnT8PGxgaPHz/GG2+8gV27dpFr164hKiqKHjhwIAeAAECWSqX6DsB5AGP5fP7uadOmCdq3b4+6ujrU19dDJpOhqqoKe/fuRYcOHZrp7IVCIRwcHNTFxcX+AP70uf26xK5sbGzkTOcSwzAICgpCUFCQucFPP/2EgoICnUqlStdqtVEVFRVbf/vtt1A+n+/o6+vLeXp6mrdvNzc3REZGwsfnmZ2+Y8eOaNu2rdmrxKSZmjRpUpNB2Nvbo1evXggNDUVYWBgSEhIwc+ZMXLlyBQaDAXq9HnK5vFXZu6Ki4nfl8pb4CZFIhDfffBN79uwxGI1GtkuXLhAKhXj48CGysrLg7OxMAZDnhBK5u7vb5+XlcQzDZAIwchxXxHGcdWBgIGfiZaysrMwfpY2NDTp37owtW7ZgzJgx8Pf3b/J8pVIpwnNr4Z/FaxGbUponFouz0tLSOnXp0qXZ9bq6Ojx69AgMwxzRarVzLCwsttva2g5ydnamaWlpVVlZWRzDMHZRUVGUEEIaGhqaeJiMGDEC3bp1Q05ODmJiYmi7du3IyZMncf/+fchkMri7u6N79+6or6/H/fv34e/vD0tLSyxcuBCWlpbQ6/XYtWsXOnTogMTExBZlXo1Gg9u3b9MRI0aQZhefvSNSUlIwYsSIZtcCAgKg0+nYFStWNHGYqKiowPfff08CAwNRVVUFX19fYmNjYxEQEIDa2tq+t27d6j9+/Hhy/fp1hIWFtTq/3bp1Q0ZGBg4fPow+ffogMjISAFBQUACFQsEBuNc6dVrHa8vZKpVqw4EDB7RXrlxp8nt1dTU2b95sJISs02g0UwCItVrt2Pv370sTEhJk3bt3t1Cr1R8qlUrHc+fO/SMmJqbUz89P/6JYRgiBVCpFbGyssrGx8dKJEyeMTk5O6NWrF8LCwlBXV4e1a9dix44dGDBgAMaMGYOuXbuaJ97X1xcsy6K4uBgXL17EuXPnYDT+27DV0NCArVu3glJacunSJWVDQ/P4gqtXrwJ4xvm3BEJIM4OFg4MDJk+ejOLiYjQ0NEAgEEAoFKKkpATJyck8e3t7UlNTA71e34SDfxkCgQAsy2LhwoW4evUqiouL0djYiMOHD4MQ0kApbd1j5HfwWqIXIWSYWCze4+7uTq5cuYLLly9r27dvb6iurkZBQQFDCFmj1+vXEkI4juNWGo1GrqamBhzHoaKiQsqy7CSDwfALgFNisXjW9evXHZ88eaIPDw/nyWQyPH36VB8fH6/T6/XbOI4b98YbbyAiIsI8lk6dOqFfv37YtGkT+Hw+cnNzYWlpaXb+2759O3g8HhiGwdy5c7Fnzx5cv34dAQEBkMvlyM3NRdu2beHg4OBYXl6eunbt2g69e/cW+Pr6QqFQ4MaNGygpKYG9vT3kcjksLCyazEVGRgacnZ1bJJi/vz+MRiPmzZvXRGSsq6vDli1bEB8fD0dHR6Snp0MgEMDW1raZW9fjx4/h5uYGKysr9OzZEz///DOqq6tBCNFrtdpm3OR/VfSSSqXfv/vuu+Sbb77B1atXMXLkSF5qamoigB4AGErpco7jfAkhPm3atOni5eWFdu3aQSQSwdHREQEBARHZ2dmXhUIh3nvvvY5Tp05loqOj9cePH39KCKk1GAypGo1mIwBfNze3dyIiIprtQA4ODhg2bBh++eUX2Nvbo6amBs7OznB1dUVhYSFkMhl++eUX9O7dGyEhIVi3bh1KS0tRWVkJmUyG4OBg+Pr68tavX99Fq9VqS0tLUVhYCIVCAZlMhsTEROzcuRM7d+7Ehx9+aH6uUqnEyZMnMWjQoBbniWEYWFpaNlPxWllZYfz48dixYwcaGxtx48YNZGZmglKKOXPmmEVLuVyOq1evYvLkyQCAwMBAXLt2DTKZTOPp6anLzs5eRwi5SinVtkSb38Prntl601ctFApBKWV69OjRs2/fvnxbW1vU19fzfv3118klJSWMyZOkV69eZrnXaDSKNmzY0Km6utowZcoU1s/PDwsXLuSWLFlyt6amxvzlikSiZd27d2/ZewFAly5dcPToUbz99ttgWRZZWVk4evQo3NzcYGlpiZs3b6Jnz55ISEiAg4MDZs2ahS+//BIKhQL79+8HwzC4f/8+KSwsFJom1+T9EhwcjHnz5uHnn39Geno6bG1tkZ2dTS9fvqxlGIaEhIQ0Vxzg2cdQXV2NltS8iYmJCA4ORnFxMY4fP45hw4Zh7ty5uHLlCoYOHYoHDx4gJiYGoaGhZmZVo9GAz+fDyclJcPfuXUG7du18MjMzewOI/bN0ey1iKxSKWd999118TEwM8+TJE7Rt29Y4evRo86FraWkJg8HADBw40MzNymQy89fLsiwGDhwoPXr0aMXUqVOVkydPFq9evVpZV1d3hhAym2GYECsrK2+GYbxaUsYAz3zRLl26BB6Phx9++AFarRZ2dnawtrZG//794ezsjE2bNmHZsmXw8fExfxB9+/bFqVOncP78efj6+iIxMZGMHDnS3K+Pjw9WrVoFANi3bx9kMplx3759dXw+X0kpTVar1d/x+fxfc3JybF7Ui5sQFxcHW1vbZg6MDQ0NePjwIT799FNs2rQJDg4OAABnZ2fzMePj44ORI0c2kWru3r2LoKAgpKen4+7du6ipqWHwTJP2p/G6KzuBEOKSkpLyFp/P/9ekSZOaHV41NTW/62ni7OwMo9FIb9++/VtKSkobnU53RSwWrw0MDLQvKysTfP3111izZo3uwYMHhrCwsCbyj1arxbZt2yCVSrF48WI4OjrCaDQiIyMDP//8M5ydnWFvb4/Fixeb5XUTTAaH2bNnQ6FQQK/XIzDw3+5cJvXljz/+CDs7O7Rr144xGAyShoYGvlqtXk0pvUcImbRjx47fRo0aJerSpQs4jkNjYyOuXr2KGzduYPjw4c3et6CgAF5eXhCLxejTpw+GDBmCzp07IzExEQAwZ86clzVmePjwIR4+fIj33nsPaWlpiIqKQkNDgwCAzR8SqQW8tiGEUlpOCMl0d3dXcxzXjNhSqRQ1NTUtOvoBMGmLHNu0afPWc1egbhqNhnF1dSVTp07FhAkTEBQUxA8ODkZBQUET0SwhIQFCoRDTpk0zn40Mw6B9+/aws7NDTU0N7O3tn73gS8xPeXk5xGIxPvzwQ6jVaqxatQpqtboJZ921a1cz0ffu3YvVq1cLtFot9+mnn64FMJRSeoEQsvns2bMfHjt2jBGLxdBoNGjTpg18fX1x7tw5ODk5NSEeIQRa7bNjViaTQSKRYPjw4Rg5ciSWL1+OPXv2oE2bNggODgalFGlpacjLy8OsWbOQnJyMcePGYdeuXYiNjeW9+eabuwC0vpJawZ8hducWflM9dwluhpCQECQkJKBDhw7NoiOAZ9tdZGQkBg8eDACglLIZGRnYv38/bt26hYMHD0IgEBh5PN7jH374wXPQoEHCsLAwRigU4tq1a2Z/sJfRrVs3XLt2rZkyAnhmovz1118hFovx3XffYc6cOWjTpg0SExPNsuyLqK6uxuPHjzFu3Djs3buXGAyGfkKhsIIQclkgEHQ1GAzMu+++CxsbG8hkMjN3npWVhZ9++gkffvgh7O3toVQqkZeXh/z8fDQ2NqKyshKRkZEYOXIk3NzcsGjRIowfPx5Hjx5FY2MjampqEBAQgOXLl0MgECAmJgaBgYGglMLBwQGU0pcXV0u0aYbX0o2/8JuAz+dXLl68WObo6NikvU6nwzfffAN/f38MHz7cvMIUCgUuXbqEhw8fYvHixc3El82bN6NTp06YNm0aFixYYMjLy1tkNBofCQSC05RSgcFggMFgwIYNG1oktkajwbp169ClSxcMHTrU/FytVot169bh4MGDGDRoEKZOnYpz585Bp9NBpVLR6OhodOvWjbzo+rxz506oVCq0bdsWubm5cHBwwMCBA1FWVma4dOkSExERQYYNG9bifB04cACpqamQyWRQKpVo164dysvLYW1tDXt7eyQlJUEkEsHa2tqsKVy+fDnmzp2LuLg4zJw5EwsWLEBycjKOHz8OlmUhFAqpTqdTqdXqf6rV6i2/R5uW8Jfs2ZRSDZ/P37t79+5/LFiwAGKx2HzNaDSC4zhkZmYiLS0NAQEBKC0tRUlJCQgh6NGjR4tyam1tLdavXw8vLy8UFhayS5YsCVWpVEHdu3dno6OjQSnFRx99hOLiYhQUFMDa2hqBgYHm3UMgEKBbt26Ii4vD7du3ERAQAEopMjIywDAMXF1dQQiBp6enWS0ZGRmpO3XqVMnp06c9nZyciFKphEKhgJWVFbRaLdzd3TFmzBg8fPgQhw8fxrvvvstevXoVwcHBrc5NaGgoSkpK0KVLF7Rt2xaenp4oLy/H119/rc/Ozubl5OTA1tYWHTt2RElJiUqv16uePHli09jYiDNnzkAul2Pv3r3Izs7Gzz//jFGjRmH69OnGw4cPH3uR0H8Gf4nY5JnWJMrGxgZr165FaGgonJycUFlZidu3b5sd6SdPnmw2buTl5UEkEqFDhw7NzmIA8PT0xIoVKzB16lR8//33SrlcfocQ8kXfvn15z5+JgIAAbN++HcOHD8fVq1eRn5+PoUOHAvi3mtPDwwOjR4/G06dPQQjBkCFDkJmZiYEDB6Jr166IjY3Fd999Bz8/P/j6+rJJSUke/v7+pGvXrpBIJPDy8gLDMCgpKcG2bdvg6+uL3r17Q6lUYt++fWAYptUYN+AZr1BXVwehUIg9e/Zg/PjxsLKygtFobBQIBFY2NjaE4zhYW1sjJyeHANi8ZcuWjzdt2sQ5ODiA4zg8ePAAQqEQvr6+pvdm+S97gPwZer2miZMAGA1giI2NzcRly5aJa2pqkJycjLq6OshkMoSGhkIsFuNf//oXevXqBQsLC9y7d888+aGhoejatWsTThh4trK//vprCAQCo0KhgNFo1BNCmK+++so8s7/88gu6deuGb775xhz3tW7dOhiNRpw8eRKZmZlwdHTEjBkzmr1HQUEBqqqqcOrUKTg5ORmlUimTmZkJjuOwfPnyFo+G27dv4/r161Cr1WY1aVVVFaRSKaZMmQJ3d/dm95w4cQIymQynTp3CwYMHsX79egiFQl1ycjK8vLz4PB4PdnZ2KCwsRHR0NPbt22fk8/nM3LlzzUYRuVyOXbt2QaFQIDw8HDExMdBoNL/qdLrRrdHm9/BaK1ssFq9ydnb+wM3NTVRYWMgsW7YMhBAMHz4cL59h9vb2sLOzw+PHj1FbWwtHR0fY2NhALpebFQcvIjU1FXw+H7/99hvz7rvvYuHChdz777/fRFFhb2+PzMxMGAwGZGZmgsfj4eTJk0hNTYW9vT2mT5+O77//HkqlssnR8nzsuHHjBlQqFfR6va6kpESg1WrRs2fPVl2QO3fujCNHjmDmzJnmI8NgMODu3bvYtm0b5s6d20TqqKioQFJSklnzlp2dDY1Gg4yMDGbw4MGkf//+ePz4MbRaLYYPHw6BQIAFCxYwX375ZZPdQiqVYv78+diwYQNiY2MREhKCmzdvDiGECCilmj9Lt9ciNsuys0+ePCk5cOAAtm/fjvz8fGg0GnTu3BkBAQGQSCTgOA4sy6K6uhqnT59G37590adPHyiVSty4cQNVVVV48OABgoODQQiBXq9HUlISzp07B4Zh0KlTJ4wfPx6FhYVgWdaYkJBAR4wYwQLP9M87duyATCYDIQRisRjZ2dmYPn26+Vjo1KkTfvnlF0yePNms2Kmrq8OmTZvQp08fzJs3DzweTwAA+/fv/11TJ8dx4DgOHh4eZt7AFDxgMBhw4MABzJ49GwaDAWlpaYiLi0O/fv2wY8cO7Ny5ExqNBo2NjVqGYUivXr34DMM029FM7lv37t1r4rjAMAyGDh2K06dPo7GxERzHQa/XO+B5UOCfwWuJXgzDlMbFxdn7+fnB2dkZTk5OoJTCysoKX331lXlCnJ2dwePx8NFHHzWZzKCgINy8eRNHjhzBhQsXYGlpidLSUjg5OaFjx47IyMjA7t27sXXrVohEIhWPx7t+/fr1Hs7OzlIej0eOHTuGrl27onPnZ0N6+PAh4uPjUVtbayb2qFGjsGvXLqxevRp9+/aFra0tYmNj0alTJ/Tv37/Jizk7OyMzM7NF7xAA5kCFl3cJ4BkjduLECaxduxYsyyIwMBCzZs2Ch4cH+vXrh/j4eCQnJ5s8WvVisbjVM9fR0RH19fXNfndycoJCoQClFBqNhgeg8aUmnZvd1AJeKw9afX39jI8//jjZtE1GR0ebnf4++ugj2NnZob6+HteuXUNlZSUUCkWzldOzZ0/cunULnTp1gouLC+zs7GBnZ4eqqiqkpaVh5cqVN+vr61MAXAJwGkCno0eP/koI8Vq4cGETv7U2bdqYOWtTX/n5+aitrYW7uztKS0uRmZmJoqIivP32203GkZWVhYyMDFRUVODu3bsICQlpcp1SanZ7ammb5/F4cHR0xPr167Fhwwbk5+ebP7iCggLEx8ebAh+4W7duNVPgvIiqqqom4cAmPLd4mXaYW0qlsq412vweXpcbz2UYxiCRSFiO4xATEwNfX98mcrOlpSVGjBgBGxsbHDlyxBws8CJMDoP9+vUz/2ZnZweDwYD6+vrIFy07AFI4jjvXu3fvd11dXZuN283NDWFhYdiwYYN5a2cYhmZlZRns7OwU5eXlQpFIJNixYwfat2+PYcOGoba2FgcPHsS2bdsgl8sxf/58KJVKdOvWDRzHobS0FGfPnkVeXh7eeuutFifCYDCgsbERXbp0wfjx4/HJJ59Ao9FALpdj3759kEgkoJTC2dkZUqkUSUlJ6Nu3b7N+5HI50tPT8cYbbzS7FhsbC7VajZSUFK1Op1vyR8RpDa9FbAsLi/UjRowwfvPNN+zChQtx8uRJTJ06tUVRpEePHrh06RJKS0ubfbUvBronJSUhNTUVAoEADMOoDAaDnhDSFkA0AAmAh0KhcFhwcHCrYw4ODkZKSgpEIhEaGxvBsiyhlPJKS0tjJRKJ/4cffthu3LhxePvtt5GUlARLS0sEBweb46jWrl2LmJgYnDp1CizLwmAwwM3NDVqtFnK5vEW9wP3790EIwcWLF7F27VpzcOG2bdvQt2/fJoT18fHBpk2bzCZW005hcto0uSmboNPpcPHiRRQXF5v4Gk4mk+0jhIRQSpvv93+A1yI2n88PHDNmDN/BwQGBgYG4e/duqzIny7Jo06ZNi8S+d+8eKisrzf7Uu3fvRmxsLLKzs4U6nU5HCDGGhIQYpVIp9/jx48aSkhJJXl5ei6IO8MyEyOPxEB0dDV9fXzAMg6KiIvz222/Di4qKyDvvvANXV1eMHz8e58+fR2BgII4ePYrdu3ejoaHBHKxQUVEBhUIBR0dHzJs3D7Gxsdi8eTPmzp1r1rlTSvHkyRMcPXoUHTp0wP79++Hg4ICnT59ixYoVsLa2bsYDODs7Y86cOdi9ezdOnDih9/LyUjU0NNCSkhIxx3GstbU1+fzzz+Hv7w+GYfD48WN4eHhg8eLF2LRpExQKBUQikY9arf4ewLQ/S7fXkrM5jpvt6Oi44YMPPhCvX78ePB4PixcvbvXeLVu2IDw8HB06dDD/9uDBAxw+fBgLFizA+fPnMXDgQHz55ZdmP+0OHTpgwoQJTfzNc3JysHPnTsydO7eZMqakpASbN2/G0qVLm/EHz9WrBisrK+P8+fP5y5cvR3R0NNq3b4+nT5/i6tWrKCsrg0QiQbdu3eDo6AhPT0/s3bsXjo6OGDp0KJYsWQIejwc3NzdzNKpKpTIZbGA0GvHZZ58hJiYGwcHBCA8PR2BgYItBAFVVVfjiiy8UBoNhDgAdwzC7GYYRrVq1ChqNBtnZ2aCUok2bNuaPKzY2FpaWlujbty/ef/99dUNDgzl4/7+aLUmr1W4vKiqa9sknnxyvqKgwlJeXo6qqqsW2dXV1yM3NRUNDA2pqalBYWIhjx47h8OHDmDVrFhwcHNCjRw/s27cPv/76KyZPngyO4zBp0qRmgQU+Pj4YNmwYzp492+w5ly9fRs+ePVsUoViWxdChQ9mcnBzj4sWLDTqdDkePHsWePXtw4cIFPHnyBAqFAr1794abmxseP36MO3fuYNiwYbhx4waePn0KqVRq1mXn5uZi1qxZ6NmzJx4/fgzg2QelVqvRuXNniEQidOzYEXK5vMU5eW7rFgG4A6DUwcFBazQazf2bvGVNhAaeydxarRYjR46ESqXizIlX/gT+cgpqoVBY7ubm5lBfX49u3bpBJBIhMDAQtra2UKvV2Lx5MxoaGiASicxiRa9evdCrV68mftF3797F3bt3kZubi/Dw8GbKGRPUajU++eQT9O3b1xgSEsIAzxQxCQkJmDJlCgICAlq8T6FQYOXKlYiKikJERARKS0vNIpWvry8qKyuxdetWcByHOXPm4MqVK+a0GRzHwcfHBwqFAiUlJfjuu+8wc+ZMZGZmol+/fvj4448BPPNy0el0aN++PY4fP44PP/wQNjbNTc/379/HxYsXUV5ermIYpp4Q4mg0Gsl7773XbMcyYf/+/bC0tIRarUZ6ejrVarVHNRrNW5RSQ2u0eRl/KQU1IUQI4GxJScl0kUiE6upq6PV6nD17FlZWVmhoaIBQKMS4ceNgZ2eHuLg4cByHqKioZp2HhIQgJCQEn3766e8qOIRCIViW1V6/fv23hISEEQBEJqVJaysJeJblyGSAAZ6dny/yEO7u7hCJRPjHP/6BlStXoqKiwmwA6dq1q5mZunDhArZt24bg4GD8+OOPTfoYO3YsUlJScO/ePWi12hZNu6b5qaurQ3BwsKhHjx4ikUiEs2fP4uzZs5g9e3YzEa+qqgoPHz7E9OnTodVqkZ+fTxiGGV5WVrYSz4LyU1t98Rfw2q7EhBCxRCJJFQgE00eOHIlPP/0UEyZMwOTJk/H555/D1dUVtra20Gq1urNnz5qsWE1celuCVqulOTk5rW43lZWVMBqNBq1WO5lhGNWVK1fQ2NiI9u3b0/j4+Fb7NZkUW3NzAgBra2vExsZCLpfj8uXLkMlkCAsLazL5AwYMQHl5OYYOHYq7d+9i1KhR5msMwyAkJARGoxEeHh7YsGED7t27B71eD0opsrOzsWXLFlhaWoJSiuHDh8PLywuOjo6YPHkytFotdu3ahbKyZ8EeJo3c999/b45lO3v2rEnaEBuNxvefL7hXwl+xeo2wsrLy9vf3R/fu3ZtcMJ25a9asgUqlqlAqlUs//PDDz7VarVddXR1pKcAegOljaHz48KGkqqqKfTFk1oTLly8DwF1KqV4oFAqdnJxACEFYWJhxx44dypiYGMngwYMZtVoNpVIJqVRqPoMZhoHRaGxVBy4SiVBeXg4rKysIBAJMnTq1WRuWZTFx4kRs2rQJ4eHhzXahW7duoaSkxJxtITExEQcOHADLsrC2tkZ4eDh69uyJw4cPIykpyRzIz3Ec5s6di8OHD2PDhg3g8XjQarWwtrbGhAkTQCnF7du3MWLECAQGBuLixYuIj48XyeXy5hxgK3gtYhNCbAkhwysqKvhjxoxpsQ3DMIiIiMDZs2edtFptrMFgWDJt2jQSGxuLEydOoF27dhAKheY4LY1GgyNHjhgACAMDA8kPP/yAkSNHQq1Wo7a2FoQQ1NbWIj09HUajsYHH453k8/lyPz8/MSEEfD4fOp2Of/36dRIfHw+dTgeRSASVSgU+n48JEybg8uXLePjwYROpwITy8nJkZmZi2LBhKC0thUajadHbBXgW8xUcHIzNmzejY8eOCAwMhE6nQ1JSEi0rKyMLFy7Ed999hylTpkAikUCn05nt+6atPSgoCLdu3WrSL8dxmDx5MhYvXgyRSARKKT7++GOzWPuiI2Lfvn1x5coVFkDLfl8t4E8TmxAilEqlydHR0S4cx+HQoUP48MMPW1QBWltbQyAQsDqdLl+r1bLXr183Zz/Kz8+HSqWCVquFj48PcnNzIZfLmYEDB7I1NTUQCoU4ePAgvL290aZNG1RWViI9PR08Hs/g5+c3rLS0FMOGDTMT7v79++yZM2dYvV6PwYMHo3v37hAKhWan/0OHDsHb2xuHDx+GVCptEvxfVFSEbdu2oX379oiLi4NOpzPyeDxm7969aN++fYviE8uy6N27NyQSCdLT01FRUYHS0tIMGxsbd0dHRxkhxHxktWSCfp5wAElJSbCwsGjigMEwDHr27Im8vLxW9Rccx8HV1RVPnz59ZRq+TkRIRwcHB9t9+/ZxwDNDQH5+fouroKioCDqdDn379mUDAgKwd+9ec44Q0wTk5eVhz549pvxhJC8vDxYWFlAqlWbPUROUSiV+/PFHtqCgAB9//HGTOKuwsDDExsZi6NChZgMJ8EzMGTx4MIRCIVJTU8EwDLZs2QIHBwelp6cn8zw6VPjWW29h7969yMvLQ1BQEBMREYGJEydi2bJl4PF4TfqUy+W4e/cu/vnPf8La2hr5+fnYunUrZRjmulwuH5mRkSEzJQzo3bt3i/N5+/ZtlJaWokOHDrh27RpycnIwfPhwPHr0CK6urnB2djaLda1Bq9XqAVS9akTIKzNo5N9lhkrLysr46enpyMrKwpMnT5pMuglKpRIJCQno3r07FAoFdu3aBQ8PD0RERDT50r28vDB//nwYDAZMmjQJ77zzDvLy8jB58mS87Ndmymfysnsw8OyjMRgM6NSpU4vj79WrF2pra9G3b1/o9frTxcXFE27evPlxVlbWOwzDZJqsSmq1GpRSbN26FZMnT8ZHH32EJ0+emPspKSnBli1bzKJjTk4Otm/fjkWLFpGJEydOp5Ti4MGDyMnJwfnz51vUPzx69AhZWVmYMGECfvnlFyQmJiI+Ph4ml6SIiAi0bdsWhYWFaCkODXjm5PGckbtOWikP9TL+9DZOKS3kOG52t27d9nAcRywtLREXF4dhw4aZY6KKiopw9OhRhISE4M033wTwLMB869atuHnzZrNMwba2tmjXrh3kcjnKysqg1Wqb+VCbIJPJ4O/vj8uXL6OmpgYVFRUmcQwmZq3FF32u/bp165aWUrqdUmrWzBBCrsTExBT5+voyJtl7w4YNePvtt7Fz505kZmaioqICSqUStbW10Gq1EAqF+Pbbb1FWVobIyEisWbMGxcXF3LFjx5z69u2LTZs24e2338b69evRq1cvBAUFQavV4s6dO6r09HSjwWBgMzMzhXq9Hk+ePAGfz8e6devQvXt3BAcHo6amBizL4uDBg5g2bRoKCgpQW1sLmUwGLy8vHDhwwPSur7xgX4tB02q1+6ysrFYePny4TY8ePdC3b1+sWrUK9vb2UKlUoJQiMjKySTCelZUVJk+ejF27dqF79+7N4p7d3NxQW1sLtVoNiUTSKtGMRiOqq6tRXFyMfv36YdCgQWhsbMT169eRnZ2NysrKJpqnF/Gc25dbWFjMIIRkU0pN++QwBwcHZXV1tZRSCoZhcPDgQezduxcAMGbMGFhYWJgdGE6cOIGGhgZzJqPjx4/j0aNH+Oyzz2BnZ4cpU6bAy8sL69evx6hRo8yVB6qqqowqleqA0Wg8IxQKF9y+fTtSJpMRo9GIzp07IzU1FeHh4TAajcjOzka7du1QX1+Pzz77DDY2NnB2dkZlZSXKy8vh7u4OOzs7RWlp6X+XGweA+vr6eaNGjTru4uKiLy0tVev1+sXl5eULXVxcOn/wwQdsS1mA3N3dIRaLUVBQ0Cwfd1VVFTQaDeLi4szccEtWpoSEBADAxx9/3OS6Kb31jz/+iCVLljQTrzIyMqDVavH111/b1NfXR69evTqSEOJBKZULhcJIvV4vNRgMeOedd8zJ6E1eqQcPHsSUKVPMY46IiMDXX38NAOaUm507dwallEZERJAlS5agqqoKu3btQqdOncy+8XFxceT8+fPhHMdNHDx4sNjPz4+o1WokJycjOTkZLMti9erV4DjOnPynpqYG7777bhMXrqKiIvz000+glHIAXllt+leyEseoVKo2OTk5/ZVKpTel9Gej0bjFwcFB3RKhTZBKpVCrm9ZTUygUuH//PgQCAWbOnAlra2scPXq02b1GoxGxsbEYO3Zsix9C7969YTQakZ6e3uT3uro6HDx4EF26dMG8efOwdOlS4uDgYAEgytraOobH4w0oLy9vQmjgmSdru3btMGHCBJw6dcpsjrWzs4NGo4G3tzemTJmCyMhIU/4XkpycDF9fXxw+fBje3t5NEgHU19cTiUTiu3TpUknv3r2Jg4MDPDw8MGrUKLz99tswGo0oKyszG2bKysowceLEZr56bm5umDFjBuRyuQhA05f9HfxVv/EKABUv/JT59OlTSiltcRvWarUoKioyB7UBzwhhytW5e/duc2jvokWLcP78efTp08csqyYkJEClUjVLlWUCIQQhISE4cOAAHj9+bM6x9vDhQwgEAty/fx9nzpxBQ0MDKioqWKFQuGvZsmWiL7/8klhZWbWql27Xrh1+++03FBcXw83NDSbRsE+fPiCEIDY2Fl9++SXee+89HD58GCtWrICrqysuXryIkydPwtbWFj169EB6ejpGjhzJtOTeFBAQgLZt2+LYsWNmwkskErOuv6CgANnZ2XB2dkZgYCDc3Nxgb2+vKykp6faq9PpLKahbQKJKpapMS0uTviiqmGBKOLNt2zZ4e3ujsbERWVlZmDZtGoqKivDzzz9jzpw5OHz4MDQaDQoKCrB69WrIZDLI5XKzT9vvgRCCbt26wdbWFo2NjfDw8EBoaCgOHjyIsWPHYtGiReYcMLW1tZxUKiVisRhOTk6t9skwDOzt7c2ccXx8PEJDQ5t80CaZ2hTm4+zsjPnz58PCwgJFRUVmhvJFxcjL6NixI9avX2/csmWLms/nw8LCQkwIwePHj3HkyBFMmjQJp06dQlVVFcLDw+Ho6KgtKSlx/D9LQf0inqejnnDgwIEr1dXVkpCQEGJlZQW5XI74+HjEx8fDYDAgOjoa9fX14DgOWVlZWLduHdRqNcLDw7Fx40bzdtmrVy+8/fbbaGhogFgshkQiwb/+9S/k5+e3uLoppbh37x6USiUYhoGLiwsEAgEOHDiAN998Ex06dDArYdLS0nD48GHeN998g6CgIHOpB4VCgTt37pgtYh07doSPjw8qKyshlUoRFxeHlJQUfPDBv2uv9e3bFytWrMDu3bvx6NEjjB07FqGhoebrPj4+8PLywuLFi2E0GlvMaggA9fX1KCoqqmUYhgAwlJeXCwoKCljT8wwGA0aPHo3jx48jPDzcVBKj6P8sBfXLoJTe5vP5qy9cuPBFTEwMEQqFUCqVEIlE8PX1hVqtRrt27fDgwQNkZGSYUzu3b98eDQ0NcHNzo0+fPs2mlGoTEhLaBwUFQST6d5L98PBwnDx5EnPnzm1m1EhMTIRKpcJ7770HkUiE7OxsnD9/3qzefBHV1dVUIBAYy8rK2LFjx+LmzZu4ePEi4uLi0K5dOzg5OeHJkyfm/Okqlcq4ZcsWDaVUuGjRIvJiyi03NzcsXrwYN2/eRH5+fqsaNy8vL6SnpzdzagT+nZs9PDzcKjg4mGVZFvfv38eOHTtgb2+P9evXY9asWXj06BFUKhVycnJQW1srBHAFwLJXoc1/vFg6IUTIMMyqu3fvMvX19XBzczMHuPfu3RsGgwEZGRmm2l3w8vKCQCCAhYUF3n//fXAcR0aNGtWW47iA0tJSxMTE4Ny5c9i3bx+OHj1q3sq/+OILJCUloaioCJmZmdi9ezfOnz+Pf/zjH3B2doaVlRW6du2Kjz76CHK5HHfu3DGPUaPR4OrVq2T06NGsk5MTKioq0K1bN8TGxmL+/Pl48803kZiYiCFDhmD48OGglEKn01Wp1equAFp09zUZRIKCgloVGwcMGICTJ082U5TU1tbi8OHDmD9/PoYPH866ubnB2dkZgwYNwgcffICSkhIMGzYM33zzjTkO7McffwSepaJuWYvUAv7jKxuAiBDCtG3bFkKhEH5+fuA4zpyUrqysDOnp6Vi0aBHmzJljTh05e/ZsyOVy5OfnIy4ujsTExMDd3R2xsbHo2rUr2rVrh4aGBhw+fBhKpRKhoaE4efKkXqfT8QBAKpXqlyxZwns544FAIMDQoUMRExODsLAwlJeX49ChQ6CUGoOCghiVSoX4+HhwHIc333wTzs7OSEtLQ3BwML777jtQSuHp6QmtVivV6XROWq3WePr0aXbevHlN7AG3bt0y+6W3hnbt2sFoNOLrr79G79694e/vD7Vajd9++80QHBzMthTLbgr+u337NioqKvDw4UNQShEVFYVTp07xALR7VcL8N7bxWktLy7Pdu3cf4O3tLYyLi+OZ/NMEAgF69eplrnw3fPhwHD161Jx598CBA+jduzdiYmKgUCiMhYWF7LJly5pkKezbty9++eUX3L9/nwoEAl5ubi70ej3at2/PUyqVLdbncnR0RGFhIdasWQOj0YiOHTuioaHBSAhhunTpguvXryM3N9dcQcDa2hrnz59HdnY2qqqqUFdXhx49eoji4uLG8fn8bysrKxd99dVXJCgoCHK5HNXV1airq8PYsWNx4sSJFtW5AJCfn282nd68eRMHDhyAXC5XUUrrg4ODW+UQBw4ciE2bNsHd3R0SiQTTpk1DQEAAkpOTaV5eXvNcH63gv7Gy0dDQMCYtLW3MgwcPPh06dGj7F5PMDhkyBDU1NcjKykKvXr1gaWkJb29vunr1aoVKpXpSWloqvHnzpprH4zlNnDjR+eV0lAzDIDQ0FMnJycTb29ssxtnZ2UGhUDTTnpnO3d69e+PevXsYOHAgrl+/rmxsbCQVFRW8J0+eoLKy0uyEDwAeHh7o0aMHTCk0xo4dC7lcThiGkSkUim95PN77MpmMl5KSgi+//BLXr183p7wyadhM279AIAAhBEqlEseOHUNkZCQ8PDzg4eGB+/fv49ChQymEECmltHVxAM8sZ6tXr27ym7W1NcnLy2s9H+dL+E+LXqa2BgBHGIbxzs3NbVKphmVZTJ48GQUFBbh69SoeP35MS0pKDmo0mm0AEjQaDSWEeAqFwoyW9OPp6ek4evQooqKicPfuXbPNu6qqqsWUHjdu3MAPP/yAyZMn4/Dhw5gzZ46yvr5+FMdx4bt3716kVCqFixYtwtatW5uEEEdERDRR9+7atUup0Wgox3EpHMfxampqMG/ePMydOxezZ88Gx3HmhL0FBQW4efMmGIaBTCaDg4MDiouLERoa2sQukJycrFSr1ccZhrFJSUnx8/Pza9HrJD09vZmtgFKKwsJCDYAn/02r16u295RIJO3z8vIMaWlpL1+Du7s7NM/wuVqtnkwpvUH/7f0o5DhO/7LKU61W4/Dhw5g9ezYiIiJgOjdNeceSk5ObjYNSipqaGgDPwmie/UQvaLXalbW1tfKhQ4eac6CaMhm/jJKSEmRkZIj4fH70hx9+KOrSpYsp0hQpKSnYtm0b7OzsMGDAAGi1Wrz55pv46quvsH79ekyZMgUqlQqOjo6IiooyM28ZGRl49OiR3mg07tHr9VuTk5ONubm5zZ5dXl6Oa9eumRPom/D48WM0NjbWAIh7Vdq8bnx2F6lUekKr1ToIBIJDjY2Ns15MsUgIEYvF4qcLFiyws7a2ZlauXElMdmyO41BUVIRz586pcnJyMjQaTW9KqeqFexkA/Xg83umlS5cKX8wnlpCQgKysLEyfPr3Z+MrKysyxXj179oRIJEJOTg6uXbsGQghcXFxQWVkJtVr9UK1WB5FnBdx169atYwQCAZRKJb7//ntYWFggKioKHh4e0Gq1SE5ORkxMDPR6vSEiIgJDhgxhHz9+jJ9++gkDBgwwF2KNjIzE3r17sWjRomY50PR6Pb799lsUFxebUoTpKyoqVHq9fujzzFMEwBqhULjEYDAwfD4fbdq0gVgs1qelpfE6duyIcePGgeM4GAwGpKam4siRI0qtVjuSUnr5v5pvXCqVFm/fvt1l0KBBGDJkiPLevXvzKaU/vdC2k4eHx/X8/HwZAPTp00eRlJSUazAYAhmGMTAMIzcYDBsNBsNXLxGaWFhYHLW2th6iUCg4Nzc3/tSpU80u0seOHYODg0Or0ZbP83HDzc3N7LQnFosxffp0EEJw9uxZZUZGxlKj0biRPCvCrtuwYQPRarXYunUrxGIxbGxskJGRgcbGRhgMBnAch9GjR+PMmTOYO3cunJyccP78eajV6iZJ5c+dOwe1Wt3EAfFFPHr0CEePHsW7776LK1eu6O/fv1+p1Wq7ACgXCAR7ZTLZqKioKImfnx80Gg3u3LmDy5cv67Va7TqBQDCGz+cHqFQqiEQiaLXaYq1WO4FSeuNl2vweXsctibAs6zhixAjIZDIMHDhQ/PDhw08IIY8opUnPmxVUVFQwJ06cgKWlJVJTU4lOpxsDIN9oNAoBNFBKW3Iz9WEY5o3MzEyRVquFm5sb3blzJwYPHgx3d3dQSs1lEluCSTM1Z84cHDlyBC4uLoiMjATLsrh79y598uQJC0BmCmYXiURZmZmZ/klJSbCzs8PkyZNBCAGlFFqtFjweD8nJyYiNjQUhxLwFl5eXmwvB1NfXo6amBkVFRc0cL19E27Zt0dDQACcnJ4wdO5ZnZWVll5CQsFutVu+3tLQctWjRIonJuCMWizFo0CAEBATwNm7c+CGAXQMHDvT5+eef+Zs2baJffvlloUajufHnKPd6zgvU0tLyxLBhw6L69esn2LRpE3Q6nbdAILjJsuy/jEbjakppLSEkatq0adsAcHK5/OMXbMfq3+neRa/XcyqVChKJBDY2NmhsbMSePXtQW1sLlmUhk8kwaNCgFj1EU1NTzTU7XFxckJubC6lUijt37igLCwvFs2bNEjx58uST27dv9wIwTK1Wf3Ho0KGdKpWKt3z5cjMxCSFmq1pYWBiuXLli8jWjAwcOJPn5+ZBKpXj48CGOHj0KDw8PFBQUoGPHjq2+WEFBAXg8HpYseRaEyefz+QaDYYBQKGw7YsQISUtWPA8PD3To0AEZGRl9hwwZwpdIJBg6dChZt25dyxabP8BrBeM3NDTsT0xMHFFZWYkZM2bAx8cHBQUF5MyZM0sLCgpGEkJC6bOSvy2HZ7QAlmUncRy309HRkXh7e4PjOPD5fFRWVuLtt982O+Rt3rwZJ0+exMiRI5sQvKysDL/++iusrKyQlpaGqqoq/ePHj7OfPn2aqVKp9AsWLIj+/vvvWY1GIxKJRIOfb+OXlEolMdXGbgmmhD3JycnaK1euUBcXF4Fer8e9e/eQl5eH48ePo3///hg5ciQSExNbXN15eXnYtWsXBgwYgF69ekEikaC8vBwXL17k3bt3r21rUSwAEBQUJM7MzNR89NFHqpKSEtGhQ4cUOp3up5eadW7p3pfxpyNCCCFSHo+3f/78+U3qZnt6emLu3LnM9u3bg3Jyci4BCG+lHzMIIW4sy840GAyWHMfN+eCDD0TOzs6orq6GTqeDo6MjKSwsxLZt27B48WLY2NhgxowZ2LVrF1auXGlO38zj8ZCXl4eRI0dCKBTizJkzqK+vh16vH6TT6QoJIZOuXLkytLi4WBIfHw+RSFSlUCgMfD5/TmBgoL6srKx1Azye8QIajaYEwPo9e/Zs8Pf3Fzg4OODu3btmf7r27dvj6tWruHz5cpNi6kajEXv27MHEiRObrHxTYEBqamqrShjgWdgugKL6+vrFa9asGWQwGB4AONwSbf4Ir6NUmdi2bVv6coF04JnCY/jw4di0aVNvQsgkSunB1johhNiKxeLUadOmWV68eJF4e3szplCaF7lZDw8PdO7cGUePHsX06dPNuT+PHj2KxYsXo6GhAV988QUWL15svs/X1xdffPEFtFptB0KIFsDh/Pz8Xj4+PtM4jqtUKpVjgGdF6Lp06SI4deoUCgsLWwwF1ul0SEtLMxoMho8BnDcYDPyampqv1Wo1T6vVIioqCn5+fsjIyIBOp0NsbKy5hBPHcbhz5w74fH6LvuqEEPj5+eHu3bvN/PJMiI+PV6hUKh0AZ4PBsJy+KkfdAv40sQUCQbfAwMBWg7Hc3NxMZ95KAK0SG0CfkJAQ/ubNm3ne3t6/m0AuJCQEu3btUixZskQMgDg6OqJTp06YN28egGfZgxUKhZnYEokEbdq0YZ88eXKaEKLj8/kHGhoaZlFK33uxX0qpUqvVmlWwc+bMwc2bN3H16lVzLZD09HQwDENtbW131dfXi1mWJdXV1aRdu3aYN28eamtrUV1djejoaKSnp+PkyZNwc3NDZmYmJBIJLCwsmiTeeRkDBgzArl274Ofn10z7d/r0aVpfXy9ZvHjx6DNnzgwrLCwMBbDod+b0d/GniW0wGOQqlYriBd8njUaD6upq8Pl8mMoDG41GL0KINaW0tpWuCh89esTLy8trEgRXWFiItLQ0aDQaODo6IiQkxHQ2F1NKu4lEooOVlZX9GxoauFOnTkEul6O0tLTJblBdXY3CwkKSn59PxGKxwN/ff3xDQ8MuADdfHIBKpXp88+bNN95//31SXV2NtWvXwsrKCo8ePUJeXh4GDx4MPp+PFStWsCKRSFpfX2/y/cIbb7wBQghsbW3Nz+7Zsydyc3ORm5sLR0dHPHz4EEaj8XcdFry9vUEp1X711Ve0W7duCAwMFKjVaiQlJcmLiorYlStXij766CNmzpw5ko4dO87E/yWx9Xr98cTExPcGDRrEU6vVOHfuHO7duwcLCwtzUjhbW1tUV1cbDQZDq2chpTRZLBav9vf3/4wQwt67d48tKipiSktLERoaCjs7O2RnZ+PcuXNwc3Mz6PV6HY/HKxEKhTpCiF6n07Fz585lLS0t8fbbbzcxgGi1WrPMzLIsbG1tjSUlJR48Hm8wx3ETKKVCvV7fKBKJ/Orr68mNGzcwatQoSCQSsCwLDw8PuLm5geM4jB8/3mxPl8lkqK2txciRI1tdqeHh4cjPz8fMmTOh0Whw7Ngx88fbEsedm5sLo9FYpdPpet66dWvuvXv3wvEsAfABHo/nde3atY8WLlwojIuLozwe77XKMprwOmf2DbVaXXT8+HGvvLw8ODk5mc+nHj16oKamxlS4RAOg+vc6UiqV6wCsI4R0jo+Pv2sqoWiSlyMiIpCXl4etW7eywcHB/sOGDeNZWlqKTCkuDh06hK5duzZxEgRgyj9u7NKli14mk2lzc3OreDzejpCQEH63bt2ELMti8+bNsLa2RllZGc6fP4/k5GR06NAB8fHxWLRoEbKzs2FpadkkX5lSqYRer0dlZSUOHDgACwsLDBgwoIlzhaWlpbmQrEAgwKRJk5CdnY19+/Zh+vTpTRix+vp67N+/X6nT6ZZTSvMBNEmOQwgRJSQk9BKJRJFisbhcoVC0XG3uFfFnnBc6A8/kbI1G0+vWrVtGKysrFBcXo0OHDujYsSO2bt0KV1dXLFiwAJRSCYDmXFzLICzL6sePH9/MZScjIwPt27fHhAkTeKbIExNjs3DhQly8eLFZodOMjAxUVVUp7t+/P/fmzZvz1Gq19fTp02UTJ04Uent7o6amBhYWFli8eDEUCgVcXFzA5/ORl5cHnU6H48eP4+zZs3j77bebiHemMo13797F/Pnz4eTkhEOHDjV5dmHhs1x0pjERQhAVFYWcnBzNp59+qvrpp5+Mly9fxpEjRzRr1qxRNzY2rjMYDLtbmhRKqaqurm4ApZSTy+UulNLWPEk7v8okv1YwPqW0RCAQyDt37myRmpqKH374AcAz3XVhYSH8/f0RGhrK3r59O4UQ0pFS2iQb33P9ZxAAIYB7DMOM6d69O+9lQlNKkZiYiLlz57a4bVpZWaFTp07YtGmTduDAgRwhBHfv3pXn5uYa9Xr9dgBOAIw+Pj68F4u/1tTUgFKK/v37g+M45Obmoq6uDlKpFIsWLUJGRgY0Gk2zUs8Mw8DPzw+jR4/GlClTEB0d3YRXMBgMuHLlCtq0aYOtW7di0aJF4DgOTk5OMBgMjUKhUOri4qK/cuUKq1ar91NKV1JKi/5o4ukfl3ZK/aM+gNd0SyKECHQ6nbRt27YoKipCYmIi7t27hydPnphTZ3h5ecHFxcWS47hrhJAmH5WFhcV2Ozu7JDc3t1gLC4srPB5vulQqbTYWrVYLlUrVYiI4E9q0aQO5XJ5+8uTJX0+cOPFbVlbWB0KhMDUyMvIf8+bN+1wkEq0OCAhoIj1IJBKIRCKMHz8ey5cvh1AohCmo/+bNm5DJZKivr2+yPZvQuXNn7N27F0lJSVi/fr1ZXKuoqMDu3bvNzgnW1tZISUkBAFNyHsmxY8eE58+f5zZs2MBaWVk5vQqh/5N4XecFLSFEz+PxuOjoaERHR8NoNGLo0KFmZ4KGhgZ4enoSnU5nV1pa+gaAkwBACGljYWEx+enTpyKRSAR/f//ueXl5goyMjGbZ/E0KC4VC0aIHCgDU19cbjUbjHbVaPe95/3Y6nW7zxYsXOR6Ph/LycoNpazWhY8eO+O2339CjRw/cu3cPw4YNwzvvvIPa2lqEhYVBKBSCz+fjypUrGD16tHkcphxnpswLJpF37dq1UKlU6NmzJwYNGgSWZdGjRw8kJiYiNDQUly9fVlBKK8rKytoAQHFxscFgMLQcsfdfxOsWhKEikejXxMTEsYMGDWJf9hE3GAy4desWJk2aBAcHB9nZs2cn4Dmx8VxkYxgG169fR3l5uaBjx454/Phxs3SSpoS1iYmJGDBgABobG5GQkABKKXr06AGZTIYbN26odDrd3hce30gpNdy+fRuBgYFISUnRq9Vq+uK7isViREZGIikpCTNnzoSVlRV69eplLlK+e/dujB8/Hnfu3MHnn3+Odu3aQSAQmENox4wZg19++QUffPABjEYjKKWws7NrwnxZWFhAoVBg7969qKmpkQCQzpkzR/7Pf/4TcrlcrlQqX8kj9D+J13ZLUqvVqy5dujTc1dVV8uJ5qNPpcPDgQTQ2NiI9PR12dnZgGObFZZlLKT3i7e091mg0inbu3MlMmDABP/zwA5YvXw69Xo/g4GBwHIfy8nJzPw4ODoiJicGwYcMgFAqxefNmODs7a3Q63W0At02dU0o1LMuO7d+//yGdTicWiUT71Gr1W48ePeK1a/dv37zBgwdDp9Nh7dq18PX1hY2NDXJyclBVVYV+/fqhffv2CAoKQkVFBeLi4nDz5k1YW1tDo9Hg3LlzkEgkKC4uRkvBEMAzf7Pq6moEBARAJpMhOjradtWqVbcqKiqmAcinr5FC+q/ir1T/ySCEDNq7d+91e3t7xt/fHyqVCunp6Wjbti3+8Y9/ID4+HklJSVCr1Wkv3EcJITMaGxs3WVpaxlBKzWojNzc3pKam4tixY+ZUl56enggODsbPP/8Md3d3kwstrl+/jkePHvEppfEvqxANBsNZQoglANLY2GgkhOzevXv31S5dunDdunUjQqHQvJP06dMHQqEQcXFxWpVK9ZvRaNwTFxf3c1JSkq2trS2prq6GWq1GSEgIJkyYYObO7927h6tXr6JDhw7NJAiNRoPY2FiMHz8eMTExkEqlCAsLYwghDpTSrNed87+Kv5wHjRAyj8/n/9CpUyfi6uqKoKCgJmq/Y8eO4c6dOyUajcaTUqp/6d5vJRLJ+4GBgcjOzsasWbNgbW2NL7/8EgMGDEDPnj3NW2N1dTU2bNiA+Ph4CAQChIeHY/Lkydi/f7+yoaFhoV6v3/kH43fh8Xjz+Hz+ZEqpE6WUb2VlpTQYDMaGhgaD0WhcqdfrN3Ec95VEIpkbEREhsbGxgYWFBbRaLc6ePQsnJydMmDABphQaP/30E7RaLUaNGmUOHyooKMCRI0dQU1Nj1ibW1tZCKBRq1Wr1h1qtdvMrTfifwKvmQftL1X+e/z7MwcHh1LJly1rUljU2NmLVqlWglH75ckZdQogbj8d7MmPGDKGHh4epsi6ys7NRW1uLyspKhISEYOzYseDxeLh79y5+++03AMAbb7yBLl26oKioCNu3by/XarU+lNJXrmRHCPHAMxOsu5WV1QyGYfg1NTU/icXibz755BPxywyhVqvFt99+izfeeMOs/jQYDPjmm2+gUCjAMAy0Wi30ej26dOmC6OholJSUwGAwgFKKffv2aRUKxUK9Xr+NEOIOoCeeOfnfoJSWvOq4W3mX/7Pwn/YBAQGtfjEymQz29vYoLy9/jxCygr6QVppSWiQUCo9fu3YteubMmWIAZtXi119/jZEjR2LEiBG4efMm+vTpg5CQEFhYWODIkSO4evUqjh8/DqPRCD6f70gIaRCJRLlqtXotgD1/ZB2ilBYQQhpFItGxb775RvbcfSk4MjKSmAhdU1OD+Ph4pKSkQK1Ww8LCAmfPnkVAQAB4PB6USiXq6+uxe/duyGQyjB07FoMGDQIhBGvWrIFYLDbb5P38/LiHDx+uFwqFU/l8flcfHx8dAH1OTg5fJBLFqNXqaZTSekKIC4DuePYhJDyPlP2P4D9BbKVcLm91YunzPCUSiQQNDQ3tAaS8eF2j0czMz8/nVqxYERUWFiZoaGhgNBoNBg4cCKlUioiICNy48W8PHIlEArVajRkzZqBNmzbQ6/VITU3FmTNnmM6dO/tkZmZuqqur60MImfEK5sC2rq6u1FQ4Zt68eUzbtm0Z4JkmbPv27QgLC8O8efNgYWGB/Px8nD9/Hjt27MCoUaPMcVhjx44FAHTo0AEXLlyAlZUV7OzsUFxcDKFQaA6qZ1lW2qlTp+6jR48GpZSn1WrBsixOnTo1NCUl5YZIJMrl8XiDXFxcdM9jtXkikeikWq1+h1L6cvWAP43/BLFPP3jw4LvWFP1Pnz41nbtGtKDEec6VjiOEBNy4cWMSx3HTPT093aKiotC/f39s27atSRWfp0+fwsfHx6wP5/P5CA0NhZeXF7799lssXLhQsm3btrEajeYogPN/MPYnxcXFdOPGjVQikRCNRgOlUml2OBg7dmyTcJ6AgAD4+flh586d2LhxI/r164ebN2/i888/N3vImPzkBgwYgODgYPD5fNTW1uLy5cuora0FpdRcr4vjOFBKERISImBZtp2jo2NgfX09q1KpBHw+H3w+H5aWlmNqamr8CCHdKKW6P0+ef+MvE5tSWiAUCq/u3bt3oCmT0dGjR5GZmQlnZ2fU1taid+/eOH/+PAPg0e/0kwlgBSEkrrS09FSHDh0kKSkpmDlzpllLpVarcfXqVXMy+Bdhb29vTi4/ZMgQyYkTJ/6JPyA2pbSOENLns88++4oQwlcqlSkJCQlzAEikUmmLcVummDWNRoMBAwagvr4ev/76K2QyGXx8fPD48WO8++67TQL7ra2tMXbsWHPSvujoaEybNg0cx6GqqgoXL14EpZTRarVNitrU1dXh3LlzbENDQ2eGYR5YW1sLNRrNjyqVau3rODH8R6I4NRrNiKdPnzZ8/vnn2LlzJ7y9vVFYWIhx48ZBKpUiOztbRQj58UW34d/BVY1Gc/jGjRsKHx8fuLm5gVKKrKwsbNiwAf7+/q1mUmrbti1KS0vRtm1bGI3G5q4hLYBSml5bWzukpqamP4BVT548UcbGxlIXF5dW86zm5eWZ85vZ29ujrKwM1tbW5mS3LQX263Q65OfnY+7cuejevbs51MjOzg6TJk1Cx44d4enp2eQjsbKywsSJE+Hl5UVkMpnfpUuXPJydnZfgWU21P40/bfVqCZRStVqt7qhSqapqa2upqTZIdHQ0GhsbkZOTU6TVape+ykOeW9XeKSsre/+nn37K++ijj7QfffSRdteuXSWVlZV0yJAhzYwiWq0Wp0+fNlcAksvlIISI7ezskoVC4Qfk5Rtaf3aDQCC4WlNTQ/Lz8/Hzzz+3SHCtVmtW35oIfOrUKZw5cwYeHh7Iy8trdk9aWhrc3d3h5eUFo9GIuLg4HDlyBLdu3QKlFEOHDkVaWpq5Eq8JhBBERkZCpVJh2rRpsLKykgAIfan7zq/yfn+G2Km/d5FSmq/X631ramoOfPDBB3TdunWYPHmyKbDdhWXZ5ntv631RSulOtVrtrdPpXHU6natGo3ETCAS/3rlzp9nsnzlzBmKxGFu2bIGtrS1OnjxJfX19Bfv37w9xdXVdTQhpuZrLSyCE+LAsOzw/Px95eXlQqVTIzs5u1k4kEpkT4ZmMJkVFRaisrERJSUmLKbTLysrMfEZMTAzKysrwzjvv4M6dO0hKSoK1tTVEIlGLSe6ysrIQFhaGAwcOwM/PDzKZ7OVo/9RXeb//aDA+pbTOaDROqa6ujlq1apV69uzZuH79OpYuXSqRSqVDXqM/Simtev6PqtXqNRcuXFDn5+c3aVdQUICVK1diyJAhWLJkCaqrq8muXbvYIUOGYOHChRKJRNL3FR/JEkLA4/HAsqy5Mt+LkMvlKCwsVF2/fl3X2NgIOzs7eHt7m4urqlSqFvOwCQQCcz70J0+eYOPGjZg+fTrWrFmDnJwcGAwGqFSqFlNkl5aW4r333kOnTp2wcOFC8Hi8louk/AH+KyG7lNJzVlZWNzIzM/skJiZymzdvVioUioT/QL+phJDxmzZtOuzn58cEBQWJNBqNOTHc+++/j3Xr1kGr1RZt2bLFYeLEidzGjRsVCoXi2is+4onRaDzVtm3bKJFIRMvLy9nTp0/r6+vrZTKZDPn5+fr4+Hid0WjcTClVfvPNN4tHjx4tDgsLQ1VVFc6cOWMu+DJ9+vQmaTudnZ1x5coVjBgxAg4ODvjxxx8hk8mwd+9e2NnZ4cGDB3B0dDRniXwRGo2GLl261KBWq3mbNm1S6nS65nUzXgF/WYP2O+1tLS0tdzAME6TVao8pFIpPWwn5+dMghMgIIVMEAkFvg8HgaDQaIRQK3TmOE2g0mkNKpXK9hYXF93w+v51CoTigVqs3vCr3+vx87waAA3ALwDCRSDSLYRhrnU6XrtVqN1FKHz5vO0YkEn2m1Wr9WZZleDwewzBMrU6ne0QI6RoWFsa3srJis7OzVVlZWZRl2ayAgICA0aNHC0+dOoXi4mL4+PigR48e2Lx5M3V1daXvvvsu86JJ9c6dOzh69GgDIeRrmUzWU6lUJqjV6nUvqp7/q4F9f6MpCCEWeLZL1po+KkJIACFkPMuyNnq9/iGAQwAMQqHwOMMwfbp37y6wsLBgc3JyFBkZGYzRaFzMcdwISmlEp06dGIFAwDx48ECjVCorNRrNCErp/d95/n+W2H/jPwdCSEeWZSewLGuj1WrvAzhAKa17fs0XwFA821mSAVx7HZm6xef+vbL//4//xjZeByCVUtqXtBzp3xmtiwD/0Wu/M4bXvfb/9fH90djxHzVx/o3//+M/nvTub/x/F38T+38IfxP7fwh/E/t/CH8T+38IfxP7fwh/E/t/CH8T+38I/w8znYz+mGkyxwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "benzene_box = plams.packmol(plams.from_smiles(\"c1ccccc1\", forcefield=\"uff\"), density=0.875, n_molecules=40)\n", "plams.plot_molecule(benzene_box);" ] }, { "cell_type": "markdown", "id": "d0ce0756-ebab-4e0a-bc61-6341c34e7b61", "metadata": {}, "source": [ "## Equilibration MD simulation\n", "\n", "Short equilibration MD simulation at 300 K for 5 ps with the GAFF force field and automatic atom typing. " ] }, { "cell_type": "code", "execution_count": 3, "id": "429020cd-f9b0-493c-8ac6-6e059922c892", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[17.03|17:48:37] JOB equilibration STARTED\n", "[17.03|17:48:37] JOB equilibration RUNNING\n", "[17.03|17:49:35] JOB equilibration FINISHED\n", "[17.03|17:49:35] JOB equilibration SUCCESSFUL\n" ] } ], "source": [ "s = plams.Settings()\n", "s.input.ForceField.Type = \"GAFF\"\n", "s.input.ForceField.AnteChamberIntegration = \"Yes\"\n", "eq_job = plams.AMSNVTJob(\n", " settings=s,\n", " molecule=benzene_box,\n", " name=\"equilibration\",\n", " nsteps=5000,\n", " timestep=1.0,\n", " temperature=300,\n", " thermostat=\"Berendsen\",\n", ")\n", "eq_job.run();" ] }, { "cell_type": "markdown", "id": "ce4f92a1-cd3e-4780-9b1b-4f8dda814a74", "metadata": {}, "source": [ "## Production MD simulation\n", "\n", "Enable the \"Bin log\" for Pressure tensor and Time. This stores these quantities at every time step of the simulation, even if the structure is not stored. It is required to be able to calculate viscosity from the Green-Kubo relation.\n", "\n", "Here, we run a very short MD simulation of 200 ps. Typically, to converge the viscosity you would need longer simulations.\n", "\n", "The time step is quite large at 1 fs. You may want to decrease it (and increase the number of MD steps correspondingly)." ] }, { "cell_type": "code", "execution_count": null, "id": "dff078df-a770-4f29-a15a-eb829f8a7f97", "metadata": {}, "outputs": [], "source": [ "prod_job = plams.AMSNVTJob.restart_from(\n", " eq_job,\n", " name=\"production\",\n", " nsteps=200000,\n", " timestep=1,\n", " temperature=300,\n", " thermostat=\"NHC\",\n", " tau=100,\n", " binlog_pressuretensor=True,\n", " binlog_time=True,\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "id": "819b203e-1eda-4ce0-8192-f929c4186a45", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[17.03|17:49:35] JOB production STARTED\n", "[17.03|17:49:35] JOB production RUNNING\n", "[17.03|18:35:10] JOB production FINISHED\n", "[17.03|18:35:11] JOB production SUCCESSFUL\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "prod_job.run()" ] }, { "cell_type": "markdown", "id": "377baf66-a1b4-4030-a195-807273804656", "metadata": {}, "source": [ "## Calculate the viscosity autocorrelation integral" ] }, { "cell_type": "code", "execution_count": 6, "id": "bf7c637b-5307-4897-a5d2-80191f269e53", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[17.03|18:35:11] JOB plamsjob STARTED\n", "[17.03|18:35:11] JOB plamsjob RUNNING\n", "[17.03|18:35:29] JOB plamsjob FINISHED\n", "[17.03|18:35:30] JOB plamsjob SUCCESSFUL\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scm.plams.recipes.md.trajectoryanalysis import AMSViscosityFromBinLogJob\n", "\n", "job = AMSViscosityFromBinLogJob(previous_job=prod_job)\n", "job.run();" ] }, { "cell_type": "markdown", "id": "bcdb7459-6b2f-413d-bc53-53e7aace3044", "metadata": {}, "source": [ "## Plot the results\n", "\n", "The viscosity is the limiting value of the viscosity integral.\n", "\n", "The viscosity integral **must** converge to a near-constant value. If it is not close to converging, then you should run a longer MD simulation.\n", "\n", "To better estimate the viscosity, a double-exponential function is fitted to the viscosity integral curve. This makes it easier to get a numerical value for the viscosity.\n", "\n", "$\\eta(t) = A(\\lambda (1 - \\exp(-t/\\tau _1) + (1-\\lambda) (1- \\exp(-t/\\tau _2))$\n", "\n", "The parameters $A$, $\\lambda$, $\\tau _1$, and $\\tau _2$ are fitted and available as ``popt[0]``, ``popt[1]``, ``popt[2]``, and ``popt[3]``.\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "7a228b77-4854-4c96-a3bd-ea8be7ddcf78", "metadata": {}, "outputs": [], "source": [ "x, viscosity_integral = job.results.get_viscosity_integral()\n", "popt, x_fit, fit_viscosity_integral = job.results.get_double_exponential_fit()" ] }, { "cell_type": "code", "execution_count": 8, "id": "c9cbdb4e-bc98-4ca6-8ce3-a716d3c257bd", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEWCAYAAABFSLFOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABD8klEQVR4nO3dd3gc1dX48e/ZVbVlW5Z7l7uxARfkQrcBGxfAEJppoSUOvPjlRwgJBl5CeUliyBsSSACHAKEEMDXBBAIY04tx77jhKluWZbnJ6to9vz9mZK9kbZG0q3o+z7PP7szce+fsWJ6zc2fmjqgqxhhjTLR46jsAY4wxTYslFmOMMVFlicUYY0xUWWIxxhgTVZZYjDHGRJUlFmOMMVFlicXElIicLiLra1i3p4gcFhFviDKHRaRPzSOMPhHZKiLn1HccxtQXSywmKoLtTFX1S1UdWJM2VXW7qqaoqs9dx2ci8pNKZVJUdXPNom7cxPGwiOS6r0dEREKUP1tE1olIgYh8KiK9Im1LRNLdOgVuG+dUaruDiLwiIgdEZL+IvBywrJuIvCMi+0QkU0RuChHjWBHxuz8Y8kRkvYhcX/OtZOqDJRZjGq/pwIXAUOBE4DzgZ1UVFJH2wNvAvUAasBh4rRptvQosA9oB9wBvikiHgOVvA7uBXkBH4P8Clv0D2AJ0AqYAvxWRcSG+1y5VTQFaA3cCfxORwSHKm4ZGVe1lr1q/gK3AOVXMHwtkVir3S2AlkA88i7PD+Q+QB3wMtHXLpgMKxAG/AXxAEXAY+ItbRoF+7ufngSeA99y2vgP6Bqx7ArAeOAg8CXwO/KSKmLsChUBawLzhwF4gHugLfALkuvNeBlKr2hZuTA+F2B5dgbeAHJyd763V2ObfANMDpm8EFgQpOx34JmC6pfsdB4VrCxgAFAOtApZ/CdwUsF23At4q1pvi/ht1CJj3NPBSkDgrbB93Xg5wCU5SWgYcAnYA94fYNu2BfwMHgH1uvJ76/n/SXF52xGLqw8XAeJwd1vk4SeVunJ2BB7i1cgVVvQdn5zBDne6vGUHavgJ4AGgLbMJJSOW/2N8E7sL51b0eOKWqBlR1F/CtG2e5K4E3VbUUEOB3OEnhOKAHcH9E3zyAiHiAd4EVQDfgbOA2ETnXXX6aiBwI0cQQt265Fe68sGVVNR/4IaB8qLaGAJtVNS/I8jE42/MFtxttkYicWf41K72Xfz4+xPdyCol4ROQiIBVYhfND5Mfu9BTgZhG5MEj1XwCZQAecHy534yQ4UwcssZj68GdVzVbVnTjJ4jtVXaaqxcA/cY4OauptVV2oqmU4RxLD3PmTgTWq+ra77HGcrptgXsFJUrjnGqa581DVTao6T1WLVTUHeBQ4M2hLwY3E+SX/oKqWqHOu6G/uulDVr1Q1NUT9FJyjr3IHgZQg51kqly0v3yqCtsLV7Y5z1PIp0Bn4A/COiLR3k9HXwL0ikiQiI3ASdosQ36urm1D3AvcB16jqelX9TFVXqapfVVfidM8F2+6lQBegl6qWqnOuzxJLHbHEYupDdsDnwiqmU2rRdmCyKAhoqytO9wkA7k4mM0Q7bwIni0hX4AycX7tfAohIRxGZIyI7ReQQzjmE9jWItRfuTrT8hfPLulOE9Q/jnIco1xo4HGQHWrlsefm8IMsD2wpXtxDYqqrPujvxOTjb+lR3+VVAb3feUzgJP9S236WqqaqapqrD3PYQkdHuBQQ5InIQuIng2/33OEesH4nIZhGZGWJ9JsossZjGpDa/OLNwflkDR45CugcrrKoHgI+Ay3C6wV4N2GH/zo3lRFVtDVxNxa6eQPlU/HXeOeDzDmCLuxMtf7VS1ckRfqc1OCfbyw1154UtKyItcc4VralqeaW21gB9RKRVkOUrCfFvo6rbVPU8Ve2gqqNxuiIXhvhewbwCzAV6qGobYDZBtruq5qnqL1S1D0536+0icnYN1mlqwBKLiaZ4t7uj/BUX5fazgZres/IecIKIXOjGdQsVd/JVeQWnT/9i93O5Vji/4g+ISDecixGCWQ5MFpE0EekM3BawbCFwSETuFJFkEfGKyPEiMjLC7/Qizg6zm3tk9QuciwWq8k/geBG5WESSgF8DK1V1Xbi2VHWD+z3uc/9dL8K5cuytgLbbisi17ne4BOec0dcAInKciLQSkQQRuRqn2+zRCL9joFbAPlUtEpFROAm/SiJynoj0c39AHMK58MNXg3WaGrDEYqLpfZxukfLX/VFu/zHgEvc+icerU1FV9wKXAo/gXM01GOeS2+IQ1eYC/YFsVQ08sf0AMALnPMN7OJfaBvMSzonurThHQEcu8VXn/pzzcc4DbcE5p/AM0AaO3Fx6OETbf8U5+b8KWO3G8tfyhSKyRkSucteVg5MgfwPsB0bjnsuJpC23bIZbdxZwidsmqroPuAC4w90mM4Gp7jYHOBfY7Na9CZhYXrea/gt4UETycBLj6yHK9se5wvAwzoUYT6rqZzVYp6kBsfNZpjlyr8jKBK5S1U/rOx5jmhI7YjHNhoicKyKpIpKIc5JcgAX1HJYxTY4lFtOcnIxz78ZenC6oC1W1sH5DMqbpsa4wY4wxUWVHLMYYY6Iq2peDNirt27fX9PT0+g7DGGMalSVLluxV1Q7BljfrxJKens7ixYvrOwxjjGlURGRbqOXWFWaMMSaqLLEYY4yJKkssxhhjosoSizHGmKiyxGKMMSaqLLEYY4yJKkssxhhjoqpZ38dijDHR9H3WIfbkFbNjXwEXDe9Gy8TmuYttnt/aGGOi7EBBCZMe+/LI9P/8azUL7z6bjq2T6jGq+mGJxRjT7BWV+jhQUMrPXlrMisyDAEwY3IkHpg6hS5vksPXLfH6GPTjvmPmjfjufLb+bjPMgy+jJLy5jW24By3ccYEjX1nRNTSb7UBEHC0s5pW+7qK+vuiyxGGOaJJ9feei9tfx7ZRY5ecW0TorjUFEZU4d1ZXTvdjw6bz17D5cErf/R2mw+WpvN+MGdePKqEcR7g5+Sfu7rLUGXzf58Mzed2SdqO/tFW/fx6Ecb+HZzbtAyFwztyg2n9WZI19bEez1kHSwkc38hx3dtQ3KCNypxhBLTYfNFZCLO42S9wDOqOqvScnGXTwYKgOtUdWmouiKShvN413Scx71epqr73UewBj57/ERghKouDxZfRkaG2lhhxjRuc1fsonvbZEb0bHtk3p/nb+QP8zZEdT1L/ucc2qUkHjO/uMzHwP/5IGTdq0b35DcXnRDReopKfTz12Q+0TPRy3old6dQ6Ca9H+Plry/nnsp01ir1ciwQvK++bQFyIJBkJEVmiqhlBl8cqsYiIF9gAjMd5BOwi4ApVXRtQZjLw3ziJZTTwmKqODlVXRB4B9qnqLBGZCbRV1TsrrfsE4B1V7RMqRkssxjRO+cVlfLcllxueP/r/t3vbZF796RgWbd3H7a+viMl6v555Ft1SK3aNpc98L+L63951VtCutTKfnzvfWsVbSzMrzO/fMYU4r4fvsw5VP+AqJMd7ee/W0+jTIaXGbdRnYjkZuF9Vz3Wn7wJQ1d8FlPkr8JmqvupOrwfG4hyNVFm3vIyqZolIF7f+wErr/q1TXO8JFaMllubpmS8389B73wNw61n96JKazLSRPeq9X9qE9tC/1/LW0kz2F5TWaxzPXZfBWYM6sXrnQe6fu4bF2/ZXWH7LuL6cNagTFz/1TZX1/379SA4VljK0eyrb9xWQ3q4luw4WMu3pun1K9oaHJpEQV7Mjl3CJJZbnWLoBOwKmM3GOSsKV6RambidVzQJwk0vHKtZ9OTC1qqBEZDowHaBnz54RfRHTdJSU+Y8kFYDHP9kEwF1vr+JnZ/bhrknHVSi/bPt+nv1qC49NG47X0/ATj6ry6sId9O3QklG905pEssw6WMiVf/uOLXvz62aFqnhQRBWPKoIi6j/y+danv0Bwl6nSNuDzPZMHMmVICvjzWHBNfy5+8utj2nroj/888rm8ngflePUDIG4MAs66K02j5fMVUcAtI+5BQvl8wZ1211V5Wc6CdnQ7bUxMNmEsE0tVf9GVD4+ClYmkbtUrFRkNFKjq6qqWq+rTwNPgHLFE0qZpGnYdKOSUWZ8EXf7XzzcDMLp3GgUlPopL/fziDadL5dN1e5h8Qhd+f+nQOom1OhZv3cej8zawLbeAnQcKKyzbOmtKncSgPh9aVIS/uBgtLMRfUoKWlKKlpXy6eie92yTQLSUOKXPmaWmpu7zk6HRpKVpScbowv4gPlu/gopIyvOp3Xz48fuez58g89+V3l6seU7ZCef/Rtrzqd3a+qngj280E9xH8EDD5TO1ai6mFf9jERY0wsWQCPQKmuwO7IiyTEKJutoh0CegK21OpzWnAq7WM3TQiZT4/WQeLWLbjAIeLyrj7n6t46qoRHNelNentWx4pFyqplPvr55uPJJhA+SU+3liS2aASy/78Eob/77GXuAZKn/keZwzowHebc3ly2on0aiEc3HuAIalxaEEBvvx8Sg7l8emybXy5YhsJpUVMHZBKKy2lQwL4C4vQ4iLnvTxxFBXiLyp2pouc+VoavHuqt/se6WlnFQ+exASIi2N/iXK8x4tPvPg8HnziwS/Ou0887jwvfvFQJl5K4uKPlDtaxuvWkyNlj9Z12vOLBz+CiuAXoVVyAgeLyvCLc2xw9N2DBs4TOVLvRyf1YFSf9uDxIB4BjwfEw7x1e/jXiiznGEOkQpuB7akIx3VpzfQz+tCxdTKIsOtgIY/N38SW3AIyeqcx+YQuJMR5GdiltXM0GvjCec/cX8iGnMP84aMNbtsAgkr5sQ2oCC3btuGiCP9NqiuW51jicE7An43zN7UIuFJV1wSUmQLM4OjJ+8dVdVSouiLyeyA34OR9mqr+ym3PA2wHzlDVY/cOldg5lsYr62AhKzMPsjknn3W7D/HO8sq/WRw901rQq10Lvty4Nyrr/dmZffjF+IE17puuLX9hIb79+zmUs4/rHp9Pq9JCUkoKnPfSAlqVFJJSWkirkgJSSgtIKS2kRWkRyWXFJPrLIlpHqXgpikugNC6B1NQUklu1RJKS8CQmIslJeBKTnOmk8vfEo9OJSXiSk9hT6OfBDzdS6omj1BNHmcd75FXiicObEE++30OZx+uW8VLmlvNL5Nu2ZYKX+y4Ywq/eXFnTTVrBdaekc8u4fnRolYjfr5T5lZlvreTtCK7G2vzbyXiCdJf6/cr0lxbz8feVfwcf9afLh3H2cR1plRRfYf7h4jLyi8voVM0bLUt9fuZ/n83S7QcY3iOVzzfkMGfR0TMM3VKT+XrmWdVqs1y9nbx3Vz4Z+BPOJcPPqepvROQmAFWd7V5u/BdgIs7lxter6uJgdd357YDXgZ44SeRSVd3nLhsLzFLViI7vLLE0Lh+u2c3PXlpS32FwxaienD+0C33ap9CxVWLQnUmkfIcPU7Z7N2V7c/Hty6Usdx9luXvx5e6jLDcXX24uZfucz1pQELwd8ZAXn8zhhBYcjk8mz30viEuiIC6RQvdVEJdIYXz59NFlBXFJFMYlUuqt2JHxxS/H0bNdC1SVpdsPMKJn6jHnbt5ckskD767hqtG9GD84+Inr2vrluQMRgXOO60Spz8+Qrm0A5xLdu99eVSEBdGiVyJ0TBwHwzaa9DO+Zitfj4e5/riIhzkPvdi2574LBHNe5NXFeOWaHHij7UBG3zVke9N6RpfeOJ61lQtj4567Yxa2vLmPS8Z2ZeHxn9hwq5pzBnSgp8zOwc6vqbIoa+XhtNj950dnnXTm6J7+N8BLoyuo1sTR0llgaj/W78zj3T1/UdxhVSmuZwIxx/WiZ6OXykRUvCPHl5VG6axdlu3dTujubsmz3ffduSrOdd39+FSelPR68aWnEtWtHXLs0vGntiGvXDm+7NP6yZC8bC+RI4jickExefAsK4xLdLpHomzlpELP+sw6AVklx/PGyYTz/zVYWbtlHic8fk3VWtuK+CbRJDr7zB6d7cG3WITwijOiVSmLcsTcDlu/zqnthg6oy49VlvLcy68i8K0f35L7zB1e5noZq6958cg4XM7xHao3vZ7HEEoIllsYh+1ARo387P+rtnj+0K++uqLoLLVIe9dO+8ACd8/fROT+XLgW5XNMzDm92FqU7duA7cKBCeRUPcR3aE9+lM54OnYjr3InDrdLo3Lcnce3bs+ywsEOTWJsHB4t9vLtyF/efP4T+HVPweoT/rN7N899srVXMDdXPzxnAKwu3kX2ouML8+88fTP9OrTi1X/t6iqyiAwUlLNtxgJP7tCMpvvEklGiyxBKCJZaGb8vefMb932dRbXP21ScxvGcqnVon8dGa3WzfV1DhEuTKvvzVOLxlJVx53xv0zNtDj7xsehzeQ8+8bLofziHe7ztStkw8ZLdIw9e5C9sT2rLO04rsFmnsTU5lb3Ib9iW2wu/xcmq/dny96Wi3yuyrR/DSgm0V5lXHsB6pLN9xoEZ1Y+nG03pz73mDyckrJutgIRf85esKyxPjPLx18ym0SY6nR1oLikp9bNpzmPP+/BXp7Vrw6R1jm8Ql002NJZYQLLE0bF9t3MvVz34X9XbX/e/EY35p+v3KybPmk32wiPaFB+l7cCeTkg5xVtwBSjZupDQzE/xOl48PIbtlGttbdWJHq47sbNmB3S3TyGrRjr3JbfB76u5X7KTjO/PHy4eRFO+lzOenzK/k5pfQtU0SJz30Mfvyg4+FVdmsH51A5zZJzFubzcvfba92LKN6p7Fwyz7AOak+/xdj6dym4gnn7ENFtEqKY+f+Qvp3Cn5OIa+olHivp9keETR0llhCsMTScH26fg/X/31RVNvslprMaz8bQ/e2LQCnz7xs1y4KV6ygcPUaDq9eQ8n6dXDQGd0WERLS00kcOJDEPn1I7NeXhL59WepL4SdzVvHhbWew+1ARl87+NqpxRuqMAR148YZRQZcXlvh4e1kmfdqn8ORnmygq9ZGRnsZTnx290yLeK1w1uhf3XzCkyjYmPfZlyKFExg3swIXDuzF+cCdaJMRxsLCU/fklFS7zNk2PJZYQLLE0TJ+u28P1z9csqYzpk8Y1Y9L5YkMOry3eQa92LWidFE9xmY97xvdlZMkeCpctp3C58yrLyQFA4uNJHDCApMHHkXjccSQNOo6kgQPwtIx8B7nnUBGjYnAuqCqPXHwil43sEb5gFVSVwlIfry/awdVjeoU8gevzK1Me/5J1u/M4vX97/u/SocR7PbyzfCc/Pjm9UYxGYKLPEksIllganqJSH4PuDT5SbILXw4heqSzY7HS5nHNcR8r8yjVjejFuYMcKl/76S0ooXL6cggXfkb9gAYWrVoF7I198jx4kDxtG8rChJA8bRtKAAUh86CuOItXnrvfwR+G/1UXDu9EqKY4Xv91WYf7z149k7MCqRjIypm5YYgnBEkvDkV9cxuzPf+DP7thdVXnowuO55KTuJMV7WbvrEN3TkmkdcO+B+v0Urf2e/K+/puC7BRQsWYoWF4PHQ9KQIbQYNZIWI0aQPHQoce1jd4WRz6/sOlDIL99cwc1j+5FfXMZ/vbyUyzN60L9TCk98uqnCQIrJ8V6Ky3z4FT69YyztUhKI8wgtEpz7SfbkFbH7YBGqzg2fbSO4X8KYWLLEEoIlloYhkqFJkuO9fP+/E4+Z78/PJ//bb8n77DMOf/45vhznDvvEAQNoMWY0LceMoUVGBt7WrWMSe6SKSn0kxnkQEXx+ZdOew/Rq55zrSYr3UlTqY39BSURPKzSmvtXn6MbGhPXmkkzueCP0szMuy+jOg1OPPzJdtm8feR/NI2/ePAoWLkRLS/GkpNDy9NNIOfNMUk4/nbh27WIderUEXt3k9cgxd1knxXstqZgmwxKLqTePfrT+yLD1wTxyyYlcelJ3fAcOsH/ePPI++ID8Bd+B309Cr160vfpqUsaOpcWI4VE7R2KMqR1LLKZe/GvZzrBJ5aHJA5i4ZxU7fvog+d9+Cz4fCb160W76T2k9aRKJAwbYzXPGNECWWEydUlV63/V+qAL0P5DJL/wb6H3Hg+w8dIi4rl1od8MNtJ48icRBgyyZGNPAWWIxdSbrYCEn/67qZ6IklhVz1o6lnLflG/ocykISE0kZP57Ui39Ei9GjEU/9DFNvjKk+SyymTuQXl1WZVLrk7+W8zd8wYdtCUsqK2JbWg3b3/pp250+p9yu5jDE1Y4nFxFxRqY8h931YYd6Q3C1cuuETRmavwy/CV11PZMwvbuLcSWdYV5cxjZwlFhNTfr8efY6KKiftWc/lG+ZzQu4WDia0ZM7As3kv/WSe/cUkhvdsW7/BGmOiIqaJRUQmAo/hPAXyGVWdVWm5uMsn4zxB8jpVXRqqroikAa8B6cBW4DJV3e8uOxH4K9Aa8AMjVbUolt/RBOfzKwP+5z/4fH5OyVrNtPUf0//gTnKS2/DUCVP5sNdoiuMS+OrOcUcGhjTGNH4xSywi4gWeAMYDmcAiEZmrqmsDik0C+ruv0cBTwOgwdWcC8wOeeT8TuFNE4oB/ANeo6gr3EcalmDrn9ys/f3057yzfxbA9G7h+7fsMOJBJZsv2/HH4ZXzSYwRlHudP786JgyypGNPExPKIZRSwSVU3A4jIHGAqEJhYpgIvqjOuzAIRSRWRLjhHI8HqTgXGuvVfAD4D7gQmACtVdQWAqtbsiUkmYqpKcZmfzTn5KMrzX2+lb8cUPli9m4KVK/nt2v8wPGcj2cmp/GH45XzS8yT8cvTqrmE9Url5bN96/AbGmFiIZWLpBuwImM7EOSoJV6ZbmLqdVDULQFWzRKR8mNcBgIrIh0AHYI6qPlI5KBGZDkwH6NmzZ+XFJkIfr83m2825PPvVlgrz2xce4IY17zEucxkHEloy+4SpvJ9+MqXein9q152Szj1TjqvLkI0xdSSWiaWqS3sqj3gZrEwkdSuLA04DRuKcr5nvDpRW4QEZqvo08DQ4g1CGadNUUlTqI6+ojJ+8WHHwznhfKT/a9AXTNnyMR5VXBp7Dm/3GUhifVGU790w5jvgQzwExxjResUwsmUDgk4i6A7siLJMQom62iHRxj1a6AHsC2vpcVfcCiMj7wAigbp681ExcMvsbVu+s+ETBkbvXcvPKd+hSkMtXXU7gmePPI7tl8EEgVz9wriUVY5qwWP7vXgT0F5HeIpIATAPmViozF/ixOMYAB91urlB15wLXup+vBd5xP38InCgiLdwT+WdS8XyOiYLApNKmOI+Zi/7Bgwueo9Tr5a5TpvOb0deGTCof334GKYl2lbsxTVnM/oerapmIzMDZ4XuB51R1jYjc5C6fDbyPc6nxJpzuq+tD1XWbngW8LiI3AtuBS906+0XkUZykpMD7qvperL5fc6OqR8+nqHLWjqX8bNU7JPuKeXHQubwxYNyRK72CmffzM+jXsVXIMsaYxs8e9GUP+orIgs25THt6AW2LDnHbstcZlb2OtWm9eGzYpWxv3TmiNrbOmhLjKI0xdcEe9GVqLetgIdOeXsDorDX8fNnrJPpKeOqEqfy7z6kVLh8OZmR6W9646ZQ6iNQY0xBYYjEh5ReXMfZ/P2DG6neZsvVbNrXpxsMZV5HZqmP4ysCATik8etmw2AZpjGlQLLGYkB7887s8/tmf6H44hzf6j+XF4yaGPZdSbtrIHsy6+MQYR2iMaWgssRgAtuzNJzney+5DRSzfvp8vNu7F9/GH3LbsdQrikrjn1Oks79A/4vb+fv1Ixg2M7KjGGNO0WGJp5opKfazeeZBLZn97ZJ7X7+Mnq//NhZu/ZE1aOr8deQ37kttE1N6UE7swYXAnxg7oEKuQjTENnCWWZuyJTzcx//tslm4/cGRem+I8/mfhixyfu4V/9j2dZ4ech8/jDdnOsB6p/Pr8wYzo2RZVteepGNPMWWJphhZszmXhln08Om9Dhfk9D+3mgQXPklp8mFkZV/F59+Fh23r44hO4fOTRMdcsqRhjLLE0I6rKG0sy+dO8Dew6WPExNcP3bOCehS9S7I3nztNuZkPb8AN0vnfraQzpGlkXmTGm+bDE0ox8tWkvv3pz5THzJ25dwIwVb7O9VSfuG3MDOS3CP8nx49vtLnpjTNUssTRhz321hV7tWrA1t4BP1+3hq017KxZQZdqG+Vz7/Qcs6jiQWSOvoSDIaMSBfvejEyypGGOCssTSBGXuL2DJtv08+O8QY3Cq8pPV73LxD1/wcY+T+OPwy/CHOUkPsOGhSSTE2cjExpjgLLE0QVc98x3bcguCLvf4ffy/5W8yYfsi3ulzGn894QI0gqFZAOK9dnLeGBOaJZYm5quNe0MmFa/fx52LX+b0XSv5x8DxvDxoAkRwJZeI83x6u+rLGBOOJZYmYltuPmt3HeLml5cGLePx+/jlklc4fddKnj7+fP7Z78yI2z9rYEduOtOeT2+MCc8SSxPw+YYcrn1uYcgyHr+PO5bO4cydK6qdVG49qx/XnpJeyyiNMc2FJZZGTlXDJxX1c/vS1xiXuYxnh0ypVlIBuH3CwNqEaIxpZiyxNHKzP98cuoAqM5a/xdmZS/n74Em82X9cxG1fNbonx3ezGyCNMdUT0+tGRWSiiKwXkU0iMrOK5SIij7vLV4rIiHB1RSRNROaJyEb3va07P11ECkVkufuaHcvv1hCoKg9/sC5kmWu+/5BJ277j1QFn8/qAsyNqd+zADnzxy3H85qITuGJU+DvwjTEmUMyOWETECzwBjAcygUUiMldVA2+umAT0d1+jgaeA0WHqzgTmq+osN+HMBO502/tBVYfF6js1JOt357Fu96GQZc7f/BVXbviY//QazYvHTQzb5pQTuvDo5UOJ93jweOzqL2NMzcSyK2wUsElVNwOIyBxgKhCYWKYCL6qqAgtEJFVEugDpIepOBca69V8APuNoYmkWPlidxU3/CH71F8DpO5dz08p3+KbzEP4y9EdhLyl+9adjOLlvu2iGaYxppmLZFdYN2BEwnenOi6RMqLqdVDULwH0PfJpUbxFZJiKfi8jpVQUlItNFZLGILM7Jyanud6pXh4vL+HhtdtikMiR3C79c8ipr26Xz8Mirw95R//2DEy2pGGOiJpZHLFX9RNYIy0RSt7IsoKeq5orIScC/RGSIqlboL1LVp4GnATIyMsK12aCceP+H+MNE3Ck/l3u/e57sFmk8MPp6SrzxIcsP7ZFKckL4oVyMMSZSsUwsmUCPgOnuwK4IyySEqJstIl1UNcvtNtsDoKrFQLH7eYmI/AAMABZH5+sc64PVu5n9+Q88csmJtEqKo0ub5Ki06/crIkefbbJixwGu+/vCsEmlRWkR9y94Do/6uW/MDRxOaBGy/PcPTrSkYoyJulgmlkVAfxHpDewEpgFXViozF5jhnkMZDRx0E0ZOiLpzgWuBWe77OwAi0gHYp6o+EemDc0FAmGtxa2fGK0sp8ysT/vgFAP99Vj/++6z+tRqkcX9+CcP/dx4jeqbiV9i05zCHi8vC1vOonzsX/4Meh3O455Sfsisl9KOBX/7JaEsqxpiYiFliUdUyEZkBfAh4gedUdY2I3OQunw28D0wGNgEFwPWh6rpNzwJeF5Ebge3Ape78M4AHRaQM8AE3qeq+WH0/OPZ8+J8/2USHVon8+OT0arVTVOpj9uc/cGlGD77f5fTcBT4uOBLXr3mPUdnr+PPQH7GiQ/+QZZ+7LoNT+7WvVvvGGBMpcS7Iap4yMjJ08eKa95QNuOc/lPj8x8zv1DqR7+4+J6I25izczr+W72TB5prnwFN3ruR/Fr3Iv3ufzBNDLw5Z9rFpw5g6rPI1FMYYEzkRWaKqGcGWR3TE4t6E2BUoBLaq6rF7U3NE9qHisGVUlTK/MvPtVbVaV7fDOdy+7DXWte3J08dPDVn2pF5tLakYY2IuaGIRkTbALcAVOCfTc4AkoJOILACeVNVP6yTKhqoW9xD+8eONPD5/Y61Wn1hWwj0LX6TU4+W3I6+h1Bv8d8K4gR349flDarU+Y4yJRKgjljeBF4HTVfVA4AL3ct5rRKSPqj4bw/iarDkLt9euAVVmrHiLXod2c+8pPwn6nPqTerXlylE9ufik7rVbnzHGRChoYlHV8SGWLQGWxCSiRiTUAcsHq3dzct92tEmueB+J36/c+85q9uSF7y4L5ZztizlnxxJeGjSBpR2Djz786/MGM7RHaq3WZYwx1RH2ulgROVVEWrqfrxaRR0WkV+xDa9xu+scS/t+cZcfM33WwkJe/q93RSpfDe7l51b9Y0b4vcwYGv0hgw0OTLKkYY+pcJDdcPAUUiMhQ4FfANpwusmYv3FN6P1ufQ/rM9ygLuHIs+1BRrdbp9fv41ZJX8ImH/xtxBf4Qz6qvzf00xhhTU5FcFVamqioiU4HHVPVZEbk21oE1JQ9/sI5T+rUnr6iMW1899iimOq5cP49B+7fzm5HXsLdFapVlbh7bl9PtPhVjTD2JJLHkichdwNXAGe6Q9qEHoGomJMLLwv725Rb+9uWWWq9vSO4WLl8/n496ZvBVt6FBy11/ajodWyXVen3GGFMTkfSVXI4zBteNqrobZ5Th38c0qkYiXFdYNCWWFfOLJXPY0yKN2SdcGLJsnMe6wIwx9SfsEYubTB4NmN6OnWOpc9et/YAuBbn88rSbKYwPfTSSkmhPnDbG1B/7aVsLdXXAMjh3Cxds/oq5vU9hdfu+Qcu1T0lg66wpdtLeGFOv7KdtA5fgK+Xny14nJzmV5wdPDlru8oweXDXGnk9vjKl/llhqQergJMtV6z6i++Ec7j5lesgusIcvOTHmsRhjTCTCJhYR6Q/8DhiMM1YYAKraJ4ZxGaDvgZ1cvPEzPug1imUdB9R3OMYYE5FIOuP/jnOTZBkwDufE/UuxDKqxiOXxiqifGSve4lBiS54dcl4M12SMMdEVSWJJVtX5OM9u2aaq9wNnxTYsM2HbIgbt384zQ84P+Yjh9imJ3DHBjmaMMQ1HJImlSEQ8wEYRmSEiFwEdI2lcRCaKyHoR2SQiM6tYLiLyuLt8pYiMCFdXRNJEZJ6IbHTf21Zqs6eIHBaROyKJsVZidMjSujifG9a8x6p2ffikx4iQZccP7sSMs0I/MdIYY+pSJInlNqAFcCtwEs4d+GGHdHHv0H8CmIRzfuYKERlcqdgknGfT9wem43S5has7E5ivqv2B+e50oD8C/4nge9VarLrCrlv7Pi3Kinhi6EV1exemMcZEQagHfXUE7gb6AauA36nq9dVoexSwSVU3u+3NAaYCawPKTAVeVOf5yAtEJFVEugDpIepOBca69V8APgPudMtdCGwG8qsRZ4MyaN9WJm37jjf7ncm21l3Clm+Z4K2DqIwxJnKhjlhexNlB/xlIAR6vZtvdgB0B05nuvEjKhKrbSVWzANz3jgDu0P53Ag9UM84ai/blxqJ+bl75L/YmteHlQRPClr/tnP7cbudXjDENTKjLjTur6j3u5w9FZGk1265qr6sRlomkbmUPAH9U1cOhdvgiMh2n242ePRvWDYVjM5cx4EAmvx9xBUVxiWHL33aOJRVjTMMTKrGIe2K8fC/tDZxW1X1h2s4EegRMdwd2RVgmIUTdbBHpoqpZbrfZHnf+aOASEXkESAX8IlKkqn8JXKGqPg08DZCRkREuWYUUzQOWxLISrl/7PhtSu/Npj+HRa9gYY+pYqMTSBufxw4G7z/KjFgXC3SC5COgvIr2BncA04MpKZeYCM9xzKKOBg27CyAlRdy7OxQOz3Pd3AFT19PJGReR+4HDlpNKQXfTDF3QoPMgjJ12Jhnh4F8BPT+/NWYM61VFkxhhTPaGeeZ9em4ZVtUxEZgAfAl7gOVVdIyI3uctnA+8Dk4FNQAFwfai6btOzgNdF5EZgO3BpbeKsjWgdsLQtOsRlGz7hqy4nhBxkstyEIZ0ZmZ4WpbUbY0x0hboqLF1Vt4ZYLkA3Vc0MVkZV38dJHoHzZgd8VuCWSOu683OBs4Ot0y1zf6jlDc01339InN/Hc0OmRFTeY5cgG2MasFBdYb93b4x8B6dLLAdnrLB+OEO7nA3ch3OepFmKxlVhPfKymbBtIe/0PZ2slPCPE/7p6b0Z1iO11us1xphYCdUVdql7U+JVwA1AF5zuqu9xjiR+o6pFdRJlAxWN44Zrvv+Q4rgEXhsQ8iDsiHumVL7H1BhjGpaQoxur6lrgnlBlTM31PZDJ6btW8vLA8RxKbFnf4RhjTFTYowZrobY9Ydd+/wF58cm83e+MiMrferaNCWaMafjsQV/1ZHDuFkZmr+PZIVMoiE8OW37rrMhO7BtjTH2zI5ZaqeEhiyrXrf0P+xJb8W7vU6MbkjHG1LOwiUVE3hKRKe4VYiYKhuVs5ITczcwZeDbFcQn1HY4xxkRVJMniKZy73jeKyCwRGRTjmBqNmp5juWL9x+xNas0HvcaELXvWoI7WDWaMaVTCJhZV/VhVrwJGAFuBeSLyjYhcLyLxsQ6wqRmSu4UTczfzZv9xlHrDn+Lq3zGlDqIyxpjoiah7S0TaAdcBPwGWAY/hJJp5MYusESj1+atd5/L1H3MgoSUf9BodtuxDFx7PHecOrEloxhhTb8L+ZBaRt4FBwEvA+eXPQgFeE5HFsQyuIXvi000cKCitVp3++3cwcs96nhs8OaJzK1eP6VXT8Iwxpt5EcrnxM+64XUeISKKqFqtqRoziatDyi8v4/Yfrq13v8g3zyYtP5r3ep4Qt++WvxtUkNGOMqXeRdIU9VMW8b6MdSGNSVOqrdp1eh7I4NWs1c/ucRkF8UtjyLRPtFiNjTOMUanTjzjiPA04WkeEcvWmjNdCiDmJrsGoy+OSlGz6l0JvAO31Pi6h8m2S7LsIY0ziF+ll8Ls4J++7AowHz84C7YxhTk9O+8ABn7lzOu31OJS8h/JhgGx6ahNdjQ+MbYxqnUKMbvwC8ICIXq+pbdRhTk3PB5q8QVf7V5/SwZVskeEmIs3tRjTGNV6iusKtV9R9AuojcXnm5qj5aRbVmoTrHEsmlRUzasoCvu57Anpbhn/poRyrGmMYu1E/j8j6bFKBVFa+wRGSiiKwXkU0iMrOK5SIij7vLV4rIiHB1RSRNROaJyEb3va07f5SILHdfK0TkokhijLXx2xeRUlbE2/3OjKh8iwRvjCMyxpjYCtUV9lf3/YGaNCwiXuAJYDzOUyYXichc9xkv5SYB/d3XaJzhY0aHqTsTmK+qs9yEMxO4E1gNZKhqmYh0AVaIyLuqWlaT+EN/t8jKedTPhT98yZq0dNanhb8n5bZz+nPeiV1rGZ0xxtSvSAahfEREWotIvIjMF5G9InJ1BG2PAjap6mZVLQHmAFMrlZkKvKiOBUCqmxRC1Z0KvOB+fgG4EEBVCwKSSBKgEcRYIxJhZ9jJu1bTpWBfxEcrt50zgH42hIsxppGL5CzxBFU9BJyHc/QwAPhlBPW6ATsCpjPdeZGUCVW3U/nd/+57x/JCIjJaRNYAq4CbqjpaEZHpIrJYRBbn5ORE8DVq7kc/fE5Wi3Ys6DIkpusxxpiGJJLEUn5DxWTgVVXdF2HbVf2sr3wUEaxMJHWPLaD6naoOAUYCd4nIMXciqurTqpqhqhkdOnQI12SN9TuQyeB923in72n47YkDxphmJJI93rsisg7IAOaLSAegKIJ6mUCPgOnuwK4Iy4Sqm+12l+G+76m8YlX9HsgHjo8gzmrTCHrZpmz5hiJvPB/3CD/qzYBOKfzlyuHRCM0YY+pdJMPmzwROxjkxXoqzw658rqQqi4D+ItJbRBKAacDcSmXmAj92rw4bAxx0u7dC1Z0LXOt+vhZ4B8AtG+d+7gUMxBnmv86llBQwNnMZn3YfQX5C+McOD+/R1k7aG2OajEhGN44HrgHOcIcy+RyYHa6ee3XWDOBDwAs8p6prROQmd/ls4H2cLrZNQAFwfai6btOzgNdF5EZgO3CpO/80YKaIlAJ+4L9UdW/4TVB9GuaA5Zzti0nylfLvCAabBPBYT5kxpgmJZKTDp3DOszzpTl/jzvtJuIruqMjvV5o3O+CzArdEWtednwucXcX8l3CG9o+5UHlF1M+ULd+wJi2dzamVr1UIWisaYRljTIMQSWIZqapDA6Y/EZEVsQqosRues5Hu+Xt5ZdCEiOu0bWEDThpjmo5IEotPRPqq6g8AItIHqP648U2IhugLm7L5Gw4ktOSrridG1NYDFwzh8pE9whc0xphGIpLE8kvgUxHZjNNn0wv3XIipqH3BAUbvXssbAyJ7nj3AtaekxzYoY4ypY2H3fqo6X0T641xlJcA6VS2OeWQNWLDjlfHbF+FF+TCC59kbY0xTFcmQLrcAyaq6UlVXAC1E5L9iH1rDVVVPmKif8dsXsbx9P3a3bFf3QRljTAMRyYWuP1XVA+UTqrof+GnMImqkTtz7A10K9vFRr1H1HYoxxtSrSBKLRwKexeuOPJwQu5AavqruvD9320IOxyXxddcTIm4n0R7oZYxpgiI5w/whzg2Js3FOL9wEfBDTqBq6SnmlZUkhp+xaxUe9RlHijezS4SevGsEZA2I3VpkxxtSXSBLLncB04Gack/cfAc/EMqjG5tSsVST6y/i4Z/hxwcq1SoojJTGyK8eMMaYxieSqMD/OEC6zRSQN6K6qzfs+lkrTp+9cQVaLNDak2v0oxhgTyVVhn7kP+koDlgN/F5Fm+7z7ylqV5DMsZyNfdhsa+aMljTGmCYvk7HEb90FfPwL+rqonAefENqyGLfBy4xF7NhCnfr7pUr0R+sMNZGmMMY1VJIklzn3uyWXAv2McT6MQeFXYiXt/ID8uiY2p3SOu3zopjqE9UmMQmTHG1L9Izh4/iHNl2FequsgdK2xjbMNqPE7c+wOr2/XG7/FGXGfl/efGMCJjjKlfkZy8fwN4I2B6M3BxLINq6Mq7sRLLiul+OIdPuo+o34CMMaYBCZpYRORXqvqIiPyZKobHUtVbYxpZI9DjcA4A21p3rudIjDGm4Qh1xLLWfV9cF4E0JuVZtteh3QBsb9Wp/oIxxpgGJlRimSgi+1T1hZo2LiITgcdwHi/8jKrOqrRc3OWTcR5NfJ2qLg1V173s+TUgHeeZ9pep6n4RGY/z2OIEoAT4pap+UtPYQyl/Hkungn0A7G6ZFovVGGNMoxTqqrCNwB9EZKuIPCwiw6rTsDum2BPAJGAwcIWIDK5UbBLQ331Nx3nkcbi6M4H5qtofmO9OA+wFzlfVE4BrqYPHFKcVHeJAQkvKPJHfQT/9jD4xjMgYY+pf0D2iqj4GPCYivYBpODdGJgGvAnNUdUOYtkcBm9yT/YjIHGAqR7vYcKdfVOcQYIGIpLqXNqeHqDsVGOvWfwH4DLhTVZcFtLsGSBKRxFg8O6b85H1aUR77klpHXG/rrCnRDsUYYxqcsPexqOo2VX1YVYcDVwIXAd9H0HY3YEfAdKY7L5Iyoep2UtUsN7YsoGMV674YWFZVUhGR6SKyWEQW5+TkRPA1gmtXdJDcaiQWY4xpDiIZ0iVeRM4XkZeB/wAbiOxy46rGN6l8dVmwMpHUrXqlIkOAh4GfVbVcVZ9W1QxVzejQoXajC6cV5bHfEosxxlQQ6nLj8cAVwBRgITAHmK6q+RG2nQkEjsrYHdgVYZmEEHWzRaSLqma53WZ7AmLuDvwT+LGq/hBhnNVW3hWWUlpAXkKLWK3GGGMapVBHLHcD3wLHqer5qvpyNZIKwCKgv4j0FpEEnPM0cyuVmQv8WBxjgINu91aounNxTs7jvr8DICKpwHvAXar6dTXirBGP30eSr5SCuMRYr8oYYxqVUCfvx9WmYVUtE5EZOMPBeIHnVHWNiNzkLp8NvI9zqfEmnMuNrw9V1216Fs6Dx24EtgOXuvNnAP2Ae0XkXnfeBFU9ckQTLYrSosw5fZMfnxzt5o0xplGL6ZOmVPV9nOQROG92wGcFbom0rjs/Fzi7ivkPAQ/VMuSIqEKL0iKAiI5YkuI9XDW6V6zDMsaYBsEeul5DLcrKE0tS2LLjBnbk3vMq38JjjDFNkyWWGlCgpXvEkh8fPrHY87+MMc2JJZYaUNWjRywRJRbLLMaY5sMSSw2Vn7wvjOAci6UVY0xzYomlBhRI9JUCUOyND1vejliMMc2JJZYaincTS4knfGIZZo8hNsY0IzG93LipUoV4vw+AEm/oTfjx7WfQt0NKXYRljDENgiWWGlES/M4RS2mYIfP7dWxVFwEZY0yDYV1hNZTgKwPCJxZjjGluLLHUgNMVVkaJJ85uUjHGmEossdRQgq/MjlaMMaYKllhqQIF4f2nYE/en9mtXNwEZY0wDYj+5a6BCV1gQL904itP71+5BYsYY0xjZEUsNhesK89i5F2NMM2WJpQYUJcFfRmmIu+4trRhjmitLLDUQSVeYDeNijGmuYppYRGSiiKwXkU0iMrOK5SIij7vLV4rIiHB1RSRNROaJyEb3va07v52IfCoih0XkL7H8XuAM6VIa5OT91GFdOalX21iHYIwxDVLMEouIeIEngEnAYOAKEan8tKtJQH/3NR14KoK6M4H5qtofmO9OAxQB9wJ3xOo7lVPF6QoLcsTy2LThJMTZwaAxpnmK5d5vFLBJVTeragkwB5haqcxU4EV1LABSRaRLmLpTgRfczy8AFwKoar6qfoWTYGIuzu+j1OM9Zv71p6bXxeqNMabBiuXlxt2AHQHTmcDoCMp0C1O3k6pmAahqloh0rE5QIjId5+iInj17VqfqEYriVT8+qZhYts6aUqP2jDGmKYnlEUtVZ681wjKR1K0RVX1aVTNUNaNDh5rdZ6IKHvXjE+vuMsaYymK5Z8wEegRMdwd2RVgmVN1st7sM931PFGOOmNfvx+exxGKMMZXFcs+4COgvIr1FJAGYBsytVGYu8GP36rAxwEG3mytU3bnAte7na4F3YvgdgvKqH78dsRhjzDFido5FVctEZAbwIeAFnlPVNSJyk7t8NvA+MBnYBBQA14eq6zY9C3hdRG4EtgOXlq9TRLYCrYEEEbkQmKCqa6P/3cCrPusKM8aYKsR0rDBVfR8neQTOmx3wWYFbIq3rzs8Fzg5SJ70W4VZLVSfvjTHG2J33NaKonWMxxpggbM9YQ95KV4WdP7RrPUZjjDENhyWWGqjqcuM4j40NZowxYImlRlKS4o45YrHEYowxDkssNdC3Qwpevw9/wDmWOK8lFmOMAUssNVb5qrCkeLtCzBhjwB5NXCOqSlxAV9j0M/rw32f1q+eojDGmYbDEUhN+P8CRxHL35OPqMxpjjGlQrCusBtTnA7D7WIwxpgq2Z6yJsjIAG9LFGGOqYHvGGjhyxGJDuhhjzDEssdSAukcsNrqxMcYcy/aMNWHnWIwxJijbM9bA0a4w23zGGFOZ7Rlrwk7eG2NMULZnrAE7eW+MMcFZYqkBLbNzLMYYE0xM94wiMlFE1ovIJhGZWcVyEZHH3eUrRWREuLoikiYi80Rko/veNmDZXW759SJybsy+mM+6wowxJpiY7RlFxAs8AUwCBgNXiMjgSsUmAf3d13TgqQjqzgTmq2p/YL47jbt8GjAEmAg86bYTdeVdYXa5sTHGHCuWe8ZRwCZV3ayqJcAcYGqlMlOBF9WxAEgVkS5h6k4FXnA/vwBcGDB/jqoWq+oWYJPbTtQd6QqzxGKMMceI5Z6xG7AjYDrTnRdJmVB1O6lqFoD73rEa60NEpovIYhFZnJOTU60vVM6TlMiatHQOJbSsUX1jjGnKYjm6cVVPvtIIy0RStybrQ1WfBp4GyMjICNdmlRL79eOOM2bUpKoxxjR5sTxiyQR6BEx3B3ZFWCZU3Wy3uwz3fU811meMMSbGYplYFgH9RaS3iCTgnFifW6nMXODH7tVhY4CDbvdWqLpzgWvdz9cC7wTMnyYiiSLSG+eCgIWx+nLGGGOqFrOuMFUtE5EZwIeAF3hOVdeIyE3u8tnA+8BknBPtBcD1oeq6Tc8CXheRG4HtwKVunTUi8jqwFigDblFVX6y+nzHGmKrF9AmSqvo+TvIInDc74LMCt0Ra152fC5wdpM5vgN/UImRjjDG1ZNfLGmOMiSpLLMYYY6LKEosxxpiossRijDEmqiyx1NCp/drVdwjGGNMgWWKpoeevj8kwZMYY0+hZYqmheK9tOmOMqUpM72Np6v724wycW3GMMcaUs8RSC+MHd6rvEIwxpsGx/hxjjDFRZYnFGGNMVFliMcYYE1WWWIwxxkSVJRZjjDFRZYnFGGNMVFliMcYYE1WWWIwxxkSVNOc7x0UkB9hWiybaA3ujFE40WVzVY3FVj8VVPU0xrl6q2iHYwmadWGpLRBarakZ9x1GZxVU9Flf1WFzV0xzjsq4wY4wxUWWJxRhjTFRZYqmdp+s7gCAsruqxuKrH4qqeZheXnWMxxhgTVXbEYowxJqossRhjjIkqSyw1ICITRWS9iGwSkZl1sL4eIvKpiHwvImtE5P+58+8XkZ0istx9TQ6oc5cb33oROTdg/kkisspd9riISC1j2+q2t1xEFrvz0kRknohsdN/b1mVcIjIwYJssF5FDInJbfWwvEXlORPaIyOqAeVHbPiKSKCKvufO/E5H0WsT1exFZJyIrReSfIpLqzk8XkcKA7Ta7juOK2r9blON6LSCmrSKyvB62V7B9Q/3+jamqvarxArzAD0AfIAFYAQyO8Tq7ACPcz62ADcBg4H7gjirKD3bjSgR6u/F63WULgZMBAf4DTKplbFuB9pXmPQLMdD/PBB6u67gq/XvtBnrVx/YCzgBGAKtjsX2A/wJmu5+nAa/VIq4JQJz7+eGAuNIDy1Vqpy7iitq/WzTjqrT8D8Cv62F7Bds31OvfmB2xVN8oYJOqblbVEmAOMDWWK1TVLFVd6n7OA74HuoWoMhWYo6rFqroF2ASMEpEuQGtV/Vadv5IXgQtjEPJU4AX38wsB66iPuM4GflDVUCMsxCwuVf0C2FfF+qK1fQLbehM4O5KjqqriUtWPVLXMnVwAdA/VRl3FFUK9bq9ybv3LgFdDtRGjuILtG+r1b8wSS/V1A3YETGcSeicfVe5h6HDgO3fWDLfr4rmAw91gMXZzP1eeXxsKfCQiS0Rkujuvk6pmgfOHD3Ssh7jKTaPif/j63l4Q3e1zpI6bFA4C7aIQ4w04v1rL9RaRZSLyuYicHrDuuoorWv9usdhepwPZqroxYF6db69K+4Z6/RuzxFJ9VWXqOrlmW0RSgLeA21T1EPAU0BcYBmThHI6HijEWsZ+qqiOAScAtInJGiLJ1GRcikgBcALzhzmoI2yuUmsQR9RhF5B6gDHjZnZUF9FTV4cDtwCsi0roO44rmv1ss/k2voOKPlzrfXlXsG4IWDbKeqMZmiaX6MoEeAdPdgV2xXqmIxOP84bysqm8DqGq2qvpU1Q/8DaebLlSMmVTs3qh17Kq6y33fA/zTjSHbPbQuP/zfU9dxuSYBS1U1242x3reXK5rb50gdEYkD2hB5V9IxRORa4DzgKrdLBLfbJNf9vASnX35AXcUV5X+3aG+vOOBHwGsB8dbp9qpq30A9/41ZYqm+RUB/Eent/iKeBsyN5Qrd/sxnge9V9dGA+V0Cil0ElF+xMheY5l7N0RvoDyx0D4nzRGSM2+aPgXdqEVdLEWlV/hnn5O9qd/3XusWuDVhHncQVoMIvyfreXgGiuX0C27oE+KQ8IVSXiEwE7gQuUNWCgPkdRMTrfu7jxrW5DuOK5r9b1OJynQOsU9Uj3Uh1ub2C7Ruo77+xcGf37VXllRiTca6++AG4pw7WdxrOoedKYLn7mgy8BKxy588FugTUuceNbz0BVzIBGTj/MX8A/oI7+kIN4+qDc4XJCmBN+bbA6X+dD2x039PqMi63vRZALtAmYF6dby+cxJYFlOL88rsxmtsHSMLp6tuEc1VPn1rEtQmnL738b6z8SqCL3X/fFcBS4Pw6jitq/27RjMud/zxwU6Wydbm9gu0b6vVvzIZ0McYYE1XWFWaMMSaqLLEYY4yJKkssxhhjosoSizHGmKiyxGKMMSaqLLGYJkdEOovIHBH5QUTWisj7IjIgButJl4DRbkOUuTJgOkNEHo/Cur8TZ+Tc7SKSI0dH0j1FRN6sbftB1nmhiPza/dzBjWFZwJAllcv/n4icFYtYTMNmlxubJsW9uesb4AVVne3OGwa0UtUvI6jvVVVfsOlKZdOBf6vq8SHaG4szMu951fgaEROR64AMVZ0Ri/YrresbnJsn94rINJx7IK4NUb4X8DdVnRDr2EzDYkcspqkZB5SWJxUAVV2uql+K4/cislqc505cDs7OX5xnWrwCrKpi2uvWWyTOQIg/q7xS98jkSxFZ6r5OcRfNAk53jyZ+7rb9b7dOmoj8y21zgYic6M6/X5zBFj8Tkc0icmukXz7wKEpErnPbf1dEtojIDBG53T3KWCAiaW65viLygTgDiX4pIoOqaHcAUOwmlWE4w7JPdr9XSxF5PmC7/tzd7tuAdiLSOdL4TdMQV98BGBNlxwNLgiz7Ec5AhkOB9sAiEfnCXTYKOF5Vt7hHGYHT04GDqjpSRBKBr0XkIyoOxLcHGK+qRSLSH+dO7QycZ2EcOWJx2y73ALBMVS90u4xedOMDGISTJFsB60XkKVUtreH2GI5z9/Qm4E5VHS4if8QZtuNPwNM4d49vFJHRwJNA5S6sU3HuIkdVl7tdYhmqOkNETgK6lR+5ifuAMNdSt+5bNYjdNFKWWExzchrwqtu1lS0inwMjgUM44yVtCSgbOD0BOFFELnGn2+CMsbQhoHw88Bf317wPZ9DBSOK5GEBVPxGRdiLSxl32nqoWA8UisgfoRMVhzSP1qTrP6cgTkYPAu+78Ve53SgFOAd6Qo4/YSKyinS5ATpB1bAb6iMifgfeAjwKW7QG61iBu04hZYjFNzRqcgfKqEurhRPkhpgX4b1X9sEJjFR/R+nMgG+doyAMURRBrqOHIiwPm+aj5/9XAdvwB0363TQ9wQFWHhWmnECehHkNV94vIUOBc4Bach17d4C5OcuuaZsTOsZim5hMgUUR+Wj5DREaKyJnAF8Dl7jmTDjiPm10YQZsfAjeLMzw5IjJAnNGcA7UBstQZ2v0anEciA+ThdGdV5QvgKrfNscBeDf0sjahz17dFRC514xA3SVT2PdCvqjZEpD3gUdW3gHtxHuFbbgBHRyM2zYQlFtOkqHOZ40XAeHEuN16D88z0XTjPi1mJM+rsJ8CvVHV3BM0+A6wFlronxv/KsUcQTwLXisgCnJ1p+RHPSqBMRFaUn9QOcD+QISIrcU7yB73CKsauAm4UkfJRqqdWUeYLYLhIlY+k7QZ8JiLLcUb7vQuOPCekH7A4BjGbBswuNzbGREREHgPeVdWPIyx/ETBCVe+NbWSmobEjFmNMpH6L85ybSMVx9DHCphmxIxZjjDFRZUcsxhhjosoSizHGmKiyxGKMMSaqLLEYY4yJKkssxhhjour/A/+LPzv0at6qAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "fig, ax = plt.subplots()\n", "ax.set_title(f\"Limiting value: {popt[0]:.5f} Pa s\")\n", "ax.plot(x, viscosity_integral, color=\"tab:blue\")\n", "ax.plot(x_fit, fit_viscosity_integral, color=\"tab:red\")\n", "ax.set_xlabel(\"Correlation Time (fs)\")\n", "ax.set_ylabel(\"Viscosity (Pa s)\");" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.12" } }, "nbformat": 4, "nbformat_minor": 5 }