{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# `torch-sim`\n", "\n", "[torch-sim](https://radical-ai.github.io/torch-sim/index.html) is a tool for running (batched) molecular dynamics and structure optimisations in a GPU-accelerated manner.\n", "\n", "You can use **any** model from `graph-pes` with `torch-sim`! Doing this is extremely straightforward, as this notebook hopes to demonstrate. \n", "One caveat is that `torch-sim` expects all models to return energies and stresses in units of `eV` and forces in units of `eV/Å`. If you have a model that was trained in a different unit system, you will need to wrap it in a [graph_pes.models.UnitConverter](https://jla-gardner.github.io/graph-pes/models/root.html#unit-conversion) object so that it can be used in `torch-sim`.\n", "\n", "To get started, you must first ensure you have the `torch-sim-atomistic` package installed:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this notebook ran successfully with this torch-sim version. \n", "# future versions may require changes as the package evolves.\n", "!pip install graph-pes torch-sim-atomistic==0.2.0" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Using device: cpu" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from pathlib import Path\n", "\n", "import ase\n", "import ase.io\n", "import load_atoms\n", "import torch\n", "import torch_sim\n", "from ase.build import bulk\n", "\n", "from graph_pes import models\n", "\n", "%config InlineBackend.figure_format = 'retina'\n", "torch.manual_seed(42) # for reproducibility\n", "\n", "DEVICE = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n", "print(f\"Using device: {DEVICE}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## NVT MD\n", "\n", "`torch-sim` makes it easy to run NVT molecular dynamics using a Langevin thermostat. \n", "\n", "To keep things simple, we use a simple `LennardJones` model to run NVT for a single structure.\n", "It's worth repeating that this proceduce will be identical for any other model in `graph-pes`, be that a model you've\n", "[trained from scratch yourself](https://jla-gardner.github.io/graph-pes/quickstart/quickstart.html), a [pre-trained foundation model](https://jla-gardner.github.io/graph-pes/interfaces/mace.html) or a foundation model you have [fine-tuned on your own data](https://jla-gardner.github.io/graph-pes/quickstart/fine-tuning.html)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", "

\n", " \n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "

\n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# make a model: this could be any model from `graph-pes`\n", "model = models.LennardJones(sigma=1.5, epsilon=1.0)\n", "ts_model = model.torch_sim_model(device=DEVICE)\n", "\n", "# create a bulk Cu structure using ASE\n", "structure = bulk(\"Cu\", \"fcc\", a=3.58, cubic=True).repeat(4)\n", "\n", "# NVT at 300 K with a timestep of 0.002 fs\n", "torch_sim.integrate(\n", " system=structure,\n", " model=ts_model,\n", " n_steps=500,\n", " timestep=0.002,\n", " temperature=300,\n", " integrator=torch_sim.nvt_langevin,\n", " trajectory_reporter=torch_sim.TrajectoryReporter(\n", " \"Cu-traj.h5md\",\n", " state_frequency=10, # save every 10 steps\n", " state_kwargs=dict(save_velocities=True, save_forces=True),\n", " ),\n", ")\n", "\n", "# post-process the trajectory to save as an ASE-compatible file\n", "trajectory = torch_sim.TorchSimTrajectory(\"Cu-traj.h5md\")\n", "trajectory.write_ase_trajectory(\"Cu-traj.traj\")\n", "load_atoms.view(ase.io.read(\"Cu-traj.traj\", index=\"-1\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Of course, a Lennard-Jones model is not particularly realistic for this case, but we can see that MD has generated a collection of [Lennard-Jones clusters](https://pubs.acs.org/doi/10.1021/jp970984n) as expected.\n", "\n", "`torch-sim` saves structure snapshots in a custom format. To extract [ase.Atoms](https://wiki.fysik.dtu.dk/ase/ase/atoms.html) objects with the force and velocities properties, we can use the following helper functions:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(256, 3)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def get_atoms_with_labels(trajectory: torch_sim.TorchSimTrajectory, index: int):\n", " atoms = trajectory.get_atoms(index)\n", " for property in [\"forces\", \"velocities\"]:\n", " if property not in trajectory.array_registry:\n", " continue\n", " atoms.arrays[property] = trajectory.get_array(\n", " property, index, index + 1\n", " )[0]\n", " return atoms\n", "\n", "\n", "def trajectory_to_ase_structures(file_name: str | Path) -> list[ase.Atoms]:\n", " traj = torch_sim.TorchSimTrajectory(file_name)\n", " atoms = []\n", " i = 0\n", " while True:\n", " try:\n", " atoms.append(get_atoms_with_labels(traj, i))\n", " except IndexError:\n", " break\n", " i += 1\n", " return atoms\n", "\n", "\n", "structures = trajectory_to_ase_structures(\"Cu-traj.h5md\")\n", "structures[0].arrays[\"forces\"].shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## NPT MD\n", "\n", "`torch-sim` also supports NPT molecular dynamics using a Langevin barostat and thermostat. Below, we do just that. \n", "\n", "Note that `graph-pes` models return stresses in the energy and length-scale units that they were originally trained on: typically this is eV and Å, but not always. \n", "Use a [graph_pes.models.UnitConverter](https://jla-gardner.github.io/graph-pes/models/root.html#unit-conversion) to convert the energies to `eV` and lengths to Å if necessary.\n", "\n", "`torch-sim` expects the external pressure provided to be in units of `eV/Å^3`." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "623189fa7ae94765a126dd794c13ceb2", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "

\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Dataset:\n", " structures: '100'\n", " atoms: 3,552\n", " species:\n", " C: 100.00%\n", " properties:\n", " per atom: (forces)\n", " per structure: (config_type, detailed_ct, energy, split)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "amorphous_carbon_structures = load_atoms.load_dataset(\"C-GAP-17\").filter_by(\n", " lambda structure: len(structure) < 50,\n", " config_type=\"bulk_amo\",\n", ")[:100]\n", "amorphous_carbon_structures" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here we load in a model trained on part of the [C-GAP-17](https://jla-gardner.github.io/load-atoms/datasets/C-GAP-17.html) dataset:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsYAAAJcCAYAAAD6lWKQAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAewgAAHsIBbtB1PgAAdk5JREFUeJzt3Qd0lNXWh/GdDgkQEnqvEkA6iCBFEFEEQQRFQZQickUU26dy7dj12hXx2lGKiKiAIgIqUgURLEiXFmrogVDS5lv75M44mcykTp/nt1ZMMvNmckIc+Odkn73DLBaLRQAAAIAQF+7rBQAAAAD+gGAMAAAAEIwBAACAHARjAAAAgGAMAAAA5CAYAwAAAARjAAAAIAfBGAAAACAYAwAAADkIxgAAAADBGAAAAMhBMAYAAAAIxgAAAEAOgjEAAABAMAYAAAByEIwBAAAAgjEAAACQg2AMAAAAiEikrxcAz8nMzJQDBw6Yt6tWrSqRkXy7AQAAXGHHuAimTp0qN910k7Rs2VIqV64sUVFREh8fL+3bt5dnn31WTp065fJjFy1aJL1795aKFStK6dKlpXHjxvLQQw/l+zElpaG4Vq1a5sUakAEAAOBcmMVisbi4Dw46d+4sK1askCZNmpiwmZiYKAcPHpSVK1fKmTNnpGHDhvLTTz9J9erVc33cK6+8Ivfcc4+EhYVJly5dpEqVKrJ06VITVpOSkmTZsmUmMLvbnj17zDpVcnKy1KxZ0+2fAwAABJcffvhBrrzySpNt1Jw5c6Rv374SCgjGRbBq1So577zzTCC2d+TIEenfv78JuNdff71Mnz7ddt+6deukbdu2Eh4eLnPnzpUrrrjC3H769Gnp16+ffP/99zJw4ED5/PPP3b5egjEAAChJKO7bt6/JKNHR0RIKCMZuojvAXbt2NaFZg7LVoEGDZObMmTJq1Ch59913c33Mrl27pH79+pKdnS0bN2405RXuRDAGAACFtXPnTmnatGnIhmJFjbGbWA+2xcTE2G5LT0+Xb775xrw9ZMiQPB9Tp04d6dSpk3n7yy+/9NpaAQAAHNWtW1fGjx8fsqFYEYzd4OTJk/L444+bt7U8wmrLli2mZEK1a9fO6cdab9eSCwAAAF969NFH5dNPPw3JUKzo31UMCxYskGnTppkSCOvhOw3HvXr1kueff9523Y4dO8zr8uXLS9myZZ0+lrXUwXptUUsl8rN///4iPyYAAAgdml8cM8p1110noYpgXAwbNmyQyZMn57pNSyVefvll077N/n82FRcX5/KxypQpY16npqYWeR3WUA0AAFCcg3bXXnut2SHu2bOnr5fjFyilKIa77rpL9Myi1hBv27ZNXnrpJfn2229NwfqSJUt8vTwAAIBCdZ84evSoKQNds2aNr5fkF9gxLgEd8NGgQQPTo1gP0XXs2FGGDh0qmzdvNkM8rL+aSEtLc/kY1gEf5cqVK/Ln104TBZVS6PARAAAAVy3ZdLe4RYsWvl6WXyAYu8mFF15odoz/+usv81OXDvLQ053q+PHjTmt47MOt9dqioP0aAAAoilDvU1wQSincyFpLnJKSYl7rVLvY2FjztqtfUVhvb9OmjdfWCQAAQg+huGAEYzc5fPiw/P777+btRo0amdf6P1qfPn3M29rFwpEO+NAR0+rqq6/26noBAEDoIBQXDsG4CJ0opk6dKmfPns1zn/Yr1lOd586dkw4dOkjz5s1t92mj7LCwMPnwww9l/vz5ttu1v/HNN98sWVlZZiS0u6feAQAAKEJx4VFjXEhaHqEH6/71r39J69atTX2vdqXYvXu3rF271vQ0btKkicyYMSPXx2mJhHat0AN6vXv3losvvlgqV65sRkjr4Tgtt3j77bd99nUBAIDgpp0nNLMoQnH+wizadwwFOnTokLz77rsm0G7atMm8n5GRIYmJiWaHeMCAATJixIhcI6HtLVq0yATk1atXmy4VtWvXlmuuuUb+/e9/uxz+UVI6AMTa61gP+XFYDwCA0KRhePr06eaFUOwawTiIEYwBAICVRj4t74Rr1BgDAAAEWU2xszJNQnHBqDEGAAAIwoN2esB/7Nixvl5SQGHHGAAAIAi7TyxYsMCUT6DwCMYAAABB2JJt5syZlE8UEcEYAAAggNGn2H0IxgAAAAGKUOxeBGMAAIAARCh2P4IxAABAgDl79qzcdNNNhGI3IxgDAAAEmFKlSsns2bOlfPnyhGI3oo8xAABAAGrbtq2sXLlS6tevTyh2E3aMAQAAAsDGjRslOzs7122NGzcmFLsRwRgAACAADtrpDvHtt9+eJxyr1DMZMmHWehk6caV5re+j6AjGAAAAAdJ9YtKkSfLee+/lueaVeZtl6eZDcuDEWfNa30fREYwBAAACqCXb8OHD81y39cDJfN9H4RCMAQAAArxP8XlVy+b7PgqHrhQAAAABPrzj7t5Jtp1iDcX3/O99FA3BGAAAIMAn2pUrHSWPDWzmxVUGJ0opAAAA/MTSpUsZ8+xDBGMAAAA/0aBBA6lZs6Z5m1DsfQRjAAAAP1G9enVZvHix3HHHHYRiHwizWCwWX3xieN6ePXukVq1a5u3k5GTbT6AAAADIix1jAAAAHx6069evn5w+fdrXSwHBGAAAwLfdJ+bOnUs49hMEYwAAAB+3ZIuNjZXISLro+hrBGAAAwM/7FMM7CMYAAABeQij2bwRjAAAALyAU+z+CMQAAgIcRigMDwRgAAMDDXn31VUJxAOD4IwAAgId9+umnpiWbdp8gFPsvgjEAAICHaSCeM2eOaclGKPZflFIAAAC42dKlS2Xfvn15wjGh2L8RjAEAANx80O7yyy+X7t275wnH8G8EYwAAAA90n9iyZYs8//zzvl4SioBgDAAA4IGWbHrY7j//+Y+vl4UiIBgDAAB4IBTPnDmTmuIAQzAGAAAoAUJx8CAYAwAAFBOhOLgQjAEAAIph48aNhOIgQzAGAAAohqSkJBk+fLh5m1AcHJh8BwAAUAzh4eHy5ptvSqtWrUxAJhQHPnaMAQAACuns2bN5wvHo0aMJxUGCYAwAAFDIg3YNGzaUtWvX+nop8BCCMQAAQCG7T+zdu1cuvfRS2bRpk6+XBA8gGAMAABShJVuXLl2kfv36vl4WPIBgDAAA4AJ9ikMLwRgAAMAJQnHoIRgDAAA4IBSHJoIxAACAHUJx6CIYAwAA2NmyZQuhOEQx+Q4AAMDOrbfeKllZWbJw4UL57LPPCMUhJMxisVh8vQh4xp49e6RWrVrm7eTkZKlZs6avlwQAQMDQiBQWFubrZcCLKKUAAAAS6jXFs2bNynM7oTj0EIwBAICE+kG76667zmk4RmghGAMAAAn17hNaUzx9+nRTPoHQRTAGAAAhx1lLtmnTplE+EeIIxgAAIKTQpxiuEIwBAEDIIBQjPwRjAAAQEgjFKAjBGAAABL2TJ0/KtddeSyhGvgjGAAAg6JUtW1Y+/fRTKVWqFKEYLjESGgAAhISePXvK0qVLpUWLFoRiOMWOMQAACEo7d+7Mc1u7du0IxXCJYAwAAILyoF3Tpk3lySef9PVSEEAIxgAAIGi7Tzz66KPy2Wef+XpJCBAEYwAAENQt2fr37+/rZSFAEIwBAEBQoE8xSopgDAAAAh6hGO5AMAYAAAGNUAx3IRgDAICAtXjxYkIx3IZgDAAAAla1atUkPj7evE0oRkkRjAEAQMBKSkoyu8ajRo0iFKPEwiwWi6XkDwN/tGfPHqlVq5Z5Ozk5WWrWrOnrJQEAAPgtdowBAEBAHbS7+eabJTMz09dLQRCK9PUCAAAAitp94vTp0/LJJ59IZCRRBu7DjjEAAAi4lmwajLOzs329LAQZgjEAAPBr9CmGtxCMAQCA3yIUw5sIxgAAwC8RiuFtBGMAAOB3CMXwBYIxAADwKzpi4bHHHiMUw+vocQIAAPxKWFiYzJ49Wy699FIzqIpQDG8hGAMAAL+TmJgo33//vcTFxRGK4TWUUgAAAJ9buXKlHD9+PNdtCQkJhGJ4FcEYAAD4/KBdjx495PLLL88TjgFvIhgDAAC/6D6xevVq6XzdnTJh1npJPZPh66UhBFFjDAAAPELD7SvzNsvWAyflvKpl5e7eSVKudJTLlmwVG18klTvdIEs3HzLvPzawmc/WjtBEMAYAAB6hodgacg+cOJsr7DqG4prNO8t51zwi4ZE5wVnDNBDUwTgrK0tWrVola9eulYMHD8qxY8dMYX2VKlWkbdu20r59e4mIiPDmkgAAgIc4hlvr+86Gd7QY/Lis3P5PfbHuMANBGYyXLVsmEydOlG+++UbS0tJyNfDWXoVWZcqUkT59+sjYsWOlU6dO3lgaAADwEA231p1i6/uuJtqdzQrLVXZxT+8kH64coSrMounUQ5YuXSp33323rFu3zoTg8PBwadasmZx//vlSoUIFKVeunJw4cUKOHDki69evlw0bNkh2drYJy23atJGXX35ZunTp4qnlBb09e/aYxugqOTlZatas6eslAQBCuMb40hpnpOclFzPRDqEXjK+//nrzP3tkZKT07t1bhg8fLpdccomULev6VyOpqammmfdHH30k8+fPl8zMTBk0aJBMnz7dE0sMegRjAIA/ycjIMPngiy++KHYoLuhAH+CXwTgmJkZGjx4tDz/8sKkhLiqtQX7iiSfk/fffl7Nn//k1DAqPYAwA8Mdw/Nprr8m4ceOKtVOsrdysB/pUl6RKdK+A/wfjnTt3St26dUv8ODt27JB69eq5ZU2hhmAMAPA1/e2v/vbYXYZOXJmrbrlqfCmZMraj2x4foc1jAz7cEYoVoRgAgMCkB+30XNG2bdvc9piO3SroXoGgnnynu5wTJkzw9TIAAEAJWLtPbNmyRbp3725+c+kOWlOs5RO6U6yv6V6BgOlKUZT+xnPmzJF3331XFi5caDpT6G0oGUopAAC+4KolG90n4O98OvlOf7Xy3nvvyeTJkyUlJcW0dFMdOnTw5bIAAEAxEYoRyLwejNPT0+Xzzz83gfinn34yt2kg1prkoUOHyk033SQNGzb09rIAAEAx2LdPizi0QWa+cIdXQjFt2xAQwXjTpk1mrPN5552X6/a//vrLlEpMmTLFjILWMKxPlKioKDl9+rRs377d3UsBAAAepuFU26cd/Xut/P7JQ5Kdcc4rO8XWz6usXSpo2wa/OXx34MABufDCC83p08aNG8tVV10lR48elQ8++EA6duwoLVq0kDfeeMPc1rx5c3n11Vdl37590qpVK3ctAQAAeJnu2Ho7FFs/b37vAz7dMb799ttl48aN5vDcokWL5LnnnjODPfQgne4Oly9fXgYPHiwjR46Utm3buuvTAgAAL3BVuqBvr9i13haKk9p180pNsX5e+37GtG2DXwVjHeXcq1cvM/ZZD89pMNZQXKNGDXn++edlwIABZhoeAAAIPK5KFzQgWyz3yKwokYzDO2X+11/Js3O3eLz2Vx9XWT8PbdvgV6UUFStWlF27dtmm1SndKd67d68Z7axlFDrmOZBt3rzZfB3Dhw835SA6yScsLEyeeuoplx/z+OOPm2vye9G6bAAA/Jmr0gUNvY9f01z+mPeBrF+5UCYvSzYBWsOzvtZA7Qn6eTWY69Q7fV2Wg3fwpx3j//u//5MxY8aYHWM9SKehUcsrvvrqKxMoH3jgAXnwwQeld+/eppxCW7mEh/vdfJF8TZo0ycx3L46WLVu6rKeOj48v4coAABCvlS5oTXHN2oki8s8oZt3o0X/7qf1FIHNbMP7Xv/4lFSpUkBkzZkjlypVl7Nix0qVLF3n55ZdlwYIFpiOFDvHQl7lz55r6Y23PpofxAkWzZs3MDwCtW7eWNm3ayDPPPCOffPJJoT62f//+ZvcYAIBAZC1dWPrTj/LnlIdlc2S4DL+4nnTr1i3XddT+IpC5tV3bNddcY14cXXbZZebl0KFD8tFHH5kexlu3bpWXXnrJdo0G50GDBvn17umoUaNyvR9oO94AAJSkdKFLQoo8/94Dkpl+VjLTRd555508wZjaXwQyrya7SpUqyX333WdKKxYvXmy6VOiBPK1FvvXWW6VatWrmtm+//dabywIAAMWYaKebXY6o/UUg89lI6K5du5qX48ePy8cff2x2jHUIiJZiaJuXzMxMCSZr166V8ePHm9IR3RXXcoy+fftK2bLF/xXTnj178r1///79xX5sAACsGPOMUOGzYGyl/Y3HjRtnXlatWmV+LfPZZ59JsNG6an2xpwH59ddfN2Owi6NWrVpuWh0AAM57Fa9ZuZRQjJDhsVKKF1980da2rbB0ct77778fVDudDRo0MIf01q1bZ3aL9WXZsmXmL5kTJ07IsGHDZOrUqb5eJgAAtl7F1lZr/e6bJJf16k0oRsgIs2iBrwfowTRt3aJtygYOHGgGfDRp0kSCifYznjx5sjz55JPy8MMPF/njdZdc+yJr7bWWRRT1L5rClFK0b9/evJ2cnCw1a9Ys8hoBAMG3I7x5f6pkZVskIixMkqqXsw3hGDpxpa2jRHraCVnx0g2Sde60baLdH8u/IxQjqHlsx1h3fq+44grZsGGDPPLII6bVWdOmTeXRRx81u6fIGf4RERFhunVoGUlRadDN70UPMwIA4LgjnJJ6To6cSpeUk+dyDeGwb60WHRcvja+6WyQsXCo2vkhaD51AKEbQ81iN8YgRI8zLyZMnTW3tF198IfPnzzdT4p5++mmpW7eubSdZR0iHosTERNPzWXd2C9r9BQCguPXCdSvFmdt++dv57ADrEA77VmsZWdkiLXtIdNlEKV+7mTSuqQM9gODm8cN32nVhyJAh5kVrlLQV26xZs2TevHmmDll7GevOpgZkfbn44otNCUYoyMrKMnXGqiTdKQAAcLU7rOwHbjhj3Sk+k3rUtFjLFazjO9KPGCHDq32MS5cubcKvHjZLSUmRr7/+2tTpnjt3Tt58803p0aOHVK1aVUaPHm12l4OtZZsjnQJ4+vRp84NAu3btfL0cAEAQKWgUc4Uy0VK5XIx0SapkQq+2ZNMD4zqES9GPGKHIZ6PboqKipHfv3qYW+eDBg/L999+bIR86Z12flH369JFnn31WAtnu3btlypQpcvZs3p/Uv/rqK9skvRtuuMH8QAAAgLvkN4o5JjJcZozrJNNuv8iE3l/+15ItLS1NbrnlFvnuu++8ulYg6LtSlMTy5cvlyy+/NF0sbr75ZvGnIR233Xab7f2///5bDh8+bA661ahRw3a7rl3LQ3777TczyKNMmTLmtV6j5SR6IFFHYqvu3bubnWO9xt20btna65iuFAAQujXG5zKz5Fhahu2+Dg0ryFODWhQ4vMNZX2PdSQaClceCsT7JdEdUp7tp54VgoGOsNcgWRPs36+HCI0eOyH/+8x/55ZdfZNu2beb99PR0qVixorRt29bUXV933XWmtZ0nEIwBAMox4GrphJZGFDTRbsKs9bY6ZaVlF9YaZCAYebyPsfbovfHGG2XkyJFB18fY3xGMASA0FWantzBjnu37Gquq8aVMzTEQrDxWY3zHHXeYdmR6yO7ll182fYw7deokH3zwgalhAgAA3plgd8u7q01YLkoodlan7KpuWR9bd5c1SOtr+88FBBKP1hhnZGTI7NmzTRheuHChaU+mu8ixsbEyaNAgs4usYRmewY4xAIQmx51e+zKIwobi/EowHG/Xnsc/bzuS53MBgSbS050nrrnmGvOyb98++eijj8wIZT149uGHH5r3GzVqZALyTTfdJFWqVPHkcgAACAkaVh2DsbV9W1xcnPn3WYNxfqHYvmVbQT2Sw8OK1ioOkFBv11a9enV58MEHZfPmzbJkyRIZNmyY2TnW98ePH292Nq+++mrT2zg7O9tbywIAIOhoTbH2KXZWBtGkRRsZNH6S1L3gcmkx+HE5m1X0oVqOwTfbUvhWcYA/82m7Nq01njFjhim1WLFihW3infb03bt3r6+WFTQopQCA0OWqDMKx00RMVLi0r1+hSK3YHB/D+jgJsdG5PhcQaHw24MP66xwto1i2bJksWLDAtDHTnH7gwAFfLgsAgIBnLYMY2eSMnF79iZQpFel0t/dcRrYJuRqiS7IjreGaKXkIdB6tMS6I9vTVYRi6Y6yHAawlFBqQAQBAydgftDt37py88sorUrdSXJ7646LWBWvofveW9nl2pIFAF+mrCXIahqdPny7Hjx83u8Ta9/iyyy4zk+6uuuoqXywLAICg4dh9QodPZWZmury+qHXBrg7mAYHMa8H46NGj8sknn5huFH/++ae5TQNxnTp1ZMSIEebFWg8LAAByFGcss6uWbGcyRdbtOpbr2oiwMLmoUUV2fAFPB2MNvt9++60Jw3PnzjV9jfW2mJgY6d+/v9kd7tGjh+3QHQAAyL81mspvpza/PsXPzlpvaortaShm5xfwcDDW1mwff/yx7N+/34Rh1bx5cxOGhw4daqbiAQCA/DnW/uZXC+wYinv3udK0ZBv57q9mt3nzvtRc12v/4Vsuqe+hlQOBx2Pt2rRmWJUrV06uv/56GTVqlLRr184Tnwou0K4NAAK/jELHOR85lW67rUPDChIVEZ6ntMLZTrGG4pXbj9s+VjtJ2D+WYkod4IUd4y5dupjd4WuvvVZKly7tqU8DAEBQl1HYB1lrizTH0oqH+zeRu+66K0/5hO4U2zt5JiPP52BKHeCFYPzTTz956qEBAAgJjqUPelBu56G0XLet/vuIDHt7tVw69mU5+fy/pEXzZraaYsfR0OlZeX9JzJQ6wIft2lJTU2XKlClm0t2hQ4fM4bv777/f3LdlyxbZuXOndO3aVUqVKuXtpQEA4FeyHKodj59Olzb1EnOF3XOZ2eb9AxIhAx98T565qZMJxUrLLKy7wsdOp+c5eBcTGZ6nG0VxumAAwcKrwVin2w0ZMkSOHTtmDuRpN4oaNWrY7t+8ebPpVqH9jQcNGuTNpQEA4Hci9HScHeuOr9YFL1yyQqIS60hE1D8T6A6cK2ULxY69hp2NcW7foEKeKXVF7YIBBBOvjYTeuHGjXH311XLixAkZM2aMzJgxw9atwuryyy+X2NhYmT17treWBQCA30qqVi7PbVpK0SUhRX6eNE7+nPqoZGWk5ymL0F1fDcJDJ640r/V93fnVg3u6SxwTFS4dG1Zw2ru4KF0wgGDjtR3jZ555Rs6ePWvqngYMGGBuu+6663Jdoz/ltmrVSn7//XdvLQsAAL+lYXbD3hO5DuDt2fiLXHnfeMnKOCdHtq6W5OUzpW63G8zBPGvQdbXr+9SgFgV+Tse6ZGqQEUq8Fox//PFHadmypS0Uu6ItxTZs2OCtZQEA4HOu6nr15d1b2ttath39e638/slDkp1xznxcUrtu0rrPTdKgWnnRoosxH6xx2q9Y39ed48LUDdvXJeu1TMRDKPFaMNaDdp07dy7wOp3jnpaW+8QtAACBHnjrVoqzlUI4htOC6npPncvME4rtJ9rZ1w/rx1vbutkf4its3bB9XTIQarwWjOPj42Xv3r0FXrd9+3apXLmyV9YEAIAnOQu89m9bA2h+db36GPs3rcmzU2wNxc4+Xg/t6QE9667v5v25d5CpGwZ8HIzbtGkjS5Yskd27d0vt2rWdXrN+/XpTX6yH9AAACHT5BdDV24+YHWXdoc2vrnfpTz/mCsWVm3aS5Yu+lrNZYfLs/8ojMrKy8xzas9/11R3llNR/OlJQNwz4uCuFjoTWw3eDBw+WAwcO5Ln/8OHD5hrtVKGvAQAIdPkFUO0prLvBSssqdIe3anwp81rres046GenyY9v32cLxRUbXySjH3lDKsTH2XajNVBr/bGWT9h/vD1njw/AhzvG11xzjRkPrb/6adCggXTq1Mncvnz5clMntXjxYjl16pTccMMNpm0bAACBzv4gm9YYr91xNNf0OVc7ynqFBt9t58pLuRpJcnznH2anWEPx/Vc5L7+IigiXKWM7On086oaBwgmzODYT9qCsrCx59NFH5dVXX7XNc7fSOqk77rhDnnvuOYmIiPDWkoLanj17pFatWubt5ORk0/EDAOBeRZkU5zhkQ/sKa6DVsc46wc7KWh+su8FZ6Wdk15IZckG/kdK0VgXb59HyiZ+3Hcn1MYRfIICCsZVOvtP2bXrQLjs724Q3HQ3NoTv3IhgDgOcDsbNQ6yqgOoZox3CrLNnZUi0h1txvH6K1VMK+n7E1VNu3VXOcYgfAj0dCWyUkJBTYzxgAgEDpOFHYA3e6k6w7ytZwfCztn6CrtCXbtu/elX9NeFvu7t3O9njOOktsP3hKkqr/MxnP67tcQBDySTAGACDQuQrA+R24Mwfq/jesw5F9n+K5L46RB65amm9niaL0Jgbg464Ua9eu9avHAQDAnRwDcExUeIEdH3Sn2DEU68fJgfXy55SHbd0nsstUlRvf/U36vPCTPPzZHyZQO3aW0F7F9uhNDPhxML7gggvk+uuvl02bNhXr4//66y8ZNGiQeRwAAPyNY1D99PaLzI5tfnW+zsLroS1r5Kf/3i+Z6Tm7vtWbdZaaff8tGZYIU7usNcgaqK1lGBrI9XGysnMXT9CbGPDjUopx48bJxIkTTXu2jh07yvDhw80Bu3r16rn8GD2Mt3DhQvnoo49k9erVpjvFnXfe6aklAgBQZI4H6N4a2S5PFwpXnSocB3k4jnnWPsVJ1z4i4ZFRTgO1Y12zHsjTA3jWw3cA/DQYv/LKKzJ69Gi577775Ntvv5WVK1ea2ytVqiRNmjSRChUqSLly5SQ1NVWOHDkiGzdulEOH/nmy9+7dW1544QVzLQAA/jzmWXeK7cOwdpuwlkzYXzPqkvq2j3UWipsPfixPKFb6ePr4znoXazDXzzvmgzUFtosD4MPDdxpqv/76a9m6dau8+eabMnv2bDMSOiUlxen1Oiq6f//+MnbsWDnvvPM8uTQAAIrFMZy62s11ds17P2wvcihWGrL18R13nHVoiP1hPg7hAQHQlUJD7muvvWZetFxi3bp1cvDgQTlx4oSUL1/e9C9u06ZNvmUWAAD4A8dwaq3t3bwvdzs1expgtavE8i05wfnQhmWFDsVW2i85PjbKlE9EhIWZVm2ZdjvTVhzCAwKoXVv9+vXNCwAAgTjJzn7Ms7W2Vz/m+OncATU6IkzKx0WbQ3KOo6Ab9bldsjPSJT3teKFCsdKDeCmp53INERk6MadM0R6H8IDio48xAAAF1BDbD+WwhmXrtVrbqzXA9sFXaShOqlbOaXlFWHi4NO5/j1iys/INxfGxkRITGSEnzmTIuYzsPLvCjrvXupvMITyg+AjGAAAUUEPsGJY37D0hp85l5gqrjjQUWx/r6PZ1ElW6nJSt1iBXONaX/LSolWB2hrUMwz5gW3eFne1eMxYaKD6CMQAABdQQO4ZlZ5Pr7EVFhJmQ+vK8zbJh7Qpz0C4iKkZaj3wpVzguiPXzOgvASks8OGgHBMCADwAAgmFwh4bQotbthoeFmZ3bNrH75I//dZ/IOJ0qySs+L9LjHEtLN7vFSgPwlLEdCxwiAqD4CMYAANhxnDD3/NyNciY9U2IiC/9Pph6U63b763L1Vf0ky677ROOr7nH5MXpYz9njaAmFlnIA8DxKKQAAKMQBvKIoap9iZT28ZwJ4mDg9bAfAs9gxBgDAQX5BNDzvxm6JQ7G9hLhoaV+/Qq7baMEGeAc7xgAAOPQv1vZrxVHSUKzs28E5HrYDEALBeMeOHfLHH39InTp1pFWrVr5eDgAgBDmOdNaaX8fexCo7703G2ROHihSKnT2+bkZrCHb8FC4+JYBALaWYM2eODBgwQFavXp3r9v/85z/SqFEjc1/btm1l5MiR3loSAAAuyyeKumlcKr6SNOh5c6F3ii9sWNEM5LCXWCbadJywhnStb9bXt7y72ky50w4VurMNIMCD8ccffyzz58+XJk2a2G7btGmTjB8/XiwWi7Rs2VJiY2Nl8uTJMnfuXG8tCwAAp3W82Zai79PW7nSNtLjhyQJDcViYmE4Xrj6Fs77J1pBMhwogCILxunXrTPgtW/afv3imTp1qXr/11luydu1a+eWXXyQiIkLeeecdby0LABDidAd2/PTfZJnD6ObCxOKMM6fy3FapaacCa4o1EK/ZcUxOnc29+xvxv5N9+R22o0MFEATB+PDhw1KjRo1cty1evFhKly4tw4cPN+83btxYOnfuLH/99Ze3lgUACHG6A6sh1VKMg3YrXhwshzYsL/bnznQoWNYx0o5DRhzLLehQAQTB4buzZ8+a3WCrrKwss0vcoUMHiY7+50lfvXp1+fnnn721LABAiHehWLHlcIm6T/z56QRpc/PLUr5O0UczR0WES/sGFfId9WzfLYMOFUCQBOPKlSvL1q1bbe9r+D1z5ox06tQp13V6W1xcnLeWBQAIUY5dKIrbkq3Cee2lXI3ChVUtlLDfI25dN8EWgF2xD8kAgqSU4qKLLpLff/9dPv30Uzlx4oQ888wzEhYWJpdeemmu6zZu3Gh2jQEAcDfdfdXODtrhYfXfR7zep7hC2Wjp0LCCmW4XExVugjJdJoAQDMYPPPCAREZGyg033CCJiYny7bffSps2baRr1662a5KTk02nigsuuMBbywIAhAgNoNr2zNoG7VxmtteHdzSsklMfrJ9bRz6v3HZEXpi7sYhfCYCAL6XQEDxv3jx5+umnJSUlRdq3by/PPvtsrms+++wziY+Plx49enhrWQCAAORYd6uH1bTkIL9rdJqdtj0rjuKG4vjYSDlxOtP2/ub9qXIsLfcO8bpdx4q1JgDuF2bRJsIISnv27JFatWrZduNr1qzp6yUBgFtoOYR9fbCWJughNvuAnOeaqHCzS+vpUJwQF2UO1WmHiY37Tsjhk/mHcV37N/dfXOR1AQjgUgoAANzFsZevliZoCLYvS8jT77eY20DZmRliyc4u9E6x7gjr7rR+/tTTBdcP6wE8ACEajFNTU81Aj6FDh8rll18uL7zwgu2+zZs3y4IFC0xrNwAAXHHVy/fnbUdsY5Mdr9EAGh2RM0CjKComXSgtbnhCKje7uNA1xRqOtY45Pct1GtcdbO1V/EDffybCAgiRGmOloXfIkCFy7NgxMwZau1LYD/3YsmWL9O/fX6ZPny6DBg3y5tIAAAFESyaUs3ZretuGvSfk5Rtbm/d157ZupTjTAcJxoEZRwrG+lFRURJhUKBNj60dc1qEuGkCI7BhrG7arr77atGobM2aMzJgxw4Rje7qDHBsbK7Nnz/bWsgAAAcja21dbnzmjh+zu+WSdCcU1E0vLrzuOmg4QhcnFWlO8a9lnnvlHNyxMpoztaNZOKAZCeMdY+xZricTMmTNlwIAB5rbrrrsu1zU6Aa9Vq1am3zEAAAW5v28T03li9fYjeQ7WWTtQaElDcQ7aWbKzpG7XwUVek5ZrtKmXKDsPpUlK6tlChXEAIbZj/OOPP0rLli1todgV7Zywf/9+by0LABAEO8fTb79IKpSJLtFjOXafOLHrL9uhu6LQnWAN7LozrJ0y7HHQDvBvXgvGhw4dkkaNGhV4XWZmpqSlpXllTQCA4AnIWlMcXvSzdfm2ZAsLL/o/k7pTrbvYSgOyHrCrGl+Kg3ZAAPBaKYUO7ti7d2+B123fvl0qV67slTUBAIJnsMcb87cUq2zBHRPtHK3Ycth0x9D16Y42gMDg1cl3S5Yskd27d0vt2rWdXrN+/XpTX6yH9AAAoSm/8Ku3WztRaO3wss2HJCoy3HScKOqIZ0+FYpVlsdjWSTAGAofXSilGjRplDt8NHjxYDhw4kOf+w4cPm2u0U4W+BgCEJmv41eCrr61lCc6GdugGcXpmtl+FYnt5howA8Gte2zG+5ppr5NprrzVdKRo0aCCdOnUyty9fvlz69esnixcvllOnTskNN9xg2rYBAEKTY5jcuC9Vrnt9uRw7lV7c4XVOp9lt/PJFj4bi/AaRAPBPXp18N23aNPn3v/9t3l60aJF5vXXrVvn6668lPT1d7r33Xvnoo4+8uSQAgJ9xDJNHTp4zB9p0T9hdwVgDcMsbn5aouAS3hGLtiKEH7LSvcseGFWyH7XSIB4DAEWZxnLLhBTr5Ttu36UG77OxsqVWrlvTo0YNDd262Z88e82erkpOTTSs8AAiUGuPN+1PleFp6vmOVS+rM0X0SU65SkUKx1jOXj42S0+lZthZs2m2CgR1A4PNJMIZ3EIwBBPKBuxfmbpSftx1x2+c4dWC7xFWuI2HhESV6nJjIcJl+x0W2A4EAgodXSykAACjsgbt1u4657fH1oN0vb4+VDbNeMBPtSkIP+t3y7moT5gEEF68F4+nTp0v9+vVl/vz5Lq/R+/Sazz//3FvLAgD44YE7HfFscdMsZfvuEwd+Wyh7Vs0u8WPaD/EAEDy8GoyPHz8ul1xyictrunfvbuqPp06d6q1lAQD88MDduYxsKVOq5KUKzlqy1bigr7gDrdiA4OO1YPzHH39IixYtJDra9Sz7mJgYadmypRnyAQAIDVqSkJGVtw/xqXOZJXpcd/cpjo7IPW+aVmxA8PFaMNahHjVq1CjwOr3G2QAQAEBwcnXITgd3+EMoDg8T04atZZ0EiYkKN4fv9H1asQHBx2sDPmJjY+XIkYJPF+s1+e0qAwCCqwPFup3uO2Tn7lCse8SJZaLNWrWu2CoqIpz2bEAQ8lowPv/8882Uu6NHj0piYqLTa/S+ZcuWSbNmzJUHgGDvQKG0C4XuyLrL8V3r3Vo+ocf/Dp/8JxBbUV8MBCevlVIMHDhQ0tLSZOjQoXL69Ok89585c0ZuvPFG81rHRwMAgpMnQ2VsxVoSm1jDo2OeFfXFQHDy2oAPDbzt2rWTTZs2SbVq1WTIkCHSuHFjc5/epl0r9u3bJ0lJSbJmzRpTeoGSYcAHAH80YdZ6246xJ6SnnZBdS6ZJg56j3BqK9fBdYpkYE4q1vphSCiD4eK2UonTp0qZP8YABA+TXX3+Vl156Kdf9ms9bt24tX375JaEYAEKsA0VJ6L8fYWH/1GNEx8XLeVeMKdZjRUWESYaLEdTl46JlytiOxV4nAP/ntWCsdPdy9erVMnfuXBOSd+3aZW6vXbu29OrVS/r165frLzcAQPDVF7tzzLMetNv103RpPuRxiSwVV+LH01CcEBslJ85kiON8kaRq5Ur8+AD8m9dKKT7++GPTp/i6667zxqcDpRQA/NDQiSvNgTt3d5+Ir32+tBr+vETGuP83jtqerX2DnPZslE8Awc1rh+9GjBghH330kbc+HQDADzkeWgtzU0u2qNh4CY/wTGhNiIuWxwY2IxQDIcBrpRQVKlRw2aYNABAavYs370uVhLgoOX0uy9QGp7uo5/XmRLuC0IECCB1eC8YXXnihGQsNAAi9QR564M5+QEZxeTMU65S79vWZcAeEEq+VUtx///2yceNG+e9//+utTwkA8PGoZ23LpjXFgRaKdVe7Ra14E+pfnrfZhHwAwc9rO8b6K7Nbb71VbrvtNpk1a5YZ+FG3bl3Txs2Zrl27emtpAAA37QzXrZTTGWLnoTRJSXXPITtvhGJT3pGeJecyclrJHUvLkF+254yqth4W1DpjAMHNa8G4W7duphWbBuRFixbJ999/7/JavS4zM9NbSwMAFDEAa93t3b2TpFzpKLMzbG3B5q6OE472rfnGY6G4Qploee+W9jLmgzUu188IaCA0eC0Y6w4wPYoBIHBpKLZOrLPfRV23M2dn1ZOaDhwvWen6OcPcXj5x6lym6BFADfuugjEH8IDQ4LVgvHjxYm99KgCABzjumtre98Kehwbh5oMft73tTlo+oaFfd8DV6r+PyLnM7Fw7yhzAA0KD1w7fAQACm+OuqXaa0IEdcdHu32M5tuN3OXPsQK7bNBB7qiWbhuHbPlhj3o536FccFRFOD2MgRHh1JDQAIHBZd1SdtV+LjggrVk/i/A7aRceVlzajXpHSCVXF3RLjouVo2j/r1x1iLaPQF90htkcZBRA6vL5jrC3btDtFUlKSlClTxrzo22PGjDH3AQD8kx6005riKWM7ml1Ue+46Q2LffeLs8YOya8mn4gmREWHSJamSVI0vZfoV24sI/+c+fU0ZBRA6vLpjrCOhNRRnZGSY7hRWW7duNS8ffvih6XM8bNgwby4LAFBEjgfV7Gty3dmSrVGfseIJSdXKmR1wrS1evf1InvtozQaEpjCLfUL1oF9//VU6dOggWVlZcuWVV8rNN98sDRo0MPdt375d3n//fZk7d65ERkbKihUrpF27dt5YVlDbs2eP1KpVy7ydnJwsNWvW9PWSAARZ6zbHg2qBMLwjPEykYtkYycq25CoHiYkMl/YNcibdUVMMhCav7Rj/5z//kezsbBOAR4wYkeu+Zs2aSb9+/cyO8siRI+Wll16S6dOne2tpAIBi9jN2B2+GYpVtEUlJzflc9hLiotkpBkKc13aMa9SoIVWrVjU7x/nRneL9+/fL3r17vbGsoMaOMYCShF8d3GF6FIeJtK6TIPf3bWLuG/XO6lwH1wIpFOuusKsdbq0nJhgDoc1rO8aHDx+W7t27F3hd48aN5c8///TKmgAAzumOsHWandK3B722TLTxhLu2U04f2evVUKzBV1mHlCjtQKEHCbVmmkN2ALwWjMuXLy+7d+8u8Dq9Jj4+3itrAoBQ5mrEs3JWJuGGUuJcSidWl1odrpZdSz/1eE2xtXbYmumtXzP1xAB8EowvuOACmT9/vvzwww9yySWXOL1G71u+fLn07t3bW8sCgJDlasSzFDAe2V20xVuDy2+RuCp1pUrz7m4Pxa4O01EuAcDnwfiOO+6QefPmSd++fWXs2LGmJVu9evVsXSn04N2kSZNs1wIAPMtxV3jzvlSZMGu9ub1upTiPfM6s9LMSEV0qVziu1voyt+wKt6mXIJHh4bLzUBq7wQD8+/CdeuSRR+Tpp5922Qhel6LXTJgwwVtLCmocvgOQHw3B9vW2jtPrIsPdWz6hB+3Wz3haWgx5XMrXbS7uxuE5AAE1+e7JJ580u8Z6CC8mJsYEYX2Jjo425RV6H6EYALxDa4rtxx87jnR2dyjWg3YZacfkt8nj5eSBv8Xd3NU+DkDo8urkO9WrVy/zooM+jhzJOfFcoUIFiYiI8PZSACCk6UE7x9HOnuDYki2hfmuJq1jb7Z9HyycAoCQ89jdiamqqnD3r+uCGBuHKlSubl1AJxTNnzpRu3bpJQkKCxMXFScuWLeWFF14wI7IBwBc8HSa91ac4JiqcdmsA/DcYa/jTQ3bOfPzxx2bscyi56667ZNCgQabrRvv27c2uubame+CBB0wZyZkzZ3y9RAAhWk6htblV40tJQpx7w6o3h3e0r1+Bg3YA/DcYW+uHnRk+fLi89957Eiq++uoree2116RMmTKyatUq+e6772TWrFmydetWad68uSxbtswcOgQAX5RT6IG1KWM7Bmwo1rZsGVnZpi8zAARUjXEoeuaZZ8zr8ePHS5s2bWy3V6xYUd566y3p0qWLvPnmmyYcM9wEgKeHeWgrtszsbPkz+YToxIvWdXPGPR9PywiIUJwYFy3R/wvDR06lmxHPOplv8JsrzM6x/aASAPDbrhShaO/evfLLL7+Yt4cMGZLn/s6dO5uWaufOnTNdOQDAk8M8dGiHhsg124/JuYxsW6i85tVltqlwJXUu9bBkZ6Z7bKdYQ7HucDseHNSvR79G/VoBoDgIxh62bt068zoxMdE20MRRu3btcl0LAO6mwzvyk+3GjvY6sKPpwAekUtPOHimfsB4YdHVwkLZtAIqLUgoP27Fjh3ldu7br1kTWIRzWa4sywCM/+/fvL9LjAQjO0gkNkBnZbmxKXMhwXLVVT5cDndzRfUJLJtTqv4+YnW8r2rYBKC6CsYedPJmzc6Ht2VzRQ3nWFndFYQ3UAOCqdEJp+YQeUPMUrSnW8gnH0c7uDsWO3SesBwcdfwigbRsAvwzGkydPNi+O9C9LV/dZ78/MzPTk0gAgqOUpJ3B/Rs190E5rii0Wqdbmco98Ht0pbl0nwWnotQZkAPDrYOyqXVsoKVs251d6aWlpLq85deqUeV2uXLkiPXZycnKBpRTaMxlA6NGdU90ptj+Y5m6O3ScObVwuVVtf5pGdYl2/HrajVzGAgAzGP/74o6ceOqDUrVu3wBBrvc96bWHVrFmzhKsDEKxs9bfbj3glFGv3iWbXPeKRUGzFoToAARuML774Yk89dEBp3bq1eX3kyBFzuM5ZZ4o1a9aY1/Y9jgGgqBxrbTUc3/bBmlw7x4E2vMMeh+oAeBrt2jxMd3UvuOAC8/a0adPy3K9T73THOCYmRnr37u2DFQIIxl7F1n6+7g6TvgrFFcpEc6gOgMcRjL3gwQcfNK+fe+45Wbt2re123UW+7bbbzNu33347U+8AlMjm/al53h91SX23nbvzVShW1BcD8AaCsRf0799fxo0bZw7ZdejQQa644gq55pprpGHDhvLnn39Kp06d5Mknn/T1MgEEuCyHKR2ZWRa555N1bplol5V+VtbPeNpjobhc6fwr+yijAOANBGMvee2112TGjBnSsWNHWbFihRn/rGUWuov8ww8/SOnSpX29RAABLsLh4Nupsxly5FS6ex47upS0uGGCRESX9shOcdMa8Xl6LYeHiVSNLyVdkipRRgHAK8Is9FQLWjoZzzoEROuY6WIBBLcJs9bbhnp4yqkD2yW2Yi23l09oAK5bKU5+3nbEdluHhhXkqUEt3Pp5ACA/7BgDQBB0o3j4sz/MaGTdZY0Mc888j9NH9ubpR1+man2P1BRrqcT9fZuY3WHrLvEDfZu4/fMAQH4YCQ0AARqGX5i7UdbtOiYZmdliX16c7caDdjU79JeGl4/2SH/iqHCR8PBwaV03Z6KdHq5jgh0AXyIYA0AA0lZs9mUH7mTffWL30hlSpnJdj4x6LhcbLTPGdXL74wJAcVFKAQAByFNT4Jy1ZKvS4hKPfK7jaRkeeVwA8Lsd4yeeeKLYH6u/snvkkUfcuh4ACCZakxvoE+3iY+lLDCBEulJo3ZgGXMeHL6hOTa/Xa7KysjyxrJBCVwoguGuMr3t9uWRkWQImFOtf//b/JJSPjZJSURG28dXlGOABIFh3jB977LE8t+3YsUM+/vhjKVWqlFx22WVSr149c/vOnTtlwYIFcvbsWRk2bJjUrVvXU8sCgKCgITIhLlpSUnOCbCDsFDtuwxw/raUUGbadbw7eAQiZYKw7lm3atDFT4CZNmiRVqlTJdX9KSorceuut8s0338iaNWs8tSwACBpJ1cpJSmrJ+hYf2/G7z8Y8e6NmGgD8siuF1gxHRUXJtGnTzI6xo8qVK5v76tevb66dPHmyt5YGAAFROqGdKDRA6iAMtf3gqRI/bqnyVSQ6rrycPX7QY6FYJ9qdy8y/iRwjnwGEVDDWUomuXbs6DcVWel+XLl1k4cKF3loWAPhl+HWsu9XbrVPt3HnornRCVWkz6hXZteRTadRnrNtDsQ4c0T7FrlrL6RjrixpVZOQzgNBq13b06FE5c+ZMgddpnfGxY8e8siYA8CfW8KvBV1/r+94oNdBw3PiquzxSPnFhgwq5JtpVKBOd634NxVpbrMM9ACBkdoxr164tP/74oxw8eDBPfbHVgQMHzDXWTgoAEEocw6/9++5qz6YH7fb+8o2cf80DEh6ZO6S6m4ZgDcX2E+0cd8XZKQYQkjvGQ4YMkVOnTkmPHj2clkosWrRIevbsKWlpaeZaAAg11tphZ+9rWYXuusboHOUSdp9I+fNH+XP6BMnOTBd3SYiLksplY0wYrlwuxqz1vVva59kJLve/kDxlbEd2igGE7o7x+PHj5bvvvpNVq1ZJr169pGLFira2bNqu7fDhw6aH8YUXXmiuBYBQt3l/qgyduNLsrI66pL5kZGXLuYz8D7EVtiWbiPaUz7+vfFGESZi8PeoCehEDCGheC8Z6sO6HH36QRx99VN555x05dOiQebEqU6aMjB492kzMy++AHgAE6wG7nYfScl17LO2fPr9/JB+TE6cz/bZP8dG0dHlh7kZ5alALtz0mAARtMFalS5eW//znP/Lkk0/K2rVrzWQ2VaNGDWnbti2BGEDQc9Zdwlp/m18dsT+HYqt1Ozk4DSCweTUYW2kAvuiii3zxqQHArw7Ybd6XKhNmrbf1J+7QsIIJmAX1/fW3UGy4rzIDAIL78J2jbdu2ycqVK2XLli2+WgIAeJ3jIIssi8XWok17/UZFhEt8bJTfh2Jrf2J7berkfh8AAo1Xg3FWVpY89dRTUrVqVUlKSpLOnTvLc889Z7t/6tSpZif5r7/+8uayAMBrrN0ltKevvo7QhGlHQ/KhVOsBueLRg8y7lkz36E5xp0aV5JGrz8/1tWhrNgAIZJHeDMVXXnmlmYAXGRkpTZo0kQ0bNuS6plOnTnLjjTfKF198Ieeff763lgYA+R6KcydruzIrLaNISf3nILKylPBzhIWFSfMhE+S3jx6QqNh4t4diHfGs/Yft+xMDQDDw2o7x22+/bdq1de/eXXbs2CHr16/Pc422b2vQoIEJzwDgL1PnvLGD7G6RMbHSavjzHqkpbt+gAv2HAQQlrwXjyZMnS2JiosycOVOqV6/u8jrdSd69e7e3lgUABU6d8yTrDrLjqOSiOr5rvaSnHc8Tjt0ZivUfDA3xTKsDEKy8VkqxadMmU1OckJD/4Yz4+HhJSUnx1rIAwGmrNMdDcp4s39BBHjq8o6QH7WITa0jrm1+U6Ljy4gmdkipROgEgqHm1xjgmJqbA6/bv31+o6wDA3SUNylpj7KldUftaZg3DR06lu637xKmD282hu/OuGCPuFBEWJhc1qshOMYCg57VgXKdOHfnjjz/yvSYjI8PUHp933nneWhYAOD0U540BH55oydag5yhxNw3F7BQDCAVeqzHu1auX7Ny504yDduWNN94wY6L79OnjrWUBgFe5q3bZG8M7dKeYmmIAocRrwfi+++4z9cO33Xab3HXXXbJixQpze1pamhkPPX78ePNSsWJFuf322721LADwKnfULntrot0FDRLN6zEfrDFt5bQMBACCWZhFO8F7yZIlS2TAgAFy9OhR02fTni6jfPnyMmfOHHNIDyW3Z88eqVWrlnk7OTlZatas6eslASFNg+UzX/0la3Yc8/tQHB0RJi1ql8+1Vt09pqQCQDDz6uS7rl27mql2999/vxngUbp0aXPQrmHDhjJu3Dj5888/CcUAgpbWF5ckFJ86sN0roVilZ1nkz+QTPmlhBwAhEYxVlSpVzBhoPYh36tQpOX36tGzevFleffVVqVGjhreXAwBes3lfaok+Pq5yHanc7GK3h2LdHV70YHcz2tmeYws5b7SwA4CQ6EqhZRRVq1aVRo0a5Xvd1q1bTcs23V0GgGCSVcLKtbDwCGk64D4pV6OR1Ligr9t2isvFRjvt5Zxtt1wdQMIhPADBzms7xt26dZPnn3++wOteeOEFMzYaAAK1jlgPqg2duDLPgbWI8NxnKwojOzMjTziu1XGAW8snmlQvl2s8te4cx0Tl/uchKiKcMdAAgp7XdoyVF8/5AYBP2Pcp1t1XLUfQUKn1uVn2W7CFPGi38csXpeWNT0uZKvU8st6YyHDbTrB9L2cN9fb9limjABAKvBqMC+PYsWNSqlTuOjcACBSOB9R+3nakxN0n1r5/r7T71xsSW8H95zDaN6jgdCfYW5MAASBkgvHu3btzva+H7Rxvs8rMzDQdKxYsWCANGjTw5LIAwGMc63Td0ZItvlZTKRVf2U0rFEmIi5KYyIh8A6+3JgECQMgE47p16+bqVzxr1izzUlC5xdChQz25LABwC60f1tIJ666q7rKOuqS+rN9zXI6lZfhtn+LG1crJk4NauO3xACBYeDQY165d2xaMdac4NjbWTLZzJjo62gygGDhwoIwZM8aTywIAj9QTW/lTKNa/gR0rm3ccSiv24wFAMPNoMN65c6ft7fDwcLn22mvlgw8+8OSnBACf1ROXZACGp3aKnR33q1sprkSPCQDBymuH7z788EMz4Q4AArlUQmtvXdUTaweK4tQYe2vMMwDAT4LxsGHDvPWpAMBjpRL2B9I0KG/Ye0KOnEo37+vrkzuOSmREmGRmFb4126n9f3s1FO+klAIAfDvgY+XKlTJy5EhZsWKFy2uWL19urlm9erW3lgUAxS6V0N1j7VFsLz3LUqRQrGp3vlYaXj7abaE4OjJMOjasYAZ1dGhYQRLjcibbWdGTGAB8vGP8zjvvyPTp081kO1d0XPTUqVMlMjJS2rdv762lAYBTjmURzgKl1uuWtD2bqtP1eqmdnS1h4SXfr2hTNzFX1wnHkhB6EgOAj4Ox7hS3atXKZVcKValSJWndurUsXbrUW8sCAJfyG3JhDZvrdh4r8uNqTbGOeq6YdGGu24sbirX5j/1gUcfB0/QkBgA/C8Z79+41obcgderUkfnz53tlTQBQ0GE7x0BpvWb19iNyLiO7yJ/DetDOkp0tLW54Ik84Lg77UKz+TjllRjq7OjQIAPBxMNZ2benpOQdU8pORkWGm4AGAtzgLu84O2+l1t7y72nbYrqgcu0/sXzu/WME4KiJMwsPC5Fym82CelW3J99AgAMDHh+90Cp4ewMsv9Op9eo3uGgOAt7tPOO4AWw/baSDWHdjBb6xwWyiu2LijnH/tg0V+nOiIMPlg9IXSvkGFXLdXKBNtDtt1SaokEeFhbuuvDAChxGvB+PLLL5eUlBR5/PHHXV4zYcIEc02vXr28tSwAcBkcrYftbMHZxQ5tcUJx88GPF7n7hIbeGeM6SbWE0qY8Qt+3dp6wPxhYv3IZp18HACB/YRaLY3WaZ+zbt0+aNWsmJ06ckH79+sktt9wijRs3Nvdt2rRJ3n33XZkzZ46ULVtW1q9fb8ZDo2T27NkjtWrVMm8nJyfzZwq4oLvB1tIDFRMZbnZk9bBd2dJRMuTNFZKSmhNqfRWKNfw+ZddpIr/167XaRs7+0KB+HQAAP6kxrl69usycOVMGDhwos2fPNiHYnuZzDcV6DQEOgK+7T9gHSa3Z9WUo1jKJB/o2KfSOtw7wmDK2Y7HWDAChzGvBWPXo0cPsBr/00kvy3Xffya5du8zttWvXNuUT99xzj22HEwC8xVk7M/sOFSfPZBT5MTPOnJI/pz1W4lCsjqWly6h3V5va4aRq5QocTU3pBAD4eSkFvI9SCqD4HMsTiuPQxuXy5/QJUuG8C4odip3R2mLHbhmOAzwonQAAP98xBoBAoEFTW7eVVKUmnaTNzS9LuRpJhQrF2oZNd4KtAzt0p9hZFYez0dS0YwOAkiMYA4Ada6/i4gzvOHvikJSKr5TrtvJ1Ch9YtTfxu7e0t5VJuNq1plQCAAIsGI8cOVLCwsLkmWeekSpVqpj3C0s/7v333/fU0gDAJS1JKE6vYutBuwY9b5bana4p1ufWdnD6+a27v9ZDgZv3pUqWxWKrMbYfTQ0ACIAaY510pwF348aN0qhRI/N+oRcVFiZZWVmeWFZIocYYKLqhE1fmOshWnO4TLW54Uio17eT0Wh298fGYDjJx4Vb5eVvecg3tS0xHCQAIsh3jDz/80LyuVq1arvcBwB/LJ16Yu1HW7TomGUUc4uGsJVuFRu1d1hDr1Dod0DGmZ0P5dfsRcazYoEwCAIIwGA8bNizf9wHAX2j5grPdW3f3Kdbd4rhSOX/tvvfD9jyhWCc5UyYBAL7D4TsAIc/VSGgNqq5mexRneEd6lkUGv7HCTNXbvD/V6TX0zwQA3yl84S8ABClX5QsaiqMjdJ/XfRPt9ICddpo4npbu9PPp7jUAIMh2jJ944olif6wevnvkkUfcuh4AcEW7P2RkZcvanUclPdOSZ5fX3tHt69wy5lkfNyEuSo6lZRRq9xoAEMDB+PHHHzcB17Hphd6WH72eYAzAG6wT47SsISvbImGmCjj/Yoao2HISERVjgnFJxjyrmMgI6ZJUPlevYg7fAUAQBuPHHnssz207duyQjz/+WEqVKiWXXXaZ1KtXz9y+c+dOWbBggZw9e9Yc0qtbt66nlgUANhqKizr2uWzVBtL65pckefnn0viqewoMxVqK4bjrbB+Crb2K7cc5AwCCrI+xI+2j26ZNG+nSpYtMmjTJDP2wl5KSIrfeeqssW7ZM1qxZI7Vr1/bGsoIafYyBkvUsrlAmuljDPuw//pUbW8u7P2w3wbdupThz+85DabYQXPZ/U+4AACF0+E5LI6KiomTatGl5QrGqXLmyuS8yMpIyCgBekV/ZQoeGFeSVG9tIjfTtsvu7iWLJLri/cfnYKImJDJeYqHDp2LCCvHdLe6meEGsm2enQjqcGtZD7+zYxn1eD8svzNptyDgBAiLVr01KJrl27mjIKV/Q+3VFeuHCht5YFIIRZyxiWbzmUpy3b2h1HZeij78nSd+6XzPSzkpZ2Rhr3v0fCXEzx1JKJD/91YYE7wPblG9bdausIaABAiATjo0ePypkzZwq8TuuMjx075pU1AQhd1oN3unNbPjZajjq0Tzuw5ddc3SfS046JJTvLFoyjwiXXgI4LG1Y0x/YmzFpvqxfW4F3OISg7dp2gCwUAhGAw1prhH3/8UQ4ePOi0lEIdOHDAXGOtiwWAwgRbVyG0uAfvCupT7HigTmuJtV745ULsButa7eua6UIBACFYYzxkyBA5deqU9OjRw2mpxKJFi6Rnz56SlpZmrgWAwgRbDZn6uiiDMTRUr95+pFihWGuIy8dF5/qYqIhwU0JRmN1gDfBdkipJ1fhS5jVdKAAgBHeMx48fL999952sWrVKevXqJRUrVrS1ZdN2bYcPHzY9jC+88EJzLQDkp6AQmt+Ost5+zr4OIp9Q3HboBMkM++evSh3nrFJS8/YeLsxusK6BmmIACPFgrAfrfvjhB3n00UflnXfekUOHDpkXqzJlysjo0aPNxLz8DugBgKsQah+GdZKdtdWaY1mDs51cx1Bco1knGfHQGzLm8sby7g9/5+ozbC2icOw9bN+TWFuz6Rq0JVxxSj0AAEHcx9jxgN3atWtNn11Vo0YNadu2LYHYzehjjGDmuCPsWOPrSEsX3hrZznzM6r+PyLnMf3aMo8IssnzibXJy72ZbKN7+6w8SHZ27ZKIo9BCe/Vq0bIKdYgDwb17bMbanAfiiiy7yxacGECSclSTk1+FBw7PLA3dhYdLypqdl3fv3SunE6manuCSh2Nla6D4BAP7PJ8FYbdu2zZRSVKhQQRo1auSrZQAI4vIK7RYRER4mWdkW2bw/VU64GKah5cYxZRKlzahXpVJivNx/VbMSd8Og+wQABB6vdaVQWVlZ8tRTT0nVqlUlKSlJOnfuLM8995zt/qlTp5qd5L/++subywIQJBw7PujkuaRq5UytcUrquVwH7k4kb5Ss9Ny91aPj4uVkuthqiAvTDWPwGytM2YTjBDu6TwBA4In0Zii+8sorzQQ8HfvcpEkT2bBhQ65rOnXqJDfeeKN88cUXcv7553traQBCrLzCetAuvmaStLzpWYmILm27Tyfgafh1VQ/s+Hhaq2wtz7D/GLpPAEDg8dqO8dtvv23atXXv3l127Ngh69evz3ONtm9r0KCBCc8AUBy6c6s7uNoNQl9rdwhX3SeO7fhDdi2ZkecxCqpVdoYaYgAIfF7bMZ48ebIkJibKzJkzJSEhweV1upP8+++/e2tZAIKM/QE7LXdIjIuWDg0ryLqdx2T/5jV5+hTX7XZDnsfIrx7Y2pJNB4TYl2ZQQwwAgc9rO8abNm2S9u3b5xuKVXx8vKSkpHhrWQACYNfXsX43P5v3peZ6/2haugnFCae2yJ9THnY50U6Fh0mB9cBaIqHhuHWdBImJCjeT8DR4U0MMAIHPqzXGMTExBV63f//+Ql0HIHR2fVV+9br2nSKOn84Z6mFPd4rn2+0UN2rXTWr0/XeuUGw/1W7MB2vy7Tihn+vnbUfyjIQGAAQ2r+0Y16lTR/744498r8nIyDC1x+edd563lgXADxW1B7B9p4j0LEu+E+00FH88dXqeUKw7v8q+44Q+rjvWBwAIDF4Lxr169ZKdO3eacdCuvPHGG6a3cZ8+fby1LAB+yLFet6D6XVfBNHXvljw1xbpT/PgXG/Ncq7vFOw+lFepxi7o+AEBg8Fowvu+++0z98G233SZ33XWXrFixwtyelpZmxkOPHz/evFSsWFFuv/12by0LgB8qag9gx84TVnGV60pCvZZ5aortzszZ6OcobOClRzEABKcwi8WSXy97t1qyZIkMGDBAjh49KmFhYbnu02WUL19e5syZYwZ/oOT27NkjtWrVMm8nJydLzZo1fb0kwCMe/uyPXDW/9rIy0iV5+Uyp3XlQnvIJ+wl5M8Z1yjPVTgMvtcMAEDq8GozVwYMH5ZVXXpF58+bJ9u3bJTs724S3K664wuwq16hRw5vLCWoEY4QK7V5hPaRnyc6SsPCIXPdrtwkd3GFVPjbKhGC9LSEuSl6/qa1US/hnyAcAIDR5rSvF7t27zS6xBjUdA20/ChoASkJ3dzUY60G7rd+8JS1vekZKla9s7tPfTdmHYnX89D/t35rVLE8oBgB4t8ZYp9pdf/313vp0AEKox7HW/NZI3y5/fPKQnDq4Xda+f4+cS80prSjoV2J0lAAAeH3HuFy5clKvXj1vfToAQcZa/2s/cc5aPtElIUVmvnCHZP2v+0Ri9XoSFVuuUI9LRwkAgNd3jJs2bWrqXAGgOKy9iu3HMKu53y6Qy6/oI2fOnDHv9+3bVzqMeMrlQTsrnVpHRwkAgE+C8S233CLLly+XX375xVufEkAQcVbyoDXF6z5+UDLTz9pC8eeffy7RUdF5hndULhdjuk/oaw3En95+kZmmR9cJAIDXg/GIESNMD+PLLrtMnnnmGdm8ebOcO5fza08AKGrJg+NEuxrNOueE4uhoSaqeu4ziXGa2ZGVb5Ilrm0tStXImZL88b7MpzwAAwOvt2iIicrdPyo92r8jMzPToekIB7doQqBz7CY+6pL5MWrhN1u06Zk7TJZzaIjOeH2fbKdbhHbc+NlHuu6qZ+bjN+1NNED6ali72f8NFR4TlGhmtO8e6awwAgFcP3xUlf3u5tTIAP60nth6w27D3hBw5lW67f+efP9tCse4Uj3jodbn/qmZmF9j6cc7Yh2JFRwoAgE+CsQ7yAIDCcAysx9NylzzU6XGzXFC3nPz999+28glnH+c42MNxx5iOFAAAn9QYA0BhOQbWCIe/qRpVKyePP/28tLphgox891dbP2PHj2tTN8EcuIsICzOvXxvWxpRPVI0vRUcKAID3d4y3bdsmX3zxhezcuVNiYmKkVatWMmjQICldmklTAJzTgR3WHeCMrGzZ+vvPEhFdWuJrNTEBVwOtlk2s+Pt4rn7G9h+nIVmvc+w6QU0xAMAnh+9effVVuf/++yUrKyvX7TVq1JB58+ZJs2b8A+VJHL5DMLh03Bvy49v3SVh4pLQe8YIk1m4qFzWqKJv3pUrKyX862+gu8JSxHX26VgBAYPNYKcWyZcvk3nvvNd0lYmNjpXXr1tKgQQPTcUID28CBA6k7BpCvH374QX767/2mJVvWuTTZ8/OXkmWxmAN2+toe9cIAAL8Nxm+++abpLjFs2DA5cOCArFmzRrZs2SJr1641AVlLLObPn++pTw+ghLRmV2t3h05caavh9XYovvLKK23dJyo1vkiaXH2f7f6I8DDqhQEAgVFjvHLlSvOr+//+97+2E+OqRYsW8tprr5l/8H7++Wfp3bu3p5YAwI0t07xZn2sNxfZjnlsOmSArt+fUFCsd1EG9MAAgIILxwYMHTei1D8VWnTt3Nq9TUlI89ekBlJBj6zPr+47DN/TAWzk3jlV2Foq1JdvZrDDb561bKU4ys7LNbrYn1gAACE0eC8bp6elSvnx5p/eVK1fOdg0A/6SB07pTbH3f0zvJrkKx/oAdbfd5tLTDV7vZAIDg5bUBHwACi7PWZ9b3PTE9Tn/L1K9fP6eh2JGn1gAACG0eDcZ6wO7jjz8u1v033XSTB1cGoCBamuBsF9bVTnJJValSxZw/GDVqVL6h2JNrAACENo/1MQ4PDzet2YpDP07bvPkT7bu8evVq+fXXX83L/v37C9UfuG7durJr1y6X91944YXmEKIn0McYnuBYY+xsiEZJfPfdd9K9e3eXodgbawAAhCaP7RjXrl272MHYHw0ZMkROnDhR7I/Xvs1lypTJc7u2rgOCYSe5OA4fPiwVK1bMddvll1/u1TUAAODxYKwjoIPJgAED5LzzzpM2bdqYl8qVKxfp41988UWzewyEImedLNasXCr9+/c3LR0HDx7s6yUCAMDhu8L64IMPfL0EIGA5drLY/ucqmfnCHeag3dChQ6Vq1aqmfAIAAF8iGAPwOPuuEUf/XiuLP3lIsjLOmfcbtukqzdu09+HqAADIQTD2kg8//FCOHj1qDhVWr15dLr74YunatauvlwV4hbWLhIbi3z95SLL/F4orNu4oNfr+WyYu2kHNMADA5wjGXvLEE0/kue2CCy6QadOmScOGDYvddSI/1s4ZgK9pTbGWTyx2CMXNBz8u4ZFR9CEGAPgFgrGH9enTx4zAbt++vdSoUUMOHDggS5culYcfflh++eUX6datm6xdu7bIh/mUtRUb4O/0oJ3WFFvLJxq162Z2ijUUK/oQAwCCuo+xv7j//vtlzpw5Rf649957zwRaV6yt6IrbH1jLKtq2bWu6d9x5553y6quvFvkxitIOjz7GcEcnCW2T5o4xzx98Mt2UT9CHGADgT4J+x3jfvn2yefPmIn/cqVOnxJMSExPlrrvuMi9z584tVjDWsFtQKYXuVAPu6CShilMHHBERYfshzn6iHTXFAAB/E/TBeMqUKebFHzVp0qRQtcKusAMMT3Ks+y1uHbAeNNXJkdqv+KOPPsp3oh0AAL4U9MHYnx05csS8LluW+kr4bycJ+/eLS8OxvgAA4M/Cfb2AUPbpp5+a15Q7wB9pTXGXpEpSNb6Uea11wIWtKdYuLEF+fAEAEITYMfag2bNnm3IHPWRn7+TJk/LII4/YDgXec889Ploh4JoetCtqHbD9QbuMjAwTkItySBQAAF8iGBfSk08+Kd98802e2/v162ermWzTpo289dZbtvt+/PFHee2116R27drSvHlzKV++vDkM+Ntvv8mxY8ckMjJSXnzxRbn00ku9+rUAnuhcEXFog23Ms/r9998lKyvL/H8OAEAg4F+sQvr7779l1apVeW5ft26d7e1SpUrluq9///6mu4X2KV6zZo1p0aYhWoPyddddJ7fddpsJzECgd65wnGhn7T5BKAYABJKg72McyrTbhXUICH2M4QlDJ66UDWtXOA3FdJ8AAAQatnMAFJuWT9iHYp1oV9RQ7K5BIgAAlBTBGECx6EE7rSm2D8UrFn1d5J1idw0SAQCgpGjXBqDIMjMzZcyYMbnGPP+5/DupEB/ns0EiAACUFMEYQJHpoTqdZqd16yWtKXYcHFKSQSIAAJQEpRQAiqVBgwayYsUKqVKlSokO2mlNsbLWGBd2kAgAAO5GMAZQqENxHSudlGX7S8n2w2dth+SsXU+8PUgEAABPIBgDKPBQnLZke/KThySxUXtpdt0jHJIDAAQlaowBOGU9BGcd3pGVcU4O/bVUklfMynU/AADBgmAMwCktl3CcaFexcUep1XGA7X4AAIIJpRQAnGoTu0+emvKwLRT37NVbLhz2lOw4klNjzCE5AECwIRgDcDq847qB/SUzPaeWuEazziYU33dVM6bSAQCCFqUUAPKE4iuvvNI2vEPLJxpd+4is3H5cBr+5QibMWm86VgAAEGwIxgBsVq5cmSsU605x88GPS3hkzi7xuYxs06lCO1YAABBsCMYAbJo2bSrNmzc3b+tEuxEPvW4LxfboSAEACEYEYwA28fHxsmDBAvn3v/9txjxrTXGXpEoSE5n7rwo6UgAAglGYxWKx+HoR8Iw9e/bYJpMlJydLzZo1fb0k+CH9KyAsLKxIU/C0I0VZDuEBAIIMXSmAED9o99hjj8ns2bMlMTHR5XWMbQYAhAKCMRBEHHd27+6d5LK9mn33iUsvvVS+//57SUhI8PqaAQDwFwRjIIhoKNauEerAiZwexM52eh1bsmmZTVxcnJdXCwCAf+HwHRBEHLtFOOse4RiKtfuEHrSLjo722joBAPBHBGMgiDh2i3B8n1AMAIBrlFIAQURripV99wgrQjEAAPkjGANBxFX3CEIxAAAFo5QCCAFTp04lFAMAUAB2jIEQ8N///tcE41OnThGKAQBwgWAMhIDIyEj5+OOPJTs7m1AMAIALlFIAQWjx4sWyefPmPOGYUAwAgGsEYyDI6EG73r17S7du3fKEYwAA4BrBGAgi9t0nDhw4IC+++KKvlwQAQMAgGANBwllLtokTJ/p6WQAABAyCMRAE6FMMAEDJEYyBAEcoBgDAPQjGQAAjFAMA4D4EYyBA7dq1i1AMAIAbEYyBAFWnTh158MEHzduEYgAASo7Jd0AAe/jhh6VRo0bSv39/QjEAACXEjjEQQE6ePJnntkGDBhGKAQBwA4IxEEAH7erWrSsLFy709VIAAAhKBGMggLpPHD16VPr16ydr1qzx9ZIAAAg6BGMgwFqy9ezZU1q0aOHrZQEAEHQIxoAfo08xAADeQzAG/NTcbxfI5Vf0sYXiK3r3IRQDAOBBBGPAT3eKB/S/SjLTz5r3KzbuKG1vfJJQDACABxGMAT8tn7APxc0HPy47juS8DwAAPINgDPiZ48ePS0ZGRq5QHB4ZJedVLevrpQEAENSYfAf4mQEDBsinn34qkz+ZIi2vf9TsFGsovqd3kq+XBgBAUAuzWCwWXy8CnrFnzx6pVauWeTs5OVlq1qzp6yWhCPSpGRYW5utlAAAQMiilAPygpvjtt9/OczuhGAAA76KUAvCTPsVZWVkyduxYXy8JAICQxY4x4CfDOxYsWGDKJwAAgG8QjAE/mWg3c+ZMyicAAPAhgjHgZYx5BgDAPxGMAS8iFAMA4L8IxoCXEIoBAPBvBGPAC86ePSs33XQToRgAAD9GMAa8oFSpUjJnzhxJSEggFAMA4KfoYwx4SZs2bWTlypVSr149QjEAAH6IHWPAQzZu3CjZ2dm5bktKSiIUAwDgpwjGQAmknsmQCbPWy9CJK81rfd960K5t27Zy++235wnHAADAP1FKAZTAK/M2y9LNh8zbB06cNa+7JKTYuk9MmjRJWrVqJaNHj/bxSgEAQEHYMQZKYOuBk7neX/rTj3lasg0fPtxHqwMAAEVBMAZK4LyqZW1vH/17rfz03/tpyQYAQIAiGAMlcHfvJOmSVEnkwHr5c8rDkpmeU05BKAYAIPAQjIESKFc6ytQUr3jvAUIxAAABjmAMlMDSpUsZ8wwAQJAgGAMl0LBhQ6lVq5Z5m1AMAEBgIxgDJVCtWjVZvHixjBs3jlAMAECAC7NYLBZfLwKesWfPHttuZnJystSsWdPXSwIAAPBb7BgDRaAT7fr16yenT5/29VIAAICbEYyBIoRiPWg3d+5cwjEAAEGIYAwUIRRbu0/ExsZKZCQT1QEACCYEY6CIoZjuEwAABCeCMZAPQjEAAKGDYAy4QCgGACC0EIwBJwjFAACEHoIx4MTrr79OKAYAIMRwrB5wYvr06aYlW+nSpQnFAACECIIx4IQG4jlz5khERAShGACAEEEpBSAiS5culX379uUJx4RiAABCB8EYIU8P2l1++eXSrVu3POEYAACEDoIxQpp994mtW7fKc8895+slAQAAHyEYI2Q5a8n24osv+npZAADARwjGCEn0KQYAAI4Ixgg5hGIAAOAMwRghhVAMAABcIRgjZGzatIlQDAAAXCIYI2QkJSXJiBEjzNuEYgAA4CjMYrFY8tyKoLBnzx6pVauWeTs5OVlq1qzp6yUBAAD4LYJxEMvMzJQDBw6Yt6tWrSqRkUwABwAAcIVgDAAAAFBjDAAAAOQgGAMAAAAEYwAAACAHwRgAAAAgGAMAAAA5CMYAAAAAwRgAAADIQTAGAAAACMYAAABADoIxAAAAQDAGAAAAchCMAQAAAIIxAAAAkINgDAAAABCMAQAAgBwEYwAAAIBgDAAAAOQgGAMAAAAEYwAAACAHwRgAAAAgGAMAAAA5CMYAAACAiET6egEIDpmZmXLgwAFfLwMAAISgqlWrSmRkyWMtwRhuoaG4Vq1avl4GAAAIQcnJyVKzZs0SPw6lFAAAAICIhFksFouvF4HA545Siv3790v79u3N26tXr5Zq1aq5aXXwB3x/gxff2+DG9ze4Bcv3tyqlFPAn+j+jO36FYaVPTHc+HvwL39/gxfc2uPH9DW7V+P5SSgEAAAAogjEAAABAMAYAAAByEIwBAAAAgjEAAACQg2AMAAAAEIwBAACAHAz4AAAAANgxBgAAAHIQjAEAAACCMQAAAJCDYAwAAAAQjAEAAIAcBGMAAACAYAwAAADkIBgDAAAABGMAAAAgB8EYPpeamiqPPPKI9OnTRxo0aCDx8fESHR0t1atXl6uuukq++eabfD8+PT1dnn/+eWnZsqXExcVJQkKCdOvWTT7//HOvfQ3I37x58+Txxx+Xvn37mu9rWFiYedmzZ0++H1e3bl3btc5eOnTo4LWvAa7NnDnTPOf0uafPQX0uvvDCC5KRkeHrpSEfmzdvljfeeEOGDx8uzZs3l8jISPO8euqpp1x+jD6P83tO6sumTZu8+nXAualTp8pNN91kno+VK1eWqKgo8+9r+/bt5dlnn5VTp065/NhFixZJ7969pWLFilK6dGlp3LixPPTQQ/l+TLCI9PUCgJSUFPMXcZkyZaRZs2bmSRweHi7btm2TOXPmmJfbbrtNJk6cmOdjT58+LT179pQVK1ZI+fLlpVevXuaJ+8MPP8hPP/0k9957r7z44os++brwjyFDhsiJEyeK/fEDBw40/3840h+k4Ft33XWXvPbaayZUXXLJJeb7pM+/Bx54QObOnSsLFiww/7DC/0yaNMl874pD/55u1aqV0/s0fME/vr/6b2OTJk2kTZs2kpiYKAcPHpSVK1fKL7/8Ih988IH5d1I3K+y98sorcs8995gfcrp06SJVqlSRpUuXyjPPPCOzZs2SZcuWmcActCyAj508edKycuVKS0ZGRp77fvjhB0tsbKxF/1edP39+nvvvvPNOc1/z5s0thw4dst2+Zs0aS5kyZcx9c+fO9fjXgPyNGDHC8swzz5jvYUpKivm+6EtycnK+H1enTh1z3Y4dO7y2VhTel19+ab4/+lz79ddfbbfrc1Gfk3rfvffe69M1wrV3333X8n//93+WqVOnWjZu3Gi58cYbzffsySefdPkxjz32mLlGX8O//fzzz5YjR47kuf3w4cOWzp07m+/j9ddfn+u+tWvXWsLCwiwRERGWefPm2W5PS0uz9OjRw3zMwIEDLcGMYAy/N3LkSPNkvOeee3LdfvToUUt0dLS5b9myZXk+Tv9y1/s6dOjgxdWiMAjGweGCCy4w35+nnnoqz31Lly4198XExFiOHz/uk/WhaIYNG0YwDhFLliwx38fExMRct1977bXm9lGjRuX5mJ07d1rCw8PN/fqDVLCixhh+T39Fq2JiYvLUrWp9ce3ataVTp05Of32vfv75Z9m3b5+XVguEhr1795pfx9o/1+x17txZatWqJefOnTPPVQD+/e+q/nv6zf/O9Dh7TtepU8f2b+2XX34pwYoaY/g1/Yd3xowZptZJD27ZW7dunXndrl07px9bv359U1N19OhR+e233/LUUSFwfPjhh+b7mJmZab6PF198sXTt2tXXywpp1uefPsfq1avn9Bp9biYnJ5trBw8e7OUVwpPWrl0r48ePN89LrSlu3bq1+Tu6bNmyvl4aCnDy5ElziFL169fPdvuWLVvMuZ38/l3V27Xe2Pr8D0YEY/iVRx99VHbv3i1nzpyRHTt2mGCsHSpef/116dixY65r9X6lO8au1KxZ0/zFbb0WgemJJ57Ic9sFF1wg06ZNk4YNG/pkTaGuMM8/3TG2vxbBQw9W6os9Dcj6d7V2QoD/0AOw+ndldna27fCdhmM9rK4dnaysz9Py5cu7/AEnFJ7TBGP4Fe1A8fvvv9ve1xPuL7/8sowcOTLPtfrEVtoeyhVrJwNtCYfAoy389Ffy2l6oRo0acuDAAbNb8fDDD5sfmrRFmO5caSsieBfPv9CknWC0O8EVV1xhfrWuNmzYIM8995x8/fXXMmzYMImIiJAbbrjB10vF/+j3Z/Lkyblu01IJ/bfVvoMIz+kcBGOUyP3332/CbFG99957JvA40pIH65NOe2y++uqrMnr0aPn000/lq6++4td0Af79LSrHFn3a11hfNDC3bdtWdu7caf6R1v9PAHjejTfemOc2rTvV3eNx48aZvsh33323XHvttea3ffCPlor6on3F9Teys2fPNi1S58+fb2qFKUvLjWCMEtFDbRpgi6qgJuHlypUzvyrXBuX6a5233npLJkyYkKsnsTUkp6WlFfh59PHgP9/fktK6Vutf9voPMsHY+3j+wZHWrerf1YcOHZJVq1aZHrjwHzrgQ3f8tUex/jCj5YlDhw41f8drr3Ge0znoSoESmTJlirb8K/KL1jYV1ogRI5yegtWdQ6U/AbtinaxmvRb+9/0tLm1arwqangfPsD6n9HCdK9b7eP6FBv2B1VrWxPPSv1144YXStGlT8xxds2ZNrufp8ePHbWUVoficJhjD71nrnXRCnj2d5KOsT2pH27dvNwfvlJ6YRnA5cuSIeU15jW9Yn1P6fXB1EMf63LQ+VxHcsrKybBMueV4G3r+tSUlJEhsbm++/q6HwnCYYw+99//335nWjRo1y3a5z3LWGTXeMly9fnufj9BSu6tChA63agpDWnSs9mAfv044vWu5k/1yzp2NjdXdJ+6TqcxXBT88jaLsvba/pqt0X/MPhw4dtB92t/7bqv6d9+vRx+ZzetWuXGTGtrr76aglWBGP4nD4Bf/311zy366/kv/jiC9OBQOkhPHsJCQkyZswY8/Ztt91m20FU2qnA2obmoYce8vBXAE/QAyLO/r/QX/FpbbH1UKDWy8E3HnzwQfNaOxLoc85Kn4v6nFS33357rpPvCFy6CaHlVWfPns1znx6OHjVqlHlbO1JUrVrVByuEfScKPaPj7Hul/Yr1cKQO39GNo+bNm9vuGz9+vPnBRnvH6+E8K/2B5+abbza/FRg4cKA0btxYglWYjr/z9SIQ2oYPH25ayegOVIsWLcxhO/2HddOmTeYnVDV27Fhz2lmfsPb0yXrppZeavowalC+55BJzcEB3mfUEroaml156yUdfGayefPJJ20QlpQdzrL+Ot55c11/N6cEdKw2/r732mumTq39x6/8XehhQO5ccO3bMTG7Sw5h33nmnD74iWOmfv/au1YM9PXr0ML+e1eef1inqAZ+FCxeagz3wP/rDjPUHGPX333+bnUT9u1jbI1rp+Y5q1aqZ554+Z7Vll77Wa7TnvIawrVu3mmu7d+9ufmi1tvWCbyxevNh8L/T5qN8r/Z7qZDv94Ua/79rTWM9paPh17EX+yiuvmH879d9bHaakdePaJnP//v2m3EJ/G1SxYkUJWr6eSQ0sW7bMMm7cOEu7du0sVatWtURFRVliY2MtjRo1sgwbNsyydOnSfD/+3LlzlmeffdbSrFkzS+nSpS3x8fGWrl27Wj777DOvfQ3In34f9a+b/F4uvvjiXB/z448/Wm6++WZL69atLVWqVDH/X8TFxVmaNGliufXWWy1//PGHz74e5DZjxgzznCtXrpx5Dupz8bnnnjPPTfgvfY4V9LzUlx07dpjrDx8+bHnggQcsl1xyiaV27drm+ajPy2rVqlmuvPJKy7Rp0yxZWVm+/rJgsVhSUlIsTz/9tKVXr16WunXrmu9VdHS0+Te2Z8+elkmTJlnOnj3r8uMXLlxoPjYxMdESExNjOe+88yz//ve/LampqZZgx44xAAAAQI0xAAAAkINgDAAAABCMAQAAgBwEYwAAAIBgDAAAAOQgGAMAAAAEYwAAACAHwRgAAAAgGAMAAAA5CMYAAAAAwRgAAADIQTAGAAAACMYAAABADoIxAABB5L///a80bNhQypcvLwMHDpT9+/f7eklAwAizWCwWXy8CAACU3Oeffy7XXnuttGvXTurUqSPz58+XJk2ayKpVqyQ8nL0woCAEYwAAgsQ111wjhw8flsWLF5v3v//+e7n00kvlr7/+kqZNm/p6eYDf48dHAACCREREhGzfvl3Wr18vaWlpJhhbbwdQMIIxgJATFhZW5Jdu3bpJoNDdwsJ8TcePH/f1Uv3a3XffbcoP1qxZU6LH6dq1q/nzHj9+fKGuv/POO831vXv3dnr/hx9+aO5/4YUX8tw3duxYU1PcvHlzKVOmjDz77LPSp08fSUpKcvpYO3bskOjoaBk0aFARvyogOEX6egEA4G3Dhg3Lc9uBAwfku+++c3l/48aNJRA5+1qsNBDBuY0bN8qbb75pDq9pvW5J3HzzzbJ06VL5+OOP5emnn8539zY9PV2mTp1q+zhnZs2aZV7r2pyF8Ntvv11effVV874G+48++sjl56tXr56MHj1aJk6cKD/99JNcfPHFRf76gGBCjTEA/G+XtXv37ubtQP9rMZi+Fl+58sor5ZtvvjElCeeff36JHuv06dNSrVo1SU1Nlblz55rHdmXmzJlm97ZSpUqyd+9eiYqKynW/Pobepwfqfvvttzwfr9/v+vXry86dO6VGjRrmMd59910ZNWqUy8+pPxTWrl1bmjVrJmvXri3R1woEOkopAACws2XLFpk3b5506NChxKFYxcbGyvXXX28rg8iP9f6hQ4fmCcXq66+/NrvKAwYMcPrxCxcuNKG4UaNG8vzzz5vbNBjnp2rVqqZsY926dbJkyZJCf11AMCIYA0ABHn/8cVPTqa93795tfsVdq1YtE1yGDx9urtFfV+s11vcdaVjR++vWrev0/mPHjsljjz0mrVq1krJly5owpXWiTz31lNlx9AYtJdHdzMqVK5syi+rVq8t1113nssZW6671a9Idai0V6Nu3r9nNdPz1va5ff7XfuXNnSUhIkJiYGNNKTK+fNm2ay7ZjvXr1Mo+na9HdTw2LGzZscHr9r7/+atZas2ZNc325cuXMzqmWG8yePbtIfw5aVqA7r66+l8VZo3XHVneMtWuEM7q7u2DBgnzLKL744guXZRT2IXjkyJHmGu1lvHr1avnzzz/z/VqsX6t+7UBI01IKAAh1P/74o9YcmBdHjz32mLl9yJAhlsTEREvVqlUtAwcOtAwYMMBy7733mms+/PBDc82wYcOcPv6OHTvM/XXq1Mlz319//WWpVauWub9atWqWXr16Wfr27WupUqWKua1Vq1aW48ePu+VrceXhhx8214eFhVk6depkGTx4sPm8eltERITl/fffz/MxF198sbn/tttus4SHh1uaNm1quf766y2XXXaZZdq0aeaa3bt3m9v1utjYWEvPnj3NNV26dLHEx8fn+fPIyMiwDBo0yFwfExNjueiiiyzXXnutpWXLlua20qVLW7799ttcH7No0SJLVFSUuV+vu+aaayxXX321pX379uYxrrrqKktR1K5d2zzWtm3bnN5fnDWq5s2bm/tffvllp4/79NNPm/svvPBCp/enpaWZP8OkpCSn9x88eND8Oej3a9++fea2MWPGmMe844478v2aT5w4Yb6HcXFxlvT09HyvBYIZwRgAChmM9WXo0KGWs2fP5rmmuMH49OnTlgYNGpj7NJyeO3cuVxDSgKr3jRgxwi1fizMa4vTaUqVKWRYsWJDrvvfee8/cp4Fr/fr1ToOxvkycODHP42ZlZVnatWtn7tewnJKSkuv+M2fOWL755ptctz344IO2cLh9+/Zc982cOdOEvoSEBMuxY8dst3fv3t18zJQpU/KsQX+gWLlypaWwNAzrY1WqVMnlNcVZo3r11VfNx2lAdua8884z97/zzjtO7581a5a5Xz+/My+88IK5/8orr7Td9ssvv5jbdD36552fFi1amGuXLl2a73VAMCMYA0Ahg7HuFrvauS1uMJ40aVKeMGPv5MmTlsqVK1siIyMtR48eLfLX4upF12vVo0cPc9s999zj9PF0bXr/Lbfc4jQYX3LJJU4/7quvvrLtguvXUZAjR46Y3VYN6Hv27HF6je5O62O+8cYbttusO9KF/fPJjwZbfayuXbu6dY3q8OHDZodZ79PAam/JkiW2XfXU1FSnj6u/sdBr1qxZ4/T+Ro0amfu//PJLp4H3k08+yfdrt/4Q9tprr+V7HRDMqDEGgELSCWLx8fFufUztfKC0PtYZ7UWr7cIyMzPll19+KVa7NmcvDRs2NPfr4y5fvty87aqm1lrv+uOPP7qctuaMjiNWQ4YMMV9HQfTxz5w5I506dTL1us5Y+0mvWLHCdlv79u3N6xtuuEGWLVtmvqbiOnjwoHldoUIFt67R+pj9+/c3b3/wwQe57rO+r+OctcbckR640/9XtEa9bdu2ee7XVmt6aFDrwx27Xmi9sXrvvffy+cr/+ZqtfwZAKKKPMQAUkquDcyWhU8rUjTfeaF7yc+jQoSI/fn49bNWRI0fk7Nmztp62zjRo0MB2OKwofy67du0qUg9o65+FTmvTQ32F/bPQIRZ//PGHfPvtt+aldOnS0qZNGxNQNSxra7PCOnHihHmth/fcuUb7HzJmzJgh06dPl5dffllKlSolp06dMm3arPc7s2jRIrM2V/dbD93p/0ORkbn/adcDgffff78tPGvHCmesX7MeBAVCFcEYAApJA1dxZWdn53u7djeoUqVKvo+hnRyC7c/F2Z+F7mbrjmx+7MO2thvTzhka/DRA6g74qlWrzOtnnnnGBOcHHnigUGvQLg7WfsHuXKNVjx49zPdRf2j48ssvZfDgwfLZZ5+Z8c0aWLt06VLkoR4aZK33a5u5n3/+Oc812kFFd51119jZxDz7Hwq0cwgQqgjGAOAG1ilyJ0+ezHf31JG2fdu0aZPZCXRVkuBJ+utzbZ927tw5sxvaokULl7ukrkoHXNGhEUq/vsLQPwul44sL2ul2ZB3bbS1j0F1wfQwdkfzggw+aP1vrznd+tBTBupPu7jUqbWU3YsQI0/pPyyc0GFvLKPR2Z7KyskzLOR0S0rFjxzz3f/LJJ7Zdf53Yl5/Jkyeb6XvOeiRbv+aCfkADghk1xgDgBtbQ6CoEWmuJHV1xxRXmte4a+oL+2l37CytXQc8a3KzT9ApLd8GVlg3ojmhBdDdVf8DQvsgpKSlF+lyOtETh1ltvNUFfd3m11KIwtAQjv4DpjjVqANaA/MMPP5iBHLqzrWOiXY3v1p1wDa1XX3210/INaxnFpEmTTP9lZy/6Z6AlL7rmOXPmOP08OuVPOathBkIFwRgA3EAPgGmNpg530B08e1o/+vrrrzv9uNGjR5tfres1+ut+ZzvOOrK3oOllJXHvvffagpXWztrTsKxBSncY77zzziI9br9+/aR169ayb98+c6jMcRdWdzm1JthKdyrvuOMOE6J1+IezoRS6s63rsf8B5MUXXzSDVxzpNVu3bi1SGYoOBdGdbq0P3rZtW577i7tGe/r4PXv2NGFVa6CVTp7THeGillFo2YQGWt31d3WAU2mg1lpj5ez/JS2j0P939ZCk9TAjEJJ83RYDAAKlXZu+zs8rr7xie4yOHTuaQRPnn3++GZrxyCOPuBzwof2B69ata+4vX768aRWmrbn69+9vWpHpx+uwD3d8LYUZ8NG5c2fz+du0aVOoAR/6+VzZuXOnGUhhbUWm/Yy1LZh+ja4GfFjbkunAidatW5thKtddd50ZPKIDKPQ++wEa+jh6W+PGjc1gD/34bt26mRZ3evtNN91kKYpx48aZj3vrrbec3l+cNTr67LPPcrXP09Z2zmRnZ5t2dxUrVrRkZmbmuX/kyJHm43XASEG2bNliW/OuXbty3ffFF1+Y+3RwCRDKCMYA4KZgrCZPnmwCpfa5LVeunOnxu3Dhwnwn3yntXasDGjRQazjWgRoaiC644ALLfffdZ1mxYoVbvpb8aJDr3bu3pUKFCiZU6oQ/DVyrVq1yen1hgrHSHsbPP/+8+VrKli1revnqn0O/fv0sn376qdOPmTdvnpksWKNGDfNnoX8mTZo0MVPzdKqeDj+x0sEeOgClWbNmpte09fGvuOIK09NXw2VRbN682fyAoJPz8lOUNTrSQS4advXPT3/o0bDtzLJly8w1GoCd/T9jDeFff/11ob62Dh06mOsfffTRXLfr90Jv/+mnnwr1OECwCtP/+HrXGgAAf6K9gLUuXGuTmzdv7rN1aJmLtnXTtWi5hSdoqY6WdzRr1kzWrl3rkc8BBApqjAEAcKAtzfRg4oQJE3y6Dm35ph0sdLiMpzz55JOSkZFhAjgQ6tgxBgDAibvvvlteffVVM3FQpw8GI23Fp+FbJ/L5qjMK4E8IxgAAAAClFAAAAEAOgjEAAABAMAYAAAByEIwBAAAAgjEAAACQg2AMAAAAEIwBAACAHARjAAAAgGAMAAAA5CAYAwAAAARjAAAAIAfBGAAAACAYAwAAADkIxgAAAADBGAAAAMhBMAYAAAAIxgAAAEAOgjEAAABAMAYAAAByEIwBAAAgEPl/wn6M9bihJXIAAAAASUVORK5CYII=", "text/plain": [ "

" ] }, "metadata": { "image/png": { "height": 302, "width": 355 } }, "output_type": "display_data" } ], "source": [ "from graph_pes.utils.analysis import parity_plot\n", "\n", "carbon_model = models.load_model(\"c-gap-17-painn.pt\")\n", "parity_plot(\n", " carbon_model,\n", " amorphous_carbon_structures,\n", " \"forces\",\n", " lw=0,\n", " s=3,\n", " units=\"eV/Å\",\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that this model generates stress predictions as expected - `graph-pes` automatically calculates stresses for all models, allowing them to be used in NPT simulations." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.07337982, 0.06381519, 0.08241484, 0.00042613, -0.01667964,\n", " 0.00498147])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "carbon_model.ase_calculator().get_stress(amorphous_carbon_structures[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with a low-density amorphous carbon structure:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", "

\n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "starting_structure = amorphous_carbon_structures[0].copy()\n", "load_atoms.view(starting_structure, show_bonds=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting cell:\n", "[[9.1 0. 0. ]\n", " [0. 9.1 0. ]\n", " [0. 0. 6.83]]\n", "\n", "Final cell:\n", "[[7.86 0. 0. ]\n", " [0. 7.86 0. ]\n", " [0. 0. 5.9 ]]\n" ] } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "# negative sign shows that we want to **compress** the system\n", "desired_pressure = -1000 # in bar\n", "bar_to_eV_per_A3 = 6.2415e-07\n", "\n", "print(f\"Starting cell:\\n{starting_structure.cell.array.round(2)}\\n\")\n", "\n", "# run NPT\n", "final_state = torch_sim.integrate(\n", " system=starting_structure,\n", " model=carbon_model.torch_sim_model(device=DEVICE),\n", " n_steps=1_000,\n", " timestep=0.001,\n", " temperature=300,\n", " integrator=torch_sim.npt_langevin,\n", " external_pressure=-desired_pressure * bar_to_eV_per_A3,\n", " trajectory_reporter=torch_sim.TrajectoryReporter(\n", " \"test.h5md\", state_frequency=10\n", " ),\n", " b_tau=0.1,\n", ")\n", "trajectory = torch_sim.TorchSimTrajectory(\"test.h5md\")\n", "cells = trajectory.get_array(\"cell\")\n", "print(f\"Final cell:\\n{cells[-1].round(2)}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the NPT simulation has compressed the system as expected" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoYAAAIdCAYAAACgK7z5AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAewgAAHsIBbtB1PgAAZwJJREFUeJzt3Qd0XNXV9vGtXqxeLMmW3Hu3cQGMCQQwmBZ6SygJJJCQEAgkAZKskBdC2kdeEvJCgECA0HuvJjQbMLjbgHuRbMtNslWsLs239jEjz71TNJJmNO3/W0uW5k7R1YwlPTrn7H3iHA6HQwAAABDz4kN9AgAAAAgPBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYBEMAAAAYiQffoafa2tpk586d5uPi4mJJTOQpBQAAkYkRw17SUFhWVmbenAERAAAgEhEMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYBAMAQAAYCSG+gQQOPvqW+T9r3ZLfkayHDGqQJISyP0AAMB/BMMo0dTaLtc/tkzKqxrM5cEF6XLtvNEysSwn1KcGAAAiBENKUeKDr3Z3hkK1dW+DXPefZfLX19dIbWNrSM8NAABEBoJhlHh39S6Px19fXinfu3eRvPuF5+sBAACcCIZRoKq+WZZv3ef1+v0NrfKHl76UV5Zu79PzAgAAkYVgGAXe+3K3dDgOXU5JjJfUpAS32/3rvU3S2NLetycHAAAiBsEwCvzXNo189Nj+8uCVM+WIkQWW4wea20zVMgAAgCcEwwhXUdUg63bWWY4dN75I+melyq3nTpSZw/Ms173KdDIAAPCCYBjh3l2903I5t1+yTB1yqEXNadMGWq5fW1kn6yqtQRIAAEARDCOYw+FwqzY+dlx/SYg/9LLOHJ4vhVkpltu8soxRQwAA4I5gGMG+2lErlfubLMe+Ob7IcjkhPk5OmTLAcuy9L3ZJfVNbn5wjAACIHATDKOpdODAvTUaXZLrd7qTJJRIfF9d5uam1w20KGgAAgGAYodraO8xuJ/aikziXAOhUkJkiR46yVii/smyHmYoGAABwIhhGqKVb9pnG1fZg6M1p06zTyVv2HJAvttUE7fwAAEDkIRhGyTTymAGZMjAv3evtpw7JlQG5aZZjry7bEbTzAwAAkYdgGIF095KF6/Zajh03vtjnfXSN4alTraOGOhVd09ASlHMEAACRh2AYgdbsqJWm1nZL6DtmbP8u73fipGJJSji0BrG13SFvraQIBQAAHEQwjNDdTlwN7d9PcjOSu7xfdnqyHD3GGiDf/5It8gAAwEEEwwi0fZ81GJb5WFtoZ+9zWLm/MWDnBQAAIhvBMAJtq25061/or+LsVMvluqY2aWnrCNi5AQCAyEUwjEDbqq0jhqXdGDH0NOW87wAFKAAAgGAYkY2td9q2wSvtxohhZmqipQBFVdc3B+z8AABA5CIYRpidNU3S3mHdscRX/0I73Rklt5911LCaEUMAAEAwjPxp5Ky0JPPWHXm26eTqeoIhAAAgGEac7b0oPHHKy0ixXCYYAgAARTCMocITpzz7VDLBEAAAEAwjv1VNdwpPnFhjCAAAPCEYRpjtthHDgbndHzHMt60x3HeAqmQAAEAwjCjNre2yu9Ya4krzezBiaAuGVUwlAwAAgmFk2b7Pffu6gbndD4b5tuKTffUt4nBYW+AAAIDYQzCM4IpknRJOS07s9RrDtg6H1Da29fr8AABAZCMYxlhFsqdgqNgWDwAAEAxjMBgmJ8ZLZpp1pJFt8QAAAMEwxppbe+1lyIghAAAxj2AYgyOGit1PAACAHcEwQtQ3tcn+htZeN7f2NmLIGkMAAEAwjBDb91lHC+PjREp60KrGKY9ehgAAwIZgGKFb4RVlp0pSQnzAguE+ik8AAIh5BMMIsa0qcOsLPU0lM2IIAAAIhhG660lvKpI9FZ+wxhAAABAMY7Ai2dOIYV1Tm7S0dfTqMQEAQGQjGEYA3cc40MEw17bGUDFqCABAbCMYRgBtU9PQ3B7QqeTM1ERJSoizHGP3EwAAYhvBMALYRws10PXPSu3VY8bFxbntmUyTawAAYhvBMAJb1ZTkpEmCNjLsJbfdT5hKBgAgphEMI8B2+/rC/N6tL/TWy5ARQwAAYhvBMAJHDHuzFZ6vymSCIQAAsY1gGIFrDAfmBmnEkKlkAABiGsEwzHU4HLJjX1+NGFKVDABALCMYhrk9tc1ujacDt8aQ4hMAAHBIosvHiIBp5LTkBLeRvkBNJe+rbzHNtLWVTW/oY7y5cqe8taJSBhf0k+9/c5hkpCb18mwBAECwEQzD3HZb4cnA3LReBzcnex/Dtg6H1Da2SXZ670LcEx9vlQc/2Gw+Xr2tRlKT4+WHx4/s1WMCAIDgIxiGuVkj8uWm1HGmZY1WJxdn966xta9g6NwWrzfB8KXF2zpDodPnG6vlh8f3+CEBAEAfIRiGuaLsVPMWDMmJ8ZKZlih1jW2dx6rqm2VIYb8ePd781TvlrrfXux3fvq9R2to7JDGBJa0AAIQzgmGMy++XYgmGOmLYEwvX7ZE/v7LG43XtHQcrqwcVuAfO5z6rkDdXVMrwogz5zlFDpDQvMIU1AACg+wiGMS43I1m27D3QrSbXOvqnFcx765qlqq5FduxrkIc+3Gxa63hTXtXgFgy/2FYj98zfYD7evOeAfLhmj3zvG0PlzBllvdryr7Gl3RTpAACA7iEYxrju7n7y8Ieb5dnPKkz48iUpIU5a2w8FxXINn6MLLbf5ZP1ey2Vty/PPdzfKB2v2yM9PGeNxhLGrwHrLc6tl0YYqGVWSKb8/b5LkBKiCGwCAWMCirxjXnd1PVlfsl/8s2NJlKJw3uUROnFTiNmJot2l3vcf7f7W9Vq58YLE89clWn6OQdi8v2S6fbqgSvcfayjp5elG53/cFAAAEw5jnFgx97H6ioasrR48plGvnjTb9C7sOhoemsO1a2zvk/vc2yZ1vrPU7HOpUtKvVFTV+3Q8AABzEVHKMc9v9xMdU8qqK/W7TxfkZKZKfmWzeTyzLkdMPG2DWBw4qsBaRVFQ1WJpn1zS0mjWKrvIzkqXK9vlfX14pSQnx8uO5I332b9TH+3J7jVsYDUTDbgAAYgXBMMbZ1xh6q0puam2XtTvqLMf+59yJMmNYvsfbD7Jt26fTz3vqmqV/VqrHaWQNf/++cpY88tEWU6nsOkb40pLtkpgQJ1cdN8JryPtsY5V02AYW65vazNdjD78AAMAzppJjnD0Y1jW1SUub+xrCNdtrzc4oTlo0PH5gttfHLchMkXRbZbCOGjrZg+GQwnRJT0mUq44fIbecPcGtKvm5z7bJv97fZEYAPbEXsviawgYAAJ4RDGOcfY2h2neg1e3YSts08ojiTBPkvNGRvTLbqGH5Xu/BcFj/jM6PZ48ulJu/Nc6ET1dPfVJuRhQ9rUdcvKna43m4fk4AAOAbwTDGZaQmmmlcV7r7id2qcmswnFSW0+Vj24PhVpd+ifZgONwlGKpvjO0vvzxtrNgnjrUq+p1VOy3HVpbvlwYvldLlVd4LXAAAgBXBMMbpyJ591HCfrQBER+S+3F5rOTZxkPdpZCd7H0LntG57R4ds2dPgdcTQ6bgJxXLDqWPcjj/4wSbT89BJ+xZ6w4ghAAD+IxhCcu1Nrm0FKOt31kmzSxBTE0q7HjEcbBsxdK4x3FbVaMKmq6EegqHSfojXnjTKcmxPbbO8tbLSfKxrDr2tL3T9nAAAoGsEQ3joZdjicxp5SGE/yU5P6vZUslYI1zW2ykbbNLIWqvh6vFOmDpBJZdYRysc/3mpGDXVEsHJ/k9f7aiV0Q/OhvaABAIB3BEN0uS3eSlujaH/WF6oBuWlu1cU6nexeeNKvy+nui+cM9Thq+MmGvW69EO1FK1QmAwAQBsHwsssuM7/Ufb01NXke7VmyZImce+65UlRUJKmpqTJ06FD5yU9+Irt37/b5OXft2iU//vGPze1TUlLM/fVxli5dGqSvMhq3xTtUfNLe4XDbQcSf9YUqMSFeBuamWY7pnsluhSdFnqeRXU0ZnCMTPYwaLlhr3e3kiJEFUpJj+5wEQwAAwqfB9ezZs2XEiBEer0tIsPa6U88++6xceOGF0tbWJjNmzDAhb/HixfKPf/xDnnnmGVmwYIHHx1u3bp3MmTPHhMdhw4bJGWecIZs3bzaP9+KLL8rTTz8tZ555ZlC+xkhmbwCtawCdO4Zs3lMvB2xTsbrDib+0AMU1mB0cMbRWCg8r7DoY6rlcMmeo/Pzx5ZZRQ32zB0PdUWX7vsbOYxVUJgMAED7B8IorrjCjh/7YsWOHXHrppSYU3nvvvfKDH/zAHG9vbzeP8eijj8pFF10kixYtsuyCoUHmggsuMKHw4osvln//+9+dofO+++6TK6+8Ui655BJZv369FBcXB+krjUwDcg7uRuJUUd0g81fvkhMmFsuqcutooY4A6ppAf9l3QFm9rcZtK7xhfowYuo4arvKyB3JqUrxMHZIjK8r3WfZ1pjIZAIAIXWN45513SkNDgxx//PGdoVBpyLvnnnskOztbPv/8c3n77bct93vjjTdk2bJlkpOTI3fffbdlJFIf57jjjpP6+nr529/+1qdfTySYOChHSmzh8P7/bjQjhdoj0H7b7rDvmfyVre2N9lAszbNO/fpca3jUEK/XTxuaJ8mJCTIo39Ymh2AIAEBkBsMXXnjBvNdRQbuMjAw5/fTTzcfPP/+8x/vp9Xo7O+fj2e+Hg+HsRyeMdGtZ8+iCLbLKtuOJfZ1fV+whzW5oYT9JiPf/v+HUIbkyodTzORwxIt9jGNVpZXt7HAAAEKJg+N5778n1119vRu5uuukmE+Kam91316irq5MNGzaYj6dPn+7xsZzHdXTQlfNyV/fTqeQDB1hzZnf4iHyZOTzPcuzZRRWyv6G1x+sLVVm+79FAf6eRrWsNPY8aznIGQ9v0dYfDITtc1hwCAIAQrjF85JFH3I6VlJTIgw8+KCeddFLnsS1bDu2DO2jQII+PVVZWZt5rUYkr5+Wu7qdrEfXzjB8/3q9z37Ztm8/rKysPNlqOdBq4fnT8SFm6+TNp63CYYwf/PUTXFtqnnLuSlpwohVkpbkUiTsMKfY8o+ho11PWKTmMGZHYW0WSkJpkWPK6NunU6ebBtJxYAANCHI4aTJ082a/pWr14ttbW1ppWMrg088sgjTaDSad/333/fMmLo1K+f51/izmlifTxXzvt2dT9P9/VFA6Wvt5kzZ0q0KM1Pl7NnHQzQnug0smvBj7/sI3i9GTFUeg5XHDtMEl0aFp45w3re9ulk9kwGACDEI4bXXXed5XJmZqaccMIJprBE28a89NJLcu2118ry5YdakCC0vjN7sMxftVOqbE2u1aRuFp446Ujdks37PF7naY9kf0woy5E/XzRFFq7bI+MGZsvRYwrd1jYu33pofSQFKAAAhMlUsqcRn9/97ncmGK5YsUIqKirM6JsGRyddB6gVyHZaWayysrIsx/W+1dXVXtcPOu/n6b6+6Ln5oiOf0TRqqFO/P/jmcPnDy1+5XefvjiddbY3nVJiZIllpXW+t540GVW9hlRFDAAAiJBiqsWPHWtbxaTAcPHhw57Hy8nKZOHGi16A2ZIi1AEEvazDU+3nivJ+GUtfP05XS0lKJNd8cXySvLN1hWcOXnZbkFrZ6O5U8tIut8HrDHkYrqhpMEUp8D6bCAQCIFSFrV1NVdagBsXOkUEfynDua6E4nnjiPT5s2zXLcebmr+40cOdJjOxscouH5xyeONG1snLTZdU/WFzp3P/HEn63wespeaNLU2uG1AAYAAIQ4GD755JOdYXD06NGdx51b1j3++OMep4NfeeUV8/FZZ51luc55v5dfftnjdLLz8ez3g2cjijLlzxdOlrkTi+XSOUPk0qOH9vixctKTJDM1MWDrC/2Rn5Es6cnW7RaZTgYAIETBUAtKNKTp1nauOjo65IEHHpCbb77ZXL7mmmskKenQOjMtRklPT5f58+fL/fff33lct8T70Y9+JPv37zf7J8+dO9fyuPPmzZOpU6ea6/V2ensn3RLv3XffNSOFP/3pT4P1JUcd3eXkF6eNlYvnDJU0W8jqDh1p9DQNHcxg6OlzUoACAECI1hhqr0AdxcvNzTXTvEVFRSa0aesa5zrACy+8UH77299a7jdgwAB56KGHzHXaEFtDpK4f1G3wNm3aZB5HR//s05p6+YknnpA5c+aYvokLFiwwAVL7G3722WeSmJhojrNPcmholfAX2w61CUpO9H8rvN58zjU76izrDAEAQAhGDLWHoY7+aSPpNWvWmK3odNROnXPOOfLaa6+ZgKeBze7cc8+VRYsWmWlfDYO6U4qOAF599dWmitm5DtFOp6RXrlxpbqe31/tpMNTH0cdzTjej79lH74YUdG8rvEAUoJTvZSoZAABf4hy6FQh6zFlR7ax8jsUqZn/sqmmSS+75VNq/3lXlR8ePkLNmem+mHQja4/C3z662rHV89tqjgvo5AQCIZCFrV4PYUpSdKnddepi8+8VOGVqYIXMnBX9Kf3C+tTJZ932uaWiV7PSe904EACCaEQzRZ0aVZJq3vlKSm2q2zXPu/awqqg5IdnrPGnUDABDtQtauBgg2XcM40FbgUk4BCgAAXhEMEdW0MtnVVgpQAADwimCIqDak0BoMV5bvD9m5AAAQ7giGiGqTB1vXE67fWS9V9WyNBwCAJwRDRLUJpdluW+N9vrE6ZOcDAEA4IxgiqiUmxMthQ/Msxz7bWBWy8wEAIJwRDBH1Zo3It1xevLla2to7QnY+AACEK4Ihot6M4dYRw4bmdlm9rSZk5wMAQLgiGCLq5WekyMjiDMsxppMBAHBHMERMmDncOp382QYKUAAAsCMYIibMsgXDLXsPyK6appCdDwAA4YhgiJgwekCWZKUlWY4xnQwAgBXBEDEhIT5OZgyjbQ0AAL4QDBEzZtra1izbsk9a2to7L2+rapBXlm6XzbvrQ3B2AACEXmKoTwDoK9OH5kmciDi+vtzU2iEryvebBthPf1ouD76/STocIkkJcfLX70yVsQOzQ3zGAAD0LYIhYkZ2epKMHZglX26v7Tz23y92y8tLdsgn6/d2Hmttd8gbKyoJhgCAmEMwRMztguIaDN9ZtdPj7SqqGvrwrAAACA+sMURM9zP0Zlt1Y9DPBQCAcEMwREwZUZQh+RnJXd5u34EWOdDc1ifnBABAuCAYIqbExcXJDA+jhuMGZkm8Vqa42M6oIQAgxhAMEXPOnVUmyYmH/uufPbPUVCH3z0q13G57NesMAQCxheITxJzBBf3knu9Nl8WbqmXSoBwZWZxpjg/MS5OdLtvkbSMYAgBiDMEQMRsO9c1VaV66LNm8r/Py9n1MJQMAYgtTyYBLMHTFiCEAINYQDIGv6VSyvWWNw+HcJwUAgOhHMAS8jBjWN7VJbWNryM4HAIC+RjAEvlaUnSIJtp41tKwBAMQSgiHwtYT4eBmQY59OZp0hACB2EAyBLtYZAgAQKwiGgAsqkwEAsYxgCPgYMaSXIQAglhAMAR8jhlp8QssaAECsIBgCPkYMm1rbpaq+JWTnAwBAXyIYAi4KMlMkOdH6bbGddYYAgBhBMARcxMfFycBcKpMBALGJYAjYDLSvM9zHiCEAIDYQDAGbUnoZAgBiFMEQ6GrEkDWGAIAYQTAEuhgx3LGvUdo7aFkDAIh+BEPApjTXOmLY2u6QPbVNITsfAAD6CsEQsMnplyTpKQmWY6wzBADEAoIhYBMXF+c2asg6QwBALCAYAn7sgMKIIQAgFhAMAX/2TKaXIQAgBhAMAQ/oZQgAiEUEQ8CPXoY79zdJa3tHyM4HAIC+QDAEPLDvl9zhcJhwCABANCMYAh5kpiVJdlqS5RiVyQCAaEcwBPyuTCYYAgCiG8EQ8LMyedPuAyE7FwAA+gLBEPBiZHGm5fLyrfvE4WDPZABA9CIYAl5MHZJruby7tlkqKUABAEQxgiHgxeCCdMntl2w5tnTLvpCdDwAAwUYwBHzsmTxlcI7l2HKCIQAgihEMgW5MJ7POEAAQzQiGgA9TBluD4f6GVtmyh+pkAEB0IhgCPpTkpEpRdqrl2DKmkwEAUYpgCHSxztA+nbxs6/6QnQ8AAMFEMAS6MNVWgLKifJ+0d3SE7HwAAAgWgiHQhcm2dYYNze2yfmd9yM4HAIBgIRgCXSjITJFB+elu1ckAAEQbgiHQg+pkClAAANGIYAj4wV6AsrqiRlraWGcIAIguBEPAD5MG5Uicy+Xmtg5Zs6M2hGcEAEDgEQwBP2SnJ8nwogzLMaaTAQDRhmAI9GJ7PAAAognBEOhhAcpX22ulsaU9ZOcDAECgEQwBP00sy5aE+EMrDds6HLJ6G7ugAACiB8EQ8FN6SqKMGZBpObaqvCZk5wMAQKARDIFuGDcw23J5+77GkJ0LAACBRjAEuqE4O9VyeVdNU8jOBQCAQCMYAt1QRDAEAEQxgiHQi2C470CLNLdSmQwAiA4EQ6AXwVDtrm0OybkAABBoBEOgm5XJmWmJlmNMJwMAogXBEOimoizrqOFOgiEAIEoQDIFeVibvJhgCAKJEnwbDX/ziFxIXF2febrvtNrfrb7nlls7rvb2tWbPG6+Nv2LBBLrvsMiktLZWUlBTzXi9v2rQpyF8ZYkl/WzBkxBAAEC2si6WC6OOPP5Y77rjDhDuHw+HztpMnT5YpU6Z4vC4729pg2GnhwoUyd+5caWhokPHjx8tRRx0lq1evlocfflieffZZmT9/vhx++OEB+VoQ2+hlCACIVn0SDDWs6chdSUmJzJgxQ1588UWftz/jjDPM6GF3Hv+8884z72+66Sa5/fbbO6+7+eab5Q9/+IO5fu3atZKWltarrwWglyEAIFr1yVSyhrX169fLfffd53XErzceeugh2bFjh4waNcptilov6/GKigp55JFHAv65EXuKcqzBsKquWVrbO0J2PgAAREwwfP/99+Wuu+6SSy65RE4++eSgfI4XXnjBvL/gggskPt76Jenl888/33z8/PPPB+XzI7arknVhxB56GQIAokBQp5Lr6+vle9/7nhQVFcmdd97p9/2WLl0qN954o1RXV5sRxqlTp8ppp50mmZmZHm+/bNky83769Oker3ced94O6I2M1ERJT0mQhuZ2y3TygFyWKQAAIltQg+ENN9wgmzdvNiN6ubm5ft/vlVdeMW+uNCD+/e9/NyOPrurq6qSqqsp8PGjQII+PV1ZWZt7v2bNHDhw4IP369fP7XLZt2+bz+srKSr8fC9FBC6h01HDzngOdx6hMBgBEg6AFw7ffflvuvfdeM72rxST+GD58uCkcmTdvngwePNgc+/LLL+WPf/yjvPrqq3LppZdKQkKCfPvb37YEQydvgS8jI6Pz49ra2m4FQ2eoBFwV51iDIb0MAQDRICjBsKamRi6//HIpLCw06wv9dfHFF7sdmz17thk9vOaaa8xjXXfddXLuuedKcnJygM8a8F9/dj8BAEShoATDa6+91kzBPvXUU1JQUBCQx9T2NXfffbeZDl60aJHMmTPHHHddd6jTxN7WOjplZWV16/NqNXNXU8kzZ87s1mMiOkYMXdGyBgAQDYISDHVNYWJiogly+ubKuXPJAw88YJpOFxcXy5NPPtnlY+bl5Un//v1NEHNd96fBUK/TQpXy8nLTHNtbuNOQ2p1pZKW7pwBdVSYTDAEA0SBoawzb2trkgw8+8Hr9li1bzJtzLWFX2tvbzRS1slcnT5s2zYTMxYsXm+plOz3uvB0QjF6G2q6mvaNDEmztkgAAiCRB+S22f/9+s+2dpzctIFG33nqruazh0B8vv/yy2dlEK0LtbWnOPPNM815HHjs6rI2G9bJOaauzzjorQF8hYp19xLDD4ZC9dS0hOx8AAAIhbIY3dBr40UcflaYm9yk53ULviiuuMB9rRbJOP7vS7fYGDBgg69atk9/85jeW6/SyHtcpYXurG6CnstOTJDXJ+u3DdDIAINL1yV7J/tA1glqV/MMf/tA0tB44cKA0NjaadjW6nZ469thj5Z577nG7b3p6ujz99NMyd+5c0+5GRxcnTJggq1evNm+6rvCZZ55hn2QEjI5ca2VyeVWDpTJ5UkjPCgCAKBkx1H6Bv/zlL02F79atW02Lmtdff91UFJ966qny+OOPm3WErj0J7W1tVqxYYUYFNWQ+99xz5r1e1uOHH354n39NiK3KZHoZAgAiXZxDF/qhx7RC2tkEW6ufqWKOHXe+sVZeXbaj8/JJk0vkhlPGhPScAACIihFDINLQyxAAEG0IhkAP0csQABBtCIZAwHoZNpm2NQAARCqCIdBDRdnWYNja7pDqenoZAgAiF8EQ6KHcfsmSlEAvQwBA9CAYAj0Ur70Ms1MsxwiGAIBIRjAEeqHYNp2sTa4BAIhUBEMggOsMaXINAIhkBEMggMGQEUMAQCQjGAIBnEpmxBAAEMkIhkAv9LcFQy0+YZdJAECkIhgCARwxbG7rkP0NrSE7HwAAeoNgCPRCXkaKJMTHWY7RsgYAEKkIhkAvaCjsn0UvQwBAdCAYAr1EZTIAIFoQDIEAB8PKfY0hOxcAAHqDYAj0UmleuuVyeVVDyM4FAIDeIBgCvTSowBoMKwiGAIAIRTAEemlQfj/L5X0HWqSukZY1AIDIQzAEeqkkJ9WtZQ2jhgCASEQwBHopMSFeBuamWY6xzhAAEIkIhkAADMq3F6AcCNm5AADQUwRDIADKKEABAEQBgiEQAGW2AhSmkgEAkYhgCARhKlmbXLe0dYTsfAAA6AmCIRAAZbZg2OEQ2cEOKACACEMwBAKgX0qi5GckW45RgAIAiDQEQyBABhVY1xlSgAIAiDQEQyBYLWv2EgwBAJGFYAgEaZ0hI4YAgEhDMASC1uS6QRwOR8jOBwCA7iIYAkFaY9jU2i576ppDdj4AAHQXwRAIEK1KTktOsBxjOhkAEEkIhkCAxMXFua0zpAAFABBJCIZAUNcZ0ssQABA5CIZAEIMhU8kAgEhCMAQCqCyfJtcAgMhFMAQCaFCBdcSwqr5F6pvaQnY+AAB0B8EQCKABuWkSHxdnOVbBOkMAQIQgGAIBlJQQLwNyUy3HmE4GAEQKgiEQ5EbXugMKAACRgGAIBJhbL0OCIQAgQhAMgWC3rKHJNQAgQhAMgSCPGO7Y3yht7R0hOx8AAPxFMASCPGLY3uGQ7fsaQ3Y+AAD4i2AIBFhGapLk9Uu2HGPPZABAJCAYAn0wnfzPdzfI5t31ITsfAAD8QTAEgmBUSabl8q6aJvnpf5bK4k3VITsnAAC6QjAEguCcWWVSnGNtdN3Q3C43P7VSXl26PWTnBQCALwRDIAjyM1LkrksPk3EDsyzHOxwOufPNdXL/extDdm4AAHhDMASCJLdfsvzloilyzNj+btc99Uk508oAgLBDMASCKCUpQW4+Y5xcdORgt+ve/WJXSM4JAABvCIZAkMXHxcn3jhkmFx81xHJ8Zfn+kJ0TAACeEAyBPjJnTKFbpbK+AQAQLgiGQB8ZUthPMtMSLcdWMGoIAAgjBEOgD6eUJ5XlWI4xnQwACCcEQ6APTRpEMAQAhC+CIRDCYLhjX6PsrWsO2fkAAOCKYAj0oWH9MyQjlXWGAIDwRDAE+lBCfJxMLMu2HFu5lWAIAAgPBEOgj7HOEAAQrgiGQIiDYUV1g1TXs84QABB6BEOgj40oypD05ATLsZXlNSE7HwAAnAiGQB9LiI+XCfZ1hkwnAwDCAMEQCAHWGQIAwhHBEAiDYLhl7wGpaWgJ2fkAAKAIhkAIjCrOlNQk1hkCAMILwRAIgcSEeBlfmmU5xnQyACDUCIZAiLDOEAAQbgiGQJgEw02766WusTVk5wMAAMEQCJHRJVmSnHjoW9AhIqsqWGcIAAgdgiEQIhoKxw20rjNcV1kbsvMBAIBgCITQsP4ZlsuV+5tCdi4AABAMgRAqyUmzXK7c3xiycwEAgGAIhFBJTqrl8k5GDAEAIUQwBMJoxLD6QIs0tbaH7HwAALGNYAiEUJFtxFAxaggACBWCIRBCui1eXr9kyzHWGQIAQoVgCIRYSa511LByH8EQABADwfAXv/iFxMXFmbfbbrvN6+3mz58vJ598shQUFEhaWpqMGTNGfvWrX0l9fb3Px9+wYYNcdtllUlpaKikpKea9Xt60aVMQvhogWJXJTCUDAKI8GH788cdyxx13mFDoy//+7//KCSecIG+++aaMHz9eTjvtNKmpqZHbb79dpk+fLnv37vV4v4ULF8rkyZPl4YcflpycHDnzzDPNe708adIk+fTTT4P0lQG9Q8saAEBMBcOGhgYzcldSUiLf+ta3vN5u2bJlcv3110tCQoK89tpr8sEHH8jTTz8tGzdulOOOO07Wrl0rV111lcfHP++888z7m266SVavXi1PPvmkea+XDxw4YK5vbOQXLsK/ZQ3BEAAQ1cFQw9n69evlvvvuk+zsbK+3+8Mf/iAOh0O++93vyrx58zqPp6enywMPPCDx8fHy3HPPyZo1ayz3e+ihh2THjh0yatQotylqvazHKyoq5JFHHgnCVwcEdsRQq5L1+wAAgKgLhu+//77cddddcskll5h1g960tLSYUUJ10UUXuV0/ePBgmT17tvn4hRdesFznvHzBBReY8OhKL59//vnm4+effz4AXxEQWCW51mDY3NYh+w60hOx8AACxK6jBUItFvve970lRUZHceeedPm+7bt06MxWsdC2hJ87jOuXsynm5u/cDwkFeRrIkJVi/FSlAAQCEQmIwH/yGG26QzZs3mxG93Nxcn7fV2yktGMnMzPR4m7KyMsttVV1dnVRVVZmPBw0a5PN+e/bsMesN+/Xr5/fXsG3bNp/XV1ZW+v1YgCfxcXFSnJ0qFdUH/zByrjMcX+p92QUAABEVDN9++2259957zfTuGWec0eXtNeApX6EtIyPDvK+trXW7n6/7Ou/nvG93gqEzVALB7mVoCYb0MgQARMtUsraXufzyy6WwsNCsLwTgG70MAQBRO2J47bXXminYp556yjSp9odz+liner1xNrjOyspyu5+v+7o2xna9rz+0mrmrqeSZM2d26zEBO1rWAACiNhjqmsLExES5++67zZsrZ6sZbT+jO5wUFxebnoNDhgwxx/fv32+mhz2tM3SGNOdtld4uLy9Pqqurpby83DS59nY/DandmUZWunsKEGyMGAIAonqNYVtbm2lQ7c2WLVvMm7ahUaNHjzb9CrUyefHixXLssce63UePq2nTplmO62UNmXq97pTi7/2AcG1ZU1XXLC1t7ZKcmBCycwIAxJ6grDHUUT9t0Ovp7dJLLzW3ufXWW81lDYcqOTlZTjnlFPPx448/7vaYW7duNdvqKd3uzpXzso48dnR0WK7Tyzqlrc4666xgfLlAr2lVsittb72rpjlk5wMAiE19tleyP2688Uazl/K///1vs1eyk44iajFLe3u7nH322TJmzBjL/XS7vQEDBpheiL/5zW8s1+llPa5TwtpkGwhH6SmJkp2WZDnGOkMAQEwHQ53qveOOO0wA1F1SdDpZdy0ZMWKEvPvuu2a6+Z///Kfb/XQKWvdU1ve33367TJw4US688ELzXi/rusJnnnlG0tKs03VAuLWscUXLGgBATAdDdd1118k777wjJ554oqxcuVJeeukl04dQ91v+/PPPvVY563Z5K1asMKOCWoiieyrre72sxw8//PA+/1qA7qAABQAQanEOXeiHHtO2PM4m2Fr9TBUzeurB9zfJ4x9v7bw8e1SB/O6ciSE9JwBAbAm7EUMgVtHLEAAQagRDIIynkhnQBwD0JYIhECaKbSOGjS3tUtvYGrLzAQDEHoIhECYKs1IkIT7OcowCFABAXyIYAmEiIT5eirJ737LG2UweAICw2RIPQM8KUHa4hMHuFqDofX/3/GrZvPuAnDCxSH528hi3UUgAiAV765rlnvkbZP+BFrnwyMEyfVheqE8pIjBiCERRL8PHFm6RjbvqpcPhkLdW7pQlm6sDfIYAEP501kT/SP7gq92yony/3PrCF1LT0BLq04oIBEMgilrWrKqosVzWkAgAsWb1thr5antt5+UDzW2ydMu+kJ5TpCAYAlEyYtjc2u62JnHfAf5CBhB7Xl6y3e2Ya1CEdwRDIIwU24LhntomaWvv8Ou+5VUNYi852d9AuxsAsaWqvlk+XLPH7fiX260zKvCMYAiE8VRyh0Nkd22zX/fduveA2zFddA0AseT1ZTukXX942mzYWS8tbe0hOadIQjAEwkhmWpJkpCb2aJ3hlj0egiGLrQHEEJ1heXXZDs/XdThkXWVdn59TpCEYAmE+ariturEXI4ZMJQOIHR+v2ytV9d7/IP6SdYZdIhgCYaY0L91yefEm/1rObN3b4HZM1xhq6xoAiAUveSg6cUUw7BrBEAgzM2xNWDUYaqsFX5o8VCQrDYV17LcMIAZs3l1veha6OmxoruXyV9tr2BmqCwRDIMwcMbJAEl12K2lt75BFG6p83qfCQ0WyE9PJAGLBy0uto4U56Uly5XEjLMd0mtnfgr5YRTAEwrAAZZrtr9wPv9rd7cITp30UoACIcjqrMn/1Lsuxk6cMkKGF/SQzLbHbbWve/3K3XHLPp3Lu3xbK2ysrJZYQDIEwdPSY/pbLn22qlsYW79PJWzwUnjgxYggg2r2zaqc0thxqRaOTLqdOGyBxcXEydkCW342uW9s75B9vr5PbXvzC7D2vmwTc+eY6qY2hJTkEQyAMHTmqQBJcppNb2jrkUx/TyeU+giG7nwCIdq8vt7aoOXJkgfTPOtjhYdzAbL9GDHfXNsnPHl0mLy62Tknrz9/1O2OnzQ3BEAhDWWlJMmVwjuXYRx46+fszlUwvQwDRTHeI2rTb+jPwtMMGdn48bmBWl42ul2yulh8+sNjraGK5h64P0YpgCISpb4y1TSdvrLJMlbhWJO/0sacyU8kAotmSzfssl3WTgCmDD63THjMgy0wte2t0/eaKSrnxiRVS42O6eKuPWZloQzAEwtTsUQUSH3fop1lTa4d8vrHK41+yvpovMGIIIJrpaJ+raUNyLUtx0lMSZUhhP4/9DLdVN8jf31rn9jM03iVIdrVcJ9oQDIEwlZ2e7Dad7GljeF+FJ4pgCCBaaa9W+4jhYUOtvWDV2AH2dYa15r5/fX2tWUNobxl2wyljLMfKq5hKBhAGjh5TaLmsBSjNrdbp5K229YWufykrppKB7tHvMS1QiKVK1Eil6wXtr5O9qbUaV2qvTK6R15dXykpbQ+xTpw6Q350zwa1gRXeRqmmIjf8PBEMgjM0eXWiZ0tD1hJ/btsizr33R9TSuqEoG/FfT0CLfu+8zuebhpfLdfy6SVbbggPCy2DaNXJqXJsU5aW63swc9bXT9z/nrLceKslPlyuOGmyU8JbmpkpRg/SO7vCo2ppMJhkAYy+2XLJMG+Z5Otk8lT7bdvqGl3W2UEYBn2qpkV83BYi4tRvj1M6vMVmsIT0ttwdDTNLIzMNobXeu6bVfXzRstackHb5MQHy8Dc9NjsjKZYAhEWLPrT9fv7Wy1oFXK9orkqUPcp1F0GgRA90egdEeNm55aaXrcIbzoz7/VFdaehNNte807eWp07eqEicVu9x1UYAuGjBgCCAdHjS6QONsI4IK1ez3+oNLbjR2YZdlrWVGAAnStvqlN1u5w72O3t65ZbnpyhdT1YM1he0eH3DN/vXz//s/k3nc3uPXPQ8/p+kBtPeO6vto+Y+JrOtl1T+WrbHsqq0H51mC4lRFDAOEgLyNFJpZZf6A98P4mMz1sLzwpyU2T1KQEyemXbDlOAQrQtRXl+8QlZ7iFgt88s6rbyzJeXbpDnvtsm2zec0CeWVQht77whdl2DYFvUzN+YJZpTeONvdG109VzR0p2epLb8UEF/WKyZQ3BEIgAZ0wvtVzWNVBPf1ru9hfs4K9/kOlfwK4oQAG6ttTW9sRu9bYauf6x5XL3O+vlyU+2yturdpo1bvZ2J67e/WKX5fIn66vk9he/NCOJ6J3FtkK8w7xMI3trdK2OGJkvx9g2E/A2lby7ttnnnvXRwnu0BhA25owpNEUorq0VnvykXMpsUx2Dv/5B5jZiyFQy0O1gqK1LPlm/11SwOq3ZUWveXOX1S5Z/fPewzr15nfYfaPG4xdpHa/fIn19ZI784baxbeyn4R9d82nsLeis8cdLRxCNHFcqCtQcL+DJTE+WaE0eZ9YeelOWlm+U5roPIFVWNMqokU6IZI4ZABNAfXD+eO9Ly125zW4ds2GWtlnR297ePGDKVDHQdNCqqrUFj3pQS+cMFkyU9JcHnfasPtMhDH2x2O/7ZpmqvuxLpSOKdb6w1TZb7ksPhiIqpbHtTaw15I4u7DmxaeXz+4YPk5Ckl8rdLpkmhLcy7SklKkOKc1JgrQCEYAhFiWP8MOXXqoY3hPXFOJWubG1eMGAK+LdviHjRGFGWa77v/OXuipCT6/nWpo1D29YfaQcCXN1ZUyv+9vd6Etb6gRTRXPbhYTvnzh3L7S19GdEBcYptG1m4M/oy+6lrC739zuPzs5DFuawg9GeS2zjD6C1AIhkAEufTooW69uJziXKro7FPJrDEEujeNPGXwoaAxZUiu3P/9mXL1CSPNaJO2NtHdNezdAhZtOLSXuYYue+ub4ycUuTVNfmnJdsv9gknXRm7cVW9GKf/7xS556pNyiUTtHQ5ZsqXarzY1vTXItlwnFgpQCIZABNG/di87epjH67QiWac+VK59Kpk+hoBXOmK31DZiOM22rdqA3DQ5c0apGW365Wlj5U8XTnHby1zDlpP212toto4g6n1/c+YEt5Gt+autBSrBUF3fLAvXWUcwn/h4a0T2Z9ywq07qGtu6tb6wpwbZRwxjYM9kgiEQYU6dWiLD+vfzOo2s3NvVMGIIeLNlzwG3UfVpHhrF231zfJHl8qKNVZ29Dj/dYA1hWrCQn5EiR44qkO8fa/3jTovKgj2drBXUOtLmStcp3//fjRLp08haJKLb2QXDYNuI4bbqxoiegvcHwRCIMLpVk05p2Q0pPPQDLCfdvsawtc8XuQORwj5aqCFDRwj96RbgOjXc2u7orHj91DY9fPiI/M6PjxhZ4Fa8ooEjWDR0vrG80uN17325O+L2g15sm/Y/bFjXIb6nBtla1ujP0R1BfK3CAcEQiECTB+fKseOsvbdmDjv0iyenn3UqWUcKdFcHAF03StbRQm8tTFxlpCbJrOGHvu+c08nbqhpkuy08HD7iUBjU0FmQmWK5fkUQw5k+9vZ93sPMP95e7zaaGK70PNdV1vbJNLLzNc7PsP6hvTXKK5MJhkCE0qq6EycVy/D+GXLV8SNkostWUPYRQ8V0MuBOpwVXltd0exrZ23Ty8q375bXlOyzHNFiMKM7ovKyhU/uSulqx1Xdz7d6wjxbat8zcuLte3rCdc7iq3N8oTa3WqdwxJd73QA5OAUqDRDOCIRCh0pIT5OenjpV7r5gh58wss1yXnBgv/WxbQ1GAArhbs71WmmxtZrT1ib9mjciX9ORDfQ513O25zyost5k5Il/ibSOQ9j19g7XOsLaxVT5cc3B62+m73xgqQ2xFFQ9+sLlHe0H3tU27rb1btWdrrm1EL9AGxVgBCsEQiFL26WRa1gBdry8cXpThVrzli3YCOGp0oeWYfVbWdX2h02RbRbPuruJruren3l29y1IsoaOFJ04qMfsD2wPkwx+5N+kON5t3W6dxtc9ksA2KsZY1BEMgSrkVoBAMATdLbf3wujON7G062VVSQrzHxxyYm+a2dm3F1sCuM9QRyNdtU8SzRxea4KujovZA+/KSHbJ1z4GIGjHsk2BYYB0xrKhqiOpiPoIhEKXoZQj4dqC5Tb7aXuezf6E/pg7JcdttyEl7HaYluzel97TO0HUv9EBYW1knm21BT7eCc7rquOFm2YmThp0nPtkqkRUMu969JNCVyc1tHbKrJvL6P/qLYAhEKXY/AXzTIOY68qOtZyaUWsOavy2kjhlr7RLgaxrZ2zrDFQFeZ2gvginOTrWsnyzOSXNbn/zfL3bLjiBMaQdCQ3ObVO63BrKhfTBimNcv2W3NdjQXoBAMgSjlvl8yI4aAq09sexmPG5htirp64pvj+3stTvHVdsq+l3GgQpmGqPe+2G05dtLkErcimLNnlknq1zsmKQ3KT39aHraNyF3p1zLYNpoXDHFxcW6jhuVR3LKGYAhEKa3Wc8UaQ8DaD+9j2xZxM4f3vB/emAFZUpJj3X1jSGE/MyrnTWlemhmNCsZ08gdf7bZUW2uHGg2GnrbZPHXqAMuxt1ZWmpAa7tPIZflpkpzYsyDfXYPtlcmMGAKI+KnkBoIh4LSqYr/bKPqcMZ5H/fwdVdJqX1ffGFPY5X3c+hkGKBi+/5V1tHDm8Hy3ptpO58wqc9vB5dlF1pY74WBTCCqSvVYmM2IIIPJHDJlKBpycW9c5aaN4f7bB8+W8wweZSt+UxHiz7d35Rwzq8j72tjWBWGfY0tYuqyqsTbtPmFjs9fYaGO2h9tVl26UmzJaf2EcM+2J9oZPbVPLehqDvbx0q7qVSAKJyxFArMFvaOixViHaNLe3S3NrerT5uQKTRdXQL1u512/e4t/R765azJ5jHt6/l88Y+YrinttkUWPQmpGoo1O9112nkrtrwaKh9fXllZzGO7i7ywucVctk3hkk40BC2aU/fVyQ7Dcq3fq66pjbTe9LbKGwkY8QQiFKe2mfU2KaT9ZfA+p118vjHW+Vnjy6TM/76kZzzt4XyP8+vjuo+XYhta3bUuq2hm2Pr6dcb/oZC5xSlfXS/t9PJSzdbm3aPKsmUzDTr57DTIGovoHlxyXbzB2U42F3bLA3N7SGbSi7OSbXscKMWbaiSaEQwBKJURmqiJNj2RHWuqdKRwbvfWS/n//1j+eGDi+XB9zeZRe+6IF/pFlpLNlsb/wLhSv/f6h66e2qbOv8P+/KRbYs4DWeDC/tu9Mm+ztBendzbAhT79+60If4V1VxwxGDL5fqmNnl16XYJx2lk/flW2IejdfFxcTJ9WJ7P5QjRgqlkIErpDzIdidDpDnsvw3+9t1FeWrK9yz1kZwzz3moDCAe6nu7nj6+QL7YdXFOnfwwVZqVIUVaqGeXR/8PfGFtoAphzSvIj2y90+w4gfU2nk7WK2GnFVuuIX3do94ENu6whyh5ovNEq6qNGFcgCl2rtZz/bJmdMLzVb/4VVY+vCfp2vaV+ZM7rQsu/0si37THjWkBpNGDEEYmxbvLrGVnlzRWWX990Sxe0YED20IbMzFCodMdy5v8lMx761cqfc9uIXZqmE08bd9eb6QK8v7I0ptnWGOm26c39jQPZ+1h6FYwdm+X3/C2dbRw31j8nHFm4Nu4rkviw8cZo5It9Svd3W4ZBFG6NvOplgCESxnH7WdUX7Glrl7VU7zZZO9vVFE0qzLce2RvlG8YgO/ky7PvTBZlm8qdrjNLLuBjKiqO9Dhr3i1b7O0B7wXOnooi4FWe7hNvZpZK161v2a/TW6JEsOs20L+NSn5bJhl3XrwL62OQR7JNvp7ieuO8eohVE4nUwwBGJsxPDVpdZtso4eUyiP/PBwufxYa/XhtqoGaWu3BkggHAtJuqKrDm9/6Uuzv609GOpoYV9PSfrTz/D5z7d5LAB7ddkOufWFL8z1v3hihSUY6zS5vfDksC6qkT258rgRlvXJOgp7x2trpL0jND8PtFPCtuqGkFUk+1p28NnGKnN+0YRgCMRQZbKuraqw/YA9/bCBHjv76zRJuO6ZCihdFlFeZf3//ItTx8ivzxgvp0yx7uZR29gqv3xiudvtA1mN3BvHjity2/7tw6/2uH0ND7y3sfOyBsf7/ruhs59eRVWD7LFVWx82tPu7ueho3EVHWqeU1++sD1nT6617G8S1piju6/WQoXDEyALT/sdJ2/ossYXxSEcwBGJoKtm+tkqrMSd/PVKRlZbktj3XFqaTEcbWVlqnN3XK9NjxRXLMuP5y7bxRMntUgeX6bdXWP3TyM5JlTDfW3wXT7NEFbqNgj3y02VJl/eiCLaZ/nqs1O+rk06/bptgDivbYszdm9teFRw5224f44Y+2mJmEUBee6NKXtOTEkP2xPcG27MZTdfK7X+wK2PaGfY1gCMTQVLKd7pHqOo1mb9mx1bZpPRBOvtpu3d1jRHFG53o6/X/981PHykAfjaJ1WrA7PQeDSc/j4qOGWo7p6KazWlkDmbdOArqGUkcP7esLda1gT6fJtVn39aeMMaNzTto0+47X1/R5j9NQ7njiz3TyJxv2Wpbd6PTyn17+Sn75xAq3pQuRgGAIxNCIoSvdtsu+TZZ9hECncIBw9ZVtfeHYAdbRP20j8tuzJ5j/656EyzSy66ihbs3nSkcJddTwvvc2eu3RqJXW7325260x9rQeTCO7GjcwW86cUeq2q4p9nXLfF56EZhrZabbt/01dY1vn6KCueXVuENDa3mE+frmL1mDhhmAIxOiIoU632XdDGGJbZ8hUMsKVrqvTaVRfwdC5Xu66k0e7Hc9OS5KJg6xTguEwanjJnCFuo4Z/f3OdfOzSW/Dgba33/duba03jelddbYPnj+9+Y5ip3HZ1/3sb3XaOCebrvNHWqiYUFcmuirJTZVRxpuXYwnV7TYHMr55aadYdOmmUj7S12gRDIMa2xXM6bdrBohNX9gIUnb4KVSViuGhobpOnPtkqjy/cYprZIjzoL1stxnDlrV/f8ROK5bRpA9z+MEqID79fgUeOKpDhtvY5ry23jtBpa5vrTh5jOWbfLk4fw9f3v7/SkhPcgrUG0Mf7qLdh9YEWt9c51MFQHTXaun5Vp4xvfHKF1NjO9Zvji+QHxw2XSBJ+3xUAAsbeG81pZHGGjC6x/sXraY2hViZvty3YjzW/e3613P/eJnnwg81y81MrOitAEV5tavT/uo7kePPD40eaNbWZaYlmJ5BLj7au5wsXcR5GDe303E+cVOyz/6K9F2FvaGWzfj5Xry/fYdr/9PX6Qm3YrTvahNps23SyBlh7cZ++Bj8/dUzYrGP1F8EQiGK6jZV943fnaKGnRelUJlvpdJlrpeeX22tlZx/8MkTX9LWwjxb6KrTQYopr542WF66bI3+8YLL5vx6ujhxZYP5480SXe5w8pcSEjcu+4T3c9qRNjS+XHT3UbdePxxZukb7f8aRfWAStwQX9TFcHb/T1++1ZE7rVXDxcRN4ZA+iWHFvQS09JcOuZ5qsApTyGC1DWVro3T471EdRwHTH0tL4wUh0cNfQc+q48fnjnFPis4fkep881jNhbqvRWYVaqnDLVOh2vWw4Ge/2cW+FJYeinkV2LhTwpyUmV28+bLOkpkbmHMsEQiHIFGdZgeMKEYrNuyBv7dLI22o1V62198tT2CFtIHo1a2tpl4y5rYBgTRcFQHT4iX0bZlnvMGJYnM4blWwLkdz1MiU8alG1mCwJNexvqyKtTe5BHDXXZxlpbgVGoK5J9ta1xLmnQEelc28/dSEIwBKLc8S4taXT67JxZZT5vby9AieU9k9ft9BAMbTvHoO9t2FlvpjKddGJxdJQFQw1915w4qrPVjgaOq+eOdLud7t3rbFLvNHP4ofAYSPkZKXK6rWjtnVU7g9b0WtsR2XdqGmmrBg6lUcWZlvZCuv7x9+dPkoF5PWsqHi4ic5wTgN/mTS6RvIxks1ZH1y6V5Hhv+OupZY22YNDK5HCs4AwmHa1Y52HEMNJaT8RC/0L9Y6ZfhE7b+aKjoA9dNcv8PxxXmu2xylgD5M9OHi3XP7bcrInV+9grsAPp/CMGyavLtne2ZNF8/p8FW+Smb40L+Od6bZm1GluLTsJlpxrnc/+bs8abBuMt7R1y8VFDwiq49lT0fScBcPvhdfiIAvPmj0G2YNjarnsmN0mZj4XW0Uj3nN3fYG09oQiGofeVrfAknMJCoOnaPn3zRUeoHrv6CNlT22Qqs3u624k/NJx+a3qpPPVJeeex/36xy+ytbF+G0tt9sLVptyvd/zocCk9clealy6/PHC/RJLaGAAB0KTs9yW1kIhbXGXoaLVSV+xu97kCBvhHNhSc9lRAfJ8U5aUENhU7nzRpkWaes3w2PLAjsWsP5q3eZLfhcv76TJpUE9HPAM4IhADdD3LbGIxi6jqDqyAxCY199i1vLIG+NrRG8Px7Psm2Vp3s6r/ewJrenyzh0utrV7FEFEV3QEUkIhgC6nE6OyWDo45cclcnhs75QF/zbC6YQfOfMLHNb13n3O+sD0gB+dUWN2z7t2pwcfYNgCMDNkBhvWeOt8MSJdYbhM42sO/joNCP6lu6zfv7h1g4Hqypq5P2vrOsCe+JVW9HJwLw0mRKAfZ/hH4IhAL8rk2OFbvVl35/VFcEwdL7aXhMzhSfh7uyZZVJs24bwvv9ulKZW677N3VHT0CofrtkT9kUn0YxgCMDvyuRY4WsaWW1j95OQ0KKftbaR3HEUnoSMNtG+8rgRlmN7apstFcvd9c6qSmltP/RHqG7DZ9+nGcFFMATgV2VyLK0z9DWNrHbso8l1KFRUNUhDi3U0ihHD0DpqdIFMHWxtsP3Up+Vm1L1nRSc73HYXyU6n6KQvEQwB+LVnciytM7RXV04ss+47q6OnHQFYZI/uWVWx33K5f1aK2Y0DoaPtcX50wkjLVK+2mbn33Q3dfqwVW/e7jcZTdNL3CIYAPIrVrfE8FZ4cM7a/5bJOdekuE+hbC9Za156NL7UGdoTG0P4Zbrut6DrB5Vv3+f0Yui7x3x9ushwblJ8uk2zb/SH4CIYA/KpMjpVgWLm/Seqa2izHZo3IN21RXO1gnWGf0qKEZVusI4azRxWG7HxgdenRQyUzzdq+5g8vfSkPvL/RVJL7GmFvaWuX3z67Sr7YZq04P2XqgD5p2A0rgiEAv0YMdX1XICqTtaJ3bWVtQPqd9cU0cnZaktlmbGCudY/pbVQm96lP1u+1hIvkxHiZNSIvpOeEQ7LSkuSyo4dZjlXVt8gTH5fLjx9aIhfe9bHc+cZa+WJbjeV7X0fff/f8F7Jks3V0MT8jWU5kp5PoC4aPPfaYXHLJJTJ58mTp37+/JCUlSXZ2tsycOVP+8Ic/SH19vdt9brnlFvMXgq+3NWvWeP2cGzZskMsuu0xKS0slJSXFvNfLmzZZh6gBdC8YamWy/pCvb/LexsXXiMDbq3bKjx9aLJfc86lc/e8lcuOTK3rV1iJY7NPII0syzc8d7aXmigKUvvWRbRp5xrA8SUu2jlAhtE6dWiJDveyXrCFRC0t++shSueqBxfLq0u3mZ8nvX/hCFm2ostw2MzVRbj9/smSk8vqGQlCf9XvuuUc+/vhjGTt2rEybNk3y8vJk165d8sknn8jnn38uDz74oHzwwQcyYID74lINk1OmTPH4uBouPVm4cKHMnTtXGhoaZPz48XLUUUfJ6tWr5eGHH5Znn31W5s+fL4cffnjAv04gWiuTc9KTZH/DoSD40Ieb5clPtsoJE4vlzOmlbm1t7Hbub5RXlu2QN5dXSo2tL6COEPz19bVy0+ljw2q6aF2ldTprVHGmeT/ANmLIVHLfqW9qkyWbqi3Hjh7DNHK4SYiPl5u/NU5++9xqn70+N+6ulzvfXCd/f2u92xRzekqC/PHCyTK8KKMPzhh9HgzvuOMOGTlypAmErqqqquSMM86QBQsWyPXXXy9PPPGE2331eh099JeGwfPOO8+8v+mmm+T222/vvO7mm282I5R6/dq1ayUtzfoDHoBnU4fkyntfWncyaGrtkFeW7jBvWq173IRi80tap5KcFYkfr9srb6zYIUs37xNfE8b//WKXjCrOkHNmDZKwKTzZWe+2s4ayTyX3xbZ4OqKiu3rE+sjYog17pa3j0P+kxPg4mTWiIKTnBO+FKA9dNUvWbK+Vhev3ysdr90pFtefRdXso1HW8fzx/sowuoQVRKAX1p82sWbM8Hs/PzzfB7eijj5a33347IJ/roYcekh07dsioUaPktttus1ynl5977jlZt26dPPLII3LllVcG5HMC0e773xwuu2ubzbogT3QLLH37x1vrTIFGYVaK/PeL3T53DbHTnRKG9c+QaUNDv15MRzkONLe5TSV7HDHc12h+sQVrR4ZHPtosjy7YIkkJ8XL9KWPkm+OLJFbZd8I4bFge04xhTL8nxpVmm7fvHztcyvceMBXlr6+olJ37Pfc3TEmMl9+fN9HcBzFafJKYePCbWtcBBsILL7xg3l9wwQUSH2/9svTy+eefbz5+/vnnA/L5gFjQPytV7rx4qvztkmnyjbH9vYYgHc1ZuG6vvLh4u89QqCOQlx091HJMB4Jue/FLM+0cLNurG0yF5C3PrXbba9fX+kKdSi/MPPgzamCeta9jc1uHVNe3BOV89Rwf+WiLeW708+ii/UZbY+dA+Gxjldzw2DL5/YtfyJ5a7w2J315ZKVfc95n88onlZj1YXxYONba0yedMI0c0XXJy0ewh8sgPD5fbz58kR4zMF9ftrfWPn/85d6JMHsx+yOEgJH9y1dXVdU4Tn3766R5vs3TpUrnxxhulurrarCmcOnWqnHbaaZKZefCvd7tly5aZ99OnT/d4vfO483YA/KPr/7RfnL5peHh56XZ5bVml36OC/VISZe7EYtPnzLkmUYPFwx9t6byNPpauS9IAam8LEwh3vLZGVlbUdO61++CVs8x5dbUV3qivC0+cVZI6qqFBzTVwFnwdHAPpyY+3Wi7rbh9afKHPYyDD56+fXtU5naevwZ8udF/XvXl3vfy/19aa223Ze8CsDZ1Uli1XHDu8T0Z3Fm2oNssTnHRq/YiRTCNHIv3DcubwfPOmO6O8+8Uu2X+gxeyFPNhL0QqiNBjqdPHjjz8uHR0dncUnGg5POukk+dOf/uTxPq+88op5c6UB8e9//7updHalj6XrFtWgQZ7XKpWVlZn3e/bskQMHDki/fv79J9y2bZvP6ysrK/16HCAaFGalyuXHDJeLjxoqn2+skvmrd8mnG6ose5s6aWPak6eUyJzRhWZPVVffPmqIrN9Vb9YiOm3cVW9Gxm48fVxAz1l/8ThDobM6UlufHD+huMsRQ9e1ThoQdTp5s8sOMDqdHOhRjq17DsgCl+fFddQuUMGwubVd/vzKV5Y1Xhr4NOjaR0bfWrnTbS2YPp/XPLJUjhxVIJd/Y1hQf6nbq5GnDM7pXM+KyKUtoC46cnCoTwOhCoZffvmlqQx2ddFFF8lf//pXtwrj4cOHm/WH8+bNk8GDB3fe/49//KO8+uqrcumll0pCQoJ8+9vftgRDJ2+BLyPjUIVTbW2t38HQGSgBiKWH3OzRheZNCyQ+WrNH3v9qt+w70GJGA06aXCKltoBhHzn45Wlj5ScPLZHyqkML0zVonnf4ILPmMFBcg5zrmjV7MGxr73DrYTjy64pkJy1AcX28YPQyfOIT62ih0/Kt+810e3FO74vntLrc9Xl3ff61UbGT9q3UAiFvNNjrHwi/O2eied0DTQOsvZXJnDHWXWgARGAwvPbaa81ba2urlJeXy0svvWQKQt58802zNlCLUJwuvvhit/vPnj3bjB5ec801ctddd8l1110n5557riQns7E2EGoZqUkyb8oA89YdOpWrgeLqhxZLQ/Oh9XMaNAIaDHe790v9fGO1NDS3SbrLdPKijVVmytbV6AHWYDggL7gtayr3N5riHW+0F+Qlc6xrNLtrVfl+eXZRhcfrdGrvkjlDOqfPdaeR6gO+11Fqf8u/vLpG/vOjwz0uA9BR2LdWVpqelWnJCZKenCip5n2CKSDRt8zUJPM+KzVJsvslda5lXbyp2tLrUtelzR7FNDIQNWsMtcG1jgj+7Gc/M2HviCOOkO985zt+t5DRdYl33323mQ5etGiRzJkzxxx3XXeo08SeuDbTzsryvxS+osLzD1DXqWRt2A2g+8ry0+X48cVm3aLT0i375PwjAjfFtMnDiKFOfesUuGulr06ZuhpfmiX5Gdb1gwNyraOgvnq19cTTn5b73DpMg+F3jhpiKQLaVtUgzyyqMMFKR1u1/6Q3WsDyl9fWeG0hpF/Pl9trO/cgnr/a+pyMKMowPSwfW7jVssZUR4pfXrLdfH77tPjPHl3WrUbm+nWMG5glE8pyZGW5dQu8iWU5ktuPAQEgKquStZXNuHHjTPBavHixX/fRfoi6g4p97Z8GQ2evRB2R9BXwCgoK/J5GVrpziq+3khK27AF6Y9pQ6xo9bX+jO6UEc8TQ3gJFg419ytLTdlylHnoZBqpCt6q+Wd5csdNtPZ0rbfWhI36uI4w/fniJvLZ8hzz1abn89JElsreu2evn+Nd7G93CrBZzuHpn1c7OELlgrXWto65xPHtmmRkdnDzIem5PfVJuKoid9Hn521vrur27jTaz/mxjtTz4/iYzYuiKamQgyvdKdga03bu9T524am9vl5qag4vI7dXJurOK8hYyncedtwMQHjRguGYTrUD9Ypv3ljLd0d7hkC0eRgydrVqcQebd1TvNbZ20+ljb89jZp5I19GioDASd3nUt4tEmzlqIM8S2u4xzZFPXRGqbGQ1STtuqG+Xnjy03IdNu6eZqeWnJoZFZ56jot2dbR2c/+Gq3OY+P1+2xTePGybHjijqXAWiPS1e6s422K3KdlraP+PWWrmkFEKXBcO/evbJixQrzsTal9sfLL79sdjbR9S/2tjRnnnmmef/kk0+a6mdXevmpp54yH5911lkB+goABEJmWpJpC+Nq6RbrSFFP6Yiaa3sZVxpAtRWKjmy9aZtGnjOm0GM7G21Noz3XXG0PwDpDnZbVnWRczZ1UbD6fvrePdGqgffCDTbJmh7VYRukuExoOq78Oh1rA8djCLaYdkCsNvz8/daxbpXOdjthtqJJ3VluLTg4bmiu5GYemcccMyJLDR1gLTp5eVG4ahNc1tso/52+wXKdfy7wpJXLsuP6mGbq2vBnWv5/pE+lPiyKd3g5GayAAfbTGUCuJtWfg2WefLampqZbrdAcS3X2kubnZ7F08ceLEzmngDz/8UM455xy3+7z44otyxRVXmI+1Irm42PrD7LLLLpPf//735rF/85vfmI+d9LIe16lfe6sbAKE3bUieJeToVnqXH9P7x93kZRrZ6cM1u00LGvuooqdpZOeo2YDcVNm6t8EynTzRNq3aXS8u3uZWZHHB1+ssj59QJP96b1Pn2kO93Z1vrDMjct5oxfENjy03076PLtwie2rdRxC1D6Gzcly3NtQpfKdnP6tw2+3meA+tcrSCWddqOtU1tskLn28zjb9d99hWPz1plM/+gzpKqff7cnuNrK6okdXbamTTrnqzHlKLVq44dpjX+wKIgGCo08NaWKIBUJtTayhraWkx4U+bV+so3tixYztH8pQ2s9aq5B/+8IfmPgMHDpTGxkYTMtevX29uc+yxx8o999zj9vnS09Pl6aeflrlz55p2Nzq6OGHCBFm9erV502nrZ555hn2SgTA0bUiuPO7S1FkrWXXUSUcTe2Pz7gNuwc61uEOnkzV0uCrOTpXJtrV9rrTPn2sw7G0Biga9FxZb+6UeM66ocwu+vIwUmTk8zxLA7KFQp511NG1nTZMlHP7vG2u97kDzrekDOy9r6x7XYOj6sdIRvSM9hDpt53PUqAJL30Vd69hkq+7WQNhVU2odidXedvrmnLLWafLt+xrMa5KdTtEJENFTyePHjzejdlo5rIUi2m5G+xDqx8cdd5wJdzqi6NqQWnsG/vKXvzRVvlu3bjX3ef31101F8amnnmqaZM+fP9/Sk9CVVjrr9LSOCmrI1P2R9b1e1uM6Ogkg/IwrzTJTm06Or/v2BbqH4TfHW7f1a2rtcKtG1qpbX/sfay9DVzpi2Bs6Nawjba4uOMJa3WufTrb7wTeHy52XTHPbz9lOvyrtMfk/50ywfI1a1JGU4P1r1ql1e4B2usS2xaEWrbiW4+jrevUJI6QntEJZm4wTCoEoGDEsLCyUm2++uVv3yc/PN42se2PEiBFuzbQBhLfkxASZUJZtdt9wXWeogSSQFclTBueanU+WbTn0eey6CmH28KW7hfTGG8t3uK3ls/dxPHxEgWSmJboFSKX7zp45o9Ssvf5/354i1z+6TCr3u+97rBXOVx0/QkYUuW8rqiOz+jnsu4w46XS2N3quWqijRSue6C43gWjKDSAGqpIBwGna0IMtp1z7GfaGjlzZp3mH9u8n3/ARNrVCuqSLEGMPhhVVDbJ8a8/OVe9rn7bVfWM97TRznEvPRSct3LjhlLGdDan7Z6XK//v2VMuoZmlemtx67kT5y0VTPIZCp+Mneg5/uke0BmpfLjlqiBmNtBuUny7nzmL3KCCSEAwBhM06Q3u17y6XNXPdtXXvAcuUphZ0DC7oZ1qe2Fr3dTpxctd9Se1TyVr1rIUef3z5y85KYH+9scI6WpiTniRHeNnZw14Qo1/Dzd8a59bQWtfo3f/9GWbLwdvPnyT/+v5Ms77PGR690S3tPO1BrE3A7b0O7XSvZNdm4U4/OXGUWxU3gPDGdyyAsDC8KMMtmPRm1NA+jawjfVpEoTtnTPJQRaxbtM3xo0+eBq8yD/tA6z7D3733M1Nh7NoT0VcV7tse1jd6C1Ja6HHm9FLzsa4P1NDlrRpap+b1sTTsJfoZzPTzHuOhd+NxPqaRXV08Z4il7Yy2wdEiFwCRhWAIICxo2LEHCW3K3FObbBXJruv25oxxD0C6Ts5bgYUrHXm79byJMqrYfVpWe/j94+31cuebnquBXX2yfq9bS5d5XYxYXj13pDz6o8Pl6WuOlNOmHaoqDpQTJxe7hfXhfu5bra1v/nLRZDl5Solcfsww+dnJowN+fgCCj2AIIGynk7VIxNfewb5s3mMdMRxaeGgHkTmjC9zWxHVVdGIPQXdddpgZtfPUCPuN5ZVdhlq9jasJpdkyyLbLiSdayJETpP2CtQJYdzTRamDdcUWno7uagnY1dmC2/OzkMXLhkYP9HqkEEF74zgUQtvsm64iaty3tfNHdTOwjhkNdRr60N+BZMw9OyyrtxafBrDt03d23DhsoD105y0zb2t377kavU8q6dtK+D7COtIWD8w8fJM9fd5T86wcz3aqjAUS/oLWrAYDu0orgkpxUS7sVHXnrbkCpPtBitplzZX+MH3xzhBwxosAUj+gUdndGxlzpNnE6sjamJFPuevtgI361cXe9vL1qp8fp4bdWVloKY9JTEjxOb4eKrz6OAKIbI4YAwm57vN4WoNh3PNGiiOKcVLcRvylDcs2+vdoOprdOnTZQhrhMV6uHPthk9jV2paOIb6ywTiNrKxp/1jcCQLARDAGE9XTyyvL9poK3N3ska2AL9iiYBs0rjxtuOabNtJ/+tMJybMnmare9i0/20LsQAEKBYAggrOgOHa4RTret27Czrldb4Wlj674wY1i+TB9mHfF8+tNy2Vt3MAjWNLSay65GFGWYVjQAEA5YYwggrOi+uNpz0HUPYl1zqBWvPR0xHFbYd0UUV35zuFkX6aw70TWM/5y/QQbmpcnzn28zO7K4msdoIYAwwoghgLBjXw+4sxs7oLR3dEj53tCMGB78XBkyb7I17L3/1W55bOFWt1CYYra6C5+iEwAgGAIIO7q7iKtdLlXKXdlW1Sit7dY2MUP7cMRQXXb00C6LSXTN4w9PGCkZqe7b0AFAqBAMAYSdYlsw3FlzaFq5K5tsja3zM5Ld9hMONm1hc+ERgzxep+snj59QJA/+YKacOpVpZADhhTWGAMJOUU6aW0PonraqcW1s3ZfOnlkmb63c2blWMu7rbfd0T+HBfuxwAgChQDAEEPYjhrtqms3WeP60nLFvhTesD9cXukpJSpA/XzRFnl1UYVrZ6JZ77CQCINwRDAGE/RpD7WO470CL5GekdH/EsI/XF9q/jqvnjgzZ5weA7mKNIYCwk5eRLEkJcd0uQNFt8OwVzH1ZkQwAkY5gCCDs6JRxf7cClK6D4bpKayNs3epuCOv5AMBvBEMAEVKZ3HUwXFtZa7k8vChDEhP4MQcA/uInJoCo6WW41jZiOLqEreYAoDsIhgCippfh2h3WEcPRJVkBPy8AiGYEQwAR0cuwq6nkvXXNUlXfYjnGiCEAdA/BEEBEjBjurmkyvQz9XV+YnpIgpfnpQTs/AIhGBEMAEdLL0CHVthFBV2t3WNcXjirO9KshNgDgEIIhgDDuZRjv99Z49hFD1hcCQPcRDAGEJR3tK8pO8WudocPhkHU7qUgGgN4iGAKIoJY1niuTK/c3SV1jm+XY6AGMGAJAdxEMAYSt4mz/KpPX2NrU5KQnSf+srvdVBgBYEQwBhK2iHNuIoZdgaF9fOKokU+IoPAGAbiMYAoj4bfHsFckUngBAzxAMAUR0L8P2jg7ZsMsWDAdQeAIAPUEwBBAxU8meehlu3dsgTa0dlmOMGAJAzxAMAYSt3H5d9zJcW2kdLdSiE70fAKD7CIYAIrqX4VpbRTKjhQDQcwRDABHdy3CdbcSQ9YUA0HMEQwAR28uwpa1DNu2ut1zPiCEA9BzBEEDE9jLUUNjWYa1SHlnMiCEA9BTBEEDE9jK0ry8sy0uXjNTEPjs3AIg2BEMAEdvL0F6RzPpCAOgdgiGAiOxlqOHwi201lutYXwgAvcOcC4CI6GXY2n6oifXO/Y3y6YYq2b7PWqHMiCEA9A4jhgAirpehTiE/8N5Gt/WFYwYwYggAvUEwBBBxvQz//cFmqWtqsxz78YkjTYgEAPQcwRBAxPUybGptt1w+dlx/OWxoXh+fFQBEH4IhgIgrQHGVnpwgVx43ok/PBwCiFcEQQMS1rHF16dFDpSDTugYRANAzBEMAERsMh/XvJ2dMH9jn5wMA0YpgCCBip5J/etJoSYjnxxgABAo/UQFERC/DlETrj6t5k0tkfGl2yM4JAKIRwRBA2NM2NMeOL+q8rGsKrzh2WEjPCQCiETufAIgIP5k7UgYXpEtdY5ucNm2AZKcnh/qUACDqEAwBRISUpAQ5d9agUJ8GAEQ1ppIBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgEAwBAABgJB58h55qa2vr/LiysjKk5wIAAGJLcXGxJCYGLs4RDHtpz549nR/PnDkzpOcCAABiS0VFhZSWlgbs8ZhKBgAAgBHncDgcBz9ETzQ1NcmqVavMx4WFhQEdznWdonaORn722WdSUlIS8M8B//BahBdej/DC6xE+eC1i5/UoZio5vKSmpsqMGTP67PPpf6ZADhmj53gtwguvR3jh9QgfvBbhpSTMXw+mkgEAAGAQDAEAAGAQDAEAAGAQDAEAAGAQDAEAAGAQDAEAAGAQDAEAAGDQ4BoAAAAGI4YAAAAwCIYAAAAwCIYAAAAwCIYAAAAwCIYAAAAwCIYAAAAwCIYAAAAwCIYAAAAwCIZh7plnnpFjjjlGcnNzpV+/fjJ58mT585//LK2traE+taiiz+e7774rP//5z2XGjBmSk5MjSUlJUlxcLKeffrq89tprPu8/f/58Ofnkk6WgoEDS0tJkzJgx8qtf/Urq6+v77GuIdr/4xS8kLi7OvN12221eb8drETwtLS3y97//XY466ijJy8uT1NRUKS0tlXnz5slTTz3l8T68HsFRXl4uP/7xj2X06NHmedXXYujQoXLppZfKihUrvN6P16P71q5dK3fddZdcdtllMnHiRElMTOzy51Bvn+8NGzaYz6ffXykpKea9Xt60aZMEne58gvD005/+VHelcSQmJjrmzp3rOOussxw5OTnm2FFHHeVoaGgI9SlGjXfeecc8r/pWXFzsOOWUUxznnXeeY8KECZ3Hf/CDHzg6Ojrc7vvXv/7VXB8XF+c4+uijHeeee655DD02evRox549e0LyNUWThQsXOuLj481zrM/rrbfe6vF2vBbBU1FR4Rg3bpx5LgsKChynnnqq4/zzz3cceeSRjvT0dMfZZ5/tdh9ej+D49NNPHZmZmeZ5HDhwoOP00093nHnmmY6hQ4d2/s54+umn3e7H69G738Vie/P2c6i3z/eCBQvM95Tebvz48eb7TN/r5X79+jk++eQTRzARDMPUCy+8YP4TZGRkOJYsWdJ5XP8jTZw40Vx3/fXXh/Qco8m7775rfrF9+OGHbtc9+eSTjoSEBPOcP/zww5brli5dar7p9frXX3+98/iBAwccxx13nLmPp1+Y8J8+lyNHjjS/AM844wyvP5B5LYJH/wgdM2aMeQ5vueUWR0tLi+V6fY6XLVtmOcbrETyTJk3q/GPV9bVob293/PrXvzbX6SBCY2Nj53W8Hj13//33O2644QbHY4895vjqq68cF198cZfBsKfPt14/YMAAc/1NN91kuU4v6/GysrKgDgwRDMPUjBkzzH+A2267ze26jz76yFyXkpLi2L9/f0jOL9Zcfvnl5jnXb2hX+hegHr/iiivc7rNlyxYzyqXX6w8T9Mw111xjnsPXXnvNcemll3r9gcxrETy/+c1vOoOIv3g9gmPv3r2dI1a7d+92u76trc2RlpZmrtdw4sTrETiX+vg51Nvn+//+7//M8VGjRpmg70ov63G9/p///KcjWFhjGIa2b98un3/+ufn4oosucrte1/eUlZVJc3OzvP766yE4w9gzdepU876iosKy3sq59tDT6zR48GCZPXu2+fiFF17os3ONJu+//75Z23PJJZeYdTre8FoEd/3tPffcYz7WNbj+4PUIHl1v5i9d16Z4PfpWSy+eb+flCy64QOLjrRFNL59//vnm4+effz5o508wDEPLli0z73Vxty4m9mT69OmW2yK41q9fb96XlJR0Hlu3bp00NDRYXg87Xqee08XZ3/ve96SoqEjuvPNOn7fltQiepUuXyt69e2XAgAEyYsQIWbVqlfzud7+TK6+8Um688UbzC7Cjo8NyH16P4MnIyJA5c+aYj3/9619bChH1dbjllluksbHRFATpAILi9ehb63rxfDsvh/J1SgzaI6PHNm/ebN4PGjTI622c3/DO2yJ4du7cKQ899JD5+Oyzz+487nzutYI5MzPT4315nXruhhtuMM+b/gWtVfm+8FoEz8qVK817rYrUIKhdEXQZktOf/vQnM6L+4osvdv7M4vUIrvvvv9+MoN93330mmGtYSEhIMGFBZ5wuvvhi+cc//tF5e16PvrW5h893XV2dVFVV+fz977zfnj175MCBA6ZbSaAxYhiG9D+H8vWC61+Nqra2ts/OKxa1tbXJd77zHampqTFtCnSUxInXKXjefvttuffee810yhlnnNHl7Xktgsf5i0pDh4bAH/3oR6Z9h35PvPPOOzJq1Chz3SmnnNI5esXrEVzaouaTTz6RuXPnmiD40ksvmalFDRk6qqstzrKysjpvz+vRt+p6+Hw77+frvs772e8bSARDwIerrrrK9DfMz8+XZ599VpKTk0N9SlFPA8fll18uhYWFZn0hQss5Oqih78ILLzQjURoGNXgcf/zxJhxqD73Vq1fLk08+GerTjQkLFy40f6jqc/7444+bWY3q6mp55ZVXzOuk3z/6BvQEwTAMOYeedZjYG2dzTNe/ChFYP/3pT+WBBx4w05jOkRFXvE7Bce2118q2bdtMAHEunu8Kr0XwuE6FuY6YO+mUl44WOpv5ut6H1yPw9u/fL2eeeaaZStRRQg3rug5Xf06deuqp8uabb0p6ero8+OCD8t5775n78Hr0rcwePt+u32ve7uvaGDtYrxXBMAwNGTLErQLWznmd87YIrOuvv97s8KBrRHRa01mV7Mr53OsPatcpAFe8Tt2nawp1Z4G7777bTIm5vukvPaWBXS/rVLPitQieYcOGefzY020qKyvNe16P4NE1hRoK9TmfNWuW2/Wux51Bndejbw3p4fOtwVCLTp072/i6n/7RHIz1hYpgGIacIUTX9nhbCLx48WLzftq0aX16brGy9dpf//pXyc7ONqHQW3WYrvPRv8xdXw87Xqeer+384IMP3N527dplrt+yZYu5/Omnn5rLvBbBo8+Xbv+ltDrZE+dx5/onXo/gcQYGX6NF+rNL6fSy4vXoW6N78Xw7L4fydSIYhiGt/tP9epWuH7FbsGCB+atB+1n56u2G7tOqy7/85S/mB6tOHztfB090vaFzCs3T67R161b5+OOPzcc69QP/6F/ZXzffd3vTfWDVrbfeai5rQFS8FsGj+4Vr71TXEShXuqZNQ7qaOXOmec/rETwDBw4079esWWPW43p6PbTFkHK2O+P16FvJvXi+nZd1va69DZRedu5JftZZZwXt/Nn5JMK2xNOu92yJFxy/+tWvOreS+uyzz/y6j742zm2P3njjjc7jbDPV9zsO8FoEz/z5883zl5uba9mntbW11fGTn/zEXKd79+7cubPzOl6P4NDdTnS/XH3+dHeNurq6zuuam5sdV199tbkuKSnJsXHjxs7reD36dueTJT18vl23xLv55pst1+llPV5aWsqWeLG+FZh+g5900knmP5GGFj02e/bsoP7HiDUvvfRS5zZT06dPN9/4nt48hXHXjdKPOeYYx3nnnecoKSlhY/oQ/EDmtQgefc71eUxMTHQceeSRjrPOOssxZMgQc0y3YHv11Vfd7sPrERz/+c9/zOugz2NhYaHj5JNPdnzrW98y+4nrMd1u7Z577nG7H69HzyxZssQxa9aszreCgoLOgOZ6fMeOHQF5vhcsWOBIT083t5swYYLjggsuMO/1sv5R4PrHWTAQDMPcU0895Tj66KMdWVlZ5oev/uf44x//aP4yROD8+9//7gyGvt4GDx7s8f7vvPOOCe95eXlmD+uRI0eaDc9ra2v7/GuJ9b/UeS2C56233nLMmzfPPLf6B2tZWZnjsssu87m/Lq9HcCxfvtw898OGDTPPa3Jysvn59O1vf9uxaNEir/fj9ei+9957z6/fD5s3bw7Y871+/XrHJZdcYkYP9XtN3+vlDRs2OIItTv8J3kQ1AAAAIgXFJwAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAADAIhgAAABD1/wEXV458ursBXwAAAABJRU5ErkJggg==", "text/plain": [ "

" ] }, "metadata": { "image/png": { "height": 270, "width": 323 } }, "output_type": "display_data" } ], "source": [ "volumes = np.linalg.det(cells)\n", "plt.plot(volumes)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", "

\n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "load_atoms.view(final_state.to_atoms()[0], show_bonds=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also express a `(3, 3)` shaped `torch.Tensor` as an anisotropic external pressure.\n", "See the `torch-sim` [docs](https://radical-ai.github.io/torch-sim/index.html) for more details.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Relaxation\n", "\n", "Relaxation in `torch-sim` is provided via `torch_sim.optimize()`. We relax the final state of the NPT simulation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", "

\n", " \n", " \n", "\n" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "relaxed_state = torch_sim.optimize(\n", " system=final_state,\n", " model=carbon_model.torch_sim_model(device=DEVICE),\n", " optimizer=torch_sim.frechet_cell_fire,\n", ")\n", "load_atoms.view(relaxed_state.to_atoms()[0], show_bonds=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Batched simulations\n", "\n", "The real power of `torch-sim` comes when you want to run many simulations at once.\n", "\n", "Here's a simple example of how to run NVE MD for a collection of structures simultaneously. Crucially, you can use this batching process with any of the above examples too!" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using MACE-OFF23 MODEL for MACECalculator with /Users/john/.cache/mace/MACE-OFF23_small.model\n" ] } ], "source": [ "from graph_pes.interfaces import mace_off\n", "\n", "# use the MACE-OFF (small) model to run NVE simulations\n", "# since these structure have no unit cell, it is important to set\n", "# compute_stress=False to avoid an error\n", "mace_off_torch_sim = mace_off(\"small\").torch_sim_model(\n", " device=DEVICE, dtype=torch.float32, compute_stress=False\n", ")" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "563bc72a7fa641c99cf62c48e2b551e2", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "

\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# load the first 100 structures from the QM7 dataset\n", "structures = load_atoms.load_dataset(\"QM7\")[:100]\n", "\n", "# run batched NVE simulations for multiple structures simultaneously\n", "# letting torch-sim handle batching and memory management for us\n", "final_states = torch_sim.integrate(\n", " system=list(structures),\n", " model=mace_off_torch_sim,\n", " n_steps=100,\n", " timestep=0.001,\n", " temperature=300,\n", " integrator=torch_sim.nve,\n", " trajectory_reporter=torch_sim.TrajectoryReporter(\n", " # save each trajectory to a unique file\n", " [f\"QM7-traj-{i}.h5md\" for i in range(len(structures))],\n", " state_frequency=10,\n", " ),\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "torchsim-raw", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.11" } }, "nbformat": 4, "nbformat_minor": 2 }