[gw] Implemented gw_sigma for DLR

YannIntVeld · YannIntVeld · commit 74ca92492feb · 2024-03-22T09:54:28.000-04:00
diff --git a/c++/triqs_tprf/lattice/gw.cpp b/c++/triqs_tprf/lattice/gw.cpp
@@ -169,36 +169,29 @@ namespace triqs_tprf {
   return sigma_r;
   }
 
-  e_k_t fock_sigma(chi_k_cvt v_k, g_wk_cvt g_wk) {
-
-  if (v_k.mesh() != std::get<1>(g_wk.mesh())) TRIQS_RUNTIME_ERROR << "fock_sigma: k-space meshes are not the same.\n";
-
-  auto _     = all_t{};
-  auto kmesh = std::get<1>(g_wk.mesh());
-
-  e_k_t sigma_k(kmesh, g_wk.target_shape());
-  sigma_k() = 0.0;
-
-  auto arr = mpi_view(kmesh);
-#pragma omp parallel for
-  for (unsigned int idx = 0; idx < arr.size(); idx++) {
-    auto &k = arr[idx];
-
-    for (auto q : kmesh) {
-
-      auto g_w  = g_wk[_, k + q];
-      auto dens = density(g_w);
+  template<typename g_t>
+  e_k_t fock_sigma_impl(chi_k_cvt v_k, g_t g_wk) {
+    auto v_r = make_gf_from_fourier(v_k);
+    auto rho_k = rho_k_from_g_wk(g_wk);
+    auto rho_r = make_gf_from_fourier(rho_k);
+    auto sigma_fock_r = fock_sigma(v_r, rho_r);
+    auto sigma_fock_k = make_gf_from_fourier(sigma_fock_r);
+    return sigma_fock_k;
+  }
 
-      for (auto [a, b, c, d] : v_k.target_indices()) { sigma_k[k](a, b) += -v_k[q](a, c, d, b) * dens(d, c) / kmesh.size(); }
-    }
+  e_k_t fock_sigma(chi_k_cvt v_k, g_wk_cvt g_wk) {
+    return fock_sigma_impl(v_k, g_wk);
   }
-  sigma_k = mpi::all_reduce(sigma_k);
-  return sigma_k;
+  e_k_t fock_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk) {
+    return fock_sigma_impl(v_k, g_wk);
   }
 
   e_k_t gw_sigma(chi_k_cvt v_k, g_wk_cvt g_wk) {
     return fock_sigma(v_k, g_wk);
   }
+  e_k_t gw_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk) {
+    return fock_sigma(v_k, g_wk);
+  }
 
   template<typename W_dyn_t, typename W_const_t, typename g_t>
   auto gw_sigma_impl(W_dyn_t W_dyn_wk, W_const_t W_const_k, g_t g_wk) {
@@ -214,29 +207,20 @@ namespace triqs_tprf {
     TRIQS_RUNTIME_ERROR << "gw_sigma: k-space meshes are not the same.\n";
 
   // Dynamic GW self energy
-  //auto g_tr = make_gf_from_fourier<0, 1>(g_wk); // Fixme! Use parallell transform
   auto g_wr = fourier_wk_to_wr(g_wk);
   auto g_tr = fourier_wr_to_tr(g_wr);
-  
-  //auto W_dyn_tr = make_gf_from_fourier<0, 1>(W_dyn_wk); // Fixme! Use parallell transform
   auto W_dyn_wr = chi_wr_from_chi_wk(W_dyn_wk);
   auto W_dyn_tr = chi_tr_from_chi_wr(W_dyn_wr);
   
   auto sigma_dyn_tr = gw_dynamic_sigma(W_dyn_tr, g_tr);
 
   // Static Fock part
-
-  //auto sigma_fock_k = fock_sigma(W_const_k, g_wk); // Has N_k^2 scaling, not fast..
-
-  auto W_const_r = make_gf_from_fourier(W_const_k);
-  auto rho_k = rho_k_from_g_wk(g_wk);
-  auto rho_r = make_gf_from_fourier(rho_k);
-  auto sigma_fock_r = fock_sigma(W_const_r, rho_r);
-  auto sigma_fock_k = make_gf_from_fourier(sigma_fock_r);
+  auto sigma_fock_k = fock_sigma(W_const_k, g_wk);
 
   // Add dynamic and static parts
   auto _ = all_t{};
-  auto sigma_wk = make_gf_from_fourier<0, 1>(sigma_dyn_tr); // Fixme! Use parallell transform
+  auto sigma_wr = fourier_tr_to_wr(sigma_dyn_tr);
+  auto sigma_wk = fourier_wr_to_wk(sigma_wr);
   for (auto w : gwm) sigma_wk[w, _] += sigma_fock_k; // Single loop with no work per w, no need to parallellize over mpi
 
   return sigma_wk;
@@ -247,11 +231,14 @@ namespace triqs_tprf {
     return gw_sigma_impl(W_dyn_wk, W_const_k, g_wk);
   }
 
-  /*
-  g_Dwk_t gw_sigma(chi_Dwk_cvt W_wk, g_Dwk_cvt g_wk) {
-    return gw_sigma_impl(W_wk, g_wk);
+  g_Dwk_t gw_sigma(chi_Dwk_cvt W_wk, chi_k_cvt v_k, g_Dwk_cvt g_wk) {
+    auto wmesh = std::get<0>(W_wk.mesh());
+    auto _ = all_t{};
+    chi_Dwk_t W_dyn_wk(W_wk.mesh(), W_wk.target_shape());
+    for (auto w : wmesh) W_dyn_wk[w, _] = W_wk[w,_] - v_k;
+
+    return gw_sigma_impl(W_dyn_wk, v_k, g_wk);
   }
-  */
   
   // ----------------------------------------------------
   // g0w_sigma via spectral representation
diff --git a/c++/triqs_tprf/lattice/gw.hpp b/c++/triqs_tprf/lattice/gw.hpp
@@ -90,7 +90,57 @@ namespace triqs_tprf {
  */
 
   g_wk_t gw_sigma(chi_wk_cvt W_wk, g_wk_cvt g_wk);
-  
+
+  /** GW self energy :math:`\Sigma(i\omega_n, \mathbf{k})` calculator for dynamic interactions
+
+    Fourier transforms the dynamic part of the interaction and the 
+    single-particle Green's function to imaginary time and real space.
+
+    .. math::
+        G_{ab}(\tau, \mathbf{r}) = \mathcal{F}^{-1}
+          \left\{ G_{ab}(i\omega_n, \mathbf{k}) \right\}
+
+    .. math::
+        W^{(dyn)}_{abcd}(\tau, \mathbf{r}) = \mathcal{F}^{-1}
+          \left\{ W^{(dyn)}_{abcd}(i\omega_n, \mathbf{k}) \right\}
+
+    computes the GW self-energy as the product
+
+    .. math::
+        \Sigma^{(dyn)}_{ab}(\tau, \mathbf{r}) =
+          - \sum_{cd} W^{(dyn)}_{acdb}(\tau, \mathbf{r}) G_{cd}(\tau, \mathbf{r})
+
+    and transforms back to frequency and momentum
+
+    .. math::
+        \Sigma^{(dyn)}_{ab}(i\omega_n, \mathbf{k}) =
+          \mathcal{F} \left\{ \Sigma^{(dyn)}_{ab}(\tau, \mathbf{r}) \right\}
+
+    The self-energy of the static part of the interaction is calculated
+    as the sum
+
+    .. math::
+        \Sigma^{(stat)}_{ab}(\mathbf{k}) = -\frac{1}{N_k}
+          \sum_{\mathbf{q},cd} V_{acdb}(\mathbf{k}) \rho_{dc}(\mathbf{k} + \mathbf{q})
+
+    where :math:`\rho_{ab}(\mathbf{k}) = -G_{ba}(\beta, \mathbf{k})` is the density matrix of the
+    single particle Green's function.
+
+    The total GW self-energy is given by
+
+    .. math::
+        \Sigma_{ab}(i\omega_n, \mathbf{k}) = 
+          \Sigma^{(dyn)}_{ab}(i\omega_n, \mathbf{k})
+          + \Sigma^{(stat)}_{ab}(\mathbf{k})
+
+    @param W_wk interaction :math:`W_{abcd}(i\omega_n, \mathbf{k})`
+    @param V_k static interaction :math:`V_{abcd}(\mathbf{q})`
+    @param g_wk single particle Green's function :math:`G_{ab}(i\omega_n, \mathbf{k})`
+    @return GW self-energy :math:`\Sigma_{ab}(i\omega_n, \mathbf{k})`
+ */
+
+  g_Dwk_t gw_sigma(chi_Dwk_cvt W_wk, chi_k_cvt v_k, g_Dwk_cvt g_wk);
+
   /** Hartree self energy :math:`\Sigma_{ab}(\mathbf{k})` calculator
 
     Computes the Hartree self-energy of a static interaction as the sum
@@ -128,6 +178,7 @@ namespace triqs_tprf {
 */
 
   e_k_t fock_sigma(chi_k_cvt v_k, g_wk_cvt g_wk);
+  e_k_t fock_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk);
 
   e_r_t fock_sigma(chi_r_cvt v_r, e_r_cvt rho_r);
 
@@ -141,6 +192,7 @@ namespace triqs_tprf {
 */
 
   e_k_t gw_sigma(chi_k_cvt v_k, g_wk_cvt g_wk);
+  e_k_t gw_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk);
   
   /** Dynamic GW self energy :math:`\Sigma(\tau, \mathbf{r})` calculator 
 
diff --git a/python/triqs_tprf/lattice_desc.py b/python/triqs_tprf/lattice_desc.py
@@ -698,6 +698,9 @@
 out
      GW self-energy :math:`\Sigma_{ab}(i\omega_n, \mathbf{k})`""")
 
+
+module.add_function ("g_Dwk_t gw_sigma(chi_Dwk_cvt W_wk, chi_k_cvt v_k, g_Dwk_cvt g_wk)")
+
 module.add_function ("triqs_tprf::e_r_t triqs_tprf::hartree_sigma (triqs_tprf::chi_k_cvt v_k, triqs_tprf::e_r_cvt rho_r)")
 module.add_function ("triqs_tprf::e_r_t triqs_tprf::fock_sigma (triqs_tprf::chi_r_cvt v_r, triqs_tprf::e_r_cvt rho_r)")
 
@@ -749,6 +752,8 @@
 out
      Fock self-energy :math:`\Sigma_{ab}(\mathbf{k})`""")
 
+module.add_function ("e_k_t fock_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk)")
+
 module.add_function ("triqs_tprf::e_k_t triqs_tprf::gw_sigma (triqs_tprf::chi_k_cvt v_k, triqs_tprf::g_wk_cvt g_wk)", doc = r"""Static GW self energy :math:`\Sigma(\mathbf{k})` calculator
 
     Computes the GW self-energy (equivalent to the Fock self-energy)
@@ -766,6 +771,8 @@
 out
      Static GW self-energy (Fock) :math:`\Sigma_{ab}(\mathbf{k})`""")
 
+module.add_function ("e_k_t gw_sigma(chi_k_cvt v_k, g_Dwk_cvt g_wk)")
+
 module.add_function ("triqs_tprf::g_tr_t triqs_tprf::gw_dynamic_sigma (triqs_tprf::chi_tr_cvt W_tr, triqs_tprf::g_tr_cvt g_tr)", doc = r"""Dynamic GW self energy :math:`\Sigma(\tau, \mathbf{r})` calculator
 
     Computes the GW self-energy as the product
diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt
@@ -11,6 +11,7 @@ set(all_tests
   lattice_utility
   gf
   gw
+  gw_dlr
   gw_hubbard_dimer
   gw_hubbard_dimer_dlr
   gw_hubbard_dimer_matrix
diff --git a/test/python/gw_dlr.py b/test/python/gw_dlr.py
@@ -0,0 +1,122 @@
+# ----------------------------------------------------------------------
+
+import itertools
+import numpy as np
+
+# ----------------------------------------------------------------------
+
+from triqs_tprf.tight_binding import TBLattice
+
+from triqs_tprf.lattice import lattice_dyson_g0_wk
+from triqs_tprf.lattice import lattice_dyson_g_wk
+from triqs_tprf.lattice import dlr_on_imfreq, dlr_on_imtime
+
+from triqs_tprf.gw import bubble_PI_wk
+from triqs_tprf.gw import dynamical_screened_interaction_W
+from triqs_tprf.gw import gw_sigma, fock_sigma
+
+from triqs.gf import Gf, Idx
+from triqs.gf import MeshImFreq, MeshImTime, MeshDLRImFreq, MeshDLRImTime
+from triqs.gf.mesh_product import MeshProduct
+
+from triqs.gf.gf_factories import make_gf_dlr
+
+# ----------------------------------------------------------------------
+
+def dlr_wk_on_imfreq_wk(g_Dwk, wmesh, g_k=None):
+    DLRwmesh, kmesh = g_Dwk.mesh.components
+    g_wk =Gf(mesh=MeshProduct(wmesh, kmesh), target_shape=g_Dwk.target_shape)
+    g_wk.data[:] = 0.0
+
+    if(g_k is None):
+        g_k = Gf(mesh=kmesh, target_shape=g_Dwk.target_shape)
+        g_k.data[:] = 0.0
+
+    for k in kmesh:
+        g_Dw = Gf(mesh=DLRwmesh, target_shape=g_Dwk.target_shape)
+        g_Dw.data[:] = g_Dwk.data[:,k.data_index,:] - g_k.data[k.data_index,:]
+        g_Dc = make_gf_dlr(g_Dw)
+        g_wk[:,k] = dlr_on_imfreq(g_Dc, wmesh)
+        g_wk.data[:,k.data_index,:] += g_k.data[k.data_index,:]
+
+    return g_wk
+
+def compare_g_Dwk_and_g_wk(g_Dwk, g_wk, g_k=None, decimal=7):
+    wmesh, kmesh = g_wk.mesh.components
+    DLRwmesh = g_Dwk.mesh.components[0]
+
+    g_ref_wk = dlr_wk_on_imfreq_wk(g_Dwk, wmesh, g_k)
+
+    np.testing.assert_array_almost_equal(g_wk.data[:], g_ref_wk.data[:], decimal=decimal)
+
+
+def test_gw_sigma_functions():
+    nk = 8
+    norb = 1
+    beta = 10.0
+    V = 1.0
+    mu = 0.1
+    dlreta = 1e-10
+    dlrwcut = 100.0
+    
+    t = -1.0 * np.eye(norb)
+    
+    t_r = TBLattice(
+        units = [(1, 0, 0)],
+        hopping = {
+            (+1,) : t,
+            (-1,) : t,
+            },
+        orbital_positions = [(0,0,0)]*norb,
+        )
+    
+    kmesh = t_r.get_kmesh(n_k=(nk, 1, 1))
+    e_k = t_r.fourier(kmesh)
+    
+    kmesh = e_k.mesh
+
+    wmesh_dlr = MeshDLRImFreq(beta, 'Fermion', dlrwcut, dlreta)
+
+    nw = int(dlrwcut * beta / (2.0 * np.pi) - 0.5)
+    wmesh = MeshImFreq(beta, 'Fermion', nw)
+
+    g0_wk = lattice_dyson_g0_wk(mu=mu, e_k=e_k, mesh=wmesh)
+    g0_Dwk = lattice_dyson_g0_wk(mu=mu, e_k=e_k, mesh=wmesh_dlr)
+    
+    V_k = Gf(mesh=kmesh, target_shape=[norb]*4)
+    V_k.data[:] = V
+    
+    print('--> pi_bubble')
+    PI_wk = bubble_PI_wk(g0_wk)
+    PI_Dwk = bubble_PI_wk(g0_Dwk)
+
+    compare_g_Dwk_and_g_wk(PI_Dwk, PI_wk)
+   
+    print('--> gw_sigma (static)')
+    sigma_k = gw_sigma(V_k, g0_wk)
+    sigma_dlr_k = gw_sigma(V_k, g0_Dwk)
+
+    np.testing.assert_array_almost_equal(sigma_k.data[:], sigma_dlr_k.data[:])
+
+    print('--> dynamical_screened_interaction_W')
+    Wr_full_wk = dynamical_screened_interaction_W(PI_wk, V_k)
+    Wr_full_Dwk = dynamical_screened_interaction_W(PI_Dwk, V_k)
+
+    compare_g_Dwk_and_g_wk(Wr_full_Dwk, Wr_full_wk, V_k, decimal=7)
+
+    print('--> gw_sigma (dynamic)')
+    sigma_wk = gw_sigma(Wr_full_wk, g0_wk)
+    sigma_Dwk = gw_sigma(Wr_full_Dwk, V_k, g0_Dwk)
+
+    compare_g_Dwk_and_g_wk(sigma_Dwk, sigma_wk, sigma_k, decimal=6)
+
+    print('--> lattice_dyson_g_wk')
+    g_wk = lattice_dyson_g_wk(mu, e_k, sigma_wk)
+    g_Dwk = lattice_dyson_g_wk(mu, e_k, sigma_Dwk)
+
+    compare_g_Dwk_and_g_wk(g_Dwk, g_wk)
+
+
+if __name__ == "__main__":
+    test_gw_sigma_functions()
+
