GW hardening [do not merge]

c++/triqs_tprf/lattice/chi_imfreq.cpp

@@ -676,6 +676,149 @@ chi_kwnn_t chiq_from_chi0q_and_gamma_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_
 
 // ----------------------------------------------------
 
+chi_kw_t chiq_sum_nu_from_chi0q_and_gamma_and_L_wn_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_ph_wnn, chi_nn_cvt L_wn) {
+
+  auto _ = all_t{};
+
+  mpi::communicator comm;
+
+  auto target_shape = chi0_wnk.target_shape();
+  // auto &[bmesh, fmesh, kmesh] = chi0_wnk.mesh();
+  auto &bmesh = std::get<0>(chi0_wnk.mesh());
+  auto &fmesh = std::get<1>(chi0_wnk.mesh());
+  auto &kmesh = std::get<2>(chi0_wnk.mesh());
+
+  double beta = fmesh.domain().beta;
+
+  chi_kw_t chi_kw({kmesh, bmesh}, target_shape);
+
+  auto arr = mpi_view(chi_kw.mesh()); // FIXME Use library implementation
+  std::cout << "BSE rank " << comm.rank() << " of " << comm.size() << " has "
+	    << arr.size() << " jobs." << std::endl;
+
+  triqs::utility::timer t;
+  t.start();
+
+#pragma omp parallel for
+  for (unsigned int idx = 0; idx < arr.size(); idx++) {
+    auto &[k, w] = arr(idx);
+
+    //triqs::utility::timer t_copy_1, t_chi0_nn, t_bse, t_chi_tr, t_copy_2;
+
+    // ----------------------------------------------------
+    //t_copy_1.start();
+
+    chi_w_t chi0_n(fmesh, target_shape);
+
+#pragma omp critical
+    chi0_n = chi0_wnk[w, _, k];
+
+    //t_copy_1.stop();
+    //std::cout << "BSE: copy_1 " << double(t_copy_1) << " s" << std::endl;
+    // ----------------------------------------------------
+    //t_chi0_nn.start();
+
+    chi_nn_t chi0_nn({fmesh, fmesh}, target_shape);
+    chi0_nn *= 0.;
+
+    for (auto const &n : fmesh) chi0_nn[n, n] = chi0_n[n];
+
+    //t_chi0_nn.stop();
+    //std::cout << "BSE: chi0_nn " << double(t_chi0_nn) << " s" << std::endl;
+    // ----------------------------------------------------
+    //t_bse.start();
+
+    auto I = identity<Channel_t::PH>(chi0_nn);
+
+    chi_nn_t gamma_nn({fmesh, fmesh}, target_shape);
+
+#pragma omp critical
+    gamma_nn = gamma_ph_wnn[w, _, _];
+
+    // this step could be optimized, using the diagonality of chi0 and I
+    auto denom = chi_nn_t{I - product<Channel_t::PH>(chi0_nn, gamma_nn)};
+
+    // also the last product here
+    auto chi = chi_nn_t{product<Channel_t::PH>(inverse<Channel_t::PH>(denom), chi0_nn)};
+
+    //t_bse.stop();
+    //std::cout << "BSE: bse inv " << double(t_bse) << " s" << std::endl;
+    // ----------------------------------------------------
+    //t_chi_tr.start();
+
+    // trace out fermionic frequencies
+    array<std::complex<double>, 4> tr_chi(target_shape);
+    tr_chi = scalar_product_PH(L_wn[w, _], chi, L_wn[w, _]);
+
+    /*
+    tr_chi() = 0;
+
+    for (auto const &n1 : fmesh) {
+      auto Ll = L_wn[w, n1];
+      for (auto const &n2 : fmesh) {
+	auto Lr = L_wn[w, n2];
+	auto C = chi[n1, n2];
+	for (auto const &[a, b, c, d] : chi.target_indices()) {
+	  for (auto const &[e, f, g, h] : chi.target_indices()) {
+	    tr_chi(a,b,c,d) += Ll(e,f,b,a) * C(e,f,g,h) * Lr(h,g,c,d);
+	  }
+	}
+      }
+    }
+
+    array<std::complex<double>, 4> tr_chi_ref(target_shape);
+    tr_chi_ref = scalar_product_PH(L_wn[w, _], chi, L_wn[w, _]);
+
+    for (auto const &[a, b, c, d] : chi.target_indices()) {
+      auto diff = std::abs( tr_chi(a,b,c,d) - tr_chi_ref(a,b,c,d) );
+      std::cout << a << ", " << b << ", " << c << ", " << d << " - " << diff << "\n";
+      if(diff > 1e-6) {
+	exit(0);
+      }	
+    }
+    */
+
+    //t_chi_tr.stop();
+    //std::cout << "BSE: chi tr " << double(t_chi_tr) << " s" << std::endl;
+    // ----------------------------------------------------
+    //t_copy_2.start();
+
+#pragma omp critical
+    chi_kw[k, w] = tr_chi;
+
+    //t_copy_2.stop();
+    //std::cout << "BSE: copy_2 " << double(t_copy_2) << " s" << std::endl;
+    // ----------------------------------------------------
+
+    if(comm.rank() == 0 && omp_get_thread_num() == 0) {
+
+      int Nomp = omp_get_num_threads();
+      int N = int(floor(arr.size() / double(Nomp)));
+
+      //double t_left = double(t) * ( N / (idx + 1) - 1. );
+      int done_percent = (N == 0) ? 100 : int(floor(100 * double(idx + 1) / N));
+
+      std::cout << "BSE " << triqs::utility::timestamp() << " "
+		<< std::setfill(' ') << std::setw(3) << done_percent << "% "
+		<< "ETA " << triqs::utility::estimate_time_left(N, idx, t)
+		<< " job no "
+		<< std::setfill(' ') << std::setw(5) << int(idx)
+		<< " of approx " << N << " jobs/thread/rank." << std::endl;
+    }
+
+  }
+
+  chi_kw = mpi::all_reduce(chi_kw);
+
+  t.stop();
+  if(comm.rank() == 0 )
+    std::cout << "BSE TIME: " << double(t) << " s" << std::endl;
+
+  return chi_kw;
+}
+
+// ----------------------------------------------------
+
 chi_kw_t chiq_sum_nu_from_chi0q_and_gamma_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_ph_wnn) {
 
   auto _ = all_t{};

c++/triqs_tprf/lattice/chi_imfreq.hpp

@@ -171,6 +171,15 @@ chi_kwnn_t chiq_from_chi0q_and_gamma_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_
  */
 chi_kw_t chiq_sum_nu_from_chi0q_and_gamma_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_ph_wnn);
 
+/** Dual lattice Bethe-Salpeter equation solver for the generalized susceptibility :math:`\chi^{(0)}_{\bar{a}b\bar{c}d}(\omega, \mathbf{k})`.
+
+  @param chi0_wnk Generalized lattice bubble susceptibility :math:`\chi^{(0)}_{\bar{a}b\bar{c}d}(\omega, \mathbf{k})`.
+  @param gamma_ph_wnn Local particle-hole vertex function :math:`\Gamma^{(PH)}_{\bar{a}b\bar{c}d}(\omega, \nu, \nu')`.
+  @param L_wn Local triangular particle-hole vertex function :math:`L^{(PH)}_{\bar{a}b\bar{c}d}(\omega, \nu)`.
+  @return Generalized lattice susceptibility :math:`\chi_{\bar{a}b\bar{c}d}(\omega, \mathbf{k})`.
+ */
+chi_kw_t chiq_sum_nu_from_chi0q_and_gamma_and_L_wn_PH(chi_wnk_cvt chi0_wnk, chi_wnn_cvt gamma_ph_wnn, chi_nn_cvt L_wn);
+
 gf<prod<brzone, imfreq>, tensor_valued<4>>
 chiq_sum_nu_from_g_wk_and_gamma_PH(gk_iw_t g_wk, g2_iw_vt gamma_ph_wnn, int tail_corr_nwf=-1);
 

c++/triqs_tprf/lattice/chi_retime.cpp

@@ -94,7 +94,6 @@ std::tuple<g_Tk_t, g_Tk_t> g0_Tk_les_gtr_from_e_k(e_k_cvt e_k, mesh::retime Tmes
 
 chi_Tr_t chi0_Tr_from_g_Tr_PH(g_Tr_cvt g_Tr_les, g_Tr_cvt g_Tr_gtr) {
 
-  auto _ = all_t{};
   auto I = std::complex(0.,1.);
 
   auto Tmesh = std::get<0>(g_Tr_les.mesh());
@@ -104,9 +103,6 @@ chi_Tr_t chi0_Tr_from_g_Tr_PH(g_Tr_cvt g_Tr_les, g_Tr_cvt g_Tr_gtr) {
 
   chi_Tr_t chi0_Tr{{Tmesh, rmesh}, {nb, nb, nb, nb}};
 
-  auto g_target = g_Tr_les.target();
-  auto chi_target = chi0_Tr.target();
-
   //for (auto const &r : rmesh) {
 
   auto arr = mpi_view(Tmesh);

c++/triqs_tprf/lattice/dynamical_screened_interaction.cpp

@@ -30,7 +30,7 @@
 namespace triqs_tprf {
 
   enum SusceptibilityType { bubble, generalized };
-
+ 
   template <SusceptibilityType susType, typename chi_t, typename v_t> auto screened_interaction_from_generic_susceptibility(chi_t &chi, v_t &V) {
 
     auto const &[freqmesh, kmesh] = chi.mesh();
@@ -82,10 +82,58 @@ namespace triqs_tprf {
     return W;
   }
 
+  chi_wk_t dynamical_screened_interaction_W(chi_wk_cvt chi_wk, chi_k_cvt V_k) {
+
+    auto const &[wmesh, kmesh] = chi_wk.mesh();
+    auto const &v_kmesh = V_k.mesh();
+
+    if (kmesh != v_kmesh) TRIQS_RUNTIME_ERROR << "dynamical_screened_interaction_W: k-space meshes are not the same\n";
+
+    auto W_wk    = make_gf(chi_wk);
+    W_wk()       = 0.;
+    size_t nb = chi_wk.target_shape()[0];
+
+    using scalar_t = typename chi_wk_cvt::scalar_t;
+    auto I         = nda::eye<scalar_t>(nb * nb);
+
+    // MPI and openMP parallell loop
+    auto arr = mpi_view(W_wk.mesh());
+
+#pragma omp parallel for 
+    for (unsigned int idx = 0; idx < arr.size(); idx++) {
+      auto &[w, k] = arr(idx);
+
+      array<scalar_t, 4> denom{nb, nb, nb, nb};
+      array<scalar_t, 4> inv_denom{nb, nb, nb, nb};
+
+      auto denom_mat = make_matrix_view(group_indices_view(denom, idx_group<0, 1>, idx_group<3, 2>));
+      auto inv_denom_mat = make_matrix_view(group_indices_view(inv_denom, idx_group<0, 1>, idx_group<3, 2>));
+
+      denom_mat = I;
+
+      for (auto const &[a, b, c, d] : V_k.target_indices())
+	for( auto e : range(nb) ) 
+	  for( auto f : range(nb) ) 
+	    denom(a,b,c,d) -= chi_wk[w,k](a,b,e,f)*V_k[k](e,f,c,d);
+
+      inv_denom_mat = inverse(denom_mat);
+
+      for (auto const &[a, b, c, d] : V_k.target_indices())
+	for( auto e : range(nb) ) 
+	  for( auto f : range(nb) ) 
+	    W_wk[w, k](a,b,c,d) += V_k[k](a,b,e,f) * inv_denom(e,f,c,d);      
+    }
+
+    W_wk = mpi::all_reduce(W_wk);
+    return W_wk;
+  }
+
+  /*
   chi_wk_t dynamical_screened_interaction_W(chi_wk_cvt PI_wk, chi_k_cvt V_k) {
     return screened_interaction_from_generic_susceptibility<bubble>(PI_wk, V_k);
   }
-
+  */
+
   chi_fk_t dynamical_screened_interaction_W(chi_fk_cvt PI_fk, chi_k_cvt V_k) {
     return screened_interaction_from_generic_susceptibility<bubble>(PI_fk, V_k);
   }