diff --git a/docs/_config.yml b/docs/_config.yml index af9f464..6d6fe35 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -11,6 +11,7 @@ execute: execute_notebooks: force exclude_patterns: - 'tour.ipynb' + - 'batching.ipynb' # Define the name of the latex output file for PDF builds latex: diff --git a/docs/_toc.yml b/docs/_toc.yml index 2642e7f..3b2d8f1 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -9,6 +9,7 @@ parts: - file: tour - file: prologue - file: optim + - file: batching - file: quirks - caption: Miscellaneous Things that Might be Useful chapters: diff --git a/docs/batching.ipynb b/docs/batching.ipynb new file mode 100644 index 0000000..5086187 --- /dev/null +++ b/docs/batching.ipynb @@ -0,0 +1,196 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "tags": [ + "remove-input" + ] + }, + "outputs": [], + "source": [ + "# Copyright (c) 2024 Graphcore Ltd. All rights reserved." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/ubuntu/miniforge3/envs/jax/lib/python3.10/site-packages/pyscf/dft/libxc.py:771: UserWarning: Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, corresponding to the original definition by Stephens et al. (issue 1480) and the same as the B3LYP functional in Gaussian. To restore the VWN5 definition, you can put the setting \"B3LYP_WITH_VWN5 = True\" in pyscf_conf.py\n", + " warnings.warn('Since PySCF-2.3, B3LYP (and B3P86) are changed to the VWN-RPA variant, '\n" + ] + } + ], + "source": [ + "import jax\n", + "import jax.numpy as jnp\n", + "import numpy as np\n", + "import optax\n", + "import seaborn as sns\n", + "from tqdm.notebook import tqdm\n", + "\n", + "from mess import Hamiltonian, basisset\n", + "from mess.structure import Structure, nuclear_energy\n", + "\n", + "sns.set_theme(style=\"whitegrid\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "def h2_hamiltonian(r: float, basis_name: str = \"sto-3g\", xc_method=\"lda\"):\n", + " mol = Structure(\n", + " atomic_number=np.array([1, 1]),\n", + " position=np.array([[0.0, 0.0, 0.0], [r, 0.0, 0.0]]),\n", + " )\n", + " basis = basisset(mol, basis_name=basis_name)\n", + " return Hamiltonian(basis, xc_method=xc_method)\n", + "\n", + "\n", + "num_confs = 32\n", + "rs = np.linspace(0.6, 6, num_confs)\n", + "H = [h2_hamiltonian(r) for r in rs]\n", + "num_orbitals = H[0].basis.num_orbitals\n", + "H = jax.tree_map(lambda *xs: jnp.sqtack(xs), *H)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "@jax.jit\n", + "@jax.vmap\n", + "def energy(Z, H):\n", + " C = H.orthonormalise(Z)\n", + " P = H.basis.density_matrix(C)\n", + " return H(P)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "Z = jnp.tile(jnp.eye(num_orbitals), (num_confs, 1, 1))\n", + "optimiser = optax.adam(learning_rate=0.1)\n", + "state = optimiser.init(Z)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "0fc26eb61c9a494fae0bb13ddcde1548", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + " 0%| | 0/128 [00:00" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "history = jnp.stack(history)\n", + "ax = sns.lineplot(history)\n", + "ax.set_xlabel(\"Iteration\")\n", + "ax.set_ylabel(\"Batched Loss (Hartree)\");" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAG9CAYAAAAIgELcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABhDklEQVR4nO3deVxU5eIG8GdW9l0FWVwARREX3FBzyaVcy6WoTFMrTSu10rqJmWlZervdcin95dKuqZl5W0xNLXdJU1wBBQQBFxAEBphhtvP7A5kkUIdhZg4Dz/fz6VPMnJl5eFN5fM973iMRBEEAEREREd2VVOwARERERI6ApYmIiIjIDCxNRERERGZgaSIiIiIyA0sTERERkRlYmoiIiIjMwNJEREREZAaWJiIiIiIzyMUOUJ+cPHkSgiBAoVCIHYWIiIjMpNPpIJFIEB0dfdfjONNkRYIgoD5vsC4IArRabb3+Hu2FY2ldHE/r4VhaF8fTemw5lub+/OZMkxVVzDC1b99e5CS2UVpaisTERISHh8PV1VXsOA6NY2ldHE/r4VhaF8fTemw5lmfOnDHrOM40EREREZmBpYmIiIjIDCxNRERERGZgaSIiIiIyA0sTERERkRlYmoiIiIjMwNJEREREZAaWJiIiIiIzsDQRERERmYGliYiIiMgMLE1EREREZnD40rR37148/PDDaN++PQYPHozvv//+nq9ZsWIFIiIiqv1n/vz5dkhNREREjsahb9h7/PhxTJ8+HY8++ijmzp2Lo0eP4o033oCbmxuGDBlyx9fFxsaiT58+lR47duwYPvjgA/Tt29fWsYmIiMgBOXRpWrVqFTp06IC3334bANCjRw9kZmZi+fLldy1NAQEBCAgIqPTYxo0b4eXlVSdLU1p2IY4nXsfo+8OhkDv85CAREZFDctifwFqtFvHx8VXK0bBhw5CamoqsrCyz36usrAy//fYbBg8eDKVSae2otbZhZxK+/jURxxOviR2FiIiowXLY0nT58mXodDqEhoZWejwsLAwAkJaWZvZ7/f777yguLsaIESOsmtFa5LLy/015hRqRkxARETVcDnt6rrCwEADg6elZ6fGKryueN8fPP/8Mf39/dOvWrda5BEFAaWlprd/ndq5Ot0pTQYnV37sm1Gp1pX+T5TiW1sXxtB6OpXVxPK3HlmMpCAIkEsk9j6tTpUmlUiEnJ+eex4WEhFjtM4uKirBv3z6MHz8eUmntJ950Oh0SExOtkOxvmtLyApiRdR2JiVqrvrcl0tPTxY5Qb3AsrYvjaT0cS+vieFqPrcbSnOU5dao07dixA/Pmzbvncdu3b4eXlxeA8qJ1u6KiIgAwPX8vO3fuhFarxUMPPVTDtNVTKBQIDw+3yntVSLuZgQPnVJA7uaNt27ZWfe+aUKvVSE9PR4sWLeDi4iJajvqAY2ldHE/r4VhaF8fTemw5likpKWYdV6dKU2xsLGJjY806VqvVQqFQIC0trdL2ARVrmf651ulOfv75Z4SGhiIyMrLmgashkUjg6upqlfeq0MjHHQCgLjNa/b0t4eLiUidy1AccS+vieFoPx9K6OJ7WY4uxNOfUHODAC8GVSiViYmKwc+fOSo9v374dYWFhCA4Ovud75OTk4M8//6yzC8AreLqVTxkWlYh/ao6IiKihctjSBADPP/88EhISsGDBAsTHx2P58uX4+eefMWPGjErHRUZGYu7cuVVev337dhiNRqudmrOVv0tTmchJiIiIGq46dXquprp27YoVK1Zg6dKl2LJlCwIDA7Fo0SIMHTq00nEGgwFGo7HK63/66Sd06NABzZo1s1dki3i6OQEon2kyd4U/ERERWZdDlyYAGDhwIAYOHHjXY5KTk6t93Jz71NUFFTNNWr0RZVoDnJ0c/n8bERGRw3Ho03MNhbNSZtrgkuuaiIiIxMHS5AAkEgkXgxMREYmMpclBmEpTKUsTERGRGFiaHARnmoiIiMTF0uQguO0AERGRuFiaHARnmoiIiMTF0uQgbt+riYiIiOyPpclBeLgpALA0ERERiYWlyUFUzDSpWJqIiIhEwdLkILimiYiISFwsTQ6CpYmIiEhcLE0O4vbSJAiCyGmIiIgaHpYmB1FRmvQGI9RlepHTEBERNTwsTQ7CWSmHUiEDwFN0REREYmBpciCertx2gIiISCwsTQ6EG1wSERGJh6XJgVSsa1KVsjQRERHZG0uTA+G2A0REROJhaXIgLE1ERETiYWlyICxNRERE4mFpciB/l6YykZMQERE1PCxNDsSDM01ERESiYWlyIDw9R0REJB6WJgdSsU+TiqWJiIjI7liaHAhv2ktERCQeliYHUrGmyWAUUKrhTXuJiIjsiaXJgTgpZHBW8qa9REREYmBpcjDcdoCIiEgcLE0OhtsOEBERiYOlycF4urI0ERERiYGlycGYth0oZWkiIiKyJ5YmB+PpzpkmIiIiMbA0ORjuCk5ERCQOliYHw9JEREQkDpYmB8PSREREJA6WJgfj4cp9moiIiMTA0uRgONNEREQkDpYmB1NRmlSlOhiNvGkvERGRvbA0OZiK0mQ0CijV6EROQ0RE1HCwNDkYhVwGFyc5AJ6iIyIisieWJgfEdU1ERET2x9LkgFiaiIiI7I+lyQF5uHHbASIiIntjaXJAnGkiIiKyP5YmB8TSREREZH8sTQ6IpYmIiMj+WJockKebEwCWJiIiIntiaXJAnGkiIiKyP5YmB8TSREREZH8sTQ7I05WliYiIyN5YmhxQxUxTsVoLA2/aS0REZBcsTQ6oYnNLQQBK1LxpLxERkT2wNDkguUwKN+eKm/ZyV3AiIiJ7YGlyUNx2gIiIyL5YmhwUr6AjIiKyL5YmB+XB0kRERGRXLE0OijNNRERE9sXS5KBYmoiIiOyLpclBVZQmFUsTERGRXbA0OSjONBEREdkXS5OD+rs0cZ8mIiIie3D40rR37148/PDDaN++PQYPHozvv//erNdduHABU6dORY8ePdC1a1eMGzcOR48etXFa6+E+TURERPbl0KXp+PHjmD59Ojp16oQ1a9Zg6NCheOONN7Bjx467vi4/Px+TJk1CQUEB3n33XXz44YdwdXXFlClTkJycbKf0tcPTc0RERPYlFztAbaxatQodOnTA22+/DQDo0aMHMjMzsXz5cgwZMuSOrzty5Ajy8vKwefNmBAcHAwC6d++O7t27Y/fu3YiIiLBL/tr4+6a9OhgMRshkDt1/iYiI6jyH/Umr1WoRHx9fpRwNGzYMqampyMrKuuNrdbrym9x6eHiYHnNycoJCoYAgCLYJbGXuLgrTf6tKedNeIiIiW3PYmabLly9Dp9MhNDS00uNhYWEAgLS0NNMs0j/1798fjRo1wpIlS/DKK69ALpfjs88+g0QiwciRI2uVSxAElJaW1uo9zOXmIkeJWo+cvEIoZe42/zy1Wl3p32Q5jqV1cTyth2NpXRxP67HlWAqCAIlEcs/jHLY0FRYWAgA8PT0rPV7xdcXz1fHy8sL69esxdepU9OnTBwDg7e2NNWvWICQkpFa5dDodEhMTa/Ue5nKSAyUAziamQJXnZJfPBID09HS7fVZ9x7G0Lo6n9XAsrYvjaT22GkulUnnPY+pUaVKpVMjJybnncbUtNnl5eZg+fTqaNWuGuXPnQiaTYfPmzXj++eexfv1602yVJRQKBcLDw2uVz1yNDqiQryqEb+NAtG3bxOafp1arkZ6ejhYtWsDFxcXmn1efcSyti+NpPRxL6+J4Wo8txzIlJcWs4+pUadqxYwfmzZt3z+O2b98OLy8vAOVF63ZFRUUAYHq+OmvXrkVhYSG2bt1qapY9e/bE8OHDsXLlSvz3v/+19FuARCKBq6urxa+vCW8PFwCFKNPDbp8JAC4uLnb9vPqMY2ldHE/r4VhaF8fTemwxluacmgMsKE1qtRqHDh3CiRMnkJqaips3b0IikcDHxwehoaHo3LkzevXqZdE3FBsbi9jYWLOO1Wq1UCgUSEtLM51iA8rXMgGostbpdikpKQgNDa00FSeTyRAREYHLly/XOLdYuO0AERGR/ZhdmpKTk/H5559j165dKC0thbOzMwICAuDl5QVBEHDp0iUcOXIEn332GVxcXDB48GA8/fTTNrt8X6lUIiYmBjt37sTEiRNNj2/fvh1hYWF3XAQOAIGBgdizZw/Kysrg5FS+FshgMCApKQlt27a1SV5bYGkiIiKyH7NK08svv4xdu3YhKioKM2bMQK9evRAeHg6ZTFbpOIPBgJSUFBw6dAg7d+7E6NGjMWTIEHz44Yc2Cf/8889jwoQJWLBgAYYOHYr4+Hj8/PPP+OijjyodFxkZiVGjRuG9994DUD6jtWXLFrzwwgsYN24cZDIZNm3ahIyMDCxatMgmWW3Bg6WJiIjIbswqTVKpFN9///09Z2EqTnFFRETgmWeeQWJiItasWWOVoNXp2rUrVqxYgaVLl2LLli0IDAzEokWLMHTo0ErHGQwGGI1G09dRUVFYu3YtVq5cibi4OBiNRoSHh2P16tXo1q2bzfJaG2eaiIiI7Mes0mTpTFHbtm1tNstUYeDAgRg4cOBdj6nu1ig9e/ZEz549bRXLLipKk4qliYiIyOYcdkdw4kwTERGRPVlcmoqLi7F69Wo8++yzGDVqFE6fPg0AKCgowOeff46MjAyrhaTq/V2aykROQkREVP9ZtE/TtWvXMH78eFy7dg3NmzdHWloaSkpKAJTvrL1x40ZkZ2ebtecSWc7TrfzKvxKNHnqDEXLetJeIiMhmLPop+/7776OkpATbtm3D119/XeUmt4MGDcKRI0esEpDuzM1FAemt/bi4romIiMi2LCpNhw4dwlNPPYXw8PBqd9EMCQnB1atXax2O7k4mlcDNheuaiIiI7MGi0qTRaODr63vH5ytO1ZHtmdY1lbI0ERER2ZJFpSksLAzHjh274/O7d+9GZGSkxaHIfLyCjoiIyD4sKk0TJ07E9u3bsXr1ahQXFwMABEFARkYGXnvtNSQkJGDSpEnWzEl3wNJERERkHxZdPTdy5EhcuXIFy5Ytw9KlSwEAkydPhiAIkEqleOWVVzBo0CBr5qQ74LYDRERE9mFRaQLK7/s2cuRI7Nq1CxkZGTAajWjWrBkefPBBhISEWDMj3QVnmoiIiOzD4tIEAIGBgTwNJ7KKvZpYmoiIqD4QBAFFJVpczy/FtbwSXMsr//fVG8XQa0vxRrgBriJlq1VpSkhIQHx8PPLy8vDkk0+iRYsWUKvVSEtLQ4sWLeDm5matnHQHnm4KACxNRETkOLQ6A67nl5b/k1eCa7cVpOv5JVCXGap9nUQCFBRr4e1l58C3WFSatFotZs2ahT179kAQBEgkEvTv3x8tWrSAVCrFM888g0mTJuH555+3dl76h4qZJm5uSUREdYnRKCC3QI3snGJk5aqQlVOM7JxiZOcWI69Qc8/X+3k5I8DPDf6+rgjwdYWPhxySshto4uNih/TVs6g0LVu2DH/88QcWLFiAmJgYDBkyxPSck5MThgwZgj179rA02QHXNBERkZhKNTpk55YXoqycYmTd+u8rucXQ6o13fJ2Lkwz+vm4I8HNFgJ8bAnxd4X+rJPn7ukKpkFX+nNJSJCYW2vrbuSuLStMvv/yCJ554Ao8//jhu3rxZ5fmwsDDs2LGj1uHo3liaiIjIHtRlemRcK0L6lSKkXy1C5vXy2aP8ojvPGsllUjRt5IbgJu4IbuKOoMbuCGrijqZ+bvB0U1Z7V5G6zKLSlJeXh4iIiDs+L5PJoNHce+qNaq+iNKnL9NDpDVDIZfd4BRER0Z0ZjAKu5ZUg/UoRLl0tRMbV8pJ0La/0jq/x9nBCUGP3KuXI38cVsnp0M3mLSlPTpk2RlpZ2x+dPnDiBZs2aWRyKzOfqrIBUKoHRWH61gZ+XeOd6iYjIsRQWlyH9VilKv1KE9GtFuHxNBa2u+oXYvp5OaNHUC82beqJ5gEd5QWriAXcXhZ2Ti8Oi0jRixAh8/vnnePDBB9GiRQsAME2xbd68Gb/++itmz55ttZB0Z1KpBJ6uShQUl7E0ERHRHRWrdUjNLMCFzJu4mFmAlKwC5N5UV3usUiFD8wAPtGjqWf5PoCeaB3jCy93JzqnrFotK07Rp03Dq1CmMHz8eoaGhkEgkWLx4MQoLC3Ht2jX069eP+zfZkYebwlSaiIiINGV6pGYXlpejzAJczLyJKzdKqj22qZ+bqRS1CPREy6ae8Pdzg0zqWOuN7MGi0qRUKrF27Vr8+OOP2LlzJ4xGI7RaLSIiIvDyyy9j5MiRDre4y5GVbztQDFUpSxMRUUOj0xtw6UoRUrIKcPFyeUHKvK6CUah6rL+vK8JDvNE6xButQnwQFuwFV+eGcWrNGmpcmjQaDT766CPExMRg5MiRGDlypC1yUQ3wCjoiooajRK1DUkY+zqXl4fylfFy8fLPaS/t9PZ3RKsT71j/lBamhn16rrRqXJmdnZ2zatAnh4eG2yEMWYGkiIqq/8os0twpSHs6n5SP9amGVWSQPVwVahfigVYg3wm8VJa5xtT6LTs+1a9cOFy5csHYWshBLExFR/SAIArJzi3H+Ur6pKFV3qX+AnysiW/qhXagfIlv6IqixO5fF2IFFpWnu3Ll47rnn0Lp1a4wePRpyea1uYUe1ZCpNxSxNRESORBAEXLlRgoTkHJxKuYHzl/JQ+I8/yyUSoGVTL0SG+iKyZXlJ4iySOCxqO3PmzIFEIsH8+fOxaNEi+Pv7w8mp8nlSiUSCH3/80Soh6e7+nmkqEzkJERHdS1GJFqcu5iLhQi4SLuQg5x+X/SvkUrRu5oPIlr5oF+qHNs194dZA9kGq6ywqTd7e3vD29kbLli2tnYcs4OF6qzTx6jkiojpHbxBwNi0fSZfTcfJCLlKzCiDctiZJLpMgsqUfOrZqjPZhjRAe4sW7O9RRFpWmr7/+2to5qBYqZppUXNNERCQ6QRBw+ZoKJy/k4q/EqziXlg+dIbvSMc0CPBDdugk6tW6MqFA/ODtxmYsjsOj/0rZt29C1a1cEBwdX+3x2djaOHTuGUaNG1SYbmal8nyYuBCciEoumTI+TF3IQf+4aTibnIL+o8nIJL3clols3QXREY3Rs1ZhrkhyURaUpLi4O77///h1L06lTpxAXF8fSZCcVM00arQFlOgOcFJzWJSKytZsqDY6dv46jZ6/i1IXcSnslKeVStAv1Q7tQb7hLCnF/zw5wc3MTMS1Zg0WlSRCq2Wb0NqWlpZDJ+IPbXlyd5ZBJJTAYBahKtHDy5t9giIhsIStHhfiz1xB/7hqSMvIrrU1q4uuKHu0C0C3SH5Et/aBUyFBaWorExERuB1BPmF2akpKSkJSUZPr6+PHjMBiq3gW5qKgIGzdu5CJxO5JIJPB0U+Kmqvz+c41YmoiIrMJoFJCccRPx567i6NlryM4trvR8WLAXekQ1RUy7ALRo6slyVM+ZXZp2796Njz/+GED5D+lNmzZh06ZN1R7r6emJf//739ZJSGb5uzRx2wEiotrQ6Q04eSEX8Wev4c/z11Cg+vvPVZlUgvbhjdCjXQC6t2uKxj78S2pDYnZpeuyxx3D//fdDEATExsZi5syZ6Nu3b6VjJBIJXFxc0KxZM254aWce3BWciMhiBqOAsyk3sO9kFg6fvoISjd70nKuzHF3b+CMmKgBd2vhzz6QGzOxm06RJEzRp0gQ6nQ5xcXF48MEH0bRpU1tmoxrgtgNERDUjCAIuZhZg38ksHEzIrnTFm6+nM3pEBSAmqinahzWCQi4VMSnVFTWeDpJKpXj//fchkUgwYcIEW2QiC3DbASIi82TlqLDvRDb2nczC1RslpsfdXRS4r2Mg+nUORruWfpBKuT6JKqtxaZLJZAgMDIRWyx/OdQlv2ktEdGd5hWrsP1lelFKzCk2PKxUy9GgXgH6dgxEd0YQzSnRXFi08Gj9+PNavX49HH30U3t7eVo5ElmBpIiKqTFWqxeHTV7DvRDbOpt0wbQ8glUrQOaIJ+kUHISaqKVy4GzeZyaJfKUajEUqlEg888AAGDx6MoKAgODs7VzpGIpFg0qRJ1shIZmBpIiIqX6d0Ni0PO49k4NDpK9Ab/t5wsl2oH/pGB+G+DoHwcne6y7sQVc+i0nT7dgJbtmyp9hiWJvtiaSKihqywuAx7j2di59F0ZOf+vU6pRVNP3N85GH2ig9DEx1XEhFQfWFSa9uzZY+0cVEserhWlifs0EVHDIAgCzqTewM4jGTh85qppVsnFSYZ+nUMwOKY5wkO8xQ1J9YpFpSkoKMjaOaiWTDNNpTqRkxAR2VaBqgx7jl3GzviMSle/tQrxxuAeLdA3OojrlMgm+KuqnqgoTVqdARqtHs5K/q8lovrDaBRwOiUXO45mIP7sVegN5au6XZzkuL9LMAbHNEdYsLe4Ianes/gna1JSEr755hucP38eKpUKRqOx0vMSiQS7d++udUAyj4uTHHKZFHqDEUUlWpYmIqoXClRl+O3PDOyKz8C1vFLT4xHNfDC4R3P07sRZJbIfi36lxcfHY/LkyfDy8kJUVBTOnz+PHj16oKysDAkJCQgPD0dUVJS1s9JdVNy0N79Ig6ISLRc8EpFDy7yuwv/2p2Lv8Uzo9OV/KXd1lqN/lxAM7tEcLQO9RE5IDZFFpWn58uUICQnB5s2bodVq0atXL0ydOhU9e/bEqVOnMGXKFLz66qvWzkr3cHtpIiJyNBULu3/4IxXHE6+bHm8V4o1hvVqid8dAOHNWiURk0a++8+fPY8aMGXB3d0dhYfnOqhWn5zp27IjHH38cy5YtQ79+/ayXlO6J2w4QkSPSG4w4eOoKtu1LMe3WLZEAMe0CMKpfOCJb+kIi4S1NSHwWlSaZTAY3NzcAgKenJ+RyOfLy8kzPh4SEIDU11ToJyWwebtx2gIgcR4lah51HM/DTgVTcKNQAKL+tycBuIRjVNwyBjd1FTkhUmUWlqVmzZkhPTwdQvpYmNDQUu3fvxsMPPwwA+OOPP9CoUSOrhSTzVMw0qUq47QAR1V05+aX48UAadsVnQF2mBwB4ezhhxH0tMaRnC+7WTXWWRaWpX79++P777zF79mzI5XI8/fTTiIuLw4MPPggAuHz5MmbNmmXVoHRvnpxpIqI67GLmTWz7IxUHT1+B0Vi+ZUCIvwdG9wtDv87BUCpkIickujuLStMLL7yACRMmQCYr/wU+evRoSKVS7Nq1CzKZDNOmTcOYMWOsGpTujWuaiKiuEQQBp1NuYONvyTib+vcyjo6tGmH0/eHoHNGE65XIYVhUmhQKBXx8fCo9NnLkSIwcOdIqocgynm7lU9osTURUF5xJvYENO5NMZUkmlaBvdBBG9QtHaBC3DCDHw2s36xHONBFRXXAuLQ8bdibhdMoNAIBcJsWQns3xSP9WaOTtInI6IsuZXZoWLVpU4zefN29ejV9DlmNpIiIxJaXnY/3OJCRcyAUAyGUSPBjTHLEDW7MsUb1gdmn65ptvqjwmkUggCEK1x0skEpYmO/N0/bs0CYLAdQJEZBcXLt/E+p1JOJGUA6D8NNyg7s3w2KDWvDsB1Stml6akpKRKX+fn56NXr174/PPP0bNnT6sHo5qrmGnSG4zQaA28HxMR2VRKZgHW70wy7d4tlUowsGsIHn8gAv6+LEtU/1j8U5WzGHWPk1IGpVwKrb78pr0sTURkC2nZhdiwMwnx564BAKQSoH/XEDw+KAJNG7mJnI7IdvhTtR6puGnvjUINikrK+Dc9IrKqjGtFWL8jCUfOXAVQXpb6dQ7GEw9EcPduahBYmuoZTzenW6WJi8GJyDoKi8uwfkcSdh5Nh1Eovy9cn05BeOKBCIT4e4gdj8huWJrqGV5BR0TWotMb8NOBNGzafQGlmvLbnfRs3xTjh7RBswBPkdMR2V+tSxPXNtUtLE1EVFuCIODP8zlYv+siruWVAgDCgr0w+eEoRIXxvqLUcJldmqKjo6stSNOmTYNUKq3yuEQiwV9//VW7dGbYu3cvli5dikuXLiEwMBDPPfccHnnkkXu+LjU1FUuWLMGxY8egUChw//33Iy4uDr6+vjbPbEseLE1EVAuXrhThiz25yMjJBgD4ejrhqaGRGNA1BFIp/5JMDZvZpWnw4MF1blbp+PHjmD59Oh599FHMnTsXR48exRtvvAE3NzcMGTLkjq8rLi7GxIkT4e/vjw8++AAajQYffvghpk6dik2bNlVbAh1FxUyTiqWJiGogv0iDr7cnYs/xyxAEQCGXYkz/cDzSvxWvxCW6xezfCUuWLLFlDousWrUKHTp0wNtvvw0A6NGjBzIzM7F8+fK7lqYNGzZApVJh27ZtaNSofKq5efPmePTRR7Fnzx488MADdslvCzw9R0Q1UaYzYNu+FGzZcxEarQEA0L65C6Y92gXNAv1ETkdUtzjslIpWq0V8fHyVcjRs2DCkpqYiKyvrjq89f/482rRpYypMANC+fXt4e3tj7969NstsDyxNRGQOQRCw/2QWnv/3HnzzaxI0WgMimvvgnee64ZH7/HjbE6JqmDXTdPLkSURHR1v0AbV57d1cvnwZOp0OoaGhlR4PCwsDAKSlpSE4OLja15aVlUGpVFZ5XKlUIi0trVa5BEFAaWlprd6jNpTy8tvaFKg0Vs+hVqsr/Zssx7G0Lo5nzVzMLMRXvybjQmYhAMDPyxnjHgxHr/YB0Gg0SC/mWFoLf21ajy3H0txbj5lVmiZOnIiOHTti7Nix6N+/P1xc7v43kJKSEuzduxcbN27E2bNncerUKfNS10BhYflvdk/Pype9Vnxd8Xx1WrRoga1bt0Kj0cDZ2RkAcOXKFeTm5sLVtXYbQup0OiQmJtbqPWrjxs3yGaabRaU2y5Genm6T922IOJbWxfG8O7XWiN9OFuJEagkAQCGXoHekB3q18YBCXoCkpALTsRxL6+J4Wo+txrK6yZR/Mqs07dy5E5988gn+9a9/QaFQoEOHDoiMjERwcDC8vLwgCAKKioqQlZWFs2fP4vTp0zAYDBg5ciQ++OADswOrVCrk5OTc87iQkBCz37M6sbGx+OqrrzB//nzMnj0bGo0Gb775JqRSaa0XuysUCoSHh9fqPWqjSaEG+DUHaq2ANm3aWHXxvlqtRnp6Olq0aHHP4kx3x7G0Lo7n3QmCgKNnr+PzHckoLC7/i1W/6EA8MSgMvp7OlY7lWFoXx9N6bDmWKSkpZh1nVmlq2rQpFi1ahFmzZuHHH3/Enj178O2330Kj0VQ6ztnZGVFRUXj55ZcxcuTIGl++v2PHDsybN++ex23fvh1eXl4AyovW7YqKigDA9Hx1QkND8e677+Ldd9/F//73PwDAgw8+iL59+6KkpKRGmf9JIpHUeraqNmQKJwCAwSgAUiVcXRRW/wwXFxdRv8f6hGNpXRzPqnJuluL/tp7BsfPlN9UNbuKO6bGd0C707ou8OZbWxfG0HluMpbkTDDW6jtTX1xeTJk3CpEmToNfrcfXqVdy8eRMA4OPjg6ZNm0Iut/zS1NjYWMTGxpp1rFarhUKhQFpaGvr06WN6vGJN0j/XOv3TqFGjMGzYMKSnp8PLywv+/v4YPnw4BgwYYHH+usBJIYOTUoYyrQGqUi3cbFCaiKjuMxgF/HIwDV//mgiN1gC5TIrHBrbCowNbQSGXiR2PyCFZ3HDkcjlCQkJqfarMUkqlEjExMdi5cycmTpxoenz79u0ICwu74yLwf75H69atAQBHjhxBeno6Ro8ebbPM9uLppkSuVo2iEi0C/HjHcaKGJi27ECu+S0BKZgEAILKlL6bHduJ94ohqyaF3LHv++ecxYcIELFiwAEOHDkV8fDx+/vlnfPTRR5WOi4yMxKhRo/Dee+8BAEpLS7FixQp069YNTk5OSEhIwOrVqzF9+vR7zlA5Ak83JXJvqrntAFEDo9HqsXFXMn7YlwqjUYCbsxyTRrTDgzHNuZs3kRU4dGnq2rUrVqxYgaVLl2LLli0IDAzEokWLMHTo0ErHGQwGGI1G09dSqRQXLlzA1q1bUVpaitDQULz11lsYM2aMvb8Fm/B0rdirqUzkJERkLyeSc7Byyylczy/fauS+joF4blT7Kgu9ichyDl2aAGDgwIEYOHDgXY9JTk6u9LWzszPWrVtny1ii8nQrXwzOmSai+q+wuAxrfzyLP/4q39C3kbcLnh/TAd3bBYicjKj+cfjSRFV5unNXcKL6ThAE7D2eiXU/noWqVAeJBHiodyjGDWkDV2deAEJkCyxN9RBvpUJUv90oUGPZppNIuJALAGgZ6InpsZ3QupmPyMmI6jeLStPq1asxcuRI+Pv7WzsPWYGHK0sTUX118FQ2PvnuFIrVOijlUjw5uA1G9guDXOawtxIlchgWlaalS5di6dKl6Nq1K0aOHInBgwfD3d3d2tnIQpxpIqp/SjU6fPrDGew9ngkACA/2wuxxXRDchNsIENmLRX81+f333zFr1iwUFhbijTfeQO/evfHKK6/gjz/+gMFgsHZGqiGWJqL65VxaHmb89w/sPZ4JqQR4bFBr/GdmXxYmIjuzaKbJ398fkydPxuTJk3HhwgX89NNP+OWXX/Drr7/Cx8cHw4YNw8MPP4yOHTtaOy+ZoaI0qViaiByaTm/Et7uS8P3eizAKgL+vK2Y92RmRLe9+CxQiso1aLwRv3bo1Zs+ejdmzZ+P48eP48ssvsWHDBmzYsAHNmjXDyJEj8fjjj8PPj7/J7cU001SqhdEocFM7IgeUeV2FDzf8hZSsQgDAwG4heG5Ue14ZRyQiq1w9V1ZWht27d+Onn37CwYMHIZPJcN9990GhUGDlypVYs2YN3n//fTzwwAPW+Di6h4rSZDQKKNXo4H5rYTgR1X2CIODXI+lY9+M5aHUGeLgq8OKjnXBfx0CxoxE1eBaXJkEQcOjQIfz000/YvXs3SkpKEBkZiVdffRUPPfSQaWYpJycHs2fPxpIlS1ia7EQhl8HFSQ51mR5FJVqWJiIHcbNIg+WbE3A88ToAoFPrxnj5iWj4ebmInIyIAAtL03vvvYft27cjLy8PjRs3xhNPPIFRo0ahVatWVY5t0qQJHn30Ubz++uu1Dkvm83BTlpemUi3491Oiuu/o2atYsTkBRSVaKORSTBoeiRG9Q3l6nagOsag0fffddxg0aBBGjRqFXr16QSK5+2/qLl26YPHixRYFJMt4uimRk1/KK+iI6jh1mR7rfjyLnUczAJRvVDn7yS5o3tRT5GRE9E8WlaZDhw7B1dXV7OODg4MRHBxsyUeRhUyLwYtZmojqqpSsArz/9XFcvVECiQQY3S8c44e2gUIuEzsaEVXDotJUk8JE4uBeTUR12674DPzf1tPQ6Y1o5OWMV57sjA7hjcWORUR3YVFpmjBhwl2fl0gkcHJyQkBAAGJiYjB48GDI5bzNnT39XZrKRE5CRLcr0xnw6dbT+O3PywCA7pEBeGVsNC/YIHIAFjUZQRBw/fp1XL58GV5eXggKCgIAZGdno7CwEM2bN4e7uztOnTqFzZs3Y/Xq1fj888/h6+tr1fB0Z5xpIqp7ruWVYPEXx5B2pRBSCTB+aFs80r8VF3sTOQiLbqPy0ksvobCwEEuWLMHhw4exdetWbN26FYcPH8bixYtRWFiIN998E0ePHsV7772HlJQUfPjhh9bOTnfh6eYEgKWJqK748/w1vPzRPqRdKYSXuxJvP9cLsQNbszARORCLZpref/99jBkzBqNGjar0uEwmw+jRo3HhwgUsXrwYmzZtwpgxY5CQkIC9e/daIy+ZyfPWVL+qlKWJSEwGo4ANO5OwefcFAEBEcx/MmdANjby59xKRo7Fopik5OfmuV8MFBwcjKSnJ9HW7du1QWFhoyUeRhXh6jkh8hcVlWLD6iKkwjejdEotf6M3CROSgLCpNjRs3xo4dO2A0Gqs8ZzQa8euvv6JRo0amxwoKCuDl5WV5SqoxliYicSVn5OPlD/9AwsVcOClleHVcF0wd3QEKuUV/7BJRHWDR6bmnn34a77zzDsaOHYvY2Fg0a9YMAJCRkYHvvvsOZ86cwbx580zH79ixAx06dLBOYjJLRWkqLtXCYBQg47oJIrsQBAHbD6dj7f/OQG8QENTYDXGTuqN5ADerJHJ0FpWmcePGQSKRYPny5Zg3b55pR3BBEODt7Y158+Zh3LhxAACtVou4uDjTFXZkHx4VN+0VgBK1zlSiiMh2NGV6fLLlFP44kQUAuK9DIGY+3gmuzgqRkxGRNVi8edKTTz6J2NhYnD17FleuXAEABAYGIioqCgrF339AKJVKdO/evfZJqUbkMincnOUo0ehRVFLG0kRkY9m5xXjviz9x+ZoKUqkET4+IxMi+Yfe8zRQROY4alya1Wo37778fU6ZMweTJkxEdHY3o6GhbZKNa8nRzulWauK6JyJaOnLmCj749CXWZHj4eTnh9Qje0C/UTOxYRWVmNS5OLiwtkMhlcXHj1R13n4abA1TxAxdJEZBOCIGDzngv45tfyq4Xbhfrh9ae6wsfTWeRkRGQLFl3G8eCDD2Lnzp0QBMHaeciKuMElke3o9AYs3XjSVJge7hOKRdN6sTAR1WMWrWkaPnw4Fi5ciAkTJiA2NhZBQUFwdq76B0W7du1qHZAsx20HiGyjsLgMi788hnNpeZBKJZg2uj2G9mopdiwisjGLStNTTz1l+u/jx49XeV4QBEgkEiQmJlqejGqNpYnI+rJyVHh7bTyu5pXA1VmO1yd0Q+eIJmLHIiI7sKg0LV682No5yAZYmois69TFXCz+8hhK1Dr4+7pi/rMxaMb9l4gaDItK0+jRo62dg2yApYnIenYezcCq70/BYBTQtoUv3ni6O7zcncSORUR2ZPE+TRVycnKQn5+PZs2awdXV1RqZyEr+Lk1lIichclwGo4AvfzmPH/5IAQD0iw7GzMc7QamQiZyMiOzN4psg7d69G0OGDEG/fv0wevRonDp1CgCQn5+PUaNG4bfffrNaSLKMhytnmohqQ1Omx+Iv/jQVpicfjMDscZ1ZmIgaKItK0969ezFjxgz4+PjgxRdfrLT1gK+vL/z9/bF161arhSTLVNxJPeemGgZD1ZsrE9Gd5RWqMWflQcSfuwaFXIpXx3XB2MFtuMM3UQNmUWn65JNP0LVrV3z77beme8zdrlOnTrxyrg5o4uMKpUIGvcGIa/mlYschchgpWQWYtXQ/UrMK4eWuxLvT7kO/zsFixyIikVlUmi5evIihQ4fe8flGjRohLy/P4lBkHVKpBMFN3AEAmddVIqchcgxHz17FnE8OIr9IgxB/D3wwsy/atvQVOxYR1QEWlSYXFxeo1eo7Pp+ZmQlvb29LM5EVhTTxAMDSRHQvgiDghz9S8N4Xf6JMa0B068b4z4w+CPBzEzsaEdURFpWmmJgYbNu2DXq9vspzubm52Lx5M3r37l3rcFR7If7lM01ZOcUiJyGqu4xGAZ/+cAaf/XQOggAM7dUCb03uATcXhdjRiKgOsWjLgZdffhmPP/44Hn30UQwZMgQSiQQHDx7E0aNHsWnTJgiCgBdffNHaWckCwf6caSK6G53eiKXfnsD+hGxIJMCzD0fh4T6hXPBNRFVYNNMUGhqKDRs2wNvbG8uWLYMgCFi3bh0+/fRTtG7dGhs2bEBwMBdN1gXNbpWmrBwVb7BM9A+aMj0WfRaP/QnZkEkleHVcF4zsG8bCRETVsnhzy1atWuGLL75AYWEhMjIyIAgCQkJC4OvLBZN1SdNGbpBJJVCXGXCjQIPGPi5iRyKqE1SlWixcexTJGTfhpJQhbmI3dGnjL3YsIqrDar0juJeXFzp06GCNLGQDcpkUTRu5ISunGJk5KpYmIpTvwTR/9RFcvqaCu4sCb03ugTYt+Bc+Iro7i0uTwWDAwYMHkZmZicLCwiqnfiQSCdc11REh/h7IyilG1nUV78ZODd6V3GK8ufoIcvJL4evpjLen9kRz3nSXiMxgUWk6c+YMZs6ciWvXrt1xnQxLU91h2quJV9BRA5eaVYAFa46ioLgMTRu54Z2pveDvy3tmEpF5LCpNCxcuhEajMe0M7unJv6XVZc14BR0RzqTewKLP4lGq0SM0yAsLpvSAj4ez2LGIyIFYVJqSk5PxyiuvYMCAAdbOQzbAbQeooYs/exX//vo4dHoj2oX64c1nYrgHExHVmEWlKSAggJevO5DgxuWn54pKtCgsLoOXu5PIiYjsZ8+xy1i+OQFGo4CYdgF47amucFLIxI5FRA7Ion2apkyZgs2bN6O4mGtkHIGzkxxNbl01x53BqSHZti8FSzeehNEoYEDXEMRN7MbCREQWs2imqaSkBG5ubnjggQcwfPhwBAQEQCar/AeRRCLBpEmTrJGRrCDY3wM5N9XIvK5Cu1A/seMQ2ZQgCPj610R8t+ciAGBUvzA8PaIdpFJuWklElrOoNP373/82/fc333xT7TEsTXVLSBMPnEjKQWYO1zVR/WYwClj1/SnsPJoBAJg4PBKP9A/nLt9EVGsWlaY9e/ZYOwfZWEjF7VSu8/Qc1V96vRHLvz6OQ6evQCoBXni0Iwb3aCF2LCKqJywqTUFBQdbOQTYW4l++GPwyr6CjekpvEPDRptM4npQLuUyKV8d3wX0dAsWORUT1iNkLwU+fPo2CggKzjs3MzMS2bdssjES2UDHTdKNADXWZXuQ0RNal0xux+WAejiflQimX4s1nYliYiMjqzC5Njz/+OA4cOGD6uqCgAB07dsSff/5Z5diTJ08iLi7OOgnJKjxclfC+tdVAFtc1UT2i0xvw4bencCFbA4VcinnPxKBzG94uiIisz+zS9M99mQRBQFlZGQwGg9VDkW0E3zpFl8l1TVRPaHUGvPfFMZy4cANymQSvj++EaN5fkYhsxKJ9msgxmRaDc6aJ6gGtzoB3v/gTxxOvQ6mQ4sl+fmgfxu00iMh2LFoITo4ppEl5abp8jaWJHFuZzoB3P4vHyQu5cFLK8Pr4TpBqrosdi4jqOc40NSAVV9BxpokcmUarx6J15YXJWSnDW5N7oF1LX7FjEVEDUKOZpuzsbJw7dw4AoFKV/+DNyMiAp6dnpeOysrKsFI+sqeL03NW8Uuj0Rijk7MzkWDRaPRZ9Fo9TF2/AWSnDgik90S7UD6WlpWJHI6IGoEaladmyZVi2bFmlxxYuXFjlOEEQuPtuHeTr6QwXJznUZXpcuVGM5gGe934RUR2hKdPjnc/icTrlBlycZHhrck/eEoiI7Mrs0rR48WJb5iA7kEgkaObvgeTLN5F1naWJHIemTI+318XjTOoNuDjJsXBKT7TlKTkisjOzS9Po0aNtmYPsJNjfHcmXb/IedOQw1GV6LFx7FOfS8uDiJMfbz/VEmxYsTERkfw69qOXQoUOYPXs2Bg0ahIiICLz99ttmv1alUmHu3Lno3r07oqOjMXPmTOTk5Ngwbd1QcQVdJq+gIwdQqtFhwZojOJeWB1dnOd6ZysJEROJx6NJ04MABJCUloVu3blUWo9/Lyy+/jEOHDmHBggX44IMPcOnSJUyZMgV6ff2+xUjFYnDONFFdV16YjuL8pXy4OcvxztReiGjOwkRE4nHofZr+9a9/Yc6cOQCA+Ph4s1938uRJHDx4EOvWrUPv3r0BAC1btsSwYcOwa9cuDBs2zCZ564KKXcGzc4phMAqQSblgn+qeisKUmJ4PNxcF3pnaE61CfMSORUQNnEPPNEmllsXfv38/PD09cd9995keCw0NRdu2bbF//35rxauT/H3doJBLodUbkXuTl2lT3aO5tYYpMT0f7i4KLJrai4WJiOoEhy5NlkpLS0PLli2rbIsQGhqKtLQ0kVLZh0wqQVDjinvQ8RQd1S06ffmtUc5fujXDNK0XwkO8xY5FRATAwU/PWaqoqAgeHh5VHvfy8sLZs2dr9d6CINT5jfaa+rkg/WoRUrPy0a6F+WvB1Gp1pX+T5TiWVekNRizddBoJt26NMmd8JwT6Ks36/cTxtB6OpXVxPK3HlmNp7v6SZpWmY8eOWRSiW7duNTpepVKZdQVbSEgIlEqlRZlsTafTITExUewYd+UkKf8Fd+5iNlr51fwXX3p6upUTNVwcy3JGQcAPR27iTHopZFLg8d4+MJRcRWLi1Rq9D8fTejiW1sXxtB5bjaU5vcKs0vTUU0/VaIfvisZW0/KwY8cOzJs3757Hbd++HWFhYTV679t5enri2rVrVR4vLCyEl5eXxe8LAAqFAuHh4bV6D1sr0F/DH2fOoESnQNu2bc1+nVqtRnp6Olq0aAEXFxcbJqz/OJZ/EwQBa39KvFWYJJj9ZEd0iWhco/fgeFoPx9K6OJ7WY8uxTElJMes4s0rTV199Vasw5oqNjUVsbKzNPyc0NBRHjhypMh136dIltG7dulbvLZFI4OrqWtuINhXWrBEA4EpuCVxcXGp8yxsXF5c6/z06ioY+loIg4POfz2P3sWxIJMDsJ7ugT3SQxe/X0MfTmjiW1sXxtB5bjKW5PwfNKk3du3evVZi6pm/fvli5ciWOHDmCXr16ASgvTOfPn8fkyZNFTmd7QY3dIJUAJRo9bqrK4OvpLHYkaqA2/nYBP/xR/je86bGdalWYiIhszaGvnsvOzsaOHTuwY8cOqNVqXL582fT17SIjIzF37lzT19HR0ejduzfmzp2LX3/9FXv37sXMmTMRERGBBx980N7fht0p5DIE+LkB4BV0JJ5t+1KxYWcSAGDKyCg8GNNc5ERERHdn8dVzZWVl2LlzJ86fPw+VSgWj0VjpeYlEgvfee6/WAe8mPj4ecXFxpq8PHDiAAwcOAACSk5NNjxsMhir5li5disWLF2P+/PnQ6/Xo3bs35s2bB7m8YVxQGOLvgSs3SpB5XYWOrWq2foSotnYezcC6H8uvVB0/pA0e7mv5GkUiInuxqCFkZ2djwoQJyM7OhqenJ1QqFby8vKBSqWAwGODj42OXc7djxozBmDFj7nnc7QWqgoeHB9577z2bF7u6KriJO+LPcaaJ7G/fiSx8siUBAPBI/3A8Nqh26wiJiOzFotNz77//PoqLi7F582bs2LEDgiDgo48+wsmTJ/Hqq6/C2dkZ69ats3ZWsqKKe9Bl5RSLnIQakvizV/HhtycgCMDQXi0wcXhkjS9EICISi0Wl6ejRoxg7diw6dOhQ6VYmSqUSkydPRo8ePRrsDI6jMN24lzNNZCenLuTi318fh9EooH+XYEwb3YGFiYgcikWlSaPRICio/CoXd3d3SCQSqFR///CNjo7GX3/9ZZ2EZBPBTcpvpXJTVYbiUq3Iaai+S7yUj3c+j4dOb0TP9k3x0uPRkPJm0UTkYCwqTU2bNsX169cBAHK5HP7+/khISDA9n5KSAicnJ6sEJNtwdVagkVf5VgOZ13mKjmwnNasAC9ceQZnWgM4RTfDa+C6QyRz6wl0iaqAsWgjeo0cP7NmzB9OnTwcAjB49GqtXr0ZRURGMRiN+/PFHjBw50qpByfqC/T1wo1CDzBwV2rb0FTsO1UOZ11WYv/oISjR6tAv1Q9ykblDIZWLHIiKyiEWl6bnnnsOZM2eg1WqhVCoxbdo05OTkYOfOnZBKpRgxYgTmzJlj7axkZSH+Hki4kMt1TWQTOfmlmPd/h1FUokV4iDfmPxsDZ2XD2NKDiOoni/4ECwwMRGBgoOlrJycnvPvuu3j33XetFoxsj1fQka0UFpdh/urDyC/SoFmABxZO6QlXZ4XYsYiIasWihQVxcXE4derUHZ8/ffp0pU0nqW4KubUYnDNNZE2aMj3eXncU2bklaOzjgref6wlPt3vfPZyIqK6zqDT98MMPuHz58h2fz8rKwrZt2yzNRHZSMdOUc7MUGq1e5DRUH+gNRiz56hguXC6Ah6sCC6f0hJ8X7+xORPWDTS5hycnJgbMzbwJb13m5O8HDVQlBALJ5io5qSRAErNicgL+ScqBUyDB/cg9TMSciqg/MXtO0e/du7Nmzx/T15s2bcfjw4SrHqVQqHD58GFFRUdZJSDYV4u+O85fykZlTjLBgb7HjkAP78pfz2Hs8E1KpBHMmdEWb5rwik4jqF7NLU2pqKnbs2AGg/Ga8p06dwtmzZysdI5FI4Orqim7duvHqOQcR4u+B85fykcV1TVQL2/al4vvfUwAAM2I7oVtkgMiJiIisz+zSNHXqVEydOhUA0KZNG7z77rt46KGHbBaM7MN0O5UcliayzB9/ZWLdj+V/gZo4PBKDujcTORERkW1YtOVAUlKStXOQSEKaVNyDjmuaqOZOJOdg6caTAICH+4Tikf7hIiciIrKdWu00l5mZif379+PKlSsAyvdv6tu3L0JCQqwSjmwv2L9824ErucXQG4yQ8/YWZKaLmTex+Is/YTAK6NspCM8+HMUb8BJRvWZxaVqyZAm++uorGI3GSo9LpVJMnDgRr7/+eq3Dke019naBs1IGjdaAqzdKeLUTmeVKbjEWrj0KjdaATq0a4+WxvAEvEdV/FpWmzz77DF988QUGDx6MZ555BmFhYQDKF4t/8cUX+OKLL+Dv749JkyZZMyvZgEQiQbC/B1IyC5CVo2Jponu6WaTB/NVHUFisRViwF+8nR0QNhkXnYjZv3owBAwZg2bJl6NixI9zd3eHu7o6OHTvio48+Qv/+/bFx40ZrZyUb+XtncK5rorsrUevw1pojuJ5fiqZ+bnhrcg/eHoWIGgyLSlN2djZ69+59x+d79+6N7Oxsi0ORffEKOjKHTm/Au5//iUtXiuDt7oSFz/WEjwc3sSWihsOi0uTn53fXK+iSkpLg68uN7RxFsOkKOpYmqp7BKOC/60/gTOoNuDjJ8NaUHmjayE3sWEREdmV2aTp27Bjy8/MBAEOGDMGWLVuwevVqlJaWmo4pLS3F6tWrsWXLFgwbNsz6ackmQm5dQZeVUwyjURA5DdU1giBgzbYzOHT6CuQyCeZO6o5w7h5PRA2Q2QvBJ0yYgPfffx8PPfQQXnrpJSQmJuLDDz/E8uXL0aRJEwDl95zT6/WIiYnBzJkzbRaarKupnxvkMgnKtAbcKFCjia+r2JGoDvluz0X8cugSJBJg1tgu6NS6idiRiIhEYXZpEoS/ZyBcXFzw5ZdfYvfu3ZX2aerduzf69euHAQMGcL8WByKTSRHY2B2Xr6mQmaNiaSKTvccz8fWviQCAySOj0Cc6SORERETiqdXmloMGDcKgQYOslYVEFNLEo7w0XS9Glzb+YsehOuDUhVws31S+2/fo+8PxcJ8wkRMREYmrRgvBOXtUfwWb1jVxMTgB6VeL8N6X5bt99+4YiEnDI8WOREQkuhrNNL322mt47bXXzDpWIpHg/PnzFoUi+6u4B93layxNDV1eoRoL1xxBqUaPdqF+eGVsZ+72TUSEGpamXr16oUWLFjaKQmKq2KspK0cFQRA4q9hAlWp0WLDmKG4UahDU2B1vPN0dSgV3+yYiAmpYmkaNGoWHHnrIVllIREFN3CGRAKpSHQqLtfD2cBI7EtmZ3mDE4i+PIf1qEbw9nLBgSg94uCrFjkVEVGfwlvYEAHBSyOB/66o57gze8AiCgI+/S0DChVw4KWWY/2wMAvy4eSUR0e1YmsikYmfwLO4M3uBs3JWMPccyIZUArz/VFa1CfMSORERU57A0kUnFuqbLLE0Nyu4/M7BhVzIAYNojHdEtMkDkREREdZPZa5rudq85qh9CmtzaduB6schJyF5OJOfg4+9OAQBiB7bC0J4txA1ERFSHcaaJTEICbt24l2uaGoRLVwqx5MtjMBgF9IsOxvghbcWORERUp7E0kUnFXk15hRqUanQipyFbyr2pxoI1R6Eu06N9WCO89EQn7sVERHQPLE1k4uaigK9n+VYDWTk8RVdfFat1WLj2CPKLNAjx98Dcp7tDIedeTERE98LSRJVUXEGXycXg9ZJOb8TiL/5ExjUVfD3L92Jyd1GIHYuIyCGwNFElFVfQsTTVP4IgYPnmkzidcgMuTjLMf7YHmvi4ih2LiMhhsDRRJX+XJp6eq2++2ZGEP/7KglQqwZwJ3REW7C12JCIih8LSRJWE+JdvO8Ar6OqXnUczsHn3BQDAi492ROc2TURORETkeFiaqJKKK+iu55VAqzOInIas4a+k61j5ffleTI8/0BoPxjQXORERkWNiaaJKvD2c4OaigFEArtwoETsO1VJadiH+/dUxGI0C+ncJxrjBbcSORETksFiaqBKJRGLaGTzzGk/RObLcm2osXHsU6jIDOoQ3wozHoiGRcC8mIiJLsTRRFabF4FzX5LBKbtuLqVmAB+ImdYdCzt/uRES1wT9FqQpuO+DYdHojlnx5zLQX01uTuRcTEZE1sDRRFRWlibuCOx5BEPDxdwlIuJgLZyX3YiIisiaWJqoi+NaapuzcYhiMgshpqCY27krG3uOZkEoleH1CN+7FRERkRSxNVEUTH1coFTLo9EZcz+cVdI5i95+XsWFXMgDg+TEd0LWtv8iJiIjqF5YmqkIqlSC4Ma+gcyQJF3Lw8XcJAIBHB7TCkJ4tRM1DRFQfsTRRtZo3LV/XdO5SvshJ6F7SrxZh8ZfHYDAK6BsdhKeGthU7EhFRvcTSRNXqEdUUAHAgIRtGrmuqs/IK1Vi45ghKNXq0C/XDy09EQyrlXkxERLbA0kTV6trWH67OctwoUCMxnbNNdVGpRoeFa4/iRqEGQY3d8cbT3aGQy8SORURUb7E0UbWUCplptmn/ySyR09A/6Q1G/Pur47h0pQje7k5YMKUHPFyVYsciIqrXWJrojvpFBwMADp2+AoPBKHIaqiAIAlZ9fxonknPgpJThzWdjEODnJnYsIqJ6j6WJ7qhjq0bwcleisFiLUxdviB2Hbvluz0Xsis+AVAK8Nq4LWjfzETsSEVGDwNJEdySTSXFfh0AAwD6eoqsT/vgrE1//mggAeG5Ue8TcOoVKRES2x9JEd9X31im6o2evQqsziJymYUu4kINlm04CAEb1C8Pw3qEiJyIialhYmuiu2rbwRSMvZ5Rq9Ei4mCd2nAYrNasA733xJ/QGAb07BuLpEe3EjkRE1OCwNNFdSaUS9DEtCL8mcpqG6VpeCRasPQp1mQEdwhth1pOduRcTEZEIWJronvpGBwEA/krORZmOV9HZU4GqDPNXH0GBqgwtAz25FxMRkYhYmuiewoK8ENTYDTq9EUlZarHjNBjqMj0WrjuKqzdK0MTXFQum9ISrs0LsWEREDZZDl6ZDhw5h9uzZGDRoECIiIvD222+b9TqtVov3338f48aNQ6dOnRAREYH8fO56fScSicS0IPxsBkuTPegNRiz58hhSMgvg6abE28/1hK+ns9ixiIgaNIcuTQcOHEBSUhK6desGT09Ps1+n0Wjw3XffwcnJCV26dLFhwvqj4hRd6lUNVKVakdPUb4IgYPmmk6bNK+c/G4Ogxu5ixyIiavAcujT961//wi+//ILFixfDw8PD7Nd5enrizz//xGeffYbhw4fbMGH9EdzEAy2aesAoAEfP5Ygdp1778pfz+P2vLEilEsyZ0A0RzX3FjkRERHDw0iSVWh5fIuHVRzV1X4cAAMBhXkVnMz/uT8X3v6cAAGbEdkLXtv4iJyIiogoOXZrIvnpFlf8AT8y4ibxCrm2ytgMns7H2x7MAgAnD2mJQ92YiJyIiotvJxQ5Q3wiCgNLSUrFj2ISbExDSWInMXC32HkvH8F7NxY7ksNRqdaV/n03Lx4ffnoAgAEN6hGBYj6B6++vIFv45nmQ5jqV1cTytx5ZjKQiCWWeg6lRpUqlUyMm593qZkJAQKJVKOySqOZ1Oh8TERLFj2Ez75q7IzNVid/wlhPrwh3ptpaen4+pNLT7/LRd6g4DIZi7o3sKIpKQksaM5pPT0dLEj1BscS+vieFqPrcbSnF5Rp0rTjh07MG/evHset337doSFhdkhUc0pFAqEh4eLHcMm1Go1ijWp2PFXAa7k6+DTpDkC/FzFjuWQ1Go10tPT4eYdgE0/noZWLyCypQ/inoqGUsHNK2uqYjxbtGgBFxcXseM4NI6ldXE8rceWY5mSkmLWcXWqNMXGxiI2NlbsGLUikUjg6lp/i4S7swxRYX44nZKHY0l5ePyBRmJHclglGgM+3XQeBcVatGjqifnP9oSbCzevrA0XF5d6/fvPnjiW1sXxtB5bjKW5F4dxITjV2H3ty6+i23cyG4IgiJzGMWm0BmzYl4ereaVo7OOCBVN6sDAREdVxdWqmqaays7Nx5swZAOXTdpcvX8aOHTsAAEOGDDEdFxkZiVGjRuG9994zPbZv3z6o1WqcPVt+tdLvv/8ONzc3hIeH19vTa9bSPbIx1v4kReZ1FdKvFqFloJfYkRyK3mDE0k2nkZ2nhbuLAgun9ISfF6ftiYjqOocuTfHx8YiLizN9feDAARw4cAAAkJycbHrcYDDAaKx8o9mFCxciOzvb9PXcuXMBANOnT8eMGTNsGdvhuTor0LWtP46cuYr9J7NZmmrAYDDiow0ncPLCDchlErw+vhNC/M3fmJWIiMTj0KVpzJgxGDNmzD2Pu71AVdi7d68tIjUYfaODyktTQjYmDGvLzULNYDAKWLrpJPYnZEMmk+Cx3r5o3cxb7FhERGQmrmkii3Rt6w8XJxly8kuRfPmm2HHqPKNRwMebE/DHX1mQSSV45fEOaB3EU3JERI6EpYks4qyUI6ZdUwDA/pPZ9zi6YTMaBaz8/hR2H7sMqVSC18Z3Rbe2TcSORURENcTSRBbrGx0EADiYkA2DkVfRVUcQBKzedgY7j2ZAKgFmje2M+zoGih2LiIgswNJEFuvUugk8XBW4qSrD2ZQbYsepcwRBwNofz+KXQ5cgkQAvPRGNfp2DxY5FREQWYmkiiynkUvTqUD5rsu9klshp6hZBEPDlL+fx4/40AMD02E4Y0JU34CUicmQsTVQr/aLLZ04On7kKnd4gcpq6Y/2OJHz/e/m2/C880gEPxvDmxkREjo6liWolMtQPvp7OKFHrcCLp3jdbbgi+3ZWMTbsvAACeG9UeQ3u1FDkRERFZA0sT1YpMKkGfTuULwnkVHfDdngvYsDMJAPDsw+3wUJ9QkRMREZG1sDRRrVVcRRd//ho0ZXqR04jnhz9S8NX2RADAhGFtMaofb8dDRFSfsDRRrbUK8UZTPzeUaQ348/w1seOI4scDqfjsp3MAgCcHt0HswNYiJyIiImtjaaJak0gk6BPdcE/RbT98CWu2ld/4+fFBrTH2wQiRExERkS2wNJFVVJyi+yvpOopLtSKnsZ+dRzOw6vvTAIBH+odj3JA2IiciIiJbYWkiq2ge4IkWTT2hNwg4fOaq2HHsYvefl/HJlgQAwMN9QzFxeCRvXExEVI+xNJHV9DWdoqvfG10KgoBvdyZh2aaTEARg+H0tMfnhKBYmIqJ6jqWJrKZi64EzKTdws0gjchrb0OoM+GD9X9iwKxkAMPr+cDw3qj0LExFRA8DSRFYT4OeGiOY+MArAgVP1b0F4gaoMb6w6hP0nsyGTSjA9thOeeagdpFIWJiKihoCliayq4hTd78czYTAKIqexnoyrRZi9fD+SMm7CzUWBhc/1xOAevDUKEVFDwtJEVtWnYxAUcilSsgqx6vtTEATHL05/JV3HaysOICe/FE0bueGDmX3QsVVjsWMREZGdsTSRVfl4OmPWk50hlZRfjv/Fz+cdujj9fDANb689CnWZHlFhfvhgZl8EN/EQOxYREYlALnYAqn96dwxCqUaPFZsTsPWPFLi7Khxuh2yDwYi1/zuLnw9dAgAM7BaCFx/tBIWcf88gImqoWJrIJh6MaY5SjQ7rfjyHr7YnwtVZgeH3tRQ7lllK1Dq8/81xnEjKAQBMHB6JR/qH8wo5IqIGjqWJbGZUv3AUq3XY9NsF/N/W03B1lqN/lxCxY93V9fxSvL3uKC5fU0GpkGH2k53Rq0Og2LGIiKgOYGkimxo3uA1K1Dr8fPASlm48CVcnOWKimoodq1pJ6flY9Hk8Cou18PV0wpvP9EB4iLfYsYiIqI7gAg2yKYlEgikj22NA1xAYjQL+/fVxnLqYK3asKv44kYW5qw6hsFiL0CAv/PelfixMRERUCUsT2ZxUKsHMxzqhR1QAdHojFn0Wj+SMfLFjAQCMRgEbdibhv+v/gk5vREy7ACx5sTcaebuIHY2IiOoYliayC5lMin891RWdWjWGRmvAgjVHkX61SLQ8giDgeOJ1vLJ0H769dUuUMfeHY+6k7nBx4llrIiKqiqWJ7EYhl2Hu093RprkPitU6zP/0MK7cKLZ7jnNpeYhbeQgL1x5FWnYhXJzkmPlYJzzNW6IQEdFd8K/UZFcuTnK8NbkH4lYeQvrVIrz5f4fx7+l97HI6LDWrAN/sSMLxxOsAAIVciuH3tcSjA1rBy93J5p9PRESOjaWJ7M7dVYm3p/bE6x8fxNUbJXjz08NY8mJvmxWX7NxirN+RhAMJ5TcRlkoleKB7MzzxQATXLhERkdlYmkgUPh7OWDS1F17/+ACycorx1pojeHfafXBzUVjtM3JvqrHxt2TsPnYZRqMAiQTo2ykYTw6JQGAjd6t9DhERNQwsTSSaJr6ueHtqL8StPIjUrEK881k8FkzpAWdl7X5ZFqjK8N3eC9h+KB16gxEA0C3SH08NbYuWgV7WiE5ERA0QSxOJKsTfAwun9MTcVYdwLi0PS748hkcHtIKzkxyuTnI4O8nhrJTBWSm/5yLtErUOP+xLwY/7U6EuMwAAosL8MGFoJNq29LXHt0NERPUYSxOJLizYG/Of7YH5q4/gr6Qc/HXrnm//5KSUwcVJDhelHM5O5UXKxan8H4VCir8Sr0NVqgMAhAd74alhkYhu3Zj3jCMiIqtgaaI6oV2oH96aHINvfk1CUUkZ1GUGaLR6aMr0MArlx5RpDSjTGlCAsju+T3ATd4wf2ha92jdlWSIiIqtiaaI6o0N4Y7w/o3GlxwRBQJnOAM2tEqUuK/9HU2aA+lapqngswM8NvToEQsa9loiIyAZYmqhOk0gkcFbKby0O515KREQkHu4ITkRERGQGliYiIiIiM7A0EREREZmBpYmIiIjIDCxNRERERGZgaSIiIiIyA0sTERERkRlYmoiIiIjMwNJEREREZAaWJiIiIiIzsDQRERERmYGliYiIiMgMLE1EREREZpAIgiCIHaK+OHHiBARBgFKpFDuKTQiCAJ1OB4VCAYlEInYch8axtC6Op/VwLK2L42k9thxLrVYLiUSCzp073/U4uVU/tYGr778hJBJJvS2E9saxtC6Op/VwLK2L42k9thxLiURi1s9wzjQRERERmYFrmoiIiIjMwNJEREREZAaWJiIiIiIzsDQRERERmYGliYiIiMgMLE1EREREZmBpIiIiIjIDSxMRERGRGViaiIiIiMzA0kRERERkBpYmIiIiIjOwNBERERGZgaWJ7ikjIwPz58/HyJEjERkZiREjRogdyWH9+uuveP7559G3b1906tQJI0eOxJYtW8D7Ztfcvn37MH78ePTo0QNRUVEYOHAgFi9eDJVKJXY0h1dSUoK+ffsiIiICZ86cETuOw9m6dSsiIiKq/PPBBx+IHc1h/fDDDxg1ahTat2+PmJgYTJ48GRqNxu455Hb/RHI4Fy9exL59+9CxY0cYjUb+gK+FL774AkFBQZgzZw58fHxw+PBhvPnmm7h27RqmT58udjyHUlBQgA4dOuCpp56Ct7c3Ll68iBUrVuDixYv47LPPxI7n0FauXAmDwSB2DIe3du1aeHh4mL729/cXMY3jWrVqFdasWYNp06ahU6dOuHnzJo4cOSLKr1GWJrqnAQMGYNCgQQCAOXPm4OzZsyInclyrVq2Cr6+v6euePXuioKAAn3/+OV544QVIpZz8NdfIkSMrfR0TEwOlUok333wT169f5w8oC6WmpmLDhg14/fXX8dZbb4kdx6G1a9eu0u93qrm0tDR8/PHHWLlyJfr162d6fPDgwaLk4Z/QdE/8QW491f0B2rZtWxQXF6O0tFSERPWLt7c3AECn04kbxIEtWrQITzzxBFq2bCl2FCJs3boVwcHBlQqTmPjTkEhkf/31F/z9/eHu7i52FIdkMBhQVlaGc+fO4ZNPPsGAAQMQHBwsdiyHtGPHDly4cAEvvvii2FHqhREjRqBt27YYOHAgPv30U57ytMCpU6fQunVrrFy5Ej179kRUVBSeeOIJnDp1SpQ8PD1HJKLjx49j+/bteP3118WO4rD69++P69evAwD69OmD//73vyInckxqtRpLlizBK6+8wgJfS40bN8aMGTPQsWNHSCQS7N27F0uXLsX169cxf/58seM5lNzcXJw9exYXLlzAW2+9BRcXF/zf//0fnnnmGezatQt+fn52zcPSRCSSa9eu4ZVXXkFMTAwmTJggdhyHtXr1aqjVaqSkpGDVqlWYNm0aPv/8c8hkMrGjOZRVq1bBz88PjzzyiNhRHF6fPn3Qp08f09e9e/eGk5MTvvzyS0ybNg1NmjQRMZ1jEQQBpaWlWLZsGdq0aQMA6NixIwYMGIBvvvkGL730kl3z8PQckQiKioowZcoUeHt7Y8WKFVw3Vgtt2rRBdHQ0YmNjsXLlSsTHx+O3334TO5ZDyc7OxmeffYaZM2dCpVKhqKjItMautLQUJSUlIid0fEOHDoXBYEBiYqLYURyKp6cnvL29TYUJKF+7GBkZiZSUFLvn4UwTkZ1pNBpMnToVKpUKmzZtqnRJMtVOREQEFAoFLl++LHYUh5KVlQWdTofnnnuuynMTJkxAx44dsXnzZhGSUUMXHh5+x9/PZWVldk7D0kRkV3q9Hi+//DLS0tKwfv16XhZvZadOnYJOp+NC8Bpq27Ytvvrqq0qPJSYmYvHixVi4cCHat28vUrL6Y/v27ZDJZIiMjBQ7ikPp378/tm7disTERLRt2xYAcPPmTZw7dw6TJk2yex6WJrontVqNffv2ASifxi8uLsaOHTsAAN27d+c+JDWwcOFC/P7775gzZw6Ki4uRkJBgei4yMhJKpVK8cA5m+vTpiIqKQkREBJydnZGUlIR169YhIiLCtK8YmcfT0xMxMTHVPteuXTu0a9fOzokc27PPPouYmBhEREQAAPbs2YPNmzdjwoQJaNy4scjpHMugQYPQvn17zJw5E6+88gqcnJywevVqKJVKPPnkk3bPIxG4vTPdQ1ZWFgYOHFjtc1999dUd/7ClqgYMGIDs7Oxqn9uzZw9nSGpg9erV2L59Oy5fvgxBEBAUFIQHHngAzz77LK/+soL4+HhMmDABW7Zs4UxTDS1atAgHDhzAtWvXYDQa0aJFC8TGxuKpp56CRCIRO57Dyc/Px+LFi/H7779Dp9Oha9euiIuLQ3h4uN2zsDQRERERmYGX7BARERGZgaWJiIiIyAwsTURERERmYGkiIiIiMgNLExEREZEZWJqIiIiIzMDSRERERGQGliYiIiIiM7A0EREREZmBpYmILLJixQpEREQgPz9f7CgmFZnuZevWrYiIiEBWVpYdUtWOtcZ5zZo1GDJkCIxGo5WSWe7bb7/F/fffD61WK3YUohphaSKqpw4cOICIiAj89NNP1T4/bdo0dOrUqU78EG3oTpw4gRUrVqCoqMgm719cXIy1a9diypQpkErF/2N/zJgx0Ol02Lhxo9hRiGpE/N89RGQTSUlJAICoqKhqnz937hxatWpVJ36INnQnT57Exx9/bLPStGXLFuj1eowYMcIm719TTk5OGDVqFL744gvw9qfkSPinJVE9lZycDHd3d7Ro0aLKc7m5ucjJyUGbNm3sH4zsbuvWrRgwYACcnJzEjmIydOhQZGdn4+jRo2JHITIbSxNRPZWcnIzIyEhIJJIqz507dw4ArFKabt68iZdeegmdO3dGTEwMFi1ahLKysirHnT9/HpMnT0bnzp0RHR2NiRMnIiEhodIxFet3MjIyMGfOHHTt2hVdunRBXFwc1Gp1pWOPHz+ORx55BO3bt8egQYOscqrn+vXriIuLQ69evRAVFYXhw4djy5YtFmcEgPj4eIwZM6ZSztvXXq1YsQLvv/8+AGDgwIGIiIiodr2VSqUy6/P+KTMzE8nJyejVq1elx7Ozs7FgwQIMHjwYHTp0QExMDGbOnGn2Oq85c+ZgwIABVR43d11ZVFQUvL29sWfPHrM+j6gukIsdgIisT6vV4tKlSxg+fDgyMjKqPF/xt3tzfrjdy8svv4ygoCDMnj0bCQkJ+Prrr1FUVGQqAgBw8eJFjBs3Dm5ubpg8eTLkcjk2bdqEp556Ct988w06duxY5T2Dg4Mxa9YsnD9/Ht999x18fX3x2muvASgvhM8++yx8fX0xY8YM6PV6rFixAn5+fhZ/Hzdu3MBjjz0GiUSCcePGwdfXF/v378cbb7yB4uJiTJo0qUYZgb+LYuPGjTFjxgwYjUZ88skn8PX1NR3zwAMPID09HT///DPi4uLg4+MDAJWOMffzqnPy5EkAQGRkZKXHz5w5g5MnT2L48OEICAhAdnY2vv32W0yYMAG//PILXFxcajyGNRUZGYkTJ07Y/HOIrIWliageSk1NhU6nw7Zt27Bt27Y7HmeN0hQcHIxVq1YBAMaNGwd3d3ds2LABzzzzjGkma+nSpdDpdPj2228REhICABg1ahSGDBmC//znP/jmm28qvWfbtm3x3nvvmb4uKCjAli1bTAVh+fLlEAQB69evR2BgIABg8ODBeOihhyz+Pj766CMYDAb89NNPpuIyduxYzJo1Cx9//DGeeOIJODs7m52xIqdMJsO3334Lf39/AOWnpYYNG2Y6pk2bNoiMjMTPP/+MQYMGITg4uNp85nxeddLS0gCgyvvef//9GDJkSKXH+vfvj8cffxw7d+7EqFGj7vq+1hASEsLSRA6Fp+eI6qHk5GQAwNy5c7Fs2bIq/3h7eyMoKAgeHh61/qxx48ZV+nr8+PEAgP379wMADAYDDh06hEGDBpkKEwA0adIEI0aMwF9//YXi4uJK7/HEE09U+rpr164oKChAcXExDAYDDh48iEGDBpkKEwCEhYWhd+/eFn0PgiBg165dGDBgAARBQH5+vumf3r17Q6VSmU5pmpOx4vs+cuQIBg4caCpMANC8eXP06dOnxhnv9Xl3UlBQALlcDjc3t0qP314AdTodbt68iWbNmsHT0xPnz5+vcT5LeHp6QqPRmHWakagu4EwTUT2UlJQEmUyGsWPHQqlUVnpOo9GgqKgIXbp0AVB+Ku+tt97CkSNHUFRUhPDwcMTFxSE6Otqsz2revHmlr5s1awapVGpaG5Ofnw+1Wo2WLVtWeW1YWBiMRiOuXr2KVq1amR6/vQwB5T9cAaCwsBBqtRoajabK5wJAy5YtsW/fPrNy3y4/Px9FRUXYtGkTNm3adMdjbne3jO7u7sjLy7tjzuoeu5d7fV5NaTQafPrpp9i6dSuuX79e6So2lUpV4/ezRMVnVrfujqguYmkiqoeSk5MREhJSpTAB5afujEaj6dScXq9HUFAQNmzYgICAAPz666+YNm0a9u7dW2V2whzW+AF4p20QbHV5esVeVQ8//DBGjx5d7TH/PJVp74yWfp63tzf0ej2Ki4srlat33nkHW7duxcSJE9GpUyd4eHhAIpHglVdeMet7uNP/Z4PBcM/XVigqKoKLi0ulWS+iuoyliageSk5ORufOnat97uLFiwD+vnLO1dUV06dPNz0/fPhwLF68GJcuXbrjHk+3y8jIqHTaLSMjA0aj0bSGxtfXFy4uLrh06VKV16alpUEqlaJp06Zmf2++vr5wdnaudoF7dZ9h7nu6ubnBaDRWucrMUn5+fnBycqo25z8fs+VMS2hoKAAgKyur0tWSFeuW5syZY3qsrKzM7FkmT0/PaveVunLlitnZsrKyTPmIHAHXNBHVM7m5ucjLy0N4eHi1z6ekpAC483YD6enpKCwsNPsU0vr16yt9XbGou2/fvgAAmUyG++67D3v27Kl0OfuNGzfw888/o0uXLjU6vSSTydC7d2/s3r270g/o1NRUHDx40Oz3+ed7Dh48GDt37sSFCxeqPG/JLUxkMhl69eqFPXv24Pr166bHMzIycODAgUrHVlypZovTYhWnWc+ePVsl3z99/fXXVWaK1Go1UlNTq4xBs2bNoFKpTJuoAkBOTg5+++03s14PlF9deKdyT1QXcaaJqJ6p+CF2+xqh26WkpMDV1RXNmjWr8pxGo8Frr72GqVOnmr1IPCsrC9OmTUOfPn2QkJCAH3/8ESNGjKhUyl5++WUcPnwYTz75JJ588knIZDJs2rQJWq32nld/VWfGjBk4cOAAxo0bh7Fjx8JgMOCbb75BeHi4aRF8Tc2ePRvx8fF47LHHEBsbi/DwcBQWFuLcuXM4cuQI/vzzzxq/5/Tp03Hw4EGMHTsWY8eOhdFoxDfffINWrVohMTHRdFy7du0AlF/BN2zYMCgUCvTv3x+urq4WfS+3CwkJQevWrXHkyBE8+uijpsfvv/9+/O9//4O7uzvCw8ORkJCAw4cPw9vbu9LrT58+jQkTJmD69OmYMWOG6fFhw4bhgw8+wPTp0/HUU09Bo9Hg22+/RcuWLSstmr/T68+ePYuCggIMHDiw1t8jkb1wpomonqkoDXeaabp48SJat25d5ZSQTqfDSy+9hGbNmuHFF180+/OWLl0KpVKJ//73v9i3bx/Gjx9f6dJ4oLzArV+/Hq1atcKnn36KTz75BIGBgfjqq6+q7NFkjjZt2mDdunXw8fHB8uXL8f3332PGjBl44IEHavxeFRo1aoTvvvsOY8aMwW+//YZ33nkHX331FQoLC/Hqq69a9J5RUVFYs2YNvLy8sGzZMmzZsgUzZ85Ez549K+3O3aFDB7z00ktISkpCXFwcZs2aZdUbIT/yyCPYu3cvNBqN6bE33ngDI0eOxE8//YQlS5YgJycHn3/+udnr2Hx8fPDxxx/DxcUF//nPf/DDDz9g1qxZ6N+/v1mv37FjBwIDA9GjRw+LviciMUgE3viHqMEzGo2YPXs21Go1Pv74Y8jlnIS2pRdeeAEpKSnYtWuXXT5PpVJh0KBBePXVVxEbG2uXz7wbrVaLAQMGYMqUKZg4caLYcYjMxpkmIsL8+fORm5uLZcuWsTBZ2e2zO0D5mrH9+/eje/fudsvg4eGBZ599FuvWrTNdKSim77//HnK5HGPHjhU7ClGNcKaJqIHLzs423cz19sXBa9asQdeuXUVMVj/07t0bo0ePRkhICLKzs7Fx40ZotVr88MMP1d5MmYjqLpYmIiIbiouLQ3x8PHJzc6FUKtGpUyfMmjXLtPibiBwHSxMRERGRGbimiYiIiMgMLE1EREREZmBpIiIiIjIDSxMRERGRGViaiIiIiMzA0kRERERkBpYmIiIiIjOwNBERERGZgaWJiIiIyAwsTURERERmYGkiIiIiMsP/A7Ta+oKL0T/OAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "E_n = jax.vmap(nuclear_energy)(H.basis.structure)\n", + "E_total = energy(Z, H) + E_n\n", + "\n", + "ax = sns.lineplot(x=rs, y=E_total)\n", + "ax.set_xlabel(\"$H_2$ bond length (a.u.)\")\n", + "ax.set_ylabel(\"Total Energy (Hartree)\");" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "jax", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}