From d488ebe42f464c06f21fa399ab8656360473404e Mon Sep 17 00:00:00 2001 From: Chia-Nan Yeh Date: Wed, 12 Jun 2024 09:35:25 -0400 Subject: [PATCH] GW embedding (#80) This PR makes a couple changes on the GW embedding routine. The main changes include: * Use dlr_wmax and dlr_eps in convert_gw_output() * extracting H0_loc and Delta_time from g_weiss numerically on a DLR mesh. * customized dlr_wmax and dlr_eps for crm_dyson_solver to avoid over-fitting --------- Co-authored-by: Alexander Hampel --- python/solid_dmft/dmft_tools/solver.py | 34 ++++++-- .../solid_dmft/gw_embedding/bdft_converter.py | 47 +++++++---- python/solid_dmft/gw_embedding/gw_flow.py | 77 ++++++++++++------- python/solid_dmft/gw_embedding/iaft.py | 48 +++++++++--- python/solid_dmft/io_tools/default.toml | 4 + python/solid_dmft/io_tools/documentation.txt | 12 +++ .../svo_cthyb_basic_crm/dmft_config.toml | 4 +- test/python/svo_ctseg_dyn/dmft_config.toml | 6 +- test/python/svo_gw_emb_dyn/dmft_config.toml | 6 +- 9 files changed, 177 insertions(+), 61 deletions(-) diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py index fe5239bb..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() @@ -1111,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 @@ -1438,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']) @@ -1475,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 = "" +crm_dlr_eps = "" delta_interface = false diag_delta = false fit_max_moment = "" @@ -107,6 +109,8 @@ random_seed = "" [[solver]] type = "ctseg" crm_dyson_solver = false +crm_dlr_wmax = "" +crm_dlr_eps = "" diag_delta = false fit_max_moment = "" fit_max_n = "" 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