DEHam/notebooks/t_J_Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. t-J Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. t-J Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a. Basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $t-J$ Hamiltonian with a S=1/2 spin-model is constructed in a basis of $\\mathcal{R}^3$ given by the vector B as follows:\n",
    "\n",
    "$$\n",
    "B=\n",
    "  \\begin{bmatrix}\n",
    "     0 \\\\\n",
    "     \\uparrow \\\\\n",
    "     \\downarrow\n",
    "  \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### b.) Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hamiltonian is constructed using two types of operators.\n",
    "\n",
    "1. The particle creation and annihilation operators $\\hat{C}^\\dagger,\\hat{C}$\n",
    "2. The spin ladder operators $\\hat{S}^+, \\hat{S}^-, \\hat{S}^z$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix form of the five operators in the above basis is given as follows:\n",
    "\n",
    "$$\n",
    "\\hat{C}^{\\dagger}_{\\uparrow}= \\left( \\hat{C}_{\\uparrow} \\right)^\\dagger=\n",
    "  \\begin{bmatrix}\n",
    "     0 & 0 & 0 \\\\\n",
    "     1 & 0 & 0 \\\\\n",
    "     0 & 0 & 0\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{C}^{\\dagger}_{\\downarrow}= \\left( \\hat{C}_{\\downarrow} \\right)^\\dagger=\n",
    "  \\begin{bmatrix}\n",
    "     0 & 0 & 0 \\\\\n",
    "     0 & 0 & 0 \\\\\n",
    "     1 & 0 & 0\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{S}^{z}\\ \\ \\ \\  =\\ \\ \\ \\ \\ \n",
    "  \\begin{bmatrix}\n",
    "     0 & 0 & 0 \\\\\n",
    "     0 & \\frac{1}{2} & 0 \\\\\n",
    "     0 & 0 & -\\frac{1}{2}\n",
    "  \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\hat{S}^{+}_{\\downarrow}= \\left( \\hat{S}^- \\right)^\\dagger=\n",
    "  \\begin{bmatrix}\n",
    "     0 & 0 & 0 \\\\\n",
    "     0 & 0 & 0 \\\\\n",
    "     0 & 1 & 0\n",
    "  \\end{bmatrix}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### c.) Hamiltonian "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Nearest Neighbor coupling**\n",
    "\n",
    "The $t-J$ Hamiltonian is written in the above basis with the above five operators and the two parameters $t$ which governs the kinetic energy of the particles and the magnetic-coupling parameter $J$ which governs the coupling of the localized spins on neighboring sites.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\hat{H_{t-J}} = \\sum^{N-1}_i \\left\\{ t_{i,i+1} \\sum_{\\sigma} \\left(\\hat{C}^{\\dagger}_{\\sigma,i}\\cdot\\hat{C}_{\\sigma,i+1} + h.c. \\right) + J_{i,i+1} \\left( \\frac{\\hat{S}^+_{i}\\cdot\\hat{S}^-_{i+1} + \\hat{S}^-_{i}\\cdot\\hat{S}^+_{i+1}}{2} + \\hat{S}^z_i\\cdot\\hat{S}^z_{i+1} \\right)\\right\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$h.c.$ stands for Hermitian conjugate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. 2D Lattice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.) Labeling\n",
    "Here a 4x4 lattice is shown with periodic boundary conditions where each site and bond is labelled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOydd3gUVffHv7vZtE0DFBIISegEUulRlNASqjQRadIRxYpI8JWuiAKvNFGKSBGQJgi8NKWDkRK60gktgUBIQuqmbfb8/rim7yabPNmdyW/P53n2YXdm1v36zcyce889c6+CiAgMwzAMYyEopRbAMAzDMOaEAx/DMAxjUXDgYxiGYSwKDnwMwzCMRcGBj2EYhrEoOPAxDMMwFgUHPoZhGMai4MDHMAzDWBQc+BiGYRiLggMfwzAMY1Fw4GMYhmEsCg58DMMwjEXBgY9hGIaxKDjwMQzDMBYFBz6GYRjGouDAxzAMw1gUHPgYhmEYi4IDH8MwDGNRcOBjGIZhLAoOfAzDMIxFwYGPYRiGsSg48DEMwzAWBQc+hmEYxqLgwMcwDMNYFBz4GIZhGIuCAx/DMAxjUXDgYxiGYSwKDnwMwzCMRcGBr5Kg0WikllApYd/KB/tWPti3ygEHvkqCVquVWkKlhH0rH+xb+WDfKgcKIiKpRTAMwzCMueAeH8MwDGNRcOBjGIZhLAoOfAzDMIxFwYGPYRiGsSg48MmQK1eu5L3Pzs7G7Nmz0atXL3z++edcLl0CS5cuRVxcHADgzp07aNeuHapUqYI2bdrg77//llidfOnXrx82bNiA1NRUqaVUKu7evYtRo0Zh6tSpSE1NxdixY+Hr64s33ngD9+/fl1oeUwIc+GTIiBEj8t5/9tlnuHPnDiZOnIj09HS888470gmTOcuWLcOLL74IAPjoo48wYcIEJCYmYu7cuexbCZw5cwY7d+6Ep6cnBgwYgN9++w1ZWVlSy5I9I0aMQKtWreDo6IigoCB4e3tj//796Nq1K0aNGiW1PKYE+HEGGdKsWTNcvHgRABAYGIiIiAhYW1uDiBAQEFCoR8jk07hxY9y8eRMA0KpVK0REROTt8/f3Z98MkHu+JScnY9euXdi0aRMiIiLQs2dPDBo0CKGhoVJLlCUFr1NPT088fPhQ7z5GfnCPT4YkJSXht99+w/bt25GZmQlra2sAgEKhgEKhkFidfOnfvz9GjBiBu3fvom/fvli0aBEePHiANWvWwNPTU2p5siX3nHJ2dsZbb72Fffv24caNG2jTpg2++eYbidXJF6VSiVu3biEiIgIajQbnzp0DINLsOTk5EqtjSoJ7fDJkxIgRhQLcN998A1dXVzx58gRDhgzB4cOHJVQnb9asWYPly5cjMjISmZmZ8PDwQJ8+fTB58mS4uLhILU+WtGvXDidOnJBaRqXj8OHDGD9+PJRKJX788UcsXLgQV65cQVJSElauXIk+ffpILZExAAe+SsKwYcPw888/Sy2j0vHWW29h/fr1UsuodPD5Vj569uyJ3bt3Q6nkZJqcUUktgClOr169im07cuQIEhMTAQC7d+82t6RKgSHfcrezb/op6hsR4ejRo3y+lYK+8+3YsWN5PT32Tb5w4JMhUVFR8PHxwZgxY6BQKEBEiIiIwMSJE6WWJmuio6PRtGlT9q2M6Dvfzp07x76VAl+nlRhiZEdOTg4tWLCAOnfuTBcvXiQiorp160qsSv6wb+WDfSsf7Fvlhcf4ZEx0dDQmTJgAV1dX7N69u1C5NGMY9q18sG/lg32rfHCqU8bUrl0b27Ztw969e+Hs7Cy1nEoD+1Y+2Lfywb5VPrjHxzAMw1gUXHNbCcjKAh49SoFOJ7WSSkZmJlIePQK4bVcmMjKAx49T2Laykp6OxOhoPt8qARz4ZEpkJPDxx4CLC2BvD3h62sPaGmjUCFi1CkhLk1qhTLl2DXj7bcDREVCrYe/hAahUgL8/sHEjkJkptUJZcuECMHSoONccHIDate2hUgGtWwPbtwPZ2VIrlCFEwF9/Af36Aba2gJMTHD09AWtrIDgY2LsX4BlcZAmnOmXGkyfAoEHA6dPimtF3w3FwENfcxInAzJkAPysL4N494M03gX/+EaZptcWPcXQEFAph2oQJ4r2F888/4ny7e1e0CfTdpx0dRdthwQJg5Ejza5QlZ84AQ4aIC1aj0d/Lc3QULYnly0VwZGQDBz4ZcecO0LYtkJCg/75dFAcHoEsXYOtWwMrK9Ppky6VLQIcOQHIyjMoHq9XiprVihUUHvz//BLp1E9kDY+4CarXIQnz1lem1yZq9e4EBA0TAMwZ7e2DOHGEeIwusZs6cOVNqEQzw7BnQpg0QG2vcvRsQHZsHD8RLzyQSlsGDB8DLLwOJicaPrWRnAzdvAunpQMeOptUnU65dA9q3B8qyBF92NnD+PODkJM5Vi+TMGaBnT+ODHiBasSdPAvXqAX5+ptPGGA0nyWTCf/4DxMWJ4YFVq4D790UH5uJFoGtXcYy1NbBtm8jqEYlhBI0G+OUXcT2ag6ioKHTo0AFNmzaFj48PFi9eDADYtm0bfHx8oFQq82apNwvvvSeqMQyZ1qYN8McfQHy8aFVs3Qq4uYluzrffim62GTDk26RJk+Dt7Q1/f3/07ds3b5owUzNypAhkhmxr0gSIiBDZh4QE4OBBsU2jASZPFlaaA0O+TZs2Df7+/ggMDERoaCgeP35sejFEwODBIpAZMq4g06aJ73TqJIwbO9Ysg/OGPMvl22+/hUKhyFu02SKR7tl5JpfkZCJ7eyKASK0mmjGDyMuLSKEg6tFD7PfyIrK2JvroI6K2bYkePyYKDhbfUSqJ+vc3j9bHjx/T+fPn/9WdTA0bNqSrV6/StWvX6MaNGxQcHEwRERHmEkNka1uyaV27CnOcnITJP/1EtH+/MM7amuiDD8wkVb9vv//+O2VnZxMRUVhYGIWFhZlcy7VrwoqSbHNxEf/mnl8ffEB0+bL4bG9PNHu2yWUSkWHfkpKS8o5ZvHgxjRs3zvRi/vyTyMGhZONEqCOqV4/oyhWiR4+IOnUS2xwdiX780eQyDXlGRPTw4UMKDQ0lT09Pevbsmcm1yBUOfDLghx/E9ZR7zRR9Xb5M1K9f4W1RUfmBDxD3fynO4169etEff/yR99msgW/GDCI7O+NNA4iaNRM3qdzPjo5EGo159BagqG9ERDt27KDBgweb/LfHjSNSqYy3zcqKaPx4orS0/G0vvkik1ZpcajH0+TZnzhx65513TP/jffqIQGeMcfv3E3XrRnTvXn7gA4jq1yfS6UyvtQAFPXv99dfp0qVL5OXlZdGBj1OdMmDdOsMZkBo1xCMMV6+W/N+wtgb27694bSVx//59XLx4EW2kGvDZsEGkOYtSkmnt2hXerlQCx46ZTKI+DPm2evVqdOvWzeS/v22b/uIpfbY9fy4s/u47UZ+RS2ameATCnBT1bcqUKfDw8MDGjRvxxRdfmPbHiYA9e/SPIxc1rn9/YZC+C/LxYzEubSYKerZr1y64u7sjICDAbL8vVzjwyQBD4yUqlXj0bN06UYtREllZYozQXKSmpuL111/HokWLpJumKSGh+LaSTPPzA6ZPByZNyt+m05nVOEO+ffXVV1CpVBgyZIjJNSQnF99myLaqVcWzpO+/L4ayclEqpT/fvvrqK0RFRWHIkCFYunSpaQUYapkWNc7RUbQQPvpI//HW1mYzrqBnKpUKc+bMMX0DobIgdZdTH/PnzycAFvS6WyxzolAQbdpEtHev/rRU0VQnkEnAxArVtXbtWr1/n6ysLAoNDaVvv/222D5Dqc6kpCRq3LhxhepLKItp9esTRUcTDR1aaHsyQMMr+O9ZNBVXmm9r1qyhoKAgSktLK/YdjUZD1apVq+DzLbtM51ruMXFxRNWr5257TkCPCtV19OjRMp9vREQPHjwgHx+fQtvS0tIqVJsTQJnGnG///S/RtGn5xxRJdT4HqHUFn2///PNPqZ5duXKFqlevTl5eXuTl5UVWVlbk4eFBMTExhb4XExNTodrq1q1LycnJev9uUiLLwGdpBAYWv9GsXk105IjhIayigU+tJlq1yvRadTodvfXWW/TRRx/p3W/WMT4PD+NM8/QUN6Bx44ob6exMtGuXyaUa8m3//v3UpEkTio2NNbmGXJycynau4d9xPo0m/1x1diYKDze9VkO+3bp1K+/9kiVL6PXXXze1EFHlU5pxFy+KwfaYGPHSaoni44nCwsR+BweimzdNLLXka5SILH6MjwOfDJg1q/C1s2wZ0alT+gtebGxEIUtUFFFIiHgPiO/fvWt6rSdPniQA5OfnRwEBARQQEEB79+6lHTt2kLu7O9nY2FCNGjUoNDTU9GLeeUdUZpZkWq1aRHfuEE2cqP+ObmcnujImxpBv9evXp9q1a+dtM0d1Yt+++TUahmzr3FkEOaVSBMrFi0WBYu75Zq6aIEO+9evXj3x8fMjPz4969uxJ0dHRphfz8sulX6TVqhG5uua/Hj4UVcW5x1WvbvKqIEOeFYQDHyM5MTH5gc/TU2xLTydKScl/DR6cnzkpipeX6P1ZHDduiNr6kkybPl3sK7g9JUWYqVQSDRggnX6JCA8X9+GSbOvfn+j6dfE5NpZozx4iPz9hm7U10YcfSvv/IAm7dolWQGkXacFXwVSnvT3RV19J+r/ACHjKMpnQuzfwv/9Bb9FYaTg6Alu2AN27V7wu2RMUVP6n99VqUdHZqlWFSpI7RECDBmJ+zvJgZyfm+Kxfv2J1yZ6cHDH5QXmLU2xtgagooHr1itXFlBmu6pQJ8+aJuTfLip0d0KyZ/okjLILvvhMBrKzY24uJTi0s6AFietIVK4QFZUWtBoYNs8CgB4gJcb//vvzGhYVx0JMJHPhkQuPGYu7bstzD7ezEDWjvXgteoaFVK2DTprLdjOztRWth0ybT6ZI5nTsDP/xQNtvUajEX+A8/mE6X7BkwAPjii7JdqGo18MYbwKxZptPFlAlLvV3KknbtxFy2rq5iImBDqFTihtWpk8jylXSsRdCrF3DggHjozNHR8HE2NqK10K8fcPSoSD1ZMCNGiBS5k1PJttnaCttGjQJ27bLwlUAA4NNPgWXLREArKU1jby+M+/RTYM0ai14JRG7wGJ8MyckRvbi5c8XsGDY2QGpqMhwdnaHVAm+9JZ6PbdJEaqUyIysL2LlTGHftGmBri+SUFDjnLmA4dizwwQdA3bpSK5UV6eliNpe5c8UE6DY2QHJyMhwcnKFSiXnA330XcHeXWqnMSE0VswfNmwc8fQpYWyM5OVmcb3Z2YhmisWPFzC6MrODAJ3NiYsQsR8+epcLd3RH165dvSMviiI4GnjxBanw8HGvXFjlhOzupVcmeBw/ETEKJiRq4u6vRsKGYbIQpASJRKRQXB01qKtTu7kDDhtw1ljEc+BiGYRiLgsf4GIZhGIuCAx/DMAxjUXDgYxiGYSwKDnwMwzCMRcGBr5Lx9ttvSy1BtuTk5GDFihWYNm0awsPDC+2bPXu2RKrkj0ajwbx58zB//nxkZGRg7dq16NWrF8LCwpCamiq1vEpFo0aNpJbAGAFXdcqQBH0LrAIgIgQEBCA6OtrMiioHY8aMgUajQevWrbF+/XoEBwdjwYIFAIDmzZvjgrmXDK8kDBgwAB4eHkhPT8fNmzfRpEkTvPnmm9i9ezeePHmC9evXSy1Rljg5OUHx70PpubdRjUYDtVoNhUKBZH0r/jKygAOfDLGysoKXlxcK/mkUCgWICI8ePUJWVpaE6uSLv78/rly5AgDQarUYP3484uLisGnTJgQFBeFiwSXEmTwCAwNx6dIlEBFq1qyJmJiYvPMtICAgz1OmMB9++CESExMxf/58uLq6AgDq1q2Le/fuSayMKQ2V1AKY4tSrVw+HDx+Gp6dnsX0eHh4SKKocFGwQqFQqrFy5ErNmzULHjh05ZWcECoUC3bt3z+vFKBSKvPdMcZYsWYLz589j0KBB6NOnD95//332q5LAY3wy5OOPP8bz58/17gsLCzOzmspDy5YtceDAgULbZsyYgZEjR+L+/fvSiKoEtGzZMq9hsHr16rztkZGRcLL4iWBLpkWLFjh06BAAIDg4GBkZGRIrYoyBU52VhGHDhuHnn3+WWkalg30rH7m+ERH3YowkJiYGvr6+iI+Pl1oKUwqc6pQhvXr1KvSZiHD06FEkJiYCAHbv3i2FLNnDvpUP9q18FPUNADIzM/O2s2/yhQOfDImKioKPjw/GjBmTV2Rw7tw5TJw4UWppsoZ9Kx/sW/mIjo5G06ZNC/kWERHBvlUGiJEdOTk5tGDBAurcuTNdvHiRiIjq1q0rsSr5w76VD/atfLBvlRce45Mx0dHRmDBhAlxdXbF79248fPhQakmVAvatfLBv5YN9q3xwqlPG1K5dG9u2bcPevXvh7OwstZxKA/tWPti38sG+VT64x8cwDMNYFPwcXyUhMzNTagmVEvatfPAD/+WDfasccOCTK3FxwNy5gJcXoFZDZW8PODkBbdsCu3YBOTlSK5QlMTHA9OmAuztgbw/Y26vg4gJ07gz88Qeg00mtUKY8eABMmgS4uQF2drBzdgaqVAFeew04cQLgxJBebt0C3n8fqFEDsLUFXFzsUK0a8OabQESE1OoYQ3CqU24kJgLvvCOCm0IBpKcXP8bJCbCxAb78Enj3XfNrlCFPnwJjxgAHDwrb9E2g4egorJs/HxgyxPwaZcn9+8Do0cBff4lWQdF5YBUKQK0GqlcHvvsO6NlTEply4/p1YNQo4PJlQKsFsrML71cqATs7wMMDWL4caN9eEpmMATjwyYnHj0WP7vHj4jcgfajVwPDhwPffixuUhXLnDvDKK0BCQvEbkD7UauDTT4FZs0yvTdZcvgx06AAkJxuXQbC3F1mIDz4wvTYZEx4OdOsGpKYa1xG2twd+/JEbW3LCaubMmTOlFsFA3Hxeegl4+NC4uzcgjrt2DdBogE6dTKtPpsTGAq1bi3+Nzf5mZ4s0lKMjEBRkWn2y5f594OWXgefPjU9jarXA8eNA3bqAn59J5cmVq1dF760sQ3laLbB/P9CiBdCggcmkMWWAx/jkwvTpQFQUMG6cuCtnZABr1uTvHzwYSEnJf6WliRuWtzewcKFovVsgH3wghkPHj9dvGwC88YZoHyQnixtX796irfDZZ8Jyi2TkSCApCXjvPcPGjR4N3L4tzrf9+4GaNYVxY8aIlLwFMnCgaDitWiXaDsnJwMWLQNeu+cd07ChSoWlpwJEjgKenGLF4802Aa63kAQc+OZCeDvz0k7gqHj8GZs8GCsySDwD45RcxQJX7Gj8eiIwELlwQadGFC80iddSoUahRowZ8fX3ztr355psIDAxEYGAg6tSpg8DAQLNoiY8Hdu8WLWpDttWqBWzYAHzyCeDsLOo3fvlFDFkRiSyxOdDn26VLlxAUFITAwEC0bNkSZ8+eNY+Ye/eA06fFmJ4h44KDgTlzRCuhWjXxnU2bxD6FAli71ixS9fl2+fJlvPTSS/Dz88Nrr71mtgVfz58H7t4FrKxEgyk4GHBxAaZOBbZuFXVoL7wA7NgBTJsmbDt3DtiyRXxfpwN+/dX0OqOiotChQwc0bdoUPj4+WLx4MQCxwHVISAgaNmyIkJAQgyvAWATSTRrD5LF2LZGjI5G4F4vXl18SrVlTeFvB15EjRNOn53+2syNKTDS51OPHj9P58+fJx8dH7/5PPvmEZs2aZXIdRERz5xLZ25dsW+vWRE+fFj4mNpYoKEi8d3Ehysw0vVZ9voWEhNC+ffuIiGjv3r0UHBxseiFERJ98QmRjU7Jx8+cTLV2a/7lmTfHdevXE51q1iHQ6k0vV51vLli3p2LFjRET0008/0dSpU02ug4ho8GAiKyv9l+Ply0T9+hGNHUsUHp6/Xa0m0miIGjcWn/38TK/z8ePHdP78eSIiSk5OpoYNG9LVq1dp0qRJ9PXXXxMR0ddff01hYWGmFyNTuMcnB378sWyDBp6eQLt2QMHldqytgX37Kl5bEdq1a4dq1arp3UdE2Lp1KwYNGmRyHYDopOgrei3IuXMi7fTaa6LSrndv0bHOXVScSAxbmRp9vikUirzeSlJSEmrVqmV6IQCwfr1xxVMFC6Zy3+f2vJKTgUuXKl5bEfT5duvWLbRr1w4AEBISgu3bt5tcBxGwfbv+ceQaNYBGjUQa3cen8KiDRiMSMz4+4vPt20B0tGm11qxZE82bNwcAODk5oUmTJnj06BF27dqF4cOHAwCGDx+OnTt3mlaIjOEpy+TAkydlO37YMODkSTHIkEtmpqjwkJCTJ0/C1dUVDRs2NMvvxcWVfoxOJ9oHv/wiysuzssSYn0aTv18q2xYtWoQuXbrg008/hU6nw19//WWeHzZmfO7AAWDzZlGLf/u2GIPW6URJLCDyfRIZ5+Pjg127dqFPnz7Ytm0boswwUJuerj/oqVTAxo3AunXAzZuiYOrZs8LHJCWJ0QlAPIUUGwvUrm1yyQCA+/fv4+LFi2jTpg2ePn2KmjVrAgDc3Nzw9OlT84iQIbLs8X377bdQKBQW87ofGVk2g4YNE1daAbKysjDp448rVNe6Ir9RGps2bTLY20tOToa3t3eF6ouPL/0G3qkTMG+eqMSzsRHjMqtWAQEBYn9qagqGDh1Zobo2b95slF/Lli3DwoULERUVhYULF2L06NHFjtHpdLCysqpQfTnGVA0fPgzMmCG6Offvi1dKSl53JTEpCa917Vqhui4bWaC1evVq/PDDD2jRogVSUlJgY2Oj97iK1Obg4AKtNrPIfz+/8/z++2JbaqoYSy6Is7OwDgCSkxPRokVQhWq7evWq3v//1NRUvP7661i0aFGxOURzv1uUJ0+eVKi2+vXrIyX3f15OSJ1rZYjI37/4oIGhMb6XXyZKTS0+JqhWE/34o1nk3rt3r9gYX3Z2NtWoUYOioqLMooGIyN29dNsmTiTasaPwMb/9JrYDRM7O4rM5KOqbs7Mz6f4dJ9PpdOTk5GQeIUXPHWPGlBs2FOddlSr5g6MnT5pFrr7zLZebN29Sq1atTK4hJ4dIqSxsyerVYqjdzi5/29ixRH/+WfiyTEvLH+NzdCS6ccPkcikrK4tCQ0Pp22+/zdvWqFEjevz4MRGJccBGjRqZXohMkWWPz+Lo2lV0RwCRQrK1Ff8WfJ/L8OGiFV50TJBIPMUtEYcOHYK3tzdqmyuHA9Gby7XGkG0REcCrr+b38AIDxefcMb7MTOme5atVqxaO/zvAeOTIEbOliNG2bf57Q8bZ2uYPTHl4ACtXAosX56dJs7LyTTUzsf+mWHU6HWbPno133nnH5L+pVBb+3122DGjSRIwdF5wl6LffxDBov37CwunTxbl286bYr1IB9eqZVisRYfTo0WjSpAk++eSTvO29evXKy+KsW7cOvXv3Nq0QOSN15GWI6MGD/GbjjBnF98+YIfbZ2hI9f07UsWPxFrkZWr1ERAMHDiQ3NzdSqVTk7u5Oq1atIiKi4cOH07Jly8yiIZeLF0WLujTb3nuP6PZtouRkoshIUdQIECkURD17mkerPt9OnjxJzZs3J39/f2rdujWdO3fOPGIOHcrv9RkyzsVFlCqmphLFxBDNmZPf5VGpiEaPNotUfb4tWrSIGjZsSA0bNqTJkyfn9ZpNzebNRE5ORJ6e4nN6OlFKSv5r8GBhT6dORNevi2rOo0eJvLzyL98pU0yv8+TJkwSA/Pz8KCAggAICAmjv3r0UFxdHHTt2pAYNGlCnTp0oPj7e9GJkCk9ZJhc6dxbjKuXB0VE8V/X66xUqqTLg7w/8/Xf5vuvoCOzdKwpkLQoi0Yt79Kh837e3F13p3B6hhZCVJSo4k5LK931bW1Hh6e5esbqYssOpTrkwf35+xVxZsLERuZNevSpeUyVgyRJxHy4rdnZiCqlXX614TbJHoSi/cfb2YqJqCwt6gLjU5s0r32WqVosJbzjoyQOeq1Mu1Kwp5j/ctUtMRWIMNjZiGZnw8Px6aQujTh0xO8sffxhvm52daCscPSpa4RZJkybi/AkPN35uWHt7MUi6a5cYrLJAWrQQw5yXLhlvm1otqop//lmMFTLSw38GOdGrF7Bnjwhijo6Gj1MoAAcHMdp+8SLw4ovm0yhDRo4Uz+mp1cIWQyiVYv8rrwBnzpRssUUwaZJYasjOruTen5WV2N+jh4W3FgTz54uVPezsxMsQ1tZi/9ChYmq9gjVqjLTwGJ8cSUsT8yLOnStWVlWpkJycLJ7FycwEQkLETeuVVyx6OaKiJCWJVvX8+WLRASsrMSOKi4sLMjOBPn3EnJ2tWkmtVGbEx4u5YhcsEE/2K5X551tWlpiZecIEi12RwRAxMcCKFSJrrNWKSzE5WZxvWi0wYgTw4YdiVhdGXnDgkzNEYh6k6GhkPH8OO1dXMbbi6iq1MllDJFJRT54AaWlavPiiCv7+YtJgpgR0OjHpeWwsMpKTYVezJtCsWfEnsplCaLViary4OCArS4vq1VVo0aJ8Y4GMeeDAxzAMw1gUPMbHMAzDWBQc+BiGYRiLggMfwzAMY1Fw4GMYhmEsCg58MqRfv37YsGEDUsuyOC2Du3fvYtSoUZg6dSpSU1MxduxY+Pr64o033sD9gmsXMoXQ6XRYvXo1evTogYCAADRv3hwDBw7EsWPHpJYma7RaLVasWIGuXbvC398f/v7+6NatG5YvX45sY59uZySBqzpliLu7O1566SUcOXIEnTt3xqBBg9CjRw+D644xgnbt2mHQoEFISkrChg0bMHLkSAwYMAB//PEHNm7ciCNHjkgtUZaMHDkSXl5e6Ny5M3799Vc4Ozvj1Vdfxdy5c9G7d2988MEHUkuUJYMGDUKVKlUwfPjwvFVJoqOjsW7dOiQkJGDLli0SK2QMwYFPhjRr1gwXL15EcnIydu3ahU2bNiEiIgI9e/bEoEGDEBoaKrVEWZLrGwB4enri4cOHevcxhfH398eV3HWaAAQFBeH06dPIzMxEYGAgrl+/LqE6+dKoUSPcunWrzPsY6eFUpwzJXRnZ2dkZb731Fvbt24cbN26gTZs2+OabbyRWJ1+USiVu3bqFs2fPQqPR4Ny5cwCAO3fuICcnR2J18sXa2roQWjkAACAASURBVBqRkZEAgAsXLuRlFmxtbfWu0s0IqlWrhm3btkGn0+Vt0+l02LJlC6pWrSqhMqY0LHOmWZnjqGcSyRdeeAHvvPOOWRbdrKzMmzcPr732GpRKJXbu3Imvv/4aV65cQVJSEn788Uep5cmW+fPno0OHDrC1tYVWq8XmzZsBAM+ePUPPnj0lVidfNm/ejMmTJ2P8+PGoWrUqiAhJSUno0KFDnoeMPOFUp0w5e/YsFAoFWrVqhWvXruHAgQPw9vZG9+7dpZYma86cOQOlUolWrVrh6tWr2L9/P5o2bcq+lcKpU6egUqn4fCsn8fHxAICPPvoIGzZskFgNUxoc+GTIrFmzsH//fmi1WoSEhODMmTPo0KEDDh48iC5dumDKlClSS5QlRX07e/Ys2rdvz76VAvtWPnrpWQPzyJEj6NixIwBg9+7d5pbEGAkHPhni5+eHS5cuITMzE25uboiOjoazszPS09PRpk2bQoUITD7sW/lg38pH8+bN0bRpU4wZMwYKhQJEhEGDBuWlOYODgyVWyBiCi1tkiEqlgpWVFdRqNerXry+WhwFgb28PJa9kaRD2rXywb+Xj3LlzaNGiBb766iu4uLigffv2sLe3R3BwMAc9mcPFLTLExsYGGo0GarUa58+fz9uelJTEN6ISYN/KB/tWPpRKJSZMmIA33ngDEyZMgKurK7RardSyGCPgVKcMyczMhK2eVa7j4uIQExMDP14QVC/sW/lg3yqGvXv3Ijw8HHPmzJFaClMKHPgYhmEYi4LzGMz/a3jORMac8PlWOeDAJ1PS04G1a4GWLYGaNQG1Oh2enkD//sCpUwD30w2QkgIsWwYEBABubsh2dgbq1AHeegvgKcsM8vw5sGAB4OMDuLqK861uXWDcOIBnLCuB2Fhg9mygcWOgRg1kOTsD9eoBH38M3L0rtTrGAJzqlBnp6cDnnwM//ggoFEDRBRoUCkCtBtzcgLlzgddfl0an7EhOBiZOBDZuBJRKIC2t8H4rK8DWFqhbV9zheb5TAMCzZ+IevWOHsE2jKbxfpQKsrUVAXLQIaNtWGp2yIyoK+PBDYP9+cVFmZBTeb20tzrlWrYAlS4DAQGl0MnrhwCcjEhKAjh2BmzeLX0f6UKuBsDBgxgzTa5M1jx8D7dqJm1FWVunH29sD8+cD771nem0yJjISePVVIC4OMCZDZ28P/PQTMGiQ6bXJmitXgA4dgKQkwJg5YNVq4NdfgW7dTK+NMQoOfDIhIwN4+WXg6lXj7t25qNUi0zJhgum0yZqkJKB5c+DhQ6AspeRqNbBiBTB0qOm0yZinT0Un5OnTsqXN7e2B7dst+B5+/744354/L9v31Grg4EFxkTOSw2N8MuHLL8VYytixQESECIRr1uTvb9JEbE9IEK+DB8U2jQaYMkX0Ei2STz4BHj0Sg1H6jCvItGniLt+pkzDu7bfFnd8CGT0aiI8Hxo/Xb5uXl7AqJSX/NXWqSMW/8UbxFLzF8OaborH13nuGzzd7e+D770UeOTEROH5cnG+9ehnXtWZMDgc+GZCVJa6TjAyRtZs9G1i9uvAxjx+LwpZq1YAXXwR27wZyJ4DXaoHFi82jddSoUahRowZ8fX3ztk2bNg3+/v4IDAxEaGgoHj9+bB4xycnApk1AZqZh43KpV0/csYtqW7nS9Dqh37dcvv32WygUCsTFxZlFy+PHwKFD4h5cmm1VqgBOTuI1e3b+9o0bzSJVr28zZ86Eu7s7AgMDERgYiH379plHzLVrwN9/AzpdycatXCku1CZNxL+56ZisLOB//zOLVEPn23fffQdvb2/4+PggLCzMLFpkCTGSs2ULkZMTkWhji9eXXxKtWVN4W+7Lyopo/HiitLT8bWo1UWqq6bUeP36czp8/Tz4+PnnbkpKS8t4vXryYxo0bZ3ohRETffUfk4GCccfv3E3XrRnTvHlGnTvnbX3iBSKs1uVR9vhERPXz4kEJDQ8nT05OePXtmch1ERFOnEtnalmybl5c41spK/zlYrx6RTmd6rfp8mzFjBs2fP9/0P16UMWOIVKqSjWvcmCgpqfgFnftq3dosUvX5duTIEerUqRNlZGQQEdHTp0/NokWOcI9PBixbJlJJxvD8uegZfvcdUHCCCCsrwBwN33bt2qFatWqFtuXO7QgAaWlp5lu8dPny4tWb+ujfX/QK9+8vvi8rCzh5suK1FUGfbwAwYcIEzJs3z6wLvq5aJewwhgcPRM3Q6tXACy/kb3/6VIxHmxpDvknCxo2ljyO3bi1MmzVLpDqvXAH69cvff/kyEBNjWp3Q79uyZcvw2Wef5c3SU6NGDZPrkCsc+GRAdLTxx1atCri4AO+/X/ixtKwss1xPBpkyZQo8PDywceNGfPHFF+b50SdPSj/G0VG0ED76yPAx5krNFmHXrl1wd3dHQECAWX/336XjSiQuTjxD6uUFtGghUp0F05sqlWS2AQCWLl0Kf39/jBo1Cs/LWmhSHtLTjas6q10b8PMT44C1aokLdd06wNtb7Le1Ne68NQG3bt3CyZMn0aZNGwQHByMiIkISHXJAloFv4cKFUCgUFvO6c+dBmfzRaERn5+efgerVxbbMzGx89FFYheqaP3++0Zq++uorREVFYciQIVi6dKkezRq4ublVqL4kY+7gM2cC69eLVrgeUlNSMGrIkArVZczq2xqNBnPmzCm1kaDT6WBlZVWh+rKzSy/BT0sDzp8X1fqxseL+3aWLaEcAQFJSIrp06VWhui5fvlyqLgB49913ERkZiUuXLqFmzZqYOHGi3uMqUlsVtRqZxjy6kBsgZ88Wg6gnTgBHj+Y9N5qYnIyXmjevUG1RUVFG+abVapGQkIDTp09j/vz5GDBgAKhISW9KSkqFamvQoAFSjE1nmRFZBr4JEyaAiCzm5evrVWaPlEpRIe3uLj6r1dZYsWJeheqaNGlSmXUNGTIE27dvL7ZdrVbjyZMnFarPpVat0gV16iQeNI6JES8PD2DrVvEAJABHZ2es3r69QnUNHDiwVFmRkZG4d+8eAgICUKdOHURHR6N58+Z4UqQ3oFQqkZOTU6H61Gor4/6YBSDK1SP+dXGpguPHd1eoLmN7vq6urrCysoJSqcTYsWNx9uxZA5orTluiVgtbhRHpaH1rF1J+cKni6IhT165VqDYPDw+jfKtduzb69esHhUKB1q1bQ6lUFiuocnJyqlBtd+7cgZOTk1H6zIksA5+l0b69mOgByJ9gxMqq8PvOncVzV0qlSDstWCDG+wpOJ9W6tSTycfv27bz3u3btgnduWsfUvPxy/p3YkHGdOgG+vsK8wECRnxs3TpTRAqJ13rKlefQWwM/PD7Gxsbh//z7u37+P2rVr48KFC3BzczP5b7dokf/ekG2tWwONGolJSapVE5OPHD0qCmkBMUYo1aINMQVy+r/99pveStkKx8oqP12Z+1mfcSdOiGdK//Mf8fnll8XD7r//nv/devVMr1cPffr0wdGjRwGItGdWVhZefPFFSbRIDjGSc/s2kZ0dEUA0Y0bx/TNmEPXvT3T9OlFKClFsLNGePUR+fpRXLObnZx6tAwcOJDc3N1KpVOTu7k6rVq2ifv36kY+PD/n5+VHPnj0pOjraPGJOn86v6jRkXNGquqJVnZ06mUWqPt8K4uXlZbaqzj178osODdk2cCDR3buiUvjxY6J164hcXcV3lEqiwYPNIlWvb0OHDiVfX1/y8/Oj1157jR4/fmweMWvXEjk6ln6+NW1K9NdfwryrV4n69BHbbWyIJkwwi1R9vmVmZtKQIUPIx8eHmjVrRocPHzaLFjnCM7fIhLZtgb/+Kt93HR3FmN+QIRWrSfYQicmBC/Q4y4Sjo5iGxMLm7czJEROfP3tWvu+r1aIQtnnzitUle9LTgRo1yv/0vp2deBawbt2K1cWUGU51yoS5c8WED2XFykpci/37V7wm2aNQAN9+Wz7jbGyAhg1FDtnCsLIS55taXfbv2tmJ7J3FBT1AnGdTppTPOHt7oG9fDnoywWrmzJkzpRbBAJ6eolDljz+Mn3LSyko83nDqlBiHsUgaNxZBLDzc+OmgbGzE8hbh4YCDg2n1yZRmzcRsWpcuGW+bra0Ynjp8WFhokbRtKzIMt28bb5y9vVgma9cucdEyksOBT0YEBor7+N694vooKQA6OooCxVOnxL8WTdu2YhG5Q4eEcSWVnTs6CpP//FPM/WbBhIQIu8LDRedZpzN8rKOjqAE6elQUV1ksCoWYczMpKf9BWkPGKRSiYRUSAuzZI1oOjCzgMT4ZEhsrpvtbvFhUzxEBqalpcHJygFYrqukmTwZeey2/GpSBmAnghx/EiwjQ6ZCWmgoHJyfROg8KAiZNEg+kccs7j7t3xUxAq1aJItmcHDEDj5OTA7KzRVHipEmi+tiYin6L4do1cZFu2CCe6M/JQWpaGhydnUW1cPfuwKefivOOjZMVHPhkTE6OmNg9OhpITs7Giy9ao1kz0WFhSiA7W3RNYmKQnZoK6+rVRX1+nTpSK5M1GRnAkSOi4ZWRocMLLyjx8sv5z4oyBkhNFcbFxUGn1UL5wgtifcjc2SUY2cGBj2EYhrEouKqTYRiGsSg48DEMwzAWBQc+hmEYxqLgwMcwDMNYFBz4ZIhGo8G8efMwf/58ZGRkYO3atejVqxfCwsKQWt7pkiyURo0aSS2BYRiZwYFPhowYMQJPnz7FvXv30KNHD5w7dw6TJk0CEeHdd9+VWp5scXJygrOzM5ydneHk5AQnJydERkbmbWf0c6XAUjrZ2dmYPXs2evXqhc8//xwajUZCZfJm6dKlecv63LlzB+3atUOVKlXQpk0b/P333xKrY0qCH2eQIYGBgbh06RKICDVr1kRMTAwUCkXemmVX9K35xeDDDz9EYmIi5s+fD1dXVwBA3bp1ce/ePYmVyZvmzZvjwoULAICJEyciPj4eI0eOxM6dOxEfH4+ff/5ZYoXyxMfHB1evXgUA9OjRA2PGjEHfvn1x7NgxTJkyBeHh4RIrZAyhkloAYxiFQoHu3btD8e+sD7mrGjP6WbJkCc6fP49BgwahT58+eP/999kvIyjY9j18+DAiIiJgbW2Ndu3aGb04rCWiLTCnYGxsLPr27QsAaN++vSxXHWfy4VSnDGnZsmXeWN7q1avztuem7RjDtGjRAocOHQIABAcHIyMjQ2JF8icpKQk7duzA9u3bkZmZCet/58HjhlbJ9O/fHyNGjMDdu3fRt29fLFq0CA8ePMCaNWvg6ekptTymBDjVKVPOnj0LhUKBVq1a4dq1azhw4AAaN25cqAfIFKegbydPnsTRo0fRsmVLdO/eXWppsmXkyJGFPn/zzTdwdXXFkydPMGTIEBw+fFgiZfJn7dq1WLZsGSIjI5GZmQkPDw/06dMHkydPhouLi9TyGANw4JMhs2bNwv79+6HVahESEoIzZ86gQ4cOOHjwILp06YIpU6ZILVGWFPXt7NmzaN++PftmBGfOnIFSqSzU0PL29uYGQykUbGhdvXoVBw4cQJMmTdg3mcOBT4b4+fnh0qVLyMzMhJubG6Kjo+Hs7Iz09HS0adOGi1sMwL6VD24wlA/2rfLCgU+GNGvWDBf/Xeur4Hsgv+KTKQ77Vj64wVA+2LfKCxe3yBAbG5u856fOnz+ftz0pKQlKJf/JDMG+lQ+VSgUrKyuo1WrUr18/75lHe3t79q0E2LfKC/91ZMiJEyegVqsBoNAFlJ2djXXr1kklS/awb+WDGwzlg32rvHCqk2EsnMzMTNja2hbbHhcXh5iYGPj5+UmgSv6wb5UXDnwMwzCMRcH9cbmSkwPs2wf06AH4+iLpxReB5s2BDz4Abt+WWp1syc4Gtm8HQkMBHx+gevVEtGwJhIUBDx9KrU7GZGYCGzcCHToATZsixc0NaNUKmDEDiImRWp1s0WiAn34C2rUDvL2BmjWTERQEzJ0LxMdLrY4xBPf45IZWCyxcCPz3v+KqKroag7U1YGUFNGsGfP01EBwsjU6ZkZkJzJkDfPedsLDojFE2NoBSCbz8MvDNN+KezkAY9cUXwIoV+Z8LYmcHEAEhIcI4Hx/za5Qh8fHAtGnAunXivCp6mdrbC9t69RKXab160uhk9MOBT06kpYkr5fRpEfRKw94eWLAAeOcd02uTMc+fA507A9evA+nppR+vVgNr1gADBphem6yJiQHatwcePBAth5JQKIRxO3aI7rQFc/eu6OHFxooMQ0kolYCjI3DgAPDSS+bRx5QOBz65oNWKG8qpU0BZ5pdUq4Fly4Bhw0ynTcakp4sbyvXrQFaW8d+ztwd+/RWw2Ak2kpJE6vzhQ3HuGYtaDRw8KLrOFsjTp0BAAPDsGaDTGf89R0fgr78ArneRBzzGJxfmzwfOnCke9Ly9gcOHgcREMbbXp0/h/RqN6PHdv282qXLis8+AmzeBsWOBiAhh35o1hY/p2FEExrQ04MgRwNNTBMwBA4CEBGl0S864cUB0tPhXn3HW1sC2bcC9eyJnl5tS12jEuLOFTv49aJBIc777rn7b2rQB/vhDHBMbC2zdCri5iVRot25lC5aM6eDAJwdyckTKsmh608oK2LUL2LMHqFYNePttYMMGoGHD4t///nvz6S3A4sWL4evrCx8fHyxatMisv51bWJCRATx+DMyeDRRYzAIA8MILIjs3bZqw8Nw5YMsWsY+oeJA0FaNGjUKNGjXg6+ubty0hIQEhISFo2LAhQkJC8Pz5c/OIiYsDdu4UXWRDxgHAn38CQ4cWL27RakVQNAP6fNu2bRt8fHygVCpx7tw5s+gAgMhIkZDRag3bVrUqsHIlUKcO4OUlhkxzz7GkJOD3382jVZ9vkyZNgre3N/z9/dG3b18kJiaaR4wcIUZ6du8mcnIiEvfi/JePD1FKSuFtv/9O9MUXxY91ciLKyDCr7L///pt8fHwoLS2NsrOzqVOnTnT79m2z/f5PPxE5Oha24csvidasyf88dixReHj+Z7WaSKMhatxYfK5Zkygnx/Rajx8/TufPnycfH5+8bZMmTaKvv/6aiIi+/vprCgsLM70Q8WNE9vYlG1fwFRVFFBxc/Nw0A/p8u3btGt24cYOCg4MpIiLCLDqIiD76iMja2njbAKJmzYiSk/M/t29vHq36fPv9998pOzubiIjCwsLMd77JEO7xyYGlS4tX0xlCoQAKtOIKsX9/xWkyguvXr6NNmzZQq9VQqVQIDg7Gjh07zPb7331XvJquKD4+wOXL+Z81GtFyzy1OTE0VrXhT065dO1SrVq3Qtl27dmH48OEAgOHDh2Pnzp2mFwKI7IAxVUAlce+eyDGbGH2+NWnSBI0bNzb5bxdl9erSi1mK0q4d8O8i7QDEufbsWcXq0v+7xX0LDQ2FSiXWHg8KCkJ0dLTphcgUDnxywND43M2bYqBg0iRApRIl5cHBosCgKFlZYszGjPj6+uLkyZOIj4+HRqPBvn37EBUVZbbff/So9GMcHUWKqSBJSUDB9XzNKLkQT58+Rc2aNQEAbm5uePr0qXl+uCLuvNbW0hknAZmZYoy4LPj5AdOni8s3F1tbkSaVmtWrV6Nbt25Sy5AMWQa+hQsX5q3+bAmvh7du6TdCqxXFLD16AE+eABMnitFyPQEuOzMTYR98YDKN77zzTt68hLk0adIEkydPRmhoKLp27YrAwEBYWVnla8rONqlvz54lFbWhGKmpwL9zB+fh7JzfwU5JScWgQaMqVNfmzZtL1VWU3O8WRafTwcrKqkL16Up7dMEIEpOS0CskpEJ1XS7YNa8AKlKbnV1V6HTG+1a/vkjAfPSRGCrNJTk5EYGBL1WotrI2Nr/66iuoVCoMGTKk2L6UlJQK1dawYUOkGJvNMiMqqQXoY8KECZgwYYLUMsyHr2/hfEhB/v5bPGuVS3i4eGq2CNb29pi3cCHmjRtnGo0GGD16NEaPHg0A+Pzzz1G7du18TdbWIBM+LVOrVumTily9CvybTQQgOsv16+fb7ezsiNWrV+P11/UUd5gYV1dXxMTEoGbNmoiJiUGNGjWKHaNUKpGTk1OxP+zgYNxzoiVQxcUFu3ftkvUEChV57uXkiE6uMf9JT0/g0CHgyy9FLVpBnJyq4PTpU2jatMKklYm1a9diz549OHz4sN6GlpOTk0mvWbkgyx6fxREUJFKZ+vDzE/kRe3vR46tZE1i7tvhxSqV4wMjMxMbGAgAePnyIHTt2YPDgwWb77ebNxZAnIApgbW3FvwXf//abaFf06ye2TZ8OXLmSPzyVnS2JbQCAXr165a0asW7dOvTu3ds8P1xw9hVDxgFiupvcSZgLvgdE7k+qu7cEWFmJSs2Cn/XZVquWeGRm6dL8yXAKotUW/u+YkwMHDmDevHnYvXt33iomFovExTUMEdHffxevsst9zZtHlJAgqjv37SOqX1//cfXrE+l0Zpf+yiuvUJMmTcjf358OHTpk1t8+ciS/qnPGjOL7Z8wQ+zp1Irp+XVRzHj1K5OWVb1tQkHm0Dhw4kNzc3EilUpG7uzutWrWK4uLiqGPHjtSgQQPq1KkTxcfHm0fM1q35VcQlGXfvXvF9Xl5ECgVRr16m1fgv+nzbsWMHubu7k42NDdWoUYNCQ0PNouX774kcHEq2bfp08T4lpfALIFKpiN5+2yxS9fpWv359ql27NgUEBFBAQACNGzfOPGJkCM/cIheaNwcKrBheJhwcxPyeY8dWrCaZQySelSpvjYWTE7B+PWCujpZsyM4GatQQkyKUBwcH8UBa27YVq0vmpKQArq7lL4i1twcuXBBzUjDSwqlOufD11+LKKCsKhbiDmzHFKBcUCjExdXmyNlZW4ibWo0fF65I91tbiif7yGGdjI1KcFjhlmZOTWBylPLbZ2YnhUA568sBq5syZM6UWwQBo0EAMFISHG/+wUG7QCw8XY38WiL+/6Lhcvmy8bVZWYoaNU6eAKlVMq0+2BAUBt26JafCMNc7aWsy/dfJk+e7+/w/o1EmcN48eGW+bra1YneHgQdFuYKSHA5+cePllcUc+elQEtZIm9nNwEMeGh1t8MzIkRFh1+rT4XJJtjo6iAOHUKcDDwzz6ZIlCIXK8z56JyuGcnJJLFh0dReMsPByoXt18OmWGQiHmeL11S0yEUNr83o6OYgWxo0eLP1bDSAeP8cmRW7eAJUtE9aZSCWRkICc7G1a2tqL6s0YN8VTs0KGFn8S2cK5cEUOdW7aIXl1mJpCdnQM7OysolWI8MCwMePPN8mWV/98SESHWf9y9W5xfGRnI0WphlWuStzcweTLQty93Wf6FCDhxQth26FD++abVaqFWq6DTiWH7sDCgZ8/8QllGHnDgkzMaDfC//wHR0dBpNFBWqSKaj23b5tfxM8VIThb38CdPgPj4VNSu7Yg2bYCWLaVWJnPi48WE6LGxyEpNhY2rq5hzy9AUeQwAcZ7t3Svm/s7IyIKrqw06dSo+lzwjHzjwMQzDMBYFV3UyDMMwFgUHPoZhGMai4MDHMAzDWBQc+BiGYRiLggMfwzAMY1Fw4JMhqampmD59Onx8fODi4oLq1asjKCgIa/WtysDkodVqsWLFCnTt2hX+/v7w9/dHt27dsHz5cmSXdelsBgDw9ttvSy1BtuTk5GDFihWYNm0awsPDC+2bPXu2RKoYY+DHGWRI79690bdvX3Tu3Blbt25FWloaBg4ciNmzZ8Pd3R1z5syRWqIsGTRoEKpUqYLhw4fnrQsYHR2NdevWISEhAVu2bJFYoTxJSEjQu52IEBAQgGg9Cx8zwJgxY6DRaNC6dWusX78ewcHBWLBgAQCgefPmuHDhgsQKGUNw4JMhAQEBhVajbtWqFSIiIqDT6dC0aVPcuHFDQnXypVGjRrhlYDX7kvZZOlZWVvDy8iq0AKlCoQAR4dGjR8jKypJQnXzx9/fHlStXAIhsw/jx4xEXF4dNmzYhKCgIF8u72gpjcjjVKUMcHBzw559/AgB2796NatWqARCrcXM7xTDVqlXDtm3boCswWadOp8OWLVtQtWpVCZXJm3r16uHYsWO4d+9e3uvu3bu4d+8eXF1dpZYnWwo2CFQqFVauXImAgAB07NgRqampEipjSoMDnwxZtmwZPvnkE1StWhXz5s3DkiVLAADPnj3De++9J7E6+bJ582b8+uuvcHNzQ6NGjdCoUSO4ublhx44d2Lx5s9TyZMvHH3+M58+f690XFhZmZjWVh5YtW+LAgQOFts2YMQMjR47E/fv3pRHFGAWnOmXK9evX8ejRIwQFBcHR0TFv+4EDB9C1a1cJlcmbM2fOQKFQoH79+rhx4wZOnTqFpk2bonv37lJLkzVnz56FQqFAq1atcO3aNRw4cADe3t7sWymwb5UTDnwyZMmSJfjhhx/g7e2NS5cuYfHixej97zLhPGhumFmzZmH//v3QarUICQnB2bNn0b59exw8eBBdunTBlClTpJYoS4r6dubMGXTo0IF9KwX2rRJDjOzw9fWllJQUIiK6d+8etWjRghYtWkRERIGBgVJKkzW+vr6k1WopLS2NnJycKCkpiYiINBoN+fn5SaxOvrBv5YN9q7yopA68THF0Ol1eerNOnTo4duwY+vfvjwcPHnBxSwmoVCpYWVlBrVajfv36cP535U97e3solTycbQj2rXywb5UX/uvIEFdXV1y6dCnvs6OjI/bs2YO4uDj8/fffEiqTNzY2NtBoNACA8+fP521PSkriG1EJsG/lg32rvPAYnwyJjo6GSqWCm5tbsX3h4eFo27atBKrkT2ZmJmxtbYttj4uLQ0xMDPz8/CRQJX/Yt/LBvlVeOPAxDMMwFgWP8cmYiAhg5UogMhKIi0uDu7sDXn0VGDsWqF5danUyhQj480/gp5+ABw8Q//gxXmjYEOjUCRg5EqhSRWqFskSnAw4fBtauBaKjgfj4JNSt64IePYChQ4ECT9QwBcnJAfbuBTZuBGJikJSQAJcGDYC+fYEBAwB73GSYJQAAGQhJREFUe6kVMnrgHp/MIALWrwfmzAGiooCMDHFTysXeXnzu3h2YMQMICJBOq6zIyRGthLlzgbg4QKMRZuaiVgvj+vUDpk8HGjeWTquMyMwEli4Fvv0WSEkBik444uAgbBs6FJgyBfDykkan7EhLAxYsAJYsESampBTen9tSGDUK+M9/AD3DFox0cOCTEVlZ4gazb5+4rkpCoRBBcN06oH9/8+iTLWlpQJ8+wF9/iYBXElZWgJ0dsHMn0LmzefTJlMREIDQU+OcfID295GNVKhEEf/8daNPGPPpky5MnQIcOwP37omVaEtbWgIsLcOwY4ONjDnWMEXDgkwlEwOuvAwcOlH4TKoi9PbB1K9Czp+m0yRqtVqQxz54t/SZUELVa3MVfecV02mRMejoQFATcuCEaXMbi6AiEhwP+/qbTJmsSE4HmzUU6Rqs17jsKBeDsDFy4ANSrZ1p9jFFwza1M+P57cR8uGvSOHhXbUlLEq+jCDOnpwJtvAjEx5tMqK2bOBM6d0x/03nwTuHZN5O/u3Ckc5DQaoEeP4ikqC+GDD4BbtwoHvdxzLPel1YpMXkFSU4GQEMBilzccNgx4/Lh40PPyEmN9CQniYvzuO5FdAESrNiUF6NKlcPqdkQwOfDJApxNjeoaydO+/Dzg5iZe3t/7vL19uWo2GqFOnDvz8/BAYGIiWLVua98czM8UNRp9xnTuL8b6RI4Vx7doBd+8WPiYnB9iwwTxai5CYmIj+/fvD29sbTZo0walTp8z220lJohajaFsh9xxzchJDUunpwLZtxb+fng7s2mUerQW5efMmAgMD817Ozs5YtGiR+QRERwN//CHOu6L88AMQGwvUrAkEBgLBwcD48fn7dTqRIj1+3Hx6C7Bw4UL4+PjA19cXgwYNQkZZsiP/H5FoxhimAAcPEjk6EonmYOHX0aNEo0fr31fwVbUqUXa2+bV7eXnRs2fPzP/DRES//GLYuPBwolGjSjeuTh0inc7s0ocNG0Y//vgjERFlZmbS8+fPzfbbixcTqdUl2zJsGFFkpOH9LVuaTa5etFotubq60v379833o//5D5GtrX5Drl0j6tYt//O8eUTLlxc+RqEg6t7dfHr/JTo6murUqUMajYaIiN544w1as2aN2XXICe7xyYDFi4tX0xXk66+BZ89ElX5wsP5jtFrRGLUoFi3Sb5xSCbRsKZ75uH1bjMd8950oailKXJx4bsSMJCUl4cSJExg9ejQAMQNIFTM+ZrF4cek1QMOHAz//bHj/1avFO9Dm5PDhw6hfvz68zFlmuny5/t4eIM7FgQPFoHutWkC3bmLAviBE4pkRAyvemxKtVov09HRotVpoNBrUqlXL7BrkBAc+GXD7tuF9kyeL8XB3d1Gt/7//6R8fz84G7t0znUZDKBQKhIaGokWLFli5cqV5f9zQmmeuroCNjSh3ffVVkXpq1gyYOrX4sQqF2Y27d+8eqlevjpEjR6JZs2YYM2YM0kor461AHj0qeb+np2hgrVtn+BgbG2nOt1w2b96MQYMGme8Hs7JEjtgQJ06Iqs3kZGHwuXOicrgotrYiZWpG3N3d8emnn8LT0xM1a9aEi4sLQkNDzapBdkjd5dTHggULCIAFvR6UmpHLfe3fT/T++/r2ZREwyWQax40bR2lpacX+VtHR0URE9PTpU/L396fjx4/n7cvKyjKpb4mGTKpSRQgYNix/W79+RBcuFDs2BaCRJtQ4depUysnJKeRZREQEWVlZ0enTp4mI6MMPP6SpU6fm7Y+Pjzfx+ZZT4jk2ZQrRsWOlnYvPCXjNZBqXLVtm8P6QmZlJL7zwAj158qTQ9sjISJPpqQJQhiEzFAqi+/eJPv+cyMaGqFo1op07iebOLXbsc4CCTPi3/emnn4r5lZCQQB06dKDY2FjKysqi3r170/r16/P2HzhwwGR6GjRoQMnJyQb/llIhy8BnaTRtalzQA4j27SP64IPi2+3tiUq4V5iFGTNm0Pz58833g25uho16+JDorbfyP/ftqzfwkbMz0bZt5tNMRDExMeTl5ZX3+cSJE9TdjGM/9vYln2M3bxKNHFnyMS4uYvxZCnbu3EkhISHm/dHsbBHg9JnxwgviGGfn/G29exP9/XfxY52ciP75x6zSt27dSqNGjcr7vG7dOnr33XfNqkFucKpTBjRrJoaliuLiIh4wtrUVldGDB4vixKJDB4DY36SJ6bUWJC0tDSn/Pg6QlpaGP/74A76+vuYTUNIDwWvWiJr96tXFNGUTJgB79hQ/TqsFmjY1nUY9uLm5wcPDAzdv3gQgxquamlFDgwaG9730kkir66vmLEhmJtCoUcXqMpZNmzaZN80JiCf4DY2LxceLAc933xUXoouLGCS9cqX4sdnZIpdsRjw9PXH69GloNBoQEQ4fPowm5r5ZyA2pIy9DFBFB5OBQvHH44otEZ88SJScTPX9OdOoUUefO+hudtWubvzgxMjKS/P39yd/fn5o2bUqzZ882r4B9+wxXdapURN9/L4yLiRGljPoq8iRa2PfixYvUokUL8vPzo969e1NCQoLZfvvnnw3btny52F9a5sHcHa5cUlNTqVq1apSYmGj+H//vfw2XwwYEiC5wQgLRs2dEW7YQ1ahR+BgrK6KhQ82vm4imT59OjRs3Jh8fHxo6dChlZGRIokMu8MwtMqFJk+IPpxuLWi0qPz/8sGI1yR6dTjw3FRtbvu87OYmKoYEDK1aXzMnIEB3hkiqJS8LREfjtNwuc8e35c9HrK+8zcGq1mPYmMLBidTFlhlOdMmH2bHFdlAdbW5FZsTiUSjHhtIND2b+rUIjA169fxeuSOXZ2IvNbnvNNpQJq1wY6dqx4XbKnalVxoZXHOBsbEfA46MkC7vHJiLAwMXVZac9YFcTBQUwG0aKF6XTJGiIxO8u2bcYblxv0IiKkG6iSGJ0O6NULOHLE+LlhlUqgWjXg0iUxDmiRZGUB7dsDFy8a3/NTqURm4tIlYSAjOdzjkxFz5wKTJokGpUJR8rF2dqJm4+hRCw56gDBq9WoR/Izp+dnbizzf6dMWG/QAEcR++w3o3ds429Rq0dM7d86Cgx4gem6HDokHHY01rlEj0cjioCcbOPDJCIVCzLl8+LBYZcfOrvg6lk5OIuMSFibGBFu1kkSqvFAqxaJy//ufmAjYzk7kfwvi5CQebJ8xA7h+3fwlsDLE2hr45RexukdwsLDNxiZ/v0IhxvM8PIBvvhHLF/F6fBDBbO9esWpvq1biIlUVWNNbqRRBsUEDMU3OuXPi3GNkA6c6ZcyzZ8DmzcCDB0B8fCZq1bJFy5bAa68Vvs6YIjx6JO7mUVHIfP4cth4eok6/Sxf9z40wAMRMLL/+KuxLTs6Eh4ct2rcXmb3SMhAWzfXrYpaWmBikJiXBsV49ca61acPGyRQOfAzDMIxFwc1fhmEYxqLgwMcwDMNYFBz4GIZhGIuCAx/DMAxjUXDgY/5fEx8fL7UEhmFkBgc+GXKgwPILSUlJGD16NPz9/TF48GA8ffpUQmXy5rPPPkNcXBwA4Ny5c6hXrx7atGkDLy8vHD9+XGJ18qV58+aYPXs2IiMjpZZSqfi/9u49KKrrjgP4d9lVIOzGaBW0MmKbBseQKGrR+EqyY3EIajSVNrXRSmwsldTBRzXWYXw0TmNUUJwYHIlxGhJTZYwmYAa1JdqWahQVa4JBYgV5aCwJPlYg7OP0j2N47bLg1r33Mvf7mdkR9l64P3/cvb97zj333KKiIlitVsyePRuVlZWIjY1Fr169EBMTg7Nnz6odHnnBwqdBK1eubP566dKlGDBgAHJzcxETE4OkpCQVI9O2gwcPom/fvgCAZcuWYc+ePfjyyy9x5MgRLF26VOXotKuurg43btyA1WrF6NGjsXnzZtTU1KgdluYlJydj+fLlmDJlCsaNG4ekpCTcvHkT69evR3JystrhkRcsfBpXVFSEdevWISIiAosXL0Z5ebnaIWmWw+GAw+EAADQ0NCDm7rQ2kZGR+Pbbb9UMTdN69+6NTZs24cqVK0hLS0NZWRlGjhwJq9WKHTt2qB2eZtntdjzzzDOYNWsWDAYDEhISAACTJk1Co69PcCBFsPBp0PXr15Geno60tDTcunULrecYcLlcKkambcnJyYiPj0dBQQHi4uKQkpKCY8eOYfXq1YjmrPhdMnHiRLz55puorq7GK6+8guPHj6sdkmYFBQXh8OHDyMnJgcFgwIEDBwAAx44dg9FoVDk68oYTX2nQ/Pnzm59sPnfuXNTW1qJfv364du0aD+BeLFy4EI8//jgyMzNx8eJFOBwOlJWVYcaMGUhNTVU7PM2K9DBZt9FoRFxcHOLi4lSIqHvYvn07li9fjoCAABw6dAiZmZlITEzEwIEDkZWVpXZ45AWnLNOoL774AtXV1RgzZgzMZnPz+/n5+TwYecG8+YZ5882FCxdQU1PDvHU3Kj35nbzYunWriIyMFNOnTxcRERHiwIEDzctGjBihYmTalpGRwbz5gPubbzIyMsSQIUOYt26IXZ0atGPHDpw+fRpmsxnl5eVISEhAeXk5UlJS2lzvo7aysrKYNx9wf/NNVlYWioqKmLduiIVPg1wuV3O3yeDBg3H06FEkJCSgoqKCHygvmDffMG++Yd66L47q1KCwsDAUFxc3f282m5GXl4fa2lqcP39exci0jXnzDfPmG+at++LgFg2qqqqCyWRC//793ZYVFhZi/PjxKkSlfcybb5g33zBv3RcLHxER6Qqv8WnZf/4D/PnPwKVLaKitRXB4ODBuHPCLXwAPPKB2dJp14QLw7rtAeTnwzTcNGDQoGE89BcycCQQGqh2dhhUXA7t3A5WVaLh1C8GDBgGxscCzzwImHio8EQL49FMgJweoqgLq6m4hMvJBxMcDcXFAAC8maRJbfFqUnw+89hpw8iTgcgFNTS3LzGb5aUtMBH7/e2DwYLWi1BQhgP37Zdo+/xyw24G7s5cBACwW+W9SErB4MfD976sTp+Y4ncBf/gKsXy9PtBob5T73HYtFFr2FC+Xr7lyoetfUBLzzDvD668DVq0B9vdwHv2M2y3PTJUuABQuABx9UL1Zyx8KnJS4XkJIC7NoF3Lnjfd0ePWTz5aOPAKtVmfg0ym4H5s6VqegsbT17ygPSkSPAj3+sTHya1dAAPPcc8M9/dp64oCB59D56FBg6VJHwtOrmTdma+/e/ZcHzJjgYCA0Fjh0DIiKUiY86x8KnFUIAv/2t7KPr7NPU2gMPyBbixIn+i03DXC7ZhXn48L2lzWyWx/vhw/0Xm6bZ7cCkScCpU7KV1xUGgyx+p08DDz/s3/g0qqEBeOIJoLQU6Oq850Yj8L3vyZ7kAQP8Gx91DXugtSI723PRe/nlloPTrl3uP1dfD0yZAuj0gaubNrkXvZ49gbfektf4bt0Czp6VZ+it2WzAT34iD2S6tGIFUFTkXvSys4GaGtmsKS0Ffv3rlmVCALdvy4LpdCobr0bMnw9cvNhx0fvRj+Q+lZ3d8p7TCXzzDRAf37Y7lNTDwqcFQgBr1nhustTUAOvWAW+/3fHPOxzAzp1+C68j8+bNQ2hoKB577DHFtw3I//brr7unzWQCKiuBp54CevUCUlOBvXvdu5oaG+X7SqusrITVasWjjz6KqKgoZGRkKBvAnTvA9u2eq/5rr8nrxr16yUEt69YBI0e2LHe55FG81cOSldLY2IjRo0dj+PDhiIqKwurVqxXd/vXrwL593hvI27bJ89T2HA5ZME+e9F983jidTowYMQJTp05VJwCNYeHTgn/9S36qPNm/H/jwQ+8tuoYGID297aAEBSQmJrZ5WrzS8vJkj1179fXA2rVARYU8pzh4ELh8GRg1qu16NpssnEozmUxIS0tDSUkJTpw4gW3btqGkpES5AHbvlt2WnpSUtAymEkK+2ndr3r6tSuICAwNRUFCAc+fOobi4GPn5+Thx4oRi2+/sgQvPPw/cuAH87W+elzc2yh4KNWRkZGCozq/NtsbCpwVbttzbBSpP6uuBgoL7E08XPfnkk+jTp4+i22wtPV0egzsTGgpERsrRnu1duQKcO3f/Y/NmwIABGHm3FWWxWDB06FBUV1crF0BamvfBLNu2yeWlpXLI4scfu69z6pRMnoIMBkPzFGF2ux12ux2Gjgq4H2zd2nFrz2IB/vhHOYqzIy6XPFm7edM/8XWkqqoKBw8exEsvvaTshjWMhU8LPv/8/+/8t9uBsrL7E083cfFi5+uYTMB778nbIUtL3ZcbjV37Pf5SXl6Os2fPYsyYMcpttKLC+/KXX5ZH8gkTgA8+8HxBKzAQuHTJP/F54XQ6ER0djdDQUMTGxiqWN7sd+O9/O17+6qvyakNn5y89eshueCUtWrQIGzZsQABvKmzGTGiBzfb//w6HQ47k0JHOBqYYDHKQQVMT8LvfeV7H6VQvbTabDTNnzsSWLVvwoJI3enVlOKLLBRQWAuHh8ka09oRQJXFGoxHFxcWoqqrCyZMn8dlnnymyXZtNFi1Phg+XA6U2b+789wQEKJu2vLw8hIaGYlT7fn69U/o5SJ1xOp1i1apVAoBuXiUtV1M6fr36qhC7dnW4/A4gFvgxxqSkJHHnzh23v9fly5dFVFSUx79lU1OTn3N31WvK3n5biIICIYKCvKX2hgB+5rcYU1NThdPp9JibyZMni7S0NLdlX3/9tV/zVt+V/e27V1aWEFu2uL1fBwirH2Ncs2ZNp8eKtWvXio0bNzZ/f+bMGT/mzSQAp8cUpaQIYbMJcfWqfN2+LUR9vRCnT3e0vz3mtzh37tzZJkcrVqwQAwcOFBERESIsLEwEBweLF154oc06+fn5founf//+Ho8batNc4dOln/5UCIPB84HHaBQiMFCIP/1JiHfekV8bje7rWSxCHDqkeOjeCp+/jR/f8fE6M1OI48eFCAnxflwPCRHizBll43a5XGLOnDkiJSVF2Q1/JzLSczL69RPi+edlUgIChJg8WR7Rp01zXzcoSIjyckXDvn79uqirqxNCCFFfXy8mTJggcnNzFdt+v36e0xYcLERYWMtr40YhcnKE6NvXfd3AQCHu/hcU98knn4gpU6aos3GNYVenFixaBISEeF6WmiqvqP/hD8CcOfLr1FT39YKD5f1VCpo1axbGjh2L0tJShIeHY6fCt1QsWdIyFVlrgwbJuQCio4Fr1+QAmNu3gV/+0n3d8HBgxAj/x9paYWEhsrOzUVBQgOjoaERHR+NjTwNI/GXJEs/7mxCyW1NOOimHIC5aBOTmuq87apTiU5FcvXoVVqsVw4YNQ0xMDGJjYxUdnr9woZzApr2GBuCrr1peNpv8mNbWtl3PYJC33D70kDLxUsc4c4sWCAH84AedDzroSHAwsGqVvClZRxwOOWKzrs63nzeb5Ui9F1+8v3Fpns0mE+fr3fsWixwxNG3a/Y1L4776St7i2NWJbtoLCZFT5Y0de1/DIh+wxacFBgOwerXvT1wwmdrOsKETJhOwbJnvaQsMlPde6Y7ZDPzmN/KE6V4ZDLLJEh9//+PSuLAwYMYMz62+zphM8nbIJ564/3HRvWOLTyuEkMVrz557n6szL0+3E1W7XMD06fKm4XtpwISEAH//e9tJSXTFbgeefho4c+be5uq0WORUZ4884tfwtKq+Hhg9Wt451PqhKd4EBAB9+si5OgcO9G981DVs8WmFwSCnhpgzp+Prfa2ZTHK9Awd0W/QAeVDZt082QLqSth495Gxcf/2rjoseIBNx+DAwfnzXEhcYKGdaLizUbdED5HnmP/4BDBvWtZ6GoCBZ7D79lEVPS1j4tMRoBDIzgfffl30iQUHuNw+FhMguqhdflKeQsbHqxKohPXvKB4Hu3CkHtAQHy1S2ZjbL14IFwPnz7HICIPel/HzgjTeAIUPk9+1vcrZY5JnCkiVyogWV5mXVkt695ZM9NmyQ1/zMZvcZ4CwWeZ6wcqWcGeiHP1QlVOoAuzq1rKxMPpHh0iXk7d2Lqb/6lbwyPnu2/LSRR+fPy4eEVlQAOTm5SEychqefBn7+c98ua+nGqVPNT2DP3bcP0+bNAyZPls/s69lT7eg0SQhZBPfskfPJ79//IebPn46pU+UIzvYnYKQNLHzdhMFgAP9U94558w3z5hvmrXtgVycREekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKCx8REekKC183kZGRoXYI3RLz5pvNmzerHUK3lJ6ernYI1AUGIYRQOwgiIiKlsMVHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES6wsJHRES68j/dTBXFRDjOGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib as mplt\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "Lx = 4\n",
    "Ly = 6\n",
    "#fig = plt.figure()\n",
    "#draw2DLattice(Lx,Ly,fig)\n",
    "draw2DLattice(Lx,Ly)\n",
    "#fig.savefig(\"/tmp/graph.png\",dpi=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.) Generating input file \n",
    "\n",
    "In this work, we will be using the [DEHam](https://github.com/v1j4y/DEHam) code to solve the Hamiltonian using Exact Diagonalization method. \n",
    "\n",
    "The input file for the code is as follows:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "  1. The DEHam program requires an input file which \n",
    "   has the topology of the Hamiltonian and the various parameters\n",
    "   as explained below in a sample inputfile:\n",
    "\n",
    "```python\n",
    "8     \t\t\t\t\t\t\t\t# The number of orbitals (total)\n",
    "1     \t\t\t\t\t\t\t\t# The largest number of non-zero elements per row (Multiple of Ndet)\n",
    "1     \t\t\t\t\t\t\t\t# The total number of processors used in parallel (Multiple of Ndet)\n",
    "1     \t\t\t\t\t\t\t\t# The number of holes\n",
    "0     \t\t\t\t\t\t\t\t# The isz (ms-1/2) value\n",
    "true     \t# Restrict the hole to the 1'st (i.e. half of natom) Family of states. *false* for no restrictions\n",
    "1,2,3,1,2,3,4,5,6,7\t\t\t\t\t# (L1) The topology of the system is specified here\n",
    "2,3,4,8,7,6,5,6,7,8\t\t\t\t\t# (L2) first and second line contain the two sites linked\n",
    "1,1,1,2,2,2,2,3,3,3\t\t\t\t\t# (ltype) third line contains the type of link (1 for t or J, 2 for K and 3 for none)\n",
    ".1430,-0.20,0.0000\t\t\t\t\t# The three types of links this line gives Jz, Kz\n",
    ".1430,-0.20,0.0000\t\t\t\t\t# The three types of links this line gives Jx,y; Kx,y\n",
    "-1.00,0.0,0.00\t\t\t\t\t\t# This line gives t(i,i+1), t(i,i+2), t(i,i+3)\n",
    "0.,0.,0.,0.,0.,0.,0.,0.,0. \t\t    # Energy of each orbital + additional energy\n",
    "2                   \t\t\t\t# The total number of roots\n",
    "1                   \t\t\t\t# I   The position of the first\n",
    "1                   \t\t\t\t# I   SBox\n",
    "1                   \t\t\t\t# I\n",
    "1                   \t\t\t\t# I\n",
    "1                   \t\t\t\t# II  The positions of the second\n",
    "1                   \t\t\t\t# II  SBox\n",
    "1                   \t\t\t\t# II\n",
    "1                   \t\t\t\t# II\n",
    "1                   \t\t\t\t# III\n",
    "1                   \t\t\t\t# III The positions of the third\n",
    "1                   \t\t\t\t# III SBox\n",
    "1                   \t\t\t\t# III\n",
    "1                   \t\t\t\t# positio of the hole\n",
    "0                   \t\t\t\t# fix the position of the first hole during the CI\n",
    "0                   \t\t\t\t# fix the position of the second hole during the CI\n",
    "0                   \t\t\t\t# Print the wavefunction. It is stored in the FIL666 file after the run\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### a.) Script for $t-J$ input file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]\n",
      "[2, 3, 4, 1, 6, 7, 8, 5, 10, 11, 12, 9, 14, 15, 16, 13, 18, 19, 20, 17, 22, 23, 24, 21, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 1, 2, 3, 4]\n",
      "[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "48\n"
     ]
    }
   ],
   "source": [
    "l1, l2, ltype = getInput(Lx,Ly)\n",
    "print(l1)\n",
    "print(l2)\n",
    "print(ltype)\n",
    "print(len(l1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.) Misc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Draw Lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function that draws a 2D Lattice\n",
    "\"\"\"\n",
    "def draw2DLattice(Lx,Ly):\n",
    "    n=Lx+2\n",
    "    m=Ly+2\n",
    "\n",
    "    G = nx.grid_2d_graph(n,m)  # Lx, Ly grid\n",
    "\n",
    "    pos = dict( (n, n) for n in G.nodes() );\n",
    "    #labels = Dict( ((i, j), i + (Ly-j) * (Lx+1) ) for (i,j) in G.nodes() );\n",
    "    #edges=Dict()\n",
    "    #for (i,edge) in enumerate(G.edges)\n",
    "    #    edges[edge]=i\n",
    "    #    print(i,\"\\t\",edge,\"\\n\")\n",
    "    #end\n",
    "\n",
    "\n",
    "    for j in range(0,m-1):\n",
    "        G.remove_edge((0,j),(0,j+1))\n",
    "        G.remove_edge((n-1,j),(n-1,j+1))\n",
    "\n",
    "    for i in range(0,n-1):\n",
    "        G.remove_edge((i,0),(i+1,0))\n",
    "        G.remove_edge((i,m-1),(i+1,m-1))\n",
    "\n",
    "\n",
    "\n",
    "    edgescenter = dict()\n",
    "    countedge = 1\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            edgescenter[((i,j),(i+1,j))] = countedge\n",
    "            countedge += 1\n",
    "\n",
    "\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            edgescenter[((i,j),(i,j+1))] = countedge\n",
    "            countedge += 1\n",
    "\n",
    "\n",
    "    labelscenter = dict()\n",
    "    countlabels = 1\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            labelscenter[(i,j)]=countlabels\n",
    "            countlabels+=1\n",
    "\n",
    "\n",
    "    color_map=[]\n",
    "    countcolor=0\n",
    "    for i in range(0,n):\n",
    "        for j in range(0,m):\n",
    "            if i==0 or j==0 or i==n-1 or j==m-1:\n",
    "                color_map.append(\"white\")\n",
    "            else:\n",
    "                if(countcolor%2 == 0):\n",
    "                    color_map.append(\"blue\")\n",
    "                else:\n",
    "                    color_map.append(\"red\")\n",
    "                countcolor+=1\n",
    "        countcolor+=1\n",
    "\n",
    "\n",
    "#    nx.draw(G, pos=pos, labels=labelscenter, node_color=color_map, font_color=\"white\", ax=fig.add_subplot(111));\n",
    "#    nx.draw_networkx_edge_labels(G,pos,edge_labels=edgescenter,font_color=\"black\", ax=fig.add_subplot(111));\n",
    "    nx.draw(G, pos=pos, labels=labelscenter, node_color=color_map, font_color=\"white\");\n",
    "    nx.draw_networkx_edge_labels(G,pos,edge_labels=edgescenter,font_color=\"black\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Generate Input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Function that generates the input\n",
    "\"\"\"\n",
    "def getInput(Lx,Ly):\n",
    "    n=Lx+2\n",
    "    m=Ly+2\n",
    "\n",
    "\n",
    "    labelscenter = dict()\n",
    "    countlabels = 1\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            labelscenter[(i,j)]=countlabels\n",
    "            countlabels+=1\n",
    "#     print(labelscenter)\n",
    "    edgescenter = dict()\n",
    "    countedge = 1\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            edgescenter[((i,j),(i+1,j))] = countedge\n",
    "            countedge += 1\n",
    "    for j in range(1,m-1):\n",
    "        for i in range(1,n-1):\n",
    "            edgescenter[((i,j),(i,j+1))] = countedge\n",
    "            countedge += 1\n",
    "\n",
    "#     print(edgescenter)\n",
    "    l1 = []\n",
    "    l2 = []\n",
    "    itype = [1 for x in range(0,len(edgescenter))]\n",
    "     \n",
    "    for x in edgescenter:\n",
    "#         print(x[0],\"\\t\",x[1])\n",
    "#         print(l1)\n",
    "#         print(l2)\n",
    "        if Lx + 1 not in x[1] or Ly + 1 not in x[1]:\n",
    "#             print(\"-->\",x[1])\n",
    "            l1.append(labelscenter[x[0]])\n",
    "#             l2.append(labelscenter[x[1]])\n",
    "            if x[1][0] > Lx and x[1][1] <= Ly:\n",
    "                l2.append(labelscenter[(x[1][0]-Lx,x[1][1])])\n",
    "            elif x[1][0] <= Lx and x[1][1] > Ly:\n",
    "                l2.append(labelscenter[(x[1][0],x[1][1]-Ly)])\n",
    "            elif x[1][0] > Lx and x[1][1] > Ly:\n",
    "                l2.append(labelscenter[(x[1][0]-Lx,x[1][1]-Ly)])\n",
    "            else:\n",
    "                l2.append(labelscenter[x[1]])\n",
    "        else:\n",
    "            l1.append(labelscenter[x[0]])\n",
    "            if x[1][0] > Lx:\n",
    "                l2.append(labelscenter[(x[1][0]-Lx,x[1][1])])\n",
    "            elif x[1][1] > Ly:\n",
    "                l2.append(labelscenter[(x[1][0],x[1][1]-Ly)])\n",
    "\n",
    "#    print(l1,\"\\n\")\n",
    "#    print(l2,\"\\n\")\n",
    "#    print(itype,\"\\n\")\n",
    "    return l1,l2,itype"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6],\n",
       " [2, 1, 4, 3, 6, 5, 3, 4, 5, 6, 1, 2],\n",
       " [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getInput(2,3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "header=r\"\"\"\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}