{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Single-orbital Hubbard model\n", "============================\n", "\n", "In this notebook you will reproduce the Bethe lattice DMFT that we did earlier with IPT,\n", "but you will use the CTHYB solver to find the solution of the impurity problem. We will\n", "consider the problem at half-filling again. \n", "\n", "In general it is a good idea to develop the script in the notebook, because it is very convenient to find bugs and to quickly come to a working code. In the beginning you should use a small number of Monte Carlo iterations (say 1000) so that the code runs quickly. Your first main goal is to have a functional script. However, once the script is done, we recommend that you do longer runs (production runs) from a shell. It will be easier for you to see the progress of the Monte Carlo solver. Think about saving the relevant data to an archive and then go back to the notebook when it comes to analyzing and plotting the results. This is usually how things are done: elaboration of the scripts and analysis of the data from the notebook, production from a shell.\n", "\n", "In order to run your script from a shell, open a terminal and go in the tutorial directory.\n", "This is where you should edit your production script. Let's call it `run_dmft.py`. Use your favourite editor (e.g. `vi` or `gedit`) to create the script `run_dmft.py`.\n", "\n", "When the script is written save it and run it. You can:\n", "\n", "- run it directly from the shell to see the Monte Carlo progress:\n", "\n", "`>>> triqs run_dmft.py`\n", "\n", "- run the following command in a notebook cell.\n", "\n", "`%run run_dmft.py`\n", "\n", "That's it! When the run is done and data has been saved into an archive, you can go back to the notebook and read the archive in order to analyze or plot the results." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 1\n", "----------\n", "\n", "Write a DMFT loop, like you did earlier but using the CTHYB solver. \n", "\n", "*Hint*: It is useful to symmetrize the Green's function (make the `up` and `down` components the same) to avoid some artificial polarization of the system close to the Mott transition. You might want to enforce the `up`-`down` symmetry on `S.G` just before the self-consistency condition. In order\n", "to have reasonable data you should have at least 10000 cycles." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 1\n", "----------\n", "\n", "The solution of this exercise is in the script `scripts/one_band.py` in the tutorial directory. We have loaded this script into the cell below. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "# %load scripts/one_band.py\n", "from triqs.gf import *\n", "from triqs.operators import *\n", "from h5 import *\n", "from triqs_cthyb import Solver\n", "import numpy as np\n", "\n", "# Parameters of the model\n", "t = 1.0\n", "beta = 10\n", "n_loops = 10\n", "\n", "# Construct the impurity solver\n", "S = Solver(beta = beta, gf_struct = [('up',1), ('down',1)] )\n", "\n", "# I run for several values of U\n", "for U in np.arange(1.0, 9.0):\n", " print('U =', U)\n", "\n", " # This is a first guess for G\n", " S.G_iw << SemiCircular(2*t)\n", "\n", " # DMFT loop with self-consistency\n", " for i in range(n_loops):\n", " \n", " print(\"\\n\\nIteration = %i / %i\" % (i+1, n_loops))\n", " \n", " # Symmetrize the Green's function and use self-consistency\n", " g = 0.5 * ( S.G_iw['up'] + S.G_iw['down'] )\n", " for name, g0 in S.G0_iw:\n", " g0 << inverse( iOmega_n + U/2.0 - t**2 * g )\n", "\n", " # Solve the impurity problem\n", " S.solve(h_int = U * n('up',0) * n('down',0), # Local Hamiltonian \n", " n_cycles = 10000, # Number of QMC cycles\n", " n_warmup_cycles = 5000, # Warmup cycles\n", " )\n", " \n", " # Save iteration in archive\n", " with HDFArchive(\"data/one_band/half-U%.2f.h5\"%U) as A:\n", " A['G-%i'%i] = S.G_iw\n", " A['Sigma-%i'%i] = S.Sigma_iw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The script saves the Green's functions and self-energies in archives in the subdirectory called `results`. They will be used for the analysis below. We study several values of $U$ to see how the Green's functions will change. Note that these runs takes a bit of time!\n", "\n", "You can\n", "\n", "- run this directly from the shell to see the Monte Carlo progress or \n", "- run the command below." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "%run scripts/one_band.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 2\n", "----------\n", "\n", "Here, you will learn to analyze the output of the solver. As discussed, the Monte Carlo algorithm provide results on the\n", "Matsubara axis. This makes the analysis of the results slightly more delicate than if we had them directly on the real\n", "axis. When we used the IPT solver, we could see the Mott transition as the appearance of a gap in the spectral function.\n", "After the Monte Carlo run, we do not have the spectral function, so we will have to use some other criteria to decide, e.g.,\n", "if the system is metallic or insulating.\n", "\n", "Plot the Green's function at the end of the DMFT loops for several values of $U$ (say between 2 and 8). Focus on the extrapolation of the imaginary part of the Green's function to zero frequency. How does it change with $U$? Is there\n", "a way to see the Mott transition just by inspecting the imaginary part of the Green's function?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 2\n", "------------\n", "\n", "Start by reading the data from the archive and plot it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from triqs.gf import *\n", "from h5 import HDFArchive\n", "%matplotlib inline\n", "from triqs.plot.mpl_interface import plt,oplot\n", "import matplotlib as mpl\n", "mpl.rcParams['figure.dpi']=100 \n", "\n", "A1 = HDFArchive(\"data/one_band/half-U2.00.h5\", 'r')\n", "A2 = HDFArchive(\"data/one_band/half-U6.00.h5\", 'r')\n", "\n", "# Plot the Green's function of the last iteration\n", "oplot(A1['G-9']['up'], 'x', label='U = 2.0')\n", "oplot(A2['G-9']['down'], 'o', label='U = 6.0')\n", "\n", "plt.xlim(0,10)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the behavior of the imaginary part is very different for the two values of $U$. When\n", "$U$ is small, the system is a metal and the imaginary part extrapolated to zero goes to a finite value.\n", "Instead, for large $U$, the system is a Mott insulator and the imaginary part goes to zero. The reason\n", "is that the extrapolation to zero is directly proportional to the density of states at the chemical\n", "potential. If the system is gapped, the density is zero; if the system is a metal, there is spectral\n", "weight and the density is finite. Therefore, even on the Matsubara axis, one has a way to decide if the\n", "system is metallic or not." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 3\n", "----------\n", "\n", "Do the same exercise as above, but analyze the self-energy. The noise usually gets bigger for larger frequencies, so just focus on the first few Matsubara frequencies. There the noise should not be too important. Again, by looking at the extrapolation to zero frequency of the imaginary part of the self-energy, can you tell where the Mott transition happens?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 3\n", "------------\n", "\n", "We now do the same for the self-energy." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "A1 = HDFArchive(\"data/one_band/half-U2.00.h5\", 'r')\n", "A2 = HDFArchive(\"data/one_band/half-U6.00.h5\", 'r')\n", "\n", "# Plot the self-energy of the last iteration\n", "oplot(A1['Sigma-9']['up'].imag, 'x', name='U = 2.0')\n", "oplot(A2['Sigma-9']['up'].imag, 'o', name='U = 6.0')\n", "\n", "plt.xlim(0,5)\n", "plt.ylim(-30,0)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 4\n", "----------\n", "\n", "A very useful quantity to measure the degree of correlation of a metal is the *quasiparticle weight* $Z$. It is defined as\n", "\n", "$$\n", "Z = \\left[ 1 - \\frac{\\partial \\text{Im} \\Sigma(i\\omega)}{\\partial \\omega} \\big|_{\\omega \\rightarrow 0} \\right]^{-1}\n", "$$\n", "\n", "For a non-interacting metal $Z=1$. As correlations appear, $Z$ gradually gets smaller. It reaches 0 at the Mott transition. Make a plot of $Z$ versus $U$ for the Bethe lattice Hubbard model. \n", "\n", "*Hint*: In order to have access to the values of $\\Sigma_\\uparrow(i\\omega_n)$, you can use `S.Sigma['up'](n)`. This will be useful to numerically compute the derivative required to compute $Z$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 4\n", "------------\n", "\n", "We estimate the derivative using the following approximation\n", "\n", "$$\n", "Z = \\lim_{\\omega_n \\rightarrow 0} \\Big( 1 - \\frac{d\\mathrm{Im} \\Sigma(i\\omega_n)}{di\\omega_n} \\Big)^{-1}\n", " \\sim \\Big( 1 - \\frac{\\mathrm{Im} \\Sigma(i\\omega_0)}{i\\omega_0} \\Big)^{-1}\n", "$$\n", "\n", "with $\\omega_0 = \\pi / \\beta$ being the first Matsubara frequency." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "beta = 10\n", "U_list = []\n", "Z_list = []\n", "\n", "for U in np.arange(1.0, 9.0):\n", " \n", " A = HDFArchive(\"data/one_band/half-U%.2f.h5\"%U, 'r')\n", " Sigma = A['Sigma-9']\n", " \n", " Z = 1 / (1 - (Sigma['up'](0)[0,0].imag * beta / np.pi))\n", " U_list.append(U)\n", " Z_list.append(Z)\n", " \n", "plt.plot(U_list, Z_list, '-o')\n", "plt.xlabel('U')\n", "plt.ylabel('Z')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 5\n", "----------\n", "\n", "Go back to your IPT code and try to modify it to extract the $Z$ versus $U$ curve. Compare this to the result you found in Exercise 4. Is the critical $U$ for the Mott transition similar to the one you found using CTHYB?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 5\n", "------------\n", "\n", "This is just the same script as we had earlier. We just add a couple of lines to extract $Z$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from triqs.gf import *\n", "from triqs.plot.mpl_interface import *\n", "\n", "class IPTSolver:\n", " def __init__(self, beta):\n", " self.beta = beta\n", "\n", " # Matsubara frequency Green's functions\n", " iw_mesh = MeshImFreq(beta=beta, S='Fermion', n_iw=1001)\n", " self.G_iw = Gf(mesh=iw_mesh, target_shape=[1,1])\n", " self.G0_iw = self.G_iw.copy() # self.G0 will be set by the user after initialization\n", " self.Sigma_iw = self.G_iw.copy()\n", " \n", " # Imaginary time\n", " tau_mesh = MeshImTime(beta=beta, S='Fermion', n_tau=10001)\n", " self.G0_tau = Gf(mesh=tau_mesh, target_shape=[1,1])\n", " self.Sigma_tau = self.G0_tau.copy()\n", " \n", " def solve(self, U):\n", " self.G0_tau << Fourier(self.G0_iw)\n", " self.Sigma_tau << (U**2) * self.G0_tau * self.G0_tau * self.G0_tau\n", " self.Sigma_iw << Fourier(self.Sigma_tau)\n", " \n", " # Dyson\n", " self.G_iw << inverse(inverse(self.G0_iw) - self.Sigma_iw)\n", "\n", "t = 1.0\n", "beta = 10\n", "n_loops = 30\n", "\n", "S = IPTSolver(beta = beta)\n", "\n", "U_list2 = []\n", "Z_list2 = []\n", "\n", "for U in np.arange(0.0, 9.0):\n", " \n", " S.G_iw << SemiCircular(2*t)\n", " for i in range(n_loops):\n", "\n", " S.G0_iw << inverse( iOmega_n - t**2 * S.G_iw )\n", " S.solve(U = U)\n", " \n", "\n", " Z = 1 / (1 - (S.Sigma_iw(0)[0,0].imag * beta / np.pi))\n", " U_list2.append(U)\n", " Z_list2.append(Z)\n", " \n", " print(Z)\n", " \n", "plt.plot(U_list2, Z_list2, '-o', label='IPT')\n", "plt.plot(U_list, Z_list, '-x', label='CTQMC')\n", "plt.xlabel('U')\n", "plt.ylabel('Z')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 6\n", "----------\n", "\n", "Try to analytically continue the Green's function on the real axis using the Pade approximation. What can\n", "you say about the result?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solution 6\n", "----------" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from math import pi\n", "\n", "with HDFArchive(\"data/one_band/half-U4.00.h5\", 'r') as A:\n", " g = A['G-9']['up']\n", "\n", "g_real = GfReFreq(indices=[0], window=[-5,5])\n", "g_real.set_from_pade(g)\n", "oplot(-g_real.imag/pi)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is completely wrong. This is because of the noise in the Monte Carlo data. One would have to make much longer runs in order to reduce the error bars. The Pade approximation can be used only on very accurate data. When the noise is still quite large, one has to use different analytical continuation methods, like MaxEnt, which produces the following spectral function:\n", "\n", "\n", "\n", "Regardless of which package you use for MaxEnt, it is very important to remember that there are some important knobs with which one can play in MaxEnt that can substantially change the results, and so one must be very careful in its use!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " Exercise 7\n", "----------\n", "\n", "In the last exercise we will analyze the impurity multiplets, the many-body states and their occupancies, using the CTHYB solver. To do this, we instruct the solver to measure the many-body density matrix of the impurity during sampling. This can be done at almost no additional numerical cost in CTHYB by passing the options `measure_density_matrix = True` and `use_norm_as_weight = True` to the `solve` call of the solver.\n", "\n", "Write a small script that measures the impurity density matrix for two U-values (metallic and insulating) and analyzes/compares the two solutions by saving the solver object after convergence is reached for both values. Then you can access the density matrix as `S.density_matrix`. We provide a function to analyze the multiplets within CTHYB: `triqs_cthyb.multiplet_tools.multiplet_analysis`. Have a look at the [documentation](https://triqs.github.io/cthyb/latest/guide/multiplet_analysis_notebook.html) to better understand how to use it. Use the code from exercise one to do this by extending it. \n", "\n", "The function `multiplet_analysis` will also characterize each impurity many-body state with its quantum numbers. What is the characteristic feature in the multiplet structure of the Mott insulating state?\n", "\n", "Tip: By default, printing a panda's data frame shows all stored multiplets sorted by occupation probability." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true, "tags": [] }, "outputs": [], "source": [ "# Parameters of the model\n", "t = 1.0\n", "beta = 10\n", "n_loops = 10\n", "\n", "# Construct the impurity solver\n", "S = Solver(beta = beta, gf_struct = [('up',1), ('down',1)] )\n", "\n", "# run metallic U\n", "U = 2.0\n", "\n", "# This is a first guess for G\n", "S.G_iw << SemiCircular(2*t)\n", "\n", "# DMFT loop with self-consistency\n", "for i in range(n_loops):\n", " print(\"\\n\\nIteration = %i / %i\" % (i+1, n_loops))\n", "\n", " # Symmetrize the Green's function and use self-consistency\n", " g = 0.5 * ( S.G_iw['up'] + S.G_iw['down'] )\n", " for name, g0 in S.G0_iw:\n", " g0 << inverse( iOmega_n + U/2.0 - t**2 * g )\n", "\n", " # Solve the impurity problem\n", " S.solve(h_int = U * n('up',0) * n('down',0), # Local Hamiltonian \n", " n_cycles = 10000, # Number of QMC cycles\n", " n_warmup_cycles = 5000, # Warmup cycles\n", " measure_density_matrix = True,\n", " use_norm_as_weight = True\n", " )\n", "\n", "solver_U_met = S\n", "\n", "################################\n", "\n", "U = 8.0\n", "S = Solver(beta = beta, gf_struct = [('up',1), ('down',1)] )\n", "# This is a first guess for G\n", "S.G_iw << SemiCircular(2*t)\n", "\n", "# DMFT loop with self-consistency\n", "for i in range(n_loops):\n", " print(\"\\n\\nIteration = %i / %i\" % (i+1, n_loops))\n", "\n", " # Symmetrize the Green's function and use self-consistency\n", " g = 0.5 * ( S.G_iw['up'] + S.G_iw['down'] )\n", " for name, g0 in S.G0_iw:\n", " g0 << inverse( iOmega_n + U/2.0 - t**2 * g )\n", "\n", " # Solve the impurity problem\n", " S.solve(h_int = U * n('up',0) * n('down',0), # Local Hamiltonian \n", " n_cycles = 10000, # Number of QMC cycles\n", " n_warmup_cycles = 5000, # Warmup cycles\n", " measure_density_matrix = True,\n", " use_norm_as_weight = True\n", " )\n", "\n", "solver_U_ins = S\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# import and use multiplet analysis function to get panda data frame of the impurity multiplets\n", "from triqs_cthyb.multiplet_tools import multiplet_analysis" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# using now the multiplet analysis function to obtain a panda data frame of the occupied multiplets\n", "res_met = multiplet_analysis(rho = solver_U_met.density_matrix,\n", " h_loc_diag = solver_U_met.h_loc_diagonalization,\n", " # one orbital\n", " n_orb = 1,\n", " # two spin channels\n", " spin_names = ['up','down'])\n", "res_met" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "res_ins = multiplet_analysis(rho = solver_U_ins.density_matrix,\n", " h_loc_diag = solver_U_ins.h_loc_diagonalization,\n", " n_orb = 1,\n", " spin_names = ['up','down'])\n", "res_ins" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that in the insulating state only the high spin states with $m_s= \\pm 0.5$ are occupied, while in the metallic state also the low spin states with $m_s=0.0$ are occupied. We can also see that their relative energy eigenvalues within the impurity change when U is varied." ] } ], "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.10.8" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 4 }