SnZn_q0_capture_h_-1_1_21.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "62564eab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from FastLine import Line\n",
    "import matplotlib.ticker as tck\n",
    "\n",
    "matplotlib.rcParams['axes.linewidth']    = 1\n",
    "matplotlib.rcParams['xtick.labelsize']   = 13\n",
    "matplotlib.rcParams['ytick.labelsize']   = 13\n",
    "matplotlib.rcParams['xtick.major.width'] = 1\n",
    "matplotlib.rcParams['xtick.minor.width'] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5039d53a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/6/eph/WSWQ'), (-0.0, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/5/eph/WSWQ'), (-0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/4/eph/WSWQ')]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAFgCAYAAAAo85E9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAArc0lEQVR4nO3deZhldX3n8fe3N6p6q2pkq65WFg24TgsyQZOQuEayqkzyJAyuWRqjzMQdNWSixKiImjFm0Ta4YtRE0MQkCojJiJO4sAjIKC4ISrWNIF290FX09p0/zrlVp27fqrpV3Ft1qur9ep77VN2z3d/1ObT96e/5fX+RmUiSJEmS6mnZfA9AkiRJkjQ5Q5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZKkJSkijoiISyPizojYHRFfj4hfKvc9OiKui4gd5evzEfHoyrkRERdHxE/K18URES0+4/kRkRHxe3P53SRJi4uhTZK0VK0Afgj8AtAHXAj8fUScAGwDfgM4EjgK+Cfg45VztwDPBjYD/wX4NeC86sUjYgPweuDWLn4HSdISYGiTJC1JmXl/Zr4hM+/IzEOZ+c/A94EnZOZwuT2BAA4Cj6ic/gLgHZl5V2YOAe8AXtj0EW8B/gK4t+tfRpK0qC360BYRvx0R10bErog40OXPOiYiPlQ+KrOrfNRmYzc/U5LUGRFxLHAylcpYRAwDo8C7gTdXDn8McFPl/U3ltsZ5Pw2cDryneyOWJC0VK+Z7AHNgB/DXQC+wtVsfEhE9wDXAl4FTgPuARwF7uvWZkqTOiIiVwEeBD2XmtxrbM7M/ItZQVNburJyyFthZeb8TWFvOa1tG8f8752fmoRZT3SRJmpFFX2nLzCsz82PA7a32R8SzI+L6iBiOiG9GxLmz/KgXAP3ASzLz3vJRm1szc9csrydJmgMRsQz4CLAPOL95f2beT1Ex+3BEHFNu3gOsrxy2HthTPk75EuDmzPxyVwcuSVoyFn1om0pEPAO4FHgZxWTzFwB/GRE/P4vLPQX4DvDB8vHIb0XEyzs2WElSx5WVsUuBY4H/lpn7Jzl0GbAaGCzf30rRhKRhM+OPVT4NeE5EbI+I7cDPAO+IiL/s9PglSUvDUng8cip/CLwrM68t3381Ii4Dng98sewg9v1prrEyMw9QdBd7CkUAfBFFN7HPRcSPM/Oj3Ri8JOlB+xuKR9mfnpkjjY3lP+rdC9wMrAHeRPG4/TfLQz4MvCIi/hVI4JUU896gaEjSU/mMK4BPUoRDSZJmbKmHthOBp0TEKyrblgONEPcD4OipLlAGNoDdwFBmvqt8f10ZAJ9FMU9CklQjEXE8RZv+B4Dtlbln51E8KvluYBMwAnwVOCszR8tj3gucBNxSvv/bchuZOdz0OfuAXZlZnQMnSVLblnpouxP4YGZe0mpnZh6i/VbNX6foFHbYZWY3NElSN2XmnRTt/CfzD1Ocm8Brytd0n/PkGQ9OkqSKRT+nLSKWl50dV5Xve8pXAP8beHlEnFketyoinhARrcLXdD4IPCQiXlpeazNwLsVjMZIkSZI0Kx0PbTNdFy0iTo+Ir0bE3oj4XkQ8t2n/MRFxRUTsjoh7IuListNXu55H8WjLlRSPPo6Ur+Mz8yrg94FLKCpqPwL+nKKV84yU/2L7y8DvAbso5i+8ITM/MdNrSZIkSVJDFE94dPCCEc+k6MTYC2zNzEkfwYyIPuC7wNspql4/D3wKeEZm/md5zNUU88VeBDwE+BxwaWZe3NGBS5IkSVINdTy0jV044snA56cJbS8C3gCcUM4PICI+AhzIzBdFxIkU66s9IjO/V+7/XeDCzDyxKwOXJEmSpBqZ70Ykm4Ebc2JyvIHikcbG/p2NwFbZf0JErG9euDoitgBbyrdP6NKYJUmSJKnjMrNlg6z5Dm3rgOYWyMPA+mn2Ux4zIbRl5lZgK0BEZLeqiJIkSZLUSZWlZw4z390jdwN9Tdv6GQ9jk+1v7JMkSZKkRW2+Q9tNwOObtp1abm/s74uIk5r23+EipZIkSZKWgm60/J9qXbRmnwLWRMSryzXSngacTfmIY2Z+H/g88LaIWF82JrkAeG+nxy1JkiRJddSNStuk66KVi1jviYiHAWTmMMXaZr9JMXftfcCLG+3+S+eW4xwCvgb8I/C2LoxbkiRJkmqnay3/55uNSCRJkiQtFBExaffI+Z7TJkmSJEmagqFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQa63hoi4jlEXFJRNwTEbsj4vKIOGqSY18fEXuaXhkRf1E55o6IGG065nGdHrckSZIk1VE3Km2vBZ4FnAFsKrd9pNWBmfnmzFzbeAGnAQlc1nTo71WPy8xbujBuSZIkSaqdFV245hbgosy8HSAiXgN8NyKOz8w72zj3xsz8ahfGJUmSJEkLTkcrbRHRDzwMuL6xLTO/B+wCNk9z7hHAC4H3ttj9zoi4LyK+HhHnTXGNLRFxXURcN4vhS5IkSVLtdLrStq78ubNp+zCwfppzfwNYBfxd0/YXUITAB4AnAx+PCDLzsHCXmVuBrQARkTMZuCRJkiTVUafntO0uf/Y1be+nqLZN5Tzgo5m5p7oxM/9PZu7JzP2ZeTXwTuC5nRisJEmSJNVdR0NbZg4DP6BoKAJARJxEUWW7ebLzIuLRwJnAe9r4mENAPKiBSpIkSdIC0Y3ukVuBCyLixIhYD1wMXJmZd0xxznnAlzPzpurGiDg+Ip4SET3lUgK/ALwc+EQXxi1JkiRJtdON0PZW4DPA14AhYDnl44wRcW5ETHj8MSJ6gefTusq2huJxyHuAHcBfUXSmfHcXxi1JkiRJtROZi7NfR0TkYv1ukiRJkhaXstliy2lg3ai0SZIkSZI6xNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNdTy0RcTyiLgkIu6JiN0RcXlEHDXJsU+OiIyIPZXXfzQd84iI+HxE3B8Rd0XEKzs9ZkmSJEmqq25U2l4LPAs4A9hUbvvIFMcfzMy1ldfPNHZExHLgM8A3gaOBXwcuiIjf6sK4JUmSJKl2uhHatgAXZ+btmbkTeA1wVkQcP4tr/TxwPPC6zNybmTcA7wVe3LnhSpIkSVJ9dTS0RUQ/8DDg+sa2zPwesAvYPMlpyyPihxGxPSL+JSKqx20Gvp2ZeyrbbpjsWhGxJSKui4jrHsz3kCRJkqS66HSlbV35c2fT9mFgfYvjvwU8HjgReCRwM/CFiNhYuV671yIzt2bm6Zl5+kwHLkmSJEl11OnQtrv82de0vZ+i2jZBZm7PzJsy80BmDmfm64D7gF+qXK+ta0mSJEnSYtTR0JaZw8APgNMa2yLiJIrK2M1tXuYQEOXvNwEnR8Sayv5Ty+2SJEmStOh1oxHJVooOjydGxHrgYuDKzLyj+cCIeGrZ0n9ZRKyNiDcAxwJXlod8EbgTeHNE9EbE44HzKJqRSJIkSdKi143Q9laKNv1fA4aA5cBzASLi3IioNhXZDFxD8Rjk7cATgWdk5g8BMvMg8GvAY4GfAP8KXJKZH+/CuCVJkiSpdiIz53sMXRERuVi/myRJkqTFJSLIzGi1rxuVNkmSJElShxjaJEmSJKnGDG2SJEmSVGOGNkmSJEmqMUObJEmSJNWYoU2SJEmSaszQJkmSJEk1ZmiTJEmSpBoztEmSJElSjRnaJEmSJKnGDG2SJEmSVGOGNkmSJEmqMUObJEmSJNWYoU2SJEmSaszQJkmSJEk1ZmiTJEmSpBoztEmSJElSjRnaJEmSJKnGDG2SJEmSVGOGNkmSJEmqMUObJEmSJNWYoU2SJEmSaszQJkmSJEk1ZmiTJEmSpBoztEmSJElSjRnaJEmSJKnGDG2SJEmSVGOGNkmSJEmqsY6HtohYHhGXRMQ9EbE7Ii6PiKMmOfaXI+ILEXFvROyIiGsj4symYzIi9kbEnsqrr9PjliRJkqQ66kal7bXAs4AzgE3lto9McuwG4N3AI4Cjgb8DPhsRD2067hczc23ltbML45YkSZKk2onM7OwFI+4ELsrMS8v3Dwe+C5yQmXe2cf524CWZeUX5PoEzM/NLMxxHdvq7SZIkSVI3RASZGa32dbTSFhH9wMOA6xvbMvN7wC5gcxvnPw44Criladc/lI9QfiUizp7i/C0RcV1EXDeb8UuSJElS3XT68ch15c/mxxeHgfVTnRgRxwCXA2/PzO9Udj0dOJHiUct3Ah+NiLNaXSMzt2bm6Zl5+izGLkmSJEm10+nQtrv82dwopJ+i2tZSRGwE/g24CnhddV9mXpOZo+XrE8BlwLkdG7EkSZIk1VhHQ1tmDgM/AE5rbIuIkyiqbDe3OiciTgCuBT6bmee3MRHtENDyWU9JkiRJWmy60T1yK3BBRJwYEeuBi4ErM/OO5gMj4pHAl4CPZearWux/bET8dESsioiVEfFs4HnA33dh3JIkSZJUO90IbW8FPgN8DRgClgPPBYiIcyNiT+XYC4BB4GVN67A1Hn88GvgAsAP4MXAh8DuZ+U9dGLckSZKkRer+Bw7wnbt38++3/Zgrb90+38OZkY63/K8LW/5LkiRJS8PBQ8k9ux9gaHiEbZXX0PBo+XOEnSP7x44/dv0RfOX1T5/HER9uqpb/K+Z6MJIkSZI0E/c/cGAsfG0rg9jY+50jbN85yv6DEws263tWsLG/l8H+Xp5w/AY29veysb+HTRt62djfO0/fZHastEmSJEmaN81VsqEWlbJqlQxg+bLguPU9DJZBrAhkRUAb3NDLQF8P63pWztM3mp2pKm2GNkmSJElds2dClazxGh17v33nKAcOta6SNapi46GsCGjHrOth+bLF1VDe0CZJkiSp48arZHvHqmLVKtnQjr3sGj0w4ZxqlWygv2esOtaolC3EKlknGNokSZIkzdiDqZKNh7LVE+aSLcYqWScY2iRJkiRNcPBQ8uPdoxO7LO6YOK9s0irZht6x+WQDfb1j75dqlawTDG2SJEnSEjNVlWxoxwh37zq8StbXu3LC3LGNZRCzStZ9hjZJkiRpEalWye7a0aINfosq2YplwXF9PWOPLja6Lg5WGn2sPcIVweaLoU2SJElaQHaP7udHO0cntr/fMV4p275rlIMt5pINblg9oUpWDWhWyerN0CZJkiTVxMFDyd27Rpl0sehJqmTHVuaSjVfHGl0YrZK149M3DnHJlbexbXiEjf29vPqZp/DsUwfne1iAoU2SJEmaM7tH948FsVbzyVpVyRpzyTb29Uxof9/4efS6I6ySPUifvnGI111xCyP7D45t6125nLec/bhaBDdDmyRJktQBBw4e4se7H5hQJRsa3jshpO1uUSUb6O9hY19v0xwyq2Rz6Wff+gWGhkcO2z7Y38v/fe1T52FEE00V2rw7JEmSpNJsqmT9q1eysa+XTRtWc8aJRzJglayWtrUIbFNtrxNDmyRJkpaE5irZUCWQtVMlO+PEI8ebe2wo2uIP9PWyxirZgrCxv7dlpW1jf+88jGZmvMMkSZK0KLSqklWDWTtVso0THl1cbZVsEXn1M09pOaft1c88ZR5H1R7ntEmSJKn2Dhw8xN1llaw5jLVTJWt0XRwLZeV6ZVbJlha7R9aMoU2SJGnh2DW6vxLIRitzySavkvX1rqzMHesZe2yxMZfsqLVWybRwGNokSZI0b1pVyYqFoscrZbsfmFglW7k8GOjrZaCv57AqmXPJtBgZ2iRJktQ1u0b3V0LY4ZWyu3c/cFiVbMPqYl2ygb5eNm0o5pBt7LdKpqXL0CZJkqRZqVbJhnYc3gZ/qirZWBDr621aMLqH1auskklVhjZJkiS11DyXrFox21auS9ZUJBurkjVCWOMRxo39vWwqq2TLrJJJM2JokyRJWoIOHDzE9l2jky4WPV2VbLB/9ViDjwGrZFJXGdokSZIWmcxk1+iBCVWx5rlk7VTJGuFsoL/HKpk0jwxtkiRJC8z+g4e4u6lK1lwp29OiStaYQ1Ztg19dMNoqmVRPhjZJkqQaaVUlu6vyyOK24RHunqJKNjihUtaollklkxYyQ5skSdIcalTJhnaMsG1nEcamq5KtWr6Mgf6eoqlHZS6ZVTJpaTC0SZIkdchsq2RHrllVtMDvm1glGyzXKDtqjVUyaSmb09AWEcuBtwIvBHqAq4DzMvPeSY4/C3gHcBLwPeAVmXlVZf8jgPcATwJ2AH+eme9oYxyGNkmSNGP7Dx5i+84ygM2wSja2HlmlBX5jjlnvquXz9I0kLQRzHdr+CHgBcBbwE+D9wOrM/KUWx54EfAPYAvw98JvAVuAxmXlHGQC/AXweuAB4JPA54H9k5iemGYehTZLUMRd++hY+9pUfcjCT5RGcc8ZDedOzHzffw9IMZSa7Rg6Mh7CdjeYeM6uSbdrQy0CfVTJJnTPXoe1O4KLMvLR8/3Dgu8AJmXln07FvBJ6amWdWtl0LfD4z3xgRTwH+BTgmM/eU+/8U+LnMfMo04zC0SZI64sJP38JlX/7BYduf+8SHGdxqplWV7K6mxaLv33dwwjmHVcnKrotjocwqmaQ5MFVo6+hs1ojoBx4GXN/Ylpnfi4hdwGbgzqZTNlePLd1Qbm/s/3YjsFX2v7SDw5YkaUof+8oPJ91uaJs7U1XJhnbsZdvwKHfvHqX532wbVbKTjl7Dz/3UUZWOi8VjjHZclFR3nW5BtK78ubNp+zCwfpLjWx37mGn2t7oWEbGF4lFLSZI65uAkT25Mtl2z06iSVeePNS8WPVmVbLC/l5/7qaMOX5vMKpmkRaDToW13+bOvaXs/sGuS46c6drr9E2TmVoo5cUSE/08qSeqI5REtA9rysDrTruYqWfPPyapkD1mzioEWVbLGz4esWWWVTNKi19HQlpnDEfED4DTg6zDWbGQ9cHOLU24CmuemnQpcU9l/ckSsycz7K/tv6uS4JUmayjlnPLTlnLZzznjoPIymnvYdKNclm2GVbGN/0WWxEciqi0UPWCWTJKB73SOfz3j3yEuBdZl5VotjHw7cAvwu8EngN4C/5fDukVcBrwVOAa4E/jAzPz7NOGxEIknqmKXcPTIz2Tmy/7Aui9Vq2Y93P9CyStYIYI3qWHU+mVUySRo3H+u0XUyxTtsRwNXAlsy8NyLOBd6bmWsrx1fXabsdeHmLddreS7FO2zDwzsx8exvjMLRJktSG1lWyiZWyvc1VshXL2NhXmTtWziUb7F89FtJ6Vlolk6R2zWloqwtDmyRJk1fJ7qoEtMmqZI12983VMqtkktR5hjZJkhapfQcO77hYtMKfuko22AhifRMbezQ6MVolk6S5ZWiTJGkBaq6SDe3Yy7amgNZulWywv1goeqCvl6PWriLsfClJtWJokySphmZTJTtixbJKZaynMp+srJT19Vglk6QFyNAmSdIca1TJ7tpRDWTjAW1oxwj37Dm8SnbU2lVji0IPbhhv8FHtuGiVTJIWH0ObJEkd9mDnkg30jbfAb4Qzq2SStHQZ2iRJmoHMZHjv/qZANsrQjvG1yVpXyY6YUBUb6Oth04Zeq2SSpGkZ2iRJqth34BA/2jky6WLR24ZHGdl/+Fyy8YWhew57hNEqmSTpwTC0SZKWjOYqWTWIzaRKVm32Mdjfy5FWySRJXWRokyQtGg8cOFiZSzaDKtmGMoQ1tcEfsEomSaoBQ5skaUHITHbs3d8Uwtqrko2tRdaiUmaVTJJUd4Y2SVItzLpK1t/LQBnKGmGs8ftxVskkSYuAoU2S1HWtqmRDOya2wb9n9wOHnXf0uiPG1iJrtMG3SiZJWmoMbZKkB+2BAwf5URm+JlTKdo6HtNH9hyac07Ny2XgA6xuvllklkyRpIkObJGlK1SrZXTtaLxY9XZWs0f5+oK93bG2yDatXWiWTJKkNhjZJWuJazSUbf3SxvSpZtePi4IaiSnbECqtkkiR1gqFNkhaxzOS++/dN6LA4kyrZpupi0WVIG+hzLpkkSXPJ0CZJC9jo/qJKNuO5ZOXjio0q2eCG8UqZVTJJkurF0CZJNdWqSlatlg0Nj3LvnsOrZMesO4KBSapkziWTJGnhMbRJ0jwZ3X+QH02oko2vR9aolrWqko0tEN2olvVbJZMkaTEztElSFzSqZEOVqti2Nqtk1XXIxgJaua3fKpkkSUuOoU2SZqF5LtlQU5VsaHiEBw5MrJL1rlw+9rhiI4wN9PUwuKGXTf2rObbvCKtkkiTpMIY2SWqSmfzk/n0TqmIzqpJtaLTC72GgDGhWySRJ0mwZ2iQtOdPNJWu3Sjb+6GKPc8kkSVLXGNokLSrNVbK7dkxs7LFteIR79+ybcE4EHL32iLGmHo0q2cbKfDI7LkqSpPliaJO0oEyoku2oVMp2joezyapkgxtWM9jfM7Y2WSOgHdfXw6oVy+bpG0mSJE3N0CapNhpVsqEd1TXJpq+SNeaSjS8Y3WPHRUmStGgY2iTNmea5ZI1w1qiSDQ2PsK9FlWxsLbK+8Rb4jWYfx663SiZJkhY3Q5ukjshM7t1T7bjYfpWs2tijaIG/emyNsr5eq2SSJGlpm7PQFhGrgb8Ezi43XQ6cn5kjkxz/fODFwKOAg8DXgNdk5i3l/hOA7wN7gcZAhzNzUxtjMbRJMzS6/+BhHRYbgWxoxwjbdo4eViVbvWp50+LQPQz0WSWTJEmaibkMbe8DHgM8hyJkfRq4KTP/YJLjXwp8G/gP4ADwv4DfAR6emXsroe2hmXnXDMdiaJMqWlXJmheL/sn9h1fJjl3Xc1gb/EYL/E39q1nfu8IqmSRJ0oM0J6EtInqB+4Bfzcxrym1PAz4DHJmZo21cowcYAZ6QmTcY2qT2taqSTVifbJIq2cQ1ySY297BKJkmSNDemCm0rOvg5pwA9wPWVbTcAvcDJwM1tXONpFI9Cfqdp+1ciYhVwK/CGzPz3VidHxBZgy8yGLdXfoUPJvfc/MD5/rM0qWWMu2WMH+3jmY46bUCVzLpkkSdLC0FZoi4gPAi+Y4pA/A64sf99Z2d74fX0bn3Ey8AHglZm5u9x8L/AkivC3kuLRyc9GxBmZeVgIzMytwNbyepbZtGBUq2RDw3sZqoSzqapkg/29DJShbLDpEcbj+npYudwqmSRJ0kLX1uOREbGWooo2mb0U1bQbgQ2ZOVyet4HikcnNrUJW5fqPBq4G3pWZb5tmLFcDX87MP57mOB+PVC3Mtkp27LoeBjc0Oi2WYaxv/NFF55JJkiQtHg/68cjM3APsmeZDbgNGgdOAL5SbT6WYo/btKc47Dfgc8KeZ+e42hnMI8G+qqo2RfQfH2t0XgWxiOPvR8Cj7Dk6skq1ZNb4u2eM29RVrk20YD2VWySRJktTQje6Rj2K85f8VwDcy88WTHP+zwD9TtPl/X4v9T6QIi9+iCJjPB/4C+NnMvL75+KZzrbTpQZusSla0vy8qZfc1VcmWBRy7vmfC/LGNfUV1rBHU1vdYJZMkSdK4uWz5vwZ4N5Os0xYRrwfOzczHlO//DfgFiscrq34pM6+NiHOAPwUGKKp4/w+4KDOvbmMshjZNq1WVbGjH+Npk01XJxueQlWuTWSWTJEnSLMxZaKsTQ5sOHUru3fNAOX/s8MWi262SNeaSDbgumSRJkrrE0KZFqVola1THxuaTzaBKNlgJaMeut0omSZKkuWdo04LzYKtkrRaLdi6ZJEmS6srQptoZ2Xew0vZ+YpVsaHiE7TunrpIN9k/8aZVMkiRJC5mhTXNquirZ0I4RduzdP+GcZQHHVeaSFfPHeq2SSZIkaUkwtKmj9u47MKEF/rbhEe6qLBb9o50j7D848X/7tUesGOuyuLGpSja4oZdj1x3BCqtkkiRJWqIMbWrboUPJPWNVsvEgVn0/XZWsMZdsoG98XbK+3pXz9I0kSZKk+jO0aUyjSjY0YS7ZzKtkmzb0joUyq2SSJEnSg2NoWyIeTJWs2ga/2nVxoM8qmSRJktRthrZFoqiSVdYia7NK1lggurnz4sZ+q2SSJElSHRjaFoDZVMmWLwuOW9/DQF9P04LRPZWOi1bJJEmSpLoztNXA3n0HGNoxMqEN/lilbGexLllzlWxdz4qxqlgjmA1WKmXHWCWTJEmSFgVDWw287oqb+dhXfzj2vlElqzb32NjfO7Y22UB/j1UySZIkaYkwtNXAdXfcx9DwyNicMqtkkiRJkhoMbZIkSZJUY1OFNks9kiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo2tmO8BLAWfvnGIS668jW3l4tqvfuYpPPvUwfkeliRJkqQFwMW1u+zTNw7xuituYWT/wbFtvSuX85azH2dwkyRJkgS4uPa8uuTK2yYENoCR/Qe55Mrb5mlEkiRJkhYSQ1uXbRsemdF2SZIkSaoytHXZxv7eGW2XJEmSpCpDW5e9+pmn0Lty+YRtvSuX8+pnnjJPI5IkSZK0kNg9sssazUbsHilJkiRpNjraPTIiVgN/CZxdbrocOD8zW07giogXAu8H9lY2fyYzz6kcczrw18BjgR8Bf5KZl7Uxllp0j5QkSZKk6cxl98h3AY8ETgFOBh4FvHOac27PzLWVVzWw9QGfpQh/G4AXA++JiCd1eNySJEmSVEsdC20R0Qs8F/jjzLw7M38M/DHwgojomeVlz6aowr0tMx/IzKuBTwFbOjJoSZIkSaq5TlbaTgF6gOsr224AeimqbpN5aERsj4gfRsTHI+LEyr7NwI1NzzneUG4/TERsiYjrIuK62X0FSZIkSaqXtkJbRHwwInKK15uAdeXhOyunNn5fP8mlvwg8DtgI/FdgFLg6ItaU+9c1XQ9geLLrZebWzDw9M09v53tJkiRJUt212z3yfOBVU+zfy3g1rY8iWDV+B9jV6qTMvL3ydntE/D5FSHsicA2wGzih6bT+ya4nSZIkSYtNW6EtM/cAe6Y6JiJuo6iUnQZ8odx8KjACfLvN8WT5anRNuQl4dtMxp5bbJUmSJGnR69ictrKt/2XARRFxTEQcA1wEfDgzR1udExG/EhGbonAk8FfAvcCXy0M+BayJiFdHxKqIeBpFc5KtnRq3JEmSJNVZp1v+v4yiqtZ43Qa8vLEzIl4fEbdWjn8y8FWKKt6twEOAZ5SVPTJzGPhl4DcpHpt8H/DizPzPDo9bkiRJkmqpo4tr14mLa0uSJElaKOZycW1JkiRJUgcZ2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqjFDmyRJkiTVmKFNkiRJkmrM0CZJkiRJNWZokyRJkqQaM7RJkiRJUo0Z2iRJkiSpxgxtkiRJklRjhjZJkiRJqrGOhraIWB0R74+I4fJ1aUT0TnH8eyJiT9MrI+IVlWMyIvY2HdPXyXFLkiRJUl1FZnbuYhHvAx4DPAdI4NPATZn5B22e/wzgX4HjM3NbuS2BMzPzSzMcS3byu0mSJElSt0QEmRkt93Uq2JQVtfuAX83Ma8ptTwM+AxyZmaNtXOOTwLLMPLuyzdAmSZIkaVGbKrR18vHIU4Ae4PrKthuAXuDk6U6OiOOAZwHvabH7HyLi3oj4SkSc3WJ/4xpbIuK6iLhuZkOXJEmSpHpqK7RFxAfLuWWTvd4ErCsP31k5tfH7+jY+5neBHwBXN21/OnAisAl4J/DRiDir1QUyc2tmnp6Zp7fzvSRJkiSp7tp6PDIi1lJU0Sazl6KadiOwITOHy/M2UDwyuTkzb57i+suA24G/ycyLpxnL+4CezHzeNMf5eKQkSZKkBWGqxyNXtHOBzNwD7JnmQ24DRoHTgC+Um08FRoBvT/MRZwEDwPvbGM4hoOWXkSRJkqTFpq3Q1o7MHImIy4CLIuIb5eaLgA+30YTkPOCKzLynujEiHgusBr5O0Y3yV4DnAb/dqXFLkiRJUp11enHtl1FU1Rqv24CXN3ZGxOsj4tbqCRExSBHGWjUgORr4ALAD+DFwIfA7mflPHR63JEmSJNVSR9dpqxPntEmSJElaKOaq5b8kSZIkqcMMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1ZihTZIkSZJqzNAmSZIkSTVmaJMkSZKkGjO0SZIkSVKNGdokSZIkqcYMbZIkSZJUY4Y2SZIkSaoxQ5skSZIk1VhHQ1tE/M+I+EpE7I2I77Z5zlkRcWtEjETENyLiF5v2PyIiPh8R90fEXRHxyk6OWZIkSZLqrNOVtm3A24A/a+fgiDgJuAJ4C9BX/vxURJxQ7l8OfAb4JnA08OvABRHxWx0etyRJkiTVUkdDW2Z+MjMvB4baPOUFwPWZeVlm7svMjwI3lNsBfh44HnhdZu7NzBuA9wIv7uS4JUmSJKmuVszz528Grm/adkO5vbH/25m5p2n/S1tdLCK2AFsq7zs3UkmSJEmaB22Ftoj4IOPVr1b+LDMvnMXnrwN2Nm0bBh4zzf71rS6WmVuBrbMYh2YhIjIzTcaaF95/mi/ee5pP3n+aL95786vdStv5wKum2L93lp+/m2IuW1U/sKvN/ZIkSZK0qLUV2srHE/dMe+DM3QQ8pWnbqcA1lf0nR8SazLy/sv+mLoxFkiRJkmqn0y3/V0RED7CyeBs95fvJfBg4PSLOiYiVEXEO8ATgQ+X+LwJ3Am+OiN6IeDxwHkUzEs2/N873ALSkef9pvnjvaT55/2m+eO/No8jMzl0s4g3AnzRvbzz/GhHnAu/NzLWVc84C3gGcBNwOvDwzr6rsfwRFSHsSxXy2d2bm2zs2aEmSJEmqsY6GNkmSJElSZ3V6cW1JkiRJUgcZ2iRJkiSpxgxtaltErI6I90fEcPm6NCJ6pzj++RHxHxGxIyLujYjPRsTj5nLMWhxmce9tLu+37RGREfFzczleLWwRsTwiLomIeyJid0RcHhFHTXH8WRFxa0SMRMQ3IuIX53K8Wjxmcu9FxGBE/GNE3Fn+OffcuR6vFpcZ3n+/HBFfKP9+tyMiro2IM+d6zEuJoU0z8S7gkcApwMnAo4B3TnH8OorGNJuAQeAG4KqIWN3lcWrxmem9tw+4AvjV7g9Ni9BrgWcBZ1D8+QXwkVYHRsRJFPfaWyjWFX0L8KmIOKH7w9Qi1Pa9BxwCrgL+O3BX94emJWAm998G4N3AI4Cjgb8DPhsRD+32IJcqG5GoLWVV4z7gVzPzmnLb04DPAEdm5mgb1+gBRoAnZOYN3RyvFo8He+9FRAJnZuaXuj5YLQoRcSdwUWZeWr5/OPBd4ITMvLPp2DcCT83MMyvbrgU+n5m2x9aMzOTeazrvDuDCzLxsTgaqRWm291/l/O3ASzLziu6OdGmy0qZ2nQL0ANdXtt0A9FJUPtrxNGAv8J3ODk2LXCfuPaktEdEPPIzK/ZaZ3wN2AZtbnLKZifcmFPdnq2OlSc3i3pM65sHef+X0l6OAW7o0xCXP0CYi4oPl8/CTvd5E8agjwM7KqY3f17fxGScDHwBemZm7O/sNtFDNxb0nzVCr+w2KdUJb3W/rZnCsNJWZ3ntSJ836/ouIY4DLgbdnpv8w3yUr5nsAqoXzgVdNsX8v4xWNPor/gBu/Q/GvMJOKiEcDV1P8x/ye2Q9Ti1BX7z1pFhr/qNTXtL2f1vfb7hkcK01lpvee1Emzuv8iYiPF3/GuAl7XlZEJsNImIDP3ZOa9U7z2ArcBo8BplVNPpZij9u3Jrh0RpwH/Drw1M9/Wxa+hBaib9540G5k5DPyAyv1WNhtZD9zc4pSbmHhvQnF/3tSlIWqRmsW9J3XMbO6/suHStcBnM/P8tFFGVxna1JbMHAEuAy6KiGPKUvhFwIcnawQRET8LXAP8UWa+e+5Gq8VklvdeRERP2fwGYFX5fvkcDVsL21bggog4MSLWAxcDV2bmHS2O/TBwekScExErI+Ic4AnAh+ZuuFpEZnLvUflzLoCV5XufotJstX3/RcQjgS8BH8vMqZ6YUYcY2jQTL6OobDRetwEvb+yMiNdHxK2V499EUWb/84jYU3m5jodm6mXM7N47nqISN1K+v6b8/XlzMVgteG+l6E76NWAIWA48FyAizo2IPY0Dy4n6ZwMXUjxCdCHwnMn+ki1No+17r9T4c+5hwPvL3y+cs9FqsZnJ/XcBxXJOL2v6O965cz3opcKW/5IkSZJUY1baJEmSJKnGDG2SJEmSVGOGNkmSJEmqMUObJEmSJNWYoU2SJEmSaszQJkmSJEk1ZmiTJEmSpBoztEmSJElSjf1/ZriRteMBU+gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use neutral structure and datapoints of disp_002,004,006\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from nonrad.ccd import get_Q_from_struct\n",
    "from nonrad.elphon import get_Wif_from_WSWQ\n",
    "from pathlib import Path\n",
    "from pymatgen.core.structure import Structure\n",
    "from glob import glob\n",
    "\n",
    "cc_dir = \"/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad\" + \\\n",
    "         \"/\" + \\\n",
    "         \"/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited\"\n",
    "\n",
    "SnZn_4 = Structure.from_file(f\"{cc_dir}/4/CONTCAR\")\n",
    "SnZn_5 = Structure.from_file(f\"{cc_dir}/5/CONTCAR\")\n",
    "SnZn_6 = Structure.from_file(f\"{cc_dir}/6/CONTCAR\")\n",
    "\n",
    "\n",
    "WSWQs = []\n",
    "for d in glob(f\"{cc_dir}/*\"):\n",
    "        if d[-2:] in [\"/4\", \"/5\", \"/6\"]:\n",
    "            pd = Path(d)\n",
    "            Q = get_Q_from_struct(SnZn_5, SnZn_6, f\"{pd}/CONTCAR\")\n",
    "            path_wswq = f\"{pd}/eph/WSWQ\"\n",
    "            WSWQs.append((Q, path_wswq))\n",
    "print(WSWQs)      \n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "ground_vr_path = f\"{cc_dir}/5/vasprun.xml\" \n",
    "\n",
    "Wifs = get_Wif_from_WSWQ(WSWQs, ground_vr_path, 2305, 2304, spin=1, fig=fig)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.ylim(-1e-6,1e-6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "15c2eba9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/6/eph/WSWQ'), (-0.0, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/5/eph/WSWQ'), (-0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/4/eph/WSWQ')]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAFgCAYAAAAo85E9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA170lEQVR4nO3de3hd+V3f+/dXN+uyt++WnPFcPB5LCoR0mGTOgR4IEFKa4daEtH1aTkJCbxPack4JEEIgbSGlkElCWko5TSZNGkIooSWXkl7IDVpCC4RkwuRywJZnxp5bbPnuvXWxZOnbP9bS1pZGtmXPlrQtv1/Psx9pr/VbWz8lazz+zO+3vt/ITCRJkiRJ7aljoycgSZIkSbo8Q5skSZIktTFDmyRJkiS1MUObJEmSJLUxQ5skSZIktTFDmyRJkiS1MUObJEmSJLUxQ5sk6aYUEVsi4j0RcSwiahHxpxHxneW5r42Iz0XE2fL1qYj42qZrIyIeiIjT5euBiIjy3O6I+J/l8XMR8YcR8U0b9XtKkm58hjZJ0s2qC3gC+FZgG/Am4D9ExH7gaeCvATuB3cBvAx9suvZ+4OXA3cBfAL4XeG15rg78bWAPsAN4APhYRHSt6W8jSdq0DG2SpJtSZk5k5s9k5tHMnM/M/ww8BrwwM8+VxxMIYA442HT5a4BfzMwnM/Mp4BeBHyw/dzozD2XmfNO1OygCoCRJ12zTh7aI+JsR8ZmIuBARl9b4Zw1GxK+WW2IulFttblnLnylJao2IGAJGgK80HTsHTAO/DPx80/DnAQ83vX+4PNb8eV8sr/1t4N9m5viaTFyStOndDFs1zgL/H9AHPLhWPyQieoFPA38EjAJngK+h2CYjSWpjEdEN/Drwq5n55wvHM3N7RAxQrKwda7qkApxven8eqERElKtzZOZfKP/d8H1Az1r/DpKkzWvTr7Rl5scz8zeAR1c6HxEvj4jPlw+L/1lEvPI6f9RrgO3AP8jMU+VWm69k5oXr/DxJ0jqIiA7g14AZ4IeXn8/MCeCdwPsjYrA8XAe2Ng3bCtQXAlvTtdPlv4N+MiLuXov5S5I2v00f2q4kIr4DeA/wIxTPGrwG+NcR8S3X8XEvBsaA95XbI/88Il7XsslKklqurPj4HmAI+KuZOXuZoR1AP7CvfP8ViiIkC+6maVvlCrqBA89utpKkm9VNHdqAfwT8UmZ+plwZ+yzwAeDVABGxPyLyKq+FLaa7KYLbZ4HnAK8CfvpZrNxJktbev6HYyv69mTm1cDAiviMi7omIzojYCryDYrv9n5VD3g/8aETsK59d/jHgfeW13xgR3xwRPRHRFxFvoAiFf7x+v5YkaTO5GZ5pu5I7gRdHxI82HesEPlN+/zhFyebLysyF4iY14KnM/KXy/eci4gPAyyiek5AktZGIuIOiTP9F4HjZZo3y2AxF8ZFbgSmK/yB3X2ZOl2PeRbFy9qXy/b8tjwFsAf5VeX62HPPdmfn0Wv4+kqTN62YPbceA92Xm21Y6WZZrPrXKz/pT4N6VPub6piZJWkuZeYyiJP/l/McrXJvAT5Sv5ef+B0u3TkqS9Kxs+u2R5daWXsrKXRHRW74C+JfA6yLiReW4noh4YUSsFL6u5n3Aroj4h+Vn3Q28Evhwi34VSZIkSTehloe2a+2LFhH3RsRnI2IyIh6JiFctOz8YER+OiFpEnIyIB8pKX6v1AxRbWz5OsfVxqnzdkZmfAP4e8DaKFbWvAv+CopTzNSn/i+13AX8XuAD8FvAzmfmb1/pZkiRJkrQgllUnfvYfGPFSikqMfcCDmXnZLZgRsQ04ArydYtXrW4CPAN+RmX9YjvkkxfNifwvYBfwO8J7MfKClE5ckSZKkNtTy0Nb44IhvAz51ldD2t4CfAfYv9LaJiF8DLmXm34qIOyn6qx3MzEfK838HeFNm3rkmE5ckSZKkNrLRhUjuBr6wrBnpQxRbGhfOn18IbE3n90fE1uWNqyPifuD+8u0L12jOkiRJktRymbligayNDm1V4PyyY+eArVc5TzlmSWjLzAeBBwEiItdqFVGSJEmSWqmp9cwzbHT1yBqwbdmx7SyGscudXzgnSZIkSZvaRoe2h4GvX3bsnvL4wvltEXFg2fmjmbl8BU6SJEmSNp21KPl/pb5oy30EGIiI15c90l4CvIJyi2NmPgZ8CnhrRGwtC5O8AXhXq+ctSZIkSe1oLVbaLtsXrWxiXY+I2wEy8xxFb7O/TvHs2ruBH1oo9196ZTnPp4A/Af4T8NY1mLckSZIktZ01K/m/0SxEIkmSJOlGERGXrR650c+0SZIkSZKuwNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltzNAmSZIkSW3M0CZJkiRJbczQJkmSJEltrOWhLSI6I+JtEXEyImoR8aGI2H2ZsT8VEfVlr4yIf9U05mhETC8b8/xWz1uSJEmS2tFarLT9JPAy4BuAW8tjv7bSwMz8+cysLLyAFwAJfGDZ0L/bPC4zv7QG85YkSZKkttO1Bp95P/DmzHwUICJ+AjgSEXdk5rFVXPuFzPzsGsxLkiRJkm44LV1pi4jtwO3A5xeOZeYjwAXg7qtcuwX4QeBdK5x+R0SciYg/jYjXXuEz7o+Iz0XE565j+pIkSZLUdlq90lYtv55fdvwcsPUq1/41oAf498uOv4YiBF4Evg34YESQmc8Id5n5IPAgQETktUxckiRJktpRq59pq5Vfty07vp1ite1KXgv8embWmw9m5v/IzHpmzmbmJ4F3AK9qxWQlSZIkqd21NLRl5jngcYqCIgBExAGKVbYvXu66iPha4EXAO1fxY+aBeFYTlSRJkqQbxFpUj3wQeENE3BkRW4EHgI9n5tErXPNa4I8y8+HmgxFxR0S8OCJ6y1YC3wq8DvjNNZi3JEmSJLWdtQhtbwE+BvwJ8BTQSbmdMSJeGRFLtj9GRB/walZeZRug2A55EjgL/ApFZcpfXoN5S5IkSVLbiczNWa8jInKz/m6SJEmSNpey2OKKj4GtxUqbJEmSJKlFDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGWh7aIqIzIt4WEScjohYRH4qI3ZcZ+20RkRFRb3r9r2VjDkbEpyJiIiKejIgfa/WcJUmSJKldrcVK208CLwO+Abi1PPZrVxg/l5mVptf/tXAiIjqBjwF/BuwB/grwhoj4G2swb0mSJElqO2sR2u4HHsjMRzPzPPATwH0Rccd1fNa3AHcAb8zMycx8CHgX8EOtm64kSZIkta+WhraI2A7cDnx+4VhmPgJcAO6+zGWdEfFERByPiP8SEc3j7gYOZ2a96dhDl/usiLg/Ij4XEZ97Nr+HJEmSJLWLVq+0Vcuv55cdPwdsXWH8nwNfD9wJPBf4IvC7EXFL0+et9rPIzAcz897MvPdaJy5JkiRJ7ajVoa1Wft227Ph2itW2JTLzeGY+nJmXMvNcZr4ROAN8Z9PnreqzJEmSJGkzamloy8xzwOPACxaORcQBipWxL67yY+aBKL9/GBiJiIGm8/eUxyVJkiRp01uLQiQPUlR4vDMitgIPAB/PzKPLB0bEt5cl/TsiohIRPwMMAR8vh/w+cAz4+Yjoi4ivB15LUYxEkiRJkja9tQhtb6Eo0/8nwFNAJ/AqgIh4ZUQ0FxW5G/g0xTbIR4FvBL4jM58AyMw54HuBrwNOA/8VeFtmfnAN5i1JkiRJbScyc6PnsCYiIjfr7yZJkiRpc4kIMjNWOrcWK22SJEmSpBYxtEmSJElSGzO0SZIkSVIbM7RJkiRJUhsztEmSJElSGzO0SZIkSVIbM7RJkiRJUhsztEmSJElSG+va6AlIkiRJ0lo6PznL4fEah47XGDtR46vnp3nw1fdu9LRWLTJzo+ewJiIiN+vvJkmSJOmZatOzjI3XGTtR49DxOmNlUBuvXWyMCSCB52zr5Q33PZeX37Nvw+bbLCLIzFjx3GYNNoY2SZIkaXOampnjyHidwydqTa86T52baozp7e5geLDK8FCF0aEqp+sz/OofHuXipfnGmL7uTn7hFc9vi+BmaJMkSZJ0w7l4aY5HT040gtnC6tnjZyZZ+Kt+T2cHB/YMMLq3yshQleHBCiNDVW7b2U9nx2IG+qa3/O6SULdg3/Y+/udPfvt6/UqXdaXQ5jNtkiRJkjbU7Nw8R09NcPjE0tWzo6cnmZsv0llnR3Dn7gG+7pZtfN89+xgZKkLa/l39dHVevb7i0ysEtisdbyeGNkmSJEnrYm4+efzMZKMgyKETNcZO1Hn0VJ3ZuSKcRcD+XQMMD1b4ruc/h+GhKqNDVe7cPUBP1/UXv79le9+KK223bO+77s9cL4Y2SZIkSS01P588dW6q8azZwsrZkfH6kmfKbt3Rx8hQlRc/d5CRoWJb48HBCr3dnS2f0+tfOsobP/wlpmbnGsf6ujt5/UtHW/6zWs1n2iRJkiRdl8zk+IXpIpgdX9zWODZeZ3JmMRw9Z1tvuWJWYXho8dmzgS3ru4b00S88xds+foinz01xy/Y+Xv/S0bYoQgIWIpEkSZL0LGQmp+ozz6jWePhEjdr0pca43ZUtjO6tMDxYLQuDVDg4WGVbX/cGzv7GYGiTJEmStCpnJ2aeEcwOn6hxdnK2MWZ7f3dZCKTSKAgyMlRl50DPBs78xmZokyRJkrTEhelZxspgduh4rWxEXedUfbERdXVLV9HnbO9iMBseqrCnsoWIFfOFrpOhTZIkSbpJTc5cYmxJKf3i+6+en26M6evuZKR83my0DGaje6vs3dprOFsnhjZJkiRpk5ueneORk3XGTtTLUvpFSf0nziyWue/p6uDgnkqxrXFvlZHy2bN92/vo6DCcbSRDmyRJkrRJzFya57FTE0WVxqZeZ0dPT1D2oaarIziwZ6DpebPi2bPbd66uEbXWn6FNkiRJusFcmpvn2JnJspR+ncPjNQ4fr/HYqQkulemsI2D/7gFGBquN1bPRoSr7dw/QbTi7oRjaJEmSpDY1P588eXaKQ8vK6T9yss5M2Yg6Am7b0f+Mao0H9gysSSNqrT9DmyRJkrTBMpOnz08Xoez4YkGQI+N1pmYXG1Hv295XFAIZqpaNqCscHKzQ37O+jai1vgxtkiRJ0jrJTE7WLpYrZ/UioI0Xz53VLy42oh6sbmF0b7VsRF1UbhwerFDttRH1zcjQJkmSJK2B0/WLSxpQj5XPnp1rakS9c6Cnsa1xoaT+yFCF7f02otYiQ5skSZL0LJyfnC0KgTRtbRwbr3GqPtMYs7W3q1g5G6oyMliW1B+qsruyZQNnrhvFuoa2iOgE3gL8INALfAJ4bWaeWmHsdwE/DvwFoBP4MvBTmfmZpjEJTAHzTZfuy8zzV5mHoU2SJEnXpH7xEmPlillzYZATFy42xgz0dDaeNVsoCDK6t8pgdYuNqHXd1ju0/TTwGuA+4DTwXqA/M79zhbGvBCaB3wPqwN8D3gZ8TWY+UY5J4EWZ+QfXOA9DmyRJklY0PTvHkfFiW+NCn7NDx2s8dW6xEfWWrg6GhypFOf29iyHtlm02olbrrXdoOwa8OTPfU76/CzgC7M/MY6u4/jjwDzLzw+V7Q5skSZKuy8VLczx6cqLxvNmhsiH1sTOTLPxVsaezo6kRdaVcRSsaUXcazrROrhTaWlo3NCK2A7cDn184lpmPRMQF4G7giqEtIp4P7Aa+tOzUf4yIbuAR4IGFQLfC9fcD91/3LyBJkqQb0uzcPMdOT3C4XDEbG69x6HiNo6cnmSsbUXd2BHfuHuBrb9nKy75+H6PlM2f7d/XTZSNqtbFWN3uoll+XP292Dth6pQsjYhD4EPD2zBxrOvWXgP9Zfv8y4Ncj4vsy83eWf0ZmPgg8WH6ey2ySJEmbzNx88sSZycaK2ULlxkdPTjAzt9iI+o6d/QwPVfnOr3tOY2vjnbsH2NJlI2rdeFod2mrl123Ljm8HLlzuooi4BfgkRdGSNzafy8xPN739zYj4S8ArgWeENkmSJG0O8/PJU+emyhWzOmPls2dHxutcvLRYn+7WHX2MDFX5ttHBxjNnd+2p0NdjONPm0dLQlpnnIuJx4AXAnwJExAGKVbYvrnRNROwHPg18JDN/fBU/Zh5wc7EkSdImkJmcuHCxsXJ26HiNw+NFSJucmWuM27u1l5G9Vf7igV2NUvrDgxUGtrR6DUJqP2tVPfLVLFaPfA9Qzcz7Vhj7XOBTwPsy800rnP86oJ8iACbw3cC/B/5mZv72VeZhIRJJkqQ2kZmcqs80VswOn6iX2xtrXJi+1Bi3u9LD8GC17He2WFJ/W1/3Bs5eWnsb0aftAYo+bVsotj3en5mnyhL/78rMSjn235XjJpZ9zGsz89cj4sXAvwb2AzMUhUjenpkfXMU8DG2SJEkb4NzkTGPFrGhEXWNsvM6ZicVG1Nv6uhkdKoLZ6N4qw4PFc2e7bEStm9S6hrZ2YWiTJElaW7Xp2caK2aGmkvona4uNqCtbuopgNlRleKjKaFlWf4+NqKUlDG2SJEm6bpMzl8pG1EWlxsMnahw+XuPp89ONMb3dHeVzZlVG91YaAe0523oNZ9IqGNokSZJ0VdOzczxyss5Yczg7UeeJs02NqLs6uGtPhdGmJtQjQxVu29FPh42opetmaJMkSVLD7Nw8j52aaKyYLaygHT09QdmHmq6O4MCegSVbGkeGqty+00bU0lowtEmSJN2E5uaTY6cnGitmC6tnj52aYHau+HtSR8D+XQNLnzvbW2X/rgF6ugxn0noxtEmSJG1i8/PJk2enONwoCFLj0Ik6j5ysM1M2oo6A23b0N1bMRsrKjXftqdDbbSNqaaMZ2iRJkjaBzOSr56cbK2aHjtcZGy+qNk7NLjaivmVbb6MB9cIzZwcHK/T32IhaaleGNkmSpBtIZnKydpHDZQn9hSbUYyfq1C4uNqLeU92y2Ous3No4PFRha6+NqKUbjaFNkiSpTZ2ZKBpRj43Xiq9lUDs/NdsYs6O/e8mq2cL3OwZ6NnDmklrJ0CZJkrTBzk/NlitmzeX0a5yqzzTGVHu7loSzhdWz3ZUee51Jm5yhTZIkaZ1MXLzE2Hh9sZz+eJ3Dx2scv7DYiLq/p7PocTZYYXRvtVFWf2jrFsOZdJMytEmSJLXY9OwcRxbCWdPq2ZNnpxpjtnR1cHBwccVsYWvjvu19NqKWtIShTZIk6TrNXJrn0VP1IpgdX9zW+PiZyUYj6u7O4K49lXLFrFIGtKIRdafhTNIqGNokSZKu4tLcPEdPTzZC2UJBkKOnJrhUprPOjmD/rv6yx1mxpXF0b4U7dg3Q3WkjaknXz9AmSZJUmptPnjgz2VQMpNja+OjJCWbmFhtR37GzvxHMhsttjQf2DLCly0bUklrP0CZJkm46mclT56YaK2YLIe3IeJ3p2fnGuH3b+4pnzfZWGRmsMrq3yl17KvT1GM4krR9DmyRJ2rQykxMXLi4po3/4RJ2xEzUmZuYa44a2bnlGr7PhoSqVLV0bOHtJKhjaJEnSpnCqfrFRDORQGcwOn6hxYfpSY8zuSg/Dg0UwWygIMjpUZVt/9wbOXJKuzNAmSZJuKOcmZ57RhPrwiTpnJhYbUW/r617yvNnCCtquypYNnLkkXR9DmyRJaku16VnGxosVs0PH64yN1zh0vMZ47WJjzEDZiHp0qFo8dzZU9D3bU7URtaTNw9AmSZI21NRM0Yj60IlaEdDKkvpPnVtsRN3b3VFua1x85mxkb5VbtvUaziRteoY2SZK0Li5emuOR8YnGitnCFscnzk6y8K/lns4O7hqsLAazMqTdtqOfDhtRS7pJGdokSVJLzc7Nc/TURFlKv95YPTt2epK5shF1V0dw5+6BpRUb91a5Y2c/XTailqQlDG2SJOm6zM0nx05PPKMoyGOnJpidK/492xGwf9cAw+WzZgsVG+/cPUBPl+FMklbD0CZJkq5ofr5oRH246XmzQ8drPHKyzsVLi42ob93Rt6QgyMhQ0Yi6t9tG1JL0bBjaJEkSUDSiPn5hulg5O764cjY2XmeyqRH1c7b1lhUbK43KjQcHKwzYiFqS1oShTZKkm0xmcrJ+kbEl2xqL72tLGlFvYXRvZUnVxuGhKtv6bEQtSevJ0CZJ0iZ2dmKmUUr/8InFsvpnJ2cbY3b0dy8tCFJ+v2OgZwNnLklasK6hLSI6gbcAPwj0Ap8AXpuZpy4z/j7gF4EDwCPAj2bmJ5rOHwTeCfxF4CzwLzLzF1cxD0ObJGlTuTA92whmC6tnh47XOVVfbERd3dJVFATZW20KaVV2V3rsdSZJbWy9Q9tPA68B7gNOA+8F+jPzO1cYewD4MnA/8B+Avw48CDwvM4+WAfDLwKeANwDPBX4H+H8y8zevMg9DmySpZd700S/xG3/8BHOZdEbw/d9wGz/38uevyc+auHiJsfEimBWl9IuS+l89P90Y09/TyfBgU5+zsjDI3q02opakG9F6h7ZjwJsz8z3l+7uAI8D+zDy2bOzPAt+emS9qOvYZ4FOZ+bMR8WLgvwCDmVkvz/8z4Jsz88VXmYehTZLUEm/66Jf4wB89/ozjr/rG259VcJuenePIeL1sRL3Y6+zJs1ONMT1dHUvDWbm1cd/2PhtRS9ImcqXQ1tISUBGxHbgd+PzCscx8JCIuAHcDx5Zdcnfz2NJD5fGF84cXAlvT+X/YwmlLknRFv/HHT1z2+GpC28yleR4rG1GPNRUFOXZ6grIPNd2dwYHdFb7+tu38jXtvKyo27q1y+85+Og1nknRTa3Xd3mr59fyy4+eArZcZv9LY513l/EqfRUTcT7HVUpKklpm7zM6N5ccvzc1z9PRkY8VsrCwKcvTUBJfKdNbZEezf1c9z91b53rtvKXqeDVXYv3uA7k4bUUuSnqnVoa1Wft227Ph24MJlxl9p7NXOL5GZD1I8E0dEuDdSktQSnRErBrcO4Fd+70hZEKTGoycnmJkrGlFHwG07+hkZqvLS5w01tjce2DPAli4bUUuSVq+loS0zz0XE48ALgD+FRrGRrcAXV7jkYWD5s2n3AJ9uOj8SEQOZOdF0/uFWzluSpJVkJk+fn+ZbRnbze4dOPuP8PPC2jx9i3/Y+hocqfOvIniWNqPt6DGeSpGdvrapHvprF6pHvAaqZed8KY+8CvgT8HeC3gL8G/FueWT3yE8BPAqPAx4F/lJkfvMo8LEQiSVqVzGS8dnGxAfXxGofHi+2N9YuXnjE+gP/zwE7ecN9zGR6sUO21EbUk6dnZiD5tD1D0adsCfBK4PzNPRcQrgXdlZqVpfHOftkeB163Qp+1dFH3azgHvyMy3r2IehjZJ0jOcrl9c0udsIaidn1psRL1roKfodTZUbRQEGRmssq3fcCZJWhvrGtrahaFNkm5u5ydnOTxeK3ud1Tl0vMbYeI1T9ZnGmK29XYwsBLOhStnrrMruypYNnLkk6WZkaJMkbVr1i5caTaibV9BOXLjYGDPQ09l41my47HM2urfKYHWLjaglSW3B0CZJuuFNzRSNqJdva3zq3GIj6t7uDoYHq42tjSNDVUb2VrllW6/hTJLU1gxtkqQbxsVLczx6cmJJMDt8osbjZyZZ+GO9p7ODA3sGGitmw4PF6tltNqKWJN2gDG2SpLYzOzfP0VMTHC4bUC80pD52epK5pkbUB3YPlM+dLRYG2b+rny4bUUuSNhFDmyRpw8zNJ4+fmSxWzY7XODxelNR/9FSd2bniz+kI2L9rgOHBSrFyVj5/dufuAXq6DGeSpM3P0CZJWnPz88lT56YYG69x6Hi9sXJ2ZLzOxUvzjXG37ugrnjUbqjJSFgU5OFiht9tG1JKkm5ehTZLUMpnJ8QvTHD5RBrNy9ezIiRoTM3ONcXu39hYl9MvnzUb2FuGssqVrA2cvSVJ7MrRJkq5ZZnKqPtNYMWsup1+bvtQYt7uypbFitrB6NjxUZVufjaglSVotQ5sk6YrOTswUgWx8cfVsbLzOmYnFRtTb+7sZWSinv7faCGk7B3o2cOaSJG0OhjZJEgAXpmcZO/HMXmcna4uNqKtbuhrBbHhwoddZhT0VG1FLkrRWDG2SdJOZnLnEkfF6Y8Xs0PGipP7T56cbY/q6OxtbGUeHFlfQ9m61EbUkSevN0CZJm9T07ByPnKw3VswWnj978uzUYiPqrg4O7qkUz53trTJSrp7duqOPDhtRS5LUFgxtknSDm7k0z2OnJp6xrfHY6QnKPtR0dQR37aksaUI9MlThjl0DdBrOJElqa4Y2SbpBXJqb59iZybIYSJ3D40VD6sdOTXCpTGcdAft3DxQrZnuLYDY6VGX/7gG6O21ELUnSjcjQJkltZn4+efLsFIfL7YzFtsY6j5ysM9PUiPq2nX2MDlWbyulXObBnwEbUkiRtMoY2SdogmcnT56c5vBDMjtcZG68xdqLO1OxiI+pbtpWNqIeqDA8WBUEODlbo77ERtSRJNwNDmyStsczkZO0ih0/Um1bOahw5Uad2cbER9WB1yzOaUI8MVaj22ohakqSbmaFNklrozMRMWUq/LAhSPnt2bnK2MWbnQE9jxWyhpP7IUIXt/TailiRJz2Rok6TrcH5qtrFiNnai3ghqp+ozjTFbe7vK5tNVRgYrje93V7Zs4MwlSdKNxtAmSVcwcfESY+N1Dh+vNRUGqXP8wmIj6oGeTg4OVRkdqiwpCjK0dYuNqCVJ0rNmaJMkikbUR8aLFbPD44urZ0+dm2qM6e3u4OBgZUk5/eHBKvu224hakiStHUObpJvKxUtzPHpyYkkT6rETNY6dmWThj4Wezg4O7BlYUhBkdKjKbTv7bUQtSZLWnaFN0qZ0aW6eo6cnioqN5fNmh47XOHp6krmyEXVnR3Dn7oHGitloWVZ//65+umxELUmS2oShTdINbW4+eeLM5JKVs8Mnajx6coKZuaIRdQTcvrOfkXLFbHioqNx45+4BtnTZiFqSJLU3Q5ukG8L8fPLUuamylH5ZGGS8xpHxOtOz841x+7b3laX0i2fPRvdWuWtPhb4ew5kkSboxGdoktZXM5MSFi00rZzUOnahz5ESNiZm5xri9W3uLFbOFao17qxwcrFDZ0rWBs5ckSWo9Q5ukDXOqfrFRSv9wU1n9C9OXGmN2V3oaz5sthLThwSrb+rs3cOaSJEnrx9Amac2dm5wpCoKcqBUNqY/XGBuvc2ZisRH1tr7uJc+bDQ8WlRt32YhakiTd5NYttEVEP/CvgVeUhz4E/HBmTl1m/KuBHwK+BpgD/gT4icz8Unl+P/AYMAksTPRcZt66irkY2qQ1UJuebWpEXW9sbxyvXWyMqWzpWlwxKwuDjAxV2FO1EbUkSdJKrhTaWv1gyC8BzwVGKULWR4F3AH//MuOrwD8F/hdwCfgnwCci4q7MnGwaN5qZT7Z4rpKuYHLmUqMR9dh4Gc6O13j6/HRjTF93J8NDFV40vIfRvUWvs5GhKrds6zWcSZIktUjLVtoiog84A3xPZn66PPYS4GPAzsycvtL15fheYAp4YWY+1LTSdtu1hjZX2qTVmZ5d3oi6WEF74mxTI+quDu7aU2G0qQn1yFCVW3f00WEjakmSpGdtvVbaRoFe4PNNxx4C+oAR4Iur+IyXUGyFHFt2/I8jogf4CvAzmfnfV7o4Iu4H7r+2aUs3h9m5eR47NbHY56wsp3/01ARlH2q6OoIDewZ4/q3b+KsvuLWxenbHThtRS5IkbZRVrbRFxPuA11xhyD8HPg78PtCxsMQVER0Uz6q9KDP/4Co/YwT4A+CfZOY7y2MV4Osowl838LeBtwLfkJlXDIGutOlmNTefHDs9seR5s7ETdR49VWd2rvhnoiNg/66Bpc+d7a2yf9cAPV2GM0mSpPX2rAuRlOGp9wpDJilW074A7MjMc+V1Oyi2TN59pZAVEV8LfBL4pcx861Xm8kngjzLzH19lnKFNm9r8fPLk2amylH6tURjkyMk6M5cWG1HftrNvSUGQ4aEKd+2p0NttI2pJkqR28ay3R2ZmHahf5YccAqaBFwC/Wx6+h+IZtcNXuO4FwO8A/ywzf3kV05kHfIhGN43M5Kvnp5c8b7awejY1u9iI+pZtvQwPVfnm4d0MDxYl9Q8OVujvsRG1JEnSjazVJf/fTVG+f6Hk/4eBL2fmD11m/DcB/5mizP+7Vzj/jRRh8c8pAuargX8FfFNmfn75+GXXutKmG0pmcrJ2ccm2xoVwVru42Ih6T3ULI0MVRspiICNlcZCtvTailiRJulGtZ5+2AeCXuUyftoj4KeCVmfm88v3vAd9Ksb2y2Xdm5mci4vuBfwY8h2IV7/8H3pyZn1zFXAxtaltnJmaWVms8XufweI1zk7ONMTv6uxeD2d4qI4NFUNsx0LOBM5ckSdJaWLfQ1k4MbWoHF6ZnGTtR49Dx+pLtjafqi42oq71dS1bNFp4/213psdeZJEnSTcLQJq2xiYuXigbUx8tgVn5//MJie8KBnk4ODhUrZqN7q2Uj6gp7t9qIWpIk6WZnaJNaZHp2jiPj9SUFQQ6fqPHk2anGmC1dHRwcbC6lX2F4sMq+7TailiRJ0soMbdI1mrlUNKI+dKJWbm+sMTZe59jpxUbU3Z3Bgd2VxvNmC73Obt/ZT6fhTJIkSdfA0CZdxqW5eY6enmRs2crZY6cmuFSms86OYP+ufkaaep2NDFXYv3uA7k4bUUuSJOnZM7Tppjc3nzxxZrIooT9e51D57NmjJyeYmSsaUUfA7Tv7GwVBFoqDHNgzwJYuG1FLkiRp7RjadNPITJ4+P83h4zUONZXUPzJeZ3p2vjFu3/a+IpjtrTIyWGxrvGtPhb4ew5kkSZLWn6FNm05mMl672FgxWygMcmS8Tr2pEfXQ1i3PaEI9MlSlsqVrA2cvSZIkLWVo0w3tVP1i2YB6sZT+4RM1LkwvhrPdlR6GyxWz4YWtjYNVtvV3b+DMJUmSpNUxtOmGcH5ylsPjZaXGE7WycmOd0xMzjTHb+rqXPG+2sIK2q7JlA2cuSZIkPTuGNrWV2vRsUyPqOmNlUBuvXWyMqWzpKlbMBqvFc2dlUBusbrERtSRJkjYdQ5s2xNRM0Yi60eusXDl76txiI+re7g6GBxdXzIqAVuWWbb2GM0mSJN00DG1aUxcvzfHoyYlGQZBDx4vVs8fPTLLwf0FPZwcH9gwwWoay4cEKo3ur3Lajnw4bUUuSJOkmZ2hTS8zOzXP01MSSJtSHT9Q4enqSubIRdVdHsH/3AKNDRUGQ0bIh9f5d/XTZiFqSJElakaFN12RuPnn8zOQzCoI8eqrO7Fzxv2lHwB27BhrPmg0PVRkdqnLn7gF6ugxnkiRJ0rUwtGlF8/PJU+emGj3OmhtRX7y02Ij6tp19SwqCDA9WOThYobfbRtSSJElSKxjabnKZyfEL0+XK2WI4GxuvMzkz1xj3nG295YrZYhPq4cEKAzailiRJktaUoe0mkZmcrF9k7ES9CGgLPc/G69SaGlHvqW5Z1uusCGlbe21ELUmSJG0EQ9smdHZipqkYyGJZ/bOTs40xO/q7yxWzoiDIQkjbMdCzgTOXJEmStJyh7QZ2YXq2KAZyvF5uaSy+P1VfbERd3dK1pAH1SFm5cU/FRtSSJEnSjcDQdgOYnLlUbGtsNKKuM3aixlfPTzfG9Pd0MjzYtK2xDGp7t9qIWpIkSbqRGdrayPTsHI+crC9WbDxe4/B4jSfOTDXG9HR1cHBP0Xx6uGlr477tfTailiRJkjYhQ1sbePD3H+GDn32Co6cnKPtQ090Z3Ll7oKkgSLFydseuAToNZ5IkSdJN40qhzVru6+CjX3iKX/m9Rzg/NUt1SxeveME+XvWNd7B/9wDdnTailiRJknR5JoY19tEvPMUbP/wlzk8VVR1rFy/xHz73JF95+oKBTZIkSdJVmRrW2Ns+foip2bklx6Zm53jbxw9t0IwkSZIk3UgMbWvs6XNT13RckiRJkpoZ2tbYLdv7rum4JEmSJDUztK2x1790lL7uziXH+ro7ef1LRzdoRpIkSZJuJC0NbRHRHxHvjYhz5es9EXHZJaWI+MGImI+IetPrN5aNuTciPhsRkxHxSES8qpVzXmsvv2cfv/CK57Nvex8B7Nvexy+84vm8/J59Gz01SZIkSTeAlvZpi4h3A88Dvg9I4KPAw5n59y8z/geBN2Xmwcuc3wYcAd4O/EvgW4CPAN+RmX94lbm0VZ82SZIkSbqcK/Vpa9lKW7mi9irgH2fmicwcB/4x8JqI6L3Oj30FMAm8NTMvZuYnKULb/S2ZtCRJkiS1uVZujxwFeoHPNx17COgDRq5w3W0RcTwinoiID0bEnU3n7ga+sGzJ7KHy+DNExP0R8bmI+Nz1/QqSJEmS1F5WFdoi4n0RkVd4/RxQLYefb7p04futl/no3weeD9wC/B/ANPDJiBgoz1eXfR7Auct9XmY+mJn3Zua9q/m9JEmSJKndda1y3A8DP36F85MsrqZtowhWC98DXFjposx8tOnt8Yj4exQh7RuBTwM1YP+yy7Zf7vMkSZIkabNZVWjLzDpQv9KYiDhEsVL2AuB3y8P3AFPA4VXOJ8vXwgN4DwMvXzbmnvK4JEmSJG16LXumLTOngA8Ab46IwYgYBN4MvD8zp1e6JiK+OyJujcJO4FeAU8AflUM+AgxExOsjoiciXkJRnOTBVs1bkiRJktpZq5tr/wjFqtrC6xDwuoWTEfFTEfGVpvHfBnyWYhXvK8AuinL+dYDMPAd8F/DXKbZNvhv4oauV+5ckSZKkzaKlfdraiX3aJEmSJN0o1qVPmyRJkiSp9QxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxgxtkiRJktTGDG2SJEmS1MYMbZIkSZLUxloa2iKiPyLeGxHnytd7IqLvCuPfGRH1Za+MiB9tGpMRMblszLZWzluSJEmS2lVkZus+LOLdwPOA7wMS+CjwcGb+/VVe/x3AfwXuyMyny2MJvCgz/+Aa55Kt/N0kSZIkaa1EBJkZK55rVbApV9TOAN+TmZ8uj70E+BiwMzOnV/EZvwV0ZOYrmo4Z2iRJkiRtalcKba3cHjkK9AKfbzr2ENAHjFzt4ojYC7wMeOcKp/9jRJyKiD+OiFescH7hM+6PiM9FxOeubeqSJEmS1J5WFdoi4n3ls2WXe/0cUC2Hn2+6dOH7rav4MX8HeBz45LLjfwm4E7gVeAfw6xFx30ofkJkPZua9mXnvan4vSZIkSWp3q9oeGREVilW0y5mkWE37ArAjM8+V1+2g2DJ5d2Z+8Qqf3wE8CvybzHzgKnN5N9CbmT9wlXFuj5QkSZJ0Q7jS9siu1XxAZtaB+lV+yCFgGngB8Lvl4XuAKeDwVX7EfcBzgPeuYjrzwIq/jCRJkiRtNqsKbauRmVMR8QHgzRHx5fLwm4H3r6IIyWuBD2fmyeaDEfF1QD/wpxTVKL8b+AHgb7Zq3pIkSZLUzlrdXPtHKFbVFl6HgNctnIyIn4qIrzRfEBH7KMLYSgVI9gD/DjgLjANvAv52Zv52i+ctSZIkSW2ppX3a2onPtEmSJEm6UaxXyX9JkiRJUosZ2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY0Z2iRJkiSpjRnaJEmSJKmNGdokSZIkqY21NLRFxP8bEX8cEZMRcWSV19wXEV+JiKmI+HJE/OVl5w9GxKciYiIinoyIH2vlnCVJkiSpnbV6pe1p4K3AP1/N4Ig4AHwY+AVgW/n1IxGxvzzfCXwM+DNgD/BXgDdExN9o8bwlSZIkqS21NLRl5m9l5oeAp1Z5yWuAz2fmBzJzJjN/HXioPA7wLcAdwBszczIzHwLeBfxQK+ctSZIkSe2qa4N//t3A55cde6g8vnD+cGbWl53/hyt9WETcD9zf9L51M5UkSZKkDbCq0BYR72Nx9Wsl/zwz33QdP78KnF927BzwvKuc37rSh2Xmg8CD1zEPXYeIyMw0GWtDeP9po3jvaSN5/2mjeO9trNWutP0w8ONXOD95nT+/RvEsW7PtwIVVnpckSZKkTW1Voa3cnli/6sBr9zDw4mXH7gE+3XR+JCIGMnOi6fzDazAXSZIkSWo7rS753xURvUB38TZ6y/eX837g3oj4/ojojojvB14I/Gp5/veBY8DPR0RfRHw98FqKYiTaeD+70RPQTc37TxvFe08byftPG8V7bwNFZrbuwyJ+Bviny48v7H+NiFcC78rMStM19wG/CBwAHgVel5mfaDp/kCKk/UWK59nekZlvb9mkJUmSJKmNtTS0SZIkSZJaq9XNtSVJkiRJLWRokyRJkqQ2ZmjTqkVEf0S8NyLOla/3RETfFca/OiL+V0ScjYhTEfHfIuL56zlnbQ7Xce/dXd5vxyMiI+Kb13O+urFFRGdEvC0iTkZELSI+FBG7rzD+voj4SkRMRcSXI+Ivr+d8tXlcy70XEfsi4j9FxLHyz7lXrfd8tblc4/33XRHxu+Xf785GxGci4kXrPeebiaFN1+KXgOcCo8AI8DXAO64wvkpRmOZWYB/wEPCJiOhf43lq87nWe28G+DDwPWs/NW1CPwm8DPgGij+/AH5tpYERcYDiXvsFir6ivwB8JCL2r/00tQmt+t4D5oFPAP838OTaT003gWu5/3YAvwwcBPYA/x74bxFx21pP8mZlIRKtSrmqcQb4nsz8dHnsJcDHgJ2ZOb2Kz+gFpoAXZuZDazlfbR7P9t6LiARelJl/sOaT1aYQEceAN2fme8r3dwFHgP2ZeWzZ2J8Fvj0zX9R07DPApzLT8ti6Jtdy7y277ijwpsz8wLpMVJvS9d5/TdcfB/5BZn54bWd6c3KlTas1CvQCn2869hDQR7HysRovASaBsdZOTZtcK+49aVUiYjtwO033W2Y+AlwA7l7hkrtZem9CcX+uNFa6rOu496SWebb3X/n4y27gS2s0xZueoU1ExPvK/fCXe/0cxVZHgPNNly58v3UVP2ME+HfAj2VmrbW/gW5U63HvSddopfsNij6hK91v1WsYK13Jtd57Uitd9/0XEYPAh4C3Z6b/YX6NdG30BNQWfhj48Sucn2RxRWMbxT/AC99D8V9hLisivhb4JMU/zO+8/mlqE1rTe0+6Dgv/UWnbsuPbWfl+q13DWOlKrvXek1rpuu6/iLiF4u94nwDeuCYzE+BKm4DMrGfmqSu8JoFDwDTwgqZL76F4Ru3w5T47Il4A/HfgLZn51jX8NXQDWst7T7oemXkOeJym+60sNrIV+OIKlzzM0nsTivvz4TWaojap67j3pJa5nvuvLLj0GeC/ZeYPp4Uy1pShTauSmVPAB4A3R8RguRT+ZuD9lysEERHfBHwa+OnM/OX1m602k+u89yIiesviNwA95fvOdZq2bmwPAm+IiDsjYivwAPDxzDy6wtj3A/dGxPdHRHdEfD/wQuBX12+62kSu5d6j6c+5ALrL9+6i0vVa9f0XEc8F/gD4jcy80o4ZtYihTdfiRyhWNhZeh4DXLZyMiJ+KiK80jf85imX2fxER9aaXfTx0rX6Ea7v37qBYiZsq33+6/P4H1mOyuuG9haI66Z8ATwGdwKsAIuKVEVFfGFg+qP8K4E0UW4jeBHzf5f6SLV3Fqu+90sKfc7cD7y2/f9O6zVabzbXcf2+gaOf0I8v+jvfK9Z70zcKS/5IkSZLUxlxpkyRJkqQ2ZmiTJEmSpDZmaJMkSZKkNmZokyRJkqQ2ZmiTJEmSpDZmaJMkSZKkNmZokyRJkqQ2ZmiTJEmSpDb2vwFmLXwyWZ0B+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use neutral structure and datapoints of disp_002,004,006\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from nonrad.ccd import get_Q_from_struct\n",
    "from nonrad.elphon import get_Wif_from_WSWQ\n",
    "from pathlib import Path\n",
    "from pymatgen.core.structure import Structure\n",
    "from glob import glob\n",
    "\n",
    "cc_dir = \"/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad\" + \\\n",
    "         \"/\" + \\\n",
    "         \"/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited\"\n",
    "\n",
    "SnZn_4 = Structure.from_file(f\"{cc_dir}/4/CONTCAR\")\n",
    "SnZn_5 = Structure.from_file(f\"{cc_dir}/5/CONTCAR\")\n",
    "SnZn_6 = Structure.from_file(f\"{cc_dir}/6/CONTCAR\")\n",
    "\n",
    "\n",
    "WSWQs = []\n",
    "for d in glob(f\"{cc_dir}/*\"):\n",
    "        if d[-2:] in [\"/4\", \"/5\", \"/6\"]:\n",
    "            pd = Path(d)\n",
    "            Q = get_Q_from_struct(SnZn_5, SnZn_6, f\"{pd}/CONTCAR\")\n",
    "            path_wswq = f\"{pd}/eph/WSWQ\"\n",
    "            WSWQs.append((Q, path_wswq))\n",
    "print(WSWQs)      \n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "ground_vr_path = f\"{cc_dir}/5/vasprun.xml\" \n",
    "\n",
    "Wifs = get_Wif_from_WSWQ(WSWQs, ground_vr_path, 2305, 2303, spin=1, fig=fig)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.ylim(-1e-6,1e-6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6e659587",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/6/eph/WSWQ'), (-0.0, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/5/eph/WSWQ'), (-0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/4/eph/WSWQ')]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAFgCAYAAAAo85E9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA31ElEQVR4nO3dfXzd913f/dfH8r0tnWPHdmJbchwnjXNnS27LVe7KKB1ruRvQwWPr1dKWMVI2uK61jK6ldBsUGL2BMgZsbVhKKWXARm8Y2yC9gV0r26C0iRQnzU2TNLHkm8R2fI4k2/KN9Ln++P0kH6myLDtH0pH0ej4e52Gf3+93zvnqkZPEb3+/v/c3MhNJkiRJUmtasdADkCRJkiRdnqFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkrQsRcSaiLg3Ip6JiKGI6I2I7yjP3RERX4yIU+XjsxFxR8NrIyLeGxEny8d7IyLKc7dGxB9HxPGIeD4i7ouIvQv1c0qSFj9DmyRpuVoJ9AN/C6gA7wL+U0TsBo4APwBsBrYA/wX4g4bX3g18H9AN7Ae+B3hzea5aXr8XuB74AvDHc/mDSJKWtsjMhR6DJEktISIeBH4uMz/ecGwlRSB7f2auL4/9b+AjmXlP+fxHgB/NzK+f5j03AyeBLZl5ch5+DEnSErPkZ9oi4h9ExOcjYjAiLs7xZ22LiN8pl8oMlkttdszlZ0qSmiMirgduBR5uOFYDRoBfB/51w+V3An0Nz/vKY9P5FuCYgU2SdK1WLvQA5sEp4N8B64B75upDImIt8DngryiWxDwP3A4Mz9VnSpKaIyJWAb8H/E5mPjp+PDOrEbEBeCPwTMNLNgL1hud1YGNERDYsYYmITuA3gZ+cy/FLkpa2JT/Tlpn3ZebvA09Ndz4ivi8ivhQRtYh4JCJed40f9UaK+xj+SWaeyMyxzHw4Mwev8f0kSfMgIlYAvwucB35i6vnMPA18EPhoRGwrDw8DHQ2XdQDDUwLbVuDTwL8r/z8kSdI1WfKhbSYR8e3AvcBbKG42fyPwGxHxLdfwdq8AvgJ8pFwe+WhEvLVpg5UkNV3Z+HgvRWHI38vMC5e5dAWwHthZPn+YooRkXDeTl1Vuoghs/yUzf7HZ45YkLS/LOrQB/xT4tcz8fDkz9gXgY8AbACJid0TkFR7jS0y3UAS3LwDbgdcDP/MCZu4kSXPv31MsZf+ezDw7fjAivj0iDkREW0R0AB+gWG7/SHnJR4GfjIid5b3L/wz4SPnaDuA+4H9l5jvm70eRJC1Vy+GetpncBLwiIhrvNWgDPl/+/hCwdaY3yMzxcpMh4HBm/lr5/IsR8THgeynuk5AktZCIuJGiFfIccKzcZo3y2HmK8pFO4CzFX8i9OjNHyms+BOwBDpbP/0N5DOD7ga8D7oyINzV85B2ZeWhOfhhJ0pK2bCr/I+Jbgc9m5sqGY/8d+IvMfH8T3v9fAf8oM7sajv0bYHtm/v0X+v6SJEmSlqclvzyyXNqyFlhdPl9bPgL4N8BbI+Ll5XWrI+IlEfHSa/iojwDXRcSPl+/VDbwO+ESTfhRJkiRJy1DTQ9vV7osWES+NiC9ExJmIeDIiXj/l/LaI+EREDEXE8Yh4b9n0NVs/RLG05T6KpY9ny8eNmflp4EeB9wMngKPAr1JUOV+VzHwG+E7gHwGDwB8BP5uZf3i17yVJkiRJ45q+PDIiXkXRxLgOuKdxOeI011aAJ4Bfppj1+hbgk8C3Z+b/Ka/5DMX9Yj8MXAf8GXBvZr63qQOXJEmSpBY0Z/e0TXcP2TTX/DDws8Du8b1tIuJ3gYuZ+cMRcRPF/mq3ZOaT5fkfAd6VmTfNycAlSZIkqYUsdHtkN/BATk6O91MsaRw/Xx8PbA3nd0dEx9SNqyPibuDu8ulL5mjMkiRJktR0mRnTHV/o0NYO1KccqwEdVzhPec2k0JaZ9wD3AERELpdmTEmSJEmLW8PWM19jodsjh4DKlGNVLoWxy50fPydJkiRJS9pCh7Y+oGfKsQPl8fHzlYjYM+X805k5dQZOkiRJkpacuaj8n2lftKk+CWyIiLeVe6S9EngN5RLHzPwq8FngfRHRURaTvB34ULPHLUmSJEmtaC5m2i67L1q5ifVwROwCyMwaxd5mP0hx79pvAT82Xvdfel05zsPA3wB/DLxvDsYtSZIkSS1nzir/F5pFJJIkSZIWi4i4bHvkQt/TJkmSJEmagaFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkiRJklqYoU2SJEmSWpihTZIkSZJamKFNkiRJklqYoU2SJEmSWpihTZIkSZJa2MqFHoAkSZIkzYdPPXCY99/3GEdqZ9lRXcfbXrWX7zuwc6GHdUWRmQs9hjkREblUfzZJkiRJV+dTDxzmpz9xkLMXRieOrVvVxi+9Zl9LBLeIIDNjunMuj5QkSZK05L3/vscmBTaAsxdGef99jy3QiGbP0CZJkiRpyTtSO3tVx1uJoU2SJEnSkrejuu6qjrcSQ5skSZKkJe9tr9rLulVtk46tW9XG2161d4FGNHu2R0qSJEla8sbLRmyPbCG2R0qSJElaLGyPlCRJkqRFytAmSZIkSS3M0CZJkiRJLczQJkmSJEktrOmhLSLaIuL9EXE8IoYi4uMRseUy174zIoanPDIi/m3DNU9HxMiUa/Y1e9ySJEmS1IrmYqbtHcD3Ai8DOstjvzvdhZn5rzNz4/gDeDGQwMemXPqPGq/LzINzMG5JkiRJajlzsU/b3cC7M/MpgIj458ATEXFjZj4zi9c+kJlfmINxSZIkSdKi09SZtoioAruAL40fy8wngUGg+wqvXQO8CfjQNKc/EBHPR0RvRLx5hve4OyK+GBFfvIbhS5IkSVLLafZMW3v5a33K8RrQcYXX/gCwGviPU46/kSIEngO+FfiDcuO5rwl3mXkPcA8Um2tfzcAlSZIkqRU1+562ofLXypTjVYrZtpm8Gfi9zBxuPJiZ/19mDmfmhcz8DPAB4PXNGKwkSZIktbqmhrbMrAGHKApFAIiIPRSzbA9e7nURcQfwcuCDs/iYMSBe0EAlSZIkaZGYi/bIe4C3R8RNEdEBvBe4LzOfnuE1bwb+KjP7Gg9GxI0R8YqIWFtuJfC3gLcCfzgH45YkSZKkljMXoe09wJ8AfwMcBtoolzNGxOsiYtLyx4hYB7yB6WfZNlAshzwOnAJ+k6KZ8tfnYNySJEmS1HIic2n2dURELtWfTZIkSdLSUpYtTnsb2FzMtEmSJEmSmsTQJkmSJEktzNAmSZIkSS3M0CZJkiRp2RgbS/qfP7PQw7gqKxd6AJIkSZI0V47VR+jtr9E3UKP3UI2Dh+tcGB3joZ97FavaFscclqFNkiRJ0pIwNHKBg4frRUjrr9HXX+fY4AgAq9qC27d38P0HdtLTVWV0LFnVtsADniUr/yVJkiQtOhdGx3js2NBEQOvtr/HE8WHGI8BNWzbQ3Vmhu6tKT1eV27d3sLaFU9pMlf+GNkmSJEktLTPpf/4sD/Sfoq+/Tt9AjYcO1zl3cQyAzRtW09NVpbuzSs+uKt2dFarrVy/wqK+OoU2SJEnSovH86fP0DVyaQevrr3HqzAUA1qxcwb6dl2bQerqqdG5aR8S0eWfRMLRJkiRJakkjF0Z5+MjgpYA2UOOZk0W7YwS8aNvGYhatnEnbe0P7oikQuRqGNkmSJEkLbmwseerEML39dXrLpY6PHB3k4ljx5/btlbUTSxz3d1bY31ll45rl0Z1oaJMkSZI0754dHJlUFHJwoM7QuYsAtK9Zyb7OysQsWk9Xles71i7wiBeOoU2SJEnSnBo+d5GDAw11+wM1jtaLuv2VK4q6/e6uCj1dm+jpqrBny0ZWrFjc96E100yhbXnMNUqSJElqmvG6/caykK88d6lu/8br1vN1uzdPzKDduaO16/ZbnTNtkiRJki4rMxk4dXbSMseHjtQZuVDU7W9av2oinI2XhWzesLjq9luByyMlSZIkzUrtzHn6Bur0HqpNzKSdPH0eKOr279pZmSgL6ems0rV58dfttwJDmyRJkqSvMXJhlC8fLer2x2fRnm6o279l68ZJ+6Et1br9VmBokyRJkpa5om7/9KSikEeODnJhtPgz8/UdayY1Oe7bWaF97aoFHvXyYWiTJEmSlpnnhkYmljj29td4sP9S3f6G1W3s7xwPaBW6u6psr6xb4BEvb4Y2SZIkaQk7fe4iBw/XJ5Y49vXXOFLW7betCG67oX3SLNrNWzfSZt1+SzG0SZIkSUvExdExHn92uJhBK2fSHn92iLHyj767Nq9vuA+twp07KtbtLwKGNkmSJGkRGq/bH29x7Ouvc/BwnbMXRgGorl9VNDk2VO5bt784GdokSZKkRaB+5sLEPWjjZSEnhou6/dUrV3Dnjo6JgNbTVWXX5vXW7S8RhjZJkiSpxZy7OMojR4foPXSq2Betv8ZXT5yeOH/Lto2T9kPbe0M7q1dat79UGdokSZKkBTQ2ljx98lLdfm9/jS831O1va59St99ZocO6/WXF0CZJkiTNo+ND5y41OZb3ow2OXKrb39dZ1OwfKIPaDR1rXea4zBnaJEmSpDly5vxFDg7Uy3BWLHM8XDsLFHX7e69vnxTQbtlm3b6+lqFNkiRJaoKLo2N85bnhiVm03v7Jdfudm9ZNCmh37aiwbrV1+7qyeQ1tEdEGvAd4E7AW+DTw5sw8Mc213wr8BXC64fCDmfmNDdfcAnwQ+AbgFPCrmfkrsxiHoU2SJEnXLDM5Uh+Z2Autt7/GwYFLdfuVdauKe9DKpY7dXVW2bFyzwKPWYjVTaFs5B5/3DuB7gZcBJ4EPA78LfMdlrh/NzI3TnSgD4J8AnwX+LnAb8GcRMZCZf9jsgUuSJGn5qp+9wIMD40UhxTLHE8PnAFjdtoI7dnTw97+ui+6uCj1dm9h9nXX7mh9zMdP2DPDuzLy3fH4z8ASwOzOfmXLttwKfzcxpw2NEvAL4b8C2zBwuj/088M2Z+YorjMOZNkmSJE3r3MVRHj06VMygHarRO1DjqeOXFn/t2bph0n5ot93QYd2+5tS8zbRFRBXYBXxp/FhmPhkRg0A38Mw0L2uLiH5gVfm6d2ZmX3muG3h8PLCV7gd+/DKffzdw9wv9OSRJkrR0ZCZfPXF6oijkgf4ajxwZ5PzoGABbNhZ1+3/vxZ10dxZ1+5V11u2rdTR7eWR7+Wt9yvEa0DHN9Y8CPcDDwEbg7cCfR8S+zDxSvt9s34vMvAe4B4qZtqsevSRJkha9E8PnJhWFNNbtr1tV1O2/6Zt2T+yLtqNi3b5aW7ND21D5a2XK8SowOPXizDwGHCuf1oCfjogfoLj/7d7y/Wb1XpIkSVp+zp4f5aEj9Ykljr2HLtXtrwjYe0MH37V/O92dVXp2Vbll60ZWtrnMUYtLU0NbZtYi4hDwYqAXICL2UMyMPTjLtxkDxv+qow+4NSI2ZOb4IuMD5XFJkiQtI6NjyRPPDdPbf4re/jp9/TUee3aI0bJvf2d1HT27qrzpG3cXdfs7O1i/ei5696T5NRdFJD8DvAF4NUV75L1Ae2a+epprvw04BDwFrAd+CngLsC8z+8v2yIcotg14B7AXuA/4p5n5B1cYh0UkkiRJi1RmcrQ+UixzLGfQDh6uc+Z8UbffsXZlUbdfPvZ3Vtnabt2+Fq/5rvx/D7AJ+BtgDfAZ4PXlQF4HfKih4r8b+G1gC8VebfcD356Z/QCZORoR3wN8iCIA1oD3XymwSZIkaXEZHLnAg/31if3QevtrHB+6VLd/+44OfvAlnfTsqtLdWeWmLRu8D03LRtNn2lqFM22SJEmt6fzFMR49Nkhff40HyqKQJxvr9rdsmJhF299Z4Y4dHaxZ2baAI5bm3kwzbYY2SZIkzZnM5JmTZ+gbqPHAoRp9AzUePjLI+YvjdfurixbHsihk/84qlfXW7Wv5MbRJkiRpXpwcPseDA/WJGbS+gRq1MxeAsm5/Z2ViiWN3V4Wd1XUuc5QwtEmSJGkOnD0/ysNH6pf2Qxuo0f/8pbr9W69vn9gLraeryou2WbcvXY6hTZIkSS/I6Fjy5PHhif3Q+vprPHpsSt1+VzF71t1Z5a6dFTassW5fmi1DmyRJkq7K0frZom6/v05v/ykODtQ5Xdbtt69dOXEfWncZ1La1r13gEUuLm6FNkiRJlzU0coGDA/WJ/dD6Bmo8O1jU7a9qC+7Y3sH+znJPtF1VbrpuAytWeB+a1EyGNkmSJAFwYXSMx44NXSoK6a/xxPFhxv/YdNOWDeUsWoXurqp1+9I8MbRJkiQtQ5nJoefPFCUh5TLHh48Mcq6s279uw+pJRSH7OytU169e4FFLy5OhTZIkaRl4/vT58j60YoljX3+NU2Xd/tpVK4q6/TKkdXdW6dxk3b7UKgxtkiRJS8zIhaJuv9iwuk5ff41Dz58BIAJu3Ta5bv/W663bl1qZoU2SJGkRGyvr9sfvQ+vtr/HYsSEulnX7OypryxbHYgZtX2eFjdbtS4uKoU2SJGkReXZwpJxBK9ocDx6uM3zuIgDta1ayv6syUbnf01VlW4d1+9JiZ2iTJElqUcPnLvLgQFEUMj6LdmxwBCjq9m/f3jERzrq7quzZYt2+tBTNFNqcN5ckSZon43X7veN1+wM1vvLcpbr93det5+v3bC72RNtV5Y7tHaxdZd2+tNw50yZJkjQHMpP+589O2rD6ocP1ibr9zeN1+51VursqdHdW2bTBun1puXJ5pCRJ0hw7dfo8vQOXNqzuG6jz/OnzAKxZWdTtj5eF9HRW6dps3b6kSwxtkiRJTTRyYZQvHx2cmEHr7a/xzMlLdfsv2raxuA9tVzGTtveGdlZZty9pBoY2SZKkazQ2ljx1Ypje/jq9/afo66/zyNHBibr97ZW15RLHoizEun1J18LQJkmSNEvPDY5M7IfWN1Djwf46Q2Xd/sY1K9nfWdTt7++scmBXleut25fUBIY2SZKkaQyfu8jBgfqkNsej9aJuf+WKYO8N7fSUM2g9XVVu3rrRun1Jc8LKf0mStOxdHB3jsWcv1e339k+u27/xuvV83e7N5TLHCnfuqFi3L6klONMmSZKWnMxk4NTZSTNoBw/XGblQ1O1vWr+qaHJsKAvZbN2+pAXk8khJkrSk1c6cp2+gPjGD1tdf42RD3f5dOysT+6Ed6Npk3b6klmNokyRJS8Z43X7jfmhfPXEaKOr2b9m6caLJsafLun1Ji4OhTZIkLUpF3f7pSzNoAzUeOTrIhdHi//Hb2tcU4WxXsWH1vs4K7WtXLfCoJenqGdokSdKi8NzQCH0N+6H1DdQYGinq9jesbmP/xH5oFXq6NnFDxbp9SUuDoU2SJLWc0+cu8tDh+sQMWl9/ncO1swC0rQhuK+v2uxvq9tus25e0RBnaJEnSgro4Osbjzw6X4axY6vj4s0OMlf+r7tq8jp6uTXR3Vjiwq2rdvqRlZ15DW0S0Ae8B3gSsBT4NvDkzT0xz7XcCPwXsB9qAh4B3ZubnG65J4Cww1vDSnZlZv8I4DG2SJC2AzORwraFuv7/OwcN1zl4YBaC6flXZ5FjlQFeV/Z0Vrtu4ZoFHLUkLa75D288AbwReDZwEPgysz8zvmOba1wFngL8AhoEfBd4P3J6Z/eU1Cbw8M//yKsdhaJMkaR7Uz1yYNIPWN1DjxHBRt7965Qru3NFR7IdWLnO88br11u1L0hTzHdqeAd6dmfeWz28GngB2Z+Yzs3j9MeCfZOYnyueGNkmSWsS5i6M8cnSI3kOnJvZFe6qs2we4eesGero20dNVoburym03dLB6pXX7knQlM4W2lU3+oCqwC/jS+LHMfDIiBoFuYMbQFhH7gC3AwSmn/nNErAKeBN47Huimef3dwN3X/ANIkqQJY2PJ0ydP0zdQo/dQjd6BOo8cGeT8aHHHwtaybv/vvaSTnq6ibr/Dun1Jarqmhjagvfx16v1mNaBjphdGxDbg48AvZ+ZXGk79beB/lb//XuD3IuL7M/PPpr5HZt4D3FO+n9NskiRdheND58rNqmsT96MNlnX761e3sb+zwg9/0+6JfdFu6FjrMkdJmgfNDm1D5a+VKcerwODlXhQRO4DPUJSW/HTjucz8XMPTP4yIvw28Dvia0CZJkmbnzPmLHByoT1Tt9/bXJtXt772+ne/av2NiP7Rbtlm3L0kLpamhLTNrEXEIeDHQCxAReyhm2R6c7jURsRv4HPDJzPypWXzMGOD/NSRJmqXRseTxZ4cmikKm1u13blpHz64qb/rG3XR3VblrZwfrVzf773UlSddqLv6LfA/w9oj4C4r2yPcC92Xm01MvjIjbgM8CH8nMd01z/i5gPUUATOC7gB8C/sEcjFuSpEUvMzlSHymr9ms80F/jocN1zpwv6vYr61axv7PC37njerrLjau3WLcvSS1trvZpey/FPm1rKJY93p2ZJ8qK/w9l5sby2t8urzs95W3enJm/FxGvAH4D2A2cpygi+eXM/INZjMP2SEnSklc/e4GDA3V6+0/R218sdzw+dA6A1W0ruGNHx0TVfndXld3W7UtSS5rXyv9WYWiTJC015y+O8cjRwYmikN7+Gk8dv/T3nnu2bpgIaD3W7UvSomJokyRpkclMnj55ht7+UxNFIV9uqNvfsnFNGc6KopB9nRUq66zbl6TFytAmSVKLOzF8btJ9aA8O1KmfvQDAulVt7OusTFrmuKNi3b4kLSWGNkmSWsjZ86M8dKQ+qc1x4FRRt78iYO8NHXSXIa27q8qLtm1kZZvLHCVpKTO0SZK0QEbHkieeG56YQevrr/HYs0OMln37O6vrJs2gWbcvScuToU2SpHmQmRwt6/Z7B4qAdnCgzumybr997cpLAa2zCGlb263blyQZ2iRJmhODI+N1+8USx77+Gs+Vdfur2oI7tndMzKD1dFXZfd0GVqzwPjRJ0tcytEmS9AKdvzjGo8cGy/vQin3Rnmys29+yodisurNCz65N3L69nTUr2xZwxJKkxcTQJknSVchMnjl5ZtJ+aA8fGeT8xfG6/dUTSxz3l0Gtun71Ao9akrSYGdokSZrByeFzPDhQnygK6RuoUTvTULe/s0J3uR9ad1eFndV11u1LkprK0CZJUuns+VEePlLch9Y3UCxz7H/+Ut3+rde3091ZpWdXMZN26/XW7UuS5p6hTZK0LI2OJU8eH6b30KU2x0ePXarb31FZW9yHVhaF7NtZYcMa6/YlSfPP0CZJWhaO1Ufo7T81URQytW6/qNkvlzl2VtjWsXaBRyxJUsHQJklacobG6/bLGbTe/hrPDk6u2+9u2A9tzxbr9iVJrcvQJkla1C6MjvHYsaFJ+6E9cXyY8f/M37RlQ1G1Xy51vH17B2tXWbcvSVo8DG2SpEUjMzn0/JkynNXpG6jx0OE658q6/c0birr98YBm3b4kaSkwtEmSWtbzp88X+6EdKqr2+/prnCrr9tesXMG+ncUM2v6uKge6qnRusm5fkrT0GNokSS1h5MIoDx8ZnFji2DdQ45mTZwCIgBdt2zgxg9bTVeXW69tZZd2+JGkZMLRJkubd2Hjdfv/4DFqdR44OcrGs299eWduwxLHKvs4KG63blyQtU4Y2SdKce3ZwZFJRyIMDdYbPXQSgfc1K9ndVik2ry6B2vXX7kiRNMLRJkppq+NxFHixnz3r7T9HXX+fY4AgAK1cEd+zomKja7+mqsGfLRuv2JUmawUyhzXUokqQZjdftN5aFfOW5S3X7u69bz8v2bC5m0XZVucO6fUmSmsqZNknShMxk4NTZScscHzpSZ+TCpbr9/eV+aD3lvWibNli3L0nSC+XySEnStGpnzk/aD62vv8bJ0+eBS3X73eNtjp1VujZbty9J0lwwtEmSGLkwypePDtLXMIv29JS6/fEljt2dVfbeYN2+JEnzxdAmScvM2Fjy1IlhevvrEyGtsW7/ho6Guv2uCvs7q9btS5K0gAxtkrTEPTc4wgMNG1Y/2F9nqKzb37hmJfs7KxMbVvdYty9JUssxtEnSEjJ87iIHB+qT2hyP1i/V7d+2vX2ibv9AV5U9WzfSZt2+JEktbV4r/yOiDXgP8CZgLfBp4M2ZeeIy178a+BVgD/Ak8JOZ+emG87cAHwS+ATgF/Gpm/kqzxy1Jreji6BiPPTs0aT+0rzw3RLnKkV2b1/N1uzdP7Id2546KdfuSJC0xTZ9pi4ifAd4IvBo4CXwYWJ+Z3zHNtXuAh4C7gf8E/CBwD3BnZj5dBsCHgM8CbwduA/4M+H8y8w+vMA5n2iQtKuN1++Mtjr39NQ4evlS3v2n9quIetIaykM3W7c+bd33qIL//1/2MZtIWwWtf1sUvfN++hR6WJGmJmNflkRHxDPDuzLy3fH4z8ASwOzOfmXLtzwHflpkvbzj2eeCzmflzEfEK4L8B2zJzuDz/88A3Z+YrrjAOQ5ukllY7c56+gaIoZPxetBPDl+r279pZmdgT7UDXJuv2F9C7PnWQj/3Voa85/vqv32VwkyQ1xbwtj4yIKrAL+NL4scx8MiIGgW7gmSkv6W68tnR/eXz8/OPjga3h/I83cdiSNOdGLozyyNHBiar9voE6Xz1xGijq9m/eupFv3btt4j406/Zby+//df9ljxvaJElzrdn3tLWXv9anHK8BHZe5frpr77zC+enei4i4m2KppSQtmKJu//TE7Nl43f6F0WL2f1v7Gnq6qvzASzrp6aqyr7NCx9pVCzxqzWT0Mis3LndckqRmanZoGyp/rUw5XgUGL3P9TNde6fwkmXkPxT1xRIT/J5U0L44PnZuYQestg9rQSFG3v2F1G/s6K/zIN++hp6uo3d9eWbfAI9bVaouYNqC1uVxVkjQPmhraMrMWEYeAFwO9MFE20gE8OM1L+oCp96YdAD7XcP7WiNiQmacbzvc1c9ySNFunz13kocP1iRm0vv46h2tnAWhbEdx2Qzvf071jYj+0m63bXxJe+7Kuae9pe+3LuhZgNJKk5Wau2iPfwKX2yHuB9sx89TTX3gwcBH4E+CPgB4D/wNe2R34aeAewF7gP+KeZ+QdXGIdFJJJekIujYzz+7PCkNsfHn71Ut9+1eV3R5FgGtDt3VFi32rr9pcr2SEnSXJrv9sg24L0U+7StAT4D3J2ZJyLidcCHMnNjw/WN+7Q9Bbx1mn3aPkSxT1sN+EBm/vIsxmFokzRrmcnh2tlJ+6EdPFzn7IVRACrrVpV7oRX7oXV3Vrlu45oFHrUkSVoq5jW0tQpDm6SZ1M9cmJhBG1/qOF63v3rlCu7c0TFpFu3G69Zbty9JkuaMoU3Ssnbu4iiPHh2aVBby1InTE+dv3rphomq/u6vKbTd0sHqldfuSJGn+GNokLRtjY8nTJ08Xs2eHavQO1HnkyCDnR8cA2FrW7Y8/rNuXJEmtwNAmack6PnRu0hLHvv4ag2Xd/vrVbezbWaFnV5WezmpZt7/WZY6SJKnlGNokLQlnzl/kocODE0sce/trk+r2X7RtIwd2FTNo3V1VXrSt3bp9SZK0KBjaJC06o2PJV54bovfQ+CxancefHWK07Nvv3LSuaHPsrNKzq8qdOzpYv7qpW09KkiTNG0ObpJaWmRytj0yaQTt4uM6Z81Pq9jsr7C+XOW5tt25fkiQtHYY2SS2lfvYCBwfq9A3UeKCcSTs+dA6A1W0ruGNHx0RRSHdXld3W7UuSpCXO0CZpwZy/OMYjRwcnFYU8efxS3f6erRsmljh2d1a5fbt1+5IkafkxtEmaF5nJ0yfPTFrm+OWGuv0tG1cXs2dlSNu/s0plvXX7kiRJhjZJc+Lk8LlJ+6H19deon70AwLpVbezrrEwKaTus25ckSZqWoU3SC3b2/CgPH6lPzKD1DdTof76o218RcOv17RP3oe3vrHLr9RtZ2eYyR0mSpNkwtEm6KqNjyRPPDRfLHMuZtMca6vZ3VtfR3VUpQ9om7tpp3b4kSdILYWiTdFmZybHBkXKJY1EUcnCgzumybr997cpLSxy7quzvqrCtfe0Cj1qSJGlpMbRJmjA4UtTtTyxz7K/xXFm3v6otuGN7B90N96HddN0GVqzwPjRJkqS5ZGiTlqnzF8d47NgQvf2n6O0v9kV78vgw4/9q3LRlQzmLVqFn1yZu397OmpVtCztoSZKkZcjQJi0Dmcmh589MmkF76Mgg5y9+bd1+d1eV/Z0VqutXL/CoJUmSBIY2aUk6OXyOBwcmtznWzjTU7e+slGUhm+juqrCzus66fUmSpBZlaJMWuZELozx0uF6Gszq9/ae+pm5/fAatp8u6fUmSpMXG0CYtIqNjyZPHhyeWOPb213j02KW6/R2VtRPhrLuryl07K2xcY92+JEnSYmZok1rYs4MjPHCoWN7Ye6jGwcN1hs9dBKB9zUr2l/uhjVfub+uwbl+SJGmpMbRJLWJo5AIHx5c59tfo669zbHAEKOr2b9/eMaksZM8W6/YlSZKWg5lCm2uqpDlyYXS8bv/SMscnGur2d1+3nq/fs7nYE62ryh3bO1i7yrp9SZIkTeZMm9QEmUn/82d5oP8UfeV+aA8drnOurNvfvGF1sRda2eTY01W1bl+SJEkTXB65wD71wGHef99jHKmdZUd1HW971V6+78DOhR6WXoBTp8/TO3BpBq2vv8apsm5/zcoVZd1+cQ9aT1eVzk3W7UuSJOnyDG0L6FMPHOanP3GQsxdGJ46tW9XGL71mn8FtkRi5MMrDRwaLe9AGipD2zMkzAETArdva6e6qsL8sCtl7QzurrNuXJEnSVTC0LaBves+fc7h29muO76yu43+949sWYESaydhY8tSJ4Yk2x77+Oo8cHeRiWbe/vbK2aHHcVZSF7Ou0bl+SJEkvnEUkC+jINIFtpuOaX88OjkwqCjk4UGeorNvfuGYl+zsr3P0teyaWOl5v3b4kSZLmmaFtju2orpt2pm1Hdd0CjGZ5Gz53kQfL2bPxkDZet79yRXDb9nb+bs8OerqqHNhVZc+WjdbtS5IkacE1NbRFxHrgN4DXlIc+DvxEZk47rRQRbwB+DLgdGAX+BvjnmXmwPL8b+CpwBhhf61jLzM5mjnsuve1Ve6e9p+1tr9q7gKNa+i6OjvHosaFyiWMR0L7y3KW6/RuvW8//ddPmYk+0rip37rBuX5IkSa2p2TNtvwbcBuylCFmfAj4A/OPLXN8O/CvgfwMXgX8JfDoibs7MMw3X7c3MgSaPdV6Ml43YHjl3MpOBU2cnLXN86EidkQtF3f6m9avo6arynfu2F8scO6ts2mDdviRJkhaHphWRRMQ64HnguzPzc+WxVwJ/AmzOzJFZvMda4Czwksy8v2GmretqQ1urFJGo+WpnztM3UKd3oiykxsnT54Gibv+unZWJGbSezipdm63blyRJUmubryKSvcBa4EsNx+4H1gG3Ag/O4j1eSbEU8itTjv91RKwGHgZ+NjP/x3Qvjoi7gbuvbthqZSMXRvny0bJuv5xFe7qhbv+WrRv5ttu2TRSFWLcvSZKkpWZWM20R8RHgjTNc8ovAfcD/BFaMT3FFxAqKe9Venpl/eYXPuBX4S+BfZuYHy2Mbgbsowt8q4B8C7wNelpkzhkBn2hafom7/9MQyx76BGo8cHeTCaPHP8fqONfR0VdnfWeVAV1G337521QKPWpIkSXrhXvA+bWV4mqnr/AzFbNoDwKbMrJWv20SxZLJ7ppAVEXcAnwF+LTPfd4WxfAb4q8z8F1e4ztDW4p4bGmlY4linb6DG0EhRt79hdRv7O6sTM2g9XVVuqFi3L0mSpKXpBS+PzMxhYPgKH/IYMAK8GPjz8vABinvUHp/hdS8G/gz4+cz89VkMZwzwBqVF5vS5ixw8XJ+YQes9VONIfUrdfvcOuruKWbQ9WzfSZt2+JEmS1LwiEoCI+C2K+v7xyv9PAA9l5o9d5vpvAv4rRc3/b01z/uspwuKjFAHzDcC/Bb4pM7809fopr3WmbYFcHB3j8WeHJ8JZ30CNx58dYqz8x7Fr8/qGGbQKd+6oWLcvSZKkZe0FL4+8ig/aAPw6l9mnLSLeCbwuM+8sn/8F8Lcollc2+o7M/HxEvBb4eWA7xSzel4F3Z+ZnZjEWQ9s8GK/bb9wP7eDhS3X71fWr6O68tMSxu6vKZuv2JUmSpEnmLbS1EkPb3KifuVDMoDWUhZwYLur2V69cwZ07Oi4FtM4qN1633rp9SZIk6QoMbbom5y6O8sjRIXoPnSr2Reuv8dUTpyfO37x1Az1dm+jpqtDTtYm9N7SzeqV1+5IkSdLVMrTpisbGkq+ePD1pP7QvN9Ttb2tfc2nD6rJuv8O6fUmSJKkpDG36GseHzk2Es/H70QYb6vb3dVaKgNZZpWdXlRs61rrMUZIkSZojhrZl7vS5izx0uN5wL1qdw7WzALStCG67oX1SQLvZun1JkiRpXhnalpGLo2N85bnhiaKQ3v7Jdftdm9dNanO8c0eFdaut25ckSZIWkqFticpMDtfO0tdfp7f/FH39dQ4ernP2wigAlXWryhm0Cj27quzvrLJl45oFHrUkSZKkqQxtS0T97AUebNiwure/zonhc8Cluv3xWbTuriq7rduXJEmSFgVD2yJ07uIojx4dKsLZoRq9AzWeOj593f7+ziq3b++wbl+SJElapAxtLS4z+eqJ02WLY50H+ms8cmSQ86NjAGxtX0N3Z5UDu6zblyRJkpYiQ1uLOTF8qW5/vDBkvG5//eo29u2sTBSFdHdV2V6xbl+SJElaygxtLeBPDx7lvx48Su+h2kTd/oqAvTd00NNVKe5F21XlRdvarduXJEmSlpmZQtvK+R7McnXwcJ3eQzV6dlV50zfupruryl07O1i/2n8EkiRJki7PmbZ5MjqWzqBJkiRJmtZMM23WDc4TA5skSZKka2FokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBbW1NAWEesj4sMRUSsf90bEuhmuf1NEjEXEcMPj96dc89KI+EJEnImIJyPi9c0csyRJkiS1smbPtP0acBuwF7gVuB34wBVe81Rmbmx4vHb8RERUgD8FPg5sAn4M+GBEfEOTxy1JkiRJLalpoa2cUXs98C8y89nMfA74F8AbI2LtNb7ta4AzwPsy81xmfgb4JHB3UwYtSZIkSS2umTNte4G1wJcajt0PrKOYdbucrog4FhH9EfEHEXFTw7lu4IHMzCnv2T3dG0XE3RHxxYj44rX9CJIkSZLUWmYV2iLiIxGRMzx+AWgvL683vHT89x2Xeev/CewDdgBfB4wAn4mIDeX59invB1C73Ptl5j2Z+dLMfOlsfi5JkiRJanUrZ3ndTwA/NcP5M1yaTatQBKvx3wMMTveizHyq4emxiPhRipD29cDngCFg95SXVS/3fpIkSZK01MwqtGXmMDA80zUR8RjFTNmLgT8vDx8AzgKPz3I8WT6ifN4HfN+Uaw6UxyVJkiRpyWvaPW2ZeRb4GPDuiNgWEduAdwMfzcyR6V4TEd8VEZ1R2Az8JnAC+Kvykk8CGyLibRGxOiJeSVFOck+zxi1JkiRJrazZlf9voZhVG388Brx1/GREvDMiHm64/luBL1DM4j0MXAd8ezmzR2bWgO8EfpBi2eRvAT+Wmf+nyeOWJEmSpJYUk4sZl46IyKX6s0mSJElaWiKCzIzpzjV7pk2SJEmS1ESGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamGGNkmSJElqYYY2SZIkSWphhjZJkiRJamFNDW0RsT4iPhwRtfJxb0Ssm+H6D0bE8JRHRsRPNlyTEXFmyjWVZo5bkiRJklpVZGbz3izit4A7ge8HEvgU0JeZ/3iWr/924L8DN2bmkfJYAi/PzL+8yrFkM382SZIkSZorEUFmxrTnmhVsyhm154HvzszPlcdeCfwJsDkzR2bxHn8ErMjM1zQcM7RJkiRJWtJmCm3NXB65F1gLfKnh2P3AOuDWK704Im4Avhf44DSn/3NEnIiIv46I10xzfvw97o6IL0bEF69u6JIkSZLUmmYV2iLiI+W9ZZd7/ALQXl5eb3jp+O87ZvExPwIcAj4z5fjfBm4COoEPAL8XEa+e7g0y857MfGlmvnQ2P5ckSZIktbpZLY+MiI0Us2iXc4ZiNu0BYFNm1srXbaJYMtmdmQ/O8P4rgKeAf5+Z773CWH4LWJuZP3SF61weKUmSJGlRmGl55MrZvEFmDgPDV/iQx4AR4MXAn5eHDwBngcev8BGvBrYDH57FcMaAaX8YSZIkSVpqZhXaZiMzz0bEx4B3R8RD5eF3Ax+dRQnJm4FPZObxxoMRcRewHuilaKP8LuCHgH/QrHFLkiRJUitr9ubab6GYVRt/PAa8dfxkRLwzIh5ufEFE7KQIY9MVkGwFfhs4BTwHvAv4h5n5X5o8bkmSJElqSU3dp62VeE+bJEmSpMVivir/JUmSJElNZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQWZmiTJEmSpBZmaJMkSZKkFmZokyRJkqQW1tTQFhH/b0T8dUSciYgnZvmaV0fEwxFxNiIeioi/M+X8LRHx2Yg4HREDEfHPmjlmSZIkSWplzZ5pOwK8D/jF2VwcEXuATwC/BFTKXz8ZEbvL823AnwCPAFuBvwu8PSL+fpPHLUmSJEktqamhLTP/KDM/Dhye5UveCHwpMz+Wmecz8/eA+8vjAN8C3Aj8dGaeycz7gQ8BP9bMcUuSJElSq1q5wJ/fDXxpyrH7y+Pj5x/PzOEp5398ujeLiLuBuxueN2+kkiRJkrQAZhXaIuIjXJr9ms4vZua7ruHz24H6lGM14M4rnO+Y7s0y8x7gnmsYh65BRGRmmoy1IPz+aaH43dNC8vunheJ3b2HNdqbtJ4CfmuH8mWv8/CGKe9kaVYHBWZ6XJEmSpCVtVqGtXJ44fMULr14f8Iopxw4An2s4f2tEbMjM0w3n++ZgLJIkSZLUcppd+b8yItYCq4qnsbZ8fjkfBV4aEa+NiFUR8VrgJcDvlOf/J/AM8K8jYl1E9ABvpigj0cL7uYUegJY1v39aKH73tJD8/mmh+N1bQJGZzXuziJ8F/tXU4+PrXyPidcCHMnNjw2teDfwKsAd4CnhrZn664fwtFCHtGyjuZ/tAZv5y0wYtSZIkSS2sqaFNkiRJktRczd5cW5IkSZLURIY2SZIkSWphhjbNWkSsj4gPR0StfNwbEetmuP4NEfG/I+JURJyIiD+NiH3zOWYtDdfw3esuv2/HIiIj4pvnc7xa3CKiLSLeHxHHI2IoIj4eEVtmuP7VEfFwRJyNiIci4u/M53i1dFzNdy8idkbEH0fEM+V/514/3+PV0nKV37/vjIg/L/98dyoiPh8RL5/vMS8nhjZdjV8DbgP2ArcCtwMfmOH6dopimk5gJ3A/8OmIWD/H49TSc7XfvfPAJ4DvnvuhaQl6B/C9wMso/vsF8LvTXRgReyi+a79Esa/oLwGfjIjdcz9MLUGz/u4BY8Cngf8bGJj7oWkZuJrv3ybg14FbgK3AfwT+NCK65nqQy5VFJJqVclbjeeC7M/Nz5bFXAn8CbM7MkVm8x1rgLPCSzLx/LserpeOFfvciIoGXZ+ZfzvlgtSRExDPAuzPz3vL5zcATwO7MfGbKtT8HfFtmvrzh2OeBz2am9di6Klfz3ZvyuqeBd2Xmx+ZloFqSrvX71/D6Y8A/ycxPzO1Ilydn2jRbe4G1wJcajt0PrKOY+ZiNVwJngK80d2ha4prx3ZNmJSKqwC4avm+Z+SQwCHRP85JuJn83ofh+TnetdFnX8N2TmuaFfv/K21+2AAfnaIjLnqFNRMRHyvXwl3v8AsVSR4B6w0vHf98xi8+4Ffht4J9l5lBzfwItVvPx3ZOu0nTfNyj2CZ3u+9Z+FddKM7na757UTNf8/YuIbcDHgV/OTP9ifo6sXOgBqCX8BPBTM5w/w6UZjQrFv8Djv4fib2EuKyLuAD5D8S/zB699mFqC5vS7J12D8b9Uqkw5XmX679vQVVwrzeRqv3tSM13T9y8idlD8Ge/TwE/PycgEONMmIDOHM/PEDI8zwGPACPDihpceoLhH7fHLvXdEvBj4H8B7MvN9c/hjaBGay++edC0yswYcouH7VpaNdAAPTvOSPiZ/N6H4fvbN0RC1RF3Dd09qmmv5/pWFS58H/jQzfyItyphThjbNSmaeBT4GvDsitpVT4e8GPnq5IoiI+Cbgc8DPZOavz99otZRc43cvImJtWX4DsLp83jZPw9bidg/w9oi4KSI6gPcC92Xm09Nc+1HgpRHx2ohYFRGvBV4C/M78DVdLyNV892j471wAq8rnrqLStZr19y8ibgP+Evj9zJxpxYyaxNCmq/EWipmN8cdjwFvHT0bEOyPi4Ybrf4Fimv1XI2K44eE+Hrpab+Hqvns3UszEnS2ff678/Q/Nx2C16L2Hop30b4DDQBvweoCIeF1EDI9fWN6o/xrgXRRLiN4FfP/l/pAtXcGsv3ul8f/O7QI+XP7+XfM2Wi01V/P9ezvFdk5vmfJnvNfN96CXCyv/JUmSJKmFOdMmSZIkSS3M0CZJkiRJLczQJkmSJEktzNAmSZIkSS3M0CZJkiRJLczQJkmSJEktzNAmSZIkSS3M0CZJkiRJLez/B1849Nra1WvGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use neutral structure and datapoints of disp_002,004,006\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from nonrad.ccd import get_Q_from_struct\n",
    "from nonrad.elphon import get_Wif_from_WSWQ\n",
    "from pathlib import Path\n",
    "from pymatgen.core.structure import Structure\n",
    "from glob import glob\n",
    "\n",
    "cc_dir = \"/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad\" + \\\n",
    "         \"/\" + \\\n",
    "         \"/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited\"\n",
    "\n",
    "SnZn_4 = Structure.from_file(f\"{cc_dir}/4/CONTCAR\")\n",
    "SnZn_5 = Structure.from_file(f\"{cc_dir}/5/CONTCAR\")\n",
    "SnZn_6 = Structure.from_file(f\"{cc_dir}/6/CONTCAR\")\n",
    "\n",
    "\n",
    "WSWQs = []\n",
    "for d in glob(f\"{cc_dir}/*\"):\n",
    "        if d[-2:] in [\"/4\", \"/5\", \"/6\"]:\n",
    "            pd = Path(d)\n",
    "            Q = get_Q_from_struct(SnZn_5, SnZn_6, f\"{pd}/CONTCAR\")\n",
    "            path_wswq = f\"{pd}/eph/WSWQ\"\n",
    "            WSWQs.append((Q, path_wswq))\n",
    "print(WSWQs)      \n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "ground_vr_path = f\"{cc_dir}/5/vasprun.xml\" \n",
    "\n",
    "Wifs = get_Wif_from_WSWQ(WSWQs, ground_vr_path, 2305, 2302, spin=1, fig=fig)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.ylim(-1e-6,1e-6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "40cab6a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/6/eph/WSWQ'), (-0.0, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/5/eph/WSWQ'), (-0.23268478223937827, '/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited/4/eph/WSWQ')]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAFgCAYAAACmKdhBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA94ElEQVR4nO3deXRc53nn+e+LHSQBgjuIutopUiIlLlVeY1leZFuLrYUsps94knYynRmn00kn6XTcSTrJpO2kE6czSU9P+kx33HHGcdLLZFikFlubJXmRl9hOFTdROyVRugWA4A4uALHUO38AlGmKuwAUlu/nHByibr1V+OH4iuZTz3PfG2KMSJIkSZLevppqB5AkSZKk6cICS5IkSZLGiAWWJEmSJI0RCyxJkiRJGiMWWJIkSZI0RiywJEmSJGmMWGBJkiRJ0hixwJIkTSshhMYQwhdDCHtCCEdDCNtCCHeOPrcyhPAPIYRDo19PhBBWnvbaEEL44xDCgdGvPw4hhNOerw0h/EEIoXP0vbeGENqq8GtKkiapumoHkCRpjNUBbwAfAF4H7gL+LoRwM9AJbAT2MPIh4y8C/wNYPfraTwP3AWuACHwNeBX4z6PPfxb4CeC9o++9Cugf719IkjR1hBhjtTNIkjSuQgg7gM/GGAunHasDfh74kxjjrNFj3wW+FGP8wujjnwP+txjje0II8xgp3NbEGHdP+C8hSZoSHBGUJE1rIYQlwHJg12nHDjPSefpz4A9PW74K2H7a4+2jxwBuBoaAjSGE7hDCiyGEXxzH6JKkKcgRQUnStBVCqAf+K/DXMcbnTx2PMbaFEGYDP8PIuOApc4Ajpz0+AswZvQ4rAeYyUqxdA1wPPBlCeDHG+LXx/U0kSVOFHSxJ0rQUQqgB/gYYAH7pzOdjjMcZubbqyyGExaOHjwGtpy1rBY7FkXn6vtFjn4sx9sUYdzBy/dZd4/QrSJKmIAssSdK0M9px+iKwBMjHGAfPsbQGmAVkRh/vYmSDi1PW8KPRwh2jf55+8bIXMkuSfowFliRpOvpPwI3A3THGU50nQggfDSGsG91uvRX4M+AQ8Nzoki8DvxZCyIQQOoB/CXwJYHRji6eB3x7dCv5G4H8CvjJRv5QkafLzGixJ0rQSQriKkd0BTwLdp93G6ucZGRf8c0aup+oDfgDcEWM8tdX6XwDXAjtHH//l6LFTPslIZ+wA0AP8bozxyXH7ZSRJU47btEuSJEnSGHFEUJIkSZLGiAWWJEmSJI0RCyxJkiRJGiMWWJIkSZI0RqbtLoIhBHfvkCRJkjRuYozhzGPTtsACcIdESZIkSePhtNuA/BhHBCVJkiRpjFhgSZIkSdIYscCSJEmSpDFigSVJkiRJY8QCS5IkSZLGiAWWJEmSJI0RCyxJkiRJGiMWWJIkSZI0RiywJEmSJGmMWGBJkiRJ0hixwJIkSZKkMWKBJUmSJGlSGhiq8Piubn7xv5U4fnKo2nEuSl21A0iSJEnSKTFGdnX2sqmY8uD2Tg4eH2DhnEZ27zvG6qSt2vEuyAJLkiRJUtX19PZz/7YyhWKZF/YepaG2ho+uXEI+l+HW6xdRVzs1hu9CjLHaGcZFCCFO199NkiRJmg76B4d54rm9FIop33ppP8OVyLor28hnE+5e3cHcWfXVjnhOIQRijOEtx6drEWKBJUmSJE0+MUZKrx+mUEr5yvZOevuHWDq3ifXrMuRzCdctmvPm2vu3lvmTx16g83AfHW3NfOb2Fdy3LlPF9D9igSVJkiSpajoP97Fla5lCMeWV/cdpqq/hzpuWks8mvPe6BdTW/Hitcv/WMr+1eSd9g8NvHmuur+WPNtw8KYosCyxJkiRJE+rEwBCP7Oxm89aU7+4+QIzwrmvmszGXcNfNS5nTeO4tId73+acoH+57y/FMWzPf+c0Pj2fsi3KuAstNLiRJkiSNmUol8v1XD7KpmPLIM12cGBjmyvmz+JXbrmfDuoQrF8y6qPfpPEtxdb7jk4UFliRJkqS37bX9x9lcSimUypQP9zGnsY67V3eQzyW88+p5hPCWZs95dbQ1n7WD1dHWPFaRx4UFliRJkqTL0ts/yFd3dFEopvzDnkOEALcsW8i/umMFH1vZTnND7WW/92duX3HWa7A+c/uKsYg+bsb8GqwQQi3weeBngSbgceDnY4z7z7H+DuBPgWuB3cCvxRgfH31uOfCHwHuBVuB14N/HGP/yInJ4DZYkSZI0xoYrkadf2kehVObxXd2cHKpw3aLZ5HMJ69dlWDp37DpM7iI48oN+G/gZ4A7gAPBXwKwY451nWXst8AzwaeDvgJ8EvgCsijG+FkJ4N/AOYAvQBbwP+ArwT2KMmy+QwwJLkiRJGiMv7j1KoZiyZWuZnqMnmdtczz1rRkYA1yRzL3kEcKqbyAJrD/C5GOMXRx9fB7wMXB1j3HPG2s8CH44xvv+0Y08DT8QYP3uO9/8fQE+M8ZcvkMMCS5IkSXobDh0f4MHtnRRKKTvSI9TWBD60YhH5bMKHb1xMY93ljwBOdROyi2AIoQ24EiieOhZj3B1C6AXWAHvOeMma09eOKo0eP9v7zwLeA/z+GEWWJEmSdJrB4Qpff76HQinlqed7GByOrFzayu9+YiX3ru1g4ZzGakec1MZ6k4uW0T+PnHH8MCPXUJ1t/dnWrjpz4ei1XX8DvAF8+Ww/PITwaUbGDSVJkiRdpBgjuzp7KZRSHtzWyYHjAyyc08Cn3ns1+WzCyo6z/VNeZzPWBdbR0T/nnnG8Deg9x/oLrg0h1AP/FVgK3BljHDzbD48xfoGRa7gIITgfKEmSJJ1Hz9F+Htg6MgL4fPdRGmpr+MjKxeSzCbcuX0R9bU21I045Y1pgxRgPhxBeB7LANnhzI4tWYMdZXrId+NAZx9YBT556EEJoAjYBc4CPxRiPjWVmSZIkaSbpHxzmyedGRgC/+eI+hiuRtVe08fv33cTdq5fSNquh2hGntPHaRfBT/GgXwS8CLTHGO86y9jpgJ/BzjBRRG4G/5Ee7CM4BHgIGgXtjjBd922Y3uZAkSZJGxBgpvX6YQinlK9s76e0for21iQ3ZDBuyCcsWz6l2xClnQja5GPV5YB7wQ6AR+Brw06Mhfgr4ixjjHHhzA4wNjNwH66+AV4D1McbXRt8rD3wQ6AP2nbb149/GGP/pOGSXJEmSpo3Ow31s2VqmUEx5Zf9xmupruGNVO/lcwk9ct5Dampm1tfpEGPMO1mRhB0uSJEkz0YmBIR59pptCKeW7uw8QI7zrmvlszCbceXM7LU311Y44LUzYfbAmCwssSZIkzRSVSuT7rx6kUEp5ZGcXxweGuXL+rJERwHUJVy6YVe2I085EjghKkiRJmgCv7T/O5lLK5q1l0kN9zGms4xOrO8jnEt559TxOu8RGE8QOliRJkjSF9PYP8vCOLgqllB++dogQ4JZlC8lnE25f1U5zQ221I84IjghKkiRJU9RwJfLtl/dTKKY8tqubk0MVrls0m3wuYf26DEvnNlc74ozjiKAkSZI0xbzcc5RNxTJbtqbs7T3J3OZ6/tE7riCfS1iTzHUEcBKygyVJkiRNIoeOD/DQjk4KxZTt6RFqawIfXL6IfC7hthsX01jnCOBk4IigJEmSNEkNDlf4xgv7KBRTnnx+L4PDkRuXtpLPZrh3bYZFLY3VjqgzWGBJkiRJk0iMkV2dvRRKKQ9u6+TA8QEWzmng3rUZ8tmElR2t1Y6o8/AaLEmSJGkS6DnazwNbOymUUp7vPkpDbQ0fWbmYfDbh1uWLqK+tqXZEvQ12sCRJkqRx1j84zBPP7aVQTPnWS/sZrkTWXNHGxmyGu9d00DarodoRdYkcEZQkSZImUIyRrW8cplBMeWh7J739Q7S3NrE+OzICuGzxnGpH1NvgiKAkSZI0AToP97Fla5lCMeWV/cdpqq/hjlXt5HMJP3HdQmpr3Fp9OrODJUmSJL1NJwaGePSZbgqllO/uPkCM8K5r5rMxm3Dnze20NNVXO6LGmCOCkiRJ0hiqVCI/eO0ghWLKwzu7OD4wzBXzm9mwLiGfTbhywaxqR9Q4ckRQkiRJGgN7DhynUCqzuZSSHupjTmMdH1+9lHw24Z1Xz6fGEcAZzQ6WJEmSdAG9/YM8srOLTcWUH752iBDglmULyWcTbl/VTnNDbbUjaoI5IihJkiRdguFK5Dsv76dQSnn0mW5ODlW4dtFs8tmE9esydLQ1VzuiqsgRQUmSJOkivNxzlE3FMvdvLdPd28/c5np+8h0j11WtvaKNEBwB1LnZwZIkSdKMd+j4AF/Z0cmmUpntbxymtibwgeWL2JhLuO3GxTTWOQKoH+eIoCRJknSaweEK33xhH5uKKU8+v5fB4cgN7S1szCXcs7aDxS1N1Y6oScwCS5IkSQJ2dR6hUCzzwLYyB44PsGB2A/euzZDPZVjVMbfa8TRFeA2WJEmSZqx9R0/ywLYym4opz3cfpaG2httuXEw+m/CBFYuor62pdkRNE3awJEmSNC31Dw7z5HM9FEop33xxH8OVyJor2tiYzXD3mg7aZjVUO6KmMEcEJUmSNO3FGNn2xmEKpZQHt3XS2z/EktZG1q9L2JjLsGxxS7UjappwRFCSJEnTVteRPjaXyhRKKa/sO05TfQ23r2onn01437KF1Na4tbomhh0sSZIkTUl9A8M8uquLQrHMd3bvJ0Z419Xzyecy3HXzUlqa6qsdUdOYI4KSJEma8iqVyA9fO0ihlPLwzm6OnRwimdfMhmxCPpvhqgWzqx1RM4QjgpIkSZqyXj9wgkIpZfPWlDcO9jG7oZa7bl5KPpfwrqvnU+MIoCYJO1iSJEmalI72D/LwzpERwB+8dpAQ4H3XLSSfy3D7qnZmNdgrUPU4IihJkqRJb7gS+e7u/RSKKY/u6qZ/sMK1i2aTzyasX5eho6252hElwBFBSZIkTWIv9xyjUErZUirT3dtPa1Md+WzCxlzC2ivaCMERQE0NdrAkSZJUFYdPDPDQ9k42lcpsf+MwtTWBDyxfRD6bcNuNi2mqr612ROmcHBGUJElS1Q0OV/jmC/solFKefK6HgeEKN7S3sDGXcO/aDItaGqsdUbooFliSJEmqmmc7e9lUTHlwe5n9xwZYMLuBe9dmyOcyrOqYW+140iXzGixJkiRNqH1HT/LAtjKFUpnnunppqK3hthsXk88mfGDFIupra6odURpzdrAkSZI0Zk4ODfPkcz0UiinfeHEfw5XImmQu+VzC3as7mDe7odoRpTHhiKAkSZLGRYyR7ekRNhXf4KHtXRzpG2RJayPr1yXksxmuX9JS7YjSmHNEUJIkSWOq60gfW7aWKRRTdu87TmNdDbevaiefS7hl2UJqa9xaXTOPHSxJkiRdtL6BYR7b1U2hlPLtl/cTI7zz6nnkswl3rV5Ka1N9tSNKE8IRQUmSJF2WSiXyg9cOUiimPLyzi+MDwyTzmtmQHRkBvGrB7GpHlCacI4KSJEm6JK8fOEGhlLJ5a8obB/uY3VDLXTcvJZ9LeNfV86lxBFB6CztYkiRJetPR/kEe3tlFoVjmB68dJAR433ULyecy3L6qnVkNfj4vgSOCkiRJOofhSuS7u/ezqZjy2K5u+gcrXLtwNvlcwvp1GTramqsdUZp0HBGUJEnSj3m55xiFUsqWUpnu3n5am+rIZxPyuYR1V7QRgiOA0qWygyVJkjSDHD4xwEPbO9lUKrP9jcPU1gRuvX4h+VzCR25cQlN9bbUjSlOCI4KSJEkz1OBwhW+9uI9CKeWJZ3sYGK5wQ3sL+WzCves6WNzSVO2I0pRjgSVJkjTDPNvZS6GU8sC2MvuPDTB/dgP3ru0gn01Y1dHqCKD0NngNliRJ0gyw/9hJHtjWyaZiynNdvdTXBm67YQn5XMIHVyyivram2hGlaW3MO1ghhFrg88DPAk3A48DPxxj3n2P9HcCfAtcCu4FfizE+ftrzfwm8F1gBfCnG+L9eZA47WJIkaUY4OTTMU8/1UCilfOOFfQxVIquTueSzCfes6WDe7IZqR5SmnYnsYP0mcC/wbuAA8FfA3wB3niXUtcBm4NPA3wE/CWwJIayKMb42umwH8P8BPz8OWSVJkqakGCPb0yMUiikPbu/kSN8gi1sa+blbrmFjLuH6JS3VjijNSOPRwdoDfC7G+MXRx9cBLwNXxxj3nLH2s8CHY4zvP+3Y08ATMcbPnrH2S8CQHSxJkjSTdR3pY8vWMoViyu59x2msq+H2Ve3kcwm3LFtIbY3XVUkTYUI6WCGENuBKoHjqWIxxdwihF1gD7DnjJWtOXzuqNHr8cn7+pxnphkmSJE0bfQPDPLarm0Ip5dsv7ydGeMdV8/ijDdfy8dVLaW2qr3ZESaPGekTwVC/6yBnHDwOt51h/trWrLueHxxi/AHwBRjpYl/MekiRJk0GMkR+8epBCKeXhnd0cOzlEpq2Zf/6hZWzIJly9cHa1I0o6i7EusI6O/jn3jONtQO851l/sWkmSpGnvjYMnKJRSNpfKvH7wBLMaarnr5qXkswnvvmY+NY4ASpPamBZYMcbDIYTXgSywDd7cyKKVkc0qzrQd+NAZx9YBT45lLkmSpMns2MkhHt7RxaZSyg9ePUgI8BPXLeBXP3I9d9zUzqwG76wjTRXj8V/rF4DfCCF8nZFdBP8YeOy0XQFP92XgMyGETwKbgI1ADvjUqQUhhAagBqgFYgihCajEGAfGIbskSdKEGK5Evrt7P4ViyqO7uukfrHDNwtn8+seWsz6bkGlrrnZESZdhvO6D9ceM3AerEfga8OkY4/4Qwk8BfxFjnHPa+tPvg/UK8C/OuA/WN4APnPFjvhlj/OAFcriLoCRJmnR27ztGoZiyZWuZriP9tDTVcfeaDvLZhOyVbYTgCKA0FZxrF8ExL7AmCwssSZI0WRw5MciDOzopFFO2vXGYmgC3Ll/ExlzCR25cQlN9bbUjSrpEFliSJEkTaGi4wrde2kehWOZrz+5lYLjCiiUtbMwl3Luug8UtTdWOKOltsMCSJEmaAM919VIopty/rZP9x04yf3YD96zpYGMuYVVHqyOA0jQxITcaliRJmon2HzvJA9tGRgCf7eqlvjbw4RsWk88mfHDFYhrqaqodUdIEsYMlSZJ0GU4ODfPUcz0USinfeGEfQ5XI6mQu+WzC3Ws6mD+7odoRJY0jRwQlSZLephgj29MjFIopD+3o5PCJQZa0NnLfugz5bMLyJS3VjihpgjgiKEmSdJm6j/SzZWuZQinl5Z5jNNbVcPuqdvK5hFuWLaS2xuuqJI2wgyVJknQWfQPDPP5sN5uKKd95eT+VCO+4ah75XMLHVy+ltam+2hElVZEjgpIkSRcQY+SHrx2iUEz56s4ujp0cItPWTD6bYUM24eqFs6sdUdIk4YigJEnSObxx8ASFUsrmUpnXD55gVkMtd960lI25hHdfM58aRwAlXSQ7WJIkaUY6dnKIh3d2USimfP/Vg4QA7712Aflswh03tTO70c+hJZ2bI4KSJGnGG65Evrf7AIVSyqPPdNM3OMzVC2axMZewPpuQaWuudkRJU4QjgpIkacZ6Zd+xN0cAu47009JUx33rMmzMZcheOY8QHAGUNDbsYEmSpGnpyIlBHtrRSaGUsvX1w9QEuHX5IvLZhI+uXEJTfW21I0qawhwRlCRJ097QcIVvvbSPQrHM157by8BQheVL5pDPJqxfl2Fxa1O1I0qaJiywJEnStPVcVy+FYsr92zrZf+wk82c3cM+aDjbmElZ1tDoCKGnMeQ2WJEmaVvYfO8kD2zopFFOe7eqlvjbw4RsWk88mfHDFYhrqaqodUdIMZAdLkiRNGSeHhvn68z1sKpb5xgs9DFUiN2fmks9muGdthvmzG6odUdIM4YigJEmakmKM7EiPUCilPLi9k8MnBlnc0sj6dRnyuYTlS1qqHVHSDOSIoCRJmlK6j/SzZWuZQinl5Z5jNNTV8LGVS9iYS7hl2ULqah0BlDT52MGSJEmTRv/gMI/t6qZQKvPtl/ZRiZC7ah75bMLHVy9lbnN9tSNKEuCIoCRJmqRijPzDnkMUiilf3dHF0ZNDZNqa2ZDNsCGbcM3C2dWOKElv4YigJEmaVN44eOLNEcA9B04wq6GWO25qZ2Mu4T3XLKCmxq3VJU09drAkSdKEOX5yiId3dlEopfz9KwcBeO+1C8jnEu68qZ3ZjX72K2lqcERQkiRVRaUS+d4rBygUUx55ppu+wWGuWTibDesyrM9mSObNqnZESbpkjghKkqQJ9cq+YxRKKVtKZTqP9NPSVMd96zJszGXIXjmPEBwBlDT92MGSJElj5siJQR7a0UmhlLL19cPUBLh1+SLy2YSPrlxCU31ttSNK0phwRFCSJI2LoeEK33ppH4Vima89t5eBoQrLl8xhYy7hvrUZFrc2VTuiJI05CyxJkjSmnu/upVBMuX9bJ/uOnmTerHruXZthYy5hVUerI4CSpjWvwZIkSW/bgWMneWDbyAjgrs5e6moCH75hMflcwodWLKahrqbaESWpquxgSZKk8xoYqvDU8z0USilff76HoUrkpkwrG7MJ96zNMH92Q7UjStKEc0RQkiRdtBgjO8tHKBRTHtzeyaETgyxqaWT9ugz5bMKK9pZqR5SkqnJEUJIkXdDe3n62bC1TKKa81HOMhroaPrZyCflcwvuXLaSu1hFASTofO1iSJM1w/YPDPLarm0KpzLdf2kclQu6qeWzIZvjE6g7mNtdXO6IkTTqOCEqSpDfFGPmHPYcoFFO+uqOLoyeHyLQ1syGbYUM24ZqFs6sdUZImNUcEJUkSbxw8weZSmc1bU/YcOMGshlruvGkp+VyG91yzgJoat1aXpLfDDpYkSdPcsZNDPLKzi03FlO+/ehCA9167gHwu4c6b2pnd6OetknSpHBGUJGkGqVQi33vlAIViyiPPdNM3OMzVC2aRzyasz2ZI5s2qdkRJmtIcEZQkaQZ4Zd8xCqWULaUynUf6aWms4751HWzMJWSvnEcIjgBK0niygyVJ0hR3pG+Qr+zopFBMKb1+mJoA779+EflcwsdWLqGpvrbaESVp2nFEUJKkaWRouMLTL+1nUynla8/uZWCowvIlc8hnE+5bl2FJa1O1I0rStGaBJUnSNPB8dy+FYsr92zrZd/Qk82bVc8+aDvK5hJszcx0BlKQJ4jVYkiRNUQePD/DAtjKFUsoz5V7qagIfumEx+WzCh29YTENdTbUjSpJG2cGSJGkSGhiq8NTzPRRKKV9/voehSuSmTCv5bMI9azpYMKex2hElaUZzRFCSpEkuxsjO8hEKxZQHt3dy6MQgC+c0sn7dyAjgDe2t1Y4oSRrliKAkSZPU3t5+tmwtUyimvNRzjIa6Gj66cgkbswnvv34hdbWOAErSVGEHS5KkKugfHOaxXd0USmW+/dI+KhGyV7aRzyV84uYO5s6qr3ZESdJ5OCIoSVKVxRgp7jlEoZTyle1dHD05RMfcJjZkEzZkM1y7aE61I0qSLpIjgpIkVUl66ASbS2U2l1JeO3CC5vpa7ry5nY3ZhPdcu4CaGrdWl6TpYsw7WCGEWuDzwM8CTcDjwM/HGPefY/0dwJ8C1wK7gV+LMT5+2vPLgP8MvBc4BPz7GOOfXkQOO1iSpAv6nft38t+//wbDMVIbAp989xX8wX03v+33PX5yiEee6aZQTPneKwcAeM+188lnE+68eSlzGv2MU5KmsgkbEQwh/DbwM8AdwAHgr4BZMcY7z7L2WuAZ4NPA3wE/CXwBWBVjfG20WHsGeAL4DeAG4FHgn8cY/98L5LDAkiSd1+/cv5O//fvX33L8p99z5WUVWZVK5O9fOcCmUsqjz3RzYmCYqxbMIp9NWL8uwxXzZ41FbEnSJDCRBdYe4HMxxi+OPr4OeBm4Osa454y1nwU+HGN8/2nHngaeiDF+NoTwIeCrwOIY47HR538fuCXG+KEL5LDAkiSd13W/9TDDZ/n/itoQ2P1Hd130+7y6/ziFYsqWrWXKh/toaazj46uXsjGXkLtqHiE4AihJ082EXIMVQmgDrgSKp47FGHeHEHqBNcCeM16y5vS1o0qjx089/+Kp4uq053/xHD//04x0wyRJuqCzFVfnO366I32DfHVHF4VSSnHPIWoC3HL9Iv7VHSu4fVU7TfW1Yx1XkjQFjPUAeMvon0fOOH4YONvdEVvOsXbVBZ4/650WY4xfYGTEkBCC7StJ0nnVhnDODtbZDA1XePrl/RSKKY8/u5eBoQrXL57Db955A+vXZVjS2jTekSVJk9xYF1hHR/+ce8bxNqD3HOvPt/ZCz0uSdNk++e4rznoN1ifffcWPPX6h+yiF0sgI4L6jJ2mbVc8n33kF+VzCzZm5jgBKkt40pgVWjPFwCOF1IAtsgzc3smgFdpzlJduBM6+lWgc8edrzy0MIs2OMx097fvtY5pYkzUynNrI42y6CB48P8MC2MoVSyjPlXupqAh9csZiNuQwfumExjXWOAEqS3mq8dhH8FD/aRfCLQEuM8Y6zrL0O2An8HLAJ2Aj8JW/dRfBx4DeBFcBjwK/EGP/HBXK4yYUk6ZIMDFV46vkeCqWUrz/fw1AlsqqjlXw24d61HSyY01jtiJKkSWIibzT8eWAe8EOgEfga8NOjIX4K+IsY4xx4cwOMDYzcB+uvgFeA9THG10afHw4h3A38BSPF2mHgTy5UXEmSdLFijDxT7mVT8Q0e3N7JoRODLJzTyP/yvqvJ5xJuaD/rZb+SJJ3VmHewJgs7WJKk8+np7WfL1pERwBf3HqOhroaPrlzCxmzC+69fSF1tTbUjSpImsQm7D9ZkYYElSTpT/+Awjz+7l0Ix5emX9lGJkL2yjXwu4RM3dzB3Vn21I0qSpggLLEnSjBRjpLjnEIVSyld2dHG0f4iOuU2sz2bYkE24btGcakeUJE1BE3kNliRJVZceOsGWUpnNW8u8uv84zfW13HlTO/lcwnuvXUBNjVurS5LGnh0sSdK0cfzkEI88002hmPK9Vw4A8J5r55PPJtx581LmNPq5oiRpbDgiKEmaliqVyN+/eoBCscwjz3RxYmCYK+fPIp9N2JDNcMX8WdWOKEmahhwRlCRNK6/tP06hlLK5VKZ8uI+WxjruWdNBPpfwjqvmEYIjgJKkiWcHS5I0ZRzpG+SrO7oolFKKew4RAtyybCEbcwkfW9lOc0NttSNKkmYIRwQlSVPS0HCFp1/eT6GY8vizexkYqrBs8Rzy2YT16zK0z22qdkRJ0gxkgSVJmlJe6D5KoZRy/9YyPUdP0jarfmQEMJuwOpnrCKAkqaq8BkuSNOkdPD7Ag9vKbCqlPFPupa4m8MEVi9mYy/ChGxbTWOcIoCRpcrODJUmqqoGhCl9/oYdCMeXrL/QwOBxZ1dFKPptwz9oOFs5prHZESZLewhFBSdKkEWNkV2cvm4opD27v5ODxARbOaWT9upFdAG9ob612REmSzssRQUlS1fX09nP/tjKFYpkX9h6lobaGj65cQj6X4dbrF1FXW1PtiJIkvS12sCRJ46p/cJivPbuXQinlWy/uoxJh3ZVt5LMJd6/uYO6s+mpHlCTpkjkiKEmaMDFGSq8fZlMx5Ss7OjnaP8TSuU2sX5chn0u4btGcakeUJOltcURQkjTuyof72FJKKZTKvLr/OE31Ndx501Ly2YT3XreA2hq3VpckTW92sCRJb8uJgSEe2dlNoZTyvVcOECO8+5r55HMJd928lDmNfpYnSZp+HBGUJI2ZSiXy968eoFAs88gzXZwYGObK+bPYkM2QzyZcMX9WtSNKkjSuHBGUJL1tr+0/TqGUsrlUpny4jzmNddy9emRr9XdePY8QHAGUJM1sdrAkSed1pG+Qr+7oolBKKe45RE2A9y1byMZcwsdWttPcUFvtiJIkTThHBCVJF224Enn6pX0USmUe39XNyaEKyxbPIZ9NWL8uQ/vcpmpHlCSpqiywJEkX9NLeo2wqpmzZWqbn6EnaZtVzz5oO8tmE1clcRwAlSRrlNViSpLM6dHyAB7d3Uiil7EiPUFcT+OCKRWzMJXzohsU01jkCKEnSxbKDJUkz0OBwha8/30OhlPLU8z0MDkdWLm0ln0u4d20HC+c0VjuiJEmTmiOCkjTDxRjZ1dlLoZTywLZODh4fYOGcRu5bO7IL4I1LW6sdUZKkKcMRQUmaoXp6+7l/W5lCscwLe4/SUFvDR1YuZmMu4dbrF1FXW1PtiJIkTRt2sCRpGuofHOZrz+6lUEr51ov7qERYe0Ub+VzC3auX0jarodoRJUma0hwRlKRpLsZI6fVDbCqW+cqOTo72D7F0bhPr12XI5xKuWzSn2hElSZo2HBGUpGmqfLiPLaWUQqnMq/uP01xfyx03tZPPJrz3ugXU1ri1uiRJE8UOliRNQScGhnhkZzeFUsr3XjlAjPDua+aTzyXcdfNS5jT6+ZkkSePJEUFJmuIqlcj3Xz3IpmLKI890cWJgmCvnzyKfTdiQzXDF/FnVjihJ0ozhiKAkTVGv7T/O5tERwPLhPuY01nH36pGt1d959TxCcARQkqTJwg6WJE1Cvf2DfHVHF4Viyj/sOUQIcMuyhWzMJXxsZTvNDbXVjihJ0ozmiKAkTXLDlci3X95PoZjy2K5uTg5VuG7RbPK5hPXrMiyd21ztiJIkaZQjgpI0Sb209yibSin3by2zt/ckc5vr+UfvuIJ8LmFNMtcRQEmSphA7WJJUBYeOD/Dg9k4KpZQd6RFqawIfWrGIfDbhwzcuprHOEUBJkiYzRwQlqcoGhyt8/fkeCqWUp57vYXA4cuPSVvLZDPeuzbCopbHaESVJ0kWywJKkKogxsquzl03FlAe3d3Lw+AAL5zRw79oM+WzCyo7WakeUJEmXwWuwJGkC9Rzt54GtIyOAz3cfpaG2ho+sXEw+m3Dr8kXU19ZUO6IkSRoHdrAkaYz0Dw7zxHN7KRRTvvXSfoYrkbVXtJHPJdy9eiltsxqqHVGSJI0RRwQlaRzEGNn6xmE2FVO+sr2T3v4h2lub2JDNsCGbsGzxnGpHlCRJ48ARQUkaQ52H+9iytUyhmPLK/uM01ddwx6p28rmEn7huIbU1bq0uSdJMZAdLki7SiYEhHn2mm0Ip5bu7DxAjvOua+WzMJtx5czstTfXVjihJkiaII4KSdBkqlcgPXjtIoZjy8M4ujg8Mc8X8ZvLZhA3rEq5cMKvaESVJUhU4IihJl2DPgeMUSmU2l1LSQ33Mbqjl46uXks8mvPPq+dQ4AihJks7CDpYkjertH+ThHV0USik/fO0QIcAtyxaSzybcvqqd5obaakeUJEmThCOCknQWw5XId17eT6GU8tiubvoHK1y3aDb5XML6dRmWzm2udkRJkjQJOSIoSad5uecom4pl7t9apru3n7nN9WzMJeSzCWuvaCMERwAlSdKlG9MOVghhFvAfgQ2jhwrAL8UY+87zmk8BvwcsBXYC/yzGWBx9rhn4MrAWuA7432OMf3CRWexgSfoxh44P8NCOTgrFlO3pEWprAh9cvoh8LuG2GxfTWOcIoCRJujgT1cH6D8ANwAogAvcDfwb8wjlC3QL8J2A98E3gV4CHQwjXxxh7R9/ju8D/DfzRGGeVNAMMDlf4xgv7KBRTnnx+L4PDkRuXtvI7H7+Re9dmWNTSWO2IkiRpGhmzDtZot+kg8IkY45Ojx24DHgLmxxj7z/KavwZqYoz/ePRxAPYAvxtj/Osz1n4DeMIOlqQLiTGyq7OXQinlwW2dHDg+wMI5Ddy7NkM+m7Cyo7XaESVJ0hQ3ER2sFUATUDztWAloBpYDO87ymjXAl049iDHGEMLW0eOXLITwaeDTl/NaSVNfz9F+HtjaSaGU8nz3URpqa/jIysXkswm3Ll9EfW1NtSNKkqRp7qIKrBDCl4CfOc+Sfws8Nvr9kdOOn/r+XB8Xt5yxHuDwedafV4zxC8AXYKSDdTnvIWlq6R8c5snneiiUUr754j6GK5E1V7Tx+/fdxN2rl9I2q6HaESVJ0gxysR2sXwJ+/TzPn2CkSwUwl5Ei6dT3AL3neN3R09ac0gbsvshckmagGCNb3zhMoZjy0PZOevuHaG9t4tO3Xks+m2HZ4pZqR5QkSTPURRVYMcZjwLHzrQkhvAD0A1ngqdHD64A+4MVzvGz76PpT7xEY2TFw88XkkjSzdB7uY8vWMoVSyiv7jtNUX8Mdq9rJ5xJ+4rqF1Na4tbokSaquMbsGK8bYF0L4W+BzIYRnRg9/Dvjy2Ta4GPVfgEdHN7t4GvhlRq7j2nJqQQihEQhADVAXQmgChmOMg2OVXdLkdWJgiMd2dVMolvnO7v3ECO+6Zj7/9NbruPPmdlqa6qsdUZIk6U1jfR+s2cCfc477YIUQ/jXwUzHGVae95lPAv+FH98H6hVP3wRp9/jXgqjN+1F/HGH/2AlncRVCaoiqVyA9eO0ihmPLwzi6ODwxzxfxmNqwbuRHwlQtmVTuiJEma4c61i+CYFliTiQWWNPXsOXCcQqnM5lJKeqiPOY113HVzO/lswjuvnk+NI4CSJGmSsMCSNCn19g/yyM4uCsUyP3jtICHALcsWks8m3L6qneaG2mpHlCRJegsLLEmTxnAl8p2X97OpmPLYrm5ODlW4dtFs8tmE9esydLQ1VzuiJEnSeU3EjYYl6bxe7jnKpmKZLVtT9vaeZG5zPT/5jpHrqtZe0cbIRqKSJElTlx0sSePq0PEBHtrRSaGYsj09Qm1N4IPLF5HPJdx242Ia6xwBlCRJU48jgpImzOBwhW++sI9CKeWJ5/YyOBy5ob2FjbmEe9dmWNTSWO2IkiRJb4sFlqRxt6vzCIVimQe2lTlwfIAFsxu4d22GfC7Dqo651Y4nSZI0ZrwGS9K42Hf0JA9sK7OpmPJ891Eaamu47cbF5LMJH1ixiPrammpHlCRJmjB2sCRdspNDwzz5XA+FYso3XtzHcCWyJplLPpdw9+oO5s1uqHZESZKkceWIoKS3JcbItjcOUyilPLS9iyN9gyxpbWT9uoSNuQzLFrdUO6IkSdKEcURQ0mXpOtLHlq1lCsWU3fuO01hXwx03tZPPJrxv2UJqa9xaXZIk6RQ7WJLeom9gmEd3dbG5VObbL+8nRnjn1fPIZxPuWr2U1qb6akeUJEmqKkcEJZ1XpRL5wWsHKRRTHt7ZxfGBYZJ5zWzIJuSzGa5aMLvaESVJkiYNRwQlndXrB05QKKVs3pryxsE+ZjfUctfNS8nnEt519XxqHAGUJEm6aHawpBnoaP8gD+/solAs84PXDhICvO+6heRzGW5f1c6sBj97kSRJOh9HBKUZbrgS+c7L+ymUUh7b1U3/YIVrF80mn01Yvy5DR1tztSNKkiRNGY4ISjPUyz3HKJRStpTKdPf209pUx8ZcQj6bsPaKNkJwBFCSJGms2MGSpqHDJwZ4aHsnm0pltr9xmNqawAeWLyKfTbjtxsU01ddWO6IkSdKU5oigNM0NDlf41ov7KJRSnni2h4HhCje0t7Axl3DP2g4WtzRVO6IkSdK0YYElTVPPdvZSKKU8sK3M/mMDLJjdwD1rO8hnE1Z1tDoCKEmSNA68BkuaRvYfO8n9W8sUSmWe6+qlvjZw2w1LyOcSPrhiEfW1NdWOKEmSNCPZwZKmiJNDwzz1XA+biinfeHEfw5XImmQu+VzC3as7mDe7odoRJUmSZgxHBKUpKMbI9vQIm4pv8ND2Lo70DbKktZH16xLy2QzXL2mpdkRJkqQZyRFBaQrpOtLHlq1lCsWU3fuO01hXwx03tZPPJrxv2UJqa7yuSpIkaTKygyVNEn0Dwzy2q5tCKeXbL+8nRnjn1fPIZxPuWr2U1qb6akeUJEnSKEcEpUkoxsgPXj1IoZTy8M5ujp0cIpnXzIbsyAjgVQtmVzuiJEmSzsIRQWkSef3ACQqllM1bU9442Mfshlruunkp+VzCu66eT40jgJIkSVOSHSxpghztH+SRnd1sKqX84NWDhAA/cd0C8tmEO25qZ1aDn3dIkiRNFY4ISlUwXIl8d/d+CsWUR3d10z9Y4dqFs8nnEu5blyHT1lztiJIkSboMjghKE2j3vmMUiilbtpbpOtJPa1MdG7IJG3MJ665oIwRHACVJkqYjO1jSGDl8YoCHdnRRKKZse+MwtTWBW69fSD6X8JEbl9BUX1vtiJIkSRojjghK42BwuMK3XtxHoZTyxLM9DAxXuKG9hXw24d51HSxuaap2REmSJI0DCyxpDD3b2UuhlPLAtjL7jw0wf3YD967tIJ9NWNXR6gigJEnSNOc1WNLbtO/oSR7YVqZQKvNcVy/1tYHbblhCPpfwwRWLqK+tqXZESZIkVZkdLOk8Tg4N89RzPWwqpnzjxX0MVyKrk7nkswn3rOlg3uyGakeUJElSFTgiKF2kGCPb0yMUiikPbu/kSN8gS1obuW9dho3ZhOuXtFQ7oiRJkqrMEUHpArqP9LN5a0qhmLJ733Ea62r42Kp2NuYSblm2kNoar6uSJEnS+dnB0ozWNzDM4892s6mY8u2X9xMjvOOqeeRzCR9fvZTWpvpqR5QkSdIk5IigNCrGyA9fO0ShmPLVnV0cOzlEpq2ZfDbDhmzC1QtnVzuiJEmSJjlHBDXjvXHwBIVSyuZSmdcPnmBWQy133rSUjbmEd18znxpHACVJkvQ22cHStHbs5BAP7+yiUEz5/qsHCQHee+0C8tmEO25qZ3ajnzFIkiTp0jkiqBljuBL53u4DFEopjz7TTd/gMNcsnE0+m2F9NiHT1lztiJIkSZriHBHUtLd73zEKxZQtW8t0HemnpamO9dkM+WxC9so2QnAEUJIkSePLDpamtCMnBnlwRyeFYsq2Nw5TWxO49fqF5HMJH7lxCU31tdWOKEmSpGnIEUFNG0PDFb754j4KpZQnnu1hYLjCDe0t5LMJ967rYHFLU7UjSpIkaZqzwNKU91xXL4Viyv3byuw/NsD82Q3cu7aDfDZhVUerI4CSJEmaMF6DpSlp/7GTPLBtZATw2a5e6msDH75hMflswgdXLKahrqbaESVJkqQ32cHSpHNyaJinnuuhUEr5xgv7GKpEVidzyWcT7l7TwfzZDdWOKEmSpBnOEUFNajFGdqRHKJRSHtzeyeETgyxuaXxzF8DlS1qqHVGSJEl604SMCIYQZgH/EdgweqgA/FKMse88r/kU8HvAUmAn8M9ijMXR594D/C7wDqAJeBn4/Rjj/WOZW9XTfaSfLVvLFEopL/cco7Guho+taiefzXDLsoXU1ToCKEmSpKljTDtYIYT/AqwC1gMRuB/YHmP8hXOsvwV4bHT9N4FfAf4lcH2MsTeEcBewEHgYOAjcA/x34NYY4w8vkMUO1iTVNzDM4892s6mY8p2X91OJkLtqHvlswifWLKW1qb7aESVJkqTzGvcRwRBCMyNF0CdijE+OHrsNeAiYH2PsP8tr/hqoiTH+49HHAdgD/G6M8a/P8XP+Hvi7GOOfXSCPBdYkEmPkH/YcolBM+eqOLo6eHCLT1syGbIYN2YRrFs6udkRJkiTpok3EiOAKRsb4iqcdKwHNwHJgx1leswb40qkHMcYYQtg6evwtQgjtjHTItp/j+U8Dn76M7Bonbxw8weZSmc1bU/YcOMGshlruvGkp+VyG91yzgJoat1aXJEnS9HFRBVYI4UvAz5xnyb9lZNQP4Mhpx09933qO17WcsR7g8NnWhxBmM3JN11dPdcjOFGP8AvCF0fW2r6rk2MkhHt7ZRaGY8v1XDxICvPfaBfzyh6/njpvamd3o3QEkSZI0PV3sv3R/Cfj18zx/gpEuFcBcRoqkU98D9J7jdUdPW3NKG7D79AMhhBbgq0AP8KmLCayJNVyJfG/3AQqllEef6aZvcJhrFs7m1z+2nPXZhExbc7UjSpIkSePuogqsGOMx4Nj51oQQXgD6gSzw1OjhdUAf8OI5XrZ9dP2p9wjAWmDzaccWAI8ArwA/HWMcupjMmhi79x2jUEzZsrVM15F+WprquG9dho25DNkr5zHyP6kkSZI0M4zHLoI38qNt2jcDz8QY/+k51t8CPArcCzwN/DLwGX60i2A78AQj13X9kxjj8CVkcZOLcXLkxCAP7ehkUzFl2xuHqQlw6/JF5LMJH125hKb62mpHlCRJksbVhNxoePQ6qT/nHPfBCiH8a+CnYoyrTnvNp4B/w4/ug/ULp90H6/dGnzvByLbvp/xhjPEPL5DFAmsMDQ1X+NZL+ygUy3ztub0MDFVYsaSFfC7DfWszLG5tqnZESZIkacJMSIE1mVhgjY3nunopFFPu39bJ/mMnmTernnvXZtiYS1jV0eoIoCRJkmakidimXdPEgWMneWBbJ4VSyq7OXupqAh++YTH5XMKHViymoa6m2hElSZKkSckOlgAYGKrw1PN72VQs840XehiqRG7OzCWfzXDP2gzzZzdUO6IkSZI0aTgiqLeIMbKzfIRCMeWB7Z0cPjHIopZG1q/LkM8mrGhvqXZESZIkaVJyRFBv2tvbz5atZQrFlJd6jtFQV8PHVi4hn0t4/7KF1NU6AihJkiRdDjtYM0T/4DCP7eqmUCrz7Zf2UYmQvbKNjbkr+Pjqpcxtrq92REmSJGnKcERwBoox8g97DlEopnx1RxdHTw6RaWtmQzbDhmzCNQtnVzuiJEmSNCU5IjiDvHHwBJtLZTZvTdlz4ASzGmq586al5HMZ3nPNAmpq3FpdkiRJGg92sKaJYyeHeGRnF5uKKd9/9SAhwHuvXUA+m3DHTe3MbrSWliRJksaKI4LTUKUS+d4rBygUUx55ppu+wWGuXjCLjbmE9dmETFtztSNKkiRJ05IjgtPIK/uOUSilbCmV6TzST0tTHfety7AxlyF75TxCcARQkiRJqgY7WFPEkb5BvrKjk03FlK2vH6YmwK3LF5HPJnx05RKa6murHVGSJEmaMRwRnIKGhis8/dJ+NpVSvvbsXgaGKixfMod8NuG+dRmWtDZVO6IkSZI0I1lgTSHPd/dSKKbcv62TfUdPMm9WPfeuzZDPJtyUaXUEUJIkSaoyr8Ga5A4cO8kD2zoplFJ2dfZSVxP48A2LyecSPrRiMQ11NdWOKEmSJOkC7GBV0cBQhaee76FQSvn68z0MVSI3Z+aSz2a4Z22G+bMbqh1RkiRJ0lk4IjhJxBjZWT5CoZjy4PZODp0YZFFLIxvWZdiQTVjR3lLtiJIkSZIuwBHBKtvb28+WrWUKxZSXeo7RUFfDx1YuYWMu4ZZlC6mrdQRQkiRJmursYE2A/+vJl/g/n3iRSoTcVfPIZxM+vnopc5vrqx1NkiRJ0mVwRLCKnn5pH99/5SAbshmuXTSn2nEkSZIkvU0WWJIkSZI0Rs5VYHnhjyRJkiSNETe5GGf3by3zJ4+9QOfhPjramvnM7Su4b12m2rEkSZIkjQNHBMfR/VvL/NbmnfQNDr95rLm+lj/acLNFliRJkjSFOSJYBX/y2As/VlwB9A0O8yePvVClRJIkSZLGkwXWOOo83HdJxyVJkiRNbRZY46ijrfmSjkuSJEma2iywxtFnbl9Bc33tjx1rrq/lM7evqFIiSZIkSePJXQTH0amNLNxFUJIkSZoZ3EVQkiRJki6RuwhKkiRJ0jizwJIkSZKkMWKBJUmSJEljxAJLkiRJksaIBZYkSZIkjRELLEmSJEkaIxZYkiRJkjRGLLAkSZIkaYxYYEmSJEnSGKmrdoDxFMJbbqwsSZIkSeMmxBirnUFTTAghxhitXjXuPNc0kTzfNFE81zSRPN8mniOCkiRJkjRGLLAkSZIkaYxYYOlyfLbaATRjeK5pInm+aaJ4rmkieb5NMK/BkiRJkqQxYgdLkiRJksaIBZYkSZIkjRELLEmSJEkaIxZYuqAQwqwQwl+FEA6Pfn0xhNB8nvWfCiF8N4RwKISwP4TwSAjh5onMrKnpMs61NaPnV3cIIYYQbpnIvJo6Qgi1IYQ/CSHsCyEcDSEUQggLz7P+jhDCrhBCXwjhmRDCxyYyr6a2SznfQgiZEMIDIYQ9o3+P/fRE59XUdYnn2l0hhKdG/212KITwdAjh/ROdeSawwNLF+A/ADcAKYDlwI/Bn51nfAvwekAAZoAQ8HkKYNc45NfVd6rk2AGwGPjH+0TTF/SZwL/BuRv5uAvibsy0MIVzLyHn1R8Dc0T+3hBCuHv+YmiYu+nwDKsDjwP8MpOMfTdPMpZxr84A/B5YBi4D/BjwSQrhivEPONO4iqPMa7R4cBD4RY3xy9NhtwEPA/Bhj/0W8RxPQB+RijKXxzKup6+2eayGECLw/xvjtcQ+rKSeEsAf4XIzxi6OPrwNeBq6OMe45Y+1ngQ/HGN9/2rGngSdijG53rAu6lPPtjNe9BvxOjPFvJySoprzLPddOe3038M9ijJvHN+nMYgdLF7ICaAKKpx0rAc2MdBguxm3ACeClsY2maWYszjXpLUIIbcCVnHZuxRh3A73AmrO8ZA0/fh7CyLl4trXSj7mM8026LG/3XBu9fGMhsHOcIs5YFlgzWAjhS6Pz3uf6+gNGxv0Ajpz20lPft17Ez1gO/D/Av4wxHh3b30BTxUSca9J5nO3cAjjM2c+tlktYK53pUs836XJd9rkWQlgMFID/I8boB+BjrK7aAVRVvwT8+nmeP8GPOgdzGfkP9tT3MPIJyTmFEFYCX2PkP97/fPkxNQ2M67kmXcCpD3fmnnG8jbOfW0cvYa10pks936TLdVnnWgihg5F/nz0O/Na4JJvh7GDNYDHGYzHG/ef5OgG8APQD2dNeuo6Ra6pePNd7hxCywDeAz8cY/904/hqaAsbzXJMuJMZ4GHid086t0Y0sWoEdZ3nJdn78PISRc3H7OEXUNHIZ55t0WS7nXBvdrOdp4JEY4y9FN2MYFxZYOq8YYx/wt8DnQgiLR1vKnwO+fK5NB0II7wOeBH47xvjnE5dWU9llnmshhNA0upEKQMPo49oJiq2p4wvAb4QQrgkhtAJ/DDwWY3ztLGu/DLwjhPDJEEJ9COGTQA7464mLqynuUs43Tvt7LAD1o4+dMtLFuOhzLYRwA/Bt4L/HGM83VaK3yQJLF+NXGekgnPp6AfgXp54MIfzrEMKu09b/ASPt6n8fQjh22pf3WtCF/CqXdq5dxUiHq2/08ZOj3//jiQirKeXzjOxI+UOgDNQCPw0QQvipEMKxUwtHLxLfAPwOI2M2vwOsP9c/jqWzuOjzbdSpv8euBP5q9PvfmbC0msou5Vz7DUZun/OrZ/z77KcmOvR05zbtkiRJkjRG7GBJkiRJ0hixwJIkSZKkMWKBJUmSJEljxAJLkiRJksaIBZYkSZIkjRELLEmSJEkaIxZYkiRJkjRGLLAkSZIkaYz8/1Po2p60KsZ1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# use neutral structure and datapoints of disp_002,004,006\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from nonrad.ccd import get_Q_from_struct\n",
    "from nonrad.elphon import get_Wif_from_WSWQ\n",
    "from pathlib import Path\n",
    "from pymatgen.core.structure import Structure\n",
    "from glob import glob\n",
    "\n",
    "cc_dir = \"/pscratch/sd/z/zkyuan/todos/kesterites/CZTS/defect/ZnCu2SnS4/static/nonrad\" + \\\n",
    "         \"/\" + \\\n",
    "         \"/as_1_Sn_on_Zn/q0_capture_h_-1_1_21/excited\"\n",
    "\n",
    "SnZn_4 = Structure.from_file(f\"{cc_dir}/4/CONTCAR\")\n",
    "SnZn_5 = Structure.from_file(f\"{cc_dir}/5/CONTCAR\")\n",
    "SnZn_6 = Structure.from_file(f\"{cc_dir}/6/CONTCAR\")\n",
    "\n",
    "\n",
    "WSWQs = []\n",
    "for d in glob(f\"{cc_dir}/*\"):\n",
    "        if d[-2:] in [\"/4\", \"/5\", \"/6\"]:\n",
    "            pd = Path(d)\n",
    "            Q = get_Q_from_struct(SnZn_5, SnZn_6, f\"{pd}/CONTCAR\")\n",
    "            path_wswq = f\"{pd}/eph/WSWQ\"\n",
    "            WSWQs.append((Q, path_wswq))\n",
    "print(WSWQs)      \n",
    "fig = plt.figure(figsize=(12, 5))\n",
    "ground_vr_path = f\"{cc_dir}/5/vasprun.xml\" \n",
    "\n",
    "Wifs = get_Wif_from_WSWQ(WSWQs, ground_vr_path, 2305, 2306, spin=1, fig=fig)\n",
    "\n",
    "plt.tight_layout()\n",
    "# plt.ylim(-1e-6,1e-6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3c47332",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}