{ "cells": [ { "cell_type": "markdown", "id": "6e2952a9-913c-4eef-8dfa-0507309eb90d", "metadata": {}, "source": [ "# Hubbard Model by UCCSD" ] }, { "cell_type": "markdown", "id": "1344f070-630d-49ab-999d-ff6dc1c31f02", "metadata": {}, "source": [ "## Overview\n", "\n", "In this notebook, we will demonstrate how to use the `from_integral` method of the `UCC` class to solve half-filled Hubbard model\n", "$$\n", "\\hat H = -t \\sum_{j, \\sigma} (\\hat c^\\dagger_{j+1, \\sigma} \\hat c_{j, \\sigma} + \\hat c^\\dagger_{j, \\sigma} \\hat c_{j+1, \\sigma}) + U \\sum_j \\hat n_{j↑} \\hat n_{j↓}\n", "$$\n", "\n", "using UCCSD. The results show that UCCSD is better than CCSD in capturing strong correlation, yet the accuracy is still not satisfactory.\n", "\n", "We also highlight that the `from_integral` method is a flexible interface for cutomized Hamiltonian." ] }, { "cell_type": "markdown", "id": "d99ad0f6-1280-4853-93cb-5f4f2b472630", "metadata": {}, "source": [ "## Setup" ] }, { "cell_type": "code", "execution_count": 1, "id": "1119cfc3-0e58-409c-852b-fadb87246e3a", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5])" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "\n", "from tencirchem import UCCSD\n", "\n", "# number of sites\n", "n = 6\n", "# half filled\n", "n_elec = n\n", "\n", "# scan U/t from 0.5 to 8.5\n", "t = 1\n", "U_array = np.linspace(0.5, 8.5, 9)\n", "U_array / t" ] }, { "cell_type": "markdown", "id": "20a04916-a709-49aa-bb88-d3f4463011e9", "metadata": {}, "source": [ "## Calculate" ] }, { "cell_type": "code", "execution_count": 2, "id": "f87e468c-a22e-4bbc-b7ed-b742eb246f3a", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Calculating U = 0.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 1.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 2.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 3.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 4.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 5.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 6.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 7.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n", "Calculating U = 8.5\n", "CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH\n" ] } ], "source": [ "# stores results\n", "e_hf_list = []\n", "e_ccsd_list = []\n", "e_uccsd_list = []\n", "e_fci_list = []\n", "\n", "for U in U_array:\n", " print(f\"Calculating U = {U}\")\n", "\n", " # set the integrals\n", " int1e = np.zeros((n, n))\n", " for i in range(n - 1):\n", " int1e[i, i + 1] = int1e[i + 1, i] = -t\n", " int1e[n - 1, 0] = int1e[0, n - 1] = -t\n", " int2e = np.zeros((n, n, n, n))\n", " for i in range(n):\n", " int2e[i, i, i, i] = U\n", "\n", " # do the calculation\n", " uccsd = UCCSD.from_integral(int1e, int2e, n_elec)\n", " uccsd.kernel()\n", " print(uccsd.opt_res.message)\n", "\n", " # record result\n", " e_hf_list.append(uccsd.e_hf)\n", " e_ccsd_list.append(uccsd.e_ccsd)\n", " e_uccsd_list.append(uccsd.e_uccsd)\n", " e_fci_list.append(uccsd.e_fci)" ] }, { "cell_type": "markdown", "id": "c38ebbf8-407d-41c9-bedb-a8f0ce454f45", "metadata": {}, "source": [ "## Plot" ] }, { "cell_type": "code", "execution_count": 3, "id": "20fc2c4a-93a1-4fd5-95d2-9d7e2fe0e4a8", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEeCAYAAACDq8KMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOlklEQVR4nO3de1xb9f348VfCnd4C9E6vQe3NWstFZ1utTujmrX5VqKvOubmV6Drnfttsxve7uZvf1bC5OevUoH7VbVoLUdduzm1Ep9bWC5B6ma2tkl5pKS1wKC1QIMnvj0MCIQnXQBJ4Px8PHpKczzl8iPTzPudzeX80LpfLhRBCiFFHG+oKCCGECA0JAEIIMUpJABBCiFFKAoAQQoxS0aGuQDhatGgRaWlp/T6vqqqK1NTUQf/8YFwnnOoSrOtIXYb2OlKXob1OqOtSWVnJJ5984v2mK8KUlpa6cnNzXWaz2VVaWurasGGDq6SkpNfzTCaTy2w2u8xms8tkMvVY9rrrrhtQ3QZ63lBcJ5zqEqzrSF2G9jpSl6G9Tqjr4u+8iOsCUhQFm82GwWDAYDCQlpZGbm5uj+cUFhai0+nIz88nPz+f9PR0DAZD0Ou2du3asLlOsOoSLOH0O4XTZxNOv5N8LkN/nWAIal0GG42GW0lJiau+vr5f5+h0Op9zevrVgxWpRwL5LPyTzyUw+Wz8C/XnMiKeAPrLbrejKAo6nc7nmNVqHf4KRZhwuvMJJ/K5BCafjX/h+LlE5CBwcXExycnJAJSVlWEymQKWtdvtft/X6XQoijIU1RtRwvGPNhzI5xKYfDb+DfZzcTid7Nx7gmqlmam6BJbNm0SUdnD38BEXAPR6Penp6ej1egDq6urIy8ujpKSkX9dJTk6mrq7O77GqqipWr17teb127Vr5oxZChMzWssMYn6ugqq7J815qciKmWzO4Pmum33M2b97M5s2bPa+rqqp8ymhcrsjOBaQoCklJSdTX1wfs5snJyaH7r5mUlITJZCI/P9/nnNWrV7Nt27ahqrIQQvTZ1rLD3LZpO90bak3Hf/9096UBg0BX/tq1iBsDsFgsXq/djX6grh73k0J3iqIEPCaEEOHA4XRifK7Cp/EHPO8Zn6vA4XQO6PoRFQAURSEvL8+rsXf34wdqzPV6PTqdzm+AyM7OHpJ6CiFEMOzce8Kr26c7F1BV18TOvScGdP2ICgA6nQ6z2ezV2BcVFZGbm+v1JFBYWOh1XkFBgdeMH4vF4rfrRwghwkm10hzUct1F3CDwmjVrvBr42tparwFgq9WK2Wxmw4YNnvc2bNhAYWEhRUVFgLok2mw2D1+lhRBiAKbqEoJarruIHwQeCjIILIQItebWdn74p3L++Kb/8U1QB4KnJyfyyW9X9zol1F+7FnFPAEIIMdJ9WtXA7X94m91HGjzvacBrMNg9C8h0a8aA1wNIABBCiDCyZed+vvt/79PU6mDyhHieNCzjVHObzzqA6b2sA+gLCQB+uBeCjbYFYHa7HbPZTGFhIXq93pMwr7a2FkVRMBqNXgPw7vJFRUUkJyeTm5tLQUEBOp0Oo9GI1WrFZrORn5+PwWDwDOK7y3dNyFdZWUlxcTH5+fk9ruwWYqRzuqCp1cHlC6fw5J3LmNLRv39tRuqAVgK7F4T5WwgWccnghkOokja1Oxyut3ZXu4p37ne9tbva1e5whKQe2dnZrvz8fK/3KisrXTqdzlVRUeFTPj093ae8y6Wm4NbpdH0uX1FR4fd9IUa6s23tXq//VnE46P/+/bVr8gQQJgay1Hs46fV6MjMzWbduHRUVFV7H3HmZuvO3Mrun8unp6QPaiEeISOVyuXj6jUp+/8purPetYtL4eACuSZ8xLD8/otYBjFTupd7dF3wcrWvitk3b2Vp2OEQ18zYcCfQCBQ0hRppTzW3c8dhO7nn6few1p3nq9c+GvQ7yBBBkZ862BzwWpdEQHxvlVdbhdHLvn8sDLvXWoC71vuL8KQH7+7QaSIgd2v+ViqJgtVp54okngn5tq9WKXq9Hr9fLAj0xKnxwoI6v/+FtKo+fJkqr4Wd5S/juVQuGvR4SAIJs6rrigMdWLZnOiz+43PNav/5FmlodPV7PvdR7wfe2cqq5zW+Z9LnJvPnzLw+kugHZ7XZP3iW73U5lZSWvvfYa6enpfsuXl5f7rMAuLS3t088qKSnBaDQOrsJCRACXy0WRdR//vXkXre1OZqYk8vS3l3PxuZNCUh8JABHCOczr9fR6vWerTZvNxpYtW8jLywtYPjMz02v1NajdOeXl5X7LuwNGbW0txcXFEgDEqPDoP/fyo+dtAFy9NJXH1n2B5LFxIauPBIAgq35iTcBjURqN12v7H25ix94abvrNG71e99nvrGD5vMl+j2k1ft8OmvT0dAoKCsjLy6O+vj4o1+waMLKysoJyTSHC3Vcv0/P0G59zxxXncNeqeWg0Q/yPtxcyCOyHex1A180U+mpMXHTAr679/+6yV54/ldTkRAL9GWhQZwNdef7UgNcd6v5/UJ8IFEUJmHZ7MLKzswPODBIikrlcLv5acdizH8mExFh23n8V3/7S/GFr/Ddv3szq1av9rgOQAOBHamoq27ZtG5ZFYFFaLaZbMwB8gkAwlnoHi3t2js1mG5Jry+wfMdLUNp5lze/e5Jbfb+eJ1zpn+MRGR/VwVvCtXbuWbdu2kZqa6nNMAkAYuD5rJn+6+1KmJyd6vT89ObHPu/0EU11dnc92me59FcrKygDvjXkCba0ZaMpooPJCjBQ799aw/Cev8o8PjhIXoyUmKjybWhkDCBPXZ80c8FLvYHGndnCnezAajZ7UDqDO1jGZTFgsFs8mO2azGbvd7kkV0TUVhMViQVEUDAaDVyoIu92O3W5Hp9ORk5MjG/OIEcPpdPHbV3Zz/4sf4XC6OGfqOP74nRUsnpUU6qr5Jemg/ZB00EKI/qppaCbf/A6v/acagJuXzeF3t2cxLiEmxDVTSTpoIYQYIvaa07yx+zgJsVE8+LVMvnqpPuSzfHojAUAIIYLgC+dOYtMdF5GpT2HBDF2oq9Mn4TkyIYQQYe5YfRO5D77BniOK573bLkuLmMYf5AlACCH6zfrRUdaZ3+Fk41lqT5/l9ftWhX13jz/yBODHYBaCCSFGDofTyfY9xyl55wDb9xznbFs7Py3+gBt+8wYnG8+yeJaOovxLwrrx72khmMwC8kNmAQkh/O3RERutpbXdCcC3vnguG29J91nhH65kFpAQQvSBe4+O7nfH7sb/7qvm86u1/jPjRhLpAhJCiC4cTifG5yr87tEBaoqWl947hMPpHM5qDQkJAEII0cXOvSd8dufryr1Hx869J4avUkNEAoAQQnRRrTQHtVw4kzGAMFNc/Sb3fPooD89fT97Uy0JaF/cmLSkpKeh0OpKTk8nNzcVoNGIymQZc3r1zmHuP4dzcXCwWCxs2bMBut2MymSgqKkKv12MwGACora0FIC0tTbaNFEOmpdXBy+8d6lPZqbqEIa7N0JMAEEZqztZj2P0QSvsZ8nc/xMqkxUyOG/4kUjabjXXr1mEymbwStdntdvLy8rDZbF4Nen/KuxPDdd1asutOY3q93pMwTq/X++wyZjAYyMvLo6SkJOi/txjdKo83cvsjb/PhwZ43PdKgZupdNi802zgGk3QBhQmXy8Wdex6m0aE+VjY6mrhrz6aQ1CUvL8+nMQe87sgHWr64uNhnX+H+bDRvNptRFIWioqI+nyNEbyzvHuDSn7zKhwfrSRkXx72rF6EhvPfoCIbI/w2GQCgWghUff5OXa3bgcKkzCxwuJy/VvE1x9ZvDVgfo7MYJlKI5Ozvbky56IOX97SrmTgvdV3l5ebKHsAia2sazfO+ZMhpb2lk2bxI7f3kV9+UuCas9OgZDFoL102AWgp1pDzwwFKWJIj4q1qdsTavC0ne/zan2JlxdJp9p0DAhegx7lz/FmKj4gNfVarQkRAVnY+m0tDSys7Mxm80By1itVk+D39/yGRkZKIqC2WzucR+AnJwcT3dQd4qikJSUREVFhc/ThBADsbXsEB8drKfghsVEd9m8xeF0hnSPjmCShWDDYOzr1wc8dvXEi3gl/X7P68lvrKHJeTZgeRcuT1fQW/Ufc7KtwW+5zPHnUfaFRwZe6S7sdjtpaWk9lunez9+f8iUlJeTk5Hju+LOzszEajf3aFMa9QU15ebkEADEgz223M1WXwJWLpwFwfdYsrs+a5VMuSqvl0gVThrt6wyYyQ9ko4u4Kanc5Ql2VoNDr9VRWVlJaWuqZ9ZOTk+O1xaQQQ+V0Sxv55ne484l3+dbjOzlxqiXUVQqpiHwCcE8jrKysBOix+wHULgiz2ezpVigtLSUrK4vc3Nyg1+30F7cGPBal8c4ZUnN5MS6Xi1s/foBXTr6PA9+VhVEaLddPWsYfz7834HW1muDFcXcD3RP3DJ2BlHfLzs4mOzsbk8mE0Whk3bp1ff7/4d5ruPs1hejJfw7Vc/sfdrDv2Cm0Gg13rZpH8tjY3k8cwSLuCcBoNLJhwwY2bNjgafh7G0BUFAWbzeaZgpiWljYkjT/AmOiEgF9d+//dZcfGJPLEov/HuOgENN3mHGjQMD4qkccW3N3jdYPV/w+Qm5uL1WrtsUzX4/0pryiK3zt9k8mEoigBN5Hvrry8HIDMzMw+lRejm8vl4v/+/TlX/Pxf7Dt2imlJCbxScCUbrj8/YvvzgyWifnt3Q961oTAYDFitVp+ZJd1VVFTgcrmorKwMu4VEk+OSeHzhPV4DwKCOATy+8J5hXQvgnq8fqFFXFIXk5OQBly8rK/NbTq/Xe/r2e2M2mzGZTH0uL0avtnYndzy2k3uefp+WNgerLpjGjl9exYr5k0NdtbAQUQEA1Lu/ro29uxugr3eP4WrNlJXcMHk5UR3dOVEaLTdOXsGaqSuHvS4lJSUYjUafRt09/77701N/yhcVFfmUs1qtfX4iKywsRFEUnwViQvgTE60lLiaKKK2GX958ISXfv5xJ4wPPqBttIn4aqMViIS8vj/r6+oB3hBaLhbq6Os+daFlZmd9UBm6h2g+g5mw983bcgdJ+hqTosXy6/KmQrAR2657aAejx6am38oqiUFxcjF6vx2azeZ3rbtDtdjtms5nCwkKfVBCKopCWliaNv+iRy+WiudVBYpw6xHnmbDt7jihkpk0Mcc1Cy1+7FvEBICMjA4PB0GPDZLPZ0Ol0nqeFoqIiSktLA6YTyMjIIDU11fN67dq1rF27NrgVDyCccgEJEWmUM62sf+o9WlrbKfn+5Wi14btT11DbvHmz12LWqqoqKioqvMpEdAAwGo0DSg7mXkgU6KlBdgQTIvKUV57k63/YwcGTZ4iJ0vLafatYOje59xNHCX/tWsSNAbhZLJY+N/7dZ564G/3eBo6FEOHP5XKx6dU95NxfysGTZ5g7eSzWn+RI498HEbkOwD2I6G78FUWhrq7O77xwRVHIy8ujsrLSZ8BY5pELEdlqG89y5xPv8I8PjgJww0Wz2HTHRUxIHN3z+/sq4p4AbDYbNpuN9PR07HY7NpuNjRs3egZ47Xa7Z6EYqHf7ZrPZq7F3z0yRaYRCRLavbtrOPz44SlyMloe+nsWz65dL498PETUGoCgKc+fO9Tvl0/1rFBUVYTKZvFandk8fXFtbG5azgIQQ/WOz1/Ltp97jCcMlLJ4VuhlzkWBEzgIaChIAhAhPNQ3NlNtruXrpDM97TqdrVM/26asRNQgshBhd3txdzbIfv8ptm97mwwN1nvel8R84CQB+hGJDGCGEmn9/+57jlLxzgO17juNwOnE4nfzqpY+4zvQ6xxta0E8ZR1xMVO8XE4BsCNNv0gUkxPDbWnYY43MVVNU1ed6bqosnaUwse6pOAXDbZXp+c1umZ5Wv6DvZEEYIEZa2lh3mtk3b6X43Wq20UK20EBet5ZFvXsxXls8NSf1GKgkAwsOdh6eoqIjk5GRyc3MpKCjwTJc1GAwUFxeTmZmJyWTy2Y2rey4g9zWMRqPfWVd9Le+e1qvT6VAUhdzcXCwWi2dDGZPJRFFRkU/uIGBAK8XF8HI4nRifq/Bp/LvSjYkl75LZw1an0UICQDhoPAQtJwMfj58I43y3qws2vV6PyWTCarV6GvmuzGYzOp3O532bzca6deswmUw+20Xm5eVhs9m8zulPefceDl2DTV5enledzWazZ9OZ7oniDAYDeXl5AfM+idDbufeEV7ePP8cbWti598SI3p4xFGQQONQaD8Hz86AkI/DX8/PUcsOka/7+7vzt/5uXl+fTmANed+QDLV9cXOzzpPHEE0/0+ju4mc1mn3UgIrxUK81BLSf6TgJAqLWcBEcv+5I6Wnp+QgghdzdOoE3ds7OzvVZh97e8oig+OZt0Ol2vu8B1lZeX5/m5IvykjOvbjnZTdQlDXJPRRwJAsLWdCfzV3uJbtr2PdzXtzT1cN3R3RhaLJWBj7ta18e1v+fT0dHJycnw2kelPv/6aNWs8u8mJ8HLgxGl+UfJhj2U0QGpyIsvmTRqeSo0iMgYQbE+MDXxs1tVw7Sudr5+eDO099316vHIttCr+j03KhDz/Wy0ONbvd7rdbqKvu/fz9KV9SUkJOTo7njj87Oxuj0dhrEOnKPYhdXl7u050kQmdr2WHWP/UuDU1tJMZG0dTqQANeg8HuJV6mWzNG/f69Q0E+UT9kIVj40Ov1VFZWUlpa6pn1k5OT43dzeRE5iqz7+Oqm7TQ0tXHRORMpe+Ba/nz3pUxPTvQqNz05kT/dfSnXZ80MUU0jX08LweQJwI/U1NSBLwRbdzrwMU231YvfqIGTH8DLK3q/7jV/g4kXBrju8MXx7on43A10T9wzdAZS3i07O5vs7GxMJhNGo5F169b1eR9hSf8dfq5Jn8GvXv6Y2y7Tc99NS4iJ1jJr4hiuzUhl594TVCvNTNUlsGzeJLnzHyT3joarV6/2OSafbLDFjAn8FR3vWza6jwNb0Qk9XDe4g2N6vb7Pm+Xk5ub69M931/V4f8oriuL3Tt9kMqEoit+ssP6Ul5cDkJmZ2afyYmh0zd+TmpyIzXQtv7x5KTHRnc1QlFbLpQumkHfJHC5dMEUa/yEmn67wkZOT0+cA4J6vH6hRVxTFa1ppf8uXlfkf29Dr9X3ez8FsNmMymWT/hxBpbm3nu0+/z4r7/sHfKo543k8e27fZP2LoSAAItfiJEBXfc5moeLXcMMnNzUWv1/vMnS8sLPQ7+6akpASj0ejTqLvn33fvqulP+aKiIp9yVqu1z90/hYWFKIris0BMDI9Pqxq4/Gf/5Ol/f45GA/uOnQp1lUQXkgzOj2FPBhcmK4G7KywspLa2lpSUFKAzMATSPbUD9Dxds7fyiqJQXFyMXq/3mcLpbtDd6SsKCwt9UkEoikJaWpo0/iHy3HY733+2jKZWB5MnxPOkYRlXnD811NUatWRDmD6SbKBCDNzplja+/2w5m3fsB+DyhVN48s5lTJGFXCEl2UCFEENu+54aNu/Yj1aj4X9uXMwPrlsog7lhSgKAECKorlqaSsF/nc9lC6ewYr4kbwtnEpb9kIVgQvTdqeY2vvfM+xyr71zV/t83XiCNf5iQHcH6ScYAhOibDw7U8fU/vE3l8dNcef5U/rLhi6GukghAxgCEEEHhcrkosu7jvzfvorXdycyURApuWBzqaol+kgAghOgX5Uwr6596j23lhwE1rcOj37pYFnZFIAkAQog+23fsFDf++t8cPHmGmCgt93/lQu5aNQ+NRtP7ySLsSAAQQvTZ9KQE4mKimDt5LM98eznp+pRQV0kMggQAIUSPlDOtjE+IQavVMDY+huLvr2TiuDgmJMaGumpikGQaqBAioJ17a/jC//ydTf/41PNe2pRx0viPEPIEIDzceXWKiopITk722qC9srKS4uJi8vPzPRk93brn9ElOTiY3Nxej0YjJZMJms7FlyxbPdXNzczEYDJKfP0w4nE6fHPwaNPz2ld3c/+JHOJwunttu565V5xEbHdX7BUXEkHUAfmRkZJCamurZSGG0ycjIIDMzE7PZ7PW+zWbDbDZ73rfZbKxbtw6TyeSz7aPRaMRms3lt/hLouiJ0tpYdxvhcBVV1nYu4puoSmDQujo8PKwDcvGwOv7s9i3EJMSGqpRiMzZs3s3nzZqqqqqioqPA6Jk8AfgxqR7ARoGs+/q7S09O99vPNy8vDbDb77M/rzsrZ9Qmip+uK0NhadpjbNm2n+x1gtdJMtdJMbLSWh76exVcv1cssnwgmO4KJoHGnbXZ3+wTanD07O1u6eMKYw+nE+FyFT+PfVdKYWG5ZMVca/xEsIp8ACgsLPQ1RXzf7GMg5A9J0JvCxqCiIi+9bWa0W4hMGVjbIrFYrer0evV7vyddvsVgCNv5u7iAhws/OvSe8un38Od7Qws69J7h0geT0GakiLgC4G3J3Q2S1WjEYDD32Kw/knAG7eGzgY5deDY++0vn68snQHOAfYeZKePqNztdfngP1ATaNWZQJL/jfOjEY3Dt4dWW32726g/zpLUCI0KlWmoNaTkSmiOsC2rhxI2vWrPG8zs7O9tm6MBjnjHbl5eUUFhZiNBopLi4OdXVEkE0a37e0DVNlE5cRLaKeAOx2O4qi+N3c22q1+r3jHMg5g/Le6cDHorpNoXujJnDZ7hto/ONA38sGQWZmpqebLCsry+e4Xq/3muHjj91ul3GAMHSk9gz3v/hxj2U0wPTkRJbNmzQ8lRIhEXEBwB+dToeiKEE7Z1ASx4S+bJD5C5K5ublYLJYez7NarT3uCSyG3993HeHOonepP9NKQmwUza0ONOA1GOwe8jXdmiE7eY1wERUAAklOTqauri5o57g3hHEbresB3Pw9PZlMJiwWS8CnKEVRZNpnGGlrd3Jf8Qc80rGid+mcZJ5Zv5yPDyk+6wCmJydiujWD67Nmhqq6Igjc8//d/G0IMyICQH8b/97OGe3rAOrq6vrUdVNSUuJ3IZiiKBQVFfnMtOrrdUXwabXwScfCrvVfmsfP11xIXEwU+injuDYj1WclsNz5R77uN67+1gFEVAAI1HgoihLw2EDOGa3cqSDsdjt2ux2dTkdOTk7AcZL09HQqKiowGo2UlpZ6UkEAXo2/OxWEezzGaDRSUFDg98lCBJfT6UKr1RCl1fLknZdQbq/l6qUzvMpEabUy1XOUirhUEElJSVRUVHg13hqNhp5+jf6eI1tCikjX0urgxy/s4my7g013XBzq6ogw4K9di7jnvIKCAqxWq+e1xWLxGmi02+0UFhb26xwhRpLPq0+R/ct/Ybbu45k3Kvn4UH2oqyTCVMQ9AYD3qt7Kykqv7JRFRUWYTCafKYo9ndOdPAGISFXyzgG++/T7nG5pJ2VcHEX5l7BqyfRQV0uEAX/tWkQGgKEmAUBEmubWdjb8uYJn3lBvfJbNm8TTdy1nenJiiGsmwoW/dq1fg8AHDhxgzpw5wayTEGKQXC4XNz34Btv31KDRwIbV5/Oj/zqf6KiI6+EVw6xffyH+uk0aGhr49a9/zQcffBCsOgkh+kGj0XD3lxcwZUI8W+/9Ij++6QJp/EWf9NgFdM4555CWlkZ6ejo5OTmUlJTw2GOP+S372muvodFo+OIXvzhklR0uo31DGBH+Tre0se/oKa9N2c+cbWdMXETN7BbDoKcNYXoMALt27WLLli1YrVZsNhsajYb09HSys7NZtWoVV1xxhVf5l156iRtvvHFofothJGMAIpx9cljha4+8zYlTLbz9y6uYNTF0aUJE5Oj3GMDSpUtZunSp5/WqVas8i39MJpNXQEhLS6O0tHREBAAhwpHL5eKZNyrZ8OcKWtocTEtKoKahWQKAGLB+PS/q9XoeeOABz2ubzcZrr71GaWkpdrudH/3oR0GvoBACTjW3cc/T72N59yAAORdMw5x/CZPGx/dyphCB9SsA5OXleb1OT08nPT2de++9N6iVEkJ0+vBAHbf/4W0qj58mSqvhp3lLuOeqBWi1slWjGJx+BYArr7xyqOohhAjgT2/ZqTx+mpkpiTz97eVcfK7k6BfB4TcA/OhHP2LVqlUjYkaPEJHul1+5kJhoLfeuXkTy2L7t5CVEX/idLGwwGPjXv/5FZmYmd911l8zxF2IYlVee5M4n3sXhdAKQEBvNxlvSpfEXQef3CWDu3Lmewd5du3bx+OOPU1FRQXZ2NgaDQVYDCzFIDqfTJwe/VqPhkX98yn3FH9DucLF4po71X54f6qqKEaxfuYBee+01zGYzDQ0N5OXlsWbNGsaPHz+U9QsJWQgmhtLWssM+u3BN0yUwRRfPBwfUzJ3/lTWTR755MRMSY0NVTTFCDHghWE9efPFFzGYzGo0Gg8Ewoub/y0IwMVS2lh3mtk3bCfSPLjpKw29uy+SOK85Bo5FZPiJ4Bp0MrqubbrqJm266iYaGBoqLi1m1ahVJSUkYDAYZPBbCD4fTifG5ioCNP0Dy2Di+fnmaNP6iU+MhaDkZ+Hj8RBg3a0CXHnTikAkTJrBu3TrWrVvH/v37sVgsbNiwgaysLAwGAxdeeOFgf4QQI8LOvSe8un38qWloYefeE7JFo1A1HoLn54GjJXCZqHi4Ze+AgkBQUwbOnTuXe++9l/LycvLz83nhhRe4+eabg/kjhIhY1UpzUMuJUaDlZM+NP6jHe3pC6MGQpQ7snkdIiNFu0vi+TeOcqksY4poIoQoYAFJSUkhOTvYkesvOzpbuHCEG6Fh9E6at/+mxjAaYnpzIsnmy0ld0OGLtvcwgBAwALpeL4uLigHfxDQ0NFBUVkZaWNqJmAAkRbK99fIxvPb6Tk41niYvWcrbdiQa8BoPdQ76mWzOI0spmLqNG2xmo2w11H0Pdf6DuE7h6G0R1PC0e2zGkPz5gAMjMzOyxC2fChAnce++97Nq1izVr1pCWlkZBQcGIXBcgxEC0O5zc/9JHPPjX3QAsnqXj2fUr2H2kwWcdwPTkREy3ZnB91sxQVVd0NYQzb6i0wL4/Q+1/4JQdus8LU/ZBymL1++mXw4Ghm5IeMADo9fo+XWDp0qUUFxeTl5dHcnIy7e3tQatcqFRVVbF69WpZCCYG5aX3Dnka/2998Vw23pJOfGwU504bz7UZqT4rgeXOP0wMZuaNywWNBzvu5v8DtR139l9+ESaco5ZRPoP9WzvPSZisNvjJ56tfiVM7j6WuHPSv03UhWHcBA0D3ecgvvvgipaWlJCUlkZOT4zPX/4knnuDFF18cdGXDQWpqqiwEE4OWd8lsSj86ylVLU7nx4tlex6K0WpnqGa76OvOm+URnANi/DWwb1S6ctkbf8rUfdwaA2VdDzBi10U9aBImTg1v/btw3sqtXr/Y51uMYQFfuhV/nnHMOiqJgt9vJzs725AXS6XRkZ2cHt+ZCRJDWdgcPv/ophuzzGJcQg0aj4Yk7l4W6WmKoVO+EyRnq946zcPxd9XttDOjmd97VpyyGqZd0njdxifrVF/ET1aeN3p5G4icO6Ffo8xOAW3Z2NiaTyW9ff1+7jYQYaQ6cOM03/rCDcnste6sapOEfDRo+7/x++qWQ8wKknA8TzoOomOD8jHGz1K6m4V4JbLVaefnll7nyyiu9GvukpKSAA706nW5AlRAikm0tO8z6p96loakNXWIM12cNcHBQDC9nu9plU1PW+XX+XTApo2/nn7Om8/vEqXDuEC16HTdr4APOvQgYACorK8nNzQXUO/ucnByys7Ox2+0BLyb5S8Ro0tLq4Mcv7MJs3QdAVloKz6xfIZu0h7OmGqj4X7WxP7nLt2ulemffA0B05C/YCxgA3F09paWlWK1WHn/8cR5//HFA3Qw+OzubVatWeT0hKIoyLJUWItQOnDjNVx/ezocH1fTN37tmAffdtISYaJnJMyDBnHbpnolzolxt6MfrYZFBPRYVBx8/3Fk2djxMyoTJWerXlIuhuWbgv0cwtLVBTRVMmQHRHU30X/8Mr25Wf7fH/h60HxUwAOTl5XnSOWzYsAFQN4exWq2UlpayZcsWTzpo9wCwzWYLWsWECGcJMVEcU5pJGRdHUf4lrFoyPdRVilyDTXjmcsLBVzq6cToa/a7BZNqlnQEgbgJc9AsYN1dt8HXngqZb0B6uALBnF9jehupDcOwQVB9W/3viqNrQ/3UvzDlPLXvEDtv/DjGx4HRCkKYMD3g/AID9+/d7nhCsVisNDQ04HI6gVCyUZEMY4U9bu9PrDv/9z08yIzmR6cmJIazVCHDCBiV96HbJq1Ab7hMV0KpAmtpFjcsFz0z1bri10ZByAUzKUgdoz7u17/UZTEBqPdvZkLu/ujbwm7bB7HPVso/9HB79mf/rx8TCk69B+gr19Z5dsLsCps2Ci6+EqKg+/zpDsiGMP5mZmZSXlwfrciEjG8KI7j6tauDrj+7g3usWcdMXZvd+gui7vgaAsTPg9BH1+zHT4fYuC5u23w2tjTC5ozsnZQlExw+8Tv66pFwuqK+DmmqoPwNKs9qw3/4DmDpDLfPYL+DRnwa+rvlfsCxH/f7Nv8HWZ2DqLLVhnzqz47+zIHlS0O7y3YK6IYw/sg5AjETPbbfz/WfLaGp18MsXP+T6rJlER0lf/7BzN/7j9Woj397S2chfuik4P6OtDar2w+RUmNRxd//y/8GTG9U7+Nazvudcdk1nAJg2C+ITYNrszkbd3cBPmwUL0jvPW3mt+hVCQQ0A7o3khRgJTre08YM/lvP82/sBWLlwCk/duUwa/2BwueBUJRx9C+wv9+2cSzfBuWshPmXwP1+pVbtVDu7r+PpM/W/VfnA44AkrfOFKtaw2Cg51zPnXaGDS9M679WmzYHKX8Z/rboPrb1fLRYAh2w9gqBQWFgLqNFUAs9ncY3mr1YrZbCYnJwe9Xk9paSlZWVmeKa5C+PPJYYWvPfI2+46dQqvR8D83LuYH1y2UfD2Ddeif8OnTasPfdKx/505d1vfG3+WCuhPeDfw1t8K556vHX/8L/PRb/s9NSASlS/fP8i/D0292NvYxsYF/bj/65sNBRAUAo9GIyWTyvDYYDOTk5FBaWhrwHEVRsNlsWCwW9Ho9RqNRGn/Ro6N1TVzx83/S3OpgWlIC/3fXclbMH9p8LSOOsx1O7IJjb0FaXudgqbIPPt+ifq+NhSkXwYRz1aAwEC5X5932JxXwx992NvqnT3mXnaHvDAD6BTB3Psw+Tx2UnX2eOuNm9nkwaZr3HfzEKerXCBQxAcDdkCuK4llxbDAYyMjIwG6395iGoqKiQlYpiz6bnpzIt648l71VDTyefwmTxg9iMHG0aG9Rp18efUtt9Kt3Qttp9VisDhZ+U/1+1pfVaZjTLlMb/+gEdRC4pwDQDjQAb70OSqn3Xf09G+GGb6jlmhrh7893nqfRwPTZaqM+61y10Xe7cBls2xPEDyAyRUwAACgvL8dut5Oerg6kuBt9WYAmBuvDA3Ukj41jZscq3l+suRCtRoNWGxl9uQMymMVXXe+8j70N27LVhGhdxenUOfhjpnW+pzsXMn/i+3Oi4uFsC9QCCYA728xRwJM5+V7fehzc1/n9uYvh/5k67+pnpkGcBO+eREwA0Ol01NfXe71ntarbpfWWhK64uJjk5GQAysrKvLqRxOjhcDp9cvBrNRqeeO0zCp63sXRuMq8WZBMTrR35A739neveUqs29EffUr/mXAtZHdMdkxaBo1XNhzPtMnXe/bTL1MRo3RdZuSm18OkH8Oku9b+7Z8LBSnWR023fgry71HLHj8HWa2HseJgzr6Nx7+iumXVu50IpAF0K3LEhGJ/OqBExAcCfjRs3Yjabe+ze0ev1pKene4JEXV0deXl5lJSUBDzHvSGMmywIi3xbyw777MI1TZdAakoi5ZW1AKSMjaOlzTE60jn0Nef9zh9C/R51U5OuYhI7A0B8EnzVDuNm+85+cbnU6ZPtbeodOcDhSrj6HP8/M3myOsd/Usd0yRQnvFkDSRMjZmZNuHAvAHPztyFMUBeC9YfFYmHLli29lisoKPB0+XRlNBpJS0sjPz+/Xz9XURSSkpKor68PGDhkIdjIsrXsMLdt2t594z2PKK2GX61dyl2r5o2ehIZ9XXzVVdKCjjv8y9SunXHdtq9sb4cDe9XplZ/ugr0fqHf3DXVw1Voo7OifdzrhkgmQMgXmXwjzl3Z8Xeg7ACuCZsgXgvVHbm7ugGfjWCyWPjf+FovF6+e4G/2uYwli5HI4nRifqwjY+IN652/IOW/0NP6gztLpi3PWqikXpq3w3rmquQmqDkDqHPV1ezusSIYzfnbDio6G1i5PG1otvHFcnW4pQiriuoDc/f7uxl9RFOrq6vyOAyiKQl5eHpWVlT4DxrJ5zeiwc+8Jr24ff2pOtbBz74nRsUXjgVfgs+fhwN/6Vn7pDyF6FnzQcUe/p+O/B/bCwgzY/L5aLjoaps9RF1LNvxDmXQgLlqr/PWcRxMZ5X1ca/7AQUQHAZrNhs9nIzc3FbrejKApbtmyhoKAAUO/qLRaLJ3upTqfDbDZ7NfZFRUXk5ubKtNBRolppDmq5iOJyqXvRJi3o3KHq0KtqAPBbHmgGurbN93wTPv7Af/n6E96ZKZ98TR2IlcVyESNiAoCiKFx55ZUoioLRaPQ65p7V41716w4AAGvWrPGsHgaora3tcQBYjCyT+ziHf6ou8jf3AKDtDBx5TU2PfPDvcOYIXP8GpK5Uj5+7FmLGqhuUv74OTqBOtTwKVHdc4w7A3RumU2fPMeuczn76+UvVu/uJU71/dvKkof3dRNBFTADwNw20u/z8fJ9xAZ1O5xUQxOhx4lQLv/3bJz2W0aAu/Fo2L4Ibr+YT8NkLaqN/9A3v+fjRCWrOHXcAmLYcbAfhyUdhF9DW7VpaoAlwb2r2nR9C4cvqNEwx4kRMABCiP97afZw7HtvB8YYWYqO1tLY70YDXYLD7Jtd0a0Zk5fhxtEFrAyRMVF83HoK3v9t5fNwcmH0NzPgynBoP77wLsxphzDj1+L6P4P0d6vdxwPSOr2lAMuBOZxMVDzMXSeM/gkkAECOKw+mkcOsnPPCX/+B0uZifOoE/rl/OvmONPusApicnYro1g+uzZvZwxSAb6OrbpuNq//3BV+Dwv0B/I3yxI33CpKUw9waYdBGcnQOfHoTNb8KuWzpn5Zx3Aaz4svr9qjx1uuXCebimp3Dv50/yZv3HOFxOojRaViZdwG/m5fdvG8YgKq5+k3s+fZSH568nb+plw/7zRxMJAH64F4LJArDI4nK5WPO7t/jXh0cBuO0yPb+5LZPEuGgWzNBxbUaqz0rgYb3z7+/q25pydbbOwVfU/W27qunyWqOFxDvgO7f4TsMcnwQZl0Hi2M73FmWoX0Bx9Rs82LQf4jqPlzfZucjRyJpxwz9NuuZsPYbdD6G0nyF/90OsTFrM5LikYa/HSNJ1R7DuJAD4kZqaKgvBIpBGo+Hqpans+LSGh76exVeWz/U6HqXVhnaqZ19X37acVAPA61+Hui5jGCnpwFI4ngBvVcKYP6n550HNdHmmEcbpIHMlZF2ufp13QcBZOTVn67lz9+/RoMHVpXNMgwbD7t9zedIFw9r4ulwu7tzzMI0OdUZWo6OJu/Zs4sUL7xu2OoxE7hvZrtkN3CQAiIjW7nBSVdfE7EnqHewdV5zDly9MJXUk7NOrz4PGSVA3Cex18OH7cMbWeXysrjMA6BdAyS41IVofctIfbTnJf33wc061N3k1/gAuXJxqP+NpfCubjvJ01T8ZExVPYlQ8Yzq+EqPiGBMVz/wxM5kRrw6itznbaXaeZUxUPFGa/uXGLz7+Ji/X7PC8dricvFTzNsXVb7Jm6sp+XUv0jQQAEbGO1J7hG4/u5Fh9E2//8ip0Y2LRaDSR2/g7Uefhu82/B9b9Qp1r79b1Dv8LXbZg1WjUKZpAs+MsB5qPY28+xv7mavY3V2NvPsbdM6/niylLAdhS/SZlp/b2UBUXL9W8zSenD3CouYb/3b85YNmH5t3FPbNvAOC9hk+5tOz7AMRqYjqCRZwneNwz67+4PXUVAAebj/Or/ZsZE6VO1X3ssO/itFA9jYwWEgBERHp1VxWGoneoP9PK+IQYPjmisHxehG3a4gROAlWo8/CPASnA1zuOj9fBRV+EhDGdXTrnLsahhaqWWvY3VzO3uYZZCerv/XrtLm792ER1a53fH7dCd74nAGQnLyU5ehz17Y1+02REoeX6yctYNHYOLpeL78y8niZnC2ccLZxpb6HJeVb93tHCtLhkz3lNXbq4Wl1ttLa3Ud/eOS5R1+X7qrMnKTry9x4/IhcuGh1NrNv9O3KnXEbG+HOZN2ZGv58uhH8SAEREaW138NPiD3nkH58CsHROMs+sX45+yrgQ18wPx1k19cKYaTD1ks73jwCfAQeA7kMCCmrmzA6HNz3H88f+rd7Fn3mB/e88xMHmGtpcai6f384z8P9m3wTAhJgxnsZ/fHQicxOmMjdhKvqEacxNmMrlSRd4rrt4vJ49y59k3o47aOjWDaRBw/joRB5bcDcA54+by6YF6/v0K+ekZNB85d88wcEdNJocasA4L3GGp2xq3ER+nvY1DjQf5+mj/wz8MbqcbDvxLttOvAvAmKh4LhyXRsb4cz1f88fMlKAwABIARMQ4cOI03/jDDsrtavrm9V+ax8/XXEhcTBj9w3c54dgO2PdnqCyGswrMuh6u/UtnmY+Agx3fx9I5D386kAIv1r7LTdMuBuD4WYUfffaUz4+J1kQxO34KMZrOf8KLxsyh7OJHmJswleSYcb0mt5scl8TjC+/hKx/9yvtXwMXjC+8ZUJeLRqMhPiqW+KhYUuh5/cDshCncl/ZVXC4XSvtptp14B4fL6VMuSqNl2YSFuADbqc8442hhh/IJO5TOAfLHF3wXw8xrATjRqnD8bD3zx8wiWjv4v42RPC1VAoCIGL+0fEi5vRZdYgyPrbuEazNm9H7ScKn/VG309z0HjQfAgXqnfzAR9v8TLtyPK9alLj6bj7rSNg210e82Scfe3LlZelriNG6d9kXPXbz7jj41PsXnjjc+KpbMCefRH2umrGTL5Dc9jW+URsv1k5YN66CrRqPh8QXf5d91H/h/GolKxLLkJ0yOS8LhcrD3zBEqTn3m+drV+DkZ4zt/75eOv82dex4mQRvn86SwoJ9BYaRPSw3ZfgDhLCMjg9TUVFkHEGZqG8/yvWfe53/XpjNr4pjeTxhOJVlQXQ6Hgf3RcFALza2ewx9++27ib1jDrJdXkuDnLtetWaOl/OptXDr7mmGotKrmbD3zdtyB0n6GpOixfLr8qZA0cluq3/B5GgHYcsH/9BiQHC4HGjRoO3Yf+/X+Yn5hf47TDt8EfwnaON7I+jUXTZgPQGN7EwnaOL9BweVycdOHv/AJjpE2LbXrOoCKigqvYxIA/JANYcJD5fFGLO8cYMP154dPrv62M7D/L2runew/Q9wE9f2/GOH+38LZzjz7J8cn8sLiFEqWTKb9wkt4+wsPY3j/++yqeS9gV0f65C9gvvi3w/TLdAqHbo5gNrhOl5PPmqqoOPUZ5af2qU8KpyppdDRRe4WF5Bi1e2rDvifYdGgrS8bpyRh/LpnjzyNj/LksHDObF2u2DygghSt/7ZoEAD8kAISe5d0DfPf/3qexpR1z/he4ZUUI929wOqDqddj7J7C/BC1n1Dv989fD7Y+oZZpO03ppMnXxURRfMImSJVPYOUeHU6sha/w8rp10MT/R38qJViXgwKsuekzI7r7DxVA+jThdTvY3V5OWON3z3rW2n/DKyfd8ysZpYmh3OXDi8vn/NCF6DHsj8P9TWO0IJoQ/za3tGJ+z8fS/Pwdg2bxJrFw4tZezhsiZY/DBg2r+/FPH4BBgB9dBDZpWF7W7/0Ly1zapTyeJY7n3/jvZpNnDhNhxfGliBk9PvIgvpWQypUtDMRQDryPJ5LgkzAu/53kaCebnodVovRp/gG1Lf87nTUe9xhQqTu3zrEbuzoWLxvYm7trzMC9e+NOg1S1U5AnAD3kCCI1Pqxr4+qM7+OSwgkYD9163iIIbFhMdNch8PV0SsP3rZAW/OVDCvXPWkDOxI9dN16RnjrbOzVPOHIVfzIB9LlwHQdMldfLhCXFYlkzh6l/9m3kT1KeTilP7aHa08oUJC3ocaBwpfcsj1UeNdpa8c2ev5W6flsO6GVdziW6BZ/whnMkTgAhbL713kLueeJemVgeTJ8TzpGEZV5wfhDv/bgnYVnV8sb+0s0xUHFz0SzXbZrsWbixVV9aOmc6pmjTGf/45GuBgUjyWCyZTsmQKe+ZMI2dSBjnazrGJrjNRetJ11ovSfobxUZ1z7kXoLR47lxsmLw84LdXt2WOlPHuslCmxSayfeR0/SfvqMNYyOCQAiLCQNCaO5jYHly+cwpN3LmNKsHbo6ksCtpaz8KcNUAmuQ3BwRglzLl4DwO6bv8v2sQ9iuWAKTQsv4OpJF/PAxCyW6xYRox34P5+h7OoQg9PbtNQJ0Yk8eJ6B1+s/4G8n3uN4az1numzC0+JoZeuJnVw1MYvx0WE2W60bCQAiZJrOtpMYp/4JXnH+VF750ZXDl6K5DXUxViVq337H5B0NsOeNYk8AuCBnHR8umotl0kXMjA9uqok1U1dG5GyS0aCnsRrzwu+xZupK7pjxZVqdbbxZ/xFzEzqfVl+r28VXPvoVsZoYrky5kBsmL2f1pEu8xoLCRfh3XIkRx+Vy8fS/P2fxD7dhP96ZG+bSBVOGp/E/ATwNlAJ2oB0OJ8fxwBVzyC1Yzd6bO9d+JEbFY5h5bdAbfxH+1kxZyQ2TlxPV0b8fpdFy4+QVXkE7VhtDTkoG5ySmet5rcbZyXuIMWl1tvHqyjPzdDzHtza9w6fvf58EDFk60KsP9qwQkg8B+yEKw4HE4nV6bsJw/K4nvP1uG5V01F8LdV83nV2uHaOOR1kZc7/8Mzbu/hVOA+ybNATyLuh1imvpVfMODZJ2zlrmJ04amLiIiDXRaqsvlYs+ZQ/ylZicv1+yg/NQ+z7FPlj3BwrGzAXUh2tiohD6vcxnIeg1ZCNZPMgsoOLaWHfbZhjFKq8HhdBGl1fCzvCV896oFaLVBXuTVeBg+fhjXvx9DYzsDn6OmXriFzmfe0x3vdfzoymv+Strsa4NbDzEiBGOR3OGWGrbWvENZw16eOf9eT4N/84f/y/sNn/Jfk5dxw+TlLE9aFDCpXddgpIse2++1CDILSAybrWWHuW3Tdp9Uww6n+s6Pb1zM965ZGNwfWlMBtl/j+lcxmo9daKq7HEsEmgD3zodjvU/tPj9cCLdgjNXMjJ/Md2Zd7/Wew+XgjfoPqWlVeOjQyzx06GUmxUxg9eRLuGHycq5MXkp8VCwwdLulyRiACDqH04nxuQq/eeZBvel+8vXPcTgDT7EbiKanvwk/24Km1AXV0BalwZo+DW5E/Rrb2xWEGD5RmijsK57lpSX3cdu0bHTRYznR1sBTVf/g2l0/YVXFjzxl3buluaeldt0tbTAkAIig27n3hFe3T3cuoKquiZ17Twz8h7Q1wX8eg7rPPG9VL7wezsDJsTE8vfoLvPqchZUPvw3T43u+VlS8uhhMiGE2JjqBG6as4I+LN1BzeTGlGQ/w7ZnXMT0uhS9NzATcGUl/73Oue7e0mrP1A/750gUkgu5YfeDGv6tqxf9y+x41VcMHm3Bs+z1RtjMcnT2D6U8cBmDuNT/lhdONLMj5Ft+Y2KV76Za9uJpPcO++It6s/9iz+nZl0gX8Zl6+90pgIUIkRhtNdko62SnpbJq/nlZne2fXT7vvvyn3bmmD6QqSACCCqrbxLI+X7uu9IDC1P4u9aj/GtfNXuP5ajPZjJ1Gn1bd19cdoOF3LhLEpaLRavrLWTybNcbPQjJvFhgkP8FSXGR2vZG2EMJybLYRWoyU+Kpb/NO7n5ZodAcu5u4I+OX2ARWPn9P/nDKKOQnh5Z98Jlv/kVcoqa3sspwFSkxNZNm9Sn67b9nQ23HUBmp+8gPYdJ5yGhjFRvHDdMt7541bGjelbI+5efTs1Vv2vrL4V4W7R2DleaxG6c69NGEjjD/IEIILA6XTx21d2c/+LH+FwukibMo47rjiHH7+wC8BrMNg94dN0a0bgRV+Os6CNVfPxANW79jNzj3ro82ljKLvxJi66+T6+kpTW77rK6lsRSfqyW9pg8kjJE4AfVVVVrF69ms2bN4e6KmHvxKkWbvzNv/l5yYc4nC5uXjaH7b/4Mt+9egF/uvtSpicnepWfnpzIn+6+lOuzZvperPkkjjcKaLs7CfvLBZ63Y9YXsX3pLLb+2sS0V2pYe+ezpA2g8RciErnTUri6zavrawrxzZs3s3r1aqqqqnyOyUIwP2QhWN/tPdrAZff9Axfwm9syue0yvdeqRofTyf27/sZDdX/k/yV/jf9Zeq3vnX/9Xppf/RGx27YRtdcJ7XAybSwT/9KIECI4KcRlIZgICpfL5Wnk502fwFN3LSdtylgWzNB1FurIwd/QeorXj/+Oue3NvHb8Ib5zbDzJseMhLgWX8jmnN3+PcW/8h4QjnacemZqI7Uu5XOd0ohmO3EBChLmhSiEuAUD0S7XSjKHoHe5dfT4r5qsJ0q7NmOFdqEsO/mTAa6nKgX+r/42K5+zfHYyzq7usuDTwweLp1N72Qy7LWc/qjhWQQgjVUKQQj6gAYLVaMZvN5OTkoNfrKS0tJSsri9zc3B7PKywsRKfTAaAoChs2bBiG2o48r318jG89vpOTjWc5dPIM5Q9c438gN1AO/nrU1bgxgKOF6qyLmVb1Pu9dkUXytx5g6bwrhvg3ECKyBXsSQ0QFAEVRsNlsWCwW9Ho9RqOxz41/fn4+oAYRg8GA2WwejiqPCO0OJ/e/9BEP/nU3AOfP1PHsd1b0LXWzCzXf/seoG6lfCpyvHpr6tUJO330ul02QDJxChEJEBQCAiooKz918X2zcuJH9+/d7XmdnZ5OTkyMBoI+O1J7hG4/u5N3P1LQN3/ziOWy8JZ2E2MB/Oi6XS53uWQW8C9R0Oah0fhufMJZ4afyFCJmICwD9YbfbURTFb8CwWq1kZ2cPf6UiyKGTZ1jxk1epP9PK+IQYNt1xETdePLvnk1xODv79x8z5G+odP6h/ZQtR7/wndBatbDpKGkO0F4AQolcRFwCKi4tJTk4GoKysDJPJFLCs3W73+75Op0NRlKGo3ogyMyWRKxZNZX/NaZ5Zvxz9lHG9n6TRMvHFf6uNvxa14c9ATcfcjaRgFiK0IioA6PV60tPT0ev1ANTV1ZGXl0dJSUm/rpOcnExdXV3A4+6FYG6jaWewAydOo0uMRTcmFo1Gwx++dTExUVriYvxvUgFw6sPncKYsQDdDvZtvWPcdxj79IFyE1x2/EGL4uHcCc/O3ECxkAcBisbBly5ZeyxUUFJCerjYs7v+6rVmzBoPBELCbJ5CeGn+A1NTUUbkQbGvZYdY/9S6XL5rKn76zAo1Gw9j4mIDlmw6+RfNvbiZlRzVHMlPRFamT+VMvuQWOPDhc1RZC+NH9xrXrTa1byAJAbm5urzN4urNYLF7nuBt9u93uExwAz5NCd4qiBDw2GrW0OvjxC7swW9Usnsfqm2lsaWd8gv/Gv63uM2ofXM2U0k9J7MjoHH+4nvazLUTHdeTWj4r3PxXUTXLwCxFyEdMFpCgKeXl5VFZWehpvdz9+oMZcr9ej0+mw2+0+ZWQAWFV5vJHbH3mbDw+qm0p875oF3HfTEmKifad4uprqqPrDaqZv3cHUBvW9Jl00n35jHRfe/jDaqI4/p3Gz4Ja90HKSutZT3PDhz2lsb2Z8dCIvLblPXQksOfiFCLmIWWev0+kwm81eDXlRURG5ubleTwKFhYVe5xUUFGC1Wj2vLRaLZ03AaGd59wCX/uRVPjxYT8q4OF78weX88ualfht/gCM/W86MP+5A2wCtiRp23X490VaF9Dse7Wz83cbNgknpJKdezvr0n3Fs/GzWp/+M5NTLYVK6NP5ChIGISganKApFRUWe17W1tV6zgIqKijCZTFRWVnqd13UlcGVlZY8zh2B0JIM7c7adDOPfqKprYtm8STx913KfzJ24nLSdriZmnDpb5+T+dxi/djmfXnkxczf8hXETpoSg5kKIgfDXrkVUABguoyEAAOzcW8NrHx+j4IbFREd53/Uf3/EIY36/AUcLTNjWuR1d69kzxMaNGe6qCiEGSbKBjnLPv21Hq9HwleVzAVg2bzLL5k32KlO3+6/w69uZYqsHp/regXf+xJxLbgOQxl+IEUQCgB/udQCROP/f4XSyc+8JqpVmpuoSWDZvEi1tTr7/bBnPv72fxNgoLjpnos+irtNHKzi98Qam7DiMRk3QyfG08Zy6dxPndjT+QojI414PEFbrAMJZpK4D2Fp2GONzFVTVdXbZTB4fT5RWwzGlGa1Gww+uW8TsSd538Ydfe5AZBT9kbMeUzoZpcRz5zo9ZtPrHSC+/EJHNfSMbVusARHBtLTvMbZu2031Ap+aUOhdfNyaWzfdcyor5vk36lEvuoD32XtrjtXx2+x0s/sbjLJKNWIQY8eRf+QjgcDoxPlfhafzbZu3nzE1baJ91wFMmISaKS86bhMvpYP+zt6PkjqWt6RQAsYlJVP9hK9H/OsUF3yySXbiEGCXkCWAE2Ln3BJrTh1iS2Igr9izN6a+Dph2WvkrCqSvQtMZR2zSOd/78Q9K3PMLcQ2on/76Hrue8/1Z36Jq55LpQ/gpCiBCQADACNNR8zq4Lf0i8tmP0trrLwflvQwM43wftr9W3XFrYv1xPytcfGfa6CiHChwSAESBJe6qz8e/KAewA9oDWqW7OdfTCSVDwHPqFOcNcSyFEuJHO3ghXXnmS372y2//BKKABdT7/LDhi+jGpf6ohVRp/IQTyBBCxXC4Xj/zjU+4r/oCFcc3gTq1zEJgKxHW8XgacBabDzKwbQlFVIUSYkicAP9wLwbpuphBOahvPsuZ3b/Lfm3fR7nCRdWEcNAGlwN+B97oUTgE6Nt6qbDo67HUVQoTW5s2bWb16tSwE66twXwh29/+9xz8+OEpcjJYHvrKEnLLfw0uod/oaIAa1w1/jfZ5swSjE6CMLwUaY/12bzrH6ZjauOMai36cz7kDHxispwBXApFDWTggRKaQLKAKcONXC5h37Pa/nTh7L4/oSLr73G4w70IIrCrgYuAlp/IUQfSYBIMxt33OcZT9+FUPRO7z28THP+ymXfhNXFJyaGcexRx6BrHh11k8gsgWjEKIb6QIKUw6nk8Ktn/DAX/6D0+UiM9VJyqvfhMV/B2Di3Es4Yt7C9CU3MD46BpZcJ1swCiH6RQJAGKpWmvnmYzt5a89xAH533j+5/V9/JqbeyV7dLczLfx6AGRlrOk8aNwvGzSIZWB+l4Z5PH6Vg/nqSp14Wgt9ACBEJJACEmdf/c4xvPraTk41nSUuoYWvLr5i95QQAjkQgbmyv11gzdSVrpq4c4poKISKdBIAwU600c7LxLL9OeYZ175YSdVp9/2jGZBI3vs68aYtCW0EhxIghg8B+DPdCMKezM4v/LSv07Iq9lzutauPfNkHDvvt+yPRnjqOTxl8I0U89LQSTAOCHeyHYcGwH+equKlbc9yonG1s877kuX41LAwdX6ml/pYrz8n495PUQQoxMa9euZdu2baSmpvockwAQIq3tDgqet7Hmd28yvsbK67+9xXPs3K8+xrHntjH7kUoSJkwLYS2FECOZjAGEwIETp/nGH3bwgf0om6N/xTVln0MFHL+5lCnz1Uyd0xfLBi1CiKElAWCYbS07zPqn3mUl/+BvR59lzHEHAI0z43C0NYe4dkKI0UQCwDD683Y7G578Jxbn/Vzy0VE0TnDFwr6bszn3+39nXHRMqKsohBhFJAAMo1XzE8g98B3i65wA1M4bh+NX25h33uWhrZgQYlSSQeAgKa5+k2lv3ExJ9Vte7+/YW4PLpU7znDxpGscXzsCRCHvvvoUUyykmS+MvhAgReQIYjMZD0HKSE2cbeHjXz5nmauH3Zfdx+dKfMk47ho2vVhP//u859IUc1n7jRwBM2riD0031zJu+OMSVF0KMdhqX+/ZUeGRkZJCamurZSMGvxkPw/DxwtPg/3gTOHaD9HE5PiyHxlUa0MXH+ywohxBDZvHkzmzdvpqqqioqKCq9j8gTgR592BGs56b/xdwF7gZ2gPQsuDdSdN4PYtiZiJQAIIYaZ7Ag2BBxOp2/6/VPAm8CRjtcT4ch372HWDQ8NZ9WEEKJPIioA5OXlkZOTQ2ZmJjqdzuuYXq/3e47VasVsNpOTk4Ner6e0tJSsrCxyc3MHVZcPD9SR3vWNGmAr0I66MUsWcAGcSL2GmYP6SUIIMTQiKgDYbDYsFovP+7m5uZSUlPg9R1EUz3l6vR6j0Tjoxh/gk4Yj3gFgIqADYoGVHd8Du5Vu5YQQIkxEVADIzc3FZDJ5vVdUVER+fn6P51VUVPg8MQzWogkzvN/QAtcC8YCm8+2Fum7lhBAiTETUOgCDweD12mq1kpmZGZK6LJmT7PtmAl6NP8CS2UnDUh8hhOiviHoC6NrPrygKdrud7OzsXs8rLi4mOVltsMvKynyeIgYiStu32NnXckIIMdwiKgB0ZTQaMZvNvZbT6/Wkp6d7gkddXR15eXkBxwygc0MYN7/rAeInQlR84HUAoB6Pn9hrHYUQItjc8//d/G0IE7KFYBaLhS1btvRarqCggPR072FUm82G2WzuUwDoTlEUkpKSqK+vDzgusHr16t7XAYDXSuAbdv2cJlcLiZp4Xl76UybFTVAb/3Gz+l1HIYQINn/tWsieAHJzcwc8G8dsNpOWltanshaLxevnuBt9u93uE1j6bdwsGDeLScB3o7Xc8+mjPDx/PZOmXja46wohxDCIyA5qq9Xap1k9iqKQl5eH3W73eg8CrxsYqDVTV3Ls8i3kjbDGf7j2RY408rkEJp+Nf+H4uURkALDb7X4bcLvdTmFhoee1TqfDbDZ7lS0qKiI3Nzfo00JHqnD8ow0H8rkEJp+Nf+H4uURkANDr9Z5ZPV25V/12tWbNGgoLCz1ftbW1PQ4AD0aw/gcH4zrh9scWTr9TOH024fQ7yecy9NcJhmDWJSIDQGVlpd/++/z8fCorK73e0+l0bNiwwfMVjCmggYTTH1s4/cFCeP1O4fTZhNPvJJ/L0F8nGIJZF0kH7ceiRYv6PMjcVVVVFampqYP++cG4TjjVJVjXkboM7XWkLkN7nVDXpbKykk8++cTrPQkAQggxSkVkF5AQQojBkwAghBCjlAQAIYQYpSI2F5AYWu71FO5ZVQNJuzHS5eTkUFpaGupqhBWj0eiZQJGcnByUvTciXVFREYqioNPpqKyspKCgIHzWIbmE6GbDhg1er/Pz813Z2dkhqk14Kikpcck/n0719fWu9PR0V319vcvlcrkqKirk83G5XCaTyfOZuFzq55Sbmxu6CnUjXUDCi3sHNXfKDFD3YbBarV4pNUYzRVGoq6sLdTXCysaNG7n55ps9d7bp6enydARs2bLF625fp9N5/dsKNQkAwkd5eblXY+9OpRFOf7ihVFxczJo1a0JdjbBSWFhIdnY2drsdq9UK0Ke9OkY6vV5PTk6O599OoDQ2oSIBQHjR6XTU19d7rbR2/4MOpz/cULFardKwdeO+WbDb7SiKgl6v9zw1jnZPPPEEdXV1JCUlYTQa/aarCSUJAKJXGzduxGw2h8/AVQi5GzjRyR0AdDqdZ/Mlk8lEXl5eiGsWejqdDoPBQG5uLoWFhZSUlITVk7QEANEjo9GIwWAgPz8/1FUJOXcmWeHNnZix6/7c7r7u0f4UYDQa0ev1lJSUUFlZSV1dHRkZGaGulocEABGQxWIhLS1NGn/UXei6NnCiU9ftVrvS6XSjeuKAu0vM3WWo1+upqKhAr9djsVhCXDuVrAMQfrnv3NyNv3vmy2jt/qirq8Nms3k+F/f6iMLCQvR6/ah+MtDpdOj1ep8BTkVRRnXQtNvtfrtNw6lrTJLBCR/uhs7dqCmKwpYtW8JrAUuI2Ww2MjIykH8+KovFQllZmSfdusViwWw2j/qpoDk5OZSUlHj9uzEYDGEzECwBQHhRFIW5c+f6HaiSPxWVxWJhy5YtWCwW8vPzycvLk5lBdK54BaitrR3SvTcihaIobNy4kZSUFM+4SH5+ftjcSEkAEEKIUUoGgYUQYpSSACCEEKOUBAAhhBilJAAIIcQoJQFACCFGKQkAQggxSkkAEEKIUUoCgBBCjFISAIQIM3a7HaPRGOpqiFFAAoAQg6QoCgaDgbS0NDQaDTk5ORgMBq8y7syq7uM9ZYM0m83k5OQMdbWFkFQQQgSLwWCguLiY+vp6v8fdd/YlJSU9XicjI4OKigqf9202G3l5eZ5MpEIMljwBCBEkVqu1x/THdrvd58mgu572HdiyZcuoTccthoYEACGCQFEU7HZ7j103paWlvWYNNZvNAYOE1Wr12qtZiMGSACBEELg3ihlsWujy8vKAjbzNZpOxARFUEgCECILS0lLPpuiBpKSk9HgNq9XqE0CsVisGg8HT8JeUlGAwGEb1VosieGQQWIggSEtLQ6/XB9wBqy9PCHl5eZhMJr/9/EajEavV6ndwWIiBkicAIQYpWP3/iqIEHOS12Wyy65gIOgkAQgxSMPr/LRZLj5uFW61W6f8XQRcd6goIEenKysoAAvb/2+32Xvv/zWZzwPUBNpsNGPwAsxDdyROAEIPUU9cN4Nk8vqfzgYAbhXef/ukuL8RgSQAQYpAyMjKoq6vze8xqtaLT6QI27gDFxcU9LhDrPn5QVFQ04LoK0ZUEACEGKT8/n+TkZJ/8PkVFRdhsth7v/kGd2pmbmxvwuE6nIy0tDfA/VVSIgZIxACGCoKKiAqPR6FkP4E4Q19vKXbvd3uPTAUBBQQEbN24kOTm517UGQvSHrAMQIoQKCwtJT0+Xu3oREhIAhAihQJk/hRgOMgYgRIj0lPlTiOEgAUCIENmyZUuv6aGFGEoSAIQIEbvdLgO6IqRkDEAIIUYpeQIQQohRSgKAEEKMUhIAhBBilJIAIIQQo9T/B3BxhenNyGOTAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from matplotlib import pyplot as plt\n", "import mpl_config\n", "\n", "plt.plot(U_array, e_hf_list, marker=\"o\", label=\"RHF\", linestyle=\"--\")\n", "plt.plot(U_array, e_ccsd_list, marker=\"d\", label=\"CCSD\", linestyle=\"--\")\n", "plt.plot(U_array, e_uccsd_list, marker=\"s\", label=\"UCCSD\", linestyle=\"--\")\n", "plt.plot(U_array, e_fci_list, marker=\" \", label=\"FCI\", linestyle=\"--\")\n", "plt.legend()\n", "plt.xlabel(\"$U/t$\")\n", "plt.ylabel(\"$E/t$\")\n", "plt.savefig(\"hubbard.pdf\")" ] }, { "cell_type": "code", "execution_count": null, "id": "9640bf24-b832-4b23-8f81-066b1a090c1c", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.13" } }, "nbformat": 4, "nbformat_minor": 5 }