{ "cells": [ { "cell_type": "markdown", "id": "670ad1a9", "metadata": {}, "source": [ "# 1.6.1 Stochastic Simulation Using ORCA " ] }, { "cell_type": "markdown", "id": "8703e0d0", "metadata": {}, "source": [ "In this example we will learn to do a stochastic Monte Carlo calculation for the band gap renormalization (similar to exercise 1.4.1). However, here our choice of ab initio code would be ORCA. As ORCA has routines that can efficiently calculate phonons; we will use it. \n", "\n", "As, usual the first step is to start with an guessed atomic cordinates and perform geometry optimizations. Then we take the geometry-optimized structure and perform frequency calculations that will give us the normal mode. These two steps can be performed at once using the following orca-input file: orca-opt-freq.inp (see orca manual: https://www.faccts.de/docs/orca/6.0/manual/ for mode details). " ] }, { "cell_type": "code", "execution_count": 1, "id": "afa4e7ab", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "!B3LYP DEF2-TZVP OPT FREQ TIGHTSCF\n", "\n", "%maxcore 28000\n", "\n", "%pal nprocs 12 end\n", "\n", "*xyzfile 0 1 co2-init.xyz\n" ] } ], "source": [ "%%bash\n", "cat orca-opt-freq.inp" ] }, { "cell_type": "markdown", "id": "61365afb", "metadata": {}, "source": [ "The calculation was submitted using orca-opt-freq.sbatch file to the midway3 cluster of RCC (University of Chicago). After finishing, orca would create several files." ] }, { "cell_type": "markdown", "id": "2433c861", "metadata": {}, "source": [ "The orca-opt-freq.hess file contains all information that we need. We would convert this file to PyEPFD restart xml file using a command line tool avaiable at utils/orca. This command line tools takes 2 mandatory arguments in the following order: (1) source path, i.e., the full path to the .hess file we would like to convert and (2) destination path, i.e., full path to the .xml file which we would like to obtain. Therefore we would use the fillowing command:" ] }, { "cell_type": "code", "execution_count": 2, "id": "635741bc", "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-19 16:23:35.058617\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", "-----------------------------------------------------\n", "Converting orca-opt-freq.hess to co2.xml\n", "Frequency scale factor = 1.0\n", "-----------------------------------------------------\n", "#Mode Frequency(cm-1)\n", " 1 -0.02\n", " 2 -0.00\n", " 3 0.00\n", " 4 0.25\n", " 5 0.25\n", " 6 676.43\n", " 7 676.43\n", " 8 1371.00\n", " 9 2408.06\n", "Time spent at write_info 0.0004 s\n", "Successfully written co2.xml.\n", "Change the values of , , , \n", " and if needed.\n", "Total time required (s): 0.0018100738525390625\n" ] } ], "source": [ "%run ../../../utils/orca/orca2pyepfd.py orca-opt-freq.hess co2.xml" ] }, { "cell_type": "markdown", "id": "f43c046e", "metadata": {}, "source": [ "The newly created co2.xml file is similar to that we saw before. As orca is a non-periodic code, while writing the cell values, PyEPFD puts a cubic cell with a large value (100 angstrom) of the cell parameter. If required, you can change it accordingly. " ] }, { "cell_type": "code", "execution_count": 3, "id": "ad2c596b", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " \n", "[ 1.88972599e+02, 1.88972599e+02, 1.88972599e+02, 9.00000000e+01, 9.00000000e+01, \n", " 9.00000000e+01 ]\n" ] } ], "source": [ "%%bash\n", "grep -A2 \" prep_orca.log\n", "cd ../" ] }, { "cell_type": "markdown", "id": "523eed52", "metadata": {}, "source": [ "Let us have a look at the contents of the following directory." ] }, { "cell_type": "code", "execution_count": 9, "id": "14de4bf2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "coord.xyz\n", "orca-sp.inp\n" ] } ], "source": [ "%%bash\n", "ls MCAP_100K/frame-1/" ] }, { "cell_type": "code", "execution_count": 12, "id": "96a5f258", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "========coord.xyz=========\n", "3\n", "# CELL(abcABC): 100.00000 100.00000 100.00000 90.00000 90.00000 90.00000 PyEPFD-Step: 0 positions{angstrom} cell{angstrom}\n", "O -0.0060215693 0.024954011 1.1515462\n", "O -0.0060215691 0.024954015 -1.1539117\n", "C 0.016042212 -0.066480597 0.0031511087\n", "=======orca-sp.inp========\n", "!B3LYP DEF2-TZVP EnGrad\n", "\n", "%maxcore 28000\n", "\n", "%pal nprocs 12 end\n", "\n", "*xyzfile 0 1 coord.xyz\n" ] } ], "source": [ "%%bash\n", "echo \"========coord.xyz=========\"\n", "cat MCAP_100K/frame-1/coord.xyz\n", "echo \"=======orca-sp.inp========\"\n", "cat MCAP_100K/frame-1/orca-sp.inp" ] }, { "cell_type": "markdown", "id": "1f134a27", "metadata": {}, "source": [ "We see that geometry of the first frame and dummy orca input file has been copied within the frame-1 directory. From the MCAP_100K directory, we can submit the jobs using a job script similar to below:" ] }, { "cell_type": "code", "execution_count": null, "id": "2bfd7e90", "metadata": {}, "outputs": [], "source": [ "# %load MCAP_100K/orca.sbatch\n", "#!/bin/bash\n", "#SBATCH --account=pi-gagalli\n", "#SBATCH -J co2\n", "#SBATCH --time=24:00:00\n", "#SBATCH --partition=gagalli-csl2\n", "#SBATCH --qos=gagalli-csl2\n", "#SBATCH --nodes=1\n", "#SBATCH --ntasks-per-node=12\n", "#SBATCH --mem-per-cpu=28GB\n", "#SBATCH --cpus-per-task=1 #increase this comensurate with number of threads\n", "#SBATCH --array=0-3\n", "export nframe=100\n", "\n", "#\n", "# SET OMP_NUM_THREADS \n", "export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK\n", "#\n", "module load orca/5.0.4\n", "orca=/software/orca-5.0.4-el8-x86_64/bin/orca\n", "#\n", "# NTASKS BELOW IS SET FOR GNU PARALLEL TO KNOW HOW MANY CORES IT HAS TO WORK WITH\n", "NTASKS=$(($SLURM_NTASKS_PER_NODE * $SLURM_JOB_NUM_NODES))\n", "echo \"NTASKS IS $NTASKS\"\n", "\n", "start_time=$(date +%s)\n", "\n", "# Calculating start and end frame indices from job array id.\n", "let begin=${nframe}*${SLURM_ARRAY_TASK_ID}+1\n", "let end=${nframe}*${SLURM_ARRAY_TASK_ID}+${nframe}\n", "\n", "for (( i=$begin; i<=$end; i++ )); do\n", " cd frame-${i}\n", " ${orca} orca-sp.inp > orca-sp.out\n", " cd $SLURM_SUBMIT_DIR\n", "done \n", "\n", "end_time=$(date +%s)\n", "total_time=$((end_time - start_time))\n", "echo \"Total execution time: $total_time seconds\"\n" ] }, { "cell_type": "markdown", "id": "213c893e", "metadata": {}, "source": [ "The above script submit 4 array jobs handling frames: 1-100, 101-200, 201-300 and 301-400 respectively. Once all jobs are finished, we would extract the data using another command line tool extract_properties.py. The command line options can be obtained by using -h (or --help) options like below. " ] }, { "cell_type": "code", "execution_count": 1, "id": "7b67be18", "metadata": { "scrolled": true }, "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-19 18:59:59.687544\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", "================================================\n", " ORCA EXTRACT_PROPERTIES HELP MESSAGES \n", "=================================================\n", "\u001b[91mUsage (optional arguments in [...]):\n", "extract_properties.py -p \u001b[92m\u001b[91m -n \u001b[92m\u001b[91m [-f \u001b[92m\u001b[91m] [-o \u001b[92m\u001b[91m] [-a \u001b[92m\u001b[91m] \n", "\u001b[00m\n", "Command line argunemts: \n", "\n", "\u001b[91m -h, --help =\u001b[00m prints help message.\n", "\n", "\u001b[91m -p, --path =\u001b[00m \u001b[92m\u001b[00m to the orca calculations.\n", " For a single frame \u001b[92m\u001b[00m of the directory where .out file \n", " is present.\n", " For many frames, \u001b[92m\u001b[00m of the directory where individual \n", " frame directories exist.\n", "\n", "\u001b[91m -o, --orca =\u001b[00m A string defining \u001b[92m\u001b[00m of orca outputs.\n", " The default is 'orca-sp'.\n", "\n", "\u001b[91m -n, --n_orb =\u001b[00m The range of orbital indices(\u001b[92m\u001b[00m) to be extracted.\n", " The indices are the same as with orca-output (starts with 0)\n", " The argument \u001b[92m\u001b[00m may take 3 integers separated by ':'.\n", " These integers define \u001b[92mstart:end:step\u001b[00m of the orbital indices.\n", " [a] only providing \u001b[92mstart\u001b[00m would extract start orbital only.\n", " [b] if only \u001b[92mstart:end\u001b[00m provided then \u001b[92mstep = 1\u001b[00m.\n", "\n", "\u001b[91m -f, --frames =\u001b[00m The range of \u001b[92m\u001b[00m.\n", " The argument \u001b[92m\u001b[00m may take 3 integers separated by ':'.\n", " These integers define \u001b[92mstart:end:step\u001b[00m of the trajectory frames.\n", " [a] only providing \u001b[92mstart\u001b[00m would extract start frame only.\n", " [b] if only \u001b[92mstart:end\u001b[00m provided then \u001b[92mstep = 1\u001b[00m.\n", "\n", "\u001b[91m -a, --atoms =\u001b[00m The list of \u001b[92m\u001b[00m (starting with 0 as orca output)\n", " for which we want to extract (charge/spin) populations.\n", " The argument \u001b[92m\u001b[00m either takes\n", " (1) A string either \n", " [a] 'none' (default): populations wont be extracted, \n", " or [b] 'all': populations of all atoms would be extracted. \n", " or\n", " (2) many integers separated by ',' defining the list of atom indices,\n", " or\n", " (3) Integers (max. 3) separated by ':' defining ranges.\n", " See the frames above for definition of ranges.\n" ] } ], "source": [ "%run ../../../utils/orca/extract_properties.py --help" ] }, { "cell_type": "markdown", "id": "7c5ef7b0", "metadata": {}, "source": [ "Now we would extract various properies for frames 1 to 400. We extract degenerate HOMOs (orbital indices 9,10) and LUMO (orbital index 11). Therefore we use the following command." ] }, { "cell_type": "code", "execution_count": 2, "id": "b9f098b2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Extracting properties from orca calculations\n", "Frame range: \u001b[94m(1, 400, 1)\u001b[00m\n", "\u001b[94m40/400\u001b[00m frames processed in \u001b[94m 1.638\u001b[00m s.\n", "\u001b[94m80/400\u001b[00m frames processed in \u001b[94m 3.223\u001b[00m s.\n", "\u001b[94m120/400\u001b[00m frames processed in \u001b[94m 5.089\u001b[00m s.\n", "\u001b[94m160/400\u001b[00m frames processed in \u001b[94m 6.934\u001b[00m s.\n", "\u001b[94m200/400\u001b[00m frames processed in \u001b[94m 8.545\u001b[00m s.\n", "\u001b[94m240/400\u001b[00m frames processed in \u001b[94m 9.672\u001b[00m s.\n", "\u001b[94m280/400\u001b[00m frames processed in \u001b[94m 10.186\u001b[00m s.\n", "\u001b[94m320/400\u001b[00m frames processed in \u001b[94m 11.825\u001b[00m s.\n", "\u001b[94m360/400\u001b[00m frames processed in \u001b[94m 13.378\u001b[00m s.\n", "\u001b[92m\n", "Property extraction finished in 14.846 s.\u001b[00m\n", "*************************************************\n", "Total single-point energies are written at: \n", "\u001b[92m MCAP_100K//etotals.dat\u001b[00m\n", "Orbital energies are written at : \n", "\u001b[92m MCAP_100K//is_0_orb-9-11.dat\u001b[00m\n", "*************************************************\n", "EnGrad informations are written at: \n", "\u001b[92m MCAP_100K//engrad.npz\u001b[00m\n", "Keywords for the engrad.npz file are:\n", " atoms = a list of atoms symbols,\n", " frames = a list of original frame indices,\n", " coords = a 2D numpy of coordinates (au),\n", " forces = a 2D numpy of forces (au),\n", " etotals = a 1D numpy of energies (au)\n", "In the last 3 arrays, \n", " 0-th dimension corresponds to frames\n", "*************************************************\n" ] } ], "source": [ "%run ../../../utils/orca/extract_properties.py -p MCAP_100K/ -n 9:11 -f 1:400 -o orca-sp" ] }, { "cell_type": "markdown", "id": "5e25b704", "metadata": {}, "source": [ "Like in example 1.4.1, we compute the band gap using an awk script." ] }, { "cell_type": "code", "execution_count": 3, "id": "05d914fe", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10.0558\n", "9.883\n", "10.1485\n", "10.313\n", "10.4468\n", "10.0803\n", "10.1678\n", "10.4725\n", "10.3487\n", "10.5374\n", "10.5014\n", "9.85465\n", "9.73425\n", "10.5914\n", "10.4057\n", "10.0914\n", "10.6387\n", "10.756\n", "9.55225\n", "10.3976\n", "10.3003\n", "10.2356\n", "9.7635\n", "10.2627\n", "10.2915\n" ] } ], "source": [ "%%bash\n", "awk '{print $3-($1+$2)/2}' MCAP_100K/is_0_orb-9-11.dat | sed '1d' > gap.dat\n", "head -n 25 gap.dat" ] }, { "cell_type": "markdown", "id": "b96b39eb", "metadata": {}, "source": [ "As usual, to obtain the renormalized gap, the next step is to average over all 400 configurations." ] }, { "cell_type": "code", "execution_count": 4, "id": "c4edaa8e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAYAAABB4NqyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAACP00lEQVR4nOzdd3ib5dU/8O+jvSXvveI4e++QATSBJC+FEPq2rEICFFrKTiiQFgiltKGhTSmFsloIvECh/Aih0DaMQBJIQshy9rLjvZcka4/n/v0hW7Zs2ZZsyZKl87muXCDpkXTLsvUcnfvc5+YYYwyEEEIIIXFEEOkBEEIIIYQMNwqACCGEEBJ3KAAihBBCSNyhAIgQQgghcYcCIEIIIYTEHQqACCGEEBJ3KAAihBBCSNwRRXoA0YjnedTW1kKtVoPjuEgPhxBCCCEBYIyhvb0dmZmZEAj6z/FQAORHbW0tcnJyIj0MQgghhAxCVVUVsrOz+z2GAiA/1Go1AM8PUKPRRHg0hBBCCAmE0WhETk6O9zzeHwqA/Oic9tJoNBQAEUIIISNMIOUrVARNCCGEkLhDARAhhBBC4g4FQIQQQgiJOxQAEUIIISTuUABECCGEkLhDARAhhBBC4g4FQIQQQgiJOxQAEUIIISTuUABECCGEkLhDARAhhBBC4g4FQIQQQgiJOxQAEUIIISTuUABECCGEkLhDARAhhBBC4g4FQFGEMTbo+/I8Q2O7LYSjIYQQQmIXBUBRxMUPPgBqaLfBYHGGcDSEEEJI7KIAKIq43IMPgCpaLLA63SEcDSGEEBK7KACKIk6eH9T9Ws0OmGwu2JyDuz8hhBASbygAiiKDzQBVtJgBADbKABFCCCEBoQAoirjcwWdwHC4eLSaH9//dQ6gjIoQQQuIFBUBRxDmI4KXV7PC5TFkgQgghZGAUAEWRwWSAmk12n8sUABFCCCEDowAoijiDrAFijPXKANFKMEIIIWRgFABFEVeQq8BMdhccLt/70EowQgghZGAUAEURl5sF1Q26s/i5O5oCI4QQQgZGAVAUYSy4btAtZgqACCGEkMGgACiKMLCAl7G73DwM1t4BENUAEUIIIQOjACiKBJMBMlid8Fcy5HDx4KkXECGEENIvCoCiCAPgDnAlmN7qf+NTxgC7iwqhCSGEkP5QABRFGGMBrwQz9BEAAVQHRAghhAyEAqAowoCAaoAYYzD2EwBRHRAhhBDSPwqAoognAzRwAGR2uPvdOJUyQIQQQkj/KACKIowFlgHqb/oLoAwQIYQQMhAKgKIIQ2CrwAyW/gMgs50CIEIIIaQ/FABFEU8GaOAiaL2f/j/dmezOoDpKE0IIIfGGAqAoEkgNkNPNwzJAhofnAYuDskCEEEJIXygAiiIM6Le4GRi4/qdTu80VghERQgghsYkCoCgSSBG02R5YYGMK8DhCCCEkHlEAFEUYBp4CC3SXi3abb6bI6abu0IQQQkgnCoCiSQAZoEA3S+2eATLbXShvNg9paIQQQkgsoQAoiniWwfefqeEDXN1ld/JwdOwJVmewoUZvhYuyQIQQQggACoCiCs/YgBmeQAMgoCsL1GC0weVmqDPYhjQ+QgghJFZQABRFGBu4EWKgU2AAYLK5oLc4YO1YEl/VaqH+QIQQQggoAIoqDIB7gGXwwcQv7XYn6o1dWR+Lw41mU/9NFAkhhJB4QAFQFOnMzvRXqxNMBshodaG+x7RXRQsVQxNCCCEUAEWRztCmv2mwYGqAzHZXr8aKeotzwL3ECCGEkFhHAVA06YhV+svyBBMA9aWcskCEEELiHAVAUYR1RED9Z4CG/jxN7faAO0oTQgghsYgCoCjCAsgABVMD1J+KFktIHocQQggZiSgAiiKdAVB/zRD5EAVAegutBiOEEBK/KACKEt378/RfAxSa57M43LQ/GCGEkLhFAVCU6F7b3HPlVnfuEDYybLdRHRAhhJD4RAFQlOge1oR7FVinnjvGE0IIIfGCAqAo0X0KrL9VYKHcyoIyQIQQQuIVBUBRIpAMEGMMA2wWHxSjlTJAhBBC4lNEA6Ddu3fjyiuvRGZmJjiOw7Zt23xuZ4zh8ccfR0ZGBuRyOZYuXYrz588H/PhPP/00OI7D/fffH9qBhwHvkwHyH+WEqgC6k8Xh7nfbDUIIISRWRTQAMpvNmDp1Kl544QW/t2/atAnPPfccXnrpJezfvx9KpRLLli2DzWbze3x3Bw4cwMsvv4wpU6aEethh0X1mq68MUKh6AHVH02CEEELiUUQDoBUrVuCpp57CqlWret3GGMOzzz6LRx99FCtXrsSUKVPw5ptvora2tlemqCeTyYQbb7wRr776KhISEsI0+vDpqwYolAXQnSgAIoQQEo+itgaorKwM9fX1WLp0qfc6rVaLuXPnYt++ff3e96677sIVV1zhc9/+2O12GI1Gn3/DLZBl8OEIgIy0EowQQkgcEkV6AH2pr68HAKSlpflcn5aW5r3Nn3fffReHDx/GgQMHAn6ujRs34te//vXgBhoirFsZtN3l9ntMOKbAKAAihBASj6I2AzQYVVVVuO+++/D2229DJpMFfL/169fDYDB4/1VVVYVxlP51T+44XLzf5e5hiH9gsbtxqKIVNXprWAIsQgghJBpFbQYoPT0dANDQ0ICMjAzv9Q0NDZg2bZrf+xw6dAiNjY2YMWOG9zq3243du3fj+eefh91uh1Ao7HU/qVQKqVQa2hcQpO6hB2OA080gEXG+x4RhCgwA2sxOtJmdEAk4pGkCDxwJIYSQkSpqM0AFBQVIT0/Hjh07vNcZjUbs378f8+fP93ufJUuW4Pjx4yguLvb+mzVrFm688UYUFxf7DX6iRc/gxt80WLgzNFQQTQghJF5ENANkMplQUlLivVxWVobi4mIkJiYiNzcX999/P5566ikUFRWhoKAAjz32GDIzM3H11Vd777NkyRKsWrUKd999N9RqNSZNmuTzHEqlEklJSb2ujzY9Qxu7i4e6x3Wh3AfMH9oagxBCSLyIaAB08OBBXHrppd7La9euBQCsXr0aW7ZswUMPPQSz2Yw77rgDer0eCxcuxPbt233qe0pLS9Hc3DzsYw+1nrGNw9W7QWGY4x+Y7JQBIoQQEh84Fq7CkhHMaDRCq9XCYDBAo9EMy3PqLQ4cLG/zXi5MVaEgWelzTK3eilO14V2iv2hMMqSi6J0qJIQQQvoSzPk7amuA4p2/DFA4+gD1ZKI6IEIIIXGAAqAo0bO+2V8RdCg3Qu1L90JoSg4SQgiJVRQARYneq8AilAHqVgdU0WIJ+/MRQgghkUABUJToGdr4mwIL9yowoKsztNPNo6zZ3GdXakIIIWQkowAoSvSMbfwFHsMxJWWxu+HmGWo7OkMbrVQTRAghJPZQABQlWI8cEM97sjDduYehBgjw9AOqarV6/58QQgiJNRQARQs/yZ2edUDDtVdXaZMZNqcnA2WkVWGEEEJiEAVAUcJfaNOzDmg4iqABoM3s8P4/ZYAIIYTEIgqAooS/2KZnHdBwBUA+Y3DyVAhNCCEk5lAAFCV61gABnuCju2GaAeuFCqEJIYTEGgqAooS/4MbhjswUWE80DUYIISTWUAAUxXplgCKUAqJCaEIIIbGGAqAo4S+4cbh9a2+GaxVYT5QBIoQQEmsoAIpi0VIDZHfyMNspC0QIISR2UAAUJfyuAouSGiAAqDPYIvbchBBCSKhRABQl/K0Cc7sZXN2CoEgGQA1GCoAIIYTEDgqAokRfsU33lWCRqgECAKvDDb3FMfCBhBBCyAhAAVCU6Cu06d4NOoIJIABAPWWBCCGExAgKgKJEXzu9dw+AIpkBAoB6gy1iS/EJIYSQUKIAKEr0FVZ0bogaDYGHy83QYqZpMEIIISMfBUBRYqAaIHek57861NNqMEIIITGAAqCo0f8UWCRXgHXXbLL7rEwjhBBCRiIKgKJEXzNc3gAoSmION8/QZLJHehiEEELIkFAAFCX6SvA43dGVAQKoKSIhhJCRjwKgKOGvESLQVQQdLTVAANBmdsDucsPNM5Q2mfpcwUYIIYREK1GkB0A8+iyC7giAWJRMgQGesZY1m9FmdsJsdyFZJYVWLo70sAghhJCAUQYoSvQVALl55vkXZVmW6lard4PUpnaqCSKEEDKyUAAUJfqaAgM8WaBoqgHqqZmKogkhhIwwFABFif7iG4eLj4pGiH0x2VywOtyRHgYhhBASMAqARgC7293nMvlo0X0aLNJbdhBCCCEDoQAoSgyUAYq2GqCemkx2MMZwvqEdJY2mSA+HEEII6RcFQFGivxqfaJ8CAwC9xYHDlW2oaLGgsd1GS+MJIYRENQqARgCHO7qLoAFPBqvN7AQA2J08DFZnhEdECCGE9I0CoCgxUAZopNXVNIZ5abyT9iMjhBAyBBQARYn+whvPMvhhG0pINBrDFwAZbU5caDKH7fEJIYTEPgqAosSAy+CjfAqsJ5vTHbZpsAtNZjQYqc6IEELI4FEAFDX6PpnbR0ANkD8NxtBvmmqwONHcbofDxaPNQnVGhBBCBocCoCjRX3zjdjM4XSMvAKpqtYQ8C1Ta3LXEvp52pSeEEDJIQW+GWlZWhq+//hoVFRWwWCxISUnB9OnTMX/+fMhksnCMMS4MFN7YXCOv0zJjwMkaA+YUJEIkHHqs7XDxaDU5vJcb220Yx6shEHBDfmxCCCHxJeAA6O2338af//xnHDx4EGlpacjMzIRcLkdraytKS0shk8lw44034uGHH0ZeXl44xxyTBprhsjlHXgAEABaHG2fq2zEpSzvkx+rcfLWTy83QanEgWSUd8mMTQgiJLwF9LZ8+fTqee+45rFmzBhUVFairq8OhQ4fwzTff4NSpUzAajfjoo4/A8zxmzZqF999/P9zjjjkDFfTanSN32XeD0QZXCJatm3oEQEB4V5sRQgiJXQFlgJ5++mksW7asz9ulUikuueQSXHLJJfjtb3+L8vLyUI0vboy8Cp/AMQborc4hZ2osfjZcbbM4/BxJCCGE9C+gDNCyZcvQ2toa0AMmJSVh5syZQxpUPBqBi7yCog/Bii1/GSCrw0070RNCCAlawJWpmZmZuO666/D555+Hczxxa6g9bRqMNnx2qh72KC2W1ocgU9OzBqhTK2WBCCGEBCngAOjVV19FU1MTli9fjvz8fDzxxBM01RVCQ00AvX+wGv88WI1DFW0hGU+oGW3OIW3n4XDxcLj81xF1XxlGCCGEBCLgAOimm27Cjh07UFJSgtWrV+ONN97A6NGjcdlll+G9996Dw0EnoaEY6hRYWYtna4ho3YSU54c2tr6yP4AnA0RdoQkhhAQj6OYsBQUF+PWvf42ysjJs374dqampuPXWW5GRkYF77703HGOMC2wIOSCD1ekNLqK5HmYo02D+6n86OV18v7cTQgghPQ2pO93SpUvx9ttv48033wQAvPDCCyEZVDwaSgKjqtXi/X9/K6Wihb+tKwLt5mx29B/gtJmjM/NFCCEkOg06AKqoqMATTzyBgoICXHvttZgxYwbefvvtUI4trgxlAqdyhARARqsTfI86oAvNJjgD6BHU3xQYQIXQhBBCghPUVhh2ux0ffPABXnvtNezcuRNZWVlYs2YNbrnlFuTn54dpiPFhKDUsPgGQM3qngtw8Q7PZjlS1Z8sUs90Fi92NNovDe11fzPb+A7s2swM8z2hbDEIIIQEJOAD6+c9/jnfffRcWiwUrV67Ef/7zH1x22WXgODrhhMJQMkDdp8CiuQYIACpaLN5gp7Hd08XZYHH2GwD1twKsk5tnaLM4kETbYhBCCAlAwAHQN998gw0bNuDHP/4xkpKSwjmm+DTICMjmdHsDCSC6p8AAT7CjtzigU0jQaPTU//irDequ3RZYfU+ziQIgQgghgQk4ADp27JjPZYfDgbKyMhQWFkIkCnpTedLDYFeBVbVZfO4Z7QEQAJQ1mzE+Q4h2m2e6rt3mhMvN97ljfEW3DFd/mk12jIU6ZOMkhBASu4IugrZYLLjtttugUCgwceJEVFZWAgDuuecePP300yEfYLwYbAlQVasVAJCp9UwhRfsUGAC0mBwobTJ5LzPWd48gvcURcKNDq8M9YLE0IYQQAgwiAFq/fj2OHj2KnTt3QibrqttYunQp3nvvvZAOLp4MNgDqLIAem+7JfDjcfEh2Xg+3Or3v8ve+psFKm8xBPW6ziXaHJ4QQMrCgA6Bt27bh+eefx8KFC30KoCdOnIjS0tKQDi5ehGIF2Ji0rqkfqzP6s0A9Gay9szytZgfazMEtb6cAiBBCSCCCDoCampqQmpra63qz2UwrwgZpsPEPzzPU6j1TYHlJCsjFQgAjow6oJ4PV6TN9xfMM5xrag34cvcU54IoxQgghJOgAaNasWfj3v//tvdwZ9Pztb3/D/PnzQzeyOMIPMgJqMTvg4hlEAg7JSinkkpEbAPE8cLzG4G2UeKHZDJMt+HoexuANCgkhhJC+BL1863e/+x1WrFiBU6dOweVy4c9//jNOnTqFvXv3YteuXeEYY8wb7ARYQ8cy8lSNFAIBB4VEiFYzYBlg24hoZbK5cL7RhHSNDBUtwdX+dFfdZkVekoIykoQQQvoUdAZo4cKFKC4uhsvlwuTJk/HZZ58hNTUV+/btw8yZM8Mxxpg32CmwzgAoTeMpRu+cAhsJK8H6UtVqwdFq/ZD2RrM53WiiWiBCCCH9GFQDn8LCQrz66quhHkvcGmwPoAaj5ySf1tFFWTGCp8C6C0UNT1WrdcDtNQghhMSvgDJAZnNw0xHBHh/vhpoBStd0BkCeeHakB0Ch0GZ2wEQ9gQghhPQhoABo9OjRePrpp1FXV9fnMYwxfP7551ixYgWee+65gJ589+7duPLKK5GZmQmO47Bt27Zej/n4448jIyMDcrkcS5cuxfnz5/t9zBdffBFTpkyBRqOBRqPB/Pnz8d///jeg8Yw0De2dU2Ce7R+8RdBRvCHqcOrcaoMQQgjpKaApsJ07d+KXv/wlnnjiCUydOhWzZs1CZmYmZDIZ2tracOrUKezbtw8ikQjr16/HT3/604Ce3Gw2Y+rUqbj11ltxzTXX9Lp906ZNeO655/DGG2+goKAAjz32GJYtW4ZTp075NGHsLjs7G08//TSKiorAGMMbb7yBlStX4siRI5g4cWJA4xpug8kAOVw8Wjo6JKdpfKfARnINUCi1mh0YldL37S43j1N1RoxNV0MqEg7fwAghhERcQAHQ2LFj8cEHH6CyshLvv/8+vv76a+zduxdWqxXJycmYPn06Xn31VaxYsQJCYeAnkhUrVmDFihV+b2OM4dlnn8Wjjz6KlStXAgDefPNNpKWlYdu2bbjuuuv83u/KK6/0ufzb3/4WL774Ir799tvoDYAGUQPU1G4HgyfoUcs8b2Os1ACFisHa/x5jZ+rb0Wi0o93mwvRcnXcKkRBCSOwL6hM/NzcX69atw7p168I1Hq+ysjLU19dj6dKl3uu0Wi3mzp2Lffv29RkAded2u/H+++/DbDb326PIbrfDbu9aNWQ0Goc2+CANJgNU320FWOdyb4XY83ZSBsiDMc8WGynq3jvE1xmsqDd4foZWhxtHqwyYX5g03EMkhBASIUEvgx8u9fX1AIC0tDSf69PS0ry39eX48eNQqVSQSqX42c9+hg8//BATJkzo8/iNGzdCq9V6/+Xk5Az9BQRhMDXQXUvgu07uI7kRYri0WXpvpcHzDGfrfbtMm+0utAa57QYhhJCRK2oDoKEYO3YsiouLsX//ftx5551YvXo1Tp061efx69evh8Fg8P6rqqoaxtEObi+wnj2AgG5TYFQE7dXiZyf5ZpMdLnfvn3l1m2U4hkQIISQKRG3RQ3p6OgCgoaEBGRkZ3usbGhowbdq0fu8rkUgwevRoAMDMmTNx4MAB/PnPf8bLL7/s93ipVAqptPc0yXAZXAbItwcQ0JUBoimwLma7CzanGzJxV21ancH/6rCmdnuvYwkhhMSmqM0AFRQUID09HTt27PBeZzQasX///qD3HON53qfGJ9qwQfT961wCn+4vA9RPAGRzuvHfE3XY8K+T2HWuKfgnHoG6T4M53TxazP5/FxjzbKNBCCEk9kU0A2QymVBSUuK9XFZWhuLiYiQmJiI3Nxf3338/nnrqKRQVFXmXwWdmZuLqq6/23mfJkiVYtWoV7r77bgCe6awVK1YgNzcX7e3teOedd7Bz5058+umnw/3yAhbsKjCz3YX2jo1CU7vVAHUWQdtdPFw8D5HAN77dW9qMfx6s9jYI/PhoLRYXJcf8nlm1eivS1DIIBBwa2+3g+wk4a/RWFCQrIRTE9s+EEELi3aACoLa2Nvz973/H6dOnAQDjx4/HrbfeisTExKAe5+DBg7j00ku9l9euXQsAWL16NbZs2YKHHnoIZrMZd9xxB/R6PRYuXIjt27f79AAqLS1Fc3Oz93JjYyNuvvlm1NXVQavVYsqUKfj0009x2WWXDealDotgS4A6sz9audhnuqZzCgzwTIOpZZ4AyO504+3vKrG3tAUAkKqWQm91Qm91oqLVgvwk5RBfQXRrMztxtFqPqdk61Bv6z/A4XTyqWi3IT47tnwkhhMQ7jgVZgbt7925cddVV0Gg0mDVrFgDg0KFD0Ov1+Pjjj7F48eKwDHQ4GY1GaLVaGAwGaDSasD9fq9mBwxVtAR9/oLwVL+++gKJUFR5ePs7ntrveOQy7i8dvr57kLZB+7svzOFZtAMcBV03NxP9MysAruy/gUGUbrpySgZXTssDzDG7GIO6jZ04s0CnE0FucAx4nEnJYMDo5pn8WhBASi4I5fwf9CX/XXXfhRz/6EcrKyrB161Zs3boVFy5cwHXXXYe77rpr0IOOZ8GuAmvu2Ok8SSXpdZu3G7TTUwfkcPE4UWMAANy/pAhXTsmEUMBhao4WAFBcpQdjDM9/VYL73iv2PnYsCiT4AQCXm6GihVaEEUJILAs6ACopKcG6det8Oj4LhUKsXbvWp56HBC7YVWDNHUu7k5W9V651djPuXAlW3WYBzwCNTIQJGV3R8OQsLTgOqGqz4l9Ha3GsxuATLMW7qjYL7C5aTUcIIbEq6ABoxowZ3tqf7k6fPo2pU6eGZFDxJtgaoJYAMkCdK8E6Mxm5SQqfYme1TIzRKSoAwMfHuja5LafMBwDA7e7dLJEQQkjsCLoI+t5778V9992HkpISzJs3DwDw7bff4oUXXsDTTz+NY8eOeY+dMmVK6EYaw4JdBdbc0bE4WdU7A9TVDdqz0qui1RPQ5Cf2LuqdlqPD+UYTAEAqEsDu4lHWbPbe/7VvyjE9V4cFo5ODGl+saDTaUWewIkMrj/RQCCGEhFjQAdD1118PAHjooYf83sZxHBhj4DgObjdNIQQkiPiHMdZvBkgu7pkB8gQ0uUmKXsdOzdHh/UPV4AD8ZGEBXthZilqDFXanG3tLW1BcrUdlmyVuAyDAs2GqTi7xWWFHCCFk5As6ACorKwvHOOIaH0QAZLS54HQzcAASFf0UQTvccLp51Oo9S+bzEnsHQOkaGW5fWACxSIDpuQnQycXepfFHq/UAPCvUjFYnNHJx0K8rFrjdDGcb2jEtRxfpoRBCCAmhoAOgvLy8cIyDBKgz+5OgkEDkZ5l2ZxG0xeFGjd4KN2NQSUVIVPYOlgBg7qiuHdALkpU4UqXH6TojzjWYvNeXt5gxJVsXwlcxsjS326G3OKDzE3ASQggZmQbdCfrUqVOorKyEw+G72eRVV1015EHFGz6IKujOFWD+pr+ArgxQm9XRVQCdqAio23N+RwD05ZlGuLulpcqa4zsAAoCSRhNm5QfX6JMQQkj0CjoAunDhAlatWoXjx497630AeE+wVPcTvGBKoDv3sfJXAA0AhR0ru4or9WjrKJbO9TP95U9BR0doc0f9kEoqgsnuQllHHVE801ucaGq3I0UduU1zCSGEhE7Qy+Dvu+8+FBQUoLGxEQqFAidPnsTu3bsxa9Ys7Ny5MwxDjH3BNEIcKAM0OlWFiwqTwNC1pD3fTwG0P/nJvsctm5gGAChvtgTdrDEWXWgyDXwQIYSQESHoAGjfvn148sknkZycDIFAAIFAgIULF2Ljxo249957wzHGmBdMbNFZA+SvCWKnH83KgUbWldzztwLMH4VEhLSODIdcLMSlY1MhFHAw2V3ewCuetdtcMFgD6yZNCIktNqebvgjGmKADILfbDbVaDQBITk5GbW0tAE9x9NmzZ0M7OtLLQBkgwDN1deNcT7G6WiZCSh/TZf50bgI6OUsLmViInARPD5xymgYD4NlZnoSWw8WDD2YpJCHDyM0zlDSasLe0GYcr9bA5qcwjVgRdAzRp0iQcPXoUBQUFmDt3LjZt2gSJRIJXXnkFo0aNCscYY16gXyoYYwPWAHWamZeAuy4pRIJCElABdKdlE9JhsrlwxeQMAJ6VYeUtFpQ1mzGbioBRb7ShKFXldwUeGZwz9UY43TwmZ+kgEQngcPHQWx1QS8XUf4lEVIvJjtN17d6gp83swL4LLUhUSMAACDkOapkIiSoJNLL4bBUykgUdAD366KMwmz3ZgCeffBLf//73sWjRIiQlJeG9994L+QDjQaCdoL09gDggQTHwH9v03ISgx5KbpMADl43xXs5PUgJoogxQB7ebobHdjkwddYcOhWaTHY1GT1B/oLwVWrkYje028LzndqVUhJl5CZCIKOAkw8fl5nG2oR11HX3UunO7GZrauzaNbjACaATGpKkDLjcg0SHoAGjZsmXe/x89ejTOnDmD1tZWJCQkBJVpIF0CzQB17tSeIPffAygcCjqmxCpaLLA53ZCJ6Rt5jd5KAVAI8DzDuW77rVkdbu8mvp3MdheO1xgwPUcHgaD/zxeeZwMeQ0h3DUYb5BIh1FKR9/xltDlxotrg7aYfqHMN7bA4XRiTqvb7e+h08xBT5jiqDLoPUHeJiTQ1MhSBVj+0BFD/E2rpGhmSVRI0mxx4e38lbltYMGzPHa0MFiesDjdNzwzRhWZTQCeZNrMD5xrbMS5d0+cxjDGcqDUgSSVFFgWnJAAlje0ob/aslBUJOUhEAogEApjsTm8GMljVrVYYLE5MztZ6m9ICQHWbBecbTEjXyjAmTQ0hBepRIeAAqK6uDs8//zx++9vfAgAWLlwIi6Vr53ChUIht27YhKysr9KOMcYGuLOis/xnOAEgg4HDbggJs+uws9l1owYRMDeZ36x59odmEBIUECXHWJbnZZEdOgP2VRjKeZ2izOGC2uyEQABlaeUg+vKvbLN6TT0DHt1ohFQm9GcnuGGM4VWdEo9GOVrMDySoJpCIKTuMZz3t+J1w8g4DzLAzRyMXeqdQGg83bKBYAXG4Gl9sNYOgFzu02F/ZfaEWySgqBALC7eLR2fHmtabOizezAlBwdVNKQ5B/6ZbK74HTx0MrFlB31I+B34K9//Sva2tq8l48ePYpbb73Vm/3573//iz/96U/4wx/+EPpRxrhAF8CY7V0NCodTUZoaV03JxEdHa/HWtxUYnaJCilqKEzUGPLvjPDK0Mvz6qokQxNEUaFOPAMhgcUItE8XUh4ybZyiu6mqoCQClTWZkaGXQKcTQyMSDmhJtarfjbLepr0CVNpog5LhedRbnGkzeWg2Xm+FcvQmTs7VBPz6JHReazag3dNXvNMLez9Gh5+YZGoy964cAzzZFB8paMTFTg1SNDE43D4vDDbU0dJ8frWYHzjW0w2RzAQCEQg5pahmK0lQ0DddNwGfSTz75BM8995zPdffdd5935de8efOwdu1aCoAGJbAIyO7yBECRqMO5YnIGTtUZcb7RhP/7tgL3LSnCPw9WAQDqDDacrW/H+AwN6o02HK3S43vjUmP6D63N7PDO6fM8w9FqPSQiASZlaYc9QA0Hf8EPADhdPCpbLKhs8Xyozh+VFNTvo8XhwolaQ1C9r7o719AOm8uNUclKiIQClDSaUNXqm0lqMNqQYpAiXSsb3JOQEaHFZEd1mxU8YxALBUhRS5GqlqLN4kR5c3Qv2nDzDMeqDZBLTN66N4EA0MrFSFXLkKqRDjqLabA4cbRK77OdkdvNUKu3otlkx7gMNVLV9LcBBBEAlZeXo6Cgq/7jsssug1LZlY4eO3Ys7RQ/SIGeDGxOz8S0LALpfYGAw5qL8rHhXydxqs6Iv3x1HrXdvmHtOteE0akq/GXHeTS02+HmGf6nYyl9LGLM8y0rTSNDvdEGh4uHw8XjUEUbFo1OHpGZIMYYqtusaDU7YLA64XD1XwjhdjOca2gPeJ84nmc4UWOE2z20nj+VLRY0GG1IVEr8rtIBgJO1Bgg4IFUT+Ad9z219SHRyunkcqdTD2KMpab3BBplYGNTeipHWveif54E2sxNtZifONbSjMEXl7cvWye5yw2BxApynWa1cLPQuiGGMwWh14UhVm0/w053DxeNYlQE5iU4UpapG5OdUKAUcADmdTjQ1NSE7OxsAsHXrVp/b29raIBDE7jf+cAr0z9Xa0YsiUsW3aRoZrpqaia1HanCixggAWDQ6GV+XNONIpR4fHqlBQ8fy0N3nm7B8YnpM/4E1tduRppH51BI4XTxazI6I7xlmd7lhsDrh5hkYA4QCT7+S7oWZ3bncPE7WGn2W9wai0WgPeI+0C82mXietwbI7+T6DH8AToJ6oNWAsz6CRiSATC/vNSDrdPI5V66G3OCERCTAxU4tEZXzVtY0EPM9wrLp38NMpVpoUMubZgNlkd6EgWYnGdjvqDTaY7a5ex4qEHMRCQUen6sAev6rVgjaL53NKIhQgUSmBMgYy18EK+BWPHTsWe/fuxfTp0/3e/vXXX2PMmDF+byP9CzwD1DEFFsGeKJdPTMN35a2obrMiTSPFjfNyUaO34kKzGZ+dagAAcJynY/WJWkNYdpHfU9KMgxVtuHVBPtQRbD7WbLKj2WTv9aHUYLSFPABijKHBaA94Wud8g8mnBgIApGIBZucn+kxZeZprOjwftrbeH66BONfQjiSlpFewy/OelVmtZgcEHDdgRinUeB44XWv0XtYqxEjXyCAQcGhqt8PmdCNLJ0eSSoLj1Qa0d7x+u5PH2fp2zC1IjOkAfqTheYaTtUa0meNnO5p6g63X33FPXQXcwTHZXN6/eYHA08coOyH2F3Z0F/CZ9LrrrsPjjz+OY8eO9brt6NGjePLJJ3H99deHdHDxItBGiJ0BkDSCvXhEAgHuWDQK03J0uG1BAUQCARYXpXhvH5WsxJJxqQA802KhZnW48Y8DlTheYwjL4wfD5WY4XWfsdX1TxxRgKF1oNuNEjQFl3WobXG4e9QYbTtUaUdGtUWW7zen3Q9Pu5HG0Sg+Xm4fd5UZpkwnflDSjuFI/6OAH8LwndT0KPnme4ViNAY1GO1xuNuzBjz8GixNn69txutaI5nY7TDYXzta3Y29Jizf46WS2u1AzDNueON180Fm3eMPzDNVtFuwtbemzsJgMDc8DZ+racbzaEPLPrmgWcAbo/vvvxyeffIKZM2fisssuw9ixYwEAZ8+exeeff4758+fj/vvvD9c4Y1qwNUDyCDcjzNTJcfelo72XZ+cn4P8drobZ4cINc3IhEwvxxelGHKsxoMVkR1IQe5ENZG9ps/fn8O2FVlwxOSOiNRt2Z+8Tu5v3dIoNVRFug9GGsiZPgNO5Esrh5lHdZoGrWz2NXCJEqlqG0qa+C0DbbS58V94Km9M96F4n/lQ0m5GplXnfixO1BjSP8BN7aZMJaRpZ2LpQ25xuHKnUw2x3IS9JgdGpKqo/6mB3uXGu3gSD1Qm7K/CpHTI0DUYbbC43pmbr4qL7esABkFgsxueff47Nmzfj3Xffxc6dOwEARUVF+M1vfoMHHngAYjHthTIYQU+BiaPrF1MqFuKRFeNgc7o7ts4AxqWrcaa+Ha98fQHLJqZjarZuyP1jeMbw5ZlG7+V6ow3lLRa/vWEird5oC0kAxPO9s0znGvwvIT9ZY4Q9lR8w8LDYQ18nYXG40dRuR6pGhhq91bu9xUjmcjNUtpoxOlUd0se1OT0/q/IWszeA9nRa5zEpSxP3QVCLyY6TtcaoyBrGI4PFie/KWlGQokRGx5RxIHieweJ0j6hVsEGNVCKR4JFHHsEjjzwSrvHEpYCnwCK4DH4g6T1W21wxOQPnG0wobTLjrztLMSFDg7WXDa1G7ESNAQ3tdsjFQoxJU+FotQH7LrREZQDUYrLD7nIPuSGf3ur0yfL0x82zQfXXCZWyZjOUUpHP9hYjXVWbFflJyn63njHanKjVW9FicmBmXkKff596iwOlTaY+a1g82zIIQh5wRQuz3QW7i/dbXN7UbkdlqxntNlfAv+8kfGxON07XGlHebEZ2ghwZWnmfGSGn27MgoaLVE9DrFGJkJcghEXo6a0dzf7SRE6rFsEAyQDxjXcvgozAA6ml8hga/uXoidp1rwqcnG3C6zgi70z2k+qUdpz3Zn4VFyRifrsbRagMOlLfiR7OyIYqyFYiMAVWtVoxOVQV0fFWrBelaWa+VSp37v40E7TYXDlf2vQR3JHK7GWr0VuQl+Q+yK1ssON/Y7v0bPlFjwMw8330RzXYXzjW0e7ey6U95swVKqQgZWs92HjzPUNlqQYvZgSnZWu/vR4PRhgSFZERMU/A8Q0WrBWXNJgDApCyttw8NzzOUNJlQ2RJ4V3AyfKwOd8cXWROSlFIkqSTQysVwuDzNG1vMDrSa7T7T6XqLE3pLV5AvFQuQk6CATiGGm2cQcBw0cnFUbAdCAdAI0T0dHG1TYH1JVcvww5k52FfaAqPNU1Q6KiWwgKAnk82Fkx1TQd8bm4pEpQRqmQjtNhdO1RrDstpsqKrbLMhPUgy4cW2D0dNIsqzZjHHpap/eNSOtQNZfTdRIV9lqQU6Cwvst1s0zmOwuVLdZei3F11ucKGs2Iz9JiXabC43tNlS1WYKqtzpZY0RZkxlqmRhmR9dKnUMVbZieq0NFiwWVLRakaWRR3/G61ezAmXqjz7Tr8WoD8pNdsDrcaDU7aKprBOB5z2fRYD6P7E4eJY0mn+sEAiA/STno80GoUAAUBQJp3NVZ/8NxgGSEdVjOTlDgVJ0R1W2DD4AqWj2FvWlqqXeJ+ez8RHx5phHFVfqoDIBcboZava3X1g2+x/Demh6Hi8eJWgPmdfTrabc5e+2OToaf3cmjzmiDSMChus0CvcXZb9b2QpMZ5S3mIRWZWxzuXhvFmmwu7Clp9j5utHa8NttdaDE50Gy2e/fA6o4xeIv6SXzi+a5FPZE0ss6kMSqQKTBvE0SxcMQVSWYneNL51W2DX1Zc3pEi7z4VMTbNUytRFsVt7ytazeD7mRIqbTL7ZE14HjjTUUPTHMCUCRkep2uNOF5tQJu5/+CnUyhX2PX3uGfqjQMGyQ4Xjxq9FS53eE84FocLx6sN2FfagnMN7X6DH0KiCWWAokAgFROR3AZjqHI6mmtVtQ1+nr+8o89NfnJXNqWz+LlGb4XDxUdlPYTdyeN0vRF5ScpeqyMMFieq/fxMWk0ONBhtI276iww/l5vhSGUbZuYn+C24tzrcOFLVBovdjXNCDhlaGbITFENaqdNisqOqzQpnR0Dl5hlcbkbL1cmIE9Bfwdq1awN+wM2bNw96MPGKBTEFNlLqf7rrngFijA0qg1XR7AkU8rtlgBIUYmhkIhhtLlS1WVAY4fnkvtTpbajT25CglHgLWXme4WRd35uCnqlvh5NqI0gALA5PP6EZuQmeuhqLA26eB8ChtuPLAeAp6K5utaK61YoEpQQZWhlS1NKANy12uHgUV/W9DQUhI01AAdCRI0d8Lh8+fBgul8vbDPHcuXMQCoWYOXNm6EcYBwLLAEXvEviBZGhlEHIcrE5P0WOwjRENVidaLQ5wAHITuzJAHMchP0mJYzUGlDebozYA6tRmduBwRRum5yagstXcbz8eCn5IMEw2F74+3xRwBqbN7ECb2QGBwJOhHZWi8lmVw5hn41qlVIiCZCVcvCfT1LNjNiEjWUAB0FdffeX9/82bN0OtVuONN95AQkICAM9GqLfccgsWLVoUnlHGuEA+tEbSEvieREIB0rWeBnlVbVZvAKS3OPD6nnIsHpOCmXkJfd6/c5uHdK2s1+vPT+4IgEbIMtp2mwsHy1u9NV2EhMpgpp943tOEscFoR2GqEqlqGQQccKrO6N12os3iAGOg4IfEnKAngv/4xz/is88+8wY/AJCQkICnnnoKl19+OdatWxfSAcaHwKfAIr0NxmBlJ8hRo7eius2CaTk6AMBXZ5twss4Ii9M9QADUe/qrU37HCquylugthO6p5+oeQiLN5nTjZI0RZ4TtUEtFPn1c4mnzURJfgi4oMRqNaGrqvQllU1MT2ttjpwPscAooA+Tq3Ah15NUAAV2F0N1Xgh2pagMA1Oqt/bYC6CyAzvOznLwzKGow2GjJOCFD5HYzn+CHkFgW9Nl01apVuOWWW7B161ZUV1ejuroaH3zwAW677TZcc8014RhjzHMF0DnXOoJrgIDeS+EbjDbUdjSRs7v4frvklveTAdLIxUhUSsDQ1SuIEEIIGUjQU2AvvfQSHnzwQdxwww1wOj3fFEQiEW677TY888wzIR9gPHAG0J8jWnaCH6ycjuLlhnYb7C43iqv0PrdXt1m8DQ6701scMFid4DggJ1Hu97ELkpRoNTtQ3mzBuHRNyMdOCCEk9gSdAVIoFPjrX/+KlpYWHDlyBEeOHEFrayv++te/QqmMvk0pR4LAAqCRuwweADQyEbRyMRjz7OnVGQCJhZ6VJzV6/00ST9d5plUztfI+Nxbt7A1UPoLqgAghhETWoM+mdXV1qKurQ1FREZRKZUC9bEhvLjcfUNdYbwA0AhshAp4l61dPywQAbCuu8e4Ns2h0CoC+A6Avz3o2QJ2V33eRdOfUWGmTiX4PCSGEBCToAKilpQVLlizBmDFj8D//8z+oq6sDANx22220AmwQnO7ATtjeZfCSkRkAAcDC0cmYNyoRPPOse8tNVHg3c6zxs01GaZMJZc1miAQcLi5K6fNxC1NUEAk4tFmcqDPY+jyOEEII6RR0APTAAw9ALBajsrISCkXXqpxrr70W27dvD+ng4oEjwP15ujJAI3MKDPBkgX48N8+7eeOMXB2ydJ66ngajvddU4I7TnuzPnIJEaOTiPh9XIhJgTMe+YCdrjeEYOiGEkBgT9Nn0s88+w+9//3tkZ2f7XF9UVISKioqQDSxedJ706402vPr1BW/zsZ5Gcifo7mRiIdZdNgbXzc7B5RPSkaAQQy4Wws0Y6ru9dr3FgUMVnmXyS8alDvi4EzM9xc8n6wzhGTghhJCYEnQAZDabfTI/nVpbWyGVBrfFAekKgL4+34T9Za344HC13+NG+iqw7hIUEiwdnwaJSACO47xL5LtPg315phFuxlCUqvLZAb4vnQHQuXpTQEXlhBASjAajzftFlMSGoAOgRYsW4c033/Re5jgOPM9j06ZNuPTSS0M6uHjgdHlqgOwdAc6xagMsjt4t5zsbIY70DJA/mR3TYJ2F0CabCzvOeKa/LpuQFtBjZOnk0MrFcLh5b4E1IYQMldnuwut7yvCrbSfwq20nerXwICNX0H2ANm3ahCVLluDgwYNwOBx46KGHcPLkSbS2tmLPnj3hGGNM66wB6myG6OIZDlfqsXB0svcYxli3RogjtwaoL9k63wzQZ6frYXfxyEmQe7fNGAjHcZiYqcHe0hacrDVifAb1AyKEDM2Ryja8tb8SBqun553B6sTzX5UgUyeDxe6G3cVDpxAjRS3F8onp3lpEMjIEfTadNGkSzp07h4ULF2LlypUwm8245pprcOTIERQWFoZjjDGtc7qm+7TN/gstPsc43Lx3u4xYzABldUyBlbWYcabe6C1+vmpqJgQc199dfUzsCHpO1lIdECFk8NptTryy+wJe2FkKg9WJNI0U6y4bg+UT08FxQK3eBr3VCavTjTqDDceqDfjj5+ewr8dnN4luQWeAAECr1eJXv/pVqMcSlzoDH1e35fBn6tuhtzggEgggFHDeLBEHQDqCV4H1JTdRAY1MBKPNhT98dg4Agsr+dOrM+lS1WVFvtCFdIwv1UAkhMe5cQzte3n3B24F+2YR0XDU1ExKRAOMzNFg4Ohn1Rht0cjGkYgHazE7sOteEQ5Vt+Ps3ZWgw2HDl1EwIBV1f3r4ra8XHx2oxJVuLq6dlQSwM/+c4YwxlLWYYrS7kJSmQoJCE/TlHmqADoFGjRuHiiy/GSy+95FP03NzcjDlz5uDChQshHWCs85cBYgD+8Pk5NHT8kd27pAiAZyNULoiMyEghEwvx0PJx2HakBgc7Vn6tnJYV9GvVyMWYkqXFsRoD/lVcizsWjwrHcAkhMeBQRRu2n6xHklKCLJ0cHAfoLU7sPt8EngHpWhl+sqAA+cm+izDStTJvKw8AyNDKMS5Dja2Ha7D9ZD0+OV6Hk3VG/HBmNoQCDl+fb8Y3Jc0A4M0W3bawwO/ehqHAGMOekhZ8cabBZ/PpJKUEV0/PwryCxJg8jwxG0AFQeXk5RCIRFi1ahH/9619IT08HALjdbloGPwiOjiJoZ0c76MIUJUqbzKjvaOjXZnGipMFT1BsLK8D6kq6R4WcXF6KmzQqjzTnoGp6rp2fhWI0B35W3YsWkdO8eZIQQ0mlfaQte21sGxoCyZrP3i1enuQWJuGleXsAlBwKOw//OzEZOghxv7a9EWbMZmz49672dA7CoKBnFVXrUGWx4+r9ncNP8PCwoTIaL92wGnaiUDDkzxBjDewer8EVHGYFYyCFFLUWdwYYWswN//6YM+8tacPO8fCQqKSMUdADEcRy2b9+OBx98EDNnzsS2bdswe/bscIwtLrh43ymwi8ekIEsnh0IiwvnGdpQ2mXGmwbMfljSGA6BOWQlyZMH/pqeByE1UYE5+Ir4rb8WHxTW493tFIRwdIWQkMFid2FPS7KnTcbihlAqRqJSA54Emkx27zzWBAZg/KglZOjnqDFZwHAeZWIDRqSrMzE0YVJZk7qgkFKWp8c7+SlxoNkEqEkIjF2HV9CyMS9fgmukubNlbjuJqPV7fU45vzjejus0Kq9MNjgNSVFIUpqgwMVODyVlaKKWBn6LdPMP/O1ztDX5WTsvE98amQikVwe5044szjfj4aC1O1Bjx649P4taFBZiarQv6NcaSoAMgxhhUKhW2bt2K9evX4+KLL8Yrr7yCyy67LBzji3k9p8AUEhFunp8PAPiouAalTWacrfcEQCO5C/RwWjktEwcrWnGs2oALzSaMSlb5PU5vceBgRRtGJSuRn6wMquCaEBKdLA4XNn16Bg1Ge7/HXTImBTfMzQ35332iUoK7vzfa720qmQg/v7QQHx+txcfH6nC+o2WHUMDBzTM0ttvR2G7HvgstUElFuOuSQhSlqcEzhvMNJhyubMOJWgOcbgaZSACpWAipSACnm0dVq9VbL3rTvDxcPKZr+yCpWIgrJmdgRq4Of/umDBUtFvzlyxKMTlFBIACSVVJcMz0LujirExpUBqjTxo0bMXHiRNx+++24/vrrQzqweNB9I9TOZfCdu6MD8DYANNk9fYFieQoslNI0MszMS8CB8jYcrzb4DYAYY3hxVylKmzw7yKtlItw8Lw/Tc/vedJUQEj2cbh6n6oww2VxwuHika2UYk6bGa3vK0WC0I0EhxvzCJMjFQpjsLrSZPUXNOoUYeYlKzM4fXJZnqAQch5XTslCYokJlqwXjMzTIS1Kg3eZCdZsFp+vacaSyDQ3tdvzx83O4ZGwKjlYZ0GTqP6ADAIVEiB/OzMaiPvZOzNDK8cjycfh/h6qx40wjSpo8Adi5BhOOVRtwy4L8uMoKDSoD1N2Pf/xjFBYWYtWqVSEbVLzovhFqZwao+xxwXpJv/UosLoEPlzGpahwob8OFZrPf2w9WtKG0ybPRqlgoQLvNhTe/rcDETC0klGkjJGq5eYZ9F1rwr6O1aDU7fG5TSISwONwQCTj8/JLRKEgOT6FxKEzK0mJSltZ7WSsXQyvXYmKmFldOzcDfvinDkUq9d0pLLhZieq4O03N00CkksLvcsDl52F1ucOCQm6hAqkY6YEZLLBTg+jm5mF+YhKZ2O3jG8OnJBlS2erJCK6dm4vtTMuKiUDroAIjne28zMH/+fBw9ehRnzpwJyaDiRfeNUDtrgETdMkA6udi7PBygACgYo1I8H3xlzWbwjPl8KDjdvHfLkf+ZnIEVk9Lx2Ecn0GxyYE9JMy4NYO8xQsjwc7l5PL+zBCdqPJsea+Vi5CTKIRIIcK6hHRaHp2HsjXNzozr4GYhUJMSdiwvxr2O1ONfQjnmjkjA3PzGkdaD5SUrvSrQZuQn4oKN+6KOjtTDanLh+di4EgsCDIIvDhVazA5k6+YgpJxhUHyB/0tLSkJYW2LYFxKP70ndvBkjQlX3gOA55SUocr/E09ovFLtDhkpUgh0QogMXhRoPRhgxtV2H1l2ca0WxyQCcXY9mENIiFAlw+IR3vfFeJT0/VY/GYFJ8eHoSQyON5hr/vKcOJGiMkIgFWTs3EpWNTvRlbp5vHiY7PyliYyhYIOFw9LWtYnkssFOC62blIVcvwj+8q8dXZJpQ1m3H5hHTMzEvo9/PQ6nDj89MN+PxUA6xON5KUEswvTIJaKoKLZ8hJUGBchjoqg6KAAqAZM2Zgx44dSEhIwPTp0/tNjR0+fDhkg4t13QOgrhog3yAnL0nRLQCiDFCgRAIB8pMVONdgwoUmszcAcvE8/nuiHoBnyXznN6oFo5Pw8bFaNJscOFDeinmjkiI29nhisrvwt68vQKeQ4JrpWVBIhdh5tgmHKtrAMwaJSIBlE9J9pgr6orc4IBRwUMvEwzByMpwsDhfePVCFA+VtEAo4/Pziwl6/E2KhICYCn0j63rhUqGUivL6nHOUtFrzy9QVoD4oxMy8Bs/ISkJ+khKSj6PpCkxn7LrTgQHkr7C7PuUzIcWgxO/DJsTqfx01WSVCUqobR6oTdxSNTJ8PkLC14xgL62w6XgAKglStXepseXn311eEcT1zp3AgV6AqGuk+BAUBetz42VAQdnFHJKpxrMKG0yYQFHXurnao1wmR3QS0TYX63IEcqEmLJuFRsK67Ff47XYVZ+AkQCyriF2weHqnGi1jOdcaSyDRq5GHUdPbA6nW8wYe1lY/rdZ6m0yYQ/fn4OMpEAT141CSpZyJLbJEIYY2gy2XGs2oB/H69Du80FDsBtCwoietKMdbPzEzE2TY2d55rw1dlGGKxOfHmmEV+eaYSQ45CskqDZ5IC7Wz1wulaGlVMzMSVbiyOVehyt1nu3bzpZa0SzyYFmU9c2ISVNJuw+34wz9e34+5rItdEJ6FNiw4YNfv+fDE1nDRBjrKsGqEeqMa9bt1DKAAWnsw6oeyH0d+WtADx/5D3Tut8bl4rPTzWg1mDDZycb8D+TM4ZvsHHofGM7vu7okJuulaHeYIPZ4YZKKsL3p2QgUSnBNyXNOFZtwPNfleCR5eOQqevdI6reYMNfviyBw8XD4eLxYXENbpqXN9wvh4TQ3tJmbD1cA33HJqSAp1nqdbNzKPgZBhq5GFdNzcT/TErHyTojDpS34lStEUabCw3tntVoapkIk7O0WDg6GUWpKu/M0LxRST4ZdLvLjSOVerSaHdAqxBALBKjWW9BotEc8005fkyKoM+vjZgydsXTPKbAEhRhqmQjtNhekVAMUlFEdRZA1eitsHY3GjlTqAXg6vfakkIhw7ewcvLanHP86WosZuQk+Le9J6Lh4Hm99WwkAWDg6GT+el4udZ5tgtDlx+YR0qDoawE3M1OCPn53DhWYz/vTFOay9bIxPPZfB6sSzO87BZHchTS1FQ7unyd2i0cm9tjAIBGMMDhcPkVBAdWAR4HTz+Md3ldh93hMYCwUc8pMUmFuQhMVjkikrO8xEQgGmZuswNVsHxhhazQ7vPouJSklAK8WkImGvQGcOEpGpk2NC5uA6/odKQAFQQkLg/RJaW1uHNKB44m8j1J5TYBzHYVKmFvsutCCDTsZB0SkkSFRK0Gp2oLzFjHabC3YXj2SVxBsc9TR/VBL2l7XiZK0RW/aW46FlY4NaCUEC8/HROtTorVBJRfjfGdkQCQRYOr73IgqpSIh7vjcamz49izqDDb/ffhb3LylCfrISNqcbf95xHs0mB1LVUjy8fBzeO1iF/WWtePu7SqxfPi7g967OYMWrX5ehps0KN2NIUIix4fsTaSothBhjcLh5SEXCXtcfrtTj27IWnK4zwubkwQG4amomlk1Mp7YUUYLjOCSppEhSSQc+eIQI6K/72WefDfMw4pO/jVDFfr7hrJ6fh6unZcbUL95wGZWsRKvZgW9KmmG0etoJzMnvezNAjuNw87w8PP6vkyhpMuHVby7gtoUF9M0zhL4ra8W/j3uKJK+fkzNgkKGWifHQsrH4847zKG+x4OntZzCnIBGtZgcqWy1Qy0S4f2kRNHIxfjgzG0er9ShrNuNv35Th1gX5EPnZX8nmdKPJZEeGRoYmkx1/+OwcDN2mW9osTnx2uh7XTM8O7YuPUxeaTXh9TzlazQ7c2a2A2eHi8c53ld7NQgFP1nvNRfmYmElTXSS8AgqAVq9eHZYn3717N5555hkcOnQIdXV1+PDDD32KrBlj2LBhA1599VXo9XosWLAAL774IoqK+t7faePGjdi6dSvOnDkDuVyOiy66CL///e8xduzYsLyGofBuhNqRARJynN9vrCKhgIKfQRqTpsbBijZ8e6ErMznHz/RXd0kqKX6ysAAv7b6AA+VtcLoZfrp41JA3KhypDFYnvilpRqPRhkydHBlaGdw8g9PNMD5DHdSqq/JmM17fWwYAWDYxDXMLAqsBUMvEePDysXhpVylO1Bqxt9RTUCkVCXDf94qQqvZkR3UKCW65qACv7L6A78pbYbK7MC5d7S18z9LJUdpkxpdnGmF1uiERCSDkOFidbmQnyHHnxYWo6Fj9suN0Iy4fnx62LBBjDP89UY8TtQZcMz0bo1P9b9kykrncPD4+Vof/nKjzFsU+/1UJbl80Ci6ex/YT9ahqs4LjgMvHp2F2fiJykxRRuWSaxJ4h/WXbbDY4HL6dODWawOf0zGYzpk6diltvvRXXXHNNr9s3bdqE5557Dm+88QYKCgrw2GOPYdmyZTh16hRkMv/TQbt27cJdd92F2bNnw+Vy4Ze//CUuv/xynDp1CkpldDXG8k6B8f5XgJGhW1zkWf11sKIV5xtMGJOmRnbCwDvET89NwF2XFOKvO0tRXKXH+4eqccOc3HAPN6q43J5v53tKW+Dmmd9jchMV+NX/jA+oXsbqcOOl3aVwuhmmZGnxgyCzKzKxEPctKcKFZjN2nWvChSYzrp+T06vWZ2ZeAu753mj8dVcpTtUZcarO6PfxRAIOjo7luzkJcqy9bAzUMjFS1VL894QcVW1WfHaqHtfMCH0WiOcZ3v6uErvONQEANn16pqPoNCNmplyrWi34+54yVLdZAXgyr06ex5FKPV7cVeo9TiUV4Y5FoyJeD0LiD8d67m0xALPZjIcffhj//Oc/0dLS0ut2t9s9uIFwnE8GiDGGzMxMrFu3Dg8++CAAwGAwIC0tDVu2bMF1110X0OM2NTUhNTUVu3btwuLFi/0eY7fbYbd37bNiNBqRk5MDg8EQVEAXrB2nG8AYUKu34vF/nYRKKsKz104L2/PFO5vTDXGQxa3FVXo8/1UJAODhZWNR1M9S7FjCmKfpXGfmbFSyEhMyNKgz2tBotEEsFKBGb4XdxeO62Tl+63d6+vs3Zdh3oQXJKgke//4EKCThra+50GzCF6caIRRwUMlEaDM7UKu3QikV4bIJaZiWrUOdwYbqNgsmZ2t9xnOksg0v7CyFVCTA09dMDmlvocpWC7YdqcGxGgM4AGPT1TjTseHxwtHJWD0/r9cULWNsRGxN4OJ5HK0yYF9pC47V6MEzT4Bz07w8zMxLgIvn8do35fiuvBXJKgnmjUrCJWNS4m4TToKwFUEbjUZotdqAzt9BfwI99NBD+Oqrr/Diiy/ipptuwgsvvICamhq8/PLLePrppwc96J7KyspQX1+PpUuXeq/TarWYO3cu9u3bF3AAZDB4mggmJvY97bFx40b8+te/HtqAg+TmmTcl7O0BFCPf/KLVYNoITMvRYeHoZHxT0owt+8qx4fsT46Ioc1txLb690AoBB/z8ktGYlqPrdczOs414a38lthXXYFZeQp8nMZ4xfHuhBfsutIDjgNsWFoQ9+AE8faDuWNz/tFJWghxZCb2X1k/L0SE3UYHKVgv2lLRg+aT0AZ+v3eaEUirqc/qmVm/FPw5U4nSdJ9gRCjjcvqgAM3MTsKe0BW/sK8c3Jc3I1Mlw+QTP85ntLry1vwIna424cU4u5o5KQo3eijf2liM7QY6b5vUOliLlXEM73tpfgVp9Vx+n6bk63DQ3Dxq5J4AUCQS4fVEBfjAjK+BVRISES9CfQh9//DHefPNNXHLJJbjllluwaNEijB49Gnl5eXj77bdx4403hmRg9fWebr09t9dIS0vz3jYQnudx//33Y8GCBZg0aVKfx61fvx5r1671Xu7MAA2XvrpAk+jwo1nZOFFjQIPRjk+O1YZlSiRanG9sx3+O13u7j988L99v8AMAi8ekYG9pCy40m/HewSr8dHGhz+1GqxP/920FTtd7VvYAwBWTMlCUGv1ZNI7jsLgoGW/tr8SBitZ+A6CqVgs+Kq5FcbUes/MTcMeiUT4ndp5n+PRUPT4qroWLZxBwwKy8RCyflI7cjkanC0cnw+pw472DVXj/UDUsdjcUUiG+ONWIVounzODVb8pQ2tF91+p040KzGWPS1BHppdJisqPBaEeL2Y46gw1VrRac7shiqaQiLBydjPmFScjy07epczURIZEWdADU2tqKUaNGAfDU+3Que1+4cCHuvPPO0I5uiO666y6cOHEC33zzTb/HSaVSb6frSOirCzSJDgqJCDfMzcVfd5biq7NNuGJKRq+lvCMdYwxv76/Ezo6aFI4DVk3LwsKOGip/BByHH8/Lw2/+fQoHytuwcpqnPwjgqfd5dsd5VLZaAAASkQAzcxPw/akjp7nkzLwEvPNdJSpaLGgw2pCm8a075BnD1sM12H6y6wvZgfI2TMxswcKOzuN1Bite31PubcY5KVODH8/LQ7KfAGDp+FTU6q34uqQZnxzv2kogTS3F6FQV9pS24Muznp3BOzdJfvdAFSZlavst1LY73dBbnUhVS4eccanVW/HB4WocrTb4vX1xUTKumZHt7eNESDQL+rd01KhRKCsrQ25uLsaNG4d//vOfmDNnDj7++GPodLqQDSw93fONq6GhARkZXR+aDQ0NmDZt2oD3v/vuu/HJJ59g9+7dyM6O7m/snX2AhpoBEgk5n55CJHSm5ei8LeAPVbThosK+A4ORhjGGfxyows5zTeA4YGFhMpZPSu91wvcnN1GByZlaHKsx4OtzTfjhrBw43Txe2FniXaJ+1yWjUZCsHHGNBdUyMcala3CqzoiDFW24oltncBfPY8vecm+N1Oz8BOjkEnx+ugHvHqhEgkKMU3VG7DjdCBfPIBcLce3sHCwoTOq3BcONc3ORqpGiVm+Dw80jVS3FFZMzIBUJkK6V4aPiWszKT8CP5+bhd/89jVq9De8fqsItCwr8PmZTux2bvziHpnY7EpUSTMnSYkq2FuPSNQNO5fasO/rkWC0+OloLxgABB6RpZEhSSpCqkSFbJ8foVJXfTt2ERKugA6BbbrkFR48excUXX4xHHnkEV155JZ5//nk4nU5s3rw5ZAMrKChAeno6duzY4Q14jEYj9u/f32+miTGGe+65Bx9++CF27tyJggL/HwzRJFQ1QCqpCHqLc+ADSdAEHIdFRSn48EgNvj7fHDMBkMvN45+HqvHlGU9mYc1F+VgQ5GtbPCYFx2oM2FPagqunZ+GfB6twpr7ds0R9SRHyk6Jr9WUwZucn4FTHVgBXTM6A3uLAd+Wt2Fvaguo2KwQcsLrjZ8bzDGXNZpQ0mfCnL857H2NSpgY3z89HonLgQl+RUIAVk/xnyVZMysDS8WneL0o3z8vH09vPYE9pCwpTVVhclOI9ljGG6jYrnt1x3tvfqNXswM5zTdh5rgkSoQAFyUoUJCsxJVvba5+18mYz/rqzFGkaKX48Lw+HKtqwrbgWgKeu5wfTs6lLOhnxgg6AHnjgAe//L126FGfOnMGhQ4cwevRoTJkyJajHMplMKCkp8V4uKytDcXExEhMTkZubi/vvvx9PPfUUioqKvMvgMzMzfXoFLVmyBKtWrcLdd98NwDPt9c477+Cjjz6CWq321gtptVrI5dH57cQZogyQQkIBUDgtKEzCtuIanG80od5gGxEnAIvDhZO1RpysNcJgdeLaWTnecVe0mPHannLU6D3LlG+alxd08AMAk7O00MnF0FudeHt/V1O7Oy8uHNHBD+Bph/DWt5WobrN2ZHxavDV7UpEAP7u4EJM7mvoJBBxuXZiPp/59GnYnj8lZWiwYnYRpObqQFft2/4wYnarCFZMz8O/jdfi/bysg4DhYHC58V9aKBqMdVqdnRW6WTo67Li1EvcGGo9UGHK82oNXiwNmGdpxtaMf2k/UYl67GymmZKExRoU5vw5++OAezw41WiwMb/nXS+5p/MCOrzwCNkJFmyBO1eXl5yMsb3MaDBw8exKWXXuq93FmIvHr1amzZsgUPPfQQzGYz7rjjDuj1eixcuBDbt2/36QFUWlqK5uauLqIvvvgiAOCSSy7xea7XX38da9asGdQ4wy0UfYCEQo72CgszncIzhXC02oCvS5rww5nDVygfrFO1Ruw614Sj1XrvyQsAylvMuOfS0ThWY8B/jtf1WqY8GEIBh4VFyfjkWJ03+FkyLjUmNq1USUWYkKnB8RqD97WNSlZi/qgkzMpP6LU8PlUtw+9WTQYHQDkMdTBXT8tEu82J3eebsWVvea/bJ2RocMfiUVBJRUhVyzClY0+nWoMNpU0mnG8w4UB5K87Ut+PM9rNQSUVgjMHscCM/SQG5ROhdtXbllAwKfkhMCboPEAAcOHAAX331FRobG8HzvM9toZwGi5Rg+ggMlptn+Kpj2mHXuSb837cVmJatw93fGz2ox1PJREjXyFDSaArlMEkPnX2B1DIR/vi/U6Oyad2JGgOe3dE1BZOulWFyphZnG9q9RcmdZuUl4Ma5uUPuc9NisuORrcfBOp7v8SsmxEy7gGPVevzlyxJkaGX44awcTMrURNXybZ5n+Ns3ZfiuvBU5CXIsHpOCsWlqJKkkARXrt5js+NfRWhwob4OjYzo+O0GOBy8fC6VEiEOVbXC5GeYW9L2FDCHBGpF9gH73u9/h0UcfxdixY5GWlubzB0F/HIPjCsEqMJlYOOKKTEeiyVlayMVCtNtcqGqzIC8KpnjsLjdMNheSVFLwjOGDw9UAgOk5Olw1NRM5HUutLQ4X/vJlCc43mqCSivDjubmYld//tiCBSlJJMb8wCcVVevxkYUHMBD8AMCVbh80/mgqlRBSVAa+go5/Qj2ZlQysXB/05nKSS4pYFBbhpXh7Kms2o0VsxKy/Ru5JrVl5ofkcIiTZBB0B//vOf8dprr0XtdNJIFIoaIJk4uA7HZHCEAg6jU1U4XmPA+UZTRAOgzgaDHxyugdHqxPVzcqGUClHVZoVcLMTqi/J9liMrJCKsvWwMjlUbMCZNFdLuxgBw64ICuHkWk7+Hof5ZhRrHcUPupiwSClCUpo6bbueEBB0ACQQCLFiwIBxjiVudNUDioWSARJQBGi5FHQHQuYb2gLaACIWj1Xr8+1gdlo5P8+6E/urXF3C+25TnO99VQt7R7XrZxDS/vVjEQsGga30CQb+DhJCRYlCrwF544QU8++yzYRhOfOrs3SMSDD4DJJcIaQflYdK5ZPhcgwmMMZypb8fre8qhU4iRl6TAotEpyE0aeMPVQDWb7Pjb12WwOt145esLOF5jwIlaA9ptLkhFAnx/SgasTjf+c7weVqcbaplo2AIzQggZqYIOgB588EFcccUVKCwsxIQJEyAW+6aGt27dGrLBxTKLwwUXz0MkEHj7AA01A+QOvp6dDEJ+kgJiIQeT3YU6gw0fHK5Gq8WBVosDF5rNOFFrxO+unhSSmjieZ/j7N57gJ0EhRpvFiX0XPJsQ5yTI8fNLRiNF7ekqLBcL8enJBlw3O2dQ+54RQkg8CToAuvfee/HVV1/h0ksvRVJS311NSd/Mdhf+589fgwHY8P0JcHYsUxYNoQZIKhbA7uIHPpAMmUgoQGGKCmfqPT1UylssEAk820K8s78STe12VLdZvcXHQ7H9ZD3ON5ogEwvw0LJxqG6z4L2DVRiXrsH1c3J8VvmsmJSB5RPT6W+SkDAQCACePmJjStAB0BtvvIEPPvgAV1xxRTjGExd2nGlEVZun+ZzDzQ95FZhA4GnK1r3fSyjQH3zfilI9AdDeUk82Zk5BIhaOTsbRKj2OVOlxuLJtyAGQ083jvyc8jTyvm52LFLUUKWoppuf2XcNDwQ8hoZekkmB8hgZGqxNnG9phd/LgOM/qW54xuNwM7hB//pLwCzoASkxMRGFh4cAHkj7951jXRoc2J9+1CmyQNUAykRAcx0E4xJNfokqCVpPDezldI0dtR5dg4stTB9T1Pi4d56m5mZ6r6wiA9Fg5LWtIz3Gy1gir0w2dXIyLCod/x29C4p1QwGFMutq7q71MLESCUgKr0w1Vt7YITjePY9UGtJkd/T3csBMJOaikIlidbtid9G22p6DPuE888QQ2bNgAi8Uy8MGkF7Pdha86dnQGAJvTPeTd4KUd9R5DqKGGViFGSo8dqpNVEiikVEviz6hkpTfgLEpVeYuep2TrIOCAGr0VDUZbQI91vrEdr+y+gI+Ka3CsWg++45vkd2WejTZn5SdQgTshw0wjF2PuqERv8NNJLBRAIxP79IQSCwWYnqPrd3ucdK1sSL3egiEWCTA2XY1FRSmYlZ+IRUUpuGh0EhIC2I8ungSdAXruuedQWlqKtLQ05Ofn9yqCPnz4cMgGF4u+PNPoU6tjdbq9U1eD7QMk69gCYygZoESlBDqF73upkYuhlYthsbuDeqwUtRRN7fZBj2UkkIqFGJOmwun6dlw+oWvFlUoqwth0NU7XteNIpR7LJ6X3+zh2lxuv7i5Dq6Xrm+O0bB1uX1SAo9V6AMCcEDUrJIQEJidRgaJUVVCNLwUCDpOytEhUSnCuob1rda+Qw8RMLVLUUticbpytbw/r56NMLMSMPB0UEt/Tu0Iiwsy8BNTqrThb305TdhhEANR9I1ISvP+eqPO5bHO6h1wD1Nn7ZSg9WJKVUqikIgiFHNxuBrlECJlYCK1cjDp9YJkMAFBIhZiYqcHe0hY4wliUzXFApBe9/WTRKDQYbb120p6Rk+AJgKraBgyAPjvZgFaLAwkKMcala3CgvBXF1Xq8uKsUdhePZJUEBcmR7zZNSCwSCjjwjIExTwZdLhahIFk5pI2OM3VyJColaDE7IBcLoZKKvJ3RZWIhpuboUNliwfnG9pB/himlIkzP1fW7CjRTJ4dWLsbxGgNMNldoBzDCBBUAuVwucByHW2+9FdnZ2eEaU8yyOFz4smP/L5VUBJPdFZoaoI5fdo7jIBRwAUX2EpHAG6CIhBw0chE4joNWLkaryQGt3JMN6vxvXxKUEpjsLjg7HitLJ4dIKEBBshJn69sH9XoGIhRwyElUoLzZHJbHD5S2I0PW0/RcHd7+rhKlTWbUG21I1/j/MG01O/Dfk54i5x/OzMGcgkTkJSnw7oEqnKg1AvBsQ0CFzYSEnk4hxrQcHYQCDi6eQSTgQva3JhMLe02ddZebpIBSKkSt3oY0rRTJSincjMHh4mGyu9Buc6K6zerNInWnVYghF3v6vgkFHIQCQCIUIlEl8dv81B+lVIQ5+Ylw8jzEAgGMNidO1BhhcwaX7R/pgjrjikQiPPPMM3C54jtqHKyvzjTB5uSRkyjH1BzPTtmeKbCh9QHqHu0HmrItTFUhuaN/TJJS6v3D1/UIfFRSUb+ZpewEOUZ1ZCg4Dt5vTlk6ORQSIRQSYcg3vEtRS5GdIEe0xgU6hQRTsj3v7yfHavs87oPD1XC4eIxOUWF2vmdl15JxqZjY7edF01+EhF6iSoLpuQkQCQXgOA7ijv8OpySVFJOztUhVyyAQeMaglIqQppFhdKoas/ITvZ/tQiGHvCQFLhqdhNn5iZiUpcWETA3GpqsxOlWN3CRFwMFPJ4GAg1QkhEDg2UZl7qhEb0+xeBH0FNj3vvc97Nq1C/n5+WEYTmzTyEW4qDAJ03N13l3bPUXQQ+sD1FkDBHjqgJwDHM9xQIpKiiSlBPssDiSpugrjPPsJmaHtqAfiOE92qM3c+1GFAg7JKikEHFDdZoVCIvT2pREIOEzN0Xm+qQg4nznxocrUySETC5GkkqI5SmuNrpqaiWPVBuwva8UVkzOQofX9NljaZML+jiLn62bneD98OY7DLRflY9OnZ5GskiInse9vkYSQvsklQmhkYiilnmwJzxhkYiHUMhFUUlHUZ1ZVUhFm5Seg0WhHhk42pL0iAyEWCjAlW4uTtUbUGwIve+hOJORC9jk/HIIOgFasWIFHHnkEx48fx8yZM6FU+tYnXHXVVSEbXKxZVJSCRUUpYIzhwfePAvCtARp0BkjUPQM08PE6hdg7Jz0mTY3EbisDtHIxREIO6m7fJrRysd8AKEUt9WaHxqSpenWiVnZ7DKVUBINloNBsYDKxEAkdwVmmVha1AVB+khLTsnUortbj46N1uH1RAdyMQSQQgGcM7x2oAgBcVJiE/B41PjqFBL8NUSdpQuKRRi7GzLyEEb83nUwsDOm2OgPhOA4TMzWQiASoarUEXKOkkYtRlKqCTiGGwepEU7sdAgEHidAzvdbUbofLzSAWCSAVCaJmqi3oAOjnP/85AGDz5s29buM4Dm53dLywaMZxHFRSz0m8ew3QYPYCE4sEPtNegawES1F11aT0nKcWCriO6aWux9H0UQfUvVAwSSUF6+evRS4WwjBgbmpg6VqZd2zJKinEIoG3/ijaXDU1E8XVenxX3orDlW1w8wyTs7UoSFbiQrMZUpEA10z33yuIgh9CAicRebYUYsyT+Zmaox3xwU+kcByHMWlq5CYqvO08+loJnKSSICdRgSSlpKuMQiHpmEnowvMMTp736VzPR8EqtKADIJ5aA4eEsqO/zlBrgKQi36ApkJVkA83zZif4fuPQysUQCAABx4ExwM17IvnEHr/k/Z20FZLQ9BPK1HUFXQIBhwytDJUt0dmTKjdJgTn5ifiuvNXb6uBYtQHHqg0AgCsmZ/T6oCCEBEYlE6EoVdWRtfYEQG0WB5QSkc+JlgyOTCxEYYoKhSkqWB1u6K0OmO1uOFw8tAoxkpSSgPccFAg4SAXCXtdFWtABEAkNlczzox9qDZCkRwA0UMM8tUwE+QDBSM9faqlIiO91dDrmeYZWiwNONx/UL3DPnhSDoZAKez3OQKvUIu3WBfm4YnIGZGIBbE4eHx2tweFKPdI0Ulw2gXZsJ2QwxqSpkZPom6kWCwVIVQ9++Trpm1wihFwSe/WIgzor7dq1C3/4wx9w+vRpAMCECRPwi1/8AosWLQrp4GJZZ41N907QocgADZT2HWqVv6Cj8DlYAwVdgUj008VUGeTKh+EmEgqQldD1wfHzS0ajqd0OhUQY9qJGQmLR+ExNv0vMCQlU0J/Ab731FpYuXQqFQoF7770X9957L+RyOZYsWYJ33nknHGOMSZ0n7u6doAdTA9Qz1TtQBihpEMFLKIRiCqznlBsAKMTCIW0BEgkpamnUB26ERKOx3fblImSogv4U/u1vf4tNmzbhgQce8F537733YvPmzfjNb36DG264IaQDjFWdGSCrw+1tXBjuGiCRkINGFpkTr1goGHLBsr99bAQCDnKxCGY79aYiJFYppEJMyNBQzRwJqaC/O1+4cAFXXnllr+uvuuoqlJWVhWRQ8UAl89SumLqduAczJSIV95gC6ycDlKCQRHR1kdJPFkgYYNCnlon6/PkE2wCMEBL9OM7zpWdchhrzCpIo+CEhF/QZNycnBzt27Oh1/RdffIGcnJyQDCoeqDsyMe3d9mIRDaIqvtcUWD+P4a+GZjj5qwMKNJ3d39iVtGM9ITFFIRVi3qgkzMxLQHaCIipWDJHYE/RX53Xr1uHee+9FcXExLrroIgDAnj17sGXLFvz5z38O+QBjVWfWonNneA6D28y0VxF0PxmeSAdAPVdwcRyQk6AIaBm7v+mvTpQBIiR2pKilmJipGXRnfEICFfSZ484770R6ejr++Mc/4p///CcAYPz48XjvvfewcuXKkA8wVvUsghUJB7cRn0QY2CowqVgQ8cLbnoXQWrkYcolwwNoggcAzfdcXVYTqmgghoSEUcMhKkCNDK4NaFt2tLUjsGNSZY9WqVVi1alWoxxJX1D1O2oOp/5H06AIN9B0ARTr7A/QOgDqX0yskQhj6CYC0cnG/2TF5x0ow6tFJSPTy9JIRotXk8LleJhZiSo4WGgp8yDAb9Fdnh8OBxsbGXp2hc3NzhzyoeCAVCSDggM5u4IOr/+kdNEVzACTv0WCxsyfRQNtkpGn6b27GcRwUEhFMNloJRkg0StVIMT5DAyHH4XiNAU0de/ilqKUYl6Gmzs0kIoIOgM6fP49bb70Ve/fu9bmeMUZ7gQWB4zjIxUKYHZ6f12Dmu6V+2pD31QdIJ498ACQSCqCQCGFxuKGQCL1Tcv01SRQIBg6AAE8dEAVA8Y2ygNFHIACKUtXISezaXmdylhYVrRakqKVUv0ciKujfvjVr1kAkEuGTTz5BRkYGbdo4BHJJVwA0mB5APet/AP8ZIIEAkImjo6BwTkEiKlst6L4NXs/MUHcpKllA04ORrm8ioSUVC6CReWrEAt2VemKmFk3tdtQbbEN+bqebp2BqiNQyESZmaXsFOQIBh4JkZYRGRUiXoM8axcXFOHToEMaNGxeO8cSV7if+QXWB9hPU+AuAZCJh1ASqIqEAo1JUPtf11yU6QxfY3j70TXLkkogEYIC3ED47UY4xqWpvfZtKKsKpWmO/jzEqRYk0jQxauRhN7XZvc9FgpGtlKExRQS4R4nSdETVt1qAfgwA6hRj5ycpBbZlDyHAK+qwxYcIENDc3h2Mscaf71E8oukADfQRAIdqJPVz62lFYKhYgKcDapZ5F5WTkGJ+hQaJSglq9FSIhhwytb2+oTJ0cbp7hbH273/unaWTeoFomFiIvSYELTeagxpCmkWFipsb7RSE3UUEBUJCEAg6jU1U+012ERLOg0w6///3v8dBDD2Hnzp1oaWmB0Wj0+UcCp/AJgIa+Dxjgvw9QKPbhCieZWOg3cEvXyALOXMnEQiioIeKIk6aRIUUthVDAISdR0Sv46ZSTqMCUbG2v3xOtQowJmRqf6/KSlANuvisSep4vUydHdqLcJ/gBPFOqyUPcODgQaRoZFoxOHpbnCheF1BN0zh2VSMEPGVGC/tq8dOlSAMCSJUt8rqci6ODJxV0//kGtAvMzBeZvJk0hjv7siEws7LWfV/dd1AORqJTAYqdv7SOFSMihKE018IEdUjUyyCVCnG80wWB1QiwQ+A2KhAIOM/MScKzaAKO1a3Vh599GikqGMemqAVce5SYq0NyxWikc0rVdWadpOTrU6K04PcBUXyeOQ0B1UeHSGbjqFOI+M7iERLugz4xfffVVOMYRl4aeAep9H3+1RDJJdBRA90ch8Q2AEpTiXp2jB5KokKC6lQKgkWJsujrok6daJsaM3AQwxuDmWZ+rJ2ViIWblJaC8xQylVIREpSTov7FEpQQqWehXFwqFHEYlK5GbqPDJOmXp5LA63Chv9p2+S1FLkZuowIlaA+xOHiIhh+k5CShpakebue/2EeEglwgxLl2NJKrvITEg6ADo4osvDsc44lL3NH1/u7j7w3H+V4H5SyQFG0hEQs8piyxd8Kl02ixx5PB0/Q0uw9cdx3ED/s0IBFyvgvtgpWlkMNlMQ3qMno/XX/apMEUJs92FpnY7JCIBshPkKEhWguM4zM5PxJn6doxOVUElFWFMmhrflbUOWyZIIxdjeq5uUF/WCIlGgzozfv3113j55Zdx4cIFvP/++8jKysL//d//oaCgAAsXLgz1GGPWUDJAEpHAb30Mx3EQCjifVTD9LTOPFt3HKBYJkDqImgiJSAC1TOSzwSyJPmqZCGPT1JEeRkASFRKUBnisUMjB7fYfjfjrh+MPx3GYnKXtuI/v37dMLMS0HJ33slomRqZOPmCxtlomQlGaGlanG2a7Cya7CyabC45+uq8DwPhMDVxuHheazZCLhRT8kJgTdAD0wQcf4KabbsKNN96Iw4cPw273zJEbDAb87ne/w3/+85+QDzJWdc/MBFsD5C/700nQLQCSiASD2mR1uHXPAGVoZYPe/TlRKaEAaBho5GJo5WJUtQ68kW13QiGHKdm6EbO7t0YugkjIwdVHYAN4srG5iQrkJytxqKKt15SZVCzAlGwdtPLAtnoI5mdTmKJCY7u9z730tAoxpuX4D1x4nsHh5tFmcaCyxeLzdzM2XY0snSdDl66VgQNHwQ+JOUEHQE899RReeukl3HzzzXj33Xe91y9YsABPPfVUSAcX64aaAeqLkOO8G0tE+wqwToqOTVFTVJ56h8FKUEpQEcDu8iQ4iSoJ8pOUEAs5SEVC7++f080H1XhwTJp6wBVa0YTjOCQqJWg0+i+GFgo5TMvWIaGjXcOETA0OdJuW0inEmJytDdtWDxKRANNydDhc2ebNPgkEnv3zEhQS5CYq+qyTEgg4yARCZGg905FmuwuWjsasKd0ysLRNBYlVQQdAZ8+exeLFi3tdr9VqodfrQzGmuKHsngEKsgaovwCoex30SFmhIRcLsbgoecgNG3VyMW2JEEICgSdoyU7wH5QWpanQZLL3OfXTXbJa6s0qjCR9BUBCIYcZOQnQKroyOxqZGHlJStQbbChIUSJTG3grh8HSysWYmq3D0Wo9MrQy5CcpB/V3r5SKqKM6iStB5zTT09NRUlLS6/pvvvkGo0aNCsmg4kX3vjXBdoL2twLM32ONlAwQx3EhOVGIhALkJVGb/VAZl67pM/gBPNmB0QEUGieqJBiXPjLqfnpKUvauR1NKRZiR6xv8dBqVrMRFhUnI0smHrQN7olKCxUUpGJeuGTFfegiJtKDD/dtvvx333XcfXnvtNXAch9raWuzbtw8PPvggHnvssXCMMWYphtAJWiLs+0Oue8Z7JKwAC7WCJCWa2+1UCzREWQlyZAaQsclOkKO8xQy7s3faTacQY1yGZkRvVSKXCCGXCGF1uCETCzE6VYU0jbTP4CZS9U0jodaPkGgS9KfSI488Ap7nsWTJElgsFixevBhSqRQPPvgg7rnnnnCMMWZ1nwILtgZILOr7w677jvAjYQVYqAkEHCZmafFdWcuAU2FCgafGQykV9eq/Es+0CnHAK7U4jkOaRobKHrVXcokQU7J1/U7XjhSJSgnccoax6WoqBiYkRgQdAHEch1/96lf4xS9+gZKSEphMJkyYMAEqlQpWqxVy+cib448UpWzwfYD6WwXW/ZvgSGiCGA4qqQhpGhnq9P0X6M7KT4Ba5pnGMNqcaDU5hmN4Ua3nZqSBSFP7BkBCAYcp2dqYCH4AoChV1WcxMSFkZBr0X7REIsGECRMwZ84ciMVibN68GQUFBaEcW8xTSrrqB8RB1gD1WwTdkQESdqzYiVcaWf/LjhNVEm/wAwCjU4fWNC8WTMrSYly6JuhpHK1C7LO6a3yGxudnO9JR8ENI7An4r9put2P9+vWYNWsWLrroImzbtg0A8Prrr6OgoAB/+tOf8MADD4RrnDFJOYRl8P0d35lNSonzdvUD1Z3k9yiW1sjESNPIwjmkqFaYqkK6dvCvv/Nnl6KWDulxCCFkOAQ8Bfb444/j5ZdfxtKlS7F371788Ic/xC233IJvv/0Wmzdvxg9/+EMI+ynMJb1JxUKIBBxcPAvtFFhHBig7yM1EY41K1vevt1rm2R+qp8JUJZpMtrhbRp+mkaEgeWir59K1MtTorRiXMTJXexFC4kvAAdD777+PN998E1dddRVOnDiBKVOmwOVy4ejRo8O21DPWcJynT4/J7goqABKLBP1OUQgEHFQyUdzvjSUWCryrd3rqa6m8QiLCpEwtjlUbwj28qJGokmBCpmbIj6OSijAtWxfX066EkJEj4HmX6upqzJw5EwAwadIkSKVSPPDAAxT8DIGA47yrtIKZAhtoybyQ4+I++9PJ3zSYSMj1u9dYqkaGorT4qAdK18owLVsXsiXU/vriEEJINAr4rOt2uyGRdGUURCIRVKr4OEmECwdPnxRg4HqV7vprggh49h5Kj+Nalu7UfqbBUtUD7zWWl6TEhEzNiNq2YSDCHoFzbpICk7K0I2ZfLkIICaWAz7qMMaxZswZSqeebs81mw89+9jMolb5TCVu3bg3tCGOYgONw07w8lLeYMSqI+ov+miACnhM8NUXz8FcHFGiBbqZOjnSNDBWtFpQ2mkI9tGGjkAgxNl0NjVyMs/Xt3m0aCgPo4EwIIbEq4ABo9erVPpd//OMfh3ww8YbjPCfZQLrtdjdQbxUKfrqopb5TMhKRAAlBTNMIBBwKkpVwuPigdz6PNI7zZLJGJSu9WZ5JWVrkJikGbBFACCGxLuAA6PXXXw/nOOLSYMungt02I57JJUKIhBxcHZt1pvazhUF/xqSpYHG40DJCGiVKxQJMydZBK+8d6FDwQwghQ2iESIZusAXksdJdd7h0rwMabG0Ux3EYn6EZdNA63IpS1X6DH0IIIR50Jo2gwc5UUQAUHJVUDI7zNOgbSlAgEwuRNAKaSyqlIqRpon+chBASSSN3i+YYwGGQGSBqyx+UnEQ5CpKVIQkcM3UyNLfbQzCq0BEKOSQpJWg0esZVmKqk9hSEEDIACoAiiDJAw0MhCd2vebJSColIAIcrOlpFS0QCTM/VQS0To9FoQ43eilQ1tUAghJCBUAAUSYMNgCgDFDECAYdMnQzlzZFfESYTCzEzL8HbqyhVI0Mq9X8ihJCA0Jk0ggSDmKYQCGhn6kgLtm1BOIg7Mj+x1KiREEKGE2WAImgwCaCBmiCS8FNIRJCJhbA5e+8xFi4CATArPxF6sxM1eismZGqgDKJ7OCGEEF/0CRpBg8kAUf1PdFDLRMMaAOUlKaGRiaGRiZGbpBi25yWEkFhFZ9MIGsxCHQqAooNmGHvsyCVC5Pexez0hhJDBobNpBA1mqTJ1gY4OGj97jIXLmDQ1bW9CCCEhFtEAaPfu3bjyyiuRmZkJjuOwbds2n9sZY3j88ceRkZEBuVyOpUuX4vz580N6zGgymHOaVEQ1QNFAPUzbSSgkQqSoqakhIYSEWkQDILPZjKlTp+KFF17we/umTZvw3HPP4aWXXsL+/fuhVCqxbNky2Gy2QT9mNBlMBohW/UQHiUgAmbjrvVDLRBiXocacUYkhfZ6MKFhxRgghsSiiRdArVqzAihUr/N7GGMOzzz6LRx99FCtXrgQAvPnmm0hLS8O2bdtw3XXXBf2YfbHb7bDbu7r7Go3GoO4/WIPJAMnFFABFC428qxB6QqbGmxXSyMUwWp2DekyhgIOb92zcynFAhpb6+hBCSDhEbQ1QWVkZ6uvrsXTpUu91Wq0Wc+fOxb59+0L6XBs3boRWq/X+y8nJCenj92UwW2HIxFH7lsWdzoBHpxD7TIklqyRBP9bYdDUWFiVj8ZgUb4F1olLik2UihBASOlF7Nq2vrwcApKWl+VyflpbmvS1U1q9fD4PB4P1XVVUV0sfvS7AzYBwHyKgGKGp0FkJnJ/guSw92w1StQoycRAVkYiGEAg5Tc7SQiYXIoukvQggJG+oDBEAqlUIqHf5C02ADIKlICAGtBooaapkYYpEAqT2KlDUyUVD7hRWmqHwuS0VCTMvVQUHZH0IICZuozQClp6cDABoaGnyub2ho8N420gXbCJGmv6KLRCTA6FRVr6CU4zgkBTgNlqAUI1HZ+1iVVETBLiGEhFHUnlELCgqQnp6OHTt2eK8zGo3Yv38/5s+fH8GRhU6wpzeqB4k+mX0UKScHOA3WM/tDCCFkeER0CsxkMqGkpMR7uaysDMXFxUhMTERubi7uv/9+PPXUUygqKkJBQQEee+wxZGZm4uqrr/beZ8mSJVi1ahXuvvvugB4zmgSbAaIl8NGnr1YGiUoJOA5grO/7Jqkk0CmCL5gmhBAydBENgA4ePIhLL73Ue3nt2rUAgNWrV2PLli146KGHYDabcccdd0Cv12PhwoXYvn07ZLKub92lpaVobm4O+DGjSbA1QLQEfuQQCwXISpCjutXa5zGjUyn7QwghkcIx1t931PhkNBqh1WphMBig0WjC+lxfnGoY+KAOM/MSkOCnXoREJ7vLjb2lLXC7e/+JpWtlmJSljcCoCCEkdgVz/o7aGqB4IQjiHaAaoJFFKvLdxFQhFUIh8Sx1H5VCm5sSQkgk0TL4CPPUkAychOM4WgU2EuUmKtBmcSA7QY5UNXV1JoSQaEEBUIQFWgYkEwsHtXcYiSyhgMOM3IRID4MQQkgPlFKIsEBXgtH0FyGEEBI6FABFWKBJHZr+IoQQQkKHzqoRFmgGiJbAE0IIIaFDAVCEBVrVQ00QCSGEkNChACjCAilsVslEAW+tQAghhJCBUQAUYQPFPyqZCDNyEyAW0ltFCCGEhAqdVSNsoBqg6bk6SET0NhFCCCGhRGfWCOsv/hEJOUhFVPtDCCGEhBoFQBEm6CcAouCHEEIICQ8KgCKu7wiIpr4IIYSQ8KAzbIT1nwGit4cQQggJBzrDRlh/y+ApA0QIIYSEB51hI6y/DJCElr4TQgghYUFn2Ajrbxk8ZYAIIYSQ8KAzbBSjAIgQQggJDzrDRhhlgAghhJDhR2fYCOuvESLVABFCCCHhQWfYCOs3A0QBECGEEBIWdIaNsL7iH7FIAEF/S8QIIYQQMmgUAEVYXzEOZX8IIYSQ8KGzbMT5j4CoAJoQQggJHzrLRlhfGSDaBoMQQggJHzrLRlhfW2FQBogQQggJHzrLRhjVABFCCCHDj86yEcZRDRAhhBAy7OgsG2F9LYOnAIgQQggJHzrLRhgFQIQQQsjwo7NshPXVCZpqgAghhJDwobNshPUVANEyeEIIISR86CwbYf7iH7FI0OfyeEIIIYQMHQVAEeYvzqHpL0IIISS86EwbYf6WwVMBNCGEEBJedKaNMH+NECkDRAghhIQXnWkjzF+tj0hI9T+EEEJIOFEAFGH+MkBiygARQgghYUVn2gjzWwNEARAhhBASVnSmjTDOzztAU2CEEEJIeFEAFGH+Qh2aAiOEEELCi860EeavEzRNgRFCCCHhRWfaCPPfCZqmwAghhJBwogAowvxlgEQCelsIIYSQcKIzbRQSUxE0IYQQElYUAEVYzwyQSMjRRqiEEEJImFEAFGE9Yx0qgCaEEELCj862EdYzAySmjVAJIYSQsKOzbYT1nOwS+dsbgxBCCCEhRQFQhAl6BDzUBJEQQggJPzrbRoHus2ASmgIjhBBCwo7OtlGgex0QZYAIIYSQ8KOzbTTolgGiGiBCCCEk/CgAigLdM0A0BUYIIYSEH51to0D3nA9lgAghhJDwowAoCvjUAFEGiBBCCAk7OttGAZ9VYFQETQghhIQdnW2jQPcAiFaBEUIIIeFHZ9sowHVUAQkEgJBqgAghhJCwowAoCnTGPJT9IYQQQoYHnXGjQOd2GBQAEUIIIcODzrhRIFklBQCIhTT9RQghhAyHiAZAu3fvxpVXXonMzExwHIdt27b53M4Yw+OPP46MjAzI5XIsXboU58+fH/BxX3jhBeTn50Mmk2Hu3Ln47rvvwvQKQiM3UQG5REgZIEIIIWSYRPSMazabMXXqVLzwwgt+b9+0aROee+45vPTSS9i/fz+USiWWLVsGm83W52O+9957WLt2LTZs2IDDhw9j6tSpWLZsGRobG8P1MoZMKOBQlKqCSEABECGEEDIcOMYYi/QgAIDjOHz44Ye4+uqrAXiyP5mZmVi3bh0efPBBAIDBYEBaWhq2bNmC6667zu/jzJ07F7Nnz8bzzz8PAOB5Hjk5ObjnnnvwyCOP+L2P3W6H3W73XjYajcjJyYHBYIBGownhq+xfi8mOpI7pMEIIIYQEx2g0QqvVBnT+jtqUQ1lZGerr67F06VLvdVqtFnPnzsW+ffv83sfhcODQoUM+9xEIBFi6dGmf9wGAjRs3QqvVev/l5OSE7oUEgYIfQgghZHhEbQBUX18PAEhLS/O5Pi0tzXtbT83NzXC73UHdBwDWr18Pg8Hg/VdVVTXE0RNCCCEkmokiPYBoIJVKIZVS9oUQQgiJF1GbAUpPTwcANDQ0+Fzf0NDgva2n5ORkCIXCoO5DCCGEkPgTtQFQQUEB0tPTsWPHDu91RqMR+/fvx/z58/3eRyKRYObMmT734XkeO3bs6PM+hBBCCIk/EZ0CM5lMKCkp8V4uKytDcXExEhMTkZubi/vvvx9PPfUUioqKUFBQgMceewyZmZnelWIAsGTJEqxatQp33303AGDt2rVYvXo1Zs2ahTlz5uDZZ5+F2WzGLbfcMtwvjxBCCCFRKqIB0MGDB3HppZd6L69duxYAsHr1amzZsgUPPfQQzGYz7rjjDuj1eixcuBDbt2+HTCbz3qe0tBTNzc3ey9deey2amprw+OOPo76+HtOmTcP27dt7FUYTQgghJH5FTR+gaBJMHwFCCCGERIeY6ANECCGEEBIuFAARQgghJO5QAEQIIYSQuEMBECGEEELiDgVAhBBCCIk7FAARQgghJO5QAEQIIYSQuEMBECGEEELiDu0G70dnb0ij0RjhkRBCCCEkUJ3n7UB6PFMA5Ed7ezsAICcnJ8IjIYQQQkiw2tvbodVq+z2GtsLwg+d51NbWQq1Wg+O4IT+e0WhETk4OqqqqYnZrDXqNI1+svz6AXmMsiPXXB9BrHArGGNrb25GZmQmBoP8qH8oA+SEQCJCdnR3yx9VoNDH7y9yJXuPIF+uvD6DXGAti/fUB9BoHa6DMTycqgiaEEEJI3KEAiBBCCCFxhwKgYSCVSrFhwwZIpdJIDyVs6DWOfLH++gB6jbEg1l8fQK9xuFARNCGEEELiDmWACCGEEBJ3KAAihBBCSNyhAIgQQgghcYcCIEIIIYTEHQqAhsELL7yA/Px8yGQyzJ07F999912khzQoGzduxOzZs6FWq5Gamoqrr74aZ8+e9TnmkksuAcdxPv9+9rOfRWjEwXviiSd6jX/cuHHe2202G+666y4kJSVBpVLhBz/4ARoaGiI44uDl5+f3eo0cx+Guu+4CMPLew927d+PKK69EZmYmOI7Dtm3bfG5njOHxxx9HRkYG5HI5li5divPnz/sc09raihtvvBEajQY6nQ633XYbTCbTML6K/vX3Gp1OJx5++GFMnjwZSqUSmZmZuPnmm1FbW+vzGP7e96effnqYX0nfBnof16xZ02v8y5cv9zkmmt/HgV6fv79JjuPwzDPPeI+J5vcwkPNDIJ+flZWVuOKKK6BQKJCamopf/OIXcLlcYRkzBUBh9t5772Ht2rXYsGEDDh8+jKlTp2LZsmVobGyM9NCCtmvXLtx111349ttv8fnnn8PpdOLyyy+H2Wz2Oe72229HXV2d99+mTZsiNOLBmThxos/4v/nmG+9tDzzwAD7++GO8//772LVrF2pra3HNNddEcLTBO3DggM/r+/zzzwEAP/zhD73HjKT30Gw2Y+rUqXjhhRf83r5p0yY899xzeOmll7B//34olUosW7YMNpvNe8yNN96IkydP4vPPP8cnn3yC3bt344477hiulzCg/l6jxWLB4cOH8dhjj+Hw4cPYunUrzp49i6uuuqrXsU8++aTP+3rPPfcMx/ADMtD7CADLly/3Gf8//vEPn9uj+X0c6PV1f111dXV47bXXwHEcfvCDH/gcF63vYSDnh4E+P91uN6644go4HA7s3bsXb7zxBrZs2YLHH388PINmJKzmzJnD7rrrLu9lt9vNMjMz2caNGyM4qtBobGxkANiuXbu811188cXsvvvui9yghmjDhg1s6tSpfm/T6/VMLBaz999/33vd6dOnGQC2b9++YRph6N13332ssLCQ8TzPGBvZ7yEA9uGHH3ov8zzP0tPT2TPPPOO9Tq/XM6lUyv7xj38wxhg7deoUA8AOHDjgPea///0v4ziO1dTUDNvYA9XzNfrz3XffMQCsoqLCe11eXh7705/+FN7BhYi/17h69Wq2cuXKPu8zkt7HQN7DlStXsu9973s+142k97Dn+SGQz8///Oc/TCAQsPr6eu8xL774ItNoNMxut4d8jJQBCiOHw4FDhw5h6dKl3usEAgGWLl2Kffv2RXBkoWEwGAAAiYmJPte//fbbSE5OxqRJk7B+/XpYLJZIDG/Qzp8/j8zMTIwaNQo33ngjKisrAQCHDh2C0+n0eT/HjRuH3NzcEft+OhwOvPXWW7j11lt9Nv4d6e9hp7KyMtTX1/u8Z1qtFnPnzvW+Z/v27YNOp8OsWbO8xyxduhQCgQD79+8f9jGHgsFgAMdx0Ol0Ptc//fTTSEpKwvTp0/HMM8+EbWohXHbu3InU1FSMHTsWd955J1paWry3xdL72NDQgH//+9+47bbbet02Ut7DnueHQD4/9+3bh8mTJyMtLc17zLJly2A0GnHy5MmQj5E2Qw2j5uZmuN1unzcTANLS0nDmzJkIjSo0eJ7H/fffjwULFmDSpEne62+44Qbk5eUhMzMTx44dw8MPP4yzZ89i69atERxt4ObOnYstW7Zg7NixqKurw69//WssWrQIJ06cQH19PSQSSa+TSlpaGurr6yMz4CHatm0b9Ho91qxZ471upL+H3XW+L/7+Bjtvq6+vR2pqqs/tIpEIiYmJI/J9tdlsePjhh3H99df7bDJ57733YsaMGUhMTMTevXuxfv161NXVYfPmzREcbeCWL1+Oa665BgUFBSgtLcUvf/lLrFixAvv27YNQKIyp9/GNN96AWq3uNb0+Ut5Df+eHQD4/6+vr/f6tdt4WahQAkUG56667cOLECZ/6GAA+8+2TJ09GRkYGlixZgtLSUhQWFg73MIO2YsUK7/9PmTIFc+fORV5eHv75z39CLpdHcGTh8fe//x0rVqxAZmam97qR/h7GM6fTiR/96EdgjOHFF1/0uW3t2rXe/58yZQokEgl++tOfYuPGjSNiy4XrrrvO+/+TJ0/GlClTUFhYiJ07d2LJkiURHFnovfbaa7jxxhshk8l8rh8p72Ff54doQ1NgYZScnAyhUNiryr2hoQHp6ekRGtXQ3X333fjkk0/w1VdfITs7u99j586dCwAoKSkZjqGFnE6nw5gxY1BSUoL09HQ4HA7o9XqfY0bq+1lRUYEvvvgCP/nJT/o9biS/h53vS39/g+np6b0WJbhcLrS2to6o97Uz+KmoqMDnn3/uk/3xZ+7cuXC5XCgvLx+eAYbYqFGjkJyc7P29jJX38euvv8bZs2cH/LsEovM97Ov8EMjnZ3p6ut+/1c7bQo0CoDCSSCSYOXMmduzY4b2O53ns2LED8+fPj+DIBocxhrvvvhsffvghvvzySxQUFAx4n+LiYgBARkZGmEcXHiaTCaWlpcjIyMDMmTMhFot93s+zZ8+isrJyRL6fr7/+OlJTU3HFFVf0e9xIfg8LCgqQnp7u854ZjUbs37/f+57Nnz8fer0ehw4d8h7z5Zdfgud5b/AX7TqDn/Pnz+OLL75AUlLSgPcpLi6GQCDoNW00UlRXV6OlpcX7exkL7yPgycrOnDkTU6dOHfDYaHoPBzo/BPL5OX/+fBw/ftwnkO0M5idMmBCWQZMwevfdd5lUKmVbtmxhp06dYnfccQfT6XQ+Ve4jxZ133sm0Wi3buXMnq6ur8/6zWCyMMcZKSkrYk08+yQ4ePMjKysrYRx99xEaNGsUWL14c4ZEHbt26dWznzp2srKyM7dmzhy1dupQlJyezxsZGxhhjP/vZz1hubi778ssv2cGDB9n8+fPZ/PnzIzzq4Lndbpabm8sefvhhn+tH4nvY3t7Ojhw5wo4cOcIAsM2bN7MjR454V0A9/fTTTKfTsY8++ogdO3aMrVy5khUUFDCr1ep9jOXLl7Pp06ez/fv3s2+++YYVFRWx66+/PlIvqZf+XqPD4WBXXXUVy87OZsXFxT5/m50rZ/bu3cv+9Kc/seLiYlZaWsreeustlpKSwm6++eYIv7Iu/b3G9vZ29uCDD7J9+/axsrIy9sUXX7AZM2awoqIiZrPZvI8Rze/jQL+njDFmMBiYQqFgL774Yq/7R/t7OND5gbGBPz9dLhebNGkSu/zyy1lxcTHbvn07S0lJYevXrw/LmCkAGgZ/+ctfWG5uLpNIJGzOnDns22+/jfSQBgWA33+vv/46Y4yxyspKtnjxYpaYmMikUikbPXo0+8UvfsEMBkNkBx6Ea6+9lmVkZDCJRMKysrLYtddey0pKSry3W61W9vOf/5wlJCQwhULBVq1axerq6iI44sH59NNPGQB29uxZn+tH4nv41Vdf+f29XL16NWPMsxT+scceY2lpaUwqlbIlS5b0et0tLS3s+uuvZyqVimk0GnbLLbew9vb2CLwa//p7jWVlZX3+bX711VeMMcYOHTrE5s6dy7RaLZPJZGz8+PHsd7/7nU/wEGn9vUaLxcIuv/xylpKSwsRiMcvLy2O33357ry+S0fw+DvR7yhhjL7/8MpPL5Uyv1/e6f7S/hwOdHxgL7POzvLycrVixgsnlcpacnMzWrVvHnE5nWMbMdQycEEIIISRuUA0QIYQQQuIOBUCEEEIIiTsUABFCCCEk7lAARAghhJC4QwEQIYQQQuIOBUCEEEIIiTsUABFCCCEk7lAARAghhJC4QwEQIXHkzJkzmDdvHmQyGaZNm+b3mEsuuQT333//kJ6nvLwcHMd59xGLdmvWrMHVV18d6WFElSeeeKLP3xFCYoEo0gMghPTW1NSErKwstLW1QSKRQKfT4fTp08jNzR3S427YsAFKpRJnz56FSqXye8zWrVshFouH9DzEE0Bs27ZtxASBhMQbCoAIiUL79u3D1KlToVQqsX//fiQmJg45+AGA0tJSXHHFFcjLy+vzmMTExCE/DyGERDuaAiMkCu3duxcLFiwAAHzzzTfe/+8Pz/N48sknkZ2dDalUimnTpmH79u3e2zmOw6FDh/Dkk0+C4zg88cQTfh+n5xRYfn4+fve73+HWW2+FWq1Gbm4uXnnlFZ/7fPfdd5g+fTpkMhlmzZqFI0eO9HrcEydOYMWKFVCpVEhLS8NNN92E5uZmn+e9++67cffdd0Or1SI5ORmPPfYYum9XaLfb8eCDDyIrKwtKpRJz587Fzp07vbdv2bIFOp0On376KcaPHw+VSoXly5ejrq7Oe4zb7cbatWuh0+mQlJSEhx56CD23ROR5Hhs3bkRBQQHkcjmmTp2K//f//p/39p07d4LjOOzYsQOzZs2CQqHARRddhLNnz3rH8etf/xpHjx4Fx3HgOA5btmzx+/PeuXMn5syZA6VSCZ1OhwULFqCiogKAJ2BduXIl0tLSoFKpMHv2bHzxxRc+98/Pz8dTTz2Fm2++GSqVCnl5efjXv/6FpqYmrFy5EiqVClOmTMHBgwd7/Zy2bduGoqIiyGQyLFu2DFVVVX7H2Olvf/sbxo8fD5lMhnHjxuGvf/2r9zaHw4G7774bGRkZkMlkyMvLw8aNG/t9PEIiKixbrBJCglZRUcG0Wi3TarVMLBYzmUzGtFotk0gkTCqVMq1Wy+68884+779582am0WjYP/7xD3bmzBn20EMPMbFYzM6dO8cYY6yuro5NnDiRrVu3jtXV1fW5S/bFF1/M7rvvPu/lvLw8lpiYyF544QV2/vx5tnHjRiYQCNiZM2cYY4y1t7ezlJQUdsMNN7ATJ06wjz/+mI0aNYoBYEeOHGGMMdbW1sZSUlLY+vXr2enTp9nhw4fZZZddxi699FKf51WpVOy+++5jZ86cYW+99RZTKBTslVde8R7zk5/8hF100UVs9+7drKSkhD3zzDNMKpV6X+Prr7/OxGIxW7p0KTtw4AA7dOgQGz9+PLvhhhu8j/H73/+eJSQksA8++ICdOnWK3XbbbUytVrOVK1d6j3nqqafYuHHj2Pbt21lpaSl7/fXXmVQqZTt37mSMde3sPXfuXLZz50528uRJtmjRInbRRRcxxhizWCxs3bp1bOLEiayuro7V1dUxi8XS62ftdDqZVqtlDz74ICspKWGnTp1iW7ZsYRUVFYwxxoqLi9lLL73Ejh8/zs6dO8ceffRRJpPJvLd3f39eeukldu7cOXbnnXcyjUbDli9fzv75z3+ys2fPsquvvpqNHz+e8Tzv83OaNWsW27t3Lzt48CCbM2eOd/yMMbZhwwY2depU7+W33nqLZWRksA8++IBduHCBffDBBywxMZFt2bKFMcbYM888w3Jyctju3btZeXk5+/rrr9k777zj93eMkGhAARAhUcLpdLKysjJ29OhRJhaL2dGjR1lJSQlTqVRs165drKysjDU1NfV5/8zMTPbb3/7W57rZs2ezn//8597LU6dOZRs2bOh3HP4CoB//+MfeyzzPs9TUVPbiiy8yxhh7+eWXWVJSErNard5jXnzxRZ8A6De/+Q27/PLLfZ6nqqqKAWBnz571Pm/3kzRjjD388MNs/PjxjDFPgCgUCllNTY3P4yxZsoStX7+eMeY5sQNgJSUl3ttfeOEFlpaW5r2ckZHBNm3a5L3sdDpZdna2NwCy2WxMoVCwvXv3+jzPbbfdxq6//nrGWFcA9MUXX3hv//e//80AeH8OPQMIf1paWhgAb2AViIkTJ7K//OUv3ss935+6ujoGgD322GPe6/bt28cAsLq6OsZY18/p22+/9R5z+vRpBoDt37/f7/gLCwt7BTS/+c1v2Pz58xljjN1zzz3se9/7ns/7R0g0oykwQqKESCRCfn4+zpw5g9mzZ2PKlCmor69HWloaFi9ejPz8fCQnJ/u9r9FoRG1tba+psgULFuD06dNDHtuUKVO8/89xHNLT09HY2AgAOH36NKZMmQKZTOY9Zv78+T73P3r0KL766iuoVCrvv3HjxgHwTPN0mjdvHjiO83mc8+fPw+124/jx43C73RgzZozP4+zatcvnMRQKBQoLC72XMzIyvGM1GAyoq6vD3LlzvbeLRCLMmjXLe7mkpAQWiwWXXXaZz/O8+eabPs/T8+eSkZEBAN7nCkRiYiLWrFmDZcuW4corr8Sf//xnn+k6k8mEBx98EOPHj4dOp4NKpcLp06dRWVnZ5zjS0tIAAJMnT+51XfexiUQizJ4923t53Lhx3mL7nsxmM0pLS3Hbbbf5/Eyeeuop789kzZo1KC4uxtixY3Hvvffis88+C/jnQEgkUBE0IVFi4sSJqKiogNPpBM/zUKlUcLlccLlc3tqOkydPRmRsPVeFcRwHnucDvr/JZMKVV16J3//+971u6wwcAnkMoVCIQ4cOQSgU+tzWfUWbv7GyHjU+Az0PAPz73/9GVlaWz21SqdTncvfn6gzcgvm5AMDrr7+Oe++9F9u3b8d7772HRx99FJ9//jnmzZuHBx98EJ9//jn+8Ic/YPTo0ZDL5fjf//1fOByOAccRirF16vyZvPrqqz7BIwDvezFjxgyUlZXhv//9L7744gv86Ec/wtKlS31qpwiJJhQAERIl/vOf/8DpdGLJkiXYtGkTZs6cieuuuw5r1qzB8uXL+12artFokJmZiT179uDiiy/2Xr9nzx7MmTMnrOMeP348/u///g82m82bBfr22299jpkxYwY++OAD5OfnQyTq+2Nn//79Ppe//fZbFBUVQSgUYvr06XC73WhsbMSiRYsGNVatVouMjAzs378fixcvBgC4XC4cOnQIM2bMAABMmDABUqkUlZWVPj/LYEkkErjd7oCOnT59OqZPn47169dj/vz5eOeddzBv3jzs2bMHa9aswapVqwB4ApHy8vJBj6k7l8uFgwcPen8/zp49C71ej/Hjx/c6Ni0tDZmZmbhw4QJuvPHGPh9To9Hg2muvxbXXXov//d//xfLly9Ha2korC0lUogCIkCiRl5eH+vp6NDQ0YOXKleA4DidPnsQPfvCDgLIkv/jFL7BhwwYUFhZi2rRpeP3111FcXIy33347rOO+4YYb8Ktf/Qq333471q9fj/LycvzhD3/wOeauu+7Cq6++iuuvvx4PPfQQEhMTUVJSgnfffRd/+9vfvFmEyspKrF27Fj/96U9x+PBh/OUvf8Ef//hHAMCYMWNw44034uabb8Yf//hHTJ8+HU1NTdixYwemTJmCK664IqDx3nfffXj66adRVFSEcePGYfPmzdDr9d7b1Wo1HnzwQTzwwAPgeR4LFy6EwWDAnj17oNFosHr16oCeJz8/H2VlZSguLkZ2djbUanWvDFJZWRleeeUVXHXVVcjMzMTZs2dx/vx53HzzzQCAoqIibN26FVdeeSU4jsNjjz026CxOT2KxGPfccw+ee+45iEQi3H333Zg3b16fAfOvf/1r3HvvvdBqtVi+fDnsdjsOHjyItrY2rF27Fps3b0ZGRgamT58OgUCA999/H+np6dDpdCEZLyGhRgEQIVFk586dmD17NmQyGb7++mtkZ2cHPEV07733wmAwYN26dWhsbMSECRPwr3/9C0VFRWEds0qlwscff4yf/exnmD59OiZMmIDf//73+MEPfuA9pjM79fDDD+Pyyy+H3W5HXl4eli9fDoGgqxTx5ptvhtVqxZw5cyAUCnHffffhjjvu8N7++uuv46mnnsK6detQU1OD5ORkzJs3D9///vcDHu+6detQV1eH1atXQyAQ4NZbb8WqVatgMBi8x/zmN79BSkoKNm7ciAsXLkCn02HGjBn45S9/GfDz/OAHP8DWrVtx6aWXQq/X4/XXX8eaNWt8jlEoFDhz5gzeeOMNtLS0ICMjA3fddRd++tOfAgA2b96MW2+9FRdddBGSk5Px8MMPw2g0BjyG/igUCjz88MO44YYbUFNTg0WLFuHvf/97n8f/5Cc/gUKhwDPPPINf/OIXUCqVmDx5srdlglqtxqZNm3D+/HkIhULMnj0b//nPf3zeX0KiCceCmRwnhJAwueSSSzBt2jQ8++yzkR5KzNuyZQvuv/9+n8wXIfGGQnNCCCGExB0KgAghhBASd2gKjBBCCCFxhzJAhBBCCIk7FAARQgghJO5QAEQIIYSQuEMBECGEEELizv9vtw4EAAAAAAT5Ww9yUSRAAMCOAAEAOwIEAOwIEACwEzZu7s/4Ij/zAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from pyepfd.epfd import mc_convergence\n", "import matplotlib.pyplot as plt\n", "co2_100K = mc_convergence(file_path='gap.dat',algo='mcap', usecols=(0))\n", "mc = co2_100K.avg[:,0] #Our gap was in column 0 of the file\n", "emc = 2*co2_100K.sd_mean[:,0] # Error is 2 times standard deviation of mean\n", "plt.plot( [i+1 for i in range(len(mc))], mc, label='MC')\n", "plt.fill_between([i+1 for i in range(len(mc))], mc + emc, mc - emc, alpha = 0.3)\n", "plt.xlabel(\"# of independent samples\")\n", "plt.ylabel(\"Renormalized Gap (eV)\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "40b4b070", "metadata": {}, "source": [ "Note, the result is different compared to what we obtained in exercise 1.4.1 because the band gap is dependent on the used basis set. Therefore, it should be repeated using larger and more accurate basis sets. " ] }, { "cell_type": "markdown", "id": "a429b431", "metadata": {}, "source": [ "Trajectory analysis, for example anharmonic measure can be computed using the following script" ] }, { "cell_type": "code", "execution_count": 1, "id": "8834e299", "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-20 13:42:21.620006\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", "Process-id0: Time spent on read_pyepfd_info class: 0.0003857612609863281 s.\n", "anharm_measure: Execution time (s): 0.08104157447814941\n" ] } ], "source": [ "from pyepfd.anharm import anharm_measure\n", "#We need pyepfd_io module to parse and store restart file\n", "from pyepfd.pyepfd_io import *\n", "#We need numpy to parse and process the energy.npz\n", "import numpy as np\n", "\n", "sd_inp = read_pyepfd_info(file_path='co2.xml')\n", "if sd_inp.ref_dynmatrix is not None:\n", " inp_dynmat = sd_inp.ref_dynmatrix\n", "else:\n", " inp_dynmat = sd_inp.dynmatrix\n", "\n", "trajectory = np.load('MCAP_100K/engrad.npz')\n", "coords = trajectory['coords']\n", "forces = trajectory['forces']\n", "\n", "\n", "# Now we calculate anharm measure\n", "anharm = anharm_measure(dynmat = inp_dynmat, \\\n", " mass = sd_inp.mass, \\\n", " forces = forces, \\\n", " disp_coords = coords, \\\n", " opt_coord = sd_inp.coord, \\\n", " remove_rot_trans = True,\\\n", " asr = sd_inp.asr)\n", "\n", "# We write the informations using write method in anharm object.\n", "anharm.write(file_path='co2-100K-mcap',atoms=sd_inp.atoms)\n", "\n", "# Task complete, we are deleting the anharm object to finish file writing\n", "del anharm\n" ] }, { "cell_type": "markdown", "id": "acb3591f", "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": 2, "id": "90809699", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "# Mode Freq(cm-1) Var(Anh) Var(Tot) Anh_Measure \n", "#--------------------------------------------------------------\n", " 1 -0.0120 0.000307268 5.89657e-05 2.28275 \n", " 2 0.0001 0.000142152 2.67528e-06 7.28941 \n", " 3 0.0006 0.0168978 3.86368e-07 209.129 \n", " 4 0.0070 0.000306387 1.16596e-06 16.2104 \n", " 5 0.0070 0.000141753 1.8682e-05 2.75457 \n", " 6 676.4269 0.0176427 0.000183404 9.80794 \n", " 7 676.4269 0.00120172 0.000475885 1.5891 \n", " 8 1371.0021 0.000555802 0.00415696 0.365655 \n", " 9 2408.0634 0.0373526 0.0142452 1.6193 \n", "#--------------------------------------------------------------\n", "# SUM 0.00828314 0.00212703 1.97338\n" ] } ], "source": [ "%%bash\n", "cat co2-100K-mcap_mode_res_anh_mes.out" ] }, { "cell_type": "markdown", "id": "262795c2", "metadata": {}, "source": [ "For the modes 6-9, we see similar trends as we have seen in exercise 1.4.1, i.e., higher anharmonic measure for bending modes and symmetric stretching has the lowest anharmonicity. However the obtained anharmonic measure is much larger. This is because, the forces are much more sensitive to the basis set and therefore to obtain force-based anharmonic measure, we need to repeat the whole calculation (starting from begining) with a much better basis set. \n", "\n", "We can also obtaine an energy based anharmonic measure (which is less sensitive); however here we will only get the total anharmonic measure and not the mode(atom) resolved anharmonic measure. For that purpose we can use the following script. " ] }, { "cell_type": "code", "execution_count": 1, "id": "03dd0ed1", "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-20 16:06:20.634063\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", "Process-id0: Time spent on read_pyepfd_info class: 0.0005810260772705078 s.\n", "Time spent on xyz class: 0.005024433135986328 s.\n", "Total anharmonic measure (Energy based) = 0.18793227018070038\n" ] } ], "source": [ "from pyepfd.anharm import boltzmann_reweighting\n", "#We need pyepfd_io module to parse and store restart file\n", "from pyepfd.pyepfd_io import *\n", "from pyepfd.coord_util import xyz\n", "#We need numpy to parse and process the energy.npz\n", "import numpy as np\n", "\n", "sd_inp = read_pyepfd_info(file_path='co2.xml')\n", "if sd_inp.ref_dynmatrix is not None:\n", " inp_dynmat = sd_inp.ref_dynmatrix\n", "else:\n", " inp_dynmat = sd_inp.dynmatrix\n", "\n", "# This can be calculated even when engrad.npz file is not present\n", "disp_energy = np.genfromtxt('MCAP_100K/etotals.dat')\n", "# We are loading the original trajectory\n", "trajectory = xyz(file_path='100K.xyz',io='r')\n", "n_frames = len(disp_energy)\n", "# below we assume only first n_frames are extracted in the previous step\n", "coords = trajectory.coords[:n_frames]\n", "\n", "#We need the energy of the optimized geometry in Hatree\n", "# which we obtain from the orca-opt-freq.engrad file\n", "opt_energy = -188.587352098676\n", "\n", "# Now we calculate anharm measure\n", "anharm = boltzmann_reweighting(dynmat = inp_dynmat, \\\n", " mass = sd_inp.mass, \\\n", " disp_coords = coords, \\\n", " opt_coord = sd_inp.coord, \\\n", " opt_energy = opt_energy,\n", " disp_energy = disp_energy, \n", " remove_rot_trans = True,\\\n", " asr = sd_inp.asr,\n", " temperature = 100)\n", "\n", "# We write the informations using write method in anharm object.\n", "print(f\"Total anharmonic measure (Energy based) = {anharm.anharm_measure}\")\n", "\n", "# Task complete, we are deleting the anharm object to finish file writing\n", "del anharm" ] }, { "cell_type": "markdown", "id": "88cfa666", "metadata": {}, "source": [ "We can see that the energy-based total anharmonic measure is much more consistent with the one we obtained in exercise 1.4.1. Note, for larges sytems and basis sets, force-calculations have significant computational overload. Therefore, while using localized basis set, it is recommended to use energy-based anharmonic measure. " ] } ], "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 }