Unstable: Add version of imtime_bubble_chi0_wk that allows for two Green's function inputs.

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<imtime>({beta, Boson, ntau}, chi_target);
+    auto g_pr_t = make_gf<imtime>(tmesh, g_target);
+    auto g_mr_t = make_gf<imtime>(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<imtime>({beta, Boson, ntau}, chi_target);
+    auto g_pr_t = make_gf<imtime>(tmesh, g_target);
+    auto g_mr_t = make_gf<imtime>(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<imtime>({beta, Boson, ntau}, chi_target);
+    auto g_pr_t = make_gf<imtime>(tmesh, g_target);
+    auto g_mr_t = make_gf<imtime>(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;
   }

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})`
 

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