examples/Relaxation of LiFePO4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "556686fe",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "This notebook shows an example of relaxing the LiFePO4 crystal.\n",
    "\n",
    "If you are running this notebook from Google Colab, uncomment the next code box to install m3gnet first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5818a2b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install m3gnet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e008446",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "from pymatgen.ext.matproj import MPRester\n",
    "from pymatgen.util.coord import pbc_diff\n",
    "\n",
    "from m3gnet.models import Relaxer\n",
    "\n",
    "for category in (UserWarning, DeprecationWarning):\n",
    "    warnings.filterwarnings(\"ignore\", category=category, module=\"tensorflow\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cc2e942",
   "metadata": {},
   "outputs": [],
   "source": [
    "mpr = MPRester()\n",
    "lfp = mpr.get_structure_by_material_id(\"mp-19017\")  # This is LiFePO4.\n",
    "\n",
    "lfp_strained = lfp.copy()  # We create a copy.\n",
    "# Create a random strain between -5% and 5% for each direction\n",
    "strains = np.random.uniform(low=-0.05, high=0.05, size=3)\n",
    "lfp_strained.apply_strain(strains)\n",
    "# In addition to the lattice strains, we also perturb the atoms by a distance of 0.1 angstrom.\n",
    "lfp_strained.perturb(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8897c0e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2022-08-16 15:10:14.833037: W tensorflow/core/platform/profile_utils/cpu_utils.cc:128] Failed to get CPU frequency: 0 Hz\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 13.1 s, sys: 2.22 s, total: 15.3 s\n",
      "Wall time: 9.48 s\n"
     ]
    }
   ],
   "source": [
    "relaxer = Relaxer()\n",
    "relax_results: dict\n",
    "%time relax_results = relaxer.relax(lfp_strained)\n",
    "relaxed_struct = relax_results[\"final_structure\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9ede3fc",
   "metadata": {},
   "source": [
    "Note that the relaxation only took < 20s."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1b355f7",
   "metadata": {},
   "source": [
    "The relaxation using the M3GNet universal IAP has brought the lattice parameters much closer to the original DFT ones and the coordinates are also within $10^{-3}$ of the original fractional coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2abb1e07",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original lattice parameters are [4.746, 10.444, 6.09]\n",
      "Strained lattice parameters are [4.655, 10.116, 6.169]\n",
      "Relaxed lattice parameters are [4.749, 10.447, 6.105]\n"
     ]
    }
   ],
   "source": [
    "print(f\"Original lattice parameters are {[round(x, 3) for x in lfp.lattice.abc]}\")\n",
    "print(f\"Strained lattice parameters are {[round(x, 3) for x in lfp_strained.lattice.abc]}\")\n",
    "print(f\"Relaxed lattice parameters are {[round(x, 3) for x in relaxed_struct.lattice.abc]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1320a18b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[9.95585866e-01, 5.09294457e-04, 9.84189370e-01],\n",
       "       [5.20894009e-01, 5.01713435e-01, 9.97474321e-01],\n",
       "       [5.04853356e-01, 4.91625924e-01, 5.07732446e-01],\n",
       "       [1.96716573e-02, 9.98423822e-01, 4.93991898e-01],\n",
       "       [4.66972193e-01, 2.26221410e-01, 7.56598929e-01],\n",
       "       [1.75367824e-02, 7.14677040e-01, 7.64243428e-01],\n",
       "       [9.74532693e-01, 2.74308700e-01, 2.60366885e-01],\n",
       "       [5.05152825e-01, 7.85728905e-01, 2.50998223e-01],\n",
       "       [4.11383909e-01, 8.64375884e-02, 2.42736461e-01],\n",
       "       [9.22660218e-01, 4.05067751e-01, 7.34201221e-01],\n",
       "       [6.64774830e-02, 5.88037118e-01, 2.49265519e-01],\n",
       "       [5.62188874e-01, 9.05333524e-01, 7.43938118e-01],\n",
       "       [7.51256294e-01, 1.04193044e-01, 2.42070362e-01],\n",
       "       [2.47872494e-01, 4.03677933e-01, 7.65556546e-01],\n",
       "       [7.64247175e-01, 6.05980547e-01, 2.46644955e-01],\n",
       "       [2.37052521e-01, 9.01273281e-01, 7.49090937e-01],\n",
       "       [1.98641403e-01, 4.49672882e-01, 2.58717668e-01],\n",
       "       [7.12504542e-01, 4.95576766e-02, 7.61240801e-01],\n",
       "       [3.04487116e-01, 9.50210774e-01, 2.57841988e-01],\n",
       "       [7.77291106e-01, 5.39120222e-01, 7.41063167e-01],\n",
       "       [2.92830014e-01, 1.70029033e-01, 6.00037370e-02],\n",
       "       [7.98603824e-01, 3.40237031e-01, 9.60862028e-01],\n",
       "       [7.93814946e-01, 3.25874071e-01, 5.43655140e-01],\n",
       "       [2.68565981e-01, 1.69799027e-01, 4.44632500e-01],\n",
       "       [2.25509920e-01, 6.64324615e-01, 4.38710904e-01],\n",
       "       [7.30509030e-01, 8.41605769e-01, 5.46991622e-01],\n",
       "       [7.04560623e-01, 8.41692824e-01, 9.60516171e-01],\n",
       "       [2.08225028e-01, 6.62225248e-01, 3.27568511e-02]])"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lfp_strained.frac_coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15e61f60",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean diff in frac. coords. between strained and original LFP structure: 0.00798\n"
     ]
    }
   ],
   "source": [
    "strained_to_orig_frac_coords_mean_diff = abs(pbc_diff(lfp_strained.frac_coords, lfp.frac_coords)).mean()\n",
    "\n",
    "print(\n",
    "    \"mean diff in frac. coords. between strained and original \"\n",
    "    f\"LFP structure: {strained_to_orig_frac_coords_mean_diff:.3}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96bd6e24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean diff in frac. coords. between original and M3GNet-relaxed LFP structure: 0.0026\n"
     ]
    }
   ],
   "source": [
    "orig_to_m3gnet_frac_coords_mean_diff = abs(pbc_diff(lfp.frac_coords, relaxed_struct.frac_coords)).mean()\n",
    "\n",
    "print(\n",
    "    \"mean diff in frac. coords. between original and M3GNet-relaxed \"\n",
    "    f\"LFP structure: {orig_to_m3gnet_frac_coords_mean_diff:.3}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1e46851",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "M3GNet reduces the mean difference in frac. coords. from the strained to the original structure by a factor of: 3.07\n"
     ]
    }
   ],
   "source": [
    "mean_pbc_diff_ratio = strained_to_orig_frac_coords_mean_diff / orig_to_m3gnet_frac_coords_mean_diff\n",
    "\n",
    "print(\n",
    "    f\"M3GNet reduces the mean difference in frac. coords. from the strained \"\n",
    "    f\"to the original structure by a factor of: {mean_pbc_diff_ratio:.3}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a064bd17",
   "metadata": {},
   "source": [
    "# Phonon Band Structure and DOS\n",
    "\n",
    "The M3GNet IAP can also be used to compute the lattice vibrations in a super cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71339aa0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING, 1 imaginary frequencies at q = ( 0.00,  0.00,  0.00) ; (omega_q = 4.480e-09*i)\n",
      "WARNING, 1 imaginary frequencies at q = ( 0.00,  0.00,  0.00) ; (omega_q = 4.480e-09*i)\n"
     ]
    }
   ],
   "source": [
    "from ase.phonons import Phonons\n",
    "from pymatgen.io.ase import AseAtomsAdaptor\n",
    "\n",
    "from m3gnet.models import M3GNet, M3GNetCalculator, Potential\n",
    "\n",
    "mp_129 = mpr.get_structure_by_material_id(mp_id := \"mp-129\")\n",
    "\n",
    "# Setup crystal and EMT calculator\n",
    "atoms = AseAtomsAdaptor().get_atoms(mp_129)\n",
    "\n",
    "potential = Potential(M3GNet.load())\n",
    "\n",
    "calculator = M3GNetCalculator(potential=potential, stress_weight=0.01)\n",
    "\n",
    "# Phonon calculator\n",
    "N = 7\n",
    "supercell = (N, N, N)\n",
    "ph = Phonons(atoms, calculator, supercell=supercell, delta=0.05)\n",
    "ph.run()\n",
    "\n",
    "# Read forces and assemble the dynamical matrix\n",
    "ph.read(acoustic=True)\n",
    "ph.clean()\n",
    "\n",
    "path = atoms.cell.bandpath(\"GHPGN\", npoints=100)\n",
    "bs = ph.get_band_structure(path)\n",
    "\n",
    "dos = ph.get_dos(kpts=(20, 20, 20)).sample_grid(npts=100, width=1e-3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "025d5dd9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.02, 'Phonon band structure and DOS of LiFePO4 Mo1 (mp-129) with (7, 7, 7) supercell')"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfYAAAE1CAYAAAAVl8KmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAB16UlEQVR4nO3ddXgUV9vA4d8TJzjBgwV3Ce7B3d3bUqhAjdK31HjrpU4F+qKlWIEixYK7huDuGjRIsED0fH/MhG+bJpDAZifZnPu6ciU7O+fMM5PdeWbOnDkjSik0TdM0TXMOLlYHoGmapmma/ejErmmapmlORCd2TdM0TXMiOrFrmqZpmhPRiV3TNE3TnIhO7JqmaZrmRJ4psYvIehF50V7BpAQR+VhEpjvbstILEQkQkRCr40gOEckjIhtF5K6IfJ/MsvdEpGhKxeYsRMRTRA6LSD6rY3kW5nocFZFcdqjrfyLy0WPeT/b+SUT+FJGOzxpbShKRiiKy1eo4UoqIFBERJSJu5usn5t0nJnYROSsiD8wdzlURmSIimewVtPZoGzdNwfqViBRPqfqtZK7bffPzeUNE1ohIjwTmaysiO8x5b4jIDBEpYPO+h4h8LyIhZl1nRWT0U4Y1GLgOZFFKvZ1ALFNE5POECiqlMimlTtvMF2nGE/fzr3WLV3fcTuCezXqMsHlfROQdETlhfq/Pi8hXIuKZQF0eInLkcQdW5oGXEpEF8aZXMqevf1y8NvN3F5GtIhKexDKDgY1KqctJqf9ZichQEdkpIhEiMiXee7VEZJWI3BSRUBH5y/aAQ0SyicgfInLN/Pk47j2lVAQwGRjBM1JKvayU+sxc5jMfEItIRaASsNB8/X68z+IDEYkVkZxJqKt+vLL3zM9HlyTGEr9sjIj8AqCU2g+EiUi7Z1hdp5LUM/Z2SqlMgD9QDfgw5ULS4os7UrNw+a5WLj8JKpmfz1LAFOBXEflv3Jsi0hWYCYwGcgLlgAhgs4hkN2d7D+OzXQPIDAQAu58ynsLAYWWf0Z++MZN93M/sJJbLZm6TXsBIEWlpTv8ZIyn2x1jPVkATYE4CdbwDhCZhWaFAbRHxsZk2ADiexFgBbmL8f0Ylcf6XgWnJqP9ZXQI+x0jC8WUHxgNFMP73d4Hfbd7/EfA2368B9BOR523enwkMSOjgymIvATPiPsdKqS9tP4vA18B6pdT1J1WklNoUr2xb4B6wPCmBxCubF3gA/GUzywwz3lTJPKB23KVvpdRjf4CzQFOb198CS8y/1wOfAVswPswrgZw287YHDgFh5rxl4tU7HNgP3AZmA1427w8CTmJ84RcB+W3eUxhf7BNm3WMASST+j4G5Zv13MXbWlWzeHwGcMt87DHSyee85YDPwHXALOAO0snnfD9hgll0F/ApMTySOnMASM96bwCaMA6tpQCzGB/Ue8B+MHYACBgLngY0YiSYksf8N4Aq8b7Muu4CCZlkF3Dfr7xG3XvHqUkBx8+8pwG9AoFmuKZAfmIexEz8DvP6Yz0wbYA9wB7gAfGzzXty6DTDX7Trwgc37Gczl3zL/H+/EX+/E4raZ1hV4CPgAApwD/hNvHhfgIPCp+XoJ8OaTvg825esAwRif3WCgjs22iwIize3dNIGyU4DPn7Q+ic1nxh73ub2BkZRzxNu+bjbzB2N810oAMUCNePUVxDjQaRzvs30EI/E/bvsHACHA/4AhNp/Fi8BIjB3/Y7dZvPpetC2TyDILYXxfbNdxCjAWWGZu9y0YCWC0+Vk6ClSJ9915z/yM3cJIxF6PW65Z7nNgyhPm8Qfu2ry+DlS3ef0+sClemRNAwwTq8jLXNaf5+gMgGqM1CIz972jbzwuQ0SwTa26Lexjf34/Nz8pUjH3EIaDaY9bjNFAvkffEfH9AUr8z8cr/Dvz+lGUHmMsWm2m+5jp7JlLmObPMXYz9Vx9z+sfY7LOJ9/3ByFtfATsw9mcLMb9r5vu1gK0Y+/V9QIDNe+uBL8zP4gOgOMZJxSqMHHAVeD+532mz3hcft42SdQQhIgWB1hg77Ti9geeB3IAHxg4EESkJ/Am8CeTCSBKLRcTDpmx3oCXGTqQixsZHRBqbG7M7kA9jxzwrXjhtgepmue5Ai8eE3gHj6C4HxtHx3yLibr53CqgPZAU+AabLP6/b1QSOYSTmb4BJIiLmezMxEmhOjC/YgMfE8DbGDjAXkAfjy62UUv0wElw7ZRyRfmNTpiFQ5gnrFmcYxtlZayAL8AIQrpRqYL5fSSXvjK83xocyM8YHdzHGB9cX4wzvTRFJLK77GGeE2TCS/Cvy7+t09TDOsJtgnFGWMaf/Fyhm/rTg8ds0MQsBN4yzo1IYicD26B6lVCzGgUozc9J2YJiIvCoiFWz+x/8iIjmApRhnvz7AD8BSEfFRSj2HcfYQd6a9+inif5LXgI4Yn4/8GIlpTAJxiojUxdiZ7MHY1iFKqR228ymlLmCsfzObyb9gfEYfJDGmqRj/czD+bwcxznLjYkl0myWxflsVgNNKqeh407tjtCbmxDhQ2YZxIJ8T4+D+h3jz9zFjLQaUxH4tkQ0wkqYtifd3+XjvH8Fo9v4HpdRDjIOghuakhhj7w7o2rzfEK3Mf44Dskvr/M924/0V7jH1pNowTpl8TWgERyYixXz6W4Boa+8zcGN+hZDHr7gr8kdyypgHAVGVmOQCl1EWMA+pSiSzvZ4yTsswYB5h7k7G8/hj703wYB1U/m/X6YnymP8fILcOBefLP/hL9MFrIMmMk8tUYrRT5MRL9GnO+JH2nkyqpif1vEQnDOHvdAHxp897vSqnjSqkHGEcZlc3pPYClSqlVSqkojLPeDBgbNc7PSqlLSqmbGIkjrmwfYLJSarcyrkG9h9HUV8Sm7CilVJhS6jywzqZsQnYppeaacfyAcRRcC0Ap9ZcZQ6yZ9E5gJIQ455RSE5RSMRgfxHxAHhEphHFg8ZFSKkIptdFch8REmWULK6WilNE09aSm2o+VUvfNbfskLwIfKqWOKcM+pdSNJJRLzEKl1BYzAVYAcimlPlVKRSrjGvAEoGdCBZVS65VSB8xtuh/jAK9hvNk+UUo9UErtwzhgiNupdQe+UErdNBPOz8kN3Pw/X8f4ssVd/0voWuxlm/e/wmha7APsBC6KSGIHFW2AE0qpaUqpaKXUnxhnhClxjW+4iISZP3FNni9jtHKEmN+Pj4Gu8S7ZXMc4K5gIjFBKrcFY18SuST/aFiLSCXBVSi1IZN5/UUptBXKISCmMHeHUeLPYc5tlwzjzim+BUmqXmQwXAA+VUlPN7+5soEq8+X9VSl0w9z9fYBwYPxPzuvRIjJamOMuBESKSWYy+Li9gNM3bumuuV0I2AA3N/29FjO9EQxHxwtgHbUxGiJuVUoHmNplGAgcTprhYEtrOYCTXuUqpe8lYdpzOGJ/PDU+aMT4RKYyxL0nooOBx2zAWKC8iGZRSl5VS8Q+8HmeaUuqgecD0EdDdvDzZFwg0t2esUmoVxr6jtU3ZKUqpQ+ZBaFvgilLqe6XUQ6XUXaVUkDlfUr7TSZbUxN5RKZVNKVVYKfVqvERzxebvcCCuY11+jCNL4NEZ0gWMM77klr2H0TyRlLIJuRAvjhBzGYhIfxHZG7fzxDiStu0McsWmbLj5Zyaz/C3znx3nHIn7FuPSwkoROS02HZqSEncSFMRofbAX22UXBvLbJJgwjLO5PAkVFJGaIrLO7Eh0G+NDG7+DzeP+97bLftw2TZDZGpMLI7HFJcOEek/ni3tfKRWjlBqjlKqLsXP4Aphs05Jg6x+fT5s4fROY91l9Z373siml4rZhYWCBzf/iCEYTu+3/I6dSKrtSqoxSKu7g6DoJbwfM6dfNs5tvgNefItZpwFCgEUZitWXPbXYL4wwovqs2fz9I4HX8fUT8z1ncPmGZTSetPkkNykzay4A3lFKbbN563Vz+CYzWpD8x9kG2MmM05yZkA8YlD3/gAEZTbkOMk5OTyTyAj/+980okecTF8q/tLCLeQDfseMadDP0wDk7OJPBegtvQ3Ef3wNgPXRaRpSJSOhnLjP85ccfYnxUGusXbL9bjn98x27KP20cn5TudZCl5Mf8SRrCA0SyIsWIXn6JsRozmu6SUTUhBm7pcgALAJfPobwLGzshHKZUNowkx0WZYG5eB7GZscQolNrN5dPa2UqooRnPYMBFpEvd2YsVs/r6PzVG+ecRo2+RzAaNJMSni15X3Ccu+AJyxSTDZlFKZlVKtEygHxiWKRUBBpVRWjOuvSdmmYGzXgjavE92mj9EBo8lsB0ZTYgjGjugR83PQhf9vCnvEbEkYg5FAyiZQ/z8+nzZxPu3nM7kuYDQr2v4/vJTRHPk4a4GCImLbIhV3ia0WxrYogXFNb5OIXAHmA/lE5Eq8FrOETANexTiLCY/3nj232X7Azw6dSuN/zi4BKKVaqf9vwp6RlIrMfclq4DOl1D869ZmtT32UUnmVUuUw9rs74lVRBqPlKiFbMZqYOwEblFKHzXhbk/hZ79MkTduY72MkoZIJvN0J46B5fXLrNT9rAfy7RSep+pPAAYXZLO5BIpcOlFIrlFLNMJLuUYz9PsTbF2L0y4gv/uckrkXwAsbZvO33MKNSyrYDaPz9aGK3sj7tdzpBKZnY5wBtRKSJeQb1NsZ1r6Tcb/gn8LyIVDZ7in4JBCmlzj5lLFVFpLO5I3jTjGM7RicThdnz1+ypGv/aV4KUUucwml0+EeO2oHo8pllRjNutipsHOLcxjsZizbevkvg/PM5xjKPrNub2/BCw7UU7EfhMREqY11Yr2ly/jF//PqCcuX29MJp9HmcHcFdE3hWRDCLiKiLlRaR6IvNnBm4qpR6aSaT3E+q3NQd4T0Syi3E72mtJLSgiOcwzrDHA10qpG+ZZwXDgQxHpLSJe5oHMRIy+CD+aZd8U4xahDCLiZjbDZ+af/UniBAIlzfrcxLgFrSxGB7ykcjVjifvxeHKRR/4HfGEmE0Qkl4h0eFIhpdRxs+wMMW7RchWRchjXSVcroz/AQYwdWWXz50WMz09lntCCZJ5FNcTo4BXfY7eZGYsXRt8IF3ObuCdQD0qpEIzWrxoJvZ8MQ0SkgBjX/z/AaK5PkBmzF0bHwLj/Xdx9xb4YB02/KqX+l0DZYiLiY65jK4xrrp/bvO+Lcdloe0LLNg+SdgFD+P9EvhXjDDSxxH4V8BGRrImv/hMF8u9LaJDIGbcY98mvf0Kd/YCtSql/nLma373HHoyISB2MFp6/Eni7IbDWbMaOXy6PiHQwT8IiMDoTxu179wINRKSQua3eS6DuviJS1myp+BTjEkQMMB1oJyIt4j6/5noUSKAOMD7r+cx9jacYl2Zqmu891Xc6MSmW2JVSxzCuQfyCcXTTDqODWGQSyq7GuJYxD+MMrhiJXM9NooUYTTG3MD5YnZVxnfsw8D1GJ5urGNeStySj3t4YnetuYnT6etxRaAmMI/p75vLGKqXWme99hZF4wkRkeEKFlVK3Mc6GJmKc5dznn815P2AkxZUYvTcnYfRpACNx/2HW393cwX9qxnMCo+9EoswPcVuMnfsZjP/nRIwOhwl5FfhURO5iXG9M6FaqxHyC0dx1xlyXpNzStE9E7mHs7F8E3lJKjbSJfzbG//0tjEs6hzG2TV2bZsxwjM/CFXP9hgBdlHlPuS2zTFuMg9UbGHcytFVJuO3HxgiM5tm4n7XJKPsTRovISnMbb8f4HCbFUIz/3XT+/3aj9RitF5jXv6/E/WB8tmPN1zFPqlwptVn9f0ct2+lP2mb9MLbDbxgdsx7w/2dVCRlnlnkWMzE+Y6cxzk4THFvA9KEZ0wiM/doD/r+z3YsYB84fi8291jZlq2I0od/F+K73iXeNtzfwR0JJycYGjCbgHTavM5PI9XWl1FGME6TT5vc+/2PqTsx4oI95MgI8OghpTML7uoI8ef+Z4Bm3WfZJJ30DgPlKqYSu+/fBSI4JccHoXHwJ4/PcEHgFwLwuPhujFWgXCR+cT8O44+AKRv+s182yFzBaB9/HODm8gNG3IsG8asbdDCMXXsHY9zYy336W7/S/yNNd5tA0TbOO2ZK3B2iinmKQGhE5i3HLUErctZCcODwxWtAaKKWuWRlLQkRkJjBHKfV3Eubdi/H/SHanXRGZCPyllFrxFGUrAuOUUrWTWzYJda/HuB1uor3rTkmWDnyiaZr2NMyz24T6P6Qp5nokpyOXQymlknwZTSlV+RmW89RDkyvjzhu7J/W0TD8ERtM0TdOciG6K1zRN0zQnos/YNU3TNM2J6MSuaZqmaU5EJ3ZN0zRNcyI6sWuapmmaE9GJXdM0TdOciE7smqZpmuZEdGLXNE3TNCeiE7umaZqmORGd2DVN0zTNiejErmmapmlORCd2TdM0TXMiDkvsItJSRI6JyEkRGZHA+54iMtt8P0hEipjTa4jIXvNnn4h0silzVkQOmO/tdNS6aJqmaVpq5ZCHwIiIK3Ac4yHzIUAw0EspddhmnleBikqpl0WkJ9BJKdVDRLyBSKVUtIjkw3h2cX7z9VmgmlLqeoqvhKZpmqalAY46Y68BnFRKnVZKRQKzgA7x5ukA/GH+PRdoIiKilApXSkWb070A/Tg6TdM0TUuEoxK7L3DB5nWIOS3BecxEfhvwARCRmiJyCDgAvGyT6BWwUkR2icjgFIxf0zRN09IEN6sDSAqlVBBQTkTKAH+IyDKl1EOgnlLqoojkBlaJyFGl1Mb45c2kPxggY8aMVUuXLu3Q+LV/unHjBgA+Pj4WR+I89Da1P71NE7dr167rSqlcttNy5sypihQpYlFEzimh7ZwUjkrsF4GCNq8LmNMSmidERNyArMAN2xmUUkdE5B5QHtiplLpoTr8mIgswmvz/ldiVUuOB8QDVqlVTO3fqfnZWmjJlCgDPPfecpXE4E71N7U9v08SJyLn404oUKYLet9pXQts5KRzVFB8MlBARPxHxAHoCi+LNswgYYP7dFVirlFJmGTcAESkMlAbOikhGEclsTs8INAcOOmBdNE3TNC3VcsgZu9mDfSiwAnAFJiulDonIpxhn3ouAScA0ETkJ3MRI/gD1gBEiEgXEAq8qpa6LSFFggYjErcdMpdRyR6yPpmmapqVWDrvGrpQKBALjTRtp8/dDoFsC5aYB0xKYfhqoZP9INU3TNC3t0iPPaZqmaZoT0Yld0zRN05yITuyapmma5kR0Ytc0TdM0J6ITu6ZpmqY5EZ3YNU3TNM2J6MSuaZqmaU5EJ3ZN0zRNcyI6sWuapmmaE9GJXdM0TdOciE7smqZpmuZEdGLXNE3TNCeiE7umaZqmORGd2DVN0zTNiejErmmapmlORCd2TdM0TXMiOrFrmqZpmhPRiV3TNE3TnIhO7JqmaZrmRHRi1zRN0zQnohO7pmmapjkRndg1TdM0zYnoxK5pmqZpTkQndk3TNE1zIjqxa5qmaZoT0Yld0zRN05yIwxK7iLQUkWMiclJERiTwvqeIzDbfDxKRIub0GiKy1/zZJyKdklqnpmmapqU3DknsIuIKjAFaAWWBXiJSNt5sA4FbSqniwI/A1+b0g0A1pVRloCUwTkTcklinpmmapqUrjjpjrwGcVEqdVkpFArOADvHm6QD8Yf49F2giIqKUCldKRZvTvQCVjDo1TdM0LV1xVGL3BS7YvA4xpyU4j5nIbwM+ACJSU0QOAQeAl833k1InZvnBIrJTRHaGhobaYXU0TdM0LXVKE53nlFJBSqlyQHXgPRHxSmb58Uqpakqparly5UqZIDVN0zQtFXBUYr8IFLR5XcCcluA8IuIGZAVu2M6glDoC3APKJ7FOTdM0TUtXHJXYg4ESIuInIh5AT2BRvHkWAQPMv7sCa5VSyizjBiAihYHSwNkk1qlpmqZp6YqbIxailIoWkaHACsAVmKyUOiQinwI7lVKLgEnANBE5CdzESNQA9YARIhIFxAKvKqWuAyRUpyPWR9M0TdNSK4ckdgClVCAQGG/aSJu/HwLdEig3DZiW1Do1TdM0LT1LE53nNE3TNE1LGp3YNU3TNM2J6MSuaZqmaU5EJ3ZN0zRNcyI6sWuapmmaE9GJXdM0TdOciE7smqZpmuZEdGLXNE3TNCeiE7umaZqmORGd2DVN0zTNiejErmmapmlORCd2TdM0TXMiOrFrmqZpmhPRiV3TNE3TnIhO7JqmaZrmRHRi1zRN0zQnohO7pmmapjkRndg1TdM0zYnoxK5pmqZpTkQndk3TNE1zIjqxa5qmaZoT0Yld0zRN05yITuyapmma5kR0Ytc0TdM0J6ITu6ZpmqY5EZ3YNU3TNM2JOCyxi0hLETkmIidFZEQC73uKyGzz/SARKWJObyYiu0TkgPm7sU2Z9Wade82f3I5aH03TNE1LjdwcsRARcQXGAM2AECBYRBYppQ7bzDYQuKWUKi4iPYGvgR7AdaCdUuqSiJQHVgC+NuX6KKV2OmI9NE3TNC21c9QZew3gpFLqtFIqEpgFdIg3TwfgD/PvuUATERGl1B6l1CVz+iEgg4h4OiRqTdM0TUtjHJXYfYELNq9D+OdZ9z/mUUpFA7cBn3jzdAF2K6UibKb9bjbDfyQiktDCRWSwiOwUkZ2hoaHPsh6apmmalqqlmc5zIlIOo3n+JZvJfZRSFYD65k+/hMoqpcYrpaopparlypUr5YPVNE3TNIs4KrFfBAravC5gTktwHhFxA7ICN8zXBYAFQH+l1Km4Akqpi+bvu8BMjCZ/TdM0TUu3HJXYg4ESIuInIh5AT2BRvHkWAQPMv7sCa5VSSkSyAUuBEUqpLXEzi4ibiOQ0/3YH2gIHU3Y1NE3TNC11c0hiN6+ZD8Xo0X4EmKOUOiQin4pIe3O2SYCPiJwEhgFxt8QNBYoDI+Pd1uYJrBCR/cBejDP+CY5YH03TNE1LrRxyuxuAUioQCIw3baTN3w+BbgmU+xz4PJFqq9ozRk3TNO3ZTZw4kZCQED7++GOrQ0mX0kznOU3TNC1t+OSTT/jkk0+sDiPd0old0zRNsytXV1erQ0jXdGLXNE3T7Mrd3d3qENI1ndg1TdM0u/Lw8LA6hHRNJ3ZN0zTNrjJmzGh1COmaTuyapmmaXWXJksXqENI1ndg1TdM0u8qdWz9B20o6sWuapml2VaBAAQCUUhZHkj7pxK5pmqbZVZ48eQC4du2axZGkTzqxa5qmaXYV9wTtDRs2WBxJ+qQTu6ZpmmZ3rq6uzJs3z+ow0iWd2DVN0zS7q1ixIsuWLeP8+fNWh5Lu6MSuaZqm2Z23tzf+/v68++67VoeS7ujErmmapqWIOnXqsHHjRsaNG2d1KOmKTuyapmlaivD09KRz586MGDGCP//80+pw0o0kPY9dRIYlYbb7Sil9WKZpmqY9kjNnTnr27Mkbb7zBkSNH+OSTTx71mtdSRlLP2N8BMgGZH/PzdkoEqGmapqVtefPmpX///syYMYPatWuzb98+q0Nyakk6YwemKaU+fdwMIqJH/dc0TdMSlDlzZvr27cvu3bsJCAigR48efP755+TMmdPq0JxOUs/Yf3jSDEqp/zxjLJqmaZoTc3FxoVq1agwePJj9+/fj5+fHwIEDOXLkiNWhOZWknrHvFZGDwJ/APKVUWMqFpCVF2MMwNp3bxM5LO8nimYU8mfKQN1NeKuSuQJ5MeawOT9M0LVHe3t60aNGCunXrsmvXLurUqUPVqlUZNmwYzZo1w93d3eoQ07SkJnZfoCnQE/hSRLZjJPmFSqkHKRWc9k+RMZH8HPQzfx78kz2X96D49wMW3Fzc6Fq2K6/VeI3aBWrrTiqapqVamTJlomHDhtSpU4cDBw7w+uuvc/PmTbp27Ur//v2pU6cOLi765q3kSlJiV0rFACuAFSLiAbTCSPKjRWSNUqpPCsaoAWtOr2FI4BCO3ThG3YJ1+W/D/9LIrxE1fGsQGRPJlXtXuHT3EouOLWLynsnMOjiLqvmqMrbNWGr41rA6fE3TtES5u7vj7++Pv78/N2/e5NChQ/Ts2ZOYmBh69+5Np06dqFWrFm5uST0XTd+SvZWUUpEichg4AlQFytg9Ku2R2w9v89KSl5h9aDbFshcjsHcgrUq0+sc8Xm5eZPHMQkmfkgQUCeDTRp8yff90vtr8FXUn1+Xrpl/zVq239Nm7pmmpXo4cOahfvz716tXj6tWr7Nq1i7lz53L79m1atGhBp06daNGiBdmzZ7c61FQryW0cIlJQRN4Rkd3AErNse6WUf4pFl87dfHCTptOaMv/IfD4J+ISDrx78V1JPSCaPTLxc7WX2vrSXtiXb8vbKt+kwqwM3H9x0QNSapmnPTkTImzcvjRo14oUXXuCFF17g/v37fPHFFxQoUIDatWvz5ZdfsmvXLmJjY60ON1VJ6gA1WzGus88BBimldqVoVBqh90NpOq0px64fY0GPBbQp2SbZdWTPkJ353efzy45fGL5yOLUm1mLDcxvIlzlfCkSsaZqWcrJmzUr16tWpXr06kZGRnDlzhqVLlzJmzBju379Po0aNaNOmDc2aNaNw4cJWh2uppDbFjwA2KaX+3VtLs7vLdy/TZGoTzoadZXGvxTQr1uyp6xIRXq/5Ov75/Gk5vSXNpjVjw3Mb8PH2sWPEmqZpjuPh4UGpUqUoVaoUALdv3+bUqVOMHTuWt99+mxw5ctC8eXNatmxJo0aNyJYtm7UBO1iSmuKVUhuVUkpESorIGvPWN0Skooh8mJQ6RKSliBwTkZMiMiKB9z1FZLb5fpCIFDGnNxORXSJywPzd2KZMVXP6SRH5WZzgIvKDqAe0nNGSC3cusKzPsmdK6rbqFarHol6LOHnzJC2mt+D2w9t2qVfTNM1qWbNmxd/fnw4dOvDmm2/SvHlzzpw5w4cffkj+/PmpXLky7733HuvWrSMiIsLqcFNccu8jmAC8B0QBKKX2Y/SOfywRcQXGYPSmLwv0EpGy8WYbCNxSShUHfgS+NqdfB9oppSoAA4BpNmV+AwYBJcyflslcn1TnrRVvsf/qfuZ0nUPDIg3tWndjv8bM7zGf/Vf303pma+5H3rdr/ZqmaVZzcXEhf/781K1bl+7duzNs2DAqV67Mtm3bGDRoEDly5KBp06b8+uuvnDp1yupwU0RyE7u3UmpHvGnRSShXAziplDqtlIoEZgEd4s3TAfjD/Hsu0ERERCm1Ryl1yZx+CMhgnt3nA7IopbablwimAh2TuT6pyuyDsxm3axz/qfOfJHWSexqtS7Tmzy5/su3CNl5e+jL66oqmObf0/h13d3fHz8+PRo0a0a9fP1577TWyZ8/O1KlTqV69OkWKFGHIkCEsW7aMhw8fWh2uXST3drfrIlIMjJFRRKQrcDkJ5XyBCzavQ4Caic2jlIoWkduAD8YZe5wuwG6lVISI+Jr12Nbpm9DCRWQwMBigUKFCSQjX8U7ePMmgxYOoXaA2nzf+PEWX1aVsFz4J+ISR60dSv1B9BlcdnKLL0zQt5Sml2H15N3/s+4PAE4GEPQwjPCqcB9EPKJy1MDV8a1DDtwYBRQKomq9qur39NUOGDJQvX57y5csTGxvL1atXOX78OG+88QZXrlyhVatW9O3bl+bNm+Pp6Wl1uE8luYl9CDAeKC0iF4EzgEMGpxGRchjN882TW1YpNR4jbqpVq5bqDl8joiPoMbcHbi5u/NnlT9xdU344xQ8afMDmC5t5fdnrVMtfDf98+q5F7dkduHqA5SeXc/LmSU7cPEHInRBK+pSkpm9NavjWoF6hemT00M+Lsre5h+fy8fqPORR6CA9XD1oWb4lvZl8yumfE082TEzdPEHwxmL8O/wVAKZ9S9K/Un74V+1Ioa+o82XEEFxcX8uXLR758+ahXrx53797lyJEjvP3221y9epW2bdvSv39/mjVrlqZGwEtWYldKnQaamk9yc1FK3U1i0YtAQZvXBcxpCc0TIiJuQFbgBoCIFAAWAP2VUqds5i/whDrThB+3/8juy7tZ0GMBhbM55jYNF3FheqfpVBlXhW5/dWPX4F1k88rmkGVrziVWxRJ4IpAft//I2jNrAcjpnZPiOYpTMU9FDoceJvBEIApFvkz5GN1yNN3Kdku3Z4z2FBEdwfCVw/k1+Fcq5anEb21+o0e5HmTPkPDgLaH3Q1l0bBHT9k/jg7Uf8OHaD2lbsi1DawyladGmuEjaSV4pIXPmzNSoUYMaNWpw584dDh8+zODBg/H09OSjjz6iZ8+eaWIc+yT9F0Wkre1rpdT9+Ek9/jzxBAMlRMTPHJK2J7Ao3jyLMDrHAXQF1po98bMBS4ERSqktNjFcBu6ISC2zN3x/YGFS1ic1uXz3Mp9v/JwOpTrQsXRHhy47V8ZczOk2h/O3z/PiohfT/bU4Lfm2h2ynwm8VaPdnO45dP8aoJqO4Ovwqoe+Esm3gNuZ2n8vhIYe59e4tAnsHki9zPnrM7UHrma05feu01eGnaefCzlH/9/r8Gvwrw2oNI3hQMC9XeznRpA7Gd36g/0DWP7eeM2+c4YP6HxB0MYgW01tQZkwZftj2A6H3Qx24FqlXlixZqFWrFs8//zw1a9bkiy++oHDhwowZMybV7yuTenj2rYhUERH/xH6ALxMrrJSKBoZijDd/BJijlDokIp+KSHtztkmAj4icBIZh3DuPWa44MFJE9po/uc33XgUmAieBU8CypK966vD+2veJio3iu+bfWbL8OgXr8Hmjz5l3ZN6jZjpNe5KI6AjeX/M+dSfX5V7kPWZ2nsmZN87wbr13yZ0x97/mz+qVlVYlWrHjxR381PIntpzfQoXfKrDx3EYLok/7Lty+QK1JtTh24xjzus/j+xbfJ/sSXpFsRfis8Wecf/M80ztNxyeDD2+vfBvfH3zpOqcr84/MJ+xhWMqsQBoiIhQvXpxevXrRpk0bvvvuO7p3756qO9oltSn+Kk9+JvuJx72plAoEAuNNG2nz90OgWwLlPgcS7E2mlNoJlH9CXKlW8MVgpuydwn/q/IfiOYpbFsfbdd7mr8N/MTRwKI39GpPTO6dlsWip35HQI/SY24MD1w4wsMpAfmjxA1k8sySprKuLK6/XfJ3OZTrTbFoz2s5sy7oB66iav2oKR+087kfep8OsDtyPvM+WF7ZQIU+FZ6rP082TPhX70KdiHw5dO8TkPZOZun8q847Mw0VcqJ6/Og0KN6BAlgLky5SPvJny2mlN0p4CBQrQp08fli5dSr169diyZUuq7GCX1Ke7BaRwHOmOUoo3lr9Bnox5+KDBB5bG4ubixuQOk6k6vipvrXiLaZ2mPbmQli6tOb2GLnO64OnmyZJeS55qqGOAAlkKsKrfKupNrkfLGS3Z+NxGyuTSz5N6klgVy3MLn2Pvlb0s7rX4mZN6fOVyl+P7Ft8zqukotodsZ/Xp1aw6vYoft/9IdGxS7mx2fu7u7nTo0IEpU6YQHBxMvXr1rA7pX9J3TwkL/XnwT7aFbOPLJl8m+WwnJVXMU5H3673P9P3TCTwR+OQCWrozcfdEWs5oScGsBdnx4o6nTupxCmQpwOr+q3EVV5pNa8a5sHN2itR5fbrhU+Yensu3zb595u3/OO6u7tQvXJ9PGn3C1oFbifgwgmvDr7H/5f2s6LsixZabVogIuXLlYufOnVaHkiCd2C0QHRvNh2s/pEreKjxX+Tmrw3nkgwYfUC5XOV5a8pIeclZ7RCnFB2s+YNDiQTTxa8Lm5zfb7e6N4jmKs7LfSu5F3qPfgn7EKv2UrsSsPbOWTzZ8wnOVn2NY7WEOXbaLuJArYy4q5KlA82LJvuPY6dy+fZsTJ07QrJl9hvy2N53YLfDXob84E3aG/zb8b6q6vcTD1YPJHSZz8c5FPl7/sdXhaKmAUophK4bx5eYvGeQ/iCW9l5DVK6tdl1ExT0V+avkTm85vYtzOcXat21lExkQyNHAoftn8GNt6rL5V0EJ37txh4cKFvPHGG5QrV87qcBKUrKwiIt1EJLP594ciMt/sEa8lkVKKr7d8TZmcZWhXqp3V4fxLDd8aDPIfxC87fuHQtUNWh6NZKFbF8tqy1xgdNJo3ar7BuLbjcHNJ7phWSdO/Un+aF2vOf1b/h/O3z6fIMtKyn4N+5sj1I/zU8icyuGewOpx06+TJk/z+++/07duXkSNHPrmARZJ7uviRUuquiNQDmmLcovab/cNyXitOrWDf1X38p+5/UtXZuq0vmnxBFs8svL789VR/v6aWMmJVLK8ufZUxwWMYXns4P7b4MUXPEkWEcW3HoZTi5SX6GQa2Lt29xCcbPqFNiTap8mQgPQgLC2Px4sWsXLmS+fPn88knn+Dq6mp1WIlKbmaJMX+3AcYrpZYCHvYNybmN2jyKAlkK0LtCb6tDSVRO75x81ugz1p5Zy7wj86wOR3MwpRRvLHuDcbvG8V699/im2TcOafotkq0IXzb5kmUnlzHzwMwUX15aMXzlcKJiovip5U9Wh5LuPHjwgLVr1zJ58mSaN2/OiRMnCAgIsDqsJ0puYr8oIuOAHkCgiHg+RR3p1rYL29hwbgNv134bD9fUfTz0UrWXqJinIsNWDNOPd01nPlz7Ib8G/8rbtd/mi8ZfOPR67pDqQ6hdoDZvrniTe5H3HLbc1GrD2Q38efBP/lP3PxTLUczqcNKNe/fusXbtWsaOHYufnx+HDx/myy+/JEsW6+9gSorkJuXuGKPHtVBKhQE5gHfsHZSz+nrL1+TIkIMX/V+0OpQncnNx45dWv3DhzgVGbR5ldTiag4zaPOpRR7lvm33r8E5ari6u/NDiB66HX+e3YH2V77/r/4tvZl9G1Bvx5Jm1Z3br1i1WrFjBuHHjKF26NPv27WPKlCn4+ib44NBUK1mJXSkVDlwD4u7Ij+YJI85phqPXj7Lw2EJeq/EamTwyWR1OkjQo3IBe5Xvx7dZv9T3G6cBvwb/x3pr36F2hN7+1+c2ynte1CtSiebHmfLv1W8Kjwi2JITUICgliw7kNDKs9DG93b6vDcWohISH8/fffTJkyhXr16nHixAnGjRuHn5+f1aE9leT2iv8v8C7wnjnJHZhu76Cc0fhd43F3cefV6q9aHUqyjGo6ChFhxBp9xuDM/jr0F0MCh9C2ZFumdJiCq4u1HYNGNhhJaHhour797estX5PNKxuD/AdZHYpTio2N5ejRo0yfPp3AwEAGDBjAhQsX+Oabb8iTJ4/V4T2T5DbFdwLaA/cBlFKXgMz2DsrZRERHMHXfVDqW7pjgAzJSs0JZCzG89nBmHZzFtgvbrA5HSwFrz6yl74K+1ClYh9ldZyf7YSIpoW6hujT2a8w3W7/hQdQDq8NxuGPXj/H30b8ZUn0ImT31LtaeoqOj2b17NxMmTODQoUN8/vnnnDt3jjfffJPMmZ1jWyc3sUcq4z4UBWA+l117gr+P/s2NBzfSxLX1hLxb713yZcrHWyve0iODOZk9l/fQcVZHSuQowaJei1JVk+/IBiO5cu8KE3ZPsDoUh/tu63d4unnyes3XrQ7FaURGRrJ161bGjh3LvXv3mDFjBvv27aNHjx64uaXM+AxWSW5in2P2is8mIoOA1UD6+9Yl08Q9EymctTBNiza1OpSnkskjE182+ZKgi0HMOjjL6nA0Ozl96zStZrQim1c2lvddTo4MOawO6R8aFmlIw8IN+XrL1zyMTr2PyLS3y3cvM3X/VJ6v/Hyaa+FLjWwTesaMGVm7di1r1qyhUaNGTjuCX3I7z30HzAXmAaWAkUqpX1IiMGdx+tZpVp9ezcAqA1PtgDRJ0b9Sf/zz+TNi9Yh02TTqbK6HX6fl9JZExkSyou8KCmQpYHVICfqowUdcunuJGftnWB2Kw4zePpro2GiG1xludShpWkxMDDt27HiU0Ddt2sT8+fOpVKmS1aGluGRnGqXUKqXUO0qp4UqpVSkRlDOZvGcyLuLC81WetzqUZ+IiLvzY4kcu3LnAd1u/szoc7RmER4XT7s92XLhzgcW9Fqfqx6U29mtM2VxlGb97vNWhOMTD6IdM2D2BzmU6UzR7UavDSbNOnz7N77//zp07d9iwYQPz58+nfPnyVoflMElK7CKy2fx9V0Tu2PzcFZE7KRti2hUdG83ve3+nVfFWqfaMKDkaFG5A17JdGbVlFCF3QqwOR3sK0bHR9JrXi6CQIGZ2nkndQnWtDumxRITB/oPZcXEHe6/stTqcFLfo2CJuPbzFYP/BVoeSJkVGRrJkyRJWr17N6NGj2bhxIxUrVrQ6LIdLUmJXStUzf2dWSmWx+cmslEobQ/FYYNmJZVy6eynNdppLyLfNviUmNoYRq/Xtb2mNUorXAl9j0bFF/NLqFzqV6WR1SEnSr1I/PF09mbDL+bvzTN4zmUJZC9HYr7HVoaQ5V69eZcqUKZQsWZLjx4/TpUsXp72G/iTJvY89RwI/1t8bk0pN2jOJvJny0qZEG6tDsZsi2YrwTp13mHFgBlsvbLU6HC0Zvtr8Ff/b9T/erfsuQ2oMsTqcJMuRIQfdynVj+oHpTj288YXbF1h5aiUDKg2wfByBtObSpUvMnDmTzz//nBkzZpAxY/q+YSu519h3A6HAcYwR50KBsyKyW0Sq2ju4tCzsYRiBJwLpXb53qrgv2J5G1BuBb2Zf3lj+hr79LY2Yum8qH6z9gD4V+vBlky+tDifZXqr6Enci7jD70GyrQ0kxU/dNRaF4rvJzVoeSpoSFhTF37lwmT57M88+n7b5M9pLcxL4KaK2UyqmU8gFaAUuAV4Gx9g4uLVt4dCFRsVH0KN/D6lDsLqNHRr5u+jU7L+1k6r6pVoejPcGqU6sYuGggTfyaMLnD5DR5d0bdgnUpk7MM43c5Zye6WBXL5L2TaVSkke40l0wbN25k6NChdO7c2epQUo3k3pVfSyn1aHxDpdRKEflOKfWS+aQ3zTT70GyKZCtC9fzVrQ4lRfSu0JsxwWN4d/W7tC/VPtXdA+0oYQ/CCLoUxP4r+zl24xjnws5xO+I29yLvER4VToyKwQUXRAQ3FzeyembFx9uHXBlzUcanDHUL1aV2gdp4uXulSHx7Lu+h85zOlM1Vlnnd56X6pwomRkQYXHUwb614i31X9lEpr3PdsrTp3CZO3zrNxw0/tjqUNCUyMpJjx46xfPlyq0NJVZKb2C+LyLtA3CglPYCrIuIK6DZZ080HN1l1ehXDag1z2s4bIsJvbX6j6viqvLPyHSZ1mGR1SCnu4p2LzD08l8ATgRy5foSr968SGROZ4LyCkchFBGOwRuOsLEbFJDi/p6snhbIWooZvDTqW6kjHMh1xc3m20bDOhp2l9czWZPfKzrI+y8jqlfWZ6rNav4r9GLF6BBN2T+DX1r9aHY5d/b73dzJ7ZKZL2S5Wh5Km3L9/nwwZMpArVy6rQ0lVkrvn6A38F/gbY1jZLeY0V4xHumrAgiMLiI6Npns5594klfJW4p067zBqyyj6VOzjdD15b4TfYNq+acw6NIv9V/fzIPr/B+ZxERd8MvhQLHsxyuUuR9lcZamUpxJV8lYhe4bsjz2gC3sYxokbJwi6GMSuy7s4GnqUU7dOcfLmSU7cPMGMAzNwEReKZi9Ku5LtGFZ7WLJvl7wRfoOW01vyMPoha15YQ/7M+Z96O6QWPt4+dCzdkdmHZjO65ehnPvBJLe5G3OWvw3/Rt0LfVDWkb1qQLVs2IiIiuHjxYpp7tGpKSvI3wzwr/0kp1SeRWU7aJ6S0b87hORTNXhT/fP5Wh5LiRjYcydwjcxm8eDAHXjlABvcMVof0TI5eP8oP235g8bHFXLl/5dH0jO4ZqZG/Bs2KNaNHuR6Uy13uqa9VZ/PKRnXf6lT3/edlmtjYWNadXcesg7NYeXolp26e4sftP/Lj9h8pmKUgz1d+nnfrvfvEnf+DqAe0n9Wes2FnWdVvFWVzlX2qOFOjHuV6MPvQbNadWUezYs2sDsculp5YSnhUOH0r9rU6lDRHRKhYsSKvvvoqf//9t9O2kCZXkvdMSqkYoLCIpM2LdA5yPfw6a06voUe5HuniQ5bBPQPj247n1K1TfLLhE6vDeSrHrh+j3/x+ZP86O2XGlGHC7glcu3+NUj6leKfOOxwfepx7798jaFAQnzf+nAp5KqRIBzQXFxeaFG3ChPYTOPfmOcI/CGdMqzFUyF2BkDshfLrxUzJ9mYkGvzdg35V9CdYRHRtNz3k92XZhG9M7T6d+4fp2j9NKLYu3JJNHJuYcmmN1KHYz78g88mbKm+oHC0qtAgIC2Lt3L1988cWjy17pXXLbsk4DW0RkEeajWwGUUj88qaCItAR+wmi2n6iUGhXvfU9gKlAVuAH0UEqdFREfjPHpqwNTlFJDbcqsB/IBcW2kzZVS15K5TnY1/8h8YlSM0zfD22rk14gXKr/Ad1u/o2vZrlTLX82hy78RfoMlx5dw4NoBDoce5uj1o0THRpM7Y25yZ8xNoayF6Fi6I038mjy69fDEjROM2jyKhccWcuPBDQDcXNyo6VuTQf6D6FOxD15uKdOhLam83Lx4tcarvFrjVe5H3mfU5lFM2D2BTec3UXlcZYpmL8q3zb6lcxmjN7BSiiFLhzwagKZr2a6Wxp8SMrhnoH2p9sw/Op+xbcam+VtJw6PCCTwRyIBKA9Lk3Qqpgbu7O507d2bcuHEcO3aMCRMm4OVl7XfXaslN7KfMHxeS8Rx2sxl/DNAMCAGCRWSRUuqwzWwDgVtKqeIi0hP4GqNz3kPgI6C8+RNfH6XUzmSuR4qZc2gOJX1KUimPc/XafZLvmn/HqtOr6DKnC7sG7yKnd84UXV5kTCSBJwKZum8qS44vISo2Ci83L0r5lKJmgZp4unoSGh7KtfvX2HJhC+N2jSNHhhwUzlqYs2FnufXwFgCu4koN3xq8XuN1epbvmWoHBsnokZHPGn/GZ40/Y/Wp1QxfNZx9V/fRZU4XCmUpRGfPzpyJPMPC8wt5v977DK0x9MmVplHdy3Zn5oGZrD2zlhbFW1gdzjNZfnI54VHhdCmjO809i2zZstG/f38CAwOpUaMGf/zxB1WqVLE6LMskK7ErpT4BEBFvpVR4MorWAE4qpU6b5WcBHQDbxN4B+Nj8ey7wq4iIUuo+sFlEiicnVitcu3+NdWfX8X6999NFM7yt7BmyM7/HfOpNrkfPuT1Z3nd5inRuioyJZMreKXyx6QvO3z5Pnox5eL3m6/Sp0IeKeSr+KzFHxUQxevtofgr6iYt3L3LzwU0ACmQpwEf1P2Kg/8BUm8wT07RYU/YW28upm6fot6Af20K2MZrRAHQs1ZHPG39ubYAprEXxFmTxzMKcQ3PSfGKfd2QePhl8aFikodWhpHkeHh506NCBPXv20LhxY3r27MlXX31FtmzZrA7N4ZI7pGxtETkMHDVfVxKRpAxM4wtcsHkdYk5LcB6lVDRwG/BJQt2/i8heEflIEsmmIjJYRHaKyM7Q0NAkVPl0Fh9bTKyKdcom0KSolr8av7X5jTVn1vD+mvftWnd0bDQTd0+k5C8leWnJS+TLlI9FPRcRMiyE75p/R5V8VR4l6NjYWMbvGk+V/1XB6wsv/rP6P1y8exG/bH68V+892pdsT8idEEauH8mkPZPS7HW5YjmKsXXgVn5u+fOjaX8f+5uOszpyL/KehZGlLC83LzqU6sD8o/MTvd0wLYiIjmDxscV0LP3stzZqBhHB39+fwYMHs3v3booXL873339PeHhyzkPTvuRe1BkNtMC4Bo5Sah/QwM4xJUcfpVQFoL750y+hmZRS45VS1ZRS1VLyfsclJ5ZQKGshKuZJf08TivN8led5pdorfLv1W2YffPbhP5VSLDiygPJjyzNo8SDyZMrDsj7L2DZwG+1KtXu0Q4yKieLXHb/iP84fj889eGnJS+y9upc8GfPwRs03CH0nlNNvnObLJl+ysNdCtr6wlVI5S/HSkpfo9lc3bj+8/cyxWmHtmbUMXzWcYh7FGJprKLm8c7Ho+CJ8vvHhl6BfrA4vxXQv152wh2GsPr3a6lCe2qrTq7gbeVc3w6cAb29vWrduTbdu3Zg6dSqFCxfmp59+4uHDh1aH5hBP8zz2C/EmJTzixj9dBAravC5gTktwHhFxA7JiHkA8JpaL5u+7wEyMJn9LPIx+yKpTq2hbom26a4aPb3TL0dQtWJe+C/ryW/BvT1WHUoq1Z9ZS7/d6dJ7TGRHh7x5/s33gdloWb4mIEHInhHdWvkPJX0ri9YUXry17jT1X9pDTOyeD/Adx7s1zXHr7EqNbjv7XNf/aBWuzfsB6vm32LX8f/ZtqE6qluceCBl8MpsOsDpT0Kcmbud+kqndVrr1zjf82+C+xKpbXl79Oxd8qcuXelSdXlsY0L9acrJ5Z03Tv+HlH5pHVMytNijaxOhSnlTdvXjp37kyXLl2YOHEiBQsW5PPPP+fmzZtWh5aikpvYL4hIHUCJiLuIDAeOJKFcMFBCRPzM2+V6AovizbMIGGD+3RVYqx7TRioibiKS0/zbHWgLHEze6tjPhrMbuB91n7Yl21oVQqrh4erB0t5LaVGsBa8GvsorS15JcpOpUorAE4HUnVyXJlObcObWGca1HceBVw5Qv1B9ftj2Aw1+b0DWUVkp+GNBvtv2HSdunsA3sy+v13idi8MucmX4Fca3G0+hrIUeuywRYXid4Wx4bgMPoh5Qa2ItZh2c9dgyqcXBawdpNaMVubxzsaLvCjK5Znr03seNPubiWxepnLcyB64doOCPBfl+6/cWRmt/Hq4edCrTib+P/k1EdITV4SRbVEwUC48upF2pdml2mN+0JF++fHTt2pVu3bqxcOFCihQpwpAhQzhz5ozVoaWI5Cb2l4EhGNfDLwKVzdePZV4zHwqswDgQmKOUOiQin4pIe3O2SYCPiJwEhgGPHvgtImeBH4DnRCRERMoCnsAKEdkP7DXjseyBzUuOL8Hb3ZtGfo2sCiFVyeqVlYU9F/Ju3Xf5367/0WxaM1adWpXgtV+lFAeuHuC/6/5LubHlaDOzDWfCztC1TFfqFazH91u/J8tXWfD51ofhq4az6fwmomOjqelbk59b/kz4++Gcf+s8P7X66alGWKtbqC57XtpDzQI16TWvF19v/jpVX3c/fuM4Tac2xcPVg1X9ViW4zrkz5WbPS3sY23osgjB81XDqTKpDeKTzXGvsWqYrtyNus+HcBqtDSbZ1Z9dx6+EtupZJn/1xrJInTx7atm3LoEGDOHjwIJUrV6Z169asXbs2VX/nkyu5veKvA4mNPPeksoFAYLxpI23+fgh0S6RskUSqTRWPilVKseTEEpoWbWr5vc+piauLK581+oxMHpn4dMOnNJ/eHEHI7JIZD/FgxHcjeBD9gAdRD4iKjQKMMdYBrty7wtwjcx9Ny5EhB5XzVqZ5sea8UOWFJ56NJ1eujLlY2Xclzy18jhFrRnA27Cy/tP4l1XVqOn3rNI3/aEysimXdgHUUy1HssfO/Uv0VupTpQoMpDdgWso083+chsHegUwxc09ivMRncMrD0+FKaF2tudTjJsvjYYjK4ZUhzcTuLLFmy0LhxY+rVq8f+/fvp378/3t7eDBs2jH79+qX557kna68lIrmAQUAR27JKqRfsG1bacjj0MGfDzvJ+Pfv2BE+Ldl/azeS9k1l7Zi3nbp8jPOqfZ4gKxZ3YO8aL+/8/XRAyuGUgq1dWcmfMTSmfUlTOW5mGRRpSq0Athwze4enmyYzOMyiStQijtowi5G4Is7vOTjXjd5+/fZ4mU5vwIPoB6waso0yuMkkqlztTbo4OPcoby97g5x0/02BKA0Y2GMknjdLmSIFxMrhnoLFfY5acWMLolqPTVN+WZSeX0civUZofgjmt8/DwoFq1alStWpUzZ84wduxY3n33Xfr168frr79OyZIlrQ7xqST3dGQhsAlYTdI6zaULS44vAaB1idYWR+J4sbGxzDk0h593/MzOSzsfnXkDZPfKTrHsxSjpU5IqeauQO2NuMntmJmhjEBlcMjC4x2CyZchGBrcMeLqljqf+uogLXzX9ikJZCzEkcAjNpzVnca/FZM+Q3dK44s7Ubz28xZr+a57qzoufWv1ElzJdaDmjJZ9u/JRN5zexst/KVNcqkRxtSrRh6YmlHLtxjNI5S1sdTpKcuHGCU7dO8Vatt6wORTOJCEWLFqVo0aKEhYWxe/duatSoQeXKlXnrrbdo27Ytrq5pZ7yL5H6jvZVS76ZIJGnYkhNL8M/nj2+W9PN0of1X9/PW8rfYdH7To2SeI0MO6hasS/dy3elWtluiyfrhTuOWkyLZizgq3GR7pfor5MqYiz7z+1D/9/qs6LvCsv/vsevHHp2pr+2/lqr5n/4KVIMiDQh5K4RqE6qx7uw6Cv5YkOBBwcl+elxq0aZkGwg0Dq7TSmIPPGFckWxVopXFkWgJyZYtG40bN6Z+/focPnyYYcOG8corrzBkyBAGDx6cJh4Rm9z2zSUikv5OSx/jRvgNtl7YStsSzt8bXinFuJ3jKDK6CJX+V4m1Z9eS0SMjfSv05eiQo9z4zw0W9VpE34p9U80Z+LPoWrYry/os49ztc9SZXIdD1w45PIZD1w7RcEpDImMiWTdg3TMl9Tg5vHNw8rWTdCrdiSv3rlD85+Jsu7DNDtE6XqGshaiQuwJLTyy1OpQkW3ZyGSV9SlI0e1GrQ9Eew93dnUqVKtG/f3/atWvHwoUL8fPzo2fPngQFBVkd3mMlN7G/gZHcH4rIHRG5KyJ3UiKwtGL5yeXEqlinv83t1x2/kuObHLy89GXO3T5H5TyVWd5nObfevcW0ztMolbOU1SGmiMZ+jdnw3AYiYyKpO7kua06vcdiy159dT/3f6+Pq4srG5zfadeAjFxcX5veYz1dNviIiJoJ6v9djxoEZdqvfkdqWbMumc5sIexhmdShPFB4Vzvqz62ldXJ8fpSX58+enTZs2DBkyhFu3btGuXTsqV67MtGnTiIhIfbdbJiuxK6UyK6VclFJeSqks5ussKRVcWrDkxBLyZMxjlzOp1GjK3in4fO3Da8te407EHTqU6sDlYZfZ8/KeND9Od1L55/Nn+8DtFMhSgJYzWjJl75QUX+Yfe/+g+bTm5M2Ul83Pb06xZuYR9UYwv/t8BKHv/L6MXDfyyYVSmTYl2hCjYlh5aqXVoTzRujPriIiJ0M3waZS3tzd16tThlVdeoUyZMnz11Vf4+vry4YcfcvXqVavDeyS5Y8WLiPQVkY/M1wVFxLLR3qwWq2JZeWolLYu3dLpHLh68epASv5Tg+YXPc+vhLdqVbMfV4Vf5u+ff5M2c1+rwHK5wtsJseWELAUUCeH7h87y94u0UGac8VsXy4doPeW7hczQo3ICtA7fil93P7sux1alMJ3YO2om3uzefbfyMQYsGpejy7K1WgVrkyJAjTTTHLzu5DG93bxoUtnIkbu1Zubi4ULp0aXr06EHPnj1ZvXo1xYsXp2/fvuzbt8/q8JLdFD8WqA30Nl/fw3gca7q098pebj64SbOizawOxW4ioiPoNKsTFf5XgZM3T1Ijfw1ChoWwqNeiFH8Ua2qX1Ssrgb0DGVJ9CD9s/4F6k+tx+tZpu9UfdzvbF5u+4MUqL7KszzKyeWWzW/2PUzlfZU69foocGXIwcc9EOs3q5JDl2oOriyutirci8EQgMbGp92YdpRTLTi6jsV9jPd6FE8mVKxetWrXi1Vdf5erVqzRu3JiGDRtaOuhNchN7TaXUEIxnpKOUugWk2/EQV51aBUDTok0tjsQ+lp1YRs5vc/L3sb/xyeBDYO9AggYFPdVobs7K3dWdX1v/ytxuczl+4zhVxlVh5oGZz/QFVkoxff90KvxWgZ2XdjKp/STGtxuPu6u7HSN/sryZ8nLm9TP4Zvbl72N/U39yfWJjYx0aw9NqU6IN18OvE3wp2OpQEnX8xnFO3zpNq+K6Gd4ZeXt7U69ePV599VWyZ89O37598ff3Z8GCBQ7/HiU3sUeJiCug4NGANWnjm58CVp9ZTYXcFciTKY/VoTyTyJhI2v/ZntYzW3M/8j5Dqg8h9J1QfR3wMbqU7cLel/dSNldZ+szvQ53Jddh8fnOy69l2YRutZrSi34J+VMxTkf0v7+eFKi9YNthKFq8snH79NCV9SrL5wmaqT6ieJpJ7y+ItcRVXlh5Pvc3xy04uA9CJ3cm5ublRpUoVBg0aRKlSpXj77bcpW7YsCxcudNgZfHIT+8/AAiC3iHwBbAa+tHtUacCDqAdsOrcpzZ+t7768m9zf5mbx8cXkyZiHPS/t4dfWv6apUbysUiRbETY9v4mJ7SZy/vZ56v9enw6zOrDw6ELuR95PtFxEdASrTq2i2bRm1Jlch12Xd/F98+9ZP2B9il9PTwoPNw+OvHqEynkqs/vKbqpOqJrqk3v2DNmp7lud1WdS72Ncl59cTimfUqnif6ylPBcXF8qUKUP//v3x9/fntddeo0qVKqxatSrFl53cseJniMguoAkgQEelVFKe7uZ0tl7YSkRMRJpO7F9t+ooP1n6AQvFC5ReY2H6iTujJ5ObixkD/gfSq0IvR20fz7dZvWXRsEZ6unjT2a0ylPJXwdPPEy82Le5H32Hx+M9tDthMRE0GejHn4rtl3vFztZTJ6pK6xqV1cXNg1eBfVJ1Zn9+XdVB5Xmb0v7cXFJfV2Em3q15QvN3/J7Ye3yeqV1epw/iEyJpKN5zbyov+LVoeiOZiIUKpUKUqUKMHhw4fp168f1atXZ8yYMRQqZN9nXsRJ9liSSqmjwNEUiCVNWXV6Fe4u7mmyd2tkTCRN/mjC5gub8XT15K9uf9GuVDurw0rTvN29eb/++wyvM5zN5zez+NhilpxYwpozax71nncRF/zz+TOk+hAaFG5A82LNU/VY4S4uLgS/GEzNiTXZeXlnqk/uTYs25fNNn7Ph3Abal2r/5AIOFBQSxIPoBzT2a2x1KJpFXFxcKF++PKVKlWLbtm1UrFiRESNG8M4779h9uNq0O0i0xVafXk3tgrXJ5JHpyTOnIudun6Pa+GpcD79OyRwl2TJwS7rv7W5PHq4eNPZrTGO/xvzY8kfAuIUtMiYSQdLciHwuLi4EvRhErUm1CL4UTK1Jtdg+cHuqTO61CtTC292b1adXp7rEvvbMWgShYeGGVoeiWczd3Z0GDRpQvnx5Jk+ezKpVq5gzZw4+Pj52W0bq+3amATfCb7D78m6a+qWtZvjlJ5dT8ueSXA+/Tp8KfTg69KhO6g7gIi54uXmluaQex8XFhe0Dt1M+V3mCLwXTdFrq/Nx7unnSoHADVp9OfdfZ155di38+f8sfJqSlHjly5KBnz55ERUVRqVIlDh8+bLe6dWJ/CuvOrkOh0tT19VGbR9F6RmuiYqP4qcVPTO88XV9P15LMxcWFPS/voVj2Yqw7u452M1PnpZsmfk04cv0IF+9ctDqUR+5H3mfbhW26GV77F1dXV5o0aUK1atVo0aIF169ft0u9OrE/hVWnVpHFMwvVfatbHUqSvLjoRd5b8x6ebp5sfG4jr9d63eqQtDTIzcWNg68epECWAiw5sYQXFr5gdUj/EnewveaM48b0f5ItF7YQFRtFE78mVoeipVJVqlShaNGidOvWzS716cT+FFafWU2jIo1S/XOsY2NjaTa1GZP2TCKHVw6ODz1OvcL1rA5LS8O83Lw4MuQIPhl8+H3v76lubPmKeSqS0ztnqkrsa8+sxc3FjXqF9HdPS1zDhg3ZvXs3Fy5ceOa6dGJPptO3TnP61ulU3wwfGR1Jhd8qsPrMavyy+XHurXMUzFrQ6rA0J5DJIxP7X97/aGz5cTvHWR3SIy7iQhO/Jqw+vdqy4TzjW3tmLbUK1Ep1tzRqqYurqytFixZl3bp1z1yXTuzJtO6MsdFT8/Wy+5H3KfZLMQ5fP0wt31qceO1Emuu9r6Vu+bPkZ8eLO3B3ceeVpa+w6Ngiq0N6pGnRply6e4mj162/KzfsYRi7Lu+icZHUu7/QUo8bN25QsmTJZ65HJ/Zk2nh+I7m8c1EmZxmrQ0lQ2MMwiv1cjJA7IbQt0ZZtL27D1cW+90hqGkC53OVY3W81IkLn2Z05ePWg1SEBPLqWnRp6x284u4FYFUuTovr6uvZ4p06dIjw8nOrVn73vlk7sybTh7AYaFG6QKnuUX7t3jaI/FeXq/av0LNeTxb0XWx2S5uQaFGnA7x1+J0bFUGtSLa6H26dX77Pwy+5H0exFU8V19rVn1pLBLQM1fWtaHYqWil2+fJlFixYxf/58uwxWoxN7MpwLO8e52+dS5WhzV+5docSvJbj18BaD/AfxZ9c/rQ5JSyf6V+rPiLojuB91n0r/q0R0bLTVIdHErwnrzq6z/DGua8+upV6heml2DAMtZSml2Lt3L7NmzWLixIk0aGCf3KITezJsPLcRINWNHnU9/DplxpThTsQdXq/xOuPbjbc6JC2d+arpV3Qo1YFLdy9Re2Jtq8MhoEgAdyLusP/qfstiCL0fysFrB2lUpJFlMWip1507d5g3bx7Hjh1jw4YNdO3a1W5168SeDBvPbSSbVzbK5y5vdSiP3Ay/SclfShL2MIxXq7/KT61+sjokLZ2a330+5XOVZ+flnby4yNqHncS1qm04t8GyGOIe49uwSOo6EdCsdffuXVatWsWECRNo3749e/fupVKlSnZdRuq+ETuV2XBuA/UL1U81ndHuPLxDyV9LPmp+H9N6jNUhaelY3LjyBX4swKQ9k6jhW4PBVQdbEkuBLAUomr0oG89t5M1ab1oSw6bzm/By86Ja/mqWLN9e7kXe4/iN49x8cJObD25y5d4Vq0NKk27dukVwcDAHDhxgwIAB/PXXX+TLly9FluWwxC4iLYGfAFdgolJqVLz3PYGpQFXgBtBDKXVWRHyAuUB1YIpSaqhNmarAFCADEAi8oVLo5tXLdy9z4uYJy3ZU8UVER1BmbBluPLjBgEoDdPO7lip4e3gT9GIQZceW5ZWlr1AlbxXLRmhsULgBi48tRillSWfXjec2UtO3Jh6uHg5f9rO49eAWk/ZMYvHxxRwOPZwqOkSmVbGxsZw6dYp9+/Zx/vx5nnvuOebOnZtiCT2OQxK7iLgCY4BmQAgQLCKLlFK2o94PBG4ppYqLSE/ga6AH8BD4CChv/tj6DRgEBGEk9pbAspRYh03nNwGp4/p6bGwsFX+ryKW7l+hUuhNTOk6xOiRNe6SETwnmdptLx9kdCfgjgHNvnrPkYUMNCjVgyt4pHLl+hLK5yjp02Xcj7rLnyh7er/e+Q5f7tCJjIvl267dM2TuFUzdPoTDOjwQhd8bclMhRgpzeOcnlnYtc3rn4iq8sjjh1u3v3Lvv372f//v3kzJmTt956i549e5Ixo2MGKXLUGXsN4KRS6jSAiMwCOgC2ib0D8LH591zgVxERpdR9YLOIFLetUETyAVmUUtvN11OBjqRQYt9wdgOZPDJRJV+VlKg+yZRS1J5Um+M3j1OvYD3m95hvaTyalpAOpTvwYf0P+XzT59ScWJMTQ084/FGvcde2N5zd4PDEvi1kG7EqlvqF6zt0ucl1M/wmry17jb8O/0VUbBRgXMboXLozAyoNoFzucgn26NeJ/d+ioqI4duwYR44c4fz583Ts2JFPP/2UmjVrOrzFyFGJ3RewHQA3BIh/Y+ejeZRS0SJyG/ABEmsH8jXrsa3TN6EZRWQwMBigUKFCyY0dMAamqVuwruXjw7ed2ZYdl3ZQPld5NjxnXccgTXuSzxp/xubzm1l/bj295/dmVtdZDl2+XzY/fDP7svH8Rl6p/opDl73p3CZcxIXaBay/QyAh9yLv0Xd+XxYdW4RCkcEtAy9WeZFPGn1Croy5rA4vzYiNjeX8+fMcOXKEw4cPU6VKFd577z06deqEt7e3ZXGli85zSqnxwHiAatWqJfsa/PXw6xy8dpBe5XvZPbbkeC3wNQJPBlIwS0F2v7Tb4WdAmpZcq/qvwvd7X2Yfmk1jv8YO7aMiIjQo3IAN5zY4/Dr7pvObqJK3Cpk9MztsmUkRExvDR+s+4rut3xEVG0WODDn4qMFHvFHzjVQ56FZqdeXKFQ4dOsSRI0fw8fHhueee46+//qJgwdTxPA5HJfaLgO0aFzCnJTRPiIi4AVkxOtE9rs4CT6jTLuJuW7FyYJrR20fza/CvZPXMyoFXDuDu6m5ZLJqWVG4ubmx/cTslfy3JK0tfoVaBWlTMU9Fhy29YuCF/HvyTU7dOUTxH8ScXsIOI6AiCLgbxSjXHthI8SVBIEG3/bMv18Ot4uHrwccOPGdlwpE7oSRQWFsbBgwc5evQo0dHR9O3blzFjxlChQgWrQ/sXRyX2YKCEiPhhJN+eQO948ywCBgDbgK7A2sf1cFdKXRaROyJSC6PzXH/gl5QIfsPZDXi5eVE9vzW9excfW8xbK97Cw9WDPS/tIatXVkvi0LSn4Zfdjz+7/Em3v7rR4PcGXBl+BS83L4csO+5gfOO5jQ5L7Lsu7+Jh9EPqF0od19cjYyIZGjiUCbsnANCuZDtmdpmpHwyVBA8ePODQoUMcP36cq1ev0rlzZz766CPq1auXqltMHZLYzWvmQ4EVGLe7TVZKHRKRT4GdSqlFwCRgmoicBG5iJH8AROQskAXwEJGOQHOzR/2r/P/tbstIoY5zG89vpFaBWpYMC7n/6n46ze6Ei7iwpv8a/LL7OTwGTXtWXct25SX/lxi3exxNpzZl8wubHbLc0jlLk8s7FxvPbeSFKi84ZJmbzhl30KSG568fu36MRn804vK9y3i7ezOr6yzalWxndVipWkxMDCdPnuTw4cOcPHmSpk2b8s0339CyZUs8PdPG0MAOu8aulArEuCXNdtpIm78fAt0SKVskkek7+fctcHZ1L/Ie+67s471676XkYhJ068Et6k6qS4yKYVqnaaliR6FpT+t/7f7H2rNr2XJhC19s/IIPGnyQ4suMu84eNxy0I2w6v8k4oLC4E9rMAzMZ8PcAomOjqVuwLkt6LyGbVzZLY0rNrly5woEDBzh06BDFihVjyJAhdO/enaxZ014LaeptS0glgi8GE6NiqFOwjkOXGxMbg/94f+5F3eM/df5D34p9Hbp8TUsJWwduxdPVk4/WfUTwxWCHLLNB4QacCTvDhdsXnjzzM4qJjWHz+c2WNsNHxUQxaPEg+szvQ3RsNO/Xe59Nz2/SST0BkZGR7Nmzh6lTp7JgwQLq169PUFAQO3fuZNCgQWkyqUM66RX/LLZe2ApArQK1HLrcVjNacTbsLE39mvJ1s68dumxNSyk5vXMyv8d82sxsQ9NpTbk6/GqKX2+3vc7ep2KfFF3WwWsHuR1x27LEfjfiLi2mt2BbyDbcXdyZ03UOHct0tCSW1Cw0NJQ9e/Zw4MABatWqxejRo2nVqhVubs6REvUZ+xNsDdlK2VxlyZ4hu8OW+e6qd1l1ehWFsxZmed/lDluupjlC6xKteaXaK9yJuEPzac1TfHkVclcgs0fmRwfpKWnLhS2ANdfXL929RNXxVdkWso0cGXKw56U9OqnbUEpx+vRp5syZw6xZs2jYsCEHDhxg5cqVtGvXzmmSOugz9seKVbFsu7CNLmW6OGyZfx/9m2+2fkNG94zseWlPqnngjKbZ09g2Y1l1ehWbzm9i9PbRKfqgFlcXV2oVqMXWkJRP7FsvbCVfpnwUyVYkxZdl69C1QzT6oxGh4aEUyFyArQO3UjBr6rin2mqxsbEcPHiQnTt34uHhwbvvvkvfvn3x8nLMnRlW0Gfsj3H0+lFuPbxF3UJ1HbK8s7fO0v2v7riICxue2+DQVgJNc7RNz2/C3cWdt1e+zYkbJ1J0WXUK1mH/1f3cjbibosvZFrKN2gVrO/Te8D2X91B7Um1Cw0MpnbM0u1/arZM6RkLfv38/EyZMICQkhP/9738cO3aMF1980amTOujE/lhxTXeO6DgXFRNFzUk1iYqN4ueWP1M1f9UUX6amWSlvprz83vF3YlUsDaY0IDY2NsWWVbdgXWJVLEEXg1JsGVfvXeX0rdPUKeC4jrY7Lu6g/u/1uRt5F/98/mwfuN3y3vhWU0px5MgRJk6cyPnz55k6dSrbt2+nVatWqfrec3tKH2v5lLZe2IpPBh9K5CiR4stqMrUJ1+5fo3u57gypMSTFl6dpqUGfCn1oX7I9V+5doff8+GNW2U/NAjURhC3nt6TYMraFbAOgdkHHjA+/9cJWAqYEcD/qPlXyVmH9gPXpfvCqq1evMmvWLHbv3s2kSZPYsWMHTZs2TXej6+nE/hhbL2ylTsE6Kf6h+GjtR2w6v4niOYozq4tjH5ShaVab12MePhl8mH1oNoEnAp9c4Clk8cxChTwVHnVuSwlbL2zFw9UD/3z+KbaMODsv7aTZtGY8jH5I+dzlWTtgbaobl96RIiIiWLlyJbNnz+a1117j8OHDtGzZMt0l9Dg6sSfievh1jt04luLN8JvObeKLTV+QwS0DQS8GpdsPopZ+ubm4sbLfSgSh25xu3Iu8lyLLqVuwLttDthMTG5Mi9W8L2YZ/Pv8Uv33vSOgRmk5tyoOoB5TOWZr1A9an63vUT506xcSJE/Hz8+PEiRMMGTLEqXq4Pw2d2BOxPWQ7kLLX1+88vEPLGS1RKJb0WkKODDlSbFmalpr55/NneJ3hhEeH02J6ixRZRp2CdbgbeZeD1w7ave7ImEiCLwan+PX1s2FnafRHI+5G3MU3sy9rB6zFx9snRZeZWkVFRbFixQpWr17NtGnTmD59Ojly6H0o6MSeqC3nt+Dm4paiD36p93s9wqPC+U+d/9C4aOMUW46mpQXfNPuGEjlKsPXCVsYGj7V7/XULGne3pMT97Huv7CUiJiJFr69fu3+Nxn80JjQ8lMyemVkzYA15M+VNseWlZtevX+ePP/4gb968HDlyhBYtUuZgMK3SiT0RW0O24p/PnwzuGVKk/jeXv8mBawfwz+uvR5bTNNP6Aetxc3HjjeVvcOnOJbvWXSRbEfJmypsi19lT+g6aB1EPaDezHedun8PdxZ0VfVdQ0qdkiiwrtTt69CjTpk3jgw8+YN68eWl22NeUpBN7AqJiothxcUeKNautPbOWn4J+IpNHJjY+77iHU2haapc/S35+afkL0bHRNJnaxK51iwh1C9ZNkTP2bSHbKJS1EPkz57d73bEqln4L+rHj0g6UUvzV7S9qFqhp9+WkdkoptmzZwtq1a1m5ciUvv/yy7pOUCJ3YE7D3yl4eRj9MkaPve5H3aP9newRheZ/lZPTIaPdlaFpa9nL1l6mRvwZHbxzl0w2f2rXuOgXrcCbsDJfvXrZrvXF30KSE91a/x7wj8wD4vvn3tCuV/h67Ghsby/Lly7l8+TK7d++mRo0aVoeUqunEnoCUvB+18R+NuR91n7drv+2wEe00La1Z2W8lnq6efLLhE07dPGW3elPiOvuF2xcIuRNC7QL2319M2j2Jb7Z+A8DAygNTdOjd1Co6OpqFCxfi4eHB1q1b8fX1tTqkVE8n9gQEXQwif+b8FMhSwK71frHxC4IvBVM2Z1m+bf6tXevWNGeS1Ssrv3cwRqWzZ5N8lXxV8HLzsut19rgTAXufsW+9sJWXl76MIDQo1ICxbcemu6bn6Oho/v77b3x9fVm5ciWZM6ffe/WTQyf2BASFBFHT177XsPZf3c9H6z7C09WTTc9vsmvdmuaMelXoReMijTl3+xxvr3zbLnV6uHpQLX81u56xb7uwjQxuGaiUp5Ld6rx09xKdZndCKUXBLAWZ12MeHq4edqs/LYiJiWHRokUULFiQBQsWOP347vakE3s818Ovc+rWKbsm9pjYGJpMbYJCMbvrbHJ463stNS0plvZZire7Nz9u+5EjoUfsUmftArXZc2UPEdERdqlv+8XtVMtfDXdXd7vUFxEdQafZnbgefh1XF1cW9lpITu+cdqk7rVBKsWLFCnx8fFiwYAEeHunroOZZ6cQez46LOwCoVaCW3ers8VcProdfp3vZ7nQo3cFu9Wqas/Ny82J6p+kolN0GrqnpW5PImEj2Xtn7zHVFxkSy5/Ieu54IDA0cyo6LO4hVsUxsN5HKeSvbre60YtOmTURERLBo0SI8PT2tDifN0Yk9nqCQIFzExW5PV1twZAHzjs4jd8bc/NnlT7vUqWnpSacynWjq15QLdy4wfOXwZ64v7qA9bnTJZ7Hvyj4iYiLsdvvZH3v/YOKeiQAMrT6UfpX62aXetOTgwYMcPXqUFStWkClTJqvDSZN0Yo8n6GIQ5XOXJ5PHs3+gwh6G0WteLwRhbf+16eaRgZpmb4t7L8bb3Zsftv3wzE3yvll8KZClANsvPntij3sMrD3O2A9eO8jLS1/GRVyoXaA237f4/pnrTGuuXLnCqlWrWLp0KXny5LE6nDRLZxobcc9rtlezWqM/GhERE8HHDT+mXO5ydqlT09IjLzcvpnacarcm+VoFahEU8uzPZg+6GES+TPme+Q6ae5H36DqnK9Gx0WTzysacbnPSXWe5hw8fsmDBAsaOHUvlypWtDidN04ndxokbJwh7GGaXxD5q8yj2XtlLhdwVGBkw0g7RaVr61qVsF5r4NeHCnQu8t/q9Z6qrlm8tzoSd4eq9q89UT1BIkPGs92e4DU0pxStLX+H4jeNEx0Yzo/MMu99qmxasWLGCNm3a0Lt3b6tDSfN0YrcR16z2rB3nTt44yQdrP8DD1YN1A9bZIzRN04BFPRfh5ebFN1u/4cytM09dT9w18bjv/NO4+eAmJ26eeOYTgT/2/cH0/UYHwffqvUfL4i2fqb606ODBg9y+fZtffvnF6lCcgk7sNoJCgsjskZnSOUs/dR1KKRr90YhYFcvkDpPT7SMVNS0leHt4M7H9RGJV7DM1yfvn88fNxe2ZOtDF3UHzLIn9xI0TDAkcgou4ULdgXT5tZN8hdNOCe/fusXr1ambNmoW3t7fV4TgFndhtBF0MorpvdVxdXJ+6jpeWvETI3RBaFGtBnwp97BidpmkAfSr0oXaB2py4eYIvN335VHV4u3tTKU+lZ0rsQSFBCEK1/NWeqnxkTCS95vUiIjqCLJ5ZmNV1Fm4ubk8dT1q1Zs0aBg4cSPXqKfeI7PTGYYldRFqKyDEROSkiIxJ431NEZpvvB4lIEZv33jOnHxORFjbTz4rIARHZKyI7nyW+B1EP2Hd13zMdfW86t4kJuyeQxTMLi3otepZwNE17jMA+gXi4ejBy3cinfrxrrQK12HFxBzGxMU9VPuhiEOVylyOz59MNczpy3Uh2Xd5FjIrh9w6/p8vr6qdOnSI0NJRPP01/LRUpySGJXURcgTFAK6As0EtEysabbSBwSylVHPgR+NosWxboCZQDWgJjzfriNFJKVVZKPd1hs2n35d1Ex0Y/dWKPiI6g7Z9tAVjSe0m669GqaY6UzSsbv7T6hRgVQ6uZrZ6qjloFanE/6j6HQg8lu6xS6pnuoFl7Zi3fbDEe7vJy1ZfpWLrjU9WTlsXExLBmzRrGjBmjm+DtzFFn7DWAk0qp00qpSGAWEH8Itg7AH+bfc4EmYnQ17QDMUkpFKKXOACfN+uzq0f2oTznQRPs/23Mn4g6D/AdRv1B9e4amaVoCBlcdTMU8Fdl/dT/jdo5Ldvm4TrJPc9vbyZsnufng5lMl9rCHYfRb0A8XcaF0ztLp8n51gF27dlG8eHHatUt/j6FNaY5K7L7ABZvXIea0BOdRSkUDtwGfJ5RVwEoR2SUigxNbuIgMFpGdIrIzNDQ0wXmCLgZROGth8mbKm/S1Mk3dO5WVp1fim9mXcW2Tv4PRNO3prOi7Aldx5fXlrxP2MCxZZYtlL4ZPBp+nus4edyJQwzf55xhvLHuDy3cvIyLM7jobb/f0d7YaERHBtm3bGD16dLp7Yp0jpPXOc/WUUv4YTfxDRKRBQjMppcYrpaopparlypUrwYp2XNzxVF/SK3evMGjJIFzEhbX91+oPqaY5UN5Mefkk4BMiYyJpNzN5Z34iQs0CNZ9qBLqgkCC83b2TPfDU/CPzmbrfGGjnm6bfUDFPxWQv2xns3LmTxo0bU6VKFatDcUqOSuwXgYI2rwuY0xKcR0TcgKzAjceVVUrF/b4GLOApm+hD74dyNuxsshO7UoomU5sQGRPJxw0/pmTOkk+zeE3TnsEHDT6gaLaibL6wmbmH5yarbC3fWhwOPZzss/2gi0FUy18tWb3Yr9y7wouLXsRFXAgoEsAbtd5I1jKdRUREBDt27OCzzz6zOhSn5ajEHgyUEBE/EfHA6AwXv9v4ImCA+XdXYK1SSpnTe5q95v2AEsAOEckoIpkBRCQj0Bw4+FTBXQoGoHr+5N1u8cGaDzh8/TDlc5fnwwYfPs2iNU2zg+V9lyMIA/4eQGR0ZJLLxfWp2Xkp6TfVRERHJPsOGqUUgxcPJuxhGBndMzK141RcJK03mD6dPXv2EBAQQJkyZawOxWk55JNlXjMfCqwAjgBzlFKHRORTEWlvzjYJ8BGRk8AwYIRZ9hAwBzgMLAeGKKVigDzAZhHZB+wAliqllj9NfMEXgxEE/3z+SS4TdDGIUVtG4e7izqq+q3QTvKZZqIRPCYbWGEp4VDjd53ZPcrm4g/m4wWaSYv/V/UTGRCarhW/mgZksPr4YhWJc23EUzFrwyYWcUExMDDt37uSDDz6wOhSn5rDREJRSgUBgvGkjbf5+CHRLpOwXwBfxpp0GKtkjtuBLwZTNVTbJ96PejbhLmxltUCjGtB5D3szJ73CnaZp9jW4xmtkHZ7Pw2EI2n99MvUL1nlgme4bslPQpmazEHjdvUhP7lXtXeDXwVQC6l+tOrwq9krwsZ3P06FH8/Pz0YDQpLH22BdlQShF8KZjqvkn7oCml6D2vNzce3KBB4Qa86P9iCkeoaVpSuLi48HfPvwHoOKsjsbGxSSpXw7cGQReDMK78PVnwpWByZ8xNwSxPPuuOa4K/G3GXXN65+K3Nb0lahrPat28fb7/9ttVhOL10n9jP3z7PtfvXknx9/fe9v7PkxBK83b35q9tfugle01KR2gVr07VsV248uMGQwCFJKlMjfw2u3LvCxbvx+/MmbMfFHVTPXz1J3/1ZB2c9aoKf0nEKOTLkSNIynNGNGzcIDQ2lU6dOVofi9NJ9Yk9Ox7njN47zytJXAJjUfhK5M+ZO0dg0TUu+GZ1mkMkjE+N2jeNI6JEnzh/XpJ6UgWruRtzl6PWjSWqGv3b/Gi8vfRmAFyq/QOsSrZ9Yxpnt37+f/v374+GhR+VMaTqxXwzG3cX9ifeTRkRH0HFWRyJjImldojU9y/d0UISapiWHh5sH0zsZj0FtM7PNE+evnLcy7i7uSbrOvuvyLhQqSScCQ5YO4W7EXfJnys+PLX9MUuzOKjY2lkOHDvH8889bHUq6oBP7pWAq562Mp5vnY+d7d/W7HLl+hMwemfm9w+8Oik7TtKfRoXQH6heqz5mwM3yx8YvHzuvp5knlvJXZcenJiT0u+T+pT87iY4uZe2QuCsX0ztPJ4pkl6cE7oZCQELJnz06lSnbp76w9QbpO7LEqlp2Xdj7x6DvwRCA/Bf0EwLi243QTvKalAUt6GQ9j+njDx1y7d+2x89bwrcHOSzuf+KS34EvB+GXzI6d3zkTnuf3wNgMXDQSMB7w08muU/OCdzNGjR+nVK/3eDeBo6TqxH7t+jLuRdx979B16P5T+C/ojCG1LtNVN8JqWRmTxysKPLX4kOjaa1jMff327hm8N7kXe4+j1o4+dLylDT7+98m1Cw0PJnyk/3zb/NtlxOxulFMePH6dr165Wh5JupOvE/qSOc0opBi0exM0HN8nskZnx7cbrXvCaloa8Wv1VyuQsw67Lu5i+b3qi88WNIhf3cJeEXLt/jfO3zz+2hW/juY1M2jMJgBldjE586d3Vq1fx8PCgfPnyVoeSbqTvxH4xmEwemSids3SC7/+x7w8WHluIQvFL61/IlzmfgyPUNO1ZLe+7HBdxYfCSwTyMfpjgPCV8SpDVM+tjO9AFXzRPBBJp4YuIjqDfgn4ADPYfTECRgGcL3EmcOnWKtm3b6pMiB0rXiX3HpR1UzVcVVxfXf713LuwcrwW+hou40LJYS/pV7GdBhJqmPatCWQvxdu23eRD9gC6zuyQ4j4u4UN23+mMT+46LO3ARl0SHnv5s42ecv32e3Blzp9tnrCckJCSE1q3T961+jpZuE3tkTCR7r+xNsFktVsXy3N/P8SD6ARncMugmeE1L475p9g35MuUj8GQg68+uT3CeGvlrsP/qfh5EPUjw/bihpxNqXj96/SijNo8CYGbnmboJ3hQdHc3Zs2cJCAiwOpR0Jd0m9gNXDxAZE5lgs9rY4LGsP7eeGBXDDy1+SLcPbNA0Z7Kol/FAyS5zuiQ43GzNAjWJUTHsvrz7X+8ppYyOc/n/3XEuVsXSe15vYlQMfSv0pUnRJvYPPo26dOkSfn5+ZM2a1epQ0pV0m9gT6zh34fYFRqwegau4ElA4gEH+g6wIT9M0O6uWvxo9yvXg5oObj0aEsxW3L0ioA93ZsLPceHAjwROB8bvGs+fKHrJ7ZefX1r/aP/A0LCQkhAYNGlgdRrqTbhP7zks78cngQ5FsRR5NU0oxNHAoD6If4O7qzsT2E3UTvKY5kemdp5PZIzMTd0/813Cz+TLno2CWgo8O+m0ldiIQej+UYSuGGXV3mk5WL31mais0NJS6detaHUa6k64Te7X81f6RuOcfmc+i44uIVbF81eQriuUoZmGEmqbZm5uL26PhZhO6t726b/VHvd9tBV8MxsPVgwp5Kvxj+guLXuBB9APal2pP65K6g1h8ly9fxt8/4c6GWspJl4k9PCqcg9cOUi1/tUfTwh6G8Wrgq7iKKzV9a/JajdcsjFDTtJTSvnR76heqz9mws/8abrZ6/uqcunWKmw9u/mP6zss7qZy3Mh6u//8Ak9WnV7PkuPGkRz3M9L9FRkZy69YtSpUqZXUo6U66TOz7ruwjRsX8o1ltxOoRhN4PRUSY3GFygrfAaZrmHBb3WpzgcLNxo8rZnrXHqlh2XdpFtXz/fyIQER1B73m9ARjfbny6fhxrYq5du0bRokVxddX7UkdLl4l956WdAI/O2Hdd2sX4XeNRKP7b8L+UzVXWyvA0TUthWb2y8kPzH/413GzVfFUB/nGdPaGhp99d/S6h4aHUKVCH3uV7Oy7wNCQ0NJQyZcpYHUa6lC4Te/ClYPJmykv+zPmJVbG8svQVRITyucrzbt13rQ5P0zQHGFJjyKPhZqfumwoYCb90ztL/GKgm7kQgroXv+I3j/Bz0M+4u7sztPld3sE1EeHg4pUsnPKqnlrLSZWK37Tg3Y/8Mgi8Fo5Tij05/4O7qbnV4mqY5yLK+y3ARF15e8jLhkeGAkcDj9glgnAhkdM9I6ZylUUrRaXYnFIqvmn6lh5l+gmLFdAdkK6S7xB6rYjl6/SjV81fnbsRd3lzxJgDv1Hkn0aEiNU1zToWzFn403GzXv4ynj1XPX50r965wK+YWYCR2/3z+uLq4Mn7XeA6HHqZEjhIMqzXMytDThEKFClkdQrqU7hL7/aj7KBTV8lfjv+v/y80HNymYpSAfB3xsdWiaplkgbrjZZSeXseb0mkcd6E5HniZaRT8aevrOwzu8ueJNBGFJ7yW6CT4J8ufPb3UI6VK6S+xxzW3ZvbLzU9BPgDFoRQb3DFaGpWmaheKGm+32Vzcq5K6Am4sbZyLOcCnqEg+jH1Ldtzq95vXiYfRDXqvxGiV9SloccdqQN29eq0NIl9JdYr8fdZ9CWQsxYvUIYlUsAyoNoEFhPeShpqVn1fJXo1f5Xtx6eIvXl79OpTyVOBt5ljMRZwCIiY0h8GQgubxz8UOLHyyONvWLiooCIHv27BZHkj6lu8QeHhVOoSyF2Hh+I1k8s/Bzq5+tDknTtFRgaqepZPHMwuQ9kymWvRhnIs5wOvI02b2y8+byNwGY132eHuMiCW7eNAb40ZcrrJHuEntEdAT7ru4DYEqHKWTxzGJxRJqmpQZuLm7M6DQDhWLNmTU8UA/YFb6LLJ5ZuP7gOu1KtqN+4fpWh5kmhIWFWR1CuuawxC4iLUXkmIicFJERCbzvKSKzzfeDRKSIzXvvmdOPiUiLpNaZmLuRd6lbsC6dynR65vXSNM15tC3VloDCAdx4cAOA+7H3OXf7HBncMjC762yLo0s77ty5Y3UI6ZpDEruIuAJjgFZAWaCXiMQf3m0gcEspVRz4EfjaLFsW6AmUA1oCY0XENYl1JshVXJnfY/6zr5imaU5nYc+FuLv8czyL39r8pjvYJsO9e/esDiFdc9QZew3gpFLqtFIqEpgFdIg3TwfgD/PvuUATMS7QdABmKaUilFJngJNmfUmpM0Fts7Qld8bcREVFERAQwPTp0wFjpKSAgABmzzaOzG/fvk1AQADz5xsHAdevXycgIIDFixcDcOXKFQICAli+fDkAFy5cICAggNWrVwNw+vRpAgIC2LBhAwDHjh0jICCArVu3AnDw4EECAgIIDjaGr9y7dy8BAQHs3bsXgODgYAICAjh48CAAW7duJSAggGPHjgGwYcMGAgICOH36NACrV68mICCACxcuALB8+XICAgK4cuUKAIsXLyYgIIDr168DMH/+fAICArh9+zYAs2fPJiAggPBw486B6dOnExAQ8KgjzJQpUwgICHi0HSdMmEDTpk0fvR47diytWrV69Pqnn36iffv2j15/9913dOnS5dHrUaNG0bNnz0evP/vsM/r27fvo9ciRI3n++ecfvX7vvfcYPHjwo9fDhw9nyJAhj16/+eabvPnmm49eDxkyhOHDhz96PXjwYN57771Hr59//nlGjhz56HXfvn357LPPHr3u2bMno0aNevS6S5cufPfdd49et2/fnp9++unR61atWjF27NhHr5s2bcqECRMevQ4ICGDKlCkAdv/s3bhxg1GjRunPnulZPntZvLLQLKoZGOPTUCZnGY7OPqo/ezafvSeJ+z9q1pC40ZVSdCEiXYGWSqkXzdf9gJpKqaE28xw05wkxX58CagIfA9uVUtPN6ZOAZWaxx9ZpU/dgIO5bWR44aPeV1JIrJ3Dd6iCcjN6m9qe3acJKKaUy204QkVDgnPmyquNDchrngVDz78JKqVzJrcDNvvGkTkqp8cB4ABHZqZSq9oQiWgrT/wf709vU/vQ2TZiI7Iw/7WkSkJYyHNUUfxEoaPO6gDktwXlExA3ICtx4TNmk1KlpmqZp6YqjEnswUEJE/ETEA6Mz3KJ48ywCBph/dwXWKuM6wSKgp9lr3g8oAexIYp2apmmalq44pCleKRUtIkOBFYArMFkpdUhEPgV2KqUWAZOAaSJyEriJkagx55sDHAaigSFKqRiAhOpMQjjj7bx62tPR/wf709vU/vQ2TZjeLqmYQzrPaZqmaZrmGOlu5DlN0zRNc2Y6sWuapmmaE9GJXXMYEbkX7/VzIvKrVfE4CxGJEZG9InJQRP4SEW+rY9Kcj4goEfne5vVwEfnYwpC0ROjErmlp3wOlVGWlVHkgEnjZ6oA0pxQBdBaRnFYHAiAiU8yDDSUiUSJyTUTWicgQEXGPN29pEflTRK6KSISInBGR70Uke7z5/ERkuoiEmPNdEpGlIlLFsWv3bNJNYrc5q4n70Ts/zRltAopbHURap/cXCYrG6A3/ltWB2FgN5AOKAM2BxcAnwCYRyQggIjUwbo/ODHTEuGX6NYznjGwVkWzmfO7AKiAX0B0oCXTBuL06h4PWxy7SxchzpgdKqcpWB5HOZRCRvTavc6DHHrAbc2CnVsByq2NxAnp/kbAxwH4R+cbqQEwRSqm4wesvAntFZCWwG/iPealgMnAcaK+UijXnPS8iuzGePfIFMATjQWPFMIYqP2nOdw7Y5pA1saN0c8aupQpxTcaVzZ3myCcV0JIk7oBpJ8Y405OsDUdzVkqpO8BU4HWrY0mMUuogxsFtF6AyRsL+3iapx813CZiB8WRQwRifPRboYh4kp1k6sWta2md7wPSa+bRDTUspozEes53R4jge5zBQFKM5HeDIY+bLDuRSSl3EOGAZCYSJyAYR+UxEyqV4tHamE7umaZqWZEqpm8AcjOSeWgmPHrybdEqpMUBeoDewGeNR4HvNp4emGTqxa5qmacn1PcYjbVOrssBpjGvrca8Tm+8W//+YVJRSd5VSi5RSHwCVgHXAZykYq93pxK45jFIqU7zXU5RSQ62Kx1nE366alhJsP2dKqatKKW+l1McWhpQgESkPtATmAnsxmuHfFhGXePPlB/oAf6pExlY3px8F0tR3TCd2TdM0La3yFJG8IpJfRCqJyDBgPbAL+M5MzC9gXGtfKCK1RaSgiLTBuFXuHPAhgIhUFpGFItJVRMqKSHERGWiWX2DFyj0t/RAYTdM0Lc0RkSn8/6O+Y4Aw4CDGmfp4206kIlIWo1NcYyAbcAmYD3ymlLplzpMT+ABognFfvAvGXSZzgFFKqYcpvEp2oxO7pmmapjkR3RSvaZqmaU5EJ3ZN0zRNcyI6sWuapmmaE9GJXdM0TdOciE7smqZpmuZEdGLXNE3TNCeiE7umaZqmORGd2DVN0zTNiejErmmapmlO5P8AH/kRpKp9jhUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 504x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the band structure and DOS:\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rc(\"figure\", dpi=150)\n",
    "\n",
    "fig = plt.figure(1, figsize=(7, 4))\n",
    "bs_ax = fig.add_axes([0.12, 0.07, 0.67, 0.85])\n",
    "\n",
    "emax = 0.035\n",
    "bs.plot(ax=bs_ax, emin=0.0, emax=emax)\n",
    "\n",
    "dos_ax = fig.add_axes([0.8, 0.07, 0.17, 0.85])\n",
    "dos_ax.fill_between(dos.get_weights(), dos.get_energies(), y2=0, color=\"grey\", edgecolor=\"black\", lw=1)\n",
    "\n",
    "dos_ax.set_ylim(0, emax)\n",
    "dos_ax.set_yticks([])\n",
    "dos_ax.set_xticks([])\n",
    "dos_ax.set_xlabel(\"DOS\", fontsize=14)\n",
    "\n",
    "fig.suptitle(\n",
    "    f\"Phonon band structure and DOS of LiFePO4 {mp_129.formula} ({mp_id}) with {supercell} supercell\",\n",
    "    fontsize=12,\n",
    "    y=1.02,\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "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.9.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "95e5cb0dae99d6911ce19b14750141ca4ad955f9302739260772daefb06054c0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}