diff --git a/python/solid_dmft/dmft_cycle.py b/python/solid_dmft/dmft_cycle.py index 9ae99dab..8d39d392 100755 --- a/python/solid_dmft/dmft_cycle.py +++ b/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) diff --git a/python/solid_dmft/dmft_tools/results_to_archive.py b/python/solid_dmft/dmft_tools/results_to_archive.py index 5bdc89ec..54cc83b0 100644 --- a/python/solid_dmft/dmft_tools/results_to_archive.py +++ b/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 diff --git a/python/solid_dmft/dmft_tools/solver.py b/python/solid_dmft/dmft_tools/solver.py index cfd6b769..497488fc 100755 --- a/python/solid_dmft/dmft_tools/solver.py +++ b/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') diff --git a/python/solid_dmft/gw_embedding/gw_flow.py b/python/solid_dmft/gw_embedding/gw_flow.py index afe02777..c456a55a 100644 --- a/python/solid_dmft/gw_embedding/gw_flow.py +++ b/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()) diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt index 27eb08e2..860a36fa 100644 --- a/test/python/CMakeLists.txt +++ b/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 ) diff --git a/test/python/svo_cthyb_basic_crm/dmft_config.toml b/test/python/svo_cthyb_basic_crm/dmft_config.toml index fdb2fd94..0ac8e447 100644 --- a/test/python/svo_cthyb_basic_crm/dmft_config.toml +++ b/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 diff --git a/test/python/svo_ctseg_crm/dmft_config.toml b/test/python/svo_ctseg_crm/dmft_config.toml new file mode 100644 index 00000000..98d0689c --- /dev/null +++ b/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 diff --git a/test/python/svo_ctseg_crm/inp.h5 b/test/python/svo_ctseg_crm/inp.h5 new file mode 100644 index 00000000..1bbda5dd Binary files /dev/null and b/test/python/svo_ctseg_crm/inp.h5 differ diff --git a/test/python/svo_ctseg_crm/ref.h5 b/test/python/svo_ctseg_crm/ref.h5 new file mode 100644 index 00000000..67562663 Binary files /dev/null and b/test/python/svo_ctseg_crm/ref.h5 differ diff --git a/test/python/svo_ctseg_crm/test.py b/test/python/svo_ctseg_crm/test.py new file mode 100644 index 00000000..7bbf2c9b --- /dev/null +++ b/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'])