GW embedding

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

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,20 +355,21 @@ 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'
 
         prec = ar['imaginary_fourier_transform']['prec']
         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

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