[feat] measure G_iw_direct test, cleanup, expose

* fix indices in direct measurement of G_iw * fill negative frequencies with conj(G_iw) * cleanup * integrate into single site bethe test the new measurement and check against Fourier trafo of G(tau)
TRIQS · Jun 20, 2023 · f8c236e · f8c236e
1 parent b814100
commit f8c236e
Show file tree

Hide file tree

Showing 3 changed files with 30 additions and 8 deletions.
diff --git a/c++/triqs_cthyb/measures/G_iw.cpp b/c++/triqs_cthyb/measures/G_iw.cpp
@@ -43,12 +43,10 @@ namespace triqs_cthyb {
         auto val      = (y.first >= x.first ? s : -s) * M;
         double dtau   = double(y.first - x.first);
         auto const &m = G_iw[0].mesh();
-        long n_iw     = m.last_index();
         auto ex_dt    = std::exp(M_PI / m.domain().beta * dtau * 1i);
         auto ex_2dt   = ex_dt * ex_dt;
         auto ex       = ex_dt;
-        for (long n = 0; n < n_iw; ++n, ex *= ex_2dt) { G_iw[block_idx][n](y.second, x.second) += val * ex; }
-        //this->G_tau[block_idx][closest_mesh_pt(dtau)](y.second, x.second) += val;
+        for (long w = 0; w <= m.last_index(); ++w, ex *= ex_2dt) { G_iw[block_idx][w](y.second, x.second) += val * ex; }
       })
         ;
     }
@@ -58,11 +56,12 @@ namespace triqs_cthyb {
 
     G_iw         = mpi::all_reduce(G_iw, c);
     average_sign = mpi::all_reduce(average_sign, c);
-    G_iw /= -real(average_sign);
-    // for (auto &gbl : G_iw) {
-    //   //double beta = gbl.mesh().domain().beta;
-    //   gbl /= -real(average_sign); //  * beta * gbl.mesh().delta();
-    // }
+    G_iw /= -real(average_sign) * G_iw[0].mesh().domain().beta;
+
+    // fill negative frequencies
+    for (auto block_idx : range(G_iw.size())) {
+      for (long w = 0; w <= G_iw[0].mesh().last_index(); ++w) { G_iw[block_idx][-w - 1] = conj(G_iw[block_idx][w]); }
+    }
   }
 
 } // namespace triqs_cthyb
diff --git a/python/triqs_cthyb/solver_core_desc.py b/python/triqs_cthyb/solver_core_desc.py
@@ -54,6 +54,11 @@
              read_only= True,
              doc = r"""Violation of the fundamental Green function property G(tau)[i,j] = G(tau)*[j,i] after the measurement""")
 
+c.add_member(c_name = "G_iw_direct",
+             c_type = "std::optional<G_iw_t>",
+             read_only= True,
+             doc = r"""Single-particle Green's function :math:`G_iw` in Matsubara Green funcitons.""")
+
 c.add_member(c_name = "G_l",
              c_type = "std::optional<G_l_t>",
              read_only= True,
@@ -420,6 +425,16 @@
              initializer = """ true """,
              doc = r"""Measure G(tau)? :math:`G_{ij}(\tau)=G_{ji}^*(\tau)` is enforced for the resulting G(tau)""")
 
+c.add_member(c_name = "measure_G_iw",
+             c_type = "bool",
+             initializer = """ false """,
+             doc = r"""Measure G_iw directly?""")
+
+c.add_member(c_name = "measure_G_n_iw",
+             c_type = "int",
+             initializer = """ 20 """,
+             doc = r"""Number of Matsubara frequencies for the G_iw measurement.""")
+
 c.add_member(c_name = "measure_G_l",
              c_type = "bool",
              initializer = """ false """,

diff --git a/test/python/single_site_bethe.py b/test/python/single_site_bethe.py
@@ -26,6 +26,8 @@
 p["n_cycles"] = 5000
 p["measure_G_l"] = True
 p["move_double"] = False
+p["measure_G_iw"] = True
+p["measure_G_n_iw"] = 10
 p["perform_tail_fit"] = True
 p["fit_max_moment"] = 3
 p["fit_min_w"] = 1.2
@@ -63,6 +65,7 @@
         Results["G_l"] = S.G_l
 
         Results["G_iw"] = S.G_iw
+        Results["G_iw_direct"] = S.G_iw_direct
         Results["G_iw_raw"] = S.G_iw_raw
 
         # we store Sigma_iw_arw here, but comparing noisy Sigma from Dyson after 2 DMFT
@@ -84,3 +87,8 @@
 
         assert_block_gfs_are_close(Results["Delta_tau"], S.Delta_tau)
 
+# compare direct measured G_iw
+iw0 = len(S.G_iw.mesh)//2
+n_iw = p["measure_G_n_iw"]
+for block, gf in S.G_iw_direct:
+    assert_arrays_are_close(S.G_iw[block].data[iw0-n_iw:iw0+n_iw, :, :], S.G_iw_direct[block].data, precision=1e-4)