{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce2efe0c-f86c-4cbd-ab3f-2aae9c5a574d",
   "metadata": {},
   "source": "## Thermal Expansion \nCalculate the thermal expansion for a Morse Pair potential using the [LAMMPS](https://www.lammps.org/) molecular dynamics\nsimulation code. In the following three methods to calculate the thermal expansion are introduced and compared for a \nMorse Pair Potential for Aluminium. \n\nAs a first step the potential is defined for the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code \nby specifying the `pair_style` and `pair_coeff` commands for the [Morse Pair Potential](https://docs.lammps.org/pair_morse.html)\nas well as the Aluminium bulk structure: "
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a96bac5-5afe-4924-90e5-04f6e6b2bedb",
   "metadata": {
    "trusted": true
   },
   "outputs": [],
   "source": [
    "from ase.build import bulk\n",
    "import pandas\n",
    "\n",
    "potential_dataframe = pandas.DataFrame(\n",
    "    {\n",
    "        \"Config\": [\n",
    "            [\"pair_style morse/smooth/linear 9.0\", \"pair_coeff * * 0.5 1.8 2.95\"]\n",
    "        ],\n",
    "        \"Filename\": [[]],\n",
    "        \"Model\": [\"Morse\"],\n",
    "        \"Name\": [\"Morse\"],\n",
    "        \"Species\": [[\"Al\"]],\n",
    "    }\n",
    ")\n",
    "\n",
    "structure = bulk(\"Al\", cubic=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe980834-e174-464a-8c07-04db9889a8c6",
   "metadata": {},
   "source": "The `pandas.DataFrame` based format to specify interatomic potentials is the same `pylammpsmpi` uses to interface with \nthe [NIST database for interatomic potentials](https://www.ctcms.nist.gov/potentials). In comparison to just providing\nthe `pair_style` and `pair_coeff` commands, this extended format enables referencing specific files for the interatomic\npotentials `\"Filename\": [[]],` as well as the atomic species `\"Species\": [[\"Al\"]],` to enable consistency checks if the \ninteratomic potential implements all the interactions to simulate a given atomic structure. "
  },
  {
   "cell_type": "markdown",
   "id": "b825e555-5e3d-43c2-8c51-50207ead60b5",
   "metadata": {},
   "source": "Finally, the last step of the preparation before starting the actual calculation is optimizing the interatomic structure. \nWhile for the Morse potential used in this example this is not necessary, it is essential for extending this example to\nother interactomic potentials. For the structure optimization the `optimize_positions_and_volume()` function is imported\nand applied on the `ase.atoms.Atoms` bulk structure for Aluminium:"
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7bd4ed1a-bffe-40f9-99f8-706797418877",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "{'optimize_positions_and_volume': Atoms(symbols='Al4', pbc=True, cell=[4.05, 4.05, 4.05])}"
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from atomistics.workflows import optimize_positions_and_volume\n",
    "\n",
    "task_dict = optimize_positions_and_volume(structure=structure)\n",
    "task_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed380b51-efa8-4ebb-bb27-b4aca90b21ec",
   "metadata": {},
   "source": "It returns a `task_dict` with a single task, the optimization of the positions and the volume of the Aluminium structure.\nThis task is executed with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code using the \n`evaluate_with_lammps()` function:"
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "991eaf43-0c4b-4ca2-8494-2b11684f4a79",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": "[jupyter-pyiron-2datomistics-2djetdbyfr:00796] mca_base_component_repository_open: unable to open mca_btl_openib: librdmacm.so.1: cannot open shared object file: No such file or directory (ignored)\n/srv/conda/envs/notebook/lib/python3.10/site-packages/pylammpsmpi/wrapper/ase.py:165: UserWarning: Warning: setting upper trangular matrix might slow down the calculation\n  warnings.warn(\n"
    },
    {
     "data": {
      "text/plain": "Atoms(symbols='Al4', pbc=True, cell=[[4.047310585424964, 2.478262976797941e-16, 2.478262976797941e-16], [0.0, 4.047310585424964, 2.478262976797941e-16], [0.0, 0.0, 4.047310585424964]])"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from atomistics.calculators import evaluate_with_lammps\n",
    "\n",
    "result_dict = evaluate_with_lammps(\n",
    "    task_dict=task_dict,\n",
    "    potential_dataframe=potential_dataframe,\n",
    ")\n",
    "structure_opt = result_dict[\"structure_with_optimized_positions_and_volume\"]\n",
    "structure_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6750cc4c-106f-4066-80ad-860bf6980732",
   "metadata": {},
   "source": "The `result_dict` just contains a single element, the `ase.atoms.Atoms` structure object with optimized positions and \nvolume. After this step the preparation is completed and the three different approximations can be compared in the following."
  },
  {
   "cell_type": "markdown",
   "id": "6c120581-efd6-4204-8413-75ee81065db1",
   "metadata": {},
   "source": "### Equation of State \nThe first approximation to calculate the thermal expansion is based on the Equation of State derived by [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790).\nSo in analogy to the previous example of calculating the elastic properties from the Equation of State, the `EnergyVolumeCurveWorkflow`\nis initialized with the default parameters: "
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b69b6ee2-b526-4913-a6c7-36018e8960af",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "{'calc_energy': OrderedDict([(0.95,\n               Atoms(symbols='Al4', pbc=True, cell=[[3.9786988461213992, 2.43625040333692e-16, 2.43625040333692e-16], [0.0, 3.9786988461213992, 2.43625040333692e-16], [0.0, 0.0, 3.9786988461213992]])),\n              (0.96,\n               Atoms(symbols='Al4', pbc=True, cell=[[3.992610493736228, 2.4447688306981026e-16, 2.4447688306981026e-16], [0.0, 3.992610493736228, 2.4447688306981026e-16], [0.0, 0.0, 3.992610493736228]])),\n              (0.97,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.00642586504517, 2.4532283058243666e-16, 2.4532283058243666e-16], [0.0, 4.00642586504517, 2.4532283058243666e-16], [0.0, 0.0, 4.00642586504517]])),\n              (0.98,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.020146608667117, 2.461629838203636e-16, 2.461629838203636e-16], [0.0, 4.020146608667117, 2.461629838203636e-16], [0.0, 0.0, 4.020146608667117]])),\n              (0.99,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.033774328510742, 2.469974409946722e-16, 2.469974409946722e-16], [0.0, 4.033774328510742, 2.469974409946722e-16], [0.0, 0.0, 4.033774328510742]])),\n              (1.0,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.047310585424964, 2.478262976797941e-16, 2.478262976797941e-16], [0.0, 4.047310585424964, 2.478262976797941e-16], [0.0, 0.0, 4.047310585424964]])),\n              (1.01,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.060756898772644, 2.486496469098726e-16, 2.486496469098726e-16], [0.0, 4.060756898772644, 2.486496469098726e-16], [0.0, 0.0, 4.060756898772644]])),\n              (1.02,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.074114747931804, 2.494675792706855e-16, 2.494675792706855e-16], [0.0, 4.074114747931804, 2.494675792706855e-16], [0.0, 0.0, 4.074114747931804]])),\n              (1.03,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.087385573728375, 2.5028018298737613e-16, 2.5028018298737613e-16], [0.0, 4.087385573728375, 2.5028018298737613e-16], [0.0, 0.0, 4.087385573728375]])),\n              (1.04,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.100570779804249, 2.51087544008222e-16, 2.51087544008222e-16], [0.0, 4.100570779804249, 2.51087544008222e-16], [0.0, 0.0, 4.100570779804249]])),\n              (1.05,\n               Atoms(symbols='Al4', pbc=True, cell=[[4.113671733924125, 2.518897460846561e-16, 2.518897460846561e-16], [0.0, 4.113671733924125, 2.518897460846561e-16], [0.0, 0.0, 4.113671733924125]]))])}"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from atomistics.workflows import EnergyVolumeCurveWorkflow\n",
    "\n",
    "workflow_ev = EnergyVolumeCurveWorkflow(\n",
    "    structure=structure_opt.copy(),\n",
    "    num_points=11,\n",
    "    fit_type=\"birchmurnaghan\",\n",
    "    vol_range=0.05,\n",
    "    axes=[\"x\", \"y\", \"z\"],\n",
    "    strains=None,\n",
    ")\n",
    "structure_dict = workflow_ev.generate_structures()\n",
    "structure_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f5d1e8d-0204-4dca-9298-878b9b2f6406",
   "metadata": {},
   "source": "After the initialization the `generate_structures()` function is called to generate the atomistic structures which are\nthen in the second step evaluated with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code to derive\nthe equilibrium properties:"
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7d1f126e-4fd0-41c5-986b-91d3b5910e3e",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "{'energy': {0.95: -14.609207927145926,\n  0.96: -14.656740101454448,\n  0.97: -14.692359030099395,\n  0.98: -14.716883724875528,\n  0.99: -14.731079276327009,\n  1.0: -14.735659820057942,\n  1.01: -14.731295089579728,\n  1.02: -14.718611862249286,\n  1.03: -14.698196715842329,\n  1.04: -14.670598736769112,\n  1.05: -14.636332030744796}}"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_dict = evaluate_with_lammps(\n",
    "    task_dict=structure_dict, potential_dataframe=potential_dataframe\n",
    ")\n",
    "result_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb679eb1-338f-4485-a953-791e147fe632",
   "metadata": {},
   "source": "While in the previous example the fit of the energy volume curve was used directly, here the output of the fit, in\nparticular the derived equilibrium properties are the input for the Debye model as defined by [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790):"
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2e0f7aab-6744-4b6f-a454-38c28833a3ac",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "66.29787349319821"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "structure_opt.get_volume()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "11c8b18d-64ff-4c93-b646-668b00eb1cf8",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "{'b_prime_eq': 6.2365371733275845,\n 'bulkmodul_eq': 216.057292780608,\n 'volume_eq': 66.29790137569191,\n 'energy_eq': -14.735658078942949,\n 'fit_dict': {'fit_type': 'birchmurnaghan',\n  'least_square_error': array([8.12779273e-07, 2.83453476e-03, 1.45091623e-03, 3.00518393e-05])},\n 'energy': [-14.609207927145926,\n  -14.656740101454448,\n  -14.692359030099395,\n  -14.716883724875528,\n  -14.731079276327009,\n  -14.735659820057942,\n  -14.731295089579728,\n  -14.718611862249286,\n  -14.698196715842329,\n  -14.670598736769112,\n  -14.636332030744796],\n 'volume': [62.98297981853827,\n  63.645958553470244,\n  64.30893728840229,\n  64.97191602333424,\n  65.63489475826624,\n  66.29787349319821,\n  66.96085222813018,\n  67.62383096306218,\n  68.28680969799419,\n  68.94978843292616,\n  69.61276716785807]}"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit_dict = workflow_ev.analyse_structures(output_dict=result_dict)\n",
    "fit_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9c7cd51c-6058-4d1d-8948-56d29c3b13e7",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": "/srv/conda/envs/notebook/lib/python3.10/site-packages/atomistics/workflows/evcurve/debye.py:80: RuntimeWarning: overflow encountered in exp\n  return xi**3 / (np.exp(xi) - 1)\n"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "workflow_ev.analyse_structures(output_dict=result_dict)\n",
    "thermal_properties_dict = workflow_ev.get_thermal_properties(\n",
    "    temperatures=np.arange(1, 1500, 50),\n",
    "    output_keys=[\"temperatures\", \"volumes\"],\n",
    ")\n",
    "temperatures_ev, volume_ev = (\n",
    "    thermal_properties_dict[\"temperatures\"],\n",
    "    thermal_properties_dict[\"volumes\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35ab7b86-0688-4520-ad47-ea54b4bfde86",
   "metadata": {},
   "source": "The output of the Debye model provides the change of the temperature specific optimal volume `volume_ev`\nwhich can be plotted over the temperature `temperatures_ev` to determine the thermal expansion. "
  },
  {
   "cell_type": "markdown",
   "id": "88ccd1f0-98c5-4e13-ab2c-febe5d3f235b",
   "metadata": {},
   "source": "### Quasi-Harmonic Approximation \nWhile the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) approach based on the Einstein crystal\nis limited to a single frequency, the quasi-harmonic model includes the volume dependent free energy. Inside the \n`atomistics` package the harmonic and quasi-harmonic model are implemented based on an interface to the [Phonopy](https://phonopy.github.io/phonopy/)\nframework. Still the user interface is still structured in the same three steps of (1) generating structures, (2) evaluating \nthese structures and (3) fitting the corresponding model. Starting with the initialization of the `QuasiHarmonicWorkflow`\nwhich combines the `PhonopyWorkflow` with the `EnergyVolumeCurveWorkflow`:"
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "493663b9-ea0c-4234-87ef-8f70774794f4",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "{'calc_energy': {0.9: Atoms(symbols='Al108', pbc=True, cell=[[11.7229062192894, 7.178209789078681e-16, 7.178209789078681e-16], [0.0, 11.7229062192894, 7.178209789078681e-16], [0.0, 0.0, 11.7229062192894]]),\n  0.92: Atoms(symbols='Al108', pbc=True, cell=[[11.80910715486485, 7.230992638996672e-16, 7.230992638996672e-16], [0.0, 11.80910715486485, 7.230992638996672e-16], [0.0, 0.0, 11.80910715486485]]),\n  0.94: Atoms(symbols='Al108', pbc=True, cell=[[11.894067681419225, 7.283015957446018e-16, 7.283015957446018e-16], [0.0, 11.894067681419225, 7.283015957446018e-16], [0.0, 0.0, 11.894067681419225]]),\n  0.96: Atoms(symbols='Al108', pbc=True, cell=[[11.977831481208684, 7.334306492094308e-16, 7.334306492094308e-16], [0.0, 11.977831481208684, 7.334306492094308e-16], [0.0, 0.0, 11.977831481208684]]),\n  0.98: Atoms(symbols='Al108', pbc=True, cell=[[12.060439826001351, 7.384889514610908e-16, 7.384889514610908e-16], [0.0, 12.060439826001351, 7.384889514610908e-16], [0.0, 0.0, 12.060439826001351]]),\n  1.0: Atoms(symbols='Al108', pbc=True, cell=[[12.141931756274893, 7.434788930393824e-16, 7.434788930393824e-16], [0.0, 12.141931756274893, 7.434788930393824e-16], [0.0, 0.0, 12.141931756274893]]),\n  1.02: Atoms(symbols='Al108', pbc=True, cell=[[12.222344243795412, 7.484027378120565e-16, 7.484027378120565e-16], [0.0, 12.222344243795412, 7.484027378120565e-16], [0.0, 0.0, 12.222344243795412]]),\n  1.04: Atoms(symbols='Al108', pbc=True, cell=[[12.301712339412747, 7.53262632024666e-16, 7.53262632024666e-16], [0.0, 12.301712339412747, 7.53262632024666e-16], [0.0, 0.0, 12.301712339412747]]),\n  1.06: Atoms(symbols='Al108', pbc=True, cell=[[12.38006930767338, 7.580606125432298e-16, 7.580606125432298e-16], [0.0, 12.38006930767338, 7.580606125432298e-16], [0.0, 0.0, 12.38006930767338]]),\n  1.08: Atoms(symbols='Al108', pbc=True, cell=[[12.457446749652004, 7.627986143754963e-16, 7.627986143754963e-16], [0.0, 12.457446749652004, 7.627986143754963e-16], [0.0, 0.0, 12.457446749652004]]),\n  1.1: Atoms(symbols='Al108', pbc=True, cell=[[12.533874715230777, 7.674784775460657e-16, 7.674784775460657e-16], [0.0, 12.533874715230777, 7.674784775460657e-16], [0.0, 0.0, 12.533874715230777]])},\n 'calc_forces': {(0.9,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[11.7229062192894, 7.178209789078681e-16, 7.178209789078681e-16], [0.0, 11.7229062192894, 7.178209789078681e-16], [0.0, 0.0, 11.7229062192894]]),\n  (0.92,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[11.80910715486485, 7.230992638996672e-16, 7.230992638996672e-16], [0.0, 11.80910715486485, 7.230992638996672e-16], [0.0, 0.0, 11.80910715486485]]),\n  (0.94,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[11.894067681419225, 7.283015957446018e-16, 7.283015957446018e-16], [0.0, 11.894067681419225, 7.283015957446018e-16], [0.0, 0.0, 11.894067681419225]]),\n  (0.96,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[11.977831481208684, 7.334306492094308e-16, 7.334306492094308e-16], [0.0, 11.977831481208684, 7.334306492094308e-16], [0.0, 0.0, 11.977831481208684]]),\n  (0.98,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.060439826001351, 7.384889514610908e-16, 7.384889514610908e-16], [0.0, 12.060439826001351, 7.384889514610908e-16], [0.0, 0.0, 12.060439826001351]]),\n  (1.0,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.141931756274893, 7.434788930393824e-16, 7.434788930393824e-16], [0.0, 12.141931756274893, 7.434788930393824e-16], [0.0, 0.0, 12.141931756274893]]),\n  (1.02,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.222344243795412, 7.484027378120565e-16, 7.484027378120565e-16], [0.0, 12.222344243795412, 7.484027378120565e-16], [0.0, 0.0, 12.222344243795412]]),\n  (1.04,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.301712339412747, 7.53262632024666e-16, 7.53262632024666e-16], [0.0, 12.301712339412747, 7.53262632024666e-16], [0.0, 0.0, 12.301712339412747]]),\n  (1.06,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.38006930767338, 7.580606125432298e-16, 7.580606125432298e-16], [0.0, 12.38006930767338, 7.580606125432298e-16], [0.0, 0.0, 12.38006930767338]]),\n  (1.08,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.457446749652004, 7.627986143754963e-16, 7.627986143754963e-16], [0.0, 12.457446749652004, 7.627986143754963e-16], [0.0, 0.0, 12.457446749652004]]),\n  (1.1,\n   0): Atoms(symbols='Al108', pbc=True, cell=[[12.533874715230777, 7.674784775460657e-16, 7.674784775460657e-16], [0.0, 12.533874715230777, 7.674784775460657e-16], [0.0, 0.0, 12.533874715230777]])}}"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from atomistics.workflows import QuasiHarmonicWorkflow\n",
    "from phonopy.units import VaspToTHz\n",
    "\n",
    "workflow_qh = QuasiHarmonicWorkflow(\n",
    "    structure=structure_opt.copy(),\n",
    "    num_points=11,\n",
    "    vol_range=0.10,\n",
    "    # fit_type='birchmurnaghan',\n",
    "    interaction_range=10,\n",
    "    factor=VaspToTHz,\n",
    "    displacement=0.01,\n",
    "    dos_mesh=20,\n",
    "    primitive_matrix=None,\n",
    "    number_of_snapshots=None,\n",
    ")\n",
    "structure_dict = workflow_qh.generate_structures()\n",
    "structure_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dcd4a1e-7122-4f57-93c1-bd9267084f70",
   "metadata": {},
   "source": "In contrast to the previous workflows which only used the `calc_energy` function of the simulation codes the `PhonopyWorkflow`\nand correspondingly also the `QuasiHarmonicWorkflow` require the calculation of the forces `calc_forces` in addition to\nthe calculation of the energy. Still the general steps of the workflow remain the same: "
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2e96e588-e279-4d6c-8f40-2eafa982933b",
   "metadata": {
    "trusted": true
   },
   "outputs": [],
   "source": [
    "result_dict = evaluate_with_lammps(\n",
    "    task_dict=structure_dict,\n",
    "    potential_dataframe=potential_dataframe,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fa40f79-f919-47df-ab07-c9a9dcd04b3d",
   "metadata": {},
   "source": "The `structure_dict` is evaluated with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code to \ncalculate the corresponding energies and forces. The output is not plotted here as the forces for the 108 atom cells \nresult in 3x108 outputs per cell. Still the structure of the `result_dict` again follows the labels of the `structure_dict`\nas explained before. Finally, in the third step the individual free energy curves at the different temperatures are \nfitted to determine the equilibrium volume at the given temperature using the `analyse_structures()` \nand `get_thermal_properties()` functions:"
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "371977fd-cd4e-469f-8955-17a2946c8629",
   "metadata": {
    "trusted": true
   },
   "outputs": [],
   "source": [
    "workflow_qh.analyse_structures(output_dict=result_dict)\n",
    "thermal_properties_dict_qm = workflow_qh.get_thermal_properties(\n",
    "    temperatures=np.arange(1, 1500, 50),\n",
    "    output_keys=[\"free_energy\", \"temperatures\", \"volumes\"],\n",
    "    quantum_mechanical=True,\n",
    ")\n",
    "temperatures_qh_qm, volume_qh_qm = (\n",
    "    thermal_properties_dict_qm[\"temperatures\"],\n",
    "    thermal_properties_dict_qm[\"volumes\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c7145b4-9a55-4212-a34f-a1300f7b440f",
   "metadata": {},
   "source": "The optimal volume at the different `temperatures` is stored in the `volume_qh_qm` in analogy to the previous section. Here the extension `_qm` indicates that the quantum-mechanical harmonic oszillator is used. "
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "70002fc3-2436-43ed-8a4c-0b3c1f9a3812",
   "metadata": {
    "trusted": true
   },
   "outputs": [],
   "source": [
    "thermal_properties_dict_cl = workflow_qh.get_thermal_properties(\n",
    "    temperatures=np.arange(1, 1500, 50),\n",
    "    output_keys=[\"free_energy\", \"temperatures\", \"volumes\"],\n",
    "    quantum_mechanical=False,\n",
    ")\n",
    "temperatures_qh_cl, volume_qh_cl = (\n",
    "    thermal_properties_dict_cl[\"temperatures\"],\n",
    "    thermal_properties_dict_cl[\"volumes\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb0db978-365f-43af-9a20-6ebb58fb8da9",
   "metadata": {},
   "source": "For the classical harmonic oszillator the resulting volumes are stored as `volume_qh_cl`. "
  },
  {
   "cell_type": "markdown",
   "id": "eb795fbd-0477-492a-b883-9cb31b58d3e2",
   "metadata": {},
   "source": "### Molecular Dynamics\nFinally, the third and most commonly used method to determine the volume expansion is using a molecular dynamics \ncalculation. While the `atomistics` package already includes a `LangevinWorkflow` at this point we use the [Nose-Hoover\nthermostat implemented in LAMMPS](https://docs.lammps.org/fix_nh.html) directly via the LAMMPS calculator interface. "
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a41d36c9-34eb-46e0-b713-57941dfb0296",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": "100%|██████████| 298/298 [06:54<00:00,  1.39s/it]\n"
    }
   ],
   "source": [
    "from atomistics.calculators import calc_molecular_dynamics_thermal_expansion_with_lammps\n",
    "\n",
    "structure_md = structure_opt.copy().repeat(11)\n",
    "result_dict = calc_molecular_dynamics_thermal_expansion_with_lammps(\n",
    "    structure=structure_md,  # atomistic structure\n",
    "    potential_dataframe=potential_dataframe,  # interatomic potential defined as pandas.DataFrame\n",
    "    Tstart=15,  # temperature to for initial velocity distribution\n",
    "    Tstop=1500,  # final temperature\n",
    "    Tstep=5,  # temperature step\n",
    "    Tdamp=0.1,  # temperature damping of the thermostat\n",
    "    run=100,  # number of MD steps for each temperature\n",
    "    thermo=100,  # print out from the thermostat\n",
    "    timestep=0.001,  # time step for molecular dynamics\n",
    "    Pstart=0.0,  # initial pressure\n",
    "    Pstop=0.0,  # final pressure\n",
    "    Pdamp=1.0,  # barostat damping\n",
    "    seed=4928459,  # random seed\n",
    "    dist=\"gaussian\",  # Gaussian velocity distribution\n",
    ")\n",
    "temperature_md_lst, volume_md_lst = result_dict[\"temperatures\"], result_dict[\"volumes\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2efeb52-ee54-4eb0-878a-184f353941bf",
   "metadata": {},
   "source": "The `calc_molecular_dynamics_thermal_expansion_with_lammps()` function defines a loop over a vector of temperatures in \n5K steps. For each step 100 molecular dynamics steps are executed before the temperature is again increased by 5K. For \n~280 steps with the Morse Pair Potential this takes approximately 5 minutes on a single core. These simulations can be \nfurther accelerated by adding the `cores` parameter. The increase in computational cost is on the one hand related to \nthe large number of force and energy calls and on the other hand to the size of the atomistic structure, as these \nsimulations are typically executed with >5000 atoms rather than the 4 or 108 atoms in the other approximations. The \nvolume for the individual temperatures is stored in the `volume_md_lst` list. "
  },
  {
   "cell_type": "markdown",
   "id": "eff137a4-61fc-4cbe-8c60-b0dc534a5f3f",
   "metadata": {},
   "source": "### Summary\nTo visually compare the thermal expansion predicted by the three different approximations, the [matplotlib](https://matplotlib.org)\nis used to plot the volume over the temperature:"
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "da8f641d-c5e6-4c10-8aeb-c891109e2e6d",
   "metadata": {
    "trusted": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": "Text(0, 0.5, 'Volume ($\\\\AA^3$)')"
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGwCAYAAACnyRH2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACywUlEQVR4nOzddVzU9x/A8dcdDRIGCCKhAiqKiqIzZgd2t5vtwjmdzp7tzNlus2M6e+bsmN2FgZ0YIEpKx31/f9y835gxVPBOeD8fj3vo9/v93Ofe31Pu3nxSpSiKghBCCCFENqfWdwBCCCGEEIZAkiIhhBBCCCQpEkIIIYQAJCkSQgghhAAkKRJCCCGEACQpEkIIIYQAJCkSQgghhADAWN8BfCw0Gg2PHz/G2toalUql73CEEEIIkQ6KovD8+XPy5cuHWv3mtiBJitLp8ePHuLi46DsMIYQQQryDBw8ekD9//jeWkaQonaytrQHtm2pjY6PnaIQQQgiRHtHR0bi4uOi+x99EkqJ0etFlZmNjI0mREEII8ZFJz9AXGWgthBBCCIEkRUIIIYQQgCRFQgghhBCAjCnKcKmpqSQnJ+s7DCGyPFNT0/+cXiuEEG9DkqIMoigKISEhREZG6jsUIbIFtVpNgQIFMDU11XcoQogsQpKiDPIiIXJwcMDS0lIWeBQiE71YTDU4OBhXV1f5eRNCZAhJijJAamqqLiHKnTu3vsMRIluwt7fn8ePHpKSkYGJiou9whBBZgHTIZ4AXY4gsLS31HIkQ2ceLbrPU1FQ9RyKEyCokKcpA0oQvxIcjP29CiIwmSZEQQgghBJIUCSGEEEIAkhSJDHbgwAFUKtUHXZpg1KhRlCpV6oO93sdIpVKxadMmfYchhBAGTZKibK5z586oVCq++uqrl6717NkTlUpF586dP3xgBqJatWqoVCpUKhVmZmY4OzvTqFEjNmzYoO/Q3kpwcDD16tXTdxhCiKzk1j5Ijtd3FBlKkiKBi4sLq1evJj7+//+5ExISWLVqFa6urnqM7MN50yrkPXr0IDg4mFu3brF+/Xq8vb1p27YtX3zxxQeM8P04OjpiZmam7zCEEFnFqQXwewv4oyukpug7mgwjSVEmURSFuOQ4vTwURXmrWEuXLo2rq2ua1o8NGzbg4uKCr69vmrKJiYn07t0bBwcHzM3N+fTTTzl9+vQb6z927BhVqlTBwsICFxcXevfuTWxsbJo6Bw4ciIuLC2ZmZnh6erJo0SIAli5dip2dXZr6Nm3a9MaZR6dPn6Z27drkyZMHW1tbqlatyrlz59KUUalUzJ07lyZNmmBlZcWPP/742vosLS1xdHTExcWF8uXLM2nSJObNm8eCBQvYu3cvADVq1KBXr15pnhcWFoaZmRl//fUXAO7u7owfP56uXbtibW2Nq6sr8+fPT/OcQYMG4eXlhaWlJQULFmT48OFpErYXXYWLFy/G1dWVHDly8PXXX5OamsrkyZNxdHTEwcGBcePGvXS//+w+e/jwIW3btiVXrlxYWVnh5+fHyZMnAbhw4QLVq1fH2toaGxsbypQpw5kzZ177/gghshFFgf3jYXt/QAEbZ8hCM0Fl8cZMEp8SzycrP9HLa59sfxJLk7dbM6lLly4sWbKEDh06ALB48WK6du3KgQMH0pQbOHAg69ev57fffsPNzY3Jkyfj7+/PrVu3yJUr10v1Xrp0CX9/f8aOHcuiRYt4+vQpvXr1olevXixZsgSAjh07cvz4cWbNmkXJkiW5e/cuz549e7ebB54/f06nTp2YNWsWAFOnTqV+/frcvHkTa2trXbmRI0cyYcIEpk+fjpGR0Vu9RqdOnfj+++/ZsGEDtWrVonv37vTq1YupU6fqWmRWrFhBvnz5qF69uu55U6dOZezYsQwdOpQ//viDr7/+mipVqlCkSBEArK2tWbp0Kfny5ePSpUv06NEDa2trBg4cqKvj9u3b7Nixg507d3L79m1atmzJ3bt38fLy4uDBgxw7doyuXbtSs2ZNypcv/1LsMTExVK1aFWdnZ7Zs2YKjoyPnzp1Do9EA0KFDB3x9fZkzZw5GRkYEBATI4ohCCNCkwvYBcEb7SyvVhkLVgZIUiazn888/Z8iQIdy7dw+VSsXRo0dZvXp1mqQoNjaWOXPmsHTpUt34lAULFrBnzx4WLVrEgAEDXqr3p59+on379nz33XcAeHp6MmvWLKpWrcqcOXMICgpi7dq17Nmzh1q1agFQsGDB97qXGjVqpDmeN28eOXPm5ODBgzRs2FB3vn379nTt2vWdXkOtVuPl5cW9e/cAaNGiBd9++y2bN2+mdevWACxZskQ3ZuuF+vXr07NnT0DbKjR9+nQOHDigS4qGDRumK+vu7s7333/PmjVr0iRFGo2GxYsXY21tjbe3N9WrV+f69ets374dtVpN4cKFmTRpEgcOHHhlUrRy5UqePn3K6dOndYmsh4eH7npQUBADBgzQxeTp6flO75EQIgtJSYQNX8CVTYAKGkyBst31HVWGk6Qok1gYW3Cy/Um9vfbbypMnDw0aNOC3335DURQaNGhAnjx50pS5ffs2ycnJVKpUSXfOxMSEcuXKcfXq1VfWe/bsWW7dusWKFSt05xRFQaPRcPfuXS5duoSRkRFVq1Z965hfJzQ0lBEjRvDXX3/x5MkTUlNTiYuLIygoKE05Pz+/93odRVF0CY+ZmRmfffYZixcvpnXr1gQEBHDhwoWXZnyVKFFC93eVSoWjoyOhoaG6c3/88QczZszg1q1bxMTEkJKSgo2NTZo63N3d07R45c2bFyMjozQ7xufNmzdNvf8UEBCAr6/vK1v2APr160f37t1Zvnw5tWrVolWrVhQqVCh9b4oQIutJfA6rO8Ddg6A2gebzoXhzfUeVKSQpyiQqleqtu7D0rWvXrrpxMb/88stL11+MVfr3eJ5/Jgf/ptFo+PLLL+ndu/dL11xdXbl169YbY1Kr1S+NkXrToGjQzqh7+vQpM2bMwM3NDTMzMypUqEBSUlKaclZWVm+s501SU1O5efMmZcuW1Z3r3r07pUqV4uHDhyxevJiaNWvi5uaW5nn/7oZSqVS6bqsTJ07Qtm1bRo8ejb+/P7a2tqxevZqpU6f+Zx1vqvffLCzenDSPGjWK9u3bs23bNnbs2MHIkSNZvXo1zZo1e+PzhBBZUOwz7YDq4AAwsYK2K6BQ9f982sdKBloLnbp165KUlERSUhL+/v4vXffw8MDU1JQjR47oziUnJ3PmzBmKFi36yjpLly5NYGAgHh4eLz1MTU3x8fFBo9Fw8ODBVz7f3t6e58+fpxmYHRAQ8Mb7OHz4ML1796Z+/foUK1YMMzOz9xqj9Cq//fYbERERtGjRQnfOx8cHPz8/FixYwMqVK9+6a+7o0aO4ubnxww8/4Ofnh6enJ/fv38/QuEHbWhUQEEB4ePhry3h5edG3b192795N8+bNdeO/hBDZSGQQLPbXJkSWuaHzn1k6IQJJisQ/GBkZcfXqVa5evfrKgcdWVlZ8/fXXDBgwgJ07d3LlyhV69OhBXFwc3bp1e2WdgwYN4vjx43zzzTcEBARw8+ZNtmzZwrfffgtou4I6depE165d2bRpE3fv3uXAgQOsXbsWgE8++QRLS0uGDh3KrVu3WLlyJUuXLn3jfXh4eLB8+XKuXr3KyZMn6dChw3+2jrxJXFwcISEhPHz4kJMnTzJo0CC++uorvv766zSDqEHbWjRx4kRSU1PfumXFw8ODoKAgVq9eze3bt5k1axYbN25857hfp127djg6OtK0aVOOHj3KnTt3WL9+PcePHyc+Pp5evXpx4MAB7t+/z9GjRzl9+vRrk14hRBYVehUW1YGwW2DrAl13gXMZfUeV6SQpEmnY2Ni8NIblnyZOnEiLFi34/PPPKV26NLdu3WLXrl3kzJnzleVLlCjBwYMHuXnzJpUrV8bX15fhw4fj5OSkKzNnzhxatmxJz549KVKkCD169NC1DOXKlYvff/+d7du34+Pjw6pVqxg1atQb72Hx4sVERETg6+vL559/rltC4F0tWLAAJycnChUqRLNmzbhy5Qpr1qzh119/falsu3btMDY2pn379pibm7/V6zRp0oS+ffvSq1cvSpUqxbFjxxg+fPg7x/06pqam7N69GwcHB+rXr4+Pjw8TJ07EyMgIIyMjwsLC6NixI15eXrRu3Zp69eoxevToDI9DCGGggk7C4rrwPBjsi0K33ZAne0y4UClvu6hNNhUdHY2trS1RUVEvJQ0JCQncvXuXAgUKvPUXochaHjx4gLu7O6dPn6Z06dL6DidLk587ITLBjd2wtiOkxEP+ctB+DVi+elLGx+JN39//JgOthcgAycnJBAcHM3jwYMqXLy8JkRDi43NhDWz6GpRU8KgNrZeB6cc1Yeh9SfeZEBngxSDps2fPMnfuXH2HI4QQb+f4r7DxC21CVKINtFuV7RIikJYiITJEtWrV3np7FSGE0DtFgX1j4Mg07XH5nlBnHKizZ5uJJEVCCCFEdpSaAtv6wrll2uOaI+HTvllq2463JUmREEIIkd0kxcH6bnB9O6jU0HAGlOmk76j0TpIiIYQQIjuJC4eVbeDhKTA2hxaLoGjD/35eNiBJkRBCCJFdRAZpt+14dgPMbaHdGnCroO+oDIYkRUIIIUR2EHIZVrTULspo4wyfrQcHWa3+nwxyePmjR4/47LPPyJ07N5aWlpQqVYqzZ8/qrj958oTOnTuTL18+LC0tqVu3Ljdv3nxjnUuXLkWlUr30SEhIyOzbEem0dOlS7Ozs9B2G3mTU/VepUoWVK1e+f0D/oX///q/c6FcIYYDuHoYl9bQJkYM3dNsjCdErGFxSFBERQaVKlTAxMWHHjh1cuXKFqVOn6r4sFEWhadOm3Llzh82bN3P+/Hnc3NyoVatWmk1DX8XGxobg4OA0D1kJV7sKc7du3ciXLx+mpqa4ubnRp08fwsLCPmgcbdq04caNG28sM2rUKEqVKvXS+Xv37qFSqf5zs1hDlp77/y9bt24lJCSEtm3bZlBUrzdw4ECWLFnC3bt3M/21hBDvIXAj/N4cEqPBtSJ02Q62zvqOyiAZXPfZpEmTcHFxSbMrt7u7u+7vN2/e5MSJE1y+fJlixYoB8Ouvv+Lg4MCqVavo3r37a+tWqVQ4OjpmWuwfozt37lChQgW8vLxYtWoVBQoUIDAwkAEDBrBjxw5OnDhBrlwfZol3CwuL99q49X0lJSVhamqqt9fPiPufNWsWXbp0Qf0B1hhxcHCgTp06zJ07l0mTJmX66wkh3sHJebBjEKBA0UbQfCGYSGPA6xhcS9GWLVvw8/OjVatWODg44Ovry4IFC3TXExMTAdK08BgZGWFqasqRI0feWHdMTAxubm7kz5+fhg0bcv78+deWTUxMJDo6Os0jK/rmm290G4RWrVoVV1dX6tWrx969e3n06BE//PCDrqxKpWLTpk1pnm9nZ5dm1/pBgwbh5eWFpaUlBQsWZPjw4SQnJ+uuX7hwgerVq2NtbY2NjQ1lypThzJkzQMZ2n6WmptKtWzcKFCiAhYUFhQsXZubMmWnKdO7cmaZNmzJhwgTy5cuHl5eXrsVp7dq1VK5cGQsLC8qWLcuNGzc4ffo0fn5+5MiRg7p16/L06VNdXRqNhjFjxpA/f37MzMwoVaoUO3fu1F1/Ue+GDRuoXr06lpaWlCxZkuPHj+vKvOr+X/w8mJubkydPHpo3b/7ae3727Bl79+6lcePGac7fvHmTKlWqYG5ujre3N3v27Enzb/mu9wzQuHFjVq1ala5/EyHEB6QosHcU7BgIKFC2O7T6TRKi/2BwSdGdO3eYM2cOnp6e7Nq1i6+++orevXuzbJl2cakiRYrg5ubGkCFDiIiIICkpiYkTJxISEkJwcPBr6y1SpAhLly5ly5YtrFq1CnNzcypVqvTasUgTJkzA1tZW93BxcXmr+1AUBU1cnF4e6V1ZOTw8nF27dtGzZ8+XWigcHR3p0KEDa9aseauVmq2trVm6dClXrlxh5syZLFiwgOnTp+uud+jQgfz583P69GnOnj3L4MGDMTExSXf96aXRaMifPz9r167lypUrjBgxgqFDh7J27do05fbt28fVq1fZs2cPW7du1Z0fOXIkw4YN49y5cxgbG9OuXTsGDhzIzJkzOXz4MLdv32bEiBG68jNnzmTq1KlMmTKFixcv4u/vT+PGjV/6//XDDz/Qv39/AgIC8PLyol27dqSkpLzyHrZt20bz5s1p0KAB58+fZ9++ffj5+b32no8cOYKlpSVFi/5/nIBGo6F58+YYGRlx4sQJ5s6dy6BBg175/Le9Z4By5crx4MED7t+//9q4hBAfWGqydg+zI39/9tYYDvWngNpIv3F9DBQDY2JiolSoUCHNuW+//VYpX7687vjMmTNKyZIlFUAxMjJS/P39lXr16in16tVL9+ukpqYqJUuWVL799ttXXk9ISFCioqJ0jwcPHiiAEhUV9VLZ+Ph45cqVK0p8fPz/64+NVa4ULqKXR2psbLregxMnTiiAsnHjxldenzZtmgIoT548URRFeWVZW1tbZcmSJa99jcmTJytlypTRHVtbWytLly59ZdklS5Yotra2b4x55MiRilqtVqysrNI8LC0tFUA5f/78a5/bs2dPpUWLFrrjTp06KXnz5lUSExN15+7evasAysKFC3XnVq1apQDKvn37dOcmTJigFC5cWHecL18+Zdy4cWler2zZskrPnj1fW29gYKACKFevXn3l/VeoUEHp0KHDG9+Pf5o+fbpSsGDBNOd27dqlGBkZKQ8ePNCd27FjR5p/y3e9Z0VRlKioKAVQDhw4kO44M8qrfu6EyPYSnivKsmaKMtJGUUblVJRzy/Udkd69+Jx61ff3vxlcS5GTkxPe3t5pzhUtWpSgoCDdcZkyZQgICCAyMpLg4GB27txJWFgYBQoUSPfrqNVqypYt+9qWIjMzM2xsbNI8shvl7xaitxln88cff/Dpp5/i6OhIjhw5GD58eJp/u379+tG9e3dq1arFxIkTuX379ivrCQoKIkeOHLrH+PHjddcKFy5MQEBAmsf27dtfqmPu3Ln4+flhb29Pjhw5WLBgQZpYAHx8fF55fyVKlND9PW/evLqy/zwXGhoKQHR0NI8fP6ZSpUpp6qhUqRJXr159bb1OTk4Aunr+LSAggJo1a77y2qvEx8e/NHHg6tWruLq6kj9/ft25ChVevSbJ29zzCy9aGOPi4tIdpxAik8Q8hd8awu19YGKp3dTV9zN9R/VRMbiB1pUqVeL69etpzt24cQM3N7eXytra2gLaMRNnzpxh7Nix6X4dRVEICAhI86GfkVQWFhQ+d/a/C2bSa6eHh4cHKpWKK1eu0LRp05euX7t2DXt7e904F5VK9VJX2j/HC504cYK2bdsyevRo/P39sbW1ZfXq1UydOlVXZtSoUbRv355t27axY8cORo4cyerVq2nWrFmaevPly5dmJtk/B3ubmpri4eGRpryxcdr/ymvXrqVv375MnTqVChUqYG1tzU8//cTJkyfTlLOysnrle/PPLj3V3/sA/fucRqNJ8xzVv/YLUhTlpXOvqvff9bzwtoOu8+TJQ0RExEsx/Nu/Y3pTbP91z+Hh4QDY29u/VaxCiAwWfke7KGP4HbDIBR3WQf7Xd7eLVzO4pKhv375UrFiR8ePH07p1a06dOsX8+fOZP3++rsy6deuwt7fH1dWVS5cu0adPH5o2bUqdOnV0ZTp27IizszMTJkwAYPTo0ZQvXx5PT0+io6OZNWsWAQEB/PLLL5lyHyqVCpWlZabUnVFy585N7dq1+fXXX+nbt2+aL+GQkBBWrFjBN998oztnb2+fZtzWzZs307QQHD16FDc3tzSDs1811sTLywsvLy/69u1Lu3btWLJkyUtJkbGx8UuJz9s4fPgwFStWpGfPnrpzr2uVel82Njbky5ePI0eOUKVKFd35Y8eOUa5cuXeut0SJEuzbt48uXbqkq7yvry8hISFERESQM2dOALy9vQkKCuLx48fky5cPIM3g7vd1+fJlTExMdDNBhRB68Pg8rGgFsU/BzhU+2wh53v3zMzszuO6zsmXLsnHjRlatWkXx4sUZO3YsM2bMoEOHDroywcHBfP755xQpUoTevXvz+eefvzQDJigoKM0XeGRkJF988QVFixalTp06PHr0iEOHDr3Xl1ZW8PPPP5OYmIi/vz+HDh3iwYMH7Ny5k9q1a+Pl5ZVmYG2NGjX4+eefOXfuHGfOnOGrr75K05Lg4eFBUFAQq1ev5vbt28yaNYuNGzfqrsfHx9OrVy8OHDjA/fv3OXr0KKdPn04zMDijeHh4cObMGXbt2sWNGzcYPnw4p0+fzvDXeWHAgAFMmjSJNWvWcP36dQYPHkxAQAB9+vR55zpHjhzJqlWrGDlyJFevXuXSpUtMnjz5teV9fX2xt7fn6NGjunO1atWicOHCdOzYkQsXLnD48OE0Sev7Onz4sG7GmhBCD27thaUNtQmRo492UUZJiN6ZwbUUATRs2JCGDV+/OV3v3r3/cyXdAwcOpDmePn16mllQQsvT05PTp08zatQoWrduTWhoKIqi0Lx5c5YvX47lP1q7pk6dSpcuXahSpQr58uVj5syZaVYab9KkCX379qVXr14kJibSoEEDhg8fzqhRowDt0glhYWF07NiRJ0+e6KaYjx49OsPv66uvviIgIIA2bdqgUqlo164dPXv2ZMeOHRn+WqD9PxkdHc33339PaGgo3t7ebNmyBU9Pz3eus1q1aqxbt46xY8cyceJEbGxs0rRE/ZuRkRFdu3ZlxYoVup8ftVrNxo0b6datG+XKlcPd3Z1Zs2ZRt27dd47rn1atWpUp/35CiHQ4txz+7ANKKhSoCm1+B/PsN/41I6mUVw06EC+Jjo7G1taWqKiolwZdJyQkcPfuXQoUKJAlVsgeOXIk06ZNY/fu3a8dlCsM05MnTyhWrBhnz5595Ti8F1QqFRs3bnzlWLL02rZtGwMGDODixYsvjen6ELLaz50Q6aYocGAiHJyoPS7RFhrPBmP9LT5ryN70/f1vBtlSJPRr9OjRuLu7c/LkST755JMPsjqyyBh58+Zl0aJFBAUFvTEpygixsbEsWbJELwmRENlWajL8+R0E/K49rtwfagyD10ygEG9HPs3EK6V3cK8wPE2aNPkgr9O6desP8jpCiL8lPoe1HeH2X6BSQ4Np4Cef1RlJkiIhsiHpNRfiIxMdDCtbQcgl7RpErZaCl7++o8pyJCkSQgghDFnoVe2U+6gHYGUP7deCc2l9R5UlSVIkhBBCGKq7h2F1B0iMgtwe0OEPyJX+3RvE25GkSAghhDBEl/7QbuyamgQu5bXbdljm+u/niXcmSZEQQghhSBQFjs2CPX8vnlu0MTSfDyaySGpmk6RICCGEMBSaVNgxCE4v0B6X7wl1fgS1kX7jyiYkKRJCCCEMQVIcrO8O17cBKvAfBxW++c+niYwjSZEQQgihb7HPYGUbeHQGjMy03WXFmuo7qmxHlioWBmPp0qXY2dnpOwy9yaj7r1KlCitXrvygr5kenTt3Tve2IqGhodjb2/Po0aPMDUoIQxB2GxbV1iZEFjmh42ZJiPREkiLBgwcP6NatG/ny5cPU1BQ3Nzf69OlDWFjYB42jTZs23Lhx441lRo0aRalSpV46f+/ePVQqFQEBAZkT3AeQnvv/L1u3biUkJIS2bdtmUFQZZ+bMmSxdujRdZR0cHPj8888ZOXJk5gYlhL7dPw4La0L4HbBzha67wU32nNQXSYqyuTt37uDn58eNGzdYtWoVt27dYu7cuezbt48KFSoQHh7+wWKxsLDAwcHhg73evyUlJenttSFj7n/WrFl06dLFIPers7W1fatWqS5durBixQoiIiIyLygh9OnSH7CsMcRHgHMZ6L4P7L30HVW2ZnifnFmEoijEJaXo5fE2Wzh88803mJqasnv3bqpWrYqrqyv16tVj7969PHr0iB9++EFXVqVSsWnTpjTPt7OzS/Pb/6BBg/Dy8sLS0pKCBQsyfPhwkpOTddcvXLhA9erVsba2xsbGhjJlynDmzBkgY7tyUlNT6datGwUKFMDCwoLChQszc+bMNGVedOdMmDCBfPny4eXlpWtxWrt2LZUrV8bCwoKyZcty48YNTp8+jZ+fHzly5KBu3bo8ffpUV5dGo2HMmDHkz58fMzMzSpUqxc6dO3XXX9S7YcMGqlevjqWlJSVLluT48eO6Mq+6/y1btuDn54e5uTl58uShefPmr73nZ8+esXfvXho3bpzmfGRkJF988QV58+bF3Nyc4sWLs3Xr1lfWcfv2bZo0aULevHnJkSMHZcuWZe/evWnK/Prrr3h6emJubk7evHlp2bKl7toff/yBj48PFhYW5M6dm1q1ahEbG5vm/f7nezZp0iQ8PDwwMzPD1dWVcePG6a77+Pjg6OjIxo0bX3vPQnyUFAUOTYH13bRrEBVpCJ22Qg79/VIotGSgdSaJT07Fe8Quvbz2lTH+WJr+9z9teHg4u3btYty4cVhYpF3/wtHRkQ4dOrBmzRp+/fVXVOncgdna2pqlS5eSL18+Ll26RI8ePbC2tmbgwIEAdOjQAV9fX+bMmYORkREBAQGYmJi8/U3+B41GQ/78+Vm7di158uTh2LFjfPHFFzg5OaXZyHTfvn3Y2NiwZ8+eNMnkyJEjmTFjBq6urnTt2pV27dphY2PDzJkzsbS0pHXr1owYMYI5c+YA2q6hqVOnMm/ePHx9fVm8eDGNGzcmMDAQT09PXb0//PADU6ZMwdPTkx9++IF27dpx69atV+40v23bNpo3b84PP/zA8uXLSUpKYtu2ba+95yNHjmBpaUnRokXTvA/16tXj+fPn/P777xQqVIgrV65gZPTq6b0xMTHUr1+fH3/8EXNzc3777TcaNWrE9evXcXV15cyZM/Tu3Zvly5dTsWJFwsPDOXz4MADBwcG0a9eOyZMn06xZM54/f87hw4dfm6QPGTKEBQsWMH36dD799FOCg4O5du1amjLlypXj8OHDdO3a9bX3LcRHJTUZtvaF88u1xxV6Qe0xMuXeQEhSlI3dvHkTRVHSfIn+U9GiRYmIiODp06fp7tYZNmyY7u/u7u58//33rFmzRpcUBQUFMWDAAIoUKQKQJmFIr0uXLpEjR4405/79xWtiYsLo0aN1xwUKFODYsWOsXbs2TVJkZWXFwoULMTU1BbQtOgD9+/fH31+72WKfPn1o164d+/bto1KlSgB069YtTQvZlClTGDRokG4sz6RJk9i/fz8zZszgl19+0ZXr378/DRo0AGD06NEUK1aMW7du6d6Pfxo3bhxt27ZNcx8lS5Z87fty79498ubNm6brbO/evZw6dYqrV6/i5aVtli9YsOBr6yhZsmSa1/jxxx/ZuHEjW7ZsoVevXgQFBWFlZUXDhg2xtrbGzc0NX19fQJsUpaSk0Lx5c9zc3ABta8+rPH/+nJkzZ/Lzzz/TqVMnAAoVKsSnn36appyzszPnz59/bbxCfFQSomBtJ7izX7vLfb3JUK6HvqMS/yBJUSaxMDHiyhj97GBsYZIxv3G8SDReJAzp8ccffzBjxgxu3bpFTEwMKSkp2NjY6K7369eP7t27s3z5cmrVqkWrVq0oVKjQS/UEBQXh7e2tOx46dChDhw4FoHDhwmzZsiVN+UePHlGtWrU05+bOncvChQu5f/8+8fHxJCUlvTRI28fH55X3V6JECd3f8+bNqyv7z3OhoaEAREdH8/jxY13C9EKlSpW4cOHCa+t1cnICtDOtXpUUBQQE0KNH+j8w4+PjMTc3f6mO/Pnz6xKi/xIbG8vo0aPZunUrjx8/JiUlhfj4eIKCggCoXbs2bm5uFCxYkLp161K3bl2aNWum6w6sWbMmPj4++Pv7U6dOHVq2bEnOnDlfep2rV6+SmJhIzZo13xiPhYUFcXFx6XwHhDBgkQ9gZWsIvQImVtByMRSuq++oxL/ImKJMolKpsDQ11ssjvV1dHh4eqFQqrly58srr165dw97eXjfORaVSvdQi88/xQidOnKBt27bUq1ePrVu3cv78eX744Yc0A5hHjRpFYGAgDRo04K+//sLb2/uVY0by5ctHQECA7vHVV1/prpmamuLh4ZHm8aJl4oW1a9fSt29funbtyu7duwkICKBLly4vDaa2srJ65b3/s0vvxfv573MajSbNc/79viuK8tK5V9X773pe+HeX5n/JkyfPS4OS37aOAQMGsH79esaNG8fhw4cJCAjAx8dH975ZW1tz7tw5Vq1ahZOTEyNGjKBkyZJERkZiZGTEnj172LFjB97e3syePZvChQtz9+7dd7638PBw7O3t3+oehDA4j89rZ5iFXoEcjtBluyREBkqSomwsd+7c1K5dm19//ZX4+Pg010JCQlixYgWdO3fWnbO3tyc4OFh3fPPmzTS/xR89ehQ3Nzd++OEH/Pz88PT05P79+y+9rpeXF3379mX37t00b96cJUuWvFTG2Ng4TdKTK9fbbYJ4+PBhKlasSM+ePfH19cXDw4Pbt2+/VR3pZWNjQ758+Thy5Eia88eOHXtt12R6lChRgn379qW7vK+vLyEhIWkSoxIlSvDw4cN0T/U/fPgwnTt3plmzZrqBzi+6FF8wNjamVq1aTJ48mYsXL3Lv3j3++usvQJvoVapUidGjR3P+/HlMTU1fmfR6enpiYWHxn/d3+fJlXfecEB+l6ztgSX2IeQIOxaDHPshXSt9RideQ7rNs7ueff6ZixYr4+/vz448/UqBAAQIDAxkwYABeXl6MGDFCV7ZGjRr8/PPPlC9fHo1Gw6BBg9K0fHh4eBAUFMTq1aspW7Ys27ZtS/OFGB8fz4ABA2jZsiUFChTg4cOHnD59mhYtWmT4fXl4eLBs2TJ27dpFgQIFWL58OadPn6ZAgQIZ/lqgbWEZOXIkhQoVolSpUixZsoSAgABWrFjxznWOHDmSmjVrUqhQIdq2bUtKSgo7duzQjc/6N19fX+zt7Tl69CgNGzYEoGrVqlSpUoUWLVowbdo0PDw8uHbtGiqVirp1X/5N1cPDgw0bNtCoUSNUKhXDhw9P05K1detW7ty5Q5UqVciZMyfbt29Ho9FQuHBhTp48yb59+6hTpw4ODg6cPHmSp0+fvjIxNDc3Z9CgQQwcOBBTU1MqVarE06dPCQwMpFu3bgDExcVx9uxZxo8f/87voRB6dXIe7BwMigYK1YRWS8Hc5j+fJvRHWoqyOU9PT06fPk3BggVp3bo1bm5u1KtXDy8vL44ePZpmQPPUqVNxcXGhSpUqtG/fnv79+2Npaam73qRJE/r27UuvXr0oVaoUx44dY/jw4brrRkZGhIWF0bFjR7y8vGjdujX16tVLM5A4o3z11Vc0b96cNm3a8MknnxAWFkbPnj0z/HVe6N27N99//z3ff/89Pj4+7Ny5ky1btrzTQPIXqlWrxrp169iyZQulSpWiRo0anDx58rXljYyM6Nq160uJ2Pr16ylbtizt2rXD29ubgQMHkpqa+so6pk+fTs6cOalYsSKNGjXC39+f0qVL667b2dmxYcMGatSoQdGiRZk7dy6rVq2iWLFi2NjYcOjQIerXr4+XlxfDhg1j6tSp1KtX75WvNXz4cL7//ntGjBhB0aJFadOmjW6cFsDmzZtxdXWlcuXKb/O2CaF/mlTYOQR2DNQmRKU7Qfs1khB9BFTK2yxqk41FR0dja2tLVFRUmoHDAAkJCdy9e5cCBQq8NND1YzRy5EimTZvG7t27qVBBVlb9mDx58oRixYpx9uzZl8ZZfWzKlSvHd999R/v27V95Pav93IksIikW1vf4e1NXoNYoqPQdpHOsp8h4b/r+/jfpPhMvGT16NO7u7pw8eZJPPvnEIFdHFq+WN29eFi1aRFBQ0EedFIWGhtKyZUvatWun71CESL/nT2BVG+3AaiMzaDYHimf88ACReaSlKJ2yU0uREB8D+bkTBuVJoHaX+6gHYJEL2q0C1/L6jkogLUVCCCHEh3NzL6zrDEnPIVch6LAOcr+8/powfJIUCSGEEO/q5HzYOUg7oNq9MrReBpZvt4SIMBySFAkhhBBvKzUFdg2FU/O0x6U+g4bTwTj9OwAIwyNJkRBCCPE2Ep/DH13h5m7tcc2R8GlfmWGWBUhSJIQQQqRX5APtgOrQQDC2gObzwLuJvqMSGUSSIiGEECI9Hp6FVW0hNhRy5NXOMHMuo++oRAaSBWhElqBSqdi0aZNeYzhw4AAqlYrIyMh0P8fd3Z0ZM2ZkWkxCiAwSuAmW1tcmRA7FoPs+SYiyIEmKsrnOnTujUqnS7EL/Qs+ePVGpVGk2hTVUwcHBr91OArLOfQohPjBFgcNTYV0nSEkAT3/otgvsXPQdmcgEkhQJXFxcWL16NfHx8bpzCQkJrFq1CldX1/eqW1EUUlJS3jfE/+To6IiZmdkby2TmfQohsqCUJNjUE/aN0R5/8rW2y8zMWr9xiUwjSZGgdOnSuLq6smHDBt25DRs24OLigq+vb5qyiYmJ9O7dGwcHB8zNzfn00085ffq07vqLLqRdu3bh5+eHmZkZhw8fpnPnzjRt2jRNXd999x3VqlUD4N69e6hUqpceL65Xq1btldfv3bsHpK/7LCPvE2D79u14eXlhYWFB9erVdbH807Fjx6hSpQoWFha4uLjQu3dvYmNj3xinEMIAxIXD8qZwYSWo1FB/CtSbCGojfUcmMpEkRZlFUbQbA+rj8Q47t3Tp0oUlS5bojhcvXkzXrl1fKjdw4EDWr1/Pb7/9xrlz5/Dw8MDf35/w8PCXyk2YMIGrV69SokSJ/3x9FxcXgoODdY/z58+TO3duqlSpAmiTl39eb968OYULFyZv3rx6uc8HDx7QvHlz6tevT0BAAN27d2fw4MFp6rh06RL+/v40b96cixcvsmbNGo4cOUKvXr3eKmYhxAf27BYsrAn3j4KpNbRfB+V66Dsq8QHI7LPMkhwH4/Pp57WHPgZTq7d6yueff86QIUN0LTZHjx5l9erVHDhwQFcmNjaWOXPmsHTpUt34nQULFrBnzx4WLVrEgAEDdGXHjBlD7dq10/36RkZGODo6AtouraZNm1KhQgVGjRoFQK5c/18hdvr06fz111+cPHkSCwsLvdznnDlzKFiwINOnT0elUlG4cGEuXbrEpEmTdPX89NNPtG/fnu+++w4AT09PZs2aRdWqVZkzZ47s1yWEIbp7GNZ8BgmRYOsK7ddAXm99R5WlJacm8yDmAfei7hGVGEUzz2Z6i0WSIgFAnjx5aNCgAb/99huKotCgQQPy5MmTpszt27dJTk6mUqVKunMmJiaUK1eOq1evpinr5+f3zrF069aN58+fs2fPHtTqtI2ZO3bsYPDgwfz55594eXm9dd0ZdZ9Xr16lfPnyqP6xWFuFChXS1HP27Flu3brFihUrdOcURUGj0XD37l2KFi361vELITLR2aWw7XvQpICzn3b8UA4HfUeVJSiKQlhCGPei7nEv+t7//4y+x8PnD0lVUgGwMLagqUfTNJ+tH5JBJkWPHj1i0KBB7Nixg/j4eLy8vFi0aBFlyminPz558oRBgwaxe/duIiMjqVKlCrNnz8bT0/ON9a5fv57hw4dz+/ZtChUqxLhx42jWLJMyUhNLbYuNPphYvtPTunbtquva+eWXX166rvzdLffv/6yKorx0zsoqbUuVWq3WPf+F5OTkl17jxx9/ZOfOnZw6dQpr67SDGa9cuULbtm2ZOHEiderUSeddvSwj7vPf9/IqGo2GL7/8kt69e790TQZ2C2FAUlNgz3A48av2uHgLaPILmLxdS7TQSkhJYP7F+ZwPPc+D5w94EvfkP59jYWyBu4077jbuxKfEY/mO32Pvy+CSooiICCpVqkT16tXZsWMHDg4O3L59Gzs7O0D7ZdS0aVNMTEzYvHkzNjY2TJs2jVq1anHlypWXvoxfOH78OG3atGHs2LE0a9aMjRs30rp1a44cOcInn3yS8TeiUr11F5a+1a1bl6SkJAD8/f1fuu7h4YGpqSlHjhyhffv2gDaxOXPmjK6L6HXs7e25fPlymnMBAQGYmJjojtevX8+YMWPYsWMHhQql3WE6LCyMRo0a0bx5c/r27fsut6eTEffp7e390sDuEydOpDkuXbo0gYGBeHh4vFe8QohMlBCl3bLj1l7tcfUfoMoA2bLjPYw/OZ6Ntzamq2y/Mv2oV6AeeS3z6q116J8MLimaNGkSLi4uaQbDuru76/5+8+ZNTpw4weXLlylWrBgAv/76Kw4ODqxatYru3bu/st4ZM2ZQu3ZthgwZAsCQIUM4ePAgM2bMYNWqVZl3Qx8RIyMjXfeQkdHLMyysrKz4+uuvGTBgALly5cLV1ZXJkycTFxdHt27d3lh3jRo1+Omnn1i2bBkVKlTg999/5/Lly7pZX5cvX6Zjx44MGjSIYsWKERISAoCpqSm5cuWiefPmWFhYMGrUKN010CZbr4o1s+/zq6++YurUqfTr148vv/ySs2fPsnTp0jT1DBo0iPLly/PNN9/Qo0cPrKysuHr1Knv27GH27NlvFbMQIhOE34GVbeHZde2WHc3mQrGm+o7qo1fNpVq6k6JabrVwtHLM5IjSz+Bmn23ZsgU/Pz9atWqFg4MDvr6+LFiwQHc9MTERIM0gVSMjI91v9q9z/Pjxl7pc/P39OXbs2CvLJyYmEh0dneaRHdjY2GBjY/Pa6xMnTqRFixZ8/vnnlC5dmlu3brFr1y5y5sz5xnr9/f0ZPnw4AwcOpGzZsjx//pyOHTvqrp85c4a4uDh+/PFHnJycdI/mzZsDcOjQIQIDA3F3d09z/cGDB3q5T1dXV9avX8+ff/5JyZIlmTt3LuPHj09TR4kSJTh48CA3b96kcuXK+Pr6Mnz4cJycnN4pZiFEBrp3BBbU1CZE1k7QdYckRBlAo2hwsnKimks1AIzVxhTJVYQ6bnXo4dODHyv9yPJ6yznc5jCXOl3CxdqwFsFUKekZHPEBvUh2+vXrR6tWrTh16hTfffcd8+bNo2PHjiQnJ+Pp6Um5cuWYN28eVlZWTJs2jSFDhlCnTh127dr1ynpNTU1ZunSprjsEYOXKlXTp0kWXaP3TqFGjGD169Evno6KiXvoyTUhI4O7duxQoUEBmFAnxgcjPnXhn55bB1r7aAdX5fKHtKrCRX1beVVRiFMceH+PIoyMcfXSUsIQwAEzUJvxc82cq5quo1/iio6OxtbV95ff3vxlc95lGo8HPz0/3W7evry+BgYHMmTOHjh07YmJiwvr16+nWrRu5cuXCyMiIWrVqvXGLhxfSM0D4hSFDhtCvXz/dcXR0NC4uhpXRCiGEeAuaVNgzAo7/rD0u1gya/Aqm+hnU+7G7EHyOz3Z1fOX4q3xW+fjx0x8p61hWD5G9O4PrPnNycsLbO+2aEEWLFiUoKEh3XKZMGQICAoiMjCQ4OJidO3cSFhZGgQIFXluvo6NjmrEoAKGhoa9d/M/MzEzXxfJfXS1CCCEMXEK0dof7FwlRtSHQcokkRO8o8c4dYjt/S+XLr+5s+qH8Dx9dQgQGmBRVqlSJ69evpzl348YN3NzcXipra2uLvb09N2/e5MyZMzRp0uS19VaoUIE9e/akObd7924qVtRvs54QQohMFn4XFtWBm7vB2FybDFUbLDPM3kFcchynfx3LjaaNyXk/nNZHNBilahMjY5Uxtd1qs6jOIqrkr6LnSN+NwXWf9e3bl4oVKzJ+/Hhat27NqVOnmD9/PvPnz9eVWbduHfb29ri6unLp0iX69OlD06ZN0wyk7tixI87OzkyYMAGAPn36UKVKFSZNmkSTJk3YvHkze/fufePgbCGEEB+5+8e0K1THhUEOR2i3EpzL6Duqj45G0TBm+/d4ztlN6VsaAC64q5jXyBTvvN4srLNQb2sLZSSDS4rKli3Lxo0bGTJkCGPGjKFAgQLMmDGDDh066MoEBwfTr18/njx5gpOTEx07dmT48OFp6gkKCkqzGnLFihVZvXo1w4YNY/jw4RQqVIg1a9Zk6BpFBjZmXYgsTX7exH86/zv8+R1oksGpJLRbDTZ62n7pI3dj6yr8R+7ELg6SjeFyK1/sO3bmT+eK5DDNoe/wMozBzT4zVG8avZ6amsqNGzdwcHAgd+7ceopQiOwlKiqKx48f4+HhkWYRUCHQpMLekXDs7/XAvJtA07kyfugdhIY94NGkCZhv2Q9AkD0k/PA1ZT7pTr81F+hSyZ16PoY9c++jnn32MTIyMsLOzo7Q0FAALC0tDWJlTiGyKo1Gw9OnT7G0tMTYWD7GxD8kRMH67trxQwBVB0HVwaA2uCG0Bi/ozAFu9OmJc5i27WRbWRUrq6lR3wsk8cRhniek8CgynlreeTExyhrvr3yaZJAXO7y/SIyEEJlLrVbj6uoqv4CI/wu7rZ1h9uyGdkB1k1/Ap6W+ozJ4scmxRCdGE5McQ2xyLM/jozBasxXb37bjnKoQngMOdSpJlE9RbC+5cveRA5CCr6sdM9v4ZpmECCQpyjAqlQonJyccHBxeudGpECJjmZqaphk3KLK52/thXWdIiATrfNB2BTiX1ndUBm9Z4DKmnJmCgrY1KHe0wjd/aigepD0+6aViWSNLptT5le/WXCAoPA61CnrV8KR3DQ+Ms1BCBJIUZTgjI6O33otLCCHEO1IUODkPdg0FJRWc/bQJkbXh7KdlyIzVxrqEqMIVDT12aciRAAkmsKS2mr981LhpOtNy7jE0igpnO3NmtPWlrHsuPUeeOSQpEkII8XFKSYRt38P55drjku2g4QwwkW1f0qt90fbUylOJB6NHYLn3FAC38qmY1UjNY+ucJAS1ITBeuzCysU0An9XJT1n3mvoMOVNJUiSEEOLjE/NUu/7QgxOgUkPtsVDhG1mQ8S3FnTtH9ICBWD56hKJWcb5uQX4u9YSISC8S7jQDjQWoEzB33ISxTQBhSW30HXKmkqRICCHExyX4IqxqB9EPwcwWWi4Gz1r6jsrgKYpCTHIMz+Kf8Sw6hJSFK7Bd9xcqjUKoLcxurOaaUwgJj5qQEqUdj+Vmr+Gn1r645fqU21G3KZGnhJ7vInNJUiSEEOLjEbgJNn0NyXGQ20O7IGMeT31HZTAexzzmZPBJnsY/1SY/8c94Gvf/vyekJuD8TOHbLakUfKJ9zsHiKhbXUROT6obqQRdS4i1eOZg6r9Wr9wrNSiQpEkIIYfg0Gjg4EQ5O0h4XqgktF4FFTv3GZWC67urKo5hHr7ymUhTqnVHosF+DaSo8N4dNzR2x8a+H9203Dl3KgUYBZzsLZrQtlWUHU7+JJEVCCCEMW2IMbPoKrv6pPa7QC2qNBiP5Cvu3dkXasf3udp7FPeNZwjM0inafstzRCl9v01Dinnam2bmCKuY2UBNumorrSV+uPk4CoHHJfIxtWhxbi+y5Srxs85FOb7NMuBBCiAwScR9Wt4cnl8HIVDu7zLfDfz5NQIomhYtPLzJ/6ud036XBKhESjWFZTTW7S6lIee5LQkhT0JiTw8yYsU2L0cw3v77DznCyzYcQQoiP372jsPZz7Q73Vg7Q5ndwzbhNvLO6mQfGYTt7DX2uats+bjrBz42MaFtzGMlXCvPn9ccAlHa1Y0YbX1xzy95wkhQJIYQwPKcXwY6BoEnR7nDfdiXYZr1WjIyWrEnm+OPjpJ48R6VJ67CNVkhVwR+fqtlUUU0tp2+Zty03T6IfY6RW0aemJz2rFcpyK1O/K0mKhBBCGI6UJNgxAM4u1R4Xa67dw0x2uP9Pp4JP8dPhH6mw+Rb1zmpbhx7ngtmNjLjlaExqWAPW33BGUZJwz23J9Dal8HWVger/JEmREEIIw/D8iba77MFJQAW1RkKl72RBxnT4LfA31v/5E9/+mUr+MO25c5XzcrCxGymxZsRdLYcm0QmAduVcGdagKFZmkgL8m7wjQggh9O/RWVj9GTx/LAsyviUlJYWUJasZtz0VYw2o8+TGecJECleqRMKxe0w6eQ1NioZcVqZMalGC2t5Zf72hdyVJkRBCCP0KWAl/fgepiZCnMLRbBbkL6Tsqg5eQksDdy0dJGDWZiteCADheRMUC/0i+z53ApsWnOHLrGQDVC9szqWUJHKxlX7g3kaRICCGEfqQmw+7hcHKO9rhwfWg2D8xl2ZPX0SgaRh0bxaabG6hzVuGz/RrMUyDWDBbVUXOkmIrk5z6MWJNIXGIi5iZqhjXwpsMnrqikG/I/SVIkhBDiw4sNg3Wd4N5h7XHVwVB1EKhlFtSbPI17yqGzG/hh+/8XYrzormJOfTXPcpiTENyYlKgyAKjNH/Jzh4rU8nTTZ8gfFUmKhBBCfFghl7QLMkYGgWkOaDYXijbSd1QGTaNo+PPWFo4vnsDU7alY/r0Q4+811OwurSI53p2Eu21QknMBGkxzH8DUfh8pJm6Aj56j/3hIUiSEEOLDubweNn0DKfGQs4B2/JBDUX1HZdBuRNyg66rmfLFDQ6eb2tahO/lN2NjWhRRnR5we+HLjvjsKKuxypNKuShJuDmXxytmB0nlL6zn6j4skRUIIITKfJhX+GgtHpmuPC9XQzjCTDV1fSaNomHthLvsf7Mf22BWm7tRgEw8pathbx57vpu7H81kcfdcEcP1xNAAtSudnVGNvrM2z575lGUGSIiGEEJkrPhLWd4dbe7THFXtDrVGgNtJnVAbt5/M/s/LUfLrs0VAlUNs6dM9Bu03HxC6/sujYfSbvuk5SigY7SxPGN/Ohvo+TnqP++ElSJIQQIvM8vQ6r2kH4bTA2h8Y/Q4lW+o7KoG24uYGTm+cxZbuG3M8BtZrcPbpT5Jtv8IlNof+6C5y4cx+AaoXtmdyiBA42MtU+I0hSJIQQInNc3Qobv4Kk52Drot3QNV8pfUdlkJI1ycwJmMPBG7uovPkuw85rW4dM3dzIN2ki5iVLsuHcI0ZtCeR5YgoWJkYMa1iU9uVkqn1GkqRICCFExtKkwoEJcOgn7bFbJWj1G+Sw129cBuxu1F0ObZ9Pr62pOEZqzyU1r03h4ZOI1BjR9/dz7AwMAbS72k9rXQr3PFb6CziLkqRICCFExomPgPU9/j9+6JOvoc5YMJLBv68TEfWEoAk/MnpLKmogKY8NHj/NwKpCBf669oSBf1ziWUwixmoVfWt78WWVgrKrfSaRpEgIIUTGeBIIqztAxF3t+KFGs6BkG31HZbACnwWya8cvlJh3AOdn2u6y/T4q+K4VHqXLMmTDRVadegCAp0MOprcpRXFnW32GnOVJUiSEEOL9Xd4Am7+B5DiwdYW2v4NTSX1HZXCiEqPYfnc7f17bQLEtgTQ5rmCkwPMcxoR804QSdepinORFvZmHCQqPA6DbpwUY4F8YcxOZrZfZJCkSQgjx7lJTYN9oODZLe1ywGrRcApa59BqWIUnRpLD1zlbW3VjHxacXKRCi8M3WVFyfaq+HVS5GkbE/USqPGzP23mDuwdNoFMhna86U1iWpWCiPfm8gG5GkSAghxLuJDYM/usDdg9rjSt9BzRGy/tC/jDw2ki23t2CUqtDqqIbmx7StQ1GWsKCumlOFr2O09RscYwZyLSQWgOalnRnVuBg2shDjByVJkRBCiLf3OADWfA5RQWBiBU1+huLN9R2VwVh4aSHLryzHxtSGe9H3cHui0HNbKgWeaK8fL6Jiob+aaAs1SWGVSX7qT6QSS86/F2KsJwsx6oUkRUIIId7OhdXwZx9ISYBcBaHNCsjrre+oDMqW21sITwgnKjaM5icUWh7RYKyBaAtY5K/meFE1mqTcpDxoR1JsfgBqFHFgYgsfHKxlIUZ9kaRICCFE+qQmw64f4NQ87bFnHWi+ACzs9BqWIYpNisXlqcKIfTmxvfsMgKdlC3KiQ0mszBLIdzc3d+6VIDlVTQ4zY0Y09KaVX35ZiFHPJCkSQgjx32JCYW0nCDqmPa46CKoOBrWsl/NvSkoKFf4KofVhDSapz1Db2uI4bBhFGjbAIyqBQesvcv2mNlGqUDA3k1uWwCWXpZ6jFiBJkRBCiP/y8Ix2/NDzx2BqDc3nQZEG+o7KICXevs3t/t/R4aoGgKtFrWg0bwvG9vasP/eI0X9v02FuomZw3SJ0rOCOWi2tQ4ZCkiIhhBCvpihwZjHsHAypSZDHC9quhDye+o7M4CipqYQv/Y3QmTNQJSUTawZLaqvx7tCZSAtbhi4/y54r2lHWpVzsmNa6JAXtc+g5avFvkhQJIYR4WXI8bO0HF1Zqj4s0hKZzwNxGv3EZoMRbt3g0dCiJFy8BcK6ginn11URYqygWWwb/GYcIj03CxEjFd7Vkmw5DZpD/Ko8ePeKzzz4jd+7cWFpaUqpUKc6ePau7HhMTQ69evcifPz8WFhYULVqUOXPmvLHOpUuXolKpXnokJCRk9u0IIcTHJfwuLKqtTYhUaqg9RrvDvSREACSnJnMl7Ap/XF3D2h/acKNJYxIvXiLWDObUVzOxtZpwS0uMn3Zn8tZwwmOTKOpkw5Zen/JNdQ9JiAyYwbUURUREUKlSJapXr86OHTtwcHDg9u3b2NnZ6cr07duX/fv38/vvv+Pu7s7u3bvp2bMn+fLlo0mTJq+t28bGhuvXr6c5Z24uUx+FEELnxm7Y0AMSIsEyD7RcDAWr6jsqvbvw9AIbb27kStgVbkbexOlJMl9vTcVDu3E9ZwupmF9P2zrkoFQjOKgGEQmmqFXQs5oHvWt6YmosyZChM7ikaNKkSbi4uLBkyRLdOXd39zRljh8/TqdOnahWrRoAX3zxBfPmzePMmTNvTIpUKhWOjo6ZEbYQQnzcNBo4OEn7QAFnP2i9DGyd9R2ZQRhzfAw3Im5glKrQ9LhCi6PadYdizOH3OqZYNqzHWt9BzNzzkFWnggAoaG/F1FYl8XXNqefoRXoZXNq6ZcsW/Pz8aNWqFQ4ODvj6+rJgwYI0ZT799FO2bNnCo0ePUBSF/fv3c+PGDfz9/d9Yd0xMDG5ubuTPn5+GDRty/vz515ZNTEwkOjo6zUMIIbKkuHBY2RoOTgQUKNsdumyXhOgfvvf7niaakoz/TUObw9qE6LSnin49jPirmIZNl69Sb+YRXULUtVIBtn1bWRKij4xKURRF30H804vurH79+tGqVStOnTrFd999x7x58+jYsSMASUlJ9OjRg2XLlmFsbIxarWbhwoV8/vnnr633xIkT3Lp1Cx8fH6Kjo5k5cybbt2/nwoULeHq+PJNi1KhRjB49+qXzUVFR2NhIv7oQIosIvgBrPoPIIDA2h4YzoFQ7fUdlME6HnKbH9i40P6ah2TEFYw08N4fFddQc9VZhrLLCLLINIcGFAXC2s2BKq5JUKJRbz5GLF6Kjo7G1tU3X97fBJUWmpqb4+flx7Ngx3bnevXtz+vRpjh8/DsCUKVNYsGABU6ZMwc3NjUOHDjFkyBA2btxIrVq10vU6Go2G0qVLU6VKFWbNmvXS9cTERBITE3XH0dHRuLi4SFIkhMg6zq+Abf2023XkdIfWy8GphL6j0psUTQqnQ06TlJqEUw4n8ljkYfSSTtRddQv3UG2Zk4VVPPiiHh6F/FAlePLL7mgehMcD0OETV4bUL0oOM4MbmZKtvU1SZHD/ck5OTnh7p91Dp2jRoqxfvx6A+Ph4hg4dysaNG2nQQLt4WIkSJQgICGDKlCnpTorUajVly5bl5s2br7xuZmaGmZnZe9yJEEIYqJRE2DEIzv49dtPTX7sgo0X27urZc38PAw8NBMA4RaHlUQ09jmt3tH+xZ1n33otpk8ePn3ZdZ8mxuygK5LM1Z1LLElT2tNfzHYj3ZXBJUaVKlV6aIXbjxg3c3NwASE5OJjk5GfW/lpY3MjJCo9Gk+3UURSEgIAAfH5/3D1oIIT4WkQ9gbUd4fA5QQfWhULm/bNcB5LHIA0Chx9od7V20O3Fw0tsY+vVgdoUvufQglvqzDnP3WSwAbcu68EODolibm+grbJGBDC4p6tu3LxUrVmT8+PG0bt2aU6dOMX/+fObPnw9op9VXrVqVAQMGYGFhgZubGwcPHmTZsmVMmzZNV0/Hjh1xdnZmwoQJAIwePZry5cvj6elJdHQ0s2bNIiAggF9++UUv9ymEEB/cnQPwR1eICwNzO2ixCDzT17qeHXhauNLhr1QanVJQKxBpCac/96V7r4WoMWPKzussPKJtHXK0MWdiCx+qFXbQd9giAxlcUlS2bFk2btzIkCFDGDNmDAUKFGDGjBl06NBBV2b16tUMGTKEDh06EB4ejpubG+PGjeOrr77SlQkKCkrTmhQZGckXX3xBSEgItra2+Pr6cujQIcqVK/dB708IIT44jQaOToe/fgRFA04ltdPtc7rrOzKDkJiayO5NM7CZ9jtNwrTDbA8XU7GklppBtVtzPTiR79ed4s5TbetQyzL5Gd7QG1sLaR3KagxuoLWhepuBWkIIYTDiI2DjV3Bjp/bY9zOoPxVMZOHahJQENl5YSfTMX6hyKg6ACBs18+vAWU81zlbuVLP6ifmHbqNRwMHajAnNfahZNK+eIxdv46MeaC2EECKDPA7Qjh+KvA9GZlD/JyjTSd9R6V1yajKrr6/mxIY5tN0SQZ6/l6E7WNqMyhPmUyb8JKdObufRg8+YG3MbgGa+zoxs5I2dpakeIxeZ7b2TouTkZEJCQoiLi8Pe3p5cuXJlRFxCCCHelaLA2aXaGWapiWDnpu0uy1dK35EZhKXHZpM6cyG9Lms7SkLsYF49NYHuqaw9Poz6OWeTeL8gqQrkyWHKuGY++BeT3RCyg3dKimJiYlixYgWrVq3i1KlTadbzyZ8/P3Xq1OGLL76gbNmyGRaoEEKIdEiK0649dGGV9rhwfWj6a7afbg/aWcfPd+2m4ug1qCMUNCrYVlbF2spqEk1VpMY7ExP8NT8H3AKgUcl8jG5cjFxW0jqUXbx1UjR9+nTGjRuHu7s7jRs3ZvDgwTg7O2NhYUF4eDiXL1/m8OHD1K5dm/LlyzN79uxXrhgthBAigz27pe0uCw3U7m5fcwRU7JOtp9vHJsdyMvgkZwL3UGjRX3hfjkYNPHO0YFqdJO7kN6Ja/lqkPvNnx/UU4jQKua1MGdu0OPV9nPQdvvjA3nqgdatWrRgxYsR/ru+TmJjIokWLMDU1pXv37u8VpCGQgdZCCIMWuAk294Kk52DloN3dvkBlfUelFwGhAewL2sexx8e4EX6dapcUOu7TkCMBUtSwqYKKDRXVpBqr6eE1ji3HbbgZGgNoW4dGNfImdw5ZvDer+Ki3+TBUkhQJIQxSajLsGQkn/l5zza2SNiGyzp5jYFZcXcHEUxMBsI9U6LFTQ6m72q+5246gGdKTxAKOPIuN4totD/48F4tGgTw5zPixaXHqFs+e71tWJrPPhBAiO4h+DOs6w4OT2uOKvaHmSDDKvh/tOUxyoFIU6pxV6HBAg3kyJBnD2spqDn9qx8Yq7bkXqmbOnxd06w41LZWPkY2KkVPGDmV7b/WTExERgaIo5MqVi6dPn3Lo0CEKFy5M8eLFMys+IYQQr3LnAPzRDeKegZmtdjB10Yb6jkpvDj44SP+D/cn9JJ5Fh5zIcf0hAFdcYF49I0r51Wd5ie+Ys+8Ji45qV6W2tzZjfDMfanvLukNCK91J0cKFC5kwYQIajYaBAweyYsUKSpQowciRI+nduzdffPFFZsYphBACtKtTH5kK+8drV6fO6wNtlkGugvqOTK8uBZ+n4YE4mh/TYKx5SLwprKim5m51T8ZV+AElvhCdFlzgXph2kcbmpZ0Z0VDWHRJppTspmj17NoGBgcTFxeHq6srdu3ext7cnOjqaKlWqSFIkhBCZLS5cuzr1zV3aY9/PtQsymljoNy49UBSFdTfWcSbkDBbXgqi+8ip5n2g3BT/roWKhvxoTJyf+qLeK6bvv8Nvx47o9y8Y3L06NItI6JF6W7qTIyMgIc3NzzM3N8fDwwN7eHtBu0KpSqTItQCGEEMDDM9rxQ1EPwNgcGkzVbtmRTa25voaph36k3UEN/mcV1Gg3cF1SW83xoipQqeictxeNZp0gKFzbOtTaLz8/NJA9y8TrpTspMjY2JiEhAXNzcw4ePKg7//z580wJTAghBNrVqU/MgT0jQJMMOQtoV6d2KqHvyPTiedJzttzewq6V45m2S6PboiPBvwIW33ZliL0zJuRg7v5gZm97AMThZGvOhOayo734b+lOiv766y/MzLTrNtja2urOx8fHs2jRooyPTAghsrv4SNj8DVzbqj32bgKNZ4O57RufltWcDD7JyqsreRDzgNCHN+m8J5UhV7XT7NXOTuQf+yNWFSsCcPjmUwavv8SjyHgA2pVzZWj9IlibS+uQ+G/pTopy5MjxyvMODg44OEj2LYQQGepxAKzrBBH3QG0C/uOhXA/IZsMVIhMi6b67OygKVS8pDN6nwToBNCqIb1mL0kMno7awICoumR+3XWHdWe2sM2c7Cya28KGyp72e70B8TLLvYhZCCGGIFAVOL4RdQyE1CexcodVv4Fxa35Hpha2ZLV/kaYr97PWUvKdtHdJ4ulFwwlQsihcDYHdgCMM2XSb0eSIqFXSq4M4A/8JYmclXnHg7Gf4/JiwsjAsXLhAQEEC/fv0yunohhMi6EqLhzz4QuEF7XLgBNP0lW2zmGhAawOVnlzEzNsPcyBxzY3PMMMZ602FqLvkTVYJCkjGs+1RNpe97UcyzGGExiYzcEsjWi8EAFMxjxaSWJSjrnkvPdyM+VulOim7dusXw4cOxs7Nj/Pjx5MyZk5s3bxIQEKBLgi5cuMDjx49RFAUrKytJioQQIr1CLsHaThB+G9TGUGs0VPgmW3SXaRQNX+/9mpjkGN05tycKX25PxT5Ee3zZVcX8empCcqloYePM5oBHjNoSSERcMkZqFT0qF+S7Wp6Ymxjp6S5EVpDupKhDhw589tlnFChQgGLFivH8+XNiY2OxtbXF29ub4sWLs2PHDhYtWkTNmjVxcXHJzLiFECJrUBQ49xvsGAQpCWCTH1otAZdy+o7sg1Gr1NR2q83GWxsxS1L44rQdlY6EodZAvLmazfXsOFTKhMTUJAqZejFzeyIHrgcAUMTRmp9alsQnf/YafC4yR7qTomfPnlG8eHEKFixIaGgogwYNomfPnjg7O+vKLF68mHLlyklCJIQQ6ZEYA9v6wcU12mPPOtBsHlhmv+6fZp7NuLNrPT12aXCICgPAum5dPIYOobSDA4qisOb0A8Ztu8rzxDBMjFR8W8OTr6oWwtRYrefoRVaR7qRo5syZfPXVV9jb2zN37lxmzpxJYGAgkydPxsvLKzNjFEKIrCf0qra77Nl1UBlBzeFQsQ+os9cX/M57O/lp5w+03R3PD1e0A6mf2kChsZNw9m8MwIPwOAZvuMjRW9pkqaSLHT+1LIFXXmu9xS2ypnQnRQ0bNqRhw/9vNtilSxfmzJlDlSpVaNGiBSNHjsyUAIUQIssJWAnbvofkOMjhCC0Xg3slfUf1Qd2Luse8C3OJ3fQnE//SkOPvafZxzWpg0b0tzgUro9Eo/Hb8HpN3Xic+ORVzEzX96xSmS6UCGKmz/lgr8eG98+wzIyMjevXqRYcOHRg5ciRFihRBo9GQmpqakfEJIUTWkRQL2wdCwO/a44LVofkCyJF91tKJT4ln1LFRXDyzg+47kvF+oD2f4uGKx8Rpumn2t0JjGLz+ImfuRwDwSYFcTGpRAvc8VvoKXWQD791OmzNnTmbNmsWRI0eoVasWNWvWZMqUKcTHx2dEfEIIkTU8uQLzq2sTIpUaqg2Fz9Znq4QI4OjdA1gu28qkRdqEKMEEltVQs2lwBSyKFyMpRcPsfTepP/MwZ+5HkMPMmB+bFmdVj/KSEIlMp1IURcnICrdu3Ur//v2JiooiODg4I6vWq+joaGxtbYmKisLGxkbf4QghPhaKAueWwY6B2tllORyhxUIoUFnfkX1QUYlRBO5bh/FPC7EOjgLgXEEVi/zVPLVTUc6xHN96T2Xw+otcC9HuqVm9sD0/NvPB2c5Cn6GLj9zbfH+/dfdZUFAQrq6ur73esGFD/P39+fnnnwF49OhRmhlqQgiRbSQ+hz+/g8t/aI8L1dTOLssGrUNxyXHMODeDPff3EB/+lM/2a6h5Qfs7+Ivd7AN8LGnh1RIv2+IEXHeh+a9H0SiQy8qUkY28aVwyH6pssE6TMBxvnRSVLVuWxo0b06NHD8qVe/U6GnFxcVhZWVG8eHG+/PJLvv322/cOVAghPirBF2BdZwi/ky1nlw08NJCDDw5Q4apCl70a7GK15/eWUrGimppYCxVfFOuIb452DNl4kQfh2sFFzXydGd7Qm1xWpvoLXmRbb50UXb16lfHjx1O3bl1MTEzw8/MjX758mJubExERwZUrVwgMDMTPz4+ffvqJevXqZUbcQghhmP69d5lNfu3sMtdP9B3ZB6MoCoVirai0RkOpu9rWoYe5YX49I665aFt+XK2KcvN6eaYGnAS0G7j+2Kw41QvLBuNCf955TFFCQgLbt2/n8OHD3Lt3j/j4ePLkyYOvry/+/v4UL148o2PVKxlTJIT4T/GRsKUXXP1Te1y4PjT5JVssxqhRNFx8epG/7uxGs2IjtfdHYpoCyUawqYKKPVVsaFy0BbXdavMwxJ5RW67wLCZJt4Frf//C5JANXEUmeJvv7wwfaJ1VSVIkhHijh2fgjy4QGQRqE6g9Bsp/neX3Ljv08BADDw0kNjmWIg8UeuxMxeWZ9tolNxUL/dVUr9SBweUG8yQ6keGbL7PnyhMAPBxyMKmFD2Xcsn7SKPQnUwdaCyGE+AeNBk78AntHgSYFcrpDyyXgXFrfkWW6h88f8s2+b8gRp/DVfg01Lmp/x46yhGU11RwupgKVigr5KrLyVBATt1/jeWIKxmoVPasV4psaHpgZywauwnBIUiSEEO8qNgw2fQ03d2mPvZtC41lgnj02J7W3sOfb4OKUXBOAzd9L0/1zIHVph9I0dunOnO3GnLx7GdBu0TGphQ9FHKXFXRgeSYqEEOJd3D8Gf3SD54/ByAzqTgC/rlm+uywyIZL4lHhS7t0nfsIMKp+9AECQPSyoa8T1/Nr7P9shgAWH7zJk9U2SUjRYmBjR378wnSu6yxYdwmBJUiSEEG9DkwqHp8GBCaCkQm4PaLUUHH30HVmm23pnKyP3D6HpsVSanFAwSYVEY/jjUzVby6lINVJRzaUaLVy+p/HPR3WLMFb2zMP4Zj645LLU8x0I8WaSFAkhRHpFP4YNX8C9w9rjEm2gwTQwy6HfuDLYo5hHLL+yHCsTK3zy+GBvaY+9hT3hhw/w08IUnLTbkRHgYcTvdS1wKVyG6YXb4GtfgRl7btNlcSCKAjktTRje0Jtmvs6yCKP4KLx3UnT48GHmzZvH7du3+eOPP3B2dmb58uUUKFCATz/9NCNiFEII/bu+Azb1hPhwMLGCBlOhVDt9R5UpFlxcwPqb63XHtrEKHfdpqByoHUgdY2tK4dGTKeJfh3Z/Jzv7rj6h/oqjPI5KAKC5rzM/NChK7hxmH/4GhHhH77W06vr16/H398fCwoLz58+TmJgIwPPnzxk/fnyGBCiEEHqVnAA7BsGqttqEyLEEfHkoyyZEAB29OwKg0ijUPqdh+vxUKgcqaIAdZVT8ObY2NnX9UalUhD5P4JuV5+j22xkeRyXgksuCZV3LMa1NKUmIxEfnvdYp8vX1pW/fvnTs2BFra2suXLhAwYIFCQgIoG7duoSEhGRkrHol6xQJkQ09u6ldeyjkkva4fE+oNQqMs/6XffzlQO4MH4T66m0A7jjC/LpG3HFSEfB5AGqVmrVnHjBu21WiE1JQq6B75YJ8V8sTS1MZmSEMxwdbp+j69etUqVLlpfM2NjZERka+T9VCCKE/igIBK2H7AEiOBcvc0HQOePnrO7JMl/r8OUdHfkOeHWdQKwpxZrCqiprdpVUoahVeOb24HxbPkA2XOHk3HIBi+WyY1KIExZ2zx1IEIut6r6TIycmJW7du4e7unub8kSNHKFiw4PtULYQQ+pEQDVv7/n9n+wJVoNl8sHHSb1yZTFEUordu4/GEcdiHRwJw2FvFH/5WjGo4kynOFUlITmH+oTvUnXmYpBQN5iZqvq9dmC6V3DE2yh4b3Yqs7b3+F3/55Zf06dOHkydPolKpePz4MStWrKB///707Nnznet99OgRn332Gblz58bS0pJSpUpx9uxZ3fWYmBh69epF/vz5sbCwoGjRosyZM+c/612/fj3e3t6YmZnh7e3Nxo0b3zlGIUQW9PAszKusTYhURlBjOHy+KUsnROdDzzN1fT9Ot27A4wEDIDySR7lgTDs1mz5zZ1ar5VR0rsj5oAia/HyMaXu06w5V9szDnr5V6VGloCREIst4r5aigQMHEhUVRfXq1UlISKBKlSqYmZnRv39/evXq9U51RkREUKlSJapXr86OHTtwcHDg9u3b2NnZ6cr07duX/fv38/vvv+Pu7s7u3bvp2bMn+fLlo0mTJq+s9/jx47Rp04axY8fSrFkzNm7cSOvWrTly5AiffJJ9dq8WQryCRgPHZsFfY7Vbddi6QstF4FJO35FlquioZ+wY/DmNTqRirIEkY1hfUc2fn6gYV30y9QrUIyYxhVFbAvnt+D0UBXJZmTK8YVGalpJp9iLryZANYePi4rhy5QoajQZvb29y5Hj3NTsGDx7M0aNHOXz48GvLFC9enDZt2jB8+HDduTJlylC/fn3Gjh37yue0adOG6OhoduzYoTtXt25dcubMyapVq14qn5iYqJtNB9qBWi4uLjLQWois5vkT2PQV3P5Le+zdFBrNBAs7fUaV6cL37iZk3I+ogp8CcLaQisV11Dy1+3tF6s/Osv9aOCM3BxIS/fc0+9LODGvgTS4rU73FLcTbepuB1hnS5mlpaYmfnx/lypV7r4QIYMuWLfj5+dGqVSscHBzw9fVlwYIFacp8+umnbNmyhUePHqEoCvv37+fGjRv4+79+EOTx48epU6dOmnP+/v4cO3bsleUnTJiAra2t7uHi4vJe9yWEMEC39sHcStqEyNhCmwy1WpplE6LnSc/pvLQhvzX25kmvPqiCn/LMBn5qoWZSq/8nRC0L9uCbFRf5cvlZQqITcMttyfJu5ZjWupQkRCJLe+95kwkJCVy8eJHQ0FA0Gk2aa40bN37r+u7cucOcOXPo168fQ4cO5dSpU/Tu3RszMzM6dtSunTFr1ix69OhB/vz5MTY2Rq1Ws3DhwjcuFhkSEkLevHnTnMubN+9rlw0YMmQI/fr10x2/aCkSQmQBKYmwbwwc/1l77OCt3dneoYh+48pESnIyRyb3p8+a25gnQ4oatpVT8UclNZ38vqRDLm8K5yzK7ouJTN19ndikJxirVXxZtSDf1vDE3ER2sxdZ33slRTt37qRjx448e/bspWsqlYrU1NS3rlOj0eDn56db/NHX15fAwEDmzJmTJik6ceIEW7Zswc3NjUOHDtGzZ0+cnJyoVavWa+v+d/+3oiiv7RM3MzPDzCzrr0UiRLbz9Aas7/r/tYf8uoH/ODCx0G9cmSRVk8rqlT9gP3cz7n9/VF9xgUX+RjSo9TWnSn6NWqXm8qMoei67xMWHUQCUccvJ+GY+FHa01mP0QnxY75UU9erVi1atWjFixIiXWmHelZOTE97e3mnOFS1alPXrtUvOx8fHM3ToUDZu3EiDBg0AKFGiBAEBAUyZMuW1SZGjo+NLrUKhoaEZFrcQwsApCpxdCjuHQEo8WOSCJr9Akfr6jizDhcaFsvjyYj63b0DAyL6UPvUIgChL+L26moM+KgI6XsBIbURcUgrT91xj8dF7pGoUrM2NGVyvCO3KuqKW3exFNvNeSVFoaCj9+vXL0MSiUqVKXL9+Pc25Gzdu4ObmBkBycjLJycmo1WmHQxkZGb3UffdPFSpUYM+ePfTt21d3bvfu3VSsWDHDYhdCGKi4cNjyLVzbqj0uWA2azs2yU+2nnfgJ1m8n5PBvFEwCjQp2+6pYU0XNpAa/8LNLVQD2Xwtl2KbLPIqMB6BBCSdGNvTGwcZcn+ELoTfvlRS1bNmSAwcOUKhQoYyKh759+1KxYkXGjx9P69atOXXqFPPnz2f+/PmAdrXsqlWrMmDAACwsLHBzc+PgwYMsW7aMadOm6erp2LEjzs7OTJgwAYA+ffpQpUoVJk2aRJMmTdi8eTN79+7lyJEjGRa7EMIA3TkIG7+E58GgNoFaI6H8N6DOmmvrxJ0+Tb2x+3AI1v6SeCMfrG5oQ6N6fThSuDXGamNCoxMY/ecVtl0KBsDZzoIfmxanehEHfYYuhN6915T8uLg4WrVqhb29PT4+PpiYmKS53rt373eqd+vWrQwZMoSbN29SoEAB+vXrR48ePXTXQ0JCGDJkCLt37yY8PBw3Nze++OIL+vbtqxsjVK1aNdzd3Vm6dKnueX/88QfDhg3jzp07FCpUiHHjxtG8efN0xSR7nwnxkUlJgv3j4OhMQIHcntBiIeQrpe/IMkXK06c8+eknorf8CUC0BayorsahZVt6lf4WO3M7NBqFlaeCmLTzGs8TUjBSq+j2aQHZr0xkaW/z/f1eSdHChQv56quvsLCwIHfu3GkGLatUKu7cufOuVRscSYqE+Ig8uwXru0FwgPa4TGfwHw+mVvqMKlMoKSk8/X0ZYbN/gdg4NCrYW0rFqqpqOlX4hq9Lfg3A9ZDnDN14ibP3IwAokd+WCc19KJZP9isTWdsH2xB22LBhjBkzhsGDB780xkcIIT44RYHzv8OOgZAcBxY5ofFsKNpI35FlirizZ7k3Yiiq20EA3HKChf5GhLhY0dSjKR29OxKflMqsv26y4NAdUjQKVqZG9PcvTMcK7hjJQGoh0nivpCgpKYk2bdpIQiSE0L/4CPizD1zZrD12rwzN5oGts37jygQpz54R+tMUojZvRgU8N4c1NUy4XjE/ZZ0/oV+ZflibWrP/WijDN5/kYYR2IHUd77yMblIMJ9usufyAEO/rvZKiTp06sWbNGoYOHZpR8QghxNu7dxQ2fAHRD0FtDDWGQcXeoM5aCw4qKSlErFrN01mz0Dx/DioVZz/JyS+fRNGh/BdM9dXuORkSlcCgdWfZfkm7DImznQWjGhejtrcsQSLEm7xXUpSamsrkyZPZtWsXJUqUeGmg9T9ngwkhRIZLTYYDE+HwVECBXAW1g6mdy+g7sgwVnxLP9b82woxFmN15DMATlxzMrBHPrXzRqFVGlHUsS6pGYdnxe0zdfYOYxP8PpO5T0xMrMxlILcR/ea+fkkuXLuHr6wvA5cuX01yT3ZOFEJnq2U3Y0AMen9ce+34GdSeB2fvtv2hokp+EsqlPI0oFRAMQYw6rqqrZWyoeRa2ilmstvin1DXGxeWj6y1EuPdKuSO3rase4pj5455OJIUKk13vNPstOZPaZEAZCUeDMItg1TLsytbkdNJoBxZrpO7IMpSQlEb58Oc9++RVNXBwaYF8pFZGd6pPT0Q0nKyd88viQz7IgU3ffYNnxe2gUsDE3ZpCsSC2EzgebfSaEEB9UTChs7gU3d2mPC1aDpnPAJp9ew8poMUeO8mTcOJLu3gUgwsOeiZ+GY+1TkhX1pgDavRt3XA6h458HeRKdCECTUvkY1sAbe2vZt1GId/FeSdGYMWPeeH3EiBHvU70QQvzfte3arTrinoGRGdQeDeW+zFIrUyc9fMidMcNRDp0AICaHEcurqTlQPBxFpcIzVTuL7EF4HCM2X2b/9acAuOe2ZGzT4lT2tNdb7EJkBe/VffZiPNELycnJ3L17F2NjYwoVKsS5c+feO0BDId1nQuhJYgzsGgrnftMe5y0OzRdAXu83P+8jkhIfx91fp5GwdDXGyamkqmCHn4p1n6qJN1dhZ2ZHkVxFaOvVgev38jNr300SkjWYGqn5qlohelYrhLlJ1pppJ0RG+WDdZ+fPn3/li3fu3JlmzbJW/74QQg8entEOpg6/A6igYi+oMRyMs0b3kKIoxOzbR+CI/tiGJ2IMXHZTcf6zMriUqMjkXEUpkqsIeS3zcupuOMM2XOZmqHbD7PIFc/FjUx88HLLWwHIh9ClTBlpfvnyZhg0bcu/evYyuWm+kpUiIDyg1RTvN/uAkUFLBxhmazYUCVfQd2Xs59+Qcp0JO4evgi92TWIxmLiX5+GkAnlnDsppqevZdjm/e0rrnhMUkMn77NdafewhALitTfqhflOalnWWWrxDpoPeB1pGRkURFRWVG1UKIrC78jnYhxofaZIHiLaHBFO2WHR+xZ/HP6LSzE+aJCi2OaWhwSsFYA8lGsLW8ERvKQ6KpipoRN/DNWxqNRmH16QdM2nmNqPhkVCpoV86Vgf6FsbM01fftCJElvVdSNGvWrDTHiqIQHBzM8uXLqVu37nsFJoTIZhQFzi+HHYMhORbMbKHBVCjRSt+RZQhbExu+fORF6fVXyBmrPXe2kIqltdQ8yfX/Fp9zT87hY1OPYZsucz4oEgBvJxt+bFac0q4fd2IohKF7r+6zAgUKpDlWq9XY29tTo0YNhgwZgrW19XsHaCik+0yITBT7TLtv2bWt2mO3T7XdZXYu+o0rg8RfvEjIuHEkXLgIQIgdLK2l5pzn/2fOqVVqJlaawdFLdqw8+RiNAjnMjOlX24uOFdwwNso6s+yE+JA+WPfZ3b/X0BBCiHd2fad2qn1sKKhNoOZwqNArS+xblvL0KaHTphO1cSMAqeYm/FUjN7s/MeN+wiNdOUWBirbfMGpNKk+itdt4NCjhxPAG3jjamusldiGyI1m8UQihHwnR2qn255drj+2LaKfaO5XQb1wZQLca9a9z0MRq+8qulHNgZrkwIqyfQcL/y2qScpM37lt2XDMHEnHPbcmYJsWp4iVrDgnxob11UtSvX790l5UNYYUQr3TvKGz6CiKDABVU+EY71d7k428ViTl4kJujhmIeHA7AM3c7fqtjwsnc4ZgbWTC4zHdYGlsSFB3M3gsqAu86c0djhKmxmp7VCvFVVVlzSAh9eeuk6FVrE72KTBUVQrwkOQH+GgvHfwEUsHOFpnPBvZK+I3tviXfv8mTCRGIPHcIciLSCFdXUHPJ5jqJS4W7jztRqU/HK6cXhm0+Z+ddl7oXFAVDZMw9jmhSnQB4r/d6EENncWydF+/fvz4w4hBBZ3eMA2PglPL2mPS7dEfzHg9nHPSEjNSaGZ7/OIXzZMkhJIUUN28qquNLImx7l+1A8LBAjlRHtirQjJt6Yb1aeY9vFYAAcrM0Y0cibBj5O8oukEAbgvccURUZGsmjRIq5evYpKpcLb25uuXbtia2ubEfEJIT52qSlwZDocnAiaFLBygMazofDHvWyHotEQvmE9IVOnoI6IBrRT7JfVVBOcW0WfwnWpnL8ylfNXJjlVw9Kj95ix9waxSamoVdCpojv9anthbW6i5zsRQrzwXknRmTNn8Pf3x8LCgnLlyqEoCtOmTWPcuHHs3r2b0qVL/3clQois69lNbevQo7PaY+8m0GA6WOXWb1zv6ebhrVwbMQCPYFADj3PBbzXVnPf4/7R5F2vtcgIn74QxfPNlbjyJAaC0qx1jmxanWD75xVEIQ/Ne6xRVrlwZDw8PFixYgLGxNr9KSUmhe/fu3Llzh0OHDmVYoPom6xQJ8RY0Gji9APaMhJR4MLeF+lPBpyV8xN1EycHBhE6ZSvS2bQDEmcL6T9Vs91ORapT2vpbX3sSyw8/ZcF479T6XlSmD6xahZZn8qNUf73sgxMfmbb6/3yspsrCw4Pz58xQpUiTN+StXruDn50dcXNy7Vm1wJCkSIp2iHsKmnnD3oPa4YHVo8gvYOus3rvegiYsjbOEiwhYvRklIQFGp+KsErK6iJirH/xOcai7V6FD4M67dc2Dqnhs8T0iR7TmE0LMPtnijjY0NQUFBLyVFDx48yFKrWQsh0kFR4OIa2D4QEqPA2ALqjIWy3T/a1iFFoyF661ZCp04j5ckTAEK8cnO6TXEC7KKJenYpTXkX48qMXZ9E4OMrAPg42/Jj0+KUdLH70KELId7BeyVFbdq0oVu3bkyZMoWKFSuiUqk4cuQIAwYMoF27dhkVoxDC0MWEwta+/9+mI39ZaDYPchfSb1zvIT4ggJAJE3Rbc6ic8jKl/FNOFo6ElKPw7P9lNSmWJD2ty69XrYBobMyNGVC3CO3LuWIkXWVCfDTeKymaMmUKKpWKjh07kpKSAoCJiQlff/01EydOzJAAhRAG7vIG2PY9xIeD2hiqDYFK34HRx7lgfnJICKFTpxH9558AJJqp2VPVhht13DkVFg4olLIvRW232iioOHfTjN3nLEhO0g6yblkmP4PrFSFPDjM93oUQ4l2805iigIAASpUqpTuOi4vj9u3bKIqCh4cHlpaWGRmjQZAxRUL8S+wz2NYPrmzWHjv6QNM52j8/Qpr4eMIWLSJs4SKUhARQqbhYLg8/lwsnMkfa1p4iuYowwnchwzZf5sKDSO05R2vGNi1OWfdceoheCPE6mT7QWq1W4+vrS/fu3Wnfvn22WJNIkiIh/iFwk7Z1KO6ZtnWocn+o/D0Yf3wDiRVFIXrrNkKnTiUlJASAEM9cBHYox7zEvVgYWzCx8kRik2N5HPOYoMhnhDwoy/7AZN1O9t/V8qRzRXfZyV4IA5TpSdHx48dZvHgxa9euJTk5mebNm9OtWzeqV6/+zkEbOkmKhABiw2B7fwjcoD12KAbN5oBTSf3G9Y7iL17kybjxxF+4AMAzWzXLqsOJIird4PAJlSfQsGBDNBqFP84+ZOLOa4THJgHQqGQ+hjUoSl6bj3/PNiGyqg82JT8+Pp61a9eyZMkSDh8+jLu7O127dqVTp07kz5//Xas1SJIUiWzv6p/awdSxT0FlBJX7QZWBH2XrUPLjx4ROn6EbN5RsZsyB6rlY6hOGdY7clMlbhr3399K2SFuGfjKUSw+jGL75MgF/d5V5OuRgdJNiVCyUR493IYRIjw+WFP3T7du3WbJkCcuWLSM4OJjatWuzffv2jKjaIEhSJLKtuHDYMRAurdMe2xeFpr+C88e3Yn1qTCxhCxYQvnQpSmIiqFTE1ipH3yJndOOGuhTrQj+/fiSmJhKfqOKnXddZeSoIRQErUyO+q+VF50rumEhXmRAfBb0kRQAxMTGsWLGCoUOHEhkZSWpqakZVrXeSFIls6dp22PodxDwBlVo7q6zaYDD+uGZWKampRK5fz9NZs0l9pp1Lf6uAOQuqJXPXUZsM5TDJwdzacymWuxhqjFh75gGTdl4jIi4ZgMYl8/GDdJUJ8dH5YIs3vnDw4EEWL17M+vXrMTIyonXr1nTr1i0jqhZC6EN8BOwYDBdXa4/zFNbOLMtfRr9xvYOYo0cJnTSZxBs3ADBxc2Ozvy1LcwWCSoWFsQU2pjY092xOSfuSXHwYyfDNgbpZZV55czC6cXEqFPq492sTQvy3d06KHjx4wNKlS1m6dCl3796lYsWKzJ49m9atW2NlZZWRMQohPqTrO+HPPhATom0dqvgtVBsKJh9XC0nirVs8+eknYg9q92BU29pi/803XKjkyNIj/QAVPUv15OuSXwMQEZvE0I2XWPV3V9mLWWWdKkpXmRDZxTslRbVr12b//v3Y29vTsWNHunbtSuHChTM6NiHEhxQXDruGwoVV2uPcntrWIZey+o3rLaWEhfH055+JXLsOUlPRGKm5VdOTgPqeRJsHcOfyHQBymuWkc7HOaDQKa848YPI/usqalsrH0PpFcZCuMiGylXdKiiwsLFi/fj0NGzbEyMgoo2MSQnxoVzbDtv4QGwqooMI3UGMYmFjoO7J00yQmErF8Oc/mzkMTEwPASS8VK6qrCMl1G57cTlP+29LfciM4kRGbz3LhYRQAhfNaM6ZJMT4pKF1lQmRH75QUbdmyJaPjEELow/Mn2nWHrv79M52nMDT5GVzK6Teut6AoCs937CB06jSSHz3SnixckGU1jdhqe5cyecvQ3PETzIzNMDcyx9zYHBMlF4cv2jHw7FEUBazNjPmuthcdK7hJV5kQ2djHuTmREOL9KIq2m2znEEiI1K5K/WlfqDLgo5pZFnfuPKGTJxMfEACAYp+LvfUcWZjvOopKhYnahEFlB1E0d1EAUlI1rDwVxJRd14lO0CZQzXydGVK/CA7W0lUmRHZnkL8SPXr0iM8++4zcuXNjaWlJqVKlOHv2rO66SqV65eOnn356bZ1Lly595XMSEhI+xC0JYTgiH8CKlrDpa21C5FQSeuzXdpd9JAlR0r17POzdh/vt22sTInMz9tTOw+edoljgfAOV2oh6BeqxpuEaXUJ06m44DWcfYcTmQKITUvB2suGPryowvU0pSYiEEIABthRFRERQqVIlqlevzo4dO3BwcOD27dvY2dnpygQHB6d5zo4dO+jWrRstWrR4Y902NjZcv349zTlzc/kwFNmERgNnFsHeUZAUA0Zm2jWHKvb+aHa0T4mI4NkvvxKxejWkpKCoVPxVAtZWTiHCOhJQUc6xHCMqjMDNxg2AJ9EJjN9+lc0BjwGwtTChv39h2pdzxUitev2LCSGyHYP7JJw0aRIuLi4sWbJEd87d3T1NGUdHxzTHmzdvpnr16hQsWPCNdatUqpee+zqJiYkkJibqjqOjo9P1PCEM0rNbsOVbCDqmPXYprx07lMdTv3GlkyYhgfBlywmbP183iNqsckV6FT7JA3ttYlM0V1E6FutIgwINUKlUJKVoWHz0LrP33SQ2KRWVCtqVc6V/ncLksvr4tiYRQmQ+g+s+27JlC35+frRq1QoHBwd8fX1ZsGDBa8s/efKEbdu2pWuxyJiYGNzc3MifPz8NGzbk/Pnzry07YcIEbG1tdQ8XF5d3uh8h9Co1BY7OhLmVtAmRiRXU+wm67PgoEiJFoyFy0yZu16vP02nT0MTEYOZdFNelS1jZrYAuIQIonbc0DQs2RKVScfDGU+rOOMTEHdeITUqltKsdW775lPHNfCQhEkK8VoZu85ERXnRn9evXj1atWnHq1Cm+++475s2bR8eOHV8qP3nyZCZOnMjjx4/f2BV24sQJbt26hY+PD9HR0cycOZPt27dz4cIFPD1f/nJ4VUuRi4uLbPMhPh4hl2FLL3j8d/JfsDo0mgk53fQbVzrFHjvGk5+mkHj1KgDG+Zww/boLl3xtOfc0gHU31qFCRV6rvHjn8qaXby/MFWfGbr3C7itPAMiTw4wh9YrQzNcZtXSVCZEt6W3vs4xgamqKn58fx44d053r3bs3p0+f5vjx4y+VL1KkCLVr12b27Nlv9ToajYbSpUtTpUoVZs2a9Z/lZe8z8dFISYTDU7UPTQqY24L/eCjVAVSGnxgkXL9B6JQpxB4+DIA6Rw5se3RlTYkYltxcQary/z0VP/f+nIFlB5KQnMqcA7eZe/A2iSkajNQqOld0p08tT2zMTfR1K0IIA/DB9z7LSE5OTnh7e6c5V7RoUdavX/9S2cOHD3P9+nXWrFnz1q+jVqspW7YsN2/efOdYhTA4949pt+h4pt3niyINocFUsE7fWDp9Sn4SytPZs4jasFE7KNzYmJzt2pGn59d8eao/p26cAqBY7mL4Ovjil9ePai7V2Hk5hB+3XeFhRDwAFQvlZlTjYnjltdbn7QghPkIGlxRVqlTppRliN27cwM3t5Sb/RYsWUaZMGUqWLPnWr6MoCgEBAfj4+LxzrEIYjPgI2DMSzv2mPbayh3qToVgzg28dSo2JIXzxYsKWLEWJ1yY21v7+OPTrS2guI36+8RunQrQJkZ2ZHasbajepvfnkOZ0Wn+HILe2u9/lszRnW0Jt6xR1RGfg9CyEMk8ElRX379qVixYqMHz+e1q1bc+rUKebPn8/8+fPTlIuOjmbdunVMnTr1lfV07NgRZ2dnJkyYAMDo0aMpX748np6eREdHM2vWLAICAvjll18y/Z6EyDSKAoEbtDvax4Zqz5XuBLVHg0VO/cb2HzRJSUSuXsOzOXNIjYgAwMLXF4eBA4gv4srkC3NYf3A9KUoKADlMcjCo3CCi4pOZufcmvx2/R6pGwdRYzReVC9KzeiEsTQ3uI00I8RExuE+QsmXLsnHjRoYMGcKYMWMoUKAAM2bMoEOHDmnKrV69GkVRaNeu3SvrCQoKQq3+/+S6yMhIvvjiC0JCQrC1tcXX15dDhw5RrtzHs52BEGlE3Idt38OtPdrjPF7agdRuFfUb139QNBqit23n6cyZJD98CICpuzv2/fpiXbs2V8Kv0OvPljyL17YAlXcqTyuvVlR2rsKWgFBqrDhAWGwSAHW88zKsgTeuuS31dj9CiKzD4AZaGyoZaC0MRmoKnJwD+8dDchwYmULl77XbdBj4itQxR44SOnXq/2eU2duTp1cvVA1rcj/uIfei7jHs6DAAXKxdGF1xNGUdy3L2fgSj/wzk4t8btxayt2Jko2JU8bLX270IIT4OH/VAayHEGzw+D1t6Q8hF7bFbJWg4A+y99BrWf4m/HEjo1CnEHT8BaGeU5ezWlSXFQtn6eDaRf4x96TmTKk8ir5kn/dYEsOG8dp8yazNj+tTypFNFd9m4VQiR4SQpEuJjkBgD+8fBybmgaMDcDuqMhVKfgdpwk4Ok+/d5OnMm0dt3AKAyMcGufXseNvuE6Q/WcPTuUV3ZvJZ5cbdxx83GjeJ5fDl6xZLZ+w7oVqNuVSY/A/yLYG9t2K1hQoiPlyRFQhi66zthe3+IeqA9Lt4S6k6AHA76jesNUp4949mvc4hYuxZSUkClwrZxI/J825sx9+aw+WRvAIxVxoz9dCw1XGpgaaIdF7T/WihjNl/h7rNrAJRysWN042KUdLHT1+0IIbIJSYqEMFTPQ2DHILiySXts5woNpoNnLb2G9SapMbGEL1lC2JIlKHFxAFhVrozD9/0wLezF1fCrbL69WVd+baO1eObUrih/91ksY7de4a9r2ll0shq1EOJDk6RICEOjSYXTi+CvHyExClRGUOEb7Y72plb6ju6VNElJRK5dp51eHxYGgHnx4lj1/opDjhGcDF7EqQunCE8IB0CFiqGfDMUzpycxiSn8/NctFh25Q3KqgomRiq6VCtCrhgfWshq1EOIDkqRICEPy4DRs6/f/gdT5fKHRLHAqod+4XkNJTSVqy588mz2b5MePATBxc8Whb1+s/f1pt60dgXcDdeUtjC341PlTvi75NYVsPVh35gGTd13n6XPtPoNVvewZ0cibQvY59HI/QojsTZIiIQxBbBjsHQnnl2uPzW2h5ggo0wXURvqN7RUUReH5nj08nTmLpNu3gb+n1/f8GruWLVGZmJCiSeFe9D3dcxbUWUAZhzKYGJlwLiiCZiuOcuHvKfbuuS0Z1sCbmkUdZDVqIYTeSFIkhD5pNHBuKewdDQmR2nOlOkCt0ZDDMNfgiT12jNDpM0i4dAkAta0teXp0J2eHDqjMzTn86DC/Bf5GQGgASRrtIottCrehvFN5QqISmLQzkI1/T7HPYWbMtzU86FzJHTNjw0v+hBDZiyRFQujLo3PaFakfn9Me5y2u3bzVtbx+43qN+AsXCJ0+g7gT2rWGVJaW5Or4Obm7dsXIxoZzT87R/8/+PI1/qnuOlYkVn3t/TjfvL/ll/y1+2X+LuH9Mse/vXxgHa3N93ZIQQqQhSZEQH1pcOPw1Fs4sARQwtYYaP0DZHmBkeD+SCTdu8HTmLGL27QP+XmuobVvyfPkFxnny6Mr9dPqnNAnRYv/FlHYozZ4rodSafki3i30Zt5yMbORNifx2H/Q+hBDivxjeJ7AQWZVGAxdWwp4REKedoYVPa+0ijNaO+o3tFZIePuTZ7NlEbflTu/GsWo1t06bYf9MTE2dnXbkUTQpnn5zVdZW98PuF/UwPSuXYbe29OtqYM6R+ERqXzCfjhoQQBkmSIiE+hOCL2gUYH5zUHtsXgfpToEBl/cb1CsmhoYTNnUfEunWQnAyAdZ062PfpjVmhQoB2oPWTuCecCD7BvAvzeBjzUPd8JcUSIhrw57ViaJQwTI3VfFmlIF9Xk13shRCGTT6hhMhMCVHajVtPzdduz2FipV1vqPzXYGRYa/CkREQQvngx4ct/R0lIAMCqUiXsv/sOC5/iANyPvs/MczMJCA1I01Vma2ZLVefqEF2ZrWcUouJTAKhX3JGh9Yvikkt2sRdCGD5JioTIDJpUCFgJ+8ZArHaFZoo1gzrjwNb5zc/9wFKjowlfupTwpb+h+XsVaouSJbHv1w+rT8qRqknlWfwzAH4L/I099/cAYKQywiunF7XcalHQpD6Td97mxpMYAIo4WjOikTcVC+V59YsKIYQBkqRIiIx27yjsHPz/BRhze0D9n6BQDf3G9S+pMbFELF9G2JKlaKKjATDzLop9r2/JUb0aKpUKRVHosL0DgWGBaZ7b2qs1/cv2JzgilXHbrjL+2gUA7CxN+L5OYdqVdcFYdrEXQnxkJCkSIqNE3NMOor7y995eZrZQdSCU+wKMTfUa2j9p4uKIWLmSsIWLSI2MBMDM04M8336Lda1aPIh5yM6bf2BuZE5CasJLCZG1iTVVnPyZuvMuvx2/R3KqgrFaRccK7vSp6YmtpWF1CwohRHpJUiTE+0p8DoenwfFfIDURVGoo0xmq/wBWhtN9pElIIGL1asIWLNTtT2bq7k6eXr2wqVcXlZF28cSe+3pyP/p+mud62HmwsclGUjUKq04F0Xf5DcJjtWOKqhe254cG3ng4yNYcQoiPmyRFQryrF1Ps942BmCfacwWqgP8EcCyu39j+QZOUROQffxA2dx4podrxTSYuLuTp2RPbRg1RGWs/BhRFYdf9XTyJfZLm+cVyF6OHTw+O3XrGmK1XuBbyHIBC9lYMb+hNtcIOH/aGhBAik0hSJMS7uH9MO24oWDuWhpwFwH8cFK4PBrIGj5KcTOTGjTybO5eUx8EAGDs5kefrr7Br1gyVyf+7uRRFYcqZKSy7suyleqo7tWbNQVt2X9EuJ2BrYULfWp50KO+GiYwbEkJkIZIUCfE2Iu5rN24N3Kg9NrP5x7ghM/3G9jclJYWorVt59suvJD94AGg3a8391ZfYtWqF2vT/45sePn/I4suLOfvkLHei7gDQqGAjLIwteBD9jLDHfkzZYE5y6hOM1Co++8SV72p5kdPKcMZICSFERpGkSIj0SIyBI9Ph2Oz/jxsq3RGqDzOYjVuVlBSi/txK2Ny5JN3Xjgkyyp2b3D26k7NtW9Tm2j3GNIqGs0/OsvPuTjbc3ECKol1TSK1SM7TcUFp6tWbdmQdsPn2dZzFJgEJlzzwMb+iNV15rfd2eEEJkOkmKhHgTTSpcWP33uKEQ7Tn3ylB3Ajj66De2vynJyURt+ZNn8+aRHBQEgJGtLbm6dyNXhw6oLdMunDjx1ERWXVulO67gVIF2Rdrhndubu0+MaTT7CFeCtVP0C+axYljDolQv7CBbcwghsjxJioR4FUWB69th31h4elV7LmcBqPMjFGlgEOOGlKQkIjdvJmzefJIfarfZMMqZk1xdu5CzXXuMclilKR8SG8Kue7vSJESL/RdT1rEs98NiGbHhKrsCtYOsrc2N+a6WF5+Xd8PUWMYNCSGyB0mKhPi3u4dh32h4eFp7bG4HlfvBJ18ZxLghJel/7d15fIzX/sDxz0z2dbKJrCJBhIqWULVTSxRV9Npq64beLrbW1SpFF7T3d3vp7YYqbdXSVrmtfbmWqFQ09n0LCZEiiezrzPn9MUyMLKLCJPF9v17zysx5vvPM+U5k5ut5znNOPtdWriJ57lwKEhOB66fJnn8e90EDix0ZAtiTtIfXt71Oal6qqW1b/23YaFyZsfYYC3+Lo0Cv0GrgmRa1GN+lPh4ybkgI8YCRokiIGxL3G0+TndlifGzjaFyjrNVocHCzZM8A46X1aStWcHXefAovGa8ms/LywvOFF3AfOACtg0Ox52TkZ/Dt0W+Zf3A+eqWnlkst9EpPG792rD2Qyb837SUly7i6fbvQGkzu0UDGDQkhHlhSFAlx9TRsfb/oijKtNUQ8B+0mgEtNy/YNMOTlGecZmv8VhUnGcU3WNWrgOeJF49VkJRRDSim2xG9hWvQ00vLSAOgZ0pOpLacSczaD99cc5eSfhwGo6+3M2z2M44aEEOJBJkWReHClJ8L2D2Hvd6D0gAbC+0HHSeARbOneYcjN5doPP5L81VemSRetvb3xHDECt35/M11NdrNzaefYcG4Dq8+u5lz6OQCCdcG88sgr1HFoxd+/O8jWE8aZqN0cbRjXOZRnWtSS+YaEEAIpisSDKDsFfpsNu+dCYa6xLbQbPD6lUsxErc/M4try5aQsWkThFWMBY+3jg+fIEbg9/TRau+Ljmi5kXGDctnEcTzluarO3smdQ2CCG1B/F51vP8fLvUegNxnXKhreqzejHZZ0yIYS4mRRF4sGRnwW/fwG/fQLXTylRqyV0mgpBLS3bN6AwNZXUxd+TsngxhjRj/6z9fPEaORJd375mky7ekJqbyqGrh5i5eyYXMi9grbWmec3mdA/pTnu/x/l5bzJdPv6NtJwCADo3qMmk7mGE1JB1yoQQ4lZSFInqLz8b9n4LUf+CLONpKGo2MhZD9bpY/PL6gj//JGXhIlJ/+AGVnQ2AbVAQniNHoHvySTQlFEM3DFoziIuZFwEIcA5gUbdFeDt6s+XYZZ7+fC9nr2YBEObjwpSeDWldt/IsUCuEEJWNFEWi+sq5Bnu+Mh4dyr5qbHOvbZyFutHToLXsOJr88+dJ/moBaatWoQqMR3LsGjTAa9RIXLp0Ma1afyulFH/8+Qe7L+02FUS2Wlvmd51PcpoD45fu5rfTyQB4OdvyRtf69GsWiJXW8nMrCSFEZSZFkah+Mi/D75/DngWQZ5yZGbda0HosNBkK1padfyf3xAmS584jff16MBgAcGzWDM9RI3Fq06bUmaOVUvx+6Xc+3f8pB68cNLVbaayY3e5rPtmQzI+xB1AKbK21PN86mFc61sHFXsYNCSFEeUhRJKqPa/HG8UL7visaQF2jgXHixYf6gpVl/7ln791H8ty5ZG7fbmpzat8Or5EjcYyIKPO5aXlpvP/7+6w/tx4wDqLuFNSJRh5NOXs+mFELLpGdrwfgyYf9+EdkfQI9ik/iKIQQonRSFImq78oJ42Kth34Eg3FxU/wjoO3rEPqERU+TKaXI2vkbyXPnkv3HH8ZGrRbXbpF4jhiBfYMGt93HvIPzmHtgLvmGfKy11gysP5DnGj3PzuN5fPTfEySlG8dJNanlxuQeDYkIcr+XKQkhRLUlRZGoui7GQtTHcHwNoIxtIR2gzXgIbmfRAdSqoID09RtIXvg1eUevr51mY4Nb76fwfOEFbGvXvv0+lCI5N5nvj31PviGfum51mfLYFPKzgnhhwTEOXTReoebv5sCbT4TRs7GvLNoqhBB3QYoiUbUoBeeijMXQ2a1F7WE9jafJ/Ms+DXWv6TMyuPbDj6QsXmxaikPj4IB7//54PPcsNj4+t91HgaGAtWfX8tWhr0wTMAKMbvQ+czfmsOHI7wA421nzSse6PNe6NvY2JQ/KFkIIUX5SFImqoTAfTqyB6M+KFmrVWBlnoG4zFrxvfxrqXipITCTl2++49uOPGLKMl8FbeXriPvgZ3AcNwtr99qe0EjMTWXV6Fb+c+cV0VZlWo6WmXTBWaU8y8uvzZou2ju0cipez5ReoFUKI6kKKIlG5JZ+Bvd/Avu+LLqu3soOmQ40LtboHWbR7OYcOk7JwIekbNoDeONDZtm4dPJ99Ftcnnyxx9umSxKfHM3jtYK7lXQPAw96DIWHDKbzWii+3xV+ffFHRoX4NJnWXRVuFEOJeqJRF0cWLF5k4cSLr1q0jJyeH0NBQFixYQMT1K3RKGzfx0UcfMWHChFL3u2LFCqZMmcKZM2eoU6cOH3zwAX369LknOYi7UJgHx1dD7CKI21HU7uwDTYbAoyMtulCrMhjI3LadlK+/Lho8DTi2fAzP554zXlZfzsHd6fnpfHf0O7488CUAdXR1GP7Qs1jnNOXjjWeJu3oGgNCazrzdoyHtQ2tUfEJCCCGASlgUpaam0rp1azp27Mi6devw9vbmzJkzuLm5mWIuXR+rccO6det44YUXePrpp0vdb3R0NAMGDOC9996jT58+rFy5kv79+7Nz505atGhxr9IRd+Lqadi7CPYvgezk640a46zTTYdDaCRYWW7OHUNuLmmr/kvKN9+QHxdnbLS2RtejOx7PPluuK8nAOIB67+W9/HzqZzae20iuPte0bcxDH/PF//4kJu4QYJx8cXyX+vRvFoC1LNoqhBD3lEYppSzdiZu9+eab/Pbbb0RFRZX7Ob179yYjI4MtW7aUGjNgwADS09NZt26dqa1bt264u7uzdOnSYvF5eXnk5eWZHqenpxMYGEhaWhqurq7l7pu4jcI8OPar8ajQuZt+5y5+xlNkTYYYJ160oMLkZFKXLCV1yRL0qakAaF1ccB/QH/chQ8o1ePqGSVGT+PXsr2ZtHvYeZGbbEqreIPqkcUoBO2stL7YN5qX2MvmiEELcjfT0dHQ6Xbm+vyvdkaJffvmFyMhI+vXrx/bt2/H39+fll19mxIgRJcb/+eefrFmzhm+++abM/UZHRzNu3DiztsjISGbPnl1i/MyZM5k+ffpfykGUw5WTxrFC+5dAToqxTaOFel0h4lmo28Wiky0qpcg9cICUJUvIWLfetAyHjb8/HsOHoev7NFbOTuXaV25hLrmFuWQVZpkVRH3r9SWyVm92HLLn60PniC40FkR9m/jzRmR9/NwcKj4xIYQQpap0RdHZs2f54osvGD9+PJMmTSImJobRo0djZ2fHsGHDisV/8803uLi40Ldv3zL3m5SURM2a5uNQatasSVJSUonxb731FuPHjzc9vnGkSNyF1PNwaiMcWQnnfytqd/WHpsOMR4V0AZbrH2DIySF97VpSv19C7tGjpnb7xo3xfP45XDp3RmNd/j+bXYm7GLVpVLH2KS2mkZsawasLT5GSlQ/AYyEeTO7RkEb+urtPRAghxB2rdEWRwWCgWbNmzJgxA4AmTZpw5MgRvvjiixKLoq+//prBgwdjb29/233fOkBbKVXqoG07OzvsynnlkCiFvhASdsPJ9cZi6Mrxom0aLYR2u35UqDNoLTvPTn58PKlLl3Ht558xpBknRdTY2uLavTvug5/BITz8zvepz2dP0h7zRqWhqfMwvlzjydkrRwAIqeHEpCca0KmBt0y+KIQQFlTpiiJfX18aNmxo1tagQQNWrFhRLDYqKooTJ06wfPny2+7Xx8en2FGhy5cvFzt6JO5SVjKc3gQnN8CZLZCbVrRNo4XAFsYB0+H9QedvuX4CSq8nMyqK1CVLyIraaZwYEuMpMvdBA9E9/XS55he6ITU3leUnlhOXFsfpa6c5m3aWwhvLjgA61Qj/vNfYfjwVyMLTyZaxXUIZ2DwQGxlELYQQFlfpiqLWrVtz4sQJs7aTJ08SFFR8Ppobl+k//PDDt91vy5Yt2bRpk9m4oo0bN9KqVau77/SDTClIOgSnNsDJjdcnVrxp7L6Du3F8UGgk1HkcHD0s1tUbClNTSfv5Z1KXLqPgwgVTu1Pbtrg/Mwjndu3QWN35kav3fn+PTec3mbV52nvSxONxEuObsfuU4gKp2FlreaFNMC91qIOrDKIWQohKo9IVRePGjaNVq1bMmDGD/v37ExMTw7x585g3b55ZXHp6Oj/++CP/+te/StzPsGHD8Pf3Z+bMmQCMGTOGdu3a8eGHH/LUU0/x3//+l82bN7Nz5857nlO1ohSkXYBLB64fEdoIGYnmMTUbGQdMh3aDgGYWPzV2Q86hQ6R+v4T0tWtR+cZxPFpXV9z69sV90EBsSyi8y+to8lGOJh81a1vRcw2r9mSzYHsceYUGAPpcH0TtL4OohRCi0ql0RVHz5s1ZuXIlb731Fu+++y7BwcHMnj2bwYMHm8UtW7YMpRSDBg0qcT/x8fFob5pAr1WrVixbtozJkyczZcoU6tSpw/Lly2WOotIYDHDtvHEF+ivHi35ePQn5meax1g4Q0t54NKheV4sPlr5ZYWoq6WvWkrZyJblHjpja7Ro2wOOZZ3Dt0QOtw18rUJRSrD+3noWHF3Is5dhN7Vp6es1k4OcnTIOoWwQbB1GHB8ggaiGEqKwq3TxFldWdzHNQpegLIfXc9cLn5uLnFBTmlPwcrTV41oXabY2FUO02YFN5jnyowkIyo6JIW7mKjK1b4frl9BobG1y6dcNj8DPYP/zwXQ1qPp5ynPd/f58DVw4AYKO1oVNgZwKterBytyLuqnH9s5AaTrzZLYwuDWvKIGohhLCAKj1P0QNHKcjLANT1gb7KNODX/HEZPwvzjPvIzzIexcnPhLzrP/OzzLeZ2jMhN914NEifX3LfrGzBKxRq1IcaYUU/PUIsOrN0aXJPniRt5SrSfv0V/dWrpna7hg1w690H1549sPa4+zFNe5L28PLml8nV5+Jg7cDzjZ6noWMPPtmcwI/njZM7ejnbMqazDKIWQoiqRIoiS9MXwCwLz39k7QA1Qs0Lnxph4BZk0QkUy6O002NWHh7onnwSXZ/e2IeFVchrbUvYxraEbaw4VXQl5Jftf2LB9qvMOnQQAHsbLSPahjCqfR2c7Sr3eyeEEMKcfGpb2h2dUtFcj7/lp7U92DqBrbPxp53LTfedr9+/+bFL0X23WqCrBeVcwLQyKO30GNbWuHTsgK5PH5zbtkVjU3FHs5KykhizdQwGZbjeB0cC9SMZ8MVhCvQKjQb6RQQwvkt9fHS3nzNLCCFE5SNFkaVprWHy5esPSil6ZCyKcdmNI0dJX7Pmnp8eu1mBvoAtCVtYfHQxBmUgxLU+ITzPxv0ajuUZAEW70Bq89UQYDXyr0VgzIYR4AElRZGkaDVjLzNklUYWFZMfuJWPzZjK2bKYw8ZJp2704PXaz+PR4fjr1E/89/V9SclNQSoM+vSkXEgZwINM4lquBryuTuofRtl6NCn99IYQQ958URaJSMeTlkfXbLjI2byZz61bTqvQAGgcHnNu2Rdf7qQo/PQZQYChgS/wWfjr5E7sv7Ta1OxdGYLjaiz9T7chB4auz5/Wu9enTxB8rrRzFE0KI6kKKImFx+owMMrfvMBZCO3agsrNN26x0Opw7dsSlS2ecWrX6y3MK3c7VnKsMWD2Ay9nGU5kaNDTWdSMtsSOH4o3jiJztrPl7hzq80CYYe5vKMSGlEEKIiiNFkbCIwqtXydjyPzI2bybr99+LBksD1jVr4tKpEy5du+AYEVHhR4RuFpcWx7vR7xL7Zyzq+vIkrWv2xCq1F2t3p6CUAWuthmda1GJMp3p4OsupTiGEqK6kKBL3hSE7m5wDB8j+I5as6Ghy9u0rmo8JsA0OxqVzZ1y6dMa+USM09+hqOKUUJ1NPEpcex8WMi8zeO9u0raFbBA6ZT7F1pz15hSkA9Aj35Y3I+gR7Od2T/gghhKg8pCgS90Rhaio5sbFkx+4lOzaW3KNHobDQLMa+USNTIWRXp84964tSigNXDrDx/EY2nd9EUlaS+XaDFU7ZkZw834nU7ALAwKO1PXizexhNa7nfs34JIYSoXKQoEndNKUXBxURy9saS/Ucs2bGx5J85UyzO2scHx4gIHJtF4NyhAza+vvelf5/s+4SvDn1leuxg7UCYRxi+Tn5kpoYSfdibPzM0QAF1vZ2Z2C2Mzg28ZVkOIYR4wEhRJO6YIS+P/HPnyNm711QEFSYlFYuzrVsHx6bGIsgxIgIbf38L9BbWnF0DgIe9B++0fIfWfq3Zez6TWeuOc/BCGgDeLnaM6xJKv4gArGVZDiGEeCBJUSRKpE9PJz8+gYKEePLjE8hPiKfgfDz5CQkU/vmn2XggAKytsX+oIY4RzXCMaIpD06ZYu1v+1NPFzItcyjLOb5SSm8KlFA1/33KQrSeuAOBka8VL7evwQttgHG3lz0EIIR5k8i3wgFFKoQoKUNnZGLKzKUhMLCp64hPIT0ig4Px59GlpZe5H6+yMQ+NwHCKMR4EcGjdG6+h4n7Io256kPSw4tIDjKcdJzk0GwFDgSsHVSN45noNB5ZiuKBvdqR5eckWZEEIIpCiyOKXXk/Ldd6AwHn0xHYEx3lc3HpexXeXno3JyMGTnYMgx3lRONoac3OuPs1E3bUOvL1ffrLy8sA0MxLZWIDa1amFbqxa2gcb7Vu7ulXLMzYZzG5i4YyJ6dT1HvQMeuQO5eDGUQr2xv93DfZgQGSZXlAkhhDAjRZGlGQxcnvWhZV7bxgYbb29sg2phE1jLWPwEBmIbFIRtQABap6pRNBQaCtn75162Jmzlp5M/oVd6Imv1QJf7JD/FZHI+x3jVW/Pa7rzVvYFcUSaEEKJEUhRZmlaL65NPGu9ruOnoy00Lwd68KKwG033N9UVjNba2aB0d0Dg4oHVwROtgj9bhpseODkWPHR3ROlx/fA8nRbwfkrKSWHJ8CT+f+pm0POPpPqU0hGgHsGt3SxKvXQOg3vUryjrJFWVCCCHKIEWRhWmsrPD/50eW7kaVkluYy/ht44m6GGVqc7N1p65dL06dacjBZAXk4uNqz/guoTwdESBrlAkhhLgtKYpElVKgL+Drw1+bFURP1xrP8VNhbD2QAihc7K15uUNdnmtdW9YoE0IIUW5SFIlK7Ur2FQ5cOWC6HU0+Sp4+DwBDvieO6QNYdMwbSMHWSsvwVkG80rEubo62lu24EEKIKkeKIlEpLTq8iH/F/qvEbc5aX+zS+5BwsRZZBuMQqz5N/BnfJZQA98oxLYAQQoiqR4oiUSktO7GsWNvkR9/j2Bk/foxJISvfeMl9h/o1mNgtjAa+rve7i0IIIaoZKYpEpaGU4kLmBY4lHyPcK5yLmRevt1vhnB3Jhz85kpJlnIn64QAdE58Io1UdL0t2WQghRDUiRZGwqN2XdrP9wnaOpxznePJxMgoyTNuU0kBGU7TXniIpyxYoIMjTkQmR9ekR7iuX1wshhKhQUhQJi0nLS+PFjS+atdlobajnFoqrvgVHTtYhMcV49ZiXsx1jOtdjYPNAbGTBViGEEPeAFEXCYm6cHrvhpyd/IiPDk//bcIpNZ1MAcLGzZlT7EJ5vIwu2CiGEuLfkW0bcVwZlYN/lfaw+u5qVp1YC4GrryvjGHzJ7bSbrDp8GwNZKy7CWQbzcsS4eTnJ5vRBCiHtPiiJxX5xKPcWvZ35l3bl1JGUlmdrb+TyJbVpv/vH9VfSGNDQa6NskgHFd6snl9UIIIe4rKYrEPbfk2BJmxcxCoQBwtnGmjW8XCpI7sH5nHrkFxivKOjfwZkJkGPV9XCzZXSGEEA8oKYrEPZFdkM32C9tZc3YN2y9sB6Ctf1t6BvchLqEW87afJy0nB4CIIHfefCKM5rU9LNllIYQQDzgpikSFMigD/9zzT1acWkFOYY6p/aXGL+NR0IN3l58iKd04bqietzP/6BZGZ1m9XgghRCUgRZGoMGl5aSw8vJDFxxYDEOAcQLfaT+CS35bF/0vj7JVDAPjp7BnXJZS+TWX1eiGEEJWHFEXirh1LPsayE8tYe3YtufpcAF5+5GUaOf6Nf244wcELiQC4O9rwSse6DHksSFavF0IIUelIUST+sujEaD7b/xkHrhwwtYW6h9LWawg7Y3z58EwMAI62VrzYNoQRbYNxsbexVHeFEEKIMklRJO5Ynj6POXvn8N3R7wCw1lrTpVYX2nj/jbWxVszelQQkY2OlYXCLIF59vC5eznaW7bQQQghxG1IUiXJTSrHx/Eb+Hftv02zUA+oPoHft5/huZzJj1ydgUJjmGhrbuR6BHjLXkBBCiKpBiiJRbrNiZrHk+BIAvB28Gd9kCvtPetN39UHyCw0AdG5QkwmR9WWuISGEEFWOFEXitvQGPctPLGfp8aUAPN/w76hrHXhzcQIZeXEAPFrbg4lP1CciSOYaEkIIUTVVyuXGL168yJAhQ/D09MTR0ZFHHnmE2NhYs5hjx47Rq1cvdDodLi4uPPbYY8THx5e6z0WLFqHRaIrdcnNz73U6VdrxlOMMWjOImTEzMRi01NeO4vsN9fhkSxwZeYU08HVl4XPNWT7qMSmIhBBCVGmV7khRamoqrVu3pmPHjqxbtw5vb2/OnDmDm5ubKebMmTO0adOGF154genTp6PT6Th27Bj29vZl7tvV1ZUTJ06Ytd3uOQ+yqzlXGblxJCm5aVhntYaUHvyRpQXyCfJ0ZHyXUJ5s7IdW5hoSQghRDVS6oujDDz8kMDCQhQsXmtpq165tFvP222/TvXt3PvroI1NbSEjIbfet0Wjw8fGpsL5WV6m5qSw/sZwlx5Zy+Yo/pIwiM8cNgJqudozuVI/+zQKxsaqUBxqFEEKIv6TSfav98ssvNGvWjH79+uHt7U2TJk2YP3++abvBYGDNmjWEhoYSGRmJt7c3LVq0YNWqVbfdd2ZmJkFBQQQEBNCzZ0/27dtXamxeXh7p6elmt+qswFBA1IUoJkVNosuPXZmzcwMXjg8i9+IQcnPccHO04a0nwtg+oSODWwRJQSSEEKLa0SillKU7cbMbp7PGjx9Pv379iImJYezYscydO5dhw4aRlJSEr68vjo6OvP/++3Ts2JH169czadIktm7dSvv27Uvc7++//87p06cJDw8nPT2dOXPmsHbtWg4cOEC9evWKxU+bNo3p06cXa09LS8PV1bVik7awjec28sHuD0jJTUGfE0je5W7os+sA1ydebBPMi+1CcJWJF4UQQlQx6enp6HS6cn1/V7qiyNbWlmbNmrFr1y5T2+jRo9mzZw/R0dEkJibi7+/PoEGDWLJkiSmmV69eODk5sXTp0nK9jsFgoGnTprRr145PPvmk2Pa8vDzy8vJMj9PT0wkMDKxWRZHeoOeTfZ/w9eGv0efWRKX0JCfNWCDaWmkZ/FgtXukoEy8KIYSouu6kKKp0Y4p8fX1p2LChWVuDBg1YsWIFAF5eXlhbW5cYs3PnznK/jlarpXnz5pw6darE7XZ2dtjZVb9iQCnFnqQ9rDi1ggNXDpCQkkPelf4UpjcBNGg18LeIAEZ3qkeAu0y8KIQQ4sFR6Yqi1q1bF7tC7OTJkwQFBQHGI0nNmzcvM6Y8lFLs37+f8PDwu+90FXHk6hFmxMzg4JWDGApcyE9+nILURwHj4qzdw30Y36U+db2dLdtRIYQQwgIqXVE0btw4WrVqxYwZM+jfvz8xMTHMmzePefPmmWImTJjAgAEDaNeunWlM0a+//sq2bdtMMcOGDcPf35+ZM2cCMH36dB577DHq1atHeno6n3zyCfv37+ezzz673yned3qDnoVHFvLZvs8oKLBFn9qT/NRW6PXGwdLtQmswoWt9wgN0Fu6pEEIIYTmVrihq3rw5K1eu5K233uLdd98lODiY2bNnM3jwYFNMnz59+PLLL5k5cyajR4+mfv36rFixgjZt2phi4uPj0WqLrpC6du0aI0eOJCkpCZ1OR5MmTdixYwePPvrofc3vfivQFzAxaiIbz+4gP6U9hmsdKSw0/tojgtx5o2t9WtbxtHAvhRBCCMurdAOtK6s7GahVGWTkZ/Bb4m8sPfoTvx+3JT+5A0rvBEBDX1cmRNanQ/0aaDQy8aIQQojqq0oPtBZ3p0BfwNLjS/niwDxSLoeRf/VxVKHxtFiIlxPju4bSvZGvzEIthBBC3EKKomriQsYF1sWtY9XpXzhzwZ28KyNRBcbTYjVcrJgQ+RB9m/hjLZMuCiGEECWSoqiK0xv0fHf0O/6z71OyrtUl/0ofDPk1AfBytuW1x+sx8NFA7KytLNxTIYQQonKToqgKy8zP5Nn1z3EkQU/elRcx5AYC4GpvzUsd6vBsq9o42sqvWAghhCgP+casgvQGPTsu7GDixrkkJ7ZDn21cDNfR1ooX2gTzYtsQdA6yJIcQQghxJ6QoqkKSspI4dPUQ7/zvKy5faIE+aygA1lYwvGUwf+9QR5bkEEIIIf4iKYqqgCvZV5gZM5P1Jw6Sf7ULhRnDANBqFH2a+vJ6l4b4uTlYuJdCCCFE1SZFUSV1Lfcatla2WGmtGLV2MgdPBlKYPgbQAoreTfwY17k+QZ5Olu6qEEIIUS1IUVTJKKWIvhTN2K1jycqxpeBqZ/Kv9eDG+mTt6rsw+YlHCPWp/BNICiGEEFWJFEWVzKf7P+XLvUvIv9qZgmstQBl/RU1q2/Lek4/SyF/WJxNCCCHuBSmKKpGUzDy27HMh6/Q/QNkCEFwzn7e7N6Zz/ToW7p0QQghRvUlRdKeyssCqhIkQrazA3t48rjRaLTgUDYxOT77GN7+d49vo82TmueCAASeXBCZ1e4S+jzxsFkt2NpS2XJ1GA46Ofy02JwcMhtL77OT012Jzc0Gvr5hYR0djvwHy8qCwsGJiHRyMvxOA/HwoKKiYWHv7on8rdxJbUGCML42dHVhb33lsYaHxvSiNrS3Y2Nx5rF5v/N2VxsbGGH+nsQaD8d9aRcRaWxvfCzD+TWRnV0zsnfzd38VnxB3FymeE8b58Rtx5bHX+jCgvJcolLS1NASrN+BFS/Na9u/kTHB1LjgOl2rc3hZ1ISlfJjrrSY5s1M99vUFDpsQ0bmsc2bFh6bFCQeWyzZqXHenmZx7ZvX3qso6N5bPfupcfe+s/vb38rOzYzsyh2+PCyYy9fLop9+eWyY+PiimLfeKPs2MOHi2KnTi07NiamKPajj8qO3bq1KPbTT8uOXb26KHbhwrJjf/ihKPaHH8qOXbiwKHb16rJjP/20KHbr1rJjP/qoKDYmpuzYqVOLYg8fLjv2jTeKYuPiyo59+eWi2MuXy44dPrwoNjOz7Ni//U2ZKSv2L35GKKWMf4OlxcpnRNFNPiOMN/mMMN6uf0aYvr/T0tTtyEJYFlanhjOyNqsQQghheRqllLJ0J6qC9PR0dDodaYmJuLqWcOXXXRwav3jhCj46B6xKqo7k0HgROTRuJIfG7zxWTp8ZyWfEX4uVzwijKvoZYfr+Tksr+fv7JlIUldOdvKlCCCGEqBzu5PtbTp8JIYQQQiBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBADWlu5AVaGUAoyr7QohhBCiarjxvX3je7wsUhSVU0ZGBgCBgYEW7okQQggh7lRGRgY6na7MGI0qT+kkMBgMJCYm4uLigkajqdB9p6enExgYSEJCAq6urhW678rmQcoVHqx8Jdfq60HK90HKFR6MfJVSZGRk4Ofnh1Zb9qghOVJUTlqtloCAgHv6Gq6urtX2H+WtHqRc4cHKV3Ktvh6kfB+kXKH653u7I0Q3yEBrIYQQQgikKBJCCCGEAKQoqhTs7OyYOnUqdnZ2lu7KPfcg5QoPVr6Sa/X1IOX7IOUKD16+tyMDrYUQQgghkCNFQgghhBCAFEVCCCGEEIAURUIIIYQQgBRFQgghhBCAFEUW9/nnnxMcHIy9vT0RERFERUVZukt3bObMmTRv3hwXFxe8vb3p3bs3J06cMItRSjFt2jT8/PxwcHCgQ4cOHDlyxCwmLy+P1157DS8vL5ycnOjVqxcXLly4n6ncsZkzZ6LRaBg7dqyprbrlevHiRYYMGYKnpyeOjo488sgjxMbGmrZXl3wLCwuZPHkywcHBODg4EBISwrvvvovBYDDFVOVcd+zYwZNPPomfnx8ajYZVq1aZba+o3FJTUxk6dCg6nQ6dTsfQoUO5du3aPc7OXFm5FhQUMHHiRMLDw3FycsLPz49hw4aRmJhoto/qkOutRo0ahUajYfbs2WbtVSXX+0IJi1m2bJmysbFR8+fPV0ePHlVjxoxRTk5O6vz585bu2h2JjIxUCxcuVIcPH1b79+9XPXr0ULVq1VKZmZmmmFmzZikXFxe1YsUKdejQITVgwADl6+ur0tPTTTEvvfSS8vf3V5s2bVJ79+5VHTt2VA8//LAqLCy0RFq3FRMTo2rXrq0aN26sxowZY2qvTrmmpKSooKAg9eyzz6rdu3eruLg4tXnzZnX69GlTTHXJ9/3331eenp5q9erVKi4uTv3444/K2dlZzZ492xRTlXNdu3atevvtt9WKFSsUoFauXGm2vaJy69atm2rUqJHatWuX2rVrl2rUqJHq2bPn/UpTKVV2rteuXVOdO3dWy5cvV8ePH1fR0dGqRYsWKiIiwmwf1SHXm61cuVI9/PDDys/PT/373/8221ZVcr0fpCiyoEcffVS99NJLZm1hYWHqzTfftFCPKsbly5cVoLZv366UUspgMCgfHx81a9YsU0xubq7S6XTqyy+/VEoZP6hsbGzUsmXLTDEXL15UWq1WrV+//v4mUA4ZGRmqXr16atOmTap9+/amoqi65Tpx4kTVpk2bUrdXp3x79Oihnn/+ebO2vn37qiFDhiilqleut355VlRuR48eVYD6/fffTTHR0dEKUMePH7/HWZWsrELhhpiYGAWY/kNa3XK9cOGC8vf3V4cPH1ZBQUFmRVFVzfVekdNnFpKfn09sbCxdu3Y1a+/atSu7du2yUK8qRlpaGgAeHh4AxMXFkZSUZJarnZ0d7du3N+UaGxtLQUGBWYyfnx+NGjWqlO/HK6+8Qo8ePejcubNZe3XL9ZdffqFZs2b069cPb29vmjRpwvz5803bq1O+bdq0YcuWLZw8eRKAAwcOsHPnTrp37w5Ur1xvVVG5RUdHo9PpaNGihSnmscceQ6fTVer809LS0Gg0uLm5AdUrV4PBwNChQ5kwYQIPPfRQse3VKdeKIAvCWsjVq1fR6/XUrFnTrL1mzZokJSVZqFd3TynF+PHjadOmDY0aNQIw5VNSrufPnzfF2Nra4u7uXiymsr0fy5YtY+/evezZs6fYtuqW69mzZ/niiy8YP348kyZNIiYmhtGjR2NnZ8ewYcOqVb4TJ04kLS2NsLAwrKys0Ov1fPDBBwwaNAiofr/bm1VUbklJSXh7exfbv7e3d6XNPzc3lzfffJNnnnnGtCBqdcr1ww8/xNramtGjR5e4vTrlWhGkKLIwjUZj9lgpVaytKnn11Vc5ePAgO3fuLLbtr+Ra2d6PhIQExowZw8aNG7G3ty81rjrkCsb/ZTZr1owZM2YA0KRJE44cOcIXX3zBsGHDTHHVId/ly5ezePFilixZwkMPPcT+/fsZO3Ysfn5+DB8+3BRXHXItTUXkVlJ8Zc2/oKCAgQMHYjAY+Pzzz28bX9VyjY2NZc6cOezdu/eO+1TVcq0ocvrMQry8vLCysipWZV++fLnY/9aqitdee41ffvmFrVu3EhAQYGr38fEBKDNXHx8f8vPzSU1NLTWmMoiNjeXy5ctERERgbW2NtbU127dv55NPPsHa2trU1+qQK4Cvry8NGzY0a2vQoAHx8fFA9frdTpgwgTfffJOBAwcSHh7O0KFDGTduHDNnzgSqV663qqjcfHx8+PPPP4vt/8qVK5Uu/4KCAvr3709cXBybNm0yHSWC6pNrVFQUly9fplatWqbPq/Pnz/P6669Tu3ZtoPrkWlGkKLIQW1tbIiIi2LRpk1n7pk2baNWqlYV69dcopXj11Vf5+eef+d///kdwcLDZ9uDgYHx8fMxyzc/PZ/v27aZcIyIisLGxMYu5dOkShw8frlTvR6dOnTh06BD79+833Zo1a8bgwYPZv38/ISEh1SZXgNatWxebXuHkyZMEBQUB1et3m52djVZr/pFoZWVluiS/OuV6q4rKrWXLlqSlpRETE2OK2b17N2lpaZUq/xsF0alTp9i8eTOenp5m26tLrkOHDuXgwYNmn1d+fn5MmDCBDRs2ANUn1wpzv0d2iyI3LslfsGCBOnr0qBo7dqxycnJS586ds3TX7sjf//53pdPp1LZt29SlS5dMt+zsbFPMrFmzlE6nUz///LM6dOiQGjRoUImX+wYEBKjNmzervXv3qscff7xSXMp8OzdffaZU9co1JiZGWVtbqw8++ECdOnVKff/998rR0VEtXrzYFFNd8h0+fLjy9/c3XZL/888/Ky8vL/WPf/zDFFOVc83IyFD79u1T+/btU4D6+OOP1b59+0xXXFVUbt26dVONGzdW0dHRKjo6WoWHh9/3S7fLyrWgoED16tVLBQQEqP3795t9ZuXl5VWrXEty69VnSlWdXO8HKYos7LPPPlNBQUHK1tZWNW3a1HQZe1UClHhbuHChKcZgMKipU6cqHx8fZWdnp9q1a6cOHTpktp+cnBz16quvKg8PD+Xg4KB69uyp4uPj73M2d+7Woqi65frrr7+qRo0aKTs7OxUWFqbmzZtntr265Juenq7GjBmjatWqpezt7VVISIh6++23zb4oq3KuW7duLfHvdPjw4UqpisstOTlZDR48WLm4uCgXFxc1ePBglZqaep+yNCor17i4uFI/s7Zu3Vqtci1JSUVRVcn1ftAopdT9OCIlhBBCCFGZyZgiIYQQQgikKBJCCCGEAKQoEkIIIYQApCgSQgghhACkKBJCCCGEAKQoEkIIIYQApCgSQgghhACkKBJCCCGEAKQoEkKISm/BggV07dr1jp7z6aef0qtXr3vUIyGqJymKhBDlotFoyrw9++yzlu5ihevQoQNjx461aB/y8vJ45513mDJliqlt2rRpPPLII2ZxUVFRuLm58dprr6GUYsSIEezZs4edO3fe5x4LUXVJUSSEKJdLly6ZbrNnz8bV1dWsbc6cOZbuYrkVFBRUmddbsWIFzs7OtG3bttSYNWvWEBkZyZgxY/jPf/6DRqPBzs6OZ555hv/85z9/+bWFeNBIUSSEKBcfHx/TTafTodFozNp27NhBREQE9vb2hISEMH36dAoLC03P12g0zJ07l549e+Lo6EiDBg2Ijo7m9OnTdOjQAScnJ1q2bMmZM2dMz7lxRGTu3LkEBgbi6OhIv379uHbtmlnfFi5cSIMGDbC3tycsLIzPP//ctO3cuXNoNBp++OEHOnTogL29PYsXLyY5OZlBgwYREBCAo6Mj4eHhLF261PS8Z599lu3btzNnzhzT0bBz586xaNEi3NzczF5/1apVaDSaYv3++uuvCQkJwc7ODqUUaWlpjBw5Em9vb1xdXXn88cc5cOBAme/7smXLyjwNtmTJEvr27cusWbOYPn262bZevXqxatUqcnJyynwNIYSRFEVCiLu2YcMGhgwZwujRozl69Chz585l0aJFfPDBB2Zx7733HsOGDWP//v2EhYXxzDPPMGrUKN566y3++OMPAF599VWz55w+fZoffviBX3/9lfXr17N//35eeeUV0/b58+fz9ttv88EHH3Ds2DFmzJjBlClT+Oabb8z2M3HiREaPHs2xY8eIjIwkNzeXiIgIVq9ezeHDhxk5ciRDhw5l9+7dAMyZM4eWLVsyYsQI09GwwMDAcr8nN/q9YsUK9u/fD0CPHj1ISkpi7dq1xMbG0rRpUzp16kRKSkqp+4mKiqJZs2Ylbvvss8947rnnWLBgAaNHjy62vVmzZhQUFBATE1PufgvxQFNCCHGHFi5cqHQ6nelx27Zt1YwZM8xivvvuO+Xr62t6DKjJkyebHkdHRytALViwwNS2dOlSZW9vb3o8depUZWVlpRISEkxt69atU1qtVl26dEkppVRgYKBasmSJ2Wu/9957qmXLlkoppeLi4hSgZs+efdu8unfvrl5//XXT4/bt26sxY8aUmbtSSq1cuVLd/HE6depUZWNjoy5fvmxq27Jli3J1dVW5ublmz61Tp46aO3duif1JTU1VgNqxY4dZ+9SpU5WtrW2x968k7u7uatGiRWXGCCGMrC1ZkAkhqofY2Fj27NljdmRIr9eTm5tLdnY2jo6OADRu3Ni0vWbNmgCEh4ebteXm5pKeno6rqysAtWrVIiAgwBTTsmVLDAYDJ06cwMrKioSEBF544QVGjBhhiiksLESn05n18dajLXq9nlmzZrF8+XIuXrxIXl4eeXl5ODk53e3bAUBQUBA1atQwPY6NjSUzMxNPT0+zuJycHLNThrduA7C3ty+2LSAgADc3Nz766COeeOIJfH19S9yHg4MD2dnZfzUNIR4oUhQJIe6awWBg+vTp9O3bt9i2m7/QbWxsTPdvjMEpqc1gMJT6WjdiNBqNKW7+/Pm0aNHCLM7Kysrs8a3Fzr/+9S/+/e9/M3v2bMLDw3FycmLs2LHk5+eXniig1WpRSpm1lTSQ+tbXMxgM+Pr6sm3btmKxt45RusHT0xONRkNqamqxbS4uLmzevJmuXbvSoUMHtm7dip+fX7G4lJQUs+JMCFE6KYqEEHetadOmnDhxgrp161b4vuPj40lMTDR94UdHR6PVagkNDaVmzZr4+/tz9uxZBg8efEf7jYqK4qmnnmLIkCGAsWg5deoUDRo0MMXY2tqi1+vNnlejRg0yMjLIysoyFT43xgyVpWnTpiQlJWFtbU3t2rXL1UdbW1saNmzI0aNHS5ynyN3dnc2bNxMZGWkqjPz9/U3bz5w5Q25uLk2aNCnX6wnxoJOB1kKIu/bOO+/w7bffMm3aNI4cOcKxY8dYvnw5kydPvut929vbM3z4cA4cOEBUVBSjR4+mf//++Pj4AMYrvWbOnMmcOXM4efIkhw4dYuHChXz88cdl7rdu3bps2rSJXbt2cezYMUaNGkVSUpJZTO3atdm9ezfnzp3j6tWrGAwGWrRogaOjI5MmTeL06dMsWbKERYsW3TaPzp0707JlS3r37s2GDRs4d+4cu3btYvLkyaZB5iWJjIwsc64hnU7Hxo0b8fLyokOHDly4cMG0LSoqipCQEOrUqXPb/gkhpCgSQlSAyMhIVq9ezaZNm2jevDmPPfYYH3/8MUFBQXe977p169K3b1+6d+9O165dadSokdkl9y+++CJfffUVixYtIjw8nPbt27No0SKCg4PL3O+UKVNo2rSp6SiLj48PvXv3Not54403sLKyomHDhtSoUYP4+Hg8PDxYvHgxa9euNV3GP23atNvmodFoWLt2Le3ateP5558nNDSUgQMHcu7cOdP4qpKMGDGCtWvXkpaWVmqMq6srGzZsoGbNmnTo0IGEhAQAli5dajbWSghRNo269eS4EEJUEtOmTWPVqlXlOj1VnfXv358mTZrw1ltvlfs5hw8fplOnTpw8ebLYoHMhRMnkSJEQQlRy//znP3F2dr6j5yQmJvLtt99KQSTEHZCB1kIIUckFBQXx2muv3dFz7nQBWSGEnD4TQgghhADk9JkQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEAP8PzeQ7Sj7BrFcAAAAASUVORK5CYII=",
      "text/plain": "<Figure size 640x480 with 1 Axes>"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(\n",
    "    temperature_md_lst,\n",
    "    np.array(volume_md_lst) / len(structure_md) * len(structure_opt),\n",
    "    label=\"Molecular Dynamics\",\n",
    "    color=\"C2\",\n",
    ")\n",
    "plt.plot(temperatures_qh_qm, volume_qh_qm, label=\"Quasi-Harmonic (qm)\", color=\"C3\")\n",
    "plt.plot(temperatures_qh_cl, volume_qh_cl, label=\"Quasi-Harmonic (classic)\", color=\"C0\")\n",
    "plt.plot(temperatures_ev, volume_ev, label=\"Moruzzi Model\", color=\"C1\")\n",
    "plt.axhline(structure_opt.get_volume(), linestyle=\"--\", color=\"red\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Temperature (K)\")\n",
    "plt.ylabel(\"Volume ($\\AA^3$)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03887d0e-da24-49ca-9bd8-9452fd666b3c",
   "metadata": {},
   "source": "Both the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) and the quantum mechanical version of the quasi-harmonic approach start at a larger equilibrium volume as they include the zero point vibrations, resulting in an over-prediction of the volume expansion with increasing temperature. The equilibrium volume is indicated by the dashed red line. Finally, the quasi-harmonic approach with the classical harmonic oszillator agrees very well with the thermal expansion calculated from molecular dynamics for this example of using the Morse Pair Potential. "
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7aeb2f2-12c5-492a-a821-384b030b4a68",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}