legacy/PYTHON/real_ipt-text_v3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP1: The Interaction-driven Metal-Insulator Transition\n",
    "\n",
    "## GOAL:\t\n",
    "\n",
    "1) To understand the difference between the electronic spectral functions of a metal and of a Mott insulator. \n",
    "\n",
    "2) To understand the interaction-driven metal-insulator transition in the Hubbard Model (the reference model for this problem), known as the **Mott Metal-Insulator Transition** (MIT).   \n",
    "\n",
    "To achieve these goals you will have to understand how to read Green's function and connect them to physical properties. We will use the (exact!)  Dynamical-Mean-Field-Theory (DMFT) solution of the Hubbard Model on the infinite-dimensional Bethe lattice. The DMFT is numerically implemented using the Iterative Perturbation Theory (IPT), which provides (approximate but good) Green's functions. To this purpose you will have to run a PYTHON code, manipulating the input/output in order to simulate different physical situations and interpret the outputs. \n",
    "\n",
    "It would be interesting to study also the details of the DMFT implementation and the IPT impurity solver, but it will not be possible on this occasion. Use it then as a black box to perform virtual theoretical experiments.\n",
    "\n",
    "You will have to provide us at the end of the Tutorial a complete report, answering to the questions of the section EXERCISES. Include graphs and commentaries whenever asked (or possible) at your best convenience.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Real frequency IPT solver single band Hubbard model\n",
    "===================================================\n",
    "Here it is the IPT code.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Author: Óscar Nájera\n",
    "# License: 3-clause BSD\n",
    "\n",
    "from __future__ import division, absolute_import, print_function\n",
    "import scipy.signal as signal\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.matplotlib.rcParams.update({'axes.labelsize': 22,\n",
    "                                'axes.titlesize': 22, })\n",
    "\n",
    "\n",
    "def fermi_dist(energy, beta):\n",
    "    \"\"\" Fermi Dirac distribution\"\"\"\n",
    "    exponent = np.asarray(beta * energy).clip(-600, 600)\n",
    "    return 1. / (np.exp(exponent) + 1)\n",
    "\n",
    "\n",
    "def semi_circle_hiltrans(zeta, D):\n",
    "    \"\"\"Calculate the Hilbert transform with a semicircular DOS \"\"\"\n",
    "    sqr = np.sqrt(zeta**2 - D**2)\n",
    "    sqr = np.sign(sqr.imag) * sqr\n",
    "    return 2 * (zeta - sqr) / D**2\n",
    "\n",
    "\n",
    "def pertth(Aw, nf, U):\n",
    "    \"\"\"Imaginary part of the second order diagram\"\"\"\n",
    "    # because of ph and half-fill in the Single band one can work with\n",
    "    # A^+ only\n",
    "    Ap = Aw * nf\n",
    "    # convolution A^+ * A^+\n",
    "    App = signal.fftconvolve(Ap, Ap, mode='same')\n",
    "    # convolution A^-(-w) * App\n",
    "    Appp = signal.fftconvolve(Ap, App, mode='same')\n",
    "    return -np.pi * U**2 * (Appp + Appp[::-1])\n",
    "\n",
    "\n",
    "def dmft_loop(gloc, w, U, beta, D, loops):\n",
    "    \"\"\"DMFT Loop for the single band Hubbard Model at Half-Filling\n",
    "\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    gloc : complex 1D ndarray\n",
    "        local Green's function to use as seed\n",
    "    w : real 1D ndarray\n",
    "        real frequency points\n",
    "    U : float\n",
    "        On site interaction, Hubbard U\n",
    "    beta : float\n",
    "        Inverse temperature\n",
    "    D : float\n",
    "        Half-bandwidth of the noninteracting system\n",
    "    loops : int\n",
    "        Amount of DMFT loops to perform\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gloc : complex 1D ndarray\n",
    "        DMFT iterated local Green's function\n",
    "    sigma : complex 1D ndarray\n",
    "        DMFT iterated self-energy\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "    dw = w[1] - w[0]\n",
    "    eta = 2j * dw\n",
    "    nf = fermi_dist(w, beta)\n",
    "\n",
    "    for i in range(loops):\n",
    "        # Self-consistency\n",
    "        g0 = 1 / (w + eta - .25 * D**2 * gloc)\n",
    "        # Spectral-function of Weiss field\n",
    "        A0 = -g0.imag / np.pi\n",
    "        # Clean for PH and Half-fill\n",
    "        A0 = 0.5 * (A0 + A0[::-1])\n",
    "\n",
    "        # Second order diagram\n",
    "        isi = pertth(A0, nf, U) * dw * dw\n",
    "\n",
    "        # Kramers-Kronig relation, uses Fourier Transform to speed convolution\n",
    "        hsi = -signal.hilbert(isi, len(isi) * 4)[:len(isi)].imag\n",
    "        sigma = hsi + 1j * isi\n",
    "\n",
    "        # Semi-circle Hilbert Transform\n",
    "        gloc = semi_circle_hiltrans(w - sigma, D)\n",
    "\n",
    "    return gloc, sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x128459970>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEQCAYAAACz0c/rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABGRElEQVR4nO3dd3xUVdrA8d9JnfRKQiCk0HsJSFEEREQEleKioq4NF/TVdX3XLbLW9V131XV1XXWxICtid0VFRSxUQVoCgZAQSAglgUAaCQnpmfP+cScYMCFtJjOTeb6fz3xm5pZzn5vyzJ1zT1Faa4QQQnR+bvYOQAghRMeQhC+EEC5CEr4QQrgISfhCCOEiJOELIYSL8LB3AE0JDw/XcXFx9g5DCCGcSlJSUoHWuktj6xw24cfFxZGYmGjvMIQQwqkopY40tc4qVTpKqaVKqTyl1N4m1iul1L+UUplKqT1KqQRrHFcIIUTLWasO/y1g2gXWXwX0sTwWAIutdFwhhBAtZJUqHa31RqVU3AU2mQm8rY1uvVuVUsFKqSitda41ji9ER8gtqeBfazLZmlWIr5c7Vw/txh2XxGHydLd3aEK0SEfV4XcHshu8z7EsOyfhK6UWYHwDICYmpoNCE6J5SUeKuP0/O6iuNTOxbxeKy2t4ZnU6q1NPsPS2UYT5e9s7RJdSU1NDTk4OlZWV9g7FbkwmE9HR0Xh6erZ4H4e6aau1fh14HWDUqFEyyI9wCFn5Zdy2dAddArx5646LiA3zA2D13hP85oNdLFyexHu/GouXh7Ry7ig5OTkEBAQQFxeHUsre4XQ4rTWFhYXk5OQQHx/f4v066i/0GNCjwftoyzIhHFp1rZl739uFp7vi3bvGnE32ANMGd+W5ucNIPHKKF74/YMcoXU9lZSVhYWEumewBlFKEhYW1+htORyX8lcCtltY6Y4ESqb8XzuDtLYfZl3uaZ64bSrdgn5+tv2ZYN+aOjOaNjVnsP1Fqhwhdl6sm+3ptOX9rNct8H9gC9FNK5Sil5iul7lZK3W3ZZBWQBWQCbwD/Y43jCmFLxeXVvLgmg4l9u3DFwMgmt1s0fQB+3h78/Zv9HRidEK1nrVY685pZr4F7rXEsITrK8i1HKK2sZdH0/he8mgr182L++Hie/+4Ae4+VMLh7UAdGKUTLyV0mIRpRWVPHsi2HmdSvC/27Bja7/e2XxBFg8uC1jVkdEJ0QbSMJX4hGfJ58jIKyahZM6Nmi7QNNnswd2YPVe3PJL62ycXTCWfz6178mNjb27PuKigomTpxIXV1di/avrq5mwoQJ1NbWWiUeSfhCNOLjxBx6dfFjXM+wFu9z05gYauo0HydlN7+x6PQOHz7MunXrqK6uprTUuKG/dOlS5syZg7t7yzrreXl5cfnll/Phhx9aJSaHaocvhCM4XHCGxCOn+MO0fq1qCdE7wp8x8aF8tCObeyb2cvlWJB3lz1+kknb8tFXLHNgtkMevGXTBbQ4ePMjYsWPx8/MjODiYo0ePEhISwq5duwgMDOTxxx/nkUce4Y033iA1NZWxY8fy7rvv8t57750tY+XKlSxbtoxPPvnk7LLFixeTlpbGSy+9BMCsWbNYtGgRN998c7vPS67whTjPip05KAWzR3Rv9b6zRnTncGE5qVZOQMLx9OrVi/Hjx7N8+XKSk5MZOnQon332GYGBgaSmprJ3715uuOEGBgwYwN69e6muriYrK4uGw74//PDD/PnPf/5Zufv27Tv7fvDgwezYscMqMcsVvhANaK35Yk8uF/cKIyro5+3umzNtUFce/WwvX+w5Lq11OkhzV+K2lJqayuDBgwHYt28f/fr1A+CRRx7hySefRCnFgAEDSE1NpaCggODg4LP77t69G7PZzODBgzly5AirVq3innvuoaam5pxvh+7u7nh5eVFaWkpAQEC74pWEL0QDB/PLOFRwhjsviWvT/iF+XozvE86Xu3N5aNqFm3MK51ZRUUFlZSUhISFkZ2cTHh6Ol5cX27ZtY/Xq1ezatYt7772XyspKhgwZgo+Pzzk9Y5OTkxk5ciQA3333HRkZGQCkpaUxbNiwc45VVVWFyWRqd8xSpSNEA9+mnQRgygU6WjVn+uAojhVXsC9Xet52ZmlpaQwYMAAwru7rX//pT3/iiy++4PDhwxw+fJjdu3eTmppKSEgIdXV1Z5O+2WymrKyMuro6VqxYQWlpKRUVFbz11lvcdNNNZ49TWFhIeHh4qwZJa4okfCEa+C7tJEOjg9pUnVNvUj9jdrl1+/OsFZZwQA2rc3x8fNi5cycvv/wy1dXVTJky5ex2kZGRlJWVUVRUxNSpU9m0aRMA06dPJysri+HDh3P33XeTmprKqFGjWLBgAQkJP80RtW7dOmbMmGGVmJXRCdbxjBo1SssUh6Ij5ZVWMvqpNTx4RV9+fXmfdpV19Us/4OPpzsd3X2yl6ERDDa+oncnOnTt54YUXWL58eYv3mTNnDk8//TR9+/b92brGfg5KqSSt9ajGypIrfCEsfjhQAMDkARHtLmtyvwiSjpyiuLy63WWJziMhIYHLLrusVR2vZs2a1WiybwtJ+EJYbD5YQKifFwNaMJRCcyb1j8CsYcOBfCtEJjqTO++8s1Udr2699VarHVsSvhAYzTE3ZxYwrlcYbm7tb1kzLDqYQJMHWw4WWiE6IaxDEr4QwMH8M5w8XcUlvcKtUp67m2J0fBhbsiTh24qj3n/sKG05f0n4QgA/HjTq7y/p3fKxc5ozrlcYRwrLOV5cYbUyhcFkMlFYWOiySb9+isPWts2XjldCAJszC+ge7ENMqK/VyqwfeG1rViFzEqKtVq6A6OhocnJyyM933Xsk9ZOYt4YkfOHyzGbNtkNFXDEg0qo9Y/t3DSDY11MSvg14enq2avJuYZAqHeHysgrOUFxew0VxoVYt181NMSY+VOrxhcOQhC9c3s4jpwBIiA2xetlj4sPILqrgREll8xsLYWOS8IXLSzpyimBfT3qG+1m97JGWD5GdR09ZvWwhWksSvnB5SUdPkRATYpX29+cb2C0Qk6cbSUck4Qv7k4QvXFpxeTWZeWVnr8StzdPdjaHdgyXhC4cgCV+4tPqqFlslfDDuDaQeL6GypmXjpwhhK5LwhUtLOnIKdzfFsOhgmx0jISaYmjrN3mMlNjuGEC0hCV+4tKQjpxjULRAfr5YNZtUWCXLjVjgISfjCZdXUmdmdXUJCjO2qcwDC/b2JDfOVenxhd5LwhctKzy2loqaOUXG2TfgAI2NCSDpS7LJjvwjHIAlfuKykI0WAbW/Y1kuIDaGgrIrsIhlITdiPJHzhspKOFhMVZGrX/LUtJR2whCOQhC9c1s4jp2wynEJj+kYG4O/tIfX4wq4k4QuXdPJ0JceKK2x+w7aeu5tiRIx0wBL2JQlfuKT6AdM6ov6+XkJMCOknTlNWVdthxxSiIUn4wiUlHTmFt4cbA6PaP2F5S42MDcGsYXd2cYcdU4iGJOELl7Tz6CmGdA/Cy6Pj/gWGxwSjFFKtI+xGEr5wOVW1dew9drrDbtjWCzR50i8yQBK+sBtJ+MLlpB4/TXWducNu2DaUEBvCzqOnMJulA5boeFZJ+EqpaUqp/UqpTKXUQ42sv10pla+USrY87rLGcYVoi59muAru8GOPjAmhtLKWzPyyDj+2EO1O+Eopd+AV4CpgIDBPKTWwkU0/1FoPtzyWtPe4QrRV0pFTRIf4EBFg6vBj17cKkmodYQ/WuMIfDWRqrbO01tXAB8BMK5QrhNWZzZpth4oYHW/dCctbKjbMlzA/L0n4wi6skfC7A9kN3udYlp3vOqXUHqXUf5VSPRorSCm1QCmVqJRKzM/Pt0JoQpwrI6+MojPVjOsZZpfjK6VIiA0h8XCRXY4vXFtH3bT9AojTWg8FvgOWNbaR1vp1rfUorfWoLl26dFBowpVsOVgAwFg7JXyAMfGhHC4s51ixDKQmOpY1Ev4xoOEVe7Rl2Vla60KtdZXl7RJgpBWOK0SrbckqJDrEhx6hvnaLYUJf42JmU4Z8ixUdyxoJfwfQRykVr5TyAm4EVjbcQCkV1eDttcA+KxxXiFapr7+3V3VOvT4R/kQGevNDRoFd4xCux6O9BWita5VS9wHfAO7AUq11qlLqSSBRa70SuF8pdS1QCxQBt7f3uEK0VnJOMcXlNYzvE27XOJRSjO/dhbXpJzGbNW5uyq7xCNfR7oQPoLVeBaw6b9ljDV4vAhZZ41hCtNWafSdxd1NM6hth71C4tE84n+zMYe/xEobacAJ1IRqSnrbCZazZl8eo2BCCfD3tHQqX9gnHTcF3aSftHYpwIZLwhUvILion/UQpUwZE2jsUAML8vRkTH8ZXKbkyz63oMJLwhUtYufs4AFcO6mrnSH4yfWgUWflnOHBShlkQHUMSvuj0tNas2JnDRXEhxITZrznm+aYN6oqbgpW7jzW/sRBWIAlfdHp7cko4mH+GOQnR9g7lHF0CvJnUL4KPEnOoqTPbOxzhAiThi07vrR8P4+flzvQhUc1v3MFuGRtDfmkV36bKzVthe5LwRaeWW1LBF7uPc/1FPQjysX/rnPNN7BtBdIgPb27Kkpu3wuYk4YtObfH6g2jgzkvi7R1Ko9zdFPdM6sXOo8WsPyBDLQjbkoQvOq3MvDLe3XaUeaN72HXsnObMHdmDHqE+/H31fmqlLl/YkCR80SmZzZrHV+7Fx9OdB6b0tXc4F+Tl4cYfp/UnLfc0/9l82N7hiE5MEr7olJZtOczmzEIWTe9PuL+3vcNp1owhUUwZEME/vttPxslSe4cjOilJ+KLT2X6oiL+tSmdy/whuGh1j73BaRCnFU7OH4O/tyYLlSZyurLF3SKITkoQvOpWs/DIWLk8kOsSH568fhlLOMxJlZKCJf9+cQHZROQvfTqKyps7eIYlORhK+6DQyTpZyw+tbcVOKpbdfRLCvl71DarXR8aH8fe5Qth4q5O53kqiqlaQvrEcSvugU9uWe5sbXtwLwwYKxxIX72Tmitps9Ipq/zh7C+v35/M87OyXpC6uRhC+c3u7sYua9sRUvDzc+WjiOPpEB9g6p3eaNjuEvswazJj2PhculekdYhyR84dS2ZRVy85JtBJg8+HDBOOKd+Mr+fLeMjeXpOUPYcCCfX72dKElftJskfOG01u/P49al2+kaZOLjhRc71EiY1nLj6BieuW4omzILmL9sBxXVkvRF20nCF05pc2YBC95Ook+kPx8tHEfXIJO9Q7KZ60f14B9zh7HlYCF3vLWdM1W19g5JOClJ+MLpJGcX86u3E+nZxY93548l1M/5WuO01pyEaF64YTjbDxVxx392UCZJX7SBJHzhVLKLyrnjP9sJ9/fm7TtHO8T8tB1l5vDu/GveCJKOnmLB24kyhr5oNUn4wmlU1tRxz7tJ1Jk1b985mojAzluN05Srh3bjmeuG8uPBQp76ap+9wxFOxsPeAQjRElprHv1sL3uPnWbp7aOcup19e/1iZDT7ck/z5qZDDIwK5PqLetg7JOEk5ApfOIUPdmTzcVIO91/eh8n9I+0djt0tuqo/l/YJ55HP9pKcXWzvcISTkIQvHN7u7GIe/zyViX278JvL+9g7HIfg4e7GS/NGEBnkzT3vJFFQVmXvkIQTkIQvHFrRmWrueSeJiEBvXrxxOO5uzjMYmq0F+3rx6i0jKTpTzX3v7ZTJU0SzJOELh1VbZ+b+93dRcKaaV28Z6ZSDodnaoG5BPH3dELZmFfH01+n2Dkc4OLlpKxzW37/Zz6bMAp79xVAGdw+ydzgOa/aIaHZnl7Bk0yGG9gjm2mHd7B2ScFByhS8c0pd7jvPaxixuGRvD9aOkFUpz/jR9ABfFhfDH/+4h/cRpe4cjHJQkfOFwkrOL+f3HexgZG8JjVw+ydzhOwcvDjVduTiDA5MHC5UmUlMuMWeLnJOELh5KZV2b0pA3wYvEtCXh5yJ9oS0UEmFh8y0iOF1dw73s7ZXRN8TPy3yQcRsbJUm5Zsg03pVh+5xgiAlyvJ217jYwN4W9zhrL5YAELZBx9cR5J+MIhJB4u4vrXtlBr1rxz1xiX7knbXr8YGc0zc4byQ0Y+Ny/ZRn6ptNEXBkn4wq7qzJpXNxzkhte3EuTjySf3jGNAVKC9w3J611/Ug5fnJZB2/DTXvryJzZkF9g5JOABJ+MJuth8q4tqXN/H01+lMHRjJyl+PJzZMruytZcbQKP57zzhMnu7cvGQbD360m5xT5fYOS9iR0lrbO4ZGjRo1SicmJto7DGFlZVW1fJ92kre3HGbn0WK6Bpp45OoBzBgShVLSi9YWKmvqeHFNBm/+cAiz1lw7rBu/GBXN2Pgw3KTncqejlErSWo9qdJ01Er5SahrwIuAOLNFaP33eem/gbWAkUAjcoLU+fKEyJeE7P601BWXVZJwsZefRU+w4fIotWYVU15rpEerDXeN7MndUNL5e0v+vI+SWVPDq+oOs2HmM0qpawv29ubhXGON6hTEgKpC+kf7yu+gEbJrwlVLuwAHgCiAH2AHM01qnNdjmf4ChWuu7lVI3ArO11jdcqNy2JnyzWVNuaZlw/rk1fHfOKt3wpW50m3P31U0sb7ycJl42uX1Ljnv+dk0tt9b5gKayxkxVrZmq2jqqaoznyhozlTV1lFTUUFRezakz1RSdqebE6UoOF5SfMzNTnwh/xvcJZ8aQKBJiQuTq0k4qquv4Nu0E69Lz2JRZeHbgNaUgKtBEZJCJroEmIgNNBPl4EmDywN/bA3+TB37eHni5u+Hp7oanu7I8//RaKVBKoSzluVleo0ChjPX8tI2bZUH9ckfgCN803ZXCx8u9TfteKOFb4+N8NJCptc6yHOwDYCaQ1mCbmcATltf/BV5WSiltg/qkU+XVjPzL99YuVrSAl7sboX5ehPh50SXAm5ExIcSF+9Gziz/Do4NdanYqR+bj5c7M4d2ZObw7WmuOFpWTfqKU/SdKOVJYzonTFRw4WcqmjAJKZSpFuxjeI5jP7r3E6uVaI+F3B7IbvM8BxjS1jda6VilVAoQB5zQdUEotABYAxMTEtCkYXy8PHp4+oEGZTW/b8JNcnbOcJparZrehRWU2X05T259/GdRkfE2V1YI4zim/wQpvDzdMnu54e7id89rk6U6Qjye+Xu4OcXUkWk4pRWyYH7Fhflw5qOvP1pvNmjPVtZRW1lJWZTxqas3UmjXVdWZqas3U1GlqzWaqa81obXyrNJ6Nb4lmbfmeaXnW2vhW2XC9o3CUULoEeNukXIeqsNNavw68DkaVTlvK8PFy51cTelo1LiFclZubIsDkSYBJvp11BtZolnkMaDi6VbRlWaPbKKU8gCCMm7dCCCE6iDUS/g6gj1IqXinlBdwIrDxvm5XAbZbXvwDW2qL+XgghRNOs1SxzOvBPjGaZS7XWTymlngQStdYrlVImYDkwAigCbqy/yXuBMvOBI+0OzrbCOe8+hAtx5XMH1z5/Vz53cPzzj9Vad2lshcN2vHIGSqnEppo/dXaufO7g2ufvyucOzn3+MrSCEEK4CEn4QgjhIiTht8/r9g7Ajlz53MG1z9+Vzx2c+PylDl8IIVyEXOELIYSLkIQvhBAuQhJ+Gyil5iqlUpVSZqXUqPPWLVJKZSql9iulrrRXjB1BKfWEUuqYUirZ8phu75hsTSk1zfK7zVRKPWTveDqaUuqwUirF8vvu1OOXK6WWKqXylFJ7GywLVUp9p5TKsDyH2DPG1pKE3zZ7gTnAxoYLlVIDMXoaDwKmAf+2DB/dmb2gtR5ueayydzC2ZPldvgJcBQwE5ll+567mMsvv2ynborfCWxj/xw09BKzRWvcB1ljeOw1J+G2gtd6ntd7fyKqZwAda6yqt9SEgE2P4aNE5nB0KXGtdDdQPBS46Ia31RoyRARqaCSyzvF4GzOrImNpLEr51NTZUdHc7xdJR7lNK7bF8/XWqr7dt4Iq/3/Np4FulVJJlOHNXE6m1zrW8PgFE2jOY1nKo4ZEdiVLqe+DnA4TDw1rrzzs6Hnu50M8BWAz8H0YS+D/gH8CdHRedsIPxWutjSqkI4DulVLrlStjlaK21Usqp2rVLwm+C1npKG3ZryVDRTqWlPwel1BvAlzYOx9463e+3tbTWxyzPeUqpTzGquVwp4Z9USkVprXOVUlFAnr0Dag2p0rGulcCNSilvpVQ80AfYbueYbMbyB19vNsbN7M6sJUOBd1pKKT+lVED9a2Aqnf93fr6GQ73fBjjVt325wm8DpdRs4CWgC/CVUipZa32l1jpVKfURxny+tcC9Wus6e8ZqY88qpYZjVOkcBhbaNRobs0zPeR/wDT8NBZ5q57A6UiTwqWUaSw/gPa31avuGZDtKqfeBSUC4UioHeBx4GvhIKTUfY/j26+0XYevJ0ApCCOEipEpHCCFchCR8IYRwEZLwhRDCRTjsTdvw8HAdFxdn7zCEEMKpJCUlFTQ1p63DJvy4uDgSEzv12ExCCGF1SqkjTa2TKh0hhHARDnuFL4Rov/zSKjLzytBa0y3Yh9gwXyzt6IULkoQvRCeUdKSIv3+zn61Z5w72GBnozawR3bllTCw9Qn3tFJ2wF0n4QnQiWmteWpvJ898dIDLQm99e0ZcRMcF4uruRlX+Gtel5LPnhEG/+cIibx8Tw68v7EO7vbe+wRQeRhC9EJ/LsN/tZvP4gs0d056nZg/H1+ulffGzPMG4aE0NuSQUvr83knW1H+Xz3cR6dMZA5Cd2lqscFyE1bITqJj3Zks3j9QeaNjuH564edk+wbigry4anZQ/jmgUvp3cWfBz/ezW3/2cGJksoOjlh0NEn4QnQCRwrP8MQXqVzcK4y/zBrcoqv13hEBfLRwHE/OHMSOQ0VMe3Ejq/ee6IBohb1IwhfCyWmteeiTFDzcFP+4fhjubi2vmnFzU9w6Lo6v7h9PjxBf7n4niUUrUiivrrVhxMJeJOEL4eTWpuexJauQB6f2IyrIp01l9Ozizyf3XMzdE3vxwY6jXP3SJvYeK7FypMLeJOEL4cTqzJq/fZ1Oz3A/bhoT066yvDzceOiq/rw7fwxnqmqZ/e/N/Ht9JnVmGUK9s5CEL4QT+zb1BJl5ZTw4tR+e7tb5d764dzirfzOBKwZG8uzq/dzw2haOFpZbpWxhX5LwhXBSWmte25hFbJgv0wY3Ns9824X4efHKTQm8cMMw9p8sZdqLG3l/+1FkwiTnJglfuJT80ioe/Wwv0/65kVuWbOObVOdtlZJ45BTJ2cXcNT6+VTdqW0opxewR0XzzwASG9whm0YoU7lqWKM03nZgkfOEyDhWc4ZqXNvHhjmy6BpnIPlXOwuVJPPb5XsxOWE/9ztYjBJo8+MXIHjY9TrdgH96ZP4bHrh7IpswCpjy/gf9sPiR1+05IetoKl1BaWcP8ZTuoqq3j03svZlC3IGrrzDz7zX5e35iFt4cbD88YaO8wW6ykvIav957gxot64OPlbvPjubkp7hwfz+UDInj081T+/EUaK3Ye46+zhzAkOsjmxxfWIVf4wiW88F0GhwrOsPiWkQzqZiQoD3c3Fl3Vn9vGxfLGD4dYlZJr5yhbbuWe41TXmrl+lG2v7s8XG+bHsjsu4qV5IzhxupJrX9nEbz9K5nhxRYfGIdpGEr7o9PafKGXZlsPMGx3D2J5h56xTSvHI1QMZFh3EohUpTlM//XFiNgOiAhnULbDDj62U4pph3fj+txNZcGlPvtyTy2XPreeZ1emUVNR0eDyi5SThi07vpbUZ+Hi684cr+zW63tPdjX/eOILKmjr+76u0Do6u9Q7ml7Enp4Tr7DzgWZCPJ4umD2DtgxO5anBXFq8/yPin1/Ls6nQKyqrsFpdomiR80akdKjjDqpRcbhkbS7CvV5PbxYf78T+TevPVnlw2ZxZ0YIStVz/ezYyhUXaOxBAd4ss/bxzBV/ePZ0LfLizecJDxz6zliZWpZOWX2Ts80YAkfNGpvbkpCw93N+4cH9fstgsn9iQm1JcnVqY6dAuUr/bkkhAT3OZhFGxlULcgXrk5ge9/O5Grh3bjna1HmPyPDdy8ZCtfp+RSU2e2d4guTxK+6LQqquv4fNdxrh4SRUSAqdntTZ7uLLqqPxl5ZXy661gHRNh6hwvOkJZ7mulDHOPqvjG9uvjz3Nxh/LhoMr+b2pfDBeXc8+5Oxv1tDU+sTCXpyCnpwGUn0ixTdFqrUnIprarlhota3pJl2uCuDOkexAvfHeCaYVF4e9i+yWNrfG2pzrnKgRN+vYgAE/dN7sM9k3qzLj2PFbtyeG/7Ud768TDRIT7MGBLF5P4RjIwNwcNKw0KIC5OELzqtD3dkEx/ux+j40Bbvo5Ti91f249al2/lgeza3XRxnuwDb4Ou9uQzrEUz3YMeqzrkQdzfFlIGRTBkYSWllDd+mnmTl7uMs3XyI1zZmEWjyYELfLlzWL4KxvcKc6tycjSR80SkdLSxn++Ei/jCtX6tbslzaJ5yxPUN5eV0mN1zUA5OnY1zl552uZE9OCb9vorWRMwgweXLdyGiuGxlNaWUNmzMLWLMvj3X78/lyj9EPIjrEhzHxYYyJD+Wi+FDiwnxl+kUrkYQvOqVVe43kcc3Qbq3eVynFby7vy7w3tvJxYja/HBdn5ejaZv2BfAAu6xdh50isI8DkybTBUUwbHIXZrNl34jTbDxWxLauIdfvz+GRnDgCBJg+GRAcxuHsQQ7sHMzQ6iOgQH/kQaANJ+KJT+joll6HRQfQI9W3T/mN7hpIQE8yrG7K4cXSM1YYebo/1+/OIDPRmQFSAvUOxOjc3xaBuQQzqFsQdl8SjtSYzr4zEI6dIOVZCSk4JSzcdoqbOuNnr7+1B7wh/+kb60ycigD6R/vSJDKBbkEk+CC5AEr7odHJOlbM7p4Q/Tuvf5jKUUtw3uTd3vpXI58nH+cXIaCtG2Ho1dWZ+yChgxpAol0hoSin6RAbQJzKAeZZlVbV1HDhRxp5jxew/UUrGyTLWpufxUWLO2f18vdyJCfUlJtSX2DBfYsL8jNehvnQP8XGID257koQvOp2vU4yWLNOHtG+M+Mv6RTAgKpB/r89k9ojuNhmCuKV2HjlFaWUtk/p1sVsM9ubt4c6Q6KCfDdZWdKaajJOlZOSVcTC/jKOF5RwqOMOGA/lU1f7U9t9NQWSgia5BJroF+dA1yERUkKnBsw+RAd6dusWQJHzR6azam8ugboHEhvm1qxylFPde1ov73tvFN6kn7Nr2fd3+fDzcFJf0DrdbDI4q1M+LMT3DGHPeOElmsyavtIqjReUcKTzD0aJyjhdXcuJ0BftOnGZteh4VNXXn7OOmjPLC/b0tDy/CGrw+uzzAizA/b7w8nOvDQRK+6FTyS6vYdbSYB6/oa5XyrhocRc/wAyxef5CrBne1W3XK+v15XBQXSoDJ0y7Hd0Zuboquliv4xprmaq05XVFL7ukKcksqyS2uJLekgoKyKvJLqykoq+JI0RkKSqt/9sFQz8/LnWBfL4J9PY2HjxdBvp4E+3gS4vvT67Pb+HgSYPLE5Olml78lSfiiU9lQ35Klv3Vasri7Ke66tCd/+jSFbYeKfjbaZkc4XlxB+olS/jS97fckxM8ppQjy9STI15P+XS886uiZqloKy6rJL6uioKyKwrJqCsuqOFVeQ3FFNSXlNRRX1JBecpqSihqKy2uovcDwHO5uCj8vdwJMnvh7e+Dn7Y6/yZMAy+veEf4smNDL2qcsCV90Luv25xER4G3VYYPnJHTnH9/uZ8kPWXZJ+Ov3d67mmM7Iz9sDP28PYsJa1upLa01ZVS3F5TVnPwCKK6o5VV5DWWUtZ6pqKauqpbSylrKqGs5U1VFSUcOxU+WUVdWSW1IpCV+IC6mtM7PxQL7Vq15Mnu7cMjaWF9dkkJlXRu8If6uV3RLr9+fRPdinw48r2k4pRYDJqL7p2ClqLsy57jgIcQE7jxZTWllrkyvhX46LxdvDjTc3HbJ62RdSVVvH5swCJvXr4hLNMYVtScIXnca6/XlGS5Y+1m/JEu7vzZyEaFbszOnQyT0SD5/iTHWdVOcIq7BKwldKTVNK7VdKZSqlHmpkvbdS6kPL+m1KqThrHFeIhtal5zEqLoRAG7VkmT8+nqpaM8u3HLFJ+Y1Zl56Hl7sbF/fu+HsHovNpd8JXSrkDrwBXAQOBeUqpgedtNh84pbXuDbwAPNPe4wrRUG6J0ZLFllfCvSP8mTIgguVbj1DZRDM9a1t/IJ8xPUPx9ZLbbaL9rPFXNBrI1FpnASilPgBmAg0nB50JPGF5/V/gZaWU0raYBaHyNHz1W8sbS52nUue9bmodP6274Haqme1aU0ZzMbVwnbun5eFleTR8bXl4eIF3EJjqH4Hg4U1ncLYli5WaYzaq8CCPhXzLjswt5Cx7n94jJsKgWcbP0gayi8rJzCtj3ugYm5Tf4cx1UFVqPKrLoKoMqk5D9Rmoq4a6GjDX/PS6rsFrbTYeaMtr3eC9Pu+9ufFtzmok7fwsFekLvv35+ub2b+X60Hi47E8/C7O9rJHwuwPZDd7nAGOa2kZrXauUKgHCgHMmD1VKLQAWAMTEtPGP3FwLx5Ia/AD1ua8bPDW+Tp/3upXb/WydtY7VgjLawsMEvmEQ2A0Cu0NQNATHQMQAiBgEfs5RlbAuPY9uQSb62KIlS+Vp+GYR7HqXGDQ+nuF45SRDzifw7SNw+WNw0V3nfohbQf3omE4znILZDEVZUJgBBRlQmAmnj0PZSSjLgzP5oNv6zchyAaXcLK/dmnjfyLL61+cU19jvSjWzTWvXt6P82spG4ms/h/qeqLV+HXgdYNSoUW3LYr6hcP8ua4blHMx1lqshyxVRbVWDKyXL8tpKI3lVnYbKEqgsNt6fKYDTOXAiBQ6sPvePLSAKeoyB+EshfiKE97HbKTalviXLzBHdrd+SpTgb3p4Jpw7BuHth3L1sOQT3v7+TD6/2ZkzWK7Dqd5C9DWb+2/gWZSXr0/OICfWlZ3j7hoiwGXMd5CTCwTWQvd240Ko6/dN63zDjAiKgK0QNA/9I4//Tyx+8A4yHlz94+RnfNOu/lbqd923VzQPcpH2JNVgj4R+Dc5qaRluWNbZNjlLKAwgCCq1wbFHPzR3cfMCznbMFaQ2lJyAvDU6mGh8CRzZD2mfG+i4DYPB1MHQuhMS1N2qrsFlLltO58NYMqCiG276EuEsAmD7YzDPBvjy318THCz+FTc/DmieN6oM5S6ySnCpr6vjxYCFzR0U7VnNMrSFnB+x8G/avgvJC4+o5chAM+QV0S4Au/SGsl5HchUOxRsLfAfRRSsVjJPYbgZvO22YlcBuwBfgFsNYm9fei/ZSCwCjj0ftyY5nWxlf1g2th7wpY9xdY9xQMuBrG/Rpizq/B61hnW7L0smL1U201fHSr8e3n9i+g+8izqzzc3Zg/Pp4nv0xjZ3YxCZc+CModvn8cwnpbpe51x+EiKmrqHKc6p64W9v4XNr9oXAx4+kH/6dB3GvSeAj7B9o5QtEC7E76lTv4+4BvAHViqtU5VSj0JJGqtVwJvAsuVUplAEcaHgnAWShlXbGG9YPSvoCQHEpcaj31fQN+r4IonoYt1BixrrfUH8hkdH4qftxVrKNc+CTnbYe5b5yT7ejdc1IMX12Tw+oYsXv3lSLjkN5C/Hzb+HeInQNz4dh1+XXo+Xh5ujOtp59ExtYb0L+G7x6HoIEQOhmteNL7leXe+iVg6O6tUjGmtV2mt+2qte2mtn7Ise8yS7NFaV2qt52qte2utR9e36BFOKijauFH5v2kw5QmjyuffY2HtU8aVcQfKOWW0ZLHqlfDxZNjyCiTcBoNmN7qJn7cHvxwbyzdpJ8jKLzM+FKf/3ajm+vxeqGnfTbf1B/IY2zMMHy87zqdbnA3v3wgf3mLUsd/wDiz8AUbeLsneSTnUTVtrKCmvYf6yHXi4K/y9PYyHyQN/b08CTB4E+3oS5PPzR4DJ064TXDglL18Y/78w4pdGa5WNzxo3fee+ZXwb6AD1zTGtlvDNZvjyAfDrYnxruYDbLo7j9R+yWLLpEH+dPQS8/WHG87B8Fvz4Ekz8fZtCOFpYTlb+GW4ZE9um/a0ifRV8drdRlTP1KRhzN7h3unThcjrlb9DLw43aOs3x4krKqoyR6Uora6muMze5j1IQ4O1hDJd63odBoGVs6zDLxAhhlkkRwvy8MHna8QrMUfiFw+xXYcA18Pl98MZkI+n3uszmh16/P5/uwT706mKl5pj7VsLxXTDr1WbrpbsEeHNdQjT/Tcrhf6f0pUuAt3HOA2fCD/+AkbeBf+tvJH+37yQAk23Zp6ApWsPav8APzxkta+a+BaE9Oz4OYROdLuEH+Xry3q/GNrqussYYgrT+UT90af3j9Nnl1ZRU1HCipJKSilpOV9Q0+WHh7+1hfAD41c+MY8yEExnoTdcgH6KCTEQGmgjz88Kts3+D6D/DaK3x/jx45zrjQ2Do9TY7XFVtHT8eLGC2tZpj1tUaN6O79G9x3L+6NJ4Pdhzl7S2HeXBqP2Ph5MeMexs//gum/qXVYXybeoJ+kQHEdXRzzLpa+OI3kPyO8a1t+nPgaerYGIRNdbqEfyEmT3dMnu5EBrbuj1hrzZnqOorKqik489PkB4VnjFlxis5UU1hWTXZROcnZxRSdqabuvMkPPN0VEQE/zaHZNdBEVLDP2cmWe4T42re+1lpC4mD+t0bSX7HA6A+Q8EubHCrx8CnKq+uYZK3mmGmfQcEBuP5to5lrC/Ts4s/UgZG8veUId0/sZdw4Du8NQ+bC9iVw8W/Av+XVTYVlVew4XMR9l/Vu40m0kbkOPplv/AwmPgSTHrJ6RzJhfy6V8NtKqZ/uB7RkAoQ6s6awrIoTpyvJLankREklJ05bnksqST1+mu/3naSy5txvDREB3sSE+hqPMOODoE9EAL26+DvXh4F3ANz0EXx4M6z8tTGEw8CZVj/M+v1WbI6ptXGjNqw39L+mVbsumNCLb1JP8lFiNndcEm8snPB72PMhJL5pJM8WWrMvD7OGqYPaNwF7q2htXNmnfQZX/B9ccn/HHVt0KEn4NuDupogINBERaGJodOPbaK05VV7D0aJy42GZZPloUTlbswr5NPnY2VEUlIIeIb70jfSnT2QAfSL86dc1gL6RAXi6O2gPRC9fuPE9WHaNcaUfEAU9Rlv1EOv3W7E5Zs4OOL7TqMZoZcepkbEhjIoNYckPh7hlbKzxOwnvY7RPT1wK43/b4h6436adoHuwj1Vn7GrW+qdh13LjQ0qSfacmCd9OlFKE+nkR6ufF8B7BP1tfVVtHdlE5GSfLOHCyjIy8UjJOlrHhQD41dcYngZeHGwOjAhkWHcSQ6GCGRQfRs4u/47Q28vSBeR/Am1cYVTx3bzI6dFlBzqlyMvLKuOEiK80ntHWxMQjasHlt2v1/LuvFnW8lsmJnDjdcZBkHavRCeG+ucSN4yC+aLaOsqpaNGQXcPCam43rXpn8FG56GYTfBZQ93zDGF3UjCd1DeHu70jgigd0QAVw35aXlNnZkjhWdIyy0lJaeY3TklfJyUwzLLGO0BJg/GxIcytmcYY3uGMSAq0L4fAH7hRtJ/fRJ8chfc+rlVmvd9n2bFlizlRUbnolHzjaaVbXBZvwiGRQfx0tpM5iREG1f5vacYLVwSl7Yo4a/ee4LqWjNXD+3WphharfAgrFgI3UbA1S9Inb0LkITvZDzd3c5+EFw7zEgMdWZNVn4Zu3NKSDxcxNasQr7flwdAoMmDi3uFc/mACCb3jyDM3w7DIXfpZ7RP/+xuo62+FYYe+Cb1JH0i/OlpjeaYqSuMweWGnz8iSMsppXhgSl/ueGsHnyTlcOPoGKNqaPhNRjPHU0cg5MLt6j/dlUNsmC8JMcFtjqPF6mrh04VGjDe8I61xXIRTJfyamhpycnKorLTN0KHObmSoiZlDB+Dp6UluSQXbsorYcrCQDQfyWZ16AqUgISaEKQMiuWZYFNEhzd+Atprh8yBrPWx8DvpfDVFD21zUqTPVbD9cxN0TrdQ+PPl9Y8iAdsQERuevYT2CeXmdcZXv5eEGQ28wEv6ejy7YEetESSU/Hizk/sl9OqY6Z/M/jfsW171p9JwWLsGpEn5OTg4BAQHExcU51giCDkBrTWFhITk5OcTHxxMV5MOsEd2ZNaI7WuuzLYO+33eSZ1an88zqdEbHhzJ7RHemD4kiyMc20wKeY9rfjAHYPr8XfrWuzVU7a9PzqDNrpg60QkuWggw4lmj0Jm0n4yq/D3f8Zwef7MwxJi4JjoHY8bDnA5jwuyarTT6z3KSfNaJ7u+NoVl66caN20GxjTBzhMhy0iUfjKisrCQsLk2TfCKUUYWFhjX77UUoxuHsQD0zpy5e/vpQf/nAZD17Rl4LSKhatSGH0U9/z+493s/dYiW2D9A2FGc/BiT2w9ZU2F/NN6gm6BpoYGm2Fmab2fGQM7ztkbvvLAib17cKImGBe+O4AZ6pqjYXDbjQmAzm2s9F96syad7cdYXRcKPG27mylNXz9e6MV1fTnpN7exThVwgck2V9AS382PUJ9+fXlfVjz4ERW3ncJ142M5ss9uVz90iauW/wjX+3JxWy20ejVA2caQ+pu+LsxC1IrlVfXsjEjn6mDIq3zt5D+FcRcDAGR7S8L43fwyIyB5JVW8eqGg8bCAdcYk3jUzylwnnXpeWQXVXDbxXFWieGCUj+FQxth8qPGDXXhUpwu4QvrUUoxNDqYv84ewtY/Xc6jVw+koKyKe9/byZX/3Mjnycd+1mPYKqY+BbUVRt12K32XZnRYmz7ECs07Tx2GvFRjXHcrGhkbwszh3Xh9YxY5p8qNMXniJxrDLTQyDcTSzYfoGmhi6iDrfOg0qabCGOSu6xAYdadtjyUckiR8AUCQjyfzx8ez9sFJvDRvBErBbz5I5ornN/B1Si5Wna8mvDdc9Cujs8+JlFbt+tmuY3QLMjE6zgqzKe3/2njud1X7yzrPH6f1Ryl49LO9xs9uwDXGNIkn956z3dasQn48WMj88fG270S3/Q04fQymPdPioSNE5yIJX5zD3U1xzbBurP7NBF69JQFPdzfueXcnN7y2lT05xdY70MQ/gHegMYZ+CxWUVbExo4Brh3e3zkB0+1cZUzbaYDTIbsE+/OHK/qzbn89HidlGyySUcZVvYTZrnl2dTkSAN78cZ+OhkCtLjKkYe085O1WjcD2S8Nvgtddeo2vXrgwbNoxevXrx9ttvN7uPh4cHffv2Zfjw4QwfPpwnnzTGWq+oqGDixInU1dU1W0Z1dTUTJkygtra23efQHDc3xbTBUXx1/3j+NmcIWQVlXPvyZh78aDeFZVXtP4BvKIy7Dw58bQxH3AJf7cmlzqyZNcIKHZMqTsHhzTa5uq93+8VxjOsZxp+/SGNviRfEXnxOwn97y2F2Hi3mD9P6236Y7S2vGOc8+RHbHkc4NEn4bZCSksITTzzB7t27ef/99/ntb3/b7D7PPPMMgwcPZteuXSQnJ/PYY48BsHTpUubMmYO7e/P/8F5eXlx++eV8+OGH7T6HlvJwd2Pe6BjW/W4S90zqxcrdx5jy/AY+ScppfzXPmIVgCob1zzS7qdaa97cfZUBUIP27WmGcmYzvQdcZQzrbiJub4sUbhxPi68Wdb+3gRPcrjPlgCzLZlFHAU6v2MalfF65LsHFTzPIiI+EPnGn0qhUuy6na4Tf05y9SSTt+2qplDuwWyOPXDGp2uz179nDddUb75fj4eLy8jIGxDh06xAMPPMCxY8dwc3Nj+fLl9OtnjJH+4IMPsnPnTv7yl7/w6KOPni3r3Xff5b333jv7fuXKlSxbtoxPPvnk7LLFixeTlpbGSy+9xKxZs1i0aBE333yzVc65pQJMnvxxWn9mj+jOohUpPPjxblbsyuGvs4cQG9bGpoSmQLj4PuPm7fFdF0xGO4+eIv1EqTGzlDXs/wr8I6FbgnXKa0JEoIm37riIX765nXkbw1jnAV9+vIQHsifQs4sfL94wwvYtz5LeguoymPAH2x5HODy5wm+DlJQU+vXrh9aal19+maeeeoqamhruuusunn/+eRITE3niiSd4+umnz9lvyZIlLF++nB9++AEwqmiysrKIi4s7u83DDz/Mn//853P269WrF/v27QNg8ODB7Nixw7YneAF9IwP4eOE4/jJrMHuyS7jqxR/4YPvRtl/tj7Zc5W987oKbvb3lCAHeHswcboXqnNoq4wq/77RWj4zZFn0iA1h53yWMGDKUfTqO7ifWMndUD/57z8UE+dq4w1tdDexYYrQS6jrYtscSDs9pr/BbciVuC9nZ2ZSWljJ9+nSOHTvG0KFDeeKJJ/jvf/9Lamrq2Sv/2tpaLr300nP23bdvH3V1dQwbNgyAgoICgoODz67fvXs3ZrOZwYMHc+TIEVatWsU999xDTU3N2atAd3d3vLy8KC0tJSDAPhNJu7kpbhkby+T+Efzu4908tCKF7/fl8fR1Qwhv7Vg9pkC46C5jSsDCg43OhZtbUsGqlFxuHhNrnaGQD2+C6lLoZ93mmBcSEWji+RuGo9fdCBueYcTUSDB1QO/mfSuNljkznrf9sYTDkyv8VkpJSWHChAkkJydz4MAB0tPT2bJlC7t37+app54iOTmZ5ORk9u7dy+LFi8/uV1VVxe23387ixYsJDDTqoH18fM7pGZucnMzIkSMB+O6778jIyAAgLS3t7IdEfVkmk/0Hu+oW7MM788fwyIwBbMzI58oXNp4dxbJVRi8Ad0/Y+u9GV7+2IQutYf74+HZGbLF/FXj6Qs+J1imvFdSAq1Hon5qE2trWxUYrpD5TO+Z4wqFJwm+lPXv2MGKEUdccEhLCTTfdxFdffUVUVBTffPMNZrMxi1VKSso51RyPPfYYY8eOZerUn/7xQkJCqKurO5v0zWYzZWVl1NXVsWLFCkpLS6moqOCtt97ippuMkRwLCwsJDw/H07MDrg5bwM1NcdelPfnivvFEBJq46+1EFq3Y89OwAi0REGnMIbvrXeMGYwMnT1fy3vajzEnoTo9QKwz2pi3JttdkY7z+jhY5GIJijB6+tpa9wxggbczdHVJ1JRyf/BW0UkpKytmED3DNNdewatUq7rzzTsxmMwMGDGD48OE888wz59yMe+6559iwYcPZZpn332/MLDR16lQ2bdoEwPTp08nKymL48OHcfffdpKamMmrUKBYsWEBCgnFzcd26dcyYYbuWJW3Vr2sAn917MXdP7MUHO7KZ8a8f2HX0VMsLGHef0ft2x5vnLH5mdTpaa+611hyvubuNKo4OrM45h1JGy6Cs9VBVZttjbVts9HVox7DPopPRWjvkY+TIkfp8aWlpP1vm7JKSkvQtt9zS4u1nz56t9+/f3+R6R/gZbTlYoC/+2xrdc9FX+vlv9+ua2rqW7bj8Oq2f7aV1dYXWWuttWYU69o9f6me+3me94Nb+VesngrUuy7dema2VtVHrxwO1Tv3MdscoztH6z6Far/6T7Y4hHBKQqJvIq3KFb2cJCQlcdtllLe54NWvWLPr27dsBkbXd2J5hfP3ApVw7rBsvrsngule3cKjgTPM7jrsXzuTD3k8oKa/hfz9MJjrEh/smW+nqHozmmD3G2HfgsJhx4BMC6atsd4wdS0CbYfSvbHcM4XQk4TuAO++8s8Udr2699dYOiKj9Ak2evHDDcF6aN4JD+WVMf/EHlm85fOHB2HpOgoiBmLe8wn3vJXHydCUvzRuBr5eVGpMVZxtj99iwd22LuHsYTUIPrDaaTVpbdTkk/ceoOgqJs375wmlJwhc2dc2wbnzzvxMYGRvCo5+nMuNfP7Apo6DxdvtKUZ6wALe8VGqzNvLX2UMYERNivWDODpbmAPdA+s+AymI48qP1y97zoTGMwph7rF+2cGqS8IXNRQX5sHz+aP59cwJlVbXc8uY2ZvxrE//ZfIh9uacpOlPNwfwy3tx0iCu/j6BIB/D36M1cf1EP6waS/iWE9TFG67S3XpPBw2Q0EbUmrWHbq9B1qDF2jxANOG3HK+FclFJMHxLF5P4RrNh5jLe3HObPX6T9bLuRsSFUR91BdPJLTXbEapOKU0aHq0vut0557eXlBz0vM5pnTnvaejNPZa2D/HSY9arMZiV+RhK+6FAmT3duGhPDTWNiyC4qZ+fRUxSWVRPo48nwHsH0jvCH0p6wZ7FxpTr979Y58IFvLYOlXWOd8qyh/wxjtNATKe2eQP2srYvBLwIGz7FOeaJTkYQv7KZHqG/jnakCusKQXxgdsS572Jgxqr3Sv4SAKMcaLbLvNEAZV/nWSPgFGZDxLUz6E3i0cogL4RLaVYevlApVSn2nlMqwPDd6h00pVaeUSrY8VrbnmMJFjL0Has7AzubnGmhWTQVkrjE6WzlSj1P/LhAz1nq9bre9Cu5eMOoO65QnOp32/vU/BKzRWvcB1ljeN6ZCaz3c8ri2nccUriBqGMSOh+2vQ107J3zJ2mB8eNhw7Ps26z8DTqbAqSPtK6fiFCS/B0Pmgn+EdWITnU57E/5MYJnl9TJgVjvLcwqvvfYaCxcubPH2Tc12BY4745VDGPc/UJIN6V80v+2FpH4KpiCIu7T5bTta/RAP7W2ts3M51JQb4+YI0YT2JvxIrXWu5fUJILKJ7UxKqUSl1Fal1KymClNKLbBsl5ifn9/O0GwnJSWFoUNbXufa1GxX4PgzXtlV32lGx6Gti5vdtEnV5Ub9/cCZ4OFltdCsJqyXMa9ue6p16mqMb0Kx461381d0Ss3etFVKfQ90bWTVww3faK21UqqpbpSxWutjSqmewFqlVIrW+uD5G2mtXwdeBxg1atSFZ9T4+iGjdYM1dR0CVz3d7GZ79uzh+uuvZ+7cuURGRpKcnEx2djbvvvsur732Gtu2bePSSy/lzTeNgcCamu0Kzp3xqrnZrgC7zXhlF27uRueh1X+EnCSIHtn6Mg58bcz2NGSu9eOzlv4zYNMLxkihvqGt3z/tc+ObkLVaNIlOq9krfK31FK314EYenwMnlVJRAJbnvCbKOGZ5zgLWAw7UVKL19u7dy9ChQ0lJSaFnz55s2rSJhQsXMn/+fJ599lnS0tL46quvqKr6abLv82e7gp/PeNXcbFdg/xmvOtyIm40RH5sYK79Zez6GgG4Qe4l147Km/jOMJqPpX7Z+X63hx5eMDmV9rrR+bKJTaW+zzJXAbcDTlufPz9/A0nKnXGtdpZQKBy4Bnm3ncVt0JW4L2dnZ+Pn5YTKZKC4u5oEHHgCMjkXz588nKioK+Glmqnrnz3YF58541ZLZrhqWa88ZrzqUdwAk3Gq0QCl5EoJaMeF3WT5kfm9Mlu7WfJWZ3XQbAeF9jXr4hFaOlXRkM+Qmw9UvOFYLJOGQ2vsX8jRwhVIqA5hieY9SapRSaollmwFAolJqN7AOeFpr/fMulk6ivv4+NTWVhIQE3Cz/ZLt372bMmDEA5OTk0K1bt7OJurHZruDcGa9aOttVfXmOMONVhxm9wBj5cfvrrdtv13Iw17Q+iXY0pWDk7ZCzHU6mtm7fH18C3zAYNs8moYnOpV0JX2tdqLW+XGvdx1L1U2RZnqi1vsvy+ket9RCt9TDL85sXLtWx7dmzhyFDhpCSknJOIt6zZ8/ZG7m7d+8+56ZuY7NdwbkzXrVktitwvBmvOkRILPS/GpLeguoWDLMMYK6DxP8YLXO69LNpeFYxbJ7Rhj5pWfPb1ju+yxhxc8zd9pm9Szgd+Q7YSvVX+CkpKQwfPhyAyspKKioqCAkx+p01TP7Q9GxX8NOMVy2Z7Qocd8Yrmxt3rzG65O73W7b9gW+g5KgxQboz8A2FgbOMtvQVLZwpbN3fwBQsTTFFyzU1M4q9HzLjVeOcYcYrmzCbtX59stbP9de6qqz5bV+bpPULg7Wure6Y+KwhN8WYCWvd081ve3Sbse3G52wfl3AqyIxXjqszznhlE0rBlU9B6XHY9M8Lb5vxLRzfCZf+DtydqOqr62CjI9bWf0NlSdPbmetg1e+MsYFGt7wDoBCS8B1AZ5zxyiZixsLgX8DmFyEvvfFtaqvgmz9BSLxz3sic9BBUnYa1f2l6mx1LjMnYr3wKvP07Ljbh9CThC+dy5V+NppqfzDcGRTvf+r9BYSbMeM4xe9Y2J2qYcd9hx5LGZ8M6kQLfPgq9p8AgGQJZtI7TJXzd2NR4AnCRn01AJMxabDRf/PAWY+iEerveMXqsJtxqJERnNfkRCO0JH90KeT91uqMgA96da9zglQlORBs41Xj4JpOJwsJCwsLCzumMJIxkX1hY6Brt8/tOhWv/BSt/Da9eYlTz5KfDvpXGROjTn7N3hO1jCoIb34dlV8OSKTDiFmP5zuXgaYLbvjSGVhailZSjXhWOGjVKJyYmnrOspqaGnJycs52VxLlMJhPR0dGu00Y/awN8/7jRHt0nBEbNh4l/6DyTfxRnG/cj6idf7z/dqNIKirZvXMKhKaWStNajGl3nTAlfiEbVVoObR+cdWqCuBpSbYw8PIRzGhRK+U1XpCNEoZ7w52xrO1LRUOLROekkkhBDifJLwhRDCRThsHb5SKh9o50SfNhcOFNg7CDtx5XMH1z5/Vz53cPzzj9VaN9qMy2ETvjNQSiU2dXOks3PlcwfXPn9XPndw7vOXKh0hhHARkvCFEMJFSMJvn1ZOwdSpuPK5g2ufvyufOzjx+UsdvhBCuAi5whdCCBchCV8IIVyEJPw2UErNVUqlKqXMSqlR561bpJTKVErtV0pdaa8YO4JS6gml1DGlVLLlMd3eMdmaUmqa5XebqZR6yN7xdDSl1GGlVIrl992pB7tSSi1VSuUppfY2WBaqlPpOKZVheQ6xZ4ytJQm/bfYCc4CNDRcqpQYCNwKDgGnAv5VSnX3Eqxe01sMtj1X2DsaWLL/LV4CrgIHAPMvv3NVcZvl9O2Vb9FZ4C+P/uKGHgDVa6z7AGst7pyEJvw201vu01vsbWTUT+EBrXaW1PgRkAqM7NjphQ6OBTK11lta6GvgA43cuOiGt9Uag6LzFM4FlltfLgFkdGVN7ScK3ru5AdoP3OZZlndl9Sqk9lq+/TvX1tg1c8fd7Pg18q5RKUkotsHcwdhCptc61vD4BRNozmNaS4ZGboJT6HujayKqHtdafd3Q89nKhnwOwGPg/jCTwf8A/gDs7LjphB+O11seUUhHAd0qpdMuVsMvRWmullFO1a5eE3wStdVsmRT0G9GjwPtqyzGm19OeglHoD+NLG4dhbp/v9tpbW+pjlOU8p9SlGNZcrJfyTSqkorXWuUioKyLN3QK0hVTrWtRK4USnlrZSKB/oA2+0ck81Y/uDrzca4md2Z7QD6KKXilVJeGDfoV9o5pg6jlPJTSgXUvwam0vl/5+dbCdxmeX0b4FTf9uUKvw2UUrOBl4AuwFdKqWSt9ZVa61Sl1EdAGlAL3Ku1rrNnrDb2rFJqOEaVzmFgoV2jsTGtda1S6j7gG8AdWKq1TrVzWB0pEvhUKQVG7nhPa73aviHZjlLqfWASEK6UygEeB54GPlJKzccYvv16+0XYejK0ghBCuAip0hFCCBchCV8IIVyEJHwhhHARkvCFEMJFSMIXQggXIQlfCCFchCR8IYRwEZLwhRDCRUjCF0IIFyEJX4gLUEpFKKUKlFLVzU12opR6USmllVIvd1R8QrSGJHwhLuyfQBjwvNY6rZltd1ieJ9g0IiHaSMbSEaIJSqkxwFbgFBCvtS5pZvuLgc1ANeDbyQfOE05IrvCFaNpjlud/N5fsLfItz14Y3wqEcCiS8IVohFIqFmOycoA3zlsXo5Qar5Q6P6k3/H/q7JPXCyckCV+Ixk0HFJCqtT5y3rq/AT8A8ect72Z5rsXJZkISrkESvhCNG2l53tHIukstz/vPWz7C8pws9ffCEUnCF6JxsZbn7IYLlVIJGPPaVmitS8/bZ7rl+WsbxyZEm0jCF6Jxnpbn8/9H7rE8n1NHr5TqD0wGzMDbtg1NiLaRhC9E4+rr7S9XSrkBKKUuA+4APgW8lFKjLcv9gaUYdf5LtdaZdohXiGZJO3whGqGUmgyssbxNwqjamQ6sB/4EJAIngQ3AeIwbttuByVrrMx0drxAtIVf4QjRCa70WuBVIAwZj3MT9F3Cd1joJ+D1GB6uZQClGm/1JkuyFI5MrfCGEcBFyhS+EEC5CEr4QQrgISfhCCOEiJOELIYSLkIQvhBAuQhK+EEK4CEn4QgjhIiThCyGEi5CEL4QQLuL/AbZh0yoHOiPmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "U = 2.0\n",
    "D = 2.0\n",
    "beta = 10**6\n",
    "w = np.linspace(-6*D, 6*D, 2**12)\n",
    "gloc = semi_circle_hiltrans(w + 1e-3j,D)\n",
    "#gloc= 0.0\n",
    "gloc, sigma_loc = dmft_loop(gloc, w, U, beta, D, 100)\n",
    "fig,ax = plt.subplots(2,1)\n",
    "ax[0].plot(w, -gloc.imag, label=r'$\\pi A(\\omega)$')\n",
    "ax[0].legend(loc=0)\n",
    "ax[1].plot(w, sigma_loc.real, '-', label=r'$Re \\Sigma(\\omega)$')\n",
    "ax[1].plot(w, sigma_loc.imag, '-', label=r'$Im \\Sigma(\\omega)$')\n",
    "ax[1].set_xlabel(r'$\\omega$')\n",
    "ax[1].legend(loc=0)\n",
    "#plt.ylim([0, 0.1])\n",
    "#plt.xlim([-1, 1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The effect of the self-consistency condition\n",
    " \n",
    "1. The DMFT implementation is an iterative process which requires a certain number of loops before convergence is reached. But one needs a starting point, i.e. one needs to input an initial guess of the spectral function. Can you say which is the starting guess in the DMFT implementation considered above?  \n",
    "\n",
    "2. For $U < Uc_2$, run the code with *just one dmft iteration*, using the same values of $U$ that you have employed above. This one-loop run corresponds to the solution of the Single Impurity Anderson Model SIAM (or **Kondo Model**, see Marc's lecture notes) for a semi-circular conduction band. Compare the spectral function $A(\\omega)$ of the one-iteration loop with the fully converged one. Do you see the MIT?  \n",
    "\n",
    "3. Repeat the same excercise of point 2., this time inputing no bath (i.e. the intial guess is zero), and always one iteration loop. This means that you are essentially solving the isolated atom problem with on-site interaction $U$. Draw a conclusions about the effect of the self-conistency in solving the MIT problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Let's start now\n",
    "\n",
    "## Excercises\n",
    "\n",
    "We shall first work at low temperature, choose $\\beta=100$. The imaginary part of the Green's function $A(\\omega)= -\\Im G(\\omega)/\\pi$ can be directly connected to physical observables (e.g. spectral functions) or transport (e.g. specific heat).\n",
    "\n",
    "## The Metal-to-Insulator Transition\n",
    "\n",
    "1. Plot first $A(\\omega)$ in the case $U=0$. This is the semicircular density of states used in the DMFT code, i.e. $-\\Im G(\\omega)/\\pi= D(\\omega)$. The system is in this case \"halffilled\", i.e. there is one elctron per atomic orbital. Can you indicate then where the occupied states are on the plot of the density of states $A(\\omega)$? Why the system is a metal according to the band theory of solids?\n",
    "\n",
    "2. Run the code for several values of $U$ (the half-bandwidth is set $D=1$) to check out the metal-insulator transition. Display $A(\\omega)$ for some rapresentative values of $U$, in the metallic and insulating sides.\n",
    "\n",
    "3. Approach the MIT at $Uc_2$ (metal to insulator) from below ($U< Uc_2$) and describe how the different contributions to the $A(\\omega)$ (quasiparticle peak and Hubbard bands) evolve as you get close to the critical value. Why can we state that for $U> Uc_2$ the system is finally insulating? Give your extimation of  $Uc_2$. \n",
    "\n",
    "4. At the MIT point, the quaiparticle peak has completely disappered. But where is the spectral weight of the peak gone? Explain why we expect that the total spectral weight [i.e. the area under the curve $A(\\omega)$] is conserved.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selfenergy\n",
    "\n",
    "1. Plot now the real part and the imaginary part of the self-energy $\\Sigma(\\omega)$, first for $U< Uc_2$, then for $U \\ge Uc_2$. Use the same values of $U$ that you have used above.\n",
    "2. We are going now to concentrate on the low energy part $\\omega\\to 0$, where a quasiparticle peak is observable in the metallic phase, even for values of $U$ very close to the MIT. The Fermi liquid theory of metals states that the selfenergy is a regular function of $\\omega$, which can be Taylor expanded: $$\n",
    "\\Re\\Sigma \\simeq \\mu_0 + \\alpha \\omega+ O(\\omega^2) \\quad \\quad\n",
    "\\Im\\Sigma \\sim \\omega^2+ O(\\omega^3)\n",
    "$$\n",
    "For which values of $U$ can we state that the system is in a Fermi liquid state?\n",
    "3. We shall now see how the Fermi liquid theory may be useful to understand the beahvior of spectral functions. The local Green's function (entering the DMFT self-consistency condition) is the Hilbert transform of the lattice Green's function $G(\\omega,\\varepsilon)$ \n",
    "$$\n",
    "G_{loc}(\\omega)= \\,\\int D(\\varepsilon) \\, G(\\omega, \\varepsilon)\\, d\\varepsilon=\n",
    "\\,\\int D(\\varepsilon)/\\left[ \\omega-\\varepsilon- \\Sigma(\\omega) \\right] \\, d\\varepsilon $$\n",
    "where $D(\\varepsilon)= -(1/\\pi) \\Im G_{loc}(\\varepsilon,U=0)$ and in our case (infinite dimension) $\\Sigma(\\omega)$ is independent of $\\varepsilon$. Show that the lattice Green's function can be written within the Fermi liquid theory ($\\omega\\to0$) as a almost free-particle Green's function $$\n",
    "G(\\omega, \\xi)\\sim \\frac{Z}{\\omega-\\xi} $$\n",
    "where the main difference is the factor $Z$, known as *quasi-particle renormalization factor* .\n",
    "    * Evaluate $Z$ for different values of $U$ (use the real part of the self-energy) and plot it as a function of $U$. Explain why $Z$ is useful to describe the MIT. *Hint: in order to extract a part of the Green's function array (let's call it e.g. $gx$) between frequencies, e.g. $[-2,2]$, define a logical array of frequencies $w_r= (w< 2 )*(w> -2)$, then the desired reduced Green's function array is simply given by $gx[w_r]$.*\n",
    "    * Calculate (for a couple of $U$ values used above) the area under the quasiparticle peaks and plot them as a function of $U$. Relate this area with the quasiparticle residue $Z$ (you should plot them as a function of $U$). \n",
    "    * Determine the exact relation between the quasiparticle peak and $Z$ also analytically. Explain finally why it is said that *a Mott transition is a breaking of the Fermi Liquid Theory*.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First order transition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MIT is a first order transiton, i.e. there is a region of coexistence between the insulator and the metallic phase.\n",
    "1. Show that there is an interval of $U$ where DMFT self-consistency provides two solutions. Determine the $Uc_1$, i.e. the critical value of $U$ where the insulating state changes into a metal by reducing $U$. *(Hint: you have to use an appropriate starting guess to obtain the insulating solution, start your investigations with values of $U$ in the insulating side).*\n",
    "2. Is $Z$ the right parameter to describe the insulator-to-metal transition too? Argument your answer. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $T-U \\,$ phase diagram and the $T$-driven transition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Run now the code for higher temperatures $T$, i.e. lower inverse temperatures $\\beta= 75, 50, 25$ and determine for each temperature $T$ the $Uc_1$ and the $Uc_2$. Sketch then the coexistence region in the $T-U$ plane.\n",
    "2. Observe the shape of the coexistence region in the $T-U$ space. Fix the interaction $U$ inside this co-existence region (e.g. $U=2.9$) and systematically increase the temperature. Compare the different density of states. Can you now explain the phase-diagram of VO$_2$, in particular the unusual not-metallic behavior at high temperature? \n",
    "3. Can you give an interpretation of this physical behaviour with temperature in terms of the associated SIAM model, i.e. a spin impurity fluctuating with temperature in a bath of free electrons?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}