diff --git a/Notebooks/testing_against_ED-Copy1.ipynb b/Notebooks/testing_against_ED-Copy1.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..6d8233b33e800ca48b465b70c0db8308254cd601 --- /dev/null +++ b/Notebooks/testing_against_ED-Copy1.ipynb @@ -0,0 +1,480 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing against ED - Hubbard on a ring" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we use the [pyALF](https://git.physik.uni-wuerzburg.de/ALF/pyALF) interface to run ALF with the Mz choice of Hubbard Stratonovitch transformation on a four site ring, at $U/t=4$ and inverse temperature $\\beta t = 2$. For this set of parameters, the exact internal energy reads: \n", + "$$\n", + " \\left\\langle -t \\sum_{\\langle i,j\\rangle, \\sigma} \\hat{c}_{i,\\sigma}^{\\dagger} \\hat{c}_{j,\\sigma}^{\\phantom\\dagger} + U \\sum_{i=1}^{N} \\hat{n}_{i,\\uparrow}\\hat{n}_{j,\\downarrow} \\right\\rangle = -1.47261997 t \n", + "$$\n", + " \n", + "To reproduce this result we will have to carry out a systematic $\\Delta \\tau t$ extrapolation keeping $\\Delta \\tau t L_\\text{Trotter} = 2$ constant. \n", + "Recall that the formulation of the auxiliary field QMC approach is based on the Trotter decomposition \n", + "$$\n", + "e^{ -\\Delta \\tau \\left( \\hat{A} + \\hat{B} \\right) } = e^{ -\\Delta \\tau \\hat{A}/2} e^{ -\\Delta \\tau \\hat{B} } e^{ -\\Delta \\tau \\hat{A}/2} + \\mathcal{O} \\left (\\Delta \\tau^3\\right)\n", + "$$\n", + "The overall error produced by this approximation is of the order $\\Delta \\tau^2$. \n", + "\n", + "Bellow we go through the steps for performing this extrapolation: setting the simulation parameters, running it and analyzing the data.\n", + "\n", + "---\n", + "\n", + "**1.** Import Simulation class from the py_alf python module, which provides the interface with ALF, as well as mathematics and plotting packages:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "from py_alf import Simulation # Interface with ALF\n", + "\n", + "import numpy as np # Numerical library\n", + "from scipy.optimize import curve_fit # Numerical library\n", + "import matplotlib.pyplot as plt # Plotting library" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "V_{U}(i) = U \\hat{n}_{i,\\uparrow}\\hat{n}_{j,\\downarrow}\n", + "$$\n", + "$$\n", + "V_{U,{\\rm SUN}}(i) = \\frac{U}{2}\\left[\n", + " \\sum_{\\sigma=\\uparrow,\\downarrow} \\left( \\hat{n}_{i,\\sigma}-\\frac{1}{2} \\right)\n", + " \\right]^2 \n", + " = \\frac{U}{2} \\left( 2\\hat{n}_{i,\\uparrow}\\hat{n}_{j,\\downarrow}\n", + " - \\hat{n}_{i,\\uparrow}\n", + " - \\hat{n}_{j,\\downarrow}\n", + " + 1 \\right)\n", + "$$\n", + "With $\\langle \\hat{n}_{i,\\uparrow} \\rangle = \\langle \\hat{n}_{i,\\downarrow} \\rangle = \\frac{1}{2}$\n", + "$$\n", + "\\left\\langle V_{U,{\\rm SUN}}(i) \\right\\rangle = \\left\\langle V_{U}(i) \\right\\rangle\n", + "$$\n", + "\n", + "\n", + "$$\n", + "V_{U,{\\rm MZ}}(i) = \\frac{U}{2} \\left( \\hat{n}_{i,\\uparrow} - \\hat{n}_{i,\\downarrow} \\right)^2\n", + " = \\frac{U}{2} \\left( 2\\hat{n}_{i,\\uparrow}\\hat{n}_{i,\\downarrow} \n", + " - \\hat{n}_{i,\\uparrow} - \\hat{n}_{i,\\downarrow} \\right)\n", + "$$\n", + "With $\\langle \\hat{n}_{i,\\uparrow} \\rangle = \\langle \\hat{n}_{i,\\downarrow} \\rangle = \\frac{1}{2}$\n", + "$$\n", + "\\left\\langle V_{U,{\\rm MZ}}(i) \\right\\rangle = \\left\\langle V_{U}(i) \\right\\rangle - \\frac{U}{2}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def extrapolate_dtau(ham_name, sim_dict, dtaus, do_ED):\n", + " sims = []\n", + " for i, dtau in enumerate(dtaus):\n", + " print(dtau)\n", + " sim_dict['Dtau'] = dtau\n", + " if i == 0 and do_ED:\n", + " sim_dict['do_ED'] = True\n", + " else:\n", + " sim_dict['do_ED'] = False\n", + " sim = Simulation(ham_name, sim_dict, \n", + " alf_dir=os.getenv('ALF_DIR', './ALF'),\n", + " branch='136-automatic-ed-form-op_t-and-op_v'\n", + " )\n", + " if i == 0:\n", + " sim.compile(target=ham_name) \n", + " sim.run()\n", + " sims.append(sim)\n", + " \n", + " if do_ED:\n", + " with open('{}/ED_Energy'.format(sims[0].sim_dir), 'r') as f:\n", + " E = float(f.read())\n", + " else:\n", + " E = None\n", + " \n", + " ener = np.empty((len(sims), 2))\n", + " for i, sim in enumerate(sims):\n", + " print(sim.sim_dir) # Directory containing the simulation output\n", + " sim.analysis() # Perform default analysis\n", + " ener[i] = sim.get_obs(['Ener_scalJ'])['Ener_scalJ']['obs'] # Store internal energy\n", + " \n", + " return ener, E" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.05\n", + "Checking out branch 136-automatic-ed-form-op_t-and-op_v\n", + "Compiling ALF... Done.\n", + "Prepare directory \"/home/jonas/Programs/pyALF/Notebooks/ALF_data/Hubbard_N_leg_ladder_L1=4_L2=1_Checkerboard=False_Symm=True_T=1.0_U=4.0_Tperp=0.0_beta=2.0_Mz=True_Dtau=0.05\" for Monte Carlo run.\n", + "Resuming previous run.\n", + "Run /home/jonas/Programs/ALF/Prog/Hubbard.out\n" + ] + } + ], + "source": [ + "dtaus = np.array([0.05, 0.1, 0.2])\n", + "ener, E = extrapolate_dtau(\n", + " 'Hubbard', # Hamiltonian\n", + " { # Model and simulation parameters for each Simulation instance\n", + " 'Model': 'Hubbard', # Base model\n", + " 'Lattice_type': 'N_leg_ladder', # Lattice type\n", + " 'L1': 4, # Lattice length in the first unit vector direction\n", + " 'L2': 1, # Lattice length in the second unit vector direction\n", + " 'Checkerboard': False, # Whether checkerboard decomposition is used or not\n", + " 'Symm': True, # Whether symmetrization takes place\n", + " 'ham_T': 1.0, # Hopping parameter\n", + " 'ham_U': 4.0, # Hubbard interaction\n", + " 'ham_Tperp': 0.0, # For bilayer systems\n", + " 'beta': 2.0, # Inverse temperature\n", + " 'Ltau': 0, # '1' for time-displaced Green functions; '0' otherwise \n", + " 'NSweep': 1000, # Number of sweeps per bin\n", + " 'NBin': 100, # Number of bins\n", + " 'Mz': True, # If true, sets the M_z-Hubbard model: Nf=2, N_sum=1,\n", + " 'do_ED': True,\n", + " }, # HS field couples to z-component of magnetizatio\n", + " dtaus,\n", + " do_ED=True\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dB3wUxRfHXyBAUAihBCkJhEgJHSUiiBSREpKQ0EGxYKEo/lWUKghYUZodBELvJQUIoaqIBcSg9FBCAAlFkBqQEkj+70324nFckjvucvX35fM+s7e7t/eYu8xv38zOG4/MzEwCAAAAjFHA2E4AAAAAIgEAACBXEEkAAACASAAAADAfRBIAAAByxDPHI05ImTJlMgMCAuztBgAAOBXbt2//h5909XV5kRCBSExMtLcbAADgVHh4eBzL6Ri6mwAAAOQIRAIAAABEAgAAgPm41JgEAOBu0tPTKTU1la5fv47qcXO8vLzIz8+PChUqZPJ7IBIAuDgiEMWLF1cPdvAApb3dAXZC8vSdO3dO/R6qVKli8vswJgGAiyMRROnSpSEQbo4H3yDI78DciBIiAYCbNBAAeNzD7wAiAQAAACJhS3pM3aIMAJBFwYIFqUGDBtn2ySefqP0tW7akGjVqUL169SgoKIhee+01unjxok2rrVixYqo8evQoLVy4EF+ZAYgkAAD/8fPnREc231kj8lr2W0DRokVpx44d2TZs2LDsYwsWLKBdu3YpK1KkCEVGRuZ5vdu3b1v9W4NIGAciAQD4j4oPEy3r/Z9QSCmvZX8+U7hwYRo3bhz99ddftHPnTqN3/KNGjaJHH32UtmzZIvmGqEWLFtSwYUNq164dnTp1Sp335ZdfUq1atVR00rNnT7VvzJgxNGHChOxr1alTR4mCPiJcP/30k4p0Pvvss3z8nzoXnhYOgnTjYgxbTbZG/IhVYi7nFuRCjp/g88K1fUu4qKGd4sN2kY810I4N5+IlNrlleJ33r7PEVwCACVRpTtRtdpYwBPOfX+KMrNey3wKuXbumGl8dw4cPpx49ehjtlqpfvz7t379flfpcvXpVNe7vv/++mvshArFixQry9fWlJUuW0IgRI2jmzJmqK+vIkSMqKjGn60reJ0ISHx9/7/9RF8TSeRJ72DqzTTXh3DfYkti8dTu44c/+lbAoTOTikrZdiwu5BajNVoFtI++rzudbP8YEANyJCIIIxOZxRM2HWCwQ+t1Npj7PbwwRkC5duqjtAwcO0J49e6hNmzbZ3U/ly5dX2xJB9OrVizp27KgM2LG7ib/MJLYDeZ3HDbwfF2FsUTkcl+eyurMt0nZJp+RivvYNtiO8nczWyBJfAQAmIl1MEkGIQEhpOEaRj0hjv3v3bqpZUzon7p4tLEKhE5LatWtnj3HIe9avX6+OrV69mgYMGKC6o6Qr6tatW+Tp6UkZGRnZ18Lsc8cbk5BRL/7F0X/f0p00Y/ubv/hD2uuKbMf1jqdq++6C9aUvW6LY2bNnreUvAO6JbgxCuphajfiv68kGQiFdSNIN5e/vr6KB3JAnouTvXcYmdO/du3evEoLjx4/TE088ocY3pLvpypUrarb5H3/8oc6VUrqjDJFZ6Wlpadb/j7m6SHDjK109e4xY3o8gZL1fxh/OsABsz+W0p/SiCPU2I+cYjUH5utPYgsWkbxIAYAEn/rhzDEI3RiH7rTAmoTP9p5uka0hEQcYbZNxBxhlMGeRevnw5DR06VI1dyDV//fVXFYk888wzVLduXXrooYdo4MCB5OPjo7qpzp8/r86bMmUKVa9e/a5rig8Sccj1MHBtxpgEN76t8zonD5qyRbBYhHLpxebN2/P5us/IQd721MY1GhpEDv56r6W76qSFfgAA8uLxN+/eJ0Jh4bhETo+sbtq0yeRrSESgjzT4mzffHeH8/PPPRsdEdN1ROV1Xkt599913JvvjLuR7dxOLwXA2P7YAbTD6e51AaIgI7ed9Igw6Vsq5LCBF2CQTVTW2bfntKwAAACuKBDfgndikcW/Ctpq31WOqXFZgSzDxMj0NuppEWPZysZRtH9tatgF4sgkAAJzsEVhuuGO5iDWyX7qGQo3sl9hyk8G+3jlc+yMuxAAAANgJzLjWeGHtC8oAAABAJAAAAJgAIgkAwF0gkzGASAAAbIYkVXj22WezX8ssaJnXFB6u0riZjUySmzx5srXcM4mjR4+quRz2Rh4b1tWbbMv8kPwEkQQAIN+5//77Va4lmVQnbNiwgSpWNJpEwSFF4haLmiXkR2pzASIBALALadfT6cTFa7T92AWrXbN9+/Yqr5KwaNEieuopSbSQhcyGlmR8Muu5cePGam0JXYrvF198US1OFBgYqNKACzJj+/Dhw2pC3eDBg9W+8ePH0yOPPKKuMXr0aKM+rF27lh5++GE1q/rJJ59U+7Zt20aPPfaYmqEtpSQPFGbPnk3dunWjDh06UNu2be+4juR+euGFF7Jndv/www9GG3BJD/L000+r84T58+dTo0aNlN/9+vVT4iHWu3dvFaXIebrZ3vJ/TkzMSqz9zz//qNQihpHNt99+q86X60ma82XLlqnryP+veXPLEzNaIwssAMDFEGHYfzqNMjKJekVtpQUvN6aGlUtafF1Z20HSfEtXiYiANP7SsAnSqEtjGxcXR99//z0999xz2VljJW24NMKSV0lyNr3yyisqrbdEJrpzZDb1oUOHVIMvyf8iIiLUbGz9hlJyPfXp00ftr1KlihImQVbEk32SkmPjxo30zjvvUHR0tDq2ZcsW5WupUqXuWH/im2++UaUkFhT/REQOHjyokhDqI/6In/J5SUlJKqX5L7/8omZ3v/rqq2rBJUlUeOLECXWeYGp6cxGN/v37q3U2Bg0apPaJyKxbt05FadZa4Q8ikU93YZev31J/bNb44wLAlmxNOacEQki/laFeW+N3LHf40tBKFBEaGnpXKg1dw9yqVSs6d+4cXbqkVg6gsLAwtTaEWNmyZenvv/++69oiEmIiNLpUGyIa+iKxdetW9VoabEEafkE+5/nnn1fny9iJJAvU0aZNm+zzDP393//+ly0ylStXViJhmJhQogbd50nKD8lMK9GOIF1v8v+RSCUlJUVdT/6vhlGLOTRt2lRFJd27d6fOnSXbkeVAJJzkLgwAW9E4sDQV8CD1Gy7kWUC9thZyhy93vdIVI0KQ2xoSWSsIkBIHHZIq3Nj4gLxfMshKF05OyDm6a+rz7rvvqm6h2NhYJWLSzaM/lpLTtUxB//3yHhGjsWPH3nWerMQnEYBEKEuXLlWLJ+mnNzc1tbl0P/3222+qW0+6oCTSKl3asu8PA9c2uAsDwJmQm5qgcsXJr2RRq9/kSBeTLEGq66PXIXf40vUiiICUKVOGvL2z1yfLM623LF8qDasuWZ9035w5c+aO9zRp0oR+/PHH7DThuu4miSR0g+gyDmEKzfX8lQhCllyVrrDckDEQyVyr80s+/9ixY2q8QcRAMtV+8MEH2SnNpTtJIg9B3mdKPcg4jSzvKt16UoeSNt1SEEk40V0YALaiuFchZdaOgv38/OiNN2SRyjuRAWoZCJbumvvuu4/mzJmT63Xk7li6VmSQVgbEZdBa+vxFCATpp5dBYunO0SGP3E6bNk11w0ijLMfkKashQ4aoO/xJkyapri5TeJXHE2Q8QMRO7vhFXPQjHmPIutsffvih6k6Sz5dxCYkcJEOt/N91UYMu0pCIS7qN5s2bl6Nf0lXVtWtXlV79q6++UoPY0m0mUYuIkuESsPeCh6lhkzMQHBycqXsawFx0KTlmhcyy2I/QLzarMYkvej6EriZgd6TxNLbSW16T6YQl/bIaXeDavwfuhtsua/IYOx+RxM+fE1V8+M5akVW4ZJEVY7n17XgXBoCtgDgAHRiTEIFY1puCLp29c/lGQ+EAAAA3BCKhLc/46qHfqOPxff+t72vhSlwAAOAKQCQEFoRlDwRQ5In9PLDxEgQCAAAgEv+xa+cc+qzwdepeuQod2DErq8sJAAAAIgkRhID171HdgiXoYEGirqWL0sD1/ejA7oX4eQC3BYtwAR3obuKnmLy78jPOZapR3TJ16ZX6r9DW+4tR1z/G0sAfBtKB81nJvgAA947MlJYZwDqT3Eu58fHHH9ukumV+xoQJE3I9Jy4ujvbt4/FKM5G5GtbAWte5VyAS8pirNkjtWcCTXm3wKq3ttpH61+9PW09tpa6rutJbm96igxcO2vWLAsCZkQljkiJCZ5LF9V5EQuZ16Sad2Yq4exQJVwEiYYQSRUrQgAYDaG2XtdS3Xl/69eSv1GVlF3p709t06MIhW39HALgkkg5DUlnoUnNL6vDp06crAZHkdxJx9OrVS+VTkslfMstZ0nxLqgnJBBscHKwyqOqnBZdUFkOHDlWJ9cSSk5PVfkl/ITOQZUa3lJJGwxD5bEm+J7OUJUXGv//+qxb0WblypUpHLv5I2guxkJAQatiwITVr1kxlgRUk3YfM+JZrSD4oY4hv+utgSCQzceJElU5E/JL/n8zilhnUuS02JLz22mvZaUQkfUeLFi2UT5Ki5NSpU+Z+HTmCyXR5iMX/HvofPVfrOZq7by4tSFpAG45toLYBbal/vf5UtWRVq30RANiCT7d9SvvPZzVquaE7R5eJIDeCSgXR0EZDcz1H1+jrkGR8PXr0oK+//lplLZVUHRcuXFCpvAXZr0sDLiIhQjJr1qzsBvajjz5S2VllLQZpXCWdty4Dq+R8khTdc+fOpTfffJPi4+NVgyrpxyX9huR4ev3111WEoI+k69B9/siRI2nGjBkqM6skJZTGWdJfCPJ5kkivWrVqKpmeiJekN5f/g4iXfI4ulbixdOnik7xHkGR+ssaFpBiXBIPiu+RykjU15HONJSQ0RLLWip8iLJJ6RNKRjxgxQv0/rQFEwgyxeLbms9lisf7oemoX0E51Sz3o86BVvgwAXL27yRBJxS0L5QwYMEBlQs0JScUtDacOaVwlD5NkhJW7ZukO0omEbjEjKQcOHJi9LkRMTIzalmVUJV+TIbKeg4iDrMMgd/ZyR26I7JfoQhYj0nHjxg1VyjoRunTn8hkSNRgiqcwlwd/JkyfV+hYlS5akSpUqqYZe1rGQdS0KFCigEhRKSvRy5crddQ1DREDFd6lLQYSzfPnyeb7PVCASZuDj5UOvP/y6iizm7JtDC5MW0rqj6ygkIESJRaBPoNW+GADyg7zu+PMjl1luyPiC5BISEZGsqJIAMK+U29KtI4PNv//+u2pkJRLRT6Wtf/ed0524sf1yHYkupLtJunGke8eYvz4+PkYFTzDlzl8iEsnqevr0aRVZCJJRVkRDuo0k8Z90mxmmB9dPHS7ojss4jXS7iRDmBxiTuEexeOPhN9SYxYt1XqRNqZuo44qONGTzEEq5lGLt7wgAl0Wylsp4gyxEJGnEdQv+SEOpv/iPPpcvX1aiUaJECXW3vWbNmjuOS3eLrtRlhZVlSRcvXpzdID/++ON3XVdSbssduHyuLg24YTpu6Q6SRYQk+tE10LoISLLS6n9GTogwyHkiFLouLBmfkay08v+WVfhkDMVYNCURk0Qucr4sYiTIuI4IjE4kxP+9e/fm+PnmApGwgJJeJenNhm/Sui7r6IU6L9Cm45uo04pOdKLQDLrhcdpa3xEATo9uTEJnMjgt6zBERUWpgVsZAJY1GiSVttC3b1/VfSQD14bInb5028jdswiLNM76SCMqayp88cUX2etFy9rYMqYh15TU23LMEFnLQd4n3Tay2px+oz5+/Hj1mTJoLQIg4xXih/igG2SWa8pYhAxc61bVM4a8R0RH1rDQdQvJ/1MyWMtgvFxf//N1+Pv7q9ThunrRrcJXuHBhJTjSvSU+Sf1Kl5i1QKpwK4bX56+fp9l7Z9Ps3Qsok9Ip/MEw9XRUlRJZyxcC4Cypwm3V3WRtpJtGGltZcAcYB6nC7Ugpr1L0VsO36JfEOnTOcx1999d3lHAkgUKrhFK/ev0ooESAPd0DwGScTRxA/oHupnzAk4rTA7e60prOa9Qg98ZjGylyRSS989M7dOzy3X2NAADrII/LIoqwLhCJfKR00dL0dvDbtKbLGvX4rMyxiIyLpBE/j6C/Lt89mQeA/MKVVqAEtv0d4BFYG6zqVaZoGRr0yCDqXac3zdwzk5YeWEqrU1ZTeGC46oby9/bPDzcAUMhErXPnzql1oU15RBO4rkCc49+B/B5sJhL8g5MZJWPYZFSsETuRmMu5BbmQ4yf4vHBtnzyrVkM7xYftIh9rwPul8z6JTZddbyvv72+Jr46AiMWQR4aox2Zn7J5Byw4uo/iUeOrwYAc1wO1fHGIBrI/MPUhNTVWPSQL3xosFIqe5KPkVSexh68w21YRz39Aafm/dDm74e+i2WRgmcqH/3NhhEQwL/XNYsZBJTSIWushi1eFVFPFgBPWp1wdiAayKPHsvz/YDYPMxCW7Ek9jyzKXNAiDSFcYWlcNxiYG7sy2yxB9nw/c+XyUWMmbRM6in6oKKiI2g0b+OphNXTtjbPQAAsM48CW7jZf76oJy6m/j4ci7GshXXzgs3OC65uifx/mDttXQ3yZRByc99mW0kH/spLz+Cg4Mz5RlpZ+XMv2dUN9Tyg8spIzODIqtGqsiiYrGK9nYNAODCcJu7Xdf+mh1J8Js3su0xYpEmfrgIwhl2YHsupz1lEEVInttK/B6ZUvgW20K+jncO1+/Llijm7H2uZe8rS8MfHU4JnROoa/WutPLwSgqPCaf3trxHJ6+ctLd7AAA3JN8jCT4mEcSzbLfYZFhdGvsYPvcZ7biMi0jfSkPel2ru9V0pkjDk9NXTKrKIPhRNmfyvY9WO1KduH6pQrIK9XQMAuBAWRRKWwh88nM2PTbqQJOXh9zqB0GjNtl9fINhhX+1pKNmW1KrV2Nwuc165+8vRiMYjVGTRpVoXWpG8gsJiw+j9Le/TqSvWW1QEAABywiKR4Aa8E5s07jIxYDVvr9P2V2BLMPEyPY0MWMsYxS6+hqRXlPGM/iwi5y3x1dnFYmTjkdliEZscS6GxofTh1g9VtAEAAPkFEvw5IRJFRO2OopjkGPLgf52rdaaX676sxAQAAKzZ3QSRcGJkMFvEQiILEQuJMkQsHrj/AXu7BgBwIiASbiAW03dPp7hDcSrtgjwZ9VKdlyAWAACTgEi4CTIBb/qu6WqAu4BHgSyxqPuSerQWAAByAiLhZqSmparIQsTCs4AndaveTaUAkRneAABgCETCTTmedlxFFjIpD2IBAMgJiISbI2Ixbdc0lURQxKJ7je4qspBEgwAA4IGnm4AgCx2JWEh68kIFCimxeKHOCxALANwcD4gEMBSLqbumKrEoXKAw9ajRQ4mFrKQHAHA/PCASwBiy3rYusihSsIgSi961e0MsAHAzPCASIDeOXjqqxGL1kdVKLHrW6EnP134eYgGAmwCRACZx5NIRJRYJRxKyxCKop4osSnmVQg0C4MJAJIDZYiFjFgkpCeTl6UVPBT2lxKKkV0nUJAAuCEQC3BMpl1Lo253f0toja5VYPB30tOqGglgA4FpAJIBFpFzUxOLoWirqWZServk0bf69FnlSMVrST7LEAwBcVSTyfdEh4PwE+gTSuBbjKCYihpr7NVer5SUXGU5nPGPp0o1L9nYPAJCPQCSAyVQtWZXGtxhP0RHRVCyjLv3juYbaRbejr/78CmIBgIsCkQBmU61kNfJL70uBN0ZT0wpN1RNRIdEh9PWfX0MsAHAxIBLgnvHKrEgTW05UkUWTCk3UE1EiFt/s+AZiAYCLAJEAFlO9ZHWa1HISLe+wnBqXb6wGudtHt6fJOybT5ZuXUcMAODEQCWA1apSqQZ898ZkSi0blG9GUnVMoZHkITdkxhdJupqGmAXBCIBIgX8Ti8yc+p2UdltEj5R6hyTsnqwFuEQ2IBQDOBUQC5BtBpYLoi1Zf0NLwpRT8QLDqfhKxkO6oKzevoOYBcAIgEiDfqVm6Jn3Z6ktaEr6EGj7QUA1si1hM3TkVYgGAgwORADajVula9FWrr2hx+GJ6qOxD9PWOr5VYyCO0iCwAcEwgEsDm1C5dm75+8mtaHJYlFjIZLyQmRK3HfTX9Kr4RABwIiASwG7XLZInForBFVN+3Pn3555cqsojaHQWxAMBBgEgAu1OnTB365slvlFjUK1OPvvjjCzUpT8Ti3/R/7e0eAG4NRAI4lFhMbj2ZFoYuVNsiFhJZSEJBiAUA9gEiARyOur51aUrrKbQgdIHqkvr8j89VZDFrzyyIBQA2BiIBHJZ6vvXo29bf0rz289STUZO2T6L2Me1p9p7ZEAsAbAREAtwTadfT6cTFa7T92IV8r8EGZRvQt22yxKJGyRo0cftEJRZz9s6ha7eu5fvnA+DOQCSA2Ygw7D+dRqkXrlGvqK02EQqdWExrO43mtp+rxGJC4gTVDQWxAMBBRcLDw6Mb2162DLbgPM4tyPYnW7zevgZsW9l2sCWyNdI7Npwtme0AWztL/ATWZWvKOcrIzNpOv5WhXtsSmVshYjEnZI5a20LEQrLOzt07l67fum5TXwBwdSyNJPawdWbbbMK5b7AlGewbx/ZeZmZmAy5Haa9FIGpx0ZOtNlsI22QRGQt9BVaicWBpKuCRtV3Is4B6bQ8efuBhimobRbNDZlNVn6o0PnG86oaav28+xAIARxAJbtyT2A7kdR438H5chLFFGV6CzVvbLsF2UtuOZFvM177BdoS3k9myowxgXxpWLklB5YqTX8mitODlxuq1Xf15oCFFtYuiWe1mUWCJQPr0908pNCaUFiQtgFgA4CRjEp+zDWHLMNj/Jtt4FpHjXE5gG67tr8gm+3Skavvugt/bV+uqSjx79qx1vQY5UtyrEFX0KWp3gdAnuFwwzWg3g2a2m0mVvSvTJ9s+yRaLG7dv2Ns9AFxTJDw8PDay7TFicrefJ3xeOBdnOCLYbuTwK2wD+Zi/lGwzdG8zcq7WC26wMzNzGluwmK+vrykuARdH1rCYFTJLiUUl70pZYhEdSguTFkIsALC2SHDj25qtjhFbYeJnNGWLYLE4yuVitla8PV879jxbjLa9TK9LSSIHEQ4dfnpdUQCYLhbcBTWj7QzyK+5HY7eNVZHFov2LIBYAOEp3E4vJcDY/tgBtMPp73n5GOywNfwttuxXbIW17pZzLYlKErQpvV2Pblt++AteDfz9qKVUZ3JZBbr9ifvTxbx8rsVi8fzHdvH3T3i4C4NKPwHZik7v+JmyreXudtr8CW4IJl+jDNpHP3cnlx2x9ZSeLyF4ulrLtY1vLNoD33bbEV+DeiFg8Wv5RJRbT2kyjCvdXoI9++0iJxZL9SyAWAOT0t8ONbw6HnI/g4ODMxMREe7vhFvSYukWVS/rJ/YHzIb/7Lae2qCVVd57dSeXuL0d96vahjlU7UuGChe3tHgC2vonaLuO6xo5hxjVw28jisQqPqVQfU1tPpbL3laUPtn5AYbFhtPTAUkq/nW5vFwFwCCASwK1RYlHxMZrffn6WWBT9TyyWHVwGsQBuD0QCAH2xCJ2v0pSXKVqG3t/yPoXHhlP0wWhKz0BkAdwTiAQABmLxeMXH1VoWk5+cTKW8StGYLWOoQ2wHijkUA7EAbgdEAoAcxKKZXzNaGLZQLa3qU8SHRv86WolF7KFYiAVwGyASAOQhFs39mqv1t0UsShQpQaN+HUURsRFKLG5l3EL9AZcGIgGAGWKxOGwxfd3qa/Iu4p0lFnERFJccB7EALgtEAgAzxaKFfwslFl+1+oqKFSpG7/7yrhKLFckrIBbA5YBIAHCPYtHSvyUtCV9CXz7xpRKLkb+MpMi4SFp5eCXEArgMmHENgJVmcP9w/AeasnMK7T+/X6Uq71evH7Wv0p48C3iijoHTzriGSABgZbH4/vj3NGXHFDpw4QAFeAdQ33p9KbRKKBUsgMUVgWOCtBwA2O6PjZ6s9CQt7bCUPmv5mcoD9c7P71DHFR1pdcpqup2BPJXAucCYBAD58YflUYBaV25Nyzoso0ktJ1GhgoVo2E/DqNPKTpSQkgCxAE4DRAKAfBaLNpXb0PIOy2lii4lU0KMgDf1pKHVe2ZnWHFkDsQAOD0QCABuJRduAthQdEU0TWkxQr4dsHkJdVnahtUfWUkam4fLvADgGEAkAbPkHx+LQLqCdEovxLcarfYM3D6bOKzrT2qMQC+B4QCQAsJNYhASEUExkDI1vPp4y+d/gHweryGLd0XWILIDDAJEAwN5iUYXFIiKGxjUfR7czb9OgHwcpsVh/dD3EAtgdiAQADoDMoZCJd7ERsfRps0+VWLz949vUdVVX2nBsA8QC2A2IBAAOJhahgaFKLD5p9olaGe+tTW9Rt1XdaOOxjRALYHMgEgA4qFiEBYZRXGQcjW02lm7evkkDNw2k7qu603fHvoNYAJuBtBwAOAGyboXMq5i6ayodu3yMgkoFUf/6/amVfys1yxsAS0DuJgBcSCwSjiTQ1J1T6a+0v5RYvFL/FXrC/wmIBbhnIBIAuKBYSC4oiSyOpx2nmqVqKrGQ9OWILIC5QCQAcGGxiE+JV5FF6pVUJRavNniVWvi1gFgAk4FIAODipGekU/zheBVZnLhygmqVrkWv1n9VLbmKyALkBUQCADcVi9qla6vIolnFZhALkCMQCQDcUCxWHV5F03ZNU2JRp3QdeqXBKxALYBSIBABuikzGkzW3RSxOXj1J9crUU2LRtEJTRBYgG4gEAG6OiEXc4Tiavms6nbp6iur51lNjFo9VeAxiAQgiAQDIFovY5Fiavns6nb56mur71ldi0aRCE2koUEtuioeHx/bMzMxgo8dk4XZXITg4ODMxMdHebgDg8Eiaj7jkONUN9fe/f1MD3waqG6pJeYiFO+KRi0gUsPDC3dj2smWwBedxbkG2P9ni9fY1YNvKtoMtka2Rtj+A7Zq2X+xbS/wEANxJ4YKFqXuN7pTQOYFGPjpSdUH129CPnl/7PG05uYVc6eYR2DfB3x62zmybTTj3DbYkg33j2N7jH2QDLkdpr3Uclv2a9bfQTwBADmLRI6iHEosRj45QT0L13dCXeq/tTVtPbYVYAMtEghvvJLYDeZ3HkYAfF2FsUYaXYPPWtkuwncR3ArJurG4AAA4lSURBVIB9xKJnUE8lFu88+g6lpqVSn/V9lFj8duo3iIUbY6tU4Z+zDWEzXO39TbbxLCLHuZzANlzvWBWte+pHtmY5XZiP9dW6qhLPnj1rdccBcCeKFCxCTwU9RQldEmh4o+FKLF5e/zK9sO4F+v307/Z2DziiSHDju5FtjxGLNOUD+LxwLs5wxLHdyOFX2AbyMX8p2WZo+0+xVeL9D3H5FttCvo4u4rgDPmeaDLiI+fr6muISAMAEsXi65tNKLIY1GkZ/Xf6LXlz3ojKIhXthlaebuAHfxMUgvlaikWNjuXiW7RabF5s09jF87jN87BJv+/B2pkfW83eXeNPbnOvrg6ebAMgfrt+6TtGHoilqdxT9c+0falSukco6G1wu1+dVgLs/3WQK/MHD2fzYAvhlT7bvRSC0wzIG0ULbbsV2SHPYV56G0rYDuajGlpLfvgIAjOPl6UW9avaiNZ3X0JBHhtDhi4dVF9TL616mP/7+A9Xmwlj6CGwntlTebMK2mrfXafsrsCWYcIk+bBP53J1cfszWV9vfnG2Xtn85W38WlvOW+AoAsI5YPFvrWVrTZQ0NDh5MyReT1WOzMm7x55k/UcUuCCbTAQDumWu3rtHSA0tp5p6ZdP76eWpcvjENaDCAGpSVp9qBK3Q3QSQAAFYXC5m5LSnKIRbOAUQCAGAT/k3/V4nFrL2zlFhIAkEZ4IZYODYQCQCAzcViyYElNGvPLLpw44JKTS6RhWSfBY4HRAIAYDexWHxgMc3eM1uJxeMVH1dZZ+v61sU34kBAJAAAdheLhfsX0py9c+jijYtqhTzphoJYOAYQCQCAQ3A1/Sot2r+IZu+dTZduXKLmfs2VWNQpU8ferrk1Hni6CQDgaGKxMIkji31zlFi08GuhxKJ2mdr2ds0t8YBIAAAckSs3r2R3Q12+eZla+rWk/g36U+3SEAtbApEAADi8WCxIWqAii7SbadTSv6Ua4K5Zuqa9XXMLIBIAAKdABELEYu6+uWr7Cf8n1KOzQaWC7O2aSwORAAA4FSIQ85Pm07y98ygtPY1a+bdSYlGjVA17u+aSQCQAAE6JjFMs2LeA5u3LEovWlVpT//r9IRZWBiIBAHB6sZi/b74SiyvpV6hN5TbUr14/iIWVgEgAAFwCeVxWuqFEMHRiIZFF9ZLV7e2aUwORAAC4nFhIVCGCIXMu2lZuq8SiWklZnwyYC0QCAOCyYiFPQskTUZL6o20Ai0W9/lS1ZFV7u+ZUQCQAAC4vFjIhT8RC1rZoF9BORRYP+jxob9ecAogEAMAtuHj9YnZkIWIREhBC/er3g1jkAUQCAOB2YiGztyU/lBKLKiGqGyrQJ9DerjkkEAkAgFty4foF1Q0l+aGu37pO7au0V5FFYAmIhT4QCQCAWyNLqYpYSJryG7dvZIlFvX5UpUQVe7vmEEAkAABAEwtZy2Lx/sVKLEKrhCqxCCgR4Nb144FU4QAA8B/nrp3LFoubGTcprEqY6oaq7F3ZLavJAyIBAAB388+1f9T620sOLFFiER4YriKLSt6VnKq6ekzdosol/ZpYXSQK3LtbAADg3JQpWoYGPTKI1nRZQ8/UfIbWH11PEXERNOLnEXT88nF7u+cQQCQAAG6PiMXgRwYrsXi65tO07ug66hDXgd795V23FwuIBAAA6InFkEeG0JrOa+ipoKdozZE1SixG/TKKjqe5Z2QBkQAAAAN87/OloY2GKrHoGdSTVqespojYCBr962hKTUt1q/qCSAAAQC5iMazRMNUN1b1Gd4o/HE8dYjvQmF/H0IkrJ9yi3iASAACQB2XvK0vDHx1OCZ0TqFuNbrTy8EoKjwlXYnHyykmXrj+IBAAAmMgD9z9A7zz6jhKLrtW7KrEIiw2j97a857JiAZEAAAAzKXd/ORrReIQSiy7VutCK5BVKLN7f8j6dunLKperTIpHw8PDoxraXLYMtOI9zC7L9yRavt68+2xa23Wyr2Lz1jg1nS2Y7wNbOEj8BACC/xGJk45HZYhGbHEuhsaH0wZYP6PTV0y5R6ZZGEnvYOrNtNuHcN9iSDPZFsQ3LzMysy2Us22DZyaJQi4uebLXZQtgmi8hY6CsAAOSvWHRKoM5VO1NMcgy1j2lPH2790OnFwiKR4MY9ie1AXudxA+/HRZgmCvrU0BOYDWxdtO1ItsV87RtsR3g7ma2RJb4CAEB+U75YeXq3ybu0utNq6lS1E0UfiqbQmFCnFgtbjUl8zjaELcNIJBKhbXdj89e2K7Lpz1xJ1fYZE6C+bIliZ8+etZ7HAABwj1QoVoFGNRmlxCKyaiRFH8wSi4+2fkR/X/37Hq/qoCLBje9Gtj1GTO7284TPC+fiDEcE240cfpFtgCSX4rI4203d24ycm2ns+nzdaZKYSszX19cUlwAAwGZiMbrJaIrvHE8RD0bQ8oPLlViM/W0snfn3jFN8C555ncCNb2sLP6MpWwQLQSiXXmzevD2fr/sM235+3VZO4n3VtS4pXeSgiyoE6a5yzefLAAAuT8ViFWnMY2Po5bovU9TuKFp6YKkSDJlz8VKdl9SkPbftbmIhGM7mxxagDUZ/LwKhCUNZrRQ/RrJ9q71tpZzL+4uwydJR1di25bevAACQn/gV91NisbLTSgoLDFPrWcgA96fbPqWz/551yUdgO7HJXb8kMV/N2+u0/RXYEky4xFN83kEu92uRwizZySKyl4ulbPvY1rIN4H23LfEVAAAcBf/i/vR+0/dpVadVailVWVZVJxayxoUj4cGNr719sBrBwcGZiYmJ9nYDAADMQtaumLprKsWnxJNnAU+VJ+rFOi+qrLSmgEWHAADAhfH39qcPH/+QVnZcSe0C2tGCpAXUPro9Tfh9gkmRRdr1dDpx8RptP3bB6r4hLQcAADgIlbwr0UePf6TEom1AW5qXNE+JxcTEiWpdbmOIMOw/nUapF65Rr6itVhcKiAQAADgYlb0rK7FYEbmC2lRuQ3P3zVVjFpMSJ90lFltTzlGGNmqQfitDvbYmEAkAAHBQAkoE0MfNPqa4yDh6stKTNGffnCyx2D6Jzl8/r85pHFiaCmgzywp5FlCvrQkGrgEAwElIuZRC03ZNo4SUBPLy9FJLrPau3Zt6Td1Nl6/foi96PkQNK5c0+7oyoVkmJBs9hqebAADA+cRi6s6pag1uEYv7rreg0rfaUHQ/NTfZqiKB7iYAAHAyAksE0qfNP1XdUC39W9K5guvoRKGZ9knLAQAAwDEJ9Amkcc3H0cH9jSmTbuXLZ0AkAADAySmSWT7fro3uJgAAABAJAAAA5oNIAgAAAEQCAACA+SCSAAAAAJEAAABgPogkAAAAQCQAAACYDyIJAAAAOQKRAAAAAJEAAABgPogkAAAAQCQAAACYDyIJAAAAOYJU4QAA4OQs6dck366NSAIAAABEAgAAgPkgkgAAAACRAAAAYD6IJAAAAEAkAAAAmA8iCQAAABAJAAAA5oNIAgAAQI54ZGZm5njQ2fDw8DjLxTELLlGG7R8ruWNN4BfqC78vx6GMC7YTlVkLfF1eJKwgMolcH8H29sMQ+IX6wu/LcfBws3YC3U0AAAAgEgAAAMwHkcSdTDO/Cm0C/EJ94fflOEyztwO29MsDYxIAAAByApEEAACAHIFIAAAAcD+R4MfBQtgOsCWzDTNyXPhSO76L7eG83svbpdg2sB3SypIO4tcYthNsOzQLtbFfM9nOsO0xeI+96ysnv+xWX1z6s/3AlsS2l+0NR6ivPPyyZ315sW1j26n59Z6D1JdXLn7Z9e9RO16Q7U+2eIvrS8YkXM2YgmyH2QLZCrPtZKtlcI58cWuk7tgas/2W13uZcWzDtG354j51EL/GsA2yR31px5qzyY90j8F77FZfefhlt/piyotP2nZxtoMO8vvKzS971pe8LqZtF5L9ctwB6ssjF7/s+veoHX+LbSFbvKV/j64aSTRiS+b/YArbTd5ezBZpcI68nsvHha287cPKWj6P90o5R9uWsqOD+GUplvglP6TNXJw3cl171lduflnKPfvF26fY/tD8S+Miia2ivesrD78sxRK/hCvaOdIYi2U6QH3l5pelWPS759KPizC2KCPvMbu+XFUk5Md9XO91qrbPlHNye+8D/IWckg2tLOsgfgmvaWHnzHsIuy3xKzfsWV95Yff64s8N4OIh7S7UYerLiF92rS+t62QHb55h28B14xD15ZGzX/b+fX3ONoQtw+A991RfrioSEoIZYqjyOZ1jynsdza8pbA+yNWCTL3+iDf3KT/LLL7vXFzccxbiIZnuT/2Avm/n5tvbLrvXFftxmk8+WO+RG7GMdMz/f1n5NsVd9sQ/hXJ5hv7ab+ZluJxKiqv56r+VLPGniObm992+9kK68dgdhd7/4B/G39oOVO4fpWrhqK79yw571lSP2ri+ui0JaQ7yAfYhxlPrKyS9715eeHxe52MQW4ki/r0wDv+xcX03ZIrg+jmrdVK14e75F9cX/EZczxpMtha2K3sBPbYNzwgwGfrbl9V5mvMHAzzgH8au83vsHyo/DVn7pHQ8wMkBst/rKwy+71Zf2ei7b50aua8/fV25+2bO+JDOpj7ZdlO0ntnAHqC/fXPyy+9+jdk5Lg4Hre6ovkx13NtNG/w9qTwmM0Pb1F9O2pXK/0Y7vZgvO7b3a/tJs37Ed0spSDuLXPO3cXWwr9X+kNvJrkRZWp2t3OC85SH3l5Jfd6ot5nC1T+2zpzxYLtXd95eGXPeurHtuf2mfLo8yjHOHvMQ+/7Pr3mItI3FN9IS0HAAAAtxuTAAAAYAUgEgAAACASAAAAzAeRBAAAAIgEAAAA80EkAQAAACIBAADAfP4PAh1xdRifYEUAAAAASUVORK5CYII=\n", + "text/plain": [ + "