[doc] documentation and cleanup

README.md

@@ -1,12 +1,12 @@
-[![build](https://github.com/TRIQS/cthyb/workflows/build/badge.svg)](https://github.com/TRIQS/cthyb/actions?query=workflow%3Abuild)
+[![build](https://github.com/TRIQS/ctseg_J/workflows/build/badge.svg)](https://github.com/TRIQS/ctseg_J/actions?query=workflow%3Abuild)
 
-## TRIQS segment picture impurity solver with spin-spin interactions (CTSEG-J)
+## TRIQS segment picture solver with spin-spin interactions (CTSEG-J)
 
-Copyright (C) 2024, N. Kavokine, H. Lu and O. Parcollet \
+Copyright (C) 2024, N. Kavokine, H. Lu, N. Wentzell and O. Parcollet.
 
 ### Documentation
 
-For documentation and installation instructions see [triqs.github.io/ctseg-j](https://triqs.github.io/ctseg-j).
+For documentation and installation instructions see [triqs.github.io/ctseg-J](https://triqs.github.io/ctseg-J).
 
 ### License
 

c++/triqs_ctseg/configuration.hpp

@@ -41,13 +41,6 @@ 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;
 
@@ -59,6 +52,14 @@ struct jperp_line_t {
   tau_t tau_Sminus, tau_Splus; // times of the S-, S+
 };
 
+// ----------------- Operator -------------------
+// (used by measure state_hist)
+struct colored_ops_t {
+  tau_t tau;
+  int color;
+  bool is_cdag;
+};
+
 // --------------- Configuration ----------------------
 // The configuration is a list of of segments for each color, 
 // and a list of J_perp lines.

c++/triqs_ctseg/dets.hpp

@@ -7,15 +7,13 @@
 using namespace triqs::gfs;
 using namespace triqs::mesh;
 
-// -----------     dets    ----------------------------------------
-
-/// A lambda to adapt the Delta function for the call of the det.
+/// A lambda to adapt Delta(tau) for the call by det_manip.
 struct delta_block_adaptor {
-  gf<imtime, matrix_real_valued> delta; // make a copy. Needed in the real case anyway.
+  gf<imtime, matrix_real_valued> delta; 
 
   double operator()(std::pair<tau_t, int> const &x, std::pair<tau_t, int> const &y) const {
     double res = delta(double(x.first - y.first))(x.second, y.second);
-    return (x.first >= y.first ? res : -res); // x,y first are time_pt, wrapping is automatic in
+    return (x.first >= y.first ? res : -res); // x,y first are tau_t, wrapping is automatic in
                                               // the - operation, but need to compute the sign
   }
 };

c++/triqs_ctseg/measures/nn_static.cpp

@@ -1,8 +1,6 @@
 #include "./nn_static.hpp"
 #include "../logs.hpp"
 
-// TODO : add non regression test.
-
 namespace measures {
 
   nn_static::nn_static(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results)
@@ -21,10 +19,6 @@ namespace measures {
 
     Z += s;
 
-    // FIXME : this is a simple O(n_color^2) algo
-    // If profiling says it is too slow,
-    // then improve by progressing on the 2 seglists in //
-    // but it is a bit more complex to write
     for (int a = 0; a < n_color; ++a)
       for (int b = 0; b < n_color; ++b) {
         for (auto const &sa : config.seglists[a]) {

c++/triqs_ctseg/measures/nn_static.hpp

@@ -2,7 +2,6 @@
 #include "../configuration.hpp"
 #include "../work_data.hpp"
 #include "../results.hpp"
-//#include "../precompute_fprefactor.hpp"
 
 namespace measures {
 

c++/triqs_ctseg/measures/perturbation_order_histo.cpp

@@ -3,12 +3,12 @@
 
 namespace measures {
 
-  perturbation_order_histo::perturbation_order_histo(params_t const &, work_data_t const &wdata,
+  perturbation_order_histo::perturbation_order_histo(params_t const &p, work_data_t const &wdata,
                                                      configuration_t const &config, results_t &results)
      : wdata{wdata}, config{config}, results{results} {
 
-    histo_delta = triqs::stat::histogram{0, 1000}; // FIXME : put as params ?
-    histo_Jperp = triqs::stat::histogram{0, 1000}; // FIXME : put as params ?
+    histo_delta = triqs::stat::histogram{0, p.histogram_max_order}; 
+    histo_Jperp = triqs::stat::histogram{0, p.histogram_max_order}; 
   }
 
   // -------------------------------------

c++/triqs_ctseg/measures/state_hist.cpp

@@ -14,7 +14,7 @@ namespace measures {
 
   void state_hist::accumulate(double s) {
 
-    LOG("\n =================== MEASURE <state histogram>  ================ \n");
+    LOG("\n =================== MEASURE STATE HISTOGRAM  ================ \n");
 
     /// Measure the state histogram
     /**

c++/triqs_ctseg/params.cpp

@@ -29,6 +29,7 @@ void h5_write(h5::group h5group, std::string subgroup_name, solve_params_t const
   h5::group grp = subgroup_name.empty() ? h5group : h5group.create_group(subgroup_name);
 
   h5_write(grp, "h_int", c.h_int);
+  h5_write(grp, "hartree_shift", c.hartree_shift);
   h5_write(grp, "n_cycles", c.n_cycles);
   h5_write(grp, "length_cycle", c.length_cycle);
   h5_write(grp, "n_warmup_cycles", c.n_warmup_cycles);
@@ -55,12 +56,12 @@ void h5_write(h5::group h5group, std::string subgroup_name, solve_params_t const
   h5_write(grp, "measure_nnt", c.measure_nnt);
   h5_write(grp, "measure_sperpt", c.measure_sperpt);
   h5_write(grp, "measure_statehist", c.measure_statehist);
-  h5_write(grp, "hartree_shift", c.hartree_shift);
   h5_write(grp, "det_init_size", c.det_init_size);
   h5_write(grp, "det_n_operations_before_check", c.det_n_operations_before_check);
   h5_write(grp, "det_precision_warning", c.det_precision_warning);
   h5_write(grp, "det_precision_error", c.det_precision_error);
   h5_write(grp, "det_singular_threshold", c.det_singular_threshold);
+  h5_write(grp, "histogram_max_order", c.histogram_max_order);
 
 }
 
@@ -71,6 +72,7 @@ void h5_read(h5::group h5group, std::string subgroup_name, solve_params_t &c) {
   h5::group grp = subgroup_name.empty() ? h5group : h5group.open_group(subgroup_name);
 
   h5_read(grp, "h_int", c.h_int);
+  h5_read(grp, "hartree_shift", c.hartree_shift);
   h5_read(grp, "n_cycles", c.n_cycles);
   h5_read(grp, "length_cycle", c.length_cycle);
   h5_read(grp, "n_warmup_cycles", c.n_warmup_cycles);
@@ -97,12 +99,12 @@ void h5_read(h5::group h5group, std::string subgroup_name, solve_params_t &c) {
   h5_read(grp, "measure_nnt", c.measure_nnt);
   h5_read(grp, "measure_sperpt", c.measure_sperpt);
   h5_read(grp, "measure_statehist", c.measure_statehist);
-  h5_read(grp, "hartree_shift", c.hartree_shift);
   h5_read(grp, "det_init_size", c.det_init_size);
   h5_read(grp, "det_n_operations_before_check", c.det_n_operations_before_check);
   h5_read(grp, "det_precision_warning", c.det_precision_warning);
   h5_read(grp, "det_precision_error", c.det_precision_error);
   h5_read(grp, "det_singular_threshold", c.det_singular_threshold);
+  h5_read(grp, "histogram_max_order", c.histogram_max_order);
 
 }
 

c++/triqs_ctseg/params.hpp

@@ -29,6 +29,9 @@ struct solve_params_t {
   /// Local Hamiltonian
   triqs::operators::many_body_operator h_int;
 
+  /// Chemical potential (high frequency limit of :math:`G_0^{-1}(i\omega) - i \omega`)
+  nda::vector<double> hartree_shift = nda::vector<double>{};
+
   /// Number of QMC cycles
   int n_cycles;
 
@@ -111,8 +114,7 @@ struct solve_params_t {
   /// Whether to measure state histograms (see [[measure_statehist]])
   bool measure_statehist = false;
 
-  /// Chemical potential (high frequency limit of :math:`G_0^{-1}(i\omega) - i \omega`)
-  nda::vector<double> hartree_shift = nda::vector<double>{};
+  // -------- Misc parameters --------------
 
   /// The maximum size of the determinant matrix before a resize
   int det_init_size = 100;
@@ -128,6 +130,9 @@ struct solve_params_t {
 
   /// Bound for the determinant matrix being singular, abs(det) > singular_threshold. If <0, it is !isnormal(abs(det))
   double det_singular_threshold = -1;
+
+  /// Maximum order for the perturbation order histograms 
+  int histogram_max_order = 1000; 
 };
 
 /// A struct combining both constr_params_t and solve_params_t

c++/triqs_ctseg/solver_core.cpp

@@ -90,9 +90,9 @@ void solver_core::solve(solve_params_t const &solve_params) {
   if (p.measure_gt) CTQMC.add_measure(measures::g_f_tau{p, wdata, config, results}, "G(tau)");
   if (p.measure_n) CTQMC.add_measure(measures::density{p, wdata, config, results}, "Density");
   if (p.measure_sign) CTQMC.add_measure(measures::sign{p, wdata, config, results}, "Sign");
-  if (p.measure_nn) CTQMC.add_measure(measures::nn_static{p, wdata, config, results}, "nn(0)");
-  if (p.measure_nnt) CTQMC.add_measure(measures::nn_tau{p, wdata, config, results}, "nn(tau)");
-  if (p.measure_sperpt) CTQMC.add_measure(measures::sperp_tau{p, wdata, config, results}, "<s_x s_x>(tau)");
+  if (p.measure_nn) CTQMC.add_measure(measures::nn_static{p, wdata, config, results}, "<nn>");
+  if (p.measure_nnt) CTQMC.add_measure(measures::nn_tau{p, wdata, config, results}, "<n(tau)n(0)>");
+  if (p.measure_sperpt) CTQMC.add_measure(measures::sperp_tau{p, wdata, config, results}, "<s_x(tau)s_x(0)>");
   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");

c++/triqs_ctseg/work_data.cpp

@@ -131,24 +131,25 @@ work_data_t::work_data_t(params_t const &p, inputs_t const &inputs, mpi::communi
   // Is there a non-zero Delta(tau)?
   for (auto const &bl : range(inputs.delta.size())) {
     if (max_element(abs(inputs.delta[bl].data())) > 1.e-13) has_delta = true;
+    // Report if Delta(tau) has imaginary part. 
+    if (!is_gf_real(inputs.delta[bl], 1e-10)) {
+      if (c.rank() == 0) {
+        spdlog::info("WARNING: The Delta(tau) block number {} is not real in tau space", bl);
+        spdlog::info("WARNING: max(Im[Delta(tau)]) = {}", max_element(abs(imag(inputs.delta[bl].data()))));
+        spdlog::info("WARNING: Disregarding the imaginary component in the calculation.");
+      }
+    }
   }
   if (not has_delta) {
     ALWAYS_EXPECTS(has_jperp, "Error : both J_perp(tau) and Delta(tau) are 0: there is nothing to expand.");
     if (c.rank() == 0) { spdlog::info("Delta(tau) is 0, running only spin moves."); }
   }
 
-  // FIXME: real??
+  // Take the real part of Delta(tau)
   delta = map([](gf_const_view<imtime> d) { return real(d); }, inputs.delta);
   for (auto const &bl : range(delta.size())) {
-    if (!is_gf_real(delta[bl], 1e-10)) {
-      if (c.rank() == 0) {
-        spdlog::info("WARNING: The Delta(tau) block number {} is not real in tau space", bl);
-        spdlog::info("WARNING: max(Im[Delta(tau)]) = {}", max_element(abs(imag(delta[bl].data()))));
-        spdlog::info("WARNING: Disregarding the imaginary component in the calculation.");
-      }
-    }
     // Construct the detmanip object for block bl
-    dets.emplace_back(delta_block_adaptor{real(delta[bl])}, p.det_init_size);
+    dets.emplace_back(delta_block_adaptor{delta[bl]}, p.det_init_size);
     // Set parameters
     dets.back().set_singular_threshold(p.det_singular_threshold);
     dets.back().set_n_operations_before_check(p.det_n_operations_before_check);
@@ -161,10 +162,13 @@ work_data_t::work_data_t(params_t const &p, inputs_t const &inputs, mpi::communi
 double trace_sign(work_data_t const &wdata) {
   double sign      = 1.0;
   auto const &dets = wdata.dets;
-  // For every block we compute the sign of the permutation that takes
+  // For every block, we compute the sign of the permutation that takes
+  // [(c c_dag) (c c_dag) (c c_dag) ...] with the cdag and c in increasing time order
+  // (the reference order of the det) to the  decreasing-time-and-color-ordered list
+  //  of operators (the order that makes the trace positive). 
+  // This is equivalent to computing the sign of the permutation that takes 
   // [(c_dag c) (c_dag c) (c_dag c) ...] with the cdag and c in increasing time order
-  // (the reference order of the det) to the completely time-and-color-ordered list
-  //  of operators (the order that makes the trace positive)
+  // to the increasing time-and-color-ordered list of operators. 
   for (auto bl : range(dets.size())) {
     auto s              = long(dets[bl].size());
     auto n_colors_in_bl = wdata.gf_struct[bl].second;

c++/triqs_ctseg/work_data.hpp

@@ -18,21 +18,20 @@ struct work_data_t {
   nda::vector<double> mu; // Chemical potential per color
 
   bool has_delta     = false; // There is a non-zero hybridization term
-  bool has_Dt        = false; // There is a non-zero dynamicam nn interaction
+  bool has_Dt        = false; // There is a non-zero dynamical nn interaction
   bool has_jperp     = false; // There is a non-zero Jperp interaction
-  bool rot_inv       = true;  // The spin-spin interaction is rotationally invariant
+  bool rot_inv       = true;  // The spin-spin interaction is rotationally invariant (matters for F(tau) measure)
   bool minus_sign    = false; // Has a move ever produced a negative sign?
   bool offdiag_delta = false; // Does Delta(tau) have blocks of size larger than 1?
 
   // Dynamical and spin-spin interaction kernels
   gf<imtime> K, Kprime, Jperp, Kprime_spin;
 
   // Hybridization function
-  using delta_target_t = matrix_real_valued;
-  block_gf<imtime, delta_target_t> delta;
+  block_gf<imtime, matrix_real_valued> delta;
 
   // The determinants
-  // Vector of the det_manip objects, one per block of the input Delta(tau)
+  // Vector of the det_manip objects, one per block of the input Delta(tau). See dets.hpp
   std::vector<det_t> dets;
 
   // Color to (block, idx) conversion tables

doc/ChangeLog.md

@@ -2,34 +2,6 @@
 
 # Changelog
 
-## Version 3.1.0
+## Version 3.2.0
 
-app4triqs version 3.1.0 is a compatibility
-release for TRIQS version 3.1.0 that
-* moves to cmake 3.12.4 and c++20
-* improves ghactions and jenkins configuration
-* switches documentation to read-the-docs theme
-* uses googletest main branch
-* fixes several skeleton issues
-
-We thank all contributors: Alexander Hampel, Dylan Simon, Nils Wentzell
-
-
-## Version 3.0.0
-
-app4triqs version 3.0.0 is a compatibility
-release for TRIQS version 3.0.0 that
-* introduces compatibility with Python 3 (Python 2 no longer supported)
-* adds a cmake-based dependency management
-* fixes several application issues
-
-
-## Version 2.2.0
-
-app4triqs Version 2.2.0 provides a project
-skeleton for TRIQS applications based on
-the TRIQS Library Version 2.2.0.
-It is intended for applications with both
-Python and C++ components.
-
-This is the initial release for this project.
+This is the initial release for this project, compatible with TRIQS 3.2.x.

doc/about.rst

@@ -1,8 +1,10 @@
 .. _about:
 
-About app4triqs
-***************
+About CTSEG-J
+*************
 
-An example application using ``cpp2py`` and TRIQS.
+The CTSEG-J solver has been written by N. Kavokine and O. Parcollet, with contriubtions from H. Lu and N. Wentzell. The algorithm has been 
+introduced in [J. Otsuki, Phys. Rev. B 87, 125102 (2013)]. 
 
-Written and maintained by N. Wentzell with contributions from H. U.R. Strand.
+The authors thank T. Ayral and M. Ferrero for sharing their unpublished implementation of CTSEG, which served as an inspiration 
+and a benchmark for the present implementation.