[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

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

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,52 +1551,39 @@ 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')
 
 
             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:

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]

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