Skip to content

Commit

Permalink
Add a method for G2 (not implemented)
Browse files Browse the repository at this point in the history
  • Loading branch information
@j-otsuki
j-otsuki committed Feb 27, 2018
1 parent a713bdb commit b2b9c42
Show file tree
Hide file tree
Showing 3 changed files with 70 additions and 56 deletions.
105 changes: 54 additions & 51 deletions c++/pomerol_ed.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/constants/constants.hpp>
#include <triqs/arrays.hpp>
// #include <triqs/arrays.hpp>
#include <fstream>
#include <algorithm>

Expand Down Expand Up @@ -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<w_nu_nup_t, tensor_valued<4>> {
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<w_nu_nup_t, tensor_valued<4>> {
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<w_nu_nup_t, scalar_valued> 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<imfreq> mesh_b{p.beta, Boson, p.n_iw};
gf_mesh<imfreq> mesh_f{p.beta, Fermion, p.n_inu};
gf_mesh<w_nu_nup_t> mesh_bff{mesh_b, mesh_f, mesh_f};
gf_mesh<w_nu_nup_t> mesh_fff{mesh_f, mesh_f, mesh_f};
block2_gf<w_nu_nup_t, tensor_valued<4>> g2;
if (p.channel == AllFermionic)
g2 = compute_g2<w_nu_nup_t>(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler);
else
g2 = compute_g2<w_nu_nup_t>(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<w_nu_nup_t, scalar_valued> 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<imfreq> mesh_b{p.beta, Boson, p.n_iw};
// gf_mesh<imfreq> mesh_f{p.beta, Fermion, p.n_inu};
//
// gf_mesh<w_nu_nup_t> mesh_bff{mesh_b, mesh_f, mesh_f};
// gf_mesh<w_nu_nup_t> mesh_fff{mesh_f, mesh_f, mesh_f};
//
// block2_gf<w_nu_nup_t, tensor_valued<4>> g2;
//
// if (p.channel == AllFermionic)
// g2 = compute_g2<w_nu_nup_t>(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler);
// else
// g2 = compute_g2<w_nu_nup_t>(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<w_l_lp_t, tensor_valued<4>> {
if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: no Hamiltonian has been diagonalized";
compute_rho(p.beta);
Expand Down
12 changes: 9 additions & 3 deletions c++/pomerol_ed.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,9 @@
#include <triqs/hilbert_space/fundamental_operator_set.hpp>
// #include <triqs/utility/optional_compat.hpp>
#include <triqs/mpi/boost.hpp>
#include <triqs/arrays.hpp>

// #include "g2_parameters.hpp"
#include "g2_parameters.hpp"

namespace pomerol2triqs {

Expand Down Expand Up @@ -60,8 +61,9 @@ namespace pomerol2triqs {
void compute_field_operators(gf_struct_t const &gf_struct);
template <typename Mesh, typename Filler> block_gf<Mesh> compute_gf(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, Filler filler) const;

using w_nu_nup_t = cartesian_product<imfreq, imfreq, imfreq>;
using w_l_lp_t = cartesian_product<imfreq, legendre, legendre>;
// using w_nu_nup_t = cartesian_product<imfreq, imfreq, imfreq>;
// using w_l_lp_t = cartesian_product<imfreq, legendre, legendre>;
using g2_t = array<std::complex<double>, 3>;
// template <typename Mesh, typename Filler>
// block2_gf<Mesh, tensor_valued<4>> compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order,
// g2_blocks_t const &g2_blocks, Filler filler) const;
Expand Down Expand Up @@ -96,6 +98,10 @@ namespace pomerol2triqs {
/// Retarded Green's function on real energy axis
block_gf<refreq> G_w(gf_struct_t const &gf_struct, double beta, std::pair<double, double> 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<w_nu_nup_t, tensor_valued<4>> G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p);
Expand Down
9 changes: 7 additions & 2 deletions python/pomerol2triqs_desc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

# Add here anything to add in the C++ code at the start, e.g. namespace using
module.add_preamble("""
#include <triqs/python_tools/converters/arrays.hpp>
#include <triqs/python_tools/converters/pair.hpp>
#include <triqs/python_tools/converters/map.hpp>
#include <triqs/python_tools/converters/set.hpp>
Expand All @@ -23,13 +24,14 @@
#include <triqs/python_tools/converters/tuple.hpp>
#include <triqs/python_tools/converters/operators_real_complex.hpp>
#include <triqs/python_tools/converters/gf.hpp>
#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_(
Expand Down Expand Up @@ -65,6 +67,9 @@
c.add_method("""block_gf<refreq> G_w (gf_struct_t gf_struct, double beta, std::pair<double,double> energy_window, int n_w, double im_shift = 0)""",
doc = """Retarded Green\'s function on real energy axis """)

c.add_method("""array<std::complex<double>, 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()

0 comments on commit b2b9c42

Please sign in to comment.