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": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "L1 = 4\n", + "L2 = 1\n", + "U = 4.0\n", + "E += L1*L2*U/2\n", + "\n", + "plt.errorbar(dtaus**2, ener[:, 0], ener[:, 1], fmt='.', label=\"Monte carlo results\")\n", + "plt.plot(0, E, 'x', label='ED result')\n", + "\n", + "def func(x, y0, a):\n", + " return y0 + a*x**2\n", + "popt, pcov = curve_fit(func, dtaus, ener[:, 0], sigma=ener[:, 1], absolute_sigma=True)\n", + "perr = np.sqrt(np.diag(pcov))\n", + "xs = np.linspace(0., dtaus.max())\n", + "p = plt.plot(xs**2, func(xs, *popt))\n", + "\n", + "plt.errorbar(0, popt[0], perr[0], label=\"Extrapolated value\", color=p[0].get_color())\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.05\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=False_Dtau=0.05\" for Monte Carlo run.\n", + "Resuming previous run.\n", + "Run /home/jonas/Programs/ALF/Prog/Hubbard.out\n", + "0.1\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=False_Dtau=0.1\" for Monte Carlo run.\n", + "Resuming previous run.\n", + "Run /home/jonas/Programs/ALF/Prog/Hubbard.out\n", + "0.2\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=False_Dtau=0.2\" for Monte Carlo run.\n", + "Resuming previous run.\n", + "Run /home/jonas/Programs/ALF/Prog/Hubbard.out\n", + "/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=False_Dtau=0.05\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/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=False_Dtau=0.05/Ener_scalJ 1\n", + "/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=False_Dtau=0.1\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/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=False_Dtau=0.1/Ener_scalJ 1\n", + "/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=False_Dtau=0.2\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/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=False_Dtau=0.2/Ener_scalJ 1\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': 300, # Number of bins\n", + " 'Mz': False, # 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": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dB3hU1daGVwiB0EIgBClJgECQJlICUqRIN/SugmAFvCiCYISLP/YG6IWrXkVRr16UXqUIhq6CmEgRpAeQ3psKCEn+9W3OhCFOGjOTmSTfy/M9+8w++5zZOTPsNWuXtX2SkpKEEEIIcUQeR5mEEEIIjQQhhJA0oSdBCCGERoIQQkjmoSdBCCEkVfKmeiYbUqJEiaTy5ct7uhqEEJKtiIuLO6UzXYNzvJGAgYiNjfV0NQghJFvh4+NzILVz7G4ihBCSKjQShBBCaCQIIYRknhw1JkEI+TtXr16VQ4cOyeXLl/l4cjn+/v4SEhIifn5+Gb6GRoKQHA4MRJEiRczEDh2g9HR1iIdAnL7Tp0+b70OFChUyfB3HJAjJ4cCDCAoKooHI5fjoDwR8DzLrUdJIEJJLGghCfG7he0AjQQghhEYiPR7+5mEjQojr8fX1lVq1aiXrzTffNPnNmzeX22+/XWrWrClVqlSRJ598Us6dO5elH0HhwoVNun//fvnqq6+y9L2zA/QkCCE3+G6CyL41Nz8RvEa+ExQoUEA2bdqUrJEjRyaf+/LLL2XLli1G+fPnl86dO6d7v4SEBJd/ajQSjqGRIITcoGwdkZkP3TAUSPEa+W4mX758MnbsWPntt99k8+bNDn/xjxkzRu666y5Zt24d4g1Js2bNpG7dutK2bVs5evSoKffvf/9bqlWrZryT++67z+S9+OKLMn78+OR71ahRwxgFe2C41q5dazydf/3rX278S7MXnAJLCLlBhaYiPf973TBEPioS+8n118h3gkuXLpnG18aoUaOkd+/eDrul7rzzTtmxY4dJ7fnjjz9M4/7yyy+btR8wEPPnz5fg4GCZPn26jB49Wj799FPTlbVv3z7jlWSm6wrXwZAsXLjw1v/QHEheJ0fKe2ryoqqqqr7Ow41No6yvJjh/WMt1sPJw7eOqk1axf+q5xZqPUK7bVTut/PWaP8iZuhJCMggMAgzEmrEiTaOdNhD23U0Znc/vCBiQ7t27m+OdO3fK1q1bpXXr1sndT6VLlzbH8CD69OkjXbp0MSKe9SS2qrqpJmWg7NNWwx+QIv9f+qW44QfeYK/m3/jpQQjJGtDFBA8CBgJphSYuMRQZAY39L7/8IlWr4nfn31cLw1DYDEn16tVNt1NKFi1aJGvWrJEFCxbIK6+8Itu2bZO8efNKYmJichmuPs+iMQn9oLarbL/2U0U9gxBN2qsmO/N+hBA3YxuDQBdTi9E3up5SDma7AXQhoRsqNDTUeANpgRlRJ0+eTDYSuBbGAIbg4MGDcs8995jxDXQ3/f7772a1+c8//2zKIkV3VEqwKv3ixYuu/8OyOVk1cI2pEfqzRG6Y8hs8qUZki+pTVTG7/Ar6eqNqtUp/yjhGzw1QxUL40hBCnODwzzePQdjGKJDvgjEJm+xnN6FrCEYB4w0Yd8A4Q0YGuWfNmiXPPfecGbvAPX/44QfjifTt21fuuOMOqV27tgwbNkwCAwNNN9WZM2dMuQ8++EAqV678t3uiDvA4cD8OXN/AJ7X+v+QCPj4xmpRycGq0XjvfKrNKkxGOxiT0HMYfovTcP/S4uVXONiZxmyanVKjEK6rSeu4Rzc+vx4X1+LQe19Xjearq+vpCWnWNjIxMutVNh2xrJD5r99ktXU+It7J9+3aH3Tckd7LdwfdB29k4bV8jb2lMQi9s5WSdGqs6aSWiNPVXBejxFL1vX9Vxu0p+rImZVqD5VzS5Yh3H6bm9egjTz23nCCEkJ3U3aSM/ShWiwowlTFpeAQOBc9r4X5+OcJ2u1kA48oOt2VA4DtckQhXv7roSQghxoZHQBryr6pAeNlQt0uOlVn4Z1eIM3GKslvsFYxJ6fI9qmJWPDlGMU2BFzSzVIDUsZ5ypa1bSe9I6I0IIydVTYLXhnqvJXAf5RzSJcpCPsYtVdq8fTOW+szWBCCGEeBCG5SCEEEIjQQjJOOwyJTboSZBbgo0IyexmNw8+eKN3+dq1aybmUocOZjZ8psEiuf/85z9Z+iHs37/frOXwNKtWrUp+bjjG+hB3QiNBCHE7hQoVMrGWsKgOfPvtt1K2bNlbvl9WG4lratScwR2hzQGNBCHEI1y8fFUOn7skcQfOuuye9957r4mrBKZOnSr3339/8jmshkYwPqx6btCggdlbwhbi+5FHHjGbE4WHh5sw4AArtvfu3WtWUD/77LMmb9y4cVKvXj1zjxdeeMFhHb755hupU6eOWVXdsmVLk7dhwwZp1KiRWaGNFMEDwX//+1/p2bOndOzYUdq0aXPTfRD76eGHH05e2b1y5UqHDTjCgzzwwAOmHJgyZYrUr1/f1HvgwIHGeEAPPfSQ8VJQzrbaG3+zbXHwqVOnTGiRlJ7Nhx9+aMrjfghzPnPmTHMf/H1Nm7om3hZDhRNCbgKGYcexi5KYJNJn8nr58rEGUrecfcScWwN7OyDMN7pKYATQ+KNhA2jU0djOmzdPVqxYIf369UuOGouw4WiEEVcJMZueeOIJE9YbnomtzLJly2T37t2mwUcUiU6dOpkgf/YNJcL2PP744ya/QoUKxjAB7IiHPITkiImJkX/+858ye/b1yZXr1q0zdS1evPhN+0+8//77JkUwQtQPRmTXrl0mCKE9qA/qiffDSmeENP/+++/Fz89P/vGPf5gNlxCo8PDhw6YcyGh4cxiNQYMGmX02RowYYfJgZJYuXWq8NFft8EcjQQi5ifXxp42BAFevJZrXrjAS+IWPhhZeRFTUzTPkv/vuu+SGuUWLFnL69Gk5f/68ed2+fXuzNwRUsmRJOX48OVBDMjASEAwNQFA/GA17I7F+/XrzGg02QMMP8D79+/c35TF2gmCBNlq3bp1cLmV9n3rqqWQjU65cOWMkUgYmhNdge7/ly5ebjZLg7QB0veHvgacSHx9v7oe/NaXXkhkaN25svJJevXpJt24I0O08NBKEkJtoEB4keXzEGAq/vHnMa1eBX/j41YuuGBgCG45iyKHBBjAONhAq3NH4AK5HBFl04aQGytjuac///d//mW6huXPnGiOGbh77sZTU7pUR7K/HNTBGb7zxxt/KYSc+eADwUGbMmGE2T7IPb57R0Obofvrxxx9Ntx66oOBpBQU59/lx4JoQchPwGqqUKiIhxQq4rKvJBrqYsAWprY/eBn7ho+sFwICUKFFCAgJSbj2TelhvbF+KhhUeBED3zYkTJ266pmHDhrJ69erkMOG27iZ4ErZBdIxDZISmdvWFB4EtV9EVlhYYA0HkWlu98P4HDhww4w0wBohUi/0vbCHN0Z0EzwPguow8B4zTYHtXdOvhGSJsurPQkyCE/L3x8fczcqWBACEhIfL009h/7GYwQI2BYHTXFCxYUD7//PM074Nfx+hawSAtBsQxaI0+fxgCgH56DBKjO8cGptx+9NFHphsGjTLOYZZVdHS0+YX/zjvvmK6ujPAPHU/AeACMHX7xw7jYezyOwL7br776qulOwvtjXAKeA3btw99u8xpsngY8LnQb/e9//0u1Xuiq6tGjhwmv/u6775pBbHSbwWuBUUq5BaxbQoVnJ7wlVLgtbtP0gde/sDmR3PA35uTQ0OnBzzfnst3VocIJIbkPGn9ig2MShBBCUoVGghBCCI1EeiQkumfZPCGEZGfoSSj7zu+TTSc3Sfz5eIk9FpvhOdCEEJLToZFQ8vvmlxIFSsj5K+fl4aUPS6d5neTzbZ/L2cuui1tDSHYCs/1sM/5I7sbZ7Ut7qrapElWR6ZT1VW1ULUyR/5Rqp3WfsXb5o1R7rHNtnalnepQpXEbKBZSTmsE15ZXGr0jR/EVlfOx4aTmzpUSvjpYNR6/HgyGE3BpYKY0VwDYh9lJavP7661nyqLE+Y/z48WmWmTdvnvz666+ZvjfWargCV93nVnF2CiwiUiFAyKQMlMUKmu2q5GWU2vhjX+vOqpraCF/R12bli6bVNLlPVV1VRhWjeZW1jFsHDnx9fKVLpS5Gu8/ultm7Z8vXe7+WJfuXGCPSLaKbdK7YWYIKuC5MASG5ASwYswXjywgwEgi0lxL8WIPy5Mm6TpB5aiQQlBCL4XIjTj1p/bC2q67H1U0DbeBDNGmvmpzi1BOqN2EgrPvZ1tHDcExDvgpr6Peo6jtT18wSUSxCRtYfKct7LpfX735dgvyD5F9x/5JWs1rJ8FXD5YcjP0hi0vUVkoSQzINwGAhlYQvNjdDhH3/8sQkDjuB38Dj69Olj4ilh8RdWOSPMN0JNIBJsZGSkiaBqHxYcoSyee+45E1gP2rMHTYeY8BdYgYwV3UgRRiMleG8E38MqZYTI+PPPP82GPgsWLDDhyFEfhL2A2rVrJ3Xr1pUmTZqYKLAA4T6w4hv3QDwoR6Bu9vtgwJN5++23TTgR1At/H1ZxYwV1WpsNgSeffDI5jAjCdzRr1szUCSFKjh49mtmPI1WyajHdBFW0qkiK/MqqJmpEXtMUEaxGqFH4SVMEUllvV+6Qlfc39NoBmkASFhbm4mqL+Of1l44VOxrFn4uXWbtnGe9i2YFlUrZwWelRuYfxPDCmQYi389aGt2THmeuNWlrYymRkXKJK8SryXP3n0ixja/RtIBhf79695b333jNRSxGq4+zZsyaUN0C+zfOAkYAh+eyzz5Ib2Ndee81EZ8VeDGhcEc7bFoEVMZ8QovuLL76QoUOHysKFC02DivDjCL+BGE9DhgwxHoI9CNdhe//nn39ePvnkExOZFUEJ0Tgj/AXA+yGQXkREhAmmB+OF8Ob4G2C88D62UOKOwqWjTrgGIJgf9rhAiHEEGETdEcsJe2rgfR0FJEwJotainjAsCD2CcOSjR482f2eWGAmtZIwmpRycGq0N+vwMXA/Td0LLxulxcwfvj+AwDVSInztDy4TjMge3cjgooPf9SJOPbGE50quPM4QHhkt0vWh5us7TsvzActMdNfHnifL+xvelWWgzYzAalmaYCkIy2t2EUNzYKGfw4MEmEmpqIBQ3Gk4baFwRhwkRYfGrGWMGNiNh28wI6bBhw5L3hZgzZ445xjaqiNeUEuznAOOAfRjwyx6/yFOCfHgX2IzIxpUrpiPE7BNhC3eO94DXkBKEMkeAvyNHjpj9LYoVK2Z+3KKhR/ca9rVAVxoCFCIkeqlSjprem4EBRd3xLAEMZ+nSpdO9zmVGQhvhVk6+R2NVJ238EUAeO3IE6PEUvW9fy0OYo8do3DdgAFzTElZ+qN090F11xMl6uHQ2VFR4lNH+8/uNsZi/Z74s/225lClURq761pPAhEaeriYhfyO9X/zuiGWWFghqh1hCMCKIiooAgOmF3Ea3Dgabf/rpJ9PIwhOxD6Vt/+s7tV/ijvJxH3gX6G5CNw66dxzVNzAwMNXxlYz88odHgqiux44dM54FQERZGA10GyHwH7rNUoYHtw8dDmzn0Xyi2w2G0B24ffRH/4BRqhAV9t7DE1lhGQgAf8+EN8TAtCb5VKdUC1BW8/KrsGNHhGqDu+t6K5QvWl6GRw43Yxfjm42X0IBQOek3X3bnHylPrXhK1hxaw4V6hKQCopZivAEbESGMuG3DHzSU9pv/2HPhwgVjNIoWLWp+bS9ZsuSm8+husaW2qLDYlnTatGnJDfLdd9/9t/si5DZ+geN9bWHAU4bjRncQNhGC92NroG0eEKLS2r9HasAwoBwMha0LC+MziEqLvxu78GEMxZE3BY8JngvKYxMjgHEdGBibkUD9t23blur7Z/UU2K4q/OrHJ7FIj5da+WVUizNwC3SahWtZzJLC0+0Pr0LBXzhDhXln36gGu3tmk7P4+fpJ2/JtZXKbyVLpyqsSlNBOfjn5iwxePljazWknH2z6QI79cczT1STEI9jGJGzC4DT2YZg8ebIZuMUAMPZoQChtMGDAANN9hIHrlOCXPrpt8OsZhgWNsz1oRLGnwsSJE5P3i8be2BjTwD0RehvnUoK9HHAdum2w25x9oz5u3Djznhi0hgHAeAXqgTrYBplxT4xFYODatqueI3ANjA72sLB1C+HvRARrDMbj/vbvbyM0NNSEDrc9F9sufPny5TMGB91bqBOeL7rEXAVDhbsxVPiUxyNl9cHVMmvXLDMbCq7o3WXvlh4RPaRJSBPJmyf7BuFlKOnsHRraW7qbXA26adDYYsMd4hiGCvci/PL4SatyrYwO/35YZu+aLfP2zJMhh4ZIyQIlpUtEF7P2ArOkCPEmsptxIO6DYTmyCBiCIXWGyLIey2TiPROlcvHK8vGWj+Xe2ffKoJhBZrbU1UTHfbCEkIyB6bL0IlxL9u3vyKagi6lFWAujI78fkbl75sqc3XNk6KqhZq0F1lzAuwgtYj+5ixDnwABrRmbekJxN0i2EF6In4QYuXr4qh89dkrgDZ9ONGTW41mBZ2n2pvNviXakRVEM+3fqpRM2JkgHLBsiy/cvkagK9C+IcWKh1+vRpxh/L5SSpgcD3AN+HzEBPwsXAMOw4dlES1WD3mbxevnysQbqbycO7aB7a3AgzoOBdzN09V4avHi7F/YtL50qdzWB3WIDrV5STnA/WHhw6dMhMkyS5G381EKmtRUkNGgkXsz7+tDEQ4Oq1RPM6PSNhT6lCpeSJO5+QAXcMMDOiMDPqi21fyGdbP5O7St0l3St3l5ZhLSWfL5aUEJI+mHuPuf2E3Ao0Ei6mQXiQ5NGuXxgKv7x5zOtbwTePr5kmC53484RZ0Y2V3dFroiUwf6CJRguDUaEo//MTQtwHjYSLgddQpVQRuXD5mky8r3amvIjUKFmwpDxe83F59I5HZf2R9SbI4Jfbv5TPf/1c6t5W18SMal2utQkXQgghroRGwg0U8fczcoWBsCePTx5pVLaR0alLp5K9i1FrR8kbP74hnSp2ku4R3aVSsUoufd/UBudhCDEG4+q/kxDiPXB2UzYF02XhWSzsutCEAmlUppFM2zlNui7oKv2W9JMFexfIpWuX3Do4f+jsJTM4n94sLkJI9oVGIpsD7+Ku0nfJuGbjTJDB4XWHm725R383WlrOaCmv//i67Dq7y+2D84SQnAm7m3IQmC77UI2HpH/1/hJ7PNbMjIKm7phq9u/GNFoEISzoV9ArBucJId4PPYkcCFbW1itVT95q+pbxLp6NfFYu/nVRxvwwRlrObCmvrn81Q7uTpTc4H1KsQIbWgRBCsi/0JHI4xfyLSb/q/eTBag/KxhMbjWeBhXrTd06X6kHVzcyoeyvcK4X8bmzq4snBeUKId0FPIhd5F3VuqyOvN3ldVvRaISPrj5QrCVfkpXUvSYsZLUy67bTrNiohhOQM6EnkQormLyp9qvaRB6o8IJtPbjbexcK9C01atXhV411EVYiSwvkKe7qqhBAP4+zOdD1V27A3tSoynbK+qo2qhSnyn1LttO4z1sorr7qk2mTpQ2fqSVL9TKRWyVry6t2vyvJey2X0XaN1MDpRXln/irSY2UJe+OEF2XJyCwPDEZKLcdaTwLaj3VSTMlD2adV2VYBdI3WPJp1VNZOSkq7o65J25fdqXi0n60cySEC+ALmvyn3S+/besvXUVrNIb/G+xSaMeeVilY130T68vSlHCMk9OGUktBFHo59unHo9j7CD7VWvqZ6xO/WE6k0YCOt+J5ypjzNwJ67r4LO8I/gOoxGRI4yhQDcU1lu8E/uOtCnfRnpW7ilJ+s9H/xFCcjZZNSYxQRWtKpIiv7KqiTZMMB6XVSPUUPxknauA7ilNL6ie1/y1jm6sZQZoAklYGENpuxKMSfS6vZcRBrVhLBbHLzarufPnKyOBCU3k/JVqZoyDEJJLxyS0EY5RbXUgdBOli5broMkJbeTjUjFSmEPZQPWsaoaWx8/To6owvaa25Xl8pdkO+zm0zEeqSCg4ODgjVSK3AKbLvtDwBVnZa6W82PBF9SHyyXG/6WZmFGJHxR2P49gFIbnRk9DGt5WT79FY1Ukb+ShNsSVSgB5P0fv21eNDqjl6jCAPGzAArmkJfYndUWxdUHGav9fyOmKdrAtxEqzWRojyGSvLyGWfg1Kv5i5ZFL9IFsYvNGHLsaobgQYD/QP5rAnJAbh9nYQ28qNUIary+vI+1QrLQIB5qhY4UEMAI4CddE7pcTBmQ1n54ZpEqOLdXVeSOfyTQuX5Bs+bVd0vN3rZDGqPix1nZkZh34ufjv1E74KQ3DwmoQ14V03eVaGfZxGmq6oBaKtpGX09WY/hPaTFpxC6rzT9S9UfXoW+bqrHL2t6TdME1SDNPuNMXYl7vYuuEV2NEExw9q7Z8nX817Jk3xIpH1BeukV0M95FUAHGeCIku+FzvacnZxAZGZkUG+v5Hqnek9aZdPrAhh6uief+xsvXLsuyA8uMwfj5xM9mH29su4qptPVL1TfRawkh3oH+II/DuK6jc1xxTdyCf15/4z1Ae8/tNTOj4F0s3b9UQgqHmHGNLpW6mH0xCCHeC3/OEbdTMbCiPFf/OTN28WaTN6VUoVIy8eeJ0npmaxm2cph8f/h7s9KbEOJ90JMgWQb24MaqbWjf+X1mNTe2YI35LUbKFi5rxi7gXWBPb0KId0BPgngETJcdHjlcYnrGyLim40wX1Lsb35U2s9rI0yuelrWH1kpCIuYsEEI8CT0J4lHy+eaTdhXaGf124TcTM2rennmy4uAKKV2o9PVZU5W6mi4qQkjWQ0+CeA1hAWEyrO4wiekRI283e9tMn/3Ppv9I29lt5cnlT8qqg6vkWiJmRRNCsgp6EsTr8PP1M4EEoYMXD5qd9ObumSurV6w24xUYu+hWqZuULlza01UlJMdDT4J4NaFFQmVInSGyrMcymdB8gkQUi5BJmycZ7+KJmCdk+W/L5WriVU9Xk5AcCz0Jki3wy+MnLcu1NDry+xEzMwoextCVQyW4QLCZFQUPI6QIotITQlwFPQmS7ShTuIw8WftJWdpjqfz7nn9LtaBq8snWTyRqTpQM/HagfHvgW3oXhLgIehIk24JQH/eE3WN07I9jxrOYs2eOPLPqGSnuX9x4F90jupsBcULIrUFPguQIMEX2iVpPyDfdvpH3W74vdwbfKZ9v+1zaz20vjy17TL7Z/438lYAYkoSQzEBPwg3k5MB+3o5vHl9pGtLU6Pgfx82aC4xfPLv6WSmWv5h0rtTZeBfliyJyPSEkPehJkBzLbYVuk4F3DpTF3RbLh60+lLq31ZX//fo/6Tivozyy9BGzFeuVBLO3FSEkFehJkFzhXTQu29jo5J8nZf7e+SaE+XNrn5PADYHSsWJHs6NeeCD2tyKE2EMjQXIVwQWD5bE7HpNHajwiPx790YQwn7pjqvEw6pSsY/a7aF2utQl1TghxsrvJx8enp2ob9qZWRaZT1le1UbXQLm86drOztB+p3blRqj2qnaq2/LCIK8GmRw3LNJS3m79twoA8U/cZOXXplPzzu3+a7Vff3PCm7D67mw+d5Hqc9SSw7Wg31aQMlH1atV0VYMtISkrqbTtWQ/C2Juet42rWftjVVdgKNQZ7YGt5hgUlLgfbqj5c42F5qPpDZl/uWbtnyYydM+TL7V+aWVLwLtqWbysF8hbg0ye5Dqc8CW20t6t2pldOG3gsg22vmpzKeR9NeqmmWlmdVdP03ldU+/R4j6q+M3UlJAPfU6lfur6MbTrWbJA0InKEXPjrgvzf9/8nLWe0lNfWvyY7z6T7dSckR5FVYxITVNGqIqmcb6I6rgbB5t+XVa23O3/IynP0H3uAJpCEhXHRVFaR06f5FvMvJv2r95d+1fqZPboxdoGptNN2TpM7StxhvIt25dtJQb+Cnq4qIZ71JLQRRlfPVgfCr/100XIdNDmhBiAujWL323kR5jIHZZIcXaj3/QgbeEPBwcEZqRIhmfIuMHX2jSZvyIpeKyS6XrT8efVPeeGHF8zYxcvrXpZfT//KJ0pyryehjW8rJ9+jsaqT/meL0hRTRgL0eIrety9O6nFea1yjbgrPIdTuNbqrjjhZD0Kcomj+ovJgtQelb9W+svnkZpm5a6Ys2LvApIgfBe8iqkKUFPIrxCdNcgw+2lg7fxMfn1WajNB7xaZTrrlVroNdXjtNRmleM7s8DFh/ZY1DYOB6uSoivYHryMjIpNjYNKtAiEs5f+W8LIpfZAa7MRsKg9swFDAY1YOqG0+EEG9Hv6dx6I1x+ZiE3rirJu+q0M+zCFNY9Y3aaoqGfbIew3tIj/tSdDXBe8G02hl6CD8eW5EN5swm4q3exQNVH5D7q9wvv5z6xXgVi/ctNtuw3l7sdmMs2oe3lyL5UhuOIyQXeBLeAj0J4g1c/OuiCfkB72LHmR3i7+tv9vBGzChMqaV3QbKTJ0EjQYibwA8wDGrbvItL1y5JpcBKxrvoEN7BeCGEeAM0EoR4mD+u/iFL9i0xU2m3nd4m+X3zS5tybYzBqF2yNr0L4lFoJAjxIuBdIMDgon2LjPEILxpujEXH8I4S6B/o6eqRXIgPu5sI8T6w3mLp/qXGu9hyaovky5NPWpdvbcYuIm+LpHdBsgwaCUK8HIT7wIyohXsXysWrF6V8QHnjXXSq2Mms/ibEndBIpMV3E0TK1hGp0PRG3r41Iod/Frl7qDs/F0L+Bga3l+1fZga7sWDPL4+ftAxraQxGvVL1TPRaQlwNjURawCDMfEik53+vG4qUrwnxEHvO7jHeBVZ1I9BgWJEw6V65u3Su2NlEriXEVdBIpIfNMEQ+KhL7CQ0E8SouX7ssMb/FmLGLuONxktcnr9wTdo/xLhqUbkDvgjgNjURGWPGayJqxIk2jRVqMdvqhE+IO4s/Hm5lR8C7OXTknZQuXNQPdXSp1MbvuEXIr0EikBz0Jks34K+EvWf7bcmMwfjz2o/j6+Erz0ObGu2hYuqHZ15uQjEIjkRYckyDZnAMXDpixi/l75suZy2ekTKEy0jWiq3St1FVuK3Sbp6tHsgE0EmnB2U0kh3A14aqsOLjCeBfrjq4zYxVNQ5pKz8o9pXGZxlm9n4MAABObSURBVPQuSKrQSBCSyzh44aDM2TNH5u6eK6cvn5ZShUoZz6JbRDdzTIg9NBKE5FKuJl6VVQdXGe/ihyM/mFXcd5e92wx2w8vImyerdjAm3gyNBCFEDv9+2OzTDe/i5KWTUrJASekS0cV4F5glRXIvPozdRAixcS3xmqw9tNbsd/Hd4e9MSPNGZRqZmVHNQpuZVd4kd+HjLiOhN+6pyYuqqqr6aW1fqmUxJw/nD9u2L9W86ZrcbhVB+Mtzeq6W5pfX4+2qnda59Zo/KL36cNMhQjLH0d+Pytw9c42HcfzP4xLkH2TWXKA7KjTAfpt5kpNxp5GAcUhUTUpvj2st+4wmqESA/R7Xduff1uS8nnvZMhIL9bhGZupDI0HIrZGQmGC8CngX8DISkhLMam54Fy1CW4ifL72L3GoknBq10pvi1366IY31fIgm7VWvqZ5xcB436KVq4Ux9CCG3BhbfoasJOv7HcZm3Z57xLkasHiHF/YtL50qdjXdRLqAcH3EuI6tCSk5QRVtehyOaqI6r0dltl1dBbcdG1WoVzhNCsgAswBt450BZ3G2xfNDqA6lTso58se0L6TC3gzy69FGzwx5WfBPvofekdUbuIF1PQhvoGE0cTawerY36/Axcj66lE1o2To+bp1LsftVUu9dHVWF6zWm9pq4ez9O0ur6+4OD+AzSBJCwsLL3qEEIy4V1guix08s+TMn/vfBNkMHpNtATmDzR7XSAqLXbWIzkXp8Ykkm/i47MqtTEJPfeGJg+qrqn8VQGqOVq2r3Uehuqwqq7mHcrs/e3hmAQh7iUxKVHWH11v1l2s+G2FXEu6JnVvq2vGLlqXa2327iZZj82LmD6wocvHJNze3aRvPEoVosJg9H2qFTYDYdFKtcPeQGiFg63ZUDjGz5QIVby760oISRuE+sB02bebvy3f9vxWhtUdZryMUWtHSYsZLeStDW+ZfTBIzsEpI6ENeFcVGneYr0V6vNTKL6NanMHb3Jeiqwlgt58teo/Nms5SDVIjcsaZuhJCXEuJAiXkkRqPyNddv5bJbSYb4zFt5zTpuqCr9FvSzwQcxE57JHvjku4mb4HdTYR4FkSh/Xrv12bsYv+F/VLEr4h0qNjBdEdVLlaZH0827G5i4BZCiMvAdNn+1ftLv2r9JPZ4rDEWGL+YumOq1CxR0xiLtuXbSkG/gnzq2QQaCUKIy8HSp3ql6hmdu3xOvo6/7l2M+WGMvPXTW9IhvINZd1E1COtxiTdDI0EIcSuB/oHyYLUHpW/VvrLp5CZjLLBYb/rO6VI9qLrxLu6tcK8U8ivETyIXL6YjhORy4F3ULllbXrv7NVnec7mMrD9SriRckZfWvWRmRr34w4uy7dQ2E3CQeA/0JAghWU7R/EWlT9U+8kCVB2Tzyc1m+9XF+xabtGrxqsa7iKoQJYXzFean42HoSRBCPOpd1CpZS15p/IrxLp6/63mzYO+V9a9Ii5ktZMz3Y2TLyS30LjwIPQlCiFdQJF8R6V2lt/S6vZdsO73NjF3Au0Aoc0yfxUA3ptMG5EPQBpJV0JMghHidd1GjRA15sdGLsrLXShnTcIzZZvWNDW9IyxktZfR3o2XTiU30LrIIehKEEK8FM556Vu5plOxdxC+WBXsXSKXASsa76FixoxnjIO6BngQhJFuA6bIvNHzBeBcvNXpJCuQtYNZcYGYUYkfFHY+jd+EG6EkQQrIVWK3dLaKb0Y4zO4x3sSh+kSyMXygVilYw3gXCmBfzL+bpquYI6EkQQrItVYpXkecbPG9mRmGGFAa1x8eOl5YzW0r06mjZcHQDvQsnoSdBCMkR3kWXSl2Mdp/dbdZbYNxiyf4lZstVm3cRVCDI01XNdtCTIITkKCKKRZjV3Ct6rpDX735dgvyD5J24d6TVrFYyfNVwWXdknVmLQTIGPQlCSI7EP6+/mfkExZ+Ll1m7ZxnvYtmBZRJSOMRsvQrPA/tikNShJ0EIyfGEB4ZLdL1oM3bxVpO3pHTh0jLx54nSemZrGbZymHx/+Ht6F6lAT4IQkmvAHtxR4VFG+87vkzm755gd9GJ+i5GyhctK10pdpWtEVylZsKSnq5pjti/tqdqmSlRFplPWV7VRtdAur5ZqvWqTKlZV3+7cKNUe1U5VW2fqSQghKcF02eGRwyWmZ4yMazZOQoqEyHub3pM2s9rIkBVDZM2hNZKQmJDrH5yznsRWVTfVpAyUfVq1XWUfeGWs6qWkpKQlagiirNfN9biatfd1dVUZVYzmVdZy/MQIIS4ln28+aVe+ndHBCwfNzCjEi1p5cKWUKlRKulXqZrwLHOdGnPIktNHertqZXjlt4EM0aa+anPIWdkYD6+qPWMedVdP03ldU+/R4jyrZyyCEEHcQGhAqQ+sOlZgeMfJO83ekYtGK8sHmD6Tt7LYyePlgWfnbSrmWeC1XPfysGpOYoIpWFUmRP1S1VI3IeMtgNbLyy6rW25U7ZOX9Db12gCaQhIWFubDKhJDcip+vn7Qu19ro0MVDZuwCu+kNOTREShYoKV0iupi1F2UKo6Mjl3sS2gijq2erA+HXfrpouQ6anFCPIM7B6SdUw/RcKFLVJ7bLHJR1uF2VXvuRKhIKDg7OSJUIISTDhOhYxZA6Q2Rpj6Uy4Z4Jcnvx2+XjLR9Lu9ntZFDMIFl+YLlcTbyaez0JbXxbOfkejVWdrDEHf1WAHk/R+/bV4/7WWAWYadcdBc8BhsNGiF1XFCGEZDl+efykZVhLo6O/H5U5e+YYD2PoqqFmrQXWXCCeVGgR+6Yr++P2dRJqDEapQlTlrcHoFZaBEKvhb2Ydt1Dtto4XoKwak/yqCnocodrg7roSQkhGKF24tAyuNViWdl8q77V4T2oE1ZBPt34qUXOiZMCyAbJs/zK5mpAzvAunxiS0Ae+qybsq9PMswlRWNQBtNUVH3WQ9hveQFo+rJmp51OOybWxBr8O02hl6+KsKo0SDObOJEOJt5M2TV5qFNjM69scxMysK3sXw1cOluH9x6Vyps/SI6CFhAdl3vNRHG19P18FlREZGJsXGxnq6GoSQXExCYoJ8f+R7mb1rtqw+tFoSkhLkrlJ3SY/KPaRFWAsz5dbV9J60zqTTBza8pev1R3kcxnUdneOKa0IIcSG+eXylaUhToxN/njCzouBdPLvmWSmWv5iJRou4UVjMlx1g7CZCCHETJQuWlAE1B8jibotlUqtJElkqUr7c/qV0mtdJHv7mYbNZ0pWEK179/OlJEEKIm8njk0calW1kdOrSKeNdoDtq5NqRUnRDUekY3tF0R1UMrOh1nwWNBCGEZCEldLrsY3c8Jo/UeEQ2HNtgtl+dtnOaTNk+RWqXrG2MRZtybUyoc2+ARoIQQjzkXTQo3cDo9KXT8vXer82eF6O/Gy1vbnhTOoR3MAajcrHKHv18aCQIIcTDBBUIkodqPCT9q/eX2OOxMnPnTONhTN0xVWoG1zTTaNuWb2u2ac1qaCQIIcRL8PHxkXql6hmdvXw22bsY88MYGfvTWGkf3t54F1WKV8myOtFIEEKIF1LMv5j0q95PHqz2oPx84mfjWczdPVem75xuVnjDWNxb4V63exdcTEcIIdmE81fOy8L4hcZg7Dm3RwrmLWi8i9hfKkuBpHJcTEcIIbmZovmLSp+qfeSBKg/I5pObjbFAl9Tl/JclIAFb7tzaiuu0YHcTIYRkw7GLWiVrGUXXj5buX7wvvkmF3fJeXHFNCCHZmIB8AeL3R2P5/Ux1iTtw1uX3p5EghJBsTJwahh3HLsqhs5ekz+T1LjcUNBKEEJKNWR9/WhKtYN5XryWa166ERoIQQrIxDcKDJI+14bNf3jzmtSvhwDUhhGRj6pYrJlVKFZELl6/JxPtqm9de40noCHtPFXaRS1RFplPWV7VRtdAur5ZqvWqTKlZV38ovr7pk5UMfOlNPQgjJyRTx95OygQVcbiBc4UlsVXVTTcpA2adV21UBdnljVS8lJSUtUUMQZb1ubp3bq/m1nKwfIYQQT3kS2ohvV+1Mr5wagBBN2qsmp7yFndEoqjriTH0IIYRIthyTmKCKVhVJkT9UtVSNyHjLYDWyO1cB3VOaXlA9r8ZoraMba5kBmkASFpZ9NxsnhJBs6UloIxyj2upAnTPyBlqugyYntJGPc3D6CdUwPReKVPWJlX9UFab5tTV9RvWV3se+myoZLfMRNvCGgoODM1IlQgghrvIktPFtlcF7pUZjVSdrzAFbLQXo8RS9b1897m+NVYCZtu4oPYdNX83GrzAuWn6vHmLnjVgn60IIIcSb1kloIz9KFaIqry/vU62wDIRYYxDNrOMWqt04UKMQjNlQ1nG4JhGqeHfXlRBCiAvHJLQB76rJuyr08yzCdFU1AG01LaOvJ+sxvIe0eFw1UcujHpdtYwtKU9XLmn9N0wTVIL3XGWfqSgghJIuNhDbcczWZ6yAfHkKUg/xVmqyye/2dJnUdlJutCUQIIcSDMCwHIYQQGglCCCGZh54EIYQQGglCCCGZh54EIYQQGglCCCGZh54EIYSQVKGRIIQQQiNBCCEk89CTIIQQQiNBCCEk89CTIIQQQiNBCCEk89CTIIQQQiNBCCEk89CTIIQQ4h4j4ePj01O1TZWoikynrK9qo2qhXd6dqnWqX1RfqwLszo1S7VHtVLV1pp6EEEI840lsVXVTrclA2adV21PkTVaNTEpKusPa4e5ZZKpRqGbth11d1U71H9ue14QQQrKJkdDGfbtqZ3rltIEP0aS9ZRTsud3OwHyr6m4dd1ZN03tfUe3T4z2q+s7UlRBCiPeOSUxQRasSHXginazjnqpQ67is6qBduUNWniMDNEAVC508edJ1NSaEEJK+kdDGN0a11YHwaz9dtFwHTU6oRxDn4PQjqsFaBueKqP6yXeagbJKj++t9P1JFQsHBwfxICSHEheRNr4A2vq2cfI/Gqk5qCKI09VcF6PEUvW9f1Q593QaFNK+y1SVl8xxsXgVAd9URJ+tBCCHE27qb1BCMUoWoyluD0StgICzDUNJKUY/nVR9aly1AWc3Pr6qgxxGqDe6uKyGEENdOge2qwq/+hqpFerzUyi+jWpyBW9yv5XZpusPyFD5DphqRbZrMUP2q+kY1WPMSnKkrIYQQN3Q3pYU23Ji2OtdBPhr8KAf5qzRZZfd6oiYTU7n3a5pAhBBCPARXXBNCCKGRIIQQknnoSRBCCEkVGglCCCE0EoQQQjIPPQlCCCE0EoQQQjIPPQlCCCHuWUxHCCHE80wfiKAX7oGeBCGEEBoJQgghmYeeBCGEEBoJQgghmYeeBCGEEBoJQgghmYeeBCGEEBoJQgghmYeeBCGEkFTxSUpKSvVkdsPHx+ekJgecuEUJ1SkXVceVsF58Xvx+eQ8lcmA7UU5tQXCONxIuMDKx+jwiPV2PlLBefF78fnkPPrmsnWB3EyGEEBoJQgghmYeexM18lPlHmCWwXnxe/H55Dx95ugJZWS8fjkkQQghJDXoShBBCUoVGghBCSO4zEjodrJ1qp2qPaqSD8+Df1vktqjrpXavHxVXfqnZbaTEvqdeLqsOqTZaisrhen6pOqLamuMbTzyu1ennseWkaqlqp2q7apnraG55XOvXy5PPyV21Qbbbq9ZKXPC//NOrl0f+P1nlf1UbVQqefF8YkcpoUX9VeVbgqn2qzqlqKMvjgluDZqRqofkzvWmWsaqR1jA/uLS+p14uqEZ54Xta5pip8SbemuMZjzyudennseSmlUSfruIhql5d8v9KqlyefF14Xto79kI/zXvC8fNKol0f/P1rnn1F9pVro7P/HnOpJ1Fft0T8wXvWXHk9TdU5RBq+/0PNgvR4HqmUtnc61SD+3jpF28ZJ6OYsz9cIXaY0mZxzc15PPK616Ocst10uPj6p+tup3UZPtqrKefl7p1MtZnKkX+N0qg8YYSvKC55VWvZzFqe+9piGatFdNdnBNpp9XTjUS+HIftHt9yMrLSJm0rr1NP5CjOLDSkl5SL/Ck5XZ+egtutzP1SgtPPq/08Pjz0vctr0lt61eo1zwvB/Xy6POyuk426eEJ1bf6bLziefmkXi9Pf78mqKJViSmuuaXnlVONBFywlKS08qmVyci13lavD1QVVbVU+PDfzsJ6uRN31cvjz0sbjsKazFYN1f+wFzL5/lldL48+L61HggrvjV/I9bWONTL5/lldrw889by0Dh00PaH1isvke+Y6IwGrGmr3Gh/ikQyWSeva43YuXWnrF4TH66VfiOPWFxa/HD623NWsqldaePJ5pYqnn5c+Cz+rIf5S6zDHW55XavXy9POyq8c5TVap2nnT9yspRb08/Lwaqzrp89hvdVO10OMpTj0v/UNynJS8qnhVBbuBn+opyrRPMfCzIb1rlXEpBn7Gekm9SttdPwxfjqyql9358g4GiD32vNKpl8eel/X6C9UEB/f15PcrrXp58nkhMmmgdVxAtVbVwQueV3Aa9fL4/0erTPMUA9e39LwyXPHsJmv0f5c1S2C0lTcIso7xcN+3zv+iikzrWis/SLVctdtKi3tJvf5nld2iWmD/Jc2iek213Oqr1i+cR73keaVWL489L+VuVZL13ujPhqI8/bzSqZcnn1dN1UbrvTGVeYw3/H9Mp14e/f+YhpG4pefFsByEEEJy3ZgEIYQQF0AjQQghhEaCEEJI5qEnQQghhEaCEEJI5qEnQQghhEaCEEJI5vl/p1ItElt6HaIAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.errorbar(dtaus**2, ener[:, 0], ener[:, 1], fmt='.', label=\"Monte carlo results\")\n", + "plt.plot(0, E, 'x', label='ED result')\n", + "\n", + "def func(x, y0, a):\n", + " return y0 + a*x**2\n", + "popt, pcov = curve_fit(func, dtaus, ener[:, 0], sigma=ener[:, 1], absolute_sigma=True)\n", + "perr = np.sqrt(np.diag(pcov))\n", + "xs = np.linspace(0., dtaus.max())\n", + "p = plt.plot(xs**2, func(xs, *popt))\n", + "\n", + "plt.errorbar(0, popt[0], perr[0], label=\"Extrapolated value\", color=p[0].get_color())\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.05\n", + "Compiling ALF... Done.\n", + "Prepare directory \"/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.05\" for Monte Carlo run.\n", + "Create new directory.\n", + "Run /home/jonas/Programs/ALF/Prog/tV.out\n", + "0.1\n", + "Prepare directory \"/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.1\" for Monte Carlo run.\n", + "Create new directory.\n", + "Run /home/jonas/Programs/ALF/Prog/tV.out\n", + "0.2\n", + "Prepare directory \"/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.2\" for Monte Carlo run.\n", + "Create new directory.\n", + "Run /home/jonas/Programs/ALF/Prog/tV.out\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.05\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.05/Ener_scalJ 1\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.1\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.1/Ener_scalJ 1\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.2\n", + "Analysing Part_scal\n", + "Analysing Kin_scal\n", + "Analysing Ener_scal\n", + "Analysing Pot_scal\n", + "Analysing SpinZ_eq\n", + "Analysing Green_eq\n", + "Analysing Den_eq\n", + "/home/jonas/Programs/pyALF/Notebooks/ALF_data/tV_L1=4_L2=1_Dtau=0.2/Ener_scalJ 1\n" + ] + } + ], + "source": [ + "dtaus = np.array([0.05, 0.1, 0.2])\n", + "ener, E = extrapolate_dtau(\n", + " 'tV', # Hamiltonian\n", + " { # Model and simulation parameters for each Simulation instance\n", + " 'Model': 'tV', # Base model\n", + " 'L1': 4, # Lattice length in the first unit vector direction\n", + " 'L2': 1, # Lattice length in the second unit vector direction\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", + " }, # HS field couples to z-component of magnetizatio\n", + " dtaus,\n", + " do_ED=True\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2dB3hUZdr+n1RCCYQaIIQSem8RUAQpoiyrYsOGiq6uqCuW/6eiIrv+FXex7Vo+QRFBQFBwpVlQEQRsAUG6gEgE6U1KKIEA+Z77yZlxCJM6M5lJcv+u677eM2fOOfPOYTh3nrc8b1hmZqYQQgghBSG8IAcTQgghNA9CCCGFgpEHIYQQmgchhJDAw8iDEEJIgYks8BnFlGrVqmXWr18/2NUghJBixbJly/bpqNzqpdY8YBxLly4NdjUIIaRYERYWtsXbfjZbEUIIKTA0D0IIITQPQgghgafU9HkQQs4lIyNDtm3bJunp6bw9pZyYmBipU6eOREVF5et4mgchpRgYR2xsrA0o0Y7RYFeHBAnkONy/f7/9Hho0aJCvc9jnQUgpBhFH1apVaRylnDD9wwG/g4JEoDQPQko5eHAQElbA3wHNgxBCSIGheRQh17/5vYkQ8gcRERHSrl07t0aOHGn7e/ToIU2bNpU2bdpIs2bN5L777pODBw8W6a2rUKGClZs3b5YpU6bwn80DmgchJH9887LIr4vO3ofX2O8DZcuWlRUrVrj12GOPud+bPHmyrFq1ylSmTBnp379/ntc7ffq0T/XxBs3jXGgehJD8kdBB5IPb/jAQlHiN/QEmOjpann/+efntt99k5cqVXiOEv//979K5c2f5/vvvkY9JLrroIunYsaNceumlsnPnTjvu1VdflRYtWlg0c8MNN9i+p556Sl588UX3tVq1amVm4QkM7euvv7bI6D//+U8Av2nxgUN1CSH5o0F3kQHvZBlG8h0iS9/Oeo39PnD8+HF7KLt4/PHH5frrr/favNW2bVtZv369lZ4cPXrUHvpPP/20zV2BccyaNUuqV68uU6dOlWHDhsm4ceOsSezXX3+1KKYgTWA4Dwbz8ccfF/6LljBoHoSQ/AOjgHEsel6k+6M+G4dns1V+5yN4A8ZyzTXX2PaGDRtkzZo10qdPH3czVq1atWwbEcfAgQPlyiuvNJHCQ/MghOQfNFUh4oBxoGzQzS8Gkh9gAqtXr5bmzZt7nR0NA3EZTMuWLa35KjuffPKJLFq0SGbPni3PPPOMrF27ViIjI+XMmTPuYzjbPn+wz4MQkj9cfRxoquo17I8mrOyd6AEATVFozkpMTLToITcwQmvv3r1u88C5MAkYxNatW6Vnz57Wf4JmqyNHjtjs+h9//NGORYlmrexgFn5aWpr/v1gxhuZBCMkf2388u4/D1QeC/X7o83DJc7QVmphgFujPQL8G+jHy07n+3//+V4YOHWp9I7jmd999Z5HLzTffLK1bt5b27dvLQw89JHFxcdbc9fvvv9txo0ePliZNmpxzTdQBEQquxw7zLMJyakMsaSQnJ2cGezEo1xyPqYPPD2o9CHGxbt06r81ApHSyzsvvISwsbJn6RHL2Yxl5EEIIKTA0D0IIITQPQgghgYeRByGEEJoHIYSQwMPIgxBSIJgdmtA88sHtn91uIoQEbhGiW265xf361KlTlpPqsssuK9T1MPlv1KhR/qpevti8ebPNRQk2CxYscN83bGN+S6Bg5EEICSrly5e3XFSYLAjmzp0rCQkJhb5eUZvHKTU7XwhECnlA8yCEhBRp6Rmy/eBxWbblgN+u+ac//cnyToH33ntPbrzxRvd7mP2NJIaY5d2lSxdb28OVSv0vf/mLLRqVlJRk6dYBZqhv2rTJZow/8sgjtu+FF16Q8847z67xj3/8w2sdPvvsM+nQoYPNIu/du7ftW7JkiVxwwQU2Ix0lki6Cd955RwYMGCCXX365XHLJJWddB7mxbr/9dvdM9q+++srrgx1pUm666SY7Drz77rvSqVMnq/fgwYPNVKDbbrvNohoc55rdju/smvS8b98+S7GSPRJ644037HhcD+nkP/jgA7sOvl/37r7nI2NiREJIvoFhrN+VJmcyRQaOTZHJd3aRjvUq+3wHsbYG0qmjyQXmAFPAAw/gYY+H8MyZM2X+/Ply6623urPwIj07Hs7IO4WcVvfcc4+lT0ck4zrmiy++kI0bN5oRIKPGFVdcYckRPR+gyIX117/+1fY3aNDADAtgBUPsQ2qSL7/8Up544gn58MMP7T3kzkJdq1Spctb6H6+//rqVSOKI+sFcfv75Z0ve6Anqg3ri8zCzG6njv/32W4mKipJ7773XFsJCgsft27fbcSC/aeRhJnfffbetc/Lwww/bPpjP559/blGdP1ZkpHkQQvJNSup+Mw6QceqMvfaHeSAiwAMYUUe/fv3Oeu+bb75xP7B79eol+/fvl0OHDtnrP//5z7Y2B1SjRg3ZvXv3OdeGeUAwIIBkiDATT/NISUmx13iQAxgCwOcMGjTIjkffDJIsuujTp4/7uOz1HTJkiNt86tWrZ+aRPaEjogzX582bN88WsEJ0BNCEh++DyCY1NdWuh++aPcopCF27drUo5rrrrpOrr7660NdxQfMghOSbLklVJTxMzECiIsPttb9ARIC/ktGkA4Nw4S3/Hh7kAKbhAinZvfU/4Hxk5EVTUE7gGNc1PRk+fLg1L82YMcPMDc1Fnn01OV0rP3iej3NgUv/617/OOQ4rJyJiQEQzbdo0W9TKM418flPIoxlr8eLF1jyIpixEZlWrFv7fjx3mhJB8gyijWc1YqVO5rN+arFygqQpLybr6AFwgIkATDoCxVKtWTSpWrJjjdbKnT8cytHjgIuIAaAbas2fPWeecf/75snDhQnc6dlezFSIPV+c9+jnyQ3eP+iLiwNK5aFLLDfSxIBOwq174/C1btlh/BkwCmX+x/ogrdTyapRCpAJyXn/uAfiAs04vmQdxDpKf3BUYehJACERsTZfKncYA6derIAw88cM5+dIyjAxrNPuXKlZMJEybkeh38NY0mGnQOoyMeneXoU4BBAPQDoHMazUIuMDR4zJgx1pyDhzXew6ivRx991CKCf//739Zklh/u1f4K9DfABBEhwHQ8IyRvYF31ESNGWLMUPh/9Hog0sMoivrsrynBFJojQ0Pw0adKkHOuFJq9rr73W0ti/9tpr1nmO5jdEOTCr7Ev5FhSmZM8D1xyP8X3H+3SjAVOyk+KQgjsv+DsuuRQkJTsjD0JIgeB6NMQvfR7qSkNUG1RrVc/ncMxm1WrVCpV7RSbdrqKaq9rolBYHa9kHbuecg9Idl+l2R2f/L6pXw7z1chFCCAld89Dndk8t+qvaaFjTUssXczm8px7TLlv4g/Um5+m+xiid12Cf6nLdj56zQapJHueMVt2lwjlQX1++AyGEkKKPPO5RjdSH/Am80PLsIQx5A+Nx9X6hvNK5znLVDmf/WlWMGlUZVS3drqjvfZ+ZNR5uouscQgghxcc8sFJ8N32oL1YtVGXNcDkXPOi/cJqgEDW4iFcP2GkHZJV/DH/4g2tUyx2Dwpi5bR7vYTvHJDj4LDSTQZhBSgghxD/k2WGuD94vtajp5a1hzvnop+iignFM0+OTnKjAk66IJPQ9mAP6Ntbr60X5+Gw0hT2nck2r9Na/keOMHP2MMVpAkpycnL+ZO4SQIhuBSEpw5KEP4ItVrbxolvOX/3SYhbJEtzEYuZqXa+zwaNaaoerkvLXbaYqCUaB0N3vp6zrOsbfqeZuc3fg87HeBbVfzFiGkGIKZ4Zjx7BJyU+XGP//5zyKpF+aXvPhibt24Yvm2fvrppwJfG3NN/IG/rhOMZquZql7Owx5NWNFOZ7cb3V9eFevadqKIrCxfIrOdDnGAcpZzXJwWSLH5uBrHt65rOU1bafp+F2eU1a2ucwghxRNMhEOqDJeQFbcw5oEGD9dkuqJiZiHNoyTgq3mMUyXpcxxm8L5qEEIQfV1b9alzTLzqG329UktEJ5/oIZ857+FPDAzL3YjSeQ3uUzVSDXeG90I1PDrpx6p+USEimePjdyCEhBhIC4KUHq4U6EjR/tZbb5mxIGkgIpSBAwdavilMasOsbqRTR8oNZNbVZmrLSOuZfh0pPYYOHWoJCaFffsEjRCwNCGZcYwY7SqQTyQ4+G0kLMSsbqUKOHTtmCy3Nnj3b0r6jPkj/AfXt21c6duwo3bp1s6y6AGlPMMMd10C+LG+gbp7rkCDyeemllyytCuqF74dZ65gxntsiUOC+++5zp1NBGpOLLrrI6oRULTt3Wjezz/g0SVBN4KQWN3vZj6akfs52qhZe58Hre8h+1tvL/hFajMjhHMwTCf6SXYSUMJ5b8pys/z3rYZcbrmPys8JmsyrNZGinobke4zIDF0hieP3118v//u//WhZYpCw5cOCApUwH2O9Ktw7zgMGMHz/e/eB99tlnLdst1sLAQxdp010ZbZETC6nQJ06cKA8++KB8/PHH9qBFmnekIUEOrPvvv98iCk+QtsT1+U8++aS8/fbblukWyRzx0EYaEIDPQwLCxo0bWxJCmBrSyOM7wNTwOa6U7d7S0qNOOAcgCSLWGEEqdyRmRN2R6wprmuBz8zPFDVmAUU8YDlKwIO37sGHD7Hv6CmeYE0JCotkqO0h5jgWM/va3v1lm2ZxAynM8UF3goYs8Vciwi7+y0azkMg/XIlMoH3roIfe6HNOnT7dtLIeLfFbZwXoaMA2sg4FI4FL9Cz472I9oBItEuThxwmYx2DodrrTy+AxEGdlByngkRtyxY4etL1K5cmWpW7euGQDWEcG6IuHh4ZbYEanna9b0No7pbGCsqDvuJYCh1qpl3cw+Q/MghBh5RQhFPdoK/RfItQRzQZZZJE7MK7U5mofQyf3DDz/YwxeRi2fKcs+/1nP6y93bflwH0QiardActECbibzVNy4uzqsR5vZ5niCCQZbcXbt2WSQCkKEXZoLmJyRMRPNb9jTsninaget99AOh+Q4G6W+Ykp0QEpIgCyz6M7BAFNK1uxZiwgPUc1EmTw4fPmxmUqlSJfvrfM6cs7tE0WzjKl1ZdrG87Pvvv+9+UF944YXnXBepzfEXOz7XlW49e9pzNCthcSdES64HtytiQpZfz8/ICRgGjoOBuJrC0P+DLL/43lg1EX003qIvRFiIdHA8FpcC6DeC8bjMA/Vfuxbzrn2H5kEICSquPg+X0CmOdTDGjh1rHcboeMYaGUhZDu666y5rhkKHeXYQGaD5B39tw3Dw0PYED1esafHKK6+41wPH2ufoM8E1keIc72UHa2ngPDT/YHVAz4c9Ur7jM9FZDmNAfwjqgTq4OrdxTfR1oMPctQqiN3AOzAhriLial/A9sV45BgHg+p6f7yIxMdFStLvui2vVxOjoaDMiNJOhTri/aFrzB0zJnhPfvCyS0EFu3zD+jxD910Ui238UufDBQt1sprImJSEle3GdJIjmHjyEsRAS8Q5TsvsDNQ754DZpVq+FrK9UPcs49LUMyN9qYoSUVIqbaZDAwA7znGjQ3Yzi3snXyPz4JJFVjnFgPyGk2IFhvcR/sM8jN9QoYBz9t+u49uQ7aBykRHJuKjpSGsks4O+A5pEb2lTVa3eqzErQDqqlb2c1XRFSgsAEtP3799NASjmZahz4HeD3kF/YbJUTTh/HqMadrc+j/8W3/9HnUcimq7T0DDmcfkqWbTkgHevZoomEBBXMndi2bZsN5ySlmxg1jpzm0niD5pETGFWlRrHeGW3l6gOx/YUwDxjG+l1pckYjw4FjU2TynV1oICToYO4A5iYQUlDYbJUTGI6b3STwupDDdFNS95txgIxTZ+w1IYQUV2geRUSXpKoS7mQniIoMt9eEEFJcYbNVEYE+jmY1Y63P45Ub2rPJihBSrKF5FCGxMVEmdpYTQoo7bLYihBBC8yCEEBJ4GHkQQgiheRBCCAk8jDwIIYTQPPzNjiM7ZPfR3VYSQgjJgkN180gWdujkITmacVQu/fBSaV6luVxc72K5uO7FkhSXlNuphBBSoqF55AIWrIdhpJ9Klz71+8i8LfPkteWvmRpUamAm0rtub2lRtUW+FrcvDXC1REJKBzSPfBATGSN/afUXE5qw5m+db0Yybs04eWv1W1KrfC0zEah9jfYSER4R6H83QggJKjSPAhJfPl5ubHaj6WD6Qflq61cy/7f5Mm3DNHl33btSJaaK9Ezsac1bnWt2lqiIqED8uxFCSFChefhAXEycXNX4KhP6Rb7Z/o1FJJ9t/kw+3PihxEbFSvfE7haRdK3d1V//ZoQQEnRoHn6ifFR5ubT+paaTp09Kys4UmffbPItKPkn9RMpElJHoqGYSe7qDHDrRQiqVqeSvjyaEkCKH5hEAoiOipXud7qbhXYbL8j3LZe6WufLBujmSFrFSekydJOfVPM+attDEVb1c9UBUgxBCAgbNI8BEhkeaUUArVlwk6WFbpGfH3RaVPJPyjIxIGSFtq7c1I+lVt5ckxiYGukqEEOIzNI8iJEzCpWxmA3mo403yYIcHZdPBTfLlb1/Kl1u+lBeXvmhqVqWZ9ZFgGHDDuIYcAkwICUloHkEC80IaVW5kurvt3bI1bav1jyAiGbVilLy+4nWpX7G+ewhwq2qtaCSEkJCB5hEioLlqUMtBpr3H9toQYEQk76x9R95e87bEl4vPiki0eQtzSdAcRgghwYJPoBAEHejXNb3OdOjEIVmwdYFFJBj+O2X9FKlcprL0SOxhRtKlVhfroCeEkKKE5hHiYEhv/0b9Tccyjsm3O761kVvQjF9m2BDh7gndpXe93tItoZuUiyoX7CoTQkoBNI9iBIyhT70+JswlWbxzsUUkaOKas3mORIdHywUJF1jzVo86PWwSIyGEBAKaRzEFTVXd6nQzDT+TNZcERoLRW2jmigiLkOSayTZqC0OAa5SrEewqE0JKEJF+GDU0RIv7VKdUn2RmZj7q5ZjNWqSpTuM4PSbZ2V9Fi6mq+iocc52+d0D399HtkSo05p9UPaL75zvnLNCiluq4c/lL9L09vn6P4gwSMcIooEfPe1R+2v+Tewjws4ufNdlcEicLcGJFziUhhATRPPRB3lOL/qo2+gA/oa9z+/O2px6zL9u+x1TzdP9IPRfb0FAVjrtc9+/Q/a10+3NVgsd5A/W9pb7UvSQPAW5ZraXpgQ4PSOrBVLeRvLTsJVOTyk2yjET7SRrHNeYQYEJIkUce96hGwjjwohARAIynh7M9QYWoYqheZ7nHMWtVMfpQLOP6HJJ/sGjVXXF3yV1t7pLtR7bbXBIYyeiVo2XUylE2RNhlJK2rtZbwMK5MTAgJvHk0UXXTB/uzWqarHtYH/A9ejstUfaHHoXxTjxnj7I/X7Z12gJY5RC7XqJZnM47xeiyawD5UjdD3cN1z0GPu0gKSunXrFuLrlSwSKiTILS1uMe07vs89l2TST5Nk/NrxUqNsDesfwRDgjvEdOZeEEFJ489AH8Jda1PTy1jDn/MqqLqrzVNP0+CQvD/OuThMUzGGuluv19aJ8fHZLLZ5TXZKtyWq7vhfrmMctqonezndMyowqOTnZq8GUVqqVrSYDmgwwHT55WBZuXWhRycxfZsr7G963IcIYsQUjOb/2+ZYVmBBC8m0e+gC+OKf39AGOZqvpjlks0ddntKym2pvtGjucco8eM0M3O6lgHrv1dS0n6kAnuLvZS1/X0QLH3qrvb/K41nanTNNjpjjX8moeJH9UjK4olze83HT81HH5bvt31k8CM5m1aZaUiyxno7rQvIUSc0sIIaUbX5utZqp6qRbog7yJMzrqrE5x3Y8nTbjzsC/vRBFPO2/PVg1yRlahnOWcgwkKn6ge1/O+9bgW6huHjnfdxhJ9l6kQGRE/UTayrPV/QBmnM+SHXT+4jeTzzZ9LVHiURSIwEsxyrxyDwJMQUtrw1TzGQfogX+MMqR2EKERf19btsbrZT8t41QyMAnI+b4ru/8w5f6TT1HWHlr+pBjj7MfS3kWq4vjfc2QfTOar63DGOCMc43vLxO5AcwBK6mHQIDes8TFbuXWlzSaBF2xZZ5zr6RlzJG2uW99a6SQgpifhkHmoCMIybvexHM1U/ZztVi7Y5nL9fi95e9o/QAvJGx8LWl/g2l6RDfAfTw8kPy/rf11tEgmV3Ry4ZacJorX0RjSX2THveakJKOJxhTgoMosjmVZubhrQfIqmHUrPSyauR7Imarh1X0+WqWROssx3NW5hX4kSehJASAs2D+ExSpSRJap0kd7a+U65+82M5HLFC4spskjGrxsgbK9+wIcIwEZhJm+ptOJeEkBIAzYP4lSipKlVP95bxfZ+U/cf3y8JtC20uyeT1k2XCTxNsiHCvxF7WIY+ledEBTwgpftA8SMCoWraqXN34alPayTT5etvX1k/yUepHMu3naRIbHSs9E3vaxMQLal9gI70IIcUDmgcpEmAU/ZL6mdJPpct3O76zUVvIADx702wzjgsTLrRRW93rdLfjCSGhC82DFDkxkTEWbUAZZzJk6a6lZiTodMciV1hiFyskuuaSIIIhhIQWNA8SVFyTDqEnOj8hq/ausj4SmMlT3z8l4Snhtma7K518rQpIREAICTY0DxIyYNJhuxrtTP+T/D/y84Gf3enkn/vhOVPLqi1t1BaiFozyIoQEB5oHCUkwL6Rplaamv7X7m2w+tNk9u/2VH18xwTwQjcBMmldpzrkkhBQhNA9SLKhfqb7c0foO066ju7ImJaqRvL3mbXlr9VtSu3xtdzr5dtXb2Yx4QkjgoHkUIVMHn1+UH1diQQ6tm5rfZDqQfsBGbMFIpm6YKu+ue1eqxFTJMhLtJ+lUs5Pl6CKE+BeaBynWIKvvVY2vMh05eUS+2f6N9ZN8mvqp/Pfn/0psVKx0T+xuRtI1oSvnkhDiJ2gepMRQIbqC9G3Q13Ti9AlJ2ZFiRoIVEz9J/URiImLMQNBPclHiRbaOCSGkcNA8SIkEKx/CIKBTZ07J0t1LLXGjq68kMixSOtfq7J5vgrQphJD8Q/MgJR7XpEPo8c6Py+p9q81AMAT4mZRnZETKCJtLYuuS1OttiRwJIXn8v8r9bUJK3lySttXbmh7q8JBsPLjRIhI0b72w9AUThv26hgBjODDTyRNyLmFZy4+XfJKTkzOXLl0a7GqUeK5/8/tiO7Js6+GtWRGJGglWTQT1K9Y3E4GZYIIijYSUNsLCwpapTyRn38/IgxCHxIqJclur20x7ju1x94+MXzNexq4ea0OEXUvudqjRgXNJSKn+Y47mQYgXapSrITc0u8F06MQh91wSDP+dvG6yzSVB0kYYCfpSoiOieR9JqYLmQUgeVCpTSfo36m86lnHMPZfk882fy/SN06V8VHlLIw8j6ZbQTcpFleM9JSUemgchBQDGcEn9S0wnT5+UlJ0pFpF89dtXMufXOTZEGBmCXenkYTyElERoHoQUEjRVIeKAhncZLsv3LHcPAUYzV0RYhC21CyPBXJLq5arzXpMSA82DED/NJYFRQEPPGypr9691G8mIxSPk2cXP2vBgVzr5xNhE3ndSrKF5EOL/oY3Sqlor0/3t75fUQ6lmIugneXHpi6ZmVZq5kzc2imvEIcCk2EHzICTARtIwrqFpcNvBsi1tm3sI8OgVo2XUilFSr2I9t5HAcDCRkZBQh+ZBSBFSJ7aO3NryVtO+4/vcRjJp7SSbT4Ihwja7XY2kQ3wHaw4jJBThL5OQIIFkjNc1vc6EuSSLti0yI5mxcYa8t/49iSsTJz0Te1o/CZI4YiQXIaECzYOQEABDei9veLkJc0m+2/Gd9ZHM3TJXZvwyQ8pFlsuaS1Ivay4J5pYQEkwig/nhpOSRlp4hh9NPybItB6RjvcrBrk6xnUuCaAPKOJ0hi3cttogETVyfbf5MosOj5YLaF5iR9KjTQ+Ji4oJdZVIKoXkQvwHDWL8rTc5kigwcmyKT7+xCA/ERLKF7YcKFpic7Pykr9q4wI0Em4AXbsuaSJMcnm5H0Suwl8eXj/fFPSUie0DyI30hJ3W/GATJOnbHXjD78R0R4hHSM72h6JPkRWff7OhsCDDP55+J/mtpUb+PucK9bsa4fP52Qs6F5EL/RJamqDjMVM5CoyHB7TQI3BLhF1Ram+zvcL6kHU93p5P+z7D+mxpUbm4nATJpUbsK5JMSv0DyI30CU0axmrPV5vHJDe0YdRUhSXJLpr23+KjuO7HBHJG+sfENGrxxtM9pd6eQRnXAuCfEVmgfxK7ExUSY2VwWP2hVqnzWXBHm2EJG8u+5deWftO1K9bPWsSYnaIY8msKjwqCDWlhRXaB6ElPC5JNc2udZ0+ORhm0uCUVuzN82WqRumSsXoipb9F81byAYcExkT7CqTYgLNg5BSAozisqTLTMdPHZfvtn+XlU5+61dmJmUjy9ocEjRtYU5JhegKwa4yKcnmoR13Q7S4T3VK9UlmZuajXo7ZrEWa6jSOc62Hq/uraDFVVV+FY67T9w7o/k66PcZ1uuop3T/DOaejFu+oyqo+VT2QWVoWYifET8AoMLwXwlySH3b9kDUEWPXFli+sKQsrJKJpC5EJVk4kxG/moQ/ynlr0V7XR5/cJfV0jl8N76jH7su17TDVP94/Uc7ENDVWtUSXr/lO6v5Zur9TyI7zW7dGqu1Qpjnn0Vc3x5XsQUtrnklyQcIHpic5PyKp9q9wd7v/47h/WuY4122EkiEqwljshvkYe96hGwjjwQss9BTwfxtPD2Z6gWqAaqtc55nEMGmEtsnCMpKK+/73zeqIWV6poHoT4aS5J+xrtTQ8nPywbDmywFCnoJxm5ZKSpVdVWWVGLGkmDSg1430spvppHE1U3fYg/q2W66mF9sP/g5Tg8/L/Q41C+qce4mqTidXunHaClZ+Si2521GKeqp7rFiUISdHubx3WxjX2EkADMJcG6I9CQ9kPk10O/ume3v/LjK6aGlRqakaDDHcfhHFI6yNM89MfwpRbe4tRhzvlIYNRFdZ5qmh6f5KUPoqvu2uGYw1wt1+vrRbl9rr6/WIuWemxzLSdoiejC2y8zx/4OPQfNW5DUrcvZtoT4AqKMO1vfadp1dJe7j2Ts6rEyZtUYSaiQkP03VvAAAA6KSURBVDW7XZu3sGoi55KUcvPQh/jFuTyc0Ww13TGLJfr6jJbVVHuzXWOHU+7RY9DxjQ5xmMduNEU5UQeapM5p9tL31ul7R3WzlRNp1PF4G9s7cqk7IhyLcpKTk9mpToifQL/HwOYDTb+n/541l0T7SZBKfuJPE6VqTFX3AldYmhf9KqRk4Wuz1UxVL9UCfcCjCStadVanuO5H7uhwfZCnOduXqJ523p6tGqQa6ZSznHPQkLrVaapCs1VT1WZ0uOtrXAeRDiKTW1Wv+fgdCCE+gJFYVze+2nTk5BH5evvXZiQfp34sH/z8gcRGx1r2XzRvIRswRnqR4o+v5oE+iXH6MMfoqJOqQYhC9HVt3R6rm/20RJrPGU5bKD5viu7/zDl/pNPUdYeWv6kGOPsvVD2m+zO0RDRzr8dIrXs8huqiKYud5YSECJgb8qcGfzKln0qX73d8b01byAD8UepHEhMRYxmCYSSYS4K5J6QUmoc+0GEYN3vZj6akfs52qhZtczh/vxa9veyfpMWkHM5Z6jRhEUJCGMxW71m3pynjTIYs273MOtsxcgvpUrDELlZIRD8JVkzEbHhSfOAMc0JIwHFNOoQe7/y4rN632owEw4Cf/v5peeb7Z2x4sGsuCfJzkdCG5kEIKVIwCgujsaCHOj4kPx/42Z1O/vkfnjch1bwrnTyyBZPQg+ZBCAka6AttWqWp6d5298pvh39zG8mry181YYiwGYn2k7So0oJzSUIEmgchJGTA6oe3t7rdtPvobpm/db6Zybg14+St1W9JrfK13OuSoJkLM+JJcKB5EEJCEqzHfmOzG00H0w/aiC30k0zbMM3WJsEQYXS0o5+kc83OnEtSxNA8CCEhT1xMnFzZ6ErT0YyjNpdk/pb5MufXOfLhxg+lQlQFG/oLI+lau6uUiyoX7CqXeGgehJBiRfmo8tK3fl/TidMnZPHOxda0hSHAn/76qZSJKGMGAiOBoVQqUynYVS6R0DwIIcUWGAUMAhreZbgs37PcZrejwx39JZFhkZYeBUaCdCmcS+I/aB6EkBIBJh3CKKChnYbK2n1rzUQQlTyT8oyMSBkh7Wq0c3e414n1TJNHCny/C3oCIYQUh7kkrau3Nj3Y4UHZdHCT20heXPqiCSnkLQuwDgNuGNeQQ4ALCM2DEFLi55I0qtzIdHfbu2Vr2lbrH4GRjFoxSl5f8brUr1jfnQW4VbVWNJJ8QPMghJQqEmMTZVDLQaa9x/bKV1u/sn6SiWsn2nyS+HLx7qatDvEdrDmMnAvvCiGk1FK9XHW5rul1pkMnDsnCbQttLgmG/05ZP0XiysS555IgL1d0BFadIIDmQQghCob0XtHwCtOxjGPy7Y5vLSJB8sYZv8ywIcLdE7pLr3q9pFtCN3tdmqF5EEJINjDJsE+9PqaTp0/Kkl1LzEjQxDVn8xyJDo+2ha2QbwsLXWESY2mD5kEIIbmApiosYAUNP5M1l2Sek7wRKVMiwiIkuWaydbaj071GuRql4n7SPAghJJ8gESOMIln16HmPyk/7fzITQVTy7OJnTW2qt3Gnk0eix5IKzYMQQgo5BLhltZamBzo8IKkHU61/BFHJv5f929SkchN3OvnGcY1L1BBgmgchhPiBpLgkGRw3WAa3HSzbj2zPWm5XI5LRK0fLqJWjbIiwy0haV2ttExmLMzQPQgjxMwkVEuSWFreY9h3fZx3tGAI86adJMn7teKlRtkbWpEQdAtwxvmOxnEtS/GpMCCHFiGplq8mAJgNMh08eloVbF1pUMvOXmfL+hvdtiDDmkqCP5Pza51uyx+IAzYMQQoqIitEV5fKGl5uOnzou327/1vpIEJXATMpFlpNudbpZ8xZKX+eSpKVnyOH0U7JsywHpWK+yn75FFjQPQggJAmUjy1qzFZRxOsPmkpiRqD7f/LlEhUdZJAIj6ZHYQyrHFOzhD8NYvytNzmSKDBybIpPv7OJXA6F5EEJIkImKiJKuCV1NwzoPk5V7V2ZlAdaIZNG2Rda5jr4RV86tmuVr5nnNlNT9Zhwg49QZe03zIISQEjyXpEN8B9MjyY/I+t/Xu41k5JKRJozWchlJ/Ur1vV6nS1JVNR0xA4mKDLfX/oSRByGEhChhYWHSvGpz05D2Q+TXQ7+6+0he/vFlU6O4Rtb0BSNpWrmpey4JooxmNWOtz+OVG9qzz4MQQkorDSo1kDtb32naeWSnLbWLuSRjVo2RN1a+YUOE0UcCM8FM99iYKJO/O8sBIw9CCCmG1KpQSwY2H2jaf3y/LNi6wKKSyesny4SfJtgQ4dORLSX2THvJOJNsHfD+hOZBCCHFnKplq8o1Ta4xpZ1Mk6+3fZ2Vc+vYQjkQ8bUcy7jW5pP4E5oHIYSUIGKjY6VfUj/TgDcXSnrYFr8bByjeyVUIIYTkSLhES7nMxjm+7ws0D0IIITQPQgghgYeRByGEEJoHIYSQwMPIgxBCCM2DEEJIMYg8wsLChqg2qNaqns/hmM2q1aoVqqUe+6uo5qo2OqXNodeyk3MstFJ1lcc5C5zPc71fw9fvQAghpAgnCeqDu6cW/VVtMjMzT+TxIO+px+zLtu8x1TzdP1LPxTY0VLVGlaz7T+n+WroNA/kIr53zBuq224RI6DB18PnBrgIhpBhEHveoRsI48ELLPQU8H8YzwdlGeaVznWMeRhGDXT7WkxBCSAiZRxNVN40KFqsWqs7L4Tg8/L/Q95ep7vLYH68msdMOyCrdkYse1xlNYbq5WnW3h5mA8U6T1fAwV/5hL+Cz0EwG7d27t7DfkRBCSEGbrfTB+6UW3patGuacj36KLioYxzQ9Pkkf9Nkjha66a4fTrIW+jfX6elFun6vvL9aipR7bXMsJWs7RfelOk9V2fR2r2x+qblFNzOEaY7SAJDk5mdELIYQUVeShD+CLVa28aJa+vU01HWahLNHtM6pqXq6xw6NZa4aqk/PWbqdPAyaF8pxmLz1nnRZHVa2c19udMk2LKR7XIoQQUkyarWaqejkPfzRhRavO6hTX/eWdKMG2tbjE6RAHs1WDnG2Us5zjGqgsKtKynhZNVRixFakyc9ISyekv87gWIYSQIsLXlOzjIH2Q4wF+UjUIIYi+rq3bY3Wzn5bxqhlO1wQ+b4ru/8w5f6TT1HWHlr+pBjj7L1Q9pvsznGjmXozUcsznc8c4IlRoUnvLx+9ACCGkKM1DH+gwjJtzaKbq52ynatE2h/P3a9Hby/5JWkzysh/NVx19qTMhhBDfYXoSQgghNA9CCCGBh5EHIYQQmgchhJDAw8iDEEIIzYMQQkjgYeRBCCGE5kEIISTwMPIghBBC8yCEEBJ4GHkQQgiheRBCCAk8jDwIIYTQPAghhAQeRh6EEEJoHoQQQgIPIw9CCCE0D0IIIYGHkQchhBCaByGEkMDDyIMQQgjNgxBCSOBh5EEIIYTmQQghJPAw8iCEEELzIIQQEngYeRBCCKF5EEIICTyRgf8IQgghwWDq4PMDdm02WxFCCKF5EEIICTyMPAghhNA8CCGEBB5GHoQQQmgehBBCAg8jD0IIITQPQgghgYeRByGEkAITlpmZWeCTiiNhYWF7tdhSyNOrqfb5sTr+gvXi/eLvK3SoVkKfE/XUJ6qXWvPw0XiW6n1KDnY9ssN68X7x9xU6hJWy5wSbrQghhNA8CCGEBB5GHvljTED/FQoP68X7xd9X6DAm2BUoynqxz4MQQkiBYeRBCCGE5kEIISTwlLrIQ4et9VVtUP2ieszL++BV5/1Vqg55navbVVRzVRudsnKI1Osp1XbVCkf9irhe41R7VGuynRPs+5VTvYJ2v7RMVH2lWqdaq3ogFO5XHvUK5v2KUS1RrXTq9f9D5H7F5FIvn++XL3XzeD9CtVz1sU/3DPM8SouUCNUmVZIqWrVS1SLbMfgHnYP7qeqiWpzXucrzqsecbfxjPhci9XpK9XAw7pfzXncVfrhrsp0TtPuVR72Cdr+UWqiTsx2r+jlEfl+51SuY9wuvKzjbUdiP90PgfoXlUi+f7pc/fvvO+/9PNUX1sS//J0tb5NFJ9Yt+8VTVSd1+X9U/2zF4PVHfBym6HacuXCuPc1FOcLZRXhki9fIVX+qFH9giLX73ct1g3q/c6uUrha6Xbu9U/ejUL02LdaqEYN+vPOrlK77UCxxxjsFDGsoMgfuVW738gU+/fS3raPFn1Vgv5xTonpU288CPfqvH623Ovvwck9u58fqPtBMbTlkjROoF7nNC13GFCN99qVduBPN+5UXQ75d+bn0t2jt/tYbM/fJSr6DeL6f5ZYVu7lHN1XsTEvcrLOd6+Xq/fK6b8rLqUdWZbOcU+J6VNvNAGJed7H8V5HRMfs4NtXqNVjVUtVPhB/FSEdYrkASqXkG/X/pAqaDFh6oH9T/x4QJ+flHXK6j3S+txWoXPxl/TnbSOrQr4+UVdr9E+3i+f6qb1uEzLPVq3ZYX4XCnt5gEHTvR4jX/cHfk8Jrdzd3uEhbWcvziCXi/9kex2fsj4K+MtJ+QtqnrlRjDvV44E+37pvYhyHtCTtQ7TQ+V+5VSvYN8vj3oc1GKBqm8o/b4ys9XLD/fL17p1VV2h92Sz09zVS7ffLfQ90y9SaqREqlJVDTw6m1pmO+bP2TqbluR1rvJCts6m50OkXrU8zn8IP5iiqpfH+/W9dEwH7X7lUa+g3S/n9UTVy16uG8zfV271Cub9QpbXOGe7rOpr1WUhcL+q51Ivn+6Xv377zjE9snWYF/ieFajiJUHOSISfnRELw5x9d0PONm746877q1XJuZ3r7K+qmqfa6JRVQqRek5xjV6lme/54i6he7znheYbz19AdIXK/cqpX0O6XcqEq0/lstJdD/YJ9v/KoVzDvVxvVcuezMeT676Hw/zGPevl8v3z97ediHgW+Z0xPQgghpMCUtj4PQgghfoDmQQghhOZBCCEk8DDyIIQQQvMghBASeBh5EEIIKTA0D0IIIQXm/wC9VTAC1oZh6wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.errorbar(dtaus**2, ener[:, 0], ener[:, 1], fmt='.', label=\"Monte carlo results\")\n", + "plt.plot(0, E, 'x', label='ED result')\n", + "\n", + "def func(x, y0, a):\n", + " return y0 + a*x**2\n", + "popt, pcov = curve_fit(func, dtaus, ener[:, 0], sigma=ener[:, 1], absolute_sigma=True)\n", + "perr = np.sqrt(np.diag(pcov))\n", + "xs = np.linspace(0., dtaus.max())\n", + "p = plt.plot(xs**2, func(xs, *popt))\n", + "\n", + "plt.errorbar(0, popt[0], perr[0], label=\"Extrapolated value\", color=p[0].get_color())\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/default_variables.py b/default_variables.py index 7590ede5eb1129232eef1f44fb01345e9abee400..6fe9b4d3f299efe275e81c021afd6863572786d8 100644 --- a/default_variables.py +++ b/default_variables.py @@ -73,6 +73,8 @@ PARAMS_GENERIC["VAR_QMC"] = { "Max_Force" : [1.5, "Max Force for Langevin" ], "HMC" : [False, "HMC update"], "Leapfrog_steps" : [0, "Number of leapfrog steps"], + "do_ED" : [False, "Calculate eigenenergies through exact diagonalization"], + "ED_N_Part" : [-1, "Fixed fermion number, for which to do ED. -1 means all possible values."], } PARAMS_GENERIC["VAR_errors"] = {