1. fix a bug when both improved_estimator and tail_fit are activated;…

… 2. fix missing conflicts with the latest unstable branch; 3. check adapol installation
TRIQS · Aug 30, 2024 · 9570cab · 9570cab
1 parent f081638
commit 9570cab
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 44 deletions.
diff --git 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
@@ -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
@@ -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