04_Calculations_for_Crystals.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "05aafea4",
   "metadata": {},
   "source": [
    "# Crystals\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01b04a84",
   "metadata": {},
   "source": [
    "Crystalline solids have multiple peculiarities connected with their fundamental property: discrete translational symmetry. Additionally, they exhibit additional symmetries which influence their properties (e.g. electronic band structure) as well as procedures used in the calculation. Dealing with crystals require use of Periodic Boundary Conditions and, often, reciprocal space, Brillouin Zone etc. ASE provides multiple tools for working with crystals. We will use some of them in the following examples - which is modeled after one of the AbInit tutorials.\n",
    "\n",
    "In the following exercise we will:\n",
    "* build a crystal of Silicon,\n",
    "* calculate its lattice parameters,\n",
    "* derive its electronic band structure.\n",
    "\n",
    "As usual the detailed technical documentation regarding ASE toolbox can be found [ASE documentation site](https://wiki.fysik.dtu.dk/ase/index.html). The documentation concerning the DFT tool used in the exercise is located at [AbInit documentation site](https://docs.abinit.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "13859ea9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "%pylab is deprecated, use %matplotlib inline and import the required libraries.\n",
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from ase.build import bulk\n",
    "from ase.visualize import view\n",
    "from ase.visualize.plot import plot_atoms\n",
    "from ase.calculators.abinit import Abinit\n",
    "import ase.io\n",
    "from ase import units as un\n",
    "import os\n",
    "\n",
    "# Patch bug in ASE AbInit interface\n",
    "import ase_patch\n",
    "ase.io.abinit.read_abinit_out = ase_patch.read_abinit_out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "297d0b00",
   "metadata": {},
   "source": [
    "### Building the Silicon crystal\n",
    "\n",
    "The silicon has a simple crystal structure with two atoms per primitive unit cell - in fact it is the same structure as diamond. Thus, building the structure is almost trivial - the only free parameter is the length of the lattice vector $a$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7d1d692f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.81837662  3.81837662  3.81837662 60.         60.         60.        ]\n"
     ]
    }
   ],
   "source": [
    "Si = bulk('Si', crystalstructure='diamond', a=5.4)\n",
    "print(Si.cell.cellpar())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4216bde8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAGFCAYAAAD3kehhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABRvElEQVR4nO3dd1RU19oG8IfeBFFAiiCDhWIDBLFHEzVqVCxRY9TEEjX2XmOPDXshltiNLcaG2Dv2ygCiSFHpvTPUYcr3h5EvXEFBz8ye8v7WustcmDnnQfRxc87eZ2tIpVIpCCGEqBxN1gEIIYTIBhU8IYSoKCp4QghRUVTwhBCioqjgCSFERVHBE0KIiqKCJ4QQFUUFTwghKooKnhBCVBQVPCGEqCgqeEIIUVFU8IQQoqKo4AkhREVRwRNCiIqigieEEBVFBU8IISpKm3UAQghRZxKJBElJSYiJicGbN29gY2ODrl27cnJsKnhCCJGzVatWISAgADExMYiNjYVQKCz7nKamJp4+fYoWLVp88Xmo4Akh5AtIpVJkZWUhJiYG0dHRiImJQU5ODlasWFHpe7Kzs2FsbAxvb2/weDw4ODigdu3a6NatGwQCARo0aMBJNg3ak5UQQqrv7NmzWLx4MaKjoyEQCMo+XqNGDTRs2BCBgYHQ1Kz6bc7p06dj27ZtsLGxQUxMDCcZaQRPCFF7hYWFiImJKTcKj46OxoIFC+Du7l7heywsLNCuXTsMGzasbBT+fiSuoaFRrfOHhYXB19cXdnZ2aNmyJRdfEgAqeEKIGsvLy0OjRo2QlpZW9jEdHR3Y29uDx+OhuLi40ve2bdsWbdu2/eIMUqkUU6ZMgYODA1JTUzm59v4eFTwhROmJRCIkJCQgOjq63Ag8JiYG7dq1g4+PT4XvMzY2xpQpU2BnZwcHBwfweDzY2NhAS0tLbtnPnDmDGzdu4OzZs6hRowYcHBw4OzZdgyeEKDypVPrRyx49e/bExYsXy/6/jY1N2WWTLl26YMSIEXJIWX0SiQQuLi5o2LAhLly4wPnxqeAJIQohLy8PYWFhH1wHj4mJQVpaGrKzsyst+bt376KwsBAODg6oV68e9PX15Zz+80VGRkJHR4fTkft7VPCEEJmTSqXIzs6Gvr4+DA0NK3zNsWPHMGTIEABArVq1yi6ZvP/1119/hY6OjjxjKz0qeEIIZ0QiES5cuPDBKPz9VMJ//vkHAwcOrPC9mZmZSEhIAI/HQ82aNeWcXDVRwRNCquT9VEIAaNy4cYWvEYvFMDAwgJaWVrnR9/tf27VrBxsbGzmmVg4JCQmYPn06fHx8OFvkBFDBE0IqcOHCBdy7d6/cKPz9VMJevXrh3Llzlb43PT0d5ubm1Z4Lrs78/PzQr18/JCYmcvoPIE2TJERN/HcqYXx8PH7++edKX+vn54dr166Bx+PBxcUFPXr0KFvI86kRpoWFBdfRVR6fz0edOnVgbW3N6XFpBE+IigoODsbWrVvLRuHx8fEQi8UAAA0NDWRnZ1d6rftT0xLJ58nNzYWuri4MDAzKfbxXr14Qi8W4dOkSp+ejETwhSkIqlSItLa3cZZPvvvsOrq6uFb6+sLAQL1++hIODA1q1alU2AufxeLC3t4eenl6l56Jyl43Zs2fj2bNnCAwMLPd7HBQUJJO5+lTwhCgwiUSCPn364M2bN4iJiUFRUVHZ52rVqgV7e/tKC75t27Z4/PixvKKSTwgMDMSePXuwZcuWcuUukUiwcOFCTp9B8x5doiFEjvLy8j6YQli3bl3Mnj270veMGDECJiYm5UbgDg4ONJVQiUilUrRr1w4CgQBBQUHQ1pbP2JpG8ITIwW+//YY///wTWVlZZR/T19cHj8dDt27dPvreAwcOyDgdkbXDhw/j4cOHuHnzptzKHaARPCHVUlJSgri4uHILeN7/961btz64efbemTNn8OrVq3IjcEtLS7rWrQYEAgEcHR3RoUMH/PPPP3I9NxU8IVV048YNdO3aFe//ymhqasLOzq6ssNevXw8zMzPGKYmimTt3Lnx9fREeHo569erJ9dxU8ETtJCQkVDgCj46OxuLFizFq1KgK35eamooLFy6UFbqtrS09G4V8VF5eHuzt7TFjxgwsWrTog89nZ2dj06ZNGDt2LGxtbTk/PxU8USvjx4/H4cOHkZ+fD+Ddopz/3rzs378/vLy8GKckqiQuLg4WFhYVXr67ceMGunTpgvDwcDg5OXF+brrJStTGo0ePsHPnTri4uODEiRPg8XgwMjJiHYuouI9dluHz+ahRowYaNWokk3NTwRO1IJFIMGnSJGhoaGDUqFFo0qQJ60iEgM/nw83NrVqbc1eHbI5KiII5cOAAAgMDIZVK4enpyToOIQDeFTyXe7D+Lyp4ovJycnIwb968sg2S3d3dGSci5N3ipx49eqBHjx4yOwddoiEqb9myZSgsLMSwYcNgYWFBK0CJQtDQ0MDmzZtlew6aRUNUWVhYGJo3b46VK1di7ty5rOMQFRcVFQVLS0uYmJiwjgKALtEQFZefn48ePXpg2rRprKMQFScWizFo0CAMGzaMdZQyNIInhBAO7Ny5E+PHj8fjx48VZi0FjeAJIeQLZWVlYcGCBRg5cmSVyr24uBhHjhxBZmamTHNRwRO18PLlS8THx7OOQVTUokWLIBKJsHr16iq9/vnz5xg2bBiio6Nlmotm0RC1MH36dBgaGsLPz491FKJiQkJCsHPnTqxbtw6WlpZVeg+fz4e2tjaaNm0q02w0gicqTyqVynxBCVFPUqkUU6ZMgaOjIyZNmlTl9/H5fDRp0gT6+voyTEcjeKIG4uPjkZmZSQVPOPf69WsEBwfj5MmT0NXVrfL75DXgoBE8URkikQhjxozBq1evyn2cz+cDABU84VyjRo0QGxuLrl27Vvk9UqkUNjY2+Oqrr2SY7B0qeKIyduzYgb1795bbmBoA0tPT0aBBA1hbWzNKRlSZqalptV6voaEBf39/jBgxQiZ5yp2L5sETVZCeng5HR0f88MMP2Llz5wefl0qltD0eUTtU8EQljB07FidOnEBUVBTMzc1ZxyFEIdBNVqL0AgMDsWfPHmzdupXKnZD/oBE8UWoSiQTt27eHQCBAUFAQtLVpzEIUl1gsRnBwMJo2bQo9PT2Zn49ushKlduTIETx8+BC+vr4Vlnt2djbEYjGDZEQV3bhx44seLxAZGQlPT088ePCAw1SVo4InSq1hw4ZYuHAhOnXqVOHnR48ejV69esk3FFFJycnJ6NevX5UfR1CR91N25bXpDP08S5RamzZt0KZNm0o/z+fzMWDAADkmIqpq/vz50NXVxW+//fbZx+Dz+ahfv361p1Z+Lip4orKysrIQExNDC5zIF3v06BEOHjyInTt3onbt2p99HHk/MoMu0RCVFRQUBIBWsJIvI5FIMHnyZLi7u2P06NFfdKysrCy5/nmkETxRWaGhoahRowYaNWrEOgpRYgcOHMCzZ89w7949aGlpfdGxQkJCIJFIOEr2aTRNkqgsqVSKjIwMWFhYsI5ClFROTg4cHR3x7bff4vDhw6zjVBtdoiEqS0NDg8qdfJGwsDAYGRlhzZo1rKN8Fip4ojQSExPx9ddf4/Xr16yjEDXRtm1bREVFoW7duqyjfBYqeKI05syZg7CwMBqVE7niYnW0VCqFQCDgIE31UMETpXDv3j0cPXoUPj4+qFmz5idfT7eWiCJJTEyEiYkJrl27Jtfz0k1WovDEYjE8PT2ho6ODR48eQVPz0+OS0aNHo6SkBIcOHZJDQkI+zt/fH3369EF8fDxsbW3ldl6aJkkU3u7duxEcHFzlcgeAx48fo127djJORkjV8Pl8WFhYyP1aPl2iIQotKysLCxYswIgRI9CqVasqvaeoqAivXr2iBU5EYbxfwSrvTWeo4IlCW7x4MUQiEXx8fKr8ntDQUIjFYip4Um3btm1DcnIy58cNDg5m8ueRLtEQhTZq1Ch07NgRlpaWVX4Pn8+HtrY2mjZtKsNkRNXcvHkTkyZNgqmpKYYOHcrpscPCwiAUCjk9ZlXQTVaicnJychAWFoa2bduyjkKUhEgkgpubG0xNTXH37l2V2b+XRvBE5ZiamlK5k2rZsWMHwsLCEBgYqDLlDtAInhCi5tLT0+Ho6IgffvgBO3fuZB2HU1TwhBC1NnbsWJw4cQJRUVEqt2k7zaIhhKitwMBA7NmzB8uXL5dJuWdnZ8POzg4BAQGcH7sqqOCJQvnSDbJXrlyJOXPmcJSGqLqsrCx069YN48aNk8nxg4ODkZCQUK1ZYFyigicK48qVK3B1dUV6evpnH+PChQtITEzkMBVRZV27dsWlS5c4eaBYRfh8PgwNDeHo6CiT438KFTxRCEKhEFOnToWFhcVn/6gsFosREhJCC5yIwuDz+XBzc/vinaA+F02TJArB19cXUVFROHHixGdPU4uMjERhYSEVPFEYfD4fXbp0YXZ+KnjCXEpKCpYtW4YJEyagWbNmn30cPp8PAHB3d+cqGiFf5NixYzA0NGR2fpomSZgbOXIkzp07h8jISNSuXfuzj/P27VsEBARg1KhRHKYjRHlRwROmHj16hDZt2uDPP//E2LFjWcchRKVQwROm2rVrh6KiIjx9+pTZjSiiHoqLizFz5kzMnTsX9erVYx1HLmgWDWFqz5492LdvH5U7kblNmzZh165dKCgoYB1FbqjgCVMuLi5wc3NjHYOouISEBKxYsQJTpkyBi4uLzM9XXFyMrl274tGjRzI/18dQwROVcOjQIWzdupV1DKKg5syZgxo1amDx4sVyOV9oaCiuX7/O/CdTKniiEg4fPowbN26wjkEU0N27d3Hs2DH4+PigZs2acjknn8+HlpbWF0375QIVPFF6UqkUfD6f5r+TD4jFYkyePBleXl4YPny43M7L5/PRuHFj6Ovry+2cFaGFTkTpJSQkICMjg1awkg/s3r0bISEhePz4MTQ15Teefb/JNmtU8ERuiouLoaGhAT09PU6PGxQUBAAK8ReKKBZ7e3ssXLgQXl5ecj3vvHnzmD1B8r9oHjyRm2XLluGff/5BcHAwdHR0ODtuUFAQjh49irVr16rUdmuEfCkqeCIXsbGxcHZ2xtSpU+Hj48M6DiFqgQqeyMXAgQNx//59REREwNjYmHUcQtQCXYMnMnfjxg2cPHkShw8fpnInRI5oBE9kqrS0FO7u7qhZsybu3btH18iJShOLxZgxYwZGjx7NfA48QPPgiYzt2LEDYWFh8PX1lUm537x5E35+fpwfl5DPERUVha1bt37RtpNcooInMiMUCrF69WqMGTNGZlMYt2/fTo8oIACAN2/eoHfv3khKSmKWQdE2naFr8ERmdHV18eDBA5iYmMjsHHw+H/369ZPZ8YnymDFjBkJCQmBqasosA5/Ph4ODA2rVqsUsw39RwROZcnBwkNmxs7OzER0dTQucCC5fvgx/f38cP36c6RZ5ivbIDCp4orSCg4MB0ApWdScUCjFt2jR06tQJAwcOZJqlZ8+eCrWZCBU8UVpaWlro0aMHHB0dWUchDPn6+iIqKgonTpxgPktr5syZTM//v2iaJCFEaSUnJ8PJyQnDhw+Hr68v6zgKhwqeEKK0RowYgfPnzyMqKkphbmwqEip4wpnk5GSYmJjAyMiIdRSiJm7cuIGcnBx8//33rKMoJCp4wgmpVIrvvvsO+fn5uHv3Lus4hMiVVCrFrl270L17d9jb27OOU4YWOhFOXLhwAZcvX8aMGTPkcr6wsDA8f/5cLuci5FOSkpIwbty4spldioJG8OSLFRcXo2nTpqhfvz6uXLkil5kMI0eOxIsXL/D06VOZn4uQTzl37hy8vb0RFxcHOzs71nHK0DRJ8sU2bdqE2NhYnDt3Tm7T1Ph8Plq3bi2XcxHyKXw+H+bm5rC1tWUdpRy6REO+SEJCAlasWIEpU6bAxcVFLucsLi7Gy5cvaYETURjvV7Cynof/v2gET77I3LlzYWxsjMWLF8vtnC9evIBYLKaCVzPvryYrWokCAI/HA4/HYx3jAzSCJ58tJSUF586dg4+PD2rWrCm382ZkZMDe3l4hnrdN5OfMmTPo0KEDBAIB6ygf2LJlC6ZPn846xgfoJiv5IqmpqbCwsICmJo0ViOwUFRXBxcUFTZs2xfnz51nHURp0iYZ8EUtLS9YRiBpYu3YtkpOTcf36ddZRlAoNuwghCi02NhY+Pj6YMWMGGjZsyDqOUqFLNIQQhTZw4EA8ePAAERERqFGjBus4H7hz5w5cXV3leh+qqmgET5RKeno68vPzWccgcnLjxg2cPHkSa9euVchyz8nJQceOHXHu3DnWUSpEBU+UyvLly9GyZUvWMYgcSCQSTJs2DW3btsWQIUNYx6mQom86QzdZSZUFBgaiQYMGzPe8VKQt0YjsaGpqYteuXahRo4ZCzn0H3v15NDAwgJOTE+soFaIRPKmSgoIC9O3bF9OmTWOWQSwWIzg4WGFHS4R7bdq0Uej1Dnw+H25ubtDS0mIdpUJU8KRKVq9ejfT0dCxZsoRZhqioKBQUFFDBE4VRXFyMVq1asY5RKbpEQz7pzZs3WLduHebOnQsHBwdmOcLDw6GhoUGXaIjCOHnyJBR5IiJNkySf1KdPHwQFBSE8PByGhoZMswgEAhgbGzPNQIiyoBE8+ajLly/D398f//zzD/NyB0DlTkg10AieVEooFKJZs2awsbHBzZs3FXYmA1ENxcXF0NHRUdgblsqIbrKSSoWFhSE7Oxtbt26lcicyt3jxYnTq1AkSiYR1lCpJSUlR6OvvABU8+Qg3NzfExcUpxDQ1sVjMOgKRoYiICGzevBndunVTiieTFhcXw87ODnv27GEd5aMU/3eSMKWvr886AoB3ozvaok81SaVSTJs2Dba2tpg1axbrOFXy4sULiEQiNG/enHWUj6KbrEQp8Pl8WFhYsI5BZOD8+fO4fPkyzpw5ozADik/h8/nQ0tJS+IKnETxReFKpFIGBgbTASQUVFxdj2rRp6Nq1K/r06cM6TpXx+Xy4uLjAwMCAdZSPohE8UXhJSUlIT0+ngldBmzZtQlxcHM6fP6+wN/JFIhESEhIQExOD6OhoREdH4/z58+jcuTPraJ9EBU8UHp/PB6C4T+wjn0cgEMDHxwdTpkyBi4sL6zgf8PX1xcaNGxEfH1/uJn/t2rWRlZWFkSNHMkxXNTQPnpQ5fvw4vv76a9SpU4d1lHJKSkrw8uVLuLu7K+woj3yekJAQ8Hg8mW6WIZVKkZaWhpiYmLJR+Ptf9+3bh7p161b4vgsXLuDOnTtwcHAAj8eDg4MDbGxs0KpVK5ibm+P27dsK/+eRCp4AAF6+fAlXV1esXr0as2fPZh2HEE68ePECXl5eKCoqKvuYqakpHBwc4ODggHXr1qF+/fpVPt7mzZsxc+ZMBAYGws3NTQaJuUUFTyCVStG1a1fExsbixYsX0NPTYx2JkHIEAkG5a+D/HYkPHz4c06dPr/B9eXl52Lt3b9kInMfjffZ+BmlpaWjUqBGGDh2K7du3f8FXIz90DZ7gzJkzuHHjBs6fP0/lTpgoLS2Fjo5OpZ/38vJCeHg4AEBPT6+ssNu0afPR6/cmJiaVln91zZ8/H1paWli+fDknx5MHGsGrucLCQjRu3BhNmzbF+fPnWcchKiwlJQWhoaHlRuDv/9vU1LSswCsSEBBQVuyWlpZyX+365MkTtGrVCtu2bcOECRPkeu4vQSN4Nbdu3TokJSXh2rVrrKNUaPv27bh16xZOnDjBOgr5CJFIhMTERFhZWVX6U+D+/fvx22+/QVNTE7a2tnBwcICzszN69OgBR0fHjx6/U6dOMkhddTt37oSrqyt+/fVXpjmqi0bwaiw2NhbOzs6YOnUqfHx8WMep0KBBg5CWloaAgADWUQjeXdP29/f/YAQeHx8PkUiEx48fw8vLq8L3hoaGQldXF/Xr1//o5RhFJBaLkZKSUumMG0VFI3g1Fh0djSZNmmDhwoWso1SKz+fD29ubdQyV99+phKamppVuIp2Xl4effvoJFhYWZdfBW7ZsWXYDs7L3SaVSTJw4Eebm5jh9+rQsvxSZ0NLSUrpyB6jg1VqnTp3w9OlThZ3Lm5OTgzdv3tACJxnYt28f+Hx+2Qg8JiYGhYWFAIBJkybB19e3wvfVrVsXAoEANWrUqNb5jh8/jrt37+Lq1atfnJ1UHV2iIQorICAAX3/9NV6+fInGjRuzjqPw/juVMC8vD8OGDav0td7e3oiOji4bhf93MU/9+vVhYmLCWa6CggI4OTnBy8tLKUfvyowKniis9PR0XL16FYMHD6Zdfipw9epV7N69u2wUnpmZWfY5S0tLpKSkMEz3/xYuXIgNGzYgLCyM6abt6ogu0RCFZWFhgaFDh7KOITdCoRBxcXHlphGOHTsWPB6vwtfn5+cjOzsb7u7u6N+/f7nRuKI8buLNmzdYt24d5s2bR+XOAI3gCWEoPT0d33//PaKjo5GYmFi2Bdz7qYRHjhxB+/btGaf8fH369EFQUBDCw8MVYtP2qiguLsaLFy/g6enJOsoXoxE8IRyRSCRITk7+4IFWrVq1wtixYyt8T82aNWFnZ4evvvqq3Ajc1tYWurq6cv4KuPX69WtcunQJR44cUZpyB4ANGzbg999/R1xcHCwtLVnH+SJU8GpCIpFg6dKlGDNmDOzs7FjHUTk//PAD/Pz8IBQKyz72fiph06ZNK32frq4ujhw5Io+IctewYUNERkbC3t6edZQqi4+Px6pVqzBlyhSlL3eALtGojUOHDuHnn3/GrVu3mK8KrIozZ87g/v37WL9+vdzOKZVKkZOTU+EDrfLy8nDnzp1K33v48GFkZ2eXjcDt7e2rPZWQsDd48GAEBAQgMjKS05lErNAIXg0IBALMmTMHgwYNUopyB97t0xkcHMz5caVSaaXz/nfv3l1uKbqRkVHZ9EE3NzdIJJJKn4HysSmJRDncvn0bx48fx4EDB1Si3AEawauFOXPm4I8//kBERITSXJ5xd3eHp6cndu/eXa33icViREZGfjACfz8qP3r0KLp161bhe9++fYvAwMCya+FmZmYKuwiMcEskEsHDwwOGhoa4f/++3B9mJis0gldxERER2Lx5MxYvXqw05V5SUoIXL15UeGNSKBSisLCw0md6FxUVlS2K0tbWRr169eDg4AB3d3f069fvo5s71K9fv1qbPxDV8eeffyI0NBSPHz9WmXIHaASv0qRSKb777jtEREQgLCwM+vr6rCNVSWBgIDw9PbFgwQJoaWmVG4UnJiZiyJAhOHz4cKXvv3fvHurVq4e6devSAilSJb1794alpSX27NnDOgqnqOBV2Llz5+Dt7Q0/Pz/06dOHdZwq++2337Bv3z6kpqbCxsam3G48Dg4OcHNzg4eHB+uY5H/w+Xy4uLjAwMCAdZRqk0qlKCoqUqrpnFVBBa/Crly5ghMnTmD37t1Kcy05OjoaLi4uGDduHHx8fJTmpw51l5OTg0aNGmHYsGHYtGkT6zjkX1TwRKH0798fT548QXh4OE0zVCLTpk3D3r17ERERARsbG9ZxyL/oJitRGNeuXcOZM2dw9OhRKncl8vLlS/zxxx9YtWoVlbuCoRE8UQilpaVwdXWFmZkZ7ty5ozSXlNSdVCpFly5dEB8fj9DQUNq0XcHQCJ4ohG3btiEiIgKBgYFU7krk9OnTuHnzJi5cuEDlroBoBE+YS01NhaOjI4YOHYrt27ezjkOqqLCwEC4uLmjWrBnOnz/POk6VvXnzBnFxcfj6669ZR5E5GsET5vT19TF+/HjMnj2bdRRSDU+fPkVeXp7SzZqZNm0anj9/jqioKKV/Yuen0AheRWRmZmLKlCnw8fFRmhWrRPnl5+cr1Q3xixcvomfPnjhx4gQGDBjAOo7MUcGriAkTJuDIkSOIjIxUicecEsI1oVCIpk2bws7ODtevX1eLez10iUYFBAcH488//8T69eup3AmpxObNm/H27VucPn1aLcodoBG80pNKpejYsSMyMjIQEhICHR0d1pEIUThJSUlwcnLCqFGjsGXLFtZx5IZG8Eru+PHjuHv3Lq5evUrlTkgl5s2bB319fSxbtox1FLlSnediqqH8/HzMmjUL/fr1Q9euXVnHqTL6oZHIk1AoREpKClavXl3pY6ZVFRW8Elu9ejUyMjKwYcMG1lGqTCqVYuDAgdi3bx/rKKSaTp48CYFAwDpGtenq6uLKlSv45ZdfWEeROyp4JdaxY0ds2rQJDg4OrKNUmb+/P06dOgULCwvWUUg1PHv2DIMGDcLRo0dZR/ksGhoaanNj9b/oJiuRm+LiYjRu3BhOTk64ePGiWv6FU0YSiQTt2rVDQUEB+Hw+tLXp1p2yoO8UkZsNGzYgISEBly5donJXIocPH8ajR48QEBBA5a5kaARP5CI+Ph7Ozs6YOHEi1q5dyzoOqaK8vDw4OTmhY8eO+Pvvv1nHIdVEBU/kYvDgwbh9+zYiIiJgYmLCOg6pojlz5uCPP/5AREQEPQJDCdHPW0Tm7t69i+PHj+PgwYNU7kokIiICmzdvxuLFi5Wq3C9dugSJRIKePXuyjsIcjeCJzBUUFODgwYMYN24cNDVp4pYiKS4uxtu3bxETE4Po6GjExMSU/ffbt29hZWUFPp+vNHvjFhYWwtnZGW5ubvD392cdhzn626YkHj58iG+++QYZGRmso1SbkZERJkyYQOXOQElJCYRCYaWfP3r0KJo0aYKePXti+vTpOHPmDHJyctCiRQuIRCJ4e3srTbkDwJo1a5CamoqNGzeyjqIQaASvBCQSCby8vCCRSPD06VNoaWmxjkQUSGlpKe7fv4/o6OiyUfj7/05KSoKfnx+8vb0rfG9ycjKioqLg4OAAGxubsj9bKSkpsLa2xsmTJ/H999/L88v5bNHR0XBxccHMmTOxcuVK1nEUAl2DVwL79+9HYGAg7t27R+WuZiQSCZKSkiASicDj8Sp9zTfffAOpVAobGxs4ODiAx+OhU6dO4PF4cHNzq/T41tbWsLa2/uDjQUFBAIAWLVpw8WXIxcyZM2Fubo758+ezjqIwqOAVXE5ODubPn49hw4ahXbt2rOMQGbp+/TqePn1abgQeFxcHoVCIAQMG4MSJExW+T09PD1FRUahbty5nl1P4fD5MTU0r/UdF0Vy7dg1nzpzB0aNHlWoDElmjgldwS5cuRVFREdasWcM6CvkMUqkU2dnZiImJQWJiInr37l3pa/fu3YvLly/DwcEBDg4O8Pb2LhuNu7i4fPQ8DRo04DS3hoYGunfvrhQL0kpLSzF16lR06NABgwcPZh1HodA1eAX28uVLuLq6YtWqVZgzZw7rOFUmFApVfq/Lyjx//hz79+8vdz08Ly+v7POFhYUwMDCo8L3q/Pv2JbKysjBmzBgsWrToo5ej1BEVvALr2rUrYmNjERoaCj09PdZxqiQ/Px9ubm5YtWoVBg0axDoOJwoLC8tNH+zSpQucnJwqfO21a9cwZcoU8Hi8stH3f381MzNTilExUQ10iUaBLVmyBBKJRGnKHQBWrlyJxMREeHl5sY7y2UQiEX766aeyEXhqamrZ57S1tXHw4MFKC75r16549eqVvKIS8lFU8Aqsffv2rCOUyc/PR3JyMpKTk5GZmQmRSITS0lJoaGhAW1sburq6KC4uxoYNGzBnzhyFujknEomQkJBQbiGPjY0Nfv311wpfr62tjfz8fDg7O6NHjx7lRuD/nUpIiKKjSzSkjFQqRXx8PAIDA9/979kzvH37BskpKRAI8qt1rFqmprCxsYajkxM8PDzh4eEBDw8PuT4HftGiRTh8+DDi4+MhFovLPm5jY4M+ffpg+/btcsuiTF6/fg0bGxsYGhqyjkK+EBW8mktISMC5c+dw4fx5PH78GBmZmQAAS/NacHOqhwZ2VrC2MIWVuSkszU1hbW6K2jVrQEdHG1r/rkwVicUQCkVIz85DakYOkjNykJqZg6S0bETEJON5ZCxy8goAAHa2ddG2XXv07t0bPXr0QO3atT+ZUSqVIi0trdxy+ujoaMTGxuL8+fOVPsL2r7/+QlhYWLkReL169ZRqZSYLzs7O6Nq1K3x9fVlHIV+ICl4NBQcHw8/PD+f8/cEPCoK2thbauDqhrVsjuDnz4ObMg7VFLc7OJ5VKEZOYjuDwGAS9isbdoAgEhb2FlpYW2rdrB+8+fdCvX78Kd6a6evUq+vXrh8LCwrKP1apVCzweDzweD/v27VO7fTZlSSAQoGbNmti3bx9GjBjBOg75QlTwaqKoqAh///03tm/bhmeBgahpbISubZqhRwc3dG7dDKbGRjI9f15BEeKSMhCbnI645AyEvUlAUFg0YpLSUVxSCpFYjG7ffosJEyeiZ8+eZde54+LicPr06bIRuL29PRW6DN27dw8dOnRASEgImjdvzjpOhdasWYOGDRsqzSMUWKKCVxBisRgaGhqcP5DrzZs32LFjB/bt24ucnFx0ad0Mo/p/jS5tmkGHw915xGIJtLQqzi6VSmH7zTgUFJUAAPT1dFDPyhz1rM1hb2OB3l97IiktC3tP30LgyzeoZ2eHX8eNw+jRo1GnTh3OMpJP27p1K+bMmQOBQAAdHR3WcT4QHh6OZs2aYenSpViwYAHrOAqPCl5BbN26FX///TcCAgI4WeySkJCAJUuW4MCBA6hpbISferXHyH5fw8H28wszMS0LkdFJiE3OQFxyBmKT0st+bdLQDme2zq70vZfvBaN2zRqwtzGHRS2TSv8hC3oVjb2nb+Lk1cfQ1NLCjBkzMWvWLNSsWfOzc5OqmzVrFu7evYvHjx+zjvIBqVSK7t274/Xr13j58iXdS6kCKngFkJaWBkdHRwwePBg7d+78omNlZWVh9erV8PX1RQ1DPcwa0Qsj+nSCgf7H/9EQicRITMtCPWvzShfiTF9zAPvPBEBTUwM2dWrD3vr/R+FNG9mhV0ePL8pe7uvIzcfWwxex85/rMDKqgd8WLMCECRPoL7UcKOqK2rNnz6Jv3744e/ZspU/HJOVRwSuAMWPG4NSpU4iMjIS5uflnHUMsFmPr1q1YtmwpREIhJg7phklDusPEqPyy+PSsPNx4HPrvyPv/R+KJaVkQiyWIuLAZlmamFZ4jPiUTYrEYdS1rc3p552OS0rKxdt9ZHDp3BzbWNti4aRMGDBggl3MTxVFcXIzGjRvD0dGRNm2vBip4xp49ewYvLy/4+vpi4sSJn3WMyMhIjBg+HA8fPULvTh6YNaI3XJ15Fb728fModBu7Eua1jN+Nvq0tYG/z/yPxtm5Onxzts/A6LgWL//gHF+/wMWjQQPzxxza5zqknbK1YsQLLli3DixcvKl1FTD5EBc+QRCJBu3btUFBQAD6fX+l87ops2bIFERERCAgIQHh4OAAp3n8nl00chKk/fVfh+0pFIpQIRahhqHyXOqRSKU5ff4zZG45AS1sP23fsoNG8GoiPj4eTkxMmTpyIdevWsY6jVKjgGfrrr78wfPhwBAQEoGPHjhAIBGVPIdTQ0PjodcYmTZogJiYGhYWFaOZYD993bYWG9axgb2MBh7p1lLLAqyotMxcz1v2F8wGB+PHHwdi9ew+MjGQ7zVMdlJSUKORzj1JTU/H7779j9erVtGl7NVHBMyIUCmFlZQUjIyNYWloiOjoaWVlZZZ9v0aIFAgMDK3zv48eP0a9vX2hKS7F72Vi0dVO/H1mlUilOXXuMaT4H0aBhI5z194e9vT3rWErN09MTHTt2xIYNG1hHIRyhXZAZuXLlCjQ1NWFrawt3d3fMmDEDR44cwYMHD5CUlISnT59W+L5Dhw6hY8eOsLMwxs19i9Sy3IF3G1IM+LY1ruz6DTmZKWjp6Ym7d++yjqW0hEIhQkNDK1xNTJQXPU2SkfebZz948KBKMwIkEgnmzZuHdevWYWivDtg452fo6SreQhR5a9LQDjf3LsLwBdvRuXNnbNu2DWPGjGEdS+mEhYVBKBTC3d2ddRTCISp4Rvh8Plq0aFGlcheJRBg1aiSOHDmK1dOHYNygrjRN7D/MTI1xZstMzNt4FGPHjkVmZibmzZvHOpZS4fP50NDQgKurK+sohENU8IwEBQVV6WFOIpEIw4YNxcmTJ7Fn2a/o37WV7MMpIR1tbayf/RPMaxlj/vz5KC4uxtKlS1nHUhpBQUFwcnKiDatVDBU8AykpKUhKSkKLFi0++jqxWIyff/4Jp06dwoGVE9G7E3crRVWRhoYG5o/pB309HSxbtgza2tpYuHAh61hKYe3atUhOTmYdg3CMCp6BWrVq4datWx99Wp9UKsXo0aPxzz8nsG/5OCr3apj+cy+IRBIsWrQIhoaGmDFjButICs/AwAD169dnHQNSqRQDBgzATz/9hL59+7KOo/RoFg0Denp66NSp00c3u/Dx8cGBAwewY9Ev6PNNSzmmUw2zR3lj+s89MXPmTJw9e5Z1HFJFx44dw+nTp2FsbMw6ikqgefAKyN/fH3379sXsUd74bUw/1nGUlkQiwc+/bUPAs3A8fPgITZs2ZR2JfER+fj6cnJzQpk0bnDx5knUclUAFr2BevnyJ1q1boZOnM/5aNZHz58Orm/zCYnQbuwpFYi08efoUZmZmrCORSsyfPx+bN2/Gq1evFGrTdmVG7aFAsrKy4N27N+ytzLBz8Rgqdw7UMNTH0bWTkZuTiYEDBkAkErGOpHAGDRqElStXMs0QFRWFjRs3Yt68eVTuHKIGkbOHDx/i+++/R05OzgefmzBhPLKzMnBk7WSVfpaMvNnbWODgygm4fecO1q5dyzqOQpFIJGWrqlmaPn06rK2tMWfOHKY5VA0VvJzdu3cPV65c+eChSadOncLx4/9g3cyh4NnQY3C51r6FM6YO64GlS5fixYsXrOMojOjoaOTl5X1yyq4sXbx4ERcuXMCGDRtgYGDw6TeQKqOClzM+nw83N7dyI6aMjAyMHzcOPTt64HtayCQz80b3RQM7S4wYPhylpaWs4ygEPp8PAEwfUeDl5YWNGzeif//+zDKoKip4OQsKCvpgtDRp0kSISkuwcc5P9AgCGdLT1cG2BaMQFBxMzxX/F5/PR926dZlubm5ubo7p06fTn30ZoIKXI4FAgMjIyHIFf/78eRw//g/WzhhS6VZ5hDseTeqXXaqJiopiHYe5MWPG4K+//mIdg8gITZOUI4FAgKNHj6J79+6wt7eHWCyGa/PmqG0A+P8xm0YwclJcIoTHoPlo36kz/v77OOs4hMgMFTxDBw8exIgRI3Bj72J4NGG/TFydHDp3B5NX7sOzZ8/g4UGPgSCqiQqekZKSEjg2agTXBpY45DOJdRy1IxKJ0XbYItg3bIyr166xjkOITNA1eEZ27NiBhMRELBpHMwdY0NbWwqJx3+Pa9eu4ceMG6zhqRSwWs46gNqjgGSgtLcUaHx8M6dkOjjwb1nHUVq+OLeDRpAFWLF/OOgoTCxYswK5du+R6zuzsbDg7O+Pq1atyPa+6ooKXk/DwcCxYsAA5OTnw8/NDSmoqxv/wLetYak1DQwPjf+iKgNu38erVK9Zx5O7gwYOIiYmR6zmXLFmClJQUNGvWTK7nVVdU8HISEBCAtWvXQl9fH9u3bUMbNyc0aWjHOpba693JAxa1a2LHjh2so8hVamoqEhMT5bqCNTQ0FNu3b8fixYthbW0tt/OqMyp4OeHz+WjSpAnevHmDgNu3Meb7b1hHIni3+Onn3h1w8OAB5Ofns44jN0FBQQDkt4JVKpViypQpaNCgAaZOnSqXcxIqeLl5v8n2jh07UMfMFL1ohyaFMbxvJ+TnF+Do0aOso8gNn89HzZo15baL08mTJxEQEIAtW7ZAV1dXLuckVPByIRQKERoaCldXVxw7ehRDe7aDrg7tlqgo6lmbo0vrZjikRis6v/76a6xbt04ui+sKCwsxc+ZM9O7dG927d5f5+cj/o4KXg8LCQkyYMAEmJibIys6m/VUVUK9OLfDg4UNkZGSwjiIXbdq0wZgxY+Ryrs2bNyM1NRWbNm2Sy/nI/6OClwNTU1Ns2rQJL1++hKV5Lbg581hHIv+jezs3SKVSXLx4kXUUlTNp0iT4+fmhQYMGrKOoHSp4OZFKpTjr54fu7Zoz31yBfKiOWU14Nm0Af3/aoJtrJiYm6NGjB+sYaomaRk4iIiLw+s0b9OjA7rnb5ON6tHfD5cuXUVxczDoKIZyggpeT69evQ1dHBx09G7OOQirRtW1zFBQU4tGjR6yjyNTevXvh7+/POgaRAyp4GUtKSsJff/2Fhw8fonFDOxjo0xQxReXiUBcG+noIDAxkHUWm1q9fj2v0gDW1QAUvY7du3cLw4cMRxA+Eu3M91nHIR2hra6G5o71KF3x+fj4iIiKY7sFK5IcKXsb4fD4cHBwQERkFNyce6zjkE1yd6iHw2VPWMWQmJCQEUqlUpitYc3NzZXZsUj1U8DLG5/PB4/EgkUjgStMjFZ6bMw+RUa+Rl5fHOopM8Pl86OrqonFj2dwLevLkCWxtbcs28yZsUcHLkEQiAZ/Ph4mJCbS1tdC4gS3rSOQTXJ3sAbx7MJYqcnBwwJQpU2TyuACJRILJkyejQYMGcHV15fz4pPpovbwMCQQCtG3bFsbGxrAyr0WPJ1ACtlZmAICEhATGSWSjV69e6NWrl0yO/ddff+HJkye4ffs2tLS0ZHIOUj00gpehmjVr4tKlS9DV1YWVuSnrOKQKTIwMYKCvh+TkZNZRlEpubi7mzZuHH3/8EV999RXrOORfVPBykJSYCEszE9YxSBVoaGjAyqIWkpKSWEdRKsuXL4dAIMDatWtZRyH/QQUvB8nJSbCmEbzSsDKrSSP4anj16hW2bNmCBQsWwNaW7jMpEip4OUhOToalmSnrGKSKLM1MkJSUyDoG565fv45nz55xflwfHx/Y29tjxowZnB+bfBkqeBnJzc3FkydPUFpaivyCAtQw0mcdiVRRDUN9FBYUsI7Bufnz5+OPP/7g/Lg7duzAuXPnoK9Pf8YVDRW8jNy+fRutWrVCamoqRCIxtLXot1pZ6GhrQSQSsY7BqdLSUjx//lwmK1gNDQ3h4uLC+XHJl6PWkRE+nw8LCwvUrVsXIpEIWvSIYKWhpampcgX/6tUrCIVCekSBmqHWkZH3e7BqaGhAS0sLYomEdSRSRWKJROXmcfP5fGhoaNACJDVDBS8j7wseAHR0tCESU8Eri1KRGNraqrUoTVNTE9999x2MjY1ZRyFyRAUvAwKBALq6umUFb2hgiMKiEsapSFUVFpXAwMCQdQxO/fzzzzh//jzrGETOVGuYoiCMjY3x9u1bSKVSAICVlSVSMnLYhiJVlpKZB3tn2j+0IiKRCHFxcahfvz7rKKQKaAQvQxoaGgAAG5u6SM2kR6gqi9TMXNjY2LCOoZD+/PNPNGnSBCkpKayjkCqggpcDaxsbpFDBK43UjBxYW1uzjqFwMjIysGjRIgwbNgxWVlas45AqoIKXAxsbG6RmUMErA0FBEfILi1RqBB8REYHU1NQvPs7ChQshkUiwcuVKDlIReaBr8BwTCoUoKipCzZo1yz5Wt25dJKVnvVvwpK1a0+9UTWJaFgCoVMGPHz8etWrVwqlTpz77GEFBQdi1axc2b96MOnXqcJiOyBKN4Dl2//59mJqaIioqquxjzZo1Q2mpCBEx9IRCRRcaGQfg3fdMFUil0nJTdj/3GJMnT0bjxo0xfvx4DtMRWaMRPMf4fD4MDQ3LzTJwd3eHhoYGgsJj0KShHcN05FOCXsXAgWeP2rVrs47CiejoaOTm5n5RwR89ehT379/H9evXoaOjw2E6Ims0gucYn8+Hm5tbuZWQxsbGcHJshJDwGHbBSJUER8TCs6UX6xiceb836pcU/N69e/H999+jc+fOXMUickIjeI7x+Xx06dLlg497eLZE8PPHDBKRqpJIJAiNjIX3DyNYR+FMSkoK6tevD0tLy88+xuXLlyEQCDhMReSFRvAcKioqwps3byocLXl4eCA0Mg7CUtV6iJUqiYpNgaCgCB4eHqyjcGbSpEnl7gd9Dl1dXZiZmXGUiMgTjeA5ZGBggJycnLIVrP/1zTffoLhEiPv8cHzdqimDdORTrj18Dn19PbRp04Z1FE5p0pNM1RZ95zlmaGgIIyOjDz7evHlz1LOzxcW7QQxSkaq4dC8EXTp3qfD7R4gyooKXEw0NDXj36YvL90MqHOETtrJy8/EoJALeffqwjkIIZ6jg5cjb2xvxyRl48TqedRTyP649eA6xWIJevXqxjsIZgUBAgwk1RwXPkar8RerYsSOMjWvgwm2+HBKR6jh/mw+vlp4q9QyaYcOGoV+/ftV6T3JyMl68eCGjRETeqOA5EhgYCBsbG0RGRlb6Gl1dXQwYMBCHz9+DmDYAURipmTm4dC8IPw4ZyjoKp/h8Ppydnav1njlz5uDbb7+FUCiUUSoiT2o7i0YqlSI+Ph5BQUGIj49HUlISkpOTkZyUhOTkJOTl5UEkEkEkEkFTUxPa2trQ19eHlZU1rG1sYG1tDWtra9jY2KBJkyZ48uQJ0tLSYGf38ZWqEydOxP79+3Hlfgi++8pdTl8t+ZhD/nego6ODESNGsI7CmbS0NCQkJFRrgdP9+/dx+PBh7N69G7q6ujJMR+RFQ6omF+lyc3Nx69YtPHv2DM+ePkVgYCAyMjMBvNtSz8q8FqzMTWFZ2wTWFqYwMTKAtrYWtLQ0IZFIIRKLUVxSirTMXCRn5CI1Kxcp6dnIyy8E8G4qmp6eHoYPHw4PDw907NgRjRo1qjBLKy8vGGsV49TmmXL7+knFRCIx3AbMRbeefbBnzx7WcThz5coVdO/eHVFRUWjYsOEnXy8Wi9GyZUtoaWnh8ePHNLVSRaj0CD4mJgbnzp2D/9mzuH3nDkpLS1HHzBTuzvYY1bc93Jx4cHW2h41FrbLNOapLUFCEl68TMHrxDmhrayHg6nns2rULEokEzk6O8O7TF97e3mjdunXZ4wsmTJyIESNG4G18Kurbff4KQ/LlrjwIQUJKBiZMmMA6CqdCQkJgYmJS5Z2X9u7di6CgIDx48IDKXYWo3Ag+OzsbBw8exP59+/A8NBQ6Otro0MIF3du74tt2rrC3Nv/sMq9MqUgE22/GYdnEHzDuh67ILyzG7WdhuHQ3CJfvPUdGdi7MzcwwYOBATJgwAQ0bNoStbV0M6OyBtTOHcZqFVE+viWtQqmWMR49V6zESEokEycnJqFu37idfm52djUaNGqFnz544ePCgHNIReVGZgufz+di2bRuOHTsKkUiEXl+1QJ9vPPFN62YwMTKQ6bklEgnC3iTCorYxLM1My31OLJYgMOwtLt7h4+/LD5GSno327drB2sYGZ/388Oyf1ahnbS7TfKRitx6/QL+p6+Hn54c+ajz/fcqUKThw4AAiIiJUahYRUfKCl0qluHLlCpYtXYpHjx/D1socI/p8hZ+8v/qgaBVBqUiEi3eCsOfULdwNDIOOtjYGfNsKOxaPYR1N7UgkEnzzy3IY1LTE/QcPOP+pTlmUlpaiTZs2GDx4MGbNmsU6DuGY0hb8o0ePMG/uXNy+cwetXR0xeWh3dGvrqjQ7JkVEJ2H6mgN4GByJ+0dWoHEDW9aR1MqZ608wcuF23L59G1999RXrOEyJxWJIJBJ61rsKUrqCDw8Px/x58+B39iyaNKyHxeP749u2rko5AhOWitBi4Fw0bWiHv9dPYx1HbZSKRGg9ZBGcmrjhwsWLrOMQIjNKc7u8tLQUK1asQPPmzcF/+hB/LhmLOweXols7N6UsdwDQ1dHG0gkDcfleMK7cD2YdR21sO3oF0QmpWO3jwzqKTIwZMwZjxtBlP6IkBR8aGorWrVph6dKlmDKkG57+vRI/9GgLLS328d/EpaD3xDWITkj7rPd/37UVurRuhmk+B5AjKOA4Hflf4dGJWLXbDzNnzkTz5s1Zx5GJ+/fv00IlAkDBC14kEmHFihXw8PBAUV46ru9ZiEXjB0BPV3GuFQaFx+Bu4CvUMvm8R8xqaGhgy28jUVBYgt82H+M4HfkvkUiMiSv2oX79+vj9999Zx5GJgoIChIeHf9EWfUR1KOxCp+zsbPwwaBBu3rqFacN6YM4vfRSq2N8LiYhBPWtzmH5mwQNA3Tq1sXLaj5i8ch/6dm6Jb9u6cpiQvPfH0csIehWN+/fvQ19fn3UcmQgJefc4aip4AijoCP7Vq1fwatkSz54+xuktMxVu1P5fIRGxcHWy/+LjDOvVAV1aN8PU1fuRlpnLQTLyXyERsWWXZlq3bs06jszw+Xzo6OigSZMmFX7+4cOHePr0qZxTEVYUruAvXLiAVq28oIMS3Nq3CB09G7OOVCmpVMpZwWtoaMB3wShIJFL8NM8XJcJSDhISAEjLzMXQub5o1qyZyl6aee+HH37AtWvXKrwGX1pail9++QVz5sxhkIywoFAFv337dvTu3Rvt3Rrh6q4F4NWtwzrSR4nFEmyeOwK9O3lycjxri1o4vGYygsJjMHPdX7RZAwdKhKX4af4fEEm14Hf2rMpemnnPwsICHTt2rPBzf/zxByIiIrB582b5hiLMKMw8+I0bN2LmzJkYP/hbrJwyWK0feHT0wj1MWL4HPtOHYtwPXVnHUVpSqRRTVu3H8SuPcPv2bZW+NPMpqampcHR0xLBhw7Bt2zbWcYicKMRN1vXr12P27NmYObwXFo77XmnntXNlSM/2CHuTgAVbj6GBnSW6tlXN6Xyy9sfRyzh07g4OHDig1uUOAL/99hu0tbVV/hIVKY/5MNnX1xezZ8/GrJG9sWj8ALUv9/eWTRyEb9u64qf5vrj9LIx1HKWz7/RNLPI9jnnz5mH48OGs4zD15MkT7Nu3DytXroSZmRnrOESOmF6iOX78OAYPHoxJQ7pj+eQfqNz/R4mwFEPmbMHD4Ej8s3EG2reo3vZr6uov/9uYsmo/Jk+ejC1btqjNn6vVq1ejqKio3ChdIpGgTZs2KCkpQWBgYNmeBEQ9MBvB8/l8jBw5AgO/baOU5Z6amYOZa/9CXHKGzM6hp6uDQz6T4dWsEQZO34ibj2kz5E/ZffI6pqzaj19//VWtyh0A/Pz8EBMTU+5jWVlZ0NbWhq+vL5W7GmJS8Kmpqejj7Q1nng22/jZSKf8S8sOisff0TWhqyja7ob4e/l4/FR08nDF41macuPJQpudTVhKJBD57/DB7/WFMmzYNO3bsUMo/V5+rtLQUISEhHyxwMjc3x71799ChQwdGyQhLci/4kpIS9O/XD6UlhTiyZjIM9JXzmRkhEbEwMzVG3Tq1ZX4ufT1dHPKZjH5dvDBmyZ9Yuu0fiMUSmZ9XWRQUlWDEgh3w2eOH5cuXY+PGjWpV7sC7p6yWlJRUuIJV3X4vyP+T+yyaGTNm4NmzZ7iwYy5s6tSS9+k5ExIRi+aO9eT2l0dPVwc7F49B04Z2WLLtH7x6m4jdv4+T+W5Vii4uOQND5/oiOikDZ86cQd++fVlHYoLP5wMA3Nzc2AYhCkWuI/ibN29i+/btWD55EFo2/fRO74rsOUcrWKtDQ0MDk4f2wD8bpuNhSCS6/rIcr94myjWDIrn9NAzfjFqOfKEGHjx4qLblDgCtW7fGzp07YWJiwjoKUSBym0UjEAjQrGlT2FkYwd93tlIvZBKLJdhw8Bw6tHBBGzdHJhmiYpMxbK4vohPTMG90X0wZ2kNpdrP6UoKCIiz54x/sO3MLXTp3xrG//4a5Oe1rS8j/klvBjx8/Hn8d3I8Hh5cr/CMIlEVxiRCrd/vB9+gluDnzsG3haLjUr8s6lkzdfhaGySv3IzOvAGvXrsO4ceOUerBAiCzJpeBv3ryJzp07Y+3MYRg7sIusT6d2nr14gwnL9yAmKR0zR/TGpCHdYWSgxzoWp9Iyc7Fq9xkc8AtAp44dsW//fjg4OLCOpXD27duHBg0aVPo8GqJeZF7wYrEY7m5uMNIS4uKOeTTakpHiEiHW7D2LbceuwNTYEHNG9cHwvh2ho60QT6P4bHkFRfA9cgnbj12Fjp4eVq5cRaP2SsTFxcHZ2RmTJk3C2rVrWcchCkDmBX/48GH89NNPuLZnodLfWFUGsUnpWL3HD8cvPYCdlRmWTBiAfp29lK4Qi0uE2Hv6FjYcPI/CYiGmTJmKefPmoVYt5Z15JStHjhxBQUEBbty4gTt37iAiIoJuthIAMi74kpISODs5oQnPHEfWTJbVaUgFXr6Ox6QVexEUHgNe3Tr4pf/XGNqrA2rXrME62kfFJqXjgF8ADp27h+y8fIwcORJLliyBra0t62gKq0uXLhAKhbh79y4OHjyIn3/+mXUkoiBkWvC+vr6YNm0aHhxZDmcH1b75p4ikUinaDFmI9LwiCAT50IAU/bt44Zf+38CjSX2FWQAjFktw8/EL7Dl9E1fvh8DExBjDh4/AxIkT4ejIZpaSspBKpTAzM4Ouri7q16+Pe/fuKd1Pa0R2ZFbwhYWF4Nnb49vWLti28BdZnIJUwZ3AV/CeuAYHDhxASkoKdu7YjpjYOPDq1kH39q7o0cEdbd0c5X6tvqCoBLeevMClu8G4cv85MrJz4ebqiomTJuHHH3+EkdHn73GrTmJiYspuNj99+hSentxsPkNUg8z+Vh8/fhwZmZmYNaK3rE5BquArDxe0dnXEwQMHcPPWLcyaNQs3btyAn58f/M+exc7j11DT2Ahd2zRDG1dHuDnz0KShLfT1uH2ERF5BEZ5HxCI4PAZ3A18h4GkYSoSlcHF2wi9jx6F///5o2bKlwvxUoQykUikuXLgAABgyZAiVO/mAzEbwLT09UUu3FCc2zZDF4Uk1nLz6CKMX70RYWBhcXFzKPi6VShEUFAR/f39cOH8OwSHPIRKJoK2tBZf6tnB1qodG9axhaV4TVmamsLIwhZWZKWoaG35QxBKJBJk5+UjJzEFqRg5SMnORkp6N8OgkBEfE4XVsEgDAwEAfXi294N2nD3r37o1GjRrJ9fdCGQkEAty6dQsxMTGIjo5GdHR02X/n5eXB1tYWgYGBqFOH1peQ8mRS8E+fPoWXlxf+Xj8N3du7cX14Uk0lwlI07TsLg4f+jK1bt1b6uuLiYoSGhiIwMBCBgYF49vQJomNikJubV+51mpqa0NHWgva/j58tFYlRKhJ9sIesuZkZGjVqBM+WLeHh4QEPDw84OztDW8mnbnKpqKgIMTExMDIyQr169Sp8TVRUFBwdHaGnpwcejwcHB4dyv3bt2pVmF5EKyaTgR44ciZtXLiDo5BpoadENH0WwfMdJ7D5zG4mJSahRo3ozaQoLC5GcnFz2v8zMTIhEIpSWlkJTUxPa2trQ0dGBhYUFbGxsYG1tDSsrK+jqKueTQmXl2LFjCA0NLTcCT01NBQDMmjUL69atq/B9YrEYaWlpsLS0pBuopFo4L/jc3FxYWVlizsjemDG8F5eHJl8gPiUTrv1nY9euXfjlF7rpzRWRSITExERER0dDIBCgd+/K7zm1b98e8fHxH4zCHRwc4OTkBEtLSzkmJ+qA85+Vr1y5guLiEgzs1obrQ5MvYGdlhjZuTvA7c4YK/gvcuHEDR48eLRuFx8fHQyQSAQAsLS2RkpJS6Xvv3LlDI3AiV5wXvL+/P5o52sPOijb3VTTd2zXHyt1+KCgooGmI/5JKpUhPTy932WT48OGwtrau8PWJiYkIDQ2Fg4MDWrZsWTYC5/F4sLf/+OOjqdyJvHF6iaa0tBR16lhgTP+OWDC2P1eHJRx5E5cCj0Hz4Ofnhz59+rCOw0xKSgpGjx5dVuqFhYVlnzM1NcWFCxfQtm1bhgkJ4QanI/j79+8jJycX33Vw5/KwhCMN6lnBkVcX/v7+KlXwAoGgbPT9/lcvLy/8+OOPFb6+Ro0a0NDQQOfOnT+4Hm5qairf8ITIEKcFf+HCBVhZ1JL7Tkek6nq0d8Xf589DKpUq/aKiYcOG4fLly8jMzCz72PuphJVdYgHeFfy5c+fkEZEQpjgt+IcPHqCtayO61qjA2ro7Ysvhi+WWuLMmFAoRFxf3wUKeoqIinDlzptL3tWnTBi4uLuVG4DSVkJD/x1nBi8ViBAUHYf4v3lwdksiAmzMPABAYGKgQBb9r1y6MHz8eEokEwLsbkba2tuDxeGjYsOFHf9KYOHGiPKMSonQ4K/iIiAgUFhbB9d8CIYrJ0swUNnXMEBgYiAEDBnByTIlEguTk5AqX0kdHR+PQoUNo3759he9t164ddu3aVTYCt7W1pQVShHCEs4IPDAwEALr+rgRcneoh8NmzKr9eKpWipKQE+vr6FX4+Kyur3PPazc3Ny6YPtmzZEhYWFpUeu0mTJmjSpEnVwxNCqozTgq9vZ42aNQy5OiSRETdnHv48GfDB5Y9Xr17h1atXH4zAY2JiMGTIEOzevbvC45mZmeHcuXPg8Xjg8XjVfhQCIUQ2OCv4N29eo1E9WmqtqAQFRUjJyEEje2s48ayRlZ2N3NzcctMCZ8yYgcuXL8PIyKjsxuX7qYReXl6VHltDQwO9etFjKQhRNJwVfFJiIprbm3J1OPKZohPScONRKOKSMxCbnI7YpAzEJWcgKzcf+no6SA7YBUszUwBAUlJSuYL/888/YWhoCDMzM6WfQkkI4bDgk5OT8a0Hj6vDkc8UGPYW8zYdha1VbdhbW6BZo3ro1bEF7G0sUM/aHFKpFNYWpgDefc8aN25c9t7KHldLCFFOnBS8WCxGalo6rMxNuTgc+QLeX3ui7zctoa2tVelr/juCJ4SoLk4KPi0tDRKJhApeAejqfPpbaqCvC1OTGkhOTpZDIkIIK5ws+cvIyAAAmJnS7AllYWZqXPZ9I4SoJk4K/v3zsHVoKzaloaOtVfZ9I4SoJk4KvrS0FADK9ugkik9bS6vs+0YIUU2cFDxNqVNO9H0jRLVxUvDa/16aEYnFXByOyIFILC77vhFCVBOnBV9K13SVhrCUCp4QVcdJwZubmwMAMrIFXByOyEFmTl7Z940Qopo4Kfg6depAU1MTqZm5XByOyFhhcQlyBQUf3fWIEKL8OCl4LS0tWNaxQHJ6NheHIzKWmvHuH2IbGxvGSQghssTZ3mbW1tY0glcSKRk5AEAjeEJUHGcFb1O3LpLTc7g6HJGh9wVPI3hCVBtnBd+gQUNExaVydTgiQ5GxyahdqxZq1qzJOgohRIY4K3gPDw+8jU9GjqCAq0MSGQkKj4aHhwctdCJExXFW8J6engCAkIhYrg5JZCQkIg6eLVuyjkEIkTHOCt7R0RFGRoYICY/h6pBEBlIzc5CclgUPDw/WUQghMsZZwWtpacHdzR3BVPAK7f33hwqeENXHWcEDQOs2bXA/OAoSiYTLwxIO3edHwLJOHdjb27OOQgiRMU4LvlevXkjNyKZRvAK7dD8EvXr3phushKgBTgu+Xbt2qFXLFBfvBnF5WMKRqNhkRMUkwdvbm3UUQogccFrw2tra6NmzFy7fC+HysIQjl+8FQ19fD126dGEdhRAiB5wWPAB4e3vjRVQc4pJpv09Fc+leCLp26QpDQ0PWUQghcsB5wXfr1g36+vo4ceUh14cmXyAuOQOPQiLRt18/1lEIIXLCecGbmJhg8ODBOOB3G2IxzaZRFAf8AlCjhhF++OEH1lEIIXLCecEDwMSJExGfkoGrD+havCIoEZbiL/+7GDFiJIyMjFjHIYTIiUwK3tPTEy09PbH39C1ZHJ5Uk/+tZ8jIzsX48eNZRyGEyJFMCh4AJk6ahOsPn+NtPD1hkrU9p26i8zffwNnZmXUUQogcyazgBw0ahDoWFlh/4JysTkGq4PbTMDx+HoXJU6awjkIIkTOZFbyBgQEWLV6Mvy89wKu3ibI6DfkIqVSKZTtOopWXFy1uIkQNyazgAWDs2LGwr1cPv+84KcvTkEr433oGfthb+KxZQ48mIEQNybTgdXV1sXzFCly6G4THz6NkeSryP0QiMZb/eRrdvv0WnTp1Yh2HEMKAhlQqlcryBBKJBO5ubjDQKMalnfOhqSnTf1PIv/acvIFZ6w+Bz+fD3d2ddRxCCAMyb1tNTU1s2boVj59HYdeJG7I+HcG7VatLt5/EqFGjqNwJUWMyH8G/N3HiROzftxf3D/2O+naW8jilWpJKpeg3dT1eJ+bgxcuXtLE2IWpMbtdL1qxZAysra0xatZ82BJGhg2dvI+DJS+zZu5fKnRA1J7eCr1GjBvbt348HQeH488R1eZ1WrcQlZ2Dh1uMYNWoUunXrxjoOIYQxuV2ieW/KlCn4c+dOnNs2B62aN5LnqVVaYXEJvhvvg6x8EZ6HhtLonRAi/4IXCoXo/M03iAx/iZv7FsHW0kyep1dJUqkUoxfvxKX7z3H//gO4ubmxjkQIUQByn7Ooq6uLU6dPQ9/QGEPn/oHC4hJ5R1A5m/66gFPXHuPAgYNU7oSQMkwmpdepUwdn/f0RFZeCySv3Qc4/RKiUS3eDsHznKSxatAgDBw5kHYcQokDkfonmv06cOIFBgwZh/OBvsWrqj7ScvpoeBEfg+2kb0a1bd5w6fZoWkRFCymHaCAMHDsS2bduw4++rWL7zFI3kq+HZizcYNGMzWrdug6PHjlG5E0I+oM06wIQJE1BcXIyZM2dCLJZg6cSBNJL/hEchURg4YyNc3dxx7vx5GBgYsI5ECFFAzAseAGbMmAFNTU1Mnz4dRSVCrJ42BFpaNCKtyO2nYRgyZys8W3rh3PnztAUfIaRSClHwADBt2jQYGBhgwoQJiE1Kx+7fx8HEiEam/7Xn5A3M3XQEXTp3wekzZ2BoaMg6EiFEgTG9yVqRy5cvY/DgH2BV2xhH10xGg3pWrCMxJywVYe7Gw9h/JgBTp07F+vXroa2tMP82E0IUlMIVPABERETAu3dvpKUkYf+K8fi6VVPWkZjJyM7Dz79tx9MXb7Bjxw788ssvrCMRQpSEQhY8AOTk5ODHwYNx9do1TB7SHfPH9IW+ni7rWHJ15X4wpvn8BTG0cPrMGbRr1451JEKIElHYggcAsViMdevWYcmSJXCoa4FtC0bBs2kD1rFkLievAPM3H8Oxi/fQvVs37Nm7F3Xr1mUdixCiZBS64N97+fIlRo4YgUA+X+VH8+9H7QUlpdi8eQtGjBhB00YJIZ9FKQoeAEQiEdavX48lS5bAsrYJ5o3ug8E92qnMdMrnkbH4fccpXH/4HN27dcPuPXtga2vLOhYhRIkpTcG/FxkZiQW//YaTp07Bub4tFo3rj+86uCvtKDc6IQ0rd53GyauP4NioEVauWoXvv/9eab8eQojiULqCf+/p06eYN3cubt66Ba9mjTDxx2/x3Vfu0FGS6YNhbxKw++R1HDp3F3Us6mDpsmUYOXIkTX8khHBGaQseePcc9OvXr2PpkiV48PAhrC1qY0TfjvjZ+ytYW9RiHe8DwlIRzgU8w55Tt/AwOAJWlpaYNn06Jk+eTIuWCCGcU+qC/6/g4GDs2LEDR44cRklJCXq0d0efbzzRpU0zmBqzW84vEonx5MVrXLobhL8vPUR6Vi46deyICRMnom/fvtDR0WGWjRCi2lSm4N/Lzc3FoUOHsH/fPvCDgqCtrYW2bk7o3t4V3du5wcG2jsyvb+fmFyLgyUtcuheMqw+eIytHAMs6dTBw0CCMHz8ejRs3lun5CSEEUMGC/6/4+HicP38e/mfP4uatWxAKhTCrZQI3J/t3/3NxgJszD3Xr1Pqsx+1KpVLk5hfiRVQ8gsNjEBIeg6CIOLyOTQIANG3SGN59+sLb2xstW7akR/oSQuRKpQv+vwQCAQICAvDs2TMEPnuGwMBnSElNAwBoa2uhTm1TWJrXhGVtE1hbmMLEyBA62lrQ1NKERCKBWCxBUYkQaZm5SM7MRWpGHlIzslFQVAwA0NfXh5urKzxbtoSHhwe++uor1K9fn+WXTAhRc2pT8BVJSkpCUFAQ4uPjkZSUhOTkZCQnJyEpMRECgQAikQilpaXQ0tKCtrY29PX1YWVlDWsbG1hbW8Pa2ho2NjZo0qQJXFxcaAYMIUShqHXBE0KIKqOLwoQQoqKo4AkhREVRwRNCiIqigieEEBVFBU8IISqKCp4QQlQUFTwhhKgoKnhCCFFRVPCEEKKiqOAJIURFUcETQoiKooInhBAVRQVPCCEqigqeEEJUFBU8IYSoKCp4QghRUVTwhBCiov4Pe8J23pOZaLIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plot_atoms(Si, rotation='10x,20y,15z')\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a74ff2b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy:  -216.001048427314\n",
      "Stress tensor:\n",
      " [[1.83981555 0.         0.        ]\n",
      " [0.         1.83981555 0.        ]\n",
      " [0.         0.         1.83981555]]\n"
     ]
    }
   ],
   "source": [
    "calc = Abinit(directory=\"work_04\",\n",
    "              ecut= 270,\n",
    "              xc='PBE',\n",
    "              pps='hgh',\n",
    "              toldfe=1.0e-6,\n",
    "              diemac=12,\n",
    "              nstep=20,\n",
    "              kptopt=1,\n",
    "              ngkpt=[2,2,2],\n",
    "              nshiftk=4,\n",
    "              shiftk=[[0.5, 0.5, 0.5],\n",
    "                      [0.5, 0.0, 0.0],\n",
    "                      [0.0, 0.5, 0.0],\n",
    "                      [0.0, 0.0, 0.5]]\n",
    "              )\n",
    "\n",
    "Si.calc = calc\n",
    "print('Energy: ', Si.get_potential_energy())\n",
    "print('Stress tensor:\\n', Si.get_stress(voigt=False)/un.GPa)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "382b1177",
   "metadata": {},
   "source": [
    "### Convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "251feb65",
   "metadata": {},
   "source": [
    "We will be scanning multiple parameters in our study. To avoid repeating the same code over and over again, we will write (somewhat) universal scanning function. This function uses some more advanced Python techniques like doc-strings inside the function, basic functional programming (`lambda` expression defining anonymous function and function assignment) and passing arguments as dictionaries. You can read about these topics in the python documentation. You can experiment with this function and use it as a basis for your own tools.\n",
    "\n",
    "To avoid cluttering the examples, all auxiliary functions are collected in one file: `abilib.py` which we import here and in other notebooks. The relevant functions are printed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "381003d7-d12c-4576-b10f-6fe594ae2670",
   "metadata": {},
   "outputs": [],
   "source": [
    "from abilib import set_params, scan_param"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "170978ed-9a46-4503-945c-3e881e9a88b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m  \u001b[0mset_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0mset_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m'''\u001b[0m\n",
       "\u001b[0;34m    Create new calculator with modified parameters based on the `calc` argument.\u001b[0m\n",
       "\u001b[0;34m    '''\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mparams\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparameters\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mparams\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0mAbinit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdirectory\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdirectory\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Project/abinitio-methods/abilib.py\n",
       "\u001b[0;31mType:\u001b[0m      function"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "?? set_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "84c94a8e-3dcd-4eb5-b63a-13f121f68d7b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m\n",
       " \u001b[0mscan_param\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mcryst\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mparam\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mstart\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mstop\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0msteps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mlista\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mlog_scale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mint_par\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mSource:\u001b[0m   \n",
       "\u001b[0;32mdef\u001b[0m \u001b[0mscan_param\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcryst\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparam\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstart\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstop\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m               \u001b[0mlista\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlog_scale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mint_par\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m'''\u001b[0m\n",
       "\u001b[0;34m    Scan a range of parameter param and plot \u001b[0m\n",
       "\u001b[0;34m    '''\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mfi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mfo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mif\u001b[0m \u001b[0mlista\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mlog_scale\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mfi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mfo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mlst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32melse\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mint_par\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mlst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlista\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mdat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mvs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mfor\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlst\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'{param}='\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m''\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mlista\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'{int(fo(v)):<8d}'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mint_par\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34mf'{fo(v):<8.2f}'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m': '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'{fo(v)}'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mend\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m': '\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32mif\u001b[0m \u001b[0mint_par\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mset_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mparam\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0;32melse\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m            \u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mset_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mparam\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mvs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfo\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mdat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_potential_energy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m                    \u001b[0mcryst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_stress\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m        \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'{dat[-1][0]:.3f} eV ; {dat[-1][1]/un.GPa:6.3f} GPa'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mdat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0;32mreturn\u001b[0m \u001b[0mvs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Project/abinitio-methods/abilib.py\n",
       "\u001b[0;31mType:\u001b[0m      function"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "?? scan_param"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c0fa09-5243-4bd1-83be-ebdc3af8944a",
   "metadata": {},
   "source": [
    "### Energy cutoff convergence\n",
    "\n",
    "The energy cut-off should be the first parameter to check for convergence. It has a strong influence on both: calculation efficiency and accuracy. It should be selected at the level which provides stable, converged energy differences (not absolute energy) and forces/stresses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d929d25c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ecut=250.00  : -215.933 eV ;  2.410 GPa\n",
      "ecut=280.62  : -216.026 eV ;  1.606 GPa\n",
      "ecut=314.98  : -216.102 eV ;  1.005 GPa\n",
      "ecut=353.55  : -216.148 eV ;  0.578 GPa\n",
      "ecut=396.85  : -216.172 eV ;  0.346 GPa\n",
      "ecut=445.45  : -216.185 eV ;  0.235 GPa\n",
      "ecut=500.00  : -216.192 eV ;  0.187 GPa\n"
     ]
    }
   ],
   "source": [
    "ecuts, eres = scan_param(Si, 'ecut', 250, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "cc743864",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAINCAYAAAA+3HWvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbZUlEQVR4nO3dfVwU9d4//tdyt6DCKiGwCCI34Ykwb1ATyxu8xbvq68mjx1tCKT1ilnZneQ5qKZXk5ZV1UMuDcjxll6klWYYlYpaKiFyKN6iIQrCEii6KciP7+f3hz71cWXAXFvZmXs/HYx6HnfnMzHs/Z2pfzXxmRiaEECAiIiKSADtzF0BERETUWhh8iIiISDIYfIiIiEgyGHyIiIhIMhh8iIiISDIYfIiIiEgyGHyIiIhIMhh8iIiISDIczF0A/R+NRoOSkhK4urpCJpOZuxwiIiKrIYTAjRs34OPjAzu7hs/rMPhYkJKSEvj5+Zm7DCIiIqtVVFQEX1/fBpcz+FgQV1dXAHf/T3NzczNzNURERNajoqICfn5+2t/ShjD4WJB7l7fc3NwYfIiIiJrgYUNFOLiZiIiIJINnfGxYnUYgs6AcZTeq4OnqjL4B7rC346BpIiKSLgYfG7U7V4WlqaegUldp5ykVzogfF4qoMKUZKyMiIjIfXuqyQbtzVZizOVsn9ABAqboKczZnY3euykyVERERmReDj42p0wgsTT0FoWfZvXlLU0+hTqOvBRERkW1j8LExmQXl9c703E8AUKmrkFlQ3npFERERWQgGHxtTdqPh0NOUdkRERLaEwcfGeLo6m7QdERGRLWHwsTF9A9yhVDijoZvWZbh7d1ffAPfWLIuIiMgiMPjYGHs7GeLHhQJAvfBz73P8uFA+z4eIiCSJwccGRYUpkTS1F7wVupezvBXOSJrai8/xISIiyeIDDG1UVJgSw0O9+eRmIiKi+zD42DB7Oxkigh4xdxlEREQWg5e6iIiISDIYfIiIiEgyrCL4XLx4ETNnzkRAQABcXFwQFBSE+Ph41NTU6LSbP38+wsPDIZfL0aNHD73bkclk9abdu3c3uv/s7GwMHz4c7du3xyOPPIIXX3wRN2/e1Gmjb7tr165t9ncnIiIi07GKMT5nzpyBRqPBunXrEBwcjNzcXMTGxqKyshKJiYnadkIIxMTE4PDhwzh+/HiD2/vpp5/w+OOPaz+7uzf8TJuSkhIMGzYMEydOxCeffIKKigq88soriI6Oxtdff63TNjk5GVFRUdrPCoWiKV+XiIiIWohVBJ+oqCidQBEYGIi8vDwkJSXpBJ+PP/4YAHD58uVGg88jjzwCb29vg/b93XffwdHREZ9++ins7O6eIPv000/Rs2dPnD9/HsHBwdq27du3N3i7RERE1Pqs4lKXPmq1utEzNY155pln4OnpiaeeeqreWZsHVVdXw8nJSRt6AMDFxQUAcODAAZ22cXFx8PDwQJ8+fbB27VpoNJqHbruiokJnIiIiopZjlcEnPz8fa9aswezZs41ar127dli1ahW+/vprfP/99xg6dCgmTpyIzZs3N7jOkCFDUFpaipUrV6KmpgbXrl3D22+/DQBQqVTadu+++y62bt2Kn376CZMmTcLChQuxYsWKRutJSEiAQqHQTn5+fkZ9HyIiIjKSMKP4+HgBoNHpyJEjOusUFxeL4OBgMXPmzEa32717d4NqiIuLE926dWu0zX/+8x/h5eUl7O3thZOTk3jttdeEl5eX+OCDDxpcJzExUbi5uTW63aqqKqFWq7VTUVGRACDUarVBtRMREdFdarXaoN9Qs47xiYuLw6RJkxpt06VLF+3fJSUliIyMREREBNavX2+SGvr164fPP/+80TaTJ0/G5MmT8ccff6Bt27aQyWRYtWoVAgICGt1uRUUF/vjjD3h5eeltI5fLIZfLm1U/ERERGc6swcfDwwMeHh4GtS0uLkZkZCTCw8ORnJysM+amOY4dOwal0rB3V90LMP/617/g7OyM4cOHN7pdZ2dntG/f3hRlEhERkQlYxV1dJSUlGDx4MDp37ozExERcvnxZu+z+u6jOnz+PmzdvorS0FLdv30ZOTg4AIDQ0FE5OTti0aRMcHR3Rs2dP2NnZITU1FR9//DE++OAD7TYyMzMxffp0/Pzzz+jUqRMA4JNPPkH//v3Rrl077NmzB6+//jref/99bahJTU1FaWkpIiIi4OLigvT0dLzzzjt48cUXeUaHiIjIglhF8ElLS8P58+dx/vx5+Pr66iwTQmj/njVrFjIyMrSfe/bsCQAoKCjQXjJ77733cOnSJdjb2yMkJAT/+te/MHXqVO06t27dQl5eHmpra7XzMjMzER8fj5s3b+JPf/oT1q1bh2nTpmmXOzo64p///CcWLFgAjUaDwMBALFu2DHPnzjVpPxAREVHzyMT9yYHMqqKiAgqFAmq1Gm5ubuYuh4iIyGoY+htqlbezExERETUFgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJhlUEn4sXL2LmzJkICAiAi4sLgoKCEB8fj5qaGp128+fPR3h4OORyOXr06KF3W0IIJCYmIiQkBHK5HH5+flixYkWj+7927RqmTZsGhUIBhUKBadOm4fr16zptCgsLMW7cOLRt2xYeHh54+eWX69VHRERE5uVg7gIMcebMGWg0Gqxbtw7BwcHIzc1FbGwsKisrkZiYqG0nhEBMTAwOHz6M48eP693W/PnzkZaWhsTERHTr1g1qtRpXrlxpdP+TJ0/G77//jt27dwMAXnzxRUybNg2pqakAgLq6OowZMwYdO3bEgQMHcPXqVcyYMQNCCKxZs8ZEvUBERETNJRNCCHMX0RQrV65EUlISLly4UG/ZkiVL8M033yAnJ0dn/unTp/HEE08gNzcXXbt2NWg/p0+fRmhoKA4dOoQnn3wSAHDo0CFERETgzJkz6Nq1K3744QeMHTsWRUVF8PHxAQBs2bIF0dHRKCsrg5ubm0H7qqiogEKhgFqtNngdIiIiMvw31CoudemjVqvh7u5u1DqpqakIDAzEd999h4CAAHTp0gWzZs1CeXl5g+scPHgQCoVCG3oAoF+/flAoFPjtt9+0bcLCwrShBwBGjhyJ6upqHD16tMFtV1dXo6KiQmciIiKilmOVwSc/Px9r1qzB7NmzjVrvwoULuHTpErZu3YqUlBRs3LgRR48exfPPP9/gOqWlpfD09Kw339PTE6Wlpdo2Xl5eOss7dOgAJycnbRt9EhIStOOGFAoF/Pz8jPo+REREZByzBp8lS5ZAJpM1OmVlZemsU1JSgqioKEyYMAGzZs0yan8ajQbV1dVISUnBgAEDMHjwYGzYsAHp6enIy8trcD2ZTFZvnhBCZ74hbR60aNEiqNVq7VRUVGTU9yEiIiLjmHVwc1xcHCZNmtRomy5dumj/LikpQWRkJCIiIrB+/Xqj96dUKuHg4ICQkBDtvMceewzA3buy9I378fb2xh9//FFv/uXLl7Vneby9vXH48GGd5deuXUNtbW29M0H3k8vlkMvlRn8PIiIiahqzBh8PDw94eHgY1La4uBiRkZEIDw9HcnIy7OyMP1n11FNP4c6dO8jPz0dQUBAA4OzZswAAf39/vetERERArVYjMzMTffv2BQAcPnwYarUa/fv317ZZvnw5VCoVlEolACAtLQ1yuRzh4eFG10lEREQtwyru6iopKcGgQYPQuXNnpKSkwN7eXrvM29tb+/f58+dx8+ZNrF27Funp6fjqq68AAKGhoXBycoJGo0GfPn3Qrl07rF69GhqNBnPnzoWbmxvS0tIAAJmZmZg+fTp+/vlndOrUCQAwatQolJSUYN26dQDu3s7u7++vczt7jx494OXlhZUrV6K8vBzR0dF47rnnjLqdnXd1ERERNY3Bv6HCCiQnJwsAeqf7DRo0SG+bgoICbZvi4mIxfvx40a5dO+Hl5SWio6PF1atXtcvT09PrrXP16lUxZcoU4erqKlxdXcWUKVPEtWvXdPZ96dIlMWbMGOHi4iLc3d1FXFycqKqqMup7qtVqAUCo1Wqj1iMiIpI6Q39DreKMj1TwjA8REVHTGPobatQYHyEEMjIy8Msvv+DixYu4desWOnbsiJ49e2LYsGG8HZuIiIgsmkEjhG/fvo0VK1bAz88Po0aNwq5du3D9+nXY29vj/PnziI+PR0BAAEaPHo1Dhw61dM1ERERETWLQGZ+QkBA8+eSTWLt2LUaOHAlHR8d6bS5duoQvvvgCEydOxOLFixEbG2vyYomIiIiaw6AxPrm5uQgLCzNogzU1Nbh06RIeffTRZhcnNRzjQ0RE1DQmfVdXWFhYvRd+NsTJyYmhh4iIiCySwU8B7NWrF8LDw5GUlAS1Wt2SNRERERG1CIODz6+//opevXrhrbfeglKpxNSpU5Gent6StRERERGZlMHBJyIiAp999hlKS0uRlJSE33//HcOGDUNQUBCWL1+O33//vSXrJEKdRuBg/lV8m1OMg/lXUafhI6iIiMg4zXqAYX5+PpKTk5GSkgKVSoXhw4fj+++/N2V9ksLBzQ3bnavC0tRTUKmrtPOUCmfEjwtFVJjSjJUREZElMPQ3tNlPbr558yb+85//4O2338b169dRV1fXnM1JGoOPfrtzVZizORsPHqiy//9/k6b2YvghIpI4k97VpU9GRgZmzJgBb29vvPHGGxg/fjx+/fXXpm6OSK86jcDS1FP1Qg8A7bylqad42YuIiAxi1CsrioqKsHHjRmzcuBEFBQXo378/1qxZg7/85S9o27ZtS9VIEpZZUK5zeetBAoBKXYXMgnJEBD3SeoUREZFVMjj4DB8+HOnp6ejYsSOmT5+OmJgYdO3atSVrI0LZjYZDT1PaERGRtBkcfFxcXLBt2zaMHTsW9vb2LVkTkZanq7NJ2xERkbQZHHx27typ8/n8+fPIz8/HwIED4eLiAiEEZDJZA2sTNU3fAHcoFc4oVVfpHecjA+CtcEbfAPfWLo2IiKyQ0YObr169iqFDhyIkJASjR4+GSqUCAMyaNQsLFy40eYEkbfZ2MsSPCwXwf3dx3XPvc/y4UNjbMXQTEdHDGR18Xn31VTg6OqKwsBBt2rTRzp84cSJ2795t0uKIACAqTImkqb3grdC9nOWtcOat7EREZBSj7uoCgLS0NPz444/w9fXVmf/oo4/i0qVLJiuM6H5RYUoMD/VGZkE5ym5UwdP17uUtnukhIiJjGB18Kisrdc703HPlyhXI5XKTFEWkj72djLesExFRsxh9qWvgwIFISUnRfpbJZNBoNFi5ciUiIyNNWhwRERGRKRl9xmflypUYPHgwsrKyUFNTgzfeeAMnT55EeXk5n9xMREREFs3oMz6hoaE4fvw4+vbti+HDh6OyshLjx4/HsWPHEBQU1BI1EhEREZlEs19SSqbDl5QSERE1jUlfUlpYWGjUzouLi41qT0RERNQaDAo+ffr0QWxsLDIzMxtso1ar8dlnnyEsLAzbt283WYFEREREpmJQ8Dl9+jQUCgWioqLg5eWFMWPGIDY2FvPmzcPUqVPRq1cveHp6YuPGjVi5ciXmzZtn0iIvXryImTNnIiAgAC4uLggKCkJ8fDxqamp02s2fPx/h4eGQy+Xo0aOH3m0JIZCYmIiQkBDI5XL4+flhxYoVje7/2rVrmDZtGhQKBRQKBaZNm4br16/rtJHJZPWmtWvXNudrExERkYkZdFeXu7s7EhMT8d577+H777/HL7/8gosXL+L27dvw8PDAlClTMHLkSISFhbVIkWfOnIFGo8G6desQHByM3NxcxMbGorKyEomJidp2QgjExMTg8OHDOH78uN5tzZ8/H2lpaUhMTES3bt2gVqtx5cqVRvc/efJk/P7779onU7/44ouYNm0aUlNTddolJycjKipK+1mhUDT1KxMREVELsNrBzStXrkRSUhIuXLhQb9mSJUvwzTffICcnR2f+6dOn8cQTTyA3Nxddu3Y1aD+nT59GaGgoDh06hCeffBIAcOjQIURERODMmTPa7chkMuzYsQPPPfdck78TBzcTERE1jUkHN1sitVoNd3fj3sidmpqKwMBAfPfddwgICECXLl0wa9YslJeXN7jOwYMHoVAotKEHAPr16weFQoHffvtNp21cXBw8PDzQp08frF27FhqNxrgvRURERC3K6AcYWoL8/HysWbMGH330kVHrXbhwAZcuXcLWrVuRkpKCuro6vPrqq3j++eexd+9eveuUlpbC09Oz3nxPT0+UlpZqP7/77rsYOnQoXFxc8PPPP2PhwoW4cuUKFi9e3GA91dXVqK6u1n6uqKgw6vsQERGRccwafJYsWYKlS5c22ubIkSPo3bu39nNJSQmioqIwYcIEzJo1y6j9aTQaVFdXIyUlBSEhIQCADRs2IDw8HHl5eQ1e/pLJ6r8IUwihM//+gHNvYPWyZcsaDT4JCQkP/f5ERERkOmYNPnFxcZg0aVKjbbp06aL9u6SkBJGRkYiIiMD69euN3p9SqYSDg4M29ADAY489BuDus4r0BR9vb2/88ccf9eZfvnwZXl5eDe6rX79+qKiowB9//NFgu0WLFmHBggXazxUVFfDz8zP4+xAREZFxmvR29rZt25pk5x4eHvDw8DCobXFxMSIjIxEeHo7k5GTY2Rk/POmpp57CnTt3kJ+fr329xtmzZwEA/v7+eteJiIiAWq1GZmYm+vbtCwA4fPgw1Go1+vfv3+C+jh07BmdnZ7Rv377BNnK5nG+0JyIiakVGpwcvLy/ExMTgwIEDLVGPXiUlJRg8eDD8/PyQmJiIy5cvo7S0VGeMDQCcP38eOTk5KC0txe3bt5GTk4OcnBzt836GDRuGXr16ISYmBseOHcPRo0fx0ksvYfjw4dqzQJmZmfjTn/6kffr0Y489hqioKMTGxuLQoUM4dOgQYmNjMXbsWO0ZotTUVHz22WfIzc1Ffn4+Pv/8c7zzzjt48cUXGWyIiIgsiTDSzp07xfjx44WTk5N49NFHRUJCgiguLjZ2M0ZJTk4WAPRO9xs0aJDeNgUFBdo2xcXFYvz48aJdu3bCy8tLREdHi6tXr2qXp6en11vn6tWrYsqUKcLV1VW4urqKKVOmiGvXrmmX//DDD6JHjx6iXbt2ok2bNiIsLEysXr1a1NbWGvU91Wq1ACDUarVR6xEREUmdob+hTX6Oz9WrV5GSkoKNGzfi1KlTGDlyJGJiYvDMM8/AwcEqbxYzOz7Hh4iIqGkM/Q01yQMM16xZg9dffx01NTXw8PDA7Nmz8dZbb6FNmzbN3bSkMPgQERE1jaG/oU0+NVNaWoqUlBQkJyejsLAQzz//PGbOnImSkhK8//77OHToENLS0pq6eSIiIiKTMzr4bN++HcnJyfjxxx8RGhqKuXPnYurUqTp3L/Xo0QM9e/Y0ZZ1EREREzWZ08HnhhRcwadIk/Prrr+jTp4/eNoGBgXjnnXeaXRwRERGRKRk9xufWrVscu9NCOMaHiIioaVpsjM+dO3f0vlNKJpNBLpfDycnJ2E0SERERtQqjg0/79u31vrvqHl9fX0RHRyM+Pr5JT1cmIiIiailGB5+NGzfinXfeQXR0NPr27QshBI4cOYJNmzZh8eLFuHz5MhITEyGXy/H222+3RM1ERERETWJ08Nm0aRM++ugj/OUvf9HOe+aZZ9CtWzesW7cOP//8Mzp37ozly5cz+BAREZFFMfpa1MGDB/Xeqt6zZ08cPHgQAPD000+jsLCw+dURERERmZDRwcfX1xcbNmyoN3/Dhg3w8/MDcPd1Fh06dGh+dUREREQmZPSlrsTEREyYMAE//PAD+vTpA5lMhiNHjuDMmTP4+uuvAQBHjhzBxIkTTV4sERERUXM06V1dly5dwtq1a5GXlwchBP70pz/hpZdeQpcuXVqgROngc3yIiIiapkWe41NbW4sRI0Zg3bp1SEhIaHaRRERERK3JqDE+jo6OyM3NbfQ5PkRERESWyujBzdOnT9c7uJmIiIjI0hk9uLmmpgaff/459uzZg969e6Nt27Y6y1etWmWy4oiIiIhMyejgk5ubi169egEAzp49q7OMl8CIiIjIkhkdfNLT01uiDiIiIqIW1+S3iJ4/fx4//vgjbt++DQBowl3xRERERK3K6OBz9epVDB06FCEhIRg9ejRUKhUAYNasWVi4cKHJCyQiIiIyFaODz6uvvgpHR0cUFhaiTZs22vkTJ07E7t27TVocERERkSkZPcYnLS0NP/74I3x9fXXmP/roo7h06ZLJCiMiIiIyNaPP+FRWVuqc6bnnypUrkMvlJimKiIiIqCUYHXwGDhyIlJQU7WeZTAaNRoOVK1ciMjLSpMURERERmZLRl7pWrlyJwYMHIysrCzU1NXjjjTdw8uRJlJeX49dff22JGomIiIhMwugzPqGhoTh+/Dj69u2L4cOHo7KyEuPHj8exY8cQFBTUEjXi4sWLmDlzJgICAuDi4oKgoCDEx8ejpqZGp938+fMRHh4OuVyOHj166N2WEAKJiYkICQmBXC6Hn58fVqxY0ej+ly9fjv79+6NNmzZo37693jaFhYUYN24c2rZtCw8PD7z88sv16iMiIiLzMvqMDwB4e3tj6dKlpq6lQWfOnIFGo8G6desQHByM3NxcxMbGorKyEomJidp2QgjExMTg8OHDOH78uN5tzZ8/H2lpaUhMTES3bt2gVqtx5cqVRvdfU1ODCRMmICIiQu97yurq6jBmzBh07NgRBw4cwNWrVzFjxgwIIbBmzZrmfXmSnDqNQGZBOcpuVMHT1Rl9A9xhb8enohMRmYJMNOHJg9evX0dmZibKysqg0Wh0lk2fPt1kxTVm5cqVSEpKwoULF+otW7JkCb755hvk5OTozD99+jSeeOIJ5ObmomvXrkbvc+PGjXjllVdw/fp1nfk//PADxo4di6KiIvj4+AAAtmzZgujoaJSVlcHNzc2g7VdUVEChUECtVhu8DtmW3bkqLE09BZW6SjtPqXBG/LhQRIUpzVgZEZFlM/Q31OgzPqmpqZgyZQoqKyvh6uqq834umUzWasFHrVbD3d3dqHVSU1MRGBiI7777DlFRURBCYNiwYfjwww+N3tb9Dh48iLCwMG3oAYCRI0eiuroaR48e5aBvMsjuXBXmbM7Gg/8lUqquwpzN2Uia2ovhh4iomYwe47Nw4ULExMTgxo0buH79Oq5du6adysvLW6LGevLz87FmzRrMnj3bqPUuXLiAS5cuYevWrUhJScHGjRtx9OhRPP/8882qp7S0FF5eXjrzOnToACcnJ5SWlja4XnV1NSoqKnQmkqY6jcDS1FP1Qg8A7bylqadQp+GrYYiImsPo4FNcXIyXX35Z77N8jLVkyRLIZLJGp6ysLJ11SkpKEBUVhQkTJmDWrFlG7U+j0aC6uhopKSkYMGAABg8ejA0bNiA9PR15eXnN+i763kwvhGj0jfUJCQlQKBTayc/Pr1k1kPXKLCjXubz1IAFApa5CZkHr/McFEZGtMvpS18iRI5GVlYXAwMBm7zwuLg6TJk1qtE2XLl20f5eUlCAyMhIRERFYv3690ftTKpVwcHBASEiIdt5jjz0G4O5dWU0Z9wPcHex9+PBhnXnXrl1DbW1tvTNB91u0aBEWLFig/VxRUcHwI1FlNxoOPU1pR0RE+hkdfMaMGYPXX38dp06dQrdu3eDo6Kiz/JlnnjF4Wx4eHvDw8DCobXFxMSIjIxEeHo7k5GTY2Rn/YvmnnnoKd+7cQX5+vvbW+7NnzwIA/P39jd7ePREREVi+fDlUKhWUyrtjMNLS0iCXyxEeHt7genK5nE+7JgCAp6uzSdsREZF+Rgef2NhYAMCyZcvqLZPJZKirq2t+VQ8oKSnB4MGD0blzZyQmJuLy5cvaZd7e3tq/z58/j5s3b6K0tBS3b9/W3tUVGhoKJycnDBs2DL169UJMTAxWr14NjUaDuXPnYvjw4dqzQJmZmZg+fTp+/vlndOrUCcDds0Hl5eUoLCxEXV2ddrvBwcFo164dRowYgdDQUEybNg0rV65EeXk5XnvtNcTGxvLuLDJI3wB3KBXOKFVX6R3nIwPgrbh7azsRETWd0cHnwdvXW0NaWhrOnz+P8+fP13s56v1348+aNQsZGRnazz179gQAFBQUoEuXLrCzs0NqairmzZuHgQMHom3bthg1ahQ++ugj7Tq3bt1CXl4eamtrtfP+8Y9/YNOmTfW2m56ejsGDB8Pe3h67du3C3/72Nzz11FNwcXHB5MmTdZ4xRNQYezsZ4seFYs7mbMgAnfBzb5RY/LhQPs+HiKiZmvQcH2oZfI4P8Tk+RERNY/Ln+IwePRpffvklFAoFgLuvcZg7d672FQ5Xr17FgAEDcOrUqeZVTiRhUWFKDA/15pObiYhaiMFnfOzt7aFSqeDp6QkAcHNzQ05Ojvburj/++AM+Pj4tMsZHKnjGh4iIqGkM/Q01+NaoB/MRr5ARERGRtTH+nnAiIiIiK2Vw8Ln3JOUH5xERERFZC4MHNwshEB0drX3gXlVVFWbPno22bdsCuPveKSIiIiJLZnDwmTFjhs7nqVOn1mvTWm9mJyIiImoKg4NPcnJyS9ZBRERE1OI4uJmIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkw8HcBRARPahOI5BZUI6yG1XwdHVG3wB32NvJzF0WEdkABh8isii7c1VYmnoKKnWVdp5S4Yz4caGIClOasTIisgVWcanr4sWLmDlzJgICAuDi4oKgoCDEx8ejpqZGp938+fMRHh4OuVyOHj166N2WEAKJiYkICQmBXC6Hn58fVqxY0ej+ly9fjv79+6NNmzZo37693jYymazetHbt2qZ8XSLJ2p2rwpzN2TqhBwBK1VWYszkbu3NVZqqMiGyFVZzxOXPmDDQaDdatW4fg4GDk5uYiNjYWlZWVSExM1LYTQiAmJgaHDx/G8ePH9W5r/vz5SEtLQ2JiIrp16wa1Wo0rV640uv+amhpMmDABERER2LBhQ4PtkpOTERUVpf2sUCiM/KZE0lWnEViaegpCzzIBQAZgaeopDA/15mUvImoyqwg+UVFROoEiMDAQeXl5SEpK0gk+H3/8MQDg8uXLeoPP6dOnkZSUhNzcXHTt2tXg/S9duhQAsHHjxkbbtW/fHt7e3gZvl4j+T2ZBeb0zPfcTAFTqKmQWlCMi6JHWK4yIbIpVXOrSR61Ww93d3ah1UlNTERgYiO+++w4BAQHo0qULZs2ahfLycpPUFBcXBw8PD/Tp0wdr166FRqNptH11dTUqKip0JiKpKrvRcOhpSjsiIn2sMvjk5+djzZo1mD17tlHrXbhwAZcuXcLWrVuRkpKCjRs34ujRo3j++eebXdO7776LrVu34qeffsKkSZOwcOHCh44dSkhIgEKh0E5+fn7NroPIWnm6Opu0HRGRPmYNPkuWLNE7KPj+KSsrS2edkpISREVFYcKECZg1a5ZR+9NoNKiurkZKSgoGDBiAwYMHY8OGDUhPT0deXl6zvsvixYsRERGBHj16YOHChVi2bBlWrlzZ6DqLFi2CWq3WTkVFRc2qgcia9Q1wh1LhjIZG78hw9+6uvgHGneklIrqfWcf4xMXFYdKkSY226dKli/bvkpISREZGIiIiAuvXrzd6f0qlEg4ODggJCdHOe+yxxwAAhYWFRo37eZh+/fqhoqICf/zxB7y8vPS2kcvlkMvlJtsnkTWzt5Mhflwo5mzOhgzQGeR8LwzFjwvlwGYiahazBh8PDw94eHgY1La4uBiRkZEIDw9HcnIy7OyMP1n11FNP4c6dO8jPz0dQUBAA4OzZswAAf39/o7fXmGPHjsHZ2bnB29+JqL6oMCWSpvaq9xwfbz7Hh4hMxCru6iopKcHgwYPRuXNnJCYm4vLly9pl999Fdf78edy8eROlpaW4ffs2cnJyAAChoaFwcnLCsGHD0KtXL8TExGD16tXQaDSYO3cuhg8frj0LlJmZienTp+Pnn39Gp06dANw9G1ReXo7CwkLU1dVptxscHIx27dohNTUVpaWliIiIgIuLC9LT0/HOO+/gxRdf5BkdIiNFhSkxPNSbT24mopYhrEBycrLA3TPf9ab7DRo0SG+bgoICbZvi4mIxfvx40a5dO+Hl5SWio6PF1atXtcvT09PrrTNjxgy9201PTxdCCPHDDz+IHj16iHbt2ok2bdqIsLAwsXr1alFbW2vU91Sr1QKAUKvVRvcRERGRlBn6GyoTQuh7XhiZQUVFBRQKBdRqNdzc3MxdDhERkdUw9DfUKm9nJyIiImoKBh8iIiKSDAYfIiIikgwGHyIiIpIMq7idnYiIdNVpBG/5J2oCBh8iIiuzO1dV7yGPSj7kkcggvNRFRGRFdueqMGdztk7oAYBSdRXmbM7G7lyVmSojsg4MPkREVqJOI7A09RT0PXzt3rylqadQp+Hj2YgawuBDRGQlMgvK653puZ8AoFJXIbOgvPWKIrIyDD5ERFai7EbDoacp7YikiMGHiMhKeLo6m7QdkRQx+BARWYm+Ae5QKpzR0E3rMty9u6tvgHtrlkVkVRh8iIishL2dDPHjQgGgXvi59zl+XCif50PUCAYfIiIrEhWmRNLUXvBW6F7O8lY4I2lqLz7Hh+gh+ABDIiIrExWmxPBQbz652QT4BGzpYfAhIrJC9nYyRAQ9Yu4yrBqfgC1NvNRFRESSwydgSxeDDxERSQqfgC1tDD5ERCQpfAK2tHGMDxERSQqfgG1e5h5QzuBDRESSwidgm48lDCjnpS4iIpIUPgHbPCxlQDmDDxERSQqfgN36LGlAOYMPERFJDp+A3bosaUA5x/gQEZEk8QnYrceSBpQz+BARkWTxCditw5IGlFvFpa6LFy9i5syZCAgIgIuLC4KCghAfH4+amhqddvPnz0d4eDjkcjl69Oihd1tCCCQmJiIkJARyuRx+fn5YsWJFs/ddWFiIcePGoW3btvDw8MDLL79crw0REZEUWdKAcqs443PmzBloNBqsW7cOwcHByM3NRWxsLCorK5GYmKhtJ4RATEwMDh8+jOPHj+vd1vz585GWlobExER069YNarUaV65cada+6+rqMGbMGHTs2BEHDhzA1atXMWPGDAghsGbNGtN2BhERkZW5N6B8zuZsyACdQc6tPaBcJoSwymdyr1y5EklJSbhw4UK9ZUuWLME333yDnJwcnfmnT5/GE088gdzcXHTt2tVk+/7hhx8wduxYFBUVwcfHBwCwZcsWREdHo6ysDG5ubgZtt6KiAgqFAmq12uB1iIiIrEVLPsfH0N9Qqzjjo49arYa7u3GnxFJTUxEYGIjvvvsOUVFREEJg2LBh+PDDD43a1oP7PnjwIMLCwrShBwBGjhyJ6upqHD16FJGRkXq3U11djerqau3niooKo74PERGRNbGEAeVWMcbnQfn5+VizZg1mz55t1HoXLlzApUuXsHXrVqSkpGDjxo04evQonn/++Wbtu7S0FF5eXjrtOnToACcnJ5SWlja4rYSEBCgUCu3k5+dn1PchIiKyNvcGlD/boxMigh5p9bvozBp8lixZAplM1uiUlZWls05JSQmioqIwYcIEzJo1y6j9aTQaVFdXIyUlBQMGDMDgwYOxYcMGpKenIy8v76HrN7Zvmaz+/3FCCL3z71m0aBHUarV2KioqMur7EBERkXHMeqkrLi4OkyZNarRNly5dtH+XlJQgMjISERERWL9+vdH7UyqVcHBwQEhIiHbeY489BuDuXVmNjftpbN/e3t44fPiwzrxr166htra23pmg+8nlcsjlcqO/BxERETWNWYOPh4cHPDw8DGpbXFyMyMhIhIeHIzk5GXZ2xp+seuqpp3Dnzh3k5+cjKCgIAHD27FkAgL+/f5P3HRERgeXLl0OlUkGpvDs4Ky0tDXK5HOHh4QbXd2+cOcf6EBERGefeb+dD79kSVqC4uFgEBweLIUOGiN9//12oVCrtdL9z586JY8eOiZdeekmEhISIY8eOiWPHjonq6mohhBB1dXWiV69eYuDAgSI7O1tkZWWJJ598UgwfPly7jcOHD4uuXbuK33//3eB937lzR4SFhYmhQ4eK7Oxs8dNPPwlfX18RFxdn1PcsKioSuHuXHydOnDhx4sSpCVNRUVGjv7VWcVdXWloazp8/j/Pnz8PX11dnmbgv2c2aNQsZGRnazz179gQAFBQUoEuXLrCzs0NqairmzZuHgQMHom3bthg1ahQ++ugj7Tq3bt1CXl4eamtrDd63vb09du3ahb/97W946qmn4OLigsmTJ+s8Y8gQPj4+KCoqgqura6Njg4xRUVEBPz8/FBUV8RZ5E2B/mh771LTYn6bHPjW9luhTIQRu3Lihc4e1Plb7HB8yDJ8NZFrsT9Njn5oW+9P02KemZ84+tcrb2YmIiIiagsGHiIiIJIPBx8bJ5XLEx8fztnkTYX+aHvvUtNifpsc+NT1z9inH+BAREZFk8IwPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDj5VJSEhAnz594OrqCk9PTzz33HP13iwfHR1d7y33/fr102lTXV2NefPmwcPDA23btsUzzzyD33//vTW/isVISkrCE088ATc3N7i5uSEiIgI//PCDdrkQAkuWLIGPjw9cXFwwePBgnDx5Umcb7E9dD+tTHqPNk5CQAJlMhldeeUU7j8dp0+nrTx6jxlmyZEm9/vL29tYut6Tjk8HHymRkZGDu3Lk4dOgQ9uzZgzt37mDEiBGorKzUaRcVFQWVSqWdvv/+e53lr7zyCnbs2IEtW7bgwIEDuHnzJsaOHYu6urrW/DoWwdfXF++//z6ysrKQlZWFIUOG4Nlnn9X+Q/nhhx9i1apV+OSTT3DkyBF4e3tj+PDhuHHjhnYb7E9dD+tTgMdoUx05cgTr16/HE088oTOfx2nTNNSfAI9RYz3++OM6/XXixAntMos6Po16iyZZnLKyMgFAZGRkaOfNmDFDPPvssw2uc/36deHo6Ci2bNminVdcXCzs7OzE7t27W7Jcq9GhQwfx+eefC41GI7y9vcX777+vXVZVVSUUCoVYu3atEIL9aah7fSoEj9GmunHjhnj00UfFnj17xKBBg8T8+fOFEILHaRM11J9C8Bg1Vnx8vOjevbveZZZ2fPKMj5VTq9UAAHd3d535+/btg6enJ0JCQhAbG4uysjLtsqNHj6K2thYjRozQzvPx8UFYWBh+++231incQtXV1WHLli2orKxEREQECgoKUFpaqtNXcrkcgwYN0vYV+7NxD/bpPTxGjTd37lyMGTMGw4YN05nP47RpGurPe3iMGufcuXPw8fFBQEAAJk2ahAsXLgCwvOPTKt7OTvoJIbBgwQI8/fTTCAsL084fNWoUJkyYAH9/fxQUFODvf/87hgwZgqNHj0Iul6O0tBROTk7o0KGDzva8vLxQWlra2l/DIpw4cQIRERGoqqpCu3btsGPHDoSGhmr/gfPy8tJp7+XlhUuXLgEA+7MBDfUpwGO0KbZs2YLs7GwcOXKk3rJ7fcLj1HCN9SfAY9RYTz75JFJSUhASEoI//vgD7733Hvr374+TJ09a3PHJ4GPF4uLicPz4cRw4cEBn/sSJE7V/h4WFoXfv3vD398euXbswfvz4BrcnhIBMJmuxei1Z165dkZOTg+vXr2Pbtm2YMWMGMjIytMsf7BdD+krK/Qk03KehoaE8Ro1UVFSE+fPnIy0tDc7Ozg2243FqGEP6k8eocUaNGqX9u1u3boiIiEBQUBA2bdqkHRRuKccnL3VZqXnz5mHnzp1IT0+Hr69vo22VSiX8/f1x7tw5AIC3tzdqampw7do1nXZlZWX1ErlUODk5ITg4GL1790ZCQgK6d++O//7v/9belfDgf3Hc31fsT/0a6lN9eIw27ujRoygrK0N4eDgcHBzg4OCAjIwMfPzxx3BwcND2CY9TwzysP/UNpuUxapy2bduiW7duOHfunMX9e5TBx8oIIRAXF4ft27dj7969CAgIeOg6V69eRVFREZRKJQAgPDwcjo6O2LNnj7aNSqVCbm4u+vfv32K1WxMhBKqrqxEQEABvb2+dvqqpqUFGRoa2r9ifhrnXp/rwGG3c0KFDceLECeTk5Gin3r17Y8qUKcjJyUFgYCCPUyM8rD/t7e3rrcNj1DjV1dU4ffo0lEql5f171KRDpanFzZkzRygUCrFv3z6hUqm0061bt4QQd+9SWLhwofjtt99EQUGBSE9PFxEREaJTp06ioqJCu53Zs2cLX19f8dNPP4ns7GwxZMgQ0b17d3Hnzh1zfTWzWbRokdi/f78oKCgQx48fF2+//baws7MTaWlpQggh3n//faFQKMT27dvFiRMnxF//+lehVCrZn41orE95jJrGg3ch8Thtnvv7k8eo8RYuXCj27dsnLly4IA4dOiTGjh0rXF1dxcWLF4UQlnV8MvhYGQB6p+TkZCGEELdu3RIjRowQHTt2FI6OjqJz585ixowZorCwUGc7t2/fFnFxccLd3V24uLiIsWPH1msjFTExMcLf3184OTmJjh07iqFDh2pDjxB3b8WMj48X3t7eQi6Xi4EDB4oTJ07obIP9qauxPuUxahoPBh8ep81zf3/yGDXexIkThVKpFI6OjsLHx0eMHz9enDx5Urvcko5PmRBCmPYcEhEREZFl4hgfIiIikgwGHyIiIpIMBh8iIiKSDAYfIiIikgwGHyIiIpIMBh8iIiKSDAYfIiIikgwGHyIiIpIMBh8iIiKSDAYfIiIikgwGHyIiIpIMBh8iIiKSDAYfIiIikgwGHyIiIpIMBh8iIiKSDAYfIiIikgwGHyIiIpIMBh89EhIS0KdPH7i6usLT0xPPPfcc8vLyGl1n3759kMlk9aYzZ860UtVERET0MAw+emRkZGDu3Lk4dOgQ9uzZgzt37mDEiBGorKx86Lp5eXlQqVTa6dFHH22FiomIiMgQMiGEMHcRlu7y5cvw9PRERkYGBg4cqLfNvn37EBkZiWvXrqF9+/ZN2o9Go0FJSQlcXV0hk8maUTEREZG0CCFw48YN+Pj4wM6u4fM6Dq1Yk9VSq9UAAHd394e27dmzJ6qqqhAaGorFixcjMjLS4P2UlJTAz8+vyXUSERFJXVFREXx9fRtczuDzEEIILFiwAE8//TTCwsIabKdUKrF+/XqEh4ejuroa//73vzF06FDs27evwbNE1dXVqK6u1tkXcPf/NDc3N9N+ESIiIhtWUVEBPz8/uLq6NtqOl7oeYu7cudi1axcOHDjQaILUZ9y4cZDJZNi5c6fe5UuWLMHSpUvrzVer1Qw+RERERqioqIBCoXjobygHNzdi3rx52LlzJ9LT040OPQDQr18/nDt3rsHlixYtglqt1k5FRUXNKZeIiIgegpe69BBCYN68edixYwf27duHgICAJm3n2LFjUCqVDS6Xy+WQy+VNLfOh6jQCmQXlKLtRBU9XZ/QNcIe9HQdNExGRdDH46DF37lx88cUX+Pbbb+Hq6orS0lIAgEKhgIuLC4C7Z2uKi4uRkpICAFi9ejW6dOmCxx9/HDU1Ndi8eTO2bduGbdu2meU77M5VYWnqKajUVdp5SoUz4seFIiqs4TBGRERkyxh89EhKSgIADB48WGd+cnIyoqOjAQAqlQqFhYXaZTU1NXjttddQXFwMFxcXPP7449i1axdGjx7dWmVr7c5VYc7mbDw4eKtUXYU5m7ORNLUXww8REUkSBzdbEEMHZjWmTiPw9Ad7dc703E8GwFvhjANvDuFlLyIishkc3CxRmQXlDYYeABAAVOoqZBaUt15RREREFoLBx8aU3Wg49DSlHRERkS1h8LExnq7OJm1HRERkSxh8bEzfAHcoFc5oaPSODHfv7uob8PDXbxAREdkaBh8bY28nQ/y4UACoF37ufY4fF8qBzUREJEkMPjYoKkyJpKm94K3QvZzlrXDmrexERCRpfI6PjYoKU2J4qDef3ExERHQfBh8bZm8nQ0TQI+Yug4iIyGLwUhcRERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGgw8RERFJBoMPERERSQaDDxEREUkGg48eCQkJ6NOnD1xdXeHp6YnnnnsOeXl5D10vIyMD4eHhcHZ2RmBgINauXdsK1RIREZGhGHz0yMjIwNy5c3Ho0CHs2bMHd+7cwYgRI1BZWdngOgUFBRg9ejQGDBiAY8eO4e2338bLL7+Mbdu2tWLlRERE1BiZEEKYuwhLd/nyZXh6eiIjIwMDBw7U2+bNN9/Ezp07cfr0ae282bNn43//939x8OBBg/ZTUVEBhUIBtVoNNzc3k9ROREQkBYb+hvKMjwHUajUAwN3dvcE2Bw8exIgRI3TmjRw5EllZWaitrW3R+oiIiMgwDuYuwNIJIbBgwQI8/fTTCAsLa7BdaWkpvLy8dOZ5eXnhzp07uHLlCpRKZb11qqurUV1drf1cUVFhusKJiIioHp7xeYi4uDgcP34cX3755UPbymQync/3riI+OP+ehIQEKBQK7eTn59f8gomIiKhBDD6NmDdvHnbu3In09HT4+vo22tbb2xulpaU688rKyuDg4IBHHnlE7zqLFi2CWq3WTkVFRSarnYiIiOrjpS49hBCYN28eduzYgX379iEgIOCh60RERCA1NVVnXlpaGnr37g1HR0e968jlcsjlcpPUTERERA/HMz56zJ07F5s3b8YXX3wBV1dXlJaWorS0FLdv39a2WbRoEaZPn679PHv2bFy6dAkLFizA6dOn8a9//QsbNmzAa6+9Zo6vQERERHow+OiRlJQEtVqNwYMHQ6lUaqevvvpK20alUqGwsFD7OSAgAN9//z327duHHj164N1338XHH3+MP//5z+b4CkRERKQHn+NjQfgcHyIioqbhc3yIiIiIHsDgQ0RERJLB4ENERESSweBDREREksHgQ0RERJJhUw8wvHjxIn755RdcvHgRt27dQseOHdGzZ09ERETA2dnZ3OURERGRmdlE8Pniiy/w8ccfIzMzE56enujUqRNcXFxQXl6O/Px8ODs7Y8qUKXjzzTfh7+9v7nKJiIjITKw++PTq1Qt2dnaIjo7G//zP/6Bz5846y6urq3Hw4EFs2bIFvXv3xj//+U9MmDDBTNUSERGROVn9Awx37dqFMWPGGNT2ypUrKCgoQJ8+fVq4qqbhAwyJiIiaxtDfUKs/42No6AEADw8PeHh4tGA1REREZMmsPvg05Pbt26itrdWZx7MoRERE0mZTt7PfunULcXFx8PT0RLt27dChQwediYiIiKTNpoLP66+/jr179+Kf//wn5HI5Pv/8cyxduhQ+Pj5ISUkxd3lERERkZjZ1qSs1NRUpKSkYPHgwYmJiMGDAAAQHB8Pf3x//+c9/MGXKFHOXSERERGZkU2d8ysvLERAQAODueJ7y8nIAwNNPP439+/ebszQiIiKyADYVfAIDA3Hx4kUAQGhoKP7nf/4HwN0zQe3btzdfYURERGQRbCr4vPDCC/jf//1fAMCiRYu0Y31effVVvP7662aujoiIiMzN6h9g2JjCwkJkZWUhKCgI3bt3N3c5D8UHGBIRETWNZB5gCAAajQYfffQRvvnmG9TW1mLYsGH4xz/+gc6dO9d7hQURERFJl01c6vrggw/w1ltvoW3btlAqlVi1ahVefvllc5dFREREFsYmgs/GjRuxZs0apKWl4dtvv8U333yDlJQU2PBVPCIiImoCmwg+ly5dwtixY7WfR44cCSEESkpKzFgVERERWRqbGONTU1MDFxcX7WeZTAYnJydUV1ebsSoytTqNQGZBOcpuVMHT1Rl9A9xhbyczd1lERGRFbCL4AMDf//53tGnTRvu5pqYGy5cvh0Kh0M5btWqVOUojE9idq8LS1FNQqau085QKZ8SPC0VUmNKMlRERkTWxidvZBw8eDJms8f/yl8lk2Lt3bytV1DS8nV2/3bkqzNmcjQcP1Hv/jydN7cXwQ0QkcZK6nX3fvn3mLoFaSJ1GYGnqqXqhBwAE7oafpamnMDzUm5e9iIjooWxicDPZrsyCcp3LWw8SAFTqKmQWlLdeUUREZLVs4ozP9evX8eWXX2LOnDkAgClTpuD27dva5fb29vjss8/4vi4rVHaj4dDTlHZERCRtNnHG57PPPsOvv/6q/bxz507Y2dlBoVBAoVDgxIkTWL16tfkKpCbzdHU2aTsiIpI2mwg+X3/9NSZPnqwz78MPP0RycjKSk5ORkJCAb7/91kzVUXP0DXCHUuGMhkbvyHD37q6+Ae6tWRYREVkpmwg++fn5CA4O1n7u2rUrnJyctJ+7d++Oc+fOmaM0aiZ7Oxnix4UCQL3wc+9z/LhQDmwmIiKD2ETwuXXrFmpqarSfs7Ky4Ovrq/1cWVkJjUZjjtLIBKLClEia2gveCt3LWd4KZ97KTkRERrGJ4BMYGIjs7OwGl2dlZSEgIMDg7e3fvx/jxo2Dj48PZDIZvvnmm0bb79u3DzKZrN505swZg/dJjYsKU+LAm0PwZWw//PekHvgyth8OvDmEoYeIiIxiE3d1/b//9/+wePFijBgxAt7e3jrLVCoV4uPjMX36dIO3V1lZie7du+OFF17An//8Z4PXy8vL03loUseOHQ1elx7O3k6GiKBHzF0GERFZMZsIPm+88Qa2bduGkJAQTJs2DSEhIdozLps3b0anTp3w5ptvGry9UaNGYdSoUUbX4enpyVvmiYiILJhNBB9XV1f8+uuvWLRoEb788ktcv34dANC+fXtMnjwZK1asgKura4vX0bNnT1RVVSE0NBSLFy9GZGRko+2rq6t1XqRaUVHR0iUSERFJmk2M8QGADh06YO3atbh69SpKS0tRWlqKq1evYu3atXB3b9lbnZVKJdavX49t27Zh+/bt6Nq1K4YOHYr9+/c3ul5CQoL2WUMKhQJ+fn4tWicREZHU2cRLSluSTCbDjh078Nxzzxm13rhx4yCTybBz584G2+g74+Pn58eXlBIRERnJ0JeUWv0Zn6ioKPz2228PbXfjxg188MEH+PTTT1uhKqBfv34PfXaQXC6Hm5ubzkREREQtx+rH+EyYMAF/+ctf4OrqimeeeQa9e/eGj48PnJ2dce3aNZw6dQoHDhzA999/j7Fjx2LlypWtUtexY8egVPJWayIiIkti9cFn5syZmDZtGr7++mt89dVX+Oyzz7SDm2UyGUJDQzFy5EgcPXoUXbt2NWibN2/exPnz57WfCwoKkJOTA3d3d3Tu3BmLFi1CcXExUlJSAACrV69Gly5d8Pjjj6OmpgabN2/Gtm3bsG3bNpN/XyIiImo6qw8+AODk5ITJkydr39elVqtx+/ZtPPLII3B0dDR6e1lZWTp3ZC1YsAAAMGPGDGzcuBEqlQqFhYXa5TU1NXjttddQXFwMFxcXPP7449i1axdGjx7dzG9GREREpsTBzRbE0IFZREREpEsyg5uJiIiIDMXgQ0RERJLB4ENERESSweBDREREkmFzwef69ev4/PPPsWjRIpSXlwMAsrOzUVxcbObKiIiIyNxs4nb2e44fP45hw4ZBoVDg4sWLiI2Nhbu7O3bs2IFLly5pn7tDRERE0mRTZ3wWLFiA6OhonDt3Ds7Oztr5o0aNeugLQ4mIiMj22VTwOXLkCF566aV68zt16oTS0lIzVERERESWxKaCj7OzMyoqKurNz8vLQ8eOHc1QEREREVkSmwo+zz77LJYtW4ba2loAd9/VVVhYiLfeegt//vOfzVwdERERmZtNBZ/ExERcvnwZnp6euH37NgYNGoTg4GC4urpi+fLl5i6PiIiIzMym7upyc3PDgQMHsHfvXmRnZ0Oj0aBXr14YNmyYuUsjIiIiC2AzwefOnTtwdnZGTk4OhgwZgiFDhpi7JCIiIrIwNnOpy8HBAf7+/qirqzN3KURERGShbCb4AMDixYt1nthMREREdD+budQFAB9//DHOnz8PHx8f+Pv7o23btjrLs7OzzVQZERERWQKbCj7PPfecuUsgIiIiCyYTQghzF0F3VVRUQKFQQK1Ww83NzdzlEBERWQ1Df0NtaowPERERUWNs6lKXnZ0dZDJZg8t5xxcREZG02VTw2bFjh87n2tpaHDt2DJs2bcLSpUvNVBURERFZCkmM8fniiy/w1Vdf4dtvvzV3KY3iGB8iIqKm4Rif+zz55JP46aefzF0GERERmZnNB5/bt29jzZo18PX1NXcpREREZGY2NcanQ4cOOoObhRC4ceMG2rRpg82bN5uxMiIiIrIENhV8/uu//ksn+NjZ2aFjx4548skn0aFDBzNWRkRERJbApoJPdHS0uUsgIiIiC2ZTY3x2796NAwcOaD9/+umn6NGjByZPnoxr166ZsTIiIiKyBDYVfF5//XVUVFQAAE6cOIEFCxZg9OjRuHDhAhYsWGDm6oiIiMjcbOpSV0FBAUJDQwEA27Ztw7hx47BixQpkZ2dj9OjRZq6OiIiIzM2mzvg4OTnh1q1bAICffvoJI0aMAAC4u7trzwQRERGRdNnUGZ+nn34aCxYswFNPPYXMzEx89dVXAICzZ8/yOT5kNeo0ApkF5Si7UQVPV2f0DXCHvV3D76AjIiLD2dQZn08++QQODg74+uuvkZSUhE6dOgEAfvjhB0RFRRm8nf3792PcuHHw8fGBTCbDN99889B1MjIyEB4eDmdnZwQGBmLt2rVN/RokYbtzVXj6g73462eHMH9LDv762SE8/cFe7M5Vmbs0IiKbYFNnfDp37ozvvvuu3vz/+q//Mmo7lZWV6N69O1544QX8+c9/fmj7goICjB49GrGxsdi8eTN+/fVX/O1vf0PHjh0NWp8IuBt65mzOxoMvzytVV2HO5mwkTe2FqDClWWojIrIVNhV8srOz4ejoiG7dugEAvv32WyQnJyM0NBRLliyBk5OTQdsZNWoURo0aZfB+165di86dO2P16tUAgMceewxZWVlITExk8CGD1GkElqaeqhd6AEAAkAFYmnoKw0O9edmLiKgZbOpS10svvYSzZ88CAC5cuIBJkyahTZs22Lp1K954440W2+/Bgwe1A6nvGTlyJLKyslBbW9vgetXV1aioqNCZSJoyC8qhUlc1uFwAUKmrkFlQ3npFERHZIJsKPmfPnkWPHj0AAFu3bsXAgQPxxRdfYOPGjdi2bVuL7be0tBReXl4687y8vHDnzh1cuXKlwfUSEhKgUCi0k5+fX4vVSJat7EbDoacp7YiISD+bCj5CCGg0GgB3b2e/9+wePz+/RgOIKdz/jrB7teibf79FixZBrVZrp6KiohatkSyXp6uzSdsREZF+NjXGp3fv3njvvfcwbNgwZGRkICkpCcDdwccPnpExJW9vb5SWlurMKysrg4ODAx555JEG15PL5ZDL5S1WF1mPvgHuUCqcUaqu0jvORwbAW3H31nYiImo6mzrjs3r1amRnZyMuLg7vvPMOgoODAQBff/01+vfv32L7jYiIwJ49e3TmpaWloXfv3nB0dGyx/ZLtsLeTIX7c3aeOP3iO8N7n+HGhHNhMRNRMMnHvmowNq6qqgr29vcEh5ObNmzh//jwAoGfPnli1ahUiIyPh7u6Ozp07Y9GiRSguLkZKSgqAu2eUwsLC8NJLLyE2NhYHDx7E7Nmz8eWXXxp1V1dFRQUUCgXUajXc3NyM/6Jk9XbnqrA09ZTOQGelwhnx40J5KzsRUSMM/Q21ueBz/fp1fP3118jPz8frr78Od3d3ZGdnw8vLS/tAw4fZt28fIiMj682fMWMGNm7ciOjoaFy8eBH79u3TLsvIyMCrr76KkydPwsfHB2+++SZmz55tVO0MPgTwyc1ERE0hyeBz/PhxDB06FO3bt8fFixeRl5eHwMBA/P3vf8elS5e0Z2gsFYMPERFR0xj6G2pTY3wWLFiAF154AefOnYOz8//d/TJq1Cjs37/fjJURERGRJbCp4HPkyBG89NJL9eZ36tSp3l1XREREJD02FXycnZ31Pv04Ly8PHTt2NENFREREZElsKvg8++yzWLZsmfY1ETKZDIWFhXjrrbf4ziwiIiKyreCTmJiIy5cvw9PTE7dv38agQYMQHBwMV1dXLF++3NzlERERkZnZ1JOb3dzccODAAezduxfZ2dnQaDTo1asXhg0bZu7SiIiIyALYTPC5c+cOnJ2dkZOTgyFDhmDIkCHmLomIiIgsjM1c6nJwcIC/vz/q6urMXQoRERFZKJsJPgCwePFiLFq0COXl5eYuhYiIiCyQzVzqAoCPP/4Y58+fh4+PD/z9/dG2bVud5dnZ2WaqjIiIiCyBTQWfZ599FjIZ32lERERE+tnUu7qsHd/VRURE1DSSelfXrVu3MHfuXHTq1Amenp6YPHkyrly5Yu6yiIiIyMLYRPCJj4/Hxo0bMWbMGEyaNAl79uzBnDlzzF0WERERWRibGOOzfft2bNiwAZMmTQIATJ06FU899RTq6upgb29v5uqIiIjIUtjEGZ+ioiIMGDBA+7lv375wcHBASUmJGasiIiIiS2MTwaeurg5OTk468xwcHHDnzh0zVURERESWyCYudQkhEB0dDblcrp1XVVWF2bNn6zzLZ/v27eYoj4iIiCyETQSfGTNm1Js3depUM1RCRERElswmgk9ycrK5SyAiIiIrYBNjfIiIiIgMweBDREREksHgQ0RERJLB4ENERESSweBDREREkmETd3URkW2p0whkFpSj7EYVPF2d0TfAHfZ2MnOXRUQ2gMGHiCzK7lwVlqaegkpdpZ2nVDgjflwoosKUZqyMiGwBL3URkcXYnavCnM3ZOqEHAErVVZizORu7c1VmqoyIbAWDDxFZhDqNwNLUUxB6lt2btzT1FOo0+loQERmGwYeILEJmQXm9Mz33EwBU6ipkFpS3XlFEZHMYfIjIIpTdaDj0NKUdEZE+DD5EZBE8XZ1N2o6ISB8Gn0b885//REBAAJydnREeHo5ffvmlwbb79u2DTCarN505c6YVKyayXn0D3KFUOKOhm9ZluHt3V98A99Ysi4hsDINPA7766iu88soreOedd3Ds2DEMGDAAo0aNQmFhYaPr5eXlQaVSaadHH320lSomsm72djLEjwsFgHrh597n+HGhfJ4PETULg08DVq1ahZkzZ2LWrFl47LHHsHr1avj5+SEpKanR9Tw9PeHt7a2d7O3tW6liIusXFaZE0tRe8FboXs7yVjgjaWovPseHiJqNDzDUo6amBkePHsVbb72lM3/EiBH47bffGl23Z8+eqKqqQmhoKBYvXozIyMgG21ZXV6O6ulr7uaKionmFE9mAqDAlhod688nNRNQiGHz0uHLlCurq6uDl5aUz38vLC6WlpXrXUSqVWL9+PcLDw1FdXY1///vfGDp0KPbt24eBAwfqXSchIQFLly41ef1E1s7eToaIoEfMXQYR2SAGn0bIZLr/hSmEqDfvnq5du6Jr167azxERESgqKkJiYmKDwWfRokVYsGCB9nNFRQX8/PxMUDkRERHpwzE+enh4eMDe3r7e2Z2ysrJ6Z4Ea069fP5w7d67B5XK5HG5ubjoTERERtRwGHz2cnJwQHh6OPXv26Mzfs2cP+vfvb/B2jh07BqWSgzGJiIgsBS91NWDBggWYNm0aevfujYiICKxfvx6FhYWYPXs2gLuXqYqLi5GSkgIAWL16Nbp06YLHH38cNTU12Lx5M7Zt24Zt27aZ82sQERHRfRh8GjBx4kRcvXoVy5Ytg0qlQlhYGL7//nv4+/sDAFQqlc4zfWpqavDaa6+huLgYLi4uePzxx7Fr1y6MHj3aXF+BiIiIHiATQvBVxxaioqICCoUCarWa432IiIiMYOhvKMf4EBERkWQw+BAREZFkcIwPEZEVqtMIPt2aqAkYfIiIrMzuXBWWpp6CSl2lnadUOCN+XCjfZ0b0ELzURURkRXbnqjBnc7ZO6AGAUnUV5mzOxu5clZkqI7IODD5ERFaiTiOwNPUU9N2Ke2/e0tRTqNPwZl2ihjD4EBFZicyC8npneu4nAKjUVcgsKG+9ooisDIMPEZGVKLvRcOhpSjsiKWLwISKyEp6uziZtRyRFDD5ERFaib4A7lApnNHTTugx37+7qG+DemmURWRUGHyIiK2FvJ0P8uFAAqBd+7n2OHxfK5/kQNYLBh4jIikSFKZE0tRe8FbqXs7wVzkia2ovP8TFSnUbgYP5VfJtTjIP5V3lHnATwAYZERFYmKkyJ4aHefHJzM/FBkNLEt7NbEL6dnYioddx7EOSDP4D3oiPPnlkfvp2diIhIDz4IUtoYfIiISFL4IEjzMve4Ko7xISIiSeGDIM3HEsZV8YwPERFJCh8EaR6W8oJdBh8iIpIUPgiy9VnSuCoGHyIikhQ+CLL1WdK4KgYfIiKSHD4IsnVZ0rgqDm4mIiJJ4oMgW48ljati8CEiIsmyt5MhIugRc5dh8+6NqypVV+kd5yPD3bNtrTGuipe6iIiIqEVZ0rgqBh8iIiJqcZYyroqXuoiIiKhVWMK4KgYfIiIiajXmHlfFS11EREQkGQw+REREJBm81GVBhLh7k19FRYWZKyEiIrIu93477/2WNoTBx4LcuHEDAODn52fmSoiIiKzTjRs3oFAoGlwuEw+LRtRqNBoNSkpK4OrqCpnMNCPcKyoq4Ofnh6KiIri5uZlkm1LG/jQ99qlpsT9Nj31qei3Rp0II3LhxAz4+PrCza3gkD8/4WBA7Ozv4+vq2yLbd3Nz4D6wJsT9Nj31qWuxP02Ofmp6p+7SxMz33cHAzERERSQaDDxEREUkGg4+Nk8vliI+Ph1wuN3cpNoH9aXrsU9Nif5oe+9T0zNmnHNxMREREksEzPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5WJiEhAX369IGrqys8PT3x3HPPIS8vT6dNdHQ0ZDKZztSvXz+dNtXV1Zg3bx48PDzQtm1bPPPMM/j9999b86tYjKSkJDzxxBPaB2lFRETghx9+0C4XQmDJkiXw8fGBi4sLBg8ejJMnT+psg/2p62F9ymO0eRISEiCTyfDKK69o5/E4bTp9/clj1DhLliyp11/e3t7a5ZZ0fDL4WJmMjAzMnTsXhw4dwp49e3Dnzh2MGDEClZWVOu2ioqKgUqm00/fff6+z/JVXXsGOHTuwZcsWHDhwADdv3sTYsWNRV1fXml/HIvj6+uL9999HVlYWsrKyMGTIEDz77LPafyg//PBDrFq1Cp988gmOHDkCb29vDB8+XPtuNYD9+aCH9SnAY7Spjhw5gvXr1+OJJ57Qmc/jtGka6k+Ax6ixHn/8cZ3+OnHihHaZRR2fgqxaWVmZACAyMjK082bMmCGeffbZBte5fv26cHR0FFu2bNHOKy4uFnZ2dmL37t0tWa7V6NChg/j888+FRqMR3t7e4v3339cuq6qqEgqFQqxdu1YIwf401L0+FYLHaFPduHFDPProo2LPnj1i0KBBYv78+UIIweO0iRrqTyF4jBorPj5edO/eXe8ySzs+ecbHyqnVagCAu7u7zvx9+/bB09MTISEhiI2NRVlZmXbZ0aNHUVtbixEjRmjn+fj4ICwsDL/99lvrFG6h6urqsGXLFlRWViIiIgIFBQUoLS3V6Su5XI5BgwZp+4r92bgH+/QeHqPGmzt3LsaMGYNhw4bpzOdx2jQN9ec9PEaNc+7cOfj4+CAgIACTJk3ChQsXAFje8cmXlFoxIQQWLFiAp59+GmFhYdr5o0aNwoQJE+Dv74+CggL8/e9/x5AhQ3D06FHI5XKUlpbCyckJHTp00Nmel5cXSktLW/trWIQTJ04gIiICVVVVaNeuHXbs2IHQ0FDtP3BeXl467b28vHDp0iUAYH82oKE+BXiMNsWWLVuQnZ2NI0eO1Ft2r094nBqusf4EeIwa68knn0RKSgpCQkLwxx9/4L333kP//v1x8uRJizs+GXysWFxcHI4fP44DBw7ozJ84caL277CwMPTu3Rv+/v7YtWsXxo8f3+D2hBCQyWQtVq8l69q1K3JycnD9+nVs27YNM2bMQEZGhnb5g/1iSF9JuT+Bhvs0NDSUx6iRioqKMH/+fKSlpcHZ2bnBdjxODWNIf/IYNc6oUaO0f3fr1g0REREICgrCpk2btIPCLeX45KUuKzVv3jzs3LkT6enp8PX1bbStUqmEv78/zp07BwDw9vZGTU0Nrl27ptOurKysXiKXCicnJwQHB6N3795ISEhA9+7d8d///d/auxIe/C+O+/uK/alfQ32qD4/Rxh09ehRlZWUIDw+Hg4MDHBwckJGRgY8//hgODg7aPuFxapiH9ae+wbQ8Ro3Ttm1bdOvWDefOnbO4f48y+FgZIQTi4uKwfft27N27FwEBAQ9d5+rVqygqKoJSqQQAhIeHw9HREXv27NG2UalUyM3NRf/+/VusdmsihEB1dTUCAgLg7e2t01c1NTXIyMjQ9hX70zD3+lQfHqONGzp0KE6cOIGcnBzt1Lt3b0yZMgU5OTkIDAzkcWqEh/Wnvb19vXV4jBqnuroap0+fhlKptLx/j5p0qDS1uDlz5giFQiH27dsnVCqVdrp165YQ4u5dCgsXLhS//fabKCgoEOnp6SIiIkJ06tRJVFRUaLcze/Zs4evrK3766SeRnZ0thgwZIrp37y7u3Lljrq9mNosWLRL79+8XBQUF4vjx4+Ltt98WdnZ2Ii0tTQghxPvvvy8UCoXYvn27OHHihPjrX/8qlEol+7MRjfUpj1HTePAuJB6nzXN/f/IYNd7ChQvFvn37xIULF8ShQ4fE2LFjhaurq7h48aIQwrKOTwYfKwNA75ScnCyEEOLWrVtixIgRomPHjsLR0VF07txZzJgxQxQWFups5/bt2yIuLk64u7sLFxcXMXbs2HptpCImJkb4+/sLJycn0bFjRzF06FBt6BHi7q2Y8fHxwtvbW8jlcjFw4EBx4sQJnW2wP3U11qc8Rk3jweDD47R57u9PHqPGmzhxolAqlcLR0VH4+PiI8ePHi5MnT2qXW9LxKRNCCNOeQyIiIiKyTBzjQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDRNQMe/fuxZ/+9CdoNBqTbveTTz7BM888Y9JtEhGDDxFRPV26dMHq1asNavvGG2/gnXfegZ3dw/91WlNTAw8PD7z33nt6lyckJMDDwwM1NTWIjY3FkSNHcODAAWNKJ6KHYPAhImqi3377DefOncOECRMMau/k5ISpU6di48aN0Pe2oOTkZEybNg1OTk6Qy+WYPHky1qxZY+qyiSSNwYeIrIoQAh9++CECAwPh4uKC7t274+uvv9Zpc/LkSYwZMwZubm5wdXXFgAEDkJ+fDwAYPHgwXnnlFZ32zz33HKKjo7XLL126hFdffRUymQwymazBWrZs2YIRI0bA2dlZZ35qairCw8Ph7OyMwMBALF26FHfu3AEAzJw5E/n5+di/f7/OOr/88gvOnTuHmTNnauc988wz+Oabb3D79m2j+oiIGsbgQ0RWZfHixUhOTkZSUhJOnjyJV199FVOnTkVGRgYAoLi4GAMHDoSzszP27t2Lo0ePIiYmRhs8Hmb79u3w9fXFsmXLoFKpoFKpGmy7f/9+9O7dW2fejz/+iKlTp+Lll1/GqVOnsG7dOmzcuBHLly8HAHTr1g19+vRBcnKyznr/+te/0LdvX4SFhWnn9e7dG7W1tcjMzDSodiJ6OAdzF0BEZKjKykqsWrUKe/fuRUREBAAgMDAQBw4cwLp16zBo0CB8+umnUCgU2LJlCxwdHQEAISEhBu/D3d0d9vb2cHV1hbe3d6NtL168CB8fH515y5cvx1tvvYUZM2Zo63v33XfxxhtvID4+HgAQExOD1157DZ988gnatWuHmzdvYuvWrVi1apXOttq2bYv27dvj4sWLGDRokMHfgYgaxjM+RGQ1Tp06haqqKgwfPhzt2rXTTikpKdpLWTk5ORgwYIA29LSk27dv17vMdfToUSxbtkynvtjYWKhUKty6dQsA8Ne//hUajQZfffUVAOCrr76CEAKTJk2qtw8XFxftekTUfDzjQ0RW494t47t27UKnTp10lsnlcgB3g0Jj7Ozs6g0srq2tbVI9Hh4euHbtWr0aly5divHjx9drfy8kKRQKPP/880hOTsbMmTORnJyM559/Hm5ubvXWKS8vR8eOHZtUHxHVx+BDRFYjNDQUcrkchYWFDV76eeKJJ7Bp0ybU1tbqPevTsWNHnXE7dXV1yM3NRWRkpHaek5MT6urqHlpPz549cerUKZ15vXr1Ql5eHoKDgxtdd+bMmRg8eDC+++47/Prrr1ixYkW9Nvn5+aiqqkLPnj0fWgsRGYaXuojIari6uuK1117Dq6++ik2bNiE/Px/Hjh3Dp59+ik2bNgEA4uLiUFFRgUmTJiErKwvnzp3Dv//9b+Tl5QEAhgwZgl27dmHXrl04c+YM/va3v+H69es6++nSpQv279+P4uJiXLlypcF6Ro4cWe85O//4xz+QkpKCJUuW4OTJkzh9+jS++uorLF68WKfdoEGDEBwcjOnTpyM4OBgDBw6st/1ffvkFgYGBCAoKakp3EZEeDD5EZFXeffdd/OMf/0BCQgIee+wxjBw5EqmpqQgICAAAPPLII9i7dy9u3ryJQYMGITw8HJ999pn27E9MTAxmzJiB6dOnY9CgQQgICNA52wMAy5Ytw8WLFxEUFNToZaapU6fi1KlT2lAF3A1D3333Hfbs2YM+ffqgX79+WLVqFfz9/eutHxMTg2vXriEmJkbv9r/88kvExsYa3UdE1DCZ0PcULSIiMsgbb7wBtVqNdevWmXS7ubm5GDp0KM6ePQuFQmHSbRNJGc/4EBE1wzvvvAN/f3+DxgQZo6SkBCkpKQw9RCbGMz5EREQkGTzjQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREkvH/AcRpX118GArxAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 600x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figsize(6,6)\n",
    "subplot(211)\n",
    "plot(ecuts, eres[0], 'o')\n",
    "ylabel('Energy (eV)')\n",
    "subplot(212)\n",
    "plot(ecuts, eres[1]/un.GPa, 'o')\n",
    "ylabel('Pressure (GPa)')\n",
    "xlabel('ecut (eV)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd31201-2aca-40da-8078-1d1ddff6abf7",
   "metadata": {},
   "source": [
    "Now we set the cutoff at the selected level for future calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "07dd9450",
   "metadata": {},
   "outputs": [],
   "source": [
    "Si.calc = set_params(Si.calc, ecut = 350)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5667f92e-0d50-4f33-8c6a-bdd955642d83",
   "metadata": {},
   "source": [
    "### Reciprocal space sampling grid\n",
    "\n",
    "The grid over which the momentum space integrations are executed plays the similar role to energy cutoff for momenta. It determines how accurate the summations over reciprocal space will be. For insulators and semiconductors, the sparse grid of just few points usually suffices. Sometimes just a single point at zone-center or off-center is enough. For metals, where Fermi surface enters the Brillouin zone, it needs to be much denser. Since this is a discrete parameter, we provide a list of grid sizes to the scanning function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7dcfd392",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ngkpt=[1, 1, 1]: -211.779 eV ; -6.245 GPa\n",
      "ngkpt=[2, 2, 2]: -214.031 eV ; -0.333 GPa\n",
      "ngkpt=[3, 3, 3]: -214.159 eV ; -0.024 GPa\n",
      "ngkpt=[4, 4, 4]: -214.172 eV ; -0.001 GPa\n",
      "ngkpt=[5, 5, 5]: -214.173 eV ; -0.003 GPa\n"
     ]
    }
   ],
   "source": [
    "kv, kres = scan_param(Si, 'ngkpt', lista=[[1,1,1], [2,2,2], [3,3,3], [4,4,4], [5,5,5]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "02ce084a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAINCAYAAADGPI5YAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABSV0lEQVR4nO3deXQUZd7+/6uz0AFMmiVAAgTComgEWcIAYRCIKKJsjooyKotAUEcEwVGJ4xjwUSJDXEYcEVwCjDqiICpfHhEdSVwAQyQRAQmGRZAQAQOdGEKApH5/+KR/NNm6Q5Luat6vc/ocuuquqs/NfaQvq+6qshiGYQgAAMDk/DxdAAAAQG0g1AAAAJ9AqAEAAD6BUAMAAHwCoQYAAPgEQg0AAPAJhBoAAOATCDUAAMAnBHi6gItFaWmpcnJyFBwcLIvF4ulyAAAwDcMwVFBQoNatW8vPr/LzMYSaepKTk6OIiAhPlwEAgGkdPHhQbdu2rXQ9oaaeBAcHS/p9QEJCQjxcDQAA5pGfn6+IiAjHb2llCDX1pOySU0hICKEGAIAaqG76BhOFAQCAT+BMjUmVlBpK25enIwWn1DI4SH06NJO/HxOQAQAXL0KNCa3bflhz1+zUYfspx7JwW5ASRkZpWNdwD1YGAIDncPnJZNZtP6z73tzqFGgkKdd+Sve9uVXrth/2UGUAAHgWocZESkoNzV2zU0YF68qWzV2zUyWlFbUAAMC3EWpMJG1fXrkzNOcyJB22n1Lavrz6KwoAAC9BqDGRIwWVB5qatAMAwJcQakykZXBQrbYDAMCXEGpMpE+HZgq3BamyG7ct+v0uqD4dmtVnWQAAeAVCjYn4+1mUMDJKksoFm7LvCSOjeF4NAOCiRKgxmWFdw7Xorl4KszlfYgqzBWnRXb14Tg0A4KLFw/dMaFjXcF0XFcYThQEAOAehxqT8/SyK6dTc02UAAOA1uPwEAAB8AqEGAAD4BEINAADwCYQaAADgEwg1AADAJxBqAACATyDUAAAAn2CKULN//35NnjxZHTp0UMOGDdWpUyclJCTo9OnTTu1mzJih6OhoWa1W9ejRo9x+UlJSNHr0aIWHh6tx48bq0aOH3nrrrWqPf/z4cY0bN042m002m03jxo3TiRMnaql3AACgNpji4Xu7du1SaWmpFi9erM6dO2v79u2Ki4tTYWGhkpKSHO0Mw9CkSZP0zTffaNu2beX2s3HjRl111VV69NFH1apVK61du1bjx49XSEiIRo4cWenx77jjDv38889at26dJGnq1KkaN26c1qxZU/udBQAANWIxDMPwdBE1sWDBAi1atEh79+4tt27OnDn64IMPlJmZWe1+hg8frlatWumNN96ocP0PP/ygqKgobd68WX379pUkbd68WTExMdq1a5e6dOniUr35+fmy2Wyy2+0KCQlxaRsAAOD6b6gpLj9VxG63q1mzZnW+n02bNslmszkCjST169dPNptNGzdurHS74uJi5efnO30AAEDdMWWo2bNnjxYuXKh77733gvazcuVKbdmyRXfffXelbXJzc9WyZctyy1u2bKnc3NxKt0tMTHTMwbHZbIqIiLigWgEAQNU8GmrmzJkji8VS5Sc9Pd1pm5ycHA0bNkxjxozRlClTanzslJQUTZw4Ua+++qquvPLKKttaLOXffm0YRoXLy8THx8tutzs+Bw8erHGtAACgeh6dKDxt2jSNHTu2yjaRkZGOP+fk5Cg2NlYxMTFasmRJjY+bmpqqkSNH6rnnntP48eOrbBsWFqZffvml3PKjR4+qVatWlW5ntVpltVprXCMAAHCPR0NNaGioQkNDXWp76NAhxcbGKjo6WsnJyfLzq9lJppSUFI0YMULz58/X1KlTq20fExMju92utLQ09enTR5L0zTffyG63q3///jWqAQAA1D5TzKnJycnR4MGDFRERoaSkJB09elS5ubnl5rRkZ2crMzNTubm5KioqUmZmpjIzMx3Ps0lJSdHw4cM1ffp03XLLLY595OXlOfaRlpamyy+/XIcOHZIkXXHFFRo2bJji4uK0efNmbd68WXFxcRoxYoTLdz4BAIC6Z4pbupcuXVrpZN5zyx88eLBSU1PLtdm3b58iIyM1ceJELVu2rNz6QYMGKSUlRdLvwSc2NtaxjSTl5eVp+vTp+uijjyRJo0aN0ksvvaQmTZq43Adu6QYAoGZc/Q01RajxBYQaAABqxuefUwMAAHAuQg0AAPAJhBoAAOATCDUAAMAnEGoAAIBPINQAAACfQKgBAAA+gVADAAB8AqEGAAD4BEINAADwCYQaAADgEwg1AADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHwCoQYAAPgEQg0AAPAJhBoAAOATCDUAAMAnEGoAAIBPMEWo2b9/vyZPnqwOHTqoYcOG6tSpkxISEnT69GmndjNmzFB0dLSsVqt69OhRbj8pKSkaPXq0wsPD1bhxY/Xo0UNvvfVWtcePjIyUxWJx+syePbu2ugcAAGpBgKcLcMWuXbtUWlqqxYsXq3Pnztq+fbvi4uJUWFiopKQkRzvDMDRp0iR988032rZtW7n9bNy4UVdddZUeffRRtWrVSmvXrtX48eMVEhKikSNHVlnDk08+qbi4OMf3Sy65pPY6CAAALpjFMAzD00XUxIIFC7Ro0SLt3bu33Lo5c+bogw8+UGZmZrX7GT58uFq1aqU33nij0jaRkZF68MEH9eCDD9a43vz8fNlsNtntdoWEhNR4PwAAXGxc/Q1160yNYRhKTU3Vl19+qf379+vkyZNq0aKFevbsqWuvvVYREREXXLir7Ha7mjVrViv7ueKKK6ptN3/+fP3P//yPIiIiNGbMGD388MNq0KBBpe2Li4tVXFzs+J6fn3/BtQIAgMq5NKemqKhI8+bNU0REhG644QatXbtWJ06ckL+/v7Kzs5WQkKAOHTroxhtv1ObNm+u6Zu3Zs0cLFy7Uvffee0H7WblypbZs2aK77767ynYzZszQO++8ow0bNmjatGl64YUX9Je//KXKbRITE2Wz2Ryf+gx8AABclAwXtG3b1rjllluMNWvWGKdPn66wzf79+4158+YZ7dq1M5YsWeLKbo2EhARDUpWfLVu2OG1z6NAho3PnzsbkyZOr3G/37t2rPPaGDRuMxo0bG8uWLXOp1nOtXLnSkGQcO3as0janTp0y7Ha743Pw4EFDkmG3290+HgAAFzO73e7Sb6hLl58+/vhjde3atco27du3V3x8vB566CH99NNPLgWqadOmaezYsVW2iYyMdPw5JydHsbGxiomJ0ZIlS1w6RkVSU1M1cuRIPffccxo/frzb2/fr10+SlJ2drebNm1fYxmq1ymq11rhGAADgHpdCTdeuXZWZmVnhbdLna9CggS699FKXDh4aGqrQ0FCX2h46dEixsbGKjo5WcnKy/Pxqdjd6SkqKRowYofnz52vq1Kk12kdGRoYkKTw8vEbbAwCA2udyMujVq5eio6O1aNEi2e32uqypnJycHA0ePFgRERFKSkrS0aNHlZubq9zcXKd22dnZyszMVG5uroqKipSZmanMzEzH82xSUlI0fPhwTZ8+XbfccotjH3l5eY59pKWl6fLLL9ehQ4ckSZs2bdLzzz+vzMxM7du3T++++67uuecejRo1Su3atau/vwQAAFA1V69nbdy40ZgyZYoREhJiNGzY0LjzzjuNzz///IKvk7kiOTm50jk35xo0aFCFbfbt22cYhmFMmDChwvWDBg1y7GPDhg1O23z77bdG3759DZvNZgQFBRldunQxEhISjMLCQrf64Or1QAAA4MzV31C3n1NTVFSkd999V8nJyfryyy8VGRmpSZMmacKECWrbtm2tBC1fxHNqAACoGVd/Q92emNKwYUNNmDBBKSkp2r17t/785z9r8eLFjlu6AQAAPOGCnyj822+/6a233tJjjz2mEydOqKSkpLZq8ymcqQEAoGbq5InC50pNTdUbb7yhVatWyd/fX7fddpsmT55c090BAABcELdCzcGDB7V06VItXbpU+/btU//+/bVw4ULddtttaty4cV3VCAAAUC2XQ811112nDRs2qEWLFho/frwmTZqkLl261GVtAAAALnM51DRs2FCrVq3SiBEj5O/vX5c1AQAAuM3lUPPRRx85fc/OztaePXs0cOBANWzYUIZhyGKx1HqBAAAArnD7lu5ff/1VQ4YM0WWXXaYbb7xRhw8fliRNmTJFDz30UK0XCAAA4Aq3Q83MmTMVGBioAwcOqFGjRo7lt99+u9atW1erxQEAALjK7Vu6169fr08++aTc04MvvfRSl9/ODQAAUNvcPlNTWFjodIamzLFjx2S1WmulKAAAAHe5HWoGDhyo5cuXO75bLBaVlpZqwYIFio2NrdXiAAAAXOX25acFCxZo8ODBSk9P1+nTp/XII49ox44dysvL09dff10XNQIAAFTL7TM1UVFR2rZtm/r06aPrrrtOhYWFuvnmm5WRkaFOnTrVRY0AAADVuuAXWsI1vNASAICacfU31KUzNQcOHHDr4IcOHXKrPQAAwIVyKdT84Q9/UFxcnNLS0iptY7fb9eqrr6pr1656//33a61AAAAAV7g0UfiHH37QvHnzNGzYMAUGBqp3795q3bq1goKCdPz4ce3cuVM7duxQ7969tWDBAt1www11XTcAAIATt+bUnDp1Sv/7v/+rL7/8Uvv371dRUZFCQ0PVs2dPXX/99eratWtd1mpqzKkBAKBmXP0NZaJwPSHUAABQM7U6URgAAMDbEWoAAIBPINQAAACfQKgBAAA+oUZv6QYAAPA2boeaVq1aadKkSfrqq6/qoh4AAIAacTvU/Oc//5HdbteQIUN02WWX6ZlnnlFOTk5d1AYAAOAyt0PNyJEjtWrVKuXk5Oi+++7Tf/7zH7Vv314jRozQ+++/r7Nnz9ZFnQAAAFWq8UTh5s2ba+bMmfruu+/03HPP6bPPPtOtt96q1q1b64knntDJkydrrcj9+/dr8uTJ6tChgxo2bKhOnTopISFBp0+fdmo3Y8YMRUdHy2q1qkePHuX2k5WVpdjYWLVq1UpBQUHq2LGjHn/8cZ05c6bK4x8/flzjxo2TzWaTzWbTuHHjdOLEiVrrHwAAuHAuvfupIrm5uVq+fLmSk5N14MAB3XrrrZo8ebJycnL0zDPPaPPmzVq/fn2tFLlr1y6VlpZq8eLF6ty5s7Zv3664uDgVFhYqKSnJ0c4wDE2aNEnffPONtm3bVm4/gYGBGj9+vHr16qUmTZrou+++U1xcnEpLSzVv3rxKj3/HHXfo559/1rp16yRJU6dO1bhx47RmzZpa6R8AALhwbr8m4f3331dycrI++eQTRUVFacqUKbrrrrvUpEkTR5sdO3aoZ8+e5c6k1KYFCxZo0aJF2rt3b7l1c+bM0QcffKDMzMxq9zNr1ixt2bJFX375ZYXrf/jhB0VFRWnz5s3q27evJGnz5s2KiYnRrl271KVLF5fq5TUJAADUTJ29JuHuu+9W69at9fXXXyszM1PTpk1zCjSS1LFjR/3tb39zu2h32O12NWvW7IL2kZ2drXXr1mnQoEGVttm0aZNsNpsj0EhSv379ZLPZtHHjxgs6PgAAqD1uX346fPiwGjVqVGWbhg0bKiEhocZFVWfPnj1auHChnn322Rpt379/f23dulXFxcWaOnWqnnzyyUrb5ubmqmXLluWWt2zZUrm5uZVuV1xcrOLiYsf3/Pz8GtUKAABc4/aZmrNnzyo/P7/cp6CgwO3LTXPmzJHFYqnyk56e7rRNTk6Ohg0bpjFjxmjKlCnuli9JWrFihbZu3aq3335ba9eudZqXUxGLxVJumWEYFS4vk5iY6JhYbLPZFBERUaNaAQCAa9w+U9OkSZMqf8zbtm2riRMnKiEhQX5+VWemadOmaezYsVW2iYyMdPw5JydHsbGxiomJ0ZIlS9yq+1xlASMqKkolJSWaOnWqHnroIfn7+5drGxYWpl9++aXc8qNHj6pVq1aVHiM+Pl6zZs1yfM/PzyfYAABQh9wONUuXLtXf/vY3TZw4UX369JFhGNqyZYuWLVumxx9/XEePHlVSUpKsVqsee+yxKvcVGhqq0NBQl4576NAhxcbGKjo6WsnJydUGJlcZhqEzZ86osvnSMTExstvtSktLU58+fSRJ33zzjex2u/r371/pfq1Wq6xWa63UCAAAqud2qFm2bJmeffZZ3XbbbY5lo0aNUrdu3bR48WL997//Vbt27fT0009XG2pclZOTo8GDB6tdu3ZKSkrS0aNHHevCwsIcf87OztZvv/2m3NxcFRUVOe5+ioqKUoMGDfTWW28pMDBQ3bp1k9Vq1bfffqv4+HjdfvvtCgj4/a8iLS1N48eP13//+1+1adNGV1xxhYYNG6a4uDgtXrxY0u+3dI8YMcLlO58AAEDdczvUbNq0Sa+88kq55T179tSmTZskSQMGDNCBAwcuvLr/s379emVnZys7O1tt27Z1WnfuGZYpU6YoNTXVqSZJ2rdvnyIjIxUQEKD58+dr9+7dMgxD7du31/3336+ZM2c6tjl58qSysrKcHsj31ltvafr06Ro6dKik30PcSy+9VGv9AwAAF87t59Rcdtlluvnmm/XMM884LZ89e7ZWr16trKwspaena/To0Tp06FCtFmtmPKcGAICacfU31O0zNUlJSRozZow+/vhj/eEPf5DFYtGWLVu0a9curVy5UpK0ZcsW3X777TWvHgAAwE1un6mRpJ9++kmvvPKKsrKyZBiGLr/8ct1zzz1OdyrBGWdqAAComTo5U3PmzBkNHTpUixcvVmJi4gUXCQAAUFvcui86MDBQ27dvr/I5NQAAAJ7g9sNexo8fr9dff70uagEAAKgxtycKnz59Wq+99po+/fRT9e7dW40bN3Za/9xzz9VacQAAAK5yO9Rs375dvXr1kiTt3r3baR2XpQAAgKe4HWo2bNhQF3UAAABckBq/QCk7O1uffPKJioqKJKnSdycBAADUB7dDza+//qohQ4bosssu04033qjDhw9L+v0VBQ899FCtFwgAAOAKt0PNzJkzFRgYqAMHDqhRo0aO5bfffrvWrVtXq8UBAAC4yu05NevXr9cnn3xS7sWSl156qX766adaKwwAAMAdbp+pKSwsdDpDU+bYsWOyWq21UhQAAIC73A41AwcO1PLlyx3fLRaLSktLtWDBAsXGxtZqcQAAAK5y+/LTggULNHjwYKWnp+v06dN65JFHtGPHDuXl5enrr7+uixoBAACq5faZmqioKG3btk19+vTRddddp8LCQt18883KyMhQp06d6qJGAACAalkMHjBTL1x9bToAAHDm6m+o25efJOnEiRNKS0vTkSNHVFpa6rRu/PjxNdklAADABXE71KxZs0Z33nmnCgsLFRwc7PS+J4vFQqgBAAAe4facmoceekiTJk1SQUGBTpw4oePHjzs+eXl5dVEjAABAtdwONYcOHdL06dMrfFYNAACAp7gdaq6//nqlp6fXRS0AAAA15vacmuHDh+vhhx/Wzp071a1bNwUGBjqtHzVqVK0VBwAA4Cq3b+n286v85I7FYlFJSckFF+WLuKUbAICaqbNbus+/hRsAAMAbuD2nBgAAwBu5HGpuvPFG2e12x/enn35aJ06ccHz/9ddfFRUVVavFAQAAuMrlUPPJJ5+ouLjY8X3+/PlOz6U5e/assrKyarc6AAAAF7kcas6fT8wrowAAgDcxxZya/fv3a/LkyerQoYMaNmyoTp06KSEhQadPn3ZqN2PGDEVHR8tqtapHjx7l9pOVlaXY2Fi1atVKQUFB6tixox5//HGdOXOmyuNHRkbKYrE4fWbPnl2bXQQAABfI5bufyn7Mz19WH3bt2qXS0lItXrxYnTt31vbt2xUXF6fCwkIlJSU52hmGoUmTJumbb77Rtm3byu0nMDBQ48ePV69evdSkSRN99913iouLU2lpqebNm1dlDU8++aTi4uIc3y+55JLa6yAAALhgLocawzA0ceJEWa1WSdKpU6d07733qnHjxpLkNN+mtg0bNkzDhg1zfO/YsaOysrK0aNEip1Dz4osvSpKOHj1aYajp2LGjOnbs6Pjevn17paSk6Msvv6y2huDgYIWFhV1INwAAQB1y+fLThAkT1LJlS9lsNtlsNt11111q3bq143vLli3r9Q3ddrtdzZo1u6B9ZGdna926dRo0aFC1befPn6/mzZurR48eevrpp8td+jpfcXGx8vPznT4AAKDuuHymJjk5uS7rcMuePXu0cOFCPfvsszXavn///tq6dauKi4s1depUPfnkk1W2nzFjhnr16qWmTZsqLS1N8fHx2rdvn1577bVKt0lMTNTcuXNrVB8AAHCfRycKz5kzp9wE3PM/5788MycnR8OGDdOYMWM0ZcqUGh13xYoV2rp1q95++22tXbvW6RJWRWbOnKlBgwbpqquu0pQpU/TKK6/o9ddf16+//lrpNvHx8bLb7Y7PwYMHa1QrAABwjduvSahN06ZN09ixY6tsExkZ6fhzTk6OYmNjFRMToyVLltT4uBEREZKkqKgolZSUaOrUqXrooYfk7+/v0vb9+vWT9Pvlq+bNm1fYxmq1OuYfAQCAuufRUBMaGqrQ0FCX2h46dEixsbGKjo5WcnJylS/WdIdhGDpz5oxbz93JyMiQJIWHh9dKDQAA4MJ5NNS4KicnR4MHD1a7du2UlJSko0ePOtade0dSdna2fvvtN+Xm5qqoqEiZmZmSfj8j06BBA7311lsKDAxUt27dZLVa9e233yo+Pl633367AgJ+/6tIS0vT+PHj9d///ldt2rTRpk2btHnzZsXGxspms2nLli2aOXOmRo0apXbt2tXr3wMAAKicKULN+vXrlZ2drezsbLVt29Zp3blnWKZMmaLU1FTH9549e0qS9u3bp8jISAUEBGj+/PnavXu3DMNQ+/btdf/992vmzJmObU6ePKmsrCzHA/msVqtWrFihuXPnqri4WO3bt1dcXJweeeSRuuwyAABwk8XgfQf1Ij8/XzabTXa7XSEhIZ4uBwAA03D1N9QUr0kAAACoDqEGAAD4BEINAADwCYQaAADgEwg1AADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHwCoQYAAPgEQg0AAPAJhBoAAOATCDUAAMAnEGoAAIBPINQAAACfQKgBAAA+gVADAAB8AqEGAAD4BEINAADwCYQaAADgEwg1AADAJ5gi1Ozfv1+TJ09Whw4d1LBhQ3Xq1EkJCQk6ffq0U7sZM2YoOjpaVqtVPXr0qHKf2dnZCg4OVpMmTao9/vHjxzVu3DjZbDbZbDaNGzdOJ06cqHmHAABArTNFqNm1a5dKS0u1ePFi7dixQ88//7xeeeUVPfbYY07tDMPQpEmTdPvtt1e5vzNnzujPf/6zrr76apeOf8cddygzM1Pr1q3TunXrlJmZqXHjxtW4PwAAoPZZDMMwPF1ETSxYsECLFi3S3r17y62bM2eOPvjgA2VmZla47aOPPqqcnBwNGTJEDz74YJVnXX744QdFRUVp8+bN6tu3ryRp8+bNiomJ0a5du9SlSxeX6s3Pz5fNZpPdbldISIhL28D8SkoNpe3L05GCU2oZHKQ+HZrJ38/i6bIAwFRc/Q0NqMeaapXdblezZs3c3u7zzz/Xe++9p8zMTL3//vvVtt+0aZNsNpsj0EhSv379ZLPZtHHjxkpDTXFxsYqLix3f8/Pz3a4V5rZu+2HNXbNTh+2nHMvCbUFKGBmlYV3DPVgZAPgmU1x+Ot+ePXu0cOFC3XvvvW5t9+uvv2rixIlaunSpy2dLcnNz1bJly3LLW7Zsqdzc3Eq3S0xMdMzBsdlsioiIcKtWmNu67Yd135tbnQKNJOXaT+m+N7dq3fbDHqoMAHyXR0PNnDlzZLFYqvykp6c7bZOTk6Nhw4ZpzJgxmjJlilvHi4uL0x133KGBAwe6tZ3FUv5ygWEYFS4vEx8fL7vd7vgcPHjQrWPCvEpKDc1ds1MVXdctWzZ3zU6VlJryyi8AeC2PXn6aNm2axo4dW2WbyMhIx59zcnIUGxurmJgYLVmyxO3jff755/roo4+UlJQk6fdgUlpaqoCAAC1ZskSTJk0qt01YWJh++eWXcsuPHj2qVq1aVXosq9Uqq9Xqdo0wv7R9eeXO0JzLkHTYfkpp+/IU06l5/RUGAD7Oo6EmNDRUoaGhLrU9dOiQYmNjFR0dreTkZPn5uX+SadOmTSopKXF8//DDDzV//nxt3LhRbdq0qXCbmJgY2e12paWlqU+fPpKkb775Rna7Xf3793e7Bvi+IwWVB5qatAMAuMYUE4VzcnI0ePBgtWvXTklJSTp69KhjXVhYmOPP2dnZ+u2335Sbm6uioiLH3U9RUVFq0KCBrrjiCqf9pqeny8/PT127dnUsS0tL0/jx4/Xf//5Xbdq00RVXXKFhw4YpLi5OixcvliRNnTpVI0aMcPnOJ1xcWgYH1Wo7AIBrTBFq1q9fr+zsbGVnZ6tt27ZO6869I33KlClKTU11fO/Zs6ckad++fU6Xsapy8uRJZWVl6cyZM45lb731lqZPn66hQ4dKkkaNGqWXXnqppt2Bj+vToZnCbUHKtZ+qcF6NRVKY7ffbuwEAtce0z6kxG55Tc3Epu/tJklOwKZtavuiuXtzWDQAucvU31JS3dAPebljXcC26q5fCbM6XmMJsQQQaAKgjprj8BJjRsK7hui4qjCcKezme+gz4DkINUIf8/Szctu3FeOqzeRA+zcHT40SoAXBRKpv3dP6kwrKnPnOZ0HsQPs3BG8aJOTUALjo89dk8eOWIOXjLOBFqAFx03HnqMzyH8GkO3jROhBoAFx2e+mwOhE9z8KZxItQAuOjw1GdzIHyagzeNE6EGwEWn7KnPld2TYdHvExx56rNnET7NwZvGiVAD4KLj72dRwsgoSSoXbMq+J4yM4pZhDyN8moM3jROhBsBFiac+ez/Cpzl40zjx7qd6wrufAO/k6YeFoXre8PwTVK8ux8nV31BCTT0h1ABAzRE+zaGuxsnV31CeKAwA8Hq8csQcPD1OzKkBAAA+gVADAAB8Apef6knZ1KX8/HwPVwIAgLmU/XZWNw2YUFNPCgoKJEkREREergQAAHMqKCiQzWardD13P9WT0tJS5eTkKDg4WBZL7czYz8/PV0REhA4ePOgzd1TRJ+/na/2R6JNZ0CdzqIs+GYahgoICtW7dWn5+lc+c4UxNPfHz81Pbtm3rZN8hISE+8x9DGfrk/XytPxJ9Mgv6ZA613aeqztCUYaIwAADwCYQaAADgEwg1Jma1WpWQkCCr1erpUmoNffJ+vtYfiT6ZBX0yB0/2iYnCAADAJ3CmBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQarzUF198oZEjR6p169ayWCz64IMPqt0mNTVV0dHRCgoKUseOHfXKK6/UfaFucLdPKSkpslgs5T67du2qn4JdkJiYqD/84Q8KDg5Wy5YtddNNNykrK6va7bx1rGrSH28fp0WLFumqq65yPAgsJiZGH3/8cZXbeOv4lHG3T94+RhVJTEyUxWLRgw8+WGU7bx+rMq70xwzjNGfOnHL1hYWFVblNfY4RocZLFRYWqnv37nrppZdcar9v3z7deOONuvrqq5WRkaHHHntM06dP16pVq+q4Ute526cyWVlZOnz4sONz6aWX1lGF7ktNTdX999+vzZs369NPP9XZs2c1dOhQFRYWVrqNN49VTfpTxlvHqW3btnrmmWeUnp6u9PR0XXPNNRo9erR27NhRYXtvHp8y7vapjLeO0fm2bNmiJUuW6KqrrqqynRnGSnK9P2W8fZyuvPJKp/q+//77StvW+xgZ8HqSjNWrV1fZ5pFHHjEuv/xyp2X33HOP0a9fvzqsrOZc6dOGDRsMScbx48frpabacOTIEUOSkZqaWmkbM42VK/0x4zg1bdrUeO211ypcZ6bxOVdVfTLTGBUUFBiXXnqp8emnnxqDBg0yZsyYUWlbM4yVO/0xwzglJCQY3bt3d7l9fY8RZ2p8xKZNmzR06FCnZddff73S09N15swZD1VVO3r27Knw8HANGTJEGzZs8HQ5VbLb7ZKkZs2aVdrGTGPlSn/KmGGcSkpK9M4776iwsFAxMTEVtjHT+Eiu9amMGcbo/vvv1/Dhw3XttddW29YMY+VOf8p4+zj9+OOPat26tTp06KCxY8dq7969lbat7zHihZY+Ijc3V61atXJa1qpVK509e1bHjh1TeHi4hyqrufDwcC1ZskTR0dEqLi7Wv//9bw0ZMkQpKSkaOHCgp8srxzAMzZo1SwMGDFDXrl0rbWeWsXK1P2YYp++//14xMTE6deqULrnkEq1evVpRUVEVtjXL+LjTJzOMkSS988472rp1q7Zs2eJSe28fK3f7Y4Zx6tu3r5YvX67LLrtMv/zyi5566in1799fO3bsUPPmzcu1r+8xItT4EIvF4vTd+L+HRZ+/3Cy6dOmiLl26OL7HxMTo4MGDSkpK8pr/wM81bdo0bdu2TV999VW1bc0wVq72xwzj1KVLF2VmZurEiRNatWqVJkyYoNTU1EpDgBnGx50+mWGMDh48qBkzZmj9+vUKCgpyeTtvHaua9McM43TDDTc4/tytWzfFxMSoU6dOWrZsmWbNmlXhNvU5Rlx+8hFhYWHKzc11WnbkyBEFBARUmJ7Nql+/fvrxxx89XUY5DzzwgD766CNt2LBBbdu2rbKtGcbKnf5UxNvGqUGDBurcubN69+6txMREde/eXf/85z8rbGuG8ZHc61NFvG2Mvv32Wx05ckTR0dEKCAhQQECAUlNT9eKLLyogIEAlJSXltvHmsapJfyribeN0vsaNG6tbt26V1ljfY8SZGh8RExOjNWvWOC1bv369evfurcDAQA9VVfsyMjI8fkr5XIZh6IEHHtDq1auVkpKiDh06VLuNN49VTfpTEW8bp/MZhqHi4uIK13nz+FSlqj5VxNvGaMiQIeXuorn77rt1+eWX69FHH5W/v3+5bbx5rGrSn4p42zidr7i4WD/88IOuvvrqCtfX+xjVyfRjXLCCggIjIyPDyMjIMCQZzz33nJGRkWH89NNPhmEYxuzZs41x48Y52u/du9do1KiRMXPmTGPnzp3G66+/bgQGBhorV670VBfKcbdPzz//vLF69Wpj9+7dxvbt243Zs2cbkoxVq1Z5qgvl3HfffYbNZjNSUlKMw4cPOz4nT550tDHTWNWkP94+TvHx8cYXX3xh7Nu3z9i2bZvx2GOPGX5+fsb69esNwzDX+JRxt0/ePkaVOf9uITOO1bmq648Zxumhhx4yUlJSjL179xqbN282RowYYQQHBxv79+83DMPzY0So8VJlt/ad/5kwYYJhGIYxYcIEY9CgQU7bpKSkGD179jQaNGhgREZGGosWLar/wqvgbp/mz59vdOrUyQgKCjKaNm1qDBgwwFi7dq1niq9ERf2RZCQnJzvamGmsatIfbx+nSZMmGe3btzcaNGhgtGjRwhgyZIjjx98wzDU+Zdztk7ePUWXODwFmHKtzVdcfM4zT7bffboSHhxuBgYFG69atjZtvvtnYsWOHY72nx8hiGP83YwcAAMDEmCgMAAB8AqEGAAD4BEINAADwCYQaAADgEwg1AADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHwCoQYAAPgEQg0AAPAJAZ4u4GJRWlqqnJwcBQcHy2KxeLocAABMwzAMFRQUqHXr1vLzq/x8DKGmnuTk5CgiIsLTZQAAYFoHDx5U27ZtK11PqKknwcHBkn4fkJCQEA9XAwCAeeTn5ysiIsLxW1oZQk09KbvkFBISQqgBAKAGqpu+wURhAADgEzhTAwDweiWlhtL25elIwSm1DA5Snw7N5O/HTRfextPjRKhxw8svv6wFCxbo8OHDuvLKK/XCCy/o6quv9nRZAC6Ap/8RRvXWbT+suWt26rD9lGNZuC1ICSOjNKxruAcrw7m8YZwshmEY9XIkk1uxYoXGjRunl19+WX/84x+1ePFivfbaa9q5c6fatWtX7fb5+fmy2Wyy2+3MqQG8hDf8I4yqrdt+WPe9uVXn/1CVxc5Fd/VirLxAXY+Tq7+hhBoX9e3bV7169dKiRYscy6644grddNNNSkxMrHZ7Qg3gXfix9H4lpYYGzP/cKXSeyyIpzBakrx69hrNrHlQf4+TqbygThV1w+vRpffvttxo6dKjT8qFDh2rjxo0VblNcXKz8/HynDwDvUFJqaO6aneUCjSTHsrlrdqqklP/n86S0fXmV/lBKv4/VYfsppe3Lq7+iUI43jROhxgXHjh1TSUmJWrVq5bS8VatWys3NrXCbxMRE2Ww2x4cH7wHew5v+EUbljhRUPkY1aYe64U3jRKhxw/n3xxuGUek98/Hx8bLb7Y7PwYMH66NEAC7wpn+EUbmWwUG12g51w5vGibufXBAaGip/f/9yZ2WOHDlS7uxNGavVKqvVWh/lAXCTN/0jjMr16dBM4bYg5dpPVXipsGyuRp8Ozeq7NJzDm8aJMzUuaNCggaKjo/Xpp586Lf/000/Vv39/D1UFoKbK/hGubMqiRb/fBcWPpWf5+1mUMDJKksqNVdn3hJFRTBL2MG8aJ0KNi2bNmqXXXntNb7zxhn744QfNnDlTBw4c0L333uvp0gC4yZv+EUbVhnUN16K7einM5nzWLMwWxB1qXsRbxolbut3w8ssv6x//+IcOHz6srl276vnnn9fAgQNd2pZbui9OPNjNu/GcGvPgvyVzqKtx4jk1XoZQc/HhB9Mc+LEEvB+hxssQai4uPNgNAGoPD98DPIQHuwGAZxBqgFrGg90AwDMINUAt48FuAOAZhBqglvFgNwDwDEINUMt4sBsAeAahBqhlPNgNADyDUAPUAW95uiYAXEx4oSVQR4Z1Ddd1UWE82A0A6gmhBqhD/n4WxXRq7ukyAOCiwOUnAADgEwg1AADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQalzw9NNPq3///mrUqJGaNGni6XIAAEAFCDUuOH36tMaMGaP77rvP06UAAIBK8EJLF8ydO1eStHTpUs8WAgAAKkWoqSPFxcUqLi52fM/Pz/dgNQAA+D4uP9WRxMRE2Ww2xyciIsLTJQEA4NMu2lAzZ84cWSyWKj/p6ek13n98fLzsdrvjc/DgwVqsHgAAnO+ivfw0bdo0jR07tso2kZGRNd6/1WqV1Wqt8fYAAMA9F22oCQ0NVWhoqKfLAAAAteSiDTXuOHDggPLy8nTgwAGVlJQoMzNTktS5c2ddcsklni0OAABIItS45IknntCyZcsc33v27ClJ2rBhgwYPHuyhqgAAwLkshmEYni7iYpCfny+bzSa73a6QkBBPlwMAgGm4+ht60d79BAAAfAuhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHyC6V6TsH//fn355Zfav3+/Tp48qRYtWqhnz56KiYlRUFCQp8sDAAAeYppQ8/bbb+vFF19UWlqaWrZsqTZt2qhhw4bKy8vTnj17FBQUpDvvvFOPPvqo2rdv7+lyAQBAPTNFqOnVq5f8/Pw0ceJEvfvuu2rXrp3T+uLiYm3atEnvvPOOevfurZdfflljxozxULUAAMATTPFCy7Vr12r48OEutT127Jj27dunP/zhD3VclXt4oSUAADXj6m+oKc7UuBpoJCk0NFShoaF1WA0AAPBGpgg1lSkqKtKZM2eclnEWBACAi5Ppbuk+efKkpk2bppYtW+qSSy5R06ZNnT4AAODiZLpQ8/DDD+vzzz/Xyy+/LKvVqtdee01z585V69attXz5ck+XBwAAPMR0l5/WrFmj5cuXa/DgwZo0aZKuvvpqde7cWe3bt9dbb72lO++809MlAgAADzDdmZq8vDx16NBB0u/zZ/Ly8iRJAwYM0BdffOHJ0gAAgAeZLtR07NhR+/fvlyRFRUXp3XfflfT7GZwmTZrU+vH279+vyZMnq0OHDmrYsKE6deqkhIQEnT59utaPBQAAas50l5/uvvtufffddxo0aJDi4+M1fPhwLVy4UGfPntVzzz1X68fbtWuXSktLtXjxYnXu3Fnbt29XXFycCgsLlZSUVOvHAwAANWOKh+9V5cCBA0pPT1enTp3UvXv3ejnmggULtGjRIu3du9flbXj4HgAANeNTD9+TpNLSUj377LP64IMPdObMGV177bV64okn1K5du3KvTahrdrtdzZo1q7JNcXGxiouLHd/z8/PruiwAAC5qpplTM3/+fM2ePVuNGzdWeHi4nnvuOU2fPr3e69izZ48WLlyoe++9t8p2iYmJstlsjk9EREQ9VQgAwMXJNJefunTpohkzZugvf/mLJGndunW66aabVFRUJIvF4vb+5syZo7lz51bZZsuWLerdu7fje05OjgYNGqRBgwbptddeq3Lbis7UREREcPkJAAA3uXr5yTShJigoSLt373ZcajIMQ0FBQdq7d6/atGnj9v6OHTumY8eOVdkmMjJSQUFBkn4PNLGxserbt6+WLl0qPz/3TnIxpwYAgJrxuTk1p0+fVsOGDR3fLRaLGjRo4HQ2xB3uvPjy0KFDio2NVXR0tJKTk90ONAAAoO6ZJtRI0t///nc1atTI8f306dN6+umnZbPZHMtq+7bunJwcDR48WO3atVNSUpKOHj3qWBcWFlarxwIAADVnmlAzcOBAZWVlOS3r37+/023VNZlbU53169crOztb2dnZatu2rdM6k1y5AwDgomCaOTVmx5waAABqxtXfUCaHAAAAn2Cay08nTpzQf/7zH913332SpDvvvFNFRUWO9f7+/nr11Vfr5P1PAADA+5nmTM2rr76qr7/+2vH9o48+kp+fn+Phdt9//71eeOEFzxUIAAA8yjShZuXKlbrjjjuclv3jH/9QcnKykpOTlZiYqA8//NBD1QEAAE8zTajZs2ePOnfu7PjepUsXNWjQwPG9e/fu+vHHHz1RGgAA8AKmmVNz8uRJnT592vE9PT3daX1hYaFKS0vruywAAOAlTHOmpmPHjtq6dWul69PT09WhQ4d6rAgAAHgT04SaP/3pT3r88ceVm5tbbt3hw4eVkJCgP/3pTx6oDAAAeAPTPHyvoKBAffv21c8//6xx48bpsssuk8Vi0a5du/Tmm2+qTZs2SktLU3BwsKdLrRAP3wMAoGZ87oWWwcHB+vrrrxUfH6///Oc/OnHihCSpSZMmuuOOOzRv3jyvDTQAAKDumeZMzbkMw3C8WLJFixZ18s6n2saZGgAAasbnztScy2KxqGXLlp4uAwAAeBFTTBQeNmyYNm7cWG27goICzZ8/X//617/qoSoAAOBNTHGmZsyYMbrtttsUHBysUaNGqXfv3mrdurWCgoJ0/Phx7dy5U1999ZX+93//VyNGjNCCBQs8XTIAAKhnpplTc/r0aa1cuVIrVqzQl19+6ZgobLFYFBUVpeuvv15xcXHq0qWLZwutBHNqAACoGVd/Q00Tas5nt9tVVFSk5s2bKzAw0NPlVItQAwBAzfj0RGFJjrdzAwAASCaZKAwAAFAdQg0AAPAJhBoAAOATCDUuGDVqlNq1a6egoCCFh4dr3LhxysnJ8XRZAADgHKYMNSdOnNBrr72m+Ph45eXlSZK2bt2qQ4cO1cnxYmNj9e677yorK0urVq3Snj17dOutt9bJsQAAQM2Y7pbubdu26dprr5XNZtP+/fuVlZWljh076u9//7t++uknLV++vM5r+Oijj3TTTTepuLjY5dvJuaUbAICacfU31HRnambNmqWJEyfqxx9/VFBQkGP5DTfcoC+++KLOj5+Xl6e33npL/fv3rzLQFBcXKz8/3+kDAADqjulCzZYtW3TPPfeUW96mTRvl5ubW2XEfffRRNW7cWM2bN9eBAwf04YcfVtk+MTHR8Swdm82miIiIOqsNAACYMNQEBQVVeNYjKytLLVq0cHk/c+bMkcViqfKTnp7uaP/www8rIyND69evl7+/v8aPH6+qrtzFx8fLbrc7PgcPHnSvowAAwC2mm1MzdepUHT16VO+++66aNWumbdu2yd/fXzfddJMGDhyoF154waX9HDt2TMeOHauyTWRkpNMlrjI///yzIiIitHHjRsXExLh0PObUAABQMz77moSkpCTdeOONatmypYqKijRo0CDl5uYqJiZGTz/9tMv7CQ0NVWhoaI1qKMuBxcXFNdoeAADUPtOFmpCQEH311Vf6/PPPtXXrVpWWlqpXr1669tpr6+R4aWlpSktL04ABA9S0aVPt3btXTzzxhDp16uTyWRoAAFD3TBVqzp49q6CgIGVmZuqaa67RNddcU+fHbNiwod5//30lJCSosLBQ4eHhGjZsmN555x1ZrdY6Pz4AAHCNqUJNQECA2rdvr5KSkno7Zrdu3fT555/X2/EAAEDNmO7up8cff9zpScIAAACSyc7USNKLL76o7OxstW7dWu3bt1fjxo2d1m/dutVDlQEAAE8yXai56aabPF0CAADwQqZ7To1Z8ZwaAABqxmff/QQAAFAR011+8vPzk8ViqXR9fd4ZBQAAvIfpQs3q1audvp85c0YZGRlatmyZ5s6d66GqAACAp/nMnJq3335bK1asqPbt2Z7CnBoAAGrmoptT07dvX3322WeeLgMAAHiIT4SaoqIiLVy4UG3btvV0KQAAwENMN6emadOmThOFDcNQQUGBGjVqpDfffNODlQEAAE8yXah5/vnnnUKNn5+fWrRoob59+6pp06YerAwAAHiS6ULNxIkTPV0CAADwQqabU7Nu3Tp99dVXju//+te/1KNHD91xxx06fvy4BysDAACeZLpQ8/DDDys/P1+S9P3332vWrFm68cYbtXfvXs2aNcvD1QEAAE8x3eWnffv2KSoqSpK0atUqjRw5UvPmzdPWrVt14403erg6AADgKaY7U9OgQQOdPHlSkvTZZ59p6NChkqRmzZo5zuAAAICLj+nO1AwYMECzZs3SH//4R6WlpWnFihWSpN27d/OcGgAALmKmO1Pz0ksvKSAgQCtXrtSiRYvUpk0bSdLHH3+sYcOGebg6AADgKT7z7qf6UFxcrL59++q7775TRkaGevTo4fK2vPsJAICa8dl3P23dulXff/+94/uHH36om266SY899phOnz5dp8d+5JFH1Lp16zo9BgAAqBnThZp77rlHu3fvliTt3btXY8eOVaNGjfTee+/pkUceqbPjfvzxx1q/fr2SkpLq7BgAAKDmTBdqdu/e7bjs895772ngwIF6++23tXTpUq1atapOjvnLL78oLi5O//73v9WoUaM6OQYAALgwprv7yTAMlZaWSvr9lu4RI0ZIkiIiInTs2LE6Od7EiRN17733qnfv3tq/f79L2xUXF6u4uNjxndvNAQCoW6Y7U9O7d2899dRT+ve//63U1FQNHz5c0u8P5WvVqpXL+5kzZ44sFkuVn/T0dC1cuFD5+fmKj493q87ExETZbDbHJyIiwq3tAQCAe0x399O2bdt055136sCBA5o1a5YSEhIkSQ888IB+/fVXvf322y7t59ixY9We2YmMjNTYsWO1Zs0apzeDl5SUyN/fX3feeaeWLVtW4bYVnamJiIjg7icAANzk6t1Ppgs1lTl16pT8/f0VGBhYq/s9cOCA06WjnJwcXX/99Vq5cqX69u3r8gP/uKUbAICacfU31HRzaiTpxIkTWrlypfbs2aOHH35YzZo1086dO9WqVSvHw/hqS7t27Zy+X3LJJZKkTp068QRjAAC8iOlCzbZt2zRkyBA1adJE+/fvV1xcnJo1a6bVq1frp59+0vLlyz1dIgAA8ADTTRSeNWuW7r77bv34448KCgpyLL/hhhv0xRdf1PnxIyMjZRiGW08TBgAAdc90oWbLli265557yi1v06aNcnNzPVARAADwBqYLNUFBQRU+8yUrK0stWrTwQEUAAMAbmC7UjB49Wk8++aTOnDkjSbJYLDpw4IBmz56tW265xcPVAQAATzFdqElKStLRo0fVsmVLFRUVadCgQercubOCg4P19NNPe7o8AADgIaa7+ykkJERfffWVPv/8c23dulWlpaXq1auXrr32Wk+XBgAAPMhUoebs2bMKCgpSZmamrrnmGl1zzTWeLgkAAHgJU11+CggIUPv27VVSUuLpUgAAgJcxVaiRpMcff1zx8fHKy8vzdCkAAMCLmOrykyS9+OKLys7OVuvWrdW+fXs1btzYaf3WrVs9VBkAAPAk04Wa0aNHO70xGwAAQPKht3R7O97SDQBAzbj6G2qaOTUnT57U/fffrzZt2qhly5a64447dOzYMU+XBQAAvIRpQk1CQoKWLl2q4cOHa+zYsfr000913333ebosAADgJUwzp+b999/X66+/rrFjx0qS7rrrLv3xj39USUmJ/P39PVwdAADwNNOcqTl48KCuvvpqx/c+ffooICBAOTk5HqwKAAB4C9OEmpKSEjVo0MBpWUBAgM6ePeuhigAAgDcxzeUnwzA0ceJEWa1Wx7JTp07p3nvvdXpWzfvvv++J8gAAgIeZJtRMmDCh3LK77rrLA5UAAABvZJpQk5yc7OkSAACAFzPNnBoAAICqEGpcEBkZKYvF4vSZPXu2p8sCAADnMM3lJ0978sknFRcX5/h+ySWXeLAaAABwPkKNi4KDgxUWFubpMgAAQCW4/OSi+fPnq3nz5urRo4eefvppnT59usr2xcXFys/Pd/oAAIC6w5kaF8yYMUO9evVS06ZNlZaWpvj4eO3bt0+vvfZapdskJiZq7ty59VglAAAXN4thGIani/CEOXPmVBs6tmzZot69e5dbvmrVKt166606duyYmjdvXuG2xcXFKi4udnzPz89XREREta9NBwAAzvLz82Wz2ar9Db1oz9RMmzbN8XLMykRGRla4vF+/fpKk7OzsSkON1Wp1evoxAACoWxdtqAkNDVVoaGiNts3IyJAkhYeH12ZJAADgAly0ocZVmzZt0ubNmxUbGyubzaYtW7Zo5syZGjVqlNq1a+fp8gAAwP8h1FTDarVqxYoVmjt3roqLi9W+fXvFxcXpkUce8XRpAADgHISaavTq1UubN2/2dBkAAKAaPKcGAAD4BEINAADwCYQaAADgEwg1AADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHwCoQYAAPgEQg0AAPAJhBoAAOATCDUuWrt2rfr27auGDRsqNDRUN998s6dLAgAA5wjwdAFmsGrVKsXFxWnevHm65pprZBiGvv/+e0+XBQAAzkGoqcbZs2c1Y8YMLViwQJMnT3Ys79KliwerAgAA5+PyUzW2bt2qQ4cOyc/PTz179lR4eLhuuOEG7dixw9OlAQCAcxBqqrF3715J0pw5c/T444/r//2//6emTZtq0KBBysvLq3S74uJi5efnO30AAEDduWhDzZw5c2SxWKr8pKenq7S0VJL0t7/9Tbfccouio6OVnJwsi8Wi9957r9L9JyYmymazOT4RERH11TUAAC5KF+2cmmnTpmns2LFVtomMjFRBQYEkKSoqyrHcarWqY8eOOnDgQKXbxsfHa9asWY7v+fn5BBsAAOrQRRtqQkNDFRoaWm276OhoWa1WZWVlacCAAZKkM2fOaP/+/Wrfvn2l21mtVlmt1lqr93wlpYbS9uXpSMEptQwOUp8OzeTvZ6mz4wEA4O0u2lDjqpCQEN17771KSEhQRESE2rdvrwULFkiSxowZ45Ga1m0/rLlrduqw/ZRjWbgtSAkjozSsa7hHagIAwNMINS5YsGCBAgICNG7cOBUVFalv3776/PPP1bRp03qvZd32w7rvza0yzlueaz+l+97cqkV39SLYAAAuShbDMM7/fUQdyM/Pl81mk91uV0hISI32UVJqaMD8z53O0JzLIinMFqSvHr2GS1EAAJ/h6m/oRXv3kxml7curNNBIkiHpsP2U0vZVfqs5AAC+ilBjIkcKKg80NWkHAIAvIdSYSMvgoFptBwCALyHUmEifDs0UbgtSZbNlLPr9Lqg+HZrVZ1kAAHgFQo2J+PtZlDDy94cAnh9syr4njIxikjAA4KJEqDGZYV3DteiuXgqzOV9iCrMFcTs3AOCixnNqTGhY13BdFxXGE4UBADgHocak/P0siunU3NNlAADgNbj8BAAAfAKhBgAA+AQuP9WTsrdR5Ofne7gSAADMpey3s7o3OxFq6klBQYEkKSIiwsOVAABgTgUFBbLZbJWu54WW9aS0tFQ5OTkKDg6WxVI7dynl5+crIiJCBw8erPFLMr0NffJ+vtYfiT6ZBX0yh7rok2EYKigoUOvWreXnV/nMGc7U1BM/Pz+1bdu2TvYdEhLiM/8xlKFP3s/X+iPRJ7OgT+ZQ232q6gxNGSYKAwAAn0CoAQAAPoFQY2JWq1UJCQmyWq2eLqXW0Cfv52v9keiTWdAnc/Bkn5goDAAAfAJnagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHwCocZLffHFFxo5cqRat24ti8WiDz74oNptUlNTFR0draCgIHXs2FGvvPJK3RfqBnf7lJKSIovFUu6za9eu+inYBYmJifrDH/6g4OBgtWzZUjfddJOysrKq3c5bx6om/fH2cVq0aJGuuuoqx4PAYmJi9PHHH1e5jbeOTxl3++TtY1SRxMREWSwWPfjgg1W28/axKuNKf8wwTnPmzClXX1hYWJXb1OcYEWq8VGFhobp3766XXnrJpfb79u3TjTfeqKuvvloZGRl67LHHNH36dK1ataqOK3Wdu30qk5WVpcOHDzs+l156aR1V6L7U1FTdf//92rx5sz799FOdPXtWQ4cOVWFhYaXbePNY1aQ/Zbx1nNq2batnnnlG6enpSk9P1zXXXKPRo0drx44dFbb35vEp426fynjrGJ1vy5YtWrJkia666qoq25lhrCTX+1PG28fpyiuvdKrv+++/r7RtvY+RAa8nyVi9enWVbR555BHj8ssvd1p2zz33GP369avDymrOlT5t2LDBkGQcP368XmqqDUeOHDEkGampqZW2MdNYudIfM45T06ZNjddee63CdWYan3NV1SczjVFBQYFx6aWXGp9++qkxaNAgY8aMGZW2NcNYudMfM4xTQkKC0b17d5fb1/cYcabGR2zatElDhw51Wnb99dcrPT1dZ86c8VBVtaNnz54KDw/XkCFDtGHDBk+XUyW73S5JatasWaVtzDRWrvSnjBnGqaSkRO+8844KCwsVExNTYRszjY/kWp/KmGGM7r//fg0fPlzXXntttW3NMFbu9KeMt4/Tjz/+qNatW6tDhw4aO3as9u7dW2nb+h4jXmjpI3Jzc9WqVSunZa1atdLZs2d17NgxhYeHe6iymgsPD9eSJUsUHR2t4uJi/fvf/9aQIUOUkpKigQMHerq8cgzD0KxZszRgwAB17dq10nZmGStX+2OGcfr+++8VExOjU6dO6ZJLLtHq1asVFRVVYVuzjI87fTLDGEnSO++8o61bt2rLli0utff2sXK3P2YYp759+2r58uW67LLL9Msvv+ipp55S//79tWPHDjVv3rxc+/oeI0KND7FYLE7fjf97WPT5y82iS5cu6tKli+N7TEyMDh48qKSkJK/5D/xc06ZN07Zt2/TVV19V29YMY+Vqf8wwTl26dFFmZqZOnDihVatWacKECUpNTa00BJhhfNzpkxnG6ODBg5oxY4bWr1+voKAgl7fz1rGqSX/MME433HCD48/dunVTTEyMOnXqpGXLlmnWrFkVblOfY8TlJx8RFham3Nxcp2VHjhxRQEBAhenZrPr166cff/zR02WU88ADD+ijjz7Shg0b1LZt2yrbmmGs3OlPRbxtnBo0aKDOnTurd+/eSkxMVPfu3fXPf/6zwrZmGB/JvT5VxNvG6Ntvv9WRI0cUHR2tgIAABQQEKDU1VS+++KICAgJUUlJSbhtvHqua9Kci3jZO52vcuLG6detWaY31PUacqfERMTExWrNmjdOy9evXq3fv3goMDPRQVbUvIyPD46eUz2UYhh544AGtXr1aKSkp6tChQ7XbePNY1aQ/FfG2cTqfYRgqLi6ucJ03j09VqupTRbxtjIYMGVLuLpq7775bl19+uR599FH5+/uX28abx6om/amIt43T+YqLi/XDDz/o6quvrnB9vY9RnUw/xgUrKCgwMjIyjIyMDEOS8dxzzxkZGRnGTz/9ZBiGYcyePdsYN26co/3evXuNRo0aGTNnzjR27txpvP7660ZgYKCxcuVKT3WhHHf79PzzzxurV682du/ebWzfvt2YPXu2IclYtWqVp7pQzn333WfYbDYjJSXFOHz4sONz8uRJRxszjVVN+uPt4xQfH2988cUXxr59+4xt27YZjz32mOHn52esX7/eMAxzjU8Zd/vk7WNUmfPvFjLjWJ2ruv6YYZweeughIyUlxdi7d6+xefNmY8SIEUZwcLCxf/9+wzA8P0aEGi9Vdmvf+Z8JEyYYhmEYEyZMMAYNGuS0TUpKitGzZ0+jQYMGRmRkpLFo0aL6L7wK7vZp/vz5RqdOnYygoCCjadOmxoABA4y1a9d6pvhKVNQfSUZycrKjjZnGqib98fZxmjRpktG+fXujQYMGRosWLYwhQ4Y4fvwNw1zjU8bdPnn7GFXm/BBgxrE6V3X9McM43X777UZ4eLgRGBhotG7d2rj55puNHTt2ONZ7eowshvF/M3YAAABMjInCAADAJxBqAACATyDUAAAAn0CoAQAAPoFQAwAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQDqREpKiiwWi06cOFFpm6VLl6pJkyY1PobFYtEHH3xQ4+1dkZWVpbCwMBUUFNTpcerLrbfequeee87TZQB1glADoFK5ubmaMWOGOnfurKCgILVq1UoDBgzQK6+8opMnT1a5bf/+/XX48GHZbLY6q+/w4cO64YYb6mz/kvS3v/1N999/v4KDgyX9/2Gta9eu5d603KRJEy1dutTlfU+cOFEWi0XPPPOM0/IPPvhAFovFrTotFouCgoL0008/OS2/6aabNHHiRMf3J554Qk8//bTy8/Pd2j9gBoQaABXau3evevbsqfXr12vevHnKyMjQZ599ppkzZ2rNmjX67LPPKt32zJkzatCggcLCwtz+cXZHWFiYrFZrne3/559/1kcffaS777673Lo9e/Zo+fLlF3yMoKAgzZ8/X8ePH7/gfVksFj3xxBNVtrnqqqsUGRmpt95664KPB3gbQg2ACv3lL39RQECA0tPTddttt+mKK65Qt27ddMstt2jt2rUaOXKko63FYtErr7yi0aNHq3HjxnrqqacqvPy0dOlStWvXTo0aNdKf/vQn/frrr1XWcPr0aU2bNk3h4eEKCgpSZGSkEhMTnY5bdvlpzpw5slgs5T5lZ04Mw9A//vEPdezYUQ0bNlT37t21cuXKKo//7rvvqnv37mrbtm25dQ888IASEhJ06tSpCrdNSUlRgwYN9OWXXzqWPfvsswoNDdXhw4cdy6699lqFhYU59et8kyZN0lVXXaXi4mJJv4fG6Oho3XnnneVqevPNN/X9999X2a9Ro0bpP//5T5VtADMi1AAo59dff9X69et1//33q3HjxhW2Of8MTEJCgkaPHq3vv/9ekyZNKtf+m2++0aRJk/SXv/xFmZmZio2N1VNPPVVlHS+++KI++ugjvfvuu8rKytKbb76pyMjICtv+9a9/1eHDhx2fpKQkNWrUSL1795YkPf7440pOTtaiRYu0Y8cOzZw5U3fddZdSU1MrPf4XX3zh2P58Dz74oM6ePauXXnqpwvWDBw/Wgw8+qHHjxslut+u7777T3/72N7366qsKDw93tPP399e8efO0cOFC/fzzz5X+PRQWFmr27NmSpL///e86duyYXn75Zad2/fv314gRIxQfH19pnySpT58+SktLc4QkwGfU2fu/AZjW5s2bDUnG+++/77S8efPmRuPGjY3GjRsbjzzyiGO5JOPBBx90arthwwZDknH8+HHDMAzjz3/+szFs2DCnNrfffrths9kqreOBBx4wrrnmGqO0tLTC9ZKM1atXl1u+adMmIygoyFixYoVhGIbx22+/GUFBQcbGjRud2k2ePNn485//XOnxu3fvbjz55JOV9uuVV14xmjVrZpw4ccIwDMOw2WxGcnKyo21xcbHRs2dP47bbbjOuvPJKY8qUKU77mjBhgjF69GjDMAyjX79+xqRJkwzDMIzVq1cb5//zvHHjRiMwMND4+9//bgQEBBipqakV/l3s2LHD8Pf3N7744gvDMAxj9OjRxoQJE5zafvfdd4YkY//+/ZX2HTAjztQAqNT5Z2PS0tKUmZmpK6+8stz/5Vd2RqPMDz/8oJiYGKdl538/38SJE5WZmakuXbpo+vTpWr9+fbU1HzhwQDfddJP++te/6rbbbpMk7dy5U6dOndJ1112nSy65xPFZvny59uzZU+m+ioqKFBQUVOn6yZMnKzQ0VPPnz69wfYMGDfTmm29q1apVKioq0gsvvFDpvubPn69ly5Zp586dFa6PiYnRX//6V/3P//yPHnroIQ0cOLDCdlFRURo/frweffTRSo/VsGFDSap2sjdgNgGeLgCA9+ncubMsFot27drltLxjx46S/v8fxXNVdpmqjGEYbtfRq1cv7du3Tx9//LE+++wz3Xbbbbr22msrnQtTWFioUaNGKSYmRk8++aRjeWlpqSRp7dq1atOmjdM2VU00Dg0NrXICb0BAgJ566ilNnDhR06ZNq7DNxo0bJUl5eXnKy8ur9O9p4MCBuv766/XYY4853a10bh++/vpr+fv768cff6y0JkmaO3euLrvsskpvd8/Ly5MktWjRosr9AGbDmRoA5TRv3lzXXXedXnrpJRUWFtbKPqOiorR582anZed/r0hISIhuv/12vfrqq1qxYoVWrVrl+FE+l2EYuuuuu1RaWqp///vfTmeZoqKiZLVadeDAAXXu3NnpExERUemxe/bsWemZkzJjxozRlVdeqblz55Zbt2fPHs2cOVOvvvqq+vXrp/HjxzsCVkWeeeYZrVmzxhGEzrVgwQL98MMPSk1N1SeffKLk5ORK9xMREaFp06bpscceK3fbuSRt375dbdu2VWhoaJV9A8yGUAOgQi+//LLOnj2r3r17a8WKFfrhhx8ck3V37dolf39/t/Y3ffp0rVu3Tv/4xz+0e/duvfTSS1q3bl2V2zz//PN65513tGvXLu3evVvvvfeewsLCKnxg35w5c/TZZ59p8eLF+u2335Sbm6vc3FwVFRUpODhYf/3rXzVz5kwtW7ZMe/bsUUZGhv71r39p2bJllR7/+uuv16ZNmyoMBud65pln9MYbbzgFwJKSEo0bN05Dhw7V3XffreTkZG3fvl3PPvtspfvp1q2b7rzzTi1cuNBpeWZmpp544gm9/vrr+uMf/6h//vOfmjFjhvbu3VvpvuLj45WTk1Phrfdffvmlhg4dWmWfAFPy9KQeAN4rJyfHmDZtmtGhQwcjMDDQuOSSS4w+ffoYCxYsMAoLCx3tVMGE3fMnChuGYbz++utG27ZtjYYNGxojR440kpKSqpwovGTJEqNHjx5G48aNjZCQEGPIkCHG1q1bKzzuoEGDDEnlPmUTd0tLS41//vOfRpcuXYzAwECjRYsWxvXXX19uwu25zp49a7Rp08ZYt25dlf0yDMMYOnSo0/Hmzp1rhIeHG8eOHXO0+eCDD4wGDRoYGRkZhmE4TxQus3//fsNqtTomChcVFRlRUVHG1KlTndr96U9/Mvr372+cPXu23N9FmXnz5hmSnCYKFxUVGSEhIcamTZsq7TdgVhbDqMGFbgC4SLz88sv68MMP9cknn3i6lFrxr3/9Sx9++KFLk64Bs2GiMABUYerUqTp+/LgKCgocr0ows8DAwHKXtwBfwZkaAADgE5goDAAAfAKhBgAA+ARCDQAA8AmEGgAA4BMINQAAwCcQagAAgE8g1AAAAJ9AqAEAAD6BUAMAAHzC/we9/705E6fHUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "subplot(211)\n",
    "plot([k[0] for k in kv], kres[0], 'o')\n",
    "ylabel('Energy (eV)')\n",
    "subplot(212)\n",
    "plot([k[0] for k in kv], kres[1]/un.GPa, 'o')\n",
    "ylabel('Pressure (GPa)')\n",
    "xlabel('Grid size (NxNxN)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e14770c9-24c3-443a-9928-64c5f266ad75",
   "metadata": {},
   "source": [
    "Now we fix both parameters for future calculations. In real production calculations, the energy convergence should be checked at this point, to make sure that the resulting calculation is indeed well converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b6e471b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "ngp=4\n",
    "Si.calc = set_params(Si.calc, ecut=350, ngkpt=[ngp, ngp, ngp])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "698b459b",
   "metadata": {},
   "source": [
    "### Cell size optimization\n",
    "\n",
    "To find optimum size of the cell we run the calculator with `optcell=1` parameter (plus some other connected parameters). After the run is finished we read-in the results and update the crystal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "2233d51c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial cell:\n",
      "   [ 3.84627441  3.84627441  3.84627441 60.         60.         60.        ]\n",
      "Initial stress:\n",
      "   [-0.00080608 -0.00080608 -0.00080608  0.          0.          0.        ] GPa\n"
     ]
    }
   ],
   "source": [
    "print('Initial cell:\\n  ', Si.cell.cellpar())\n",
    "print('Initial stress:\\n  ', Si.get_stress()/un.GPa, 'GPa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "0a0042c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-214.171809130701 [-0.00080608 -0.00080608 -0.00080608]\n"
     ]
    }
   ],
   "source": [
    "# Set up the calculation\n",
    "Si.calc = set_params(Si.calc, \n",
    "                     optcell=1,              # Optimize cell volume\n",
    "                     ionmov=2,               # Update atomic positions\n",
    "                     ntime=20,               # Max. number of steps\n",
    "                     ecutsm=0.5*un.Ha,       # Energy buffer for changing basis\n",
    "                     chkdilatmx=0,\n",
    "                    )\n",
    "\n",
    "# Run the calculation and check the result\n",
    "print(Si.get_potential_energy(), Si.get_stress()[:3]/un.GPa)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0cd06c5",
   "metadata": {},
   "source": [
    "Update the crystal from the `Si.calc.results` dictionary and check if we have acceptable convergence by switching off optimization mode (`optcell=0`) and recalculating the stress tensor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "f80c32ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update the cell matrix\n",
    "Si.set_cell(Si.calc.results['atoms'].cell, scale_atoms=True)\n",
    "\n",
    "# Reset the calculator to the *single point* mode\n",
    "Si.calc = set_params(Si.calc, optcell=0, ionmov=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "621d4572",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final cell:\n",
      "   [ 3.84627443  3.84627443  3.84627443 60.         60.         60.        ]\n",
      "Final stress:\n",
      "   [-0.00080502 -0.00080502 -0.00080502  0.          0.          0.        ] GPa\n"
     ]
    }
   ],
   "source": [
    "print('Final cell:\\n  ', Si.cell.cellpar())\n",
    "print('Final stress:\\n  ', Si.get_stress()/un.GPa, 'GPa')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "da618293",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final lattice constant a = 5.4395 A\n"
     ]
    }
   ],
   "source": [
    "print(f'Final lattice constant a = {2*Si.cell[0,2]:.4f} A')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d1112b0-b206-4289-89ee-ef4d6fbfeb58",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}