{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6845c633",
   "metadata": {},
   "source": [
    "# DFT-D4 dispersion, block by block: reproducing the reference `dftd4` with `xnn`\n",
    "\n",
    "[DFT-D4 (Caldeweyher *et al.*, *J. Chem. Phys.* **150**, 154122, 2019)](https://doi.org/10.1063/1.5090222)\n",
    "is the charge-dependent successor of the D3 London-dispersion correction: atomic\n",
    "dynamic polarizabilities of tabulated reference systems are **scaled by atomic\n",
    "partial charges** (eq 2, classical EEQ charges by default), **interpolated in the\n",
    "coordination number** (eq 8), integrated to pairwise $C_6$ coefficients\n",
    "(Casimir–Polder, eq 9), and summed with BJ damping (eqs 18–21) plus an\n",
    "Axilrod–Teller–Muto three-body term (eqs 22–27).\n",
    "\n",
    "In xnn this is `xnn.common.models.D4Dispersion`, a **model-agnostic add-on** (like\n",
    "`LatentEwald` / `ForceStressOutput`): it stands alone as the model `\"d4\"` or wraps any\n",
    "short-range model (`extra: {dispersion: {...}}`), and deploys through every channel\n",
    "(PyTorch, TorchScript, ASE, LAMMPS). The implementation is an independent one written\n",
    "from the paper; this notebook checks every block against the **reference\n",
    "implementation** [`dftd4`](https://github.com/dftd4/dftd4) (v4.2.0, through its Python\n",
    "package), the same protocol as the other fidelity notebooks.\n",
    "\n",
    "**Reference code.** xnn does not copy, vendor, link or import any of `dftd4`'s code: the\n",
    "package is used here purely as an external oracle whose numbers we compare against,\n",
    "exactly as one would use a compiled reference program. The *reference data* of the\n",
    "method (TD-DFT polarizabilities, reference CNs and charges, element constants) are\n",
    "numerical values extracted from the published sources by `tools/build_d4_reference.py`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71794685",
   "metadata": {},
   "source": [
    "## 0. Setup: `float64`, and `dftd4` *before* `torch`\n",
    "\n",
    "The `dftd4` wheel bundles its own OpenMP runtime. If `torch` is imported first, the\n",
    "two runtimes clash and `dftd4` returns **wrong EEQ charges** (silently: the water\n",
    "molecule comes out with $q_\\mathrm{O}\\approx -0.009$ instead of $-0.586$). Importing\n",
    "`dftd4` first avoids it (so does `torch.set_num_threads(1)`, which the test suite\n",
    "uses); we assert the known water charges as a guard.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8b276907",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:33.118796Z",
     "iopub.status.busy": "2026-09-23T19:19:33.118609Z",
     "iopub.status.idle": "2026-09-23T19:19:36.313442Z",
     "shell.execute_reply": "2026-09-23T19:19:36.312912Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xnn: 0.2.1 | dftd4 (reference) 4.2.0 | PBE0-D4 parameters: {'s6': 1.0, 's8': 1.20065498, 'a1': 0.40085597, 'a2': 5.02928789, 's9': 1.0, 'alp': 16.0}\n"
     ]
    }
   ],
   "source": [
    "# dftd4 (the reference) must be imported before torch, see above\n",
    "import numpy as np\n",
    "from dftd4.interface import DampingParam, DispersionModel\n",
    "\n",
    "import logging, warnings\n",
    "logging.disable(logging.WARNING)\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "from ase.build import molecule, bulk\n",
    "\n",
    "torch.set_default_dtype(torch.float64)\n",
    "\n",
    "import xnn\n",
    "from xnn.common.data import structure_to_graph, collate\n",
    "from xnn.common.models import D4Dispersion, DFTD4, ForceStressOutput, c6_matrix\n",
    "from xnn.common.models.d4 import BOHR, HARTREE, PBE0_D4\n",
    "\n",
    "probe = DispersionModel(np.array([8, 1, 1]), molecule(\"H2O\").positions / BOHR).get_properties()\n",
    "assert abs(probe[\"partial charges\"][0] + 0.58639069) < 1e-6, \"dftd4/torch OpenMP clash: restart the kernel\"\n",
    "print(\"xnn:\", xnn.__version__, \"| dftd4 (reference) 4.2.0 | PBE0-D4 parameters:\", PBE0_D4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c63d50",
   "metadata": {},
   "source": [
    "## Helpers\n",
    "\n",
    "Both codes are driven in **atomic units**: positions in bohr go to `dftd4`, and xnn's\n",
    "Angstrom/eV results are converted back with the same CODATA-2018 factors it uses\n",
    "internally. `dftd4` reports the gradient $\\partial E/\\partial \\mathbf r$ and the strain\n",
    "derivative $\\partial E/\\partial\\epsilon$ (its \"virial\"); xnn gives forces and the stress\n",
    "$\\sigma = V^{-1}\\,\\partial E/\\partial\\epsilon$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0ecc15e6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:36.317904Z",
     "iopub.status.busy": "2026-09-23T19:19:36.317813Z",
     "iopub.status.idle": "2026-09-23T19:19:36.325128Z",
     "shell.execute_reply": "2026-09-23T19:19:36.324543Z"
    }
   },
   "outputs": [],
   "source": [
    "def reference(atoms, charge=0.0, params=PBE0_D4, cutoffs=None):\n",
    "    periodic = bool(atoms.pbc.any())\n",
    "    model = DispersionModel(atoms.numbers, atoms.positions / BOHR, charge=charge,\n",
    "                            lattice=atoms.cell[:] / BOHR if periodic else None,\n",
    "                            periodic=np.array(atoms.pbc) if periodic else None)\n",
    "    if cutoffs:\n",
    "        model.set_realspace_cutoff(**cutoffs)\n",
    "    res = model.get_dispersion(DampingParam(**params), grad=True)\n",
    "    res.update(model.get_properties())\n",
    "    return res\n",
    "\n",
    "def xnn_d4(atoms, charge=0.0, model=None, **options):\n",
    "    model = model or D4Dispersion(**options)\n",
    "    s = {\"pos\": atoms.positions, \"atomic_numbers\": atoms.numbers, \"total_charge\": charge}\n",
    "    if atoms.pbc.any():\n",
    "        s[\"cell\"], s[\"pbc\"] = atoms.cell[:], atoms.pbc\n",
    "    graph = structure_to_graph(s, model.cutoff)\n",
    "    out = ForceStressOutput(model, compute_stress=bool(atoms.pbc.any()))(graph)\n",
    "    out = {k: v.detach() for k, v in out.items()}\n",
    "    out[\"energy_au\"] = float(out[\"energy\"]) / HARTREE\n",
    "    out[\"gradient_au\"] = -out[\"forces\"].numpy() / HARTREE * BOHR\n",
    "    if atoms.pbc.any():\n",
    "        out[\"virial_au\"] = out[\"stress\"][0].numpy() / HARTREE * BOHR**3 * (atoms.get_volume() / BOHR**3)\n",
    "    out[\"c6\"] = c6_matrix(out[\"dynamic_polarizabilities\"]).numpy()\n",
    "    return out\n",
    "\n",
    "def report(name, atoms, charge=0.0, **options):\n",
    "    ref, mine = reference(atoms, charge), xnn_d4(atoms, charge, **options)\n",
    "    row = {\"system\": name, \"N\": len(atoms), \"E [Eh]\": ref[\"energy\"],\n",
    "           \"|dE|\": abs(mine[\"energy_au\"] - ref[\"energy\"]),\n",
    "           \"|dgrad|\": np.abs(mine[\"gradient_au\"] - ref[\"gradient\"]).max(),\n",
    "           \"|dCN|\": np.abs(mine[\"coordination_numbers\"].numpy() - ref[\"coordination numbers\"]).max(),\n",
    "           \"|dq|\": np.abs(mine[\"eeq_charges\"].numpy() - ref[\"partial charges\"]).max(),\n",
    "           \"|dalpha|/alpha\": np.abs(mine[\"polarizabilities\"].numpy() / ref[\"polarizabilities\"] - 1).max(),\n",
    "           \"|dC6|/C6\": np.abs(mine[\"c6\"] / ref[\"c6 coefficients\"] - 1).max()}\n",
    "    if atoms.pbc.any():\n",
    "        row[\"|dvirial|\"] = np.abs(mine[\"virial_au\"] - ref[\"virial\"]).max()\n",
    "    return row\n",
    "\n",
    "def show(rows):\n",
    "    keys = list(rows[0])\n",
    "    for extra in (\"|dvirial|\",):\n",
    "        if any(extra in r for r in rows) and extra not in keys:\n",
    "            keys.append(extra)\n",
    "    print(\"\".join(f\"{k:>16s}\" for k in keys))\n",
    "    for r in rows:\n",
    "        print(\"\".join(f\"{r.get(k, ''):>16.2e}\" if isinstance(r.get(k, ''), float) else f\"{str(r.get(k, '')):>16s}\" for k in keys))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c0a153f",
   "metadata": {},
   "source": [
    "## Block 1: Coordination numbers and EEQ charges · eqs 6, 11–16\n",
    "\n",
    "The D4 coordination number weights each error-function count by an\n",
    "electronegativity factor (eq 6); the EEQ charges solve the constrained linear\n",
    "system of eq 16 with the plain (electronegativity-free) CN of eq 14, softly capped at\n",
    "8 in the reference code. `dftd4` exposes both through `get_properties()`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dc08b9a9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:36.326427Z",
     "iopub.status.busy": "2026-09-23T19:19:36.326361Z",
     "iopub.status.idle": "2026-09-23T19:19:39.005989Z",
     "shell.execute_reply": "2026-09-23T19:19:39.005347Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          system               N          E [Eh]            |dE|         |dgrad|           |dCN|            |dq|  |dalpha|/alpha        |dC6|/C6\n",
      "             H2O               3-0.0001576174317030107        5.42e-20        5.42e-20        0.00e+00        3.33e-16        2.22e-16        2.22e-16\n",
      "           CH3OH               6-0.0012335328614920178        2.17e-19        1.29e-19        0.00e+00        3.75e-16        6.66e-16        1.55e-15\n",
      "            C6H6              12-0.00898814444745347        3.47e-18        3.52e-19        4.44e-16        8.33e-17        6.66e-16        1.33e-15\n",
      "        CH3CH2OH               9-0.0030075001576794188        4.34e-19        9.15e-20        4.44e-16        1.67e-16        6.66e-16        1.11e-15\n",
      "             CO2               3-0.0007974225976500424        4.34e-19        1.02e-20        0.00e+00        1.67e-16        3.33e-16        9.99e-16\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "water CN     dftd4: [1.60884227 0.80442113 0.80442113] \n",
      "             xnn  : [1.60884227 0.80442113 0.80442113]\n",
      "water q      dftd4: [-0.58639069  0.29319534  0.29319534] \n",
      "             xnn  : [-0.58639069  0.29319534  0.29319534]\n"
     ]
    }
   ],
   "source": [
    "rows = []\n",
    "for name, charge in [(\"H2O\", 0.0), (\"CH3OH\", 0.0), (\"C6H6\", 0.0), (\"CH3CH2OH\", 0.0), (\"CO2\", 0.0)]:\n",
    "    rows.append(report(name, molecule(name), charge))\n",
    "show(rows)\n",
    "\n",
    "water = molecule(\"H2O\")\n",
    "ref, mine = reference(water), xnn_d4(water)\n",
    "print(\"\\nwater CN     dftd4:\", ref[\"coordination numbers\"], \"\\n             xnn  :\", mine[\"coordination_numbers\"].numpy())\n",
    "print(\"water q      dftd4:\", ref[\"partial charges\"], \"\\n             xnn  :\", mine[\"eeq_charges\"].numpy())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "037ef80d",
   "metadata": {},
   "source": [
    "## Block 2: Charge scaling, reference weighting and $C_6$ · eqs 2–9\n",
    "\n",
    "The charge-dependence is the point of D4. Removing or adding an electron to the\n",
    "ammonium/amide pair changes the EEQ charges, the scaled polarizabilities (eq 4) and\n",
    "hence the $C_6$ coefficients; both codes agree for every total charge. The static\n",
    "polarizabilities $\\alpha(0)$ and the full $C_6$ matrix are compared element by element.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "787e1809",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:39.007872Z",
     "iopub.status.busy": "2026-09-23T19:19:39.007785Z",
     "iopub.status.idle": "2026-09-23T19:19:42.949154Z",
     "shell.execute_reply": "2026-09-23T19:19:42.948511Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          system               N          E [Eh]            |dE|         |dgrad|           |dCN|            |dq|  |dalpha|/alpha        |dC6|/C6\n",
      "             NH3               4-0.000336297649224461        5.42e-20        1.02e-19        0.00e+00        1.11e-16        0.00e+00        3.33e-16\n",
      "     NH4+ (q=+1)               5-0.0003545463401048189        2.71e-19        8.13e-20        0.00e+00        8.88e-16        6.66e-16        1.11e-15\n",
      "     NH4- (q=-1)               5-0.0013716042669924796        1.30e-18        8.67e-19        0.00e+00        8.88e-16        6.66e-16        1.33e-15\n",
      "     NH3+ (q=+1)               4-0.0001816518458468447        2.71e-20        1.63e-20        0.00e+00        1.67e-16        0.00e+00        3.33e-16\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6sAAAFSCAYAAAAU4El9AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAApHlJREFUeJzs3XdYVEfbBvB7lw7SwYIiKCKIiAWxYe89xh57Yo0ae+yx9957bzH2EkuMxt6xK1YEGyIiTTrszvcHL/u5UgQpZ4H7d11cyZ4zM/vsLjL7nJkzIxNCCBARERERERFpELnUARARERERERF9jckqERERERERaRwmq0RERERERKRxmKwSERERERGRxmGySkRERERERBqHySoRERERERFpHCarREREREREpHGYrBIREREREZHGYbJKRJQHvH//Hl5eXlAqlVKHkms8ffoUT58+zZK2Hj16hBcvXmRJW5T33L59G2/evEl2XAiBd+/e4datW/Dz8/uutj9//gwvLy+Eh4dnMkoiIs3DZJWIco3g4GB4eXnhzp07KZ5/9eoVvLy8UqwTGhqaYp3Xr1/Dy8sLCQkJqT5vXFwcvLy84OXlBYVC8d3xf+3Dhw/ffO70lt+4cSM8PDwQFRWVZfFll8ePH+PZs2eSt9G7d2/0798/U20k6dSpE4YNG5YlbVHeU6NGDSxcuFDt2Nu3b+Hu7o4KFSqgf//+2Lp1K4KCguDl5YXY2Nh0t33r1i14eHjgypUrqZZRKpW4desWvLy8EBcX992vg4gopzFZJaJc4/jx4/Dw8EClSpVw6tSpZOenT58ODw8PtWQuqc65c+dSbHPWrFnw8PBINZkFgEmTJsHDwwMeHh74/PlzZl+Gyvbt27/53Jkpr6m6d++OgQMHSt4GkZRmzJiBDx8+qC6YTZ48GYcOHYKHh0eKo7CZMW/ePFSuXBkeHh7w9/fP0raJiLITk1UiynW0tLQwevToHJnyevnyZSxcuBCurq7Z/lxElH/cvXsXrq6uMDAwyNbnuXfvHiZPngw3N7dsfR4iouygLXUAREQZ1bt3b6xbtw7bt29Hz549s+15IiIi0KNHD/Tv3x8xMTF4+PBhuuqFhITAx8cHACCTyWBkZISSJUtCV1dXVeb169d4+/YtgMQvk6ampgAAZ2dnFChQIFmbGS3/6tUrREdHo3Tp0pDL1a9Lvn//Hu/evUOlSpUgl8vx5s0bBAUFoWLFiqoyCQkJePnyJeLj41GqVCno6eml+FrTWy7JvXv3EBUVBS0tLdWUbV1d3WRfpCMjI+Hr6wt9fX04ODhAJpOlu42XL18iODgYQOKFDQsLC9jZ2aUZV3oolUq8fPkSQgg4ODgke1+TxMXF4cWLFzAzM4ONjU2y8+mN79GjR9DT00OpUqUQFxeHly9fwsDAQK1s0n2ySe/R7du3YW1tDVtb22TthYeH4/Xr1zA0NETJkiVTjD0iIgJv376FiYlJirH7+PggJCQEFStWhJaWVoptAIC3tzdiYmJQqVKlZO0/efIEhQoVgq2tLWJjY/HgwQPY29vDysoKoaGhePPmDezt7WFsbKxWNyNlv5ZUt1ixYihcuHCy848fP4YQAi4uLul+L9Ly9efypbCwMDx//hwBAQHQ0dFR/Q6bmZnh1atXAICHDx+qZlA4Ojqq/r0nvZYXL17A0tIyxdfy9evu3r07OnTogMKFC+P+/fsZeh1ERJITRES5xPbt2wUA8e+//4qqVasKW1tbER0drTrfu3dvAUDEx8cnq3Pw4MEU2+zfv78AID5+/JjsXJ8+fYStra0IDw9XtR0SEvLNOM+cOSPc3d2Fu7u7qFSpkrCxsREGBgZixIgRQqlUCiGEmDdvnihWrJgAIMqXL68qf/fu3RTb/Fb56dOnCwDCx8dH1K1bV5QqVUoYGRmJ4sWLCy8vL7W2kso+efJE1KxZU5QsWVKYmpoKIYRQKBRi2rRpwtTUVBQsWFDY29uLAgUKiJkzZ6q1kd5yX2vQoIEwNDQUxsbGqtfQvHlz1fnIyEjRp08foaurK+zt7YWpqakoUqSI2L59e7rbGD16tOp4uXLlhLGxsShevLj4+++/1WLx9PQUderUSTNeIYRISEgQ06ZNE+bm5sLMzEw4OjoKKysrsWzZMlWZsmXLihYtWogTJ06IkiVLCkdHRyGXy0XLli1FZGSkWnvpjS+pzUOHDgl7e3tha2srevfuLYQQ4vr168LR0VHo6ekJR0dH4ebmJp4+fSr09PTE0KFD1doJCAgQbdu2Fbq6usLBwUFYWFgIBwcHcfr0aVWZ6Oho0a1bN6GrqytKly6ter5t27aptfXDDz+k69/Bzp07BQCxcOFC1TGlUilat24tjI2NxZMnT4QQQjx//lwAEOvXrxfjx48Xtra2omjRokJXV1fMnz9frc2MlP1aQkKCKFKkiNrvSZLQ0FBhYGAgBg4cmKH3IiVeXl7Cyckpzc/l4sWLwt3dXejr6wtTU1PV78LEiROFnZ2dACBcXV1Vxy9evKhqf8OGDcLU1FSYm5sLOzs70aFDB3H8+HEBQJw4cSJZPKNGjRJWVlbi48ePYuTIkQKA8PX1/ebrICLSFExWiSjX+DJZvXjxogAg5syZozqfVrK6YMECcfPmzWQ/bdu2TTFZ/fvvvwUAcezYMbW205OspuTEiRNCX19frF69WnVs/vz5qSbKKUmrfFIC2qlTJ1UiEBISIlxdXUWlSpVSLNumTRvh7e0thEhMfoQQYsiQIUJXV1ccOHBAVf6ff/4Rurq6YvHixapj6S2XEnd3d9GgQYMUz7Vr104YGxuLs2fPCiESk+JRo0YJAGrPlVYbX4uOjhbDhg0TBQoUEK9fv1YdT2+yOnDgQKGrq6uWMIeEhIhp06apHpctW1a4urqKPn36qC6gXLhwQcjlcjF37tzviq9s2bLCxcVFdO/eXdXm9evXxYcPH4SFhYWoW7euCA0NFUII8ebNG9GxY0ehq6urlqxGRkaKMmXKCBcXF/Hs2TMhROJ7OmLECKGvr6/6/KdMmSKMjY3F06dPVXX9/f3FrFmz1GIdOXKkcHd3F58/f/7m+/brr78KbW1tVbI1e/ZsAUDs2bNHVSYpAa1WrZrYsmWL6vioUaOEXC5XxZfRsikZO3as0NLSEm/fvlU7vnr1agFA3Lp1K0PvxdcCAwOFlZWVqF27turvhJ+fn+pCwdcXEZycnESLFi3Ujq1fv14AEM+fP0/WflJSOmbMGNVFrxMnTogaNWqkmKwm/f7t2rVLCCGYrBJRrsRklYhyjS+TVSESR3lMTU1FUFCQECLtZNXBwUE1UvHlj7W1dbIE8OPHj6Jw4cKia9euqmPfk6xGRkaKR48eCS8vL3Hz5k1RuXJl0bRpU9X57EhWd+7cqXZ8zZo1AoB48eJFsrLr169XK/vq1SuhpaUlRo0alaz9fv36iUKFCgmlUpnucqlJLdF8/PixACDGjx+vdjwhIUE4ODiIihUrfrONLwUEBIi7d++KmzdvilOnTgkAaklOepLVly9fCrlcLoYNG5ZmubJly4oCBQqoksckdevWFZUrV/6u+MqWLSv09fWTfd6zZs0SAMTjx4/VjiddYPkyKVqxYoUAIC5duqRWNi4uThQtWlT069dPCCFE586dRcmSJdP83DIqNjZWeHh4CBsbG7F7926hpaWVLGFLSkBbtWqldjwkJETIZDIxY8aM7yqbkufPnwuZTJZs9L9y5cqiQoUKqsff+14kJeNfJ82HDh1K9rkIkfFktU6dOqJEiRIiISFB7fiAAQOSJavh4eGiRIkSomXLlqpjTFaJKDfiPatElGvNnTsXrq6umD59OpYsWZJm2QULFqBNmzbJjg8YMABr165NdkyhUHyzzdQEBQVhwIABOHLkCAoXLgxLS0toaWnBx8cHkZGR39Vmenl6eqo9dnBwAJB4D2vS/yepXbu22uMLFy5AoVDAzs4Od+/ehUi8oAkAsLS0xIcPH/D27dt0l0vpvsm0XL9+HQBQt25dteNaWlqoXbs2tmzZgpiYGOjr66fZzrlz5/Dbb7/h2bNnsLe3V7unN+mewPS6dOkSlEolmjZt+s2yrq6uavcWAonv/5EjR747vjJlysDKykrt2PXr12FlZQVnZ2e1419/9gBw9uxZ6OnpwcDAQO2zEkKgePHiqvslmzVrht27d6NOnTro2rUr6tSpk6z9jNLV1cXevXtRqVIldO7cGdWrV8f8+fNTLPt17GZmZrCwsEjx88pI2S+VKlUKtWvXxubNmzFu3DjIZDI8ePAAXl5eWL58uarc974X165dg5WVFcqUKaN2/Ot/Z9/r+vXraN++fbJ7hWvXro01a9aoHRs2bBg+ffqE1atXZ8lzExFJhckqEeVaTk5O6NOnD1avXo0hQ4ZkSZsPHjzA/v37MWbMGPj5+cHPzw9AYgIKJK7gaW1tjbJly6baRu/evXHt2jXcvXtXbcGWhg0bqhZJyi6WlpZqj5MSu5T2X/16cZbw8HAAwOrVq7Fly5Zk5d3d3REfH5/uchkVEREBADA3N092ztzcHEIIREVFpZmsfvz4Ea1atUKjRo1w9epVVSIYFBQEa2vrDK8gnbRVkYWFxTfLfv3eA4nv/5fvfUbjS2kBncjIyGRJMQCYmJgkW/Qp6bPq169fijEnLdbUo0cPWFpaYvPmzZg4cSKCgoJgZ2eHGTNmoFu3bmm97DSZm5vDwsICwcHBqF69OnR0dFIsl5737nvKfq1Pnz7o3r07Lly4gDp16mDjxo3Q09ND165dVWW+973IyOeSUQqFAjExMSm2b2Zmpvb49evX2LRpEwYNGoSAgAAEBAQASNynGUj8GxceHs7VgYkoV2CySkS52pQpU7Bjxw6MGzfumyuCpoe2tjbc3d1x+vRpnD59WnU8adRm+PDhcHJywu7du1Nt4+TJk+jTp49aoiqEwLNnz2BoaJjpGLPK11+gk0ZCZ8+ejdatW6daL73lMqpYsWIAEt/rypUrq53z8/NDgQIFkn0x/9qlS5cQERGB4cOHq41YPnny5LtiKl68OIDEVXA9PDy+q43MxJdSklO0aFFcvnwZSqVS7fy7d++SJbtJn9WVK1fUVqNOSYsWLdCiRQsAiQnN6NGj0bNnT1SqVEntdzkjevfuDX9/f/Ts2ROLFy9GvXr10LJly+9qKyu0a9cOv/32GzZt2oTq1atj586daNu2bbILJN/zXhQtWhSXLl2CQqFQG/18+/ZtprfZ0tLSQqFChfD69etk574+JpfL4e7ujmvXruHatWtqcQDAhAkTYGtri2PHjmUqJiKinMB9VokoVytUqBBGjRqFPXv2qKY0ZkaZMmXg5eWV7OeHH34AkDitMq1EFUgc6UjadiLJwYMH4e/vn6wckPKoZ2rtZqR8RjVo0ACFChXCqlWrUjyfNIU5veVSY2ZmluJrqFu3LoyNjbFp0ya14+/fv8fx48fRsmVLVXKWWhtJ79HX7//KlSvTjCk19evXR8GCBbFs2bJkCUd0dHSG28uK+Fq2bIno6Gjs27dP7fiOHTuSbZPStWtXxMbGYuPGjSm2lfRZff2ZlStXDoMHD4ZSqVRtwwQkJu1eXl5QKBTfjHPx4sXYt28f1q5di02bNqFBgwbo0aOHaraCFAwMDNClSxfs27cP27ZtQ1BQEHr37q1WJr3vxddatWqFmJgY7N27V+341q1bk30uqUnr33irVq3w77//qkZKgcSLYDt27FArV6xYsRT/hiWNCh85coSJKhHlGhxZJaJcb9SoUVi7di3u3bsndSgAEqcazp07F25ubqhWrRpu3ryJQ4cOoXHjxnj58qWqXPXq1SGXy7FgwQL89NNP0NHRSXXf1O8pn1GGhob4888/0aZNGzRo0AD9+vWDjY0N3r59i3PnzuHZs2c4e/ZsusulxtPTE3PnzsXu3bvh4OAAPT09uLm5wdTUFMuXL8cvv/yCbt26oUePHvj06RMmT54MKysrtfsdU2ujRo0aKFOmDIYNG4bY2FgYGxtj+/btye7Xzch7smPHDrRp0wb16tXDwIEDUbBgQdy6dQt79uzBjRs3MtReVsTXrl07NG/eHH369MGHDx/g6uqKixcvIiwsDLq6umqJUf369TFx4kQMHToUjx8/RuPGjWFkZITnz59jz549aNSoEcaMGYOffvoJ5ubmaNy4Mezt7REQEICZM2eiZMmSavcQjxw5EocPH0ZISEiao9xXrlzBmDFjMGDAAFWStGvXLlSsWBEdOnTA5cuXvznSm1369OmDVatWYdiwYShRogTq16+vdj6978XX2rZtixYtWqBv37748OED3NzccO7cOYSHh6f7tXp4eEBXVxdLlixBnz59oKurq9pnddq0aThx4gTq1auHKVOmwNTUFJs2bUKVKlVw6dKlzLwlREQai8kqEeUalpaWcHd3h4mJidpxIyMjzJs3T7Ug0pdf1pPqpHQfJJB4z567u3uq99Ilsbe3h7u7O7S1v/1nc9q0abC1tcXRo0dx4sQJ1KhRA8eOHcP06dPVpgGXLVsWf/75J3bs2IFhw4ZBoVBg48aNKF++fIrtplXexsYG7u7uyRZfMTY2hru7u1pikVpZAKhXrx4ePXqEdevWYdOmTYiKioKdnR3q16+PZcuWZbhcSsaMGQOFQoG1a9ciIiICBQsWVI309OzZE05OTli3bh2mTp0KfX19dO7cGUOGDFFbaCi1NvT09HDu3DnMmTMHK1asQIECBdC5c2d07NgRJ0+ehI2NjaoNZ2fnFN+DrzVq1Aj379/H6tWrsWbNGtVU8b///ltVxtXVFUWKFElW187ODpUqVVI9zkh8qbUpk8lw8OBBLF++HIcOHcKhQ4fQunVrjBkzBosXL4aBgYFa+enTp6Np06bYtm0b5s+fDx0dHTg6OmLixImq5Gv//v3Yu3cvTpw4AT8/P1hYWKBz587o27ev2vT6UqVKffPfQVxcHObNm4fmzZurLVJmbW2Nv/76C8OGDcPmzZvRv39/6Ovrw93dHdbW1snaKV++PEqWLKl6nJGyaalYsSLatWsHPz8//Pzzz8lGPdP7XnxNJpPhwIEDWL58OY4cOYIjR46gZcuWmDx5Mm7cuKGaUp7E1dU12SJkdnZ22LdvHzZt2oThw4erFnqrWbMmihQpgps3b2LOnDlYvnw5LC0t8euvv8La2hrnz59P8X7WL9na2sLd3R16enrpep+IiDSBTCQt4UhERES51osXL+Do6Ig1a9agf//+UodDRESUabxnlYiIKJeJjY1Ndmz9+vWQy+Vo3LixBBERERFlPU4DJiIiymXGjh0LLS0t1KtXD1paWjh+/DhWrFiBsWPHokSJElKHR0RElCU4DZiIiCiXiYmJwZo1a3Dx4kUEBASgWLFi6NSpE9q2bSt1aERERFmGySoRERERERFpHN6zSkRERERERBqHySoRERERERFpHCarREREREREpHGYrBIREREREZHGYbJKREREREREGofJKhEREREREWkcJqtERERERESkcZisEhERERERkcZhskpEREREREQah8kqERERERERaRwmq0RERERERKRxmKwSERERERGRxmGySkRERERERBqHySoRERERERFpHCarREREREREpHGYrBIREREREZHGYbJKREREREREGofJKhEREREREWkcJqtERERERESkcZisEhERERERkcZhskpEREREREQaR1vqALKbUqmEv78/jI2NIZPJpA6HiIhyOSEEAMDExCRf9SvsT4mIKKsIIfD582fY2NhALk99/DTPJ6v+/v6wtbWVOgwiIspjwsLCYGJiInUYOYb9KRERZbU3b96gWLFiqZ7P88mqsbExgMQ3Ij99qSAiouwRHh6eL5M29qdERJRVkvrSpL4lNXk+WU2aqmRiYsLOlYiI6DuxPyUioqz2rdtKuMASERERERERaRwmq0RERERERKRxmKwSERERERGRxsnz96yml0KhQHx8vNRhUDbT0dGBlpaW1GEQEeVJSqUScXFxUodBOYD9KRHlhHyfrAohEBAQgNDQUKlDoRxiZmaGwoULc59AIqIsFBcXB19fXyiVSqlDoRzC/pSIslu+T1aTEtWCBQvC0NCQf3DzMCEEoqKiEBgYCAAoUqSIxBEREeUNQgi8f/8eWlpasLW1TXODd8r92J8SUU7J18mqQqFQJaqWlpZSh0M5wMDAAAAQGBiIggULcgoTEVEWSEhIQFRUFGxsbGBoaCh1OJQD2J8SUU7I15c+k+5RZceavyR93rxHmSh/ePPmDfr164eQkBCpQ8mzFAoFAEBXV1fiSCgnsT8lyj+USiV27NiBVatW5ejz5utkNQmn/uYv/LyJ8oekeyePHTuG6OhoJCQkSBxR3se/r/kLP2+i/EGpVCI2NhZr165FwYIFc/S5mawSEVGeIoTAnj174OzsDF9fXwwYMADbt2+HtbW11KERERHlGu/fv0evXr0wcOBAGBgY4OLFi2jfvn2OxsBkNR8KDg7GlClTMrwC8vfWyyoBAQFYvHgx5s2bJ8nzE5Hmu3PnDurUqYOhQ4diwoQJsLOzkzokysMOHDiAQ4cO5Vi9rCCEwD///IOpU6fixo0bksRARJotNjYWc+bMQenSpREdHY1x48ZJFku+XmApK7wLjUZIZOp7ypkb6aKomUEORvRtwcHBmDp1Knr16gUzM7Nsr5cVfv75Z5w+fRpFixbFixcvMHr06Bx9fiLSbPHx8dDR0cGZM2fg6emJY8eOwdjYWOqwKJ1yY18KJCad2traaNOmTY7Uy6xLly7h559/hr29PS5evAhLS0tUqVIlR2MgIs0lhEBCQgKio6Nx5swZHDt2DLVr15Y0JiarmfAuNBr1F5xDbELqe8rpacvx36i6Wd7JCiFw+PBhPHv2DI6OjqhXrx6WLFmCgQMHpjmXPCYmBsuWLQMALFmyBGZmZihVqhS6desGADh37hyuX78OQ0NDtGjRAiVLlkxXvezWsWNHrFu3DitXrsSMGTNy5DmJSDrpTV7i4uKwYsUKLF68GPfu3cOoUaNyMErKClL2pQDw+vVrHDhwAPr6+qhTpw58fX0RERGBjh07plnv5MmTuH//PuRyOaZMmQIA6NOnD4oVK4agoCAcOnQIHz9+hLOzM1q3bq1aLTetetnN3NwcJ0+ehIODQ45fdCYiaaS3P3348CGGDRuGatWqYcaMGfj3339zMMrUMVnNhJDIuDQ7VwCITVAiJDIuyzvYLl264L///kPHjh1x6dIlTJw4Ed7e3mjfvv133/jcp08fHDx4EJ07d0ZgYCBGjRqF3bt348cff8xQO7NmzUJcXOr/KAoXLowBAwZkqM1mzZplqDwR5V7pTV7+qBCHGRPHQFtbG5s3b4aFhUUORklZRcq+9NatW6hbty48PT3h5OSEuXPnQkdHB25ubt9MVlNz//591KtXD+XLl0f58uWxfv16LFiwAGfPns3QaslXrlzBqVOn0izTuXNnODs7p7vNsmXLprssEeV+6elPteMjUevzWezevgUDBw7UuIu+TFZzobNnz2Lv3r14+PChqpP69ddf4e3t/c26+vr6GDJkCJYvX45hw4bB3t4eAHD+/Hls2rQJt2/fRoUKFQAAkyZNwuDBg9G0aVMYGBikWI+IKKt9K3lRxscgFvq4dv0GBg8ejAEDBkBbm90ZZdyIESPwww8/YMeOHarHTk5OcHNz+2bdpk2bYseOHdDW1laNkAJAt27dULduXezfvx8AMGHCBDg5OWH16tUYOnRoqvWIiLJaWv2pUCoAoURMbCzevfPH3bt34eTklMMRfht791zoxIkTqF69utrV1N69e2PNmjXf3ebx48dRpUoVVaIKAAMGDMD06dNx7949VKtWLd1tjR8//rvjICJKjTI2EmGXdyPy6SUU7bsOA4aPhmtRU6nDolwqKioKFy5cwOnTp1XH7Ozs0LBhw+9uMzo6GhcuXMDx48dVx6ysrNCuXTucOHECQ4cOTXdbNWrUQI0aNb47FiKi1ET73UXImfUoUKEpTNxbYdG6bXDS0P6UyWou9P79exQuXFjtWJEiRTLVpr+/f7I2ihQpArlcjnfv3mWorcxMA96yZQv8/PwAJI4Cjx07NkPPTUR5U8T9Uwg5vxW61nYo2H4yZNo6UodEudz79+8BIMX+NCQk5Lva9Pf3hxAiWX9qY2ODq1evZqitzEwD9vPzw5YtW1SPmzZtmqGLzkSUNyWEByL4zHrE+N2DmedPMK7QVOqQvolb1+RCRYoUQUBAgNqxpE43PVLaxNvGxiZZGwEBAVAqlShatGiq9YiIspMyJgIAkBAWCMsmg1Gw00zoWttLGxTlCUkJZVb2pzY2NpDJZMnaeP/+vaovTakeEVF2UsZFQygVEIoEaBmYomi/dTCp8iNkWpp/4Zcjq7lQ06ZNsWjRIjx9+lQ1t3zjxo3prm9qmjjM//nzZ9WxZs2aYcGCBbh//77qXp1169ahSJEiKF++fKr1UpKZacC9evX67rpElHckhAci5OxmxAe9RpFfVsCsVs6sPE75h6GhIWrVqoXNmzejQYMGAIBXr17h9OnTaNmyZbraMDU1VUtMDQwMULt2bWzcuBFNmyaOWAQHB2P//v34448/Uq2XksxMA7a3t+f9sEQEpVKJiIdnEHp+K8wb9IORc01YNh0sdVgZwmQ1F6pfvz7atWuH2rVro1OnTvDz84OPj0+661tZWaFcuXIYNGgQ6tati9KlS6Nbt274+eefUa9ePfz000/4+PEjDh06hN27d8PAwCDNejlh165dePbsGa5du4aoqChVJ/ytrXqIKHdRKpVYtXA2/DcsgaFTTRTsOI2jUJRtFi5ciHr16qFp06ZwdnbG4cOHYWtrm+769evXR/fu3TF8+HCYmpqiT58+WLZsGerVq4cGDRrAzc0NR44cgZOTE3799dc06+XE1jX+/v5Yt24dgMQt6Y4fP46goCC4uLh89+rHRKSZ7t27hx4/90Hoi1cwr98Xhk6eUof0XZisZoK5kS70tOXf3F7B3Cj9S9Wn1+7du3Ho0CE8e/YMdevWRbVq1dSmGH3L+fPnsX//frUruxs2bEC3bt1w/fp1ODk5Yfbs2ap9VtOql5OqVavG+26I8iAhBAIDA1GoUCEkKBQo1HkW9Gw0b1VCynpS9qUeHh549OiRap/V48ePY+XKlcmmBqembdu2sLKywo0bNxARkThl3c3NDU+fPsXBgwfx8eNHLFiwQG2f1dTq5TSuCUGUN3348AEFCxaEQqFAzXoNEVKvMuQ6+lKH9d1kQgghdRDZKTw8HKampggLC4OJiYnauZiYGPj6+qJEiRLQ1/++DzG9G+1mt6CgIFhbW+PBgwdwdXXN9ufLzbLicyeirHP79m0MHToUBgYGOHXqVLr3Wf1vVN0c+fv6tbT6lbwsO/tTTelLAWDw4MEICAjAvn37cuT5cjP2p0SaIyYmBosXL8asWbNw8uRJeHp6anR/mt6+lCOrmVTUzECSL0tpOX78OG7cuJHiudKlS6NLly45HBERUXIxMTH47bffsHPnTgwfPhzjxo0DkPh39b9RdTUmedF0GzZsSJZYFSpUCFu3blU79vDhQ6xatQr+/v4oU6YMRowYAWtr65wMNVWa2JdGRUVh3rx5qZ7v168fbGxscjAiIqKUnTlzBn379oWZmRlOnDgBT8/EKb95oT9lsppHGBoaYvLkyShYsCBev34tdThERKmKi4vDx48fYWNjAxsbGzx8+DDZLQeamLxoqidPniAoKAgzZsxQHUtaayDJvXv3UKNGDXTp0gUdO3bExo0bUa1aNdy+fVu1eB4lat68uWRTc4mIMsLPzw/29vbQ1dXFhAkT0KtXL7VbDoDc359yGjCnr+Q7/NyJpHPs2DEMHz4cnp6e2Lx5s9ThfBdNmwY8atQoPHz4ECdPnky1TOvWrREdHY1///0XQOKooa2tLUaNGqUa0f4W9qf0NX7uRNL49OkTJk2ahE2bNuH+/ftwdHSUOqQMS29fyn1WiYgo2wUHB6N58+bo0aMHhg8fjvXr10sdUp7y+PFjtGvXDr169cLatWuRkJCgOqdQKHD69Gm0adNGdczQ0BDNmjVLM8ElIiLNs3PnTjg6OuL169e4d+9erkxUM4LTgImIKNuEhoYiMjIShQsXRq1atbBjxw5YWFhIHVaeoq2tjSZNmqBBgwb48OEDZs2ahR07duDcuXPQ0tLChw8fEB0dnWxLlmLFiuHixYupthsbG4vY2FjV4/Dw8Gx7DURElLZHjx6hbNmysLa2xq5du1R7Oed1HFklIqIsp1AosG7dOjg6OmLt2rXQ0tLCuHHjmKhmg4kTJ2LdunXo1KkThgwZgv/++w83btzAjh07ACTeIwwkv4/V0NBQdS4ls2fPhqmpqeonI/uPEhFR1njx4gXatGmDWrVq4dOnT2jcuHG+SVQBJqtERJTF/Pz84O7ujnnz5mHjxo2YOnWq1CHlaQUKFFB77ODgACcnJ9y+fRsAYGZmBiBxKvaXPn36BHNz81TbHTduHMLCwlQ/b968ydrAiYgoTfPmzYObmxsKFy6Mp0+fwtLSUuqQcpzk04BfvXqFffv2wdzcHL/88kuKZe7evYuzZ89CX18fLVu25NVdIiKpKBXAqytAxAegQCHArgYgT1x50M/PD3p6eihSpAj69OmDvn37Qk9PT+KA86fg4GDo6uoCSExW7ezscOfOHXTq1ElV5s6dOyhfvnyqbejp6fHzIyLKYUqlEl5eXqhSpQpcXV1x5coVVKhQQeqwJCPZyKpCoUCrVq1Qt25dbNmyBatWrUqx3KJFi+Dp6Yl79+7h5MmTcHZ2xpkzZ3I4WiIigvcRYIkrsLUlsL934n+XuCLSaw/++OMPlC1bFseOHYOenh4GDx7MRCeHrFy5EgqFQvV4wYIF8Pf3x48//qg61rNnT2zbtg3v378HAFy8eBGXL19Gz549czxeIiJK2ZUrV1C1alV06dIFMTExaN68eb5OVAEJk1UhBPr164cXL16gQYMGKZZ5/fo1xo4di7Vr12LLli04fPgwunfvjr59+yKP77iTrWJjY3Ht2jW1hTOys15OEELg8ePH8Pb2hlKplDocorzH+wiwpwcQ7q92+PbTt3Cq3xmnj+7FuXPn0KdPH4kCzL9evXoFW1tb1KxZEw4ODpg3bx62bNmCGjVqqMqMHz8eHh4ecHZ2RtWqVdG4cWOMGzcuX933lB18fHzg4+OTY/VyQnBwMO7cuYNPnz5JHQpRvjJw4EA0atQIrVu3xoMHD7gd1P9INg1YW1sbrVq1SrPM4cOHoa+vj44dO6qO9e3bF2vXrsWtW7dQuXLl7A4z/dKYGqdp3rx5g+rVq8PX1xf29vbZXi8zYmNjcefOHVSrVi3VMt7e3vjxxx8REhICuVwOIyMjHDhwIM3pbUSUAUoFcHIMgP+/SOjlr0BJczkcLWSY30gfnWvIIHOvJF2M+di8efMwceJEPHz4ECYmJnB0dEw2qq2np4fDhw/j6dOn8Pf3h5OTE2xsbCSKOA25qC8FgMmTJ0NbWxtbtmzJkXqZ4ePjAz09PRQrVizVMuPGjcOiRYtgb28PPz8/DBo0CIsWLcqxGInym+joaNy8eRO1a9dGixYtMG7cON7u+BXJ71lNy5MnT2BnZ6e67wYASpcurTqXUrIqyVL73kcSv8h9OeJgYgM0nQu4tM62pw0NDcXr16/h4OAAPT09eHl5oXz58slWfPySQqHA3bt3ASTerxQQEAAzMzM4OzsDSHz/njx5AkNDQ7V9m75VL7skJcipjaQrlUq0b98e5cuXx+7duyGTydC1a1e0b98ejx8/hra2Rv+KE+UOr66o/r4FRCgx/kwsdj+Mx94OBmhRWgc/ldMBPvsnlitRS+Jg8ycTExO1kdTUODk5wcnJKQci+g4S9aVA4ncKfX192NnZ4dWrV4iPj//m3oWvXr1CUFAQtLS0cO3aNQCAq6urasGrN2/eICgoCA4ODmob3n+rXnaZPHkySpUqhSlTpqR4fvfu3ViyZAkuXboEDw8P3LlzB56ennBzc0OvXr2yNTai/EYIgf3792PUqFGws7PDuXPn0KJFC6nD0kga/U0+IiICpqamaseMjY2hpaWFiIiIFOvMnj07Z1eeTJoah6+SqfD3icc7bsuWTnb9+vUYMmQIihUrhuDgYHTo0AFr167FgwcP4Orqmmq96OhozJw5EwAwbdo06OnpoXr16li8eDF27tyJwYMHw8LCAmFhYShcuDD2798PJyenNOt97ebNm2r3T33NyMgI5cqVy+Q7kOjKlSt4/Pgx/vrrL8jlibPaJ0+eDGdnZ5w/fz7VKeZElAERHwAA/7xIQIe9UWjsoI1HAwughLk8xXJEGSZRXxoUFIQWLVrg8ePHKFKkCGQyGYoVKwYzMzPs27cvzboHDx7ErVu3AADDhg0DAGzcuBE2Njbo1KkTrly5AltbW/j5+WH06NGq7yap1Stbtqxa++/fv8erV6/SjMHJySnNFZ0zYsOGDWjevDk8PDwAABUrVkTr1q2xceNGJqtEWUgIgWbNmuHOnTuYPXs2evXqBZlMJnVYGkujk1UjI6NkI6MRERFQKBQwMjJKsc64ceMwYsQI1ePw8PDsG05PYWrc/xMAZMDJsYBziyydxuTr64tBgwZhy5Yt6NKlC6KiotCkSZN01S1QoAD27t0LR0dHHDx4UDWd18/PD71798by5cvRt29fxMXFoV27dujZsyeuXbuWar2UjBkzBlFRUameL126NLZt25aRl5yq27dvQ09PTy35dXJygomJCe7cucNklSiThBA4fuMl6sYJVCoix+HOhqhXIpWuo0ChnA2O8gaJ+lIgsb9SKpV4+/YtTExMcOTIEfzwww9o167dN+sOGzYMXl5eyabz9uvXDwEBAXj16hUsLS1x8eJF1K9fH9WqVUOzZs1Srfe1s2fPYtmyZWnGMG/ePNSuXTu9LzdNt2/fxujRo9WOVa1aFX/88UeWtE+U3338+BEPHz5EvXr1MGzYMNSoUUNt1gWlTKOTVScnJ+zcuRPx8fHQ0dEBADx//lx1LiU5utT+F1PjUiaA8HdZPjVu+/btKF26NLp06QIgcWP3P/74I90Ja2ptFi9eHH379gUA6OrqYtasWXBzc8PDhw/THK392n///ffdcQDqI7Nv374FANVUKUB9ZDY4OBgWFhbJ2rC0tOTiEESZ5O3tjeHDh+PWrVs43s0CVcxCUa9ESuvyyRKna9p9exoqUTIS9aUJCQnYvn07du3apfrC2Lp1a1SvXv2721QoFNi6dSs2b96s2g+xVq1aqhHKZs2apbutLl26qPr57/H1yGxQUBD09fXV+tOkkVkhBEJDQ5P1p5aWloiMjERsbCxX9yb6TvHx8Vi1ahWmTJmC1q1bo169elzcLgM0Ollt3bo1Ro0ahQMHDqj2htu0aROKFy+uGYsrpXfKWxZPjfPx8UmWrGf23lEfHx+UKVNG7ViZMmUgl8vh4+OToWQ1s9OAvxyZTbr/OGmqFKA+Mqujo5Pi6sTR0dGqCxxElHFbt27FgAED0L9/f+zevRvm7y/+b5qmDOojYP+butR0jkYvhEMaTKK+9M2bN4iPj0+xP/3e9S7evHmDuLi4ZP1p2bJlcfTo0Qy1ldlpwF+PzPr4+ODhw4d4+PCh6ljSyKxMJoOWllay/jQ6OhoAuP4D0XeKjY2Fu7s7lEoldu/enamBpfxK0r8+69atQ2BgIG7cuIH3799jxowZABKTFR0dHZQoUQJTp05F7969ce7cOQQHB+Pw4cM4dOiQ6v5ESaV3ylsWT40zNjbGhw/qnXZkZGSm23zz5o3asaioKCiVShgbG2eorcxOA/5yZPbFixdwdHRUuxL8JTs7O4SEhCA6Olq1sFRsbCw+ffoEOzu7DMVNlN8lJCRg165d+Omnn1C/fn3cunULLi4uiSfNWyfeN5jiAjhzsn0BHMrDJOxLgeT9Z2b606Q2v15XIyIiIsN9aWanAX89MtutW7c0F1iys7ODv7/6CLe/vz+KFSsGLS1eiCLKiGfPnuHdu3eoV68eli1bhlq1anEQ5TtJmqzGxsYiJiYG9evXBwDExMQAgNrKr+PHj0eDBg3w33//wdHREbNnz0bJkiUliTcZuxqJX9TC3yPle22yZ2pc5cqVsWPHDkRERKhWDzx16lS66ycldfHx8apjHh4e2LJlC0JDQ2FmZgYAOHnypNr9oCnVS0lmpwFnRN26dSGTyXDs2DG0b98eAHDixAkkJCSofq+I6NvOnj2LoUOHIjo6Gp6ennBwcEheyKV14n2DuWhrEcoFJOpLraysYG9vj9OnT6u2R4uNjcWFCxfg6emZrjYMDAzULs5aWlqiRIkSOHXqFGrVSpyyLITAqVOn0Lhx41TrpSSz04AzqkGDBjh27BhmzZoFmUwGIQSOHj3KtR+IMiAsLAwzZszAihUrMHr0aNSrV4/fRzNJ0mT1t99+S1e5qlWromrVqtkczXeQayUuqZ/DU+O6du2K2bNno02bNhgxYgT8/PxUo9LpUbhwYVhaWmLdunX48ccfYWFhgS5dumDhwoVo1aoVxowZg48fP+L333/HyJEjYW1tnWq97N665luKFSuGwYMHY9CgQYiJiYGWlhZGjhyJvn37okSJEpLGRpRbzJw5E3PnzsWkSZMwZMgQte3CkpFrcXsayloS9aUAMGXKFAwYMAAFChSAs7MzVqxYgdDQ0HTXd3V1xaJFi/DPP//A1NQUrq6umDNnDnr06AFjY2OUL18eW7Zswfv37zFmzJg062X31jXfMnbsWFSqVAk///wzOnXqhP379+PFixfYvXu3pHER5RafPn1C2bJl4eLiguvXr8PNzU3qkPIEDZhLm8u5/G9qnEkR9eMmNtm21L6Ojg7Onj0LBwcHzJ8/H3fv3sXhw4fTXV9LSwv79++Hr68vfv/9d6xduxba2to4e/YsPD09sXDhQuzbtw+zZ89WS4JTqpfd9PX1v3mhYtGiRfjjjz+wbds2bNq0CaNHj8bKlSuzPTai3CwiIgIrVqyAUqlEly5d8Pz5c4waNSrtRJUou0jQlwJAz549sWHDBpw4cQKrV69G06ZN8csvv6S7fr9+/dC1a1fMnz8fw4YNw6tXr9CxY0fs378f169fx8yZM2FiYoIbN26gUKFCadbLbqVKlUKxYsVSPV+iRAlcuXIFMpkMs2fPRnx8PK5cuaLa354oX1MqAN+LwIN9if9V/v/aLJcvX8aFCxdgaWmJQ4cO4cyZM0xUs5BMfDnnNg8KDw+HqakpwsLCki0PHRMTA19fX5QoUQL6+vqZeyKlQtKpcUFBQbC2tv7mPquUxZ87US6jVCqxa9cujBkzBiVKlMChQ4dgZWUldVi5Slr9Sl6WI/2pxH0pAAwePBgBAQHf3GeV2J9SPuF9JMX1Gt6UH4UxG8/gyJEjWLJkCfr06SNdjLlQevtSLu+WVTRoapyfnx8CAgJSPGdubp7qtj9ElPf9+uuvOH78OBYuXIhOnTpxI3LSLBrUlyYkJMDLyyvV8+XKlUt1z3ciyiO8j/zvFgX1sb0Xfm9RfmwftGtWF0+fPkXRokWliS8fYLKaR+jo6KBq1aowNDTE/v37sXfv3hTL1axZEwsWLMjh6Igop1z1v4o5N+ZgbJWxqG6TuF/k+/fvsXPnTowcORJjxozBokWL+CWbKBUlSpSAiYkJIiMj1bZN+9qWLVskX7eBiLKRUpE4ovq/RFUIgb3eCbAzlaFKUS1c71MArqU+AEUKSxtnHsdkNY8wNTVVbe8ycuRIjBw5UuKIiCinCSGw9PZSvAx7iaW3l6KCeQUsXboUM2fORPPmzRETE6M5q6kTaagv+8/Utk0jonzg1RXV1N877xUYejIGz4OV2NBKHzKZDK4FZUD4u8RyGjIjJC9iskpElEdc8b+CR58eAQAefXqEes3rITY0FseOHUt1L0YiIiJKQcQHAMCVNwlouC0KQ6vq4lgXQxjryVIsR9mDySoRUR4ghMDyO8sR9y4OYXfDUKhlIVh3t8bBHgehrc0/9UREROkVFxeHlX+dR6NgBaoV04L3oAKwN0tlE5UChVI+TlmCW9cQEeUBJx+dxL9L/sXzKc+REJ4AhUIBPy0/XP9wXerQiIiIco0TJ07Azc0NGw+dRZSuFeQyeSqJqgwwKZq4ajllGyarRES5nFKpxE+tfkL8x3iUmlYKRX4qAplcBrlMjuV3liOP71BGRESUJfbv34+uXbvit99+w927d1Gl37L/nfl65fz/PW46J8e318pvmKwSEeVSZ86cwbx583D1/VUUHlYYdiPsoFdET3VeKZR49OkRrvhfkTBKIiIizRUWFobRo0fDz88PrVu3xvPnzzFo0KDEW2hcWgMdtwEmRdQrmdgkHndpLU3Q+QiTVSKiXObly5do27Yt2rZtCx0dHSy/sxy6ZroplpVBxtFVIiKirygUCmzYsAGOjo64desWEhISoKOjA0tLS/WCLq2BYQ+Bnn8D7TYm/nfYAyaqOYSrbhAR5SJxcXGoXbs2mjdvjufPn8PM0gyH9x2GQMrJqIBAQGQA4pXx0NVKOaElIiLKb1avXo3Fixdj3bp1+OGHHyCTfT3V9wtyLW5PIxEmq/mQr68vKlasiPv376N48eLZXi+7vXnzBuvWrcO2bdugra0NHx8fqUMiylJKpRI7duzAx48fMXLkSDx8+BBmZmaq87tb7kZwTHCq9S30LZioEmWD/v37Q0tLC6tWrcqRetlJoVDg2LFjWL16Na5du4Z58+ahb9++UodFlKVev36NyZMnY+7cuejduzf69OkDfX19qcOiNHAacBa66n8VPxz6AVf9r0odSpoUCgXCwsKgVCpzpF5m+Pr6qn0pT0n79u2hra2Nzp07IywsLGcCI8oh169fR/Xq1TFx4kTY2dkBQLJ/E4WNCsPF0iXVn8JGhSWInOj75Ja+FAAiIyMRFRWVY/Uyo3///pgzZ06q57ds2YINGzZg8ODBUCqViI2NzcHoiLJXVFQUpkyZgjJlyqi+xxoYGDBRzQWYrGYRIQSW3l6Kl2EvsfT20my/P+zevXuoW7cubGxsUKdOHezduxdmZmZ4/PhxmvXCwsJQqVIlAICbmxvMzMzQvn17AMDTp0/Rpk0bFClSBA4ODhg1ahSio6O/WS87JSXIabl27RomT56MIkWKpFmOKLcJCQlBy5Yt0bx5czx58iRH/s0RSSmn+1KlUonJkyfD3t4ezs7O+O233zBo0CD07Nnzm3XHjRuHPXv2YNeuXTAzM4OZmRm8vLwghMCiRYvg6uqKQoUKoU6dOjh//vw362W3yMhIxMTEpHr+l19+wZEjR9CiRYu0p0MS5UIzZ87EyZMn8d9//2Hr1q0oWLCg1CFROnEacBa54n8Fjz49AgDV6pueRT2z5bkiIiLQuHFjtGrVChs3bsSrV6/Qo0cPhIWFQaFQpFnXxMQE58+fR6VKlXDp0iUUL14cOjo6iIyMRIMGDVCnTh1cuHABgYGB6NWrF4KCgrBly5ZU66WkWLFiiIiISDWGcuXK4eLFi5l6D77ETpXykpiYGCxevBgmJiYYNGgQXr16BUNDQ6nDIsoROdmXAsDChQuxevVqbN26Fc7Ozli8eDFWrVqFdu3afbPupEmT8PLlS7XpvMbGxli6dClmzJiBzZs3o0KFCti4cSMaN26Mhw8fwtHRMdV6X1uxYgUmTpyYZgw7d+5EixYtvuOVJ8e+lPKa27dvY+bMmdiyZQsmTpyI6dOnQy7nOF1uw2Q1CwghsPzOcshlciiFUrW3YQ2bGtnyx3/Lli2Qy+VYvXo1dHR04ODggFmzZqXrSrBMJlN1iiYmJqrphCtXrkRsbCw2btwIfX19ODo6YuXKlWjWrBmmT58OW1vbFOul5NGjR2leDdfS4n5URF8TQuDQoUMYOXIkzMzMsHz5cgBgokr5Rk73pQAwb948TJ48Gc2aNQMALFmyBMePH09XXQMDA+jo6EBbW1utT5w5cyYmTJiAH374AQAwbdo0nD59GvPmzcP69etTrfe1vn37olu3bmnGYGRklK5YifKTwMBATJgwATt37sTw4cMhl8thYGAgdVj0nZisZoEvrwQD6nsbZscV4QcPHqBy5cpqI5vVq1fPdJvu7u5qc/dr1qwJpVIJb29v2NraprstU1PTTMXy5chs0n0FX3boWT0yS6QJ/Pz8MHjwYEybNg29evXiRR3Kd3K6L33//j2CgoLU+k+5XI6qVat+9/2aAQEBCAoKgqenery1atXChQsXMtSWnp4e9PT0vl0wFV+PzEZFRUEul2PJkiWqY1k5MkukKcaPH4+QkBA8fPgQJUuWlDocyiQmq5n09ZXgJNl5RVihUCT7IqutnbmPMiEhIVkbWlpakMlkSEhIyFBbmZ0G/OXIrK+vLypVqgQ/Pz+1uIjygk+fPmHy5MmoXLkyevXqBT8/v1Sn1xOlh5eXFzZs2ICaNWsmG5X78OEDtm7dCn9/f5QpUwY9e/bUmMVFpOpLgeR9ira29ncnq0n95df9qba2dob70sxOA/56ZLZ///4oUaIExo4dqzrGkVnKK44fP44VK1bg4MGDqpmHlDcwWc2kr68EJ8nOK8KlS5fGf//9ByGEqvO+c+dOuusndaJfrurr5OSEf/75Ry0Rvn37NoQQcHJySrVeSjI7DfjLkdmkqcffWhGYKDeJj4/HmjVrMHnyZNSoUQNDhgwBAHaulClhYWH46aefEBoaioSEBLVExdfXF9WqVUPFihVRu3ZtrFixAhs3bsSFCxc0ImGVoi+1sbFBgQIFcPfuXVSoUEF1/M6dOyhdunS62tDW1lbrE4sWLYoCBQrg1q1bqFy5suq4l5eXqi9NqV5KMjsN+OuRWR0dHejr67M/pTzlyZMnGD58OG7cuIHp06dDS0sr0wM4pFl4l3EmJF0JliHlq70yyLD8zvIsX83w559/RmBgIKZNm4bY2Fi8fPnym1dfv1S4cGFoa2vjwYMHqmO//PILIiIiMGHCBMTExMDf3x/Dhw9Hy5YtUapUqVTrpcTU1FS1wmFKPyktJEGUn9y/fx+rV6/Gzp078ffff6f7izFRWvr164fOnTurJUVJJk6cCHt7exw7dgzjx4/H2bNn8fjxY2zYsEGCSNVJ1ZfK5XIMHDgQ06dPx5MnTxAfH4+FCxd+s4/7kq2traoukLguxLBhwzBr1izcv38fCoUCW7ZswenTpzFs2LBU66VET08vzb7UzMyMF7go3xs/fjwcHR3x/PlzDBw4kIlqHsRkNRPilfEIiAyAQModqIBAQGQA4pWpd0bfw9raGvv378fmzZthbGyMunXr4pdffkl3fX19fcyYMQPdunWDiYkJ2rdvD0tLSxw9ehQnT56EsbExSpQogaJFi2Lz5s1p1tMEnTt3hpmZGcaNG4fg4GBVJ/6tbXyIctKLFy/Qpk0bHDx4EO7u7nj48KFqUReizFq/fj18fX0xefLkZOeEEDh69Cg6deqkmtliZWWFJk2a4MiRIzkdajJS9aUAMGXKFNSoUQPlypWDlZUVzp49m66VgJP07dsXCQkJMDY2Vm1B88cff+CHH36Ap6cn9PX1MXnyZGzfvl1tpDWlelK7e/euqv8MCwvD6NGjYWZmlq7FG4lyikKhwLp169CkSRMoFArs27cPy5Ytg4WFhdShUTaRiezexExi4eHhMDU1RVhYGExMTNTOxcTEwNfXFyVKlPjuaVABkQEIjglO9byFvgUKGxX+rrbTIzY2Fnp6eggKCoK1tTUePHgAV1fXdNf//Pkz5HK52lSimJgY6OjopDldN6V62UGpVCI8PDzNaUuRkZEpXp02MTFJcYnyrPjcidLr8+fPmDlzJpYtW4YePXpg+vTpsLa2ljosyoS0+hUpeHt7o27durh8+TIcHR1Rs2ZNODs7q0ZNAwICUKRIERw4cAA//vijqt7vv/+O/fv34+XLlym2Gxsbq3bvZnh4OGxtbbOlP5W6L42Pj4eWlhbkcjkGDx6MgIAA7Nu3L0P1IyMjYWxsrOo7hRCIjo5Oc0XvlOpll6QFllL7fBQKBT5//pzsuK6ubqqvgf0p5aTz589j6NChiIiIwKJFi9CqVStuuZSLpbcv5Vh5JhU2KpytHei3ZGalQCDlvd3S0+Hk1FReuVz+zftruEAEaaKke8qvXr2Kq1ev4sqVK2r3xRFlhZiYGHTq1AmzZ8+Go6NjimWio6MBJP+7bWJigqioqFTbnj17NqZOnZp1waZB6r40s9NpdXR0kvVVMpnsm1tPpVQvu3wrFi0tLd7PShopaVxt1qxZ6NKlC4YOHZrp77+Ue3AacB6UNHUnpZ/OnTtLHR5Rnnf16lVUrVoVFy5cQOPGjXHu3DkmqpQt9u3bB19fX9y8eRMDBgzAgAED4OPjg8uXL2PAgAGIjo5WXbEODQ1VqxsSEpLm1exx48YhLCxM9fPmzZvsfCka58vbSlL6uXv3rtQhEuVpkZGRmDRpkmoE9eTJkxg9ejQT1XyGI6t5hKWlpeqLx9SpUzF+/PgUy3ExBqLs8+7dO4wdOxYHDx7EmDFj4OHhAQCcpkTZxsPDAwsWLFA7dvr0aVhaWqJChQrQ0tKCpaUlChUqBG9vb7Vy3t7eKFu2bKptZ3afz9xq/vz5EELAwMBAbdu0r3GxQKKsc9X/KubcmIOxVcaiWpFq+PPPPzFmzBgUK1YMy5YtA8C+NL9isppHyGQy1fQdAwMDGBgYSBsQUT6StOXTxYsXoVQq8fjxY9ja2kodFuUDTk5OyVb/3bFjB5ydnTFgwADVsc6dO2P79u0YPnw4jI2N4e3tjf/++w/bt2/P6ZA13pf9J6fFEmU/IQSW3l6Kl2EvsfjmYuxsvhNbt27FnDlz0KVLFyap+VyuSFYTEhLw7t076OrqonDhwvylJSKNIITA/v37MWbMGBw+fBidO3fmVHvSSFOmTMHFixdRoUIFVK5cGWfOnEGnTp3QsWNHqUMjonzuiv8V3PW5iw/7P+BN3BvcqHwD//zzj9RhkYbQ+GR10aJFmDZtGkxNTVWr6q1atQrNmzfPsufI4wsi01f4eVNWuHfvHoYOHYonT55g9uzZcHFxkTokIgCJq/x+PSJoZmaG69ev48yZM6p9tKtVq5blz82/r/kLP2/KrNjYWAycOBDPdz9HAdcCsOlsg+V3lqOGTQ0OThEADU9Wvby8MHLkSPz111/o2LEjlEolRo0ahY4dOyI4OBi6urqZaj/p/s2oqChOm81Hkla/5P279D2Stou6fv06qlSpgiNHjmjE9iVESX744YcUj2tra6NJkybZ8pxJ267ExcWxP81H2J/S9xJCIC4uDpfeXILvHV8UH1ocBcoUAAA8+vQIV/yvwLOop8RRkibQ6GT17du3AIBmzZoBSNzGpGnTpli8eDGCg4NRuHDmlrlPWqY9MDAQQOKy7ryKk3cJIRAVFYXAwECYmZll+552lLfEx8dj1apVmDVrFu7cuYN+/fpJHRKRxtDW1oahoSE+fvwIHR2dFPe4pryD/Sllhre3N4YPH45SpUrhU7NPKDGiBJRCqTovl8k5ukoqGp2sNmnSBJUrV8agQYMwcOBA1RLWffv2zXSimiSpnaSElfI+MzOzLPv9ofzh33//xZAhQyCTybBt2zbY2NhIHRKRRpHJZChSpAh8fX3x6tUrqcOhHML+lDIiLCwMkyZNwrp169C/f3806dsEo71GJyunFEqOrpKKRierBgYGmDlzJnr27ImzZ88iOjoadnZ2GDVqVKp1YmNjERsbq3ocHh6e5nMkdbAFCxZEfHx8lsVOmklHR4dXgCndIiIiUKBAAXh7e2PAgAEYOHAgp7sRpUJXVxeOjo6Ii4uTOhTKAexPKb0SEhIQHx8PIQQCAgJw69YtlClTBj8d+wkyyCCQ/N5nGWQcXSUAGp6sXrp0CS1atMBff/2Ftm3bQqlUYuLEifD09MTTp09hYWGRrM7s2bMxderUDD+XlpYW/+gSEYDEq78zZszA1q1b8fLlSwwdOlTqkIhyBblcDn19fanDICINcfbsWQwdOhSdOnXChAkT8NdffwEA4hRxCIgMSDFRBQABgYDIAMQr46Grlbk1aih3kwkNXsptxIgROHHiBB4/fqw6FhYWBjMzM+zZswcdOnRIVielkVVbW1uEhYVxERQi+qatW7di9OjRcHFxwdKlS+Hm5iZ1SKRhwsPDYWpqmu/6lfz6uoko4969e4chQ4bg33//xaRJkzBkyJBkC6MGRAYgOCY41TYs9C1Q2IjTzPOq9PYpGj2yam5ujtDQUCgUCtWo58ePH1XnUqKnpwc9Pb0ci5GI8oZPnz7B0tISgYGBWLVqFdq2bcupR0RERBkQEREBbW1tyGQyWFtb49mzZ6ne11zYqDCTUfomjV6u76effkJERAR69OiBy5cv49SpU+jWrRvKlCkDT0/ecE1EmffmzRt06dIFHh4eSEhIwO+//4527doxUSUiIkonpVKJHTt2wMnJCbt27YKNjQ3WrFnDBbgo0zQ6WS1VqhSuX78OfX19/P7775g1axZq166N8+fPcx83IsoUIQSmT58OZ2dnaGtr4+LFi9DW1ujJJkRERBrn0aNH8PT0xLhx47BgwQL8/PPPUodEeYjGfzNzcXHBxo0bpQ6DiPIIIQT8/f1RtGhR6Orq4vTp06hevbrUYREREeUqAQEBsLKygq6uLho3bozRo0fDyMhI6rAoj9HokVUioqx0584d1KlTB+3bt4cQAmPGjGGiSkRElAExMTGYM2cOHB0dcfr0aTg6OmLq1KlMVClbMFklojwvISEB/fr1g6enJzw9PXHq1Cnek0pERJRBFy5cQNmyZfHXX3/h2LFjaNq0qdQhUR6n8dOAiYi+V1xcHPz9/WFvb49SpUrhwYMHcHBwkDosIiKiXOX58+coVaoUjI2NMXbsWPzyyy+qnTqIshNHVokoTzpx4gTc3Nzw+++/AwBGjx7NRJWIiCgDgoOD8dtvv8HNzQ337t1DxYoV0bdvXyaqlGOYrBJRnvL582e0aNECXbt2xeDBg/Hnn39KHRIREVGus3fvXjg6OsLHxwd37txBhQoVpA6J8iFOAyaiPCE0NBTBwcEoUaIEGjZsiG3btsHS0lLqsIiIiHKVu3fvokKFCihatCi2b9+O5s2bSx0S5WMcWSWiXE2hUGD9+vUoXbo0Vq5cCZlMhuHDhzNRJSIiygAfHx/8+OOPqFOnDt6+fYsaNWowUSXJMVklolzL398fHh4emD17NtatW4cFCxZIHRIREVGus3TpUpQrVw7W1tZ4/vw5ihUrJnVIRAA4DZiIcqHXr19DqVTC1tYW/fv3R8+ePaGvry91WERERLmGUqnElStXULNmTZQvXx6XL19GxYoVpQ6LSA1HVoko14iKisKUKVNQpkwZHD58GFpaWujfvz8TVSIiogy4fv06qlevjm7duiE8PBx169ZlokoaickqEeUK3t7ecHZ2xsmTJ/Hff/9h6NChUodERESU64wcORL169dH8+bN4e3tDRMTE6lDIkoVpwETkUa7ffs2ChcuDAcHB8yfPx8dOnSAXM7rbEREROkVExODS5cuoWHDhmjRogWGDh2K4sWLSx0W0TfxGx8RaaTAwED07dsXNWvWxJUrV6Cnp4dOnToxUSUiIkonIQQOHjwIFxcXTJw4EQkJCahfvz4TVco1+K2PiDTOhQsX4OjoiODgYDx8+BDt27eXOiQiIqJcp127dhgwYAAmTJiAy5cvQ1ubkyopd2GySkQa4/jx4wgODkaFChVw8OBB7N+/HyVLlpQ6LCIiolzj06dPOHHiBABgyJAhePbsGXr37g0tLS2JIyPKOCarRCS5J0+eoFmzZujevTsePXoEExMT1K9fX+qwiIiIco34+HgsX74cjo6O2Lx5M4QQqFu3LkxNTaUOjei7cS4AEUlq79696NGjB/r27YudO3fCwsJC6pCIcp3AwECcPXsWYWFhKFOmDGrVqpWsTExMDI4fPw5/f3+UKVMGDRo0kCBSIsoOCoUC1apVQ2RkJHbu3IlmzZpJHRJRluDIKhHlOIVCgc2bNyMyMhJ169aFl5cXli1bxkSV6DssW7YMNWvWxNGjR3Hz5k106tQJtWvXRmRkpKpMUFAQ3N3dMWnSJNy+fRvdunVD27ZtoVQqJYyciDLLx8cHx48fh5aWFpYtW4YHDx4wUaU8hSOrRJSjzp8/j6FDhyIiIgJVq1aFi4sLrK2tpQ6LKNeqVKkSvL29VQunfPjwAba2tti7dy969eoFAJg0aRKUSiVu3boFQ0NDvHjxAmXLlsWuXbvQrVs3CaMnou/x+fNnzJw5E8uWLcNvv/2G5s2bw9PTU+qwiLIcR1aJKMcsWbIErVq1QpcuXfDo0SO4uLhIHRJRrlezZk21FT6NjIygpaWltpjK3r170a1bNxgaGgIASpUqhYYNG2LPnj05Hi8RZU5ERATKlCmDq1ev4sqVK5g7d67UIRFlG46sElG2ioyMxLp16zB48GB06NABnTp1QpEiRaQOiyhPeffuHf788098/vwZx48fR5cuXdC5c2cAiVOAg4KC4OTkpFbH2dkZR44cSbXN2NhYxMbGqh6Hh4dnT/BElC5Xr15FeHg4mjRpgsOHD6NSpUqQyWRSh0WUrTiySkTZQgiBP//8E87OztizZw8+fvyIokWLMlElygbx8fEICAjAmzdvEBAQAIVCgYSEBACJ0wUBwMzMTK2OmZmZ6lxKZs+eDVNTU9WPra1ttsVPRKl79+4dunfvjkaNGuHFixcAAHd3dyaqlC9818hqQkIC/v33X/j6+iIuLk7tnLu7e4qrEBJR/jJq1Cjs3r0b8+bNQ5cuXdipEqXAx8cHV65cQXBwMIQQaucGDRoEHR2ddLVjb2+PBQsWAEgcSXV1dcWCBQvwxx9/qKb+fj0yGhYWBiMjo1TbHDduHEaMGKF6HB4ezoSVKIf5+/vD2dkZrVu3xuPHj/lvkPKdDCerMTEx8PT0xKNHj2BnZ5esI1UoFExWifKpgIAArF+/HhMmTMDw4cMxdepUFChQQOqwiDTSrl270LNnT1hZWcHS0jLZ+f79+6c7Wf2SlZUVKleujNu3bwMAChUqBFNTU/j4+KiV8/HxgaOjY6rt6OnpQU9PL8PPT0SZI4TA/v37YWFhgfr168PLyyvZNH6i/CLDyerRo0cRERGBd+/epdi5ElH+Exsbi2XLlmH69Olo0qQJIiIiUKxYManDItJoM2fOxOLFizF48OBMtfPo0SOULVtW9Tg0NBS3b99Gjx49VMd++OEH/Pnnnxg+fDi0tbXx/v17nDp1CosWLcrUcxNR1rp//z6GDh2Kx48fY/Xq1QDARJXytQwnq6GhoWjQoAETVSJS6d69O54+fYojR46gbt26UodDlCuEhoaiY8eOmW7n559/RvHixVGuXDl8/vwZ+/btg42NDX7//XdVmZkzZ6JatWqoX78+atWqhb1798Ld3R0///xzpp+fiLLG/fv3Ua1aNfz22284fPgwTExMpA6JSHIZXmCpUqVKuHPnTnbEQkS5iLe3N/744w8IIbBkyRLcunWLiSpRBlSqVEk1VTczrl+/jl9++QUymQyFChXC2rVrcfPmTbWLysWKFcP9+/fRqVMnAIn7rp45cwa6urqZfn4i+n7x8fFYunQprl69inLlyuHp06eYO3cuE1Wi/8nwyKqlpSUUCgU6duyIDh06wNjYWO18iRIlOF2BKA8LCQnBlClTsG7dOvTv3x8JCQmwsbGROiyiXKdp06bo1asXxowZAycnJ8jl6tePGzVqpLZXampkMhmaN2+O5s2bp1nOwsICgwYNylTMRJR1/vnnHwwbNgwymQyVK1eGTCbjAkpEX8lwsnrgwAHcvn0bt2/fxoEDB5KdHzlyZJZvTnzq1CmsXLkSvr6+qFixImbMmMF/zEQSadGiBYyNjXHr1i24uLhIHQ5RrjVhwgRERESoTdf90rdW6yWi3OvUqVPo3Lkzpk6dil9//fW7FlMjyg9k4uu18jXMhg0b8Ntvv2HWrFlo1KgRnjx5gn379mH37t3pqh8eHg5TU1OEhYVxSgVRCt6FRiMkMi7V8+ZGunh25xr+++8/TJ8+HR8/foSVlRW3oqF8K7/2K/n1dRNllbCwMMyYMQM9e/aEi4sLQkJCuAYM5Vvp7VO+a5/VnBIUFIQhQ4ZgxowZGD58OADA1dUVP/74o8SREeUN70KjUX/BOcQmKFM8Hx8agPBzmwD/B5j0v/tTra2tczhKIiKi3EuhUGDLli0YP348XFxc0Lt3b8jlciaqROnwXclqQkIC/v33X/j6+iIuTn1Ext3dPcv2WT1y5AhiYmLwyy+/qB1Pzz08RPRtIZFxqSaqQqlA4J4/oFfMFX+f90LdiqVzODqivM/HxwdXrlxBcHAwvp7oNGjQIE4NJMoDduzYgenTp2PVqlVo27YtZyYRZUCGk9WYmBh4enri0aNHsLOzS9aRKhSKLEtWnz59Cnt7e1y4cAELFixAVFQUKlasiD/++AN2dnYp1omNjUVsbKzqcXh4eJbEQpQfCKFEpPd5JIT4w6xmVxTpuRRyPUNYFSwkdWhEec6uXbvQs2dPWFlZpTjC0r9/fyarRLnUmzdvMGHCBEybNg1dunRBx44dYWBgIHVYRLlOhpPVo0ePIiIiAu/evcv26QtxcXF4//495s6dizlz5kBfXx/Tpk2Dp6cnHjx4AHNz82R1Zs+ejalTp2ZrXER5Uaz/UwSfWQfF508wr5u496Jcz1DiqIjyrpkzZ2Lx4sUYPHiw1KEQURaJiorCggULMHfuXLRt2xb6+vrQ0dHhhSei75ThZDU0NBQNGjTIkXn21tbWiImJwbp16+Dq6goA2LlzJywsLHDixAl06dIlWZ1x48ZhxIgRqsfh4eFcOZjoG5RxMfh4cBYKuDWGSdV2kOvqSx0SUZ4XGhqKjh07Sh0GEWWh5cuX4++//8bp06dRvXp1qcMhyvXk3y6irlKlSrhz5052xJJM1apVAQBmZmaqYwUKFICOjg4iIyNTrKOnpwcTExO1HyJKLiYmBhtWLEbY1T2Q6+qjaP/1MKvVlYkqUQ6pVKkSbt++LXUYRJRJd+7cQatWrfDx40cMHz4c165dY6JKlEUyPLJqaWkJhUKBjh07okOHDjA2NlY7X6JECTg5OWVJcHXr1kXZsmUxZ84cLF26FFpaWpg/fz7kcjkaNGiQJc9BlN8IIXD48GGMHDkSugZG0CvfFQAg09aVODKi/KVp06bo1asXxowZAycnJ8jl6tePGzVqxAUFiTRYYGAgJk6ciB07dmDo0KHQ19eHri77UqKslOFk9cCBA7h9+zZu376NAwcOJDs/cuRIzJ07N0uC09LSwuHDh9G1a1eYm5tDW1sbZmZm2Lt3L0qWLJklz0GU33z48AHDhg3DhAkTUKVJO/yw6qrUIRHlSxMmTEBERAR+//33FM+HhYXByMgoh6MiovSaO3cuPn78iAcPHsDBwUHqcIjyJJn4eq18DfXx48fv2pOKm5gTAcHBwZg8eTIcHBwwbNgwJCQkQFtb+5v7rAKAnrYc/42qi6JmXMWQCMi//Up+fd1EXzpx4gTmzZuHI0eOwMDAANra37ULJFG+l94+JcP/wl69egVbW9tk05WAxNV77927Bw8Pj4w2+03W1tZZ3iZRXpeQkIC1a9di0qRJqFq1KgYNGgQAqs61qJkB/htVFyGRcam2YW6ky0SVKBu8fPky1VlCN2/eRPny5TmlkCib3XoVgjfBUamet7UwhLudOZ4+fYoRI0bg6tWrmDZtGhNVohyS4X9le/fuxbFjx7Bt2za1VXa9vb3RtWtXNGrUKFuSVSLKGCEEfHx8sGrVKmzfvh3NmzdPsVxRMwMmo0QSqF69OqZNm4b+/furjikUCsyZMwdTp05FSEgIk1WibHTrVQjarb6SZhkhBA4M9MSKOXNQokQJbNu2LUd2xCCiRBleDbhbt27Q1dVF+fLlsWfPHgDAypUrUblyZZQpUwYTJkzI8iCJKP1evnyJtm3bYsuWLXBycsLDhw9TTVSJSDrz58/H77//jjZt2iAoKAh+fn6oW7culi5digMHDvB+VaJsltaIqlAq8PnuSXzYORovP4Ri06ZNWLFiBRNVohyW4WS1cOHCOHnyJP744w/06NEDTk5OGD9+PNauXYtdu3bB1NQ0O+Ikom/4/Pkzxo8fD1dXV1hZWaFFixYAAJlMJnFkRJSSHj164M6dOwgICEC5cuVQvnx5GBsb48GDB2jZsqXU4RHlWzFvHiJg2wiEX9sLkyptoa2tw76USCLfNdleJpPB2dkZxsbG8PX1haurKzw9PbM6NiJKB6VSCZlMBm9vb1y+fBmXL19GxYoVpQ6LiNLB2toaJUuWxO3bt5GQkIBGjRqhYMGCUodFlC8JoYRMJsfn28dg6FwTJpV/gExbl4kqkYQyPLIaExODoUOH4ocffsCQIUPw9u1b2NnZoUKFCtiyZUs2hEhEqbl+/Tpq1KiBo0ePomrVqjh37hwTVaJc4sqVKyhfvjy8vb1x9+5d7NmzB9OnT0eTJk3w/v17qcMjyjeU8TEIvbQLH3aNhRACVq1Hw7RaB+4/TqQBMpysrl69GseOHcPFixfxxx9/oGDBgjh48CAWLlyIwYMHY8GCBdkRJxF9wd/fHz179kT9+vXRrFkzNGzYEACn/BLlJs2bN0ebNm1w/fp1uLi4oH379rh//z4UCgXKlSuH2NhYqUMkytOEEIh8fAH+639F9MtbMK/XGzKZjH0pkQbJ8DTgOnXqoE+fPjA2NlY73rdvX9StWxc3btzIsuCISF18fDx0dHTw4MEDxMXF4fHjxyhevLjUYRHRdzh8+DDq1KmjdqxYsWI4ffo0FixYAKUy9f2PiShz4uPjAQBRTy7BrHZ3GJWtC5ksw2M4RJTNMvyvslKlSskS1SSOjo7o2rVrpoMiInVCCBw8eBDOzs64evUqmjRpgj///JOJKlEu9nWimkQmk+H333+HgQG3lCLKaoGBgejbty9++OEHyGQyWP84HgVc6zNRJdJQ/JdJpOEePnyIhg0bYsCAARg/fjyqVKkidUhERES5Snx8PBYuXAhHR0cEBwdjxYoVsLUwTFfd9JYjoqz3XasBE1H2i46OhoGBAby9vVGxYkUcOHCAW0MRERFlUHR0NLS1tXHt2jUcOHAADRo0UJ3b/2uNNPdbtbUwhLudeU6ESUQpkAkhhNRBZKfw8HCYmpoiLCwMJiYmUodD9E0JCQlYs2YNpk6diqtXr6JUqVJSh0REX8iv/Up+fd2Uez158gQjRoyAtbU1tm7dKnU4RPSF9PYpnAZMpEHOnj2LChUqYMWKFdi2bRsTVSIiogyKiIjAiBEjULFiRZQqVQqLFy+WOiQi+k6cBkykAcLCwmBqago/Pz/06dMHgwYNgo6OjtRhERER5RoKhQLR0dGQy+UICgqCl5cXypYtK3VYRJQJHFklktDnz58xduxY2Nvb49OnT/j5558xbNgwJqpElCHh4eE4f/48zp8/j5CQkFTL3b59G3///TdevHiRg9ERZb/z58/D3d0d06dPh6GhIbZt28ZElSgPYLJKJJFdu3ahdOnSuHr1Ks6ePQtLS0upQyKiXOj333+Hs7MzJk+ejPHjx6N48eJYs2aNWpmoqCg0bNgQTZs2xaJFi1ChQgUMHz5cooiJss6HDx/QsWNHtGrVCl26dMG0adOkDomIshCnARPlsMDAQBQsWBDh4eFYtmwZ2rdvD5lMJnVYRJRL2dnZ4cWLFzA0TNxeY8uWLejduzfq1asHJycnAMCMGTPw7NkzPHr0CNbW1rhx4wZq1KiB+vXro1WrVlKGT/RdIiMjIYSArq4urK2t8fTpUxQpUkTqsIgoi3FklSiHvHv3Dt27d0f58uURGxuLAQMGoEOHDkxUiShTBg8erEpUAaBjx45QKpW4d++e6tj27dvRo0cPWFtbAwCqVKmC2rVrY/v27TkeL1FmCCHw559/wtnZGRs3boS5uTlWrlzJRJUoj2KySpQD5s2bBycnJyiVSty4cQN6enpSh0REedR///0HAHBxcQEAhISE4O3bt3Bzc1Mr5+bmhvv376faTmxsLMLDw9V+iKT07Nkz1KpVC6NGjcKcOXMwZMgQqUMiomzGZJUomwgh8OrVKwCAoaEh/vnnH+zcuRO2trYSR0ZEeVVAQAB+/fVXdO3aFa6urgASVxsHAHNzc7WylpaWCA0NTbWt2bNnw9TUVPXDv10klQ8fPiAmJgYFChRAgwYN8PTpU3Tt2pUzk4jyASarRNng/v37qF+/Plq2bAkhBAYPHgxPT0+pwyKiPCwoKAiNGzdGqVKlsH79etXxpJkcUVFRauUjIiKgr6+fanvjxo1DWFiY6ufNmzfZEzhRKmJjYzF//nw4Ojri77//ho2NDaZOnYoCBQpIHRoR5RAmq0RZSKlUYtCgQahWrRqqVKmCy5cv88ovEWW7T58+oWHDhjA3N8fff/8NAwMD1blChQrB0NAQr1+/Vqvz+vVrlCxZMtU29fT0YGJiovZDlFOuXbsGV1dX7NixA0eOHEH79u2lDomIJMBklSgLxMfH49mzZ5DL5ShTpgzu3buHuXPn8ssdEWW74OBgNGzYEKampjh+/DiMjIzUzsvlcjRp0gT79u1THfv8+TNOnjyJZs2a5XS4RGl6+vQpEhISYGFhgZEjR+LWrVuoW7eu1GERkURkQgghdRDZKTw8HKampggLC2PiQNnin3/+wbBhw2Bvb48TJ05IHQ4RZTNN6leUSiU8PDzg6+uLxYsXqyWqFSpUQKlSpQAA3t7eqF69Olq0aIF69ephy5YtCA4Oxs2bN9M9pVKTXjflPSEhIZgyZQrWrVuH06dP89YZojwuvX0K91kl+k4xMTHo0KEDLl26hKlTp+LXX3+VOiQiymcUCgVKlCiBEiVK4OjRo2rnDAwMVMmqi4sLbt26hbVr1+LcuXNo1qwZBg8ezHv/SCMcPnwYvXv3hru7O27duqVayZqIiCOrRBkUFhaG9+/fw9nZGWvWrEH79u1hZWUldVhElEPya7+SX183ZZ8bN27Aw8MD9+7dw9u3b9GiRQuu80CUT6S3T+E9q0TppFQqsXHjRpQuXRpLliwBAAwYMICJKhERUQb4+vqiXbt2aNiwIZ4/f44KFSqgZcuWTFSJKBkmq0Tp8OnTJ1SpUgXTp0/HqlWrsHr1aqlDIiIiynXWrFmDsmXLwszMDM+ePUPp0qWlDomINFiuumd17ty5eP/+PaZOnQpTU1Opw6F84M2bN4iKikLp0qUxaNAgdO7cWW1LCCIiIkqbUqnE2bNn0aBBA1SqVAkXLlxA5cqVpQ6LiHKBXDOyumbNGixcuBBLly7F58+fpQ6H8rioqChMmzYNzs7O2L9/P2QyGX7++WcmqkRERBlw48YNeHp6omfPnvj48SOqVKnCRJWI0i1XJKsPHz7EjBkzMH/+fKlDoXzAx8cHZcqUwdGjR/Hvv/9i/PjxUodERESU60yYMAF169ZFo0aN8PTpU1hbW0sdEhHlMho/DTg6OhqdO3fGkiVLYGZmJnU4lIfdvXsXxsbGsLe3x8KFC9G2bVvI5bnieg4REZFGiI2Nxb///ouWLVuiZcuW6Nu3L+zt7aUOi4hyKY3/Jj5kyBBUqVIF7du3T1f52NhYhIeHq/0QpeXjx4/o168fatSogStXrkBLSwvt27dnokpERJROQggcPnwYLi4umDRpEqKjo1G9enUmqkSUKRr9bXzPnj04e/Ysli1blu46s2fPhqmpqerH1tY2GyOk3O7mzZtwdHTEx48f8eDBA3Tv3l3qkIiIiHKdbt26oU+fPhg7dixu3rzJNR6IKEvIhBBC6iBSU7p0adja2qJcuXIAEldmPXDgAHr37o0WLVrgxx9/TFYnNjYWsbGxqsfh4eGwtbXlJuak5sSJE3Bzc4OlpSWuXr2KevXqSR0SEeUS6d3IPK/Jr6+bUhccHIxz586hbdu2uHbtGpydnXnLFhGlS3r7FI0eWR07dixatWoFe3t72Nvbo1ChQgAAW1tbWFhYpFhHT08PJiYmaj9ESZ4+fYoWLVqga9eu8Pb2hr6+PhNVIiKiDEhISMDKlSvh6OiIrVu3QqlUolq1akxUiSjLafQCS7/88ova49OnT2P16tXo3bs3ihUrJlFUlFv9/fff6NChA3r37o1t27bB0tJS6pCIiIhyFSEE6tSpg6CgIGzfvh3NmzeXOiQiysM0OlklyiyFQoGtW7eiTZs2qF27Nm7cuKGaVk5ERETp8/LlS9y+fRvt27fH0qVL4ebmBl1dXanDIqI8TqOnAX/NyckJixcv5jQTSpeLFy/Cw8MDM2bMwKtXr2BiYsJElYiIKAM+f/6M8ePHw9XVFZcuXQIAVK5cmYkqEeWIXJWs2traYtiwYShQoIDUoZCGW7duHZo3b46OHTvC29sbFStWlDokIiKiXCU2NhZubm64dOkSLl++jCVLlkgdEhHlM5wGTHlGVFQUVqxYgcGDB6Nt27Zo2bIlbGxspA6LiIgoV7l+/Trevn2Ldu3a4ciRI3B1dYVMJpM6LCLKh3LVyCpRSoQQ2L17N5ydnXHgwAEEBgbCysqKiSoREVEG+Pv7o0ePHqhfvz5evXoFAChXrhwTVSKSDEdWKdebPHkyNmzYgLlz56Jr166Qy3kNhoiIKCOCg4Ph4uKCZs2a4fHjxyhevLjUIRERcWSVcqfAwEBMnDgR8fHxGDRoEJ49e4bu3bszUSUiIkonIQQOHjyII0eOwMLCArdu3cKff/7JRJWINAa/2VOuEhcXh4ULF8LR0RHe3t4IDw9HoUKFuOgWERFRBjx48AANGzbEgAEDEB0dDQBwcHCQOCoiInWcBky5yq+//oobN27g4MGDqF+/vtThEBER5TovX75E1apVMWjQIBw4cACmpqZSh0RElCImq6Txnjx5gvXr12P+/PmYPXs2LCwsoK3NX10iIqL0SkhIwJo1a+Dk5IRGjRrBx8cHRYoUkTosIqI08Rs/aazQ0FBMmzYNq1evRt++fREbG4uCBQtKHRYREZHk5hx/jJt+wYhXihTPO1oXwIgmTihqZoDTp09j2LBhSEhIwNq1awGAiSoR5QpMVkkaSgXw6goQ8QEoUAiwqwHItdSKdO7cGQqFAl5eXihbtqxEgRIREWmWOccfY82FlwAAbSSgh9YpFJcF4rUoiG2KxkiANu6/DcOhu/5YXM8QPTq2x5QpUzBo0CDo6OhIHD0RUfoxWaWc530EODkGCPf//2MmNkDTuTj/0RSHDh3CokWL8Oeff8LMzIz7uxEREX3h/rswAMBYrV3oq30MWrL/H12doL0TyyMbY/IlAUPHqijh0ge+vr4wNzeXKlwiou/G1YApZ3kfAfb0UE9UAbx6/Q4dO7RHy+ZNUbhwYQghYG5uzkSViIgoBWO1dqG/9t+Q4/8TVaUQ2HY3FnPX7kMh/wuQ6RpCJpMxUSWiXIvJKuUcpSJxRBXq99cIIdDmr0gU0AWe/V4MY34fxf1SiYgyKCQkBF5eXggLC0u1TFBQEB4+fIiIiIgcjIyymlwZh77axwAAX17TPfwkAVPPx2JZU3086hULQ+uiEkVIRJQ1mBFQznl1RTWiKoTAnw/iMfKfGMhkMlz62QibfjBAEQQmliMionTx9vbGzz//DGdnZ3h4eODixYvJyigUCvTu3RtFixZF69atUbBgQSxdulSCaCkrNI78G1oyAZkMeBeuRI+D0XgUqMAPztp4PKgAOpTVgY5coIfWKalDJSLKFCarlHMiPgAAbvkrUGtzFEb9G4NKRRIXVTLSlSUrR0RE33br1i3Url0bd+/eTbXMwoULcfjwYTx69AgvX77Enj17MHz4cFy4cCHnAqUsU1DxHtHxAjMvxMJpRQQUQsBMXwa5TAYDnf/vT4vLAiWMkogo87jAEuWcAoUQrxDovD8anctq42RNQxTQTeGe1AKFcj42IqJcqnv37gASt/tKzfr169GjRw+UKlUKANCyZUt4eHhgw4YNqF27dk6ESVkoUKsIttyNx8En8finmyE8i6f8de61KIjKORwbEVFWYrJK2S42NhbLli3Dp6CPmFOoKB4N9IeuVkolZYmrAtvVyOkQiYjyrPDwcLx48QIeHh5qx6tWrYr//vtPoqjoe9y/fx9jx46Fbs1e+KuSLvq560BLnvyirxCAAnJsUzRGWwniJCLKKpwGTNlGCIGjR4/C1dUVO3bsQJOmzYCmc6GrJQPwdef6v8dN5yTbb5WIiL5fcHAwAMDS0lLtuJWVFT59+pRqvdjYWISHh6v9kDSCgoLw66+/omrVqnB1dUWF0iWwWdkScpkMQn3NQtXjDQnNIWQ6MDfSzfmAiYiyCEdWKduEhYVhxIgRGDlyJPr06QNt7f/9unXclso+q3MAl9bSBEtElEfp6OgASEw+vxQdHa06l5LZs2dj6tSp2Robpc/KlSvx7t073L9/H46OjgCAOYYTceiBHlpHH4QWlKqySpkch/Tb4HmJobjQxAlFzQykCpuIKNOYrFKWCgkJwdSpU2Fubo7JkyfjyZMn0NL6aqTUpTXg3CJx1d+ID4n3qNrV4IgqEVE2KFSoEHR1deHvr76/tb+/P4oXL55qvXHjxmHEiBGqx+Hh4bC1tc22OEndqVOnMHXqVBw6dAgTJ05M1peObV4GaL4JSFgD3FwPhPgB5vbQ8uiLdtq6aCdN2EREWYrJKqXp5MMAPP/wOdXzjoWM0dS1MBQKBdavX4+JEyfC3d0dixcvBoDkiWoSuRZQolZ2hExERF/Q1tZGnTp1cOzYMfTv3x8AEB8fj3/++Uf1OCV6enrQ09PLqTDpf54/f46RI0fi4sWLmDp1KszMzFLvSwFAWxeoPijnAiQiykFMVilVJx8GYMCOW98st7prJbiaJWD16tXYsmULWrRoAZkshVV+iYgoy4WEhMDHxwcREREAgBcvXsDLywuFCxdGsWLFAABTp05FnTp1MG7cONSvXx9r1qyBXC7Hb7/9JmXo9AUhBGQyGVauXIlixYrh+fPnsLKykjosIiJJyYT4+tb8vCU8PBympqYICwuDiYmJ1OHkKsvPPMfCf5+lej4+NAChZzehSf06OLB6lqqjJSLKyzStX/nvv/8wevToZMe7d++OoUOHqh5fvnwZS5Ysgb+/P8qUKYMJEyagRIkS6X4eTXvdeYVSqcTmzZuxYsUKXLp0CYaGhuxLiSjPS2+fwpFVyjBlXDTCru3D55uHYFimNtzrtwQAdq5ERBKoX78+vLy8vlnO09MTnp6eORARpdfly5cxdOhQfPr0CQsWLGCiSkT0FSarlG5CJK42mBD+EbFvH6FQlznQK+IIEwtriSMjIiLKPRQKBbS0tLB161b8+OOPGDFiBAwMuGovEdHXmKxSusT6P0XwmXUwrtgCBVzro3CXOVKHRERElKtERUVhwYIFOHDgAG7duoV169ZJHRIRkUZjskppSogIRuj5rYh6cgkmVX6EYekaUodERESkeSKCga3NgYgAoEBhoOdxoIAFgMTFk/bu3Yvff/8dBQsWxKpVq9Je4ZeIiAAwWaVUxMXFQQgBRVggRHwsbPqsgrZpIanDIiIi0jzzHYHIwP9/HB0CLCgBGBVE3NBH0NHRwYkTJzB16lT06NEDcrlculiJiHKRXPHX8tOnT7h58yYCAgKkDiXPE0Lg8OHDKFOmDJ54XYJeUWdYtxnLRJWIiCglXyeq//MxUol+u1+hvpMpZDIZNm/ejF69ejFRJSLKAI3+i/no0SM0bdoUzs7OGDhwIEqXLo1WrVohLCxM6tDypMePH6Nx48bo06cPxo4dixbNmqSrnmMh42yOjIiISANFBAORgdhoYoxKdsWw0cQYCUqBxVdj4bg8Ah+jBLa20k4sR0REGabR04B9fX0xfPhwNGmSmDR9/PgRNWrUwMiRI7FhwwaJo8s7IiMjYWRkhNevX6NcuXLYu3cvzMzMAABaWlp4/uFzqnUdCxmjqWvhHIqUiIhIg2xtDiWANeamiJfLsdLQGD3DwnE/UIl9HQ3RsKS2qhwGXZM0VCKi3Eijk9WWLVuqPba2tkabNm1w7NgxiSLKWxISErB27VpMmjQJZ86cQZMmTVQXBpI0dS3MZJSIiPKtOccf4/67lGd0rQl6i92mJgj7EI/3f76HTFuGTd1MsfmHr/ZKjeBtTERE30Ojk9WU3LhxA46Ojqmej42NRWxsrOpxeHh4ToSV61y4cAEDBw5EfHw8tm/fjgoVKkgdEhERkUaZc/wx1lx4mer5ZyiASVc+IvDsO5jXNkfBNtZYb6SFPmHh6vdZFeBFXyKi76HR96x+beXKlbhy5QrGjx+fapnZs2fD1NRU9WNra5uDEWq+4ODE+2Y+fvyIX375BQ8ePEDz5s0ljoqIiEjzpDaiKpQKKGMi0NOsPuJiBRwmOcCmuw20jXUQI5djg6mJeoWex3MgWiKivCfXJKt79uzB8OHDsWHDBlStWjXVcuPGjUNYWJjq582bNzkYpeb6/Pkzxo8fDzs7O7x9+xbt2rXDiBEjoKurK3VoREREuUbMm4cI2DYCIRe2QlboGor2soG+rf7/FxAC681MoEx6bFRQtd8qERFlTK5IVvft24fu3btjzZo16NmzZ5pl9fT0YGJiovaT31z1v4ofDv2Aq/5XAQB79+6Fk5MTLl68iPPnz6NYsWISR0hERJS7KKLC8PHwXATumwpDJ08UbO8ImTwekH11f6pM9v+jq0YFgd+fSxMwEVEeoPHJ6oEDB9C1a1esWrUKv/zyi9ThaDwhBJbeXoqXYS8x5985EEIgLi4OixcvxoULF1CpUiWpQyQiItJYrQ60Qrmt5dDqQCsAgDI+BsqYCMh09KBlZA6bPmtgWr09DIpchBCpNCKA9daFoBz5NOcCJyLKgzR6gaVTp06hc+fO6NKlCxwcHHDu3DkAgLa2NmrWrCltcBrqiv8V3H1xFx/2fsDje49xutppdO3aVeqwiIiINF5cXBz8PvsBAHzDfWF04wT8/1qNAhWawKxGZ1g07JdYUB4DyBKSDaqqyIBYRRyiEqJQQLdAjsRORJQXaXSy6ufnhxo1asDPzw9TpkxRHS9QoAD+/vtv6QLTUEIIDJ48GM93PodxeWM4TnfExucb0bB0Q8hS7VGJiIgIAH44+gMAIC4wDm83vEXchzkwqz0YRmXrqhdU6iPK71fIdT8BABwLFsCg+qXUitgZ2zFRJSLKJI1OVvv164d+/fpJHYbmiIsG/p0IBL8ELEoCjWYAugYQQuDly5cI0A/AR3yE3XA7GDkZAQAefXqEK/5X4FnUU+LgiYiIJKZUAK+uABEfgAKFALsagFwLwP9GVd/5QaYrg5aRFozKGMFuuBWi/GoipbumlDHFoYwpDgAwt7ZEy5LVcvKVEBHlCxqdrNIX/vwJePrF0vc+/wE3N+CBQQ0MOxWHV69eodLCSrCsbQmlUK1BCLlMjuV3lqOGTQ2OrhIRUf7lfQQ4OQYI9///YyY2QNO5iCvVFG593fB873MU6VIE5rXMUejHQgAAw5JLEPVydJpNuxU1zc7IiYjyLSarucHXiSoSp/wO/ycWa2+dxKAmThh9dCVG3RiVrKpSKDm6SkRE+U6/rTfhHfAZAFBbcRUz4+ZB4Ksx0vD3uLWkKzr/Y4w3ymAU/604CrioT93V0gtGZQdj6CLlrd7cippibPMy2fMiiIjyOSarmi4uWi1RTVAKPA1SomxBLZQvJMfd/kYobfUeP73YCBlkEEi+NKEMMo6uEhFRvtFv602cehwIAJBDicF6GxIT1S+6wKdBCtiZyVHEWAbdWkqUql8KMq2U+8hIi9nY1+5EDkRORERf0vita/K9fyeq/vf0ywRUWBOJX4/FQAiBnyvqwslKC/EAAkJ8UkxUAUBAICAyAPHK+BwKmoiISDpJI6oAUEX+BDayYFWiGhojMOKfGFRYG4kLrxSwMpZB3tAq1UQVAN5GvEVcXFx2h01ERF/hyKqmC36JBKVAh73ROOubgCl19TDIQ1dthFQXwG4UQXDLeak2Y6FvAV2tlKcwERER5QkRwcDW5vg7+i3e65igQ/xkFESo6vQ/LxLQ7WA0KhSWw6uvEcoW1IJPOpt+E/EGDhYO2RI2ERGljMmqBvv8+TNehhqjvFyGNk7aWNdSH9ZGKQ+GF7Z0QmFLlxyOkIiISEPMdwQiE6f+mgEw0/qMB/J+CFEa4drbBHjYaKGkuQwbWumjtZO26qKvA4C+IWG4X6IqYGCWYtPlrcozUSUikgCTVQ2kVCqxbds2jBs3Dg3r18N2R6BnhW+MijaakTPBERERaYrgd8C6GkBMaIqnX4UqMfr0R5x4kYDzvYxQqYgWHC21violwxBhDLT+S7WNDRERaQbes6phIiIiUK1aNUyePBnLli3Dth07AafmaVdyag7oGuRMgERERJpgmjWwzAVXEY0fihbGVX09tdNb78ah7OoIGOnI8HRQAVQorAXlV0s7qDZ6azqHiSoRkQbiyKqGePfuHYKCglC+fHkMHz4cbdq0gYHB/xLQn/5McfsaAImJ6k9/5mywREREUppmjVaFLeCnqwNdIRAnl2OphRmqvgvAKR8FmjhooVIRLZzraQSPoolJ6FuFBeRywAbBqmY+wgqFOi4GXFpL9UqIiCgNTFZz0OTLk3HgxQG0LdUWUz2nAgCio6OxaNEizJ49G8OHD0f58uXx008/Ja/805+J29j8OxEIfglYlEyc+ssRVSIiyoCYmBjo6+tLHcb3C36HykWtEaslB2QyxP3v3lOv90pU2BqDT8EJuNbbCOUKqY+UGshi4RG7FlXkT1AQoQiEGfxNKuCCSyMpXgUREaUDk9UcolAocPDFQQDAwRcHManaJHz48AGenp6wtLTEP//8A09Pz7Qb0TUAWizMgWiJiCivmT17NhYsWICwsDDY2tpi8eLFaNOmjdRhZVj0uhqILWwMfLEqfuCRQHz8+yPKNDTHk/LxMNZNvg1NoDCDEnJcU/7/YoSNi5jlRMhERPSdeM9qDpl0dZJqH9ToN9H4ddevKFKkCJYsWYIbN258O1ElIiL6Tps3b8aMGTOwZ88eREdHY9iwYejQoQPu378vdWgZ5lnQCJDJoIxXIuxmGADAuLwxHGc6Qt7RBvdN1Wccif/9DNedhWLmBqqfxmUKYl1Pj5x/AURElG4yIYT4drHcKzw8HKampggLC4OJiYkkMSgUClTcURHxn+MReCAQIZdCYNPVBq82vIKWFhd0ICLKTTShX8mo8uXLo1q1ali7dq3qWNmyZVGrVi2sWbMmXW1owuuOjo6Gx+7K+HwvAgG7AyDXlaPE+BLQMkjsS+VCoExcHP70/wC1sVWjgsDvzyWJmYiIkktvn8KR1Rww6eokRL+NxrMxzxAfEo9S00vBvK45Jl2dJHVoRESUx0VFReHBgweoVauW2vE6derg+vXrEkX1fTz3eMJ/63u82/gOVk2s4DDFQZWoAoBSJsMjPT1cMfjinlwmqkREuRbvWc1mJ06cwL77+6Bno4fivxVHgTIFVOeO+hzFtOrTOLpKRETZJigoCEIIWFlZqR0vWLAgAgMDU60XGxuL2NhY1ePw8PBsi/FbQkJCcPjwYcTL4mFR3wKFOxaGllHKfadMCCw3N0WNAvaQ9TwBFLDI4WiJiCircGQ1mzx//hytW7dG245tEfM+BjK5TC1RBQABwdFVIiLKEUqlUu1xQkICZLLkCxElmT17NkxNTVU/tra22R1iMgqFAmvWrIGjoyOGLBoCkSBgYGeQaqIKAEImQ4CpDeIHXGSiSkSUyzFZzQZnz56Fm5sbihYtihKzS8CkYurzsI/6HIVCocjB6IiIKD8pVKgQtLS08OHDB7XjgYGBsLGxSbXeuHHjEBYWpvp58+ZNdoeaTLNmzbBw4UKsXbsWxUcUh0w79eQ6ycJaC7G75W7oaunmQIRERJSdmKxmEaVSiU2bNuHdu3eoXr06rl+/jpmLZkLLOO0pvgICYXFhORQlERHlN3p6evDw8MCZM2fUjv/777+oWbNmmvVMTEzUfnKCr68vtm3bBgCYP38+Hj58iBLVS6Q5CpyktX1rNC7ZGIWNCmd3mERElAN4z2oWuHz5MoYOHYpPnz6hbNmyKFq0KNzc3AAAy+otw7PQZ6nWdTJzgoUBpykREVH2GTduHNq1a4d69eqhfv36WL58OYKCgvDbb79JHZpKREQE5syZg0WLFqFHjx7o3r07ypcvDwCoZFMJjWwb4XHI41Trl7Uoi5l1ZuZUuERElAO4dU0m7dq1C/369cO4ceMwYsQIGBgYfLsSERHlWpqwhcv32LVrF+bPnw9/f3+UKVMGc+bMQbVq1dJdPztft0KhgKurKywsLLB06VJUrlw5S9snIiLNkt4+hcnqd4iKisLSpUsxYMAAaGlp4fPnzyhatGiWtE1ERJottyarmZUdr/vGjRt4+vQpunfvjidPnsDJySld032JiCh34z6r2UAIgT179qBMmTI4cOAAgoKCYGJiwkSViIgoA96/f49evXqhbt26ePfuHQDA2dmZiSoREanhPasZMH/+fCxevBhz5sxB9+7dIZcz1yciIsqIyMhIuLm5oX79+nj8+DHs7OykDomIiDQUk9UM6NOnD3799VcYGxtLHQoREVGuZGRkBC8vLyapRET0TUxWM8DCgqv2EhERZRYTVSIiSg/OYyUiIiIiIiKNw2SViIiIiIiINA6TVSIiIiIiItI4Gp+sBgYGokePHihatCgcHBwwfvx4xMfHSx0WERERERERZSONXmBJqVSiZcuWMDAwwD///IOQkBB07twZ4eHhWLFihdThERERERERUTbR6GT19OnTuHnzJp4/f45SpUoBAGbNmoU+ffpg2rRpXJ2XiIiIiIgoj9LoacAXL15E8eLFVYkqADRq1AgJCQm4du2ahJERERERERFRdtLokdV3796hUKFCascKFiwImUwGf3//FOvExsYiNjZW9TgsLAwAEB4enn2BEhFRvpFf+xMhBID8+/qJiCjrJPUlSX1LajQ6WQUAuVye7LFMJkv1hc2ePRtTp05NdtzW1jZb4iMiIsoPPn/+DID9KRERZZ3Pnz/D1NQ01fManawWLFgQFy5cUDsWFBQEpVKJggULplhn3LhxGDFihOqxUqnEq1evUKFCBbx58wYmJibZGnNuEB4eDltbW74f/8P3Qx3fD3V8P9Tx/fj/q8DGxsYSR5KzbGxs8ObNGxgbG0Mmk+X48+f33738/voBvgd8/Xz9een1CyHw+fNn2NjYpFlOo5PVqlWrYt68efD391e9kPPnz0Mmk8HDwyPFOnp6etDT01M7ljQ6a2Jikic+3KzC90Md3w91fD/U8f1Qx/cj/5HL5ShWrJjUYeT73738/voBvgd8/Xz9eeX1pzWimkSjF1hq0aIFHBwcMHToUISHh+PNmzeYMmUKOnTo8M0snIiIiIiIiHIvjU5WdXV1cezYMQQEBMDS0hIODg5wc3PD+vXrpQ6NiIiIiIiIspFGTwMGgNKlS+PixYuIi4uDlpYWtLS0MtyGnp4eJk+enGx6cH7F90Md3w91fD/U8f1Qx/eDpJLff/fy++sH+B7w9fP158fXLxPfWi+YiIiIiIiIKIdp9DRgIiIiIiIiyp+YrBIREREREZHGYbJKREREREREGkfjF1jKCm/fvkVAQAAcHR3TtZ9PXpaQkIBnz55BT08P9vb237VgVV708eNHPH36FPb29hqxj6CUYmJi8PjxYxQuXBhFihSROhxJxcXF4eXLl4iJiUHJkiXzzL5m6RUREYF79+6hePHisLW1TbFMcHAwfHx8UKxYsXz/+0LZ6/nz5/j8+TNcXFygr68vdTjZSgiB27dvQ0dHB25ubimWiYuLw6NHj2BgYABnZ+ccjjB7RUVF4dmzZ7C2tkbRokVTLCOEgLe3NxISEuDq6pqnvs8olUq8ePECCQkJKFmyZKq/7z4+PggNDUWZMmVgaGiYw1Fmv6ioKNy+fRsFCxZE6dKlk51///493r59CwcHB1hYWEgQYdb79OkTHj9+nOx41apVoaOjo3YsJCQEPj4+KFKkSKr/TvIEkYfFxMSI9u3bCwMDA1GmTBmhr68vFi9eLHVYkkhISBCTJk0S1tbWwsXFRdja2gp7e3vx77//Sh2a5OLi4kTVqlWFXC4Xs2fPljocSS1ZskSYmJgIFxcXUapUKdG9e3cRGxsrdViSOHDggChUqJAoWbKkKF++vDAwMBBjxoyROqwc8e7dOzFo0CBRpEgRoaurK6ZPn55iualTpwo9PT3h4uIi9PT0RK9evURCQkIOR0t53YcPH0S1atWEmZmZKFWqlDA3NxeHDx+WOqxsoVAoxPz580WpUqWEmZmZqFq1aorlTp8+LaytrUWJEiWEpaWlKF++vHj9+nUOR5v13r9/L3r27ClMTU1FhQoVhLm5ufD09BS+vr5q5R4/fixKly4tChUqJIoVKyaKFSsmrl+/Lk3QWWzTpk3Czs5OuLi4CEdHR2FqaiqWL1+uViY4OFjUqVNHGBsbC0dHR2FiYiJ2794tUcTZp3v37kIul4uuXbuqHU9ISBC9evUS+vr6qv5n6tSpEkWZtfbu3Su0tLSEp6en2s+nT5/Uys2ZM0fo6+ur8psuXbqIuLg4iaLOXnk6WZ04caKwsbERb9++FUIIcfToUQFAXL58WeLIct7nz5/F1KlTRWhoqBBCCKVSKUaNGiVMTExESEiItMFJ7Pfffxc9e/YUhQoVytfJ6oYNG4Senp44c+aM6tiuXbtEcHCwhFFJIy4uThgZGaklp6dPnxYAxNmzZ6ULLIdcuHBBLFu2TISGhgo7O7sUk9W///5baGtriwsXLgghhHj27JkwMzMTCxcuzOlwKY9r06aNqFy5soiMjBRCCDF79mxhaGgo3r9/L3FkWS8qKkqMHDlSPH/+XAwaNCjFZDUkJESYm5uL8ePHCyGEiI2NFbVq1RL16tXL6XCz3LVr18TWrVtFfHy8ECLxu0vt2rVFrVq1VGWUSqVwc3MTbdq0EQqFQgghRO/evYWtra2IiYmRJO6stGTJEhEYGKh6vHXrVgFAeHt7q45169ZNuLq6irCwMCGEEMuXLxe6urrJkvrcbNu2baJq1aqiTp06yZLVBQsWCAsLC/HixQshhBDnzp0TWlpa4tixY1KEmqX27t0rTE1N0yxz+vRpIZfLxenTp4UQQvj6+gorKysxc+bMHIgw5+XpZNXGxkZMnDhR7ViFChVE7969JYpIs7x48UIAEOfPn5c6FMmcPHlSODg4iPDw8HydrCoUClG0aFExZMgQqUPRCKGhoUImk4lDhw6pjkVFRQmZTCb27NkjYWQ5L7VktW3btqJhw4ZqxwYMGCDKli2bU6FRPvDx40chl8vVRo2io6OFsbFxnp8plVqyumnTJqGrq6tKVIRIvHgEIE8lK0k2bNggtLW1VYnpjRs3BADh5eWlKuPn5ycAiKNHj0oVZrZ5/fq1AKBKTD5//ix0dXXFhg0bVGUSEhKElZVVqrNgcptnz56JwoULi+fPn4sGDRokS1ZdXFzE4MGD1Y7VrVtXtGvXLifDzBZ79+4VJiYm4tGjR+LBgwciOjo6WZkuXbqImjVrqh0bNmyYcHBwyKkwc1SeXWApMDAQ/v7+cHd3VztepUoV3LlzR6KoNMvNmzchk8lQsmRJqUORxIcPH/DLL79g27ZtMDY2ljocST19+hTv3r1Dq1b/196dR0VxpX0A/jWoyAiyiCyCdIgLqERHMYACyrC5iyxCXFqjiJoZ9WiCymQ06rjGxJioMXHABSdRQhgElcW4gNqiBAQGRRDEiEJAbRBE9qW+P/i6xoLupjFAQ/s+53AOfWt7b1V1vXWrb1XNgkgkQmpqKkpKShQdlsJoaWkhMDAQmzZtQnh4OOLi4rBw4ULY29tj1qxZig6vW0hLS5N4fM3KykJtba2CoiLKJiMjA01NTZx9rW/fvnjvvffe2lyelpaGYcOGce6ht7a2Zocpm+TkZJiZmUFFpfmUNS0tDSoqKhg7diw7Dp/Ph76+vtLU/+nTpxAKhYiKisKHH34INzc3TJ48GQBw79491NXVcb4TqqqqGDdunFLUv66uDr6+vti5cyeGDh3aarj4uRrKfH7/8uVLuLu7w93dHbq6uti1axdnuLT8m5eXh4qKiq4MtUso7QOWSktLAQADBgzglA8YMIAd9jb7/fff8fHHH2PJkiVv5QOFGIaBQCCAv78/Jk6cqOhwFO73338HAJw7dw4LFy6EkZERsrOz4ePjg+Dg4FY39b8N5s2bh0uXLiEgIACampoQiUQ4ePCg0j/YRV6lpaUSj69NTU0oKyuDgYGBgiIjyoRyeWuSvnvih8so2zq5fPkygoODcfz4cbastLQU2trabONVTJn2ieTkZOzZswdFRUWoqanB4cOH0atX8ym7rO9EUVFRl8fa0davXw8zMzMsXbpU4vCysjIwDKO0xwQ+n4/U1FT2Ysz58+cxZ84cDB48GAKBAID0/Csepmw/wCjtL6vik+uamhpOeXV1Nfr06aOIkLoNkUiEKVOmwNzcHIcOHVJ0OAoRFBSEe/fuwdHREUKhEEKhEPX19cjPz0dSUpKiw+ty4u9LWloaHj58iLS0NGRkZCAqKgpfffWVgqPreiUlJZg0aRKmT5+OR48e4c6dOwgLC8O8efMQGxur6PC6hd69e0s8vgJ464+xpONQLm9N0ndP/FmZ1klycjI8PT0REBDAnqQDkusPKNc+MXPmTAiFQuTl5WH//v3w9PSEUCgEoNzfiRs3biAoKAiLFy9mz83Ky8vx/Plz9jxNmesPAO+//z6n18DMmTMxc+ZMhIaGsmVvW/5V2saqsbExVFVVUVhYyCkvLCyEqampgqJSvJKSEri4uEBXVxfnz5+Hurq6okNSCB6Ph3feeQebNm1CYGAgAgMDUVFRgdjYWGzfvl3R4XW5d955BwAgEAjYx98PGzYMzs7OuH79ugIjUwyhUIiysjL87W9/Y8scHBwwevRonDt3ToGRdR98Pl/i8VVTUxM6OjoKioooGz6fDwCUy18j7bsHQGnWSUpKCtzc3ODv7489e/ZwhvH5fFRVVaGsrIwta2howLNnz5Sm/q/z8fGBqakp4uLiACj3d6Kmpgbjxo3D3r172XOz3NxcpKWlsedpurq60NTUVMr6S2NgYMCpr7RjgJqaGvT19bs6vE6ntI3Vvn37wsHBAWfPnmXLqqurcfHiRbi6uiowMsUpLS2Fi4sLtLS0EBMTg379+ik6JIXx9/dnr9qJ/3R1dbFy5UqcP39e0eF1OVNTU1hYWEg8+A0cOFBBUSmOuM4FBQVsWUNDA4qLi9/K9SGJq6srYmNj0dDQwJZFRUXBxcVFgVERZfPee+/BwMCAk8tzcnKQlZX11uZyV1dXFBYWIjU1lS2LiopC//79YWNjo8DIOsbt27fh6uqKpUuX4ssvv2w13NHREb179+bsExcvXkRVVVWPP/7U1tZyjqlA8/2Lz58/Z7t5Dh06FGZmZpz6FxQUsOutJ3N2dm51bjZ+/Hi4ubmx52k8Hg/Ozs6c+tfX1yMmJqbH1x8AKisrOZ/r6+tx9epVWFpasmWurq64cOEC6urq2LKoqCg4OTkp1fuGxZT2nlUA2LFjB/7yl79gw4YNsLe3x7fffgttbW2sXLlS0aF1uerqari6uuL58+fYs2cP5yZ0CwsL6OnpKTA60h188cUXWLBgAXR0dDBixAhERUUhPT0d33//vaJD63LW1tawsbHBggUL8Nlnn6F///4IDg5GRUUFFi9erOjwOl1NTQ1SUlIANJ88PX78GEKhENra2mzCXLt2LU6cOAFfX198+OGHiI2NRXJyMm7evKnI0ImSUVVVxa5du/DRRx9BR0cHpqam2LZtGyZPnoxp06YpOrxOkZqaiqqqKhQVFaGiooLt/mlnZwcejwc7OzvMmjUL8+bNw/bt2yESibBlyxbs3Lmzx99Tn52dDTc3N4wdOxYeHh5s3QHAxsYGvXv3hp6eHgICArB27Vo0NDRATU0NGzduxNKlS2Fubq7A6P+4oqIieHl5wc/PD8OGDcPTp09x4MAB6Ovrc3LP7t27sXDhQujr68Pc3By7du2ClZUVvLy8FBh919myZQsmTpyIVatWYcqUKTh27BgaGxuxdu1aRYf2hy1evBhDhw7FxIkTUVdXh8OHD+P58+fYtGkTO87q1asRHBwMb29v+Pv74/Lly7h69arS9oTjMQzDKDqIznTz5k0cOnQIxcXFsLS0RGBgIIyMjBQdVpcrLi6Gt7e3xGHbtm2Ds7NzF0fU/bi7u8PHxwcLFixQdCgKk5CQgCNHjkAkEmHYsGFYs2YNLCwsFB2WQrx69Qrffvstbt26hZqaGlhYWGD16tVvxdOzCwsL4evr26r8/fffx/79+9nP+fn5+Pzzz3H//n2YmJhg7dq1nHttCOkokZGROHnyJCoqKmBnZ4f169crbe+gBQsWID8/v1X5lStX2PvRampqsH//fsTHx0NdXR0ffPAB5s2b19WhdrjY2Fjs3LlT4rCzZ8+yD5JiGAZHjx5FZGQkGhoaMHXqVKxatYp9CFFP9uDBAxw+fBiZmZnQ1taGra0t/P39oaGhwRkvNjYWR48eRVlZGWxsbLBhwwZoaWkpKOrOs27dOujq6mLz5s2c8rS0NHz99dcoKCiAubk5Nm7cyHaR7slqa2tx5MgRxMfHg2EYjB49GmvWrGn1o1JBQQH27NmDrKwsDBo0CGvWrMH777+voKg7l9I3VgkhhBBCCCGE9DxKe88qIYQQQgghhJCeixqrhBBCCCGEEEK6HWqsEkIIIYQQQgjpdqixSgghhBBCCCGk26HGKiGEEEIIIYSQbocaq4QQQgghhBBCuh1qrBJCCCGEEEII6XaosUpkKisrQ2hoKOrr6znl6enpCA8Px82bN+WeV05ODi5cuNDRIXaKe/fu4dKlS4oOo0Pj6Og6ZWRkICEhoV3LbOtzV7t58yays7PlHv/27du4c+eOzHFSUlIQGhqK0NBQFBYWApBvXcmjo+bzupycHDberKysDp03IUTxlDWPd0b+UHRO+iMk5R55vEleoZxGulIvRQdAurdHjx5h3rx5ePHiBbS1tQEAgYGBCAkJgb29PRwcHKCpqYmnT5/C2dlZ5rxiYmJw4sQJTJkypdWwpqYmREREQE1NDbNmzXqjWO/evStXHPJMFxERgfPnz8PFxeWNYukoHRlHR9fp1KlTSElJgaOjo9zLbOvzm27DN1FcXIzZs2cjJSWFU/7gwQPcvXsXBgYGsLGxgYrK/67p1dTUwMvLC1lZWVBXV5c43+DgYERHR8POzg5Dhw6FsbGxXOtKHn9kPmfOnEFdXR309PTg4OCAPn36AADy8vIQGRmJX375BRs2bMCIESP+UIyEkO6lJ+Xx9uiMPN2VuT87Oxs5OTno168fLCwsYGxs3GqchoYGJCUloby8HLa2ttDV1ZU6P0m5Jzc3F7dv35Y4vo+PD1RUVN4or3SHnFZTU4MvvvgCmZmZ0NPTw6FDhwBQTlNG1Fgl7XbkyBEcO3YMHh4eAICtW7fi0qVLf6iBsXv3bvzzn/+EkZHRGye58PDwN4pD0nSjRo1CbW3tG8XRXSmiTm0ts+XwN92Gb2LHjh2YPXs2+Hw+WxYYGIhDhw5h4sSJyMzMhKmpKeLi4qClpQUAsLOzg6mpKb7//nusW7dO6rxtbGwQGhra6XVoj9jYWJSXlyMvLw99+/aFUCgEAEybNg3Tpk3D+PHjFRwhIaSrdNc83h6dkdO6Ik+mp6dj6dKlePLkCaytrdHY2Ijs7GyMGTMGQUFB0NfXB9D8S6mPjw/69OmDESNGYN26ddi3bx9mzpwpdd4tc8+FCxewbt06eHl5tRrXy8uLczG2pxEIBFi0aBE2b96Mly9fsuWU05QPNVZJK2lpaSgsLGx1NaqkpAQXLlxAWVkZ0tLSUFtbC2NjY9y9excikYg9QE6YMAF8Ph8NDQ24fv066uvrMXbsWKnLS0xMRHBwMFauXImoqKg24ysqKkJ6ejo0NTVhZWUFdXV15ObmSoyjsbERv/76KwBAQ0MDo0aNgpmZGTsvadOZm5tDU1Oz1bIfPHjANmJk1SkjIwOlpaWwsbFBeno6RCIR7O3toaOjwxmPYRikpKTgyZMn4PP5sLKykln3hw8fyqzP68u2tbVFYmIiysrK4OnpyalTU1MTwsLCWs2/d+/e8PLykms5YtXV1VLrKG09ShouaVuYmZnh0aNHmDt3LiepPnnyBImJia3KxRoaGnDjxg1UVlbCzs4OWlpauHHjBvT19TFs2DC8evUKISEhiImJYae5dOkS9u7dC6FQiIkTJ6K8vBxWVlbYtGkTDh48yI4nEAiwd+9emY1VaWStK0D+/aG2thbp6ekoLS3F+PHjMXDgQHaYtO1/5MgR/PTTT7hz5w4yMjLaHTshpOfoznlcnvwoTw5qmV8kHftcXFwQExMDDw8PqKmpAQCio6OhqamJSZMmAQBEIhEuX76MuXPnSsxZks45Xnf//n1kZWXByMgIY8eOZXutSJKXlwdHR0d4e3sjMTERffv2ZYdFR0ejsrISQPN2mj59OubOnYtDhw6Bx+Ph1atX7DppDzU1NbkunsrKK9IoKqc5OzsjMzMTaWlpiIqKwtSpU+Ht7d1mvKRnosYq4fDz80NYWBjs7OyQk5MDCwsLdtiLFy/YJCQUCpGTkwN7e3tkZ2ejpKQEkZGRAAATExNoaWnB2dkZxcXFGD16NDIyMjB8+PBWyysrK8PChQtx7NgxXL16tc34jhw5goCAALYh+vvvv+PUqVN4/vy5xDjq6+vZzy9fvsT169exYsUKfPnllwCA3377TeJ0CQkJnK5AdXV1+PDDD3Hu3DlMmDABJSUlMDIywtmzZyU2lk6dOoWIiAioqKhg0KBBKC0txePHjxEdHY0JEyYAACoqKjBz5kzcv38f48aNQ0pKCsaMGYOoqCj86U9/klj//Px8mfURL/vMmTPo1asXjIyMYGJiAk9PT073psbGRnY+YsnJyRCJRPDy8pJrOQDw+PFjjB07Vmod2+pS9fpwSdtCIBBg8eLF0NHRgZubGzvdrl27cO/ePfj6+raa57NnzzBlyhTU1NTA0NAQS5Yswa1btzBv3jyEh4cDABISEtDY2AhbW1t2uh9++AG2traYOHEiAEBLSwt+fn744osvcODAAfB4PACAk5MTli9fjuzsbM73oy1trSt59wfx1XgDAwOUl5cjOzsbcXFx7HykbX+GYRAREYHMzEz4+PjIHTchpGfp7nlcnvwoTw5qmV8kHftmzJiBJUuWQF9fH05OTqisrISHhwcMDQ3x+PFjAMB//vMffPXVV/D19W01T2nnHGPHjmXPCy5fvgxra2vk5+ejvr4eZ8+exbBhwyTWffv27dDQ0MDBgwc5DVUAmDFjBvt/UFAQ281VnHs0NDTg5OTU5vp9E23lFUkUmdP+/Oc/48GDB/Dx8cG7774LV1dXDB8+HKNHj+6U9UMUjCHk/0VHRzN9+vRh7t69yzAMw1RXVzO2trYMAObFixcMwzBMfX09A4CJj49np9uyZQtjZ2fHmdf69euZESNGMOXl5QzDMMyjR48YLS0tZsyYMZzxPD09mfXr17Pz4fP5MmM0NDRkTp48yX4uKChgkpKSpMbRUl5eHqOhocHcunVLZvzbt29nbGxs2M9///vfGUNDQ+a3335jy6Kjo5m6ujqJy9m4cSMDgDl9+jRbtmzZMsbS0pJpbGxkGIZhAgMDmSFDhjAikYhhGIYpLi5mTExMmG3btkmNQ576iJd98eJFmXV6XWpqKqOurs4cOHCg3cuRVceWy2zrs6Rt8cEHHzC+vr7s5+rqakZbW5sJCQmRGKtAIGDee+89pra2lmEYhvH392cmTZrEDB48mGlqamIYhmG2bdvGjB49mjPduHHjGH9/f07ZuXPnGADMkydPOOXq6uqc/fB1K1asYLy8vDhlHbU/bNy4kVFRUWFu3LjBWT/Tp09vtayW21+sqqqK6dWrF1NWVsYpt7KyYnbv3i1xGkJIz9AT8rg8x8OWJOWglvlD2rHP0dGR2bx5M8MwDBMXF8eYm5szmpqaTF5eHsMwzcfQ5cuXS5ynrHOOzz77jLG2tmYqKyvZ4StXrmScnZ2l1t3Q0JBZtmyZ1OFi06ZNY6ZOncq8fPmSOX/+PHPlyhV2+0kjKfccPHiQUVNTY06fPs35i4uLY8eRJ6+0pOicVltby5iamrI5/aOPPmLOnTvHiZFymvLouZ3VSYcLCwvD9OnTMWrUKABA3759sWbNmjee14oVK9C/f38AAJ/Pb/Ur2HfffYeHDx9ix44dcs9XXV0d9+7dQ0NDAwDA2NgY1tbWMqeprq6GUChEeHg4fv31VxgaGra7K01ISAhWrVqFd955hy2bPn06evfuLXUaExMTfPDBB+znDRs24O7du+zT6UJDQ7FixQoMGDAAAGBgYIBly5a12V1HnvqMGDFC7gdEPHv2DO7u7pg/fz5Wr17druW0VceOsGzZMkRGRuLFixcAmq+CNzU1Se3yEx0djcWLF7Ndsby9vXHt2jV4enqyV6hFIlGr7krl5eWtHl4h3jZlZWWcch0dHYhEonbVo6P2hzFjxrC//gKAo6Mj7t+/zxmn5favrKxEaGgofvzxR3zyyScYNGgQ+90khCiPnpDHAflyx5vkbkm5z9HREfHx8QCae9W4uLjAzs6OLbt69arUB/zIOuc4fvw4Ro8ejZiYGPz8888ICwuDvr4+rl+/zo7f0tOnT2FqaiqzDkBzXq6oqMC4cePw7bffYuPGjTAzM2vVI0oe4p5Ur/+J6y4mT15pSZE5rU+fPggICIBAIMDWrVuRm5vbab86E8WjbsCE9fjxY4wZM4ZTJu0+RVkYhkFBQQGnYSeeV1JSEoDmrkgff/wxNmzYgIiICADNT4IVn1Tb29vDxMSk1byPHz+Ov/71r/j+++/h4OAALy8vCAQCqQ8JiI+Ph4+PDwwNDfHuu+9CXV0d5eXlePbsmdz1EXf9kdT9SRZJ9QeauzeNHDkST548wbvvvssZZ8iQIcjPz5c6T3nrY2RkJFeM9fX18Pb2xuDBg3H48OF2L0dWHcUnS3+Uk5MTjI2N8eOPP2LVqlU4duwYfH19JXaVbmxsRGlpKScu8T7t6enJlmloaLD3Bompqanh1atXnDLx55bdtSorK2XeiytJR+0PLRvUampqqKmp4ZS13P5VVVWIjIyEqqoqBg8ejISEBLbhTghRHj0hjwNt5443zd2Scp+joyN27dqFqqoqJCQk4JNPPsHgwYMRHx8Pe3t7FBUVSW2sSjvnYBgGT548QXZ2NioqKjjTeHh4oLq6WmKO0NDQQElJicw6AM0559atW0hOTmbvFf7000+xePFiFBcXS30avSTy3LMqT15pSZE5DQBWr14NgUCA0tJSbNmyhXKaEqPGKmENGDCA/fVKrOVnefB4POjo6MicF8MwcHd3x/3799kraNnZ2exJtZmZmcQkN3nyZGRmZiI/Px8XLlxAYGAgMjMzsXfvXomxfPLJJ1i+fDl27tzJlllaWoJhGLnro6qqKneCeZ20+uvp6YHH40FXVxelpaWccUpLS6Gnpyd1nvLWR96D9urVq/Hbb78hJSWF81AIeZcjq44dhcfjwc/PD8eOHcPMmTMRHx/Piet1qqqq0NbW5pw8PHz4EADYJ/oCwPDhwxEUFMSZdsiQIew9TGL5+fno1asX50r4ixcvUF5eDnNz83bVozP2B2labv+BAwd2u6cTE0I6Xk/I45Jiapk73jR3S8p9tra2UFFRwS+//ILU1FRMnjwZgwcPxsGDB+Hg4IDhw4dLvcAr65yjX79+mDt3brt+uba2tsatW7faHG/IkCF4+vQp56FWXl5e2L17N3Jzc7vFvZmKzGli2tra7OuYiPKibsCEZW9vjwsXLnCuaImvlsqioaHR6iqYvb09p7tKY2Mj5wmBurq67EubxX/e3t7sSbWNjU2r5TQ1NaG4uBhAc3ek5cuXY9GiReyBX1IcxcXFnEZFZmYmsrOz24y/JTc3N/zwww+cRFlSUoKmpiap02RlZSEnJ4f9HBERAW1tbYwcORJA8zpquX7Dw8Nhb28vdZ7y1Edehw8fxsmTJxEZGQkDA4M3Wk5bdWwvadtiyZIlyMjIwLp16zBixAjOg5FasrKy4rzU/YcffgAAzrvmnJycUFJSwulqNH36dFy+fJmTXMPCwuDi4sJpyN+4cQOampoS91FZOmN/IISQ13X3PC7W1vGwI3OdmpoabG1tsXPnTgwfPhwDBw7E+PHj8fLlSwQHB0v9VbWtc46pU6ciODi4VZffwsJCqbGsW7cOSUlJOH36dKthFRUV7C0ns2bNQlFREcrLy9nhWVlZ4PF4Ui8AdDXKaaSr0C+rhLV8+XIcOnQIzs7OWLRoEdLT0+W6P2L8+PH4xz/+ga+//hqGhoaYMGECtm3bBltbWwgEAjg4OCA8PBwlJSUSX3otr8bGRjg4OMDNzQ3jxo1DSUkJgoODsXXrVqlxzJkzB5s2bUJFRQWqqqrw9ddfo1+/fm3G39LevXthb28PJycn+Pj4QCQS4eeff8bt27eldkHu378/ZsyYgTVr1qCkpASff/45du/eDQ0NDQDNT7S1sbGBr68vnJ2dERsbi7t37+LEiRNS14E89ZFHfn4+1q5dCw8PD+Tm5iI3NxfA/15dI+9y2qpje0naFnw+H0ZGRpgxYwYiIyNbPZG4pU8//RRubm7Q09ODkZERTp48CT8/P+zbtw+mpqZwcXGBqakp3NzccPr0aXb/WbJkCYKCguDq6go/Pz8kJiYiMTGRfR+p2E8//QSBQCDz9QSSdMb+QAghr+vueVysreNhR+U6MUdHR2zdupV9LoOqqirs7OwQFxeHjz/+WOI0bZ1z7Nu3Dw4ODpgwYQIEAgF4PB6uXbsGHo8n8dVwQPMTf/fu3YslS5YgJiYGkydPRkNDA7KysnDu3DnExsZCW1sbXl5e+Ne//gU3Nzf4+flBJBLhyy+/xPr161t1m21LQ0ODxJ41M2bMaPftLK+jnEa6CjVWCUtdXR03b97EgQMH8Ouvv2LkyJFISEjA9u3b2RNzFRUV+Pr6si+tBpqTwL///W9cuXIFSUlJMDExgb29PZKSkhAUFIT09HQsWrQIa9euxbVr16Qu39LSUuaLxHv37o3//ve/CAkJwa1bt9CvXz+EhobC1dVVahzffPMNLC0tkZycjP79++Onn37ClStXOI/ylzRdyxeDDxkyBHfu3MHRo0eRlJQEMzMz/PLLLzIfsGRlZYXt27fjzJkzEIlE+PHHHzn3TVpYWCA9PR3BwcEQCoWwtLTEN998w+ly2jIOeeozZswYiQmo5bzEsbx+ItO3b194eXnJvZyAgAA4OTlJrWPLZbb1WdK24PP5AJpPXmJjYyEQCKSuc6D5V9OEhAScPn0a2dnZOHPmDCZPnozPP/8cQqGQfUjD5s2b4ePjg40bN0JdXR1qamq4du0avvvuOyQlJWHQoEFIS0vjvIKgqKgI0dHRnF9pJSkoKEBoaCgcHBxgbGws17qSZ3+QtG3NzMzg7u4ucxxZcnJykJqa+kZdBQkh3Ut3z+NibeVHeXJQy/wh69jn7u6OrKwszJ07ly3z8/ODlpYW58E8r8+zrXMOPp+PO3fuICQkBOnp6ejfvz/mz5+POXPmyKz7+vXr4eHhgbCwMNy4cQP9+vXDyJEjsW3bNrZLK4/HQ3R0NI4fP46bN29CW1sbp06dwtSpU2XOu2XuGT58OObMmSPxgsWkSZOgqakpV15piXIa6Uo8pj037xFC5BIYGIiUlBROd1Tyx/j4+KCpqYl9V2pH2LFjByZNmsS+HL4t4l8WVqxYIXWco0eP4uLFiwCAgIAAjB8/vkNi7SyxsbEICQkBAMyfPx+zZ89WcESEEGVG+bFz9LTc01kopykfaqwS0gkoGXccoVCI69evY8uWLUhMTHxrEzAhhCgDyo+EkPagbsCEdIL2dlsh0iUlJeH+/fuIiIighiohhPRwlB8JIe1Bv6wSQgghhBBCCOl26NU1hBBCCCGEEEK6HWqsEkIIIYQQQgjpdqixSgghhBBCCCGk26HGKiGEEEIIIYSQbocaq4QQQgghhBBCuh1qrBJCCCGEEEII6XaosUoIIYQQQgghpNuhxiohhBBCCCGEkG6HGquEEEIIIYQQQrqd/wO9dbuG+1HsOgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 950x340 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nh3 = molecule(\"NH3\")\n",
    "nh4 = nh3.copy(); nh4 += molecule(\"H\"); nh4.positions[-1] = [0.0, 0.0, 1.5]\n",
    "rows = [report(\"NH3\", nh3, 0.0), report(\"NH4+ (q=+1)\", nh4, 1.0), report(\"NH4- (q=-1)\", nh4, -1.0),\n",
    "        report(\"NH3+ (q=+1)\", nh3, 1.0)]\n",
    "show(rows)\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(9.5, 3.4))\n",
    "for q, style in [(-1.0, \"s\"), (0.0, \"o\"), (1.0, \"^\")]:\n",
    "    ref, mine = reference(nh4, q), xnn_d4(nh4, q)\n",
    "    ax[0].plot(ref[\"polarizabilities\"], mine[\"polarizabilities\"].numpy(), style, label=f\"q_tot = {q:+.0f}\")\n",
    "    ax[1].plot(ref[\"c6 coefficients\"].ravel(), mine[\"c6\"].ravel(), style, label=f\"q_tot = {q:+.0f}\")\n",
    "for a, lab in zip(ax, [\"static polarizability α(0) [bohr³]\", \"pairwise C6 [Eh bohr⁶]\"]):\n",
    "    lim = a.get_xlim(); a.plot(lim, lim, \"k--\", lw=0.8); a.set_xlabel(f\"dftd4 {lab}\"); a.set_ylabel(\"xnn\"); a.legend()\n",
    "plt.suptitle(\"NH4 at three total charges: xnn vs dftd4\"); plt.tight_layout(); plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ffdc26b",
   "metadata": {},
   "source": [
    "## Block 3: Two-body (BJ) and three-body (ATM) energies · eqs 18–27\n",
    "\n",
    "Energies, gradients and virials for molecules of increasing size. The ATM term uses\n",
    "$C_6$ coefficients from *neutral* polarizabilities (paper sec II.C); `s9 = 0` switches it\n",
    "off in both codes and must give identical two-body-only energies.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c34be3bd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:42.951068Z",
     "iopub.status.busy": "2026-09-23T19:19:42.950988Z",
     "iopub.status.idle": "2026-09-23T19:19:44.638087Z",
     "shell.execute_reply": "2026-09-23T19:19:44.637258Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          system               N          E [Eh]            |dE|         |dgrad|           |dCN|            |dq|  |dalpha|/alpha        |dC6|/C6\n",
      "            C2H6               8-0.0021913693974552222        8.67e-19        6.78e-20        0.00e+00        1.11e-16        8.88e-16        1.33e-15\n",
      "            C6H6              12-0.00898814444745347        3.47e-18        3.25e-19        4.44e-16        8.33e-17        6.66e-16        1.33e-15\n",
      "             C60              60-0.1967102894599138        5.55e-17        7.05e-18        1.33e-15        2.03e-15        8.88e-16        2.00e-15\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "C60 two-body only:   dftd4 -2.042758827134e-01   xnn -2.042758827134e-01   |diff| 1.1e-16 Eh\n",
      "C60 ATM three-body: xnn E3 = 7.565593e-03 Eh  (-3.70 % of E2)\n"
     ]
    }
   ],
   "source": [
    "rows = [report(n, molecule(n)) for n in (\"C2H6\", \"C6H6\", \"C60\")]\n",
    "show(rows)\n",
    "\n",
    "c60 = molecule(\"C60\")\n",
    "two_body_ref = reference(c60, params={**PBE0_D4, \"s9\": 0.0})[\"energy\"]\n",
    "two_body_xnn = xnn_d4(c60, s9=0.0)[\"energy_au\"]\n",
    "full = xnn_d4(c60)\n",
    "print(f\"\\nC60 two-body only:   dftd4 {two_body_ref:.12e}   xnn {two_body_xnn:.12e}   |diff| {abs(two_body_ref - two_body_xnn):.1e} Eh\")\n",
    "print(f\"C60 ATM three-body: xnn E3 = {float(full['energy_3body']) / HARTREE:.6e} Eh  ({float(full['energy_3body']) / float(full['energy_2body']) * 100:+.2f} % of E2)\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6eae636a",
   "metadata": {},
   "source": [
    "## Block 4: Periodic systems · Ewald-summed EEQ\n",
    "\n",
    "For a periodic cell the $1/r$ Coulomb matrix of the EEQ model is Ewald-summed. The\n",
    "reference code chooses the splitting parameter automatically, sums fixed $\\pm 2$\n",
    "real- and reciprocal-space windows and averages over equivalent Wigner–Seitz images;\n",
    "xnn follows the same conventions, so crystals agree to the same precision as\n",
    "molecules -- energies, gradients **and** the strain derivative (virial), on\n",
    "primitive cells, sheared cells and a small periodic water box.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "be521f67",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:19:44.639707Z",
     "iopub.status.busy": "2026-09-23T19:19:44.639619Z",
     "iopub.status.idle": "2026-09-23T19:20:12.518867Z",
     "shell.execute_reply": "2026-09-23T19:20:12.518182Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          system               N          E [Eh]            |dE|         |dgrad|           |dCN|            |dq|  |dalpha|/alpha        |dC6|/C6       |dvirial|\n",
      "   NaCl rocksalt               2-0.006890954411354103        2.90e-16        3.89e-19        1.15e-14        4.44e-16        2.22e-16        4.44e-16        5.12e-17\n",
      "      Si diamond               2-0.023251471901675425        2.91e-16        8.73e-19        1.33e-15        0.00e+00        2.22e-16        4.44e-16        2.43e-16\n",
      " Si sheared cell               2-0.022196428758707833        3.47e-18        1.21e-17        8.88e-16        0.00e+00        4.44e-16        8.88e-16        4.52e-16\n",
      "water box (9 atoms)               9-0.002941339168787877        0.00e+00        1.18e-18        2.22e-16        1.11e-15        1.55e-15        3.11e-15        1.54e-16\n"
     ]
    }
   ],
   "source": [
    "nacl = bulk(\"NaCl\", \"rocksalt\", a=5.64)\n",
    "si = bulk(\"Si\", \"diamond\", a=5.43)\n",
    "sheared = si.copy(); c = sheared.cell[:].copy(); c[1, 0] += 0.5; sheared.set_cell(c, scale_atoms=False)\n",
    "rng = np.random.default_rng(0)\n",
    "box = molecule(\"H2O\") * 1\n",
    "from ase import Atoms\n",
    "box = Atoms(cell=[6.2, 6.2, 6.2], pbc=True)\n",
    "for _ in range(3):\n",
    "    w = molecule(\"H2O\"); w.rotate(rng.uniform(0, 360), rng.normal(size=3)); w.translate(rng.uniform(0, 6.2, 3)); box += w\n",
    "rows = [report(\"NaCl rocksalt\", nacl), report(\"Si diamond\", si), report(\"Si sheared cell\", sheared), report(\"water box (9 atoms)\", box)]\n",
    "show(rows)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d5d1f85",
   "metadata": {},
   "source": [
    "## Block 5: Matching non-default real-space cutoffs\n",
    "\n",
    "The upstream defaults (60 / 40 / 30 bohr for pairs / triples / CN) are far longer\n",
    "than an MLIP needs. Both codes accept shorter cutoffs; with the same values they\n",
    "still agree exactly. xnn additionally offers quintic **switching windows** at the\n",
    "cutoffs (`switch_width_pair`, `switch_width_triple`) so energy and forces stay\n",
    "continuous in MD -- a deliberate extension, off by default.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2e71cb08",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:20:12.520489Z",
     "iopub.status.busy": "2026-09-23T19:20:12.520385Z",
     "iopub.status.idle": "2026-09-23T19:20:17.195066Z",
     "shell.execute_reply": "2026-09-23T19:20:17.194166Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cutoffs 12/9/10 bohr: dftd4 -8.814949146389e-03   xnn -8.814949146389e-03   |diff| 3.5e-18 Eh\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhwAAAE1CAYAAACstTkfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdvBJREFUeJzt3XdYU9cbB/BvgLA3speiIiBORMFR916469aqVVuprdVfta1t1VZaV+uoWsWBe4t7K04cKCgqLhQERIYgYSaQnN8fKbdGVgKJIfh+nofH5N5zbt6bmOTNOeeew2OMMRBCCCGEqJCWugMghBBCSM1HCQchhBBCVI4SDkIIIYSoHCUchBBCCFE5SjgIIYQQonKUcBBCCCFE5SjhIIQQQojKUcJBCCGEEJWjhIMQQjREYWEhPvvsM5iZmaFjx4549eqVukMiRG6UcBBCiIb4559/kJCQgPj4eLRt2xbffPONukMiRG6UcBBCiIaIiorCwIEDYW5ujtGjR+PWrVvqDokQuVHCQchHZtq0aTA3N1d3GCpx5swZNG/eHAYGBuDxeIiLiwMALFmyBG5ubtDR0YGHh4d6g6yCZs2aYf/+/Xj79i22bt2Kpk2bqjskQuRGCQchRO0mTpyIWrVqVekYWVlZGDp0KNq2bYvMzEwwxlC7dm2cO3cOs2bNwrJlyyAUCvHo0SO1xHvixAl07twZNjY2MDIygoeHB37++WdkZ2fLfYwpU6bg7du3sLCwwB9//IE1a9YoGjohakMJByEfmVWrVuHt27fqDkPpIiMj8fbtWwwZMgT6+vrc9gsXLsDQ0BABAQHQ1tZWS2z79+9Hr1694OTkhDt37iA1NRVBQUFYsmQJ+vXrJ/dxRCIRnj17Bg8PDxQWFuLJkycqjJoQ5aKEgxBSI6SmpgIADAwMSmx/f9uHFhISAh0dHfzzzz9wcnKCkZERBgwYgHHjxiEsLAwvX76U6zj79u1DVlYWduzYATs7O2zYsEHFkROiPJRwEFKNjRo1Ck5OTkhNTcWgQYNgamqKWrVq4fPPP4dAIJAp++uvv4LH43F/ZmZm6Ny5M86fPy9TrrQxHMWPk56ejhEjRsDCwgK+vr7lxhYZGYlBgwbBxsYGhoaG8PX1xa5du7j9HTp0gJ+fX4l6xXHm5OQAALp06YINGzbgzZs3MvGnp6dzdU6cOIF27drB2NgYRkZGaNu2LY4dO8btb9GiBYYNGwYA8PX1lTnO+vXrZY79119/lXlO8jyH8sT7PjMzszL3aWlpwcTEpMz97woODkabNm3QrFkzfPbZZ9i3b59CXTKEqBUjhFRbI0eOZPb29iwgIICdOXOGZWVlsVOnTjFbW1vWrl07JhaLS61XVFTEnj9/zj7//HOmr6/PHj58yO378ssvmZmZWamPM2DAAHbq1CmWkZHBNm/eXGZc586dY7q6uqxv377s3r17LCcnh92+fZsNHz6cJSYmMsYYa9++PWvVqlWJugsWLGAAWHZ2NrdtwoQJzMrKqtTH2rVrF+PxeGzq1KksISGBJSYmssDAQMbj8di2bdu4cjt37mQA2K1bt2Tql3fs8pT3HCp6zHv37jErKys2ZswYlpCQwHJyctiBAweYoaEhmzNnjlzHePr0KQPAtmzZwhhj7MWLF0xLS4utW7dOsRMjRE0o4SCkGhs5ciQDwA4ePCizfdeuXQwAO3DgQIXHsLa2ZrNnz+bul5VwAGD79++XK6769eszT09PVlRUVGYZZSQcEomEOTk5sebNm5fY5+fnx2xtbbkYlJ1wvOv957Ayx7x69SpzdHRkALi/r7/+usyk8X2zZ89mFhYWLD8/n9vWo0ePUp9jQqoj6lIhpJrT0tJCnz59ZLb1798fPB5PpqlfIBBg1qxZcHd3h76+PtfUn5aWhmfPnlX4ODwer8TjlObp06d4+vQphg0bpvJBmI8fP0ZiYiIGDhxYYt+gQYOQkpKC+/fvK+3xqvocluXYsWPo3Lkzunfvjvj4eAgEAhw8eBBbt27FoEGDKqwvFosREhKCsWPHygyInTx5Mm7cuIEHDx5UOjZCPhRKOAip5iwtLaGjoyOzTV9fH6ampjLjBgYOHIgtW7bgr7/+QnJyMiQSCXdpaGFhYYWPY21tDV1d3QrLFQ/OdHR0VPBMpBhjcpd98+YNAMDOzq7EvuJt5Y2dUFRVn8OyzJgxA46Ojli3bh1cXFxgYmKCgIAAzJs3D6GhoTh58mS59Y8dO4bk5GT89ddfMuNGBgwYAAA0eJRoBEo4CKnmMjIyUFRUJLNNKBRCIBDAysoKAPDq1SucO3cO3377LXr16gULCwvweDwUFRUhMTFRrsfh8/lylbO2tgYAJCUllVvOzMys1AGNFdV7l6WlJQAgJSWlxL7ibVWdv6OYMp7D0kgkEjx9+hReXl4lWoS8vb0BoMK5QTZs2IA+ffqASbvBZf7+/PNPbN26FSKRqNIxEvIhUMJBSDUnkUhkrsgAgEOHDoExhs6dO8ts19PTk7m/Y8eOEslKVbm7u6N+/frYvXs3xGJxmeXq1q2LuLg47moUQLr42OnTp0uUNTIyglAoLLG9QYMGcHR0RGhoaIl9Bw4cgK2tLRo2bFi5EymDPM9hWfGWRktLC05OToiJiYFEIpHZFx0dDQBwdXUts/7r169x/Phx9OjRo9T9PXv2RHp6Og4fPixXPISoCyUchFRz9vb2CAkJwblz55CdnY0zZ85g+vTpaN26NTdplIODA3x8fLB8+XJERkZCIBBg37592LRpE+rWrav0mNauXYvY2FgMHDgQ0dHRyM3NxZ07dzBy5EiuBWPChAnIz8/HzJkz8ebNG7x48QITJkxAs2bNShzP29sbOTk5uHDhgkwSo6WlhUWLFiEiIgLTpk1DUlISXr16ha+//hrXrl3DokWLSnQ3VZYiz2FZ8Zblu+++w7NnzzBlyhQkJCQgJycHoaGh+OWXX9CoUSP07t27zLohISEoKioqM+Fo0KAB6tSpQ90qpPpT23BVQkiFRo4cyRwdHdnr169Z//79mbGxMbO0tGQTJkxgmZmZMmXj4+NZQEAAMzc3Z+bm5mzo0KEsJSWFNWjQgPXv358rV9ZVKo6OjgrFFhERwfr3788sLS2ZoaEh8/X1ZTt37pQps337dla/fn2mq6vLmjRpwk6fPl3qVSpCoZCNHj2aWVpaMh6PxwCwtLQ0bv+RI0dY69atmaGhITMwMGD+/v7s0KFDMo+ljKtU5H0OK4q3NIcPH2bt27dnNjY2zNDQkHl6erI5c+aUeB3f5+7uzurVq1dumS+++IJpaWmxly9fynWehKgDjzEFRnARQj6oUaNGISwsrEpjCAghpDqgLhVCCCGEqBwlHIQQQghROUo4CCGEEKJyNIaDEEIIISpHLRyEEEIIUTlKOAghhBCicsqZMaeGk0gkePXqFUxMTMDj8dQdDiGEEFItMMaQnZ0NBwcHaGmV34ZBCYccXr16BWdnZ3WHQQghhFRLCQkJcHJyKrcMJRxyMDExASB9Qk1NTdUcDSEfN6FQiGXLlmHGjBkl1j0hhHxYAoEAzs7O3PdkeSjhkENxN4qpqSklHISoWX5+PtLS0mBsbAwDAwN1h0MIAeQabkAJByFEoxgYGCA4OFjdYRBCFERXqRBCNEpBQQFmzJiBgoICdYdCCFGAxrVwFBYWIiwsDCkpKWjUqBGaNGmikjqEEEIIUR6Nmmn0zZs36Ny5M3JyctCoUSOcP38eo0ePxqpVq5Ra530CgQBmZmbIysqiMRyEEELIvxT5ftSoFo45c+agsLAQd+/ehZGREW7duoVWrVqhd+/e6Nmzp9LqEEKqr/z8fAQGBmLlypU0aJQQDaIxYzgkEgl2796Nzz77DEZGRgAAX19f+Pn5YefOnUqrQwip3rS0tODk5FThJEOEkOpFY1o4EhISIBAI4OXlJbO9YcOGuH37ttLqANLr/IVCIXdfIBBUIXJCajDGAKEAyH8LFLwFhNmAKA8Q5QCiXOlfUT5QJAKKCgCxCCgSSv9lYkAiBiRF//0LJj0mY//elsg8nIQxxKXkIMBAjJi/+qnhhAnRfIZd58C9eYcP/rgak3AUf+lbWFjIbLe0tCwzIahMHQAICgrCvHnzqhIuIZqtSAhkJQJv46X/5qQAOan//psm/Tc/U5pkvJcUqJIWAOdChrGh+QgJMIAhn5YaIERRUZmv1fK4GpNwFPfVZmdny2zPzs4usx+3MnUA6biPGTNmcPeLZ1IjpEYpEgIZL4D0J//+PQUyYoG3CUCOgh9IOvqAvjmgbwrwDQFdY0DXCNA1lN7X1gV09P79Vx/Q5gNaOv/+aUv/5Wn9+8cDwJO9/a/7r7Kw4XIs8l2jEdm4J3T5GvMRRki14VzfRy2PqzHvVhcXF/D5fMTFxclsf/HiBerVq6e0OgCgp6dHUyaTmiUvA3h9D0i+9++/d4E3sdJujbLwjQBzZ8DMGTCxBYyL/2wAIxvA0AowMJcmGnz9D3Iaj1kiDvLuot2gT9FmWKsP8piEEOXQmIRDV1cX3bt3x65duzBx4kTweDwkJyfjwoULWL16NVfuypUrSE1NxcCBA+WuQ0iNIpEAaY+AhOvAy3//3saXXlbXBKhVH6jlLv3Xqh5g4QqYuQCGlv+2MFQf2QWFkIgKcHH518j99Cw3GJwQUv1p1DwcDx8+ROvWrdGpUyf4+fkhJCQE5ubmCAsLA5/PBwBMnDgR169fx/379+WuUxGah4NUa4xJu0NizwPPLwAvw4GCrJLlLGoDdo0B+8aAXRPAzhswsa92SUV5Vp57iiUnH8ArJwqHln8PXV1ddYdEyEetxs7D4eXlhXv37mHz5s2Ij49HYGAgxo8fL5M4tGvXDo6OjgrVIUTjCLOBZ+eA2HNA7AUgK0F2P98QcGoBuPgDzq0ARx9p94eGyxYWgafNR9s+QynZIETDaFQLh7pQCwepFvIygMfHgZgj0iRD/N+l29DWA1z9gbqdgNrtALtG0oGZNcycA/ew/cpTaJ/8FfdvXaEuFULUrMa2cBDy0SnIAh6EAvf3AXFXZQd5WtYF3HtIkwzX1tIrQmo4QUEReNo66Dl8IrVwEKJhKOEgpLoRF0nHYtzdCTw6Jp0wq5htI8CzL+DVD7D20KjxF8qQ/W/C0alXf+oWJUTDUMJBSHXx9iUQsRGI2iGdWKtYrQZAk0+BhgGApZvawqsOpFep5OPbYZ3RI+o2jI2N1R0SIUROlHAQok4SCfAiDLi5Hnhy8r9ZOw2tAO/BQNPhgH3Tj64loyzZBUXg6eji6x8WQF//w8z9QQhRDko4CFEHUR4QuQ24+Q/w5tl/2906AL4TgfrdAR0ao/C+7IJC8LS00bVbd+jo0McXIZqE3rGEfEj5b4Fb64Hra4G8dOk2XROg6QhpomHtrtbwqrvsgiJIhHno6uuFJ48fwcTERN0hEULkRAkHIR9C9mvg+mrg1kZA9O/aPuauQOtA6fgMPfrirEiRWII8kRg8vh42b9tZ7npIhJDqhxIOQlQpLwO4vBS4Ffzf1SY2DYG23wANBwDa9BaUV46wCADA09JGx0/aQkdHS80REUIUQZ92hKiCMBsIXw1cW/lfi4ZTS6Ddt4B7dxoEWgnZBdKEQ0dcgFqW5khMTKSJ+AjRIJRwEKJMRSJpa8blpf+N0bBrBHT+GajXhRKNKhAUFAIAzEyMcSg8nGYZJUTDUMJBiLI8OQWcnANkxErvW9UDOv4AeAUAWtT8X1XFLRxmhnpo2LChmqMhhCiKEg5Cqir9qTTReHZGet/IBuj0A9B0FI3RUKLihEOfCcHj8WhtI0I0DH0aElJZwmzg4h/SS1wlhYAWH/CbCnwyC9CnL0Jly/63S8XCzBQJCQk0yyghGoYSDkIq4/FJ4NgMQJAkvV+/G9A9CKhVT71x1WDFLRwmBnyYmpqCR+NhCNEolHAQooicNODkd8D9/dL7FrWBnoukV54QlSpu4dBjQrmXwyaEVB+UcBAiD8akq7ee+h7IzwR4WoD/NKDDnI9iWfjqoLiFw8rCHFlZWTTLKCEahhIOQiqS/Ro4NO2/QaF2jYB+KwGHZuqN6yMj+DfhMNbThkAggLGxMXWrEKJBKOEgpDwPDwNHpgP5GYCOvrRFw/9LQJuv7sg+OsVdKrpMBGdnZ+pSIUTDUMJBSGkKBMCJ74C7O6T37RoDA9cDNh7qjesjVtylYmNlAcaYmqMhhCiKZiMi5H3x4cDaNtJkg6cFtJ0BTDxHyYaaFbdwGPF5ePDgAcRisZojIoQoglo4CCkmkQBX/wLOLwCYBDB3AQasA1z91R0ZwX8tHHxJIfz9/WktFUI0DCUchADSVV0PTgGenpLebzwM6LWEJvCqRooTDjtrCwgEAjVHQwhRFHWpEJIYAfzziTTZ0NYD+q4ABvxDyUY1U9ylYqDDQ3h4OIqKitQcESFEEZRwkI8XY8CNf4CNPYCsBMDSDZh4FvAZS6u6VjNiCUOuSDpmg88KMWTIEOTn56s5KkKIIqhLhXycCguAo9/8dxWKZz+g/ypA30y9cZFS5RT815phb22BxMRENUZDCKkMauEgH5/s18Dm3v9ehaINdF8IDN1CyUY1li38dw4OHS1og+HUqVPUpUKIhqGEg3xckm4D6zoASRGAvjkwar90Ii/qQqnWigeMmurroKCgADNmzEBBQYGaoyKEKIK6VMjH494e6RTlYiFg7QF8ugOwqqvuqIgcuJVi9fkwNjbGgwcP1BwRIURR1MJBaj7GgAsLgQOTpMmGe09gwhlKNjRI8RUqJvo6KCwsxN69e1FYWKjmqAghiqCEg9Rs4kLg0JfAxT+k99t+I23ZoEteNcp/LRw6EIlEWLZsGUQikZqjIoQogrpUSM1VIAD2jAGeX5AODu2zDPAZp+6oSCVwLRx6fBgZGSE8PFzNERFCFEUtHKRmErwCNvWUJht8I2D4Lko2NJjgvRaO4OBgauEgRMNQwkFqntRHQHAXIOU+YGQDjD8GuHdTd1SkCt4dNEpjOAjRTNSlQmqWpNvAtkFAfiZgVR8YtQ+wqK3uqEgVvTto1MjICKdOnVJzRIQQRVELB6k5XlwGQvpJkw1HH2DCaUo2aoh3B40KhUIsW7YMQqFQzVERQhRBCQepGR6flLZsiHKAOp8AYw4BhpbqjoooSXELh6k+H2KxGOHh4RCLxWqOihCiCOpSIZoveh9wcDIgKQIa9AIGbwL4+uqOiijRuy0choaG2Lt3r5ojIoQoilo4iGa7HQLsnyhNNhoNla6JQslGjfPuoFGhUIhffvmFulQI0TAamXA8ffoUV65cQXp6utx14uPjcfbsWWRmZqowMvJBRWwCjnwFgAEtJgAD/gG0+eqOiqjAu4NGJRIJEhMTIZFI1BwVIUQRGpVw5OXloVevXmjRogUCAwPh7OyMJUuWlFvn+vXr6NOnD/z8/NC1a1fcvXv3A0VLVCpiE3D0a+ltvy+B3ksBLY3670wU8G6XioGBAYKDg2FgYKDmqAghitCoT+iff/4ZDx48wNOnTxEZGYn9+/dj1qxZuHbtWpl1YmNjMXnyZNy4ceMDRkpUKmKjbLLR/Tda7bUGk0gYckT/danQarGEaCaNSjhCQkIwceJE2NjYAAB69eqFpk2bYvPmzWXWGTlyJPr27Qst+vVbM0RsBI5+I73tP42SjY9AjqgIjElvm+jTOHdCNJXGvHuTkpKQlpaG5s2by2xv3rw5oqKilPpYQqFQZkCaQCBQ6vFJJUVskk02uv1KycZHoLg7RVdbC/p8bYCvjWXLlqk5KkKIotSacDx8+BCvXr0qt0ybNm1gYGDADfa0tJSdW8HKykrpA0GDgoIwb948pR6TVNHdXf91o1Cy8VF5d8AoAOTn5yMwMBArV66kcRyEaBC1JhzHjh2rcIriLVu2wMDAALq6ugCkHzbvysvL4/Ypy5w5czBjxgzuvkAggLOzs1Ifgyjg4WEgdKr0dsvJlGx8ZN4dMAoAWlpacHJyom5SQjSMWhOOWbNmYdasWXKVdXZ2hpaWFhITE2W2JyYmonbt2kqNS09PD3p6eko9JqmkZ2eBfZ8BTAI0HQX0+J2SjY/Mfy0c0kue9fT08Msvv6gxIkJIZWjMTwQDAwO0a9cOoaGh3Lbs7GycO3cO3br9txJoTEwMrl+/roYIidLFXwN2jQIkhYBXANBvBV36+hF6v4UjLy8PQ4YMQV5enjrDIoQoSGMGjQLAwoUL0bFjR3zzzTfw9/fH33//DUdHR0ycOJErs3TpUly/fh33798HACQnJ+PBgwfcJGG3b99GUVER3Nzc4ObmppbzIHJIugNsHwoU5QP1uwED1wNa2uqOiqiB4L2EQ1tbG/7+/tDWpv8PhGgSuRMODw8PhQ786NEjhYOpSOvWrXH16lWsXr0aISEh8Pf3x8yZM2FkZMSV8fLyAu+dJvfo6GgsWrQIANC5c2ecOHECJ06cwJgxYyjhqK7Sn/67EFs2ULuddLpyHeWO0yGao7QulXfHWBFCNAOPseIr3CsoyONh5cqVch00MDAQch5WIwgEApiZmSErKwumpqbqDqdmEyQDG7oBWS8Bh+bA2MOAnom6oyJq9MfJR1gTFovxbWrj574NkZubi4EDB+LAgQMyPzYIIR+eIt+PCnWpTJs2Ta5ygYGBihyWEKmCLGD7EGmyYVkXGLmXkg1SooWDz+djyJAh4PNp3RxCNIncCYcic13QAmlEYUVCYNdIICUaMLYFRh8AjGqpOypSDRQPGjX9dwyHrq6uzLgtQohmkHvIv7m5udwHVaQsIZBIgIOTgbjLgK6JtGXDora6oyLVxPtXqeTm5sLf3x+5ubnqDIsQoiCFrjFs0KABFi1ahJSUFFXFQz42jAGnvgceHAS0+MCn2wD7JuqOilQjxV0qxnrSLhRdXV3MmDFD6RP+EUJUS6GEo0uXLggKCoKTkxMGDBiAY8eOQSwWqyo28jG4vga4sUZ6e8BawK2DWsMh1c/7LRw0hoMQzaRQwvH3338jOTkZmzdvRlZWFvr27QtXV1f8+OOPePHihapiJDXVo+PS1g0A6LoAaDRYvfGQaun9hCMnJwcNGzZETk6OOsMihChI4Wkb9fX1MXLkSJw/fx7Pnj3D+PHjERISgrp166Jz587YuXOnKuIkNc2rKGD/BAAM8BkPtKYrm0jpBO9dpaKvr49ly5ZBX19fnWERQhQk9zwc5ZFIJNi3bx+mTJmCzMzMGjUHB0DzcChdVhIQ3BnITgbcOkoHiWpT8zgpSSJhqPvDcTAG3Py+M2xMKckgpDpR5PuxygtTXLlyBRMnTsRnn30GkUiE8ePHV/WQpCYT5gA7h0mTDWsPYGgIJRukTLmiIhT/filu4cjOzoaTkxOys7PVGBkhRFGVSjiSk5Pxxx9/oEGDBmjXrh0ePHiAP//8E8nJydi4caOyYyQ1hUQs7UZ5HQ0YWQMj9gD6ZuqOilRjxeM3dLR40OdLP64MDAywd+9eGBgYqDM0QoiCFJpp9ODBg9i4cSNOnDgBc3NzjBo1Cvv374e3t7eq4iM1yZmfgCcnAR194NOdgIWruiMi1dy7A0aL10jS0dGBv7+/OsMihFSCQi0cgwYNQkFBAbZv346kpCT89ddflGwQ+UTtBMJXSW8HrAacfdUbD9EI709rDkj7jE1NTSEQCNQVFiGkEhRq4Xj+/Dlq166tolBIjZUYARyZLr3dbibgPUi98RCN8f4lsQBgZGSE8PBwWriNEA2jUAvHu8lGbm4uDh06hD///JPb9uTJkxp3hQqpIkGydI0UsRBo0Bvo+IO6IyIa5L9LYv9LOLS1tdGwYUNoa2urKyxCSCVUatDokydP0LBhQ0yePBkzZszgti9YsAC7d+9WWnBEwxUWALtHAjmvAWtPYOA/gFaVL4wiH5H/Wjhku1R4PB51qRCiYSr16f/1119j6NChSE5Oltk+ffp0LFmyRCmBEQ3HmLQbJek2YGABDN9BS80ThZXWpWJsbIyEhAQYGxurKyxCSCUoNIajWHh4OHbs2MGNGi/m4eGB6OhopQRGNNz11cC9XQBPGxiyGbB0U3dERAMVDxo1faeFg8fjwdTUtMTnDyGkeqtUC4dEIoFIJAIAmTd9XFwczMxoXoWP3ovLwOm50tvdf6MF2Uil5QhLtnBkZ2fDzMyMJv4iRMNUKuHo2rUr13VSnHCkp6fjq6++Qo8ePZQXHdE8glfAvvEAEwONhwGtpqg7IqLBSutSMTExQVZWFkxMqIuOEE1SqYRj6dKl2Lt3Lzw9PcEYQ8eOHVGnTh3ExcXh999/V3aMRFMUiYA9Y4HcNMDWG+jzF0DN3qQKSpuHgzEGgUBAV8QRomEqNYbD1dUV9+7dw7Zt2xAREQGJRIKBAwdi7NixtLjZx+zU90DiTUDPDBi2FdA1VHdERMMJSmnhyMnJgbOzMy2mSIiGqVTCAUibNadOnarMWIgmu7sLuLVeenvgOhokSpSitMtiTU1NqXWDEA1U5UkRCgoKSvyRj8zraODI19Lb7b8DGtA4HqIc2aVM/CUWi/HgwQOIxWJ1hUUIqYRKJRxPnz5F165dYWhoCAMDgxJ/5CNSkAXsHg0U5QP1ukgTDkKUpLiFw/SdhCM3Nxf+/v7Izc1VV1iEkEqoVJfK2LFjYWVlhT179sDc3FzJIRGNwRhwOBDIfAGYuQAD1wNaNN00UQ7G2DuXxcp2qdAso4RonkolHJGRkUhOTqZk42N3cz3w8BCgxQeGbAIMLdUdEalB8kRiiCXSsRrvdqkUFRXh1q1b8PX1hY5OpYehEUI+sEp1qbi6uiI9PV3ZsRBNknRHelUKAHSdDzi1UG88pMYp7k7R1uLBgP9fy1l+fj6GDBmC/Px8dYVGCKmESiUc8+bNw/jx43H9+nWkpaUhPT1d5o/UcPlvgb3jAEkh4NEH8KOrlYjyvTtg9N0ZjU1MTJCYmEgTfxGiYSrVHmlra4vo6Gj4+/uXup8uWavBGAMOfQm8jQfMXYH+f9PkXkQlSpuDA5B2qZw7dw6dO3emLhVCNEil3q1Tp05Ft27d8MUXX9A4jo/NjbXAo6OAtq50UTYDc3VHRGooroVDjy+zvaCgADNmzMCNGzdoxVhCNEilEo74+Hhcv36dFmr72LyK/G9Rtm6/AY7N1RsPqdFKW0cFkC5P/+DBA3WERAipgkqN4ahXrx6Sk5OVHQupzoTZwL7PpOM2PPsCLSepOyJSQ2XkirD9RjzWXXoOQPaSWAAoLCzE3r17UVhYqI7wCCGVVOl5OEaNGoXFixejXr16MgO6AMDJyUkpwZFq5PgsIOM5YOoE9FtJ4zaIUmXlFeLkg2QcvZeMa7FvuMthAeAT91oyZUUiEZYtW4ZevXqBz+e/fyhCSDXFY5UY4fl+gvG+mjZoVCAQwMzM7ONdLOrubuDg5wBPCxh3HHAtfbAwIYrIFRbhbEwKDke9wqWnaSgU//e50cjRDL0b26N3I3s4W9IigIRUV4p8P1aqhSMmJqZSgREN9CYWODZDerv9bEo2SJUIi8S4+DgNh+6+wrmYFBQUSrh9HnYm6NvEAb0b2aN2LaMyjyESibBlyxaMGTMGurq6HyJsQogSVCrh8PDwUHYcpDoqEgH7JwCiHMC1DfDJTHVHRDSQRMJwMy4Dh6KScDz6NbLy/xt7UdvKEP2aOKBvEwfUt5VvXo3iMRzDhw+nhIMQDSJ3whEQEIDQ0FCllyXV2Pn50itT9M2lS87TOilEAY9eC3AwMgmHo14hOeu/VaRtTfXQt7ED+jd1hLejaYVdtO8zMjLCqVOnlB0uIUTF5B7DwePxkJ2dLddBTUxMatQ4jo9yDEfseWDrAOntYdsBzz7qjYdohFRBAQ5FvcKByCTEJP+3wJqJvg56edujf1MHtHKzgrZW5QcdC4VC/P333/jyyy+hp6enjLAJIZWksjEc1WUq4dzcXGRmZsLe3h7a2vL96s7MzISOjk61OYdqLS8DCP1CervFZ5RskHLli8Q49eA19t9JxNVn6Si+wISvzUPHBjYY0MwRHT1soM9XTguZWCxGeHg4pkyZopTjEUI+DLkTjvDwcFXGIRexWIyvvvoKwcHBMDIygra2NlauXIlPP/20zDobNmzAkiVLkJqaCpFIBDc3N6xevRpt2rT5gJFrEMaAI9OB7GTAqr50gi9C3sMYw80XGdh/JxHHo19zy8gDgI+rBQY0c0TvRvawMFL+GAtDQ0Ps3btX6cclhKiW3AmHn5+fKuOQy6JFi7Bnzx5ERUXB09MT69atw6hRo+Dt7Q1vb+8S5Yt/CYWGhqJBgwYoLCzEN998g759++Lp06ewsrJSw1lUc1E7gJjDgJYOMGg9oEuXJJL/JGbmYf/tJOy/k4iXGXncdmdLAwxs5oQBzRzLvcJEGYRCIYKCgjBnzhzqUiFEg1RqHg51cXFxwciRIxEUFMRtc3d3R48ePbBixQq5jpGcnAwHBwecOHECPXr0kKvORzOGI+MFsLat9KqUzj8B7b5Vd0SkGigolHaZ7IlIwLXYNyj+xDDW00GvRnYY1NwJvrUtoVWFcRmKyM/PR2BgIFauXAkDA4MP8piEkNKpfB4OdUhJSUFCQkKJFWrbtGmDiIgIuY/z6NEjAICjo6NS49N44iLg4GRpsuHSGmjztbojImrEGMP9JAF2R7zEoahX3LomANC6rhWGtHBCj4b2MND98FcuGRgYIDg4+IM/LiGkatSacKSmpkIgEJRbxtXVFXw+H+np6QBQohukVq1auHr1qlyPl5OTg8DAQHTr1g2NGjUqs5xQKIRQKOTuVxRjjXBlGZBwA9AzBQaspUtgP1JZeYUIjUrCrlsJMleZOJobYLCPEwb7OKl95s+CggJ8//33WLhwIfT19dUaCyFEfmpNOJYuXYr9+/eXW+b8+fNwcXGBlpZ0nbn3F2wSiURyXalSUFCAgIAASCQSbN++vdyyQUFBmDdvXoXHrDGSbgNhv0tv91oCWLiqNx7yQRUPAN158yWO338NUZF09k9dbS1097bDsBbOaF3X6oN1mRBCaqZKj+FIS0tDcHAwN825l5cXJk6ciFq1alVQs3Kys7NhamqKnTt3ylyVMmzYMGRkZODMmTNl1hUKhejfvz/i4uIQFhYGOzu7ch+rtBYOZ2fnmjmGozAfWNsOePMUaDgQGLyRFmb7SGTkirD/diJ23nqJ52m53HYPOxMM83VGQFNHlVxlQgipORQZw1Gp5enDwsLg5uaG9evXIz8/H/n5+Vi3bh3c3Nxw6dKlSgVdERMTE/j4+MjMMCgSiXD27Fl06NCB25aWloaEhATuvlAoREBAAOLi4nDhwoUKkw0A0NPTg6mpqcxfjXV2njTZMLYDei+lZKOGY4zhxvM3+GpnJPwWnsNvx2PwPC0Xhrra+NTXGaFftsGJ6e0wvk2dapts5OfnY+LEicjPz1d3KIQQBVSqSyUwMBBffPEFgoKCuK4OiUSCOXPmYNq0abh3755Sgyz2888/Y+DAgWjatCn8/f2xbNky6OrqykwANGfOHFy/fh3379+HWCzGoEGDcOfOHezfvx+5ubl49uwZAMDGxqZmJxLyeHEJuLFGerv/KsDQUr3xEJXJyi/EgTuJ2H7jJZ6l5nDbGzmaYXhLF/Rr6gBjPc0YQ66lpQUnJyfus4cQohkq1aViaGiIpKQkWFhYyGzPyMiAk5MT8vLyyqhZdaGhoVixYgVSUlLQqFEjLFiwAPXr1+f2f//994iKisLx48eRmZkJX1/fUo8TFBSEIUOGyPWYNfKy2IIsYE0bICsB8BkP9P1L3RERFbiflIUt4XE4fPcVtzKrAV8bAc0cMKKlKxo5mak5QkKIJlP5ZbHu7u549uxZiS/z2NhYuLu7V+aQcgsICEBAQECZ+xcuXMjdtrCw4Fo0yHtOfi9NNixqA91+VXc0RIkKCsU4Hp2MLeHxiEp4y21vYGuCUX4u6N/MEab6fPUFWEV5eXkYO3YsQkJCYGhIE9MRoikqlXBMmTIFQ4YMwfz58+Hr6wvGGCIiIvDTTz9hzpw5SExM5Mo6OTkpLViiJI+OA1HbAPCAgDWAnrG6IyJKkPQ2H9uux2P3rQRk5IoASNcz6eltj9H+rmjhaqHwyqzVkba2Nvz9/eVeR4kQUj1UqktFkQ8tDZrItEw1qkslNx1Y7QfkpgGtvwK6LVB3RKQKGGO4/jwDm6+9wJmHKdzCaQ5m+hjRygXDfF1gbULTfxNCVEPlXSrFl8ISDXRshjTZsPYEOv6g7mhIJeWLxDgYmYSQa3F4nJLNbW9Tzwpj/Gujs4cNdLRr5qDK3NxcDBw4EAcOHICRkWrXbSGEKE+lEg4PD48S2+Li4mBjY0N9qtXZ/QPAw0MAT1s6myifZmnUNMlZ+dgSHo+dN1/ibZ50EjwDvjYGNnfE2Na14W5rouYIVY/P52PIkCHg8zV3HAohH6NKJRy3b9/G5s2bsXLlSgDAqFGjsH37dpiYmODkyZNo3bq1UoMkSpCTBhyfKb3d7lvAoalawyGKufMyE5uuxuF4dDLE//abOFsaYKx/bQxp4Qwzg4/ny1dXVxcTJ05UdxiEEAVVKuGYOXMmFiyQ9v3fv38fR44cwfXr13H69GnMmTMHFy9eVGqQRAmOfwvkvQFsvYFPZqk7GiKHIrEEpx+mYP3l54h8+Zbb7udmifFt6qCLpy20P8LpxnNzc9GlSxecPXuWulQI0SCVSjgiIiLg4+MDADh9+jQGDBiAVq1aoWHDhliyZIlSAyRKUNyVoqUDBKwGdKrnDJJEKkdYhN23ErDp6gskZkpn09TV1kK/pg4Y36Y2Gjp83HNn6OrqYsaMGdDVpf/HhGiSSiUcxsbGiIuLg6enJ44ePYpx48YBALKysmBsTJdYVivvd6XYN1FvPKRMyVn52Hw1DjtuvES2ULocvIUhH6P9XDHK3xU2JjTmBvhvDAchRLNUKuEYPHgwevfuDW9vb0RFRaFv374AgFOnTqFnz55KDZBUAWPSq1Ly3gC2jYB2M9UdESnF49fZWHfpOQ7fTUKhWDo+w83aCBPa1sHAZk4w0KX5Jt6Vk5ODVq1a4caNG/QDhxANUqmEY9myZahXrx7i4+Mxb948borzR48e4eeff1ZqgKQKHhwAYg5TV0o1VDx/xj+XYhH2OI3b3rKOJSZ/4oaODWxoOfgy6OvrY9myZdDXpxYfQjRJpZen/5ho5MRfuenA3y2lrRvtZwMd56g7IgJAImE4E5OCNWGx3LTjPB7Qo6EdPv/EDc1cLMo/ACGEVCMqmfjr0aNHAKRzcBTfLktp83SQD+zEd9Jkw6ahdOwGUatCsQSHol5h7cVYbrVWXR0tDPFxwqR2bqhdi662kFd2djY8PT0RExMDE5OaP+8IITWF3AmHp6cnAGlTcPHtslCjiZo9Og7c3wfwtKTLzlNXitoUFIqx+1YC/rkYi1dZBQAAEz0djPZ3xfg2dWja8UowMDDA3r17YWBgoO5QCCEKkDvhSEhIKPU2qWby3wJHv5Hebh0IODZXazgfq+yCQmy7/hIbrjxHeo50ITVrEz1MaFsHI1q5aPRqreqmo6MDf39/dYdBCFGQ3AnHu6u+0gqw1djpH4Gc14BlXaADjdv40N7mibDpahw2XX0BQYH00lZHcwNM6VAXQ3ycoM+nK06qSiAQwMnJCYmJiZozpooQIn/CERUVJfdBmzZtWolQSJXFXgAit0pv918F8KnJ+UPJyBVhw5XnCLkWj5x/59BwszbCFx3qoX9TB/Br6EJq6mBkZITw8HCaZZQQDSN3wtGsWTO5D0pjONRAmAMc+Up623cS4Err2XwIadlCBF9+jq3X45EnEgMAPOxMMK1TPfT0tv8opx5XNW1tbTRs2FDdYRBCFCR3wpGZmcnd3rdvH37//Xf8/vvv8PX1BQDcunULs2fPxpw51IyvFucXAG9fAmYuQBeaC0XV0rKFWHcpFluvx6OgUAIA8HY0xVed6qOLpy3NoaFCGnmZOiFE/oTD3Nycu718+XLs27dPpuvE1dUV9erVw5gxYzBhwgRlxkgqknATuPGP9HbfvwA9ulRQVdJzhFh36Tm2hscjv1DaotHE2RzTO9dDxwY24PEo0VA1Y2NjJCQk0CyjhGiYSs00+uzZMzg7O5fY7uzsjGfPnlU5KKKAIhFwOBAAA5qMAOp1VndENVJGrgj/XIrFlmvvJBpOZvi6qzs6uFtTovEB8Xg8mJqa0nNOiIapVMLh4eGB+fPnY+nSpdDRkR6iqKgI8+fPp0m/PrQrfwJpjwDDWkD339QdTY2TlV+I4MvPsfHKC+T+O0ajsZMZvunijg4NKNFQh+zsbOpSIUQDVSrhWL16Nfr06cN1qzDGEBUVBZFIhGPHjik7RlKWtMfA5SXS2z3/AAwt1RtPDZIrLMLma3H452Isd3lrQwdTzOjqjk4e1HWiTiYmJsjKyqJZRgnRMJVKOPz9/fH8+XNs2bIFDx8+BAD07NkTY8eOpV8cH4pEAhz+ChCLgPrdAe9B6o6oRigoFGPb9XisCYvFm1zphF3utsaY0dUd3RvaUaJRDTDGIBAIYGxsTK8HIRqkUgkHAJiZmSEwMLDcMj169MDJkycr+xCkPLc3AgnXAV1joPdS6QpgpNKKxBLsv5OI5WefclOQ17YyxDdd3dGnsQNd3lqN5OTkwNnZmbpUCNEwlU445HHq1ClVHv7jlZUEnPlFervzz4B5yQG8RD6MMZy4/xpLTj/G87RcAIC9mT6+7lIfg5o7QYcm7Kp2TE1Naa4fQjSQShMOogKMAcdnAqJswKkl4EuXIFfWtWfp+P3kI9xLzAIAWBrp4osOdTHKz5WmIK/GxGIxHj16BA8PD2hr0+tEiKaghEPTPDwEPD4OaPGBfisALfrAVdTDVwL8fvIRLj1JAwAY6WpjYjs3TGxXBya0qFq1l5ubC39/f1pLhRANQwmHJsl/C5z4n/R2228AG0+1hqNpEjPzsOz0ExyMSgJjAF+bh5GtXDGtUz3UMqZl4jWFqakpBAKBusMghCiIEg5Ncm4ekJMCWNUD2n2r7mg0RlZ+IVZfeIZN1+IgKpJOQ963iQNmdnOHqxUtAKZpioqKcOvWLfj6+nLzABFCqj96t2qKl9eBiI3S233+Avj6ag1HE4iKJNh+Ix4rzj1FZl4hAKB1XSvM7umBxk7m6g2OVFp+fj6GDBmCmJgYmouDEA2i0oRj+vTpqjz8x6NIBBz597lsNgqo00698VRzjDGcevAav594hLg3eQCA+jbG+L63J01DXgOYmJggMTFR3WEQQhSk0DV/a9eulbkfFRVVoszs2bO523/99VelgiLvubb8v+nLuy5QdzTVWnRiFob9cx1Ttt1B3Js81DLWw8IBjXBiejtaXK2GKCoqwqlTp1BUVKTuUAghCuAxBS5o5/F4Mte/v3+/rG2aTq3LYb+JBVb7A2IhMHA90Hjoh318DZEiKMDiU4+x/04iGAP0+VqY1M4Nk9vXhbEe9RzWJDk5OWjVqhVu3LhBK8YSomaKfD/SJ3F1xpi0K0UsBOp2AhoNUXdE1U5BoRjrLz3HmouxyPt3cbWApg74Xw8POJgbqDk6ogrGxsZ48OCBusMghCiIEo7q7O4uIO4yoGMA9F5G05e/o3iG0N+OxSDpbT4AoJmLOX7q44VmLhZqjo6oUmFhIUJDQxEQEAA+n+ZNIURTUMJRXeVlAKd/kN7u8B1gWUe98VQjMckCzDvyANefZwAAHMz0MbuXJ/o2tqcxGh8BkUiEZcuWoVevXpRwEKJBFE449u3bV+59oiRn5gJ5bwAbL8B/mrqjqRYyc0VYduYJtt+Ih4QBejpamNy+Lqa2rwsDXZpx9WNhZGSE8PBwdYdBCFGQwoNG5UGDRqso/hqwqaf09menAZdWqn/MakwsYdh16yUWn3qMt//Op9GrkR2+7+UJJwtDNUdHPjSRSIQtW7ZgzJgx0NXVVXc4hHzUVDZoNDs7u0qBETkUiYAjX0tvNx/70ScbUQlv8dOh+9wCaw1sTfBzPy+0rltLzZERdSksLMTevXsxfPhwSjgI0SAKzcNhbGws158qLV68GI6OjuDz+WjWrBkuXrxYbvmnT59izJgxsLW1hZGREVq2bIkjR46oNMYqubYCSH8snXOjyy/qjkZt3uQIMXv/PQxYfRX3ErNgoqeDn/p44dhXbSnZ+MgZGRnh1KlTMDKiaekJ0SQKJRzqtm7dOsybNw8bNmxAeno6evXqhV69eiEuLq7MOqtXr0a/fv0QExOD169fY9CgQRgwYECpk5apXcZz4NJi6e0eQYChpXrjUQOJhGHHjZfotPQidt1KAGPAoOZOODezPT5rWwc62hr1X5aogFAoxLJlyyAUCtUdCiFEAQqN4VC3Bg0aoEePHli+fDkA6VgRFxcXjBw5Er///rtcx5BIJNDT08PatWsxYcIEuep8kDEcjAHbBgGx5wC3DsDo0I/uMtj7SVn4MfQ+ohLeAgA87U2xoH9DtKj98SVepGx5eXkYO3YsQkJCYGhIY3gIUacaOfFXRkYGnjx5gqCgIG4bj8dDhw4dcO3atXLrMsYgFoshEAjwzz//wNTUFF27dlV1yIp5cECabGjrfXRzbmQXFGLZmScIuRYHCQOM9XQwo6s7xvi7UosGKcHQ0BB79+5VdxiEEAWpNeGQSCSQSCTllilefvr169cAAGtra5n9NjY2uHnzZrnHOHPmDHr16gWxWAxzc3Ps2rULLi4uZZYXCoUyzbUCgaDc41dZQRZwco70drsZgFVd1T5eNVE8ede8Iw+QIpA+370b22Nuby/YmdFquKR0QqEQQUFBmDNnDvT09NQdDiFETmr9+Thp0iTo6+uX+xcbG1vuMSQSSYWX63br1g1FRUXIzMzE7Nmz0a9fv3Kv4w8KCoKZmRn35+zsXKnzk9v534CcFMCqHtD2G9U+VjWR9DYfE0Mi8MX2O0gRCOFqZYgtn7XE3yOaU7JByiWRSJCYmFjhjxVCSPWiMWM4MjMzYWlpiX379mHQoEHc9lGjRiExMRFhYWFyH8vX1xfNmjXDunXrSt1fWguHs7OzasZwvIoE1ncCmAQYc0g6fqMGKxJLsPlaHJadeYI8kRh8bR6mtq+LLzrWgz6fJu8ihBBNosgYDo3pILewsICnpycuXLjAbWOM4cKFC2jdujW3TSKRQCwWl3ssoVAILa2yT11PTw+mpqYyfyohEQNHv5EmG42G1Phk435SFgJWX8Wvx2KQJxLDt7YFjn/VDjO6NaBkg8itoKAAM2bMQEFBgbpDIYQoQGMGjQLArFmz8OWXX6Jr167w9/fHH3/8AYFAgKlTp3JlPv/8c1y/fh33799Hbm4uRo0ahe+++w5eXl54+/YtVq1ahZiYGKxdu1aNZ/KviI3SFg49M6Dbb+qORmXyRWL8de4Jgi+/gFjCYKqvg+97eWJoC2doaX08g2OJ6onFYhQWFqo7DEJqDD6fD21t5fwg1KiEY/z48cjOzsY333yDlJQUNGrUCKdPn5YZY6Gtrc0NNDUyMsLkyZPx/fffIzIyEkZGRtxkYe+2iqhFdgpwbr70due5gImteuNRkWux6ZhzIBrxb/IAAH0a2+Pnvg1hbUKD/Ujl6OvrY9myZTLbGGN4/fo13r59q56gCKnBzM3NYWdnV+XFMTVmDIc6qWQejn0TgPv7AIdmwMRzgFbN6lLIyi9E0PEY7LqVAACwM9XHrwHe6OJVMxMr8uHk5+cjMDAQK1euhIGBAQAgOTkZb9++hY2NDQwNDWnVYEKUgDGGvLw8pKamwtzcHPb29iXK1Mh5OGqU2AvSZIOnBfT5s8YlG+diUvD9wWjuUtdRfi74rocHTPRpKXFSdVpaWnBycuLGYYnFYi7ZsLKyUnN0hNQsxUl9amoqbGxsqtS9QgnHhyaRACe+k95u+bm0haOGeJsnwrwjD3EwMgkA4FbLCL8PaoyWdWimUKI8enp6+OWXX7j7xWM2aNZRQlSj+L1VWFhYpYRDY65SqTG0tIBB64EGvYCOP6g7GqU5ef81uiy7hIORSdDiAZM/ccPx6e0o2SBKl5eXhyFDhiAvL09mO3WjEKIaynpvUcKhDvZNgOE7AX0VXW77AWXmijBtxx1M2XYb6TlC1Lcxxv6prTGnlydd6kpUQltbG/7+/kobOV9dBQYG4saNG+oOo1rZs2cPAgMDMW7cOEgkEuTk5GDp0qX4/PPPsXTp0jLrpaenY+rUqXQFUymeP3+OGTNm4EMM56SEg1Ta6Qev0fXPSzh6LxnaWjx82bEujn7VFs1cLNQdGqnB9PT0MGPGjBo/rfnOnTvx4sULdYehdDNmzMCVK1cUrrd9+3Z88cUXcHd3R4cOHcDj8TBp0iTs378fLVu2ROPGjcus+8MPP8DAwAB8/ocbR1bZ8/zQ3NzccPXqVYSEhKj8sWgMB1FYVl4h5h15gAP/jtWob2OMpUOboLGTuXoDIx+F3NxcDBw4EAcOHICRkZG6wyEK2rNnDxo3boy2bdsqVO/ChQvo1asXAgMDZbatWLECQ4cOLbNeUlISNm/ejCdPnlQ65sqo7Hmqw1dffYVffvkFY8eOVWnXJCUcRCFhj1Px3f57SBEIocUDJn3ihm+6uFP3Cflg+Hw+hgwZ8kF/rapKUlISQkJCkJiYiPr162P8+PEwNzfn9jPGsGvXLty4cQPGxsaYMGECateuze1fsGABYmNjwePxYG9vj27duqFDhw4yjxEYGIiRI0fi0aNHuHHjBtq1a4cRI0YgMDAQI0aMwLNnzxAREQFTU1N8/vnnFa4ddffuXezbtw+ZmZlo164dhg0bxu0bN24c/ve//8HLy4vbNmfOHHTp0gWdO3fGTz/9hMzMTAQHByMsLAx8Ph/r168HIF1k8/jx4xAKhdxxi69E+u6773Du3Dnw+XyMGzeOO/abN2+wdu1aHD9+HJ999hk++eSTEvFu3LgRLVu2hKurq9xxFj9vFT0/Zb1+ZZ1nWa8FYwwHDx7EhQsXoKenh/bt26Nv377c42RlZWH69OkApP//69Spg9GjR8vE8vTpU/z222+YN28e9u/fj2fPnsHDwwOTJ0/GmzdvsGHDBiQnJ8Pf3x+jR4+WeY4GDBiAzz77DGFhYejYsWO5r39VUJcKkUueqAg/hkZj3KZbSBEIUaeWEfZOaY05PWmsBvmwdHV1MXHiROjq6qo7lCpJSUlB06ZN8eDBAzRu3BivXr1Cu3btZBal+/bbb3H48GHUq1cP0dHR8PHxQUZGBre/efPm6NChAz755BNoa2tj4MCBWLVqlczj7Ny5EwMGDMCxY8fQtGlT1KtXT2b7nj17ULduXdy8eRMtWrTAq1evyow5ODgYvr6+SE1NRYMGDRAaGorvvvuO2x8SElKi/sGDBxETEwMAaNGiBfT09ODh4cHFDQB//PEHAgICYGRkBBcXF3z77bcYNWoUd4xWrVrBxsYG9vb26NChA/enra0Nb29vdOjQodQ5IgDg7NmzaNOmjcy2iuKU5/kp7/Ur6zxLey0YYxg8eDDmzZuH2rVrw87ODtOnT+cSDED6f774nH19ffHkyRM0bNgQjx8/5sqkpKQgJCQEHTt2xNu3b1GvXj0sXLgQffv2RYcOHSAWi1GnTh18/fXX+OOPP2TO3dDQED4+Pjh9+nSpz6HSMFKhrKwsBoBlZWWpOxS1uBOfwTosvsBcvzvKXL87yn4+dJ/lCYvUHRb5SOXk5DA/Pz+Wk5PDGGMsPz+fPXz4kOXn53NlJBIJyxUWquVPIpHIdR779u1j9vb2MttevXrF1beysmKjRo3i9hUWFjJbW1u2cePGco9Zq1YtmW1WVlYsICCgRFkrKyvWtm1b7vHEYjFr1qwZmzZtWqnHTk1NZfr6+mz16tUy25OSkrjbANiZM2dk9jdo0ICtXLmSu+/o6Mg2bdrE3X/9+jUzMDBg27dv57ZFRUUxHo/HLl68yG0bNmwYmzBhgsyxjYyM2MGDB0uNt5iNjQ1bu3atzDZ54qzo+ano9Xv/PIuP+f5rsWPHDubk5MT9f2aMsdjYWKalpcWePHlS5nmNGjWKTZw4kbt/+fJlBoAdOHCA27Z69WoGgB0/fpzbtmjRIubp6VnieCNGjGD9+/cv9bFKe48VU+T7kbpUSJkKxRKsPPcUf4fFQixhsDfTx+LBTdC2fi11h0Y+Yrq6upgxY0a5LRz5hWJ4/XTqA0b1n4fzu8NQt+KPVk9PT6SlpWHhwoUYPnw46tSpU+JXenHzPgDo6OjAzc0NiYmJ3La8vDzs3LkT9+/fx9u3b5GVlYX09HS8efNGZhK07t27lxrDwIEDuT57LS0tDB48GLt37y617IULFyAUCjF+/HiZ7Q4ODhWea3lu3rwJoVCIIUOGcNuaNGkCLy8vXLp0qdRuEkVkZ2dXeqxPec+PPK9fad5/LY4dOwYdHR1Mnz4djDHuj8/n4+7du6hfvz4AID4+Hjt37kR8fDzy8/MRExPDTcr1rne7ROrWrVvqtqSkpBL1jI2N8fr16wrjrwrqUiGlik3LwaA117Di/DOIJQz9mzrg5PRPKNkgaldTxnB4eXnh1KlTuH37Nnx8fODm5lbi0k59fX2Z+1paWtxq2AUFBfD398e6detgZ2eHtm3bomXLlgCkX7LvendcyLssLS1L3E9LSyu1bGZmJkxNTUvEVFVpaWkwNTUt8XrWqlULqampVT6+lZUVMjMzK1W3vOdHntevNO+/Fm/evIGDgwPatm2Ldu3a4ZNPPkH79u2xdu1aNG3aFABw584deHp64sGDB/D09ESHDh3g5uZW4nUGZP/PFI+BeX9baSuqv337VuUz9VILB5HBGMOuWwmYf+Qh8gvFMDPg49cAb/RtUrVfMYQoS05ODlq1asUNpCyNAV8bD+eX/qte1QwUGNPUqVMndOrUCRKJBCdPnkT//v3RqFEjdOvWrcK64eHhePz4MTIzM7lfumfPnlUo1pcvX8rcj4+PL3PQqJOTE9eCUqtW6T88dHV1IRKJZLa9/2X//lUQLi4uXOuMmZkZtz0uLk5m4GRlFY+zUDROoOLnp7zXT96rPZycnHDnzh2ZwbDvCwkJQdeuXbF161ZuW0REhFKvvLl//z7GjBmjtOOVhlo4CCcjV4TJW29jzoFo5BeK0aaeFU59/QklG6RaKV4ttrxf2jweD4a6Omr5k/eL5saNG3j27BkA6a/ODh06wMDAQO4Vb3k8HsRiMfcrVyQSYfHixXLVLbZp0ybu8dLT07FlyxYMHjy41LKdOnWCvb095s6dy00SVVhYiAsXLnBl6tevj8uXL3P3jxw5UqKVwtLSUmbga+vWreHg4CDTOrBnzx68evUK/fv3V+h8StO7d2+ZGOWNEyj/+ano9Xv/PMsyevRo3LlzB9u2bZPZfuTIEeTk5ACQvtbv/r+Ii4vD9u3bKzy2vFJSUhATE4PevXsr7ZiloRYOAgC49CQNM/feRWq2EHxtHv7X3QMT2taBlhZNF02qFx0dnTLHJGgSbW1t9O/fH+bm5nBxccGNGzfQrFkz9OnTR6767dq1Q5s2bdC8eXO0bt0akZGRsLa2VigGFxcXNGnSBD4+Prh+/Trq1KmDL7/8stSyBgYG2L9/PwYOHIiLFy/C09MT9+7dw/fff8+VmTdvHkaOHImoqCgwxpCeng4bGxuZ4wwaNAgLFy5EREQEjIyMsH79egQHB2PYsGG4cOECTExMcOHCBSxevJi7oqYqRo0ahdmzZ+P69evw8/OTO86Knp+KXr/SzrM0n3zyCVasWIFJkyZhzZo1sLe3R3R0NDw9PbmWrqlTp6JNmzZo2bIlHB0dceXKFdSrV09pM6du374dbdu2hbe3t1KOVxZanl4OKlmevpoQFUmw5PRjrLv0HABQz8YYyz9tioYOZhXUJEQ9srOz4enpiZiYGJiYmKCgoAAvXrxAnTp1lD6+QNWKiopw8+ZNJCYmws3NDS1atOD27dq1C61atUKdOnW4bSdOnIC9vT3Xty+RSBAWFobXr1+jQYMGcHd3x/79+zF48GCuu6m04wDSMRKrVq1Cq1atEBUVBRMTE3Ts2LHCKePz8/Nx6dIlrmvLyclJZv+LFy9w584d2NnZwdfXF8eOHYOnpyc8PDy4MhEREXj69CmKioq4OSEyMjJw6dIliEQi+Pn5wcXFRea4Fy9eBJ/PR+vWrblt27dvR7t27UqUfd+SJUtw+fJlHDp0SO445Xl+ynv9SjvPsl4LQDqW4+rVqygsLESTJk1KJFsZGRm4fPkyCgsL4efnh6ysLMTGxqJfv34ApKu5Hj9+HKNHj+ZiTE5OxqlTp2S6a16+fInLly9j5MiRAAChUAh3d3cu6ShNee8xRb4fKeGQQ01NOOLSc/HVrkjcS8wCAIz2c8X3vTxhoEvzapDqq6ioCLdu3YKvry90dHQ0OuFQp+Iv1E8//VTdoahcYWEhdu7cieHDh8s92PhjeX5evXqFW7duldt9payEg7pUPlKhkUn44WA0ckVimBvysWhQY3RraKfusAipkI6ODvz9/dUdBtEgfD5f5QMiNZWDg4NSxsrIgxKOj0yusAhzD93HgTvS67Bb1rHE8k+bwt6s5PXchFRHAoEATk5OSExMrFEtjh9acXcBKR09P8pHCcdHJCZZgC+338Hz9Fxo8YDpnd0xrVM9aNPAUKJBjIyMEB4eTgu3VVFN7yqoKnp+lI8Sjo8AYww7bybglyMPICqSwN5MH8s/bYaWdSwrrkxINaOtrY2GDRuqOwxCiIIo4ajhsgsK8f3B+zhyV7rgUCcPGywZ0gSWRpq98BX5eNXUQdyE1HSUcNRgD15lYdqOSLxIz4W2Fg/f9WiAiW3daG4NotGMjY2RkJBQ5iyjhJDqiRKOGogxhj0RCZh7SNqF4mCmj5UjmsPH1ULdoRFSZTweD6ampnLP6EkIqR4o4ahh8kVizD10H/tuS1eU7ORhg2VDm8DckLpQSM2QnZ1NXSqEaCBKOGqQF+m5mLrtNh69zoYWD/i2WwNMbV+XulBIjWJiYoKsrCyYmJioOxSNIxAIcOXKFTRt2rTKy8oToihKOGqIk/eTMXPvPeQIi1DLWBcrhjdD67q0lDypeRhjEAgEMDY2pm4VBeTm5nLTat+7dw8XLlxA/fr11R0W+YjQarEarkgsQdCJGEzZdgc5wiK0rG2JY1+1o2SD1Fg5OTlwdnbmVtIk8omKioKvry+OHz+OWbNm4cSJEx/ssY8cOYLnz59X6RgHDx4ssVy8Kh5H2ZQV08mTJ/H06VMlRKQ+lHBosDc5QozZeBP/XJT+Z57Urg62T2oFW1NaT4LUXKampmCM1ajxGzdu3MCuXbuQlZVVYdmrV69i165dMn/nz5+vsJ6Xlxfu37+P9evXY+/evfjkk0+UEbpcAgMDZWI8fPgwXrx4odAxJk2ahGvXrin0ONWBsmKaOXMmTp06pYSI1Ie6VDTU3YS3mLrtNl5lFcBQVxuLBjdGn8bUJ0tqPrFYjEePHsHDw6PClU2ru9DQUPzyyy/Izc3Fs2fPEB0dDTOz8ldq/vPPPxEVFSWzMqmXlxc6depUbj0LCwv4+flhypQpaNGiBbfi7IfQt29f1K1bl7v/xRdf4Ndffy111VRlPk51UB1jUhdKODTQ7lsvMTf0AURiCdxqGWHtaB+429IAOvJxyM3Nhb+/f41YS0UoFGLz5s0AgGbNmsldr0uXLli7dq1CjyUSibBnzx7Mnz8fc+fOxePHj9GgQYMK6925cweJiYmoX78+PD09AUi7Z7Kzs9GuXTsAQFpaGs6dO4e2bdtyy9XfvHkTjDG0atUK3bp1g6urKwDg1KlTyM/Px40bN6Cvrw9tbW0MGTIEgHRV15s3byIzMxO+vr6wtbUtEU9KSgqioqJgbGwMPz8/maTz3ccBpN0wPj4+0NPTK7MOALx+/Ro3btyAg4MDmjRpgvPnz6Nu3bqljnGpyrkrElNKSgquX7/OxVSa7OxsXLt2DUKhEK1atZJ5vhSN84NgpEJZWVkMAMvKylJrHMJCMfvh4D3m+t1R5vrdUTYx5BbLyhepNSZC1C0/P589fPiQ5efnqzuUSouMjGQAWHR0dIVlBw0axAYMGMAOHjzIrly5wgQCgVyPsXv3bmZlZcWEQiFr164d+9///ldu+YKCAtaxY0fm7OzM+vfvz7y9vdmAAQOYRCJhwcHBrE6dOlzZVatWMQBs3rx53LYWLVqwRYsWMcYYc3V1ZevXr2eMMfbdd98xAwMD1rJlSzZs2DA2evRoxhhjd+7cYbVr12Zubm6sV69erHbt2mzr1q3c8aysrFj37t2Zm5sb69OnD7O3t2ft2rVjRUVFXJl3H0feOgcOHGD6+vrMz8+PtWvXjjVr1ow5ODiwlStXlvq8VOXc5Y3p6NGjzNDQkPn5+bG2bduypk2blojp7NmzzNzcnPn4+LD27dszAwMDtmbNmkrHWZ7y3mOKfD9SwiGH6pBwpAoK2JA115jrd0dZ7dlH2cpzT5hYLFFbPISoS2FhIbt27RorLCxkjJXxYSiRMCbMUc+fRPH3paIJh62tLevVqxfz9PRkVlZWbP/+/RXW6969O5s+fTpjjLHNmzczOzs77jkszdGjR5mFhQXLycnhtoWGhjKxWMyePXvGALD4+HjGGGODBw9mLVq0YB07dmSMST8ztbW12c2bNxljJb90HR0d2aZNm7j7BQUFzNnZmY0bN4774s3Pz2cnT57kylhZWbGWLVuy3NxcxhhjKSkpzMjIiB04cIArU9qXe3l18vLymJ2dHfv555+5On/++ScDUGbCUdVzryim/Px85ujoKBPT4sWLZWLKz89nzs7ObObMmVyZzZs3Mz09Pfb8+fNKxVkeZSUc1KWiAe4lvsXkrbeRnFUAEz0dLB/eFJ08SjY1EvIxyM/Px5AhQxATE1P2XByFecBCNY1p+v4VoKu6lWynTJmCnTt3gs/nAwB+/PFHjB49Gj4+PjJN9+9KSEjAmTNnsGjRIgDAkCFD8NVXX+H48ePo169fqXUMDAwgFArx+PFjNG/eHADQv39/AEDdunXh7OyMCxcuYOzYsbh06RI2btyIwYMHQygU4vLlyzAyMuLqVeT8+fNISEjA77//znUt6Ovro3v37jLlxo4dC0NDQwCAjY0NPD098fjx43KPXV6dK1euICUlBd9++y1X/osvvsAPP/xQ5vGUce4VxfTq1SuZmAIDA/HTTz9x969evYqEhATMmTOH2zZmzBjMnTsXBw8exIwZM5T+GikDXaVSzR2MTMSQteFIziqAm7URQqe1oWSDfNRMTEyQmJj40U781aVLFy7ZAICffvoJhYWFOHfuXJl1Nm/eDGtrazx8+BC7du3C4cOH0bBhQ2zcuLHMOp06dcJXX32FTp06oUGDBvjiiy8QFRXF7W/fvj0uXLiABw8eoLCwED179oSTkxPCw8MRFhaGtm3byj2o9+XLlzA1NS11zMa7LC1lV7jW09NDQUFBpeskJCTAyspK5v+Srq5uhZOiVfXcy4vp5cuXJWLS09OTiSk+Ph7m5uYyx+HxeHBzc0N8fLzS4lQ2auGopsQShj9OPsK6S9JLXjt72ODPT5vCVJ9fQU1CaraioiKcO3cOnTt3ho5OGR9hfENpS4M68A0/6MPp6upCX18faWlppe5njGHz5s3w8vJCaGgot93R0RGhoaFITU2FjY1NqXWDgoIwf/583L59Gzt27EDLli0RGRmJhg0bokOHDliwYAF8fX3Rvn17aGlpoUOHDggLC0NYWBiGDh0q9zmYm5sjJycHQqEQenp6Cp1/VVhaWkIgEEAikUBL67/f35mZmeXWU+a5v8/KyqrCmGrVqoXs7GwUFRXJvAcyMjJQq9Z/czCpMs7KoBaOaii7oBATQ25xyca0jvWwfkwLSjYIAVBQUIAZM2aU/8uWx5N2a6jjT8mzn165coWbx6GgoABv376V2X/mzBlkZ2fD19e31PphYWF4+fIl9u/fLzN3x969e+Hu7o4tW7aUWu/169cQi8Xg8/nw8/PDihUrYGJigtu3bwOQfpnFx8dj8+bN6NChA7ft8OHDiIyM5LaVxtjYWOb1a9++PXR0dLB9+3aZcmUlUcri6+sLxhhOnjzJbbt69SrevHlTbr2qnLs8MQHA8ePHuW0XL15ERkaGTBkdHR0cOnSI2/bo0SM8ePAAbdu2/SBxVga1cFQz8W9yMTEkAk9Tc6Cno4UlQ5qgbxOaX4OQYsbGxnjw4IG6w1CKx48fIzIykmsGP3nyJO7fvw8fHx/ukswlS5bg7du36NSpE3Jzc9G2bVv069cPDRo0wNOnT7Fq1SqMHj26zHk4Nm7ciPbt28PCouRq0QEBAdi0aRNmzpxZYt/Vq1cxf/58DB48GC4uLrhw4QJ4PB73OMVjBCIiIrBhwwYAQMeOHTFq1CiYmpqWOzagRYsW2LhxIwwNDWFgYIAhQ4Zg8eLFmDp1KqKjo+Hl5YVr166hVq1aWLx4sWJPqgIcHR3x1VdfYfTo0fjuu++go6ODFStWVDhtflXOvSIODg745ptvMHr0aMyePRs6OjpYvnw5jI2NuTL29vb4/vvvMX78eDx69AgmJiZYunQpAgIC0LFjxw8SZ2VQwlGNXItNxxfb7+BtXiFsTfWwfkwLNHYyV3dYhFQrhYWFCA0NRUBAgMxYBk309OlTrptj2LBhiIiIQEREBCwtLbmEo127dsjNzQUgbW4PDw/H5s2buS/kgwcPokuXLqUeXyKRgM/nY9q0aaXuHzFiBGJjY5GcnAx7e3uZfYMGDYKnpyd27dqFixcvol69eoiKiuLmcACAb775BlFRUWjUqBEA6Zfl1KlTYWdnJzM24P3Jr1asWIHVq1fj3LlzYIxxg1hbtGiBvXv34ubNm+jYsSNGjx7N1Rk4cGCJQbFdunSBt7d3mY8jT51FixbB29sbFy9ehL29PQ4fPozu3btXOEaosucuT0xBQUHw8vLChQsXYG9vj6NHj2Lr1q1wd3fnyvz0009o1qwZjh07BpFIhJ9//hljxoypdJwfAo8xxj7oI2oggUCg8uWwd958ibmh91EkYWjiZIZ1Y1rQFOWElCI3NxddunTB2bNnYWRkhIKCArx48QJ16tSBvj69Z4hiMjIyZAZfRkVFoVmzZoiOjpZJAj5m5b3HFPl+1LgxHHv27EGbNm1Qr149DBgwADExMXLXXbZsGZycnPDLL7+oLkAFiSUMvx17iDkHolEkYejXxAG7J/tTskFIGYyMjBAeHg4jI9Vdeko+Hnv27MHQoUOxbt06LFy4ED169MDw4cMp2VABjUo4Dh48iJEjR2L06NHYv38/TExM0L59e7kGFt28eRMrV64En88vMehKXXKFRZi89TbWX5YuYvRNF3cs/7Qp9PmavT4EIaokEokQHBwMkUik7lBIDTBlyhQMHz4c9+7dQ0pKClasWIEdO3aoO6waSaMSjgULFmDs2LGYMmUKmjRpgo0bN4IxhjVr1pRbTyAQYMSIEQgODq5wYaQPJTkrH0PWhuNsTAp0dbSwYngzTO9Sv9yBSoQQ6RiOvXv3orCwUN2hkBpiwIABWLVqFZYvX/7BLxX9mGhMwiEQCBAZGYmuXbty23R0dNC5c2dcunSp3Lqff/45AgIC0LlzZ1WHKZf7SVkI+PsqHiYLYGWki52T/NCPrkQhRC5GRkY4deoUdakQomE0JuFISkoCANjZ2clst7W15faVJjg4GI8ePcJvv/0m92MJhUIIBAKZP2UpEksQuDMSKQIh3G2NEfplG/i4lrxcjRBSOqFQiGXLlkEoFKo7FEKIAtSacMyYMQNOTk7l/sXFxQGQXt4FoMTMgnw+H2KxuNTjx8TEYPbs2di+fbtCs9cFBQXBzMyM+3N2dq7cCZZCR1sLKz5thu4NbbFvams4W37YWQkJ0XRisRjh4eEl3vd0wR0hqqGs95ZaL4vNzMzkri8vi729PbS1tZGamgpbW1uEhoZyCwgBwLhx4/D06VNcvXq1RN1Vq1Zh5syZMlO9pqamQk9PD2ZmZoiPjy/1OmShUCjz60kgEMDZ2Vmll8USQipHLBbjyZMnsLGxgZWVlbrDIaTGefPmDVJTU+Hu7l7iO1ORy2LVOvGXhYVFqbPflcbGxgaurq64evWqTMJx5coV9O3bt9Q648ePR0BAgMy27t27w9/fH7/88kuZk57o6el90Pn8CSHyEwqFCAoKwpw5c6CnpwdtbW2Ym5sjNTUVAGBoaEiDrwlRAsYY8vLykJqaCnNz8ypPFKZRM41++eWX+P333zFy5Eg0atQIq1atwsuXL/H5559zZf73v//hzp073KRA7w8s4/P5MDY2lpktjxCiOSQSCRITE7luVuC/sV3FSQchRHnMzc1LjJ+sDI1KOL799lskJyfDz88P2traMDU1xa5du+Dp6cmVycjIwOvXr9UYJSFElQwMDBAcHCyzjcfjwd7eHjY2NnS5LCFKxOfzlTYFukZObS4SiZCVlYVatWqVaDrNzMyESCSCra1tqXXfHcMhrw8xtTkhRD4FBQX4/vvvsXDhQprKnBA105gxHJWlq6sLa2vrUvdVNCbExsZGFSERQgghpBwa2cLxoVELByGEEFJSjW/h+NCKczJlTgBGCKmc/Px8zJo1C4sXL4aBgYG6wyHko1b8vShP2wW1cMghMTFRqZN/EUIIITVJQkJChVd/UsIhB4lEglevXsHExEQp1/cXTySWkJBQI7po6HyqNzqf6o3Op/qqSecCqOZ8GGPIzs6Gg4MDtLTKn7yculTkoKWlpZJ5O0xNTWvEf+JidD7VG51P9UbnU33VpHMBlH8+8l71qTGLtxFCCCFEc1HCQQghhBCVo4RDDfT09PDzzz/XmPVa6HyqNzqf6o3Op/qqSecCqP98aNAoIYQQQlSOWjgIIYQQonKUcBBCCCFE5SjhIIQQQojK0TwcVXT16tVSp3R1cHCAm5tbuXXj4+ORlpaGBg0awMTEpNJllEUoFOLWrVul7nNzc4ODg0O5dR8/fgwzMzO4uLiUmCAtOjoaWVlZMtusra3RoEGDqgdehqysLERHR5e6r2HDhmUu9Hfz5k2IRCKZbS4uLnBxcZHZJpFI8PDhQ4jFYnh7eyttCeeyvH79Gs+ePSt1n4+PT6nTfMtTRyKR4Nq1ayX2u7u7f5DFDjMzMxEXFwdjY2PUqVMHOjoVfyzl5eUhJiYG5ubmqFu3bqXLqEJKSgoSExNhaWmJ2rVryzVZYGJiIjIyMlC3bl0YGRnJ7MvIyMDDhw9L1GnZsiV0dXWVFndZEhISkJKSAjs7uwrnI0pKSsKLFy9ktmlpaaF169YlyqampiI+Ph61a9cuczFOZWOM4cWLF8jMzISzs3OF/7/Dw8MhFotLbLezs0O9evUAALGxsUhOTpbZb2RkhGbNmikv8DIUFhbi+fPnKCgoQMOGDeV67wDAs2fPkJWVBS8vrzKXB5CnjMIYqZL27duzNm3acH8+Pj4MAJs3b16ZdXJzc1mfPn2YoaEh8/DwYIaGhuyff/5RuIyyJScny5xLmzZtWMOGDRkAtn379lLrCAQC9tVXXzELCwvWuHFjZmtry7y9vdmdO3dkyrVp04Y5OzvLHPvnn39W6fncunWrxPm4ubkxAOzGjRtl1rO1tWXu7u4y9datWydT5sGDB6xevXrMzs6OOTo6MhcXF3br1i2Vnk9oaGiJ83FwcGDa2tosNTW10nWys7MZANa4cWOZcsePH1fp+YjFYjZp0iRmYGDAmjVrxpycnJiLiwu7cOFCufV27drFTE1NWf369ZmJiQlr3749y8jIULiMsuXk5LC+ffsyIyMj1rx5c2Ztbc28vb3Zw4cPy6xz+PBh5u3tzRwdHVmjRo2YkZER++mnn2TKHDx4kGlpaZV4Hct6zZUlKSmJ+fn5MQsLC+bj48PMzMxY27ZtWUpKSpl1Fi9ezIyMjGTi7NSpk0wZiUTCAgMDmZ6eHvPy8mJ6enrs22+/Vem5MMbY3bt3mYeHB7O1tWXNmzdnhoaGbODAgSwvL6/MOl27dpU5F19fXwaAzZkzhyszefJkZmVlJVNu9OjRKj+fw4cPM3t7e1a3bl3m7u7OHBwc2MWLF8utk5aWxlq3bs3MzMxYvXr1mJmZGTtw4IDCZSqLEg4l27hxI+PxeOzFixdllvn6669Z7dq1uTfuzp07GY/HY5GRkQqV+RC+//57ZmZmxnJzc0vdHxsby5YvX869aUUiEfv000+Zi4sLk0gkXLkPkWDIY+jQoczT07PcMra2tmzr1q1l7heLxaxhw4Zs8ODB3DmOHTuWubq6MqFQqNR4K9K8eXPWt2/fKtUpTjjCw8OVHV65duzYwbS1tVl0dDRjTPq8jhkzhrm6upZZ5/nz50xXV5etXr2aMcZYVlYW8/LykvmAl6eMKsyfP59ZWVmxV69eMcYYEwqFrH379qxz585l1lm+fDm7f/8+d//ixYuMz+ez3bt3c9sOHjzIjIyMVBd4GcaMGcO8vLxYTk4OY4yxt2/fsnr16rHJkyeXWWfx4sWsSZMm5R43ODiYGRsbs3v37jHGpD8M9PT0yvxRoyytW7dm3bt3Z4WFhYwxxuLj45mZmRlbvHix3MfYsWMHA8AePXrEbZs8eTIbNGiQ0uMtz6tXr5iRkRFbsGABt23lypXMysqKvX37tsx6gwcPZs2bN+de08WLFzMDAwOWmJioUJnKooRDydq0acO6detW5n6xWMwsLCzY77//LrO9fv36bPr06XKX+RCKioqYo6Mj++KLLxSqd/bsWQaAxcfHc9vatGnDvvnmG3bz5k2WkJCg7FDlkp6ezvT09NiyZcvKLWdra8uWLVvGbt26VeqvuWvXrjEALCoqitv27NkzBoCdOHFC6XGXJTIykgFghw8frlKd4oRj7969LCIigmVmZqog2pKWL1/OLCwsZLatWbOGGRkZySSr75o/fz6zsbFhYrGY27Z27Vqmp6fHfUDKU0YVAgMDWYsWLWS2ffvtt6xx48YKHadp06YsMDCQu1+ccDx8+JDdu3eP5efnKyXeivTs2ZN9+umnMtv69+/PBg4cWGadxYsXM29vb3b37l0WExPDRCJRiTKtW7dmo0aNktkWEBBQbmKmDPXr12c//vijzLaGDRuy//3vf3Ifo3PnzuyTTz6R2TZ58mTWu3dvdvv2bfbs2TOZ/3eqcuDAAQZAptVOJBIxbW1ttmHDhlLrZGRkMG1tbbZt2zZum1AolEm65ClTFZRwKNGjR4+4D+6yxMbGMgDszJkzMttHjBjB/UeWp8yHcPToUQZA4VaVoKAgZmhoKPNh06ZNG2ZgYMCaNGnCTE1NWePGjT94a81ff/3FdHV1WXp6ernlbG1tmampKWvSpAkzNDRkHTt2ZHFxcdz+v//+m+no6JT4UrS0tGQLFy5USeylmTZtGnNwcGBFRUVVqlOccNSqVYt5e3szXV1dNnToUJUnHunp6axRo0Zs4sSJ7PTp02zr1q2sTp06bMWKFWXWGTBgAOvevbvMtjt37jAA7ObNm3KXUYXHjx8zZ2dn9uOPP7IzZ86wtWvXMjs7OxYaGir3MTIyMpiJiQlbvnw5t+3gwYMMAKtbty6rW7cu09fXZ/Pnz1fFKci4dOkSs7W1ZX/++Sc7e/Ys++OPP5idnV25z+HixYsZj8djnp6ezMnJiZmbm5f4AjQwMGB//vmnzLaFCxcyS0tLVZwGZ+fOnczR0ZEFBwezM2fOsFmzZrE6derIvLfL8+LFC8bj8Uq0fk6ePJnp6OiwRo0aMWtra+bk5MSOHTumilPghIWFMQBc6yBjjCUkJDAAZf5AvHjxIgPAYmJiZLZ/8sknbMSIEXKXqQpKOJRo1qxZzNrautSsvtitW7cYgBJjHAIDA5m3t7fcZT6EAQMGlPjFVpHo6GhmbGxc4gNx69atXLdMTk4OGzBgAHN1dWXZ2dlKi7cijRo1KvGLrTQbNmzgml1TU1NZq1atmJ+fH5dgLFiwgNnY2JSoV79+fTZz5kzlBl2GgoICZmFhUeIXW2Xq5Ofnsy1btnD3Y2NjWZ06dZTyAVORf/75h1laWrImTZowBwcH1qlTJ/b8+fMyy7dv375EXPHx8QwAO3nypNxlVKGoqIjNmzePmZmZsWbNmrFatWqxYcOGsTdv3shVXyKRsMGDBzMXFxeZZvGIiAgWERHB3T9+/DjT0dFhGzduVPo5vCs3N5dNnTqVmZubs+bNmzMzMzM2ffr0cltYzpw5w549e8bdX716NePxeNzYgvz8fAZA5v8bY9L/Bzwer8yWLWVISUlhAwYMYNbW1qxZs2bMzMyMBQUFyf2Yc+fOZRYWFiXO/8CBA1xLaFFREZs1axYzNDSUeR6UrbCwkPn6+rLGjRuzvXv3soMHD7JWrVoxMzOzMj/jihPX98f+DBgwgPXo0UPuMlVBl8UqSVFREbZs2YJx48aBz+eXWa54X0FBgcz2/Px8bsS5PGVULTU1FUePHsWkSZPkrhMbG4sePXqgb9+++OGHH2T2jRo1CoaGhgCkI7iXLVuG+Ph4hIeHKzXusty6dQvR0dFync9nn33Gjfa2trbGr7/+iuvXryMuLg6A9PV5/7UBPuzrc+DAAbx9+xYTJkyoch19fX2MHj2au+/m5ob//e9/2LdvH4qKipQW8/s2bNiAGTNm4OLFi4iKikJCQgIaNGiAjh07lvr8AqU/9/n5+QAg8/6pqIwq/Pjjj1i/fj0ePnyIO3fuICEhATk5OejXr59c9b/66iucP38ehw8flll908fHBz4+Ptz9nj17on///ti1a5fSz+FdY8eORUREBOLj43H79m3ExcUhLCwMU6ZMKbNOly5dZK4Imjp1Kho3bozdu3cDKP+zTUdHR64reiqDMYaePXsCkF4RdOfOHURHR2P58uX47bffKqwvkUiwefNmjBo1Cvr6+jL7BgwYwF3toq2tjaCgIPD5fBw+fFj5J/IvHR0dnDt3DoMHD8bmzZuxceNG/Pjjj3B3dy/zipLq8N1DCYeSHD16FCkpKZg4cWK55VxdXQFILx97V1JSEnfZpTxlVG3Lli3Q1dXF8OHD5Sr//PlzdOzYEW3btsXWrVuhpVX+fy1bW1sAJc9RVTZs2IC6deuiY8eOCtd9P1ZXV1cIBAJkZ2dzZUQiEdLS0j7Y67NhwwZ07doVtWvXVkkdW1tbiEQipKenVz7IChw9ehTt27eHt7c3AOnlk1OnTkV8fDzu3btXah1XV9dS3xcAZN4/FZVRhaNHj2LgwIHc5eP6+vqYMGECrl69ioyMjHLrfvPNN9i+fTvOnDmDJk2aVPhYtra2Kn/vHD16FKNHj+aWMTc3N8fIkSMV/iJ9N1ZtbW04OjqW+voUf+6pwqtXr3Dnzh1MmjSJ++J0dnZGv3795Dqf06dPIyEhQa4fLNra2rCyslL562NiYoK5c+fi6NGjOHz4MLp06YIHDx5w76f3VYfvHko4lCQ4OBjt27eHu7t7iX1PnjzB3bt3AUjftC1atJD5T56VlYWwsDB07dpV7jKqtmHDBnz66aelzv1x7949PH78mLsfFxeHjh07wt/fH9u2bSsxH0V+fj4kEonMttOnTwNAmW8OZcrLy8POnTsxceLEUn9B3bhxAy9fvgQA5Obmlth/+vRp6OjowMPDAwDQsWNH6Ojo4MiRI1yZU6dOQSQSoUuXLio6i/+8ePEC58+fL/XDTywW48qVK0hNTZW7TlnnbG1trdJ5OKytrUt8sCUkJHD7imO7cuUKN4dL165dERERITPvwaFDh1CnTh3ul7U8ZVR1PomJiSXOR1dXl/vSfvPmDa5cuSIzz8uMGTMQEhKCM2fOoHnz5iWO+/7rU1RUhLCwMJW/d8o6n3fnzEhMTJSZw+X9WDMzMxERESETa9euXXHkyBFu/iKJRIIjR46o9LPN0tISWlpaFZ7Ps2fPEBkZWaJ+cHAwWrVqhUaNGslsZ4xxrWfFnj59ivj4eJW/Pu8/bkhICMRiMYYNG8Zti46ORkxMDADAy8sLDg4OMt8rz58/R3R0NPfcy1OmSqrcKUNYUlJSiZG97xo5ciTz8fHh7p8+fZrp6OiwuXPnstDQUNa+fXvWoEEDmevB5SmjKleuXGEA2PXr10vd36pVKzZs2DDGGGOvX79mtWvXZt7e3iwsLIxdvnyZ+xMIBIwx6ZwVvr6+bO3atezUqVNsyZIlzMLCgg0fPlzl58IYY5s3b2Y6OjosOTm51P1WVlbshx9+YIxJB8p26dKFbdy4kZ08eZLNnTuX6evrs7lz58rUmTVrFrOysmIbN25kW7duZfb29uzzzz9X+bkwxtiPP/5Y5lihzMxMBoBt2rRJ7jrLly9nw4YNY9u3b2fHjh1jX375JdPR0SlxDGW7ffs24/P5bPz48ezEiRNs8+bNzMXFhfXu3ZsrU3xVTfHcHEVFRczX15f5+vqyAwcOsKCgIKajoyNzGak8ZVThwIEDjMfjsdmzZ7NTp06xVatWMQsLC5kry/bu3csAcFdqzZ07l2lpabHly5fLvHfeHbQ3dOhQNmvWLHbo0CG2b98+1qVLF2Zubi5zOa0qLF68mOnr67M//viDnT59mv32229MV1eXrVy5kisTFBTE9PT0uPtt2rRh8+fPZ8eOHWPbtm1jTZs2lbm8nzHGnjx5wkxNTdm4cePY4cOH2YgRI5iFhYXcgzcra/LkyaxWrVpszZo17NSpU+zbb79lPB5PZr6ZCRMmsIYNG8rUS0tLY7q6uiw4OLjEMYVCIfP09GRLly5lJ0+eZOvWrWO1a9dmvr6+rKCgQKXnM27cOPb777+zM2fOsAULFjADA4MS8wW1b9+e9e/fn7sfEhLC+Hw+W7JkCdu/fz9r2rQpa9u2rcyVNfKUqSxaLVYJ9u3bhzVr1uDYsWMl+vcAYMGCBYiPj0dwcDC37eLFi1izZg3S0tLQtGlTzJ49u8Rse/KUUYW//voL165dw549e0rdP3nyZNjb2+OXX35BZGQkAgMDSy23bt06eHl5AQAePnyINWvW4MmTJ7C3t0efPn0wePBglZ3Du2bNmgWxWIxly5aVur9Pnz7o378/9+v/2rVr2LRpE16+fAkXFxeMGDGiRFeMRCJBcHAwDh06BIlEgl69emHq1Klyz/RXFcOHD0eLFi3w7bffltiXk5ODHj164IcffuD6rCuqAwBHjhzB3r17kZaWBjc3N3z++edyNe1X1YMHD/D3338jNjYWJiYmaN++PT7//HNu+exnz55h3LhxWLlyJTdzY1ZWFhYtWoQbN27A3NwcEyZMkDlXecuowpUrV7j/O7Vq1UKvXr0wcuRIrosxLCwMP/74Iw4ePAhra2tMnDgRjx49KnGcjh07YsGCBQCks/iuX78e58+fh0Qigbe3N6ZPn/5BPgsOHTqEffv24fXr17C3t8ewYcPQu3dvbv/27duxceNGnDt3DoD0ef/7779x7do1GBgYoEWLFpg2bVqJ2VMfPnyIpUuX4sWLF6hbty5mzZpVauuwMonFYmzZsgWnTp3Cmzdv4OrqiokTJ8LPz48r88cff+DBgwfYsmULt+3w4cNYtmwZjh49CmNj4xLHTU5OxsqVKxEZGQkLCwu0a9cOEydOLHcsnzLk5ORg0aJFuH79OhwcHDBx4kS0bdtWpkxgYCBMTEywcOFCmfMJCQmBQCCAv78/Zs2aVaIlW54ylUEJByGEEEJUjsZwEEIIIUTlKOEghBBCiMpRwkEIIYQQlaOEgxBCCCEqRwkHIYQQQlSOEg5CCCGEqBwlHIQQQghROUo4CCFEze7cuYMLFy6ApkUiNZnqp0UkhBBSpnXr1mHJkiXcGkqrV69Wd0iEqAQlHIQQtUhMTMStW7cwYMAAANIF5qKjo+Vezl2TPHz4kFsB18fHB/Xr1+f27d27FydPnoSrqys8PT1L1E1MTMSVK1cAAO7u7qUu8EaIJqAuFUJqgPPnz2Pfvn0adezr16/LrF578eJFfPHFF3LVff78ucxqvdXdgQMHMG3aNISGhuLFixcy+zp27Ijff/8dQUFBpSYTSUlJCA0NxcyZM7Fx48YPFTIhSkctHIRoOIFAgL59+yIvLw/Xr19Hq1atNOLY76tTpw769+8vV9nz58/j119/Rd++fVUWj7LVq1cPu3btKrF9woQJsLe3B4/Hw9GjR0vsb9WqFXbt2oWAgIAPECUhqkMJByEabseOHXBwcECjRo0QHBxcIil4+vQpnj59iq5du+L69etITk7GgAED5FrNUpnHZowhPDwcWVlZaNy4cYn9Li4u6N69u8y2N2/e4Pbt29DR0UGLFi1gamqKpKQk3Lp1C7m5udwXeJMmTWBhYYGwsDAAgIGBARo0aAAPD49S4+3WrRvu3buH5ORkNG3aFI6OjiXiSUhIQGRkJGxtbdGiRQtoa2vL7I+Pj0dUVBSsrKzQvHlzGBoaVvh8lmbr1q3w9PRE06ZNsWXLlg+yqi0h6kAJByEaLjg4GBMnTkTTpk0xePBg/PnnnzLLaJ86dQq//vorXFxcoK+vDwcHB/Tp00euhENZxxaJROjTpw8iIyPh6+uLe/fuwcvLS6bMxYsX8eOPP3JjOEJDQzFmzBj4+PhAV1cXsbGx+Oeff2BmZobIyEjk5eUhNDQUAGBoaAhXV1fufm5uLq5evYqePXti27Zt4PF4XLy//fYbateuDT09PQDAjRs3sHv3bu5xGWOYPn06l2AVFBRAR0cHx44dg6mpKRhj+Prrr7Ft2zb4+fkhPT0dycnJ2L9/P3x9fRV45aQ2bdqECRMmoGnTpujduzcyMzNhYWGh8HEIqfYYIURj3b17l/H5fPb69WsmkUhYnTp12IYNG2TKrFy5kgFgmzZtUtuxly9fzuzs7FhycjJjjLH09HTm6urKrKysuDKbNm1ijo6O3P0WLVqwhQsXcvffvHnDwsLCGGOMrV+/nrm6upb7mGlpacze3p7t27evRLzvbps5cyZr3Lgxd//vv/9mhoaG7O7du9y2a9eusdTUVMYYY+vWrWN169Zl6enp3P5ff/2VeXh4lBnLggULWKtWrUpsDw8PZ3w+n6WmpjKJRMLc3NzYqlWrSj1G//792ZdfflnuORNSndGgUUI02Pr169GvXz/Y2tqCx+NhwoQJCA4OLlHO1NQUY8eOVdux9+zZg1GjRsHOzg4AYGVlhQkTJpRbx8DAAE+ePEFBQQEAwNLSEu3bty+3TmFhIa5fv44DBw7g7NmzcHZ2xs2bN2XKWFtbY9CgQdz9Dh064PHjx9z9kJAQjBkzRqbbx9/fH9bW1gCkLRKNGzfGhQsXsHfvXuzZswfGxsZ49OgRkpOTy43vfRs2bEDfvn1hbW0NHo+Hzz77DBs2bFDoGIRoCupSIURDCYVCbN++HaNGjeLGMhgZGSE8PBwxMTEyl1ja2dlx3QrFShvAWKdOHbRq1arKx37fy5cvMXz48BKPVZ6VK1di0qRJsLa2RuvWrdGvXz9MnDiR6wp5X3R0NHr16sWN3zAyMkJKSgpSU1NlyllaWsrc19PTg1AolIn1008/LTOuuLg4iESiElfuDBs2DCKRqNxzeldeXh727NmD0aNHc8+xgYEBIiMjcffuXTRp0kTuYxGiCSjhIERD7d+/Hzo6OkhNTeXGLgCAh4cHgoODsXTpUm5baQnBu3WKtW/fHq1atarysd9nZWWFzMxMmW3v339fkyZNcPPmTSQnJ+PcuXOYN28erl69ih07dpRa/ocffkCnTp0QEhLCbevRo4fCs3eam5vjzZs3Ze43NTVFp06dsGjRIoWO+769e/dCT08P6enpJZ7jDRs2YMWKFVU6PiHVDSUchGio4OBgjBkzBkuWLJHZvnfvXnz55ZcICgqCrq5umfVLa+FQ1rHf17ZtWxw6dAg//PADl6AcOHCg3DpJSUlwdHSEvb09Ro0ahczMTPz5558AAGNjY66rpdjr16/h5+cnc//q1asy3Sfy6NatG3bv3o2ffvqJO8fc3FzweDwYGhqiR48e2L59O+bOnQsTE5MS8cpr48aNZT7HU6ZMweLFi8tszSFEE1HCQYgGev78OcLCwjB//vwS+3r27AmBQIDDhw9j8ODB1eLYs2fPRuPGjdGvXz/07dsXZ86cwf3798utM3DgQDRs2BCtWrVCfn4+li1bxnXLNG/eHGlpafj1119Rr149NGnSBAEBAVi6dCn09fWho6ODlStXQkdH8Y+4uXPn4uTJk/D398fYsWNRUFCAnTt34vjx4zA0NMRPP/2Ec+fOwdfXF5MmTYKhoSFu3LiB58+f49KlS3I9RmxsLC5fvozffvutxL6ePXsiNzcXhw4dwtChQxWOn5DqihIOQjRQTEwMRo8ejdatW5fYZ2xsjDlz5nADGN3d3dG7d2+1HtvR0RERERH4+++/ERERgS5duuDrr7+W6f54f+KvK1euYPv27QgPD4euri6WL1/O7Xd3d8fx48dx+PBhPHjwAIaGhpg9ezacnJxw8eJF6OnpYfny5Xjx4gW0tP4bG19avHZ2dhg2bBh3v1atWoiIiMDGjRsRGRkJBwcH7Nu3D/b29gCkY0Bu3ryJbdu24caNG9DX10eXLl3KHfcBSOcU2bVrF3x8fPD48eMKn+OUlBQA/01tnpSUBCcnpwqfa0KqKx5TtIOTEEKIQg4ePIjdu3cDAD777DN069ZN7ro3btzgupK6d++O8ePHqyRGQlSNEg5CCCGEqBzNw0EIIYQQlaOEgxBCCCEqRwkHIYQQQlSOEg5CCCGEqBwlHIQQQghROUo4CCGEEKJylHAQQgghROUo4SCEEEKIylHCQQghhBCVo4SDEEIIISpHCQchhBBCVO7/cv2pUC1dQMoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 550x320 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cluster = molecule(\"CH3CH2OH\") + molecule(\"CH3CH2OH\"); cluster.positions[9:] += [4.0, 0.5, 0.3]   # two ethanols\n",
    "ref = reference(cluster, cutoffs=dict(disp2=12.0, disp3=9.0, cn=10.0))\n",
    "mine = xnn_d4(cluster, cutoff_pair=12 * BOHR, cutoff_triple=9 * BOHR, cutoff_cn=10 * BOHR)\n",
    "print(f\"cutoffs 12/9/10 bohr: dftd4 {ref['energy']:.12e}   xnn {mine['energy_au']:.12e}   |diff| {abs(mine['energy_au'] - ref['energy']):.1e} Eh\")\n",
    "\n",
    "# smooth switching: a dimer scanned through the pair cutoff\n",
    "d = np.linspace(7.0, 9.0, 81)\n",
    "sharp = D4Dispersion(s9=0.0, cutoff_pair=8.0, cutoff_cn=8.0, cutoff_eeq_cn=8.0, cutoff_triple=8.0)\n",
    "smooth = D4Dispersion(s9=0.0, cutoff_pair=8.0, switch_width_pair=1.5, cutoff_cn=8.0, cutoff_eeq_cn=8.0, cutoff_triple=8.0)\n",
    "E = {m: [float(xnn_d4(Atoms(\"Ar2\", positions=[[0, 0, 0], [x, 0, 0]]), model=mod)[\"energy\"]) * 1e3 for x in d]\n",
    "     for m, mod in [(\"sharp cutoff (upstream)\", sharp), (\"1.5 Å switching window\", smooth)]}\n",
    "plt.figure(figsize=(5.5, 3.2))\n",
    "for k, v in E.items(): plt.plot(d, v, label=k)\n",
    "plt.axvline(8.0, color=\"k\", ls=\":\", lw=0.8); plt.xlabel(\"Ar–Ar distance [Å]\"); plt.ylabel(\"E_disp [meV]\"); plt.legend(); plt.title(\"pair cutoff at 8 Å\"); plt.tight_layout(); plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cb6a763",
   "metadata": {},
   "source": [
    "### The one place the two codes differ: isolated atoms, at $10^{-9}$ e\n",
    "\n",
    "An atom with no neighbor inside about two covalent radii has an EEQ coordination number\n",
    "of exactly zero. Upstream evaluates the soft cap $\\mathrm{CN}' = \\ln(1+e^{8}) - \\ln(1+e^{8-\\mathrm{CN}})$\n",
    "with the two logarithms coming out one ulp apart (they take different paths through the\n",
    "Fortran runtime), leaving $1.8\\\\cdot10^{-15}$\n",
    "instead of $0$; the EEQ right-hand side $\\kappa\\,\\mathrm{CN}/\\sqrt{\\mathrm{CN} + 10^{-14}}$ then amplifies\n",
    "that by $10^{7}$, so the charges of such systems differ by $\\sim10^{-9}$ e and the energies by\n",
    "$\\sim10^{-13}$ hartree. xnn evaluates the cap exactly and gets zero. The effect is only\n",
    "visible for genuinely isolated atoms (a stretched cluster, a line of far-apart atoms),\n",
    "never inside a molecule or a solid, and it is well below any physical relevance.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "74bb0189",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:20:17.198317Z",
     "iopub.status.busy": "2026-09-23T19:20:17.198217Z",
     "iopub.status.idle": "2026-09-23T19:20:17.676081Z",
     "shell.execute_reply": "2026-09-23T19:20:17.675016Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "isolated atoms:  dftd4 CN = [0. 0. 0.]  xnn CN = [0. 0. 0.]\n",
      "                 |dq| = 1.4e-09 e   |dE| = 2.7e-15 Eh\n",
      "ln(1+e^8) in double: 8.000335406372896;  one ulp at this magnitude is 1.8e-15, the residual upstream's cap leaves for CN = 0\n"
     ]
    }
   ],
   "source": [
    "line = Atoms(numbers=[1, 7, 6], positions=[[i * 12 * BOHR, 0.0, 0.0] for i in range(3)])   # three atoms 12 bohr apart\n",
    "ref, mine = reference(line), xnn_d4(line)\n",
    "print(\"isolated atoms:  dftd4 CN =\", ref[\"coordination numbers\"], \" xnn CN =\", mine[\"coordination_numbers\"].numpy())\n",
    "print(\"                 |dq| = %.1e e   |dE| = %.1e Eh\" % (np.abs(mine[\"eeq_charges\"].numpy() - ref[\"partial charges\"]).max(),\n",
    "                                                        abs(mine[\"energy_au\"] - ref[\"energy\"])))\n",
    "import math\n",
    "print(\"ln(1+e^8) in double: %.17g;  one ulp at this magnitude is %.1e, the residual upstream's cap leaves for CN = 0\" % (math.log(1 + math.exp(8.0)), np.spacing(8.0)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d22758b4",
   "metadata": {},
   "source": [
    "## Block 6: One implementation, every channel\n",
    "\n",
    "The same `DFTD4` module serves the batched `AtomicGraph` path (training), the\n",
    "scripted whole-system and pair-style ABIs (`TorchScriptPotential`, LAMMPS) and the\n",
    "ASE calculator; and it wraps a short-range model without touching that model's own\n",
    "neighborhood (the core sees only the edges within its cutoff).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "20b6d3a2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-09-23T19:20:17.677669Z",
     "iopub.status.busy": "2026-09-23T19:20:17.677584Z",
     "iopub.status.idle": "2026-09-23T19:20:19.413886Z",
     "shell.execute_reply": "2026-09-23T19:20:19.413155Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eager           -0.244579870193917 eV\n",
      "TorchScript     -0.244579870193917 eV   |dF| = 0.0e+00\n",
      "pair-style ABI  -0.244579870193917 eV\n",
      "ASE calculator  -0.244579870193917 eV\n",
      "batch vs single |dE| = 0.0\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MACE+D4: cutoff 12.0 Å, E = E_sr + E_disp: 6.332386 = 6.576966 + -0.244580  (|E_sr - MACE alone| = 0.0e+00)\n"
     ]
    }
   ],
   "source": [
    "from xnn.common.config import from_dict\n",
    "from xnn.common.data import build_neighbor_list\n",
    "from xnn.common.deploy import TorchScriptPotential, XNNCalculator\n",
    "from xnn.common.models import build_model\n",
    "\n",
    "atoms = molecule(\"C6H6\")\n",
    "d4 = D4Dispersion(cutoff_pair=12.0, cutoff_triple=9.0, cutoff_cn=10.0, cutoff_eeq_cn=10.0)\n",
    "eager = xnn_d4(atoms, model=d4)\n",
    "scripted = torch.jit.script(TorchScriptPotential(d4, d4.cutoff).eval())\n",
    "ts = scripted(torch.tensor(atoms.positions), torch.tensor(atoms.numbers))\n",
    "ei, cs = build_neighbor_list(torch.tensor(atoms.positions), d4.cutoff)\n",
    "ts_pair = scripted.forward_lammps(torch.tensor(atoms.positions), ei, cs, torch.tensor(atoms.numbers), torch.zeros(3, 3))\n",
    "atoms.calc = XNNCalculator(ForceStressOutput(d4), cutoff=d4.cutoff)\n",
    "print(f\"eager           {float(eager['energy']):.15f} eV\")\n",
    "print(f\"TorchScript     {float(ts['energy']):.15f} eV   |dF| = {float((ts['forces'] - eager['forces']).abs().max()):.1e}\")\n",
    "print(f\"pair-style ABI  {float(ts_pair['energy']):.15f} eV\")\n",
    "print(f\"ASE calculator  {atoms.get_potential_energy():.15f} eV\")\n",
    "\n",
    "# batching: a batch of three molecules equals the three single evaluations\n",
    "graphs = [structure_to_graph({\"pos\": molecule(n).positions, \"atomic_numbers\": molecule(n).numbers}, d4.cutoff) for n in (\"H2O\", \"CH4\", \"C6H6\")]\n",
    "e_batch = d4(collate(graphs))[\"energy\"]\n",
    "e_single = torch.stack([d4(g)[\"energy\"][0] for g in graphs])\n",
    "print(\"batch vs single |dE| =\", float((e_batch - e_single).abs().max()))\n",
    "\n",
    "# wrapping a short-range model: (untrained) MACE + D4, the core keeps its 4.5 Å neighborhood\n",
    "cfg = from_dict({\"model\": {\"name\": \"mace\", \"cutoff\": 4.5, \"n_interactions\": 1, \"n_rbf\": 6, \"n_features\": 8,\n",
    "                           \"extra\": {\"species\": [1, 6], \"l_max\": 2,\n",
    "                                     \"dispersion\": {\"cutoff_pair\": 12.0, \"cutoff_triple\": 9.0, \"cutoff_cn\": 10.0, \"cutoff_eeq_cn\": 10.0}}}})\n",
    "torch.manual_seed(0)\n",
    "wrapped = build_model(cfg.model)\n",
    "out = wrapped(structure_to_graph({\"pos\": atoms.positions, \"atomic_numbers\": atoms.numbers}, wrapped.cutoff))\n",
    "core = wrapped.model(structure_to_graph({\"pos\": atoms.positions, \"atomic_numbers\": atoms.numbers}, 4.5))\n",
    "print(f\"MACE+D4: cutoff {wrapped.cutoff} Å, E = E_sr + E_disp: {float(out['energy']):.6f} = {float(core['energy']):.6f} + {float(out['energy_disp']):.6f}  (|E_sr - MACE alone| = {abs(float(out['energy_sr']) - float(core['energy'])):.1e})\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0866cc51",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "| block | quantity | xnn vs `dftd4` |\n",
    "|---|---|---|\n",
    "| 1 | D4 and EEQ coordination numbers, EEQ charges | ~1e-13 |\n",
    "| 2 | charge-scaled polarizabilities, $C_6$ matrix (any total charge) | ~1e-12 relative |\n",
    "| 3 | two-body + ATM energies, gradients (molecules up to C60) | ~1e-15 Eh, ~1e-14 Eh/bohr |\n",
    "| 4 | periodic energies, gradients, virials (Ewald EEQ) | ~1e-15 Eh, ~1e-13 |\n",
    "| 5 | shorter, matching real-space cutoffs | ~1e-16 Eh |\n",
    "| 5 | isolated atoms (upstream's rounded CN cap) | 1e-9 e in q, 1e-13 Eh (documented) |\n",
    "| 6 | eager / TorchScript / pair-style / ASE / batch | identical |\n",
    "\n",
    "Unit conversions matter at this level: `dftd4` derives the bohr radius from the\n",
    "CODATA-2018 constants ($a_0 = \\hbar / m_e c \\alpha$), which differs from the tabulated\n",
    "0.529177210903 Å in the 12th digit; xnn uses the same derived value so the covalent\n",
    "radii, and with them the coordination numbers, agree to round-off.\n"
   ]
  }
 ],
 "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
}