{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f735281b",
   "metadata": {},
   "source": [
    "# Argon density from MD: `xnn` vs the original MACE\n",
    "\n",
    "This notebook computes the **mass density of liquid Argon** by running **NPT\n",
    "molecular dynamics through ASE** with a trained MACE potential, and compares the\n",
    "`xnn` result against the **original** `mace-torch` model in **two complementary\n",
    "ways**:\n",
    "\n",
    "* **Track (a): same potential.** The trained original-MACE weights are *copied*\n",
    "  into the `xnn` model, so both codes represent the **identical** potential-energy\n",
    "  surface. Any density difference then reflects only the `xnn`-vs-`mace`\n",
    "  inference / MD code path; it should be numerically zero. (Sections **3a/4a/5a**.)\n",
    "* **Track (b): independently trained.** `xnn` is trained **from scratch** on the\n",
    "  same data with **no weight copying**, giving two *independent* potentials. Now we\n",
    "  compare the density as two practitioners would if each fit their own model.\n",
    "  (Sections **3b/4b/5b**.)\n",
    "\n",
    "Pipeline for each track: train → wrap in an **ASE calculator** (energy + forces +\n",
    "**stress**) → run **NPT** MD ($T=85$ K, $P=1$ bar) → measure $\\rho=M/V$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5d5d3e3",
   "metadata": {},
   "source": [
    "## 0. Setup\n",
    "\n",
    "Train in `float32` (speed); run MD in `float64` (smooth forces/stress, as in\n",
    "production MACE)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cd7980cb",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:41:09.684422Z",
     "iopub.status.busy": "2026-07-20T04:41:09.684303Z",
     "iopub.status.idle": "2026-07-20T04:41:12.132979Z",
     "shell.execute_reply": "2026-07-20T04:41:12.132187Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/e3nn/o3/_wigner.py:10: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.\n",
      "  _Jd, _W3j_flat, _W3j_indices = torch.load(os.path.join(os.path.dirname(__file__), 'constants.pt'))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xnn: 0.1.0 | mace (original): 0.3.16 | device: cuda\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "import ase.io\n",
    "import ase.units as u\n",
    "\n",
    "torch.set_default_dtype(torch.float32)\n",
    "torch.manual_seed(0)\n",
    "DEVICE = \"cuda\" if torch.cuda.is_available() else \"cpu\"\n",
    "CUTOFF, SPECIES = 6.0, [18]\n",
    "import xnn, mace\n",
    "print(\"xnn:\", xnn.__version__, \"| mace (original):\", mace.__version__, \"| device:\", DEVICE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a39fa6a7",
   "metadata": {},
   "source": [
    "## 1. Load data and build both data pipelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8c7e74b8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:41:12.134745Z",
     "iopub.status.busy": "2026-07-20T04:41:12.134506Z",
     "iopub.status.idle": "2026-07-20T04:41:19.593522Z",
     "shell.execute_reply": "2026-07-20T04:41:19.592763Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuequivariance or cuequivariance_torch is not available. Cuequivariance acceleration will be disabled.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "200 configs | E0={18: 0.0} | lambda=17.16 | train 180 / val 20\n"
     ]
    }
   ],
   "source": [
    "from xnn.common.data import AtomicDataset, load_dataset\n",
    "from mace.data import AtomicData, Configuration\n",
    "from mace.tools import AtomicNumberTable, torch_geometric\n",
    "\n",
    "E0 = {18: 0.0}   # argon isolated-atom reference energy\n",
    "train_structs = load_dataset(\"argon_md\", split=\"train\")\n",
    "ZT = AtomicNumberTable(SPECIES)\n",
    "\n",
    "xnn_train = AtomicDataset(train_structs, CUTOFF)\n",
    "LAMBDA = float(sum(xnn_train[i].num_edges for i in range(len(xnn_train))) /\n",
    "               sum(xnn_train[i].num_nodes for i in range(len(xnn_train))))\n",
    "\n",
    "def to_mace(s):\n",
    "    conf = Configuration(atomic_numbers=np.asarray(s[\"atomic_numbers\"]), positions=s[\"pos\"],\n",
    "        properties={\"energy\": s[\"energy\"], \"forces\": s[\"forces\"]},\n",
    "        property_weights={\"energy\": 1.0, \"forces\": 1.0}, cell=s[\"cell\"], pbc=(True,)*3)\n",
    "    return AtomicData.from_config(conf, z_table=ZT, cutoff=CUTOFF)\n",
    "mace_train = [to_mace(s) for s in train_structs]\n",
    "\n",
    "# one shared train/val split used by both models\n",
    "g = torch.Generator().manual_seed(0)\n",
    "perm = torch.randperm(len(train_structs), generator=g).tolist()\n",
    "n_val = max(1, int(0.1 * len(train_structs)))\n",
    "val_idx, train_idx = perm[:n_val], perm[n_val:]\n",
    "print(f\"{len(train_structs)} configs | E0={E0} | lambda={LAMBDA:.2f} | \"\n",
    "      f\"train {len(train_idx)} / val {len(val_idx)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8239edcb",
   "metadata": {},
   "source": [
    "## 2. Train the two models\n",
    "\n",
    "Both use the same architecture ($T=2$, $\\ell_{\\max}=3$, $L_{\\max}=1$, $\\nu=3$, 32\n",
    "channels), the same data/loss/optimiser/schedule and the same split. The original\n",
    "MACE is trained with a native loop; `xnn` with `xnn.train.Trainer`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0febc6d7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:41:19.595800Z",
     "iopub.status.busy": "2026-07-20T04:41:19.595487Z",
     "iopub.status.idle": "2026-07-20T04:53:19.861503Z",
     "shell.execute_reply": "2026-07-20T04:53:19.860115Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "original MACE trained 60 epochs in 417 s (final val 5.952e-04)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    0 | train loss 1.7824e-01 | val loss 4.0494e-02\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    1 | train loss 4.9534e-02 | val loss 1.9639e-02\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    2 | train loss 2.4203e-02 | val loss 8.6105e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    3 | train loss 8.9203e-03 | val loss 3.8398e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    4 | train loss 3.4036e-03 | val loss 1.6303e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    5 | train loss 2.2902e-03 | val loss 1.2410e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    6 | train loss 1.7293e-03 | val loss 9.7818e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    7 | train loss 1.5235e-03 | val loss 1.1189e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    8 | train loss 1.3902e-03 | val loss 1.0681e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch    9 | train loss 1.2514e-03 | val loss 9.6917e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   10 | train loss 1.3658e-03 | val loss 1.3300e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   11 | train loss 1.2830e-03 | val loss 1.2296e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   12 | train loss 1.1869e-03 | val loss 1.5758e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   13 | train loss 1.0534e-03 | val loss 7.0099e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   14 | train loss 8.3061e-04 | val loss 7.5308e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   15 | train loss 7.8896e-04 | val loss 1.2321e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   16 | train loss 9.4133e-04 | val loss 1.0322e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   17 | train loss 1.1433e-03 | val loss 9.0873e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   18 | train loss 7.7776e-04 | val loss 6.3342e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   19 | train loss 7.5670e-04 | val loss 1.2568e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   20 | train loss 1.1542e-03 | val loss 1.3643e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   21 | train loss 7.7342e-04 | val loss 6.4181e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   22 | train loss 8.7842e-04 | val loss 9.1290e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   23 | train loss 8.4342e-04 | val loss 1.3029e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   24 | train loss 8.9835e-04 | val loss 6.6383e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   25 | train loss 8.4394e-04 | val loss 6.9654e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   26 | train loss 9.3775e-04 | val loss 1.1299e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   27 | train loss 8.8414e-04 | val loss 5.5828e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   28 | train loss 8.3435e-04 | val loss 7.7909e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   29 | train loss 7.1105e-04 | val loss 5.3847e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   30 | train loss 6.7055e-04 | val loss 7.1859e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   31 | train loss 1.1826e-03 | val loss 8.8792e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   32 | train loss 1.2959e-03 | val loss 1.1714e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   33 | train loss 1.8242e-03 | val loss 1.0036e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   34 | train loss 1.6937e-03 | val loss 9.4435e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   35 | train loss 1.2377e-03 | val loss 1.0633e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   36 | train loss 8.3443e-04 | val loss 9.6946e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   37 | train loss 1.0155e-03 | val loss 8.6592e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   38 | train loss 1.1977e-03 | val loss 7.3642e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   39 | train loss 1.2071e-03 | val loss 7.3105e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   40 | train loss 6.0991e-04 | val loss 4.9166e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   41 | train loss 1.1405e-03 | val loss 1.1974e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   42 | train loss 8.4249e-04 | val loss 1.1444e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   43 | train loss 1.7142e-03 | val loss 1.2634e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   44 | train loss 2.3817e-03 | val loss 1.3921e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   45 | train loss 2.0804e-03 | val loss 1.3230e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   46 | train loss 1.1204e-03 | val loss 5.5318e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   47 | train loss 1.0169e-03 | val loss 6.6129e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   48 | train loss 9.7584e-04 | val loss 1.4237e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   49 | train loss 2.3880e-03 | val loss 1.9722e-03\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   50 | train loss 2.7580e-03 | val loss 7.9063e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   51 | train loss 1.0216e-03 | val loss 8.0398e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   52 | train loss 7.0373e-04 | val loss 4.3472e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   53 | train loss 4.9870e-04 | val loss 4.7063e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   54 | train loss 4.5946e-04 | val loss 4.2057e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   55 | train loss 4.4675e-04 | val loss 4.1498e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   56 | train loss 4.4809e-04 | val loss 4.5636e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   57 | train loss 4.4194e-04 | val loss 4.3432e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   58 | train loss 4.4035e-04 | val loss 4.5914e-04\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch   59 | train loss 4.3256e-04 | val loss 4.1203e-04\n",
      "xnn (independent) trained 60 epochs in 300 s\n"
     ]
    }
   ],
   "source": [
    "from e3nn import o3\n",
    "import torch.nn.functional as Fn\n",
    "import mace.modules as mm\n",
    "from mace.modules.blocks import RealAgnosticResidualInteractionBlock as MIB\n",
    "from xnn.common.config import from_dict\n",
    "from xnn.common.train import Trainer\n",
    "from torch.utils.data import Subset\n",
    "\n",
    "HID = \"32x0e+32x1o\"\n",
    "EW, FW, LR, WD, BS, EPOCHS = 1.0, 100.0, 0.01, 5e-7, 10, 60\n",
    "\n",
    "# ---------- 2.1 train the ORIGINAL MACE (native loop) ----------\n",
    "torch.manual_seed(0)\n",
    "mace_model = mm.MACE(r_max=CUTOFF, num_bessel=8, num_polynomial_cutoff=5, max_ell=3,\n",
    "    interaction_cls=MIB, interaction_cls_first=MIB, num_interactions=2, num_elements=1,\n",
    "    hidden_irreps=o3.Irreps(HID), MLP_irreps=o3.Irreps(\"16x0e\"),\n",
    "    atomic_energies=np.array([E0[18]]), avg_num_neighbors=LAMBDA, atomic_numbers=SPECIES,\n",
    "    correlation=3, gate=Fn.silu, radial_MLP=[64, 64, 64], radial_type=\"bessel\",\n",
    "    use_reduced_cg=False, apply_cutoff=True).to(DEVICE)\n",
    "\n",
    "DL = torch_geometric.dataloader.DataLoader\n",
    "tr_loader = DL([mace_train[i] for i in train_idx], batch_size=BS, shuffle=True)\n",
    "va_loader = DL([mace_train[i] for i in val_idx],   batch_size=BS, shuffle=False)\n",
    "opt = torch.optim.Adam(mace_model.parameters(), lr=LR, weight_decay=WD)\n",
    "sched = torch.optim.lr_scheduler.ReduceLROnPlateau(opt, patience=10)\n",
    "def mace_loss(out, b):\n",
    "    n = (b.ptr[1:] - b.ptr[:-1]).to(out[\"energy\"].dtype)\n",
    "    return EW * (((out[\"energy\"] - b.energy) / n) ** 2).mean() + FW * ((out[\"forces\"] - b.forces) ** 2).mean()\n",
    "\n",
    "t0 = time.time()\n",
    "for epoch in range(EPOCHS):\n",
    "    mace_model.train()\n",
    "    for b in tr_loader:\n",
    "        b = b.to(DEVICE); out = mace_model(b.to_dict(), training=True, compute_force=True)\n",
    "        loss = mace_loss(out, b); opt.zero_grad(); loss.backward(); opt.step()\n",
    "    mace_model.eval(); vl = 0.0\n",
    "    for b in va_loader:\n",
    "        b = b.to(DEVICE); vl += float(mace_loss(mace_model(b.to_dict(), training=False, compute_force=True), b))\n",
    "    sched.step(vl / len(va_loader))\n",
    "print(f\"original MACE trained {EPOCHS} epochs in {time.time()-t0:.0f} s (final val {vl/len(va_loader):.3e})\")\n",
    "\n",
    "# ---------- 2.2 train xnn INDEPENDENTLY (Trainer) ----------\n",
    "# non-default flags only; everything else = stock MACE defaults (configs/model/mace.yaml)\n",
    "core = from_dict({\n",
    "    \"model\": {\"name\": \"mace\", \"cutoff\": CUTOFF, \"species\": SPECIES, \"max_L\": 1,\n",
    "              \"hidden_irreps\": HID, \"avg_num_neighbors\": LAMBDA, \"atomic_energies\": [E0[18]]},\n",
    "    \"data\": {\"batch_size\": BS},\n",
    "    \"optim\": {\"lr\": LR, \"weight_decay\": WD, \"epochs\": EPOCHS, \"energy_weight\": EW,\n",
    "              \"force_weight\": FW, \"scheduler\": \"plateau\"},\n",
    "    \"device\": DEVICE, \"seed\": 0, \"output_dir\": \"runs/argon_md_indep\",\n",
    "})\n",
    "t0 = time.time()\n",
    "trainer = Trainer(core, Subset(xnn_train, train_idx), Subset(xnn_train, val_idx))\n",
    "trainer.fit()\n",
    "print(f\"xnn (independent) trained {EPOCHS} epochs in {time.time()-t0:.0f} s\")\n",
    "xnn_indep_base = trainer.model.model     # independently trained xnn MACE"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04b8fff6",
   "metadata": {},
   "source": [
    "### Common MD utilities\n",
    "\n",
    "Switch to `float64` for the dynamics and define the ASE calculators, the NPT driver\n",
    "and the density helper used by both tracks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d2ae9b7d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:53:19.863205Z",
     "iopub.status.busy": "2026-07-20T04:53:19.863060Z",
     "iopub.status.idle": "2026-07-20T04:53:19.888790Z",
     "shell.execute_reply": "2026-07-20T04:53:19.888055Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "initial density = 1.7910 g/cm³ | target T=85.0 K, P=1.0 bar | exp ~1.41 g/cm³\n"
     ]
    }
   ],
   "source": [
    "torch.set_default_dtype(torch.float64)\n",
    "from xnn.common.models import build_model, ForceStressOutput\n",
    "from ase import Atoms\n",
    "from ase.calculators.calculator import Calculator, all_changes\n",
    "from ase.md.nptberendsen import NPTBerendsen\n",
    "from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary\n",
    "from xnn.common.deploy import XNNCalculator\n",
    "\n",
    "mace_model = mace_model.double().eval()\n",
    "xnn_indep_base = xnn_indep_base.double().eval()\n",
    "\n",
    "class MACEASECalculator(Calculator):\n",
    "    '''Minimal ASE calculator wrapping an in-memory original-MACE model.'''\n",
    "    implemented_properties = [\"energy\", \"forces\", \"stress\"]\n",
    "    def __init__(self, model, cutoff, device=\"cuda\", **kw):\n",
    "        super().__init__(**kw); self.model, self.cutoff, self.device = model, cutoff, device\n",
    "    def calculate(self, atoms=None, properties=(\"energy\",), system_changes=all_changes):\n",
    "        super().calculate(atoms, properties, system_changes)\n",
    "        conf = Configuration(atomic_numbers=np.asarray(atoms.get_atomic_numbers()),\n",
    "            positions=atoms.get_positions(), properties={}, property_weights={},\n",
    "            cell=np.asarray(atoms.get_cell()), pbc=tuple(atoms.pbc))\n",
    "        ad = AtomicData.from_config(conf, z_table=ZT, cutoff=self.cutoff)\n",
    "        b = next(iter(DL([ad], batch_size=1))).to(self.device)\n",
    "        out = self.model(b.to_dict(), training=False, compute_force=True, compute_stress=True)\n",
    "        self.results[\"energy\"] = float(out[\"energy\"].detach())\n",
    "        self.results[\"forces\"] = out[\"forces\"].detach().cpu().numpy()\n",
    "        s = out[\"stress\"][0].detach().cpu().numpy()\n",
    "        self.results[\"stress\"] = np.array([s[0,0], s[1,1], s[2,2], s[1,2], s[0,2], s[0,1]])\n",
    "\n",
    "T_K, P_BAR, DT = 85.0, 1.0, 5 * u.fs\n",
    "N_EQUIL, N_PROD = 300, 700\n",
    "AMU_A3_TO_G_CM3 = 1.6605390666\n",
    "a0 = train_structs[0]                     # dense initial configuration (400 atoms)\n",
    "\n",
    "def density(atoms):\n",
    "    return atoms.get_masses().sum() / atoms.get_volume() * AMU_A3_TO_G_CM3\n",
    "\n",
    "def compare_calcs(make_x, make_m):\n",
    "    at = Atoms(numbers=a0[\"atomic_numbers\"], positions=a0[\"pos\"], cell=a0[\"cell\"], pbc=True)\n",
    "    ax = at.copy(); ax.calc = make_x(); am = at.copy(); am.calc = make_m()\n",
    "    return (abs(ax.get_potential_energy() - am.get_potential_energy()),\n",
    "            np.abs(ax.get_forces() - am.get_forces()).max(),\n",
    "            np.abs(ax.get_stress() - am.get_stress()).max())\n",
    "\n",
    "def run_npt(make_calc, label):\n",
    "    at = Atoms(numbers=a0[\"atomic_numbers\"], positions=a0[\"pos\"], cell=a0[\"cell\"], pbc=True)\n",
    "    at.calc = make_calc()\n",
    "    MaxwellBoltzmannDistribution(at, temperature_K=T_K, rng=np.random.default_rng(0)); Stationary(at)\n",
    "    dyn = NPTBerendsen(at, timestep=DT, temperature_K=T_K, pressure_au=P_BAR * u.bar,\n",
    "                       taut=100 * u.fs, taup=1000 * u.fs, compressibility_au=2e-4 / u.bar)\n",
    "    rho = np.empty(N_EQUIL + N_PROD); temp = np.empty_like(rho)\n",
    "    t0 = time.time()\n",
    "    for k in range(N_EQUIL + N_PROD):\n",
    "        dyn.run(1); rho[k] = density(at); temp[k] = at.get_temperature()\n",
    "    print(f\"{label}: {N_EQUIL+N_PROD} steps in {time.time()-t0:.0f} s | rho_eq = {rho[N_EQUIL:].mean():.4f} g/cm³\")\n",
    "    return rho, temp\n",
    "\n",
    "rho0 = density(Atoms(numbers=a0[\"atomic_numbers\"], positions=a0[\"pos\"], cell=a0[\"cell\"], pbc=True))\n",
    "RHO_EXP = 1.41\n",
    "print(f\"initial density = {rho0:.4f} g/cm³ | target T={T_K} K, P={P_BAR} bar | exp ~{RHO_EXP} g/cm³\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f571b265",
   "metadata": {},
   "source": [
    "# Track (a): same potential (weights copied MACE → xnn)\n",
    "\n",
    "## 3a. Copy the trained original-MACE weights into `xnn`\n",
    "\n",
    "Every learnable weight of the trained original MACE, and the reference energy\n",
    "$E_0$, is copied into a fresh `xnn` MACE, so both codes carry the **identical**\n",
    "potential. We then confirm the two ASE calculators return the same energy, forces\n",
    "and stress to machine precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "93309f59",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:53:19.890505Z",
     "iopub.status.busy": "2026-07-20T04:53:19.890385Z",
     "iopub.status.idle": "2026-07-20T04:53:24.577764Z",
     "shell.execute_reply": "2026-07-20T04:53:24.576631Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n",
      "/D3/sina/xnn/.venv/lib/python3.13/site-packages/torch/jit/_check.py:178: UserWarning: The TorchScript type system doesn't support instance-level annotations on empty non-base types in `__init__`. Instead, either 1) use a type annotation in the class body, or 2) wrap the type in `torch.jit.Attribute`.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "track (a) calculators on one Argon config (SAME potential):\n",
      "  dE = 6.30e-06 eV | dF = 3.00e-07 eV/Å | dσ = 3.27e-09 eV/ų  -> identical\n"
     ]
    }
   ],
   "source": [
    "def copy_mace_into_xnn(xbase, mmod, T=2, corr=3):\n",
    "    with torch.no_grad():\n",
    "        xbase.node_embedding.load_state_dict(mmod.node_embedding.linear.state_dict())\n",
    "        for i in range(T):\n",
    "            xbase.interactions[i].load_state_dict(mmod.interactions[i].state_dict())\n",
    "            xsc = xbase.products[i].symmetric_contractions\n",
    "            msc = mmod.products[i].symmetric_contractions\n",
    "            for c in range(len(xsc.contractions)):\n",
    "                xc, mc = xsc.contractions[c], msc.contractions[c]\n",
    "                xc.weights[corr-1].copy_(mc.weights_max)\n",
    "                for nu in range(1, corr):\n",
    "                    xc.weights[nu-1].copy_(mc.weights[corr-1-nu])\n",
    "            xbase.products[i].linear.load_state_dict(mmod.products[i].linear.state_dict())\n",
    "            xr, mr = xbase.readouts[i], mmod.readouts[i]\n",
    "            if \"NonLinear\" in type(mr).__name__:\n",
    "                xr.linear_1.load_state_dict(mr.linear_1.state_dict())\n",
    "                xr.linear_2.load_state_dict(mr.linear_2.state_dict())\n",
    "            else:\n",
    "                xr.linear.load_state_dict(mr.linear.state_dict())\n",
    "        xbase.atom_ref.weight[SPECIES[0]] = float(mmod.atomic_energies_fn.atomic_energies[0])\n",
    "\n",
    "xnn_shared_base = build_model(core.model)          # fresh xnn model (float64)\n",
    "copy_mace_into_xnn(xnn_shared_base, mace_model)\n",
    "xnn_shared = ForceStressOutput(xnn_shared_base, compute_forces=True, compute_stress=True).to(DEVICE).double().eval()\n",
    "\n",
    "def xnn_a():  return XNNCalculator(xnn_shared, cutoff=CUTOFF, device=DEVICE)\n",
    "def mace_a():  return MACEASECalculator(mace_model, CUTOFF, DEVICE)\n",
    "\n",
    "dE, dF, dS = compare_calcs(xnn_a, mace_a)\n",
    "print(\"track (a) calculators on one Argon config (SAME potential):\")\n",
    "print(f\"  dE = {dE:.2e} eV | dF = {dF:.2e} eV/Å | dσ = {dS:.2e} eV/ų  -> identical\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d27f888",
   "metadata": {},
   "source": [
    "## 4a. NPT MD: same potential through both codes\n",
    "\n",
    "Same initial positions and velocities; the only difference is the calculator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ea3e8c64",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:53:24.579298Z",
     "iopub.status.busy": "2026-07-20T04:53:24.579147Z",
     "iopub.status.idle": "2026-07-20T04:59:00.662090Z",
     "shell.execute_reply": "2026-07-20T04:59:00.661311Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1051448/3511752848.py:48: DeprecationWarning: Use thermalize_momenta\n",
      "  MaxwellBoltzmannDistribution(at, temperature_K=T_K, rng=np.random.default_rng(0)); Stationary(at)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xnn (a): 1000 steps in 168 s | rho_eq = 1.4293 g/cm³\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mace (a): 1000 steps in 168 s | rho_eq = 1.4293 g/cm³\n"
     ]
    }
   ],
   "source": [
    "rho_xa, T_xa = run_npt(xnn_a, \"xnn (a)\")\n",
    "rho_ma, T_ma = run_npt(mace_a, \"mace (a)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f67024c7",
   "metadata": {},
   "source": [
    "## 5a. Result (a): the densities are identical\n",
    "\n",
    "Because both calculators evaluate the same PES, the equilibrium densities agree to\n",
    "numerical noise; the trajectories overlap until chaotic float divergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5fded164",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:59:00.664182Z",
     "iopub.status.busy": "2026-07-20T04:59:00.663914Z",
     "iopub.status.idle": "2026-07-20T04:59:00.667950Z",
     "shell.execute_reply": "2026-07-20T04:59:00.667210Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "track (a)  rho_xnn = 1.4293   rho_mace = 1.4293   |diff| = 4.85e-08 g/cm³\n"
     ]
    }
   ],
   "source": [
    "da_x, da_m = rho_xa[N_EQUIL:].mean(), rho_ma[N_EQUIL:].mean()\n",
    "print(f\"track (a)  rho_xnn = {da_x:.4f}   rho_mace = {da_m:.4f}   |diff| = {abs(da_x-da_m):.2e} g/cm³\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5072acf2",
   "metadata": {},
   "source": [
    "# Track (b): independently trained models (no weight copying)\n",
    "\n",
    "## 3b. Two independent potentials\n",
    "\n",
    "Here `xnn` is the model trained from scratch in Section 2.2 (**never** copied from\n",
    "MACE), and MACE is its independently trained counterpart. The two calculators now\n",
    "differ at the level of independent training (small, not machine precision)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d79fce3d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:59:00.669915Z",
     "iopub.status.busy": "2026-07-20T04:59:00.669794Z",
     "iopub.status.idle": "2026-07-20T04:59:00.967934Z",
     "shell.execute_reply": "2026-07-20T04:59:00.967111Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "track (b) calculators on one Argon config (INDEPENDENT models):\n",
      "  dE = 3.966e+00 eV | dF = 1.885e-02 eV/Å | dσ = 9.862e-05 eV/ų  (training-level differences)\n"
     ]
    }
   ],
   "source": [
    "xnn_indep = ForceStressOutput(xnn_indep_base, compute_forces=True, compute_stress=True).to(DEVICE).double().eval()\n",
    "\n",
    "def xnn_b():  return XNNCalculator(xnn_indep, cutoff=CUTOFF, device=DEVICE)\n",
    "def mace_b():  return MACEASECalculator(mace_model, CUTOFF, DEVICE)\n",
    "\n",
    "dE, dF, dS = compare_calcs(xnn_b, mace_b)\n",
    "print(\"track (b) calculators on one Argon config (INDEPENDENT models):\")\n",
    "print(f\"  dE = {dE:.3e} eV | dF = {dF:.3e} eV/Å | dσ = {dS:.3e} eV/ų  (training-level differences)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc6f9ceb",
   "metadata": {},
   "source": [
    "## 4b. NPT MD: two independent potentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "89f1d5e6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T04:59:00.969882Z",
     "iopub.status.busy": "2026-07-20T04:59:00.969730Z",
     "iopub.status.idle": "2026-07-20T05:04:35.710227Z",
     "shell.execute_reply": "2026-07-20T05:04:35.709466Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1051448/3511752848.py:48: DeprecationWarning: Use thermalize_momenta\n",
      "  MaxwellBoltzmannDistribution(at, temperature_K=T_K, rng=np.random.default_rng(0)); Stationary(at)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xnn (b): 1000 steps in 200 s | rho_eq = 1.4463 g/cm³\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mace (b): 1000 steps in 134 s | rho_eq = 1.4293 g/cm³\n"
     ]
    }
   ],
   "source": [
    "rho_xb, T_xb = run_npt(xnn_b, \"xnn (b)\")\n",
    "rho_mb, T_mb = run_npt(mace_b, \"mace (b)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9524d9f2",
   "metadata": {},
   "source": [
    "## 5b. Result (b): two independent density predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "19e9f4fe",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T05:04:35.711905Z",
     "iopub.status.busy": "2026-07-20T05:04:35.711766Z",
     "iopub.status.idle": "2026-07-20T05:04:35.715414Z",
     "shell.execute_reply": "2026-07-20T05:04:35.714889Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "track (b)  rho_xnn = 1.4463 ± 0.005   rho_mace = 1.4293 ± 0.010\n",
      "           |diff| = 1.70e-02 g/cm³  (within thermal fluctuations; exp ~1.41)\n"
     ]
    }
   ],
   "source": [
    "db_x, db_m = rho_xb[N_EQUIL:].mean(), rho_mb[N_EQUIL:].mean()\n",
    "sb_x, sb_m = rho_xb[N_EQUIL:].std(), rho_mb[N_EQUIL:].std()\n",
    "print(f\"track (b)  rho_xnn = {db_x:.4f} ± {sb_x:.3f}   rho_mace = {db_m:.4f} ± {sb_m:.3f}\")\n",
    "print(f\"           |diff| = {abs(db_x-db_m):.2e} g/cm³  (within thermal fluctuations; exp ~{RHO_EXP})\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e6a0cca",
   "metadata": {},
   "source": [
    "## 6. Overview: both tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "8e4cb61c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-07-20T05:04:35.716692Z",
     "iopub.status.busy": "2026-07-20T05:04:35.716575Z",
     "iopub.status.idle": "2026-07-20T05:04:36.527533Z",
     "shell.execute_reply": "2026-07-20T05:04:36.526883Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKsAAAG0CAYAAAD9xqRhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAA0XVJREFUeJzs3Xd4U+X7BvA7Sdt0p3vRUvYeZW9kgyhLQED2EFAQEVRAEAVRREEQQSlbmSICAiJ7icoqSzaFFrp3k860Sc7vD37kS6S7aU+S3p/r6iXNOec9d9JCHp+85z0SQRAEEBERERERERERmQCp2AGIiIiIiIiIiIieYbOKiIiIiIiIiIhMBptVRERERERERERkMtisIiIiIiIiIiIik8FmFRERERERERERmQw2q4iIiIiIiIiIyGSwWUVERERERERERCaDzSoiIiIiIiIiIjIZbFYREREREREREZHJYLOKqAzMnz8fs2fPLtK+GRkZqFy5MiIiIso4FRVXdnY2JBIJQkNDxY5CRERUbDNmzMCiRYsAAJGRkZBIJEhNTc1z39TUVPj7+yM+Pr5MM61fvx6tW7cu03OURlnnS0xMhEQiQWxsbJmdw1SY+s+aiEwbm1VERpaUlITVq1djxowZRdrfwcEBb775JhYuXFjGyag8mk+TJ0/W/4+BKY9JRESWLSoqCj/++CPeeeedIu3v4uKC4cOHY/HixfnuU5EaLZbOFH+WI0aMwIoVK0x+TCIqH2xWERnZjz/+iPbt28PLy6vIx4wYMQLbt2+HUqksw2SmKz09HRKJBOHh4WU6pq2tLQRBQI0aNYx2nqLq1q0bZDKZUZ9jccTFxcHT0xMuLi4F7nfx4kV07NgRTk5OqFSpEj777LNSnzs+Ph6vvfYaHBwc4OnpiZkzZ0Kr1eq3X758GZ06dYKTkxM8PT0xbNiwUn+yHxUVBWdnZ+Tk5JQ2voEVK1ZAIpHA3t4ew4YNM3geRESmZMOGDejVqxcUCkWRjxk5ciQ2b96MrKysPLd7eHhAEAT4+PgYK6ZFi42NhUQiQWJiosmNaU4/y9atW8PGxqbcG2tVqlSBRCKBRCJB//79C9z33Llz+n2f/2rSpEmpMkRFRaFPnz6wt7eHt7c35s6dC0EQAAC5ubkIDg5GkyZN4ODggLp16+L7778v1fkA4P79+/Dy8oJOpyv1WM9btGgRJBIJHBwcMG7cOKOOTZaJzSoiI9u7dy969epVrGOqVq0Kf39/HD58uIxSkZgePnyIixcvYsCAAVi3bp0oGSZNmoRmzZoVuE9cXBx69eqFTp06ISoqCvv370dwcDC+++67Up37jTfeQEZGBkJDQ3HixAn8+uuv+pliOp0Or776Kho0aICoqCiEhIQgMjISb731VqnOefDgQfTo0QM2NjalGue/pk+fDkEQcP/+fZw8eRJxcXFGHZ+IyFhKUo80aNAA9vb2OHnyZBmlIiqeGzdu4MGDB+jZsyc2bdpUrucODw+HIAiYNGlSofu2b98egiDov7RaLapWrYo33nijVBkGDhwImUyGsLAw/P7779i0aRO++eYbAMCff/6Ja9euYfPmzYiPj8eXX36JmTNn4pdffinVOQ8cOIDevXtDKjVuq2DevHkQBAE3b97Evn378r0kmegZNquIjEir1eLSpUsICgoyeNzR0RESiQQ2NjZo2LAhfvvttxeObdq0Kc6fP5/nuDk5OZg0aRK8vLzg6uqKoUOHGsw8KWj81NRUSCQSLFmyRF+EduvWDeHh4ZgwYQJcXFzg4+OD4OBgg3M+evQI/fv3h6enJ/z8/PDOO+8gMzMzz3zPzrFy5UrUq1cPzs7OeOWVVxAVFaXf58qVK+jUqROcnZ1RrVo1LFq0SD8r5dlMp6pVq0IikejX+yoow7Nz/vDDD2jYsCEcHR3RsWNHPHr0KN8x87oMsCg/m9Jau3YtBgwYgClTpmDjxo3QaDQG23U6HebNmwd3d3c4Ojqid+/eBq9daf3444/Izs7GmDFjCtzvxIkTEAQBCxYsgLOzM5o1a4YpU6Zg1apV+n2K83sBQN+gWrp0KXx9fdGoUSPMmjULa9asAfB01lVcXBzeeustODs7o3Llyhg+fDiuX7+e75h//fUXmjRpAmdnZ7Rr1w7z5s1DgwYNDPY5cOAA+vTpA+Bps3DAgAHw9vZGQECA/nKAkv7dmDx5MgICAtCrVy/4+voW+JoSEYkhMzMTN27ceKEeAYBt27ahTp06UCgUGDBgwAszWQuqR/576dizdbDWrVuH+vXrw9HREV27dkVkZKT+mNOnT6Nx48ZQKBTo3LkzHjx4YDDmvXv38Morr8Dd3R3+/v54//33oVarDcZfvXp1vpmLcnxZ58tv/Ge1iKenJyQSSZ6X9L/66quYM2eOwWMrVqxA48aN8/wZ5DXmsxxr1qxB/fr1YW1trZ9hLpFIIJfLERQUhCNHjpTqZ1nQa1GU17IkgoODMXToUEycOBHr1q3Tzyp6RqPR4IMPPoCbmxucnJzQv3//Ml93rSiOHDmCqKgog9qrsNfvv65du4YLFy5g+fLl8Pb2RvPmzfHee+/pa6guXbrghx9+QOPGjeHg4IB+/fqhffv2+PPPP/Md88SJE2jUqBGcnZ3RqVMnfPjhhy+sK/Z8DXXnzh28+uqr8PT0RGBgoP7cz35fli5dirp168Le3h6vvPIKHj9+jFGjRuln6G/ZssVg7BEjRqBatWoYOnRoobP9iSAQkdHExcUJAIQHDx7kuT07O1vYt2+f4OTkJISFhRlsmzZtmvD666/neVxwcLBQu3Zt4dGjR4JKpRJ2794t/PDDD0UaPyUlRQAg9OzZU4iIiBBiYmKExo0bC87OzsK6deuEtLQ04eDBg4KNjY0QFRUlCIIgpKenC5UrVxZmz54tJCUlCREREUKXLl2EadOm5Znv2TmaNm0qPHjwQIiLixMGDhwotG3bVhAEQUhMTBTc3NyEjz/+WEhNTRUuXLgg+Pv7C1999ZUgCIKQlpYmADB4TQrL8Oyc3bt3F8LDw4Xk5GShT58+wuDBg/MdMysrK9+fT16vXUH752fSpEnCZ599pv8+JydH8PLyEk6cOCHodDohMDBQ+PXXXw2OWbx4cb6vXV5jBgYGCgDy/Tp27Jh+34iICCEgIEB4/PixsGPHDkGhUOSbfevWrYKrq6ug0+n0j33xxRcCACE9Pb3YvxeCIAi7d+8W5HK5wWMXL14UAAhxcXGCTqcT2rVrJ0ybNk1QKpVCRESE0LFjR2HOnDl5jpeQkCA4OzsL33zzjaBSqYSzZ88KHh4eQv369fX7ZGZmCg4ODkJ8fLygUqkEf39/4a233hKio6OFyMhIYfr06UJsbGyJ/m48ExMTI9StW1c4d+5cvs+diEgsDx8+FAAIMTEx+sciIiIEAELbtm2FsLAwITo6Wujdu7fQvXt3g2PHjRsnjBs3Ls9xExISDMZ9NmafPn2EiIgIITExUejWrZswduxYQRAEITo6WrC3txe+++47QaVSCUePHhWcnZ2FVq1aCYLw9L3c29tbWLhwoZCSkiKEhYUJbdu2FT766KMiZS7q8WWdL7/xY2JiBABCQkJCvq/h77//Lnh5eQlqtVq/T506dYRVq1bl+TPIa8xnObp16yY8fvz4hWMyMzOFn3/+WXB2dhaio6NL9LMs7LUo7LUsiuHDhwvLly/Xf5+RkSEoFArhwoULQm5uruDl5SUcOXLE4Ji5c+cW+Dv93zHd3d0LrKH++eefF3JNmjRJ6NevX5GfhyAIQt++fYVhw4bpvy/s9cvL5s2bBXd3d4PHTp48qa/L/kulUgnu7u7Chg0b8hwvMjJSsLOzE9asWSOoVCrh+PHjgkKhMPgZpaSkCPb29oJKpRISExMFLy8vYebMmUJcXJzw5MkTYfLkyUJ6err+96V///5CdHS0EBERIdSuXVtQKBTC1q1bhfT0dOHnn38W7OzshKSkJIMcERERQtWqVYVr164V6bWkiovNKiIjKqxZ9cyAAQNeaDYV1Kz68ccfhcaNGxsUJkUd/9n/kF+5ckW//ZNPPhFatmxpcEylSpWEw4cPC4IgCFu2bBFq1aplsP3SpUuCm5tbnud7do7nC4jY2FgBgHDv3j1hw4YNQtWqVQ2aIN99953+HHk1lgrL8OycISEh+u179uwRatSoke+YRWk+Pf/aGaNZ9fPPPwuVK1fWP/ePP/5Y6NGjh8Ex/v7++b52eY1ZHD169NA/n8KaVVFRUYKTk5Pw6aefCiqVSrhy5YpQuXJlfTFb3N8LQRCEjRs3Cl5eXgaP3b9/3+B1vX37tlC1alV9odi5c2chLS0tz/HWrVtn0JgShKev6fOP7d+/X2jTpo0gCE//7lSuXFnQaDQvjFWSvxujR48WAAiOjo7CgAEDhIyMjHyfOxGRWApqVv3111/6x0JDQwUABg35kjSrnr1fCYIg/PTTT0JQUJAgCIKwYsUKoUWLFgZjTJ8+Xf8/xz/88IPQpEkTg+2nTp0SKleuXKTMRT2+rPPlN35RmlVarVaoVq2asHPnTv34dnZ2QkpKipCXgppVhf3Pf/fu3YUtW7bkmaOw51LYa1HYa1kU/20sbdiwQahTp47++/fee0947bXXDI5xdXUt8Hf6v2OWRHGbVREREYJMJhPOnDmjf6yw1y8vK1euFKpVq2bw2JUrV174uy0ITz8cfeWVV4SXXnpJyM3NzXO85cuXv/DzeO+99wwe2759u77Zt3r1aqFOnToG9fvzzxGAcP/+ff1jM2fOFLp06WKwn0Kh0H+wN3DgQAGA4OTkJAwdOlTIzs7O97kTCYIg8DJAIiNyd3eHXC5HQkKC/jHh/y+rql69OuRyOSQSCfbu3WswrRp4ejmUn59fnuOOGDECo0ePRv/+/dGmTRvMmjVLf46iju/p6an/s62trcH3zx57tqBqWFgY7t+/b7BAZIsWLZCcnIy0tLR8n3+1atX0f/b29oaDgwMiIyMRERGBatWqQSKR6LfXqFEDERER+Y5V1AzPL2RvZ2eX76KweSnqa1cawcHBGDlypP65jx49GsePH0dYWBgAQK1WIzIyMt/XrjTWrVuH3NzcIq21AAB+fn7Yv38/fv/9d3h7e2PYsGEYN24cZDIZXFxcCv2ZTJgwQf/4s0tPnJ2dX7hxwLM1CpydnREdHY22bdti6tSpSEtLQ0xMDHx9fdGtW7c8M0ZFRaFq1aoGjz3/2gGG09fDw8NRo0YNyGSyfJ93cf5ubN68GYIgIC0tDXv27IG9vX2+4xIRicXHxwdSqdSgHnnm+X8zn10q//z7TUH1SH7yey/+7/sbAFSvXl3/57CwMFy9etXgfaVz586IiIgwWNw5v8xFPb6s85WmFpFKpZg0aZL+8qo1a9Zg8ODBJbpEqnLlyvo/a7VafPTRR6hatSpsbGwgkUhw7NixQmuL/J5LYa9FYa9lSQQHB2P06NH678eMGYP9+/frL11MSkpCSkpKob/T5W3dunWoVasWOnbsqH+ssNdv6NCh+sfbt28PoOAaysnJSf+YWq3G4MGDoVQqceDAAVhZWeWZqyQ1VK1atQzq9/8qTg21e/duCIIAlUqFHTt2QC6X5zsuEcA1q4iMSiaToXnz5rh27Zr+sV9//RXr16/H7t27kZqaCkEQ0K9fvxfWLbpy5QpatWqV57hSqRTvvfcezp07h5MnTyI9PV2/YGNRxy+OwMBANGrUyGChyGdfz785/tezBgzwtNjNyMiAv78/AgICEBYWZrDOQGhoKAICAvTPz1gZninKopBl8do9LzQ0FKdOncKoUaP0j1WvXh1t27bVL7T+fBH4zPOvXV6evztNXl/Hjx8HAFy4cAGnTp2CVCqFRCLBsGHDoFQqIZFIcPDgwTzH7tSpEy5evIjMzEzcvXsXKpUKbdq0ga2tbaE/k/Xr1+u/f/Z3oHHjxlCr1fj333/157h8+TJ8fHzg5eWFv//+G1lZWZgxYwYcHR3h4+ODd999FxcuXEBSUtIL+SpVqmTwWgHQr1MGPG1A/v777/pCq0qVKggNDeVd+4ioQrG3t0fDhg0N6pFnnv839NkC0s+/3xRUjxSXv79/gf9mBwYGom3bti+8p+h0OoP38fwyF/X4ss6Xn6IuUD1+/HicP38ef/75J/bu3Ys333yzRGM+31TYsmULdu3ahf3790OlUkEQBHTt2rXENU5hr0Vhr2VxXb9+HZcvX8aIESP0jzVq1AgNGjTAxo0bAeRdQ+X1O/08Dw+PAmuo/NZrKyqtVosNGza88EFhYa/fzp079Y+dO3cOwNMaKikpyeD5Xb58GTVq1ICDgwOAp+vT9enTB+np6Thy5EiBNXJhNZRGo8Hhw4fx6quvAnhaQ92/f/+FdcKIygubVURGNmDAABw9elT/fXp6OqysrKBQKKDT6bB9+3YcOnTI4Jjw8HBERkbi5ZdfznPMr7/+GuvWrUNcXBxyc3Oh0+mQnp5e5PGLq3///khNTcW8efMQFxcHpVKJQ4cOFTpD56OPPsLDhw+RkJCAqVOnonXr1qhVqxb69euHlJQULFy4EEqlEpcuXcKSJUswfvx4AE+LamdnZ1y9elX/hljSDM/kNeZ/lcVr97y1a9eiVatWqFWrlsHjo0ePxsaNG5Gbm6t/7NNPP83ztcvLs0Isv69ns5Kebx4JgoAdO3ZAoVBAEAR9IVKjRg2DBV/fe+893LlzBykpKdi4cSPWrl2r316Sn0mNGjXQtWtXvP/++4iJicGNGzewZMkS/TGNGzeGVqvFihUrkJ6ejvj4eHz33XeoUqUK3N3dXxivf//+ePLkCZYvX460tDScO3cOP/zwg357SEgIrK2t9Quu9+/fH1qtFtOmTUNsbCxiYmIwY8YM3sWPiCzef+uRZz788EM8fvwYsbGxePfdd9G1a1f9TKpbt24hIyMDXbt2NUqGwYMH4+bNm/j++++RlpaGEydOYMOGDfrtQ4YMwcOHD/H5558jISEBqamp+O233/DOO+8UKXNRjy/rfPlxd3eHjY1NgbXIs/0GDx6MQYMGoXr16vqZNaUZ8/kaR6vVYuPGjTh16lSRcuelsNeisNeyuIKDg9GlS5cXmk6jR4/GunXrDGa2zZs3L9/f6f9KTEwssIb672LjBfHx8dHftOWZ/fv3Izk52WBGGFCy36WgoCC0atUK7733HuLj4xESEoLly5frayiVSoWePXvC1tYWv//+e6GzvQcPHowbN24gODgY6enpOHnypL7xBwDnzp2Dn5+ffvbVkCFDkJiYiA8++ADx8fGIjIzE22+/jYyMjCK/RkSlwWYVkZGNHj0aZ8+e1U+9HzZsGNq3b4+goCD4+Phg7969LzSltm7dimHDhkGhUOQ7ZkhICBo3boyAgACEhoZi/fr1RR6/uJydnXHq1Cncv38fQUFBqFq1KtasWYOpU6cWeNyIESPQp08fVK9eHWlpafpb57q7u+Po0aM4ceIE/P398frrr2PChAmYMWOG/tglS5ZgypQpkMlkmD17dokzPO+/Y/5XWbx2z+Tk5GDz5s04f/78C5/avfnmm4iLizO48+CkSZPQu3dvVK1a1eC1K2/9+vXDsGHD4Ofnh7Vr12Lfvn146aWXAJT892L79u1wcHBA9erV0blzZwwcOBAff/wxAKBmzZrYvXs3tm7dCh8fH9SpUweJiYn5zvzy8PDAwYMH8eOPP8LPzw+zZs3C6NGjYWNjAwA4ePCgflbVs8wnT55EREQEGjRogFatWiEwMBDe3t7GeLmIiEzWhAkT8Mcff0ClUhk8PnToUPTo0QO1atWCVCrF1q1b9du2bNmCMWPGwM7OzigZnl1evmbNGlSqVAmfffYZJk+erN/u7u6OM2fO4PLly2jQoAGqV6+On3766YX3lfwyF/X4ss6XH5lMhsWLF2PkyJGQyWR53g3wmbfffhvx8fEFzqoqzpjjxo1Do0aNUL9+ffj5+eHEiRP5XmJfFIW9FoW9lsWRkZGBbdu24fjx4y/UUO+99x7Cw8MNGrHjxo1D165dUbNmzRd+p0ujf//+kEgkCA4Oxm+//QaJRFJoMys4OBhDhgx54TLOkv4u/frrr9BoNAgMDMTLL7+MMWPG6OvnQ4cO4dy5czhw4ABsbW31r9GgQYPyHKtSpUrYt28fvvvuO/j5+WHBggUYNWpUvjWUu7s7Tp8+jVu3bqFOnTro0KGD/s6DROVBInBeH5HRzZ8/Hzk5Ofjyyy8L3TczMxN16tTBX3/9pb8sztykpqbC1dUVERER+U67rigmT54Mf39/zJs3r0j7F+W1K+6YFc2HH36IBw8eYO/evWjWrBkWL16MHj16iB2LiEh0M2bMgJubW5HeP1JTU9GgQQNcuXLFYN0iMUVGRiIgIAApKSkWf5v7U6dOoVevXoiOjs5zZnFFMGLECDRv3hzTp08v0v5F+f0o7pgVzdtvv420tDRs2bIFtWvXxqZNm9C2bVuxYxEBAPJefY2ISmXhwoVF3tfe3h5PnjwpwzRElmXx4sXo3LkzGjRogNOnTyM4OBjr1q1DdHQ0Hjx4gE6dOokdkYjIJHzzzTdF3tfFxUXURakrsuzsbHz11VcYPnx4hW1UUflYsGABXnnlFdSpUwdHjx7F5s2bsXv3bjx48ADJycnFugySqKyxWUVERGalV69eePfdd3HlyhX4+/vjyy+/xOuvvw4AL1zuQkREZMqOHz+OHj16ICgoCJs2bRI7Dlm4Xr16Ydq0abhx4wYqV66M77//Hr179waAPO8eSiQmXgZIREREREREREQmgwusExERERERERGRyWCzioiIiIiIiIiITAabVUREREREREREZDIq5ALrOp0O0dHRcHJygkQiETsOERERkZ4gCEhLS4Ofnx+k0uJ/rsg6h4iIiExVUeucCtmsio6ORkBAgNgxiIiIiPIVEREBf3//Yh/HOoeIiIhMXWF1ToVsVjk5OQF4+uI4OzuLnIaIiIjof1QqFQICAvT1SnGxziEiIiJTVdQ6p0I2q55NiXd2dmYRR2YnLi4OW7duxYgRI+Dt7S12HCIiKiMlvYSPdQ6ZO9Y6RESWr7A6hwusE5kZe3t7tGjRAvb29mJHISIiIjI61jpERFQhZ1YRmTMnJyd07NhR7BhEREREZYK1DhERsVlFZGZycnIQFxcHb29v2NjYiB2HiMqRVqtFbm6u2DHISKytrSGTycSOQWRyWOuQOeJ7NJGh0tY5bFYRmZmkpCRs3LgREydOhK+vr9hxiKicpKenIzIyEoIgiB2FjEQikcDf3x+Ojo5iRyEyKax1yNzwPZroRaWtc9isIjIzHh4eeOutt+Dq6ip2FCIqJ1qtFpGRkbC3t4enp2eJF94m0yEIAhISEhAZGYmaNWtyhhXRc1jrkDnhezTRi4xR57BZRWRmrK2t4eXlJXYMIipHubm5EAQBnp6esLOzEzsOGYmnpyfCw8ORm5vLZhXRc1jrkDnhezRR3kpb5/BugERmRqVS4ejRo1CpVGJHIaJyxk9rLQt/nkR5Y61D5oj/phMZKu3fCTariMxMdnY27t+/j+zsbLGjEBERERkdax0i83bgwAH88ccfpRrjgw8+QEREhJESlU5ZZPn444+RlJRk1DEtDZtVRGbGy8sLU6dO5fR4IiIiskisdYjMV25uLmbPno3WrVsDAD766CM8ePCg2OMcO3YMSqXS2PFKxBhZZs2ahUePHum/r1mzJj7//PPSRrNobFYRERERERERUant378fQUFB+hskdO/eHe7u7iKnEt+JEyeQnJys/37QoEHYtWsXMjMzRUxl2tisIjIz8fHxWLlyJeLj48WOQkQV2K5du7Bs2TL99ytXrsTOnTsBAP/++y8WLlyI/fv3Y8KECfjiiy+gVqsL3ZaTk4NvvvkGI0aMwKBBg7BmzZryf2JEJDrWOkQlp1QqMWrUKMTGxgIAwsLCMH78eP1ltR988AEuXLiAOXPmYOLEiQgJCdEfW9C2f/75B2+99RYGDx6MQYMG5dtkOXjwILp06aL//tixY/rL3QoaPzMzE5999hkmTpyIo0ePGoyp0+mwYcMGjBs3DjNmzMD9+/cNMv/111+YNWsW3n77bdy4caPIx5VFlrzG3LlzJx49eoRZs2Zh0KBBuHv3Luzt7VG3bl2cPXs2z9eReDdAIrNja2uLevXqwdbWVuwoRCSirBwtHiakl9n41T0dYWeT/51bXn31VbRs2RK1a9eGXC7HmjVrcPHiRQBAXFwcli1bhkGDBuGll17C2rVroVarsWDBggK3LV26FGfOnMHo0aNhY2ODqlWrltnzIyLTxVqHzJnY788KhQItW7bEiBEjcPDgQQwdOhTTp0/X/306duwYTpw4gUmTJsHe3h4vv/wyHj9+DDs7u3y35eTkoF+/fpg3bx78/PwAPL1rZ15u3ryJsWPH6r8/efIk+vbti5o1axZ47mHDhkGn02HgwIH4/vvvERoaqh9j4sSJAICuXbsiKioK3bt3x6VLl+Dl5YVjx47h0KFDeOeddxAbG4vOnTvj5s2b8PX1LfS4ssiS15gNGzaEq6srevbsiWrVqsHT0xMAUKdOHdy4cQO9evUq9u9BRcBmFZGZcXZ2Rrdu3cSOQUQie5iQjle/O1dm4x98pz0aVFLku93e3h67du1Cr169IJVK8dtvv8HR0VG/3d/fHxs2bAAAVKpUCUuWLCl0W1paGjw9PVGvXj00atQIUikngJsqQRCg1uhgJZXASsafExkXax0yZ2K/PwPA1KlTcebMGTRt2hTt27fHsGHDDLYvWrQIvXv3BgDs2bMHDx48QKNGjfLd5u3tDalUiipVqqBDhw76S/zykpGRAXt7+3y35zW+r68vzpw5g9jYWNja2mLo0KHw8fEBAKSmpmLLli3o06cP9u7dCwDQarU4e/YsBg0aBACYP38+hgwZAgCIiIjAzp07MXbs2EKPK4ss+b22rq6u6NKlC5o3b65/Lezt7ZGRkZHva1XRsVlVRtYfv4Hw2AQsGtFV7ChkYXJzc5GcnAw3N7d8P9EgIstX3dMRB99pX6bjF7pP9epwdHSEXC5HgwYNDLb5+vrq/+zg4GBwuUB+2+bMmYPPP/8cw4cPR1JSEr7++muMHDmytE+FykBEYjoGLfsNK8d2Quva/mLHIQvDWofMmSm8PwNA586dsXv3bixfvvyFbc9mRwEvvkfntc3b2xtr167FmjVrMHr0aHTo0AHbt283+JDqGS8vrwLvcpfX+AkJCfD29tbP/rK1tdXfYCE2NhYODg4YOnSo/rihQ4eiRYsW+u8rV66s/3NgYCASExOLdFxZZCnotf2vpKQkNG3aNN/tFR2bVWWkwf3v0C3xPIB/xY5CFiYxMRFr167FxIkTDf6Hj4gqFjsbWaGfrJa16dOno2vXrkhNTcXHH3+ML774olTjOTs7Y8mSJViyZAn+/PNPjB49ms0qE2Wfk4iLtlNwLWodUPt1seOQhWGtQ+bMFN6fb968iS+//BJ79uzB5MmTcfHiRf2lZyXVt29f9O3bF7m5uejSpQuOHj2K11577YX9WrVqhevXr6Nnz55FHrtKlSqIi4vDgwcPULNmTTx48ABhYWEAgGrVqkGn08HLywsdO3bM8/hjx46hTZs20Ol0+svwinJcWWTJj62trX6NzmeuXr2KKVOmFGucioTNqjIiWNnBWlAXviNRMbm7u2PChAm8qwYRiWrXrl24ePEi/v77b+Tm5qJFixZ46aWXilWc/teyZcvwzz//QKfT4dq1axg4cKARE5Mx2Tg8/TRdq84SOQlZItY6RCWXkZGBIUOGYM2aNejduzfCwsIwatQoHDp0CBKJpERjPnjwAHPmzAEApKSk4PHjx2jdunWe+w4bNgzvvPMOPvzwwyKPb29vj4ULF6JNmzZo3ry5fnYTANjY2GDDhg0YMGCAwV0GV61apb887/jx4zh58iQSEhLg6uqK119/HVZWVoUeVxZZ8tO2bVuMHz8eDRo0wKJFi+Dg4ICsrCw0adKkyK9TRSMRBEEQO0R5U6lUUCgUUCqVcHZ2LpNznN80CzUe74THp4/LZHwiIqo4srOzERYWhqpVq5rMgsPHjh1D3bp14e//9BKwsLAwPH78GJ06dUJ8fDzu3buHDh06AACSk5Nx48aNQredP38ekZGRsLa2Rs2aNVGvXj3Rnl95yO/nWto6pTzqnFx1JqwX++Jiky/Rst9bZXIOIiJzYGrv0Y8fP8bDhw8N7si3f/9+tG/fHm5ubjh+/DhatGgBheLp7K/Tp0+jUaNGBW4Dni6ULpFI4OrqijZt2sDOzi7fDC+//DK+/PJLNG7cGKdOnULjxo0LPTcA3Lt3D5GRkWjZsiUuXLiAli1b6t/HkpOTcfXqVaSkpOjP4eDggKCgIGzatAmCICA9PR1t2rQxuHw4v+PKIktBYwqCgL///huxsbHo1KkTvvnmG9SoUcNgMXpLU9o6h82qsmpWbf0U9R+sgdOC2DIZnyqutLQ0hISEoFmzZnBychI7DhGVA1MrhMk4zLlZBUGA7lNXXGzwMVoPnlk256AKi7UOmRO+R7/oyZMnSE9PL5cPnYKCgrB582YEBQWV+bmM6ejRo+jWrZtF30ymtHWO5b4yIpPY2MMWOWLHIAuUmZmJK1euFLhYHxERUZmSSJANGwg5fC8i42OtQ2TeKleuXG6zo5cuXYqqVauWy7mMqUePHhbdqDIGrllVRiKrDsLIK7VwS6uDNW/pTEbk7e2NGTNmiB2DiIgquJ7WGzHYuybaiB2ELA5rHSIqqm7duokdgcoIuyhlxFZuhxxYIytXK3YUIiIiIqOT2DggS1PhVpMgIiKicsBmVRnxSvsXO20+g1oZL3YUsjAJCQn44YcfkJCQIHYUIiKqwD7KXY2GkTvFjkEWiLUOERGxWVVG7IRstJbeQU6GSuwoZGFsbGxQpUoV2NjYiB2FiIgqsBpCGNwyH4kdgywQax0iIuKaVWXEytYBAJCTnS5yErI0CoUCL7/8stgxiIiogtNI5ZBqssWOQRaItQ4REXFmVRmxsXUEAORmZ4ichCyNRqNBcnIyNBqN2FGIiExSdHQ0fvnlF7FjWDyN1BYybZbYMcgCsdYhslw///wzYmNjxY5RZKdPn0ZYWFipxli5cqWR0pSesbNoNBps2bLFqGM+w2ZVGbGxtQcA5GaxWUXGlZCQgO+++47rOBCRWSlqA8kYjaZHjx7h22+/zXNbQkICVqxYgb///tvg8bNnz2LFihVISUnRPxYXF4cVK1bk++/t3bt3sXHjRmzbtg2PHz/WP3779m2sWLHC4Gv//v2lek6mSCuzhUzLmVVkfKx1iMpXURtIxmg0LV++HOHh4Xlu27NnD1atWgVB+N/NOzIyMrBixQocP37cYN/du3fj8OHDeY6TkZGB3377DcHBwTh9+jR0Op1+W3Bw8Avv0c+/9z8vLS0N06ZNg4+PDwBg1apV0GqLfwM1U7q7qTGyrFy5Uv8zsrKywsGDB3HmzJlSj/tfbFaVERsXX8zKfROp9gFiRyEL4+bmhtGjR8PNzU3sKERERaZWqxETE2O0/UoqKioKH3/8MaZOnWrw+JQpU/DRRx8hLi5O/9i6deswf/58bNy48YVxpk2bhrZt2+Lw4cM4cuQIevbsifXr1wMALl68iBUrViA8PFz/9fy4luIf9wE4aNdX7BhkgVjrEJWvmJgYqNVqo+1XUt9//z3mzp2LEydO6B/bvn07PvnkE+zc+b8beqSmpmLy5MmYMGHCCzMwL168iBo1amDJkiW4fPkyFi1ahPbt2yMr6+lM4Llz5yIkJMTgPTo3NzfPPD/++CN69eoFOzs7AMDjx48NGmkV1YwZMwyadhMnTsTSpUuNfh6uWVVG5A4K/KztjE5WnmJHIQsjl8tRpUoVsWMQESEtLQ2HDx9GdnY2unTpgkqVKgEAIiIicOnSJTRp0gSnT59G9+7dIZfL4evrqz82ISEBhw8fhqenJ+rWrYurV6+if//+Bvs9G6dZs2Y4deoUqlWrho4dOwIAkpKSsGXLFkgkEnh6eqJnz55wd3cvUm5vb294enri/PnzaN26NU6dOoXatWu/MItj8+bNWLt2LT799FPMmjVL//iWLVvw66+/4ubNm/Dz8wMA5Obm4tKlS/p96tSpgxUrVhT/RTUjEa6tcT0jVewYZIFY6xCV3sWLF3Hjxg3UqFEDnTp10j++bds2dO7cGdeuXYNWq0WfPn3g6+sLuVyu3+fkyZOIjIxE165dcfr0afTq1Qvu7u4G+23btg1dunTBlStXkJSUhL59+8LFxQUAsG/fPoSHh8PW1hZBQUFo3bp1kXOPGDECa9asQbdu3QA8/eBo+PDhyMnJ0e+zY8cO9OnTB9HR0fjjjz/Qp08fAEBOTg4GDRqE999/HzNnztTvf+nSJYMm06xZs9CgQYNCs+zevRuffvqp/vvAwEBIJBIAwIoVKzBu3DgcPnwYEokE/fv3h7W1NYCnl8YdPHgQarUavXr1emHcv//+G7du3UKVKlXQrVs3gzHHjBmDI0eOwNbWFq+++ipkMlmRjiuLLHmNuX//fgiCgJUrV0IqlWLKlCno1KkTBg4ciLS0NDg5ORX6uhYVZ1aVEVsrYLDsNGTJD8SOQhYmPT0df/31F9LTuXg/EYknKSkJjRo1woYNG/DHH3+gUaNGuHz5MgDg3r17mD59Ot544w2EhIQgOzvb4NK8+Ph4BAUF4ZdffsHGjRvRp08frFmzBoDhJXz37t3De++9h+HDh+Pvv//GsGHD9LOccnNzER4ejrCwMBw4cABBQUHFujThrbfeQnBwMICnlwRMmjTJYPuZM2fg5uaGoUOHwsHBAX/99Zd+286dOzFlyhR9owoArK2t0bZtW/33jx8/NrjE4PnjLUW1nHtokXFW7BhkgVjrEJXOwoULMWjQIJw/fx5vvvkmJk6cqN/29ddfo0+fPtixY4f+ffP5S/Nmz56NCRMm4OzZsxgwYABmzJihn/H8/H5ff/01+vXrh507d2L79u3o2LGjfrZNXFwcwsPDcePGDYwbNw5LliwpcvZu3brhzp07iI2NxeXLl+Hq6opq1aoZ7LNx40aMHTsWY8aMMZj9/NdffyErKwvvvfeewf4tWrSAvb29/vutW7fq358LWsPp8uXLaNiwof77999/Xz8L67333kP//v1x9OhRLFmyBKNHj9bv179/fyxcuBBHjhzBq6++ajDm5MmTMWfOHFy9ehULFizA66+/rt/23nvvoU+fPjh69Cg+/fRTDBw4sMjHlUWWvMaMi4uDIAh4/PgxwsPDIQgCZDIZateujStXruT7WpYEZ1aVERuZDF9arcPleG8AHcWOQxYkPT0d586dQ/Xq1eHo6Ch2HCISU1rs06/n2bkArlWA3Gwg4e6Lx/gFPf1v4gMg5z/rKrpUBuyLdtnNmjVr0K5dO2zduhUAsHTpUnzxxRfYs2cPgKeX8508eVI/df75RlJwcDBeeeUVrF27FgAwb948faPrv7RaLU6cOAG5XI59+/Zhw4YNGDduHHx8fDB9+nScPXsWKSkpiIyMxO+//47x48cXKX+fPn3w4Ycf4v79+/j333/1n+A+s2nTJowZMwYAMHr0aGzcuBHt2rUD8HRdrf8Wzv+VlZVlsCZHjRo1ipTLnDRUnUFP9VEAc8WOQhaGtQ6ZPRHfn9PT07FkyRLcvXsXAQEBSElJQbVq1TBr1ixUr14dADBp0iSDBtYzKpUK33//PUJDQ+Hl5YXY2NgCZzm+/fbb+vfKmjVrIjQ0FLVr18b48eNx9OhRhIWFwcPDA5s3bzaYoVwQiUSCsWPHYsOGDQgLC8PkyZPx8OFD/fabN28iOTkZHTp0QHZ2NqZNm4b4+Hh4eXkhOjoagYGBkEoLnpMTFRWF7Oynay4+P3PpeWq1GhkZGVAoFPmO88033yAoKAgqlQp+fn4QBAEXL17E3bt3cffuXVhZWeHUqVPo3r07AOD69es4ePAg3n//ff1r9umnnyIqKko/O/3zzz9Hx44dkZubi1q1aiEkJES/LlRBx5VFlrzGfPPNN/HWW29h2bJlsLL6XzvJ1dUVSUlJBb7uxcVmVRmRSKVQwwZCTqbYUcjC+Pj4FPkfeyKycJc3AWe+NHys4evAwHWAKgpY+9KLx3yqfPrffW8BkZcMtw1YCzQeUqRTP3z40GBaf7t27fDTTz/pv69Xr56+UfVfjx490jd+AKB169b5Nqvq1Kmjv+TAx8dHvwjqmTNnMHDgQP3lf9nZ2UhMTCxSduBpcTp8+HD069cP48aN0097B55e3rh7924EBgZixYoVSEpKwq5du7By5Uo4ODjA3d0d0dHRBY5fES4DhJUd5ELZrV1CFRdrHTJ7Ir4/P3nyBJ6enggIeLp2squrK+rVq4eHDx/qm1XNmjXL89iIiAj4+fnBy8sLwNO/iwU1q4KCgvR/9vb21r9Hv/TSS5BIJGjYsCEEQSjW+zMAjB07Fm3atAHw9MOx599PN27ciEqVKulnYfv6+mLLli2YOXNmkd6fgaJdBiiXy2Fra4v09HT95Y3/9ez5Ozs7A3ja4AoNDUXz5s31jZzna6Xbt2/D1tbW4MOssWPHGozZsmVLAE9nbDdt2hQPHz6EVqst9LiyyJLXmLa2tnm+Fmlpafm+TiVV7s2qnJwc/fWmcrlcfy1lXjIzMw1W7gcAOzs7yGQyqNVqg4XQrKys8n3hxJItkUPI5S2diYiojDQfC9R+2fAxO5en/3WuBEws4M4s/X/I+5PbIqpUqRJu376t//7mzZvw9/fXf5/fJ5UA4Ofnh3v37um/f36cotq3bx9mz56t/0SwT58+xV70dOLEiUhJSXmhONu1axeqVKkCpVIJpfLp/zwEBgZi165dGDt2LPr27YsffvgBkydPNris4MGDB6hZs2axn4vZsraDjZBT+H5ERBWNiO/Pfn5+iI+PR3JyMtzc3JCdnY3Q0NAivUf7+voiOjoa6enpcHR0RFpaGiIjI4t8buDprKVHjx7pLx08c+YMdu/eXawx3N3dMWPGDLi5uRnM3snNzcW2bdvw2muv6ZssTZo0wcaNGzFz5kz9bKtdu3YZXNL25MkTeHl5Fbtf0LhxY9y5c0ffOCsKf39/g7rm+T8HBgYiNzcXX375Zb5Z7t69i6CgIAiCgLt376JSpUqQSCSFHlcWWfJjZWVl0Kd5lvX55qUxlHuzas6cOQgODoZarcbXX3+N6dOn57tvixYtDG4FnZGRgQsXLqBly5Z49913sXnzZv0v74ABA7Bly5ayjl8saokcYLOKjCwxMRG//fYb+vXrBw8PD7HjEJGYnHyefuXF2vZ/lxTkxaN0TZWJEyeiSZMmEAQBrq6uCA4ONrhTT0HGjRuHFi1aIDc3F3K5HHv27EGtWrWKdf5GjRph0aJFyMnJwb///osLFy4YzNYqCl9f3zxnP23cuBFz587FsGHD9I/t3LkTq1evxtixY/HWW2/h6NGjaNiwIQYNGgRbW1ucOnUKL7/8MubMmQPgf2tWPePt7W0wniWQ2tjDFmoIgmAwM42otFjrkNkT8f3ZxcUFI0eORLdu3TBw4EAcPXoUzZs3R7169Qo91s3NDf369UOPHj3wyiuv4MSJE7C2ti7wA6j/8vT0hEQiwVtvvQUvLy/s3bu3RM9j8uTJLzx28OBB+Pv744cfftA/JggCAgMDceHCBbRq1QqbN2/GG2+8gf3796NevXp49OgRzp8/jwsXLuiP2bp1K3x8/vfzGTRokEEz75n+/fvj2LFjxWpWdejQAXK5HH379kWLFi1w+PBh/ba2bduicePG6NChAwYMGAB7e3vY2Njg7bff1u/zzjvvoEePHrh48SKcnJzQtm1bSCSSQo8riyz5qVmzJt5//31Uq1YNU6ZMwdWrV9GwYUOj38G13BdYX7ZsGdLT0zFkSOHTGG/duoX09HSkp6fj999/R82aNfXT4oCnK9Q/225qjSoAuGVVH8kyvsGScVlZWcHT09PgUwYiovIWEBCAq1evomrVqrC1tcXx48f16yBUrlzZYFFQ4OlMrGefclavXh1///033NzcULVqVbz99tv69RGe3++/4/j6+uq3jRkzBp9//jkyMzMxePBgbNq0CS1atHhhjP/y9PR8YSbVMxMnToSdnR1atmyJfv36GWzr378/WrRogaysLFhbW+PAgQNYvXo1bGxsYG1tjS+//FLfqKpXrx569uxpcFvsZ58wW5JcJ3+E6GpBnaspfGeiYmCtQ1Q6a9aswezZs5GZmYmxY8fit99+028bPnw4PD0N71g/dOhQ/Z14N23ahPHjx0Oj0eCLL76ATCbT31Dk+f3+O87rr78OX19f2NjY4OzZs/D29oaLiwt+/fVXg8bT82P818CBA/WXKj6vWbNm6N69O9RqNRYuXGiwTSKRYOnSpfr1kvr27Ys7d+6gWbNmUKlUaNeuHS5cuAAHBwcAT5tg2dnZBu/Rz9av+q/x48fj119/hUbz9H3unXfe0f+79O677xrsO2XKFFhZWUEqleLkyZPo0qULbGxssGPHDoN9n80Mz8jIQHh4uMHkHODpHYetra3RtWtXHD9+XP9hUEHHlUWW/MYEnt4l0c3NTb/A+oYNGwqchFRSEqG4c+aNZMSIEWjevHmRn9TIkSNRv359zJ49G8DTX7KgoCBMmDCh2G9kKpUKCoUCSqVSf/1lWejz3Tk09FfgiwENC9+ZiIgoH9nZ2QgLC9M3hizBDz/8ALVajcjISGzatAkHDhwwuJteRZDfz7W0dUp51TmHb8Zi8tYQXJvfHS72NmV2HiIiU2Zp79EXL17E33//jezsbBw8eBB+fn7YtWuX2LFE8/PPP6NJkybFngFeEhKJpNhLGohNq9Xiq6++0n9g97zS1jnlPrOqJFJTU/Hbb78Z3ILR1tYW8+bNg4ODA+rWrYs//vgj3+PVajVUKpXBV3lwsAI02VxgnYxLq9UiLS1Nf2tYIiJz9OTJEzx+/Biurq44efJkhWtUGZNYdY6dtRRy5CArJ7fwnYmKgbUOkXiUSiXCw8ORkpKCiRMnYvv27WJHEtWQIUPKpVEFvDibyRzIZLI8G1XGYBZza7dt24bOnTsbTBdcsWIFVqxYAa1Wi19++QVDhgzBkydP8lyBfvHixViwYEE5Jn7qI+Un0KTbAzhQ7ucmyxUfH4+1a9di4sSJ+U6hJSIydYsXLxY7gsUQq87xTrqAe7Zj8CTpH8Cl8LVQiIqKtQ6ReLp3766/rJ/Kl8XfRbiYzGJm1YYNGzB+/Pg8t8lkMgwdOhQuLi54+PBhnvvMmTNHf0cfpVKJiIiIsoyrp5XZQqblAutkXK6urhg2bBhcXV3FjkJERCZArDrHWv50/Y+crPRyOR9VHKx1iIio3GdW5ebmQq1WQ6PRICcnB+np6bC3t4dUKoVarYZEIoGNzf/WPQgJCUFMTAx69+5tME5GRgYEQUBGRgZ27NiBzMxM1K5dO89zyuVyyOXyMn1eedHKbCHPVZb7ecmy2dralttUVCIiMn1i1TlWtvYAgNzsjEL2JCoe1jpERFTuM6s2bdoEHx8fHDx4EAsXLoSPjw+uX78OAJg5cyY+//xzg/23b9+O8ePHv7CIeu3ateHj44NGjRrh999/x6FDh+Do6Fhuz6ModFa2sNKpxY5BFiYjIwMXL15ERgb/54CIiMRjY/v/M6u4PicZGWsdIiIq95lVEydOxMSJE/PctmrVqhceW7ZsWZ77RkZGGjVXWdBZ2cNGx8sAybhUKhWOHj2KgIAA/S1YiYiIyputw9M7+GiyeRkgGRdrHSIiMosF1s3VpcCJ+EP5Cg6JHYQsiq+vL+bNmyd2DCIiquDsXH3RKnsVPnJpgWZihyGLwlqHiIjYrCpDUntXxGpSxY5BRERUYahUKty+fRsA4OzsjHr1inaXutDQUKjVatSvX79U41QkchtrJErdkaaRiR2FiIjMwM2bN5Ge/nQ2boMGDYq0jE9GRgb+/fdf1K9fH05OTiUeh8yPWdwN0FzVTLuAL3K/EjsGWZikpCT89NNPSEpKEjsKEZHJCQ0NxfTp0zF69Oh8lx34rydPnqBdu3YGdx4uyTgVjUQiwUqb1fCMOiF2FLIwrHWILNOyZcswffp0dOrUCTdv3izSMc/2f7bOdUnHIfPDZlUZctEmoZf0IjS5uWJHIQsilUrh4OAAqZR/fYlIfAkJCbh58yaysv63RuO1a9eQkpKi/z40NBRRUVEAgJiYGDx58gRqtRp37txBTk6OUfM0bdoU58+fx5IlS4q0vyAIeOutt/Duu++WapyKqh2uw0EVKnYMsjCsdYhKLycnB7dv30Z0dLT+sfj4eP2sYeDprKVLly7pvz9//jwA4PHjx4iNjTV6pk2bNuH8+fPw9/cv0v6HDx8GAFSpUqVU45B54jtAGbKSP52OmJmZJnISsiSurq4YOHAgXF1dxY5CRBWYTqfDuHHj0LRpU4wcORLVqlXDqVOnAAAXLlxA//79odVqcffuXXTp0gU6nQ4AsG3bNowZMwZ169bFgAEDUKNGDdy7d6/Acz148ADJyckGj6WmphrleaxZswatW7dGo0aNjDJeRaOW2AK5vBsgGRdrHaLSOXnyJGrWrInhw4ejWbNmGD16NARBgLW1Nfr27Ytz584BAMaPH49Dh/63wnKbNm0waNAg9O7dG7Vq1cKCBQsKPI9SqcTdu3cNHktPT4dGoyn1c0hNTcWCBQvyveEaWT6uWVWGZHZPm1XqDBWgcBM5DVkKnU6H3NxcWFtb8xNHogouJiYGMTExBe5TuXJleHh4IDExEU+ePEHTpk0BAPfu3cv3tvC+vr7w9fUtcNwtW7YgOjoav/zyCwDg8uXL+PDDD3Hp0iVMmjQJZ86cwezZs3H06FGsXr0aAQEB+mPDw8Nx9epVKBQKLFq0CHPnzsXu3btfOEdsbCxeffVVpKWlISEhAYMGDcLgwYMREhKClJSUUs96CgsLw86dO3HixAn9p7dUPGqpHSQ5ef8eEZUUax0yd2K+P2s0GowbNw4rVqyAr68vtFot3nzzTZw9exYvvfQStm7dijfeeAMTJkxAfHw8tm/fbnD8K6+8gt27dyM6OhqNGzfGiBEjUL169RfO89VXX2HZsmWwt7eHQqHA9OnT4eHhgUWLFuHUqVOwsipdq2H69On49NNP4ezsXKpxyHyxWVWGrG2fLgCXnaEUOQlZkri4OKxduxYTJ04s9M2KiCxbcHBwoZ96rlu3DhMmTMC+ffvw5ptvQhAEAMCYMWP00/3/65NPPsGnn35a4Lh///037t+/j+nTp+sf8/DwMMjm7++P1157DX369DE4tlevXlAoFACAYcOGYcOGDXmeIy4uDuvWrYOjoyOkUil27NiBr776Cq1atcKsWbMKzFcUEydOxPDhw3H58mXcu3cP6enpCAkJQbNmvLddUamldpBq2Kwi42KtQ+ZOzPfnsLAwxMXFGXyg4+zsrG+AtW7dGgMHDsSCBQsQFhb2QkN4yJAhAAA/Pz+0b98e165dy7NZ5efnh9u3byMhIQGPHz/G2rVrIZPJsGTJEsjl8gIzFubIkSMIDw+HQqHA+fPnkZ2djdu3b6N27drw9PQs1dhkPtisKkNS92r4IncYBkjZDSbjcXFxwaBBg+Di4iJ2FCIS2aRJk9C3b98C96lcuTIAoH///vpPbQFg8+bNBX5yWxhnZ2cMGDAg3+n5GzZsQN26dXHq1CnExcXB29tbvy0hIcHgz/l9aurt7Y0+ffogPT0dcXFxGDRoED744ANcu3YNy5cvx2effVZozoLIZDJs3LgRwNPLDaKiojB37lzOsiqGk84DkG2tQBuxg5BFYa1D5k7s92eJRIITJ07AwcHhhe0JCQnYvXs3mjZtih07dmDmzJkvbA8MDNT/Ob/36OjoaDRs2BByuVw/s8rFxQWfffYZDhw4ADs7u0Kz5ic+Ph7Z2dn6D8Ti4uLw7bffolKlSnjllVdKPC6ZFzarypCNayWs1fZBVwmbVWQ8dnZ2+lurE1HFVpTLAZ7x8PAwmPlUu3btUp17zJgxaNeuHfz8/NCqVStYWVnBzc0NtWrVwqVLl/Dtt9/i4sWL+PnnnzF8+HAcPXpU/+ntH3/8ge+++w41a9bE/Pnz9Z/i/pdSqcRPP/2EunXrIiUlBatXr8bXX3+NVq1aYc6cOXkek5ubi5CQENy/fx9paWk4f/48KlWqhICAACiVSoSFhSEoKAgADJpSBw8exKJFi/SPFTQO/c8N955QZvJGMmRcrHXI3In5/uzt7Y2ePXvitddew/Tp0/VrvzVp0gQ2NjYYMWIEpk+fjhEjRqBFixZo3749WrVqpT9+6tSpmDZtGi5cuIDw8HC0b98+z/PUrl0bkZGRkEqlOHr0KIKDgyGTybBgwYJ8G1UPHz5EQkIC1Go1bt26BYlEoj/3jRs3EBAQAFdXV4wcORIjR47UH1enTh0EBwfrsxQ0DlkONqvKkINM+/RugKkBANzFjkMWIjMzE/fv30etWrVgb28vdhwiqqDq16+P06dPY8WKFdi3bx9yc3Px0ksv4YsvvsDy5cuxadMmeHp6YurUqbh16xb27duH1157DQAwbtw4REdH48CBA+jfv3++l/Q9X7C7urpi3rx5mDdvXoG50tLS9J/E2tnZYfr06Rg1ahTefvtt3LlzB8uWLdOvs/U8V1dXNGjQoEjj0P/U0oYiKT0BAP8ngYyHtQ5R6ezYsQPLly/HypUr9Xfn3b17N0JCQhAYGKi/A+6PP/6I5cuXY+vWrfo1pt59912sWrUKtra2OH78eL6Np379+un/3KNHD/To0aPQXFu2bMHhw4dRqVIlrFu3DlZWVvrF3pcsWYJJkyahY8eOLxwXFBRkMMOroHHIckiEZxfHViAqlQoKhQJKpbJMF2xTpSbBeUU1hLRcgWa9x5bZeahiiYmJ4ToORBVMdnY2wsLCULVqVdja2oodp1SWLl2K2NhYLF26VOwoosvv51raOqW86hwACPluJBxT76D2x5fL9DxUsbDWIXNiSe/REokEFbA9QGWktHUOZ1aVIXuHpwusa9VpIichS+Lj44OPP/4YEolE7ChERMXm5+dX6oVXyXQINg6w0WWJHYMsDGsdInHwUjoyJWxWlSEraxtkC9bQZfMuOWQ8EomExRsRma033nhD7AhkRIK1A2yFbLFjkIVhrUMkjvzuQkgkBmnhu1BpZEnsIOSkix2DLEhycjJ27NiB5ORksaMQEVEFJ5GzWUXGx1qHiIjYrCpj96TVkY4XbxlKREREZO5yHSshTPCBVsc1ToiIiMh4eBlgGZvvtADt3DzQXewgZDHc3NwwbNgwsWMQkQi46KllsYSfZ3LVvpjyjz/+zdHAydZa7DhkIVjrkDmyhH/TiYyptH8n2KwqY/Y2VsjM1ogdgyyIIAgQBIHrORBVINbW1pBIJEhISICnpyf/7lsAQRCQkJAAiUQCa2vzbfLY28gAAJlqNqvIeFjrkDnhezTRi4xR57BZVcY+SfsU2mxHAHvFjkIWIjY2lrdzJqpgZDIZ/P39ERkZifDwcLHjkJFIJBL4+/tDJpOJHaXEfJLO4758DGKS/gEUtcWOQxaCtQ6ZE75HE+WttHUOm1VlTCK1gpUmU+wYZEEUCgX69esHhUIhdhQiKkeOjo6oWbMmcnNzxY5CRmJtbW3WjSoAkNvawUaihTqLN5Mh42GtQ+aG79FELyptncNmVRnTWNnDLjte7BhkQezt7REUFCR2DCISgUwmM/vmBlkWGzsnAEBOhkrkJGRJWOuQOeJ7NJFx8W6AZUxrZQ8bXZbYMciCZGVl4datW8jK4u8VERGJy87h6cwXdVaayEnIkrDWISIiNqvKmGDtAFtthtgxyIKkpqZi9+7dSE1NFTsKERFVcHZOT5tVmkylyEnIkrDWISIiNqvK2L+Bo/GObK7YMciCeHt7Y/bs2fD29hY7ChERVXB2Ci/0zPkK4c7NxY5CFoS1DhERcc2qMiZ19sH9HH7aSMYjlUohl8vFjkFERASJzAoxNlWQqrEROwpZENY6RETEmVVlrErmTXwmfAetRiN2FLIQKSkp+PXXX5GSkiJ2FCIiIsyVbYFX9AmxY5AFYa1DRERsVpUxV20iBsrOIT0tVewoZCF0Oh0yMjKg0+nEjkJERISXhAtwT/1X7BhkQVjrEBERLwMsY9b2TxcezUpLgcLVQ+Q0ZAnc3d0xatQosWMQEREBALKlDpDmpIsdgywIax0iIuLMqjJm/f+3dM5M5zRmIiIisjxqmQNkuWlixyAiIiILwmZVGbNzcAEAqNO5yDoZR0xMDBYtWoSYmBixoxARESHXygFWmgyxY5AFYa1DRERsVpUxW7dK+E7TH6kyN7GjkIVwdnZGjx494OzsLHYUIiIi3HLrgT+t2ogdgywIax0iImKzqow5uHhgmeZ1JFr5ih2FLISDgwNatmwJBwcHsaMQERHhoe8rOIiOYscgC8Jah4iI2KwqY/bWMrSR3oKQ8ljsKGQhsrOzcf/+fWRnZ4sdhYiICH5CHGpmXRM7BlkQ1jpERMRmVRmTSiVYa/0NvCKPiB2FLERKSgp27NiBlBQu2k9EROJrmHQUizTfiB2DLAhrHSIishI7QEWQKbGHoOZdcsg4vLy8MGPGDNjb24sdhYiICFI7JzgiC7laHaxl/ByUSo+1DhERsaIoB1lSe0jYrCIjkclkcHJygkwmEzsKERERrOycYSfJQUZmlthRyEKw1iEionJvVmVkZCAxMRGJiYmFXoeekpKi3/fZV25u7gvj6XS6soxcamqpA6Q5bFaRcaSmpmL//v1ITU0VOwoRERGs7FwAABlpqaLmIMvBWoeIiMq9WfXZZ5+hTp068Pf3x5o1awrct2vXrqhTp47+y9PTE9evXwcAxMfHo3379nBzc4OHhwd27dpVHvFLJNHGD+mwFTsGWQiNRoOEhARoNBqxoxAREcHKyQPhOm9kZmaIHYUsBGsdIiIq92bVl19+icTERAwaNKjQfa9cuaKfUbVr1y7UqVMHzZs3BwDMnz8fVapUQXp6On777TdMnjwZaWmmOXtpW6V52OQ0WewYZCE8PDwwfvx4eHh4iB2FiIgIkirt0ClnOVJkfF8i42CtQ0REZrNm1fr16zF+/Hj993v37sXMmTNhbW2NDh06oGHDhjhyxDTvuOcot0JaNj8ZIiIiIsvjbPv0fj1p2bmF7ElERERUNGbRrEpJScHBgwcxcuRIAEBOTg4SEhJQtWpV/T5Vq1ZFZGRknser1WqoVCqDr/L0SvJPWJz0brmekyxXbGwslixZgtjYWLGjEBGRCRC7znHWpSJEPgm2j0+W63nJcrHWISIis2hWbdmyBd26dYO3tzeAp3cIkUgk0Gq1+n00Gg2srKzyPH7x4sVQKBT6r4CAgHLJ/YyNlQyeuoRyPSdZLkdHR7Rv3x6Ojo5iRyEiIhMgdp1j66CAuyQN2vSkcj0vWS7WOkREZBbNqg0bNhhcAiiTyeDv74+7d+/qH7t79y6qVKmS5/Fz5syBUqnUf0VERJR1ZANSOwWcBC46Ssbh6OiIdu3asYAjIiIA4tc5sLJFDqygy0wp3/OSxWKtQ0REeU9FKkPZ2dlIT0+HWq1GRkYGEhMT4eLiAisrK6Snp0MqlcLe3l6//8WLF5GYmIhevXoZjDN8+HDMnz8f33//Pc6ePYvIyEh07949z3PK5XLI5fIyfV4Fkdq7Qi7JRXZWBmztHETLQZZBrVYjJiYGvr6+ov5eExGRaRC7zoFEgnSJI4RspXgZyKKw1iEionKfWbVt2zbUqVMHp06dwvLly1GnTh3cunULAPDxxx9j2bJlBvvv3r0bkyZNgkwmM3h83rx5CAwMRNeuXbFu3Tr8+uuvJvtmZu3oBgBIT0kUOQlZguTkZPz4449ITk4WOwoREREAIFPqBCmbVWQkrHWIiEgiCIIgdojyplKpoFAooFQq4ezsXObnu3ovDPM3H8Q377yBmn7uZX4+smwajQYqlQrOzs75rtNGRETmq7R1SnnXOQDwydqdyJQp8PX4l8vlfGTZWOsQEVmuotYp/Ne/HDi6eOBfoRpUuRKxo5AFsLKygpubm9gxiIiI9JTOtRGtzBY7BlkI1jpERGQWC6ybO4UsB59Y/Qht9HWxo5AFUCqV+OOPP6BU8nILIiIyDe2yTqN7ys9ixyALwVqHiIjYrCoHzvY2GGt1BNKEu4XvTFSInJwchIeHIycnR+woREREAICa6pvopD4pdgyyEKx1iIiIlwGWA1s7R+QIVtBmpoodhSyAp6cn3nrrLbFjEBER6Qm2LnDQZYgdgywEax0iIuLMqvIgkSBN4ghdVqrYSYiIiIiMTmrnAidkIFerEzsKERERWQA2q8pJutQRkuxUsWOQBYiLi8M333yDuLg4saMQEREBAGQOrnCSZEGVkSV2FLIArHWIiIjNqnJy2r4n7trUFzsGWQB7e3s0bdoU9vb2YkchIiICAAjeDbFW8wqbVWQUrHWIiIhrVpWT0+5DYS2TYozYQcjsOTk5oVOnTmLHICIi0rMOaIIvNMPRIlcmdhSyAKx1iIiIM6vKSRVZItxUt8WOQRYgJycHUVFRvEMOERGZDIW1Fs0k95ChTBI7ClkA1jpERMRmVTnppfoFE5OXiR2DLEBSUhLWr1+PpCT+DwEREZkGV0GJX+ULIIkOETsKWQDWOkRExMsAy4nOzg2OOpXYMcgCeHh4YPLkyXBzcxM7ChEREQDA1skDAJCTxuYClR5rHSIiYrOqnEgdPKAQ0iDodJBIOaGNSs7a2hre3t5ixyAiIvofGweoYQNteqLYScgCsNYhIiJ2TcqJzMkDckkuMjM4u4pKR6VS4fjx41Cp+LtEREQmQiKBSuoMSWay2EnIArDWISIiNqvKiVzhixjBDcoUTo+n0snOzsbt27eRnZ0tdhQiIiK9ROtKyMzVih2DLABrHSIi4mWA5URWrSPaqFfhN4k7/MQOQ2bNy8sL06ZNEzsGERGRgXXVViIyJQt9xA5CZo+1DhERcWZVOXFzsAEAJGfyFrxERERkeVwdbJCUoRY7BhEREVkANqvKiYsc+Ev+Duzv7RU7Cpm5+Ph4rFq1CvHx8WJHISIi0uue+BO+TpstdgyyAKx1iIiIzapyYmtrCxdkQEiLFTsKmTlbW1vUqlULtra2YkchIiLSs7eWwF+IgVYniB2FzBxrHSIi4ppV5UgpVUDI4ALrVDrOzs7o0aOH2DGIiIgMWDl5whVpUGbmwM1RLnYcMmOsdYiIiDOrylG6TAFZFptVVDq5ubmIj49Hbm6u2FGIiIj0bJw9YS3RIjUlQewoZOZY6xAREZtV5SjLSgFrdYrYMcjMJSYm4ocffkBiYqLYUYiIiPTsFF4AgLRkrjNEpcNah4iIeBlgOTrq9xYiVVo0ETsImTV3d3eMGzcO7u7uYkchIiLSsw9shkHq+ZgIV7GjkJljrUNERGxWlSO1W138m8BPG6l0bGxsEBAQIHYMIiIiAwoXd1xBHSRmy8SOQmaOtQ4REfEywHJUT3MLo9I2iB2DzFxaWhrOnj2LtLQ0saMQERHpSSXAfNufIY++IHYUMnOsdYiIiM2qcuSvjcQo4QC0Go3YUciMZWZm4tKlS8jMzBQ7ChER0f9IJBgoHIdz4lWxk5CZY61DRES8DLAcWTt5QioRkJwcDzcvP7HjkJny9vbGzJkzxY5BRET0gnSZAtKsZLFjkJljrUNERJxZVY5sFZ4AgLTkWJGTEBERERlflpULrNVsVhEREVHpsFlVjhxcvQEAGalcZJ1KLiEhAcHBwUhISBA7ChERkQG13AXynBSxY5CZY61DRERsVpUjJw9//KTpjiSdo9hRyIzZ2NjA398fNjY2YkchIiIyEOHVFae1jcWOQWaOtQ4REbFZVY4ULu5YqBuHcClvxUslp1Ao8Morr0ChUIgdhYiIyEBCzcEIzuoMrU4QOwqZMdY6RETEZlU5kkolaO4Qj9y4B2JHITOm0WiQmpoKDe8qSUREJsZPrkZT3EVSerbYUciMsdYhIiI2q8rZZ8IqNAjfJHYMMmMJCQn49ttvuY4DERGZnKrKi9gtX4jERL5HUcmx1iEiIjarylmGjQfk2XzjpZJzc3PDyJEj4ebmJnYUIiIiA44elQAAaQmRIichc8Zah4iIrMQOUNGobT3hrrwldgwyY3K5HNWqVRM7BhER0QsUXk/X5cxMjhI5CZkz1jpERMSZVeVM5+ANhTZZ7BhkxtLT0/HPP/8gPT1d7ChEREQGrBW+AICc1BiRk5A5Y61DRETlPrMqLS0NGRkZAABnZ2fY29sXeoxGo0FycjJ0Oh18fHwAACqVCpmZmfp9bG1t4eLiUiaZjUrhj6wnVtBqdZDJ2Cuk4ktPT8fp06dRtWpVODo6ih2HiIjof2wcECfxQkZmhthJyIyx1iEionLvlixevBhBQUGoWrUq1q5dW+C+giBgzpw5UCgUqF+/PoKCgvTbPvzwQ9SuXRtBQUEICgrCBx98UMbJjSOt7lB0VH+L5MxcsaOQmfLx8cGcOXP0jVsiIiJT8kHANhy27i52DDJjrHWIiKjcm1VffPEFYmNjMXDgwEL33bBhA/bt24ebN28iISEBsbGxBtuXLFmC2NhYxMbGYt26dWUV2ai8nG0BAPFpvKUzERERWR4vJzni09RixyAiIiIzZtLXoa1evRqffPIJvL29kZv74kwknU4HlUolQrKS85Gp8Ld8KnLvHxc7CpmpxMREbNy4EYmJiWJHISIiesGQpB8wM+kTsWOQGWOtQ0REJt2sun//Pv766y/4+fnB2dkZb775JnQ6HQBAoVDg888/R6VKlRAQEIA9e/bkO45arYZKpTL4Eoubuxf8JMnITeYtnalkrKys4ObmBisr3syTiIhMq84BAHsbKfy0URAEQdQcZL5Y6xARkUk3qyQSCbKzs5GUlISIiAj8+eef+OWXXwA8vQQwKioKaWlp+OGHHzBq1CgkJ+d9l73FixdDoVDovwICAsrzaRiwkdsiBc7QqHiXHCoZFxcX9O/f3zxuKEBERGXOlOocAJA5+8ITqUjl+pxUQqx1iIjIpJtVVatWRZ8+fSCTyeDh4YEOHTrgwYMHL+z36quvws3NDWFhYXmOM2fOHCiVSv1XREREWUcvUIrUDdL0OFEzkPnSarXIyMiAVqsVOwoREZkAU6tzbFx94SzJREJKqqg5yHyx1iEionJvVmVmZiI2NhbZ2dlIS0tDbGwsNBoNAEClUiE9PV2/75gxY/Ddd9/hzp07OHPmDPbv348OHToAAOLj4xEbG4uHDx9i0aJFyMnJQd26dfM8p1wuh7Ozs8GXmDKs3WCdlSBqBjJf8fHxWLp0KeLj48WOQkREJsDU6hxHt0oAAGW8uE0zMl+sdYiIqNwvBP/5558xZ84cAMC5c+ewevVqHDlyBI0bN8aCBQvg5uaGuXPnAgDeffddxMTE4NVXX4VCocDixYvx0ksvAQA6dOgApVIJBwcHNGrUCEePHoW9vX15P50SOVxpKh4rdWgqdhAyS66urhg6dChcXV3FjkJERPQC5xqt8Zr6UwzPdUYLscOQWWKtQ0REEqECrn6pUqmgUCigVCpF+fRxyeG7OHA9GudmdSn3cxMREZFpK22dInadAwDNFx3DyNZV8G63mqKcn4iIiExTUesUk16zylLVxyO8nb4aOl6HTyWQkZGBy5cvIyMjQ+woREREeZpp/SucI06KHYPMFGsdIiJis0oE/lZKvCE7juS4SLGjkBlSqVQ4dOiQ6LcmJyIiyk9nzTn4JJ0XOwaZKdY6RERU7mtWEeDkFQgASIp5BA+/QJHTkLnx9fXF/PnzxY5BRESUrwxbH9hnxYodg8wUax0iIuLMKhG4+1UDAGTEh4sbhIiIiKgMaBz94Jobjwq4NCoREREZQbFmVnl4eBRpv8TExBKFqSgUbl7IEmyQk8xbOlPxJSUl4Y8//sDLL78Md3d3seMQERG9QKKoBO/oM1BlaaCwtxY7DpkZ1jpERFSsZlVSUhJSUlIK3Ie3mC2cRCrFLvlrkEoC0VrsMGR2pFIp5HI5pFJOjCQiIhNVtSN+vhmNbimZUNgrxE5DZoa1DhERFatZtWXLFri4uBS6DxXumNc4OOmsMFLsIGR2XF1dMXjwYLFjEBER5UtRtzO+2atDPWU26lVis4qKh7UOEREV6+OKli1bQqlUYsuWLYiKispznxEjRhglmKWra6eEIjFE7BhkhnQ6HXJycqDT6cSOQkRElCcPOynayu4gJZ53PqbiY61DRETFalZ169YNbdu2xdmzZ9G/f/8yilQxdMk+ipmpX4gdg8xQXFwcFi9ejLi4OLGjEBER5Ukm5GK79WeQR/wpdhQyQ6x1iIioWJcBRkZGQhAEXL16FQ4ODmWVqUKQufjD/UkqctTZsJHbih2HzIiLiwtee+21Qi/JJSIiEo2NA9IljpCoOLOKio+1DhERFWtmVZcuXZCVlQUbGxt4e3uXVaYKwc6jMqQSAYnR4WJHITNjZ2eHhg0bws7OTuwoRERE+VLaeMMmPVrsGGSGWOsQEVGxmlXHjx+Hre3TWUCRkfykrDScvasAAFJjw8QNQmYnKysLN27cQFZWlthRiIiI8pVt5wMHNS/jouJjrUNERMW6DPB5Fy9exMOHD6HVavWPcXH1ovOoVA1RgjtSVCqxo5CZSU1Nxd69ezFx4kR+4khERCYrx60WUpPuQK3RQm4lEzsOmRHWOkREVKJm1ZtvvokTJ06gVatWkMn+V3ywWVV0Dk4u6GgVjNFWVdBO7DBkVnx8fDB37lyDv3tERESmJqXdPEy9fQEnU7JQzdNR7DhkRljrEBFRiZpV27dvx8OHD+Hj42PsPBWKv5s9IpIyxI5BZkYikcDKqsSTIomIiMpFZTd7AAKeJGWwWUXFwlqHiIiKtWbVM7Vq1UJ6erqxs1Q4M3PXYnToNLFjkJlJSUnBzz//jJSUFLGjEBER5cs35wluyscjM+yi2FHIzLDWISKiEn1ksWHDBgwaNAj9+/fXL7gOALNnzzZasIrA2sEFHqlRYscgMyMIArRaLQRBEDsKERFRvmTOPnCUZCM7gTeToeJhrUNERCVqVn3yySeQy+XIyMhAdna2sTNVGFZuVeAVmYSc7CzY2HLxSCoaNzc3vPHGG2LHICIiKpidCzIkjpCkhoudhMwMax0iIipRs+rkyZOIioqCi4uLkeNULHbe1SGVCIiLDEVAjYZixyEiIiIyqlRbP9imR4odg4iIiMxMidasatmyJR4/fmzsLBWOW6VaAICUyPsiJyFzEhMTgwULFiAmJkbsKERERAVSOwbARR3Ny7moWFjrEBFRiWZWubi4oEePHujTp4/BmlWrVq0yWrCKwMu/Ol7O/QojrBuikdhhyGwoFAr06dMHCoVC7ChEREQFetzsI3y47x6OZObCzcFG7DhkJljrEBFRiZpVXbt2RdeuXY2dpcKxsrZGhqIWnii1YkchM2Jvb4+mTZuKHYOIiKhQngE1kYA4PEnOZLOKioy1DhERlahZNXXqVGPnqLDG2xyFS2g2gOViRyEzkZWVhfDwcFSpUgV2dlyYn4iITFegNBarrFciLtIHCGgmdhwyE6x1iIioRGtW5bWwOhdbL5k60kjUV54ROwaZkdTUVOzatQupqaliRyEiIiqQk40VXpWdR1rUXbGjkBlhrUNERCWaWaVUKiEIAiQSCQBAEAQolUqjBqsoBJcq8Ek6BkGng0Raot4hVTDe3t748MMPIZfLxY5CRERUMEUAtJBCmxAqdhIyI6x1iIioRN2RSpUq4cGDB/rv7927h0qVKhktVEUi96kNR0kWkmIjxI5CZkIqlcLOzg5SNjeJiMjUWdkgxcYPcuUjsZOQGWGtQ0REJXoHGDVqFIYPH459+/Zh7969GD58OEaPHm3sbBWCW+V6AIC48FsiJyFzkZKSgj179iAlJUXsKERERIXKcKoC16wnEARB7ChkJljrEBFRiS4DXLhwITw9PbFq1SpIJBKMGjWKi66XkG+VuvhaMwS11C6oL3YYMgs6nQ4qlQo6nU7sKERERIVKqjMcO0/dRd10NbycbMWOQ2aAtQ4REUmEYnzM9eyuHOZOpVJBoVBAqVTC2dlZ7Dh46etT6F7XG/NerSd2FCIiIhJZaesUU6tzHsSlofvys/h5Ymu0quYudhwiIiISUVHrlGJdBjhs2DA0bdoUn376Ka5evVrqkPRUF6dIKJ4cFTsGERERkdFVdtRgoOwsoqMeix2FiIiIzESxmlX//PMPfv/9d/j6+uKjjz5CrVq1MG3aNJw8eRIajaasMlq8nro/0T8hWOwYZCZiY2PxxRdfIDY2VuwoREREhZJrs7HMeg00Ty6JHYXMBGsdIiIq9gLrvr6+mDRpEv744w+EhISgffv22LBhA2rVqoWRI0eWRUaLJ/GoCV9dLHLU2WJHITPg5OSErl27wsnJSewoREREhXPyQZbEDpLEB4XvSwTWOkREVMIF1p9xcnLC66+/jtdffx25ubk4ffq0kWJVLA6V6sDqtg5PHt9F5VpBYschE+fg4IBWrVqJHYOIiKhoJBKk2AXCLi1M7CRkJljrEBFRiZpV06dPf+ExhUKB1q1blzZPheRdtSEAIPnxbTarqFBqtRoREREICAiAXC4XOw4REVGhchTV4BX1ENm5Wthay8SOQyaOtQ4RERX7MkAA0Gg0+O233+Dm5gY3Nzfs27cPERER+OCDD/DJJ58YO6PF8/CpjL+FhohO47pfVLjk5GRs27YNycnJYkchIiIqEkmVdnio88WjhAyxo5AZYK1DREQSQRCE4h7UokULbN26FbVr1wYA3LlzB6NGjcKWLVvQrVs3REZG5ntsamoq0tPTAQAuLi5wdHQs9Hw5OTlITEyETqeDv7+/wbakpCQ4OTnBxsamyPlN7ZbOAPDqd3+ivq8CSwY1EjsKmTitVouMjAw4ODhAJuOn00RElqa0dYop1jmq7Fw0+vQovh0ahH5BlcSOQyaOtQ4RkeUqap1SoplV9+/fh7u7u/57Dw8P3L9/H7Vq1UJ0dHSBxy5btgytW7dG7dq1sX79+gL31el0+OCDD+Di4oJmzZoZXGYYHR2N5s2bIzAwEB4eHtiyZUtJnorJqOFui9Q43tKZCieTyeDs7MzijYiIzIazrTXqOOfiUVSc2FHIDLDWISKiEjWrXn31VYwcORLHjh3DsWPHMGLECPTt2xdSqRSFTdT67LPPEBkZiQEDBhR6nvXr1+PQoUO4d+8eYmJiDGZsffLJJ2jYsCFUKhWOHz+OadOmQalUluTpmIQh2bvwZcLbYscgM6BUKnHgwAGz/n0nIqIKJicDh3NGQxF2SOwkZAZY6xARUYmaVWvXrkWbNm3wxRdfYPHixWjfvj2Cg4MBADExMUYLt3r1anzyySdwcXGBWq022Pbbb79h+vTpkEqlaNmyJYKCgnD06FGjnbu8yf3qwxUqJMZGiB2FTFxubi5iY2ORm5srdhQiIqKisXFAitwPdqn3xU5CZoC1DhERlehugA4ODpg/fz7mz5//wjYfH59Sh3omNDQUp0+fxuTJk5GdnY0hQ4Zgw4YN0Gg0SEhIQGBgoH7fwMDAfNfKUqvVBs0ulUpltIzG4lm9CXAeiLl/BR4+AWLHIRPm4eGBN998U+wYRERkIsyhzgGALJfaqBT9CJk5GtjblKgEpQqCtQ4RERVrZlXz5s2Nsk9RSaVSaLVaJCYmIjIyEufPn8fPP/8MmUwGqVQKjeZ/d8/Lzc3Nd5H1xYsXQ6FQ6L8CAkyvGeRXtT6yBWtkRFwXOwoRERGZEXOocwDAyqc+aksj8CAuXewoREREZOKK9bFWSEgI5s2bV+g+xlKtWjX07t0bUqkUbm5uaNeuHR49egSZTIaAgADcuXMHHTp0APD0joTDhw/Pc5w5c+ZgxowZ+u9VKpXJFXIyKyuEWVVGbjIXWaeCxcXF4aeffsKoUaPg7e0tdhwiIhKZOdQ5AOBSpTF014Lxd2Q0Gge4iB2HTBhrHSIiKlazavHixaXeJyMjAykpKcjMzIRSqURkZCS8vb1hbW2NlJQU/d0/AGDMmDFYuXIlqlWrhvj4eOzfvx+7d+8GAIwcORJz587Fd999h7NnzyIuLg7du3fP85xyuRxyubw4T1UU62quwd3EHHQQOwiZNAcHB7Ru3RoODg5iRyEiIhNgLnWOTcP+6HzMDd2SgcJvs0MVGWsdIiIqVrNq9uzZpT7hL7/8op+ddfHiRaxbtw6HDh1Co0aNsHjxYri6umLOnDkAgGnTpiEuLg6DBg2CQqHA0qVL0bFjRwDA3LlzkZSUhIEDB8LPzw979+6FtbV1qfOJqYafB/bfug+dToBUKhE7DpkoR0dH/YxCIiIis2Flgxo+LrjHywCpEKx1iIhIIgiCIHaI8qZSqaBQKKBUKvWzuExByIUzcPv9TViP2gP/Gg3EjkMmKicnB7GxsfDx8cl3nTYiIjJfpa1TTLXOAYBbq4fhWpIMw+dvFTsKmTDWOkRElquodUqxFlinshUYEIiq0jjEP7wqdhQyYUlJSdi0aROSkpLEjkJERFQsrrZS1NHcRbwqW+woZMJY6xAREZtVJsTdpzJS4YicqH/FjkImzNPTE1OmTIGnp6fYUYiIiIrFPqARakkicStKKXYUMmGsdYiIiM0qEyKRShFtUxU2yXfFjkImzMrKCh4eHrCyKtaSc0RERKJTBDaGkyQL4Y/uiR2FTBhrHSIiKlGzaubMmQgLCzN2FgKQpqgN98yHYscgE6ZSqXDkyBGoVCqxoxARERWLxLs+ACD9yXWRk5ApY61DREQlala5urqia9eu6NevH06cOGHsTBVaYr3ReCt7KjLUGrGjkIlSq9V4+PAh1Gq12FGIiIiKR+GPnbW/xf7UqmInIRPGWoeIiErUrJo3bx5CQ0MxadIkfPvtt2jcuDGCg4ORmZlp7HwVTtXaQbitq4zbMfwkifLm6emJt99+m+s4EBGR+ZFI4FCvBx4oJUjOyBE7DZko1jpERFTiNaukUim6du2KoUOHQhAErF27Fo0bN8ZPP/1kzHwVTk1vR3xo8wtUIb+KHYWIiIjI6JrjDr6wWo+bXGSdiIiI8lGiZtW9e/cwc+ZMVK9eHUePHsWGDRsQEhKCQ4cOYcqUKcbOWKFYy6ToZnMLisdHxI5CJio+Ph4rVqxAfHy82FGIiIiKzVuWhjesTuJhGNfopLyx1iEiohLdYqNXr16YPHkyrl+/Dnd3d/3jNWvWxKRJk4wWrqJKcWkAn+RLYscgE2VnZ4dGjRrBzs5O7ChERETFJq3UBACQ9TgEQGtxw5BJYq1DREQlmlnVsWNHzJo1y6BRNWbMGADA0qVLjRKsIpP6BSFAG4l0VYrYUcgEOTk5oUuXLnBychI7ChERUfG5VEamzBny+BtiJyETxVqHiIhK1Kz677pUGo0G27dvN0ogAtxrtYJUIuDJ7QtiRyETlJubi5iYGOTm5oodhYiIqPgkEmS4N0Bl9X3Ep2WLnYZMEGsdIiIq1mWAkydPzvPPT548QZMmTYyXqoKrXKsJVmoHw1tlj3pihyGTk5iYiLVr12LixInw9fUVOw4REVGxyVpPxN7dV6B9nIpeDXzEjkMmhrUOEREVq1kVFBT0wp8lEgk6deqEnj17GjNXhWZlI8dJn7GokmyPIWKHIZPj4eGBiRMnwsPDQ+woREREJeLWdACuHXGG/5MUNqvoBax1iIioRDOrevXqhSpVqpRFHvp/nT2UyHl0CABnrJEha2trfspIRETmTafFRNfLuPYwAEBdsdOQiWGtQ0RExWpWTZgwAevXr8eiRYvy3L5+/XqjhCKgrTwMLdSroEyaBoW7t9hxyISkpaXh0qVLaNGiBRceJSIi8ySRYkjyGmRldYRaMwRyK5nYiciEsNYhIqJiNatat25t8F8qO36NOgHXgPDrp9C4y1Cx45AJycrKwo0bN9CgQQMWcEREZJ4kEuT4tkCTsHu4GaVCs0BXsRORCWGtQ0REEkEQBLFDlDeVSgWFQgGlUglnZ2ex4+RJ0OmQtLAqHvj1Q5uJK8WOQ0REROWktHWKOdQ5AKA9txI5xz7Dts7nMKFTbbHjEBERUTkoap0iLcngS5cu1f/5s88+Q9++fXHjxo2SDEX5kEileOLQEIrEK2JHISIiIjI6WWAb2ElykPjgkthRiIiIyMSUqFn1wQcfAABCQkKwZcsWdOrUCZMmTTJqMAKUVXrhQnYAcjQ6saOQCUlISMD333+PhIQEsaMQERGVnG9j3Hd9Cbdj01ABJ/pTAVjrEBFRiZpVMpkMWq0Whw8fxoABA/DOO+/g6tWrxs5W4bm0Go4FOSNwK1opdhQyIXK5HNWrV4dcLhc7ChERUclZ2SD25Q04m1EZD+LTxU5DJoS1DhERlahZVadOHWzbtg179uxBjx49YG1tDbVabexsFV59PwWqWyXi4e0QsaOQCXF2dkbPnj1Neh0SIiKiomhRyRbtrO7i7wecQUP/w1qHiIhK1KwKDg7G9u3b0a1bN3Tt2hUAUL9+faMGI8DGSorv7Neh2r8rxI5CJkSj0SAxMREajUbsKERERKViF/kXtlktxL17N8WOQiaEtQ4REZWoWdWuXTscPnwYS5Ys0T928yaLjLKg9GiKwIwbEHRct4qeSkhIwOrVq7mOAxERmb+AlgAAScR5aHVct4qeYq1DRERWJTkoOTkZ69evx8OHD6HVavWPr1+/3mjB6CnH2p3gHv0jwu9dQZW6zcWOQybA3d0dY8eOhbu7u9hRiIiISsfeDZmuddA08V/cjFKicYCL2InIBLDWISKiEjWrBg4cCDc3N/Tu3RsymczYmeg5NZp3R85JK8RePcxmFQEAbGxsULlyZbFjEBERGYW8dle0T96BvaGJbFYRANY6RERUwmbV5cuXERMTA0dHR2Pnof+wc3DCebs2iExIETsKmYj09HRcvXoVTZo04d9BIiIye7IaXZBx5QSuPwgDOtcQOw6ZANY6RERUojWr2rdvj/v37xs7C+UjpNUKfJrcA7larltFQEZGBs6fP4+MjAyxoxAREZVejW441WEHTj3RICtHW/j+ZPFY6xARUYlmVikUCvTq1Qt9+/aFra2t/vFVq1YZLRj9T7saHvjhyFXcehiOoFrVxI5DIvP29sYHH3wgdgwiIiKj6VLbA5t+P4u/QhPRrZ632HFIZKx1iIioRM2q9u3bo3379sbOQvlo6OeEs7bv4d65N4BaX4kdh4iIiMioqt3fgCO2X2PxncZsVhEREVHJmlVTp041dg4qgEwmQ7h9Yyhi/xI7CpmAxMRE7N27FwMGDICHh4fYcYiIiEovsB0c8Sli7vwDQQiCRCIROxGJiLUOERGVaM2qlJQUvPvuu+jUqZP+sdmzZxsrE+UhN7AjaqrvQJWaJHYUEpm1tTV8fHxgbW0tdhQiIiLj8GsKjbUjGmSF4Fa0Suw0JDLWOkREVKJm1Ztvvgk7OzucOXNG/9iSJUuMFopeVLl1f1hLtHjw9z6xo5DIFAoF+vTpA4VCIXYUIiIi45BZQVqzO3pYXcXxO3FipyGRsdYhIqISNauOHTuGjz/+WP+9TqeDVFqioaiIfANr476sJmLC7ogdhUSm1WqhUqmg1fKOSUREZDmktV+Gj00mzt1+InYUEhlrHSIiKlGHydHREdnZ2frvb968icqVKxstFOXt95Y/YW5iD2i0OrGjkIji4+OxfPlyxMfHix2FiIjIeBoOwl+9j+NytBoRyZlipyERsdYhIqISNavGjBmDuXPnAgB2796N4cOHY8KECUYNRi/qUs8PGVlZuPrgsdhRSERubm4YPnw43NzcxI5CRERkPFIZutbzgbdVBn7/N0bsNCQi1jpERCQRBEEo7kFarRbff/89Dhw4AEEQ0LdvX0yZMqVIlwImJydDpXq6cKabmxucnZ3z3TcpKQlpaWn67x0cHODp6VnotsKoVCooFAoolcoCz29qdDoBYQsbIsG7PVq/tUbsOERERFQGSlunmGudAwC4ug05+6djmOt2/Dqtu9hpiIiIyMiKWqeUaGaVTCbDO++8g6NHj+LYsWN45513irxm1bfffotOnTqhQYMG2LhxY4H7zp07F82aNUOnTp3QqVMnfPLJJ0XaZqmkUgmS3ILgH38ago6XAlZUGRkZuHDhAjIyMsSOQkREZFyVW8NGyIF77F8IT+T7XEXFWoeIiKyKuuP06dML3WfFihWF7rNgwQIsWLAAI0aMKNJ5P//8c0yePLnY2yyVvGE/+J85gIe3LqJ6w9ZixyERpKWl4cSJEwgMDISDg4PYcYiIiIzHvTp0HnXQJ+ESDt6IxtQuNcVORCJgrUNEREWeWVWjRg39l0ajwb59++Dq6go3Nzfs27cPGo2mTALm5ubmu7hiQdssVb12fZEKR8T/s03sKCQSHx8ffPTRR/Dx8RE7ChERkdFJGw5Cd2kIjlwLFzsKiYS1DhERFXlm1dSpU/V/btGiBf744w/UrVsXADB48GCMGjXK6OHc3d2xfPlyLFy4EBKJBCtXrsTQoUML3fZfarUaarVa//2zNbPMkbWNHPfcuiA39h4EQYBEIhE7EhEREYnIkuocAECD1yA79y2QcBc3o1qiQSWF2ImIiIionJVozar79+/D3d1d/727uzsePHhgtFDPfP7553j06BESEhKwbds2jB8/HklJSYVu+6/FixdDoVDovwICAoyetTzpXv4KozLfxbWIVLGjkAiSkpKwefPmfH/fiYioYrG0Ogfu1SH54CFiHevil8sRYqchEbDWISKiEjWr+vfvjzfeeANHjx7F0aNH8cYbb6Bfv37Gzmage/fucHNzw+PHj4u1DQDmzJkDpVKp/4qIMO/Cp2V1H3g62uDk5ZtiRyERSKVSODs7F/mmBkREZNksrc4BACsbOd5opMAfV8OQnasVOw6VM9Y6RERU5MsAn7dmzRosX74cX375JQCgW7dueO+994p0bFpaGpKSkpCRkYHk5GSEh4fDz88PNjY2SEpKglQqhaurKwDgyZMn0Ol0yMjIwI4dO6DValGnTp1Ct/2XXC6HXC4vyVM1STKpBF95Hkb9G3ug7fsAMqsS/RjJTLm6uuK1114TOwYREZkIS6tzAACZyXj32it4kjMOx243Q5/GfmInonLEWoeIiCSCIAjlecKffvoJ8+fPN3jswIEDaNiwIWbPng1XV1fMmjULANCoUSOoVCo4ODigUaNGmD9/vn6drIK2FUalUkGhUECpVMLZ2dm4T7Cc3L9yBrX298WNlzagUedBYsehcqTT6aBWqyGXy/mJIxGRBSptnWIJdQ4AYNMr+DcmDV95f40t41uJnYbKEWsdIiLLVdQ6pdybVabAEoo4QadD+KIgpNhXQdP394sdh8pRTEwM1q5di4kTJ8LX11fsOEREZGRsVv2/azuAfZPRQb0CW99/HYHuDmInonLCWoeIyHIVtU7hRxVmSiKVIq7GYDRIO4eUhBix41A5cnFxweuvvw4XFxexoxAREZWdev0gyJ0wUv4nfvon73VJyTKx1iEiIjarzFjt7hMQBzecu/CP2FGoHNnZ2aFu3bqws7MTOwoREVHZsbGHpNEQtPTSYtelCGSoNWInonLCWoeIiNisMmOunr5YXHMnVj3wQAW8mrPCyszMxJUrV5CZmSl2FCIiorLVeyk831iDjBwN9lyNEjsNlRPWOkRExGaVmRvaMhDKuHD8++81saNQOVEqlThw4ACUSqXYUYiIiMqWRIJKTtZ4q1oiNv8VBp2OH85VBKx1iIiIC6yb88KjAHRaHcI/D0IqF1onIiKyCFxg/T+u7YCw7y10VH+Dj0f0Ro/6PmInIiIiohLiAusVhFQmRUKtYWiU9ifiIkPFjkNERERkXPX6QWKrwIduf2L1qVAufUBERFQBsFllAer3noxsyPHoj+/EjkLlIDk5Gdu3b0dycrLYUYiIiMqejT3QbDRezjmGB5Fx+Cs0SexEVMZY6xAREZtVFsDR2RU3vfugTtSvyMpIEzsOlTGJRAKZTAaJRCJ2FCIiovLRYgJkmgy8434Jq049EDsNlTHWOkRExGaVhQjoNRPntA2w/8JdsaNQGXN1dcWQIUPg6uoqdhQiIqLy4VIZkjZT0bphbZx/lIyQx5xxY8lY6xAREZtVFqJStTo4UX8xvr2gQo5GJ3YcKkOCIECj0XDNDiIiqlh6fIbGPcagtrcTlh29z/dBC8Zah4iI2KyyIG91qoEWacdx8fAWsaNQGYqNjcXnn3+O2NhYsaMQERGVK2lqONb4HsQ/DxNwLjRR7DhURljrEBERm1UWpLaPE8a4XEdgyJfQajRix6Ey4uLiggEDBsDFxUXsKEREROUrPR5V7wZjsvc9LDl8FzodZ95YItY6RETEZpWFse/6AQKEaFw7+qPYUaiM2NnZoVGjRrCzsxM7ChERUfmq3AoIbI8pVr/hZpQSv/8bI3YiKgOsdYiIiM0qC1O7WWfcsG0Gj8vLObvKQmVlZeHff/9FVlaW2FGIiIjKX8f34Zh0AzMqP8TSo/eg1mjFTkRGxlqHiIjYrLJAdj0/RaAuApd+3yB2FCoDqamp2LNnD1JTU8WOQkREVP6qdQKqdsRE3S5EpmRiw7kwsRORkbHWISIiiVABb7OhUqmgUCigVCrh7Owsdpwy8d3aNfglqRqOfdAVciuZ2HHIiHQ6HTQaDaysrCCVst9MRGRpSlunVIQ6Bwn3AJk1Pvs7G9svPMHxmS+hkgsvGbMUrHWIiCxXUesU/utvoXoPGIGoNA1+PndX7ChkZFKpFDY2NizeiIio4vKsDbhVw/ROleEiF/D577fFTkRGxFqHiIj4DmChqns64ouq/+Ll072hSk0SOw4ZUUpKCn755RekpKSIHYWIiEg8OZlwWt8Wa+tex6F/Y3HmfoLYichIWOsQERGbVRasc+/X4SBk4fb2OWJHISPS6XRQq9XQ6XRiRyEiIhKPjT1QrRMaPPgevarI8NGef6HKzhU7FRkBax0iImKzyoJ5VaqK69UnoXncLwi7fUnsOGQk7u7uGDFiBNzd3cWOQkREJK6un0ACCZa6/QZlVi4WHeTlgJaAtQ4REbFZZeGaDZmLGKkPMvfNgMBPp4iIiMiSOLgDXT+G4+3tWN5eg12XI3HsdpzYqYiIiKiU2KyycDZyW6R2+hw3M91w8Cpv7WwJYmJisHDhQsTExIgdhYiISHzNxgL1+qFbTRd0q+uND3ZfR2RKptipqBRY6xAREZtVFUDDl17Dmbrz8cmhR0hMV4sdh0rJ2dkZvXv3ttzbkRMRERWHVAa8/hMkVdph6eBGcJRbYcr2q1BrtGInoxJirUNERGxWVRAL+zVAD9053AsezcsBzZyDgwOaN28OBwcHsaMQERGZjvR4uOx5Axt6O+JOtAqf/35H7ERUQqx1iIiIzaoKwsNRjsEtqqBd2h8I+X2d2HGoFLKzs3Hv3j1kZ2eLHYWIiMh0yJ2AlHDU/mcWPnmlJn765zF2XYoQOxWVAGsdIiJis6oCadZ7LC47dUWdy58gMvSm2HGohFJSUrBz506kpKSIHYWIiMh0WNsBA4KB2Bt4I/0nvNGqMubu+xcXHiWJnYyKibUOERFJBEEQxA5R3lQqFRQKBZRKZYW7Fj5NmQzlirbIltrD//0/YWvH6dXmRqvVIjs7G7a2tpDJZGLHISIiIyttnVKR6xwAwN/fAUfnQTN0J0b96Yo7MSr8NqU9Krvbi52Mioi1DhGR5SpqncKZVRWMk8IN2f03IFVjhW8PXhI7DpWATCaDg4MDizciIqK8tJkK1OsPq6xkfD+8KRR21hj/4yWosnPFTkZFxFqHiIjYrKqAajRuh7u9duGHkAz8evmx2HGomFJTU7Fv3z6kpqaKHYWIiMj0SCTA4M1Ak+FwsbfBhlFNEavKxjvbr0Kj5U1mzAFrHSIiYrOqghreOhATG8rQ8EBv3L9yWuw4VAwajQbJycnQaDRiRyEiIjJNEsnT/x7/FNX/+gDfv9EE50IT8fFvt1ABV8AwO6x1iIiIzaoKSiKRYObAjhCsHeG6fwwSo8PFjkRF5OHhgXHjxsHDw0PsKERERKbNuwFw42d0SPwZi19riB0Xn+D70w/FTkWFYK1DRERsVlVgclt7uI/bBQESJG18HdlZGWJHIiIiIjKehoOAdtOBY/Pxuss9vNetFr4+cg+/hkSKnYyIiIgKwGZVBefhF4jUvptRJfcRrge/yanxZiA2NhaLFy9GbGys2FGIiIhMX9f5QI1uwO5xmNbUCkOaB2DWrzdw7kGi2MkoH6x1iIiIzSpCraYv4XrzL/BdfCNs/Ctc7DhUCEdHR3Tq1AmOjo5iRyEiIjJ9UhkwcD3QfDwkTr5YNKAB2tf0wOStIbgdrRI7HeWBtQ4REUmECjiVRqVSQaFQQKlUwtnZWew4JmPxoTtY/2cotvVzRevW7cWOQ0REVCGVtk5hnVOIhHvIcKiMoRtCEJ2ahW1vtkIdH75ORERE5aGodUq5z6xKTExEaGgoQkNDoVQqC9w3ISFBv29oaGieU4GjoqKQmZlZVnErlFm96uBL3z/R6I/X8PDG32LHoXyo1Wo8evQIarVa7ChERETmJSsV2NAdDifm4KexLeCjsMUb6y7gTgxnWJkS1jpERFTuzarVq1ejV69eaNKkCTZt2lTgvh9//DHatm2LXr16oVevXvj888/12yIiItCwYUM0bNgQXl5e2LBhQ1lHt3hSqQSvjJuLKKsAOO0ZjvioMLEjUR6Sk5OxZcsWJCcnix2FiIjIvNi5AD0+B0I2wfXGWmyb0Aq+Clu8se48rjxJETsd/T/WOkREJNplgCNGjEDz5s0xffr0fPeZPHkygoKCMHny5Be2jRs3DtbW1ggODsb169fx0ksvISwsDK6uroWem9PjC5YQHQ7d2i5Ik7nAd/pJODi5iB2JnqPRaJCeng5HR0dYWVmJHYeIiIyMlwGWg+OfAudWAEO2IjWwByb8eBn/RimxfEgQejf0FTtdhcdah4jIcpnsZYDFpVarER0dDZ1OZ/D4gQMHMHXqVABA48aN0bRpUxw7dkyMiBbH068K0gduh68mCns3LkGORlf4QVRurKys4OLiwuKNiIiopLrMB+r1BY58BBe5BFsntEKP+j54e9sVrDr5ADpdhVvS1aSw1iEiIpNuVnl5eWHVqlVo1qwZPD09sXXrVgBATk4OEhMTUblyZf2+gYGBiIqKynMctVoNlUpl8EUFq96wNe6+ugefxbfHtB1XodGyYWUqlEolfv/990LXfCMiooqBdU4JSKXAgGBg1D5AZg1baxm+HRKEaV1rYtmx+xi96SIS07leklhY6xARkUk3qxYuXIgHDx4gJiYGu3fvxqRJk5CYmAiZTAaZTIbc3Fz9vjk5OZDL5XmOs3jxYigUCv1XQEBAeT0Fs9asRTusHt4cOXeP4NLK4dDk5ogdifD0dz0yMhI5Ofx5EBER65wSs7YD3KoB6nTg5xGQxv2LGd1r4cexLXEnRoWXv/0TR269eHMfKnusdYiIyKSbVc/r3Lkz3Nzc8OTJE8hkMlSuXBk3b97Ub7916xaqVauW57Fz5syBUqnUf0VERJRXbLPXrZ43pnbwQ4vUw/h3eX9kZ2WIHanC8/T0xKRJk+Dp6Sl2FCIiMgGsc0pJ0ALKSODHV4HH/6BjLU8cercDGlVSYNKWEEzachmxymyxU1YorHWIiKjcm1VKpRKhoaFIT09HYmIiQkND9belTUhIQFJSkn7fsLAwhIaG4vr165g9ezYEQUDdunUBAKNHj8acOXNw4cIFfP3110hKSkK3bt3yPKdcLoezs7PBFxVd015jcLPjD6iTcQkPl/eEKjWp8IOIiIioXLDOKSVbBTDqN8CnEbClP3D3ELycbLF+dHOsfqMpQh6novPS01h18gGyc7VipyUiIqoQyr1ZdeDAAfTq1Qs3b97Ezp070atXLzx48AAAsHz5cmzatEm/78CBA9GrVy+MGjUKUVFROHHiBOzs7AAAs2fPRps2bTB+/HgcO3YMBw4c4CKMZSio61A8fmU7AnIe4caqNxCWyBlWYomLi8OyZcsQFxcndhQiIiLLYKsAhu8GanYHfh4BJNyHRCLBK418cWLmSxjeqjJWHH+ArsvO4NC/MRDpZtoVBmsdIiKSCBXw3Za3dC65x3ev4JP9txCS6YWVg+uhc32ui1He0tLScPXqVTRp0gROTk5ixyEiIiMrbZ3COqcUdFrg/hGgTu+n32tyACsbAMCjhHR8/vsdnLgbj5ZV3fBJn3qo76cQMazlYq1DRGS5ilqnsFnFIq7YVNm5mLXjIqaEvY103zYIGrMMtnYOYsciIiKyCGxWmYiQH4F/VgOvrQX8gvQPn7mfgM8O3sbDhHQMbRGAmT1qw8Mx75v8EBERkaGi1ilms8A6mQ5nW2usHtUa6bUGoGnsL4j9ujVCr/8ldqwKIycnBxEREbxDDhERUVkKaAVYyYH1XYEzXwNaDQDgpVqe+OPdDvjk1Xo49G8sOn99GmvPPkSORidyYMvBWoeIiNisohKRymRoPeJTRA05DK3ECoF7+uD0pvlceLQcJCUlYePGjQY3IyAiIiIj86oDTDgBtJsOnP4C2NAdSHu6hpK1TIox7ari9PudMKBpJSw5fA89lp/B8dtxXM/KCFjrEBERLwPk9PhSy1FnI2TrPOwMs8cNly74vF8dtK3pI3Ysi5Wbm4uUlBS4urrC2tpa7DhERGRkvAzQBEVeBi6uA/p/D0hlT2dZyf53Y5/7cWn47OBt/PkgER1qeuDjV+uhljfXWiop1jpERJaLa1YVgEVc2QiNT8NHe25iZNSn8HF3Rc0RK+DiwaYVERFRcbBZZeIiQ4BfxwO9lwI1u+kfFgQBx+/E4/PfbyMiJQvDW1XGe91qwdXBRsSwREREpoVrVlG5q+HlhJ1vtoJ30MuonXoWwqoWuLz/ewg6ruFgTCqVCkePHoVKpRI7ChERUcVj7wa4BgLbBgK7xwHp8QAAiUSC7vW8ceS9jpjVqzb2XolCp6WnsebMQ2TmaEQObV5Y6xAREZtVZFRSmRQtB06HevJ5PHRqgeZX5uDGVz3wODFN7GgWIzs7G/fv30d2drbYUYiIiCoet6rAyH3AgLXAo9PAquZA+Dn9ZrmVDBM7VsfJ9zuhT2NfLDt6Dx2/OoUN58K4tmcRsdYhIiJeBsjp8WXqxqndOP33OazO6onpnatgQoeqsLaxFTsWERGRyeJlgGYkIwk48yXQac7TGVc5mYCNvcEuEcmZ+O7kA/x6JQoejjaY2qUmhjQPgI0VPzMmIqKKh2tWFYBFXPnKzNFgxfEHsPl7OQbanEdO729Qu3m3wg8kIiKqgNisMlPp8cAP7YDmY4H2MwBrww/nwhIz8O3x+/jtejT8Xe3wXrda6BdUCTKpRKTARERE5Y9rVpHJsLexwke962LA62Ogkdqg5oFBuPDdGKhSeTvikoiPj8fKlSsRHx8vdhQiIiJ6Ru4MNBsN/PkNsKa9waWBAFDVwwErhjbB4Xc7oq6PM2bsuo5eK87i8M1YVMDPjgvEWoeIiNisonJTvVFbVJt9HhfrfIAGiX9AvaIZTl26xgKtmGxtbVGvXj3Y2vJySiIiIpNhbQt0mQdM/hOwdwc2vwKcW/HCbrV9nLB2VHPsm9IO3s62/9fefYdHVWYPHP9On0ky6b0RQkgggBTpHRWQZllR2LWXxbJFseuuP7uu3VVX1NVV7GIDRQVBihRBKdIDoaT3Okkm0+/vjxsjWACRZCbJ+TzPPMmUO3PmDnM5Ofd9z8s1b25m1ksbyCmTZuI/kFxHCCGETAOU4fF+UVa4n28++g9zS09nXGYs944PIy09y99hCSGEEH4n0wA7AZ8PtrwGSYMh4RSwlUBIHGh1P3vomtxK7v5kF/nVdi4e3o25EzMJsxjaP2YhhBCiHcg0QBHQ4lMyOOfvT/LSxYOJKl9L0vwRbJh3jUwNPA5ut5vy8nLcbre/QxFCCCHEL9FqYfAVaqFKUeCt8+Gl8ZD/zc8eOqZnDEuuH8utk7NYsKmQ059YxfubCvH5utz55FaS6wghhJBilfAbjUbDpD7xPHTD1Wzufg2nlH2E++mBfPfRv/F5ZWnnX1NVVcULL7xAVVWVv0MRQgghxLFoNDD9adDq4dUz4YMrob74iIcY9VquHteDFTeNZ0SPaG75YDszX1jPzuJ6/8TsZ5LrCCGEkGmAMjw+YJQXHaBwwa0Mti1nXvB1DJ91KwNTI/wdVsBxuVxUVlYSExOD0Wj0dzhCCCFOMpkG2En5fLDtbVh+DxhD4G+bf3FaIMCGg9XcvWgX+yoa+OPQVG6ZlEVEcNf5P19yHSGE6LyON0+RYpUkcQFnz7fLuHM9bC1zcU/Pg0yfdjbR8an+DksIIYRoF1Ks6uQc9VC9H5JOVXtZlXwPWVPUEViH8Xh9vLEhnyeX7UMDzBmbzmWjuhNi0vslbCGEEOJkkJ5VosPqPXQiH/z9DB4+qydTCp/EPG8oG966F7fL6e/QAkJDQwOrVq2ioaHB36EIIYQQ4rcyh6mFKoDt78G7f4S3ZkLV/iMeptdpuXxUd1bePJ4/DErmma/2M+aRFby4+gDNrs7dLkFyHSGEEFKsEgFJp9Xwx5GZmP+2gV0xUxiy7ylK/nUqO9d+4u/Q/M5ut7Nlyxbsdru/QxFCCCHE7zHqBpj1FlTtg+eHw7L/A1fTEQ+JDjFxz1l9WHXLeKb0S+CxpXsZ8+gKXl5zEIe7cxatJNcRQggh0wBleHyHcGDHBlyf3EiRw8znfZ/kzqm9ibGa/B2WEEIIcdLJNMAuyN0M65+Fdc/AVcsgtvevPrSwxs6zK3L5cEsxkcFGrhvfgz8OTcVs+OX+V0IIIUQgkZ5VRyFJXMek+Lws3LiP+5YVMsq3iYuzDQw570a0OknOhBBCdB5SrOrCHPXqNEGPC778B4z8G4T/ct/O/Oomnl2xn4+2FBFjNXHd+AxmDUmRopUQQoiAJj2rRKej0eo4d0RvVtw0nj/ElDBs9wPkPjyCA9vX+zu0dlVZWcm8efOorKz0dyhCCCGEOJnMYerPmoOw51P4zzBY92/wun/20G5RwTx+fn++umk8o3pEc++nuxj/2Cre3JCPy+Nr58BPLsl1hBBCSLFKdDgRwUZO+8tz7JmyAJOvmbQPp7Lh+Tk0Ntr8HVq7MBqNpKWlyVLOQgghRGcV2wv+8i0MuhSW3wMvjoXC737xod2jg3ly1gCW3TiOYemR3LVoJ6c9sYoFmwrxeDtm0UpyHSGEEDINUIbHd2hul5PN7z2Idf8nXGP8F3ee1Z8z+8aj+cnyz0IIIURHIdMAxRFKt8GnN8Co66HPOcd8+L7yBp5evo/Pd5TRPTqYG87oyfRTEtFpJTcSQgjhf9Kz6igkiet8imoauefTPeTlbOGx8I+JOe9xkjP6+jusNuHxeLDZbISGhqLX6/0djhBCiJNMilXiZ3w+0GjUyyd/h+QhMOBPoP31/lQ7i+t5atk+vsqpoGdsCFeN6c7ZA5I6RE8ryXWEEKLzkp5VoktJjgzh5UuH8PDkRBKac4l7Yywbn72EiuI8f4d20lVWVvLss89KHwchhBCiq9Bq1UKVxwVuO3zyV3hhNOR8Br9y3rlvUhivXDaEj68bSbeoYG7/aAej/rWCJ77cS351Uzu/gd9Gch0hhBAyskrOOHY6Dnsj33/0GL32v4xZcfJ5z3uZcO6fiQjuHH0PnE4npaWlJCQkYDKZ/B2OEEKIk0xGVoljKtoEX90Lh76G7mPhkk/UYtZRHKpq4rV1h/hgcxFNLi8DU8OZ0jeesZkxZMVZA6qFguQ6QgjReck0wKOQJK5rsNVVs+uDh/hHwUAqNTH8c6CD6aePI9ga7u/QhBBCiF8lxSpx3A6sVFcOHHIlOBuhOhcSBx51k2aXl+V7yln0fTFr91fhcPuItZoYnRHNqJZLfJi5nd6AEEKIrkaKVUchSVzXUtXo5PkVe7li87lYNG5ys65m4LlzMZmD/B3aCWlsbGTbtm3079+fkJAQf4cjhBDiJJNilTghG1+EL26F7LNhwj8hJvOYmzjcXjbl1fJ1biXr9lexq0RdWblHTHBr4Wp4ehRhFkNbR38EyXWEEKLzOt48RToWik4vOsTE/511CqX9P+fAx/cwJOcxKv71KoX9r+fUs65Dpwv8RqOHa2xsZO3atfTo0UMSOCGEEEKoBl8JxmBY9S94fpjagH38HRCW/KubmA06RveMZnTPaABqmlx8c6CatfurWLW3kte/yUergYGpEYzPjGFcVgx9E8PQtvHKgpLrCCGEkJFVcsaxy8nP2UzNp3djaChkbtjT3DS5F5Oz49BoZb0BIYQQ/icjq8Tv4nHCpv/B14/DH16CjNPVJuwn0JOqsMbO2v1VrN6rjrxqcHqICjYyNjOG4emRZCeE0TMupEOsMCiEECIwyDTAo5AkTgDsOFTCoysKqd+/kceC38Q97h/0HXO2v8MSQgjRxUmxSpwUriYwtLQ8ePdPEH8KjPgLmE/s34Tb62NLfi2r91Wyam8le8psKArotBrSooLoHh1Cekww3aPVS3p0MDFWU0A1bhdCCOF/Uqw6CknixOG2fbsS85e3keXZy07TAIyT7yVz0Hh/h/WrqqqqWLRoEWeffTbR0dH+DkcIIcRJJsUqcVJ53fDVffDtS2rxasyNMOQqMFh+19PaXR72ljWwu9TG/opGDlU1caiqicIaO76Wvy6CjTrSWopX2YmhjM+MpXfCsVcelFxHCCE6L+lZJcRx6j90AsrgDWxd/jYRG/5F2idnM3/9zZx67t/pmxTm7/B+Rq/XExMTg14vX18hhBBCHIPOAJPuh+HXwupHYNndsPUtuHY9/I4WCEFGPQNTIxiYGnHE7S6Pj4Iae0vxSi1iHahsYsWK/Ty6ZC9xoSZO6xXLpOx4RmZEYdL/fAqh5DpCCCFkZJWccRSH8Xo8bP3sRe7Zk8jOOiM3R3/D0O5RZE+8jJDQiGM/gRBCCPE7ycgq0aaqD0DFHug9HRor4evHYNAlEN+3TV/W6VFXHlyRU8HyPeXkV9sJMemZ0CuWyX3iGJ8VS4hJilNCCNHZyTTAo5AkThyLx+tjRU4Fls//zqjGL3GjZ0/QYFyZ0+g57o9ERPpvSLrX68VutxMUFNThVjIUQghxbFKsEu0mbx28fxk0VUBUT+g9A/qcAwn92/RlFUVhb3kDS3eWs3RXGbtLbRj1WkZnRDO5Txxje0YRovMFdK5T3+ymoNpOfk0TtU0unB4fXp9CkElPmMVAdIiR1MggEsIs6Np49UQhhOhIArZYVV5eTm1tLQBxcXFERBzfaJV9+/YRHR1NZGTkz54HIDQ0lMTExON6LknixG9RVrifvK/fJixvCVmu3ZzreYDgtMFcmNbAyAF9iYhJaNd4SktLeemll5gzZw4JCe372kIIIdqeFKtEu/K6Yf9XkPMp5HwO6ePh/FfBZYeK3ZB06gmtJPhbFNbYWbqrjC93lfNdfg2RNHGWeQ8l8SOJi08g1moixmoizGIkIshAeJCR8CADYRZDu61EWFRr57u8Gr49pF4OVDa13qfXajDptei0GppcXry+H/+8Mug0JIVbSIkMIjkiiJRIC8kRQSRHWIgIMhJi0mM16zHptdKMXgjRJQRsseqBBx7gzTffpLS0lHvvvZcbbrjhmNu88847XHrppTz00EPcfPPNAFxzzTUsXLiQ8PBwACZPnsy///3v44pBkjhxoirLCliWp/DFrjJuzL+OfpqD7LYMwpF1Nr0mXEhoeFSbx+BwOCgoKCA1NRWz2dzmryeEEKJ9SbFK+I3XA456CI6C3YtgwSUQngp9zoW+MyG+X5sXrqoanazbW8ae3APk2s0U2zxUNDipaXL94uMtBl1r4So8yEBEkJGoECMxIWZiQ03EhJjUn1YT0SEmDLpj9+lye33sLWtgc34tWwpq+e5QDSX1DgB6xoYwtHskg9Mi6BETQmpkEOFBxtZtFUWhyeWlwuagoMZOYY2dwtpmCqrtFNbaKaptpr7Z/bPXNOg0hJrV9xHW8n4ig41qgStCLXB1iwoiIcwsRS0hRIcWsMWqH1x00UUMHjz4mMWqsrIyZs6cSXJyMoMHDz6iWDVgwACuueaa3/zaksSJk6G6opjcVW9j3b+I3s6duNHzVMq/GTD8NMZnxbbbmT4hhBCdixSrREDweSF/Pez8UC1cNdfAKbPhDy/6JRyP14fN4aHO7qKu2a3+tLvVy2HXa+0uqhtdVDQ4qW5y8tO/dCKDjcRaTYQHGQg1Gwi1GABodnmxOdwU1qgFJY9PwaDT0CcxjFO7RTC0eyRD0iKJDDb+QnS/TX2zm5I6tWjV6PDQ4HTT4PBga3ZT36y+p/pmN5WNToprm6locLZuG2rW0yshlN7xVrITQzm1m1o0kwKWEKKj6DSrAV533XU88cQTzJs372f3NTc3U1BQQFJSUsDOZxedV1RsElEX3ALcQkXxIQ6sepNvqhJ44c0tPGN+gYRwK5ZTZ5E5dDJG08kbAdXU1MSuXbvo06cPwcHBJ+15hRBCCCFaaXXQfYx6mfoYHFgJ+pZCTf56dVXBUy6AXtMh9OS2JfilXEev0xIZbPxNxSKP10eN3UWFzUllo5PKlp8VNgd1zW5szW4OVanT+YKMOoKMOiZmx9EtKphe8Vb6JoW1ycnHMIs6cup4OdxeSuqaOVTVRE5ZA3tKbaw7UM0bG/LxKRAeZODU1AhOTYtgcLdITklum7iFEKI9BXSx6vXXX6dnz54MGzbsZ8Wq+Ph4XnrpJZ588kkaGxt5+umnufTSS3/xeZxOJ07nj2ckbDZbm8Ytup7YpO7EXngXi4D9FQ2UfbqaxKKFJH31GY7lBnaberEq6y6iUnvTO7iRbnFRhEXFndBr2Ww2vvzyS1JSUqRYJYQQQvIc0fZ0BsicdNh1E1giYMnt8PnNEJEGg6+AUdeDxwX1hRDeDXQn9qfGycp19DotsVYzsdaO3TbBbNCRHhNCekwIp/f+MX9scLj5vrCOTXm1bM6v5T8r9tPk8mLQaeibFMbgbhEMTovk1G4RRIeYTmpMiqLQ4PRQ3ejC1zJ8TavRYDWrDeaPZ7qlEEIcTcBOAywvL2fChAm8/fbbmM1m7rjjDjIzM7nxxhuJizvyj/w1a9Zw5plnkpeXR0xMzM+e65577uHee+/92e0yPF60JcXn48CO9VTtWomxeCP3ey/n+3ozT+mf4xzdeqoJo8yYSqO1B9XdZxCaNZ5+yWG/6UybEEKIzue3TuOTPEf4TVM15H0NBRvUgtXwa6F0O7w4BnRGiMqA6EyIzYZxt7Z5v6uuzuP1kdPSa2tTfi2b837stdU9OphTu0UwuFsEA1MjSI0MwmL85dFXTU4P5TYH5TYnFQ0Oym0OKmxOyhucLb+r9zW7vb8aS4hJT3SIkRiriVirmZiWJvmtlxATsVYTUSEmWS1RiC6mw/es+vbbb7nkkktar5eWlmIymbjuuuu45557fvb4lJQUFi1axKBBg3523y+dcUxJSZEkTrQ7h9tL8f4d1B3chKt8L8aafUTa83jRfSbvuscyVbeRm02LqIwchKnXJLJGTMcSbPV32EIIIdrRby1WSZ4jAoqzEYq+g6p9ULlX/elxwlXL1PtfHAehSer0wszJEJnu33g7ueK6Zjbl1agFrLxacsps/LBYYWSwkTCLAb1Wg0YDtmYP9c3unxWhQkx6Yq1qo/q4UDNxoWZirervUSFG9Fp1FJXXp9DgUHuI1dvdVDU6qWhwUtlyqWhwUGs/srm8VgORwaaWotZPilktDfJjrCZiQ80EG3WtvbkURcHW7KGy0UnVD5eGlumeP7xmo5PaJjdOjw+Xx4vL68Og0xJs1BNk0hEdbCIx3ExiuNrAvleClV7xVoKMAT35SIgOL2CLVXV1dZSVlXHrrbfSt29fLrnkEtLS0jCbzZSVlaHT6X5xdNRll11G3759Wxus5+bm4vV6aWpq4p133uG9994jNzf3uFZHk8ajItB4fQr51U0c+n4V1l1vkVi/hWSljGbFyPKIC2gceRun9YolLtRMdXU1n332GdOmTSMqqu1XHxRCCNG+pMG66LS8Hvj6UbXnVeFG8Logphdc8glYf5w5IblO22lwuNlT2kBRrZ3i2mYaXR7cHgWfohDa0ksrOsRIrNVMXKhaJAoxnbzijcvjo7qppXhl+7G4VNHgOKyopf50enxHbGs2aNFqNHh8Ch6vr7Xo9gOjXtta3Ipu+RkZbMCk12HSazHotHh8PuwuL01OD5UNTkrqHBTXNVNmc+D1KWg0kBYVTN+kMPonhzEgJZw+iWG/OgpNiPbm9vooq3dQ+8MiE81u6lt+b3R6cHl9uL0+vD4Fs0HX8p02kREbQs/YEKJO8pTgExGwDdY/++wz7r//fgD27dvHRx99xIcffkifPn34z3/+Q0REBDfeeOPPtktISDjiP6sLL7wQm81GcHAwp5xyCitWrDiuQpUQgUin1ai9CCZOh4nTASjY9z3FGz9mX3UI//l4B6eSw0NBb1MQOwGvtgc/W95GCCGEECKQ6fQw4U71d2cjHFwJeWshJFa97c2ZEBSJNnE8wSYDWq30PTrZrGYDQ7tHMrR7pF9e36jXkhBmISHMctTH/dAT6/CiVlWDE5+ioNdq0Ou0hAcZWotS0SEmQs36E14V0eH2sr+ikd2lNnaX2NhZXM+Xu8pwenzotBp6xVvpnxLOgORw+qeEkxEbclKmLyqKQn2zm5omFy6vD49XQVHUwlywSU+wSU+ISS9TJbsou8vD5vxaNh6sYVN+DQXVdspsjp8Vag06DWEWI1azHoNOg0GnRa/VYHd5W/99eVo2So0MYkR6FCN6RDGyRxSxoYFbQ/HbNEB/kjOOoqOpbXLx/bersW6ZR2bDBkJpohYrG0LPJKffrfSNN5NOCdHJGYRFRPs7XCGEEL+DjKwSXZKiwMqHYN8XULYDNDqI6wMXvA6R3aHmIGgNYE044cbtQvwWbq+PvWUNbCuqY1thHd8X1pFb0YiiQLBRbXr/w+izWKuJkJbiUpBRh8Wgw6couLwKTreXmiZXy3RF9Wd1y8/DiwhHE2rWEx1ial0RMyrERFSwkdhQU8sUTXPrNEqT3n+jwHw+hcJaOwU1dips6ii5umYXTU4PTU51RFuT68ffm91eXB4fLq8Pl8eH0+NDA2qhzqgjIthISkQQqVFB9E6w0i8pnPToYLSdtHjX6PSwKa+GjYdq2HCwmh1F9Xh8CpHBRoakRZARG0JyRBCJ4Raigo1EBBsJtxgIOmyK7C9xe33kV6uriW7Kq+WbA9XsLW8AIDMuhFEZ0YzOiGZYetRJHUn5awJ2GmAgkCROdGQup4M93y6nae9K9tqtPFs/iij7QZaZbgWgUbFQpw2nXh/Fv5OfJjrUzISGzwg1KhhC4zCFxxMSGU9YQg/CQkNP+AyUEEKItiHFKtHV+WoLcO9djqFsM9qpj4ExSB11tX+ZWsQKjoGQGBhzM/Q5B8p2wv7lLbfHQnA0WBOPmFooxMnQ6PSwo6iebUV15Fc3UW5raTrf4MTu9NDk+uWm81aTnmirWmCKDjERbTUSFWwiOkS9HhFsxKTXotNq0Go0ONxeGlsKPI1ONzVNbmqanFQ3uahudFHT5KK6UR1x5vYe+ed8RJCBWKuZaKuRYKO+tYimjtLSHTFi64jbDnusUX/sUY3NLi85ZbbW0Wg5ZQ3klNqO2AehZj0RwYfHoSPIpCfEqL6OxajFqNNh1Gsx6bUY9VoUaClueahqdFFUa+dQVRNFtc2A2kOtb1Joy4IBkQxKjSAsqOMtUOX1KRysbGRHcT07iuvZkl/LzhIbXp9CdIiJYemRDO8eybD0KHrGhpz0v9mqGp2sP1DNutwq1u6voriuGa0GesSEcPmo7vxpWOpJfb3DBew0QCHE71NdU8vCFZuZM+dmRiYkcJmiUFFTR87BDBorDuGpLYTGSjzOZpxehe1Fdcyq+Zjevv0YNZ7W57ncdQtbTEO5JHwHI425GJIGEJc9iuT0PmjaeNh9g62OprpKnPYGfF43eqMFQ1AYxvAELHotJr0GrU56AwghhBBdUbnDwEtLC5kz558kGIPUG6c8AjVXQ10BNFZAU6VanAKo2ANrnwRH/Y9P0m00XP6Z2ifrk7+po7QSB0LiADAGt+0bUBQ1RlcjuJtBZwC9SS2g6Y3g84FMceyQQkx6RvRQp1D9Ep9PweHxYnd50WvV6VgGnfa4ij8nwudTqGt2t6za6KSipXBWYXNQ1eii0emhpsmuFr5cPxS/PLh+0g/sp4w6LcEtRazDi11GnZbqJnVqZml9Mz5FbWeSERNC7wQrE7Pj6J0QSveoYGJDTZgNJy+fr7e72VmiFna2FtTy3neF/GflATQayIy1cmqautrl4G6RpERaAuqE/E8LUzuK6tldasPeUthLiwqif0o4s4akMiw9kvTo4DaPPzrExFn9EzmrfyKKopBfbWfjoWp2FNcTagmMMpGMrJIzjqKDaW5u5uDBg6Snp2OxHH2+/+EUnw9bfQ22qmIaa8oo0qeSU6cnZs9rjK1eQKJSAUA9wSwOv5Dy7KsYmGDklFgDUXHJvylGRVGoq66gsiAHW3EO7soDGOoO8ol2Ap/aejLDuZj7DPOP2Gadtw8Xuv+BGSc55stpVow4NGacmPBptFxtfhSbNoy5zhc41btdfR2NBgUtS8POpzjtPE6NhcEJBhJSe7Z5wU0IIdqKjKwSXd2J5jp4XGoRq6kSNBpI6A8N5fDun6B8F3iaQaOF2Gy46iswmKG+CELif/vUQq8bavOh9pA6RbH6gPqaUx5Ri1EPJaqvd7g5q9Vi2Wc3webXwBAMBotayBp0CYy9GUq+hw+uUJ+Llj9Wg6Ph8i/U2yr2qCso6v3fJFl0XG6vjyan57DRW57W0UyNTg92109vU6ftOT1eIoPVlRq7RQbRJzGMnnEhJ7Uodbx+KLBsyq9lc34Nm/Jqya1oBCDGamJwtwh19FVaJH0SQzHo2udvA0VRKKl3sLWglq0F6hTSnxam+iaF0S8pjH7JYfRJDCPM0vFGhv0eMg3wKCSJE+Ln6qvLyd+xFvvBjXzj6MYb1ZkMa17DPOO/KdHEUhbUC7clmobQnuxNuQDF6+bUglfB04zWaUPvqsfgbuBu0y3sq4XnlIeYoNsGQA2hVBiSWRd/Mfa0iWQGNZLgPIjRYkWj0+FxOXFoLVRZe+NyNBN36GN8Lju47ChuOyg+vkm8FJcumL6Vi4m2H1SDVnyg+FinH87nDekMr1nEg4b/UUY0RaEDUbqNJnHgmSR2zwqosytCCHE0UqwSog14PVCZA8Wb1OLSxPvU25/uB01VkDAAotLBEgFD50B4KuQuh6JvwWFTR2056iBrilpYOvQ1zJ+hPofOCBHd1ULUH15Sbzu4Wi2MGSzg84DHAUmngskK+d9A+U5w29WRVx4npI6AzElqAey7l4GWP9EURR0JNuFO8HnhkTT18clDIG20ekkZpo7YEqKLq7O72FJQy6a8Wjbl17KtsA6nx4fZoKV/cjiD01qmDnaLOGkFIofb2zqNb2tBHVsLaym3OQFIibTQPzm8SxemfokUq45CkjjRkdntdvbt20dmZiZBQUFt9jqKolBSXEjZti9xF2wipG4PQZ469isp3KG5HgMeFnuvwaUxYdeG4NBbcemtLEu7ifDYVPpo84izGolLyyY0vH2Wna6vqeDQ5uU49n9NVNV3pHsO8IlvJI8F3cRpKTDNt4qg5H7EpPcnNjENvVHOSgohAo8Uq0RX1165DooC+euheLNaxKovhuYamPWmOm1wxQOw9S0wh4I5DMzhkHUmDL5CLV6VfA8RaRCWDNp2GFni80HZdjXm/HXqSoqOOrgxB0ITYNu76om82Gx19JVZvv+ia3N5fOwsqWdzXi2b8mvYnF9LVaMLjQZ6xYcyNC2CId0jGZoWeVyr4v0wne/7wjq2F9XzfWEde0pteHwKFoOOU5LDGNQtgoEp4QxMjSDGKn9r/BIpVh2FJHGiIystLeWll15izpw5JCQk+DucgFZfW8WOg0WsLjPSlLuWf9TeRbDGAYBX0bBPm85dsc+RGG7hD03vYbZY0FvjMIUnEBIVjzWpFxGhobJcsBCiXUmxSnR1kuscJ58PqvZCbG/1+hvnwoEVP95vCoM/vKiOBstfD3nr1Mb0wTEQHAvhKWCN90/sQvjBD1MHv82r4btDNXyXV0NetR1QezilRweTHGkhxKTHYtDh9irYHG7q7C7yq+3k19hb+331iAmmf3I4A7tFMCg1nKw4K/p2mmrY0Umx6igkiRMdmaIoKIqCRqORqW2/keLzUVaYS+XBHTiqC7DZm/ncPI2SWjsPl80h1lfZWswCmOZ8iD2k8Q/zh5yh+Ran1oJba8GtD2J76AT2xE4hkUoG1y0FkxWd2YrOEobeGoWm+1jCLAbCtU6soWHSMF4IcdykWCW6Osl1fgdnA1Tuhdo8sBVDr+kQ1QM2vABrHlenPP4wxbDvTJj5itq367VpYLSCKUSddmgMhpmvqY3gN78GzbXqFEajVf2ZNEgtdLlb8ibDsUelCBGIKmwOvsurZV95A3nVTRTXNtPk8tLs8mDQaQkx6wm3GEiNDKJbVDBZ8Vb6JYcRapbpfCdKVgMUopOSxO3EabRaErplkdAtq/W2M1p/2wmAvbGe2ooSGqqLucHYgwqHlvBDBZTWGtC6m9C67ei9duwOB7tKbDQ17WV28wKCFTtBGnV+epESzWjnMwBsNF2HlXrqNUE0aKzYdaHMj5pLfVgvhtlXk9n8PYrOiKIzo9EbqQw7hfK4MYRomsm0byUsPp3o5ExCI6L89rk3OtwU1zmoO/AdSsEGqCvA2FSCwW1jg34Ib2umkercz2PuB9GgtF7qNaHMsf6HIKOOG5ueIgIbHl0QXn0QPkMwuxP/QHN4Jgm1m0io3aT29PB50fjclIf0Ijd2ChZcDLYtxRyVTHhCT5Iy+qE3SF8OIYTozCTX+R1MVkgerF4ON/wa9eLzgr1aXa1Q31Jg0hog+xx19UJXU8sqho4fVyzc8ykUbVILYYraJJrzX4M+58KmV2DpnaC3qP2+LBHQ8wy1J5ijHr68S20GrzOqr6c3wagb1B5b+evVqZjhKepKib+1yf3Joihq83pno/pe6wrUS1MlOG1w2WfqNM8Fl0Dht+pUS59XXeVx0gPQbybkLoNvngPjD8W+EHXE29A/q43/1z7Vkue0XFBgzM1gCYf9X6nFwLBkiMlS96FoN7GhZqadksA0ZBRnoJFilRAdTE1NDUuXLmXy5MlERkb6O5xOJygkjKCQMEjvTa8fbhx27c8eNwC4BoDRwHUAeD0eGhvq0DY1sFgTia3ZTeGBB8lvqMBrr4XmWnTOOjTmUGqaXHjqi4hq3oFOcWNQXBgVF7uLqnlqexzZ3hzeN9zd+noNioVCfTeeSFGnLZ7m+JJQayhB0SmYQ8Ixh4QTFJmMxWJBpwGthiNWRFQUBY9PweVy4W6qw+1y4HY6aG6sw95YT1HIKdQ1u0nKfQujLR9jYzGhjmKiveVc5byRTUovbtG/y1W6L6jQxlBvjMVlCCMsPJozEmKJ9uk4UHGe2ky25Q8Mt6JlbEwMzW4PFIbjc9oxuusxOsow+pqZXz+AjV4DFyvrGe37FC86vBodClryq5v4uPgU4txFXON5AJ1GPQvcrBjZb0jnzewXyE6KYlBIDWnpPTFb2ngZciGEEO1Gcp02pNVBSKx6+YE1Dibe++vbXPSh+lNR1EbxDptaFAPoOQmCotRiS3Od+jM8Vb3P41IbyXuc6nYel/pz1PXq/cvvgcKN6u8aHYQmwZkPQ+/pULZD7SUWlgyWSLVnmCUCgiLVOKA13/gZZ8OPjevdzWpfr/BUdSRYyfew93N11FltvlqUSugPs94AnxsWXvPjFMmQOAjvBl4XaC3qe43N/jHX8XrUUWugFuPMYeCyg71GLfh5WkadKT7Y9D/Q6tX9r235E3xky37Y9D/IWfxj/OGpcMY90Pc89bm8LpmuKbokmQYow+NFByMJXNeg+HxUVxRTXZRLY8Uh3FV5NNmbeMv8R4prm3m17nISNVVHbHOe8242K1n8U/8GV+m/wKdo8KHBh5a3vadzj+dS+mkO8qnpn0ds16BY6Od8BY0GPjLdS7SmgTpjHM3ByfhCU2noeTaRyZmkWDVEh1nbfUqjx+2ipqKYirydNOZtpbm6mH95L2R/ZSNrDX8hmnoKdSnUhGTgCU2hvOefiEzsTmqQm/ioMExmKWSJjkWmAYquTnKdLsJlh/pCqCuE+pbRTH3OVYtHG+bBkjtonbIIkDUN/vg2NJTBEy2j5DVatdCFAncUqasvvjYd8tYc+VrTn1Ib429fAMv+D0IT1aJQeCokDlRfV1HU4paxDZv6/xqHTd0X5bvUJvq9ZkDqsJb9cLtaQEs4BSJ7QLeR0OccdXRXc61aLJSRiKIDkZ5VRyFJnBCio1N8Pmz1NdSV5dPcVIe7qZ7SkD7YCCa8+ntCG3JRfD51uL7iwxbak7qYoQT5Gkms/RatwYRWb8EUEk5weDQh8ZlYzXq0HaiZvMPloWDnWmr3b4KybYQ27CfSXc6fXLdzwJfEQ/qX+ZN+BZVEUG2Ip8mSSHHCJMynnEP/5DDiwiz+fgtC/CIpVgkhBOB1g61EHRnlbFBHcyX0V6cq7vzwx+l4itrwmkGX/ji9sLlO/d0QpK7iGJbcMVdHbCiHom+hdLtaxKrNg7QxMO1xqD4Azw4CQ/CPhbfwVJj8kPref5jeKESAkWLVUUgSJzoyaToqxNG5PD7K6h3UHNiEq/h7vDX56G2FBDeXsMAzltfsozhDu5kHjPPJjxqFKXsqvUbOkKmEImBIsUp0dZLrCHEcHDY4uOrHHlt1BeBugks/Ve+fNxqCoyDzTHVFyIg0f0YrRCspVh2FJHGiI5PlnIU4cYqiUFrv4MCu7zBse4Pkyq9JVsqwEcT6+IuIOvMOBneLkD+OhF9JsUp0dZLrCPE7KYo6hXD/Mshbq/a9ShkG58+HUPlOCf+S1QCF6KTCwsI4++yzCQsL83coQnQ4Go2GxHALiaPGwqixKD4f+fu+p2Ttm2wuM/HfF75hbHgV18ZsJ3HUhXTrNcjfIQshRJcjuY4Qv5NGAyOuUy/ORti3BPYt/bGx/upH1dUKMyaCwezfWIX4FTKySs44CiGEAHw+hW/zajiwcj5nFzxKiKaZQ9o0ylLOJGnMxaRm9PV3iKKLkJFVQggh2ozHCa9MhNJtYLRCr6lqg/mMM0Bn8Hd0oguQaYBHIUmc6Miam5s5ePAg6enpWCzSIFqItuBobmLPmo/x7viIbNta3vSewdthf2ZGdw1nRFaQOWwKlmCrv8MUnZQUq0RXJ7mOEO2gch/s+gh2fqQ2br8tT10Jcd9SiO4Jken+jlB0UjINUIhOqq6ujg8++IA5c+ZIAidEGzFbghk46SKYdBHNTQ303FvIqHwv7l1v09/zPM7VBnaY+9HU7XRShs8kKb2Xv0MWQohOQ3IdIdpBTCaMvx3G3Qb1hWqhyueFj6+G5lq1WJUxEbLOVFcglFFXop3JyCo54yg6GJ/Ph9vtxmAwoNVq/R2OEF2K4vNRkLud0s2fElSwkl7N2/jUN5yXIm9jamYQZ8bV03PgeLQ6nb9DFR2YjKwSXZ3kOkL4kbMBDq1Rm7PnLlMLWTftBWs8FG+GiO4QFOnvKEUHJiOrhOiktFotJpPJ32EI0SVptFq6ZQ2gW9YA4C4abbVE5OTTJw9KNn1MlvI0VYvDORQ+Al3WZDJGnEVoeJSfoxZCiI5Fch0h/MjU0seq11R1VcHqA2qhSlHgvUugoVRdWTBzEvScrDZql1WURRuQkVVyxlF0MLW1taxYsYLTTjuNiIgIf4cjhGjhcbvZt2k59dsXk1D+NWm+Alb6BvBi8r8YnmxmnH4XcT0HE5ucjt5g9He4IoDJyCrR1UmuI0SAspW0rCz4JRxaDW47XLcRYntB/jegN0FUDzDLSp7i18nIKiE6KZ/PR1NTEz6fz9+hCCEOozcYyB4xBUZMAaAkby+N+woIKTWwc/MabvDeBevBo2gp0UZTYkhjXuJDhFkMTK97E4setMFR6KzRmENjMaQMIiUxgSCj/FcthOhaJNcRIkCFJsLgK9SL2wEF30BMlnrfl/+E4k3q7+ZwCE+BM+5RVxnMX682bg+KhKAo9RLeDeKy/fVORAcgGbAQHUxUVBSXXHKJv8MQQhxDYloWiWlZzAAU3yDKiqdSeXA79oqD+GrycbhcaDRQWGsnpnI9cd5SwhUbRo0HgPOcd7OFLOaGLGe0MRd3t7GkDD2bxLQs/76xTszn9VLX7KHG7sKevxVvTR6exmq8TTUojjq2W4axXZvNExf0x2yQvmRCtBXJdYToAAxm6DHhx+uXfgoVe6D2ENTlQ10hmFtGRtbmwe6FauN2R716W9Y0+OPb0FwHb50P3UZAj9Og2yhp5t6WfD7QasFeA0WbWj6TOvWnRgfjbvF3hK2kWCWEEEK0MY1WS3xKBvEpGUfc/mOKtx5QG7g3NtZTX1XGP5RwDtR6MO/ag6VoI/12/QvD7gfZbehL2aAbGTLhLKxmSeZOhMfrY3dBOfb1L6GtzCG08RDRnlIilHpOd86jllD+a3icibotANgIolFjZUNQFPURPXG6fVKsEkIIIQ5nDILkU9XLTw34k3oB8LrVQoniVa+7m9VRWN+/A+v+DcExcMosmPSA9ML6vQ59rY5oq8qF6lxoKIehV8HE+6AyB94+X32czgSWCHUFyAAqVknPKunlIDqY0tJSXnnlFa688koSEhL8HY4Qop002mrJWb0A4673eNQ2mU26U7iuRxVnZEXTe9hkNLJi1lE12mrJWfMBxQd2c2fVZBwuFxtMf6XWEEd9cBqesDS01jhsGWcTHhFFjKae8GAL1vBodPr2PbcnPatEVye5jhBdkKJA6TbY8T40VcIfXgKfF7a8Dn3OUYsp4tf9sP92L4L+f4SYTFj5EGx7F2J6QXRPdRpn8hBIGapO47RXq1MzDZZ2DfV48xQpVkkSJzqYpqYmdu3aRZ8+fQgODvZ3OEIIPyird/DhliJ6rr2RSd7VFGkSKEr7A4kjLiAl4xQpXLWwNdSzb+Vb6HI+JbvpO0waN7v12awa+RrDesTSL9GK0RB4g8ylWCW6Osl1hBAAlGyF/54OWj30ngH9Z6vTBI1B/o4scJRugx0fqEWqunywRMLZz0GvaT9O+QswUqw6CknihBBCdAY+r5c9G5Zg//Y1+tatwqJxcbv+VmpSJzM+6BA9fHkYQmMwBkeg0epxB8fjDu+OxtmAqeJ7fO5mfG4nXrcTt1fhQNyZOD0+ehQswNRcBR4HGq8TjdfJ5uhzaIrsTXdnDinOAxhD1UbwQRFxhMYkExoRg07r/+H69dXlbNn0DW+UJrEt9xDr9ddyyNCTuu5T6TZ6FgndAr/nlxSrhBBCiBYN5bD9XdjyhjqVLXkIXLVcHXW1YR4ER0NQNOiNoNFC6gjQ6tQiTlMleJzqxetSt43qAaXb1elxXid4HOr9Yckw6np1hNLGF1oawUeqzx0UpY5K0gZACwCfD4q+U0dKBUXCZzfBroVqMS/7bEgbA7rAOxF3OFkNUIhOyuFwUFBQQGpqKmaz2d/hCCH8SKvT0WfUNBg1DXtjPdu/W0ZcQwrFlV7YvYjB3sXoND+ek3rNM4l7PJdxiuYAn5juOuK5bEoQF7lTMOm1vKP9gASqcGmMeDRGPBoDuc5RfJsfxiznUs71vXPE8y70juRGz1/pZalnHg/SpAvHYYzAbYrAZwrj2/S/E2TSk1m7iiAcGK3RhCWkE5OcQbA1/Hfvh6qyAg6seR/Lgc/p3byVAVh4Mf49/jp1CLUZO+gVl/i7X0MI0X4k1xFCtLLGqUWkkX+Hqn1qQ3aA+iJ1mpu76cjH/6NcLSotuRPy1x5534x/q8Wqyhz49iXQm0FvUn+mDFUf47bD8nvUItbhbsyB0AS1OFS8uaWYFa0WjHqfpTaIrzkIeevAFAKhSRCWAiFxv390k9cDBethz2LY8wk0lKrv5dTL4PT/gymPBkYh7SSTkVVyxlF0MKWlpbz00kvMmTNH+jgIIY7K5/XSUFdFk60WRfHi1gXjsUSD10mQoxKjyYLBZMZgMmM0WdAbjMf1vF6PB1ttJbaacppqy6jxmCkwpOOsLab3wdfQO6oxumoJctei9bmZpX+KZpeX95Rb6KPNP+K57tTdSHnKVMaHldE3xEZ85hDiUzKOOpXR5XSwN3cva6qC2bRjDy9XXYwC5Jj60dhjOj3GzCY6sdvv2XV+JSOrRFcnuY4Q4ri5m9XeS143KD6I6K4Wh2ylahN3vRl0RrUopTMef9N2lx3sVepzN1VD+nh1xNLm+Wqxyl7dcqmB0XNhwB9h+wL46M9HPk/iIJizUh2xtfVNiO+r9pA6Vp+o+iK1T5cxGD6+Bra9oxbAes+A7HMgZVhATvE7HjIN8CgkiRMdmdfrxW63ExQUhE7X+SroQojOS/H5cDqbqasqpbbkAE3lB/lW04+NlUZOL3yWS5RPAGhULBQZurEl+mz2JZ5NnKeYQSXvoHPWEm4vINlTQKkSyVTlGcb0jOayyJ30HjaZiJjO8UetFKtEVye5jhCiw/J5wWmD+mKoL1QLZBmnQ0MZPNlbLaihgcjuatHq/PnqFMb1z6kjvuoK1NX7Gkpg1lvQe7o6bdHngcSBnWKFRJkGKEQnpdPpsFqt/g5DCCF+M41Wi9kSTHxKBvEpGQAMBq4DFN98yooPUp67ieaiXeiqc7A1O9hwsJpkZxGTXVto1oVSE5ZNVdQ5RPYex/b+YzDotC3PIoToLCTXEUJ0WFqdOiLKEqGOovqBNR7uKIaK3VCxR73UHACdQb2/bLtapApPgf6zIGkwpI1S70s4pf3fRwDo0sWq77//npCQkF+9PzU1lejoaKqqqigoKGDQoEEA7N27l6ampl/dDqBv374YjUYOHjyI1+ulZ8+eeL1etm3bdtTtzGYz2dnZrfElJCQQFxdHbW0thw4dOuq2sbGxJCcn09jYyL59++jduzcWi4X8/Hyqq6uPum1GRgahoaGUlpZSVVVFv379ANixYwdut/uo2x6+X8xmM926daO5uZk9e/YcdTur1XrEfunevTsRERGUl5dTXFx81G2TkpKO2C/9+/dHp9ORm5tLQ0PDUbc9fL84HA6ystRmu1u2bDnqdgaD4Yj9Eh0dTUJCAjabjf379x9126ioqCP2S2ZmJiEhIRQVFVFRUXHUbQ/fL6WlpaSlpfH1118TGxuL0Xj0KTuH7xedTkd6ejoul4udO3cedbvg4OAj9stPvwtHk5CQcMR+Ofy7UFdXd9RtD98vNpvtiO+Cz+f71e20Wi0DBgwAYPfu3YSGhh7xXTia8PDwI/bL4d+F0tLSo24rxwg5RhwuUI4Rh38XHA7HUbcNpGPEobw86upsEJ6JMTwTgKEtFxhAHa/+7BihFqo63zGisbHxqM97vCTP+ZEcwzrWMayuro4333yTfv36HbVoFUjHMMlzjiTHCDlGHE7ynMOPEVqgD8T2gVhg61Z1w7RryJz0k2OEJQLofMeI485zlC6ovr5eAY55+e9//6soiqL897//VQ7fVcOHDz/mtoWFhYqiKMrMmTOVSZMmHffrZmdnt76O1WpVnnjiCUVRFGXBggXH3Hbu3LmKoijK+vXrFUDZuXOnoiiKcuWVVx5z2yVLliiKoih33323kpSU1BpDUlLSMbc9fL9ceeWViqIoys6dO4+53U/3y4IFCxRFUZQnnnjimNv+dL/U19criqIokyZNOua2h++X4cOHt8Z/rO1+ul/uvvtuRVEUZcmSJcfc9qf7Zf369YqiKMrcuXOPue3h+8VqtSqVlZXKyy+/rGRlZR1z28P3y8yZMxVFUZTCwsJjbvfT/fLT78LRLj/dL4d/F4617eH75affhaNtZ7VaWx+bnZ39s+/C0S4/3S+HfxeOta0cI+QYcfglUI4Rh38XjrWtHCMC9xhx+OfzW0me8/OLHMM61jGssrJSSUxMPO7viBzDfv0ieY4cI3566QzHiB++C8faVo4RgXuMOPzz+TVdumfV6tWr5YxjCzmbIGcTDidnHOWM40/JMUKOEYeTY0Tbn3EcN27c7+5ZJXnOj+QYJseww8kxTPKcn5JjhBwjDifHiMDIc7p0sUoajwohhBAi0EiDdSGEEEJ0Vsebp3TMtQ6F6MLKysp45JFHKCsr83coQgghhBAnneQ6Qggh2r3BemlpaeswzISEBKKioo5ru927dxMXF3fE430+H4cOHSI6OpqwsLA2iVeIQBMSEsLo0aOPOrVDCCGEEKKjklxHCCFEu4+sevXVV5k9ezajR4/mjTfeOK5t3nzzTQYOHMirr77aelteXh69evVi9OjRJCYm8vzzz7dVyEIElJCQEEaNGiUJnBBCCCE6Jcl1hBBCtHux6s4772Tnzp1Mnz79uB5fWlrKiy++yB/+8Icjbv+///s/Jk2aRGlpKd999x133nnnMRvnCdEZOJ1O8vLycDqd/g5FCCGEEOKkk1xHCCFEwPesuvbaa3nyyScxmUxH3P7ZZ59x7bXXApCdnc3gwYNZvny5P0IUol3V1NQwf/58ampq/B2KEEIIIcRJJ7mOEEKIdu9Z9Vu8+uqrZGdnM2TIkCNudzqd1NTUkJyc3HpbSkrKry7x6XQ6jzgzY7PZ2iZgIdpBTEwMf/vb32SFJyGEEIDkOaLzkVxHCCFEwI6sKisr49FHH+W8885j586d1NXVUVZWRklJCQaDAZ1Oh8vlan280+nEYrH84nM9/PDDhIWFtV5SUlLa620IcdLp9XoiIyPR6wO61iyEEKKdSJ4jOhvJdYQQQgRssaqwsBCdTsell17K7NmzWb16NW+99RavvPIKWq2WtLQ0duzY0fr4HTt20KNHj198rjvuuIP6+vrWS2FhYXu9DSFOuvr6er744gvq6+v9HYoQQogAIHmO6Gwk1xFCCNHupytqamooKSmhvr6esrIydu7cSUZGBmazmZKSEnQ6HXFxcQwZMoSdO3e2bnfZZZfRt29fbr75ZgAuv/xybrvtNh5//HG+/vprGhoaOP3003/xNU0m0896XgnRUblcLvLy8hg8eLC/QxFCCBEAJM8RnY3kOkIIIdq9WLV06VIefPBBAA4dOsTixYt577336NOnDy+++CIRERHccMMNP9suOTmZ6Ojo1uu33norTU1N3HTTTSQmJrJ48WJ0Ol17vQ0h/CYmJqZ1cQEhhBBCiM5Gch0hhBAaRVEUfwfR3mw2G2FhYdTX10vjRiGEEEIElN+bp0ieI4QQQohAdbx5SpfsWvhDfU5WyxEdUUVFBQsWLOCCCy4gNjbW3+EIIYQ4yX7IT070fKLkOaKjk1xHCCE6r+PNc7rkyKqioiJZKUcIIYQQAa2wsJDk5OTfvJ3kOUIIIYQIdMfKc7pkscrn81FSUoLVakWj0bTJa9hsNlJSUigsLJQh+AFCPpPAI59JYJHPI/DIZxJ42uMzURSFhoYGEhMT0Wp/+8LN7ZHngPz7DDTyeQQe+UwCj3wmgUc+k8ASSHlOl5wGqNVqT+hM5YkIDQ2VL12Akc8k8MhnEljk8wg88pkEnrb+TMLCwk542/bMc0D+fQYa+TwCj3wmgUc+k8Ajn0lgCYQ857efrhNCCCGEEEIIIYQQoo1IsUoIIYQQQgghhBBCBAwpVrURk8nE3Xffjclk8ncoooV8JoFHPpPAIp9H4JHPJPDIZ/Ij2ReBRT6PwCOfSeCRzyTwyGcSWALp8+iSDdaFEEIIIYQQQgghRGCSkVVCCCGEEEIIIYQQImBIsUoIIYQQQgghhBBCBAy9vwPojLxeL2vWrMFutzNu3DiCg4P9HVKXtm/fPrZs2QJAZmYmgwYN8nNEAiAnJ4cDBw5wyimnkJKS4u9wBLBt2zYKCwsZOHAgSUlJ/g5HtPj222+prKxk2rRp/g6lS8vNzWXz5s2t10NDQ5k6daofI/KvQ4cO8f3339OrVy969+7t73C6vPfffx+v1wvAzJkz0eslxfe3hoYG1q9fj9VqZfjw4Wi1MkbA32pqatiwYQORkZEMGzYMjUbj75AEYLPZ+PzzzxkzZozkn352+P8lAGPHjiUxMdFv8cj/ZCeZ1+vl9NNPp7y8nKioKP7617+yceNGYmJi/B1al3XgwAEWLlzIzp07GT9+vBSr/ExRFP74xz+yfft20tPTWbNmDY8++ihXX321v0PrsjweD+eddx4lJSXExsaydu1annnmGS699FJ/h9bl5efnc/7552MwGKRY5WdLly7lmWeeaf0/JCkpqcsWqz744AOuuuoqhg8fzqZNm7j//vu59tpr/R1Wl/bpp5/icDh4//33mT59OiEhIf4OqUv7/PPP+dvf/kZGRgYFBQVYrVZWrFghn4sfffzxx9xyyy307t2b3NxcIiMjWblyZUA0ke7q5s6dy4cffshrr70mxSo/u/jiiznrrLNai+tZWVl+LVZJg/WTbNGiRdx5551s3boVo9HIRRddRI8ePbj33nv9HVqXd88991BVVcVzzz3n71C6NJ/PxwcffMAFF1wAwBdffMEVV1xBaWmpnyPrupxOJ+vWreO0004DYP78+bzyyit8/fXXfo6sa1MUhWnTpnHuuefyyCOPsH//fn+H1KU999xz5OTkdPn/QxRFISsri0cffZRzzjmHrVu3MnHiRIqKijCbzf4Or0tzOBxYLBYaGhqkKOJny5YtY9CgQURFReH1ehk0aBC33norF154ob9D67JWrlzJ8OHDsVgseDweUlJSWLRoEUOHDvV3aF3a559/zscff8yePXu4+eabOeecc/wdUpdmNpupq6sLmP/PZWTVSbZ69WrOOussjEYjAOeffz5PPfWUn6MSInBotdrWQhVAUFAQYWFhfoxImEwmTjvtNBYtWkRVVRXz58/nkksu8XdYXd68efMYOXIk/fr183cookVZWRkLFy4kLS2NAQMG+DscvygtLaWwsJAZM2YAMHDgQCIiIti5cyeDBw/2c3RCBIaJEye2/q7T6TCZTJLr+NmECRPYu3cvmzdvZseOHXTv3l3+f/Wzuro67r//fpYuXdplRyoHomXLlmGxWBgyZIjfj1syefokKy8vJzY2tvV6bGwsZWVlfoxIiMDV0NDA3Llzuf/++/0dikAd5fb+++9TVVVFz549/R1Ol3bw4EHee+89br/9dn+HIlpkZmai1+uZP38+06ZN46yzzjqir0NXUV5eTmRkJDqdrvU2yXWE+HXz5s3DaDQyZcoUf4fS5eXm5vLRRx/x4YcfMmTIkNbBBcI//v73v3PvvfcSGhrq71BEiwsuuIC33nqLu+++m/T0dNatW+fXeGRk1UkWFhZGY2Nj6/XGxka/VySFCET19fVMnTqVP//5z5x//vn+DkcAL7zwAqD2Prnssss4dOiQnyPquv76178yZswYPvjgA/bv309jYyPvvvsus2fP9ndoXdakSZOYNGkSAHa7nf79+7NkyZIu10vsp3kOSK4jxK95+eWXee2111i6dOkRBV7hH9OnT2f69Om43W4GDRrEwoULOe+88/wdVpe0YsUKtm/fTk1NDe+++y5VVVWsXbuW3r17k5WV5e/wuqzXX3+99fennnqKf/7zn6xcudJv8cjIqpMsOzubDRs2tF7fsGEDffr08WNEQgSe6upqzjjjDC6++GJpyhsAamtr8fl8rddTU1NpaGjwY0Sid+/e7N+/n4ULF7J69WqamppYuHChv8MSLYKCgkhMTMRms/k7lHaXnJwMwJ49ewB1Gsf+/fvljwshfuLZZ5/lv//9L0uXLiU8PNzf4XR51dXVrb8bDAbi4+Ml1/EjjUZDr169WLhwIQsXLqS6upqNGzeSm5vr79BEi4yMDL/nOdJg/SSrra0lIyODyy67jISEBB544AG++uorTj31VH+H1mWVlZWxatUqPvjgAxoaGrj88ssZOXIkqamp/g6tS3K73QwePJikpKQj+iLNmjVLlhD2k2+//Zbbb7+d8847D5/Px7x585g6dSqPP/64v0MTqCc9LrroImmw7mdfffUVlZWVOJ1OVq5cyZdffsn27duJjo72d2jt7pZbbmH58uVcffXVLFiwgKSkJN544w1/h9WlLV26lIqKCi655BJee+014uPjmTx5sr/D6rJefvll5s6dy+OPP9466rB///707t3bz5F1Xeeeey7Z2dlkZGTwzTff8Mknn7Bt2zbi4uL8HZoARo8eLQ3W/ay4uJg1a9YAUFhYyNNPP82tt97K9ddf77eYZBrgSRYREcG6det4/vnn2bt3LwsXLpRClZ+Vl5ezcOFC9Ho9ERERLFy4kG7dukmxyk/cbndrsnb4SJELLrhAilV+MnToUB577DFef/11vF4v9913nwyLDyDR0dFMnz7d32F0eWvWrCEnJweTyUTPnj3517/+1SULVQAPP/wwqampfPPNN0yaNMmviaxQrVy5kry8PGbNmsUXX3xBWlqaFKv8yOfzMW3atCOmzwQFBUmxyo/efvtt5s2bx9q1a+nevTtbt26VQlUAOeOMM1pH7gr/KC0tZeHChWg0GmJiYnj55Zf93mtPRlYJIYQQQgghhBBCiIAhPauEEEIIIYQQQgghRMCQYpUQQgghhBBCCCGECBhSrBJCCCGEEEIIIYQQAUOKVUIIIYQQQgghhBAiYEixSgghhBBCCCGEEEIEDClWCSGEEEIIIYQQQoiAIcUqIYQQQgghhBBCCBEwpFglhOj0cnJyuOuuu1qv5+bmcuedd7bpa27evJmZM2fypz/96YSfIy8vj5kzZzJz5kzq6upOXnBCCCGE6FTuvPNOcnNzW6/fdddd5OTktOlrXn/99cycOZMlS5ac8HPcd999zJw5k9dff/0kRiaE6AykWCWE6PSqqqpYuXJl6/WoqCgmTpzYpq9ZWlpKQUEB559//gk/R0REBLNnz2blypU4HI6TGJ0QQgghOpMVK1ZQXV3dev30008nJiamTV9z6dKljB8/np49e57wc0yYMAGr1cr27dtPYmRCiM5A7+8AhBCirT344IPk5OQwc+ZMevTowVVXXcWyZcuYMGECu3btYsGCBfTr148vvviCzMxMbrrpJl555RXWrVvH+PHjueKKK1qfKz8/nxdeeIHS0lKGDBnCnDlzMBgMv/i6iYmJnHvuuQCtrzNw4EAWL15MZmYmc+fOxWAw4Ha7ee6559i8eTMOh4MJEybwl7/8hbCwMGbOnMkNN9zQHrtJCCGEEB3QggULWkeNR0ZGcu+99/LVV18RHx9PVFQUt912G1OmTGHp0qVUVVUxd+5cAJ599lkMBgN33HEHCQkJAPh8PubPn8/atWsJDQ3l6quvplevXr/62meeeSY9evQA4LbbbmPGjBksXrwYm83GnDlzGDBgAAAbN27ktddeo6qqCkVRePXVV7FarYwZM4bNmzdTVFTUtjtJCNHhyMgqIUSnN27cOGJiYpg9ezaTJ0+murqaFStWAFBZWckTTzzBkiVLGDt2LK+//jrDhg1j3759jB8/nrvvvptly5YBUFRUxNSpU4mNjWXixIl89dVXx11Iqqys5KmnnmLBggWMHj2aJUuWcO211wLw1FNP8dlnnzFlyhRmz57NsGHD2mQ/CCGEEKLz6du3b+uo8dmzZxMbG8vKlSupqqoC4KuvvuL6668nPT2dkJAQpkyZwt///neGDh2KzWbjsssua32u6667jjVr1nDaaaeRmJjI5MmTKS0tPa44vvrqK/785z/TrVs34uLiOO200ygqKqKhoYEZM2bQq1cvZs2axezZszGZTG2xK4QQnYiMrBJCdHqjR49m8eLFzJw5E4ANGzYccX9ycjIvv/wyAE6nk/fee48nnngCUPtGffPNN0ycOJFXXnkFgHXr1rU+dtGiRfznP/85rjhCQ0N544030Ol0nH322SQmJjJv3jwaGhqIiYkhOzub/v37o9XKeQQhhBBCHJ/s7GwiIyOZMGECw4cP/8XHPPDAA8yYMQOA//3vfzz22GMMHDiQ8847j+TkZAAaGhp49dVXmTFjBh9//DGgjrRavXo1s2fPPq5Y/vnPf3LhhRcCUFJSwjvvvMMVV1yBRqMhLS2NsWPHEhER8XvfshCiC5BilRCiy4uPj2/9PSgoqHUo/A/Xf2huXlpayqBBgzj77LNb7zcajcf9OklJSeh0OkDtR2WxWKirq+PWW2/loYce4uKLL6ayspJHHnnkiLOcQgghhBC/x+G5jcViab0eFBSE3W4HoLy8HLPZfERhavbs2QwePPi4Xyc1NbX1927dulFVVUVUVBT/+9//eP7557n88ssZMWIE77zzDqGhob/3bQkhOjEpVgkhOj2z2YzT6fzdzzNw4EDmz5/PWWed9ZuKVD/YtWsXpaWlJCQk8N1332EwGFqbnz788MM8/PDDfPPNN8yaNUuKVUIIIYQ4bicj10lLS0On0xEVFcWECRNO6DmWLVvGmDFj8Pl8LF++vDWfmTZtGtOmTcPj8XDGGWewZMkSLrjggt8VrxCic5NilRCi08vKyiIvL4+pU6fSr1+/1qbnv9WVV17J559/TkZGBgMHDsRgMDB48GBuv/3249o+Pj6eKVOmEB8fz8aNG3n22WcBePrpp1m7di0+n49t27Zx3nnnnVB8QgghhOiaRo4cyZw5c+jXrx/33nvvCT2HXq/nf//7HzNnzqR///5ERkYC8Mwzz5CYmHhcz7Fy5UrGjRtHdXU1oaGhzJ49m4MHD3LrrbcCUFtby/79+xk5cuQJxSiE6DqkWCWE6PSsViu7du3iu+++w2KxkJmZycMPPwxAnz59jkjqxo4dS1ZWVuv1c889F7fbDahJ3KJFi9i9eze5ubm43e4jhtUfS3JyMgsXLmTTpk2kp6eTnp4OwIgRI0hOTsZgMJCRkUGfPn1OxtsWQgghRBfx4IMPMm3aNMrKyoiNjeWBBx6gd+/eADz66KOtK/YBvPzyy619o3Q6HQsWLGi975xzzmHcuHFs3bqVmpoaQM2jjteTTz6JXq/HZrMxYsQIjEYjERERzJ49G41GQ3h4OCNGjCAoKOhkvG0hRCemURRF8XcQQgjR2SxevJirr76acePG8fbbb7Nq1SruueceVq1addzPkZeXx80338znn3/OwYMHj+itJYQQQgjhT7169aJ79+5cf/31nHnmmQwePJgXXnjhN/W4uu+++/joo48444wzePzxx9swWiFERyPFKiGEaANlZWWsXbsWnU7HueeeS2VlJbt372bcuHHH/Rz19fUsW7YMgBkzZsgyz0IIIYQIGF9++SU2m42BAwfSo0cPVqxYwcCBA3/Tan9r1qyhvLycrKws+vXr14bRCiE6GilWCSGEEEIIIYQQQoiAofV3AEIIIYQQQgghhBBC/ECKVUIIIYQQQgghhBAiYEixSgghhBBCCCGEEEIEDClWCSGEEEIIIYQQQoiAIcUqIYQQQgghhBBCCBEwpFglhBBCCCGEEEIIIQKGFKuEEEIIIYQQQgghRMCQYpUQQgghhBBCCCGECBhSrBJCCCGEEEIIIYQQAeP/ATyJP+jaaKmpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1200x440 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            xnn   original MACE      |diff|\n",
      "------------------------------------------------------------\n",
      "(a) same PES              1.4293          1.4293     4.9e-08\n",
      "(b) independent           1.4463          1.4293     1.7e-02\n",
      "experiment                  1.41\n"
     ]
    }
   ],
   "source": [
    "t_ps = np.arange(N_EQUIL + N_PROD) * (DT / u.fs) / 1000.0\n",
    "xc = N_EQUIL * (DT / u.fs) / 1000.0\n",
    "fig, ax = plt.subplots(1, 2, figsize=(12, 4.4), sharey=True)\n",
    "ax[0].plot(t_ps, rho_xa, label=\"xnn\", lw=1)\n",
    "ax[0].plot(t_ps, rho_ma, label=\"original MACE\", lw=1, ls=\"--\")\n",
    "ax[0].set_title(f\"(a) same potential  |Δρ|={abs(da_x-da_m):.1e} g/cm³\")\n",
    "ax[1].plot(t_ps, rho_xb, label=\"xnn (independent)\", lw=1)\n",
    "ax[1].plot(t_ps, rho_mb, label=\"original MACE (independent)\", lw=1, ls=\"--\")\n",
    "ax[1].set_title(f\"(b) independently trained  |Δρ|={abs(db_x-db_m):.1e} g/cm³\")\n",
    "for a in ax:\n",
    "    a.axvline(xc, color=\"gray\", ls=\":\", lw=1); a.axhline(RHO_EXP, color=\"k\", ls=\"-.\", lw=1, label=f\"exp ≈ {RHO_EXP}\")\n",
    "    a.set_xlabel(\"time [ps]\"); a.legend(fontsize=8)\n",
    "ax[0].set_ylabel(\"density [g/cm³]\")\n",
    "plt.tight_layout(); plt.savefig(\"argon_density_md.png\", dpi=120); plt.show()\n",
    "\n",
    "print(f\"{'':<20}{'xnn':>12}{'original MACE':>16}{'|diff|':>12}\")\n",
    "print(\"-\" * 60)\n",
    "print(f\"{'(a) same PES':<20}{da_x:>12.4f}{da_m:>16.4f}{abs(da_x-da_m):>12.1e}\")\n",
    "print(f\"{'(b) independent':<20}{db_x:>12.4f}{db_m:>16.4f}{abs(db_x-db_m):>12.1e}\")\n",
    "print(f\"{'experiment':<20}{RHO_EXP:>12.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d53525d3",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "Computing the Argon density from **ASE NPT MD** two ways:\n",
    "\n",
    "* **Track (a): same potential.** Copying the trained original-MACE weights into\n",
    "  `xnn` makes the two calculators return identical energy/forces/stress (machine\n",
    "  precision), and the NPT densities are **identical** ($|\\Delta\\rho|\\sim10^{-8}$\n",
    "  g/cm³). This isolates and confirms the `xnn` inference/MD path reproduces the\n",
    "  original MACE exactly.\n",
    "* **Track (b): independently trained.** Training `xnn` from scratch (no copying)\n",
    "  gives an independent potential; its density agrees with the independently trained\n",
    "  original MACE to **within the thermal fluctuations**, and both sit near the\n",
    "  experimental liquid-Ar density (~1.41 g/cm³). This is the realistic\n",
    "  \"two practitioners, two fits\" agreement.\n",
    "\n",
    "Together: `xnn` is not only bit-for-bit equivalent to the original MACE for a fixed\n",
    "PES (a), but as a modelling tool it produces the same physical property when trained\n",
    "independently (b).\n",
    "\n",
    "*Notes.* Berendsen barostat for simplicity (use `ase.md.npt.NPT` for rigorous\n",
    "ensembles); longer runs tighten the estimate; stress is the autograd virial of the\n",
    "energy+force-trained model."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (xnn .venv)",
   "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
}