{ "cells": [ { "cell_type": "markdown", "id": "3ea3ba7b", "metadata": {}, "source": [ "# BAMBOO charges, energy decomposition, and deployment\n", "\n", "A quick tour of what BAMBOO predicts beyond the energy: **per-atom partial\n", "charges** (conserved to the total charge), the **energy split** into the\n", "semi-local NN term and the charge-equilibrium electrostatic term, the molecular\n", "**dipole**, and how to run it as an ASE calculator. No training; this is about\n", "the model's outputs and symmetries." ] }, { "cell_type": "markdown", "id": "e9837591", "metadata": {}, "source": [ "## 0. Build a BAMBOO" ] }, { "cell_type": "code", "execution_count": 1, "id": "220e1161", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:40.580995Z", "iopub.status.busy": "2026-07-20T04:24:40.580781Z", "iopub.status.idle": "2026-07-20T04:24:42.417133Z", "shell.execute_reply": "2026-07-20T04:24:42.416294Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BAMBOO | 103620 params | 64 invariant features\n" ] } ], "source": [ "import logging, warnings\n", "logging.disable(logging.WARNING); warnings.filterwarnings(\"ignore\")\n", "import numpy as np, torch\n", "torch.set_default_dtype(torch.float64)\n", "import xnn\n", "from xnn.common.config import from_dict\n", "from xnn.common.data import structure_to_graph\n", "from xnn.common.models import build_model, ForceStressOutput\n", "\n", "cfg = from_dict({\"model\": {\"name\": \"bamboo\", \"cutoff\": 5.0, \"n_features\": 64,\n", " \"n_rbf\": 32, \"n_interactions\": 3, \"extra\": {\"num_heads\": 16}}})\n", "model = build_model(cfg.model)\n", "print(model.__class__.__name__, \"|\", sum(p.numel() for p in model.parameters()),\n", " \"params |\", model.node_feature_dim, \"invariant features\")" ] }, { "cell_type": "markdown", "id": "ead95b23", "metadata": {}, "source": [ "## 1. Outputs on a small molecule\n", "\n", "A single water-like cluster. BAMBOO returns the per-structure energy, its\n", "`energy_nn` / `energy_elec` components, per-atom `charges`, the `dipole`, and\n", "the invariant `node_features` (consumed e.g. by the LES wrapper)." ] }, { "cell_type": "code", "execution_count": 2, "id": "572ed16d", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:42.420099Z", "iopub.status.busy": "2026-07-20T04:24:42.419932Z", "iopub.status.idle": "2026-07-20T04:24:42.445811Z", "shell.execute_reply": "2026-07-20T04:24:42.444833Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "node_energy shape (6,)\n", "energy shape (1,)\n", "node_features shape (6, 64)\n", "charges shape (6,)\n", "dipole shape (1, 3)\n", "energy_nn shape (1,)\n", "energy_elec shape (1,)\n", "\n", "energy == energy_nn + energy_elec : True\n", "energy == sum(node_energy) : True\n" ] } ], "source": [ "rng = np.random.default_rng(0)\n", "s = {\"pos\": rng.uniform(0, 3, (6, 3)), \"atomic_numbers\": [8, 1, 1, 6, 7, 9]}\n", "g = structure_to_graph(s, 5.0)\n", "out = model(g)\n", "for k, v in out.items():\n", " print(f\"{k:14s} shape {tuple(v.shape)}\")\n", "print(\"\\nenergy == energy_nn + energy_elec :\",\n", " torch.allclose(out[\"energy\"], out[\"energy_nn\"] + out[\"energy_elec\"]))\n", "print(\"energy == sum(node_energy) :\",\n", " torch.allclose(out[\"energy\"], out[\"node_energy\"].sum()))" ] }, { "cell_type": "markdown", "id": "f59920ce", "metadata": {}, "source": [ "## 2. Charge conservation: neutral and charged\n", "\n", "BAMBOO squashes the raw charge to `[-charge_ub, charge_ub]` and then shifts every\n", "atom so the total exactly matches the requested charge (0 by default, or a\n", "`total_charge` attribute on the graph)." ] }, { "cell_type": "code", "execution_count": 3, "id": "370930a0", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:42.447383Z", "iopub.status.busy": "2026-07-20T04:24:42.447201Z", "iopub.status.idle": "2026-07-20T04:24:42.457479Z", "shell.execute_reply": "2026-07-20T04:24:42.456890Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "neutral total charge: 1.3877787807814457e-17\n", "anion total charge : -1.0\n" ] } ], "source": [ "print(\"neutral total charge:\", float(model(g)[\"charges\"].sum()))\n", "g.total_charge = torch.tensor([-1.0]) # e.g. an anionic cluster\n", "print(\"anion total charge :\", float(model(g)[\"charges\"].sum()))" ] }, { "cell_type": "markdown", "id": "cac29356", "metadata": {}, "source": [ "## 3. Symmetries: energy invariant, forces equivariant, dipole invariant in magnitude" ] }, { "cell_type": "code", "execution_count": 4, "id": "fd5d11b7", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:42.459067Z", "iopub.status.busy": "2026-07-20T04:24:42.458894Z", "iopub.status.idle": "2026-07-20T04:24:42.561899Z", "shell.execute_reply": "2026-07-20T04:24:42.561311Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dE (rot+trans) : 1.1102230246251565e-16\n", "dF (rotate forces) : 9.71445146547012e-17\n", "d|dipole| : 1.6653345369377348e-15\n" ] } ], "source": [ "fs = ForceStressOutput(model)\n", "R, _ = np.linalg.qr(rng.standard_normal((3, 3)))\n", "if np.linalg.det(R) < 0: R[:, 0] *= -1\n", "o1 = fs(structure_to_graph(s, 5.0))\n", "s2 = {\"pos\": s[\"pos\"] @ R.T + 4.0, \"atomic_numbers\": s[\"atomic_numbers\"]}\n", "o2 = fs(structure_to_graph(s2, 5.0))\n", "print(\"dE (rot+trans) :\", float((o1[\"energy\"] - o2[\"energy\"]).abs().max()))\n", "print(\"dF (rotate forces) :\", float(np.abs(o1[\"forces\"].detach().numpy() @ R.T\n", " - o2[\"forces\"].detach().numpy()).max()))\n", "print(\"d|dipole| :\", float((o1[\"dipole\"].norm() - o2[\"dipole\"].norm()).abs()))" ] }, { "cell_type": "markdown", "id": "0c4888d0", "metadata": {}, "source": [ "## 4. The electrostatic term is long-range\n", "\n", "Unlike the 5 Å semi-local GET, `energy_elec` sums a damped Coulomb over **all**\n", "pairs. We slide two fluoride-like atoms apart (freezing the predicted charges)\n", "and watch the electrostatic energy follow the expected ``~ q_i q_j / r`` tail\n", "well beyond the cutoff." ] }, { "cell_type": "code", "execution_count": 5, "id": "20c11fe3", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:42.563708Z", "iopub.status.busy": "2026-07-20T04:24:42.563635Z", "iopub.status.idle": "2026-07-20T04:24:43.286132Z", "shell.execute_reply": "2026-07-20T04:24:43.285498Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAekAAAE2CAYAAABFpmptAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZipJREFUeJzt3XdcU9f7B/BPAiHsDTJkiLgHbrHuvavFXbditbVua8VarXV+q3XUtrbuuidSV+teCIoDVFQcCLJBVtgEkvP7gx+3xgQIIZCAz/v1ykty77nnPgmRJ+feM3iMMQZCCCGEaB2+pgMghBBCiGKUpAkhhBAtRUmaEEII0VKUpAkhhBAtRUmaEEII0VKUpAkhhBAtRUmaEEII0VKUpAkhhBAtRUmaEEII0VKUpEmNsXDhQujr62s6jGrNzs4O3t7emg5Da+jr62Pu3LmaDqPStWnTBr169apQHR/Le6UIYwwrVqyAq6srdHR00KZNGwBASkoKPv/8c9ja2oLH42Hp0qXlrpuSdCUKCQkBj8eTeZiYmKBp06ZYs2YNcnNzSzz24cOH4PF40NXVRUxMTJn1r1ixQmGZb7/9livz+vXrCscWGBiIESNGwMHBAXp6eqhVqxaGDBmCq1evlvhaVDlGG3399dcwNjbWaAzm5uaYMWOGxuvQRjX1dWmTmvgeq+M1+fr64ocffsCOHTtQUFCA+/fvAwCWLl2KgIAA3L9/H4wxrFq1qtx1U5KuAitXrgRjDIwxREdHY9GiRVi+fDnGjx9f4jE7duyAlZUVBAIBdu/eXWr9RkZG+Ouvv/DhNOwSiQQHDhyAkZGRWmLbsGEDOnbsCEtLS1y9ehVZWVkIDAyEu7s7evXqhWXLlqnlGKI5CQkJ2Llzp6bDIKRauXbtGmxtbdG7d2/w+XyZ7V27doWzs7PKdVOSrmLm5uaYMGECBg0ahFOnTilssebk5ODQoUPw9vbGyJEjsXv3bkil0hLr9PLyQkREBG7evCmz/cKFC4iLi8OwYcMqHNvFixexaNEizJkzB3/++ScaNmwIPT09uLm54eeff8by5cuxcuVKHD9+vELHEEJIdZOUlAQDAwOlt5cLI5UmODiYAWArV66U2zd06FAmEAhYQUGB3L7du3czPp/PIiIi2J07dxgA9u+//5ZY/6ZNm1jz5s3ZpEmTZPYPHz6cderUia1cuZIBYK9evVI5to4dOzJjY2OWkZGh8LXm5+czGxsb1qRJkwodU5LHjx8zLy8vZmVlxfT09FjDhg3Z1q1bZcosWLCACYVClY4tLjdixAhma2vLDAwMWKtWrdi+ffuYVCplAwcOZADkHtHR0YwxxqZOncqsrKxYamoqGz9+PLOwsGBNmzZljDFWUFDAVq5cyerXr8/09PSYtbU1GzNmDHvz5o3M+Y8fP87at2/PzMzMmKWlJevevTv3e09LS1N4/q5du3LHe3h4cNt1dHSYg4MDmzhxIouLi1O6jlq1arGpU6fKvTd//fUXa9++PTMyMmJWVlZsyJAh7MmTJ0rFXpqyYlYmBmVel6rxCYVCNmfOHHb16lXWsmVLJhQKmbu7O/v9998Vlt+zZw9r06YNMzAwYMbGxqx3797s3r173P6WLVsyDw8Phce2a9eONWrUiHuu7OemOMbAwEDWvn17JhQKWZ06ddj27dvlzpGTk8Pmzp3LatWqxYyMjFjfvn1ZREQEa926NevZs2eJ74My73F54lDmvSpNaZ/Hd+/eMQBs/fr1csd5eHiwvn37Kv2ayvodREREKKyjpEfx34vyoCRdiRQlQpFIxPbv388EAgH7+uuvFR73ySefsIEDB3LPW7duzYYNG1Zi/Zs2bWIbN25kxsbGLCsrizHGWEpKCtPT02M7d+5UOkmXFJtIJGJ8Pp/179+/1Nc7evRoBoDFxMSodExJAgICmIGBAfvss8/Y8+fPWUZGBjt69CgzMTFh3333HVdOUZJW9thbt24xfX191q9fPxYcHMwyMzPZw4cP2fjx41l4eDhjjLGZM2cyIyMjhTFOnTqVWVpashEjRrBz586x1NRUtnv3bsYYYyNGjGCGhoZs3759LD09nT18+JC1bNmS2dracv9pb968yXg8Hlu3bh1LTk5mIpGI3bx5k/Xv358VFhZy5zEzM2PTp08v9T1lrOiPcUBAAGvRogVr166d0nUoStJz585lAoGArVmzhkVFRbHU1FT2999/M29v73LFXpGYy4qhtNdVkfiEQiHr168fGzp0KHv58iV79+4dW7NmDePxeGzjxo0yZRcuXMiEQiHbsmULS0xMZLGxseyrr75i+vr67OHDh4wxxnbs2MEAsNu3b8sc++DBAwaAbdiwgdumzOemOMaBAweyUaNGsRcvXrCUlBQ2f/58BoAFBgbKnGfIkCHMzMyMnTx5kolEInb79m3Wr18/1qRJk1KTdFnvcXnjUOa9KklZnwVlk7Qyr0nZ38GwYcOYi4tLud4vZVGSrkTFiVDRY/To0Sw3N1fumKdPnzIA7Ny5c9y2nTt3MoFAwBITExXWv2nTJpaUlMQEAgHbu3cvY4yxX375hRkaGrKMjIxSk7QysRXHNGPGjFJf77fffsv9AVLlmJK0bNmSNWzYkInFYpnt69atY3p6etz7oihJK3tss2bNWN26deXKva+sJA2AHTx4UGZ7QEAAA8DWrFkjsz0iIoLp6upy/4FXrVrF+Hy+wisr7yvvf/rr168zAOz+/ftK1fFhkg4JCWEA2JIlS0o8h7KxqxqzMjEwVvLrqkh8QqGQmZuby10NGj58ODM1NeW+FIeGhjIej8eWL18uU04qlbLmzZtzX7qzs7OZubk5GzdunEy5qVOnyvwfV/ZzUxyjjY0NFwtjRVeprKys2OTJk7ltt2/fZgDkriL5+/szAGpJ0srEoex7pYgynwV1Jeny/A4qM0nTPekq8H7nrKysLFy+fBmBgYHo0aMHxGKxTNkdO3agTp066NevH7ft888/5zqHlcTGxgb9+/fH3r17AQB79uzB8OHDYWJiUuHY2Acd0kpSXI7H46l0jCKxsbEIDg7G0KFDIRAIZPb16tULYrEYgYGBFTo2OjoaT548wYgRI+TKldenn34q8/zKlSsAivoNvM/V1RWtW7fm9nt4eEAqlWLMmDG4efOm3OdCGY8fP8awYcNgZ2cHXV1d8Hg8dOvWDQBkevaXx/nz5wEUfQZLUpHYlYlZmRhKU9H3tnv37nL/j4YOHYqMjAyuF+/Zs2fBGMOIESNkyvF4PPTo0QM3btwAABgaGmLSpEk4fvw4kpOTAQDp6ek4fPgwBg8eDFtbWwDKf26KdevWTaaDqJ6eHurXr483b95w24qP+fAz2rFjR1hbW5fjHSmZMnEo+14pUtHPQnmU93dQWShJVzEjIyP07NkTq1evRmBgIA4cOMDtE4vF2L9/PyIiIqCjo8MNjTI0NER6enqZvW4nTZqEGzdu4O+//0ZwcDAmTZqklticnJzA5/Px9u3bUo+PioriyqtyjCIJCQkAgJ9++gm6urrQ0dGBjo4O+Hy+zFjEihyblJQEAHB0dCw11rKYmprKDdEqjs3Ozk6uvJ2dHfeHetCgQdi+fTueP3+Orl27wszMDH369MG///6r1Lmjo6PRqVMnZGVl4cKFC8jMzARjjEsiBQUFKr0mZd4bVWNXNuaK/n4q+t7WqlWrxG3Fv7/iz5qHhwf3WePz+eDz+di8eTOysrK4LwdfffUVxGIxN2rjr7/+Qk5ODqZOncrVr+znppi9vb1cOVNTU6Snp8vVWdrrqShl4ijPe/Whin4WlG08AOX/HVQWStIaUr9+fQBAaGgot83X1xcZGRnIyMjgWrfFj+joaLx8+bLUb5mDBg2ClZUVJk+eDFdXV65FUtHYTE1N0b59e9y6dQtZWVkKjxGLxbhy5QqaNGmC2rVrq3SMIsXf8H/88UcUFhZCIpFAIpFAKpVy782UKVMqdKyNjQ2AopZ3RShqhVtaWgIAEhMT5fYlJibKtGCmTZuG0NBQJCYmYv/+/cjNzcWAAQNw69atMs996tQpZGZmYtu2bfDw8OB6lEZERKj6cgBA6fdGldiVjVkdv5+KvLcl/e4AwMrKCsB/n7Xw8HDusyaVSmU+a3p6egCAevXqoVevXvjjjz8gkUjwxx9/wNHREX379uXqL8/nBij5StT7imMt7fVUlDJxlOe9+pAynwVTU1PweDxkZmbK7SvPZ6i8v4PKQklaQ16+fAlA9lvajh070LFjR4WXqGvXro2mTZuW2poWCAQYO3Ys0tLSMHHiRKX+wygb2/fff4+srCx8//33Co9Zu3Yt3r17J7NflWM+5OLigmbNmsHX1xeFhYXleh3KHuvs7IxmzZrh+PHjpbY4jYyMyn2ptGfPngCKEtL73r59iwcPHnD732dra4vhw4dj//79YIzJJBIjIyPk5+eXeD6hUCjzfN++fQpfR2l1vG/gwIEAgEOHDilVvrTYS1JWzMrGoMzrUiW+69evy33R/Pvvv2FiYoK2bdsCKPqCDABHjx4tsz6gqDUdERGBRYsWISwsDJMmTYKOjg63X5XPTVl69OgBADhz5ozM9oCAAKVbheX57JSkvO/V+5T5LOjp6aF27doyDSAAuHPnjsKrbiW9psr4HaikQne0SakU9aDOyspiV65cYS4uLszW1pYlJCQwxhgLDw/neqCW5JtvvmH6+vosNTVVpv5NmzaVGoeyvbtLiq3YunXrGI/HY9OnT2dhYWEsPz+fvXnzhi1YsIDxeDzm4+Mjd25VjvlQQEAAMzQ0ZEOHDmUPHz5k2dnZLCoqivn6+rJu3bpxHVVK6t2tzLHFvbv79+/PQkJCWFZWFgsODmYTJkzgenfv27ePAWDnz5+X6xVcPARLES8vL2ZsbMwOHjzI0tPTWXBwMGvTpg2ztrZmUVFRjDHGli5dypYtW8YePXrEsrOzWWJiIvvmm28Yj8dj/v7+XF19+vRhjRo1khuiFBYWxvT09Njo0aO53rKLFi1iw4cPZwDY/v37y6yDMcW9u+fMmcMEAgFbu3Yti46OZqmpqez06dNs2rRp5Yr9Q+WJuawYSntdqsbHWFFnqL59+zIvLy/26tUrlpyczP73v/8xPp8v1zHpm2++YUKhkP3000/s7du3LCcnh4WGhrL169ezmTNnypQtLCxkzs7ODADj8XjcZ+x9ynxuimOcM2eO3PF9+/aVG+41aNAgZmFhwU6dOsUyMjJYQEAA69+/v9K9u0v77JQnjvK8Vx9S5rPw448/MoFAwI4dO8YyMzPZzZs3mZeXF2vUqJFcx7HSXpOyvwPq3V1NKepBra+vz+rVq8e++uormS78Pj4+DAALCQkpsb6rV68yAOyXX36Rqb8iSVqZ2N5369YtNmzYMFarVi0mEAiYjY0NGzx4MLt48WKJ51flmA89f/6cjRs3jjk4ODCBQMBcXFzY8OHD2fXr17kyJY2TVuZYxop6jnp5eTFLS0tmaGjIWrduzY2TZqxozOTkyZOZlZUV4/F4CsdJKyIWi9mKFSuYu7s7EwgEzMrKio0aNYq9fv2aK/Pu3Tu2atUq1qJFC2ZoaMhsbGxY79695d6jJ0+esA4dOjADAwO5MZ3nz59nLVu2ZAYGBszZ2ZmtXr2aPX78WC7hlVZHSeOkd+3axVq3bs0MDAyYjY0N++yzz1hoaGi5YldE2ZjLiqG011WR+IoTz6VLl5iHhwcTCoWsbt267Ndff1VY/siRI6xLly7MxMSEGRoasqZNm7LFixcrHGK4evVqBoB1795dYV3KfG7ej/FDipJjTk4Omz17NrOxsWGGhoasd+/eLDw8vMxx0sVK++yUJw7Gyvdefaisz0J+fj77+uuvmZWVFTMyMmJDhgxh8fHxCnt3l/aalP0dVGaS5jFWjjvphBBC1OLXX3/FrFmzcODAAYwdO1bT4RAtRUmaEEI0oHPnznj27BliY2Np9TZSIuo4RgghVYgxhtOnT8Pf3x/z58+nBE1KRS1pQgipIqGhoWjWrBlMTU0xfPhw/PHHHxWeQIfUbJSkCSGEEC1Fl7sJIYQQLUVJmhBCCNFSupoOoKaQSqWIi4uDiYmJyjN9EUII+TgwxpCZmQkHBwfw+SW3lylJq0lcXFyJi0QQQgghikRHR5e4dgFASVptiufbjo6OhqmpqYajIeUlFovx888/AwAWLFhQ4gT/hBCiDhkZGXBycipzOWHq3a0mGRkZMDMzg0gkoiRdDTHG8O7dOwBFK+3QLQtCSGVSNmdQS5oQFC2xZ2trq+kwCCFEhtJJeu7cueWqePPmzeUMhRBCCCHvU/pyN4/H49byLMu5c+fwsV1Fp8vd1ZtEIuHWFu7cubPM2r6EEKJulXK5++zZs0qVo/t5pLqRSCS4ceMGAOCTTz6hJE0I0QpKT2Zy+PBhpSstT1lCCCGkuogX5SIgPBnxotwqOZ/SLenRo0crXWl5yhJCCCHVwdF7UfDxfQIpA/g8YK1XM4xq61yp56RpQQkhhHxUytMaziuQ4HVSFk4Fx2LxyaIEDQBSBizxDa30FrXSLWlXV1elK42MjFQhFEIIIaRyfdgaXj64Mdq7WSEmNRex6UWPmLQcxKYV/ZycJS6xLgljiEzOgb2ZQaXFq3SSXrp0aaUFQQghhJRHvCgXEcnZqGNtVGqSzCuQIDo1B9FpOQiNzcCmSy9RPPZIyoDlp5+VeS5joS5sTYR4k5wts12Hx4OrtWFFXkaZlE7S3t7elRkHIYQQopQPW8OL+jVECydzRKXmICY1B1GpOYhOy0VUag7eZeaXWZ+RUAd1rI3gaG4AR3ND1LYwgKOFARzNDeBkYQhTA13weDwcvReFJb6hkDAGHR4Pa7yaVmorGlDDtKD5+flgjEFfX19dMVVLNE66epNKpYiPjwcA2Nvbl7oqDSFEPcrTGo5Jy0Fkcg4ex6Tjl6uvy3UeE6EunCwNYWOih5svk/F+0uPzgNuLeyidbONFuYhMzoGrtWGFEnSlTwu6Y8cO/O9//8ObN28AAG5ubvj2228xbdo0VaskRGP4fD4cHR01HQYhH40PW8MrPm2KtnUsEJmcg7cp2YhMKfr3bUoO4kS5KKs5WctUiPq1TOBkaQgnC0M4WxrCydIAzpaGMDMQcPN3VLQ1bG9mUOmt5/eplKTXr1+PlStXYvbs2fD09ASPx0NgYCAWLFiA9PR0fPPNN+qOkxBCSDUnLpQiKjUb9yPT4OP7RObe8Pd/h5Z6rLFQFy5Whqhlqo9rYUlyrWG/mR2VSp6j2jqjS30btbSGq4JKl7tdXFzw+++/y00Teu7cOcycOfOj7N1Nl7urN4lEgjt37gAAPD09acYxQpT04SVrxhiSMvMR/i4Lb95lIyI5G2/eZeFNcjaiU3O4IUyKGAl14G5rAhdLQ7haGcLFygiu1kX/WhnpldoaruzxyupWqZe74+Pj0blzZ7ntnTp1QlxcnCpVEqJREokEly9fBgC0bduWkjQhZcgrkGDb9XD8cuUV16p1MNeHKKcA2WJJiccZ6enAydIQYQmZMtv5PODy/K41sjVcESol6bp168LX1xeTJk2S2X7ixAnUrVtXHXERQgipIqV14MrOL8TrpCy8TsrCq6QsvE7KxKukLLxNyZGrJy49DwCgw+fBycIAbjbGcLM2Qh0bI7hZG6OujRFsTIRq6Sld1feGNUWlJP39999j4sSJ+Pfff9GuXTsAwN27d+Hr64t9+/apNUBCCCGV5/0OXDwe4NXSEeaGelxijk0v34xaG0Z44FMPB+jplj5C4mNqDVeESkn6888/h6OjI9avX4+tW7eCx+OhcePGuHLlCrp06aLuGAkhhKhJ8TSXYQmZeBCZisP3orl9jAEnH8bKHWNtLEQ9W2O42xqjXq2if031Bfj0V3+Ze8w6PB46uluVmaCLfSyt4YpQeQhW165d0bVrV3XGQgghREUfXrKWShmi03IQlpCJFwmZCEvIQFhCJiKTs0vtvAUA/ZvaoUt9G7jbGsPdxhgWRnoKy631alblk3t8bFRO0oQQQrTDntsR+PHsM24scW0LA6Rmi5FTQgcuC0MBGtiZwMnCECcexMgMZ9Lh8bBscGPqwKUlVErSubm5+P333+Hv74+0tDS5/devX69oXIQQQj7AGENiRj6exonwNC4DT+NEeBQjQoIoT6ZcTFrRfWQ9XT7q2RqjgZ0JGtqZoKGdKRramXCdtwCgjasFdeDSYiol6S+++ALXrl3D0KFD0aRJE3XHREiV09XVxcSJE7mfCdGE9y9Z1zLRR2RK9v8n46KE/CwuAynZJa/K9L4NIzwwtIUDdHWoA1d1ptJfo7///hv3799H/fr11R1PmRITE3Ho0CEkJiaiWbNmGDVqVJl/VDMzM3Hw4EGEhYXh66+/hru7u1rqJTUHn88v13KshKiTRMqw9eorbLn835hjPR0exBL5m8c6fB7q2hihqYMZGjuYws5UH7OPBCvswFVWgi5GrWHtpVIWMjQ0hLW1tbpjKdOrV6/QsWNHtGzZEu3atcPy5cuxZ88eXLhwocTJJ3bu3Illy5ahc+fOOHbsGIYOHSqXpFWplxBCVMEYQ3RqLh7FpONRdDoex4jwODYdeQVSmXJiCYOeDg+NHMzQxMH0/x9maGhnAn2B7N+lbHEhdeCqoVSaFvTbb79FXl4eNm7cWKVJzMvLC8nJybh+/Tr4fD6ioqLg7u6OXbt2Yfz48QqPefToEdzc3CASieDk5IRr166hW7duFa73QzQtaPUmkUjw4MEDAEDr1q3pyxlR2Ye9rN9l5uPx/yfkRzEiPI5JR1pOgVJ1HZjaDp3q2Sh9XrpkXX1U6rSgs2bNQvPmzXHw4EG4uLhwHRCK3b9/X5VqS1VQUIDz589j8+bN3DKCzs7O6NatG/z8/EpMph4eHgAAkUik1npJzSKRSPDPP/8AAFq0aEFJmqhkX0Aklp95yvWyNjcQID1XPiELdHhobG+K5rXN0by2GRzNDTBu1125S9Z1bY2VPjddsq6ZVErSkydPhoWFBUaMGAFzc3M1h6RYVFQU8vPz5aYdrVu3Lm7fvl3l9ebn5yM//7/FxDMyMlSOgRBSPcWm5+LB2zQ8iEzFnYgUvEjIktlfnKDr2RqjeW1zeDiZwaO2ORram0CoK/tFkMYcE0VUStL+/v54+vQp3NzcKnTyQ4cOISgoqNQyS5cuhbW1NXJyiuaJNTExkdlvamrK7VOFqvWuXbsWK1asUPm8hBDtVNI81gUSKZ7HZ+DB2zTcf5uGh2/TEP/B0CdFdk9sgx6NapVZjnpZE0VUStJ2dnZque9qbW1dZo9agUAA4L8kmp6eLrM/LS1NLsGWh6r1+vj4YP78+dzzjIwMODk5qRwHIUTz3p/Hms8DJnZwhZFQF/ffpuJRtAi5BbKTg+jweWjiYIpWzhZwszHCD6efyl2ybuSg/N9KumRNPqRSkh4xYgS+++47bN26FXp6iqeLU0afPn3Qp08fpco6OzvD2NgYYWFh6NevH7c9LCwMjRs3VjkGVesVCoUQCoUqn5cQol1CY0VYfPIJNwRKyoA9AZEyZUz1ddHaxeL/H5bwcDKDod5/f0aFuny6ZE3USqUkfebMGYSFheHw4cNwcnKS6zgWGhqqluDex+fzMWzYMOzduxczZsyAvr4+Hj9+jNu3b8PX15crd/ToUURGRuLbb79Va72EkJolOSsfd9+kIvBNMu68ScXrpCyF5brUt0b/pvZo42KBujbG4PN5CssBdMmaqJ9KSXr69OnqjkMp69atQ9euXdG2bVu0bNkS58+fx7hx4zBkyBCuzKVLl3Dnzh0uSQcFBeHQoUPIyir6D/jrr7/Cz88P/fr141rOytRLCKk+FN1XLk7Kd96k4M6bFLwqISm/T4fHw/+GNadpMonGqDROWpNyc3Pxzz//cDODderUSWb/5cuXkZCQgHHjxgEAnj17hosXL8rV4+npCU9PT6XrLQuNk67epFIpXr9+DQBwd3fnhuOR6kdmfWQAnm6WSMkW42WifFJuaGcCTzcreLpZoX0dS1x8liB3uXpUW+eqfxGkxlM2Z5QrSdetWxeDBg3CoEGD0LVr1wrdj65pKEkToln5hRJcfJqA2YdDUNIftQ+TsqIlGGlSEFIVKmUyk61bt+LMmTOYOnUq0tPT0bt3bwwaNAgDBgxArVplDzEghBB1YYwhIjkbN1++w81XybjzJqXEpRnn9aqH8R1cYVnCusjvo8vVRJuofLn70aNHOHv2LM6ePYv79++jVatWXCu7ZcuW6o5T61FLunqTSCR48uQJAKBZs2Y045gWUHRfOSOvAAGvU3Dz1TvcfPmOW5KxmIWhHtJyZFeJ0uHx4L+4OyVeolUqdVpQoGi6TQ8PD3z33Xd49+4dzp8/j7Nnz2LDhg0wMTHBwIED8eeff6paPSFVSiKR4O+//wYANG7cmJK0hsncV+YBvRrVQnqOGA+j0iF5byCyQIeHNi6W6FLfBl3qW6ORnSmOP4imYVCkxlDLWow2NjaYOHEiJk6ciIKCAty8eRNnz55VR9WEkI/M66QsmfHKjAGXniVy+92sjbik3L6OFYyEsn/GaBgUqUnUvmCyQCBAz5490bNnT3VXTQipoeJFubj8PAmXniUi4HWywo5fUzu5YtIndeBkaVhmfXRfmdQUSifpsqbvfF9kZKQKoRBCPhaMMTyPz8Tl54m49CwRT2IVr1JXTIfHg3dnN0q85KOjdJJeunRpZcZBCKmB3u/8ZW0sRFBEKi49K0rMsen/dfri8YCWTubo3dgOvRvXwv3IVHx3iu4rE6J0kvb29q7MOAghNcz7nb+Aonmt8wul3H6hLh+d61mjd+Na6NGwFmxM/psL393WGF0b0H1lQtR+T5oQ8nHLFUtw8mEMlvrJzuGfXyiFuYEAfZrUQq9GtdC5ng0M9EruRU/3lQmpQJJ+9uwZ/Pz8EBUVhcLCQpl9O3furHBghFQlXV1dDB8+nPuZlE9+oQQ3Xybj7OM4XHqWWOKkIr+NbYWO7tZVHB0h1ZdKf41Onz6NUaNGoWvXrrhw4QKGDBmChw8fIjo6GgMGDFB3jIRUOj6fjyZNmmg6jGqlQCJFQHgKzj6Kw79PE5CZ99+XdTtTfSRm5Mn00tbh8eBmY1T1gRJSjamUpJcvX45du3bh888/B4/Hg5+fHwoKCvDVV19BKpWWXQEhpNp4v/OXrYk+giJScfZxHP4JTUBq9n+ze9UyFWJgMwcM9rBHCydzHLtPk4oQUlEqTQtqYGCA5ORkGBkZQSAQICMjAwYGBkhMTETTpk3x7t27yohVq9G0oNWbVCrF8+fPAQCNGjWiVbD+34crShkLdZGZ/1+L2cpID/2b2WFwcwe0dbWUW2uZFqsgRLFKnRY0Ly8PRkZFl63s7e3x6tUrNG/eHBKJBLm5uWUcTYj2KSwsxIkTJwAAPj4+tMIbgMfR6bIzfwHIzC+EsVAHA5rZY7CHAzq4WUFXp+QvNNT5i5CKqXAPmcGDB2Pq1KkYO3YsfH190blzZ3XERQjRAHGhFFeeJ+LY/Whcf/FO4cxf28a2Ruf6NlUeGyEfI5WS9JkzZ7if161bh0WLFmHv3r1o2LAhfv75Z7UFRwipGi8SMnHsfjROBcfK3Gf+kA6PB/daxlUYGSEfN5WS9KBBg7ifTUxMsG3bNrUFRAipGhl5BTjzKA7H7sfgUXQ6t93WRIjhrWtjRBsnBEWkUOcvQjRIpSTNGEN0dDScnZ1ltkdFRcHJyQk8Hq+EIwkhmlDcQ9vVyhBvU3Jx/H40zofGI6+gaDSGLp+HXo1qYWTb2uhSz4a7z1zn/1ecos5fhGiGSkn6l19+QUREBDZv3iyzfePGjahbty5mzZqljtgIIWrw4fSc76tna4xRbZ0wtKUjrI2F8gVAnb8I0SSVkvTmzZtx/fp1ue1z585F7969KUkToiUuP0/EtyefyG0f0sIekz6pgxZO5nTlixAtplKSTkhIgKGh/JquhoaGiImJqXBQhFQ1HR0dDBkyhPu5OiuQSPFPaAL23I5AcFS6wjKj27qgpbNF1QZGCCk3lZJ0q1atsHPnTvj4+Mhs37FjB1q0aKGOuAipUjo6OtX+s5uWLcahoCjsD3yLhIw8AICAz0PBB9e5dXg8uFrLf8kmhGgflZL0ihUr0L9/fwQGBqJLly5gjOHmzZv4999/8c8//6g7RkJIKV4kZGLP7QicCo7lloK0NtbD2PYuGOvpjGthSdRDm5BqSqVpQQHg1q1bWLduHR4+fAgej4eWLVvCx8cHnTp1UneM1QJNC1q9SaVSvH79GgDg7u6uldOCvj+Hdi0TfVwNS8KegAjcfp3ClWniYIopHetgkIc9hLo6MsdSD21CtEelTgtaWFiIzp07K5xdLD09Hebm5qpUS4jGFBYW4vDhwwC0c1rQD+fQtjTSQ8r/TzrC5wF9m9hhcsc6aOtqobAjGPXQJqR6Uqm54O3tDUUN8PT0dPTu3bvCQZXm8uXLGDJkCDw9PTFt2jRER0eXecz169cxevRotGjRAvfv35fbv3r1arRo0ULmUby2MCGaFi/KlRlCxQCkZIthLNTBF13ccHNRd2wb1xrt6lhST21CahiVWtJ37tzBggULsHHjRm5bWloaevfuDTMzM7UF96FLly5hwIABWL58OTp06IAtW7agY8eOePLkSYnnXbx4MQIDA/HZZ5/h6NGjyMrKkisTHR0NGxsbrF+/nttmYECtDqJ5SRl5+PHMM4VjnH/9vBW6NbCt+qAIIVVGpSR98eJFdOzYEdbW1liyZAnS0tLQq1cvWFhYyMzrrW7Lli3DqFGjsHTpUgBAx44dYW9vjz///BOLFi1SeMz3338PIyMjxMTEYN68eSXWbWZmVu1795KaIzEjD9uuh+NwUBTXGex9OjweGtiZaCAyQkhVUilJOzs74+LFi+jcuTMEAgGOHDnCJejKaoFmZWXh7t27MhOl6Ovro2fPnrhy5UqJSbp4Sc2y3LlzB5988gnMzMzQuXNnzJs3j1rTpMrFi3Lxx/VwHL4XDfH/J+dWzuZo4WyBv25HQMJAPbQJ+YiovFRlo0aNcO7cOfTs2RMdOnTA6dOnKzWpxcbGgjEGe3t7me0ODg64fPlyheo2MTHB9OnT0bVrV8TGxmLFihXw8/NDQEAAdHUVv0X5+fnIz8/nnmdkZFQoBvJxi0vPxbbr4Th6LxpiSVFybutqgTk966OjuxV4PB6mda5DPbQJ+cgonaSbNm2qcLtAIEBkZCTatm3LbQsNDVWqzmXLluH06dOlljlz5gycnJxQUFAAABAKZecXFgqF3D5VrVmzBgKBgHvu6emJevXq4ejRoxg7dqzCY9auXYsVK1ZU6LyExKTl4Pfr4Th+PxoFkqIbz+3qWGJuz3roUNdKpiMY9dAm5OOjdJL29vZW+8m9vb3h5eVVaplatWoBAKysrAAAKSkpMvtTUlK4fap6P0EDQJ06deDq6oonT+TnPC7m4+OD+fPnc88zMjLg5ORUoTiI5ujo6KB///7cz5WleKyzUJePEw9icOJBDJecPd0sMadnfXSoW7HPMyGk5lA6Sc+dO1ftJ3d2dpZb7rIk9vb2cHBwQFBQEAYPHsxtv3PnDrp3767WuAoKCpCUlAQTk5I75giFQrlWPam+dHR00K5du0o9R0mrUXV0t8LsHvXQ3o2SMyFEltLjpP/44w+lKy1P2fKYNm0adu7cicjISADAoUOHEBYWhqlTp3Jlfvzxx3KNcRaLxVixYgVycnK45/PmzUN+fj5GjBih1vjJx+tVYiYWn5RP0H+Oa4WD3p6UoAkhCimdpL/88kulKy1P2fL47rvv0Lt3bzRo0ADOzs6YMWMGduzYgVatWnFloqKiEBYWxj0/ffo0WrRogb59+wIousTeokUL7ouEQCCArq4uXFxcUK9ePVhbW+PatWv4559/UL9+/Up5HUT7SKVSREZGIjIyElKp/JAnVUmkDIeDojBsWwAUzb9raqBdM5sRQrSL0nN383g8dOzYUalKb9++rXBGMnVJTk7Gu3fv4OrqKtejPDo6Gjk5OWjQoAEAIDU1FVFRUXJ12NnZwc7OjnsulUoREREBCwsLWFpaljsmmru7ehOLxVi7di0A9U0LGvA6GT+efYawhEyF+3V4PPgv7k6dwQj5CKl97u7ly5crffJevXopXVYV1tbWsLa2Vrjvw85blpaWSiVdPp+PunXrqiU+8nGLSM7GmvPPcelZIgDAzECAub3qQajLx/d+T2k1KkKI0pRO0j/88EMlhkFI9SfKLcDWK6/wV2AkCiQMOnwexnu6YE7PerAwKmqZd29oS2OdCSFKU3kyE0JIkUKJFIfvRWPTpZdI/f+Vqbo1sMHSgY3gbis7QoDGOhNCyoOSNCHl9P66zq8Ss7Dq3DO8TCxauMXd1hhLBzaihS8IIWpBSZqQcihprLO5oQDze9fHmHbOEOiotAIsIYTIUSlJ5+bm0uIT5KPz4brOxUa1ccKSAY1gZihQfCAhhKhIpa/89vb2mDFjBoKCgtQdDyEaoaOjg169eqFXr14lTgt6/km8wnWdh7Z0pARNCKkUKiXprVu34sWLF/D09ESzZs2wadMmvHv3Tt2xEVJldHR00LFjR3Ts2FEuSWfnF2LZ36FYefa5/HE8HlytDasqTELIR0alJD1+/Hhcu3YNr1+/xtChQ7Fp0yY4Ojpi2LBhOHfuHCQSibrjJEQjAl4no+/mm9gX+BYA0K6OBfj/vzAVjXUmhFQ2pWccK41UKsVvv/2GhQsXQiwWw8HBAbNnz8bcuXM/mkUoaMax6k0qlSI+Ph5A0e2cnAIp1p5/joN3i2arczQ3wLphzdC5ng3iRbk01pkQUiFqn3FMkaysLBw7dgy7d+9GYGAgunfvDm9vbyQlJWHz5s24c+cOTp06VZFTEFIlCgsLsXPnTgBAl5FfYOnfzxGbngsAGNveGT4DGsFYWPTfhcY6E0KqikpJ2t/fH7t378axY8dgZmaGSZMmYf/+/ahTpw5XZtSoUUovQ0mINhAzHdwrqI09fz0EANS2MMBPw5rjE3fFU9ASQkhlUylJd+vWDQMGDMDhw4cxYMAAhb1ha9WqhcmTJ1c4QEKqws1XyTiV3wQ5rGj6zgkdXPBtv4YwEtJUAoQQzVHpL1B0dDTs7e3LLFdZ60oTUlHFs4ZZGwux4+YbHH8QA0APJrw8/D6pIzo3sCuzDkIIqWwqJWllEjQh2krRrGE8HtCIn4hWgli0rzNYc8ERQsh7VErS76/D/CGhUAg3NzdMnjwZEyZMUDkwQipDSbOG/TKqOR7+vVszQRFCSAlUStLTp0/HTz/9hAkTJqBVq1bg8Xi4f/8+9u/fj9mzZ0MsFuPLL79EYWEhpkyZou6YCVFZYHiKwlnDzA1oxjBCiPZRKUnfunULBw8ehJeXF7ftiy++QJ8+fbBt2zZcuXIFnp6eWLVqFSVpojUuPk3A936hctt1eDzUsTFB165di56XMC0oIYRUNZUmMzEzM0NMTAxMTGTXys3IyICzszPS09ORkZEBBwcHZGVlqS1YbUaTmWivQokUP196iW3XwwEALpaGiE7LgZT9N2vYqLY0XJAQUnUqdTITY2NjnD9/HqNGjZLZfu7cORgbGwMA4uLi4Obmpkr1hKhNSlY+Zh8Jxu3XKQCAyR1dsWRAIyRn5dOsYYQQradSkl68eDEmTJiA8+fPo02bNmCM4cGDBzhy5Ag2bNgAAFi/fj0WLFig1mAJKY+Q6HR8eeAB4kV5MBDoYN2wZhjSwhGA/KxhjDFukRgbGxvweDyNxEwIIe9Tee7uf//9F1u2bMHz58/B4/HQsGFDzJkzB/369VN3jNUCXe7WHowxHLwbhRVnnqJAwuBmbYRt41qjgZ1JiceIxWKsXbsWAODj4wM9Pb2qCpcQ8hGq1MvdZ8+exaBBgz7ahEy0V65Ygu/8nsD3YSwAoG+TWtgwwgMm+tR7mxBS/aiUpIcOHYqCggK6JEi0ytuUbMw48BDP4zPA5wGL+jXE9C5u9DklhFRbKiXpunXr4tmzZ2jSpIm64yFEJVeeJ2Lu0RBk5hXCykgPW8e0pIUxCCHVnkpJetGiRRg/fjx++uknNG7cWO7+nbV15f9xLCgogECg/CVMqVQKPp+v9nqJ5sSLchGelI0rzxOxJyASANDS2Ry/j21FPbYJITVC2VlLAW9vbwQHB6N3795wdHSEjY2NzKMyrVy5ElZWVtDX10ejRo1w+fLlUsufO3cO3bt3h6mpKYyNjdG3b1+EhspPaFHeeolmHb0XhY7rrmLcrrtcgp7QwQVHv+hACZoQUmOolKSDg4NLfVSW3377DevXr8fJkyeRlZWFzz//HIMHD0Z4eLjC8hKJBNu2bcMPP/yAd+/eITo6GpaWlujTpw9EIpHK9RLNUjT/No8HfNmtLvR0VfpIE0KIVlJ5CJYm1KtXD4MGDcKmTZu4bS4uLhg5ciTWr1+vVB0xMTFwcnLChQsX0KdPH7XVS0Owqs6hu2+x5JT81ZDD0zzRoa6VSnVKJBJcuXIFANCzZ0+aGpQQUqkqdQgWAOTl5eHGjRt48+YNvvzySwBAVFQUnJycKqU3bUpKCl6/fo0uXbrIbO/atSvu3LmjdD2xsUVDc6ysrNRaL6ka9yJTsfrcc7ntOjweXK0NVa5XR0eH+9JGCCHaQqVrg+Hh4WjWrBlGjx6Nr776itu+aNEinDx5Uul68vPzkZWVVeqjuKGfmJgIAHL3vG1sbJCUlKTU+cRiMebMmQNPT0+0atWqQvXm5+cjIyND5kEq17WwJIzfdRfZYglcrQzB///vgsXzb9O9aEJITaNSS3revHkYOHAgNm7cKHNZcP78+fj6668xfPhwpeqZM2cODhw4UGqZx48fy8wBLpVKZfZLpVKlWu4SiQTjxo1DXFwc/P395Y4pb71r167FihUryjwvUY/Tj+Iw/2gICqUM3RrYYNvY1kjPFatt/m3GGNdPwczMjMZWE0K0gkpJ2t/fH3v37pUb0tS4cWM8evRI6Xr++OMP/PHHH0qVdXBwAAC51m1SUhLs7e1LPVYikWDChAkICAjAjRs34Oz834pHqtbr4+OD+fPnc88zMjLg5OSk1Gsh5XPgzlt8/3coGAMGezjg5xEe0NPlw0DPQG2t54KCAmzZsgUATQtKCNEeKl3uLiws5Fqe77c4oqOj5ZavVBdzc3M0adKE69wDFLV2r169ik6dOnHb8vPzkZubK1Nm4sSJuHHjBq5fv466deuqVO+HhEIhTE1NZR5EvRhj+O3aayz1K0rQ4zydsXlUC+rBTQj5aKj0165Hjx7YvHkzgP+SdHp6OubOnYu+ffuqLbgPLV68GHv27MHRo0cRFRWF2bNnIzc3l+u4BgAzZ85E27ZtART9kZ8yZQouXbqEs2fPws7OjrvXXVhYWK56SdVijGHdP2FYf+EFAODr7u5YOaQpdPh0GZoQ8vFQ6XL3zz//jK5du+Ls2bNgjKF///64c+cOTExM4O/vr+4YOePGjUNOTg5WrFiBxMRENGvWDJcvX+YuWQOAvr4+DA2LevmmpqbixIkTACDXKv7tt98wceJEpeslVUciZVji+wRH70cDAL4b0AjTutDa5ISQj4/K46TT0tKwd+9e3L9/H1KpFK1atYK3tzcsLCzUHWO1QOOk1SO/UIJ5R0Nw/kkC+DxgnVdzjGxb+ff6aalKQkhVqvRx0hYWFpg3b56qhxMiJ0dciOn7H+DWq2To6fCxZXQL9G9WeqdAQgipyVRO0vn5+QgPD0dqaqrcvtI6XBGiiCinAJP3BuFhVDoM9XTw5/jW6FyvcueBJ4QQbadSkr527RrGjBnDTQTyoWo00yjRsHhRLoLfpuHnSy8R/i4bZgYC7JncFq2cq/a2CZ/PR5s2bbifCSFEG6h0T7pZs2bo1asXFi1apPAetL6+vlqCq07onnT5Hb0XJbNQholQFye+/AQN7CpnGB8hhGiLSr0n/ebNG6xcuRLGxsYqB0g+bopWssoWF8LUQOU7MIQQUuOodF2vfv36iImJUXcs5CPyJEYkk6ABQMqAyOQcjcTDGEN2djays7Ppdg0hRGuolKRnz56NSZMmITAwEO/evUNycrLMg5DS5BVIsPXqa7ntFV3JqiIKCgqwYcMGbNiwAQUFBRqJgRBCPqTStcUpU6YAAD755BOF+6klQkoikTLMPRKCJ7Ei6OvyIZZIIWW0khUhhCiiUpIODg5WdxzkI8AYw8qzz/Dv0wTo6fCxd0o7uFgZqm0lK0IIqWlUStItWrRQcxjkY7Dj1hvsDYgEAPw80gOeblYAQMmZEEJKoPKA0Ly8PFy4cAHbtm3jtkVFRdGlbqLQ6UdxWHM+DEDRXNyDPWhedEIIKYtKSTo8PBzNmjXD6NGj8dVXX3HbFy1ahJMnT6otOFIzBIanYOGxonXGJ3d0hXfnOhqOiBBCqgeVkvS8efMwcOBApKSkyGyfP38+fvrpJ7UERmqGFwmZ+GL/fYglUvRvaoelAxvLrEFOCCGkZCrdk/b398fevXvlpk9s3LgxHj16pJbASPWXIMrDpD1ByMwrRFtXC2wa1UJr14Pm8/nw8PDgfiaEEG2gUpIuLCyEVCoFAJlWUXR0NExMaEpHAmTkFWDSniDEi/JQ18YIOya0gb5AR9NhlUhXVxdDhw7VdBiEECJDpSZDjx49sHnzZgD/Jen09HTMnTsXffv2VVtwpHoSF0oxY/8DhCVkwsZEiL2T28HckNZnJoSQ8lKpJf3zzz+ja9euOHv2LBhj6N+/P+7cuQMTExP4+/urO0ZSjUilDItOPEJAeAqM9HSwZ1JbOFlqZhax8mCMcTONCQQCum9OCNEKKiXpunXr4smTJ9i7dy/u378PqVSKJUuWwNvbW+GqWOTjsf7iC/iFxEGXz8Pv41qjqaOZpkNSSkFBAdauXQsA8PHxgZ4etfwJIZqn8pJDFhYWmDdvnjpjIdXc/sBIbLseDgBY69UMXevbaDgiQgip3mhdQFJh8aJcHL8fg02XXgIA5veujxFtnDQcFSGEVH+UpEmFHL0XhcW+T1A80Vy7OhaY1cNds0ERQkgNQQNCicriRbnweS9BA8D9yDQkZORpLihCCKlBKEkTlYUnZUP6wVTtUgZEJudoJiBCCKlhVErSWVlZCufoPnnyJLKysiocFKke/F+/k9umw+PB1Vr7h1wRQkh1oFKSXrRoEZKSkuS2JyUlYfHixRUOimi/kOh07LwVAQAoHlKsw+NhjVfTarn0JJ/PR+PGjdG4cWOaFpQQojV4TIW1Ja2trfHixQtYWVnJbE9JSUGTJk2QkJCgtgAVef78ORITE9GoUSPUqlWrzPK5ubkIDQ2Frq4uGjZsCAMD2STy7NkzREVFyWwzMTFBx44dlY4pIyMDZmZmEIlEMDU1Vfq46igrvxADf7mFtyk5GNTcHksGNMTblFy4WhtWywRNCCFVTdmcofLc3SKRSC5Jp6enIzc3V5UqlZKdnQ0vLy8EBQXB3d0doaGhWL58eYmtd8YYlixZgt27d8PFxQU5OTlITEzE77//jhEjRnDlfvnlF/j5+aFFixbcNldX13Il6Y/JD6ef4m1KDhzNDbD6s2YwMxDAwZwucRNCiLqplKR79OiBxYsX46+//uJapbm5uVi0aBF69Oih1gDft2zZMrx8+RKvXr2CtbU1Lly4gH79+qFTp07o1KmTXHnGGMzMzBAeHg5jY2MAwE8//YTx48fjk08+gaOjI1e2U6dOOHHiRKXFXlOcfRyHEw9iwOcBm0a1gJmBQNMhEUJIjaXSzbf169fD398fderUwaBBgzBw4EC4uroiMDAQGzZsUHeMnH379mHq1KmwtrYGAPTt2xctW7bEX3/9pbA8n8/H4sWLuQQNABMnTkR+fr7ckprZ2dm4efMmHj16hLw8GkKkSGx60ZArAPiqmzva1bHUcETqIxaLsWLFCqxYsQJisVjT4RBCCIAKzt29e/duPHz4EDweDwsXLsTUqVNhaVk5f7hjYmKQnJyMli1bymxv2bJludawDggIAADUr19fZvvNmzchEokQFxeHnJwcbNu2DcOGDSuxnvz8fOTn53PPMzIylI6hOpJIGeYdCUFmXiFaOJljTq96mg6JEEJqPJVnHLOyssI333xToZM/efIEsbGxpZbp0qULDA0NkZ6eDgByXwKsrKyQlpam1Pni4+Mxa9YsTJw4Ee7u/82K9emnn+J///sfzMzMwBjDjz/+iLFjx6Jp06Zo0KCBwrrWrl2LFStWKHXemuCPG+EIikyFkZ4OtoxuAYEO9YAmhJDKpnSSjoyMBFDUoar455K4uroqVefFixdx6dKlUst4eHjA0NCQW5Xow45pOTk5Sq1YlJycjD59+qBevXrYtm2bzL4BAwZwP/N4PHz//ffYsmULzp49W2KS9vHxwfz587nnGRkZcHKqmfNVh0Snc/NyrxjSFC5WRhqOiBBCPg5KJ+k6deoAKOqMVfxzSZQd1bVgwQIsWLBAqbJOTk7Q0dFBdHS0zPaYmJgyvxSkpKSgV69esLS0xNmzZ+WGYH2Iz+fD3Nwc8fHxJZYRCoUQCoVKxV6dZeUXYs6RYBRKGQY1t8ewVo5lH0QIIUQtlE7Sr169UvhzVTEwMECXLl3g5+eHiRMnAgBEIhEuX77MrQMMAE+fPkV6ejo3fCo1NRW9evWCmZkZzp8/DyMj2VagVCpFVlaWzDi1p0+f4u3btzJDsj5WHw634hXPXEIIIaTSKZ2k37+HO3fuXJw9e1ZhuUGDBpW4r6LWrFmDbt26YdasWejQoQO2bdsGJycnTJ06lSuzadMm3LlzB6GhoRCLxejduzfi4uKwdOlS3Lp1iyvXtGlT1K5dGxKJBO3atcPnn3+OJk2aICoqChs2bECnTp0watSoSnkd1QUNtyKEEM1SqePYuXPnFG6XSqX4559/KhRQaTw9PREYGIht27bh6NGj6Nq1K+bPnw9Dw/8m0mjatCl0dYteVl5eHmxsbGBjY4MdO3bI1DVv3jzUrl0bAoEAAQEB+OOPP3DkyBFYWFhg3bp1GDt27Ec9PWRNHm6lCJ/PR7169bifCSFEG5RrWtCYmBgARfeHP7w3LJVKcfv2bSxatEhu38egJk0LKpEyjNl+B0GRqWjhZI7jMzpQb25CCFGjSpkW9P3ey4p6Muvp6WHTpk3lqZJooW3XX9NwK0II0QLlStLPnz8HADRq1Ij7uZhAIICDg0OZPaeJdguJTsemy0UdA2m4FSGEaFa5knTDhg0BAO/eveOm5iQ1R/FwK8lHONxKLBZzU9ouXLhQqbH3hBBS2VTqOEYJuuaJF+Xi25OPP+rhVgUFBZoOgRBCZKiUpKVSKXbt2oXjx48jKioKhYWFMvtfv36tluBI1Th6LwqLTz5BcQ/CwR72NNyKEEK0gEo9gtatW4cff/wRPXr0wIsXL/D111+jXbt2iIiIwKBBg9QdI6lE8aKioVbvd/HfcTMC8aLKWxecEEKIclRqSRe3oj09PeHj44O5c+cCKFoM49SpU+qMj1SyiORsSD8YhCdhDJHJObA3o06AhBCiSSq1pN++fYvWrVsDKJquMzMzEwDw+eef4/bt2+qLjlS6/AKp3DYdHg+u1oYKShNCCKlKKiVpiUQCgaDonqWLiwsCAwMBAG/evPkoFp2oKRhj2HY9HABQ3EVMh8fDGq+m1IomhBAtoPJ60sWmTZuG0aNHo0uXLggICMDIkSPVERepAheeJiAoMhX6Aj6OTPNEboEUrtaGH2WC5vF4cHFx4X4mhBBtoFKSLr68DQDz58+Ho6MjAgMD0a9fP3h7e6stOFJ58gslWHM+DADwRWc3tHC20HBEmiUQCDBp0iRNh0EIITJUutw9evRomeejRo3C5s2bMWPGDAwdOlQdcZFK9ldAJKJSc2BrIsT0rnU1HQ4hhBAFVErSmloFi6hHSlY+tl4pGsu+sG8DGAkrfNeDEEJIJSjXX+fiVbA+/Bn4bxUsBwcH9URGKs2myy+RmV+IJg6mGN6qtqbD0QpisRhbtmwBAMyZM4emBSWEaAVaBesj8zIxE4fuRgEAvh/UGHw+dZIqlpOTo+kQCCFEBq2C9ZFZfe45pAzo26QWPN2sNB0OIYSQUqhtFayEhATo6OioLzKidtdfJOHGy3cQ6PDg07+RpsMhhBBSBpU6jsXGxmLevHncc29vb9jb28POzg5BQUFqC46oT6FEitXniq5+TOzgCldrWieaEEK0nUpJesGCBdxQq2fPnuHYsWO4fv06vvrqKyxevFid8RE1OXwvGq+SsmBhKMCsnvU0HQ4hhBAlqDT2JigoCG3btgUAXLx4EUOHDkXXrl3RqlUr1K5NvYW1jSi3AJsuvQQAzOtdn5ahJISQakKllrShoSGio6MBFI2Z7tGjBwAgKysLRkZ0GVXb/HbtNVKzxXC3Ncbn7Zw1HY5W4vF4cHBwgIODA00LSgjRGiq1pD/77DN8+umn8PDwwL1793DkyBEARa3qvn37qjVAUjFvU7Kx53YEAOC7AY2gq6PS97IaTyAQYNq0aZoOgxBCZKj0F3vz5s2YMmUKrK2tceHCBVhZFQ3lCQkJwfLly9UaIKmYtefDUCBh6FzPGt0a2Gg6HEIIIeXAY4wxTQdRE2RkZMDMzAwikQimpqaaDgcAcPdNCkZtvwM+D/hnThc0sDPRdEiEEEKgfM6o0KTN0dHRCAsLA2MMjRo1UjgLGdEMqZRh1f8PuRrTzpkSdBkKCgrw22+/AQBmzpzJrZdOCCGapPJSlTNmzMDhw4dR3BDn8XgYM2YM/vzzTxgbG6s1yPelpKTg2LFjSExMRLNmzfDZZ5+Bzy/9qn1YWBguXbqErKwsNGnSBAMHDpSbeEWVerWZb3AsnsSKYCLUxbze9TUdjtZjjEEkEnE/E0KINlApC82fPx8hISG4ePEisrKykJWVhYsXLyIkJATz589Xd4yciIgINGvWDEePHkV2djYWLFiAwYMHQyqVlnjMkiVL8PnnnyMyMhLp6emYO3cuPD09kZWVVaF6tVmOuBDrLxStFT2zhzusjYUajogQQohKmAosLS3ZkydP5LY/efKEWVlZqVKlUkaMGME6dOjACgsLGWOMhYeHM11dXXbo0KESjwkJCZF5npCQwPh8Pjt8+HCF6v2QSCRiAJhIJCrPS6oUGy++YC7fnmWd/neF5RUUajqcaiE/P5/98MMP7IcffmD5+fmaDocQUsMpmzNUaklnZ2fD3t5ebrudnZ1MC1WdCgsLcebMGYwbN467VO3m5oauXbvC19e3xOM8PDxknkskEgCAiYlJherVVvGiXPx5MxwA4NO/EYS6NJ86IYRUVyrdk27Xrh1WrFiBjRs3Qle3qIrCwkL88MMPaN++vVoDLBYVFYW8vDy4u7vLbHd3d0dgYGCpx758+RK7d+9GRkYGbty4gaVLl2LgwIEVqjc/Px/5+fnc84yMjPK+pEqx/t8XyCuQoq2rBfo3tdN0OIQQQipApSS9adMm9OnTB76+vmjVqhUA4OHDh8jLy8OFCxeUruf48eN48OBBqWW++eYbWFlZITs7GwDkuqqbmZlx+0oiEAhgbm6OwsJC5OTk4OnTp8jJyYGhoaHK9a5duxYrVqwo9bxV7crzRPgGxwIAlg5sTDNnEUJINafS5e7WrVvj9evXmDt3LqysrGBtbY25c+fi9evXaN26tdL1GBoawtzcvNRH8SXo4h7jxT1wi6Wnp5fZm7xOnTpYvHgxNmzYgAcPHuDGjRvYsmVLher18fGBSCTiHsXTpGrKkaAoTP3rPvc8LEE7WvbVBY/Hg42NDWxsbOjLDSFEa6jUkvb29sbOnTuxcOHCCp184MCB3GXnsjg7O8PIyAgvXryQmXr0xYsXaNRI+bWRLS0t0bhxYzx9+rRC9QqFQgiF2tFrOl6UC59TT2S2LfENRZf6NrA3M9BQVNWLQCDAV199pekwCCFEhkot6UOHDiEvL0/dsZRKR0cHn332Gfbt2wexWAygaJlMf39/jBgxgivn6+uLzZs3AwDEYjHu3bsnU090dDRCQkK4DmXK1qvNIpKz8eHQXgljiEzO0UxAhBBC1EKlaUF79+6Nr7/+GkOGDKmMmEoUGxuLzp07w8LCAq1bt8bp06fRrVs3HD58mLtE6e3tjTt37iA0NBQFBQXo1asXBAIBGjdujPT0dJw+fRpdunTB0aNHYWBgoHS9ZdHktKBXwxIxZe99mW06PB78F3enljQhhGihSp0WtG3bthg7diwmTZqExo0bQ09PT2a/t7e3KtWWydHREY8fP8bp06eRmJiIkSNHolevXjJlhg0bhg4dOgAouoR548YNBAQEIDg4GEZGRliwYIHcsCxl6tVmfsFxAAAeAIaiBL3Gqykl6HIoKCjAjh07AADTpk2jaUEJIVpBpZa0q6trqfsjIyNVDKf60lRLOjY9F11+ugaJlGHv5LYQ6urA1dqQEnQ5icVirF27FkBRp8APv3gSQog6VWpL+vXr19z46A8VFhaqUiVR0d7bEZBIGTq4WaFbA1tNh0MIIUSNVOo4VtqlQLpMWHUy8wpwJKho6Ne0LnU0HA0hhBB1U+syT/n5+VozLOljcPReNDLzC+FmY4Ru9akVTQghNU25Lnf/+uuvCn8GAKlUiqCgIDRu3Fg9kZFSFUqk2HM7EgDg3ckNfD5NwEEIITWN2pK0QCCAq6srtm/frp7ISKn+fZqA2PRcWBrpwauVo6bDIYQQUgnKlaTDworWKO7UqRP8/f0rJSBSNsYYdtyKAACM83SBvoBWuqooHo8HMzMz7mdCCNEGKvXupgStWfffpuFRdDr0dPmY0MFF0+HUCAKBAHPnztV0GOQjIpFIUFBQoOkwSCURCATc2hMVoVKSfvbsGQ4ePIjVq1fLbF+yZAnGjx9frrm0SfntvPUGAODV0hHWxtRRj5DqhDGGhIQEpKenazoUUsnMzc1hZ2dXoatzKiXpWbNm4bvvvpPb3qtXL8yZMwcXL15UOSBSusjkbFx8lggAmNqJhl0RUt0UJ2hbW1sYGhrS7ZUaiDGGnJwcJCUlAQDs7e1VrkulJB0YGIh27drJbW/bti0CAgJUDoaUbc/tCDAGdGtgg3q1TDQdTo1RUFCAvXv3AgAmTZpE4/1JpZBIJFyCtrKy0nQ4pBIVrw2RlJQEW1tblS99qzRO2srKCiEhIXLbg4ODYW5urlIgpGzpOWIcux8DAJjW2U3D0dQsjDHExcUhLi4OKsyUS4hSiu9BGxoaajgSUhWKf88V6XugUpIeO3YspkyZgqtXr0IsFiM/Px9XrlzBlClT8Pnnn6scDCndwbtRyC2QoKGdCT6pS9/CCamu6BL3x0Edv2eVLnevWLECERER6NmzJ/j8ojwvlUoxcuRIrFy5ssJBEXniQin+CogEUNSKpv/khBBStUJDQ/HixQvY2tqic+fOKCwsREBAAOLj49GmTRvUrVtX7edUKUkLhUIcPXoUK1euRHBwMHg8Hlq2bIl69eqpOz7y/848ikNSZj5sTYQY7OGg6XAIIR+p5ORkBAcHIycnB87OzmjYsCF3/xUAsrOzcebMGYXHdurUCcnJydycGyUZOnQo9PX11Ro3AKSlpeHChQsYNmxYufudrFq1Cr/88gu6dOmCdu3aoWPHjujatStEIhGaNm0KKysr7UnSxerXr4/69eurKxZSgqLJS4qGXU38xBV6umqdcp0QQsokEokwa9YsHD9+HK1bt4aNjQ1iYmIQHx+PSZMmYdWqVQCAd+/eYcyYMejevTtsbWXXFKhTpw5ev34tk8RPnTqFhg0bygzd7d+/f6Uk6fDwcIwZMwZpaWnl7j+1fft2/Pzzzxg/fjwA4PHjxwgICEBqaiosLCzUHmsxlZP03bt3sXfvXrx58wYXLlwAAOzevRsjR46EsbGx2gIkQEB4CsISMmEg0MHY9s6aDocQ8pGRSCQYMGAAMjIyEBYWBheX/yZRysjIwMGDB+WOWbZsGbp16ya3vX379hg7diz33NraGmPGjMHixYuViiU1NRVBQUEQCoXw9PTkWvGxsbG4c+cOhg0bxpXNzMzEuXPnMGTIEEilUm54sK+vLwwNDeHi4oIOHToAACIjIxESEgJTU1N88skn3JcEkUiEf/75BwkJCQgNDcWRI0cAFF36BsDlPy8vr0pZh16lJtmZM2fQvXt3ZGdny4yJjouLw8aNG9UWHClS3Ioe2aY2zA3V/yEgRQwNDanXLSEK+Pr6IiAgANu3b5dJ0ABgamqKL7/8skri+P333+Hk5IQffvgBq1atQtu2bfHy5UsAwL179zBx4kSZ8rGxsRgzZgxSUlKQk5ODGzduAADOnTsHPz8/3Lt3DwDw3XffoUmTJvjtt98wc+ZM1KtXD48fPwZQlOj9/PwgkUgQGBgIPz8/+Pn54fbt2wDAPc/Pz6+cF81U0LJlS3by5EnGisaqcNtfvHjBXFxcVKmy2hOJRAwAE4lEaq33ZUIGc/n2LHNdfJZFvMtSa92EkKqVm5vLnj17xnJzc+X25efnl/goKChQuqxYLC6zbHl5e3szOzs7pcpGREQwAOz7779nhw8flnnk5OTIlbeysmJr164ts15/f3/G4/HYiRMnuG1v3rxhjx8/ZowxdurUKWZkZCRzzPPnzxkAFh0dzRhj7N69ewwAS0tLk6v31q1bjDHGJBIJGzZsGGvXrh2TSqVcOaFQyM6cOcM9v3TpEisrhZb2+1Y2Z6h0uTssLAz9+vUDINvF3MHBAXFxcRX82kDet8u/aCGNPo1rwdXaSMPREEIqy9q1a0vcV69ePZnhrRs2bChx7K2LiwsmTZrEPd+yZQtycnJkyixfvrxcsSUkJMDZWfZW26tXr/DgwQPued++fWXuzfr7+3Ot3GJ9+vSR6WRWHvv27UP79u1lLmfXqVPxWRePHDmCzp07o1OnTgAAPp8PHx8ftGnTBhEREXBz0+ycFColaQsLC0RFRaFhw4YySdrf3x9OTk5qC+5j9y4zH77BsQBo8hJCiOYYGxvj+fPnMtsiIiLg5+cHkUiEf//9F/fu3UObNm24/SXdk1ZVVFRUpXRUfvv2rVwiLu6lrWhfVVMpSY8dOxZfffUVdu/eDQDIzc3FxYsXMXPmTHh7e6s1wI/Z/jtvIS6UwsPJHK1dKq/3ICmaEai488vYsWNpWlBS5Xx8fErcVzwfRbGFCxeWWPbDORTmzJlTscAAtGvXDidOnEBiYiJq1aoFoKhV3KdPH4SGhuLff/+t8DnKYm5ujpSUlBL38/l8SKVSmW15eXll1mttbS1Xb2pqKrdP01TqOLZy5UrY29ujTp06kEqlMDY2xtChQ9GtWzeFC2+Q8ssrkODAnbcAgGmd69DkJZWMMYa3b9/i7du3NC0o0Qg9Pb0SH7q6ukqX/fALpqIy5TVx4kSYmppiwYIFKCwsrNDrVFWfPn1w9epVxMTEcNukUimXYB0dHZGbmytzy/XatWsydRSPPHo/eXfq1AnXrl1DWloat+348eOwtbXViiHGKk9mcvDgQaxcuRIPHjyAVCpFy5YtteIF1RS+D2ORmi2Go7kB+jWx03Q4hJCPmKWlJU6fPg0vLy+0bdsWn332GVxcXCASiXDmzBnY2NjIjRW+evUqEhISZLa1bNkSDRo0UCmG8ePH48iRI+jQoQNmzZoFQ0NDblKtbt26oUWLFmjWrBlGjhyJyZMn4+XLlzh8+LBMHa6urrCyssKyZcvQvXt3uLq6YuLEifjjjz/QrVs3TJ8+HVFRUdi0aRN27NgBoVDzSwFXaDITNzc3jV+vr4mkUoad/kXDriZ3dIWuDk1eQgjRrI4dOyI8PBzHjh3D/fv38ebNGzg6OmLKlCnw8vLiEpqRkRFGjRqFly9fynUcMzY2lkvSXl5eaNy4cZnn19XVxblz53D48GHcuHEDJiYmWL9+PTw9PQEAOjo6uH79On755RcEBgaiSZMmuHLlCr7//ntuaKW+vj4uX76Mffv24cyZM/D09ESHDh1w/fp17Ny5E0FBQTAxMcHly5fRuXNnmfOPHDkSjo6O3HM7OzuMGjWq/G9kOfGYktf2Nm/erHSlc+fOVTGc6isjIwNmZmYQiUQwNTWtUF1XwxIxZe99mAh1EeDTAyb6dH+0sonFYq53rY+PT6VMSkBIXl4eIiIiUKdOnUqZUYtol9J+38rmDKVb0jt37lQ6sMpM0gEBAfjtt9+QmJiIZs2aYfHixVxHBkUKCwuxb98+XLhwAVlZWWjSpAlmzpwpMyB/w4YN8PPzkznOzc0N+/btq6yXUaodN4uGXY1u50QJmhBCPmJKJ+niKdA06ebNm+jVqxfmz5+PkSNH4tdff0XHjh0REhJS4lSkEydOhIWFBUaPHg1dXV3s2LEDbdq0QXBwMGrXrg0AeP36NQQCgcwKXkZGmhmTfP1FEgLfpIDPAyZ1rPgYQEIIIdVXhe5JV7XvvvsOQ4cOxbp16wAAvXr1gr29PbZv34758+crPGb79u0yCbdv374wMjLC5cuXZQb8W1lZcYPZNeXovSh8e/IJAEDKAP9X7zCqLc3VXVVo2BUhRNuo3CPp7t27+PLLL9G3b19u2+7du5GVlaWWwD6Uk5ODgIAAfPrpp9w2IyMj9OrVC5cuXSrxuA9bxEFBQZBIJGjatKnM9vv376NPnz4YMWIEfvnllxJn86ks8aJc+Pg+kdm2xDcU8aLcKo3jY6Wnp4clS5ZgyZIldD+aEKI1qs0CG9HR0ZBKpTK964CisXFRUVGlHnv37l106tQJzZs3x5AhQ3Dq1CmZmXEMDAwwatQozJ07F3369MGmTZvQs2dPSCSSEuvMz89HRkaGzKMiIpKzIf2gC5+EMUQm5yg+gBBCSI2n0uXu5cuX48CBA/Dy8sL+/fu57SNHjkSfPn2wbNkypepZs2YNzp8/X2qZo0ePwtHRkWvZfthDzsDAAGKxuNQ6GjRogHXr1iE5ORm7du3C3Llz0bp1a+6e9Nq1a2Xq7dq1Kxo3bowTJ06U2MV+7dq1WLFiRZmvUVl1rI3A50EmUevweHC1plWZCKlpaMKcj4M6fs8aXWBjxIgR6NKlS6llrKysABQNpgcgN31bSkoKt68k5ubm3P3mwYMHo379+ti8eTM2bNgAQD7x169fHy4uLggJCSkxSfv4+MjcB8/IyKjQvOX2ZgZY69UMS3xDIWEMOjwe1ng1hb2ZapPRk/IpLCzEsWPHABR92fxwhidC1KG430NOTo7KC02Q6qN4YZOK9HfR6AIb9erVQ7169ZQq6+DggFq1auHBgwcYNGgQtz0oKAgdO3ZU+pw6Ojqws7NDUlJSiWUkEglSUlJKXVtYKBSqfTaaUW2d0aW+DSKTc+BqbUgJugpJpVK8evWK+5mQyqCjowNzc3Pu74+hoSFN+VsDMcaQk5ODpKQkmJubQ0dHR+W6qtUCG5MnT8bOnTvxxRdfwN7eHn5+fggNDZUZw/2///0PT58+xb59+5Cbm4vt27dj5syZXMvo/PnzCAoKwsyZMwEUTWKxdetWzJ49GwKBABKJBEuWLEF2drbMkmhVxd7MgJIzITWYnV3RNL+lNRRIzWBubs79vlWlUpJeuXIlpkyZwq3laWxsDKlUirFjx1bqAhvLly/Hixcv4O7uDhcXF0RGRmLLli1o3749V+bVq1d4+PAhgKLWbnx8POzt7eHk5IT09HSIRCKsXr2aW5tVIBAgOTkZDg4OqF27NuLi4mBoaIhTp04pNVUdIYSUB4/Hg729PWxtbat8FAmpOgKBoEIt6GJKTwuqyJs3bzSywEZUVBQSExNRv359mJmZyex7/fo1MjMz0bJlS25bXl4eXrx4ASMjI7i4uCi8P1BcxsLCArVr15ZbGq4s6pwWlFQ9mhaUEFKV1D4tqCKaWmDD2dkZzs6KJ/lwd3eX26avrw8PD49S61SmDCGEEFKVaHklQgghREvROBM1Kb5rUNFJTYhmiMVibiH4jIwMutxNCKlUxbmirDvOFbonTf4TExNToXHShBBCPj7R0dHcxFqKUJJWE6lUiri4OJiYmNC4xzIUT/wSHR1NneyUQO9X+dD7VT70fpWPut4vxhgyMzPh4OBQakdlutytJnw+v9RvQ0Seqakp/VEoB3q/yofer/Kh96t81PF+fTg6SRHqOEYIIYRoKUrShBBCiJaiJE2qnFAoxPLly9U+93lNRe9X+dD7VT70fpVPVb9f1HGMEEII0VLUkiaEEEK0FCVpQgghREtRkiaEEEK0FCVpohF3797FvXv3NB1GtRATE4OQkBCIxWJNh6L10tLSEBwcjJcvX6KwsFDT4WiduLg4+Pv7lzp9cVRUFO7fv09THKNozW9/f3+kpaUp3C+VSvHixQu8ePGi8pYdZYRUsV27djE+n88cHR01HYpWi4mJYd27d2dmZmasTZs2zM3Njf3777+aDksrSSQSNmPGDGZgYMBatGjBnJycmJOTE7t8+bKmQ9MKd+/eZZ999hmzsbFhANitW7fkyuTk5LAhQ4YwQ0ND1rBhQ2ZgYMB+++03DUSreY8ePWKjR49mtWrVYgDYmTNn5MqsX7+eOTg4sAYNGjA3NzdWq1Ytdvz4cbXHQi1pUqXCwsKwbNkyeHt7azoUrZafn4++fftCX18fsbGxuHfvHoKCgpCenq7p0LTSyZMnsX37dgQGBiI4OBhv375F7969MWHCBE2HphVCQ0MxduxYBAQElFhm2bJlCA4Oxps3b/D8+XPs27cPX3/9NR48eFCFkWqHp0+f4tNPP8Xjx49LLCMSiRAcHIywsDCEh4dj4cKFGDt2LCIiItQbjNrTPiElyM3NZc2bN2dHjhxhK1eupJZ0Kf766y/G5/NZdHS0pkOpFn777TdmYmLCpFIpt23Hjh1MX1+fSSQSDUamXSIiIhS2pKVSKbOysmKrVq2S2d6oUSM2c+bMqgxRq6SlpZXYkv5QVlYWA8COHj2q1hioJU2qzIIFC+Dh4YFRo0ZpOhStd+XKFbRs2RKOjo4IDQ2t3HteNcCYMWPg5uYGb29vXLx4EYcOHcK6deuwdu3aUhcvIEWio6ORkpKC1q1by2xv27YtgoODNRRV9VLcx8bd3V2t9dICG6RKnDp1Cv/88w9CQkI0HUq1EBcXBwMDA3h6eiI7OxtZWVkoKCjArl270K9fP02Hp3UsLCwwa9YsfPvtt7h//z5SUlJQr149DB48WNOhVQupqakAACsrK5ntVlZWCAoK0kRI1YpIJMKMGTMwYMAAtGrVSq1101dMUumysrIwdepUzJgxA48fP4a/vz+ioqIgFovh7++P5ORkTYeodQQCAfz9/bFw4UKEhoYiIiICY8aMwZgxYyASiTQdntb566+/MHv2bFy5cgWPHj1CVFQUmjVrhu7duyM3N1fT4Wk9gUAAAMjLy5PZnpubCz09PU2EVG1kZ2dj0KBBMDAwwIEDB9RePyVpUunEYjEaN26M06dPY/HixVi8eDEuXLgAkUiExYsX49GjR5oOUeu4urrC0tISI0aMAADweDxMnz4d6enpePLkiYaj0z5nz55Fp06d4OHhAaBo6dgvv/wS0dHRdLlWCc7OzuDxeIiNjZXZHhsbC2dnZw1Fpf2ys7MxcOBAZGRk4PLly7CwsFD7OShJk0pnaWkJf39/mce0adNgY2MDf39/9OzZU9Mhap2+ffsiKytLZqxqTEwMAMDGxkZTYWktGxsbxMbGgr23FEF0dDS3j5TOxMQE7du3x+nTp7ltmZmZuHr1Knr37q3ByLRXTk4OBg4ciLS0NFy5ckXuVoG60D1pQrTQp59+ig4dOuCzzz7DvHnzkJmZiWXLlmHIkCFo0KCBpsPTOtOnT8fu3bsxadIkjBkzBklJSVi+fDn69eun9o481VFiYiJevXqFhIQEAOCuxri4uMDJyQkAsHr1avTt2xd16tRBu3bt8Msvv8DBwQFTp07VWNyakpKSgufPnyMrKwsA8OzZM5ibm6N27dpwdXWFVCrlhmjt2bMHYWFh3LFubm5wcHBQWyy0ChbRiH379sHPzw++vr6aDkVr5eTkYOPGjbh58yZMTEzQvXt3TJ8+nbt/SGQ9e/YMv//+O169egUTExN07twZ06dPh76+vqZD07izZ89i3bp1ctu/+OILmbHkt27dwu+//46kpCR4eHhg8eLFsLW1rcpQtcL169exdOlSue1jxozBzJkzIRaL0aNHD4XHzpkzh7tNpQ6UpAkhhBAtRfekCSGEEC1FSZoQQgjRUpSkCSGEEC1FSZoQQgjRUpSkCSGEEC1FSZoQQgjRUpSkCSGEEC1FSZoQ8tF6+PAhrl27BpougmgrmhaUEPJR2r59OzZs2ABzc3O0adMGv//+u6ZDIkQOJWlCSJWLiIjAkydP8Omnn1bqeZ49e4bHjx8DAFq3bo169epx+44fP45///0XLi4uaNSokdyxMTEx8Pf3BwDUr19f7esEE6IMutxNCKlUb968wZkzZ2S23bhxA1999VWln9vX1xdff/01/Pz8EBERIbOve/fuWLduHdauXaswAcfGxsLPzw8LFy7E7t27Kz1WQhShljQhpFJdvXoVq1atwuDBg7ltderUwZAhQ6rk/O7u7jhy5Ijc9qlTp8Le3h48Hg9nz56V29++fXscOXIEQ4cOrYIoCVGMkjQhWi4lJQUPHjyArq4u2rRpA1NTU5n9b9++RUhICKysrNCqVSsYGhpy+169eoVXr16hd+/eePToEWJjY9G+fXvY2dlxZRISEnD9+nUAgIGBARo0aICGDRvKnOP9eu7cuYP4+Hh89tlnSElJKfXY2NhY3Lt3D9nZ2Vyi9PDwgLOzM/r27Sv3Wh8/fozXr1/D3t4e7du3B5//38W+4hj69OmDx48fIz4+Hi1atICjo6NK7+v+/fvRqFEjtGjRAvv27UP//v1VqoeQykRJmhAt5ufnhwkTJqB169bQ09NDeHg4/vzzT/Ts2ROMMcydOxcHDhyAp6cnkpOTER8fj5MnT6Jt27YAgAsXLmDVqlVwdHSEgYEBpFIpQkJCcOTIEe5+cGJiIvz8/AAA2dnZuH37Nvr3748DBw6Ax+PJ1OPs7Ax9fX04ODhg0KBBZR4bHx+P4OBg5OTkcOUMDQ2RmpqKpUuXcjEUFhZi5MiRuH79Ojw9PfH48WPY2dnhn3/+gY2NDRfD6tWr4erqCqFQCAC4e/cujh49qtK97T179mDq1Klo0aIFBg4ciLS0NFhYWKj0eyKk0jBCiNZq06YNW7NmDfc8JSWFXb9+nTHG2Pbt21ndunVZcnIyt3/VqlWsYcOG3POtW7cyAGz9+vXctuXLlzM7OzuWk5Oj8Jzv3r1j9vb27MSJE3L17Nmzp9R4FR27Y8cO5uLiIlNuz549zNHRkXv+66+/MktLS/b27VvGGGMZGRmsefPmbNq0aXIxvF/3woULWfPmzUuMZ+XKlax9+/Zy2wMDA5lAIGBJSUlMKpUyNzc39uuvvyqsY8iQIWzmzJmlvm5CKgt1HCNEixkYGODly5fIy8sDAFhaWqJr164AilqCzZs3x7Vr13D8+HEcO3YMxsbGCAsLQ3x8PFeHUCjErFmzuOcLFy5EUlISbt68yW0rKCjAnTt34Ovri8uXL8PJyQlBQUEysZiammLixIlyMSpzbFmOHDmCcePGwdnZGQBgYmKC2bNn4+jRozLlbGxsMGzYMO55t27d8OLFi3KdCwB27dqFwYMHw8bGBjweD1OmTMGuXbvKXQ8hlY0udxOixbZu3Ypp06bBxsYGn3zyCT799FN4e3tDKBQiMjISYrEYJ06ckDlm1KhREIvF3HN7e3vu8jAAGBsbw8bGBm/fvgUAPHnyBAMGDODuKRsZGSExMRFJSUky9drZ2XGXv4spe2xZ3r59i+HDh8tsq1u3LjIyMmQuQ1taWsqUEQqFyM/PL9e5cnJycOzYMYwfP567T25gYIDg4GA8evQIHh4e5aqPkMpESZoQLebh4YGgoCDEx8fjypUrWLFiBW7fvo1Dhw7B1NQUPXr0wE8//VRqHenp6TLPGWNIT0+HtbU1AOC7775Djx498Ndff3Fl+vXrJzcL14cJujzHlsXa2hqpqaky21JTUyEQCOQ6ylXU8ePHIRQKkZyczN0nB4CGDRti165d+OWXX9R6PkIqgi53E6LFYmNjARS1hseNG4fZs2fjzp07AIqS4cGDB5GZmanwmGLp6elcD2wAOH/+PKRSKde5LCEhAQ0aNOD2JyQk4Pbt20rFp8yxxsbG3OX6knTq1Al+fn6QSqXctuPHj6NDhw7Q0dFRKhZl7d69GxMmTMCRI0dkHj/++CMOHjxY7pY5IZWJWtKEaDEvLy80adIE7du3R25uLjZu3IgxY8YAAJYtW4YrV66gbdu2mDZtGgwNDXH37l28efNG5n6zkZERJkyYgNmzZ0MqlWLdunWYPXs2nJycAABDhw7Fzz//DH19fejq6mLr1q3Q1VXuT4Myx7Zq1Qrv3r3DqlWr4O7urvBy8tKlS3H8+HEMGDAAXl5eCAwMxN9//y3z5UIdwsPDcevWLaxevVpuX//+/ZGdnY2///4bI0eOVOt5CVEVJWlCtJi/vz8OHjyIwMBA6OnpYcuWLdwkIJaWlggKCsKBAwdw9+5d6Ovro1evXhg9erRMHbVr18bx48dx5MgRxMXFYePGjTIdwBYvXozatWvjxo0bEAqF2LJlCyIiImTGKNevXx8DBw6Ui0/ZY8+fP4/Tp0/j6dOnMDQ0lJvMxNbWFiEhIfjzzz/h7+8POzs7PHz4UGbMtaIY7OzsMGrUqFLfw5SUFBw5cgStW7fGixcvMH78eHzyySdy5YyNjeHj44PExEQA/00LGhsbi9q1a5d6DkIqC4+V9+YRIaTa+PXXX/Hrr78iLCxM06FoxKlTp7ge4lOmTEGfPn2UPvbu3bvYtGkTAKBv376YPHlypcRISGkoSRNSg33sSZqQ6o46jhFSg5V0mZoQUj1QS5oQQgjRUtSSJoQQQrQUJWlCCCFES1GSJoQQQrQUJWlCCCFES1GSJoQQQrQUJWlCCCFES1GSJoQQQrQUJWlCCCFES1GSJoQQQrTU/wFbGxlvbADi/wAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m = model\n", "seps = np.linspace(2.5, 12.0, 25)\n", "e_elec = []\n", "for d in seps:\n", " pair = {\"pos\": np.array([[0, 0, 0], [d, 0, 0]]), \"atomic_numbers\": [9, 3]} # F, Li\n", " gp = structure_to_graph(pair, 5.0)\n", " e_elec.append(float(m(gp)[\"energy_elec\"][0]))\n", "import matplotlib.pyplot as plt\n", "plt.figure(figsize=(5, 3.2))\n", "plt.axvline(5.0, ls=\"--\", c=\"gray\", label=\"GET cutoff\")\n", "plt.plot(seps, e_elec, \"o-\", ms=3)\n", "plt.xlabel(\"separation [Å]\"); plt.ylabel(\"electrostatic energy [kcal/mol]\")\n", "plt.title(\"BAMBOO electrostatics acts beyond the cutoff\"); plt.legend()\n", "plt.tight_layout(); plt.savefig(\"bamboo_electrostatic_tail.png\", dpi=110); plt.show()" ] }, { "cell_type": "markdown", "id": "c7c94cd3", "metadata": {}, "source": [ "## 5. Deploy as an ASE calculator\n", "\n", "`XNNCalculator` wraps any xnn model for ASE (energies, forces, and, for\n", "periodic cells, stress). BAMBOO works in kcal/mol and Å." ] }, { "cell_type": "code", "execution_count": 6, "id": "7d3ba5b6", "metadata": { "execution": { "iopub.execute_input": "2026-07-20T04:24:43.287765Z", "iopub.status.busy": "2026-07-20T04:24:43.287686Z", "iopub.status.idle": "2026-07-20T04:24:43.298674Z", "shell.execute_reply": "2026-07-20T04:24:43.298083Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "energy [kcal/mol]: -0.3304\n", "forces shape : (3, 3)\n" ] } ], "source": [ "from ase import Atoms\n", "from xnn.common.deploy import XNNCalculator\n", "\n", "atoms = Atoms(\"OH2\", positions=[[0, 0, 0], [0.96, 0, 0], [-0.24, 0.93, 0]])\n", "atoms.calc = XNNCalculator(ForceStressOutput(model), cutoff=5.0)\n", "print(\"energy [kcal/mol]:\", round(atoms.get_potential_energy(), 4))\n", "print(\"forces shape :\", atoms.get_forces().shape)" ] }, { "cell_type": "markdown", "id": "4ccf986c", "metadata": {}, "source": [ "## Summary\n", "\n", "BAMBOO is a single model that predicts a physically-decomposed energy\n", "(`energy_nn` + `energy_elec`), conserved per-atom partial charges, and the\n", "molecular dipole, with correct E(3) symmetries and a genuinely long-range\n", "electrostatic term, all through the standard xnn `AtomicGraph` /\n", "`ForceStressOutput` / `XNNCalculator` interfaces." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.12" } }, "nbformat": 4, "nbformat_minor": 5 }