Measuring static observables / impurity density matrix

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,

\[\hat\rho_\mathrm{imp} = \mathrm{Tr}_\mathrm{bath}[e^{-\beta\hat H}/Z].\]

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):

\[\langle\hat O\rangle = \mathrm{Tr}_\mathrm{imp}[\hat\rho_\mathrm{imp} \hat O].\]

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.

Measuring \(\hat{\rho}_\mathrm{imp}\) in cthyb

Let us consider a single-orbital Anderson impurity problem in a weak external magnetic field (script: static.py).

At first, we import all necessary modules, define some input parameters and construct a ctyhb Solver object in the usual way:

from h5 import HDFArchive
from triqs.gf.descriptors import Fourier
from triqs.gf import Gf, MeshImFreq, iOmega_n, inverse, GfImTime, BlockGf, Wilson
from triqs.operators import c, c_dag, n

# plotting interface
from triqs.plot.mpl_interface import plt,oplot

# redirect ctyhb c++ output to notebook
from triqs.utility.redirect import start_redirect
start_redirect()

from triqs_cthyb.solver import Solver

Warning: could not identify MPI environment!

Starting serial run at: 2023-08-14 16:13:15.854417

# Parameters
D = 1.0         # Half-bandwidth of the bath
V = 0.2         # Hybridisation amplitude
U = 4.0         # Coulomb interaction
e_f = -U/2      # Local energy level
h = 0.01        # External field
beta = 50       # Inverse temperature

# Construct the impurity solver with the inverse temperature
# and the structure of the Green's functions
S = Solver(beta = beta, gf_struct = [ ('up',1), ('down',1) ])

# Initialize the non-interacting Green's function S.G0_iw
# External magnetic field introduces Zeeman energy splitting between
# different spin components
S.G0_iw['up']   << inverse(iOmega_n - e_f + h/2 - V**2 * Wilson(D))
S.G0_iw['down'] << inverse(iOmega_n - e_f - h/2 - V**2 * Wilson(D))

# setting up a local interaction Hamiltonian
h_int = U * n('up',0) * n('down',0)

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 for more information.

Now we run the solver:

# Run the solver
S.solve(h_int = h_int,                  # interaction Hamiltonian
        n_cycles  = 500000,             # Number of QMC cycles
        length_cycle = 200,             # Length of one cycle
        n_warmup_cycles = 10000,        # Warmup cycles
        measure_density_matrix = True,  # Measure the reduced density matrix
        use_norm_as_weight = True)      # Required to measure the density matrix

Insert error : recovering ...

The results of the density matrix accumulation are accessible via the density_matrix attribute of the solver:

# Extract accumulated density matrix
rho = S.density_matrix

Compute expactation values of static observables

Now we can use the function trace_rho_op() from triqs.atom_diag (see doc) to compute expectation values \(\langle\hat O\rangle = \mathrm{Tr}_\mathrm{imp}[\hat\rho_\mathrm{imp} \hat O]\) of any time independent operator.

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:

# Object containing eigensystem of the local Hamiltonian
h_loc_diag = S.h_loc_diagonalization

This allows now to compute for example the expectation values of the impurity occupations, or the double occupancy:

from triqs.atom_diag import trace_rho_op

# Evaluate impurity occupations
print("<N_up> =", trace_rho_op(rho, n('up',0), h_loc_diag))
print("<N_down> = ", trace_rho_op(rho, n('down',0), h_loc_diag))

# Evaluate double occupancy
print("<N_up*N_down> =", trace_rho_op(rho, n('up',0)*n('down',0), h_loc_diag))

<N_up> = 0.6200042407162714
<N_down> =  0.3801439064667559
<N_up*N_down> = 0.003665497313727422

Measuring moments of Green’s function and self-energy

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.

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,

# accessing the high frequency moments stored in the Solver class
S.Sigma_moments, S.G_moments

({'up': array([[[1.52057563+0.j]],

         [[3.77015227+0.j]]]),
  'down': array([[[2.48001696+0.j]],

         [[3.76958372+0.j]]])},
 {'up': array([[[ 0.        +0.j]],

         [[ 1.        +0.j]],

         [[-0.48442437+0.j]]]),
  'down': array([[[0.        +0.j]],

         [[1.        +0.j]],

         [[0.48501696+0.j]]])})

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.

Evaluation of the impurity interaction energy

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}\):

Eint = trace_rho_op(rho, h_int, h_loc_diag)
print("E_int=<H_int>=", '{:5.4f}'.format(Eint))

E_int=<H_int>= 0.0147

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.

In comparison to the widely used Galitski-Migdal Formula: \(\frac{1}{2} Tr [ \Sigma G]\) (Ref: Galitski & Migdal paper), 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.

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 for more information. We first used the pyED exact diagonalization benchmark 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)\).

The results are depicted below:

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.