diff --git a/aiida_base.ipynb b/aiida_base.ipynb new file mode 100644 index 0000000..6f0f4fd --- /dev/null +++ b/aiida_base.ipynb @@ -0,0 +1,409 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# AiiDA \n", + "This implementation uses AiiDA to implement the equation of state workflow with Quantum ESPRESSO.\n", + "It intentionally provides a minimal implementation for the sake of simplicity:\n", + "\n", + "* Quantum ESPRESSO is run locally, but could be run on a remote cluster\n", + "* Not all steps of the procedure are captured by provenance, just the Quantum ESPRESSO calculations\n", + "* The entire procedure is not yet wrapped up in a single workflow\n", + "* The SCF calculations are run serially but could be run in parallel\n", + "\n", + "## Installation and setup\n", + "```\n", + "pip install aiida-core[atomic-tools]\n", + "pip install aiida-shell\n", + "verdi presto\n", + "```\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Profile" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from aiida import load_profile\n", + "\n", + "load_profile()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Implementation of a new simulation code\n", + "\n", + "AiiDA allows developers to implement complex interfaces for any external code, with features such as input validation, error handling etc.\n", + "However, if you just want to run an external code without too much faff, you can use `aiida-shell` to do so." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "from ase.build import bulk\n", + "from aiida.orm import SinglefileData\n", + "\n", + "from functions import write_input\n", + "\n", + "FILEPATH_PSEUDOS = Path.cwd() / 'espresso' / 'pseudo'\n", + "pseudopotentials = {\"Al\": \"Al.pbe-n-kjpaw_psl.1.0.0.UPF\"}\n", + "\n", + "input_string = write_input(\n", + " input_dict={\n", + " \"structure\": bulk('Al'), \n", + " \"pseudopotentials\": pseudopotentials, \n", + " \"kpts\": (3, 3, 3),\n", + " \"calculation\": \"vc-relax\",\n", + " \"smearing\": 0.02,\n", + " },\n", + " return_string=True\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Although the `pw.x` executable can be run without any parsing, let's add a parser that converts the outputs into AiiDA nodes to improve the richness of the provenance:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'xmlschema'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[2], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01maiida\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m orm\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01madis_tools\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mparsers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m parse_pw\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mpw_parser\u001b[39m(dirpath):\n\u001b[1;32m 6\u001b[0m parsed_data \u001b[38;5;241m=\u001b[39m parse_pw(dirpath \u001b[38;5;241m/\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpwscf.xml\u001b[39m\u001b[38;5;124m'\u001b[39m)\n", + "File \u001b[0;32m~/aiida_projects/adis/git-repos/ADIS2023/adis_tools/parsers.py:3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mnumpy\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mxmlschema\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m XMLSchema\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mqe_tools\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m CONSTANTS\n\u001b[1;32m 6\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mase\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Atoms\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'xmlschema'" + ] + } + ], + "source": [ + "from aiida import orm\n", + "from adis_tools.parsers import parse_pw\n", + "\n", + "def pw_parser(dirpath):\n", + " \n", + " parsed_data = parse_pw(dirpath / 'pwscf.xml')\n", + "\n", + " return {\n", + " 'structure': orm.StructureData(ase=parsed_data['ase_structure']),\n", + " 'energy': orm.Float(parsed_data['energy']),\n", + " 'volume': orm.Float(parsed_data['ase_structure'].get_volume()),\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now it's time to run the relaxation:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "from aiida_shell import launch_shell_job\n", + "\n", + "results_relax, node_relax = launch_shell_job(\n", + " 'pw.x',\n", + " arguments='-in {input_file}',\n", + " parser=pw_parser,\n", + " nodes={\n", + " 'input_file': SinglefileData.from_string(input_string, filename='input.pwi'),\n", + " },\n", + " outputs=['pwscf.xml', ],\n", + " metadata={\n", + " 'options': {\n", + " 'prepend_text': f'export ESPRESSO_PSEUDO {FILEPATH_PSEUDOS.as_posix()}',\n", + " 'redirect_stderr': True\n", + " }\n", + " }\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## EOS" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "from functions import generate_structures\n", + "\n", + "scaled_structures = generate_structures(\n", + " results_relax['structure'].get_ase(), [0.9, 0.95, 1.0, 1.05, 1.1])\n", + "\n", + "structures = []\n", + "energies = []\n", + "\n", + "scf_inputs = {\n", + " \"pseudopotentials\": pseudopotentials, \n", + " \"kpts\": (3, 3, 3),\n", + " \"calculation\": \"scf\",\n", + " \"smearing\": 0.02,\n", + "}\n", + "\n", + "for scaled_structure in scaled_structures:\n", + "\n", + " scf_inputs['structure'] = scaled_structure\n", + " input_string = write_input(scf_inputs, return_string=True)\n", + "\n", + " results_scf, node_scf = launch_shell_job(\n", + " 'pw.x',\n", + " arguments='-in {input_file}',\n", + " parser=pw_parser,\n", + " nodes={\n", + " 'input_file': SinglefileData.from_string(input_string, filename='input.pwi'),\n", + " },\n", + " outputs=['pwscf.xml', ],\n", + " metadata={\n", + " 'options': {\n", + " 'prepend_text': f'export ESPRESSO_PSEUDO {FILEPATH_PSEUDOS.as_posix()}',\n", + " 'redirect_stderr': True\n", + " }\n", + " }\n", + " )\n", + "\n", + " structures.append(results_scf['structure'].get_ase())\n", + " energies.append(results_scf['energy'].value)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlwAAAHACAYAAAB3ULYVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABfr0lEQVR4nO3deVxVdeLG8c9lB4XrgoC4AK64L7mhTtriVmplZmVROuXYNFba4mirNuXW5kxO21RoTbtKY5tLaWYCLiXivoIbIi54QVG2e35/oPxCFgWBw70879frviYu59z73DNHePjec79fi2EYBiIiIiJSaVzMDiAiIiLi7FS4RERERCqZCpeIiIhIJVPhEhEREalkKlwiIiIilUyFS0RERKSSqXCJiIiIVDIVLhEREZFKpsIlIiIiUslUuEREREQqmQpXNZWUlMQDDzxAWFgY3t7eNG/enBdeeIHs7OzL7rtjxw6GDx+O1WrF19eXXr16cfDgwYLvp6SkEBkZSVBQELVq1aJr164sXLiw2MfKysqic+fOWCwW4uPjy/QaZs6cSffu3fH19SUgIIBbb72VXbt2lekxREREnIEKVzW1c+dO7HY77777Ltu2beONN97gnXfe4emnny51v3379tG3b1/Cw8P5+eef2bx5M8899xxeXl4F20RGRrJr1y6WLFnCli1bGDFiBHfeeSebNm0q8niTJ08mODi4XK9h9erV/O1vfyMuLo4VK1aQm5vLwIEDOXv2bLkeT0RExGEZ4jDmzJljhIWFlbrNnXfeadx7772lblOrVi3jo48+KnRfvXr1jPfff7/Qfd9//70RHh5ubNu2zQCMTZs2Ffr+tm3bjCFDhhi1atUyAgICjHvvvdc4fvx4ic+bmppqAMbq1atLzSciIuJsNMLlQGw2G/Xq1Svx+3a7ne+++45WrVoxaNAgAgIC6NmzJ19//XWh7fr27csXX3zBqVOnsNvtfP7552RlZdG/f/+CbY4dO8a4ceP4+OOP8fHxKfJcR48epV+/fnTu3JmNGzeydOlSjh07xqhRo0rND5T6GkRERJyRCpeD2LdvH2+++SYPPfRQidukpqZy5swZZs2axeDBg1m+fDm33XYbI0aMYPXq1QXbffHFF+Tm5lK/fn08PT0ZP3480dHRNG/eHADDMBgzZgwPPfQQ3bp1K/a53n77bbp27cqMGTMIDw+nS5cufPjhh6xatYrdu3cX2d4wDB5//HH69u1L+/btr/JoiIiIOBYVrio2bdo0LBZLqbeNGzcW2ic5OZnBgwdzxx138OCDD5b42Ha7HYBbbrmFSZMm0blzZ6ZMmcLQoUN55513CrZ79tlnSUtL48cff2Tjxo08/vjj3HHHHWzZsgWAN998k/T0dKZOnVric/3222+sWrWK2rVrF9zCw8OB/HJ4qQkTJpCQkMBnn3125QdLRETESbiZHaCmmTBhAnfddVep24SGhhb8d3JyMtdddx0RERG89957pe7n7++Pm5sbbdu2LXR/mzZt+PXXX4H8MjRv3jy2bt1Ku3btAOjUqRNr1qzh3//+N++88w4rV64kLi4OT0/PQo/TrVs37rnnHhYsWIDdbmfYsGHMnj27SI6GDRsW+vqRRx5hyZIl/PLLLzRu3LjU1yAiIuKMVLiqmL+/P/7+/le07ZEjR7juuuu45ppriIqKwsWl9AFJDw8PunfvXmTqhd27dxMSEgJAZmYmQJHHcnV1LRgh+9e//sVLL71U8L3k5GQGDRrEF198Qc+ePQHo2rUrixYtIjQ0FDe34k8jwzB45JFHiI6O5ueffyYsLOyKXreIiIjTMfmifSnBkSNHjBYtWhjXX3+9cfjwYePo0aMFtz9q3bq1sXjx4oKvFy9ebLi7uxvvvfeesWfPHuPNN980XF1djTVr1hiGYRjZ2dlGixYtjD/96U/GunXrjL179xqvvvqqYbFYjO+++67YLImJiUU+pXjkyBGjQYMGxsiRI41169YZ+/btM5YtW2aMHTvWyM3NNQzDMP76178aVqvV+Pnnnwvlz8zMrOCjJSIiUr2pcFVTUVFRBlDs7Y8AIyoqqtB9H3zwgdGiRQvDy8vL6NSpk/H1118X+v7u3buNESNGGAEBAYaPj4/RsWPHItNE/FFxhevi49x2221GnTp1DG9vbyM8PNyYOHGiYbfbC7IVd7s0r4iIiLOzGIZhmDCwVmZpaWk8+uijLFmyBIDhw4fz5ptvUqdOnRL3MQyD6dOn895775GWlkbPnj3597//XXDtEsD48eP58ccfSU5Opnbt2vTu3ZvZs2cXXABe3ucWERERuchhPqU4evRo4uPjWbp0KUuXLiU+Pp7IyMhS95kzZw6vv/468+bNY8OGDQQFBTFgwAAyMjIKtrl4fdSOHTtYtmwZhmEwcOBA8vLyruq5RURERC5yiBGuHTt20LZtW+Li4gou2o6LiyMiIoKdO3fSunXrIvsYhkFwcDATJ07k73//O5C/LmBgYCCzZ89m/PjxxT5XQkICnTp1Yu/evTRv3rxczy0iIiLyRw7xKcXY2FisVmtB4QHo1asXVquVmJiYYktPYmIiKSkpDBw4sOA+T09P+vXrR0xMTLGF6+zZs0RFRREWFkaTJk3K/dzFsdvtJCcn4+vri8ViueLXLiIiIuYxDIOMjAyCg4MvO1tAaRyicKWkpBAQEFDk/oCAAFJSUkrcByAwMLDQ/YGBgRw4cKDQfW+99RaTJ0/m7NmzhIeHs2LFCjw8PMr93JA/mpaVlVXw9ZEjR4rMjyUiIiKO4dChQ1c1l6SphWvatGlMnz691G02bNgAUOyokGEYlx0tuvT7xe1zzz33MGDAAI4ePcqrr77KqFGjWLt2LV5eXuV+7pkzZxb72g4dOoSfn1+pmUVERKR6SE9Pp0mTJvj6+l7V45hauK501vWEhASOHTtW5HvHjx8vMoJ1UVBQEJA/QvXHmc9TU1OL7GO1WrFarbRs2ZJevXpRt25doqOjufvuuwkKCirzcwNMnTqVxx9/vODri/+H+fn5qXCJiIg4mKu9HMjUwnWls65HRERgs9lYv349PXr0AGDdunXYbDZ69+5d7D5hYWEEBQWxYsUKunTpAkB2djarV68udjmaPzIMo+DtwPI8N+RfL3bp0jgiIiJSMznEtBBt2rRh8ODBjBs3jri4OOLi4hg3bhxDhw4tdNF6eHg40dHRQH4TnThxIjNmzCA6OpqtW7cyZswYfHx8GD16NAD79+9n5syZ/Pbbbxw8eJDY2FhGjRqFt7c3N910U5meW0RERKQkDnHRPMAnn3zCo48+WvCpw+HDhzNv3rxC2+zatQubzVbw9eTJkzl37hwPP/xwwcSny5cvL3gf1svLizVr1jB37lzS0tIIDAzk2muvJSYmptCF8lfy3CIiIiIlcYh5uJxBeno6VqsVm82ma7hEREQcREX9/naItxRFREREHJkKl4iIiEglU+ESERERqWQqXCIiIiKVTIVLREREpJKpcImIiIhUMhUuERERkUqmwiUiIiJSyVS4HJzdbrBmz3E0f62IiEj1pcLlwAzDYPi/fyXyg/Ws2XPC7DgiIiJSAhUuB2axWOgWUg+AqLWJJqcRERGRkqhwObgxvUOxWGDVruPsP37G7DgiIiJSDBUuBxfqX4vrWwcAsCAmydwwIiIiUiwVLicwtk8YAAt/O0z6+RyT04iIiMilVLicQJ8W9WkVWJuz2Xl8ueGQ2XFERETkEipcTsBisTCmd/4o14LYJPLsmiJCRESkOlHhchK3dWlEHR93Dp06x087jpkdR0RERP5AhctJeHu4clf3pgDM18XzIiIi1YoKlxO5LyIEVxcLMftOsjMl3ew4IiIicoEKlxMJruPN4HZBAMxfm2RuGBERESmgwuVkxvQJBSB60xFOnc02N4yIiIgAKlxOp1tIXdo38iMr185n6w+aHUdERERQ4XI6FouFsRemiPg49gA5eXaTE4mIiIgKlxMa2qkh/rU9SUk/z9KtKWbHERERqfFUuJyQp5sr9/TMnyIiam2iyWlEREREhctJ3dOrKe6uFn4/eJrNh06bHUdERKRGU+FyUgG+XgzrGAxolEtERMRsKlxObGyf/Ivnv9tylNT08yanERERqblUuJxYh8ZWuoXUJSfP4L/rNEWEiIiIWVS4nNzFUa5P1x0gKzfP5DQiIiI1kwqXkxvYLpCGVi9OnMnmm81HzY4jIiJSI6lwOTl3VxciI0KA/IvnDcMwOZGIiEjNo8JVA9zdvSle7i5sS05nQ1Ka2XFERERqHBWuGqBuLQ9u69II0BQRIiIiZlDhqiHGXFhfcdm2FA6nZZqcRkREpGZR4aohWgf50qdFfexG/qLWIiIiUnVUuGqQsRdGuT5bf5DM7FyT04iIiNQcKlw1yPXhAYTU9yH9fC6Lfz9idhwREZEaQ4WrBnFxsXB/RCgA82OSNEWEiIhIFVHhqmHu6NaYWh6u7E09w697T5gdR0REpEZQ4aphfL3cuaNbEwCi1iaZG0ZERKSGUOGqge7vHYrFAit3ppJ44qzZcURERJyeClcNFOZfi+taBwCwICbJ3DAiIiI1gApXDTW2TygAX208RPr5HHPDiIiIODkVrhqqbwt/WgbU5mx2Hl9tPGx2HBEREaemwlVDWSwWxlwY5VoQk0SeXVNEiIiIVBYVrhpsRJfGWL3dOXgqk5U7U82OIyIi4rRUuGowbw9X7upxcYqIRJPTiIiIOC8VrhruvohQXF0sxOw7yc6UdLPjiIiIOCUVrhquUR1vBrULBGC+JkIVERGpFCpcwtg+YQBEbzpC2tlsk9OIiIg4HxUuoVtIXdoF+5GVa+ezDQfNjiMiIuJ0VLgEi8VSMMr1cewBcvLsJicSERFxLipcAsCwTg3xr+3BUdt5lm1LMTuOiIiIU3GYwpWWlkZkZCRWqxWr1UpkZCSnT58udR/DMJg2bRrBwcF4e3vTv39/tm3bVmib8ePH07x5c7y9vWnQoAG33HILO3fuLLRNaGgoFoul0G3KlCkV/RJN5enmyuieIQBE6eJ5ERGRCuUwhWv06NHEx8ezdOlSli5dSnx8PJGRkaXuM2fOHF5//XXmzZvHhg0bCAoKYsCAAWRkZBRsc8011xAVFcWOHTtYtmwZhmEwcOBA8vLyCj3Wiy++yNGjRwtuzz77bKW8TjPd26sp7q4WfjuQRsLh02bHERERcRoWwzCq/ZouO3bsoG3btsTFxdGzZ08A4uLiiIiIYOfOnbRu3brIPoZhEBwczMSJE/n73/8OQFZWFoGBgcyePZvx48cX+1wJCQl06tSJvXv30rx5cyB/hGvixIlMnDix3K8hPT0dq9WKzWbDz8+v3I9T2SZ9EU/0piPc1qURb9zZ2ew4IiIipqqo398OMcIVGxuL1WotKFsAvXr1wmq1EhMTU+w+iYmJpKSkMHDgwIL7PD096devX4n7nD17lqioKMLCwmjSpEmh782ePZv69evTuXNnXn75ZbKznXP6hLEX1lf8NiGZ1PTz5oYRERFxEg5RuFJSUggICChyf0BAACkpxV/gffH+wMDAQvcHBgYW2eett96idu3a1K5dm6VLl7JixQo8PDwKvv/YY4/x+eefs2rVKiZMmMDcuXN5+OGHS82clZVFenp6oZsj6Ni4DteE1CUnz+C/6zRFhIiISEUwtXBNmzatyMXol942btwI5E9dcCnDMIq9/48u/X5x+9xzzz1s2rSJ1atX07JlS0aNGsX58/8/ujNp0iT69etHx44defDBB3nnnXf44IMPOHnyZInPO3PmzIIL/K1Wa5ERs+rs4ijXp+sOkJWbV/rGIiIiclluZj75hAkTuOuuu0rdJjQ0lISEBI4dO1bke8ePHy8ygnVRUFAQkD/S1bBhw4L7U1NTi+xzsRS1bNmSXr16UbduXaKjo7n77ruLfexevXoBsHfvXurXr1/sNlOnTuXxxx8v+Do9Pd1hStegdkE0tHpx1HaebzYfZeQ1jc2OJCIi4tBMLVz+/v74+/tfdruIiAhsNhvr16+nR48eAKxbtw6bzUbv3r2L3ScsLIygoCBWrFhBly5dAMjOzmb16tXMnj271OczDIOsrKwSv79p0yaAQkXuUp6ennh6epb6PNWVu6sL9/YK4ZVlu4ham8jtXRtddiRRRERESuYQ13C1adOGwYMHM27cOOLi4oiLi2PcuHEMHTq00CcUw8PDiY6OBvLfSpw4cSIzZswgOjqarVu3MmbMGHx8fBg9ejQA+/fvZ+bMmfz2228cPHiQ2NhYRo0ahbe3NzfddBOQf8H+G2+8QXx8PImJiXz55ZeMHz+e4cOH07Rp06o/GFVkdI+meLq5sC05nY0H0syOIyIi4tBMHeEqi08++YRHH3204FOHw4cPZ968eYW22bVrFzabreDryZMnc+7cOR5++GHS0tLo2bMny5cvx9fXFwAvLy/WrFnD3LlzSUtLIzAwkGuvvZaYmJiCi/Q9PT354osvmD59OllZWYSEhDBu3DgmT55cRa/cHHVreXBbl0Z8vuEQUWsT6R5az+xIIiIiDssh5uFyBo4yD9cf7UxJZ/DcNbi6WPhl8nU0quNtdiQREZEqVaPm4RJzhAf50bt5ffLsBh/FJpkdR0RExGGpcEmpxvYJA+Dz9YfIzM41OY2IiIhjUuGSUl0fHkDTej7YzuUQvemI2XFEREQckgqXlMrVxcL9vUMBmL82CV3yJyIiUnYqXHJZd3RrTC0PV/aknuHXvSfMjiMiIuJwVLjksvy83LmjW/4s+VFrk8wNIyIi4oBUuOSK3BcRAsDKnakknjhrchoRERHHosIlV6RZg9pc17oBAAtikswNIyIi4mBUuOSKXZwiYuFvh8k4n2NyGhEREcehwiVX7E8t/WkRUJszWbl8tfGw2XFEREQchgqXXDGLxcKYC1NELIhNIs+uKSJERESuhAqXlMmIro3w83LjwMlMVu1MNTuOiIiIQ1DhkjLx8XDj7h5NAYiKSTQ5jYiIiGNQ4ZIyi4wIwcUCa/eeZFdKhtlxREREqj0VLimzxnV9GNQuCID5GuUSERG5LBUuKZeLU0Qs/v0IaWezTU4jIiJSvalwSbl0D61Lu2A/snLtfLbhoNlxREREqjUVLimXP04R8XHsAXLy7OYGEhERqcZUuKTchnUKpn4tD47azrN82zGz44iIiFRbKlxSbl7urtzT88IUEWt18byIiEhJVLjkqtzbKwR3VwsbD6Sx5bDN7DgiIiLVkgqXXJUAPy9u7tAQ0CiXiIhISVS45KpdnCLim4RkUjPOm5xGRESk+lHhkqvWqUkdujatQ06ewSdxmiJCRETkUipcUiEujnJ9su4AWbl5JqcRERGpXlS4pEIMbh9EkJ8XJ85k8+3mo2bHERERqVZUuKRCuLu6EBkRAkBUTCKGYZicSEREpPpQ4ZIKc3ePpni6ubD1SDobD6SZHUdERKTaUOGSClOvlge3dm4EaIoIERGRP1Lhkgo1tm8oAMu2HePI6XPmhhEREakmVLikQoUH+RHRrD55doOPYw+YHUdERKRaUOGSCje2TygAn60/yLlsTREhIiKiwiUV7oY2gTSp543tXA7Rm46YHUdERMR0KlxS4VxdLNwfEQrAfE0RISIiosIllWNU9ybU8nBl97EzrN170uw4IiIiplLhkkrh5+XOyGsaA5oiQkRERIVLKs39vUMBWLkrlaQTZ80NIyIiYiIVLqk0zRrU5rrWDTAMmB+TZHYcERER06hwSaUa0ycMgIW/HSbjfI7JaURERMyhwiWV6tqW/jRvUIszWbks/O2w2XFERERMocIllcpisRSMci2IScJu1xQRIiJS86hwSaW7vWsj/LzcSDqZyapdqWbHERERJ7Nmz3E+/DWRvGr8R70Kl1Q6Hw837urRFICotUnmhhEREadyJiuXKYu28OK326v1NEQqXFIl7osIwcUCv+49we5jGWbHERERJ/HK0p0cOX2OxnW9ufvCH/fVkQqXVInGdX0Y2DYI0CiXiIhUjA1Jp1gQewCAWSM6UsvTzeREJVPhkioztk8oANGbDnM6M9vcMCIi4tDO5+Tx94UJANzZrQl9W/qbnKh0KlxSZXqE1aNtQz/O59j5bP0hs+OIiIgDm/vjHvafOEuArydP39zG7DiXpcIlVcZisRSMcn0cm0Runt3cQCIi4pC2HLbxnzX7AXjp1vZYvd1NTnR5KlxSpYZ1CqZ+LQ+SbedZtu2Y2XFERMTB5OTZmbwogTy7wdCODRnYLsjsSFdEhUuqlJe7K6N7Xpwiovp+fFdERKqnd37ex46j6dT1cWf68HZmx7liKlxS5e7tFYKbi4WNB9LYcthmdhwREXEQe45l8ObKvQBMG96O+rU9TU505VS4pMoF+nlxc8eGAETFaJRLREQuL89u8NTCBLLz7NwQHsDwTsFmRyoTFS4xxdgL6yt+u/koxzOyTE4jIiLV3fyYJOIPncbX042XbmuPxWIxO1KZqHCJKTo3qUOXpnXIzrPzyboDZscREZFq7ODJTF5dtguAqTe1oaHV2+REZecwhSstLY3IyEisVitWq5XIyEhOnz5d6j6GYTBt2jSCg4Px9vamf//+bNu2rcRthwwZgsVi4euvv77q55bLuzjK9d+4g2Tl5pmcRkREqiPDMJiyOIFzOXlENKvP3T2amB2pXBymcI0ePZr4+HiWLl3K0qVLiY+PJzIystR95syZw+uvv868efPYsGEDQUFBDBgwgIyMomv5zZ07t8ThyfI8t1zekPZBBPl5ceJMFt8lHDU7joiIVEOfbzhEzL6TeLm7MOv2Dg73VmIBwwFs377dAIy4uLiC+2JjYw3A2LlzZ7H72O12IygoyJg1a1bBfefPnzesVqvxzjvvFNo2Pj7eaNy4sXH06FEDMKKjo6/quYtjs9kMwLDZbFe8T00wb+UeI+Tv3xpD/7XGsNvtZscREZFqJPl0ptH++aVGyN+/Nf7zyz5TMlTU72+HGOGKjY3FarXSs2fPgvt69eqF1WolJiam2H0SExNJSUlh4MCBBfd5enrSr1+/QvtkZmZy9913M2/ePIKCik6eVp7nlit3d4+meLi5sOWIjd8OpJkdR0REqgnDMHg2eisZWbl0blKn4DIUR+UQhSslJYWAgIAi9wcEBJCSklLiPgCBgYGF7g8MDCy0z6RJk+jduze33HJLhT03QFZWFunp6YVuUlS9Wh7c2jn/o71Ra5PMDSMiItXGks3J/LQzFXdXC3NGdsTVxUHfSrzA1MI1bdo0LBZLqbeNGzcCFPuerWEYl30v99Lv/3GfJUuWsHLlSubOnVumx7iS5545c2bBRfZWq5UmTRzzIr+qcPGvlqXbUkg+fc7kNCIiYraTZ7KY/s12AB65viWtAn1NTnT1TC1cEyZMYMeOHaXe2rdvT1BQEMeOFV137/jx40VGsC66+PbgpaNQqampBfusXLmSffv2UadOHdzc3HBzcwPg9ttvp3///gWPU9bnBpg6dSo2m63gdujQocsfkBqqTUM/ejWrR57d4KNYTREhIlLTTftmO6fOZhMe5Mtf+zc3O06FcDPzyf39/fH397/sdhEREdhsNtavX0+PHj0AWLduHTabjd69exe7T1hYGEFBQaxYsYIuXboAkJ2dzerVq5k9ezYAU6ZM4cEHHyy0X4cOHXjjjTcYNmxYuZ8b8q8X8/R0nCUHzDa2Txhx+0/x+YaDPHZDS7w9XM2OJCIiJlix/RjfbE7GxQKvjOyEu6tDXP10WaYWrivVpk0bBg8ezLhx43j33XcB+Mtf/sLQoUNp3bp1wXbh4eHMnDmT2267DYvFwsSJE5kxYwYtW7akZcuWzJgxAx8fH0aPHg3kj14Vd6F806ZNCQsLK9Nzy9W5sU0gTep5c+jUOb6OP8LdPZqaHUlERKqY7VwOz369BYBx1zajQ2OryYkqjsPUxk8++YQOHTowcOBABg4cSMeOHfn4448LbbNr1y5stv9fDHny5MlMnDiRhx9+mG7dunHkyBGWL1+Or2/Z3gu+kueWq+PqYuH+iFAAotYmYhiGuYFERKTKzfx+B8fSswjzr8WkG1uZHadCWQz9ZqsS6enpWK1WbDYbfn5+ZseplmzncoiY+ROZ2Xl88mBP+rS4/NvNIiLiHNbuPcE9768D4MvxEfQIq2dyonwV9fvbYUa4xPlZvd0ZeU1jIH+US0REaobM7FymLE4A4L6IkGpTtiqSCpdUK/f3DgXgp52pHDh51twwIiJSJV5ZtotDp87RqI43kweHmx2nUqhwSbXSvEFt+rdugGHA/Jgks+OIiEgl++1AWsHP+xkjOlDb0yE+z1dmKlxS7VycCPWrjYfJOJ9jchoREaksWbl5/H1RAoYBt3dtTL9WDcyOVGlUuKTa+VMLf5o1qMWZrFwW/nbY7DgiIlJJ3vxpL3tTz+Bf25PnhrYxO06lUuGSasfFxcLYC9dyLYhJwm7XB2lFRJzNtmQbb6/eB8BLt7ajjo+HyYkqlwqXVEsjujbG18uNpJOZrNqVanYcERGpQLl5diYvTCDPbjCkfRCD2zc0O1KlU+GSaqmWpxt3dc9f8FsXz4uIOJf31uxnW3I6Vm93pt/Szuw4VUKFS6qt+yJCcbHAmj0n2HMsw+w4IiJSAfYdP8PcH/cA8PzQtgT4epmcqGqocEm11aSeDwPaBgIQpVEuERGHZ7cb/H1hAtm5dvq1asCIro3MjlRlVLikWrs4RcTi3w9zOjPb5DQiInI1PopNYuOBNGp5uDJjRAcsFovZkaqMCpdUaz3D6tGmoR/nc+x8vuGQ2XFERKScDp3KZM6yXQBMGRJOozreJieqWipcUq1ZLBbG9gkF4KOYJHLz7OYGEhGRMjMMg6ejt5CZnUePsHrc0zPE7EhVToVLqr3hnYKpV8uDZNt5lm8/ZnYcEREpo69+O8yaPSfwdHNh9u0dcXGpOW8lXqTCJdWel7sro3s0BSBqbaLJaUREpCxS08/z0rfbAXh8QCvC/GuZnMgcKlziECIjQnBzsbAhKY2tR2xmxxERkStgGAbPfr2V9PO5dGxs5YG+YWZHMo0KlziEQD8vbuqQPxPxhxrlEhFxCN9vSWH59mO4uViYfXtH3Fxrbu2oua9cHM7Fi+e/3XyU4xlZ5oYREZFSpZ3N5oUlWwF4+LoWtGnoZ3Iic6lwicPo0rQunZvUITvPzqfrDpodR0RESvHit9s5cSabVoG1mXBdC7PjmE6FSxzKxVGu/647QHaupogQEamOVu1MJXrTEVwsMGdkJzzcVDd0BMSh3NShIYF+nhzPyOK7LclmxxERkUtknM/h6egtADzQN4zOTeqYG6iaUOESh+Lu6kJkr/wJ86LWJmEYhsmJRETkj2b9sJOjtvOE1Pfh8QGtzY5TbahwicO5u0dTPNxcSDhs4/eDaWbHERGRC2L3neSTC9fYzhrREW8PV5MTVR8qXOJw6tf25NbOwQB8uDbJ3DAiIgLAuew8pixOAGB0z6ZENK9vcqLqRYVLHNKY3vmT5y3dmkLy6XMmpxERkddX7OLAyUwaWr2YOiTc7DjVjgqXOKS2wX70DKtHnt3g47gDZscREanR4g+d5oNf8yelfvm29vh6uZucqPopV+E6e/ZsRecQKbOxffJHuT5bf5Bz2XkmpxERqZmyc+1MXrgZuwG3dg7m+vBAsyNVS+UqXIGBgfz5z3/m119/reg8IldsQNtAGtf15nRmDl/HHzE7johIjfTvVXvZfewM9Wt58PywdmbHqbbKVbg+++wzbDYbN9xwA61atWLWrFkkJ2tOJKlari4W7o8IBSBqbaKmiBARqWI7U9L596q9AEy/pR31anmYnKj6KlfhGjZsGIsWLSI5OZm//vWvfPbZZ4SEhDB06FAWL15Mbm5uRecUKdao7k3w8XBl97EzxO47aXYcEZEaIzfPzuSFCeTaDQa2DeTmDg3NjlStXdVF8/Xr12fSpEls3ryZ119/nR9//JGRI0cSHBzM888/T2ZmZkXlFCmW1dud27s2BjRFhIhIVfpwbSIJh234ernxj1vbY7FYzI5UrV1V4UpJSWHOnDm0adOGKVOmMHLkSH766SfeeOMNoqOjufXWWysopkjJxlxYX/Gnncc4cFIf6BARqWyJJ87y2vLdADx3c1sC/bxMTlT9uZVnp8WLFxMVFcWyZcto27Ytf/vb37j33nupU6dOwTadO3emS5cuFZVTpETNG9SmX6sGrN59nAUxB3h+WFuzI4mIOC273eDvixLIyrXzp5b+3NGtsdmRHEK5RrjGjh1LcHAwa9euJT4+ngkTJhQqWwDNmjXjmWeeqYiMIpc19sIo11cbD3EmS9cQiohUlk/WH2R94il8PFyZcVsHvZV4hco1wnX06FF8fHxK3cbb25sXXnihXKFEyuralg1o1qAW+4+fZeHGQ4y5MEeXiIhUnCOnzzHr+x0ATB7Umib1Su8C8v/KNcKVm5tLenp6kVtGRgbZ2dkVnVHkslxcLIzpHQrAgtgD2O2aIkJEpCIZhsEz0Vs4m53HNSF1ue/CtDxyZcpVuOrUqUPdunWL3OrUqYO3tzchISG88MIL2O32is4rUqLbuzbG18uNxBNn+Xl3qtlxREScSvSmI/y86zgebi7Mvr0jLi56K7EsylW45s+fT3BwME8//TRff/010dHRPP300zRq1Ii3336bv/zlL/zrX/9i1qxZFZ1XpES1PN24s1sTAKI0RYSISIU5npHFi99uB+CxG1rSIqC2yYkcT7mu4VqwYAGvvfYao0aNKrhv+PDhdOjQgXfffZeffvqJpk2b8vLLL/P0009XWFiRy7m/dygfrk1kzZ4T7DmWQctAX7MjiYg4vBeWbOV0Zg7tgv34y7XNzI7jkMo1whUbG1vslA9dunQhNjYWgL59+3Lw4MGrSydSRk3q+XBjm/yFU6NikswNIyLiBJZuPcr3W1JwdbEwZ2RH3F2vagrPGqtcR61x48Z88MEHRe7/4IMPaNIk/y2dkydPUrdu3atLJ1IOYy98QnHx74exZeaYnEZExHHZMnN49uttADzUrxntgq0mJ3Jc5XpL8dVXX+WOO+7ghx9+oHv37lgsFjZs2MDOnTtZuHAhABs2bODOO++s0LAiV6JXs3qEB/myMyWDzzccZHy/5mZHEhFxSP/4bjsnzmTRvEEtHrm+pdlxHJrFMIxyfX7+wIEDvPPOO+zatQvDMAgPD2f8+PGEhoZWcETnkJ6ejtVqxWaz4efnZ3Ycp/flhkNMXpRAozrerH6qP24aAhcRKZPVu49z/4frsVhg4UO9uSakZr5rVVG/v8s8wpWTk8PAgQN59913mTlzZrmfWKQyDe8czKylOzly+hwrth9jiFaxFxG5Ymeycnl68RYAxvQOrbFlqyKV+c9+d3d3tm7dqqn8pVrzcndldI+mgKaIEBEpqzkX/mBtXNebpwa1NjuOUyjX+yz33XdfsRfNi1Qn9/YKwc3FwvqkU2w9YjM7joiIQ1ifeIqPYg8AMGtER3w8ynW5t1yiXEcxOzub999/nxUrVtCtWzdq1apV6Puvv/56hYQTuRpBVi+GdGjIN5uTiVqbxGujOpkdSUSkWjufk8eURQkA3NmtCX1b+pucyHmUq3Bt3bqVrl27ArB79+5C39NbjVKdjO0Tyjebk/lmczJThoTTwNfT7EgiItXW3B/3sP/EWQL9PHn65jZmx3Eq5Spcq1atqugcIpWia9O6dGpSh82HTvPpuoM8dqM+1iwiUpwth238Z81+AF66tQNWb3eTEzmXq/qs/N69e1m2bBnnzp0D8lcSF6lu/twnFID/rjtAdq4WVBcRuVR2rp2nFm4mz24wrFMwA9oGmh3J6ZSrcJ08eZIbbriBVq1acdNNN3H06FEAHnzwQZ544okKDShytYa0b0iAryfHM7L4bkuy2XFERKqdd1fvY2dKBnV93Jk2rK3ZcZxSuQrXpEmTcHd35+DBg/j4+BTcf+edd7J06dIKCydSETzcXIjsFQLkTxGhkVgRkf+351gGb67cC8C04e2oX1vXulaGchWu5cuXM3v2bBo3blzo/pYtW3LgwIEKCSZSkUb3bIqHmwsJh238fvC02XFERKqFPLvBUwsTyM6zc0N4AMM7BZsdyWmVq3CdPXu20MjWRSdOnMDTU81Yqp/6tT255cIPkqi1iSanERGpHqLWJhJ/6DS+nm68dFt7zTRQicpVuK699lo++uijgq8tFgt2u51XXnmF6667rsLC/VFaWhqRkZFYrVasViuRkZGcPn261H0Mw2DatGkEBwfj7e1N//792bZtW4nbDhkyBIvFwtdff13oe6GhoVgslkK3KVOmVNArk6oy5sLF8z9sTeGo7Zy5YURETHbwZCavLt8FwNSb2tDQ6m1yIudWrsL1yiuv8O677zJkyBCys7OZPHky7du355dffmH27NkVnRGA0aNHEx8fz9KlS1m6dCnx8fFERkaWus+cOXN4/fXXmTdvHhs2bCAoKIgBAwaQkZFRZNu5c+eW2uxffPFFjh49WnB79tlnr/o1SdVqF2ylR1g98uwGH8fqrW8RqbkMw2DK4gTO59iJaFafu3s0MTuS0ytX4Wrbti0JCQn06NGDAQMGcPbsWUaMGMGmTZto3rx5RWdkx44dLF26lPfff5+IiAgiIiL4z3/+w7fffsuuXbuK3ccwDObOncszzzzDiBEjaN++PQsWLCAzM5NPP/200LabN2/m9ddf58MPPywxg6+vL0FBQQW32rVrV+hrlKpxcYqIz9Yf5HxOnrlhRERM8vmGQ8TsO4mXuwuzbu+gtxKrQLnn4QoKCmL69Ol8++23fP/997z00ks0bNiwIrMViI2NxWq10rNnz4L7evXqhdVqJSYmpth9EhMTSUlJYeDAgQX3eXp60q9fv0L7ZGZmcvfddzNv3jyCgoJKzDB79mzq169P586defnll8nOzq6AVyZVbUDbIBrV8SYtM4evNx0xO46ISJU7ajvHjO92APDkwNaE1K91mT2kIpR7RcrTp0+zfv16UlNTsdsLTyZ53333XXWwP0pJSSEgIKDI/QEBAaSkpJS4D0BgYOHJ2wIDAwt9knLSpEn07t2bW265pcTnf+yxx+jatSt169Zl/fr1TJ06lcTERN5///0S98nKyiIrK6vg6/T09BK3larj6mLh/t4hzPh+J1Frk7izexP9ZSciNYZhGDwbvZWMrFw6N6nD2D5hZkeqMcpVuL755hvuuecezp49i6+vb6FfWBaL5YoL17Rp05g+fXqp22zYsKHgcS9lGMZlf1le+v0/7rNkyRJWrlzJpk2bSn2MSZMmFfx3x44dqVu3LiNHjiwY9SrOzJkzL/vaxBx3dmvKGyv2sOtYBrH7TtK7hRZnFZGaYcnmZH7amYq7q4U5Izvi6qI/OKtKud5SfOKJJ/jzn/9MRkYGp0+fJi0treB26tSpK36cCRMmsGPHjlJv7du3JygoiGPHjhXZ//jx40VGsC66+PbgpSNgqampBfusXLmSffv2UadOHdzc3HBzy++ft99+O/379y8xd69evYD8pY1KMnXqVGw2W8Ht0KFDJR8IqVJWH3duv6YRAB+uTTI3jIhIFTl5Jovp32wH4JHrW9Iq0NfkRDVLuUa4jhw5wqOPPlrsXFxl4e/vj7//5UcXIiIisNlsrF+/nh49egCwbt06bDYbvXv3LnafsLAwgoKCWLFiBV26dAEgOzub1atXF3yScsqUKTz44IOF9uvQoQNvvPEGw4YNKzHPxRGx0q5Z8/T01Jxk1diY3mH8N+4gP+08xsGTmTStf3XnsohIdTftm+2cOptNeJAvf+1f8R9wk9KVa4Rr0KBBbNy4saKzlKhNmzYMHjyYcePGERcXR1xcHOPGjWPo0KG0bt26YLvw8HCio6OB/LcSJ06cyIwZM4iOjmbr1q2MGTMGHx8fRo8eDeSPgrVv377QDaBp06aEheW/rx0bG8sbb7xBfHw8iYmJfPnll4wfP57hw4fTtGnTKjsGUrFaBNTm2lYNMAxYEJtkdhwRkUq1YvsxvtmcjKuLhVdGdsLdtdyfmZNyKtcI180338xTTz3F9u3b6dChA+7u7oW+P3z48AoJ90effPIJjz76aMGnDocPH868efMKbbNr1y5sNlvB15MnT+bcuXM8/PDDpKWl0bNnT5YvX46v75UPo3p6evLFF18wffp0srKyCAkJYdy4cUyePLliXpiYZmyfUH7ZfZwvNxxi0oBW1PYs92dIRESqLdu5HJ6J3gLAuD81o0Njq8mJaiaLUY6VfF1cSm7GFouFvDzNb3Sp9PR0rFYrNpsNPz8/s+MIYLcb3Pj6avafOMv04e24v3eo2ZFERCrclEUJfL7hEGH+tfjhsT/h5e5qdiSHUlG/v8s1pmi320u8qWyJo3BxsRQs9zM/Jgm7vcx/e4iIVGtr957g8w35H9qafXtHlS0Tlalw3XTTTYXesnv55ZcLrWd48uRJ2rZtW2HhRCrbiK6N8fV0I/HEWVbvPm52HBGRCpOZncuUxQkA3BcRQo+weiYnqtnKVLiWLVtWaDLP2bNnF5oGIjc3t8SldkSqo9qebozqnr+G2IdrE01OIyJScV5ZtotDp87RqI43kweHmx2nxitT4br0cq9yXP4lUu3cHxGKxQJr9pxgb2rRhc1FRBzNbwdOMT8mCYAZIzroQ0HVgD4XKjVe0/o+3NgmfzLcKE2EKiIO7nxOHpMXJmAYcHvXxvRr1cDsSEIZC5fFYimyVI7WoRNnMPbCxfOLfz+CLTPH3DAiIldh3sq97Dt+Fv/anjw3tI3ZceSCMo0xGobBmDFjCmZQP3/+PA899BC1auWvNP7H67tEHElEs/qEB/myMyWDzzccZHw/zcIsIo5nW7KNt1fvA+ClW9tRx8fD5ERyUZlGuO6//34CAgKwWq1YrVbuvfdegoODC74OCAi44oWrRaoTi8VSMMr1UewBcvPs5gYSESmjnDw7kxcmkGc3uKlDEIPbl7z8nFS9Mo1wRUVFVVYOEdPd0rkRs37YyZHT5/hxxzH9sBIRh/LeL/vZlpyO1dud6cPbmx1HLqGL5kUu8HJ3ZXTP/PUxP9TF8yLiQPamnuGfP+0B4PmhbWng62lyIrmUCpfIH0T2CsXVxcL6xFNsS7ZdfgcREZPZ7QZTFiWQnWunX6sGjOjayOxIUgwVLpE/CLJ6MaR9EKApIkTEMXwUm8TGA2nU8nBlxogOmj2gmlLhErnE2D5hACyJT+bEGX3yVkSqr0OnMpmzLH+Flyk3taFRHW+TE0lJVLhELtG1aR06NbaSnWfn03UHzY4jIlIswzB4OnoLmdl59Airxz09mpodSUqhwiVyifwpIvJHuT6OO0B2rqaIEJHq56vfDrNmzwk83VyYfXtHXFz0VmJ1psIlUoybOjQkwNeT4xlZfL/lqNlxREQKSU0/z0vfbgfg8QGtCPOvZXIiuRwVLpFieLi5cG+vEACi1iZqoXYRqTYMw+DZr7eSfj6Xjo2tPNA3zOxIcgVUuERKMLpnUzxcXdh82MbvB0+bHUdEBIDvthxl+fZjuLlYmDOyI26u+lXuCPT/kkgJ/Gt7MrxzMADzY5LMDSMiAqSdzeaF/20D4OHrWhAe5GdyIrlSKlwipbi4vuIPW46SYjtvbhgRqfFe/HY7J89m0yqwNhOua2F2HCkDFS6RUrQLttIjrB65doOP45LMjiMiNdiqnalEbzqCiwXmjOyEh5t+hTsS/b8lchlje4cC8Om6g5zPyTM3jIjUSBnnc3g6egsAD/QNo3OTOuYGkjJT4RK5jAFtA2lUx5u0zBz+F3/E7DgiUgPN/GEnR23nCanvw+MDWpsdR8pBhUvkMtxcXbgv4uIUEUmaIkJEqlTsvpMFq17MGtERbw9XkxNJeahwiVyBu7o3xdvdlZ0pGcTuP2l2HBGpIc5l5zFlcQKQP1VNRPP6JieS8lLhErkCVh93RnRtBOSPcomIVIXXV+ziwMlMGlq9mDok3Ow4chVUuESu0MUpIn7ccYyDJzPNDSMiTi/+0Gk++DURgBm3dcDXy93kRHI1VLhErlCLAF/+1NIfw4AFsUlmxxERJ5aVm8fkhZuxG3Bbl0ZcFx5gdiS5SipcImXw5z75a5Z9ueEQZ7JyTU4jIs7qrVX72H3sDPVrefDc0LZmx5EKoMIlUgb9WjWgmX8tMrJyWfz7YbPjiIgT2pmSzr9X7QVg+i3tqFfLw+REUhFUuETKwMXFwv0XJkKdvzYJu11TRIhIxcnNszN5YQK5doOBbQO5uUNDsyNJBVHhEimj269pjK+nG/tPnGX1nuNmxxERJ/LBr4kkHLbh5+XGS7e2x2KxmB1JKogKl0gZ1fZ0445uTQBNESEiFWf/8TO8vmI3AM8ObUuAn5fJiaQiqXCJlMOY3qFYLPDL7uPsTT1jdhwRcXB2u8GUxVvIyrXzp5b+3HFNY7MjSQVT4RIph6b1fbghPBCA+TGJJqcREUf3yfqDrE88hY+HKzNu66C3Ep2QCpdIOf35wkSoi347gi0zx9wwIuKwjpw+x6zvdwAweVBrmtTzMTmRVAYVLpFyimhen9aBvpzLyeOLjQfNjiMiDsgwDJ5evIWz2Xl0C6nLfRGhZkeSSqLCJVJOFoulYLmfBTEHyM2zmxtIRBxO9KYjrN59HA83F2aP7IiLi95KdFYqXCJX4dYujajr486R0+f4cccxs+OIiAM5npHFi99uB+CxG1rSvEFtkxNJZVLhErkKXu6u3N2jKQAfaooIESmDF5Zs5XRmDu2C/fjLtc3MjiOVTIVL5CpFRoTg6mJhfeIptiXbzI4jIg5g6dajfL8lBTcXC3NGdsTdVb+OnZ3+Hxa5Sg2t3gxpHwTkL/cjIlKa05nZPPv1NgAe6tecdsFWkxNJVVDhEqkAFy+e/9/mZE6eyTI3jIhUay99t4MTZ7Jo3qAWE65vYXYcqSIqXCIVoGvTunRsbCU7186n6zRFhIgUb/Xu4yz87TAWC8wZ2Qkvd1ezI0kVUeESqQB/nCLi47gDZOdqiggRKexMVi5PL94C5C8Pdk1IXZMTSVVS4RKpIDd3CKaBryepGVn8sPWo2XFEpJqZs3QnR06fo0k9b54a1NrsOFLFVLhEKoiHmwv39gwBNEWEiBS2PvEUH8UeAGDWiI74eLiZnEiqmgqXSAUa3bMpHq4ubD50mt8PppkdR0SqgfM5efx9UQIAd3ZrQp8W/iYnEjOocIlUoAa+ngzrFAxAlEa5RASY++MeEk+cJdDPk6dvbmN2HDGJCpdIBbt48fwPW46SYjtvbhgRMdWWwzb+s2Y/AC/d2gGrt7vJicQsKlwiFax9Iys9QuuRazf4OC7J7DgiYpLsXDtPLdxMnt1gWKdgBrQNNDuSmEiFS6QSXBzl+nTdQc7n5JkbRkRM8c7qfexMyaCujzvThrU1O46YTIVLpBIMaBtIozrepGXmsCQ+2ew4IlLF9hzL4M2VewCYNrwd9Wt7mpxIzKbCJVIJ3FxdiIy4OEVEIoZhmJxIRKpKnt3gqYUJ5OQZ3BAewPALH6SRms1hCldaWhqRkZFYrVasViuRkZGcPn261H0Mw2DatGkEBwfj7e1N//792bZtW6Ft+vfvj8ViKXS76667rvq5Re7q3gQvdxd2pmQQt/+U2XFEpIpErU0k/tBpfD3dePm2DlgsFrMjSTXgMIVr9OjRxMfHs3TpUpYuXUp8fDyRkZGl7jNnzhxef/115s2bx4YNGwgKCmLAgAFkZGQU2m7cuHEcPXq04Pbuu+9e9XOL1PHxYETXxkD+D2ARcX4HTp7l1eW7AHj65jYEWb1MTiTVhUNMdbtjxw6WLl1KXFwcPXv2BOA///kPERER7Nq1i9atiy6RYBgGc+fO5ZlnnmHEiBEALFiwgMDAQD799FPGjx9fsK2Pjw9BQUEV9twiF43tHcqn6w6yYscxDp3KpEk9H7MjiUglMQyDKYu2cD7HTkSz+tzVvYnZkaQacYgRrtjYWKxWa0HhAejVqxdWq5WYmJhi90lMTCQlJYWBAwcW3Ofp6Um/fv2K7PPJJ5/g7+9Pu3btePLJJwuNgJXnuUUuahnoy59a+mMYsCAmyew4IlKJPt9wiNj9J/Fyd2HW7XorUQpziBGulJQUAgICitwfEBBASkpKifsABAYWnvckMDCQAwcOFHx9zz33EBYWRlBQEFu3bmXq1Kls3ryZFStWlPu5AbKyssjKyir4Oj09vZRXKM5sbJ9Q1uw5wRcbDzFpQCtqeTrEPzsRKYOjtnPM+G4HAE8ObE1I/VomJ5LqxtQRrmnTphW5YP3S28aNGwGK/UvBMIzL/gVx6fcv3WfcuHHceOONtG/fnrvuuouFCxfy448/8vvvv5f4GFfy3DNnziy4yN5qtdKkiYaWa6r+rQII869FxvlcFv1+2Ow4IlLBDMPg2eitZGTl0qVpHcb2CTM7klRDphauCRMmsGPHjlJv7du3JygoiGPHjhXZ//jx40VGsC66eE3WpaNQqampJe4D0LVrV9zd3dmzZ0/B45T1uQGmTp2KzWYruB06dKjEbcW5ubhYuP/CFBHz1yZht2uKCBFnsmRzMj/tTMXD1YU5t3fE1UVvJUpRpr634e/vj7//5VdNj4iIwGazsX79enr06AHAunXrsNls9O7du9h9Lr5NuGLFCrp06QJAdnY2q1evZvbs2SU+17Zt28jJyaFhw4blfm7Iv17M01MT3Um+kd2a8Nry3ew/cZbVe45zXeuib1OLiOM5eSaLaUvypxt65PoWtAz0NTmRVFcOcdF8mzZtGDx4MOPGjSMuLo64uDjGjRvH0KFDC31KMDw8nOjoaCD/bcCJEycyY8YMoqOj2bp1K2PGjMHHx4fRo0cDsG/fPl588UU2btxIUlIS33//PXfccQddunShT58+ZXpukdLU9nTjjm75byvPX5tkbhgRqTDTvtlOWmYO4UG+PNS/udlxpBpziMIF+Z8k7NChAwMHDmTgwIF07NiRjz/+uNA2u3btwmazFXw9efJkJk6cyMMPP0y3bt04cuQIy5cvx9c3/y8QDw8PfvrpJwYNGkTr1q159NFHGThwID/++COurq5lem6Ry7m/dwgWC6zefZy9qWfMjiMiV2nF9mN8szkZVxcLr4zshLurw/xKFRNYDK05UiXS09OxWq3YbDb8/PzMjiMmeXDBBn7ckUpkrxD+cWt7s+OISDnZzuUw4PXVpGZk8VC/5kwZEm52JKkkFfX7W3VcpApd/PTSot8PYzuXY3IaESmvGd/tIDUji2b+tZh4Y0uz44gDUOESqUK9m9endaAvmdl5fLlBn1wVcURr9+bPqwcwe2RHvNxdL7OHiAqXSJWyWCyM6RMKwILYJPI0RYSIQ8nMzmXK4gQA7osIoXtoPZMTiaNQ4RKpYrd2bkQdH3cOp51jxfaic7yJSPX1yrJdHDp1jkZ1vJk8WNdtyZVT4RKpYt4ertzdoykAUWsTTU4jIlfqtwOnmH9hTdSZIzpQW8t0SRmocImYILJXCK4uFtYlnmJbsu3yO4iIqc7n5DF5YQKGASOvacy1rRqYHUkcjAqXiAmC63gzuH3+8lOaCFWk+pu3ci/7jp/Fv7Ynz97cxuw44oBUuERM8ucLF8//b3MyJ89kmRtGREq0LdnG26v3AfDSre2o4+NhciJxRCpcIibp2rQuHRpZyc6189n6g2bHEZFi5OTZmbwwgTy7wU0dghjcvqHZkcRBqXCJmMRisTD2wijX/JgD7Diabm4gESnivV/2sy05nTo+7kwfrtUhpPxUuERMdHPHhjSt58OJM1kMn/cr81buITfPbnYsEQH2pp7hnz/tAeD5oW1p4OtpciJxZCpcIibydHNl4V8jGNA2kJw8g1eX7+a2t2LYfSzD7GgiNZrdbvD3RQlk59rp16oBt3VpZHYkcXAqXCImC/D14r3Ia3jjzk74ebmx5YiNof/6lbd/3qfRLhGTfBSbxG8H0qjl4cqMER2wWCxmRxIHp8IlUg1YLBZu69KYFY/34/rwALLz7MxeupOR78SyN/WM2fFEapRDpzKZs2wXAFNuakOjOt4mJxJnoMIlUo0E+nnxwf3deGVkR3w93Yg/dJqb/rWG//yyX+suilQBwzCYungLmdl59Airxz0XVoUQuVoqXCLVjMVi4Y5uTVg26VqubdWA7Fw7L3+/g1HvxpJ44qzZ8USc2lcbD/Pr3hN4urkw+/aOuLjorUSpGCpcItVUcB1vFoztzqwLa7b9diCNIf/8hQ9/TcSu0S6RCncs/Tz/+G47AI8PaEWYfy2TE4kzUeESqcYsFgt39WjK0ol/ok+L+pzPsfPit9u56z9xHDyZaXY8EadhGAbPfb2VjPO5dGxs5YG+YWZHEiejwiXiABrX9eG/D/TkH7e2x8fDlfWJpxj8z1/4ODZJo10iFeC7LUdZvv0Y7q4W5ozsiJurfj1KxdIZJeIgLBYLkb1CWPrYtfQMq0dmdh7P/W8b936wjkOnNNolUl4nz2Txwv+2AfBw/xaEB/mZnEickQqXiINpWt+Hz8b1Ytqwtni5uxCz7ySD5/7Cp+sOYhga7RK5Unl2g8/XH2TgG79w8mw2rQN9+dt1LcyOJU7KYugndJVIT0/HarVis9nw89NfT1Ixkk6c5cmvNrPxQBoAf2rpz+zbOxKseYNESrUh6RTTv9nG1iP5a5g2a1CLt++5htZBviYnk+qmon5/q3BVERUuqSx5doOotYm8smwXWbl2fD3deG5YW+64prFmxxa5xJHT55j1w06+2ZwMgK+XGxNvbMV9ESG467otKYYKl4NR4ZLKtjf1DE9+tZn4Q6cBuK51A2bd3pFAPy9zg4lUA+ey83j3l328s3of53PsWCxwd4+mPDGgFfVra1FqKZkKl4NR4ZKqkGc3+M+a/by+fDfZeXb8vNyYNrwdt3VppNEuqZEMw+C7LUeZ+f1Ojpw+B0CPsHq8MKwt7YKtJqcTR6DC5WBUuKQq7TmWwRNfbSbhsA2AG9sEMmNEewJ8NdolNcfWIzZe/GY765NOAdCojjdP39SGmzoE6Q8QuWIqXA5GhUuqWm6enXd/2c/cH3eTk2dQx8ed6cPbMbxTsH7ZiFM7eSaLV5fv4vMNhzAM8HJ34eH+LfjLtc3wcnc1O544GBUuB6PCJWbZcTSdJ7/azLbk/E9jDWkfxD9ubY+/rlsRJ5Oda+ej2CT++dMeMs7nAjC8UzBThoTrk7tSbipcDkaFS8yUk2fnrVX7eHPlHnLtBvVqefDSre25qUNDs6OJVIifd6Xy4rfb2X88f4H39o38eGFYO7qH1jM5mTg6FS4Ho8Il1cHWIzae/GozO1MyABjasSEv3tKeerU8TE4mUj77j5/hpe92sHJnKgD+tT14alBrRl7TBFcXvXUuV0+Fy8GocEl1kZ1r582Ve3jr533k2Q38a3vw8m0dGNQuyOxoIlcs/XwOb/60h6i1SeTaDdxdLYztE8aE61vg5+VudjxxIipcDkaFS6qbhMOneeLLzexJPQPArZ2DmTa8HXV8NNol1Vee3eCrjYd4ZdkuTp7NBuCG8ACeubkNzRrUNjmdOCMVLgejwiXV0fmcPP750x7eXb0PuwEBvp7MHNGBG9oEmh1NpIhLl+Np3qAWzw1tS//WASYnE2emwuVgVLikOtt0MI0nvtpccMHxyGsa89zQtli99daMmC/59DlmajkeMYkKl4NR4ZLq7nxOHq8t38X7vyZiGBDk58Ws2zto9EBMcy47j/d+2c/bq/cWLMdzV/emPDlQy/FI1VHhcjAqXOIoNiad4smvNpN0MhOAu7o34Zmb2+CrC5Glimg5HqlOVLgcjAqXOJJz2XnMWbaTqLVJQP6SKLNv70jflv7mBhOnp+V4pLpR4XIwKlziiOL2n2TywgQOnsof7bq3V1OmDmlDLU83k5OJsyluOZ6/9stfjsfbQ8vxiHlUuByMCpc4qrNZucz6YScfxx0AoHFdb14Z2YmI5vVNTibOQMvxSHWnwuVgVLjE0a3de4LJCxMKrqkZ0zuUyYNb4+Oh0S4pHy3HI45AhcvBqHCJM8g4n8OM73fy2fqDAITU9+HVOzrpF6SUiZbjEUeiwuVgVLjEmfyy+zh/X5TAUdt5LBb4c58wnhrUGi93XWsjJbu4HM/8mCRy8rQcjzgGFS4Ho8Ilzib9fA4vfbudLzceBqCZfy1euaMT14TUNTmZVDd5doOFv+Uvx3PiTP5yPNeHB/CsluMRB6DC5WBUuMRZrdqZypTFCRxLz8LFAuOubcakG1tptEuAosvxNLuwHM91mlBXHIQKl4NR4RJnZsvMYfq321j8+xEAWgTU5rU7OtGpSR1zg4lptByPOAsVLgejwiU1wYrtx5i6eAsnzmTh6mLhoX7NePSGlni6abSrptByPOJsVLgcjAqX1BRpZ7N5Yck2llwY2Wgd6MtrozrRvpGWZHFmWo5HnJUKl4NR4ZKa5octR3n2662cPJuNm4uFv13Xgr9d1wIPN72d5Gy2JduY/s121idqOR5xPipcDkaFS2qik2eyeO5/W/l+SwoAbRv68eodnWgbrH8DziB/OZ7dfL7hoJbjEaelwuVgVLikJvs2IZnnvt5KWmYO7q4WHr2+JQ/1b66Lpx2UluORmkSFy8GocElNdzwji2eit7B8+zEAOjSy8tqoTrQK9DU5mZTFpcvxtAv2Y9pwLccjzkuFy8GocInkX1j9v/hkXliyDdu5HDxcXZg0oBXj/hSGm0a7qjUtxyM1lQqXg1HhEvl/x9LP8/TiLfx04Zd35yZ1ePWOTrQI0Kzj1Y2W45GaToXLwahwiRRmGAaLfj/C9G+2kXE+Fw83F54a2Jo/9w3TiEk1oOV4RPKpcDkYFS6R4h21nePvi7bwy+7jAFwTUpdX7+hEmH8tk5PVXFqOR+T/VdTvb4e5aCItLY3IyEisVitWq5XIyEhOnz5d6j6GYTBt2jSCg4Px9vamf//+bNu2rdA2/fv3x2KxFLrdddddhbYJDQ0tss2UKVMq+iWK1EgNrd4sGNudWSM6UNvTjd8OpDHkn7/w4a+J2O36e7AqJZ8+xyOfbeKOd2LZeiQdXy83nr25DcsmXquyJXKVHGaEa8iQIRw+fJj33nsPgL/85S+EhobyzTfflLjP7Nmzefnll5k/fz6tWrXipZde4pdffmHXrl34+uZ/Mqp///60atWKF198sWA/b29vrNb/nxk5NDSUBx54gHHjxhXcV7t2bWrXvvJhdY1wiVze4bRM/r4ogbV7TwL5M5W/MrIjIfU12lWZSlqO54mBrfDXcjxSw1XU72+3CsxUaXbs2MHSpUuJi4ujZ8+eAPznP/8hIiKCXbt20bp16yL7GIbB3LlzeeaZZxgxYgQACxYsIDAwkE8//ZTx48cXbOvj40NQUFCpGXx9fS+7jYhcncZ1ffjvAz35ZN1BZny/g/WJpxg8dw1Tbwrn3p4huOjargpV7HI8ofV4flhbLcUkUsEc4i3F2NhYrFZrQdkC6NWrF1arlZiYmGL3SUxMJCUlhYEDBxbc5+npSb9+/Yrs88knn+Dv70+7du148sknycjIKPJ4s2fPpn79+nTu3JmXX36Z7OzsUjNnZWWRnp5e6CYil2exWLi3VwjLJl5Lr2b1OJeTx/P/28a9H6zj0KlMs+M5jW3JNu58L44Jn27iyOlzBFu9mDe6C1+M76WyJVIJHGKEKyUlhYCAotcPBAQEkJKSUuI+AIGBgYXuDwwM5MCBAwVf33PPPYSFhREUFMTWrVuZOnUqmzdvZsWKFQXbPPbYY3Tt2pW6deuyfv16pk6dSmJiIu+//36JmWfOnMn06dPL9DpF5P81qefDpw/24uO4A8z6YScx+04yeO4vPHNzW+7u0URr9JWTluMRMYephWvatGmXLSUbNmwAKPaHq2EYl/2he+n3L93nj9dltW/fnpYtW9KtWzd+//13unbtCsCkSZMKtunYsSN169Zl5MiRBaNexZk6dSqPP/54wdfp6ek0adKk1KwiUpiLi4X7e4fSr1UDnlq4mQ1JaTwdvYUfth5l9u0dtYxMGeTk2fko9gBzf9yt5XhETGBq4ZowYUKRTwReKjQ0lISEBI4dO1bke8ePHy8ygnXRxeutUlJSaNiwYcH9qampJe4D0LVrV9zd3dmzZ09B4bpUr169ANi7d2+JhcvT0xNPT11sKlIRQv1r8flfIpgfk8ScpTtZs+cEg974heeGtuWObo012nUZP+9K5R/fbmefluMRMY2phcvf3x9/f//LbhcREYHNZmP9+vX06NEDgHXr1mGz2ejdu3ex+1x8m3DFihV06dIFgOzsbFavXs3s2bNLfK5t27aRk5NTqKRdatOmTQClbiMiFcvVxcIDfcPo37oBT361mU0HTzN5UQI/bD3KzBEdCbJ6mR2x2tFyPCLVh0NNC5GcnMy7774L5E8LERISUmhaiPDwcGbOnMltt90G5F/oPnPmTKKiomjZsiUzZszg559/LpgWYt++fXzyySfcdNNN+Pv7s337dp544gm8vb3ZsGEDrq6uxMbGEhcXx3XXXYfVamXDhg1MmjSJbt268b///e+K82taCJGKk2c3eH/Nfl5bsZvsXDt+Xm5MG96O27o00mgXRZfjcXOxMLZPKI/c0FLL8YiUUY2aFgLyP0n46KOPFnzqcPjw4cybN6/QNrt27cJmsxV8PXnyZM6dO8fDDz9MWloaPXv2ZPny5QVzcHl4ePDTTz/xz3/+kzNnztCkSRNuvvlmXnjhBVxd8y8e9fT05IsvvmD69OlkZWUREhLCuHHjmDx5chW9chG5lKuLhfH9mnN9eABPfrWZzYdtPP7lZr7fksKMEe0J8K2Zo10lLcfzzM1taK7leERM5TAjXI5OI1wilSM3z867v+xn7o+7yckzqOPjzvTh7RjeKbhGjXZpOR6RyqG1FB2MCpdI5dqZks4TX25mW3J+4RjcLoiXbmvv9DOlJ58+x6wfdrJkczIAvl5uPHZDS+6LCMXDzSGmWhSp1lS4HIwKl0jly8mz89aqfby5cg+5doN6tTz4xy3tubmj833ARcvxiFQNFS4Ho8IlUnW2Jdt44svN7EzJXzViaMeGvHhLe+rV8jA52dXTcjwiVUuFy8GocIlUrexcO/NW7uHfP+8jz27gX9uDl2/rwKB2jrsm6rZkG9O/2c76xFMABFu9ePrmNtzcoWGNul5NpCqpcDkYFS4Rc2w5bOOJr+LZfewMALd2Dmba8HbU8XGc0S4txyNiHhUuB6PCJWKerNw85v64h3dX78NuQANfT2aN6MANbUpedaI6KG45nmEXluNppOV4RKqECpeDUeESMd+mg2k8+dXmgiVubu/amOeHtcXqXf0mAy1uOZ4XhrWjR5iW4xGpSipcDkaFS6R6OJ+Tx+srdvOfNfsxDAjy82LW7R3oX03mq7p0OZ76tfKX47mjm5bjETGDCpeDUeESqV5+O3CKJ79KIPFE/gjSXd2b8MzNbfA1aemb9PM5zFu5l6i1iVqOR6QaUeFyMCpcItXPuew8Xlm2i6iYRAwj/1N/c0Z2om9L/yrLoOV4RKo3FS4Ho8IlUn2t23+SpxYmcPBUJgD39GzK1JvaUNuzcpeb1XI8ItWfCpeDUeESqd4ys3OZ/cNOFsQeAKBxXW/mjOxI7+YVP9pVZDkeTzceu1HL8YhURypcDkaFS8QxxOw9wVMLEwpmcR/TO5TJg1vj43H1o13nc/J4d/Wly/E04YmBrbUcj0g1pcLlYFS4RBzHmaxcZny/g0/XHQQgpL4Pr4zsVO4pGQzD4PstKcz4foeW4xFxMCpcDkaFS8Tx/LL7OFMWJZBsO4/FAn/uE8aTA1uXaXb34pbjmXpTG4Z21HI8Io5AhcvBqHCJOKb08zm8/O0Ovth4CIBm/rV45Y5OXBNSt9T9iluO56F+zRl/bXMtxyPiQFS4HIwKl4hjW7UrlSmLEjiWnoWLBcb9qRmTBrTCy71wedJyPCLORYXLwahwiTg+W2YO07/dxuLfjwDQIqA2r97Ric5N6gBajkfEGalwORgVLhHnsWL7MZ6O3sLxjPzRrgf/1Iy9qWe0HI+IE1LhcjAqXCLOJe1sNtO+2cb/4pML7tNyPCLOp6J+f1fuNMoiIk6qbi0P/nlXF4a0D2LmDztpGVCbqTdpOR4RKZ4Kl4jIVRjcviGD2zc0O4aIVHNaQ0JERESkkqlwiYiIiFQyFS4RERGRSqbCJSIiIlLJVLhEREREKpkKl4iIiEglU+ESERERqWQqXCIiIiKVTIVLREREpJKpcImIiIhUMhUuERERkUqmwiUiIiJSyVS4RERERCqZCpeIiIhIJXMzO0BNYRgGAOnp6SYnERERkSt18ff2xd/j5aXCVUUyMjIAaNKkiclJREREpKwyMjKwWq3l3t9iXG1lkytit9tJTk7G19eXjIwMmjRpwqFDh/Dz8zM7WrWQnp6uY3IJHZPi6bgUpWNSlI5J8XRcirrcMTEMg4yMDIKDg3FxKf+VWBrhqiIuLi40btwYAIvFAoCfn59O+EvomBSlY1I8HZeidEyK0jEpno5LUaUdk6sZ2bpIF82LiIiIVDIVLhEREZFKpsJlAk9PT1544QU8PT3NjlJt6JgUpWNSPB2XonRMitIxKZ6OS1FVdUx00byIiIhIJdMIl4iIiEglU+ESERERqWQqXCIiIiKVTIXrKv3yyy8MGzaM4OBgLBYLX3/9dYnbjh8/HovFwty5c0t9zPnz52OxWIrczp8/X7HhK8nljsmYMWOKvLZevXpd9nEXLVpE27Zt8fT0pG3btkRHR1fSK6h4lXFMHP08gSv797Njxw6GDx+O1WrF19eXXr16cfDgwVIf15nPFSj7MXH0c+Vyx6S412axWHjllVdKfVxHPk+gco6Ls58rZ86cYcKECTRu3Bhvb2/atGnD22+/fdnHrYhzRYXrKp09e5ZOnToxb968Urf7+uuvWbduHcHBwVf0uH5+fhw9erTQzcvLqyIiV7orOSaDBw8u9Nq+//77Uh8zNjaWO++8k8jISDZv3kxkZCSjRo1i3bp1FR2/UlTGMQHHPk/g8sdl37599O3bl/DwcH7++Wc2b97Mc889V+prdPZzpTzHBBz7XLncMbn0dX344YdYLBZuv/32Eh/T0c8TqJzjAs59rkyaNImlS5fy3//+lx07djBp0iQeeeQR/ve//5X4mBV2rhhSYQAjOjq6yP2HDx82GjVqZGzdutUICQkx3njjjVIfJyoqyrBarZWSsaoVd0zuv/9+45ZbbinT44waNcoYPHhwofsGDRpk3HXXXVeZsOpV1DFxpvPEMIo/Lnfeeadx7733lulxnP1cKc8xcaZzpaSfs390yy23GNdff32p2zjTeWIYFXdcnP1cadeunfHiiy8Wuq9r167Gs88+W+LjVNS5ohGuSma324mMjOSpp56iXbt2V7zfmTNnCAkJoXHjxgwdOpRNmzZVYsqq9/PPPxMQEECrVq0YN24cqamppW4fGxvLwIEDC903aNAgYmJiKjNmlSrrMQHnPk/sdjvfffcdrVq1YtCgQQQEBNCzZ89S37YH5z5XyntMwLnPlT86duwY3333HQ888ECp2znzeVKcKz0u4NznSt++fVmyZAlHjhzBMAxWrVrF7t27GTRoUIn7VNS5osJVyWbPno2bmxuPPvroFe8THh7O/PnzWbJkCZ999hleXl706dOHPXv2VGLSqjNkyBA++eQTVq5cyWuvvcaGDRu4/vrrycrKKnGflJQUAgMDC90XGBhISkpKZcetEuU5Js5+nqSmpnLmzBlmzZrF4MGDWb58ObfddhsjRoxg9erVJe7nzOdKeY+Js58rf7RgwQJ8fX0ZMWJEqds583lSnCs9Ls5+rvzrX/+ibdu2NG7cGA8PDwYPHsxbb71F3759S9ynws6VMo2HSam4ZPhy48aNRmBgoHHkyJGC+67kLcVL5eXlGZ06dTIeeeSRCkpadS49JsVJTk423N3djUWLFpW4jbu7u/Hpp58Wuu+///2v4enpWRExq1RFHZNLOfJ5YhhFj8uRI0cMwLj77rsLbTds2LBSh/Kd+Vwp7zG5lCOfK5f799O6dWtjwoQJl30cZzpPDKPijsulnO1ceeWVV4xWrVoZS5YsMTZv3my8+eabRu3atY0VK1aU+DgVda64la2eSVmsWbOG1NRUmjZtWnBfXl4eTzzxBHPnziUpKemKHsfFxYXu3bs7zV8Yl2rYsCEhISGlvr6goKAif02kpqYW+avDWVzJMbmUs50n/v7+uLm50bZt20L3t2nThl9//bXE/Zz5XCnvMbmUs50rF61Zs4Zdu3bxxRdfXHZbZz5PLlWW43IpZzpXzp07x9NPP010dDQ333wzAB07diQ+Pp5XX32VG2+8sdj9Kupc0VuKlSgyMpKEhATi4+MLbsHBwTz11FMsW7bsih/HMAzi4+Np2LBhJaY1z8mTJzl06FCpry8iIoIVK1YUum/58uX07t27suOZ4kqOyaWc7Tzx8PCge/fu7Nq1q9D9u3fvJiQkpMT9nPlcKe8xuZSznSsXffDBB1xzzTV06tTpsts683lyqbIcl0s507mSk5NDTk4OLi6Fq4+rqyt2u73E/SrsXCnTeJgUkZGRYWzatMnYtGmTARivv/66sWnTJuPAgQPFbl/cW4qRkZHGlClTCr6eNm2asXTpUmPfvn3Gpk2bjLFjxxpubm7GunXrKvOlVJjSjklGRobxxBNPGDExMUZiYqKxatUqIyIiwmjUqJGRnp5e8BiXHpO1a9carq6uxqxZs4wdO3YYs2bNMtzc3Iy4uDgzXmKZVcYxcfTzxDAu/+9n8eLFhru7u/Hee+8Ze/bsMd58803D1dXVWLNmTcFj1KRzxTDKd0wc/Vy5kp+zNpvN8PHxMd5+++1iH8PZzhPDqJzj4uznSr9+/Yx27doZq1atMvbv329ERUUZXl5exltvvVXwGJV1rqhwXaVVq1YZQJHb/fffX+z2xRWufv36Fdp+4sSJRtOmTQ0PDw+jQYMGxsCBA42YmJjKexEVrLRjkpmZaQwcONBo0KCB4e7ubjRt2tS4//77jYMHDxZ6jEuPiWEYxldffWW0bt3acHd3N8LDw8t0fZPZKuOYOPp5YhhX9u/ngw8+MFq0aGF4eXkZnTp1Mr7++utCj1GTzpWLynpMHP1cuZJj8u677xre3t7G6dOni30MZztPDKNyjouznytHjx41xowZYwQHBxteXl5G69atjddee82w2+0Fj1FZ54rFMAyjbGNiIiIiIlIWuoZLREREpJKpcImIiIhUMhUuERERkUqmwiUiIiJSyVS4RERERCqZCpeIiIhIJVPhEhEREalkKlwiIiIilUyFS0TkMkJDQ5k7d67ZMUTEgalwiYhTGzZsGDfeeGOx34uNjcVisfD7779XcSoRqWlUuETEqT3wwAOsXLmSAwcOFPnehx9+SOfOnenatasJyUSkJlHhEhGnNnToUAICApg/f36h+zMzM/niiy944IEHWLRoEe3atcPT05PQ0FBee+21Eh8vKSkJi8VCfHx8wX2nT5/GYrHw888/A/Dzzz9jsVhYtmwZXbp0wdvbm+uvv57U1FR++OEH2rRpg5+fH3fffTeZmZkFj2MYBnPmzKFZs2Z4e3vTqVMnFi5cWJGHQ0RMosIlIk7Nzc2N++67j/nz52MYRsH9X331FdnZ2URERDBq1CjuuusutmzZwrRp03juueeKFLTymDZtGvPmzSMmJoZDhw4xatQo5s6dy6effsp3333HihUrePPNNwu2f/bZZ4mKiuLtt99m27ZtTJo0iXvvvZfVq1dfdRYRMZfF+ONPIBERJ7Rz507atGnDypUrue666wDo168fjRo1wmKxcPz4cZYvX16w/eTJk/nuu+/Ytm0bkH/R/MSJE5k4cSJJSUmEhYWxadMmOnfuDOSPcNWtW5dVq1bRv39/fv75Z6677jp+/PFHbrjhBgBmzZrF1KlT2bdvH82aNQPgoYceIikpiaVLl3L27Fn8/f1ZuXIlERERBVkefPBBMjMz+fTTT6viUIlIJdEIl4g4vfDwcHr37s2HH34IwL59+1izZg1//vOf2bFjB3369Cm0fZ8+fdizZw95eXlX9bwdO3Ys+O/AwEB8fHwKytbF+1JTUwHYvn0758+fZ8CAAdSuXbvg9tFHH7Fv376ryiEi5nMzO4CISFV44IEHmDBhAv/+97+JiooiJCSEG264AcMwsFgshbYtbeDfxcWlyDY5OTnFbuvu7l7w3xaLpdDXF++z2+0ABf/73Xff0ahRo0LbeXp6Xu7liUg1pxEuEakRRo0ahaurK59++ikLFixg7NixWCwW2rZty6+//lpo25iYGFq1aoWrq2uRx2nQoAEAR48eLbjvjxfQl1fbtm3x9PTk4MGDtGjRotCtSZMmV/34ImIujXCJSI1Qu3Zt7rzzTp5++mlsNhtjxowB4IknnqB79+784x//4M477yQ2NpZ58+bx1ltvFfs43t7e9OrVi1mzZhEaGsqJEyd49tlnrzqfr68vTz75JJMmTcJut9O3b1/S09OJiYmhdu3a3H///Vf9HCJiHo1wiUiN8cADD5CWlsaNN95I06ZNAejatStffvkln3/+Oe3bt+f555/nxRdfLChkxfnwww/JycmhW7duPPbYY7z00ksVku8f//gHzz//PDNnzqRNmzYMGjSIb775hrCwsAp5fBExjz6lKCIiIlLJNMIlIiIiUslUuEREREQqmQqXiIiISCVT4RIRERGpZCpcIiIiIpVMhUtERESkkqlwiYiIiFQyFS4RERGRSqbCJSIiIlLJVLhEREREKpkKl4iIiEglU+ESERERqWT/B4iAHI6NCodDAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "volumes = [structure.get_volume() for structure in structures]\n", + "\n", + "plt.plot(volumes, energies)\n", + "plt.xlabel('Volume')\n", + "plt.ylabel('Energy')\n", + "plt.savefig('evcurve.png')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Submission to HPC - up-scaling for high throughput screening " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "AiiDA allows registering remote HPC compute resources in its database, which can then be used for submitting simulations. Registration can be achieved via the CLI with `verdi computer setup`, followed by `verdi computer configure`, or via the Python API (see below). With the first command, the database entry for the `Computer` is created (and the provided information immutably stored in the DB), while the second command _configures_ the compute resource, which are, among others, the SSH connection parameters.\n", + "\n", + "AiiDA provides interfaces for various _transport technologies_, most notably SSH, implemented, as well as the most common schedulers, such as SLURM or SGE. This means that after successful registration in the database, using such HPC resources is as simple as changing a single string in the submission script that specifies the computer where the calculation should be run. When a simulation is started, AiiDA takes care to communicate with the HPC and upload all required files, as well as generate the submission script. In addition to the remote compute resource, also the codes (executables) have to be registered to such that they can then be used to run simulations.\n", + "\n", + "One advantage of this design, among others, is the fact that utilizing different HPCs, even within one workflow, can be easily achieved by using the corresponding identifiers of the compute resources. In addition, the software environment does _not_ have to be duplicated on the HPC, but AiiDA instead runs locally. However, this means that compute resources first have to be correctly registered in AiiDA, which can present an initial barrier for new users. The following code snippets show how a computer and a code can be registered via the Python API:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The `get_or_create` method allows keeping the `Computer` generation atomic within this notebook\n", + "from aiida import orm\n", + "\n", + "created, computer = orm.Computer.collection.get_or_create(\n", + " label=\"Todi\",\n", + " description=\"New Todi HPC at CSCS\",\n", + " hostname=\"todi\",\n", + " workdir=\"/scratch//aiida\",\n", + " transport_type=\"core.ssh\",\n", + " scheduler_type=\"core.slurm\",\n", + ")\n", + "\n", + "if created:\n", + " computer.store()\n", + " computer.set_minimum_job_poll_interval(10.0)\n", + " computer.set_default_mpiprocs_per_machine(128)\n", + " computer.set_append_text('')\n", + " computer.set_prepend_text('')\n", + " computer.configure()\n", + "\n", + "print(computer)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Register a code for later use \n", + "try:\n", + " code = aiida.orm.load_code(\"pw.x@todi\")\n", + "except aiida.common.NotExistent:\n", + " code = aiida.orm.Code(\n", + " computer=load_computer('todi'),\n", + " remote_computer_exec=[computer, data.pwx_path],\n", + " )\n", + " code.label = \"pw.x\"\n", + " code.description = \"Quantum ESPRESSO pw.x code\"\n", + " code.set_prepend_text(\"export OMP_NUM_THREADS=1\")\n", + " code.store()\n", + "code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# For `aiida-shell` the computer can directly be passed to the metadata options of the simulation\n", + "# This is because the executable is directly passed as the first argument to the `launch_shell_job` function\n", + "# And the path is resolved internally and a `Code` entity created\n", + "metadata={\n", + " 'options': {\n", + " 'prepend_text': f'export ESPRESSO_PSEUDO {FILEPATH_PSEUDOS.as_posix()}',\n", + " 'redirect_stderr': True\n", + " 'computer': load_computer(''),\n", + " }\n", + " }" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Cache Workflow Steps\n", + "\n", + "AiiDA provides caching functionality out of the box. \n", + "By default, caching is turned off, however, it can be globally enabled via the CLI with the command:\n", + "`verdi config set caching.default_enabled True` \n", + "as well as for individual simulation types with, e.g.\n", + "`verdi config set caching.enabled_for aiida.calculations:quantumespresso.pw`.\n", + "\n", + "To achieve this, all nodes are automatically hashed, and the hashes saved in the _extras_ of the respective Nodes once they are stored in the database:\n", + "```python\n", + "from aiida import orm\n", + "node = orm.Int(1)\n", + "node.store()\n", + "\n", + "print(node.base.caching.get_hash())\n", + "```\n", + "\n", + "Caching is then made possible by comparing the hashes of the `ProcessNode` of a given simulation that is to be executed, including all its input `Data` nodes, and comparing those with the existing entries in the DB. If the exact same simulation has already been carried out successfully, the corresponding nodes are duplicated and inserted into the workflow, rather than running the same simulation again. While this saves computational time, it does not necessarily save space, due to the duplication of data.\n", + "\n", + "There are a few details to keep in mind for caching, because, as the saying goes _\"There are only two hard things in Computer Science: cache invalidation and naming things.\"_ For one, a choice has to be made which data should be included when calculating the hash. E.g., which input nodes should be considered for calculating the hash of a Process? Or to which precision should the AiiDA, as well as the plugin version be included for evaluating the hash? We can assume that results between different major versions actually _are_ different, but what about minor and patch versions? Overall, a balance has to be found between minimizing the chance of false negatives, where two nodes should have the same hash but do not, and false positives, where two different nodes get the same hash by mistake. For instance, if every patch version would lead to different hashes, despite the fact that results can be expected to be the same, the feature would not be very useful. Furthermore, it is not immediately clear how far a workflow should be traced to generate the hash of a given Data node that it produced." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Keep Track of the Execution\n", + "\n", + "As provenance is one of the principal features behind AiiDA, great care has been taken during its development to ensure that the execution of a workflow is being tracked and that data can always be traced back. In addition, AiiDA implements checkpointing, meaning that each successful step completed during the execution of a workflow is recorded in the database. This allows restarting a running workflow, even after the workstation that runs AiiDA has been shut down or an SSH connection to an HPC has dropped." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Data Organisation\n", + "\n", + "AiiDA stores data in an SQL database, as well as an internal file repository. The former is used for input/parsed output Data nodes of well-defined types, such as `str` or `int` (using the corresponding AiiDA data types), while raw files are stored in the file repository (as they are often too large to store the contents in the SQL database).\n", + "\n", + "To organize data, AiiDA provides the concepts of _Groups_ that can be used to create collections of specific data nodes, such as, for example, all crystal structures used for a study, or all band structure simulations (but groups with mixed types are of course possible). As AiiDA further stores plenty of useful metadata (such as the creation and last modification time), one can easily query and/or filter data even years after its creation. Without a computational infrastructure like AiiDA, it can be otherwise difficult to navigate files and folders. Of course, also with AiiDA, proper organization using groups, and, possibly, adding additional metadata to the Nodes is just as helpful as keeping a logical directory structure when working with files and folders." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Sharing Workflows\n", + "\n", + "Lastly, once all required simulations have been conducted AiiDA, thus concluding a research project, all data can be exported using the `verdi archive create` CLI endpoint. Here, collections of data nodes, groups, or other entities can be selected for exporting (in addition to _all_ profile data). These archives can then be shared with other researchers, who can import them into their AiiDA profiles.\n", + "\n", + "While the data of a concrete, executed workflow is contained in the exported archives, the workflow logic cannot be shared in the same way. Instead, however, it is defined by the Python source code that was used to set up the workflow. These workflows are usually grouped together (with other required code infrastructure), in AiiDA-plugins, of which an extensive collection is available to the general public in the [AiiDA plugin registry](https://aiidateam.github.io/aiida-registry/)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "b" + ] + } + ], + "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.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/functions.py b/functions.py new file mode 100644 index 0000000..061ea40 --- /dev/null +++ b/functions.py @@ -0,0 +1,74 @@ +import os + +from ase import Atoms +from ase.io import write +from adis_tools.parsers import parse_pw +import matplotlib.pyplot as plt + + +def generate_structures( + structure: Atoms, + strain_lst: list + ) -> list[Atoms]: + + structure_lst = [] + + for strain in strain_lst: + structure_strain = structure.copy() + structure_strain.set_cell( + structure_strain.cell * strain**(1/3), + scale_atoms=True + ) + structure_lst.append(structure_strain) + + return structure_lst + + +def append_to_list(lst: list, item: float): + lst.append(item) + + +def split_string(string: str, character: str) -> list: + return string.split(character) + + +def write_input(input_dict, working_directory=".", return_string=False): + + filename = os.path.join(working_directory, 'input.pwi') + + os.makedirs(working_directory, exist_ok=True) + + write( + filename=filename, + images=input_dict["structure"], + Crystal=True, + kpts=input_dict["kpts"], + input_data={ + 'calculation': input_dict["calculation"], + 'occupations': 'smearing', + 'degauss': input_dict["smearing"], + }, + pseudopotentials=input_dict["pseudopotentials"], + tstress=True, + tprnfor=True + ) + + if return_string: + with open(filename) as f: + return f.read() + + +def collect_output(working_directory="."): + output = parse_pw(os.path.join(working_directory, 'pwscf.xml')) + return { + "structure": output['ase_structure'], + "energy": output["energy"], + "volume": output['ase_structure'].get_volume(), + } + + +def plot_energy_volume_curve(volume_lst, energy_lst): + plt.plot(volume_lst, energy_lst) + plt.xlabel("Volume") + plt.ylabel("Energy") + plt.savefig("evcurve.png")