postprocess-lpg/DOS/d_band_moments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The system used as an example is PdML/Pt(111).\n",
    "Here we are looking at the projected density of states of the dband of Pd.\n",
    "Vaspkit was used to obtained the postprocessed data as PDOS_USER.dat\n",
    "This script can be used on other other similar file formats.\n",
    "\n",
    "OUTLINE\n",
    "1. First visualize the total pdos.\n",
    "2. Choose the range of integration.\n",
    "3. Calcualte the momemnts (center of the band and the width as the sqrt of the second moment)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAGFCAYAAABDgaHaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XecXFX9//HXZ2Z3s5uQTk+AJIjU0IIgQUKLIs0KYqFZflaUSFGsKCAdQdCvIiqINFERBFEktIAgSDCAlCCEQAgE0nuy7fz+uHNn79y5907JzE7Z9/PxiLNz61lM5n7mnM/5HHPOISIiIlKqVK0bICIiIo1JQYSIiIiURUGEiIiIlEVBhIiIiJRFQYSIiIiURUGEiIiIlKWl1g3oL2bmACZNmlTw2JkzZ9Ksx9ZLO+rh2HppRz0cWy/tqIdj66Ud9XBsvbSj0Y6tl3b4xzrnrODBZbKBUifCDyKK+X3NvP/ezXhsvbSjHo6tl3bUw7H10o56OLZe2lEPx9ZLOxrt2HppR+DYqgURGs4QERGRsiiIEBERkbIoiBAREZGyDJjESuec7bXXXkUNZpWSJ9Jox9ZLO+rh2HppRz0cW81rH3DAATVvQ738t2i0Y+ulHY12bL20wzmHmc0s+oQyqCdCREREyqIgQkRERMpSd8MZZnY1sJ1z7sAijh0PXAr4x94JnOacW1i1BopISR544IFaN0FEqqSueiLM7LPA54o8djRwP/Bu4EK8YOIDwD1m1la1RoqIiBRp3LhxjBs3riLXOvDAA7O1H+pFXQQRZpY2s+8DV5dw2qnAWOAQ59yFzrkfAUcDuwEnVqGZIlKGSy65hEsuuaTWzRApmpnl/Emn02y88cYcfPDB3HDDDbVuXl2p+XCGmbUDjwG7AtcBhxR56seBB5xzz/sbnHPTzWx2Zl8pAYmIVMmdd94JwOmnn17jloiU5qyzzgKgq6uL2bNnc9ttt3H//fczc+ZMfvzjH9e4dfWh5kEE0A4MA451zt1iZnMLnWBmI4EJwB8jdj8JHFHRFoqIyIDzgx/8IOf9vffey3vf+14uv/xyvva1r1VsmKKR1cNwxgq8RMpbSjhnTOZ1fsS+N4FhZjZ8g1smIiKSccghh7DDDjvgnOPf//53drtzjp/+9KfsvPPOtLe3M2bMGE4++WSWL19e1n1uvvlmJk2aREdHB5tuuinHH388b7zxRqV+jYqqeU+Ec64X6C3xtKGZ1zUR+9ZmXocAOf8Pzpw5My8pZaAsQCaN45Z/z2P52i7+35QJtW6KiIT4z4zgs2TatGlcccUVbLHFFnz+85+ntbWV22+/nccee4zOzk7a2orP9b/ssss49dRTGTFiBCeccAIjRozg7rvvZvLkyQwfnvzduBZJlzUPIsrk96AkRQClBiYideEbf3oaoGmCiI6Ojlo3QaQipk+fzuzZszEz3vWudwHwyCOPcMUVV7Dtttvy+OOPM2rUKAB+9KMfcdBBB/Hmm2+yzTbbFHX9uXPncuaZZzJy5EiefPLJ7HDJ+eefzzHHHMOtt95ald9rQzRqELEy8xr16dQROiZr0qRJPPHEE1VrlIjk+9vf/lbrJoiUxc+JCCZWOuf4+te/ng0MrrnmGgC+853vZAMIgPb2ds4//3wOOuigou93ww030NnZyZlnnpmTb5FKpbj44ou57bbb6O2N/34c1bNe7d6JRg0iXsu8bhGxb0tgmXNudT+2R0REmswPf/hDwHsQjxgxgv3335/PfvazHHfccdljnnzySSB6jZj999+flpbiH7NJ15owYQJbbbUVr776akm/Q7U1ZBDhnFtmZq8Ae0bs3gNQd4NInTjnnHMA+N73vlfjloiUppicOT95crPNNsvbl06nGT16dNH3S7oWwOabb153QUQ9zM4o15+AqWa2g7/BzKYC2wM316xVIpLj3nvv5d577611M0Sqwk92fOutt/L29fT0sHjx4opcC2DBggVltLC6GqInwswmAJOBR5xzczKbLwJOAO41s0vx6k18A5gJXF+ThoqISEX88I5nee6NFRW95k5bDuOso3au6DX33HNPnnzySR588EEmTMhNhn7ooYfo7u4u6Vq33norDz74IAcffHDOvjlz5jBv3ryKtLmSGqUnYgrwu8wrAJlFtqYATwFnA9OA24DDnHPra9FIEREZWE466STAm42xZMmS7PZ169bxrW99q6RrfepTn6K1tZUrr7ySuXPnZrf39vZyxhlnJCZV1krd9UQ458ZFbLsWuDZi+2zg8Ko3SkRE+lWlewyqZb/99uOrX/0qV155JbvssgtHH310tk7EyJEj2WKLqPz/aOPGjeOCCy7gtNNOY4899uDYY49l+PDh3H333Sxbtoxdd92Vp59+uoq/TekapSdCRBrU6NGjS0ouE2k0P/nJT7jyyisZPnw4V111FTfddBOHHnoo06dPL6nQFMCpp57KjTfeyPjx47n22mv5zW9+wy677MIjjzzCyJEjq/QblM8GUsXGvfbay6lOhNS7cWf+FYC5F2gJGBHZMGY20zm3V7Wur54IERERKYuCCBGpqm9961slJ5iJSGOou8RKEWkujz76aK2bICJVop4IERERKYuCCBERESmLgggREREpi3IiRKSqxo4dW+smiEiVKIgQkaq6/notZSPSrDScISIiImVRECEiVTVt2jSmTZtW62aISBVoOENEqmrWrFm1boKIVIl6IkRERKQsCiJERESkLAoiREREpCzKiRCRqvnv/OUsSo9i3wmja90UEakCBREiUjVfvH4mq971Wb7/jYNq3RQRqQINZ4iIiEhZFESISFUt/vuVfPPrJ9e6GSJSBRrOEJGq6loyn1deXlPrZohIFagnQkRERMqiIEJEqs/VugEiUg0KIkSkasxq3QIRqSblRIhIVbVtOoGdJo6tdTNEpAoURIhIVY2a+nnOOkN1IkSakYYzREREpCwKIkSkqhbdcQmnfPEztW6GiFSBhjNEpKq6Vy5iwRs9tW6GiFSBeiJERESkLAoiREREpCwKIkRERKQsyokQkaoaNGYH9njXNrVuhohUgYIIEakawxh5wEl884wDa90UEakCDWeIiIhIWRREiEhVLfzzeXzxpE/WuhkiUgUazhCRqupZu4KlS/R9RaQZ6V+2iFSN0xrgIk1NQYSIiIiURUGEiFSNYbVugohUkXIiRKSq2rfZjX33G1/rZohIFSiIEJGqGrHfJzj5tANq3QwRqQINZ4hI1TnlV4o0JQURIlJVb91yFp/9xIdr3QwRqQIFESJSNWbgutezbu26WjdFRKpAQYSIiIiURUGEiFSdUiJEmlNdBBFmNt7MbjWzJZk/15nZJkWcN8nM7jGz1Wa2wszuMLPt+6PNIiIiA13Np3ia2WjgfqANuBCvTWcAu5rZ3s65zpjztgceANYAZ2c2nwY8bGa7OefeqHbbRWrhv/OX8685i/nc/hNq3ZSidGy7NwdOUZ0IkWZU8yACOBUYC0x0zj0PYGaPAfcAJwJXx5w3DdgImOKc+0/mvPuAx4Gv4wUiIg3LOYdZfsXHI698GKBhgojh+3yET39x/1o3Q0SqoB6GMz4OPOAHEADOuenA7My+OBOARX4AkTnv38BiYGKV2irSb5qptkIz/S4i0qemQYSZjcQLBmZG7H4S2Cvh9P8Bo4K5E2Y2ChgBvFnJdorUQrM8dxfceCYnHn14rZshIlVQ656IMZnX+RH73gSGmdnwmHMvAl4HbjKzXc1sInAT0AlcUfGWivQz10xf35voVxGRPrXOiRiaeV0TsW9t5nUIsDy80zn3mpmdB/wUeCqzuQc4OjjEETRz5sy8Meam+qCWpqK/mSJSiqgcqmqrdU+Ef/+kz8veqI1mdg7wC+CfwKeAE/CSKm8xs6Mq2UiRWmiG+FYLgYs0t1oHESszrx0R+zpCx2SZ2Qi82RdPAIc45250zv0OOAB4DrjazAaFz5s0aRLOuZw/IvXKqS9CREoQfr71xzOu1kHEa5nXLSL2bQksc86tjti3HTAIuMk51+NvdM51ATcAmwE7VLitIv2qWWLcITvsz3uP/FCtmyEiVVDTnAjn3DIzewXYM2L3Hng9DVHWZ17TEfv8bbUOkEQEGLrnEXzshP1q3QwRqYJ6eND+CZhqZtmeAzObCmwP3BxzzrPAG8BJZtYeOK8dLzdiEfDfqrVYpB80S09Eb9c61q6Nyp0WkUZXD0HERcAS4F4zO9XMvg38Ea92xPUAZjbBzI4zswkAmSGMk/GGLB43s1PM7FS8nosdgGmZoQ2RhtUsORFv/+EHfPXEj9W6GSJSBTUPIpxzC4EpeNM0z8YrZ30bcJhzzh+2mAL8LvPqn/dn4L14Ach5wDnAUuBw59wN/fYLiIiIDFC1rhMBgHNuNhBb0s45dy1wbcT2+4D7qtYwkRpqluEMaK7fRUT6lBVEmFkLsDEwEq8o1HwNH4hUlp67IlLvig4izGwf4BPAgcAu5NaR6TWz2cA/8KZd/ruSjRQZiArN8Y5b5VNEpL8UDCLM7HDg28C+eIHDEuARvBkQq/EWvNoY2BEvn+EUM3sI+JFz7p4qtVuk6RXqiXAOGiGG2GjiVI48aNtaN0NEqiA2sdLMNjOzPwF34n2efQXY1jm3sXNuinPuI865451zRznn9sUb2tgdOBMvqPi7mf3FzKIKSYlIAY2YR/DDO55l3Jl/zb43MzaaOJUjPvrJGrZKRKolaXbG83hVIfd0zu3vnPuFc+6VuIOd52nn3MXOuV3wZlIMwavpICKlKhBE1GOMcc0/5+Zt61mznKWLF/V/Y0Sk6pKCiGOcc0c652aVc2Hn3D+dc4cAHy+vaSIDW6E6EfW89su/5ixm+Rov13rhbefzrZM/XeMWiUg1xOZEOOfurcQNnHP/qMR1RAaaOo4RCvr4L//F3uNG1XWgIyIbLikn4iUz+56ZbdOfDRIRT8HEyn5pRfn+93ZgAd56b6yIlCVpOGNr4AfAy2Z2r5mdYGaD+6dZIlJ4imc/NaRMremUpqCKNLmkIGJzvBkZj+LVhrgGWGBmvzGzKQnniUgFFO6JqO8oojVd86r6IlJlsf/KnXNLMjMy9gfGAd8BXgVOAu43szlmdpaZje+XlooMMPXe01BIW4v38TJ0j8P50CdPqm1jRKQqivqq4Jyb55y7wDk3EdgNuCRz7lnAS2Z2v5mdZGZDqthWkQGl8OyMfmpImVrT3lDGkB2ncPDhH65xa0SkGkrub3TOPeOc+6ZzbhzeMMfVwE7Ab/CGO66tZANFpDH5wxndKxay4M3Xa9waEamGDRq0dM7NcM59ERiDV6nSgOMr0TCRAc/B4lXrWbRqfa1bUpbWdAoDFt15KT864yu1bo6IVMEGBRFmdqCZ/RyYB5wPpIE/VaJhIgOdAyadO529zp0evb+I4YxXF6/mO39+hp7e/h/7aAskVtb5yIuIlKnkIMLM9jWzn5jZfOBe4AvAXLyZHFs45z5W2SaKDEyFgoRiZmecfON/uOGx1/jv/OUValXxWlsC0zsVRYg0paKWAjezSXjlqz8GjMUbtngNr/fht865/1WthSIDVCWncNaiXMOC5euYs2h1/99YRPpNbBBhZhOBYzN/JuAFDquB3+EFDvf3SwtFBqiCPRF1/u3+5YUKIESaXVJPxFN4nZAOuB/4LfAn59ya/miYyEDX6GWvfcP2/jDHHLJdrZshIlWQFES8iBc4/M45p/lZIv2scNnrxggjBr9jH/Y9aK9aN0NEqiBpFc8dCp1sZkOcc+qzFKmCBokRYrWmja4eR9fi13ltzgjYYbNaN0lEKqyk2Rnm+aKZPWZm64Dlme0nZ9bU0KeESD8pJsaISs68fdZ8PvfbJyrfoJBBLWkAFt/9Uy77wRlVv5+I9L+iZmcAmFkLcDvwfqALWAmMyuwej7emxv5mNtk5t7DC7RQZcCqZWGn0Tc845eZZZbaoNL2N3pUiIgWV0hNxOnAYcBle8PCzwL5v4q2jsS3wrYq1TmSACeY5FDPF8+0V65j56pJqNqls3cECV4onRJpSKUHECcA/nXOnZ2ZoZD8WnHPdzrlzgPuAIyvcRpGm8uwby3l7xbrIfcEv7wW/yDt472Uz+OjPH61c4yooWCVTMYRIcyoliJgAPFTgmCfwilGJSIwjrniY91wUXWbFxfwcfaxj+dquirWr0mpRaltE+lcpQcQyYJsCx2xLJtlSROJ1dvdGbq/XaZvjzvwr3/nzM2WdO3zyx/nE50+pcItEpB6UEkRMBz5iZrtH7TSzdwMfxFtPQ0TKkNMTUbBORBHXq2BMcsNjr5V1Xse43dlj3ymVa4iI1I1SgoizgPXAP83sCmAfADM70cx+ilfVch1wbsVbKTIAVbJiZS3WzvB1vjWHl57/b+0aICJVU/QUT+fcy2Z2MF4Vy5MDu36Dt67GK8AJzrkXKttEkYEjN7GyOSpWLrn3l/zyv0M5+Zj31ropIlJhRQcRAM65J4GJZrYPMAkYAawCngZmOOeiB3pFpCjBaZ0NEiMUZV1XT62bICJVkLSK5y3A34F7nHPzgvucc48Bj1W5bSIDzo2BvINCkxsaKcZ4bYnW7RNpRkk9EX5BqTYzmw3cjRdUzHDOre2PxokMNJf+48Xsz4UqPjZTT4WINKbYxErn3FS8QOJIvADiUOBvwBIzu8fMTjezXfunmSLN7aH/LWT52i5a0n0ZkNUOEholp0JE6ldiTkSmx+FvmT+Y2VZ4a2e8D/g2cKGZvQX8Ay/QmK51M0RKs2R1J8f/+nH2e8doWlJ9QUTBnogiBjSSLuFc/8zaGDHlxOrfRERqotTEynnA1cDVZpYC3o3XQ/E+4DgAM5vlnNur0g0VaVbru72kw5feXkW6hCBiQ5Miep0jxYZFEXMWrmL+srXsv90msce0j91xg+4hIvWrpCAiKDMT45HMn7PMbAReMPG+CrVNZEAxLCeI6K5y2ehKXP3gSx8EYO4FR8Qes+715zM/xR8jIo2p7CDCZ2Y7AZOBuc65W4BbNrhVIgNIsMMhJ4joKTScUbn75m53rO7sYaNBLdn3hfT0upy2By2b8dvMT6eX00wRqWOlVKzEzL5rZm+aWUfm/THAU8BVwN1mdq+ZtVehnSJNz4ycnIjunuSyKxuaFxmXU/Hrh19hl7Pu5s3l3iSsYjpElq3p3LDGiEhDKjqIMLMvAmcDvcBIMzPgx0A3Xknsa4GDgG9UvpkiA0Pw23xXgad3MYmV2WMjDo0LQv723wUAvL7UDyLi7+MnZi5VECEyIJXSE/EZ4CVge+fcG8B7gDHAjc65c51znwXuAz5R+WaKNK/gI7ol3fdPslBPRM41Yh70Lvuavz8uNggPSiQt6T2ioxWAxasURIgMRKUEETsCtzvnVmXeH473GfWXwDFPUHi5cBGJYOQOZ3QVyolw0T8XOja7rUBPhn9O0rXbWryPkHUxS5uLSHMrJbGyi9yg4zCgB2/1Tt8oYFkF2iUyYAR7EXJnZxTIiYj5udCxffeNPtYfovDblTSc0etyj40y6pDPF2idiDSqUnoingOOMLPBZjYF2BV4yDm3AsDMxgPHAFrzV6QMZhZKrCwl5yFmOCMxAIjeZ6EBjeBxh142gxkv9tWT86+f1NK2zSbQttmEhCNEpFGVEkT8DHgnsAC4F+9z4woAM5sGzAKGAhdWuI0iA0ZOYmXB2RmBFT8LXDcqmCg2RAmmRMx+ayXfva3ve0I2XyLhYmvnzmLt3FlF3k1EGknRwxnOuZvMrBc4A2/49irn3O2Z3YOBecC3nHP3Vr6ZIs0r+HxvDSZWFpqdUUpORIHzk84JByDBUtl9MUT8xZY/cnPmp+8k31BEGk6pZa9/D/w+YtfFzrnzKtMkkYErVUJPRFCxSZKhk6KFpmeEYxnL2efiry8iTa+kYlM+M9vRzD5hZl/JbNrCzDYqtxFmNt7MbjWzJZk/15lZfDH+vvM2MbOrzewtM1thZg+Y2bvLbYdIrQVzIjoDMx7WdfUwb8ma2PMKP8TzDyh2qfHwcWb5K40qiBAZmEqtWLmTmf0LL3nyeuAnmV2fBl43s2NLbYCZjcab4fFuvHyKS4EPAPeYWVvCeUOBGcDHgJ8D38OrW3Gfme1SajtEas2MbKlpgM5AT8QXfjeT/S+6P+f4Uh7c0VM8Y9qR3V/M7IzCiZUi0ryKHs7IzL6YgZc8eSOwOXBwZvcreAHJDWY23zn3cAltOBUYC0x0zj2fuddjwD3AiXirhkY5E9geONA5NyNz3u+BOXhVM08ooQ0idWHrUYOzP3d19z2aH8zMiMhNpgz8XFZORMzsDAsfF9of+NlPrPSv1Zo2unocbelUThAkIs2plJ6Ic/ASKPd1zh0PZAMF59x1eD0Ja4BvldiGjwMP+AFE5nrTgdmZfXkyJbdPBP7qBxCZ8xbgrfLzUIltEKmZuAAgKicidtGsMnIiih0ByeuJCEQR2eGMzPvWdIrPvWc8W43qyB4z+tCTGX3oyYXuJiINqJTEyvcCtzjnnoza6Zx7zsz+gFfJsihmNhKYAPwxYveTxK8dPA5v6OKizHUMGOKcW+Wc+79i7y9SD/wAwCw3GOiKKDaVU2CqpIqVpZS9DteJCO8P7stNrHQu83sEZ5yMHpvcOBFpWKX0RAwD3i5wzDJgeAnXHJN5nR+x701gmJlFXW+7zOvbZnZx5r4rzewlMzuqhPuL1Jz/wDUs5+EbHM7w9eYMZxD5c+Q9Iu9b3FLjvQlTTfva05c/kTLLaeealx5jzUuPFWihiDSiUnoi5gD7x+3M9AYcmDmuWEMzr1Fp52szr0OA5aF9IzKv5+CV4z4FrwT3GcBtZnZoZkgkx8yZM3Myy6HwB6lItcUFA509PfnHxvx1jUt+TJo9EZtYmS17HX1u8N9Qb+gYB2DQEzhpxeN/zvx0dswdRaQSws+3/lBKT8T1wD5mdp6ZpYM7zGwQ3rLguwM3R51c4P5JT/Ko7KxBmdcRwH7OuWudc78DpuD1SpxfQhtEaioukC3YE5Hzc4F7lLCKZ9I9IX+Vz9CNMIwCy36ISJMoJYi4BG8q5pnAQuDzAGb2APA6Xm/A48DFJVxzZea1I2JfR+iYoNWZ11udc0v9jc65ZXirik6KqlsxadIknHM5f0RqLS7PoWDZ69g3hQ72NxX39z+/TkT85R0Os8I1KESk8sLPt/54xhUdRDjnuoBDgW/i5UZsgfelZAreQ/1s4CDn3PoS7v9a5nWLiH1bAsucc6sj9vk5FFE5Gm9n2lV28SuR/uT/Q89PrCy+QFTB2RkR2+JSHbLDGdk8h9D+iL6InMRKAmtqiEhTK7XsdTdeT8PFZjYEL4lylb+SZ6mcc8vM7BVgz4jdewBPxJz6X2A9sHPEvvHAOrzeEpG615dYmauzOzknorTZGVHbYupEUHrekB9wOMj0RBQ8RUSaQNE9EWZ2n5llCzg551Y7594IBhBm9lUze6HENvwJmGpmOwSuMxWvkFRkfkWmd+IvwJFmtnPgvPF41S5vd87lfwKL1KG4YYmuiKXAXcy7wqMZpedE9JW9zt0eOZyR7YlwXk5E4OIbH3kaGx95WoEWikgjiu2JMLMOoNV/izfz4lEzGxZzShvwPmCbEttwEV51yXvN7FKgHa/i5Ey8ZE7MbAIwGXjEOefP/vhGpk33m9lPgE68vIy1wLdLbINIzeT0KAS2R+VE5CZWBq+RPMxRytBo33BG/j3jBCd6hnMiWoYVXAZHRBpUUk/EZ4ClmT9LMtvODGwL/3kLr9DUf0ppgHNuIV5exVN4eRXTgNuAwwL5FVMAf/aFf95cvCqZD+JN7fwuMAuYHAg0ROpe3IyL4AJcffujr1FOnYhikx+TFuDqa1dfsBLOiVj9/AxWPz8j7xwRaXxJORE/x3tob5p5PwUvEXJuxLEOr17DfOBHpTbCOTebhEqXzrlrgWsjts8Bjin1fiK+u555ky1HdLD7ViMKH1wBSZUjww/n6LLXMcWmKlixMnxOyQnells0a+V/7sr8dGGJFxKRehcbRDjneoHsqpxm1gtc45xTxRhpGl++waviPveCuArrlRVd9CkzOyO0PzInIi6xsozZGfHFpixnfzF1IpwLzDJBszNEBoqiZ2c450paNlxEilNuTkTuRQrdJOq+xQ5n5L6PrhPhAj0q+e1M1aCSnohUX0lTPAHM7B14Qxxp+r6UGF4S5mjgcOfciRVroUgTKfisDxwQtZR2zhDGBs7OKLazoKjEStfXhvDsDIgOPESk8RUdRJjZxsBfgb2KOFxBhEiEqG//2Qdu6EFb2uyMQveN3Bp5rIV2h9scN8UzWDQrHKCogKVIcyqlJ+Ic4F3AM8DDwKeA5/FmY+yEl3j5FvD+CrdRpKnF5TZErZ0R1+VQ1gJcEdt+9NfnePDFhTltKapiJcGeiNz2bPKhbzGoVaOhIs2olCDiMGA2sIdzrtfMNgM6nHNfBsgUoroWb9rlU5VuqEgzSEpwDCdWdkeWvQ6cF5NLUep9g65+6JW+e2U6QsJJknFDE8GciGDb0oOHk25RECHSjEr5l70F8I/MrA3weiDe7e90zl2HV7NBQxkiMUopPx3Vu5CbBxFdXyL6vgnDKDH8/R//5b9ytkfPznB9s0wCUcavTtiLnVbOZPlT9yTeS0QaUylBxFq8NSl8LwMjzGxMYNvjwIRKNEykGUWWnw6995/B3ZF1Igr/HH3fiG1lnANkGxiuWRF1vS1HdPDCQ3ey4unpyTcTkYZUShDxPIGeB7yhDcNbKMs3DBhcgXaJ1MSTry3l+F8/VnAZ7krqS0g0nHOkMw/pqFoL5S6xXWxORFS74uQ0L/BzcLijrcXCu0WkiZQSRNwM7G9mvzWzccDTwBvA2Wa2g5kdCHwCeLHSjRTpL6ff8hQP/W8Rry9dW5XrF3qYOyCV8h68UTkRsWUiCj6lk4dGosQuFZ7dnzu00jfJpC+KaE2nIoc/RKQ5lBJE/Ay4FTge2D+zSuY5wO7As8C9eD0Rqm0rUoJgYiWQ2BNRasXKpDAhHHj0hu4XF5j4PQ3h6aZ9ORF9x7amlVAp0sxKqVjZDRxtZvsAr2e2XWVmS/B6INYB1zvn7kq4jEhDKLaaY+nXTd7mHKT8nIioICKBdVL2AAAgAElEQVQuWIjtoXBF3Regqzd3CCdu6MSPEcI9KFHlLrJBhApFiDSlUopNbQ0sc849FtzunPsD8IfMMZub2cHOufsq20yR5uVCwwJJJaJ7Y3ofypvimbs1vGposbM3IFSxMpgTkU7xuXN/yS9mvFyghSLSiErpa3wFOKXAMafgLeMt0piqPICfVH7af/j6ORGR58dWrCw0xbPwtqgFvwDGbzwk570/hTMc0PQtwBXIiWgxBrV3YC3tie0TkcYU2xNhZu8FdgxuAvY1s6/FnNKGt+pnT+WaJ9JcklfxtJzhjChxyY6FeyIKTy2N64lY35X7TzoysTKmJ6I1neKfd9zIytkLgf5ZKVVE+k/ScMZS4Md4nxmG9xlxKIXLWv+0Mk0TqZ1qjeBHXjeUW5BOiCJeX7om8lrlrJ0RHq4IT2v1d0ctBAbgApvj6kS0pIynZvyd1UvX4JzLKUQlIo0vNohwzj1hZh8ANsELIn6DN1Rxe9ThQBcw3zk3oxoNFekPtXjE5QcD8a046Zp/B44tXLHShV5z9oU2rs/riYjeHjU7I9gVEQwUgj87p9U8RZpNYmJlcKaFmR0A/Nk595eqt0qkxqo1mSDqYR/eVOyDNqbWU9H3DZ8VN5wR3u7nPIR7MvqGZWLaUKCNItJ4Spni+emo7WaWBsYDC5xzqyrVMJHaqtIUz4ht/sPYCx4qV5Gy1HPyqnQ6L/jIG87I9kTkHJqzAFfE4ZlARl0RIs2kpEowZjbFzH6fCRwws93wZm3MBt42s7Oq0EaRflPtMfvoxMrc/cW2oJhiU0n3DSdphoOIXufo6nH5PSXZa+YmVvrBUNwUVfVEiDSfooMIMzsYryrl0cBWmc1XA2OB+4G5wPfN7LgKt1GkqYUfxsXHMdHTPaNM+/0szr/r+dj7Qn6FzF4X0TuB1751XT3sfd69OdeKmp0BcMpl17P5Jy9QvSmRJlRKT8Q3gJXA3s65uWa2I7AXcLdzbipe+esXgK9Uvpki/atqD7wCPRGQW2eh6MsW0d6rZsxJvG+4Z8Lh6I6pHfHGsrWhY6MrVkJf706h3hIRaTxF50QA7wJuds7NzLw/Eu+z4xYA51ynmf0d+H+VbaJI84h8kJYwLJFzWpnnRZ3vvY/oiejN74loTafyhiyCa2dgxvRTp7BoVScA03//K5bPXohzhWaHi0ijKSWIaAeWB94flnm9J7AtBXRvaKNEaqXaaX9J9Rq8pcDLnJ0RF0MkxBbhoCF8qHMuchGwQS0RQUTgAga8Y9OhvGNT7/2z/7qftUvWICLNp5Qg4iVgHwAzGwvsBzzrnHs9s60NrySdiuRLw+vPjvfwMEAlApkXFqzg/Zc/lHzf0Pu8KZuxORGWF+gk5UQErycizaWUnIhbgQPN7H7gQbwA5BoAMzsCeBTYFi/ZUkQiFKrWUMpzNm5ZcIC7nllQ0vmQnxPR66JzIrzKk/HXi8vpUE6ESPMppSfiXGBzvJwHA34PXJHZNxnYDa9MtoIIaVjVrqgYXWwqPDujuEaEV/8suS2Eex7yeyK6I3Iiel3+NM5gToR6IkQGjlKKTfUAXzKzbwAp51wwP+Jq4Arn3FuVbqBILVStYmXEtrjlvUu5ViWKTeX3TLjIlT2jeiK8VTy9n8MxRNugdqylR/0QIk0oaRXPwc65vGwo59zKiG1zE64zxDm3uuwWitRA/3a9536DL6c3JJzPUM45UTkRUYmVvS46AImtE3HxNfzorucLLlcuIo0nKSfiJTP7tJmVVNXSZ2ZpM/sKXkKmSEPwx/Mr8bzr6XX09uY/mMNyykWXcN/cKZ6ly5+NEd7vIhMrXdS59A2HhHMi/KBCIYRI80kKEM4ALgRmm9k0M9u0mAua2QQzOwevHPYPM9cRaSiVCCK2/fZdfOBnD+deN+JRGk6sLH6KZ/EVK2MukCO/ZwK6I3oinMsPjnLaEGr/nb+9kmX/vEk5ESJNKGkp8BvM7G7gYuAi4BIzexp4AngOWASsAYYDGwM7A+/BW4yrG7gO+J5zrnCauEidqdRwxn/nrwhfOP9eoW3lVKws53t++HeMmp0R2RNR4HcIt/75mY+ybvFqdUWINKFCS4EvAj5tZt8FTgGOBT7n7w4c6n9uzMYLOn7q148QaUT9+a25NzAMUFLeQMIUzyit6dzHe3jiReTsjIjEyl6XvyhXTmJlqCvFAseISHMpanaGc24+3toZ3zCzCcCewKZ4vRBLgDeBZ5xzr1SroSLV5M046Cui1J+zM8ofzki+blg6FVFlMum9c5FTPHNKXEdsi2u+hjNEmk8pdSIAcM7NAeYUPFCkgVz36KucOHlc9n21vjVHDwX0zc4oZSnwQtcNa03lpkDlr5URkRMR0xMRkRKRmyAa1cbCTRSRBlPWzAuRZhB8iP7judzUnaiHZEXuWcSjtPhiU8GfC1+3JZ3cExFZsTIqgTLifkllrzcaPpJ0xzBN8RRpQiX3RIg0i9xkwHAFxv574Pm3evr15XlDDonnBWdnFHF8Szq5JyJyFc/IxMr8UMi5+CmeXzv/53z/9mfVEyHShNQTIQNWTmZwXgXGKt0z4sLBYYT/vLas6OGM3z36auJ1w1rDORGFKlb25q+dkTI/YAidS3yxqWxipaIIkaajIEIGrKTehmr1REQmVoY3FhlF/OO5virzxQyTpMOzMwrWichPrGxJpTLHRSRWxjThlp9fxNIHr9XsDJEmpOEMGbByeyIs57VqszOiFuAKvS9rDbAi2tsSSqzsiRi+CO9/dXFu5ftUyjsufKwL9EWEczpe+u+TrF+kOhEizUg9ETJgBZ+h4VSE/nzeVaLXo6iciMwv6f+u+SW5Q0FFr+NfcxYDfTUmWlKpTGJl6P6BnojYKZ5FtFFEGktRQYSZbWNmQ0LbtjCzs8zsejO70Mx2qU4TRarjP68tzf4cfvBFlXWuhMgpnqH3xc7OCCpmAS4/sdJfxju8uFb4Ej29sK6rl6k7bsqglnTmXD+xMnQuSTkR1e3dEZHaSQwizGwXM3sCry7ElMD23YGnge8Dn8RbH+M/mQW3ROpeb6/j2F/+K/s+/OCuZU9EtepE9PVEZIKIgnUiHOu7exjUks62MZ2y2PyHvp6I6MxK5USINJ/YIMLMNgYewKtO+S9gYWZ7Cm9djNHA48BkYD/gIeByM3tXdZsssuFWrOvKeW+h1/781lyJezngNw+/wom/eTz2GL9OhJ8akTecETreCyJ6GdSSyu5LpyxTbCp/PCNbsTIUQ4zadHNahm6sngiRJpSUWHkaMAo43jl3Q2D7VGAXYB3wEefcmwBm9kHgZbw1No6rTnNFKmPZmlAQEc6JqNbsjKKGM8q5ruPsO59LPMavWFlsT0RPr6Ozu5e2llS23SmL7okI5kmEm//Vs3/CN//0jPohRJpQ0nDG4cCMUAABcFTm9W4/gABwzq0E7iIw7CFSr5at7Qpt6Z/hjMilwEObylnFs7hiU7nDGeGeiKgpn309EX3DGVH5FzmJlTE5EY++vLiIVopII0kKIsYBsyK2H4T3mfWPiH2v4y3MJVLXlq7pzHnvz1io+gJckT0RoZyIKq0Eng7NzsgraV1ET4R/jchVPPtSK3P2XfvjH7Bk+i85/Q9PFW6kiDSUpCAiDeR8XTOzTYGdMm/vizhnBLCq1EaY2Xgzu9XMlmT+XGdmm5R4jV3NrNPMflDq/WXgWbE2eTijmNkOlVKJiSClJC2mUtGzM/J7IsgmVvr8nojIVTxjeiJeffE5Ot/Wmn0izSgpJ+JV4J2hbUf4+5xzsyPOORCYW0oDzGw0cD/QBlyYadMZwK5mtrdzrjPp/Mw1WoBrgdZS7i0D1/ru/DUhgqo3nBG1MXmqZVHXLeKcYF4DRM/GCOrs7qXX4fVEZLalMzkR+cWm8q8vIs0vqSfir8Bh/mwLMxuElzTpgJvCB5vZiXi9FHeX2IZTgbHAIc65C51zPwKOBnYDTizyGt8Cdi7xvjKAhb+Fp7IVK7331UusLFyxMpzwWNx1izgGx7quHpas9uLy8Npa4cBgbVcPAINaUtlGpvyeiIjAJzs7I+++ItKskoKIi4DlwINmdh/wIrAr8DZwqX+Qmb3HzC4FrgaWAT8psQ0fBx5wzj3vb3DOTQdmZ/YlMrOJwHeBc0q8rwxg3aEnaL8twBW1LTyM0Ov42F5jN/i6Ufe58r7/Zd9f+PcX+PINMwP7c6+ytrMviPB7KVpSllNYqu/+LnY4Q0uAizSv2CDCObcIr/7D43jDFFsBM4H3OueWBA69Bfg6sBL4kHPu7WJvbmYjgQmZ64Y9CexV4PwW4BpgOnB9sfcV6QqtTplXbCrhwfe3Z97kmdeXV6wt4Xv1OMdFR++2QdeI0utc3tTWu55ZELhG7vF+T0RbICciboonLr5i5aZjx9M6akzB9olI40lcgMs59yJwYKbkdatzblnEYVfg9Vjc4JxbUeL9/U+W+RH73gSGmdlw51zcJ/Y3ge2AD6HFxKQE4eGMvC74hGfyl254EoC5FxwRf1CM6KXAw++TA4KtRnUwb8naxGtE3ci55B6L8H2DPRHhYlMu0DPR3ZtZfsv5wxm5/zU/dfq5vHrXC4m/k4g0pqLWznDOrY4JIHDOXeCc+3kZAQTA0Mzrmoh9/qfkkIh9mNnOeGW3T3fOvV7MzWbOnImZ5fyRgamrNzqxspR1HtZ395Rx5/wL+9/4fTFNy9p17Ii8bcXMJgkmPxbTsjWd3QAMak3ll73OHHPNp99FR2s6s55GRuifVbjXR0SqI/x8649nXNHf3jNTLo/EG37YGK8M9qvAHaUMYYT4QUziF6SItqTxhjEeds5dXea9ZQDrCT3YwjMKinko73XOdJ754aEl3Td82e6eXn58z4u5bSsw5zMd8cGQV/MhYNexw2lvSWf+lcUfF/6dV633gojBbRFTPLPVKQ2z5FU8b7jkOyyev5zR7/9q7L1FpDEVDCIywcOVwEfwakcEPyMc0GNmtwMnO+feKvH+KzOvHRH7OkLHBJ2BN3vjPZk1PgBGZl4HZ7Ytcc7lBCCTJk3iiSeeKLGJ0oy6wsMZZSRWrsw8ZDfEX595My9oKDQ7I7xsOeRWnwznR2w2rJ3V67vp7O5N7okI7Vu5zvv9Olpb8qZ4ZocuLPiB4G/LbeBbr79C19Lc4RcRqbyo3Khq90YkBhFm9k68hbU2wet5uAt4HliCN8ywE3AY8FFgipkd5JxLLuCf67XM6xYR+7YEljnnVkfsez9eXYmo1YbOyPwZT4k1K2Tg6AmNGWQX4Kp2xcrQ+7mL8kfyCi1DHlWHIRiIRBWiNCs8nBG+r19LY3Bbuq8GRIqcoYvsfy/ieyI0OUOkecUGEWbWDtyOt1rn94CLowo/ZYYWTsErFPVHM5vknCvqa4dzbpmZvYK3UmjYHkBct8Fp9PU8+DbDm6HxO7xVRheETxLxdW/A7IwNEb7ssrX5tdT8nohP7rM1f3zidTrzpqNGBBEuGETk1ZPEsMzDP/73itsTHM5oSXlJlv7QhzeckZsnkT/FM/aWItLgknoiPgVsD0xzzl0Rd5Bzrgf4sZmtBX6WOe9XJbThT8A0M9vBOfcCgJlNzdz74ph75k0JNbNxmR/nZOpMiMQK5xAUWyCpUL5CIVHlouPucd6HJ9LRmubXD7+Ss7/QcEZPr8u7blE9ETE7Bw/q+5jwi035v4Y/nJFTJyK8mFnouivXddHRmqYlXVRet4jUsaR/xScAL+HlQxTknPs58AKlLwN+Ed7wyL1mdqqZfRv4I17tiOsBzGyCmR1nZhNKvLZIjtP/8BQPvrgwr9jUrf+Zz5LVndnHX9zDtitc5nED9PY6HnxxYd724L2j2pGOiCKCPRFeWerc9+Hegsj2xOwc3BpIrPSTKDPvLfM/4TyJoE3GbU/bpn3/dCf+4B+ccnPU2n4i0miSgojtgPtcaf269+L1IBTNObcQb/nwp4CzgWnAbcBhzrn1mcOm4A1TaJlx2SB/nPk6J/7m8bzESoA3lvWNwsV9K0+aBVGM4GV/9fAcXlmUn/IT7O2IGn5IRQQRvTk5ES4nIOh1zustcPk9FLGNC+jIm50RXGyrr98hJ7AIOOikMxg19fM52/76zJsJDRGRRpE0nDESr4egFCvpq/1QtMxiXocn7L8Wb4GtpGvMJf/zSyQr+HAOT/EEaA10r8c9a8M9GKUKPqefiql6Ge5VCIsazugJBRHh2D87nJE4xbPv521GD+bVxWtoS6e8tTOy9/YGL/wgK2V+L0cgQAm1b/PhfZOvNnQ4SETqS1JPxAJK7FXAW/Xz1fKbI1I9waGIqGJTrenAZMWYb+U/u/+lDWpD8CG+bE30ArVJ0zUhZnZGsOehF557s6/2myOTt+CIjI78+wV7X/whk/EbD8lJ5Mz2RGTem+UHKOGciCeu/SGL7riEfcaPquhwkIjUXlIQ8RDwfjOLmn6Zx8zG4k33fLgSDROptO5Q8mFYyn8iEp8TcfVDr0TvKMPiVTFBRLAnImJ/VNsefXlxzvkP/W9RzvEpM+LmZvgBVXBfSyaI2HbT3IKx4bLXfsjwwoKVnPDrx7PHBL0xfz6Du5fT6xz3vVBuXToRqUdJQcTP8Ao+3Whmg5MuYmZDgT8AgygyEVOkv/lDESnLn+IJuQ/RpG7/YsT1ZAQ3v7Agqo5a7rBC1GWWRvRgTH++r85buFiVw4uNenujczp6Inoi/N6O4R1tOceGy15bZjjj8VeWZK89ZFCaMANeW7KGL2fWHRGR5pC0iudjwOXAAcBMMzvRzDYPHmNmo8zs03hJkXsDP3TO/beaDRYpl19vIZ0yuiOGM1xCLsLVM+Yw7sy/Fn2vStVGiErw3GXM8MRzwr0s3u/lLeEdldPhr20RvFVLZmgnOMQDXsXK4BTPlIUHL2DooNa8e5gZq9eXs9aIiNSzQhO1TwMuxct1+A0w38yWmdlrZrYQr4rlr4CxwNnOuXOq2lqRDeD3PqRTFtkT0euIneJ5+fQX845PEje7o9TgIurwER2t/OPr8ROVovIOvPUtXHJPRGBfOuV9NLSkcj8i0ikvAaKv2FT+lM6N2rWgrshAkRhEOM8ZwF54UyzfAIbhBQ2j8ZIvfwHs6Zz7YZXbKrJB/Ifruq5eXnw7fygh2BMRDgJKndpZqTkI0bMzLDK50hcOkFwgOIrKBfF7J6JyIvJ6IsILcPmFIgLCwxn77rsvm71jYpmrnopIPSt2KfD/OOdOcs5tBWwEjAEGO+fGOOe+oiEMaQTBJannLcmvzJ6Ti5C3r7SwYLvv/C1vm3OOk66JWu4lScR9LXqapy88VPOdI3b0EiudFwztttWI0PHxszNaQkFEKmUEJ3nkD2bAoJbcIOL8889n74+drCXBRZpQyf2Ozrk1QP6qQSJ1LioPIuhn97/ErHnLvDeh553/oE2nrOxaB+u6elm8OnpGRpyoJhfqiQg/rHfcYpiXWOkcPb292V4GX3dEToS/3HjecIblzs4ITGhJFLV8uYg0vqKCCDMbhjd9cz+8FTeHAauAOcAjwF2B6pIidamrO/nh/5en3sj+HO558N/uNnY4T762rKz7d3ZXpkZCyqJrRfii8j38Wg49vS5vCmZ2imfO1FLv56jhjLzZGQXa+9GPftQLzg4+tcCRItJoCgYRZvZN4JvAcKI/L04FlpnZOc65yyvcPpGKiSowFScu3Ijrkn9l0Wq2GN5Oe2s6tqeinJyAqKmmZt6S3HGifk9/Fc+eXkdrOsVGg1pYtb4bCE7x7Du+bzgjP7EyZ50MLNsTsf92G/PTT+YvyLt48WLWrVpF/pwNEWl0iTkRZnYNcB5e/Yfrgc8DRwBTgQ8DX8db52IIcKmZ/bqqrRXZAFHf0AGOmJhfT62UBbjWdvZw0CUPcOots2KPAVhfRk9EVDxSamIl9PVEdGd6InYZMyy7z29vVJ2I8NBHOmW5K3ZaX15ER2ua4R3RoYIGM0SaU2xPhJkdA5wIPAEc45yLK2f9EzPbGrgZOMnMbnPO3VH5popsmLiHe9SqmKUswOXXn3joxUU578PKCSLigpnEICKqJyKTWOkPZ/z0E3txwd+f56bH59HT67j1yddZtqYre3w6Ozsj93tGyqLLXgfPiaKUCJHmlNQT8RlgKXBkQgABgHPuNeAovAW4Pl255olUTlwQEfXsix/OyL+Gf74feHTFBAsbMpzxib23YuTg1sz9LHE4I9gT8bWD3wGQXcWzp9fRkjKGD27l0J292nEvvrWKU295KmdlzZaY2RktqdzFtoLFpqJWFxWR5pYUROwO3O2cK6rYvXNuMfB34F2VaJhIpcXlM0Q+/OJ6IhKmKfp7wvd5e8U6IL8n4sb/t0/stcLN2GubUWw1yqs+X3h2hnefCz86kVPf562hZ+bNDgkmVvq9DAtX5udE+9cP3yeVWYDrnDuf864L2QW64mZgHHLIIYybWPh3FZHGkxREjATmlXi9ucAmZbdGpIrilvGOevjFzeKMGqrwj832RISO8deLWN+Vuz0pEPAFp1L6rMDsDD+ICR4zb8kaFqxYxwsLVmanbfrBxIV/fyHvGn4PRHhYx7/igkxgFGxG3HDG9773PfY/9oux7RWRxpUURLQBq0u83lpQErbUp66YyCDqgRy3gFY4EPntI3Ozx/qnhAONJZkFs8LDGUk5BNl2BNqYHTYoUGyqJ1DTwjd/WV9xrVS2JyL+Iv5/k/B/svz/Vn2zM5ICm3CgphwJkeZQVMVKkWYQl6sQ9UCLG7QID2ec9Zdnsw9oP4gI90T49SHCwxnF9ET4D/HcnghLzD/wEyuD1w/OmvDzHdIJiRX+MeFgKtxks+D1oq912GGHcfM5X8q+/8TeW7FRm9bXEGkGCiJkwChUsTIobjijq7eXIW25ZZ1fWeR12PlJkOGiVn5QkR9EeK+Ttx0d246oHhEjOQDxg5ZgoJFTjdJPmkwIRPwAI1zzInxKyix7vbg2rV27lu5OL+/CDAa3teQtVy4ijanQ14GTzOzAEq43rvymiFRXZwlrN+QuC973c3eP413bjmLOwtXZHoi1nd4wRW/McEZXj2PFui4em7M4Z3s6ZTx/9vtpSVvkWhsQnErZ94BOmSUGANkS3YFzgg/tuJkXQf4x4WDKwsMS9CVoJvWOpFPQDbS3pGnZgNLhIlJfCgUR4yg9MNCng9SluMTKQlas7e67Rq8jZcbx+27DBX/zEhLXdXlBRE+vY/Gq9XnDGV3dvXzmmn/zxKtLc7anzOgI9Wrk8YczApvMoL01za5jh/P068vzTulb5yNwmcieiPiOyFRoOGPqjpsy/fm3o4cz0smzMyA428Of4aGPCZFmkBREjO+3Voj0g6TpmWHBZ9wby3NX/EynLOeBed8LfbOgD718Bpcfu0fO8Z09vX0LewUUNTuDwOyM0LTLzYe18zTLaU1bzrRSP1gK9hoEv/mni0is9AMQ/7yfHzeJNZ093PjYaznHGZY32yP6en370qaeCJFmERtEFCowJVLPnnxtKQ+9uIhTpm6X3RZXSTLqWR78prymM3dWhf9t2nf/7IXZnxet6sxPrOzpjSz7XNTsjGxPhGXf+O31Xw/eYVPufvat7Dn+/dMFgoik+5+47zgeeXkxH9pjDOANWQzvSOXlRJj1BSNxQdGRRx7J068vYwZex0rK4nNORKSxFLMA11bAocBoYD7wt0xhKZG69ZH/ewQgJ4gopScimEMQXn3TzBKnWIaDleCql0FxsxnC50JuQmM4LyG8SJbfK5HOSazMz4kIl7QO2mrUYO477cC87dGzM/yeiOhrnX766dz1zJvMyNTL8AOw3l6nKpciDS4xiDCzs4EzgeDA7VozO805d1VVWyZSYaXMzugJBBzhoCAdmJEQZVmmLkQh4WAgSq+LGs7InJ/p32gNtcXvdQhePneFzsLDD3HCvQ1mls2JSAoIBmdyP5zr6yHpcY6UluYSaWixX0XM7FPAd4H1wI3AxcCf8YpQ/Z+ZHdIvLRSpkLjhjCg/e+AlHn9liXdeeGpmKjmf4dXFa4q6R1Iioq/v2d83nBG+d35PRGY4I/BQj5qd0ZqUWFlkNajg7Iy43+fAAw/ktBM/0nftzP2VFyHS+JI6VD8HLAN2d84d75w70zl3NDAZL7D4Sn80UKQS/j13SV4wkGRdVy8fu+pRIHo4I+lb/Lyla2P3AbS35vcEbLvJkOiFwKKGM7Lt8F7DCZJ/++8C7/qBh3pv4IHtP8TTCYmVcTFEfk9EsNhU/PXaW1OMGdHB+R+ZmD1OMzREGl/ScMZE4Bbn3MvBjc65J8zsTrxgQqSuvfT2Su5+9i0uvnt22dcIJ0qmzRJ7Ee546o3E62WnOwYeuvd8/QAcsO2378o5tm/tDMs+2S37mmlPttiTN2zx0tur8q4fWSciaenu2Lbn/y7+cEbS8ExrOsUDZx4MwNUz5nhtUk+ESMNLCiKGAXEreL4IfKDyzRGprKk/nrFB5zvnmPb7WTnbUgatLeWP5adCuQ0Qn0+QLTYV2BZ+VvuJjaOGDGLRqr4VOVMxPRHFVKyMCwiiik359ylmeAaCiZVFHS4idSxpOKMFr8hclC600JY0jfiH38sL89egS6WMtnSBIlEJ/G78onIi/DyIwL9U/6GdTazM9AS0pY2h7X3fCwoVm0oafii2JwILBBFFFtH3R1FU+lqk8WkVHJEEh//kobxtKTPaWspfdqY3GxgUMzvDe7XAYz18mp9YaWaMHtLGynXd2Xb6gg/s9pa+4+PE7srrieib7hr3+3zsYx/LeZ9WYqVI01AQIZIgakZHypKrPRbSlyxZ+Bp7jx/Fgy8uZMsRHdlt4dOCUzyDq3XmzM4IPLAHD4r/Zz9ycCtL13TFBhj5ORGFhzO+/OUv556jxEqRplEoiNjdzE6I2g5gZscT0fPpnLuuAm0TqUvpVGV6IooZzvjSAdtyxMQtGLfxkOy27APect+bwdD2vlRlrdMAACAASURBVCAiJyci8MAekrAM94Uf3ZVDdtwsdr+F/rmbWfb+ccMja9Z4U14HDx7sHWfqiRBpFoWCiA9m/oT5nxbXRmx3gIIIaVpmxqANCiK814QyDVmplOUEEBCY4hlx/LCOvn/SuUFE3zEdbX03/vHHduPUW57Kvt9oUEtirkRe2evAtrielcMPPxyABx54wDtOwxkiTSMpiPhhv7VCpEIWrlzP7bPmV/Ue3nBG+UGEiykaVfT5MdtTZgwd1NcTEex9CP4cTAptb81NEB3UmpwwGlX2OlWgJyLM74nQcIZI40tagEtBhDScr970JP+as6Sq90hvcGJl5jollp3OmxgRqhdhltsTEcznCD6vg/kc+YFM8oM9nCJiWLZHpdh1MJRYKdI8lFgpTWX52rhZyfGK6RAY1t7CisysBzOjbQN6InxRPRHfeP/2bDa0vaTrjN5oEAAf3XNszvaumAqdwQDIrxWx+1Yj+PAeY9hz65GJ9wqvP2IpGNTi9V6sC612GkeJlSLNQ0GENJXy50wk+9OXJvPey7zCVRs6xdMX9cX9ywe+o+jz/dOHDmrhhXPez6CWFNc/9lp2f9wjevSQQXnbRg5u5cTJ4wreM68EOPCVg97B60vX8p7tNi6q3f7vrY4IkcanIEKaSjlpBsV8IQ4GDekUFemJKGcVzTh+bsPGQ9qy2/YeNyrvuD9+cV+2Hj04+7478yQPL+IVJzwE0dsLmwwdxK9O3Cv2nJNOOinnvWZniDQPBRHSVMrMVSwoGERUqieimKXAk8/3Xl2gz2HMSK+exEf2GJOTo7DH1iP4z2vL2CsUWPjDE8XWvegOPfiLqToZDiI0O0OkeWz4J6FIHQnXMUgydmQH+4wfxYmTt8lu++b7d4g8NtjzYBsYRHxmv/FlnxsU9ZvuOnYElx+7O2d9YOec7dd9Zm/+8fUpecd392R6IoqZb0ruYmSD29I5ZbbjLFq0iEWLFmXfa3aGSPNQT4Q0lVK+3I/eaBC//8K+LA4sWrXD5kOzVRuDwsMZGzLF8/tH7cT3j9qp7PML+dAeY/K2DW1vzSlE5fODgpZieyIyQceglhTPnf3+os45+uijgb46EZqdIdI81BMhTaWcAYLgLIm4aYrBoCFlRkvKOHavrcq424YLf4HfkC/0/oO8tcSeiGlT31n2PTU7Q6R5KIiQ5lJGnkFOEGHRuQrBY1pSKcyMC4/ele03G1peOyvAb+eGPIv9aph7bjOiqOO7Mj0RG7J2SF9iZdmXEJE6oeEMaSolPdoyT1+LWGY7LDiTIji04QoUZ6omv0Ub0oJ3TxjNfacdwPhQae04fiJmywbMLPE7PTScIdL4FERIUylnwkNuT4Rly1LnHtP3c04Q0Y/PwbzfrUIzUSZsslHRx/o9EcVOCY2ixEqR5qEgQppKOetRBAOEuC/YwSGOYBCR9G36l8dP4rZZ87nrmQUAHLbL5owY3BZ7fCH18MwdMdhLzhxZwu/xpS99Kee9EitFmoeCCGkq4RggZYUrIwanhaZTVjAQCa7g2RkzsL/VqA7et/Pm3BZYDOznx01KbkiZonpOquWUQ7Zjq5GDOXzi5kWfc+yxx+a8z9aJqIeoSEQ2SF0kVprZeDO71cyWZP5cZ2abFHHeoWb2kJmtMbNVZjbdzN7dH22W+hR+/vvrOkTxH2HBc1IpY3hH/lTI3Gv2/bPxpzzmXdvlvlZC3gqaVSvyHa+9Nc0n99m6pEJZ8+bNY968edn3fs2NcAltEWk8NQ8izGw0cD/wbuBC4FLgA8A9ZhbbZ2pmBwB/A0YA38Fbunxb4EEz27va7Zb6FH6wtrfG/xX3H/DBnodBLSlGDknuqg8WnuoqcorB/31qz6KOS3L0JG+Bra1GdeRsr/fv88cffzzHH3989r1fontdV3ELdolI/aqH4YxTgbHAROfc8wBm9hhwD3AicHXMeZcD84B9nHNrMuddBzwP/Ah4b5XbLfUo9AXZe2B1RR7qC+ZBDGpJM3Jwck9EMCciHETsvOUwnn1jBbuOHZ65tkU1qyyf3HtrPrn31nlLgNd9FBHS0eYFEeu71BMh0ujqIYj4OPCAH0AAOOemm9nszL68IMLMRgK7AZf6AUTmvLfM7EHgfdVvttSj8MN6UKg8dWvasjMM/OmZwemb7a0pPrj7GKY//3bsPXKDiNwn+An7bsM7NxvKDpsPA+CsD+zE8MGtHLLjZiX/LmHhIYT+H8yojPbMf7+16okQaXg1Hc7IBAMTgJkRu58E4pYGXAFsD1wWsW9joLsiDZSGUygn4sEzDuLsD+auK2E5wxlpjtpty8R7BK/p103wjRkxmD22Hpn9tr3p0HbO+/DEiizYFaeWtSrKoeEMkeZR654Iv8j//Ih9bwLDzGy4c255cIdzrgf4X/gEM9sV2A+4u9INlcYQzokYFMqJ2HJEB3tsNTL2/KQcCl9cT8Rlx+7Ge7bbuNimbrBqrVhabX4QoZ4IkcZX6yDCrxm8JmLf2szrEGB5xP4cZrYRcF3m7QVRx8ycOTOvS7g/p8dJ9eX3RPQ98PfPPOD92Rc7bTEs7/yk2Ry+tphCS7tsObzYZlZUvf8VPu2003Lep1NGWzrFOuVEiFRUKbOmKqXWQYT/aZz0MVjwk8bMBgN/wcuTON8592AF2iYNKG8aZGDD7z67DwBbjx7Mn740mZ23zA8iotaEuOAjE3PeV3NoohR+r0udxxAcddRRedvaW1MazhBpArX+NFyZee2I2NcROiaSmY0A/gEcBPwGb7pnpEmTJuGcy/kjzWXB8nU57/2cyU2GDsrZPmmbkdlu9aCoSP7je2+d8z7Yu7HvhNHZn/v7b1OjDGfMnj2b2bNn52xrb00riBCpsPDzrT+ecbUOIl7LvG4RsW9LYJlzbnXcyWa2KV6Nif2AXwKfc4oMBqz/vbWSlxfm/nVJmfH3afsz/esHVOw+wSDi1yftVVQeRTUc9+5tADjgnQXrstXUF77wBb7whS/kbOtoSysnQqQJ1DSIcM4tA14Boirx7AE8EXeumQ3FS6DcHbjMOfcFBRAD2wsL8jutUmbssPkwhheo/RB251ffE7svOJwxuK2FrUYOLunalbLLmOHMveAIthwR1ZFX39pb1BMh0gxq3RMB8Cdgqpnt4G8ws6l4UzhvTjjvZ3gBxE+cc6dWt4nSCFavz5/ZW26X/3abxa9sGZcToRC2eO1taSVWijSBWidWAlwEnADca2aXAu3AN/BqR1wPYGYTgMnAI865OWa2I3A83qyNWWZ2XPiizrnr+6n9UidWRQYRxUURt31lPxatXJ99n044Lzw7wz+00eo11FJ7S0rDGSJNoOZBhHNuoZlNwSscdTbedM/bgDOcc/6n+hTgGuDTwBzAH+AentkeRUHEALNsTX5567ilvcN232pEzvt0wokt4SCiYWtH1k57a5plazpr3QwR2UA1DyIAnHOzgcMT9l8LXBt4/wvgF1VvmDSUqDH2Qst6xyllvvXOY4Yx+62VbDSoLv451Z3vfve7eds6WtO8qZ4IkYanTz1pGp0RK2oW2xOxIc778EQ+tc/WjK1RgmW9mzp1at42r06EciJEGl09JFaKVERnd/5DqT8quLW3ppm0zaiq36dRzZo1i1mzZuVs0xRPkeagnghpGlFBRH/0REiyadOmAfDAAw9ktw3SFE+RpqCeCGka6yOHMxRF1KOONgURIs1AQYQ0jc7uXnbYfCjfPWLH7DYFEfWpvSVNV4+jOyLwizLuzL9y/l3PV7lVIlIqBRHSNDq7exnUkuKz7xnP/9t/vLdRMURd6mjzPnrWRQxBxblqxpxqNUdEyqQgQppGZ3cvbS0pzIz37+ItxzJxTOWW5/7+kTtxxqHbV+x6A5m/+JmGNEQamxIrpaEtXd1JW0uKIYNa6OzpzS6GNWmbkdz51few0xb5y32X6zPvGV+xaw0k5513Xt629pbigwgtiSNSv9QTIQ1tj3Pu4eBLHwAyPRGBapK7jBlOStMzam7y5MlMnjw5Z1t7mxdE/P7f8wqe36sYQqRuKYiQhvfWCq86uj+cIfXlkUce4ZFHHsnZ1p75/+nK+14qeH6PogiRuqXhDGkK85asYfZbK3nn5kNr3RQJ+fa3vw3k1onoyPREhD38v0WMGdnB+I2HZLcpiBCpXwoipCnsf9H9QP4KmxvinA/uzKghgyp2PenjJ1aGHffrxwCYe8ER2W09yokQqVsKIqSpDI75hluO4/cdV7FrSa6OmCAiSk+PggiReqUBZGkqQ7SSZkMoJXdFPREi9UtBhDSVtrRmYzSCUvIcunvLW+1zyepOjrjiIV5bvKas80WkMH1tk+aiMtd15/LLL8/bts3o4pdNLzOG4M6n3+DZN1Zw9UNzOOdDu5R3ERFJpCBCmorKQtSf3XffPW/b4LYWTth3G/7y1BvZbXFFpcrtiRCR6tNwhjQV02IZdWf69OlMnz49b3trOkV3IGmyO2aIo1AMsWp9Nx/82T95YcGKnO1KpRCpPgUR0lQ0mlF/zj33XM4999y87S1pozOwimd3zCyMQomV/3p5MU/NW8ZFf5+ds13lskWqT0GENKzeiG+uB++waQ1aIuVoS6foCgYRMV0OPQW6InozwULcUJYCS5HqUU6ENKzwN9RggSKpfy2pFM55MzXSKYvviSgwnOHHkhYTLSiGEKke9URIwwpOExyq+hANp7XFe7z7vRFdMT0OhRIrC/VEiEj1KIiQhhVMxNtkqMpTNxq/RLkfRMT1RBRKrOwLIhRFiPQ3fX2ThhXsibj4mF1r2BJJctVVV0Vub0n5PRHe/49xQUShngg/CNGy7yL9Tz0R0rD8IOKHH9iZSduMqnFrJM7222/P9ttvn7e9tSW3JyJuOKO3wCyLdV2ZIKIGPRHL13Zxzp3Psb67p9/vLVIPFERIw1q6phPQN9B6d8cdd3DHHXfkbW9NFTecEbfdt67Le4CH/xr0xwTPy+55kV8//AqfufbfPPy/Rf1wR5H6ouEMaViHXjYD6OsWl/p06aWXAnDUUUflbO/IrLi6an03QM50T4D5y9byysLVpAp81VmbDSJy/x74PVVxszYqwa9z8c+XFvPPlxZnZwgde9WjpMy46fPvrtq9ReqBeiKkYfmJlWkFEQ1p61He+hn+AlnhipVfv3kWx/36MZ5+fXnidf5/e2ceJlV1JfDf6Y1egG5oFqFBNokigiCoiQiikmR0xi2O0TiJEpcxMZMJIXEyMYxZjGMckxhHIiaOBiWJYlyzkChEcUEDKqKskVWg7aaBht73PvPHe1VdXdRe1V2v2vP7vvpe133n3Hdun6p659177r2+4YzgT0FbgluIf3i0id+u3RvTYlXhPnlrd1fzxq7DCV3fMDIJ64kwMp5sy8rPSMaWFgGw6NlN1DS1UVZS4D/X2NrOuj3VALyxM/LNuLKmCaDb6pcA7dEWmAjBi9sOcN3StwCYdUIpY1wbwxHc+9He0UljW2z5EfUt7by3/yhnTRgSt52G4RWsJ8LIeOLZVtrwDsWFuZx30jCq6lq45cn3WLB8g//cC5sP+P/eWlEbSt3P5g+d8y3t3YOGtgQ+F74AAuCdvUejygfHr7XN7Ww/UO9/7xuqCcXXl2/g6gfXUlXXHLedhuEVLIgwMp6WBJ44DW+w+Orp/M/lU7njslOobW7zl9/5560cP7iQT548nKq6loh1+IZDapraupX7eiKCcy3C0RoUhKzfeySqTnAfWH1zO9sP1Pnflx9pCqu7qdwZpmlps8+vkbnYcIaR8QT/+BveYtmyZWHPFebl8NnTRwMwpayYixevAZxhgl98YQaNre2s3HIgrH5reyd17tP+ut3VbNxfQ2NrO8vf2sfT68uB6LM7AL72+Ds8t+HDbmXvRsnFgGOTNutb2tlb3eh/v6OqnhOPGxBS19eD1hTj8EeyNLa28/6BeqaWFduMJiNlWBBhZDw2R9/bjB49Oia5qaNK2H3nhew/0kRp/zwK85yfpyH9+3GoPnRvxFF3mu+8ScNZtfUAFy1+DYDc7K6bZLj1J3w0t3X4A4gl/3IaF0wZwXee2ciKjRVRbQ4ezqhvaedgXQuDCnM52tTGjqr60Ip0rX/REDDkcbCuheqG1rCBRyg+ONxA+dEmzpowhJb2Dpau2cP8WWPpl5PdrY3n/+RlKmqaWXb9GcyeODTm+g0jEhZEGBmPdQd7m+XLlwNw5ZVXRpUVEUa7szZ8rFo4h2k/WMnwgf14Z+8RTh450H+DPNLoDGFcNr2Mr50/kYVPbODS6WVcMm0kZ9/1EnBsT4Rvw69tlbXc9uxmmt0g9K7Lp3DBlBGAM3PkSGMbtc1tDMzP7ab/l02VbKusZcG8jx2zk2xDSztVdS2MHlyIHGmKmO/g64loaHGur6qcfscqAAYV5nLmuFKuPvN4DtQ2c9GpI8nP7QoKOjvV35twzt2rAdh2+z/wyOt7uPPP28jNzuK6s8f55W99ZiMVNY4t+yMMsRhGvFgQYWQkRxpa/X8HJ9QZ3mLJkiVAbEFEKEoK8/jcGaN5bN0+Lrv/db7xyY/x1fMnArDRzSso7Z/HlFHFrFx4jl9v1cI5zPvpK2ytqGVrRS2lRXms2XmIhU+8y8jiAipqmvw7gN44exyfndnVY+KbfrqvupHJI4u72XPL796lrqWdT5183DGfvS/9+m1aOzo578RhNLV2cLi+tdv5to5Oct09Q3zX/vxDa5k6qpjJIwf65Y40tvGXzZX8ZXMlAHsON9DRCQvmTWRTeQ1X/fJv/PiKUxk9uGtGy/oPjviDqsDejfKjTfzpva5elYNRckwioao9uu6GkXlYEGFkJLsOdXUTX/OJMWm0xOgNAqdSrttTzWPr9vL85kpe236IU0cVc9rxg47ROWHYAO658lQWPbOJC+591V9+3MB8ZowZRGVNATnZwp2fmXLMVM7RAWtYBAcRvp6L7VV1/iDi+QVz+PTPXqGlvZMThw/gurPHcd+L2znc0HXD3lfdyKfueYXZE4cwsqSgWyLoe/trjlkPY9SgAoYPzGfXwXp+/tJOAM4cP5jdBxto71QWLN/AvEnD/fJbKmr9QySNbR1s3F9DTraw+MUdALz4jXM47ycv89etB7hk2kh+v+FDPqxpoigvhytmjuZgXQuDi/I4OSCYAXh1+0EeW7eXuy6fyjUPr6N/vxx+Nf90crItL9+wIMLIUN7d5/zgrvnP8xgZsL6A0TcZNiDf//er2w/x6vZDjB5cwPVnj+OG2ePJywl9Q7ts+ijOPmEoK7ccQFGyRZgxZhATh0fOOTi+1Aki/vfFHTy7oZycrCyKC3P5p6kjyMnKoq2jg1+8vIstFbWcMKw/Y4cUMqI4n+tmjePGOeMB+O26vWxxp582tLTz2Lq9NLV18IKbKDphaBE7DzZ0u25ZSQGXzxjFjbPHMSA/F1VlW2UdNzzyFuVHm3h7z5FuOUCrth7gK+dO4Mm393PPyvdpaHXOLVm9kyWrd/rl/nXOeMYP7c/Y0kLe3V/jHwLx8X+v7QagMC+bZdefyYwxXUHZFx5aB8CKjZX+sh/+aStlJQWUDSrgglOOs96JjzASy6psfYWZM2fqW2+9FV3Q8CxPr9/PB4cb+cN7TiLci9+Ym16DjKjMnTsXgNWrVydcR01TGz994e9cMr2Mg3UtTBjanwlDi3r05nX90jfZfbiBvOws2jo6j7nhZ2cJk0YMYP5Z4/jnGaOO0b/rL9t48JVdLL56Orf/cSvlR5uYUlbMnZ+ZQmNrB2eMG8yyN/bwD6eMID83i3N//DJ3XT6F8wN6FwK5ePFr/t6K0qI8Zo4dxMiSAm759IlsrajlybfLeXNPNTuq6hk/tIjPTC9jwtD+DMjP5awJpWRlCe8fqOO+F3cwcVh/ph9fwtjSIvJysnhuQzkFudk88PIuKmqa/EHZgPxc//DH3BOHcsWM0Sx9fTdv7uma/vqdCyf5AyfDe4jI26o6s8fqtyDC8AptHZ1U1bVQVlJAZ6dyqL6Fv26rorW9k1+t2c35k4bzkPvElJedxb1XTfMnwhneJRVBhBfYUVXHjqp68nOz+fj40m6JjqF4bkM5X3vcWUCrX04WM8YM4vsXT47aCxKOptYOnn5nP0+9vZ9rPjGWS6eXJVRPJA7WtfDI63to6+hEgbrmNvJzs/n38yYyqCgPcHpVtlTUMmFof/7zqfd4YcsBhvTPY+KwAcz+2BBOHzuYUYMKOFzfyoSh/f17pBjpwYKIFOK1IGL3oQb/nPJAP3TziAb+GSATWB7kwsC34erVMBfpVm8Y+bB2hLluMDVNbeRmZ1GYl01bh9La3klLuzPNbsO+o+TlZNHRqSFXohw1qIBVbvJctB9xwxscOuTsbjlkyEdreeeaxja+/4fNnHPiUP5xyog+mUNQ39LOw6/t5sOjTby7v+aY1UVHDSrg/JOGkZ2VRU62kCVCTpaQ7b6q6ppRhaJ+ORTmZVOUl0NRvxyK+mVTmJdDbrZk1FDJ1LJif7DlFSyISCFeCyJ+8sLfuc9NejK6+KepIxhTWsiwAflMHVVMZU0za3dX09LeyU1zxjN2SOT9DAzDSA9Vtc1srqil/EgTLe2dPPPOfvZVN9HRqbR3dtLZiXN0bzv5uVkU5eXQ0Nru30gtk/nNDWcy6wRvBcsWRKQQrwUR5UebqKzpmkceGHAHxt6BkXj38u71ScDZcMF792uElo9JJqwd3S8cSse3iFB9Szt52Vn0y80iLzuLgrxs613ogyxduhSA+fPnp9UOwzuoOj2NWSL+9S46OpXG1nYaWztoaHGOmTZ9e+Lw/sesK5JuejqIsNkZaaSspKDbzoWG0RexIMIIRkTIye7+wJGdJQzIz2WAx27CRmT63iCdYRiGYRi9ggURhmEYhmEkhAURhmEYhmEkhAURhmEYhmEkhCeCCBEZJyJPi0i1+3pURKLuVZuonmEYvceKFStYsWJFus0wDKMHSPvsDBEpBV4C8oC7cGy6BZgqImeoamsq9QzD6F0KCwujCxmGkZGkPYgAFgKjgCmquhVARNYCK4FrgQdTrGcYRi9y//33A3DzzTen2RLDMFJN2hebEpGdwG5VnRdUvg0oV9XzU6EnIgqRl2MOkKWvynrFDi/IesUOL8h6xQ4vyHrFDi/IesWOTJP1ih0Bsj22dnhacyJEZBAwHng7xOn1QMhVthLVMwzDMAwjdaQ7sdK3DV15iHMVwEARKU6hnmEYhmEYKSLdQYRvT9zGEOea3GOo3ZYS1UNEur0MwzAMoy+QjvtbuoMI3/UjDe6E2oElUT3DMAzDMFJEumdn1LnHULtQFQTJJKWnqiIiB4EPAssjRWvxRHKZJusVO7wg6xU7vCDrFTu8IOsVO7wg6xU7Mk3WI3aMibnSBEh3ELHXPY4IcW4kcFRVG1Klp6q2EJVhGIZhpIi0Dmeo6lFgN3BaiNPTgbdSqWcYhmEYRupId04EwFPAPBE5yVcgIvOAE4HHe0DPMAzDMIwU4IUg4n+AauCvIrJQRG4FnsRZA+LXACIyXkQ+LyLjo+g9jZNsuZcQiMiXRWSbiDSJyEYRuSpWI5PRTSUislRENMJrdRT9CyPontJLzYhIsjZ6xVehEJGBInKviOwXkVYR2SMid4hIXgy6nvOd9PF9b0Tk0yLyqog0iki9iKwSkY/HoLcujJ+e7A27YyVRO73uPxEZG+V3UkVkbgR9T/tPRB4M9VufjF8S1U13TgSqelBE5gD3AD/Ambb5LHCLqra4YnOAXwFfBHaF0WsC2gEB2oKvIyLfBO4GfufqXAY8JiKqqssj2ZiMbg/wC2BViPLLgUuBP0TRn4wTaH0R6Ag6ty9p61JDwjZ6zFfBtglOoDsX+CXwHvAJ4NvAyTi2RsJTvpM+vu+NiJwD/BnYDHwHx86bgZdFZLaqrgujJ8AknN+xp4JOf3CsRnpI1M4M8d9B4AshyguA+4Aq4N1Qil73n4hcD9wAvBxUnrBfkvKpqmb8C8fhW3F+YBWYG3S+BKgHfhNQlgW8gtNrkR2h7oR1e7H9xwM1wPO4S5lHkH0E2JNum3vCRq/7CrjE/Xx+Jaj8Drd8Vib5zrW7HZgUUDbPbcuNqdZLQ/vewblpFAaUDcfpAV0ZQW+c25b56W5DlPYlZGem+C+M7T/DCcBnZ5r/gGzgNpzlCxRYnSq/JKPrheGMpBCRa4ENQCnwf2HELsJZfGqJr0BVO4H7gdHAWREukYxub/ETIB+4WV3vR2AyTsDlZRK10eu+museHwkqf8I9RrPPa767CueHzG+Tqq4C/u6eS7VeryHO0vqnAk+oqn9RO1U9gPMEGMlXk92jl3wVikTt9Lz/QiEiU4CvAktV9dUIop7zn4jk42zp8H1gGaFXa07GLwnrZnwQAUzFSaScAqwJI+PbS2N9UPn6oPOp1u1xROQ0nKGMxaq6M4qsr5tui/s+X0TSPqQVSJI2etpXwO3AdFWtDyof4h7bwyl6zXfS9/e9qcVJ0r4nxLkhRPAVQTchEQm5eq4HiNvODPJfKO7AGfZeFEXOi/7LBwYCV6rqtQR9/pLxS7I+7QtBxLdV9Vr3CSEcZcCRwCcKlwr3eHwP6fYGi4AW4EcxyI4HCoEJIrIRJ/+kQUQeF5EhkVV7jWRs9LSvVLVaVTeEOPVl9/h6BHWv+a5P73ujqh2qul1VPwwsF5GpwCwi++oUnMXufioidUC9iOwUDyX4uiRiZ0b4LxjXbxcBD6hqRRRxL/qvFpioqk+EOZ+MX5LyqaeeQn2IyHFRROp9T3MaWxLPABLYZyMFujETT5sDdMpwvhhLVfVgDJfxRdgfxwk6dgNnAwuAyW4CTVM45WSIo33J2NgrvgomEd8F6H4BpyfpJVVdG6GOtPkuDLHuX1OTIr20IyL9gUfdt5GC9sk47SwBrnGPX8NJ8M1V1WU9amjsJGJnpvrvyzi5EPfFIOs5/7nDspG2ckjGL0n5kcl1bAAACXtJREFU1JNBBF1PjuG4g+hdUoFkEXqfDV9ZJOckoxsPibT5izg+jOWLAc7Mlh/gJB6+75Y9JyI7gAeA64HFMdYVL7G2Lxkbe8tXwST0eRWRi4GHgEpgfpQ60um7UHyk9r0RkULg9zh5Eneq6ssRxH+Jk8T78wD9x4FNwN0i8ltVDZ5dkw4SsTPj/CciBcDngd+raiyzKzLFf4Ek45ekfOrVIOLGKOffibO+OkLvs1EYcL4ndOMhkTZfDLyvqu/FcgFV3YTzRQjmYZwb0Ln03I0opvYlaWNv+SqYuH0nIp/DSbCsAy5Q1ZBrm/hIs+9C0Wv73qQbESkB/ogzjPEwznTPsKjqAyHKmkRkGfBdnOm8G3vA1LhI0M6M8x/Od6M/zrTvqGSK/4JIxi9J+dSTQYSqhptlkSh7gcEi0k+71p4AZ58NCD0WlArdmIm3zSIyDCfh5a4UXLtNRI7gfNF6hGR9GqONveKrELbF67ubcGaMVAPzVDXkfPUYr93jvgtDr+57ky7c79nzwDScJ9QvxTADKhxV7rG3fRUvkezMKP+5XAi0An9Ksh4v+y8ZvyTl076QWBkL63EWoZoWVD7dPb7ZQ7o9ySdw7PprrAoicruI7BKRgUHlg4GhOOPsaSVJG73qKz8icg3O8EMlMCfWAMJrvtOPwL43IjKArgDiHlW9KVoAISJlIrJZRG4Lcdq3RL8XvmcJ2ZlJ/gtgFvCmqtZGE8wU/wWTjF+S9mm0BS4y6YUzrhxqsalSnASRXweU+RYh2kPkxaYS1u3htn7XbeuQOHRudHW+HlT+M7d8jgd8mLCNXvVVgC2TgGacJ5qJvfV/6cH23I2zOuxJAWW+BWquT7VeGtr3qGvTz+LU24+zgujAgLLROBn2L6a7XcnamSn+c+3KxZm9dl9f8p/7e7Y6VX5JSjfd/4wU/2PnEyKIcM/5brq/xVky9M/u+yuC5C4FLk1Et5fbugxoiCLTrS3uF2qt+2G5Hydj+Sm3LQ+k23/x2pgpvgqwbblry6M4iV7Br6mZ5DucHpADOMNEC4FbgaM4Ty79XJnxbtvGx6OX7hdOwKeuXfND+StC+y51dTfhZPUvwul5qiZgRcB0v2KxM1P9F2DrBLeNt4Q5n5H+I3QQEZNfUu3TtP8zUvyPnU/4IEKAb+J02zThrJv+z2GcsycR3V5u6wqgIoYPWnBbBtPVnd6Cs6DK14GsdPsvXhszxVcBtlXStTR7qNf3Ms13OAsyrcBZbrwKWAoMDTjv+07Oj0cv3S/gS1F8pVHadwlO0NcMHAGeIeApzyuvaHZmqv8C7DzDtf+mMOcz0n+ECCJi9UuqfSqusmEYhmEYRlx8VBIrDcMwDMNIMRZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIZhGIaREBZEGIaRNkTkQRH5fQJ6z4uIisglUeSyReSAiNSLSJFb9kMReUNE7PfPMJLEvkSG0ccQkZXuDTbSa1Fv1xWi7nNx1vH/jwTUH3GPV0aR+yQwDHhSVRvcsrtwNiH6twSuaxhGALZ3hmH0MUTkMDAQuCOC2BOquqU36wqqNwdnA7HXVfXaeHRd/QKcXQezgGGq2hhGbhnOjoXnqurqgPIFwO04W7JXxnt9wzAcLIgwjD6EiEwAdgDvqOppXqkrRN1XAY8BM1X17QTreAi4Dvisqv4uxPlCnEDjEM62xxpwbgDObqj3quqtiVzfMAwbzjCMvsZM9/imx+oKZiGwLVQAISIDReRHIrJTRFpEpFxElojIsCBR35DGVWGucTHQH3hUg56WVLUO+CPwJV+uhGEY8WNBhGH0LTwfRLg9HKcDz4c4VwysAb4F7AbuBd4A/hVYJyIjAsRfdWUudHsWgvkXQOkKNoJ5HhgEfDqxlhiGkZNuAwzDSCm+G/9ZIjIqjMzdAUmGvVVXIHPd41shzv03cArwFVW931coIhcDz+EEFZ8FUFV1cx5uw+l1+E2AfClOcPCqqu4KY4fv+nOBp+Nsg2EYWE6EYfQZRESAI0BxBLFqVS3tqbpE5Ebg28DxwEOqelOIun8O3AzMUNX1AeU57jU/UNVTQui9BnwcGKyqtW6ZL2/jj6p6UYDsl4H7getU9Vdh2pgHNAF/U9VZEdppGEYYrCfCMPoOH8O56b+mqrN7uy4ROQlYAlwB/A2oCyPqy204FFR+Ik4OQ7aIfC+EXj6QDUzBGfJAVXeKyBrgUyJSoqpHXdmrgUbgyXD2qmqriNQCQ6M0zTCMMFgQYRh9B9/ww7vRBEVkGnBpcLmqfi/eugK4GNikqs9EkfP1bgRPyyxxjycB342gPzjo/SPALOAzwMMiMsZ9/2s3gTISDTh5EYZhJIAFEYbRd4jnxj+N0Dfq7yVQFyLyPjDR/VuBZ1X1sjDi1e6xmO69EfXucZmqXhPLdV2eAP4XZ+Gph4HPAQIsjUG3JMAewzDixGZnGEbfIeYbv6ouVVUJfiVSl8vZwPvAImAEEGkBqQr3OCSo/O9ACzDDzcnohogsEJFFbtJkYFtqcJIuzxOREpzhlL3AS5EMFpF8oAjYF0nOMIzwWBBhGH0Adx+IaUAnsDENddXiLCW9RlUrfYmPYdjkHicHFqpqM7AcOBlnHYlAm+YCP8ZZXOpIiDofwelZvQk4jRBrQ4TAl7wZz5CNYRgB2HCGYfQNJuEkJdYC3wrxIO/jl6r6YQ/UdQrO78mGGGxdgbN+w9k4ww+BfBM4C/ixu7nWWmAUTr5DG85si84Qdb6A08Nxm/s+3NoQgfhmZKyMQdYwjBDYFE/D6AOIyDVEv3F2AiXRkg0TqUtEbgAWqerYGO1dizPsMTY4KBCRwcCtwGU4AcRBnNket6tq2F4DEbkbJwiJaUaJiPwJJ2AZ4faCGIYRJxZEGIaRNCKyGBilqsfM+Agj79s741Oq2us9ASIyEidv4k5V/a/evr5h9BUsJ8IwjFQwjdiGMnw8gZOIeWPPmBOVLwLNOCtgGoaRIBZEGIaRFO5MiqnEEUS4QxgLgMtFZGpP2RYKdwbHQuD7qhq84JVhGHFgwxmGYaQNdzvvEap6YS9e807gPOAsVe3oresaRl/EggjDMAzDMBLChjMMwzAMw0gICyIMwzAMw0gICyIMwzAMw0gICyIMwzAMw0gICyIMwzAMw0gICyIMwzAMw0gICyIMwzAMw0gICyIMwzAMw0iI/weaeKTO6MxxvQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "#mpl.use('Agg') #silent mode\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "import sys\n",
    "\n",
    "#This is a modified script form the vaspkit mypdosplot.py. \n",
    "#This version contains the normalization\n",
    "#------------------ Functions ----------------------\n",
    "\n",
    "#use to normalized the y axis using  y' = a + (y-ymin)(b-a)/y.max-y.min  where a , b is the range you want the scale to be.\n",
    "def normalize(datarange):\n",
    "    \"\"\"this takes the passed data and normalize it in the range [a,b]\n",
    "    \"\"\"\n",
    "    #y_norm = (datarange - np.min(datarange)) / np.ptp(datarange) #(np.max(datarange)- np.min(datarange))\n",
    "    y_norm =  (datarange-np.min(datarange)) / (np.max(datarange)-np.min(datarange)) #np.ptp(datarange)\n",
    "    return np.nan_to_num(y_norm) # this convert NAN values to zero. If there are any.\n",
    "\n",
    "\n",
    "#------------------ FONT_setup ----------------------\n",
    "font = {'family' : 'arial',\n",
    "    'color'  : 'black',\n",
    "    'weight' : 'normal',\n",
    "    'size' : 20.0,\n",
    "    }\n",
    "\n",
    "#------------------- Data Read ----------------------\n",
    "pdosfile= \"PDOS_USER.dat\"\n",
    "\n",
    "with open(pdosfile,\"r\") as reader:\n",
    "    legend = reader.readline()\n",
    "legends=legend.split()[1:] #omits the 0th label\n",
    "legends=[i.replace(\"_\",\" \") for i in legends]\n",
    "legend_s=tuple(legends)\n",
    "datas = np.loadtxt(pdosfile,dtype=np.float64,skiprows=1)\n",
    "\n",
    "datanorm = normalize(datas) # normalize the Y axis the orbitals\n",
    "\n",
    "#--------------------- PLOTs ------------------------\n",
    "axe = plt.subplot(111)\n",
    "# Color methods! choose only one of the two methods\n",
    "\n",
    "\n",
    "axe.plot(datas[:,0], datanorm[:,1:],linewidth=1.5) #auto colors\n",
    "#axe.fill_between(datas[:,0], datanorm[:,8].min(), datanorm[:,8], facecolor='orange', alpha=0.5,linestyle='-' )\n",
    "#axe.plot(datas[:,0],datas[:,1:],linewidth=1.0) #auto colors\n",
    "\n",
    "\n",
    "\n",
    "axe.set_xlabel(r'${E}$-$E_{f}$ (eV)',fontdict=font)\n",
    "axe.set_ylabel(r'PDOS (states/eV)',fontdict=font)\n",
    "plt.yticks(fontsize=font['size']-2,fontname=font['family'])\n",
    "plt.xticks(fontsize=font['size']-2,fontname=font['family'])\n",
    "plt.legend(legend_s,loc='best', frameon=False)\n",
    "plt.xlim(( -10,  10)) # set y limits manually\n",
    "axe.axvline(x=0, color='k',linestyle='--',linewidth=1.5)\n",
    "span= 0.5\n",
    "axe.xaxis.set_ticks(np.arange(-10,10, span), minor=True)\n",
    "axe.tick_params(axis='both',  right=True, top= True, which='both', direction='in', width=2, size=5)\n",
    "\n",
    "#axe.axhline(y=0, color='k',linestyle='--',linewidth=1.5)\n",
    "leg = plt.gca().get_legend()\n",
    "ltext = leg.get_texts()\n",
    "plt.setp(ltext, fontsize=font['size'])\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches( 8, 6)\n",
    "name = pdosfile\n",
    "plt.show()\n",
    "#plt.savefig(name+'.png', dpi= 300, format='png', transparent=True, style='presentation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Choose the range of integration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   #Energy  Pd_d\n",
      "0 -8.31815   0.0\n",
      "1 -8.30115   0.0\n",
      "2 -8.28415   0.0\n",
      "3 -8.26815   0.0\n",
      "4 -8.25115   0.0\n",
      "New data frame\n",
      "   #Energy  Pd_d     E-Ef\n",
      "0 -8.31815   0.0 -8.31815\n",
      "1 -8.30115   0.0 -8.30115\n",
      "2 -8.28415   0.0 -8.28415\n",
      "3 -8.26815   0.0 -8.26815\n",
      "4 -8.25115   0.0 -8.25115\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "\"\"\"\n",
    "Here the data is obtained from the  PDOS_USER.dat file. \n",
    "I rewrote the code so that it uses pandas to do the filtering, instead of np.genfromtxt. \n",
    "\n",
    "The data for the projected density of states is usually in a format where the first column shows the energies and the other columns contains the projected density of states (pdos).\n",
    "column1name = Energy, column2name = sum of density of states, columnNthname =  bands(orbitals) s p d.\n",
    "Always check the documentation of the program used to get the density of states.\n",
    "\n",
    "The Energy column from PDOS_USER.dat it is aalready referenced to the fermi energy as (E-Ef). I used vaspkit and chose the option to have the Energy referenced to the Ef. However, not all the pdos data files will have the Energy referenced to the  fermi level. \n",
    "\n",
    "This script can be modified to accept fermi_energy as an argument.\n",
    "\"\"\"\n",
    "\n",
    "first_limit_energy_i = -4.8\n",
    "second_limit_energy_j = 1.0\n",
    "\n",
    "#load the data and convert it to a dataframe object. \n",
    "dataPDOS = pd.read_csv(\"PDOS_USER.dat\", sep=\" \", skipinitialspace=True) #skipinitialspace->Skip spaces after delimiter.\n",
    "\n",
    "print dataPDOS.head()\n",
    "\n",
    "def energy_minus_fermi_energy (df, fermi_energy=None):\n",
    "    \n",
    "    \"\"\"\n",
    "    Retrun the energy referenced to the fermi level.\n",
    "    Transforms the pandas column energy to energy minus the fermi level. E- Ef/  \n",
    "    \"\"\"\n",
    "    if fermi_energy == None:\n",
    "        return df\n",
    "    if fermi_energy != None:\n",
    "        return df - float(fermi_energy)\n",
    "\n",
    "\n",
    "#Next a new column is created using the apply method and the function to reference the  energy to the fermi energy. apply takes as argumnet itself, and then the other arguments. I prefred to pass the second argument of the main function as a keyword. If it is passed as argument, you have to pass it as a tuple like args=(fermi_energy,) ( see documentation). \n",
    "\n",
    "dataPDOS[\"E-Ef\"]= dataPDOS.iloc[:,0].apply(energy_minus_fermi_energy, fermi_energy= 0)\n",
    "\n",
    "#Here we filtered the data in the range of the energies we want to integrate for the d band and width calcualtion.\n",
    "\n",
    "pdos_filtered_data =  dataPDOS.loc[(dataPDOS[\"E-Ef\"] >= first_limit_energy_i) & (dataPDOS[\"E-Ef\"] <= second_limit_energy_j)]\n",
    "\n",
    "pdos = pdos_filtered_data[\"Pd_d\"]\n",
    "energy = pdos_filtered_data[\"E-Ef\"]\n",
    "#dataPDOS.iloc[:,0].apply(energy_minus_fermi_energy, args=(0.35,)) #passed as argument\n",
    "#dataPDOS.iloc[:,0].apply(energy_minus_fermi_energy, fermi_energy= 0.35) #  passed as a keyword\n",
    "\n",
    "print \"New data frame\\n\",dataPDOS.head()\n",
    "\n",
    "#print(pdos_filtered_data.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Calcualte the first and second momemnts ( center of the band and the width of the band resp.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d-band center = -1.87375847851 eV\n",
      "d-band center width = 1.37351172548 eV\n",
      "integral=303.37017444\n"
     ]
    }
   ],
   "source": [
    "#J. Chem. Phys.137, 084703 (2012)  \n",
    "# dband center  first moment of the projectedd-band density of states about the Fermi level (Ef)\n",
    "\n",
    "# d-band width (Wd) computed as the square-root of the second moment of projected d-band density of states about thed-band center (Ed).\n",
    "# You can obtain the p band  center and also calcualte the pband-width, in the same way.\n",
    "\n",
    "I = np.trapz(pdos, energy)\n",
    "center = np.trapz(pdos * energy, energy) / I\n",
    "W_centered_width =  np.sqrt(np.trapz(pdos * (energy - center)**2, energy) / I)\n",
    "\n",
    "print 'd-band center = %s eV\\nd-band center width = %s eV\\nintegral=%s' % (center, W_centered_width, I)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2.7.15 ('EnvPy2.7')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "19b9ec579bbeb4879300fb1fc3a220ab550e61f3ee6f5799a60444e61b93c703"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}