From 602626708f3053083fadd76bde4ce9a1b04dcb0a Mon Sep 17 00:00:00 2001 From: Alexander Hampel Date: Mon, 8 Jul 2024 15:50:13 -0400 Subject: [PATCH] [fix] align ctseg interface --- .../dmft_tools/results_to_archive.py | 4 +- python/solid_dmft/dmft_tools/solver.py | 66 ++++++++++--------- python/solid_dmft/gw_embedding/gw_flow.py | 3 +- python/solid_dmft/io_tools/default.toml | 6 +- test/python/svo_ctseg_dyn/dmft_config.toml | 2 +- test/python/svo_gw_emb_dyn/dmft_config.toml | 4 +- 6 files changed, 44 insertions(+), 41 deletions(-) diff --git a/python/solid_dmft/dmft_tools/results_to_archive.py b/python/solid_dmft/dmft_tools/results_to_archive.py index 9b36c73b..8d3b1234 100644 --- a/python/solid_dmft/dmft_tools/results_to_archive.py +++ b/python/solid_dmft/dmft_tools/results_to_archive.py @@ -116,9 +116,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 diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py index 80f6647d..8649d434 100755 --- a/python/solid_dmft/dmft_tools/solver.py +++ b/python/solid_dmft/dmft_tools/solver.py @@ -698,14 +698,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') @@ -728,17 +720,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 @@ -764,7 +769,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 @@ -988,7 +993,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] @@ -997,23 +1003,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): @@ -1410,7 +1412,7 @@ 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 @@ -1420,10 +1422,10 @@ def set_Gs_from_G_l(): # 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()): + 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[i] + 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 = {} @@ -1591,11 +1593,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(): diff --git a/python/solid_dmft/gw_embedding/gw_flow.py b/python/solid_dmft/gw_embedding/gw_flow.py index c456a55a..9be85944 100644 --- a/python/solid_dmft/gw_embedding/gw_flow.py +++ b/python/solid_dmft/gw_embedding/gw_flow.py @@ -419,7 +419,8 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params): # some printout of the obtained density matrices and some basic checks from the unsymmetrized solver output if solvers[ish].solver_params['type'] == 'ctseg': - density_shell[ish] = np.sum(solvers[ish].triqs_solver.results.densities) + 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] = solvers[ish].orbital_occupations density_mat[ish] = density_mat_unsym[ish] diff --git a/python/solid_dmft/io_tools/default.toml b/python/solid_dmft/io_tools/default.toml index d74d498b..428ba0ab 100644 --- a/python/solid_dmft/io_tools/default.toml +++ b/python/solid_dmft/io_tools/default.toml @@ -123,12 +123,12 @@ legendre_fit = false length_cycle = "" 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 = "" n_l = "" -n_tau_k = 10001 +n_tau_bosonic = 10001 n_warmup_cycles = "" off_diag_threshold = 0.0 perform_tail_fit = false diff --git a/test/python/svo_ctseg_dyn/dmft_config.toml b/test/python/svo_ctseg_dyn/dmft_config.toml index 590d0a8c..e4f683cd 100644 --- a/test/python/svo_ctseg_dyn/dmft_config.toml +++ b/test/python/svo_ctseg_dyn/dmft_config.toml @@ -33,7 +33,7 @@ off_diag_threshold = 1e-4 crm_dyson_solver = true crm_dlr_wmax = 2 crm_dlr_eps = 1e-6 -n_tau_k = 10001 +n_tau_bosonic = 10001 [gw] code = "aimbes" diff --git a/test/python/svo_gw_emb_dyn/dmft_config.toml b/test/python/svo_gw_emb_dyn/dmft_config.toml index e2ae5fcf..3ab3f82f 100644 --- a/test/python/svo_gw_emb_dyn/dmft_config.toml +++ b/test/python/svo_gw_emb_dyn/dmft_config.toml @@ -21,7 +21,7 @@ it_2e = 1 [solver] type = "ctseg" -n_tau_k = 10001 +n_tau_bosonic = 10001 length_cycle = 600 n_warmup_cycles = 1e+3 n_cycles_tot = 1e+5 @@ -31,4 +31,4 @@ crm_dyson_solver = true crm_dlr_wmax = 0.2 crm_dlr_eps = 1e-6 -measure_nnt = true +measure_nn_tau = true