Merge pull request #1 from j-otsuki/specify_three_freqs

New interface: Enable to specify arbitrary combinations of three freqs

c++/g2_parameters.hpp

@@ -3,6 +3,8 @@
 #include <triqs/hilbert_space/fundamental_operator_set.hpp>
 using triqs::hilbert_space::gf_struct_t;
 using indices_t = triqs::hilbert_space::fundamental_operator_set::indices_t;
+using four_indices_t = std::tuple<indices_t, indices_t, indices_t, indices_t>;
+using three_freqs_t = std::tuple<int, int, int>;
 
 namespace pomerol2triqs {
 
@@ -11,29 +13,67 @@ namespace pomerol2triqs {
 
   // using g2_blocks_t = std::set<std::pair<std::string, std::string>>;
 
-  struct g2_iw_inu_inup_params_t {
+  struct g2_iw_legacy_params_t {
 
     /// Block structure of GF
     gf_struct_t gf_struct;
 
     /// 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 Matsubara frequencies.
+    /// Number of non-negative bosonic Matsubara frequencies
     int n_b;
-
-    /// Number of fermionic Matsubara frequencies.
+    /// 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;
 
-    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_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 non-negative bosonic Matsubara frequencies
+    int n_b;
+    /// Number of positive fermionic Matsubara frequencies
+    int n_f;
+
+    /// four operator indices in TRIQS convention: (block_name, inner_index)*4
+    std::vector<four_indices_t> four_indices;
   };
+
+  struct g2_iw_freq_fix_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).
+    std::vector<three_freqs_t> three_freqs;
+
+    /// four operator indices in TRIQS convention: (block_name, inner_index)*4
+    std::vector<four_indices_t> four_indices;
+  };
+
 /*
   struct g2_iw_l_lp_params_t {
 

c++/pomerol_ed.cpp

@@ -398,17 +398,19 @@ namespace pomerol2triqs {
   // Two-particle Green's function //
   ///////////////////////////////////
 
-  auto pomerol_ed::G2_iw(g2_iw_inu_inup_params_t const &p) -> g2_t {
+  // core function for computing G2
+  triqs::arrays::array<std::complex<double>, 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<three_freqs_t> const &three_freqs){
+
     if (!matrix_h) TRIQS_RUNTIME_ERROR << "G2_iw: no Hamiltonian has been diagonalized";
-    compute_rho(p.beta);
-    compute_field_operators(p.gf_struct);
+    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(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::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),
@@ -419,54 +421,105 @@ namespace pomerol2triqs {
     pom_g2.prepare();
     pom_g2.compute(false, {}, comm);
 
-    // indices of fermionic Matsubara frequencies [-n_f:n_f)
-    std::vector<int> 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<int> 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<index_wb.size(); ib++)
-    //   for(int if1=0; if1<index_wf.size(); if1++)
-    //     for(int if2=0; if2<index_wf.size(); if2++)
-    //       g2(ib, if1, if2) = -pom_g2(index_wb[ib]+index_wf[if1], index_wf[if2], index_wf[if1]);
-
-    // create a list of three frequency-indices
-    std::vector< std::tuple<int, int, int> > three_freqs;
-    for(int ib=0; ib<index_wb.size(); ib++)
-      for(int if1=0; if1<index_wf.size(); if1++)
-        for(int if2=0; if2<index_wf.size(); if2++)
-          three_freqs.push_back( std::make_tuple(ib, if1, if2) );
-    // std::cout << three_freqs.size() << std::endl;
-
     // compute g2 value using MPI
-    g2_t g2(p.n_b, 2*p.n_f, 2*p.n_f);
+    g2_iw_freq_fix_t g2( three_freqs.size() );
     for(int i=comm.rank(); i<three_freqs.size(); i+=comm.size()){
-      int ib = std::get<0>(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>(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<three_freqs.size(); i++){
-      int ib = std::get<0>(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);
+      boost::mpi::broadcast(comm, g2(i), sender);
     }
 
     // std::cout << "End freq loop: rank" << comm.rank() << std::endl;
     return g2;
   }
 
+
+  std::vector< triqs::arrays::array<std::complex<double>, 1> > pomerol_ed::compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &four_indices, std::vector<three_freqs_t> const &three_freqs){
+
+    std::vector<g2_iw_freq_fix_t> vec_g2;
+    // TODO: MPI
+    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;
+  }
+
+
+  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 ind4 = std::make_tuple(p.index1, p.index2, p.index3, p.index4);
+    p2.four_indices.push_back( ind4 );
+
+    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)
+    std::vector<three_freqs_t> three_freqs;
+    std::vector< std::tuple<int, int, int> > three_indices;  // (ib, if1, if2)
+    {
+      // indices of fermionic Matsubara frequencies [-n_f:n_f)
+      std::vector<int> 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<int> 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<index_wb.size(); ib++)
+        for(int if1=0; if1<index_wf.size(); if1++)
+          for(int if2=0; if2<index_wf.size(); if2++){
+            three_freqs.push_back( std::make_tuple(index_wb[ib], index_wf[if1], index_wf[if2]) );
+            three_indices.push_back( std::make_tuple(ib, if1, if2) );
+          }
+      // std::cout << three_freqs.size() << std::endl;
+    }
+
+    // compute g2 values
+    std::vector<g2_iw_freq_fix_t> 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<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 vec_g2_freq_box;
+  }
+
+
+  auto pomerol_ed::G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p) -> std::vector<g2_iw_freq_fix_t> {
+    return compute_g2(p.gf_struct, p.beta, p.channel, p.four_indices, p.three_freqs);
+  }
+
 /*
   template <typename Mesh, typename Filler>
   auto pomerol_ed::compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order, g2_blocks_t const &g2_blocks,

c++/pomerol_ed.hpp

@@ -82,11 +82,16 @@ namespace pomerol2triqs {
 
     // using w_nu_nup_t = cartesian_product<imfreq, imfreq, imfreq>;
     // using w_l_lp_t   = cartesian_product<imfreq, legendre, legendre>;
-    using g2_t = triqs::arrays::array<std::complex<double>, 3>;
+    using g2_iw_freq_box_t = triqs::arrays::array<std::complex<double>, 3>;
+    using g2_iw_freq_fix_t = triqs::arrays::array<std::complex<double>, 1>;
     // template <typename Mesh, typename Filler>
     // block2_gf<Mesh, tensor_valued<4>> compute_g2(gf_struct_t const &gf_struct, gf_mesh<Mesh> const &mesh, block_order_t block_order,
     //                                              g2_blocks_t const &g2_blocks, Filler filler) const;
 
+    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<three_freqs_t> const &three_freqs);
+
+    std::vector<g2_iw_freq_fix_t> compute_g2(gf_struct_t const &gf_struct, double beta, channel_t channel, std::vector<four_indices_t> const &four_indices, std::vector<three_freqs_t> const &three_freqs);
+
     double density_matrix_cutoff = 1e-15;
 
     public:
@@ -117,9 +122,17 @@ namespace pomerol2triqs {
     /// Retarded Green's function on real energy axis
     block_gf<refreq> G_w(gf_struct_t const &gf_struct, double beta, std::pair<double, double> const &energy_window, int n_w, double im_shift = 0);
 
-    /// Two-particle Green's function, Matsubara frequencies
+    /// 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, in a low-frequency box with 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, for fixed frequencies (wb, wf1, wf2)
     TRIQS_WRAP_ARG_AS_DICT
-    g2_t G2_iw(g2_iw_inu_inup_params_t const &p);
+    std::vector<g2_iw_freq_fix_t> G2_iw_freq_fix(g2_iw_freq_fix_params_t const &p);
 
     /// Two-particle Green's function, Matsubara frequencies
     // TRIQS_WRAP_ARG_AS_DICT

example/2band.atom.py

@@ -95,10 +95,13 @@
 # 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
+# G2_iw = ed.G2_iw_legacy( index1=('up',0), index2=('dn',0), index3=('dn',1), index4=('up',1), **common_g2_params )
 
+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",
@@ -125,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]

example/slater.py

@@ -137,18 +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( 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 )
+
+# 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