Matrix-valued Green's functions (Anderson impurity model) giving strange results in Solver

[bfield1]

This issue may be related to #72

Description

I attempted to formulate the Anderson impurity model from a matrix perspective. The Anderson impurity model has an impurity with energy level $\epsilon_f$, a substrate with isolated Green's function $G_{sub}(i\omega_n)$ (typically assumed to be flat), and coupling between the impurity and substrate $V$. The non-interacting Green's function in the basis of the impurity and the substrate is

$$ G_0(i\omega_n) = \begin{pmatrix} i\omega_n - \epsilon_f & -V \\ -V & G_{sub}^{-1}(i\omega_n) \end{pmatrix}^{-1} $$

In the tutorial https://triqs.github.io.cthyb/latest/guide/aim.html, we take just the impurity element of this matrix (after inversion) and apply CTHYB's Solver. But why not put the entire matrix into CTHYB (besides potentially being inefficient)? CTHYB can handle matrix-valued Green's functions, after all.

So I did that, using the example script below. However, the matrix-valued Green's function gave garbage results. I would expect the Green's functions (in the impurity matrix element) to be the same between the impurity-only and full-matrix solver.

Plots

The non-interacting Green's functions in the impurity element (being fed into CTHYB) agree between the two methods.

But the interacting Green's functions (returned by CTHYB) disagree strongly.

Here's a detailed plot of the interacting Green's function at all the distinct matrix elements.



The impurity element has the very strange behaviour of being effectively constant except at the tail (where it seems to agree with the impurity-only solution).
The off-diagonal element is zero. This is, as far as I can tell, due to issue #72

Example Script

Note: this script took about 6 minutes to execute on a 4-core machine. Fewer cores may produce noisier results.

import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

from triqs.gf import *
from triqs.operators import *
from triqs_cthyb import Solver
from h5 import HDFArchive
import triqs.utility.mpi as mpi
from triqs.plot.mpl_interface import oplot

# Solve single impurity
# Adapted from https://triqs.github.io.cthyb/latest/guide/aim.html

def aim1(D, V, U, e_f, beta):
    # Construct the impurity solver with the inverse temperature
    # and the structure of the Green's functions
    S = Solver(beta = beta, gf_struct = [ ('up',[0]), ('down',[0]) ])

    # Initialize the non-interacting Green's function S.G0_iw
    for name, g0 in S.G0_iw: g0 << inverse(iOmega_n - e_f - V**2 * Wilson(D))

    # Run the solver. The results will be in S.G_tau and S.G_iw
    S.solve(h_int = U * n('up',0) * n('down',0),     # Local Hamiltonian
            n_cycles  = 500000,                      # Number of QMC cycles
            length_cycle = 200,                      # Length of one cycle
            n_warmup_cycles = 10000)                 # Warmup cycles

    # Save the results in an HDF5 file (only on the master node)
    archive = f"aim_solution_D{D}V{V}U{U}.h5"
    if mpi.is_master_node():
        with HDFArchive(archive,'w') as Results:
            Results["G_tau"] = S.G_tau
            Results["G_iw"] = S.G_iw
            Results["G0_iw"] = S.G0_iw
            Results["Sigma_iw"] = S.Sigma_iw
    return S
# Two orbital case

def aim2(D, V, U, e_f, beta):
    # Construct the impurity solver with the inverse temperature
    # and the structure of the Green's functions
    S2 = Solver(beta = beta, gf_struct = [ ('up',[0,1]), ('down',[0,1]) ])

    # Initialize the non-interacting Green's function S.G0_iw
    for name, g0 in S2.G0_iw:
        g0[0,0] << iOmega_n - e_f
        g0[1,0] << -V
        g0[0,1] << -V
        g0[1,1] << inverse(Wilson(D))
    S2.G0_iw.invert()

    # Run the solver. The results will be in S.G_tau and S.G_iw
    S2.solve(h_int = U * n('up',0) * n('down',0),     # Local Hamiltonian
            n_cycles  = 500000,                      # Number of QMC cycles
            length_cycle = 200,                      # Length of one cycle
            n_warmup_cycles = 10000)                 # Warmup cycles

    # Save the results in an HDF5 file (only on the master node)
    archive = f"aim2_solution_D{D}V{V}U{U}.h5"
    if mpi.is_master_node():
        with HDFArchive(archive,'w') as Results:
            Results["G_tau"] = S2.G_tau
            Results["G_iw"] = S2.G_iw
            Results["G0_iw"] = S2.G0_iw
            Results["Sigma_iw"] = S2.Sigma_iw
    return S2

def plot_G0(G1, G2):
    oplot(G1['up'])
    oplot(G2['up'][0,0].real, ls=':')
    oplot(G2['up'][0,0].imag, ls=':')
    plt.xlim(0,10)
    plt.legend(['1*1 Real','1*1 Imag','2*2 Real','2*2 Imag'])
    plt.title('G0_iw for Anderson Impurity')
    plt.savefig('G0_iw_aimvsaim2.png')
    plt.show()
    plt.close()

def plot_results(G1, G2):
    oplot(G1['up'])
    oplot(G2['up'][0,0])
    plt.xlim(0,10)
    plt.legend(['1*1 Real','1*1 Imag','2*2 Real','2*2 Imag'])
    plt.title('G_iw for Anderson Impurity')
    plt.savefig('G_iw_aimvsaim2.png')
    plt.show()
    plt.close()

    oplot(G2['up'][0,0])
    plt.xlim(0,None)
    plt.title('G_iw[0,0] for 2*2 matrix form of Anderson Impurity')
    plt.savefig('G_iw00_aim2.png')
    plt.show()
    plt.close()

    oplot(G2['up'][0,1])
    plt.xlim(0,None)
    plt.title('G_iw[0,1] for 2*2 matrix form of Anderson Impurity')
    plt.savefig('G_iw01_aim2.png')
    plt.show()
    plt.close()

    oplot(G2['up'][1,1])
    plt.xlim(0,10)
    plt.title('G_iw[1,1] for 2*2 matrix form of Anderson Impurity')
    plt.savefig('G_iw11_aim2.png')
    plt.show()
    plt.close()

if __name__ == "__main__":
    # Parameters
    D, V, U = 1.0, 0.2, 4.0
    e_f, beta = -U/2.0, 20

    S = aim1(D, V, U, e_f, beta)
    S2 = aim2(D, V, U, e_f, beta)

    plot_G0(S.G0_iw, S2.G0_iw)
    plot_results(S.G_iw, S2.G_iw)

Versions

TRIQS library version: 3.0.1
TRIQS git hash: fea808e79cd2e66627ddd8eed535317bc1561466
triqs_cthyb version: 3.0.0
triqs_cthyb git hash: ba23ff9

I used the Docker distribution of TRIQS, running through Singularity. I run on Linux on a supercomputing cluster.
Python 3.8.5.
OpenMPI 4.0.3.
The plots above were made with 4 CPUs. Running time was just over 6 minutes.

[Wentzell]

Dear @bfield1,

Thank you for pointing out this issue.

What you are seeing is indeed related to #72. By initializing your Weiss-field as

for name, g0 in S2.G0_iw:
        g0[0,0] << iOmega_n - e_f
        g0[1,0] << -V
        g0[0,1] << -V
        g0[1,1] << inverse(Wilson(D))

you will generate a hybridization function (see also S.Delta_tau after your S.solve(..) call),
that vanishes for all components but [1,1]. The sampling of the Green function is however based
on the removal of hybridization lines from all the diagrams that the MC Markov Chain explores.

Given that Delta[0,0] vanishes, you will never see any diagrams that contain the line,
incorrectly resulting in a vanishing G[0,0]. The fact that you are still seeing a constant background
in your G[0,0] is attributed to the fact that the code manually enforces the discontinuity that
is to be expected for the diagonal entries of G_tau, c.f.

cthyb/c++/triqs_cthyb/measures/G_tau.cpp

Line 73 in 2b155eb

G_tau_block[last] = -1. - G_tau_block[0]; // Enforce 1/iw discontinuity (nb. matrix eq.)

The proper course of action here would be to add a warning in cases where the incorrect G sampling will show up.
This could be done for example by comparing the hybridization function that is extracted from G0_iw against the
atomic Green function that can be calculated from h_loc directly. If the atomic Green function has finite entries where
the Hybridization function vanishes, we can warn the user. @krivenko Would you agree?

Another important route will be the addition of Worm sampling to the TRIQS cthyb code.
This sampling method does not have the issue described above. The implementation of this sampling
approach is currently on-going.

[bfield1]

Dear @Wentzell ,

Thank you for your reply. That makes sense. My hybridisation function indeed has zero elements (see plots below, and as we'd expect from the definition of the hybridisation function), which would explain my anomalous results.
I'll leave you to figure out the best course of action for cthyb. Meanwhile, I'll reformulate my problem to avoid vanishing hybridisation functions.

And just for comparison: the impurity-only solution.