From 74ca92492feb3108eed7f322a60b49c65a6ae2f1 Mon Sep 17 00:00:00 2001 From: Yann in 't Veld Date: Fri, 22 Mar 2024 09:54:28 -0400 Subject: [PATCH] [gw] Implemented gw_sigma for DLR --- c++/triqs_tprf/lattice/gw.cpp | 65 +++++++--------- c++/triqs_tprf/lattice/gw.hpp | 54 ++++++++++++- python/triqs_tprf/lattice_desc.py | 7 ++ test/python/CMakeLists.txt | 1 + test/python/gw_dlr.py | 122 ++++++++++++++++++++++++++++++ 5 files changed, 209 insertions(+), 40 deletions(-) create mode 100644 test/python/gw_dlr.py diff --git a/c++/triqs_tprf/lattice/gw.cpp b/c++/triqs_tprf/lattice/gw.cpp index e8ccc6628..8c3272c4e 100644 --- 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 + 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 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 index 574c342db..d0d166065 100644 --- 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 index a466c4f07..794a149f5 100644 --- 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 index a0e31db3b..7838bd917 100644 --- 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 new file mode 100644 index 000000000..a70adfb34 --- /dev/null +++ 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() +