{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "226ac273",
   "metadata": {},
   "source": [
    "# 1.5 Analysizing a Monte Carlo or Molecular Dynamics Trajectory. \n",
    "\n",
    "Here we would learn to analyze a molecular dynamics trajectory. We would compute:\n",
    "\n",
    "(1.5.1) Vibrational density along normal modes\n",
    "\n",
    "(1.5.2) Anharmonic Measure proposed by Carbogno and coworkers.\n",
    "\n",
    "1.5.1. Vibrational Density along Normal Mode\n",
    "====\n",
    "The file co2-b3lyp.pos_0.xyz contains a 10000 step trajectory from a Path-Integral PIGLET simulation. We would discard the first 2000 step (1ps) for equilibration and compute the vibrational density along normal modes. We would use xyz2nmode class of pyepfd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "45e16bf4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                 \n",
      "          ███████████                            \n",
      "         ░░███░░░░░███                           \n",
      "          ░███    ░███ █████ ████                \n",
      "          ░██████████ ░░███ ░███                 \n",
      "          ░███░░░░░░   ░███ ░███                 \n",
      "          ░███         ░███ ░███                 \n",
      "          █████        ░░███████                 \n",
      "         ░░░░░          ░░░░░███                 \n",
      "                        ███ ░███                 \n",
      "                       ░░██████                  \n",
      "                        ░░░░░░                   \n",
      " ██████████ ███████████  ███████████ ██████████  \n",
      "░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███ \n",
      " ░███  █ ░  ░███    ░███ ░███   █ ░  ░███   ░░███\n",
      " ░██████    ░██████████  ░███████    ░███    ░███\n",
      " ░███░░█    ░███░░░░░░   ░███░░░█    ░███    ░███\n",
      " ░███ ░   █ ░███         ░███  ░     ░███    ███ \n",
      " ██████████ █████        █████       ██████████  \n",
      "░░░░░░░░░░ ░░░░░        ░░░░░       ░░░░░░░░░░   \n",
      "PyEPFD version     :  1.1\n",
      "Author             : Arpan Kundu\n",
      "Author Email       : arpan.kundu@gmail.com\n",
      "Today              :  2024-12-18 17:47:54.804572\n",
      "*************************************************\n",
      "                   CITATIONS                     \n",
      "=================================================\n",
      "Please cite the following 3 references:          \n",
      "(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,  \n",
      "L070801, \n",
      "(2) A. Kundu and G Galli, \n",
      "J. Chem. Theory. Comput. (2023), 19, 4011\n",
      "(3) https://pyepfd.readthedocs.io/en/latest/\n",
      "*************************************************\n"
     ]
    }
   ],
   "source": [
    "from pyepfd.xyz2nmode import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e9f4fcd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Process-id0: Time spent on read_pyepfd_info class: 0.0005273818969726562 s.\n",
      "Total execution time (s): 1.6892316341400146\n"
     ]
    }
   ],
   "source": [
    "nm_density = prob_dist_on_nm(\n",
    "            phonon_info_file='../2_normal_mode_phonon/enmfdphonon.xml', #File containing phonon information.\n",
    "            ref_struc='../1_cartesian_phonon/co2_b3lyp_opt.xyz', #xyz file containing optimized geometry\n",
    "            traj_file='co2-b3lyp.pos_0.xyz', #trajectory file\n",
    "            prefix='co2-100K', # The prefix that would be used for outputs\n",
    "            par_exec = False, #If true then run using mpirun\n",
    "            nmode_traj = False, #If true a normal mode trajectory would be written\n",
    "            start_frame = 2001, #First 2000 steps discarded forequilibration\n",
    "            end_frame = 10000, #Default is the last frame\n",
    "            mode_list = [6,7,8,9]\n",
    "            )\n",
    "## deleting the object to finish the process\n",
    "del nm_density"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eb036d7",
   "metadata": {},
   "source": [
    "This created two files: (1) co2-100K.stat.dat & (2) co2-100K.prob-density.dat. The first file prints some statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bbe5a94f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#-----------------------------------------Information-----------------------------------\n",
      "#        Skewness = Fischer - Pearson coefficient of skewness (moment based)         \n",
      "#       Skew-Test Null Hypothesis: The distribution is drawn from a Gaussian         \n",
      "#   Therefore if p-value is less than 0.05 that means the distribution is skewed     \n",
      "#                      with (1 - 0.05)*100 = 95% confidence                          \n",
      "#               Mean, Median and Stdev are in Freq scaled coordinate                 \n",
      "#---------------------------------------------------------------------------------------\n",
      "# Mode No      Mean          Median          Stdev        Skewness    Skew-Test(p-value)\n",
      "#---------------------------------------------------------------------------------------\n",
      "      6      0.0150425      0.0139688       0.697128    -0.00828963       0.761954\n",
      "      7     0.00316282   -1.45686e-05       0.730729      0.0254228       0.352964\n",
      "      8      0.0972618      0.0554929       0.721932       0.137939    5.24893e-07\n",
      "      9     0.00258965    -0.00323408       0.719454      0.0350835       0.199971\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat co2-100K.stat.dat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f847ce1",
   "metadata": {},
   "source": [
    "The second file prints the the probability density of the 7,8 & 9th normal modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "112e64cd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#Column-1 => Grid points(Freq Scaled)\n",
      "#Other columns => probability densities for Modes:6 7 8 9\n",
      "-4.995  0.0  0.0  0.0  0.0\n",
      "-4.985  0.0  0.0  0.0  0.0\n",
      "-4.975  0.0  0.0  0.0  0.0\n",
      "-4.965  0.0  0.0  0.0  0.0\n",
      "-4.955  0.0  0.0  0.0  0.0\n",
      "-4.945  0.0  0.0  0.0  0.0\n",
      "-4.9350000000000005  0.0  0.0  0.0  0.0\n",
      "-4.925  0.0  0.0  0.0  0.0\n",
      "-4.915  0.0  0.0  0.0  0.0\n",
      "-4.905  0.0  0.0  0.0  0.0\n",
      "-4.895  0.0  0.0  0.0  0.0\n",
      "-4.885  0.0  0.0  0.0  0.0\n",
      "-4.875  0.0  0.0  0.0  0.0\n",
      "-4.865  0.0  0.0  0.0  0.0\n",
      "-4.855  0.0  0.0  0.0  0.0\n",
      "-4.845  0.0  0.0  0.0  0.0\n",
      "-4.835  0.0  0.0  0.0  0.0\n",
      "-4.825  0.0  0.0  0.0  0.0\n",
      "-4.8149999999999995  0.0  0.0  0.0  0.0\n",
      "-4.805  0.0  0.0  0.0  0.0\n",
      "-4.795  0.0  0.0  0.0  0.0\n",
      "-4.785  0.0  0.0  0.0  0.0\n",
      "-4.775  0.0  0.0  0.0  0.0\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "head -n 25 co2-100K.prob-density.dat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85848f3d",
   "metadata": {},
   "source": [
    "This file could be used to plot the probability density using matplotlib. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6842c1a4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp4AAAHOCAYAAAAxEAn5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAACnoUlEQVR4nOzdd3xN9//A8ddNZCIRYoVI7E3UiFlVKaW0tlKltH6tWpVOrTaNVqkus3zNTrUatGrv2jWrJGqPIBIjiYSse35/HG5dGe6Ne3PueD8fjzy4557xTsj7vO/5LJ2iKApCCCGEEEJYmYvWAQghhBBCCOcghacQQgghhCgQUngKIYQQQogCIYWnEEIIIYQoEFJ4CiGEEEKIAiGFpxBCCCGEKBBSeAohhBBCiAIhhacQQgghhCgQUngKIYQQQogCIYWn0JxOp+Pjjz/WOgwhhAOS/CKEbZHCUwDw3XffodPp0Ol0bN++Pdv7iqIQGBiITqejU6dOGkSYf127dqVPnz6A+n34+fnx3XffmXy8Xq9nxowZhISE4OXlRYkSJXjyySc5fPiwlSIWwrFIfsnZvZ9JTl9PPfWUFaMWQjuFtA5A2BZPT08WLFhAy5YtjbZv3bqVixcv4uHhoVFk+bd3717eeecdAKKjo7l58yZNmzY1+fhBgwbx888/079/f4YNG0ZKSgoHDx7k6tWr1gpZCIck+cXYjz/+mG3bvn37mDx5Mu3atbNonELYCik8hZGOHTuyZMkSpkyZQqFC//33WLBgAQ0bNiQhIUHD6Mx38eJFLl26ZLgR7Nq1C19fX6pXr27S8YsXL+b7778nKiqKrl27WjNUIRye5Bdj/fr1y7Zty5Yt6HQ6w1NUIRyNNLULI3369OHatWusX7/esC09PZ2lS5fSt2/fHI9JSUnhzTffJDAwEA8PD6pXr86XX36JoihG+6WlpTFq1ChKlixJ0aJFefbZZ7l48WKO54yNjWXQoEGULl0aDw8Pateuzbx580z6HtLS0khISCAhIYHNmzfj5uZGYGAgCQkJbNu2jXr16nHt2jUSEhLQ6/V5nuvrr7+mSZMmdO3aFb1eT0pKikkxCCGyk/zy8HP/+uuvtG7dmvLly5t1rBB2QxFCUZT58+crgPLXX38pzZs3V1588UXDe8uXL1dcXFyU2NhYJSgoSHnmmWcM7+n1euXJJ59UdDqd8sorryjTpk1TOnfurADKG2+8YXSNfv36KYDSt29fZdq0aUq3bt2UevXqKYASERFh2O/KlStK+fLllcDAQGXs2LHKjBkzlGeffVYBlG+++cbk78WUrzNnzuR6nsTEREWn0ylDhw5VRo8erRQpUkQBlIoVKyqLFi0y+WcrhLOT/GKaqKgoBVBmz55t1nFC2BMpPIWiKMY3hmnTpilFixZVUlNTFUVRlJ49eypt2rRRFEXJdmNYvny5Aiiffvqp0fl69Oih6HQ65eTJk4qiKMqhQ4cUQHn99deN9uvbt2+2G8PLL7+slC1bVklISDDa9/nnn1d8fX0NceXm0qVLyvr165X169crQUFBSv/+/ZX169crv/zyiwIoU6ZMMbx/+/btXM9z4MABBVBKlCihlC5dWvn222+Vn3/+WWnSpImi0+mU1atX5xmHEEIl+cU03bt3Vzw8PJQbN26YdZwQ9kQKT6EoivGN4erVq0qhQoWUxYsXK0lJSYqXl5fhE/iDN4b/+7//U1xdXZWkpCSj8+3atUsBlKlTpyqKoiifffaZAigxMTFG++3du9foxqDX65VixYop//d//6fEx8cbfd2Lcfv27SZ9Tzdu3FBcXFyUtWvXKoqiKEuWLFE8PT2VO3fumHT8tm3bDE8udu/ebdienJys+Pv7Ky1atDDpPEI4O8kvD5eYmKh4enoqXbt2zdfxQtgLGVwksilZsiRhYWEsWLCA1NRUsrKy6NGjR477njt3joCAAIoWLWq0vWbNmob37/3p4uJC5cqVjfZ7sBN+fHw8N2/eZNasWcyaNSvHa+Y1mjwjI4PExEQA1q5di4uLCzVq1CAhIYG1a9fSoEEDkpOTSU5OxtfXFzc3t1zP5eXlBUDFihUJDQ01bC9SpAidO3fmp59+IjMz02iQhBAib5Jfcvbrr79y584dXnjhBZOPEcIeyR1T5Khv374MHjyYK1eu0KFDB4oVK1Yg173XGb9fv34MGDAgx33q1auX6/E7duygTZs2RtuCgoKMXpcsWRKAzZs388QTT+R6roCAAABKly6d7b1SpUqRkZFBSkoKvr6+uZ5DCJGd5Jfsfv75Z3x9fe1uHlMhzCWFp8hR165defXVV9m9ezeLFi3Kdb+goCA2bNhAcnKy0VOJmJgYw/v3/tTr9Zw6dcroKcTx48eNzndvRGpWVhZhYWFmx12/fn3DiNkhQ4bQtGlTBgwYQGJiIj169GDy5MnUqlXLsG9eAgICKFOmDLGxsdneu3TpEp6entmexAghHk7yi7HLly+zefNmXnrpJbucy1QIc8h0SiJHRYoUYcaMGXz88cd07tw51/06duxIVlYW06ZNM9r+zTffoNPp6NChA4DhzylTphjtN2nSJKPXrq6udO/enV9//ZV//vkn2/Xi4+PzjNvPz4+wsDBatmzJ+fPn6d69O2FhYRQuXBhXV1defvllwsLCCAsLw8/PL89zAfTu3ZsLFy4YTf+SkJDAihUrePLJJ3FxkV8hIcwl+cXYwoUL0ev10swunII88RS5yq0p6n6dO3emTZs2fPDBB5w9e5b69euzbt06VqxYwRtvvGHocxUSEkKfPn349ttvSUxMpHnz5mzcuJGTJ09mO+eECRPYvHkzoaGhDB48mFq1anH9+nUOHDjAhg0buH79+kPj2rdvH+np6TRv3hyAnTt3Uq9ePQoXLmzWz2D06NEsXryY7t27Ex4ejq+vLzNnziQjI4PPPvvMrHMJIf4j+eU/P//8MwEBAWY1zQtht7Qe3SRsw/2jTvPy4KhTRVFHeY8aNUoJCAhQ3NzclKpVqypffPGFotfrjfa7ffu2MmLECKVEiRJK4cKFlc6dOysXLlzINt2JoihKXFycMnToUCUwMFBxc3NTypQpo7Rt21aZNWuWSd/PhAkTlMqVKxteh4WFKUOHDjXp2AedOnVK6dq1q+Lj46N4eXkpTz75pLJ37958nUsIZyT5JXcxMTEKoISHh+freCHsjU5RHlj+QQghhBBCCCuQDmpCCCGEEKJASOEphBBCCCEKhBSeQgghhBCiQGheeE6fPp3g4GA8PT0JDQ1l7969ee5/8+ZNhg4dStmyZfHw8KBatWqsWrWqgKIVQtgbyTFCCGE7NJ1OadGiRYSHhzNz5kxCQ0OZNGkS7du35/jx45QqVSrb/unp6Tz11FOUKlWKpUuXUq5cOc6dO1dgq14IIeyL5BghhLAtmo5qDw0NpXHjxobJgfV6PYGBgQwfPpz33nsv2/4zZ87kiy++ICYmxqw1cIUQzklyjBBC2BbNCs/09HS8vb1ZunQpXbp0MWwfMGAAN2/eZMWKFdmO6dixI8WLF8fb25sVK1ZQsmRJ+vbty7vvvourq2uO10lLSyMtLc3wWq/Xc/36dUqUKIFOp7P49yWEM1MUheTkZAICAjRf1UlyjBCOxZbyi8g/zZraExISyMrKonTp0kbbS5cubViH90GnT59m06ZNvPDCC6xatYqTJ0/y+uuvk5GRQURERI7HjB8/nsjISIvHL4TI3YULFyhfvrymMUiOEcIx2UJ+EflnV0tm6vV6SpUqxaxZs3B1daVhw4bExsbyxRdf5HpTGD16NOHh4YbXiYmJVKhQgQsXLuDj41NQoYt8GD/eF4DRoxM1jkSYKikpicDAQIoWLap1KPkiOcax3cspeZF8Y7vsPb8IlWaFp7+/P66ursTFxRltj4uLo0yZMjkeU7ZsWdzc3IyavGrWrMmVK1dIT0/H3d092zEeHh54eHhk2+7j4yM3BRvn6an+mde/U2Sk2pQZESELcNkSW2hilhwjHnQvp+RF/s1sny3kF5F/mnWScHd3p2HDhmzcuNGwTa/Xs3HjRpo1a5bjMS1atODkyZPo9XrDtn///ZeyZcvmeEMQQjgvyTFCCGF7NG1qDw8PZ8CAATRq1IgmTZowadIkUlJSGDhwIAD9+/enXLlyjB8/HoAhQ4Ywbdo0Ro4cyfDhwzlx4gSfffYZI0aM0PLbEBq496RTiLxIjnFu0iLiXPR6PXfu3NE6DKfk7u5OoUKmlZSaFp69e/cmPj6ejz76iCtXrhASEsKaNWsMgwHOnz9vNHItMDCQtWvXMmrUKOrVq0e5cuUYOXIk7777rlbfghDChkmOEcI5pKWlcezYMaPWClGw/P39qVChwkO7Qmg+uGjYsGEMGzYsx/e2bNmSbVuzZs3YvXu3laMSQjgKyTFCODZFUTh79iyFChWiYsWKMtVSAdPr9dy6dYvY2FgAgoKC8txf88JTiBxlZhJ4DspcAb74AurV47OtT5PhIc1mQohHV/walLsIxW7C9RJwsTwkFtM6KpEfGRkZ3Lp1i4oVK1KkSBGtw3FK937usbGxFC5cGH9//1z3NbvwjIiIYNCgQQ+taIXIt1Wr4O23GXQMMl2BP8fCrVuMKAxb2sDYCJ2Gw+KENUl+EdbmewPaboS6/6ivU73A+zbodXDwMeC1OHhg7ldh2zIzMwFynF1CFJx7xeeaNWt44okncp1r1ezb94oVK6hcuTJt27ZlwYIFRit2CGGqyEhd9gFCigIffADPPAOlSjFvEIx/H0hMhGPHOFUZOq2EXougUIYmYQsrk/wirKnCOXhtJgSfhd+ehQnvwRfvwsR3YH07qHUUaNAAjhzROlSRDzLNkrbudXFISkpi3bp1ZGTkfKM2u/A8dOgQf/31F7Vr12bkyJGUKVOGIUOG8Ndffz1axMK5KQoMGQKffaY2rW/axIUKoHeFyE9ciVxci+Xd4OcXoPIp6PcjuElN4nAkvwirWb2afj/C5bIwfZj6dDPt7ryet71hdzOYPhT1aefjj8OePZqGK4S98vHxITk5mZSUlBzfz1eDZYMGDZgyZQqXLl1i7ty5XLx4kRYtWlCvXj0mT55MYqKs/CDM9OWX8L//seJZiEx5G3L55HqyKvwwAMpehmd/A6S7p8OR/CIs7tgx6NWL05XUD69puUwkn1IU2LIFatXiVtumfP2mPEETwlwuLi4oioKi5HyDfqTBRYqikJGRQXp6Ooqi4Ofnx7Rp0/jwww+ZPXs2vXv3fpTTC2excSO89x689x6HPCc8dPeLgbC8C/RaApcCYFcL64coCp7kF2EJHneArl0hKIhfnz1KlttDDvD1hagosqqVoedi4LM0kL6Ddsma8z2bOsg1ODiYq1ev4uLiQtGiRenVqxdfffUVVapU4aeffqJly5YAHDlyhA8//JCtW7eiKAqBgYH06dOH8PBwPD09+fjjj7l48SJz5swxOv/Zs2epWLEihQsXNtq+a9cuo4UyUlJS8Pb2NnRHOHbsGBUqVHiUH0G+5euJ5/79+xk2bBhly5Zl1KhRNGjQgOjoaLZu3cqJEycYN26cTLgsTHPzJvTrB23bwqefmnxYdG3Y0QLCNkDpy9YLTxQ8yS/CktqvBa5cgWXLyDC1fixdmsW91JYVPvnEitEJZ7Bu3Tpu3brFtm3bWLx4MbNnzzZ6Pzo6mhYtWlC7dm2io6O5efMmy5cvJz4+ngsXLjz0/K6urty6dcvoq27dukavPTw8OHr0qOG1VkUn5KPwrFu3Lk2bNuXMmTPMnTuXCxcuMGHCBKpUqWLYp0+fPsTHx1s0UOGgRo+GlBSYNw/uWx/bFJvaQIK/OuCIrCzrxCcKlOQXYUlBZ6DBQWDiRKha1axjL5WH7S3vHnv0qFXiE86latWqtGrViqMP/H+KjIwkLCyMcePGUaZMGQAqV67MN998Q1Uz/9/aA7MLz169enH27Fn++OMPunTpgmsOxYK/v7+sHiAeqvx5YOZMVrdMJnJuoNnNIvpCsLIzlI8FZsywSoyiYEl+ERaTlkanlXA+EBg8OF+n2N4KqFQJXn0V5P+ceETHjx/nzz//JCQkxGj75s2bee6557QJSgNmF573+lo96Pbt24wdO9YiQQknoMDTayA2AP5qnP/TXKgA+xsCH36oNtsLuyb5RVjMt99S/Lr64TTyE9d89ffLKgTfNT0OO3bAkiVWCFI4gw4dOlCsWDE6dOjASy+9xKBBg4zev3btmuFJJ0B4eDjFihWjcOHC/Pjjjw89f1ZWFsWKFTP6yrLhVkCzC8/IyEhu3bqVbXtqaiqRkZEWCUo4vhoxUO4SbHgKlEecDH7zE0B6ujoNk7Brkl+ERSQnw2efcbABxJd6tFOdqwj/VoVrrz8PucxLKEReVq9ezc2bNzl9+jTjx4/PtqRn8eLFuXLliuH1119/zc2bN2nVqpVJBaSrqys3b940+sqptchW5OuJZ06TtB4+fJjixYtbJCjh4LKyaLMJTlWCsxUf/XQpRYGRI2HSJHUQgbBbkl+ERXz9NSQns7W1+YfmtLjFprZQ4jrw3XcWCU+I+7Vp04bffvtN6zAKjMmFp5+fH8WLF0en01GtWjWKFy9u+PL19eWpp56iV69e1oxVOIpffqFUPGxsm/suZjeLvf02uLuzu0tZq06hIaxD8ouwmBs34KuvYNgwkn0tc8q4MnCkDhAZqbauCGFBH330EevWreOjjz4iLi4OUKdJio2NNdovKyuLO3fuGL5yWxnI1pk8j+ekSZNQFIVBgwYRGRmJr+9/v9Hu7u4EBwcbzRklRI4UBSZO5N+qcLmcBc/r5wcjRvDY+LFsy8dTDqEtyS/CYmbMUIvDt9+GmV9Z7LTbHoe638bCggXw0ksWO6+wHlPn2tRa7dq12b59Ox9++CHVq1cHoEKFCvTt25cePXoY9vvuu+/47r6n7m3btmXOnDlkZWUZ1km/Z+nSpTz99NMFEr+5TC48BwwYAEDFihVp3rw5bm4Pm4VXiBysXQtHjrDzJSuce9gwdJ+NpZGsrmh3JL8Ii7hzB6ZMUQvD0qUteuqEUkDnzuoqa/37g8sjdk4XTuHs2bMmba9fv36eze0ff/wxH3/8cY7v5bZC0P3u3Lnz0H0KikmFZ1JSEj4+PoC6nN3t27e5fft2jvve20+IHE2cCE2acC5or+XPXbIkhxpA6B7g9m3w8rL8NYTFSX4RFvPjj3D1Krz5pnXO/8470KoVrF4NzzxjnWsI4eBMKjz9/Py4fPkypUqVolixYjl2/r83KMCWh/ALjR04AJs3s7gnYKVumLuaQcN9qDeg//s/61xEWJTkF2ERiqL27eza1ezJ4k3WogU0bao+9ZTCU4h8Manw3LRpk2FE6ebNm60akHBg334LgYHE1Hj4EmD5daM4/FsNanz7rTppdA5FjLAtkl+ERWzaBMePw6xZ1ruGTqfOoNGnD9OH6hg63T76EAphS0wqPFu3bp3j34Uw2c2baqf8Dz5AyRxj1Uvtaww1fjoMu3eDDEixeZJfhEXMmAG1a6tN4dbUrRu3CkOjfda9jBCOyuze0WvWrGH79u2G19OnTyckJIS+ffty48YNiwYnHMj336uTL7/8stUvdaoSULmy+oRV2BXJLyJfYmNh+XIYMsT6rRzu7hx8DOofBnJY7EAIkTezC8+3336bpKQkAI4cOUJ4eDgdO3bkzJkzhIeHWzxA4QAURX0a0b073LcsmNW4AK+9BosXQ3y89a8nLEbyi8iXOXPA0xNefNGql7k3ufz+huCejtqKI4Qwi9mF55kzZ6hVqxYAv/76K507d+azzz5j+vTprF692uIBCgewYwccP873nosKbHL3iVffVgven3/OdR+ZaN72SH4RZtPrYd486NsXCmjWg8RicLIKMHdugVxPCEdiduHp7u5OamoqABs2bKBdu3aAutbovScVQhiZPx+CgzkbZL1LPLjM3e3CwLPPqtc2YY4zYRskvwizbdoE58/DwIEFetlDDYC9e+HYsQK9rhD2zuzCs2XLloSHh/PJJ5+wd+9enrk7pcS///5L+fLlLR6gsHMpKaT9NI8twWfz8b/tEQ0cCH//DQcPFvCFRX5JfhFmmz8fqldXpzkqQMerASVKqNcXwoG89tprTJw40WrnN7sUmDZtGoUKFWLp0qXMmDGDcuXUdQ9Xr15ts8szCQ0tXYpHOhyqr8G127eHsmXlxmBHJL8Is9y8CVFR6ofMAp46TV8IeOEFdc5gO10zW1hfcHCw0YBJUFcheuWVVzSK6OFmzpzJO++8Y7Xzm7xk5j0VKlRg5cqV2bZ/8803FglIOJjvv+d0RUj00+DahQqpgw3mzFEnlnZ31yAIYQ7JL8IsS5ao67K/+KKhq02Brs89cKC6ROfatdCpU8FdVzidzMxMChUyu2SzSflq/NTr9fz7779s376dbdu2GX0JYRAbC1u2cKSehjH06wfXr8O6dRoGIcwh+UWYbMECaNsWAgK0uX5IiDp36C+/aHN9Yfe6detGqVKlKF68OD179uT69euAupZ7oUKFmDlzJuXKleOll17ipZdeYsSIEbRp04YiRYrQuXNn4uPj6dGjBz4+PrRp04Zr164Zzh0VFUXNmjXx8/OjU6dOxMbGGp179uzZlC1bljJlyvD9998bjnvppZf49NNPDa8XLVpEnTp1KFq0KHXr1uX48eOP9D2bXT7v3r2bvn37cu7cuWwL08uSdsLIkiXg5kZ0jXTtYqhbV70xLFwoTyTsgOQXYbJLl2DrVu1HlvfpA+PHQ2oqeHtrG4v4T2oqxMRY9xo1ajzyv3m3bt348ccfyczMpHfv3owdO5ZJkyYBkJWVxaFDhzh16hSKojBkyBCWLFnChg0bqFChAs2bN6dly5bMmTOHBQsW0KlTJ6ZMmUJkZCQxMTEMHDiQP/74g8aNG/P222/Tr18/w+pwWVlZ/PPPP5w7d44tW7bQrVs3unXrRtGiRY3i27FjB0OHDmXFihU0a9aMf//9F59HnD3C7MLztddeo1GjRvzxxx+ULVs2x3WVhQDUpwAdOpDmtULbOJ5/HiZMkBuDHZD8Iky2eDG4ualrs2upd28YM4YlAwpzrE4BN/WL3MXEQMOG1r3G/v3w2GMP3a1Dhw64uroaXt+5c4d+/foBGP4EGDVqFB988IHRsREREXh6ehpe9+zZk9q1awPQsWNHoqOjaXV3ta6uXbuydu1aAJYsWUKXLl1o2bIlAJ999hl+fn5cvnzZcK6PPvoId3d32rVrh7e3N6dOnSIkJMTo+t999x2vvvoqLVq0AKBGjRoP/X4fxuzC88SJEyxdupQqVao88sWFAzt9Wp1qZOFCiLGBwvPDD+GPP6BnT21jEXmS/CJMtnAhdOgAxYppFoKhX2njxtT55y+O1dEsFPGgGjXUwtDa1zDB6tWrDQUgqIOLLl68SGZmJm+99RbLli3jxo0bKIqCv7+/YT8XFxfKli1rdK5SpUoZ/u7l5ZXt9a27q2ldunSJChUqGN4rUqQIJUqU4NKlS5QoUQJXV1dKlChheN/b29tw7P0uXrxIaGioSd+nqcwuPENDQzl58qTcGETeFi5Uny526gRWbu14qCpVoHFj9QmsFJ42TfKLMMnp07Bnj+30rezTh6pv/4XHHa0DEQbe3iY9jdTSzz//zJYtW9i5cyflypVj7dq1vPrqq4b3H6XFJyAggJMnTxpep6SkcO3aNQICAkhLSzP5PIGBgZw9ezbfceTE7MJz+PDhvPnmm1y5coW6devi5uZm9H69elqOJBE245df4LnnoHBhrSNRPf88vP8+JCaCr6/W0YhcSH4RJrn3wbZz51x3KdCVyXr1wjU8nBrRBXdJYf+Sk5Px9PTEz8+PhIQEvvzyS4udu0ePHjRt2pSdO3fSqFEjxowZQ/PmzSlbtqxZheSAAQPo0qULzzzzDE2bNjX08XzwSaw5zC48u3fvDsCgQYMM23Q6HYqiSOd/ofrnH/Vr3DitI/lP797w1luwfDkMGKB1NCIXkl+ESRYuVFcms5UPtuXKcS4I6vyjdSDCnvTv358//viD0qVLU758eV555RVOnDhhkXPXrFmTOXPmMHDgQK5evUqzZs346aefzD5PixYtmDx5MoMGDeLixYtUrFiRJUuWPFLhqVMeHDr6EOfOncvz/aAgK66LaAFJSUn4+vqSmJj4yCOzRC7GjIHp0+HKFfDw0GxN9Gyd/Fu3Bi8vWLMGUJ+IyEAAy3rU3y97zy8gOcbqjh6FOnVgxQp49tls+eXe73RB5p2ICIWVnXR0XAUucVehZMkCu7Yzye13KzU1lejoaGrWrIm3DCDVzL1/h+PHj3P16lUGDBiAn1/2SbzNfuJpD4lfaEhR1KcR3bqBh4fW0Rjr0weGDYP4eLkx2CjJL+KhFi5UBxS1b691JAaRkTq8akGH1cDSpTBkiNYhCWGz8jWB/I8//kiLFi0ICAgwPKGYNGkSK1ZoPHpZaG/fPjh1Si3ybE2PHuqfS5dqG4fIk+QXkSsTPthGRuo0aWW5XRhOV8J2BjwJYaPMLjxnzJhBeHg4HTt25ObNm4Y+V8WKFTNMemqu6dOnExwcjKenJ6Ghoezdu9ek4xYuXIhOp6NLly75uq6wgqVLwd8fnnhC60iy8/eHsDCjwlOrm5TImeQXkafDh+HkSejVS+tIcnS0DrB9O9w3V6IQwpjZhefUqVOZPXs2H3zwgdGEqI0aNeLIkSNmB7Bo0SLCw8OJiIjgwIED1K9fn/bt23P16tU8jzt79ixvvfWWYeJUYQMUBaKioEsXdZ10W9Stm7rayX3LignbIflF5CkqCooV45PtT9vkh8bj1QAXF7X/qShwZg5ZERam1+uBh/87mF14njlzhgYNGmTb7uHhQUpKirmn4+uvv2bw4MEMHDiQWrVqMXPmTLy9vZk3b16ux2RlZfHCCy8QGRlJpUqVzL6msByj5H/0qPo0QuuVRPLy3HOg18Pvv2sdiciB5BeRp2XLoFMn9Db6ufaON2prz7JlWofiVArdfdBhzvyUwvLuTUCfkZGR535m//pWrFiRQ4cOZRsEsGbNGmrWrGnWudLT09m/fz+jR482bHNxcSEsLIxdu3bletzYsWMpVaoUL7/8Mn/++Wee10hLSzP6z5iUlGRWjMIMUVGkuYNH27ZaR5IjwyojLVqqT06svJqaMJ+95ReQHFNg/v0X/vmHRbVte86iPzw28vQacL1xA3IY0Sssz83NjSJFihAbG4u7uzsuLvkaviLySa/Xc+vWLWJjY7l58yaZmZl57m924RkeHs7QoUO5c+cOiqKwd+9efvnlF8aPH8+cOXPMOldCQgJZWVmULl3aaHvp0qWJicl5uZvt27czd+5cDh06ZNI1xo8fT2RkpFlxiXyKiuLfalDX1kazP6hbNxg9Gvc6kG7joTobe8svIDmmwCxbRkYhOFlZ60DyFlMDnlkFrFwJL76odThOQafTERwczLFjxzh+/LjW4TitmzdvEhcXR0ZGBoUKFcLd3T3H/cwuPF955RW8vLwYM2YMqamp9O3bl4CAACZPnszzzz//yIHnJTk5mRdffJHZs2cbrWeal9GjRxMeHm54nZSURGBgoLVCdF6nT8Phw0T3hKh7TxZtdY7Mrl0hPJwqJzCsrRxp6zE7CXvLLyA5psBERXGiKmTmfC+zGbd84EJ5CIyKksKzAHl4eFC/fn1iYmLYsWMHmZmZeHp6ah2WU1AUhbS0NLKyssjKyiI1NZWQkJBc51TNV0+ZF154gRdeeIHU1FRu3bpltEi9Ofz9/XF1dSUuLs5oe1xcHGXKlMm2/6lTpzh79iyd71sm7V5n1kKFCnH8+HEqVzb+OOzh4YGHrT+BcwTLloGHByer2E4fm1wHHgQHQ4MG1Ig5aCg8he2wp/wCkmMKxMWLsHcvMTbcffx+MTWhzMrlfPGBjvfHyYfZguLi4kKtWrVQFIXDhw+TkpIiA44KiKurK66urri7u1O3bl1atWqV61rzj9RF29vb+5FWCXB3d6dhw4Zs3LjRMGWJXq9n48aNDBs2LNv+NWrUyDaydcyYMSQnJzN58mR5yqClqCho354Mj9+0jsQ03bpR7ZODuGZAltvDdxcFT/KLMFi+HAoV4t9qefcdsxXRNeCp9VD5pNaROKfatWtTu3ZtrcMQuTCp8GzQoEGuleuDDhw4YFYA4eHhDBgwgEaNGtGkSRMmTZpESkoKAwcOBNS1TMuVK8f48ePx9PSkTh3jR1TFihUDyLZdFJwiycDOnfDdd3DWTgrPrl3x+PBDKp2BE9W0Dsa5SX4RDxUVBW3bkua1VutITHKjBMSVgprRWkcihO0xqfC8fwLlO3fu8O2331KrVi2aNWsGwO7duzl69Civv/662QH07t2b+Ph4PvroI65cuUJISAhr1qwxDAg4f/68jFCzcTViQK+DL469BPayTG6tWiSUUG8M9xee0tez4El+EXlKSFDn3p0xAy7bR+EJEF0Tmu4G0tMhl0EWQjgjnWJmB4hXXnmFsmXL8sknnxhtj4iI4MKFC3nOj2cLkpKS8PX1JTExER8fH63DsXuRkTr6/aD+/af+xu/dK95sbZLne3Ftb6WjwQH46k1QXHPeR5jnUX+/7D2/gOQYi5s/H2XQIL56C1KKaB2M6UpdgSEzgTVrbGpdeXsmv1uOweyP+kuWLKF///7Ztvfr149ff/3VIkEJ++GZCsFn1SlE7E10TSicChXOax2JuEfyi8gmKorzFeyr6AS4Whqu+6F2ExBCGJhdeHp5ebFjx45s23fs2CFTFzihav+Cq94+C89LAZDoI/2wbInkF2EkORnWrSPGvLUDbINOHd3O8uWQlaV1NELYDLNHtb/xxhsMGTKEAwcO0KRJEwD27NnDvHnz+PDDDy0eoLBtNaPVOetu5dDqYWtN7Nno1IK5Rgys6aC+FtqS/CKMrFoF6elE2+EHW1BbVZrvvKoOvmzVSutwhLAJZhee7733HpUqVWLy5Mn89NNPANSsWZP58+fTq1cviwcobFhKCpVPweY2WgeSf9E1IXQvBFyCS+W0jkZIfhFGli2Dxx4j0c+82QxsxcVykFwE/nnncZrtkn7jQkA+5/Hs1auX3ASclNGo7zVrcMvEPpvB7jpfAVK91Ce3UnjaBskvAoA7d+CPP+DddyHLPgtPXNRWlZrRgKKAidOGCeHIZB4RkX9RUVwpDTeKax1I/imucLwG1IgG5IGEELZjwwa4dQu6ddM6kkcSXROKJQIHD2odihA2QQpPkT/p6bBypV0OKnpQdA3wvwb+8VpHIoQwWLYMqlWDmnbcpAKcC4bbnqjfjxBCCk+RPz8N8oCkJKLt+54AwOlKkOYuo9uFsBmZmbBihfq0086bp/WucLw6Mq2SEHdJ4SnypWaMOkfd1dJaR2K+yEid0Yj7LDc4UVUd3S6E0N53r7jBtWvQtavWoVhETA3g2DH491+tQxFCc2YXnps3b7ZGHMKO6PRQPeZuMrXvhxEGMTUg4DL43tQ6Eucm+UWA+sE2qSjQqJHWoVjEqcqAl5c0twtBPgrPp59+msqVK/Ppp59y4cIFa8QkbFz5i1Akxb5Hsz/oRFXIdJWnnlqT/CJQlP8+2Lo4RqNcpjvQoYMUnkKQj8IzNjaWYcOGsXTpUipVqkT79u1ZvHgx6enp1ohP2KAa0XCrsDpxvKNI91T7etaQfp6akvwiOHCAYomO9cEWULsN7NkDsbFaRyKEpswuPP39/Rk1ahSHDh1iz549VKtWjddff52AgABGjBjB4cOHrRGnsBWKOghHfRqhdTCWFVNDXbfdOyV7P1BRMCS/CKKiSPWCc0FaB2JhzzwDhQqpS2gK4cQeqXR47LHHGD16NMOGDePWrVvMmzePhg0b0qpVK44ePWqpGIUNKR0Hfjcd8GkEd0eeAtWPaxuHUEl+cVLLlvFvNXU0uEN9APTzgzZtpLldOL18FZ4ZGRksXbqUjh07EhQUxNq1a5k2bRpxcXGcPHmSoKAgevbsaelYhQ2oEQ13POBMsNaRWF5qEXUlI2lu15bkF+cUGalj2nAdREc75AdbQJ0eassWuH5d60iE0IzZhefw4cMpW7Ysr776KtWqVePgwYPs2rWLV155hcKFCxMcHMyXX35JTIyM0nBENWJQn0bka7FV2xdTAyqdBvc7WkfinCS/OLca0YC3tzoK3BE99xzo9fD771pHIoRmzC48jx07xtSpU7l06RKTJk2iTp062fbx9/eXaVEc0alTlIlzzGb2e2JqQKEsqHpS60ick+QX51YjBujQgUw3rSOxvMhIHZQtC02bSnO7cGpmF54RERH07NkTDw8Po+2ZmZls27YNgEKFCtG6dWvLRChsx7JlZBSCk1W0DsR6Ev3gchlpbteK5BfnVTQRysdClPKr1qFYV9eusHYtpKRoHYkQmjC78GzTpg3Xc+ifkpiYSJs2bSwSlLBRy5ZxqjJkuGsdiHVF14SqJ8A1Q+tInI/kF+dVIwayXODfqlpHYmVdu8KdO2rxKYQTMrvwVBQFXQ5r5167do3ChQtbJChhg65cgV271GmUHFxMTfBIh0pntI7E+Uh+cV41YuBMRUjz0joSK6tSBerWlbXbhdMyeYhIt27dANDpdLz00ktGTWFZWVn8/fffNG/e3PIRCtuwYgW4uPBv9SytI7G6+JJwrbg0txckyS9O7to1gs/Cqo5aB1JAunaFyZMhPR3cHbwJSYgHmFx4+vr6AuoTiaJFi+Ll9d/HUnd3d5o2bcrgwYMtH6GwDcuWQevW3PbepHUk1qdTm9sbHASyssDVVeuIHJ7kFye3ciU6BY47QYsKoE6rNHasOrVSu3ZaRyNEgTK58Jw/fz4AwcHBvPXWW9Ls5Uxu3oSNG1n1VKbWkRSYmJrQcgewfTvIQBark/zi5KKiuBAIt4pqHUgBqVcPKlZUm9ul8BROJl+j2uWm4GR++w0yMx16GqUHxQZAog+wdKnWoTgVyS9OKCkJ1q51qvyCTgfdu6stSZnO84FeCDDxiedjjz3Gxo0b8fPzo0GDBjl2/r/nwIEDFgtO2IjFi6FlS5J9tmsdScFxgWO1oNnSpTBpkjS3W5HkFyf3+++QlsbRWloHUnAiI3UEJMLgq8C2bfDkk1qHJESBManwfO655wyd/bt06WLNeIStuXED1q2Dr76C605UeAJHa0Oz3Vekud3KJL84ucWLoVkzkort0jqSAnUpAAgOVr9/KTyFE9EpiqJoHURBSkpKwtfXl8TERHx8fLQOx/Z99x0MGgQXLxI5u5zW0RQsBUZOghPVYNUzEBHhVL8q+SK/X/IzyEtkpPo02/C7lJgIpUrB558TmThKw8i0EZH6DsybB5cvQyEHXYfYguR3yzGY3cdTOJnFi6FVKwgI0DqSgqdTm9trHgOdXutghHBAv/2mTinUo4fWkWijVy9ISFBHtwvhJEz6iOXn55dnv6v75bTqiLBT167B+vXqfHNO6mhtaL4Lgs5pHYnjkvzixBYtgpYtoXx5rSPRxmOPQaVK6gf8sDCtoxGiQJhUeE6aNMnKYQibtHw56PXqnHNO6lI5uOkLtY9qHYnjkvzipO7vP+6sdDr1qefs2TB9Ori5aR2REFZnUuE5YMAAa8chbNHixeqgmjJltI5EOzr1qWfIIdRpT6QflsVJfnFSK1aov1Pdu2sdibZ69YIJE2DzZpnTUzgFk+6iSUlJho68SUlJee4rHX4dREICbNwI06ZpHYnmjtaGFjuBrVuhbVutw3E4kl+c1N3+4043aPFBISHq+u2LF0vhKZyCSYOL/Pz8uHr1KgDFihXDz88v29e97cJBLFsGiuLUzez3XA6AG8VQbwzC4iS/OKHr19X+4716aR2J9u41t0dFQUaG1tEIYXUmPfHctGkTxYsXB2Dz5s1WDUjYiEWLOB2k58cZpWUaobvN7S1//RWmTgV3d60jciiSX5zQsmWQlaU2s/9vmNbRaMYwvVSvQ/DZZ2ox3rGjtkEJYWUmFZ6t75s8u7VMpO34Ll6ETZs48qzWgdiOI/Wg5Y5rsGYNPCs/GEuS/OKEfvxRHcXtzP3H71evHtSurf5cpPAUDi5fIyVu3LjB3LlziY6OBqBWrVoMHDjQ8NRC2LmffwZPT47VvK11JDbjamnUvlg//iiFp5VJfnFsvjdQ+0v/+KPWodgOnQ7694eICHXteunLLByY2RPIb9u2jeDgYKZMmcKNGze4ceMGU6ZMoWLFimzbts0aMYqCpChc/eo9jlS+Tbqn1sHYmBdfVCe8vnFD60gcluQXx1fvb0h3g8+iXzQ0NQugb19IS4OlS7WORAirMrvwHDp0KL179+bMmTNERUURFRXF6dOnef755xk6dGi+gpg+fTrBwcF4enoSGhrK3r17c9139uzZtGrVyjDoICwsLM/9hZkOHqRUPByur3UgNqhvX3X6FxlkZDWSXxycohaex2pBhnSVNla+vDprxg8/aB2JEFZlduF58uRJ3nzzTVxdXQ3bXF1dCQ8P5+TJk2YHsGjRIsLDw4mIiODAgQPUr1+f9u3bG0a5PmjLli306dOHzZs3s2vXLgIDA2nXrh2xsbFmX1vk4IcfuFUYTlf6b1NkpE6eTIDaH61dO2kitCLJL46tXCz4X4O/62kdiY168UW1G8I5WSpNOC6zC8/HHnvM0PfqftHR0dSvb/5jsq+//prBgwczcOBAatWqxcyZM/H29mbevHk57v/zzz/z+uuvExISQo0aNZgzZw56vZ6NGzfmuH9aWhpJSUlGXyIXGRmwYAFH6oHi+vDdnVL//rBjB5w6pXUkDsne8gtIjjFH/cOQVBTOVtQ6EhvVrRt4e8NPP2kdiRBWY9Lgor///tvw9xEjRjBy5EhOnjxJ06ZNAdi9ezfTp09nwoQJZl08PT2d/fv3M3r0aMM2FxcXwsLC2LVrl0nnSE1NJSMjI9eBB+PHjycyMtKsuJzWunUQH89hJ19IJE/PPQdFi6pPPT/+WOtoHII95xeQHGOy9HRq/wMHHgPF7EceTqJIEXWKqR9+gPffVwcdCeFgdIqiPHSSRhcXF3Q6HQ/bVafTkZWVZfLFL126RLly5di5cyfNmjUzbH/nnXfYunUre/bseeg5Xn/9ddauXcvRo0fx9Mw+GiYtLY20tDTD66SkJAIDA0lMTJRVUB7UuzdERxPZ/YjWkdgkw3ymgwapzWEnT8qN4QFJSUn4+vqa9ftlz/kFJMeYbPly6NqVb1+H+FJaB2ObIiIUdS7Pdu1g924IDdU6JJuSn/wibI9JTzzPnDlj7TjyZcKECSxcuJAtW7bkelPw8PDAw8OjgCOzQzdvqmsnf/oppLytdTS2rX9/mD8fdu6EFi20jsbu2XN+AckxJvvxRy6VlaLzoZ58EgIC1FYVKTyFAzKp8AwKCrLKxf39/XF1dSUuLs5oe1xcHGUeMrHwl19+yYQJE9iwYQP16klP9Ue2cKHax7NvX5gthWeeHn8cgoI4OLwlvz2HrOz0iCS/OIGrV+H33/m7rdaB2AFXV+jXD2bPhi++AC8vrSMSwqLyNYE8wLFjxzh//jzp6elG2581Y3Jtd3d3GjZsyMaNG+nSpQuAoSP/sGG5L6M2ceJExo0bx9q1a2nUqFG+4hf3URSYMQM6d1Y/aYu8ubjAK69QJ/JD1rXTOhjHJPnFsWzoW5rWChyWGv6hIiN1FL8Nw2+gTt02YIDWIQlhUWYXnqdPn6Zr164cOXLEqF+W7m5fN3P6YAGEh4czYMAAGjVqRJMmTZg0aRIpKSkMHDgQgP79+1OuXDnGjx8PwOeff85HH33EggULCA4O5sqVKwAUKVKEIkWKmPvtCIBdu+Dvv2HiRK0jsWmGdZUjFHjlFVwiPqTeYY2DcjCSXxxQVhYN98PR2nDHW+tgbNu9HHO9BJyqBO5jXiJQCk/hYMweWzhy5EgqVqzI1atX8fb25ujRo2zbto1GjRqxZcsWswPo3bs3X375JR999BEhISEcOnSINWvWULp0aQDOnz/P5cuXDfvPmDGD9PR0evToQdmyZQ1fX375pdnXFmqiOzykBVSuDE89pXU4diEyUkfk/8oSUwMa7UN9YiwsQvKLA1q7Fr+bsK+x1oHYl32NIfAicOiQ1qEIYVEmjWq/n7+/P5s2baJevXr4+vqyd+9eqlevzqZNm3jzzTc5ePCgtWK1CBkVZ2ziOzrCv4ZC4yfC22rfTpks3jTBZ2DA98DmzfDEE1qHYxMe9ffL3vMLSI7JpnNnLu9byaxXAUktJtNlwRuTwKfP/8H//qd1ODZBfrccg9lPPLOysihatCig3iQuXboEqAMEjh8/btnohNU1OHT3L3ebHoXpzgZDvD9q/1hhEZJfHMy5c/DHH+xrhBSdZlJc4UBD4OefQRYlEA7E7MKzTp06HD6sdmwLDQ1l4sSJ7Nixg7Fjx1KpUqWHHC1sil5Pw31q3yv8/bWOxv7oYH8jICoK7vYFFI9G8otj+fPFYNLcFI7U1ToS+3TgMdCnprCqj6/WoQhhMWYXnmPGjEGv1wMwduxYzpw5Q6tWrVi1ahVTpkyxeIDCitavp/gN1KcRIl8O1Qfc3GDuXK1DcQiSXxxIejoNDsDh+pAh05zmS7IPal/yv5C+5MJhmD2qvX379oa/V6lShZiYGK5fv46fn59h5KmwEzNmcKU0XAzUOhD7leYF9OkDs2bBe++pc/CJfJP84hgiI3XU/gd6pMigoke1rxH0/xHYvh1atdI6HCEe2SOtmHvhwgUuXLhA8eLF5aZgb06fht9/V28KOhlQ9EiGDoXz59UlAYXFSH6xb6G74WyQrFT0qM5UhIQSwOTJWocihEWYXXhmZmby4Ycf4uvrS3BwMMHBwfj6+jJmzBgyMjKsEaOwhm++geLFOVxf60AcwGOPcSYYLo7sIc1hj0jyi2Mof16dCmhXc60jcQAusKsZal/ykye1jkaIR2Z24Tl8+HBmzZrFxIkTOXjwIAcPHmTixInMnTuXESNGWCNGYUGRkTomvquDefNg2DAy3bSOyDHsag7lY4EdO7QOxa5JfnEMzXeqMz78W1XrSBzD3/VRB4B+843WoQjxyMzu47lgwQIWLlxIhw4dDNvq1atHYGAgffr0YYZMLWPzGv0F6PVqE/H0jw3bpbk9/05Ugasl4forrVjUR9Zvzy/JLw7g33+pEQO/d+YRO3OJezLdYHOdeFrO+pZJhb/l7YmSX4T9MjsteHh4EBwcnG17xYoVcXd3t0RMworc0iB0D/xV9w6R00tqHY5DiIzUqc1hzaHGcSgVl8s+4qEkvziAzz/nVhH4W9Zlt6i/GoGiU/vOCmHPzC48hw0bxieffEJaWpphW1paGuPGjWPYsGEWDU48mshIXbaCp9E+8LwD21toFJQD+7se3CgGrbZpHYn9kvxivyIjdUx+Qwc//MDOFpAl3Xgs6nZhdYaA0D3AjRtahyNEvpnU1N6tWzej1xs2bKB8+fLUr6+OTDl8+DDp6em0bdvW8hEKy0lNpflOONQAkoppHYzj0bvC9lbQ6Xfg2DGoVUvrkOyC5BfH0XI7ULw4+xpe1ToUh7SzOTTeC1u6FmdrG+nSI+yTSYWnr6/xqgndu3c3eh0YKBNB2oWZM/G6DX+21DoQx3WoPjy+FXw/+QR++UXrcOyC5BfH4HsTQg4Cn79NZsrbWofjkFKKqPN6Nt0Ne5pqHY0Q+WNS4Tl//nxrxyGsLTERPvuMQyGQ6Kd1MI5LXwi2tYbOCxfCu+9CSIjWIdk8yS+O4YnNcNsLptx4G6Q7rtXsaAEN90MLmUBD2Kl8jzmMj49n+/btbN++nfj4eEvGJKzhyy8hJYWtT2gdiOM7GAJUqwajR2sdit2S/GJnjhyh/mHY2hoypOi0qpSisLupOsjoqzez9+MXwtaZXXimpKQwaNAgypYty+OPP87jjz9OQEAAL7/8MqmpqdaIUTyiL9/Swddfw8iRJPtoHY3jU1xhcf1/Yc0a2LxZ63DsiuQXO/X++9zwgwMNtQ7EOexsARlu8MQWrSMRwnxmF57h4eFs3bqV33//nZs3b3Lz5k1WrFjB1q1befPNN60Ro3hEbTdCqpLK53yudShOI7oWXCwHV/o8ydgIeSJhKskvdmjdOli5ko1t1QF2wvrSPGHb49DgIJS+onU0QpjH7MLz119/Ze7cuXTo0AEfHx98fHzo2LEjs2fPZunSpdaIUTyCgIvQ4BBsehLueGkdjRPRweoOUCYOHtuvdTD2Q/KLncnIgJEj4fHHOVZb62Ccy19N4FoJ6LAKWapX2BWzC8/U1FRKly6dbXupUqWkKczW6KHDarhSWprAtHCpvNrf88lN4CW/GiaR/GJnpk2Df/+FKVNAHuwXKL0rrOkAQeeBhQu1DkcIk5ldeDZr1oyIiAju3Llj2Hb79m0iIyNp1qyZRYMT+XNv4vhG+9T1w1d1BEWWrtPExjBw0cNT67SOxD5IfrEj587BRx/BkCFwd85VUbBOV4boGsCoUTKpvLAbZq/VPmnSJJ5++ulsEzx7enqydu1aiwco8sfnJoRtUOd8uxCkdTTOK6UIrGsPz/4GrF8PTz2ldUg2TfKLnVAUeO01KFYMPvtM62ic2qqOUHPuHXjzTZg3T+twhHgoswvPunXrcuLECX7++WdiYmIA6NOnDy+88AJeXtKJ0CYo0GklpHnA+jCtgxEHG0CdI1Bp8GA4coTIr9WpBWTVkewkv9iJH35QZ21YuRJ8ZKoMLd3ygd8fT6Tz/Pn8mDGf01UktwjbZlbhmZGRQY0aNVi5ciWDBw+2Vkwin+7N59Z4L1Q9CT/3hXRPjYMSoIPfO8NrM88R3coHumodkG2S/GInTp+G4cPhxRfhmWe0jkYABx6DWkehy3KYMUTraITIm1k9/9zc3Iz6XgnbUyoO2q2DPU3gZDWtoxH33CyuNomFHFafforsJL/Yvk8+0kHfvuDvrw4sErZBB8u7gGuW2q0n8mMZ6SVsl9lDToYOHcrnn39OZmamNeIRj8DjDvRcrE6xsV66Etqcv+vDkTrQ6XcoIYvx5Ejyi21ruwH0e/cw58kzhib2e4MZhbZu+cCKLlDjODTbJf8uwnaZ3cfzr7/+YuPGjaxbt466detSuHBho/ejoqIsFpwwg15P1ygocgtmD4YsN60DEtnoYGVneHkOPL8QGJcIvr6Gm4P0y5L8YtMWLKD5LljzNMQGah2MyMm/1dW13MPWq9PonamsdURCZGd24VmsWDG6d+9ujVjEoxgzhmr/wi994bq/1sGI3KR7wKLnYfAsoE8fWLFC65BsiuQXG7VnD7zyCofqw55QdZM8TbNNG9uqqxn1WArzXtY6GiGyM7vwnD9/vjXiEI9gdUcdHVbD+nZwQvp12rzrJWBJL3jxl/UweDAEI5Nv3yX5xQYdP05qm6Yk+MMfnZD/qzZOcYFfe8DAedDvR/jKXcctH2lREbbD5D6eer2ezz//nBYtWtC4cWPee+89bt++bc3YhCl+/JGnV8OuZrCrudbBCFOdrgx89x18/z3t1wKKcz9BkvxS8EzqA3j6NLRrx60i8EsfyJQuPHbhjhf83A90CvT7CbxvaR2REP8xufAcN24c77//PkWKFKFcuXJMnjyZoUOHWjM28TDz58OAARxsAOtkMJH9eeEFmD6dprvh6dWAEz+QkPxig06ehCeeAHd3fuoHd7y1DkiYI8kXfnoRvFNhwPfAlStahyQEYEbh+cMPP/Dtt9+ydu1ali9fzu+//87PP/+MXq+3ZnwiJ4oCn38OgwbBq6/ye2fyMT+B0FpkpI7I+KH83glC90KXZfDph845ElXyi/bu/383+/903AqpCt7esHUryb4aBibyLaEkfP8SeN6B67XKMm248+UWYXtMLlfOnz9Px44dDa/DwsLQ6XRcunTJKoGJXNy5o/YLfO89+PBDIkvPlKLTzh1oBEt6QO2jap+swk7YLCb5xYYsWcJL8+GGH/DnnxAQoHVE4hFc84d5gyDLVZ1Rg40bgf+6WjjjB12hLZNLlszMTDw9jZfBcXNzIyMjw+JBiVycOAHNm8NPP6n9A8eOlY7+DuJYHfhhAPgnwKszgc2btQ6pQEl+0Z5rBuqKRL16EVND/f9IyZJahyUsINEP5r4MlwKAp56Cjz9Gl6V1VMJZmTyqXVEUXnrpJTw8PAzb7ty5w2uvvWY0157Ms2cFWVkwdSq8/z7XvG6z5CV4bcAAraMSFnahAsx8Dbr/CkXbtoVhw+Czz6BIEYef61Pyi8Z27eLV/wFJs2D6dKKuDpUPtQ4m7e6Ao1Z/KrQeG8nLZWHFcxBfWuvIhLMxufAckEOh069fP4sGI3KwdSuMGgWHDsGwYfyvyFQyPJx7BLQjSykKP/aHj/y+hvffh6goGD8e9Dh0lwrJL9ookozaV/y770gvC9++nE58/H9Fp+QZx6K4wLbWcLoSPLsCXv0f/NUEGHED/Py0Dk84CZ2iKI75CCUXSUlJ+Pr6kpiYiM/dJd9sjqLAtm0wbhysXw+hofDNN9CsmdwInEREhAJnzsDbb8OvvxJXCv58HHr8kgGFzJ5+t8DYxe+XldnFzyA2lt09ytNwH7j5loDISMbGDUNx1TowUVBcM6DZbmj5J3h4+cKIEXxx5xNSC9tuy4pd/G6Jh3LgZyh2KC4Ovv0WGjRQpzGJi4PFi2HnTmjWTOvoRAGKjNQR+UMlWLoUtm8nuai6EgmBgfDuuxATo3WIwgblOVgkLQ3++AN69oSgIOofUpdX5ORJGDpUik4nk+UG21vB1OGwq1Yi6RM+Ifwr6LEYWL1a7eJ1HxmIJCzFJgrP6dOnExwcjKenJ6GhoezduzfP/ZcsWUKNGjXw9PSkbt26rFq1qoAitTBFUQuIyZPh8cdRypRBP2woVKgAa9aozes9e4KLi/zSO7MWLfj5RbX/Jz17wuzZULOm+iQ8IgJ27IDMTK2jtFlOm18A4uPh11/VOWNLlYJOnbi6dSmr2mUx6Q3Y2gYoVkzjIIWWUorCuvYwaRSsf0od4EjHjiT5FWJ/Ix2Le+v4/D259wjL0bypfdGiRfTv35+ZM2cSGhrKpEmTWLJkCcePH6dUqVLZ9t+5cyePP/4448ePp1OnTixYsIDPP/+cAwcOUKdOnYdeT7NH9RkZ6iog0dFw7Bjs2qU+ybx+Hdzc4KmnWFFoFcerw+3/xlIYmjyk6BT3RLx7W13jfelSdWqUGzfAxweaNFGfljdoAPXrQ6VK8MBIcWuztaawgs4voM3PIPJjHUVuqUWDfzwEXILAC+B/7e4OdeuypcQRjtWE+FLIwCGROwXKXoL6f0Plk+r/Ib0OrpSBK2Xh8t0/r5WA214Q8XHBlRC2ll9E/mheeIaGhtK4cWOmTZsGqEvnBQYGMnz4cN57771s+/fu3ZuUlBRWrlxp2Na0aVNCQkKYOXPmQ6+X7/+4WVmQnq42V93/lZ4Ot27BzZvq140b6p/XrsGlS3DxIsTGwvnzavEJ4OurFgnNm0OLFtC0KRQtKsWlMIlR/6usLNi/X+0L/Ndf6lPyc+f+e79sWahYEcqXB39/46/ixaFwYfDyUicKv/9PLy9wdQWdef8nbe3GUND5BfL5M1AUo/zy9YRSuGbCyNeO8b9JtfBMgwHP/vpfjrlxQ12JJjaWuANr8E0EzzT1VFkucLUUXAhUZ0q4EAiJxUwLQ4gH+dyEyqfUDzJlrkCpq+B6d12HNHfwqFEPgoKgVCm2x8wl1Rva9Z2v5pgiRQw5ZcqcemQUgjc/vKHmF3d3u88vIn80HaWQnp7O/v37GT16tGGbi4sLYWFh7Nq1K8djdu3aRXh4uNG29u3bs3z58hz3T0tLIy0tzfA6MTERgKTAQPU//YN194OvMzPVm4GpK6i4uKiFpZ+fetMPCICQELX5vHp1pq7tzK0iiYwevZTx431hE+qXECYaPTrnZD36B/X/Ntev89MHFSl2EzrV7a8OUoqLU5+0X7sGCQn/fQh6GJ1OLUAf/HJxMf773RtI0t3fE1sYs1gQ+QXMzDEP/lz0erXgfODf45W7fyZNqUWfe3//vrt6PTe44wkphSG5KCQHQlItuF5CfQp1sxjoH+xEdSfX8IXI0x1PuFobdtVWX7tmQYkE8LsBxW5CscS/8Tn6N177oPxtdYnOpPUDs53npbt/Jn39wOj5vPLL/X9iW/lF5J+mhWdCQgJZWVmULm08kVjp0qWJyWXwxJUrV3Lc/0ou69COHz+eyMjIbNsDk5LyGfVD6PX/PZE4fTrX3SZMkDXohGXl+H9q1eePdlJFUT98mdmHNDk5GV9fbf+PF0R+AQ1yTMbdr2RAlt8W9i4rK9tApoexhfwi8s9252WxkNGjRxs9wdDr9Vy/fp0SJUqgM/MxvzUkJSURGBjIhQsXpOkgB/LzyZut/XwURSE5OZkAJ1pmUXKM/ZKfTd5s7efjjPnFEWlaePr7++Pq6kpcXJzR9ri4OMqUKZPjMWXKlDFrfw8PD6PVUACK2eAoTh8fH5v4xbZV8vPJmy39fGzlSURB5BeQHOMI5GeTN1v6+dhKfhH5p+l0Su7u7jRs2JCNGzcatun1ejZu3EizXOatbNasmdH+AOvXr891fyGEc5L8IoQQtkfzpvbw8HAGDBhAo0aNaNKkCZMmTSIlJYWBA9XOyf3796dcuXKMHz8egJEjR9K6dWu++uornnnmGRYuXMi+ffuYNWuWlt+GEMIGSX4RQgjbonnh2bt3b+Lj4/noo4+4cuUKISEhrFmzxtDB//z587i4/Pdgtnnz5ixYsIAxY8bw/vvvU7VqVZYvX27yHHu2xsPDg4iIiGxNdUIlP5+8yc8nb86eX0D+j+RFfjZ5k5+PsAbN5/EUQgghhBDOwSaWzBRCCCGEEI5PCk8hhBBCCFEgpPAUQgghhBAFQgpPIYQQQghRIKTwtCFnz57l5ZdfpmLFinh5eVG5cmUiIiJIT0/XOjTNTJ8+neDgYDw9PQkNDWXv3r1ah2QTxo8fT+PGjSlatCilSpWiS5cuHD9+XOuwhA2T/JKd5JecSX4R1iSFpw2JiYlBr9fzv//9j6NHj/LNN98wc+ZM3n//fa1D08SiRYsIDw8nIiKCAwcOUL9+fdq3b8/Vq1e1Dk1zW7duZejQoezevZv169eTkZFBu3btSElJ0To0YaMkvxiT/JI7yS/CmmQ6JRv3xRdfMGPGDE6fPq11KAUuNDSUxo0bM23aNEBddSYwMJDhw4fz3nvvaRydbYmPj6dUqVJs3bqVxx9/XOtwhJ2Q/CL5xRSSX4QlyRNPG5eYmEjx4sW1DqPApaens3//fsLCwgzbXFxcCAsLY9euXRpGZpsSExMBnPL/isg/yS+SX0wh+UVYkhSeNuzkyZNMnTqVV199VetQClxCQgJZWVmGFWbuKV26NFeuXNEoKtuk1+t54403aNGihV2vsCMKluQXyS+mkPwiLE0KzwLw3nvvodPp8vyKiYkxOiY2Npann36anj17MnjwYI0iF/Zg6NCh/PPPPyxcuFDrUIQGJL8Ia5L8IixN87XancGbb77JSy+9lOc+lSpVMvz90qVLtGnThubNmzNr1iwrR2eb/P39cXV1JS4uzmh7XFwcZcqU0Sgq2zNs2DBWrlzJtm3bKF++vNbhCA1IfjGf5BfTSH4R1iCFZwEoWbIkJUuWNGnf2NhY2rRpQ8OGDZk/fz4uLs75UNrd3Z2GDRuyceNGunTpAqhNPhs3bmTYsGHaBmcDFEVh+PDhLFu2jC1btlCxYkWtQxIakfxiPskveZP8IqxJCk8bEhsbyxNPPEFQUBBffvkl8fHxhvec8VN4eHg4AwYMoFGjRjRp0oRJkyaRkpLCwIEDtQ5Nc0OHDmXBggWsWLGCokWLGvql+fr64uXlpXF0whZJfjEm+SV3kl+ENcl0Sjbku+++yzXpOes/07Rp0/jiiy+4cuUKISEhTJkyhdDQUK3D0pxOp8tx+/z58x/a7Cqck+SX7CS/5Ezyi7AmKTyFEEIIIUSBcM4OPkIIIYQQosBJ4SmEEEIIIQqEFJ5CCCGEEKJASOEphBBCCCEKhBSeQgghhBCiQEjhKYQQQgghCoQUnkIIIYQQokBI4SmEEEIIIQqEFJ5CCCGEEKJAyFrtQhSAhIQESpYsabXzywJkQjg3yTHCXkjhKUQBWLVqFUlJSRQtWlTrUIQQDkhyjLAX0tQuRAGQG4IQwpokxwh7IYWnjXviiSd44403bPZ8juLatWuUKlWKs2fPWvzc6enpeHh4WPy89zz//PN89dVXVju/cGySYwqG5BghVFJ4auSll15Cp9Oh0+lwc3OjdOnSPPXUU8ybNw+9Xm/YLyoqik8++UTDSO2XOTfAcePG8dxzzxEcHGzYtm3bNjp37kxAQAA6nY7ly5fnevzAgQMZM2ZMju/9+eeftGzZ0ozIzTNmzBjGjRtHYmKi1a4h7I/kGOuTHCOE+aTw1NDTTz/N5cuXOXv2LKtXr6ZNmzaMHDmSTp06kZmZCUDx4sWl+cTKUlNTmTt3Li+//LLR9pSUFOrXr8/06dPzPD4rK4uVK1fy7LPP5vh+TEwMNWvWtFi8D6pTpw6VK1fmp59+sto1hH2SHGMbJMcI8R8pPDXk4eFBmTJlKFeuHI899hjvv/8+K1asYPXq1Xz33XdA9k/US5cupW7dunh5eVGiRAnCwsJISUkx7Dts2DCGDRuGr68v/v7+fPjhh3mORlyzZg0tW7akWLFilChRgk6dOnHq1CnD+3q9nokTJ1KlShU8PDyoUKEC48aNM3p//PjxVKxYES8vL+rXr8/SpUuNrvHEE08wfPhw3njjDfz8/ChdujSzZ88mJSWFgQMHUrRoUapUqcLq1avNPu+IESN45513KF68OGXKlOHjjz8G1Kc9W7duZfLkyYanPrk1ca1atQoPDw+aNm1qtL1Dhw58+umndO3aNdefH8DOnTtxc3OjcePGee73oCFDhuT6lKJ8+fJMmDDB5HN17tyZhQsXmnV94fgkx0iOkRwjbI0UnjbmySefpH79+kRFRWV77/Lly/Tp04dBgwYRHR3Nli1b6Natm1HS//777ylUqBB79+5l8uTJfP3118yZMyfX66WkpBAeHs6+ffvYuHEjLi4udO3a1dAUN3r0aCZMmMCHH37IsWPHWLBgAaVLlzYcP378eH744QdmzpzJ0aNHGTVqFP369WPr1q1G1/n+++/x9/dn7969DB8+nCFDhtCzZ0+aN2/OgQMHaNeuHS+++CKpqalmn7dw4cLs2bOHiRMnMnbsWNavX8/kyZNp1qwZgwcP5vLly1y+fJnAwMAcfwZ//vknDRs2fMi/TO5+++03OnfujE6nAzB8DwBnzpyhUqVK2Y45evQos2bNYuLEiTmes2bNmhw6dMjkGJo0acLevXtJS0szL3jhdCTHSI4ByTFCQ4rQxIABA5Tnnnsux/d69+6t1KxZU1EURWndurUycuRIRVEUZf/+/QqgnD17NsfjWrdurdSsWVPR6/WGbe+++67hXA+eLyfx8fEKoBw5ckRJSkpSPDw8lNmzZ+e47507dxRvb29l586dRttffvllpU+fPkbXbNmypeF1ZmamUrhwYeXFF180bLt8+bICKLt27cr3eRVFURo3bqy8++67Jn2v9zz33HPKoEGD8twHUJYtW5bje1WrVlVWrlyppKenKyNGjFA+/fRTw3tz585VUlNTsx0zYMAAJTQ0NNfr9erVS2nduvVDY7/n8OHDef7fEM5HcozkGMkxwhbJE08bpCiK4ZPt/erXr0/btm2pW7cuPXv2ZPbs2dy4ccNon6ZNmxod26xZM06cOEFWVlaO1zpx4gR9+vShUqVK+Pj4GDq+nz9/nujoaNLS0mjbtm2Ox548eZLU1FSeeuopihQpYvj64YcfjJrSAOrVq2f4u6urKyVKlKBu3bqGbfeecFy9ejXf5wUoW7YsV69ezTHe3Ny+fRtPT0+zjrknOjqaS5cu0bZtW9zc3Pjkk09Yv3694f3U1FS8vLyMjsnMzCQqKoru3bsbtr366qvMnTvX8Do5OTnbcXm5t+/9T0KEyI3kGMkxkmOEVmQCeRsUHR1NxYoVs213dXVl/fr17Ny5k3Xr1jF16lQ++OAD9uzZk+P+pujcuTNBQUHMnj2bgIAA9Ho9derUIT09/aFJ6datWwD88ccflCtXzui9B6f2cHNzM3p9b6Tt/a9B7Xf1qOe9f8SuKfz9/bPdXE3122+/8dRTTxluKj4+Pnh6epKQkEDhwoXx9vbOdsypU6dITk423BT1ej1Lliwxuvn+/fff9O7d2+Q4rl+/DmDVlUuE45AcIzlGcozQijzxtDGbNm3iyJEjRp9U76fT6WjRogWRkZEcPHgQd3d3li1bZnh/z549Rvvv3r2bqlWr4urqmu1c165d4/jx44wZM4a2bdtSs2ZNo+RYtWpVvLy82LhxY46x1KpVCw8PD86fP0+VKlWMvnLr62QKS53X3d0916cw92vQoAHHjh3LV6wrVqzgueeeM9rWoUMHVq1axZYtW3jiiSeyHXPz5k0AihQpAsDatWu5ceOG4caye/duYmNjDQMOrl27RpUqVQD1SUZISEi2wRz//PMP5cuXx9/fP1/fh3AekmMkx0iOEVqSJ54aSktL48qVK2RlZREXF8eaNWsYP348nTp1on///tn237NnDxs3bqRdu3aUKlWKPXv2EB8fbzSNxvnz5wkPD+fVV1/lwIEDTJ06NdeJf/38/ChRogSzZs2ibNmynD9/nvfee8/wvqenJ++++y7vvPMO7u7utGjRgvj4eI4ePcrLL79M0aJFeeuttxg1ahR6vZ6WLVuSmJjIjh078PHxYcCAAfn6uVjqvMHBwezZs4ezZ89SpEgRihcvjotL9s9a7du3Z/To0dy4cQM/Pz/D9lu3bnHy5EnD6zNnznDo0CGKFy9OhQoVuHr1Kvv27eO3334zOt8zzzzDBx98wOOPP06HDh2yXS8oKAidTscvv/xC4cKFeeutt3jmmWdYsWIFgYGBvPbaa4SFhRlGo5YoUcLwtKFQoUKUKlWKixcvGt0g//zzT9q1a2fSz0U4D8kxOZMcIzlGaEjbLqbOa8CAAQqgAEqhQoWUkiVLKmFhYcq8efOUrKwsw373d14/duyY0r59e6VkyZKKh4eHUq1aNWXq1KlG+77++uvKa6+9pvj4+Ch+fn7K+++/bzQQ4MHO8OvXr1dq1qypeHh4KPXq1VO2bNli1Mk9KytL+fTTT5WgoCDFzc1NqVChgvLZZ58Zjtfr9cqkSZOU6tWrK25ubkrJkiWV9u3bK1u3bs31moqiKEFBQco333xjtO3+6+b3vM8995wyYMAARVEU5fjx40rTpk0VLy8vBVDOnDmT679HkyZNlJkzZxpt27x5s+Hf6P6ve+efM2eO0qJFixzPFxoaavRv86DPPvtM8fHxUUqXLq3MnTtXOXTokBIUFKQULlxYef7555Xr168b7V+hQgXl5s2bSlZWllK9enUlKSnJ8N7t27cVX19fZdeuXbleTzgfyTGSYyTHCFskhacDMXWEpchu5cqVSs2aNY1uyA/TuXNn5fPPP8/xvfDwcGXDhg2WCk9p2rSpEh0drUydOlUZOnSo0Xvffvut8tRTT1nsWkLkRnJM/kmOEUIlTe1CoDZdnThxgtjYWJP7eLVs2ZI+ffrk+N7w4cMJCAiwWHxly5Zl8eLFbNiwgTVr1hi95+bmxtSpUy12LSGE5UmOEUKlU5Q8lpwQduWJJ54gJCSESZMmaR2KsLChQ4eyfft2Nm7cKJ37hWYkxzguyTGioEjhKYQQQgghCoRMpySEEEIIIQqEFJ5CCCGEEKJASOEphBBCCCEKhBSeQgghhBCiQEjhKYQQQgghCoQUnkIIIYQQokBI4SmEEEIIIQqEFJ5CCCGEEKJASOEphBBCCCEKhBSeQgghhBCiQEjhKYQQQgghCoQUnkIIIYQQokBI4SmEEEIIIQqEFJ5CCCGEEKJASOEpNKfT6fj444+1DkMI4YAkvwhhW6TwFAB899136HQ6dDod27dvz/a+oigEBgai0+no1KmTBhHmX9euXenTpw+gfh9+fn589913Jh+/ePFimjZtSrFixShRogStW7fmjz/+sFK0QjgeyS+5mzZtGjVr1sTDw4Ny5coRHh5OSkqKlaIVQntSeAojnp6eLFiwINv2rVu3cvHiRTw8PDSI6tHs3buXpk2bAhAdHc3NmzcNrx9m6tSp9O7dG39/fyZMmMCHH35IYmIinTp1IioqypphC+FwJL8Ye/fddxk+fDh16tRh8uTJdO/enalTp9KtWzdrhiyEpgppHYCwLR07dmTJkiVMmTKFQoX++++xYMECGjZsSEJCgobRme/ixYtcunTJcCPYtWsXvr6+VK9e3aTjp06dSuPGjfn999/R6XQADBo0iHLlyvH999/LDUIIM0h++c/ly5f5+uuvefHFF/nhhx8M26tVq8bw4cP5/fff6dy5s9ViF0Ir8sRTGOnTpw/Xrl1j/fr1hm3p6eksXbqUvn375nhMSkoKb775JoGBgXh4eFC9enW+/PJLFEUx2i8tLY1Ro0ZRsmRJihYtyrPPPsvFixdzPGdsbCyDBg2idOnSeHh4ULt2bebNm2fS95CWlkZCQgIJCQls3rwZNzc3AgMDSUhIYNu2bdSrV49r166RkJCAXq/P81xJSUmUKlXKUHQC+Pj4UKRIEby8vEyKRwihkvzyn127dpGZmcnzzz9vtP3e64ULF5oUjxD2Rp54CiPBwcE0a9aMX375hQ4dOgCwevVqEhMTef7555kyZYrR/oqi8Oyzz7J582ZefvllQkJCWLt2LW+//TaxsbF88803hn1feeUVfvrpJ/r27Uvz5s3ZtGkTzzzzTLYY4uLiaNq0KTqdjmHDhlGyZElWr17Nyy+/TFJSEm+88Uae38Mvv/zCwIEDjbaVK1fO6HXJkiUBOHPmDMHBwbme64knnmDp0qVMnTqVzp07c+fOHaZOnUpiYiIjR47MMw4hhDHJL/9JS0sDyPYB1tvbG4D9+/fnGYcQdksRQlGU+fPnK4Dy119/KdOmTVOKFi2qpKamKoqiKD179lTatGmjKIqiBAUFKc8884zhuOXLlyuA8umnnxqdr0ePHopOp1NOnjypKIqiHDp0SAGU119/3Wi/vn37KoASERFh2Pbyyy8rZcuWVRISEoz2ff755xVfX19DXLm5dOmSsn79emX9+vVKUFCQ0r9/f2X9+vXKL7/8ogDKlClTDO/fvn07z3PFxcUpbdu2VQDDl7+/v7Jz5848jxNC/EfyS3b79+9XAOWTTz4x2r5mzRoFUIoUKZJnHELYK2lqF9n06tWL27dvs3LlSpKTk1m5cmWuzWCrVq3C1dWVESNGGG1/8803URSF1atXG/YDsu334NMFRVH49ddf6dy5M4qiGJq0EhISaN++PYmJiRw4cCDP+MuWLUtYWBiNGjXiwoULvPDCC4SFhVGoUCE8PT35v//7P8LCwggLC8PT0zPPc3l7e1O9enUGDBjAkiVLmDdvHmXLlqVbt26cPHkyz2OFENlJflE99thjhIaG8vnnnzN//nzOnj3L6tWrefXVV3Fzc+P27dt5xiGEvZKmdpFNyZIlCQsLY8GCBaSmppKVlUWPHj1y3PfcuXMEBARQtGhRo+01a9Y0vH/vTxcXFypXrmy034Od8OPj47l58yazZs1i1qxZOV7z6tWrucaekZFBYmIiAGvXrsXFxYUaNWqQkJDA2rVradCgAcnJySQnJ+Pr64ubm1sePwno2bMnhQoV4vfffzdse+6556hatSoffPABixYtyvN4IYQxyS//+fXXX+nduzeDBg0CwNXVlfDwcLZu3crx48fzPFYIeyWFp8hR3759GTx4MFeuXKFDhw4UK1asQK57rzN+v379GDBgQI771KtXL9fjd+zYQZs2bYy2BQUFGb2+1/9q8+bNPPHEE7me6/Tp06xZsybbDap48eK0bNmSHTt25HqsECJ3kl9U5cqVY/v27Zw4cYIrV65QtWpVypQpQ0BAANWqVXvYtyOEXZLCU+Soa9euvPrqq+zevTvPp3pBQUFs2LCB5ORko6cSMTExhvfv/anX6zl16pTRU4gHP9XfG5GalZVFWFiY2XHXr1/fMGJ2yJAhNG3alAEDBpCYmEiPHj2YPHkytWrVMuybl7i4OACysrKyvZeRkUFmZqbZ8QkhJL88qGrVqlStWhWAY8eOcfnyZV566SWz4xPCHkgfT5GjIkWKMGPGDD7++OM855Lr2LEjWVlZTJs2zWj7N998g06nM4xcvffng6NWJ02aZPTa1dWV7t278+uvv/LPP/9ku158fHyecfv5+REWFkbLli05f/483bt3JywsjMKFC+Pq6srLL79s6H/l5+eX57mqVKmCi4sLixYtMpq65eLFi/z55580aNAgz+OFEDmT/JIzvV7PO++8g7e3N6+99prZxwthD+SJp8hVbk1R9+vcuTNt2rThgw8+4OzZs9SvX59169axYsUK3njjDUOfq5CQEPr06cO3335LYmIizZs3Z+PGjTkO0JkwYQKbN28mNDSUwYMHU6tWLa5fv86BAwfYsGED169ff2hc+/btIz09nebNmwOwc+dO6tWrR+HChU3+/kuWLMmgQYOYM2cObdu2pVu3biQnJ/Ptt99y+/ZtRo8ebfK5hBDGnD2/AIwcOZI7d+4QEhJCRkYGCxYsYO/evXz//fdUqFDBrHMJYTc0HFEvbMj9053k5cHpThRFUZKTk5VRo0YpAQEBipubm1K1alXliy++UPR6vdF+t2/fVkaMGKGUKFFCKVy4sNK5c2flwoUL2aY7URR1GqOhQ4cqgYGBipubm1KmTBmlbdu2yqxZs0z6fiZMmKBUrlzZ8DosLEwZOnSoScfeLyMjQ5k6daoSEhKiFClSRClSpIjSpk0bZdOmTWafSwhnJfklZ/Pnz1fq16+vFC5cWClatKjStm1byS3C4ekU5YHlH4QQQgghhLAC6eMphBBCCCEKhBSeQgghhBCiQEjhKYQQQgghCoTmhef06dMJDg7G09OT0NBQ9u7dm+f+N2/eZOjQoZQtWxYPDw+qVatmWC5NCCEeJDlGCCFsh6bTKS1atIjw8HBmzpxJaGgokyZNon379hw/fpxSpUpl2z89PZ2nnnqKUqVKsXTpUsqVK8e5c+cKbNULIYR9kRwjhBC2RdNR7aGhoTRu3NgwObBerycwMJDhw4fz3nvvZdt/5syZfPHFF8TExDx0DVwhhJAcI4QQtkWzwjM9PR1vb2+WLl1Kly5dDNsHDBjAzZs3WbFiRbZjOnbsSPHixfH29mbFihWULFmSvn378u677+Lq6prjddLS0khLSzO81uv1XL9+nRIlSqDT6Sz+fQnhzBRFITk5mYCAAFxctO3JIzlGCMdiS/lF5J9mTe0JCQlkZWVRunRpo+2lS5c2rMP7oNOnT7Np0yZeeOEFVq1axcmTJ3n99dfJyMggIiIix2PGjx9PZGSkxeMXQuTuwoULlC9fXtMYJMcI4ZhsIb+I/LOrJTP1ej2lSpVi1qxZuLq60rBhQ2JjY/niiy9yvSmMHj2a8PBww+vExEQqVKjAhQsX8PHxKajQxUOMH++b63ujRycWYCTiUSQlJREYGEjRokW1DiVfJMc4l5zyjuQb22Xv+UWoNCs8/f39cXV1JS4uzmh7XFwcZcqUyfGYsmXL4ubmZtTkVbNmTa5cuUJ6ejru7u7ZjvHw8MDDwyPbdh8fH7kp2BBPz9zfk38n+2MLTcySY8TD5JR35N/M9tlCfhH5p1knCXd3dxo2bMjGjRsN2/R6PRs3bqRZs2Y5HtOiRQtOnjyJXq83bPv3338pW7ZsjjcE4RgiI3VERkqiEeaRHCMkdwhhezTtnRseHs7s2bP5/vvviY6OZsiQIaSkpDBw4EAA+vfvz+jRow37DxkyhOvXrzNy5Ej+/fdf/vjjDz777DOGDh2q1bcghLBhkmOcV14FpxSkQmhH0z6evXv3Jj4+no8++ogrV64QEhLCmjVrDIMBzp8/bzRyLTAwkLVr1zJq1Cjq1atHuXLlGDlyJO+++65W34LIp3tJPyJCs9m8hBOQHCPyS3KUENah+eCiYcOGMWzYsBzf27JlS7ZtzZo1Y/fu3VaOSgjhKCTHCCGE7ZCJsIRN02WBa6bWUQghHJIChdK1DkII52L2E8+IiAgGDRpEUFCQNeIRAl0WNDgEj+2HUlfBNQsS/AF9BLz1FshUGg5L8osoCAEX4fFtUP4iFE6FpKJwvgJsba11ZEI4PrOfeK5YsYLKlSvTtm1bFixYYLRihxCPqkQ8vDYTOv8OST6wsS2segYulgcmToSqVWHtWq3DFFYi+UVYVVoavPIKg+dAsZvwV2NY3gUOh0DAJRgyA3jnHbhvVgMhhGWZXXgeOnSIv/76i9q1azNy5EjKlCnDkCFD+Ouvv6wRn3AiARdh0Dz177P+DxY/D3uawf5G8PtzwPHj8Nhj0KkTLFigaazCOiS/CKtJTlZzx08/8Xsn+N9rsLWNWnRuagvfDlX/5KuvoG9fXKSLjxBWka8+ng0aNGDKlClcunSJuXPncvHiRVq0aEG9evWYPHkyiYmy8oMwj3889P8BrpWA+QPhckAOO1WoAL/9Bv36wQsvwPLlBR2mKACSX4Sljf1Ix5l6PrB3L6xdy4FGoDxw98sqBDtaAkuWwLJldFkOyIB2ISzukQYXKYpCRkYG6enpKIqCn58f06ZNIzAwkEWLFlkqRuHokpLovRASfeGnF+GOdx77FioE8+ZB9+7Qvz/ksua2sH+SX0R+5DRHZ9gGCDoH33VLInLLE3mfoFs3+Okn6v4DTXdZL04hnFW+Cs/9+/czbNgwypYty6hRo2jQoAHR0dFs3bqVEydOMG7cOEaMGGHpWIWjeuUVityCRb0hPfvKg9lEjnWB+fOhfHn1JnHnjvVjFAVG8ouwpBrR0HwXrGsP54JNPKhnT3a0gKfWAzt2yITzQliQ2YVn3bp1adq0KWfOnGHu3LlcuHCBCRMmUKVKFcM+ffr0IT4+3qKBCse08HkdLFnCyk5w3d+MA4sWhaVL4eRJGDfOavGJgiX5RViSx23o+AfEVIc9oeYdu7EtxJYDBg+WKd2EsCCzC89evXpx9uxZ/vjjD7p06YKrq2u2ffz9/Y3WOhYiJ+53oOMq+LcqHK2TjxPUqgXvvw8TJsDRoxaPTxQ8yS/CksI2gHs6rOoImPnAUnGB3zsDJ07QYrs1ohPCOZldeN7ra/Wg27dvM3bsWIsEJZxD663geUedLsmUm0KOzV2jR0PlyvD666DISAB7J/lFWEpALDTarz65TPbN3zniSwNvv02rP6HYdYuGJ4TTMrvwjIyM5NatW9m2p6amEhkZaZGghOPzSYQme2FHC0gs9ggn8vCAb76BbdtgzRpLhSc0IvlFWErbDXC1JOxrbN5x2T7gjhnDbS94YotFwxPCaeXriadOl/3x1OHDhylevLhFghKOr/UWSPOA3c3yd7zRzeHpp6FVK7XZXZpg7ZrkF2EJFU9DpTPqvJwPTptkNm9vtrWGen8D//xjifCEcGomL5np5+eHTqdDp9NRrVo1o5tDVlYWt27d4rXXXrNKkMLBHD9OyCF1lKkpo9gfSqeDzz5Ti88lS6B3bwucVBQkyS/CYhR4ciNcLAfHq1vmlAcaQLOdUHzMGJk/WIhHZHLhOWnSJBRFYdCgQURGRuLr+1+nGXd3d4KDg2nWLJ+Pr4Rz+fxzbhWBfQ0f/VT3nnpGRCjqk89x46BXL7UYFXZD8ouwlEqnoXws/PgiZg8oyo2+kLqOe9flK5jxuo4h30p/ciHyy+TCc8CAAQBUrFiR5s2b4+bmZrWghAOLjYWffmJ3G8iy9H+h996DJ55Q13J/+mkLn1xYk+QXYSnNd8ClsnC6kmXP+09deHKT+uRTCJF/JvV+SUpKMvy9QYMG3L59m6SkpBy/hMjT5Mng5cX+x6xw7scfh8aNYeJEK5xcWIvkF2ExBw5Q+TTsbIHFnnbeo3dV+6TXPQJcuGDZkwvhREx64unn58fly5cpVaoUxYoVy7Hz/71BAVlZWRYPUjiIpCSYORNef510z88tempDk/s7S6BnT9i3Dxo1sug1hHVIfhEW88UX3CgGx2pa5/QHHoPHt8Kh3hVY1/5uFx8hhFlMKjw3bdpkGFG6efNmqwYkHNgPP8Dt2zB8OMyxbOFp0LUrVKwI06bBd99Z5xrCoiS/CIu4fBmWLmV3GCjZ1x2wiHQP2N8QGu2DzW2scw0hHJ1JhWfr1q1z/LsQJlMU+PZb6NIFypWz3nVcXeHVVyEiAr76CkqUsN61hEVIfhEWMWcOuLtzOMS661vuawQtdkAdmVlJiHwxe4azNWvWsH37f+uHTZ8+nZCQEPr27cuNGzcsGpxwIFu3QnS0usKQtQ0apBa68+db/1rCoiS/iHzJzIT//Q9eeIE0T+teKtEPTlSFRn8hq6UJkQ9mF55vv/22oZP/kSNHCA8Pp2PHjpw5c4bw8HCLBygcxIwZUKOGOurc2kqWVKdUmjlTJpS3M5JfRL6sXKnOmFEQH2yBvxpDwGWY/X8u2ZfxFULkyezC88yZM9SqVQuAX3/9lc6dO/PZZ58xffp0Vq9ebfEAhQOIj4eoKHjttYKbX3PIEDh1CjZtKpjrCYuQ/CLy5X//g9BQCAmx+KlzKixPVYEbxdS14IUQ5jG78HR3dyc1NRWADRs20K5dOwCKFy8u052InP30E7i4QL9+BXfNZs2genVpbrczkl+E2WJjYd06ePnlAruk4gKHQ6DWUXBLL7DLCuEQzC48W7ZsSXh4OJ988gl79+7lmWeeAeDff/+lfPnyFg9Q2Ll7fS2ffbZAB/pEjnVR+3pGRTHhPZ00h9kJyS/CbD/8AB4eBb5U7qEQcE+HWscK9LJC2D2zC89p06ZRqFAhli5dyowZMyh3d4Ty6tWreVpWixEPOnAAjhxRi8CC9uKLkJFBnaMFf2mRP5JfhFnufbDt0QN8fAr00onF4ExFCDlYoJcVwu6ZvGTmPRUqVGDlypXZtn/zzTcWCUg4mPnzISAA7jaZFqiyZeHppwk5+Af7ZS55uyD5RZhl5044cULt46mBQw2gWxRqf/LKlTWJQQh7Y3bhCaDX6zl58iRXr15F/8Co4ccff9wigQkHkJEBCxeqfa9c/5vRuUCbvQcMoPwff1D8WsFdUjwayS/CZD/9xE1fmLzlSSLaFPzURtE1IM0dPBYsgA8/LPDrC2GPzC48d+/eTd++fTl37hzKA3OYyZJ2wsiGDXDtGvTpo10MnTqR5i6TPdsLyS/CZBkZsGQJR+sALgX8gfauTHeIqQH1f/kFxowpuFk7hLBjZvfxfO2112jUqBH//PMP169f58aNG4av69evWyNGYa8WLlTn7qxfX7sYvLyIqQF1jiCTPdsByS/CZBs3wrVrHKmjbRj/1EFdHOPIEW0DEcJOmP3E88SJEyxdupQqVapYIx7hKG7fhmXL4M03DU8BtBpZfrQO1P8b9cZQr54mMQjTSH4RJvvlF6henbgyxzUN43QloHhxNR7JL0I8lNlPPENDQzl58qQ1YhGOZNUqSE5mWsLHWkfCqUqQ6oV6YxA2TfKLMMm9D7Z9+oDGrdv6Qqij6hculFYVIUxg9hPP4cOH8+abb3LlyhXq1q2Lm5ub0fv15BOfAFi4kMtl4Jq/1oGoN4boWtBw4UL47DN1jk8gIkJuErZG8oswyerVkJwMzz8PCz/WOhq1AJ41C/bsgaZNtY5GCJtmduHZvXt3AAbdNy+jTqdDURTp/C9USUmwciX/tNI6kP/8Uwcafn9WvTEImyX5RZjkl1+gQQN1dTJb0KqVOn3bL79I4SnEQ5hdeJ45c8YacQhH8ttvcOcO/9TWOpD/nAtCvTEsXAh+WkcjciP5RTxUcjKsXAmRkVpHYhD5aSEier+h5pevvzaaPk4IYczswjMoKMgacQhH8ssv0KIFScV2aB2JgeKCuqTewoXo/k99HRmpk+Z2GyP5RTzUihVw506BL5H5UH36wKRJsHUrPPmk1tEIYbPMHlwE8OOPP9KiRQsCAgI4d+4cAJMmTWLFihUWDU7YoWvXYN06te+VxiIjH1ij/fnn4coVgs5qFpIwgeQXkadffoHmzcHWPqQ0bgyVKskgRiEewuzCc8aMGYSHh9OxY0du3rxp6HNVrFgxJk2aZOn4hL1Ztgz0eujZU+tIsmvSBCpWzHHt9mxFqtCE5BeRpxs3YN06Vvvu1DqS7HQ69cPtr7+qk9sLIXJkduE5depUZs+ezQcffIDrff1YGjVqxJF8TqA7ffp0goOD8fT0JDQ0lL1795p03MKFC9HpdHTp0iVf1xVWEBWldrQvXVrrSLLT6aB7d6rHgE7/8N1FwZP8IvK0ciVkZhJdU+tActG9u1ocb9umdSRC2CyzC88zZ87QoEGDbNs9PDxISUkxO4BFixYRHh5OREQEBw4coH79+rRv356rV6/medzZs2d56623aNXKhoZOO7vERHWZzG7dtI4kd926USQFyl/QOhCRE8kvIk9RUdC0Kck+WgeSiwYN1C4AUVFaRyKEzTK78KxYsSKHDh3Ktn3NmjXUrGn+x9Cvv/6awYMHM3DgQGrVqsXMmTPx9vZm3rx5uR6TlZXFCy+8QGRkJJUqVTL7msJK/vhDbWLq2lWzEB7aZB4aSnIRqBldcDEJ00l+EblKSYE1awwfbG2ye4xOp8Z3r8uRECIbswvP8PBwhg4dyqJFi1AUhb179zJu3DhGjx7NO++8Y9a50tPT2b9/P2FhYf8F5OJCWFgYu3btyvW4sWPHUqpUKV5++eWHXiMtLY2kpCSjL2ElUVFqB/vAQKPNNnWDcHEhpsbdwlMGtNsce8svIDmmwKxZo45m1/CDrUm6doXLl2XOYCFyYfZ0Sq+88gpeXl6MGTOG1NRU+vbtS0BAAJMnT+Z5M0cyJyQkkJWVRekH+gOWLl2amJiYHI/Zvn07c+fOzfGpSE7Gjx9PpA3N9+awbt9WVxMZM0brSB4quiY03gdlrmgdiXiQveUXkBxTYJYtg7p1oUoVrSPJW/PmUKqUGm+zZlpHI4TNydd0Si+88AInTpzg1q1bXLlyhYsXL5r8dOBRJCcn8+KLLzJ79mz8/U1bi3H06NEkJiYavi5ckM59VrFuHaSm2nT/zntPXs8Fw21PaW63VfaUX0ByTIFIT4fff7fp/GLg6gpduqgtQLJ2uxDZmP3E837e3t54e3vn+3h/f39cXV2Ji4sz2h4XF0eZMmWy7X/q1CnOnj1L586dDdv0d/vRFCpUiOPHj1O5cmWjYzw8PPDw8Mh3jMJEUVFQq5btLGGXB70rHK8ONaTwtGn2kF9AckyB2LRJXYrXhgvPe92JIiIUtbl91iw4cgTq1dM4MiFsi0mFZ4MGDdDpTOujd+DAAZMv7u7uTsOGDdm4caNhyhK9Xs/GjRsZNmxYtv1r1KiRbUqVMWPGkJyczOTJkwl8oG+hKCAZGdxe/AN/NYHHtY7FRDE1IeQwcPy4XRTLjkzyi3ioZcugcmW1qd0ePPkk+PiocUvhKYQRkwrP++exu3PnDt9++y21atWi2d3+K7t37+bo0aO8/vrrZgcQHh7OgAEDaNSoEU2aNGHSpEmkpKQwcOBAAPr370+5cuUYP348np6e1KlTx+j4YsWKAWTbLgrQ1q143YHoGrDZTpahPFUZ0t3AfdkyeO89rcNxapJfRJ6ysmD5cujfXx01bg/c3aFzZ7UlKCJC62iEsCkmFZ4R9/3ivPLKK4wYMYJPPvkk2z756dvUu3dv4uPj+eijj7hy5QohISGsWbPGMCDg/PnzuLjkqyuqKChRUdz0hStltQ7EdJlucKIqFJs6mjlpo7UOx6lJfhF52rkTrl616Wb2HHXrBj//DCdP2v6AKCEKkE5RzOv97Ovry759+6hatarR9hMnTtCoUSMSExMtGqClJSUl4evrS2JiIj4+tjoLsR3R66F8eXYFXWbd0+qme088bWYKpVzU+Ru6R8E3oyDJV91mD09rbdmj/n7Ze34ByTEWN2oULFoEFy/C3Q8Jtp5bACLeugX+/jB2LLz9ttbhOAT53XIMZn/U9/LyYseOHdm279ixA09PT4sEJezInj1w+bLtLmGXhxPVIMsFauQ8s47QgOQXYURR1H6SXbsaik67UbgwPP20Gr8QwsDsUe1vvPEGQ4YM4cCBAzRp0gSAPXv2MG/ePD788EOLByhsXFQUlCrFxcC8lyC0RWmecLqSOrp9b6jW0QiQ/CIecPAgnDtnmDTeHp50GunWTe2beukSBARoHY0QNsHswvO9996jUqVKTJ48mZ9++gmAmjVrMn/+fHr16mXxAIUNUxS18OzSBcVlltbR5EtMTXhmJXinQGphraMRkl+Ekago8POD1q21jiR/OnWCQoXUwVH5GBwnhCPK1zyevXr1kpuAgL//htOn1U/1u+208KwOnX6H6sfh4GNaRyNA8ou4T1QUPPssuLlpHUn++PlxqkImui+HUkkKTyGAfK5cJASg9l3y9YU2bYw229Ta7A+RWgTOBclk8kLYnJgYiI62v9HsD4iuCcFngevXtQ5FCJsghafIv6gotSnJ3V3rSB5JTA2odBrc72gdiRDCYNky8PaGp57SOpJHcrwG6BTgt9+0DkUImyCFp8ifkyfhyBEWZ/xsN083cxNTEwplQdUTWkcihAC11SR22vvQoQN4eWkdziO5VRQuBCKj24W4SwpPkT/LloGnJycdYF7kxGJwqaxMqySErfBJhHKXsPtm9ntiagDr1kFKitahCKE5swvPzZs3WyMOYW+WLYP27cmw71Z2g5gad5943pH2di1JfhGRkTpqxKhz7PLMM1qHYxExNVFzy5o1WocihObMLjyffvppKleuzKeffpqvJeyEA7h8GXbtMsyt5whiaoJHOrBxo9ahODXJLwLUwX5nKqIOXrRT9w+yvFEc4kohze1CkI/CMzY2lmHDhrF06VIqVapE+/btWbx4Menp6daIT9iiFSvA1RU6d9Y6EouJLwkJJVAHTAnNSH4RXqkQdO7uU0IHEl0TWLkS5P+ycHJmF57+/v6MGjWKQ4cOsWfPHqpVq8brr79OQEAAI0aM4PDhw9aIU9iSZcvgiSegeHGtI7EcndrcnrJwHmM/su/BUvZM8ovzuveEsNpxdRR4THX7mprtYWJqAomJIN1JhJN7pMFFjz32GKNHj2bYsGHcunWLefPm0bBhQ1q1asXRo0ctFaOwJTduwKZNDtXMfk9MTSicChWkhdcmSH5xTjWj1VHgKUW1jsSy4koDFStKc7twevkqPDMyMli6dCkdO3YkKCiItWvXMm3aNOLi4jh58iRBQUH07NnT0rEKW/DHH5CZCV26aB2JxcUGQFJRmUxea5JfnJdbGlQ+5XjN7ADoUD+wL18OWVlaRyOEZsxeMnP48OH88ssvKIrCiy++yMSJE6lTp47h/cKFC/Pll18SEBBg0UCFjVi2DJo0gXLltI7E8lzU5vYaMajr0Osco4nPnkh+cW5VTqpz6kbX0DoSK+nWDb7+GnbvhhYttI5GCE2YXXgeO3aMqVOn0q1bNzw8PHLcx9/fX6ZFcUS3b6vTgXz4odaRWE1MTWjyF3DgADRsqHU4Tkfyi3OrEQNXSsPNB7qPO0o/T5o1g9Kl1Q/wUngKJ2V2U3tERAQ9e/bMdlPIzMxk27ZtABQqVIjWrVtbJkJhO9atg9RUh5nUOSfnguC2J9IPSyOSX5yXSyZU+9dBm9nvcXGB555T84uiaB2NEJowu/Bs06YN169fz7Y9MTGRNm3aWCQoYaOioqBWLahWTetIrEbvCserI4WnRiS/OK+KZ8AzzYGb2e/p1g1On4YjR7SORAhNmF14KoqCLoe+b9euXaNw4cIWCUrYoIwM+P13hxzN/qCYmsCxY3D8uNahOB3JL86rZgxc94OrpbWOxMratFEnxpc5g4WTMrmPZ7e7zas6nY6XXnrJqCksKyuLv//+m+bNm1s+QmEbtm1Tp1JygsLzVGXA21t96vnee1qH4xQkvzi5rCyqx8Df9VBHfzsyd3f+rpBIqZmRlPn4Y62jEaLAmVx4+t5dukxRFIoWLYqXl5fhPXd3d5o2bcrgwYMtH6GwDcuWQYUK8NhjWkdidZluwNNPS+FZgCS/OLnduymS4uD9O+8TUwPqHUFtcq9USetwhChQJhee8+fPByA4OJi33npLmr2ciV6vzj3Xo4fzTDHUrRv06wcXL0L58lpH4/Akvzi5ZctILgIXnORX7WQVyCgEbsuWwZtvah2OEAUqX6Pa5abgZPbuhdhYo2Z2R1rKLkfPPANubtIPq4BJfnFCigK//qoO6nuktfTsR4bH3S49S5dqHYoQBc6kJ56PPfYYGzduxM/PjwYNGuTY+f+eAwcOWCw4YSMWL4YyZaBlS60jKTjFikG7drBkCYwYoXU0Dk3yi5Pbtw/OnuXo41oHUrCO1YYaUbvh/Hm1G5MQTsKkwvO5554zdPbv4oBLJYo86PVq8dWjB7i6ah1NgYmM1FHPDbpuR33a64grNdkIyS9ObvFiKFWKc0FXtY6kQB2vBnh4qE89w8O1DkeIAmNS4RkREZHj34UT2L1b7efYq5fWkRS449Uh0xUKLV0KI0dqHY7DkvzixBRFLTy7d0dxnaF1NAUq3RPo0EH9/qXwFE7ESXrUiHxbvBjKljUs7+bwfTvvk+Z1tx/W4sVahyKEY9q7V21qdsIPtoD6fe/ZA2fPah2JEAXGpCeefn5+efa7ul9Oq44IO3Wvmb1nT3WpNyd0tDZUX7YTLlyAwECtw3FIkl+c2OLF6trlrVrBVq2D0UCnTuDpqTa3v/WW1tEIUSBMKjwnTZpk5TCETdq5Ey5dYt6tKVyInEJEhPOtLXy8Ov/1wxo1SutwHJLkFyflpP3HjRQtCh07qgW4FJ7CSZhUeA4YMMDacQhbtHgxlCvHhfKxWkeimXRP1MnkFy+WwtNKJL84qT171JYEZ21mv6dXL3j+eThzBipWNHRlcsYP+sI5mFR4JiUl4ePjY/h7Xu7tJ+xcVpb6lO/558HlG62j0dSvrivovhs4dw6CgrQOx+FIfnFOu8ObU7sIFHWmadoeEBmpwy0N3vfyUp/+vvOO1iEJYXUmddzz8/Pj6lV1qotixYrh5+eX7eveduEgtm+Hy5eNnkY4y6CiB/17d3R7bpM9O+vPxVIkvzghvZ5aR9W5LJ21//g9GR6oC1YsXiy5RDgFk554btq0ieLFiwOwefNmqwYkbMTixeqkxqGhsFbrYLSV7gEnqkLNRYtkeTsrkPzihHbuxCdZHbwXqnUstqBXL+jVC78WcKO41sEIYV0mFZ6tW7fO8e/CQWVkqIXnSy85z9rsD3G0DtRc+hecOAFVq2odjkOR/OKEfv6ZRB/nWZv9nlyfaD7zDBQpQt0jt9gmvwLCwZlUeD7oxo0bzJ07l+joaABq1arFwIEDDU8thJ1bswYSEqB/f60jsRnHqwM+PvDjjzB2rNbhODTJLw4uLQ0WLeLveshM0ndFflGY5ypDvcOw7XFAPu8LB2b2r/22bdsIDg5mypQp3Lhxgxs3bjBlyhQqVqzItm3brBGjKGg//AD160PdulpHYjMy3VDnM/3xR3UaGGEVkl+cwB9/wI0b/F1f60Bsy9/1ocR1KHdR60iEsC6zC8+hQ4fSu3dvzpw5Q1RUFFFRUZw+fZrnn3+eoUOH5iuI6dOnExwcjKenJ6GhoezduzfXfWfPnk2rVq0Mgw7CwsLy3F+Y6cYNMpctZW3pw1pHYnPmZ81VVxjZsUPrUByW5Bcn8MMP0KgRCSXVl860GlpezgZDog/Ul9QrHJzZhefJkyd58803cb1vwl9XV1fCw8M5efKk2QEsWrSI8PBwIiIiOHDgAPXr16d9+/aGUa4P2rJlC3369GHz5s3s2rWLwMBA2rVrR2ys8841aVFLluCihyPysDOb8xXgRjHUG6ewCskvDi4hAVatYnXJfVpHYnMUFzhSD+r8A66ZWkcjhPWYXXg+9thjhr5X94uOjqZ+ffPbTr7++msGDx7MwIEDqVWrFjNnzsTb25t58+bluP/PP//M66+/TkhICDVq1GDOnDno9Xo2btxo9rVFDn74gVOVIaWo1oHYoLs3BpYsgdu3tY7GIUl+cXCLFoGi8E8drQOxTYfrgdcdqHpC60iEsB6TBhf9/fffhr+PGDGCkSNHcvLkSZo2bQrA7t27mT59OhMmTDDr4unp6ezfv5/Ro0cbtrm4uBAWFsauXbtMOkdqaioZGRm5DjxIS0sjLS3N8PphE1Q7tVOnYMcO/u6udSC263A9eHxbIkv6e3PsvpunrDaSf/acX0ByjFl++AE6dCC18O9aR2KTEkrBpbLqICMhHJVJhWdISAg6nQ5F+e+m+k4OKyz07duX3r17m3zxhIQEsrKyKF26tNH20qVLExMTY9I53n33XQICAggLC8vx/fHjxxMZGWlyTE7tp5+gaFFiqidrHYnNuu6vTgFT/zBGhafIP3vOLyA5xmTHj8PevWqLwVEpPHNzuD60W4faLcHfX+twhLA4kwrPM2fOWDuOfJkwYQILFy5ky5YteHp65rjP6NGjCQ8PN7xOSkoiMDCwoEK0H4qijtju0YNM9/laR2PT/q4HHVZD4VuQUkTraOyfPecXkBxjsh9/BF9f6NQJjmodjO36pw60X4vaLSGfA+qEsGUmFZ5BVlqf2t/fH1dXV+Li4oy2x8XFUaZMmTyP/fLLL5kwYQIbNmygXr16ue7n4eGBh4eHReJ1aDt3qk3tc+bAVik883K0Djy9Buoegd3NtI7G/tlzfgHJMSbR69UWld69IY8iXkBqEXWltOrffy+Fp3BI+ZpAHuDYsWOcP3+e9PR0o+3PPvusyedwd3enYcOGbNy4kS5dugAYOvIPGzYs1+MmTpzIuHHjWLt2LY0aNcpX/OIB//sfVKoEjz8OW9VNMsVJzm57Q0wNeGw/7G6KTPZsBZJfHMzatXDuHAwcqHUkduFQCFRf/BccOgQhIRpHI4RlmV14nj59mq5du3LkyBGjflm6u0srZmVlmXW+8PBwBgwYQKNGjWjSpAmTJk0iJSWFgXcTVP/+/SlXrhzjx48H4PPPP+ejjz5iwYIFBAcHc+XKFQCKFClCkSLS7pkv166pS2SOHQsuspSIKfY1ggE/QNBZOFdR62gch+QXx3T8jY74lIFZa5rBWq2jsX3HqwMBATBjhvpQQAgHYnaVMXLkSCpWrMjVq1fx9vbm6NGjbNu2jUaNGrFlyxazA+jduzdffvklH330ESEhIRw6dIg1a9YYBgScP3+ey5cvG/afMWMG6enp9OjRg7Jlyxq+vvzyS7OvLe6aP1/t4ylPI0x2tiIklIBGD0xHKJNhPxrJLw7o/HmqnlA/rEnrgGkUV9hS9RLp82eBzJIgHIxOuX8oqQn8/f3ZtGkT9erVw9fXl71791K9enU2bdrEm2++ycGDB60Vq0UkJSXh6+tLYmIiPj4+WoejPb0eqlWDpk3VPlhIE7upQnfBU+vhm1HZ5z111mmVHvX3y97zC0iOyWbMGNK+GMdXb0KGdIU1WdEkeOMbcJk6Tfp63iW/W47B7CeeWVlZFC2q3mX9/f25dOkSoA4QOH78uGWjE9a3YYM6qGjIEK0jsTuHQ0DvAg1svxayG5JfHEx6OsyZw+H6UnSaK9nnbpP7zJlqi5QQDsLswrNOnTocPqzObhsaGsrEiRPZsWMHY8eOpVKlShYPUFjZjBlQty40b651JHbnjpc69UnD/aDTax2NY5D84mCWL4e4OPY11joQ+7SvMfDPP7Bjh9ahCGExZheeY8aMQa9X77Jjx47lzJkztGrVilWrVjFlyhSLByis6OJF+O039WmnTprX82NfYyiWCFUeWOJO+nrmj+QXBzNjBrRqRXwprQOxT6crAlWqqD9HIRyE2aPa27dvb/h7lSpViImJ4fr16/j5+RlGngo7MXs2aYX0fH3hdUYjTe35camcusRdo31worrW0dg/yS8OJDoatmyBn3+GE39qHY19cgFeew3efx8mTYKSJbWOSIhH9khz51y4cIELFy5QvHhxuSnYm/R0mD2bI/UgXfpePZJ9jaHqCfC7rnUkjkXyi5379ltSvOHT6Be0jsS+vfSSOs3dnDlaRyKERZhdeGZmZvLhhx/i6+tLcHAwwcHB+Pr6MmbMGDIyMqwRo7CGX36By5fZE6p1IPbvSF1I9Yamu7SOxP5JfnEQ16/DvHnsawxZ+V6mRABQogS8+CJMmQJpaVpHI8QjM7vwHD58OLNmzWLixIkcPHiQgwcPMnHiRObOncuIESOsEaOwNEWBL7+ETp1IuNtyI/0R8y/TDf5qoo5u90o1fk/6eppH8ouDmDED9Hr2yqAiywgPhytXYMECrSMR4pGZ/Vl0wYIFLFy4kA4dOhi21atXj8DAQPr06cMM6QRt+9auVUdKTpsGW1ZqHY1D+KsxtNgOjf+Cba21jsZ+SX6xf5+O0TFmTmkYMIDUIrLqjkXUqAGdO6sPDAYMkBXmhF0z+3+vh4cHwcHB2bZXrFgRd3d3S8QkrO3zz6FxY3Vd9vvI07n8Sy2srq/cZA+4pT90d5ELyS/2r/5hUOLimOomRadFvf02HDsGf/yhdSRCPBKzC89hw4bxySefkHZfX5O0tDTGjRvHsGHDLBqcsII//1RHmr7/vkyhZGE7W4DXbXVeT5E/kl/sXEYGLbfD0dpw3V/rYBxMy5bQsiWxrz0rE8oLu2ZSU3u3bt2MXm/YsIHy5ctTv359AA4fPkx6ejpt27a1fITCsj75BOrVg2ef1ToSh3PTDw7Xh+Y71HWpM93+ey8yUue0y2g+jOQXB/LTT/jdhIV9tA7EcdxrhYqIUOCjjyjXrp3aXerppzWOTIj8Manw9PX1NXrdvXt3o9eBgYGWi0hYz65dsH49LFkifYSs5M9WalPjYwdgr8wYYBLJLw4iMxPGjeNYTbhaWutgHFRYGBfKQ+DYsdC+vbRaCbtkUuE5f/58a8chrE1RYMwYqFMHHnjCJCznRgn4ux60/BMONoCMHLolGj3BEJJfHMV338GpU2x7VetAHNO9vFG5NfT7eResXg0dO2oclRDmy/djr/j4eLZv38727duJj4+3ZEzCGtavh02bYNw4edppZVufAO9UCN2tdST2S/KLnbl9Gz7+GJ5/nriyWgfj2E5VQR0YOno03F1eVgh7YvZ0SikpKQwfPpwffvjBsKayq6sr/fv3Z+rUqXh7e1s8SPGI9Hp47z1o0YLIA8/BQa0Dcmw3/dTVjFrsgP2N4PbdXwmZMeDhJL/YqalTybocy/SSC7WOxPHpYE7VbbwyF3Vez379tI5ICLOY/egrPDycrVu38vvvv3Pz5k1u3rzJihX/396dx0VVrw8c/wyIA6Ggsrog4pJimnYNEfWnmCiuqZWZtxTN2y8NSaVuLtk1SqMsCy9uaWr5K9Nyw8glpERLxKtm5QI3TaKQzQ0MRGTm/P44OTkyEBjMGeR5v17z0jnrw2F45jnnfM/3G0dSUhLPPfdcTcQo/qr16+Hbb+H110FqH6vY1wd0CvRJ0jqS2kXySy10/jxER3O0m9rURNS8TB841QG1+dTVq1qHI0SVVLnw3Lx5M6tXr2bw4MG4uLjg4uLCkCFDWLVqFZs2baqJGMVfceUKvPACPPKI2h2HsIoiZ9jfB7ofAvdcraOpPSS/1EJz54KisDdY60Dqlj0DgHPnYOFCrUMRokqqXHgWFRXh5VX2kUVPT0+KioosrCE09dprcPmyOuKFsKqDPeBSYxi8E7DwHJF02F+W5Jda5ttvYeVKeOUVipy1DqZuuegGX3e/zvX5L8PPP2sdjhCVVuXCMygoiHnz5lFcXGyadvXqVaKioggKCqrW4MRflJoKb7+ttu/09dU6mjrHUA92DYLWZ6HjCa2jqR0kv9QiBgOEh5PrrvDK+WlaR1Mn7esDV52AadNAUeRkVtQKVX64KCYmhkGDBpXp4NnR0ZHdu3dXe4DiNhmNMGkStGqlDrUmNHH6bjjpr171/Kk1FMuzMRWS/FKLLFsGycnETwTFXutg6qbrevXk9tFP42DzZq3DEaJSdIpS9bG3ioqK+Oijj0hNTQXA39+fxx9/HCcnp2oPsLoVFBTg6upKfn4+Li4uWodTc5YsgYgISEqCPn3kLFhDDa7AM0shrT3EjSo7/07qz7M6/r5qc36BOpJj0tPVPoHDwojyXKZ1NHWbAvO+GwnJybwRlmM6ub2T8soNdeJvqw6o0hXP69ev06FDB+Lj43nqqadqKibxV/33vzBzJkyZovb3JjT1W0P4YiCM2A6p/pDWwXz+zScFd+KXRWVJfqklDAaYOJF8+0KWuUjRqTkdsHQpdOzI0M9h8yNI7yXCplWpjaeDg4NZ2ythg0pKYOxYaN5cnna0Icfug9T28GAcNCzQOhrbJPmllli4EJKS2DoKShy1DkYA0KwZrFhBpxPqkL1C2LIqP1wUHh7OG2+8QWlpaU3EI/6q2bPhhx/g44+hQQOtoxE36GD7g1BaD0ZtAZ1B64Bsk+QXG5ecDC+9BLNm8bOf1sGIG6KidESljeXbrjDkc3A7r3VEQpSvyg8X/ec//yExMZEvvviCzp074+xs3ofGli1bqi04UUUbNqhPsS9aBN26aR2NuMVVZ9jyEIxfByGJkDBQ64hsj+QXG5aVxZWBPWnYPQiiouC1aK0jErfYNRha/ApjNgALCkDaQQobVOXCs1GjRjz88MM1EYv4K44dgyefhL//HWbM0DoaUY6f/dT2noN2Q5Y3HL9X64hsi+QXG3XtmjoIBahPTzs4aBuPsKhEDxsfg3+sAsaNg61bwU69sXmjLXldbkcubEOVC8+1a9fWRBzir0hPhyFDwN8fVq0CnXnLcnmi3bak9ADvbBgRpz54lC63LE0kv9ggg0EtYo4cYeMT8I+mTbWOSFTggjtseRj+viEepk+HxYvLfCcIoaVKt/E0Go288cYb9OrVi4CAAGbNmsVVGSNWe3l5EBoKTk6wYwfcJR1F2jwdfDYcfvZVb4l5Z2kdkPYkv9goRVE7J9+8GTZsINNHRtyqDX68G7Wf1dhYdfQ6IWxIpa94LliwgJdffpmQkBCcnJxYvHgxubm5rFmzpibjExXJyYH+/SE/H77+Gn4falC+FGyfsR58MgbGfwDj1sH/jYPsZmVvh9WV22OSX2yQopASZEdginqidPQ7C53QCpsVlT2ZPsHQb+5c9nw1F/5H64iEUFX6iue6detYtmwZu3fvZtu2bXz22Wd89NFHGI3GmoxPlCcjA4KD4eJF2LsX2rbVOiJRRSV6teC82ER94KhFhtYRaUfyi40pLYUpUwhMgfihcFSeVayV9vWFvX3VhxmDvwTu7PNXUUtUuvDMyMhgyJAhpvchISHodDrOnTtXI4GJChw5AoGBUFysFp0dOvzpKsI2XXOCD8dBjheEfQD3HFen17XbmZJfbMiVK/y3owPGle8S9yAcCdA6IHHbdJDUD/b0h777YORW1AfFhNBQpQvP0tJSHB3Newt2cHDg+vXr1R6UKIeiwHvvQe/e4OsLBw/C3XdrHZX4i645qsXnyY7wyCYY8AXY1bF+PiW/1Iwqn8CcOAGBgfj+DB89Dsf+VnOxCev55n9g88Nwzwn4tY0jMTPqzkmtsD2VbuOpKAoTJkxAr9ebphUXFzN58mSzvvakn70acv68Ovb6hg3wv/8LMTHqA0XijmCoB1sfgqymELIHfDIgbgRc8NA6MuuQ/KIxoxGWL4d//hPatGHVU3Xns1dXHO8MFxvD6E/h6RVAwHp1lDt54l1YWaULz7CwsDLTnnjiiWoNRlhgNMKHH8Lzz6vtrtb/nizEnUcHB3vCLz4waitMXgH7+kByT60Dq3mSXzT0/fcwdSrs3w+TJ8OiRVx40/nP1xO1zrkW8O5kGBoPnR5/nLSox9k9CJ5dLI0/hfVUuvCU/vWsTFFg9251eLrDh+HRR9X+2Ly9tY5M1LBMH1gxBYL3Qt8k6HYE4n7Q8X0XMNrfmU+8S36xHtPnZvxPsGABxjWrudgE4sNgwvLlGkcnalqxE2weDZ3mb8U7bBThSwC7GfDCCyB9tAorqPJY7aKGFRfDunUQEACDB/NL9mHWTAQ2brRYdN7cfquuPZByJyt1gD0DYFk4ZDaHEdsh4t/Qez+Qmal1eMKGVOlvXlFo8QuM2gy0aweffcbuUFj+DDL2eh0T9d0olkyFpL7A2rXg5wf/+AckJ6sXPoSoITZReC5dupRWrVrh6OhIYGAghw4dqnD5Tz/9lA4dOuDo6Ejnzp3ZsWOHlSKtIVevQlwcjB+vFpdhYZy5dIR142DNJPjFt+KiUgrOO9dFN/h0DCyfonY43ycJaNkShg7lnuPgWPTH718+B5bV+fyiAEePwosvQocOTFqtjufNW2/B2bMc6qFeSQfJJXVNaX3Y3xd19LuXXoKEBOjZkzxPOxIG6FgSoSPqZflMiOpV5SEzq9vGjRuJjIxkxYoVBAYGEhMTQ2hoKGlpaXh6epZZ/sCBA4wdO5bo6GiGDRvG+vXrGTlyJEePHqVTp04a/AS3obAQDh2CAwfU17598NtvcM89MG0asRdf4aK75VUlAdRNuV6w7SHYOQRm+S2H1at55JBaU2Q2h5/aqG1Ds6Qlhpk6mV+MRjh5Er75hhFboVU6ENWNIidIaw8nnoAzrYH8GfDmDI2DFTahUSN48UWiSubidxbu+xb6fQUD9sBlVzW/0G493Hcfr2zoiGJ3ZzTxEdrQKYq219QDAwMJCAhgyZIlgDp0no+PDxEREcyaNavM8mPGjKGwsJD4+HjTtB49etC1a1dWrFjxp/srKCjA1dWV/Px8XFxcqu8HuVlJCVy4AOfOwa+/QmYmBz4Mx/083F3qp55dKgq4ukJQEPTpAw89BO3bA1JciordSPjvROpofQbanIHWP8FdN0aY9PbmtHM2l5pAwKNvqrfQWrQAd3f15eJS6SdZq9qO1Cp/X1Vg7fwCVjgGBoM6Wll2NmRmsm3pQEbe9zKkpqqvtDT1Loq9Pec8DfzSEtLuhp9b/XFlU4g/41ACvulqfmlzBjzOq9NLHNR+h30eGK/mFj8/tXs/T081v7i5gX3NfNBsLb+I26PpFc+SkhKOHDnC7NmzTdPs7OwICQkhOTnZ4jrJyclERkaaTQsNDWXbtm0Wl7927RrXbuowNz8/H4CCV14Bvb5sW5Zb3xsMaoe7167B9evqvyUlZtN+SfsSx2LwqOelfiEUF5tvw96eZg3VEWoKhg5VC8xu3YiO6wm6XVCwC96fU8GREuIPs2f/XjTqIbcjHOyovm10CTxzwTsnG49caHwW8ub8E/0tXWEadGDv6UVeaQ4l9aF5615q11y/v75N+wRDPbi/RzhB/wHFDnbu1qHYQe8+s9UvlZtfdnamQrbg98++xuezgHXyC1Qxx1jKNzfyyS15xTTt6lUoKODyzz/gWAz6Erj5tOEB4FzCy1xsAhfc4EJPyPWAc80MXK9/04LXf38JUQnFwPGW6ot+4HgVvHLBMwe8cuDKnnU0ugwNisquW+QIV53ArWUXcHLip+yDtO44xCzP4OBgOZdYem+ntgq0pfwibp+mhef58+cxGAx4/T7G+A1eXl6kpqZaXCc7O9vi8tnZ2RaXj46OJioqqsx0n0WLbjPqiuRYnmwwwGXU1+9XXoTQjALk3PRZzfrG8nJHlpad9k10pXZx5coVXF1dqx5bNbJGfgFr55hyFP3++tV6uxSiXMW/vy5998e09OprK20L+UXcPs3beNa02bNnm13BMBqNXLx4ETc3N3Q20HFuQUEBPj4+/PLLL3LrwAI5PhWzteOjKApXrlyhWbNmWodiNZJjai85NhWzteNTF/PLnUjTwtPd3R17e3tycsyvFObk5OBdTn+V3t7eVVper9ebjYYC0KhRo9sPuoa4uLjYxB+2rZLjUzFbOj62ciXCGvkFJMfcCeTYVMyWjo+t5Bdx+zTtTql+/fp069aNxMRE0zSj0UhiYiJBQUEW1wkKCjJbHiAhIaHc5YUQdZPkFyGEsD2a32qPjIwkLCyM+++/n+7duxMTE0NhYSETJ04EYPz48TRv3pzoaLVt2bRp0+jbty+LFi1i6NChbNiwgcOHD7Ny5UotfwwhhA2S/CKEELZF88JzzJgx5OXl8a9//Yvs7Gy6du3Krl27TA38MzIysLP748Jsz549Wb9+PXPnzmXOnDm0a9eObdu21Z4+9m6h1+uZN29emVt1QiXHp2JyfCpW1/MLyGekInJsKibHR9QEzfvxFEIIIYQQdYNNDJkphBBCCCHufFJ4CiGEEEIIq5DCUwghhBBCWIUUnkIIIYQQwiqk8BRCCCGEEFYhhacNSU9PZ9KkSfj5+eHk5ESbNm2YN28eJSUlWoemmaVLl9KqVSscHR0JDAzk0KFDWodkE6KjowkICKBhw4Z4enoycuRI0tLStA5L2DDJL2VJfrFM8ouoSVJ42pDU1FSMRiPvvvsuJ06c4J133mHFihXMmTNH69A0sXHjRiIjI5k3bx5Hjx6lS5cuhIaGkpubq3VomktKSiI8PJyDBw+SkJDA9evXGThwIIWFhVqHJmyU5Bdzkl/KJ/lF1CTpx9PGvfnmmyxfvpyffvpJ61CsLjAwkICAAJYsWQKowx36+PgQERHBrFmzNI7OtuTl5eHp6UlSUhJ9+vTROhxRS0h+kfxSGZJfRHWSK542Lj8/nyZNmmgdhtWVlJRw5MgRQkJCTNPs7OwICQkhOTlZw8hsU35+PkCd/KyI2yf5RfJLZUh+EdVJCk8bdvr0aWJjY3n66ae1DsXqzp8/j8FgMA1teIOXlxfZ2dkaRWWbjEYj06dPp1evXrV6aEdhXZJfJL9UhuQXUd2k8LSCWbNmodPpKnylpqaarZOZmcmgQYMYPXo0Tz31lEaRi9ogPDyc48ePs2HDBq1DERqQ/CJqkuQXUd3qaR1AXfDcc88xYcKECpdp3bq16f/nzp2jX79+9OzZk5UrV9ZwdLbJ3d0de3t7cnJyzKbn5OTg7e2tUVS2Z+rUqcTHx7Nv3z5atGihdThCA5Jfqk7yS+VIfhE1QQpPK/Dw8MDDw6NSy2ZmZtKvXz+6devG2rVrsbOrmxel69evT7du3UhMTGTkyJGAessnMTGRqVOnahucDVAUhYiICLZu3crevXvx8/PTOiShEckvVSf5pWKSX0RNksLThmRmZhIcHIyvry9vvfUWeXl5pnl18Sw8MjKSsLAw7r//frp3705MTAyFhYVMnDhR69A0Fx4ezvr164mLi6Nhw4amdmmurq44OTlpHJ2wRZJfzEl+KZ/kF1GTpDslG/L++++Xm/Tq6q9pyZIlvPnmm2RnZ9O1a1f+/e9/ExgYqHVYmtPpdBanr1279k9vu4q6SfJLWZJfLJP8ImqSFJ5CCCGEEMIq6mYDHyGEEEIIYXVSeAohhBBCCKuQwlMIIYQQQliFFJ5CCCGEEMIqpPAUQgghhBBWIYWnEEIIIYSwCik8hRBCCCGEVUjhKYQQQgghrEIKTyGEEEIIYRUyVrsQVnD+/Hk8PDxqbPsyAJkQdZvkGFFbSOEphBXs2LGDgoICGjZsqHUoQog7kOQYUVvIrXYhrEC+EIQQNUlyjKgtpPC0ccHBwUyfPt1mt3enuHDhAp6enqSnp1f7tktKStDr9dW+3Rsee+wxFi1aVGPbF3c2yTHWITlGCJUUnhqZMGECOp0OnU6Hg4MDXl5eDBgwgDVr1mA0Gk3LbdmyhVdffVXDSGuvqnwBLliwgBEjRtCqVSvTtH379jF8+HCaNWuGTqdj27Zt5a4/ceJE5s6da3He/v376d27dxUir5q5c+eyYMEC8vPza2wfovaRHFPzJMcIUXVSeGpo0KBBZGVlkZ6ezs6dO+nXrx/Tpk1j2LBhlJaWAtCkSRO5fVLDioqKWL16NZMmTTKbXlhYSJcuXVi6dGmF6xsMBuLj43nwwQctzk9NTcXf37/a4r1Vp06daNOmDR9++GGN7UPUTpJjbIPkGCH+IIWnhvR6Pd7e3jRv3py//e1vzJkzh7i4OHbu3Mn7778PlD2j3rRpE507d8bJyQk3NzdCQkIoLCw0LTt16lSmTp2Kq6sr7u7uvPTSSxU+jbhr1y569+5No0aNcHNzY9iwYZw5c8Y032g0snDhQtq2bYter6dly5YsWLDAbH50dDR+fn44OTnRpUsXNm3aZLaP4OBgIiIimD59Oo0bN8bLy4tVq1ZRWFjIxIkTadiwIW3btmXnzp1V3u6zzz7LCy+8QJMmTfD29ubll18G1Ks9SUlJLF682HTVp7xbXDt27ECv19OjRw+z6YMHD2b+/PmMGjWq3OMHcODAARwcHAgICKhwuVtNmTKl3KsULVq04PXXX6/0toYPH86GDRuqtH9x55McIzlGcoywNVJ42pgHHniALl26sGXLljLzsrKyGDt2LE8++SSnTp1i7969PPTQQ2ZJ/4MPPqBevXocOnSIxYsX8/bbb/Pee++Vu7/CwkIiIyM5fPgwiYmJ2NnZMWrUKNOtuNmzZ/P666/z0ksvcfLkSdavX4+Xl5dp/ejoaNatW8eKFSs4ceIEM2bM4IknniApKclsPx988AHu7u4cOnSIiIgIpkyZwujRo+nZsydHjx5l4MCBjBs3jqKioipv19nZmZSUFBYuXMgrr7xCQkICixcvJigoiKeeeoqsrCyysrLw8fGxeAz2799Pt27d/uQ3U77t27czfPhwdDodgOlnADh79iytW7cus86JEydYuXIlCxcutLhNf39/jh07VukYunfvzqFDh7h27VrVghd1juQYyTEgOUZoSBGaCAsLU0aMGGFx3pgxYxR/f39FURSlb9++yrRp0xRFUZQjR44ogJKenm5xvb59+yr+/v6K0Wg0TZs5c6ZpW7duz5K8vDwFUH744QeloKBA0ev1yqpVqywuW1xcrNx1113KgQMHzKZPmjRJGTt2rNk+e/fubXpfWlqqODs7K+PGjTNNy8rKUgAlOTn5trerKIoSEBCgzJw5s1I/6w0jRoxQnnzyyQqXAZStW7danNeuXTslPj5eKSkpUZ599lll/vz5pnmrV69WioqKyqwTFhamBAYGlru/Rx99VOnbt++fxn7Dd999V+FnQ9Q9kmMkx0iOEbZIrnjaIEVRTGe2N+vSpQv9+/enc+fOjB49mlWrVnHp0iWzZXr06GG2blBQED/++CMGg8Hivn788UfGjh1L69atcXFxMTV8z8jI4NSpU1y7do3+/ftbXPf06dMUFRUxYMAAGjRoYHqtW7fO7FYawL333mv6v729PW5ubnTu3Nk07cYVjtzc3NveLkDTpk3Jzc21GG95rl69iqOjY5XWueHUqVOcO3eO/v374+DgwKuvvkpCQoJpflFREU5OTmbrlJaWsmXLFh5++GHTtKeffprVq1eb3l+5cqXMehW5sezNV0KEKI/kGMkxkmOEVqQDeRt06tQp/Pz8yky3t7cnISGBAwcO8MUXXxAbG8uLL75ISkqKxeUrY/jw4fj6+rJq1SqaNWuG0WikU6dOlJSU/GlS+u233wD4/PPPad68udm8W7v2cHBwMHt/40nbm9+D2u7qr2735id2K8Pd3b3Ml2tlbd++nQEDBpi+VFxcXHB0dOT8+fM4Oztz1113lVnnzJkzXLlyxfSlaDQa+fTTT82+fL///nvGjBlT6TguXrwIUKMjl4g7h+QYyTGSY4RW5Iqnjfnyyy/54YcfzM5Ub6bT6ejVqxdRUVF8++231K9fn61bt5rmp6SkmC1/8OBB2rVrh729fZltXbhwgbS0NObOnUv//v3x9/c3S47t2rXDycmJxMREi7F07NgRvV5PRkYGbdu2NXuV19apMqpru/Xr1y/3KszN7rvvPk6ePHlbscbFxTFixAizaYMHD2bHjh3s3buX4ODgMutcvnwZgAYNGgCwe/duLl26ZPpiOXjwIJmZmaYHDi5cuEDbtm0B9UpG165dyzzMcfz4cVq0aIG7u/tt/Ryi7pAcIzlGcozQklzx1NC1a9fIzs7GYDCQk5PDrl27iI6OZtiwYYwfP77M8ikpKSQmJjJw4EA8PT1JSUkhLy/PrBuNjIwMIiMjefrppzl69CixsbHldvzbuHFj3NzcWLlyJU2bNiUjI4NZs2aZ5js6OjJz5kxeeOEF6tevT69evcjLy+PEiRNMmjSJhg0b8vzzzzNjxgyMRiO9e/cmPz+fb775BhcXF8LCwm7ruFTXdlu1akVKSgrp6ek0aNCAJk2aYGdX9lwrNDSU2bNnc+nSJRo3bmya/ttvv3H69GnT+7Nnz3Ls2DGaNGlCy5Ytyc3N5fDhw2zfvt1se0OHDuXFF1+kT58+DB48uMz+fH190el0fPzxxzg7O/P8888zdOhQ4uLi8PHxYfLkyYSEhJieRnVzczNdbahXrx6enp78+uuvZl+Q+/fvZ+DAgZU6LqLukBxjmeQYyTFCQ9o2Ma27wsLCFEABlHr16ikeHh5KSEiIsmbNGsVgMJiWu7nx+smTJ5XQ0FDFw8ND0ev1yt13363ExsaaLfvMM88okydPVlxcXJTGjRsrc+bMMXsQ4NbG8AkJCYq/v7+i1+uVe++9V9m7d69ZI3eDwaDMnz9f8fX1VRwcHJSWLVsqr732mml9o9GoxMTEKO3bt1ccHBwUDw8PJTQ0VElKSip3n4qiKL6+vso777xjNu3m/d7udkeMGKGEhYUpiqIoaWlpSo8ePRQnJycFUM6ePVvu76N79+7KihUrzKZ99dVXpt/Rza8b23/vvfeUXr16WdxeYGCg2e/mVq+99pri4uKieHl5KatXr1aOHTum+Pr6Ks7Ozspjjz2mXLx40Wz5li1bKpcvX1YMBoPSvn17paCgwDTv6tWriqurq5KcnFzu/kTdIzlGcozkGGGLpPC8g1T2CUtRVnx8vOLv72/2hfxnhg8frrzxxhsW50VGRip79uyprvCUHj16KKdOnVJiY2OV8PBws3nLli1TBgwYUG37EqI8kmNun+QYIVRyq10I1FtXP/74I5mZmZVu49W7d2/Gjh1rcV5ERATNmjWrtviaNm3KJ598wp49e9i1a5fZPAcHB2JjY6ttX0KI6ic5RgiVTlEqGHJC1CrBwcF07dqVmJgYrUMR1Sw8PJyvv/6axMREadwvNCM55s4lOUZYixSeQgghhBDCKqQ7JSGEEEIIYRVSeAohhBBCCKuQwlMIIYQQQliFFJ5CCCGEEMIqpPAUQgghhBBWIYWnEEIIIYSwCik8hRBCCCGEVUjhKYQQQgghrEIKTyGEEEIIYRVSeAohhBBCCKuQwlMIIYQQQljF/wPMIf1q1Q7PYgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "hist_data = np.genfromtxt('co2-100K.prob-density.dat')\n",
    "fig, axs = plt.subplots(2,2)\n",
    "fig.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.3, hspace=0.6)\n",
    "plt.setp(axs,xlim=(-3,3),ylim=(0,0.6))\n",
    "mode_list = [6,7,8,9]\n",
    "imode=0\n",
    "bins=np.linspace(-5,5,201) # Histogram bins\n",
    "#Harmonic density is a Gaussian\n",
    "harmonic_density = [1/np.sqrt(np.pi)*np.exp(-1.0*x**2) for x in bins]\n",
    "for plot_row in range(2):\n",
    "    for plot_col in range(2):\n",
    "        axs[plot_row,plot_col].set_xlabel(\"Displacement $(1/\\sqrt{\\omega_{\\\\nu}})$\")\n",
    "        axs[plot_row,plot_col].set_ylabel(\"Probability density\")\n",
    "        axs[plot_row,plot_col].set_title(\"Mode # \" + str(mode_list[imode]))\n",
    "        axs[plot_row,plot_col].hist(hist_data[:,0], weights=hist_data[:,imode+1], bins=bins, density=True,\\\n",
    "                      histtype= 'stepfilled', color = 'olive', label = 'PIGLET')\n",
    "        axs[plot_row,plot_col].plot(bins,harmonic_density, linewidth = 1.0, color = 'red', label = \"Harmonic\")\n",
    "        if imode==1:\n",
    "           axs[plot_row,plot_col].legend(loc='best',bbox_to_anchor=(1,1),shadow=True, fontsize='small')\n",
    "        imode += 1\n",
    "plt.show()              "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52657aa1",
   "metadata": {},
   "source": [
    "1.5.2. Calculation of Anharmonic Measure\n",
    "===\n",
    "Now we will calculate the anharmonic measure as proposed by Carbogno and co-workers, see https://journals.aps.org/prmaterials/abstract/10.1103/PhysRevMaterials.4.083809. For that purpose we need to import anharm module available within pyEPFD. We need a trajectory including the forces. Files co2-1.r and co2-2.r contains forces, Kohn-sham eigenvalues etc. for the PIGLET trajectory sampled every 5-th step after discarding 2000 steps for equilibration.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "69d50c45",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time spent on qbox class: 0.16231036186218262 s.\n",
      "Time spent on qbox class: 0.1451880931854248 s.\n",
      "Process-id0: Time spent on read_pyepfd_info class: 0.00029730796813964844 s.\n",
      "anharm_measure: Execution time (s): 0.32349658012390137\n"
     ]
    }
   ],
   "source": [
    "from pyepfd.anharm import anharm_measure\n",
    "#We need qbox class to read qbox outputs \n",
    "from pyepfd.coord_util import qbox\n",
    "#We need pyepfd_io module to parse and store restart file\n",
    "from pyepfd.pyepfd_io import *\n",
    "\n",
    "### qbout a list of objects each of which stores qbox output information\n",
    "qbout = [] # initialize it with a empty list\n",
    "##\n",
    "\n",
    "# Initiating the qbout objects from co2-1.r and co2-2.r\n",
    "for i in range(2):\n",
    "    qbout.append( qbox(file_path='co2-'+str(i+1)+'.r', io = 'r') )\n",
    "\n",
    "## Concatenating coordinates and forces into a single array\n",
    "for i in range(len(qbout)):\n",
    "    if i == 0:\n",
    "       coords = qbout[i].coords\n",
    "       forces = qbout[i].forces\n",
    "    else:\n",
    "       coords = np.vstack((coords, qbout[i].coords))\n",
    "       forces = np.vstack((forces, qbout[i].forces))\n",
    "\n",
    "### Now we read the normal mode information from the restart file\n",
    "phonon_info = read_pyepfd_info(file_path='../2_normal_mode_phonon/enmfdphonon.xml')\n",
    "\n",
    "## Now we calculate anharm measure\n",
    "anharm = anharm_measure(dynmat = phonon_info.ref_dynmatrix, \\\n",
    "                        mass = phonon_info.mass, \\\n",
    "                        forces = forces, \\\n",
    "                        disp_coords = coords, \\\n",
    "                        opt_coord = phonon_info.coord, \\\n",
    "                        remove_rot_trans = True,\\\n",
    "                        asr = phonon_info.asr)\n",
    "\n",
    "# We write the informations using write method in anharm object.\n",
    "anharm.write(file_path='co2-100K-piglet',atoms=phonon_info.atoms)\n",
    "\n",
    "# Task complete, we are deleting the anharm object to finish file writing\n",
    "del anharm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84779c75",
   "metadata": {},
   "source": [
    "This will create two files storing the mode-resolved and atom resolved anharmonic measures.\n",
    "\n",
    "Let us first look at the mode-resolved measure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2acfee66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Mode  Freq(cm-1)       Var(Anh)        Var(Tot)     Anh_Measure \n",
      "#--------------------------------------------------------------\n",
      "    1     -0.0000    5.97692e-06     6.12601e-06       0.987757 \n",
      "    2     -0.0000    6.31072e-06     1.97373e-06        1.78812 \n",
      "    3     -0.0000    0.000206224     8.77498e-06        4.84782 \n",
      "    4      0.0000    6.15809e-06     6.91862e-06       0.943438 \n",
      "    5      0.0000    6.09401e-06     5.90723e-05       0.321188 \n",
      "    6    681.8844    0.000174674     0.000346854       0.709644 \n",
      "    7    681.8851    6.19722e-06     0.000368153       0.129743 \n",
      "    8   1345.3420    5.64072e-06      0.00337331      0.0408921 \n",
      "    9   2367.9766    0.000259688       0.0141033       0.135695 \n",
      "#--------------------------------------------------------------\n",
      "# SUM                7.52182e-05      0.0020305       0.192469\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat co2-100K-piglet_mode_res_anh_mes.out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2decfd52",
   "metadata": {},
   "source": [
    "We can neglect the first 5-modes (translation & rotation). Now 6th and 7-th modes are degenerate in-plane and out-of-plane bending. Being degenerate, any linear-combination of the obtaned eigenvectors is also a possible description of these normal modes and therefore, there anharmonic measure should not be reported separately but the average anharmonic measure of such Eg modes must be reported. \n",
    "\n",
    "The anharmonic measure is the square-root of the ratio of the varience of anharmonic forces (column 3) and the varience of total force (column 4). Therefore, the average anharmonic measure of the Eg modes can be calculated as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "55ad0aa4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Anharmonic Measure of Bending Modes = 0.5029555295307137\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "# Now we read force variences (columns 3 & 4) from the above file\n",
    "f_var = np.genfromtxt(\"co2-100K-piglet_mode_res_anh_mes.out\", usecols=(2,3)) \n",
    "anh_mes = np.sqrt( (f_var[5,0] + f_var[6,0]) / (f_var[5,1] + f_var[6,1]) )\n",
    "print(\"Anharmonic Measure of Bending Modes = \" + str(anh_mes) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4909887",
   "metadata": {},
   "source": [
    "Similarly we can also check the anharmonic measure of each atom separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "adc83dad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#  Atom         Var(Anh)        Var(Tot)     Anh_Measure \n",
      "#--------------------------------------------------------\n",
      "    1   O     7.28372e-05      0.0014313       0.225586\n",
      "    2   O     6.23086e-05     0.00135273       0.214619\n",
      "    3   C     9.05087e-05     0.00330748       0.165423\n",
      "#-------------------------------------------------------\n",
      "# SUM         7.52182e-05      0.0020305       0.192469\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat co2-100K-piglet_atom_res_anh_mes.out"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1c4dbc7",
   "metadata": {},
   "source": [
    "If we want to obtain average anharmonic measure of a group of atoms, e.g., O atoms in CO2 we can compute the average anharmonic measure just we did it for bending modes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f052e3e5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Anharmonic Measure O atoms = 0.22032527522933562\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "# Now we read force variences (columns 3 & 4) from the above file\n",
    "f_var = np.genfromtxt(\"co2-100K-piglet_atom_res_anh_mes.out\", usecols=(2,3)) \n",
    "anh_mes = np.sqrt( (f_var[0,0] + f_var[1,0]) / (f_var[0,1] + f_var[1,1]) )\n",
    "print(\"Anharmonic Measure O atoms = \" + str(anh_mes) )"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}