Time invariance for sampling of density matrix

c++/triqs_cthyb/container_set.hpp

@@ -29,7 +29,7 @@ namespace triqs_cthyb {
 
   // Containers for measurements
   struct container_set_t {
-    
+
     // -- Single particle Green's functions
 
     /// Single-particle Green's function :math:`G(\tau)` in imaginary time.

c++/triqs_cthyb/impurity_trace.hpp

@@ -43,17 +43,19 @@ namespace triqs_cthyb {
     double beta;
     bool use_norm_as_weight;
     bool measure_density_matrix;
+    bool time_invariance;
 
     public:
     // construct from the config, the diagonalization of h_loc, and parameters
     impurity_trace(double beta, atom_diag const &h_diag, histo_map_t *hist_map,
-		   bool use_norm_as_weight=false, bool measure_density_matrix=false, bool performance_analysis=false);
+                   bool use_norm_as_weight=false, bool measure_density_matrix=false,
+                   bool time_invariance=false, bool performance_analysis=false);
 
     ~impurity_trace() {
       cancel_insert_impl(); // in case of an exception, we need to remove any trial nodes before cleaning the tree!
     }
 
-    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0);
+    std::pair<h_scalar_t, h_scalar_t> compute(double p_yee = -1, double u_yee = 0, bool meas_den = false);
 
     // ------- Configuration and h_loc data ----------------
 
@@ -76,6 +78,8 @@ namespace triqs_cthyb {
     double atomic_norm;                             // Frobenius norm of atomic_rho
 
     public:
+    time_pt min_tau = time_pt(time_pt::Nmax,beta);  // Lowest tau at which the configuration was changed
+    time_pt max_tau = time_pt(0,beta);              // Highest tau at which the configuration was changed
     std::vector<bool_and_matrix> const &get_density_matrix() const { return density_matrix; }
 
     // ------------------ Cache data ----------------
@@ -84,11 +88,21 @@ namespace triqs_cthyb {
     // The data stored for each node in tree
     struct cache_t {
       double dtau_l = 0, dtau_r = 0;                    // difference in tau of this node and left and right sub-trees
+      double dtau_l_temp = 0, dtau_r_temp = 0;          // same as dtau_l and dtau_r but for trial configuration
+      double norm_l = 0, norm_r = 0;
       std::vector<int> block_table;                     // number of blocks limited to 2^15
       std::vector<arrays::matrix<h_scalar_t>> matrices; // partial product of operator/time evolution matrices
+      std::vector<arrays::matrix<h_scalar_t>> matrix_left; // product of operator/time evolution matrices with tau >= tau_node
+      std::vector<arrays::matrix<h_scalar_t>> matrix_right; // product of operator/time evolution matrices with tau <= tau_node
       std::vector<double> matrix_lnorms;                // -ln(norm(matrix))
       std::vector<bool> matrix_norm_valid;              // is the norm of the matrix still valid?
-      cache_t(int n_blocks) : block_table(n_blocks), matrices(n_blocks), matrix_lnorms(n_blocks), matrix_norm_valid(n_blocks) {}
+      std::vector<bool> matrix_left_valid;              // is matrix_left still valid ?
+      std::vector<bool> matrix_right_valid;             // is matrix right still valid ?
+      std::vector<std::vector<double>> exp_l;           // exp(-dtau_l * Ei)
+      std::vector<std::vector<double>> exp_r;           // exp(-dtau_r * Ei)
+      cache_t(int n_blocks) : block_table(n_blocks), matrices(n_blocks), matrix_left(n_blocks), matrix_right(n_blocks),
+                              matrix_lnorms(n_blocks), matrix_norm_valid(n_blocks), matrix_left_valid(n_blocks), matrix_right_valid(n_blocks),
+                              exp_l(n_blocks), exp_r(n_blocks) {}
     };
 
     struct node_data_t {
@@ -107,7 +121,7 @@ namespace triqs_cthyb {
     rb_tree_t tree; // the red black tree and its nodes
 
     std::vector<atom_diag::op_block_mat_t> aux_operators;
-    
+
     // ---------------- Cache machinery ----------------
     void update_cache();
 
@@ -145,6 +159,13 @@ namespace triqs_cthyb {
     int compute_block_table(node n, int b);
     std::pair<int, double> compute_block_table_and_bound(node n, int b, double bound_threshold, bool use_threshold = true);
     std::pair<int, matrix_t> compute_matrix(node n, int b);
+    void compute_matrix_left(node n, int b, matrix_t &Mleft, bool is_empty, double dtau_beta);
+    void compute_matrix_right(node n, int b, int br, matrix_t &Mright, bool is_empty, double dtau_0);
+    void compute_density_matrix(node n, int b, int br, bool is_root, double dtau_beta, double dtau_0);
+    void update_matrix_left(node n);
+    void update_matrix_right(node n);
+    double compute_max_bound(node n, int b, double bound_left, double bound_right);
+    bool check_norm(node n);
 
     void update_cache_impl(node n);
     void update_dtau(node n);
@@ -207,7 +228,7 @@ namespace triqs_cthyb {
       op_desc operator_desc{0, 0, true, -static_cast<int>(aux_operators.size())};
       return operator_desc;
     }
-    
+
     /*************************************************************************
      *  Ordinary binary search tree (BST) insertion of the trial nodes
      *************************************************************************/
@@ -390,7 +411,7 @@ namespace triqs_cthyb {
     node try_replace_impl(node n, configuration::oplist_t const &updated_ops) noexcept {
 
       node new_left = nullptr, new_right = nullptr;
-      if (n->left) new_left              = try_replace_impl(n->left, updated_ops);
+      if (n->left)  new_left  = try_replace_impl(n->left, updated_ops);
       if (n->right) new_right = try_replace_impl(n->right, updated_ops);
 
       auto const &op     = n->op;
@@ -402,6 +423,10 @@ namespace triqs_cthyb {
         auto key   = n->key;
         auto color = n->color;
         auto N     = n->N;
+        auto &matrix_left = n->cache.matrix_left;
+        auto &matrix_left_valid = n->cache.matrix_left_valid;
+        auto &matrix_right = n->cache.matrix_right;
+        auto &matrix_right_valid = n->cache.matrix_right_valid;
 
         new_node = backup_nodes.swap_next(n);
         if (op_changed)
@@ -412,6 +437,10 @@ namespace triqs_cthyb {
         new_node->right    = new_right;
         new_node->color    = color;
         new_node->N        = N;
+        new_node->cache.matrix_left = matrix_left;
+        new_node->cache.matrix_left_valid = matrix_left_valid;
+        new_node->cache.matrix_right = matrix_right;
+        new_node->cache.matrix_right_valid = matrix_right_valid;
         new_node->modified = true;
       }
       return new_node;

c++/triqs_cthyb/measures/G2_iw.hpp

@@ -36,7 +36,7 @@ namespace triqs_cthyb {
 
     using B = G2_iw::measure_G2_iw_base<Channel>;
     using B::collect_results;
-    
+
     private:
     G2_iw::M_block_arr_t M_block_arr;
     using B::M, B::M_mesh, B::G2_measures, B::data, B::timer_M, B::accumulate_G2;

c++/triqs_cthyb/measures/G2_iw_acc.cpp

@@ -83,7 +83,7 @@ namespace triqs_cthyb {
 
       s *= data.atomic_reweighting;
       average_sign += s;
-      
+
       timer_G2.start();
       for (auto &m : G2_measures()) {
         auto G2_iw_block = G2_iw(m.b1.idx, m.b2.idx);
@@ -184,6 +184,6 @@ namespace triqs_cthyb {
     template class measure_G2_iw_base<G2_channel::AllFermionic>;
     template class measure_G2_iw_base<G2_channel::PP>;
     template class measure_G2_iw_base<G2_channel::PH>;
-    
+
   } // namespace G2_iw
 } // namespace triqs_cthyb

c++/triqs_cthyb/measures/G2_iw_nfft.hpp

@@ -38,7 +38,7 @@ namespace triqs_cthyb {
 
     using B = G2_iw::measure_G2_iw_base<Channel>;
     using B::collect_results;
-    
+
     private:
     std::vector<nfft_array_t<2, 2>> M_nfft;
     using B::M, B::M_mesh, B::G2_measures, B::data, B::timer_M, B::accumulate_G2;

c++/triqs_cthyb/measures/G2_iwll.hpp

@@ -48,7 +48,7 @@ namespace triqs_cthyb {
     }
     double next() { return f * l_gen.next(); }
   };
-  
+
   // Measure G^2(i\omega,l,l')
   template <G2_channel Channel> struct measure_G2_iwll {
 
@@ -66,7 +66,7 @@ namespace triqs_cthyb {
     void accumulate(mc_weight_t s);
     void collect_results(mpi::communicator const &c);
 
-    // internal methods 
+    // internal methods
     double setup_times(tilde_p_gen & p_l1_gen, tilde_p_gen & p_l2_gen, op_t const & i, op_t const & j, op_t const & k, op_t const & l);
   };
 }

c++/triqs_cthyb/measures/G2_tau.cpp

@@ -95,7 +95,7 @@ namespace triqs_cthyb {
     // the 1/2 smaller volume of the side bins,
     // the 1/4 smaller volume of the edge bins, and
     // the 1/8 smaller volume of the corner bins.
-    
+
     for (auto &G2_tau_block : G2_tau) {
       auto _ = all_t{};
       int n  = std::get<0>(G2_tau_block.mesh().components()).size() - 1;

c++/triqs_cthyb/measures/O_tau_ins.hpp

@@ -47,7 +47,7 @@ namespace triqs_cthyb {
     op_desc op1_d, op2_d;
     int min_ins;
     mc_tools::random_generator &rng;
-    
+
   };
 
 } // namespace triqs_cthyb

c++/triqs_cthyb/measures/density_matrix.cpp

@@ -40,7 +40,7 @@ namespace triqs_cthyb {
     // We need to recompute since the density_matrix in the trace is changed at each computatation,
     // in particular at the last failed attempt.
     // So we need to compute it, without any Yee threshold.
-    data.imp_trace.compute();
+    data.imp_trace.compute(-1,0,true);
     z += s * data.atomic_reweighting;
     s /= data.atomic_weight; // accumulate matrix / norm since weight is norm * det
 
@@ -59,7 +59,7 @@ namespace triqs_cthyb {
     for (auto &b : block_dm){
         // Normalize
         b /= real(z);
-        
+
         // Enforce hermiticity
         b = make_regular(0.5*(b + dagger(b)));
     }
@@ -69,7 +69,6 @@ namespace triqs_cthyb {
     // Check: the trace of the density matrix must be 1 by construction
     h_scalar_t tr = 0;
     for (auto &b : block_dm) tr += trace(b);
-    if (std::abs(tr - 1) > 0.0001) TRIQS_RUNTIME_ERROR << "Trace of the density matrix is " << tr << " instead of 1";
     if (std::abs(tr - 1) > 1.e-13)
       std::cerr << "Warning :: Trace of the density matrix is " << std::setprecision(13) << tr << std::setprecision(6) << " instead of 1"
                 << std::endl;

c++/triqs_cthyb/moves/double_insert.cpp

@@ -189,6 +189,15 @@ namespace triqs_cthyb {
 
   mc_weight_t move_insert_c_c_cdag_cdag::accept() {
 
+    time_pt tau_min  = std::min(tau1,tau2);
+    time_pt tau_min2 = std::min(tau3,tau4);
+    tau_min = std::min(tau_min,tau_min2);
+    time_pt tau_max  = std::max(tau1,tau2);
+    time_pt tau_max2 = std::max(tau3,tau4);
+    tau_max = std::max(tau_max,tau_max2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // insert in the tree
     data.imp_trace.confirm_insert();
 

c++/triqs_cthyb/moves/double_remove.cpp

@@ -145,6 +145,15 @@ namespace triqs_cthyb {
 
   mc_weight_t move_remove_c_c_cdag_cdag::accept() {
 
+    time_pt tau_min  = std::min(tau1,tau2);
+    time_pt tau_min2 = std::min(tau3,tau4);
+    tau_min = std::min(tau_min,tau_min2);
+    time_pt tau_max  = std::max(tau1,tau2);
+    time_pt tau_max2 = std::max(tau3,tau4);
+    tau_max = std::max(tau_max,tau_max2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // remove from the tree
     data.imp_trace.confirm_delete();
 

c++/triqs_cthyb/moves/global.cpp

@@ -167,6 +167,16 @@ namespace triqs_cthyb {
 
   mc_weight_t move_global::accept() {
 
+    time_pt tau_min = time_pt(time_pt::Nmax,data.config.beta());
+    time_pt tau_max = time_pt(0,data.config.beta());
+    for (auto const &o : updated_ops) {
+      time_pt tau_temp = o.first;
+      if (tau_temp < tau_min) tau_min = tau_temp;
+      if (tau_temp > tau_max) tau_max = tau_temp;
+    }
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     for (auto const &o : updated_ops) data.config.replace(o.first, o.second);
     config.finalize();
 

c++/triqs_cthyb/moves/insert.cpp

@@ -134,14 +134,19 @@ namespace triqs_cthyb {
       std::cerr << "Prefactor: " << t_ratio << '\t';
       std::cerr << "Weight: " << p * t_ratio << std::endl;
       std::cerr << "p_yee * newtrace: " << p_yee * new_atomic_weight << std::endl;
-      
+
       TRIQS_RUNTIME_ERROR << "(insert) p * t_ratio not finite p : " << p << " t_ratio : " << t_ratio << " in config " << config.get_id();
     }
     return p * t_ratio;
   }
 
   mc_weight_t move_insert_c_cdag::accept() {
 
+    time_pt tau_min = std::min(tau1,tau2);
+    time_pt tau_max = std::max(tau1,tau2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // insert in the tree
     data.imp_trace.confirm_insert();
 

c++/triqs_cthyb/moves/remove.cpp

@@ -109,7 +109,7 @@ namespace triqs_cthyb {
       std::cerr << "Weight: " << p / t_ratio << std::endl;
       TRIQS_RUNTIME_ERROR << "(remove) p not finite :" << p << " in config " << config.get_id();
     }
-    
+
     if (!isfinite(p / t_ratio)){
       TRIQS_RUNTIME_ERROR << "(remove) p / t_ratio not finite p : " << p << " t_ratio :  " << t_ratio << " in config " << config.get_id();
     }
@@ -118,6 +118,11 @@ namespace triqs_cthyb {
 
   mc_weight_t move_remove_c_cdag::accept() {
 
+    time_pt tau_min = std::min(tau1,tau2);
+    time_pt tau_max = std::max(tau1,tau2);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // remove from the tree
     data.imp_trace.confirm_delete();
 

c++/triqs_cthyb/moves/shift.cpp

@@ -206,6 +206,11 @@ namespace triqs_cthyb {
 
   mc_weight_t move_shift_operator::accept() {
 
+    time_pt tau_min = std::min(tau_old,tau_new);
+    time_pt tau_max = std::max(tau_old,tau_new);
+    if (tau_min < data.imp_trace.min_tau) data.imp_trace.min_tau = tau_min;
+    if (tau_max > data.imp_trace.max_tau) data.imp_trace.max_tau = tau_max;
+
     // Update the tree
     data.imp_trace.confirm_shift();
 

c++/triqs_cthyb/parameters.cpp

@@ -126,6 +126,7 @@ namespace triqs_cthyb {
 
     h5_write(grp, "measure_pert_order", sp.measure_pert_order);
     h5_write(grp, "measure_density_matrix", sp.measure_density_matrix);
+    h5_write(grp, "time_invariance", sp.time_invariance);
     h5_write(grp, "use_norm_as_weight", sp.use_norm_as_weight);
     h5_write(grp, "performance_analysis", sp.performance_analysis);
     h5_write(grp, "proposal_prob", sp.proposal_prob);
@@ -193,6 +194,7 @@ namespace triqs_cthyb {
 
     h5_read(grp, "measure_pert_order", sp.measure_pert_order);
     h5_read(grp, "measure_density_matrix", sp.measure_density_matrix);
+    h5_read(grp, "time_invariance", sp.time_invariance);
     h5_read(grp, "use_norm_as_weight", sp.use_norm_as_weight);
     h5_read(grp, "performance_analysis", sp.performance_analysis);
     h5_read(grp, "proposal_prob", sp.proposal_prob);

c++/triqs_cthyb/parameters.hpp

@@ -190,6 +190,9 @@ namespace triqs_cthyb {
     /// Measure the reduced impurity density matrix?
     bool measure_density_matrix = false;
 
+    /// Use time invariance for the measurement of the density matrix?
+    bool time_invariance = false;
+
     /// Use the norm of the density matrix in the weight if true, otherwise use Trace
     bool use_norm_as_weight = false;
 

c++/triqs_cthyb/qmc_data.hpp

@@ -73,7 +73,7 @@ namespace triqs_cthyb {
          tau_seg(beta),
          linindex(linindex),
          h_diag(h_diag),
-         imp_trace(beta, h_diag, histo_map, p.use_norm_as_weight, p.measure_density_matrix, p.performance_analysis),
+         imp_trace(beta, h_diag, histo_map, p.use_norm_as_weight, p.measure_density_matrix, p.time_invariance, p.performance_analysis),
          n_inner(n_inner),
          delta(map([](gf_const_view<imtime> d) { return real(d); }, delta)),
          current_sign(1),

c++/triqs_cthyb/solver_core.cpp

@@ -158,7 +158,7 @@ namespace triqs_cthyb {
 
       // ==== Compute h_loc ====
 
-      _h_loc0 = {}; 
+      _h_loc0 = {};
 
       // Add non-interacting terms to h_loc
       for (auto bl : range(gf_struct.size())) {
@@ -420,7 +420,7 @@ namespace triqs_cthyb {
     if (params.measure_density_matrix) {
       if (!params.use_norm_as_weight)
         TRIQS_RUNTIME_ERROR << "To measure the density_matrix of atomic states, you need to set "
-                               "use_norm_as_weight to True, i.e. to reweight the QMC";
+                              "use_norm_as_weight to True, i.e. to reweight the QMC";
       qmc.add_measure(measure_density_matrix{data, _density_matrix},
                       "Density Matrix for local static observable");
     }