[fix] align ctseg interface

python/solid_dmft/dmft_tools/results_to_archive.py

@@ -111,9 +111,9 @@ def _compile_information(sum_k, general_params, solver_params, solvers, map_imp_
             if solver_params[isolvsec]['measure_pert_order']:
                 write_to_h5['pert_order_histo_imp_{}'.format(icrsh)] = solvers[icrsh].perturbation_order_histo
                 write_to_h5['avg_order_imp_{}.format(icrsh)'] = solvers[icrsh].avg_pert_order
-            if solver_params[isolvsec]['measure_nnt']:
+            if solver_params[isolvsec]['measure_nn_tau']:
                 write_to_h5['O_NN_{}'.format(icrsh)] = solvers[icrsh].triqs_solver.results.nn_tau
-            if solver_params[isolvsec]['measure_statehist']:
+            if solver_params[isolvsec]['measure_state_hist']:
                 write_to_h5['state_hist_{}'.format(icrsh)] = solvers[icrsh].state_histogram
             if solver_params[isolvsec]['crm_dyson_solver']:
                 write_to_h5['G_time_dlr_{}'.format(icrsh)] = solvers[icrsh].G_time_dlr

python/solid_dmft/dmft_tools/solver.py

@@ -683,14 +683,6 @@ def make_positive_definite(G):
         if self.solver_params['type'] == 'ctseg':
             # fill G0_freq from sum_k to solver
             self.triqs_solver.Delta_tau << self.Delta_time
-            # extract hartree shift per orbital for solver
-            hloc_0_bm = block_matrix_from_op(self.Hloc_0, self.sum_k.gf_struct_solver_list[self.icrsh])
-            chemical_potential = []
-            for block in hloc_0_bm:
-                for iorb in range(block.shape[0]):
-                    # minus sign here as the solver treats the term as chemical potential and not Hloc0
-                    chemical_potential.append(-block[iorb,iorb].real)
-            mpi.report('impurity levels:', chemical_potential)
 
             if self.general_params['h_int_type'][self.icrsh] == 'dyn_density_density':
                 mpi.report('add dynamic interaction from AIMBES')
@@ -713,17 +705,30 @@ def make_positive_definite(G):
                     Uloc_dlr_2idx_prime[coeff] = Uprime
 
                 # create full frequency objects
-                Uloc_tau_2idx = make_gf_imtime(Uloc_dlr_2idx, n_tau=self.solver_params['n_tau_k'])
-                Uloc_tau_2idx_prime = make_gf_imtime(Uloc_dlr_2idx_prime, n_tau=self.solver_params['n_tau_k'])
+                Uloc_tau_2idx = make_gf_imtime(Uloc_dlr_2idx, n_tau=self.solver_params['n_tau_bosonic'])
+                Uloc_tau_2idx_prime = make_gf_imtime(Uloc_dlr_2idx_prime, n_tau=self.solver_params['n_tau_bosonic'])
+
+                Uloc_tau_2idx_sumk = BlockGf(name_list=['up', 'down'], block_list=[Uloc_tau_2idx, Uloc_tau_2idx])
+                Uloc_tau_2idx_prime_sumk = BlockGf(name_list=['up', 'down'], block_list=[Uloc_tau_2idx_prime, Uloc_tau_2idx_prime])
+                Uloc_tau_2idx_solver = self.sum_k.block_structure.convert_gf(Uloc_tau_2idx_sumk,
+                                                                             ish_from=self.icrsh,
+                                                                             space_from='sumk',
+                                                                             space_to='solver')
+                Uloc_tau_2idx_prime_solver = self.sum_k.block_structure.convert_gf(Uloc_tau_2idx_prime_sumk,
+                                                                                   ish_from=self.icrsh,
+                                                                                   space_from='sumk',
+                                                                                   space_to='solver')
+
 
                 # fill D0_tau from Uloc_tau_2idx and Uloc_tau_2idx_prime
-                norb = Uloc_dlr.target_shape[0]
-                # same spin interaction
-                self.triqs_solver.D0_tau[0:norb, 0:norb] << Uloc_tau_2idx.real
-                self.triqs_solver.D0_tau[norb:2*norb, norb:2*norb] << Uloc_tau_2idx.real
-                # opposite spin interaction
-                self.triqs_solver.D0_tau[0:norb, norb:2*norb] << Uloc_tau_2idx_prime.real
-                self.triqs_solver.D0_tau[norb:2*norb, 0:norb] << Uloc_tau_2idx_prime.real
+                for iblock, Uloc_i in Uloc_tau_2idx_solver:
+                    for jblock, Uloc_j in Uloc_tau_2idx_solver:
+                        # same spin interaction
+                        if iblock == jblock:
+                            self.triqs_solver.D0_tau[iblock, jblock] << Uloc_tau_2idx_solver[iblock]
+                        # opposite spin interaction
+                        else:
+                            self.triqs_solver.D0_tau[iblock, jblock] << Uloc_tau_2idx_prime_solver[iblock]
 
                 # TODO: add Jerp_Iw to the solver
 
@@ -749,7 +754,7 @@ def make_positive_definite(G):
             self.triqs_solver_params['move_move_segment'] = False
             # Solve the impurity problem for icrsh shell
             # *************************************
-            self.triqs_solver.solve(h_int=self.h_int, chemical_potential=chemical_potential, **self.triqs_solver_params)
+            self.triqs_solver.solve(h_int=self.h_int, h_loc0=self.Hloc_0, **self.triqs_solver_params)
             # *************************************
 
             # call postprocessing
@@ -973,7 +978,8 @@ def _create_ctseg_solver(self):
 
         # Separately stores all params that go into solve() call of solver
         self.triqs_solver_params = {}
-        keys_to_pass = ('length_cycle', 'max_time', 'measure_statehist', 'measure_nnt')
+        keys_to_pass = ('length_cycle', 'max_time', 'measure_state_hist', 'measure_nn_tau', 'measure_G_tau',
+                        'measure_pert_order',)
         for key in keys_to_pass:
             self.triqs_solver_params[key] = self.solver_params[key]
 
@@ -982,23 +988,19 @@ def _create_ctseg_solver(self):
         # cast warmup cycles to int in case given in scientific notation
         self.triqs_solver_params['n_warmup_cycles'] = int(self.solver_params['n_warmup_cycles'])
 
-        # Rename parameters that are differentin ctseg than cthyb
-        self.triqs_solver_params['measure_gt'] = self.solver_params['measure_G_tau']
-        self.triqs_solver_params['measure_perturbation_order_histograms'] = self.solver_params['measure_pert_order']
-
         # Makes sure measure_gw is true if improved estimators are used
         if self.solver_params['improved_estimator']:
-            self.triqs_solver_params['measure_gt'] = True
-            self.triqs_solver_params['measure_ft'] = True
+            self.triqs_solver_params['measure_G_tau'] = True
+            self.triqs_solver_params['measure_F_tau'] = True
         else:
-            self.triqs_solver_params['measure_ft'] = False
+            self.triqs_solver_params['measure_F_tau'] = False
 
 
         gf_struct = self.sum_k.gf_struct_solver_list[self.icrsh]
         # Construct the triqs_solver instances
         triqs_solver = ctseg_solver(beta=self.general_params['beta'], gf_struct=gf_struct,
                         n_tau=self.general_params['n_tau'],
-                        n_tau_k=int(self.solver_params['n_tau_k']))
+                        n_tau_bosonic=int(self.solver_params['n_tau_bosonic']))
         return triqs_solver
 
     def _create_ftps_solver(self):
@@ -1384,14 +1386,58 @@ def set_Gs_from_G_l():
 
         # first print average sign
         if mpi.is_master_node():
-            print('\nAverage sign: {}'.format(self.triqs_solver.results.sign))
+            print('\nAverage sign: {}'.format(self.triqs_solver.results.average_sign))
 
         # get Delta_time from solver
         self.Delta_time << self.triqs_solver.Delta_tau
 
         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 block, norb in 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[block][iorb]
+
+        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)
@@ -1510,11 +1556,11 @@ def set_Gs_from_G_l():
             self.Sigma_freq << inverse(self.G0_freq) - inverse(self.G_freq)
 
 
-        if self.solver_params['measure_statehist']:
+        if self.solver_params['measure_state_hist']:
             self.state_histogram = self.triqs_solver.results.state_hist
 
         if self.solver_params['measure_pert_order']:
-            self.perturbation_order_histo  = self.triqs_solver.results.perturbation_order_histo_Delta
+            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():

python/solid_dmft/gw_embedding/gw_flow.py

@@ -404,8 +404,9 @@ 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':
-            density_shell[ish] = np.sum(solvers[ish].triqs_solver.results.densities)
+        if solvers[ish].solver_params['type'] == 'ctseg':
+            for block, occ_mat in solvers[ish].orbital_occupations.items():
+                density_shell[ish] += np.trace(occ_mat)
             density_tot += density_shell[ish]
             density_mat_unsym[ish] = {}
             for i, (block, norb) in enumerate(sumk.gf_struct_solver[ish].items()):

python/solid_dmft/io_tools/default.toml

@@ -119,12 +119,12 @@ legendre_fit = false
 length_cycle = "<no default>"
 max_time = -1
 measure_G_tau = true
-measure_nnt = false
+measure_nn_tau = false
 measure_pert_order = true
-measure_statehist = true
+measure_state_hist = true
 n_cycles_tot = "<no default>"
 n_l = "<none>"
-n_tau_k = 10001
+n_tau_bosonic = 10001
 n_warmup_cycles = "<no default>"
 off_diag_threshold = 0.0
 perform_tail_fit = false

test/python/svo_ctseg_dyn/dmft_config.toml

@@ -34,7 +34,7 @@ crm_dyson_solver = true
 fit_min_n = 10
 fit_max_n = 60
 fit_max_moment = 4
-n_tau_k = 10001
+n_tau_bosonic = 10001
 
 [gw]
 code = "aimbes"

test/python/svo_gw_emb_dyn/dmft_config.toml

@@ -21,10 +21,14 @@ it_2e = 1
 
 [solver]
 type = "ctseg"
-n_tau_k = 10001
-length_cycle = 100
+n_tau_bosonic = 10001
+length_cycle = 600
 n_warmup_cycles = 1e+3
 n_cycles_tot = 1e+5
 off_diag_threshold = 1e-5
 diag_delta = false
 crm_dyson_solver = true
+crm_dlr_wmax = 0.2
+crm_dlr_eps = 1e-6
+
+measure_nn_tau = true