Add possibility of starting the run from a user provided configuration

c++/triqs_cthyb/configuration.hpp

@@ -51,11 +51,23 @@ namespace triqs_cthyb {
       h5_write(g, "inner", op.inner_index);
       h5_write(g, "dagger", op.dagger);
     }
+
+    friend void h5_read(h5::group g, op_desc &op) {
+      h5_read(g, "block", op.block_index);
+      h5_read(g, "inner", op.inner_index);
+      h5_read(g, "dagger", op.dagger);
+    }
+
+    bool operator==(op_desc const &op) const { return (block_index == op.block_index && inner_index == op.inner_index
+      && dagger == op.dagger && linear_index == op.linear_index);
+    }
   };
 
   // The configuration of the Monte Carlo
   struct configuration {
 
+    bool operator==(configuration const &config) const { return (beta_ == config.beta_ && oplist == config.oplist); }
+
     // a map associating an operator to an imaginary time
     using oplist_t = std::map<time_pt, op_desc, std::greater<time_pt>>;
 
@@ -66,6 +78,7 @@ namespace triqs_cthyb {
     ~configuration() { configs_hfile.close(); }
 #else
     configuration(double beta) : beta_(beta), id(0) {}
+    configuration() : beta_(double(0.)) {}
 #endif
 
     double beta() const { return beta_; }
@@ -86,15 +99,36 @@ namespace triqs_cthyb {
       return out;
     }
 
+    static std::string hdf5_format() { return "Configuration"; }
+
     // Writing of configuration out to a h5 for e.g. plotting
     friend void h5_write(h5::group conf, std::string conf_group_name, configuration const &c) {
+      auto beta = c.beta();
+      auto t1 = time_pt(1,beta);
       h5::group conf_group = conf.create_group(conf_group_name);
+      write_hdf5_format(conf_group, c);
+      h5_write(conf_group, "beta", beta);
       for (auto const &op : c) {
         // create group for given tau
         auto tau_group_name        = std::to_string(double(op.first));
         h5::group tau_group        = conf_group.create_group(tau_group_name);
         // in tau subgroup, write operator info
         h5_write(tau_group, op.second);
+        h5_write(tau_group, "n", floor_div(op.first, t1));
+      }
+    }
+
+    friend void h5_read(h5::group conf, std::string conf_group_name, configuration &c) {
+      h5::group conf_group = conf.open_group(conf_group_name);
+      op_desc op;
+      uint64_t n;
+      double beta;
+      h5_read(conf_group, "beta", beta);
+      for (auto &sgrp : conf_group.get_all_subgroup_names()) {
+        auto tau_group = conf_group.open_group(sgrp);
+        h5_read(tau_group, op);
+        h5_read(tau_group, "n", n);
+        c.insert(time_pt(n,beta), op);
       }
     }
 

c++/triqs_cthyb/parameters.hpp

@@ -25,6 +25,7 @@
 
 #include "./config.hpp"
 #include "./types.hpp"
+#include "./configuration.hpp"
 
 namespace triqs_cthyb {
 
@@ -193,6 +194,9 @@ namespace triqs_cthyb {
     /// Use the norm of the density matrix in the weight if true, otherwise use Trace
     bool use_norm_as_weight = false;
 
+    /// Initial configuration of the run
+    configuration initial_configuration;
+
     /// Analyse performance of trace computation with histograms (developers only)?
     bool performance_analysis = false;
 

c++/triqs_cthyb/qmc_data.hpp

@@ -34,7 +34,7 @@ namespace triqs_cthyb {
  ***********************/
   struct qmc_data {
 
-    configuration config; // Configuration
+    configuration &config; // Configuration
     time_segment tau_seg;
     std::map<std::pair<int, int>, int> linindex; // Linear index constructed from block and inner indices
     atom_diag const &h_diag;                     // Diagonalization of the atomic problem
@@ -68,8 +68,8 @@ namespace triqs_cthyb {
 
     // Construction
     qmc_data(double beta, solve_parameters_t const &p, atom_diag const &h_diag, std::map<std::pair<int, int>, int> linindex,
-             block_gf_const_view<imtime> delta, std::vector<int> n_inner, histo_map_t *histo_map)
-       : config(beta),
+             block_gf_const_view<imtime> delta, std::vector<int> n_inner, histo_map_t *histo_map, configuration &c)
+       : config(c),
          tau_seg(beta),
          linindex(linindex),
          h_diag(h_diag),
@@ -78,11 +78,36 @@ namespace triqs_cthyb {
          delta(map([](gf_const_view<imtime> d) { return real(d); }, delta)),
          current_sign(1),
          old_sign(1) {
+      config.clear();
+      std::vector<std::vector<std::pair<time_pt, int>>> X(delta.size()), Y(delta.size());
+      for (auto const &o : p.initial_configuration) {
+        auto tau = o.first;
+        auto op  = o.second;
+        auto block_index = op.block_index;
+        auto inner_index = op.inner_index;
+        if (block_index >= delta.size() || inner_index >= n_inner[block_index])
+          TRIQS_RUNTIME_ERROR << "Inconsistency in the block structure of the initial configuration";
+
+        op.linear_index = linindex[std::make_pair(block_index, inner_index)];
+
+        imp_trace.try_insert(tau, op);
+        imp_trace.confirm_insert();
+
+        if (op.dagger) X[block_index].push_back(std::make_pair(tau,inner_index));
+        else Y[block_index].push_back(std::make_pair(tau,inner_index));
+
+        config.insert(tau, op);
+        config.finalize();
+      }
       std::tie(atomic_weight, atomic_reweighting) = imp_trace.compute();
+
       dets.clear();
+      dets.reserve(delta.size());
       for (auto const &bl : range(delta.size())) {
+        if (X[bl].size() != Y[bl].size())
+          TRIQS_RUNTIME_ERROR << "In the initial config, the number of c and c_dag operators should be equal";
 #ifdef HYBRIDISATION_IS_COMPLEX
-        dets.emplace_back(delta_block_adaptor(delta[bl]), p.det_init_size);
+        dets.emplace_back(delta_block_adaptor(delta[bl]), X[bl], Y[bl]);
 #else
         if (!is_gf_real(delta[bl], 1e-10)) {
           //TRIQS_RUNTIME_ERROR << "The Delta(tau) block number " << bl << " is not real in tau space";
@@ -92,13 +117,14 @@ namespace triqs_cthyb {
             std::cerr << "WARNING: Dissregarding the imaginary component in the calculation.\n";
           }
         }
-        dets.emplace_back(delta_block_adaptor(real(delta[bl])), p.det_init_size);
+        dets.emplace_back(delta_block_adaptor(real(delta[bl])), X[bl], Y[bl]);
 #endif
         dets.back().set_singular_threshold(p.det_singular_threshold);
         dets.back().set_n_operations_before_check(p.det_n_operations_before_check);
         dets.back().set_precision_warning(p.det_precision_warning);
         dets.back().set_precision_error(p.det_precision_error);
       }
+      update_sign();
     }
 
     qmc_data(qmc_data const &) = delete; // Member imp_trace is not copyable

c++/triqs_cthyb/solver_core.cpp

@@ -61,7 +61,7 @@ namespace triqs_cthyb {
   };
 
   solver_core::solver_core(constr_parameters_t const &p)
-     : beta(p.beta), gf_struct(p.gf_struct), n_iw(p.n_iw), n_tau(p.n_tau), n_l(p.n_l), delta_interface(p.delta_interface), constr_parameters(p) {
+     : beta(p.beta), gf_struct(p.gf_struct), n_iw(p.n_iw), n_tau(p.n_tau), n_l(p.n_l), delta_interface(p.delta_interface), _configuration(p.beta), constr_parameters(p) {
 
     if (p.n_tau < 2 * p.n_iw)
       TRIQS_RUNTIME_ERROR
@@ -262,7 +262,7 @@ namespace triqs_cthyb {
     }
 
     // Initialise Monte Carlo quantities
-    qmc_data data(beta, params, h_diag, linindex, _Delta_tau, n_inner, histo_map);
+    qmc_data data(beta, params, h_diag, linindex, _Delta_tau, n_inner, histo_map, _configuration);
     auto qmc =
        mc_tools::mc_generic<mc_weight_t>(params.random_name, params.random_seed, params.verbosity);
 
@@ -431,10 +431,17 @@ namespace triqs_cthyb {
 
     // --------------------------------------------------------------------------
 
+    mc_weight_t sign = data.current_sign * data.atomic_weight / std::abs(data.atomic_weight);
+
+    for (size_t block = 0; block < _Delta_tau.size(); ++block) {
+      auto det = data.dets[block].determinant();
+      sign *= det / std::abs(det);
+    }
+
     // Run! The empty (starting) configuration has sign = 1
     _solve_status =
        qmc.warmup_and_accumulate(params.n_warmup_cycles, params.n_cycles, params.length_cycle,
-                                 triqs::utility::clock_callback(params.max_time));
+                                 triqs::utility::clock_callback(params.max_time), sign);
     qmc.collect_results(_comm);
 
     if (params.verbosity >= 2) {

c++/triqs_cthyb/solver_core.hpp

@@ -32,6 +32,7 @@
 #include "types.hpp"
 #include "container_set.hpp"
 #include "parameters.hpp"
+#include "configuration.hpp"
 
 namespace triqs_cthyb {
 
@@ -53,6 +54,7 @@ namespace triqs_cthyb {
     double _average_order;                 // average perturbation order
     double _auto_corr_time;                // Auto-correlation time
     int _solve_status;                     // Status of the solve upon exit: 0 for clean termination, > 0 otherwise.
+    configuration _configuration;          // Final configuration of the run
 
     // Single-particle Green's function containers
     std::optional<G_iw_t> _G0_iw; // Non-interacting Matsubara Green's function
@@ -153,6 +155,9 @@ namespace triqs_cthyb {
     /// Status of the ``solve()`` on exit.
     int solve_status() const { return _solve_status; }
 
+    /// Configuration
+    configuration const &get_configuration() const { return _configuration; }
+
     /// is cthyb compiled with support for complex hybridization?
     bool hybridisation_is_complex() const {
 #ifdef HYBRIDISATION_IS_COMPLEX

python/triqs_cthyb/solver_core_desc.py

@@ -227,6 +227,8 @@
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | use_norm_as_weight            | bool                                                     | false                         | Use the norm of the density matrix in the weight if true, otherwise use Trace                                     |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
+| initial_configuration         | Configuration                                            | {}                            | Initial configuration for the run                                                                                  |
++-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | performance_analysis          | bool                                                     | false                         | Analyse performance of trace computation with histograms (developers only)?                                       |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | proposal_prob                 | dict(str:float)                                          | {}                            | Operator insertion/removal probabilities for different blocks                                                     |
@@ -313,6 +315,10 @@
                getter = cfunction("int solve_status ()"),
                doc = r"""status of the ``solve()`` on exit.""")
 
+c.add_property(name = "configuration",
+               getter = cfunction("configuration get_configuration()"),
+               doc = r"""Configuration""")
+
 c.add_property(name = "hybridisation_is_complex",
                getter = cfunction("bool hybridisation_is_complex ()"),
                doc = r"""cthyb compiled with support for complex hybridization?""")
@@ -537,6 +543,11 @@
              initializer = """ false """,
              doc = r"""Use the norm of the density matrix in the weight if true, otherwise use Trace""")
 
+c.add_member(c_name = "initial_configuration",
+             c_type = "triqs_cthyb::configuration",
+             initializer = """ {} """,
+             doc = r"""Initial configuration for the run""")
+
 c.add_member(c_name = "performance_analysis",
              c_type = "bool",
              initializer = """ false """,
@@ -641,5 +652,30 @@
 
 module.add_converter(c)
 
+# The class configuration
+c = class_(
+        py_type = "Configuration",  # name of the python class
+        c_type = "triqs_cthyb::configuration",   # name of the C++ class
+        doc = r"""Core class of the cthyb configuration""",   # doc of the C++ class
+        comparisons = "==",
+        is_printable = True,
+        hdf5 = True
+)
+
+c.add_constructor(signature = "()",
+                  doc = "Create empty configuration")
+
+c.add_constructor(signature = "(double beta)",
+                  doc = "Initialize beta")
+
+c.add_property(name = "beta",
+               getter = cfunction("double beta()"),
+               doc = r"""Value of beta""")
+
+c.add_property(name = "size",
+               getter = cfunction("int size()"),
+               doc = r"""Size of the configuration""")
+
+module.add_class(c)
 
 module.generate_code()