{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "similar-reply",
   "metadata": {},
   "source": [
    "# Quick Start - From Python\n",
    "Here we present a broad overview of using the PyBigDFT library to drive BigDFT calculations using Python. If you have installed from source, you should make sure you have setup the proper environment variables using the following command:\n",
    "\n",
    "```\n",
    "source install/bin/bigdftvars.sh\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c59d05",
   "metadata": {},
   "source": [
    "## System Manipulation\n",
    "Here we define a system which is compsed of two fragments: H2 and Helium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ede8e9ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "jewish-country",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.Systems import System\n",
    "from BigDFT.Fragments import Fragment\n",
    "from BigDFT.Atoms import Atom\n",
    "from BigDFT.Visualization import InlineVisualizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bored-scottish",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H2:1 H [0.0, 0.0, 0.0]\n",
      "H2:1 H [0.0, 0.0, 1.4]\n",
      "He:2 He [10.0, 0.0, 0.0]\n"
     ]
    }
   ],
   "source": [
    "# Create Three Atoms\n",
    "at1 = Atom({\"r\": [0, 0, 0], \"sym\": \"H\", \"units\": \"bohr\"})\n",
    "at2 = Atom({\"r\": [0, 0, 1.4], \"sym\": \"H\", \"units\": \"bohr\"})\n",
    "at3 = Atom({\"r\": [10, 0, 0], \"sym\": \"He\", \"units\": \"bohr\"})\n",
    "\n",
    "# Construct a System from Two Fragments (H2, He)\n",
    "sys = System()\n",
    "sys[\"H2:1\"] = Fragment([at1, at2])\n",
    "sys[\"He:2\"] = Fragment([at3])\n",
    "\n",
    "# Iterate Over The System\n",
    "for fragid, frag in sys.items():\n",
    "    for at in frag:\n",
    "        print(fragid, at.sym, at.get_position())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "colonial-bottom",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_1684308313740075\"  style=\"position: relative; width: 400px; height: 300px\">\n        <p id=\"3dmolwarning_1684308313740075\" style=\"background-color:#ffcccc;color:black\">You appear to be running in JupyterLab (or JavaScript failed to load for some other reason).  You need to install the 3dmol extension: <br>\n        <tt>jupyter labextension install jupyterlab_3dmol</tt></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://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.0.1/3Dmol-min.js');\n}\n\nvar viewer_1684308313740075 = null;\nvar warn = document.getElementById(\"3dmolwarning_1684308313740075\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nvar viewergrid_1684308313740075 = null;\nviewergrid_1684308313740075 = $3Dmol.createViewerGrid(document.getElementById(\"3dmolviewer_1684308313740075\"),{rows: 1, cols: 1, control_all: false},{backgroundColor:\"white\"});\nviewer_1684308313740075 = viewergrid_1684308313740075[0][0];\nviewergrid_1684308313740075[0][0].zoomTo();\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 1},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 2},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 3},{\"sphere\": {\"scale\": 0.2, \"color\": \"#FF0000\"}});\n\tviewergrid_1684308313740075[0][0].addUnitCell({\"model\": -1},{\"box\": {\"color\": \"black\"}, \"alabel\": \"\", \"blabel\": \"\", \"clabel\": \"\"});\n\tviewergrid_1684308313740075[0][0].zoomTo();\nviewergrid_1684308313740075[0][0].render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_1684308313740075\"  style=\"position: relative; width: 400px; height: 300px\">\n",
       "        <p id=\"3dmolwarning_1684308313740075\" style=\"background-color:#ffcccc;color:black\">You appear to be running in JupyterLab (or JavaScript failed to load for some other reason).  You need to install the 3dmol extension: <br>\n",
       "        <tt>jupyter labextension install jupyterlab_3dmol</tt></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://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.0.1/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_1684308313740075 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_1684308313740075\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "var viewergrid_1684308313740075 = null;\n",
       "viewergrid_1684308313740075 = $3Dmol.createViewerGrid(document.getElementById(\"3dmolviewer_1684308313740075\"),{rows: 1, cols: 1, control_all: false},{backgroundColor:\"white\"});\n",
       "viewer_1684308313740075 = viewergrid_1684308313740075[0][0];\n",
       "viewergrid_1684308313740075[0][0].zoomTo();\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n",
       "\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 1},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n",
       "\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 2},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n",
       "\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n",
       "\tviewergrid_1684308313740075[0][0].addModelsAsFrames(\"MODEL 0\\nHETATM    1 H    H2  A   1       0.000   0.000   0.000  1.00  0.00       B   H  \\nHETATM    2 H    H2  A   1       0.000   0.000   0.741  1.00  0.00       B   H  \\nHETATM    3 He   He  A   2       5.292   0.000   0.000  1.00  0.00       B  He  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308313740075[0][0].setStyle({\"model\": -1, \"serial\": 3},{\"sphere\": {\"scale\": 0.2, \"color\": \"#FF0000\"}});\n",
       "\tviewergrid_1684308313740075[0][0].addUnitCell({\"model\": -1},{\"box\": {\"color\": \"black\"}, \"alabel\": \"\", \"blabel\": \"\", \"clabel\": \"\"});\n",
       "\tviewergrid_1684308313740075[0][0].zoomTo();\n",
       "viewergrid_1684308313740075[0][0].render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "viz = InlineVisualizer(400, 300)\n",
    "viz.display_system(sys)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caring-timeline",
   "metadata": {},
   "source": [
    "## Calculation\n",
    "Calculate the created system using a grid spacing of $0.4$ and the PBE functional. A logfile is generated from which we can access the computed properties. This logfile has built in properties and can be accessed like a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "acquired-headline",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.Inputfiles import Inputfile\n",
    "inp = Inputfile()\n",
    "inp.set_hgrid(0.4)\n",
    "inp.set_xc(\"PBE\")\n",
    "inp[\"perf\"] = {\"calculate_forces\": False,\n",
    "               \"multipole_preserving\": True}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "decent-irish",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.Calculators import SystemCalculator\n",
    "calc = SystemCalculator(skip=True, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "responsible-smoke",
   "metadata": {},
   "outputs": [],
   "source": [
    "log = calc.run(sys=sys, input=inp, name=\"quick\", run_dir=\"scratch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "modified-romance",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-4.053822554528733\n",
      "{'Peak Value (MB)': 123.207, 'for the array': 'f_i', 'in the routine': 'vxcpostprocessing', 'Memory Peak of process': 'unknown'}\n"
     ]
    }
   ],
   "source": [
    "print(log.energy)\n",
    "print(log.log[\"Memory Consumption Report\"]\n",
    "             [\"Memory occupation\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "worldwide-situation",
   "metadata": {},
   "source": [
    "## Periodic Systems\n",
    "We setup a BCC unit cell of iron and perform the calculation using a 2x2x2 k-point grid with a Monkhorst-Pack grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "champion-parliament",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.UnitCells import UnitCell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "official-yellow",
   "metadata": {},
   "outputs": [],
   "source": [
    "pat = Atom({\"Fe\": [0, 0, 0], \"units\": \"angstroem\"})\n",
    "psys = System({\"CEL:0\": Fragment([pat])})\n",
    "psys.cell = UnitCell([2.867, 2.867, 2.867], units=\"angstroem\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "documented-recorder",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_1684308322316214\"  style=\"position: relative; width: 400px; height: 300px\">\n        <p id=\"3dmolwarning_1684308322316214\" style=\"background-color:#ffcccc;color:black\">You appear to be running in JupyterLab (or JavaScript failed to load for some other reason).  You need to install the 3dmol extension: <br>\n        <tt>jupyter labextension install jupyterlab_3dmol</tt></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://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.0.1/3Dmol-min.js');\n}\n\nvar viewer_1684308322316214 = null;\nvar warn = document.getElementById(\"3dmolwarning_1684308322316214\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nvar viewergrid_1684308322316214 = null;\nviewergrid_1684308322316214 = $3Dmol.createViewerGrid(document.getElementById(\"3dmolviewer_1684308322316214\"),{rows: 1, cols: 1, control_all: false},{backgroundColor:\"white\"});\nviewer_1684308322316214 = viewergrid_1684308322316214[0][0];\nviewergrid_1684308322316214[0][0].zoomTo();\tviewergrid_1684308322316214[0][0].addModelsAsFrames(\"MODEL 0\\nCRYST1    2.867    2.867    2.867     90     90     90                          \\nHETATM    1 Fe   CEL A   0       0.000   0.000   0.000  1.00  0.00       B  Fe  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308322316214[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n\tviewergrid_1684308322316214[0][0].addModelsAsFrames(\"MODEL 0\\nCRYST1    2.867    2.867    2.867     90     90     90                          \\nHETATM    1 Fe   CEL A   0       0.000   0.000   0.000  1.00  0.00       B  Fe  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n\tviewergrid_1684308322316214[0][0].setStyle({\"model\": -1, \"serial\": 1},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n\tviewergrid_1684308322316214[0][0].addUnitCell({\"model\": -1},{\"box\": {\"color\": \"black\"}, \"alabel\": \"\", \"blabel\": \"\", \"clabel\": \"\"});\n\tviewergrid_1684308322316214[0][0].zoomTo();\nviewergrid_1684308322316214[0][0].render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_1684308322316214\"  style=\"position: relative; width: 400px; height: 300px\">\n",
       "        <p id=\"3dmolwarning_1684308322316214\" style=\"background-color:#ffcccc;color:black\">You appear to be running in JupyterLab (or JavaScript failed to load for some other reason).  You need to install the 3dmol extension: <br>\n",
       "        <tt>jupyter labextension install jupyterlab_3dmol</tt></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://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.0.1/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_1684308322316214 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_1684308322316214\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "var viewergrid_1684308322316214 = null;\n",
       "viewergrid_1684308322316214 = $3Dmol.createViewerGrid(document.getElementById(\"3dmolviewer_1684308322316214\"),{rows: 1, cols: 1, control_all: false},{backgroundColor:\"white\"});\n",
       "viewer_1684308322316214 = viewergrid_1684308322316214[0][0];\n",
       "viewergrid_1684308322316214[0][0].zoomTo();\tviewergrid_1684308322316214[0][0].addModelsAsFrames(\"MODEL 0\\nCRYST1    2.867    2.867    2.867     90     90     90                          \\nHETATM    1 Fe   CEL A   0       0.000   0.000   0.000  1.00  0.00       B  Fe  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308322316214[0][0].setStyle({\"model\": -1},{\"line\": {\"color\": \"black\"}});\n",
       "\tviewergrid_1684308322316214[0][0].addModelsAsFrames(\"MODEL 0\\nCRYST1    2.867    2.867    2.867     90     90     90                          \\nHETATM    1 Fe   CEL A   0       0.000   0.000   0.000  1.00  0.00       B  Fe  \\nENDMDL\\n\",\"pdb\",{\"keepH\": \"true\"});\n",
       "\tviewergrid_1684308322316214[0][0].setStyle({\"model\": -1, \"serial\": 1},{\"sphere\": {\"scale\": 0.2, \"color\": \"#7F00FF\"}});\n",
       "\tviewergrid_1684308322316214[0][0].addUnitCell({\"model\": -1},{\"box\": {\"color\": \"black\"}, \"alabel\": \"\", \"blabel\": \"\", \"clabel\": \"\"});\n",
       "\tviewergrid_1684308322316214[0][0].zoomTo();\n",
       "viewergrid_1684308322316214[0][0].render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "viz = InlineVisualizer(400, 300)\n",
    "viz.display_system(psys)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "tracked-honey",
   "metadata": {},
   "outputs": [],
   "source": [
    "inp = Inputfile()\n",
    "inp.set_hgrid(0.3)\n",
    "inp.set_xc(\"LDA\")\n",
    "inp[\"kpt\"] = {\"method\": \"mpgrid\", \"ngkpt\": [2, 2, 2]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "buried-tolerance",
   "metadata": {},
   "outputs": [],
   "source": [
    "log = calc.run(sys=psys, input=inp, name=\"psys\", run_dir=\"scratch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "sunrise-interim",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAG6CAYAAAAPuZLqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABuxUlEQVR4nO3deXxU5dUH8N+dPTPZ94SEJCyyKrsIiqBUFLWKdWmrL6K1KorayqtY9G3V2harVq0rWhFU6o6iFlTQAiqLEgRRdiEkQBKyZ5JJMut9/5i5k0SSkGVm7jK/7+eTj2FyZ+aZGJiTc85zHkEURRFEREREUUYn9wKIiIiI5MAgiIiIiKISgyAiIiKKSgyCiIiIKCoxCCIiIqKoxCCIiIiIohKDICIiIopKBrkXoFQ+nw+lpaWIi4uDIAhyL4eIiIi6QRRFNDQ0IDs7Gzpd17keBkGdKC0tRW5urtzLICIiol44cuQIcnJyuryGQVAn4uLiAPi/ifHx8TKvhoiISD4OhwPZ2dkA/EkCm80m84o6Z7fbkZubG3wf7wqDoE5IJbD4+HgGQUREFNX0en3w8/j4eEUHQZLutLKwMZqIiIiiEoMgIiIiikoshxEREVGXDAYD5syZE/xcK7TzSoiIiCgszGYzli1bJvcyQo7lMCIiIopKzAQRERFRl0RRRFNTEwDAarVqZogwM0FERETUpaamJsTGxiI2NjYYDGkBgyAiIiKKSgyCiIiIKCoxCCIiIqKoxCCIiIiIohKDICIiIopKDIKIiIh+wuXxwe31yb0MCjPOCSIiImrjh2P1uPblb9Di9uLZa8binCHpci9Jdnq9HldccUXwc60QRFEU5V6EEtntdiQkJKC+vh7x8fFyL4eIiCJAFEVc/PRX2FVqBwBkJVjw5YJzYNCzcKIWPXn/5v9VIiKigB+O2YMBkCAAZfUt+O/eCplXReHCIIiIiCjgsz3HAQAXnpqJ6ybnAwCDIA1jEERERBTwdVE1AODMQamYFugF+vJAlZxLUgSHwwFBECAIAhwOh9zLCRk2RhMREQFwerzYXlIHAJhYkIKsBAsEAThW14yKhhakx1nkXSCFHDNBREREAA4cb4TT40OS1YiBaTbYzAYMTIsF4N8xRtrDIIiIiAjAvvIGAMCQzDgIggAAOLVfAgDg+6N22dZF4cMgiIiICMD+Cn8QdEpGXPC2Edn+Lda7y5gJ0iIGQURERAD2l58YBA1K95fDDlVqpxmYWjEIIiIiArD/eCOA9kGQ1BNUXN0Er4+zhbWGu8OIiCjqNbk8OFbXDAAYHMj+AEB2YgxMBh1cHh+O1jYhL8Um1xJlpdfrceGFFwY/1woGQUREFPWO1foDoDiLAUk2U/B2vU7AgFQb9pY34GBlY9QGQRaLBatWrZJ7GSHHchgREUW9o4EgKCfJesLX8gOBT0l1U0TXROHHIIiIiKLe0Vp/gJObFHPC1/oFbpMCJdIOBkFERBT1usoE5QSCIKlnKBo5HA7YbDbYbDYem0FERKQlrUFQB5mgRAZBANDUpL1yIDNBREQU9Y4EymEdBUFSdojlMO1hEERERFGvq3KY1BNU43ChyeWJ6LoovBgEERFRVHM4PahxuAC0BjxtJcQYEWf2d4+URnlJTGsYBBERUVQrq28BAMSaDUiIMXZ4DXeIaRODICIiimoVDf4gKCPe3Ok1GfGWwLXOiKyJIoO7w4iIKKpV2P2BTXqcpdNr0uPMgWtbIrImpdHpdJg6dWrwc61gEERERFFNygSlMxPUqZiYGKxfv17uZYScdsI5IiKiXpAyQVKg0xEpQDoepZkgrWIQREREUe14g1QO6zwTJJXKojUTpFUMgoiIKKpJfT7p3cgESVmjaONwOJCWloa0tDQem0FERKQVFd3IBLX2BLVAFEUIghCRtSlJVVWV3EsIOWaCiIgoqkmZoK56gtJi/QGS2yuitskdkXVR+DEIIiKiqNXo9MDh8gLoOhNkMuiQbDMBaN1NRurHIIiIiKKWlAWKNRtgM3fdISIFScejtC9IixgEERFR1Kpq9J8ZlhprOum1aVE+MFGLGAQREVHUqnH4szpSqasrqYG+IOmwVVI/7g4jIqKoVR0IaJJtnfcDSaRAKRqDIJ1Oh/Hjxwc/1woGQUREFLVqAuWwlG5kglICJTOphBZNYmJisHXrVrmXEXLaCeeIiIh6KJgJ6kZPUEowE8TGaK1gEERERFFLKm11KxNkY0+Q1jAIIiKiqFXb5A9okqwnD4KSo7gc1tTUhPz8fOTn56OpqUnu5YSMKoKgRYsWYcKECYiLi0N6ejpmzZqFffv2dXmf9evXQxCEEz727t0boVUTEZHSVTd2vxyWGsWZIFEUUVxcjOLiYoiiKPdyQkYVQdCGDRswb948bNmyBWvXroXH48GMGTO6dYjbvn37UFZWFvwYPHhwBFZMRERq0JNymBQoNbu9aHJ5wrouigxV7A775JNP2v156dKlSE9Px7Zt23D22Wd3ed/09HQkJiae9DmcTiecztZmN7vd3qu1EhGROoiiGAyCujMnyGbSw2TQweXxobrRBWuyKt5CqQuqyAT9VH19PQAgOTn5pNeOGTMGWVlZmD59OtatW9fpdYsWLUJCQkLwIzc3N2TrJSIi5Wl0euDy+gC0Nj13RRAEpAaCpeooLIlpkeqCIFEUMX/+fJx11lkYOXJkp9dlZWXhxRdfxIoVK/Dee+9hyJAhmD59Or744osOr1+4cCHq6+uDH0eOHAnXSyAiIgWQskAxRj1iTPpu3UcqiXGbvDaoLpd32223YefOnfjqq6+6vG7IkCEYMmRI8M+TJk3CkSNH8Nhjj3VYQjObzTCbT/6bABERaUN1D0phEiljFI07xLRIVZmg22+/HR9++CHWrVuHnJycHt//jDPOwIEDB8KwMiIiUpvgtOhu7AyTpETp0RmCIGD48OEYPnw4BEGQezkho4pMkCiKuP322/H+++9j/fr1KCgo6NXjbN++HVlZWSFeHRERqVFPmqIlUsBU3Rhd5TCr1Ypdu3bJvYyQU0UQNG/ePLz++uv44IMPEBcXh/LycgBAQkICYmJiAPh7eo4dO4ZXX30VAPDkk08iPz8fI0aMgMvlwvLly7FixQqsWLFCttdBRETK0ZtymHTQKhujtUEVQdDzzz8PAJg2bVq725cuXYrrrrsOAFBWVoaSkpLg11wuF+666y4cO3YMMTExGDFiBFatWoULL7wwUssmIiIFk5qbuzMjSJJkNQIA6pvcYVkTRZYqgqDuTKdctmxZuz8vWLAACxYsCNOKiIhI7Woc/kAmqQdBUGIgCJKO24gWTU1NmDBhAgBg69atsFqtMq8oNFQRBBEREYVafXP3zw2TJAaurWuOrkyQKIrYvXt38HOtUNXuMCIiCp/1+yrwm2Vb8fz6g/D5tPNG15m6QEkrMcbY7ftImaA6lsM0gZkgIiLCzqN1uOGVQnh9Iv67twImgw43nNW7nbhqIWVzEqzdD4KkrFFdkws+nwidTjvbxaMRM0FERITH1uyH1yfCGpic/M/P9qPZ5ZV5VeHVmgnqfjksIZA18olAg5OHqKodgyAioih33N6CLw9UAgBW3TEFOUkxsLd48OmucplXFj6iKMLei0yQxahHjNEfKHKHmPoxCCIiinKrvy+DKAJj+yeiINWGX4zpBwBYu/u4zCsLn2a3N3h4ak96goDo3SGmRQyCiIii3FcHqgAAF4zMBABMHZLuv/3HKng12iAtlcKMeiFYAuyuaNwhJggC8vLykJeXx2MziIhIG7w+Ed8U1QAAJg9MBQCMyklAnMWA+mY3fjhWj1G5iTKuMDzqpVJYjLHHb+pS5qguijJBVqsVhw8flnsZIcdMEBFRFNtTZkeD04M4iwHDsuIBAAa9DuPzkgAAO47Uybi68JEyQQk9LIUBQJKN2+S1gkEQEVEU23m0HgAwOjcR+jbbvU/LSQQAfKfRIEgalJjYg0GJkoTAbjL2BKkfgyAioii2q9QfBA3Pjm93++hACey7o3URXlFktC2H9VRSFA5MbG5uxoQJEzBhwgQ0NzfLvZyQYU8QEVEU21VqBwCMyE5od/vIfv4/H6pyoNnlRUwPm4eVrjfToiVtByZGC5/Ph8LCwuDnWsFMEBFRlPL6ROwtl4Kg9pmgtDgzkm0miCJwsLJRjuWFVW+mRUuk+0TT7jCtYhBERBSliqocaHH7EGPUIz/FdsLXB6XHAgAOVDREemlh15tp0RIpE1QbReUwrWIQREQUpQ4FMjwD023tmqIlgwNB0P7j2ssEBadFx/S8K0QallgfReUwrWIQREQUpQ5VOQAAA1JjO/z6KRlxAIADGgyC6vqwOywpODGamSC1YxBERBSlpExQQeqJpTCgNRP0o4bLYb3qCQqU0Owtbs1O1I4W3B1GRBSliqRMUFrHQVB+IDg6WtsMj9cHg147vzf3ZYu8VA4TRX9ZLcnW82ySGqWmpsq9hJDTzk80ERH1yKFKfxA0MK3jclhmvAUmvQ4en4iy+pZILi3s6vuwRd6o1yHW7M8hRMsOMZvNhsrKSlRWVsJm6zhoViMGQUREUai+yY1qh78vJr+TcphOJyAnOQYAUFLTFLG1hZvH60OD0wOgdz1B/vvxJHktYBBERBSFDlX5+4HS48zBrEZH+idbAWgrCLK3eIKfx1t61xUildHqoyQTpFUMgoiIotDRWv/RB3kp1i6vy9NgECRNeo4zG3rd5yQFQfYoCYKam5sxbdo0TJs2jcdmEBGRuklBUE5S10FQrhaDoD5Mi5bEW6IrCPL5fNiwYUPwc61gJoiIKAodq/MHNf0SY7q8Li8wSbqkWjtBULApug9BEMth2sAgiIgoCrVmgroOgrTYE9SX7fGS+MCk6bb9RaQ+DIKIiKLQsUAQ1O8kQZD09fpmNxxObbzhSz1BfQmCgpkgTo1WNQZBRERRRhTFbvcExZoNiAvsHtPKrCApexOKIMjewiBIzRgEERFFmdomN5rdXgBAVoLlpNdnBq4p10gQ1BAIXKTm5t6IZ0+QJjAIIiKKMkdr/f096XFmWIz6k16fFWieLq3XxtZoe7M/ExTXyxlBQHQGQVarFVZr15lDteEWeSKiKNPdfiBJVrzGMkHOQCaoL43Rlugqh9lsNjgcDrmXEXLMBBERRZnu9gNJshL9QVAZM0FBbIzWBgZBRERRRmpwzu5GPxDQ2jeklcbo0PQE+QOoBqcHPp8YknVR5DEIIiKKMsft/mAmI757QVBmgr9sVlanjSBI2h0W15cgKHBfUUTwMFYta2lpwUUXXYSLLroILS3a+DkA2BNERBR1ehoEZSdoqxwWzATF9P4t0GLUw2zQwenxwd7s7tN2ezXwer1YvXp18HOtYCaIiCjKHG/wB0GZCeZuXS9tkbe3eFQ/MFEUxTY9QX0LXHh0hvoxCCIiiiKiKOK43QkASI/rXiYozmJErEYGJjo9Pri8/gNA4/vQGA1wYKIWMAgiIooidU1uuDz+ICA9vnuZIKA1GySV0tRKClgEAbCZ+hYESVvso+UkeS1iEEREFEXKA0FMss0Es+HkgxIlabH+gKmywRmWdUVKsBRmNkCnE/r0WCyHqR+DICKiKNLTpmiJlDVSexAkNUX3tR8IaC2nSYEVqQ+DICKiKNIaBHW/FAa0yQQ1qjsIkrbH92VatISZIPXjFnkioigiNUVndLMpWpIW5w+CKlTeE9SaCer72198FDVG22w2iKL2hkIyE0REFEWknqCMbk6LlgTLYWrPBAVKV32ZFi1hJkj9GAQREUWRil6Xw/xBk1Z6gvq6PR7g7jAtYBBERBRFpHJYZg8bo4PlMJUHQfaWvp8gL5GySdGQCWppacGVV16JK6+8UlPHZjAIIiKKIuW93R0WCILqmtxwetR7bEJDS99PkJdEUznM6/Xi3XffxbvvvstjM4iISH08Xh+qAj09PQ2CEmKMMOr9c3WqG10hX1ukSKWrUPQESWePSTvOSH0YBBERRYmqRhdEEdDrBKTYTD26r04nIDVW/SUxZoKoLQZBRERRQsoCpdhMvZqWLJXE1NwcHdKeoMBjuDw+tLi1UyKKJqoIghYtWoQJEyYgLi4O6enpmDVrFvbt23fS+23YsAHjxo2DxWLBgAEDsHjx4gislohImYJBUGzPdoZJ0jQQBIUyExRrMkCKJblDTJ1UEQRt2LAB8+bNw5YtW7B27Vp4PB7MmDEDDoej0/sUFRXhwgsvxJQpU7B9+3bce++9uOOOO7BixYoIrpyISDmkXp7U2J6VwiRpcerfJh/KniCdTggevxENAxO1SBUToz/55JN2f166dCnS09Oxbds2nH322R3eZ/Hixejfvz+efPJJAMCwYcNQWFiIxx57DJdffvkJ1zudTjidrX+x7XZ76F4AEZECVDtay2G90bpNXr1bpEOZCQL8fUH1zW72BamUKjJBP1VfXw8ASE5O7vSazZs3Y8aMGe1uO//881FYWAi3+8Qf1kWLFiEhISH4kZubG9pFExHJrNrhzwRFaznM6xPR4Azd2WFAa3O01g9RtVqtaGxsRGNjI6xWq9zLCRnVBUGiKGL+/Pk466yzMHLkyE6vKy8vR0ZGRrvbMjIy4PF4UFVVdcL1CxcuRH19ffDjyJEjIV87EZGcpHJYSi/LYVIGqcahzi3yjc7WQCVUmSBpm7zWM0GCIMBms8Fms0EQet5Ur1SqKIe1ddttt2Hnzp346quvTnrtT/9HSYe/dfQ/0Gw2w2zu3W9HRERqUB1ojE619e7fOrUHQVI/kNmgg9mgD8ljcpu8uqkqE3T77bfjww8/xLp165CTk9PltZmZmSgvL293W0VFBQwGA1JSUsK5TCIiRWoth/UyExS4X5VKD1Ft7QcKTSkMaG2w1vruMKfTieuuuw7XXXddu/5ZtVNFECSKIm677Ta89957+O9//4uCgoKT3mfSpElYu3Ztu9vWrFmD8ePHw2gM3V8AIiK1aC2H9S4TlBzIINlbPHB7fSFbV6TYQ3h4qiR4iKrGd4d5PB688soreOWVV+DxaKf/SRVB0Lx587B8+XK8/vrriIuLQ3l5OcrLy9Hc3By8ZuHChbj22muDf547dy6Ki4sxf/587NmzBy+//DKWLFmCu+66S46XQEQkK1EU2w1L7I3EGGNwLk6tCktiwUxQiJqigehpjNYqVQRBzz//POrr6zFt2jRkZWUFP956663gNWVlZSgpKQn+uaCgAKtXr8b69esxevRoPPTQQ3jqqac63B5PRKR1DpcXTo8/e9PbcphOJyDZJpXE1BcEtc4ICmEmyBIdjdFapYrGaKmhuSvLli074bapU6fi22+/DcOKiIjURWqKjjHqYTX1/p/+ZJsJVY0uVTZHN7SEblCiJFrKYVqlikwQERH1TVUft8dLUgJ9QdLgRTWxh3hQIsAgSO0YBBERRYGaPg5KlCQHgqhqFZbDGkJ4eKpEyiqxHKZODIKIiKJA64ygvmaC1DsrSGpejjOHLhOUEBiWyMZodVJFTxAREfVNX2cESdRcDmtwhiETFHishhY3fD4ROp12pim3ZbVaUVFREfxcKxgEERFFgeD2+CguhwUzQSHdHeYPgnwi0OjyhLTpWkkEQUBaWprcywg5lsOIiKJAcFBiH8thqSouh4Vjd5jFqIfJ4H8r1frUaC1iEEREFAWk8lVqXzNBgSCoWoVBUDh2hwHRMTDR6XRi3rx5mDdvHo/NICIidenrCfKSlGA5TH1vhOHYHQa0DkzU8jZ5j8eD5557Ds899xyPzSAiInWR5gQl93l3WOv5YS6Pus4PC0dPENAaVHGbvPowCCIi0jifT0Rtkz8I6ms5LCHGCH1gB5T0mGrQ4vbCFTj0NdSZoNZyGIMgtWEQRESkcfXNbnh9/uOHkqx9ywTpdAKSrP43fTXtEJNKVYIAxPbh2JCOSI3WUs8RqQeDICIijZOaohNijMGdTH2hxllB0gnysWZDyGf5xMfwEFW1YhBERKRxoTo3TJKswm3yrSfIh36OTzATxCBIdRgEERFpnFS2SrX1rR9IIg1MrFJROawhTNvjgTY9QRreHaZVnBhNRKRxUtkqVJkgaeBirZoyQWEYlCiJj4LG6JiYGBQVFQU/1woGQUREGhfqclhioLlaTbvDpEyQ1L8TSq3lMO02Rut0OuTn58u9jJBjOYyISOOkwYYpoSqHBXaHqSkIkrI0cWHIBLEcpl4MgoiINC7YExSiTFBSsBymnjf9YCYoDD1B0bA7zOVy4e6778bdd98Nl0s9we/JMAgiItI4qScoOUSZoCQVlsOkLE04MkHRsDvM7Xbjsccew2OPPQa3Wzuvk0EQEZHGhercMIkag6Bw9gRJ5TCHywuPV11HiUQ7BkFERBonnfgeunKY1BPkhiiKIXnMcAtnT1DbbfecGq0uDIKIiDTM5fEFe1VC1RgtZYJcHh+aXN6QPGa4tfYEhT4IMuh1sJn0ALRdEtMiBkFERBomlaz0OiFYtukrq0kfPH5DLSWx1p6g8EyG4Q4xdWIQRESkYVWNUlO0KWRnZglC6yGqatkhFs6J0UDrwEQt7xDTIgZBREQaFmyKtoWmH0iitubo4MToEGXDfioaBiZqESdGExFpmLQ9PjU2NP1AEjUFQT6fiEZnZDJBWi2HxcTE4Icffgh+rhUMgoiINCzU2+MlySo6P6zR5YG0iS0cjdFA69Z7rTZG63Q6jBgxQu5lhBzLYUREGhY8NyxEO8MkiYGeoJom5b/pS/1AJr0OZkN43vak4Io9QerCTBARkYYFzw0LUyaoTgXlsIY2O8MEITTN4T+l9d1hLpcLf/vb3wAA9957L0ym0P48yYVBEBGRhtU4wtMYLZ0kX6OCcpjUrByufiCgTU+QRhuj3W43HnzwQQDA3XffrZkgiOUwIiINq5KCoJA3Rvvf9OtUUQ4L37RoiXQwK8th6sIgiIhIw8JVDpNOkldDJiic54ZJtL47TKsYBBERaZi0Oyw1xI3R0hZ5NfQEBadFm8OXCQr2BDETpCoMgoiINKrJ5UGz23+2V8gbo6WeIBUEQeGeFg203R2mzZ4grWIQRESkUVIWyGLUwRo44DNUEgMnybe4fWhxK/sQ1XBPi/Y/tqHdc5E6MAgiItIo6dywFJs55FvD48wGGAJnkSl9anQkMkFSOczlUX5QSK24RZ6ISKOC/UAhLoUB/kNUE60mVDU6UeNwIStBuUcpSH064dwdZjMZoBMAn+h/PosxtJk3uVksFnzzzTfBz7WCQRARkUZJ54Ylh3hGkCTZZkRVo1Px2+QjkQnS6QTEWYyob3bD3uJGerx2AgUA0Ov1mDBhgtzLCDmWw4iINCp4ZEaIZwRJ1DIwUZoTFK5zwyRSSYzN0erBTBARkUYFp0WHoRwGtO4QU/o2+eCcoDBmggBtH6Lqcrnwz3/+EwDwu9/9TjMToxkEERFplDQoMdQzgiRJgR1iNQ5lv+nbIzAxGmjNNGlxh5jb7caCBQsAALfeeqtmgqCIlMMqKytht9sj8VRERBRQHeZMkDQwUS27w8I5MRrgwEQ1ClsQ5HK5sGDBAqSmpiIzMxNJSUkYNGgQFi9eHK6nJCKiNsLdE6SGIMjj9aHJ5d+yHqlMEM8PU49eBUGbN2+GXq9HWloanE7nCV8XRRGXXHIJ/vGPf6CmpgaiKEIURRw6dAjz5s3Dvffe2+eFExFR14LnhoVpd5h0flitgneHNTpbm5TDuTsMaDswkY3RatGrIOjLL7+EKIr49a9/DbP5xN8wXn31VaxZswYAkJ6ejhtvvBF33nkn8vLyIIoiHn30Ufzwww99WzkREXXK5xODjdGpYcsE+TMftQreHWYP7NSyGHUw6sPbAcJymPr06ifiq6++giAIuOSSSzr8+jPPPAMAyMvLw86dO/HCCy/gH//4B77//nuceuqp8Pl8WLZsWa8XTUREXbO3uOHxiQDCNyeoNROk4CAoQtvjAZ4kr0a9CoIOHToEAJg4ceIJX6usrMS2bdsgCAIWLlyI9PT04NdiY2Nx3333QRRFfPXVV71cMhERnYzUDxRvMcBkCE8GJNgTpOBMUCQGJUrYE6Q+vfqpOH78OOLj4xEXF3fC1zZv3hz8/NJLLz3h6zNnzgQA/Pjjj715aiIi6oZgP1CYSmFAaznM4fLC6fHCbFDeURENEdoeD7Qth2mvJ8hisWDdunXBz7WiV78e1NfXw+vt+IC4bdu2AQD69+/fLgskiYuLQ2xsLBoaGrr9fF988QV+/vOfIzs7G4IgYOXKlV1ev379egiCcMLH3r17u/2cRERqFhyUGKZSGODPfATOUEW9Qpuj7ZHMBGn4JHm9Xo9p06Zh2rRp0OuVF+z2Vq+CoISEBDgcjg5n/2zduhUAMGbMmE7vLwhCj76JDocDo0aNCvYadde+fftQVlYW/Bg8eHCP7k9EpFZVYZ4RBPjPywoenaHQvqDgkRkxEegJYjlMdXoVGg8dOhSbNm3CypUrce211wZvb2pqwpdffglBEDBp0qQO79vQ0ICGhgbk5uZ2+/lmzpwZLKP1RHp6OhITE3t8PyIitYtEOQwAEq1G1DhcqFXo1OhIHZkBtGmMbnZDFEUIghD254wUt9uNF198EQBw0003wWgMf1AZCb3KBF1wwQUQRREPPvggjh07Frz9T3/6ExwOBwB0unPsm2++AQAMGTKkN0/dI2PGjEFWVhamT58erGV2xul0wm63t/sgIlKr6kBjdGoYy2FA6/lhSt0hJkdPkE/090lpicvlwm233YbbbrsNLpcy/1/3Rq+CoJtvvhnJyck4fPgwBg0ahEmTJiEvLw9PPPEEBEHAeeed12mQ88EHH0AQBJx++ul9WnhXsrKy8OKLL2LFihV47733MGTIEEyfPh1ffPFFp/dZtGgREhISgh89yVQRESlNtSNSmSBlB0FSk3IkMkFmgw6mwCwilsTUoVc/FampqXjrrbdw2WWXobGxEV9//XXwa9nZ2XjhhRc6vF9zczPeeOMNAMB5553Xm6fuliFDhrQLwiZNmoQjR47gsccew9lnn93hfRYuXIj58+cH/2y32xkIEZFqtR6ZEeZMUOAQ1TqFNkY3OCOXCRIEAfExBlQ1umBvdqNfYkzYn5P6pteh8fTp07Fr1y688MIL2LFjBwDg9NNPx7x585CSktLhfbZt24Zp06bBaDRiypQpvX3qXjnjjDOwfPnyTr9uNps7nH5NRKRGrUdmhPffNWlWUI1CZwVFck4Q4O8LkoIgUr4+/VTk5ubiL3/5S7evP+uss3DWWWf15Sl7bfv27cjKypLluYmIIi3cJ8hLFF8OCwZBkWnk5Q4xdYlMaNxHjY2N7YYrFhUVYceOHUhOTkb//v2xcOFCHDt2DK+++ioA4Mknn0R+fj5GjBgBl8uF5cuXY8WKFVixYoVcL4GIKGLcXl+wPBWuc8Mkii+HNUvHZkQuEwTwEFW1CNtPRXFxMSoqKiAIAtLS0pCXl9frxyosLMQ555wT/LPUuzNnzhwsW7YMZWVlKCkpCX7d5XLhrrvuwrFjxxATE4MRI0Zg1apVuPDCC3v/goiIVEIqTel1AhLDPB8nUeHlsEhngniIqrqENAgqKyvDokWL8Oabb6K6urrd11JSUnD11Vfjnnvu6XFZatq0aRBFsdOv//Qw1gULFmDBggU9eg4iIq2oCvQDJdtM0OnCO6tG6gmqU2g5rHWLfIQyQYHn0Vo5zGw24z//+U/wc60I2U/Fxo0bMWvWLNTU1HQYsFRVVeHpp5/G66+/jpUrV2Ly5MmhemoiImojuDMszDOCgNZyWK0Cy2FOjxdOjw9AZCZGt30erR2dYTAYcNFFF8m9jJALSRBUUVGBSy65BLW1tYiPj8fcuXNx3nnnIScnBwBw9OhRfPbZZ3jhhRdQVVWFSy65BLt37+7wbDEiIuobaWdYWlz4f2OXymH1zW54vD4Y9OE5sb43Gtr05cSaI5MJ0vIhqloUkp+Kf/zjH6itrcXQoUOxdu1a9OvXr93XpWGFt99+O372s59h3759ePzxx/Hwww+H4umJiKiN6ghmgtr2HNU3u8M+nLEnpCAo1myAPsxlQYlWd4e53W78+9//BgBcc8010X1sxk+tWrUKgiDgX//61wkBUFvZ2dn417/+BVEUg7VFIiIKraoInRsGAAa9LtgHo7SSWKT7gQDtniTvcrlw/fXX4/rrr+exGT91+PBh2Gw2nHnmmSe99swzz4TNZkNxcXEonpqIiH5C6gkK9/Z4SZJNmbOCWo/MiFzWgrvD1CUkQZAgCF3u3upIT68nIqLuac0Ehb8cBrQZmKiwbfKyZIICAVcD5wSpQkiCoLy8PDQ1NWHLli0nvXbz5s1wOBzIz88PxVMTEdFPSIenpkYoCEq2KnNgYqSPzABad4dprSdIq0ISBM2cOROiKOKmm25CZWVlp9dVVFTgpptugiAIHFxIRBQm1ZEuh0kDE5VWDmuJ3OGpEqk/qtHpgcfri9jzUu+EJDy+6667sGTJEuzatQvDhg3DLbfcgunTp6Nfv34QBAFHjhzB559/jhdeeAHV1dVITEzEXXfdFYqnJiKiNkRRbN0dFqEgSKnnh0nToqVm5UhoO4+o0ekJfm9ImULyk5GRkYH3338fl112GWpqavC3v/0Nf/vb3064ThRFJCYmYuXKlZwRREQUBvYWD1yBDEQktsgDbc4PcyirBNQgQybIqNfBatKjyeVFfbObQZDChWyq1dSpU7Fz507cfPPNSEpKgiiK7T6SkpJwyy234Pvvv8fZZ58dqqclIqI2pKboOLMBFqM+Is+ZqNBymBw9QUBrc7SWBiaazWa8/fbbePvtt3lsRmdycnLw/PPP4/nnn0dRUREqKioAAOnp6SgoKAjlUxERUQdaS2GRy0Ao9fwwKRMUyS3ygH+bfLm9RVOzggwGA6688kq5lxFyYQuPCwoKGPgQEUWYdGRGpJqiASApUA5T2knyUiYm4pmgGG0eoqpFIfnJqKurw8qVK7FhwwYcPHgQNTU1APwnxw8cOBDTpk3DrFmzEB8fH4qnIyKiTkR6RhDQNhOkrDf9Bqc8maDWcpiyvh994fF48P777wMALrvsMhgMkQ0sw6XPr+Lvf/87Hn74Ydjt9uBt0iBEQRDw1Vdf4ZVXXsHvf/973HvvvdwVRkQURlUR3hkGtAmCmt3w+UToInRO18nI1ROUoMGT5J1OJ6666ioAQGNjI4MgAJg9ezZef/31YNCj1+sxYMAAJCcnQxRF1NbW4tChQ/B6vairq8M999yDXbt2YenSpSFZPBERtdc6KDFyQVBiYFii1yeiocWDBKsyDtdsDYIinAniwETV6PXusMWLF+Pf//43RFHEmDFj8M4776Curg779u3D5s2bsWXLFuzbtw91dXV4++23MWbMGIiiiFdffRUvvfRSKF8DEREFVDVIgxIjVw6zGPWwmvw70ZQyK0gUxWAQkhAT6XJY4BBVDe0O06peBUFutxt//OMfIQgCfv3rX2PLli24/PLLYbPZTrjWZrPhiiuuwJYtW/CrX/0Koijivvvug8fDHw4iolCTIxMEtJbElBIENbm88Pr8VYpIDkv0Px8zQWrRqyDoww8/RHV1NQoKCrBkyRIYjSePso1GI15++WUUFBSgqqoKH330UW+emoiIuhDsCYrQoESJtENMKUGQFIAY9QJiIjQvSaLUCdp0ol4FQevWrYMgCLjttttgsVi6fT+LxYJ58+ZBFEV8/vnnvXlqIiLqQuvuMJkyQQqZGi01JSfEGCEIkW3UTlLogbJ0ol4FQdu3bwcAnHfeeT2+7/nnn9/uMYiIKDScHm+wGTgtysth9U3ybI8HmAlSk14VSktKSiAIAoYPH97j+w4fPhw6nQ4lJSW9eWoiIupEZYM/C2TS6yLeByNlP5Tyxt96eGrkgyDpe1GvoUyQyWQK7uw2mbRzHlqv/pbY7XbExcX1KsUoCALi4+PbzRUiIqK+qwgEQWlx5oiXgFqzH8p445d6guQJgvzfiwanB26vD0Z9yI7plI3RaMR1110n9zJCrlf/ZxobGxETE9PrJzWbzXA4HL2+PxERnajC3hoERVqyTVnnh8m1PR7wB15SDMq+IGXrVRAkDUfsi1A8BhERtapsaAEApMsQBEkDE5Vyfph0ZEV8hKdFA4BeJwR7kZQSFPaVx+PBqlWrsGrVKk2NuNHG3GsiIgqWw9LjIx8EKe38MDkzQYC/L6i+2a2Y8mBfOZ1OXHzxxQB4bAYA4Pjx49Drezd7QRTFiNeriYi0TiqHpcd1f3RJqEjlMOU0RsvXEwQEeqSqmxTz/aCO9ToIYjmLiEhZKhRQDqt1uBXxi65dAZkgQDvlMK3qVRB0//33h3odRETUR0ooh7m8PjS5vLCZ5S2XSOd2yRcEKWu3HHWMQRARkUZIc4LSYiNfDrOa9DAZdHB5fKhtcskeBAW3yMswLBHgwES1UP/wAiIigtcnBo/MkCMTJAhC68BEBRydoYTGaACoU8D3gjrHIIiISAOqHU74REAQIn94qkRJR2e0NkbLk5FKVFijOHVMG3vciIiinLQzLMVmhkGmCcVKCYLcgb4kQAGZoGZtZIJMJhOeeeaZ4OdawSCIiEgDpH4gOXaGSZJsUjlM3iDI3ibwiJOrJyhGWRO0+8poNGLevHlyLyPkWA4jItKA4PZ4GfqBJEo5P0zqB4ozG6DXybNVPzgygLvDFI2ZICIiDWgdlChfEJSskHKYnCfIS5LanKWmhLlJfeX1evHll18CAKZMmdLrYclKwyCIiEgDgjOCZJgWLVFK9kPOE+QlUk+Q2yvC4fIiVuaRAX3V0tKCc845B4D/2AybzSbzikKD5TAiIg1QQjms9fwweTNBrdvj5Qs8Yoz+uUmA/D1S1DkGQUREGlBe7w+CMuPlywRJ54fJfZK8XeZBiUD7uUlKOVSWTsQgiIhIA0oDQVBWQoxsa0hUyJu+3IMSJUoZGUCdYxBERKRyLo8vOC06K1G+TJBS3vTlPkFe0tojxSBIqRgEERGp3HF7C0QRMOl1wR1acpB2RDW5vGhxe2Vbh9wnyEtae6RYDlMqBkFERCpXbveXwjISzNDJNBcHAOItrXN55HzjV0o5TCnlQeqcuvfsERERyhTQDwT4m4ETY4yodrhQ2+RCZoI8pTl7szQnSN63OC2dJG80GvHII48EP9cKBkFERCpXVtcMAMiSKehoK8lm8gdBMu4QU0omSCnDI0PBZDLh7rvvlnsZIcdyGBGRyiklEwS0DgmUc2BisDFaxi3yAJASq4yRAdQ5ZoKIiFSurF45mSAllIAUkwkKNIpXNao/CPJ6vfj2228BAGPHjuWxGUREpAzlwUyQ/EFQsAQkU/bD5xNbgyCrvEFQaqx/end1YHyBmrW0tOD0008HwGMzIu6LL77Az3/+c2RnZ0MQBKxcufKk99mwYQPGjRsHi8WCAQMGYPHixeFfKBGRDJQwKFGSaJO3HGZvcUMUA2uJkW9cANC+HCZKiyJFUUUQ5HA4MGrUKDzzzDPdur6oqAgXXnghpkyZgu3bt+Pee+/FHXfcgRUrVoR5pUREkaWUQYmSZJnPD5OCL5up9ewuuUjlMI9PDO5YI2VRRTls5syZmDlzZrevX7x4Mfr3748nn3wSADBs2DAUFhbisccew+WXXx6mVRIRRZ5SBiVKpAGBNTIFQVLwlaiA74XZoEec2YAGpwdVDqfs5Tk6kSoyQT21efNmzJgxo91t559/PgoLC+F2d5yidTqdsNvt7T6IiJTuSG0TAKBfUoysgxIliTLvDpMGEyYqJOCQSmLVGmiO1iJNBkHl5eXIyMhod1tGRgY8Hg+qqqo6vM+iRYuQkJAQ/MjNzY3EUomI+uRojX9nWE6S/P1AQGsJSK7G6Lpm//MmKSATBAApGmqO1iJNBkGAf3JpW1JT2k9vlyxcuBD19fXBjyNHjoR9jUREfVVS488E5SZbZV6JnxQEyTUbp9ahjJ1hkpTA96Oas4IUSRU9QT2VmZmJ8vLydrdVVFTAYDAgJSWlw/uYzWaYzeZILI+IKGSkclh/hQRBUuaj0elBi9sLizGy82TqAtvjk5QSBGmkHGY0GnH//fcHP9cKTQZBkyZNwkcffdTutjVr1mD8+PGa+p9HRHREygQlKSMIircYYNQLcHtFVDtc6JcY2TKd1BitmHKYLVAOc6i7HGYymfDAAw/IvYyQU0U5rLGxETt27MCOHTsA+LfA79ixAyUlJQD8paxrr702eP3cuXNRXFyM+fPnY8+ePXj55ZexZMkS3HXXXXIsn4gobI7U+nuCcpOV0RMkCELrG78MfTBSY7Tc06IlwUwQy2GKpIogqLCwEGPGjMGYMWMAAPPnz8eYMWPwpz/9CQBQVlYWDIgAoKCgAKtXr8b69esxevRoPPTQQ3jqqae4PZ6INKXZ5UVlgz/QUEomCJD3jb9WaZkgjTRG+3w+7Nq1C7t27YLP55N7OSGjinLYtGnTupy2uWzZshNumzp1avCcEyIiLToa6AeKNRsUsyUcaPvGH/kgSDoyQynfj2BjtMp7gpqbmzFy5EgAPDaDiIgUQGqKzk22drrzVQ6tb/yRz37UKmhYIsBymNIxCCIiUqniaqkpWhn9QBI5t4XXOZSWCfJnxWqbXPD6eH6Y0jAIIiJSqUOVDgDAgLRYmVfSnlQOq4pwJsjt9aHB6T+jSyk9QUlWIwQBEMXWLBUpB4MgIiKVOljZCAAYmKas/gy5ZuNI/UCAf6u+Ehj0OiQGdqqpvS9IixgEERGplFIzQamx8kyNlrbHx1sMMOiV8/YmV2aMTk45PyVERNRtjU4Pyu0tABSYCZJpTpCSTpBvKz3O//2oaGiReSX0U8rIFxIRUY8UBbJAKTaT4t70pXJYlcMFURQjtnNNygQp5cgMSTAIsqs3E2Q0GoMDh7V08gKDICIiFTpUJfUDKasUBrRmglweHxqdHsRZIvOmKTUeJygsKMyItwAAKhrUGwSZTCY8+uijci8j5FgOIyJSoQPHA0FQurJKYQAQY9LDavIfnBrJZuB6hR2eKkkLlsPUGwRpFYMgIiIV2lNmBwAMzYyXeSUdax0SGLk3/uCgRIWcGyZJD2SCjtvV2xPk8/lw+PBhHD58mMdmEBGRvHYHgqDh2QoNgmxmHKlpRlUEM0FST5DSeqSknqBKFWeCmpubUVBQAIDHZhARkYxqHS6U1fuzCkMz42ReTcfk2CbfGgQpKxMU7AlScSZIqxgEERGpjFQK659sjVjTcU/JsU1eaSfIS6RMkMPlRWNgojUpA4MgIiKVCZbCspRZCgPabJOPYDlMyjol25QVBNnMBtgCjeLMBikLgyAiIpX57mg9AGCEQvuBgNYpyZE8RLVaoUEQoI1t8lrEIIiISGW2Ha4BAIzLS5J5JZ1LDZ4fFpk3fZ9PDGaCpCyUknCbvDIxCCIiUpHSumaU1rdArxMwun+i3MvplJSNidScIHuLG16f2O65lSSdzdGKxC3yREQy6c2REoXFtQD8/UBWk3L/CQ82RkdoTpBUCoszG2A26CPynD2RrvJMkMFgwK233hr8XCu080qIiFTicJUDD360C18eqEKsxYDLxvTD/PNO6dZOr40HqgAAE/KTw73MPmm7Rd7rE6HXhff8sGBTtAJLYQCQES+dH6bOTJDZbMazzz4r9zJCjuUwIqIIOnC8AbOe24h1+yrh8Ymoa3Jj6cbDuPSZjSita+7yvqIoYt2+CgDAOUPTIrHcXku2mSAIgE9s3boeTlLZTYmlMKC1MVqa70TKwCCIiChCnB4vbv33t6hrcuO0nAR8/LspePm68chOsOBQlQO/enFLl43Eu0rtqGhwwmrS4/QCZWeCDHodUgIBSSROTw82RSs0COqXGANAvUGQKIqorKxEZWUlRFGUezkhwyCIiFRjy6FqPPSf3VjyVREcKhw699KXRThQ0YjUWBOWXX86hmXF49yhGXjnlsnITY5BSU0T5i7fBqfH2+H9P9pZCgCYMjhVkX0vP5UWJ20LD/8bvxQ8Sr1ISpMVDIKa4fOpL4hoampCeno60tPT0dTUJPdyQoZBEBGpwj8/O4BfvbgFS74qwkP/2Y1Ln92IqghOI+6rRqcHL35xCABw30XD2pVt+iXGYOl1ExBnNmDr4VosXPH9Cb9tu70+vP/tMQDAZWNyIrfwPojktvBqpfcExZmhEwC3V1TVz63WMQgiIsX7bPdxPPHZfgDAxadlIT3OjB8rGvG7N7erJjX/+tfFqG92Y0CqDZeM6nfC1welx+HZa8ZCrxPw3vZjePq/P7b7+vvfHkNFgxOpsWacOzQ9Usvuk0geHKr0cphBr0NmoC/o2El6vyhyGAQRkaK1uL344wc/AABunFKAZ64ei9dvnAizQYeNP1bjw+9KZV7hyYmiiNe/LgEA3Dx1QKc7pc4+JQ0PXToSAPD42v1YvqUYgL/U8+iafQCAuVMHwGRQxz/dcgRBSm2MBoDsQEmstE6dfUFapI6/SUQUtd4pPIKy+hb0S4zB/POGAPBnTW6dNggA8Mx/f1R8j8W24locrm6C1aTHxadld3nt1RP74+apAwAA/7fyB1z9ry34xfObUNngxKD0WPzPGXmRWHJItM7GiUBPkAqCoLZ9QaQMDIKISLG8PhEvfunvo7np7AGIMbU2A19/Vj7iLAYcqGjE+v0Vci2xW97ddhQAMHNkFmzmk49n+8MFQzH/vFOg1wnYdLAaxdVNyIg3Y/H/jIXFqPyGaEnrlOQI9AQF+mxSY5XZGA0A2YkshykNhyUSkWJtOVSNIzXNiLcYcNX43HZfi7cYcdX4XCz5qghvbz2Kc4dmyLTKrrm9Pqz6vgwAcPm4E3uBOiIIAu6YPhgXn5aFdfsqYTHq8PNR2YjvxjBFJYnUlGRRFIOziJScCeoXLIcxCFIKBkFEpFjvBXZDXTwqu10WSCIFQZ/tOY7qRmfw5HIl2VpUg4YWD1JsJkwsSOnRfQekxWJAWmyYVhZ+aW3KYb05IqS77C0euL3KPTdMkp2g3p4gg8GAOXPmBD/XCu28EiLSlGaXFx//EMigjO04gzIkMw6n5SRg59F6rP6+DLMn5Udwhd3z+V5/qW7akPSwHx2hNOmBOUEtbh8anJ6wZbKkpmibSa/ocmG2ijNBZrMZy5Ytk3sZIceeICJSpI0/VqHJ5UW/xBiM7Z/U6XUXnpoFAFiz+3ikltZtoiji8z3+df1smDq2tYdSjEmPuEAPVDh3iNUEDmlV6owgidQTVO1wocXd8UBMiiwGQUSkSFIGZfqw9C7LKOcN9/cCbTlUjYYWd0TW1l2Hqhw4XN0Ek16HKaco+6yvcEkLHhwaviBIOjdMqdOiJQkxRlgDZV21NUeLogiHwwGHw6Ga2VzdwSCIiBRHFEX8d68/g3KywYAD02IxINUGt1fEhv2VkVhet23Y51/PxAHJiO3GrjAtisQ2+apAEJSq8EyQIAjon2wFAJRUq+voiaamJsTGxiI2NpbHZhARhdOuUjuO252IMepxxoCTNxNL2aDP9yhrq/zmQ9UAgDMHpcq8EvlIfUHhLIdJAZZ0VpmS5aX4g6DiaofMKyGAQRARKdC6QCnsrMGp3Wp0nRooNW38sUoxqXqfT8Q3RTUA0K1ATqsicX6Y9NjScylZfooNAHBYZZkgrWIQRESKs+mgP4Nydjf7aMbmJcFk0KGiwYmDlY3hXFq37S6zo77ZjVizASOz4+VejmwicXSG9NjpKgiC+gcyQSU1DIKUgEEQESlKi9uLb0tqAQCTuplBsRj1GJ/n30G28cfqsK2tJ7YESmET8pNg0EfvP7Xp8eHvCapQURDUmgliOUwJovdvJhEp0ndH6uD0+JAaa8bANFu37yf13Ww6WBWupfWIFARFcykMaO0JOh7G3WFVKiqHSY3RR2qa4FX4mXfRgEEQESnKlkP+PpqJA5J7NGF48kB/sLH5YLXsby5en4iv2Q8EAMhM8AdB5fXhyQSJothaDotXfmN0dmIMjHoBbq/Ig1QVIDr3bBKRYvU2g3JqvwTEmg2wt3iwr7wBw2Xsw9lTZkdDiwexZgNGRHE/EABkBYKgRqcH9hZ3yKdG1ze74fL6ACh/izwA6HUCcpOsOFTlQHF1E3KSrHIvqVv0ej2uuOKK4OdawUwQESlG+36g5B7d16DXYUz/RADAtuKaUC+tR7YfqQMAjOmfGNX9QABgNRmQEOMPfMKRDZL6gRKtRpgN6nhzlrbJF1Wppy/IYrHgnXfewTvvvAOLRfkZt+6K7r+dRKQoPxyrh9PjQ4rNhIG9ODh0XKA5urC4NtRL65HtgUBuTBfHfUQTKRsUjjOzpFJYmgIPz+2M9LP9Y4UydjJGMwZBRKQY20vqAPi3vPfmxPHxef7sUeFheYOgHYHXIWWmop10cGhZWDJB/seUdqGpwSkZcQAYBCkBgyAiUoztR6QMSmKv7j+6fyJ0gv9cpnA14p5MrcOFQ4Eyx+icRFnWoDRSc3RYgiC7tD1ePSWaQRn+TND+4w0yr6T7HA4HBEGAIAhwONRTxjsZBkFEpBhSJmhMbu/KSLFmA4Zl+RuRC2XqC9pxtA4AUJBqQ5JN+Y26kZAtBUHhLIepYHu8ZHC6PwiqaHCivklZh/5GGwZBRKQIZfXNKKtvgU4ATstJ6PXjSEMT5SqJtQZyibI8vxJlJoSzHKaeQYmSOIsx2Cd1oEI92SAtYhBERIog9dEMyYyHrQ8nro8NBEHbZGqObm2KTpTl+ZVIygSVhmEujhozQQAwONAXdIB9QbJiEEREiiBtKx/bx+BhfL6/OXp3mR1NLk8fV9UzPp+I74Lb47kzTJIVaIwur28J+QG3x+2BxmgV9QQBrSUxNfUFaZFqgqDnnnsOBQUFsFgsGDduHL788stOr12/fn2wgavtx969eyO4YiLqiW+LQ7OtvF9iDDLjLfD6RHx3pD4US+u2Q1UO2Fs8sBh1GJIZF9HnVrLMwCTnJpcX9ubQBaaiKAazS9mJ6gqChgQyQXvK7DKvJLqpIgh666238Pvf/x733Xcftm/fjilTpmDmzJkoKSnp8n779u1DWVlZ8GPw4MERWjER9YTL48P3x/wBSyjKSNK8IGnwYqRIpbDT+iXCGOVDEtuKMemRZPUPTAxlSay+2Y0Wt39adIYKjsxoa2Q/f9/brmN2+HiGmGxU8bf08ccfxw033IDf/va3GDZsGJ588knk5ubi+eef7/J+6enpyMzMDH5oadQ3kZbsLbfD6fEhIcaIgpTuH5raGSmQ2h7pIKjNpGhqLyuhtSQWKqV1/sdKsZlgMarr3/fBGbEwGXRocHpQXNMk93JOSq/X48ILL8SFF16oqfdSxQdBLpcL27Ztw4wZM9rdPmPGDGzatKnL+44ZMwZZWVmYPn061q1b1+W1TqcTdru93QcRRYa0o2p0biJ0up4PSfypscFMUF3Ie1C6sp1DEjuVFYbm6HK7/7GyVFYKAwCjXhcc5yBlQZXMYrFg1apVWLVqFY/NiKSqqip4vV5kZGS0uz0jIwPl5eUd3icrKwsvvvgiVqxYgffeew9DhgzB9OnT8cUXX3T6PIsWLUJCQkLwIzc3N6Svg4g6F+odVSOy42HS61DjcKG4OjK/ZTucHuwr9//yxKboE/VL8meCjtWGLgiSMkGZ8TEhe8xIOrWfPwj6QQVBkFap5hT5n47QF0Wx07H6Q4YMwZAhQ4J/njRpEo4cOYLHHnsMZ599dof3WbhwIebPnx/8s91uZyBEFCHbQ7yjymzQY2S/eHxbUodvS2qRn9r3EtvJ7DxaD5/o3w6utv6USMgNnJZeEsLST5lKm6Ilp/VLBFCCnYEBmxR5is8EpaamQq/Xn5D1qaioOCE71JUzzjgDBw4c6PTrZrMZ8fHx7T6IKPyqG53BbE0oj5kY2z+yzdHSkR+jWQrrUG6yP1tzJISZIGn4onQsh9qcGhgK+sMxO7wKb452OByw2Wyw2Ww8NiOSTCYTxo0bh7Vr17a7fe3atZg8eXK3H2f79u3IysoK9fKIqI++C/wWPCDNhoTADqJQaB2aWBeyx+zKjj4e+aF1ucn+TNCRUGaCAuWw7AR1lsNOyYhDnMWARqcHu0uV34fa1NSEpiblN3H3hCrKYfPnz8fs2bMxfvx4TJo0CS+++CJKSkowd+5cAP5S1rFjx/Dqq68CAJ588knk5+djxIgRcLlcWL58OVasWIEVK1bI+TKIqAM72jRFh5KUCdpXbkej04PYPkyhPhlRFLkz7CSkIKjG4YLD6enTVHBJeWBQYpZKM0F6nYAJ+cn4794KfF1UHcwMUeSoIgj65S9/ierqavz5z39GWVkZRo4cidWrVyMvLw8AUFZW1m5mkMvlwl133YVjx44hJiYGI0aMwKpVq3DhhRfK9RKIqBM7jgbmA4U4CMpMsKBfYgyO1TVj55E6TB6UGtLHb+tYXTMqG5ww6ITg/BdqL95iRKLViLomN47UNmFoZt9aDkRRRGngQNYslWaCAOD0AikIqsFvpwyQezlRRxVBEADceuutuPXWWzv82rJly9r9ecGCBViwYEEEVkVEfSGKrcdMjA5DGWlM/0Qcq2vGtyW1YQ2CpK3xw7PjVTevJpJyk6yoa6pHSXXfg6AahwtOjw+CAGQkqOvcsLYmFviPedl6uAY+nxiSERHUfYrvCSIi7SqqcqC+2Q2TITzHTLQ2R9eF/LHb4snx3RPK5mjpMTLjLTAb1Bt4juyXAKtJj7omN/aUK78vSGsYBBGRbHYEskAjs+NhMoT+nyOpOXp7SW1YhyZKO8M4H6hroWyOlrbaS1vv1cqo12HywBQAwLq9FTKvJvowCCIi2YSzFAYAw7PiYTboUNvkRlFVeLb1Oj1e7DomDUlMDMtzaIUUsIQiCJIeQwqs1Gz6MP+4l88VHATpdDpMnToVU6dOhU6nndBBNT1BRKQ9UiZoVG54molNBh1O7ZeAwuJabCuuxYC02JA/x+5SO1xeH5JtJvTXwBtyOEnfnyO1oQyC1NsULTl3aDoA/9+HqkYnUmOV1+MUExOD9evXy72MkNNOOEdEqtLi9mJ3WSCDEsbZOm3PEQuHtv1AnU2xJ7+8FH8QVFzd1OeT06VymBYCz4x4C07tlwBRBD7fc1zu5UQVBkFEJIs9ZXa4vSKSbaaw/jYvNUeH60R5zgfqvpwkK0x6HZweH47V9a05WktBEACcP8JfEvvwu1KZVxJdGAQRkSx2BPuBwptBGZuXCADYd7wBDS3ukD9+6+GvbIo+Gb1OQH6qP2g5WNnY68dxe33BIzO00BMEAJeO7gcA2HSwOngmmpI4HA6kpaUhLS2Nx2YQEfVVsB8ohOeFdSQ9zoKcpBiIIvDdkdCe1l3R0IKjtc0QBOA0TvvtloGBvqyDlb1/Iy2ra4HXJ8Js0CFNgf0zvZGbbMXp+ckQReCDHcrMBlVVVaGqqkruZYQUgyAikkVwZ1gEykjhOkxVOvLjlPQ4xFlCd+6Zlg1IswEADvUhE1TSZmeYloYL/mKsPxv05jclfe6Zou5hEEREEVfrcOFw4OT4URHIoIwNBFqhDoK2BR4v1OeeaVlrJqj3QVBRtT+LlKeRUpjk56OyEW8x4HB1E/6r4O3yWsIgiIgibkfg5PiCVBsSraawP19wh1hxbUh/w95aVAMAmBA4+oBOLhTlsIMV/gBqYHroRx7IyWY24NcT+wMAXt5YJPNqogODICKKuG2H/RmUsRFqJh6WFQ+LUQd7iweHqnqfgWir2eXF98f8PUan5zMI6i6pHFbZ4IS9l43qUhZpYOCxtOTaSfnQ6wRsOlgdLBlT+DAIIqKI23rYn0EZnx+ZIMio1+G0fokAgG+L60LymNuP1MLtFZEZb9HEwL5IibMYkR7nb2Y+1MtskHS/gWEYfim3fokxmBXYKfbYmn0yr0b7GAQRUUS5PD58FyiHTYhQEAS0HZoYmr6grUX+x5lQkMwhiT0kBS/7jzf0+L5NLk9wxpAWgyAA+P3PBsOoF/DlgSpsPlgt93IA+I/NGD9+PMaPH6+pYzO080qISBV2ldajxe1DotWIAamRexMLdXP0N4f9b06nsx+ox4ZlxQPwD8zsKSkLlGwzIckW/n4yOeQmW/GrCf7eoIc/3qOInWIxMTHYunUrtm7dipgY7WQ+GQQRUURtK/YHIePzkiK6vVnKBB2oaER9c9+GJrq9vmBZjf1APTcsKw5A74IgqR9oQKr2+oHaun36IMSaDfjuaD3eKjwi93I0i0EQEUWU1A80Li+ywUNqrBn5KVaIYuuurt764Vg9mt1eJMQYMVhjO5QiYXi2PxO0u9QOUexZlmNfub+ENjhD29/39DgL7jzvFADA3z/ZixqHS+YVaRODICKKGFEUg5mgSPYDSSYPSgUAfPVj36bebgr0aUwsSNbUsL5IGZweB6NegL3F0+MzxKRDd4dna39C95xJeRiaGYe6Jjce+WSvrGtpampCfn4+8vPz0dTUJOtaQolBEBFFTHF1E6oaXTDpdRjZL/JvYlNCFAR9sb/S/3inpPV5TdHIZNBhULpUEutZc7RUQhseKKlpmUGvw0OzRgIA3tx6BJsOyndkhSiKKC4uRnFxcY+zd0rGIIiIIubrIn8G5bScBFiM+og//+SBqdAJwI8Vjb0+pNLh9ASbq6WginpueKA5eldp989zq2504rjdCUEAhmTGh2tpijIhPxlXBwYo/mHF92hyeWRekbYwCCKiiNn4oz8ImjwwRZbnT7AacWrgwNavDvTut+qvi6rh9orITY5BXoq2jm2IpBGBvqDvj3Y/CJKyRnnJVsSaDWFZlxItnDkU2QkWlNQ04ZFPODsolBgEqZi/v6IGj6/dj/9b+T2e+vxAj36rIookURSDvTSTZcyg9LUk9sV+//2mDE7jfKA+GBfYrbetpPtHmUj/vklb7KNFnMWIRZefBgB4ZfPh4OaCSCita8ZTnx/Aja8WBm9b/X0ZPF5fxNYQTtETSmvM/uMNuPe971FY3H7myeNr92PqKWn4y6yRyNXY4YKkbgcqGlHV6ITFqMOYCJwc35kpg1PxzLofsWF/JTxeHwz67v8uKIoi1u/zH2x59mCWwvpieLb/KJO6JjcOVTUGe4S6sr2kDgAwKgoPrJ16ShquGp+DtwuPYsG7O7H6jimIMYWvpOz0ePHsf3/E8xsOwu0V4XO1BL/2v29/h8Ubj+Lxq0ar/v8FM0Eq9MkP5Zj17EYUFtfCYtThklHZuP3cQZg5MhNGvYAN+yvx82e+wpZDypg0SgQAGwOZlwn5yTAbIt8PJBmXl4Rkmwl1TW5808Ot8vuPN+JwdRNMBh2mDGZTdF8Y9TqMCpQmCw+ffIClKIrBX/rG50V+Z6ES3HfRcGTEm1FU5cA/wnikRnWjE1f/62s89d8f4faKOL0gGX/6+fDg1xOtRhysdOCKxZvw0XelYVtHJDAIUplPfijHvNe/RZPLizMHpWD9XefgqV+Pwf/OGILn/2cc1tw5FaNyE1HX5Ma1L3/DQIgUo7UfSN4MikGvw8+GpQMAPt1V3qP7rglcP2VQKmxR1JMSLsGSWPHJg6AjNc2oanTCqBdk2VmoBAkxRiz6xakAgCUbi7Cpj7scO1LZ4MQVizdjW3Et4iwGPH/NWLx10xm4emIehg8fjuHDh+PT30/FjOEZcHtF3PHmdnyw41jI1xEpDIJU5JuiGtz+xrfw+kT8Ymw/vHL96chMsLS7piDVhrduOgPTh6bD5fHhxlcKsbe851NZiULJ4/UFd4adOUiepui2zh+RCQD4dNfxHh1J8OlufxA0Y0RGWNYVbaQDdLvT41JY7L9mZD95dhYqxblDM/Dr0/tDFIH/fec71Df1bfp5Ww0tbly39BsUVTnQLzEG7996JmaemgVBEGC1WrFr1y7s2rULWakJeP5/xgXXcfc7O3ucVVUKBkEqUWFvwbzXv4XbK2LmyEw8cvlpnfYyWIx6PHvNWJyen4wGpwdzX9sGe0vo/qIQ9dT2I3VoaPEgIcaIEQoYcnfmoFTYTHqU21uw/Uhdt+5TVOXAD8fs0AnA9GEMgkJhQn4yDDoBh6ubUFzd9Yny0pvsuP7RWQpr648XD0NBqg1l9S24b+X3IZnb4/R4MXf5NuwqtSPFZsLy307EoC6moet1Av46ayRmjsyEy+vDza8Vory+pdPrlYpBkAp4vD7c9vp2VDY4MSQjDv+4atRJmzktRj1emD0O/RJjcLi6Cfe8u1NTA65IXT7bcxwAcM6QNOgVMGHZYtRjRiAbtOLbo926z4pt/uvOPiUNqbHmsK0tmsRZjMFs0Pp9lZ1eJ4picEDlWWxIh9VkwBO/HA29TsB/dpbhgx1968vx+UTMf/s7bPyxGjaTHsuuPx0F3TibTacT8PhVozEiOx61TW7c+dYOeBVw2GtPMAhSgZc3FuGbwzWINRvw/P+MhdXUvV6EJJsJz1w9Bka9gI9/KMe727r3jz1RqH2+x7+jSkkZlCvH5wAAPtpRimaXt8trfT4R7wWCpcvH5oR9bdFk2hB/f9a6wK67jvxY0YjS+haYDTqcMUD+cqoSjM5NxO+mDwYA/HHlDzhc1XUmrTOiKOLBj3Zh1c4yGPUCFs8eh1NzTszWNjU1YcSIERgxYkS7YzNiTHo8/esxsJr02HyoGos3HOzdC5IJgyCF8+8C2A/AnwIdkNazQwPH9E/C/POGAAD+/NHuXk/JJeqt4moHfqxohEEn4GwFHTNxRkEKcpNj0OD0YPX3ZV1e+8WBSpTWtyDOYsB5w5UTyGnBOYEgaNPBatQ3d1y2/3yvP0A6Y0BKVPcD/dSt0wZifF4SGpwe3PRaIRqdPZ8m/ey6H/HK5mIAwD+uGt3prkdRFLF7927s3r37hKrCgLRYPHjJCADAk5/tx/7jPTsKRU4MghTM5xNxz4qdcHp8OGtQKq4an9urx7lxSgFG5yaiwenBH1aEpn5M1F2fBbJAE/KTkRBjlHk1rXQ6Ab+a4D+O4F9fHury78WLXxwCAFwxLodvwiF2SkYsTsmIhcvj6zQYlco9bEhvz6DX4blrxiI9zoz9xxvxv2/3rBz12ubDeCzwS/b9Px+OS0Zl93otV4zLwc+GpcPtFbHg3Z2qKYsxCFKwf39djG+KahBj1GPRL07t9XRag16Hx64cBZNBhw37K/FOIctiFDlrAzuqpge2pSvJ/0zMg82kx97yBvx3b8flmB1H6rDpYDUMOgG/nTIgwivUPkEQgiXGFR2U7PeVN2BPmR1GvYCLTs2K9PIULz3egsWzx8Gk1+HTXcdx3/vfd2vH47vbjuJPH+4CANwxfTCuP7OgT+sQBAF/mXUq4swG7DhSh6Ubi/r0eJHCIEihjtY24eGP9wIAFlwwpM/Tnwelx+KuGacAAB5atRvH7err4if1Ka9vwdeBXT0zFfgGlmA14n8m5QEA/rZ6D1ye9kcBiKKIv/xnNwDg0tH90C8xJuJrjAaXjekHvU5AYXEtvvvJbr1lmw4DAM4dmo5Eqynyi1OBsf2T8PgvR0En+E+bv2fFzhN+liWiKOL59Qdx1zvfQRSBOZPycOfPBodkHZkJFtx70TAAwGNr9p10x58SMAhSIFEUce/7P8Dh8mJcXhLmTMoPyeP+5swCjMpJQEOLB/e9z7IYhd9/dpZCFP0TfpUaQNw6dRBSY004WOnA0/890O5rr2w6jMLiWsQY9bjr/FNkWqH2pcdbMGt0PwDAP9buD/7bdKSmKbh774azmIXrysWnZeOxK0dBEIB3th3FlYs3YefRunbXHKpsxPXLtuLvn/h/wb5xSgHu//mIkJ6B96sJuZg0IAUtbh8Wvqf89xkGQQq04ttj+GJ/JUwGHf5++WnQhWhLsUGvwyNXjIJRL+CzPRX4UOXjzkn5pJ+xS0b3vtcg3BKsRvzxYv+RAE//90cs3VgEj9eHFduO4i+r9gDwZ2OzEpQZxGnF7ecOgkEn4Iv9lVj+dQmcHi/ufvc7uDw+TB6Yggn5nA90Mr8Ym4Ol101AnMWA747W45JnNmLmP7/Erf/ehlnPbsT0xzdg/T7/e8tDl47AfRcND9n7i0QQBDx8+amwGHXYdLAab249EtLHDzVBVHqYJhO73Y6EhATU19cjPj5yJxZX2Fvws8c3wN7iwYILhuDWaYNC/hxPfX4Aj6/dj2SbCWvvPBspnHlCYfBjRSN+9vgG6HUCvr53uuJn6/x11W7860t/H4MgANK/jJeN6YfHrxrFE+Mj4Nl1P+LRT/1nYiVZjahtcsNq0uOj28/CwB7ujI1m5fUt+Psne/Hhd6UnNChPH5qOe2YOxSkZJz+wtq2mpiYMH+7/ZWH37t2wWrtu0Xjpy0P4y6o9iDMbsGb+2RH9JaIn798MgjohRxAkiiJufm0b1uw+jpH94rHy1jN7dMJ1d7m9Pvz86a+wt7wBF5+WhWeuHhvy5yB68KNdWLrxMH42LB0vzZkg93JOShRFvLzxMJ76/ADqm/1vvjedPQC3nztYEQMeo4HPJ+Lvn+zFi18egij6A6Fnrx6LyYM4ILE3qhqd+PpQDSoaWpBkNeH0gmRkR6gs7fWJuPz5TdhxpA7Th6bjpTnjI/aLBIOgEJAjCProu1Lc/sZ2GPUCPrztLAzLCt/zfn+0HrOe2wivT8QLs8cFz1IiCoUmlwcT//Y5Glo8WHb9hOBAPDVweXyoaGhBepwFJgM7BuRwpKYJR2ubcWpOAmJ5UK1q7T/egIue+hJur4h//mo0Lg30fYVbT96/+TdcIaobnbg/sF1x3jmDwhoAAcCpOQm4+Wx/o+H/rfwhpIfwEX2woxQNLR70T7bi7E6GrymVyaBDTpKVAZCMcpOtmDQwhQGQyp2SEYfbzvHvPHvgw12obnTKvKIT8W+5AoiifyhijcOFoZlxYekD6sgd0wdjYJoNlQ1OPLRqd0Sek7TP6xPxr8Bwwf85o3/IGy+JKPKam5sxYcIETJgwAc3N3T954JZpAzE0Mw61TW7cq8BdyQyCFGDZpsP4bE8FTAYdHr9qdMR+A7UY9XjkCv+Wyne3HcXHJzk6gKg7/rOzFIeqHEi0GnH1xDy5l0NEIeDz+VBYWIjCwkL4fB3PIOqIyeAf1mvUC/h013G8Epj7pBQMgmS282gdFq32z2z4v4uGYXh25HaiAcC4vCTcFCiL3f3uThysbIzo85O2uDw+/PNz/6yd355VwHIGEWFkvwQsnOkfovi31Xvx/dF6mVfUikGQjMrrW3Djq4VweX2YMTwDs8+Q57fmu2cMwekFyWh0ejD3tW1oaGF/EPXOq5sP41ClAyk2E66dnC/3cohIIa4/Mx8zhmfA5fVh7vJtqGhQxqkFDIJk4nB6cOOrhThud2Jweiwek3EOiUGvwzNXj0F6nBkHKhoxd/k2OD1eWdZC6lVa14x/fubPAi24YAjiLco5LJWI5CUIAh69YhQKUm04VteM375SiCZXz0+9DzUGQTJocnlw/bKt+P5YPZJtJiyZM0H2N4z0OAuWzJkAm0mPjT9W4863dsDt7X7dl6Kb1yfi92/tQIPTg9G5ibhyXK7cSyIihUmwGrH0uglIshqx82g9bn5tG1rc8v7CzSAowppcHly/dCu+KapBnNmAl6+bgP4pfTscNVROzUnA4tnjYNQLWP19OW5+bRuaXcwI0cn9/ZO9+KaoBlaTHk/+cjR3hBFRh/JTbXhpzgRYTXp8eaAKN75aKGsgxCAownaV2rG9pA5xZgNeveF0jM5NlHtJ7UwZnIYXZo+D2aDDf/dW4FcvbsaRmia5l0UK9sKGg3gxsCV+0S9ORX6qTeYVEVE4pKamIjW179O7x+UlYdn1p8Nq0qOywYkmGX/Z5sToToRzYvS6vRVIsBoxtr9yDwTcergGv32lEPXNbsRZDPjjRcNxxbgc/oZPQS6PD3//ZC+WfOU/b+uuGafgtnMHy7wqIlKL7SW16J9sDfn5lZqcGP3cc8+hoKAAFosF48aNw5dfftnl9Rs2bMC4ceNgsVgwYMAALF68OEIrPblzhqYrOgACgAn5yVj9uykYl5eEhhYPFqzYiUuf3YhPd5WfcCAfRZ/CwzX4xfMbgwHQgguGYN45kRnySUTaMKZ/kuwHeKsiE/TWW29h9uzZeO6553DmmWfihRdewEsvvYTdu3ejf//+J1xfVFSEkSNH4sYbb8TNN9+MjRs34tZbb8Ubb7yByy+/vFvPKdcp8krj9vqwLHCoZIPT38nfLzEGF4zMxPRh6RiVkwgbZ8FEhcoGJzbsr8Q7hUfwdVENACDRasTDvzgVF4zMknl1RER+mjtAdeLEiRg7diyef/754G3Dhg3DrFmzsGjRohOuv+eee/Dhhx9iz549wdvmzp2L7777Dps3b+7WczIIaq+q0YklXxXhjW9KUNfmnDGdAAxOj0NeihU5SVZkJpgRbzEizmJEnMUAq0kPvU6AQafz/1cvQK8ToBcEdDYRQEDHXwjFBIG2P+0ixE5ul24TT7jtp9fiJI/R2XO2u60ba2r3jCe5vv1z9+w1uL0iGlo8sDe7Ud/sxnF7C4qqHPixohGHqhzB6/Q6AVeNz8GdPzsF6fGWjhdKRJrR3NyMmTNnAgA+/vhjxMRE5jT63ujJ+7fif4V3uVzYtm0b/vCHP7S7fcaMGdi0aVOH99m8eTNmzJjR7rbzzz8fS5YsgdvthtF44nZ0p9MJp7P1cDe73R6C1WtHaqwZ91wwFHecOxgb9ldiza5ybDpYjXJ7C/Ydb8C+4w1yL5HCTBCAoZnxuHBkJi4fl4PsROX+I0hEoeXz+bBhw4bg51qh+CCoqqoKXq8XGRkZ7W7PyMhAeXl5h/cpLy/v8HqPx4OqqipkZZ2Yul+0aBEefPDB0C1co2JMelwwMhMXjMwEAFTYW7Cr1I4jtU04WtuM4/YWNLR40NDiRkOLBy1uLzw+EV6f2Ppfr6/TvqLO0pKdZkU6uYcots8ctc0utb9duk044baf/qHDazv4ettr2t/W8QN3vJ7erb3dMwidfB64Z0ePpdMJiLcYkRBjRHyMESk2EwpSbchPtWFUTgISraYOn4uISI0UHwRJfvoPvSiKXU5Y7uj6jm6XLFy4EPPnzw/+2W63IzeXA99OJj3ewnIIERGpkuKDoNTUVOj1+hOyPhUVFSdkeySZmZkdXm8wGJCSktLhfcxmM8xmebvUiYiIKHIUv0XeZDJh3LhxWLt2bbvb165di8mTJ3d4n0mTJp1w/Zo1azB+/PgO+4GIiIgo+ig+CAKA+fPn46WXXsLLL7+MPXv24M4770RJSQnmzp0LwF/Kuvbaa4PXz507F8XFxZg/fz727NmDl19+GUuWLMFdd90l10sgIiIihVF8OQwAfvnLX6K6uhp//vOfUVZWhpEjR2L16tXIy8sDAJSVlaGkpCR4fUFBAVavXo0777wTzz77LLKzs/HUU091e0YQERERtWe1KuOcy1BSxZwgOXBOEBERkfpo8tgMIiIiolBiEERERERRiUEQERERdamlpQUXXXQRLrroIrS0tMi9nJBRRWM0ERERycfr9WL16tXBz7WCmSAiIiKKSgyCiIiIKCoxCCIiIqKoxCCIiIiIohKDICIiIopK3B3WCWmQtt1ul3klRERE8nI4HMHP7Xa7oneISe/b3TkQg0FQJxoaGgAAubm5Mq+EiIhIObKzs+VeQrc0NDQgISGhy2t4dlgnfD4fSktLERcXB0EQ5F7OSdntduTm5uLIkSOaPetM669R668P0P5r1PrrA7T/GrX++gDtv0ZRFNHQ0IDs7GzodF13/TAT1AmdToecnBy5l9Fj8fHxmvyhbkvrr1Hrrw/Q/mvU+usDtP8atf76AG2/xpNlgCRsjCYiIqKoxCCIiIiIohKDII0wm824//77YTab5V5K2Gj9NWr99QHaf41af32A9l+j1l8fEB2vsbvYGE1ERERRiZkgIiIiikoMgoiIiCgqMQgiIiKiqMQgiIiIiKISgyCNczqdGD16NARBwI4dO+ReTshccskl6N+/PywWC7KysjB79myUlpbKvayQOXz4MG644QYUFBQgJiYGAwcOxP333w+XyyX30kLmr3/9KyZPngyr1YrExES5lxMSzz33HAoKCmCxWDBu3Dh8+eWXci8pZL744gv8/Oc/R3Z2NgRBwMqVK+VeUkgtWrQIEyZMQFxcHNLT0zFr1izs27dP7mWF1PPPP4/TTjstOCRx0qRJ+Pjjj+VelqwYBGncggULVHPOS0+cc845ePvtt7Fv3z6sWLECBw8exBVXXCH3skJm79698Pl8eOGFF7Br1y488cQTWLx4Me699165lxYyLpcLV155JW655Ra5lxISb731Fn7/+9/jvvvuw/bt2zFlyhTMnDkTJSUlci8tJBwOB0aNGoVnnnlG7qWExYYNGzBv3jxs2bIFa9euhcfjwYwZM9odHKp2OTk5ePjhh1FYWIjCwkKce+65uPTSS7Fr1y65lyYfkTRr9erV4tChQ8Vdu3aJAMTt27fLvaSw+eCDD0RBEESXyyX3UsLmkUceEQsKCuReRsgtXbpUTEhIkHsZfXb66aeLc+fObXfb0KFDxT/84Q8yrSh8AIjvv/++3MsIq4qKChGAuGHDBrmXElZJSUniSy+9JPcyZMNMkEYdP34cN954I1577TVYrVa5lxNWNTU1+Pe//43JkyfDaDTKvZywqa+vR3JystzLoA64XC5s27YNM2bMaHf7jBkzsGnTJplWRX1RX18PAJr9O+f1evHmm2/C4XBg0qRJci9HNgyCNEgURVx33XWYO3cuxo8fL/dywuaee+6BzWZDSkoKSkpK8MEHH8i9pLA5ePAgnn76acydO1fupVAHqqqq4PV6kZGR0e72jIwMlJeXy7Qq6i1RFDF//nycddZZGDlypNzLCanvv/8esbGxMJvNmDt3Lt5//30MHz5c7mXJhkGQijzwwAMQBKHLj8LCQjz99NOw2+1YuHCh3Evuke6+Psndd9+N7du3Y82aNdDr9bj22mshKnwAek9fIwCUlpbiggsuwJVXXonf/va3Mq28e3rz+rREEIR2fxZF8YTbSPluu+027Ny5E2+88YbcSwm5IUOGYMeOHdiyZQtuueUWzJkzB7t375Z7WbLhsRkqUlVVhaqqqi6vyc/Px69+9St89NFH7f7x9Xq90Ov1uOaaa/DKK6+Ee6m90t3XZ7FYTrj96NGjyM3NxaZNmxSd2u3paywtLcU555yDiRMnYtmyZdDplP17S2/+Hy5btgy///3vUVdXF+bVhY/L5YLVasU777yDyy67LHj77373O+zYsQMbNmyQcXWhJwgC3n//fcyaNUvupYTc7bffjpUrV+KLL75AQUGB3MsJu5/97GcYOHAgXnjhBbmXIguD3Aug7ktNTUVqaupJr3vqqafwl7/8Jfjn0tJSnH/++XjrrbcwceLEcC6xT7r7+joixfJOpzOUSwq5nrzGY8eO4ZxzzsG4ceOwdOlSxQdAQN/+H6qZyWTCuHHjsHbt2nZB0Nq1a3HppZfKuDLqLlEUcfvtt+P999/H+vXroyIAAvyvW+n/boYTgyAN6t+/f7s/x8bGAgAGDhyInJwcOZYUUt988w2++eYbnHXWWUhKSsKhQ4fwpz/9CQMHDlR0FqgnSktLMW3aNPTv3x+PPfYYKisrg1/LzMyUcWWhU1JSgpqaGpSUlMDr9QbnWA0aNCj4M6sm8+fPx+zZszF+/HhMmjQJL774IkpKSjTTx9XY2Igff/wx+OeioiLs2LEDycnJJ/ybo0bz5s3D66+/jg8++ABxcXHBXq6EhATExMTIvLrQuPfeezFz5kzk5uaioaEBb775JtavX49PPvlE7qXJR7Z9aRQxRUVFmtoiv3PnTvGcc84Rk5OTRbPZLObn54tz584Vjx49KvfSQmbp0qUigA4/tGLOnDkdvr5169bJvbRee/bZZ8W8vDzRZDKJY8eO1dT26nXr1nX4/2vOnDlyLy0kOvv7tnTpUrmXFjK/+c1vgj+faWlp4vTp08U1a9bIvSxZsSeIiIiIopLymwyIiIiIwoBBEBEREUUlBkFEREQUlRgEERERUVRiEERERERRiUEQERERRSUGQURERBSVGAQRERFRVGIQREQUYocPH4YgCCd8KPGQ2CeffPKEdU6bNk3uZRFFBIMgIo164IEHOnwj7uyDwiM1NRUZGRnIyMgI2SG4N954IwRBQEpKSo8Ovxw0aBAEQcAll1wSvM1mswXXZ7PZQrI+IrVgEEQUBaQ3ua4+KDy2bt2K8vJylJeXIz4+PiSPecMNNwAAampq8MEHH3TrPhs2bMDBgwfb3R/wB1TS+u66666QrI9ILXiKPFEUkE7EJm0444wzMHz4cOzevRtLly7FVVddddL7LF26FIA/IL7ooovCvUQiVWAmiIhIhaRszpo1a3D06NEur21oaMC7774LALj22mthMPD3XyKAQRARdSI/Px+CIGDZsmVwuVx49NFHMWrUKNhsNiQkJODcc8/FJ598ctLH2b59O37zm99g4MCBsFqtiI2NxahRo/B///d/qKqq6vA+Uj+T1KC7YsUKzJgxA+np6dDpdHjggQfaXf/BBx9g+vTpSExMDD7+I488ArfbfcJjAUBtbS2sVisEQcDbb7/d5fr/+Mc/QhAEDBgwAKIonvT19tbKlSsxa9YsZGdnw2QyISkpCWeffTYWL14Mt9t9wvWzZ8+G0WiEz+fDK6+80uVjv/XWW3A4HACA3/zmN2FZP5EqiUSkSffff78IQOztX/O8vDwRgPj000+LEydOFAGIRqNRjI2NDT6uIAjikiVLOn2MP/3pT6IgCMHrrVaraDKZgn/OysoSv/32207XPnXqVHH+/PnB50pKShL1er14//33B6/93//93+DjARATExNFg8EgAhDPPvts8d577w0+Vltz5swRAYjTp0/vdP0ej0fs16+fCED861//2u3vXVFRUXA9RUVFXV7b0NAgXnzxxe1eQ3x8fLvv26RJk8SampoT7nv55ZeLAMRBgwZ1+RyTJ08WAYhnnnlml9e1/b4TRQMGQUQaFaogKCkpSezXr5+4cuVK0eVyiaIoinv37hXPOOMMEYAYGxsr1tXVnXD/J554QgQgxsXFiYsWLRLLyspEUfQHFoWFheK5554rAhBzcnLEhoaGDtcuBVwLFiwQKyoqRFEUxZaWFvHw4cOiKIriG2+8EXyNV199tXj06FFRFEWxublZfPHFF0WLxSImJSV1+Ma+ZcuWYHB18ODBDr8HH374oQhANBgMwfV3R0+CoFmzZgUDmddff1202+3B1/DBBx+IAwYMEAGIs2bNOuG+q1evDj7Phg0bOnz8vXv3Bq/pKmAVRQZBFH0YBBFpVNsgKCMjo8uPO+6444T7S0GQ2WwW9+zZc8LXKyoqRIvFIgIQly9f3u5rlZWVotVqFQVBED/77LMO1+d2u8Vx48aJAMQnnnii07XPnz+/w/v7fD5x8ODBIgDxvPPOE30+3wnXLF26NPg4Hb2xjxkzRgQg/uEPf+jwOaQMzS9+8YsOv96Z7gZB//nPf0QAYmZmZjCA+6kjR46INptNBCBu37693de8Xq+Yk5MjAhDnzJnT4f0XLFgQDCh/Gmz+FIMgijbsCSKKAsePH+/yo76+vtP7XnHFFRg6dOgJt6elpWHSpEkAgJ07d7b72r///W80NTVh/PjxmD59eoePazAY8Otf/xoA8Omnn3Z4jU6nwz333NPh13bs2IEDBw4AAO69994OZx3NmTMH/fv37+SVAXPnzgXg3zn1076bY8eO4eOPPwYA3HzzzZ0+Rl+89NJLAPz9Pf369evwmpycHJxzzjkATvw+6XQ6zJkzBwDw7rvvorGxsd3XvV4vXnvtNQDAL3/5S8TGxoZ0/URqxyCIKAqI/qxvpx/Lli3r9L4TJ07s9GvZ2dkA/PNq2vrqq68AAD/88AMyMzM7/fjzn/8MACguLu7w8QcNGoT09PQOv/btt98CAIxGIyZPntzhNYIgYOrUqZ2u/+qrr0Z8fDyOHz+Ojz76qN3XXn75ZXi9XhQUFOC8887r9DH6Qvo+vfjii11+nz777DMAHX+ffvOb30AQBDgcDrz11lvtvvbxxx+jrKwseB0RtccgiIi6FBcX1+nXpK3WP82ilJaWAgCam5u7zEDZ7XYAQFNTU4eP31kABACVlZUAgJSUFJhMpk6v6yzDAgCxsbG45pprAPgDEYnP58OSJUsAtE5nDjW32x3cHVdfX9/l96mlpQVAx9+nAQMGBHe+vfzyy+2+Jv156NChnQaKRNGMQRARhZzX6wXgLzedLAsliiIOHz7c4ePo9fpOn0MMbFc/WYAiXdeZW265BQCwdu3a4DrWrFmD4uJiGAwGXH/99V3ev7ek7xEAvPnmm936PnWWsZNmBm3atAn79u0DAFRVVeE///lPu68TUXsMgogo5DIzMwEA33//fdieQ8oSVVVVweVydXqdlJXqzKmnnorJkye3y/7861//AgBceumlwdcSahaLBQkJCQD6/n26/PLLkZiYCKB1MvRrr70Gt9sNg8GA2bNn9+nxibSKQRARhdyZZ54JANiyZUun/T59NXbsWAD+stKmTZs6vEYURXzxxRcnfSwpG/Tyyy/j2LFjwf6gm266KUSr7Zj0fXrnnXfg8/l6/TgWiwVXX301AODVV1+F1+sNBkMXX3wxz4Yj6gSDICIKudmzZyMmJgZerxfz5s1rV/r5KZ/Ph7q6uh4/x+jRozFo0CAAwMMPP9xh2Wv58uXdCsKuvPJKpKSkoLS0FFdffTXcbndYG6IlUpC1f/9+PProo11e63A4usx4SSWvsrIyPPTQQ8HsEkthRJ1jEEREIZeZmYmHH34YALBq1Sqcd9552LhxYzAYEkURe/fuxeOPP46RI0cGe1d6QhAEPPjggwD8W8fnzJkTLH21tLRgyZIluPnmm5GUlHTSxzKbzbjuuusAIJg5CldDdFuXXnopLrvsMgDAH/7wB9xyyy3Yv39/8Osulwtff/017rnnHuTl5aGioqLTxxo7dixGjx4NAHjooYcAAFlZWZg5c2b4XgCRyvEUPaIo0J2+lvfeey+kO4juuOMOOJ1OLFy4EOvWrcNZZ50Fk8mEuLg42O32djvKehtsXH311di6dSuefPJJvPbaa1i+fDkSExPR2NgIt9uNc889FxMnTsSiRYtgsVi6fKy5c+fi8ccfhyiKYW2I/qnly5fjhhtuwJtvvonFixdj8eLFsNlsMJlMqK+vb1cmO9n36YYbbsDtt98evM+cOXO6bC4ninbMBBFFgZMNSzx+/HiXpZbeuvvuu7F3717ceeedOO2002CxWFBXV4fY2FhMmDABCxYswKZNm4L9LL3xxBNP4L333sO0adMQFxcHp9OJYcOG4dFHH8Wnn34aPDhUahzuzKBBg4KZlHA2RP+U1WrFG2+8gXXr1mH27NkYMGAAfD4fGhsbkZ6ejnPPPRePPPIIDhw40OV2fwC45ppr2gV7nA1E1DVBPNn+USIiFTvzzDOxadMm/PnPf8Yf//jHTq8rLy9Hbm4uPB4PPv30U8yYMaPXz3n48GEUFBQAAIqKipCfn9/rx4qkBx54AA8++CCmTp2K9evXy70corBjJoiINGvDhg3BnWMXXHBBl9cuXrwYHo8HgwYNCntDNBEpA4MgIlK1efPmYdmyZSgvLw/uEKurq8MLL7yASy+9FABw7rnnYsKECZ0+RmFhIf7xj38AAObPnx/ShuiCggIIggBBEHq1Cy7cnnzyyeD6pEZzomjBxmgiUrWNGzfiueeeA+Df5WW1WlFXVxcMiIYPH45XX321w/vm5+fD6XSivLwcADBmzBj89re/7fOa9Hp9h7N5dDrl/d5ps9lOWGtycrJMqyGKLPYEEZGqffjhh3j//ffxzTff4Pjx46ivr0d8fDxGjBiBX/ziF7jppptgtVo7vK+U8cnMzMQFF1yAhx9+mIMFiaIIgyAiIiKKSsrLzRIRERFFAIMgIiIiikoMgoiIiCgqMQgiIiKiqMQgiIiIiKISgyAiIiKKSgyCiIiIKCoxCCIiIqKo9P9ctim2yr9fWQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "_ = log.get_dos().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aging-extraction",
   "metadata": {},
   "source": [
    "## File I/O\n",
    "Read and write a PDB file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "informed-joyce",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.IO import read_pdb, write_pdb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "uniform-theory",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"scratch/temp.pdb\", \"w\") as ofile:\n",
    "    write_pdb(sys, ofile)\n",
    "with open(\"scratch/temp.pdb\", \"r\") as ifile:\n",
    "    sys = read_pdb(ifile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tracked-dialogue",
   "metadata": {},
   "source": [
    "## Linear Scaling\n",
    "Activate the linear scaling mode of BigDFT and compute some per fragment quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "intended-yugoslavia",
   "metadata": {},
   "outputs": [],
   "source": [
    "from BigDFT.PostProcessing import BigDFTool\n",
    "from BigDFT.IO import read_mol2\n",
    "from io import StringIO\n",
    "from matplotlib import pyplot as plt\n",
    "from BigDFT.Systems import plot_fragment_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "varied-passion",
   "metadata": {},
   "outputs": [],
   "source": [
    "inp = Inputfile()\n",
    "inp.set_hgrid(0.55)\n",
    "inp.set_psp_nlcc()  # Soft pseudopotentials\n",
    "inp.set_xc(\"PBE\")\n",
    "inp[\"import\"] = \"linear\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "expected-panel",
   "metadata": {},
   "outputs": [],
   "source": [
    "istr = \"\"\"@<TRIPOS>MOLECULE\n",
    "N3\n",
    " 48 0 0 0 0\n",
    "SMALL\n",
    "GASTEIGER\n",
    "\n",
    "@<TRIPOS>ATOM\n",
    "      1  N        -15.9520   11.4820   75.1020 N.4     1  ALA         0.3855\n",
    "      2  H        -16.8590   11.1760   75.4220 H       1  ALA        -0.0890\n",
    "      3  H2       -15.9190   12.4910   75.1350 H       1  ALA        -0.0890\n",
    "      4  H3       -15.3270   11.2070   75.8460 H       1  ALA        -0.0890\n",
    "      5  CA       -15.6750   10.9380   73.6710 C.3     1  ALA        -0.0080\n",
    "      6  HA       -15.5000    9.8640   73.7340 H       1  ALA         0.0175\n",
    "      7  CB       -16.9370   11.1850   72.9780 C.3     1  ALA        -0.0702\n",
    "      8  HB1      -16.8760   10.8330   71.9480 H       1  ALA         0.0218\n",
    "      9  HB2      -17.7950   10.7240   73.4660 H       1  ALA         0.0218\n",
    "     10  HB3      -17.0660   12.2620   72.8690 H       1  ALA         0.0218\n",
    "     11  C        -14.3150   11.5090   73.1160 C.2     1  ALA         0.1989\n",
    "     12  O        -13.5240   11.9470   73.9260 O.2     1  ALA        -0.2776\n",
    "     13  N        -14.0770   11.2840   71.7960 N.am    2  VAL        -0.3051\n",
    "     14  H        -14.7740   10.8610   71.1990 H       2  VAL         0.1494\n",
    "     15  CA       -12.9310   11.8890   71.0690 C.3     2  VAL         0.1016\n",
    "     16  HA       -12.1370   12.1190   71.7800 H       2  VAL         0.0597\n",
    "     17  CB       -12.2680   10.8300   70.1710 C.3     2  VAL        -0.0197\n",
    "     18  HB       -13.0570   10.2210   69.7300 H       2  VAL         0.0318\n",
    "     19  CG1      -11.2730   11.3540   69.0870 C.3     2  VAL        -0.0605\n",
    "     20 HG11      -10.7560   10.5590   68.5490 H       2  VAL         0.0233\n",
    "     21 HG12      -11.7530   11.9770   68.3320 H       2  VAL         0.0233\n",
    "     22 HG13      -10.5100   11.8940   69.6470 H       2  VAL         0.0233\n",
    "     23  CG2      -11.4480    9.8090   70.9840 C.3     2  VAL        -0.0605\n",
    "     24 HG21      -11.9590    9.4330   71.8700 H       2  VAL         0.0233\n",
    "     25 HG22      -11.3140    8.9080   70.3860 H       2  VAL         0.0233\n",
    "     26 HG23      -10.4970   10.2560   71.2750 H       2  VAL         0.0233\n",
    "     27  C        -13.2980   13.1490   70.4060 C.2     2  VAL         0.2342\n",
    "     28  O        -14.1940   13.1280   69.5820 O.2     2  VAL        -0.2738\n",
    "     29  N        -12.5920   14.2640   70.6930 N.am    3  LEU        -0.3007\n",
    "     30  H        -11.8060   14.2210   71.3260 H       3  LEU         0.1496\n",
    "     31  CA       -12.8570   15.7250   70.2250 C.3     3  LEU         0.1249\n",
    "     32  HA       -13.7770   15.8070   69.6470 H       3  LEU         0.0616\n",
    "     33  CB       -12.9890   16.4610   71.5050 C.3     3  LEU        -0.0216\n",
    "     34  HB2      -12.3120   16.1130   72.2850 H       3  LEU         0.0291\n",
    "     35  HB3      -12.7420   17.5000   71.2880 H       3  LEU         0.0291\n",
    "     36  CG       -14.4560   16.4840   72.0620 C.3     3  LEU        -0.0445\n",
    "     37  HG       -15.0910   16.5870   71.1830 H       3  LEU         0.0296\n",
    "     38  CD1      -14.9030   15.1850   72.7940 C.3     3  LEU        -0.0626\n",
    "     39 HD11      -15.9610   15.3610   72.9910 H       3  LEU         0.0232\n",
    "     40 HD12      -14.8720   14.2770   72.1920 H       3  LEU         0.0232\n",
    "     41 HD13      -14.2220   15.0500   73.6340 H       3  LEU         0.0232\n",
    "     42  CD2      -14.5970   17.6630   72.9710 C.3     3  LEU        -0.0626\n",
    "     43 HD21      -15.6450   17.5910   73.2620 H       3  LEU         0.0232\n",
    "     44 HD22      -13.9320   17.5610   73.8280 H       3  LEU         0.0232\n",
    "     45 HD23      -14.3580   18.5620   72.4040 H       3  LEU         0.0232\n",
    "     46  C        -11.7070   16.2810   69.3360 C.2     3  LEU         0.3758\n",
    "     47  O        -10.5210   15.8410   69.5920 O.co2   3  LEU        -0.2443\n",
    "     48  OXT      -11.9460   17.0960   68.3940 O.co2   3  LEU        -0.2443\n",
    "\"\"\"\n",
    "Lsys = read_mol2(StringIO(istr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "asian-soundtrack",
   "metadata": {},
   "outputs": [],
   "source": [
    "log = calc.run(sys=Lsys, input=inp, name=\"Lsys\", run_dir=\"scratch\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "native-giving",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHVCAYAAAB8NLYkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABL4UlEQVR4nO3de1xVVf7/8fcB5KAJpwtyS1LGDDWtvHwFbPIyBWKZZfVV06gmM8saQ2s0LAmhJJq+Zo63Ll6mycpmzK5GUqbVCN5GskzJyryCt/Acr6C4f3/44Pw6cj3E4cjm9Xw8zuMxe+219vrsZsz3rHP22hbDMAwBAACg0fPxdgEAAACoHwQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkzBtsJs9e7aioqIUEBCg7t2766uvvqqy77333iuLxVLhc+WVVzr7LFy4sNI+J0+ebIjbAQAAqJEpg93ixYuVnJysJ598Uhs3btR1112nAQMGaOfOnZX2f+mll1RYWOj87Nq1SxdffLH+93//16VfUFCQS7/CwkIFBAQ0xC0BAADUyGIYhuHtIupbTEyMunXrpjlz5jjbOnbsqFtvvVWZmZk1jn/vvfd02223afv27WrTpo2ksyt2ycnJOnz4cJ3rOnPmjPbu3avAwEBZLJY6XwcAADQdhmHoyJEjioiIkI9P9Wtyfg1UU4MpLS3Vhg0b9MQTT7i0JyQkaPXq1bW6xrx583TDDTc4Q125o0ePqk2bNiorK9M111yjjIwMde3atcrrlJSUqKSkxHm8Z88ederUyY27AQAAOGvXrl1q3bp1tX1MF+wOHjyosrIyhYaGurSHhoaqqKioxvGFhYX65JNP9Oabb7q0d+jQQQsXLlSXLl3kcDj00ksv6dprr9U333yj9u3bV3qtzMxMTZkypUL7rl27FBQU5MZd/X9ZWVmaOnWqJk2apIkTJ3p8HAAA8C6Hw6HIyEgFBgbW2Nd0wa7cuV91GoZRq68/Fy5cqAsvvFC33nqrS3tsbKxiY2Odx9dee626deumv//975oxY0al10pJSdH48eOdx+X/xQQFBdU52D377LMKCAhQamqqAgICNHny5BrHZGRkaOrUqUpPT69VfwAAcP6pTY4xXbALDg6Wr69vhdW5/fv3V1jFO5dhGJo/f76SkpLk7+9fbV8fHx/9z//8j7Zt21ZlH6vVKqvVWvvia6k8nKWmprocVyYjI0OpqamEOgAAmgDTBTt/f391795dOTk5Gjx4sLM9JydHt9xyS7VjV61apR9//FEjR46scR7DMJSfn68uXbr87prrojbhjlAHAEDTYrpgJ0njx49XUlKSevToobi4OL3yyivauXOnHnzwQUlnvyLds2ePXn/9dZdx8+bNU0xMjDp37lzhmlOmTFFsbKzat28vh8OhGTNmKD8/X7NmzWqQe6pMdeGOUAcAQNNjymA3dOhQHTp0SOnp6SosLFTnzp21bNky51OuhYWFFfa0s9vtWrJkiV566aVKr3n48GE98MADKioqks1mU9euXfXll1+qZ8+eHr+f6lQW7gh1AAA0Tabcx+585XA4ZLPZZLfb6/zwRFXKw5y/v79KS0sJdQAAmIQ7+YFg14A8Geyksw9rlJaWyt/f32X/PAAA0Hi5kx9M+UqxpigjI8MZ6kpLS5WRkeHtkgAAQAMj2JnAb39TV1JSovT0dKWmphLuAABoYkz58ERTUtmDEu7scwcAAMyDYNeIVff0K+EOAICmh2DXSNVmSxPCHQAATQvBrhFyZ586wh0AAE0Hwa6Rqcvmw4Q7AACaBoJdI1NWVlanzYfL+5eVlXmiLAAAcB5gg+IG5OkNigEAgPmwQTEAAEATRLADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJMg2AEAAJgEwQ4AAMAkCHYAAAAmYdpgN3v2bEVFRSkgIEDdu3fXV199VWXflStXymKxVPhs3brVpd+SJUvUqVMnWa1WderUSUuXLvX0bQAAANSaKYPd4sWLlZycrCeffFIbN27UddddpwEDBmjnzp3VjisoKFBhYaHz0759e+e53NxcDR06VElJSfrmm2+UlJSkIUOGaM2aNZ6+HQAAgFox5btiY2Ji1K1bN82ZM8fZ1rFjR916663KzMys0H/lypXq16+fiouLdeGFF1Z6zaFDh8rhcOiTTz5xtiUmJuqiiy7SW2+9Vau6eFcsAABwV5N+V2xpaak2bNighIQEl/aEhAStXr262rFdu3ZVeHi4rr/+en3xxRcu53Jzcytcs3///tVes6SkRA6Hw+UDAADgKaYLdgcPHlRZWZlCQ0Nd2kNDQ1VUVFTpmPDwcL3yyitasmSJ3n33XUVHR+v666/Xl19+6exTVFTk1jUlKTMzUzabzfmJjIz8HXcGAABQPT9vF+ApFovF5dgwjApt5aKjoxUdHe08jouL065du/TCCy+od+/edbqmJKWkpGj8+PHOY4fDQbgDAAAeY7oVu+DgYPn6+lZYSdu/f3+FFbfqxMbGatu2bc7jsLAwt69ptVoVFBTk8gEAAPAU0wU7f39/de/eXTk5OS7tOTk56tWrV62vs3HjRoWHhzuP4+LiKlxz+fLlbl0TAADAk0z5Vez48eOVlJSkHj16KC4uTq+88op27typBx98UNLZr0j37Nmj119/XZI0ffp0tW3bVldeeaVKS0v1xhtvaMmSJVqyZInzmo8++qh69+6trKws3XLLLXr//ff12Wef6euvv/bKPQIAAJzLlMFu6NChOnTokNLT01VYWKjOnTtr2bJlatOmjSSpsLDQZU+70tJSPf7449qzZ4+aN2+uK6+8Uh9//LFuvPFGZ59evXrp7bff1lNPPaXJkyerXbt2Wrx4sWJiYhr8/gAAACpjyn3szlfsYwcAANzVpPexAwAAaKoIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJgh0AAIBJmDbYzZ49W1FRUQoICFD37t311VdfVdn33XffVXx8vFq1aqWgoCDFxcXp008/demzcOFCWSyWCp+TJ096+lYAAABqxZTBbvHixUpOTtaTTz6pjRs36rrrrtOAAQO0c+fOSvt/+eWXio+P17Jly7Rhwwb169dPN998szZu3OjSLygoSIWFhS6fgICAhrglAACAGlkMwzC8XUR9i4mJUbdu3TRnzhxnW8eOHXXrrbcqMzOzVte48sorNXToUKWmpko6u2KXnJysw4cP17qOkpISlZSUOI8dDociIyNlt9sVFBRU6+sAAICmy+FwyGaz1So/mG7FrrS0VBs2bFBCQoJLe0JCglavXl2ra5w5c0ZHjhzRxRdf7NJ+9OhRtWnTRq1bt9bAgQMrrOidKzMzUzabzfmJjIx072YAAADcYLpgd/DgQZWVlSk0NNSlPTQ0VEVFRbW6xv/93//p2LFjGjJkiLOtQ4cOWrhwoT744AO99dZbCggI0LXXXqtt27ZVeZ2UlBTZ7XbnZ9euXXW7KQAAgFrw83YBnmKxWFyODcOo0FaZt956S2lpaXr//fcVEhLibI+NjVVsbKzz+Nprr1W3bt3097//XTNmzKj0WlarVVartY53AAAA4B6PrtgdOHBAKSkpiouL0xVXXKHNmzdLkl5++eUav8asq+DgYPn6+lZYndu/f3+FVbxzLV68WCNHjtQ777yjG264odq+Pj4++p//+Z9qV+wAAAAakseC3fbt23X11VdrxowZslgs+umnn5wPEmzatKnKVa7fy9/fX927d1dOTo5Le05Ojnr16lXluLfeekv33nuv3nzzTd100001zmMYhvLz8xUeHv67awYAAKgPHvsqdsKECbrwwgu1fv16hYSEyN/f33nuj3/8o55++mlPTa3x48crKSlJPXr0UFxcnF555RXt3LlTDz74oKSzv33bs2ePXn/9dUlnQ93dd9+tl156SbGxsc7VvubNm8tms0mSpkyZotjYWLVv314Oh0MzZsxQfn6+Zs2a5bH7AAAAcIfHgt3nn3+uOXPmKCIiQmVlZS7nwsPDtXfvXk9NraFDh+rQoUNKT09XYWGhOnfurGXLlqlNmzaSpMLCQpc97V5++WWdPn1aDz/8sB5++GFn+z333KOFCxdKkg4fPqwHHnhARUVFstls6tq1q7788kv17NnTY/cBAADgDo/tY9eiRQu9//77io+PV1lZmZo1a6b169erW7du+vjjj3XnnXfK4XB4Yurzljv70AAAAEjnyT520dHR+uyzzyo99+WXX6pz586emhoAAKBJ8thXsaNGjdL48eMVERGhESNGSDq7efC///1vzZ49WzNnzvTU1AAAAE2SR18p9sADD+i1116Tj4+Pzpw5Ix8fHxmGoVGjRmnu3Lmemva8xVexAADAXe7kB4+/KzYvL08ff/yx9u3bp+DgYA0cOLDabUfMjGAHAADc5U5+8PibJ859YwMAAAA8w3TvigUAAGiqPBbsfHx85OvrW+nHz89PwcHBSkxM1BdffOGpEgAAAJoUjwW71NRUtWnTRhdffLHuueceTZgwQUlJSbr44ot12WWX6a677tLu3bsVHx9f4fVfAAAAcJ/HfmN38cUXKywsTN9++60uuOACZ/vRo0cVHx+vSy+9VPn5+YqPj9ezzz6r+Ph4T5UCAADQJHhsxW7GjBl6/PHHXUKdJLVs2VKPP/64Zs+eLT8/Pz344IP673//66kyAAAAmgyPBbvdu3erWbNmlZ7z8/NTUVGRpLPvjT116pSnygAAAGgyPPpKsZdeekmnT592aT99+rReeuklRUdHS5IKCwvVqlUrT5UBAADQZHjsN3bp6em6/fbbdfnll+vWW29VaGio9u3bp/fee0979uzRkiVLJEk5OTmKi4vzVBkAAABNhkffPJGdna3U1FRt2LBBhmHIYrGoR48eSk9PV//+/T017XmLN08AAAB3ef3NE6WlpVq5cqU6deqktWvX6vjx4youLtZFF12kFi1aeGJKAACAJs8jv7Hz8/PTwIEDtW3bNklSixYtdOmllxLqAAAAPMgjwc7Hx0etW7eWw+HwxOUBAABQCY89FTty5EjNmjVLZWVlnpoCAAAAv+Gxp2L9/f1VUFCgjh07atCgQQoPD5fFYnGet1gsGjdunKemBwAAaHI89lSsj0/1i4EWi8Wjq3mzZ8/W3/72NxUWFurKK6/U9OnTdd1111XZf9WqVRo/frw2b96siIgITZgwQQ8++KBLnyVLlmjy5Mn66aef1K5dOz377LMaPHhwrWviqVigaUtLS5Ovr68mT57s9tiMjAyVlZUpLS2t/gsDcF5zJz947KvY7du3V/v5+eefPTW1Fi9erOTkZD355JPauHGjrrvuOg0YMEA7d+6sstYbb7xR1113nTZu3KhJkyZp7Nixzr32JCk3N1dDhw5VUlKSvvnmGyUlJWnIkCFas2aNx+4DgLn4+voqNTVVGRkZbo3LyMhQamqqfH19PVQZANMwTKhnz57Ggw8+6NLWoUMH44knnqi0/4QJE4wOHTq4tI0ePdqIjY11Hg8ZMsRITEx06dO/f39j2LBhta7Lbrcbkgy73V7rMQDMJT093ZBkpKene6Q/APNxJz94bMXOW0pLS7VhwwYlJCS4tCckJGj16tWVjsnNza3Qv3///lq/fr3zPbZV9anqmpJUUlIih8Ph8gHQtE2ePFnp6em1WrkrX6lLT0+v09e3AJoejz08IUlffvmlZsyYoS1btujEiRMu5ywWi3766ad6n/PgwYMqKytTaGioS3toaKiKiooqHVNUVFRp/9OnT+vgwYMKDw+vsk9V15SkzMxMTZkypY53AsCsykNaamqqy/FvEeoA1IXHVuy+/vprXX/99bLb7dqyZYs6dOigSy+9VDt37pSfn5969+7tqaklyeUJXEnOV5q50//cdnevmZKSIrvd7vzs2rWr1vUDMLfqVu4IdQDqymMrdk8//bT+/Oc/a86cOWrWrJmeeeYZdevWTZs2bVJiYqJuu+02j8wbHBwsX1/fCitp+/fvr7DiVi4sLKzS/n5+frrkkkuq7VPVNSXJarXKarXW5TYANAGVrdwR6gD8Hh5bsfvuu+80ePBg54pW+dYmV111lfP/qXqCv7+/unfvrpycHJf2nJwc9erVq9IxcXFxFfovX75cPXr0ULNmzartU9U1AaA2frtyZ7VaCXUAfhePBbvjx4+rZcuW8vHxkdVq1cGDB53nOnTooO+//95TU2v8+PF67bXXNH/+fG3ZskXjxo3Tzp07nfvSpaSk6O6773b2f/DBB7Vjxw6NHz9eW7Zs0fz58zVv3jw9/vjjzj6PPvqoli9frqysLG3dulVZWVn67LPPlJyc7LH7ANA0TJ48Wf7+/iotLZW/vz+hDkCdeSzYXXbZZdq3b58kqVOnTvr444+d51atWuX8itMThg4dqunTpys9PV3XXHONvvzySy1btkxt2rSRJBUWFrrsaRcVFaVly5Zp5cqVuuaaa5SRkaEZM2bo9ttvd/bp1auX3n77bS1YsEBXXXWVFi5cqMWLFysmJsZj9wGgacjIyHCGutLSUrf3uQOAch5788TDDz8si8WimTNnau7cuRozZoz69esnq9Wq5cuX67HHHlNWVpYnpj5v8eYJAOc69zd1/MYOwLncyQ8ee3hiypQp+vXXXyWd/arz+PHjWrRokSwWi5566ik9+eSTnpoaABqFykJcbbZCAYCqeGzFDhWxYgegXE0rc6zcASh3XqzYAQAqV5vQxsodgLrwaLD7+uuv9eabb2rHjh2Vvnni888/9+T0AHDecWcljnAHwF0eC3YLFizQyJEjdfHFF+uKK66osFEv3wADaGrq8vUq4Q6AOzwW7J5//nkNGTJE//jHP3j7AgDo7EbtdfnNXHn/8o3eAaAqHnt4okWLFvrggw90ww03eOLyjRIPTwAAAHe5kx88tkFxx44dnRsUAwAAwPM8FuymTp2q5557Tnv27PHUFAAAAPiNev2N3aBBg1yO7Xa7rrjiCl1zzTUVXiFmsVj0/vvv1+f0AAAATVq9BrtNmzbJYrE4j319fRUSEqK9e/dq7969Ln1/2w8AAAC/X70Gu19++aU+LwcAAAA31Otv7IqLi3X77bfro48+qrLPRx99pNtvv12HDh2qz6kBAACavHoNdq+99pq++eYbJSYmVtknMTFR3377rWbNmlWfUwMAADR59Rrs3n77bY0aNUp+flV/w+vn56dRo0bpgw8+qM+pAQAAmrx6DXY//PCDevToUWO/bt266YcffqjPqQEAAJq8eg12p0+fVrNmzWrs16xZM506dao+pwYAAGjy6jXYhYeH6/vvv6+x3+bNmxUWFlafUwMAADR59Rrs+vTpo9mzZ1e7Gnfq1CnNmTNH/fr1q8+pAQAAmrx6DXbjxo3T1q1bNXjw4AobEkvS3r17deutt6qgoEDjxo2rz6kBAACavHoNdldddZVmzZqlTz/9VFFRUerVq5dGjBihESNGqFevXoqKitLy5cs1a9YsdenSpT6ndiouLlZSUpJsNptsNpuSkpJ0+PDhKvufOnVKEydOVJcuXXTBBRcoIiJCd999d4Vg2rdvX1ksFpfPsGHDPHIPAAAAdWExDMOo74vm5uZq6tSp+uKLL3T8+HFJUosWLXT99dcrJSVFsbGx9T2l04ABA7R792698sorkqQHHnhAbdu21Ycfflhpf7vdrjvuuEOjRo3S1VdfreLiYiUnJ+v06dNav369s1/fvn11xRVXKD093dnWvHlz2Wy2WtfmcDhks9lkt9sVFBRUxzsEAABNiTv5wSPBrtyZM2d08OBBSVJwcLB8fOp1gbCCLVu2qFOnTsrLy1NMTIwkKS8vT3Fxcdq6dauio6NrdZ1169apZ8+e2rFjhy677DJJZ4PdNddco+nTp9e6npKSEpWUlDiPHQ6HIiMjCXYAAKDW3Al2Hk1aPj4+CgkJUUhIiMdDnXR2pdBmszlDnSTFxsbKZrNp9erVtb6O3W6XxWLRhRde6NK+aNEiBQcH68orr9Tjjz+uI0eOVHudzMxM51fCNptNkZGRbt0PAACAO6p+RUQjVFRUpJCQkArtISEhKioqqtU1Tp48qSeeeELDhw93ScUjRoxQVFSUwsLC9N133yklJUXffPONcnJyqrxWSkqKxo8f7zwuX7EDAADwhEYR7NLS0jRlypRq+6xbt06SZLFYKpwzDKPS9nOdOnVKw4YN05kzZzR79myXc6NGjXL+586dO6t9+/bq0aOH/vvf/6pbt26VXs9qtcpqtdY4LwAAQH1oFMHukUceqfEJ1LZt22rTpk3at29fhXMHDhxQaGhoteNPnTqlIUOGaPv27VqxYkWN32F369ZNzZo107Zt26oMdgAAAA2pUQS74OBgBQcH19gvLi5Odrtda9euVc+ePSVJa9askd1uV69evaocVx7qtm3bpi+++EKXXHJJjXNt3rxZp06dUnh4eO1vBAAAwIM8/0RDA+rYsaMSExM1atQo5eXlKS8vT6NGjdLAgQNdnojt0KGDli5dKuns+23vuOMOrV+/XosWLVJZWZmKiopUVFSk0tJSSdJPP/2k9PR0rV+/Xr/88ouWLVum//3f/1XXrl117bXXeuVeAQAAzmWqYCedfXK1S5cuSkhIUEJCgq666ir985//dOlTUFAgu90uSdq9e7c++OAD7d69W9dcc43Cw8Odn/Inaf39/fX555+rf//+io6O1tixY5WQkKDPPvtMvr6+DX6PAAAAlfHoPnZwxQbFAADAXefNPnYAAABoOAQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkyDYAQAAmITpgl1xcbGSkpJks9lks9mUlJSkw4cPVzvm3nvvlcVicfnExsa69CkpKdFf/vIXBQcH64ILLtCgQYO0e/duD94JAACAe0wX7IYPH678/HxlZ2crOztb+fn5SkpKqnFcYmKiCgsLnZ9ly5a5nE9OTtbSpUv19ttv6+uvv9bRo0c1cOBAlZWVeepWAAAA3OLn7QLq05YtW5Sdna28vDzFxMRIkl599VXFxcWpoKBA0dHRVY61Wq0KCwur9Jzdbte8efP0z3/+UzfccIMk6Y033lBkZKQ+++wz9e/fv/5vBgAAwE2mWrHLzc2VzWZzhjpJio2Nlc1m0+rVq6sdu3LlSoWEhOiKK67QqFGjtH//fue5DRs26NSpU0pISHC2RUREqHPnztVet6SkRA6Hw+UDAADgKaYKdkVFRQoJCanQHhISoqKioirHDRgwQIsWLdKKFSv0f//3f1q3bp3+9Kc/qaSkxHldf39/XXTRRS7jQkNDq71uZmam87d+NptNkZGRdbwzAACAmjWKYJeWllbh4YZzP+vXr5ckWSyWCuMNw6i0vdzQoUN10003qXPnzrr55pv1ySef6IcfftDHH39cbV01XTclJUV2u9352bVrVy3vGAAAwH2N4jd2jzzyiIYNG1Ztn7Zt22rTpk3at29fhXMHDhxQaGhorecLDw9XmzZttG3bNklSWFiYSktLVVxc7LJqt3//fvXq1avK61itVlmt1lrPCwAA8Hs0imAXHBys4ODgGvvFxcXJbrdr7dq16tmzpyRpzZo1stvt1Qawcx06dEi7du1SeHi4JKl79+5q1qyZcnJyNGTIEElSYWGhvvvuOz3//PN1uCMAAID61yi+iq2tjh07KjExUaNGjVJeXp7y8vI0atQoDRw40OWJ2A4dOmjp0qWSpKNHj+rxxx9Xbm6ufvnlF61cuVI333yzgoODNXjwYEmSzWbTyJEj9dhjj+nzzz/Xxo0bddddd6lLly7Op2QBAAC8rVGs2Llj0aJFGjt2rPMJ1kGDBmnmzJkufQoKCmS32yVJvr6++vbbb/X666/r8OHDCg8PV79+/bR48WIFBgY6x7z44ovy8/PTkCFDdOLECV1//fVauHChfH19G+7mAAAAqmExDMPwdhFNhcPhkM1mk91uV1BQkLfLAQAAjYA7+cFUX8UCAAA0ZQQ7AAAAkyDYAQAAmATBDgAAoAZpaWnKyMio09iMjAylpaXVb0FVINgBAADUwNfXV6mpqW6Hu4yMDKWmpjbYLhqm2+4EAACgvk2ePFmSlJqa6nJcnfJQl56eXqv+9YFgBwAAUAvuhDtvhDqJYAcAAFBrtQl33gp1EsEOAADALdWFO2+GOolgBwAA4LbKwp23Q53EK8UaFK8UAwDAXMrDnL+/v0pLSz0S6tzJDwS7BkSwAwDAfKxWq0pLS+Xv76+SkpJ6vz7vigUAAGgAGRkZzlBXWlpa502M6wvBDgAAoA5++5u6kpISpaen12kT4/rEwxMAAABuquxBibpsYlzfCHYAAABuqO7pV2+HO4IdAABALdVmSxNvhjuCHQAAQC24s0+dt8IdwQ4AAKAGddl82BvhjmDXgMq3DHQ4HF6uBAAAuOP48eOaNGmSHn30Ubf+Hn/00Ud18uRJHT9+vM5//5ePq83Ww2xQ3IB2796tyMhIb5cBAAAaoV27dql169bV9iHYNaAzZ85o7969CgwMlMVi8XY5AACgETAMQ0eOHFFERIR8fKrfgphgBwAAYBK8eQIAAMAkCHYAAAAmQbADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJMg2AEAAJgEwQ4AAMAkCHYAAAAmQbADAAAwCT9vF9CUnDlzRnv37lVgYKAsFou3ywEAAI2AYRg6cuSIIiIi5ONT/Zocwa4B7d27V5GRkd4uAwAANEK7du1S69atq+1DsGtAgYGBks7+FxMUFOTlagAAQG1lZWVp6tSpmjRpkiZOnOjxcb/lcDgUGRnpzBHVIdg1oPKvX4OCggh2AAA0Is8++6wCAgKUmpqqgIAATZ48ucYxGRkZmjp1qtLT02vVvya1+RkXwQ4AAKAWysNZamqqy3FlMjIylJqaWm+hrrYIdgAAALVUm3DnrVAnEewAAADcUl2482aokwh2AAAAbqss3Hk71EmSxTAMwyszN0EOh0M2m012u52HJwAAMIHyMOfv76/S0lKPhDp38gPBrgER7AAAMB+r1arS0lL5+/urpKSk3q/vTn7glWIAAAB1lJGR4Qx1paWlysjI8Go9BDsAAIA6+O1v6kpKSpSenq7U1FSvhjsengAAAHBTZQ9KuLPPnacQ7AAAANxQ3dOv3g53BDsAAIBaqs2WJt4MdwQ7AACAWnBnnzpvhTuCHQAAQA3qsvmwN8IdwQ4AAKAGZWVlddp8uLx/WVmZJ8qqoFFtd1JcXKykpCTZbDbZbDYlJSXp8OHD1Y4xDENpaWmKiIhQ8+bN1bdvX23evNl5/tdff9Vf/vIXRUdHq0WLFrrssss0duxY2e323z03AAAwh7S0tDqvuE2ePFlpaWn1W1AVGlWwGz58uPLz85Wdna3s7Gzl5+crKSmp2jHPP/+8pk2bppkzZ2rdunUKCwtTfHy8jhw5Iknau3ev9u7dqxdeeEHffvutFi5cqOzsbI0cOfJ3zw0AANCQGs0rxbZs2aJOnTopLy9PMTExkqS8vDzFxcVp69atio6OrjDGMAxFREQoOTlZEydOlCSVlJQoNDRUWVlZGj16dKVz/etf/9Jdd92lY8eOyc/Pr05zl8/121eLOBwORUZG8koxAABQa6Z8pVhubq5sNpszWElSbGysbDabVq9eXemY7du3q6ioSAkJCc42q9WqPn36VDlGkvMfnJ+fX53nlqTMzEznV7c2m02RkZG1vl8AAAB3NZpgV1RUpJCQkArtISEhKioqqnKMJIWGhrq0h4aGVjnm0KFDysjIcFnNq8vckpSSkiK73e787Nq1q8q+AAAAv5fXg11aWposFku1n/Xr10uSLBZLhfGGYVTa/lvnnq9qjMPh0E033aROnTrp6aefrvYatZnbarUqKCjI5QMAAOApXt/u5JFHHtGwYcOq7dO2bVtt2rRJ+/btq3DuwIEDFVbkyoWFhUk6u+IWHh7ubN+/f3+FMUeOHFFiYqJatmyppUuXqlmzZi7XcXduAACAhub1YBccHKzg4OAa+8XFxclut2vt2rXq2bOnJGnNmjWy2+3q1atXpWOioqIUFhamnJwcde3aVZJUWlqqVatWKSsry9nP4XCof//+slqt+uCDDxQQEPC75wYAAGhojeapWEkaMGCA9u7dq5dfflmS9MADD6hNmzb68MMPnX06dOigzMxMDR48WJKUlZWlzMxMLViwQO3bt9fUqVO1cuVKFRQUKDAwUEeOHFF8fLyOHz+upUuX6oILLnBeq1WrVvL19a313DVx56kWAAAAyb384PUVO3csWrRIY8eOdT7lOmjQIM2cOdOlT0FBgcvmwhMmTNCJEyc0ZswYFRcXKyYmRsuXL1dgYKAkacOGDVqzZo0k6fLLL3e51vbt29W2bdtazw0AAOBNjWrFrrFjxQ4AALjLlPvYAQAAoHoEOwAAAJMg2AEAAJgEwQ4AAMAkCHYAAAAmQbADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJMg2AEAAJgEwQ4AAMAkCHYAAAAmQbADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJMg2AEAAJgEwQ4AAMAkCHYAAAAmQbADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJNoVMGuuLhYSUlJstlsstlsSkpK0uHDh6sdYxiG0tLSFBERoebNm6tv377avHmz8/yvv/6qv/zlL4qOjlaLFi102WWXaezYsbLb7S7Xadu2rSwWi8vniSee8MRtAgAA1EmjCnbDhw9Xfn6+srOzlZ2drfz8fCUlJVU75vnnn9e0adM0c+ZMrVu3TmFhYYqPj9eRI0ckSXv37tXevXv1wgsv6Ntvv9XChQuVnZ2tkSNHVrhWenq6CgsLnZ+nnnrKI/cJAABQFxbDMAxvF1EbW7ZsUadOnZSXl6eYmBhJUl5enuLi4rR161ZFR0dXGGMYhiIiIpScnKyJEydKkkpKShQaGqqsrCyNHj260rn+9a9/6a677tKxY8fk5+cn6eyKXXJyspKTk2tdc0lJiUpKSpzHDodDkZGRstvtCgoKqvV1AABA0+VwOGSz2WqVHxrNil1ubq5sNpsz1ElSbGysbDabVq9eXemY7du3q6ioSAkJCc42q9WqPn36VDlGkvMfXHmoK5eVlaVLLrlE11xzjZ599lmVlpZWW3NmZqbza2ObzabIyMja3CoAAECd+NXc5fxQVFSkkJCQCu0hISEqKiqqcowkhYaGurSHhoZqx44dlY45dOiQMjIyKqzmPfroo+rWrZsuuugirV27VikpKdq+fbtee+21KmtOSUnR+PHjncflK3YAAACe4PVgl5aWpilTplTbZ926dZIki8VS4ZxhGJW2/9a556sa43A4dNNNN6lTp056+umnXc6NGzfO+Z+vuuoqXXTRRbrjjjucq3iVsVqtslqt1dYGAABQX7we7B555BENGzas2j5t27bVpk2btG/fvgrnDhw4UGFFrlxYWJiksyt34eHhzvb9+/dXGHPkyBElJiaqZcuWWrp0qZo1a1ZtTbGxsZKkH3/8scpgBwAA0JC8HuyCg4MVHBxcY7+4uDjZ7XatXbtWPXv2lCStWbNGdrtdvXr1qnRMVFSUwsLClJOTo65du0qSSktLtWrVKmVlZTn7ORwO9e/fX1arVR988IECAgJqrGfjxo2S5BIYAQAAvMnrwa62OnbsqMTERI0aNUovv/yyJOmBBx7QwIEDXZ6I7dChgzIzMzV48GBZLBYlJydr6tSpat++vdq3b6+pU6eqRYsWGj58uKSzK3UJCQk6fvy43njjDTkcDjkcDklSq1at5Ovrq9zcXOXl5alfv36y2Wxat26dxo0bp0GDBumyyy5r+H8YAAAAlWg0wU6SFi1apLFjxzqfch00aJBmzpzp0qegoMBlc+EJEyboxIkTGjNmjIqLixUTE6Ply5crMDBQkrRhwwatWbNGknT55Ze7XGv79u1q27atrFarFi9erClTpqikpERt2rTRqFGjNGHCBE/eLgAAgFsazT52ZuDOPjQAAACSSfexAwAAQPUIdgAAACZBsAMAADAJgh0AAIBJEOwAAABMos7BrkuXLnr55Zd1/Pjx+qwHAAAAdVTnYNeqVSs99NBDuvTSSzVu3Dht27atPusCAACAm+oc7FasWKHvvvtOd955p+bNm+d8M8RHH31Un/UBAACgln7Xb+w6deqk2bNna8+ePZo2bZp27NihW265RX/4wx/0wgsvqLi4uL7qBAAAQA3q9c0TRUVFGjFihL744gtJUosWLTR69GhlZGSoRYsW9TVNo8WbJwAAgLsa/M0Tubm5uuuuu9S2bVutXbtWDz30kFauXKnRo0drzpw5uu++++pjGgAAAFTDr64DT548qTfffFOzZs1Sfn6+2rRpo2effVb333+/bDabJKl37966+uqr9fDDD9dbwQAAAKhcnYPdpZdeqsOHD6t3795asmSJbrnlFlkslgr9oqOjdezYsd9VJAAAAGpW52A3ePBgPfroo+rSpUu1/WJiYnTmzJm6TgMAAIBaqvNv7O6++25FRUVVeu7o0aP68ssv61wUAAAA3FfnYNevXz99//33lZ4rKChQv3796lwUAAAA3FfnYFfdLimnTp2Sjw+vofWEtLQ0ZWRk1GlsRkaG0tLS6rcgAABw3nDrN3YOh0OHDx92HhcVFWnnzp0ufU6cOKF//OMfCgsLq5cC4crX11epqamSpMmTJ9d6XEZGhlJTU5Wenu6p0gAAgJe5FexefPFFZzCwWCwaPHhwpf0Mw9CkSZN+f3WooDzMuRPufhvq3AmDAACgcXEr2CUkJKhly5YyDEMTJkzQX/7yF1122WUufaxWq7p06aI+ffrUa6H4/9wJd4Q6AACaDreCXVxcnOLi4iRJx44d06hRoxQREeGRwlC92oQ7Qh0AAE1Lnfexe/rpp+uzDtRBdeGOUAcAQNPjVrB7/fXXddNNN+mSSy7R66+/XmP/u+++u86FoXYqC3eEOgAAmiaLUd2+Jefw8fFRXl6eevbsWeN2JhaLRWVlZb+7QDNxOByy2Wyy2+0KCgqq12uXhzl/f3+VlpYS6gAAMAl38oNbwW7Hjh0KDw+Xv7+/duzYUWP/Nm3a1PbSTYIng5109sGV0tJS+fv7q6SkpN6vDwAAGp47+cGtr2LLg1ppaakKCgrUoUOHCk/FwjsyMjKcoa60tFQZGRms2AEA0MTU6fUQfn5+GjhwoLZt21bf9VSruLhYSUlJstlsstlsSkpKctkwuTKGYSgtLU0RERFq3ry5+vbtq82bN7v0GT16tNq1a6fmzZurVatWuuWWW7R169bfPXdD+e1v6kpKSpSenq7U1NQ6v6ECAAA0TnUKdj4+PmrdurUcDkd911Ot4cOHKz8/X9nZ2crOzlZ+fr6SkpKqHfP8889r2rRpmjlzptatW6ewsDDFx8fryJEjzj7du3fXggULtGXLFn366acyDEMJCQkuvxGsy9wNobIHJSZPnky4AwCgKTLq6JlnnjGuv/564/Tp03W9hFu+//57Q5KRl5fnbMvNzTUkGVu3bq10zJkzZ4ywsDDjueeec7adPHnSsNlsxty5c6uc65tvvjEkGT/++GOd5y6fy263Oz+7du0yJBl2u73W912d9PR0Q5KRnp5ep/MAAOD8Z7fba50f6ryPnb+/vwoKCtSxY0cNGjRI4eHhslgszvMWi0Xjxo37PZnTRW5urmw2m2JiYpxtsbGxstlsWr16taKjoyuM2b59u4qKipSQkOBss1qt6tOnj1avXq3Ro0dXGHPs2DEtWLBAUVFRioyMrPPckpSZmakpU6bU+Z6rU5stTery+jEAANB41TnYTZw40fmfp02bVuF8fQe7oqIihYSEVGgPCQlRUVFRlWMkKTQ01KU9NDS0wlO9s2fP1oQJE3Ts2DF16NBBOTk58vf3r/PckpSSkqLx48c7jx0OhzMs/h7u7FNHuAMAoOmoc7Dbvn17vRSQlpZW46rWunXrJMllRbCcYRiVtv/WuecrGzNixAjFx8ersLBQL7zwgoYMGaL//Oc/CggIqPPcVqtVVqu12trcVZfNhwl3AAA0DXUOdvW1R90jjzyiYcOGVdunbdu22rRpk/bt21fh3IEDByqsyJULCwuTdHbFLTw83Nm+f//+CmPKn3Zt3769YmNjddFFF2np0qW68847FRYW5vbcnlJWVlanzYfL+7NpNAAA5lXnYFdfgoODFRwcXGO/uLg42e12rV27Vj179pQkrVmzRna7Xb169ap0TFRUlMLCwpSTk6OuXbtKOrsH36pVq5SVlVXtfIZhODf5rcvcnpKWllbnsazUAQBgbm69eeK3oqKiqv0a0mKx6KeffqpzYZUZMGCA9u7dq5dfflmS9MADD6hNmzb68MMPnX06dOigzMxMDR48WJKUlZWlzMxMLViwQO3bt9fUqVO1cuVKFRQUKDAwUD///LMWL16shIQEtWrVSnv27FFWVpa++uorbdmyxfnbutrMXRNPv3kCAACYj8fePPFbffr0qRDsDh48qNWrVysoKEh9+vSp66WrtGjRIo0dO9b5lOugQYM0c+ZMlz4FBQWy2+3O4wkTJujEiRMaM2aMiouLFRMTo+XLlyswMFCSFBAQoK+++krTp09XcXGxQkND1bt3b61evdrlgYnazA0AAOBNdV6xq8qhQ4cUHx+vSZMm6Y477qjPSzd6rNgBAAB3uZMf6vTmiepccskl+utf/+qx/dsAAABQuXoPdtLZByJ+/vlnT1waAAAAVaj3YHfq1Cm9+uqrioqKqu9LAwAAoBp1fnjiT3/6U4W2kpIS/fDDD/r111/1j3/843cVBgAAAPfUOdhV9sxFUFCQ7rjjDiUlJTX4/m4AAABNndvB7sSJE3rvvffUv39/tWrVSoMGDVKrVq08URsAAADc4Faw27t3r3r37u3yntjHH39cn3zyiWJjY+u9OAAAANSeWw9PPPXUU9qzZ4+eeuopffTRR3rxxRfl7++vhx56yFP1AQAAoJbcWrHLycnRpEmTnO8cHTBggNq1a6dBgwZp3759Cg0N9UiRAAAAqJlbK3ZFRUXq3bu3S1vfvn1lGIb27dtXr4UBAADAPW4Fu7KyMjVv3tylLSAgQJJ0+vTp+qsKAAAAbnP7qdiCggL5+f3/YWVlZZKkrVu3VujbrVu331EaAAAA3GExKtuQrgo+Pj6yWCwV2g3DcGkvPy4PfTjLnZf4AgAASO7lB7dW7BYsWPC7CgMAAIDnuBXs7rnnHk/VAQAAgN/JrYcnAAAAcP4i2AEAAJgEwQ4AAMAkCHYAAAAmQbADAAAwCYIdAACASRDsAAAATIJgBwAAYBIEOwAAAJMg2AEAAJgEwQ4AAMAkGlWwKy4uVlJSkmw2m2w2m5KSknT48OFqxxiGobS0NEVERKh58+bq27evNm/e7NJn9OjRateunZo3b65WrVrplltu0datW136tG3bVhaLxeXzxBNP1PctAgAA1FmjCnbDhw9Xfn6+srOzlZ2drfz8fCUlJVU75vnnn9e0adM0c+ZMrVu3TmFhYYqPj9eRI0ecfbp3764FCxZoy5Yt+vTTT2UYhhISElRWVuZyrfT0dBUWFjo/Tz31lEfuEwAAoC4shmEY3i6iNrZs2aJOnTopLy9PMTExkqS8vDzFxcVp69atio6OrjDGMAxFREQoOTlZEydOlCSVlJQoNDRUWVlZGj16dKVzbdq0SVdffbV+/PFHtWvXTtLZFbvk5GQlJyfXuuaSkhKVlJQ4jx0OhyIjI2W32xUUFFTr6wAAgKbL4XDIZrPVKj80mhW73Nxc2Ww2Z6iTpNjYWNlsNq1evbrSMdu3b1dRUZESEhKcbVarVX369KlyzLFjx7RgwQJFRUUpMjLS5VxWVpYuueQSXXPNNXr22WdVWlpabc2ZmZnOr41tNluF6wEAANSnRhPsioqKFBISUqE9JCRERUVFVY6RpNDQUJf20NDQCmNmz56tli1bqmXLlsrOzlZOTo78/f2d5x999FG9/fbb+uKLL/TII49o+vTpGjNmTLU1p6SkyG63Oz+7du2q1b0CAADUhdeDXVpaWoWHEs79rF+/XpJksVgqjDcMo9L23zr3fGVjRowYoY0bN2rVqlVq3769hgwZopMnTzrPjxs3Tn369NFVV12l+++/X3PnztW8efN06NChKue1Wq0KCgpy+QAAAHiKn7cLeOSRRzRs2LBq+7Rt21abNm3Svn37Kpw7cOBAhRW5cmFhYZLOrtyFh4c72/fv319hTPnXpe3bt1dsbKwuuugiLV26VHfeeWel146NjZUk/fjjj7rkkkuqrR8AAKAheD3YBQcHKzg4uMZ+cXFxstvtWrt2rXr27ClJWrNmjex2u3r16lXpmKioKIWFhSknJ0ddu3aVJJWWlmrVqlXKysqqdj7DMFwefDjXxo0bJcklMAIAAHiT14NdbXXs2FGJiYkaNWqUXn75ZUnSAw88oIEDB7o8EduhQwdlZmZq8ODBslgsSk5O1tSpU9W+fXu1b99eU6dOVYsWLTR8+HBJ0s8//6zFixcrISFBrVq10p49e5SVlaXmzZvrxhtvlHT2wY28vDz169dPNptN69at07hx4zRo0CBddtllDf8PAwAAoBKNJthJ0qJFizR27FjnU66DBg3SzJkzXfoUFBTIbrc7jydMmKATJ05ozJgxKi4uVkxMjJYvX67AwEBJUkBAgL766itNnz5dxcXFCg0NVe/evbV69WrnwxpWq1WLFy/WlClTVFJSojZt2mjUqFGaMGFCA905AABAzRrNPnZm4M4+NAAAAJJJ97EDAABA9Qh2AAAAJkGwAwAAMAmCHQAAgEkQ7AAAAEyCYAcAAGASBDsAAACTINgBAACYBMEOAADAJAh2AAAAJkGwAwAAMAmCHQAAgEkQ7AAAAEyCYAcAAGASBDsAAACTINgBAACYBMEOAADAJAh2AAAAJkGwAwAAMAmCHQAAgEkQ7AAAAEyCYAcAAGASBDsAAACTINgBAACYBMEOABpIWlqaMjIy6jQ2IyNDaWlp9VsQANMh2AFAA/H19VVqaqrb4S4jI0Opqany9fX1UGUAzKJRBbvi4mIlJSXJZrPJZrMpKSlJhw8frnaMYRhKS0tTRESEmjdvrr59+2rz5s1V9h0wYIAsFovee++93z03APzW5MmTlZ6e7la4Kw916enpmjx5socrBNDYNapgN3z4cOXn5ys7O1vZ2dnKz89XUlJStWOef/55TZs2TTNnztS6desUFham+Ph4HTlypELf6dOny2Kx1NvcAHAud8IdoQ6A24xG4vvvvzckGXl5ec623NxcQ5KxdevWSsecOXPGCAsLM5577jln28mTJw2bzWbMnTvXpW9+fr7RunVro7Cw0JBkLF269HfNXT6X3W53fnbt2mVIMux2u7u3D8Bk0tPTDUlGenp6nc4DaDrsdnut80OjWbHLzc2VzWZTTEyMsy02NlY2m02rV6+udMz27dtVVFSkhIQEZ5vValWfPn1cxhw/flx33nmnZs6cqbCwsHqZW5IyMzOdX93abDZFRka6dc8AzKu6lTtW6gDUlZ+3C6itoqIihYSEVGgPCQlRUVFRlWMkKTQ01KU9NDRUO3bscB6PGzdOvXr10i233FJvc0tSSkqKxo8f7zx2OByEOwBO5aEtNTXVeUyoA/B7eD3YpaWlacqUKdX2WbdunSRV+vs3wzCq/F1cuXPP/3bMBx98oBUrVmjjxo1uXaM2c1utVlmt1mqvC6Bp+224e+aZZ1RaWkqoA1BnXg92jzzyiIYNG1Ztn7Zt22rTpk3at29fhXMHDhyosCJXrvxr1aKiIoWHhzvb9+/f7xyzYsUK/fTTT7rwwgtdxt5+++267rrrtHLlSoWFhbk9NwDU1uTJk52hzt/fn1AHoM68HuyCg4MVHBxcY7+4uDjZ7XatXbtWPXv2lCStWbNGdrtdvXr1qnRMVFSUwsLClJOTo65du0qSSktLtWrVKmVlZUmSnnjiCd1///0u47p06aIXX3xRN998c53nBoDaysjIcIa60tJSZWRkEO4A1I2nn+SoT4mJicZVV11l5ObmGrm5uUaXLl2MgQMHuvSJjo423n33Xefxc889Z9hsNuPdd981vv32W+POO+80wsPDDYfDUeU8Ouep2NrOXRN3nmoB0DSc+/QrT8MCOJc7+cHrK3buWLRokcaOHet8ynXQoEGaOXOmS5+CggLZ7Xbn8YQJE3TixAmNGTNGxcXFiomJ0fLlyxUYGFjvcwOAOyp7UKKyByoAoLYshmEY3i6iqXA4HLLZbLLb7QoKCvJ2OQC8qKanX3k6FkA5d/JDo1qxAwAzqE1oY+UOQF0Q7ACgAbmzEke4A+Augl0DKv/W2+FweLkSAN6QlZWlqVOnatKkSXr00Udr9e+CRx99VCdPnlRqaqpOnjypiRMnNkClAM4n5f+uqM2v5/iNXQPavXs3b54AAAB1smvXLrVu3braPgS7BnTmzBnt3btXgYGBNb4toy7KX1m2a9cuHs4AGhn+/AKNl6f//BqGoSNHjigiIkI+Pj7V9uWr2Abk4+NTY9KuD0FBQfzFADRS/PkFGi9P/vm12Wy16ld97AMAAECjQbADAAAwCYKdiVitVj399NOyWq3eLgWAm/jzCzRe59OfXx6eAAAAMAlW7AAAAEyCYAcAAGASBDsAAACTINgBAACYBMEOAADAJAh2AAAAJkGwM6HTp09r586d3i4DAAA0MIKdCW3evFlRUVHeLgNAFWbPnq0bbrhBQ4YM0YoVK1zOHTx4UH/4wx+8VBmAujIMQ2fOnPF2GQQ7AGhIM2bM0F//+ld16NBBVqtVN954ozIzM53ny8rKtGPHDi9WCKA6p0+f1lNPPaU+ffro6aefliT97W9/U8uWLdW8eXPdc889Ki0t9Vp9fl6bGXXWrVu3as+fOHGigSoB4K6XX35Zr776qoYPHy5JGjNmjG699VadOHFC6enpXq4OQE2mTJmi1157TSNGjNC///1v7d+/Xx9//LFeeeUVnTlzRpMmTdL06dM1YcIEr9THK8UaoYCAAA0bNqzKr1sLCwv16quvqqysrIErA1CTFi1a6Pvvv1fbtm2dbZs3b9b111+vP//5z0pOTlZERAR/foHzVLt27fTSSy9p4MCB+vHHHxUdHa0333xTQ4cOlST961//Unp6ur799luv1MeKXSPUuXNnxcTE6KGHHqr0fH5+vl599dUGrgpAbQQHB2vXrl0uwe7KK6/UihUr9Kc//Ul79uzxXnEAarR3715dffXVkqTLL79c/v7+zmNJ6tGjh1d/TsFv7BqhP/7xjyooKKjyfGBgoHr37t2AFQGorT/+8Y9asmRJhfZOnTrp888/V3Z2theqAlBbNptNhw8fdh5369ZNgYGBzuOSkhJZLBYvVHYWK3aN0PTp06s9365dO33xxRcNUwwAtzzxxBPasGFDpeeuvPJKffHFF/r3v//dwFUBqK1OnTrpv//9r7p06SJJ+s9//uNy/ttvv1X79u29UZokfmMHAABQaz/88IOaNWtW5e/c33zzTfn5+WnIkCENXNlZBDsAAACT4Dd2jVxUVJTi4+Nd2m644QY2OAUaAf78Aqhv/MaukbvnnnvUqlUrl7Zbb71Vhw4d8lJFAGqrsj+/gwcP1sGDB71UEQB39OvXT23atNHChQudbffcc4927dpV4a0yDYWvYgEAAOrgz3/+s8LDwzV16lRn26RJk1RYWKgFCxZ4pSaCnckUFxfrjTfe0Lx585Sfn+/tcgC46aefftKoUaO89v/2ATRufBVrEp999pnmzZun9957T8HBwbrtttu8XRKAOjh69KhWrVrl7TIANFIEu0Zs586dWrBggRYsWKCjR4+quLhY77zzjm6//XZvlwYAgClV9U5nm82m6OhoJSQkyMfHe8+mEuwaoXfeeUevvfaa/vOf/+jGG2/USy+9pAEDBuiCCy5Qx44dvV0eAACmtXTp0krbDx8+rD179ujKK6/Up59+qpCQkAau7CyCXSM0fPhwTZgwQUuWLHF5jQkAAPCsjRs3VnmusLBQw4cP16RJk/Taa681YFX/H8GuEbrvvvs0e/ZsrVq1SklJSRo6dKguuugib5cFoBa6du1a7Xskjx8/3oDVAKhP4eHheuaZZ5SUlOS1Ggh2jdArr7yil156Se+8847mz5+v5ORk9e/fX4Zh6MyZM94uD0A1br31Vm+XAMCDLr30Uu3fv99r87PdiQls27ZN8+fP1+uvv66jR4/qpptu0h133MGTsQAANLD3339fTz75pL777juvzE+wM5EzZ87o448/1rx58/TJJ5+opKTE2yUBcAP7UALnP4fDUWm73W7XunXr9Nhjj+n+++/Xk08+2cCVncVXsSbi4+Ojm2++WYmJidX+uBPA+YV9KIHG48ILL6zyd7IWi0WjR4/WhAkTGriq39TAip35fPPNN+rWrZvKysq8XQqAKrAPJdA4VbWBeFBQkNq3b6+WLVs2cEWuWLEDgAbEPpRA49anTx9vl1At722NDABN0PDhw9WjRw8VFRXpX//6l2655Rb5+/t7uywAtfT888/rxIkTzuMvv/zS5TftR44c0ZgxY7xRmiSCHQA0qPJ9KBMTEzV37lwVFxd7uyQAbkhJSdGRI0ecxwMHDtSePXucx8ePH9fLL7/sjdIk8VVso7Rp06ZqzxcUFDRQJQDcxT6UQON27qMJ59ujCjw80Qj5+PjIYrFU+z8mi8XCwxPAeSg/P1/XXHON83jbtm2aN2+e/vnPf7IPJdAI+Pj4qKioyPku2MDAQH3zzTf6wx/+IEnat2+fIiIivPZ3MMGuEdqxY0eNfYqLi13+8gBwfvDx8VHXrl11//33a8SIEQoKCpLEPpRAY0GwQ4Ox2+1atGiRc3NTVuyA809ubq7mz5+vd955R6dOndJtt92mkSNHql+/fs4++/fvd/6lAeD84uPjo2eeeca5rcnEiRP117/+VcHBwZLOPjyRmppKsEPdrVixQvPnz9e7776rNm3a6Pbbb9ftt9+url27ers0AFU4ceKE3nnnHS1YsEBfffWV2rZtq/vuu0/33HOPWrdu7e3yAFShbdu2VW5QXM5isejnn39uoIrOmZtg1zjt3r1bCxcu1Pz583Xs2DENGTJEc+fO1TfffKNOnTp5uzwAbvjpp5+0YMECvf766yosLFR8fLyWLVvm7bIA1MHOnTuVlpam+fPne2V+gl0jdOONN+rrr7/WwIEDNWLECCUmJsrX11fNmjUj2AGN1NGjR7Vo0SJNmjRJhw8f5qcUQCPl7bc/sd1JI7R8+XKNHTtWDz30kNq3b+/tcgD8DqtWrdL8+fO1ZMkS+fr6asiQIRo5cqS3ywLQSLFBcSP01Vdf6ciRI+rRo4diYmI0c+ZMHThwwNtlAailXbt2KSMjQ+3atVO/fv30008/6e9//7v27t2rV199VbGxsd4uEUAjxVexjdjx48f19ttva/78+Vq7dq3Kyso0bdo03XfffQoMDPR2eQAqER8fry+++EKtWrXS3Xffrfvuu0/R0dHeLgtAPfH2V7EEO5MoKChwbnJ6+PBhxcfH64MPPvB2WQDOMWjQII0cOVIDBw6Ur6+vt8sB4KaaNg8/fPiwVq1aRbBD/SgrK9OHH36o+fPnE+wAAKhnf/7zn2vVb8GCBR6upHIEOwAAAJPg4QkAAACTINgBAACYBMEOAADAJAh2AAAAJkGwA9DkLFy4UBaLpdLP448/7u3yvG7ZsmVKS0vzdhkA6oBXigFoshYsWKAOHTq4tEVERHipmvPHsmXLNGvWLMId0AgR7AA0WZ07d1aPHj1q7Hfq1ClZLBb5+fGvTADnN76KBYDfWLlypSwWi/75z3/qscce06WXXiqr1aoff/xRBw4c0JgxY9SpUye1bNlSISEh+tOf/qSvvvqqwnV2796tO+64Q4GBgbrwwgs1YsQIrVu3ThaLRQsXLnT2u/fee9WyZUtt3bpV/fv31wUXXKDw8HA999xzkqS8vDz98Y9/1AUXXKArrrhC//jHPyrMVVRUpNGjR6t169by9/dXVFSUpkyZotOnTzv7/PLLL7JYLHrhhRc0bdo0RUVFqWXLloqLi1NeXp5LPbNmzZIkl6+of/nll3r6JwzAk/i/nwCarLKyMpfw81spKSmKi4vT3Llz5ePjo5CQEB04cECS9PTTTyssLExHjx7V0qVL1bdvX33++efq27evJOnYsWPq16+ffv31V2VlZenyyy9Xdna2hg4dWulcp06d0m233aYHH3xQf/3rX/Xmm28qJSVFDodDS5Ys0cSJE9W6dWv9/e9/17333qvOnTure/fuks6Gup49e8rHx0epqalq166dcnNz9cwzz+iXX36psPv9rFmz1KFDB02fPl2SNHnyZN14443avn27bDabJk+erGPHjunf//63cnNznePCw8N/zz9qAA3FAIAmZsGCBYakSj85OTmGJKN37941Xuf06dPGqVOnjOuvv94YPHiws33WrFmGJOOTTz5x6T969GhDkrFgwQJn2z333GNIMpYsWeJsO3XqlNGqVStDkvHf//7X2X7o0CHD19fXGD9+vMs1W7ZsaezYscNlrhdeeMGQZGzevNkwDMPYvn27Icno0qWLcfr0aWe/tWvXGpKMt956y9n28MMPG/z1ADROfBULoMl6/fXXtW7dOpdP+e/obr/99krHzJ07V926dVNAQID8/PzUrFkzff7559qyZYuzz6pVqxQYGKjExESXsXfeeWel17RYLLrxxhudx35+frr88ssVHh6url27OtsvvvhihYSEaMeOHc62jz76SP369VNERIROnz7t/AwYMMBZy2/ddNNN8vX1dR5fddVVkuRyTQCNF1/FAmiyOnbsWOHhiZUrV0qq/KvHadOm6bHHHtODDz6ojIwMBQcHy9fXV5MnT3YJdocOHVJoaGiF8ZW1SVKLFi0UEBDg0ubv76+LL764Ql9/f3+dPHnSebxv3z59+OGHatasWaXXPnjwoMvxJZdc4nJstVolSSdOnKh0PIDGhWAHAJWwWCwV2t544w317dtXc+bMcWk/cuSIy/Ell1yitWvXVhhfVFRUv0VKCg4O1lVXXaVnn3220vNs3wI0LQQ7AKgli8XiXOEqt2nTJuXm5ioyMtLZ1qdPH73zzjv65JNPnF+JStLbb79d7zUNHDhQy5YtU7t27XTRRRfVyzV/u4rXvHnzerkmgIbBb+wAoJYGDhyo5cuX6+mnn9aKFSs0Z84c9e/fX1FRUS797rnnHl1++eW66667NGfOHOXk5Gj8+PH69NNPJUk+PvX3r9709HQ1a9ZMvXr10pw5c7RixQotW7ZMs2fP1sCBA7V79263r9mlSxdJUlZWltasWaP169ertLS03moG4Dms2AFALT355JM6fvy45s2bp+eff16dOnXS3LlztXTpUudv8yTpggsu0IoVK5ScnKwJEybIYrEoISFBs2fP1o033qgLL7yw3moKDw/X+vXrlZGRob/97W/avXu3AgMDFRUVpcTExDqt4g0fPlz/+c9/NHv2bKWnp8swDG3fvl1t27att7oBeIbFMAzD20UAQFMwdepUPfXUU9q5c6dat27t7XIAmBArdgDgATNnzpQkdejQQadOndKKFSs0Y8YM3XXXXYQ6AB5DsAMAD2jRooVefPFF/fLLLyopKdFll12miRMn6qmnnvJ2aQBMjK9iAQAATIKnYgEAAEyCYAcAAGASBDsAAACTINgBAACYBMEOAADAJAh2AAAAJkGwAwAAMAmCHQAAgEn8P1KCO+Qj13QYAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "fig, axs = plt.subplots(2, 1)\n",
    "btool = BigDFTool()\n",
    "plot_fragment_information(axs[0], btool.fragment_population(Lsys, log))\n",
    "axs[0].set_ylabel(\"Charge\", fontsize=12)\n",
    "axs[0].set_xticks([])\n",
    "axs[0].set_xlabel(\"\")\n",
    "plot_fragment_information(axs[1], btool.run_compute_purity(Lsys, log))\n",
    "axs[1].set_ylabel(\"Purity\", fontsize=12)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d9191bd",
   "metadata": {},
   "source": [
    "## Remote Calculations\n",
    "Calculations of large systems require the use of powerful supercomputers. One way to run on a supercomputer is to write a python script (or convert a Jupyter notebook), copy the script and related data to the machine, and run the calculation through the job system. As an alternative, we have developed the [remotemanager](https://l_sim.gitlab.io/remotemanager/) framework for automatically launching your jobs. It works by wrapping arbitrary python functions that are then sent and run on the remote machine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "a4a93037",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate(hgrid):\n",
    "    from BigDFT.Calculators import SystemCalculator\n",
    "    from BigDFT.Inputfiles import Inputfile\n",
    "    from BigDFT.Database.Molecules import get_molecule\n",
    "    \n",
    "    sys = get_molecule(\"N2\")\n",
    "    inp = Inputfile()\n",
    "    inp.set_hgrid(hgrid)\n",
    "    calc = SystemCalculator()\n",
    "    log = calc.run(sys=sys, input=inp, name=str(hgrid))\n",
    "    \n",
    "    return log.energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "e8eb72d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "appended run runner-0\n",
      "appended run runner-1\n",
      "appended run runner-2\n",
      "assessing run for runner dataset-54fed8e7-runner-0... checks passed, running\n",
      "assessing run for runner dataset-54fed8e7-runner-1... checks passed, running\n",
      "assessing run for runner dataset-54fed8e7-runner-2... checks passed, running\n"
     ]
    }
   ],
   "source": [
    "from remotemanager import URL  # Replace this with the computer of your choice\n",
    "from remotemanager import Dataset  # Store a set of remote calculations\n",
    "url = URL()\n",
    "ds = Dataset(function=calculate, url=url)\n",
    "for h in [0.3, 0.35, 0.4]:\n",
    "    ds.append_run({\"hgrid\": h})\n",
    "ds.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "fc1583a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "checking remotely for finished runs\n",
      "checking remotely for finished runs\n",
      "[-19.888866344196096, -19.888459196244888, -19.88727486068828]\n"
     ]
    }
   ],
   "source": [
    "from time import sleep\n",
    "while not all(ds.is_finished):\n",
    "    sleep(10)\n",
    "ds.fetch_results()\n",
    "print(ds.results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5a35752",
   "metadata": {},
   "source": [
    "For simpler runs, you can use Jupyter magics for the same effect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "3d587caa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext remotemanager"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "dd339cca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "runner runner-1 already exists\n",
      "assessing run for runner dataset-469c350d-runner-0... skipping already completed run\n",
      "checking remotely for finished runs\n"
     ]
    }
   ],
   "source": [
    "%%sanzu url=url\n",
    "%%sargs hgrid = 0.35\n",
    "from BigDFT.Logfiles import Logfile\n",
    "log = Logfile(\"log-\" + str(hgrid) + \".yaml\")\n",
    "log.log[\"Memory Consumption Report\"][\"Memory occupation\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "fc80d1a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{'Peak Value (MB)': 101.461, 'for the array': 'zt', 'in the routine': 'G_PoissonSolver', 'Memory Peak of process': 'unknown'}]\n"
     ]
    }
   ],
   "source": [
    "print(magic_dataset.results)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}