From 8112865fc72bf31a83d3e4714185722f3823892b Mon Sep 17 00:00:00 2001 From: Alexander Hampel Date: Wed, 19 Jun 2024 19:19:46 -0400 Subject: [PATCH] [feat] extract Sigma moments when using ctseg * extract Sigma Hartree from measured occupations with high accuracy and store to archive * no more tail fit parameters are necessary for ctseg usage when using the CRM solver --- .../dmft_tools/results_to_archive.py | 3 +- python/solid_dmft/dmft_tools/solver.py | 65 +++++++++++++------ test/python/svo_ctseg_dyn/dmft_config.toml | 3 - test/python/svo_gw_emb_dyn/dmft_config.toml | 3 - 4 files changed, 46 insertions(+), 28 deletions(-) diff --git a/python/solid_dmft/dmft_tools/results_to_archive.py b/python/solid_dmft/dmft_tools/results_to_archive.py index 54cc83b0..2d414466 100644 --- a/python/solid_dmft/dmft_tools/results_to_archive.py +++ b/python/solid_dmft/dmft_tools/results_to_archive.py @@ -105,6 +105,8 @@ def _compile_information(sum_k, general_params, solver_params, solvers, map_imp_ if solver_type_per_imp[icrsh] == 'ctseg': # if legendre was set, that we have both now! write_to_h5['orbital_occupations_{}'.format(icrsh)] = solvers[icrsh].orbital_occupations + write_to_h5['Sigma_Hartree_{}'.format(icrsh)] = solvers[icrsh].Sigma_Hartree + write_to_h5['Sigma_moments_{}'.format(icrsh)] = solvers[icrsh].Sigma_moments if (solver_params[isolvsec]['legendre_fit']): write_to_h5['G_time_orig_{}'.format(icrsh)] = solvers[icrsh].G_time_orig write_to_h5['Gimp_l_{}'.format(icrsh)] = solvers[icrsh].G_l @@ -121,7 +123,6 @@ def _compile_information(sum_k, general_params, solver_params, solvers, map_imp_ if solver_params[isolvsec]['crm_dyson_solver']: write_to_h5['G_time_dlr_{}'.format(icrsh)] = solvers[icrsh].G_time_dlr write_to_h5['Sigma_dlr_{}'.format(icrsh)] = solvers[icrsh].Sigma_dlr - write_to_h5['Sigma_Hartree_{}'.format(icrsh)] = solvers[icrsh].Sigma_Hartree if general_params['h_int_type'][icrsh] == 'dyn_density_density': write_to_h5['D0_time_{}'.format(icrsh)] = solvers[icrsh].triqs_solver.D0_tau write_to_h5['Jperp_time_{}'.format(icrsh)] = solvers[icrsh].triqs_solver.Jperp_tau diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py index 497488fc..80f6647d 100755 --- a/python/solid_dmft/dmft_tools/solver.py +++ b/python/solid_dmft/dmft_tools/solver.py @@ -1384,6 +1384,8 @@ def _ctseg_postprocessing(self): r''' Organize G_freq, G_time, Sigma_freq and G_l from ctseg solver ''' + from triqs.operators.util.extractors import extract_U_dict2, dict_to_matrix + from solid_dmft.postprocessing.eval_U_cRPA_RESPACK import construct_Uijkl def set_Gs_from_G_l(): @@ -1423,6 +1425,43 @@ def set_Gs_from_G_l(): for iorb in range(norb): self.orbital_occupations[block][iorb, iorb] = self.triqs_solver.results.densities[i] + self.orbital_occupations_sumk = self.sum_k.block_structure.convert_matrix(self.orbital_occupations, ish_from=self.icrsh, space_from='solver', space_to='sumk') + self.Sigma_Hartree = {} + self.Sigma_Hartree_sumk = {} + self.Sigma_moments = {} + if mpi.is_master_node(): + # get density density U tensor from solver + U_dict = extract_U_dict2(self.h_int) + norb = get_n_orbitals(self.sum_k)[self.icrsh]['up'] + U_dd = dict_to_matrix(U_dict, gf_struct=self.sum_k.gf_struct_solver_list[self.icrsh]) + # extract Uijij and Uijji terms + Uijij = U_dd[0:norb, norb:2*norb] + Uijji = Uijij - U_dd[0:norb, 0:norb] + # and construct full Uijkl tensor + Uijkl = construct_Uijkl(Uijij, Uijji) + + # now calculated Hartree shift via + # \Sigma^0_{\alpha \beta} = \sum_{i j} n_{i j} \left( 2 U_{\alpha i \beta j} - U_{\alpha i j \beta} \right) + for block, norb in self.sum_k.gf_struct_sumk[self.icrsh]: + self.Sigma_Hartree_sumk[block] = np.zeros((norb, norb),dtype=float) + for iorb, jorb in product(range(norb), repeat=2): + for inner in range(norb): + self.Sigma_Hartree_sumk[block][iorb,jorb] += self.orbital_occupations_sumk[block][inner, inner].real * ( 2*Uijkl[iorb, inner, jorb, inner].real - Uijkl[iorb, inner, inner, jorb].real ) + + # convert to solver block structure + self.Sigma_Hartree = self.sum_k.block_structure.convert_matrix(self.Sigma_Hartree_sumk, ish_from=self.icrsh, space_from='sumk', space_to='solver') + + # use degenerate shell information to symmetrize + self.sum_k.symm_deg_gf(self.Sigma_Hartree, ish=self.icrsh) + + # create moments array from this + for block, hf_val in self.Sigma_Hartree.items(): + self.Sigma_moments[block] = np.array([hf_val]) + + self.Sigma_Hartree = mpi.bcast(self.Sigma_Hartree) + self.Sigma_moments = mpi.bcast(self.Sigma_moments) + self.Sigma_Hartree_sumk = mpi.bcast(self.Sigma_Hartree_sumk) + if mpi.is_master_node(): # create empty moment container (list of np.arrays) Gf_known_moments = make_zero_tail(self.G_freq,n_moments=2) @@ -1450,10 +1489,8 @@ def set_Gs_from_G_l(): fit_max_n=self.solver_params['fit_max_n'], fit_min_w=self.solver_params['fit_min_w'], fit_max_w=self.solver_params['fit_max_w'], - fit_max_moment=self.solver_params['fit_max_moment'],) - self.Sigma_Hartree = {} - for block in tail.keys(): - self.Sigma_Hartree[block] = tail[block][0] + fit_max_moment=self.solver_params['fit_max_moment'], + fit_known_moments=self.Sigma_moments) # recompute G_freq from Sigma with fitted tail self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq) @@ -1482,7 +1519,6 @@ def set_Gs_from_G_l(): from triqs.gf.dlr_crm_dyson_solver import minimize_dyson mpi.report('\nCRM Dyson solver to extract Σ impurity\n') - self.Sigma_Hartree = {} # fit QMC G_tau to DLR if mpi.is_master_node(): if self.solver_params['crm_dlr_wmax'] is not None: @@ -1515,44 +1551,32 @@ def set_Gs_from_G_l(): # minimize dyson for the first entry of each deg shell self.Sigma_dlr = self.sum_k.block_structure.create_gf(ish=self.icrsh, gf_function=Gf, mesh=mesh_dlr_iw, space='solver') # without any degenerate shells we run the minimization for all blocks - _, tail = _fit_tail_window(Sigma_iw, - fit_min_n=self.solver_params['fit_min_n'], - fit_max_n=self.solver_params['fit_max_n'], - fit_min_w=self.solver_params['fit_min_w'], - fit_max_w=self.solver_params['fit_max_w'], - fit_max_moment=self.solver_params['fit_max_moment'],) if self.sum_k.deg_shells[self.icrsh] == []: for block, gf in self.Sigma_dlr: - # tail, err = Sigma_iw[block].fit_hermitian_tail() - self.Sigma_Hartree[block] = tail[block][0] np.random.seed(85281) print('Minimizing Dyson via CRM for Σ[block]:', block) gf, _, _ = minimize_dyson(G0_dlr=G0_dlr_iw[block], G_dlr=G_dlr[block], - Sigma_moments=tail[block][0:1] + Sigma_moments=self.Sigma_moments[block] ) else: for deg_shell in self.sum_k.deg_shells[self.icrsh]: for i, block in enumerate(deg_shell): if i == 0: - # tail, err = Sigma_iw[block].fit_hermitian_tail() - self.Sigma_Hartree[block] = tail[block][0] np.random.seed(85281) print('Minimizing Dyson via CRM for Σ[block]:', block) self.Sigma_dlr[block], _, _ = minimize_dyson(G0_dlr=G0_dlr_iw[block], G_dlr=G_dlr[block], - Sigma_moments=tail[block][0:1] + Sigma_moments=self.Sigma_moments[block] ) sol_block = block else: print(f'Copying result from first block of deg list to {block}') self.Sigma_dlr[block] << self.Sigma_dlr[sol_block] - self.Sigma_Hartree[block] = tail[block][0] self.Sigma_freq[block] = make_gf_imfreq(self.Sigma_dlr[block], n_iw=self.general_params['n_iw']) - self.Sigma_freq[block] += tail[block][0] + self.Sigma_freq[block] += self.Sigma_moments[block][0] - # recompute G_freq from Sigma with fitted tail self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq) print('\n') @@ -1560,7 +1584,6 @@ def set_Gs_from_G_l(): mpi.barrier() self.Sigma_freq = mpi.bcast(self.Sigma_freq) self.Sigma_dlr = mpi.bcast(self.Sigma_dlr) - self.Sigma_Hartree = mpi.bcast(self.Sigma_Hartree) self.G_time_dlr = mpi.bcast(self.G_time_dlr) self.G_freq = mpi.bcast(self.G_freq) else: diff --git a/test/python/svo_ctseg_dyn/dmft_config.toml b/test/python/svo_ctseg_dyn/dmft_config.toml index cd54184f..590d0a8c 100644 --- a/test/python/svo_ctseg_dyn/dmft_config.toml +++ b/test/python/svo_ctseg_dyn/dmft_config.toml @@ -33,9 +33,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 n_tau_k = 10001 [gw] diff --git a/test/python/svo_gw_emb_dyn/dmft_config.toml b/test/python/svo_gw_emb_dyn/dmft_config.toml index cd16e90c..e2ae5fcf 100644 --- a/test/python/svo_gw_emb_dyn/dmft_config.toml +++ b/test/python/svo_gw_emb_dyn/dmft_config.toml @@ -30,8 +30,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 measure_nnt = true