From a39818f7f1508e5cee61882e147c7076c524e5e0 Mon Sep 17 00:00:00 2001 From: Yuichi Motoyama Date: Tue, 10 Sep 2024 15:56:33 +0900 Subject: [PATCH] make dcore_anacont and spectrum work with TRIQS --- src/dcore/anacont/spm.py | 10 ++++++---- src/dcore/dcore_anacont.py | 24 +++++++++++++----------- src/dcore/dcore_spectrum.py | 7 ++++++- 3 files changed, 25 insertions(+), 16 deletions(-) diff --git a/src/dcore/anacont/spm.py b/src/dcore/anacont/spm.py index 431ee6ac..86557f15 100644 --- a/src/dcore/anacont/spm.py +++ b/src/dcore/anacont/spm.py @@ -334,9 +334,11 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver): mesh_iw = MeshImFreq(beta, "Fermion", n_matsubara_retain) gf_iw = GfImFreq(data=data, beta=beta, mesh=mesh_iw) n_orbitals = gf_iw.data.shape[1] - matsubara_frequencies = np.imag(gf_iw.mesh.values()[n_matsubara_retain:]) + iws = list(gf_iw.mesh.values()) + ws = list(mesh_w.values()) + matsubara_frequencies = np.imag(iws[n_matsubara_retain:]) sigma_w_data = np.zeros( - (mesh_w.size, n_orbitals, n_orbitals), dtype=np.complex128 + (len(ws), n_orbitals, n_orbitals), dtype=np.complex128 ) for i_orb in range(n_orbitals): print( @@ -351,7 +353,7 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver): ) density, gf_tau_fit, energies_extract, sum_rule_const, chi2 = get_single_continuation( - tau_grid, gf_tau, n_sv, beta, mesh_w.values()[0], mesh_w.values()[-1], mesh_w.size, sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, solver=solver, solver_opts=params_spm_solver + tau_grid, gf_tau, n_sv, beta, ws[0], ws[-1], len(ws), sum_rule_const=const_imag_tail, lambd=L1_coeff, verbose=verbose_opt, solver=solver, solver_opts=params_spm_solver ) energies, gf_real, gf_imag = get_kramers_kronig_realpart( energies_extract, dos_to_gf_imag(density) @@ -373,4 +375,4 @@ def anacont(sigma_iw_npz, beta, mesh_w, params_spm, params_spm_solver): plt.close() sigma_w = GfReFreq(mesh=mesh_w, data=sigma_w_data) data_w[key] = sigma_w.data - return data_w \ No newline at end of file + return data_w diff --git a/src/dcore/dcore_anacont.py b/src/dcore/dcore_anacont.py index 4104ab7b..8aef0ab2 100644 --- a/src/dcore/dcore_anacont.py +++ b/src/dcore/dcore_anacont.py @@ -61,9 +61,17 @@ def dcore_anacont(inifile): params_ac["post.anacont.spm.solver"] = {} data_w = spm.anacont(sigma_iw_npz, beta, mesh_w, params_ac["post.anacont.spm"], params_ac["post.anacont.spm.solver"]) else: - assert False, "Unknown solver: " + solver + sys.exit("ERROR: Unknown anacont solver " + solver) + omega = numpy.array(list(mesh_w.values())) + data_w["omega"] = omega + + dir_post = params["post"]["dir_post"] + if not os.path.exists(dir_post): + os.makedirs(dir_post) + file_sigma_w = os.path.join(dir_post, "sigma_w.npz") + print("Writing to", file_sigma_w + "...") + numpy.savez(file_sigma_w, **data_w) - data_w["omega"] = mesh_w.values() if params["post.anacont"]["show_result"] or params["post.anacont"]["save_result"]: import matplotlib.pyplot as plt @@ -75,8 +83,8 @@ def dcore_anacont(inifile): sigma_w = data_w[f"data{idata}"] for iorb, jorb in itertools.product(range(sigma_w.shape[1]), range(sigma_w.shape[2])): plt.axhline(y=0, xmin=omega_min, xmax=omega_max, color="lightgrey") - plt.plot(mesh_w.points, sigma_w[:,iorb,jorb].real, label="Real") - plt.plot(mesh_w.points, sigma_w[:,iorb,jorb].imag, label="Imag") + plt.plot(omega, sigma_w[:,iorb,jorb].real, label="Real") + plt.plot(omega, sigma_w[:,iorb,jorb].imag, label="Imag") plt.xlim(omega_min, omega_max) plt.xlabel(r"$\omega$") plt.legend() @@ -90,12 +98,6 @@ def dcore_anacont(inifile): plt.show() plt.close() - dir_post = params["post"]["dir_post"] - if not os.path.exists(dir_post): - os.makedirs(dir_post) - file_sigma_w = os.path.join(dir_post, "sigma_w.npz") - print("Writing to", file_sigma_w + "...") - numpy.savez(file_sigma_w, **data_w) def run(): @@ -120,4 +122,4 @@ def run(): args = parser.parse_args() - dcore_anacont(args.path_input_files) \ No newline at end of file + dcore_anacont(args.path_input_files) diff --git a/src/dcore/dcore_spectrum.py b/src/dcore/dcore_spectrum.py index ccf74a0e..e280e55c 100644 --- a/src/dcore/dcore_spectrum.py +++ b/src/dcore/dcore_spectrum.py @@ -167,7 +167,12 @@ def __init__(self, seedname, params, n_k, xk, nkdiv_mesh, kvec_mesh, dir_post): self._sigma_w_file = os.path.join(dir_post, "sigma_w.npz") ws = _read_ws(self._sigma_w_file) self.mesh_w = MeshReFreq(ws[0], ws[-1], len(ws)) - self.mesh_w._points = ws + if TRIQS_COMPAT: + self.mesh_w._points = ws + else: + ws_linear = np.linspace(ws[0], ws[-1], num=len(ws)) + if not np.allclose(ws, ws_linear): + sys.exit("Error: omega points must be linearly spaced when using TRIQS (DCORE_TRIQS_COMPAT=0)") # Construct a SumKDFT object self._omega_min = ws[0]