{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce2efe0c-f86c-4cbd-ab3f-2aae9c5a574d",
   "metadata": {},
   "source": [
    "## Thermal Expansion \n",
    "Calculate the thermal expansion for a Morse Pair potential using the [LAMMPS](https://www.lammps.org/) molecular dynamics\n",
    "simulation code. In the following three methods to calculate the thermal expansion are introduced and compared for a \n",
    "Morse Pair Potential for Aluminium. \n",
    "\n",
    "As a first step the potential is defined for the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code \n",
    "by specifying the `pair_style` and `pair_coeff` commands for the [Morse Pair Potential](https://docs.lammps.org/pair_morse.html)\n",
    "as well as the Aluminium bulk structure: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a96bac5-5afe-4924-90e5-04f6e6b2bedb",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.build import bulk\n",
    "import pandas\n",
    "\n",
    "potential_dataframe = pandas.DataFrame({\n",
    "    \"Config\": [[\n",
    "        \"pair_style morse/smooth/linear 9.0\",\n",
    "        \"pair_coeff * * 0.5 1.8 2.95\"\n",
    "    ]],\n",
    "    \"Filename\": [[]],\n",
    "    \"Model\": [\"Morse\"],\n",
    "    \"Name\": [\"Morse\"],\n",
    "    \"Species\": [[\"Al\"]],\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 \n",
    "the [NIST database for interatomic potentials](https://www.ctcms.nist.gov/potentials). In comparison to just providing\n",
    "the `pair_style` and `pair_coeff` commands, this extended format enables referencing specific files for the interatomic\n",
    "potentials `\"Filename\": [[]],` as well as the atomic species `\"Species\": [[\"Al\"]],` to enable consistency checks if the \n",
    "interatomic 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. \n",
    "While for the Morse potential used in this example this is not necessary, it is essential for extending this example to\n",
    "other interactomic potentials. For the structure optimization the `optimize_positions_and_volume()` function is imported\n",
    "and applied on the `ase.atoms.Atoms` bulk structure for Aluminium:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7bd4ed1a-bffe-40f9-99f8-706797418877",
   "metadata": {},
   "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.\n",
    "This 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": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/janssen/projects/pylammpsmpi/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 \n",
    "volume. 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 \n",
    "The 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).\n",
    "So in analogy to the previous example of calculating the elastic properties from the Equation of State, the `EnergyVolumeCurveWorkflow`\n",
    "is initialized with the default parameters: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b69b6ee2-b526-4913-a6c7-36018e8960af",
   "metadata": {},
   "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\n",
    "then in the second step evaluated with the [LAMMPS](https://www.lammps.org/) molecular dynamics simulation code to derive\n",
    "the equilibrium properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7d1f126e-4fd0-41c5-986b-91d3b5910e3e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'energy': {0.95: -14.609207927145926,\n",
       "  0.96: -14.656740101454448,\n",
       "  0.97: -14.692359030099395,\n",
       "  0.98: -14.716883724875531,\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, \n",
    "    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\n",
    "particular 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": {},
   "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": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'fit_type': 'birchmurnaghan',\n",
       " 'volume_eq': 66.29790137569186,\n",
       " 'energy_eq': -14.735658078942949,\n",
       " 'bulkmodul_eq': 216.05729278060872,\n",
       " 'b_prime_eq': 6.236537173329183,\n",
       " 'least_square_error': array([8.12779273e-07, 2.83453476e-03, 1.45091623e-03, 3.00518393e-05]),\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],\n",
       " 'energy': [-14.609207927145926,\n",
       "  -14.656740101454448,\n",
       "  -14.692359030099395,\n",
       "  -14.716883724875531,\n",
       "  -14.731079276327009,\n",
       "  -14.735659820057942,\n",
       "  -14.731295089579728,\n",
       "  -14.718611862249286,\n",
       "  -14.698196715842329,\n",
       "  -14.670598736769112,\n",
       "  -14.636332030744796]}"
      ]
     },
     "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": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/janssen/projects/atomistics/atomistics/workflows/evcurve/debye.py:10: RuntimeWarning: overflow encountered in exp\n",
      "  return xi**3 / (np.exp(xi) - 1)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "temperatures_ev, volume_ev = workflow_ev.get_thermal_expansion(\n",
    "    output_dict=result_dict, \n",
    "    temperatures=np.arange(1, 1500, 50),\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`\n",
    "which 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 \n",
    "While the [Moruzzi, V. L. et al.](https://link.aps.org/doi/10.1103/PhysRevB.37.790) approach based on the Einstein crystal\n",
    "is 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/)\n",
    "framework. Still the user interface is still structured in the same three steps of (1) generating structures, (2) evaluating \n",
    "these structures and (3) fitting the corresponding model. Starting with the initialization of the `QuasiHarmonicWorkflow`\n",
    "which combines the `PhonopyWorkflow` with the `EnergyVolumeCurveWorkflow`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "493663b9-ea0c-4234-87ef-8f70774794f4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'calc_energy': OrderedDict([(0.9,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[3.9076354064297996, 2.3927365963595605e-16, 2.3927365963595605e-16], [0.0, 3.9076354064297996, 2.3927365963595605e-16], [0.0, 0.0, 3.9076354064297996]])),\n",
       "              (0.92,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[3.9363690516216168, 2.4103308796655574e-16, 2.4103308796655574e-16], [0.0, 3.9363690516216168, 2.4103308796655574e-16], [0.0, 0.0, 3.9363690516216168]])),\n",
       "              (0.94,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[3.9646892271397416, 2.4276719858153393e-16, 2.4276719858153393e-16], [0.0, 3.9646892271397416, 2.4276719858153393e-16], [0.0, 0.0, 3.9646892271397416]])),\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.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",
       "              (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.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.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.06,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[4.12668976922446, 2.5268687084774325e-16, 2.5268687084774325e-16], [0.0, 4.12668976922446, 2.5268687084774325e-16], [0.0, 0.0, 4.12668976922446]])),\n",
       "              (1.08,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[4.1524822498840015, 2.542662047918321e-16, 2.542662047918321e-16], [0.0, 4.1524822498840015, 2.542662047918321e-16], [0.0, 0.0, 4.1524822498840015]])),\n",
       "              (1.1,\n",
       "               Atoms(symbols='Al4', pbc=True, cell=[[4.177958238410259, 2.558261591820219e-16, 2.558261591820219e-16], [0.0, 4.177958238410259, 2.558261591820219e-16], [0.0, 0.0, 4.177958238410259]]))]),\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`\n",
    "and correspondingly also the `QuasiHarmonicWorkflow` require the calculation of the forces `calc_forces` in addition to\n",
    "the 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": {},
   "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 \n",
    "calculate the corresponding energies and forces. The output is not plotted here as the forces for the 108 atom cells \n",
    "result in 3x108 outputs per cell. Still the structure of the `result_dict` again follows the labels of the `structure_dict`\n",
    "as explained before. Finally, in the third step the individual free energy curves at the different temperatures are \n",
    "fitted to determine the equilibrium volume at the given temperature using the `get_thermal_expansion()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "371977fd-cd4e-469f-8955-17a2946c8629",
   "metadata": {},
   "outputs": [],
   "source": [
    "temperatures_qh_qm, volume_qh_qm = workflow_qh.get_thermal_expansion(\n",
    "    output_dict=result_dict, \n",
    "    temperatures=np.arange(1, 1500, 50),\n",
    "    quantum_mechanical=True,\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": {},
   "outputs": [],
   "source": [
    "temperatures_qh_cl, volume_qh_cl = workflow_qh.get_thermal_expansion(\n",
    "    output_dict=result_dict, \n",
    "    temperatures=np.arange(1, 1500, 50),\n",
    "    quantum_mechanical=False,\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\n",
    "Finally, the third and most commonly used method to determine the volume expansion is using a molecular dynamics \n",
    "calculation. While the `atomistics` package already includes a `LangevinWorkflow` at this point we use the [Nose-Hoover\n",
    "thermostat 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": {},
   "outputs": [],
   "source": [
    "from atomistics.calculators import calc_molecular_dynamics_thermal_expansion_with_lammps\n",
    "\n",
    "structure_md = structure_opt.copy().repeat(11)\n",
    "temperature_md_lst, volume_md_lst = 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",
    ")"
   ]
  },
  {
   "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 \n",
    "5K 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 \n",
    "further accelerated by adding the `cores` parameter. The increase in computational cost is on the one hand related to \n",
    "the large number of force and energy calls and on the other hand to the size of the atomistic structure, as these \n",
    "simulations are typically executed with >5000 atoms rather than the 4 or 108 atoms in the other approximations. The \n",
    "volume for the individual temperatures is stored in the `volume_md_lst` list. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eff137a4-61fc-4cbe-8c60-b0dc534a5f3f",
   "metadata": {},
   "source": [
    "### Summary\n",
    "To visually compare the thermal expansion predicted by the three different approximations, the [matplotlib](https://matplotlib.org)\n",
    "is used to plot the volume over the temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "da8f641d-c5e6-4c10-8aeb-c891109e2e6d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Volume ($\\\\AA^3$)')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkUAAAGwCAYAAACnyRH2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACy+klEQVR4nOzdd3yN1x/A8c+9N3sbGRJZJEEIQqi9ia323qOqSmlR1GzN2m2prVSN1qw9au8VK6lNjJCIDNnjPr8/LvfX1KgRbhLf9+uVV/s8z7nnfp+Q3K/zfM85KkVRFIQQQgghPnBqQwcghBBCCJEVSFIkhBBCCIEkRUIIIYQQgCRFQgghhBCAJEVCCCGEEIAkRUIIIYQQgCRFQgghhBAAGBk6gOxCq9Vy7949rK2tUalUhg5HCCGEEK9AURQeP36Ms7MzavXLx4IkKXpF9+7dw9XV1dBhCCGEEOIN3L59m/z587+0jSRFr8ja2hrQfVNtbGwMHI0QQgghXkVsbCyurq76z/GXkaToFT19ZGZjYyNJkRBCCJHNvErpixRaCyGEEEIgSZEQQgghBCBJkRBCCCEEIDVFmS49PZ3U1FRDhyFEjmdiYvKf02uFEOJ1SFKUSRRF4f79+0RHRxs6FCE+CGq1Gk9PT0xMTAwdihAih5CkKJM8TYgcHBywsLCQBR6FeIeeLqYaFhaGm5ub/LwJITKFJEWZID09XZ8Q5cmTx9DhCPFBsLe35969e6SlpWFsbGzocIQQOYA8kM8ET2uILCwsDByJEB+Op4/N0tPTDRyJECKnkKQoE8kQvhDvj/y8CSEymyRFQgghhBBIUiSEEEIIAUhSJDLZ3r17UalU73VpgtGjR1OyZMn39n7ZkUqlYv369YYOQwghsjRJij5wXbp0QaVS0bt372eu9enTB5VKRZcuXd5/YFlEtWrVUKlUqFQqTE1NcXFxoVGjRqxdu9bQob2WsLAw6tWrZ+gwhBA5SUoCXPvL0FFkKkmKBK6urqxcuZLExET9uaSkJFasWIGbm5sBI3t/XrYKec+ePQkLC+Pq1ausWbMGX19f2rRpQ69evd5jhG/HyckJU1NTQ4chhMgpEqNgWVP4tQVc2mroaDKNJEXviKIoJKQmGORLUZTXirVUqVK4ubllGP1Yu3Ytrq6u+Pv7Z2ibnJxMv379cHBwwMzMjEqVKnHixImX9n/48GGqVKmCubk5rq6u9OvXj/j4+Ax9Dh48GFdXV0xNTfH29mbhwoUALFmyBDs7uwz9rV+//qUzj06cOEHt2rXJmzcvtra2VK1aldOnT2doo1Kp+Pnnn2nSpAmWlpZ89913L+zPwsICJycnXF1dKVeuHJMmTWLu3LnMnz+fXbt2AVCjRg369u2b4XWRkZGYmpry11+6f0l5eHgwfvx4unXrhrW1NW5ubsybNy/Da4YMGYKPjw8WFhYUKFCAESNGZEjYnj4qXLRoEW5ublhZWfHpp5+Snp7O5MmTcXJywsHBgXHjxj1zv/98fHbnzh3atGlD7ty5sbS0JCAggGPHjgFw9uxZqlevjrW1NTY2NpQuXZqTJ0++8PsjhPjAxIbB4vpw+yiYWIGZnaEjyjSyeOM7kpiWyEe/fWSQ9z7W7hgWxq+3ZlLXrl1ZvHgx7du3B2DRokV069aNvXv3Zmg3ePBg1qxZwy+//IK7uzuTJ08mMDCQq1evkjt37mf6PX/+PIGBgXz77bcsXLiQiIgI+vbtS9++fVm8eDEAnTp14siRI8yaNYsSJUpw48YNHj58+GY3Dzx+/JjOnTsza9YsAKZOnUr9+vW5cuUK1tbW+najRo1iwoQJTJ8+HY1G81rv0blzZ7788kvWrl1LrVq16NGjB3379mXq1Kn6EZnly5fj7OxM9erV9a+bOnUq3377LcOGDeOPP/7g008/pUqVKhQuXBgAa2trlixZgrOzM+fPn6dnz55YW1szePBgfR/Xrl1j69atbNu2jWvXrtGiRQtu3LiBj48P+/bt4/Dhw3Tr1o2aNWtSrly5Z2KPi4ujatWquLi4sHHjRpycnDh9+jRarRaA9u3b4+/vz5w5c9BoNAQFBcniiEIInYdXdSNEMaFg5QQd1oBTMUNHlWkkKRIAdOzYkaFDh3Lz5k1UKhWHDh1i5cqVGZKi+Ph45syZw5IlS/T1KfPnz2fnzp0sXLiQQYMGPdPv999/T7t27fjiiy8A8Pb2ZtasWVStWpU5c+YQGhrK6tWr2blzJ7Vq1QKgQIECb3UvNWrUyHA8d+5ccuXKxb59+2jYsKH+fLt27ejWrdsbvYdarcbHx4ebN28C0Lx5cz7//HM2bNhAq1atAFi8eLG+Zuup+vXr06dPH0A3KjR9+nT27t2rT4q++eYbfVsPDw++/PJLVq1alSEp0mq1LFq0CGtra3x9falevTqXLl1iy5YtqNVqChUqxKRJk9i7d+9zk6LffvuNiIgITpw4oU9kvby89NdDQ0MZNGiQPiZvb+83+h4JIXKYu6dheQtIiITcBaHjWsjlYeioMpUkRe+IuZE5x9odM9h7v668efPSoEEDfvnlFxRFoUGDBuTNmzdDm2vXrpGamkrFihX154yNjSlbtiwhISHP7ffUqVNcvXqV5cuX688pioJWq+XGjRucP38ejUZD1apVXzvmFwkPD2fkyJH89ddfPHjwgPT0dBISEggNDc3QLiAg4K3eR1EUfcJjampKhw4dWLRoEa1atSIoKIizZ88+M+OrePHi+v9XqVQ4OTkRHh6uP/fHH38wY8YMrl69SlxcHGlpadjY2GTow8PDI8OIl6OjIxqNJsOO8Y6Ojhn6/aegoCD8/f2fO7IHMHDgQHr06MGyZcuoVasWLVu2pGDBgq/2TRFC5EzX98LK9pASB/lKQPs1YGVv6KgynSRF74hKpXrtR1iG1q1bN31dzE8//fTM9ae1Sv+u5/lncvBvWq2WTz75hH79+j1zzc3NjatXr740JrVa/UyN1MuKokE3oy4iIoIZM2bg7u6Oqakp5cuXJyUlJUM7S0vLl/bzMunp6Vy5coUyZcroz/Xo0YOSJUty584dFi1aRM2aNXF3d8/wun8/hlKpVPrHVkePHqVNmzaMGTOGwMBAbG1tWblyJVOnTv3PPl7W77+Zm788aR49ejTt2rVj8+bNbN26lVGjRrFy5UqaNm360tcJIXKoi+tgbS9ITwHPKtB6OZjZ/PfrsiEptBZ6devWJSUlhZSUFAIDA5+57uXlhYmJCQcPHtSfS01N5eTJkxQpUuS5fZYqVYqLFy/i5eX1zJeJiQl+fn5otVr27dv33Nfb29vz+PHjDIXZQUFBL72PAwcO0K9fP+rXr0/RokUxNTV9qxql5/nll1+IioqiefPm+nN+fn4EBAQwf/58fvvtt9d+NHfo0CHc3d0ZPnw4AQEBeHt7c+vWrUyNG3SjVUFBQTx69OiFbXx8fBgwYAA7duygWbNm+vovIcQH5sQC+L2rLiHybQLt/8ixCRFIUiT+QaPREBISQkhIyHMLjy0tLfn0008ZNGgQ27ZtIzg4mJ49e5KQkED37t2f2+eQIUM4cuQIn332GUFBQVy5coWNGzfy+eefA7pHQZ07d6Zbt26sX7+eGzdusHfvXlavXg3ARx99hIWFBcOGDePq1av89ttvLFmy5KX34eXlxbJlywgJCeHYsWO0b9/+P0dHXiYhIYH79+9z584djh07xpAhQ+jduzeffvpphiJq0I0WTZw4kfT09NceWfHy8iI0NJSVK1dy7do1Zs2axbp169447hdp27YtTk5OfPzxxxw6dIjr16+zZs0ajhw5QmJiIn379mXv3r3cunWLQ4cOceLEiRcmvUKIHEpRYO9E2PwloEBAN2ixGIxy9tIekhSJDGxsbJ6pYfmniRMn0rx5czp27EipUqW4evUq27dvJ1euXM9tX7x4cfbt28eVK1eoXLky/v7+jBgxgnz58unbzJkzhxYtWtCnTx8KFy5Mz5499SNDuXPn5tdff2XLli34+fmxYsUKRo8e/dJ7WLRoEVFRUfj7+9OxY0f9EgJvav78+eTLl4+CBQvStGlTgoODWbVqFbNnz36mbdu2bTEyMqJdu3aYmZm91vs0adKEAQMG0LdvX0qWLMnhw4cZMWLEG8f9IiYmJuzYsQMHBwfq16+Pn58fEydORKPRoNFoiIyMpFOnTvj4+NCqVSvq1avHmDFjMj0OIUQWpU2HLV/B3gm646pfQ4NpoH69WbrZkUp53UVtPlCxsbHY2toSExPzTNKQlJTEjRs38PT0fO0PQpGz3L59Gw8PD06cOEGpUqUMHU6OJj93QrwDacmw7hNdHREqqDcZPso+C9U+z8s+v/9NCq2FyASpqamEhYXx9ddfU65cOUmIhBDZT/JjWNVBN9NMbQzN5kKx5v/5spxEkiIhMsGhQ4eoXr06Pj4+/PHHH4YORwghXk/8Q1jeEu6dBmNLaPMrFKzx36/LYSQpEiITVKtW7bW3VxFCiCwhOlS3SnXkVTDPDR3+AJfSho7KICQpEkIIIT5U4SG6hOhxGNi6Qoe1YO9j6KgMRpIiIYQQ4kN06zCsaANJMWBfRLePma2LoaMyKEmKhBBCiA9N8EZY0wPSk8H1I2i7Eiyev/XPh0SSIiGEEOJDcnw+bBkEKFC4ITRfAMZvvsBtTiJJkRBCCPEhUBT461s48GQ/xdJdocHUD2JRxleVJVe0vnv3Lh06dCBPnjxYWFhQsmRJTp06pb/+4MEDunTpgrOzMxYWFtStW5crV668tM8lS5agUqme+UpKSnrXtyNe0ZIlS7CzszN0GAaTWfdfpUoVfvvtt7cP6D989dVXz93oVwiRBaWnwoa+/0+Iqg+HhtMlIfqXLJcURUVFUbFiRYyNjdm6dSvBwcFMnTpV/2GhKAoff/wx169fZ8OGDZw5cwZ3d3dq1aqVYdPQ57GxsSEsLCzDl6yEq1uFuXv37jg7O2NiYoK7uzv9+/cnMjLyvcbRunVrLl++/NI2o0ePpmTJks+cv3nzJiqV6j83i83KXuX+/8umTZu4f/8+bdq0yaSoXmzw4MEsXryYGzduvPP3EkK8hZR4WNEWgn4FlQYazYKqg0GlMnRkWU6We3w2adIkXF1dM+zK7eHhof//K1eucPToUS5cuEDRokUBmD17Ng4ODqxYsYIePXq8sG+VSoWTk9M7iz07un79OuXLl8fHx4cVK1bg6enJxYsXGTRoEFu3buXo0aPkzv1+iu/Mzc3fauPWt5WSkoKJiYnB3j8z7n/WrFl07doVtfrd/3vHwcGBOnXq8PPPPzNp0qR3/n5CiDfwz0UZjcyh5RIoVNfQUWVZWW6kaOPGjQQEBNCyZUscHBzw9/dn/vz5+uvJyckAGUZ4NBoNJiYmHDx48KV9x8XF4e7uTv78+WnYsCFnzpx5Ydvk5GRiY2MzfOVEn332mX6D0KpVq+Lm5ka9evXYtWsXd+/eZfjw4fq2KpWK9evXZ3i9nZ1dhl3rhwwZgo+PDxYWFhQoUIARI0aQmpqqv3727FmqV6+OtbU1NjY2lC5dmpMnTwKZ+/gsPT2d7t274+npibm5OYUKFWLmzJkZ2nTp0oWPP/6YCRMm4OzsjI+Pj37EafXq1VSuXBlzc3PKlCnD5cuXOXHiBAEBAVhZWVG3bl0iIiL0fWm1WsaOHUv+/PkxNTWlZMmSbNu2TX/9ab9r166levXqWFhYUKJECY4cOaJv87z7f/rzYGZmRt68eWnWrNkL7/nhw4fs2rWLxo0bZzh/5coVqlSpgpmZGb6+vuzcuTPDn+Wb3jNA48aNWbFixSv9mQgh3rNHN2BhbV1CZJ4bOv8pCdF/yHJJ0fXr15kzZw7e3t5s376d3r17069fP5YuXQpA4cKFcXd3Z+jQoURFRZGSksLEiRO5f/8+YWFhL+y3cOHCLFmyhI0bN7JixQrMzMyoWLHiC2uRJkyYgK2trf7L1dX1te5DURS0CQkG+XrVlZUfPXrE9u3b6dOnzzMjFE5OTrRv355Vq1a91krN1tbWLFmyhODgYGbOnMn8+fOZPn26/nr79u3Jnz8/J06c4NSpU3z99dcYGxu/cv+vSqvVkj9/flavXk1wcDAjR45k2LBhrF69OkO73bt3ExISws6dO9m0aZP+/KhRo/jmm284ffo0RkZGtG3blsGDBzNz5kwOHDjAtWvXGDlypL79zJkzmTp1KlOmTOHcuXMEBgbSuHHjZ/5+DR8+nK+++oqgoCB8fHxo27YtaWlpz72HzZs306xZMxo0aMCZM2fYvXs3AQEBL7zngwcPYmFhQZEiRTJ8H5o1a4ZGo+Ho0aP8/PPPDBky5Lmvf917Bihbtiy3b9/m1q1bL4xLCGEA94JgYR14dB1s3aD7DnAtY+iosj4lizE2NlbKly+f4dznn3+ulCtXTn988uRJpUSJEgqgaDQaJTAwUKlXr55Sr169V36f9PR0pUSJEsrnn3/+3OtJSUlKTEyM/uv27dsKoMTExDzTNjExUQkODlYSExP/3398vBJcqLBBvtLj41/pe3D06FEFUNatW/fc69OmTVMA5cGDB4qiKM9ta2trqyxevPiF7zF58mSldOnS+mNra2tlyZIlz227ePFixdbW9qUxjxo1SlGr1YqlpWWGLwsLCwVQzpw588LX9unTR2nevLn+uHPnzoqjo6OSnJysP3fjxg0FUBYsWKA/t2LFCgVQdu/erT83YcIEpVChQvpjZ2dnZdy4cRner0yZMkqfPn1e2O/FixcVQAkJCXnu/ZcvX15p3779S78f/zR9+nSlQIECGc5t375d0Wg0yu3bt/Xntm7dmuHP8k3vWVEUJSYmRgGUvXv3vnKcmeV5P3dCCEVRru5WlHHOijLKRlFmV1SU2DBDR2RQT39PPe/z+9+y3EhRvnz58PX1zXCuSJEihIaG6o9Lly5NUFAQ0dHRhIWFsW3bNiIjI/H09Hzl91Gr1ZQpU+aFI0WmpqbY2Nhk+PrQKE9GiF6nzuaPP/6gUqVKODk5YWVlxYgRIzL82Q0cOJAePXpQq1YtJk6cyLVr157bT2hoKFZWVvqv8ePH668VKlSIoKCgDF9btmx5po+ff/6ZgIAA7O3tsbKyYv78+RliAfDz83vu/RUvXlz//46Ojvq2/zwXHh4OQGxsLPfu3aNixYoZ+qhYsSIhISEv7DdfvnwA+n7+LSgoiJo1az732vMkJiY+M3EgJCQENzc38ufPrz9Xvnz5577+de75qacjjAkJCa8cpxDiHTq7SldDlBIHnlWg6xawllraV5XlCq0rVqzIpUuXMpy7fPky7u7uz7S1tbUFdDUTJ0+e5Ntvv33l91EUhaCgoAy/9DOTytycQqdP/XfDd/Ter8LLywuVSkVwcDAff/zxM9f//vtv7O3t9XUuKpXqmUdp/6wXOnr0KG3atGHMmDEEBgZia2vLypUrmTp1qr7N6NGjadeuHZs3b2br1q2MGjWKlStX0rRp0wz9Ojs7Z5hJ9s9ibxMTE7y8vDK0NzLK+Fd59erVDBgwgKlTp1K+fHmsra35/vvvOXbsWIZ2lpaWz/3e/PORnurJDI1/n9NqtRleo/rXTA5FUZ4597x+/93PU69bdJ03b16ioqKeieHf/h3Ty2L7r3t+9OgRAPb29q8VqxAikykKHP4Bdo7QHRdrDh/PASNTw8aVzWS5pGjAgAFUqFCB8ePH06pVK44fP868efOYN2+evs3vv/+Ovb09bm5unD9/nv79+/Pxxx9Tp04dfZtOnTrh4uLChAkTABgzZgzlypXD29ub2NhYZs2aRVBQED/99NM7uQ+VSoXKwuKd9J1Z8uTJQ+3atZk9ezYDBgzI8CF8//59li9fzmeffaY/Z29vn6Fu68qVKxlGCA4dOoS7u3uG4uzn1Zr4+Pjg4+PDgAEDaNu2LYsXL34mKTIyMnom8XkdBw4coEKFCvTp00d/7kWjUm/LxsYGZ2dnDh48SJUqVfTnDx8+TNmyZd+43+LFi7N79266du36Su39/f25f/8+UVFR5MqVCwBfX19CQ0O5d+8ezs7OABmKu9/WhQsXMDY21s8EFUIYgFYLO4bD0dm64/J9ofa38B5moeY0We47VqZMGdatW8eKFSsoVqwY3377LTNmzKB9+/b6NmFhYXTs2JHChQvTr18/Onbs+MwMmNDQ0Awf4NHR0fTq1YsiRYpQp04d7t69y/79+9/qQysn+PHHH0lOTiYwMJD9+/dz+/Zttm3bRu3atfHx8clQWFujRg1+/PFHTp8+zcmTJ+ndu3eGkQQvLy9CQ0NZuXIl165dY9asWaxbt05/PTExkb59+7J3715u3brFoUOHOHHiRIbC4Mzi5eXFyZMn2b59O5cvX2bEiBGcOHEi09/nqUGDBjFp0iRWrVrFpUuX+PrrrwkKCqJ///5v3OeoUaNYsWIFo0aNIiQkhPPnzzN58uQXtvf398fe3p5Dhw7pz9WqVYtChQrRqVMnzp49y4EDBzIkrW/rwIED+hlrQggDSE2CNd3+nxDV+Q4Cx0lC9Iay3EgRQMOGDWnYsOELr/fr1+8/V9Ldu3dvhuPp06dnmAUldLy9vTlx4gSjR4+mVatWhIeHoygKzZo1Y9myZVj8Y7Rr6tSpdO3alSpVquDs7MzMmTMzrDTepEkTBgwYQN++fUlOTqZBgwaMGDGC0aNHA7qlEyIjI+nUqRMPHjzQTzEfM2ZMpt9X7969CQoKonXr1qhUKtq2bUufPn3YunVrpr8X6P5OxsbG8uWXXxIeHo6vry8bN27E29v7jfusVq0av//+O99++y0TJ07ExsYmw0jUv2k0Grp168by5cv1Pz9qtZp169bRvXt3ypYti4eHB7NmzaJu3cyZlrtixYp38ucnhHgFCY9gZTsIPQJqY/h4NhRvZeiosjWV8ryiA/GM2NhYbG1tiYmJeaboOikpiRs3buDp6ZkjVsgeNWoU06ZNY8eOHS8syhVZ04MHDyhatCinTp16bh3eUyqVinXr1j23luxVbd68mUGDBnHu3Llnarreh5z2cyfEa4m6BctbwMPLYGoDrX+FAlUNHVWW9LLP73/LkiNFwrDGjBmDh4cHx44d46OPPnovqyOLzOHo6MjChQsJDQ19aVKUGeLj41m8eLFBEiIhPmj3zsDyVhAfDjYu0P4PcPT979eJ/yS/zcRzvWpxr8h6mjRp8l7ep1UrGaYX4r27vAN+7wypCeBYDNr/DjbOho4qx5CkSIgPkDw1FyIbOrUENg0EJR0KVIdWS8Hsw1tD712SpEgIIYTIyhQF/voODkzRHZdoB41ngSbzt0j60ElSJIQQQmRVaSmw8XM4t1J3XHUIVBsKL1iEVbwdSYqEEEKIrCgpBlZ1hBv7QKWBRjOgVCdDR5WjSVIkhBBCZDUxd3V7mIVfBGNLaPULeNc2dFQ5niRFQgghRFZy/4IuIXp8D6wcod1qcC5p6Kg+CJIUCSGEEFnF9b26R2bJsZC3kG7Kfa53u+aY+D9JioQQQois4OxK2PAZaNPAvSK0WQ7muQwd1QdFlioWWcaSJUuws7MzdBgGk1n3X6VKFX777bf3+p6vokuXLq+8rUh4eDj29vbcvXv33QYlRFagKLBvMqz7RJcQFW0KHdZKQmQAkhQJbt++Tffu3XF2dsbExAR3d3f69+9PZGTke42jdevWXL58+aVtRo8eTcmSJZ85f/PmTVQqFUFBQe8muPfgVe7/v2zatIn79+/Tpk2bTIoq88ycOZMlS5a8UlsHBwc6duzIqFGj3m1QQhhaWgqs7wN7xumOK3wOzReBseznZwiSFH3grl+/TkBAAJcvX2bFihVcvXqVn3/+md27d1O+fHkePXr03mIxNzfHwcHhvb3fv6WkpBjsvSFz7n/WrFl07do1S+5XZ2tr+1qjUl27dmX58uVERUW9u6CEMKTEaPi1GZz9DVRqaDAV6nwHWfDn90Mh3/l3RFEUElLSDPL1Ols4fPbZZ5iYmLBjxw6qVq2Km5sb9erVY9euXdy9e5fhw4fr26pUKtavX5/h9XZ2dhn+9T9kyBB8fHywsLCgQIECjBgxgtTUVP31s2fPUr16daytrbGxsaF06dKcPHkSyNxHOenp6XTv3h1PT0/Mzc0pVKgQM2fOzNDm6eOcCRMm4OzsjI+Pj37EafXq1VSuXBlzc3PKlCnD5cuXOXHiBAEBAVhZWVG3bl0iIiL0fWm1WsaOHUv+/PkxNTWlZMmSbNu2TX/9ab9r166levXqWFhYUKJECY4cOaJv87z737hxIwEBAZiZmZE3b16aNWv2wnt++PAhu3btonHjxhnOR0dH06tXLxwdHTEzM6NYsWJs2rTpuX1cu3aNJk2a4OjoiJWVFWXKlGHXrl0Z2syePRtvb2/MzMxwdHSkRYsW+mt//PEHfn5+mJubkydPHmrVqkV8fHyG7/c/v2eTJk3Cy8sLU1NT3NzcGDdunP66n58fTk5OrFu37oX3LES2FXULFtaBmwfAxEo3w6xMD0NH9cGTQut3JDE1Hd+R2w3y3sFjA7Ew+e8/2kePHrF9+3bGjRuHubl5hmtOTk60b9+eVatWMXv2bFSvuHqqtbU1S5YswdnZmfPnz9OzZ0+sra0ZPHgwAO3bt8ff3585c+ag0WgICgrC2Djzl6rXarXkz5+f1atXkzdvXg4fPkyvXr3Ily9fho1Md+/ejY2NDTt37syQTI4aNYoZM2bg5uZGt27daNu2LTY2NsycORMLCwtatWrFyJEjmTNnDqB7NDR16lTmzp2Lv78/ixYtonHjxly8eBFvb299v8OHD2fKlCl4e3szfPhw2rZty9WrV5+70/zmzZtp1qwZw4cPZ9myZaSkpLB58+YX3vPBgwexsLCgSJEiGb4P9erV4/Hjx/z6668ULFiQ4OBgNBrNc/uIi4ujfv36fPfdd5iZmfHLL7/QqFEjLl26hJubGydPnqRfv34sW7aMChUq8OjRIw4cOABAWFgYbdu2ZfLkyTRt2pTHjx9z4MCBFybpQ4cOZf78+UyfPp1KlSoRFhbG33//naFN2bJlOXDgAN26dXvhfQuR7dw5BStaQ3wEWDtDu1WQr7ihoxJIUvRBu3LlCoqiZPgQ/aciRYoQFRVFRETEKz/W+eabb/T/7+HhwZdffsmqVav0SVFoaCiDBg2icOHCABkShld1/vx5rKysMpz79wevsbExY8aM0R97enpy+PBhVq9enSEpsrS0ZMGCBZiYmAC6ER2Ar776isDAQAD69+9P27Zt2b17NxUrVgSge/fuGUbIpkyZwpAhQ/S1PJMmTWLPnj3MmDGDn376Sd/uq6++okGDBgCMGTOGokWLcvXqVf3345/GjRtHmzZtMtxHiRIlXvh9uXnzJo6Ojhkene3atYvjx48TEhKCj48PAAUKFHhhHyVKlMjwHt999x3r1q1j48aN9O3bl9DQUCwtLWnYsCHW1ta4u7vj7+8P6JKitLQ0mjVrhru7bgqxn5/fc9/n8ePHzJw5kx9//JHOnTsDULBgQSpVqpShnYuLC2fOnHlhvEJkOyF/wpqekJYIjn66hMjWxdBRiSckKXpHzI01BI8NNNh7Z4anicbThOFV/PHHH8yYMYOrV68SFxdHWloaNjb/38V54MCB9OjRg2XLllGrVi1atmxJwYIFn+knNDQUX19f/fGwYcMYNmwYAIUKFWLjxo0Z2t+9e5dq1aplOPfzzz+zYMECbt26RWJiIikpKc8Uafv5+T33/ooX//+/2hwdHfVt/3kuPDwcgNjYWO7du6dPmJ6qWLEiZ8+efWG/+fLlA3QzrZ6XFAUFBdGzZ89nzr9IYmIiZmYZizODgoLInz+/PiH6L/Hx8YwZM4ZNmzZx79490tLSSExMJDQ0FIDatWvj7u5OgQIFqFu3LnXr1qVp06b6x4E1a9bEz8+PwMBA6tSpQ4sWLciV69kZNCEhISQnJ1OzZs2XxmNubk5CQsIrfgeEyMIUBY78BDu+ARTwqg0tF4OptaEjE/8gNUXviEqlwsLEyCBfr/qoy8vLC5VKRXBw8HOv//3339jb2+vrXFQq1TMjMv+sFzp69Cht2rShXr16bNq0iTNnzjB8+PAMBcyjR4/m4sWLNGjQgL/++gtfX9/n1ow4OzsTFBSk/+rdu7f+momJCV5eXhm+no5MPLV69WoGDBhAt27d2LFjB0FBQXTt2vWZYmpLS8vn3vs/H+k9/X7++5xWq83wmn9/3xVFeebc8/r9dz9P/fuR5n/JmzfvM0XJr9vHoEGDWLNmDePGjePAgQMEBQXh5+en/75ZW1tz+vRpVqxYQb58+Rg5ciQlSpQgOjoajUbDzp072bp1K76+vvzwww8UKlSIGzduvPG9PXr0CHt7+9e6ByGynPQ02PIV7BgOKBDQDdqulIQoC5Kk6AOWJ08eateuzezZs0lMTMxw7f79+yxfvpwuXbroz9nb2xMWFqY/vnLlSoZ/xR86dAh3d3eGDx9OQEAA3t7e3Lp165n39fHxYcCAAezYsYNmzZqxePHiZ9oYGRllSHpy5879Wvd24MABKlSoQJ8+ffD398fLy4tr1669Vh+vysbGBmdnZw4ePJjh/OHDh1/4aPJVFC9enN27d79ye39/f+7fv58hMSpevDh37tx55an+Bw4coEuXLjRt2lRf6Pz0keJTRkZG1KpVi8mTJ3Pu3Dlu3rzJX3/9BegSvYoVKzJmzBjOnDmDiYnJc5Neb29vzM3N//P+Lly4oH88J0S2lBwHK9vBiQWASje7rME00MiDmqxI/lQ+cD/++CMVKlQgMDCQ7777Dk9PTy5evMigQYPw8fFh5MiR+rY1atTgxx9/pFy5cmi1WoYMGZJh5MPLy4vQ0FBWrlxJmTJl2Lx5c4YPxMTERAYNGkSLFi3w9PTkzp07nDhxgubNm2f6fXl5ebF06VK2b9+Op6cny5Yt48SJE3h6emb6e4FuhGXUqFEULFiQkiVLsnjxYoKCgli+fPkb9zlq1Chq1qxJwYIFadOmDWlpaWzdulVfn/Vv/v7+2Nvbc+jQIRo2bAhA1apVqVKlCs2bN2fatGl4eXnx999/o1KpqFu37jN9eHl5sXbtWho1aoRKpWLEiBEZRrI2bdrE9evXqVKlCrly5WLLli1otVoKFSrEsWPH2L17N3Xq1MHBwYFjx44RERHx3MTQzMyMIUOGMHjwYExMTKhYsSIRERFcvHiR7t27A5CQkMCpU6cYP378G38PhTCo2HvwWyu4fx6MzKDZPPBtYuioxEvISNEHztvbmxMnTlCgQAFatWqFu7s79erVw8fHh0OHDmUoaJ46dSqurq5UqVKFdu3a8dVXX2FhYaG/3qRJEwYMGEDfvn0pWbIkhw8fZsSIEfrrGo2GyMhIOnXqhI+PD61ataJevXoZCokzS+/evWnWrBmtW7fmo48+IjIykj59+mT6+zzVr18/vvzyS7788kv8/PzYtm0bGzdufKNC8qeqVavG77//zsaNGylZsiQ1atTg2LFjL2yv0Wjo1q3bM4nYmjVrKFOmDG3btsXX15fBgweTnp7+3D6mT59Orly5qFChAo0aNSIwMJBSpUrpr9vZ2bF27Vpq1KhBkSJF+Pnnn1mxYgVFixbFxsaG/fv3U79+fXx8fPjmm2+YOnUq9erVe+57jRgxgi+//JKRI0dSpEgRWrdura/TAtiwYQNubm5Urlz5db5tQmQN98/Dglq6/1rkhS6bJSHKBlTK6yxq8wGLjY3F1taWmJiYDIXDAElJSdy4cQNPT89nCl2zo1GjRjFt2jR27NhB+fLlDR2OeA0PHjygaNGinDp16pk6q+ymbNmyfPHFF7Rr1+6513Paz53IQa7sgt87Q0oc5PV5sqmrh6Gj+mC97PP732SkSDxjzJgxzJo1i2PHjr2wCFhkTY6OjixcuFA/Wyy7Cg8Pp0WLFrRt29bQoQjxek4u1j0yS4kDj8rQfYckRNmIjBS9og9ppEiI7EB+7kSWok2HXaPg8A+64xJtodEsMHr1JU3Eu/E6I0VSaC2EEEK8jZR43YKMl56sOF9tGFQdDK+4PIrIOiQpEkIIId5U7D34rTXcPwcaU/h4Nvi1+O/XiSxJkiIhhBDiTdwLghVt4HGYboZZm9/A7SNDRyXegiRFQgghxOsK2QRre0JqAtgX1u1hJgXV2Z4kRUIIIcSrUhRdMfXOkYACBWtAyyVgZmvoyEQmkKRICCGEeBXpqbB5IJxeqjsO6A71JsuWHTmIrFMkcgSVSsX69esNGsPevXtRqVRER0e/8ms8PDyYMWPGO4tJCJFJEqPg12a6hEilhroTocFUSYhyGEmKPnBdunRBpVJl2IX+qT59+qBSqTJsCptVhYWFvXA7Ccg59ymEMIDIa7CgNtzYDyZWuh3uy30qU+5zIEmKBK6urqxcuZLExET9uaSkJFasWIGbm9tb9a0oCmlpaW8b4n9ycnLC1NT0pW3e5X0KIXKoW4d1e5hFXgEbF+i2DXwCDR2VeEckKRKUKlUKNzc31q5dqz+3du1aXF1d8ff3z9A2OTmZfv364eDggJmZGZUqVeLEiRP6608fIW3fvp2AgABMTU05cOAAXbp04eOPP87Q1xdffEG1atUAuHnzJiqV6pmvp9erVav23Os3b94EXu3xWWbeJ8CWLVvw8fHB3Nyc6tWr62P5p8OHD1OlShXMzc1xdXWlX79+xMfHvzROIUQWEbQCfmkMiY/A2R96/gVOfoaOSrxDkhS9K4qiW+XUEF9vsHNL165dWbx4sf540aJFdOvW7Zl2gwcPZs2aNfzyyy+cPn0aLy8vAgMDefTo0TPtJkyYQEhICMWLF//P93d1dSUsLEz/debMGfLkyUOVKlUAXfLyz+vNmjWjUKFCODo6GuQ+b9++TbNmzahfvz5BQUH06NGDr7/+OkMf58+fJzAwkGbNmnHu3DlWrVrFwYMH6du372vFLIR4z7Ra2P0trO8N2lQo0hi6bAFrJ0NHJt4xqRB7V1ITYLyzYd572D0wsXytl3Ts2JGhQ4fqR2wOHTrEypUr2bt3r75NfHw8c+bMYcmSJfr6nfnz57Nz504WLlzIoEGD9G3Hjh1L7dq1X/n9NRoNTk66XzhJSUl8/PHHlC9fntGjRwOQO3dufdvp06fz119/cezYMczNzQ1yn3PmzKFAgQJMnz4dlUpFoUKFOH/+PJMmTdL38/3339OuXTu++OILALy9vZk1axZVq1Zlzpw5sl+XEFlRSgJs6AMX1+mOKw2EGiNALWMI70pCagI3Ym5wLeYascmxdPDtYLBYJCkSAOTNm5cGDRrwyy+/oCgKDRo0IG/evBnaXLt2jdTUVCpWrKg/Z2xsTNmyZQkJCcnQNiAg4I1j6d69O48fP2bnzp2o//WLaOvWrXz99df8+eef+Pj4vHbfmXWfISEhlCtXDtU/Ci3Lly+foZ9Tp05x9epVli9frj+nKAparZYbN25QpEiR145fCPEOxd6Dle3g3hlQG0OjmeDf3tBR5RixKbFcj77O9ZjrXIu+xrWYa1yPvk5YfJi+jZnGjHZF2qFWGSYJzZJJ0d27dxkyZAhbt24lMTERHx8fFi5cSOnSpQF48OABQ4YMYceOHURHR1OlShV++OEHvL29X9rvmjVrGDFiBNeuXaNgwYKMGzeOpk2bvpubMLbQjdgYgrHFG72sW7du+kc7P/300zPXlSeP5VT/mnGhKMoz5ywtM45UqdVq/eufSk1NfeY9vvvuO7Zt28bx48extrbOcC04OJg2bdowceJE6tSp84p39azMuM9/38vzaLVaPvnkE/r16/fMNSnsFiKLuXtalxA9DgPz3NB6GXhUMnRU2dbt2NtMPDGR0NhQbj++TbqS/tL2uc1yU8C2AAXtCpKUloTFG36Ova0slxRFRUVRsWJFqlevztatW3FwcODatWvY2dkBug+jjz/+GGNjYzZs2ICNjQ3Tpk2jVq1aBAcHP/Nh/NSRI0do3bo13377LU2bNmXdunW0atWKgwcP8tFH72CvGpXqtR9hGVrdunVJSUkBIDDw2dkVXl5emJiYcPDgQdq1awfoEpuTJ0/qHxG9iL29PRcuXMhwLigoCGNjY/3xmjVrGDt2LFu3bqVgwYIZ2kZGRtKoUSOaNWvGgAED3uT29DLjPn19fZ8p7D569GiG41KlSnHx4kW8vLzeKl4hxDt2YS2s7wNpibotO9quhNyeho4q20rTptFqUyviUuP+s21n38509+tOLrNc7yGy/5blkqJJkybh6uqaoRjWw8ND//9Xrlzh6NGjXLhwgaJFiwIwe/ZsHBwcWLFiBT169HhuvzNmzKB27doMHToUgKFDh7Jv3z5mzJjBihUr3t0NZSMajUb/eEij0Txz3dLSkk8//ZRBgwaRO3du3NzcmDx5MgkJCXTv3v2lfdeoUYPvv/+epUuXUr58eX799VcuXLign/V14cIFOnXqxJAhQyhatCj3798HwMTEhNy5c9OsWTPMzc0ZPXq0/hrokq3nxfqu77N3795MnTqVgQMH8sknn3Dq1CmWLFmSoZ8hQ4ZQrlw5PvvsM3r27ImlpSUhISHs3LmTH3744bViFkK8A4oC+ybB3gm6Y6/a0GIRmNkYNq5szkhtRMtCLVl8YfF/tq3qWjXLJESQBZOijRs3EhgYSMuWLdm3bx8uLi706dOHnj17Arqp0kCGIlWNRqP/l/2LkqIjR448M8IQGBj4wtWEk5OT9e8FEBsb+za3lW3Y2Lz8l8HEiRPRarV07NiRx48fExAQwPbt28mV6+V/qQMDAxkxYgSDBw8mKSmJbt260alTJ86fPw/AyZMnSUhI4LvvvuO7777Tv65q1ars3buX/fv3AxkTZIAbN248c+593Kebmxtr1qxhwIABzJ49m7JlyzJ+/PgMM9mKFy/Ovn37GD58OJUrV0ZRFAoWLEjr1q1fO14hRCZLTdSNDl18skRHuc+gzregfr1/ZImMFEXhctRlTDWmaFQa0pV0nCydKJSrEK7WrrjZuOFmrftysnLCWG38352+RyrlVYoj3qOnyc7AgQNp2bIlx48f54svvmDu3Ll06tSJ1NRUvL29KVu2LHPnzsXS0pJp06YxdOhQ6tSpw/bt25/br4mJCUuWLNE/DgH47bff6Nq1a4bk56nRo0czZsyYZ87HxMQ884GalJTEjRs38PT0lBlFQrwn8nMn3lhsGKxs+6Sg2ggaTodSnQwdVbYVmRjJkbAjHL57mMP3DhOZFKm/1qBAA8ZUGIOp5uWL675LsbGx2NraPvfz+9+y3EiRVqslICCA8ePHA+Dv78/FixeZM2cOnTp1wtjYmDVr1tC9e3dy586NRqOhVq1aL93i4alXKRB+aujQoQwcOFB/HBsbi6ur61vcmRBCCIO7dwZWtJWC6kxw6dElWm1ojlb9/M/Rbz76hlaFWr3wczYrynILL+TLlw9fX98M54oUKUJoaKj+uHTp0gQFBREdHU1YWBjbtm0jMjIST88XF8Y5OTllqEUBCA8Pf+Hif6amptjY2GT4EkIIkY1dXAeL6ukSIvvCuhWqJSF6I4qicH3xbKbOT8cy8fkPnOp61s1WCRFkwaSoYsWKXLp0KcO5y5cv4+7u/kxbW1tb7O3tuXLlCidPnqRJkyYv7Ld8+fLs3Lkzw7kdO3ZQoUKFzAlcCCFE1qQosHci/N5FN8PMqzZ03yEzzN7Qg9uXONmxKR7zd+DyCGoF/T8pcrV2pbNvZ3a22Imtqa0Bo3wzWe7x2YABA6hQoQLjx4+nVatWHD9+nHnz5jFv3jx9m99//x17e3vc3Nw4f/48/fv35+OPP86wdk2nTp1wcXFhwgTdrIL+/ftTpUoVJk2aRJMmTdiwYQO7du3i4MGD7/0ehRBCvCdSUJ1prkZdZcm8vjRcdQvbBEjRwPLqag6Ut6KZZ11GlhuJJpt/X7NcUlSmTBnWrVvH0KFDGTt2LJ6ensyYMYP27f+/qmhYWBgDBw7kwYMH5MuXj06dOjFixIgM/YSGhmZYDblChQqsXLmSb775hhEjRlCwYEFWrVqVqWsUZbGadSFyNPl5E/9JCqozjTYxkUvDBtJ+zy0AHuQz43r/xrT4qCHjHEtmuVlkbyrLzT7Lql5WvZ6ens7ly5dxcHAgT548BopQiA9LTEwM9+7dw8vLK8MioEIAUlCdSRRFIfjwn6SPmorxnXAANpVR0WHmNo7d0rDg4HVW9iqPlWmWG2PRy9azz7IjjUaDnZ0d4eG6vzAWFhbZrrhMiOxEq9USERGBhYUFRkbya0z8y4U1sP4zWaH6LSnp6fw5tjuevx/DWAuPrOCnhmrOuZlz+I8znLmu+9n75fBNPqueM1bul98mmeTpDu9PEyMhxLulVqtxc3OTf4CI/9NqYc84ODBFd+xVG1osBLPsV/D7PmkVLdHJ0cSnxBOXGkdcahwJt29hPXER3hdvABBU1JxLPapjmWqP8bminIlTo1Gr6F/Tm0+qFDDwHWQeSYoyiUqlIl++fDg4ODx3o1MhROYyMTHJUDcoPnDJj2FtL7i0RXdcoR/UGi0F1f8hXZtOm81t+PvR3/pzlS5o6bFDi3kyJJrAotpq8jRtjG1iO5buv4JWAdfc5sxo7U9p96yzRUdmkKQok2k0mtfei0sIIcRbeHRDVz8UEQIaU2j8A5SQ7XRehVqlRoVutNUiSaHHdi2VgnWlxpdc4MdGGqJsXTA7UoK7kVcAaFbKhTGNi2JtlvNq+SQpEkIIkX3d2A+rO0FiFFg5QZvfIH9pQ0eVbahUKlY0WMHVPRtIGj0Zk4gY0tWwpoKaNRVUpMSVIulKE9CagDoRM6d1jGz8U45MiECSIiGEENmRosCJBbB1CCjp4FwK2iwHG2dDR5atKCkpRP7wI9oFCzBRFBIcbfitlQPbze+TdP9j0mJLAqAxv4GZyyrUxtGExYdhZ2Zn0LjfFUmKhBBCZC9pKbB1EJxaojsu3hoazQRjc4OGlR2kadOISooiIjGCqJBzWIyfh9n1MAB2l1DxS8144tIh6Xp/lLRcqFQKLT8yZ3jdztyJr05KegpF8hQx8F28O5IUCSGEyD7iH8KqjhB6GFBB7TG6omqZhQjoZpIdvXeUG7E3iEyMJCIxgojECN3/J0QQlRyFok2n3gmFdnu1mKTDYzOYW1/NMR8jUiJqkvqoBoqiwj2PBTNal8TfTVdMbWtW1MB39+5JUiSEECJ7uH8eVrSDmFAwtYHmC8Gnzn+/7gOy7/Y++u3p98LreWIU+mzW4ndLV0x9uoCKK71r425XiMvHXbgdqUsLWpTOz+jGRbP0oozvwod1t0IIIbKn4A2wrjekJkDuAroFGe0LGTqqLKdY3mJUc63GzZibRCRGEJ8ar7ugKFS6qNB9hxbLZEgyhmU11OwoqcIl0o47JwsQn5KOjZkR45v50bD4h1mbJdt8vKLXWSZcCCFEJtFqYf9k2Kvb3JsC1aHlYjDPWevjvCsPEx+y4vg8LGb8SoUQ3cf9ZWfdVPswWwuSwpqS9rg4AGU9czO9dUlc7HJWbZZs8yGEECL7S46D9Z9CyEbdcbnPoPZY0MhH16t4mPiQoVPq0WVDHLnjIF0Ff1RSs62KBZPLbWbQ7+eIe5yMkVrFgNo+9K5aEI36w67Nkr9ZQgghsp6om7CyPTy4ABoT3Q73/h0MHVW2cDfuLhdun8Ji/h8M3BKnO5cbfmiswcTXj6KPO9Jp4QkACthbMqN1SYrntzNgxFmHJEVCCCGylut74fcuugUZLR2g9a/g9pGho8ryElITWHB+Aft3LOLTjcnkfaQ7v7W0iuXV1SSmO2N2oQWPYnWP0TqWc2dY/SKYm8guDE9JUiSEECJrUBQ4Ogd2fPP/BRlb/wq2LoaOLMtLSE2g1frmlN1+i9GHFDQKPLY1YVu7glz3sSUhxJKkiDokKBryWpkwuUVxahR2NHTYWY4kRUIIIQwvNRH+/ALOrdQdl2ine2RmbGbQsLKLOxeP8cnsm3jp1mHEul49fEaNJD+mfLn6LInhkQDU9nVkYjM/8liZGjDarEuSIiGEEIYVcwdWdYB7Z0ClgcBx8FFvWZDxPyiKQnj8A+78sgCTuSvxSoF4U1gQqCay8m3a3ohl1IYQHielYW6sYVQjX1qXcUUl39cXkqRICCGE4dw6rNvQNT4CzHNDyyVQoKqho8rSzkacpcOWDuSJUfh0i5biN3U1Quc8VMxpoOahpTmqkPIMPHQegJKudkxvXRLPvJaGDDtbkKRICCHE+6cocHIRbB0M2jRw9NNt6JrL3dCRZXlrL6+h6jktXXdpsUiGZCNYXl3N9tIqUhMKknS9FUqaHZBOyUL3+aNTPYw0akOHnS1IUiSEEOL9SkuGLYPg9C+646LNoMmPYCIjGS/zIP4BP/41Du95u2hz5f8LMf7UUMO9XEYkh9ch9VFlQI3K+CHmLquIt05Fo+5p2MCzEUmKhBBCvD+P7+s2dL1zHFBBrVFQ8QupH3oJRVEYvH8wUdu20HObFptESFPD9lp2nK7pjoPiRnhwaeJirQAo65NCwzJmmJl0oa5nXakheg2SFAkhhHg/7pzUFVQ/DgMzW2i+CLxrGTqqLOtK1BVmnZ7FiSt76LZDS+Vg3ejQTQcwHT2YAdW6sPDgdaZsv0xKupbcliZMbOZHnaJOBo48+5KkSAghxLt35lfYNADSU8C+MLT5DfIUNHRUWVZEQgS9d/Ym34X7TN2iJXccaFWwvryK6HaB9C7elLbzj3L8hm6FxhqFHZjY3A8Ha1nC4G1IUiSEEOLdSU+F7cPh+FzdceGG0PRnMLU2bFxZWKo2lTa/f0yTbVHUPqMbHTLycCP/xEkMK1GC30/dockPx4hLTsPCRMOIhr60kan2mUKSIiGEEO9GXDj83hVuHdQdVxsKVQaDWmZCPU9wZDA/nPmBuBPHGbEhAcdo3flcHTviMHAAj9LV9Ft2ip3BDwAIcM/F1FYlcM8jBeqZRZIiIYQQme/OqSf1Q/fAxAqazYPCDQwdVZa29MwCPJfuo8FxBTUQn8cCr8kzsKtYmR0X7zN07Xki41Mw1qgYWLsQvaoU+OB3tc9skhQJIYTIXKd+gS1f6eqH8njr6ofsfQwdVZalKApn9v1OzdHbyP9Q97jMullTfIYNI8HIlMF/nGX1yTsAFHayZlqrkvg62xgy5BxLkiIhhBCZIy1ZtxjjqSW648IN4eM5YCYf4M+TkJrA5kvriZg9myp7IsmvQLQlLG5ozk8jx3DyVixf/n6cO1GJqFTQq3IBBtbxwdRIdrV/VyQpEkII8fZi7+m267hzAlBBjeFQ6UupH/oXRVE4G3GW9VfXE3JkM13Xx1EsQnftVlk30r/oytAClfh++1XmH7iOokD+XOZMa1WSsp65DRv8B0CSIiGEEG/n1mFY3Rniw5+sP7QQvGsbOqosJTQ2lGXBy1h7ZS3pqck0O6TlmyMKRlpIsjJB+1VP6rTuw99hcfReHMSlB48BaB3gyohGvliZysf1+yDfZSGEEG9GUeD4PNg+TLd/mUNRaPMr5C5g6MiylOsx12myvgkA7g8UPtuUjke47tqRwioW1kknJmkefsstORniSGq6Ql4rEyY0K05tX0cDRv7hkaRICCHE60tN1C3GeHaF7rhYc2j8g+xf9sTDxIf02N6DuNQ4EtMS0aQrfHxEofkhLUZaiDWHhYFqjhRRo03JQ+KtVhxJdAAU6vg6MqGZH3msTA19Gx8cSYqEEEK8nuhQ3XT7sLOg0kDtsVD+M9m/7B9OPjjJtZhrALhGKAzblE6B+7prx31UzK+rJtpCTeqjcqRFNESr1WBtasTIRr60KJ1fFmI0EEmKhBBCvLpre+CPbpD4CCzyQMsl4FnF0FFlOYmpiai1Co2PKrQ9rEaVCmlWZpzvWI7LpXJRMDqRiyH+xD3KA0BFrzxMblECFztzA0f+YZOkSAghxH9TFDg8C3aNBkUL+UpC61/BztXQkWVJ4cGn+G5pOl5hAFqsqlfHacxoitnb8/upO2zYFUxcchrmxhqG1i9Mh4/cUctCjAYnSZEQQoiXS46DjX3h4jrdccn20GAqGMuoxr8p6elELl5MuelrMEmHeFPwHPUduZs2I+JxMkN/Ocnuv3VV1qXdczGlZQk880odVlYhSZEQQogXe3hFVz8U8TeojaHeRAjoLvVDz5F8/Tr3hg4l6ew5TIDTBVTMra9mV5NGbD4fxjfrLxCdkIqJRs2XdXzoUVm26chqJCkSQgjxfCF/wrpPIeUxWDlBq1/ArZyho8pylLQ0IhcvJnzWLFSpaSSYwi811ewprkJJt2LgqgtsOhcGQFFnG6a1KkkhJ2sDRy2eJ0suNXr37l06dOhAnjx5sLCwoGTJkpw6dUp/PS4ujr59+5I/f37Mzc0pUqQIc+bMeWmfS5YsQaVSPfOVlJT0rm9HCCGyl/Q02DlSN0KU8hjcK8In+yUh+oeIhAj23d7H0j+/Y1+DikRMnYYqNY0gTxVf9tCwp4SatLgiaO6OYNO5MDRqFf1rerP+s4qSEGVhWW6kKCoqiooVK1K9enW2bt2Kg4MD165dw87OTt9mwIAB7Nmzh19//RUPDw927NhBnz59cHZ2pkmTJi/s28bGhkuXLmU4Z2Zm9q5uRQghsp+4CPijK9w8oDsu3xdqjQaNsUHDMrSktCSWhyzndPhpgiODiYqLoMlRhRYHdesOxZvCklpq9vmpMMIG25iO3LnnRiLg5WDFtFYlKJ7fztC3If5DlkuKJk2ahKurK4sXL9af8/DwyNDmyJEjdO7cmWrVqgHQq1cv5s6dy8mTJ1+aFKlUKpycnN5F2EIIkf3dPqHbv+zxPTC2hCY/QrFmho4qS9h/Zz8zTs8AdKtSf7U5nQIPdNdOeqnY1sqDbtW+pF1KCQb9cY470bpNXHtWLsDA2j6YGcsmrtlBlnt8tnHjRgICAmjZsiUODg74+/szf/78DG0qVarExo0buXv3LoqisGfPHi5fvkxgYOBL+46Li8Pd3Z38+fPTsGFDzpw588K2ycnJxMbGZvgSQogcSVHg+HxYXE+XEOX1gV57JCH6hwrOFehaqCNdjpoxYYkuIXpsBrMaqZncQs3Z9Af0X32YdguOcTc6EbfcFqzqVZ5h9YtIQpSNqBRFUQwdxD89fZw1cOBAWrZsyfHjx/niiy+YO3cunTp1AiAlJYWePXuydOlSjIyMUKvVLFiwgI4dO76w36NHj3L16lX8/PyIjY1l5syZbNmyhbNnz+Lt7f1M+9GjRzNmzJhnzsfExGBjY5NJdyuEEAaWkqDbruPcSt2xbxNo8hOYSt0LQEJqAl/u+5K7Jw/w6Zb/71l2zEfFgkA1MVYqrNNLEX27IQmJFgC0/8iNYfWLYCmbuGYJsbGx2NravtLnd5ZLikxMTAgICODw4cP6c/369ePEiRMcOXIEgClTpjB//nymTJmCu7s7+/fvZ+jQoaxbt45atWq90vtotVpKlSpFlSpVmDVr1jPXk5OTSU5O1h/Hxsbi6uoqSZEQIueIvKZ7XPbgwpPtOsboaog+4On29+Luce7hOZwsnLC3sOf8vTOcmjiYJkcUNMr/9ywr1rIX3nbF2HvWhlXHdZmSi505k5oXp5J3XgPfhfin10mKslwamy9fPnx9fTOcK1KkCGvWrAEgMTGRYcOGsW7dOho0aABA8eLFCQoKYsqUKa+cFKnVasqUKcOVK1eee93U1BRTU9mMTwiRQ13aCms/geQYsLSHFovBs7KhozK4Fhtb8Dj1MQAF7yn02ZxOs4e6a4eLqDjathizmy/nTOhjBv1xlluRuoSobVk3htUvjLXZh12Qnt1luaSoYsWKz8wQu3z5Mu7u7gCkpqaSmpqKWp2xHEqj0aDVal/5fRRFISgoCD8/v7cPWgghsgttOuwZDwem6I7zl9WtP2TjbNi4sojHqY8xTlVodUBLo+MKagWiLWBPa2+a9ZxMOytvJmy5xOLDN1AUyGdrxqTmxaniY2/o0EUmyHJJ0YABA6hQoQLjx4+nVatWHD9+nHnz5jFv3jxAN62+atWqDBo0CHNzc9zd3dm3bx9Lly5l2rRp+n46deqEi4sLEyZMAGDMmDGUK1cOb29vYmNjmTVrFkFBQfz0008GuU8hhHjv4iNhTXe4vkd3XPYTqPMdGJkYNq4spFNaGUotOoLzI93xgaIqKn2/mKGeZTkdGkX9RQe48TAegFYB+fmmoS82MjqUY2S5pKhMmTKsW7eOoUOHMnbsWDw9PZkxYwbt27fXt1m5ciVDhw6lffv2PHr0CHd3d8aNG0fv3r31bUJDQzOMJkVHR9OrVy/u37+Pra0t/v7+7N+/n7Jly77X+xNCCIO4ewpWd4aY22BkDo1nQfFWho4qy/j7dhAXvhtCw32hADyygvl11RhVLk8LJ1/GbwlhwUHd6JCjjSkTmxeneiEHA0ctMluWK7TOql6nUEsIIbIMRYETC2D7MEhPgdwFdLvbOxY1dGRZQnBkMJtXfMdHS89g/2Tlld0lVCyroSbBTMX0cusYv+kO1yN0o0MtSudnRENfbM1ldCi7yNaF1kIIITJJchxs+gLO/647LtxQN93e3M6QUWUJwZHBLDw0E69lB2h4Xjc2EG4Lt/o0oGObEfz8WzVSwmvRa8kltAo4WJsyoZkfNYs4Gjhy8S69dVKUmprK/fv3SUhIwN7enty5c2dGXEIIId5GxCXddPuIv2W6/b/EJMfw4/T2dNyahF08aIGtZVSsrKImme1En66MU9R0rkbqRoea+rswqpEvdhZSe5XTvVFSFBcXx/Lly1mxYgXHjx/PsJ5P/vz5qVOnDr169aJMmTKZFqgQQohXdP4P2NgPUuN1u9u3XAzuFQwdVZaQGh5O7Lff8vlO3Wbgd/LAnPoaruRXoWiNSAmvxYy/jdAq8eS1MmV802LUKSrbQ30oXjspmj59OuPGjcPDw4PGjRvz9ddf4+Ligrm5OY8ePeLChQscOHCA2rVrU65cOX744YfnrhgthBAik6Ulw45v4Lhuti4elaHFIrD6cAuCFUXhSvQVDt45QMz69VRdex3zRC1atYrNFU1ZUS6VfLncaG7Xmj2nXLgVqftHfuMSzoxpXJRcljI69CF57ULrli1bMnLkyP9c3yc5OZmFCxdiYmJCjx493irIrEAKrYUQWVp0KPzeRTfLDKDyl1B9OKg/vH23YpJj2HJjC4fvHebIvSPYRCbRc5uWkjd0H3fXnODn+hpuOaoomtufwnzJL4fvoCiQ18qUcU2LESijQzlGtt7mI6uSpEgIkWVd2QVre0BiFJjZQbN54PPyDbJzqscpj6m2qhop2hRUikLgKYV2e7WYpUKKERyt706+Hp8QkxbH3Qgzth7Pw63IRACa+bswUmqHchyZfSaEEB8CbTrsmwT7JgMKOPtDy18gl7uhIzMYjUqDq7UrCdev8umWdArf0Z0PdoW59TS0qNWMWh6N+H77JZYcvomiJOJoY8r4pjKzTLxmUhQVFYWiKOTOnZuIiAj2799PoUKFKFas2LuKTwghxPPEP4Q1Pf6/OnVAN6g7EYw+zD0bk9KSGLx/MPtv/sWgS974bwFVKiSawPJqas5VdOKrjwZjk16aujMOEPooAdCtSj28gaw7JHTU/91EZ8GCBQQEBFC6dGnmzJlD06ZN2b17N23atNFvwSGEEOI9CD0GP1fWJUTGFtB0HjSc/sEmRADxqfHcPvoXkxanU2rD36hS0zlTQMXgXqa4denFbw3WceisC23nHyP0UQL5bM1Y0rUMk1uUkIRI6L1yTVGJEiU4duwYCQkJuLm5cePGDezt7YmNjaVKlSoEBQW941ANS2qKhBAGpyhwdA7sHAHaNMjjDa2XgUMRQ0dmENdjrrP04lJioh/gvz6E0vsfoAZizWFxbTWHfFX80XgNDx/lYciac9yJ0tUOyY72H5Z3UlOk0WgwMzPDzMwMLy8v7O11OwLb2NigksXAhBDi3UqKgQ19IWSj7rhoM93+ZabWho3LQKKToum5oyfO5+/Tc5tWv0XHvmIqltZU89hChY9NcZbsTWbF8WMAuNiZM6l5cSp55zVg5CIre+WkyMjIiKSkJMzMzNi3b5/+/OPHj99JYEIIIZ64FwS/d4aom6A2hsBxULbXB7k6dbo2navRV+m6qjlddmmpfFH3sCPNMTcpX3ajeqVKfGxqy/nQNEatv8SpmNsAdCznzpB6hbEylflF4sVe+W/HX3/9hamp7nm1ra2t/nxiYiILFy7M/MiEEOJD9+/NXO3coOUScClt6Mjeq+ikaBZeWMiFhxe4FPk3Jc/EMn23FptE0Kogb+cu2Pf7HLWFBTEJqXy7OZg/Tummnbnm1o0OVSgoo0Piv71yUmRlZfXc8w4ODjg4fLirpQohxDuRFAt/9oeLa3XHhRrAxz+BeS7DxmUA35/8no3XNmIfrfD5PxZhfJTfhsKTZ5Gn1EcAbLtwnxEbLhDxOBmVCjqX92Bw3UJYmMjokHg18jdFCCGymvvnYXVneHQN1EZQawyU/+yDfFwG0KJgM9JXrKf1ft0ijGlGKvJ9/gWFu3VFZWxMxONkRm+8yObzYQAUsLdkcvPiBHjIBuXi9WR6UhQZGcnZs2cJCgpi4MCBmd29EELkXIoCp3+BrUMgLQls8us2c3Uta+jI3rmY5Bi239yOoiiYGplipjHDVGOKxa1wrKYupfPfWgAuusG8ump2fNILRVFYd+YOY/4MJjohFY1aRe+qBfi8hjdmxh/e9ibi7b1yUnT16lVGjBiBnZ0d48ePJ1euXFy5coWgoCB9EnT27Fnu3buHoihYWlpKUiSEEK8qOQ42DYDzq3XH3nWg6Vyw+DBGO5YGL2Xeuf+veWecqtDikJZGxxTUWog3hWU11OwpoaKWRx3uRScyfN159lyKAMA3nw2TWxSnmIvti95CiP/0yklR+/bt6dChA56enhQtWpTHjx8THx+Pra0tvr6+FCtWjK1bt7Jw4UJq1qyJq6vru4xbCCFyjgfButllDy+DSgM1R0KFfqB+5fV1s73qrtVZdGERado0yt4xo/vWVHI91O1Yf66oJasaWBNhkYaNAmbxtagzfT9xyWmYaNT0r+VNryoFMNZ8ON8v8W68clL08OFDihUrRoECBQgPD2fIkCH06dMHFxcXfZtFixZRtmxZSYiEEOJVnVkOm7+EtESwzgctFoN7eUNH9d4VsC2AS6o1jTY/pMrFOACMHB1xGvENRWrVojVwKzKeIWvO8dvZR0AapdzsmNyiOF4OH+ZaTSLzvXJSNHPmTHr37o29vT0///wzM2fO5OLFi0yePBkfH593GaMQQuQ8KQmw5SsIWq47LlgDms0Hyw9r6nhiWiLdtnYlz97zjNytxToJtMDtwGLUHrcEjZUl6VqFxYduMGXHJZJStZgbaxgUWIjOFTzQqD/M4nPxbrxyUtSwYUMaNmyoP+7atStz5syhSpUqNG/enFGjRr2TAIUQIseJuKSbXRYRAio1VBsGlb/8oB6XJacns/7KehZu+Zae27UUu6WbZp/k6UTsF+2pUL0VGhNLLj94zOA/zhF0OxqAil55mNC0OG55LAwYvcip3nj2mUajoW/fvrRv355Ro0ZRuHBhtFot6enpmRmfEELkLEErdI/LUuPByhGaLwDPKoaO6r3admMbU49MosLecL4/pMUkHVKMVeT7vD/2XbuhMjYmJU3LrN1X+OGvK6SmK1ibGjG8QRFal3GVraXEO/PKG8L+l+DgYAYMGMCZM2cYPHgwn332Gebm5pnRdZYgG8IKId5KchxsGQRnf9Mde1aB5gvB6sNb/Lbn95VouSES14e64yBPFUvqGjOv60bcbdw5ExrF12vOc+mBbhupWkUc+O5jP5xszQwYtciuXufzO9OSoqc2bdrEV199RUxMDGFhYZnZtUFJUiSEeGP3L8AfXZ/MLlNDtaFPHpd9OGvppGnTuBx6hofTpmO/4wwAMRawpJZuN3tUKhbX/o3Np9QsOXwTRYE8liaMbORL4xLOMjok3tjrfH6/9uOz0NBQ3NzcXni9YcOGBAYG8uOPPwJw9+7dDDPUhBDig6EocGoJbPtatxijdT7d4zKPSoaO7L3YcXMHSy4u4XzEOcr/rdBlpxb7eN213SVU/FpdTflCtRnsWArTZH/6Lw3nbnQiAM38XRjR0JdcliYGvAPxoXntpKhMmTI0btyYnj17Urbs81dZTUhIwNLSkmLFivHJJ5/w+eefv3WgQgiRrfx77zKv2tD05w9mdtnBuwf5ct+X5I1RGLJdS+lruocSd3PDvHoaQtxUWBhZMKLMRL7bHMLaM7cAcLEzZ1zTYlQr9OE9VhSG99pJUUhICOPHj6du3boYGxsTEBCAs7MzZmZmREVFERwczMWLFwkICOD777+nXr167yJuIYTIuu6dgd+7QtQN3d5lNUdC+c8/qNlleYxsaXRUS8uDuv3KUjWwvryKdeXVpBmpUBRo6TKJ2tP3ExmfgkoFXSp48FWdQliayracwjDeuKYoKSmJLVu2cODAAW7evEliYiJ58+bF39+fwMBAihUrltmxGpTUFAkh/pOiwLGfYccI0KaCrRu0WASuZQwd2XsRkRDBX6F/EbxvHRV+PYdbhO7jJdgV5tfV4O1fnTrudShsU46Jm2+y++9wAHwcrZjYvDil3HIZMnyRQxm00DqnkqRICPFSCY9gw2dwaYvuuHBDaPIjmOfsD/qY5Bi+3Pclx8KOYZmo0G6vltpBuo+VWHP4tYaavX4qjrQ7ioWRJcuP3WLStkvEJadhrFHRt7o3n1YriInRhzOKJt6vd1poLYQQ4l9Cj8Ef3SD2DmhMoM44KNsTPoAZU0MPDOXYvaNUuqjQabcWuwTd+T3FdYXUjy1UVHapzP1o+HrNEU7eigKglJsdk5oXx9tRtugQWYckRUII8aa0Wjg8E3Z/C0o65C6g27vMuaShI3tvGpsEUHXlXorf1I0O3cmje1QW4qbCxsSGHt5t0UbXpP7MA6Ska7E00TC4bmE6lnNHLVt0iCxGkiIhhHgTcRGw7hO4tlt3XKw5NJwBZjn78XpyejKPEh+hTUkmeckK3JesgBSFFCNYU0HNxnIq0jUqfqr5E1ZaP4auOc+lB9cBqFbInnFN/XCxyzkL+4qcRZIiIYR4Xdf3wtpPIO4+GJlDvUlQqlOOf1yWmJZIvTX1cPo7gp7btTg/0p0P8lSxMFDNg1wqXKxc6Fv8K3afys3So4dRFMhlYcyoRkVpUlIWYRRZmyRFQgjxqtJTYc94ODgdUCBvIWi5BBx9DR1ZplIUhTVX1nA87DjV3arjYOGAvbk9aZGRtF0dTpWLukdl0VYqltcx5UbpfLT2bkpT76acuZHOyD8uEBajW3eoWSkXvmngS25ZhFFkA289++zAgQPMnTuXa9eu8ccff+Di4sKyZcvw9PSkUqWcs2qrzD4T4gMXdRPW9IA7J3THpTpD3YlgkvN2az8fcZ52W9rpj1WKQo2zCu33aLFKAi1g2qIRnoO/QfPk92F4bBKj/7zIlvP3AXDLbcH4pn5U8v4wFqsUWdfrfH6/1RzINWvWEBgYiLm5OWfOnCE5ORmAx48fM378+LfpWgghso4La+HnyrqEyNRWNzrUeFaOTIgAiuQpQmWXygC4hSuMWZbOJ1t1CdENRxjdxRT7EcPR2Nig1Sr8diyUmtP2seX8fTRqFb2rFmT7F1UkIRLZzluNFPn7+zNgwAA6deqEtbU1Z8+epUCBAgQFBVG3bl3u37+fmbEalIwUCfEBSonX7Vt2eqnu2PUj3d5ldi/e/zGnSI+L5+7MqTxevgK1FpKMYVUVNVsDVMypM48KLhW4Gh7HsLXnOX5TV1xUPL8tE5r5UdTZ1sDRC/F/722dokuXLlGlSpVnztvY2BAdHf02XQshhGHdP69be+jhZUCl29W+2lDQ5OxSTEVR+Pv3hSROm415dCJq4EhhFb/UVPPIRlckrVaZMHPXFX7ac5WUdC0WJhq+rFOILhU80Mg0e5GNvdVPd758+bh69SoeHh4Zzh88eJACBQq8TddCCGEYigLH58OObyA9WbezfbN54PnsPwBzmpSbN7n/7Xdw6BDmwH07WFhHjVfdluwu9w0alYbjNx/yzcpgroQ/AHTT7L/7uBj5c+XMR4niw/JWNUWffPIJ/fv359ixY6hUKu7du8fy5cv56quv6NOnzxv3e/fuXTp06ECePHmwsLCgZMmSnDp1Sn89Li6Ovn37kj9/fszNzSlSpAhz5sz5z37XrFmDr68vpqam+Pr6sm7dujeOUQiRAyU8gpXtYOsgXULkUxd6H8rRCdHDxIf8eHQaW4d35krDhsQfOkSKBlZXUvNlTw3VWg5gVPlRJKQojNhwgdZzj3MlPI68VibMauvP4i5lJCESOcZbjRQNHjyYmJgYqlevTlJSElWqVMHU1JSvvvqKvn37vlGfUVFRVKxYkerVq7N161YcHBy4du0adnZ2+jYDBgxgz549/Prrr3h4eLBjxw769OmDs7MzTZo0eW6/R44coXXr1nz77bc0bdqUdevW0apVKw4ePMhHH330RrEKIXKQGwdgbS94fO/JVh3fQdleOX7toeWLv6L40mM4ReuOgzxVLKqjxtW3LKfqLARg+8UHjNp4gQexusk0rQLyM6x+EewsZJq9yFkyZUPYhIQEgoOD0Wq1+Pr6YmVl9cZ9ff311xw6dIgDBw68sE2xYsVo3bo1I0aM0J8rXbo09evX59tvv33ua1q3bk1sbCxbt27Vn6tbty65cuVixYoVz7RPTk7Wz6YDXaGWq6urFFoLkdOkp8G+SbD/e0CBPN66ne3zFTd0ZO9USlgYN78dSfpfBwGItIYltdQcK6QClYo+JfrQxKMLIzdcZFeI7lGZZ15LxjUtRoWCMqtMZB/vbUr+UxYWFgQEBFC2bNm3SogANm7cSEBAAC1btsTBwQF/f3/mz5+foU2lSpXYuHEjd+/eRVEU9uzZw+XLlwkMDHxhv0eOHKFOnToZzgUGBnL48OHntp8wYQK2trb6L1dX17e6LyFEFhR9G35pCPsnAwr4d4BP9uXohGhdyB8M7VOUi3VqkP7XQdJV8GdZFQN6ajhWWA0qFfmt3EiKrECtafvYFfLgyW72XmztX1kSIpGjvfU0iqSkJM6dO0d4eDharTbDtcaNG792f9evX2fOnDkMHDiQYcOGcfz4cfr164epqSmdOnUCYNasWfTs2ZP8+fNjZGSEWq1mwYIFL10s8v79+zg6OmY45+jo+MJlA4YOHcrAgQP1x09HioQQOcTFdfBnf0iKAVMbaDgd/FoYOqp3KvLoQSwGjaBThO74bxdYUFfDR5VaMSSPL755fEmJd2TEhhBmnrgDQIB7LsY388NHdrMXH4C3Soq2bdtGp06dePjw4TPXVCoV6enpr92nVqslICBAv/ijv78/Fy9eZM6cORmSoqNHj7Jx40bc3d3Zv38/ffr0IV++fNSqVeuFff97zx1FUV64D4+pqSmmpqavHb8QIotLjoOtQyDoV92xS4Bu7aHcnoaN6x06/fcejo3sS7VzWtyAWHNYXl3N/apF+KXOz+Q1z0t8chrTdl5m8aGjaBWwMTNiaP0itA5wld3sxQfjrZKivn370rJlS0aOHPnMKMybypcvH76+GfcRKlKkCGvWrAEgMTGRYcOGsW7dOho0aABA8eLFCQoKYsqUKS9MipycnJ4ZFQoPD8+0uIUQ2cDdU7qtOh5d5/9rD30NGmNDR5bpfrn4C45meXHafhbt3GVUe1IiubuEiuXV1CxouYqieYsCsCv4ASM3XOBeTBIAjUs4803DIjhYmxkqfCEM4q2SovDwcAYOHJipiUXFihW5dOlShnOXL1/G3d0dgNTUVFJTU1GrM5ZDaTSaZx7f/VP58uXZuXMnAwYM0J/bsWMHFSpUyLTYhRBZlDYdDs3QbeaqTQOb/Lq1hzwqGjqyd+J27G02rv+e7jvSMX0ApsB1J1hYR0Pxai3Y+9FQTDWm3I9JYvTGi2y7qPsHo2tuc75tUoxqhRwMewNCGMhbJUUtWrRg7969FCxYMLPiYcCAAVSoUIHx48fTqlUrjh8/zrx585g3bx6gWy27atWqDBo0CHNzc9zd3dm3bx9Lly5l2rRp+n46deqEi4sLEyZMAKB///5UqVKFSZMm0aRJEzZs2MCuXbs4ePBgpsUuhMiCYu7A2k/g1pOf9aJNdfVD5rkMG9c7khYZSdz4iXy3WVe+EGcGK6qqMWnagNllBuJk6US6VmHpkZtM3naJuOQ0NGoVPSsXoH9Nb8xNNAa+AyEM562m5CckJNCyZUvs7e3x8/PD2DjjEHS/fv3eqN9NmzYxdOhQrly5gqenJwMHDqRnz5766/fv32fo0KHs2LGDR48e4e7uTq9evRgwYIC+RqhatWp4eHiwZMkS/ev++OMPvvnmG65fv07BggUZN24czZo1e6WYZO8zIbKhi+ufFFNHg7El1P8eSrbLkWsPKenpRK1aRcSMmWhjYwHdo7KTHxemf80R+Dv4AxB8L5Zh684TdDsagJKudkxo5keRfPJ7TeRMr/P5/VZJ0YIFC+jduzfm5ubkyZMnQ9GySqXi+vXrb9p1liNJkRDZSHKcbiPXM8t0x86ldMXUeTJvVDsriT99mrtjR5P+9xUArjvqHpVFeuXhr5Z/oVFrSEhJY+auKyw4eIN0rYK1qRGD6xWmXVk32a9M5GjvbUPYb775hrFjx/L1118/U+MjhBAGcff0k2Lqa4AKKg2A6sNyZDF12qNHhE+ZSszatYDuUdnKKmp2l9JQwtGf78oOQaPWsDvkASM3XORudCIADfzyMbKRL442UkgtxD+9VVKUkpJC69atJSESQhieVguHZ8Jf3z0ppnaBpnPBs7KhI8t0Sno60atXEz59hv5R2Z7iKrbXdSCPSwGWlvqC4vbFCYtJpPeyU/pCahc7c8Y2KUrNIjLrVojneaukqHPnzqxatYphw4ZlVjxCCPH6Yu7Cuk/g5pPtgYo0hkYzwSK3YeN6BxLPnuX+mLEkBQcDEOOeh8lVookv5ML2FtsBSEvXsvDgDabtuER8SjpGahXdK3vSv6Y3FiZvvWavEDnWW/10pKenM3nyZLZv307x4sWfKbT+52wwIYR4J4I3wsbPnxRTW0C9SeDfMUcVUyuKwuXrJ4ib9TMW248AkGxuxNrqZqz3i0ZRq2jmXA6AoNvRDF93nov3dCNIpd1zMa5pMQo7SS2kEP/lrZKi8+fP4++vm9Fw4cKFDNdetFK0EEJkiqRYXTF10HLdcb6S0Hwh5PUyaFiZTUlLY+uUfjis2IPlkwUY9/qp+LW6QqxlEu42HvQt2ZdyTtUZueECy47eQlHA1tyYr+sVlhWphXgNb5UU7dmzJ7PiEEKIVxd6FNb2guhbgAoq9ofqw8HIxNCRZar4Y8d58N13eF55MqvMCc51/AhNcV96W+bD3cadcvnKsf3CQ2r/doCIx7qsqZm/C8MaFCGvlWxVJMTrkIfLQojsIz0V9k6Eg9NA0YKtGzT9OcetTJ0aFsaDyZN5vHWb7oStDXPLx/FXCRWbmo/FzcYNgFuR8fT4JYj9l3U7vBbIa8l3HxejgpfsZC/Em3irpGjs2LEvvT5y5Mi36V4IIf7v4RVY2xPundEdl2irqx8yszVsXJlIm5JC2PyfiZ63AHVyKloV7C1twrJK8cSb62b5JqYlkpyWzvz91/nhr6skp2kxMVLzWTUvelcrgKmRrEgtxJt6q6Ro3bp1GY5TU1O5ceMGRkZGFCxYUJIiIcTbUxQ4uRC2fwNpiWBmB41m6LbryEHCdmziwYQJmIY9Qg2E5IdFdTTcctRirDahiJ0X5Z3LE/koN30WH+BaRDwAlbzy8u3HxfDMa2nYGxAiB3irpOjMmTPPnIuNjaVLly40bZqzfmEJIQzg8QPY2Beu7NAdF6gGH88BG2eDhpWZUm7d4sw3/bA5cRlT4JEVbK/viGWDuvTK60uhXIUoYFeAmAQt47eEMPP0cQDyWpkwoqEvjUs4y8QWITLJW23z8SIXLlygYcOG3Lx5M7O7NhjZ5kOI9+zvzbqp9gmRoDGF2mOg7CeQjReLDYsLY/ONzXjaeOKktsP8ty0k//o7pKaSpobNZVSY9ejIZxW/wlitW+JEq1VYcSKUydsuEZOYikoF7cq6MTiwMLYWOW+VbiEy23vb5uNFoqOjiYmJeRddCyFyuuQ42D4UTi/VHTv6QbN54Ohr2LgywbCDwzh5/wTl/1bouFtL3se68+cKqFlUS8W9PCpKxf2tT4gu3I1h+PoLnH2yeWtRZxu++7gY/m65DHQHQuRsb5UUzZo1K8OxoiiEhYWxbNky6tat+1aBCSE+QLdP6Iqpo24AKqjwOdT4BoxyxtTyDsaVaPDbUYqG6o7DbeGXWmpOeKv0i01ejLxIeNxjZv8VytIjN9EqYG1qxJd1fOhQzh0jTfYdKRMiq3urx2eenp4ZjtVqNfb29tSoUYOhQ4dibW391gFmFfL4TIh3KD0V9k+B/d+Dkg42+XVT7XPIvmVpUVFEzJxJ9OrfQaslzVjN2o9gQzkVqcb/rwcaUuZrEqOKMH/vI8KfrDnUuIQz3zQogoNs3irEG3lvj89u3LjxNi8XQgiIuAzre8PdU7pjv5ZQfwqY2xk0rMygpKURtWIlET/8oN+4NaRkbrY1cOCI9mqGtl4WFdh00IOj18MA3ZpDY5sUo5K3rDkkxPsiizcKIQxDq4Xj82DXKEhL0q031GAa+LUwdGSZIv7IER6MH0/yFV3yE+duz9RKUVx0iwVtrL6dojXCKq4ZFy8HkJIejamRmr7VvehVVdYcEuJ9e+2kaODAga/cVjaEFUI8V/RtWP/p/3e1L1AdmvwEti6GjSsTpNy5w7Vvv4F9xwBItDTiYH03FhS4haJW0bhgYyq5VOJe3D2OX4/nQFA+7idaAFqqFbJnbONiuOWxMOxNCPGBeu2k6HlrEz2PrJshhHiGosDZFbB1CCTH6na1r/MtBHTP9rvaa+PjeThvPg8XLUKVmkq6CraXVvF7JYV481BMNWYM/2g4Tb2bEhaTyMZDwWy/cB+AfLZmjGrkS2BRJ/ndKYQBvXZSJJvACiHeSFwE/NkfLm3WHecvqyumzlPQsHG9JUVRiN20ifDvp5AWHo4KOOehYkktNd91/IV8jy7xIOEB9T3rU8DWm3n7rzFj1xUSUtLRqFV0q+jBF7V8sDSVagYhDO2tfwqjo6NZuHAhISEhqFQqfH196datG7a2OWc/IiHEWwr5E/78AhIegtoYqg/T7Wyvzt41M4kXLnBrzEiU8yEAPLCDpTV1U+ydrVwo7Via0o6lATh6PZL6Sw5wJTwOgAD3XHzXtBiFnWQ2qxBZxVslRSdPniQwMBBzc3PKli2LoihMmzaNcePGsWPHDkqVKpVZcQohsqPEaNj2te6RGYBjMd3okJOfQcN6W4/DbvPn120pcSwSNZBkDOsqqNlUVkWqke7xl6etbsmS8MdJjN8cwvqgewDksTTh63qFaV4qP2q1PCoTIit5q3WKKleujJeXF/Pnz8fISJdfpaWl0aNHD65fv87+/fszLVBDk3WKhHhN1/bAhs8g9i6o1LqRoWpDs/VCjNrkZB4tXcr92T+iSUwBYH9RFcurq4myzpjgjK0wjqgHfkzfeZnHyWmoVND+IzcG1ZHtOYR4n17n8/utkiJzc3POnDlD4cKFM5wPDg4mICCAhISEN+06y5GkSIhXlJIAu0bD8bm641ye0HQuuH1k0LDehqIoPN65k/DJ35N65w4AN12MWVBDy+X8/0+G3Kzd6FasG06aCny36SohYbqp9yVc7fiuSTH88ktZgRDv23tbvNHGxobQ0NBnkqLbt2/nqNWshRCv6M5JWPcJRD5ZmDCgO9QeC6ZWho3rLSQFB/NgwkQSTpwAIMHOjFPNfLlSJh+XQ7dnaKtR7DgcVIA1p08DYGtuzJC6hWlTxlUelQmRDbxVUtS6dWu6d+/OlClTqFChAiqVioMHDzJo0CDatm2bWTEKIbK6tGTYNwkOTgdFC9b5oMmP4FXL0JG9sbSHD3Vbc/yxBhQFlakpuytZsbhkNMkm5yD0nL6toqhIjS7LuUuBnNPqaodaB7gypF5hcluaGOoWhBCv6a2SoilTpqBSqejUqRNpaWkAGBsb8+mnnzJx4sRMCVAIkcXdOwPr+0B4sO7YryXU/x7Ms+dO7k/rhiJ/nos2Ph6Ak37mHG5SkEsmkSQnqrA0tuTTEp+iQsXdSCM2HbPkbqTu16lvPmu+/diP0u7Z8/6F+JC9UU1RUFAQJUuW1B8nJCRw7do1FEXBy8sLC4uctxqr1BQJ8S9pKbB/MhyYptvE1SIvNJwOvo0NHdkbeV7dUJxXPiaWD89QN/TU6vqbWHYwlhXHQ1EUsDYz4qs6hehQzh2NPCoTIst454XWarUaf39/evToQbt27T6INYkkKRLiH+4FPRkduqg7LtpUt4mrZfbcvDQpJIQH4yfo64ZibYx42LkuI6x2kko6fUv2paBdQe7G3eVeXBihd104dtGRR/GpADTzd2Fo/SLYW2ffmXVC5FTvPCk6cuQIixYtYvXq1aSmptKsWTO6d+9O9erV3zjorE6SIiHQjQ4dmAIHpoI2DSzyQIOpuqQoG/p33VCasZr1ZRU2lFOTbKIb7anjXocpVXWlAufuRDNiw0XO3o4GwMfRim+bFOOjAnkMeBdCiJd5b1PyExMTWb16NYsXL+bAgQN4eHjQrVs3OnfuTP78+d+02yxJkiLxwbt/HtZ9Cg/O6459m0D9qWBlb9i43oA2KYlHvywlcu5ctE+WDgn2z83cigmEWafRuGBj/gr9CydLJ5bWW0paqimTt19i5QndozIrUyO+qOVN5woeGGvUBr4bIcTLvLek6J+uXbvG4sWLWbp0KWFhYdSuXZstW7ZkRtdZgiRF4oOVnqobGdr/vW50yDy3bnSoWDNDR/baFK2W2M1bCJ8+jbR7YQBofAsx7qN7nHZKBHQrUa9vsh6toiVNq2XNqTC+336J6ATdo7Km/i4MrVcYBxszg92HEOLVGSQpAoiLi2P58uUMGzaM6Oho0tPTM6trg5OkSHyQ7l+A9Z/C/SfTzws31BVTWzkYNq43kHD6NA8mTiLpnO5e4nNbsLBiEoeKqlCe7Ew//KPh1C9QHxsTG4JuRzNywwXO3YkBoLCTNWObFKOsZ26D3YMQ4vW9t8Ubn9q3bx+LFi1izZo1aDQaWrVqRffu3TOjayGEIaSnwsEZurWHtKm66fX1p0Cx5qDKXjOrUkJDCZ86jcfbdQstqi0seNiyKp/b7yDVWI2R2gg7Uzu87bxp4tWExGQNX/95jlUnb+tmlZkaMaC2D53Ku2Mkj8qEyNHeOCm6ffs2S5YsYcmSJdy4cYMKFSrwww8/0KpVKywtLTMzRiHE+/QgGNb3hrCzuuNCDXSjQ9aOho3rNaXHxvJwzs88+vVXSE0FtRq7Fi2w6tODtjubkKpV4ZvHl5UNVqJSqUjXKvx2PJQp2y8Rk/hkVlkpF76uVxgHa3lUJsSH4I2Sotq1a7Nnzx7s7e3p1KkT3bp1o1ChQpkdmxDifXo6OrR/MqSngJmdbhFGv5bZanRISU0latVqHv74I+nR0QA8KJaPUy39eOCUSNS5saRqdUnP1KpTUalUnA6NYuSGC1y4q9urrEg+G75tUpQAD3lUJsSH5I2SInNzc9asWUPDhg3RaDSZHZMQ4n27dwY29IUHF3THPvWg0QywdjJoWK9DURTi9uwh/PsppNy4AcDdvGqW1ICzBSMg6S+4+f/21fJXw0xlz6Dfz/L7Kd1ijU8XYGz/kZs8KhPiA/RGSdHGjRszOw4hhCGkJsLeiXD4B92q1Oa5od5k8GuRrUaHkoKDeTBpMgnHjgGgymXHyUZefJ/vDHktHfnEuylmRmaYacwwNTLFVG3B3bue1Jiyl9gk3RZFLUvnZ0i9wuS1kgUYhfhQZUqhtRAiG7p1GDZ+/v8d7Ys1h7qTstW6Q6n37xMxcxYx69eDooCJMVfqFOb7IteJNgoCVHQq2onORTvrX3P0eiSjN17k7/u60aSizjaMbVJM9ioTQpAlx4fv3r1Lhw4dyJMnDxYWFpQsWZJTp07pr6tUqud+ff/99y/sc8mSJc99TVJS0vu4JSGyjuTHsPlLWFxPlxBZ54M2K6DFomyTEKXHxRE+fQbX6tYjZt06UBRulM3PZ921DPcLIdooGb+8fsyqPotOvp0ACItJ5PMVZ2gz7yh/33+MnYUx45oWY2PfSpIQCSGALDhSFBUVRcWKFalevTpbt27FwcGBa9euYWdnp28TFhaW4TVbt26le/fuNG/e/KV929jYcOnSpQznzMxkVon4gFzdBX9+ATG3dcelOkHtb8HczpBRvTJ9EfVPP5EeFQXAVTdjFlXTctXlPqDCTGPGzBozKZ+vPCqViuS0dBYcuMGPf10lMTUdlQraf+TGl7ULkcvSxLA3JITIUrJcUjRp0iRcXV1ZvHix/pyHh0eGNk5OGYs/N2zYQPXq1SlQoMBL+1apVM+89kWSk5NJTk7WH8fGxr7S64TIkhIewfbhcPY33bGdOzSeBQWqGTSsV/V0B/uIqdNIuXULABNPT2Z8FMlej3h9/VM//360L9IeC2MLAPb8Hc6YPy9yM1K3lUdp91yMaVyUYi45fxNrIcTry3KPzzZu3EhAQAAtW7bEwcEBf39/5s+f/8L2Dx48YPPmza+0WGRcXBzu7u7kz5+fhg0bcubMmRe2nTBhAra2tvovV1fXN7ofIQwueAP89NGThEgF5fpAnyPZJiFKOHOGW+07cLdff1Ju3UKTJw9Oo0Zy7YfP2euZoE+IPGw86Fm8JxbGFtyKjKf7khN0XXKCm5EJ2FubMq1VCf7oXV4SIiHEC2XqNh+Z4enjrIEDB9KyZUuOHz/OF198wdy5c+nUqdMz7SdPnszEiRO5d+/eSx+FHT16lKtXr+Ln50dsbCwzZ85ky5YtnD17Fm9v72faP2+kyNXVVbb5ENnH4wew5SsIeTJbNG8haPITuJYxbFyvKOXWLcKnTdevRK0yM8O6cweuN/DjTNzf/H75d6KTo8lvlR9zY3N6F+9NJecazN5zjXn7r5OSrsVIraJbJU8+r+GFtZmxge9ICGEIBtv7LDOYmJgQEBDA4cOH9ef69evHiRMnOHLkyDPtCxcuTO3atfnhhx9e6320Wi2lSpWiSpUqzJo16z/by95nIttQFDi7ArYNhaRoUBtBpQFQZRAYZf3p5mlRUTycPYeolSv1K1HbNm1KUJPCTLw+j6jkKH3bQrkKsaLBCozURmw5f59xm4O5F6ObPFHZOy+jGvni5WBtqFsRQmQB733vs8yUL18+fH19M5wrUqQIa9aseabtgQMHuHTpEqtWrXrt91Gr1ZQpU4YrV668caxCZDmR12DTALixT3ecr4RudMjJz7BxvQJtUhKPli4jct48tHFxAFhWqYzDl1+xPOUAM05PAsDZ0pkyTmUobl+c+p71ufEwiVEbLnLkeiQALnbmjGjoS2BRR1TZaK0lIYThZbmkqGLFis/MELt8+TLu7u7PtF24cCGlS5emRIkSr/0+iqIQFBSEn1/W/7AQ4j+lpcDhWbD/e0hLAiMzqPY1lP8cNFnuxzwDJT2dmD//JGLmLNKezCw1LVIEx8GDIKA4a66uZ8bpGfr2m5puwlhjTExCKlO2XWbZ0VukaxVMjdT0rlqQ3lULYm4iK+0LIV5flvttOWDAACpUqMD48eNp1aoVx48fZ968ecybNy9Du9jYWH7//XemTp363H46deqEi4sLEyZMAGDMmDGUK1cOb29vYmNjmTVrFkFBQfz000/v/J6EeKdCj8GmLyA8WHdcoDo0nAa5Xz4b09AURSFu3z4ipk4j+cmIrVG+fDh80R/z+nVZeXkVC9YMyfC4rHWh1qhVRvx2LJQpOy7xKD4FgDq+joxo6ItrbguD3IsQImfIcklRmTJlWLduHUOHDmXs2LF4enoyY8YM2rdvn6HdypUrURSFtm3bPref0NBQ1Or/T66Ljo6mV69e3L9/H1tbW/z9/dm/fz9ly5Z9p/cjxDuTGA27x8DJxYACFnkgcAIUb5Xlt+hIDAoifMpUEk6eBEBtY0PeXj3J1aEDceoUuu/sQVBEEADuNu60LdyWuh51ufFAReMfD3Lxnm6JDC8HK0Y18qWyd/ZYdFIIkbVluULrrEoKrUWWoSi6afZbh0Dcfd25kh2gzrdgkbV3dU++foOI6dN5vHMnACoTE3J36ohF147cVUVzM/Ymk45PIiIxAhO1CUM/GsrHXh/z8HEaE7eGsD7oHqDbuHVALR86lnfHWDZuFUK8RLYutBZCvET0bd00+8vbdMd5vKDhDPCsbNCw/kvqg3Ae/vQT0WvWQHr6kxllH3OqgRc/3lvJvU2/PPOaT0p8QkPPpszdd4Of9lwlIUW3GnXrAFe+CiwkG7cKITKdJEVCZAfpaXB8Lvw1DlLjQW2sm2Zf+Uswzrpb1aQ/fkzkgoU8+uUXlCf7DFrVqEF6z9b8knKI3/6epm9rY2KDh60HHjYeeNl6Y6/UIHDGfm49WY26lJsdYxoXwy+/LL4ohHg3JCkSIqu7FwR/9oOws7pjt/K60SGHwoaM6qW0yclE/baCyJ9/Jj0mBgBzf38cBn3FTpvbfHPoM33bDkU68EnxT7AzswPgangcYzcFs//yeQAcrE0ZVr8ITUo6yxR7IcQ7JUmREFlVchzsnQBHZ4OiBTNbqD0W/DuBOmvW0ein18+aRdo93fR6k4IFcfhyIFbVq/Mg4QHf/NFZ3352zdlUcqmESqUiNimVH3ZfYfGhm6RpFUw0arpX9uSz6l5YmcqvKiHEuye/aYTIahQF/t4M277+/272RZtB3Ylg7WjY2F5AURTi9u4lYvoMki9fBsDI0ZHcfftwqpQ1v4Uf4tj6adyKvaV/TaBHIJXzV0arVfj95G0mb/+bh3G6Kfa1ijjwTQNfPPJaGuR+hBAfJkmKhMhKHl2HLYPhqm52FrZuujWHvGsbNq6XiD9+nIjpM0h8ssHyP6fXTzv/I78c/H8RtVqlpoR9CXr49aCyS2VO3XrE6I3BnL+re8RWIK8lIxr5Ur2Qg0HuRQjxYZOkSIisIDURDk6HgzMgPVlXSF2xn66Q2iRrjpYkXrhIxIwZxB88COg2bM3doT15evRAY2fHv1f76FqsKz39emJtYk1YTCJfrApiw9Mp9qZGfF7Tiy4VPDExypqPBoUQOZ8kRUIY2qVtsHUwRD95tFSgOtSfAnm9DBvXCyRfv07EzFn63esxMsKuZQvy9v4UY0cHLj26xKL9Ezhw9wCPUx4D4JvHlz4l+oBizKzdV5iz9xqJqf+fYv9lnULYW8sUeyGEYUlSJIShRN3U7WR/aYvu2NoZ6k4A3yZZckXq1Lt3ifhpNjHr14NWCyoVNo0aYv/555i4unI//j6jd/bm0L1D+tcYqY2o7lqdsRXGsjs4ivFbQrgbnQhAGY9cjGpUlGIuMsVeCJE1SFIkxPuWmqTbvPXAVN3mrWojKP8ZVBkMplaGju4ZaZGRPPx5LtErV6KkpgJgVbMm9v37Yebjo2/3a/CvGRKir8t+TSufVlx+kEC3xec5fuMRAM62ZgytX4SGxfPJFHshRJYiSZEQ79PVXbBlkK6gGsCjMjSYCvaFDBvXc6Q/fkzkokU8+mUpSoJuAUWLcuVwGPAF5iVK6NspisLlqMs8Tn2c4fXLL2zgQkhxVp0IRauAmbFuF/tPqsgu9kKIrEmSIiHeh+jbsH0ohPypO7ZygsBxUKx5lntUpk1MJGr5ch7OX4D2ycKLZn5+OAz4AssKFfTtYpJjOBtxlkUXFnHqwSn9eUVRkxpVnpvXGhGcEgpAw+L5GFq/CC525u/3ZoQQ4jVIUiTEu5SWAkd+hP3fQ2oCqDTwUW+o9jWYZa2NhbUpKUT//juRP88lLSICABOvgtj37491rVqoVCriU+OZfmo6x8KOcTP2pv61JmoTKrhUIJ+qJjtO2XDzoa5uqKizDaMaFaWsZ9beqFYIIUCSIiHeDUWBKztg+zCIvKo751YBGkwBx6KGje1flNRUoteu4+HPP5MWpluF2tjFhbyf98W2USNUGg1RSVGkK+kcuHOAVZdW6V/rZu1Geefy1M7XgZ//CmfepQggkTyWJgwKLETLAFc06qw1EiaEEC8iSZEQmS08RJcMXftLd2xpD3W+g+Kts9SjMiUtjZiNf/Jw9mxS79wBdKtQ5+39CXbNm6MyMQFgwrEJ/Pb3bxle65fXj59q/oRaa8XM3VdovzGYNK2CkVpF14oefF7TGxsz4/d+T0II8TYkKRIis8RH6vYqO7kIlHTQmEC5T3ULMJplnWnnSno6sVu28vCnn0i5eRMATd685O3VE7vWrYlVEtkQqlsmwERtwvab2zO83khlRF2PBvx5JoZpO08QnaCbkVariAPD6hehgH3Wm0EnhBCvQpIiId5WeiqcWKBLiJJ0hckUbgh1voXcBQwb2z8oWi2Pd+wg4scfSbl6DQBNrlzk6dGDXO3aojbXFUFPPjCGP6//+czrD7c9jLWJNfsvR/DtpmCuhF8EwMfRihENfansbf/+bkYIId4BSYqEeFP6uqHhEHlFd87RD+qOB88qho3tHxRFIe6vv4j44UeS//4bALWtLXm6diVXhw5orP6/jciZ8DOEPg7N8Hp7c3vaFm5LRIyKLzafYPff4QDksjBmYG0f2pZ1w0gjW3MIIbI/SYqEeBP/rhuyyAs1R4B/R1BnjTV4FEUh/sABImb9QNKFCwCorazI3bkzubt0RmNtnaH9lutbGHpwKFpFm+H8R47VuB9akTor9+vrhjqV96B/TW9sLaRuSAiRc0hSJMTriI+EvePh5OIsWzekKAoJR44QMesHEoOCAFBZWJC7QwfydOuKxs5O3zY2JZZlwcs4cOcAIY9C0CpaKjpXxMXKhcjEaO7f92TLXi9iEm8AUL2QPcMb+OLlIHVDQoicR5IiIV5FWoqubmjfxCxbN6QoCvGHD/Pwp9kknj4N6Hauz9WuHXl6dMco9//XCroSdYWdt3bya8iv+k1bAZp7N2dk+ZEcufaIsX8Gc+nBYyAdLwcrvmlQhGqFHN73bQkhxHsjSZEQL6MocHkb7Pjm/+sNZbG6IUVRiD94kIc//kTi2bMAqExMsGvVijy9emLskDGR2XlrJwP3DtQfF7QtSDe/bhTLWwxVqiO9lp5mV8gDAGzNjRlQy5v25dwxlrohIUQOJ0mREC9y8xDsHgO3j+mOs1jdkKIoxO3bx8PZc0g6dw4Alakpdq1bkad7D4wdMyZDCakJrLu6jnnn5unPTag8gXoe9YhL0jLrryssPbKP1HQFjVpFx3LufFHLGzsLk/d6X0IIYSiSFAnxb2HnYPdYuLpTd2xkptuao/LALFE3pCgKcXv28PCn2SRd1E2LV5mZkatNG/J074aR/bNT42/H3qbfnn5cjb6qP/djjR+p4FyZX4+GMmPXZaKerDdUrZA9w+sXwdvR+pl+hBAiJ5OkSIinIq/BnnFwYY3uWKWBUp2g6mCwcTZsbDxZZ2j3bh7OnkNySAgAKnNzcrVrS56uXTHKm/eZ16RqU9lwdQPTTk3jccpjcpvlxtbUFof/tXff4VGW2cPHvzPpdUghhBQCCYQE6QGpUgQJTQRWQLrogv4sNGURRAFdAd13XdBFBJYFXaVKUZoQECkSCNJ7h9BigPQ2SWbu94+BgSGFRAKTcj7XNVcy93PmmfskYebwzF2cfMhOrUXnmTs4fzMdgFo+rkzqXoe2obLekBCiYpKiSIiUG7DjMzjwLRhzTW11/wLt3wevEOv2jTvFUNQWbn31FfrTpwHQOjvjMXAgnsNethhAfb/jt4/z3o73zBu31veuz7/a/4ukVCf+vv4kw7cdBMDTxZ4xz4XSv2mgrDckhKjQpCgSFVdmIuyaCXvnQq5pV3dqPmcaN1S1gVW7BneKoU2bTFeGzpoWh9S6uOAxeBCeQ4di6+GR5zG3Mm+x9fJW1l5Yy+GbpkHXno6evFr3VToF9mbWpossiYnFqMDORsOwVjV4s31NdE6y3pAQQkhRJCqe7HTY+zXsmgX6O9PrA5tBh8lQvZV1+4Zp1/rk9eu5Pf8/ZJ83bcehdXXFc8hgPIcMsVhn6K70nHTG7xjP9qvbzW1ajZaO1Toyvsn7rDmQyHPLfyNVb7oSFvlUFSZ0Cae6t0uecwkhREUlRZGoOHKz4cA3sOMfkGaaco5PHejwIYR2tvoO9sasLJJWriRhwX/JuX4dAK27O55DhuA5ZDA27u55HpOVm8Whm4eYf2Q+MXExANT1qkvnGp3pUr0LBy4a6DPnMJdvZwDwlJ87k7rVoUWI15NLTAghyggpikT5Z8iBoz+YFl5MvGRqqxRkGjNU70WrT683pKaSuGQpCd98g+H2bQBsvLzwHDoUj/4v5dmO434Td00k6rJplpyzrTPzO82nfuX6HLuWzNvfnWDvxQQAKrs5MC6yNn9pHICN1rrFnxBClFZSFInyKycTDn4Hv30ByXc2OXXxMc0mazwUbK27/k7u7dskfPs/Er//HmNaGgB2fn54/vVVKvXujdbRscDHnkk8w65ru8wFEcBXHb/Cxz6Ud1ccZuWBqygFDrZaRrQJ5vW2Ibg4yD93IYQojLxKivInK9m0JceeOZB+09TmUhlavAlPjwB7646jybl+ndv/XUjSDz+gsrIAsA8JwXvEcNy7dkVjV/Cg5zOJZ5h9cDa/XPnFov39ph+x44gb83b8SmaOAYAeDfwY3yUM/0pOjy8ZIYQoR6QoEuVH2k3Y85WpINKnmNoqVYOWI6HRILCzbnGgP3+e2/P/Q/K6dZBrGvDsWL8+3iOG4/rss2i0BU+H1xv0zDsyj/lH5qNQaNDQJqANET5NSb5Vm89XpxGfapqhFhHkwaRu4TSqlnd2mhBCiIJJUSTKvqRY2P2laZ2hXNOVFyqHQesxpvWGbKw73Tzz6DFuz5tH6pYtpr3UAOcWzfEeMQLn5s3RPGSA96ZLm/go+iNSsk2FXqegTrzZ8E3ibrnz9/UnOXHDNA4p0NOJCV3C6VLX96HnFEIIkZcURaLsunkadv0Ljq64t+iifwQ88w6EdoFCrrw8bnc3aU1YuJD03dHmdteOHfAeMQKn+vWLdJ5kfTKLTy4mJTsFH2cfxkSMIdytHdN+PMmWO6tauzna8vazNRnasjoOttbfk00IIcoqKYpE2XNtP+z8HE6tB0xXXqjR1lQM1Whj1an1xuxsUtauJWHRIvRn7+wzZmODrns3vP76Vxxq1XroOZRS7Lq2i/lH53Mw/qC5ffhTY9h3tDqj9+4g12jatHVQs2qM6hiKp4ts2iqEEI9KiiJRNhgNcG4r7JkNF3691x7W3bRRq3+E1boGkJuYSNLSpSR8vxjDrVuAaSuOSn1exGPwEOwD/B96jmR9MmvOrWH9hfWcTDhpbq/sWBXH9M78fZktafrLAHQM9+G9LuHU9HF9PAkJIUQFJEWRKN2Sr5qm1R/4H6RcNbVpbKB+X2g1GnzCrNq97EuXuP3NNySvXmOeSWbr64vn4EFU6tMn3wUX85OancrLP79s3sXeydaJPrX6EqDtzle/XOdCQiZgoE5VdyZ1C6dlzbybvwohhHg0pbIounbtGuPHj2fjxo1kZmYSGhrKggULiIgwXQ0oaBDpZ599xrhx4wo878qVK/nggw84f/48ISEhfPLJJ/Tq1eux5CAegSEXzkXB/kVwdjMoo6ndyQMaDIBmr4FHkNW6p5Qi88ABbi9cSNrWX8yDpx3qhOM1bBjunTsXOq3+fnqDns2XNjNx10QAHG0cGR0xmgC71nwZdY3Zl03bfPi4OfCuLL4ohBCPVakrihITE2nVqhXt27dn48aN+Pj4cP78eSrdt9/TjRs3LB6zceNGXn31Vf7yl78UeN7o6Gj69evHxx9/TK9evVi9ejV9+/Zl165dNGvW7HGlI4ojKdZ0Rejgd5B6/V579WdMiy2GPw92BS9o+Lip3FxSo6K4vXARWUeOmNtd27bFc9gwnJs9XeRZX2cTz7Lq7CrWXlhL8t3914D3Gv+DbQfdWXfkGACOdlpeaxPCiDbBsviiEEI8Zhql7vw3t5R47733+O2339i5c2eRH9OzZ09SU1PZunVrgTH9+vUjJSWFjRs3mts6d+6Mh4cHS5YsyROv1+vR6/Xm+ykpKQQGBpKcnIx7ET8SEUVgyIEzm0xXhc5twTxw2tkLGg4wFUPeDx+c/DgZ0tJJXvkDCd98a96TTGNvj+6FF/B8eSgOISFFPtfqs6v5cPeHFm2ejp7cTs+kkd177D3lRLbBiEYDLzYO4J1OtfHVWa8QFEKIsi4lJQWdTlek9+9S91/Pn376icjISPr06cP27dvx9/fnjTfeYPjw4fnG//HHH6xfv55vvvmm0PNGR0czZswYi7bIyEhmzpyZb/z06dOZOnXqn8pBFEHiJdO6Qge/h7S4e+012kDEy6YB1LYO1uodAPpz50hcvITkH3/EmJ4OgI2HBx4DBuAxoD+2XkXbVDXbkE1GTga5KteiIHrG/xn61HqJ87H+/HvbeXZm5ABGWtf0ZmLXcOr4SfEthBBPUqkrii5cuMCcOXMYO3YsEydOJCYmhpEjR+Lg4MCQIUPyxH/zzTe4ubnRu3fvQs8bFxdHlSpVLNqqVKlCXFxcvvETJkxg7Nix5vt3rxSJR5By3TRG6MSPcH4b964KeUOjgaarQl5Fv+ryOKicHFK3/kLi4sVkxMSY2+2Dg/EcOhTdCz0K3ZPsQZdTLtN/fX9Ss1Mt2p8P7sEzld5iyvJTXL59CoBaPq5M7BZOu9DKsviiEEJYQakrioxGI02aNGHatGkANGrUiOPHjzNnzpx8i6L//ve/DBw4EMcivFE9+EajlCrwzcfBwQEHB+teqSjzjAbTmkJnNsHZTRB31PJ4cHvTVaHaXa2+OWtOfDxJK1aQtGw5ufHxpkatFtdn2+M5YADOLVoUu1AxGA0cv3U8T0FU07Ejp451ZfHlAwB4uzrwTqdQ+kQEYGtjvQUnhRCioit1RVHVqlWpU6eORVt4eDgrV67ME7tz505Onz7NsmXLHnpeX1/fPFeF4uPj81w9Eo8oMwnOb4Uzm00zyDJu33dQY1pPKDQS6r0InsHW6iVwZxbZ77+TuGQJKZujzPuR2Xh5mdYX6tcPu6pVi3y+rNwsVp1dxZFbR7icfJlzSefIMmSZjzsY/Whg8yFbDyYBSTjaaRnxTDAj2obgKoOohRDC6krdK3GrVq04ffq0RduZM2cICso7BfvuNP0GDRo89LwtWrQgKirKYlzR5s2badmy5aN3uiJTCm6eunM1aDPE7gFluHfcwR1CnoXQzlCzI7hWtl5f7zCmp5O8di2Ji5egP3PG3O7UqBEeAwbgFtkJrX3xr1zNOzKP+UfnW7S52bnRuHJr0uNbs+uELVsNSTKIWgghSqlSVxSNGTOGli1bMm3aNPr27UtMTAzz5s1j3rx5FnEpKSmsWLGCf/7zn/meZ8iQIfj7+zN9+nQARo0aRZs2bfj000954YUX+PHHH9myZQu7du167DmVK0qZxgbFHTXNFjuzCZJjLWO8a0NoJ6gVCdWaW31D1rv058+TuGQpyWvWYExLA0Dj5ISue3c8BvTHMTz8T587NiWWUwmnLNp+eH4NO4/Dv7edJykjB1AyiFoIIUqxUlcUNW3alNWrVzNhwgQ++ugjatSowcyZMxk4cKBF3NKlS1FK0b9//3zPExsbi/a+DUFbtmzJ0qVLmTRpEh988AEhISEsW7ZM1igqiNEIyVdMm67ePGX59YExMtg4QI1nTEVQaCfwqG6VLufHmJ5OyqbNJK9eTca+feZ2+6AgPAb0R9erV5FXnc7Pnht7mH9kPjFx9wZlKwVdfSYyYkEsl29nABBaxZWJXcNpK4OohRCi1Cp16xSVVsVZ56BMMRpM0+MfLH5unYGcjPwfo7ExjQeq3spUCAW3BXuXJ9rtwiijkYyYfSSvXk3K5s2ozEzTAa0W13bt8BgwAJeWLdBo//yg5utp15keM51fr/wKgAYNrfxbUd+1Fz/vd+bwFdOCjJXdHBj7nAyiFkIIaynT6xRVOEqBPhVQd7aLUOZtIyzvF/I1Vw/ZaaBPg+x005Uc8/dp9x2706ZPvfc18RIY9Pn3TWtnWjixcm2oHHbvq2eI1WeL5Sc7NpbkNT+SvGaNeZFFAPvq1dH17InuhR7FGjhdkEvJlxj681ASshKw0djQr3Y/nvXtz4Ltt/j0xB9AMs72NoxoE8zwZ2QlaiGEKCvk1draDDkww8rrH9k63il+7it8KoeBRw2wKd1/Ioa0dFI3/Uzy6jVk/P67uV3r6op7167oevXEqWHDEvnI6mD8QTZe3MiSU/dWQF/QYTlr9mUz4McTGIwKrQZeeroaozvWwsdNBlELIURZUrrf8SqCYr1Za+7EP/DVxgEcXMHe1fQxloOb6Wt+981xd45VCoRKQaC1eUwJljzTx2Mxdz4ei7r38ZhGg0vLluh69cKtY4diLbL4MNmGbN7c+qZ5zSFltKNydj+GzrtImt40lb9juA/vdQmjpo9biT2vEEKIJ0eKImvT2sKkO4sFFlT0yMBclFJknztHysaNJK1ZQ+71e5sC21evjq5XL9PHY76+Jfq8RmUk+no0i08tJjU7FQ8HL1q7v8OmA3ZcTDMAudQP0DGhSzgtQoq27YcQQojSSYoia9NorL7HV2mljEayjhwhdcsWUqO2kH35svmY1s0N965dqdSrJ44NGpT4jK74jHjWnFvDqrOruJZ2DaXAkB6K/sYQFidpAQMBHk6Mi6zN8/X90GqlcBVCiLJOiiJRqqicHNJjYkjdsoW0rb/c23ID0NjZ4dyyBboePXDrULIfj4HpqtDu67v54cwP/HrlVwx3FqF0NNTEPqkf1+LdyATcHW15+9laDGkZhINt2fnYUQghROGkKBJWZ8zIIG3XLlMh9Ot2jCkp5mNaFxdc27bBrWNHXNq0wcbV9bH0IVmfzOtRr3Ps9jFzWx331hhvd+P38xqUAnsbLUNaBPHWszWp5Fz6Zt8JIYR4NFIUCaswJCWRuu1XUrdsIX3XLpT+3rIANp6euHV4FreOHXFu0eJPbblRVLcybzF191R2XdtFrjINmK7l1oAQ9QY//p5MtsEIwPMN/PhbZG0CPZ0fW1+EEEJYlxRF4okw6vVkHTlCxv79pO/Za1pd2nBvjzQ7f3/cOnbE7bmOODVqhMbm8XwspZTiYspFLiRd4FraNeYcnkN6TjoAwe618c7uzb4jXhzISgSgRbAXE7qGUT+g0mPpjxBCiNJDiiLxWBhSUsg4cIDM/fvJ+H0/WceOoXJyLGIcQkPNhZBDWNhj3f7iTOIZNl3axOZLm7mUcsnimFIaclMacuvGyxxOygRyCfN1Y3yXMNrJthxCCFFhSFEkSkTOH3+Q8fvv5iJIf/bsvZW577Cp7I1zRBOcGzfGtV1b7KtVeyJ9W356OR/v+dh8305rR5hnGP6u/uSm1WLP8arcSNByjUyq6hwZ+1wovRsHYCMzyoQQokKRokgUmzE7m5zLl8k4eNBcBOVcu5Ynzj4oCKcmEaZCKKIxdtWqWeWqy/oL683fT39mOu0C2nHppoEZG0+x69wtANwcbXmjXU2GtaqOo53MKBNCiIpIiiKRL0NaGjmxsWTHXiH7SqzF97k34vJcBUKrxTE8HKeIxuYiyNbb2zqdv8+tzFsciD9gvh+bkMYH+86x+qCpiLO30TK4RRBvta+Jh4vMKBNCiIpMiqIKSGVnY8zMxJiRQc6NG2THxpITe4Xs2Ng7BdAVDImJhZ5D6+yMY926ODeJwKlxBE4NG2Lj6vKEMijc6YTTzD40mxO3T/BHxh8AKIMT+lvt+fy0G7lGU0H0QkM/3u0kM8qEEEKYSFFkZcpgIOF//wOF6eqL+QqM6Xt1934hx1V2NiozE2NGJsasLIyZGaiMTFPhk5lpeT8rC3Jzi9Q3G09P7KtVw65aIPaB1bCvFojdna82Xl6lcgDyvrh9vLHlDbIMWQAooy2e+r8Qf70BOdlaAFqGeDGhSzj1AnTW7KoQQohSRooiazMaiZ/xqXWe284Ou8qVsatWDfvAQOyDqpmLHrvAwMe2UGJJMyojR28dZVvsNn44+wNZhiya+7Yi1HYQK/ZmEpuSDUCYrxvvdQmjrcwoE0IIkQ8piqxNq8X9+edN32u47836vo1g798UVoP5e82dTWM1dnZonZ3QODmhdXJG6+yE1qmA+87OaJ3u3Leze5KZlriErASWn17OitMriM80bQeiFPhpOnHh+PNExScD4F/JiXc6hfJCQ3+ZUSaEEKJAUhRZmcbGBv9/fGbtbpQpBqOB6THTWXZ6mbnN1c6VMOeuxF5qzJk4gHQqOdvxVvuaDGoeJDPKhBBCPJQURaJMyTXmsu7COouCqK1PP/TxkUTtMV0tcrTT8kqrGrzWNgSdU9m+GiaEEOLJkaJIlGqJWYkcvnmYIzePcPjmYY7eOkpmbiYAxhx3tInd2XCqPkYVj1YDfZsEMrpjKL46Ryv3XAghRFkjRZEolVadXcXk3ZPzPeas9cIlvSdXr9Yi+85EuufqVGF859rU9HF7gr0UQghRnkhRJEqlVWdX5Wkb3ehd4m+Es2xPEn9kmvZRiwjyYEKXMJpU93zSXRRCCFHOSFEkSg2lFHHpcZxMOEnNSjU5fPPwnXYNDhmt+c96f24k3wSgpo8r4zuH0THcR6bXCyGEKBFSFAmrOnzzMFsvb+VUwilOJpwkSZ9kPqYUqPSnsE9+kdspToCeKu4OjH0ulL80DsDWRmu1fgshhCh/pCgSVpNrzOWVn18h25htbrPV2BJSKQQvTVNOnQ3jcrwt6YC7oy1vtK/Jyy1lw1YhhBCPhxRFwmqupF6xKIiWdl+KUe/LrKgLbDppml7vYKtlWKsa/F/bEHTOMr1eCCHE4yNFkXiilFIcu3WM9RfXs+L0CgDstHb8rdEM/rs1l1UH96AU2Gg19G0SwKgOMr1eCCHEkyFFkXgiLqdc5qfzP7Hx4kaupF4xtzet3J7K+kF8uOQ22YarAHSp68s7nWpT06ds7L0mhBCifJCiSDx2a8+v5cPfPiRXmRYVcrJ1orXvs2hSOrB5r5FUvWlGWYtgL8Z3CaNhYCUr9lYIIURFJUWReCz0Bj07r+5kw8UNRF2OAuBp36fpGdKbm3+EMufXy9xM1QNQp6o747uE0aaWt0yvF0IIYTVSFIkSpZRi9qHZfH/ye9Jy0sztg8IHU9tuAP9v9Vku3z4DQDVPZ97pFMrz9f3Qyu71QgghrEyKIlFi0rLTWHp6KXOPzAWginMVOlfvgg/tWLI7gzk3TIsxervaM7JDLV5qWg17W1lrSAghROkgRZF4ZOcSz7H09FLWnl9LRm4GAIPCB9G+8l/5f5vOEHMpDgA3B1teaxvMsFY1cHGQPz0hhBCli7wziT/tYPxBvjjwBb//8bu5rYauBu18BnDsRDBzTu0BwN5Wy8stq/N/bUPwcLG3VneFEEKIQklRJIotx5jD14e/5j9H/4NRGbHR2NA+sD3tq77ItsMufLn2GkrdNK81NLJDLarqnKzdbSGEEKJQUhSJIlNKsf3qdv75+z+5lHIJgB4hPRhQ63WW70nmnS2x5BiSAOhWrypjO4USUlnWGhJCCFE2SFEkiuzLg18y/+h8ADwcPBjTcAIXYoPpM/s4GdkGAJ6p5c24yNrUD6hkxZ4KIYQQxSdFkXgoozKy+uxqFh5bCMCgsGG4ZHRh6vIrJGWcA6BBYCXGR9amZU1va3ZVCCGE+NNK5Xzoa9euMWjQILy8vHB2dqZhw4bs37/fIubkyZP06NEDnU6Hm5sbzZs3JzY2tsBzLlq0CI1Gk+eWlZX1uNMp0y4kXWDwxsFMiZ5CjtFIiGYoq7c24P9tukBSRg41fVz5elAEa95oKQWREEKIMq3UXSlKTEykVatWtG/fno0bN+Lj48P58+epVKmSOeb8+fO0bt2aV199lalTp6LT6Th58iSOjoVvHOru7s7p06ct2h72mIosWZ/M8Kjh/JF+E216U+ySenIo1QbIwk/nyOiOofRu7I+tTamsrYUQQohiKXVF0aeffkpgYCALFy40t1WvXt0i5v3336dr16589tln5rbg4OCHnluj0eDr61tifS2vkvXJrDizgsUnlnDjpgcqYQhZGaarQF4u9rzZviYDm1fDwdbGyj0VQgghSk6p+y/+Tz/9RJMmTejTpw8+Pj40atSI+fPnm48bjUbWr19PaGgokZGR+Pj40KxZM9asWfPQc6elpREUFERAQADdu3fn4MGDBcbq9XpSUlIsbuVZrjGX3dd38+FvH/LcD8/xz10/cflUbzKvDiMrwxs3B1veeS6UHX9rzyuta0hBJIQQotzRKKWUtTtxv7sfZ40dO5Y+ffoQExPD6NGjmTt3LkOGDCEuLo6qVavi7OzM3//+d9q3b8/PP//MxIkT2bZtG23bts33vHv27OHcuXPUq1ePlJQUZs2axYYNGzh8+DC1atXKEz9lyhSmTp2apz05ORl3d/eSTdrKtl/ZztToqdzMvIkh0w/9zUgM6bUBcLDV8nIr08KLlZxl4UUhhBBlS0pKCjqdrkjv36WuKLK3t6dJkybs3r3b3DZy5Ej27dtHdHQ0169fx9/fn/79+7N48WJzTI8ePXBxcWHJkiVFeh6j0Ujjxo1p06YNX3zxRZ7jer0evV5vvp+SkkJgYGC5KoqMysjcw3P56vBXGPXeGBO6kZkUDoCtVsNLTwfy9rO1qOIu466EEEKUTcUpikrdmKKqVatSp04di7bw8HBWrlwJgLe3N7a2tvnG7Nq1q8jPo9Vqadq0KWfPns33uIODAw4ODsXsfdmw/4/9/HDmBw7fPMzlhBSyb/6F3OQmKDRoNNCzoT+jO9YiyMvF2l0VQgghnphSVxS1atUqzwyxM2fOEBQUBJiuJDVt2rTQmKJQSnHo0CHq1av36J0uI04nnGba3mkciD+AMdeF7FvtyUlqDsr0Z/BcnSq80ymUMN/ycSVMCCGEKI5SVxSNGTOGli1bMm3aNPr27UtMTAzz5s1j3rx55phx48bRr18/2rRpYx5TtHbtWn799VdzzJAhQ/D392f69OkATJ06lebNm1OrVi1SUlL44osvOHToELNnz37SKT5xRmXkfyf+x6wDs8jOsSE3oTO5ic+QazANlm4R7MW4zrVpXM3Dyj0VQgghrKfUFUVNmzZl9erVTJgwgY8++ogaNWowc+ZMBg4caI7p1asXX3/9NdOnT2fkyJHUrl2blStX0rp1a3NMbGwsWu29yXVJSUmMGDGCuLg4dDodjRo1YseOHTz99NNPNL8nLdeYy6TfJrHu3GayE1qhEjuQk2sHQIMAHe9G1qZ1TW80Go2VeyqEEEJYV6kbaF1aFWegVmmQnpPOb9d+Y+nJlfx2SkP2rfYog2lz1tpV3BjbKZROdapIMSSEEKJcK9MDrcWjyTHmsOL0Cr469DU3/6hJ9q0OqNxKAAR5OTP2uVC61/fDRivFkBBCCHE/KYrKietp19lwcQM/nvuJs1fc0d96BZVdGQAvVxve7VSHFyMCsJMtOYQQQoh8SVFUxhmMBhafWswXB74kLak6+ps9MOqrAuB5d0uOZtVwtJMVqIUQQojCSFFUhmXkZDBs0zCOxOrRxw/DmFUNADdHG15rE8KwVjVwcZBfsRBCCFEU8o5ZBhmMBn67/ht/2/w1N6+2xJBREwAnOxuGtarOiDbBsiWHEEIIUUxSFJUhcelxHL99nA9+mccfV57GkG5apsBWqxjUvAZvtA/Bx0225BBCCCH+DCmKyoBbmbeYETODDacPkn3zOXJThwCg0Sh6NfJl7HN1CPBwtnIvhRBCiLJNiqJSKikrCXsbe2y1try2YRKHT/uTmzIa0AKKHg38GPNcbWp4y/5kQgghREmQoqiUUUqx58YeRm8bTVqmHTm3OpKd1BUwzR57JtSN97s0JKxq6V9AUgghhChLpCgqZWYfms2cA9+TfasDOUnNzJu1Nqpux0fdm1EvQGflHgohhBDlkxRFpUhCmp4tB11JP/c3UKbZY9WrZDOxSz06hdW0cu+EEEKI8k2KouJKTwebfBZCtLEBR0fLuIJoteDkZL6bcjuJb367xLfRl0nTu+GEERe3K0zo3JC/NGxgEUtGBhS0XZ1GA87Ofy42MxOMxoL77OLy52KzssBgKJlYZ2dTvwH0esjNLZlYJyfT7wQgOxtyckom1tHx3t9KcWJzckzxBXFwAFvb4sfm5pp+FgWxtwc7u+LHGgym311B7OxM8cWNNRpNf2slEWtra/pZgOnfREZGycQW59/9I7xGFCtWXiNM38trRPFjy/NrRFEpUSTJyckKUMmml5C8t65dLR/g7Jx/HCjVtq057HRcirrtrCs4tkkTy/MGBRUcW6eOZWydOgXHBgVZxjZpUnCst7dlbNu2Bcc6O1vGdu1acOyDf34vvlh4bFravdihQwuPjY+/F/vGG4XHXrx4L/bddwuPPXbsXuzkyYXHxsTci/3ss8Jjt227F/vvfxceu27dvdiFCwuPXb78Xuzy5YXHLlx4L3bdusJj//3ve7HbthUe+9ln92JjYgqPnTz5XuyxY4XHvvvuvdiLFwuPfeONe7Hx8YXHDh16LzYtrfDYF19UFgqL/ZOvEUop07/BgmLlNeLeTV4jTDd5jTDd7rxGmN+/k5PVw8hGWFYWUtkV2ZtVCCGEsD6NUkpZuxNlQUpKCjqdjuTr13F3z2fm1yNcGr9+7RZV3B3z37leLo3fI5fGTeTSePFj5eMzE3mN+HOx8hphUkZfI8zv38nJ+b9/30eKoiIqzg9VCCGEEKVDcd6/5eMzIYQQQgikKBJCCCGEAKQoEkIIIYQApCgSQgghhACkKBJCCCGEAKQoEkIIIYQApCgSQgghhACkKBJCCCGEAKQoEkIIIYQApCgSQgghhACkKBJCCCGEAKQoEkIIIYQApCgSQgghhADA1todKCuUUoBpt10hhBBClA1337fvvo8XRoqiIkpNTQUgMDDQyj0RQgghRHGlpqai0+kKjdGoopROAqPRyPXr13Fzc0Oj0ZTouVNSUggMDOTKlSu4u7uX6LlLm4qUK1SsfCtSrlCx8q1IuULFyrci5KqUIjU1FT8/P7TawkcNyZWiItJqtQQEBDzW53B3dy+3f5QPqki5QsXKtyLlChUr34qUK1SsfMt7rg+7QnSXDLQWQgghhECKIiGEEEIIQIqiUsHBwYHJkyfj4OBg7a48dhUpV6hY+VakXKFi5VuRcoWKlW9FyrUoZKC1EEIIIQRypUgIIYQQApCiSAghhBACkKJICCGEEAKQokgIIYQQApCiyOq++uoratSogaOjIxEREezcudPaXSq26dOn07RpU9zc3PDx8aFnz56cPn3aIkYpxZQpU/Dz88PJyYl27dpx/Phxixi9Xs/bb7+Nt7c3Li4u9OjRg6tXrz7JVIpt+vTpaDQaRo8ebW4rb7leu3aNQYMG4eXlhbOzMw0bNmT//v3m4+Ul39zcXCZNmkSNGjVwcnIiODiYjz76CKPRaI4py7nu2LGD559/Hj8/PzQaDWvWrLE4XlK5JSYmMnjwYHQ6HTqdjsGDB5OUlPSYs7NUWK45OTmMHz+eevXq4eLigp+fH0OGDOH69esW5ygrucLDf7f3e+2119BoNMycOdOivSzl+1gpYTVLly5VdnZ2av78+erEiRNq1KhRysXFRV2+fNnaXSuWyMhItXDhQnXs2DF16NAh1a1bN1WtWjWVlpZmjpkxY4Zyc3NTK1euVEePHlX9+vVTVatWVSkpKeaY119/Xfn7+6uoqCh14MAB1b59e9WgQQOVm5trjbQeKiYmRlWvXl3Vr19fjRo1ytxennJNSEhQQUFB6uWXX1Z79+5VFy9eVFu2bFHnzp0zx5SXfP/+978rLy8vtW7dOnXx4kW1YsUK5erqqmbOnGmOKcu5btiwQb3//vtq5cqVClCrV6+2OF5SuXXu3FnVrVtX7d69W+3evVvVrVtXde/e/UmlqZQqPNekpCTVsWNHtWzZMnXq1CkVHR2tmjVrpiIiIizOUVZyVerhv9u7Vq9erRo0aKD8/PzUv/71L4tjZSnfx0mKIit6+umn1euvv27RFhYWpt577z0r9ahkxMfHK0Bt375dKaWU0WhUvr6+asaMGeaYrKwspdPp1Ndff62UMr1Q2dnZqaVLl5pjrl27prRarfr555+fbAJFkJqaqmrVqqWioqJU27ZtzUVRect1/PjxqnXr1gUeL0/5duvWTb3yyisWbb1791aDBg1SSpWvXB984yyp3E6cOKEAtWfPHnNMdHS0AtSpU6cec1b5K6xIuCsmJkYB5v+QltVclSo436tXryp/f3917NgxFRQUZFEUleV8S5p8fGYl2dnZ7N+/n06dOlm0d+rUid27d1upVyUjOTkZAE9PTwAuXrxIXFycRa4ODg60bdvWnOv+/fvJycmxiPHz86Nu3bql8ufx5ptv0q1bNzp27GjRXt5y/emnn2jSpAl9+vTBx8eHRo0aMX/+fPPx8pRv69at2bp1K2fOnAHg8OHD7Nq1i65duwLlK9cHlVRu0dHR6HQ6mjVrZo5p3rw5Op2uVOefnJyMRqOhUqVKQPnL1Wg0MnjwYMaNG8dTTz2V53h5y/dRyIawVnLr1i0MBgNVqlSxaK9SpQpxcXFW6tWjU0oxduxYWrduTd26dQHM+eSX6+XLl80x9vb2eHh45IkpbT+PpUuXcuDAAfbt25fnWHnL9cKFC8yZM4exY8cyceJEYmJiGDlyJA4ODgwZMqRc5Tt+/HiSk5MJCwvDxsYGg8HAJ598Qv/+/YHy97u9X0nlFhcXh4+PT57z+/j4lNr8s7KyeO+99xgwYIB5Q9Tyluunn36Kra0tI0eOzPd4ecv3UUhRZGUajcbivlIqT1tZ8tZbb3HkyBF27dqV59ifybW0/TyuXLnCqFGj2Lx5M46OjgXGlYdcwfQ/zCZNmjBt2jQAGjVqxPHjx5kzZw5Dhgwxx5WHfJctW8Z3333H4sWLeeqppzh06BCjR4/Gz8+PoUOHmuPKQ64FKYnc8osvrfnn5OTw0ksvYTQa+eqrrx4aXxZz3b9/P7NmzeLAgQPF7ldZzPdRycdnVuLt7Y2NjU2eCjs+Pj7P/9bKirfffpuffvqJbdu2ERAQYG739fUFKDRXX19fsrOzSUxMLDCmNNi/fz/x8fFERERga2uLra0t27dv54svvsDW1tbc1/KQK0DVqlWpU6eORVt4eDixsbFA+frdjhs3jvfee4+XXnqJevXqMXjwYMaMGcP06dOB8pXrg0oqN19fX/744488579582apyz8nJ4e+ffty8eJFoqKizFeJoHzlunPnTuLj46lWrZr5Nevy5cu88847VK9eHShf+T4qKYqsxN7enoiICKKioizao6KiaNmypZV69ecopXjrrbdYtWoVv/zyCzVq1LA4XqNGDXx9fS1yzc7OZvv27eZcIyIisLOzs4i5ceMGx44dK1U/jw4dOnD06FEOHTpkvjVp0oSBAwdy6NAhgoODy02uAK1atcqzvMKZM2cICgoCytfvNiMjA63W8iXRxsbGPCW/POX6oJLKrUWLFiQnJxMTE2OO2bt3L8nJyaUq/7sF0dmzZ9myZQteXl4Wx8tTroMHD+bIkSMWr1l+fn6MGzeOTZs2AeUr30f2pEd2i3vuTslfsGCBOnHihBo9erRycXFRly5dsnbXiuX//u//lE6nU7/++qu6ceOG+ZaRkWGOmTFjhtLpdGrVqlXq6NGjqn///vlO9w0ICFBbtmxRBw4cUM8++2ypmMr8MPfPPlOqfOUaExOjbG1t1SeffKLOnj2rvv/+e+Xs7Ky+++47c0x5yXfo0KHK39/fPCV/1apVytvbW/3tb38zx5TlXFNTU9XBgwfVwYMHFaA+//xzdfDgQfOMq5LKrXPnzqp+/foqOjpaRUdHq3r16j3xaduF5ZqTk6N69OihAgIC1KFDhyxes/R6fZnL9WH55ufB2WdKla18Hycpiqxs9uzZKigoSNnb26vGjRubp7GXJUC+t4ULF5pjjEajmjx5svL19VUODg6qTZs26ujRoxbnyczMVG+99Zby9PRUTk5Oqnv37io2NvYJZ1N8DxZF5S3XtWvXqrp16yoHBwcVFham5s2bZ3G8vOSbkpKiRo0apapVq6YcHR1VcHCwev/99y3eKMtyrtu2bcv33+nQoUOVUiWX2+3bt9XAgQOVm5ubcnNzUwMHDlSJiYlPKEuTwnK9ePFiga9Z27ZtK3O5KvXw3+2D8iuKylK+j5NGKaWexBUpIYQQQojSTMYUCSGEEEIgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBCCCEEIEWREEIIIQQgRZEQQgghBCBFkRBClHoLFiygU6dOxXrMv//9b3r06PGYeiRE+SRFkRCiSDQaTaG3l19+2dpdLHHt2rVj9OjRVu2DXq/nww8/5IMPPjC3TZkyhYYNG1rE7dy5k0qVKvH222+jlGL48OHs27ePXbt2PeEeC1F2SVEkhCiSGzdumG8zZ87E3d3dom3WrFnW7mKR5eTklJnnW7lyJa6urjzzzDMFxqxfv57IyEhGjRrFl19+iUajwcHBgQEDBvDll1/+6ecWoqKRokgIUSS+vr7mm06nQ6PRWLTt2LGDiIgIHB0dCQ4OZurUqeTm5pofr9FomDt3Lt27d8fZ2Znw8HCio6M5d+4c7dq1w8XFhRYtWnD+/HnzY+5eEZk7dy6BgYE4OzvTp08fkpKSLPq2cOFCwsPDcXR0JCwsjK+++sp87NKlS2g0GpYvX067du1wdHTku+++4/bt2/Tv35+AgACcnZ2pV68eS5YsMT/u5ZdfZvv27cyaNct8NezSpUssWrSISpUqWTz/mjVr0Gg0efr93//+l+DgYBwcHFBKkZyczIgRI/Dx8cHd3Z1nn32Ww4cPF/pzX7p0aaEfgy1evJjevXszY8YMpk6danGsR48erFmzhszMzEKfQwhhIkWREOKRbdq0iUGDBjFy5EhOnDjB3LlzWbRoEZ988olF3Mcff8yQIUM4dOgQYWFhDBgwgNdee40JEybw+++/A/DWW29ZPObcuXMsX76ctWvX8vPPP3Po0CHefPNN8/H58+fz/vvv88knn3Dy5EmmTZvGBx98wDfffGNxnvHjxzNy5EhOnjxJZGQkWVlZREREsG7dOo4dO8aIESMYPHgwe/fuBWDWrFm0aNGC4cOHm6+GBQYGFvlncrffK1eu5NChQwB069aNuLg4NmzYwP79+2ncuDEdOnQgISGhwPPs3LmTJk2a5Hts9uzZDBs2jAULFjBy5Mg8x5s0aUJOTg4xMTFF7rcQFZoSQohiWrhwodLpdOb7zzzzjJo2bZpFzP/+9z9VtWpV831ATZo0yXw/OjpaAWrBggXmtiVLlihHR0fz/cmTJysbGxt15coVc9vGjRuVVqtVN27cUEopFRgYqBYvXmzx3B9//LFq0aKFUkqpixcvKkDNnDnzoXl17dpVvfPOO+b7bdu2VaNGjSo0d6WUWr16tbr/5XTy5MnKzs5OxcfHm9u2bt2q3N3dVVZWlsVjQ0JC1Ny5c/PtT2JiogLUjh07LNonT56s7O3t8/z88uPh4aEWLVpUaIwQwsTWmgWZEKJ82L9/P/v27bO4MmQwGMjKyiIjIwNnZ2cA6tevbz5epUoVAOrVq2fRlpWVRUpKCu7u7gBUq1aNgIAAc0yLFi0wGo2cPn0aGxsbrly5wquvvsrw4cPNMbm5ueh0Oos+Pni1xWAwMGPGDJYtW8a1a9fQ6/Xo9XpcXFwe9ccBQFBQEJUrVzbf379/P2lpaXh5eVnEZWZmWnxk+OAxAEdHxzzHAgICqFSpEp999hldunShatWq+Z7DycmJjIyMP5uGEBWKFEVCiEdmNBqZOnUqvXv3znPs/jd0Ozs78/d3x+Dk12Y0Ggt8rrsxGo3GHDd//nyaNWtmEWdjY2Nx/8Fi55///Cf/+te/mDlzJvXq1cPFxYXRo0eTnZ1dcKKAVqtFKWXRlt9A6gefz2g0UrVqVX799dc8sQ+OUbrLy8sLjUZDYmJinmNubm5s2bKFTp060a5dO7Zt24afn1+euISEBIviTAhRMCmKhBCPrHHjxpw+fZqaNWuW+LljY2O5fv26+Q0/OjoarVZLaGgoVapUwd/fnwsXLjBw4MBinXfnzp288MILDBo0CDAVLWfPniU8PNwcY29vj8FgsHhc5cqVSU1NJT093Vz43B0zVJjGjRsTFxeHra0t1atXL1If7e3tqVOnDidOnMh3nSIPDw+2bNlCZGSkuTDy9/c3Hz9//jxZWVk0atSoSM8nREUnA62FEI/sww8/5Ntvv2XKlCkcP36ckydPsmzZMiZNmvTI53Z0dGTo0KEcPnyYnTt3MnLkSPr27Yuvry9gmuk1ffp0Zs2axZkzZzh69CgLFy7k888/L/S8NWvWJCoqit27d3Py5Elee+014uLiLGKqV6/O3r17uXTpErdu3cJoNNKsWTOcnZ2ZOHEi586dY/HixSxatOiheXTs2JEWLVrQs2dPNm3axKVLl9i9ezeTJk0yDzLPT2RkZKFrDel0OjZv3oy3tzft2rXj6tWr5mM7d+4kODiYkJCQh/ZPCCFFkRCiBERGRrJu3TqioqJo2rQpzZs35/PPPycoKOiRz12zZk169+5N165d6dSpE3Xr1rWYcv/Xv/6V//znPyxatIh69erRtm1bFi1aRI0aNQo97wcffEDjxo3NV1l8fX3p2bOnRcy7776LjY0NderUoXLlysTGxuLp6cl3333Hhg0bzNP4p0yZ8tA8NBoNGzZsoE2bNrzyyiuEhoby0ksvcenSJfP4qvwMHz6cDRs2kJycXGCMu7s7mzZtokqVKrRr144rV64AsGTJEouxVkKIwmnUgx+OCyFEKTFlyhTWrFlTpI+nyrO+ffvSqFEjJkyYUOTHHDt2jA4dOnDmzJk8g86FEPmTK0VCCFHK/eMf/8DV1bVYj7l+/TrffvutFERCFIMMtBZCiFIuKCiIt99+u1iPKe4GskII+fhMCCGEEAKQj8+EEEIIIQApioQQQgghACmKhBBCCCEAKYqEEEIIIQApioQQQgghACmKhBBCCCEAKYqEEEIIIQApioQQQgghAPj/YPNvOcqlx8AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(temperature_md_lst, np.array(volume_md_lst)/len(structure_md) * len(structure_opt), label=\"Molecular Dynamics\", color=\"C2\")\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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}