Merge pull request #6 from mili1247/unstable

add measure_statehist

Author: mili1247

c++/triqs_ctseg/configuration.cpp

@@ -222,6 +222,30 @@ double K_overlap(std::vector<segment_t> const &seglist, tau_t const &tau, bool i
   return is_c ? result : -result;
 }
 
+// ---------------------------
+
+// List of operators containing all colors.
+// Time are ordered in decreasing order, in agreement with the whole physic literature.
+std::vector<colored_ops_t> colored_ordered_ops(std::vector<std::vector<segment_t>> const &seglists) {
+  int c = 0;  // index of color
+  std::vector<colored_ops_t> ops_list; // list of all the operators
+  for (auto const &seglist : seglists) {
+    for (auto const &s : seglist) {
+      ops_list.push_back({s.tau_c, c, false});
+      ops_list.push_back({s.tau_cdag, c, true});
+      if (is_cyclic(s)) {
+        ops_list.push_back({tau_t::beta(), c, false});
+        ops_list.push_back({tau_t::zero(), c, true});
+      }
+    }
+    ++c;
+  }
+  std::sort(ops_list.begin(), ops_list.end(), [](const colored_ops_t &a, const colored_ops_t &b) {
+    return b.tau < a.tau; // Note the order of b and a for descending sort
+  });
+  return ops_list;
+}
+
 // ===================  PRINTING ========================
 
 std::ostream &operator<<(std::ostream &out, std::vector<segment_t> const &sl) {
@@ -248,3 +272,11 @@ std::ostream &operator<<(std::ostream &out, configuration_t const &config) {
   return out;
 }
 
+// ---------------------------
+
+std::ostream &operator<<(std::ostream &out, std::vector<colored_ops_t> const &col) {
+  for (auto const &[i, co] : itertools::enumerate(col))
+    out << "\n" << "i: " << i << ", tau: " << co.tau << ", color: " << co.color 
+        << ", " << (co.is_cdag ? "cdag" : "c");
+  return out;
+}

c++/triqs_ctseg/configuration.hpp

@@ -41,6 +41,13 @@ struct segment_t {
   static segment_t full_line() { return {tau_t::beta(), tau_t::zero()}; }
 };
 
+// operator descriptor: time, color, is_cdag[creation/annihilation]
+struct colored_ops_t {
+  tau_t tau;
+  int color;
+  bool is_cdag;
+};
+
 // simple alias
 using vec_seg_iter_t = std::vector<segment_t>::const_iterator;
 
@@ -141,9 +148,15 @@ double K_overlap(std::vector<segment_t> const &seglist, tau_t const &tau_c, tau_
 double K_overlap(std::vector<segment_t> const &seglist, tau_t const &tau, bool is_c, gf<imtime, matrix_valued> const &K,
                  int c1, int c2);
 
+// List of operators containing all colors.
+std::vector<colored_ops_t> colored_ordered_ops(std::vector<std::vector<segment_t>> const &seglists);
+
 // ===================  PRINTING ========================
 
 std::ostream &operator<<(std::ostream &out, std::vector<segment_t> const &sl);
 
 std::ostream &operator<<(std::ostream &out, configuration_t const &config);
 template <> struct fmt::formatter<configuration_t> : ostream_formatter {};
+
+std::ostream &operator<<(std::ostream &out, std::vector<colored_ops_t> const &col);
+template <> struct fmt::formatter<std::vector<colored_ops_t>> : ostream_formatter {};

c++/triqs_ctseg/measures.hpp

@@ -7,3 +7,4 @@
 #include "./measures/density.hpp"
 #include "./measures/sign.hpp"
 #include "./measures/perturbation_order_histo.hpp"
+#include "./measures/state_hist.hpp"

c++/triqs_ctseg/measures/state_hist.cpp

@@ -0,0 +1,59 @@
+#include "./state_hist.hpp"
+#include "../logs.hpp"
+
+namespace measures {
+
+  state_hist::state_hist(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results)
+     : wdata{wdata}, config{config}, results{results} {
+
+    beta    = p.beta;
+    H       = nda::zeros<double>(ipow(2, config.n_color()));
+  }
+
+  // -------------------------------------
+
+  void state_hist::accumulate(double s) {
+
+    LOG("\n =================== MEASURE <state histogram>  ================ \n");
+
+    /// Measure the state histogram
+    /**
+    *Measures the time the impurity spends in a certain atomic state:
+    * an atomic state is characterized by occupation number for each color, e.g.
+    * $|n_1 n_2 n_3 n_4\rangle = |0 1 1 0\rangle$ for a model with 4 colors
+    *
+    * - the index of a state in the histogram is given by $\sum_i n_i 2^i$
+    *
+    * - the length of the histogram is 2^n_colors
+    */
+    
+    Z += s;
+
+    double tau_prev = beta; // time of prevous operator; start with beta
+    nda::vector<bool> state = nda::zeros<bool>(config.n_color());
+    for (auto const &op : colored_ordered_ops(config.seglists)) {
+        int state_idx = 0;
+        for (auto c : range(config.n_color())) 
+            if (state(c)) state_idx += ipow(2, c); // get the index of the impurity state
+        H(state_idx) += (tau_prev - op.tau);
+        tau_prev = (double)op.tau;
+        ALWAYS_EXPECTS((state(op.color) == op.is_cdag), "Operator error at color {}", op.color);
+        state(op.color) = !op.is_cdag;
+    }
+
+    // get edge state contribution; tau_prev has time of last operator
+    ALWAYS_EXPECTS((state == nda::zeros<bool>(config.n_color())), "Operator error");
+    H(0) += tau_prev;
+  }
+  // -------------------------------------
+
+  void state_hist::collect_results(mpi::communicator const &c) {
+
+    Z = mpi::all_reduce(Z, c);
+    H = mpi::all_reduce(H, c);
+    H = H / (Z * beta);
+
+    // store the result (not reused later, hence we can move it).
+    results.state_hist = std::move(H);
+  }
+} // namespace measures

c++/triqs_ctseg/measures/state_hist.hpp

@@ -0,0 +1,25 @@
+#pragma once
+#include "../configuration.hpp"
+#include "../work_data.hpp"
+#include "../results.hpp"
+#include "../util.hpp"
+
+namespace measures {
+
+  struct state_hist {
+
+    work_data_t const &wdata;
+    configuration_t const &config;
+    results_t &results;
+    double beta;
+
+    nda::vector<double> H;
+
+    double Z = 0;
+
+    state_hist(params_t const &params, work_data_t const &wdata, configuration_t const &config, results_t &results);
+
+    void accumulate(double s);
+    void collect_results(mpi::communicator const &c);
+  };
+} // namespace measures

c++/triqs_ctseg/params.hpp

@@ -108,6 +108,9 @@ struct solve_params_t {
   /// Whether to measure langle s_x(tau)s_x(0)rangle (see [[measure_sperp_tau]])
   bool measure_sperpt = false;
 
+  /// Whether to measure state histograms (see [[measure_statehist]])
+  bool measure_statehist = false;
+
   /// Hartree shift of the chem pot
   nda::vector<double> hartree_shift = nda::vector<double>{};
 

c++/triqs_ctseg/results.cpp

@@ -15,6 +15,7 @@ void h5_write(h5::group h5group, std::string subgroup_name, results_t const &c)
   h5_write(grp, "densities", c.densities);
   h5_write(grp, "perturbation_order_histo_Delta", c.perturbation_order_histo_Delta);
   h5_write(grp, "perturbation_order_histo_Jperp", c.perturbation_order_histo_Jperp);
+  h5_write(grp, "state_hist", c.state_hist);
 }
 
 //------------------------------------
@@ -34,4 +35,5 @@ void h5_read(h5::group h5group, std::string subgroup_name, results_t &c) {
   h5_read(grp, "densities", c.densities);
   h5_read(grp, "perturbation_order_histo_Delta", c.perturbation_order_histo_Delta);
   h5_read(grp, "perturbation_order_histo_Jperp", c.perturbation_order_histo_Jperp);
+  h5_read(grp, "state_hist", c.state_hist);
 }

c++/triqs_ctseg/results.hpp

@@ -38,6 +38,9 @@ struct results_t {
   /// Perturbation order histogram
   std::optional<triqs::stat::histogram> perturbation_order_histo_Jperp;
 
+  /// State histogram
+  std::optional<nda::vector<double>> state_hist;
+
   /// Average sign 
   double sign;
 };

c++/triqs_ctseg/solver_core.cpp

@@ -93,6 +93,7 @@ void solver_core::solve(solve_params_t const &solve_params) {
   if (p.measure_sperpt) CTQMC.add_measure(measures::sperp_tau{p, wdata, config, results}, "<s_x s_x>(tau)");
   if (p.measure_perturbation_order_histograms)
     CTQMC.add_measure(measures::perturbation_order_histo{p, wdata, config, results}, "Perturbation orders");
+  if (p.measure_statehist) CTQMC.add_measure(measures::state_hist{p, wdata, config, results}, "State histograms");
 
   // Run and collect results
   CTQMC.warmup_and_accumulate(p.n_warmup_cycles, p.n_cycles, p.length_cycle,

c++/triqs_ctseg/util.hpp

@@ -26,3 +26,8 @@ long det_lower_bound_x(auto const &d, auto const &x) {
 long det_lower_bound_y(auto const &d, auto const &y) {
   return lower_bound([&d](long i) { return d.get_y(i).first; }, d.size(), y);
 }
+
+// Integer power
+constexpr unsigned int ipow(unsigned int n, unsigned int m) {
+  return m == 0 ? 1 : m == 1 ? n : n * ipow(n, m - 1);
+}