diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py
index 45e62728..85f57a5c 100755
--- a/python/solid_dmft/dmft_tools/solver.py
+++ b/python/solid_dmft/dmft_tools/solver.py
@@ -469,6 +469,14 @@ def solve(self, **kwargs):
             self.triqs_solver.solve(h_int=self.h_int, **self.triqs_solver_params)
             # *************************************
 
+            # dump Delta_tau constructed internally from cthyb when delta_interface = False
+            if self.general_params['store_solver'] and mpi.is_master_node():
+                with HDFArchive(self.general_params['jobname'] + '/' + self.general_params['seedname'] + '.h5',
+                                'a') as archive:
+                    if not self.solver_params['delta_interface']:
+                        archive['DMFT_input/solver/it_-1'][f'Delta_time_{self.icrsh}'] = self.triqs_solver.Delta_tau
+            mpi.barrier()
+
             # call postprocessing
             self._cthyb_postprocessing()
 
@@ -744,6 +752,9 @@ def make_positive_definite(G):
                         archive['DMFT_input/solver/it_-1'][f'Uloc_dlr_2idx_{self.icrsh}'] = Uloc_dlr_2idx
                         archive['DMFT_input/solver/it_-1'][f'Uloc_dlr_2idx_prime_{self.icrsh}'] = Uloc_dlr_2idx_prime
             mpi.barrier()
+
+            # turn of problematic move in ctseg until fixed!
+            self.triqs_solver_params['move_move_segment'] = False
             # Solve the impurity problem for icrsh shell
             # *************************************
             self.triqs_solver.solve(h_int=self.h_int, chemical_potential=chemical_potential, **self.triqs_solver_params)
@@ -1108,8 +1119,16 @@ def set_Gs_from_G_l():
                 mpi.report('\nCRM Dyson solver to extract Σ impurity\n')
                 # fit QMC G_tau to DLR
                 if mpi.is_master_node():
-                    G_dlr = fit_gf_dlr(self.triqs_solver.G_tau,
-                                       w_max=self.general_params['dlr_wmax'], eps=self.general_params['dlr_eps'])
+                    if self.solver_params['crm_dlr_wmax'] is not None:
+                        dlr_wmax = self.solver_params['crm_dlr_wmax']
+                    else:
+                        dlr_wmax = self.general_params['dlr_wmax']
+                    if self.solver_params['crm_dlr_eps'] is not None:
+                        dlr_eps = self.solver_params['crm_dlr_eps']
+                    else:
+                        dlr_eps = self.general_params['dlr_eps']
+                    mpi.report(f"crm_dyson_solver with (wmax, eps) = ({dlr_wmax}, {dlr_eps}). ")
+                    G_dlr = fit_gf_dlr(self.triqs_solver.G_tau, w_max=dlr_wmax, eps=dlr_eps)
                     self.G_time_dlr = make_gf_dlr_imtime(G_dlr)
 
                     # assume little error on G0_iw and use to get G0_dlr
@@ -1435,8 +1454,16 @@ def set_Gs_from_G_l():
             self.Sigma_Hartree = {}
             # fit QMC G_tau to DLR
             if mpi.is_master_node():
-                G_dlr = fit_gf_dlr(self.triqs_solver.results.G_tau,
-                                   w_max=self.general_params['dlr_wmax'], eps=self.general_params['dlr_eps'])
+                if self.solver_params['crm_dlr_wmax'] is not None:
+                    dlr_wmax = self.solver_params['crm_dlr_wmax']
+                else:
+                    dlr_wmax = self.general_params['dlr_wmax']
+                if self.solver_params['crm_dlr_eps'] is not None:
+                    dlr_eps = self.solver_params['crm_dlr_eps']
+                else:
+                    dlr_eps = self.general_params['dlr_eps']
+                mpi.report(f"crm_dyson_solver with (wmax, eps) = ({dlr_wmax}, {dlr_eps}). ")
+                G_dlr = fit_gf_dlr(self.triqs_solver.results.G_tau, w_max=dlr_wmax, eps=dlr_eps)
                 self.G_time_dlr = make_gf_dlr_imtime(G_dlr)
                 self.G_freq = make_gf_imfreq(G_dlr, n_iw=self.general_params['n_iw'])
 
@@ -1472,7 +1499,7 @@ def set_Gs_from_G_l():
                         gf, _, _ = minimize_dyson(G0_dlr=G0_dlr_iw[block],
                                                   G_dlr=G_dlr[block],
                                                   Sigma_moments=tail[block][0:1]
-                                                    )
+                                                  )
                 else:
                     for deg_shell in self.sum_k.deg_shells[self.icrsh]:
                         for i, block in enumerate(deg_shell):
diff --git a/python/solid_dmft/gw_embedding/bdft_converter.py b/python/solid_dmft/gw_embedding/bdft_converter.py
index ca0a4572..de6bebc8 100644
--- a/python/solid_dmft/gw_embedding/bdft_converter.py
+++ b/python/solid_dmft/gw_embedding/bdft_converter.py
@@ -90,6 +90,18 @@ def _get_dlr_from_IR(Gf_ir, ir_kernel, mesh_dlr_iw, dim=2):
     return Gf_dlr
 
 
+def check_iaft_accuracy(Aw, ir_kernel, stats,
+                        beta, dlr_wmax, dlr_prec, data_name):
+    mpi.report('============== DLR mesh check ==============\n')
+    mpi.report(f'Estimating accuracy of the user-defined (wmax, eps) = '
+               f'({dlr_wmax}, {dlr_prec}) for the DLR mesh\n')
+    ir_imp_kernel = IAFT(beta=beta, lmbda=beta * dlr_wmax, prec=dlr_prec)
+    Aw_imp = ir_kernel.w_interpolate(Aw, ir_imp_kernel.wn_mesh('f'), 'f')
+
+    ir_imp_kernel.check_leakage(Aw_imp, stats, data_name, w_input=True)
+    mpi.report('=================== done ===================\n')
+
+
 def calc_Sigma_DC_gw(Wloc_dlr, Gloc_dlr, Vloc, verbose=False):
     r"""
     Calculate the double counting part of the self-energy from the screened Coulomb interaction
@@ -288,7 +300,8 @@ def calc_W_from_Gloc(Gloc_dlr, U):
     return W_dlr
 
 
-def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
+def convert_gw_output(job_h5, gw_h5, dlr_wmax=None, dlr_eps=None,
+                      it_1e=0, it_2e=0, ha_ev_conv = False):
     """
     read bdft output and convert to triqs Gf DLR objects
 
@@ -298,6 +311,10 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
         path to solid_dmft job file
     gw_h5: string
         path to GW checkpoint file for AIMBES code
+    dlr_wmax: float
+        wmax for dlr mesh, defaults to the wmax from the IR basis
+    dlr_eps: float
+        precision for dlr mesh, defaults to the precision from the IR basis
     it_1e: int, optional
         iteration to read from gw_h5 calculation for 1e downfolding, defaults to last iteration
     it_2e: int, optional
@@ -338,6 +355,7 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
         gw_data['beta'] = ar['imaginary_fourier_transform']['beta']
         gw_data['lam'] = ar['imaginary_fourier_transform']['lambda']
         gw_data['gw_wmax'] = gw_data['lam'] / gw_data['beta']
+        gw_data['gw_dlr_wmax'] = gw_data['gw_wmax'] if dlr_wmax is None else dlr_wmax
         gw_data['number_of_spins'] = ar['system/number_of_spins']
         assert gw_data['number_of_spins'] == 1, 'spin calculations not yet supported in converter'
 
@@ -345,13 +363,13 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
         if prec == 'high':
             # set to highest DLR precision possible
             gw_data['gw_ir_prec'] = 1e-15
-            gw_data['gw_dlr_prec'] = 1e-13
+            gw_data['gw_dlr_prec'] = 1e-13 if dlr_eps is None else dlr_eps
         elif prec == 'mid':
             gw_data['gw_ir_prec'] = 1e-10
-            gw_data['gw_dlr_prec'] = 1e-10
+            gw_data['gw_dlr_prec'] = 1e-10 if dlr_eps is None else dlr_eps
         elif prec == 'low':
             gw_data['gw_ir_prec'] = 1e-6
-            gw_data['gw_dlr_prec'] = 1e-6
+            gw_data['gw_dlr_prec'] = 1e-6 if dlr_eps is None else dlr_eps
 
         # 1 particle properties
         g_weiss_wsIab = ar[f'downfold_1e/iter{it_1e}']['g_weiss_wsIab']
@@ -393,19 +411,26 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
     # get IR object
     mpi.report('create IR kernel and convert to DLR')
     # create IR kernel
+    mpi.report("")
     ir_kernel = IAFT(beta=gw_data['beta'], lmbda=gw_data['lam'], prec=gw_data['gw_ir_prec'])
 
+    if dlr_wmax is not None or dlr_eps is not None:
+        # check user-defined dlr_wmax and dlr_eps for the impurity
+        check_iaft_accuracy(g_weiss_wsIab, ir_kernel, 'f',
+                            gw_data['beta'], gw_data['gw_dlr_wmax'], gw_data['gw_dlr_prec'],
+                            "fermionic Weiss field g")
+
     gw_data['mesh_dlr_iw_b'] = MeshDLRImFreq(
         beta=gw_data['beta']/conv_fac,
         statistic='Boson',
-        w_max=gw_data['gw_wmax']*conv_fac,
+        w_max=gw_data['gw_dlr_wmax']*conv_fac,
         eps=gw_data['gw_dlr_prec'],
         symmetrize=True
     )
     gw_data['mesh_dlr_iw_f'] = MeshDLRImFreq(
         beta=gw_data['beta']/conv_fac,
         statistic='Fermion',
-        w_max=gw_data['gw_wmax']*conv_fac,
+        w_max=gw_data['gw_dlr_wmax']*conv_fac,
         eps=gw_data['gw_dlr_prec'],
         symmetrize=True
     )
@@ -473,17 +498,13 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
     mpi.report(f'writing results in {job_h5}/DMFT_input')
 
     with HDFArchive(job_h5, 'a') as ar:
-        if 'DMFT_results' in ar and 'iteration_count' in ar['DMFT_results']:
-            it = ar['DMFT_results']['iteration_count'] + 1
-        else:
-            it = 1
         if 'DMFT_input' not in ar:
             ar.create_group('DMFT_input')
-        if f'iter{it}' not in ar['DMFT_input']:
-            ar['DMFT_input'].create_group(f'iter{it}')
+        if f'iter{it_1e}' not in ar['DMFT_input']:
+            ar['DMFT_input'].create_group(f'iter{it_1e}')
 
         for key, value in gw_data.items():
-            ar[f'DMFT_input/iter{it}'][key] = value
+            ar[f'DMFT_input/iter{it_1e}'][key] = value
 
     mpi.report(f'finished writing results in {job_h5}/DMFT_input')
     return gw_data, ir_kernel
diff --git a/python/solid_dmft/gw_embedding/gw_flow.py b/python/solid_dmft/gw_embedding/gw_flow.py
index acb39891..10f39347 100644
--- a/python/solid_dmft/gw_embedding/gw_flow.py
+++ b/python/solid_dmft/gw_embedding/gw_flow.py
@@ -34,6 +34,8 @@
 from triqs.utility import mpi
 from triqs.gf.tools import inverse
 from triqs.gf import (
+    iOmega_n,
+    fit_hermitian_tail,
     Gf,
     BlockGf,
     make_hermitian,
@@ -253,6 +255,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
         gw_data, ir_kernel = convert_gw_output(
             general_params['jobname'] + '/' + general_params['seedname'] + '.h5',
             gw_params['h5_file'],
+            general_params['dlr_wmax'], general_params['dlr_eps'],
             it_1e = gw_params['it_1e'],
             it_2e = gw_params['it_2e'],
         )
@@ -358,21 +361,32 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
 
             # prepare solver input
             imp_eal = sumk.block_structure.convert_matrix(gw_params['Hloc0'][ish], ish_from=ish, space_from='sumk', space_to='solver')
-            delta_dlr = sumk.block_structure.convert_gf(gw_params['delta_dlr'][ish], ish_from=ish, space_from='sumk', space_to='solver')
-            # fill Delta_time from Delta_freq sumk to solver
-            for name, g0 in delta_dlr:
-                # make non-interacting impurity Hamiltonian hermitian
-                imp_eal[name] = (imp_eal[name] + imp_eal[name].T.conj())/2
+            for name, g0 in G0_dlr:
+                # Estimate the HF shift
+                G0_iw = solvers[ish].G0_freq[name]
+                Delta_iw = Gf(mesh=G0_iw.mesh, target_shape=G0_iw.target_shape)
+                Delta_iw << iOmega_n - inverse(G0_iw)
+                tail, err = fit_hermitian_tail(Delta_iw)
+                # overwrite H0 using estimation from high-frequency tail
+                imp_eal[name] = tail[0]
+                mpi.report(f"Tail fitting for extracting the HF shift in g_weiss with error {err}")
+
                 if mpi.is_master_node():
                     print('H_loc0[{:2d}] block: {}'.format(ish, name))
                     fmt = '{:11.7f}' * imp_eal[name].shape[0]
                     for block in imp_eal[name]:
                         print((' '*11 + fmt).format(*block.real))
 
+                G0_dlr_iw = make_gf_dlr_imfreq(g0)
+                Delta_dlr_iw = Gf(mesh=G0_dlr_iw.mesh, target_shape=g0.target_shape)
+                for iw in G0_dlr_iw.mesh:
+                    Delta_dlr_iw[iw] = iw.value - inverse(G0_dlr_iw[iw]) - imp_eal[name]
+                Delta_dlr = make_gf_dlr(Delta_dlr_iw)
+                # create now full delta(tau)
+                Delta_tau = make_hermitian(make_gf_imtime(Delta_dlr, n_tau=general_params['n_tau']))
+
                 # without SOC delta_tau needs to be real
                 if not sumk.SO == 1:
-                    # create now full delta(tau)
-                    Delta_tau = make_hermitian(make_gf_imtime(delta_dlr[name], n_tau=general_params['n_tau']))
                     solvers[ish].Delta_time[name] << Delta_tau.real
                 else:
                     solvers[ish].Delta_time[name] << Delta_tau
@@ -444,37 +458,45 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
                 sumk.symm_deg_gf(Sigma_dlr_iw[ish],ish=ish)
                 Sigma_dlr[ish] = make_gf_dlr(Sigma_dlr_iw[ish])
 
-                for i, (block, gf) in enumerate(Sigma_dlr[ish]):
-                    # print Hartree shift
-                    print('Σ_HF {}'.format(block))
-                    fmt = '{:11.7f}' * solvers[ish].Sigma_Hartree[block].shape[0]
-                    for vhf in solvers[ish].Sigma_Hartree[block]:
-                        print((' '*11 + fmt).format(*vhf.real))
-
-                # average Hartree shift if not magnetic
-                if not general_params['magnetic']:
-                    if general_params['enforce_off_diag']:
-                        solvers[ish].Sigma_Hartree['up_0'] = 0.5*(solvers[ish].Sigma_Hartree['up_0']+
-                                                                  solvers[ish].Sigma_Hartree['down_0'])
-                        solvers[ish].Sigma_Hartree['down_0'] = solvers[ish].Sigma_Hartree['up_0']
-                    else:
-                        for iorb in range(gw_params['n_orb'][ish]):
-                            solvers[ish].Sigma_Hartree[f'up_{iorb}'] = 0.5*(solvers[ish].Sigma_Hartree[f'up_{iorb}']+
-                                                                          solvers[ish].Sigma_Hartree[f'down_{iorb}'])
-                            solvers[ish].Sigma_Hartree[f'down_{iorb}'] = solvers[ish].Sigma_Hartree[f'up_{iorb}']
+            for i, (block, gf) in enumerate(Sigma_dlr[ish]):
+                # print Hartree shift
+                print('Σ_HF {}'.format(block))
+                fmt = '{:11.7f}' * solvers[ish].Sigma_Hartree[block].shape[0]
+                for vhf in solvers[ish].Sigma_Hartree[block]:
+                    print((' '*11 + fmt).format(*vhf.real))
+
+            # average Hartree shift if not magnetic
+            if not general_params['magnetic']:
+                if general_params['enforce_off_diag']:
+                    solvers[ish].Sigma_Hartree['up_0'] = 0.5*(solvers[ish].Sigma_Hartree['up_0']+
+                                                              solvers[ish].Sigma_Hartree['down_0'])
+                    solvers[ish].Sigma_Hartree['down_0'] = solvers[ish].Sigma_Hartree['up_0']
+                else:
+                    for iorb in range(gw_params['n_orb'][ish]):
+                        solvers[ish].Sigma_Hartree[f'up_{iorb}'] = 0.5*(solvers[ish].Sigma_Hartree[f'up_{iorb}']+
+                                                                      solvers[ish].Sigma_Hartree[f'down_{iorb}'])
+                        solvers[ish].Sigma_Hartree[f'down_{iorb}'] = solvers[ish].Sigma_Hartree[f'up_{iorb}']
 
             iw_mesh = solvers[ish].Sigma_freq.mesh
             # convert Sigma to sumk basis
             Sigma_dlr_sumk = sumk.block_structure.convert_gf(Sigma_dlr[ish], ish_from=ish, space_from='solver', space_to='sumk')
             Sigma_Hartree_sumk = sumk.block_structure.convert_matrix(solvers[ish].Sigma_Hartree, ish_from=ish, space_from='solver', space_to='sumk')
             # store Sigma and V_HF in sumk basis on IR mesh
+            ir_nw_half = len(ir_mesh_idx)//2
             for i, (block, gf) in enumerate(Sigma_dlr_sumk):
                 Vhf_imp_sIab[i,ish] = Sigma_Hartree_sumk[block]
-                for iw in range(len(ir_mesh_idx)):
-                    Sigma_ir[iw,i,ish] = gf(iw_mesh(ir_mesh_idx[iw]))
+                # Make sure sigma_ir[iw].conj() = sigma_ir[-iw]
+                for n in range(ir_nw_half):
+                    iw_pos = ir_nw_half+n
+                    iw_neg = ir_nw_half-1-n
+                    Sigma_ir[iw_pos,i,ish] = gf(iw_mesh(ir_mesh_idx[iw_pos]))
+                    Sigma_ir[iw_neg,i,ish] = gf(iw_mesh(ir_mesh_idx[iw_pos])).conj()
 
                 if not general_params['magnetic']:
                     break
+    if mpi.is_master_node():
+        print("\nChecking impurity self-energy on the IR mesh...")
+        ir_kernel.check_leakage(Sigma_ir, stats='f', name="impurity self-energy", w_input=True)
 
     # Writes results to h5 archive
     if mpi.is_master_node():
@@ -495,7 +517,6 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
             ar[f'downfold_1e/iter{iteration}']['Sigma_imp_wsIab'] = Sigma_ir
             ar[f'downfold_1e/iter{iteration}']['Vhf_imp_sIab'] = Vhf_imp_sIab
 
-
     mpi.report('*** iteration finished ***')
     mpi.report('#'*80)
     mpi.barrier()
diff --git a/python/solid_dmft/gw_embedding/iaft.py b/python/solid_dmft/gw_embedding/iaft.py
index c00050a8..4a9ca312 100644
--- a/python/solid_dmft/gw_embedding/iaft.py
+++ b/python/solid_dmft/gw_embedding/iaft.py
@@ -1,8 +1,9 @@
+import sys
 import numpy as np
 import sparse_ir
 
 """
-Fourier transform on the imaginary axis based on IR basis and the sparse sampling technique.
+Fourier transform on the imaginary axis based on IR basis and the sparse sampling technique.  
 """
 
 
@@ -84,6 +85,7 @@ def __init__(self, beta: float, lmbda: float, prec: float = 1e-15):
         self.Twt_bb = np.dot(Twl_bb, self.Tlt_bb)
 
         print(self)
+        sys.stdout.flush()
 
     def __str__(self):
         return "*******************************\n" \
@@ -94,13 +96,12 @@ def __str__(self):
                "    - lambda = {}\n" \
                "    - nt_f, nw_f = {}, {}\n" \
                "    - nt_b, nw_b = {}, {}\n" \
-               "*******************************".format(self.prec, self.beta, self.lmbda, self.nt_f, self.nw_f,
+               "*******************************\n".format(self.prec, self.beta, self.lmbda, self.nt_f, self.nw_f,
                                                         self.nt_b, self.nw_b)
 
-    def wn_mesh(self, stats, ir_notation= True):
+    def wn_mesh(self, stats: str, ir_notation: bool = True):
         """
         Return Matsubara frequency indices.
-
         :param stats: str
             statistics: 'f' for fermions and 'b' for bosons
         :param ir_notation: bool
@@ -118,7 +119,7 @@ def wn_mesh(self, stats, ir_notation= True):
             wn_mesh = (wn_mesh-1)//2 if stats == 'f' else wn_mesh//2
         return wn_mesh
 
-    def tau_to_w(self, Ot, stats):
+    def tau_to_w(self, Ot, stats: str):
         """
         Fourier transform from imaginary-time axis to Matsubara-frequency axis
         :param Ot: numpy.ndarray
@@ -173,7 +174,7 @@ def w_to_tau(self, Ow, stats):
         Ot = Ot.reshape((Ttw.shape[0],) + Ow_shape[1:])
         return Ot
 
-    def w_interpolate(self, Ow, wn_mesh_interp, stats, ir_notation=True):
+    def w_interpolate(self, Ow, wn_mesh_interp, stats: str, ir_notation: bool = True):
         """
         Interpolate a dynamic object to arbitrary points on the Matsubara axis.
 
@@ -195,7 +196,7 @@ def w_interpolate(self, Ow, wn_mesh_interp, stats, ir_notation=True):
             raise ValueError("Unknown statistics '{}'. "
                              "Acceptable options are 'f' for fermion and 'b' for bosons.".format(stats))
         if ir_notation:
-            wn_indices = wn_mesh_interp
+            wn_indices = np.asarray(wn_mesh_interp)
         else:
             wn_indices = np.array([2*n+1 if stats == 'f' else 2*n for n in wn_mesh_interp], dtype=int)
         Tlw = self.Tlw_ff if stats == 'f' else self.Tlw_bb
@@ -214,7 +215,7 @@ def w_interpolate(self, Ow, wn_mesh_interp, stats, ir_notation=True):
         Ow_interp = Ow_interp.reshape((wn_indices.shape[0],) + Ow_shape[1:])
         return Ow_interp
 
-    def tau_interpolate(self, Ot, tau_mesh_interp, stats):
+    def tau_interpolate(self, Ot, tau_mesh_interp, stats: str):
         """
          Interpolate a dynamic object to arbitrary points on the imaginary-time axis.
 
@@ -244,9 +245,38 @@ def tau_interpolate(self, Ot, tau_mesh_interp, stats):
         Ot_interp = np.dot(Ttt, Ot)
 
         Ot = Ot.reshape(Ot_shape)
-        Ot_interp = Ot_interp.reshape((tau_mesh_interp.shape[0],) + Ot_shape[1:])
+        Ot_interp = Ot_interp.reshape((np.shape(tau_mesh_interp)[0],) + Ot_shape[1:])
         return Ot_interp
 
+    def check_leakage(self, Ot, stats: str, name: str = "", w_input: bool = False):
+        if w_input:
+            Ot_ = self.w_to_tau(Ot, stats)
+            self.check_leakage(Ot_, stats, name, w_input=False)
+            return
+
+        if stats not in self.statisics:
+            raise ValueError("Unknown statistics '{}'. "
+                             "Acceptable options are 'f' for fermion and 'b' for bosons.".format(stats))
+        nts = self.nt_f if stats == 'f' else self.nt_b
+        Tlt = self.Tlt_ff if stats == 'f' else self.Tlt_bb
+        if nts != Ot.shape[0]:
+            raise ValueError("Inconsistency between nts = {} and Ot.shape[0] = {}".format(nts, Ot.shape[0]))
+
+        # coeff_first
+        O_l0_i = np.einsum('t,ti->i', Tlt[0], Ot.reshape(nts, -1))
+        coeff_first = np.max(np.abs(O_l0_i))
+
+        # coeff_last
+        O_lm2_t = np.einsum('lt,ti->li', Tlt[-2:], Ot.reshape(nts, -1))
+        coeff_last = np.max(np.abs(O_lm2_t))
+
+        leakage = coeff_last/coeff_first
+        print("IAFT leakage of {}: {}".format(name, leakage))
+        if leakage >= 1e-8:
+            print("[WARNING] check_leakage: coeff_last/coeff_first = {} >= 1e-8; "
+                  "coeff_last = {}, coeff_first = {}".format(leakage, coeff_last, coeff_first))
+        sys.stdout.flush()
+
 
 if __name__ == '__main__':
     # Initialize IAFT object for given inverse temperature, lambda and precision
diff --git a/python/solid_dmft/io_tools/default.toml b/python/solid_dmft/io_tools/default.toml
index 89220269..4577196f 100644
--- a/python/solid_dmft/io_tools/default.toml
+++ b/python/solid_dmft/io_tools/default.toml
@@ -65,6 +65,8 @@ w_range = [-10, 10]
 [[solver]]
 type = "cthyb"
 crm_dyson_solver = false
+crm_dlr_wmax = "<none>"
+crm_dlr_eps = "<none>"
 delta_interface = false
 diag_delta = false
 fit_max_moment = "<none>"
@@ -107,6 +109,8 @@ random_seed = "<none>"
 [[solver]]
 type = "ctseg"
 crm_dyson_solver = false
+crm_dlr_wmax = "<none>"
+crm_dlr_eps = "<none>"
 diag_delta = false
 fit_max_moment = "<none>"
 fit_max_n = "<none>"
diff --git a/python/solid_dmft/io_tools/documentation.txt b/python/solid_dmft/io_tools/documentation.txt
index b16be00b..359eaa42 100644
--- a/python/solid_dmft/io_tools/documentation.txt
+++ b/python/solid_dmft/io_tools/documentation.txt
@@ -221,6 +221,12 @@ cthyb
 crm_dyson_solver : bool, default = False
         use CRM Dyson solver to extract Sigma_imp from G(tau) (conflict with legendre_fit and tail_fit)
         set dlr_wmax and dlr_eps parameters in general section to use
+crm_dlr_wmax: float, default = None
+        customized dlr_wmax for the crm_dyson_solver. Only used if crm_dyson_solver = True.
+        Set to dlr_wmax if crm_dlr_wmax = None.
+crm_dlr_eps: float, default = None
+        customized dlr_eps for the crm_dyson_solver. Only used if crm_dyson_solver = True.
+        Set to dlr_eps if crm_dlr_eps = None.
 delta_interface : bool, default = False
         use delta interface in cthyb instead of input G0
 diag_delta : bool, default = False
@@ -317,6 +323,12 @@ ctseg
 crm_dyson_solver : bool, default = False
         use CRM Dyson solver to extract Sigma_imp from G(tau) (conflict with legendre_fit and tail_fit)
         set dlr_wmax and dlr_eps parameters in general section to use
+crm_dlr_wmax: float, default = None
+        customized dlr_wmax for the crm_dyson_solver. Only used if crm_dyson_solver = True.
+        Set to dlr_wmax if crm_dlr_wmax = None.
+crm_dlr_eps: float, default = None
+        customized dlr_eps for the crm_dyson_solver. Only used if crm_dyson_solver = True.
+        Set to dlr_eps if crm_dlr_eps = None.
 diag_delta : bool, default = False
         approximate the hybridization function as diagonal when using the delta interface
 fit_max_moment : int, default = None
diff --git a/test/python/svo_cthyb_basic_crm/dmft_config.toml b/test/python/svo_cthyb_basic_crm/dmft_config.toml
index 65cd906f..0ac8e447 100644
--- a/test/python/svo_cthyb_basic_crm/dmft_config.toml
+++ b/test/python/svo_cthyb_basic_crm/dmft_config.toml
@@ -8,8 +8,6 @@ mu_initial_guess = -0.027041
 prec_mu = 0.001
 n_iw = 501
 n_tau = 10001
-dlr_wmax = 4
-dlr_eps = 1e-8
 
 h_int_type = "kanamori"
 U = 8.0
@@ -39,3 +37,5 @@ n_cycles_tot = 1e+4
 imag_threshold = 1e-5
 measure_density_matrix = true
 crm_dyson_solver = true
+crm_dlr_wmax = 4
+crm_dlr_eps = 1e-8
diff --git a/test/python/svo_ctseg_dyn/dmft_config.toml b/test/python/svo_ctseg_dyn/dmft_config.toml
index 250d6d0e..cd54184f 100644
--- a/test/python/svo_ctseg_dyn/dmft_config.toml
+++ b/test/python/svo_ctseg_dyn/dmft_config.toml
@@ -6,8 +6,8 @@ mu_initial_guess = 13.223155
 prec_mu = 0.001
 n_iw = 200
 n_tau = 20000
-dlr_wmax = 2
-dlr_eps = 1e-6
+dlr_wmax = 20
+dlr_eps = 1e-10
 
 h_int_type = "dyn_density_density"
 # h_int_type = "crpa_density_density"
@@ -31,6 +31,8 @@ n_warmup_cycles = 1e+4
 n_cycles_tot = 1e+6
 off_diag_threshold = 1e-4
 crm_dyson_solver = true
+crm_dlr_wmax = 2
+crm_dlr_eps = 1e-6
 fit_min_n = 10
 fit_max_n = 60
 fit_max_moment = 4
diff --git a/test/python/svo_gw_emb_dyn/dmft_config.toml b/test/python/svo_gw_emb_dyn/dmft_config.toml
index d023b6b2..cd16e90c 100644
--- a/test/python/svo_gw_emb_dyn/dmft_config.toml
+++ b/test/python/svo_gw_emb_dyn/dmft_config.toml
@@ -6,8 +6,8 @@ enforce_off_diag = false
 
 n_iw = 1000
 n_tau = 10001
-dlr_wmax = 0.2
-dlr_eps = 1e-6
+dlr_wmax = 10
+dlr_eps = 1e-10
 
 gw_embedding = true
 h_int_type = "dyn_density_density"
@@ -28,6 +28,8 @@ n_cycles_tot = 1e+5
 off_diag_threshold = 1e-5
 diag_delta = false
 crm_dyson_solver = true
+crm_dlr_wmax = 0.2
+crm_dlr_eps = 1e-6
 fit_min_n = 10
 fit_max_n = 100
 fit_max_moment = 4