diff --git a/python/solid_dmft/dmft_tools/solvers/ctseg_interface.py b/python/solid_dmft/dmft_tools/solvers/ctseg_interface.py index 4fddd4cd..ee83ea41 100644 --- a/python/solid_dmft/dmft_tools/solvers/ctseg_interface.py +++ b/python/solid_dmft/dmft_tools/solvers/ctseg_interface.py @@ -278,8 +278,24 @@ def set_Gs_from_G_l(): # get G_time, G_freq, Sigma_freq from G_l set_Gs_from_G_l() elif self.solver_params['perform_tail_fit']: - mpi.report('Self-energy post-processing algorithm: tail fitting with analytic static impurity self-energy') - self.Sigma_freq = inverse(self.G0_freq) - inverse(self.G_freq) + if not self.solver_params['improved_estimator']: + mpi.report('Self-energy post-processing algorithm: tail fitting with analytic static impurity self-energy') + self.Sigma_freq = inverse(self.G0_freq) - inverse(self.G_freq) + else: + mpi.report('Self-energy post-processing algorithm: ' + 'improved estimator + tail fitting with analytic static imppurity self-energy') + self.F_freq = self.G_freq.copy() + self.F_freq << 0.0 + self.F_time = self.G_time.copy() + self.F_time << self.triqs_solver.results.F_tau + F_known_moments = make_zero_tail(self.F_freq, n_moments=1) + for i, bl in enumerate(self.F_freq.indices): + self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i]) + + for block, fw in self.F_freq: + for iw in fw.mesh: + self.Sigma_freq[block][iw] = self.F_freq[block][iw] / self.G_freq[block][iw] + # without any degenerate shells we run the minimization for all blocks self.Sigma_freq, tail = self._fit_tail_window( self.Sigma_freq, @@ -295,42 +311,20 @@ def set_Gs_from_G_l(): self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq) # if improved estimators are turned on calc Sigma from F_tau, otherwise: - elif self.solver_params['improved_estimator']: - if self.solver_params['perform_tail_fit']: - mpi.report('Self-energy post-processing algorithm: ' - 'improved estimator + tail fitting with analytic static imppurity self-energy') - else: - mpi.report('Self-energy post-processing algorithm: improved estimator') + elif self.solver_params['improved_estimator'] and not self.solver_params['perform_tail_fit']: + mpi.report('Self-energy post-processing algorithm: improved estimator') self.F_freq = self.G_freq.copy() self.F_freq << 0.0 self.F_time = self.G_time.copy() self.F_time << self.triqs_solver.results.F_tau F_known_moments = make_zero_tail(self.F_freq, n_moments=1) - if mpi.is_master_node(): - for i, bl in enumerate(self.F_freq.indices): - self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i]) + for i, bl in enumerate(self.F_freq.indices): + self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i]) - self.F_freq << mpi.bcast(self.F_freq) - self.G_freq << mpi.bcast(self.G_freq) for block, fw in self.F_freq: for iw in fw.mesh: self.Sigma_freq[block][iw] = self.F_freq[block][iw] / self.G_freq[block][iw] - # Further tail fitting on top of Sigma_freq from the improved estimator - if self.solver_params['perform_tail_fit']: - # without any degenerate shells we run the minimization for all blocks - self.Sigma_freq, tail = self._fit_tail_window( - self.Sigma_freq, - 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'], - fit_known_moments=self.Sigma_moments, - ) - # recompute G_freq from Sigma_freq - self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq) - # should this be moved to abstract class? elif self.solver_params['crm_dyson_solver']: mpi.report('Self-energy post-processing algorithm: crm dyson solver') @@ -408,9 +402,6 @@ def set_Gs_from_G_l(): if self.solver_params['measure_pert_order']: self.perturbation_order_histo = self.triqs_solver.results.pert_order_Delta - bin_vec = np.arange(0, self.perturbation_order_histo.data.shape[0]) - self.avg_pert_order = np.sum(bin_vec * self.perturbation_order_histo.data[:]) - if mpi.is_master_node(): - print(f'Average perturbation order: {self.avg_pert_order}') + self.avg_pert_order = self.triqs_solver.results.average_order_Delta return diff --git a/python/solid_dmft/gw_embedding/bdft_converter.py b/python/solid_dmft/gw_embedding/bdft_converter.py index c314f591..8d40c50d 100644 --- a/python/solid_dmft/gw_embedding/bdft_converter.py +++ b/python/solid_dmft/gw_embedding/bdft_converter.py @@ -181,7 +181,12 @@ def tail_fit_g_weiss(g_weiss_wsIab, ir_kernel, gw_data, wmax_imp=None, eps_imp=N def bath_fitting(delta_wsIab, iw_mesh, Np=5): mpi.report(f"Bath fitting for hybridization with nbath/impurity orbital = {Np}") - from adapol import hybfit + try: + from adapol import hybfit + except ImportError: + raise ImportError("bath fitting requires the adapol package (https://github.com/flatironinstitute/adapol). \n" + "Please ensure that it is installed. ") + nspin, nImp = delta_wsIab.shape[1:3] error = -1 for s in np.arange(nspin): diff --git a/python/solid_dmft/gw_embedding/gw_flow.py b/python/solid_dmft/gw_embedding/gw_flow.py index 4782fcd5..46c27788 100644 --- a/python/solid_dmft/gw_embedding/gw_flow.py +++ b/python/solid_dmft/gw_embedding/gw_flow.py @@ -233,8 +233,9 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params): assert gw_params['code'] == 'aimbes', 'Only AIMBES is currently supported as gw code' # prepare output h5 archive + archive = general_params['jobname']+'/'+general_params['seedname']+'.h5' if mpi.is_master_node(): - with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar: + with HDFArchive(archive, 'a') as ar: if 'DMFT_results' not in ar: ar.create_group('DMFT_results') if 'last_iter' not in ar['DMFT_results']: @@ -336,7 +337,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params): # write block structure to h5 archive if gw_params['it_1e'] == 1 and mpi.is_master_node(): - with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar: + with HDFArchive(archive, 'a') as ar: if 'block_structure' not in ar['DMFT_input']: ar['DMFT_input']['block_structure'] = sumk.block_structure if 'deg_shells' not in ar['DMFT_input']: @@ -558,16 +559,13 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params): 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 - mpi.report('Writing iter {} results to h5 archives {}/{}.h5 and {}.'.format( - iteration, general_params['jobname'], general_params['seedname'], gw_params['h5_file'])) + mpi.report('Writing iter {} results to h5 archives {} and {}.'.format(iteration, archive, gw_params['h5_file'])) if mpi.is_master_node(): - with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar: - results_to_archive.write(ar, sumk, general_params, solver_params, solvers, - map_imp_solver, solver_type_per_imp, iteration, - False, gw_params['mu_emb'], density_mat_pre, density_mat) - + results_to_archive.write(archive, sumk, general_params, solver_params, solvers, + map_imp_solver, solver_type_per_imp, iteration, + False, gw_params['mu_emb'], density_mat_pre, density_mat) + with HDFArchive(archive, 'a') as ar: # store also IR / DLR Sigma ar['DMFT_results/it_{}'.format(iteration)]['ir_mesh'] = ir_mesh ar['DMFT_results/it_{}'.format(iteration)]['Sigma_imp_wsIab'] = Sigma_ir @@ -576,7 +574,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params): ar['DMFT_results/it_{}'.format(iteration)][f'Sigma_dlr_{ish}'] = Sigma_dlr[ish] # write results to GW h5_file - with HDFArchive(gw_params['h5_file'],'a') as ar: + with HDFArchive(gw_params['h5_file'], 'a') as ar: ar[f'downfold_1e/iter{iteration}']['Sigma_imp_wsIab'] = Sigma_ir ar[f'downfold_1e/iter{iteration}']['Vhf_imp_sIab'] = Vhf_imp_sIab