Wrap results.nn_tau as a Block2Gf

c++/triqs_ctseg/measures/nn_tau.cpp

@@ -10,9 +10,12 @@ namespace measures {
     ntau    = p.n_tau_chi2;
     dtau    = p.beta / (ntau - 1);
     n_color = config.n_color();
+    block_number = wdata.block_number;
+    index_in_block = wdata.index_in_block;
 
-    q_tau   = gf<imtime>({beta, Boson, p.n_tau_chi2}, {n_color, n_color});
+    q_tau   = gf<imtime>({beta, Boson, ntau}, {n_color, n_color});
     q_tau() = 0;
+    q_tau_block = make_block2_gf<imtime>({beta, Boson, ntau}, p.gf_struct);
   }
 
   // -------------------------------------
@@ -65,7 +68,12 @@ namespace measures {
     q_tau = mpi::all_reduce(q_tau, c);
     q_tau = q_tau / Z; //(beta * Z * q_tau.mesh().delta());
 
-    // store the result (not reused later, hence we can move it).
-    results.nn_tau = std::move(q_tau);
+    // store the result 
+    for (auto const &c1 : range(n_color)) {
+      for (auto const &c2 : range(n_color)) {
+        q_tau_block(block_number[c1], block_number[c2]).data()(range::all, index_in_block[c1], index_in_block[c2]) = q_tau.data()(range::all, c1, c2);
+      }
+    }
+    results.nn_tau = std::move(q_tau_block);
   }
 } // namespace measures

c++/triqs_ctseg/measures/nn_tau.hpp

@@ -2,7 +2,6 @@
 #include "../configuration.hpp"
 #include "../work_data.hpp"
 #include "../results.hpp"
-//#include "../precompute_fprefactor.hpp"
 
 namespace measures {
 
@@ -14,8 +13,10 @@ namespace measures {
     double beta;
     double dtau;
     int ntau;
+    std::vector<long> block_number, index_in_block;
 
     gf<imtime> q_tau;
+    block2_gf<imtime> q_tau_block;
 
     double Z = 0;
     int n_color;

c++/triqs_ctseg/results.hpp

@@ -15,7 +15,7 @@ struct results_t {
   std::optional<block_gf<imtime>> F_tau;
 
   /// Density-density time correlation function :math:`\langle n_a(\tau) n_b(0) \rangle`.
-  std::optional<gf<imtime>> nn_tau;
+  std::optional<block2_gf<imtime>> nn_tau;
 
   /// Perpendicular spin-spin correlation function :math:`\langle s_x(\tau) s_x(0) \rangle`.
   std::optional<gf<imtime>> sperp_tau;

test/c++/anderson.cpp

@@ -70,7 +70,7 @@ TEST(CTSEGJ, Anderson) {
   if (c.rank() == 0) {
     h5::file ref_file("anderson.ref.h5", 'r');
     block_gf<imtime> G_tau, F_tau;
-    gf<imtime> nn_tau;
+    block2_gf<imtime> nn_tau;
     nda::matrix<double> nn_static;
     nda::array<double, 1> densities;
     h5_read(ref_file, "G_tau", G_tau);
@@ -81,7 +81,7 @@ TEST(CTSEGJ, Anderson) {
     EXPECT_ARRAY_NEAR(densities, Solver.results.densities, precision);
     EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision);
     EXPECT_BLOCK_GF_NEAR(F_tau, Solver.results.F_tau.value(), precision);
-    EXPECT_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
+    EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
     EXPECT_ARRAY_NEAR(nn_static, Solver.results.nn_static.value(), precision);
   }
 }

test/c++/dynamical_U.cpp

@@ -81,7 +81,7 @@ TEST(CTSEGJ, Dynamical_U) {
   if (c.rank() == 0) {
     h5::file ref_file("dynamical_U.ref.h5", 'r');
     block_gf<imtime> G_tau, F_tau;
-    gf<imtime> nn_tau;
+    block2_gf<imtime> nn_tau;
     nda::matrix<double> nn_static;
     nda::array<double, 1> densities;
     h5_read(ref_file, "G_tau", G_tau);
@@ -92,7 +92,7 @@ TEST(CTSEGJ, Dynamical_U) {
     EXPECT_ARRAY_NEAR(densities, Solver.results.densities, precision);
     EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision);
     EXPECT_BLOCK_GF_NEAR(F_tau, Solver.results.F_tau.value(), precision);
-    EXPECT_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
+    EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
     EXPECT_ARRAY_NEAR(nn_static, Solver.results.nn_static.value(), precision);
   }
 }

test/c++/jperp.cpp

@@ -76,7 +76,7 @@ TEST(CTSEGJ, J_perp) {
   if (c.rank() == 0) {
     h5::file ref_file("jperp.ref.h5", 'r');
     block_gf<imtime> G_tau;
-    gf<imtime> nn_tau;
+    block2_gf<imtime> nn_tau;
     nda::matrix<double> nn_static;
     nda::array<double, 1> densities;
     h5_read(ref_file, "G_tau", G_tau);
@@ -85,7 +85,7 @@ TEST(CTSEGJ, J_perp) {
     h5_read(ref_file, "densities", densities);
     EXPECT_ARRAY_NEAR(densities, Solver.results.densities, precision);
     EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision);
-    EXPECT_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
+    EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
     EXPECT_ARRAY_NEAR(nn_static, Solver.results.nn_static.value(), precision);
   }
 }

test/c++/spin_spin.cpp

@@ -85,7 +85,7 @@ TEST(CTSEGJ, Spin_Spin) {
   if (c.rank() == 0) {
     h5::file ref_file("spin_spin.ref.h5", 'r');
     block_gf<imtime> G_tau, F_tau;
-    gf<imtime> nn_tau;
+    block2_gf<imtime> nn_tau;
     nda::matrix<double> nn_static;
     nda::array<double, 1> densities;
     h5_read(ref_file, "G_tau", G_tau);
@@ -96,7 +96,7 @@ TEST(CTSEGJ, Spin_Spin) {
     EXPECT_ARRAY_NEAR(densities, Solver.results.densities, precision);
     EXPECT_BLOCK_GF_NEAR(G_tau, Solver.results.G_tau, precision);
     EXPECT_BLOCK_GF_NEAR(F_tau, Solver.results.F_tau.value(), precision);
-    EXPECT_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
+    EXPECT_BLOCK2_GF_NEAR(nn_tau, Solver.results.nn_tau.value(), precision);
     EXPECT_ARRAY_NEAR(nn_static, Solver.results.nn_static.value(), precision);
   }
 }

test/python/multiorb.py

@@ -24,8 +24,9 @@
 n_tau_bosonic = 1001
 
 # Solver construction parameters
+gf_struct = [(f'up{(l+2)//2}', 1) if l % 2 == 0 else (f'down{(l+1)//2}', 1) for l in range(2 * n_orb)]
 constr_params = {
-    "gf_struct": [(f'up{(l+2)//2}', 1) if l % 2 == 0 else (f'down{(l+1)//2}', 1) for l in range(2 * n_orb)],
+    "gf_struct": gf_struct,
     "beta": beta,
     "n_tau": n_tau,
     "n_tau_bosonic": n_tau_bosonic
@@ -69,12 +70,18 @@
 # Solve
 S.solve(**solve_params)
 
+# Unwrap Block2Gf nn_tau into matrix Gf
+nn_tau = GfImTime(indices=range(2*n_orb), beta=beta, statistic = "Boson", n_points=n_tau_bosonic)
+for i, (bl_name_i, si) in enumerate(gf_struct): 
+    for j, (bl_name_j, sj) in enumerate(gf_struct): 
+        nn_tau[i, j] = S.results.nn_tau[bl_name_i, bl_name_j][0, 0]
+
 # Save and compare to reference
 if mpi.is_master_node():
     with HDFArchive("multiorb.out.h5", 'w') as A:
         A['G_tau'] = S.results.G_tau
         A['F_tau'] = S.results.F_tau
-        A['nn_tau'] = S.results.nn_tau
+        A["nn_tau"] = nn_tau
         A['nn'] = S.results.nn_static
         A['densities'] = S.results.densities
 

test/python/spin_spin.py

@@ -17,8 +17,9 @@
 n_tau_bosonic = 2001
 
 # Solver construction parameters
+gf_struct = [('down', 1), ('up', 1)]
 constr_params = {
-    "gf_struct": [('down', 1), ('up', 1)],
+    "gf_struct": gf_struct,
     "beta": beta,
     "n_tau": n_tau,
     "n_tau_bosonic": n_tau_bosonic
@@ -62,12 +63,20 @@
 # Solve
 S.solve(**solve_params)
 
+# Unwrap Block2Gf nn_tau into matrix Gf
+n_color = 2
+nn_tau = GfImTime(indices=[0, 1], beta=beta, statistic = "Boson", n_points=n_tau_bosonic)
+nn_tau[0, 0] = S.results.nn_tau["down", "down"][0, 0]
+nn_tau[1, 0] = S.results.nn_tau["up", "down"][0, 0]
+nn_tau[0, 1] = S.results.nn_tau["down", "up"][0, 0]
+nn_tau[1, 1] = S.results.nn_tau["up", "up"][0, 0]
+
 # Save and compare to reference
 if mpi.is_master_node():
     with h5.HDFArchive("spin_spin.out.h5", 'w') as A:
         A['G_tau'] = S.results.G_tau
         A['F_tau'] = S.results.F_tau
-        A['nn_tau'] = S.results.nn_tau
+        A['nn_tau'] = nn_tau
         A['nn'] = S.results.nn_static
         A['densities'] = S.results.densities
 

test/python/two_orbitals.py

@@ -57,7 +57,8 @@
 # Save output
 if mpi.is_master_node(): 
     with HDFArchive("two_orbitals_4b.out.h5", "a") as A: 
-        A['nn_tau'] = S.results.nn_tau
+        A['nn_tau_0'] = S.results.nn_tau["up1", "up1"][0, 0]
+        A['nn_tau_1'] = S.results.nn_tau["dn1", "dn1"][0, 0]
         A['nn'] = S.results.nn_static
         A['densities'] = S.results.densities
 
@@ -90,7 +91,8 @@
 # Save output
 if mpi.is_master_node(): 
     with HDFArchive("two_orbitals_2b.out.h5", "a") as A: 
-        A['nn_tau'] = S.results.nn_tau
+        A['nn_tau_0'] = S.results.nn_tau["up1", "up1"][0, 0]
+        A['nn_tau_1'] = S.results.nn_tau["up1", "up1"][1, 1]
         A['nn'] = S.results.nn_static
         A['densities'] = S.results.densities