From d367dd1cff7a442a2480f0dd76f93f506dcf4e3c Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:05:22 +0200 Subject: [PATCH] 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.cpp | 324 ++++++++++++++++++-- c++/triqs_cthyb/impurity_trace.hpp | 32 +- c++/triqs_cthyb/measures/density_matrix.cpp | 3 +- c++/triqs_cthyb/moves/double_insert.cpp | 9 + c++/triqs_cthyb/moves/double_remove.cpp | 9 + c++/triqs_cthyb/moves/global.cpp | 10 + c++/triqs_cthyb/moves/insert.cpp | 5 + c++/triqs_cthyb/moves/remove.cpp | 5 + c++/triqs_cthyb/moves/shift.cpp | 5 + c++/triqs_cthyb/parameters.cpp | 2 + c++/triqs_cthyb/parameters.hpp | 3 + c++/triqs_cthyb/qmc_data.hpp | 2 +- c++/triqs_cthyb/solver_core.cpp | 6 +- python/triqs_cthyb/solver_core_desc.py | 7 + 14 files changed, 383 insertions(+), 39 deletions(-) diff --git a/c++/triqs_cthyb/impurity_trace.cpp b/c++/triqs_cthyb/impurity_trace.cpp index ec84b14c..e3811145 100644 --- a/c++/triqs_cthyb/impurity_trace.cpp +++ b/c++/triqs_cthyb/impurity_trace.cpp @@ -50,22 +50,23 @@ namespace triqs_cthyb { // -------- Constructor -------- impurity_trace::impurity_trace(double beta, atom_diag const &h_diag_, histo_map_t *hist_map, bool use_norm_as_weight, bool measure_density_matrix, - bool performance_analysis) + bool time_invariance, bool performance_analysis) : beta(beta), use_norm_as_weight(use_norm_as_weight), measure_density_matrix(measure_density_matrix), + time_invariance(time_invariance), h_diag(&h_diag_), density_matrix(n_blocks), atomic_rho(n_blocks), atomic_z(partition_function(*h_diag, beta)), atomic_norm(0), - histo(performance_analysis ? new histograms_t(h_diag_.n_subspaces(), *hist_map) : nullptr) { + histo(performance_analysis ? new histograms_t(h_diag_.n_subspaces(), *hist_map) : nullptr) { // init density_matrix block + bool for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl] = bool_and_matrix{false, matrix_t(get_block_dim(bl), get_block_dim(bl))}; // prepare atomic_rho and atomic_norm - if (use_norm_as_weight) { + if (use_norm_as_weight || measure_density_matrix) { auto rho = atomic_density_matrix(h_diag_, beta); for (int bl = 0; bl < n_blocks; ++bl) { atomic_rho[bl] = bool_and_matrix{true, rho[bl] * atomic_z}; @@ -127,7 +128,7 @@ namespace triqs_cthyb { if (n->right) { std::tie(b1, lnorm) = compute_block_table_and_bound(n->right, b, lnorm_threshold, use_threshold); if (b1 < 0) return {b1, 0}; - lnorm += n->cache.dtau_r * get_block_emin(b1); + lnorm += n->cache.dtau_r_temp * get_block_emin(b1); } if (use_threshold && (lnorm > lnorm_threshold)) return {-1, 0}; @@ -135,8 +136,9 @@ namespace triqs_cthyb { if (b2 < 0) return {b2, 0}; int b3 = b2; + if (n->left) { - lnorm += n->cache.dtau_l * get_block_emin(b2); + lnorm += n->cache.dtau_l_temp * get_block_emin(b2); if (use_threshold && (lnorm > lnorm_threshold)) return {-1, 0}; double lnorm3; std::tie(b3, lnorm3) = compute_block_table_and_bound(n->left, b2, lnorm_threshold, use_threshold); @@ -177,9 +179,17 @@ namespace triqs_cthyb { matrix_t M = (!n->delete_flag ? get_op_block_matrix(n, b1) : nda::eye(get_block_dim(b1))); if (n->right) { // M <- M * exp * r[b] - dtau_r = double(n->key - tree.min_key(n->right)); + if (n->modified) dtau_r = n->cache.dtau_r_temp; + else dtau_r = n->cache.dtau_r; auto dim = M.shape()[1]; // same as get_block_dim(b2); - for (int i = 0; i < dim; ++i) M(_, i) *= std::exp(-dtau_r * get_block_eigenval(b1, i)); // Create time-evolution matrix e^-H(t'-t) + if (updating) { + if (n->cache.exp_r[b1].empty()) n->cache.exp_r[b1].resize(dim); + for (int i = 0; i < dim; ++i) { + n->cache.exp_r[b1][i] = std::exp(-dtau_r * get_block_eigenval(b1, i)); + M(_, i) *= n->cache.exp_r[b1][i]; + } + } + else for (int i = 0; i < dim; ++i) M(_, i) *= std::exp(-dtau_r * get_block_eigenval(b1, i)); if ((r.second.shape()[0] == 1) && (r.second.shape()[1] == 1)) M *= r.second(0, 0); else @@ -191,9 +201,17 @@ namespace triqs_cthyb { auto l = compute_matrix(n->left, b2); b3 = l.first; if (b3 == -1) return {-1, {}}; - dtau_l = double(tree.max_key(n->left) - n->key); + if (n->modified) dtau_l = n->cache.dtau_l_temp; + else dtau_l = n->cache.dtau_l; auto dim = M.shape()[0]; // same as get_block_dim(b1); - for (int i = 0; i < dim; ++i) M(i, _) *= std::exp(-dtau_l * get_block_eigenval(b2, i)); + if (updating) { + if (n->cache.exp_l[b2].empty()) n->cache.exp_l[b2].resize(dim); + for (int i = 0; i < dim; ++i) { + n->cache.exp_l[b2][i] = std::exp(-dtau_l * get_block_eigenval(b2, i)); + M(i, _) *= n->cache.exp_l[b2][i]; + } + } + else for (int i = 0; i < dim; ++i) M(i, _) *= std::exp(-dtau_l * get_block_eigenval(b2, i)); if ((l.second.shape()[0] == 1) && (l.second.shape()[1] == 1)) M *= l.second(0, 0); else @@ -218,9 +236,221 @@ namespace triqs_cthyb { return {b3, std::move(M)}; } + // At each node, computes recursively the full left matrix (keeping only the operators with tau >= tau_node) for block b + + void impurity_trace::compute_matrix_left(node n, int b, matrix_t &Mleft, bool is_empty, double dtau_beta) { + + auto _ = arrays::range(); + + int b1 = b; + if (n->right) b1 = (n->right)->cache.block_table[b]; + int b2 = get_op_block_map(n,b1); + + int dimb1 = get_block_dim(b1); + int dimb2 = get_block_dim(b2); + + if (!n->left && is_empty && !n->cache.matrix_left_valid[b1]) { + if (n->cache.matrix_left[b1].shape()[0] != dimb2 || n->cache.matrix_left[b1].shape()[1] != dimb1) + n->cache.matrix_left[b1] = matrix_t(dimb2,dimb1); + for (int i = 0; i < dimb2; ++i) + n->cache.matrix_left[b1](i,_) = get_op_block_matrix(n,b1)(i,_) * std::exp(- dtau_beta * get_block_eigenval(b2,i)) ; + n->cache.matrix_left_valid[b1] = true; + } + + if (!n->left && !is_empty && !n->cache.matrix_left_valid[b1]) { + if (dimb1==1 && dimb2==1) + n->cache.matrix_left[b1] = Mleft * get_op_block_matrix(n,b1)(0,0); + else + n->cache.matrix_left[b1] = Mleft * get_op_block_matrix(n,b1); + n->cache.matrix_left_valid[b1] = true; + } + + if (n->left && !n->cache.matrix_left_valid[b1]) { + compute_matrix_left(n->left, b2, Mleft, is_empty, dtau_beta); + node x = tree.max(n->left); + if (n->cache.matrix_left[b1].shape()[0] != dimb2 || n->cache.matrix_left[b1].shape()[1] != dimb1) + n->cache.matrix_left[b1] = matrix_t(dimb2,dimb1); + for (int i = 0; i < dimb2; ++i) + n->cache.matrix_left[b1](i,_) = get_op_block_matrix(n,b1)(i,_) * n->cache.exp_l[b2][i]; + if (dimb1==1 && dimb2==1) + n->cache.matrix_left[b1] = x->cache.matrix_left[b2] * n->cache.matrix_left[b1](0,0); + else + n->cache.matrix_left[b1] = x->cache.matrix_left[b2] * n->cache.matrix_left[b1]; + n->cache.matrix_left_valid[b1] = true; + } + + if (n->right) { + int dim = n->cache.matrix_left[b1].shape()[0]; + if (Mleft.shape()[0] != dim || Mleft.shape()[1] != dimb1) Mleft = matrix_t(dim,dimb1); + for (int i = 0; i < dimb1; ++i) + Mleft(_,i) = n->cache.matrix_left[b1](_,i) * n->cache.exp_r[b1][i]; + compute_matrix_left(n->right, b, Mleft, false, dtau_beta); + } + } + + // At each node, computes recursively the full right matrix (keeping only the operators with tau <= tau_node) for block b + + void impurity_trace::compute_matrix_right(node n, int b, int br, matrix_t &Mright, bool is_empty, double dtau_0) { + + auto _ = arrays::range(); + + int b1 = br ; + if (n->right) b1 = (n->right)->cache.block_table[br]; + int b2 = get_op_block_map(n,b1); + + int dimb1 = get_block_dim(b1); + int dimb2 = get_block_dim(b2); + int dimb = get_block_dim(b); + + if (!n->right && is_empty && !n->cache.matrix_right_valid[b]) { + if (n->cache.matrix_right[b].shape()[0] != dimb2 || n->cache.matrix_right[b].shape()[1] != dimb) + n->cache.matrix_right[b] = matrix_t(dimb2,dimb); + for (int i = 0; i < dimb; ++i) + n->cache.matrix_right[b](_,i) = get_op_block_matrix(n,b)(_,i) * std::exp(- dtau_0 * get_block_eigenval(b,i)) ; + n->cache.matrix_right_valid[b] = true; + } + + if (!n->right && !is_empty && !n->cache.matrix_right_valid[b]) { + if (dimb1==1 && dimb2==1) + n->cache.matrix_right[b] = get_op_block_matrix(n,b1)(0,0) * Mright; + else + n->cache.matrix_right[b] = get_op_block_matrix(n,b1) * Mright; + n->cache.matrix_right_valid[b] = true; + } + + if (n->right && !n->cache.matrix_right_valid[b]) { + compute_matrix_right(n->right, b, br, Mright, is_empty, dtau_0); + node x = tree.min(n->right); + if (n->cache.matrix_right[b].shape()[0] != dimb2 || n->cache.matrix_right[b].shape()[1] != dimb1) + n->cache.matrix_right[b] = matrix_t(dimb2,dimb1); + for (int i = 0; i < dimb1; ++i) + n->cache.matrix_right[b](_,i) = get_op_block_matrix(n,b1)(_,i) * n->cache.exp_r[b1][i]; + if (dimb1==1 && dimb2==1) + n->cache.matrix_right[b] = n->cache.matrix_right[b](0,0) * x->cache.matrix_right[b]; + else + n->cache.matrix_right[b] = n->cache.matrix_right[b] * x->cache.matrix_right[b]; + n->cache.matrix_right_valid[b] = true; + } + + if (n->left) { + if (Mright.shape()[0] != dimb2 || Mright.shape()[1] != dimb) Mright = matrix_t(dimb2,dimb); + for (int i = 0; i < dimb2; ++i) + Mright(i,_) = n->cache.matrix_right[b](i,_) * n->cache.exp_l[b2][i]; + compute_matrix_right(n->left, b, b2, Mright, false, dtau_0); + } + } + + /// Mark the matrix left that need to be recomputed + void impurity_trace::update_matrix_left(node n) { + if (n->key <= max_tau) { + for (int b = 0; b < n_blocks; ++b) n->cache.matrix_left_valid[b] = false; + if (n->left) update_matrix_left(n->left); + } + if (n->right) update_matrix_left(n->right); + } + + /// Mark the matrix right that need to be recomputed + void impurity_trace::update_matrix_right(node n) { + if (n->key >= min_tau) { + for (int b = 0; b < n_blocks; ++b) n->cache.matrix_right_valid[b] = false; + if (n->right) update_matrix_right(n->right); + } + if (n->left) update_matrix_right(n->left); + } + + void impurity_trace::compute_density_matrix(node n, int b, int br, bool is_root, double dtau_beta, double dtau_0) { + + double weight = 0; + double eu = 0; + double ev = 0; + + if (n->left || n->right || is_root) { + int b1 = br ; + if (n->right) b1 = (n->right)->cache.block_table[br]; + int b2 = get_op_block_map(n,b1); + + int dimb = get_block_dim(b); + int dimb1 = get_block_dim(b1); + int dimb2 = get_block_dim(b2); + double dtau = dtau_beta + dtau_0; + double epsilon = 1.e-15; + + // Special case when sampling the density matrix between last + // and first operators. We only do it once, at the root node. + if (is_root) { + for (int u = 0; u < dimb; ++u) { + eu = get_block_eigenval(b,u); + for (int v = 0; v < dimb; ++v) { + ev = get_block_eigenval(b,v); + if (std::abs(eu - ev) < epsilon) + weight = std::exp(-eu * dtau) * dtau / beta; + else + weight = (std::exp(- eu * dtau) - std::exp(- ev * dtau) ) / (beta * (ev - eu)) ; + density_matrix[b].mat(u,v) = density_matrix[b].mat(u,v) + n->cache.matrices[b](u,v) * weight; + } + } + density_matrix[b].is_valid = true; + } + + if (n->left) { + double dtau_l = n->cache.dtau_l; + matrix_t M = {}; + node x = tree.max(n->left); + if (dimb==1 && dimb2==1) + M = n->cache.matrix_right[b](0,0) * x->cache.matrix_left[b2]; + else + M = n->cache.matrix_right[b] * x->cache.matrix_left[b2]; + for (int u = 0; u < dimb2; ++u) { + eu = get_block_eigenval(b2,u); + for (int v = 0; v < dimb2; ++v) { + ev = get_block_eigenval(b2,v); + if (std::abs(eu - ev) < epsilon) + weight = n->cache.exp_l[b2][u] * dtau_l / beta; + else + weight = (n->cache.exp_l[b2][u] - n->cache.exp_l[b2][v] ) / (beta * (ev - eu)) ; + density_matrix[b2].mat(u,v) = density_matrix[b2].mat(u,v) + M(u,v) * weight; + } + } + density_matrix[b2].is_valid = true; + compute_density_matrix(n->left, b, b2, false, dtau_beta, dtau_0); + } + + if (n->right) { + double dtau_r = n->cache.dtau_r; + matrix_t M = {}; + node x = tree.min(n->right); + if (dimb==1 && dimb1==1) + M = x->cache.matrix_right[b](0,0) * n->cache.matrix_left[b1]; + else + M = x->cache.matrix_right[b] * n->cache.matrix_left[b1]; + for (int u = 0; u < dimb1; ++u) { + eu = get_block_eigenval(b1,u); + for (int v = 0; v < dimb1; ++v) { + ev = get_block_eigenval(b1,v); + if (std::abs(eu - ev) < epsilon) + weight = n->cache.exp_r[b1][u] * dtau_r / beta; + else + weight = (n->cache.exp_r[b1][u] - n->cache.exp_r[b1][v] ) / (beta * (ev - eu)) ; + density_matrix[b1].mat(u,v) = density_matrix[b1].mat(u,v) + M(u,v) * weight; + } + } + density_matrix[b1].is_valid = true; + compute_density_matrix(n->right, b, br, false, dtau_beta, dtau_0); + } + } + } + // ------- Update the cache ----------------------- - void impurity_trace::update_cache() { update_cache_impl(tree.get_root()); } + void impurity_trace::update_cache() { + update_cache_impl(tree.get_root()); + if (measure_density_matrix) { + if (tree.get_root() && time_invariance) { + update_matrix_left(tree.get_root()); + update_matrix_right(tree.get_root()); + } + } + } // -------------------------------- @@ -232,6 +462,8 @@ namespace triqs_cthyb { update_cache_impl(n->right); n->cache.dtau_r = (n->right ? double(n->key - tree.min_key(n->right)) : 0); n->cache.dtau_l = (n->left ? double(tree.max_key(n->left) - n->key) : 0); + n->cache.dtau_l_temp = n->cache.dtau_l; + n->cache.dtau_r_temp = n->cache.dtau_r; for (int b = 0; b < n_blocks; ++b) { auto r = compute_block_table_and_bound(n, b, double_max, false); n->cache.block_table[b] = r.first; @@ -248,24 +480,29 @@ namespace triqs_cthyb { if ((n == nullptr) || (!n->modified)) return; update_dtau(n->left); update_dtau(n->right); - n->cache.dtau_r = (n->right ? double(n->key - tree.min_key(n->right)) : 0); - n->cache.dtau_l = (n->left ? double(tree.max_key(n->left) - n->key) : 0); + n->cache.dtau_r_temp = (n->right ? double(n->key - tree.min_key(n->right)) : 0); + n->cache.dtau_l_temp = (n->left ? double(tree.max_key(n->left) - n->key) : 0); } //-------- Compute the full trace ------------------------------------------ // Returns MC atomic weight and reweighting = trace/(atomic weight) - std::pair impurity_trace::compute(double p_yee, double u_yee) { + std::pair impurity_trace::compute(double p_yee, double u_yee, bool meas_den) { double epsilon = 1.e-15; // Machine precision auto log_epsilon0 = -std::log(1.e-15); double lnorm_threshold = double_max - 100; std::vector> init_to_sort_lnorm_b, to_sort_lnorm_b; // pairs of lnorm and b to sort in order of bound + if (meas_den) { + min_tau = time_pt(time_pt::Nmax,beta); + max_tau = time_pt(0,beta); + } + // simplifies later code if (tree_size == 0) { + if (meas_den) density_matrix = atomic_rho; if (use_norm_as_weight) { - density_matrix = atomic_rho; - return {atomic_norm, atomic_z / atomic_norm}; + return {atomic_norm, atomic_z / atomic_norm}; } else return {atomic_z, 1}; } @@ -329,7 +566,19 @@ namespace triqs_cthyb { double norm_trace_sq = 0, trace_abs = 0; // Put density_matrix to "not recomputed" - for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false; + if (meas_den) { + if (time_invariance) { // reset to 0 + for (int bl = 0; bl < n_blocks; ++bl) { + int dim = get_block_dim(bl); + density_matrix[bl].mat = matrix_t(dim,dim); + density_matrix[bl].mat = h_scalar_t{0}; + density_matrix[bl].is_valid = false; + } + } + else { + for (int bl = 0; bl < n_blocks; ++bl) density_matrix[bl].is_valid = false; + } + } auto trace_contrib_block = std::vector>{}; //FIXME complex -- can histos handle this? @@ -340,24 +589,25 @@ namespace triqs_cthyb { // determine at which block we have exceeded the bound and hence can stop. // Can tighten bound on trace by using sqrt(dim(B)) in the case of Frobenius norm only. bound_cumul[n_bl] = 0; - if (!use_norm_as_weight) { + if (!use_norm_as_weight) { for (int bl = n_bl - 1; bl >= 0; --bl) - bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first) * std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second)); + bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first) * get_block_dim(to_sort_lnorm_b[bl].second); } else { - for (int bl = n_bl - 1; bl >= 0; --bl) bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first); + for (int bl = n_bl - 1; bl >= 0; --bl) bound_cumul[bl] = bound_cumul[bl + 1] + std::exp(-to_sort_lnorm_b[bl].first) * + std::sqrt(get_block_dim(to_sort_lnorm_b[bl].second)); } int bl; for (bl = 0; bl < n_bl; ++bl) { // sum over all blocks - // stopping criterion - if ((bl > 0) && (bound_cumul[bl] <= std::abs(full_trace) * epsilon)) break; + // stopping criterion + if ((bl > 0) && (bound_cumul[bl] <= std::abs(full_trace) * epsilon) && !(meas_den && time_invariance)) break; int block_index = to_sort_lnorm_b[bl].second; // index in original (unsorted) order // additionnal Yee quick return criterion if (p_yee >= 0.0) { - auto current_weight = (use_norm_as_weight ? std::sqrt(norm_trace_sq) : full_trace); + auto current_weight = (use_norm_as_weight ? std::sqrt(norm_trace_sq) : full_trace); auto pmax = std::abs(p_yee) * (std::abs(current_weight) + bound_cumul[bl]); if (pmax < u_yee) return {0, 1}; // pmax < u, we can reject } @@ -380,15 +630,23 @@ namespace triqs_cthyb { trace_abs += std::abs(x); } - if (use_norm_as_weight) { // else we are not allowed to compute this matrix, may make no sense + if (use_norm_as_weight || (meas_den && !time_invariance)) { // recompute the density matrix - density_matrix[block_index].is_valid = true; double norm_trace_sq_partial = 0; - auto &mat = density_matrix[block_index].mat; + matrix_t M = {}; + matrix_t *mat; + if (meas_den && !time_invariance) { + mat = &density_matrix[block_index].mat; + density_matrix[block_index].is_valid = true; + } + else { + M = matrix_t(dim,dim); + mat = &M; + } for (int u = 0; u < dim; ++u) { for (int v = 0; v < dim; ++v) { - mat(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v)); - double xx = std::abs(mat(u, v)); + (*mat)(u, v) = b_mat.second(u, v) * std::exp(-dtau_beta * get_block_eigenval(block_index, u) - dtau_0 * get_block_eigenval(block_index, v)); + double xx = std::abs((*mat)(u, v)); norm_trace_sq_partial += xx * xx; } } @@ -396,9 +654,15 @@ namespace triqs_cthyb { // internal check if (std::abs(trace_partial) - 1.0000001 * std::sqrt(norm_trace_sq_partial) * get_block_dim(block_index) > 1.e-15) TRIQS_RUNTIME_ERROR << "|trace| > dim * norm" << trace_partial << " " << std::sqrt(norm_trace_sq_partial) << " " << trace_abs; - auto dev = std::abs(trace_partial - trace(mat)); + auto dev = std::abs(trace_partial - trace(*mat)); if (dev > 1.e-14) TRIQS_RUNTIME_ERROR << "Internal error : trace and density mismatch. Deviation: " << dev; } + if (meas_den && time_invariance) { + matrix_t M = {}; + compute_matrix_left(root, block_index, M, true, dtau_beta); + compute_matrix_right(root, block_index, block_index, M, true, dtau_0); + compute_density_matrix(root, block_index, block_index, true, dtau_beta, dtau_0); + } #ifdef CHECK_MATRIX_BOUNDED_BY_BOUND if (std::abs(trace_partial) > 1.000001 * dim * std::exp(-to_sort_lnorm_b[bl].first)) @@ -438,9 +702,9 @@ namespace triqs_cthyb { } // return {weight, reweighting} - if (!use_norm_as_weight) return {full_trace, 1}; + if (!use_norm_as_weight) return {full_trace, 1}; // else determine reweighting - auto rw = full_trace / norm_trace; + auto rw = full_trace / norm_trace; if (!isfinite(rw)) rw = 1; //FIXME if (!isfinite(rw)) TRIQS_RUNTIME_ERROR << "Atomic correlators : reweight not finite" << full_trace << " "<< norm_trace; return {norm_trace, rw}; diff --git a/c++/triqs_cthyb/impurity_trace.hpp b/c++/triqs_cthyb/impurity_trace.hpp index 9463f77e..fe6bad86 100644 --- a/c++/triqs_cthyb/impurity_trace.hpp +++ b/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 compute(double p_yee = -1, double u_yee = 0); + std::pair 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 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 block_table; // number of blocks limited to 2^15 std::vector> matrices; // partial product of operator/time evolution matrices + std::vector> matrix_left; // product of operator/time evolution matrices with tau >= tau_node + std::vector> matrix_right; // product of operator/time evolution matrices with tau <= tau_node std::vector matrix_lnorms; // -ln(norm(matrix)) std::vector 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 matrix_left_valid; // is matrix_left still valid ? + std::vector matrix_right_valid; // is matrix right still valid ? + std::vector> exp_l; // exp(-dtau_l * Ei) + std::vector> 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 compute_block_table_and_bound(node n, int b, double bound_threshold, bool use_threshold = true); std::pair 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; diff --git a/c++/triqs_cthyb/measures/density_matrix.cpp b/c++/triqs_cthyb/measures/density_matrix.cpp index 5b9a1289..3c8ecd68 100644 --- a/c++/triqs_cthyb/measures/density_matrix.cpp +++ b/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; diff --git a/c++/triqs_cthyb/moves/double_insert.cpp b/c++/triqs_cthyb/moves/double_insert.cpp index a78ab91d..1389d1fc 100644 --- a/c++/triqs_cthyb/moves/double_insert.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/double_remove.cpp b/c++/triqs_cthyb/moves/double_remove.cpp index 4988f5ca..5bd92f62 100644 --- a/c++/triqs_cthyb/moves/double_remove.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/global.cpp b/c++/triqs_cthyb/moves/global.cpp index c0de2cfe..afcb2496 100644 --- a/c++/triqs_cthyb/moves/global.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/insert.cpp b/c++/triqs_cthyb/moves/insert.cpp index 9e2acf6e..c4cc691f 100644 --- a/c++/triqs_cthyb/moves/insert.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/remove.cpp b/c++/triqs_cthyb/moves/remove.cpp index 3c262d3d..96434f1b 100644 --- a/c++/triqs_cthyb/moves/remove.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/moves/shift.cpp b/c++/triqs_cthyb/moves/shift.cpp index 81604e9d..1f7ca02e 100644 --- a/c++/triqs_cthyb/moves/shift.cpp +++ b/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(); diff --git a/c++/triqs_cthyb/parameters.cpp b/c++/triqs_cthyb/parameters.cpp index 40b393a2..1ea39986 100644 --- a/c++/triqs_cthyb/parameters.cpp +++ b/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); diff --git a/c++/triqs_cthyb/parameters.hpp b/c++/triqs_cthyb/parameters.hpp index 65d0b006..c070c551 100644 --- a/c++/triqs_cthyb/parameters.hpp +++ b/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; diff --git a/c++/triqs_cthyb/qmc_data.hpp b/c++/triqs_cthyb/qmc_data.hpp index 12223822..2c3dec2c 100644 --- a/c++/triqs_cthyb/qmc_data.hpp +++ b/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 d) { return real(d); }, delta)), current_sign(1), diff --git a/c++/triqs_cthyb/solver_core.cpp b/c++/triqs_cthyb/solver_core.cpp index e4d951f1..25a63cd4 100644 --- a/c++/triqs_cthyb/solver_core.cpp +++ b/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"); } diff --git a/python/triqs_cthyb/solver_core_desc.py b/python/triqs_cthyb/solver_core_desc.py index 1f199ce6..c0a4d9cb 100644 --- a/python/triqs_cthyb/solver_core_desc.py +++ b/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",