From 9a9570bf4543fb88c19ebc9d0e2ab2a4672c54a0 Mon Sep 17 00:00:00 2001 From: Andrew Hardy Date: Wed, 7 Aug 2024 11:26:29 -0400 Subject: [PATCH 1/2] header --- c++/triqs_tprf/lattice/chi_imtime.cpp | 139 ++++++++++++++++++++++++++ c++/triqs_tprf/lattice/chi_imtime.hpp | 3 +- 2 files changed, 141 insertions(+), 1 deletion(-) diff --git a/c++/triqs_tprf/lattice/chi_imtime.cpp b/c++/triqs_tprf/lattice/chi_imtime.cpp index 725c6674..9e8547c6 100644 --- a/c++/triqs_tprf/lattice/chi_imtime.cpp +++ b/c++/triqs_tprf/lattice/chi_imtime.cpp @@ -134,7 +134,58 @@ chi_tr_t chi0_tr_from_grt_PH(g_tr_cvt g_tr) { return chi0_tr; } +// ---------------------------------------------------- +// chi0 bubble in imaginary time for multiple Green's functions + +chi_tr_t chi0_tr_from_grt_ab_PH(g_tr_cvt ga_tr,g_tr_cvt gb_tr) { + + auto _ = all_t{}; + + auto tmesh = std::get<0>(ga_tr.mesh()); + auto rmesh = std::get<1>(ga_tr.mesh()); + + int nb = ga_tr.target().shape()[0]; + int ntau = tmesh.size(); + double beta = tmesh.beta(); + + chi_tr_t chi0_tr{{{beta, Boson, ntau}, rmesh}, {nb, nb, nb, nb}}; + + auto g_target = ga_tr.target(); + auto chi_target = chi0_tr.target(); + + // -- This does not work on the boundaries!! The eval wraps to the other + // regime! + // -- gt(beta) == gt(beta + 0^+) + // chi0_tr(tau, r)(a, b, c, d) << g_tr(tau, r)(d, a) * g_tr(-tau, -r)(b, c); + + //for (auto r : rmesh) { + auto arr = mpi_view(rmesh); + +#pragma omp parallel for + for (unsigned int idx = 0; idx < arr.size(); idx++) { + auto &r = arr[idx]; + + auto chi0_t = make_gf({beta, Boson, ntau}, chi_target); + auto g_pr_t = make_gf(tmesh, g_target); + auto g_mr_t = make_gf(tmesh, g_target); + +#pragma omp critical + { + g_pr_t = ga_tr[_, r]; + g_mr_t = gb_tr(_, -r); + } + + for (auto t : tmesh) chi0_t[t](a, b, c, d) << g_pr_t(t)(d, a) * g_mr_t(beta - t)(b, c); + +#pragma omp critical + chi0_tr[_, r] = chi0_t; + } + + chi0_tr = mpi::all_reduce(chi0_tr); + + return chi0_tr; +} // -- memory optimized version for smaller nw chi_wr_t chi0_wr_from_grt_PH(g_tr_cvt g_tr, int nw=1) { @@ -180,7 +231,51 @@ chi_wr_t chi0_wr_from_grt_PH(g_tr_cvt g_tr, int nw=1) { chi0_wr = mpi::all_reduce(chi0_wr); return chi0_wr; } +// -- memory optimized version for smaller nw +chi_wr_t chi0_wr_from_grt_ab_PH(g_tr_cvt ga_tr, g_tr_cvt gb_tr,int nw=1) { + auto _ = all_t{}; + + auto tmesh = std::get<0>(ga_tr.mesh()); + auto rmesh = std::get<1>(ga_tr.mesh()); + + int nb = ga_tr.target().shape()[0]; + int ntau = tmesh.size(); + double beta = tmesh.beta(); + + chi_wr_t chi0_wr{{{beta, Boson, nw}, rmesh}, {nb, nb, nb, nb}}; + + auto g_target = ga_tr.target(); + auto chi_target = chi0_wr.target(); + + auto arr = mpi_view(rmesh); + +#pragma omp parallel for + for (unsigned int idx = 0; idx < arr.size(); idx++) { + auto &r = arr[idx]; + + auto chi0_t = make_gf({beta, Boson, ntau}, chi_target); + auto g_pr_t = make_gf(tmesh, g_target); + auto g_mr_t = make_gf(tmesh, g_target); + +#pragma omp critical + { + g_pr_t = ga_tr[_, r]; + g_mr_t = gb_tr(_, -r); + } + + for (auto t : tmesh) chi0_t[t](a, b, c, d) << g_pr_t(t)(d, a) * g_mr_t(beta - t)(b, c); + +#pragma omp critical + { + auto chi0_w = make_gf_from_fourier(chi0_t, nw); + chi0_wr[_, r] = chi0_w; + } + } + + chi0_wr = mpi::all_reduce(chi0_wr); + return chi0_wr; +} // -- optimized version for w=0 chi_wr_t chi0_w0r_from_grt_PH(g_tr_cvt g_tr) { @@ -219,6 +314,50 @@ chi_wr_t chi0_w0r_from_grt_PH(g_tr_cvt g_tr) { auto int_chi0 = chi_trapz_tau(chi0_t); +#pragma omp critical + chi0_wr[0, r] = int_chi0; + } + + chi0_wr = mpi::all_reduce(chi0_wr); + return chi0_wr; +} +chi_wr_t chi0_w0r_from_grt_ab_PH(g_tr_cvt ga_tr,g_tr_cvt gb_tr) { + + auto _ = all_t{}; + + auto tmesh = std::get<0>(ga_tr.mesh()); + auto rmesh = std::get<1>(ga_tr.mesh()); + + int nw = 1; + int nb = ga_tr.target().shape()[0]; + int ntau = tmesh.size(); + double beta = tmesh.beta(); + + chi_wr_t chi0_wr{{{beta, Boson, nw}, rmesh}, {nb, nb, nb, nb}}; + + auto g_target = ga_tr.target(); + auto chi_target = chi0_wr.target(); + + auto arr = mpi_view(rmesh); + +#pragma omp parallel for + for (unsigned int idx = 0; idx < arr.size(); idx++) { + auto &r = arr[idx]; + + auto chi0_t = make_gf({beta, Boson, ntau}, chi_target); + auto g_pr_t = make_gf(tmesh, g_target); + auto g_mr_t = make_gf(tmesh, g_target); + +#pragma omp critical + { + g_pr_t = ga_tr[_, r]; + g_mr_t = gb_tr(_, -r); + } + + for (auto t : tmesh) chi0_t[t](a, b, c, d) << g_pr_t(t)(d, a) * g_mr_t(beta - t)(b, c); + + auto int_chi0 = chi_trapz_tau(chi0_t); + #pragma omp critical chi0_wr[0, r] = int_chi0; } diff --git a/c++/triqs_tprf/lattice/chi_imtime.hpp b/c++/triqs_tprf/lattice/chi_imtime.hpp index 03bad53c..84097a96 100644 --- a/c++/triqs_tprf/lattice/chi_imtime.hpp +++ b/c++/triqs_tprf/lattice/chi_imtime.hpp @@ -38,7 +38,7 @@ namespace triqs_tprf { chi_tr_t chi0_tr_from_grt_PH(g_tr_cvt g_tr); chi_Dtr_t chi0_tr_from_grt_PH(g_Dtr_cvt g_tr, bool symmetrize=false); chi_wr_t chi0_wr_from_grt_PH(g_tr_cvt g_tr, int nw); - +chi_wr_t chi0_wr_from_grt_ab_PH(g_tr_cvt ga_tr, g_tr_cvt gb_tr, int nw); /** Generalized susceptibility zero imaginary frequency bubble in the particle-hole channel :math:`\chi^{(0)}_{\bar{a}b\bar{c}d}(\omega=0, \mathbf{r})` Computes @@ -52,6 +52,7 @@ chi_wr_t chi0_wr_from_grt_PH(g_tr_cvt g_tr, int nw); @return Generalized susceptibility :math:`\chi^{(0)}_{\bar{a}b\bar{c}d}(\mathbf{r})` in real-space. */ chi_wr_t chi0_w0r_from_grt_PH(g_tr_cvt g_tr); +chi_wr_t chi0_w0r_from_grt_ab_PH(g_tr_cvt ga_tr,g_tr_cvt gb_tr); /** Static susceptibility calculation :math:`\chi_{\bar{a}b\bar{c}d}(\omega=0, \mathbf{r})` From 7b8d294f13a3231fbf10002f925b45fe30164af7 Mon Sep 17 00:00:00 2001 From: Andrew Hardy Date: Wed, 7 Aug 2024 11:28:59 -0400 Subject: [PATCH 2/2] lattice_utils update --- python/triqs_tprf/lattice_utils.py | 91 ++++++++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/python/triqs_tprf/lattice_utils.py b/python/triqs_tprf/lattice_utils.py index 6af29aee..7dfd1ef7 100644 --- a/python/triqs_tprf/lattice_utils.py +++ b/python/triqs_tprf/lattice_utils.py @@ -275,7 +275,98 @@ def imtime_bubble_chi0_wk(g_wk, nw=1, save_memory=False, verbose=True): del chi0_wr return chi0_wk +# ---------------------------------------------------------------------- +def imtime_bubble_chi0_ab_wk(ga_wk,gb_wk, nw=1, save_memory=False, verbose=True): + ncores = multiprocessing.cpu_count() + + wmesh, kmesh = g_wk.mesh.components + + norb = g_wk.target_shape[0] + beta = wmesh.beta + nw_g = len(wmesh) + nk = len(kmesh) + + ntau = 2 * nw_g + + is_dlr_mesh = type(wmesh) == MeshDLRImFreq + if is_dlr_mesh: + iw_values = np.fromiter(wmesh.values(), dtype=complex) + is_symmetrized = np.allclose(iw_values, -iw_values[::-1]) + + + # -- Memory Approximation + + ng_tr = ntau * np.prod(nk) * norb**2 # storing G(tau, r) + ng_wr = nw_g * np.prod(nk) * norb**2 # storing G(w, r) + ng_t = ntau * norb**2 # storing G(tau) + + nchi_tr = ntau * np.prod(nk) * norb**4 # storing \chi(tau, r) + nchi_wr = nw * np.prod(nk) * norb**4 # storing \chi(w, r) + nchi_t = ntau * norb**4 # storing \chi(tau) + nchi_w = nw * norb**4 # storing \chi(w) + nchi_r = np.prod(nk) * norb**4 # storing \chi(r) + + if nw == 1: + ntot_case_1 = ng_tr + ng_wr + ntot_case_2 = ng_tr + nchi_wr + ncores*(nchi_t + 2*ng_t) + ntot_case_3 = 4 * nchi_wr + + ntot = max(ntot_case_1, ntot_case_2, ntot_case_3) + + else: + ntot_case_1 = ng_tr + nchi_tr + ncores*(nchi_t + 2*ng_t) + ntot_case_2 = nchi_tr + nchi_wr + ncores*(nchi_w + nchi_t) + + ntot = max(ntot_case_1, ntot_case_2) + + nbytes = ntot * np.complex128().nbytes + ngb = nbytes / 1024.**3 + + if verbose and mpi.is_master_node(): + print(tprf_banner(), "\n") + print('beta =', beta) + print('nk =', nk) + print('nw =', nw_g) + print('norb =', norb) + print() + print('Approx. Memory Utilization: %2.2f GB\n' % ngb) + + if verbose: mpi.report('--> fourier_wk_to_wr') + ga_wr = fourier_wk_to_wr(ga_wk) + gb_wr = fourier_wk_to_wr(gb_wk) + del ga_wk + del gb_wk + + if verbose: mpi.report('--> fourier_wr_to_tr') + ga_tr = fourier_wr_to_tr(ga_wr) + gb_tr = fourier_wr_to_tr(gb_wr) + del ga_wr + del gb_wr + + if nw == 1: + if verbose: mpi.report('--> chi0_w0r_from_grt_PH (bubble in tau & r)') + chi0_wr = chi0_w0r_from_grt_ab_PH(ga_tr, gb_tr) + del ga_tr + del gb_tr + else: + if not save_memory: + if verbose: mpi.report('--> chi0_tr_from_grt_PH (bubble in tau & r)') + chi0_tr = chi0_tr_from_grt_ab_PH(ga_tr, gb_tr) + del ga_tr + del gb_tr + + if verbose: mpi.report('--> chi_wr_from_chi_tr') + chi0_wr = chi_wr_from_chi_tr(chi0_tr, nw=nw) + del chi0_tr + elif save_memory: + chi0_wr = chi0_wr_from_grt_PH(g_tr, nw=nw) + + if verbose: mpi.report('--> chi_wk_from_chi_wr (r->k)') + chi0_wk = chi_wk_from_chi_wr(chi0_wr) + del chi0_wr + + return chi0_wk # ---------------------------------------------------------------------- def chi_contraction(chi, op1, op2): """Contract a susceptibility with two operators