Change interface of G2 calc: Pass sets of four operator indices as a …

…vector

c++/g2_parameters.hpp

@@ -13,7 +13,7 @@ namespace pomerol2triqs {
 
   // using g2_blocks_t = std::set<std::pair<std::string, std::string>>;
 
-  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<four_indices_t> 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<four_indices_t> 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) {}
   };
 
 /*

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<g2_iw_freq_box_t> 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<g2_iw_freq_box_t> {
 
     // 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.size(); i++){
-      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_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<g2_iw_freq_vec_t> 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<g2_iw_freq_box_t> 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.size(); i++){
+        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];
+      }
+      vec_g2_freq_box.push_back(g2);
     }
 
-    return g2;
+    return vec_g2_freq_box;
   }
 
 

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_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

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
 

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)

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 <triqs/python_tools/converters/vector.hpp>
+#include <triqs/python_tools/converters/string.hpp>
+#include <algorithm>
+
+namespace triqs { namespace py_tools {
+
+template <> struct py_converter<g2_iw_legacy_params_t> {
+ 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 <typename T, typename U> 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<T>(PyDict_GetItemString(dic, name));
+  else
+   r = init_default;
+ }
+
+ template <typename T> static void _get_optional(PyObject *dic, const char *name, T &r) {
+  if (PyDict_Contains(dic, pyref::string(name)))
+   r = convert_from_python<T>(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<gf_struct_t>(PyDict_GetItemString(dic, "gf_struct"));
+  res.beta = convert_from_python<double>(PyDict_GetItemString(dic, "beta"));
+  _get_optional(dic, "channel"  , res.channel     ,PH);
+  res.n_b = convert_from_python<int>(PyDict_GetItemString(dic, "n_b"));
+  res.n_f = convert_from_python<int>(PyDict_GetItemString(dic, "n_f"));
+  res.index1 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index1"));
+  res.index2 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index2"));
+  res.index3 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index3"));
+  res.index4 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index4"));
+  return res;
+ }
+
+ template <typename T>
+ static void _check(PyObject *dic, std::stringstream &fs, int &err, const char *name, const char *tname) {
+  if (!convertible_from_python<T>(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 <typename T>
+ 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<T>(dic,fs,err,name,tname);
+ }
+
+ template <typename T>
+ 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<T>(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<std::string> ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"};
+  pyref keys = PyDict_Keys(dic);
+  if (!convertible_from_python<std::vector<std::string>>(keys, true)) {
+   fs << "\nThe dict keys are not strings";
+   goto _error;
+  }
+  ks = convert_from_python<std::vector<std::string>>(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<gf_struct_t             >(dic, fs, err, "gf_struct", "gf_struct_t");
+  _check_mandatory<double                  >(dic, fs, err, "beta"     , "double");
+  _check_optional <pomerol2triqs::channel_t>(dic, fs, err, "channel"  , "pomerol2triqs::channel_t");
+  _check_mandatory<int                     >(dic, fs, err, "n_b"      , "int");
+  _check_mandatory<int                     >(dic, fs, err, "n_f"      , "int");
+  _check_mandatory<indices_t               >(dic, fs, err, "index1"   , "indices_t");
+  _check_mandatory<indices_t               >(dic, fs, err, "index2"   , "indices_t");
+  _check_mandatory<indices_t               >(dic, fs, err, "index3"   , "indices_t");
+  _check_mandatory<indices_t               >(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" <<fs.str();
+   if (raise_exception) PyErr_SetString(PyExc_TypeError, fs2.str().c_str());
+  return false;
+ }
+};
+
+}}
+
+
 // --- C++ Python converter for g2_iw_freq_vec_params_t
 #include <triqs/python_tools/converters/vector.hpp>
 #include <triqs/python_tools/converters/string.hpp>
@@ -112,15 +223,12 @@ namespace triqs { namespace py_tools {
 template <> struct py_converter<g2_iw_freq_box_params_t> {
  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> {
   g2_iw_freq_box_params_t res;
   res.gf_struct = convert_from_python<gf_struct_t>(PyDict_GetItemString(dic, "gf_struct"));
   res.beta = convert_from_python<double>(PyDict_GetItemString(dic, "beta"));
-  _get_optional(dic, "channel"  , res.channel     ,PH);
+  _get_optional(dic, "channel"         , res.channel            ,PH);
   res.n_b = convert_from_python<int>(PyDict_GetItemString(dic, "n_b"));
   res.n_f = convert_from_python<int>(PyDict_GetItemString(dic, "n_f"));
-  res.index1 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index1"));
-  res.index2 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index2"));
-  res.index3 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index3"));
-  res.index4 = convert_from_python<indices_t>(PyDict_GetItemString(dic, "index4"));
+  res.vec_four_indices = convert_from_python<std::vector<four_indices_t>>(PyDict_GetItemString(dic, "vec_four_indices"));
   return res;
  }
 
@@ -179,7 +284,7 @@ template <> struct py_converter<g2_iw_freq_box_params_t> {
   std::stringstream fs, fs2; int err=0;
 
 #ifndef TRIQS_ALLOW_UNUSED_PARAMETERS
-  std::vector<std::string> ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","index1","index2","index3","index4"};
+  std::vector<std::string> ks, all_keys = {"gf_struct","beta","channel","n_b","n_f","vec_four_indices"};
   pyref keys = PyDict_Keys(dic);
   if (!convertible_from_python<std::vector<std::string>>(keys, true)) {
    fs << "\nThe dict keys are not strings";
@@ -191,15 +296,12 @@ template <> struct py_converter<g2_iw_freq_box_params_t> {
     fs << "\n"<< ++err << " The parameter '" << k << "' is not recognized.";
 #endif
 
-  _check_mandatory<gf_struct_t             >(dic, fs, err, "gf_struct", "gf_struct_t");
-  _check_mandatory<double                  >(dic, fs, err, "beta"     , "double");
-  _check_optional <pomerol2triqs::channel_t>(dic, fs, err, "channel"  , "pomerol2triqs::channel_t");
-  _check_mandatory<int                     >(dic, fs, err, "n_b"      , "int");
-  _check_mandatory<int                     >(dic, fs, err, "n_f"      , "int");
-  _check_mandatory<indices_t               >(dic, fs, err, "index1"   , "indices_t");
-  _check_mandatory<indices_t               >(dic, fs, err, "index2"   , "indices_t");
-  _check_mandatory<indices_t               >(dic, fs, err, "index3"   , "indices_t");
-  _check_mandatory<indices_t               >(dic, fs, err, "index4"   , "indices_t");
+  _check_mandatory<gf_struct_t                >(dic, fs, err, "gf_struct"       , "gf_struct_t");
+  _check_mandatory<double                     >(dic, fs, err, "beta"            , "double");
+  _check_optional <pomerol2triqs::channel_t   >(dic, fs, err, "channel"         , "pomerol2triqs::channel_t");
+  _check_mandatory<int                        >(dic, fs, err, "n_b"             , "int");
+  _check_mandatory<int                        >(dic, fs, err, "n_f"             , "int");
+  _check_mandatory<std::vector<four_indices_t>>(dic, fs, err, "vec_four_indices", "std::vector<four_indices_t>");
   if (err) goto _error;
   return true;
 

python/pomerol2triqs_desc.py

@@ -67,10 +67,13 @@
 c.add_method("""block_gf<refreq> G_w (gf_struct_t gf_struct, double beta, std::pair<double,double> energy_window, int n_w, double im_shift = 0)""",
              doc = r"""Retarded Green's function on real energy axis""")
 
-c.add_method("""triqs::arrays::array<std::complex<double>, 3> G2_iw_freq_box (**pomerol2triqs::g2_iw_freq_box_params_t)""",
+c.add_method("""triqs::arrays::array<std::complex<double>, 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<std::vector<std::complex<double> > > G2_iw_freqs_vec (**pomerol2triqs::g2_iw_freq_vec_params_t)""",
+c.add_method("""std::vector< triqs::arrays::array<std::complex<double>, 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<std::complex<double> > > 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)