Use of time invariance for sampling of the density matrix

Add files via upload

Add files via upload

Add files via upload

Add files via upload

Add files via upload

Add files via upload

Delete solver_core.hpp

Delete configuration.hpp

Delete solver_core.cpp

Delete qmc_data.hpp

Delete parameters.hpp

Delete impurity_trace.hpp

Add files via upload

Update CMakeLists.txt

fixes.cpp

Update qmc_data.hpp

Update solver_core.hpp

Update configuration.hpp

Update solver_core.cpp

Update parameters.hpp

Update solver_core.cpp

Update insert.hpp

Update insert.cpp

Update remove.hpp

Update remove.cpp

Update insert.hpp

Update insert.cpp

Update remove.hpp

Update remove.cpp

Update insert.cpp

Update remove.cpp

Update insert.cpp

Update remove.cpp

Update solver_core.cpp

Update update_time.hpp

Add files via upload

Delete c++/triqs_cthyb/measures/update_time.hpp

Add files via upload

Update parameters.cpp

Update parameters.hpp

Update solver_core_desc.py

Update impurity_trace.cpp

Update impurity_trace.cpp

Update impurity_trace.cpp

Update impurity_trace.cpp

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,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;          // same as dtau_l and dtau_r but for trial configuration
       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 {
@@ -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/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
 
@@ -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

@@ -188,6 +188,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

@@ -144,6 +144,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

@@ -166,6 +166,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

@@ -141,6 +141,11 @@ namespace triqs_cthyb {
   }
 
   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

@@ -117,6 +117,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

@@ -205,6 +205,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

@@ -189,6 +189,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

@@ -418,9 +418,9 @@ namespace triqs_cthyb {
       qmc.add_measure(measure_perturbation_hist_total(data, *perturbation_order_total), "Perturbation order");
     }
     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";
+      //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";
       qmc.add_measure(measure_density_matrix{data, _density_matrix},
                       "Density Matrix for local static observable");
     }

python/triqs_cthyb/solver_core_desc.py

@@ -225,6 +225,8 @@
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | measure_density_matrix        | bool                                                     | false                         | Measure the reduced impurity density matrix? Automatically also determines high frequency moments for G and Sigma |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
+| time_invariance               | bool                                                     | false                         | Use time invariance to sample the density matrix?                                                                 |
++-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | use_norm_as_weight            | bool                                                     | false                         | Use the norm of the density matrix in the weight if true, otherwise use Trace                                     |
 +-------------------------------+----------------------------------------------------------+-------------------------------+-------------------------------------------------------------------------------------------------------------------+
 | performance_analysis          | bool                                                     | false                         | Analyse performance of trace computation with histograms (developers only)?                                       |
@@ -531,6 +533,11 @@
              c_type = "bool",
              initializer = """ false """,
              doc = r"""Measure the reduced impurity density matrix?""")
+
+c.add_member(c_name = "time_invariance",
+             c_type = "bool",
+             initializer = """ false """,
+             doc = r"""Use time invariance to measure the density matrix?""")
 
 c.add_member(c_name = "use_norm_as_weight",
              c_type = "bool",