{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [ "hide-cell" ] }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "%config InlineBackend.figure_format='retina'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Quickstart Guide\n", "\n", "This guide demonstrates the core functionality of `site_analysis` using a simple example: analysing the diffusion of a single mobile ion through a simple cubic lattice." ] }, { "cell_type": "markdown", "metadata": {}, "source": "## Example Scenario\n\nWe'll analyse a short molecular dynamics trajectory where:\n\n1. We have a 4×4×4 simple cubic structure with oxygen atoms at the corners\n2. A single mobile lithium ion moves from one interstitial site to an adjacent site\n3. We track and analyse this movement using `site_analysis`\n\nFor this quickstart, we'll use a pre-generated trajectory (available in the examples directory as `simple_cubic_li.XDATCAR`). This trajectory contains 11 frames showing a lithium ion moving between interstitial sites in a simple cubic structure." }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading the Structure and Trajectory\n", "\n", "First, let's load our structure and trajectory using pymatgen:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "from pymatgen.io.vasp import Xdatcar\nimport numpy as np\n\n# Load the trajectory from the XDATCAR file\nxdatcar = Xdatcar(\"examples/simple_cubic_li.XDATCAR\")\n\n# Get the structures from the trajectory\nstructures = xdatcar.structures\n\n# Get the initial structure\nstructure = structures[0]\n\nprint(f\"Loaded trajectory with {len(structures)} frames\")\nprint(f\"Structure contains {len(structure)} atoms\")\nprint(f\"Composition: {structure.composition}\")" }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's examine the first structure to understand what we're working with:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lattice parameters: (16.0, 16.0, 16.0)\n", "O: 64.0\n", "Li: 1.0\n", "Li ion at index 64, position: [0.375 0.375 0.375]\n" ] } ], "source": [ "# Examine the structure\n", "print(f\"Lattice parameters: {structure.lattice.abc}\")\n", "\n", "# Show the species and their counts\n", "for element, count in structure.composition.items():\n", " print(f\"{element}: {count}\")\n", "\n", "# Find the position of the Li ion\n", "for i, site in enumerate(structure):\n", " if site.species_string == \"Li\":\n", " print(f\"Li ion at index {i}, position: {site.frac_coords}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's also look at how the Li ion moves during the trajectory:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Li positions throughout trajectory:\n", "Frame 1: [0.375 0.375 0.375]\n", "Frame 2: [0.4 0.375 0.375]\n", "Frame 3: [0.425 0.375 0.375]\n", "Frame 4: [0.45 0.375 0.375]\n", "Frame 5: [0.475 0.375 0.375]\n", "Frame 6: [0.5 0.375 0.375]\n", "Frame 7: [0.525 0.375 0.375]\n", "Frame 8: [0.55 0.375 0.375]\n", "Frame 9: [0.575 0.375 0.375]\n", "Frame 10: [0.6 0.375 0.375]\n", "Frame 11: [0.625 0.375 0.375]\n" ] } ], "source": [ "# Track the Li ion position throughout the trajectory\n", "li_positions = []\n", "for struct in structures:\n", " for site in struct:\n", " if site.species_string == \"Li\":\n", " li_positions.append(site.frac_coords)\n", " break\n", " \n", "# Print the Li positions\n", "print(\"Li positions throughout trajectory:\")\n", "for i, pos in enumerate(li_positions):\n", " print(f\"Frame {i+1}: {pos}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining Sites\n", "\n", "Now we'll define the sites we want to track. In this case, we'll create spherical sites at the centers of each unit cell in our simple cubic structure:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Defined 64 interstitial sites\n", "Site 0: center=[0.125 0.125 0.125], radius=2.0, label=Site_0_0_0\n", "Site 1: center=[0.125 0.125 0.375], radius=2.0, label=Site_0_0_1\n", "Site 2: center=[0.125 0.125 0.625], radius=2.0, label=Site_0_0_2\n", "Site 3: center=[0.125 0.125 0.875], radius=2.0, label=Site_0_0_3\n", "Site 4: center=[0.125 0.375 0.125], radius=2.0, label=Site_0_1_0\n" ] } ], "source": [ "from site_analysis import TrajectoryBuilder\n", "\n", "# Define the centers of our interstitial sites (at the centers of unit cells)\n", "# For our 4×4×4 structure, we need sites at the centers of each unit cell (4×4×4 sites)\n", "site_centers = []\n", "for i in range(4):\n", " for j in range(4):\n", " for k in range(4):\n", " site_centers.append(np.array([i + 0.5, j + 0.5, k + 0.5])/4)\n", "\n", "# Set a radius for each site equal to half the unit cell length (2.0 Å)\n", "# This ensures the sites meet exactly at the boundaries between unit cells\n", "# without overlapping\n", "site_radii = [2.0] * len(site_centers)\n", "\n", "# Create labels for the sites (optional but helpful for analysis)\n", "site_labels = [f\"Site_{i}_{j}_{k}\" for i in range(4) for j in range(4) for k in range(4)]\n", "\n", "print(f\"Defined {len(site_centers)} interstitial sites\")\n", "\n", "# Print the first few sites\n", "for i in range(5):\n", " print(f\"Site {i}: center={site_centers[i]}, radius={site_radii[i]}, label={site_labels[i]}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating the Trajectory Analysis\n", "\n", "Now let's use the `site_analysis` package to analyse our trajectory. The `TrajectoryBuilder` class provides a convenient way to set up our analysis:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|████████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:00<00:00, 2866.92 steps/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Analysed trajectory of 11 frames\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# Use the builder pattern to create our trajectory analysis\n", "trajectory = (TrajectoryBuilder()\n", " .with_structure(structure) # Initial structure for reference\n", " .with_mobile_species(\"Li\") # Specify which species we're tracking\n", " .with_spherical_sites( # Define our spherical interstitial sites\n", " centres=site_centers,\n", " radii=site_radii,\n", " labels=site_labels)\n", " .build())\n", "\n", "# Process the entire trajectory\n", "# The progress=True flag shows a progress bar\n", "trajectory.trajectory_from_structures(structures, progress=True)\n", "\n", "print(f\"Analysed trajectory of {len(trajectory)} frames\")" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At the start, Site_1_1_1 (site index 21) is occupied by the Li ion (atom index 64)\n" ] } ], "source": [ "# 1. Query which site is occupied at the start of the simulation\n", "\n", "timestep = 0\n", "first_frame = trajectory.sites_trajectory[timestep]\n", "occupied_sites_start = [i for i, atoms in enumerate(first_frame) if atoms]\n", "if occupied_sites_start:\n", " start_site = trajectory.site_by_index(occupied_sites_start[0])\n", " print(f\"At the start, {start_site.label} (site index {start_site.index}) is occupied by the Li ion (atom index {start_site.trajectory[timestep][0]})\")\n", "else:\n", " print(\"No sites occupied at the start\")" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At the end, Site_2_1_1 (site index 37) is occupied by Li ion (atom index 64)\n" ] } ], "source": [ "# 2. Query which site is occupied at the end of the simulation\n", "\n", "timestep = len(trajectory) - 1\n", "last_frame = trajectory.sites_trajectory[timestep]\n", "occupied_sites_end = [i for i, atoms in enumerate(last_frame) if atoms]\n", "if occupied_sites_end:\n", " end_site = trajectory.site_by_index(occupied_sites_end[0])\n", " print(f\"At the end, {end_site.label} (site index {end_site.index}) is occupied by Li ion (atom index {end_site.trajectory[timestep][0]})\")\n", "else:\n", " print(\"No sites occupied at the end\")" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Li ion trajectory (site indices over time):\n", "[21, 21, 21, 21, 21, 21, 37, 37, 37, 37, 37]\n" ] } ], "source": [ "# 3. Query the mobile ion trajectory (site indices over time)\n", "\n", "li_atom = trajectory.atoms[0] # We only have one Li atom\n", "print(\"Li ion trajectory (site indices over time):\")\n", "print(li_atom.trajectory)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Trajectories for the occupied sites:\n", "Site_1_1_1 trajectory: [[64], [64], [64], [64], [64], [64], [], [], [], [], []]\n", "Site_2_1_1 trajectory: [[], [], [], [], [], [], [64], [64], [64], [64], [64]]\n" ] } ], "source": [ "# 4. Query the trajectories for the two sites involved\n", "# Find the two sites that were ever occupied\n", "\n", "ever_occupied_sites = []\n", "for site in trajectory.sites:\n", " if any(site.trajectory):\n", " ever_occupied_sites.append(site)\n", " \n", "print(\"Trajectories for the occupied sites:\")\n", "for site in ever_occupied_sites:\n", " print(f\"{site.label} trajectory: {site.trajectory}\")\n", "# This shows when the site was occupied (empty list = unoccupied, list with atom indices = occupied)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Site_1_1_1 transitions:\n", " -> Site_2_1_1 x 1.\n", "Site_2_1_1 transitions:\n", " No transitions\n" ] } ], "source": [ "# 5. Query the transitions between sites\n", "\n", "for site in ever_occupied_sites:\n", " print(f\"{site.label} transitions:\")\n", " if site.transitions:\n", " for index, number_of_transitions in site.transitions.items():\n", " print(f\" -> {trajectory.site_by_index(index).label} x {number_of_transitions}.\")\n", " else:\n", " print(f\" No transitions\")\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualising the Results\n", "\n", "For this simple example, a basic visualisation showing the ion movement between the two sites is sufficient:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9UAAAJICAYAAABvxcc2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAewgAAHsIBbtB1PgAAitBJREFUeJzs3Xd8VFX+//H3JJNMQgohlBASCKAgQhIBAQsIKIoNURQrKthW7HV1rbiWdW2r6Fp2XaWp648VxYKFIlUUVEpCDxBQeiihJJAyc39/8M01YWZSJjN3kpnX8/Hg8RjunLnzmTOfJPOZc+45NsMwDAEAAAAAgDqLCHYAAAAAAAA0VhTVAAAAAAD4iKIaAAAAAAAfUVQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPqKoBgAAAADARxTVAAAAAAD4iKIaAAAAAAAfUVQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPqKoBgAAAADARxTVAAAAAAD4yB7sAAAAAAArbdmyRb/88os2b96sgwcPKi4uTq1bt1b37t3VpUsX2Wy2YIcIoBFhpBpAozJ+/HjZbDbZbDYNHDjQ5/PMmTPHPM+oUaN8OsemTZvMc4TiB7BRo0ZVeX3t2rWTy+Xy6VzPPvtslXOFYn/BOk899ZRbPnn716RJE7Vp00b9+/fXo48+qk2bNgU7fASJYRj66KOPdMopp6ht27YaNmyY7r33Xj3xxBO6//77dc0116hr167q1KmTXn31VZWVlXk9lz/+hgAIHRTVAIBa+f333zVnzpw6P84wDE2YMMH/AaFRGThwoFmEWFnYHj58WNu3b9f8+fP1/PPPKzMzU1988UXAn5eiq2HZv3+/zj33XI0YMUKLFy+utu2GDRt0//336/TTT1dBQYHPz1nx/rdv397ncwBoHJj+DQCotQkTJuiss86q02MWLlyo9evXBygihLusrCwNHTrU431HjhzRvn37lJeXpwULFsgwDBUVFenKK6/UwoUL1aNHD4ujRTA4nU5ddNFFmj9/viQpOTlZd911l84991y1a9dOzZo10759+7RkyRJNmDBBU6ZMkST98ssvOvfcc/XTTz8pOjo6mC8BQANHUQ0AqLUpU6bozTffVHx8fK0fM378+MAFhLDXs2dPPfvsszW2W7VqlS688EJt2rRJR44c0SOPPKJvv/3WgggRbO+8845ZUGdkZGj+/Plq27ZtlTZNmjRRWlqaLrroIn355Ze6/PLLVVJSoqVLl2rs2LH685//HIzQATQSTP8GEJYGDhwowzBkGAZFXy307dtXklRUVGSO4tTG4cOHNXnyZElSZmZmQGIDaqNr16569913zf9/99132rNnTxAjglX+/e9/m7dfeeUVt4L6WBdddJHGjBlj/v+9996TYRhV2vA3BEBlFNUAgBpdfvnliomJkaQ6XR89depUHThwQJK4rhRBd+aZZ5p5LIlFy8JAaWmpcnJyzP/X9vKV0aNHKyoqSpK0du1a7dy5MyDxAQgNFNUAwlIwFhFas2aNHn30UfXs2VMtW7ZUdHS02rRpowEDBuiVV16pcdSsffv2VRa9cblc+vjjj3XuuecqLS1NUVFRatmypQYNGqT3339f5eXlfou9adOmGjZsmCRp9uzZ2rx5c60eVzGCExkZqREjRtTpOX/++Wfde++9yszMVHJyshwOh9q2batzzz1Xb7/9toqKijw+7rnnnjPf28suu6zG59m7d6/sdrtsNpvS0tLkdDrd2hiGoRkzZuj6669Xhw4dFBsbqxYtWqh79+569NFHq+2PyitVVyz0tnz5co0ePVqdOnUyzzVgwAC9//77VVZY37Vrl/76178qOztb8fHxio2NVefOnXXLLbcoNze3xtcmHZ32/MADD6hbt25KTExUfHy8TjzxRI0aNUo//fST18dV/hl56qmnJEk7d+7U008/re7du6tZs2aKiYlR+/btdf3112vJkiVu56i8gvzcuXPN4x06dAjKIl6RkZFq1qyZ+X+Hw1Ft+7r2XeUdAc4880zz+IQJE6qsSL5p0yatXr3a/H9iYmK1P6+dOnUy2z733HNe282bN6/G3He5XJoyZYqGDx+udu3ayeFwKCUlRb1799bzzz+vXbt2VdsnFcrKyjRu3DhdcMEFSktLk8PhUFpamvr27at//vOf5pdpnlRetK7C119/rWHDhikjI0PR0dFq3ry5Tj/9dL366qs6cuRIrWLyZPfu3VX+X9udB5o1a6YePXrI4XDI4XC4FdXe/oZU/J6u/DybN2+u8v57G9n2x3vz888/6/bbb1dmZqaaNGkiu91uvpa77rrL488pAD8wAKARGTdunCHJkGQMGDDA5/PMnj3bPM/IkSN9Okd+fr55jup+nZaWlhr33nuvERkZWaX9sf+aNWtmvP/++17Pk5GRYUgyMjIyjP379xsXXnhhtecbNGiQUVRU5NNrMwzDGDlypHmucePGGd9++635/2eeeabGx2/ZssWIiIgwJBkXXnihYRhGrfrr4MGDxogRI6p9bZKMtLQ046uvvnJ7/MaNG802sbGxNfbBxIkTzfb333+/2/27d+82zj///GpjiYmJMd59912P5x8zZozZbvbs2cYrr7xi2O12r+e66qqrDJfLZcyaNcto1aqV13Z2u9348MMPvb6u8vJy489//nONeXfHHXcYR44ccXt85Z+RMWPGGD/88IORkpLi9Tw2m814++23q5yjcg55++fLz1/lPq3L40tKSozo6Ggz3gMHDnhs52vfHfs7wdu//Px8w+VyGccff7x5bOHChR5j2b59e5XHDh482Ovre+KJJ8x2EyZMcLs/Pz/fOOWUU2r8PfTFF19U249LliwxunTpUu150tPTvb6mAQMGmO1KS0uNG2+8sdpzZWdnG7t27ao2Jm8OHz5c5VzvvfeeT+c5lre/IRW/p6v7N27cOLfz1fe9KSsrM2699dZa5d/DDz9suFwuv/QDgKNYqAwAAqi8vFyXXnqpvvrqK/NYenq6zjnnHCUnJ2vz5s2aPn26Dhw4oH379unGG29UQUGBHnroIa/nNAxDo0eP1rRp0xQVFaUzzzxTXbp0UUlJiX788UdzquOsWbM0ZswYvfTSS355LWeffbbatGmjbdu2aeLEiXrssceqHfWZNGmSOeo6cuTIWj3HwYMHNWjQIP3888/msU6dOunMM89UfHy88vLyNGPGDB05ckRbt27V0KFDNWnSJF1zzTVm+w4dOuj000/XwoULdfjwYX377be69NJLvT7n1KlTzdtXX311lfv27Nmjfv36ac2aNZKkpKQkDR48WO3atdP+/fs1e/ZsrV+/XkeOHNEtt9wip9OpW2+91etzjRs3ThMnTpQknX766erVq5dcLpdmzJihtWvXSpI+/vhjnXzyyfrrX/+qQ4cOqV27djr77LPVvHlzbdmyRV988YWKiopUXl6uW265RWeddZZat25d5XlcLpdGjBih//f//p8kKSoqSmeddZa6du2q8vJy/fLLL/rxxx8lSW+++ab27t2rDz/80Ov7WdHXe/bsUUpKis4++2y1bt1aW7Zs0fTp07Vv3z4ZhqE777xTp512mk466SRJR69NTU9Pl3Q0H3777TdJ0h133KGkpCRJsnQF7s8//1ylpaWSpEsuuUQJCQluberTd02bNtVjjz0m6ejo5AcffCDJfYXypk2bymaz6eKLL9Yrr7wiSZo5c6ZOO+00t3h++OEHt/+XlZWZU5MrmzVrliQpIiJCF1xwQZX71q9fr/79+2v79u2SpNatW+ucc85RSkqKdu/erRkzZmjr1q3at2+fLr74Yn311Vdu55Ckn376SYMHD9bBgwclHR2ZPeuss5ScnKytW7dq+vTp2rNnj7Zs2aJBgwbphx9+qPY9HjNmjN5//33ZbDb169dPJ510kgzD0K+//mrOBsjJydHtt9+u//3vf17P401MTIx69eqlX375RZJ0++23mz+vnvqwvu68804VFhZKkjmrIDExUXfddZfZJjs7u8pj/PHePPHEE/rXv/5l/v+kk07SKaecooSEBG3bts18XyTphRde0PHHH6+bb77Z768fCFtBLuoBoE4a20j1008/bd4fERFh/OMf/zDKysqqtNm3b59xww03VBnxmz17ttu5KkZAbDabIcno1auXsWHDhiptnE5nlVG8hIQE4/Dhwz69vmNHqg3DMB566CHz2A8//OD1sS6XyxzJatasmRlDTf1VecSqSZMmxgcffOA2orJt27Yqo/SxsbHG6tWrq7R58803zftHjBjhNc7i4mKjSZMmhiSjU6dOVZ7L5XIZQ4cONc9z6623uo1sOp1OY+LEieY57Ha7sWrVqiptKr8fkoykpCRj+vTpVdqUlZUZw4YNcxtRevLJJ43S0tIqbdevX28kJyebbf75z3+6va5//OMf5v2nn366kZ+f79Zm0aJFRseOHc12//nPf6rcX/lnpCLnHnnkEbdR7R07dhi9evUy295yyy0e+7ry6KSneOrCl5HqWbNmmf2WkJBg5OTkeGznj74zjNr9jpk3b16Nv8/uuece8+ehou2iRYvc2h04cMCcAXHGGWdUua+0tNQ4+eSTzcc/9dRTbu9jaWmp8corr5ij802bNjV2795dpc2+ffuMdu3aGZKMqKgo46233jLKy8urtCkqKqrye6Jjx45GSUlJlTaVc8FmsxkdO3Y0li5dWqWNy+Uy3nvvvSrtfv/9d499VJPPP//c7WerTZs2xl133WV88cUXxp49e+p8ztq8vxX3Z2RkeD2PP96b/fv3GzExMWY/ffLJJ27PU1JSYtx+++3m8xx//PGMVgN+RFENoFFpTEV1QUGB+UFHktfpwYZx9ANk5SK2V69ebh94Kk8rbN26tdcPguXl5UZWVpbZdvHixT69Pk9F9cqVK81jf/rTn7w+dtGiRWa72267zTxeXX+tWLGiyv3fffed1/OXlZUZZ511ltl2+PDhVe4vKCgwC4zExESP05sNo+qH7SeffLLKfZVz5Oqrr672A2jl81x55ZVV7ju2qJ4xY4bHcyxdurRKu6uuusrr891///1muzvvvLPKfYWFhUZiYqL5RYG3Kc6GYRibN282mjZtakgy2rZtW6X4qfz6a3q/f/75Z7Nd165dPbYJVFGdlZVlPPbYYx7//fnPfzZGjhxpdO3a1Wzftm1bY+7cuR7P66++M4za/Y4pKyszmjdvbhaphw4dcmtTUXDddtttZmH94osvurX76quvzOd76aWXqtxX+ffmI4884vU1GYZhvP7662bbhx9+uMp9lfv9X//6V7XnqZyjx14WUDkXYmNjjby8PK/nqfwF2uTJk6t9zuq8+uqr5uUonv6deOKJxi233GL873//M/bu3Vvj+fxVVPvjvVm8eLF5/OSTT/b6eJfLZWRmZpptfZ1SD8AdRTWARqUxFdUvv/yyeV/fvn1rHBXYt2+fER8fbz7ml19+qXJ/5aLa04fqykaPHm22/d///lf3F2d4LqoNwzB69+5tjpZ4GwW/7bbbzMdWHlWrrr/uvPNO877qRpcrrFu3zhxBjYyMNHbs2FHl/sofxr/++muP56g8Q+DY0e5LL73UvO+3336rMZ6KQsHhcFQpjioXIoMGDfL6eJfLZV7vK8n48ccfvbb9+OOPveZv5Q/eH3zwQY1xV45v5syZ5vHKPyORkZHGtm3bvJ7D6XQasbGxhiQjLi7OY5tAFdV1+dekSRPj3//+t1FcXOzxvP7qO8Oo/e+Yyj9n3377bZX7Dhw4YI5OfvLJJ8bZZ59tSDKGDBnidp777rvPPM/atWur3NezZ09DOjpCX9MaA+Xl5UaHDh3MLwsqfm+Vl5cbLVq0ML9wcDqd1Z7nwIED5pcA/fr1q3Jf5Vy4/fbbqz3P3//+d69fFtTVzz//bFxwwQXm7w1v/yIiIowBAwYYH3/8sdfX6a+i2h/vTeUvO5s1a1btyPvMmTONt99+23j77beNffv2Vft8AGqP1b8BIEC+//578/Ytt9xS46qzSUlJVVbsnT17tte2nq51rKxVq1bmbW+rZPuq4vro/fv36/PPP3e7v6SkRB9//LEkqUuXLurdu3etzntsf9WkU6dO6t+/vyTJ6XRq/vz5Ve6vvNr4p59+6vb48vJyffnll5KOXtfbpUuXKvfNnDlT0tFrtGva11aSudpzSUmJ5s2b57HN+eef7/XxNputyvuWlZXltW2LFi283vfdd9+Zt8844wyv7SpUXqV6+vTpHtv06NFDqampXs8RERGhli1bSvJ/vvlTcXGx/vSnP+n44493yxcpMH1Xk4svvti8XZFzFRYtWmSuRn/GGWeY+T5//ny3Veorrqfu0qWLOnfubB4vKCgwV3zu1auXmjRpUm08kZGRGjBggCTp999/N9cTWLp0qbmSdr9+/RQRUf1HyISEBPXp00eStHDhQh06dMhjOyt/l/Xq1UvTpk3T77//rvfee09XX32123oE0tHr6ufOnaurrrpK/fr105YtW+r1vN7467057rjjzJ+/ffv26bTTTtOUKVN0+PBht3MMGjRIo0eP1ujRo811DQDUH0U1AATI8uXLzdv9+vWr1WMqf5Cv/PhjdezYsdrzVC7gDcOo1XPX1lVXXWUu8ONpz+ovv/xS+/btk/THlko1KSsr06pVqyQd/eB46qmn1iqW6vpr6NChiouLk3R0MbJjtyz64YcfzCKh8kJn0tFFpiq2BMrPz6+yHY63fxXbTknSihUrPMbbvHnzal9PZGSkebsi9praHatyP2RkZNQY98CBA2uMu6Z8k2q/VZG/jRw5UsbRmXdu/1wulwoLC7V69Wq999576tq1qyRp27ZtOv/8882cqxCIvqvJOeecY27tdWxRXVH4d+nSRa1atTILqv3791fZUq2goMBcoLDygmiSquzRPHv27FrlcuUtnypeV+W+GTduXK3OU7GFnMvl0urVqz2+/mD8LktLS9ONN96ojz76SNu2bdOGDRs0adIk/elPf1JGRkaVtj/++KP69u1b45aHvvDXe+NwOPTaa6+ZX3SsW7dOw4cPV/PmzXXBBRfob3/7m+bOnWsu0gfA/yiqASBAKn8Ia9OmTa0eU3k0sLoPcbGxsb4HVk/NmzfXRRddJOnoyF7FirUVKj70RURE6Nprr63VOffu3Vvl/DXtH1yhuv6Ki4szV/3evXu3FixYUOX+yqt+X3nllVXuO3Zv27oK1MhWbdQndm9xBzPf6qNiRe4uXbroxhtv1C+//GKOLhcVFbmtsh+IvqtJfHy8zj77bEnSsmXLqsRQkbMVI9R9+vQxfzYqz4aoPMvj2KLaX7kcqJ+JYOeWzWZTx44dde211+pf//qX8vPzNXv2bJ1++ulmm99++00PP/yw35/bn316zTXX6LvvvlNmZqZ57PDhw/rmm2/02GOPaeDAgUpOTtbw4cPNldAB+A9bagFAgNR3VMXfI8z+NHLkSH366adyuVz68MMP9eCDD0qSduzYoW+//VbS0RG4tLS0Wp3PH6/V0zlGjBihSZMmSZKmTJlijiwahmEW1f379692end6enqttwSrcPLJJ9epfaA88sgjNU7TrSw5OTmA0QRfbGys3njjDbPw+Prrr3XgwAElJia6tbWy7y6++GJNmzZN0tEC+YorrlBZWZm5pVRFUR0TE6NTTjlF8+bN09y5c3X33XdL+mPqd8uWLaud5dGtWzddcskldYrN03ZYffv2rTJKXxsdOnSoU/tgqZiBMHfuXF1zzTXmNl7//e9/9frrr9c4RdtX/nhvzj77bOXk5GjlypWaNm2avv32Wy1atMicBl5UVKQpU6ZoypQp+tOf/qS33nqr2lkvAGqPohoAAqR58+batm2bJGn79u21mkK7Y8eOKo9vqM4//3y1bNlSBQUFmjBhgh544AHZbDZ9+OGH5rWeo0aNqvX5Khcke/bsUWlpqaKjo2t8XE39NWjQILVq1Uq7du3SZ599prFjxyoiIkI5OTnatGmTJPep38fG06xZMz377LO1fi3BlpycbObdgw8+GPKFcl1169ZNTZs21f79+2UYhtatW6devXpJCl7fDRkyxLw9a9YsXXHFFVq2bJmKi4slVb3MYcCAAZo3b57mzZsnwzBks9nMonrIkCFuRVLl19CxY0efc7nyeXr27NlofibOO+88rV+/XpL0zTffqFOnTrV6nN1u1xtvvGEW1cXFxdq0aZN5CYE/+Ou9qcxmsykzM1OZmZl6+OGHVVZWpqVLl+q7777Te++9p82bN0uS/v3vf+uEE07Q/fffX+/nBMD0bwAImMoLTS1cuLBWj6ncrrqFqoItKirKXAhsxYoVWrp0qQzDMK+xTkxMrLIAU02io6N1wgknSDq66NjixYtr9bia+stut+uqq66SJG3dulU///yzJOmzzz4z7x8+fLjb49q3b29OS123bp1KSkpqjOXIkSMqLCxUYWFhUK9drPyhv/J1t96Ul5ebcVcUcaGu8uUFla+1D1bfpaam6pRTTpH0x3XVFddTt2/fXu3atTPbVoxa7969W6tXr9amTZu0ceNGSe5TvyX311SbWSHFxcXm66ron7r2jSQdOnTIPI/L5arVY/ytpKREGzZs0IYNG5SXl1enx6akpJgLgEky11nwF3+9N9WJiopSnz599MQTT2jt2rW6/vrrzfveeecd3wIH4IaiGgACpPLKwP/5z39qbH/gwAFzVERSnadXWq3ylOgJEyZo2bJl5oftq666qs7XSta1vyqufZSOjs54W6258irgU6ZMkfTH9dTnnnuuxxHuqKgocxptSUmJvvjiixrjufbaa9WsWTM1a9bMXJU3GCr3w+TJk2ts/84775hx16bfG7s9e/Zo165d5v8rF6zB7LuKL6E2btyo/Px8t+upK5x22mmy249ONJw3b545Su1wOHTOOee4nTctLc2cer1p06ZaXU/bv39/NWvWTMnJyTp48KCko5c0VPxMz58/320thWOVl5fruOOOU7NmzXTcccfV+JyBkp2dbd6u+PmvLafTqf3795v/T09P91tckv/em2effdYcna7uiwOHw6EXXnjB/H/FqDWA+qOoBoAAGTVqlDmFee7cuR5Xyq5gGIbuv/9+cySke/fu5nY0DVX37t3ND6wfffSR3n33XfO+ul6DLEm33nqreXvixIlVFl86Vnl5uUaPHm2Ofl188cVeF4Pr3bu3jj/+eElHt9bKz883VzK++uqrvT7HjTfeaN5+5JFHqh2l2rx5s7766itJR4u0yosFWe366683V0x+9913tWzZMq9tS0tLqxSDNW1v5C/BXC+g8uhct27dquRNoPquNq+38ijzzJkzvRbVcXFx5nT1uXPnmkX12Wef7XXF+Mq5fN9991U7k2Lx4sX69ddfJR3dtaBZs2aSjl7PXXGphNPp1L333lvt6/r000/NLy8uvPDCOl2f7k+VFyGcMGGCeZ16bXz11VdmX3Xo0KHWC04eq7p+8sd7ExERoZUrV2rlypXm7yFvKm/FVt3WfADqhqIaAAIkJSWlyurCN910k1577TW3KXv79+/XLbfcovfee8889uKLLwZti6K6qCied+/ebRYrnTp10mmnnVbnc3Xv3t1cLdwwDA0dOlQfffSR2wfS7du369JLLzX3BXY4HHrmmWe8ntdms5mj1Rs2bNDTTz8t6eiiVdVNUb/88svVrVs383GDBg3SunXr3NqtX79eQ4YMMaeI33///UErIKSj04VvvvlmSUe3Khs8eHCV/ZcrFBYWauTIkeYXDJdccon55UMgVGzDJh3dAspqhmHogw8+qLL12bHXk/qz7+r6ert27WqO6L711lvmY44tqiWZW2vNmTPHLKqry+XbbrtNKSkpko5uJTd06FDz2vHKlixZomHDhpn/f+CBB6rc//DDD5uj1ZMnT9bIkSNVWFjodp6ZM2fqpptuknT05+++++7zGlugnXbaaeYXHk6nU+eff75mzJhR4+Nyc3OrfNF377331vnnumJGwb59+1RWVuaxjT/em8pfyDz++OOaMmWKx+n2Bw8e1F133WX+v/LsIAD1w0JlABqtnJwcnXfeeXV6zAMPPOBximSgPPnkk1q0aJFmzJghp9Op++67T6+88orOOeccNW/eXJs3b9Z3331XZRT0r3/9q6Ux1seIESP00EMPyel0msXvyJEjff5C4J///Kdyc3O1fPlyFRUVacSIEfrrX/+qgQMHKiEhQXl5eZo+fbqOHDki6egH9nfeeafGkeGK80h/bPk1dOhQxcfHe32Mw+HQhx9+qIEDB6qwsFC//PKLunbtqv79+5vPt27dOs2cOdMc/RkwYIDuuOMOn167P73yyitavHixli9froKCAp133nnKzs7Wqaeeqvj4eG3evFkzZ840p7ampKTojTfeCGhMlUf5rr32Wg0ZMkSxsbHq2bOnLrvsMp/Pu2TJEj3++ONe73e5XNq3b5/mzZtXZV/qiy66yONiev7qu8qv99tvv9U111yjdu3ayW6364EHHjBHGSvYbDZdfPHF+sc//mGOkLdu3drjFx0DBgzQCy+8UGWhvsqLnR2refPmmjRpkoYMGaLS0lJ999136tixowYNGqROnTqptLRUq1at0ty5c83HjBgxwu0a7U6dOuntt9/WDTfcIMMwNGnSJH322Wc655xz1KFDBx08eFDLli0z1y6QpMcee8zjCuJWsdls+ve//61+/fpp06ZNKiws1ODBgzV48GBdc801OuWUU5SWliaXy6WdO3dq3bp1+vzzzzVp0iTzi7LzzjvPp5/rNm3a6LffftPBgwd1zjnnqE+fPoqOjtall16qnj17SvLPe5OZmambbrpJ7733noqLizV8+HC1a9dO/fv3V0pKisrLy7VhwwbNmTNHhw4dknR0K7cnnniiPl0LoDIDABqRcePGGZJ8/jdu3DjDMAxj9uzZ5rGRI0f6FEt+fn6Vc3tz5MgRY/To0YbNZqs2toSEBOPtt9/2ep6MjIwan6vCmDFj3F5zXY0cObJW57jwwgvNdjabzdi8ebPXtrXpr8LCQuOyyy6r8b1s1aqVMWXKlFq/nt69e1d5/Oeff16rx+Xm5hrdunWrMZ5hw4YZBw4ccHt8Xd6L2r7Htcnfffv2GUOHDq0x7szMTGPNmjU+PUddYp8yZYrH5/fl569yn/ry74YbbjAOHz7s9fz17TvDMAyXy2VkZWV5fFx+fr7Hx8yZM6dKuyuuuMJju8LCQiMiIsJs16dPn1r125w5c4x27drV+LpGjx5tlJaWej3P5MmTjebNm1d7jsjISOOvf/2r4XK53B4/YMCAGvuiQuXf+WPGjKnV6/Tk999/N/r27VvnXBkxYoRRXFzsdr7a/HzcddddHs/p6fdAfd+bkpISY9SoUbV6TWlpacaCBQt87ksA7pj+DQAB5nA49Pbbb2vp0qW6//77lZ2dreTkZNntdqWkpKhv37567rnnlJeXp9GjRwc73DqrPNp31llnVVn4yRdNmzbVJ598ovnz5+u2227TiSeeqKZNmyoqKkpt2rTRWWedpVdffVV5eXm69NJLa33eiqnlkpSUlKRzzz23Vo/LzMzUsmXLNHHiRA0bNkzt2rVTTEyMHA6HOnTooBEjRmjmzJmaMmWKEhIS6vx6AyUpKUmff/65vv/+e918883q3LmzEhMTzby74IILNGHCBP3666/myuuBNGzYMI0dO1aZmZl1XsSuvmJjY9WxY0fdeOON+vHHH/X+++8rJibGa3t/9J3NZtNnn32mSy+9VK1atarV1OG+fftW2WbJ09Rv6ejPSOXRX0+rfnsyYMAArV69Wm+++abOPfdcpaWlKTo6Wk2aNNEJJ5ygW265RT///LPefvvtKtPXj3X55ZcrLy9PL774ogYMGKDU1FRFRUUpPj5eWVlZuvfee7VixQo9+eSTDeYylvT0dM2bN08fffSROX3em+joaF144YWaPXu2PvjgA5/z9fnnn9fdd9+tjIyMavtTqv97Ex0drXHjxmnRokW644471L17dyUmJioiIsLM/0suuUTvvvuu8vLy1LdvX59eEwDPbIYRxNVCAAAAAIvt2bNHv/76qzZs2GBeF56cnKwTTzxRPXr0aFBfkAFo+CiqAQAAAADwEdO/AQAAAADwEUU1AAAAAAA+oqgGAAAAAMBHFNUAAAAAAPiIohoAAAAAAB9RVAMAAAAA4COKagAAAAAAfERRDQAAAACAjyiqAQAAAADwEUU1AAAAAAA+oqgGAAAAAMBHFNUAAAAAAPiIohoAAAAAAB/Zgx0AanbkyBHl5uZKklq2bCm7nbcNAAAAAOqqvLxcBQUFkqSsrCzFxMTU+5xUZ41Abm6u+vTpE+wwAAAAACBkLF68WL179673eZj+DQAAAACAjxipbgRatmxp3l68eLFSU1ODGE1VJSUlWrlypSSpW7ducjgcQY4IoYpcg5XIN1iFXIOVyDdYpSHn2vbt281ZwJXrrPqgqG4EKl9DnZqaqvT09CBGU1VJSYl5TUJ6enqD+oFBaCHXYCXyDVYh12Al8g1WaSy55q+1qpj+DQAAAACAjyiqAQAAAADwEUU1AAAAAAA+oqgGAAAAAMBHLFSGeomIiFB8fLx5GwgUcg1WIt8QSKXlLs1YtVNz1+1Szpb92lhwSGVOQ9Fffa+OLeOVlZaoAZ1b6ZyuKYq2k38VKvdb7tYD2lhwSKVOl6IjI+i3apBvdUeu+Sacc81mGIYR7CBQvS1btqht27aSpN9//71Brf4NAABqp8zp0nsL8vWf+fnafaikxvYtExy6qV8H3dSvg6IiQ+sDaF3Qb76h3+qOPvNNY+u3QNRWYVlUFxcX65133tH//vc/rV69WkVFRWrevLlOPvlkXX311br66qsVGRlptp8zZ47OPPPMOj3HmDFj9NRTT/klXopqAAAat3U7D+r+ycu0YuuBOj82My1R/7iiuzqnJAQgsoaNfvMN/VZ39JlvGmO/BaK2CruvVH7//Xf17t1bDzzwgH766Sft379f5eXl2rlzp77++mtdd911GjRokA4ePBjsUAEAQAj4dfNeXfbWQp8+dErSiq0HdNlbC/Xr5r1+jqxho998Q7/VHX3mG/rtD2FVVLtcLl1xxRVatWqVJOnuu+/WkiVLtGHDBn3zzTc677zzJElz587Vn//8Z/Nxp59+urZv317jv3HjxkmS2rZtqz/96U/Wv8AgKCsr07p167Ru3TqVlZUFOxyEMHINViLf4C/rdh7UqPd/1sGS8nqd52BJuUa9/7PydobHl/70m2/ot7qjz3xDv1UVVkX1d999p59++kmS9Nprr2ns2LHq0aOHOnbsqPPOO0/Tpk3TBRdcIEkaP368OVodHR2t1q1bV/tvz549uvvuuxUXF6evv/5abdq0CdrrtJLL5dLevXu1d+9euVyuYIeDEEauwUrkG/yhzOnS/ZOX1ftDZ4WDJeW6b/IylTlDOyfpN9/Qb3VHn/mGfnMXVkX1tGnTJElt2rTR7bff7nZ/RESErrnmGklSSUmJ8vLyanXeoqIiXXbZZTp48KBeffVVZWZm+i9oAADQKL23IN/naZHerNh6QO8tyPfrORsa+s039Fvd0We+od/chVVRvXHjRklSjx49FBUV5bFNkyZNzNt2e+12HLvnnnu0du1anX/++br55pvrHygAAGjUSstdAfuA+N6C/EY9olMd+s039Fvd0We+od88C6uiesqUKTp48KA++eQTj/eXl5frzTfflCSlpaWpS5cuNZ5z5syZeu+992S32/Xaa6/JZrP5NWYAAND4zFi1UwUHa95axhcFB0s0feXOgJw72Og339BvdUef+YZ+86x2Q7EhIjY21u2YYRjav3+/li9frueee06zZs2S3W7Xv/71L0VHR1d7vpKSEnMa+Z133qnOnTsHJG4AANC4zF23K6Dnv/OjJbr749D7It/lCuxOr/Sbb0Kx3+gz3wS63+atK9CF2akBfY5ACKui2pPY2FiVlPzxbUuHDh00adIk9e3bt8bH/utf/1JeXp4cDocefvhhn2PYsmVLtfdv377dvF1SUlIl3sqioqIUEXF08oHL5apxxVqHw2HedjqdKi/3vthARERElSnzZWVlcrlcKi0tNRfxKS0tNe+PjIysMn2+tLRU1W2Jbrfbzb3BDcOoci5PoqOjzVkBNcVus9mqfEFSXl4up9Pptb2311qb2CV5fX8qBON98qYxvU/H5hrvU8N8n47VWN+ninyrOHdtYud98i4cf55y/Xy94bEMSc4Af7gNRfSbb+i3uqPPfJOzZV+VvymB+PtU098sX4R9UX2s/Px8Pfroo3r11VfVs2dPr+0OHz6s559/XpJ04403qnXr1j4/Z8Xm47WxcuVKFRQUeLyvW7duSkg4unl6UVGRVq5cWe25Tj31VPP2rl27tHnzZq9t4+PjqyzAlp+fb66Ku3//fklSbm6u+aEpNTVVGRkZZvsVK1ZU+0GkU6dOat68uaSjH3CWLl1abew9evQwP3QVFhZWu6hcdHR0lfdy69atVb6oOFZycnKVWQdr167VoUOHvLbPyMhQauof36jVFHsw3idvGtP7dGyutWjRgvfJC36evKvt++RyuVRUVGSeW+J9qqyhvE9Sw/152ljgvf8AAJ5tKDhU5fdsIP4+7drl/5lEYXVNtSe5ubnKz8/XnDlz9Pjjj8vhcGjevHk6++yztWbNGq+PmzBhgnbs2CHp6H7X4cpmsykmJkYxMTFcT46AItdgJZvNpri4OKWmplb5xhuordJGutgOAARTmfdJVQ2azahuzlMY+uqrr3TRRRdJkq699lpNmjTJrY1hGOrRo4eWL1+u0047TQsXLqzXc9Zm+nefPn0kSevXr1d6errHduE4va42sTMNkvfp2Ngl3qfKeJ/8Ezvvk3fh+D6d8Pg3KimnsAaAuoiJitDyx88y/x+Iv09btmzR8ccfL0n6/fffvdZWdcH072MMGTJEl19+uf73v//ps88+k2EYbqNiv/76q5YvXy5Juv766+v9nHV5Ix0OR5UPG95ERETUql2FyMjIOo3GeNuSzJuaFn2rzGazBTR2u91e6+3SpLq/1rrEzvvkHe+Td7xPnvE+ecf75F2g3qeOLeO1envgrqtOjovW0JPaBOz8wfL5sq3aV1z9ly71Qb/5JhT7jT7zTaD7rWOL+Gp/z/rj71Nd/gbVVtgU1Rs3btTrr78uSXrsscfUsmVLr2379++v//3vfyoqKtLu3bvd2k6ZMsW8PXTo0MAEDAAAGq2stMSAFtXnnJiip4Z2C9j5g6W4tFyTf6l+Bl990G++CcV+o898E+h+y0prGrBzB1LYXFNdXFyssWPHauzYsVq7dm21bStPkfP0DfbUqVMlST179lSbNqH3DVRdlJaWasmSJVqyZEmNU+KA+iDXYCXyDfU1oHOrgJ6/f2fvgwONGf3mG/qt7ugz39BvnoVNUX388ccrJiZGkvTtt99W2/a7776TdHSF0MTExCr37dq1y1zAbMCAAQGItHGpuL6spmvSgPoi12Al8g31dU7XFLVM8P8UQ0lqmeDQ4G4pATl3sNFvvqHf6o4+8w395lnYFNUxMTEaMmSIJGns2LFeV/b+8MMP9c0330iSrrnmGrfrqZcsWWLe7tGjR4CiBQAAjVm0PUI39esQkHPf1K+DoiJD8yMc/eYb+q3u6DPf0G+eNc6offTMM88oLi5Ohw4dUr9+/fTqq69qxYoV+u233zRnzhxdfvnluvbaayVJaWlpevTRR93OkZuba97Ozs62LHYAANC43NSvgzLTEmtuWAdZaU11c4A+0DYU9Jtv6Le6o898Q7+5C6uiukuXLvr000/VtGlT7dmzR/fff7+ysrKUkZGhM888U5988okk6YQTTtCMGTOUnJzsdo49e/aYt1u1Cuw1BQAAoPGKiozQP67orgSHf9aFTYix6x9XnCR7Ix3JqS36zTf0W93RZ76h39w13sh9NHjwYK1Zs0aPP/64Tj75ZCUmJioyMlKtWrXS4MGD9e6772rZsmU68cQTPT6+sLDQvN2sWTOLogYAAI1R55QEjb+xd70/fCbE2DX+ht7qlJLgp8gaNvrNN/Rb3dFnvqHfqrIZrMDS4G3ZskVt27aV5L8Nyv2lpKRES5culXT0GvNA7PsGSOQarEW+wd/ydh7UfZOXacXWum+zlZmWqFev6N7oP3T6gn7zDf1Wd/SZbxpjvwWitgq7kWoAAACrdUpJ0Ge399WDgzvX+jEtExz6y/ld9NntfcPyw7r0R7/95fwutV5xmH6j33xBn/mGfjuKkepGoCGPVDudTnNKfFJSkiIjI4MbEEIWuQYrkW8IlKW/7dOwtxZ6vC8mKkIdW8QrK62p+nduqcHdUhrtSriBUOZ0afrKnZq3rkC5W/dr4+5DKil3yWGn36pTtd8KtaGgSKXlLjnIN6/INd80llwLRG1FUd0INOSiGgAA1N7EHzfpyc9Xuh3/6JZTdPpxLYIQEQCEF6Z/AwAANGI5W/Z7PJ6Z1tTiSAAA/uKfddARtgzDUGlpqSQpOjpaNpstyBEhVJFrsBL5hkDJ9VBUt2/eRA6bS4ZhkGsIKH63wSrhlmuMVKNeSktLtXTpUi1dutT8wQECgVyDlcg3BEJxabnydh10O54WW06uwRL8boNVwi3XKKoBAAAssGrbAbk8rGTTMYmJgwDQmFFUAwAAWMDb9dQU1QDQuFFUAwAAWCB3q3tRbbNJHSiqAaBRo6gGAACwQM6WQrdjx7WIU4w9tBfwAYBQR1ENAAAQYAePlGnj7iK345ltEoMQDQDAnyiqAQAAAmzltgMyPCxSltkmwfpgAAB+RVENAAAQYJ72p5akzDRGqgGgsWNlDNRLdHS0evToYd4GAoVcg5XIN/hbjodFyiIjbMpu21yRbZMkkWsIPH63wSrhlmsU1agXm80mh8MR7DAQBsg1WIl8g7/lelikrFOreDVx2MXHMViF322wSrjlGtO/AQAAAmh/cZk27Sl2O56V1jQI0QAA/I2vRlEvTqdThYWFkqSkpCRFRkYGNyCELHINViLf4E8rtnm+njo7vSm5BkuRb7BKuOUaI9Wol/LycuXl5SkvL0/l5eXBDgchjFyDlcg3+FOOl0XKstKTyDVYinyDVcIt1yiqAQAAAih3a6HbMXuETV1as50WAIQCimoAAIAA8jRSfULrBMVEhfZ0SAAIFxTVAAAAAbK3qFRb9h12O56dziJlABAqKKoBAAACJNfD/tSSlJWWZG0gAICAoagGAAAIEE/7U0uMVANAKKGoBgAACBBP11NHR0aocwqLlAFAqGCfatSLzWZTdHS0eRsIFHINViLf4C8rPEz/PjE1QdH2o+Ma5BqsRL7BKuGWaxTVqJfo6Gj17Nkz2GEgDJBrsBL5Bn8oOFiibfuPuB3PqjT1m1yDlcg3WCXcco3p3wAAAAHgaZRakrJZpAwAQgpFNQAAQAB4up5aqjpSDQBo/Jj+jXopLy/X1q1bJUlpaWmy20kpBAa5BiuRb/CH3K2Fbscc9gh1ahVv/p9cg5XIN1gl3HKNkWrUi9Pp1Pbt27V9+3Y5nc5gh4MQRq7BSuQb/MHTSHW3NomyR/7x8Ytcg5XIN1gl3HKNohoAAMDPdh44ol0HS9yOZ6cnWR8MACCgKKoBAAD8zOv11GlcTw0AoYaiGgAAwM9ytxR6PJ7NImUAEHIoqgEAAPwsx8N2Wk2iI9WxZbyH1gCAxoyiGgAAwI8Mw1Cuh+nfmW2aKjLCFoSIAACBRFENAADgR9v2H9GeolK34+xPDQChKbQ3DEPARUREKDk52bwNBAq5BiuRb6iPulxPTa7BSuQbrBJuuUZRjXqJiopS586dgx0GwgC5BiuRb6iPuqz8Ta7BSuQbrBJuuRb6XxsAAABYKNfDImUJDrvaN48LQjQAgECjqAYAAPATwzA8jlRnpjVVBIuUAUBIYvo36qWsrExr166VJJ1wwgmKiooKckQIVeQarES+wVe/7z2s/YfL3I5725+aXIOVyDdYJdxyjaIa9eJyuXTo0CHzNhAo5BqsRL7BVzlbCz0e97byN7kGK5FvsEq45RrTvwEAAPzE0/7UkpSdlmRtIAAAy1BUAwAA+Imn66mbxkapbXJsEKIBAFiBohoAAMAPXC5DKzys/J2d3lQ2G4uUAUCooqgGAADwg017inSwpNztuKf9qQEAoYOiGgAAwA887U8teV/5GwAQGiiqAQAA/MDT9dSSlJWeZG0gAABLsaUW6sVutysjI8O8DQQKuQYrkW/whaeVv5vHRatN0xivjyHXYCXyDVYJt1wL/VeIgIqMjFRqamqww0AYINdgJfINdeV0GVqxzb2ozqphkTJyDVYi32CVcMs1pn8DAADU08aCQyoudbodz2aRMgAIeRTVAAAA9cT11AAQvpj+jXopKSnR0qVLJUk9evSQw+EIckQIVeQarES+oa58XfmbXIOVyDdYJdxyjZFqAACAesrZUuh2rFWCQymJ3hcpAwCEBopqAACAeih3urRy2wG34+xPDQDhgaIaAACgHvJ2HVJJucvteFZakvXBAAAsR1ENAABQD572p5YYqQaAcEFRDQAAUA85Wws9Hs9kOy0ACAsU1QAAAPXgaaS6TdMYtUwI7dVuAQBHUVQDAAD4qLTcpdXbD7odz2LqNwCEDfapRr1ERUWpW7du5m0gUMg1WIl8Q22t23lQpU73Rcqy05Nq9XhyDVYi32CVcMs1imrUS0REhBISEoIdBsIAuQYrkW+orRwvi5Rl1fJ6anINViLfYJVwyzWmfwMAAPgo18siZbUtqgEAjR8j1agXl8uloqIiSVJcXJwiIvieBoFBrsFK5Btqy9NIddvkWDWLi67V48k1WIl8g1XCLddC+9Uh4MrKyrRy5UqtXLlSZWVlwQ4HIYxcg5XIN9TGkTKn1u5wX6QsOy2p1ucg12Al8g1WCbdco6gGAADwwZodB1XuMtyOs/I3AIQXimoAAAAf5G4p9Hg8m+upASCsUFQDAAD4wNvK390oqgEgrFBUAwAA+CB3q3tR3aFFnJrGhv6erACAP1BUAwAA1NHhUqfW7XRfpIyttAAg/FBUAwAA1NGq7fvlYY0yZbNIGQCEHYpqAACAOvJ2PTUj1QAQfuzBDgCNm8Ph0KmnnhrsMBAGyDVYiXxDTXI9FNU2W90XKSPXYCXyDVYJt1xjpBoAAKCOcjwsUnZcy3jFOxivAIBwQ1ENAABQB4dKyrWh4JDbcfanBoDwxNepqBen06ldu3ZJklq1aqXIyMggR4RQRa7BSuQbqrNy634ZHhYpy/JhkTJyDVYi32CVcMs1imrUS3l5uTZv3ixJSk5ODvkfGAQPuQYrkW+ojqf9qSXfVv4m12Al8g1WCbdcY/o3AABAHXha+TvCJnVNZfo3AIQjimoAAIA68DRS3TklQbHRoT0SAwDwjKIaAACglvYfLlP+7iK34+xPDQDhi6IaAACgllb68XpqAEBooKgGAACoJW+LlGWlJ1kbCACgwaCoBgAAqKUcD0W1PcKmLq0TghANAKAhYEst1EtERITi4+PN20CgkGuwEvkGb3I9rPzdOSVBMVG+LVJGrsFK5BusEm65RlGNeomKilJmZmaww0AYINdgJfINnhQWl+q3vcVux+tzPTW5BiuRb7BKuOVa6H9tAAAA4Afer6dmkTIACGcU1QAAALWQ42HqtyRlpyVZGwgAoEFh+jfqpaysTPn5+ZKkDh06KCoqKsgRIVSRa7AS+QZPPF1PHR0Zoc6t430+J7kGK5FvsEq45Roj1agXl8ulvXv3au/evXK5XMEOByGMXIOVyDd44mn6d5fUBDnsvi1SJpFrsBb5BquEW65RVAMAANRg96ESbS087HY8K43rqQEg3FFUAwAA1MDbImX1WfkbABAaKKoBAABq4Ol6aknKYpEyAAh7FNUAAAA18LTyt8MeoU4pvi9SBgAIDRTVAAAANcjdWuh2rGubREVF8lEKAMIdfwkAAACqsfPAEe08UOJ2PJtFygAAYp9q1FNkZKRSU1PN20CgkGuwEvmGyrxeT52eVO9zk2uwEvkGq4RbrlFUo17sdrsyMjKCHQbCALkGK5FvqCwngCt/k2uwEvkGq4RbrjH9GwAAoBq5WwrdjsVGReq4lixSBgCgqAYAAPDKMAyPe1RnpiUqMsIWhIgAAA0N079RL6WlpVqxYoUkKTMzU9HR0UGOCKGKXIOVyDdU2L7/iHYfKnU77q/9qck1WIl8g1XCLdcoqlEvhmGotLTUvA0ECrkGK5FvqOBpf2rJP9dTS+QarEW+wSrhlmtM/wYAAPDC0/7UkpTlp6IaAND4UVQDAAB44WmkOt5hV4fmcUGIBgDQEFFUAwAAeFDdImURLFIGAPg/FNUAAAAebNl3WIXFZW7Hs9OTrA8GANBgUVQDAAB44G2Rsqw0rqcGAPyBohoAAMCDHC+LlPlr5W8AQGhgSy3Ui91uV6dOnczbQKCQa7AS+QZJyvUwUp0YY1e75CZ+ew5yDVYi32CVcMu10H+FCKjIyEg1b9482GEgDJBrsBL5BpfL8yJl2elJstn8t0gZuQYrkW+wSrjlGtO/AQAAjrF5b7EOHil3O87+1ACAYzFSjXoxDEOlpaWSpOjoaL9+ew9URq7BSuQbcrYUejye7edFysg1WIl8g1XCLdcYqUa9lJaWaunSpVq6dKn5gwMEArkGK5Fv8HQ9teT/kWpyDVYi32CVcMs1imoAAIBj5Hi4njo5LlppSbFBiAYA0JBRVAMAAFTidBla6aGozkprGvJTGAEAdUdRDQAAUEn+7kMqKnW6HWd/agCAJxTVAAAAleR4u57az4uUAQBCA0U1AABAJd6K6uz0JGsDAQA0ChTVAAAAleR6uJ66ZYJDKYmOIEQDAGjoKKoBAAD+T7nTpZXb3IvqbBYpAwB4YQ92AN4sW7ZM3bt3D3YYqEF0dLR69Ohh3gYChVyDlci38LW+4JCOlLncjvt7f+oK5BqsRL7BKuGWawEfqX7ooYdUVlZW6/Yul0vPPfecTjnllABGBX+x2WxyOBxyOBx8g4+AItdgJfItfHm/njowRTW5BiuRb7BKuOVawIvql19+Waeffrry8vJqbJuXl6e+ffvqySefVHl5eaBDAwAAqCLXS1GdycrfAAAvLLmmesmSJerRo4fef/99r23eeust9ejRQ4sXL5ZhGMrKyrIiNNST0+nUnj17tGfPHjmd7nt6Av5CrsFK5Fv4yvGwSFlq0xi1SogJyPORa7AS+QarhFuuBbyonjZtmlJTU1VcXKxbbrlFV1xxhfbt22fev3XrVp177rm66667VFxcrKioKD399NP6+eefAx0a/KC8vFx5eXnKy8tjdgECilyDlci38FRa7tLq7Qfcjgdyf2pyDVYi32CVcMu1gBfV559/vlasWKGrr75ahmFoypQp6t69u+bOnauPPvpIWVlZmjlzpgzD0Omnn65ly5bp8ccfV1RUVKBDAwAAMK3beVCl5e6LlAXqemoAQGiwZPXvpKQkffjhhxo2bJhuu+02/f777zrrrLMkSYZhKD4+Xs8//7xuv/32sLiQHQAANDye9qeWpKz0JGsDAQA0KpZuqTV8+HB17dpVvXv31pEjR2QYhqKjozVr1iz17t3bylAAAACq8LbydyCnfwMAGj9LFiqr8Mknn+iss84yC2pJKisr07nnnqvx48dbGQoAAEAVuVsL3Y6lN4tVclzo77EKAPCdJUX1zp07ddlll+nKK6/Url275HA49Oabb2rmzJlKT09XYWGhbrrpJg0ePFibNm2yIiQAAADTkTKn1u446Hac66kBADUJeFE9ceJEde3aVVOnTpVhGOrZs6eWLFmi2267TWeddZZyc3M1cuRIGYahWbNmKTMzU6+++qpcLveFQgAAAAJh7Y6DKnMabsez0pKsDwYA0KgEvKgeNWqUuYXWo48+qh9//FFdunQx709MTNS4ceM0depUtWzZUsXFxXrwwQd12mmnBTo0+IHNZlN0dLSio6NZZA4BRa7BSuRb+PG0P7UU+JFqcg1WIt9glXDLNZtRcXFzgERERKh9+/aaOHGi+vXrV23b3bt3a/To0fr0009ls9nCYqPw2tiyZYvatm0rSfr999+Vnp4e5IgAAAgtD32yXJN/2eJ2fPmTg9W0Cdt8AkCoCERtFfCR6uuuu07Lli2rsaCWpBYtWuiTTz7RpEmT1LQp1zABAABreFr5u33zJhTUAIAaBbyonjBhghITE+v0mBEjRmjFihUBiggAAOAPh0udytt1yO04+1MDAGrD0n2qJWn58uWaMWOGVq5cqV27dunDDz9UUlKSPv74Yw0aNEgtW7aUJLVp08bq0OCD8vJybd26VZKUlpYmu93ylEKYINdgJfItvKzafkBOl/vVcNkW7E9NrsFK5BusEm65Ztmr++WXX3T33Xdr0aJFkiTDMGSz2VRaWipJuu2223TkyBHdd999evrpp0O+40OF0+nU9u3bJUmtW7fmfUPAkGuwEvkWXnK3FHo8nmXBdlrkGqxEvsEq4ZZrluxT/dlnn+mMM87QokWLZBiGMjIy3No4HA6VlJTohRde0KhRo6wICwAAwOPK3zab1K1N3S5fAwCEp4AX1evXr9d1112nkpISnXTSSVq+fLk2btzosd2TTz4pwzD03//+V9OmTQt0aAAAAMr1sEhZxxZxSohhkTIAQM0CXlSPHTtWxcXF6tChgxYsWKCsrCyP7eLj4/XUU0/ppptukmEYevfddwMdGgAACHNFJeXaUOC+SFk2i5QBAGop4EX1d999J5vNpueee05NmjSpsX3F1O+Ka68BAAACZdX2A/KwRpmyLFikDAAQGgJeVG/ZskWSdOqpp9aqfcVG3Hv27AlYTAAAAJLn/aklKduCRcoAAKEh4EV1xeh0xSrfNdm9e7ck1XlvawAAgLrytPJ3hE3qyiJlAIBaCvja5pmZmZo/f76mT5+uzp0719j+yy+/lKRatUXwRUREKDk52bwNBAq5BiuRb+HD08rfnVolqEm0Ndu/kGuwEvkGq4RbrgX8L8Y111yjefPm6cknn9Qll1yi9PR0r21zcnL08ssvy2azadiwYYEODX4QFRXFFyCwBLkGK5Fv4eHgkTJtLChyO27F/tQVyDVYiXyDVcIt1wL+tcFNN92k7t27q7CwUN27d9d//vMf8zpr6ei08DVr1ui5555Tv379VFRUpDZt2uj2228PdGgAACCMrdh6wONxrqcGANRFwEeqIyMjNW3aNJ111llau3atbr31VkmSzWaTJGVkZJhtDcNQ06ZN9cUXXyguLi7QoQEAgDCWu7XQ43FW/gYA1IUlE9xTU1P166+/6qGHHlJ8fLwMw3D7FxkZqSuvvFK5ubnq0aOHFWHBD8rKyrRixQqtWLFCZWVlwQ4HIYxcg5XIt/DgaeVve4RNJ6Zat0gZuQYrkW+wSrjlmjWrcOjoKuB///vfNWbMGP3www9as2aN9u/fr4SEBB133HE67bTTzIvZ0Xi4XC4dOnTIvA0ECrkGK5Fv4SHXwyJlnVMSFBMVaVkM5BqsRL7BKuGWa5YV1RViY2N19tln6+yzz7b6qQEAACRJ+4vLtHlPsdtxrqcGANRV6K9vDgAAcAxPo9SSlMn11ACAOvLrSPWNN97ot3PZbDa99957fjsfAABAhRwvi5QxUg0AqCu/FtXjx483V/WuD8MwKKoBAEDA5HpYpCwq0qYTWicEIRoAQGPm16K6f//+XovqX375RUVFRZKkZs2aKTs7W4mJidq+fbvWr1+vwsJCSZLD4dBJJ53kl+IcAADAE08rf3dpnSiH3bpFygAAocGvRfWcOXM8Hn/22Wc1d+5ctW/fXmPHjtUFF1ygyMg//mi5XC59/fXXuu+++7Rx40adcMIJmjBhgj9DAwAAkCTtOVSirYWH3Y5nMfUbAOCDgK/+vXDhQj355JNq06aNFi1apJYtW7q1iYiI0JAhQ3TqqaeqV69e+uCDD3TOOefo2muvDXR4qCe73a6MjAzzNhAo5BqsRL6FNm+LlGUHYZEycg1WIt9glXDLtYCv/v3GG2/IZrPp8ccf91hQV9aiRQs98cQTMgxD48ePD3Ro8IPIyEilpqYqNTW1yuwDwN/INViJfAttnq6nloIzUk2uwUrkG6wSbrkW8KJ63rx5kqSBAwfWqv0ZZ5whScrJyQlUSAAAIIzleBipjrZHqHMKi5QBAOou4EX17t27JUkxMTG1al/xTcaBAwcCFhMAAAhfnkaqu6YmKioy4B+LAAAhKOAT3Fu0aKEdO3ZoxYoVat++fY3tf/31V0mqcao4GoaSkhItXbpUktSjRw85HI4gR4RQRa7BSuRb6Np14Ih2HDjidjxY+1OTa7AS+QarhFuuBfwr2dNOO02GYeill16Sy+Wqtm1JSYn+9re/yWazqV+/foEODQAAhBlvi5RlBWGRMgBAaAh4UX3LLbdIkhYsWKCLL75Yv//+u8d2GzZs0HnnnWdeS3377bcHOjQAABBmPO1PLUnZ6UnWBgIACBkBn/597rnn6oYbbtC4ceP09ddfq2PHjjrttNPUqVMntW7dWkVFRVqyZIkWLlwowzAkSffcc4+5YBkAAIC/eBqpjo2K1HEt44IQDQAgFFiyadi///1vJSUl6bXXXpPT6dSCBQv0ww8/mPdXFNN2u12PP/64nnjiCSvCAgAAYcQwDI8j1d3aJMrOImUAAB9ZUlRHRkbqlVde0S233KL3339fc+bMUV5eng4dOqTExER16dJFgwYN0k033WRuEg4AAOBPOw4c0e5DJW7Hg7E/NQAgdFhSVFfo0qWLXnzxRSufEgAAQFJ111NTVAMAfMdcJwAAEBY87U8tSVlpSdYGAgAIKZaOVBcVFWnr1q3auXOneR11dfr3729BVKiPqKgodevWzbwNBAq5BiuRb6Epx8MiZXHRkerYIniLlJFrsBL5BquEW65ZUlTn5ubq7rvv1rx582r9GJvNpvLy8gBGBX+IiIhQQkJCsMNAGCDXYCXyLfQYhqHcLYVuxzPTmioiwmZ9QP+HXIOVyDdYJdxyLeBF9fr163XGGWfo4MGDtRqdBgAA8Lct+w5rX3GZ23GupwYA1FfAi+rnnntOBw4ckN1u1/33369LL71UrVq1ks0WvG+F4T8ul0tFRUWSpLi4OEVEcJk+AoNcg5XIt9DjaX9qScpKT7I2kGOQa7AS+QarhFuuBbyonjlzpmw2mx555BH99a9/DfTTwWJlZWVauXKlJKlHjx5yOBxBjgihilyDlci30ON15e+04I5Uk2uwEvkGq4RbrgX8K4Ndu3ZJkoYPHx7opwIAAPAod2uh27GEGLsymjexPhgAQEgJeFHdrFkzSVJ8fHygnwoAAMCNYRgeR6qz05tyORoAoN4CXlSffPLJkqTly5cH+qkAAADcbN5TrINH3HcUYX9qAIA/BLyovueee2QYhv7+97+rrMx91U0AAIBA8rQ/tcTK3wAA/wh4UT148GA9//zz+vnnn3XBBRdo1apVgX5KAAAAk6f9qSUpK8iLlAEAQkPAV/++/vrrJUkdO3bU999/r6ysLLVr107p6emKiory+jibzaZZs2YFOjwAABDiPF1P3axJlNKbxQYhGgBAqAl4Uf3BBx/IZrPJMAzz2ObNm7V58+ZqH8fCIQAAoL5cLkMrPEz/zkpP4rMGAMAvLBmp5o9W6HI4HDr11FODHQbCALkGK5FvoWPj7iIVlTrdjgd7f+oK5BqsRL7BKuGWawEvqsePHx/opwAAAPDI0/7UkpTFImUAAD8J+EJlAAAAweLpemqJlb8BAP4T8JFqhDan06ldu3ZJklq1aqXIyMggR4RQRa7BSuRb6Mj1UFS3iHeodWJMEKJxR67BSuQbrBJuuebXonrixIn+PJ25cjgarvLycnPRueTk5JD/gUHwkGuwEvkWGsqdLq3cdsDteHZ60waz3gu5BiuRb7BKuOWaX4vqUaNG+e2PlM1mo6gGAAA+21BQpMNl7ouUsT81AMCf/D79u/LWWQAAAMGSs6XQ43GupwYA+JNfi+rZs2f783QAAAA+y/WwP7XESDUAwL/8WlQPGDDAn6cDAADwmaeVv1snxqhVA1mkDAAQGthSCwAAhJwyp0urtrsvUsb+1AAAf6OoBgAAIWfdzoMqLXe5Hc9m6jcAwM/Ypxr1EhERofj4ePM2ECjkGqxEvjV+nvanlhreSDW5BiuRb7BKuOUaRTXqJSoqSpmZmcEOA2GAXIOVyLfGL6eRLFJGrsFK5BusEm65FvpfGwAAgLDjaaQ6LSlWzeMdQYgGABDKKKoBAEBIKSl3as0O90XK2J8aABAITP9GvZSVlSk/P1+S1KFDB0VFRQU5IoQqcg1WIt8at7U7DqrMabgdb2jXU0vkGqxFvsEq4ZZrjFSjXlwul/bu3au9e/fK5XJfZRXwF3INViLfGjdP+1NLUnZakrWB1AK5BiuRb7BKuOVa0Irq8vJyFRQUyDDcv0kGAADw1YpGskgZACA0WFpUr169Wvfdd58yMzPVpEkTpaamqqCgQJI0dOhQ/fOf/1RxcbGVIQEAgBDjaaQ6o3kTNW0S2tMPAQDBYUlR7XK59MADDygrK0uvv/66Vq1apfLy8iqj1PPnz9c999yjk046SWvXrrUiLAAAEGKOlDm1budBt+OMUgMAAsWSovrGG2/Ua6+9JpfLpZ49e+qpp55yazN8+HDFxMRow4YNOu+883To0CErQgMAACFk9fYDKne5X1rGyt8AgEAJeFH9+eefa+LEiZKkZ599VosWLdKTTz7p1u7dd9/V4sWLlZKSot9++03//Oc/Ax0aAAAIMbler6dOsjYQAEDYCHhR/c4770iSLrvsMj366KOKiPD+lN26ddOjjz4qwzD02WefBTo0AAAQYryt/J2ZlmhxJACAcBHwfapzcnJks9l0zz331Kr9wIEDJR1d1AwNX2RkpFJTU83bQKCQa7AS+dZ45Xooqju2jFNCTMNcpIxcg5XIN1gl3HIt4EX1nj17JEnt2rWrVfuEhARJUmlpacBigv/Y7XZlZGQEOwyEAXINViLfGqfi0nLl7XJfpCy7AS9SRq7BSuQbrBJuuRbw6d+tWrWSJG3btq1W7StW/m7ZsmXAYgIAAKFn1bYD8rBGmbLSkyyPBQAQPgJeVA8YMEDSH9dW1+Sdd96RzWbTqaeeGsiwAABAiPF2PTUrfwMAAingRfVdd90lwzA0ceJEvf32217bGYahMWPG6PPPP5ck3XDDDYEODX5QWlqqJUuWaMmSJUzZR0CRa7AS+dY4eVr5O8ImdU1tuIuUkWuwEvkGq4RbrgX8muo+ffrooYce0osvvqg777xTH374oc477zzz/okTJ2r37t367LPPtH79eknSsGHDdMEFFwQ6NPiBYRjmD4pheJhzB/gJuQYrkW+NU86WQrdjx7eKV5wj4B93fEauwUrkG6wSbrlmyV+Z559/XtHR0Xruuee0cOFC/fjjj7LZbJKkhx9+WNIfnX3ppZdq0qRJVoQFAABCxMEjZdq4u8jtOPtTAwACLeDTvyXJZrPp6aefVm5urm666SalpqbKMAzzX3JysoYNG6ZvvvlGn3zyiWJjY60ICwAAhIiV2w7I02AI11MDAALN0vlQXbt21bvvvitJOnTokPbv36+EhAQlJjbca50AAEDD52l/aknKoqgGAARY0C4yio+PV3x8fLCeHgAAhJAcD4uURUbYGvQiZQCA0BDw6d8RERGy2+3atWtXrdrv3LnTfAwAAEBt5HpYpKxzSoJioiKtDwYAEFYsuaa6Liu+lZeX1/kxAAAgfO0vLtOmPcVux7PTmPoNAAg8vw8H//bbbx6Pb9myRUeOHKn2sS6XS5999pkkKTo62t+hIQDsdrs6depk3gYChVyDlci3xmXFtsZ7PTW5BiuRb7BKuOWa319hhw4d3I4ZhqHevXvX+hw2m818E9CwRUZGqnnz5sEOA2GAXIOVyLfGJcfLImWNYeVvcg1WIt9glXDLNb8X1d6mbddlOndcXJyee+45f4UEAABCWO7WQrdjUZE2ndA6wfpgAABhx+9F9ezZs83bhmHorLPOks1m0yeffKLk5OQaHx8fH68uXbooLi7O36EhAAzDUGlpqaSjU/ZtNluQI0KoItdgJfKtcfE0Ut2ldaIc9oa/SBm5BiuRb7BKuOWa34vqAQMGeDzet29ftWrVyt9PhyArLS3V0qVLJUk9evSQw+EIckQIVeQarES+NR57i0q1Zd9ht+OZjWSRMnINViLfYJVwy7WAXzV+/fXXy2azKTY2NtBPBQAAwkyuh/2ppcZxPTUAIDQEvKgeP358oJ8CAACEKU/7U0tSViMZqQYANH5+LarnzZsn6ehqb3379q1yzBf9+/f3S1wAACA0ebqeOtoeoc4pLFIGALCGX4vqgQMHymazKT4+Xvv3769yrK5sNpvKy8v9GR4AAAgxnqZ/n5iaqGh7RBCiAQCEo4BsqeVyudyOAQAA+NOug0e0ff8Rt+PZTP0GAFjIr0V1fn6+JFUZma44BgAA4E8rvCxSlsUiZQAAC/m1qM7IyKjVMQAAgPrydD21xMrfAABrBXz1b4S26Oho9ejRw7wNBAq5BiuRb41DroeiOiYqQse3jA9CNL4h12Al8g1WCbdcC/gqHgUFBR6Pr169WldccYVSUlIUFxenbt266YUXXlBJSUmgQ4If2Ww2ORwOORwOnxakA2qLXIOVyLeGzzAM5XiY/t2tTVPZIxvPImXkGqxEvsEq4ZZrAfmrs3PnTv3pT39SUlKSUlNT3RYu++mnn9S7d29NmTJFBQUFOnz4sFavXq1HH31UZ555pg4dOhSIsAAAQIjYeaBEBQfdv4hnf2oAgNX8Pv171apVOvPMM7V7924ZhuH2zURpaalGjBih4uJiSVJmZqa6du2q1atXKzc3V4sWLdIdd9yhCRMm+Ds0BIDT6VRhYaEkKSkpSZGRkcENCCGLXIOVyLeGL2dLocfjje16anINViLfYJVwyzW/jlSXlZXpsssuU0FBgQzD0Nlnn62xY8cqIuKPp5k8ebLy8/Nls9n00EMPKScnRx9//LGWL1+uN954Q4Zh6IMPPtDKlSv9GRoCpLy8XHl5ecrLy2NfcQQUuQYrkW8Nn6f9qaXGV1STa7AS+QarhFuu+bWo/vjjj7V27VrZbDa9//77mj59uu68884qbf7f//t/kqR27drp2WefrXLfHXfcoTPPPFOS9NFHH/kzNAAAEEI8rfwdFx2pDi0azyJlAIDQ4Nei+tNPP5UkDRo0SKNGjXK7v7S0VLNnz5bNZtPll18uu9199vmIESNkGIbmzZvnz9AAAECIMAzD40h1t7SmiowI/QVxAAANi1+L6mXLlslms+mWW27xeP+iRYvMa6kvuugij20qll5fu3atP0MDAAAhYmvhYe0tKnU7ns0iZQCAIPBrUb1r1y5JUufOnT3ev3jxYkmS3W5Xnz59PLZp0aKFJJkXtgMAAFTmaX9qScpqZNdTAwBCg1+L6oqL0JOSkjzeX1FUn3TSSXI4HB7bVOxTHQ6bhAMAgLrztD+1JGWnJ1kbCAAA8nNRXTHKvH37drf7XC6Xvv/+e9lsNvXr18/rOSqmfbds2dKfoQEAgBDhaaQ6IcaujOQmQYgGABDu/FpUn3TSSZKkOXPmuN33888/a8+ePZJkrvDtybRp0yRJ2dnZ/gwNAWKz2RQdHa3o6Gi3PckBfyLXYCXyreEyDMPjHtVZaU0V0QgXKSPXYCXyDVYJt1xzX367Hi688EJ9++23+ve//6277rpLcXFx5n0vvPCCJKlJkyYaOHCgx8dv2rRJ48ePl81m0/nnn+/P0BAg0dHR6tmzZ7DDQBgg12Al8q3h+m1vsQ4ccd/ztLFeT02uwUrkG6wSbrnm15Hqa6+9VsnJydq8ebNOPfVUjR8/Xt98841uvvlmTZ06VTabTSNHjlRCQoLbYzdu3KghQ4boyJEjatq0qa666ip/hgYAAEKAp/2pJSk7LcnaQAAA+D9+Halu2rSp3n77bV155ZVatWqVbrrppir3x8XF6YEHHqhybNy4cZo2bZqmTZum0tJS2Ww2vfjii14XOwMAAOHL0/7UkpTdSEeqAQCNn1+Lakm6/PLLFRERoTvvvFM7d+40j7ds2VL//e9/1aFDhyrtJ0+erOnTp8swDEVEROjpp5/WzTff7O+wECDl5eXaunWrJCktLU12u99TCpBErsFa5FvD5el66qQmUUpvFmt9MH5ArsFK5BusEm65FpBXd9lll+niiy/WwoULtW3bNrVo0UJnnHGGx220oqKi1LZtWw0cOFB33nmnevXqFYiQECBOp9Nc7b1169Yh/wOD4CHXYCXyrWFyuQyt2HrA7XhWWtNGuxAOuQYrkW+wSrjlWsBend1uV//+/Wts98UXXwQqBAAAEELy9xTpUIn7ImVM/QYABJNfFyoDAAAIFE/7U0tSFouUAQCCiKIaAAA0Cl5X/makGgAQRBTVAACgUcjdWuh2rEV8tFKbxlgfDAAA/4eiGgAANHjOEFykDAAQGiiqAQBAg7eh4JAOlzndjmelJ1kfDAAAlYT22uYIuIiICCUnJ5u3gUAh12Al8q3h8Xo9dVrjvp6aXIOVyDdYJdxyjaIa9RIVFaXOnTsHOwyEAXINViLfGp7cLYUej2c18kXKyDVYiXyDVcIt10L/awMAANDo5Wx1H6lOSXQoJZFFygAAwUVRDQAAGrQyp0urtnlapCzJ+mAAADgG079RL2VlZVq7dq0k6YQTTlBUVFSQI0KoItdgJfKtYcnbeUgl5S6346GwPzW5BiuRb7BKuOUaRTXqxeVy6dChQ+ZtIFDINViJfGtYVniY+i01/uupJXIN1iLfYJVwyzWmfwMAgAYtZ2uhx+NZjXzlbwBAaKCoBgAADVquh+200pJi1SLeEYRoAACoiqIaAAA0WKXlLq3eftDtOKPUAICGgqIaAAA0WOt2HlSp0/16vFC4nhoAEBooqgEAQIOV42HqtxQaK38DAEIDRTUAAGiwclmkDADQwLGlFurFbrcrIyPDvA0ECrkGK5FvDYenkep2yU2U1CQ6CNH4H7kGK5FvsEq45Vrov0IEVGRkpFJTU4MdBsIAuQYrkW8Nw5Eyp9bu8LBIWQhN/SbXYCXyDVYJt1xj+jcAAGiQ1uw4qHKX4XY8m6nfAIAGhKIaAAA0SLlbCj0eD6WRagBA48f0b9RLSUmJli5dKknq0aOHHA5HkCNCqCLXYCXyrWHwtvJ3ZgiNVJNrsBL5BquEW64xUg0AABqk3K3uRXXHFnFKjIkKQjQAAHhGUQ0AABqcw6VOrdsZ2ouUAQBCA0U1AABocFZt3y8Pa5SxPzUAoMGhqAYAAA2Ot+ups9OTrA0EAIAaUFQDAIAGJ9dDUW2zSd3aJAYhGgAAvKOoBgAADU6Oh0XKjm8ZrzgHG5cAABoWimoAANCgHCop14aCQ27HWaQMANAQ8XUv6iUqKkrdunUzbwOBQq7BSuRbcK3cul+Gh0XKskNwkTJyDVYi32CVcMs1imrUS0REhBISEoIdBsIAuQYrkW/B5Wl/aknKCsFFysg1WIl8g1XCLdeY/g0AABoUTyt/R0bY1DWVRcoAAA0PI9WoF5fLpaKiIklSXFycIiL4ngaBQa7BSuRbcHkaqe7UKl6x0ZFBiCawyDVYiXyDVcIt10L71SHgysrKtHLlSq1cuVJlZWXBDgchjFyDlci34Nl/uEz5u4vcjmeH6CJl5BqsRL7BKuGWaxTVAACgwVgZRtdTAwBCA0U1AABoMDztTy2F5srfAIDQQFENAAAajFwPi5RFRdrUJTV8VpEFADQuFNUAAKDByNla6HbshNYJcthDb5EyAEBooKgGAAANwr6iUv2+97Db8ay0JOuDAQCgliiqAQBAg+BpKy0pdFf+BgCEBopqAADQIHgrqrNYpAwA0IDZgx0AGjeHw6FTTz012GEgDJBrsBL5Fhw5WwrdjkVHRqhzSuguUkauwUrkG6wSbrnGSDUAAGgQPK38fWJqgqLtfFwBADRc/JUCAABBV3CwRNv2H3E7nsX11ACABo7p36gXp9OpXbt2SZJatWqlyEi2PEFgkGuwEvlmvRXeFikL8ZW/yTVYiXyDVcIt1yiqUS/l5eXavHmzJCk5OTnkf2AQPOQarES+WS/Hw9RvKfRHqsk1WIl8g1XCLdeY/g0AAIIud2uh2zGHPUKdWsVbHwwAAHVAUQ0AAILO00h1tzaJskfyUQUA0LDxlwoAAATVzgNHtOtgidvx7PQk64MBAKCOKKoBAEBQeb2eOi20r6cGAIQGimoAABBUuVsKPR7PDvFFygAAoYGiGgAABFWOh+20mkRHqmNLFikDADR8bKmFeomIiFB8fLx5GwgUcg1WIt+sYxiGcj1M/85s01SREbYgRGQtcg1WIt9glXDLNYpq1EtUVJQyMzODHQbCALkGK5Fv1tm2/4j2FJW6HQ/1/akrkGuwEvkGq4RbroX+1wYAAKDB4npqAEBjR1ENAACChpW/AQCNHdO/US9lZWXKz8+XJHXo0EFRUVFBjgihilyDlcg36+R6WKQswWFX++ZxQYjGeuQarES+wSrhlmuMVKNeXC6X9u7dq71798rlcgU7HIQwcg1WIt+sYRiGx5HqzLSmigiDRcokcg3WIt9glXDLNYpqAAAQFL/vPaz9h8vcjnM9NQCgMaGoBgAAQZGztdDj8XBZ+RsAEBooqgEAQFB42p9akrLTkqwNBACAeqCoBgAAQeHpeuqmsVFqmxwbhGgAAPANRTUAALCcy2VohYeVv7PTm8pmC49FygAAoYGiGgAAWG7TniIdLCl3O87+1ACAxoZ9qlEvkZGRSk1NNW8DgUKuwUrkW+B52p9aCr+Vv8k1WIl8g1XCLdcoqlEvdrtdGRkZwQ4DYYBcg5XIt8DzdD21JGWlJ1kbSJCRa7AS+QarhFuuMf0bAABYztPK383jotWmaUwQogEAwHcU1QAAwFJOl6GV29yL6iwWKQMANEJM/0a9lJaWasWKFZKkzMxMRUdHBzkihCpyDVYi3wIrf/chFZU63Y5nh+EiZeQarES+wSrhlmsU1agXwzBUWlpq3gYChVyDlci3wOJ66j+Qa7AS+QarhFuuMf0bAABYyltRHW4rfwMAQgNFNQAAsJSn7bRaJTiUksgiZQCAxoeiGgAAWKbc6fK4SBmj1ACAxoqiGgAAWGZ9wSEdKXO5Hc9KS7I+GAAA/ICiGgAAWIbrqQEAoYaiGgAAWCbXS1GdGYbbaQEAQgNbaqFe7Ha7OnXqZN4GAoVcg5XIt8DJ8bBIWZumMWqZ4AhCNMFHrsFK5BusEm65FvqvEAEVGRmp5s2bBzsMhAFyDVYi3wKjtNyl1dsPuB3PCuOp3+QarES+wSrhlmtM/wYAAJZYt/OgSsvdFynLTk+yPhgAAPyEkWrUi2EYKi0tlSRFR0fLZrMFOSKEKnINViLfAsPT/tSSlBXG11OTa7AS+QarhFuuMVKNeiktLdXSpUu1dOlS8wcHCARyDVYi3wLD28rf4VxUk2uwEvkGq4RbrlFUAwAAS+RuLXQ71jY5Vs3ioq0PBgAAP6GoBgAAAXekzKm1Ow66Hc9OS7I+GAAA/IiiGgAABNzaHQdV5jTcjofzyt8AgNBAUQ0AAALO0/7UkpQdxtdTAwBCA0U1AAAIuNwthR6Pd6OoBgA0chTVAAAg4Dyt/N2hRZyaxkYFIRoAAPyHohoAAATU4VKn8nYdcjsezltpAQBChz3YAaBxi46OVo8ePczbQKCQa7AS+eZfq7YfkNPlvkhZNouUkWuwFPkGq4RbrlFUo15sNpscDkeww0AYINdgJfLNv7xdT81INbkGa5FvsEq45RrTvwEAQEB5WvnbZmORMgBAaGCkGvXidDpVWFgoSUpKSlJkZGRwA0LIItdgJfLNv3I9LFJ2XMt4xTv4GEKuwUrkG6wSbrnGSDXqpby8XHl5ecrLy1N5eXmww0EII9dgJfLNf4pKyrW+wH2RMvanPopcg5XIN1gl3HKNohoAAATMym0HZLivUaYsFikDAIQIimoAABAwOV4WKWPlbwBAqKCoBgAAAZPrYZGyCJvUNZWiGgAQGiiqAQBAwHhapKxzSoJio0N70RoAQPigqAYAAAFx4EiZNu4ucjvO/tQAgFBCUQ0AAAJihYep3xLXUwMAQgsbRKJebDaboqOjzdtAoJBrsBL55h+epn5LUlZ6krWBNGDkGqxEvsEq4ZZrFNWol+joaPXs2TPYYSAMkGuwEvnmHzkeRqrtETZ1aZ0QhGgaJnINViLfYJVwyzWmfwMAgIDwNFJ9QusExUSxSBkAIHRQVAMAAL8rLC7Vb3uL3Y5zPTUAINQw/Rv1Ul5erq1bt0qS0tLSZLeTUggMcg1WIt/qz9P+1JKUlZZkbSANHLkGK5FvsEq45Roj1agXp9Op7du3a/v27XI6ncEOByGMXIOVyLf6y/GySBkj1VWRa7AS+QarhFuuUVQDAAC/83Q9dXRkhDqnsEgZACC0UFQDAAC/8zT9+8TUBEXb+egBAAgt/GUDAAB+tftQibYWHnY7nsXUbwBACKKoBgAAfuV9kTKKagBA6KGoBgAAfuXpemqJlb8BAKGJohoAAPiVp5W/HfYIdUqJD0I0AAAEVmhvGIaAi4iIUHJysnkbCBRyDVYi3+ond2uh27GubRIVFUlfHotcg5XIN1gl3HKNohr1EhUVpc6dOwc7DIQBcg1WIt98t/PAEe08UOJ2PJvrqT0i12Al8g1WCbdcC/2vDQAAgGW8Xk+dnmRtIAAAWISiGgAA+E2Ol5W/s9lOCwAQopj+jXopKyvT2rVrJUknnHCCoqKighwRQhW5BiuRb77L3VLodiw2KlLHtWSRMk/INViJfINVwi3XKKpRZ6XlLs1YtVNz1+1SzpZCbdh1SOUuKdr+uzq2jFdWWqIGdG6lc7qmKNrOZIgKlfstd+sBbSw4pFKnS9GREfSbF+Sab8g135Bvvqmab/u1evtBtzZtkmLkdBmKjLAFIcKGzeVy6dChQ+ZtIJDIN1gl3HLNZhiGEewgUL0tW7aobdu2kqTff/9d6enpQYmjzOnSewvy9Z/5+dp9yH0RmmO1THDopn4ddFO/DmG94iv9Vnf0mW/oN9/Qb76h3/yjpKRES5culST16NFDDocjyBEhlJFvsEpDzrVA1FZh+VfNMAxNnTpVQ4cOVUpKiqKiotSqVSsNGTJEkydPrtO3KRdddJFsNpvmzJkTuIAbgHU7D2rYWz/o79+sqdWHJ0kqOFiiv3+zRsPe+kHrdrqPXIQD+q3u6DPf0G++od98Q78BAPCHsCuqi4uLdemll2rYsGH68ssvtWvXLpWXl6ugoEDTpk3TlVdeqcGDB+vgwZr/4G/dulUzZsywIOrg+nXzXl321kKt2HrAp8ev2HpAl721UL9u3uvnyBo2+q3u6DPf0G++od98Q78BAFBV2BXVo0eP1tSpUyVJI0eO1E8//aRNmzbp66+/Vp8+fSRJs2bN0qhRo6o9z/79+3X11VerpKR239A3Vut2HtSo93/WwZLyep3nYEm5Rr3/s/LCZHSCfqs7+sw39Jtv6Dff0G8AALgLq6J6yZIlmjRpkiTp3nvv1fjx43XKKacoIyND559/vn788Uedf/75kqRPP/1Uv/zyi/nYkpISrV27VnPnztVzzz2nrl27av78+UF5HVYpc7p0/+Rl9f7wVOFgSbnum7xMZc7QXqyAfqs7+sw39Jtv6Dff0G8AAHgWVkX1Rx99JElKTEzUM88843Z/RESE3nnnHfP/X3/9tXn7xx9/VJcuXTRw4EA9/vjj2rZtW+ADDrL3FuT7PL3PmxVbD+i9Bfl+PWdDQ7/VHX3mG/rNN/Sbb+g3AAA8C6uiumLk+cwzz1R8vOf9Mtu1a6eUlBRJR6+ZDlel5a6AfdB5b0F+yI5M0G91R5/5hn7zDf3mG/otcOx2uzIyMpSRkSG7nZ1OEVjkG6wSbrkWVkX1zp07JUnt27f32sYwDBUVFUmSYmJizOMDBw6UYRhV/uXnh+636zNW7VTBwcBcL15wsETTV+4MyLmDjX6rO/rMN/Sbb+g339BvgRMZGanU1FSlpqYqMjIy2OEgxJFvsEq45Vrof21QyaJFi+RyuardJ23x4sXmRuXZ2dlWhdbgzF23K6Dn/9e8Dfp9X3FAnyMYvs7dHtDzh2K/0We+od98Q7/5JtD9Nm9dgS7MTg3ocwAAEChhVVQnJiZWe7/T6dRf/vIXSVKTJk108cUXWxGWtmzZUu3927f/8WGmpKTE64rjUVFRiog4OvnA5XKprKys2vNW/nLB6XSqvPyPxWdythTWFHa95GzZr5wt+wP6HKGIfqs7+sw39Jtv6Dff5GzZV+VvW3V/n44VERGhqKgo8/9lZWVyubxPJ4+MjKwyFbG0tFSGYXhtb7fbzVEWwzBUWlpa7WuJjo6WzWarVew2m03R0dHm/8vLy+V0Or22r+trrRy7pBp3LPHn54j6xs775B3vk2e8T7xPx8YueX6fArF7U1gV1dUpLS3VzTffrDlz5kiSxowZoxYtWljy3G3btq1125UrV6qgoMDjfd26dVNCQoIkqaioSCtXrqz2XKeeeqp5e9euXdq8ebP5/w27DtU6JgAA6mNDwSEtXbrU/H91f5+OFR8fr8zMTPP/+fn52rvX+x7YqampysjIMP+/YsWKaj8wdurUSc2bN5d09LNC5Tg96dGjh/nhuLCwUHl5eV7bRkdHq1u3buY5W7Rood27d3ttn5ycrM6dO5v/X7t2rTm7zpOMjAylpv4xA6Cm2P35OeJYjf196tmzp/n/rVu3VhnwOFZDfp+2bt1qtm/atKlZ9FXgfWoY71Mo/Dy5XC7t33/0S+bKudYQ3qddu/w/I5eiWtL69et17bXXatGiRZKku+++Ww8++GCQowqu8vBdMwYAYLEy74MUAAA0eDajurH/EOd0OvXPf/5Tjz76qIqLi2W32/X000/rL3/5izmFoTqbNm1Shw4dJEmzZ8/WwIEDfYqjNtO/+/TpI+noFwDp6eke2/lzmkn2M9+rhMoaAGCBmKgILX/8LPP/4TQN0jAMcyQlKyur2lVyma4avtNVK6vP+1RcXKxly5ZJOppvlV+XL7HzPnkX7j9PpaWlys3NlVQ11xrC+7RlyxYdf/zxkqTff//da21VF2E7Ur1+/XrdcMMNWrBggaSjUzTGjx+vXr16WR5LXd5Ih8NR7UJrFSIiImrVrkJkZGSVBOzYMl6rt/t3P1IAADzp2CLe69+sY/8+1aTyh6/aOLaoqI7NZqvX31ZPKn/gs9vtdTp/XV9rXc5d388RNWls71Nldru9TlsENbT3qfI03Joey/vkHT9PnlWOvTa5Foz3qS6vp7bCsqieNWuWLr74YhUVFSkmJkZPPfWU7r///jq/SaEsKy0xoEX10JPa6PELTwzY+YPlmWmr9OXywK2SG4r9Rp/5hn7zDf3mm0D3W1Za04CdGwCAQAu7onrGjBkaOnSojhw5oh49eujDDz/UiSeG3geg+hrQuZUm/1L9tPT6OLdba7VKjKm5YSNzXrfUgH7wDMV+o898Q7/5hn7zTaD7rX/nlgE7NwAAgRZRc5PQsWvXLl111VU6cuSI+vXrp/nz51NQe3FO1xS1TPD/1AhJapng0OBuKQE5d7DRb3VHn/mGfvMN/eYb+g0AAO/Cqqh+5ZVXtHfvXqWnp2vatGmKi4sLdkgNVrQ9Qjf16xCQc9/Ur4OiIkMz9ei3uqPPfEO/+YZ+8w39BgCAd2Ez/dvlcmnixImSpKuuukrFxcUqLi6u9jEJCQlhXXjf1K+DvsrZphVb/XdtdVZaU90coA9mDQX9Vnf0mW/oN9/Qb76h3wIjKipK3bp1M28DgUS+wSrhlmth89XwmjVrtGPHDknSyy+/rNTU1Br/vfvuu0GOOriiIiP0jyu6K8Hhn+9eEmLs+scVJ8ke4iMS9Fvd0We+od98Q7/5hn4LjIiICCUkJCghIcFcKRcIFPINVgm3XAv9V/h/8vLygh1Co9Q5JUHjb+xd7w9RCTF2jb+htzqlJPgpsoaNfqs7+sw39Jtv6Dff0G8AALizGdXt/I0GYcuWLWrbtq0k/21QXld5Ow/qvsnLfJr2l5mWqFev6B6WH57ot7qjz3xDv/mGfvMN/eY/LpdLRUVFkqS4uLiwGNFB8JBvsEpDzrVA1FYU1Y1AQyiqJanM6dJ7C/L13oJ8FRwsqbF9ywSHburXIewXoaHf6o4+8w395hv6zTf0m3+UlJRo6dKlkqQePXrI4QjMKuuARL7BOg051yiqw1RDKaorlDldmr5yp+atK1DOln3aUHBIZU7JERWhji3ilZXWVP07t9Tgbil8cKqkcr/lbt2vjbsPqaTcJYedfvOGXPMNueYb8s035Fv9NOQPngg95Bus0pBzjaI6TDW0orqyhvwDg9BCrsFK5BusQq7BSuQbrNKQcy0QtRVfGQMAAAAA4COKagAAAAAAfERRDQAAAACAj+q30SQsUV5ebt7evn17ECNxV1JSol27dkk6en1CQ7peAqGFXIOVyDdYhVyDlcg3WKUh51rleqpynVUfFNWNQEFBgXm7T58+QYwEAAAAAEJDQUGB2rdvX+/zMP0bAAAAAAAfsaVWI3DkyBHl5uZKklq2bCm7veFMMNi+fbs5er548WKlpqYGOSKEKnINViLfYBVyDVYi32CVhpxr5eXl5kzgrKwsxcTE1PucDac6g1cxMTHq3bt3sMOoUWpqaoPaQxuhi1yDlcg3WIVcg5XIN1ilIeaaP6Z8V8b0bwAAAAAAfERRDQAAAACAjyiqAQAAAADwEUU1AAAAAAA+oqgGAAAAAMBHFNUAAAAAAPiIohoAAAAAAB/ZDMMwgh0EAAAAAACNESPVAAAAAAD4iKIaAAAAAAAfUVQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPqKoBgAAAADARxTVAAAAAAD4iKIaPps/f76uuuoqpaenKzo6WklJSerVq5f+9re/6cCBA8EODyHCMAxNnTpVQ4cOVUpKiqKiotSqVSsNGTJEkydPlsvlCnaICAPXX3+9bDabRo0aFexQEEIOHjyoN954Q6eeeqpSUlIUHR2ttm3b6rrrrtPy5cuDHR5CyPfff68rr7xSaWlpioqKUlJSknr27KnHHntMO3bsCHZ4aMQKCgpks9nUvn37Gtvu3LlTjzzyiLKyspSQkKCYmBi1b99eN954o3JycgIfbADZDMMwgh0EGhfDMPSXv/xFL774otc2GRkZ+uabb3TiiSdaGBlCTXFxsUaMGKGpU6d6bTNo0CB99tlnSkhIsC4whJXJkyfryiuvlCSNHDlS48ePD25ACAl5eXm66KKLtHbtWo/32+12TZw4UVdffbXFkSGUOJ1O3XnnnXrnnXe8tklOTtbnn3+ufv36WRgZQsUrr7yiBx98UBkZGdq0aZPXdgsWLNDFF1+svXv3erw/IiJCY8eO1Z133hmgSAOLkWrU2fjx482Cum3btvrkk0+0ceNGLVmyRLfeeqskafPmzRo6dKhKSkqCGSoaudGjR5sF9ciRI/XTTz9p06ZN+vrrr9WnTx9J0qxZsxg9RMBs3bpVo0ePDnYYCDH79+/Xueeeq7Vr16pVq1YaO3asli9frjVr1mjcuHFKS0tTeXm5brzxxmo/pAI1efXVV82CulevXvrkk0+0du1aLVmyRH/7298UHx+vvXv36tJLL2WWIepswYIFevLJJ2tst2vXLl1yySXau3evbDabHn/8ca1YsULr16/XuHHjlJSUJJfLpbvuukszZ860IPIAMIA6cDqdRnp6uiHJSE1NNbZs2eLW5sknnzQkGZKMd999NwhRIhT8+uuvZh7de++9bvc7nU7j/PPPN9v8/PPPQYgSoczpdBpnn322Icno2LGjIckYOXJksMNCCLjzzjsNSUaLFi2MvLw8t/tXr15txMbGGpKMRx99NAgRIhQ4nU6jdevWhiQjMzPTKCoqcmvz6aefmn9H33nnnSBEicZk9+7dxrJly4ypU6cao0aNMiIjI838ycjI8Pq4Rx991Gz3/vvvu92/YsUKIyYmxpBk9O3bN4CvIHAYqUadLFq0SFu2bJEkPfjgg0pLS3Nr8/DDD6tVq1aSpE8//dTS+BA6PvroI0lSYmKinnnmGbf7IyIiqkxn+/rrry2LDeHhjTfe0MyZM3XrrbfqjDPOCHY4CBF79uzRe++9J0l64YUXdPzxx7u16dKli4YOHaq4uDgtXbrU6hARInbv3m1eL33jjTeqSZMmbm2GDh2q2NhYSdKyZcusDA+N0AMPPKDu3bvrkksu0fjx4+V0Omt8jGEY+uSTTyRJ2dnZHmcXduvWTbfccosk6YcffmiU1/lTVKNOvvnmG/P28OHDPbZp0qSJLrjgAklHp+aWlpZaEhtCyy+//CJJOvPMMxUfH++xTbt27ZSSkiLp6DRdwF9Wrlyphx9+WB07dtTLL78c7HAQQr755hsdPnxYSUlJuuaaa7y2+/jjj3Xo0CG+MITP4uLizNs2m63G9k2bNg1kOAhT+fn5WrdunSTpsssu85qLl112mXl7+vTplsTmTxTVqJPNmzdLklq0aKG2bdt6bderVy9JUmlpaaP8tgnBt3PnTkmqdjVJwzBUVFQkSYqJibEiLISBkpISjRgxQqWlpRo/frzXL3UAX3z//feSpIEDB1b5vVVWVqa9e/eyowH8Ji4uTgMGDJAkffDBBx4HOb7++msdPnxYknTuuedaGh8an/Hjx8swjCr/xowZU+1jKmoHSTr55JO9tqt8X2NcS4KiGnVSUSCnpaVV+61nmzZtzNvbtm0LeFwIPYsWLdK+ffv0/PPPe22zePFiHTp0SNLRKUWAPzz55JNavny5HnjgAaZ9w+/WrFkj6egU77KyMr3yyivKzMyUw+FQ8+bN1aRJE51//vmaPXt2kCNFKPjHP/6huLg4/frrrxoyZIgWLlyogwcPaseOHfrPf/6jkSNHSpKuueYaDRw4MLjBIiRVHlzzdNlohfj4eCUmJkpqnLWDPdgBoHGpWM27Ium9qTyF6MiRIwGNCaGpphxzOp36y1/+IunoJQcXX3yxFWEhxM2bN08vvfSSunbt6vFafqC+KkZg4uPjNWjQIM2fP7/K/SUlJfr222/17bff6tlnn9Vjjz0WhCgRKnr27Kl58+bphhtu0IwZMzRjxgy3Nk888YQef/zxWk0RB+qq8k5AtakfDhw40ChrB0aq4ZOIiOpTp6b7gfooLS3VDTfcoDlz5kiSxowZoxYtWgQ3KDR6+/fv1/XXX6/IyEhNnDiRSwoQEPv375ckvf7665o/f77OPPNMzZ8/X0VFRdq7d68mT56sjIwMSdLjjz+ur776KpjhopE7cuSIJkyYYM6Q8OS///1vlTVzgEAJ5fqh8UYOICytX79e/fv316RJkyRJd999tx588MEgR4VQcNddd2nz5s164oknqr3uC6iPihGYXbt2afDgwZo+fbr69eunJk2aqFmzZrr88su1cOFCtWzZUtLRUUTAFy6XS1dccYVef/11lZaW6qKLLtI333yj/Px85eTk6KmnnlJMTIzWr1+vSy+9VFOnTg12yECjRVENoFFwOp0aO3asTjrpJC1atEh2u11/+9vf9NprrzXqbzbRMEyePFmTJk1Sr1699MgjjwQ7HISwytMfX3nlFdnt7lfitWnTRvfff7+ko9sc5efnWxYfQsfUqVP15ZdfSpLuuOMOff755zrvvPPUvn17ZWVlacyYMfr111/VrFkzuVwuPfjgg7XaIgmAOz6Jok6ioqIkqcZfupXvj46ODmhMCH3r16/XwIEDde+996q4uFjdunXTjz/+qEceeYRrwFBvO3bs0OjRo+VwODRx4kTz9xwQCMnJyZKOFs6ZmZle2/Xu3du8nZOTE/C4EHoqZnSlpqbq5Zdf9vj3smvXruZ1+xs2bGBfdPhd5b+pta0fGmPtQFGNOmndurUkad++fdW227t3r3m78krgQF3NmjVL3bt314IFCxQTE6O///3vWrp0qbltG1Bfa9as0b59+1RSUqKuXbvKZrO5/ZswYYIkacKECeYxpkrCFxXrP9S0DkR6erp5u2KXA6Au1q9fL0k6/fTTq10jYvDgwebtxriVERq2itpBqr5+MAzDrB8aY+1AUY06adeunSQpLy+v2pX5li9fLkmy2+1KTU21JDaEnhkzZmjIkCEqKipSjx49tGTJEj388MOMJAJotLKysiQdvabaMAyv7Xbt2mXeTklJCXhcCF+V90aPjIwMYiQIRRW1gyStWLHCa7sNGzaouLhYkszFGhsTimrUyXnnnSfp6OrLFdfpHMvlcmnKlCmSpDPPPFMOh8Oy+BA6du3apauuukpHjhxRv379NH/+fJ144onBDgsh6PTTT9f27dur/XfFFVdIkq644grzWMXvQ6AuBg0aJOnoZQcLFy702u5///ufpKNFDgvnwRfHH3+8JOnHH39UWVmZ13bff/+9ebtz584Bjwvh5fjjj1fHjh0l/fF7zZPK951zzjkBj8vfKKpRJ6effro5jeOxxx4ztwap7O233za3bhg2bJil8SF0vPLKK9q7d6/S09M1bdo0xcXFBTskhKjo6Gi1bt262n+xsbGSpNjYWPMYW27BF8OGDTOnft99990e/46uWLFC77zzjiTpoosuUrNmzSyNEaHh0ksvlSRt27ZNzzzzjMc2Gzdu1PPPPy9J6tatm7p27WpZfAgPNptNl112mSTp66+/9rh92+bNm/XKK69Ikvr06VPl8pfGgqIadRIZGamnnnpK0tEp4P369dOnn36q9evX6+eff9af//xn3XPPPZKk9u3ba9SoUcELFo2Wy+XSxIkTJUlXXXWViouLtWPHjmr/FRUVBTlqAKhZTEyMXn/9dUnSkiVL1LNnT02ePFkbN27UmjVr9Prrr2vAgAEqKytTkyZN9OKLLwY5YjRW11xzjbng3TPPPKNhw4Zp+vTp2rx5s1atWqWXXnpJJ510kgoKChQREaGxY8ey+CcC4v7771dSUpIkafjw4Xr++ee1bNkyrVmzRhMnTtQZZ5yhPXv2SJJZZzQ6BlBHLpfLuPvuuw1JXv+1adPGyMnJCXaoaKRWrlxZbX55+vfqq68GO2yEsJEjRxqSjJEjRwY7FISIv//979X+TktOTjZmz54d7DDRyO3YscPo06dPtbkWGxtrTJo0KdihopEaM2aMIcnIyMiott2sWbOMxMREr3los9mMl156yZqgA4CRatSZzWbT2LFjNX36dA0bNkypqamy2+2Kj49Xdna2Hn/8ceXk5JiLsQB1lZeXF+wQACCgHn74YS1atEjXXHONjjvuOMXExKhp06bKzs7WY489ptWrV2vgwIHBDhONXEpKin744QdNmjRJF1xwgfmZLS4uTllZWXrggQe0atUqXXvttcEOFSHurLPO0sqVK3XvvffqxBNPVJMmTRQdHa309HRdffXV+umnn/Tggw8GO0yf2QyjmqUnAQAAAACAV4xUAwAAAADgI4pqAAAAAAB8RFENAAAAAICPKKoBAAAAAPARRTUAAAAAAD6iqAYAAAAAwEcU1QAAAAAA+IiiGgAAAAAAH1FUAwAAAADgI4pqAAAAAAB8RFENAAAAAICPKKoBAAAAAPARRTUAAAAAAD6iqAYAAAAAwEcU1QAAAAAA+IiiGgAAAAAAH1FUAwBgkfbt28tms/n8T5J5e/z48cF9MQAAQBJFNQAAaKAGDhwom82mUaNGBTsUAAC8sgc7AAAAwsXPP/8sp9PpdnzhwoW67LLLJEmvvfaarrzySq/niIuLkyTZ7fwJBwCgIeAvMgAAFmnZsqXH48nJyebtpk2bqnXr1l7PcejQIb/HBQAAfMf0bwAAAAAAfERRDQBAI+JtobKK4zNnztThw4f1/PPP66STTlJcXJyaNWumQYMG6dtvvzXb5+Tk6Nprr1VaWpqioqLUunVrDR8+XIsWLar2+RctWqTrrrtObdu2lcPhUJs2bXTeeedp6tSpMgzD42MOHDig559/Xqeddppat26tmJgYZWRkqH///nrxxRerjL5v2rTJfC1z586VJE2YMME89tRTT7mdv7i4WC+//LL69Omjpk2bKj4+XieeeKIeeeQRbdu2zWNMFYvG/eUvf1FxcbH+9re/qVevXkpKSlJsbKw6deqku+++Wxs3bqy2PwAAkAEAAIJq9uzZhiRDkjFu3Lhq23prV3H8448/NrKyssz/H/vvk08+Mf773/8a0dHRHu+PjIw0pk+f7va8TqfTePjhh72eV5Jx3XXXGWVlZVUet379eqN9+/bVPq5Lly7Gjh07DMMwjPz8/Grbjhkzpsr5V65caXTo0MFr+5YtWxrz5893ez0ZGRmGJOP22283unXr5vXxsbGxxn//+9+a30QAQNhipBoAgBBy9913a9WqVXrooYe0bNky5efn66233pLD4ZAk3XPPPRo1apTatm2rSZMmaf369Vq7dq3+/ve/KzIyUk6nUw8++KDbeZ966im98MILkqTrrrtOCxcu1KZNmzRv3jxzkbVJkybpz3/+c5XH3Xjjjdq0aZMcDoeeffZZLVmyRPn5+Zo9e7ZGjBghSVqzZo0efvhhSVLbtm21fft2bd++Xaeddpok6YorrjCPVY5t+/btOvPMM5Wfn682bdro/fffV15entauXat33nlHLVq0UEFBgYYOHar8/HyP/fWf//xHK1eu1BlnnKEvvvhCGzZsUG5url566SUlJSXp8OHDGjFihDlqDgCAm2BX9QAAhDt/jlTr/0arj/XXv/7VvL9Zs2bmyHBlt9xyi9mmsLDQPL5ixQrDZrMZkow33njDY1zPPvusIcmw2WzGihUrDMMwjMLCQvN8f/vb3zw+7rbbbjMkGfHx8YbL5apy34ABAwxJxsiRIz0+9sorrzQkGccdd5yxd+9et/t37NhhpKWlGZKMK6+8ssp9FSPVkozhw4cb5eXlbo/Py8szWrRoYY6mO51Oj3EAAMIbI9UAAISQgQMH6oorrnA7ft5555m37733XqWkpLi1GTx4sHl737595u3XXntNhmHo1FNP1R133OHxef/yl7/o5JNPlmEY+ve//y1JVbYP2717t8fH3XvvvRozZoweeOABlZeX1/Dq/vDbb79p8uTJkqQ333xTzZo1c2uTkpKit956S5I0efJkjzHExsbqrbfeUmRkpNt9xx9/vJ577jlJR0fT582bV+v4AADhg6IaAIAQMmTIENlsNrfjLVq0MG9nZ2d7fGzlrb0qGIahr7/+WpJ0/vnnezy3JEVGRuqiiy6SJM2YMcM8X9euXSVJ//jHPzRy5Ej9/PPPcrlc5uM6d+6sp556Sk899ZSioqJq8xIlSd99950Mw5DdbtfZZ5/ttd0FF1ygyMhIGYah77//3u3+c845x+tWZ5J0zTXXmHHNmjWr1vEBAMIHRTUAACGkefPmHo9HRPzxJz8pKanGNhUOHTpkrqA9ZswYcxVuT/8qVuZet26dSktLJR1dubuiaJ04caL69OmjlJQUXXTRRXrppZe0YsUKn17nmjVrJEnl5eWy2+1eY4qKijJHzFeuXOl2nhNOOKHa54mPjzfbrFq1yqdYAQChjaIaAAB4tX///jo/xul0mtPHe/XqpTVr1ujpp5/WySefLOnoVPCvvvpKDz30kLKysjRw4ED98ssvAY9r165dbsdat25d4+Pi4+N9fk4AQOizBzsAAADQcCUkJJi3X3jhBV1//fW1elzl6ebJycl64okn9MQTT2jHjh1asGCBfvjhB02dOlWbNm3S3Llz1b9/f/3yyy/mdPHaxtWyZUvl5OTU6jGxsbFuxw4cOFDj47Zv317lOQEAqIyiGgAAeJWYmKjmzZtrz5492r17d61GdqvTunVrDR8+XMOHD9fLL7+st956S3fffbcOHz6ssWPH6l//+letznPcccdJOjrqHRcX53PBu2zZsmrvLy4uNqe/d+zY0afnAACENqZ/AwAAr2w2m/r37y9J+uyzz6qs6H2sW265RUlJSRoyZIgk6euvv9bxxx+v448/XocOHXJrHxkZqbvuuktZWVmSpE2bNtU6roqYDMPQ1KlTvbbbsGGDkpKSlJSUpAULFrjd/91333mcFl7h888/V1lZWZXnBACgMopqAABQrT/96U+SpPXr15uLkR1r5cqV+vDDD7V//36dc845kqTU1FRt2LBBGzZs0Pvvv+/xcSUlJSooKJAktWnTxmObyquFV8jOztYpp5wi6eh2Xr/99ptbG8Mw9Mwzz2j//v1yOBzq3bu3W5sjR47o1ltv9bidV0FBgR555BEztgsuuMBjfACA8EZRDQAAqnXuuefq0ksvlSQ9++yzuuSSS/Ttt98qPz9fOTk5euONNzRo0CAdPnxYHTp00E033SRJ6t69u3r16iXp6H7Ud9xxh+bOnav8/HytX79eU6ZMUf/+/bVjxw5J0nXXXVfleSsWCFu8eLFWrlyp7du3VxnxfvPNN2W327Vt2zadfPLJeuWVV5STk6P8/HzNmjVLw4cP14QJEyRJTz31lBwOh9tra9KkiaZOnap+/fppypQpWr9+vdauXatx48apV69e2rx5s6Sje3XXZcsvAED44JpqAABQLZvNpkmTJskwDH322Wf6/PPP9fnnn7u1a9u2rb766iuzGLbZbProo4909tln67ffftNbb72lt956y+NzPPvsszrrrLOqHOvdu7emTZumtWvXKjMzU9LRbb0qRstPPvlkffnll7rqqqu0e/duPfjggx7P/cADD2j06NEe77v77rv1008/ac6cORo+fLjb/ZGRkXr55Zd1+eWXe+4cAEDYY6QaAADUqEmTJpoyZYq+/PJLXXrppWrTpo2ioqIUFxenXr166dlnn9WKFSvcVu/u1KmTcnNz9eKLL6pv375q1qyZIiMj1bx5c3Xr1k033nijFi9erMcee8ztOR944AFdc801Sk5Ols1m8xjXeeedp7y8PD3++OM66aSTlJiYKLvdrjZt2uiKK67Q3Llz9fLLL3t9vMPh0MyZM/Xee+/pzDPPVNOmTeVwONS+fXvddNNNWrp0qe6999569x8AIHTZDMMwgh0EAACAldq3b6/NmzdXGfkGAMAXjFQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPmL1bwAAAAAAfMRINQAAAAAAPqKoBgAAAADARxTVAAAAAAD4iKIaAAAAAAAfUVQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPqKoBgAAAADARxTVAAAAAAD4iKIaAAAAAAAfUVQDAAAAAOAjimoAAAAAAHxEUQ0AAAAAgI8oqgEAAAAA8BFFNQAAAAAAPqKoBgAAAADARxTV/789OCQAAAAAEPT/tR/MAAAAAEwB1X4UasJG+DIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "image/png": { "height": 292, "width": 490 } }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "\n", "# Visualise the Li ion's site occupation over time\n", "plt.figure(figsize=(5, 3))\n", "\n", "# Get the site indices the Li ion occupied at each timestep\n", "li_atom = trajectory.atoms[0] # We only have one Li atom\n", "site_indices = li_atom.trajectory\n", "\n", "# Plot the trajectory\n", "timesteps = range(len(site_indices))\n", "plt.plot(timesteps, site_indices, 'o-', linewidth=2, markersize=8)\n", "\n", "# Format the plot\n", "plt.xlabel('Timestep')\n", "plt.ylabel('Site Index')\n", "plt.title('Li Ion Movement Between Sites')\n", "plt.grid(True, linestyle='--', alpha=0.7)\n", "if 'start_site' in locals() and 'end_site' in locals():\n", " plt.yticks([start_site.index, end_site.index]) # Only show the relevant site indices\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Alternative Site Types\n", "\n", "In the example above, we used spherical sites, which are intuitive but have a key limitation: they do not completely fill space. Even with appropriately sized radii, there will be regions between spheres that are not assigned to any site, and with larger radii, spheres can overlap, resulting in some regions being assigned to multiple sites.\n", "\n", "`site_analysis` provides alternative site types that offer complete spatial discretization:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|███████████████████████████████████████████████████████████████████████████████████████████████████| 11/11 [00:00<00:00, 11943.40 steps/s]" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Li ion trajectory with Voronoi sites (site indices over time):\n", "[85, 85, 85, 85, 85, 85, 101, 101, 101, 101, 101]\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\n" ] } ], "source": [ "# Example: Using Voronoi sites instead of spherical sites\n", "# Voronoi sites divide space according to proximity to site centres,\n", "# creating a complete tessellation with no gaps or overlaps\n", "\n", "# Create a trajectory using Voronoi sites\n", "trajectory = (TrajectoryBuilder()\n", " .with_structure(structure) # Initial structure for reference\n", " .with_mobile_species(\"Li\") # Specify which species we're tracking\n", " .with_voronoi_sites( # Define our spherical interstitial sites\n", " centres=site_centers)\n", " .build())\n", "\n", "# Analyse the trajectory\n", "trajectory.trajectory_from_structures(structures, progress=True)\n", "\n", "# Compare the Li ion trajectory using Voronoi sites\n", "li_atom_voronoi = trajectory.atoms[0] # We only have one Li atom\n", "print(\"Li ion trajectory with Voronoi sites (site indices over time):\")\n", "print(li_atom_voronoi.trajectory)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": "# Example: Using polyhedral sites for more complex structures\n# Polyhedral sites are defined by a set of vertex atoms forming a polyhedron\n# and provide a complete spatial discretisation based on the crystal structure\n\n# Create a reference structure with just the framework\nreference_structure = structure.copy()\n# Remove the Li ion\nreference_structure.remove_species([\"Li\"])\n# Add dummy Li ions to every polyhedron centre\nfor c in site_centers:\n reference_structure.append(species=\"Li\", coords=c)\n\n\ntrajectory = (TrajectoryBuilder()\n .with_structure(structure) # Initial structure for reference\n .with_reference_structure(reference_structure)\n .with_mobile_species(\"Li\") # Specify which species we're tracking\n .with_polyhedral_sites( # Define our spherical interstitial sites\n centre_species=\"Li\",\n vertex_species=\"O\",\n cutoff=4.0,\n n_vertices=8,\n label=\"Cubic\")\n .with_structure_alignment(align_species=\"O\")\n .build())\n\n# Analyse the trajectory\ntrajectory.trajectory_from_structures(structures, progress=True)\n\n# Compare the Li ion trajectory using polyhedral sites\nli_atom_poly = trajectory.atoms[0] # We only have one Li atom\nprint(\"Li ion trajectory with polyhedral sites (site indices over time):\")\nprint(li_atom_poly.trajectory)" }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Summary\n", "\n", "In this quickstart guide, we've:\n", "\n", "1. Loaded an existing trajectory of a simple cubic structure with a mobile lithium ion\n", "2. Defined interstitial sites at the centers of the unit cells\n", "3. Used `site_analysis` to track site occupations throughout the trajectory\n", "4. Analyzed the results to understand ion movement patterns\n", "5. Visualized the diffusion pathway and site occupations\n", "6. Demonstrated alternative site types for different analysis scenarios\n", "\n", "This demonstrates the core functionality of `site_analysis` for tracking and analyzing ionic diffusion in crystalline materials. The same approach can be extended to more complex structures and longer trajectories with multiple mobile ions." ] } ], "metadata": { "celltoolbar": "Tags", "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.1" } }, "nbformat": 4, "nbformat_minor": 4 }