From 0254d634563d42c026442a9eba5e632b145fb733 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Mon, 7 May 2018 20:18:53 +0900 Subject: [PATCH] Enable arbitrary combinations of three freqs in G2 calc --- c++/g2_parameters.hpp | 30 +++++++-- c++/pomerol_ed.cpp | 138 +++++++++++++++++++++++++++++++++++++++++- c++/pomerol_ed.hpp | 5 ++ 3 files changed, 167 insertions(+), 6 deletions(-) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index 1127d40..82216ea 100644 --- a/c++/g2_parameters.hpp +++ b/c++/g2_parameters.hpp @@ -3,6 +3,7 @@ #include using triqs::hilbert_space::gf_struct_t; using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t; +using three_freqs_t = std::vector< std::tuple >; namespace pomerol2triqs { @@ -22,11 +23,8 @@ namespace pomerol2triqs { /// Channel in which Matsubara frequency representation is defined. channel_t channel = PH; - /// Number of bosonic Matsubara frequencies. - int n_b; - - /// Number of fermionic Matsubara frequencies. - int n_f; + /// Number of bosonic and fermionic Matsubara frequencies. + int n_b, n_f; /// indices of operators in TRIQS convention: (block_name, inner_index) indices_t index1, index2, index3, index4; @@ -34,6 +32,28 @@ namespace pomerol2triqs { g2_iw_inu_inup_params_t() {} g2_iw_inu_inup_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} }; + + struct g2_three_freqs_params_t { + + /// Block structure of GF + gf_struct_t gf_struct; + + /// Inverse temperature + double beta; + + /// Channel in which Matsubara frequency representation is defined. + channel_t channel = PH; + + /// three frequencies (wb, wf1, wf2). + three_freqs_t three_freqs; + + /// indices of operators in TRIQS convention: (block_name, inner_index) + indices_t index1, index2, index3, index4; + + g2_three_freqs_params_t() {} + g2_three_freqs_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + }; + /* struct g2_iw_l_lp_params_t { diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 06e94fe..48561c6 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -397,7 +397,7 @@ namespace pomerol2triqs { /////////////////////////////////// // Two-particle Green's function // /////////////////////////////////// - +/* auto pomerol_ed::G2_iw(g2_iw_inu_inup_params_t const &p) -> g2_t { if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized"; compute_rho(p.beta); @@ -466,6 +466,142 @@ namespace pomerol2triqs { // std::cout << "End freq loop: rank" << comm.rank() << std::endl; return g2; } +*/ + + // core function for computing G2 + std::vector > pomerol_ed::compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, three_freqs_t const &three_freqs){ + + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized"; + compute_rho(beta); + compute_field_operators(gf_struct); + + if (verbose && !comm.rank()) { std::cout << "\nPomerol: computing TwoParticleGF" << std::endl; } + + Pomerol::ParticleIndex pom_i1 = lookup_pomerol_index(index1); + Pomerol::ParticleIndex pom_i2 = lookup_pomerol_index(index2); + Pomerol::ParticleIndex pom_i3 = lookup_pomerol_index(index3); + Pomerol::ParticleIndex pom_i4 = lookup_pomerol_index(index4); + + Pomerol::TwoParticleGF pom_g2(*states_class, *matrix_h, + ops_container->getAnnihilationOperator(pom_i2), + ops_container->getAnnihilationOperator(pom_i4), + ops_container->getCreationOperator(pom_i1), + ops_container->getCreationOperator(pom_i3), + *rho); + pom_g2.prepare(); + pom_g2.compute(false, {}, comm); + + // compute g2 value using MPI + g2_three_freqs_t g2( three_freqs.size() ); + for(int i=comm.rank(); i(three_freqs[i]); + int wf1 = std::get<1>(three_freqs[i]); + int wf2 = std::get<2>(three_freqs[i]); + g2[i] = -pom_g2(wb+wf1, wf2, wf1); + } + + // broadcast results + for(int i=0; i g2_t { + + // create a list of three frequencies, (wb, wf1, wf2) + three_freqs_t three_freqs; + std::vector< std::tuple > three_indices; // (ib, if1, if2) + { + // indices of fermionic Matsubara frequencies [-n_f:n_f) + std::vector index_wf(2*p.n_f); + std::iota(index_wf.begin(), index_wf.end(), -p.n_f); + // for( auto i : iw_f ) std::cout << i << std::endl; + + // indices of bosonic Matsubara frequencies [0:n_b) + std::vector index_wb(p.n_b); + std::iota(index_wb.begin(), index_wb.end(), 0); + // for( auto i : iw_b ) std::cout << i << std::endl; + + for(int ib=0; ib(three_indices[i]); + int if1 = std::get<1>(three_indices[i]); + int if2 = std::get<2>(three_indices[i]); + g2(ib, if1, if2) = g2_three_freqs[i]; + } + + return g2; + } + + + auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_three_freqs_t { + return compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, p.three_freqs); + } +/* + auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_three_freqs_t { + if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized"; + compute_rho(p.beta); + compute_field_operators(p.gf_struct); + + if (verbose && !comm.rank()) { std::cout << "\nPomerol: computing TwoParticleGF" << std::endl; } + + Pomerol::ParticleIndex pom_i1 = lookup_pomerol_index(p.index1); + Pomerol::ParticleIndex pom_i2 = lookup_pomerol_index(p.index2); + Pomerol::ParticleIndex pom_i3 = lookup_pomerol_index(p.index3); + Pomerol::ParticleIndex pom_i4 = lookup_pomerol_index(p.index4); + + Pomerol::TwoParticleGF pom_g2(*states_class, *matrix_h, + ops_container->getAnnihilationOperator(pom_i2), + ops_container->getAnnihilationOperator(pom_i4), + ops_container->getCreationOperator(pom_i1), + ops_container->getCreationOperator(pom_i3), + *rho); + pom_g2.prepare(); + pom_g2.compute(false, {}, comm); + + // compute g2 value using MPI + g2_three_freqs_t g2( p.three_freqs.size() ); + for(int i=comm.rank(); i(three_freqs[i]); + // int if1 = std::get<1>(three_freqs[i]); + // int if2 = std::get<2>(three_freqs[i]); + // g2(ib, if1, if2) = -pom_g2(index_wb[ib]+index_wf[if1], index_wf[if2], index_wf[if1]); + int wb = std::get<0>(p.three_freqs[i]); + int wf1 = std::get<1>(p.three_freqs[i]); + int wf2 = std::get<2>(p.three_freqs[i]); + g2[i] = -pom_g2(wb+wf1, wf2, wf1); + } + + // broadcast results + for(int i=0; i diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index 6dd546c..cd58066 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -83,9 +83,11 @@ namespace pomerol2triqs { // using w_nu_nup_t = cartesian_product; // using w_l_lp_t = cartesian_product; using g2_t = triqs::arrays::array, 3>; + using g2_three_freqs_t = std::vector >; // 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; + std::vector > compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, three_freqs_t const &three_freqs); double density_matrix_cutoff = 1e-15; @@ -121,6 +123,9 @@ namespace pomerol2triqs { TRIQS_WRAP_ARG_AS_DICT g2_t G2_iw(g2_iw_inu_inup_params_t const &p); + TRIQS_WRAP_ARG_AS_DICT + g2_three_freqs_t G2_iw_three_freqs(g2_three_freqs_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);