1d_hubbard_totalenergy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cb509096-42c6-4a45-8dc4-a8eed3116e67",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pymf\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "396d935c-146e-438c-878b-04ed70aa6d63",
   "metadata": {},
   "source": [
    "To simulate infinite systems, we provide the corresponding tight-binding model.\n",
    "\n",
    "We exemplify this construction by computing the ground state of an infinite spinful chain with onsite interactions.\n",
    "\n",
    "Because the ground state is an antiferromagnet, so we must build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian in is:\n",
    "$$\n",
    "H_0 = \\sum_i c_{i, B}^{\\dagger}c_{i, A} + c_{i, A}^{\\dagger}c_{i+1, B} + h.c.\n",
    "$$\n",
    "We write down the spinful by simply taking $H_0(k) \\otimes \\mathbb{1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5529408c-fb7f-4732-9a17-97b0718dab29",
   "metadata": {},
   "outputs": [],
   "source": [
    "hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))\n",
    "h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "050f5435-6699-44bb-b31c-8ef3fa2264d4",
   "metadata": {},
   "source": [
    "To build the tight-binding model, we need to generate a Hamiltonian on a k-point and the corresponding hopping vectors to generate a guess. We then verify the spectrum and see that the bands indeed consistent of two bands due to the Brillouin zone folding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b39a2976-7c35-4670-83ef-12157bd3fc0e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAGyCAYAAAD9IyA0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABnO0lEQVR4nO3deVxUZeP+8c8MCEQIZIiIomIqILigkvtSFu5pmZaVqZW5ZmqWj9qT2mZZpmWalaVZmS3mlppabrnmhooCkqG4Ia4gIiLM/P7g53zjcQkUODNwvV+veT0PZ87MXJDMXJxzn/s2Wa1WKyIiIiKSJ2ajA4iIiIg4EpUnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB2ejAzgKi8XC8ePHKV26NCaTyeg4IiIikgdWq5ULFy7g7++P2Vwwx4xUnvLo+PHjBAQEGB1DREREbsGRI0eoWLFigTyXylMelS5dGsj54Xt6ehqcRkRERPIiNTWVgIAA2+d4QVB5yqOrp+o8PT1VnkRERBxMQQ650YBxERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB4csTxMmTCAiIoLSpUvj6+tLly5diIuL+9fHrVu3jvr16+Pm5kbVqlWZMWNGEaQVERGR4sQhy9O6desYNGgQW7ZsYdWqVWRlZREZGcnFixdv+JiEhATat29P8+bN2bVrF6NHj2bIkCHMnz+/CJOLiIiIozNZrVar0SFu16lTp/D19WXdunW0aNHiuvuMHDmSxYsXExMTY9vWv39/du/ezebNm//1NVJTU/Hy8mLdunXcfffduLm54erqioeHB56enpjNDtlDRUREig2LxcLp06dxdnamTJkywP99fqekpODp6Vkgr+NcIM9isJSUFADbD+p6Nm/eTGRkZK5tbdq04YsvvuDKlSuUKlUq132XL1/m8uXLtq9TU1MBaNmy5XWf32QyYTabcXJywsXFBXd3d7y8vChTpgzlypWjQoUKhISEEB4eTr169XB3d7+l71VERKSkOH/+PNu2bWPPnj3ExsZy/PhxkpOTOXfuHKmpqaSnp3PlyhWysrKwWCy2x7Vu3Zrffvut0HI5fHmyWq0MHz6cZs2aERYWdsP9kpKSKFeuXK5t5cqVIysri9OnT1O+fPlc902YMIHx48fnK0d2djbZ2dlkZmaSlpZGcnLyDfc3m824u7tTrlw57rnnHmrXrk2TJk2477778Pb2zvPrioiIOLKkpCR+//13tm7dyt69e0lISOD06dOkp6dzqyfHLly4UMApc3P48jR48GD27NnDhg0b/nVfk8mU6+ur/1H+dzvAqFGjGD58uO3r1NRUAgICSElJwcPDg8zMTDIyMkhNTSU5OZlTp05x6tQpTp8+zalTpzh27BhJSUmcOnXK1pDT0tK4cuUKkHNoMS0tjbS0NA4ePMjKlSttr+Xi4kK5cuUICQmhSZMmdOrUiXr16t3Sz0dERMQeWCwWNmzYwNKlS9m6dSsHDhzg1KlTZGVl/etjXVxcKF26NF5eXtx1112ULVsWPz8//P398fX1xcfHB19fX8qVK2f7ujA5dHl64YUXWLx4MevXr6dixYo33dfPz4+kpKRc25KTk3F2dubuu+++Zn9XV1dcXV2v+1xmsxk3Nzfc3Nzw9vamUqVKec6clZVFbGwsO3bsYOfOnezbt4+///6b5ORk24D3zMxMjhw5wpEjR1i5ciXjxo3DyckJPz8/6tSpQ9u2benRo0eh/+MQERG5VYcPH2bu3Ln89ttvREdHc+rUqRseSTKZTJQuXTrX2Zh69erRoEEDAgMD7W5csUMOGLdarbzwwgssWLCAtWvXUr169X99zMiRI1myZAn79++3bRswYABRUVH5GjBekAPO/pfFYiEqKoply5axceNGYmJiSEpKyjX26p88PDyoVasWHTp04Nlnn8XPz69QcomIiPybgwcPMnPmTH799Vfi4uK4dOnSdfe74447qFChAmFhYTRv3pwOHToQFBRUaLkK4/PbIcvTwIEDmTt3LosWLcr1A/fy8uKOO+4Ack67HTt2jDlz5gA5UxWEhYXRr18/+vbty+bNm+nfvz/fffcdXbt2/dfXLIrydCNJSUl8//33rFixgqioKE6ePJlrYNxVHh4eRERE0LNnT5588klcXFyKNKeIiJQcaWlpfPbZZ3z//ffs2bOHjIyMa/ZxdnamQoUK1K9fn/bt29OtW7ci/wxVefr/rjdGCWDWrFn07t0bgN69e3Po0CHWrl1ru3/dunUMGzaMffv24e/vz8iRI+nfv3+eXtPI8vS/rp43nj17NuvXr+fQoUNkZ2fn2sdkMhEQEEC7du14+eWXueeeewxKKyIixcX27dv54IMPWLNmzTVDYSBnbFK1atV44IEH6NOnD3Xr1i36kP9D5clA9lSermfDhg18/vnn/P777xw7duya+++66y7uv/9+XnrpJRo3bmxAQhERcUS//PILU6dOZePGjddMRm0ymQgMDKRt27YMGDDgple9G0XlyUD2Xp7+KT09ndmzZzN37lx27NhxzaFUd3d37r//fl599VUaNmxoUEoREbFXy5YtY+LEiWzevJnMzMxc95UuXZpGjRrRp08funXrhrOzfV97pvJkIEcqT/9r7dq1TJkyhbVr19omFL3K09OTNm3aMHbsWEJDQw1KKCIiRtuwYQNvvfUWa9euveaP7rJly9K2bVtGjBhB7dq1DUp4a1SeDOTI5emf9u3bx1tvvcWyZcuuKVJ+fn48/fTT/Pe//8XDw8OghCIiUlSSkpL473//y08//cT58+dz3efr60vXrl0ZPXr0v04HZM9UngxUXMrTP23bto0333yT3377jfT0dNt2k8lEnTp1GDNmDI8++qiBCUVEpKBZLBamT5/ORx99RHx8fK77vL296dSpE2PHji02FxqpPBmoOJanf1q8eDFvvvkmO3bsyDUNgoeHBz169OCdd9656dqBIiJi3w4fPszw4cP55Zdfco1jKlWqFC1btuT1118vlhcUFcbnt31N2SmGeeihh/jzzz+5ePEi48aNo0KFCkDOPB6ff/45Pj4+NGjQgFWrVhmcVERE8uObb76hRo0aVKlShZ9//tlWnKpXr860adPIyMhg1apVxbI4FRYdecqj4n7k6Xr27NnDiBEjWL16da55pHx9fRkxYgQvvfSS3U2ZLyIikJGRwZgxY/j8889zLZLr6urKQw89xPvvv5+vpcUcmY48SZGqXbs2K1euJD09nXHjxuHr6wvkrAn4yiuv4O7uTq9evTh79qzBSUVEBHJOzXXs2BEPDw8++OADW3GqUqUK06dPJz09nR9++KHEFKfCovIk/8rFxYWxY8dy8uRJVq5cSa1atQC4fPkyc+bMwcfHh8jISA4fPmxwUhGRkmn79u00aNCAKlWqsHTpUrKzszGZTDRv3pxdu3aRkJDAgAEDdLaggOinKPny4IMPsmfPHg4dOkTHjh1xcnLCarWyatUqqlSpwr333svOnTuNjikiUiIsX76coKAgIiIi2LFjB5Bzaq53796cPXuW9evX28USKcWNypPcksqVK7NkyRJSU1MZMGAArq6uQM70B/Xr1yckJISNGzcanFJEpHiaP38+FStWpH379hw4cAAALy8vxo0bR3p6OrNmzcLb29vYkMWYypPcFnd3d9t59LFjx9oG48XGxtKsWTNCQkLYsGGDwSlFRIqHH374gQoVKvDoo4/a1jH18/NjxowZnD9/nrFjx+rUXBHQT1gKhNlsZty4caSkpDB16lTbnFCxsbE0b96c4OBglSgRkVv0008/UaFCBR577DGOHz8OQEBAAD///DMnTpygX79+BicsWVSepMANHjyYM2fOMH36dFuJiouLo3nz5tSuXZvo6GiDE4qIOIbff/+dypUr061bN1tpqlSpEosWLSIxMZGHH37Y4IQlk8qTFJoBAwZw5swZZsyYYStRe/fupVatWjRt2lRX54mI3MD27dsJCQnhgQceIDExEcg50rRo0SIOHz7MQw89ZHDCkk3lSQpdv379OHPmDO+//z533nknAJs2bSIwMJAOHTpcsxiliEhJdfjwYRo2bEhERASxsbFAzsTE33//PYmJiSpNdkLlSYrMSy+9RGpqKv/5z39wdXXFarWybNkyypYty5AhQ3KtqSciUpKkp6fz6KOPEhgYyJ9//gmAp6cn06ZN4+TJk3Tv3t3ghPJPKk9SpMxmMxMmTCA1NZVnn30Ws9lMVlYWU6dOxdPTk+nTpxsdUUSkyFgsFkaNGoW3tzfz58/HarXi4uLCf//7X86dO8fAgQONjijXofIkhnBxcWHmzJmcPHmS1q1bA3Dx4kUGDRpExYoVNUeUiBR7P/zwA2XKlOGdd97hypUrmEwmnnjiCVJSUnj99dc15YAd038ZMZSPjw+//fYbu3fvJiQkBIBjx47RrFkzWrZsyenTpw1OKCJSsOLi4ggNDeWxxx4jJSUFgKZNm3L06FG+/fZb3NzcDE4o/0blSexC7dq12b9/Pz/++KNtVtz169fj5+fH8OHDNR5KRBxeRkYG3bp1IyQkhP379wM5V9Bt2rSJDRs24O/vb3BCySuVJ7Erjz76KGfOnGH48OE4OTmRnZ3N5MmT8fHxYdmyZUbHExG5JZ9++il33XUXP/30E1arFTc3N6ZMmUJiYiKNGzc2Op7kk8qT2B2z2cykSZNISkqiZcuWAJw7d44OHTrQtGlTncoTEYcRHx9PUFAQ/fv3JyMjA5PJxOOPP865c+d48cUXjY4nt0jlSeyWj48Pa9euZc2aNZQtWxbImR+qfPnyjBs3zthwIiI3kZWVRZ8+fQgKCrIt3Fu9enViYmL47rvvNK7Jwak8id1r1aoVSUlJvPLKKzg5OZGVlcX48eOpWLEiUVFRRscTEcll2bJl3H333cyePRur1YqrqyvTpk3jwIEDBAUFGR1PCoDKkzgEs9nMu+++y9GjR4mIiAByrsoLDw/n6aefJisry+CEIlLSpaWl0bp1azp06EBqaioAHTt25OzZs5qvqZhReRKH4ufnx59//sn333+Pu7s7AF9//TU+Pj6sWLHC4HQiUlJ99tln+Pj4sHr1agDKli3Lhg0bWLJkie29SooPlSdxSN27d+fcuXN07twZgJSUFNq2bUvbtm3JyMgwOJ2IlBTJycnUrVuXfv36cfnyZUwmE4MHDyYpKYmmTZsaHU8KicqTOCwXFxcWLlzIH3/8gY+PDwArVqzAx8eHRYsWGZxORIq7qVOnUqFCBXbv3g1AtWrViIuLY+rUqZodvJjTf11xeM2aNePkyZP07dsXk8nExYsX6dKlC23atNFRKBEpcMnJydSpU4chQ4aQlZWFk5MTb7/9NvHx8VSvXt3oeFIEVJ6kWDCbzXz22Wds374dX19fAFauXMndd9+to1AiUmCuHm3as2cPAMHBwRw6dIhRo0YZnEyKkkOWp/Xr19OpUyf8/f0xmUwsXLjwpvuvXbsWk8l0zS02NrZoAkuRqVevHidOnKBfv36YTCbS09Pp0qULnTt31hV5InLLzp8/T4MGDa452hQTE0PFihWNjidFzCHL08WLF6lTpw4ff/xxvh4XFxfHiRMnbDcdXi2ezGYzM2bMYPv27bbJNRcvXkzZsmXZuHGjwelExNF8//33+Pn5sWPHDgBq1Kiho00lnEOWp3bt2vHmm2/yyCOP5Otxvr6++Pn52W5OTk6FlFDsQb169UhKSuKJJ54Acv5ybNasGX379tVCwyLyrzIzM2nTpg2PP/647Uq6V199lbi4OB1tKuEcsjzdqvDwcMqXL0/r1q1Zs2bNTfe9fPkyqampuW7ieMxmM99++y2//fYbpUuXBmDmzJlUqlSJhIQEg9OJiL3auHEjZcuWZeXKlQCUL1+effv28cYbbxicTOxBiShP5cuX57PPPmP+/Pn8/PPPBAUF0bp1a9avX3/Dx0yYMAEvLy/bLSAgoAgTS0Fr3bo1p0+fJjIyEsiZnbx69ep8+OGHBicTEXszZMgQmjVrZvuj+bnnnuPo0aOEhIQYnEzshclqtVqNDnE7TCYTCxYsoEuXLvl6XKdOnTCZTCxevPi691++fJnLly/bvk5NTSUgIICUlBQ8PT1vJ7IYbM6cOTz33HNcuXIFgCZNmrBq1SrNAixSwh0/fpzmzZvz999/A+Dh4cGSJUto1aqVscHktqSmpuLl5VWgn98l4sjT9TRq1Ij4+Pgb3u/q6oqnp2eumxQPTz/9NImJidSoUQOATZs24evra1tWQURKnlmzZlG5cmVbcWrRogWnTp1ScZLrKrHladeuXZQvX97oGGIQPz8/4uLiGD58uG1izdatW/Piiy8aHU1EilBWVhadOnXimWeesU1BMG3aNNatW4ebm5vR8cROORsd4FakpaXx119/2b5OSEggKiqKMmXKUKlSJUaNGsWxY8eYM2cOAFOmTKFKlSqEhoaSmZnJN998w/z585k/f75R34LYiUmTJtG9e3ciIyNJTU3lo48+YuXKlWzcuJEyZcoYHU9EClFcXBwtWrQgOTkZgAoVKrBx40YqV65scDKxdw555Gn79u2Eh4cTHh4OwPDhwwkPD+e1114D4MSJEyQmJtr2z8zMZMSIEdSuXZvmzZuzYcMGli5dmu+pDqR4atiwISdPnrQt4hkbG0uFChX45ZdfDE4mIoXl448/pmbNmrbi1KNHDxITE1WcJE8cfsB4USmMAWdif9555x1Gjx7N1V+Lfv36MWPGDINTiUhBycrKon379qxatQrIWWD866+/pnv37gYnk8JSGJ/fKk95pPJUcuzZs4f77ruPs2fPAjlrV23evBlvb29jg4nIbYmPj6dp06acOnUKgKpVq7Jx40b8/PwMTiaFSVfbiRSB2rVrc+LECdtVNrGxsfj7++tqPBEH9sUXXxASEmIrTr179+bgwYMqTnJLVJ5ErsPFxYU1a9YwYcIETCYTly5donXr1lrLSsTBWCwWunXrxnPPPUd2djalSpXi+++/Z9asWUZHEwem03Z5pNN2Jde2bdto3bo1Fy5cAODee+/VZcwiDiA5OZl7772Xw4cPAzlX023atIlKlSoZnEyKkk7biRggIiKCpKQk6tWrB8Cff/6Jv78/MTExBicTkRv5/fffqVSpkq04de7cmcTERBUnKRAqTyJ54O7uzo4dOxgyZAgA586do1atWra5xETEfrz22ms88MADXL58GbPZzNSpU1m4cCFmsz7ypGDotF0e6bSdXLVo0SK6detmWxuvd+/eGj8hYgeysrJ44IEHWLduHQClS5dm7dq1tqPGUjLptJ2IHejcuTN//fWX7Sqd2bNnExYWRlpamsHJREquxMREKlSoYCtOtWrV4vjx4ypOUihUnkRuQaVKlThy5AitW7cGYN++fVSoUIF9+/YZnEyk5Fm+fDnVq1e3zRbet29f9uzZg4eHh8HJpLhSeRK5Rc7Ozvz222/897//BXIODdepU4dvv/3W4GQiJcf48eNp3749mZmZODk58dVXX/HZZ58ZHUuKOY15yiONeZKbWbZsGV26dLGNgxo4cCDTpk0zOJVI8WWxWGjfvj0rVqwAcsY3bdq0ibCwMIOTib3RmCcRO9W+fXvi4+MpW7YsANOnT6dRo0ZkZWUZnEyk+Dl79ixVq1a1FaeQkBCOHz+u4iRFRuVJpIBUrlyZo0eP0rhxYwC2bt1KQEAAx48fNziZSPGxc+fOXPM39ejRg/3792t8kxQplSeRAuTi4sKmTZsYOHAgAElJSdxzzz1s2LDB4GQijm/OnDlERERw8eJFTCYTkydPZu7cuUbHkhJI5UmkEEybNo2ZM2diNpvJyMigRYsWTJ8+3ehYIg5r2LBh9OrVC4vFgqurK6tXr2bo0KFGx5ISSuVJpJA8++yzbNmyhTvuuAOr1cqgQYN4/vnnjY4l4lAsFgv3338/U6ZMAaBs2bL89ddftGrVytBcUrKpPIkUooiICA4dOkSFChUA+Pzzz2nWrJkGkovkQWpqKtWqVWPNmjUA1K9fn6NHj1KxYkWDk0lJp/IkUsh8fX05dOgQTZo0AWDjxo1UrVqVs2fPGpxMxH7FxcUREBBAQkICAD179mT79u24uLgYnExE5UmkSDg7O7Nx40aeffZZAI4cOUKlSpWIjo42OJmI/Vm2bBlhYWGkpqYC8O6772oRbrErKk8iRWjmzJlMnjwZk8nExYsXqVu3LosXLzY6lojdmDJlCh07diQrKwtnZ2eWLFnCK6+8YnQskVxUnkSK2NChQ1m+fDmlSpUiOzubzp0788EHHxgdS8RwAwYMYNiwYVitVjw8PIiKiqJjx45GxxK5hsqTiAHatGnD3r17KV26NAAvvfQSAwYMMDiViDEsFgutW7dmxowZAFSsWJEjR44QGhpqcDKR61N5EjFIUFAQiYmJtiuHZsyYQevWrbFYLAYnEyk66enpBAcHs3r1aiDnCtWEhAS8vb2NDSZyEypPIgby9vYmISGBe++9F4DVq1cTHBxMenq6wclECt/Ro0cJCAggPj4egMcff5w///wTZ2dng5OJ3JzKk4jBnJ2d2bp1Kz169AAgPj6eSpUqkZSUZHAykcKzc+dOqlevbpuyY9y4cXz33XcGpxLJG5UnETsxd+5cxo4dC8CZM2eoWrUqe/bsMTiVSMH75ZdfuPfee8nIyMBsNvP111/b/u2LOAKVJxE7Mm7cOL788ktMJhOXLl2ifv36rFixwuhYIgVm+vTpPPTQQ2RnZ1OqVCnWrFnDU089ZXQskXxReRKxM3369GHVqlU4OzuTlZVFu3btmDlzptGxRG7byJEjGTRokG0qgr1799KiRQujY4nkm8qTiB1q3bo1UVFRuLu7Y7Va6du3L+PGjTM6lsgt69GjBxMnTgTAz8+PhIQEgoKCDE4lcmtUnkTsVGhoKAcPHqRs2bIAjB8/nr59+xqcSiR/LBYLLVu2ZN68eQCEhISQkJCAj4+PwclEbp3Kk4gd8/Pz49ChQ1StWhXIWd6lffv2mgtKHEJmZiahoaGsX78egFatWhEdHY2bm5vByURuj0OWp/Xr19OpUyf8/f0xmUwsXLjwXx+zbt066tevj5ubG1WrVrXNZCti79zd3YmPjyciIgKA5cuXExERQVZWlsHJRG7s/PnzVKlShdjYWCDntN2aNWswmx3yY0ckF4f8V3zx4kXq1KnDxx9/nKf9ExISaN++Pc2bN2fXrl2MHj2aIUOGMH/+/EJOKlIwzGYzf/75p22dr507d1KtWjVNpil2KTExkSpVqnDixAkARowYwdy5cw1OJVJwTFar1Wp0iNthMplYsGABXbp0ueE+I0eOZPHixcTExNi29e/fn927d7N58+Y8vU5qaipeXl6kpKTg6el5u7FFbtmAAQNsR07Lli3L/v37NX5E7EZ0dDT33nsvly5dAmDy5MkMHTrU2FBSohXG57dDHnnKr82bNxMZGZlrW5s2bdi+fTtXrly57mMuX75MampqrpuIPfjkk08YP348AKdOnaJq1aocPnzY4FQisGHDBurVq8elS5cwmUzMnTtXxUmKpRJRnpKSkihXrlyubeXKlSMrK4vTp09f9zETJkzAy8vLdgsICCiKqCJ58tprrzF9+nRMJhMXLlwgODiYqKgoo2NJCbZo0SJatmzJlStXcHJyYsWKFbYlh0SKmxJRniDn9N4/XT1b+b/brxo1ahQpKSm225EjRwo9o0h+DBgwgB9++AGz2UxGRgYRERGsXbvW6FhSAn3xxRc8/PDDWCwWXF1d2bJlCw8++KDRsUQKTYkoT35+ftcsspqcnIyzszN33333dR/j6uqKp6dnrpuIvXn00UdzzUbeunVrFixYYHQsKUEmTpzIc889h9Vq5c4772Tfvn00aNDA6FgihapElKfGjRuzatWqXNtWrlxJgwYNKFWqlEGpRArG/fffz7Zt23Bzc8NisdC1a1dmzZpldCwpAUaNGsXIkSMBKFOmDH/99Rf33HOPwalECp9Dlqe0tDSioqJsYzwSEhKIiooiMTERyPmFfvrpp2379+/fn8OHDzN8+HBiYmL48ssv+eKLLxgxYoQR8UUKXN26ddm/fz8eHh5YrVaeeeYZpkyZYnQsKcb69+/PO++8A4C/vz8JCQn4+fkZnEqkaDhkedq+fTvh4eGEh4cDMHz4cMLDw3nttdcAOHHihK1IAQQGBrJs2TLWrl1L3bp1eeONN/joo4/o2rWrIflFCkNgYCDx8fHcddddAAwbNsz2OyFSkLp3786nn34KQLVq1Th48KCGNkiJ4vDzPBUVzfMkjuL8+fOEhITYxvkNHjyYqVOnGpxKigOLxULbtm1twyDq1KnD9u3bcXZ2NjiZyI1pnicR+Vfe3t4cPHiQwMBAAD7++GN69eplcCpxdBaLhaZNm9qKU4sWLdi5c6eKk5RIKk8ixZC7uzsHDhwgNDQUgDlz5vDII48YnEocVVZWFuHh4WzZsgWAjh07sm7dOq1TJyWW/uWLFFPOzs7s2bPHdtn4ggULaNOmjcGpxNFkZmYSGhrKnj17gJwFfpcsWWJwKhFjqTyJFGNms5mtW7fSsmVLIGeKjqZNm2KxWAxOJo4gPT2dGjVqcODAAQD69u2rBX5FUHkSKfbMZjNr166lffv2AGzatIn69eurQMlNpaamUq1aNdu6icOHD+ezzz4zOJWIfVB5Eikhli5dymOPPQZAVFQUtWrVIisry+BUYo/Onj1LtWrVOHHiBABjx45l0qRJBqcSsR8qTyIlyLx58+jTpw8A+/fvp2bNmmRmZhqcSuxJcnIy1apV49SpU0DO8ivjxo0zNpSInVF5EilhvvzySwYOHAhAfHw8NWrUICMjw+BUYg+OHz9O9erVOXfuHAAfffQRL7/8ssGpROyPypNICTRt2jSGDx8OwOHDh6lWrRppaWkGpxIjHT58mKCgIFJTUwH4/PPPeeGFFwxOJWKfVJ5ESqhJkybx6quvAnDs2DGqVatm++CUkuXgwYPUrFmTtLQ0TCYTX3/9Nc8995zRsUTslsqTSAn2xhtv8OabbwJw8uRJqlWrxvnz540NJUUqLi6OWrVqkZ6ejslkYt68eTz11FNGxxKxaypPIiXcmDFjmDhxIgCnTp2iWrVqnD171uBUUhT27dtH3bp1uXTpEmazmYULF9K9e3ejY4nYPZUnEeHll19m8uTJAJw5c4Zq1apx+vRpg1NJYYqOjqZ+/fpkZGRgNptZsmQJDz30kNGxRByCypOIADB06FCmTp0KwLlz56hevTpJSUkGp5LCEBUVRf369bl8+TJms5lff/3VNomqiPw7lScRsRk8eDAzZswA4Pz58wQFBalAFTPbt2/n3nvvJTMzEycnJ1atWsWDDz5odCwRh6LyJCK59OvXjy+//BLIWaIjKCiI48ePG5xKCsK2bdto0qQJV65cwcnJidWrV3P//fcbHUvE4ag8icg1+vTpw+zZs4GcAhUcHKwC5eD+WZycnZ1Zt24dLVq0MDqWiENSeRKR6+rVqxdfffUVJpOJCxcuEBwczNGjR42OJbdg69atNGnShKysLJydnVm7di1NmzY1OpaIw1J5EpEbevrpp5kzZ46tQIWEhKhAOZitW7fSrFkzW3Fav369ipPIbVJ5EpGbeuqpp2wFKi0tTQXKgVyvODVu3NjoWCIOT+VJRP7VU089xddff20rUDVr1tQYKDu3bds2FSeRQqLyJCJ58uSTT+Y6hadB5PZr+/btucY4qTiJFCyVJxHJs6eeeorZs2fnGgOlAmVf/rc4rV27VsVJpICpPIlIvjz99NPMmjULyJnGICQkRBNp2omdO3fmmo5AV9WJFA6VJxHJt169euWaSDM4OJjk5GSDU5VsUVFRNGrUyDYB5po1a1ScRAqJypOI3JI+ffrYClRKSgpBQUFaTNgg0dHRNGzYMNfM4c2aNTM6lkixpfIkIresT58+fPrpp8D/rYV3/vx5Y0OVMDExMTRo0MC2Vt3KlSs1c7hIIVN5EpHb8vzzzzNt2jQAzp49S/Xq1UlNTTU4VckQHx9PvXr1uHz5MmazmeXLl2utOpEioPIkIrdt4MCBTJkyBYDTp09TvXp10tLSjA1VzCUkJFCnTh0yMjIwm8388ssvPPjgg0bHEikRVJ5EpEC8+OKLvP/++wAkJydTo0YN0tPTDU5VPCUmJhIWFsalS5cwmUwsWLCAdu3aGR1LpMRQeRKRAvPSSy8xYcIEAE6cOEFQUBAZGRkGpypejh8/TmhoKOnp6ZhMJubPn89DDz1kdCyREsVhy9P06dMJDAzEzc2N+vXr88cff9xw37Vr12Iyma65xcbGFmFikZLhP//5D+PHjwfg6NGjBAcHk5mZaXCq4iE5OZmQkBDS0tIwmUzMmzePhx9+2OhYIiWOQ5an77//nqFDhzJmzBh27dpF8+bNadeuHYmJiTd9XFxcHCdOnLDdqlevXkSJRUqW1157jdGjRwNw+PBhQkNDycrKMjiVYzt9+jRBQUG2wfhz5syhe/fuBqcSKZkcsjx98MEHPPvsszz33HOEhIQwZcoUAgIC+OSTT276OF9fX/z8/Gw3JyenIkosUvK89dZbjBgxAoC//vqLWrVqYbFYDE7lmM6fP09wcLBtGoiZM2fy1FNPGRtKpARzuPKUmZnJjh07iIyMzLU9MjKSTZs23fSx4eHhlC9fntatW7NmzZqb7nv58mVSU1Nz3UQkf9577z0GDx4MQGxsLHXr1lWByqe0tDSCgoI4c+YMANOmTePZZ581OJVIyeZw5en06dNkZ2dTrly5XNvLlSt3w/W1ypcvz2effcb8+fP5+eefCQoKonXr1qxfv/6GrzNhwgS8vLxst4CAgAL9PkRKiqlTp/Lcc88BsHfvXho2bKgClUcZGRkEBQXZlr6ZPHkyAwcONDiViDgbHeBWmUymXF9brdZrtl0VFBREUFCQ7evGjRtz5MgR3n///RvOxDtq1CiGDx9u+zo1NVUFSuQWff7551y6dIlvv/2W7du306pVq5v+8SI5R9mDg4M5fvw4AG+//TZDhw41NpSIAA545MnHxwcnJ6drjjIlJydfczTqZho1akR8fPwN73d1dcXT0zPXTURu3TfffEPXrl0B+OOPP6459S7/Jysri7CwMA4fPgzAq6++yqhRowxOJSJXOVx5cnFxoX79+qxatSrX9lWrVtGkSZM8P8+uXbsoX758QccTkZv46aefaN++PZDzO9ulSxdjA9khi8VC3bp1bX/cjRgxgjfeeMPgVCLyTw5XngCGDx/OzJkz+fLLL4mJiWHYsGEkJibSv39/IOeU29NPP23bf8qUKSxcuJD4+Hj27dvHqFGjmD9/vm0gq4gUnaVLl3LfffcBsGjRInr06GFwIvthsViIiIhg3759QM6yN++9957BqUTkfznkmKfHHnuMM2fO8Prrr3PixAnCwsJYtmwZlStXBnJmNv7nnE+ZmZmMGDGCY8eOcccddxAaGsrSpUttfwGLSNH67bffaNasGZs3b2bevHnceeedzJw50+hYhmvRogU7d+4EoHfv3rYFl0XEvpisVqvV6BCOIDU1FS8vL1JSUjT+SaQAWCwW6tevT1RUFABDhgzhww8/NDaUgSIjI23DEbp168YPP/xgcCKR4qEwPr8d8rSdiDg+s9nMjh07CA4OBuCjjz6yzUpe0jz88MO24tShQwcVJxE7p/IkIoYxm83s3buXwMBAIGd+tbfeesvgVEXrySefZOHChQC0atWKX375xdhAIvKvVJ5ExFDOzs7s37+fChUqADmX5U+ZMsXYUEWkf//+zJ07F4B7772X33//3eBEIpIXKk8iYjg3NzdiY2Px9fUFYNiwYcV+APlLL73Ep59+CkCtWrXYvHkzZrPekkUcgX5TRcQueHh4EBcXx1133QXA888/z/fff29wqsIxbtw4PvjgAwCqV6/Ozp07VZxEHIh+W0XEbnh7e7N//35Kly6N1WqlR48exW4M0KRJkxg/fjwAAQEBREdH4+zskLPGiJRYKk8iYlf8/PyIjo7mjjvuwGq10rlzZ9auXWt0rALx2WefMWLECCDn+4yNjcXFxcXgVCKSXypPImJ3KlWqxK5du3B1dcVisfDggw+ydetWo2Pdlu+++45+/foBUKZMGWJiYnB3dzc4lYjcCpUnEbFLQUFBbNmyhVKlSpGVlUWLFi1sy5Y4msWLF/Pkk08C4OnpSUxMDN7e3saGEpFbpvIkInarbt26rFu3DicnJzIzM2nQoAEJCQlGx8qXtWvX8vDDD2O1WnF3d2ffvn22qwpFxDGpPImIXWvcuDHLly/HbDaTkZFBrVq1OH78uNGx8mTbtm08+OCDWCwWXF1diYqKomLFikbHEpHbpPIkInbvwQcf5KeffsJkMnHx4kVCQ0M5e/as0bFuat++fTRr1oysrCxKlSrFn3/+SfXq1Y2OJSIFQOVJRBzCww8/zOzZswE4f/48wcHBpKWlGRvqBhISEmjQoAGZmZk4OTmxZs0aateubXQsESkgKk8i4jCefvpppk6dCsCpU6cICQkhMzPT4FS5JSUlUbt2bTIyMjCbzSxfvpymTZsaHUtECpDKk4g4lMGDB/Pmm28CcPToUcLCwsjKyjI4VY7z589Ts2ZN0tLSMJlM/PDDDzz44INGxxKRAqbyJCIOZ8yYMbz88ssAxMfH06BBAywWi6GZ0tPTCQ4O5ty5cwDMmjWLrl27GppJRAqHypOIOKSJEyfSt29fAHbv3k2rVq0My5KZmUnNmjU5efIkAJMnT6ZXr16G5RGRwnXL5enq2kwiIkb57LPP6NatGwB//PEHHTp0KPIMFouFOnXqcPjwYSBn0d+hQ4cWeQ4RKTq3XJ4WLVpk+//PPvtsgYQREcmvH374gTZt2gCwbNkynnjiiSJ7bYvFwr333ktsbCwAw4YNY+zYsUX2+iJijAI5bbdr166CeBoRkVvy66+/0rhxYyBnDblBgwYVyes+8MAD7NixA4DevXvzwQcfFMnrioixbrk8nTp1iiVLlnDo0KECjCMicms2bNhArVq1AJg+fTpjxowp1Nd75JFHWLNmDZAzB9WsWbMK9fVExH6YrFar9VYe+MEHHxAdHU10dDQHDhwgNDSUkJAQ2619+/YFndVQqampeHl5kZKSgqenp9FxROQ6srKyCAoK4u+//wZyBpVfvSqvIPXp08c2Yef999/P77//XuCvISIFozA+v2+5PP2vv//+21am9u/fzzfffFMQT2s3VJ5EHENGRgZVq1blxIkTAHz66ac8//zzBfb8w4cPZ/LkyQA0aNCArVu3YjbrwmURe2VoeerZsyeffvop7u7uBfLCjkblScRxpKamUrVqVc6cOYPJZGLevHl07979tp/39ddftw0IDw4OZt++fSpOInauMD6/8/xbP3fu3FzrSPXr1882GdxVV65cKZBQIiK3w9PTk/3791O6dGmsViuPP/44K1asuK3nnDp1qq04Va5cmd27d6s4iZRQef7N/98DVN99912u8nTy5ElKly5dcMlERG6Dr68v0dHR3HHHHVitVjp06MDmzZtv6bnmzJnDkCFDbM+7f/9+XFxcCjKuiDiQW/6z6Xpn++xtgU4RKdkqVarEtm3bcHFxITs7m5YtWxIdHZ2v51i8eDG9e/cGwNvbm5iYmBI7fEFEchToMWeTyVSQTycicttCQ0PZsGEDTk5OXLlyhYiICBISEvL02LVr1/Lwww9jtVq588472bdvH2XKlCnkxCJi7/JVnubOncvOnTttY5tUlkTEEURERLBixQrMZjMZGRnUrl2bpKSkmz5m+/btPPjgg1gsFlxdXdm9ezf+/v5FlFhE7Fmer7Zr0aIFu3fv5sKFC5QqVYqsrCy6d+9Os2bNqFevHmXLliUoKIjs7OzCzmwIXW0n4vjmz59Pt27dsFqt3HXXXfz99994e3tfs19cXBy1a9cmMzMTZ2dntm3bRt26dYs8r4jcPruY5yk+Pp4dO3awc+dOduzYwa5duzh//rztKJTKk4jYsy+++ILnnnsOAD8/Pw4ePJhrDNPRo0cJCgoiPT0ds9nMmjVraNGihVFxReQ2GTpVwVXVq1fn8ccfZ+LEifz++++cPXuWgwcPMm/ePEaOHFkgofJi+vTpBAYG4ubmRv369fnjjz9uuv+6deuoX78+bm5uVK1alRkzZhRRUhGxJ88++yyTJk0CICkpibCwMLKysgA4e/YsYWFhpKenYzKZWLRokYqTiFyjQAaMBwYG0q1bN95+++2CeLp/9f333zN06FDGjBnDrl27aN68Oe3atSMxMfG6+yckJNC+fXuaN2/Orl27GD16NEOGDGH+/PlFkldE7Mvw4cMZPXo0kPP+EB4eTlpaGsHBwaSkpADw9ddf07FjRyNjioidKrDlWYpSw4YNqVevHp988oltW0hICF26dGHChAnX7D9y5EgWL15MTEyMbVv//v3ZvXt3nud90Wk7keJn0KBBTJ8+HYBSpUrZLoaZNm0aAwcONDKaiBQQuzhtZ7TMzEx27NhBZGRkru2RkZFs2rTpuo/ZvHnzNfu3adOG7du333BW9MuXL5OamprrJiLFy7Rp03j88ceB/1sh4Y033lBxEpGbuq3ytHPnziKfGPP06dNkZ2dTrly5XNvLlSt3w0uPk5KSrrt/VlYWp0+fvu5jJkyYgJeXl+0WEBBQMN+AiNgNi8VCfHx8rm379+83KI2IOIrbKk8REREcOnSogKLkz//OMWW1Wm8679T19r/e9qtGjRpFSkqK7XbkyJHbTCwi9uaBBx5gx44dAPj4+AA5S08NGjTIyFgiYuduqzwZMVzKx8cHJyena44yJScnX3N06So/P7/r7u/s7Mzdd9993ce4urri6emZ6yYixccjjzzCmjVrAHj44Yc5efIktWrVAnKu5v3vf/9rZDwRsWMON+bJxcWF+vXrs2rVqlzbV61aRZMmTa77mMaNG1+z/8qVK2nQoAGlSpUqtKwiYp+eeeYZFixYAECrVq34+eefMZvN7Ny5k8DAQADefPNN25QGIiL/5HDlCXIuM545cyZffvklMTExDBs2jMTERPr37w/knHJ7+umnbfv379+fw4cPM3z4cGJiYvjyyy/54osvGDFihFHfgogY5OWXX2bWrFkA1KtXj99//912n7OzM/v378fPzw+AESNG8MUXXxiSU0Tsl7PRAW7FY489xpkzZ3j99dc5ceIEYWFhLFu2jMqVKwNw4sSJXHM+BQYGsmzZMoYNG8a0adPw9/fno48+omvXrkZ9CyJigLfeeov3338fyJnwd9u2bZjNuf+GdHNzIyYmhqpVq3Lu3Dn69u2Lt7e33i9ExOa25nkym83ExsZSo0aNgsxklzTPk4hjmz59um0geMWKFTl48CAuLi433D8pKYnq1auTlpaG2Wxm5cqVtG7duqjiikgB0TxPIiK34Ntvv7UVp7JlyxITE3PT4gQ5F5rs2bMHNzc3LBYLbdq0YevWrUURV0TsnMqTiBRry5Yto2fPngB4enqyf/9+PDw88vTYwMBAtm3bRqlSpcjOzqZFixa5VioQkZLptsrT2LFjbXOjiIjYmw0bNtCpUyesVivu7u7s27cv3+9ZYWFhrFu3DicnJzIzM6lfvz6HDx8upMQi4ggccm07I2jMk4hjiYqKIiIigqysLFxdXdm9ezdBQUG3/HwrVqygXbt2WK1WSpcuzV9//YWvr28BJhaRwqAxTyIieRAfH0+jRo3IysrC2dmZDRs23FZxgpz1MOfNm4fJZOLChQuEhIRozUuREkrlSUSKlaNHjxIeHs7ly5cxm82sWrWKBg0aFMhzd+/enRkzZgBw9uxZgoODycjIKJDnFhHHkefy1LNnT9LT0wszi4jIbTl79ixhYWFcvHgRk8nEggULaNWqVYG+xvPPP8/EiROBnDnlQkNDycrKKtDXEBH7lufyNHfuXNLS0mxf9+vXj3PnzuXa58qVKwWXTEQkH9LS0ggODiYlJQWAOXPm8NBDDxXKa7388suMGjUKgL///pt69ephsVgK5bVExP7kuTz977jy7777Lld5OnnyJKVLly64ZCIieZSZmUlISAinTp0CYOrUqTz11FOF+ppvv/02AwYMAGDv3r00bdq0UF9PROzHLY95ut5FepmZmbcVRkQkv7KysggLC+Po0aMAvPHGGwwePLhIXnv69On06NEDgC1bttCmTZsieV0RMVaBDhg3mUwF+XQiIjdlsViIiIggPj4eyFnI99VXXy3SDHPnzqV9+/YArFy5ku7duxfp64tI0ctXeZo7dy47d+60jW1SWRIRI913331ERUUB8Nxzz/Hee+8ZkmPp0qU0a9YMgB9//JHnn3/ekBwiUjTyPElmixYt2L17NxcuXKBUqVJkZWXRvXt3mjVrRr169ShbtixBQUFkZ2cXdmZDaJJMEfvSsWNHli5dCsCjjz7Kjz/+aGgei8VCvXr12L17N5AzqPzqVXkiYpzC+PzO9wzj8fHx7Nixg507d7Jjxw527drF+fPnbUehVJ5EpLA9+eSTzJ07F4DIyEhWrFhhcKIcWVlZ1KxZ03Ya8c0332TMmDEGpxIp2Qrj89s5vw+oXr061atX5/HHH7dtS0hIYPv27ezatatAQomI3MiAAQNsxalRo0Z2U5wAnJ2diY6O5p577uHo0aO8+uqreHp68sILLxgdTUQKkNa2yyMdeRIx3siRI22nwmrVqkVUVBRms/0tlJCWlsY999xDcnIyALNnz6ZXr14GpxIpmbS2nYiUWBMmTLAVp2rVqrFz5067LE4AHh4exMTE4O3tDUCfPn1YsGCBsaFEpMDY5zuPiMg/fPzxx4wePRqAChUqsHfvXpyd8z3qoEiVKVOGffv2ceedd2K1Wnn00UdZtWqV0bFEpACoPImIXZszZ45tzFDZsmWJjY3Fzc3N4FR54+/vz969e3Fzc8NisdCuXTs2b95sdCwRuU0qTyJit+bPn0/v3r0B8PLyIjY2Fg8PD2ND5VNgYCDbtm3DxcWF7OxsWrZsaZubSkQck8qTiNilFStW0L17d6xWK3feeSf79++nTJkyRse6JWFhYWzYsAFnZ2euXLlCo0aNiIuLMzqWiNwilScRsTsbNmygQ4cOWCwW3Nzc2Lt3L/7+/kbHui0RERGsWrUKs9nM5cuXCQ8PJzEx0ehYInILVJ5ExK7s3LmT++67j+zsbFxcXNi+fTuBgYFGxyoQrVq1YsmSJZhMJi5dukRYWBhJSUlGxxKRfFJ5EhG7ERMTQ+PGjcnKysLZ2ZmNGzcSGhpqdKwC1b59e+bNm4fJZOLChQvUrFmT8+fPGx1LRPJB5UlE7EJCQgL169cnMzMTJycn1qxZQ4MGDYyOVSi6d+/O559/DsC5c+cICgoiLS3N4FQiklcqTyJiuKNHj1KrVi0uXbqEyWRiyZIlNGvWzOhYherZZ59lypQpACQnJxMcHExGRoaxoUQkT1SeRMRQp0+fJjQ0lIsXL2Iymfjxxx9p166d0bGKxIsvvsibb74JwLFjx6hZsyZZWVkGpxKRf6PyJCKGOX/+PEFBQaSmpgI5E2J27drV4FRFa8yYMbbZ0xMSEqhVqxYWi8XgVCJyMypPImKI9PR0goODOXv2LAAzZszgqaeeMjiVMd566y2GDBkCQGxsLPXq1VOBErFjKk8iUuQyMjIICgri5MmTAEyaNIl+/foZnMpYH374Ic899xwAu3fvpkmTJipQInZK5UlEilRmZiY1a9bk6NGjAIwfP57hw4cbnMo+fP755/To0QOArVu30rp1a4MTicj1OFx5OnfuHD179sTLywsvLy969uz5r3Ok9O7dG5PJlOvWqFGjogksIjZZWVnUrl2bhIQEAEaNGsVrr71mcCr7MnfuXDp37gzA2rVrad++vcGJROR/OVx5euKJJ4iKiuLXX3/l119/JSoqip49e/7r49q2bcuJEydst2XLlhVBWhG5ymKxUK9ePduabkOHDuXtt982OJV9WrhwIZGRkQAsX768xA2iF7F3DlWeYmJi+PXXX5k5cyaNGzemcePGfP755/zyyy//usimq6srfn5+tpujLjAq4ogsFgsNGzZk7969APTr14/JkycbnMq+rVixghYtWgDw888/8+STTxqcSESucqjytHnzZry8vGjYsKFtW6NGjfDy8mLTpk03fezatWvx9fWlRo0a9O3bl+Tk5Jvuf/nyZVJTU3PdROTWtGzZku3btwPQs2dPZsyYYXAix7BmzRoiIiKAnNN5VweUi4ixHKo8JSUl4evre812X1/fmy6u2a5dO7799ltWr17NpEmT2LZtG/fffz+XL1++4WMmTJhgG1fl5eVFQEBAgXwPIiVN69at2bBhAwBdu3Zlzpw5BidyHGazmS1btlC7dm0AvvjiCwYNGmRwKhGxi/I0bty4awZ0/+/t6l+tJpPpmsdbrdbrbr/qscceo0OHDoSFhdGpUyeWL1/OgQMHWLp06Q0fM2rUKFJSUmy3I0eO3P43KlLCtGvXjtWrVwPQsWNHfvrpJ4MTOR6z2cyuXbsICQkBYPr06bo6UcRgzkYHABg8eDCPP/74TfepUqUKe/bssc0L80+nTp2iXLlyeX698uXLU7lyZeLj42+4j6urK66urnl+ThHJrUuXLvz6668AREZGsmTJEoMTOS6z2cyePXuoWbMm8fHxTJ48GTc3Nw24FzGIXZQnHx8ffHx8/nW/xo0bk5KSwp9//sm9994L5MyFkpKSQpMmTfL8emfOnOHIkSOUL1/+ljOLyI11796dRYsWATnjnVasWGFwIsfn7OxMdHQ0wcHBJCQkMGHCBFxdXRk7dqzR0URKHLs4bZdXISEhtG3blr59+7Jlyxa2bNlC37596dixI0FBQbb9goODWbBgAQBpaWmMGDGCzZs3c+jQIdauXUunTp3w8fHh4YcfNupbESm2nn76aX788Ucg5w+eq6ft5Pa5uLiwf/9+2xjMcePGMWHCBINTiZQ8DlWeAL799ltq1apFZGQkkZGR1K5dm6+//jrXPnFxcaSkpADg5OTE3r176dy5MzVq1KBXr17UqFGDzZs3U7p0aSO+BZFiq0+fPrbfx/r167NhwwbMZod7m7Frbm5uxMbG2o6cjx49mvfee8/gVCIli8lqtVqNDuEIUlNT8fLyIiUlBU9PT6PjiNid559/ns8//xyAOnXqsHPnThWnQpSWlka1atVs40AnT57M0KFDjQ0lYocK4/Nb72wictsGDRpkK05hYWEqTkXAw8ODAwcOULZsWQCGDRvGxx9/bHAqkZJB724icltefPFFpk+fDuSMN9y1a5eKUxHx9PTkwIED3H333QC88MILfPrppwanEin+9A4nIrfspZde4qOPPgKgevXq7N27F2dnu7iIt8Tw9vbmwIED3HXXXQD079+fmTNnGpxKpHhTeRKRW/Lyyy/zwQcfAFC1alWio6NVnAxSpkwZDhw4gLe3NwB9+/bliy++MDaUSDGm8iQi+TZy5Ejef/99AAIDA4mJicHFxcXgVCWbj48PcXFxeHl5AfDcc88xa9Ysg1OJFE8qTyKSL6NGjWLixIlAzsz/sbGxKk52wtfXl9jYWNsVRc888wxfffWVwalEih+VJxHJs9GjR/POO+8AULlyZeLi4lSc7Iyfnx9xcXG2AtW7d28txixSwFSeRCRPRo8ebZvNulKlSjriZMf8/PyIiYmxTQSsAiVSsFSeRORf/bM4BQQEEBcXh5ubm8Gp5Gb8/f2JjY2ldOnSWK1WFSiRAqTyJCI3NWrUqFxHnA4cOKDi5CD+t0D16tVLY6BECoDKk4jc0MiRI68Z46Ti5FiuFqh/joHSVXgit0flSUSu65VXXrFdVVe5cmViY2NVnByUv78/MTExua7CU4ESuXUqTyJyjZdeeon33nsPyJmOQKfqHJ+/v3+ueaCeeeYZPvvsM4NTiTgmlScRyWXIkCG2mcMDAwM1HUEx4ufnR2xsrK1A9evXz7YuoYjkncqTiNgMGDCAqVOnAlCtWjVNR1AM+fn55VoLb9CgQXz44YcGpxJxLCpPIgLkrIc2Y8YMAIKCgrTkSjHm6+vLgQMHKFOmDABDhw5l0qRJBqcScRwqTyJC7969mTlzJgA1a9bUIr8lgI+PDwcPHsTHxweAESNG2K6sFJGbU3kSKeF69Ohhm/unVq1a7N27V8WphPD29ubgwYP4+voCOXN6jR8/3uBUIvZP5UmkBOvSpQvz5s0DoF69ekRFRWE2622hJPH09OTgwYOUL18egHHjxjFq1CiDU4nYN71LipRQbdu2ZdGiRQA0atSIbdu2qTiVUB4eHvz1118EBAQA8M477/Diiy8anErEfumdUqSEsVgstGrVihUrVgDQqlUrNm7cqOJUwrm7u3PgwAGqVq0KwEcffUT//v0NTiVin/RuKVKCWCwWmjZtyrp16wBo06YNa9asUXESANzc3IiJiSEoKAiATz/9lF69ehmcSsT+6B1TpITIysqiXr16bNmyBcgZ7/Trr78anErsjYuLC9HR0dSqVQuAOXPm0LVrV4NTidgXlSeREiAzM5PQ0FB2794N5Fxht2DBAoNTib1ydnYmKiqK+vXrA/Dzzz/Tpk0bg1OJ2A+VJ5FiLj09nRo1anDgwAEgZzLMuXPnGpxK7J3ZbObPP/+kZcuWAKxcuZJmzZphsVgMTiZiPJUnkWIsNTWVatWqcfjwYQCGDx+uxWAlz8xmM2vXrqV9+/YAbNy4kQYNGqhASYmn8iRSTCUnJ1O1alVOnDgBwNixY7UEh9ySpUuX8thjjwGwa9cuQkNDyczMNDiViHFUnkSKocTERKpVq8aZM2cAmDhxIuPGjTM2lDi0efPm8eyzzwIQGxtLtWrVSE9PNziViDFUnkSKmZiYGIKDg7lw4QImk4kZM2bw8ssvGx1LioGZM2cybNgwAI4cOUJgYCBnz541OJVI0VN5EilGtm3bRt26dbl06RImk4l58+bRr18/o2NJMfLBBx/wxhtvAP93avjo0aMGpxIpWipPIsXE77//TuPGjcnMzMTJyYnly5fTvXt3o2NJMfTqq68ydepUAFJSUggKCiIuLs7gVCJFx+HK01tvvUWTJk1wd3fH29s7T4+xWq2MGzcOf39/7rjjDlq1asW+ffsKN6hIEfrhhx+IjIwkOzubUqVK8ccff2heHilUgwcP5ptvvsFkMpGenk7t2rXZtm2b0bFEioTDlafMzEy6devGgAED8vyYiRMn8sEHH/Dxxx+zbds2/Pz8ePDBB7lw4UIhJhUpGtOnT+exxx7DYrHg5ubGzp07ady4sdGxpAR48sknWbx4MWazmczMTBo3bmxbM1GkOHO48jR+/HiGDRtmWzrg31itVqZMmcKYMWN45JFHCAsL46uvviI9PV0TBYrDGz9+PIMGDQLA09OT2NhYwsLCDE4lJUnHjh1Zv349pUqVIjs7m3bt2vHdd98ZHUukUDlcecqvhIQEkpKSiIyMtG1zdXWlZcuWbNq06YaPu3z5MqmpqbluIvZk4MCBtukHfH19OXjwIJUrVzY2lJRITZs2ZdeuXdxxxx1YrVaeeOIJ25gokeKo2JenpKQkAMqVK5dre7ly5Wz3Xc+ECRPw8vKy3QICAgo1p0h+dOvWjU8++QSAypUrk5CQgI+Pj8GppCQLDQ0lNjYWLy8vAIYMGcKYMWMMTiVSOOyiPI0bNw6TyXTT2/bt22/rNUwmU66vrVbrNdv+adSoUaSkpNhuR44cua3XFykIFouF5s2b89NPPwFQq1Yt/vrrL9zd3Q1OJgKVKlXi77//tv2x+vbbb9OnTx+DU4kUPGejA0DOVRuPP/74TfepUqXKLT23n58fkHMEqnz58rbtycnJ1xyN+idXV1dcXV1v6TVFCkNGRgZ169a1XRLeqlUrfv/9d8xmu/gbSASAMmXK8Pfff1O7dm0OHjzI7NmzOX78OMuXL9e/VSk27KI8+fj4FNoph8DAQPz8/Fi1ahXh4eFAzhV769at49133y2U1xQpaGfPniUsLMy2Tt0TTzzBt99+a3Aqketzd3fnwIEDNGnShK1bt7Jy5Urq16/Ptm3bcHa2i48dkdvicH8GJCYmEhUVRWJiItnZ2URFRREVFUVaWpptn+DgYBYsWADknK4bOnQob7/9NgsWLCA6OprevXvj7u7OE088YdS3IZJnhw8fzrXA78svv6ziJHbPbDazZcsWOnfuDEBUVBTVqlXL9V4t4qgc7k+A1157ja+++sr29dWjSWvWrKFVq1YAxMXFkZKSYtvnlVde4dKlSwwcOJBz587RsGFDVq5cSenSpYs0u0h+bd++nebNm5ORkQHAlClTePHFFw1OJZJ3CxcuZODAgXzyySccPnyYypUrs3fvXvz9/Y2OJnLLTFar1Wp0CEeQmpqKl5cXKSkpeHp6Gh1HSoBFixbRtWtXsrOzMZvNfP/99zz66KNGxxK5JRMmTGD06NEA3HHHHWzatIm6desaG0pKhML4/Ha403YiJcHUqVN5+OGHyc7OxsXFhfXr16s4iUMbNWoUX3/9NSaTiUuXLtGgQQOWLVtmdCyRW6LyJGJnhg8fzpAhQ7BarZQuXZr9+/fTtGlTo2OJ3LannnqK1atX22Yj79ixI59++qnRsUTyTeVJxI488sgjTJ48GYDy5ctz6NAh7rnnHoNTiRScVq1asXv3bu68806sViv9+/dn5MiRRscSyReVJxE7kJmZSd26dW1XidaqVYtDhw5RpkwZg5OJFLyQkBAOHTpkm2tv4sSJPPzwwwanEsk7lScRg50+fZoqVaqwe/duANq1a0dUVBQuLi4GJxMpPD4+Phw6dIjQ0FAg56q8unXrkpmZaXAykX+n8iRioOjoaKpUqWKbw2no0KEsW7ZMMzFLieDm5saePXto164dALt376Zy5cqcPn3a4GQiN6d3aBGDLFu2jPDwcC5evIjJZGLatGm28U4iJYXZbGbZsmUMHToUyFlKq0qVKkRHRxsbTOQmVJ5EDDBp0iQ6duxIVlYWzs7OLF++nIEDBxodS8QwkydPZtq0aZhMJi5evEjdunVZtGiR0bFErkvlSaSI9enThxEjRmC1WvHw8CAqKoo2bdoYHUvEcAMHDmT58uU4OzuTnZ1Nly5dmDBhgtGxRK6h8iRSRLKysmjUqBGzZ88GoFKlShw5csQ2YFZEoE2bNkRHR9tmgh49erTWIRW7o/IkUgSuXlG3detWAFq0aEFCQgLe3t7GBhOxQ0FBQRw5coSqVasC8N1331G/fn1diSd2Q+VJpJBt376dypUrc+zYMQD69+/PunXrdEWdyE14enoSHx/PAw88AMDOnTsJCAjg6NGjBicTUXkSKVRfffUVDRs2JD09HZPJxNSpU/nkk0+MjiXiEMxmM6tWrWLIkCEAJCcnU61aNdavX29wMinpVJ5ECsmQIUPo3bs3FosFV1dX1q5dy+DBg42OJeJwPvzwQ7788kvMZjOXL1+mVatWTJ061ehYUoKpPIkUsKysLJo3b257c/f19eWvv/6iRYsWBicTcVx9+vRh69atuLu7Y7VabX+ciBhB5UmkAB0/fpzKlSuzYcMGABo0aMCRI0eoWLGiwclEHF+DBg04fPiw7ffpq6++Ijw8nIyMDIOTSUmj8iRSQFavXk3VqlU5fvw4AM8++yzbtm3TGnUiBcjHx4fDhw/TqlUrAKKiovD39yc+Pt7YYFKiqDyJFIB33nmHBx54gMuXL2M2m/n000+ZOXOm0bFEiiWz2cyaNWt4+eWXATh37hw1a9bkp59+MjiZlBQqTyK3wWKx0LlzZ0aNGoXVauXOO+9k69atPP/880ZHEyn2Jk6cyM8//4yzszNZWVl069aN4cOHGx1LSgCT1Wq1Gh3CEaSmpuLl5UVKSopt5lsp2U6fPk1ERASHDh0CoHr16vz555+a+FKkiB08eJB7772Xs2fPAtC4cWPWrl2rU+YCFM7nt448idyC1atXExAQYCtOXbt2JTY2VsVJxAD33HMPx44do379+gBs3ryZ8uXLExcXZ3AyKa5UnkTyafz48TzwwANkZGRgMpmYMmUKP/30k2YMFzGQm5sb27dv54UXXgDg7NmzhIaG8u233xqcTIojnbbLI522k6ysLCIjI1mzZg0AHh4erFmzhgYNGhicTET+af78+fTo0YMrV64AOVe+6gKOkkun7UQMEh8fj7+/v604hYaGcuzYMRUnETvUtWtX4uPj8fX1BeCLL74gJCSE8+fPGxtMig2VJ5F/MWvWLEJCQjh16hSQM9NxdHS0jkCK2LGri3Hfd999AMTGxuLv78/q1asNTibFgcqTyA1YLBa6devGM888Q3Z2Ns7OzsybN48vv/zS6GgikgfOzs6sXr2at99+G5PJxKVLl2jdujWjRo0yOpo4OI15yiONeSpZjh8/TqNGjThy5AgA/v7+bN68mUqVKhmcTERuxbZt27j//vtJS0sDcpZ6WbduHe7u7gYnk8KmMU8iReD777+nSpUqtuLUuXNnjhw5ouIk4sAiIiI4efIk9erVA2D79u2UK1eOjRs3GpxMHJHKk8j/Z7FY6N69O48//jhXrlzBycmJGTNmsHDhQk1DIFIMuLu7s2PHDtuyLmlpaTRv3pzRo0cbnEwcjU7b5ZFO2xVvhw8fpkmTJrZFfcuVK8cff/xB9erVDU4mIoVh/fr1tG/fnosXLwJQp04d1q9fr/f3Ykin7UQKwcyZM6lWrZqtOHXp0oXjx4+rOIkUYy1atCApKck2K/nu3bvx8/Nj1apVBicTR+Bw5emtt96iSZMmuLu753kpjN69e2MymXLdGjVqVLhBxe5lZmbSunVr+vbtS1ZWFs7Oznz11VcsWLBAp+lESgAPDw+2b9/O2LFjbVfjRUZG0qdPHywWi9HxxI453CdEZmYm3bp1Y8CAAfl6XNu2bTlx4oTttmzZskJKKI5g8+bNlC1b1jbnS+XKlTl48CBPP/20wclEpKiNGzeOrVu32v4gnz17NpUrVyYhIcHYYGK3HK48jR8/nmHDhlGrVq18Pc7V1RU/Pz/brUyZMoWUUOzdsGHDaNKkCampqQD079+fQ4cO6Wo6kRLs6tV4bdu2BeDo0aNUr16dqVOnGpxM7JHDladbtXbtWnx9falRowZ9+/YlOTn5pvtfvnyZ1NTUXDdxbImJidxzzz1MmTIF+L+16T755BNjg4mIXXBxcWH58uXMnj2bUqVKkZ2dzZAhQ2jUqJE+AySXElGe2rVrx7fffsvq1auZNGmSbbK0y5cv3/AxEyZMwMvLy3YLCAgowsRS0CZNmkTVqlX5+++/AWjWrBknT56kVatWxgYTEbvTq1cvEhMTbReNbN26lXLlyjF//nyDk4m9sIvyNG7cuGsGdP/vbfv27bf8/I899hgdOnQgLCyMTp06sXz5cg4cOMDSpUtv+JhRo0aRkpJiu12dMFEcy9mzZ6lXrx4jRoywLbEyY8YM/vjjD80sLCI35Ofnx4EDB3jllVcwmUxkZGTw6KOP0rFjRzIzM42OJwZzNjoAwODBg3n88cdvuk+VKlUK7PXKly9P5cqViY+Pv+E+rq6uuLq6FthrStGbM2cOzz//vO0IY0hICKtXr8bPz8/gZCLiKN5991169uzJAw88wMmTJ1m6dCm+vr4sXLhQR65LMLsoTz4+Pvj4+BTZ6505c4YjR45Qvnz5IntNKTqpqam0bduWzZs3A2A2m3nttdcYO3aswclExBGFhYVx/PhxnnnmGb766itSUlK47777ePTRR/nuu+9wdraLj1IpQnZx2i4/EhMTiYqKIjExkezsbKKiooiKirIt9ggQHBzMggULgJzp90eMGMHmzZs5dOgQa9eupVOnTvj4+PDwww8b9W1IIfnqq6/w9fW1FafKlSuzf/9+FScRuS1ms5nZs2ezbt0625QGP/30Ez4+PrYpT6TkcLjy9NprrxEeHs7YsWNJS0sjPDyc8PDwXGOi4uLiSElJAcDJyYm9e/fSuXNnatSoQa9evahRowabN2+mdOnSRn0bUsDOnz9P48aN6d27N5cvX8ZsNvPKK69w6NAhgoKCjI4nIsVEixYtOHXqlG2oSUpKCq1bt6Zr165kZWUZnE6Kita2yyOtbWe/pk6dyksvvcSVK1eAnKNNq1at0vIqIlKo1q9fT+fOnTl//jyQM/3JN998Q+fOnY0NJrlobTuRfzh8+DAhISEMGTKEK1euYDab+c9//sOhQ4dUnESk0F09CtWjRw8gZ5hIly5daNmypeaFKuZUnsThWCwWXnnlFapWrUpsbCyQcyXd33//zYQJEwxOJyIlibOzM3PnzuXPP/+0Xcm7fv16ypYty/Tp0w1OJ4VF5Ukcyvr16ylfvjzvvfceFosFFxcXPvroI/bv30/lypWNjiciJVRERATHjh1j+PDhmM1mMjMzGTRoEEFBQTedFkcck8qTOIS0tDQiIyNp2bKlbWmdpk2bcvLkSV544QWD04mI5FyRN2nSJP766y/bhSoHDhwgKCiI3r17a0B5MaLyJHbvgw8+4O6772bVqlUAeHt7s2TJEjZs2GC7ZFhExF4EBgYSGxvL9OnTcXNzw2q18tVXX1GmTBl++ukno+NJAVB5Eru1detWKlWqxEsvvURmZiYmk4n+/ftz5swZOnbsaHQ8EZGbGjBgQK73qwsXLtCtWzfq1KlDQkKCwenkdqg8id05e/YsrVu3plGjRrY1BUNDQzl48CCffPIJZrP+2YqIY3B3d2fJkiVs2bKFChUqALBnzx7uuecennzySa2T56D0KSR24+pVdOXKlbPN2Ovp6cm8efOIjo4mMDDQ4IQiIremYcOGHD16lHfffRdXV1esVitz587F29ubjz/+2Oh4kk8qT2IXro4HeO+998jKysJsNvPCCy9w7tw5HnvsMaPjiYgUiFdeeYWzZ8/yyCOPAHDp0iVeeOEFKlasyO+//25wOskrlScx1MaNG6lSpQq9e/e2LanTrFkzTpw4wUcffaRTdCJS7Li7uzN//nz2799PcHAwAMeOHeOBBx4gPDxcUxs4AH0yiSESEhJo1KgRzZo14/DhwwBUrVqVLVu28Mcff+Dr62twQhGRwhUSEkJMTAw///wzZcuWBSAqKooaNWrQsWNHzp49a3BCuRGVJylSycnJtGnThnvuuYetW7cCOVMPfPPNNxw8eJCGDRsanFBEpGg9/PDDJCcnM2HCBNzc3ABYunQpZcuW5emnnyYjI8PghPK/VJ6kSKSlpdG9e3fKly/PypUrsVqtuLq68t///pczZ87w5JNPGh1RRMRQ//nPf0hJSaFPnz44OTlhsVj4+uuv8fT0ZNiwYZpk046oPEmhSk9Pp0+fPtx11138+OOPWCwWnJ2d6d+/P6mpqbz++usa1yQi8v+5uLjw5Zdfcvr0abp06YLJZOLKlStMmTIFT09PRo8ejcViMTpmiadPLSkU6enp9OrVCy8vL2bPnm27gq5bt26cO3eOTz75BBcXF6NjiojYJW9vbxYsWMDRo0dp1aoVkHNl3oQJE/Dw8GDkyJE6EmUglScpUKmpqTz99NN4eXkxZ84csrKyMJlMdOzYkRMnTvDDDz/g4eFhdEwREYfg7+/PmjVrOHDgAE2aNAFyStTEiRMpXbo0L7/8skqUAVSepEAcP36cjh07ctddd/H111/bjjQ99NBDJCcns2TJEl1BJyJyi6pXr87GjRtzlaiMjAzef/997rzzTvr06UNaWprBKUsOlSe5Lfv27aNZs2ZUrFiRpUuXYrFYMJvNdO7cmZMnT7Jo0SJ8fHyMjikiUiz8s0Q1a9YMk8lEZmYms2fPxsvLi06dOnH8+HGjYxZ7Kk9yS+bPn09QUBBhYWFs3LgRq9WKi4sLvXv35ty5cyxcuFClSUSkkFSvXp0//viDxMREOnTogNlsxmKx8Msvv1ChQgUiIiLYuHGj0TGLLZUnybOsrCzGjRvH3XffzaOPPsqBAwcAuPPOO3nllVe4ePEis2bNwtPT0+CkIiIlQ8WKFfnll184d+4cvXr1olSpUgBs377ddlbgk08+0RV6BUzlSf5VfHw8nTp1wt3dnfHjx9tmva1QoQLTp08nNTWVd999F2dnZ4OTioiUTJ6ensyePZv09HTGjh3LXXfdBeQs+zJw4EA8PDzo1asXycnJBictHlSe5LosFguffvopVatWpUaNGvzyyy9cuXIFgAYNGrBhwwaOHj3KgAEDNE+TiIidcHZ2Zty4cZw9e5affvqJGjVqADlX6M2ZM4dy5cpRu3ZtFi1aZHBSx6ZPPcklOjqazp074+7uTv/+/UlISADgjjvuoGfPnpw4cYJt27bRtGlTg5OKiMjNdO3albi4OGJjY+nQoYPt7MDevXvp0qULnp6e9OrVi6NHjxqc1PGoPAlpaWmMHDkSPz8/atWqxeLFi7l8+TKQMyhx9uzZpKWlMWfOHPz8/AxOKyIi+REUFMQvv/xim2TT398fgAsXLjBnzhwCAgKoWrUq7733nuaMyiOT1Wq1Gh3CEaSmpuLl5UVKSkqxGBCdmZnJxx9/zOeff05cXBz//Gfg4eHBQw89xNtvv03lypUNTCkiIoUhKiqKMWPG8Ntvv5GZmWnbbjabCQ8PZ8iQITz11FPFYlhGYXx+qzzlUXEoT5mZmXzxxRd8+umn7N27N9fVF2azmYiICF599VU6duxoYEoRESkqFouFr776ikmTJrF///5cf0g7OzvTsGFDhgwZwqOPPuqwRUrlyUCOWp7Onz/Phx9+yLx58zhw4ECuwmQymQgKCqJv374MHjxYa82JiJRgqampvPPOO3z77bckJibmus/Z2ZnatWvTu3dv+vbti5ubm0Ep80/lyUCOVJ62bt3K9OnT+e23364702xgYCC9evXipZde0jpzIiJyjeTkZN566y1+/PFHTpw4kes+k8lEYGAgbdu2ZfDgwYSEhBiUMm9Ungxkz+Xp+PHjfPnllyxcuJDo6GjbYO+rzGYzoaGhPPXUUwwePBh3d3eDkoqIiKM5ffo0kyZN4scff+Tvv//mf2vDnXfeSXh4OI8++ii9evXC29vbmKA3oPJkIHsqT4mJicyePZulS5eyb98+Ll68eM0+Hh4eRERE0KtXL3r27Omw56pFRMR+ZGRkMGPGDObNm8fu3bvJyMi4Zh9vb2/q1q1L586deeqppwxfqkvlyUBGlaf09HSWLFnCokWL2LZtG4mJibmujLjK2dmZwMBAOnTowODBg7nnnnuKLKOIiJRMO3fuZPr06axcuZJjx45ddxmYO+64g8DAQBo3bswjjzxCZGRkka5IUeLL06FDh3jjjTdYvXo1SUlJ+Pv789RTTzFmzJibDna2Wq2MHz+ezz77jHPnztGwYUOmTZtGaGhonl+7sMuTxWJh//79rFixgj/++IN9+/Zx7NgxLl26dN39r5al+++/n169etG4ceMCzyQiIpJXFouFFStW8M0337Bx40aOHDlywzX1PDw8qFixInXq1KFly5ZERkYW2h/9Jb48/frrr3z//ff06NGDatWqER0dTd++fenZsyfvv//+DR/37rvv8tZbbzF79mxq1KjBm2++yfr164mLi6N06dJ5eu2C+OGnpaWxc+dOdu3aRUxMDAcOHCAhIYHk5GTS09Nv+tgyZcpQs2ZN7rvvPnr06GH3A/RERKRks1gsbNu2jR9++IH169dz4MABUlNTb7i/yWTCw8ODcuXKUbVqVYKCgggJCaFevXrUqVPnlq/wK/Hl6Xree+89PvnkE/7+++/r3m+1WvH392fo0KGMHDkSgMuXL1OuXDneffdd+vXrl6fXufrDX7RoEc7Ozly6dInMzEzS09M5c+YM586d49y5c5w/f55z585x+vRpzp07R2pqKunp6Vy+fDlPq1o7OztTtmxZgoKCuPfee2nXrh3NmjXTorsiIuLwMjIy+O2331i1ahXbtm3jr7/+4uzZs2RnZ//rY52cnHB1deXOO+/E09OTu+66i7Jly+Lt7Y23tzd33XUXd999N3fffTfh4eHUrl0bKJzy5PCfyCkpKZQpU+aG9yckJJCUlERkZKRtm6urKy1btmTTpk03LE+XL1/OddXa1bbcuXPn285cqlQpPDw88PHx4Z577qF27do0bdqUFi1a2N1VCiIiIgXFzc2Njh07XjMZc1JSEuvWrWPLli3s3buXhIQEzpw5w8WLF21LxmRnZ5Oenk56ejqnTp266es0atSIzZs3F9r34dDl6eDBg0ydOpVJkybdcJ+kpCQAypUrl2t7uXLlOHz48A0fN2HCBMaPH3/d+0wmk+1/zWYzTk5OlCpVCldXV1xdXXF3d7c1Yn9/fypWrEhgYCD169cnJCREV76JiIj8g5+fH4899hiPPfbYNfdlZmYSHR3Nrl27OHToEEeOHCEpKcl2hufSpUtkZGSQmZnJlStXyM7OvuYzv6DZRXkaN27cDYvKVdu2baNBgwa2r48fP07btm3p1q0bzz333L++xtXCc5XVar1m2z+NGjWK4cOH275OTU0lICDALqYqEBERKSlcXFyoV68e9erVMzqKjV2Up8GDB/P444/fdJ8qVarY/v/x48e57777aNy4MZ999tlNH+fn5wfkHIEqX768bXtycvJNm+nVo0giIiIi/2QX5cnHxyfPk2gdO3aM++67j/r16zNr1qx/PQUWGBiIn58fq1atIjw8HMg5BLhu3Trefffd284uIiIiJYtDDb45fvw4rVq1IiAggPfff59Tp06RlJRkG9d0VXBwMAsWLAByTtcNHTqUt99+mwULFhAdHU3v3r1xd3fniSeeMOLbEBEREQdmF0ee8mrlypX89ddf/PXXX1SsWDHXff+ccSEuLo6UlBTb16+88gqXLl1i4MCBtkkyV65cmec5nkRERESucvh5noqKPa1tJyIiInlTGJ/fDnXaTkRERMRoKk8iIiIi+aDyJCIiIpIPKk8iIiIi+aDyJCIiIpIPKk8iIiIi+aDyJCIiIpIPKk8iIiIi+aDyJCIiIpIPDrU8i5GuTsSemppqcBIRERHJq6uf2wW5oIrKUx6dOXMGgICAAIOTiIiISH6dOXMGLy+vAnkulac8KlOmDACJiYkF9sMXEfuQmppKQEAAR44c0dqVIsVMSkoKlSpVsn2OFwSVpzwym3OGh3l5eenNVaSY8vT01O+3SDF19XO8QJ6rwJ5JREREpARQeRIRERHJB5WnPHJ1dWXs2LG4uroaHUVECph+v0WKr8L4/TZZC/LaPREREZFiTkeeRERERPJB5UlEREQkH1SeRERERPJB5UlEREQkH1Se8mD69OkEBgbi5uZG/fr1+eOPP4yOJCIiIv8wYcIEIiIiKF26NL6+vnTp0oW4uLhCeS2Vp3/x/fffM3ToUMaMGcOuXbto3rw57dq1IzEx0ehoIiIi8v+tW7eOQYMGsWXLFlatWkVWVhaRkZFcvHixwF9LUxX8i4YNG1KvXj0++eQT27aQkBC6dOnChAkTDEwmIrcrODj4hn+ZfvjhhwwZMqSIE4lIQTl16hS+vr6sW7eOFi1aFOjvu4483URmZiY7duwgMjIy1/bIyEg2bdpkUCoRKSgLFiwA4Pfff+fEiRMkJibi7OzMjz/+SL9+/QxOJyK3IyUlBcC2IHBB/r6rPN3E6dOnyc7Oply5crm2lytXjqSkJINSiUhBSUpKwtnZmaZNm+Ln58eZM2fIysqiefPmmm1cxIFZrVaGDx9Os2bNCAsLAwr29925MEIXNyaTKdfXVqv1mm0i4nj27t1LjRo1bG+cUVFRlC1b9po/mETEsQwePJg9e/awYcMG27aC/H1XeboJHx8fnJycrjnKlJycrDdXkWJgz5491KpVy/Z1VFQUtWvXNjCRiNyuF154gcWLF7N+/XoqVqxo216Qv+86bXcTLi4u1K9fn1WrVuXavmrVKpo0aWJQKhEpKHv27Mn15qnyJOK4rFYrgwcP5ueff2b16tUEBgbmur8gf99Vnv7F8OHDmTlzJl9++SUxMTEMGzaMxMRE+vfvb3Q0EbkNFouFffv25Xrz/Pvvv6lcubKBqUTkVg0aNIhvvvmGuXPnUrp0aZKSkkhKSuLSpUsF/vuuqQryYPr06UycOJETJ04QFhbG5MmTadGihdGxROQ2xMfHU6NGDQ4fPkylSpUA6NSpExs2bGDRokX6HRdxMDcaizxr1iyaNm1aoL/vKk8iIiIi+aDTdiIiIiL5oPIkIiIikg8qTyIiIiL5oPIkIiIikg8qTyIiIiL5oPIkIiIikg8qTyIiIiL5oPIkIiIikg8qTyIiIiL5oPIkIiIikg8qTyJSIrz00kt06tTJ6BgiUgyoPIlIiRAVFUXdunVvuk/v3r35z3/+UzSBRMRhqTyJSImwe/duwsPDb3i/xWJh6dKldO7cuQhTiYgjUnkSkWLvyJEjnDlzxnbk6fz583Tq1IkmTZpw4sQJADZu3IjZbKZhw4YAvP7669SqVYs777yTcuXKMWDAAK5cuWLUtyAidkTlSUSKvaioKLy8vAgMDGTv3r1ERERQvnx51q5dS/ny5QFYvHgxnTp1wmw2Y7Vayc7O5tNPP2X//v3Mnj2bn376iZkzZxr8nYiIPXA2OoCISGGLioqiTp06fPfddwwaNIh33nmHfv365dpn8eLFvP/++wCYTCbGjx9vu69y5co8+OCDxMbGFmluEbFPOvIkIsVeVFQUe/fuZfDgwSxduvSa4hQTE8PRo0d54IEHADh8+DCDBw8mLCyMu+66Cw8PD3744QcqVqxoRHwRsTMqTyJS7EVFRdG1a1cyMjI4f/78NfcvXryYBx98kDvuuIPTp09z7733cvr0aT744AM2bNjA5s2bcXJy+ter9USkZNBpOxEp1i5cuEBCQgIDBw6kadOm9OjRg02bNhEaGmrbZ9GiRTz33HMALFu2jKysLL777jtMJhMA06ZNIzMzU+VJRACVJxEp5qKionBycqJmzZqEh4ezb98+OnXqxJ9//omPjw/Jycls27aNhQsXAlCmTBlSU1NZvHgxNWvWZMmSJUyYMIEKFSpQtmxZY78ZEbELOm0nIsXa7t27CQ4OxtXVFYB3332XmjVr8sgjj5CZmcmSJUto2LAhvr6+AHTo0IFnn32Wnj170qxZM44dO0b37t111ElEbExWq9VqdAgREaM89NBDNGvWjFdeecXoKCLiIHTkSURKtGbNmtGjRw+jY4iIA9GRJxEREZF80JEnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXxQeRIRERHJB5UnERERkXz4f+Vq/hDWzQe3AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set number of k-points\n",
    "nk = 100\n",
    "ks = np.linspace(0, 2*np.pi, nk, endpoint=False) \n",
    "hamiltonians_0 = pymf.tb_to_khamvector(h_0, nk, ks=ks) \n",
    "\n",
    "# Perform diagonalization\n",
    "vals, vecs = np.linalg.eigh(hamiltonians_0)\n",
    "# Plot data\n",
    "plt.plot(ks, vals, c=\"k\")\n",
    "plt.xticks([0, np.pi, 2 * np.pi], [\"$0$\", \"$\\pi$\", \"$2\\pi$\"])\n",
    "plt.xlim(0, 2 * np.pi)\n",
    "plt.ylabel(\"$E - E_F$\")\n",
    "plt.xlabel(\"$k / a$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ec53b08-053b-4aad-87a6-525dd7f61687",
   "metadata": {},
   "source": [
    "Here, in the workflow to find the ground state, we use a helper function to build the initial guess. because we don't need a dense k-point grid in the self-consistent loop, we compute the spectrum later on a denser k-point grid."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc59e440-1289-4735-9ae8-b04d0d13c94a",
   "metadata": {},
   "source": [
    "Finally, we compute the eigen0alues for a set of Ualues of $U$. For this case, since the interaction is onsite only, the interaction matrix is simply\n",
    "$$\n",
    "H_{int} =\n",
    "\\left(\\begin{array}{cccc}\n",
    "    U & U & 0 & 0\\\\\n",
    "    U & U & 0 & 0\\\\\n",
    "    0 & 0 & U & U\\\\\n",
    "    0 & 0 & U & U\n",
    "\\end{array}\\right)~.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "32b9e7c5-db12-44f9-930c-21e5494404b8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def compute_sol(U, h_0, nk, filling=2):\n",
    "    h_int = {\n",
    "        (0,): U * np.kron(np.eye(2), np.ones((2, 2))),\n",
    "    }\n",
    "    guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
    "    full_model = pymf.Model(h_0, h_int, filling)\n",
    "    mf_sol = pymf.solver(full_model, guess, nk=nk, optimizer_kwargs={\"M\":0})\n",
    "    return pymf.add_tb(h_0, mf_sol)\n",
    "\n",
    "\n",
    "def compute_gap_and_vals(full_sol, nk_dense, fermi_energy=0):\n",
    "    h_kgrid = pymf.tb_to_khamvector(full_sol, nk_dense)\n",
    "    vals = np.linalg.eigvalsh(h_kgrid)\n",
    "\n",
    "    emax = np.max(vals[vals <= fermi_energy])\n",
    "    emin = np.min(vals[vals > fermi_energy])\n",
    "    return np.abs(emin - emax), vals\n",
    "\n",
    "\n",
    "def compute_phase_diagram(\n",
    "    Us,\n",
    "    nk,\n",
    "    nk_dense,\n",
    "):\n",
    "    gaps = []\n",
    "    vals = []\n",
    "    for U in tqdm(Us):\n",
    "        full_sol = compute_sol(U, h_0, nk)\n",
    "        gap, _vals = compute_gap_and_vals(full_sol, nk_dense)\n",
    "        gaps.append(gap)\n",
    "        vals.append(_vals)\n",
    "\n",
    "    return np.asarray(gaps, dtype=float), np.asarray(vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6a8c08a9-7e31-420b-b6b4-709abfb26793",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 75%|███████▌  | 15/20 [00:02<00:00,  5.07it/s]/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.07651e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      "/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.54968e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      "/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.17567e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      "/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.88621e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      "/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.35564e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      "/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.73285e-17): result may not be accurate.\n",
      "  gamma = solve(self.a, df_f)\n",
      " 85%|████████▌ | 17/20 [00:06<00:01,  2.47it/s]\n"
     ]
    },
    {
     "ename": "NoConvergence",
     "evalue": "[ 1.86715628e+00  4.47287535e-02 -1.23356829e-04  1.20430729e-03\n  4.47287535e-02  2.63543104e+00  2.63234161e-03  3.71702031e-03\n -1.23356829e-04  2.63234161e-03  1.86002862e+00  4.51677226e-02\n  1.20430729e-03  3.71702031e-03  4.51677226e-02  2.64245378e+00\n  0.00000000e+00  1.61192759e-01 -1.51104882e-03  1.47113079e-03\n -1.61192759e-01  0.00000000e+00  5.03808388e-03  1.99828384e-03\n  1.51104882e-03 -5.03808388e-03  0.00000000e+00  1.64032886e-01\n -1.47113079e-03 -1.99828384e-03 -1.64032886e-01  0.00000000e+00]",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNoConvergence\u001b[0m                             Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/1169016709.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mUs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0.5\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mendpoint\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk_dense\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m40\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mgap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompute_phase_diagram\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mUs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mUs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk_dense\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnk_dense\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py\u001b[0m in \u001b[0;36mcompute_phase_diagram\u001b[0;34m(Us, nk, nk_dense)\u001b[0m\n\u001b[1;32m     26\u001b[0m     \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     27\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mU\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtqdm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mUs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 28\u001b[0;31m         \u001b[0mfull_sol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompute_sol\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mU\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     29\u001b[0m         \u001b[0mgap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_vals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompute_gap_and_vals\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfull_sol\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk_dense\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     30\u001b[0m         \u001b[0mgaps\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgap\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py\u001b[0m in \u001b[0;36mcompute_sol\u001b[0;34m(U, h_0, nk, filling)\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mguess\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpymf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgenerate_guess\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfrozenset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh_int\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh_0\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalues\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0mfull_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpymf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mh_int\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilling\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m     \u001b[0mmf_sol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpymf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolver\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfull_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mguess\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mpymf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_tb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mh_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmf_sol\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/kwant-scf/pymf/solvers.py\u001b[0m in \u001b[0;36msolver\u001b[0;34m(Model, mf_guess, nk, optimizer, optimizer_kwargs)\u001b[0m\n\u001b[1;32m     60\u001b[0m     \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpartial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcost\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mModel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     61\u001b[0m     result = rparams_to_tb(\n\u001b[0;32m---> 62\u001b[0;31m         \u001b[0moptimizer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmf_params\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0moptimizer_kwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mh_int\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     63\u001b[0m     )\n\u001b[1;32m     64\u001b[0m     \u001b[0mfermi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcalculate_fermi_energy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0madd_tb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mModel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mh_0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mModel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilling\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py\u001b[0m in \u001b[0;36manderson\u001b[0;34m(F, xin, iter, alpha, w0, M, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py\u001b[0m in \u001b[0;36mnonlin_solve\u001b[0;34m(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)\u001b[0m\n\u001b[1;32m    239\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    240\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mraise_exception\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 241\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mNoConvergence\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_array_like\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    242\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    243\u001b[0m             \u001b[0mstatus\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNoConvergence\u001b[0m: [ 1.86715628e+00  4.47287535e-02 -1.23356829e-04  1.20430729e-03\n  4.47287535e-02  2.63543104e+00  2.63234161e-03  3.71702031e-03\n -1.23356829e-04  2.63234161e-03  1.86002862e+00  4.51677226e-02\n  1.20430729e-03  3.71702031e-03  4.51677226e-02  2.64245378e+00\n  0.00000000e+00  1.61192759e-01 -1.51104882e-03  1.47113079e-03\n -1.61192759e-01  0.00000000e+00  5.03808388e-03  1.99828384e-03\n  1.51104882e-03 -5.03808388e-03  0.00000000e+00  1.64032886e-01\n -1.47113079e-03 -1.99828384e-03 -1.64032886e-01  0.00000000e+00]"
     ]
    }
   ],
   "source": [
    "# Interaction strengths\n",
    "Us = np.linspace(0.5, 10, 20, endpoint=True)\n",
    "nk, nk_dense = 40, 100\n",
    "gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e17fc96c-c463-4e1f-8250-c254d761b92a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "\n",
    "ds = xr.Dataset(\n",
    "    data_vars=dict(vals=([\"Us\", \"ks\", \"n\"], vals), gap=([\"Us\"], gap)),\n",
    "    coords=dict(\n",
    "        Us=Us,\n",
    "        ks=np.linspace(0, 2 * np.pi, nk_dense),\n",
    "        n=np.arange(vals.shape[-1])\n",
    "    ),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a87dcc1-208b-4602-abad-a870037ec95f",
   "metadata": {},
   "source": [
    "\n",
    "We observe that as the interaction strength increases, a gap opens due to the antiferromagnetic ordering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "50f02d06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fe0a154a130>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAGdCAYAAABO2DpVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAgGklEQVR4nO3db2xb5d3/8Y9tht0ixyyZUjs3KZgKqbgeg1AylZbCftAqAyyYJqYxMhBs0hSl0FJpKv82Y0YbFbau0iqCUk0ILSr0wYCRB0RDsLXrAKU0FAjZqBgRVOAo7A6yAyxBtc/vQRffNYnb2L2Oj/+8X1Ie+OSk/uKY+K1zfC67LMuyBAAAYIDb6QEAAEDtICwAAIAxhAUAADCGsAAAAMYQFgAAwBjCAgAAGENYAAAAYwgLAABgzBnlvsNsNquPP/5Yfr9fLper3HcPAABKYFmWpqam1NLSIre78HGJsofFxx9/rNbW1nLfLQAAMODo0aM655xzCn6/7GHh9/slHR+soaGh3HcPAABKkE6n1dramnsdL6TsYTF7+qOhoYGwAACgypzqbQy8eRMAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMCYsi+QBQBAPcpkLQ2NTWpialrNfp/aw43yuGvvM7MICwAAbDY4klRiYFTJ1HRuWyjgUzwWUUc05OBk5nEqBAAAGw2OJNXVP5wXFZI0nppWV/+wBkeSDk1mD8ICAACbZLKWEgOjsub53uy2xMCoMtn59qhOhAUAADYZGpucc6TiRJakZGpaQ2OT5RvKZoQFAAA2mZgqHBWl7FcNCAsAAGzS7PcZ3a8aEBYAANikPdyoUMCnQheVunT86pD2cGM5x7IVYQEAgE08bpfisYgkzYmL2dvxWKSm1rMgLAAAsFFHNKTezjYFA/mnO4IBn3o722puHQsWyAIAwGYd0ZDWRYKsvAkAAMzwuF1atazJ6TFsx6kQAABgDGEBAACMISwAAIAxhAUAADCGsAAAAMYQFgAAwBjCAgAAGENYAAAAY1ggCwCA/8pkrbpYHdNOhAUAAJIGR5JKDIwqmZrObQsFfIrHIjX3eR524lQIAKDuDY4k1dU/nBcVkjSemlZX/7AGR5IOTVZ9CAsAQF3LZC0lBkZlzfO92W2JgVFlsvPtga8iLAAAdW1obHLOkYoTWZKSqWkNjU2Wb6gqRlgAAOraxFThqChlv3pHWAAA6lqz32d0v3pHWAAA6lp7uFGhgE+FLip16fjVIe3hxnKOVbUICwBAXfO4XYrHIpI0Jy5mb8djEdazWCDCAgBQ9zqiIfV2tikYyD/dEQz41NvZxjoWRWCBLAAAdDwu1kWCrLx5mggLAAD+y+N2adWyJqfHqGqcCgEAAMYQFgAAwBjCAgAAGENYAAAAYwgLAABgDGEBAACMISwAAIAxhAUAADCGsAAAAMYUFRbHjh3TAw88oHA4rEWLFun888/XQw89pGw2a9d8AACgihS1pPf27dv1+OOP68knn9SKFSv0+uuv6/bbb1cgENDGjRvtmhEAAFSJosLi1Vdf1Q033KDrrrtOknTeeefpqaee0uuvv27LcAAAoLoUdSpkzZo1eumll3TkyBFJ0ptvvqkDBw7o2muvtWU4AABQXYo6YrFlyxalUiktX75cHo9HmUxGW7du1c0331zwZ2ZmZjQzM5O7nU6nS58WAABUtKKOWOzdu1f9/f3as2ePhoeH9eSTT+rXv/61nnzyyYI/09PTo0AgkPtqbW097aEBAEBlclmWZS1059bWVt1zzz3q7u7ObXv44YfV39+vf/7zn/P+zHxHLFpbW5VKpdTQ0HAaowMAgHJJp9MKBAKnfP0u6lTIF198Ibc7/yCHx+M56eWmXq9XXq+3mLsBAABVqqiwiMVi2rp1q5YuXaoVK1bojTfe0I4dO3THHXfYNR8AAKgiRZ0KmZqa0i9+8Qs9++yzmpiYUEtLi26++Wb98pe/1Jlnnrmgf2Ohh1IAAEDlWOjrd1FhYQJhAQBA9Vno6zefFQIAAIwhLAAAgDGEBQAAMIawAAAAxhAWAADAGMICAAAYQ1gAAABjCAsAAGAMYQEAAIwhLAAAgDGEBQAAMIawAAAAxhAWAADAGMICAAAYQ1gAAABjCAsAAGAMYQEAAIwhLAAAgDGEBQAAMIawAAAAxhAWAADAGMICAAAYQ1gAAABjCAsAAGAMYQEAAIwhLAAAgDFnOD0AAKB8MllLQ2OTmpiaVrPfp/Zwozxul9NjoYYQFgBQJwZHkkoMjCqZms5tCwV8isci6oiGHJwMtYRTIQBQBwZHkurqH86LCkkaT02rq39YgyNJhyZDrSEsAKDGZbKWEgOjsub53uy2xMCoMtn59gCKQ1gAQI0bGpucc6TiRJakZGpaQ2OT5RsKNYuwAIAaNzFVOCpK2Q84GcICAGpcs99ndD/gZAgLAKhx7eFGhQI+Fbqo1KXjV4e0hxvLORZqFGEBADXO43YpHotI0py4mL0dj0VYzwJGEBYAUAc6oiH1drYpGMg/3REM+NTb2cY6FjCGBbIAoE50RENaFwmy8iZsRVgAQB3xuF1atazJ6TFQwzgVAgAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAwLZAFAhclkLVbHRNUiLACgggyOJJUYGFUyNZ3bFgr4FI9F+DwPVAVOhQBAhRgcSaqrfzgvKiRpPDWtrv5hDY4kHZoMWDjCAgAqQCZrKTEwKmue781uSwyMKpOdbw+gchAWAFABhsYm5xypOJElKZma1tDYZPmGAkpAWABABZiYKhwVpewHOIWwAIAK0Oz3Gd0PcAphAQAVoD3cqFDAp0IXlbp0/OqQ9nBjOccCikZYAEAF8LhdiscikjQnLmZvx2MR1rNAxSMsAKBCdERD6u1sUzCQf7ojGPCpt7ONdSxQFVggCwAqSEc0pHWRICtvomoRFgBQYTxul1Yta3J6DKAknAoBAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYEzRYfHRRx+ps7NTTU1NWrx4sS6++GIdOnTIjtkAAECVKepy008//VSrV6/Wd77zHb3wwgtqbm7Wv/71L5199tk2jQcAAKpJUWGxfft2tba26oknnshtO++880zPBAAAqlRRp0Kef/55rVy5UjfddJOam5t1ySWXaPfu3XbNBgAVKZO19Oq//ld/OvyRXv3X/yqTtZweCagYRR2xeP/999Xb26vNmzfrvvvu09DQkO666y55vV7deuut8/7MzMyMZmZmcrfT6fTpTQwADhocSSoxMKpkajq3LRTwKR6L8FkegCSXZVkLTu0zzzxTK1eu1CuvvJLbdtddd+ngwYN69dVX5/2ZBx98UIlEYs72VCqlhoaGEkYGAGcMjiTV1T+sr/7RnP0UDz4oDLUsnU4rEAic8vW7qFMhoVBIkUgkb9uFF16oDz/8sODP3HvvvUqlUrmvo0ePFnOXAFARMllLiYHROVEhKbctMTDKaRHUvaJOhaxevVrvvvtu3rYjR47o3HPPLfgzXq9XXq+3tOkAoEIMjU3mnf74KktSMjWtobFJPkAMda2oIxZ33323XnvtNW3btk3vvfee9uzZo76+PnV3d9s1HwBUhImpwlFRyn5ArSoqLC677DI9++yzeuqppxSNRvWrX/1KO3fu1C233GLXfABQEZr9PqP7AbWqqFMhknT99dfr+uuvt2MWAKhY7eFGhQI+jaem532fhUtSMOBTe7ix3KMBFYXPCgGABfC4XYrHjr953fWV783ejsci8ri/+l2gvhAWALBAHdGQejvbFAzkn+4IBnxcagr8V9GnQgCgnnVEQ1oXCWpobFITU9Nq9h8//cGRCuA4wgIAiuRxu7ikFCiAUyEAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhzhtMDAKhMmaylobFJTUxNq9nvU3u4UR63y+mxAFQ4wgLAHIMjSSUGRpVMTee2hQI+xWMRdURDDk4GoNJxKgRAnsGRpLr6h/OiQpLGU9Pq6h/W4EjSockAVAPCAkBOJmspMTAqa57vzW5LDIwqk51vDwAgLACcYGhscs6RihNZkpKpaQ2NTZZvKABVhbAAkDMxVTgqStkPQP0hLADkNPt9RvcDUH8ICwA57eFGhQI+Fbqo1KXjV4e0hxvLORaAKkJYAMjxuF2KxyKSNCcuZm/HYxHWswBQ0GmFRU9Pj1wulzZt2mRoHABO64iG1NvZpmAg/3RHMOBTb2cb61gAOKmSF8g6ePCg+vr6dNFFF5mcB0AF6IiGtC4SZOVNAEUr6YjFZ599pltuuUW7d+/W17/+ddMzAagAHrdLq5Y16YaL/0erljURFQAWpKSw6O7u1nXXXadrrrnG9DwAAKCKFX0q5Omnn9bw8LAOHjy4oP1nZmY0MzOTu51Op4u9SwAAUCWKOmJx9OhRbdy4Uf39/fL5FnYde09PjwKBQO6rtbW1pEEBAEDlc1mWteBF/5977jl973vfk8fjyW3LZDJyuVxyu92amZnJ+540/xGL1tZWpVIpNTQ0GPhPAAAAdkun0woEAqd8/S7qVMjVV1+tt99+O2/b7bffruXLl2vLli1zokKSvF6vvF5vMXcDAACqVFFh4ff7FY1G87adddZZampqmrMdAADUH1beBAAAxpS8QNasv/71rwbGAAAAtYAjFgAAwJjTPmIBwBmZrMWS2wAqDmEBVKHBkaQSA6NKpqZz20IBn+KxCB8SBsBRnAoBqszgSFJd/cN5USFJ46lpdfUPa3Ak6dBkAEBYAFUlk7WUGBjVfKvazW5LDIwqk13wuncAYBRhAVSRobHJOUcqTmRJSqamNTQ2Wb6hAOAEhAVQRSamCkdFKfsBgGmEBVBFmv0L+/C/he4HAKYRFkAVaQ83KhTwqdBFpS4dvzqkPdxYzrEAIIewAKqIx+1SPBaRpDlxMXs7HouwngUAxxAWQJXpiIbU29mmYCD/dEcw4FNvZxvrWABwFAtkAVWoIxrSukiQlTcBVBzCAqhSHrdLq5Y1OT0GAOThVAgAADCGsAAAAMYQFgAAwBjCAgAAGENYAAAAYwgLAABgDGEBAACMISwAAIAxhAUAADCGsAAAAMawpDdgk0zW4rM8ANQdwgKwweBIUomBUSVT07ltoYBP8ViETx8FUNM4FQIYNjiSVFf/cF5USNJ4alpd/cMaHEk6NBkA2I+wAAzKZC0lBkZlzfO92W2JgVFlsvPtAQDVj7AADBoam5xzpOJElqRkalpDY5PlGwoAyoiwAAyamCocFaXsBwDVhrAADGr2+4zuBwDVhrAADGoPNyoU8KnQRaUuHb86pD3cWM6xAKBsCAvAII/bpXgsIklz4mL2djwWYT0LADWLsAAM64iG1NvZpmAg/3RHMOBTb2cb61gAqGkskAXYoCMa0rpIkJU3AdQdwgKwicft0qplTU6PAQBlxakQAABgDGEBAACMISwAAIAxhAUAADCGsAAAAMYQFgAAwBjCAgAAGENYAAAAYwgLAABgDGEBAACMYUlv1LVM1uLzPADAIMICdWtwJKnEwKiSqenctlDAp3gswieQAkCJOBWCujQ4klRX/3BeVEjSeGpaXf3DGhxJOjQZAFQ3wgJ1J5O1lBgYlTXP92a3JQZGlcnOtwcA4GQIC9SdobHJOUcqTmRJSqamNTQ2Wb6hAKBGEBaoOxNThaOilP0AAP+HsEDdafb7jO4HAPg/hAXqTnu4UaGAT4UuKnXp+NUh7eHGco4FADWBsEDd8bhdiscikjQnLmZvx2MR1rMAgBIQFqhLHdGQejvbFAzkn+4IBnzq7WxjHQsAKBELZKFudURDWhcJsvImABhEWKCuedwurVrW5PQYAFAzOBUCAACMISwAAIAxhAUAADCGsAAAAMYQFgAAwBjCAgAAGFNUWPT09Oiyyy6T3+9Xc3OzbrzxRr377rt2zQYAAKpMUWGxb98+dXd367XXXtOLL76oY8eOaf369fr888/tmg8AAFQRl2VZVqk//Mknn6i5uVn79u3T2rVrF/Qz6XRagUBAqVRKDQ0Npd41AAAoo4W+fp/WeyxSqZQkqbGRT4EEAACnsaS3ZVnavHmz1qxZo2g0WnC/mZkZzczM5G6n0+lS7xIAAFS4ksNiw4YNeuutt3TgwIGT7tfT06NEIlHq3aDOZbIWHxIGAFWkpPdY3HnnnXruuee0f/9+hcPhk+473xGL1tZW3mOBUxocSSoxMKpkajq3LRTwKR6L8LHmAFBmtrzHwrIsbdiwQc8884xefvnlU0aFJHm9XjU0NOR9AacyOJJUV/9wXlRI0nhqWl39wxocSTo0GQDgZIoKi+7ubvX392vPnj3y+/0aHx/X+Pi4/vOf/9g1H+pQJmspMTCq+Q6lzW5LDIwqky35giYAgE2KCove3l6lUildddVVCoVCua+9e/faNR/q0NDY5JwjFSeyJCVT0xoamyzfUACABSnqzZunseQFsGATU4WjopT9AADlw2eFoOI0+31G9wMAlA9hgYrTHm5UKOBToYtKXTp+dUh7mIXZAKDSEBaoOB63S/FYRJLmxMXs7XgswnoWAFCBCAtUpI5oSL2dbQoG8k93BAM+9Xa2sY4FAFSoklfeBOzWEQ1pXSTIypsAUEUIC1Q0j9ulVcuanB4DALBAnAoBAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMYQFAAAwhrAAAADGEBYAAMAYwgIAABhDWAAAAGMICwAAYAxhAQAAjCEsAACAMWc4PQDslclaGhqb1MTUtJr9PrWHG+Vxu5weCwBQowiLGjY4klRiYFTJ1HRuWyjgUzwWUUc05OBkAIBaxamQGjU4klRX/3BeVEjSeGpaXf3DGhxJOjQZAKCWERY1KJO1lBgYlTXP92a3JQZGlcnOtwcAAKUjLGrQ0NjknCMVJ7IkJVPTGhqbLN9QAIC6QFjUoImpwlFRyn4AACwUYVGDmv0+o/sBALBQhEUNag83KhTwqdBFpS4dvzqkPdxYzrEAAHWAsKhBHrdL8VhEkubExezteCzCehYAAOMIixrVEQ2pt7NNwUD+6Y5gwKfezjbWsQAA2IIFsmpYRzSkdZEgK28CAMqGsKhxHrdLq5Y1OT0GAKBOcCoEAAAYQ1gAAABjCAsAAGAMYQEAAIzhzZsVIJO1uHIDAFATCAuHDY4klRgYzfvQsFDAp3gswloTAICqw6kQBw2OJNXVPzznk0jHU9Pq6h/W4EjSockAACgNYeGQTNZSYmBU1jzfm92WGBhVJjvfHgAAVCbCwiFDY5NzjlScyJKUTE1raGyyfEMBAHCaCAuHTEwVjopS9gMAoBIQFg5p9vtOvVMR+wEAUAm4KmQB7LgctD3cqFDAp/HU9Lzvs3Dp+CeRtocbT+t+7MalsigFzxugdhEWp2DX5aAet0vxWERd/cNySXlxMfvnNR6LVPQfWy6VRSnK8bwhXGoXv9vK57Isq6yXHaTTaQUCAaVSKTU0NBj5N+16os1eDvrVB2j2X+7tbDvtP4TV+uJcjscGtYf/p2qfnS/8/G5Pzu7oWujrd0lh8dhjj+nRRx9VMpnUihUrtHPnTl1xxRVGB1sou55omaylNdtfLnjlxuypigNb/t9p/+KqrcDL+digdpTjeUPwOsvOF35+tydXjuha6Ot30W/e3Lt3rzZt2qT7779fb7zxhq644gp997vf1YcffnhaA5fCzgWmynk5qMft0qplTbrh4v/RqmVNFf9izKWyKIXdzxvWhnGWnX+P+d2eXKUttlh0WOzYsUM/+clP9NOf/lQXXnihdu7cqdbWVvX29toxX0F2P9G4HLQwHhuUwu7nDcHrHLv/HvO7LawSo6uosPjyyy916NAhrV+/Pm/7+vXr9corrxgd7FTsfqJxOWhhPDYohd3PG4LXOXb/PeZ3W1glRldRYfHvf/9bmUxGS5Ysydu+ZMkSjY+Pz/szMzMzSqfTeV8m2P1Em70ctNBJCZeOn7+q9MtB7cBjg1LY/bwheJ1j999jfreFVWJ0lbRAlsuV/6fBsqw522b19PQoEAjkvlpbW0u5yznsfqLNXg4qac4fwmq5HNQuPDYohd3PG4LXOXb/PeZ3W1glRldRYfGNb3xDHo9nztGJiYmJOUcxZt17771KpVK5r6NHj5Y+7QnK8UTriIbU29mmYCD/FxIM+Or+Hcg8NiiFnc8bgtc5dv895ndbWCVGV9GXm37729/WpZdeqsceeyy3LRKJ6IYbblBPT88pf97k5aaz74SV5l9gytQLXLVdDlpOPDYoBWsd1J5y/D3mdzu/cr0W2raOxd69e/XjH/9Yjz/+uFatWqW+vj7t3r1b77zzjs4991xjgy0UTzQAX0XwOoNVVZ1TSetYlLxA1iOPPKJkMqloNKrf/va3Wrt2rdHBisETDQAqA3+PnVPVK2+eDjvCAgAA2Mu2lTcBAAAKISwAAIAxhAUAADCGsAAAAMYQFgAAwBjCAgAAGENYAAAAYwgLAABgDGEBAACMOaPcdzi70Gc6nS73XQMAgBLNvm6fasHusofF1NSUJKm1tbXcdw0AAE7T1NSUAoFAwe+X/bNCstmsPv74Y/n9frlcfDBNOaTTabW2turo0aN8PkuZ8dg7h8feOTz2zrHzsbcsS1NTU2ppaZHbXfidFGU/YuF2u3XOOeeU+24hqaGhgf/JHcJj7xwee+fw2DvHrsf+ZEcqZvHmTQAAYAxhAQAAjCEs6oDX61U8HpfX63V6lLrDY+8cHnvn8Ng7pxIe+7K/eRMAANQujlgAAABjCAsAAGAMYQEAAIwhLAAAgDGERQ3r6enRZZddJr/fr+bmZt1444169913nR6rLvX09MjlcmnTpk1Oj1IXPvroI3V2dqqpqUmLFy/WxRdfrEOHDjk9Vs07duyYHnjgAYXDYS1atEjnn3++HnroIWWzWadHqzn79+9XLBZTS0uLXC6XnnvuubzvW5alBx98UC0tLVq0aJGuuuoqvfPOO2WZjbCoYfv27VN3d7dee+01vfjiizp27JjWr1+vzz//3OnR6srBgwfV19eniy66yOlR6sKnn36q1atX62tf+5peeOEFjY6O6je/+Y3OPvtsp0eredu3b9fjjz+uXbt26R//+IceeeQRPfroo/rd737n9Gg15/PPP9e3vvUt7dq1a97vP/LII9qxY4d27dqlgwcPKhgMat26dbnP67ITl5vWkU8++UTNzc3at2+f1q5d6/Q4deGzzz5TW1ubHnvsMT388MO6+OKLtXPnTqfHqmn33HOP/v73v+tvf/ub06PUneuvv15LlizR73//+9y273//+1q8eLH+8Ic/ODhZbXO5XHr22Wd14403Sjp+tKKlpUWbNm3Sli1bJEkzMzNasmSJtm/frp/97Ge2zsMRizqSSqUkSY2NjQ5PUj+6u7t13XXX6ZprrnF6lLrx/PPPa+XKlbrpppvU3NysSy65RLt373Z6rLqwZs0avfTSSzpy5Igk6c0339SBAwd07bXXOjxZfRkbG9P4+LjWr1+f2+b1enXllVfqlVdesf3+y/4hZHCGZVnavHmz1qxZo2g06vQ4deHpp5/W8PCwDh486PQodeX9999Xb2+vNm/erPvuu09DQ0O666675PV6deuttzo9Xk3bsmWLUqmUli9fLo/Ho0wmo61bt+rmm292erS6Mj4+LklasmRJ3vYlS5bogw8+sP3+CYs6sWHDBr311ls6cOCA06PUhaNHj2rjxo3685//LJ/P5/Q4dSWbzWrlypXatm2bJOmSSy7RO++8o97eXsLCZnv37lV/f7/27NmjFStW6PDhw9q0aZNaWlp02223OT1e3XG5XHm3Lcuas80OhEUduPPOO/X8889r//79fGR9mRw6dEgTExO69NJLc9symYz279+vXbt2aWZmRh6Px8EJa1coFFIkEsnbduGFF+qPf/yjQxPVj5///Oe655579MMf/lCS9M1vflMffPCBenp6CIsyCgaDko4fuQiFQrntExMTc45i2IH3WNQwy7K0YcMGPfPMM3r55ZcVDoedHqluXH311Xr77bd1+PDh3NfKlSt1yy236PDhw0SFjVavXj3nsuojR47o3HPPdWii+vHFF1/I7c5/WfF4PFxuWmbhcFjBYFAvvvhibtuXX36pffv26fLLL7f9/jliUcO6u7u1Z88e/elPf5Lf78+ddwsEAlq0aJHD09U2v98/570sZ511lpqamniPi83uvvtuXX755dq2bZt+8IMfaGhoSH19ferr63N6tJoXi8W0detWLV26VCtWrNAbb7yhHTt26I477nB6tJrz2Wef6b333svdHhsb0+HDh9XY2KilS5dq06ZN2rZtmy644AJdcMEF2rZtmxYvXqwf/ehH9g9noWZJmvfriSeecHq0unTllVdaGzdudHqMujAwMGBFo1HL6/Vay5cvt/r6+pweqS6k02lr48aN1tKlSy2fz2edf/751v3332/NzMw4PVrN+ctf/jLv3/fbbrvNsizLymazVjwet4LBoOX1eq21a9dab7/9dllmYx0LAABgDO+xAAAAxhAWAADAGMICAAAYQ1gAAABjCAsAAGAMYQEAAIwhLAAAgDGEBQAAMIawAAAAxhAWAADAGMICAAAYQ1gAAABj/j+4y9R0fBf9vgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(Us, gap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "868cf368-45a0-465e-b042-6182ff8b6998",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'_PlotMethods' object has no attribute 'scatter'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3492086229.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvals\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mscatter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"ks\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"Us\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mec\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxhline\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mls\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"--\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"k\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxticks\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m\"$0$\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"$\\pi$\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"$2\\pi$\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"$E - E_F$\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: '_PlotMethods' object has no attribute 'scatter'"
     ]
    }
   ],
   "source": [
    "ds.vals.plot.scatter(x=\"ks\", hue=\"Us\", ec=None, s=5)\n",
    "plt.axhline(0, ls=\"--\", c=\"k\")\n",
    "plt.xticks([0, np.pi, 2 * np.pi], [\"$0$\", \"$\\pi$\", \"$2\\pi$\"])\n",
    "plt.xlim(0, 2 * np.pi)\n",
    "plt.ylabel(\"$E - E_F$\")\n",
    "plt.xlabel(\"$k / a$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0761ed33-c1bb-4f12-be65-cb68629f58b9",
   "metadata": {},
   "source": [
    "The Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))\n",
    "$$\n",
    "\\epsilon_{HF}^{\\sigma}(\\mathbf{k}) = \\epsilon(\\mathbf{k}) + U \\left(\\frac{n}{2} + \\sigma m\\right)\n",
    "$$\n",
    "where $m=(\\langle n_{i\\uparrow} \\rangle - \\langle n_{i\\downarrow} \\rangle) / 2$ is the magnetization per atom and $n = \\sum_i \\langle n_i \\rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\\Delta=U$. And we can confirm it indeed follows the expected trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "ac2eb725-f3bd-4d5b-a509-85d0d0071958",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9+klEQVR4nO3dd3xV9eH/8dfJTiAJewTCFIgBAiGTgIgVsZVSKO4BYalVRAG31lkFcbfiAiIEkGqrRanWukGRkUXYEPY0hJlLBhn3nt8ftvzKV0UIST53vJ+PR/7ISUJfvWru+3HO4V7Ltm0bEREREQ/lZzpARERE5HxozIiIiIhH05gRERERj6YxIyIiIh5NY0ZEREQ8msaMiIiIeDSNGREREfFoAaYD6prL5eLAgQOEh4djWZbpHBERETkLtm1z4sQJoqKi8PM787kXrx8zBw4cIDo62nSGiIiI1MDevXtp27btGb/H68dMeHg48MODERERYbhGREREzobD4SA6OvrU8/iZeP2Y+e+lpYiICI0ZERERD3M2t4joBmARERHxaBozIiIi4tE0ZkRERMSjacyIiIiIR9OYEREREY+mMSMiIiIeTWNGREREPJrGjIiIiHg0jRkRERHxaBozIiIi4tGMjplvvvmGoUOHEhUVhWVZfPDBB6d93bZtHn/8caKioggNDWXgwIFs2LDBTKyIiIi4JaNjprS0lF69ejFjxoyf/Pqzzz7Liy++yIwZM8jOzqZVq1ZcdtllnDhxop5LRURExF0ZHTO/+c1veOqppxgxYsSPvmbbNi+//DIPP/wwI0aMoEePHmRmZlJWVsbChQsN1IqIiMj/sl0u8r98F5fTabTDbe+Z2blzJ4WFhQwePPjUseDgYC6++GKWL1/+sz9XUVGBw+E47UNERERqV4njGDkvXU3vb29h1duPG21x2zFTWFgIQMuWLU873rJly1Nf+ynTpk0jMjLy1Ed0dHSddoqIiPiabWtXcPSlNJJOfEG17YdlWUZ73HbM/Nf/fYBs2z7jg/bggw9SXFx86mPv3r11nSgiIuITbJeLlX97luj3h9LOPsBBmrJtyLukjnzSaFeA0f/1M2jVqhXwwxma1q1bnzpeVFT0o7M1/ys4OJjg4OA67xMREfEljuKjbJ01htSSJWDBmtBU2o/LJKZZK9Np7ntmpmPHjrRq1YrPP//81LHKykqWLl1KWlqawTIRERHfsjX/Wxwv9yWhZAlVtj+rLphM3L2f0MgNhgwYPjNTUlLCtm3bTn2+c+dO8vPzadKkCe3atWPSpElMnTqVLl260KVLF6ZOnUpYWBg33HCDwWoRERHfYLtcrHznGRK2vECQVc33VnNODJ1JSsKvTKedxuiYycnJ4ZJLLjn1+ZQpUwBIT09n7ty53HfffZSXl3P77bdz7NgxUlJS+OyzzwgPDzeVLCIi4hOKjx5i++zR9C1bBhasbtCPTuMzad24uem0H7Fs27ZNR9Qlh8NBZGQkxcXFREREmM4RERFxe5tzvybio1uIsouotP1ZHXM3ydc+iOVXf3ennMvzt9veACwiIiL1y+V0seqvfyJx658JtJzst1pSPmw2KfEDTKedkcaMiIiIcOzwQXZlpNO3fAVYkNfwYrqMn0ObRk1Np/0ijRkREREftynrcxr/6w/Ec5gKO5A13e8j6ap76vWy0vnQmBEREfFRLqeTlQseI3nHqwRYLvZaUVSOeIvkuL6m086JxoyIiIgPOlK0n31vjSLtZA5YkBMxiJjxs2kY0dh02jnTmBEREfEx65d/QovPbqcXRzlpB7Iu7mESf3+Xx1xW+r80ZkRERHyEs7qaVfMeJmX3m/hbNrv92mJfNYek2GTTaedFY0ZERMQHHCrcQ+GcUaRVrP7hslKjXxM7fiZhDSNNp503jRkREREvt+7bxbT+ciI9OU6ZHczG+EdJHH6H6axaozEjIiLipaqrqsjOfICUvRn4WTa7/NphXZNJYkwf02m1SmNGRETECxXt30VR5kj6Vq4FC7KaDCVu/OuEhHnf+xtqzIiIiHiZNUveJ3rJJHrgoNQOYVPikyQPvdV0Vp3RmBEREfESVVWVZM+5h7QDmQBs9+9I0HXzSOwSZ7isbmnMiIiIeIHCvds5Om8kaVUbAMhqNpy4ca8REtrAcFnd05gRERHxcKu/eIcOy+4hlhOU2KEUpE4l+TdjTWfVG40ZERERD1VZUUHuW5Poe3AhANv8OxN24wL6dIo1XFa/NGZEREQ80IFdW3AsGEXf6s0ArGp+NfHjXiEoJNRwWf3TmBEREfEweZ8u4IIV9xFFKQ7C2NF3OimXjzKdZYzGjIiIiIeoqChn9ew7ST30NwAKAroSftN8eneIMVxmlsaMiIiIB9i3fSNlC0eR6twKwKpWN9Bn7EsEBoUYLjNPY0ZERMTN5X4yh64rH6StVc5xGrL7oudJufR601luQ2NGRETETZ0sLyV/9gRSjywCCzYHxtJ41Dx6RXcxneZWNGZERETc0J6ta6l8J51U5w4AVkaNInH08wQEBRsucz8aMyIiIm4m+58zic15hAbWSY4Rwd6BL5E68CrTWW5LY0ZERMRNlJeWsHb2H0g59k+wYGNQT5qnzyeuTUfTaW5NY0ZERMQN7Nq8Gvtv6aS4duOyLbKix5A0+ln8AwJNp7k9jRkRERHDsj6YQY/VTxJmVXCESAoHvULqRcNMZ3kMjRkRERFDSk8Us2H2LSQX/xss2BDcm1aj59O9dTvTaR5FY0ZERMSAHRuz8XtvDMmuvThti5wOt5A4cir+AXpqPld6xEREROqR7XKRtegvxK19mlCrkkM05tDlr5KSNsR0msfSmBEREaknJY5jbJo1npQTX4AF60ISaDNmHrEt25pO82gaMyIiIvVg29oVBC0aS5J9gGrbj5xOt5N805P4+fubTvN4GjMiIiJ1yHa5WPXeC8RvmE6wVUURTTg65A1Sky83neY1NGZERETqiKP4KFtnjSG1ZAlYsCY0hfZjM4lp3tp0mlfRmBEREakDW/O/JfTD8STYhVTZ/uR1uZPkGx7B8tNlpdqmMSMiIlKLbJeLle88Q8KWFwiyqimkOcVD3yQl8VLTaV5LY0ZERKSWFB89zPbZo+lb9i1YkB+WRqdxmbRq2sJ0mlfTmBEREakFm3O/JuKjW+hjF1Fp+7M6ZgrJ1z6E5ednOs3racyIiIicB5fTxaq//onErX8m0HJywGpJ2bDZpMQPMJ3mMzRmREREauj44YPszEinb/kKsCCv4cVcMH4OUY2amk7zKRozIiIiNbAp63Ma/+sPxHOYSjuA/O73kXTVvbqsZIDGjIiIyDlwOZ2sXPAYyTteJcBysc9qTeWIt0iOSzOd5rM0ZkRERM7SkaL97HsrnbST2WBBbsSldBufQcOIxqbTfJrGjIiIyFnYsOITmn96O704ykk7kHVxD5H4+0m6rOQGNGZERETOwFldzap5D5Oy+038LZs9fm1wXTmHpO4pptPkPzRmREREfsahwj0UzhlFWsVqsCAn8nJib55JWMNGptPkf2jMiIiI/IR13y6m9ZcT6clxyu0gNsQ/SuLwiaaz5CdozIiIiPyP6qoqsjMfIGVvBn6WzS6/dljXzCUxJsF0mvwMjRkREZH/KNq/i6LMkfStXAsWZDceQs+b3yQkLNx0mpyBxoyIiAiwZsn7RC+ZRA8clNnBbEp8kqShfzCdJWdBY0ZERHxaVVUlOXPuoe+BTAB2+Hcg8NpMErr2NhsmZ01jRkREfFbh3u0cnTeSvlUbAMhqOpy4ca8SEtbQcJmcC40ZERHxSflfvkv7b+8mlhOU2KEUpEwl+YqxprOkBjRmRETEp1RWVJDz1mTSDr4NwDb/zoTeMI8+nXsYLpOa0pgRERGfcWDXFhwLRpFWvRmAVc2vove4VwgOCTNcJudDY0ZERHxC3mcL6Lz8PqIoxUEY2/tOJ+XyUaazpBZozIiIiFerqChn9ew7ST30NwAKAroSftN84jvEGC6T2qIxIyIiXmv/jo2UvT2KVOdWAFa1up4+Y18mMCjEcJnUJo0ZERHxSrmfzKXrygdoY5VTTAN2XfQCKZdebzpL6oDGjIiIeJWT5aXkZ9xB6uF/gAWbA2NpPHIevdp1MZ0mdcTPdMCZVFdX88c//pGOHTsSGhpKp06dePLJJ3G5XKbTRETEDe3ZupZ9z1/0w5ABVkaN4oJ7l9BSQ8arufWZmenTp/PGG2+QmZlJ9+7dycnJYcyYMURGRnLXXXeZzhMRETeS89EsYrIfoaFVzjEi2DvwJVIHXmU6S+qBW4+ZFStWMGzYMIYMGQJAhw4d+Otf/0pOTo7hMhERcRflpSWszbiNlKOLwYKNQT1pnj6fuDYdTadJPXHry0z9+/fnyy+/pKCgAIA1a9awbNkyrrjiip/9mYqKChwOx2kfIiLinXZvWc33L/Qj5ehiXLbFyrZj6Xbf1zTXkPEpbn1m5v7776e4uJiYmBj8/f1xOp08/fTTXH/9z9+NPm3aNJ544ol6rBQREROyPniN7qsfp4FVwREi+X7QX0i9aLjpLDHArc/MvPvuuyxYsICFCxeSl5dHZmYmzz//PJmZmT/7Mw8++CDFxcWnPvbu3VuPxSIiUtfKSorJevl6kvMfpIFVwfrg3ti3fksPDRmfZdm2bZuO+DnR0dE88MADTJgw4dSxp556igULFrB58+az+jMcDgeRkZEUFxcTERFRV6kiIlIPdmzMxu+9MXRw7cVpW2S3v4WkUVPxD3DrCw1SA+fy/O3W//TLysrw8zv95JG/v7/+araIiI+xXS6yFr1C3NqnCLUqOURjDl3+KqlpQ0yniRtw6zEzdOhQnn76adq1a0f37t1ZvXo1L774ImPHjjWdJiIi9aTkxHE2zxpHiuMLsGBdSAJtxswjtmVb02niJtz6MtOJEyd45JFHWLRoEUVFRURFRXH99dfz6KOPEhQUdFZ/hi4ziYh4rm3rVhD0j7G0sw9QbfuR0+l2km96Ej9/f9NpUsfO5fnbrcdMbdCYERHxPLbLxar3XiR+wzMEW1UU0YSjv3mDmJTLTadJPfGae2ZERMT3OIqPUjBrLKklX4MFa0JTaD82k5jmrU2niZvSmBEREbexNf9bQj8cT6JdSJXtT16XO0m+4REsP11Wkp+nMSMiIsbZLher3p1On83PE2RVU0hzioe+SUripabTxANozIiIiFHFRw+zffZoUsu+BQvyw9LoNC6TVk1bmE4TD6ExIyIixmzJ/Zrwj26hj11Epe3P6pgpJF/7EJafW79AvbgZjRkREal3LqeLVX99ioStLxNkOTlgtaRs2GxS4geYThMPpDEjIiL16vjhg+zMSKdv+QqwYHXDAXQeP4eoRs1Mp4mH0pgREZF6synrcxr/6w/Ec5hKO4D87veRdNW9uqwk50VjRkRE6pzL6WTl24+TvH0GAZaLfVZrKkZkkBzXz3SaeAGNGRERqVNHivaz76100k5mgwW54b+i280ZNIxoYjpNvITGjIiI1JkNKz6h+ae304ujnLQDWRf3EIm/n6TLSlKrNGZERKTWOZ1OVs17mJRdb+Bv2ezxa4PzyjkkdU8xnSZeSGNGRERq1aHCPRTOSSetIg8syIm8nNibZxLWsJHpNPFSGjMiIlJr1i9bTKsvJtKT45TbQayPf5Sk4RNNZ4mX05gREZHzVl1VRXbmA6TszcDPstnl1w7rmrkkxSSYThMfoDEjIiLn5dCB3RycexN9K9eCBdmNh9Bj/JuENgg3nSY+QmNGRERqbM2S94leMokeOCizg9mU+ARJQ28znSU+RmNGRETOWXVVJVlz7iV1fyZ+ls0O/w4EXptJQtfeptPEB2nMiIjIOSncu51j80aSVrUBLMhqOpy4ca8SEtbQdJr4KI0ZERE5a/lfvkv7b+/mQk5QYodSkPI0yVeMM50lPk5jRkREflFlRQW5cybTt/BtALb5dyb0hnn06dzDcJmIxoyIiPyCA7u24Fgwir7VmwFY1fwqeo97heCQMMNlIj/QmBERkZ+V99kCOi+/jyhKcRDG9r7TSbl8lOkskdNozIiIyI9UVJSTl3EXfYveBaAgoCsNb5xPfMcYw2UiP6YxIyIip9m/YxOlb4+ir7MAgJUtr6fP2JcJCg4xXCby0zRmRETklNxP5tJl5YO0scoopgG7+r9A6qDrTWeJnJHGjIiIcLK8lPyMO0g9/A+wYEvghTQaOZ9e7bqYThP5RRozIiI+bs/WdVS+k06qczsAK6NGkTj6eQKCgg2XiZwdjRkRER+W89EsYrIfoaFVzjHC2TfwZVIHXmU6S+ScaMyIiPig8tIS1mbcRsrRxWDBpqAeNEtfQM82HU2niZwzjRkRER+ze8tqnO+OJsW1C5dtkRU9hsT06QQEBplOE6kRjRkRER+S9cFrdF/9OA2sCo4QyfeD/kLqRcNNZ4mcF40ZEREfUFZSzPrZfyD5+L/Agg3BvWg5ej49Wrc3nSZy3jRmRES83M6NOVjvjSbZtRenbZHd/haSRk3FP0BPAeId9G+yiIiXsl0usha9Qtzapwi1KjlMI4ouf43UtCGm00RqlcaMiIgXKjlxnE2zxpPi+BwsWBeSQJsx84ht2dZ0mkit05gREfEy29atJOgfY0my9+O0LXI63U7STX/Cz9/fdJpIndCYERHxErbLxar3XiR+wzMEW1UU0YSjV7xOSsqvTaeJ1CmNGRERL+AoPkrBrLGklnwNFqwNTabd2HnENG9tOk2kzmnMiIh4uK353xL64c0k2t9TZfuT12UiyTc8iuWny0riGzRmREQ8lO1yserd6fTZ/DxBVjWFNKd46BukJA4ynSZSrzRmREQ8UPHRw2zLGE1q6bdgQX5YGp3GZdKqaQvTaSL1TmNGRMTDbMldQsOPbiHBPkil7c/qmCkkX/sQlp+f6TQRIzRmREQ8hO1ysfKvT5NQ8BJBlpMDVkvKhs0mJX6A6TQRozRmREQ8wPHDB9mZkU7f8hVgweqGA+g8fg5RjZqZThMxTmNGRMTNbcr6gkb/upV4DlNpB5Afex9JV9+ry0oi/6ExIyLiplxOJyvffpzk7TMIsFzss1pTMSKD5Lh+ptNE3IrGjIiIGzpadIC9b40i7WQ2WJAX/iu63pxBw4gmptNE3I7GjIiIm9mw4hOaf3o7vTjKSTuQdXEPkvj7ybqsJPIzNGZERNyE0+lk1byHSdn1Bv6WzR6/NriunENS9xTTaSJuTWNGRMQNHCrcQ+GcdNIq8sCCnMjLib15JmENG5lOE3F7GjMiIoatX/ZPWn1xBz05TrkdxIbej5I4/A6wLNNpIh5BY0ZExBBndTWrMh8gZc9s/C2b3X7RWNdkkhiTYDpNxKNozIiIGHDowG4Ozh1JWuUasCC78RB63vwmIWHhptNEPI7GjIhIPVu7dBFtvr6LHhRTZgezKfEJkobeZjpLxGNpzIiI1JPqqkqy5t5H6r65+Fk2O/w7EHhtJglde5tOE/FoGjMiIvWgcO92js4bRVrVerAgq+kw4sa9RkhYQ9NpIh5PY0ZEpI6t+epd2n1zN7GcoMQOpSD5KZKHjDedJeI1NGZEROpIVWUFOW9Npm/h2wBs9+9MyA3z6NO5h+EyEe+iMSMiUgcO7C7AMX8Ufas3AZDV/Cp6jfsLwSENDJeJeB+3f6OP/fv3c9NNN9G0aVPCwsLo3bs3ubm5prNERH5W3udv02DOQGKqN+EgjPy+fyF5QoaGjEgdceszM8eOHaNfv35ccsklfPLJJ7Ro0YLt27fTqFEj02kiIj9SUVFOXsZd9C16F4CCgK40vHEevTteaLhMxLu59ZiZPn060dHRzJkz59SxDh06nPFnKioqqKioOPW5w+GoqzwRkVP279hE6cJR9K0uAGBVy2uJH/sXgoJDDJeJeD+3vsy0ePFiEhMTufrqq2nRogXx8fHMmjXrjD8zbdo0IiMjT31ER0fXU62I+KrcTzIJz/wVXasLKKYBa/q/TsptMzVkROqJZdu2bTri54SE/PCLYMqUKVx99dVkZWUxadIk3nzzTUaNGvWTP/NTZ2aio6MpLi4mIiKiXrpFxDecLC9ldcZE+h5+H4AtgRfSaOQ8WrbrarhMxPM5HA4iIyPP6vnbrcdMUFAQiYmJLF++/NSxO++8k+zsbFasWHFWf8a5PBgiImdrz9Z1VLyTThfndgBWtR5JwpgXCAgKNlwm4h3O5fnbrS8ztW7dmtjY2NOOXXjhhezZs8dQkYgIZH80myYLLqOLczvHCWfdxbNJuXWGhoyIIW59A3C/fv3YsmXLaccKCgpo3769oSIR8WXlpSWsybid1KMfggWbArvTbPQCerbpZDpNxKe59ZiZPHkyaWlpTJ06lWuuuYasrCxmzpzJzJkzTaeJiI/ZvSUf57vppLp24bItsqLHkJg+nYDAINNpIj7Pre+ZAfjoo4948MEH2bp1Kx07dmTKlCncfPPNZ/3zumdGRM5X1oev0z3vMRpYFRwhku8v/Qs9LhpuOkvEq3nNDcC1QWNGRGqqrNTBulm3knL8XwBsCO5Fy9HzadZal7pF6tq5PH+79WUmERFTdmzMwXpvDCmuPbhsi+z2N5M4ahr+Afq1KeJu9F+liMj/sG2brEWvELfmT4RalRymEUWXv0pK2m9Np4nIz9CYERH5j5ITx9k46xZSHJ+CBetD+hA1Zh6xLfVK4iLuTGNGRATYtm4lgf8YR7K9D6dtkdPpdpJu+hN+/v6m00TkF2jMiIhPs10uVr3/Er3XTyPEqqKIJhy94nVSUn5tOk1EzpLGjIj4LEfxUQpmjSO15CuwYG1oMu3GziOmeWvTaSJyDjRmRMQnFeQvI+zD8STa31Nt+5Hb5U6Sb3gUy0+XlUQ8jcaMiPgU2+VixbvPkrj5OYKsagppTvHQN0hJHGQ6TURqSGNGRHxG8dHDbMsYQ1rpN2BBflgancZl0qppC9NpInIeajxmnE4nixYtYtOmTViWRUxMDMOHDydALyglIm5oc+4Swj+6hQT7IJW2P/ndJpN03cNYfn6m00TkPNVoeaxfv55hw4ZRWFhIt27dgB/ezbp58+YsXryYnj171mqkiEhNuZwuVr7zFIkFLxNkOfneakHpsAyS4weYThORWlKjMTN+/Hi6d+9OTk4OjRs3BuDYsWOMHj2aW265hRUrVtRqpIhITRw/fJAdGaNJK18OFqxuOIDO4+fQulEz02kiUotqNGbWrFlz2pABaNy4MU8//TRJSUm1FiciUlMbs76k8b9uoQ+HqbQDyI+9j6Sr79VlJREvVKMx061bNw4ePEj37t1PO15UVMQFF1xQK2EiIjXhcjpZ+fYTJG2fQaDlZL/VipMjMkiO6286TUTqSI3GzNSpU7nzzjt5/PHHSU1NBWDlypU8+eSTTJ8+HYfDcep7f+ltu0VEasuRogPsfWs0aSdXgQV54ZfQ9ea3aBPRxHSaiNQhy7Zt+1x/yO9/TtNalgX88E6z//dzy7JwOp210VljDoeDyMhIiouLNaxEvNiGFf+m2ae30ZKjVNiBrI17kMTfT9ZlJREPdS7P3zU6M/P111/XKExEpLY5nU5WznuElF2vE2C52OPXBteVc0jqnmI6TUTqSY3GzMUXX1zbHSIi5+zQwb18/1Y6/SpywYKcyMHE3jyLsIaNTKeJSD06r1e4KysrY8+ePVRWVp52PC4u7ryiRER+ybplH9HqizuI4xjldhDrez9K0vA74D+XukXEd9RozBw6dIgxY8bwySef/OTXTd8nIyLeq7qqiqx5D5KyZzb+ls1uv2isa+aSFJNoOk1EDKnRnXGTJk3i2LFjrFy5ktDQUP7973+TmZlJly5dWLx4cW03iogAcOj7PWx6bhBpe2fhb9lkNx5Cy3tW0E5DRsSn1ejMzFdffcWHH35IUlISfn5+tG/fnssuu4yIiAimTZvGkCFDartTRHzcmqWLaPv1XfSkmDI7mE0JT5D0u9tMZ4mIG6jRmCktLaVFix/eZbZJkyYcOnSIrl270rNnT/Ly8mo1UER8W1VVJdlz7yN131z8LJud/h0IuDaThK69TaeJiJuo8SsAb9myhQ4dOtC7d2/efPNNOnTowBtvvEHr1q1ru1FEfFThvh0cyRxJWtV6sCCr6TDixr1GSFhD02ki4kZqNGYmTZrE999/D8Bjjz3G5ZdfzoIFCwgKCiIzM7NWA0XEN63+6u+0/2Yy3TlBiR1KQfJTJA8ZbzpLRNxQjV4B+P8qKytj8+bNtGvXjmbN3OvdaPUKwCKepbKiguw5d9OvcD4A2/07E3JDJm069zRcJiL1qc5fAXjKlCk/edyyLEJCQrjgggsYNmwYTZro/VBE5Owd2F1A8fxR9KveBEBW8yvpNe4VgkMaGC4TEXdWozMzl1xyCXl5eTidTrp164Zt22zduhV/f39iYmLYsmULlmWxbNkyYmNj66L7rOnMjIhnyPt8IZ2+u5dGlOAgjB19n6H35emms0TEkHN5/q7R68wMGzaMQYMGceDAAXJzc8nLy2P//v1cdtllXH/99ezfv58BAwYwefLkGv0fEBHfUVFRzvLXbqXPd7fRiBK2BnShJP0rDRkROWs1OjPTpk0bPv/88x+dddmwYQODBw9m//795OXlMXjwYA4fPlxrsTWhMzMi7mv/jk2ULBxFt+oCAFa1vJb4sX8hKDjEcJmImFbnZ2aKi4spKir60fFDhw7hcDgAaNSo0Y/es0lE5L9yP8kkPPNXdKsuoJgGrOn/Oim3zdSQEZFzVqMbgIcNG8bYsWN54YUXSEpKwrIssrKyuOeeexg+fDgAWVlZdO3atTZbRcQLnCwvY3XGHfQ9/D5YsCUwhkYj59OrnX5fiEjN1OgyU0lJCZMnT2bevHlUV1cDEBAQQHp6Oi+99BINGjQgPz8fgN69e9dm7znTZSYR97Fn6zoq3kmni3M7AKtaj6TPmBcIDAo2XCYi7uZcnr/P63VmSkpK2LFjB7Zt07lzZxo2dL9X5dSYEXEP2R/P5sKsP9LQKuc44ey9+CV6XnK16SwRcVN1/joz/9WwYUPi4uLO548QES9XXlrCmozbST36IViwKbA7zUYvoGebTqbTRMRLnNeYERE5k90Fa6h+J51U105ctkVW9BgS06cTEBhkOk1EvIjGjIjUiVUfvkGPvEdpYFVwlAgODHqF1IuGm84SES+kMSMitaqs1MG6WX8g5fjHYMGG4F60HD2fHq3bm04TES+lMSMitWbnxlx4bzQprj24bIvs9uNJHPUM/gH6VSMidUe/YUTkvNm2TdaiV+i55inCrAoO04iiwTNI6TfUdJqI+ACNGRE5LyUnjrNx1i2kOD4FC9aH9CFqzDxiW0abThMRH6ExIyI1tm39KgLfH0uyvQ+nbZHT6XaSbvoTfv7+ptNExIdozIjIObNdLla+/zLx66cSYlVRRBOOXfE6KSm/Np0mIj5IY0ZEzomj+ChbZo+j74mvwIJ1oclEj51Ht+atTaeJiI/SmBGRs7Z1zXJCPhhHkn2AatuP3C53knzDo1h+uqwkIuZozIjIL7JdLlb+7Vn6bHqeYKuKQqs5jt++SUripabTREQ0ZkTkzIqPHWHb7NH0Lf0GLMgPS6PTuExaNW1hOk1EBNCYEZEz2Jy3lPB/3kyCfZBK25/8mCkkXfsQlp+f6TQRkVM0ZkTkR2yXixV/nUpiwYsEWU6+t1pQOiyD5PgBptNERH5EY0ZETnP8SBE7MkaTVvYdWLC64UVcMH4urRs1M50mIvKTNGZE5JRNWV/S6F+30odDVNoBrOl+H4lX3avLSiLi1jRmRASX08nKt58gafsMAi0n+61WVIx4i6S4fqbTRER+kcaMiI87cuh79makk3ZyFViQF34JXW9+i4YRTUyniYicFY0ZER+2YcW/afbp7fTmCBV2IOviHiLh95N0WUlEPIrGjIgPcjqdrJr/CMk7XyfAcrHXrw3OK+eQ2D3FdJqIyDnTmBHxMYcO7uX7t9JJq8gFC3IjB3PhzbMIa9jIdJqISI1ozIj4kHXLPqLVF3cQxzHK7SA29H6ExOETwbJMp4mI1JjGjIgPqK6qImveQ6TsmYW/ZbPbLxrrmkwSYxJMp4mInDeNGREvV3RgDwfnjiStMh8syGl8BT3Gv0lIgwjTaSIitUJjRsSLrVm6iDZfT6Inxymzg9mU8ASJv7vNdJaISK3yqL9/OW3aNCzLYtKkSaZTRNxaVVUly2dNoudXY2jGcXb6d+DIjZ+RoCEjIl7IY87MZGdnM3PmTOLi4kyniLi1wn07OJI5krSq9WBBdtNh9Bz3GiFhDU2niYjUCY84M1NSUsKNN97IrFmzaNy48Rm/t6KiAofDcdqHiK9Y/dXfCZo9gO5V6ym1Q1id/AJJE+dpyIiIV/OIMTNhwgSGDBnCoEGDfvF7p02bRmRk5KmP6OjoeigUMauyooLv3riD+G/G04QTbPfvxPFRXxB/xXjTaSIidc7tLzO988475OXlkZ2dfVbf/+CDDzJlypRTnzscDg0a8WoHdm+leP5I+lVvAiCr+ZX0GvcKwSENDJeJiNQPtx4ze/fu5a677uKzzz4jJCTkrH4mODiY4ODgOi4TcQ+5ny2k8/J7iaKEE4Syve8zJF8+2nSWiEi9smzbtk1H/JwPPviA3//+9/j7+5865nQ6sSwLPz8/KioqTvvaT3E4HERGRlJcXExEhF5XQ7xDRUU5uRmTSCt6B4CtAV1oeNN8Wne40HCZiEjtOJfnb7c+M3PppZeybt26046NGTOGmJgY7r///l8cMiLeaN/OzZS+PZK06gIAslpeR/y4PxMYdHZnL0VEvI1bj5nw8HB69Ohx2rEGDRrQtGnTHx0X8QU5n2TSZeWDtLVKcdCAXf2fI3nQjaazRESMcusxIyI/OFleRl7GRNIOvwcWbAmMofGo+cRFdzWdJiJinMeNmSVLlphOEKlXe7at5+Rf00lzbgMgq/WN9BnzEgFButFdRAQ8cMyI+JKsj2YTm/1HGlrlHCecvRe/SPIl15jOEhFxKxozIm6ovKyU/Nm30ffoh2DB5qBYmqa/Tc82nUyniYi4HY0ZETezu2AN1e+k09e1E4BVbUeTOPp5/AMCDZeJiLgnjRkRN7LqwzfokfcoDawKjhJB4a/+TMqAEaazRETcmsaMiBsoK3WwdtYfSD3+MViwMbgXLUfPJ7Z1e9NpIiJuT2NGxLAdm3Lh72NIde3GZVvktB9Pwqhn8A/Qf54iImdDvy1FDLFtm1WLZhC35k+EWRUcoRFFg18hud/vTKeJiHgUjRkRA0pOFLNx1s2kOj4FCzYEx9N67DwubNnOdJqIiMfRmBGpZ9vWryLw/bEk2/tw2hY5HW8j6aY/4afLSiIiNaLfniL1xHa5WPn+y8Svn0qIVcUhmnD0ildJSbnCdJqIiEfTmBGpB47io2yZPZ6+J74EC9aHJtJ27Hy6NY8ynSYi4vE0ZkTqWMGa7wj9YDxJ9gGqbT/yLriDpBsfx/LzN50mIuIVNGZE6ojtcrHib8+RsOk5gq0qDlpNcfz2TZITLzOdJiLiVTRmROpA8bEjbJ09hrTSpWDB2rBUOoyfR8smLU2niYh4HY0ZkVq2OW8p4f+8mUT7IFW2P/ndJpF43R+x/PxMp4mIeCWNGZFa4nK6WPHOVJIKXiTIcvK91ZzS380mqc9A02kiIl5NY0akFhw/UsSOjNH0K/sOLMhv0J/ON8+ldaPmptNERLyexozIedqY/SWNP76VPhyi0g5gTew9JF59vy4riYjUE40ZkRpyOZ2sePtJkre/QqDl5IDVipO/zyCpV3/TaSIiPkVjRqQGjhz6nj0Z6fQ7uQosWB0+kC7j5xAV2cR0moiIz9GYETlH61d8SvNPbyOeI1TYgazv+QB9RkzRZSUREUM0ZkTOktPpZOW8R0jZ9ToBlot9flE4r5xDQvdU02kiIj5NY0bkLBw6uJcDb6XTryIXLMiLHETM+NmEhTc2nSYi4vM0ZkR+wdplH9P6iwn04hgn7UA29H6EhOF3gmWZThMRETRmRH5WdVUVqzIfInXvLPwtm91+0VhXzyXhwkTTaSIi8j80ZkR+QtGBPRTOHUm/ynywIKfxFfQY/yYhDSJMp4mIyP+hMSPyf+Qv/YC2X99FHMcps4PZnPg4iUNvN50lIiI/Q2NG5D+qqirJmnMffffPxc+y2eXXnoDr5tGna2/TaSIicgYaMyLA9/t2cDRzJP2q1v9wWanJUHqMf52QsHDTaSIi8gs0ZsTnrf7y77T/djLdOUGpHcLW5D+ROOQW01kiInKWNGbEZ1VWVJA95276Fc4HYId/R0Kun0/vC3oaLhMRkXOhMSM+6cDurRTPH0m/6k0AZDcfQa9xrxIUEma4TEREzpXGjPic3M8W0nn5vURRwglC2ZE6jaRfjzGdJSIiNaQxIz6joqKc3NmTSDv0DgDbAi6gwY3z6dUx1nCZiIicD40Z8Qn7dmymZOEo0qq3AJDd6lp6j/kzgcGhhstEROR8acyI18v+ZB5dVz5AW6sUBw3Y3f85kgbdaDpLRERqicaMeK2T5WXkZUwk7fB7YEFBYDcajZxPz3bdTKeJiEgt0pgRr7R723oq/ppOmnMbANmtbyR+zIsEBIUYLhMRkdqmMSNeZ9VHGcRmP0y4Vc5xGrL/4hdJuuRa01kiIlJHNGbEa5SXlbJ69u2kHf0ALNgc1J1m6fPp3qaz6TQREalDGjPiFXZtWYPz3XTSXDsByGo7moT05/APDDJcJiIidU1jRjyabdusWvwmPfMeo4F1kmNE8P2lfyb5ohGm00REpJ5ozIjHKi1xsG72H0g9/jFYsCk4jhajFxDbur3pNBERqUcaM+KRtm/MxXpvDKmu3bhsi9z240gY9Qx+AYGm00REpJ5pzIhHsW2blYtepdeaJwmzKjhCIw4NfoWkfr8znSYiIoZozIjHOOE4zsbZt9DX8SlYsCEknqgx84lpGW06TUREDNKYEY+wdV0WQf8YQ4q9D6dtkdfpDyTc+BR+AfpXWETE1+mZQNya7XKx4v2X6bN+KiFWFYdowtErXiMp5Tem00RExE1ozIjbKj5+lC2zx5NW8iVYsD40iehx8+jWLMp0moiIuBGNGXFLW/K/I/TD8STbB6i2/cjvcgcJNzyO5edvOk1ERNyMxoy4FdvlYvnfniNx03MEW1UctJri+O2bJCZeZjpNRETclMaMuI3jxw6zdfY4+pUuAQvWhqXSYfw8WjZpaTpNRETcmMaMuIVNuUuJ+OfNJHGQKtufNTGTSLj2j1h+fqbTRETEzWnMiFEup4vl70wlueBFgiwnhVZzSofNJjF+oOk0ERHxEBozYszRwwfZkTGG/uXfgQVrGvan8/hMWjVqZjpNREQ8iMaMGLEh60sa/+tWEjlEpe3Puu730ueq+3VZSUREzpnGjNQrl9PF8refIGX7KwRaTg5YragcMZuEuItMp4mIiIfSmJF6c7joe/a+lU7/k6vAgvyIgXQZN4cGkU1Mp4mIiAfTmJF6sW75v2n+2QTiOUyFHcj6uAdIGHE3WJbpNBER8XBufYPCtGnTSEpKIjw8nBYtWjB8+HC2bNliOkvOgdPpZNmch7jw0+tpxWH2+UVReM1HJFx5j4aMiIjUCrceM0uXLmXChAmsXLmSzz//nOrqagYPHkxpaanpNDkLRYX7WPfsYPrvfpUAy8XqyMtoOnkF7bunmk4TEREvYtm2bZuOOFuHDh2iRYsWLF26lAEDBpzVzzgcDiIjIykuLiYiIqKOC+W/1iz7mNZf3EELjnLSDmRT/KPED5uoszEiInJWzuX526PumSkuLgagSZOfv2G0oqKCioqKU587HI4675L/r7qqihXzHiZtz0z8LZs9ftFY18wlPibRdJqIiHgpt77M9L9s22bKlCn079+fHj16/Oz3TZs2jcjIyFMf0dHR9Vjp2w4e2M3G5wZx0d438bds8hr/hhZ3LydaQ0ZEROqQx1xmmjBhAh9//DHLli2jbdu2P/t9P3VmJjo6WpeZ6tjqJYuIXnIXzSimzA6mIPFxeg+93XSWiIh4KK+7zDRx4kQWL17MN998c8YhAxAcHExwcHA9lUlVVRUr59xPv/1v4WfZ7PJvT+B18+jdpbfpNBER8RFuPWZs22bixIksWrSIJUuW0LFjR9NJ8j/279nBsfmjuKhqHViQ23QoPca/TnBouOk0ERHxIW49ZiZMmMDChQv58MMPCQ8Pp7CwEIDIyEhCQ0MN1/m2nC//Tqdvp9AGB6WEsC35KRKuuNl0loiI+CC3vmfG+pm/xjtnzhxGjx59Vn+G/mp27aqorCD7rbvpXzgfgB3+nQi9YR6tO/c0XCYiIt7Ea+6ZceOd5ZP27SrAsSCd/tUbAchtPoKe414lKCTMcJmIiPgytx4z4j6yP/0rXZbfQ1urhBOEsqvvNBIuH2M6S0RERGNGzuzkyZPkZNxF/0PvgAXbA7rQ8Kb59Oxwoek0ERERQGNGzmDP9k2ULUynv/OHN/fMaXUNvcb8hcBg3XwtIiLuQ2NGftKqf2USs+pB2lmlOGjAnoueJ/HSG0xniYiI/IjGjJymvKyM3IyJ9D/yHliwNTCGRqPm0yO6q+k0ERGRn6QxI6fsKlhH5buj6e/cBkBO1I3Ej3kJ/0C9orKIiLgvjRkBYMXi2fTI/SPhVjnHaciBgS+ROPAa01kiIiK/SGPGx5WWlpA/+3b6HfsQLNgS1J1m6QuIbdPJdJqIiMhZ0ZjxYds352P/bTT9XDsByI0eTe9Rz+EfGGS4TERE5OxpzPgg27ZZ8eGb9Fr9GA2skxwjgoOX/pmEi0aYThMRETlnGjM+5sSJYtbNvo204o/Bgs3BPWk5egExrTuYThMREakRjRkfsnV9Dv7vjyHN3oPLtsjrMI4+I5/BLyDQdJqIiEiNacz4ANu2+e79V+iz7inCrAqO0Igjv36FxL6/M50mIiJy3jRmvFxx8XE2zr6F/ic+BQs2hcQTNXY+XVtEm04TERGpFRozXmzzmlUEfzCWvvY+nLZFfufb6HPjn7D89Y9dRES8h57VvJDtcvHt318ieeM0QqwqDluNOX7FGyQk/9p0moiISK3TmPEyx48dZfPs8Qwo/RIs2BCWRPTYeVzQLMp0moiISJ3QmPEiG1cvo8Him0m1D1Bt+7G220Tir3sMy8/fdJqIiEid0ZjxAi6ni2/feY7UgucItqoosppS+rs36dPnMtNpIiIidU5jxsMdPXKIrRnjuLhsKViwvkFfOoyfR4vGLUyniYiI1AuNGQ+2LnsJjT++hRQOUmX7syF2Mr2ufhjLz890moiISL3RmPFATqeLb99+mr7bXybYqqbQakHF72fRu9dA02kiIiL1TmPGwxw6VMiujLEMPPkdWLAu/CI6j59LWGQz02kiIiJGaMx4kPwVX9D809tIoohK25+NPe6j91X3g2WZThMRETFGY8YDVFc7WbbgCfrtnEGg5eSAXyucI96id49+ptNERESM05hxcwcLD7BvzmgGVqwCC9ZGXkLX8XMICW9sOk1ERMQtaMy4sdxl/ybqiwkkcJgKO5AtvR8kbvgUXVYSERH5HxozbqiqupplmY9w0Z43CLBc7PeLgqvnEndhiuk0ERERt6Mx42YO7N9DYeZoLqnM/eGyUuPL6DpuNiENG5lOExERcUsaM24ka8k/ab9kIn04xkkC2ZbwKHG/najLSiIiImegMeMGKior+W7OQ1x8YDb+ls1e/2gCr5tHjy59TKeJiIi4PY0Zw/bs2cWR+en8qir/h8tKza4gZuxMgsLCTaeJiIh4BI0Zg5Z/8Q+6fjuJdlYx5QSzM+UJ4n5zm+ksERERj6IxY8DJikq+y7iPSw7Oxc+y2RPQgdAbMont1Nt0moiIiMfRmKlnO3duxfH2GC6tXvfDZaUWw4gd+xoBIQ1Np4mIiHgkjZl69O0n7xK78h46Wg5KCWFf2jTiBo81nSUiIuLRNGbqQWn5SVbMnsKgI2+DBbsDO9HwpgV0a9/ddJqIiIjH05ipY1u3bqH8ndEMcm4EYF3rq4gdMwP/oFDDZSIiIt5BY6aO2LbNko/epnfOAzS2TlBCGN8PmE7PX40ynSYiIuJVNGbqgKO0jFWzJ3HZsXfBgp1BXWg86m26tO1mOk1ERMTraMzUsk2b1uP8+1guc20BYH3b64kd9TJ+QSGGy0RERLyTxkwtsW2bLxfNIXHNH2lkleKgAYd/9QI9BlxvOk1ERMSraczUguMnSsiZNZFBjn/8cFkpOIamoxfSqXVn02kiIiJeT2PmPK1bv4aA98cyyN4GwIYOo4i96XmsgGDDZSIiIr5BY6aGXC6bz9+fRdr6Rwm3yikmnOOD/0z3tCtNp4mIiPgUjZkaWvLGnVxeNA8s2B7ag1ZjF9K+eXvTWSIiIj7Hz3SAp2ofNwCXbbGx0zg63f01DTRkREREjNCZmRrq3P9qTnToTmzbWNMpIiIiPk1nZs5DuIaMiIiIcRozIiIi4tE0ZkRERMSjacyIiIiIR9OYEREREY+mMSMiIiIeTWNGREREPJrGjIiIiHg0jRkRERHxaBozIiIi4tE0ZkRERMSjacyIiIiIR9OYEREREY+mMSMiIiIeLcB0QF2zbRsAh8NhuERERETO1n+ft//7PH4mXj9mTpw4AUB0dLThEhERETlXJ06cIDIy8ozfY9lnM3k8mMvl4sCBA4SHh2NZlukcn+FwOIiOjmbv3r1ERESYzvEpeuzN0WNvjh57M+rycbdtmxMnThAVFYWf35nvivH6MzN+fn60bdvWdIbPioiI0C8WQ/TYm6PH3hw99mbU1eP+S2dk/ks3AIuIiIhH05gRERERj6YxI3UiODiYxx57jODgYNMpPkePvTl67M3RY2+GuzzuXn8DsIiIiHg3nZkRERERj6YxIyIiIh5NY0ZEREQ8msaMiIiIeDSNGak106ZNIykpifDwcFq0aMHw4cPZsmWL6SyfNG3aNCzLYtKkSaZTfML+/fu56aabaNq0KWFhYfTu3Zvc3FzTWV6vurqaP/7xj3Ts2JHQ0FA6derEk08+icvlMp3mdb755huGDh1KVFQUlmXxwQcfnPZ127Z5/PHHiYqKIjQ0lIEDB7Jhw4Z669OYkVqzdOlSJkyYwMqVK/n888+prq5m8ODBlJaWmk7zKdnZ2cycOZO4uDjTKT7h2LFj9OvXj8DAQD755BM2btzICy+8QKNGjUyneb3p06fzxhtvMGPGDDZt2sSzzz7Lc889xyuvvGI6zeuUlpbSq1cvZsyY8ZNff/bZZ3nxxReZMWMG2dnZtGrVissuu+zU+yPWNf3VbKkzhw4dokWLFixdupQBAwaYzvEJJSUl9OnTh9dee42nnnqK3r178/LLL5vO8moPPPAA3333Hd9++63pFJ/z29/+lpYtW5KRkXHq2JVXXklYWBjz5883WObdLMti0aJFDB8+HPjhrExUVBSTJk3i/vvvB6CiooKWLVsyffp0br311jpv0pkZqTPFxcUANGnSxHCJ75gwYQJDhgxh0KBBplN8xuLFi0lMTOTqq6+mRYsWxMfHM2vWLNNZPqF///58+eWXFBQUALBmzRqWLVvGFVdcYbjMt+zcuZPCwkIGDx586lhwcDAXX3wxy5cvr5cGr3+jSTHDtm2mTJlC//796dGjh+kcn/DOO++Ql5dHdna26RSfsmPHDl5//XWmTJnCQw89RFZWFnfeeSfBwcGMGjXKdJ5Xu//++ykuLiYmJgZ/f3+cTidPP/00119/vek0n1JYWAhAy5YtTzvesmVLdu/eXS8NGjNSJ+644w7Wrl3LsmXLTKf4hL1793LXXXfx2WefERISYjrHp7hcLhITE5k6dSoA8fHxbNiwgddff11jpo69++67LFiwgIULF9K9e3fy8/OZNGkSUVFRpKenm87zOZZlnfa5bds/OlZXNGak1k2cOJHFixfzzTff0LZtW9M5PiE3N5eioiISEhJOHXM6nXzzzTfMmDGDiooK/P39DRZ6r9atWxMbG3vasQsvvJD333/fUJHvuPfee3nggQe47rrrAOjZsye7d+9m2rRpGjP1qFWrVsAPZ2hat2596nhRUdGPztbUFd0zI7XGtm3uuOMO/vGPf/DVV1/RsWNH00k+49JLL2XdunXk5+ef+khMTOTGG28kPz9fQ6YO9evX70cvQVBQUED79u0NFfmOsrIy/PxOfxrz9/fXX82uZx07dqRVq1Z8/vnnp45VVlaydOlS0tLS6qVBZ2ak1kyYMIGFCxfy4YcfEh4efuo6amRkJKGhoYbrvFt4ePiP7k1q0KABTZs21T1LdWzy5MmkpaUxdepUrrnmGrKyspg5cyYzZ840neb1hg4dytNPP027du3o3r07q1ev5sUXX2Ts2LGm07xOSUkJ27ZtO/X5zp07yc/Pp0mTJrRr145JkyYxdepUunTpQpcuXZg6dSphYWHccMMN9RNoi9QS4Cc/5syZYzrNJ1188cX2XXfdZTrDJ/zzn/+0e/ToYQcHB9sxMTH2zJkzTSf5BIfDYd911112u3bt7JCQELtTp072ww8/bFdUVJhO8zpff/31T/5+T09Pt23btl0ul/3YY4/ZrVq1soODg+0BAwbY69atq7c+vc6MiIiIeDTdMyMiIiIeTWNGREREPJrGjIiIiHg0jRkRERHxaBozIiIi4tE0ZkRERMSjacyIiIiIR9OYEREREY+mMSMiIiIeTWNGRNzewIEDmTRp0o+Of/DBB1iWVf9BIuJWNGZERETEo2nMiIhXWLNmDZdccgnh4eFERESQkJBATk6O6SwRqQcBpgNERGrDjTfeSHx8PK+//jr+/v7k5+cTGBhoOktE6oHGjIh4hT179nDvvfcSExMDQJcuXQwXiUh90WUmEfEKU6ZMYfz48QwaNIhnnnmG7du3m04SkXqiMSMibi8iIoLi4uIfHT9+/DgREREAPP7442zYsIEhQ4bw1VdfERsby6JFi+o7VUQM0JgREbcXExPzkzfzZmdn061bt1Ofd+3alcmTJ/PZZ58xYsQI5syZU5+ZImKIxoyIuL3bb7+d7du3M2HCBNasWUNBQQGvvvoqGRkZ3HvvvZSXl3PHHXewZMkSdu/ezXfffUd2djYXXnih6XQRqQeWbdu26QgRkV+Sm5vLww8/zOrVqzl58iRdu3bl7rvv5rrrrqOyspL09HS+++47Dh48SLNmzRgxYgTPPfccISEhptNFpI5pzIiIiIhH02UmERER8WgaMyIiIuLRNGZERETEo2nMiIiIiEfTmBERERGPpjEjIiIiHk1jRkRERDyaxoyIiIh4NI0ZERER8WgaMyIiIuLRNGZERETEo/0/jDGtakN4da0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds.gap.plot()\n",
    "plt.plot(ds.Us, ds.Us)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06e0d356-558e-40e3-8287-d7d2e0bee8cd",
   "metadata": {},
   "source": [
    "We can also fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5499ea62-cf1b-4a13-8191-ebb73ea38704",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(0.9997852)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.gap.polyfit(dim=\"Us\", deg=1).polyfit_coefficients[0].data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "0cb395cd-84d1-49b4-89dd-da7a2d09c8d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf(\"./data/1d_hubbard_example.nc\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce428241",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "```{code-cell} ipython3\n",
    "def compute_phase_diagram(\n",
    "    Us,\n",
    "    nk,\n",
    "    nk_dense,\n",
    "    filling=2,\n",
    "):\n",
    "    gap = []\n",
    "    vals = []\n",
    "    for U in tqdm(Us):\n",
    "        # onsite interactions\n",
    "        guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))\n",
    "        full_model = Model(h_0, h_int, filling)\n",
    "        mf_sol = solver(full_model, guess, nk=nk)\n",
    "        hkfunc = transforms.tb_to_kfunc(add_tb(h_0, mf_sol))\n",
    "        ks_dense = np.linspace(0, 2 * np.pi, nk_dense, endpoint=False)\n",
    "        hkarray = np.array([hkfunc(kx) for kx in ks_dense])\n",
    "        _vals = np.linalg.eigvalsh(hkarray)\n",
    "        _gap = (utils.compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=nk_dense))\n",
    "        gap.append(_gap)\n",
    "        vals.append(_vals)\n",
    "    return np.asarray(gap, dtype=float), np.asarray(vals)\n",
    "\n",
    "import xarray as xr\n",
    "\n",
    "ds = xr.Dataset(\n",
    "    data_vars=dict(vals=([\"Us\", \"ks\", \"n\"], vals), gap=([\"Us\"], gap)),\n",
    "    coords=dict(\n",
    "        Us=Us,\n",
    "        ks=np.linspace(0, 2 * np.pi, nk_dense),\n",
    "        n=np.arange(vals.shape[-1])\n",
    "    ),\n",
    ")\n",
    "\n",
    "# Interaction strengths\n",
    "Us = np.linspace(0.5, 10, 20, endpoint=True)\n",
    "nk, nk_dense = 40, 100\n",
    "gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)\n",
    "\n",
    "ds.vals.plot.scatter(x=\"ks\", hue=\"Us\", ec=None, s=5)\n",
    "plt.axhline(0, ls=\"--\", c=\"k\")\n",
    "plt.xticks([0, np.pi, 2 * np.pi], [\"$0$\", \"$\\pi$\", \"$2\\pi$\"])\n",
    "plt.xlim(0, 2 * np.pi)\n",
    "plt.ylabel(\"$E - E_F$\")\n",
    "plt.xlabel(\"$k / a$\")\n",
    "plt.show()\n",
    "\n",
    "```\n",
    "\n",
    "The Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))\n",
    "$$\n",
    "\\epsilon_{HF}^{\\sigma}(\\mathbf{k}) = \\epsilon(\\mathbf{k}) + U \\left(\\frac{n}{2} + \\sigma m\\right)\n",
    "$$\n",
    "where $m=(\\langle n_{i\\uparrow} \\rangle - \\langle n_{i\\downarrow} \\rangle) / 2$ is the magnetization per atom and $n = \\sum_i \\langle n_i \\rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\\Delta=U$. And we can confirm it indeed follows the expected trend.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}