From b2b9c422cf64647854ae372581562b126986fbae Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 27 Feb 2018 16:50:43 +0900 Subject: [PATCH] Add a method for G2 (not implemented) --- c++/pomerol_ed.cpp | 105 ++++++++++++++++++----------------- c++/pomerol_ed.hpp | 12 +++- python/pomerol2triqs_desc.py | 9 ++- 3 files changed, 70 insertions(+), 56 deletions(-) diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 67464ca..b9fdc74 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -2,7 +2,7 @@ #include #include -#include +// #include #include #include @@ -453,65 +453,68 @@ namespace pomerol2triqs { return make_block2_gf(block_names, block_names, std::move(gf_vecvec)); } +*/ - auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf> { + auto pomerol_ed::G2_iw(g2_iw_inu_inup_params_t const &p) -> g2_t { + // auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf> { if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_inu_inup: no Hamiltonian has been diagonalized"; compute_rho(p.beta); compute_field_operators(p.gf_struct); if (verbose && !comm.rank()) std::cout << "G2_iw_inu_inup: filling output container" << std::endl; - auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { - long mesh_index = 0; - for (auto w_nu_nup : g2_el.mesh()) { - if ((mesh_index++) % comm.size() != comm.rank()) continue; - - if (p.channel == AllFermionic) { - - int n1 = std::get<0>(w_nu_nup).index(); - int n2 = std::get<1>(w_nu_nup).index(); - int n3 = std::get<2>(w_nu_nup).index(); - - if (p.block_order == AABB) - g2_el[w_nu_nup] = -pom_g2(n2, n1 + n3 - n2, n1); - else - g2_el[w_nu_nup] = +pom_g2(n1 + n3 - n2, n2, n1); - - } else { // p.channel == PH or PP - - int w_n = std::get<0>(w_nu_nup).index(); - int nu_n = std::get<1>(w_nu_nup).index(); - int nup_n = std::get<2>(w_nu_nup).index(); - - int W_n = p.channel == PH ? w_n + nu_n : w_n - nup_n - 1; - - if (p.block_order == AABB) { - g2_el[w_nu_nup] = -pom_g2(W_n, nup_n, nu_n); - } else { - g2_el[w_nu_nup] = +pom_g2(nup_n, W_n, nu_n); - } - } - } - }; - - gf_mesh mesh_b{p.beta, Boson, p.n_iw}; - gf_mesh mesh_f{p.beta, Fermion, p.n_inu}; - - gf_mesh mesh_bff{mesh_b, mesh_f, mesh_f}; - gf_mesh mesh_fff{mesh_f, mesh_f, mesh_f}; - - block2_gf> g2; - - if (p.channel == AllFermionic) - g2 = compute_g2(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler); - else - g2 = compute_g2(p.gf_struct, mesh_bff, p.block_order, p.blocks, filler); - - g2() = mpi_all_reduce(g2(), comm); - - return g2; + // auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { + // long mesh_index = 0; + // for (auto w_nu_nup : g2_el.mesh()) { + // if ((mesh_index++) % comm.size() != comm.rank()) continue; + // + // if (p.channel == AllFermionic) { + // + // int n1 = std::get<0>(w_nu_nup).index(); + // int n2 = std::get<1>(w_nu_nup).index(); + // int n3 = std::get<2>(w_nu_nup).index(); + // + // if (p.block_order == AABB) + // g2_el[w_nu_nup] = -pom_g2(n2, n1 + n3 - n2, n1); + // else + // g2_el[w_nu_nup] = +pom_g2(n1 + n3 - n2, n2, n1); + // + // } else { // p.channel == PH or PP + // + // int w_n = std::get<0>(w_nu_nup).index(); + // int nu_n = std::get<1>(w_nu_nup).index(); + // int nup_n = std::get<2>(w_nu_nup).index(); + // + // int W_n = p.channel == PH ? w_n + nu_n : w_n - nup_n - 1; + // + // if (p.block_order == AABB) { + // g2_el[w_nu_nup] = -pom_g2(W_n, nup_n, nu_n); + // } else { + // g2_el[w_nu_nup] = +pom_g2(nup_n, W_n, nu_n); + // } + // } + // } + // }; + + // gf_mesh mesh_b{p.beta, Boson, p.n_iw}; + // gf_mesh mesh_f{p.beta, Fermion, p.n_inu}; + // + // gf_mesh mesh_bff{mesh_b, mesh_f, mesh_f}; + // gf_mesh mesh_fff{mesh_f, mesh_f, mesh_f}; + // + // block2_gf> g2; + // + // if (p.channel == AllFermionic) + // g2 = compute_g2(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler); + // else + // g2 = compute_g2(p.gf_struct, mesh_bff, p.block_order, p.blocks, filler); + // + // g2() = mpi_all_reduce(g2(), comm); + // + // return g2; } +/* auto pomerol_ed::G2_iw_l_lp(g2_iw_l_lp_params_t const &p) -> block2_gf> { if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: no Hamiltonian has been diagonalized"; compute_rho(p.beta); diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index 60e37a3..a0b8542 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -12,8 +12,9 @@ #include // #include #include +#include -// #include "g2_parameters.hpp" +#include "g2_parameters.hpp" namespace pomerol2triqs { @@ -60,8 +61,9 @@ namespace pomerol2triqs { void compute_field_operators(gf_struct_t const &gf_struct); template block_gf compute_gf(gf_struct_t const &gf_struct, gf_mesh const &mesh, Filler filler) const; - using w_nu_nup_t = cartesian_product; - using w_l_lp_t = cartesian_product; + // using w_nu_nup_t = cartesian_product; + // using w_l_lp_t = cartesian_product; + using g2_t = array, 3>; // template // block2_gf> compute_g2(gf_struct_t const &gf_struct, gf_mesh const &mesh, block_order_t block_order, // g2_blocks_t const &g2_blocks, Filler filler) const; @@ -96,6 +98,10 @@ namespace pomerol2triqs { /// Retarded Green's function on real energy axis block_gf G_w(gf_struct_t const &gf_struct, double beta, std::pair const &energy_window, int n_w, double im_shift = 0); + /// Two-particle Green's function, Matsubara frequencies + TRIQS_WRAP_ARG_AS_DICT + g2_t G2_iw(g2_iw_inu_inup_params_t const &p); + /// Two-particle Green's function, Matsubara frequencies // TRIQS_WRAP_ARG_AS_DICT // block2_gf> G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p); diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py index f71b47d..43c9b0a 100644 --- a/python/pomerol2triqs_desc.py +++ b/python/pomerol2triqs_desc.py @@ -14,6 +14,7 @@ # Add here anything to add in the C++ code at the start, e.g. namespace using module.add_preamble(""" +#include #include #include #include @@ -23,13 +24,14 @@ #include #include #include +#include "./pomerol2triqs_converters.hxx" """) module.add_using("namespace pomerol2triqs") module.add_using("namespace Pomerol") module.add_enum("spin", ["down", "up"], "Pomerol", "Spin projection") -# module.add_enum("block_order_t", ["AABB", "ABBA"], "pomerol2triqs", "G^{(2)} block order") -# module.add_enum("channel_t", ["PP", "PH", "AllFermionic"], "pomerol2triqs", "G^{(2)} channel") +module.add_enum("block_order_t", ["AABB", "ABBA"], "pomerol2triqs", "G^{(2)} block order") +module.add_enum("channel_t", ["PP", "PH", "AllFermionic"], "pomerol2triqs", "G^{(2)} channel") # The class pomerol_ed c = class_( @@ -65,6 +67,9 @@ c.add_method("""block_gf G_w (gf_struct_t gf_struct, double beta, std::pair energy_window, int n_w, double im_shift = 0)""", doc = """Retarded Green\'s function on real energy axis """) +c.add_method("""array, 3> G2_iw (**pomerol2triqs::g2_iw_inu_inup_params_t)""", + doc = r"""Two-particle Green's function, Matsubara frequencies""") + module.add_class(c) module.generate_code()