DLR interpolation leads to tail issues

[andrewkhardy]

Continuing my attempt to extend DLR for everything, if you construct a hybridization $\mathcal{G}_0$ to feed to say CTHYB even from very precise non-interacting starting point evaluated on a DLR grid, the tail of $\mathcal{G}_0$ has issues. In the below MWE, you have the warning

WARNING: max(abs(imag(S.Delta_tau[0][0, 0]))) = 1.82808e-10 setting to zero.

Feeding this into a full DMFT loop eventually leads to errors such as (S.G0_iw violates fundamental Green Function property G0(iw)[i,j] = G0(-iw)*[j,i] ) and the code eventually crashing.

Attempting at least the simplest tail-fitting (taken from here for example) does not alleviate this problem.

I believe this is a core TRIQS issue rather than a CTHYB issue, but I'm not sure the best place for the solution. Perhaps there is just some hermiticity / fitting post-processing that should be added to the function make_gf_imfreq ?

from triqs.gf import *
import numpy as np
from triqs.operators import *
from triqs_cthyb import Solver
class IPTSolver:
    def __init__(self, beta):
        self.beta = beta

        # Matsubara frequency Green's functions
        iw_mesh = MeshDLRImFreq(beta=beta, statistic='Fermion', w_max= 10, eps =10e-10, symmetrize = True)
        self.G_iw = BlockGf(mesh=iw_mesh, gf_struct = [('up',1), ('down',1)] )
        self.G0_iw = self.G_iw.copy() # self.G0 will be set by the user after initialization
        self.Sigma_iw = self.G_iw.copy()

        # Imaginary time
        tau_mesh = MeshDLRImTime(beta=beta, statistic='Fermion', w_max= 10, eps =10e-10, symmetrize = True)
        self.G0_tau = BlockGf(mesh=tau_mesh, gf_struct = [('up',1), ('down',1)] )
        self.G_tau = self.G0_tau.copy()
        self.Sigma_tau = self.G0_tau.copy()

    def solve(self, U):
        self.G0_tau << make_gf_dlr_imtime(self.G0_iw)
        self.Sigma_tau << (U**2) * self.G0_tau * self.G0_tau * self.G0_tau
        self.Sigma_iw << make_gf_dlr_imfreq(self.Sigma_tau)

        # Dyson
        self.G_iw << inverse(inverse(self.G0_iw) - self.Sigma_iw)

t = 1.0
U = 5.0
beta = 20
n_iw = 2**12+3
w_max= 10 
dlr_error =10e-10
S = IPTSolver(beta = beta)
S.G_iw << SemiCircular(2*t)
S_CTHYB = Solver(beta = beta, gf_struct = [('up',1), ('down',1)] , n_iw = n_iw)

#n_loops = 1
#for i in range(n_loops):
for block, g in S.G_iw:
    S.G0_iw[block] << inverse( iOmega_n + U/2.0 - t**2 * g)

for block, g0 in S.G0_iw:
    S_CTHYB.G0_iw[block] = make_gf_imfreq(S.G0_iw[block], n_iw = n_iw)


S_CTHYB.solve(h_int = U * n('up',0) * n('down',0),   # Local Hamiltonian
n_cycles  = 50000,                           # Number of QMC cycles
n_warmup_cycles = 5000,                      # Warmup cycles
)`

[the-hampel]

Hi @andrewkhardy ,

I had a look at your script and I am a bit confused. Right now it seems that there are 2 solver objects here: IPT and cthyb. And somehow I miss how the DMFT loop goes here. S_CTHYB.solve is called, but S_CTHYB.G_iw is not used? Maybe the important step missing in the code above (which you probably did but maybe forgot to share) is how you go from S_CTHYB.G_iw or S_CTHYB.G_tau to DLR space? This can be done via the fitting functions for DLR, is this what you tried? If you share such script with me I can try to reproduce the problem. In principle, apart from uncertainty on how to get the DLR result from cthyb the workflow should work, especially for a bethelattice where you do not need to compute Sigma.

See for example the scripts here https://github.com/flatironinstitute/DLR_DMFT_scripts/blob/master/srvo3_dlr_delta_interface.py for our paper https://arxiv.org/abs/2301.07764 , in which we looked at a generic hubbard square lattice, SVO, and SRO. But we did not try there for the Bethe Lattice.

[andrewkhardy]

Hi Alex @the-hampel,

Thanks for the detailed response. Just to be clear about the goal of this post, at Jason's suggestion, I am attempting to put together a "simple" DLR tutorial, so I was trying to use code directly from the tutorials, which perhaps made the example much more convoluted than it needed to be.
Let me attach a more detailed example that matches a more likely use-case, still trying to match the tutorial examples as closely as possible, if you are willing to look at it. I receive the same Solver warnings until the script crashes. I believe (although happy to be corrected) the whole error arises from naively extrapolating a square lattice hybridized Green's function. This error does not occur if a regular mesh is used.

I will update this if I am able to understand what goes wrong based on the SVO example you provided.

from triqs.gf import *
import numpy as np
from triqs.operators import *
from triqs_cthyb import Solver

params = {
    "beta": 10.0,
    "t": 1.0,
    "U" : 5.0,
    "beta" : 20,
    "n_iw" : 2**12+3,
    "w_max" : 10,
    "dlr_err" :10e-10,
    "mu": 0.0,
    "n_iter" : 10,
    "n_k" : 1000,
    "n_bins" : 50,
    "n_cycles" : 50000,
    "length_cycle" : 50,
    "n_warmup_cycles" : 5000,
    "alpha" : 0.75,

}

class DLR_GreensFunctions:
    def __init__(self, params):
        self.gf_struct =[('up',1), ('dn',1)]
        iw_mesh = MeshDLRImFreq(beta=params["beta"], statistic='Fermion', w_max= params["w_max"], eps = params["dlr_err"])#, symmetrize = True)
        self.G_iw   = BlockGf(mesh=iw_mesh, gf_struct=self.gf_struct)
        self.G_tau =  make_gf_dlr_imtime(self.G_iw)
        self.G0_iw = self.G_iw.copy()
        self.G0_tau = self.G_tau.copy()
        self.Sigma_iw = self.G_iw.copy()
        self.Sigma_tau = self.G_tau.copy()
        self.Delta_iw = self.G_iw.copy()
        self.Delta_tau = self.G_tau.copy()

# initializing the DLR Green's functions
GF = DLR_GreensFunctions(params)
# creating density of states
k_linear = np.linspace(-np.pi, np.pi, params["n_k"], endpoint=False)
kx, ky = np.meshgrid(k_linear, k_linear)
epsk = -2 * params["t"] * (np.cos(kx) + np.cos(ky))
h = np.histogram(epsk, bins=params["n_bins"], density=True)
epsilon = 0.5 * (h[1][0:-1] + h[1][1:])
rho = h[0]
delta = h[1][1]-h[1][0]
# interacting Hamiltonian
H_U = params["U"] * n('up',0) * n('dn',0)
# self-consistent loop
S_CTHYB = Solver(beta = params["beta"], gf_struct = [('up',1), ('dn',1)] , n_iw = params["n_iw"])

GF.G_iw.zero()

GF.Sigma_iw << params["U"] / 2
params["mu"] = params["U"] / 2

for iter_i in range(params["n_iter"]):

    if iter_i > 1:
        G_old = GF.G_iw.copy()

    GF.G_iw.zero()

    for spin in ['up', 'dn']:
        for i in range(params["n_bins"]):
            GF.G_iw['%s'%(spin)] += rho[i] * delta * \
                inverse(iOmega_n + params["mu"] - epsilon[i] - GF.Sigma_iw['%s'%(spin)])

    if iter_i > 1:
        for block, g in GF.G_iw:
            GF.G_iw[block] =   (params["alpha"]) * GF.G_iw[block] + (1- params["alpha"] )* G_old[block]


    for block, g0 in GF.G0_iw:
        g0 << inverse(inverse(GF.G_iw[block]) + GF.Sigma_iw[block])
        
    for block, g0 in S_CTHYB.G0_iw:
        g0 << make_gf_imfreq(GF.G0_iw[block], n_iw = params["n_iw"])


    n_min = int(1*params["beta"]+5)
    n_max = int(2*(1*params["beta"]+10))
    S_CTHYB.solve(h_int = H_U,
            n_cycles  = params["n_cycles"],
            length_cycle = params["length_cycle"],
            n_warmup_cycles = params["n_warmup_cycles"],
            measure_density_matrix=True, 
            use_norm_as_weight = True,
            perform_tail_fit=True,
            fit_max_moment = 5,
            fit_min_n = n_min,
            fit_max_n = n_max
            )

    # extracting Sigma implicitly 

    for block, g in GF.G_iw:
        G_dlr = fit_gf_dlr(S_CTHYB.G_tau[block], w_max=params["w_max"], eps=params["dlr_err"])
        g << make_gf_dlr_imfreq(G_dlr)
        GF.G_tau[block] << make_gf_dlr_imtime(g)
        GF.Sigma_iw << inverse(GF.G0_iw[block]) - inverse(GF.G_iw[block])

[the-hampel]

Hi Andrew,

In this example the problem arises from the way you get GF.Sigma_iw :

for block, g in GF.G_iw:
        G_dlr = fit_gf_dlr(S_CTHYB.G_tau[block], w_max=params["w_max"], eps=params["dlr_err"])
        g << make_gf_dlr_imfreq(G_dlr)
        GF.G_tau[block] << make_gf_dlr_imtime(g)
        GF.Sigma_iw << inverse(GF.G0_iw[block]) - inverse(GF.G_iw[block])

this will lead to instabilities. If the tail properties of this object is not "good" or "treated" the following will give you the problematic behavior you described:

for block, g0 in GF.G0_iw:
        g0 << inverse(inverse(GF.G_iw[block]) + GF.Sigma_iw[block])

Now the moments of Sigma and G_iw^-1 have to cancel such that their inverse g0 becomes a proper Gf again (tail goes like 1/iwn) etc. If this is not the case the following step will give you garbage:

g0 << make_gf_imfreq(GF.G0_iw[block], n_iw = params["n_iw"])

this goes horribly wrong if the moments of the DLR objects are not such they are physical. This is the only point in the script where you demand that the DLR object is indeed representing a physical Gf object that has a positive spectrum and the spectrum falls off as 1/w. Anyway, to stabilize this loop you need to treat the tail of GF.Sigma. You can probably visually see this by plotting the resulting g0. It will at some point look very funky.

In the scripts I gave you this is done by actually doing a tail fit on the GF.Sigma before it is further used - this was done because we wanted to only test if the actual DMFT step can be done on the DLR grid. However, in this paper https://arxiv.org/abs/2310.01266 we found that the best way to get Sigma on the DLR mesh is via directly minimizing dyson's equation via an objective function, i.e. https://triqs.github.io/triqs/unstable/documentation/python_api/triqs.gf.dlr_crm_dyson_solver.minimize_dyson.html#triqs.gf.dlr_crm_dyson_solver.minimize_dyson . Additionally one should measure in cthyb the moments of Sigma and then use them in the tail-fit or in the crm_dyson_solver. Trying to impose moments on Gimp and then extracting Sigma via normal Dyson equation always gave inferior results.

However, all this does not explain why a more simple Bethe lattice example has problems. Here, Sigma is never needed and this should be much more stable and doable via fit_gf_dlr(S_CTHYB.G_tau[block] . Maybe we should tackle that problem first and then work towards the version you posted above? I am happy to help in either of the two to make this work. :-)

Ah and I would always use the symmetrized version of the DLR mesh. I found the unsymmetrized DLR meshes can give a lot of trouble when it comes to physicality of Gf objects.

Best,
Alex

[andrewkhardy]

Hi Alex @the-hampel to your last point, I would need to install the unstable version to get symmetrized DLR to work as you pointed out.

I also tried the very nice constrained optimization fitting for Sigma at one point, but was never able to get sensible results out of it (which also does not take symmetrize = true by the way as an argument in v3.3 at least). Let me post a simpler Bethe lattice version to work out the kinks and perhaps these other points will become clearer.