diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..773a6df --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.dat diff --git a/README.md b/README.md index c56de3b..4f8581b 100644 --- a/README.md +++ b/README.md @@ -3,4 +3,6 @@ This is used to test the performance of TAB on Tully's model. To run TAB dynamics, simply run -'python3 name-main.py' +'python3 namd-main.py' + +Initial conditions, calcH, and diffH must be set manually. diff --git a/analysis.ipynb b/analysis.ipynb new file mode 100644 index 0000000..b06d475 --- /dev/null +++ b/analysis.ipynb @@ -0,0 +1,233 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 82, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import glob\n", + "import os\n", + "# Specify the paths to the directories containing the ene*.dat files\n", + "data_directory1 = '/home/adurden/Ari/TAB/momentum/m9_10000_split/new'\n", + "data_directory2 = '/home/adurden/Ari/TAB/momentum/m9_10000_split/old'" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuwAAAGPCAYAAAAHqIWVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACLG0lEQVR4nOzdd3gU1dfA8e9JQhol9NB7kd6LAhK6YgE7KAI27NgVK+DPgr0DdtFX7KKIKEgJWFFApCq9956EloT7/nEnsISU3ZSd3c35PM88uzs7M/dsZpOcvXvnXDHGoJRSSimllApMYW4HoJRSSimllMqeJuxKKaWUUkoFME3YlVJKKaWUCmCasCullFJKKRXANGFXSimllFIqgGnCrpRSSimlVADThF0pFdBEpJaIGBEZ5XYsBU1EuonIHyKS5LzGoW7HpFQgEpFRzu9ILbdjUcoNmrArFUBEJFZE7hSRn0Vkr4ikisgOEZkqIkNFJMLtGFXBEJEywNdAceAe4Gpgrpf79nWSl3QRqVGIYYYMj4Qvu2VGPo/dvwDDDXkisj6X8+G5JLgdr1Ju03/+SgUIEakHfA80AGYATwO7gYpAT+B9oDFwv1sxumQDEAOkuR1IAWsHlAauM8Z87eO+1wKbgHjgGmB0wYYW0h4D1mWxfls+jjkSmAB8k49jFDV3AiU8HjcCHgImYT/IeloB/AKMAY76IzilAo0m7EoFABGJAaYAdYBLskjgnhGRdtgkr0gQkZLGmCRjp2M+4nY8haCSc7vXl51EpAJwIfA/oBUwVEQeNy5NWy0ixYBwY0ywnKMfjDHz3Qwg473tZgxuM8Z84/nY6UV/CFhsjPm/bHYLtQ/tSnlNh8QoFRiuBxoCL2TX22qM+csYM9ZznYj0F5FfRSTZWX4VkX6Z93W+fk4UkRYiMsPZdqeIPC8iESIS7dzfIiJHRGSuiDTKdIyhztfTPZ0hABtE5KiILBaRAVm02VtEPhORtSJyWET2i8h0EemaxbaJTox1RORLEdkLHHSey3IMu4gMFpE/neOmOO187CS0ntudLSI/icgBJ46FInJdDjFUEZFPRGSfc9xpItIgq3OSFRFpLiKTRGSP87NcLiL3i0i45/nA9sgCzM746t/LJq7GdrZ8BHwA1AJ6eBy7kXO8F7OJ7xMROeb5cxKRyiIyTkQ2Os9tFZG3RKRipn0zhpU0EZEXRWQz9sNUR+f5K0RksnOcoyKyW0S+EZHm2cRys4j85/ycVorIbR7vs4RM28aJyDMisto59i7ntdTx8ufmNY8YuovIvSKyxmlzpYgM8diulsd5GyIewzg8tjEi8oGI9BCRX0QkGfjO43lff4dbi8gsZ9u9IjLB8zyJyMVOm9dn89qWOT9Dyeb5cLF/BxZm8/yNzvH7O4+jnffFfyJyyPl9XCIiz+X4Q/aRZDGG3WNdYxF5WUS2Ob+zM0WkobPNxc7v/GHnZzgsm+P3FPv3ab/zflwsIjcV5GtQKj+0h12pwHCpc/uWtzuIyC3AG8C/wBOAAYYC34jIjcaYzMeqBvwEfAZ8CfTGjp1OB5pgh52MAcoD9zrHaWSMOZ7pOM9gx12Pc9q8BvhERKKNMR94bDcUKAt8CGwGqmI/mMwUkW7GmJ8zHbcEMAf4FXgYOxQou9c+CJvw/owd4nAYqAGc6+y3y9nuAuxX7NuBF4AkYADwjojUMcY8nOnQxbHjyP/A9vbVBu4AvhWRpsaY9Oxictpr67yGVOy52Q5c4PzMWgBXOZve6cQ6DHgK+5W/t64F5hhj1jsJ805n3QwAY8wKEfkLuFJE7vOMWURKAf2wvcwZP6MawO9AJPAusAaoB9wMdBORtsaYA5li+Bj7M38B+x7IGE5yG/Ybg7ec117XeY2/ikhrY8wqj1gewL7fFmJ/1rHAfTjnzpOIxAG/Yc/xe8AyoDJwCzDPiXGDlz+/OBEpn8X6FGPM4UzrnsL+XryJHYpxM/CBiKw2xvzqxHo19sPTz2T/+9sWuAR4m5Mf1PL6OzwT+Ar7O9wae+7bikg7Y8whYDL2Z38d8I7nziLSETus7uHsvpExxqSLyMfAfc57fmmmTQZjh+p97zx+w4nhQ+AlIByoD3TP5mdRGCYAydjzVQH7d22aiDwKPIv9W/Ue9mfypogsN8b8krGzk8SPx/7ePwmkAL2AcSJS1xhznx9fi1JZM8booosuLi/AHuCgD9uXwf6DWg2U8lhfCptwJQGlPdavxyYDl2U6zgLgOPAtIB7rhzvb9/FYN9RZtwGI81gf56zbC8R4rC+eRdzx2H/2UzOtT3SO/UQW+9Rynhvlse5rbA98RA4/o3Anrv1AFY/1kdgPBelA/SxiuD/Tce7L/LPIoc1fsV/bN/dYJ8DnzjF6ZPHzTPDhvHdw9hnqse4lbPJcxmPdrc52fTPtf52z/mKPdd9ik/5qmbZt67wWz5/7KGf/xKx+9tmc80bYZHesx7qyTsyLgWiP9ZWAA5l/LsArzvYtMh27pvM++MCLn11G7Nkt92Zxbv4GIj3WV3VeyyeZjm2yi8Hj+D0L6Hf4zkzHuctZP8Jj3VPOusaZtn3bOadVsorVY7smzv7PZlpf11n/qse6vWT6fc7LAiSQ6fc8m/NXK4t135H1368koIbH+grYb4Q+8VhX2Vk3MYs2X8H+naib39eniy75XXRIjFKBoRTOEBAv9cL2Br9qjDmxn3P/NWxvdc9M+2wxxnyRad0v2ITyNWOMZ49bRu93/SzaHmc8elyd++OxCUiCx/qUjPsiUkJEymH/+c3DJp5ZeT6b9ZkdwPbInpfdV/tAG5weWWPMVo+4jgHPYYcEZh56cBx4NdO6Wc5tVj+LE5xhCWcBk40xiz3aM9gECuCinI7hheuwvX9feqx7H4gGrvRY9wlwDNsb6mkwNsGa4sQcB5yP7ZU9IiLlMxZsgrga+01MZi8bY04bT5xxzsUq5RxnF/Afp57zXk7M44zH2HdjzHZs7/0Jzvm9CvvNx5ZMMaZge0WzijE7tzrtZ14+z2Lbsc77JSO+LcBKcnkvZOEfY0zmKjR5+R0+iO0tPiVGZ73ne+ttbMJ6YuiXiBQHrsB+u7KVHBhjlmE/zF8lIp55Qsb7aYLHugNAExFpmtMxC9mr2fz9+tYYszFjpbHfKv3HqefvUiAKeNfzveW8v77D/p3ogVIu0yExSgWGg0BJH7av7dwuy+K5jK+wM4/tzaoyxr5snstYXy6LfbIavrE8c5siUhf79XIfbDUUT1l9Hb/LGLM/i/VZeQo4G1uVY4+IzAF+AD4zJy/my8vPaKs5/eLJPc5tVj8LTzm1txz7YSDP461FJBY7nCcRqOTxOeUQNrG+Djs8AWPMXhH5HugnInHGmAPO2N8unJqENsQmJNfhkdxlsjaLdSuzibEV9mLYBGwy6snzPZbxs/ovi8NkXlcB+7PvTRbDZRyZh23l5E/j/UWnWb32PdiefV9k9fPKy/tzrTHmlCopxpijIrLWc1tjzDqxZSqvFpERxphU4HLs35hThsnk4ENsD3NPYLqzbhCwzBizwGO7O7FDgpY4cczGJrrfmdOH0xWWzOcpu79rGc95nr+Ma3VyKusZn8e4lCowmrArFRiWAmc746qzShIyy65XOSc5jb/O7rms2skq2T5lOxEpge0RLQ68DCzBfj19HHiQrMe3HsohvlMDMGaViDTG9nz1ALpiexVHi8jZxpg12cSem5x+RrkdLy/t+SIj4TrPWU4PQKSlMWaR83ACttf1MmySdrUT44eeuzi3/8epvaaeMo/rhizOlTMWfi72w+f/sIl3Cvb98jKnlvDz5WeVse0M7LUA/uTL70VOsnpv5+X9kuW482yO9RbwBbai0FfYD2TbOTn2PDcTsd94DQami0gX7IeCB04JyJhvnQ+DfbG/hz2dtn4WkZ6e31AUouzOkzfnL+P+YLIv7enN32SlCpUm7EoFhq+wPcbXYy/Ay80a57YJ9iI0T42d28L6J9MYO4TCU0YvVUabPYAqwLXGmPc9NxSRJwoiCKencaqzICJ9scnI3dhhD54/o8wK42eUcays2jsD25Odn/auBbZiL4LNLBKbiF8H3O6sm4rtkR7MyYT9X2PMnx77rcYmgZFZDNnw1UXYpPxCY8xszyec4VCePcMZPZ8NOTnkCI91nnZhr0MoVQAxBpK8/A7XFZFIzyRYRKKwvfX/Zto249qE60RkKdAJeCaroUxZMcbsFpGpwEXOB/DB2A/cp5VcNMbsddb/nzOEaQx2voh+2A8NgSzjQujdIfb+UiFGx7ArFRjewfZI3itZlHQDEJE2TlUJsNVeUoDbRaSkxzYlsQlbsrNNYbjZGfuc0WYccBM2qZrjrM7o2crc896b7Mevey2bKh8ZZejKejzeCFwjIhk1zzPqhmdcSPptfmPJYIzZia1kcoHneF4ngXnQeTgpL8cWW1ayC/CVMebLLJaJ2HG7VzoJHM4wiE+AziJyJXbc7im96MaYPdjE/mKngkjmdkUylcnMQXbn/AZO1pzP8BNO1RURifbYthInK+lkxHgcO669vYhcShYkU/lJFyRz8n3nrbz8DpfCVsbxdIuz/hvPlc75/wA7JG2ks/pdH2OcgL1WZBD2m5qfPMe/iy0BWTpTuxkX64LvPxM3fI59L44WOx/GKcSWE43yf1hKnUp72JUKAMaYQyJyPraH+BsRmY79Z70HO4a3G/Yf77PO9vtF5H7smOV5IvKBc6ih2JJ8N5rTS/EVlN1Om+9hk7NrsBd3Xm9sWTmwF7NuB15wvi7fDLTE9vIuAZrlM4bpInIAOwRjE3aM/FBsEv4RnChPdxs2Sf5LRN7CDsu5Als3/CnjUWawgNyB/dDys4hklHU8H3vuJhpjMvekeuta5/arHLb5Cjt2/CLgU2fdBGzFjHFk0zuKLVX4CzBXRD7EJlth2OEP/bA996O8iPEH7NCPj0TkdexY4U7YoRJr8Ph/Y4zZIyKjsdci/Coi/4dNDIdhx3u35dThHw87x/pcRD7HXmh6DDsWuS/2AsmhXsQIcK6InJHF+hRjTJ4+UDnx9BRbqnIjNm/9NKcd8vg7vAYY6XwgXIC9sPpabO965oulwQ4Tuw8YiC0F6uv7/Xvs36BnsB8KMg+bKglsE5HJ2PfNTmxv/83Y8/8dAc4Ys1lEbsZ2mqwQkY+w1aUqYP9O9cd+47HerRiVArSsoy66BNKCTVruwiZQ+7D1vHdg/3FejZ1R0nP7i7C9uinO8hvQP4vjrgcSs1g/ikyl0pz1tTi9lOJQZ11PYDQ2MTmKHX9/ZRbHbg786LyOJOzFkl2wvX4m07aJwPpsfiZZxXID9gPNdmzitg3bU9wti/27OtsexJZv+xv74SLzdlnGkFX7uZzDFtjezr3Oz2cFdnhA5nOX8fNMyOV44dihMDszHyPTdlWxSfn0TOuXOO38lMO+5bGVc1Y6P6P9zn6v4FEaMLv3i8fzZzvv3STnGN8DTXP42d7qtHkUOzThNmzvsgHaZ/G78agT12GnjRXYpLSDF+clI/bsls3enJusXgv224vpznvM4PH+JoeSj3n5HcbWXp/lbLsP+wE1Pofjz3RiuNqb928W+7/m7H8Aj7KtznORwNPAn9jE/qgT53t4lEz1sp0E8l7WsVambWtld6wc3oudsB/ud2L/pmzFXkB7Dx6lR3XRxa1FjMnuGhallDpJRIZiSwh2M8YkuhuNClUi8ho2ca9ijMnuIsAiR+zsuOuNMQk+7jcVOBP788zqAmKlVBDQMexKKaX8znPsuse6ytiLG5dqsp5/IlIPOxzrI03WlQpuOoZdKaWUGxJE5DnsrLWbscMYbsBWmhnhYlxBT0Q6YCs3DccO73jR3YiUUvmlCbtSSik3rMZeRHkDdmKkI8B84Gmj5fXy62bsNxVrgauMMevdDUcplV86hl0ppZRSSqkApmPYlVJKKaWUCmA6JCYH5cuXN7Vq1XKl7ZSUFIoXL+5K2yprek4Ck56XwKPnJDDpeQk8ek4Cj5vnZMGCBbuNMVlOVqcJew5q1arF/PnzXWk7MTGRhIQEV9pWWdNzEpj0vAQePSeBSc9L4NFzEnjcPCcisiG753RIjFJKKaWUUgFME3allFJKKaUCmCbsSimllFJKBTBN2JVSSimllApgmrArpZRSSikVwDRhV0oppZRSKoBpwq6UUkoppVQA0zrsSimllFJZOHjwIDt37iQ1NbVQjh8XF8eKFSsK5dgqbwrynERERBAdHU2FChWIjo7O37EKJCKllFJKqRBy8OBBduzYQdWqVYmJiUFECryNpKQkSpYsWeDHVXlXUOfEGENaWhrJycls3LiR+Ph44uLi8nw8TdiVUkoppTLZuXMnVatWJTY21u1QVBASEYoVK0aZMmWIiopi+/bt+UrYdQy7UkoppVQmqampxMTEuB2GCgExMTEcPXo0X8fQHnalVJF14ACsXm3v160LpUtDcjJERtpFKVW0FcYwGFX0FMT7SHvYlVJF0oED0KwZtG1rlw8+sOs//RQuuACMcTU8pZRS6gTtYVdKFTmbN0NCAmzZAu+8A2XLQteu9rmDB2H6dJg5E3r2dDVMpZRSCtAedqVUEbJgge1JX7gQ4uJg5Ei47jq46CKbtAPceitUqwaPPaa97EoplR8JCQnUqlXL7TBCgibsSqki45134PbboU8fm7w/9tjp20RFwSOPwO+/wyef+D9GpZRS1qJFixg1ahTr1693OxTXacKulCoS0tLg22+hd2+blOfkmmugdm14883stzHGDp/JmE/l+HFISiq4eJVSKthNnz6d//77L8/7L1q0iNGjR2vCjibsSqki4rvvYNs2GDw4920jI2HOHBg1yj7euBFSUk7d5rbb7LCali0hPd0+rlkTli8v6MiVUio4RUZGEpVbD4lLjDEkJye7HYbXNGFXSoU8Y+Dii6FECTjvPO/2qV4dunWDJUtsIv7VVyefW7fO9r43bmwT9B9+sAn7vn12PLyOfVdKBZMPPvgAEWHGjBmMGjWKmjVrEhUVRfPmzfn0009P2/6bb76hU6dOlChRghIlStCpUye+/fbb07bLagx7xrqtW7cycOBAypQpQ/HixenTpw8rV648sd2oUaO45pprAOjWrRsigogwdOjQE9scPXqUp556iiZNmhAdHU3p0qW54IIL+Pvvv09pMzExERHhgw8+4I033qBx48ZER0fz/PPP5+On5l9aJUYpFbK2bIFx4+Dmm+HCC+GSSyDCx796TZtCrVo2Qb/ySrt/5crw8su2/ONZZ9mLWM8/346Rv/56WLbM7qeUUsHkgQceICUlhZtvvhkR4f3332fgwIEcOXLkRKI8duxYbr31Vs444wweeeSRE4lw//79efPNNxk2bFiu7aSkpHD22WfTsWNHnnrqKdatW8crr7xCv379WLp0KeHh4Vx88cVs27aNt956i4ceeohGjRoBULduXcBObHXOOefw22+/cfXVV3Pbbbdx4MAB3n77bTp16sTcuXNp27btKe2+/PLL7NmzhxtuuIFKlSpRvXr1gv0BFiZjjC7ZLG3atDFumT17tmttq6zpOQlMOZ2Xa681JjLSmLVr89fGe+8ZA8ZceqkxX3556nMTJxozfry9v2GD3e7FF/PXXrDT35XApOfFN8uXLy/0Ng4ePFjobXjj/fffN4CpUaOG2b9//4n1+/fvNzVq1DBlypQxhw4dMnv37jXFixc3devWNQcOHDix3YEDB0ydOnVMiRIlzL59+06s79q1q6lZs+YpbXXt2tUA5plnnjll/bPPPmsA8+OPP54WV1bv3RdffPG07TNiqV69uunateuJdbNnzzaAKVOmjNmxY0eOP4vCOifevJ+A+SabnNTVHnYROQd4BQgH3jHGjMn0vDjP9wUOAUONMQtz2ldELgNGAY2A9saY+R7HexC4DkgHhhtjphXqC1RKuWbUKHjvPbjnHnsBaX5ccw1s2gRPPQU//mh71jNmQh048OR2NWpAw4aQj2uslFKBbsGdsG9RgRwqJj0dwsPzf6AyLaHNy/k+zM0330xcXNyJx3Fxcdx000089NBDJCYmkpSUREpKCsOHD6dUqVIntitVqhS33347d911FzNmzODSSy/NsZ2wsDCGDx9+yrru3bsDsGrVKvr06ZNrrP/3f//HGWecQZs2bdi9e/cpz/Xq1YsJEyZw+PBhYmJiTqwfPHgwFStWzPXYgci1hF1EwoE3gF7AZuAvEZlsjPG8ZOtcoL6zdADGAR1y2XcpcDFwSn0HEWkMDACaAFWAGSLSwBiTXogvUynlgm3bYPRom0yPGZP79t547DG44w44fPhksp6VBQugePGCaVMppfwpY9iJp8aNGwOwdu3aExdpNmnS5LTtmjrjANeuXZtrO1WqVCE6OvqUdeXKlQNgz549XsW6YsUKDh8+TIUKFbLdZvfu3acMe2nQoIFXxw5EbvawtwdWG2PWAojIp0A/wDNh7wd86HxN8IeIlBaRykCt7PY1xqxw1mVurx/wqTHmKLBORFY7MfxeSK9PKeWSqVPt7YgRvo9Zz0lcnF1ykpGsHz4MHh07SqlQUQA92RkOJyVRsmTJAjtefmWRO2E8rqL3vJ8f4Tl8q+BtG8YYmjVrxosvvpjtNpmT+djYWO8CDEBuJuxVgU0ejzdje9Fz26aql/tm1d4fWRxLKRViEhLgpZegWTN32n/tNXj6aVtNJkArmiml1GmWL1/OhRdeeMq6FStWAFCnTp0TPezLli2jR48ep+2bsV1ByeoDRIb69euza9cuunfvTlhY6Bc9dDNhz+osZP5Yld023uybl/YQkWHAMID4+HgSExNzOWzhSE5Odq1tlTU9J4HJ87ykpQlLl5aiWbODtGxpmDPHnZgOHSrLtm3Nee65pXTuvDv3HUKM/q4EJj0vvomLiyOpkGdDS09PL/Q2vHHkyBHAVoAZNGjQiXHsBw4cYNy4cZQuXZrWrVtz7NgxihcvziuvvMKll1564tuBpKQkXnnlFUqUKEHHjh1PvKb09HSMMae8xqzWASc+DBw9evTEcxk98Vu2bDlt+yuuuIJHHnmEp59++rTx8AA7d+48MV790KFDJ15nbj/vwjonR44cydfvn5sJ+2bAs55ONWCrl9tEerFvXtrDGPMW8BZA27ZtTUJCQi6HLRyJiYm41bbKmp6TwOR5XhIT4a674OuvbT10t3TuDC+8AIsWNeWRR9yLwy36uxKY9Lz4ZsWKFYU+XCUpQIbEZIwnr1ChAj179uTaa6/FGMP777/Ppk2beOedd4iPjwfg2Wef5dZbb6Vnz54nSj1+8MEHrF27ljfffJNq1aqdOG54eDgicsprzGodQIkSJQCIioo68dzZZ59NWFgYL774IkeOHKF48eLUrl2bDh06cP/99zN37lweeeQRfv31V7p3706pUqXYuHEjM2fOJDo6mtmzZwMnh8JER0fn+vMurHMSHR1Nq1at8ry/mwn7X0B9EakNbMFeEHplpm0mA7c5Y9Q7AAeMMdtEZJcX+2Y2GZgoIi9iLzqtD/xZYK9GKeW6776zF4T26uVuHBER9oLX8ePh999tZZlhw6CqDsJTSgWwZ555hp9//pnXX3+dHTt2UL9+fT7++GOuvPJkinXLLbdQuXJlnnvuOUaPHg1AixYtmDRpEv379y/QeGrUqMF7773HM888w80330xqaipDhgyhQ4cOFCtWjO+//56xY8fy0UcfMXLkSMBe0Nq+fXuGDBlSoLG4zbWE3RiTJiK3AdOwpRnfM8YsE5GbnOfHA1OxJR1XY8s6XpPTvgAichHwGlAB+F5EFhlj+jjH/hx7UWsacKtWiFEqdKSnw+TJ0L27ndHUbVdfDa++aidWAnjjDTumPQA605RSKksRERGMHj36RCKenYsuuoiLvPgaM6shINkNC6lVq1aWF5wOGTIk2+Q7IiKC4cOHZzkkxlNCQkKBXTDrFlfrsBtjpmKTcs914z3uG+BWb/d11k8CJmWzz5PAk/kIWSkVgIyxtdJXr4bHH3c7GqtNGzsspm9fmDIF7rsPpk+3s60qpZRSvgj9y2qVUiHv4EFYs8aWcfScyMhNInD33XDGGXDnnVCmjP0GQCmllPKVqz3sSilVEOLi4Ndf3Y4iexERcN55sHmz25EopZQKRtrDrpQKWhs2wODB7Rk50g6LCWTvvgszZ9r7K1aAU0VNKaVcN3ToUIwxWkUogGnCrpQKSseP2+EvmzfH8Pjj8MorbkeUs8hIe3vkCPTpA61b29lQlVJKqdxowq6UCkr//GNLJg4fvoqPP4bBg92OyDvR0fDEE7aXfdYst6NRSikVDDRhV0oFpW3boFIlOOusPVx5JZQt63ZE3rviCoiNhamn1blSSimlTqcJu1IqKPXtC1u3QsWKR90OxWdRUdCzp03YA33svVJKKfdpwq6UCjrG2EXE7Ujyrm9f2LkTtmxxOxKllFKBThN2pVTQmTQJ6tSBefPcjiTvBg2yCXu1am5HopRSKtBpwq6UChrGwLhxMHw4lCxpZxMNVsWL20UppZTKjSbsSqmgMWcO3HILxMTA22/bCYmC2fLl0LKlfV1KKaVUdjRhV0oFNGNscl6vHjz3nK2usmgRdOjgdmT5V6sWrF9vX59SSgWrhIQEatWqles6XwwdOhRx4UKlpk2bBuQEUpqwK6UC0qpV0KQJtGsHw4ZBt25www0wenToDCWJjYWrr4YvvoCXX7aTQSmllCpco0aN4ptvvnE7DJ8E+RfKSqlQdfvtsGEDVK8OL74Id9wBYSHYxTBiBCxZAnfdBRUrwpVXuh2RUkrl3/Tp0zEBWrd29OjRDBkyhP79+5/23IIFCyhVqpT/g8qFJuxKqYA0YoSdHGngQLcjKVxVq9oZTzt1gs2b7bqrroKhQ6FXL1dDU0qpPIuMjHQ7hDyJiooKyNhDsL9KKRUKEhJCP1nPEBYGv/wC999/ciKlc88N7rKVSqngkJSUxCOPPEKHDh0oX748UVFR1KtXjxEjRnDo0KFTtt23bx833HAD5cuXp3jx4iQkJLBgwYIsj5vVGPY///yToUOH0qBBA2JjYylZsiSdOnVi0qRJ2ca3a9cuBg8eTLly5ShevDg9evTg77//Pm27sWPH0rt3b6pWrUpkZCSVK1dm0KBBrF+//sQ269evPzEufsKECYjIiSVDdmPY//77by677DLi4+OJioqievXqDBw4kDVr1mQbe0HSHnalVEA5cgT+9z+4/nqoXdvtaPwnPNzeisDYsbbX/Y474LffQnMokFIqMGzZsoV33nmHSy65hCuvvJKIiAjmzJnDs88+y99//820adMASE1NpU+fPvz1119cffXVdOzYkUWLFtGzZ0/KlSvnVVuTJk3i33//5fLLL6dmzZrs2bOHCRMmcPHFF/Pxxx9zZRZjAs855xzKli3LqFGj2L59O6+//jpnn302v//+O02bNj2x3fPPP0/Hjh0ZPnw4ZcuWZenSpbzzzjvMmjWLJUuWUK5cOSpUqMBHH33E1VdfTZcuXRg2bJhXcU+ZMoVLLrmE4sWLc/3111OvXj22b9/OtGnTWLp0KXXr1vXqOPlijNElm6VNmzbGLbNnz3atbZU1PSf+MWOGncf0+++92z5Uz8sHH9ifwyuvuB2J70L1nAQ7PS++Wb58eaG3cfDgwUJvIzdHjx41x44dO239I488YgAzb948Y4wxb775pgHMY489dsp2L730kgFMzZo1T1nftWvX09YlJyef1k5KSopp0KCBadSo0SnrhwwZYgBz0UUXmePHj59YP3/+fCMipk+fPrkee8aMGQYwzzzzzCnrATNkyJDTtjfGmBo1apiuXbueEl/58uVNhQoVzObNm0/bPj09PcvjZObN+wmYb7LJSbWHXSkVUBYtsrft27sahusGD4Yvv4T77oP+/aFGDbcjUkplyKrq3+WX23kiDh2Cvn1Pf37oULvs3g2XXmrXpafHnPh27eab4YorYNMmWz0qs3vugQsugP/+gxtvPPW5xMS8vxbP8dppaWkkJSWRnp5Oz549eeKJJ5g3bx7t27fnm2++ITw8nHvuueeU/W+++WZGjhzpVVvFPUp8HTp0iMOHD2OMoXv37owfP56DBw+edsHn/ffff8qQlTZt2tCrVy9mzJhBcnIyJUqUOOXYx48fJykpidTUVFq0aEFcXBzz8jG+cNq0aezevZsxY8ZQtWrV054P89NXoPpFq1IqoCxeDFWqQPnybkfiLhF47z1YsECTdaVU4Ro7dizNmzcnKiqKsmXLUqFChRPjuPft2wfA2rVrqVy58mkJdVRUFHXq1PGqnZ07dzJs2DDi4+MpXrw45cuXp0KFCowfPx6A/fv3n7ZPo0aNTlvXuHFj0tPT2bBhw4l1s2bNIiEhgeLFi1O6dGkqVKhAhQoVOHDgwInXkBerVq0CoFWrVnk+RkHQHnalVED55x9o0cLtKAJDhQp2UUoFlpx6tGNjc36+fPmTzyclHaZkyZKnPF+9es77N2yYvx71zF588UXuueceevfuzfDhw6lSpQqRkZFs2bKFoUOHctyZIMIYk+1ERsaL8o3GGHr37s2KFSsYPnw47dq1Iy4ujvDwcN5//30mTpx4oi1vjuXpr7/+onfv3tSrV48xY8ZQu3ZtYmJiEBEGDBjg9XFzasuNSZw8acKulAoY6em29vo557gdSeDYuRPuvhsuvNB+5a6UUgXpo48+olatWvzwww+nDO/48ccfT9mubt26TJ8+/bRhK0ePHmXdunWUKVMmx3YWL17MP//8w2OPPcbo0aNPee6dd97Jdr8VK1bQsWPH09aFh4dTs2ZNACZOnEh6ejo//PADtT2qFaSkpOSrdx2gYcOGgK0S08vFWrs6JEYpFTDCw+34zkcfdTuSwFGuHKxcCcOHQxbfFiulVL6Eh4cjIqf0WqelpTFmzJhTtuvXrx/p6em88MILp6wfN24cBw8e9KodOL13fOnSpTmWdXz22WdP2WfhwoXMmDGDHj16nBi/nt2xn3rqqSx710uUKMHevXtzjRmgd+/elC9fnhdeeIFt27ad9rw33y4UBO1hV0oFlPBw8LguqcgLD4fx46FdO3j4YXjjDbcjUkqFkksvvZQHH3yQc889l4svvpiDBw8yceJEihUrdsp211xzDW+99RaPP/4469at48wzz+Tvv//miy++oG7duqSlpeXYTqNGjWjSpAnPPvsshw4domHDhqxcuZI333yTpk2bsnDhwiz327BhA3369OHCCy9k27ZtvP7668TExPDcc8+d2Oaiiy7ipZdeom/fvgwbNozIyEh++uknFi9eTPksLojq2LEjM2bM4JlnnqFGjRonhs5kJTY2lnfffZdLL72Upk2bnijruGvXLqZNm8bdd99Nv379cvsx55v2sCulAsbo0dq7npXWreH222HcOJg+3e1olFKh5L777uOpp55i7dq13HHHHbzxxhv07t2bDz/88JTtMpLga6+9lu+//557772XlStX8tNPP1GtWrVc2wkPD+f777/nggsuYMKECdxxxx3MmTOHCRMmcMEFF2S7348//kh8fDwjR47kpZdeok2bNsyZM4fmzZuf2KZTp0589dVXFC9enEcffZRRo0YRExPDnDlzTqlMk2Hs2LF07tyZJ598kiuvvJKBuczSd+GFF/LLL7/QpUsX3n33XW677TbGjx9P1apVadasWa6vvSCIv7ryg1Hbtm3N/PnzXWk7MTExy5m2lHv0nBS++vWhaVPI4dvR0xSV85KUBJ06QefOtpfd5eufclRUzkmw0fPimxUrVmRZoaQgJSUlnXbRqXJXYZ0Tb95PIrLAGNM2q+e0h10pFRB27oTVq21Sqk5XsiRMmwZNmsDx4/Duu+DFsFGllFIhQBN2pZSr0tJg4kQ7oQhAz56uhhPQKleGW2+1teqHDbOTrKSmuh2VUkqpwqYJu1LKVeHhsG8f/PijHerRsqXbEQW+Vq1g7Fj7M/v8c7ejUUopVdg0YVdKuUoEhgyBXbvstN7KO9dfDyVKwG+/uR2JUkqpwqZlHZVSrhkzBg4fttVhnHK6ykvh4dC+Pcyb53YkSimlCpv2sCulXGGMLVO4eLHbkQSvPn2gZk37s1RKKRW6tIddKeWKxYth40YYOdLtSILX/ffb22XL7AWpZcu6G49SocYYgwRyDVUVFAqihLr2sCulXPHZZ3ZYx/nnux1JcDMG/u//YPhwtyNRKrRERETkOnunUt5ITU0lPDw8X8fwuoddRGKBzkAToCJggF3AUuBXY8yhfEWilCoyjh+Hjz+2QzoqVnQ7muAmYj/4fPyxvWj3rLPcjkip0BAdHU1ycjJlypRxOxQV5A4ePJjvyZhyTdhF5FzgJuAcZ/vM3w0ZIE1EfgDGG2N+zFdESqmQd/CgnbHz4ovdjiQ0jBgB778Pd9xhL0IN0+9Olcq3ChUqsHHjRqKiooiJidGhMconxhhSU1M5ePAg+/bto0aNGvk6XrYJu4h0AV4A2gLrgfeA34E1wB5s4l4WqAeciU3op4rIAuAuY8wv+YpMKRWySpe2PcKqYJQoAc88A1dfDV9+CZdf7nZESgW/6Oho4uPj2b59O0ePHi2UNo4cOUJ0dHShHFvlTUGek/DwcEqWLEmNGjWIiorK17Fy6mFPBL4B7jHG/JzDdr8CEwBEpCtwp7OvXtCqlDrN0aOwZg00bux2JKHlyivhoYfstQGasCtVMOLi4oiLiyu04ycmJtKqVatCO77yXaCek5yS6tbGmH98OZgxZg4wR0Ra5C8spVSomjzZJpRz50KXLm5HEzrCwuD776FePbcjUUopVdCyTdh9TdYLal+lVGh77TWoXVsvjiwMzZq5HYFSSqnC4PWlSSIyS0R65PB8NxGZVTBhKaVC0aJF8PPPcNtttrKJKniffAL9++tkSkopFUp8qSWQAMTn8HxFoGu+olFKhbTvvrNlCIcMcTuS0HX4MHz7rR1ypJRSKjQUZPGv0kDhXEatlAoJP/9sh22UK+d2JKFr4ECIjrZJu1JKqdCQYyUXEWkOtPRY1UVEstqnLHALsLzgQlNKhZovvoAtW9yOIrTFxMCZZ0Jioi2d2a0bVKnidlRKKaXyI7fSixcBI537BrjRWbKSBOjk2EqpbMXF2UUVrq5dYdQoePZZGDMG/vlHJ1NSSqlgllvC/gG2proAs4CngJ8ybWOAZGC5MeZIAcenlAoR330HCxfCww9DhM7SUKj69IElS6BFC3jsMfjwQxg61O2olFJK5VWO/zaNMRuADQAicg0wxxiz3g9xKaVCzBdfwIwZMHJk7tuq/OnY0c54aoytGvP++5qwK6VUMPP6S1JjzARN1pVSebVgAbRu7XYURYuIvQj1559h82a3o1FKKZVXPo9qFJG2InKriDwiIo9lWh718VjniMh/IrJaREZk8byIyKvO84tFpHVu+4pIWRH5SURWObdlnPXFRGSCiCwRkRUi8qCvr10plTcpKfDvv9CmjduRFD0DBtjbJUvcjUMppVTeeT2SVERigK+B3tgx7ca5xeO+Af7n5fHCgTeAXsBm4C8RmWyM8aw0cy5Q31k6AOOADrnsOwKYaYwZ4yTyI4AHgMuAKGNMMxGJBZaLyCf6rYFShW/RIjh+XBN2N9Svby86bdYMjh2zw2SiotyOSimllC986WF/DJusPwl0wyboQ7BJ9c/AX0BjH47XHlhtjFlrjDkGfAr0y7RNP+BDY/0BlBaRyrns2w+Y4NyfAPR37huguFOWMgY4Bhz0IV6lVB4lJUFsLLRv73YkRVOzZvb2s89guNbyUkqpoONLwn4p8IUx5jFgqbNuizFmGtATiASG+nC8qsAmj8ebnXXebJPTvvHGmG0Azm1FZ/2XQAqwDdgIPG+M2etDvEqpPDrnHFi/HipVcjuSou3vv+Gdd2DFCrcjUUop5QtfiqtVB1507qc7t5EAxpg0EfkEuBnwdmy4ZLHOeLmNN/tm1h4bdxWgDPCziMwwxqw9pUGRYcAwgPj4eBITE3M5bOFITk52rW2VNT0nvktPh0WLStO69X4kq9/aAqDnxXtnn12Mt97qQNeuaYwZs4Q6dVIKpR09J4FJz0vg0XMSeAL1nPiSsCd5bJ8EHMcmvxkOAL70n23GfgjIUA3Y6uU2kTnsu0NEKhtjtjnDZ3Y6668EfjTGpAI7ReRXoC1wSsJujHkLeAugbdu2JiEhwYeXVHASExNxq22VNT0nvvvmG7j3XpgyBc47r3Da0PPim19+gfPOi2D06HbMnw8VKhR8G3pOApOel8Cj5yTwBOo58WVIzBqgAYAxJh1Yhh0mg4gIcDGnDlPJzV9AfRGpLSKRwABgcqZtJgODnWoxHYEDzjCXnPadjB1bj3P7rXN/I9DdOVZxoCPwrw/xKqV89P33ULq0HRKjAkPLljB5sr2u4Icf3I5GKaWUN3zpYZ8BXCsidzoJ+5vA6yKyBjscpTbwkLcHc4bR3AZMA8KB94wxy0TkJuf58cBUoC+wGjgEXJPTvs6hxwCfi8h12CT9Mmf9G8D72PH3ArxvjFnsw+tXSvlo9mzo2hXCw92ORHlq0wa2b4fISLcjUUop5Q1fEvYxwEc448eNMWNFJBoYhB0b/jbwrC+NG2OmYpNyz3XjPe4b4FZv93XW7wF6ZLE+mZPJu1KqkG3cCGvWaFWSQBUZaUs8Llyo5TaVUirQ+TLTabIx5j9jTJrHuheNMa2NMe2MMc84CbZSSjFnjr3t1s3dOFT2Xn4ZOnSA116DPXtOrjcGJk60tfOVUkq5z+eZTpVSyhtXXgmzZkGTJm5HorJz/fXQqZP9FqR8ebjsMlt+8/vv4aqr4Lrr7GRLSiml3OXLkJiMi0t7YmceLcfp5RWNMcarmU6VUqFn+nSYNMn23EZFae96oCtZEhIT4c8/4euv4fnnoV49eOopGDkSRo+GsDB49123I1VKqaLN64RdROoD3wBnkHUddLAXn2rCrlQRs3EjLF4MF1xgH8+dC337wpgxIXLB6c5fIC0FqvRxO5ICJ2KHxXToAIMG2eEwIjBqFOzYAR98AC+8YKv9KKWUcocvPeyvAXWBB4BZwJ6cN1dKFQU7dkCtWtCli621HhcHn3wCN98cRMl68jqQMChe0z5OPwKr3oSkVRBZBlY8Bxg4ZwGUbupqqIWpWbNTH193HYwfD59/DsOGuROTUkop3xL2zsDLxpjnCysYpVTw+eEH2yv7yiu2xnd6Ojz9NNSo4XZkOdg9D2KqQPHqNjn/qRMc2QVtXobyZ8GvAyBpJRQrBakHoXQzOLwd5t8GPRPdjt5v2rSBhx+Gdu3cjkQppYo2XxL2Y8C6wgpEKRWcpkyBqlWhRQv7ODw8wJP1ZU/BPw/bhL3t67B/MRzeBsVrw9L/QURxm8Qn/GiHwKQmQ0QsLH8W/nkQUjac7IkPcSLwxBNuR6GUUsqXKjHTgE6FFYhSKvgcO2YvND3vPJvcBbx9i+GfR6DK+ZB+GH6+GJaMgjKtoeUYOLIDktdCq+dOjlcvVsIOl6l5uX382yDYNt21l+BvKSnw++9w6JDbkSilVNHlSw/73cBcEbkHeM0Yo8W+lCrifv7ZTnF/3nluR+KlxY/aYS5nfQjH0yF5DRzeCmVaQXQFCI+B8GiofvHp+5aoA2Vawq5fYHYfuOIIhEf5/SX4W2IinH++PdedO7sdjVJKFU2+JOy/AsWxs5mOEZGt2BlOPRljTN2CCk4pFdjatrUXmPY4bW7hAHJoq71wNGUDbJkMzUbZxwDR5U/dtsXTtkc9PDrrY531Mfz3Kqx+E7bPgKrB8kkl75o619guW6YJu1JKucWXhH0jtmyjUkoBtiLMgAFuR5GDIzthyhl2zHnJehAWCfVvzn77M+7I+XhxjaHNq7DhU9j0ZZFI2GvUgBIlbNlOpZRS7vA6YTfGJBRiHEqpIHPNNbZ29003uR1JFoyxF4jungfpKbZ3/cBSqHs9RFfM37HDI+2QmQ2fQNNH7VCZjDaDYiC/b0TsbKjjx0P16jBihNsRKaVU0ePLRadKKQXAunV2Qp09gTobw94FsPwZ2JkItQZB/43Qd6mtClMQmj8OEgF/XAOpSbDwHtuTn7KhYI4fYD75BC691F6AqpRSyv+y7WEXkXBjTOYx6l4RkQhjTFrew1JKBbK33rI9r0OGuB1JNtZNgLAo6DnXTnQUEQuRpQvu+LHVoN04+GMofN8YDm22639oBbE1bHudP7fbhYAyZeCzz9yOQimliq6cethXishgEfF6rkIRiRCRa4GV+Q9NKRWIduyAV1+Fyy6DaoGYj6Yfs8NVqvWH8u1t8lwYag+C7jNt/fbyZ0LvP2y5yOgKsH8pzOppJ10KMT/8AA8+6HYUSilVtOQ0hv1L4E3gWRGZCPwA/GWM2e+5kYiUBToCfYErsJVkXimUaJVSrnv5ZTh6NIAn1Nn6PRzdA3X80P0f3xV6zT35uHwHe7tjDszqDn/eZCvLhNDY9smT4e237YWopUvDrbe6HZFSSoW+bBN2Y8wDIjIeGAFcD9wBICL7gL2AAGWB0s4uycD/Ac8aY0JzIKdSim7doEIFqF/f7UiysW4CRFeCSr3ciyG+KzT/n51R9dh+KN0EWj4bEon7ww/Du+/CI4/Yx126QPPm7saklFKhLscqMcaYdcCNInIvcB5wNtAYqIAt8bgYWAokAj8aY/SSJKVCXO/edglIu/+ALVPgjLshzJeqtYWg8YN2SMzK12HbD1D1Qjt0xu248qlaNXj6adi6Ffr21WRdKaX8wav/HMaYJOBTZ1FKFVHr18P+/TZJCwu0GlOpSfDrAIitDk0ecjsa25vecowt/TipCvxyqR1f32MGlG3jdnT5cs89pz4O0YqWSikVMALtX65SKoC99Ra0aweHD7sdSRYWPwopG+2Y8YKsCJNfEcWh5pV2EieTDol9YV/ozEI0fDgMHux2FEopFdo0YVdKeW3uXGjbFooXdzuSTA6uhJWvQb0bocJZbkdzulbPQLdpcM5fdrbVGV1PloIMcqmp8O23kKaFfJVSqtBowq6U8srhw/Dnn3D22W5H4sEYe7v8aZsINxvlajjZKlYKKveGUg2hx2w4fhSmdYBvqtuKNkGse3dISoIFC9yORCmlQpcm7Eopr/zxh+1NDZiE/fB2mFSFFrvvhnUfQt1hEBPvdlS5K1nPfrA4vNX2su+Y5XZE+ZKQYG9nBffLUEqpgKYJu1LKK3Pn2gsLO3d2OxLHf6/CkR3EHVsCVftBiyfdjsh7je+HS/ZARAnYMdvtaPKlQgVo2hS++ebkFx5KKaUKVnDXF1NK+c1dd9nhD3FxbkcCpCbDqrFQ41J+PXYNXbqcE3xlSqLKQoUusH0G7Ftka8fHVHI7qjy55RZISQm+U6CUUsHC64RdRKobYzYVZjBKqcBVqpSdJCcgrP8YUg9Aw7tIX3Y0eDPFSj3g7x/gh1YgEdD1O6hyjttR+ezmm92OQCmlQpsvQ2LWi8gPInKJiGjPvFJFyNKl8PjjsGuX25EAu363FWFKt4DyHd2OJn/q3QgdP4AuX0FcY/jtKji61+2o8uyFF+DJIBqZpJRSwcKXhP1NoAPwObBVRJ4XkUaFE5ZSKpBMmQIjR7odBbD7T/jpLDiwDBrcFrw96xmKlYA6Q6D6xdBuPBzbC1t/cDuqPFu6FB59FCZPruJ2KEopFVK8TtiNMbcAlYHBwFLgLmCpiPwmIteISGwhxaiUctHBg/D669C+vb3A0FX/vgjF4uC85VD3OpeDKWDl2kNUOdg2ze1I8uyNN6BvX3jppQaMGeN2NEopFTp8qhJjjDlqjPnYGNMdqAc8DVQD3gG2i8hbItK+EOJUSrnk6adh61Z49VWXA0nZBJu+hHrDIK5R8PeuZxYWDpV6wfbpYI67HU2exMbCpEnQpcsuRo6E/fvdjkgppUJDnss6GmPWGWMeAc4APgZKANcDv4vI3yJyWQHFqJRy0bRptjpMhw4uB7L2AzDpUD+Er3Csch4c2QGzetpKOEGoWDEYMGAjxYrBP/+4HY1SSoWGPF88KiLNgeuAq4CywAbgXeAYcCPwqYg0MsY8XhCBKqX8zxho0QJat3Y5kPQjsPY9qNQTStR2OZhCVOtKSNkAix+BbT9CjUvdjihPGjVKYudO2+OulFIq/3zqYReRUiJyk4j8BfwN3AzMAfoCdYwxTxhjngUaAF8CtxZ0wEop/xGB99+H2293MYije+D7JpCy3l5oGsokzE6qFFE8qCdUErHJujFw6JDb0SilVPDzOmEXkQ+BbcBYbI/6I0B1Y8wlxpgfjTk5x50xJh34FnD7EjWlVD6kpATA7JWbvoLktdBlElTr53IwfhBWzE6otDPR7UjyJT0dzjwT7r7b7UiUUir4+dLDfgUwFehjjKlrjHnaGLMjh+1/A67JV3RKKVddfTV06uRyEJsmQYk6RSNZzxDfDQ4sh0Nb3I4kz8LD7VCqd96B335zOxqllApuviTs1Y0xlxljfvJmY2PMemPMhDzGpZRymTHw++9Qt66LQRzbBztmQrWLQq8qTE6q9bc97fNvD4CvOPLu6aehRg0YNAiOHXM7GqWUCl6+1GHfmXFfRGJFpJGz6GVFSoWgTZtg+3bo6OZkoktGw/E0qD3IxSBcUKoBNH8SNk+CXb+4HU2excXZGv7r1sEnn7gdjVJKBS9fLzptLCJTgf3YyZOWAvtFZKqINCmE+JRSLvn9d3t75pkuBbDrd1j5mi3jWKalS0G4qN4wkHDYNt3tSPLl3HOheXOYN8/tSJRSKnh5XdZRRFpiK8KUAH4ClgMCNAZ6A51EpKsxZlHBh6mU8rc//oCYGGjWzIXGj+2HXwdAbE1o8ZQLAQSAyDg7++n2GdDif25Hk2ci8Pbb0ES7dJRSKs98qcP+PHAcaGeMWej5hIi0BmYBzwG9Ci48pZRbLroIqle3E+H43T8PweHN0Ot3m7gWVZV6wrInIWUjFK/hdjR51t6Z/3rHDujRww7L//xzTeKVUspbvgyJ6Qi8njlZB3DWvQG49eW5UqqAnX22SyX5DqyAVeOhwXAo396FAAJItYuAMPiuPqz90O1o8m3xYpuk79sHnTvDr7+6HZFSSgUHXxL2I8D2HJ7fChzOXzhKqUDwyy+w8LSP5n6y6WvAQOMHXAoggJRtBeevsHXZ/xhix/UHsV694LPPbJnHChWgd29YscLtqJRSKvD5krBPBS7M4fkLgR/yF45SKhDcfz/ccINLjW+ZAmXbQUwllwIIMCXrwdnfQHg0bAiNUiu1asHs2VCuHCxb5nY0SikV+HxJ2O8GyonIFyLSTkRKOkt7EfkSO/vpXYUTplLKX7Zvtxec9u/vQuNHdsGeeVD1fBcaD2DFSkDlc+y3D+a429EUiKpVYdUquPRStyNRSqnA58tFpzsBA7QGLs70XMaMJjvl1MlNjDHGlzaUUi777jt7UaArCfv2GYCBKue60HiAq34JbP7GLtUz/wkOTlFRkJ4OKSlQqpTb0SilVODyJZn+EJuwK6VC2OefQ5060LSpC43vmAXFSkOZ1i40HuCqXwz/vgS/DYJzFkBcI7cjyjdjoHZtOO88GDfO7WiUUipweZ2wG2OGFmIcSqkAsG8f/Pkn3HuvrZ/tdztmQXwChIW70HiAi4iFhKkwuTb89zK0f9PtiPJNBBo21EmVlFIqNz7NdKqUCm1lysDmzTB8uAuNJ6+H5LUQ392FxoNETDzUHADrP4bUg25HUyA6dLDlHnfvdjsSpZQKXD4n7CLSTUReFZEpzvKqiHTLS+Mico6I/Cciq0VkRBbPi3P81SKy2JmgKcd9RaSsiPwkIquc2zIezzUXkd9FZJmILBGR6LzErVQoK1kS4tyYq2jHLHurCXvO6t0IaSlO+cvgd+WVdhz766+7HYlSSgUurxN2EQkTkY+AGcBtwDnOchswQ0Q+FPH+S3QRCcdOtnQu0BgYKCKNM212LlDfWYYB47zYdwQw0xhTH5jpPEZEIoD/A24yxjQBEoBUb+NVKtStWQOtWtka2a7YMQuiK0Jc5j8D6hTl2kNMZdjyvduRFIjGjaFfP3jvPTgeGgVwlFKqwPnSw34PcBXwJdAKiHGWlsDnznO+zIvYHlhtjFlrjDkGfAr0y7RNP+BDY/0BlBaRyrns2w+Y4NyfAPR37vcGFhtj/gEwxuwxxqT7EK9SIW3aNFi0yE5o43fGOOPXu7s0eD6IiECVvrB9OhwPjT6HZ56B116DMB2kqZRSWfKlSsxQYLox5opM6xdje7jLANcCL3h5vKrAJo/Hm4EOXmxTNZd9440x2wCMMdtEpKKzvgFgRGQaUAH41BjzbOagRGQYtjef+Ph4EhMTvXw5BSs5Odm1tlXWQv2cfPRRMypXjmXz5nls2eLftmNTN9L+8Db+21+VbT7+jEP9vGSl/OGaNE09yKKf3mB/VEu3wzlNXs5JXJydTOm77yoTF5dG1667Cie4Iqwo/q4EOj0ngSdQz4kvCXsdYGwOz38HPO/D8bLqRstcNjK7bbzZN7MIoDPQDjgEzBSRBcaYmaccxJi3gLcA2rZtaxISEnI5bOFITEzErbZV1kL5nOzeDfPnw513QrduCf4PYOVY2AUNz76ZhiXr+rRrKJ+XbKW2hi9H0bLyXmie4HY0p8nrOUlLg4cegqVLbWnRQYP0C5eCVCR/VwKcnpPAE6jnxJcvIFOA+Byer+Rs463NQHWPx9WArV5uk9O+O5xhMzi3Oz2ONccYs9sYcwiYip0ESqki77PPbLJ09dUuBbBjFsTWgBJ1XAogyBQrBWXawI7ZbkdSoCIi4IsvoHlzGDwYbrlFx7UrpRT4lrD/DNwmIk0yP+Fc8HkrMNeH4/0F1BeR2iISCQwAJmfaZjIw2KkW0xE44Ax3yWnfycAQ5/4Q4Fvn/jSguYjEOhegdgWW+xCvUiGrWjW47z6bKPmdOW4Tz0o6ft0n8QmwZx6kHXI7kgJVrRr8/LN9P44fDy+95HZESinlPl+GxDwG/AH8LSLfcjLZbQJcABwDRnp7MGNMmojchk2kw4H3jDHLROQm5/nx2F7wvsBq7DCWa3La1zn0GOBzEbkO2Ahc5uyzT0RexCb7BphqjAmNMgtK5VO/fnZxxf7FcGyvlnP0VXw3WPEcbJ4MtQa4HU2BCguzF6Ju3QqlSrkdjVJKuc+XmU6XiEhX4BXgEmfJ8BtwhzFmiS+NG2OmYpNyz3XjPe4bbM+9V/s66/cAPbLZ5/+wpR2VUo79++1Ss6ZLHdzbM+qv52k6h6IrvhuUaQnzroO4RlCmhdsRFSgR+D/9a62UUoCPEycZY+YbYzphx7J3BM7EVmXpbIxZUBgBKqUK1+TJULs2rFjhUgA7ZkHJBhBbzaUAglR4NCT8AMVKwJ832aFFIejQIfjpJ7ejUEopd3mVsItICRFJF5FHAYwxu4wxfxpj5hljtPaWUkHsr7+geHFo2NCFxo+nws45Ohwmr2IqQcvnYM8fsGmS29EUirFjoXdvWL/e7UiUUso9XiXsxphkYD8nK64opUJAejr88AN06ADh4S4EsHMOpCVDpSxHsSlv1LoKIsvClu/cjqRQnHuuvZ05M+ftlFIqlPkyJGY2trKKUipEfP45rFljy+f5XWqSHcoRWwMq93EhgBARFg6Ve8O2aXbG2BDTuDFUrGgnVVJKqaLKl4T9PqCziIwWEb1uX6kQMHYsnHEGXHSRC40vfwaS18BZH0Oxki4EEEIqnwNHttuKOyFGBLp1swl7CH4eUUopr/iSsM8EooFHgH0isl1E1mZa1hROmEqpwjBlCnzzjS2j51eHd8B/L0PNAVCxs58bD0GVe9vbbT+6G0ch6dYNtm2DDRvcjkQppdzhSx32jdj65UqpEBEXZxe/W/mqnfCn2WgXGg9BMZWhdAvY+iM0fsDtaArcgAF2aEytWm5HopRS7vClDntCIcahlPKzsWPhyBG4+24/N5x2CFaNh2r9oFQDPzcewqqcAytesNcGhNgQo7g46NLF3t+wwc4ZoJRSRYnXX4SLyNkiUiGH58uLyNkFE5ZSqjClpdmZJF2pb73xSzuzacM7XWg8hFU+B0wa7AjdqzO//hrq1IF777XvYaWUKip8rRLTK4fnezjbKKUC3LffwsaNcOONLjS+eRLEVIWK+vm+QJU/E8KiYNfPbkdSaM45B4YOhRdegA8/dDsapZTyH18S9twmLQ8HQnOqPaVCzNixdnbTCy7wc8Nph2z5wWr9bfkPVXDCo6Bsa9j9u9uRFJrYWHjnHWjeHF58UavGKKWKDl9rQ+T05/EsYHc+YlFK+UFSEsydC1dc4cJkSavfgvTDUL2/nxsuIsqfCXvmQ/oxtyMpNCL2uotly2D6dLejUUop/8gxYReROzJKNjqrXs6ilONaEdkH3AxMKfSIlVL5smcP9OlzcgZJv9n6Iyy8CyqfCxW7+bnxIqL8mXD8KOyd73YkhWrgQOjd24VypEop5ZLcqsTsBzIq39YC9gA7Mm1jgKXAH8DLBReaUqow1Kpl66/73fIxULwmnD3Jzs6pCl6FzhAeDbN6QNcpUKmH2xEVishImDbN3jdGR1cppUJfjgm7MWYCMAFARNYBI4wxk/0RmFKqcOzZA+XK+bnRff/AzjnQ6jk71loVjphK0PsPSOwL/74Usgl7hvR0GDIE+vWDyy5zOxqllCo8Xn+haIyprcm6UsFt+3aoUAHeftvPDW+aBBIGda71c8NFUJkWUOsqe3HvkdC+rOjoUVuXfeBAW/JRKaVCVZ5GAIpIrIhUF5EamZeCDlApVXBmzrRDCFq39nPDO+dA6ZYQVdbPDRdRta6yNdk3feF2JIUqNhamToX27e1F1L+HboEcpVQR58vESWEiMkJEtgBJwHpgXRaLUipAzZgBZctCq1Z+bDT9KOz5Ayp29WOjRVzp5hDXBNZ/7HYkha5kSfjhB4iPhxtugONaXFgpFYJyu+jU0xjgXmAZ8BX2AlSlVJBITbW9kb16+bm6xp6/IP2ITpTkTyJQ60r452FIXg8larkdUaGKi4Onn4ZHH7Xv8yi9TEIpFWJ8SdgHAT8aY/oWVjBKqcIzdSrs3AmDBvm54Q2fQlikJuz+VtNJ2Dd9DY3udjuaQjdoEFSpYj+rvPsudO0K9eq5HZVSShUMX/rZygDfFlYgSqnC1aULjBtnp3f3myO7Ye17UGuQjl/3txK1oHgt2DPP7Uj8QgR69IB9++COO2DECLcjUkqpguNLwr4EqFxYgSilClfZsnDTTRDhy/dq+bXuQzuzaaN7/NioOqFcO9jzp9tR+FV8PNx/P3z1FXz5pdvRKKVUwfAlYR8N3CQi1QsrGKVU4Zg3D958Ew4f9nPDW6dCXGO7KP8r2w5S1sORXW5H4lf33GMrx1x+OUzWYsRKqRDgS19bG+ysp8tFZBK2Ikx6pm2MMeZ/BRWcUqpgvPsufPYZXHONHxtNTYZdP0OD2/3YqDpFuXb2du98qHKuu7H4UfHikJgIHTrAnXdC375+/mZJKaUKmC9/wkZ53M/usjUDaMKuVABJT4fvvoM+feyU7n6zYxYcP1akEsWAU7YNILZSTxE7DzEx8M47dmy7JutKqWDny5+x2oUWhVKq0MyaZWc4veIKPze8IxHCo6FCZz83rE4oVhJKnVHkxrFnaN/+5P30dAgPdy8WpZTKD6/HsBtjNnizFGawSinfTZwIpUrBeef5ueFdP0O5DhCuRbFdVa497P3LTnFbBO3ebavHRETAx6E/j5RSKkQV2PQpIhIrInUK6nhKqYKxeTNcfDFER/ux0dRk2Pe39q4HgnLt4MhOOLTJ7UhcUbIk1KoF1arByy+7HY1SSuVNjgm7iBwTkQEej0uKyGQRaZbF5hcBqwo6QKVU/vz0k73o1K/2zAOTrgl7ICjrXHi65y9343BJVJR9/99/P8yfDwsXuh2RUkr5Lrce9ohM20QC5wMVCi0ipVSBCyuw79K8cDwdlj0N4TFQ4Sw/NqyyVKaFvZZg5xy3I3HVoEFQurROqKSUCk7+/DeulPKzK66AYcP83OjqN2HHTGj7BhQr5efG1WnCoyC+O2yZUmTHsQOUKQNffw1vv+12JEop5TtN2JUKUceOwdSpLlTG2PgZlG4Gdf1Z9F3lqOoFkLIODix3OxJXdesGNWvaCcS+/rpIf35RSgUZTdiVClG//ALJyXbSGL85sgt2/QLVLvJjoypXVc+3t1t02k+ADz6ASy6Bxx93OxKllPKOJuxKhaipU+1ESd27+7HRzd+COQ7V+vuxUZWr2Gq2xObGz92OJCAMGwYDB8KTT8KKFW5Ho5RSufNm4qS+IlLJuR+Lnc30MhFpmWm7NgUZmFIq74yBKVOga1c7TbvfGl09Hko1gjIt/dSo8lrNAbDwLjj4H5Rq6HY0rgoPh1desR9qH3kEvvrK7YiUUipn3iTsVzqLpxuz2VZHBCoVANLSYPBgaNDAj43u/gP2LoB2Y+188Cqw1LgMFt4NaydAy6fcjsZ1FSrA8OHwv//B0qXQtKnbESmlVPZyS9i7+SUKpVSBKlYMHnrIz42ufBWKxUGtq/3csPJKbFWocSmsfB0a3QtRZd2OyHV33AHz5kFcnNuRKKVUznJM2I0xRbtwr1JBKD0dvv/eDofxWyJyaCts/BIa3A7FSvipUeWzpo/Bxi9s6c0mD7odjevKlYNp0+x9Y/SLIaVU4NKLTpUKMX//Df362fG5frPmXTuzaYNb/dio8lnpprbk5s65bkcSUA4cgAsv9PPvjFJK+UATdqVCzKxZ9rabPwe0bf7azmpasq4fG1V5Uq4j7Jlnq/koAKKiYO1aO6b92DG3o1FKqdNpwq5UiJk9Gxo1gkqVct+2QKRsgH2LtJRjsCjfEY7tg6RVbkcSMKKj4bnnYM0aePBBOK6fZZRSAUYTdqVCyOHDMHeuC7XXAar282OjKs/Kd7S3u393N44Ac+65cMMN8OKL0Lo1/PGH2xEppdRJmrArFUISE+HQIejf34+Nbp8JJepCqfp+bFTlWakzIKIE7F3odiQBRQTefBM++ghatIAaNdyOSCmlTvKmDrtSKkiccw78848dEuMX5jjs+hmqX+ynBlW+SZidOOngf25HEnBEYNAguwCsXAl79sCZZ7obl1JK5amHXUTqiUgnEdHqtUoFEBFo3tzWYfeL/YvteOiKXf3UoCoQJRtA0kq3owh4r7xih5dNn+52JEqpos6nhF1EzheRNcB/wFygjbO+ooisFpFLCyFGpZQX5syBa66BHTv82OgOZ6oGTdiDS6mG9mLh9CNuRxLQRo+2swUPGgS7drkdjVKqKPM6YReRBGASsBcYDZyYYsIYsxNYAwwo2PCUUt765BP44gsoVcqPje6cA8VrQ3Ed8BtUSjYADCStdjuSgFa+PEycaOu036pTDCilXORLD/tjwD9AB+CNLJ7/HWhdEEEppXyTng6TJsF550FMjJ8aNcdtwh6vvetBp1QDe6vDYnLVpAmMGmU/DH/+udvRKKWKKl8S9rbAx8ZkO9vGZsCnys8ico6I/OcMpxmRxfMiIq86zy8Wkda57SsiZUXkJxFZ5dyWyXTMGiKSLCL3+hKrUoHst99g50642J/Xfh5YBsf26nCYYFTSSdj1wlOv3HefTdrbtHE7EqVUUeVLwh4OHM3h+fKA13PEiUg4tqf+XKAxMFBEGmfa7FygvrMMA8Z5se8IYKYxpj4w03ns6SXgB2/jVCoYfP21na2xb18/Nqrj14NXsZJ2HPu2aW5HEhQiImDkSKirE/kqpVziS8K+AuiSw/PnY4fMeKs9sNoYs9YYcwz4FMg880o/4ENj/QGUFpHKuezbD5jg3J8A9M84mIj0B9YCy3yIU6mAV6ECDBkCJUv6sdFtPzrj12v5sVFVYGoNskOakte5HUnQ2LQJunWD1Tr0XynlZ77UYX8XeFVEZgCTnXVGRGKBMcCZwGAfjlcV2OTxeDN2fHxu21TNZd94Y8w2AGPMNhGpCCAixYEHgF5AtsNhRGQYtjef+Ph4EhMTfXhJBSc5Odm1tlXWAvmcnHWWXfwVXtjxw3Ta/hPbip/P6jlz/NNoNgL5vASyqLT6dETYOOsR1pW6oUCPHarnZMOGWBYubEXPnscYP34B0dHZjRANTKF6XoKZnpPAE6jnxOuE3RgzTkQ6AW8DLwAG+AQohx0u874x5mMf2pYs1hkvt/Fm38xGAy8ZY5JFstrdOYgxbwFvAbRt29YkJCTkctjCkZiYiFttq6wF6jnZtAmqVIHwcD82uvlb2H6Mah1uplqlBD82fLpAPS9B4ddvqLnpK2omjCrQmWpD+ZxUrQq9exfjq6/O5u233Y7GN6F8XoKVnpPAE6jnxKc67MaYQcAl2LHh/2JLPE4FLjPGXOdj25uB6h6PqwFbvdwmp313OMNmcG53Ous7AM+KyHrgTuAhEbnNx5iVCijGQM+eMMDfBVW3TIGIklDhbD83rApU6xchPAr+utm+mVSuevaEBx+Ed96xJR+VUsoffJ7p1BgzyRhziTGmiTGmsTGmnzHmqzy0/RdQX0Rqi0gktob75EzbTAYGO9ViOgIHnOEuOe07GRji3B8CfOvE3cUYU8sYUwt4GXjKGPN6HuJWKmAsX26nT+/e3Y+NGgNbf4DKvSA80o8NqwIXUxlaPAU7ZsLq8W5HEzRGj4bOnWH8eP2co5TyD1/GsBcoY0ya08M9DTuk5j1jzDIRucl5fjy2974vsBo4BFyT077OoccAn4vIdcBG4DI/viyl/Orrr0EE+vf3Y6MHlsHhLVD5XD82qgpNvZtg0yT46xZAoP5NbkcU8CIi4NtvoVgx+/unlFKFzeuEXUQey2UTAxzGJsmJzuynOe9gzFRsUu65brzHfQNkOb9cVvs66/cAPXJpd1RusSkVDL7+Gs48EypX9mOjW52qqFXO8WOjqtCEhUO3H2Bae1j/f5qwe6lsWXt74IBN2v06w7BSqsjxpYd9FCcv7Mzcp5B5faqIPG+MeTgfsSmlcvDvv7BoEbzwgp8b3vYjxDWF2Gp+blgVmrBiEN8DVr4O6cd0qJOXtm2DmjXh+edh+HC3o1FKhTJfxrA3BRYCvwNXAC2dZQDwBzAf6IgdgjIfGCEiNxZgrEopD/XqwY8/wqBBfmw0NQl2/QxVdDhMyCl/Jhw/Cvv+djuSoFG5MjRoAF9+6XYkSqlQ50vCfgNwBOhqjPnCGLPYWT4HugKpwADnAtSuwBJAE3alCklEBPTpAxUr+rHRHbPgeKom7KGo/Jn2dvfv7sYRZK66Cn7+WSdTUkoVLl8S9gHA58aY9MxPGGPSgM+BgZkeNyyIIJVSp5o+He6/H5KS/Nzwlu8hogSU7+TnhlWhi60CJerB2vftsBjllcGDISyMoKvJrpQKLr4k7HHO4u3zu8l9MiOlVB489xx88gnExPix0eNpsPlrqHq+jnEOVa2fh/2L4V9/XxgRvKpWhUsvhTffhL173Y5GKRWqfEnY/wFuEZGamZ8QkVrALcAij9UNgW35CU4pdbpff4UZM+Dmm+2wGL/ZMRuO7oEal/uxUeVX1fpBxa6w8XO3IwkqEybArFknK8copVRB8+Xf/Qhs3fMVIvINsNJZ3xDoh03+BwKISBRwFTClwCJVSnH8ONx5p+3Vu+MOPze+5m07HKaylnMMaZV6weJH7IezqHJuRxMUoqOhdWs7idLhwxAb63ZESqlQ43XCboyZIyI9gRex49k9zQfuNcbMdbY96vTEpxZYpEoppk6F+fPhgw+geHE/NrzzF9j4BTR9DCL8OQ5H+V28M23ujkSocYmroQQTY6BNG2je3P5+KqVUQfJlSAzGmF+MMe2BSsCZwFlAJWNM+4xk3WPbo8aY4wUXqlKqbl247Ta48ko/N/zv8xBdCRrf7+eGld+Va2u/Sdk+3e1IgooIdOpkh8e8957b0SilQo1PCXsGY8xOY8w8Y8wf3sxoqpQqGI0awWuv2SnR/Sb9KGyfAdUvggh/dusrV4QVs2PZ139i6+4rr73wAvTuDTfeaMe0K6VUQclTwi4iJUSkmojUyLwUdIBKKTt2/a67YOXK3LctcLt+hrQUqNLXhcaVKxoMh7QkW+JReS0yEj7/3E6mdOmlsGqV2xEppUKFTwm7iAwQkaXAAWADsC6LRSlVwGbOhJdfhj//9HPDxsDaDyAsCuK7+blx5Zry7aFCF1gyGg5vdzuaoBIXB1OmQNu29mJUpZQqCF4n7CLSH5iIvVD1TUCAT4AvsBeXLgQeL/gQlVJvv21Lxl12mZ8b/vcFWP8xNLxDh8MUNe3fgvRD8PsQW4Nfea12bTu5WfXqkJIC/ftrjXalVP740sN+L7ACaAk85qx7zxgzAGgLNODUOuxKqQKwbh18842dUTEqyo8NH9sHSx6HqhdAyzF+bFgFhLgzoM1r9uLTxY+4HU3QWrMGfvgBevWCNP3co5TKI18S9ubABGPMESCj+ks4gDFmKfAW8GDBhqdU0WYM3HSTTdTvvtvPjf/3mh3H3PwJWwJDFT31roe618OK52DPfLejCUrNm9vKMQsXwr332vuHD7sdlVIq2PiSsIcDe5z7GX9u4jye/w9oWhBBKaWso0ehVSt49VX79bpfbf0eKnSGMs393LAKKK2eh+h4WOjvT4yh44orbMnHV16xZVkz7NwJ2/USAaWUF3xJ2DcDNQGMMYeBndihMBkaAikFF5pSKjoaxoyBa67xc8PpR2HfIih/pp8bVgEnMg7OuNdWC9q3yO1ogpKIHdb2xx+wZAnExEBSEnTtCpUrw/nn63AZpVTOfEnYfwN6ejyeDNwhIo+JyCjgViCx4EJTqmgzBmbPhlQ35gve9w8cPwbl2rvQuAo4da+B8FhY8bzbkQSt8uWhQweoVcs+3r8f7rgD7rsPvv8ePvzQzeiUUoHOl4R9LJAoIhnzkj+MHQYzCnsR6hrshalKqQLwySfQvTt8/LELje9x6kdqwq4AIsvYSkHrP4ZNX7sdTUioXt1en/LMMzaRf+QR2LbN7aiUUoHK64TdGPOXMeYhZzgMxphdxpiW2KoxzYAWxphNhRKlUkXMrl1w661w5pkwaJALAez+3Y5bjvX3wHkVsJqNgtItYPFjuW6qvCcC48ZBRARs0v+gSqlseJWwi0hxZ+hLn8zPGWMWG2OWGWOOZ7WvUsp3L7wABw7Au+/af+R+lX4Mtk6FSr21Oow6KTwSag+GA8sgZYPb0YSUVq3srKjt29shcDpDqlIqM68SdmNMCvAQoN1tShWyrVvh9ddhwABo1MiFALbPgNT9UPNyFxpXAa3qefZ2y/fuxhGCMuZYeP99aN0aduxwNx6lVGDxZQz7GqBSYQWilLJ27rQXpo0e7VIAm76EYqWgUi+XAlABq2QDKFEXtnzndiQh68wz7eyor77qdiRKqUDi60WnN4hIucIKRikFLVvC4sVQv74LjRsD23+yw2HC/TmtqgoKIlDjMvseOaxdwIWhWTO4+GJ44w29CFUpdZIvCXsSsBf4T0ReEJGbRGRw5qWQ4lQq5H3+OVx2mS33FubLb2ZBSl4DhzZDpe4uBaACXu2rwaTDholuRxKynngCjh2DwYPtZ2illPLlcrYPPO7flc02BtBqskrlwccfw99/Q1xc7tsWmh2z7G3Fbi4GoQJaXGMo2w7WvAMN79QLkwvBGWfAiy/CqFGwfj3Uru12REopt/mSsOt/cKUKSVISTJsGN9/scv6zfSZEV4JSDV0MQgW8BrfBH0Ps0JjKvd2OJiQNGwYDB7r8AV4pFTC8TtiNMXMKMxClirLPPoOjR+FyNwuzpB2Crd9DzYHaa6pyVnMA/DMCFo+Eil31eodCEBZmk/XUVNi82fayHzoEyclQsaLb0Sml/C1PI2VFJEpEqopIZEEHpFRR9NZb0LQpdOzoYhBbpkBaik3YlcpJeCS0ehH2/AGzesGO2W5HFLIuugj69LGzopYqBfHxMHmy21EppfzNp4RdRFqLyCzsBagbgc7O+ooiMlNEehZCjEqFtLQ0OPdcuO8+lzu2N0yEmMq2x1Sp3NQaAO3GQfJqmHsxpB9xO6KQNGwYrFtnP9Rffz2MHAm9nVFIW7fqRalKFRVeD4kRkZbAz8Bu7IWl12Q8Z4zZKSIxwBBgRgHHqFRIi4hwseZ6hkNbbQ/7GfdAWLjLwaigUf8mKFkfZvWETZOgln47U9AuvNDOepySAhUqnFyfkgJnnWWHyiQl2Z74G2+E8uXdi1UpVXh86WF/HNgKNAFGAJn7AmcC7QsoLqWKhB9/hK++cjsKYO17tlRfvRvcjkQFm/huULwWLP0f7P7T7WhCUmzsqck6QEwMXHstLF1qHz/yiJ1w7Z9//B6eUsoPfEnYuwBvG2OSseUbM9sIVCmQqJQqAtavh6uvhv/9zw6LcU36UVg11s5sWrKei4GooCRh0OZVOLYHZibA3gVuR1QkhIXBY4/Brl0wf75N1IcPhxo14Phxt6NTShU0XxL2aOBADs+XymcsShUZe/ZAt26Qng6ffmqHxbhm3YdweBs0vt/FIFRQq3YBnLsYoivC3IsIP57idkRFTvPm8NRTUKYMPPQQtGgBc7S2m1Ihw5eEfQ3QJofnuwPL8xeOUkXDO+/YHvapU+0kKa5aNQ7KtIb4Hi4HooJaTDx0+gwOb6HOwbfcjqZI69gRDh+GHj3g9tth1Sq3I1JK5ZcvCftE4OpMlWAMgIjcA5wDfFSAsSkVkoyBd9+Fs892uYwjwKHNsO9vqHmF1l5X+Ve+A9S9gcqHfoDUg25HU2T172+HyVx/PYwbBw0awJtvuh2VUio/fEnYnwf+AKYBc7HJ+ksisgV4FvgJGFvgESoVYoyBV1+Fxx93OxJgy/f2tur57sahQkftwYSRaqsOKdeUKgXjx8OGDbaGezedq1ypoOZ1wm6MOQb0Au4FDgNHgAbYMo/3A+cbY/RSF6VyERYG55wDXQOh3PmWKVC8NpRq5HYkKlSU78jRsPKw8Qu3I1FA1aone9lTU4VLLoGffjr5vF6gqlRw8GniJGNMmjHmJWNMW2NMcWNMrDGmhTHmBWOMm3UulAoK06fDXXfZusmuO54GO+dA5d46HEYVHAljR0wP2Pwt7PzF7WiUh5SUCNasgb59oWFDWyoyOjpASssqpXLky8RJFwLfG2PSCzEepUKWMfDww7B7NzzzjNvRAPsWQVoSVExwOxIVYjaUHEwN5sGsHlChEzR/Asq2gfAot0Mr0kqXTmXOHBg1ys6SWrasrSrTr5993hj97K5UoPKlmNw3wC4RmQh8aIz5u3BCUio0vf66vRDsvfcgMtLtaLC96wDxgTA2R4WS9LBY6DYNVr8F6ybAT50gsiy0eQVqD3I7vCItLg5eeun09cuXw5lnQufOsGmTHbZ3661Qs6b/Y1RKnc6XITE3A6uBO4D5IrJYRO4RkUqFE5pSoWPGDLjjDjvN+ODBbkfj2DEbSjaAmMpuR6JCUakG0Pp5OP9fOGsixDWC36+GBXfa4VgqoHz6KVx0ESQmQkoKPPccPP+821EppTL4ctHpm8aYTkA94AkgFngO2CQi34vI5SKi33cqlYkxdtx6/frwyScQHu52RMDRPbB9OlTp63YkKtRFlYNaA6HHHGh4J/z3CsztD8dT3Y5MeXj8cfjgA0hOhtWr4c8/4cYb7XMzZ9pStEop9/g8v6IxZi0wEhgpIl2AwcCl2DrsB4CyBRqhUkFOBKZMgb17ITbW7Wgc6yfahKnONW5HooqKsHBo8xKUqAsLbof5t0HtIVCuHYQVczs65cgYw96u3cl1EybAxx/DwYO2utVvv9k5JNq2dSdGpYoin6rEZGaM+Rk7RGYEkATEFURQSoWCvXth+HDbY1WzJrRq5XZEHtZ+AGVaQZnmbkeiipqGt0G9m+z49p86wVcV4K9bQasCB6yxY+1Eb3ffDW3a2NlTN2ywz23ZoqUhlfKHPCfsItJTRD4EdmAnTEoF3iiowJQKZsbYWQbHj4clS9yOJpN9i2HfQqgz1O1IVFHVbiyctwI6fwmVesCqsfaaChWQSpSA2bPhn3/giy9g6VK45BJITbUXqYaH29sJE2DqVNi1y+2IlQo9Pg2JEZHG2CEwVwFVgDRgKjABW/JRByUqha1rPGkSPPusrbwQUNZ+YIcg1LzS7UhUUSUCcWfYpUpf2F4Z1k6wybsKWM2b2yWDCDzxBCxbZqtfDR1q148caUtHzpwJX34Jl1+uM60qlV9e97CLyHxgCXZW0+3AnUAVY8xFxphvNFlXyjpyBO67D5o1s18hB5T9S2D1OKh2EUSXdzsapSAiBmoOgE1fwqGtbkejfBARAVddBU89ZYfIrFwJP/8M115rn1+1yo59794dbrrJDg9USuWNL0NiKgHPA02NMe2MMa8ZY/Z4buBrlRgROUdE/hOR1SIyIovnRURedZ5fLCKtc9tXRMqKyE8issq5LeOs7yUiC0RkiXPb3ZdYlfLWyy/D+vXwwgsBUhEmgzHw+1AoFmfrYSsVKM64x97Ou96+T1XQiYqylbA6d4YaNey6m26CnTttB8Zbb9kLWdescTdOpYKVLwl7DWPMA8aY5ZmfEJE2IjIW8Lp7RETCsWPezwUaAwOdITeezgXqO8swYJwX+44AZhpj6gMznccAu4ELjDHNgCHAR97GqpQvrr0WXn0VevVyO5JMds61Y9ebPwExOn2CCiCl6kOLp2HbD7B1qtvRqAIUHW2HBs6caZP36dPdjkip4ORLHfZTrgN3erKHi8g/wJ/ATYAvl5q0B1YbY9YaY44BnwL9Mm3TDzurqjHG/AGUFpHKuezbDzumHue2vxP/38aYjA8Uy4BorRuvCooxto7x+vVQsaKtohBw/n0RIstALR27rgJQg1tsycd/HoL0I25HowpYt26waBEMG2Yfjxlj70+fDvv3a6UZpXLjcx12EekDXAtcCEQCK4HRwFfGmGU+HKoqsMnj8WaggxfbVM1l33hjzDYAY8w2EamYRduXAH8bY45mfkJEhmF784mPjycxMdHb11OgkpOTXWtbZS2nczJ1aiWee+4Mtm5dw4ABm7Lcxk2VU6bQ8MBk1pa8jo2//Ol2OAVKf1cCT17PSfliV9N03yiSvmrGphID2BnT/WRhcJVvgfC7smaN7eD48896TJlShbffPtlv2K7dXp59drGL0flfIJwTdapAPSdeJewiUhu4BjuUpBq2J/1L4ErgYWPM13loO6u/wpkHL2a3jTf7Zt2oSBPgGaB3Vs8bY94C3gJo27atSUhI8OawBS4xMRG32lZZy+6cLF8Ob7xhL6waO7YuYWF1/R9cTvYthmmvQ6Xe1El4kzphgTSwPv/0dyXw5P2cJMCmFpRceBeN9z9B41olofH9BRxd0RVIvyvdusGhQzBnjv0bmpQE1aqVJSEhgeRkaN0aqlWDhx6CHj1C93NbIJ0TZQXqOckxYReRK4HrgK7YEo7fA7c7t7Wx5R3zajNQ3eNxNU4fA5/dNpE57LtDRCo7veuVgZ0er6caMAkYbIzRS19Uvs2fDxdcYOsUf/ghhOVrKrJCkJYCvw6wQ2HO+sjONqlUIKveH6pdCL9dBYtGQHx3KKdTaoai2Fg491y7eNq2zV68OmOGvRYoIcEOoemQ+Tt4pYqQ3NKL/wNqcrKE4yXGmMnGmHS87NHOwV9AfRGpLSKRwABgcqZtJgODnWoxHYEDznCXnPadjP0mAOf2WwARKY39oPGgMebXfMauFGBrDYeH20lFqlZ1O5osLLwHDv4LZ/0fRGc1OkypACRh0P5NiCqrY9qLoPr1bV33VavglVdsD3znzvYaIaWKqtwS9mNALeyFnOeKSExBNWyMSQNuA6YBK4DPjTHLROQmEbnJ2WwqsBZYDbwN3JLTvs4+Y4BeIrIK6OU8xtm+HvCoiCxyFs1gVL68/bad9a9x5vpGgeDIbljzLtS/RSekUcGnWClo/BBs/wm+qQEpgXdtiCpcUVEwfDj8+6+djK5WrdO32bgRfvvN76Ep5Xe5jWGvBAzCXmT6ETBORL7AVl/J9wwXxpip2KTcc914j/sGuNXbfZ31e4DTshNjzBPAE/kMWSmMsWXKBg+GypXdjiYHm74Ckwb1rnc7EqXy5oy7oHhN+OVSWDcBmj7idkTKBWXKwIUX2vuffQbvvw+NGtmJmqZOhZYt4e+/7fOvvQbVq0O/fqE77l0VTTn2sBtj9htjXjfGtAbaYpP2/sBs4BfssJi4wg5SqUDyxhswYgR8/rnbkeRiw6dQqiGUbuF2JErljQjUuAQqdoV1H+qkSopdu2yi/s47tmf9f/+zw2cAtm+HRx+Fiy6CO+/Ut4sKLb7UYV9ojLkVqAJcja1lDvCOM7zkEacCi1Ih6fhxm6jffjv07Qu33eZ2RDk4tAV2zoGaA7WbSQW/OkMhaRX8cpm9kFoVWbfdBmvX2qoye/fCI49Aq1b2uUqVYPdum6y/+qq91fruKlT4XIfdqV0+EZgoIrWww2WGAI8Do/JyTKUC3eHD8PjjjZkzB268EV5/3V5sGrA2fgEYqDnA7UiUyr/ag+HwNvjnYQiPhTMn6AdRleVbICICXnzR3n/5Zdi3z1bwUirY5Su5NsasBx4TkZFAxoRKSoWcbdtgx45onn8e7r47CHKF9ROhTCs7JEapYCdh0ORBOJ4KS0ZCsZLQ5lUtU6qyJGKT9iZNTl6ounAhPPUU9OwJQ4ZATIGV0FDKPwqkarSxfjTGXF4Qx1Mq0NSpA2PHLuSee4IgWd/5M+z9C+pc43YkShWspo9Ao/th1VibuCuVDRG4/nqboANs3Qp//gk332wnZRo7Fo4dczdGpXwRaNO8KBVQNmyAa66BAweCIFEHSDsEix+BqApQ9zq3o1GqYEkYtHoG6l4Py56EZU/Dsf16daHK1fnn27/nP/5oJ7i79VZ44AG3o1LKe5qwK5WNzZuhd2/4+mtbmSDgpaXA9DNh51xo8SRExLodkVKFo+3rUOMKO6nSl2XgsxiYVBVm9oTktW5HpwKUCPTpA8uW2QICLVuefO7ff211mapV7TeqjzgVRI3Rz4MqMGjCrlQW3n/fjn/cutXW+a1Xz+2IvLDwHti/BM6eDPVucDsapQpPeBR0+gS6z4CWz0DD4VC5D+yZBwvucjs6FQSeftqOZQe4/35b133uXDuE5qyzoEYN+9yKFVCxou2NT052L16ltKKLUpk88YSt5ZuQAG++CQ0auB2RFzZPhtVvQqP7oNoFbkejVOETsTP4es7iW7IB/PMgTGlsv2WqfpF78amgcPw4REbaC1KvucaWhvRUrJj9X/DsszBlil1q13YlVFXEacKuVCYDBtivQB96KMBLN2Y4sgvmXQ9lWkLz/7kdjVLuOeNOSD8Mm7+Fny+x5SArdIboeIiuCDGVILZGkFyQovwhLMx20mSnfn344guYORMuvRQGDoRffw2S/w0qpOiQGKUcS5bY3pZ69WwPe9D8QV4wHFL3w5kf2aECShVV4dHQfDT0/h3OuAs2fgZ/3gBzL4TpHeHbWvDLpXDsgNuRqiDTo4etLDNvnq02k+HoUZ2cSfmHJuxKAdOnQ7t2dlxjUNn0DWz4FJo8CqWbuh2NUoEhIgZavwCX7od+G6D3POg6BZo+Znvfp7WDxSNh32K3I1VBZMAA+L//g44d7eOBA6FECWjYEC6/HLp2hccfP7m9jnlXBUkTdlXkzZgB/frBGWfATTe5HY0PDq6Ev26C0i2gyQi3o1Eq8IRHQfEaUL49VD3P9r73mA3mOCz9H/zYClaNcztKFSRE4KqrTo6oio+H4cOhenX7DW16OsTF2ed27LAVZzwTeKXyQ8ewqyJt1iy44AI7TnHGDChXzu2IvJR6EGb1sIlHp4kQVsztiJQKDhW7wIWr4dg++G0QzB8OhEGdoTqkTPnk5Zdzfv7cc2HkSFtG8rbboEsXv4SlQpT2sKsiKzkZrrgC6ta1FxSVL+92RD5Y8QIc2gxdJ0NcY7ejUSr4RJaBsz6GMi3sN1W/X60Ft1WBiY+Hjz+2ifqMGdCrl70F+zZbscIm/E89dQaXXWYncjp40NWQVYDThF0VWSVK2EmRfvwRKlRwOxofHFgBK56HGpdB+Y5uR6NU8IosDX3+stWVNn4Bq99yOyIVQsLD4bXXYOVKW8xgrTOn13vvQePGcNdd8PffZViyBD79FCKcMQ+HD7sXswpcmrCrIiU1FcaNg2eesY+7dIFq1dyNySfpR+Dni6BYCWj9otvRKBX8RKDxg3bipfm3wtZpbkekQky5cvDHH3D99fZxs2Z2jo916+CLL37n339hzx6IjYUDB6B5c3jsMUhLczduFVg0YVdFxvz50KoV3HKLHbselN9+//cKHPzPlnCMDaZPGkoFsLBw6Pw5xDW1ZR+XPwNbf4SkNXaeg6D8Y6ECSYkStuY7QPv2MGwY1Kp1+nbh4bYj6X//s1Vn1q8/+Zwxtid+2DBbrUbflkWLJuwq5C1fbqu/dOwI+/fDpEl2GExQzZ1ijsOmr21li6oXQuXebkekVGgpVgoSvoeYyrBoBCSeC9/Vg68rwqxesHeB/YZLqUJUooQdMjNxIixdCi1bwgcf2Oc2b7ZVaiZOhKuvhiuvhGnTbHUaFfo0YVchTwS++gquu87+AezfP8iS9QPL4bsGdubGkg2h3RtuR6RUaIqtCuf/C5fuhe4zoOMEaDYadv0MP7aF75vZce6bv4XDO9yOVoWwgQNh0SJo1Mgm6GDLR06fbjueRo60/9cuvPDk/7PUVLeiVf6gZR1VyPr1VzjrLPsHb9u2kxf0BJW9f8Ps3iARdhhMjcu09JxShUnCbAWZSj1Orqs5EPbMg38ehj9vtOsiy9iJmGpeYXvllSpgtWvb/2N79pxc18N5W44aBffcYzuhwsJsL3uNGtC3r93vwAHYu9dO6NSnDxw5AnPn2m+aS5Vy5eWofArGFEapHKWkwJgx8MQTdrzfFVcEabKesgFmdodiJaH7TChV3+2IlCqaStW3S80r4PA2SNkIC++GhXfBfy/DuX/bBF6pAhYWln0Vs5Il4cwz7f0jR2xi/sUXkJQEMTFQtuzJWVn/+ss+HxkJd9wBQ4faSjUZ1q2z20RF2frxkZGF+rJUHgRjGqNUloyxXxdef70d6zdoEFxyidtR5ZExMO8GMGl2ZsaSdd2OSCkVVszOnFq8BvSZZ4fKzOoJ39aCqApQLA4q94LmT9oLWZXyk+LF7Vj399+3Q2MyJ9xt2sBPP8GHH8Jzz9klIQG++86Om7/8cluYAaB7d9uD37ixrXCTnAzz5tkPDhUq2B76mJiTF9Eq/9CEXYWMQYPsWL/69WHOHDj7bLcjyoe178H2n6DdWE3WlQpEIlDxbOjyDWydAsf2w5HttsLMpq9tFafK50Cj+4LsohkVzESy7h2PjYWePe3y5JPw2Wd2DHzx4vb5cePsN9ELFthKamefDW+8Ye//+6/dL7NPPoEBA+yET8uW2XbXrLHDc44ds8NzKlSw5SnDw/XXIL80YVdByxjbO1CtGrRuDTfcYHsMBg2yn/6D1qEt9uv2il2h3o1uR6OUyknVvnbJsOZd2PQNHN4Mix4ABBrdq9mKChjVq8O999olQ9u29rZlS7joIvjtN3v9F0DDhpCYCLt22SUpyQ49bdrUPv/mm3D//ae3s3OnvX38cVuGsksXO3Tn2DGoWhVeecUm8gsXQlycnXXcGFi8GI4eteUv1UmasKugdOCA/Qpv+nQ7Vr11a5usJyS4HVk+GQN/3gTHU6HDO/YCOKVU8Kh7nV3Mcfjlclh0P2z5DkrWg6jyUKIO1LrKXpuiVAAqWxbOP//k45IlbU347AwbZj8EREbapDsy0o6Fzxh737mzrXgzc6ZdHxsLK1bYZB1sj//XX9vx9klJtrf+qqtskn/smO2Eu/BCe+yjR2HVKqhT5+QFuEWFJuwq6Pzzj/0F/u8/eP31k7PHhYT1E+3X661fsv/glVLBScKg02ew4lnY9BVsmwbH9tpa7osegDrXQcPhUKKW25EqlS9xcXZoTHZ697ZLdp580o6xnzwZypSBsWPhggvsc2vX2lliv/ji1H3mzrW3331nK+mkp9tx9vv22So5n3xih/ssWWI/INSpE/xfcmnCroLKd9/BxRfbC2GmTs16XF3QOrQVFgyH8mdCg9vdjkYplV9h4dDkQbtk2D3Pzli88jVY8w70mgtlWroWolJuO+MMeOghu2T13Nq1tnzl9u0QHQ2VK0ODBvb5xYvhpZfsl9Nt29pe+MOHT47Nf+AB+OGHkxVz2rSxF9NmfMAYP96WvxSBdu2gfHk4ejQwv9nWhF0FtOPH7S9jlSp2Iomzz7azlo4aZZP2kHE8DX4baHvfOr6vFSaUClXlO0D5idDiCZjRFRLPh/OWQWSc25EpFZAiIuzY+qw8/DCMGGEvbI3KYoqSp5+2nXx//AG//247+rp2PZmwv/aanQ3d0xVX1KJPnwJ9CQVCE3YVsFatsleoz5gB115rE/a4OPsLFnIWPwY759rJkUo1dDsapVRhK1EHunwN0zrA3/dA27EQrsWvlfJVePjJ8fCZtWhhl4yhs6mp9sLXDPPmQbFitld+wQI7pObgwR1AjUKP21easKuAk5Zmv6a691578cpbb4XYOPXMtv4Ay5+GujdA7UFuR6OU8pdy7aDhnfDfS7Duw5MXpYbHQngMhEc7tzF2YqbyHaFUI51ETak8KlbMLhlKlLC3UVEnL2JNTEzxf2Be0IRdBZzff4fbb4dzzrGTQFSq5HZEhShpNfx+NZRuAW1ecTsapZS/tX4eKveBnbPhyA5IXg+pB+HITkg/bIfJpR+GY/vsRGoAnT61s64qpYoMTdiVq/bts1d/jx0LvXrZ2dc6d4bZs+04s2C/qjtHKRtgZg9AoPMXEBHMxeOVUnkiYVClj11yknYI9v1jK8z8diWseAGiykJEcdsjHxFrb2MqQYPb7HqlVMjQhF254u237UQL+/fbx23a2FtjbJIe9PXUc3PgX5hznu1J6zlbv+JWSuUsIhYqnAldJ9tkffdvttf90BZIP2QT+rQUSEuCg/9Cy2chooQdVhPSPR9KFQ2asIe4DRtg0iQ7SQHYyiqDBtmyRoXt+HFbM337djvxwowZcM899v7evbYHvWPHk7dF5n/K7nkws7vtUe82TUu6KaW8F1kaWvwv++f/eRiWPQVrP7CPJRxiKkNMFVokpUFiFVuRptH9UKa5PyJWShUATdhD2IcfwpAh9n7Firb3OiUF7rzTrhs3DqZMsRMKNGliL7ioUyf7q61zkpYG//4LpUpBjRp2quGLL7YfGDKI2OEuPXvCffdBWGCWOi1cqcl2zHpUeej9O8RWcTsipVQoafY4lGkNh7faHvfUg3BkGxzaSpjZAoe3wK5fYMsU6DETyrZxO2KllBc0YQ8xR4/C5s128oC+feGZZ2ziXM+ZNNOYk9sePmy3/fVXOHDArqtUySbZkZF21rG9e6FZM6hf304RfOQIVKtmt336aUhOhq1b7YRGe/bYi0VffdW217QpPP64jWX3btuLHh9v9y2SyfrxVPh1ICSvge4zNVlXShW8sHCocUmWT/2dmEhCQgKkbIQZZ8OcC2wp2XLtIKJkEfqaU6ngowl7CFm9Gq64wo4LX77czth1//2nbuP59/juu+0Ctnf8119tD3ykUwr40UftLGKeevWC6dPt/ffeg/Xr7RCXc86xHxDatrXPlSple++V4+ge+PliW2u93TiIT3A7IqVUUVW8BnSdAjMTYFbGdNFysoxkWCSEFbO3EbE2mS9e05aUjDsD4nvYC16VUn6jCXuI+OoruOYaOyPYBx9kPeNXTs44wy6efv3V9p7/849NzCMjTx37vnKldsh4JXkd/HwJHFhue7O01rpSym2lm0K/jbD9Jzj4H6Ql2/KRaYfh+DEwqZB+zF7QmnoQdv8OGz4FjK0JX7mPR4Wa4nacfFQFCIuw4+YlAsq1t1VrlFL5pgl7CPj8czsLaPv28Nlndgx5QShRAho0sEtWNFnPRdphWDoaVjxve6vOngRVznU7KqWUsiJioVo/77dPOwT7FsHyZ2DvAqcqjbNk1Ij3FFMZ+szX4X9KFQBN2IOcMXYm0LPOgh9/hOJaejcwbP4OFgyHlPVQ51po/jjEVnU7KqWUyruIWKhwFnT99tT1xtiJnlL3w/E0MOn2otdfLoPEvnD2N3YIjRRzeuAjtMdHKR9pwh7kRODbbyE9XZP1gJCaDAvvgjXvQFxj6DFbx6srpUKbCMTE2yVDmebQ5WuYeyFMrp3FPhFQ8Ww7ll4njVMqV5qwB6nly+Hhh+Gdd2xtdeWy5HXw78uw7gM73rPJQ9BslB0Ko5RSRVHlXtD7D3ux/fFUOy7+uLMc2w8rX4OfOkOJWk6vuzP2PcyL+xJhK+Jk3Jdwj977HO4XxvH12wLlB5qwB6Fdu+DCCyEpyZZmVC7avwz+e9lOUiICNS6HBrdD+Q5uR6aUUu4r08IuWYk7A1a9aS96NenOcBpnSI1JOzm8JvN9k+7f15AbCcvbBwIJp1VSCkwv7ST9YR63YYDYb2pbPAXFSrj7GpXrNGEPMocP22R982aYPftkTXRV+MSkwp75sOdP2DMPdv1qa6qHRUH9m6HxAzpOXSmlvFX/Zrv4yhgwx53kPe30ZD+7+zl9CMj2w4E3x/fY9rgv+6WRnmLstQHmOJDxujziWPUGrHnXltw8Lan34jYj8c/zbQ7HDouERvdBxS75fScoL2jCHqBqH3wbfn3Tls8Ki4JiJTlerCJX33U+8+bV5PPXf+XMqttgd00oUceOAQyPcX7JVJ4YA6kH4PA2OLLdud0FKetgz5902T0ftqXabaMrQrmO0PBOqHkFRFdwNXSllCoyRGxPNeGAjzWMA8zijMmssrNjDmz6GjjuzHyYy61n4u/Ttjndpme9fcpGe1Fx7SGZEvzMCb/HY2/XZTz29Vg57ZfbsYqVgMrn5Od0FipN2APRkZ3UTJ4IqRWcT+KpkJbMtr2V+XP+hTw38F4uLf0i/JLFvmGRNnHPWCJindsSNvmPLG1vizm3GeuKxUFknF0fVT60LgIyxo4rP7QZDm2ElA12ObzNVjY4utPeHtkJx4+evn94DJRtw5biF1G91aW2tnBsDR23qJRSqnDFd7VLIDq0Beb2gw0TnQ8CznIisc9hXcbjQFOiHs2PxcH+96F0M7ejOYUm7IFoz1/2tsuXUPFs+8E2PZ2qaXtYfNER4kreAgyzk1scXGkTz/RDdtKLjIkvMj9OS4KkVXBsny29lZaSffsSDsVrQXiUneEutpq9eFLCs1/CoyCynL0fVc4uEsFpX6md8snY86s37OtJP+KxHLa3x7NYl7EcP+JM9HHk1PUct+2nJsGRHacn4hJhawRHx9uldDM76Ud0vLO+knNb0X6ICQtnTWIi1WskFOCJVkoppYJUbFU4Z37+j5M5ifc5+ff89iAf+x1YCqvGE350h7NtYHE1YReRc4BXsN9rvWOMGZPpeXGe7wscAoYaYxbmtK+IlAU+A2oB64HLjTH7nOceBK4D0oHhxphphfwS82bvfAyClGlNWhrceSekpYUzblxFSmceIp3XT4Dpx+zwj2P77NX6qQecx/tt73PSKjuGLvWAnaHzxBhAz+X4yfsZiXRhkjDnm4OM6bOjnfvOElHCfjsQFmW3NWkQXtzOtBcdDzFV7PTaxWvahDwsvHDjVUoppVTOThme4qK4RlDjMv5OTCQhuwulXeRawi4i4cAbQC9gM/CXiEw2xiz32OxcoL6zdADGAR1y2XcEMNMYM0ZERjiPHxCRxsAAoAlQBZghIg2MCbTLzYE9f5ESXoPFf5Xg7rth3jy4554CbiM8EsIrFNzYa2Ns0m7S4egeOLbXSei9HSN33PbSZyTkpyTjGev0CyGllFJKFT1uZkDtgdXGmLUAIvIp0A/wTNj7AR8aYwzwh4iUFpHK2N7z7PbtByQ4+08AEoEHnPWfGmOOAutEZLUTw++F+Bp9Zwy//Wa49IW/2LYXKlaEiRNh4EC3A8uFiB0vD1CsJPYUKaWUUkqp/HIzYa8KbPJ4vBnbi57bNlVz2TfeGLMNwBizTUQqehzrjyyOdQoRGQYMA4iPjycxMdH7V1QAwo4fJiquJC2bbuTKdnvp3n0nsbHp+DkMlYXk5GS/vx9U7vS8BB49J4FJz0vg0XMSeAL1nLiZsGdVYiPzJcPZbePNvnlpD2PMW8BbAG3btjU5llsqLN3P5f4miSQktAEa+r99laXE3MpvKVfoeQk8ek4Ck56XwKPnJPAE6jlxc4T/ZqC6x+NqwFYvt8lp3x3OsBmc250+tKeUUkoppVRAcTNh/wuoLyK1RSQSe0Ho5EzbTAYGi9UROOAMd8lp38nAEOf+EOBbj/UDRCRKRGpjL2T9s7BenFJKKaWUUgXBtSExxpg0EbkNmIYtzfieMWaZiNzkPD8emIot6bgaW9bxmpz2dQ49BvhcRK4DNgKXOfssE5HPsRempgG3BmSFGKWUUkoppTy4WifPGDMVm5R7rhvvcd8At3q7r7N+D9Ajm32eBJ7MR8hKKaWUUkr5lctV6pVSSimllFI50YRdKaWUUkqpAKYJu1JKKaWUUgFME3allFJKKaUCmCbsSimllFJKBTBN2JVSSimllApgmrArpZRSSikVwDRhV0oppZRSKoBpwq6UUkoppVQAEzuZqMqKiOwCNrjUfHlgt0ttq6zpOQlMel4Cj56TwKTnJfDoOQk8bp6TmsaYClk9oQl7gBKR+caYtm7HoU7ScxKY9LwEHj0ngUnPS+DRcxJ4AvWc6JAYpZRSSimlApgm7EoppZRSSgUwTdgD11tuB6BOo+ckMOl5CTx6TgKTnpfAo+ck8ATkOdEx7EoppZRSSgUw7WFXSimllFIqgGnCHmBE5BwR+U9EVovICLfjKcpEZL2ILBGRRSIy31lXVkR+EpFVzm0Zt+MMZSLynojsFJGlHuuyPQci8qDzu/OfiPRxJ+rQl815GSUiW5zfl0Ui0tfjOT0vhUxEqovIbBFZISLLROQOZ73+vrgkh3OivysuEpFoEflTRP5xzstoZ31A/67okJgAIiLhwEqgF7AZ+AsYaIxZ7mpgRZSIrAfaGmN2e6x7FthrjBnjfKAqY4x5wK0YQ52InA0kAx8aY5o667I8ByLSGPgEaA9UAWYADYwx6S6FH7KyOS+jgGRjzPOZttXz4gciUhmobIxZKCIlgQVAf2Ao+vviihzOyeXo74prRESA4saY/2/vbkMtq+o4jn9/jQ+IRZEPo6ig5CSFhFlJIeW8qMGHaKgUrZRJIiUc0AjBLDXL0l5YZg9KQ9bMVIqmo9ILMwg1QXQyzBpnIh/KhhmdUmvG1Hzo34u973Q63HPvnTszd+/r/X7gcM5ee521FyzWuf+7zv+s/WySXYG7gbOBj9DjueIKe78cBTxcVY9W1YvAdcDijvuk/7cYWN6+Xk7z4audpKruAp4eKh41BouB66rq31X1GPAwzZzSDjZiXEZxXGZAVW2sqt+2r7cAa4EDcL50ZoIxGcUxmQHVeLY93LV9FD2fKwbs/XIA8NeB4/VMPLm1cxVwe5L7k5zRls2vqo3QfBgD+3bWu7lr1Bg4f7q3NMmDbcrM2NfJjssMS3Iw8HbgXpwvvTA0JuBc6VSSeUkeADYBv6yq3s8VA/Z+yThl5ix15+iqOhI4DjirTQNQfzl/unUV8CbgCGAjcHlb7rjMoCSvBW4EzqmqzRNVHafMcdkJxhkT50rHquqVqjoCOBA4KsnhE1TvxbgYsPfLeuCggeMDgQ0d9WXOq6oN7fMmYBXNV2BPtnmJY/mJm7rr4Zw1agycPx2qqifbP4L/AZbxv6+MHZcZ0ubj3gj8pKpuaoudLx0ab0ycK/1RVf8A7gCOpedzxYC9X1YDC5IckmQ34BTg1o77NCcl2bP9kRBJ9gQWAX+gGY8lbbUlwC3d9HBOGzUGtwKnJNk9ySHAAuC+Dvo3J439oWt9mGa+gOMyI9of0v0AWFtV3xg45XzpyKgxca50K8k+Sd7Qvt4DeD+wjp7PlV1m+oIarapeTrIU+AUwD7imqtZ03K25aj6wqvm8ZRfgp1V1W5LVwPVJPgU8DpzUYR9f9ZJcCywE9k6yHrgIuIxxxqCq1iS5HngIeBk4y90Vdo4R47IwyRE0XxX/GTgTHJcZdDRwGvD7NjcX4HycL10aNSYfc650an9gebsz32uA66vq50nuocdzxW0dJUmSpB4zJUaSJEnqMQN2SZIkqccM2CVJkqQeM2CXJEmSesyAXZIkSeoxA3ZJepVL8qMkvd4SLMm+STYn+fRQ+d5JViTZkKSS3LENbV6R5I/tzWskadYyYJekWaYNXKf6OLjr/k7RJTR3FvzhUPnlwMnA1TR7Wn91G9q8jOYOhZ/ZER2UpK64D7skzTJJTh0qei9wBvB94NdD51YBLwLzquqFGejeNktyIM0NZD5XVd8aOrcB+E1VfWiabV8DHAccVFUvb29fJakL3ulUkmaZqvrx4HGSXWgC9nuGzw14aad3bPrOpLnr47XjnNsPeHo72l4JnA4sBm7cjnYkqTOmxEjSq9x4OexjZUn2al//PcmWJDcn2a+tc0aStUleSLIuyeIR7Z+c5O72/c8luTfJidvQxZNoVtE3DbT5pbbPAZYMpPh8sj1/QpI7234/n+TxJDclefNQ23cB/2qvIUmzkgG7JM1ttwGvBy4ElgEfBFYlORc4F1gOnAfsBvwsySGDb05yCXAdsAW4oK37HHBDkrMmu3iS+cBhwH1Dp26iyVmHJs3ntPZxV5JjgFvbfl8KLG37vhdw6GAjVfUKsBo4ZrK+SFJfmRIjSXPbfVW1NbBOAvBZ4ADg8Kra3Jb/CvgdTerN59uyI4EvAJdW1fkDbV6Z5Gbg0iQrqmrLBNd/a/v8yGBhVT0IPJhkJfDoYKpPkqU0C06LBlflga+MuMYjwMIke1XVUxP0RZJ6yRV2SZrbrhg6HvvR6oqxYB22BtCbgQUDdT9Bk3u+vN1+ceuDZgX8dcB7Jrn+Pu3ztuSp/7N9/mibvz+ZsSB93224hiT1hivskjS3PTp0/Ez7/Ng4dZ+hSTsZ8xaaHPN1E7Q/f5Lrj+XWZ5J6g75D8yPS7wFfT3I3TWrPtVX1t3Hqj7XttmiSZiUDdkmaw9oc7/GMKs/Q66LZNnFU/TWTdGEswH7jJPW2qqqnkryLZjvLDwDvA74JXJzk+Kq6Z+gtY22PF8xLUu8ZsEuSputPwLHA41W1dpptrKEJ+hdMVnFQ+4/GHe2DJG8D7ge+CJwwVP1Q4Anz1yXNVuawS5Kma2X7/LUk84ZPJpk0Z7xNYXkIePdUL9rmyA9bBzzP0Ep92693AndOtX1J6htX2CVJ01JVq5NcBFwMPJDkBmADsD/wDuB4mu0gJ3MDcEGS/atq4xTqL2vvjno78BdgD+Bkmh+5rhiquxDYs72GJM1KrrBLkqatqr5Ms3f7BuAc4Ls0Wz/uDpw9xWaW0aTFfHyK9VcCG4ElwLdp0mBeAk6sqquG6p4KPAHcMsW2Jal3UuWP5iVJ3UpyNbAIOKyqXtpBbe5HswvOeVV15Y5oU5K64Aq7JKkPLqTZMvL0HdjmecB6YHjVXZJmFVfYJUmSpB5zhV2SJEnqMQN2SZIkqccM2CVJkqQeM2CXJEmSesyAXZIkSeoxA3ZJkiSpxwzYJUmSpB4zYJckSZJ67L/fiBWqjGW6SAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of files processed in Dataset 1: 500\n", + "Number of files processed in Dataset 2: 1000\n", + "Time range Dataset 1: 0.0 to 300.0\n", + "Time range Dataset 2: 0.0 to 300.0\n", + "Average entropy at start (Dataset 1): 0.0\n", + "Average entropy at end (Dataset 1): 1.3913508375000012e-05\n", + "Average entropy at start (Dataset 2): 0.0\n", + "Average entropy at end (Dataset 2): 7.821582446250007e-05\n" + ] + } + ], + "source": [ + "\n", + "\n", + "def process_data(directory):\n", + " columns = [\"t\", \"Mean-Field Energy\", \"PE Difference\", \"Total Energy\", \"norm\", \"last_column\", \"total count\"]\n", + " file_paths = glob.glob(os.path.join(directory, 'ene*.dat'))\n", + " dfs = []\n", + "\n", + " for file_path in file_paths:\n", + " df = pd.read_csv(file_path, delim_whitespace=True, names=columns, skiprows=1)\n", + " dfs.append(df)\n", + "\n", + " combined_df = pd.concat(dfs)\n", + " average_df = combined_df.groupby('t').mean().reset_index()\n", + " return average_df, len(file_paths)\n", + "\n", + "\n", + "\n", + "# Process data from both directories\n", + "avg_df1, file_count1 = process_data(data_directory1)\n", + "avg_df2, file_count2 = process_data(data_directory2)\n", + "\n", + "# Plot the data\n", + "plt.figure(figsize=(12, 6))\n", + "#set legend font size to 18\n", + "\n", + "plt.plot(avg_df1['t'], (9/8)*avg_df1['last_column'], linestyle='-', label='pointer', color='orange')\n", + "plt.plot(avg_df2['t'], (9/8)*avg_df2['last_column'], linestyle='--', label='adiabatic', color='blue')\n", + "plt.xlabel('Time (fs)', fontsize=18)\n", + "plt.ylabel('Average Entropy (nat)', fontsize=18)\n", + "plt.title('Comparison of Average Entropy vs Time', fontsize=18)\n", + "plt.legend(fontsize=18)\n", + "plt.grid(True)\n", + "plt.savefig('entropy_comparison.pdf', format='pdf')\n", + "plt.show()\n", + "# Print some statistics\n", + "print(f\"Number of files processed in Dataset 1: {file_count1}\")\n", + "print(f\"Number of files processed in Dataset 2: {file_count2}\")\n", + "print(f\"Time range Dataset 1: {avg_df1['t'].min()} to {avg_df1['t'].max()}\")\n", + "print(f\"Time range Dataset 2: {avg_df2['t'].min()} to {avg_df2['t'].max()}\")\n", + "print(f\"Average entropy at start (Dataset 1): {(9/8)*avg_df1['last_column'].iloc[0]}\")\n", + "print(f\"Average entropy at end (Dataset 1): {(9/8)*avg_df1['last_column'].iloc[-1]}\")\n", + "print(f\"Average entropy at start (Dataset 2): {(9/8)*avg_df2['last_column'].iloc[0]}\")\n", + "print(f\"Average entropy at end (Dataset 2): {(9/8)*avg_df2['last_column'].iloc[-1]}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 84, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000\n", + "749\n", + " t S0 pop S1 pop S2 pop S3 pop S4 pop \\\n", + "0 0.01 1.0 2.604167e-09 2.604167e-09 2.604167e-09 2.604167e-09 \n", + "1 0.02 1.0 1.041665e-08 1.041665e-08 1.041664e-08 1.041664e-08 \n", + "2 0.03 1.0 2.343738e-08 2.343737e-08 2.343736e-08 2.343735e-08 \n", + "3 0.04 1.0 4.166625e-08 4.166622e-08 4.166619e-08 4.166615e-08 \n", + "4 0.05 1.0 6.510313e-08 6.510305e-08 6.510297e-08 6.510288e-08 \n", + "\n", + " S5 pop S6 pop S7 pop S8 pop poptot \n", + "0 2.604165e-09 2.604165e-09 2.604164e-09 2.604164e-09 NaN \n", + "1 1.041662e-08 1.041661e-08 1.041661e-08 1.041661e-08 NaN \n", + "2 2.343723e-08 2.343722e-08 2.343720e-08 2.343719e-08 NaN \n", + "3 4.166579e-08 4.166575e-08 4.166570e-08 4.166565e-08 NaN \n", + "4 6.510201e-08 6.510190e-08 6.510179e-08 6.510167e-08 NaN \n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/intel/oneapi/intelpython/latest/envs/2022.0.2/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1402: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead.\n", + " ndim = x[:, None].ndim\n", + "/opt/intel/oneapi/intelpython/latest/envs/2022.0.2/lib/python3.9/site-packages/matplotlib/axes/_base.py:276: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead.\n", + " x = x[:, np.newaxis]\n", + "/opt/intel/oneapi/intelpython/latest/envs/2022.0.2/lib/python3.9/site-packages/matplotlib/axes/_base.py:278: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead.\n", + " y = y[:, np.newaxis]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAGaCAYAAABdWLbHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdB3xT5foH8F/SPWjLKHvvjSAgIigiICKoqAgiijgQFfR6VS6oV0AU99YrKm4QcCICykb2kL0F2XtTWujO/3mOp/zTJm2TkrRJ+vvez3NP8p6RtznSPnnzDovNZgMREREREfkGa1FXgIiIiIiI/h8TdCIiIiIiH8IEnYiIiIjIhzBBJyIiIiLyIUzQiYiIiIh8CBN0IiIiIiIfwgSdiIgCnsViGSlhk/iyqOtCRJQfJuhE5Lck2QqWuFfid4nDEqkSpyW2SkyX+I9EqzzOVwMllkmckTgnsVbiaYnQwvxZyKXkuiDRge8tEfmb4KKuABFRQUjiFS+bGRIt7YqTdZdEPYn6Et0kzkrEOTk/RDZTzGNUqkSGxGVm9JJjOtpstkTeoSKn9+Cok3L9EFXSfHzCvH85pdrt3y5x2OO1IyLyMAtXEiUifyTJ80zZdJE4JzFa4hv5fXbE3FdCNldI9JS4UcqrOzn/VdkMNZP6QRLjJTL1eImvJEpJfCvn3uX9n4YKwmwdn28+rSH3ag/fSSIKBEzQicgfEzNtHd9qPu0lidkPeRwbIfsv5CgrLxtN5sIkHpf97+XYf7PZum6TuEz2b/Bk/ckzmKATUaBiH3Qi8kdN7B5Py+vAnMm56TYzOdfuL584OecX2fxldpfp607FdBCi2fdZ+02HS4yS2CZxQeKYxESJuvlco5zEm+Z55yXOSqyUeFIizFuv6+SaERIJ5nW753PsNvO4x3KUXyPxg8QBc4yA/iw7JKZIPCRhLepBonb91atL1JOYYI5pOG+OSbjbybiFP80xC6ckJklUzef19drvS2w3r6vnrjbHSUR542cmIv/FBJ2I/F2lApxzrbldKMm4dnFxZpa57ViA66sws/vF8xI1zL7Q2m++j4QmfVc7O0nKW8tmi8S/zb706WZfax3s+obECjmmrKdfN48PN/pNAvL6oCLXbGHWVfuAT7YrHyibBeYHIr1PaRJBErUl9FuKsebP5iv0vV9l/qwREuHmeISvzQ9H+oFtgsTHEk0lLGYf+N4Si2R3aWcXlfJbzW98BktkfUjSn1vft1ckdJByOW/9UETkf5igE5E/Wm33+ENzwKg7GprbzXkco0myamAmZu562Ezi+ktES7IbK9vmEmskIiW+k8tmDXA0mM+nmP3fN0q0lvNi9HyJXhKnJZqZSaLHXjcf35rbm+Q8Pd+ZO83tPHk9YzCneeybZvnnElVlX5SE/iyayN4gMdHs9+8r9NuUPyRqSj11YHGc+SFCvWBGDwltUdefQ8c6tJfQsQ/agv6fnBeU90E/WE2S0EHJOu6hmoS2mOv700ZihfmN0Nfe+qGIyP8wQScivyPJ0y67hOZ6Ce0+MUfiRe0/7kLCXsHcHsrjmKx9mohpuEsT44FS168l0sx6rzPre1JCW0wfzXHOYLNuZyS6yPGrzPMyzH722gquOukMMx583bzMljhmJpXa6p2N+eFFW5Dtk3nV2Hzfksz67M/aIY9PSfwu0Vcia5YVX6A/Z0+p0259ItsE873aaSbUz+lzKR+v9ZZQi83Bxup2J9d820zOn5Zjh0nsM8/Te7rC/KCi/611kbfSfkYiIirGmKATkb96UOItiVSzu8B1Es+aLdDHzD7bd+XS+p3V59dZ//Qs5+0eFyRB35sjYTVIUnbC7CLhLKHLej4ua0aaHOdqt5tl5tM7PPi6udJEUjbf52gpt9dOooqEdhX6ya5ck1uYyanTrh8+6A35ebVL0UXyXFv455lPD5iz/eQ019zWsO9PLo9ryeYq87+zrJZ45Li+fivym/m0c8GrTkSBhAk6EfklswXzSTM5HGR2l9ihu8xDWpnJ1OQ8BiJmHesNf2gzaW77zG1jqZvRB9vcaqsz7KYOdCYrWWzhidd1UVbC31XO0+439rL6pk83W5yz6L3QCDX7WD+hs+8UsLtQYdFuRbm1rKstZsKek/0c7fZz7rc1t/oe7JYf/YizsPtmRP9bJiJigk5E/k0SpmMSH5vdJeqaXUS0dT2rS4X23R6S4zTtdqFy61Odc19BFis66MI+HTCZ1R9cE1+rC+dqK66K99Dr5kve16Wy2W22huuAT4Mkl8F2rfHfOml572u+Zk3z2w4dKHlCzvte4iYfTNZzW8QoI6/95s+aRd+jnF2pgsyuRblFlAv/PRJRMcIWdCIKKDpIUWKc2cKc1bJ5Xy79yyvmcamKdsm5p1cTzS8xdTqVYiG8bl50oGPO2Vy0S0YZc7rK6TlPkPvwp2zqSPQzxwzsMj+IaFKvU1lOlxxdk9dA/xu7Vt4LXXckv7i3SGtLRD6DCToRBSSzz7UmgapuLjO0NHJhppeteXQZyUteyX9Wy6q2vGofZHVKIqv7hM70kZvK5va4h17XVVkzx1wtSXXFHH3Sf5K3KCW3qRolJkj0l6hltqa/rLvMAZLaPSlQZX1ArGN+20BE5BIm6EQUyLK6suScKSSrj3d7XdQnl3M75xgA6K5rXNi3KWsWE3O7Kcc87c5kzd6yxhOv6yo5frPZR1v/bvQx37dbzN3funGd3RLP2M2Xnld9/d0yu0HGXYqyIkTkX5igE5HfkeRQZ8uolc8xkXYJpE4zaE9nG0kxB/Q94OTcHubCOzZz8GlB6MqRDrOemIMsdQEf2M2OkkWnUlT3ynEVnJyrSd6V5tPvPPi6rspKxPX6Pcx5wI84G9TqwiDUC17uzlPk5IPINtksN5++mteKoeaqrQH7XhCRe5igE5E/0q4pumT6TxJ32CezmgSZCfYicyVN9a79yeYUhlllr+lS7ll9oWXbTTZfmPsmyrEbClhH7Zf9qVyvX1b3BtnqAkIzzQGeOjPI/3Kc84E5EFFXsfw9a15srZvEbXb9wOdIveZ58HXdSdD1Q4vWa7hZNjnHIMks3eR1dfaWByWq2X9w0jJ5eJdZpPUKZEPMD4ONzdVGO9ndF6tEIwmdX/1vCYcPZURUPLFPHBH5o6wl43uaocmOtsimmgv1ZNHE8XlJIO3n587ynJk0dTMHMGpSm2E3k8aqS+wf/ZFEB4lvJMbJtTVJ01VBs+ZY72XOgX2RPpfjtNX/dwlNqlfJ83PmzCBZXXE22CW3HnldV+kiO3K9pebc3ro6aX7dW3SlzDZ29yfZ/NYia7DqDHP1zoClA2XlZ+9pfhPT3Fz4KdW8rzE5Zn3x5rSfRORH2IJORP6Y9Mw0u6A8ZS5MpCs9ZvX1PWP2z35HopkcOyaXa2iS38NMwrUbQoqZIK0zl2xvJ8doElVQKWZf8hfMxYNCzYGd2greQq69MJd6rTQHqOoKlH+ZCZwunqMzojwtcYVOLenp1y3AYFH1t1lfZ7SF/26Jr8y+6+fNLjG6mukcif4SPXIuDBSI5Gf8zRyo/KL532bWBxWdN14/8Dwv0UCO0/tFRASd1olvAxGRh0jL6Jdm8jlKfr+OLKw3tqhel4iIPI8t6EREREREPoQJOhERERGRD2GCTkRERETkQ5igExERERH5EA4SJSIiIiLyIcVmHvS4uDhb7dq1i7oaVEBJSUmIisp1ET7yYbx3/o33z7/x/vkv3jv/tnr16hM2m00XhyuQYpOglytXDn/+qdMIkz9asGABOnTQtVfI3/De+TfeP//G++e/eO/8m8ViuaR1DdgHnYiIiIjIhzBBJyIiIiLyIUzQiYiIiIh8CBN0IiIiIiIfwgSdiIiIiMiHMEEnIiIiIvIhTNCJiIiIiHxIsZkHnYiIiKgoJScn4/jx48Y2PT09z2NjY2OxdevWQqoZ5Sc4OBjh4eGIj483tt7GBJ2IiIjIy86ePYujR48aCV758uWNhM9iseR6/Llz51CiRAneFx9gs9mMD1SJiYnYt2+fsfilfoDyJiboRERERF524sQJVK5cGZGRkXyv/YxFPkiFhISgZMmSCAsLw5EjR7yeoFt98E34XOKYxKZc9qv3JHZKbJBoUdh1JCIiInJHamoqIiIi+Kb5uQi5hykpKV5/HZ9L0MWXEl3z2H+DRB0zBkp8VBiVIiIiIrrUlljyb5ZCuodWH+zns1A2p/I45GaJr+U4tVwex8mbVSG/615IOocq5UoiOjLceHOdRdeuXZGUlOSpH4WIiIiIqFj0Qa8ksd/u+QGz7HBeJ5XNPIo/H46WR2FIs0XiUEYs9qVEY29SKHaftWDXyVT8dWgrqlWKx8mzFy6ep33Ftm/fbvQbu1Rz5sxB586dL/k6Of3444+45ZZbYLX63OctIiIiInKTRUem+hppza4um2lSt8ZO9k2Xzcuyb7H5fK5shsrz1U6O1S4wGqhUtuTlLz3aEyG2FITbLiDOdhZlpKG+Ak4g3nI223lnNYFHPI5ZyiAlKAbWoFCkWcORYglHKkKRYQlGpnz5EGTLkP/PQBDSEWyTkG2ILc3YBss2RI7W5yGQsKVefG6Vsy1G2GA13n+b8Tgr9P/1+kZYrHqU8Sr6uvoqGfoKFnlubPX5/z829htlWY+N2pnnZh2nNZGfyRJq7JM3Ktd7oYMiGjZsaIw2L0o6cjo6Wj9gkb/hvfNvvH/+jffPd+igwtq1a7t8fEaG/HUPCvJijaigdu7caczKk5drr712teSmLYtTC7q2mFexe65N24ecHShvzCey0UCVKlVs94743OkFIyLCUKdyPOpUikOdsuGoVcqKmtGpqBG+H1VsJxCckVngyiZlhOC8TSPUiASbtOBLgmyTyJTk2QhJnDVd13Rfk3N99k/aLVtLhpH8a2IvV0C4JPuhljT9HgBhlnR5nC7bjALXL11+tMQ0C85ppNiQmGrDueQMJKZk4pw+TgXWTrOZ+2CUnb5gw4nzNpzM2kqkZXq3hX/BggXo0KFDgc+nosN75994//wb75/v0DnN3Zk2kdMs+i6dB7158+ZefQ1/TNCnSgyW1vFJsr1C4qwk4nl2b8mvU/+FCynYsOOAEdkEhSC6XhuULl0KJaPDUCpKnodZEW61yadai5Fop2ZYJYBkaR/PKFUDqSExuJAOnE+z4YJEZno6bOmpEinyWLcSaSnG838e/7PPlpGmtTRasy3Sam60apthsYbAEiwRFGzUyRIcKo+ztsGwyr7wYAvCpE7hIRJyWJiGfPAOD7L981ySeCOpl4iwpCBSahwtEWlJRlSQPA66gMiIf8pipKyCfMsQBQndLxEhHwTykpBmxclkiyTtwImkTIl0nEhMN5L4WWP64JtnbTiSaMPRxEwcTbLJ+5P7tXQaoy1bthjzxBIREREVNz6XoEsiPVE22lRaRh5rxjxCIkT3SSI+VjYzJLpJ7JQ4LzHA1Wlx8urO47R/uCTNiVsWIVEe7nX7J/Ex1qB/knq7xN4SEg5LaDisoRGyjTC21tC4i4/tt6GhoYgOD0GJMAtKhmWiVGg6SoemoZQ1CSVxDiUtiSgZKdvIRFQocw4NjbJzKCHJvTPnJKE/dsFqJO1HzqXjwJlU7D5tw+4zmdh1+ixqV62AJLsk/o033tCviy4+1243S5YsQVxcnLffOSIiIqLinaBLEn1nPvs1y37U06/bqVMnhwR+165dqFOnDjIzC97FJafrr7/e6PIRFRXlket999136N27d/4HZmbAppHmPGEuKCOJDy8hEY2giH+21vByErWM8nD5OeNjwlG+RBDKyY9cNjwDZa3n5NPXWZSJ+iealD+DLpaTDsn88eRgbD1txcZDKaifugqdagbhr5OZ2H/WZrSwa0u7M9oV5tdff2WfdSIioiLw5ZdfYsCAAZg7dy7WrFmDjz76CAcOHEC1atXw7LPPon///g6NpK+99hpWrlyJ5ORk1K1bF4888ggGDRp08Zh7770XkyZNwpkzZ4wuJmrp0qW46qqrjHxAF4LK6kr722+/oVu3bpg8eTLuuOOOwvvBAzlB9yU1a9Y0Bmn4Mv0Pz1v/8ek/qssvvzzPY2ypF5ChkXBMOs84d0Ri48VnFlglkQ+KipMoCWukbKPrGc/LlIhATcm5q8fYUCM6DdWCT6BuuQO4p/wBlEj5Ddfd/c+HmgvSrWhnQgi2n0jH5kMXsPlYOubvyTD6w2f1uXTWz4+z3RARERWeZ555RroRX8BDDz1krMCpibom2jpYVhNr9cknnxiJeJs2bYzkXRswZ8+ejYcffhh///03Xn/9deO4jh074quvvjK+Pb/uuuuMsnnz5hlJ+enTp7F27dqLOYuWa9dm+2/e/Q0TdMpVixYt8uwWVLAWfhsyLyQYkXZiX7Zjz5j9lgwWK4JjyyGkVCUEl74SN3XrjKQDm1A3LhO1whNQI+YwLo89iJ61jyHIEooMmwXLTsVgyuYL+H3TKWw7kSll2ety2223OdRv48aNaNzYYbIgIiIiukS64uaqVauMbrLq9ttvNxo/P/jgAyNBP3z4MB577DH06dMH33777cXzHpHW88cffxxvvfWWkbzXqlUrW1Ju/7hHjx6YP3++0Vpvn6Dr3/b4+Hi/vYfBRV0BCmzutvBr15VGjRpJHp+J9DOHjcCuP9G8dzs89dk/i8Zaw6Ikaa8syXsVRJe5Ak0rhKFrhTO4oeR2vHH1WYloXMgMxrqEEpi9W77qWn8MK/ddQKaTzxpNmjTJ9rxt27b4/fff3RppT0REVFCjft2MLYcSfGqaxYYVYzCih/wtvkSaaGcl56pSpUpG95UdO3YYz3/44Qcjib///vuNLir2ekji/d577xmJtyboeq52O9bkW2lXmGXLluHNN9803is9bujQoUYXmHXr1mHIkCGXXP+ixASdfIoO/nTWaq/dVrLKs/qcpR7aDl33dbYZ/ylVGfXqXoZ2taLQvKwNV5Q4gucu243nm4fgdEYJLDhVBrP+TsfExXtwNinV6evrtWNiYi4+HzVqFIYPH27MB09ERESu09bynEqXLo29e/denHoyaxxgbo4ePXrxsXZz+eyzz4wpKLVlXpN0LdMk/7nnnkNqaqqRL+jYQS33Z0zQye9oK3duXW/sv/rSAawVa9RBl8al0blKGq4tvQ8948/gv63K49+rymHyrNVGS31eRowYYUTWL5XNmzejXLlynv2BiIio2MqtpToQ5kHP7RuArL/hWduvv/4aFSpUyDfJ7yhJ98cff4yFCxcarecVK1ZE/fr1jQT9/PnzWL58uZEH6Otec801Hv5pChcTdAoo+o/XWfKufeEr9O2Hay+vjbHXpWBSm+0Y3bIWZp0oi5+2pGDBqu3ITD6X57VPnjx5cW72yMhIbN++HZUr6zpZRERE5C7tsqLKlCmTZyu6/d94Hfyp3Vk0QdfnqmnTpkZ/cy3X/ug6hk5XbvVnBV/akciPaD94XRhq3ootqPPCVmxu+RLOBZfEA+W2Ym7HvzDjyctx38B7Ed/mJgSXqpTv9fSTepUqVYxfFGXLlsWxY8cK4acgIiIKrL/NOruLflOts73kdPbsWaN1PIsm8jr4c9q0afjzzz8vJuj6t1inWP7++++Nb7qzyv0ZE3QqdnShpkbdB6PFyBUIGb4Hh1v8G9eE7cBnFX7C9i5/YuIjzTFwyCMo07KbMZtMfo4fP250e9FfEG+//TbS0/NedZWIiIhgfAutUy+uXr0aDRo0wH//+1+MGzcOL7/8Mvr27Wt0e9GZXuxp8q2DTPVvrX0iro+1T7t+i84EncjP6WJKFW4agfBn9gB3/4ywBtejR9g6fFxqPObecBDNBr6GmDa9jDnbXfHvf//bGFCq/dXtB7YQERGRI13QSPuUN2/e3OhfrjO/vP/++0ZiPnr06ItdS7NkjTPTvum68FHOcv0b3K5dO79/q9kHnUgFyT+FWh0RKYGMdNi2TkWjnx/Gn+Vfwsz4y/Huda9j2YyfkLhhlrE/P6dOnbr4S2Xnzp3GFFFERETFgS5GpOGMzrKSk87MlrVwUX569OjhdKyZ9me/1LVbfAm7uBA5SdYtjW9F0KCFCG77CG4I24RZwf/Ggu5HMfT55/DsxMUIiS3r8vumK6Zp95dt27bxvSYiIiIm6EQFFl8P6PIirI+vhe26kWgQm4ZXM9/C/VsGYMzgm/Hyj8uwaOUaly+n/euYqBMREVF+2IJOlJ8S5WFp/wQin5Bk/I6vUaFydTwVPBl3r++LH3/5Cc/8sAZHEy7gxx9/dCtR37VrF997IiIiYoJOVGDWIKDhzQgfKP3QH5iHMqVK4+3g9zBoQy+Me/0/OBZdC+dT0rFx40aXLqf90nU+9Zwj1ImIiKh4Yws6UUFUvhxhQ1YAfb5FfMUaeMb6JW5b0gOvvjoSO1JjjWWGdfrFZs2a5XkZnfdVV0Jr3bo1kpKSeC+IiIiICTpRgVnl8239GxH+0GzgvlmIiK+BkRnvodzPd2D0TysRE1cK69atM5LwBx54IM9LrVq1CtHR0Vi8eDFvCBERUTHHFnQiT6h6BaIfmY/MG99Gm6BtaL3uWdz3/i9YvuskwsPD8emnnxot5LpqWl7at29vLFes0zQSERFR8cQEnchj/5qssLa6D9Yuo9E1aBXGn70Xaz7/F177bbMxN6v2N588ebKxgJGuPJqbEydOGAsdTZ8+nfeGiIioGGKCTuRpbQcDDy1CerN+eCR4Khou/TdGTFmH5LQMY3fZsmVx5MgRY0nivHTv3t1YeOHcuXO8R0RERMUIE3Qib6jQFME9P0Rmp9HoHrQc3dc+hPvf/RFbDiVcPKR+/fpGy/q0adNyvYyuQhoTE4MVK1bwPhERERUTTNCJvPkPrN1jwK2fokXYQbya+BweGDsLy/4+me2YG2+8EWfPnsWVV16Z63XatGmDTz75hPeKiIioGGCCTuRtTe9A8D1TUCnoND4OeRODxs3Di9O2IDPTdvEQbSVfunQpFi1alOtlHnroITRp0gTnz5/nPSMiIgpgTNCJCkOVVrD0/BiNbTswO2Y0Ni6ZgdHTtxhdXOy1a9fO6HOeW2v6pk2bEBUVhW3bthVGrYmIiLxqwYIFxuraX375pUvHd+jQAdWrVw/4u8IEnaiwNL4Vln4/Ij7CholhY7Bz2VSMnrY1W0u60vnQtTV9ypQpuV6qQYMGmDNnjrdrTEREREWACTpRYap5DSwPL4WlXH18Gv4+ti6bhqE/bkB6RqbDoTfffDP27NmT66U6d+6MCRMmeLO2REREVASYoBMVtrASsPT9DmGlq2B82Cs4u3YKBn+7Finp/0zDaK9atWrIyMjAo48+6vRS/fr1w3vvveftGhMREVEhYoJOVBRiK8Ny/ywEVWyOsWHvI3Prr3hs4lpk5OjuoqxWKz744INc++c9/vjjeOmll7xdYyIiIpfponvauFSlShWEhoYaW31+8mT2mcxyc/r0aTz44IMoU6aMMfZK+56vXr262NyB4KKuAFGxFR4L3PU9gsbfho8PvYNXtx3G098FYcxtTREeEuRweP/+/Y1fcNddd53Dvueeew7JyckYPXp0YdSciIgoVzp1cNu2bY21PO677z60aNECa9euxUcffYR58+Zh5cqVKFGiRK7np6Wl4frrr8eqVatw9913G1MNr1u3Dp06dTJW2i4OfLIF3WKxdJXYLrFTYpiT/SUlfpbYILFSonFR1JPokkWWAgbMgKVRTwwLmYSGm15F/89X4kKqY3cX1bFjR2zevNnpvhdffNFobcjMdOzPTkREVFhee+017Nixw/j2d9y4cXjkkUfw6aef4v333zdmIdP9efniiy+M5Pz555/H119/bZyva4GMGDECf//9dyH9FEXL51rQJdnWpsMPJTpLHJBYJWVTbTbbFrvDnpFYJ2U9ZV9983jHZkUifxASAdz+ORAVjwdWfoyD++Mx8BsrxvVvibBgx5b0hg0b4tChQ6hYsaLDPv1FqJGamoqQkJDCqD0REV2K36Qd8shGh+KIjHQgqIjStPJNgBteKfDpP//8M+Lj5W/ZwIEO63mMHDnS2J/XN75TpkxBUFAQnnzyyWzlDz/8sJGkFwe+2ILeWmKnJN+7JFLl8SSJm3Mc01Birj6QY3RC6OqSqJcr3GoSeZDFAnR9GajfHSOCv0adXd/gkfFrkJQiv6CdqFChgjFfem60vx9b0omIqCjs3r0b9erVQ3Bw9g8Y+lzLd+3alef5u2S//p3TRfzshYWFoWbNmh6vry/yuRZ0UUliv91zbUW/Iscx6yVulVgsibkm9NUkKkscLZQaEnmDVVrLb/8C+PF+PL/1G7yxMwW3/a8fJj10JeIiQx0O1/nStZ9er169nM6Zri0NH3/8Me8VEZEvy6Wl+oI0wuTVTzuQ2Ww2Y/Gi3PYVB76YoDu7Iznvhv7X/K7cvHWy1e+F1ko4NDXKfv1uxfh+Rb9q0dWqyD8lJiYWm/tnib8X9U+exVPHvkO1U0dxz/sP4LFWJRBsteQ6i4sm6YcPH3bYN2nSJJQvX97bVc5Tcbp3gYj3z7/x/vmO2NjYPL/5zEmn2HXneF+iK31qX3OdicW+FT09PR3bt2839mf9bOfPnze2OtFBVlm1atWMwaQHDx7M1oqekpJitM7HxcUV6XujdfX23zVfTNC1xbyK3XNtGT+U49NTgmwGmEm4Zi27zUCO4z6RjYZ+pWLTKXrIP+k/hGJ1/zpcC/zxCnr98SqqJx3GtJMfYOTtrXNtUVA6mMZZn753330Xjz32mDdrm6did+8CDO+ff+P98x1bt251q0X8nB+3oN96660YM2YMJk+ejEGDBl0s11lcdPpFLcv62SIjI41teHj4xbLbbrsNs2fPNgaGjho16uL5n332GRISElCyZMkifW+0rs2bNy92CfoqiTqSiNSQ7UGJPhJ97Q+QfXGyOW/2UX9AYqGZtBMFBqsVuPYZoExdXP7jg9i/4QV8UPJtDOlUN9dTXnjhBaNlYfz48Q4t7PqLbMAA4zMtERGRVw0dOhTff/+9Me/5mjVrjGRWp1nUBFv7oOv+vAyQv1eanGf9XbvyyiuN8/WatWrVMlriA521qCuQkyTa+q4PlpgpsVXiOynbLEE+RFAAACAASURBVEn5IA3zsAYSWqYDRG+QeLxoakvkZU1uB64ZiluDFsO6YDQ+W5T3wJrcFjPSeWh11DwREVFhdOdZsmSJMWvLjBkzjG9xdast54sXL8639Ts0NNRoQde/XdOnT8dTTz2Fv/76yyirXFk7VgQ+S3HpbK9dXLTfE/mnYv01bWYGMn99Ata1X+GT9BtR9rbXcUtzHUud+9eiOUe+Z1m4cCHat2/vrZo6VazvXQDg/fNvvH++1cWlQQNtXwz8Li6BbqsL91IakVdLjt0yYFrQichxdhfrTe8io+WDGBg8HTt+HoO9J5NyfZv0F7qzAaPq6quvNlohiIiIyHcxQSfyBxYLgrq9igt1bsLT1vH4+au3kZ6R+4qhOnOLruLmjPb/O3PmjLdqSkRERJeICTqRH7WkR/Qeh5OlW2Lg2fcwcYaxVleuateujRUrVjjdxy4nREREvosJOpE/CQ5D6f7jkREcjiv+/BeOHDuW5+GtW7c2BtjktH79enzzzTfeqiURERFdAiboRP4mpgKSbx6HGjiM41/fm+/h3bp1wzvvvONQfs8992D+/PneqCERERFdAiboRH4ovmkXLKv+MJokLsFvv0zM9/jcFirq2LGjMccsERER+Q4m6ER+6qq7nsPJ4LKotPpVzNigC/DmTlcgzW1Z5Jo1ayI1Vdf8IiIiIl/ABJ3ITwWFRqDEjS+iqXU3Nnz/MlbtOZXn8dHR0cbcrc4MGTLEG1UkIiKiAmCCTuTHQi+7A6m1b8C/gybjky+/wIHT5/M8vn79+pg5UxfpzU6XVNZV34iIiKjoWYu6AkR0CaTrSmjPD4HStfCe7RV88sOMfE/p0qWL0fc8p3bt2uFYPrPCEBERkfcxQSfyd1GlEXrvVGlOj8RN+1/F9PUH8z1l2rRpTsubN2/u6doRERGRm5igEwWCEuUQcsMYtLT+hTU/v4PDZy/keXhERAS2bdvmUH7o0CH89ttv3qolERERuYAJOlGACG7eFxcqt8PjtvEYPXEeMjNteR5fr149fPihdI9xMm/6vn37vFVNIiIiygcTdKIA6o8e0fM9RAWlo9vB9/Dh/J35njJo0CCn5dWqVfN07YiIiHzOlClTMHLkyKKuhgMm6ESBpHQtWK95Gt2DlmP13MmYuv5QnodbrVacOHHC6b4ZM/IfcEpEROTvCfqoUaOKuhoOmKATBRjLVf9CZpl6eC38K4z6cSV2n0jK8/jSpUtj0aJFDuU33ngjDhzIewEkIiIi8jwm6ESBJjgU1h7vomzmMTxm/QGPTVyL1PTMPE/RKRZr167tUF6nTh1v1ZKIiAJYSkoKxowZg0aNGiE8PBxxcXHo0aMH1q5dm+243r17IygoCAsWLMhWPnPmTONb3nvuuedimU5u8MgjjxjXLFGiBCIjI3H55Zfj008/dVqHhIQEPPvss2jQoIFRB22Q0r93kyZNMvZ36NABX3311cUVt7Piyy+/9ORbUSDBRV0BIvKCalcCl9+Lu9d8g+8OXYk3ZpXGM90a5HnKn3/+afwCtZecnGz80tRfYkRERK5IS0tD165dsXTpUtx9990YPHgwzp49ayTSV111FRYuXIiWLVteXChv1apV6NevH9atW4cyZcrgyJEjRmKuDUf/+9//Ll5X/x7pud27d0eNGjWQlJSE77//HgMHDjS6aw4fPvzisWfOnDGS8c2bN+P222/Hww8/jIyMDOMDgk413KdPHyN5z8zMNL5F/uabby6e27Zt26K/0TabrVhE3bp1ZUP+av78+UVdBf9z/pTN9lpt29ExTW2X/2eC7Y/tx/I9RVoSdOoXhzh+/HiBq8F75994//wb75/v2LJli1vHS+uvl2rifW+99Zbxt+P333/PVi5Juq1KlSq2a665Jlv58uXLbSEhITZJvG2SRNs6depkCw0NtUnDUbbjEhMTHV5Lj9frxcTE2FJTUy+WS0Ju1OHjjz92ek6W/v37G8d5+l7KNf+8lLyVXVyIAlVESaDXF4iXri6TIl/D0O/W4OyFtDxP0RYMZ+SXpTdqSEREHqItwUePHjWSu6I2fvx41K9f3+h+oi3bWSEJNDp37ozFixfjwoX/X6/jiiuuwIsvvmi0bF999dWYM2cOXnnlFeN8e1FRUdm+4T158iROnTplrJCt3Vmy1vfQ90K7sWjXlgcffNChftp1xtdZi7oCRORF1dvBcvOHqJ25G1dfmIO3Z/+V5+H6S0t/wee0fv16bNy40Vu1JCKiS6AJ6bXXXovKlSsbXRL1eVHaunWrkSzHx8c7xOeff250Nck5g9jTTz+N9u3bY8mSJUbC/a9//cvhutKCjqeeegpVq1Y1FtzT7jB6Te2qok6fPm1s9dr6+LLLLjP6lPsj9kEnCnSNegLLPsSzx35Gm2Vt0btVFTSoEJPr4WXLljUG4dj3+1NNmzY1Bv3I147erjEREblBuiEa/b3T09ONrT4vV65ckb2H2orfpEkTSFeXXI/RxNrenj17sGHDBuPxzp07jWRcB4La69u3r9HKrn3OtaW9VKlSCA4ONqYFfvvtty9+MMn6FsFfk3PFBJ0o0OkvqM6jEPvljRgYPgcjfimPyQ+1yfMXl/6iy5mgq5deeskn54slIirOtGFFBzZqcq5bfV6UdAYw/ZDQsWNHl7qT6AeLO++809i+9957ePzxx41BndpVxn7QpybnOuh07Nix2c7XLjE5k/+SJUsag07z46tJPLu4EBWTri6o3RmPBv2CHdJKkd8CRtpKvmPHDofyF154AceOHfNWLYmIqAA0yZw/f76xdoXOdFLUSafOwKIzseTWgp6zK+Vzzz2HFStW4IMPPsCQIUPw5JNPYsKECRenQFQ6FaPK2cf+8OHDGDduXLYy/VCgCf+WLVvw2WefOby+/TWio6ONrfZl9yVsQScqLrqMRsjY9ngnZgKG/1YWXRqWR0ToP7/wnNHprXSQT9agmyzaIrJp0yZv15aIiNygSWlRdmuxpy3gs2fPNvqVz5s3z/i7ERMTg3379mHu3LnGnOT6gSKr9fu1114zuq/ce++9RtmYMWPwxx9/GNMz6jcC2iKv3V20b7q2qmv/81atWmHv3r34+OOPjSkXdcCoPR10qq/9wAMPYNasWcaUi5qY6zSL2lKfNa1imzZtjA8G2rVTF+gLCQkxBq3qNf0yQZdPZ3qurmxSUiIkr2PlDVlY0NchIg8p2wCWDsNwzbzRqJV6Fb5cWh0Pd6iV5ynLly93mBtd55TV/oHOFjYiIiLSJHf69OlGV0lNhEeMGGG8KRUrVkTr1q3Rv39/47l+I6tdVmrWrJmt20qInD9x4kQ0b97caAnXrjv6za4m58OGDcOvv/5qtK5r4q5dL/X4AQMGZHvjtYvLsmXLjGT/p59+ws8//2wk+Q0bNjRa6bPo9TVp11lfdE517cf+xRdfFHmCbnF3Oh5JzLXGL0vcJBHmwik6HWSRt9TXq1fPtn379qKuBhUQF8vxkPRU4L3LsDWlNO5Kfx6Lhl6LqLC8/3nedddd+Pbbbx3KXf3dwXvn33j//Bvvn+/QmU102j9XnTt3zmGQJPnPvZR8ebX8nfxnNSZv90GXF9Mms5USvSTCzWLtkLovj9hf0MoRkYcFhwJth6BBygZUP78Jk1bl/8/TWf89pSu/ERERUdEPEh0tUVrioMTtEmHy6aCCRI28wt1KyQeBrhLbJXZKDHOyP1biV4n1Epslsn+vQUS5a3EPEFEKw0rMwA+rD+T7TmlfQe0LmJN+TXn+/Hm+00REREWcoHeU0O+175TE+yeJdA/XR5NvHbX2ocQNEg31taRMt/YeldB1VpvJtoPEm3IMJ2cmckVoFNDmEbROlS/DjmzE1sMJ+Z6ii0c4M3ToUL7nRERERZyga2eoC5IYL/FwPey1ltgpr7FLQjrMYpLEzTmO0Q8JJSz/zCOk8+Po3Dge/7BAFLBaPwBbSBQGBrvWiq7/1JytMPrhhx86rAZHREREhZug79NzzMTYWypJ2HeMPWCW2ftAQnvn62TOuv7445LMF+26tkT+JKIkLC3uRo+gZZi7agMSktPyPUUXvnA2KOaWW27xRg2JiIiKLXdnV9HW7P9KXCeRfdkmz3GW/OecLuJ6iXVmlxudJ262fGZYJEl6tu/qpWygbAZmrSqlo9nJP+mSv7x/nhVhu0y+rspAz4zf8eLEUuhWI/9eYu+//77Tldl0vlud5soZ3jv/xvvn33j/fEdsbKwxM4urMjIy3DqeCk9ycrLXcxJ3E/RXzOkVP5bkt5MkxLu9UCdtMa9i97yyRM5lD3VQ6Cs6f6NsdSCp1qO+OcPMRbL7E9l8kjXNYocO2l2d/BGnCvOSs1MxYMc8dNp3O57p3QFxkfkn6TqnrC7+YK9atWrYs2eP0+N57/wb759/4/3zran5dNVKVzshcJpF36Spp06eoHO0+1KCfofEFxKjJDbKf2Q/yFbnWjuXzw/ztRuvoderY863rrPF9JHo66SrjbbiL5LjdNmsehK73HgNIlJtHkbM9hm4Nu0PvDW7Ol64uXG+78vUqVMvLo2cRVdz0wS9evXqfF+JiJzQperT0tKMBXfIf6XJPdR76W3uJuhfSmirddbHv7vNyI/LCbrODCNJ92B5OFNC34HPpUynUhxk7h9rTvf4pZRtNOvyHynnSDUid1VvD5RvimGnpqL98jbo06oqGlaMyfOUqKgoYxU4XcXNXpMmTfh1LBFRLnTRoYSEBJQpU4bvkR9LkHtYGAtIuZugL5Rwb+nRApBke4ZsZuQou7gGrDzWLi9dvF0PooCnX7V2ex1xn3fFM2HfY+SvFTF5YJt8v4LVJZlzJuja13XXrl3Gks1ERJRdqVKlsG+fdgAAYmJijHE73p1zgzzZrUVbzjU5P336NKpWreqpS3smQZcKshM3UaCpKgl5i3vQe+0EvL27O37dUA03NauY5yna/27YsGF45RUdlvL/6tevj9RUnR2ViIjshYWFGYndqVOnjC6BOgg0v4GI+ruWfIN2a9GWc72Hei99rQWdiALRVY/DuuZrPFlyIcZMr4DODcohIjTvPnajRo1ySNC1heGvv/5C3bp1vVlbIiK/pIldhQoVjHBlgK+3ByJS4MyDTkSBqHQtWOp1w+2ZM3Eh4QQ+X5L/BE060EkXKsqpXj0ds01ERESFnqBLv6maEkMlJknMNWOSWcZOqET+5tpnEJyagHfKTsNHC/7GqaT8u6oMHGgsM+BAW9GJiIiokBJ0Sb4jJHRucf0L/LI59eK1Ztxhlv0lx4zVYwtWLSIqdOUbA60eRIeEXxGbegRfuNCKHhwcjIkTJzqUsxWdiIiokBJ0i8Wix/8icb95rs6mMkHiVTMmmGW670GJKXIOhygT+YsrH4EFNjxdYQO+XLoHCclp+Z7Sq1cvp+VsRSciIiqcFnRdwbOTRIrEQxJVbTbb3RLDzdA50XXuGZ2zPNU8Vs8hIn9Qsrr8C74SXTMX4Jwk598s2+vSyPbJkyc7lLMVnYiIqHAS9HvMedAfk2T8UwmHOdG1TGgXmMcktPW8f8GqRkRFolkfhJ/ZiUHVDuOzxbtxPjU931Nuu+02p+U7d+70dO2IiIgCnrsJehMJ/c47+wolzn1lHqvnEJG/aHIHEF0ej2IyTiWlYPxy11rRP/lEP5dn17BhQ2/UkIiIKKC5m6DroM/z0kKeb8dUOUa7uCSZ5xCRvwiNBK5+CiWOrsJDVQ4YM7pod5f8DBjg2JtN50XnwkVERETeTdB1AGisxWKpnd+BcoyuVBJnnkNE/qSF9GaLKotHw37D6fNp+HzxHpdmdHn++ecdyrds2eKNGhIREQUsdxP0OWa/8o8lAc91/Vlz31gJ7aM+u+DVI6IiERwGtHoAMQcW4J7ayRi3aBdOuzAv+rPPPutQpstZ79u3zxu1JCIiCkjuJug6lWKyRAeJDZKID5KoL1FCoozE5RJPyb4dEteYx77m2SoTUaFodT8QEomnwqYgMTUdYxf+7dLqoo8//rhDeZs2bbxRQyIiooDkVoIu/cp3mYsRnZfQbi66zvdmiTMSRyVWmkl8JfOY3uY5RORvosoAVw5GzN9T8Vi9BHy1dA+OJehn7ry9+qr+Csju8OHDOHLkiDdqSUREFHDcXklUEu5psmkm8YVEgtnlxT7OSnyux5jHEpG/uuoxSdTj8XDqF0jPyMQH8/OfNjEsLAx9+/Z1KL/xxhu9UUMiIqKA43aCrrRVXOJ+iZJmS/qVZtSWslISD7DlnCgAhJWQDm3DEX5oBUbU3YuJK/dh/yn9cixv48aNcyhbs2YNzpzRL9uIiIjI4wm6k2R9hRnszkIUaFr0B8rUxZ1nxyHEkoF35+oQk7xFRESgY8eODuVDhgzxRg2JiIgCyiUn6EQU4IKCgU6jEHz6b7xeYy2mrD3oUl/0X375xaFs/PjxSE7O/1wiIqLiLNcE3WKxVDWjgpMyt6JwfhQi8pp6NwBV2+L6U9/CakvDt9LVJT/R0dFOVxJ94403vFFDIiKiYtGCvtsMnfs8Z5k7wW4vRP7OYgHaPYHgpMN4suIWTFixD2kZmfmetmTJEoey//73v96oIRERUbFI0LNmZbE/JueMLa4Eu9EQBYLanf7pi54xFcfPJWPhX8fzPSUuLg4hISEO5du2bfNGDYmIiAJCrsmzzWazmtHASZlbUTg/ChF5lVX+KV/5KGLObEGXyB34cc0Bl05z1s2lQYOLv1aIiIgo55/cHM+JiHLXtDcQWQZPlpiNOVuO4VRSar7vVnBwMMqWLetQvmfPHr7TREREl5qgm4M+K7lxfEUOEiUKICERQKsHUO/sElTJ3I8Jy/e6dJqz1UVbt27t6doREREVyxZ0bfJa6cbxOkKMg0SJAokk6AgKw7OlF+CrZXuRkp6R7yn9+vVzKDt+/DiOHTvmjRoSEREVuy4uOvDTm8cTkS+Ljgea9cY1F+YgLfEkZm0+6lI3l3vuucehvFevXt6oIRERkV/zdh/0SIl0L78GERW2Vg8iKCMF/aJX4ee1B106ZezYsQ5lCxcuxPnz5z1dOyIiIr/mtQTdYrHUlk0ZiSPeeg0iKiLlmwBlG+Gu8KX446/jOH4uJd9TIiIinPY7f+mll7xRQyIiIr8VnE+SfbNsNOzFSvnneZ0mESfRznw+391KyfW7yuZdiSCJcTab7ZUc+5+WzV12P4PO2RYvx51y97WIqIALF112JyrOeg7VbQcwaeU+DLmuTr6nzZ49G7GxsdnKxowZwySdiIjI1QRdXCZxb46yCCdluflbwq1lAyX51qT8Q4nOEjrR8iopmyrJ95asY+Tx67J53Ty+h2yeYHJOVMia3QnMHY1nSi/E0KU18ODVNREeov98cxcTE4MqVapg//792cp37tyJ2rX1SzciIiLKL0FfkOP5CIlEiTfzOEfX/06Q2KznS+Lsbh90/Q58p5y3y0zAJ5mt+BcT9BwkS8BEN1+DiC5VlPRga9ILHTb9iLSkHkZf9DtbV833tJkzZzosXtSsWTMkJSXxnhAREeWXoEuS/IdsNGAmy0aCLuWjvPju6Tzr9s1r2op+hbMDpT46CFW7wwz2Yn2IKDdtBiFo3XgMjl0qCXpVlxL0+vXrO5TpQNEjR46gfPnyfK+JiKjYs0iy7fKbIAlxNdlkyDmurfFdAPIaOu/a9fIaD5jP75ZNa3k+xMmxvWXTT/b1yOVaA2Wjgfj4+Mu/++47b1WbvCwxMRHR0dF8n31Qs3XPwnbuCFonvY1Xr4lG6Qhrvvfu4MGDRkJuLyoqymnyTkWL//b8G++f/+K982/XXnvtaslPW3qri0s28kKuLRt4aTT5r2L3vLLEoVyO7ZNX9xap7yey0UC9evVsHTp08FQdqZAtWLAAvH8+qvwzwKS+6Gz9Eyej+uO2q2vle+9SUlIQHh7ucKlz587xg5iP4b89/8b7579474q3S5pmUVqoVSmJKhJVcws3L7tKoo6cV0Mi1EzCpzp5bZ0K4hqJXy7lZyCiS1RXepnFVcOQyNn4ZV1un6WzCwsLQ7t2WRM9/b+RI0fydhARUbFXoARdkuPbJOaaA0aPS+yR2J1LGIM9XWUOKtU+5TMltkp8J2Wb5fUGadgd2lNiluzjyDKiomQNAq4YhIbpW2A9vA47j51z6bQZM2Y4lL35Zl7jz4mIiIoHtxN0SZI/ko125r7WnHJR5z3PK9x+DUm6Z0jUlaglYaxiItuxGnbHfCmhretEVNSa90NmSBQGBP+OqS62opcoUcKYcjEnnXKRiIioOLO623Ium4cktNVaB2+WMncdMfuzVzLnSP9L4qREF0mivbZaKRH5iPAYWC+7E92DVuCP9dv1A7RLp+mUizm1aNHC07UjIiLyK+4mzzqziv7lHSZ/gCdInMnaIY8zJQ5LfC1PLze7tvwsST2nZSAqDlrej1Ckoc2ZGdh8SJdCyJ+zWVt0oOiJEyc8XTsiIqKATdCzmrbG53Uds1+49iOPkhhesKoRkV8p1xBpFVuhZ9BiTNtw2KVT5AM8Hn30UYfy/v37e7p2REREAZugx0mckwTcvnksVcJhgmo5RmdjSTL7qhNRMRDS5FbUt+7H6jWrkJ6hiwrn74033nA6gDQ1VX+1EBERFT/uJug6Y0vOyYtPSURIS1gZJ8cHSZQtSMWIyA81+GfNsBbnl2DetmMunaLzoTdu3Nih/OuvtbccERFR8eNugr5fIkSScfv1uNeb2+vtD5RjrjaT+dMFrx4R+ZW4KrBVaI5bQ5fjm2U6+6prFi1a5FD22GOPebJmREREAZugLzC37e3KfjCnU3xLkvJeErrI0K3y/GtzQOmsS68mEfkLS4t+qGfbjaRdy3DkbLJL58TFxaF8efvP/cCFCxdw6JBrUzYSEREV5wT9ZzMZv8eu7EuJZRLxEpMktkl8L6EriOpUDM9fejWJyG807YPM0Bj0D5qFqesPunza8OGO48lvuOEGT9aMiIgo8BJ0m822UjYlJO6wK8uQTReJ1yX0O+10cw70iRJtZP9ej9WWiHxfWDSsze9Ct6CVmLN6u8unPfzwww5lGzZsQGKiLlhMRERUfBRklc8kiQtOyv5jrvwZJlFW4i6J3Z6rKhH5jcvuRIh8Vq99Yg4OnnNtNpeQkBB07tzZoXzUqFGerh0REZFP4yqfROR55ZsivXQ93Bq0GMsO65dqrvnpp58cyt59911P1oyIiMjnMUEnIs+zWBB8WW+0tG7HvkMHkZmp48XzFx0djSpVqmQrS0tLw549rs8IQ0RE5O+Cc9thTpPoETabbaGnrkVEfqJJL2DuC7g6bQn+3HszWtco5dJpEyZMwNVXZ//1065dOxw4cMAbtSQiIvKfBN2cUtG1Zq+82fJ5HSIKRHFVkVGlLW7dtxjj1h5wOUG/6qqrHMoOHjyIhIQExMTEeLqWREREftXFZZ+HQhc3IqJiKEi6udS0HMbfG5YhNd21waJWqxU333yzQ/lzzz3n6eoRERH5V4Iu3VKqS9TwRBTmD0REPqTejfIVmgVXpi/HvG1HXT5Nu7nk9MEHH3iyZkRERD6Lg0SJyHui43E2pj5uCFmL7/50vQ95VFQUatasma1MPuxj165dnq4hERGRz2GCTkRedbLMFahn242/tm/GkbPJLp83ffp0h7I2bdp4smpEREQ+iQk6EXnV8fh/kupu1uWYsu6gy+fVq1fP8VrHj+PMmTMeqxsREZEvcmt2FYvFMq8AryHfTNuuK8B5RBQAkiMqAJVaos+RFRiy/hAGXVPL1d83uPPOOzFx4sRs5U899RTGjRvnjaoSERH5BHenP+zg4nFZ0zNa7B4TUXHV9A7UPDgUaYc34+/jzVErPtql0zQRz5mgf/XVV0zQiYgooLmboA/IZ3+sRCuJ2yTOS4yUOFeAehFRIGl0K2y/D8ctQUvw6/qO+Fenui6dFhkZidq1a2Pnzp0Xy9LT043nWk5ERITinqBLV5WvXPxqepRsZkncK9GuAPUiogCbzcVSqyN67V6OO9cdwOPX1TG6sLhixYoVKF26dLayli1bsi86EREFLK8MEpVEXpu7Bkm0kBjujdcgIv/r5hKfcQylTq7B1sOuf7FWqlQplCtXLlvZ2bNncerUKU/XkIiIKOBncZktoXOq9fHiaxCRv6h/I2whkbgjeKFbs7mo7t27O5Q9+uijnqoZERFRsZpmUdf2ruLl1yAifxAaBYu0ot8UvBzTV2zB2QtpLp/6/vvvO5T98MMPnqwdERFRsUjQ20pESiR48TWIyJ+0vA+hthR0SZ+Pb5btcfm0iIgIh0GhWYNFiYiIAo3HE3SLxRIs0VMeTjCnWJzj6dcgIj9VoZkxJ/oDEQswccU+Ha/i8qk6WDQnHSxKRERUrBN0Sbx35ROHzH7n+t1zVYmTEv91t1Jyna4S2yV2SgzL5ZgOEuskNkv84e5rEFERaXU/KqXvR+WEtViz77TLp3GwKBERFRfutqBXzyfKm9dMlfhO4gppIdvtzgtIsh0kmw8lbpBoKHGnlDXMcUycbP4ncZNcv5Fse7n5cxBRUWnUE7awWPQLmYtf1x9261QOFiUiouLA3QT92nyivUQTiRhJnPu4m5ybWkvslHN3SWiiP0ni5hzH9JX4Sfbv0yeyPVaA1yGiohASAUujm9EpeD1+X78P6Rk6ltw1HCxKRETFgVsJuiTCf+QTSyQ2S7g+PYOjShL77Z4fMMvs6TKEJaUlfYHEaol7LuH1iKiw1e2KiMwk1LiwEct3uT6fOQeLEhFRceDWSqKFxNnygjYn9b5c4jqJCIllkqQvlw8Gf2W7kMUyUDYaiI+Px4IFC7xQXSoMiYmJvH8BdO+C0q1oawnB9UFr8PHvzZB+MMzl640dOxbr1q3LVjZt2jRcdtllHqkvZcd/e/6N989/8d4Vb5eUoJt9wfWvYrxZdFxinSTKZy7hsgdyzJ1eWeKQk2NOyOskyTZJ6rFQ1pijPgAAIABJREFUts0ksiXosv8T2WigXr16tg4dOlxCtagoaYLH+xdg9+7INeixdx3eOpmJK9q2R0SoDj9xTb9+/XDkyJFsZcePH0eZMmUutbqUA//t+TfeP//Fe1e8FWiaRUmI20rMkocnJOaa/cQnmY9PyL7fJa4sYJ1WSdSR82tIhJorkU7NccwvEu3NKR11rvUrJLYW8PWIqCg0uQOl0w6jQcpmTF3v3sqivXo5jgu/7777PFUzIiIi/0rQJSH+t2wWmt1L9PxMczrFU+ZjLesisUiOfcLd60urd7psBkvMNJPu77Rfu1xrkIZ5jJb/LrFBYqXEOCnb5O5rEVERatADtrAYPBC9BF8v2+vWnOivvfaaQ9mvv/6KzEzXB5wSEREFyjzo18vmDfO8BWYiXkL+sJaV0G4u0RKdJeaZx7wh5+hzt8i1ZkjUlagl8ZJZNlbD7pjXJRpKNJZ4x93XIKIiFhoJS+NbcW3mUvx96DjW7HO9Z1x4eDjq16/vUL5s2TJP1pCIiMgvWtCfNrfjJSm+TmKOhC5MZJDHKRJzJTrpMeaAz6EeqisRBZpGtyI4IxldwzZJK/oet05dunSpQ1m3bt08VTMiIiK/SdB1Xe1Mu0Q9L3qMfmfdyt1KEVExUe0qIKIU7iu1ETM2HsaJxBSXTy1ZsiTKl9e10f5fQkICTpzQoTFERETFJ0HXFvGz0kJ+NL8DzWMuZTYXIgp0QcFA/W5olLQM1owUzN6S76+WbDhYlIiIApG7CboOzoyRfuUl8jtQjonRY81ziIica9obQann8EDMCszanH3qxPxwsCgREQUidxP0jyR0suJhLhz7H/PY/7lbKSIqRqq3Byo0w32W6Vi68zgSU3QiJ9dwsCgREaG4J+jSbeUr2eiMKcOkhXysRLWcx0hZVYmPzCT+bTnnG89UlYgCkkV6zrV9DKVT9qGNbb10c3GvFZ2DRYmIqFivJCqJt06fqBIkHtSQsn2yzVplpKJENbtjmtudY0/ydpvOo05E9M+c6BElcXfQEny68lr0bK4LCLs3WNR+ZdGswaJcWZSIiIpDFxddr1sj1hwwajET8rZmVLcrj7U73lkQEf0jOAyWxrejg20Vtu7ej13HEy95sOiAAQP47hIRUeC3oAv+xSMi77isL4JXfYoewSsweVUTDO/WwK3Bou+//362smnTpiE9PR3Bwe7+miMiIipawQXog05E5HkVmwPx9TEgYRl6r+6KJ7vUQ2iw1eXBoo0bN8amTZsckvRbbrmFd4uIiAK6iwsRkfcGi0oreu2UzYg5v9ftOdEXLVrkUHbXXXd5qnZERET+k6BbLJYIiSpmRHiiUkRUfOdEt1mCcH/kIkxapePPXRcXF4dy5cplKzt//ny2waNEREQBm6BLIl5KYqTEFnl6TmKPGee0TGKERElPVpSIioES5WFp0AO3W+Zh1Y6D2H/qvFunDx482KHs1ltv9VTtiIiIfDNBl8S7tWy0o+d/Jeqb18iaucVqlj2vx5jHEhG57opBCE9PQM+gJW63og8dOtShbNmyZUhOTuYdICKiwEzQJeHW749/kygvcUbiZYnOEjrdQgPz8SvmvgoS081ziIhcU7UNULYhBkQtxS/rDungdJffudDQULRo0cKh/L333uO7T0REAduCrs1T2nVlg0QD+cP5rMRcie1m6ONnZF9DiY0SpSSe9myViSjgB4s2vQN1U7fAcmYP1u3Xz/uumzfPcW20ESNGeKp2REREPpeg3yihzVn3SSJ+LLeDZJ9Ov3Cf2e2le8GrR0TFUpN/Fh66LXgppq4/5NapsbGxqFhRFzX+f9rF5cCBAx6rHhERkS8l6FUlzkkCvia/A+WY1XqseQ4RketiKwPV26NP+DJMWXMAF1Iz3Hr3xowZ41DWqVMn3gEiIgrIBD1VIlT6lWvLeJ7kEL12iHkOEZF7pJtL+bQDqJK8HT9Iku4OZ/Ofb9++HYmJibwLREQUcAn6NokwiZ4uHKvHhEtsd7dSRERocBNsQWEYGLcKny3ahYxM1weLBgcHo0uXLg7lI0eO5BtLREQBl6B/J6Gt559IC7nO2OKU7LtJj5HQv6gTC149Iiq2IuJgaXgTuqbOQeLJw26vLPrjjz86lL355pueqh0REZHPJOgfSKwzZ2f5XRLxFRKvSAyReErifQmd4eVnc7YXPfZ/nq0yERUb1wxDUGYKhkf/ik+lFd0d0dHRqFatmkP5zp07PVU7IiKiok/QbTab9ifX741nmi3prcxpFN+ReFXiEYnG5r7fJa43zyEicl+Z2rC0uBu3ZMzCsX3bsO1Iglunr16tY9Wza968Oe8EEREF1kqiknCfkLhBHl4toat/LJH4y4wlZtnVckw3PdaTlSWi4tmKbg0KwdMhP+CHP90bLFq6dGnEx8dnK9OBoseO5TpLLBERkf8l6Fkk+V4s8S8JTcZ10aIG5mMtW+zJShJRMRZTAZYrBqK7dSmWrV2HtIxMt07v3bu3S2VERER+n6ATERWalvfJLysbOibPxdyt7g0Wff311x3KFixYYCxeRERE5LcJusViuUziI4mVEtsklku8K9HA2xUkIkLJ6rDpwkUhi/Dt8j1uvSHh4eFO+52/844OnSEiIvLDBF2S8Mdk86fEQImWEnUlWksMllgn+x1XBCEi8jBLi3tQCUdh2TUfe08muXXu/PnzHcqGDx/uqaoREREVXoIuybcm5G+Zx6Wag0B1LvTlEpnmSqGfynE1PVkpuV5Xie0SOyWGOdnfQeKshH5A0Hjek69PRD6o4S3IiCqPh4Kn4duV+9w6NTY2FhUqVHAo37XLvakbiYiIfKEFfbB5zGaJxjabrb1EH4m28lyT9/3myqIPeapCkmwHyeZDCZ0ppqHEnVKm25wWST0uM+MFT70+Efmo4FAEXfkw2lo3Y+2qJUhJz3Dr9EmTJjmUtWypv8aIiIj8K0G/ylwN9CFJgv+23yHP18vmSXPOcz3OU7T7zE65/i5zDnX9q3qzB69PRP6qeT/5hWTBlSlL3V5ZtF27dg5lp0+f5pSLRETkdwl6RYk0s0uLMwvsjvOUSmbLfJYDZllOV0rL+nqJ3yQaefD1ichXRZUBqlyBG0LW4Nf1h9w61Wq14t5773Uo79Wrl6dqR0RE5BHB+eyPkDgiLdlOJx7WhYgkOdaH4R6pzT8szl4qx/M1EtXk9RPl9bvJ4ykSdRwuZLHowFYNY7ESnVqN/JMuLsP75588fe+qhNRDfWkz2LJ1M36bcw4Rwc5+ZTh3zz33oHFjXew4u3nz5hkJPDnivz3/xvvnv3jvirf8EvSioC3mVeyeV5bI1lQmifnF9b7l8QxJxP8nUSbnyqXy/BPZaKBevXq2Dh06eK/W5FWa4PH++SeP37sT8ivhg69wrWU1LpS6Dje00F8Rrhs2bBhWrlyZreyJJ57AW2/peHjKif/2/Bvvn//ivSvefLHJaJVEHUm4a0iEyuM+ElPtD5Dy8hJGs5lsWps/x8lCrykRFb4ytWErUxc3ha3DhBXuzeaiZs+e7VD29ttvIzPTvRVKiYiIijJBLydJcEZuYXY/yeuYdHcqJK3e6ebsMTMltkp8J2Wb5TqDNMzDbpfYpH3QZfuehM4sk7MbDBEFKEu9bmhh24wdew9gw4Ezbp0bExODqlWrOpQvWrTIU9UjIiLyeoKuLdWXGm7RbisSdSVqSbxklo3VMB9/INFIoplEG4ml7r4GEfkxSdCt8lm+a+gGfLV0r9unr1mjw1iyYxcqIiLylz7oAwqlFkRE7qjcCoitgocylqLbhqvwfPeGiI3UddNcU7p0aZQtW9ZhisU9e/agevXqvBdEROS7Cbq0TH9VWBUhInKZzrjS8j7UmjsKVTL24ae1BzDgqhpuvYGvvvoqBgzI3gbRqlUrHD9+nDeCiIiKlC8OEiUiyl+Le4CgMPw79g9jsKi7w1D69evnUHbixAkcOXKE7z4RERUpJuhE5L+LFjW+DZ3T5uGodFVZufuUW6cHBwcb86Ln1L17d0/VkIiIqECYoBOR/2r9IEIyLqBv+JICTbk4dqwx7jyb1atX48wZ92aGISIi8iQm6ETkvyq1kGiJB8Pm4reNB3HozAW3To+IiEC7du0cygcP1pleiYiIigYTdCLyb60HokzKPlxp2YTPFu92+/QZM2Y4lE2YMAHnz5/3RO2IiIjcxgSdiPxbo1uAyDJ4uuQiTFy5D2fOp7p1eokSJdCgQQOH8pEjR3qqhkRERG5hgk5E/i04DLj8XjROWopSaUcwfrn7CxctXeq41tnrr7+O5ORkT9SQiIjILUzQicj/tRwAi/xvePxSfLFkD5LTMtw6PS4uzukCRa+99pqnakhEROQyJuhE5P9iKwP1b0SXlFlITErED6sPuH0Jnb0lpxEjRiA9Pd0TNSQiIvJ+gm6xWKIkypkRVdDrEBF5arBoSMppPBK/Hp8u2oWMTPcWLipVqhQqVqzoUP7999/zBhERkW8m6JKEN5YYI7Fc4qwUJUgcMiNBy8x9L+mx3qowEZFT1dsB8Q1wb9BM7D2ZhN83ub8i6Pr16x3K+vbti8zMTL7pRETkOwm6JNuREt/Iw3US/5FoLVFCd+WIEua+YXqsnPM1W9aJqNBYLMbCRbFntqBHyf0Y+8ffsNnca0UvU6YMypYt61A+b948T9WSiIgoX8F57ZQEO0Q2+peplT6V+EtilsQWs+U8a6LgSAn9brihRGeJehJ3SdSRa7SXP5LsxElE3tesDzD3BQyNmYP2e6ti6d8ncVXtMm5dYtOmTQ5JeufOnZGRkQGrlcN2iIioiBN08bjZKn5S4n5JtKe6clFJynvI5nPzXL3Gm5dSSSIil4RGAS3vQ+Ul7+Cy6J5GK7q7CXp8fDxKly6Nkyf11172VvROnTrxRhARkdfl1xzUV0K/I+7vanKu5Nhf9Ryz1V2vQURUOKSbi8Vixahyi7FoxwlsPqRDZtyzZYt+SQiHVnT2RSciIl9I0GtLXJCE23Et7HyY52gXmDoFqRgRUYHESG+7Rrei6fGpKB2cjO9W7Xf7EtrFRVvRc5o5cyZvChERFXmCrqt9BEmXFW0Jd4ucotcOMq9BRFR42jwMS2oSnqqwHtM3HkZ6RqZHWtG7devG1UWJiKjIE/StEqHmgE93adeWMAnHv3JERN5Usbl0Jq+PrpmLcCIxFct2Ze9P7morevny5R3KX3jhBU/UkIiIqMAJ+pcS2nr+sbSIPyqhyXqe9BiJR+ThWLP/ul6DiKjw6Jd+TXuj5Mk1aBh+Cl8s2VOgy2zcuNGh7OWXX0ZiYuKl1pCIiKhgCbr0I/9ENr9JREi8J3FEku8p5mJEmrDfJzHAfKxlU+SYwxLvm1MvzpBrfJrXaxAReUWTXvJ/Foyushrzth3Dmn2n3b6EzoterVo1h/InnnjCAxUkIiJyzpVJfW+WeF0iVSJO4iZzMSJN2DX5Hmc+HmbuK2ke+5pETxeuT0TkeXFVgAbd0eL4z6gUmYkP5u0s0GXWrFnjUDZu3DicOHHiUmtIRERUsARdFxmS0BVEq0sMlvhRYpuEzl2mCxClm4+1v/oPEo/qsXLOMC5QRERFqu1jsCSfwejq6zF/+zHsOZHk9iVKlSqFxo0bO5T36KHLPRAREXmey8viSbJ9VOJ/Er0kGkmUkggzQx83lrhD4iM91vNVJSJyU5XWQOXWuPrkdwix2PD1sr0FegsXL17sULZ8+XLs27ePt4SIiDyO61YTUWBrOwTBZ/fiP9V34vs/9yMpRb/0c09sbCxuuOEGh/JGjRp5ooZERETZMEEnosBW/0agZA30TpuCc5Kc/7T2YIEu88MP2oMvO53NxdlML0REREWeoFssliESayWSJE5LzJfQwaVEREXLGgRc+Siij6/FHeUO4uule7TLntuXiYyMxKhRoxzKmzZtitRUHRdPRERUCAm6JNktJU5J/C0Rlssxk2TzjkRTczrGWIlrJH6SfTq41G1yXleJ7RI7JYblcVwriQyJ2wvyOkRUTFzWFwiPw5CIWdhxLLFAUy6q4cOHOy3XudGJiIgKqwW9ozm1os5nnpJzpyTGulroHfpQ4piEzpv+tsRus2y0HNPAnQrJ8dLchQ8ltMNnQ4k7paxhLse9KjHTnesTUTEUGgU074fKxxagYuh5TF61v0CXCQkJwa5duxzKR44ciWPH9FcgERGR9xP0qyX0u+Cfc9n/uLnVqQx0FpdBEk/qY4m1EppE3+9mnVpL7JTr7JLQ7421hd5Zd5kh5pSP/KtIRPlr1geWzDQMrbwZv64/jBOJDm0OLqlR4//auxPAKMqzD+D/HCSBhIQjhIRw33IIhEPwDFY5PAAtKtZa7YeftVRrP9uKPTxKvdpqtRatWouKFRWqFK2WQyVUhSKESwQCyBkgHOEIIQe5vueZI0yS3SSbbJLZ7P/X7/lmd2Z2d3ZeNzzzzjPv9EDv3r2rzL/kkkvYCkRE1CgJek8rQV/toQc7XiYjreWzJJnOtpfJ43yZPKKrWeUuvkiWcHZvZVrznJ+dbN0E6UUf35uIglXiYKDjYEwo/hSFxSV4Me2bOr/VmjVrqszbvn07Nm7cWJ8tJCIiMoQb/9+7RIkcSbg93d3jQmuqCfoHHpZ/4kjyfaFJfWWVr+jSmveZsl1af+79jUJC7pSJBjp06IC0tDQfN4XcQkfLYPsFJje1XXLrMeiz82V8t8NuvL5STvW1yEKbyLpdKz937twqZS0ff/wxsrOzERrafAbIclP7ke/YfoGLbRfcakrQo607hXqiveewylGOVl4o8/IkQdY7jLb2cZu0x7yL43lniYOV1hkh8baVnGtP/lXyWO94+s9K26A18Rro169fWWpqqo+bQm6hCQLbLzC5qu0KhwNPv4X7E9PxxtEe2Ckn536W2q9Ob6Ujt0RGVr12/vbbb8err75a3y11DVe1H/mM7Re42HbBraZuHi1biZLkN8HDstFWz/baal4fIeHr+GN67riPfGYPCX39NIn3KyXePSS6a8hTHZx4RuXknIioikjpL0i5FTE738f3ehfgjf/urdONi1RERAR27NhRZf5rr72Gbdu2cecTEVGDJeh2QeV3PdSf21dErfD0Qlkn0Rp20ae7gkiirf9a3m2NzrJVYr7M+1re7y4NX96LiKiKS34KtIjGT0tfxan8s3ht5Z467yS9WHTUKL2uvaLzzjsPBQUF3PlERNQgCfo7ElpH8pAkx9dpj7b2bMvzuY7ecW8jvNgJ/GZfN0oSch3Wsa9EL4nHrHkvanhY93aJqrf4IyLyJFr6F1JnIvbAZ7ir+xG8tOIbnMorqvO++uQT+3Kbiu65RweaIiIi8n+C/oZEukSshCbBOjrLTonxVnnLbEmOj3l57TRrnc993ywiogY0/HY5v9cWM6KWIFdKXB77aEud3yomJgbp6fpnsqJXXnmFo7oQEZH/E3QdJcW6YdAyqyfdGZq8e7ytnvSy68gtk6ynnkZ4ISJq2hsXDf8+YvcswUMjSzB/bSY+3Xa4zm+XkpKCCy+0B7Y6Z+jQocjJyanPlhIRURCqcSww7SGX0B7z86y7hmr0skpLvF1dVSoxRWKirKM97kRE7jLmR0DrRNy270GcH1+GRz/ciqIS/dNVN0uXLvU4f+rUqXV+TyIiCk61HqxXEu0MrfW2YncN6+6R+FDC879YRERuqEW/4XWEnNyDp3tuwK6jZ/D2mv11f7voaGzYsKHK/GXLlnlN3omIiDxpPnfTICLyVdcLgM6j0PvQvzCyWxs8/+lO4y6jdTVkyBBjHPTKxo8fj4yMDLYPERHVChN0IgpuQ6Yh5OhW/DKlGFk5BZhfj1509dJLL3mc379/f9ajExFRrTBBJ6LgNvA6ILwlhu57FcO7tcULad/Uqxddb2C0d+9ej8smTJhQ5/clIqLgwQSdiIJbq3bARfci5OuFePD8HBw6VYAFazPr9ZZdu3bF8uXLq8xftWoV5syZU6/3JiKi5o8JOhHRRT8GWidhyNe/Q0qXWDy/fCcKiurei65SU1Px+OOPV5k/ffp0LFq0iPuciIiYoBMRVTsu+hWPIOTgOvyub4bRi/7qF3vqvcNmzpzpcf6UKVOQlpbGBiEiIo/Yg05EpAbfCHRKQZ9NT2Fi39Z47pMdSN97vF77JjQ0FKdOnfK4bOzYsVizZg33PRERMUEnIvJIkmlMeBI4fQhPJ/wbiXFRmP76WuQUFNVrh8XGxmL79u0el40aNQrr169ngxARUcV/kio8IyIK9nHRh92KVmtfwDs9F+NkXlG9h11Uffr0MS4Q9SQlJQWbNm2q92cQEVHzwQSdiMjp2j8BQ29BwqYXMamLWYteVFJa7300evRoLFmyxOsNjtLT09kORETEBJ2IqIrQMODyB43pz9t9hgMn8zH7051+2VHjxo3DwoULPS4bMWIEVqxYwQYhIiL2oBMRVRGbBAyYgi573sW0IW0xe/lObDmY45cdpSO4vPPOO16HZvzHP/7BBiEiCnIscSEi8uSCu4DCHDzc9SvERIbjqaUZfttPN954I+bNm+dx2Q033IBnnnmGbUJEFMSYoBMRedJ5BNBpGFqufwV3XdoTn247gi9312/YRaebb74ZCxYs8LjsvvvuQ3JystchGomIqHljgk5E5ElICDB6BnBsO/4n5gt0iovCQ4s2++WCUdvUqVOxePFij8sOHjyINm3aYPny5WwfIqIgwwSdiMibQVOBbhch8pMH8fiVHbAt6zSeWeZ5TPO6Gj9+PNatW+d1+eWXX46kpCScOHGC7UREFCSYoBMRef0LKX8ir30OKCpA6je/x00juuCFtG/w5uq9ft1nw4YNw+HDh70uz8rKQrt27XDPPfegqKh+N04iIiL3Y4JORFSd+N7A2F8AWz/AY/13YWy/Dnjwn5ux9Ossv+63hIQElJWV4YknnvC6zuzZsxEREYFf//rXTNSJiJoxJuhERDUZcw+QeD7CF9+P56/vicHJcfjp/I3YfzzP7/vugQcewNKlS6td57HHHjMSdb1DKUtfiIiaHyboREQ1CQsHJs8GzhxDq7SHMfs7KUAI8OO31/v1olHblVdeiezsbPTo0aPa9Xbu3GmUvoSEhBjjq+fl+f+AwR9ycnKMC2J1O+sa27f7t/afiMjNmKATEdVG0hDgonuB9X9HlxOr8fh1g7F+30m/XzRq08R7165dWLVqVa3WX7RoEaKjo8sT2jlz5qCkpKRBtk0vavUluY6Li8O7775br8/s169fvRL8WbNmsSyIiAIGE3Qiotq6bCbQXmrSP7gX154Xh2kjzYtGP9h4sMH24ejRo43a9IULF/r0uunTpyM8PLxeSa1Genp6lXnDhw9voG/bcB5++GGjLKi++4OJPhE1BiboRES11SIKmPRn4OReYPEDeOTaARjRrS1+umAj0vc27DCIWsJSl0Sd3JnoN1Z4OsBiBMY+YNuFNHkb1CfEsPr8rWGCTkTki24XAhffJ3UecxG1bCZevnU4kuKicOfctTiSU9Dg+9JO1Gtb+kJERE0itNkl6HLkMUEiQ2KnxAMelk+W2CSxQWKtxMVNsZ1EFKS+9RBw4T3AmlfQ7qu/4W+3jcTpwmI8+uHWRtsEu/RFY8mSJY32uXVlX+hpb3Nt4/Tp07j55pubevOJiBpVqAv/iIfJ5HmJiRIDJG6WeTp1+kRiiPzxHirT/5F4pXG3koiCmp6+vGIW0P8aYOmv0bt4B+66rBfel1r099ZlNvrmjBs3rkJSW98LMmsjOTkZR48erXWiXVpaagwL6auYmBjMmzfP58TeGSwLIqJA47oEXYyS2Cl/VHdJnJXHb0tMdq4g83Mlyqyn0TqrkbeRiIKd3mVUh16MSQAW/hAzLu6MMT3bG/XoizcfatJNu/766+uV0DpDLwj1ND8zMxPx8fFN+j19LQuqb+hIOUREtVSvMXjD6/PiBpIssd/xXLujLqi8kvSqXycTveWe/OuIqxtn04iIHFq2Ba59Dph3A6JWPoVXv/8r3PTyfzHz3a8wuHMbJLdpyd3VjEyaNMlI1ANJWlpawG0zmdh2gU3y1PX1er3bfrjyhW6QyXjZrjus57fKZJQ8v8fL+pfK5CFZfoWHZXfKRAMdOnQYPn/+/IbbcGpQubm5xqluCjzB0Hb9tj2HxKzl2DD0UWxvcR4eXpmPFlKsd+uASIxKdGM/SO0FQ/s1Z2y/wMW2C2xjx45Nl9x0RHNK0MfI5BHZrvHW81/oVJ4/Uc1rdstkpKxzrJqbXJRlZGT4e3OpEXsSUlNTG+nTyJ+Cou0KcoC/jpXpKeAH/8G2vBj8fMEmfHXgFCYOSsSsyYPQoXVkU29lnQRF+zVjbL/AxbYLbJKb1itBd2MN+hqJPvLFekhEyONpEu87V5D5vUOsQSZlkiITXS+70beUiEhFxQI3vQmczQPmfw/94yOxcMaFuH9CP3yy7QiufGaFcQEpERFRQCbocrRRLJO7JXTcMB2zbL7M+1oS8bs0rNW+LbFZh1m0Rny5yXHRKBFR40voD0x5AciUPoa0JxEeFooZqb3x0Y8vRo/4aPz4rfWYv8Z5eQ0REZFnriyOlFz7I5l8VGnei47Hv5OJBhGRewycAmz/DrDyOWDwVKDjQPROaI137hyD6a+vwcz3NuGbY7n4+bh+RgJPRETkCf+FICLyp3G/lZKXNsCrV0myvtSYFREeipduHY5pI7vgpRW7cI/0pheV1GsELiIiasaYoBMR+VN0PHDHMiCuizH8otak48u/olUY8MT15+PBawbg35uzMOPNdSgsLuG+JyIiJuhERA2uXU8zSR95B3BgnRTs/QyYOwkoKcb0i3tg1uSBWLblMCY++xnW7DnOBiEiogrYg05E1BBatASufhr4yVfmdO8XwPq5xqLvjemO174/EiVlZbjlldWYv3Y/Skt5nTsRETFBJyJqeDoi7IjpQNcxwPIngDPm7RrF/i9cAAAdwklEQVRS+yXgnzMuwvnJcbj/H5vwnVf+i4IilrwQERF70ImIGidJn/g7oPA08KbUpZ89Y8xuGx2B+T8Yg8evG4zVu4/jrr+n41ReEVuEiCjIscSFiKgxJA0Bps4BDm0AFtwu9ehmIh4aGoLvXNAVj00ZjM93HMP4Z/+DL3Z6vSkyEREFASboRESNpf9VwDXPADuWAh/cqzd1KF+kSfp7My5Eq8gwoy591gdb2JtORBSkXHmjIiKiZmu49J7nHAJWPAm0TgK+9WD5ovM7t8GH91yC3y3ehjlf7MZbX+7D9SnJmDAoESO7t0NUi7Am3HAiImosTNCJiBpb6gPAaUnSP3sKiE0yh2O0tIwIwyOTBuKmkV3wqiTpC9Iz8ebqfejVIRqzv5OC85Ji2V5ERM0cS1yIiJriotGr/wj0nQh8+DNg7avGGOlOmoj/fuoQpP/6CvzllhScyi/GNX/+HD9bsBGrd2WzzYiImjH2oBMRNYWwcPOi0TenAv/6CbDudeC77wGt2lVYrXVUC0wcnIQLerbHsx9vx7vSo/4PiVE92uGbI7kY0CkWF/aKl2iPQclxCAuV5F/knS3GwZMF6J0Qg9zCYsREnvtzfzLvLA7nFOLAyTzsP56PyPBQRMtyXbdXhxhEyHN9/a6jZ5B1phRlZWUoLi1D1qkCY/rl7mz0lPVahIUiO7cQ2WfOonv7aGObiIio/pigExE1lYhWwG3/Aja/CyyaAbx2jZm0J/Svsmq76AjMmjwIv7zqPPxhSQYWb87CGEnKdxzONWrWVdtWLTC2XwK6tm+FBWszcfBUPi7r2wFpGUdxSZ94dIprid3ZZ5C+9wRKvNwYKVwS/G7yek3cz5aUGvN+898lxjS/hnHah3SOQ5J8Rq+EaJzIK8LI7m0xWg4sEmOjkHH4NPYcO4OMrFzsPX7GOBCIlnKeVhHh6NMxBp3btjK+o32AQUQUzJigExE1pVCpNDz/BiC6PfCu1KK/nArc9AbQ50qPq+uFog9eM8AI29HThVj5zTEjEf8044j0kBcZJTL9Elvj021HcNXgRGzKPIVtWaeR3KYl7rikBwZLb3t8TKSRKBdJIn4qvwjbJYnOkHV2SM+8Jvop3dpizYbNCG2bbCT0+n7agz5C5meeyJdkWg8KIoz4aPMh4/O3ZuVg8ddZRo/9vNX7jO3TxPu49LLb1T3toyPx3roDVb6bLmsn79U+JsJYp5Uk8HqQoAcWk4Z2Qh/p4deeft2G5fK9dkvC36lNFAZ2ikNC60hjXHnt1SciCnRM0ImI3KDX5cAPV0nJy7fNmxn1uAS44IdAv4lm5lqNDpKcTh6abIQm0hpaplIq08OnC4xe7Zp0ksTd0wWorbIzkJp67mDAVnndGam9jVCa8IfJNm8+eArr9500Dg60Vz21bwLiW0dIMh1llNAUFpXidEExthzKwRHZzmO5Z82SGZ2eKURWTpGRcP97/yG8s3Z/jd9BxUaFo3t8NHrLgUdiXBSGy8HE5f0ToCcM7N55LdnRshw9sGCPPRG5ERN0IiK3aN0RuP1DYNXzwIa3gLdvBrpLoj5gMpB4vtSYHDfvRjrwOsk2W3h8C0047aRTb4JUm+Tc3+xebB02UsMTLW2R/Njo9daSnOoUSGlNmpwZ0AT+TGGxcbwySHrNteZ93/E848yAJtwnJI5Jgq8966t2ZRtnFrS3vaWcddDyHE3eu7RrZbxGDwxay3MdvnJgp1jjAEUPbL45misHC4VGfX8bKRkaNyDRKCXq2SEakeHnhrnMKSgyyozssxpRckAUKdOtcrChBxm6D0rkQCBavmfe2RJok4SHmQdNRaWlRklPfzkjESpf5lT+WazdcwJrpfRoU+ZJ4z31dXp2o6McZGhrJsmZgq6y7UNkf+rZAr02oVjeRz9HryE4LfvlmGx3mXUAokPs6+N9OSXGe+p+KC7RawlKy6ch8r+E2Ej07diaZx6IXIYJOhGRm0TFAWN/CVx6P7B+LrDsEWDPZxXX+WQWkHIbMOwW6TLu1DTb2Yg0AZ4wKMnjMr1YVcOTYunJf3ddJrYczJFkO8Lold+bnYdhXdsYF7V+c/QMVu/OxortR8tr8rUuvp2U2PRNaI1Dpwow619bjPl6RqJz25ZGCU6yTLWOX8t8PG9vqJEEa/KtJToRmphbF9oaibqUNdn1/c7a/4FSdjRtZFfjs4wzC3IG4itJrnXTFm8uqPIam76nl0sKTCu/qGYhjAQ/rmUL43ONsJJ+473lzS/t0wEdY6MQIwc0WmbUXg6qvjpwCvvlQEc/V7+bGcDZ4lLjAMg4iyKv1e8VFmq+nx4gxcrn6IGMHjTptQ56IKL7Xl+v+0sPvpxT/W4hlaae1rOnPCtCzQUTdCIit47yMuJ/gKGShOefAHatkK7pKPmrHWX2sC9/FEh7HOg7QZL170mNSgoQk1BjOUww0R7rmyThrYkm8lk5BUavsybhmgjadhw+ja3SQ79p/0kcknW0V36DPNa6+j9MHWKsX1hcIr38pUYvfRfpGdfSGmWPfmOfUdDn+t463SMHCprgal6tyeqg5FjjrII32vOuBwSbDpw0rjHQBF6TX+2NzyssMe5AmySfayar2jdu/qewbcsWDD1/sOwLTZZDramEddBwQN9TDgL0/TS5LpR9oVMN3bac/CL86ZMdPu137eHXawU0SdfkW/eBngXR0YQcN89tEHpWRM886PULemBn7gfPib3ObykHU4lylkn3XUc5m6D7yF5HS8f0YET3lXmgYR9wmPuYqCExQScicrPwSMk6EoEhN52b13c8kP2N9LC/IfEmkPGROT9GSmSShgIlZyWzOiBJvtSQtJe68EjpYd4mpTMlRUCbbsCo/5Xk//vma9Jflx7WP0v2KL2z7XsBY+4GelxaMdE/ewY4vtt8n4JTQHw/s8QmtHnc2VQTMC058aSPlH9oTBri+5kKTeJaSELsfG5Pe8RHG1Fb2pOtpUA1lQNVFnN8O1IHyH8XXqR0bYtra/humljnSgKvFxJrCY+W03SX7dCLc/Vag5DQc73dmrw6S4EqH2TkWtceaKmSliLpwY2dFGvybvfEa4HOud75c2U7Xp/L//RgQLdPD2TypaxIy5DMdSqur29vTMvM76ZlUb4cOOh31Gsc9IyC3bpGwq/PzP+rMM88GLDmWc8r/PdQvq45T99fz2QcP1qARYc3mO9X6fX6tH1MpJyRCMdhOdNjcx44OD+z4nPH9jm2wfnA+ZqQWr6fzvC2bmMJlY+LbGGeBbK/m7a9MZU4eDLfODNWWeUDLk9bXfmYrPI6/j5oY4JORBSINJm+Qspfxv7K7F3P3glkrgGOZpgjwyScJ13DhUDWJuB0lnmxaXSCuY6Ou757hZmUb1kEJI+QxL0LsF+WzZ0k9RytgW4XSq/8UPTftloS+HRJ0k87PjzE/OcuTl4TKet2uQD41kNVxnCn5kFH5NHQMwM6kk9d6UFGbJQc2MkJBu2d1ot53UB7+vW6gyNyhsRO+DXZPyzP9aCkvHbfquM/IwcZOw/nokAOLpRR728dJNiJvnEtgPy85N0qzrMOEOxrBOw3ODfP/AwtZzqdW4rMwuPnXl8+NdfXA4siWVdLh3TfmttxLhl1PrCfO5ef2y5rXqV1z7224veq/D7O17pdhBz46EhWFXLpStvu6avY39fbOpW/v71P64MJOhFRINOe7D5XmIG7al6/VJKKT38LrJ1j/iujte6XzTRLaoqkZ2njW8DhzVLbsVRiCdq1iAMGStLec6wk/FJzHRFtHgQoPSjQi1bXzZWDhOXA7dKTH5fckN+WyO+0BEmTNg03SUtLQ2pqqtflevZBS5H0gmE3KU/eHUl9Y/Wjl8iHFlolWk725+tZj8YaijVE/qzWBxN0IqJgomUpds+7co4GozXudumLkl7DlStWIHWsJOfV2Se97H//NvDXy4GRdwC9v2WW0ugNmLQbURP50mKzVEcjJtFM5LVUJlzKcIjIZ1pK5K2cqCmdK+Vq/M8OlVTcSMClMjDQMUEnIgpGXoZprEBLZWrzr2xXKXG57X3g44fNi1c1nCJizHp4HSbSKSTMvLBVh5CM7mCW7XQZZda8n9oPtO0hZTbDzF5/O8HnxXlEFASYoBMRUf0lp0iS/gFw+rA5LOSRLcCAKZJUJ5m16dpzX3wWOHPEXOfEbimV2SaJeCZwaKNZK7/h79V/hta66/CSejFsrLxvwgCg4CSQe1QS/HigZVuzx37vSvOAIEkS/xbuKlsgIqoNJuhEROTfmy0NnurlXxxJmuM6m9F5eNXlecfNi1hDw82LXI9sBQ5/bSbZ2qu+9m/AohneP1tf16KVFOfmnOu5734x0O0icyjKlp5vmkRE5DZM0ImIyB20p12HkLTpTZi0nt120b3S0/4VUJQHHNsOnNwntaaxZo+6ls/kHjbHjNe7r2qivn0xsH+1OdXymzZdgfi+UkrTR6Z9zMc61fIaHZpyVxrw1QJzqgcLuj1aVtPrcqD/NebrQ1uYF9TqxbFHtpkXzWrir58Xm2yeOdAhLvV1+r46X0NLinQ8QpboEFEtMEEnIqLAoMmtlq2orqNrXr//Veb04AZzrPhjO4Bsid2fmSPSOO/equU3Oi+qjXnzJy2h0ST9xB5g5Wzgiz/V/Hl6E6niqmMsV6AJvvb0a8KuCXwPOZjoJ9uZl61fUA4UCs1hM/UAQNfRMwdaf6/baIeeUTASf51KhFtTLespK7HG99MoQ8eszVI6dNAcalPPLJSPB1dWcbvt0AuFYzub1wPoe+owmhpOWlKkByH25+g1Avq5+t201KmZjI9P1JSYoBMRUfPWaagZjtFpkJNpJux20q7Jbc9UczjJyiPLaEKa+aXZQ19SbPa26w2kOg40E3LttS/MNWvqdWx4rZXPlYT4zDGzt19DX6eJdmmROdXnegZg03wg/bWKn6ej3OgBgq6jCXOYfNbJ/UDBZvNGUXogodtQC1IoBMhmGTSBNnrwy++UYybselDgVYh5pkHPFOgZAE3E9ZoBTcw9iZOzDG27mUm6ccYgzDyb0HWMeUZCD0r0Zlp6FoKIvHLlLyQkJES6L6DdFXoY/kpZWdmTlZbfIhN7hEn5q4gfyjryF4OIiKgWo9NouYqGs4TGmxhJKvtf3TC7VXvptc5ey3ns5FmHqNRtrI4m75qoF+VbBwEFZtJu3NbTTo5DsfrLNbhg9GgzMa7cE27TJF0PNDT0/Y7vMg8INHHXC3o1Idf3tnvML/oJ0HmEmXzbn6VTPZDREiHtqdebZBnrS0K//0vzrrc2LUsqL/+xev/LH7cyzxLowY8eBJWfDfAQ9nUG+r30zIf2+uuoQETNgOsSdEm+NSl/XuJKCeniwBqZ974k4FLYV263xGUy74QsmyiPX5aQLgsiIqIAor3KWubiK+2BDvNQflJJfispRWnXs/r30gMDO1HWkXD0YKGuzr+h6jzjbMFes/dfp7v/I49zHAcX+eZZCON5vlnu4yxB8oVej6AXIetZB03w7WghZwAS+psHB+UcZT56QKFnSPQASUuQjAgzDxh0H+u+qXL9gIchSGtax4fl8UfljMlW7YN0ML5PlLltRrmS9R20zClxkLmcmgXXJehilMROSb53WQn72zKZLFGeoMuylY71/yshv0YiIiJyHT2Y0N5tpTXqA6+rfn07oddecvvCWuuMQHnYia3W6BdK4p8nCf6xDEBr7nUoT024dZkm+9r7rwcEG+dV/7malCstQXIBSbcBOblSa8aZDA/3N/B6YXJtDjB8XddRPmU/L1/N07IQD8s8rFfdsirrwU/vbx+sWWeKbJ6u4/A0rxkm6HqfaDm3Vi6zht7x6RL/btAtIiIiosZP6H3R54rqlxu99JUu4rWTNE1s9WyGXZevBwclRWZPvl6wqz37TuUJWYWZNazj2/K1a9ZgxIgRjsWObSq/BsBKQvNPmvcS0OXVfaYv21/tutWtZ02dPfy1XeZxPW/LnI8b4v31AmgN6/qRCgcklZN8b/OaV4Lu6Zt5/K9DetfHWgn6xV6W3ykTDXTo0AFpaWn+2kZqZLm5uWy/AMW2C2xsv8DG9gtcuSEJSMuodPfdKuxEVEp3QiUVquHSBWpM85tdgq495l0cz7V85aCH5FvH2npFYqKUvOj4VFXIfK1N10C/fv3KUlNT/b+11Cj04IrtF5jYdoGN7RfY2H6Bi20X3Nx4rLVGoo8k4D0kdKyraRLvO1eQ+V1l8p7ErZKEb2+CbSQiIiIiahCu60GXhLtYEvC75eESCa3KnyPzvpZ5d1nLX5TJQxLtJV6Q+Tq7WOY7CrWIiIiIiAKT6xJ0Jcn2RzL5qNI8Tcztx3fIRIOIiIiIqFlxY4kLEREREVHQYoJOREREROQiTNCJiIiIiFyECToRERERkYswQSciIiIichEm6ERERERELsIEnYiIiIjIRZigExERERG5CBN0IiIiIiIXCSkrK2vqbWgUISEhp2WS0dTbQXUWL3GM+y8gse0CG9svsLH9AhfbLrD1kxy7dV1fHO7PLXG5DNlRI5p6I6jOB1hr2X6BiW0X2Nh+gY3tF7jYdoHffvV5PUtciIiIiIhchAk6EREREZGLBFOC/nJTbwDVC9svcLHtAhvbL7Cx/QIX2y6I2y9oLhIlIiIiIgoEwdSDTkRERETkeqFBciXtBIkMiZ0SDzT19lD1pI32SHwlscG+Clqm7SSWSeywpm25H91B2mKOxBGJzY55XttLHv/C+i3qb3J802w11dB+j0gcsH6DGlc5lrH9XELaoovEcomtEl9L3GvN5+8vcNuOv70AEBISEiXxpcRGq/1+4+/fXrMvcZGdECaT7RJXSmRKrJG4Wb73libdMKquzfbIZIS00THHvN/L5LjMe9I6yGorj2dyNzY9aY9LZZIrMVfaZFB17SWPB8jjtyRGSXSS+FiirywraaLND3pe2u8RnSfPn3LuILafu0h7JMkkSdppnTzW8ZbTJaZI3C7B319gtt2NEvztuZy0WYhMoqX9cuVhC3n8uYQeZF3vr99eMPSg687YKTthl8RZefy2xOQm3ibynbbZ69bj160/ZOQC8rv6j0yO17K9dP7b8ppCid3yeKf1GyV3tZ83bD93td0hTfCsx3ozvq0SyRL8/QVu23nD356LSJsp7dhQmqBrlPnztxcMCbr+B7/f8Tyzhh8BNT39j3ypHHGmS9xpzeuof9CMheY0ocm2jmrDW3vx9xg47pbf3yarBMY+Tcv2cylpo+4yGSaxWoK/v8BtO8XfXoBUaEhskIdHJJbJv3V+/e0FQ4KupyEqa951PYHvIvkPO0WmEyV+ZJ2Cp+aBv8fA8BeJXhJDJfQfmaet+Ww/F5K/kTEyeVfiJ/K3M6e6VT3M47+H7mo7/vYCRJmUp0jo38jOEqOkLY0SQX/99oIhQdejlC6O57ojDzbRtlAtyH/wRvvIVI9KF1qngQ5bNXt27Z4uI/fy1l78PQYA+e0dtv7xKZWnf3WcimX7uYxV/6oJ3pvSXu9Zs/n7C9C2428v8JSVlZ2USZrEBH/+9oIhQdeLQvvIjuohESGPp0m838TbRF5IG0VbF8wYj2UyTmKz1Wa3WavpdBF3oqt5ay+dP03aNlJ/k/K4j8SXTbB9VA37HxjLdRL2CC9sP/ddqPY3ia2SJPzRsYi/vwBtO/72AkNISEgHiTbW45YyuUJimz9/e+ENseFuIv/hF8vOuFseLpHQEV3myLyvm3izyLuOEgulzez/PudJey2W53qgNV+m02W6T+IG7kR3kDbRK9NTJeLlsfYSPCzxpKf20t+ezJsvD3UUpWKJH8k8juDivvZLlcdDrVOwOqrSD3Rdtp/rXCRxq4QxLK0175cS/P0FbtvdzN9eQNBOjNelrcKszu758vfxX/J8lb/+7Wv2wywSEREREQWSYChxISIiIiIKGEzQiYiIiIhchAk6EREREZGLMEEnIiIiInIRJuhERERERC7CBJ2IiIiIyEWYoBMRuUhISEiZFd2belv8Tb7T5xJ6b4reHpapu3VMaIm8+u4Hed0yiRKJwfXfciKixsUEnYjITxxJpa+ht4lu1uQ7TrJuzvJ2WVnZTg+r6E1a/iwxRFeXOGxFXW9k9aj1b9wTdXw9EVGTafZ3EiUiakSaUHrSTqKFRIHEKQ/LjzseZ1jTIj9uV1Mn55ooPy6hd8Z7zMtq91rT+ySeLavnXfTk5Svkcz+Th1fL9GJ5/nl93o+IqDHxTqJERA39h9bsIb9M4nVJFG8Pwu8/USYfSXwm3/9SD8sTHAc3rWWdXD997m0yeU3iPXnPb/vjPYmIGgNLXIiIqKHdYU3f9rK8pf3AX8m5ZaF11uJa6yCAiCggMEEnInIRbxdHyvNHrPmvWRdU/khivcQZiUMSr0t0dqzfx5qXKVEgsVnif2v47FCJW60LLI9KnJU4KPGOxAV1/D7tZXKthJasLKi0LFW/kzzc4+H7azzimN9D4i8S2yXyrQtJ9+rZCYlfSMRX/mxJ9nNkssQqL7qlLttPRNQUWINORBR43pK4SeKsVaueKPE9iUusRLqXxL8l2lg17xESAyVeluVtJHH9Q+U3lPmtZfKexBXWLE2cT0skSdwoMVXWuVdeO9vHbR1rJcjb5bVHKy07a5W2hEnEe6jjN3rT5XNTZKJlQrqNsL7zGYmuVmj50HqJxR4+/wuJyRLjJJ7xcduJiJoEe9CJiALLFImrJb5rJawaWtedJdFD4rdWKYleFNlLkmJN0jVetF4/y+rVrmyulZxvst4/Wl4bJ9O21ggrxRJ/ktfqSCy+sNdPr7xA3n+lhB5cjHTMS3TEU9bsp6zvuVoiReZHSOh2RVuvfdbLxbdqrTW90LpYlYjI9fjHiogosGjSfLckqG9KnNXRToSOVnK/tfwHVs/0dTJ/l6PU40cSOrxhlJWAl5PE9Qor8ddSk7Gy/kcS+dZrT0roUIUPWv9m/MLH7R1lTTXxr6vR1lR78LWn3CCP8yTWSvyfxCovr91oTWMlzqvHNhARNRom6EREgSVT4g0P8z92PP6DJKza411OnpfKZLn1dFCl1+poJ+o1Wc855KPTPGs6VhJ6LUmpLS2RUcd8eE1lOZXeyxcnHGOp1+X1RESNjgk6EVFg2WIl25UdcTze7OW1dn23loc4XWhN/0+S7yxP4SgVaSXhqUTGm3hHolxXOkSjmivb8qTEaAmta6+RNZ66Xf5S5UJSIiI3YoJORBRYDnmaKXloSU3rCHudysltkqN8pmM1AUeSXluR1lTLburq5xIrrTr0mRJazpIjSfqnEj+UKB+m0QsdalHVtB4RkSswQSciIvvfgsmS6OsN7GqK8mERa8EumdELVetEPi9bJhdLXCnxnMR6a2QaHSHmBYnNziEmPbDPGOj7EBG5HhN0IiKyS18GNMCuOOalrMYn1sWwH0vohaIpVrnKD6wDgJ7ehlCUxD3S0XNenzp4IqJGwwSdiIjsEVC+3QC7IsOa6hCQfiNJ+gmJl60hIGGNhe6JfcOnMse2EBG5GhN0IiJ6zdoFI6THWW945JUs97Un/Av7veuym627m1Z3Uz1jOEhHrXtl9hjrW61SGSIi12OCTkQU5CRxXWzdRVTNkYT4NxJJzqRcYrLEInn6Rx/fXm+YpIb5ODwjHOOX75TX/kpisP0eVuL+LXn4mLXekhoSdB0rnogoIDBBJyIipT3n/5TQBPghiYOSAJ+UOGXVeeuySXXYVTo84y7rrp+pddzV3SQetW52lC/blG2NCqNjv3e23v8+L6+9ypq+U8fPJiJqdEzQiYhIe9HPSFwnu+Iaqzf9gHVxZYR1B1K9UdFUiRm+7C5rHPI51tNpdbxJkW7TsxJfShy1hls8I7FG4lcSQ+Vj9AZOFUgir2U1va0EPq0On01E1CR0uKwm+WAiIgoOkih3kokOzXhaopP8u1PYSJ/7tNWz/kv5zCca4zOJiPyBPehERNSgJDk+KJOXJNpJfL+RknO96dJ0q8d9dmN8JhGRvzBBJyKixvBbiVyJmTWMyuIvP5bQJP1xOUDQnnsiooDRGH8kiYgoyEmSfMQawnGIdWGnL3cjrYsT1sWuf2ngzyEi8jvWoBMRERERuQhLXIiIiIiIXIQJOhERERGRizBBJyIiIiJyESboREREREQuwgSdiIiIiMhFmKATEREREbnI/wOGEN5esCVLRwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Define the column names\n", + "\n", + "def process_data2(directory):\n", + " columns2 = [\"t\", \"S0 pop\", \"S1 pop\", \"S2 pop\", \"S3 pop\", \"S4 pop\", \"S5 pop\", \"S6 pop\", \"S7 pop\", \"S8 pop\", \"poptot\"]\n", + " file_paths = glob.glob(os.path.join(directory, 'dpop*.dat'))\n", + " dfs = [] \n", + " count=0\n", + " for file_path in file_paths:\n", + " df = pd.read_csv(file_path, delim_whitespace=True, names=columns2, skiprows=1)\n", + " dfs.append(df)\n", + " count+=1\n", + " print(count)\n", + " combined_df = pd.concat(dfs)\n", + " average_df = combined_df.groupby('t').mean().reset_index()\n", + " return average_df, len(file_paths)\n", + "avg_df1, file_count1 = process_data2(data_directory1)\n", + "avg_df2, file_count2 = process_data2(data_directory2)\n", + "#avg_df3, file_count3 = process_data2(data_directory3)\n", + "# Plot the data\n", + "\n", + "# Load pop.dat with 1-second interval data\n", + "columns_pop = [\"t\", \"S0 pop\", \"S1 pop\", \"S2 pop\", \"S3 pop\", \"S4 pop\", \"S5 pop\", \"S6 pop\", \"S7 pop\", \"S8 pop\", \"poptot\"]\n", + "pop_df = pd.read_csv('/home/adurden/Ari/TAB/momentum/m9_10000_split/exact/pop.dat', delim_whitespace=True, names=columns_pop, index_col=False)\n", + "pop_df['t'] = pop_df['t'] / 100\n", + "print(pop_df.head())\n", + "#Plot the dat\n", + "\n", + "\n", + "\n", + "plt.figure(figsize=(12, 6))\n", + "plt.plot(avg_df1['t'], avg_df1['S0 pop'], label='new')\n", + "plt.plot(avg_df2['t'], avg_df2['S0 pop'], label='old')\n", + "#plt.plot(avg_df3['t'], avg_df3['S0 pop'], label='rescales all directions')\n", + "plt.scatter(pop_df['t'], pop_df['S0 pop'], label='exact', color='black', marker='o', s=5) # Plot as dots with smaller size\n", + "plt.xlabel('Time (fs)', fontsize=24)\n", + "plt.ylabel('S0 Population', fontsize=24)\n", + "\n", + "#plot x axis up to 300 fs\n", + "plt.xlim(0, 300)\n", + "\n", + "plt.title('S0 pop vs Time', fontsize=24)\n", + "plt.legend(fontsize=18)\n", + "plt.grid(True)\n", + "#plt.savefig('S0_pop_m9_1000_split.pdf', format='pdf')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "2022.0.2", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..38f3a8f --- /dev/null +++ b/analysis.py @@ -0,0 +1,48 @@ +import pandas as pd +import matplotlib.pyplot as plt +import glob +import os +# Specify the paths to the directories containing the ene*.dat files +data_directory1 = '/home/adurden/Ari/TAB/tab-model/data6' +data_directory2 = '/home/adurden/Ari/TAB/tab-model/data4/old' +# Define the column names + +def process_data2(directory): + columns2 = ["t", "S0 pop", "S1 pop", "S2 pop", "poptot"] + file_paths = glob.glob(os.path.join(directory, 'dpop*.dat')) + dfs = [] + count=0 + for file_path in file_paths: + df = pd.read_csv(file_path, delim_whitespace=True, names=columns2, skiprows=1) + dfs.append(df) + count+=1 + print(count) + combined_df = pd.concat(dfs) + average_df = combined_df.groupby('t').mean().reset_index() + return average_df, len(file_paths) +avg_df1, file_count1 = process_data2(data_directory1) +avg_df2, file_count2 = process_data2(data_directory2) +#avg_df3, file_count3 = process_data2(data_directory3) +# Plot the data + +columns_pop = ["t", "S0 pop", "S1 pop", "S2 pop", "poptot"] +pop_df = pd.read_csv('/home/adurden/Ari/TAB/tab-model/data4/pop.dat', delim_whitespace=True, names=columns_pop, index_col=False) +pop_df['t'] = pop_df['t'] / 100 + + +plt.figure(figsize=(12, 6)) +plt.plot(avg_df1['t'], avg_df1['S0 pop'], label='new') +plt.plot(avg_df2['t'], avg_df2['S0 pop'], label='old') +#plt.plot(avg_df3['t'], avg_df3['S0 pop'], label='rescales all directions') +plt.scatter(pop_df['t'], pop_df['S0 pop'], label='exact', color='black', marker='o', s=7) # Plot as dots with smaller size +plt.xlabel('Time (fs)', fontsize=24) +plt.ylabel('S0 Population', fontsize=24) + +#plot x axis up to 300 fs +plt.xlim(0, 300) + +plt.title('S0 pop vs Time', fontsize=24) +plt.legend(fontsize=18) +plt.grid(True) +plt.savefig('momentum_rescale.pdf', format='pdf') +#plt.show() \ No newline at end of file diff --git a/calcH.py b/calcH.py index 283ad38..74d094a 100644 --- a/calcH.py +++ b/calcH.py @@ -1,18 +1,17 @@ -def buildH(dimH, x1, x2, w1, w2, c, delta): +def buildH(dimH, x, w1, w2, c, delta): #Build the Hamiltonian for at position (x1,x2) import numpy as np #initialize Hamiltonian matrix H = np.zeros((dimH,dimH)) - - H[0][0] = -1.0*w1*x1 #diabatic potential + H[0][0] = -1.0*w1*x[0] #diabatic potential i = 1 while i < dimH: - H[i][i] = w2*x1 - (i-1)*delta #diabatic potential - H[i][0] = c*x2 #diabatic coupling + H[i][i] = w2*x[0] - (i-1)*delta #diabatic potential + H[i][0] = c*x[1] #diabatic coupling H[0][i] = H[i][0] i = i + 1 pass diff --git a/cgauss.py b/cgauss.py index 82d85cd..76155e0 100644 --- a/cgauss.py +++ b/cgauss.py @@ -1,4 +1,4 @@ -def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale): # Determines which coherent sub-block +def gcollapse(dimH,ndof,deltatn,aforce,poparray,dcp,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale): # Determines which coherent sub-block """of the electronic density matrix the WF collapses into""" # Standard library imports ===================================== @@ -26,7 +26,10 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol, while i < dimH: j = i + 1 while j < dimH: - invtau[i][j] = ((aforce1[i]-aforce1[j])**(2.0)/(8.0*dcp1)+(aforce2[i]-aforce2[j])**(2.0)/(8.0*dcp2))**(0.50) + invtausum = 0.0 + for k in range(ndof): + invtausum+= (aforce[k, i]-aforce[k, j])**2/(8.0*dcp[k]) + invtau[i][j] = (invtausum)**(0.50) invtau[j][i] = invtau[i][j] j = j + 1 i = i + 1 @@ -370,7 +373,7 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol, # constructing npop from the density matrix if (track == 0): - return poparray + return poparray, track pass k = 0 @@ -390,5 +393,5 @@ def gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol, # print A[track] - return npop + return npop, track diff --git a/collapseanalysis.ipynb b/collapseanalysis.ipynb new file mode 100644 index 0000000..b0914f1 --- /dev/null +++ b/collapseanalysis.ipynb @@ -0,0 +1,334 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Standard library imports =====================================\n", + "import numpy as np\n", + "import math\n", + "import sys\n", + "import random\n", + "from scipy.optimize import lsq_linear\n", + "\n", + "# Custom function imports ======================================\n", + "from rtelem import minelem\n", + "from sremove import rstates\n", + "\n", + "poparray = np.sqrt(1/3)*np.array([1,1,1])\n", + "# Rules ========================================================\n", + "# i, j, k, l, m, and n are all reserved for integer incrementing\n", + "print('poparray:',poparray)\n", + "\n", + "\n", + "#\tNot super efficient, could put sooner, but here is the rank reduction \n", + "#\tto working with only populated electronic states -> add a tolerance to pass\n", + "#\tfrom namd-main.py\n", + "\n", + "vstates = []\n", + "dimH = len(poparray)\n", + "zpop = 1.0e-6\n", + "i = 0\n", + "while i < dimH:\n", + " if (poparray[i] >= zpop):\n", + " vstates.append(i)\n", + " i = i + 1\n", + "pass\n", + "\n", + "rank = len(vstates)\n", + "print('rank:',rank)\n", + "#\tprint 'odotrho'\n", + "#\tprint odotrho\n", + "#\tprint 'w'\n", + "#\tprint w\n", + "\n", + "#\tsys.exit()\n", + "# constructing vectorized target density matrix\n", + "# This target is rho**(d) in the original python manuscript\n", + "# However it is organized into a column vector of the diagonal and\n", + "# lower triangular elements so it can be used as the target in a \n", + "# library linear least squares optimization\n", + "\n", + "eseg = np.zeros((rank,rank))\n", + "nblock = 0\n", + "\n", + "vtarget=[1, 0.05, 0.3, 1, 0.9, 1]\n", + "#print('vtarget_matrix:',vtarget_matrix)\n", + "#print('precollapseentropy:',precollapseentropy)\n", + "iter = 0 # Used to track what column is sent to minelem\n", + "bcore = [] # block coordinates for fastest decaying element\n", + "listbank = []\n", + "nzthresh = 1.0e-6\n", + "eseg = ([1, 0.05, 0.3], [0.05, 1, 0.9], [0.3, 0.9, 1])\n", + "dgscale = 1.0e6\n", + "p, q, leave = minelem(eseg,rank,nzthresh)\n", + "\n", + "bcore.append(p)\n", + "bcore.append(q)\n", + "\n", + "# print 'bcore'\n", + "# print bcore\n", + "\n", + "# print 'pre-algorithm check before individual element stuff is made'\n", + "# sys.exit()\n", + "\n", + "while leave == 'no':\n", + "\n", + " # Finding the largest possible coherent sub-block\n", + " i = 0\n", + " bstates = []\n", + " while i < rank: # starting by adding all possible states\n", + " bstates.append(i)\n", + " i = i + 1\n", + " pass\n", + "\n", + " # checking for zeroes in eseg such that blocks can be removed from the coherent block\n", + " icstates = [] # list of states with zero elemnts in rows or columns shared by the\n", + " # coherence indicated by bcore\n", + "\n", + " if (nblock > 0):\n", + " icstates = rstates(rank,eseg,bcore,nzthresh) # function identifying the icoherent states\n", + " pass\n", + "\n", + "# print 'icstates', icstates\n", + "\n", + " i = 0\n", + " if (nblock == 0):\n", + " pass\n", + " else:\n", + "# print('is this getting killed by poor loop navigation?')\n", + " while i < len(icstates):\n", + " j = -(i + 1)\n", + " temp = bstates.pop(icstates[j])\n", + " i = i + 1\n", + " pass\n", + " pass\n", + "\n", + " listbank.append(bstates[:])\n", + "\n", + "# print 'bstates for current block', bstates\n", + "\n", + " # projecting out the current block from eseg\n", + " iter = len(bstates)\n", + " i = 0\n", + " temp4 = eseg[bcore[0]][bcore[1]]\n", + " while i < iter:\n", + " j = 0\n", + " while j < iter:\n", + " eseg[bstates[i]][bstates[j]] = eseg[bstates[i]][bstates[j]] - temp4\n", + " j = j + 1\n", + " pass\n", + " i = i + 1\n", + " pass\n", + "\n", + "# print 'updated eseg'\n", + "# print eseg\n", + "\n", + " nblock = nblock + 1\n", + " if (nblock > 500):\n", + " print('nblock blew up')\n", + " print('Pnc matrix')\n", + "# print(Pnc)\n", + " print('current eseg')\n", + " print(eseg)\n", + " sys.exit()\n", + " pass\n", + "\n", + " bcore = [] # block coordinates for fastest decaying elements\n", + "\n", + " p, q, leave = minelem(eseg,rank,nzthresh)\n", + "\n", + " bcore.append(p)\n", + " bcore.append(q)\n", + "\n", + "pass\n", + "\n", + "i = 0\n", + "while i < rank:\n", + " if (abs(poparray[vstates[i]]*eseg[i][i]) > nzthresh):\n", + " listbank.append([])\n", + " listbank[-1].append(i)\n", + " pass\n", + " i = i + 1\n", + "pass\n", + "\n", + "\n", + "#\tprint 'vstates', vstates\n", + "#\tprint 'listbank'\n", + "#\tprint listbank\n", + "\n", + "\n", + "\n", + "\n", + "#\tprint 'vtarget'\n", + "#\tprint vtarget\n", + "\n", + "#\tsys.exit()\n", + "\n", + "# generating the vectorized set of coherent block density matrices\n", + "# for the TAB wave function collapse step\n", + "# These are determined using the current time step electronic populations,\n", + "# and therefore must be recomputed at each collapse step. (Unlike the list of \n", + "# states populated in each block [listbank])\n", + "\n", + "A = [] \t#stores the block density matrices in form of A[block-ID][element]\n", + "aindex = 0\n", + "\n", + "\n", + "i = 0\n", + "while i < len(listbank):\n", + " k = 0\n", + " A.append([])\n", + " while k < rank:\n", + " l = k\n", + " while l < rank:\n", + " if (k in listbank[i] and l in listbank[i]):\n", + " m = 0\n", + " popsum = 0.0\n", + " while m < len(listbank[i]):\n", + " popsum = popsum + poparray[vstates[listbank[i][m]]]\n", + " m = m + 1\n", + " pass\n", + " if (popsum <= nzthresh):\n", + " aelem = 0.0\n", + " else:\n", + " if (k == l):\n", + " aelem = dgscale*((poparray[vstates[k]]*poparray[vstates[l]])**(0.5))/popsum\n", + " else:\n", + " aelem = ((poparray[vstates[k]]*poparray[vstates[l]])**(0.5))/popsum\n", + " pass\n", + " pass\n", + " else:\n", + " aelem = 0.0\n", + " A[aindex].append(aelem)\n", + " l = l + 1\n", + " pass\n", + " k = k + 1\n", + " pass\n", + " i = i + 1\n", + " aindex = aindex + 1\n", + "pass\n", + "\n", + "#\tTest to see if there is linear dependence, and if not will do the real\n", + "#\tlinear least squares optimization\n", + "At = np.transpose(A)\n", + "\n", + "#\tehrv = []\n", + "#\tehrv.append([])\n", + "\n", + "#\ti = 0\n", + "#\twhile i < len(A[-1]):\n", + "#\t\tehrv[0].append(A[-1][i])\n", + "#\t\ti = i + 1\n", + "#\tpass\n", + "\n", + "#\tehrvt = np.transpose(ehrv)\n", + "#\tehrw = lsq_linear(ehrvt, vtarget, bounds=(0.0, 2.0), method='bvls')\n", + "\n", + "#\tif (abs(ehrw.x[0]-1.0) <= pehrptol):\n", + "#\t\ti = 0\n", + "#\t\treturn poparray\n", + "#\tpass\n", + "\n", + "\n", + "#\tprint 'A[dimH-1]'\n", + "#\tprint A[dimH-1]\n", + "\n", + "# scipy linear least squares optimization\n", + "#\tAt = np.transpose(A)\n", + "optw = lsq_linear(At, vtarget, bounds=(0.0, 2.0), method='bvls')\n", + "\n", + "# Error analysis of the linear least squares target wave function\n", + "\n", + "#Convergence analysis\n", + "# Sort the weights in descending order\n", + "sorted_weights = np.sort(optw.x)[::-1]\n", + "\n", + "# Initialize cumulative sum and term counter\n", + "cumulative_sum = 0.0\n", + "term_count = 0\n", + "\n", + "# Calculate how many terms it takes to reach or exceed 0.90 cumulative weight\n", + "for weight in sorted_weights:\n", + " cumulative_sum += weight\n", + "\n", + "\n", + "\n", + "# Collapsing the wave function\n", + "\n", + "i = 0\n", + "ptotal = 0.0\n", + "while i < len(optw.x):\n", + " ptotal = ptotal + optw.x[i]\n", + " i = i + 1\n", + "pass\n", + "\n", + "#\tprint 'optw'\n", + "#\tprint optw.x\n", + "\n", + "check2 = random.random()\n", + "check = check2*ptotal\n", + "\n", + "track = 0\n", + "psum = 0.0\n", + "i = 0\n", + "j = 0\n", + "while j < 2:\n", + " psum = psum + optw.x[i]\n", + " if (check <= psum):\n", + " j = 3\n", + " track = i\n", + " else:\n", + " i = i + 1\n", + " pass\n", + "pass\n", + "\n", + "# track is the collapsed into density matrix according to A\n", + "# constructing npop from the density matrix\n", + "print('poparray:',poparray)\n", + "if (track == 0):\n", + " print('track == 0')\n", + " return poparray, precollapseentropy.real, cumulative_sum\n", + "pass\n", + "k = 0\n", + "index2 = 0\n", + "while k < rank:\n", + " l = k\n", + " while l < rank:\n", + " if (k == l):\n", + " npop[vstates[k]] = A[track][index2]/dgscale\n", + " index2 = index2 + 1\n", + " l = l + 1\n", + " k = k + 1\n", + "pass\n", + "\n", + "#\tprint 'vectorized target'\n", + "#\tprint A[track]\n", + "return npop, precollapseentropy.real, cumulative_sum" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py27", + "language": "python", + "name": "py27" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/diffH.py b/diffH.py index 7c44929..cb675df 100644 --- a/diffH.py +++ b/diffH.py @@ -1,18 +1,16 @@ # Calculate dH/dR -def dHcalc(dimH,x1, x2, w1, w2, c, delta): +def dHcalc(dimH,ndof,x, w1, w2, c, delta): import numpy as np - - dH1 = np.zeros((dimH, dimH)) - dH2 = np.zeros((dimH, dimH)) + dH=np.zeros((ndof, dimH, dimH)) - dH1[0][0] = -1.0*w1 + dH[0, 0, 0] = -1.0*w1 i = 1 while i < dimH: - dH1[i][i] = w2 - dH2[i][0] = c - dH2[0][i] = c + dH[0, i, i] = w2 + dH[1, i, 0] = c + dH[1, 0, i] = c i = i + 1 pass - return dH1, dH2 + return dH diff --git a/efF.py b/efF.py index d5b3117..057ed31 100644 --- a/efF.py +++ b/efF.py @@ -1,23 +1,29 @@ -# Calculate Ehrenfest force -def calEff(dimH, ct, x1, x2, w1, w2, c, delta): - - import numpy as np +#Calculate Ehrenfest force +import numpy as np +import sys +from diffH import dHcalc - from diffH import dHcalc - - ctbra = np.transpose(np.conjugate(ct)) - - dH1, dH2 = dHcalc(dimH, x1, x2, w1, w2, c, delta) +def calEff(dimH, ndof, ct, x, w1, w2, c, delta): + # build bra for later + ctbra = ct.conj().T # shape (dimH,) - cnorm = np.dot(ctbra, ct) + # dH has shape (ndof, dimH, dimH) + dH = dHcalc(dimH, ndof, x, w1, w2, c, delta) - Efft1 = -np.dot(ctbra, np.dot(dH1, ct))/cnorm - Eff1 = Efft1.real + # normalization scalar + cnorm = (ctbra @ ct).real # same as np.dot(ctbra, ct) - Efft2 = -np.dot(ctbra, np.dot(dH2, ct))/cnorm - Eff2 = Efft2.real + # 1) first mat‐vec: each slice dH[k] dot ct → shape (ndof, dimH) + tmp = dH @ ct - return Eff1, Eff2 + # 2) then dot each of those rows with ctbra → shape (ndof,) + # and scale / take real part + Eff = - (tmp @ ctbra) / cnorm + if np.any(np.abs(Eff.imag) > 1e-10): + print(Eff) + raise ValueError("Ehrenfest force has non-zero imaginary part, which is unexpected.") + sys.exit(1) + return Eff.real diff --git a/exact.f90 b/exact.f90 new file mode 100644 index 0000000..1066ea0 --- /dev/null +++ b/exact.f90 @@ -0,0 +1,260 @@ +program exact + implicit none + ! This program solves the TDSE of a nucleus on PESs of 3 coupled states + + !--------------- grid size and domain -------------------------------- + integer, parameter :: ngrid1 = 1000, ngrid2 = 1000 + real*8, parameter :: x1min = -4.0d0, x1max = 8.0d0 + real*8, parameter :: x2min = -4.0d0, x2max = 8.0d0 + + !--------------- model parameters ------------------------------------ + real*8, parameter :: w1 = 0.25d0 + real*8, parameter :: w2 = 0.025d0 + real*8, parameter :: delta = 0.01d0 + real*8, parameter :: epsil = 0.00d0 + real*8, parameter :: c = 0.025d0 + real*8, parameter :: pmass = 1.845d3 + + !--------------- integration parameters ------------------------------ + real*8, parameter :: deltat = 0.01d0 + real*8, parameter :: tstepmax = 4.0d4 + + !--------------- initial conditions ---------------------------------- + real*8, parameter :: R1bar = -1.0d0, R2bar = 0.0d0 + real*8, parameter :: P1bar = 10.0d0, P2bar = 10.0d0 + + !--------------- other parameters ------------------------------------ + real*8, parameter :: pi = 3.14159265359d0 + real*8, parameter :: alpha1 = 6.0d0, alpha2 = 6.0d0 + + !--------------- wavefunction arrays --------------------------------- + real*8 :: wfr(3,ngrid1,ngrid2) + real*8 :: wfi(3,ngrid1,ngrid2) + real*8 :: dwfrdt(3,ngrid1,ngrid2) + real*8 :: dwfidt(3,ngrid1,ngrid2) + real*8 :: V(3,ngrid1,ngrid2) + + !--------------- diagnostics arrays ---------------------------------- + real*8 :: pop(3) + + !--------------- workspace variables --------------------------------- + real*8 :: deltaR1, deltaR2 + real*8 :: prefac1, prefac2 + real*8 :: gauss1, gauss2 + real*8 :: R1, R2, R1disp, R2disp + real*8 :: PbarRdisp, cosPbarRdisp, sinPbarRdisp + real*8 :: kedenom1, kedenom2, halfdeltat + real*8 :: norm + real*8 :: start, finish + integer*8 :: igrid1, igrid2, istate, n + + !--------------- precompute constants -------------------------------- + deltaR1 = (x1max - x1min) / dble(ngrid1) + deltaR2 = (x2max - x2min) / dble(ngrid2) + prefac1 = sqrt(sqrt(2.0d0 * alpha1 / pi)) + prefac2 = sqrt(sqrt(2.0d0 * alpha2 / pi)) + kedenom1 = 2.0d0 * deltaR1**2 * pmass + kedenom2 = 2.0d0 * deltaR2**2 * pmass + halfdeltat = 0.5d0 * deltat + + call cpu_time(start) + + !--------------- initialize wavefunction and potentials -------------- + do igrid1 = 1, ngrid1 + do igrid2 = 1, ngrid2 + + ! boundary conditions + if (igrid1==1 .or. igrid1==ngrid1 .or. igrid2==1 .or. igrid2==ngrid2) then + wfr(:,igrid1,igrid2) = 0.0d0 + wfi(:,igrid1,igrid2) = 0.0d0 + else + R1 = x1min + (dble(igrid1)-0.5d0)*deltaR1 + R2 = x2min + (dble(igrid2)-0.5d0)*deltaR2 + + R1disp = R1 - R1bar + R2disp = R2 - R2bar + gauss1 = exp(-alpha1 * R1disp**2) + gauss2 = exp(-alpha2 * R2disp**2) + PbarRdisp = P1bar*R1disp + P2bar*R2disp + cosPbarRdisp= cos(PbarRdisp) + sinPbarRdisp= sin(PbarRdisp) + + wfr(1,igrid1,igrid2) = prefac1*prefac2*gauss1*gauss2*cosPbarRdisp + wfi(1,igrid1,igrid2) = prefac1*prefac2*gauss1*gauss2*sinPbarRdisp + wfr(2,igrid1,igrid2)= 0.0d0 + wfi(2,igrid1,igrid2)= 0.0d0 + wfr(3,igrid1,igrid2)= 0.0d0 + wfi(3,igrid1,igrid2)= 0.0d0 + end if + + V(1,igrid1,igrid2) = -w1 * R1 + do istate = 2, 3 + V(istate,igrid1,igrid2) = w2 * R1 - (istate-1) * delta + end do + + end do + end do + + !--------------- write initial wavefunction -------------------------- + open(unit=21, file="wf_init.dat", status="replace") + do igrid1 = 1, ngrid1 + do igrid2 = 1, ngrid2 + R1 = x1min + (dble(igrid1)-0.5d0)*deltaR1 + R2 = x2min + (dble(igrid2)-0.5d0)*deltaR2 + write(21,*) R1, R2, wfr(1,igrid1,igrid2), wfi(1,igrid1,igrid2) + end do + end do + close(21) + + !--------------- compute initial norm & populations ------------------ + norm = 0.0d0 + pop = 0.0d0 + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + do istate = 1, 3 + norm = norm + wfr(istate,igrid1,igrid2)**2 + wfi(istate,igrid1,igrid2)**2 + pop(istate) = pop(istate) + wfr(istate,igrid1,igrid2)**2 + wfi(istate,igrid1,igrid2)**2 + end do + end do + end do + norm = norm * deltaR1 * deltaR2 + pop = pop * deltaR1 * deltaR2 + + open(unit=30, file="norm.dat", status="replace") + open(unit=31, file="pop.dat", status="replace") + write(30,*) 0_8, norm + write(31,*) 0_8, pop + + n = 1_8 + + !--------------- main propagation loop ------------------------------- + do while (n < tstepmax) + + ! compute dwfrdt + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + do istate = 1, 3 + dwfrdt(istate,igrid1,igrid2) = -(( & + wfi(istate,igrid1-1,igrid2) - 2.d0*wfi(istate,igrid1,igrid2) + wfi(istate,igrid1+1,igrid2) ) / kedenom1 + & + ( wfi(istate,igrid1,igrid2-1) - 2.d0*wfi(istate,igrid1,igrid2) + wfi(istate,igrid1,igrid2+1) ) / kedenom2 ) + & + V(istate,igrid1,igrid2) * wfi(istate,igrid1,igrid2) + end do + end do + end do + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + R2 = x2min + (dble(igrid2)-0.5d0)*deltaR2 + do istate = 2, 3 + dwfrdt(1,igrid1,igrid2) = dwfrdt(1,igrid1,igrid2) + c * R2 * wfi(istate,igrid1,igrid2) + dwfrdt(istate,igrid1,igrid2)= dwfrdt(istate,igrid1,igrid2) + c * R2 * wfi(1,igrid1,igrid2) + end do + end do + end do + + ! propagate wfr by dt/2 + wfr = wfr + halfdeltat * dwfrdt + + ! compute dwfidt + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + do istate = 1, 3 + dwfidt(istate,igrid1,igrid2) = (( & + wfr(istate,igrid1-1,igrid2) - 2.d0*wfr(istate,igrid1,igrid2) + wfr(istate,igrid1+1,igrid2) ) / kedenom1 + & + ( wfr(istate,igrid1,igrid2-1) - 2.d0*wfr(istate,igrid1,igrid2) + wfr(istate,igrid1,igrid2+1) ) / kedenom2 ) - & + V(istate,igrid1,igrid2) * wfr(istate,igrid1,igrid2) + end do + end do + end do + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + R2 = x2min + (dble(igrid2)-0.5d0)*deltaR2 + do istate = 2, 3 + dwfidt(1,igrid1,igrid2) = dwfidt(1,igrid1,igrid2) - c * R2 * wfr(istate,igrid1,igrid2) + dwfidt(istate,igrid1,igrid2)= dwfidt(istate,igrid1,igrid2) - c * R2 * wfr(1,igrid1,igrid2) + end do + end do + end do + + ! propagate wfi by dt + wfi = wfi + deltat * dwfidt + + ! compute dwfrdt again + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + do istate = 1, 3 + dwfrdt(istate,igrid1,igrid2) = -(( & + wfi(istate,igrid1-1,igrid2) - 2.d0*wfi(istate,igrid1,igrid2) + wfi(istate,igrid1+1,igrid2) ) / kedenom1 + & + ( wfi(istate,igrid1,igrid2-1) - 2.d0*wfi(istate,igrid1,igrid2) + wfi(istate,igrid1,igrid2+1) ) / kedenom2 ) + & + V(istate,igrid1,igrid2) * wfi(istate,igrid1,igrid2) + end do + end do + end do + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + R2 = x2min + (dble(igrid2)-0.5d0)*deltaR2 + do istate = 2, 3 + dwfrdt(1,igrid1,igrid2) = dwfrdt(1,igrid1,igrid2) + c * R2 * wfi(istate,igrid1,igrid2) + dwfrdt(istate,igrid1,igrid2)= dwfrdt(istate,igrid1,igrid2) + c * R2 * wfi(1,igrid1,igrid2) + end do + end do + end do + + ! propagate wfr by dt/2 again + wfr = wfr + halfdeltat * dwfrdt + + ! diagnostics every 100 steps + if (mod(n,100)==0) then + norm = 0.0d0 + pop = 0.0d0 + do igrid1 = 2, ngrid1-1 + do igrid2 = 2, ngrid2-1 + do istate = 1, 3 + norm = norm + wfr(istate,igrid1,igrid2)**2 + wfi(istate,igrid1,igrid2)**2 + pop(istate) = pop(istate) + wfr(istate,igrid1,igrid2)**2 + wfi(istate,igrid1,igrid2)**2 + end do + end do + end do + norm = norm * deltaR1 * deltaR2 + pop = pop * deltaR1 * deltaR2 + write(30,*) n, norm + write(31,*) n, pop + if (mod(n,400)==0) call dump_density(n, wfr, wfi, ngrid1, ngrid2) + end if + + n = n + 1_8 + end do + + call cpu_time(finish) + print*, "Program completes in ", finish-start + + contains + + subroutine dump_density(nstep, wfr, wfi, ng1, ng2) + implicit none + integer*8, intent(in) :: nstep + integer, intent(in) :: ng1, ng2 + integer :: i, j, k + real*8, intent(in) :: wfr(3,ng1,ng2), wfi(3,ng1,ng2) + real*8, allocatable :: rho(:,:) + character(len=20) :: fname + + allocate(rho(ng1,ng2)) + rho = 0.0d0 + do k = 1, 3 + do i = 1, ng1 + do j = 1, ng2 + rho(i,j) = rho(i,j) + wfr(k,i,j)**2 + wfi(k,i,j)**2 + end do + end do + end do + + write(fname,'("rho_t",I8.8,".bin")') nstep + open(unit=99, file=trim(fname), status="replace", form="unformatted", access="stream") + write(99) rho + close(99) + + deallocate(rho) + end subroutine dump_density + + end program exact + \ No newline at end of file diff --git a/mainout.py b/mainout.py index 43927af..16c4816 100644 --- a/mainout.py +++ b/mainout.py @@ -1,4 +1,4 @@ -def writemain (t,dimH,x1,x2,ct,odotx1,odotx2,H,posout,eneout,popout,dpopout,outp,pmass): #writing to output files +def writemain (t,dimH,ndof,x,ct,odotx,H,posout,velout,eneout,popout,dpopout,outp,pmass): #writing to output files import numpy as np import sys from hwrsort import eigsort @@ -27,7 +27,7 @@ def writemain (t,dimH,x1,x2,ct,odotx1,odotx2,H,posout,eneout,popout,dpopout,outp EMFr = EMF.real #Real part of the MF energy #Kinetic energy - KE = 0.5*pmass*(odotx1**2.0+odotx2**2.0) + KE = 0.5*pmass*(np.inner(odotx,odotx)) #Kinetic energy of the wave function Etot = EMFr + KE #Total Energy @@ -80,10 +80,15 @@ def writemain (t,dimH,x1,x2,ct,odotx1,odotx2,H,posout,eneout,popout,dpopout,outp # Formatting and writing the outputs line1 = format(t,'.4f').rjust(8) - line2 = format(x1,'.10f').rjust(20) - line3 = format(x2,'.10f').rjust(20) - lineout = line1 + line2 + line3 + '\n' - posout.write(lineout) + lineout = line1 + for linei in range(ndof): + lineout+= format(x[linei],'.10f').rjust(20) + posout.write(lineout + '\n') + + lineout = line1 + for linei in range(ndof): + lineout+= format(odotx[linei],'.10f').rjust(20) + velout.write(lineout + '\n') line2 = format(EMFr,'.10f').rjust(20) line2p = format(dPE,'.10f').rjust(20) diff --git a/namd-main.py b/namd-main.py index 8b736c2..e925483 100644 --- a/namd-main.py +++ b/namd-main.py @@ -55,11 +55,14 @@ # number of total states in a single trajectory dimH = 9 +#number of nuclear degrees of freedom +ndof = 2 + # (absolute) slope of diabatic1 potential along x1 direction w1 = 0.25 # slope of diabatic2-dimH potential along x1 direction -w2 = 0.025 +w2 = 0.25 # linear coupling constant c = 0.025 @@ -78,15 +81,14 @@ #Number of trajectories to be run in a calculation trajnum = 1000 - #Maximum number of nuclear time steps within a simulation tstepmax = 6000 -#Decoherence correction parameter on x1 direction -dcp1 = 6.0 -#Decoherence correction parameter on x2 direction -dcp2 = 6.0 +#Decoherence correction parameter in each direction +dcp = [6.0, 6.0] +rescale='old' +reverse=int(0) #reverse the velocity for frustrated hops; 0 for no, 1 for yes #Particle mass (Nuclear mass) pmass = 1845 @@ -110,17 +112,23 @@ # Initial population only on state-0 intpop = np.zeros((dimH)) intpop[0] = 1.000 - +x= np.zeros((ndof)) #Initial position vector +odotx= np.zeros((ndof)) #Initial velocity vector #loops over trajectories-k -k = 450 +k = 101 while k <= trajnum: #-----------Initial Conditions for trajectory-k ----------------- + #Initial conditions, buildH, and diffH are still defined manually along each dof np.random.seed(k) - x1 = np.random.normal(0.0,0.204)-1.0 #Initial paritcle position on x1-direction - x2 = np.random.normal(0.0,0.204) #Initial particle position on x2-direction - odotx1 = np.random.normal(10.0,2.451)/pmass #Initial particle velocity on x1-direction - odotx2 = np.random.normal(10.0,2.451)/pmass #Initial particle velocity on x2-direction - KE = 0.5*pmass*(odotx1**2.0+odotx2**2.0) #Initial kinetic energy + random.seed(k) #Set the random seed for reproducibility + x[0] = np.random.normal(0.0,0.204)-1.0 #Initial paritcle position on x1-direction + x[1] = np.random.normal(0.0,0.204) #Initial particle position on x2-direction + odotx[0] = np.random.normal(10.0,2.451)/pmass #Initial particle velocity on x1-direction + odotx[1] = np.random.normal(10.0,2.451)/pmass #Initial particle velocity on x2-direction + KE = 0.5*pmass*(np.inner(odotx, odotx)) #Initial kinetic energy + + + t = 0.0 #Initial simulation time #-----------Initilize the Time dependent wave function---------- @@ -135,7 +143,7 @@ ct = cr + 1j*ci #Calculating the Hamiltonian matrix at initial positions - H = buildH(dimH, x1, x2, w1, w2, c, delta) + H = buildH(dimH, x, w1, w2, c, delta) #------------Creating trajectory-k specific output files-------- @@ -148,10 +156,25 @@ #Heading for position output file of each trajectory line1 = str('t').rjust(8) - line2 = str('x').rjust(20) - line3 = str('y').rjust(20) - heading = line1 + line2 + line3 + '\n' - posout.write(heading) + heading = line1 + for linei in range(ndof): + linei = str('x' + str(linei+1)).rjust(20) + heading+=linei + posout.write(heading + '\n') + + #Opening trajectory specific velocity output file + line1 = str('vel') + line2 = str(k) + line3 = str('.dat') + line = line1 + line2 + line3 + velout = (open(line, 'w')) + + #Heading for velocity output file of each trajectory + line1 = str('t').rjust(8) + for linei in range(ndof): + linei = str('v' + str(linei+1)).rjust(20) + heading+=linei + velout.write(heading + '\n') #Opening trajactory specific energy and norm output file line1 = str('ene') @@ -218,13 +241,13 @@ dpopout.write(heading) #Writing out t = 0 outputs - null = writemain(t,dimH,x1,x2,ct,odotx1,odotx2,H,posout,eneout,popout,dpopout,outp,pmass) + null = writemain(t,dimH,ndof,x,ct,odotx,H,posout,velout,eneout,popout,dpopout,outp,pmass) #----------Begin the simulation--------------------------------- #Compute the Ehrenfest forces - mfdF1, mfdF2 = calEff(dimH,ct, x1, x2, w1, w2, c, delta) - - amp = np.zeros((dimH),dtype=np.complex) + mfdF = calEff(dimH, ndof, ct, x, w1, w2, c, delta) + #----changed upto here + amp = np.zeros((dimH),dtype=complex) poparray = np.zeros((dimH)) oldpop = np.zeros((dimH)) @@ -237,7 +260,7 @@ sw, sVR = eigsort(dimH,w,VR) tsVR = np.transpose(sVR) - temp1 = np.zeros((1),dtype=np.complex) + temp1 = np.zeros((1),dtype=complex) i = 0 while i < dimH: @@ -255,17 +278,13 @@ i = i + 1 pass - oldforce1 = np.zeros((dimH)) - oldforce2 = np.zeros((dimH)) - dH1,dH2 = dHcalc(dimH,x1,x2,w1,w2,c,delta) #diratives of diabatic Hamiltonian - i = 0 - while i < dimH: - force1 = -np.dot(tsVR[i,:],np.dot(dH1,sVR[:,i])) - force2 = -np.dot(tsVR[i,:],np.dot(dH2,sVR[:,i])) - oldforce1[i] = force1 - oldforce2[i] = force2 - i = i + 1 - pass + oldforce = np.zeros((ndof, dimH)) + dH = dHcalc(dimH,ndof,x,w1,w2,c,delta) #diratives of diabatic Hamiltonian + + for j in range(ndof): + for i in range(dimH): + oldforce[j,i]= -np.dot(tsVR[i,:],np.dot(dH[j,:,:],sVR[:,i])) + #Time steps count n = 1 @@ -283,14 +302,13 @@ ct = cr+1j*ci #Step positions forward in time - acel1 = mfdF1/pmass - acel2 = mfdF2/pmass + acel = mfdF/pmass + - x1 = movex(x1, odotx1, acel1, deltatn) - x2 = movex(x2, odotx2, acel2, deltatn) + x = movex(x, odotx, acel, deltatn) #Calculte H at the new position - H = buildH(dimH, x1, x2, w1, w2, c, delta) + H = buildH(dimH, x, w1, w2, c, delta) #Propagate WF through the other half time step using H(t+dt) i = 0 @@ -306,22 +324,20 @@ norm2ct = (cnorm.real)**(0.50) #Storing forces from the last step - mfdFprev1, mfdFprev2 = mfdF1, mfdF2 + mfdFprev = mfdF #Calculate Ehrenfest forces - mfdF1, mfdF2 = calEff(dimH, ct, x1, x2, w1, w2, c, delta) - + mfdF = calEff(dimH, ndof, ct, x, w1, w2, c, delta) #Step velocities forward in time - odotx1 = vcalc(odotx1, mfdF1, mfdFprev1, deltatn, pmass) - odotx2 = vcalc(odotx2, mfdF2, mfdFprev2, deltatn, pmass) + odotx = vcalc(odotx, mfdF, mfdFprev, deltatn, pmass) #-------------- TAB Starts from Here --------------------------------------- poparray = np.zeros((dimH)) #array holding state populations - ampdir = np.zeros((dimH),dtype=np.complex) #Stores amplitude directions for each state - amp = np.zeros((dimH),dtype=np.complex) #Stores amplitudes for each state - KE = 0.5*pmass*(odotx1**2.0+odotx2**2.0) + ampdir = np.zeros((dimH),dtype=complex) #Stores amplitude directions for each state + amp = np.zeros((dimH),dtype=complex) #Stores amplitudes for each state + KE = 0.5*pmass*(np.inner(odotx, odotx)) w, VR = np.linalg.eigh(H) sw, sVR = eigsort(dimH,w,VR) @@ -341,7 +357,7 @@ temp4 = np.dot(temp3,temp1) poparray[i] = temp4.real - if (poparray[i] == 0): + if (poparray[i] < nzthresh): ampdir[i] = 1.0 else: ampdir[i] = amp[i]/(poparray[i]**(0.5)) @@ -353,45 +369,32 @@ rEMF = EMF.real roldEMF = rEMF - dH1,dH2 = dHcalc(dimH,x1,x2,w1,w2,c,delta) #diratives of diabatic Hamiltonian - newforce1=np.zeros((dimH)) #Adiabatic State Force along x1 direction - newforce2=np.zeros((dimH)) #Adiabatic State Force along x2 direction - i = 0 - while i < dimH: - force1 = -np.dot(tsVR[i,:],np.dot(dH1,sVR[:,i])) - force2 = -np.dot(tsVR[i,:],np.dot(dH2,sVR[:,i])) - newforce1[i] = force1 - newforce2[i] = force2 - i = i + 1 - pass - - aforce1 = np.zeros((dimH)) - aforce2 = np.zeros((dimH)) + dH = dHcalc(dimH,ndof,x,w1,w2,c,delta) #diratives of diabatic Hamiltonian + newforce=np.zeros((ndof,dimH)) #Adiabatic State Force + + for j in range(ndof): + for i in range(dimH): + newforce[j,i]= -np.dot(tsVR[i,:],np.dot(dH[j,:,:],sVR[:,i])) + + aforce = np.zeros((ndof,dimH)) + aforce = (oldforce + newforce)/2.0 #Average force for the current time step odotrho = np.zeros((dimH)) i = 0 while i < dimH: - aforce1[i] = (oldforce1[i] + newforce1[i])/2.0 - aforce2[i] = (oldforce2[i] + newforce2[i])/2.0 odotrho[i] = (poparray[i] - oldpop[i])/deltatn i = i + 1 pass #----------- New Collapse Routine Goes Here --------------------- npop = np.zeros((dimH)) - npop = gcollapse(dimH,deltatn,aforce1,aforce2,poparray,dcp1,dcp2,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale) + npop, track = gcollapse(dimH,ndof,deltatn,aforce,poparray,dcp,nzthresh,errortol,npthresh,pehrptol,odotrho,tolodotrho,nta,dtw,zpop,dgscale) - oldforce1 = np.zeros((dimH)) - oldforce2 = np.zeros((dimH)) - i = 0 - while i < dimH: - oldforce1[i] = newforce1[i] - oldforce2[i] = newforce2[i] - i = i + 1 - pass + oldforce = np.zeros((ndof,dimH)) + oldforce = newforce poparray = npop - + namp = np.zeros((dimH),dtype=complex) nct = np.zeros((dimH,1),dtype=complex) @@ -400,25 +403,113 @@ namp[i] = ampdir[i]*(npop[i]**(0.5))*norm2ct i = i+1 pass - + #print ('namp',namp) + #print ('namp norm', np.dot(namp, np.transpose(np.conj(namp)))) + #namp_residue= np.zeros((dimH,1),dtype=complex) i = 0 while i < dimH: nct = nct + namp[i]*np.transpose([tsVR[i]]) i = i+1 pass - - i = 0 - while i < dimH: - cr[i] = nct[i][0].real - ci[i] = nct[i][0].imag - i = i + 1 - pass - - ct = cr+1j*ci - ccon = np.conjugate(ct) - ccont = np.transpose(ccon) - oldnorm = cnorm - cnorm = np.dot(ccont,ct) + EMF = (np.dot(np.conjugate(np.transpose(nct)),np.dot(H,nct))/np.dot(np.conjugate(np.transpose(nct)),nct)).item() + if (abs(EMF.imag) > nzthresh): + print('Warning: Mean-field energy has non-zero imaginary part, which is unexpected.') + print('EMF:', EMF) + sys.exit() + rEMF = EMF.real + a_f = np.dot((nct.flatten()),ct) + vrescale = np.ones(ndof)/np.sqrt(ndof) # Effective "NAC" vector for rescaling velocities. We normalize it later. + deltav = 0.0 + if(track==0): + vrescale = np.ones(ndof) + else: + print('old pop', oldpop) + print('npop', npop) + a_r = np.sqrt(np.maximum(0.0,1-abs(a_f)**2.0)) + print('a_f', a_f) + print('a_r', a_r) + if (abs(a_r) > nzthresh): + print('nct', nct) + print('ct', ct) + ct_residue = (ct - a_f*nct.flatten()) / a_r + print('ct_residue', ct_residue) + #print('namp',namp) + #print('amp_residue',amp_residue) + F_i = calEff(dimH, ndof, ct, x, w1, w2, c, delta) + print('F_i', F_i) + F_f = calEff(dimH, ndof, nct.flatten(), x, w1, w2, c, delta) + print('F_f', F_f) + F_r = calEff(dimH, ndof, ct_residue, x, w1, w2, c, delta) + vrescale = F_i - (abs(a_f)**2)*F_f - (abs(a_r)**2)*F_r + vrescale = vrescale / np.linalg.norm(vrescale) + print('vrescale', vrescale) + + + + #--Rescaling Kinetic Energy------- + deltaEMF = rEMF - roldEMF + nKE = KE - deltaEMF + if not (track==0): + if (nKE < 0.0): + print ('frustrated hops are needed') + odotx = ((-1)**reverse)*odotx #reverse the velocity or not + #restore the old wave function + elif(rescale=='old'): + scale = (nKE/KE)**0.50 + odotx = odotx*scale + + i = 0 + while i < dimH: + cr[i] = nct[i][0].real + ci[i] = nct[i][0].imag + i = i + 1 + pass + + ct = cr+1j*ci + ccon = np.conjugate(ct) + ccont = np.transpose(ccon) + oldnorm = cnorm + cnorm = np.dot(ccont,ct) + norm2ct = (cnorm.real)**(0.50) + + else: + tempv=np.dot(odotx, np.transpose(vrescale)) + print('tempv',tempv) + print('deltaKE',-1*deltaEMF) + discriminant = tempv**2 - (2*deltaEMF/pmass) + if discriminant >= 0: #we pick the root of smaller absolute value + if (tempv<0): + deltav = -1*tempv - np.sqrt(discriminant) + else: + deltav = -1*tempv + np.sqrt(discriminant) + + print('deltav',deltav) + print('odotx',odotx) + print('vrescale',vrescale) + oldKE=0.5*pmass*(np.inner(odotx, odotx)) + odotx = odotx + deltav*vrescale + print('odotx after rescale',odotx) + newKE = 0.5*pmass*(np.inner(odotx, odotx)) + print('KE sanity check:', newKE-oldKE+ deltaEMF) + i = 0 + while i < dimH: + cr[i] = nct[i][0].real + ci[i] = nct[i][0].imag + i = i + 1 + pass + + ct = cr+1j*ci + ccon = np.conjugate(ct) + ccont = np.transpose(ccon) + oldnorm = cnorm + cnorm = np.dot(ccont,ct) + norm2ct = (cnorm.real)**(0.50) + + else: + print("Warning: Negative discriminant encountered. Assuming frustrated hop..") + print("Discriminant:", discriminant) + odotx = ((-1)**reverse)*odotx #reverse the velocity, or not + #--Norm Conservation Check--------- # if (abs(cnorm-oldnorm).real >= 1e-12 or abs(cnorm-oldnorm).imag >= 1e-12 ): @@ -428,22 +519,10 @@ # sys.exit() # pass - mfdF1, mfdF2 = calEff(dimH, ct, x1, x2, w1, w2, c, delta) - EMF = np.dot(ccont,np.dot(H,ct))/cnorm - rEMF = EMF.real - - #--Rescaling Kinetic Energy------- - deltaKE = rEMF - roldEMF - nKE = KE - deltaKE - - if (nKE < 0.0): - print ('frustrated hops are needed') - print ('aborting program') - sys.exit() - pass + mfdF = calEff(dimH, ndof, ct, x, w1, w2, c, delta) - odotx2 = math.copysign(abs((2.0*nKE/pmass)-odotx1**2.0)**0.50,odotx2) + i = 0 while i < dimH: oldpop[i] = poparray[i] @@ -451,7 +530,7 @@ pass if (n%twrite == 0): - null = writemain(t,dimH,x1,x2,ct,odotx1,odotx2,H,posout,eneout,popout,dpopout,outp,pmass) + null = writemain(t,dimH,ndof,x,ct,odotx,H,posout,velout,eneout,popout,dpopout,outp,pmass) pass n = n+1 #forward one time step diff --git a/run.exe b/run.exe new file mode 100755 index 0000000..a75a60f Binary files /dev/null and b/run.exe differ diff --git a/test.slurm b/test.slurm new file mode 100644 index 0000000..4b842fa --- /dev/null +++ b/test.slurm @@ -0,0 +1,15 @@ +#!/bin/bash +# +#SBATCH --job-name=exa +#SBATCH --output=res.txt +#SBATCH --ntasks-per-node=1 +#SBATCH --nodes=1 +#SBATCH --time=48:00:00 +#SBATCH -p long-28core +#SBATCH --mail-user=fangchun.liang@stonybrook.edu +#SBATCH --mail-type=ALL + +cd $SLURM_SUBMIT_DIR + +gfortran exact.f -o run.exe +./run.exe diff --git a/traj_animate.py b/traj_animate.py new file mode 100644 index 0000000..1a2eb7c --- /dev/null +++ b/traj_animate.py @@ -0,0 +1,103 @@ +import os +import glob +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.animation import FuncAnimation + +# —————————————————————— +# 1) Paths to your data directories +new_dir = '/home/adurden/Ari/TAB/momentum/m9_10000_split/new' # new files (Δt = 0.05) +old_dir = '/home/adurden/Ari/TAB/momentum/m9_10000_split/old' # old files (Δt = 0.5) + +# —————————————————————— +# 2) Gather file lists +new_files = sorted(glob.glob(os.path.join(new_dir, 'pos*.dat'))) +old_files = sorted(glob.glob(os.path.join(old_dir, 'pos*.dat'))) + +# —————————————————————— +# 3) Load all data into numpy arrays +# Each array has shape (n_steps, 3) for columns [t, x1, x2] +new_data = [np.loadtxt(f, skiprows=1) for f in new_files] +old_data = [np.loadtxt(f, skiprows=1) for f in old_files] + +# Extract the time axes +times_new = new_data[0][:, 0] # 0.00, 0.05, 0.10, … +times_old = old_data[0][:, 0] # 0.00, 0.50, 1.00, … + +# —————————————————————— +# 4) Build index map: for each old time, find its index in new_times +# Assumes exact matching entries exist. +idx_map = [np.where(times_new == t)[0][0] for t in times_old] + +# Stack positions into arrays of shape (n_files, n_steps, 2) +new_coords = np.stack([d[:, 1:] for d in new_data], axis=0) # (n_new_files, n_new_steps, 2) +old_coords = np.stack([d[:, 1:] for d in old_data], axis=0) # (n_old_files, n_old_steps, 2) + +# Number of animation frames equals number of old time‐steps +n_frames = len(times_old) + +# —————————————————————— +# 5) Set up the matplotlib figure +fig, ax = plt.subplots(figsize=(6,6)) +ax.set_xlabel('x1') +ax.set_ylabel('x2') +ax.set_title('New (○) vs Old (×) Trajectories') + +# Create two scatter artists, one for each data set +scatter_new = ax.scatter([], [], marker='o', label='new') +scatter_old = ax.scatter([], [], marker='x', label='old') + +# Text handle for displaying the current time +time_text = ax.text( + 0.95, 0.95, + '', + transform=ax.transAxes, + ha='right', va='top' +) + +# Auto‐scale axes to fit all points +all_x = np.concatenate([new_coords[...,0].ravel(), old_coords[...,0].ravel()]) +all_y = np.concatenate([new_coords[...,1].ravel(), old_coords[...,1].ravel()]) +ax.set_xlim(all_x.min(), all_x.max()) +ax.set_ylim(all_y.min(), all_y.max()) +ax.legend(loc='lower right') + +# —————————————————————— +# 6) Animation functions +def init(): + scatter_new.set_offsets(np.empty((0,2))) + scatter_old.set_offsets(np.empty((0,2))) + time_text.set_text('') + return scatter_new, scatter_old, time_text + +def update(frame): + # Current time from old_times + t_cur = times_old[frame] + + # Old: direct lookup + pts_old = old_coords[:, frame, :] # shape (n_old_files, 2) + + # New: subsample at idx_map + j = idx_map[frame] + pts_new = new_coords[:, j, :] # shape (n_new_files, 2) + + # Update scatters and time text + scatter_old.set_offsets(pts_old) + scatter_new.set_offsets(pts_new) + time_text.set_text(f't = {t_cur:.2f}') + return scatter_new, scatter_old, time_text + +# —————————————————————— +# 7) Build and save the animation +ani = FuncAnimation( + fig, + update, + frames=range(n_frames), + init_func=init, + interval=100, # milliseconds between frames + blit=True +) + +# Save as MP4 (requires ffmpeg installed) +ani.save('trajectories_overlay_aligned.mp4', writer='ffmpeg', dpi=200) +print("Saved aligned overlay as 'trajectories_overlay_aligned.mp4'")