New features for internal energy calculations

c++/triqs_cthyb/configuration.hpp

@@ -70,6 +70,7 @@ namespace triqs_cthyb {
 
     double beta() const { return beta_; }
     int size() const { return oplist.size(); }
+    oplist_t const &get_oplist() const { return oplist; }  
 
     void insert(time_pt tau, op_desc op) { oplist.insert({tau, op}); }
     void replace(time_pt tau, op_desc op) { oplist[tau] = op; }
@@ -97,7 +98,7 @@ namespace triqs_cthyb {
         h5_write(tau_group, op.second);
       }
     }
-
+    
     long get_id() const { return id; } // Get the id of the current configuration
     void finalize() {
       id++;

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);  
+    time_pt max_tau = time_pt(0,beta);  
     std::vector<bool_and_matrix> const &get_density_matrix() const { return density_matrix; }
 
     // ------------------ Cache data ----------------
@@ -84,11 +88,20 @@ 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;  
       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;   
+      std::vector<arrays::matrix<h_scalar_t>> matrix_right;  
       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;   
+      std::vector<bool> matrix_right_valid;  
+      std::vector<std::vector<double>> exp_l;  
+      std::vector<std::vector<double>> exp_r;  
+      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 {
@@ -145,6 +158,11 @@ 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); 
 
     void update_cache_impl(node n);
     void update_dtau(node n);
@@ -402,6 +420,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 +434,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/G_tau.cpp

@@ -31,7 +31,6 @@ namespace triqs_cthyb {
     results.G_tau_accum = block_gf<imtime, G_target_t>({data.config.beta(), Fermion, n_tau}, gf_struct);
     G_tau.rebind(*results.G_tau_accum);
     G_tau() = 0.0;
-
     results.asymmetry_G_tau = block_gf{G_tau};
     asymmetry_G_tau.rebind(*results.asymmetry_G_tau);
   }
@@ -47,7 +46,7 @@ namespace triqs_cthyb {
         double dtau = double(y.first - x.first);
         this->G_tau[block_idx][closest_mesh_pt(dtau)](y.second, x.second) += val;
       })
-        ;
+        ; 
     }
   }
 

c++/triqs_cthyb/measures/density_matrix.cpp

@@ -40,20 +40,33 @@ 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();
-    z += s * data.atomic_reweighting;
-    s /= data.atomic_weight; // accumulate matrix / norm since weight is norm * det
 
     // Careful: there is no reweighting factor here!
-    int size = block_dm.size();
-    for (int i = 0; i < size; ++i)
-      if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s * data.imp_trace.get_density_matrix()[i].mat; }
+    if (data.updated) {
+      if (!flag) {
+        z += data.n_acc * old_z;
+	mc_weight_t s_temp = data.n_acc * old_s;
+	int size = block_dm.size();
+	for (int i = 0; i < size; ++i) 
+	  if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s_temp * data.imp_trace.get_density_matrix()[i].mat; }
+      }	
+      data.imp_trace.compute(-1,0,true);
+      old_z = s * data.atomic_reweighting;
+      old_s = s / data.atomic_weight;
+    }
+    flag = false;
   }
 
   // ---------------------------------------------
 
   void measure_density_matrix::collect_results(mpi::communicator const &c) {
 
+    z += data.n_acc * old_z;
+    mc_weight_t s_temp = data.n_acc * old_s;
+    int size = block_dm.size();
+    for (int i = 0; i < size; ++i)
+      if (data.imp_trace.get_density_matrix()[i].is_valid) { block_dm[i] += s_temp * data.imp_trace.get_density_matrix()[i].mat; }
+
     z                          = mpi::all_reduce(z, c);
     block_dm                   = mpi::all_reduce(block_dm, c);
     for (auto &b : block_dm){
@@ -69,7 +82,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/measures/density_matrix.hpp

@@ -27,6 +27,9 @@ namespace triqs_cthyb {
     qmc_data const &data;
     std::vector<matrix_t> &block_dm; // density matrix of each block
     mc_weight_t z = 0;
+    bool flag = true;  
+    mc_weight_t old_z = 0; 
+    mc_weight_t old_s = 0; 
 
     measure_density_matrix(qmc_data const &data, std::vector<matrix_t> &density_matrix);
     void accumulate(mc_weight_t s);

c++/triqs_cthyb/measures/update_time.hpp

@@ -0,0 +1,72 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2021, Simons Foundation
+ *    author: N. Wentzell
+ *
+ * TRIQS is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU General Public License as published by the Free Software
+ * Foundation, either version 3 of the License, or (at your option) any later
+ * version.
+ *
+ * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU General Public License along with
+ * TRIQS. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "../qmc_data.hpp"
+
+namespace triqs_cthyb {
+
+  /// Measure of the average update time
+  struct measure_update_time {
+
+    measure_update_time(qmc_data &_data, double &_update_time) : data(_data), update_time(_update_time) { update_time = 0.0; }
+
+    void accumulate(mc_weight_t) {
+      if (data.updated && !flag) {
+        update_time += data.n_acc;
+        data.n_acc = 1.;
+	++N;
+	data.updated = false;
+      }
+      else
+        data.n_acc += 1.;
+      if (flag) data.updated = false; // need to set it back to false at the first measurement since it has been set to true during warmup
+      flag = false;
+    }
+
+    void collect_results(mpi::communicator const &comm) {
+      if (N == 0) {
+        N = 1;
+	update_time = data.n_acc;
+      }
+
+      N = mpi::all_reduce(N, comm);
+
+      // Reduce and normalize
+      update_time = mpi::all_reduce(update_time, comm);
+      update_time = update_time / N;
+    }
+
+    private:
+    // The Monte-Carlo configuration
+    qmc_data &data;
+
+    // Flag to indicate whether or not this is the first measure
+    bool flag = true;
+
+    // Reference to double for accumulation
+    double &update_time;
+
+    // Accumulation counter
+    long long N = 0;
+  };
+
+} // namespace triqs_cthyb

c++/triqs_cthyb/moves/double_insert.cpp

@@ -189,6 +189,16 @@ namespace triqs_cthyb {
 
   mc_weight_t move_insert_c_c_cdag_cdag::accept() {
 
+    data.updated = true;
+    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,16 @@ namespace triqs_cthyb {
 
   mc_weight_t move_remove_c_c_cdag_cdag::accept() {
 
+    data.updated = true;
+    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,17 @@ namespace triqs_cthyb {
 
   mc_weight_t move_global::accept() {
 
+    data.updated = true;
+    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();