fit_gf_dlr should allow a mesh to be specified.

[andrewkhardy]

I'd like to propose that the fit_gf_dlr function either have a symmetrize=true keyword, like MeshDLRImFreq() or even allow you to provide the DLR grid you're fitting over.

Otherwise it's a bit tricky to combine quantities fit to a DLR grid and quantities computed on a DLR grid directly.

The simplest example I can think of that shows the issue comes from adapting the tutorials. Here is my best example of a MWE with issues

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 - t**2 * g)

S.solve(U = U)

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

for block, g in S_CTHYB.G_iw:
    S_CTHYB.G0_iw[block] << inverse( iOmega_n + U/2.0 - t**2 * g)

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
)
G_iw_copy = S.G_iw.copy()
G_tau_copy = S.G_tau.copy()
Sigma_iw_copy = S.Sigma_iw.copy()

for block, g in G_iw_copy:
    G_dlr = fit_gf_dlr(S_CTHYB.G_tau[block], w_max=w_max, eps=dlr_error)
    g = make_gf_dlr_imfreq(G_dlr)
    G_tau_copy[block] = make_gf_dlr_imtime(g)

This arrives at the error

AssertionError: Green function meshes are not compatible:
  DLR imtime mesh of size 34 with beta = 20, statistic = Fermion, w_max = 10, eps = 1e-09
and
  DLR imtime mesh of size 33 with beta = 20, statistic = Fermion, w_max = 10, eps = 1e-09

The general idea being that if you would like to evaluate G on a smaller grid and then take Sigma or G from the solver, you can run into problems.

[andrewkhardy]

@jasonkaye might be interested.

[the-hampel]

Hi @andrewkhardy ,

did you see that the unstable branch of triqs has this option: https://triqs.github.io/triqs/unstable/documentation/python_api/triqs.gf.gf_factories.fit_gf_dlr.html . I think we discussed this at some point and I think Nils or someone else added it already. It is just not yet in any release... But handing over a mesh is not yet possible.

Best,
alex

[andrewkhardy]

Hi @the-hampel,

I will check the unstable branch next time! Thanks a lot. Would handing a specific mesh be useful, if you want to take it from a Green's function object directly for example? Perhaps specifying the error and range is effective enough, in which case I will close the issue.