{ "cells": [ { "cell_type": "markdown", "id": "402ab8fc", "metadata": {}, "source": [ "# Exact diagonalization\n", "\n", "Quantax wraps [QuSpin](https://github.com/QuSpin/QuSpin) to provide a simple interface\n", "of exact diagonalization (ED). In this tutorial, we show how to perform ED to obtain eigenstates\n", "of the Heisenberg model in Quantax. You will find it very convenient and fast!\n", "\n", "Reference:\n", "\n", "[P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems part I: spin chains, SciPost Phys. 2, 003 (2017).](https://scipost.org/SciPostPhys.2.1.003)" ] }, { "cell_type": "markdown", "id": "12b77d46", "metadata": {}, "source": [ "## System definition\n", "\n", "First, we define the lattice and Hilbert space of the problem, encoded in {py:class}`~quantax.sites.Lattice`. We will use a small 4x4 lattice with periodic boundary condition (PBC). The Heisenberg model has a conserved number of up and down spins, which should also be specified." ] }, { "cell_type": "code", "execution_count": null, "id": "c83be0c7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(1, 4, 4)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAMtCAYAAABXYgSXAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYkRJREFUeJzt3Xl0VOX9x/HPTFZCFpaQkECAsMsiIAYJqIAKFEVA24po2USoChaKdaG/CqK1cUO0SkUFiUgRZW9BEQQBZRHZFJQdwp6wJxDINnN/fyDTBJKQSTIZePJ+nTPnzNz73DvfyT1PZj53ea7NsixLAAAAAGAQu7cLAAAAAIDSRtABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADCOr7cLKAqn06kjR44oJCRENpvN2+UAAAAA8BLLsnT27FlFR0fLbi/4uM11EXSOHDmimJgYb5cBAAAA4Bpx8OBB1axZs8D510XQCQkJkXTxw4SGhnq5GgAAAADekpaWppiYGFdGKMh1EXQuna4WGhpK0AEAAABw1UtaGIwAAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdEogLS1NM2bM0FNPPaUOHTqofv36CgsLk7+/vyIiItSxY0e99tprOnnypLdLBcq1jh07ymazufVYvny5t8sGyqVjx45pwYIFGj16tLp166bw8HBXvxwwYECR1uF0OvXLL78oMTFRTzzxhOLi4hQQEED/BkpJafTTxMTEIn8nJyYmFqtO32ItBUnSunXr1KdPn3znHT9+XCtWrNCKFSv0+uuva9q0aeratWsZVwigOOx2uxo0aODtMoByKTIyssTr+OSTT4r8YwuA+0qjn5YFgk4JxcTEqFOnTmrdurViYmIUFRUlp9OpQ4cOadasWZozZ45OnDihHj16aN26dWrRooW3SwbKnSlTpig9Pb3QNr/88ot69+4tSbrzzjtVo0aNsigNQCFq1aqlxo0ba/HixW4tZ1mW67mfn5+aN2+u7OxsbdmypbRLBMq94vbT3L766itFR0cXOL9mzZrFWi9BpwQ6deqkAwcOFDj/gQce0Lx583TfffcpKytLY8eO1Zw5c8qwQgCSFBsbe9U2n3zyiet5v379PFkOgEKMHj1acXFxiouLU2RkpJKSkorUh3Nr0qSJ/vnPfyouLk4tW7ZUYGCgXnjhBYIOUEpKo5/m1rBhQ9WpU6f0CvwVQacEfHx8rtqmV69eatSokXbs2KFvv/22DKoC4C6n06l///vfkqTg4GDdf//9Xq4IKL/Gjh1b4nW0adNGbdq0KYVqAOSnNPppWWAwgjIQEhIiScrIyPByJQDys3TpUh0+fFiS9Lvf/U5BQUFerggAAJQUQcfDduzYoc2bN0uSGjdu7N1iAORr6tSpruectgYAgBkIOh5w/vx57dq1S2+++aY6dOignJwcSdKIESO8WxiAK5w7d05z586VJNWuXVsdO3b0bkEAAJQzAwcOVHR0tPz9/RUeHq62bdvqb3/7m+tsi+LiGp1SkpiYqIEDBxY4/7nnntNDDz1UhhUBKIrZs2e7RmT7wx/+IJvN5uWKAAAoX3Lf2+rkyZM6efKkvv/+e40bN05vvfWW/vjHPxZrvQQdD2vZsqU++OADxcXFebsUAPngtDUAALyjbt26uv/++xUfH6+YmBhJ0t69ezV79mzNmjVLGRkZeuyxx2Sz2TRkyBC310/QKSW9evXSzTffLEm6cOGC9uzZo88//1xz585Vnz599NZbb6l79+5erhJAbocOHXLtRWrbtq0aNmzo3YIAACgn7rvvPvXv3/+KMyni4uLUu3dvLViwQPfff7+ys7P15z//WT169FD16tXdeg+u0SkllSpVUrNmzdSsWTPFxcXpwQcf1Jw5czR16lTt3btXPXv2VGJiorfLBJDLtGnT5HQ6JUn9+/f3cjUAAJQfYWFhhZ4u3r17d40ePVrSxevfJ0+e7PZ7EHQ8rG/fvvr9738vp9OpYcOG6dSpU94uCcCvLt0kNCAgQL179/ZyNQAAILchQ4a4wtCKFSvcXp5T14royJkLmrvpsA6dPi+bzabYqhV13001FB4ccNVle/bsqc8//1zp6elatGgRgxIAHrQr5az+8+MRHUvLlJ+vTU2iwtSjZbSCA/L+u1u/fr1++eUXSRf3GlWuXNkb5QLlkmVZ2njgtBb/nKLT57NUwc9HretU0W+aVpe/L/tggWuFw2lp5c7j+nbXCZ3LzFZwgJ9ubxiu2xtUk93u+cF7IiIiVLVqVZ04caJYI7ARdK4iOTVDz8/fqq+3pcgmyf5rqnRall5ZtF09WkTrhXubKizIr8B1VKtWzfV8//79ni4ZKJe2J6fp+Xlb9UPSafnYbbr07zfHeUAvLfhFf2hbS093bez6EZV7EAJOWwPKzqrdJzT2vz9rZ8o5+eb6ofTxmv2qVMFPT3Sqp8G31WUERMDL5m46pNcW7dDR1Az52m2yJNkkfbRqn6IrBerZ3zRWz5Y1PF5HSf4XEHQKcfDUef32vdU6mZ4ly5IsXQw4Lpal/2w+os0Hz2j24+1UpaJ/vuvJnUCDg4M9XDVQ/mw8cFoPf/i9snIuXm/jcFp55l/Idmjyd/u09XCaEh+Jk91yasaMGZIu7ojo1q1bmdcMlEcLfzqqJz/dqEtfpTmX9dUzF7L1jy+2a+/xdCXc35ywA3jJe8v36NVF212vL++rR85kaPiMzUpJy9CQ2+t5rI7jx4/rxIkTkqTo6Gi3l+f4cAEcTksDE3/QyfSsK3405WlnWTpw6ryGTd9YYJuZM2e6njdv3rxU6wTKu9Tz2Ro4ZZ0ycxxyWAX3Vaclrd13Un9fsE1ffvmljh8/Lkl66KGH5OvLPh/A03alnNXwGZtcOw4LM+OHg/p4dVJZlAXgMt9sP5Yn5BTmH19s14qdxz1WywcffCDr1+/2Dh06uL28W0Hnvffe04033qjQ0FCFhoYqPj5eX375ZaHLzJw5U40bN1ZgYKCaN2+uL774wu0ivWHlzuPafexcoSHnEofT0uo9J/XzkdQr5o0fP971mWNjY3XbbbeVeq1AeTZzw0GlXchREbqqLEua8cMBfTTtM9c07p0DlI0pq5Nk6eoh55KJK/YU6TsYQOn61/Ld8iniwVQfmzRx+R633yMpKUmbNm0qtM2CBQv04osvSpIqVKiggQMHuv0+bu3GrFmzpl555RU1aNBAlmXp448/Vs+ePbVp0yY1bdr0ivarV69Wnz59lJCQoO7du2v69Onq1auXNm7cqGbNmrldbFn6eE2SfOy2Iv+TtcvS63PWaFibKjp39qy2/vyzPp0xQ6vXrJEk+fv761/vvitndrac2dmeLB0oNyzLUuKqfW4tk+OwtPJQliSpadOmat6kibIzMjxRHoBfnc3I0ewNh9wKLslpmfrrPz9RlQv/O/17186dmvzhh3na9evbN9/lp/46quIlGzdscD1fuGCB9uza5Xpdr25dtW/fvsi1AabafTxdPySdLnJ7hyWt2XtSH3627n/r2L37iluqDBgwIM/rpKQkderUSfHx8br33nvVokULRURESLp4w9BZs2Zp1qxZrqM5b7zxhmrUcP96IJtlFXKuRxFUqVJFr7/+ugYNGnTFvN69eys9PV0LFixwTWvbtq1atmypiRMnFrjOzMxMZWZmul6npaUpJiZGqampCg0NLUm5RdZy7GKdueBeIMk6vl9HPxp6xfSwCoHqHXejGlavls9SAIor0+avD+pc+b+nUJZT6du/04n/vKZ7bmysTo09d24xgIuS/SM0s8Zv3VrGcuQode3nSv1ueqHt3njgnnyn/+XzhUV+r5vr1NSDbVq4VR9gou3BDbWk2p1uL3f8P6/r/LaCh3++PG4sX75cnTp1uup6g4KCNH78eA0ZMiTP9LS0NIWFhV01GxT7xHSHw6GZM2cqPT1d8fHx+bZZs2aNRo4cmWda165dNW/evELXnZCQoLFjxxa3tFKR5XC6vYyPX4DsNpv8fHwUHOiv6EqhahIVoRYx0fL39fFAlUD55rAVp1/ZZPPxk91m0021PT9aDADJYS9OX7Vk88l/kB8AnpFTrO9VyebrXl9t3bq1pk2bpjVr1mj9+vU6evSoTpw4oZycHFWuXFlNmzbVnXfeqUcffdR1pKdYdbl7RGfLli2Kj49XRkaGgoODNX36dN199935tvX399fHH3+sPn36uKb961//0tixY5WSklLge1wLR3Tav7JMh89cKHJ7m6RWMWH6bNDNnisKQB6ZOU7d+PI3Rbo+5xJfu033t4zSyz1u8FxhAPLYfTxd3SasdWsZu016+q76erR9bQ9VBeByS7Yd1xOf/eT2cpP63ay7mkR6oKL8eeyITqNGjbR582alpqZq1qxZ6t+/v1asWKEmTZqUqODcAgICFBBw9RtxelKvVtF6b/ket35A3XdTTfkFBnquKAB5+Enq3CRSX287VuRz/3Oclnq0iqGvAmWocc0A1atWUXuPpxd5MALLkrrTV4Ey1bFplIL8f9b5LEeRlwkO8NWtDcI9WFXxuT28tL+/v+rXr6/WrVsrISFBLVq00Ntvv51v2+rVq19x5CYlJUXVq1cvXrVlqE+bWnLnWJe/r129WnEaDFDW+sXXKXLIsdmkWlWC1K5eVQ9XBSA3m82mAe1jixxyfOw2dWxUTTFVgjxaF4C8gvx91TsuRj72og275mOz6cG4GAX6XZuXaJT4PjpOpzPPaWa5xcfHa+nSpXmmLVmypMBreq4lNSsHadgd9Yvc/m/dmygk0M+DFQHIT7t6VXV38+q62v9k26+Pv/dqJnsR/4EDKD2/b11TzWuEXfUHlN0mBfjaNepuTi8FvOGJjvVVLTjgqn3Vx25TRGiAHut47Q7q41bQGTVqlFauXKmkpCRt2bJFo0aN0vLly/Xwww9Lung/ilGjRrnaDx8+XIsWLdK4ceO0fft2vfDCC1q/fr2GDRtWup/CQ0Z2bqght9eVpHw39qVpf727sfq25RxiwBtsNpvefKClOv96bnBBfdXXx6Z3H7pJtzdk9EPAGwL9fPTxI23UNPri+fT5/Yay2y6eBvPJoFvUMDKkjCsEIEnVQgL06ZC2igwNcO0kzO3StKiwQM0Y0lbhwd693KQwbg1GMGjQIC1dulRHjx5VWFiYbrzxRj377LPq3LmzJKljx46qU6dOnrGzZ86cqb/97W9KSkpSgwYN9NprrxU4eEFBinrBkaf8kHRKH69O0hdbjrqu2fH3salXqxrqF19HzWqElXlNAPJyOi0t235MH69J0re7TrimBwf4qE+bWvpD29qqXbWiFysEIEmZOQ4t/OmopqxK0pbD/7vRdniwv/rH19GDbWqpWsi1+8MJKC/OZmRr9oZDmrJqn/af+t8AXXWqBmlg+1j9tnVNBQcUewDnEilqNijxfXTKgreDziWnzpzTuGGPyWZZ+su/3lelsGCv1QKgYCkn0/TP4UNltxx69v1JCg7mPH/gWnT4eKom/HmYfK0cjfpgsipUpK8C15qsCxf0j0EDlWX31x/f+KdqRYbJZvPuKeAev49OeRQS6Ksq2WckSRW9lGABXF2Viv6qkn3xzs4BviW+FBGAh0SEBKjqr33V14e+ClyLbDabQhznJIcUXSnQ6yHHHfxXAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACM41bQSUhIUFxcnEJCQhQREaFevXppx44dhS6TmJgom82W5xEYGFiiogEAAACgMG4FnRUrVmjo0KFau3atlixZouzsbHXp0kXp6emFLhcaGqqjR4+6Hvv37y9R0QAAAABQGF93Gi9atCjP68TEREVERGjDhg26/fbbC1zOZrOpevXqxasQAAAAANxUomt0UlNTJUlVqlQptN25c+dUu3ZtxcTEqGfPnvr5558LbZ+Zmam0tLQ8DwAAAAAoqmIHHafTqREjRqh9+/Zq1qxZge0aNWqkjz76SPPnz9e0adPkdDrVrl07HTp0qMBlEhISFBYW5nrExMQUt0wAAAAA5VCxg87QoUO1detWzZgxo9B28fHx6tevn1q2bKkOHTpozpw5qlatmt5///0Clxk1apRSU1Ndj4MHDxa3TAAAAADlkFvX6FwybNgwLViwQCtXrlTNmjXdWtbPz0+tWrXS7t27C2wTEBCggICA4pQGAAAAAO4d0bEsS8OGDdPcuXO1bNkyxcbGuv2GDodDW7ZsUVRUlNvLAgAAAEBRuHVEZ+jQoZo+fbrmz5+vkJAQJScnS5LCwsJUoUIFSVK/fv1Uo0YNJSQkSJJefPFFtW3bVvXr19eZM2f0+uuva//+/Xr00UdL+aMAAAAAwEVuBZ333ntPktSxY8c806dMmaIBAwZIkg4cOCC7/X8Hik6fPq3BgwcrOTlZlStXVuvWrbV69Wo1adKkZJUDAAAAQAHcCjqWZV21zfLly/O8Hj9+vMaPH+9WUQAAAABQEiW6jw4AAAAAXIsIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoltH79er344ovq0qWLatasqYCAAAUHB6thw4YaOHCgvvvuO2+XCCAfBw4c0JgxY3TzzTerWrVqCgwMVExMjG677TaNHj1aW7du9XaJQLl07NgxLViwQKNHj1a3bt0UHh4um80mm82mAQMGuL2+L7/8Uvfdd5/rO7pmzZq677779OWXX5Z+8UA5Udr9NLfz58+rbt26rvXVqVOn2OvyLVEl5dztt9+ub7/99orpWVlZ2rVrl3bt2qXExET169dPH374ofz9/b1QJYDLvfPOOxo1apTS09PzTD906JAOHTqk7777TmlpaXrrrbe8UyBQjkVGRpbKepxOp4YMGaLJkyfnmX748GEdPnxY8+bN06OPPqr3339fdjv7fQF3lFY/zc/o0aO1b9++UlkXQacEjhw5IkmKjo7W73//e912222qVauWHA6H1qxZo3Hjxunw4cOaOnWqsrOzNX36dC9XDODvf/+7nn/+eUlSw4YNNXjwYMXFxSksLEwnT57Upk2bNHfuXH74ANeAWrVqqXHjxlq8eLHby/7f//2fK+S0atVKzzzzjOrVq6c9e/botdde06ZNmzRp0iRVq1ZN//jHP0q7dKDcKEk/vdymTZv01ltvKTAwUH5+fjp79myJ1kfQKYHGjRvrH//4h37729/Kx8cnz7y2bduqb9++at++vXbu3KlPP/1Ujz32mG6//XYvVQtg6dKlrpDTr18/TZo0SX5+fnna3HnnnfrLX/6irKwsb5QIlHujR49WXFyc4uLiFBkZqaSkJMXGxrq1jp07d+qNN96QJN18881auXKlKlSoIEmKi4tTjx491KFDB61fv16vv/66HnnkEdWvX7/UPwtgqtLop5dzOBwaPHiwHA6HxowZo8mTJ5c46LDLsgQWLFigBx544IqQc0l4eLjGjRvnej1r1qyyKg3AZZxOpx5//HFJUosWLTR58uQrQk5unGoKeMfYsWPVvXv3Ep0a89ZbbyknJ0fSxVNVL4WcS4KCgvTOO+9IknJycjR+/PjiFwyUQ6XRTy/39ttva8OGDWrUqJGeffbZUlknQcfDOnXq5Hq+Z88eL1YClG+LFy/Wrl27JEnPPvusfH05oA2YyLIszZ8/X9LFMy/atm2bb7u2bduqUaNGkqT58+fLsqwyqxFAXvv379fo0aMlSRMnTiy1nY0EHQ/LzMx0PS/oyA8Az5s5c6YkyWazqXv37q7pp06d0q5du3Tq1ClvlQagFO3bt891DW2HDh0KbXtp/uHDh5WUlOTp0gAU4IknnlB6err69u2rjh07ltp6CToetmLFCtfzG264wYuVAOXb2rVrJUl16tRRSEiIpk+frubNm6tq1apq2LChqlatqkaNGumNN97Is4MCwPXll19+cT1v3LhxoW1zz9+2bZvHagJQsBkzZuiLL75Q5cqV81zyURoIOh7kdDr1yiuvuF4/8MADXqwGKL+cTqe2b98u6eK1c8OHD9fDDz98xb1ydu7cqaefflp33HGHzpw544VKAZTUoUOHXM9r1qxZaNuYmBjX84MHD3qsJgD5O336tEaMGCFJeuWVV1StWrVSXT9Bx4PGjx+vdevWSZLuv/9+tW7d2ssVAeVTamqqnE6nJGnLli365z//qaioKE2bNk2nTp3S+fPntWLFCte5/KtXr9YjjzzizZIBFFPuUZqCg4MLbVuxYkXX83PnznmsJgD5e/rpp5WSkqL4+HgNHjy41NdP0PGQFStW6LnnnpMkRURE6L333vNyRUD5lfvGoBkZGQoKCtI333yjhx9+WJUrV1aFChV0++23a9myZWrRooUkae7cufr++++9VTKAYsrIyHA9v9oFzQEBAa7nFy5c8FhNAK60cuVKffTRR/L19dXEiRNls9lK/T0IOh7w888/67777lNOTo4CAwM1c+ZMRUREeLssoNwKDAzM8/rRRx91jbaUW4UKFfTyyy+7Xn/22Wcerw1A6crd3692P6zc1+NdPgQ1AM/JzMzUkCFDZFmWhg8frhtvvNEj70PQKWX79u1Tly5ddPr0afn4+GjGjBncJBTwspCQkDyvu3TpUmDbO++80zX09A8//ODRugCUvtz9/Wqno+U+2nu109wAlJ6XX35ZO3bsUExMjMaOHeux9+FGEm5IOnleBwJryiZLh89cUJ3qefcSHzlyRHfddZeOHDkim82mjz76SD179vRStUD5ZFmWtief1f4KMfKxHDpxLlNR4WGqVq2ajh8/LinvBciXCwwMVHh4uJKTk13tAZQ+h9PST4fTtL9CjHytHJ3NyFGVwKsvdzW5ByDIPTBBfnIPQFDY/wWgPMvMduhoQKSy7P7aeOCMWsRWU6BfyW6Z8uqrr0qS7rrrLv33v//Nt82lHRHp6emaMWOGpIuXg9xxxx1Ffh+CzlU4nZb++9MRJa5K0qaDZ6SoeyVJ895arfb1q+qR9rG6o3GETp48qc6dO2vv3r2SLt6JuV+/fl6sHChfMnMc+vyHg5qyOkl7j6dL1S/eK2f+uO/0m2bVVe+Wu3R8waeSJIfDUei6Ls3npqJA6Uu9kK1pa/dr6uokpZzNdPXVBW98q9/eVFODbq2j+hEhV1lLwZo0aeJ6fmm0xYLkns8tIIC8jqZeUOKqJE3//oDORt8vSfrPRxsUGuirPrfU0sB2saoeVry9E5dOK50yZYqmTJlSaNsTJ06oT58+ki7e+4qgU0qyHU79+bPNWvDTUdnzuT5q7Z5TWrX7pB66OUqLEga7xu5/5ZVXNHTo0DKuFii/Ui9ka+CUddp04MwV85yW9NXPKXI0fUghR87q7MYF2rt3r1q1apXvutLS0nTixAlJUo0aNTxZNlDuHDx1Xg99uFaHz1yQ08o7LyvHqc/XH9TsDYc04eGb1LlJZLHeIzY2VtHR0Tpy5Eiee9nlZ+XKlZIu9vU6deoU6/0AE/148Iz6Tv5e6VkOOS7rrGkZOZq0cp9mrDuoaYNuUfOaYV6q8uq4RqcQo+dt1cItRyXpin/IkuSwLk6cvv6o9la4eGHz//3f/+nZZ58tsxqB8s7htDRk6nr9eDBVlqR8uuqv/6RtqtL5MQXd0EFz584tcH1z586V9Wvfvu222zxSM1AepWVk6+FJ3+tIaka+36nSxb6a7XDq8WkbtGH/qWK9j81mc502vn37dtfNgi+3du1a1xGdnj17emTEJ+B6dODkefWd/L3OZeZcEXIucViWzmZk6+FJa3Xw1Hm338OyrKs+ateuLUmqXbu2a9ry5cvdeh+3gk5CQoLi4uIUEhKiiIgI9erVSzt27LjqcjNnzlTjxo0VGBio5s2b64svvnCrSG/YfeysPv3hoKwC/hlfLiy+tx4b/hf9/e9/92xhAPJYtv2Yvt93yrXjoVCWpcp3DNKnMz7T0qVLr5idnJysv/3tb5IuDks7cODA0i4XKLc+/f6ADp4+X+APp0ssSU7L0itfFn7aWWFGjBghH5+L1xA8+eSTVwwdfeHCBT355JOSLp6ieumGhQCkd7/ZpfQsR4E7JC5xWlJ6lkPvrdhTNoUVg1unrq1YsUJDhw5VXFyccnJy9Ne//lVdunTRL7/8kuemW7mtXr1affr0UUJCgrp3767p06erV69e2rhxo5o1a1YqH8ITpq09IB+77ar/kC+x2e2KbNvrijut5+bv76+GDRuWVokAJH28Okk+NlvRgo7NJt/gKgqoe7O6d++uESNG6O6771aFChW0bt06JSQkuC5efumllzh1DSglTqelxNVJRd556LSkH5JO67X3P5H9bIpr+u7du5WYmJin7YABA65YvmHDhnr66af1yiuvaP369Wrfvr2effZZ1atXT3v27NGrr76qTZs2Sbp4w8IGDRoU96MBRkk9n615m44U+fevw2np83X7VXHXTte0ovbTsmCzrKL+27nS8ePHFRERoRUrVhQ4hHLv3r2Vnp6uBQsWuKa1bdtWLVu21MSJE4v0PmlpaQoLC1NqaqpCQ0OLW65bbvnH10pJy7x6w19ZllNZR3cp+ZOnCmxTu1Yt7SrCETAARZOR7VDzl5e7tYzdJmXt/E6H57yS73ybzabnnn1WY8eMKYUKAUjS9uSzunfiOreWsZwOnVn5idK+n1Vou6wCbvTpdDr12BNPKPHjjwtcduCAAXpvwgTZ7ZzJD0jSlz+n6E8zC95pX5Bjc17WhV1rCpxfnLhRp04d7d+/X7Vr11ZSUlKeeUXNBiUajCA1NVWSVKVKlQLbrFmzRiNHjswzrWvXrpo3b16By2RmZua5iVdaWlpJyiyWtAs5brW32eyyVyg8hKWdOK5/9v9dScoCkEu6TwWp1gC3lnE6LdWtHqUmTRro5yMpOpV+XjlOp0IDA1SvWlXd2qCOqu7dQl8FStHhwCgpqpd7C1lO2StcffS1wvpqM0mDbovT2j0HdPD0GaVnZqtigJ9iKldS23q1dMP543p34APu1QUYbGvIDVJ4R7eX86lwbd6HqthBx+l0asSIEWrfvn2hp6AlJycrMjLvyCmRkZFKTk4ucJmEhASP3jyoKCr4+ehCduFD0F6uekU/jXzgHg9VBOByfk73dkhIkk2WAm1OdW3WUF2bcSopUBZ8i9FX7XYfdW5YS7dElOx79YaoCN0QFVGidQDlRXG+VyXp31MT1aNFdKnWcvlRnOIodtAZOnSotm7dqu+++67ERVxu1KhReY4CpaWllfmNvNrWq3JxSNoinqPoY7Opx+2t9Kff9PZwZQByWzZhrXYfT893tLX8WDa7Hu7dXQ/FPebRugD8T0a2Q1+88a3OZRZ9B6Jls+vJEUPUNvY5D1YGILdDpy9oyduri/ydKkk2m9S6dmWP1VQSxQo6w4YN04IFC7Ry5co8dyDOT/Xq1ZWSkpJnWkpKiqpXr17gMgEBAQoICChOaaWmX3wdfbGl4KNOl3NYlvq2ryu/wFK4rTOAIht4a139de6WIrev4Oej37WpI78AbiMGlBW/QKlPm1r6aFVSkXYg2iTVrhqkWxtHMewzUIZiowLVsVE1rdx1okh91cduU8dG1VSjUoUyqM59bl19Z1mWhg0bprlz52rZsmWKjY296jLx8fFXDOO6ZMkSxcfHu1dpGbsltopuqlVJPvndKfQydpt0d/Mo1a12bZ6fCJisV6toRYUFFqmvStLg2+uqIiEHKHP94usowNee7w24L2dJGnFXQ0IO4AVDO9Uv8uABlmXpiY71PVxR8bkVdIYOHapp06Zp+vTpCgkJUXJyspKTk/OMT9+vXz+NGjXK9Xr48OFatGiRxo0bp+3bt+uFF17Q+vXrNWzYsNL7FB5gs9k0qX+cYsMryqeQ/7N2m9SqVmW98fsby644AC5B/r76ZNAtqlTB76php2fLaI24k2FkAW+IqRKkSf1vlp+PXT4FBJhLU0fc1UC9WjG8O+ANN9epotd/10I2mwrcMWH/dd64B1pcs6etSW4OL13QnpUpU6a4xsfu2LGj6tSpk2f87JkzZ+pvf/ubkpKS1KBBA7322mu6++67i1ykN4aXdr13RrbeXLxTn60/qAtZDtkshySbLJtdoYG+6htfW0/e0UCBfj5lWheAvI6cuaDXv9qh//54cfx/m+WQ9WtfrR4WqCG31dWAdnVkL+KRHwCe8cuRNL2+eLuWbz8u2SSb8399tUFEsJ68s0GpX9QMwH2rdp/QuMU7tPHAGdksp2yyJLuPnJZ0c+3KGtmlodrVC/dKbUXNBiW6j05Z8WbQuSQ9M0cLNx3U7CmJslmW+gwZpN+0rKkAXwIOcC05eS5TCzcf1IJPpsnHcuiRYX/UHc1qFPnUNgBl4+Cp8/rqp8Na/Om/5Wvl6Mmnhqltw0hOVwOuMVv3n1DC38cry+6vO3s/qE5NotWo+tWHfvekMrmPTnlSMcBX97WM0sEzGyRJdzf7q/wIOcA1p2pwgPrcXFPH31kvSerYMJyQA1yDYqoEqX/bGKW+d7Gv3ly7EiEHuAY1igzWTWk/SpIGtXv2uhp4i1sBAwAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMI7bQWflypW69957FR0dLZvNpnnz5hXafvny5bLZbFc8kpOTi1szAAAAABTK7aCTnp6uFi1aaMKECW4tt2PHDh09etT1iIiIcPetAQAAAKBIfN1doFu3burWrZvbbxQREaFKlSoVqW1mZqYyMzNdr9PS0tx+PwAAAADlV5ldo9OyZUtFRUWpc+fOWrVqVaFtExISFBYW5nrExMSUUZUAAAAATODxoBMVFaWJEydq9uzZmj17tmJiYtSxY0dt3LixwGVGjRql1NRU1+PgwYOeLhMAAACAQdw+dc1djRo1UqNGjVyv27Vrpz179mj8+PH65JNP8l0mICBAAQEBni4NAAAAgKG8Mrx0mzZttHv3bm+8NQAAAIBywCtBZ/PmzYqKivLGWwMAAAAoB9w+de3cuXN5jsbs27dPmzdvVpUqVVSrVi2NGjVKhw8f1tSpUyVJb731lmJjY9W0aVNlZGRo0qRJWrZsmRYvXlx6nwIAAAAAcnE76Kxfv16dOnVyvR45cqQkqX///kpMTNTRo0d14MAB1/ysrCw99dRTOnz4sIKCgnTjjTfq66+/zrMOAAAAAChNbgedjh07yrKsAucnJibmef3MM8/omWeecbswAAAAACgur1yjAwAAAACeRNABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQKaFjx45pwYIFGj16tLp166bw8HDZbDbZbDYNGDDA2+UB+NWlfnm1R8eOHb1dKlDuZWRk6F//+pfuvPNOVatWTf7+/oqOjtbdd9+tGTNmeLs8oNwrjd+/27Zt07vvvqv+/fvrpptuUs2aNRUYGKiKFSuqbt266t27t+bPny/Lsopdp2+xl4QkKTIy0tslAABgjB07dqhnz57asWNHnulHjx7V0aNH9eWXX2rKlCmaPXu2goODvVQlUL6Vxu/fl19+Wf/+97/znbdv3z7t27dPn3/+uTp06KDZs2eratWqbr8HQacU1apVS40bN9bixYu9XQqAAjz++ON64oknCpxfsWLFMqwGQG7Hjh1T586ddfDgQUnS73//e/Xv31/R0dE6cuSIPv74Y82cOVOLFy/Wgw8+qAULFni5YgDF/f3r6+urW265Re3bt1fz5s1VvXp1VatWTadPn9b27dv1/vvva+vWrVqxYoXuvfdefffdd7Lb3TsZjaBTQqNHj1ZcXJzi4uIUGRmppKQkxcbGerssAAWIiIhQs2bNvF0GgHy8+OKLrpAzZswYvfDCC655rVq10j333KMxY8boxRdf1MKFCzVr1iz97ne/81K1QPlVGr9/J02aJF/f/KPIXXfdpccff1wPPPCA5syZozVr1mjBggXq0aOHW+/BNTolNHbsWHXv3p1T2AAAKAGHw6Fp06ZJkmrXrq3nn38+33ajR49WrVq1JEmvvPJKmdUH4H9K4/dvQSHnEh8fHz399NOu199++63b70HQAQAAXrdr926lpqZKkjp37iwfH5982/n4+Khz586SpA0bNmjfvn1lViOAshUSEuJ6npGR4fbyBB0AAOB1p06edD2/2l7i3POLs5cXwPUh9yiLjRs3dnt5gg6AcmXmzJlq0qSJgoKCFBISogYNGqh///765ptvvF0aUK5VzDWC2qUjOwXJPf+XX37xWE0Ayt6JEye0Zs0aDRo0SC+//LIkKTw8XA8//LDb62IwAgDlyuU/inbv3q3du3dr6tSp6tWrlxITExUWFual6oDyq369evLz81N2drZWrlxZaNvc8w8cOODp0gB4WMeOHbVixYp854WHh2vu3LmqVKmS2+vliA6AciEoKEgPPvigPvzwQ3377bfatGmTFi9erP/7v/9zjc0/b9489ezZU9nZ2V6uFih/KlasqDvuuEOS9NNPP+nTTz/Nt92nn36qLVu2uF6fPXu2TOoDUPb+9Kc/adu2bbr11luLtTxHdACUC4cPH853b1Dnzp315JNPqlu3btq0aZNWrFih9957T3/605/KvkignHvhhRe0dOlS5eTkqH///tqzZ4/69eunqKgoHT16VFOnTtWLL74of39/ZWVlSZIuXLjg5aoBlNSUKVOUnp4uy7J05swZrV+/Xu+9957effdd7d27V5MmTSrWCG8c0QFQLhR2yDsyMlKzZs2Sn5+fJOmdd94po6oA5Na2bVu9//778vX1VXZ2tp5//nnVrl1b/v7+riGnfX199eabb7qWyT0qE4DrU2xsrJo1a6bmzZvrtttu05///Gf99NNPuvvuu7VgwQLFxcXp0KFDbq+XoAMAkurWresasnb37t06cuSIlysCyqdHHnlE33//ve677z5VrFjRNd3X11c9evTQxo0bdfPNN7umV65c2RtlAvCwwMBATZkyRUFBQTp48KCeeeYZt9fBqWtFkHohW7M3HNKn3+9XUu1BkqRl769T3/g66tEyWkH+/BmBa0FKWoY+XXdAs9cf1NHaj8puObUucYP6t6+rzk0i5edT+L6dJk2a6IsvvpB08VS36OjosigbKHf2Hj+naWsPaMFPh3Wq9qPysRza9umPGnBrPbWvX1U33XST5syZo5ycHB09elRZWVmqUaOGAgMDJcl1Y1FJatq0qbc+BmC8Hw+e0cer9ujLWv2VbfPXv19bqTtviFTf+Nq6sWYlj79/eHi42rdvryVLlmj+/PnKzs52nX1RFPxCv4ovthzVyM82KzPHKUmy7P6SpF+OntVzc7bo5S+26f0/tFa7+uHeLBMo1yzL0off7tWrX+6QJUtOS5L94j/CH/af0fdJGxVdKVCJA9uoYWTBp7nYbLYyqhgon3IcTo397y/6ZO1++dhtcjgtye6nHPnpm10n9fWOE2oWHaqPBsYpIiRQvr6+iomJuWI9GzZscD1v06ZNWX4EoFxIz8zRk59u0rLtxy72VZ8gSdKp89mas+mwZm44pLtuiNA/+7Ty+A7/atWqSZLOnz+vEydOKCoqqsjLcupaIRb+dFRD/71RmTlOWZKsXPMuPU/PzFG/j9Zp7d6T+awBQFn41/I9+scX2+Wwfg05uVx6nZKWqd+9t1p7j58rcD25h57maA5QuizL0rOzf9K0tfsl6WLIyeXS623JZ/W799bodHpWvutxOByaM2eOJCkmJkbt2rXzYNVA+ZOZ49CAKeu0fMcxSQX31WXbj2nglB+U9evBAE85fPiw63lwrvttFYXbQWflypW69957FR0dLZvNpnnz5l11meXLl+umm25SQECA6tevr8TERHfftsylXsjWU59vlpQ34FzOaUlOy9LQ6RuV7fDshgZwpR3JZ/X6Vzuu2s7htJSe5dBfZv6Y7/x9+/ZpyZIlkqR69eqpRo0apVonUN4t/iVFszceLvQ7VbrYVw+fuaBXFm3Pd/7kyZNd98754x//KB8fn1KuFCjfPvouSRv2n75ix+HlnJa0LumUpqza57FaDh06pDVr1kiSateu7fbgI24HnfT0dLVo0UITJkwoUvt9+/bpnnvuUadOnbR582aNGDFCjz76qL766it337pMzdl4yHUk52qclnTyXJa+/iXF43UByGvar6fAFIXDaWnjgTPanpyWZ3pKSop++9vfuoarfeKJJ0q9TqC8S1yVpCJ2VTmcluZuPKzUC3nvabVs2TKNGDFCktSwYUM99dRTpVwlUL45nJY+Xr3vqiHnEsuSElcnyWkVcYFf7dy5U8uWLSu0TWpqqh566CHXd3O/fv3ceg+pGNfodOvWTd26dSty+4kTJyo2Nlbjxo2TJN1www367rvvNH78eHXt2tXdty8zn65z707LNll6Y85qdfHf6Zq2a+dOTf7wwzzt+vXtWyr1Abh4vv+sDQevOKxeGLssvTV3jQbdVEknT57UipUrNWnyZJ04cUKS1L5dOw0ZNEjZGRmeKhsod46cydAaN0/xznI4NH7mN+pWP1gHDh7U/P/8R5/OmCGn06kqVaro3598Ih+JvgqUotV7Tyk5LdOtZY6mZuifn37per179+4rzt4aMGBAntdHjhzRnXfeqRYtWqhXr15q3bq1qlevLl9fXyUnJ2vVqlWaPHmykpOTJUnNmjXTc8895/bnsVmWmxEs98I2m+bOnatevXoV2Ob222/XTTfdpLfeess1bcqUKRoxYoRSU1PzXSYzM1OZmf/7I6elpSkmJkapqakKDQ0tbrluaTZmkc5lOtxaJvvUIR358LFC27zxwD0lKQtALhfsgZpUe6Bby1hOp87vWqMT8xKumNe8ZnU9cPONquBf9BFdAFzdkYDqmh19n1vLWI5spa2bqzMrp+aZHhkarIfbtlJ0pbL5PQCUJz8HN9ayap3cXu7EwjeVvrXgIzSXx43ly5erU6eivc8999yjKVOmuAYlkC5mg7CwsKtmA4+PupacnHzFnUwjIyOVlpamCxcuqEKFClcsk5CQoLFjx3q6tEIVZ/SlEmRGAMVSvD7na7PJbrMpwNdHlYIqqHbVyrq5Tk3VCed+HIAn2IrVV23ytdvkY7crOMBfUWEhujEmSq1r15CPnbGUAE8oXl/VxXPY3NC+fXt99dVX+vrrr7V+/XodOnRIKSkpOn/+vEJDQxUbG6u2bduqT58+at++ffFq0jU6vPSoUaM0cuRI1+tLR3TKUmzVitp6JLXI5yj62KR72t+kf71zwbOFAXBxOC199tpKnc3IKfIyvj52Pf5IHz01/XkPVgYgt+NnMzV73Hdu/YSy+fjq3YQxur/lRI/VBSCv9fvPaOmUDVdveJmVX8zWTbWKvrPQz89PXbp0UZcuXdx+L3d4POhUr15dKSl5L9JPSUlRaGhovkdzJCkgIEABAQGeLq1Qf2hbW8/M/qnI7R2W9Id2sfL79WZmADzPT9JDt9TSpJX75Cji3iSnJfWJp68CZSk6MFB33BCh5TuOF/mauiB/H917U4z8uCk3UGbaNoxUnapB2n/yfJF2TNgk1a1WUa1iKnm4suLx+LHf+Ph4LV26NM+0JUuWKD4+3tNvXSL3tohWWAXfIo0Q42OTalUJ0m3cNBQoc3+4pXaR2/rYbbq9YTXVrlrRgxUByM+g9rFFDjl2m/RQm1oevxEhgLxsNpsG3Va3yEdfLUmDbq17zd5w2+2gc+7cOW3evFmbN2+WdHH46M2bN7vGtB81alSe4d8ee+wx7d27V88884y2b9+uf/3rX/r888/15z//uXQ+gYdU8PfRxD/cLLvNVmjY8bFJAX4+er9va9mLOm4mgFITUyVICb9tLuninqWC+NhtiggJ0Ou/v7FsCgOQR7v64XqiY72rtrPbpGY1wvRUl0ZlUBWAyz3cppZ+06x6od+p0sXv3HuaR+nBuLK9vMQdbged9evXq1WrVmrVqpUkaeTIkWrVqpVGjx4tSTp69Kgr9EhSbGysFi5cqCVLlqhFixYaN26cJk2adE0PLX1JfL2qmj64raoGXzyNLneO8fn1L1ejcpBmP95ON0Qx+gvgLQ/cHKO3H2ypoICLNw605emrF180rxGmuU+0V0QIp6wB3vJ010Z69jeN5edju+JHlM+vE+66IVLTB7dVBX9uBAp4g91u0zt9WukPbWvLZtMVO/ztv07rF19bbz/Y8pre0V+i4aXLSlGHkPOUHIdTX287pk+/T9JPP++WZCmuRSM9HB+r2+qHX9MbGChPLmQ59N8fj2j2hgPasStJdsuhjnFN1bd9XbWMqXTNHloHypsz57M0a8MhLfjxsPbtOyA/Z47uvr2l+rWvp/oRwd4uD8Cvjpy5oH+v3qvPl6xTlj1AMbVq6s4m1dWnTS1VD/PejsOiZgOCjhuyMzL0z/6/kyT96eNZXMwMXKPoq8D1gb4KXPuuxX5a1GzAQPQAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMU6ygM2HCBNWpU0eBgYG65ZZbtG7dugLbJiYmymaz5XkEBgYWu2AAAAAAuBq3g85nn32mkSNHasyYMdq4caNatGihrl276tixYwUuExoaqqNHj7oe+/fvL1HRAAAAAFAYt4POm2++qcGDB2vgwIFq0qSJJk6cqKCgIH300UcFLmOz2VS9enXXIzIyskRFAwAAAEBh3Ao6WVlZ2rBhg+66667/rcBu11133aU1a9YUuNy5c+dUu3ZtxcTEqGfPnvr5558LfZ/MzEylpaXleQAAAABAUbkVdE6cOCGHw3HFEZnIyEglJyfnu0yjRo300Ucfaf78+Zo2bZqcTqfatWunQ4cOFfg+CQkJCgsLcz1iYmLcKRMAAABAOefxUdfi4+PVr18/tWzZUh06dNCcOXNUrVo1vf/++wUuM2rUKKWmproeBw8e9HSZAAAAAAzi607j8PBw+fj4KCUlJc/0lJQUVa9evUjr8PPzU6tWrbR79+4C2wQEBCggIMCd0gAAAADAxa0jOv7+/mrdurWWLl3qmuZ0OrV06VLFx8cXaR0Oh0NbtmxRVFSUe5UCAAAAQBG5dURHkkaOHKn+/fvr5ptvVps2bfTWW28pPT1dAwcOlCT169dPNWrUUEJCgiTpxRdfVNu2bVW/fn2dOXNGr7/+uvbv369HH320dD8JAAAAAPzK7aDTu3dvHT9+XKNHj1ZycrJatmypRYsWuQYoOHDggOz2/x0oOn36tAYPHqzk5GRVrlxZrVu31urVq9WkSZPS+xQAAAAAkIvbQUeShg0bpmHDhuU7b/ny5Xlejx8/XuPHjy/O2wAAAABAsXh81DUAAAAAKGsEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BJ1StH//fj311FNq3LixKlasqCpVqiguLk6vv/66zp8/7+3yAOTj2Weflc1mcz2WL1/u7ZIASMrKytKkSZPUtWtXRUVFKSAgQMHBwWrUqJEGDhyo1atXe7tEoNw6duyYFixYoNGjR6tbt24KDw93fY8OGDCgSOs4f/685syZo8cff1xxcXGqXLmy/Pz8VLVqVcXHx+uFF15QcnJyier0LdHScPnvf/+rP/zhD0pLS3NNO3/+vNavX6/169dr0qRJWrhwoerXr+/FKgHktnnzZr355pveLgPAZfbv36977rlHP//8c57pWVlZ2rlzp3bu3KnExEQ9+eSTevvtt2Wz2bxUKVA+RUZGlmj5n376Se3bt9e5c+eumHfq1CmtXbtWa9eu1fjx4/XBBx+od+/exXofjuiUgk2bNql3795KS0tTcHCwXn75Za1evVpLly7V4MGDJUk7d+7UPffco7Nnz3q5WgCS5HQ6NWTIEOXk5CgiIsLb5QD4VXZ2dp6Qc+ONNyoxMVFr1qzR4sWLNXr0aFWsWFGS9M477+jVV1/1ZrlAuVerVi116dLFrWXS0tJcIad9+/ZKSEjQkiVLtHHjRn311Vf64x//KLvdrrS0ND388MP68ssvi1UbR3RKwfDhw3XhwgX5+vpq8eLFio+Pd82744471KBBAz3zzDPauXOnxo0bpxdeeMF7xQKQJP3zn//UDz/8oMaNG+u+++5TQkKCt0sCIGn+/PmukBMfH69vv/1WPj4+rvmdO3dWjx49FB8fr+zsbL366qv6y1/+Il9fftIAZWX06NGKi4tTXFycIiMjlZSUpNjY2CIvb7fb9cADD2jMmDFq0qTJFfO7dOmibt266b777pPD4dCTTz6pXbt2uX30liM6JbRu3Tp9++23kqRBgwblCTmXPPXUU7rhhhskSW+//bays7PLtEYAeR04cEDPP/+8JGnixIny9/f3ckUALsl97c2oUaPyhJxLWrdure7du0uSzpw5o23btpVZfQCksWPHqnv37sU+ha1du3b67LPP8g05l/Ts2VP333+/JGnPnj3atGmT2+9D0CmhefPmuZ4PHDgw3zZ2u139+vWTdPEf8jfffFMWpQEowNChQ3Xu3Dn1799fHTp08HY5AHLJyspyPa9bt26B7erVq5fvMgDM0alTJ9fzPXv2uL08QaeEvvvuO0lSxYoV1bp16wLb5f4xtWrVKo/XBSB/n3/+uRYsWKAqVarojTfe8HY5AC7TqFEj1/O9e/cW2O7Sjx6bzaYGDRp4vC4AZS8zM9P1PL+ju1dD0CmhS4fL69evX+j5wY0bN75iGQBl68yZMxo+fLgk6dVXX1V4eLiXKwJwuT59+ig0NFTSxX7qcDiuaLNp0yYtXLhQkvTQQw+52gMwy4oVK1zPL10G4g6CTglkZGToxIkTkqSaNWsW2rZy5cquUWIOHjzo8doAXOmZZ55RcnKy2rdvr0GDBnm7HAD5CA8P1yeffKKgoCCtWrVKcXFxmjp1qtauXauvv/5aY8eOVYcOHZSVlaWbbrpJ48aN83bJADzgxx9/dO3QaN68ebGCDkOUlEDuoaKDg4Ov2r5ixYpKT0/Pd8xwAJ717bffatKkSfL19dXEiRO57wZwDevRo4c2bNigcePGafLkyerfv3+e+ZGRkXrppZc0ePBgBQUFealKAJ6SmZmpRx991HVE9+WXXy7WejiiUwIZGRmu50UZtSkgIECSdOHCBY/VBOBKWVlZGjJkiCzL0p///Gc1a9bM2yUBKERWVpamTp2q+fPny7KsK+anpKRo2rRp+vrrr71QHQBPGzZsmNavXy9J6t+/v+69995irYegUwKBgYGu50UZ8eXSBVUVKlTwWE0ArvSPf/xD27dvV61atTRmzBhvlwOgEOnp6brrrruUkJCgU6dO6ZlnntG2bduUmZmp1NRULV68WLfeeqvWr1+vXr166c033/R2yQBKUUJCgiZNmiRJiouL04QJE4q9LoJOCYSEhLieF+V0tPT0dElFO80NQOnYvmOH62ag77zzjutaOQDXphdeeMF1f7rJkyfr1VdfVePGjeXv76/Q0FB17txZ33zzjTp16iTLsvT000/rxx9/9HLVAErD+++/r7/+9a+SLg7k9cUXX5Toe5trdIroyJkLmrVuv5ZV7SCbLFVcvV+/a1NHVatW1cmTJ3Xo0KFClz99+rQr6MTExJRFyUC5tCvlrOZuOKDl4R3lYzn09Qf/VbZ8VLduXZ0/f14zZsy4YpmtW7e6ni9btkzJycmSpHvvvZdgBHiAZVnaeOC0vvzxsFaHd5SvM0d1tyTr7pYx+uijjyRJDRs2vOLanEt8fX310ksv6dZbb5XT6VRiYqLGjx9flh8BKBccTkvLd57Qt1XaKcvur7OLdqrjDdV1e4NqsttL91rXTz/9VE888YQkqXbt2lqyZEmJR0cl6FxFcmqGnp+/VV9vS5FNkkIuju//y5Ldev3rPYro8bROT39Ju3fvVk5OToFDTG/fvt31vDijRgAo3PbkND0/b6t+SDotH7tNzuCGsklyyq6awz7RqU1fqM/DfSVnTqHreemll1zP9+3bR9ABStmq3Sc09r8/a2fKOfnabXIEN5JNlv48+2eNWbhdOfU7SOvmqlWrVoWuJ/e963J/xwIoHXM3HdJri3boaGqG7KHNZMmmnesOKXHtQUVXCtSzv2msni1rlMp7/ec//1G/fv3kdDoVFRWlpUuXXnVE46Lg1LVCHDx1Xj3e/U7Lth+TZUlOS3LafC4+rIsp90JEU1XvO04XnD7asGFDgevKPQ54+/bty6J8oNzYeOC07puwWhv3n5F0sW9av/ZV2Wyy+wUq9OZeinhgrOTD/h3AWxb+dFR9J3+vXSkXT/fOcVqybPaLfVVSWoZDlTsNUpWuw5SdU/hOiezsbNfzwu5jB8B97y3foz9/9qOOpl4ceMtp85FlsyvHeXFwkCNnMjR8xmZ9sHJPid9r6dKleuCBB5STk6OqVatqyZIlqlevXonXKxUz6EyYMEF16tRRYGCgbrnlFq1bt67Q9jNnzlTjxo0VGBio5s2b64svvihWsWXJ4bQ0MPEHnUzPksN55Ygvl1iyybdylMJ7PKspU6bk28bpdGrq1KmSpEqVKqlTp04eqRkoj1LPZ2vglHXKzHHIkc/oTJfY7HYF1Wmhv83aJMuy8jxyD1DwzTffuKbXqVOnDD4BUD7sSjmr4TM2ybKkgnvqRSEtf6N1pwKVU0jYyb0DMTY2tpSqBPDN9mN6dVHRjpL+44vtWrHzeLHfa/Xq1erZs6cyMzMVFhamr776Sk2bNi32+i7ndtD57LPPNHLkSI0ZM0YbN25UixYt1LVrVx07dizf9qtXr1afPn00aNAgbdq0Sb169VKvXr3ynBN/LVq587h2HztXaMi5xGb3UYU6LTT1P8u0Zs2aK+aPGzdO27ZtkyQNHz5cfn5+pV4vUF7N3HBQaRdyVISuKsuSZvxwQKfTrz5KIoDSNWV1kixdPeRIF6/hsW7orJf+nv+9M06fPq1nn33W9bp79+6lUyQA/Wv5bvkU8fIbH5s0cXnxjups3rxZ99xzj9LT01WxYkUtXLgwzymppcHtY71vvvmmBg8erIEDB0qSJk6cqIULF+qjjz7Sc889d0X7t99+W7/5zW/09NNPS7p4/vuSJUv07rvvauLEiSUs33M+XpMkH7utSEFHkiynQxVu7KIuXbro2WeeUcfbb9eFjAx9PnOmJk2eLElq0KCB/jR0qLJz3X8HQPFZlqXEVfvcWibHaemz7/dpULvarmmOXHuNc7Ky6KNAKTubkaPZGw4V+TvVZrPJNyRcr0+boA3rf1Dfhx9WbGysMjMy9P26dXrn3Xd14OBBSdIdnTqp0+2302+BUrD7eLp+SDpd5PYOS1qz96Q+/Ox/Z3ft3r1biYmJedoNGDAgz+s9e/aoa9euOnPmjCTp73//u8LCwgo9EBIREaGIiIgi1yZJNiu/O3EVICsrS0FBQZo1a5Z69erlmt6/f3+dOXNG8+fPv2KZWrVqaeTIkRoxYoRr2pgxYzRv3rwCh4PMzMx03XNGktLS0hQTE6PU1FSFhoYWtdwSaTl2sc5cyL56w1yyT+zXkclD851XLaSiBt0ap/AQLmwGSkumzV8f1Bnk1jI2y6n66Xv0m+P/u9HgV1t3askvuyRJj3Vsq/oRVUu1TqC8S/aP0Mwav3VrGcuRo9S1nyv1u+kFtqkfUVX92rVWkD9nSgClYXtwQy2pdqfbyx3/z+s6v21FgfMvjxuJiYmugyZFNWbMGL3wwguSLmaDsLCwq2YDt47onDhxQg6HQ5GRkXmmR0ZGFjjiSXJycr7tLw3fmp+EhASNHTvWndJKXZbD6fYylUPCVL9hrLYdPaYz5zPka7epanBFtYiJUvv6deTv6+OBSoHyy2Fzv09ZsinHxoXLQFly2N3vqz52m+pHRelMtSpKTjunjOxs2W02hQQGKKZKJbWqFa2m0ZGy2Up3iFugPMspxveqJNl8/Uu5ktJxTX7bjxo1SiNHjnS9vnREpyxVDvLX+awLRW5vk9Swfm199uMvnisKQB6ZOU5NefmbIl2fc4mvj11xt7XXn3o86pr2Jw/UBuB/dh9P15wJa91byO6jQX8crEen/t0zRQG4wpJtx/XNZz+5vdy8GZ/oriaRV2/4qwEDBlxxOpsnuBV0wsPD5ePjo5SUlDzTU1JSVL169XyXqV69ulvtJSkgIEABAQHulFbqerWK1nvL97j1A+q+m2rKLzDQc0UByMNPUucmkfp627Ein/uf47TUo1UMfRUoQ41rBqhetYraezy9SIMRSBcHD+lOXwXKVMemUQry/1nnsxxFXiY4wFe3NijZjT09xa1R1/z9/dW6dWstXbrUNc3pdGrp0qWKj4/Pd5n4+Pg87SVpyZIlBba/VvRpU0tFv3pJ8ve1q1er0rlpEoCi6xdfx40LnKVaVYLUrh7X4ABlyWazaUD72CKHHB+7TR0bVVNMlSCP1gUgryB/X/WOi5GPvWinhPrYbHowLkaBftfm5RluDy89cuRIffjhh/r444+1bds2Pf7440pPT3ddUNSvXz+NGjXK1X748OFatGiRxo0bp+3bt+uFF17Q+vXrNWzYsNL7FB5Qs3KQht1Rv8jt/9a9iUICuRgSKGvt6lXV3c2r62r/k22/Pv7eq5nsRfwHDqD0/L51TTWvEXbVH1B2mxTga9eou28oo8oA5PZEx/qqFhxw1b7qY7cpIjRAj3UsnZt7eoLbQad379564403NHr0aLVs2VKbN2/WokWLXAMOHDhwQEePHnW1b9eunaZPn64PPvhALVq00KxZszRv3jw1a9as9D6Fh4zs3FBDbq8rSflu7EvT/np3Y/VtW/uK+QA8z2az6c0HWqrzr+cGF9RXfX1sevehm3R7w2plXSIASYF+Pvr4kTZqGn1xhKT8fkPZbRdPg/lk0C1qGBlSxhUCkKRqIQH6dEhbRYYGuHYS5nZpWlRYoGYMaavwYO9eblIYt4aX9paiDiHnKT8kndLHq5P05ZZk153X/X3s6tUqWv3i66hZjbAyrwlAXk6npWXbj+njNUn6dtcJ1/TgAF/1aROjP7StrdpVGd4d8LbMHIcW/nRUU1YlacvhVNf0asEB6hdfWw+2qaVqIdfuDyegvDibka3ZGw4pcXWSkk6ed02vUzVIA9vH6retayo4wDvjmhU1GxB03KkjI1vH0jJlt0mRoYGq6KWNC6Bwp9KzdPJcpvx87IqqFKgAhnYHrknH0jJ05kK2Kvj5KCosUL4+bp9oAsDDLMvSkdQMncvIUXCgr6LDAr0+rLtH7qNT3oUG+imU63CAa16Viv6qUvHaHNMfwP9EhAYqIpRR1YBrmc1mU41KFbxdRrGw6wQAAACAcQg6AAAAAIxD0AEAAABgHIIOAAAAAOMQdAAAAAAYh6ADAAAAwDgEHQAAAADGIegAAAAAMA5BBwAAAIBxCDoAAAAAjEPQAQAAAGAcgg4AAAAA4xB0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxiHoAAAAADAOQQcAAACAcXy9XUBRWJYlSUpLS/NyJQAAAAC86VImuJQRCnJdBJ2zZ89KkmJiYrxcCQAAAIBrwdmzZxUWFlbgfJt1tSh0DXA6nTpy5IhCQkJks9m8WktaWppiYmJ08OBBhYaGerUWlA62qZnYruZhm5qJ7Woetql5rrVtalmWzp49q+joaNntBV+Jc10c0bHb7apZs6a3y8gjNDT0mtjQKD1sUzOxXc3DNjUT29U8bFPzXEvbtLAjOZcwGAEAAAAA4xB0AAAAABiHoOOmgIAAjRkzRgEBAd4uBaWEbWomtqt52KZmYruah21qnut1m14XgxEAAAAAgDs4ogMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9DJx4QJE1SnTh0FBgbqlltu0bp16wptP3PmTDVu3FiBgYFq3ry5vvjiizKqFEXlzjZNTEyUzWbL8wgMDCzDanE1K1eu1L333qvo6GjZbDbNmzfvqsssX75cN910kwICAlS/fn0lJiZ6vE64x93tunz58iv6qs1mU3JyctkUjKtKSEhQXFycQkJCFBERoV69emnHjh1XXY7v1WtXcbYp36vXvvfee0833nijQkNDFRoaqvj4eH355ZeFLnM99FOCzmU+++wzjRw5UmPGjNHGjRvVokULde3aVceOHcu3/erVq9WnTx8NGjRImzZtUq9evdSrVy9t3bq1jCtHQdzdppIUGhqqo0ePuh779+8vw4pxNenp6WrRooUmTJhQpPb79u3TPffco06dOmnz5s0aMWKEHn30UX311VcerhTucHe7XrJjx448/TUiIsJDFcJdK1as0NChQ7V27VotWbJE2dnZ6tKli9LT0wtchu/Va1txtqnE9+q1rmbNmnrllVe0YcMGrV+/XnfccYd69uypn3/+Od/2100/tZBHmzZtrKFDh7peOxwOKzo62kpISMi3/QMPPGDdc889eabdcsst1h//+EeP1omic3ebTpkyxQoLCyuj6lBSkqy5c+cW2uaZZ56xmjZtmmda7969ra5du3qwMpREUbbrN998Y0myTp8+XSY1oeSOHTtmSbJWrFhRYBu+V68vRdmmfK9enypXrmxNmjQp33nXSz/liE4uWVlZ2rBhg+666y7XNLvdrrvuuktr1qzJd5k1a9bkaS9JXbt2LbA9ylZxtqkknTt3TrVr11ZMTEyhezRwfaCfmq1ly5aKiopS586dtWrVKm+Xg0KkpqZKkqpUqVJgG/rr9aUo21Tie/V64nA4NGPGDKWnpys+Pj7fNtdLPyXo5HLixAk5HA5FRkbmmR4ZGVngOd/JyclutUfZKs42bdSokT766CPNnz9f06ZNk9PpVLt27XTo0KGyKBkeUFA/TUtL04ULF7xUFUoqKipKEydO1OzZszV79mzFxMSoY8eO2rhxo7dLQz6cTqdGjBih9u3bq1mzZgW243v1+lHUbcr36vVhy5YtCg4OVkBAgB577DHNnTtXTZo0ybft9dJPfb1dAHCtiY+Pz7MHo127drrhhhv0/vvv66WXXvJiZQBya9SokRo1auR63a5dO+3Zs0fjx4/XJ5984sXKkJ+hQ4dq69at+u6777xdCkpJUbcp36vXh0aNGmnz5s1KTU3VrFmz1L9/f61YsaLAsHM94IhOLuHh4fLx8VFKSkqe6SkpKapevXq+y1SvXt2t9ihbxdmml/Pz81OrVq20e/duT5SIMlBQPw0NDVWFChW8VBU8oU2bNvTVa9CwYcO0YMECffPNN6pZs2ahbflevT64s00vx/fqtcnf31/169dX69atlZCQoBYtWujtt9/Ot+310k8JOrn4+/urdevWWrp0qWua0+nU0qVLCzxHMT4+Pk97SVqyZEmB7VG2irNNL+dwOLRlyxZFRUV5qkx4GP20/Ni8eTN99RpiWZaGDRumuXPnatmyZYqNjb3qMvTXa1txtunl+F69PjidTmVmZuY777rpp94eDeFaM2PGDCsgIMBKTEy0fvnlF2vIkCFWpUqVrOTkZMuyLKtv377Wc88952q/atUqy9fX13rjjTesbdu2WWPGjLH8/PysLVu2eOsj4DLubtOxY8daX331lbVnzx5rw4YN1oMPPmgFBgZaP//8s7c+Ai5z9uxZa9OmTdamTZssSdabb75pbdq0ydq/f79lWZb13HPPWX379nW137t3rxUUFGQ9/fTT1rZt26wJEyZYPj4+1qJFi7z1EZAPd7fr+PHjrXnz5lm7du2ytmzZYg0fPtyy2+3W119/7a2PgMs8/vjjVlhYmLV8+XLr6NGjrsf58+ddbfhevb4UZ5vyvXrte+6556wVK1ZY+/bts3766Sfrueees2w2m7V48WLLsq7ffkrQycc777xj1apVy/L397fatGljrV271jWvQ4cOVv/+/fO0//zzz62GDRta/v7+VtOmTa2FCxeWccW4Gne26YgRI1xtIyMjrbvvvtvauHGjF6pGQS4NK3z549J27N+/v9WhQ4crlmnZsqXl7+9v1a1b15oyZUqZ143CubtdX331VatevXpWYGCgVaVKFatjx47WsmXLvFM88pXf9pSUp//xvXp9Kc425Xv12vfII49YtWvXtvz9/a1q1apZd955pyvkWNb1209tlmVZZXf8CAAAAAA8j2t0AAAAABiHoAMAAADAOAQdAAAAAMYh6AAAAAAwDkEHAAAAgHEIOgAAAACMQ9ABAAAAYByCDgAAAADjEHQAAAAAGIegAwAAAMA4BB0AAAAAxvl/xIMMcifZfLQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import quantax as qtx\n", "import matplotlib.pyplot as plt\n", "\n", "# 4x4 square lattice, PBC by default, 8 spin-up, 8 spin-down\n", "lattice = qtx.sites.Square(4, Nparticles=(8, 8))\n", "\n", "print(lattice.shape) # output: (1, 4, 4), 1 spin per unit cell, and Lx = Ly = 4\n", "lattice.plot()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "3b10f68a", "metadata": {}, "source": [ "Now we define the Heisenberg Hamiltonian\n", "\n", "$$\n", "H = \\sum_{\\left< i,j \\right>} (\\sigma_i^z \\sigma_j^z + \\sigma_i^x \\sigma_j^x + \\sigma_i^y \\sigma_j^y)\n", "= \\sum_{\\left< i,j \\right>} (\\sigma_i^z \\sigma_j^z + 2\\sigma_i^+ \\sigma_j^- + 2\\sigma_i^- \\sigma_j^+)\n", "$$\n", "\n", "We use the last equation to avoid complex numbers in the Hamiltonian. It's implemented an instance of {py:class}`~quantax.operator.Operator` in Quantax." ] }, { "cell_type": "code", "execution_count": 2, "id": "702fa795", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+4.0 Sᶻ₀ Sᶻ₄ +4.0 Sᶻ₀ Sᶻ₁ +4.0 Sᶻ₁ Sᶻ₅ +4.0 Sᶻ₁ Sᶻ₂ +4.0 Sᶻ₂ Sᶻ₆ +4.0 Sᶻ₂ Sᶻ₃ +4.0 Sᶻ₃ Sᶻ₇ +4.0 Sᶻ₃ Sᶻ₀ +4.0 Sᶻ₄ Sᶻ₈ +4.0 Sᶻ₄ Sᶻ₅ +4.0 Sᶻ₅ Sᶻ₉ +4.0 Sᶻ₅ Sᶻ₆ +4.0 Sᶻ₆ Sᶻ₁₀ +4.0 Sᶻ₆ Sᶻ₇ +4.0 Sᶻ₇ Sᶻ₁₁ +4.0 Sᶻ₇ Sᶻ₄ +4.0 Sᶻ₈ Sᶻ₁₂ +4.0 Sᶻ₈ Sᶻ₉ +4.0 Sᶻ₉ Sᶻ₁₃ +4.0 Sᶻ₉ Sᶻ₁₀ +4.0 Sᶻ₁₀ Sᶻ₁₄ +4.0 Sᶻ₁₀ Sᶻ₁₁ +4.0 Sᶻ₁₁ Sᶻ₁₅ +4.0 Sᶻ₁₁ Sᶻ₈ +4.0 Sᶻ₁₂ Sᶻ₀ +4.0 Sᶻ₁₂ Sᶻ₁₃ +4.0 Sᶻ₁₃ Sᶻ₁ +4.0 Sᶻ₁₃ Sᶻ₁₄ +4.0 Sᶻ₁₄ Sᶻ₂ +4.0 Sᶻ₁₄ Sᶻ₁₅ +4.0 Sᶻ₁₅ Sᶻ₃ +4.0 Sᶻ₁₅ Sᶻ₁₂ +2.0 S⁻₀ S⁺₄ +2.0 S⁻₀ S⁺₁ +2.0 S⁻₁ S⁺₅ +2.0 S⁻₁ S⁺₂ +2.0 S⁻₂ S⁺₆ +2.0 S⁻₂ S⁺₃ +2.0 S⁻₃ S⁺₇ +2.0 S⁻₃ S⁺₀ +2.0 S⁻₄ S⁺₈ +2.0 S⁻₄ S⁺₅ +2.0 S⁻₅ S⁺₉ +2.0 S⁻₅ S⁺₆ +2.0 S⁻₆ S⁺₁₀ +2.0 S⁻₆ S⁺₇ +2.0 S⁻₇ S⁺₁₁ +2.0 S⁻₇ S⁺₄ +2.0 S⁻₈ S⁺₁₂ +2.0 S⁻₈ S⁺₉ +2.0 S⁻₉ S⁺₁₃ +2.0 S⁻₉ S⁺₁₀ +2.0 S⁻₁₀ S⁺₁₄ +2.0 S⁻₁₀ S⁺₁₁ +2.0 S⁻₁₁ S⁺₁₅ +2.0 S⁻₁₁ S⁺₈ +2.0 S⁻₁₂ S⁺₀ +2.0 S⁻₁₂ S⁺₁₃ +2.0 S⁻₁₃ S⁺₁ +2.0 S⁻₁₃ S⁺₁₄ +2.0 S⁻₁₄ S⁺₂ +2.0 S⁻₁₄ S⁺₁₅ +2.0 S⁻₁₅ S⁺₃ +2.0 S⁻₁₅ S⁺₁₂ +2.0 S⁺₀ S⁻₄ +2.0 S⁺₀ S⁻₁ +2.0 S⁺₁ S⁻₅ +2.0 S⁺₁ S⁻₂ +2.0 S⁺₂ S⁻₆ +2.0 S⁺₂ S⁻₃ +2.0 S⁺₃ S⁻₇ +2.0 S⁺₃ S⁻₀ +2.0 S⁺₄ S⁻₈ +2.0 S⁺₄ S⁻₅ +2.0 S⁺₅ S⁻₉ +2.0 S⁺₅ S⁻₆ +2.0 S⁺₆ S⁻₁₀ +2.0 S⁺₆ S⁻₇ +2.0 S⁺₇ S⁻₁₁ +2.0 S⁺₇ S⁻₄ +2.0 S⁺₈ S⁻₁₂ +2.0 S⁺₈ S⁻₉ +2.0 S⁺₉ S⁻₁₃ +2.0 S⁺₉ S⁻₁₀ +2.0 S⁺₁₀ S⁻₁₄ +2.0 S⁺₁₀ S⁻₁₁ +2.0 S⁺₁₁ S⁻₁₅ +2.0 S⁺₁₁ S⁻₈ +2.0 S⁺₁₂ S⁻₀ +2.0 S⁺₁₂ S⁻₁₃ +2.0 S⁺₁₃ S⁻₁ +2.0 S⁺₁₃ S⁻₁₄ +2.0 S⁺₁₄ S⁻₂ +2.0 S⁺₁₄ S⁻₁₅ +2.0 S⁺₁₅ S⁻₃ +2.0 S⁺₁₅ S⁻₁₂\n" ] } ], "source": [ "# Define Heisenberg Hamiltonian\n", "\n", "# Method 1: Using built-in operator\n", "H = qtx.operator.Heisenberg()\n", "\n", "\n", "# Method2: Customized operator\n", "from quantax.operator import sigma_z, sigma_p, sigma_m\n", "\n", "H = 0\n", "\n", "Lx, Ly = lattice.shape[1:]\n", "for x in range(Lx):\n", " for y in range(Ly):\n", " # Periodic boundary taken into account automatically\n", " H += sigma_z(x, y) @ sigma_z(x + 1, y) + sigma_z(x, y) @ sigma_z(x, y + 1)\n", " H += 2 * (sigma_m(x, y) @ sigma_p(x + 1, y) + sigma_p(x, y) @ sigma_m(x + 1, y))\n", " H += 2 * (sigma_m(x, y) @ sigma_p(x, y + 1) + sigma_p(x, y) @ sigma_m(x, y + 1))\n", "\n", "print(H)" ] }, { "cell_type": "markdown", "id": "69aff0f6", "metadata": {}, "source": [ "## Direct ED\n", "\n", "Then we use ED to obtain 2 lowest eigenstates of H, which is internally performed by QuSpin.\n", "In Quantax, we implement {py:class}`~quantax.state.DenseState` to process dense wavefunctions." ] }, { "cell_type": "code", "execution_count": 3, "id": "fede224a", "metadata": {}, "outputs": [], "source": [ "E, wf = H.diagonalize(k=2)\n", "\n", "dense0 = qtx.state.DenseState(wf[:, 0]) # Ground state\n", "dense1 = qtx.state.DenseState(wf[:, 1]) # First excited state" ] }, { "cell_type": "markdown", "id": "62e0fdc5", "metadata": {}, "source": [ "One can easily compute the wavefunction amplitude $\\psi(s) = \\left$ from\n", "{py:class}`~quantax.state.DenseState`." ] }, { "cell_type": "code", "execution_count": 4, "id": "870d61c3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "s = [ 1 -1 -1 1 1 -1 -1 1 1 -1 -1 1 1 -1 1 -1]\n", " = [0.00286136]\n" ] } ], "source": [ "s = qtx.utils.rand_states()\n", "print(\"s =\", s) # +1/-1 represents spin-up/down\n", "print(\" =\", dense0(s))" ] }, { "cell_type": "markdown", "id": "822627cd", "metadata": {}, "source": [ "There are convenient expressions to compute $\\left<\\psi|\\phi\\right>$ and $\\left<\\psi|H|\\psi\\right>$." ] }, { "cell_type": "code", "execution_count": 5, "id": "fc6a3563", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " = 1.0000000000000018\n", " = -2.220446049250313e-16\n", " = -44.91393283371552 \t E0 = -44.913932833715506\n", " = -42.59953949065392 \t E1 = -42.599539490653896\n" ] } ], "source": [ "print(\" =\", dense0 @ dense0) # output: 1, normalization\n", "print(\" =\", dense0 @ dense1) # output: 0, orthogonality\n", "\n", "print(\" =\", dense0 @ H @ dense0, \"\\t\", \"E0 =\", E[0])\n", "print(\" =\", dense1 @ H @ dense1, \"\\t\", \"E1 =\", E[1])" ] }, { "cell_type": "markdown", "id": "c36bbdd9", "metadata": {}, "source": [ "One can also measure customized observables from the ED results. Here we show the spin-spin correlation as an example. The spin-spin correlation is defined as\n", "\n", "$$\n", "C_{ij} = \\left< \\mathbf{S}_i \\cdot \\mathbf{S}_j \\right>\n", "= \\frac{1}{4} \\left< \\sigma_i^z \\sigma_j^z + 2\\sigma_i^+ \\sigma_j^- + 2\\sigma_i^- \\sigma_j^+ \\right>\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "id": "abc992c1", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " = -0.35089010026340245\n", " = 0.21376528539942638\n" ] } ], "source": [ "def correlation(i, j):\n", " # sigma_z(x, y) and sigma_z(i) with i=4x+y are equivalent in 4x4 lattice\n", " return (\n", " sigma_z(i) @ sigma_z(j)\n", " + 2 * sigma_p(i) @ sigma_m(j)\n", " + 2 * sigma_m(i) @ sigma_p(j)\n", " ) / 4\n", "\n", "C1 = correlation(0, 1)\n", "C2 = correlation(0, 2)\n", "\n", "print(\" =\", dense0 @ C1 @ dense0)\n", "print(\" =\", dense0 @ C2 @ dense0)" ] }, { "cell_type": "markdown", "id": "05c4e46e", "metadata": {}, "source": [ "## ED with symmetries\n", "\n", "A common trick in ED is reducing the Hilbert space dimension using {py:class}`~quantax.symmetry.Symmetry`.\n", "In the square lattice Heisenberg model, the system is invariant under translation,\n", "rotation, mirror flip, and spin flip (time reversal symmetry)." ] }, { "cell_type": "code", "execution_count": 8, "id": "ea8bdcfe", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/aochen/quantax_env/lib/python3.12/site-packages/quantax/symmetry/symmetry.py:288: GeneralBasisWarning: using non-commuting symmetries can lead to unwanted behaviour of general basis, make sure that quantum numbers are invariant under non-commuting symmetries!\n", " basis = spin_basis_general(\n" ] } ], "source": [ "from quantax.symmetry import TransND, C4v, SpinInverse\n", "\n", "# The symmetry superposition can be performed by \"@\"\n", "symm = TransND() @ C4v() @ SpinInverse()\n", "\n", "# We will get \"GeneralBasisWarning\" from QuSpin, but it doesn't matter here\n", "E0_symm, wf0_symm = H.diagonalize(symm)\n", "dense0_symm = qtx.state.DenseState(wf0_symm, symm)" ] }, { "cell_type": "markdown", "id": "d1b5a9ff", "metadata": {}, "source": [ "The wavefunction is identical to the unsymmetrized one up to a possible extra minus sign,\n", "and the measurement gives the same result." ] }, { "cell_type": "code", "execution_count": 9, "id": "ba8e64bb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "E0 = [-44.91393283]\n", "Check overlap: -1.0000000000000009\n", " = [-0.00286136]\n", " = -44.91393283371545\n" ] } ], "source": [ "print(\"E0 =\", E0_symm)\n", "print(\"Check overlap:\", dense0_symm @ dense0)\n", "print(\" =\", dense0_symm(s))\n", "print(\" =\", dense0_symm @ H @ dense0_symm)" ] }, { "cell_type": "markdown", "id": "7f089366", "metadata": {}, "source": [ "However, be careful that some operators like the spin-spin correlation $C_{ij}$ are ill-defined in the symmetrized Hilbert space. One cannot directly measure them in this case." ] }, { "cell_type": "code", "execution_count": 10, "id": "1fc4d946", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wrong : 0.15074025930850307\n", "Wrong : 0.23714288042475265\n" ] } ], "source": [ "print(\"Wrong :\", dense0_symm @ C1 @ dense0_symm)\n", "print(\"Wrong :\", dense0_symm @ C2 @ dense0_symm)" ] }, { "cell_type": "markdown", "id": "8931027b", "metadata": {}, "source": [ "One can obtain the excited state as the ground state in a specific symmetry sector. Here we compute the triplet excitation energy with -1 eigenvalue under spin inversion." ] }, { "cell_type": "code", "execution_count": 11, "id": "f53096a7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Excited state energy = [-42.59953949]\n", "Check overlap: -1.0000000000000013\n" ] } ], "source": [ "spin_inv = SpinInverse(-1)\n", "\n", "E1_symm, wf1_symm = H.diagonalize(spin_inv)\n", "print(\"Excited state energy =\", E1_symm) # It's the first excited state shown above\n", "\n", "dense1_symm = qtx.state.DenseState(wf1_symm, spin_inv)\n", "print(\"Check overlap:\", dense1_symm @ dense1)" ] } ], "metadata": { "kernelspec": { "display_name": "quantax_env", "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.12.11" } }, "nbformat": 4, "nbformat_minor": 5 }