From 0254d634563d42c026442a9eba5e632b145fb733 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Mon, 7 May 2018 20:18:53 +0900 Subject: [PATCH 1/8] 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); From 429ce9c904214a7844e23e97c392aa11809aad42 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 00:40:27 +0900 Subject: [PATCH 2/8] Update python wrapper --- c++/g2_parameters.hpp | 8 +-- c++/pomerol_ed.hpp | 3 +- python/pomerol2triqs_converters.hxx | 108 ++++++++++++++++++++++++++++ python/pomerol2triqs_desc.py | 5 +- python/pomerol2triqs_parameters.rst | 26 ++++++- 5 files changed, 142 insertions(+), 8 deletions(-) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index 82216ea..d911df2 100644 --- a/c++/g2_parameters.hpp +++ b/c++/g2_parameters.hpp @@ -29,8 +29,8 @@ namespace pomerol2triqs { /// indices of operators in TRIQS convention: (block_name, inner_index) indices_t index1, index2, index3, index4; - 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) {} + // 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 { @@ -50,8 +50,8 @@ namespace pomerol2triqs { /// 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) {} + // g2_three_freqs_params_t() {} + // g2_three_freqs_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} }; /* diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index cd58066..4213915 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -119,10 +119,11 @@ 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 + /// Two-particle Green's function. Specify frequency cutoff, n_b and n_f. TRIQS_WRAP_ARG_AS_DICT g2_t G2_iw(g2_iw_inu_inup_params_t const &p); + /// Two-particle Green's function. Specify three frequencies (wb, wf1, wf2). TRIQS_WRAP_ARG_AS_DICT g2_three_freqs_t G2_iw_three_freqs(g2_three_freqs_params_t const &p); diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx index c8eae7d..b20165f 100644 --- a/python/pomerol2triqs_converters.hxx +++ b/python/pomerol2triqs_converters.hxx @@ -111,4 +111,112 @@ template <> struct py_converter { } }; +}} + + +// --- C++ Python converter for g2_three_freqs_params_t +#include +#include +#include + +namespace triqs { namespace py_tools { + +template <> struct py_converter { + static PyObject *c2py(g2_three_freqs_params_t const & x) { + PyObject * d = PyDict_New(); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "three_freqs", convert_to_python(x.three_freqs)); + PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); + PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); + PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); + PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); + return d; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r, U const &init_default) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = init_default; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = T{}; + } + + static g2_three_freqs_params_t py2c(PyObject *dic) { + g2_three_freqs_params_t res; + res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); + res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); + _get_optional(dic, "channel" , res.channel ,PH); + res.three_freqs = convert_from_python(PyDict_GetItemString(dic, "three_freqs")); + res.index1 = convert_from_python(PyDict_GetItemString(dic, "index1")); + res.index2 = convert_from_python(PyDict_GetItemString(dic, "index2")); + res.index3 = convert_from_python(PyDict_GetItemString(dic, "index3")); + res.index4 = convert_from_python(PyDict_GetItemString(dic, "index4")); + return res; + } + + template + static void _check(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!convertible_from_python(PyDict_GetItemString(dic, name), false)) + fs << "\n" << ++err << " The parameter " << name << " does not have the right type : expecting " << tname + << " in C++, but got '" << PyDict_GetItemString(dic, name)->ob_type->tp_name << "' in Python."; + } + + template + static void _check_mandatory(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!PyDict_Contains(dic, pyref::string(name))) + fs << "\n" << ++err << " Mandatory parameter " << name << " is missing."; + else _check(dic,fs,err,name,tname); + } + + template + static void _check_optional(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (PyDict_Contains(dic, pyref::string(name))) _check(dic, fs, err, name, tname); + } + + static bool is_convertible(PyObject *dic, bool raise_exception) { + if (dic == nullptr or !PyDict_Check(dic)) { + if (raise_exception) { PyErr_SetString(PyExc_TypeError, "The function must be called with named arguments");} + return false; + } + std::stringstream fs, fs2; int err=0; + +#ifndef TRIQS_ALLOW_UNUSED_PARAMETERS + std::vector ks, all_keys = {"gf_struct","beta","channel","three_freqs","index1","index2","index3","index4"}; + pyref keys = PyDict_Keys(dic); + if (!convertible_from_python>(keys, true)) { + fs << "\nThe dict keys are not strings"; + goto _error; + } + ks = convert_from_python>(keys); + for (auto & k : ks) + if (std::find(all_keys.begin(), all_keys.end(), k) == all_keys.end()) + fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; +#endif + + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "three_freqs", "three_freqs_t"); + _check_mandatory(dic, fs, err, "index1" , "indices_t"); + _check_mandatory(dic, fs, err, "index2" , "indices_t"); + _check_mandatory(dic, fs, err, "index3" , "indices_t"); + _check_mandatory(dic, fs, err, "index4" , "indices_t"); + if (err) goto _error; + return true; + + _error: + fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_three_freqs_params_t\n" <, 3> G2_iw (**pomerol2triqs::g2_iw_inu_inup_params_t)""", - doc = r"""Two-particle Green's function, Matsubara frequencies""") + doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") + +c.add_method("""std::vector > G2_iw_three_freqs (**pomerol2triqs::g2_three_freqs_params_t)""", + doc = r"""Two-particle Green's function. Specify three frequencies (wb, wf1, wf2).""") module.add_class(c) diff --git a/python/pomerol2triqs_parameters.rst b/python/pomerol2triqs_parameters.rst index be8cd2e..6dddaaf 100644 --- a/python/pomerol2triqs_parameters.rst +++ b/python/pomerol2triqs_parameters.rst @@ -7,9 +7,31 @@ +----------------+--------------------------+---------+----------------------------------------------------------------------+ | channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | +----------------+--------------------------+---------+----------------------------------------------------------------------+ -| n_b | int | -- | Number of bosonic Matsubara frequencies. | +| n_b | int | -- | Number of bosonic and fermionic Matsubara frequencies. | +----------------+--------------------------+---------+----------------------------------------------------------------------+ -| n_f | int | -- | Number of fermionic Matsubara frequencies. | +| n_f | int | -- | Number of bosonic and fermionic Matsubara frequencies. | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| index1 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| index2 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| index3 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| index4 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ + + + ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++================+==========================+=========+======================================================================+ +| gf_struct | gf_struct_t | -- | Block structure of GF | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | ++----------------+--------------------------+---------+----------------------------------------------------------------------+ +| three_freqs | three_freqs_t | -- | three frequencies (wb, wf1, wf2). | +----------------+--------------------------+---------+----------------------------------------------------------------------+ | index1 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | +----------------+--------------------------+---------+----------------------------------------------------------------------+ From 9fb8b76ac280ad3e4bc88915ff70ab41a4482175 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 02:33:42 +0900 Subject: [PATCH 3/8] Change the name convention for G2 --- c++/g2_parameters.hpp | 18 +++-- c++/pomerol_ed.cpp | 30 ++++++-- c++/pomerol_ed.hpp | 13 ++-- example/2band.atom.py | 2 +- example/slater.py | 2 +- python/pomerol2triqs_converters.hxx | 105 ++++++++++++--------------- python/pomerol2triqs_desc.py | 4 +- test/python/anderson_g2_matsubara.py | 12 +-- test/python/wick.py | 12 +-- 9 files changed, 105 insertions(+), 93 deletions(-) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index d911df2..a8a9984 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 four_indices_t = std::tuple; using three_freqs_t = std::vector< std::tuple >; namespace pomerol2triqs { @@ -12,7 +13,7 @@ namespace pomerol2triqs { // using g2_blocks_t = std::set>; - struct g2_iw_inu_inup_params_t { + struct g2_iw_freq_box_params_t { /// Block structure of GF gf_struct_t gf_struct; @@ -29,11 +30,11 @@ namespace pomerol2triqs { /// indices of operators in TRIQS convention: (block_name, inner_index) indices_t index1, index2, index3, index4; - // 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) {} + // g2_iw_freq_box_params_t() {} + // g2_iw_freq_box_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} }; - struct g2_three_freqs_params_t { + struct g2_iw_freq_vec_params_t { /// Block structure of GF gf_struct_t gf_struct; @@ -47,11 +48,12 @@ namespace pomerol2triqs { /// 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; + /// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 + // indices_t index1, index2, index3, index4; + std::vector vec_four_indices; - // g2_three_freqs_params_t() {} - // g2_three_freqs_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + // g2_iw_freq_vec_params_t() {} + // g2_iw_freq_vec_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} }; /* diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 48561c6..45fa862 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -492,7 +492,7 @@ namespace pomerol2triqs { pom_g2.compute(false, {}, comm); // compute g2 value using MPI - g2_three_freqs_t g2( three_freqs.size() ); + g2_iw_freq_vec_t g2( three_freqs.size() ); for(int i=comm.rank(); i(three_freqs[i]); int wf1 = std::get<1>(three_freqs[i]); @@ -510,8 +510,22 @@ namespace pomerol2triqs { return g2; } + std::vector< std::vector > > pomerol_ed::compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs){ + + std::vector vec_g2; + // TODO: MPI + for( auto four_indices : vec_four_indices ){ + indices_t index1 = std::get<0>(four_indices); + indices_t index2 = std::get<1>(four_indices); + indices_t index3 = std::get<2>(four_indices); + indices_t index4 = std::get<3>(four_indices); + vec_g2.push_back( compute_g2(gf_struct, beta, channel, index1, index2, index3, index4, three_freqs) ); + } + return vec_g2; + } - auto pomerol_ed::G2_iw(g2_iw_inu_inup_params_t const &p) -> g2_t { + + auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> g2_iw_freq_box_t { // create a list of three frequencies, (wb, wf1, wf2) three_freqs_t three_freqs; @@ -538,11 +552,10 @@ namespace pomerol2triqs { } // compute g2 values - // g2_three_freqs_t g2_three_freqs = G2_iw_three_freqs(p_core); - g2_three_freqs_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); + g2_iw_freq_vec_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); // reshape G2 - g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f); + g2_iw_freq_box_t g2(p.n_b, 2*p.n_f, 2*p.n_f); for(int i=0; i(three_indices[i]); int if1 = std::get<1>(three_indices[i]); @@ -554,8 +567,11 @@ namespace pomerol2triqs { } - 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_iw_freq_vec_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_freqs_vec(g2_iw_freq_vec_params_t const &p) -> std::vector { + return compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, p.three_freqs); } /* auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_three_freqs_t { diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index 4213915..c5946f1 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -82,12 +82,15 @@ 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 >; + using g2_iw_freq_box_t = triqs::arrays::array, 3>; + using g2_iw_freq_vec_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); + g2_iw_freq_vec_t 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); + + // std::vector compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs); + std::vector< g2_iw_freq_vec_t > compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs); double density_matrix_cutoff = 1e-15; @@ -121,11 +124,11 @@ namespace pomerol2triqs { /// Two-particle Green's function. Specify frequency cutoff, n_b and n_f. TRIQS_WRAP_ARG_AS_DICT - g2_t G2_iw(g2_iw_inu_inup_params_t const &p); + g2_iw_freq_box_t G2_iw_freq_box(g2_iw_freq_box_params_t const &p); /// Two-particle Green's function. Specify three frequencies (wb, wf1, wf2). TRIQS_WRAP_ARG_AS_DICT - g2_three_freqs_t G2_iw_three_freqs(g2_three_freqs_params_t const &p); + std::vector G2_iw_freqs_vec(g2_iw_freq_vec_params_t const &p); /// Two-particle Green's function, Matsubara frequencies // TRIQS_WRAP_ARG_AS_DICT diff --git a/example/2band.atom.py b/example/2band.atom.py index 415d90e..3b2b7c4 100644 --- a/example/2band.atom.py +++ b/example/2band.atom.py @@ -95,7 +95,7 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) print type(G2_iw) print G2_iw.shape diff --git a/example/slater.py b/example/slater.py index b11781e..1a34571 100644 --- a/example/slater.py +++ b/example/slater.py @@ -145,7 +145,7 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) if mpi.is_master_node(): print type(G2_iw) diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx index b20165f..c0e460e 100644 --- a/python/pomerol2triqs_converters.hxx +++ b/python/pomerol2triqs_converters.hxx @@ -3,25 +3,21 @@ // c++2py.py ../c++/pomerol_ed.hpp -I../../../local/pomerol/include -I/usr/include/eigen3 -I../c++ -p -mpytriqs.applications.impurity_solvers.pomerol2triqs -o pomerol2triqs --moduledoc "TRIQS wrapper around Pomerol ED library" -// --- C++ Python converter for g2_iw_inu_inup_params_t +// --- C++ Python converter for g2_iw_freq_vec_params_t #include #include #include namespace triqs { namespace py_tools { -template <> struct py_converter { - static PyObject *c2py(g2_iw_inu_inup_params_t const & x) { +template <> struct py_converter { + static PyObject *c2py(g2_iw_freq_vec_params_t const & x) { PyObject * d = PyDict_New(); - PyDict_SetItemString( d, "gf_struct", convert_to_python(x.gf_struct)); - PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); - PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); - PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); - PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); - PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); - PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); - PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); - PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "three_freqs" , convert_to_python(x.three_freqs)); + PyDict_SetItemString( d, "vec_four_indices", convert_to_python(x.vec_four_indices)); return d; } @@ -39,17 +35,13 @@ template <> struct py_converter { r = T{}; } - static g2_iw_inu_inup_params_t py2c(PyObject *dic) { - g2_iw_inu_inup_params_t res; + static g2_iw_freq_vec_params_t py2c(PyObject *dic) { + g2_iw_freq_vec_params_t res; res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); - _get_optional(dic, "channel" , res.channel ,PH); - res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); - res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); - res.index1 = convert_from_python(PyDict_GetItemString(dic, "index1")); - res.index2 = convert_from_python(PyDict_GetItemString(dic, "index2")); - res.index3 = convert_from_python(PyDict_GetItemString(dic, "index3")); - res.index4 = convert_from_python(PyDict_GetItemString(dic, "index4")); + _get_optional(dic, "channel" , res.channel ,PH); + res.three_freqs = convert_from_python(PyDict_GetItemString(dic, "three_freqs")); + res.vec_four_indices = convert_from_python>(PyDict_GetItemString(dic, "vec_four_indices")); return res; } @@ -80,7 +72,7 @@ template <> struct py_converter { std::stringstream fs, fs2; int err=0; #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS - std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"}; + std::vector ks, all_keys = {"gf_struct","beta","channel","three_freqs","vec_four_indices"}; pyref keys = PyDict_Keys(dic); if (!convertible_from_python>(keys, true)) { fs << "\nThe dict keys are not strings"; @@ -92,20 +84,16 @@ template <> struct py_converter { fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; #endif - _check_mandatory(dic, fs, err, "gf_struct", "gf_struct_t"); - _check_mandatory(dic, fs, err, "beta" , "double"); - _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); - _check_mandatory(dic, fs, err, "n_b" , "int"); - _check_mandatory(dic, fs, err, "n_f" , "int"); - _check_mandatory(dic, fs, err, "index1" , "indices_t"); - _check_mandatory(dic, fs, err, "index2" , "indices_t"); - _check_mandatory(dic, fs, err, "index3" , "indices_t"); - _check_mandatory(dic, fs, err, "index4" , "indices_t"); + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "three_freqs" , "three_freqs_t"); + _check_mandatory>(dic, fs, err, "vec_four_indices", "std::vector"); if (err) goto _error; return true; _error: - fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_inu_inup_params_t\n" < 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_vec_params_t\n" < struct py_converter { }} -// --- C++ Python converter for g2_three_freqs_params_t +// --- C++ Python converter for g2_iw_freq_box_params_t #include #include #include namespace triqs { namespace py_tools { -template <> struct py_converter { - static PyObject *c2py(g2_three_freqs_params_t const & x) { +template <> struct py_converter { + static PyObject *c2py(g2_iw_freq_box_params_t const & x) { PyObject * d = PyDict_New(); - PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); - PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); - PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); - PyDict_SetItemString( d, "three_freqs", convert_to_python(x.three_freqs)); - PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); - PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); - PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); - PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); + PyDict_SetItemString( d, "gf_struct", convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); + PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); + PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); + PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); + PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); + PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); return d; } @@ -149,12 +138,13 @@ template <> struct py_converter { r = T{}; } - static g2_three_freqs_params_t py2c(PyObject *dic) { - g2_three_freqs_params_t res; + static g2_iw_freq_box_params_t py2c(PyObject *dic) { + g2_iw_freq_box_params_t res; res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); - _get_optional(dic, "channel" , res.channel ,PH); - res.three_freqs = convert_from_python(PyDict_GetItemString(dic, "three_freqs")); + _get_optional(dic, "channel" , res.channel ,PH); + res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); + res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); res.index1 = convert_from_python(PyDict_GetItemString(dic, "index1")); res.index2 = convert_from_python(PyDict_GetItemString(dic, "index2")); res.index3 = convert_from_python(PyDict_GetItemString(dic, "index3")); @@ -189,7 +179,7 @@ template <> struct py_converter { std::stringstream fs, fs2; int err=0; #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS - std::vector ks, all_keys = {"gf_struct","beta","channel","three_freqs","index1","index2","index3","index4"}; + std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"}; pyref keys = PyDict_Keys(dic); if (!convertible_from_python>(keys, true)) { fs << "\nThe dict keys are not strings"; @@ -201,19 +191,20 @@ template <> struct py_converter { fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; #endif - _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); - _check_mandatory(dic, fs, err, "beta" , "double"); - _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); - _check_mandatory(dic, fs, err, "three_freqs", "three_freqs_t"); - _check_mandatory(dic, fs, err, "index1" , "indices_t"); - _check_mandatory(dic, fs, err, "index2" , "indices_t"); - _check_mandatory(dic, fs, err, "index3" , "indices_t"); - _check_mandatory(dic, fs, err, "index4" , "indices_t"); + _check_mandatory(dic, fs, err, "gf_struct", "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "n_b" , "int"); + _check_mandatory(dic, fs, err, "n_f" , "int"); + _check_mandatory(dic, fs, err, "index1" , "indices_t"); + _check_mandatory(dic, fs, err, "index2" , "indices_t"); + _check_mandatory(dic, fs, err, "index3" , "indices_t"); + _check_mandatory(dic, fs, err, "index4" , "indices_t"); if (err) goto _error; return true; _error: - fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_three_freqs_params_t\n" < 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_box_params_t\n" < G_w (gf_struct_t gf_struct, double beta, std::pair energy_window, int n_w, double im_shift = 0)""", doc = r"""Retarded Green's function on real energy axis""") -c.add_method("""triqs::arrays::array, 3> G2_iw (**pomerol2triqs::g2_iw_inu_inup_params_t)""", +c.add_method("""triqs::arrays::array, 3> G2_iw_freq_box (**pomerol2triqs::g2_iw_freq_box_params_t)""", doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") -c.add_method("""std::vector > G2_iw_three_freqs (**pomerol2triqs::g2_three_freqs_params_t)""", +c.add_method("""std::vector > > G2_iw_freqs_vec (**pomerol2triqs::g2_iw_freq_vec_params_t)""", doc = r"""Two-particle Green's function. Specify three frequencies (wb, wf1, wf2).""") module.add_class(c) diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py index 0c6baf4..b716700 100644 --- a/test/python/anderson_g2_matsubara.py +++ b/test/python/anderson_g2_matsubara.py @@ -73,14 +73,14 @@ # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order -G2_iw_ph_uuuu = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_dddd = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_uudd = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_dduu = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_uuuu = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_dddd = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_uudd = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_dduu = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order -G2_iw_ph_uddu = ed.G2_iw( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_duud = ed.G2_iw( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_uddu = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_duud = ed.G2_iw_freq_box( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) # Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order # G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", diff --git a/test/python/wick.py b/test/python/wick.py index 55b517c..3837ec2 100644 --- a/test/python/wick.py +++ b/test/python/wick.py @@ -81,12 +81,12 @@ 'n_f' : g2_n_wf, 'n_b' : g2_n_wb, } -G2_ph_uuuu = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_dddd = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_uudd = ed.G2_iw( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_dduu = ed.G2_iw( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_uddu = ed.G2_iw( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_ph_duud = ed.G2_iw( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) +G2_ph_uuuu = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_ph_dddd = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_ph_uudd = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_ph_dduu = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_ph_uddu = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +G2_ph_duud = ed.G2_iw_freq_box( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) G2_ph_uuuu_wick = G2_ph_uuuu.copy() G2_ph_dddd_wick = G2_ph_dddd.copy() From 92897ab1c90509aed5e2fab09a01eb5e69ac47dc Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 03:37:03 +0900 Subject: [PATCH 4/8] Change interface of G2 calc: Pass sets of four operator indices as a vector --- c++/g2_parameters.hpp | 27 +++-- c++/pomerol_ed.cpp | 41 ++++++-- c++/pomerol_ed.hpp | 6 +- example/2band.atom.py | 2 +- example/slater.py | 2 +- python/pomerol2triqs_converters.hxx | 150 ++++++++++++++++++++++----- python/pomerol2triqs_desc.py | 7 +- python/pomerol2triqs_parameters.rst | 50 +++++---- test/python/anderson_g2_matsubara.py | 12 +-- test/python/wick.py | 12 +-- 10 files changed, 232 insertions(+), 77 deletions(-) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index a8a9984..168d168 100644 --- a/c++/g2_parameters.hpp +++ b/c++/g2_parameters.hpp @@ -13,7 +13,7 @@ namespace pomerol2triqs { // using g2_blocks_t = std::set>; - struct g2_iw_freq_box_params_t { + struct g2_iw_legacy_params_t { /// Block structure of GF gf_struct_t gf_struct; @@ -30,8 +30,26 @@ namespace pomerol2triqs { /// indices of operators in TRIQS convention: (block_name, inner_index) indices_t index1, index2, index3, index4; - // g2_iw_freq_box_params_t() {} - // g2_iw_freq_box_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + // g2_iw_legacy_params_t() {} + // g2_iw_legacy_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} + }; + + struct g2_iw_freq_box_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; + + /// Number of bosonic and fermionic Matsubara frequencies. + int n_b, n_f; + + /// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 + std::vector vec_four_indices; }; struct g2_iw_freq_vec_params_t { @@ -51,9 +69,6 @@ namespace pomerol2triqs { /// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 // indices_t index1, index2, index3, index4; std::vector vec_four_indices; - - // g2_iw_freq_vec_params_t() {} - // g2_iw_freq_vec_params_t(gf_struct_t const &gf_struct, double beta) : gf_struct(gf_struct), beta(beta) {} }; /* diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 45fa862..8a56472 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -525,7 +525,21 @@ namespace pomerol2triqs { } - auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> g2_iw_freq_box_t { + auto pomerol_ed::G2_iw_legacy(g2_iw_legacy_params_t const &p) -> g2_iw_freq_box_t { + g2_iw_freq_box_params_t p2; + p2.gf_struct = p.gf_struct; + p2.beta = p.beta; + p2.channel = p.channel; + p2.n_b = p.n_b; + p2.n_f = p.n_f; + four_indices_t four_indices = std::make_tuple(p.index1, p.index2, p.index3, p.index4); + p2.vec_four_indices.push_back( four_indices ); + + std::vector vec_g2 = G2_iw_freq_box(p2); + return vec_g2[0]; + } + + auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> std::vector { // create a list of three frequencies, (wb, wf1, wf2) three_freqs_t three_freqs; @@ -552,18 +566,23 @@ namespace pomerol2triqs { } // compute g2 values - g2_iw_freq_vec_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); - - // reshape G2 - g2_iw_freq_box_t g2(p.n_b, 2*p.n_f, 2*p.n_f); - for(int i=0; i(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]; + // g2_iw_freq_vec_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); + std::vector vec_g2_freq_vec = compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, three_freqs); + + // reshape G2 (from freq_vec to freq_box) + std::vector vec_g2_freq_box; + for( auto g2_freq_vec : vec_g2_freq_vec ){ + g2_iw_freq_box_t g2(p.n_b, 2*p.n_f, 2*p.n_f); + for(int i=0; i(three_indices[i]); + int if1 = std::get<1>(three_indices[i]); + int if2 = std::get<2>(three_indices[i]); + g2(ib, if1, if2) = g2_freq_vec[i]; + } + vec_g2_freq_box.push_back(g2); } - return g2; + return vec_g2_freq_box; } diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index c5946f1..661e71f 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -124,7 +124,11 @@ namespace pomerol2triqs { /// Two-particle Green's function. Specify frequency cutoff, n_b and n_f. TRIQS_WRAP_ARG_AS_DICT - g2_iw_freq_box_t G2_iw_freq_box(g2_iw_freq_box_params_t const &p); + g2_iw_freq_box_t G2_iw_legacy(g2_iw_legacy_params_t const &p); + + /// Two-particle Green's function. Specify frequency cutoff, n_b and n_f. + TRIQS_WRAP_ARG_AS_DICT + std::vector G2_iw_freq_box(g2_iw_freq_box_params_t const &p); /// Two-particle Green's function. Specify three frequencies (wb, wf1, wf2). TRIQS_WRAP_ARG_AS_DICT diff --git a/example/2band.atom.py b/example/2band.atom.py index 3b2b7c4..68d7cce 100644 --- a/example/2band.atom.py +++ b/example/2band.atom.py @@ -95,7 +95,7 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) print type(G2_iw) print G2_iw.shape diff --git a/example/slater.py b/example/slater.py index 1a34571..cb975b1 100644 --- a/example/slater.py +++ b/example/slater.py @@ -145,7 +145,7 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) if mpi.is_master_node(): print type(G2_iw) diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx index c0e460e..abc54be 100644 --- a/python/pomerol2triqs_converters.hxx +++ b/python/pomerol2triqs_converters.hxx @@ -3,6 +3,117 @@ // c++2py.py ../c++/pomerol_ed.hpp -I../../../local/pomerol/include -I/usr/include/eigen3 -I../c++ -p -mpytriqs.applications.impurity_solvers.pomerol2triqs -o pomerol2triqs --moduledoc "TRIQS wrapper around Pomerol ED library" +// --- C++ Python converter for g2_iw_legacy_params_t +#include +#include +#include + +namespace triqs { namespace py_tools { + +template <> struct py_converter { + static PyObject *c2py(g2_iw_legacy_params_t const & x) { + PyObject * d = PyDict_New(); + PyDict_SetItemString( d, "gf_struct", convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); + PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); + PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); + PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); + PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); + PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); + return d; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r, U const &init_default) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = init_default; + } + + template static void _get_optional(PyObject *dic, const char *name, T &r) { + if (PyDict_Contains(dic, pyref::string(name))) + r = convert_from_python(PyDict_GetItemString(dic, name)); + else + r = T{}; + } + + static g2_iw_legacy_params_t py2c(PyObject *dic) { + g2_iw_legacy_params_t res; + res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); + res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); + _get_optional(dic, "channel" , res.channel ,PH); + res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); + res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); + res.index1 = convert_from_python(PyDict_GetItemString(dic, "index1")); + res.index2 = convert_from_python(PyDict_GetItemString(dic, "index2")); + res.index3 = convert_from_python(PyDict_GetItemString(dic, "index3")); + res.index4 = convert_from_python(PyDict_GetItemString(dic, "index4")); + return res; + } + + template + static void _check(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!convertible_from_python(PyDict_GetItemString(dic, name), false)) + fs << "\n" << ++err << " The parameter " << name << " does not have the right type : expecting " << tname + << " in C++, but got '" << PyDict_GetItemString(dic, name)->ob_type->tp_name << "' in Python."; + } + + template + static void _check_mandatory(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (!PyDict_Contains(dic, pyref::string(name))) + fs << "\n" << ++err << " Mandatory parameter " << name << " is missing."; + else _check(dic,fs,err,name,tname); + } + + template + static void _check_optional(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) { + if (PyDict_Contains(dic, pyref::string(name))) _check(dic, fs, err, name, tname); + } + + static bool is_convertible(PyObject *dic, bool raise_exception) { + if (dic == nullptr or !PyDict_Check(dic)) { + if (raise_exception) { PyErr_SetString(PyExc_TypeError, "The function must be called with named arguments");} + return false; + } + std::stringstream fs, fs2; int err=0; + +#ifndef TRIQS_ALLOW_UNUSED_PARAMETERS + std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"}; + pyref keys = PyDict_Keys(dic); + if (!convertible_from_python>(keys, true)) { + fs << "\nThe dict keys are not strings"; + goto _error; + } + ks = convert_from_python>(keys); + for (auto & k : ks) + if (std::find(all_keys.begin(), all_keys.end(), k) == all_keys.end()) + fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; +#endif + + _check_mandatory(dic, fs, err, "gf_struct", "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "n_b" , "int"); + _check_mandatory(dic, fs, err, "n_f" , "int"); + _check_mandatory(dic, fs, err, "index1" , "indices_t"); + _check_mandatory(dic, fs, err, "index2" , "indices_t"); + _check_mandatory(dic, fs, err, "index3" , "indices_t"); + _check_mandatory(dic, fs, err, "index4" , "indices_t"); + if (err) goto _error; + return true; + + _error: + fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_legacy_params_t\n" < #include @@ -112,15 +223,12 @@ namespace triqs { namespace py_tools { template <> struct py_converter { static PyObject *c2py(g2_iw_freq_box_params_t const & x) { PyObject * d = PyDict_New(); - PyDict_SetItemString( d, "gf_struct", convert_to_python(x.gf_struct)); - PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); - PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); - PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); - PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); - PyDict_SetItemString( d, "index1" , convert_to_python(x.index1)); - PyDict_SetItemString( d, "index2" , convert_to_python(x.index2)); - PyDict_SetItemString( d, "index3" , convert_to_python(x.index3)); - PyDict_SetItemString( d, "index4" , convert_to_python(x.index4)); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); + PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); + PyDict_SetItemString( d, "vec_four_indices", convert_to_python(x.vec_four_indices)); return d; } @@ -142,13 +250,10 @@ template <> struct py_converter { g2_iw_freq_box_params_t res; res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); - _get_optional(dic, "channel" , res.channel ,PH); + _get_optional(dic, "channel" , res.channel ,PH); res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); - res.index1 = convert_from_python(PyDict_GetItemString(dic, "index1")); - res.index2 = convert_from_python(PyDict_GetItemString(dic, "index2")); - res.index3 = convert_from_python(PyDict_GetItemString(dic, "index3")); - res.index4 = convert_from_python(PyDict_GetItemString(dic, "index4")); + res.vec_four_indices = convert_from_python>(PyDict_GetItemString(dic, "vec_four_indices")); return res; } @@ -179,7 +284,7 @@ template <> struct py_converter { std::stringstream fs, fs2; int err=0; #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS - std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"}; + std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","vec_four_indices"}; pyref keys = PyDict_Keys(dic); if (!convertible_from_python>(keys, true)) { fs << "\nThe dict keys are not strings"; @@ -191,15 +296,12 @@ template <> struct py_converter { fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; #endif - _check_mandatory(dic, fs, err, "gf_struct", "gf_struct_t"); - _check_mandatory(dic, fs, err, "beta" , "double"); - _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); - _check_mandatory(dic, fs, err, "n_b" , "int"); - _check_mandatory(dic, fs, err, "n_f" , "int"); - _check_mandatory(dic, fs, err, "index1" , "indices_t"); - _check_mandatory(dic, fs, err, "index2" , "indices_t"); - _check_mandatory(dic, fs, err, "index3" , "indices_t"); - _check_mandatory(dic, fs, err, "index4" , "indices_t"); + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "n_b" , "int"); + _check_mandatory(dic, fs, err, "n_f" , "int"); + _check_mandatory>(dic, fs, err, "vec_four_indices", "std::vector"); if (err) goto _error; return true; diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py index 8eb4ae6..0a5b27d 100644 --- a/python/pomerol2triqs_desc.py +++ b/python/pomerol2triqs_desc.py @@ -67,10 +67,13 @@ 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 = r"""Retarded Green's function on real energy axis""") -c.add_method("""triqs::arrays::array, 3> G2_iw_freq_box (**pomerol2triqs::g2_iw_freq_box_params_t)""", +c.add_method("""triqs::arrays::array, 3> G2_iw_legacy (**pomerol2triqs::g2_iw_legacy_params_t)""", doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") -c.add_method("""std::vector > > G2_iw_freqs_vec (**pomerol2triqs::g2_iw_freq_vec_params_t)""", +c.add_method("""std::vector< triqs::arrays::array, 3> > G2_iw_freq_box (**pomerol2triqs::g2_iw_freq_box_params_t)""", + doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") + +c.add_method("""std::vector< std::vector > > G2_iw_freqs_vec (**pomerol2triqs::g2_iw_freq_vec_params_t)""", doc = r"""Two-particle Green's function. Specify three frequencies (wb, wf1, wf2).""") module.add_class(c) diff --git a/python/pomerol2triqs_parameters.rst b/python/pomerol2triqs_parameters.rst index 6dddaaf..2e303b1 100644 --- a/python/pomerol2triqs_parameters.rst +++ b/python/pomerol2triqs_parameters.rst @@ -22,22 +22,34 @@ -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| Parameter Name | Type | Default | Documentation | -+================+==========================+=========+======================================================================+ -| gf_struct | gf_struct_t | -- | Block structure of GF | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| beta | double | -- | Inverse temperature | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| three_freqs | three_freqs_t | -- | three frequencies (wb, wf1, wf2). | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| index1 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| index2 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| index3 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ -| index4 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | -+----------------+--------------------------+---------+----------------------------------------------------------------------+ \ No newline at end of file ++------------------+-----------------------------+---------+------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++==================+=============================+=========+==================================================================+ +| gf_struct | gf_struct_t | -- | Block structure of GF | ++------------------+-----------------------------+---------+------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++------------------+-----------------------------+---------+------------------------------------------------------------------+ +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | ++------------------+-----------------------------+---------+------------------------------------------------------------------+ +| three_freqs | three_freqs_t | -- | three frequencies (wb, wf1, wf2). | ++------------------+-----------------------------+---------+------------------------------------------------------------------+ +| vec_four_indices | std::vector | -- | | ++------------------+-----------------------------+---------+------------------------------------------------------------------+ + + + ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++==================+=============================+=========+====================================================================================+ +| gf_struct | gf_struct_t | -- | Block structure of GF | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| n_b | int | -- | Number of bosonic and fermionic Matsubara frequencies. | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| n_f | int | -- | Number of bosonic and fermionic Matsubara frequencies. | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ +| vec_four_indices | std::vector | -- | set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 | ++------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ \ No newline at end of file diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py index b716700..a8ef2e5 100644 --- a/test/python/anderson_g2_matsubara.py +++ b/test/python/anderson_g2_matsubara.py @@ -73,14 +73,14 @@ # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order -G2_iw_ph_uuuu = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_dddd = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_uudd = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_dduu = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order -G2_iw_ph_uddu = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_duud = ed.G2_iw_freq_box( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) +G2_iw_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +G2_iw_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) # Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order # G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", diff --git a/test/python/wick.py b/test/python/wick.py index 3837ec2..858ff72 100644 --- a/test/python/wick.py +++ b/test/python/wick.py @@ -81,12 +81,12 @@ 'n_f' : g2_n_wf, 'n_b' : g2_n_wb, } -G2_ph_uuuu = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_dddd = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_uudd = ed.G2_iw_freq_box( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_dduu = ed.G2_iw_freq_box( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_uddu = ed.G2_iw_freq_box( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_ph_duud = ed.G2_iw_freq_box( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) +G2_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +G2_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +G2_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +G2_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) G2_ph_uuuu_wick = G2_ph_uuuu.copy() G2_ph_dddd_wick = G2_ph_dddd.copy() From 19eab2ab28fc617eb2f62e8d9d55ed4d2f463985 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 04:08:05 +0900 Subject: [PATCH 5/8] Use the new interface in test and examples --- example/2band.atom.py | 3 +- example/slater.py | 3 +- test/python/anderson_g2_matsubara.py | 36 +++++++++++++++----- test/python/wick.py | 50 ++++++++++++++++++++-------- 4 files changed, 68 insertions(+), 24 deletions(-) diff --git a/example/2band.atom.py b/example/2band.atom.py index 68d7cce..2f1628f 100644 --- a/example/2band.atom.py +++ b/example/2band.atom.py @@ -95,7 +95,8 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0] print type(G2_iw) print G2_iw.shape diff --git a/example/slater.py b/example/slater.py index cb975b1..ff227cb 100644 --- a/example/slater.py +++ b/example/slater.py @@ -145,7 +145,8 @@ # G^{(2)}(i\omega;i\nu,i\nu') # ############################### -G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) +G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0] if mpi.is_master_node(): print type(G2_iw) diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py index a8ef2e5..26c9785 100644 --- a/test/python/anderson_g2_matsubara.py +++ b/test/python/anderson_g2_matsubara.py @@ -72,15 +72,33 @@ 'n_b' : g2_n_wb, } -# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order -G2_iw_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_iw_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) - -# Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order -G2_iw_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_iw_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) +# # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order +# G2_iw_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +# G2_iw_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +# G2_iw_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +# G2_iw_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +# +# # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), ABBA block order +# G2_iw_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +# G2_iw_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) + + +four_indices = [] +four_indices.append( (('up',0), ('up',0), ('up',0), ('up',0)) ) # uuuu +four_indices.append( (('dn',0), ('dn',0), ('dn',0), ('dn',0)) ) # dddd +four_indices.append( (('up',0), ('up',0), ('dn',0), ('dn',0)) ) # uudd +four_indices.append( (('dn',0), ('dn',0), ('up',0), ('up',0)) ) # dduu +four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) ) # uddu +four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) ) # duud + +G2_iw_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params ) +G2_iw_ph_uuuu = G2_iw_ph[0] +G2_iw_ph_dddd = G2_iw_ph[1] +G2_iw_ph_uudd = G2_iw_ph[2] +G2_iw_ph_dduu = G2_iw_ph[3] +G2_iw_ph_uddu = G2_iw_ph[4] +G2_iw_ph_duud = G2_iw_ph[5] + # Compute G^{(2),pp}(i\omega;i\nu,i\nu'), AABB block order # G2_iw_inu_inup_pp_AABB = ed.G2_iw_inu_inup(channel = "PP", diff --git a/test/python/wick.py b/test/python/wick.py index 858ff72..d2a3c6d 100644 --- a/test/python/wick.py +++ b/test/python/wick.py @@ -81,19 +81,43 @@ 'n_f' : g2_n_wf, 'n_b' : g2_n_wb, } -G2_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) -G2_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) -G2_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) -G2_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) - -G2_ph_uuuu_wick = G2_ph_uuuu.copy() -G2_ph_dddd_wick = G2_ph_dddd.copy() -G2_ph_uudd_wick = G2_ph_uudd.copy() -G2_ph_dduu_wick = G2_ph_dduu.copy() -G2_ph_uddu_wick = G2_ph_uddu.copy() -G2_ph_duud_wick = G2_ph_duud.copy() +# G2_ph_uuuu = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('up',0), index4=('up',0), **common_g2_params ) +# G2_ph_dddd = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +# G2_ph_uudd = ed.G2_iw_legacy( index1=('up',0), index2=('up',0), index3=('dn',0), index4=('dn',0), **common_g2_params ) +# G2_ph_dduu = ed.G2_iw_legacy( index1=('dn',0), index2=('dn',0), index3=('up',0), index4=('up',0), **common_g2_params ) +# G2_ph_uddu = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',0), index4=('up',0), **common_g2_params ) +# G2_ph_duud = ed.G2_iw_legacy( index1=('dn',0), index2=('up',0), index3=('up',0), index4=('dn',0), **common_g2_params ) + +four_indices = [] +four_indices.append( (('up',0), ('up',0), ('up',0), ('up',0)) ) # uuuu +four_indices.append( (('dn',0), ('dn',0), ('dn',0), ('dn',0)) ) # dddd +four_indices.append( (('up',0), ('up',0), ('dn',0), ('dn',0)) ) # uudd +four_indices.append( (('dn',0), ('dn',0), ('up',0), ('up',0)) ) # dduu +four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) ) # uddu +four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) ) # duud + +G2_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params ) +G2_ph_uuuu = G2_ph[0] +G2_ph_dddd = G2_ph[1] +G2_ph_uudd = G2_ph[2] +G2_ph_dduu = G2_ph[3] +G2_ph_uddu = G2_ph[4] +G2_ph_duud = G2_ph[5] + +# G2_ph_uuuu_wick = G2_ph_uuuu.copy() +# G2_ph_dddd_wick = G2_ph_dddd.copy() +# G2_ph_uudd_wick = G2_ph_uudd.copy() +# G2_ph_dduu_wick = G2_ph_dduu.copy() +# G2_ph_uddu_wick = G2_ph_uddu.copy() +# G2_ph_duud_wick = G2_ph_duud.copy() + +g2_shape = (g2_n_wb, 2*g2_n_wf, 2*g2_n_wf) +G2_ph_uuuu_wick = np.zeros(g2_shape, dtype=np.complex) +G2_ph_dddd_wick = np.zeros(g2_shape, dtype=np.complex) +G2_ph_dduu_wick = np.zeros(g2_shape, dtype=np.complex) +G2_ph_uudd_wick = np.zeros(g2_shape, dtype=np.complex) +G2_ph_uddu_wick = np.zeros(g2_shape, dtype=np.complex) +G2_ph_duud_wick = np.zeros(g2_shape, dtype=np.complex) G = lambda s, i: G_iw[s].data[i + n_iw, 0, 0] From 70721e45de2c8f53906bb7fc8907d7dd5de9a9cd Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 09:53:05 +0900 Subject: [PATCH 6/8] Change some names and types --- c++/g2_parameters.hpp | 33 ++++++----- c++/pomerol_ed.cpp | 44 +++++++------- c++/pomerol_ed.hpp | 16 +++--- example/2band.atom.py | 10 ++-- example/slater.py | 31 +++++++--- python/pomerol2triqs_converters.hxx | 86 ++++++++++++++-------------- python/pomerol2triqs_desc.py | 8 +-- python/pomerol2triqs_parameters.rst | 62 ++++++++++---------- test/python/anderson_g2_matsubara.py | 2 +- test/python/wick.py | 2 +- 10 files changed, 157 insertions(+), 137 deletions(-) diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index 168d168..584dc32 100644 --- a/c++/g2_parameters.hpp +++ b/c++/g2_parameters.hpp @@ -4,7 +4,7 @@ using triqs::hilbert_space::gf_struct_t; using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t; using four_indices_t = std::tuple; -using three_freqs_t = std::vector< std::tuple >; +using three_freqs_t = std::tuple; namespace pomerol2triqs { @@ -21,11 +21,13 @@ namespace pomerol2triqs { /// Inverse temperature double beta; - /// Channel in which Matsubara frequency representation is defined. + /// Channel in which Matsubara frequency representation is defined channel_t channel = PH; - /// Number of bosonic and fermionic Matsubara frequencies. - int n_b, n_f; + /// Number of non-negative bosonic Matsubara frequencies + int n_b; + /// Number of positive fermionic Matsubara frequencies + int n_f; /// indices of operators in TRIQS convention: (block_name, inner_index) indices_t index1, index2, index3, index4; @@ -42,17 +44,19 @@ namespace pomerol2triqs { /// Inverse temperature double beta; - /// Channel in which Matsubara frequency representation is defined. + /// Channel in which Matsubara frequency representation is defined channel_t channel = PH; - /// Number of bosonic and fermionic Matsubara frequencies. - int n_b, n_f; + /// Number of non-negative bosonic Matsubara frequencies + int n_b; + /// Number of positive fermionic Matsubara frequencies + int n_f; - /// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 - std::vector vec_four_indices; + /// four operator indices in TRIQS convention: (block_name, inner_index)*4 + std::vector four_indices; }; - struct g2_iw_freq_vec_params_t { + struct g2_iw_freq_fix_params_t { /// Block structure of GF gf_struct_t gf_struct; @@ -60,15 +64,14 @@ namespace pomerol2triqs { /// Inverse temperature double beta; - /// Channel in which Matsubara frequency representation is defined. + /// Channel in which Matsubara frequency representation is defined channel_t channel = PH; /// three frequencies (wb, wf1, wf2). - three_freqs_t three_freqs; + std::vector three_freqs; - /// set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 - // indices_t index1, index2, index3, index4; - std::vector vec_four_indices; + /// four operator indices in TRIQS convention: (block_name, inner_index)*4 + std::vector four_indices; }; /* diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 8a56472..5f138d4 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -469,7 +469,7 @@ namespace pomerol2triqs { */ // 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){ + triqs::arrays::array, 1> pomerol_ed::compute_g2_core(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, std::vector const &three_freqs){ if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized"; compute_rho(beta); @@ -492,34 +492,35 @@ namespace pomerol2triqs { pom_g2.compute(false, {}, comm); // compute g2 value using MPI - g2_iw_freq_vec_t g2( three_freqs.size() ); + g2_iw_freq_fix_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); + g2(i) = -pom_g2(wb+wf1, wf2, wf1); } // broadcast results for(int i=0; i > > pomerol_ed::compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs){ - std::vector vec_g2; + std::vector< triqs::arrays::array, 1> > pomerol_ed::compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &four_indices, std::vector const &three_freqs){ + + std::vector vec_g2; // TODO: MPI - for( auto four_indices : vec_four_indices ){ - indices_t index1 = std::get<0>(four_indices); - indices_t index2 = std::get<1>(four_indices); - indices_t index3 = std::get<2>(four_indices); - indices_t index4 = std::get<3>(four_indices); - vec_g2.push_back( compute_g2(gf_struct, beta, channel, index1, index2, index3, index4, three_freqs) ); + for( auto ind4 : four_indices ){ + indices_t index1 = std::get<0>(ind4); + indices_t index2 = std::get<1>(ind4); + indices_t index3 = std::get<2>(ind4); + indices_t index4 = std::get<3>(ind4); + vec_g2.push_back( compute_g2_core(gf_struct, beta, channel, index1, index2, index3, index4, three_freqs) ); } return vec_g2; } @@ -532,8 +533,8 @@ namespace pomerol2triqs { p2.channel = p.channel; p2.n_b = p.n_b; p2.n_f = p.n_f; - four_indices_t four_indices = std::make_tuple(p.index1, p.index2, p.index3, p.index4); - p2.vec_four_indices.push_back( four_indices ); + four_indices_t ind4 = std::make_tuple(p.index1, p.index2, p.index3, p.index4); + p2.four_indices.push_back( ind4 ); std::vector vec_g2 = G2_iw_freq_box(p2); return vec_g2[0]; @@ -542,7 +543,7 @@ namespace pomerol2triqs { auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> std::vector { // create a list of three frequencies, (wb, wf1, wf2) - three_freqs_t three_freqs; + std::vector three_freqs; std::vector< std::tuple > three_indices; // (ib, if1, if2) { // indices of fermionic Matsubara frequencies [-n_f:n_f) @@ -566,8 +567,8 @@ namespace pomerol2triqs { } // compute g2 values - // g2_iw_freq_vec_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); - std::vector vec_g2_freq_vec = compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, three_freqs); + // g2_iw_freq_fix_t g2_three_freqs = compute_g2(p.gf_struct, p.beta, p.channel, p.index1, p.index2, p.index3, p.index4, three_freqs); + std::vector vec_g2_freq_vec = compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, three_freqs); // reshape G2 (from freq_vec to freq_box) std::vector vec_g2_freq_box; @@ -577,7 +578,7 @@ namespace pomerol2triqs { int ib = std::get<0>(three_indices[i]); int if1 = std::get<1>(three_indices[i]); int if2 = std::get<2>(three_indices[i]); - g2(ib, if1, if2) = g2_freq_vec[i]; + g2(ib, if1, if2) = g2_freq_vec(i); } vec_g2_freq_box.push_back(g2); } @@ -586,11 +587,8 @@ namespace pomerol2triqs { } - // auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_iw_freq_vec_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_freqs_vec(g2_iw_freq_vec_params_t const &p) -> std::vector { - return compute_g2_indices_loop(p.gf_struct, p.beta, p.channel, p.vec_four_indices, p.three_freqs); + auto pomerol_ed::G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p) -> std::vector { + return compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, p.three_freqs); } /* auto pomerol_ed::G2_iw_three_freqs(g2_three_freqs_params_t const &p) -> g2_three_freqs_t { diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index 661e71f..5b4f4e5 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -83,14 +83,14 @@ namespace pomerol2triqs { // using w_nu_nup_t = cartesian_product; // using w_l_lp_t = cartesian_product; using g2_iw_freq_box_t = triqs::arrays::array, 3>; - using g2_iw_freq_vec_t = std::vector >; + using g2_iw_freq_fix_t = triqs::arrays::array, 1>; // 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; - g2_iw_freq_vec_t 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); - // std::vector compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs); - std::vector< g2_iw_freq_vec_t > compute_g2_indices_loop(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &vec_four_indices, three_freqs_t const &three_freqs); + g2_iw_freq_fix_t compute_g2_core(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, std::vector const &three_freqs); + + std::vector compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector const &four_indices, std::vector const &three_freqs); double density_matrix_cutoff = 1e-15; @@ -122,17 +122,17 @@ 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. Specify frequency cutoff, n_b and n_f. + /// Two-particle Green's function, Legacy support TRIQS_WRAP_ARG_AS_DICT g2_iw_freq_box_t G2_iw_legacy(g2_iw_legacy_params_t const &p); - /// Two-particle Green's function. Specify frequency cutoff, n_b and n_f. + /// Two-particle Green's function, in a low-frequency box with cutoff n_b and n_f TRIQS_WRAP_ARG_AS_DICT std::vector G2_iw_freq_box(g2_iw_freq_box_params_t const &p); - /// Two-particle Green's function. Specify three frequencies (wb, wf1, wf2). + /// Two-particle Green's function, for fixed frequencies (wb, wf1, wf2) TRIQS_WRAP_ARG_AS_DICT - std::vector G2_iw_freqs_vec(g2_iw_freq_vec_params_t const &p); + std::vector G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p); /// Two-particle Green's function, Matsubara frequencies // TRIQS_WRAP_ARG_AS_DICT diff --git a/example/2band.atom.py b/example/2band.atom.py index 2f1628f..24a5841 100644 --- a/example/2band.atom.py +++ b/example/2band.atom.py @@ -96,10 +96,12 @@ ############################### # G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) -G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0] -print type(G2_iw) -print G2_iw.shape +four_indices = [(('up',0), ('dn',0), ('dn',1), ('up',1))] +G2_iw = ed.G2_iw_freq_box( four_indices=four_indices, **common_g2_params ) + +print "G2_iw :", type(G2_iw), "of size", len(G2_iw) +print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape # # Compute G^{(2),ph}(i\omega;i\nu,i\nu'), AABB block order # G2_iw_inu_inup_ph_AABB = ed.G2_iw_inu_inup(channel = "PH", @@ -126,4 +128,4 @@ ar['G_iw'] = G_iw ar['G_tau'] = G_tau ar['G_w'] = G_w - ar['G2_ph'] = G2_iw + ar['G2_ph'] = G2_iw[0] diff --git a/example/slater.py b/example/slater.py index ff227cb..4a41cbf 100644 --- a/example/slater.py +++ b/example/slater.py @@ -137,19 +137,36 @@ common_g2_params = {'channel' : "PH", 'gf_struct' : gf_struct, - 'beta' : beta, - 'n_f' : g2_n_wf, - 'n_b' : g2_n_wb, } + 'beta' : beta,} ############################### # G^{(2)}(i\omega;i\nu,i\nu') # ############################### # G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params ) -G2_iw = ed.G2_iw_freq_box( vec_four_indices=[(('up',0), ('dn',0), ('dn',1), ('up',1)),], **common_g2_params )[0] + +# four-operators indices +four_indices = [] +four_indices.append( (('up',0), ('dn',0), ('dn',1), ('up',1)) ) +four_indices.append( (('up',0), ('up',0), ('dn',1), ('dn',1)) ) + +# compute G2 in a low-frequency box +G2_iw = ed.G2_iw_freq_box( four_indices=four_indices, n_f=g2_n_wf, n_b=g2_n_wb, **common_g2_params ) + +if mpi.is_master_node(): + print "G2_iw :", type(G2_iw), "of size", len(G2_iw) + print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape + with HDFArchive('slater_gf.out.h5', 'a') as ar: + ar['G2_iw_freq_box'] = G2_iw + +# compute G2 for given freqs, (wb, wf1, wf2) +three_freqs = [] +three_freqs.append( (0, 1, 2) ) +three_freqs.append( (0, 1, -2) ) +G2_iw = ed.G2_iw_freq_fix( four_indices=four_indices, three_freqs=three_freqs, **common_g2_params ) if mpi.is_master_node(): - print type(G2_iw) - print G2_iw.shape + print "G2_iw :", type(G2_iw), "of size", len(G2_iw) + print "G2_iw[0] :", type(G2_iw[0]), "of size", G2_iw[0].shape with HDFArchive('slater_gf.out.h5', 'a') as ar: - ar['G2_iw'] = G2_iw + ar['G2_iw_freq_fix'] = G2_iw diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx index abc54be..aa543f1 100644 --- a/python/pomerol2triqs_converters.hxx +++ b/python/pomerol2triqs_converters.hxx @@ -114,21 +114,22 @@ template <> struct py_converter { }} -// --- C++ Python converter for g2_iw_freq_vec_params_t +// --- C++ Python converter for g2_iw_freq_box_params_t #include #include #include namespace triqs { namespace py_tools { -template <> struct py_converter { - static PyObject *c2py(g2_iw_freq_vec_params_t const & x) { +template <> struct py_converter { + static PyObject *c2py(g2_iw_freq_box_params_t const & x) { PyObject * d = PyDict_New(); - PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); - PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); - PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); - PyDict_SetItemString( d, "three_freqs" , convert_to_python(x.three_freqs)); - PyDict_SetItemString( d, "vec_four_indices", convert_to_python(x.vec_four_indices)); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); + PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); + PyDict_SetItemString( d, "four_indices", convert_to_python(x.four_indices)); return d; } @@ -146,13 +147,14 @@ template <> struct py_converter { r = T{}; } - static g2_iw_freq_vec_params_t py2c(PyObject *dic) { - g2_iw_freq_vec_params_t res; + static g2_iw_freq_box_params_t py2c(PyObject *dic) { + g2_iw_freq_box_params_t res; res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); - _get_optional(dic, "channel" , res.channel ,PH); - res.three_freqs = convert_from_python(PyDict_GetItemString(dic, "three_freqs")); - res.vec_four_indices = convert_from_python>(PyDict_GetItemString(dic, "vec_four_indices")); + _get_optional(dic, "channel" , res.channel ,PH); + res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); + res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); + res.four_indices = convert_from_python>(PyDict_GetItemString(dic, "four_indices")); return res; } @@ -183,7 +185,7 @@ template <> struct py_converter { std::stringstream fs, fs2; int err=0; #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS - std::vector ks, all_keys = {"gf_struct","beta","channel","three_freqs","vec_four_indices"}; + std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","four_indices"}; pyref keys = PyDict_Keys(dic); if (!convertible_from_python>(keys, true)) { fs << "\nThe dict keys are not strings"; @@ -195,16 +197,17 @@ template <> struct py_converter { fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; #endif - _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); - _check_mandatory(dic, fs, err, "beta" , "double"); - _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); - _check_mandatory(dic, fs, err, "three_freqs" , "three_freqs_t"); - _check_mandatory>(dic, fs, err, "vec_four_indices", "std::vector"); + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory(dic, fs, err, "n_b" , "int"); + _check_mandatory(dic, fs, err, "n_f" , "int"); + _check_mandatory>(dic, fs, err, "four_indices", "std::vector"); if (err) goto _error; return true; _error: - fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_vec_params_t\n" < 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_box_params_t\n" < struct py_converter { }} -// --- C++ Python converter for g2_iw_freq_box_params_t +// --- C++ Python converter for g2_iw_freq_fix_params_t #include #include #include namespace triqs { namespace py_tools { -template <> struct py_converter { - static PyObject *c2py(g2_iw_freq_box_params_t const & x) { +template <> struct py_converter { + static PyObject *c2py(g2_iw_freq_fix_params_t const & x) { PyObject * d = PyDict_New(); - PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); - PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); - PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); - PyDict_SetItemString( d, "n_b" , convert_to_python(x.n_b)); - PyDict_SetItemString( d, "n_f" , convert_to_python(x.n_f)); - PyDict_SetItemString( d, "vec_four_indices", convert_to_python(x.vec_four_indices)); + PyDict_SetItemString( d, "gf_struct" , convert_to_python(x.gf_struct)); + PyDict_SetItemString( d, "beta" , convert_to_python(x.beta)); + PyDict_SetItemString( d, "channel" , convert_to_python(x.channel)); + PyDict_SetItemString( d, "three_freqs" , convert_to_python(x.three_freqs)); + PyDict_SetItemString( d, "four_indices", convert_to_python(x.four_indices)); return d; } @@ -246,14 +248,13 @@ template <> struct py_converter { r = T{}; } - static g2_iw_freq_box_params_t py2c(PyObject *dic) { - g2_iw_freq_box_params_t res; + static g2_iw_freq_fix_params_t py2c(PyObject *dic) { + g2_iw_freq_fix_params_t res; res.gf_struct = convert_from_python(PyDict_GetItemString(dic, "gf_struct")); res.beta = convert_from_python(PyDict_GetItemString(dic, "beta")); - _get_optional(dic, "channel" , res.channel ,PH); - res.n_b = convert_from_python(PyDict_GetItemString(dic, "n_b")); - res.n_f = convert_from_python(PyDict_GetItemString(dic, "n_f")); - res.vec_four_indices = convert_from_python>(PyDict_GetItemString(dic, "vec_four_indices")); + _get_optional(dic, "channel" , res.channel ,PH); + res.three_freqs = convert_from_python>(PyDict_GetItemString(dic, "three_freqs")); + res.four_indices = convert_from_python>(PyDict_GetItemString(dic, "four_indices")); return res; } @@ -284,7 +285,7 @@ template <> struct py_converter { std::stringstream fs, fs2; int err=0; #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS - std::vector ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","vec_four_indices"}; + std::vector ks, all_keys = {"gf_struct","beta","channel","three_freqs","four_indices"}; pyref keys = PyDict_Keys(dic); if (!convertible_from_python>(keys, true)) { fs << "\nThe dict keys are not strings"; @@ -296,17 +297,16 @@ template <> struct py_converter { fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized."; #endif - _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); - _check_mandatory(dic, fs, err, "beta" , "double"); - _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); - _check_mandatory(dic, fs, err, "n_b" , "int"); - _check_mandatory(dic, fs, err, "n_f" , "int"); - _check_mandatory>(dic, fs, err, "vec_four_indices", "std::vector"); + _check_mandatory(dic, fs, err, "gf_struct" , "gf_struct_t"); + _check_mandatory(dic, fs, err, "beta" , "double"); + _check_optional (dic, fs, err, "channel" , "pomerol2triqs::channel_t"); + _check_mandatory >(dic, fs, err, "three_freqs" , "std::vector"); + _check_mandatory>(dic, fs, err, "four_indices", "std::vector"); if (err) goto _error; return true; _error: - fs2 << "\n---- There " << (err > 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_box_params_t\n" < 1 ? "are " : "is ") << err<< " error"<<(err >1 ?"s" : "")<< " in Python -> C++ transcription for the class g2_iw_freq_fix_params_t\n" <, 3> G2_iw_legacy (**pomerol2triqs::g2_iw_legacy_params_t)""", - doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") + doc = r"""Two-particle Green's function, Legacy support""") c.add_method("""std::vector< triqs::arrays::array, 3> > G2_iw_freq_box (**pomerol2triqs::g2_iw_freq_box_params_t)""", - doc = r"""Two-particle Green's function. Specify frequency cutoff, n_b and n_f.""") + doc = r"""Two-particle Green's function, in a low-frequency box with cutoff n_b and n_f""") -c.add_method("""std::vector< std::vector > > G2_iw_freqs_vec (**pomerol2triqs::g2_iw_freq_vec_params_t)""", - doc = r"""Two-particle Green's function. Specify three frequencies (wb, wf1, wf2).""") +c.add_method("""std::vector< triqs::arrays::array, 1> > G2_iw_freq_fix (**pomerol2triqs::g2_iw_freq_fix_params_t)""", + doc = r"""Two-particle Green's function, for fixed frequencies (wb, wf1, wf2)""") module.add_class(c) diff --git a/python/pomerol2triqs_parameters.rst b/python/pomerol2triqs_parameters.rst index 2e303b1..4030bee 100644 --- a/python/pomerol2triqs_parameters.rst +++ b/python/pomerol2triqs_parameters.rst @@ -5,11 +5,11 @@ +----------------+--------------------------+---------+----------------------------------------------------------------------+ | beta | double | -- | Inverse temperature | +----------------+--------------------------+---------+----------------------------------------------------------------------+ -| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined | +----------------+--------------------------+---------+----------------------------------------------------------------------+ -| n_b | int | -- | Number of bosonic and fermionic Matsubara frequencies. | +| n_b | int | -- | Number of non-negative bosonic Matsubara frequencies | +----------------+--------------------------+---------+----------------------------------------------------------------------+ -| n_f | int | -- | Number of bosonic and fermionic Matsubara frequencies. | +| n_f | int | -- | Number of positive fermionic Matsubara frequencies | +----------------+--------------------------+---------+----------------------------------------------------------------------+ | index1 | indices_t | -- | indices of operators in TRIQS convention: (block_name, inner_index) | +----------------+--------------------------+---------+----------------------------------------------------------------------+ @@ -22,34 +22,34 @@ -+------------------+-----------------------------+---------+------------------------------------------------------------------+ -| Parameter Name | Type | Default | Documentation | -+==================+=============================+=========+==================================================================+ -| gf_struct | gf_struct_t | -- | Block structure of GF | -+------------------+-----------------------------+---------+------------------------------------------------------------------+ -| beta | double | -- | Inverse temperature | -+------------------+-----------------------------+---------+------------------------------------------------------------------+ -| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | -+------------------+-----------------------------+---------+------------------------------------------------------------------+ -| three_freqs | three_freqs_t | -- | three frequencies (wb, wf1, wf2). | -+------------------+-----------------------------+---------+------------------------------------------------------------------+ -| vec_four_indices | std::vector | -- | | -+------------------+-----------------------------+---------+------------------------------------------------------------------+ ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++================+=============================+=========+=========================================================================+ +| gf_struct | gf_struct_t | -- | Block structure of GF | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| n_b | int | -- | Number of non-negative bosonic Matsubara frequencies | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| n_f | int | -- | Number of positive fermionic Matsubara frequencies | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| four_indices | std::vector | -- | four operator indices in TRIQS convention: (block_name, inner_index)*4 | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| Parameter Name | Type | Default | Documentation | -+==================+=============================+=========+====================================================================================+ -| gf_struct | gf_struct_t | -- | Block structure of GF | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| beta | double | -- | Inverse temperature | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined. | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| n_b | int | -- | Number of bosonic and fermionic Matsubara frequencies. | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| n_f | int | -- | Number of bosonic and fermionic Matsubara frequencies. | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ -| vec_four_indices | std::vector | -- | set of indices of four operators in TRIQS convention: (block_name, inner_index)*4 | -+------------------+-----------------------------+---------+------------------------------------------------------------------------------------+ \ No newline at end of file ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| Parameter Name | Type | Default | Documentation | ++================+=============================+=========+=========================================================================+ +| gf_struct | gf_struct_t | -- | Block structure of GF | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| beta | double | -- | Inverse temperature | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| channel | pomerol2triqs::channel_t | PH | Channel in which Matsubara frequency representation is defined | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| three_freqs | std::vector | -- | three frequencies (wb, wf1, wf2). | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ +| four_indices | std::vector | -- | four operator indices in TRIQS convention: (block_name, inner_index)*4 | ++----------------+-----------------------------+---------+-------------------------------------------------------------------------+ \ No newline at end of file diff --git a/test/python/anderson_g2_matsubara.py b/test/python/anderson_g2_matsubara.py index 26c9785..e61f8fc 100644 --- a/test/python/anderson_g2_matsubara.py +++ b/test/python/anderson_g2_matsubara.py @@ -91,7 +91,7 @@ four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) ) # uddu four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) ) # duud -G2_iw_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params ) +G2_iw_ph = ed.G2_iw_freq_box( four_indices=four_indices, **common_g2_params ) G2_iw_ph_uuuu = G2_iw_ph[0] G2_iw_ph_dddd = G2_iw_ph[1] G2_iw_ph_uudd = G2_iw_ph[2] diff --git a/test/python/wick.py b/test/python/wick.py index d2a3c6d..a5d621f 100644 --- a/test/python/wick.py +++ b/test/python/wick.py @@ -96,7 +96,7 @@ four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) ) # uddu four_indices.append( (('dn',0), ('up',0), ('up',0), ('dn',0)) ) # duud -G2_ph = ed.G2_iw_freq_box( vec_four_indices=four_indices, **common_g2_params ) +G2_ph = ed.G2_iw_freq_box( four_indices=four_indices, **common_g2_params ) G2_ph_uuuu = G2_ph[0] G2_ph_dddd = G2_ph[1] G2_ph_uudd = G2_ph[2] From dee49fa88fc44db3bfcbe282b722790e8a9fceff Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 11:08:38 +0900 Subject: [PATCH 7/8] Clean up --- c++/pomerol_ed.cpp | 118 +-------------------------------------------- 1 file changed, 1 insertion(+), 117 deletions(-) diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index 5f138d4..8fb9049 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -397,76 +397,6 @@ 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); - 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); - - // 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; - - // std::cout << "Start freq loop: rank" << comm.rank() << std::endl; - // g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f); - // // std::cout << typeid(g2).name() << std::endl; - // for(int ib=0; ib > three_freqs; - for(int ib=0; ib(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]); - } - - // broadcast results - for(int i=0; i(three_freqs[i]); - int if1 = std::get<1>(three_freqs[i]); - int if2 = std::get<2>(three_freqs[i]); - int sender = i % comm.size(); - boost::mpi::broadcast(comm, g2(ib, if1, if2), sender); - } - - // std::cout << "End freq loop: rank" << comm.rank() << std::endl; - return g2; - } -*/ // core function for computing G2 triqs::arrays::array, 1> pomerol_ed::compute_g2_core(gf_struct_t const &gf_struct, double beta, channel_t channel, indices_t index1, indices_t index2, indices_t index3, indices_t index4, std::vector const &three_freqs){ @@ -540,6 +470,7 @@ namespace pomerol2triqs { return vec_g2[0]; } + auto pomerol_ed::G2_iw_freq_box(g2_iw_freq_box_params_t const &p) -> std::vector { // create a list of three frequencies, (wb, wf1, wf2) @@ -559,7 +490,6 @@ namespace pomerol2triqs { for(int ib=0; ib vec_g2_freq_vec = compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, three_freqs); // reshape G2 (from freq_vec to freq_box) @@ -590,51 +519,6 @@ namespace pomerol2triqs { auto pomerol_ed::G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p) -> std::vector { return compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, 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 From 1e591121d5aa47eb78e45ac5e4b786ce7c728be7 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Tue, 8 May 2018 11:38:48 +0900 Subject: [PATCH 8/8] Add a test of checking consistency bw G2_iw_freq_box() and G2_iw_freq_fix() --- test/python/CMakeLists.txt | 1 + test/python/g2_freq_box_and_freq_fix.py | 87 +++++++++++++++++++++++++ 2 files changed, 88 insertions(+) create mode 100644 test/python/g2_freq_box_and_freq_fix.py diff --git a/test/python/CMakeLists.txt b/test/python/CMakeLists.txt index de8460f..1c9e136 100644 --- a/test/python/CMakeLists.txt +++ b/test/python/CMakeLists.txt @@ -7,6 +7,7 @@ file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURREN triqs_add_python_test(anderson_gf) triqs_add_python_test(slater_gf) triqs_add_python_test(wick) +triqs_add_python_test(g2_freq_box_and_freq_fix) foreach(TEST_MPI_NUMPROC 1 2 4) triqs_add_python_test(anderson_g2_matsubara) diff --git a/test/python/g2_freq_box_and_freq_fix.py b/test/python/g2_freq_box_and_freq_fix.py new file mode 100644 index 0000000..883b6c5 --- /dev/null +++ b/test/python/g2_freq_box_and_freq_fix.py @@ -0,0 +1,87 @@ +from pytriqs.archive import HDFArchive +from pytriqs.gf import * +from pytriqs.operators import Operator, c, c_dag, n +from pytriqs.utility import mpi +from pytriqs.applications.impurity_solvers.pomerol2triqs import PomerolED +from pytriqs.utility.comparison_tests import * +import numpy as np +from itertools import product + +# Single orbital Anderson model + +#################### +# Input parameters # +#################### + +beta = 10.0 # Inverse temperature +U = 4.3 # Coulomb repulsion +mu = 1.1 # Chemical potential + +# Levels of the bath sites +epsilon = [-1.5, 1.3] +# Hopping amplitudes +V = [0.52, 0.55] + +spin_names = ("up", "dn") + +# Number of bosonic Matsubara frequencies for G^2 calculations +g2_n_wb = 3 +# Number of fermionic Matsubara frequencies for G^2 calculations +g2_n_wf = 3 + +# GF structure +gf_struct = {'up' : [0], 'dn' : [0]} + +# Conversion from TRIQS to Pomerol notation for operator indices +index_converter = {} +index_converter.update({(sn, 0) : ("loc", 0, "down" if sn == "dn" else "up") for sn in spin_names}) +index_converter.update({("B%i_%s" % (k, sn), 0) : ("bath" + str(k), 0, "down" if sn == "dn" else "up") + for k, sn in product(range(len(epsilon)), spin_names)}) + +# Make PomerolED solver object +ed = PomerolED(index_converter, verbose = True) + +# Number of particles on the impurity +H_loc = -mu*(n('up', 0) + n('dn', 0)) + U * n('up', 0) * n('dn', 0) + +# Bath Hamiltonian +H_bath = sum(eps*n("B%i_%s" % (k, sn), 0) + for sn, (k, eps) in product(spin_names, enumerate(epsilon))) + +# Hybridization Hamiltonian +H_hyb = Operator() +for k, v in enumerate(V): + H_hyb += sum( v * c_dag("B%i_%s" % (k, sn), 0) * c(sn, 0) + + np.conj(v) * c_dag(sn, 0) * c("B%i_%s" % (k, sn), 0) + for sn in spin_names) + +# Complete Hamiltonian +H = H_loc + H_hyb + H_bath + +# Diagonalize H +ed.diagonalize(H) + +########### +# G^{(2)} # +########### + +common_g2_params = {'channel' : "PH", + 'gf_struct' : gf_struct, + 'beta' : beta,} + +four_indices = [] +four_indices.append( (('up',0), ('dn',0), ('dn',0), ('up',0)) ) # uddu + +# compute in a low-freq box +G2_iw_freq_box = ed.G2_iw_freq_box( four_indices=four_indices, n_f=g2_n_wf, n_b=g2_n_wb, **common_g2_params )[0] +print G2_iw_freq_box.shape + +# compute for fixed freqs +three_freqs = [ (iwb, iwf1-g2_n_wf, iwf2-g2_n_wf) for iwb, iwf1, iwf2 in product( range(g2_n_wb), range(2*g2_n_wf), range(2*g2_n_wf) ) ] +G2_iw_freq_fix = ed.G2_iw_freq_fix( four_indices=four_indices, three_freqs=three_freqs, **common_g2_params )[0] +print G2_iw_freq_fix.shape +# reshape from 1D array to 3D array +G2_iw_freq_fix = np.reshape( G2_iw_freq_fix, (g2_n_wb, 2*g2_n_wf, 2*g2_n_wf) ) +print G2_iw_freq_fix.shape + +assert( np.allclose(G2_iw_freq_box, G2_iw_freq_fix) )