diff --git a/c++/g2_parameters.hpp b/c++/g2_parameters.hpp index c828371..fd75245 100644 --- a/c++/g2_parameters.hpp +++ b/c++/g2_parameters.hpp @@ -2,13 +2,14 @@ #include using triqs::hilbert_space::gf_struct_t; +using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t; namespace pomerol2triqs { - enum block_order_t { AABB, ABBA }; + // enum block_order_t { AABB, ABBA }; enum channel_t { PP, PH, AllFermionic }; - using g2_blocks_t = std::set>; + // using g2_blocks_t = std::set>; struct g2_iw_inu_inup_params_t { @@ -21,23 +22,19 @@ namespace pomerol2triqs { /// Channel in which Matsubara frequency representation is defined. channel_t channel = PH; - /// Order of block indices in the definition of G^2. - block_order_t block_order = AABB; - - /// List of block index pairs of G^2 to measure. - /// default: measure all blocks - g2_blocks_t blocks = g2_blocks_t{}; - /// Number of bosonic Matsubara frequencies. - int n_iw = 30; + int n_b; /// Number of fermionic Matsubara frequencies. - int n_inu = 30; + int n_f; + + /// 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) {} }; - +/* struct g2_iw_l_lp_params_t { /// Structure of G^2 blocks. @@ -71,4 +68,5 @@ namespace pomerol2triqs { g2_iw_l_lp_params_t() {} g2_iw_l_lp_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 b9fdc74..91b135f 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -294,6 +294,8 @@ namespace pomerol2triqs { if (!states_class || !matrix_h || !rho || !ops_container) TRIQS_RUNTIME_ERROR << "compute_gf: internal error!"; + if (verbose && !comm.rank()) { std::cout << "\nPomerol: computing GF" << std::endl; } + struct index_visitor { std::vector indices; void operator()(int i) { indices.push_back(std::to_string(i)); } @@ -375,6 +377,46 @@ 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; + + g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f); + for(int ib=0; ib auto pomerol_ed::compute_g2(gf_struct_t const &gf_struct, gf_mesh const &mesh, block_order_t block_order, g2_blocks_t const &g2_blocks, @@ -453,68 +495,65 @@ namespace pomerol2triqs { return make_block2_gf(block_names, block_names, std::move(gf_vecvec)); } -*/ - auto pomerol_ed::G2_iw(g2_iw_inu_inup_params_t const &p) -> g2_t { - // auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf> { + auto pomerol_ed::G2_iw_inu_inup(g2_iw_inu_inup_params_t const &p) -> block2_gf> { if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_inu_inup: no Hamiltonian has been diagonalized"; compute_rho(p.beta); compute_field_operators(p.gf_struct); if (verbose && !comm.rank()) std::cout << "G2_iw_inu_inup: filling output container" << std::endl; - // auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { - // long mesh_index = 0; - // for (auto w_nu_nup : g2_el.mesh()) { - // if ((mesh_index++) % comm.size() != comm.rank()) continue; - // - // if (p.channel == AllFermionic) { - // - // int n1 = std::get<0>(w_nu_nup).index(); - // int n2 = std::get<1>(w_nu_nup).index(); - // int n3 = std::get<2>(w_nu_nup).index(); - // - // if (p.block_order == AABB) - // g2_el[w_nu_nup] = -pom_g2(n2, n1 + n3 - n2, n1); - // else - // g2_el[w_nu_nup] = +pom_g2(n1 + n3 - n2, n2, n1); - // - // } else { // p.channel == PH or PP - // - // int w_n = std::get<0>(w_nu_nup).index(); - // int nu_n = std::get<1>(w_nu_nup).index(); - // int nup_n = std::get<2>(w_nu_nup).index(); - // - // int W_n = p.channel == PH ? w_n + nu_n : w_n - nup_n - 1; - // - // if (p.block_order == AABB) { - // g2_el[w_nu_nup] = -pom_g2(W_n, nup_n, nu_n); - // } else { - // g2_el[w_nu_nup] = +pom_g2(nup_n, W_n, nu_n); - // } - // } - // } - // }; - - // gf_mesh mesh_b{p.beta, Boson, p.n_iw}; - // gf_mesh mesh_f{p.beta, Fermion, p.n_inu}; - // - // gf_mesh mesh_bff{mesh_b, mesh_f, mesh_f}; - // gf_mesh mesh_fff{mesh_f, mesh_f, mesh_f}; - // - // block2_gf> g2; - // - // if (p.channel == AllFermionic) - // g2 = compute_g2(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler); - // else - // g2 = compute_g2(p.gf_struct, mesh_bff, p.block_order, p.blocks, filler); - // - // g2() = mpi_all_reduce(g2(), comm); - // - // return g2; + auto filler = [&p, this](gf_view g2_el, auto const &pom_g2) { + long mesh_index = 0; + for (auto w_nu_nup : g2_el.mesh()) { + if ((mesh_index++) % comm.size() != comm.rank()) continue; + + if (p.channel == AllFermionic) { + + int n1 = std::get<0>(w_nu_nup).index(); + int n2 = std::get<1>(w_nu_nup).index(); + int n3 = std::get<2>(w_nu_nup).index(); + + if (p.block_order == AABB) + g2_el[w_nu_nup] = -pom_g2(n2, n1 + n3 - n2, n1); + else + g2_el[w_nu_nup] = +pom_g2(n1 + n3 - n2, n2, n1); + + } else { // p.channel == PH or PP + + int w_n = std::get<0>(w_nu_nup).index(); + int nu_n = std::get<1>(w_nu_nup).index(); + int nup_n = std::get<2>(w_nu_nup).index(); + + int W_n = p.channel == PH ? w_n + nu_n : w_n - nup_n - 1; + + if (p.block_order == AABB) { + g2_el[w_nu_nup] = -pom_g2(W_n, nup_n, nu_n); + } else { + g2_el[w_nu_nup] = +pom_g2(nup_n, W_n, nu_n); + } + } + } + }; + + gf_mesh mesh_b{p.beta, Boson, p.n_iw}; + gf_mesh mesh_f{p.beta, Fermion, p.n_inu}; + + gf_mesh mesh_bff{mesh_b, mesh_f, mesh_f}; + gf_mesh mesh_fff{mesh_f, mesh_f, mesh_f}; + + block2_gf> g2; + + if (p.channel == AllFermionic) + g2 = compute_g2(p.gf_struct, mesh_fff, p.block_order, p.blocks, filler); + else + g2 = compute_g2(p.gf_struct, mesh_bff, p.block_order, p.blocks, filler); + + g2() = mpi_all_reduce(g2(), comm); + + return g2; } -/* auto pomerol_ed::G2_iw_l_lp(g2_iw_l_lp_params_t const &p) -> block2_gf> { if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw_l_lp: no Hamiltonian has been diagonalized"; compute_rho(p.beta); diff --git a/example/2band.atom.py b/example/2band.atom.py index 0439dda..38f1c30 100644 --- a/example/2band.atom.py +++ b/example/2band.atom.py @@ -49,6 +49,7 @@ # Pomerol: site_label, orbital_index, spin_name index_converter = {(sn, o) : ("loc", o, "down" if sn == "dn" else "up") for sn, o in product(spin_names, orb_names)} +print "index_converter:", index_converter # Make PomerolED solver object ed = PomerolED(index_converter, verbose = True) @@ -86,15 +87,21 @@ # G^{(2)} # ########### -# common_g2_params = {'gf_struct' : gf_struct, -# 'beta' : beta, -# 'blocks' : g2_blocks, -# 'n_iw' : g2_n_iw} +common_g2_params = {'channel' : "PH", + 'gf_struct' : gf_struct, + 'beta' : beta, + 'n_f' : 5, + 'n_b' : 1, } ############################### # 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 ) +print type(G2_iw) +print G2_iw.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", # block_order = "AABB", diff --git a/python/pomerol2triqs_converters.hxx b/python/pomerol2triqs_converters.hxx index 17edfce..c8eae7d 100644 --- a/python/pomerol2triqs_converters.hxx +++ b/python/pomerol2triqs_converters.hxx @@ -1,117 +1,6 @@ // DO NOT EDIT // Generated automatically using libclang using the command : -// c++2py.py ../c++/pomerol_ed.hpp -I../../../pomerol/installed/include -I/usr/include/eigen3 -I../c++ --only_converters -p -mpytriqs.applications.impurity_solvers.pomerol2triqs -o pomerol2triqs --moduledoc "TRIQS wrapper around Pomerol ED library" - - -// --- C++ Python converter for g2_iw_l_lp_params_t -#include -#include -#include - -namespace triqs { namespace py_tools { - -template <> struct py_converter { - static PyObject *c2py(g2_iw_l_lp_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, "block_order", convert_to_python(x.block_order)); - PyDict_SetItemString( d, "blocks" , convert_to_python(x.blocks)); - PyDict_SetItemString( d, "n_iw" , convert_to_python(x.n_iw)); - PyDict_SetItemString( d, "n_l" , convert_to_python(x.n_l)); - PyDict_SetItemString( d, "n_inu_sum" , convert_to_python(x.n_inu_sum)); - PyDict_SetItemString( d, "inu_sum_tol", convert_to_python(x.inu_sum_tol)); - 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_l_lp_params_t py2c(PyObject *dic) { - g2_iw_l_lp_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, "block_order", res.block_order ,AABB); - _get_optional(dic, "blocks" , res.blocks ,g2_blocks_t{}); - _get_optional(dic, "n_iw" , res.n_iw ,30); - _get_optional(dic, "n_l" , res.n_l ,20); - _get_optional(dic, "n_inu_sum" , res.n_inu_sum ,500); - _get_optional(dic, "inu_sum_tol", res.inu_sum_tol ,1e-6); - 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","block_order","blocks","n_iw","n_l","n_inu_sum","inu_sum_tol"}; - 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_optional (dic, fs, err, "block_order", "pomerol2triqs::block_order_t"); - _check_optional (dic, fs, err, "blocks" , "g2_blocks_t"); - _check_optional (dic, fs, err, "n_iw" , "int"); - _check_optional (dic, fs, err, "n_l" , "int"); - _check_optional (dic, fs, err, "n_inu_sum" , "int"); - _check_optional (dic, fs, err, "inu_sum_tol", "double"); - 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_l_lp_params_t\n" < struct py_converter { static PyObject *c2py(g2_iw_inu_inup_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, "block_order", convert_to_python(x.block_order)); - PyDict_SetItemString( d, "blocks" , convert_to_python(x.blocks)); - PyDict_SetItemString( d, "n_iw" , convert_to_python(x.n_iw)); - PyDict_SetItemString( d, "n_inu" , convert_to_python(x.n_inu)); + 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; } @@ -152,11 +43,13 @@ template <> struct py_converter { g2_iw_inu_inup_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, "block_order", res.block_order ,AABB); - _get_optional(dic, "blocks" , res.blocks ,g2_blocks_t{}); - _get_optional(dic, "n_iw" , res.n_iw ,30); - _get_optional(dic, "n_inu" , res.n_inu ,30); + _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; } @@ -187,7 +80,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","block_order","blocks","n_iw","n_inu"}; + 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"; @@ -199,13 +92,15 @@ 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_optional (dic, fs, err, "block_order", "pomerol2triqs::block_order_t"); - _check_optional (dic, fs, err, "blocks" , "g2_blocks_t"); - _check_optional (dic, fs, err, "n_iw" , "int"); - _check_optional (dic, fs, err, "n_inu" , "int"); + _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; diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py index 43c9b0a..81c7025 100644 --- a/python/pomerol2triqs_desc.py +++ b/python/pomerol2triqs_desc.py @@ -30,7 +30,7 @@ module.add_using("namespace Pomerol") module.add_enum("spin", ["down", "up"], "Pomerol", "Spin projection") -module.add_enum("block_order_t", ["AABB", "ABBA"], "pomerol2triqs", "G^{(2)} block order") +# module.add_enum("block_order_t", ["AABB", "ABBA"], "pomerol2triqs", "G^{(2)} block order") module.add_enum("channel_t", ["PP", "PH", "AllFermionic"], "pomerol2triqs", "G^{(2)} channel") # The class pomerol_ed