{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "import matplotlib\n", "matplotlib.use(\"Pdf\")\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "\n", "# set matplotlib parameters\n", "params = {'backend': 'ps',\n", " 'axes.labelsize': 13,\n", " 'font.size': 13,\n", " 'legend.fontsize': 13,\n", " 'xtick.labelsize': 13,\n", " 'ytick.labelsize': 13,\n", " 'text.usetex': False,\n", " 'text.latex.preamble': \"\\\\usepackage{mathpazo}, \\\\usepackage{amsmath}\",\n", " 'font.family': 'sans-serif',\n", " 'font.sans-serif': ['Computer Modern Sans serif']}\n", "plt.rcParams.update(params)\n", "\n", "import warnings \n", "warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n", "\n", "# numpy\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Measuring static observables / impurity density matrix\n", "\n", "In addition to the interacting Green's functions one is often interested in such observables, as orbital occupation numbers and double occupancy of the impurity. Expectation values of these time-independent operators are readily expressed in terms of the reduced many-body density matrix of the system, \n", "\n", "$$ \\hat\\rho_\\mathrm{imp} = \\mathrm{Tr}_\\mathrm{bath}[e^{-\\beta\\hat H}/Z]. $$\n", "\n", "Here, $e^{-\\beta\\hat H}/Z$ is the density matrix of the full system (impurity + bath) described by the Hamiltonian $\\hat H$, and the Tr operator runs only over the bath degrees. Knowledge of $\\hat{\\rho}_\\mathrm{imp}$ allows to compute the thermodynamic expectation value of any time independent operator $\\hat{O}$ acting on the impurity degrees of freedom ([wikipedia: density matrix](https://en.wikipedia.org/wiki/Density_matrix#Measurement)):\n", "\n", "$$ \\langle\\hat O\\rangle = \\mathrm{Tr}_\\mathrm{imp}[\\hat\\rho_\\mathrm{imp} \\hat O]. $$\n", "\n", "The ctyhb solver allows to measure the impurity density matrix in the many-body basis, with only little additional calculational cost, and thus gives access to the multiplet structure of the impurity problem.\n", "\n", "## Measuring $\\hat{\\rho}_\\mathrm{imp}$ in cthyb\n", "\n", "Let us consider a single-orbital Anderson impurity problem in a weak external magnetic field ([script: static.py](static.py)).\n", "\n", "At first, we import all necessary modules, define some input parameters and construct a ctyhb ``Solver`` object in the usual way:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Warning: could not identify MPI environment!\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Starting serial run at: 2023-08-14 16:13:15.854417\n" ] } ], "source": [ "from h5 import HDFArchive\n", "from triqs.gf.descriptors import Fourier\n", "from triqs.gf import Gf, MeshImFreq, iOmega_n, inverse, GfImTime, BlockGf, Wilson\n", "from triqs.operators import c, c_dag, n\n", "\n", "# plotting interface\n", "from triqs.plot.mpl_interface import plt,oplot\n", "\n", "# redirect ctyhb c++ output to notebook\n", "from triqs.utility.redirect import start_redirect\n", "start_redirect()\n", "\n", "from triqs_cthyb.solver import Solver" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Parameters\n", "D = 1.0 # Half-bandwidth of the bath\n", "V = 0.2 # Hybridisation amplitude\n", "U = 4.0 # Coulomb interaction\n", "e_f = -U/2 # Local energy level\n", "h = 0.01 # External field\n", "beta = 50 # Inverse temperature\n", "\n", "# Construct the impurity solver with the inverse temperature\n", "# and the structure of the Green's functions\n", "S = Solver(beta = beta, gf_struct = [ ('up',1), ('down',1) ])\n", "\n", "# Initialize the non-interacting Green's function S.G0_iw\n", "# External magnetic field introduces Zeeman energy splitting between\n", "# different spin components\n", "S.G0_iw['up'] << inverse(iOmega_n - e_f + h/2 - V**2 * Wilson(D))\n", "S.G0_iw['down'] << inverse(iOmega_n - e_f - h/2 - V**2 * Wilson(D))\n", "\n", "# setting up a local interaction Hamiltonian\n", "h_int = U * n('up',0) * n('down',0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we instruct the ``solve`` function to accumulate the density matrix by passing the solver parameters ``measure_density_matrix = True`` and ``use_norm_as_weight = True``. The latter parameter tells the solver to employ a reweighting scheme, using a spherical norm instead of the trace to calculate the atomic weight. This reweighting allows to take into account important contributions to $\\hat\\rho_\\mathrm{imp}$, which otherwise would be missed. See also [cthyb solver parameters](https://triqs.github.io/cthyb/latest/reference/solve_parameters.html) for more information.\n", "\n", "Now we run the solver:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Insert error : recovering ... \n" ] } ], "source": [ "# Run the solver\n", "S.solve(h_int = h_int, # interaction Hamiltonian\n", " n_cycles = 500000, # Number of QMC cycles\n", " length_cycle = 200, # Length of one cycle\n", " n_warmup_cycles = 10000, # Warmup cycles\n", " measure_density_matrix = True, # Measure the reduced density matrix\n", " use_norm_as_weight = True) # Required to measure the density matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results of the density matrix accumulation are accessible via the `density_matrix` attribute of the solver:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Extract accumulated density matrix\n", "rho = S.density_matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute expactation values of static observables\n", " \n", "Now we can use the function `trace_rho_op()` from `triqs.atom_diag` ([see doc](https://triqs.github.io/triqs/latest/documentation/cpp_api/triqs/atom_diag/trace_rho_op.html)) to compute expectation values $\\langle\\hat O\\rangle = \\mathrm{Tr}_\\mathrm{imp}[\\hat\\rho_\\mathrm{imp} \\hat O]$ of any time independent operator.\n", "\n", "To do so, information about the structure of the local Hilbert space is also needed. This information is stored as a special object in `h_loc_diagonalization`:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Object containing eigensystem of the local Hamiltonian\n", "h_loc_diag = S.h_loc_diagonalization" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This allows now to compute for example the expectation values of the impurity occupations, or the double occupancy:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " = 0.6200042407162714\n", " = 0.3801439064667559\n", " = 0.003665497313727422\n" ] } ], "source": [ "from triqs.atom_diag import trace_rho_op\n", "\n", "# Evaluate impurity occupations\n", "print(\" =\", trace_rho_op(rho, n('up',0), h_loc_diag))\n", "print(\" = \", trace_rho_op(rho, n('down',0), h_loc_diag))\n", "\n", "# Evaluate double occupancy\n", "print(\" =\", trace_rho_op(rho, n('up',0)*n('down',0), h_loc_diag))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Measuring moments of Green's function and self-energy\n", "The high-frequency moments of the Green's function and self-energy contain crucial high-frequency information that can be used to (1) improve the Fourier transform of the Green's function from imaginary time to Matsubara frequencies and (2) improve the tail fit of the self-energy. When a user turns on the measurement of the density matrix (``measure_density_matrix = True``), these high-frequency moments are automatically calculated and used. \n", "\n", "The moments are stored as members of the ``Solver`` class. ``Solver.Sigma_moments`` are the calculated self-energy moments and ``Solver.G_moments`` are the calculated Green's function moments. For example," ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "({'up': array([[[1.52057563+0.j]],\n", " \n", " [[3.77015227+0.j]]]),\n", " 'down': array([[[2.48001696+0.j]],\n", " \n", " [[3.76958372+0.j]]])},\n", " {'up': array([[[ 0. +0.j]],\n", " \n", " [[ 1. +0.j]],\n", " \n", " [[-0.48442437+0.j]]]),\n", " 'down': array([[[0. +0.j]],\n", " \n", " [[1. +0.j]],\n", " \n", " [[0.48501696+0.j]]])})" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# accessing the high frequency moments stored in the Solver class\n", "S.Sigma_moments, S.G_moments" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The moments are stored as Python dictonaries and have the same structure as ``gf_struct``. For more technical information and details, see the [high frequency moments notebook](high_freq_moments.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Evaluation of the impurity interaction energy\n", "\n", "The function `trace_rho_op()` can also be used to evaluate the interaction energy of the impurity by computing the expectation value of $\\hat{H}_\\mathrm{int}$:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "E_int== 0.0147\n" ] } ], "source": [ "Eint = trace_rho_op(rho, h_int, h_loc_diag)\n", "print(\"E_int==\", '{:5.4f}'.format(Eint))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This provides a powerful way to calculate the impurity interaction energy, which is for example needed for calculating the total energy in a combined DFT+DMFT calculation.\n", "\n", "In comparison to the widely used Galitski-Migdal Formula: $\\frac{1}{2} Tr [ \\Sigma G]$ ([Ref: Galitski \\& Migdal paper](http://www.jetp.ac.ru/cgi-bin/e/index/e/7/1/p96?a=list)), using $\\mathrm{Tr}_\\mathrm{imp}[ \\hat\\rho_\\mathrm{imp} \\hat{H}_\\mathrm{int}]$ has the advantage that the self-energy $\\Sigma$ is not needed for calculating the interaction energy, which can be very noisy, especially for large $i \\omega_n$. Especially, the influence of the high-frequency tail of the real part of the self-energy $Re \\Sigma (i \\omega_n \\rightarrow \\infty)$, the so-called Hartree shift, critically determines the value calculated in the Galitski-Migdal formula. Therefore, calculating the interaction energy with `trace_rho_op()` provides a significantly more stable procedure when using a QMC based solver.\n", "\n", "To test and benchmark this feature of the cthyb solver, we performed calculations for a Dimer coupled to two discrete bath states, including a full Kanamori interaction Hamiltonian. See [triqs benchmark Dimer](https://github.com/TRIQS/benchmarks/blob/master/Dimer/) for more information. We first used the [pyED exact diagonalization benchmark](https://github.com/TRIQS/benchmarks/blob/master/common/pyed) to compute the exact interaction energy of the impurity. Next, we compared these results with interaction energies obtained from `trace_rho_op()`, and the Galitski-Migdal formula using the cthyb solver. We performed either a tail-fit, or used sampling directly in the Legendre basis to obtain smooth high-frequency behavior of $\\Sigma (i \\omega_n)$. \n", "\n", "The results are depicted below:\n", "\n", "![](Eint_benchmark.svg)\n", "\n", "The black dashed horizontal line represents the interaction energy calculated with the exact diagonalization solver. The green crosses show the interaction energy computed with cthyb using `trace_rho_op()` for an increasing number of total QMC cycles (`n_cycles*mpi.size`). The error in the interaction energy is smaller than 1 meV compared to the ED result in this system, even for only $48 \\times 10^6$ QMC cycles. In comparison, the interaction energy obtained from the Galitski-Migdal formula using Legendre sampling (blue circles), or using a tail-fitting procedure (orange squares), show errors as large as 70 meV. Even for incredible large numbers of QMC cylces ($4 \\times 10^9$), the error is still as large as $\\approx 10$ meV. Even though, the high-frequency tail of $\\Sigma(i \\omega_n)$ looks very smooth for this large number of QMC cycles, $Re \\Sigma (i \\omega_n \\rightarrow \\infty)$ has not exactly the same value as the $\\Sigma$ calculated with the ED solver, resulting directly in an error in the interaction energy. Hence, we highly recommend to calculate the interaction energy via `trace_rho_op()`, which converges very fast and reliably with increasing number of QMC cycles." ] } ], "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.10.7" } }, "nbformat": 4, "nbformat_minor": 2 }