[feat] use orbital occupations measured by cthyb / ctseg for reportin…

…g and add SVO CRM ctseg test

python/solid_dmft/dmft_cycle.py

@@ -740,10 +740,21 @@ def _dmft_step(sum_k, solvers, it, general_params,
             mpi.report('Actual time for solver: {:.2f} s'.format(timer() - start_time))
 
         # some printout of the obtained density matrices and some basic checks from the unsymmetrized solver output
-        density_shell[icrsh] = np.real(solvers[icrsh].G_freq_unsym.total_density())
-        density_tot += density_shell[icrsh]*shell_multiplicity[icrsh]
-        density_mat_unsym[icrsh] = solvers[icrsh].G_freq_unsym.density()
-        density_mat[icrsh] = solvers[icrsh].G_freq.density()
+        if ((solver_type_per_imp[icrsh] == 'cthyb' and solvers[icrsh].solver_params['measure_density_matrix']) or
+             solver_type_per_imp[icrsh] == 'ctseg' or
+            (solver_type_per_imp[icrsh] == 'hubbardI' and solvers[icrsh].solver_params['measure_density_matrix'])):
+            mpi.report('\nExtracting impurity occupations from measured density matrix.')
+            for block, occ_mat in solvers[icrsh].orbital_occupations.items():
+                density_shell[icrsh] += np.trace(occ_mat)
+            density_tot += density_shell[icrsh]*shell_multiplicity[icrsh]
+            density_mat_unsym[icrsh] = solvers[icrsh].orbital_occupations
+            density_mat[icrsh] = density_mat_unsym[icrsh].copy()
+            sum_k.symm_deg_gf(density_mat[icrsh], ish=icrsh)
+        else:
+            density_shell[icrsh] = np.real(solvers[icrsh].G_freq_unsym.total_density())
+            density_tot += density_shell[icrsh]*shell_multiplicity[icrsh]
+            density_mat_unsym[icrsh] = solvers[icrsh].G_freq_unsym.density()
+            density_mat[icrsh] = solvers[icrsh].G_freq.density()
         formatter.print_local_density(density_shell[icrsh], density_shell_pre[icrsh],
                                       density_mat_unsym[icrsh], sum_k.SO)
 

python/solid_dmft/dmft_tools/results_to_archive.py

@@ -59,6 +59,8 @@ def _compile_information(sum_k, general_params, solver_params, solvers, map_imp_
                     write_to_h5['Sigma_moments_{}'.format(icrsh)] = solvers[icrsh].Sigma_moments
                     write_to_h5['G_moments_{}'.format(icrsh)] = solvers[icrsh].G_moments
                     write_to_h5['Sigma_Hartree_{}'.format(icrsh)] = solvers[icrsh].Sigma_Hartree
+                    write_to_h5['orbital_occupations_{}'.format(icrsh)] = solvers[icrsh].orbital_occupations
+
 
         elif solver_type_per_imp[icrsh] == 'ftps':
             write_to_h5['Delta_freq_{}'.format(icrsh)] = solvers[icrsh].Delta_freq
@@ -102,6 +104,7 @@ 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
             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

python/solid_dmft/dmft_tools/solver.py

@@ -513,6 +513,7 @@ def solve(self, **kwargs):
                 self.h_loc_diagonalization = self.triqs_solver.ad
                 # get moments
                 from triqs_cthyb.tail_fit import sigma_high_frequency_moments, green_high_frequency_moments
+                from triqs_cthyb.util import orbital_occupations
                 self.Sigma_moments = sigma_high_frequency_moments(self.density_matrix,
                                                  self.h_loc_diagonalization,
                                                  self.sum_k.gf_struct_solver_list[self.icrsh],
@@ -524,6 +525,10 @@ def solve(self, **kwargs):
                                                  self.sum_k.gf_struct_solver_list[self.icrsh],
                                                  self.h_int
                                                  )
+                self.orbital_occupations = orbital_occupations(self.density_matrix,
+                                                 self.sum_k.gf_struct_solver_list[self.icrsh],
+                                                 self.h_loc_diagonalization
+                                                 )
 
             # *************************************
 
@@ -1200,6 +1205,7 @@ def set_Gs_from_G_l():
             self.Sigma_moments = self.triqs_solver.Sigma_moments
             self.Sigma_Hartree = self.triqs_solver.Sigma_Hartree
             self.G_moments = self.triqs_solver.G_moments
+            self.orbital_occupations = self.triqs_solver.orbital_occupations
 
         if self.solver_params['measure_pert_order']:
             self.perturbation_order = self.triqs_solver.perturbation_order
@@ -1410,6 +1416,13 @@ def set_Gs_from_G_l():
         self.G_time << self.triqs_solver.results.G_tau
         self.sum_k.symm_deg_gf(self.G_time, ish=self.icrsh)
 
+        # get occupation matrix
+        self.orbital_occupations = {bl: np.zeros((bl_size,bl_size)) for bl, bl_size in self.sum_k.gf_struct_solver_list[self.icrsh]}
+        for i, (block, norb) in enumerate(self.sum_k.gf_struct_solver[self.icrsh].items()):
+            self.orbital_occupations[block] = np.zeros((norb,norb))
+            for iorb in range(norb):
+                self.orbital_occupations[block][iorb, iorb] = self.triqs_solver.results.densities[i]
+
         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)
@@ -1442,6 +1455,9 @@ def set_Gs_from_G_l():
             for block in tail.keys():
                 self.Sigma_Hartree[block] = tail[block][0]
 
+            # recompute G_freq from Sigma with fitted tail
+            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']:
             self.F_freq = self.G_freq.copy()
@@ -1535,6 +1551,9 @@ def set_Gs_from_G_l():
 
                             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]
+
+                # recompute G_freq from Sigma with fitted tail
+                self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq)
                 print('\n')
 
 

python/solid_dmft/gw_embedding/gw_flow.py

@@ -418,14 +418,10 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
         mpi.report('Actual time for solver: {:.2f} s'.format(timer() - start_time))
 
         # some printout of the obtained density matrices and some basic checks from the unsymmetrized solver output
-        if solver_params[ish]['type'] == 'ctseg':
+        if solvers[ish].solver_params['type'] == 'ctseg':
             density_shell[ish] = np.sum(solvers[ish].triqs_solver.results.densities)
             density_tot += density_shell[ish]
-            density_mat_unsym[ish] = {}
-            for i, (block, norb) in enumerate(sumk.gf_struct_solver[ish].items()):
-                density_mat_unsym[ish][block] = np.zeros((norb,norb))
-                for iorb in range(norb):
-                    density_mat_unsym[ish][block][iorb, iorb] = solvers[ish].triqs_solver.results.densities[i]
+            density_mat_unsym[ish] = solvers[ish].orbital_occupations
             density_mat[ish] = density_mat_unsym[ish]
         else:
             density_shell[ish] = np.real(solvers[ish].G_freq_unsym.total_density())

test/python/CMakeLists.txt

@@ -41,6 +41,7 @@ set (integration_tests
      lno_hubbardI_mag
      svo_cthyb_basic_crm
      svo_cthyb_basic_tf
+     svo_ctseg_crm
      lco_ftps
      nio_cthyb_hartree
      )

test/python/svo_cthyb_basic_crm/dmft_config.toml

@@ -9,7 +9,7 @@ prec_mu = 0.001
 n_iw = 501
 n_tau = 10001
 
-h_int_type = "kanamori_den_den"
+h_int_type = "kanamori"
 U = 8.0
 J = 0.65
 

test/python/svo_ctseg_crm/dmft_config.toml

@@ -0,0 +1,34 @@
+[general]
+seedname = "inp"
+jobname = "out"
+enforce_off_diag = false
+block_threshold= 0.001
+mu_initial_guess = -0.027041
+
+prec_mu = 0.001
+n_iw = 501
+n_tau = 10001
+
+h_int_type = "kanamori_den_den"
+U = 8.0
+J = 0.65
+
+store_solver = true
+
+beta = 10
+
+n_iter_dmft = 1
+
+dc_type = 1
+dc = true
+dc_dmft = false
+
+
+[solver]
+type = "ctseg"
+length_cycle = 60
+n_warmup_cycles = 1e+4
+n_cycles_tot = 1e+6
+crm_dyson_solver = true
+crm_dlr_wmax = 4
+crm_dlr_eps = 1e-8

test/python/svo_ctseg_crm/test.py

@@ -0,0 +1,18 @@
+import shutil
+import triqs.utility.mpi as mpi
+import solid_dmft.main as solid
+
+import sys
+import importlib.util
+
+# try triqs_ctseg import
+ctseg = importlib.util.find_spec("triqs_ctseg") is not None
+if not ctseg:
+    mpi.report('ImportWarning: ctseg needs to be installed to run this test')
+    sys.exit()
+
+if mpi.is_master_node():
+    shutil.rmtree('out', ignore_errors=True)
+mpi.barrier()
+
+solid.main([None, 'dmft_config.toml'])