From f884fa4c490e2074bdcaedc35ef091f914e6dd78 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 1/7] Add files via upload --- configuration.hpp | 120 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 120 insertions(+) create mode 100644 configuration.hpp diff --git a/configuration.hpp b/configuration.hpp new file mode 100644 index 00000000..16578f97 --- /dev/null +++ b/configuration.hpp @@ -0,0 +1,120 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "./util.hpp" +#include +#include +#include +#include + +#include
+ +#include + +namespace triqs_cthyb { + + using triqs::utility::time_pt; + using triqs::utility::time_segment; + + // The description of the C operator + struct op_desc { + int block_index; // the block index of the operator + int inner_index; // the inner index inside the block + bool dagger; // is the operator a dagger + long linear_index; // the cumulative index + + friend std::ostream &operator<<(std::ostream &out, op_desc const &op) { + out << (op.dagger ? "Cdag(" : "C(") << op.block_index << "," << op.inner_index << ")"; + return out; + } + + friend void h5_write(h5::group g, op_desc const &op) { + h5_write(g, "block", op.block_index); + h5_write(g, "inner", op.inner_index); + h5_write(g, "dagger", op.dagger); + } + }; + + // The configuration of the Monte Carlo + struct configuration { + + // a map associating an operator to an imaginary time + using oplist_t = std::map>; + +#ifdef SAVE_CONFIGS + configuration(double beta) : beta_(beta), id(0), configs_hfile("configs.h5", exists("configs.h5") ? 'a' : 'w') { + if (NUM_CONFIGS_TO_SAVE > 0) h5_write(configs_hfile, "c_0", *this); + } + ~configuration() { configs_hfile.close(); } +#else + configuration(double beta) : beta_(beta), id(0) {} +#endif + + double beta() const { return beta_; } + int size() const { return oplist.size(); } + oplist_t 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; } + void erase(time_pt const &t) { oplist.erase(t); } + void clear() { oplist.clear(); } + + oplist_t::iterator begin() { return oplist.begin(); } + oplist_t::iterator end() { return oplist.end(); } + oplist_t::const_iterator begin() const { return oplist.begin(); } + oplist_t::const_iterator end() const { return oplist.end(); } + + friend std::ostream &operator<<(std::ostream &out, configuration const &c) { + for (auto const &op : c) out << "tau = " << op.first << " : " << op.second << std::endl; + return out; + } + + // 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) { + h5::group conf_group = conf.create_group(conf_group_name); + 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); + } + } + + long get_id() const { return id; } // Get the id of the current configuration + void finalize() { + id++; +#ifdef SAVE_CONFIGS + if (id < NUM_CONFIGS_TO_SAVE) h5_write(configs_hfile, "c_" + std::to_string(id), *this); +#endif + } + + private: + double beta_; + long id; // configuration id, for debug purposes + oplist_t oplist; + +#ifdef SAVE_CONFIGS + // HDF5 file to save configurations + h5::file configs_hfile; +#endif + }; +} From 93047f1503e1d3eb92ebc9bac615d4ced10e81a3 Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:33:08 +0200 Subject: [PATCH 2/7] Add files via upload --- c++/triqs_cthyb/impurity_trace.cpp | 325 ++++++++++++++++++++++++++--- 1 file changed, 296 insertions(+), 29 deletions(-) diff --git a/c++/triqs_cthyb/impurity_trace.cpp b/c++/triqs_cthyb/impurity_trace.cpp index ec84b14c..acd6b23e 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].size() != dim) n->cache.exp_r[b1] = std::vector(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].size() != dim) n->cache.exp_l[b2] = std::vector(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,224 @@ namespace triqs_cthyb { return {b3, std::move(M)}; } + + void impurity_trace::compute_matrix_left(node n, int b, matrix_t &Mleft, bool is_empty, double dtau_beta) { + + int dimb1 = 0; + int dimb2 = 0; + auto _ = arrays::range(); + + int b1 = b; + if (n->right) b1 = (n->right)->cache.block_table[b]; + int b2 = get_op_block_map(n,b1); + + dimb1 = get_block_dim(b1); + 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); + } + } + + void impurity_trace::compute_matrix_right(node n, int b, int br, matrix_t &Mright, bool is_empty, double dtau_0) { + + int dimb1 = 0; + int dimb2 = 0; + int dimb = 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); + + dimb1 = get_block_dim(b1); + dimb2 = get_block_dim(b2); + 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); + } + } + + 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); + } + + 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; + int dimb = 0; + int dimb1 = 0; + int dimb2 = 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); + + dimb = get_block_dim(b); + dimb1 = get_block_dim(b1); + 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 +465,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 +483,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 +569,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 +592,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; + //if ((bl > 0) && (bound_cumul[bl] <= std::abs(full_trace) * epsilon)) 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 +633,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 +657,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 +705,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}; From bfea411d08dfb8d086531da83adbbb39255d5765 Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:36:03 +0200 Subject: [PATCH 3/7] Add files via upload --- impurity_trace.hpp | 522 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 522 insertions(+) create mode 100644 impurity_trace.hpp diff --git a/impurity_trace.hpp b/impurity_trace.hpp new file mode 100644 index 00000000..0b892394 --- /dev/null +++ b/impurity_trace.hpp @@ -0,0 +1,522 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "./configuration.hpp" +#include "./parameters.hpp" +#include "triqs/utility/rbt.hpp" +#include +#include + +//#define PRINT_CONF_DEBUG + +using namespace triqs; +using histo_map_t = std::map; +using triqs::stat::histogram; +using triqs::utility::rb_tree; +using triqs::utility::rbt_insert_error; + +namespace triqs_cthyb { + + /******************************************** + Calculate the trace of the impurity problem. + ********************************************/ + class impurity_trace { + + 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 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, bool meas_den = false); + + // ------- Configuration and h_loc data ---------------- + + const configuration *config; // config object does exist longer (temporally) than this object. + const atom_diag *h_diag; // access to the diagonalization of h_loc + const int n_orbitals = h_diag->get_fops().size(); // total number of orbital flavours + const int n_blocks = h_diag->n_subspaces(); // + const int n_eigstates = h_diag->get_full_hilbert_space_dim(); // size of the hilbert space + + // ------- Trace data ---------------- + + private: + struct bool_and_matrix { + bool is_valid; + matrix mat; + }; + std::vector density_matrix; // density_matrix, by block, with a bool to say if it has been recomputed + std::vector atomic_rho; // atomic density matrix (non-normalized) + double atomic_z; // atomic partition function + 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 const &get_density_matrix() const { return density_matrix; } + + // ------------------ Cache data ---------------- + + private: + // 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 block_table; // number of blocks limited to 2^15 + std::vector> matrices; // partial product of operator/time evolution matrices + std::vector> matrix_left; + std::vector> matrix_right; + std::vector matrix_lnorms; // -ln(norm(matrix)) + std::vector matrix_norm_valid; // is the norm of the matrix still valid? + std::vector matrix_left_valid; + std::vector matrix_right_valid; + std::vector> exp_l; + std::vector> 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 { + op_desc op; + cache_t cache; + node_data_t(op_desc op, int n_blocks) : op(op), cache(n_blocks) {} + void reset(op_desc op_new) { op = op_new; } + }; + + using rb_tree_t = rb_tree>; + using node = rb_tree_t::node; + +#ifdef EXT_DEBUG + public: +#endif + rb_tree_t tree; // the red black tree and its nodes + + std::vector aux_operators; + + // ---------------- Cache machinery ---------------- + void update_cache(); + + private: + // The dimension of block b + int get_block_dim(int b) const { return h_diag->get_subspace_dim(b); } + + // the i-th eigenvalue of the block b + double get_block_eigenval(int b, int i) const { return h_diag->get_eigenvalue(b, i); } + + // the minimal eigenvalue of the block b + double get_block_emin(int b) const { return get_block_eigenval(b, 0); } + + // node, block -> image of the block by n->op (the operator) + int get_op_block_map(node n, int b) const { + if( n->op.linear_index >= 0 ) + return (n->op.dagger ? h_diag->cdag_connection(n->op.linear_index, b) : h_diag->c_connection(n->op.linear_index, b)); + else { + int aux_idx = -n->op.linear_index - 1; + return aux_operators[aux_idx].connection(b); + } + } + + // the matrix of n->op, from block b to its image + matrix const &get_op_block_matrix(node n, int b) const { + if( n->op.linear_index >= 0 ) + return (n->op.dagger ? h_diag->cdag_matrix(n->op.linear_index, b) : h_diag->c_matrix(n->op.linear_index, b)); + else { + int aux_idx = -n->op.linear_index - 1; + return aux_operators[aux_idx].block_mat[b]; + } + } + + // recursive function for tree traversal + 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); + + bool use_norm_of_matrices_in_cache = true; // When a matrix is computed in cache, its spectral radius replaces the norm estimate + + // integrity check + void check_cache_integrity(bool print = false); + void check_cache_integrity_one_node(node n, bool print); + int check_one_block_table_linear(node n, int b, bool print); // compare block table to that of a linear method (ie. no tree) + matrix_t check_one_block_matrix_linear(node n, int b); // compare matrix to that of a linear method (ie. no tree) + + // Pool of detached nodes + class nodes_storage { + + const int n_blocks; + std::vector nodes; + int i; + + // make a new detached black node + node make_new_node() { return new rb_tree_t::node_t(time_pt{}, node_data_t{{}, n_blocks}, false, 1); } + + public: + inline nodes_storage(int n_blocks, int size = 0) : n_blocks(n_blocks), i(-1) { + for (int j = 0; j < size; ++j) nodes.push_back(make_new_node()); + } + inline ~nodes_storage() { + for (auto &n : nodes) delete n; + } + + // Change the number of stored nodes + inline void reserve(int size) { + if (size > nodes.size()) + for (int j = nodes.size(); j < size; ++j) nodes.push_back(make_new_node()); + } + inline int index() const { return i; } + inline bool is_index_reset() const { return i == -1; } + inline int reset_index() { + int i_ = i; + i = -1; + return i_; + } + inline node swap_next(node n) { + std::swap(n, nodes[++i]); + return n; + } + inline node swap_prev(node n) { + std::swap(n, nodes[i--]); + return n; + } + inline node take_next() { return nodes[++i]; } + inline node take_prev() { return nodes[i--]; } + }; + + public: + + // attach auxiliary operators + op_desc attach_aux_operator(many_body_op_t const &op) { + aux_operators.push_back(h_diag->get_op_mat(op)); + op_desc operator_desc{0, 0, true, -static_cast(aux_operators.size())}; + return operator_desc; + } + + /************************************************************************* + * Ordinary binary search tree (BST) insertion of the trial nodes + *************************************************************************/ + // We have a set of trial nodes, which we can glue, un-glue in the tree at will. + // This avoids allocations. + + int tree_size = 0; // size of the tree +/- the added/deleted node + + // a pool of trial nodes, ready to be glued in the tree. Max 4 to allow for double insertions + nodes_storage trial_nodes = {n_blocks, 4}; + + // for each inserted node, need to know {parent_of_node,child_is_left} + std::vector> inserted_nodes = {{nullptr, false}, {nullptr, false}, {nullptr, false}, {nullptr, false}}; + + node try_insert_impl(node h, node n) { // implementation + if (h == nullptr) return n; + if (h->key == n->key) throw rbt_insert_error{}; + auto smaller = tree.get_comparator()(n->key, h->key); + if (smaller) + h->left = try_insert_impl(h->left, n); + else + h->right = try_insert_impl(h->right, n); + if (inserted_nodes[trial_nodes.index()].first == nullptr) inserted_nodes[trial_nodes.index()] = {h, smaller}; + h->modified = true; + return h; + } + + // unlink all glued trial nodes + void cancel_insert_impl() { + for (int i = 0; i <= trial_nodes.index(); ++i) { + auto &r = inserted_nodes[i]; + if (r.first != nullptr) (r.second ? r.first->left : r.first->right)= nullptr; + } + if (tree_size <= trial_nodes.index() + 1) tree.get_root() = nullptr; + } + + /************************************************************************* + * Node Insertion + *************************************************************************/ + + public: + // Put a trial node at tau for operator op using an ordinary BST insertion (ie. not red black) + void try_insert(time_pt const &tau, op_desc const &op) { + if (trial_nodes.index() > 3) TRIQS_RUNTIME_ERROR << "Error : more than 4 insertions "; + auto &root = tree.get_root(); + node n = trial_nodes.take_next(); // get the next available node + inserted_nodes[trial_nodes.index()] = {nullptr, false}; + n->reset(tau, op); // change the time and op of the node + root = try_insert_impl(root, n); // insert it using a regular BST, no red black + tree_size++; + } + + // Remove all trial nodes from the tree + void cancel_insert() { + cancel_insert_impl(); + trial_nodes.reset_index(); + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + // confirm the insertion of the nodes, with red black balance + void confirm_insert() { + cancel_insert_impl(); // remove BST inserted nodes + // then reinsert the nodes in in balanced RBT + int imax = trial_nodes.reset_index(); + for (int i = 0; i <= imax; ++i) { + node n = trial_nodes.take_next(); + tree.insert(n->key, {n->op, n_blocks}); + } + trial_nodes.reset_index(); + update_cache(); + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + /************************************************************************* + * Node Removal + *************************************************************************/ + private: + std::vector removed_nodes; + std::vector removed_keys; + + public: + // Find and mark as deleted the nth operator with fixed dagger and block_index + // n=0 : first operator, n=1, second, etc... + time_pt try_delete(int n, int block_index, bool dagger) noexcept { + // traverse the tree, looking for the nth operator of the correct dagger, block_index + int i = 0; + node x = find_if(tree, [&](node no) { + if (no->op.dagger == dagger && no->op.block_index == block_index) ++i; + return i == n + 1; + }); + removed_nodes.push_back(x); // store the node + removed_keys.push_back(x->key); // store the key + tree.set_modified_from_root_to(x->key); // mark all nodes on path from node to root as modified + x->delete_flag = true; // mark the node for deletion + tree_size--; + return x->key; + } + + // Clean all the delete flags + void cancel_delete() { + for (auto &n : removed_nodes) n->delete_flag = false; + removed_nodes.clear(); + removed_keys.clear(); + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + // Confirm deletion: the nodes flagged for deletion are truly deleted + void confirm_delete() { + for (auto &k : removed_keys) tree.delete_node(k); // CANNOT use the node here + removed_nodes.clear(); + removed_keys.clear(); + update_cache(); + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + /************************************************************************* + * Node shift (=insertion+deletion) + *************************************************************************/ + + // No try_shift implemented. Use combination of try_insert and try_delete instead. + + // Cancel the shift + void cancel_shift() { + + // Inserted nodes + cancel_insert_impl(); + trial_nodes.reset_index(); + + // Deleted nodes + for (auto &n : removed_nodes) n->delete_flag = false; + removed_nodes.clear(); + removed_keys.clear(); + + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + // Confirm the shift of the node, with red black balance + void confirm_shift() { + + // Inserted nodes + cancel_insert_impl(); // first remove BST inserted nodes + + // then reinsert the nodes used for real in rb tree + int imax = trial_nodes.reset_index(); + for (int i = 0; i <= imax; ++i) { + node n = trial_nodes.take_next(); + tree.insert(n->key, {n->op, n_blocks}); + } + trial_nodes.reset_index(); + + // Deleted nodes + for (auto &k : removed_keys) tree.delete_node(k); // CANNOT use the node here + removed_nodes.clear(); + removed_keys.clear(); + + // update cache only at the end + update_cache(); + tree_size = tree.size(); + tree.clear_modified(); + check_cache_integrity(); + } + + /************************************************************************* + * Node replacement (replace op_desc according to a substitution table) + *************************************************************************/ + private: + // Store copies of the nodes to be replaced + nodes_storage backup_nodes = {n_blocks}; + + 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->right) new_right = try_replace_impl(n->right, updated_ops); + + auto const &op = n->op; + auto it = updated_ops.find(n->key); + bool op_changed = it != updated_ops.end(); + auto const &new_op = op_changed ? it->second : op; + node new_node = n; + if (op_changed || new_left != n->left || new_right != n->right) { + 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) + new_node->reset(key, new_op); + else + new_node->reset(key, op); + new_node->left = new_left; + 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; + } + + node cancel_replace_impl(node n) { + node n_in_tree = n; + if (n_in_tree && n_in_tree->modified) n= backup_nodes.swap_prev(n); + if (n_in_tree->right) cancel_replace_impl(n_in_tree->right); + if (n_in_tree->left) cancel_replace_impl(n_in_tree->left); + return n; + } + + public: + void try_replace(configuration::oplist_t const &updated_ops) { + if (tree_size == 0) return; + + if (!backup_nodes.is_index_reset()) TRIQS_RUNTIME_ERROR << "impurity_trace: improper use of try_replace()"; + backup_nodes.reserve(tree.size()); + auto &root = tree.get_root(); + root = try_replace_impl(root, updated_ops); + } + + void confirm_replace() { + backup_nodes.reset_index(); + update_cache(); + tree.clear_modified(); + check_cache_integrity(); + } + + void cancel_replace() { + if (tree_size == 0 || backup_nodes.is_index_reset()) return; + auto &root = tree.get_root(); + root = cancel_replace_impl(root); + check_cache_integrity(); + } + + private: + // ---------------- Histograms ---------------- + struct histograms_t { + + // How many block non zero at root of the tree + histogram &n_block_at_root; + // how many block kept after the truncation with the bound + histogram &n_block_kept; + + // What is the dominant block in the trace computation ? Sorted by number or energy + histogram &dominant_block_bound; + histogram &dominant_block_trace; + histogram &dominant_block_energy_bound; + histogram &dominant_block_energy_trace; + + // Various ratios : trace/bound, trace/first term of the trace, etc.. + histogram &trace_over_norm; + histogram &trace_abs_over_norm; + histogram &trace_over_trace_abs; + histogram &trace_over_bound; + histogram &trace_first_over_sec_term; + histogram &trace_first_term_trace; + +#define ADD_HISTO(NAME, HISTO) NAME(histos.emplace((#NAME), (HISTO)).first->second) + histograms_t(int n_subspaces, histo_map_t &histos) + : ADD_HISTO(n_block_at_root, histogram(0, n_subspaces)), + ADD_HISTO(n_block_kept, histogram(0, n_subspaces)), + ADD_HISTO(dominant_block_bound, histogram(0, n_subspaces)), + ADD_HISTO(dominant_block_trace, histogram(0, n_subspaces)), + ADD_HISTO(dominant_block_energy_bound, histogram(0, 100, 100)), + ADD_HISTO(dominant_block_energy_trace, histogram(0, 100, 100)), + ADD_HISTO(trace_over_norm, histogram(0, 1.5, 100)), + ADD_HISTO(trace_abs_over_norm, histogram(0, 1.5, 100)), + ADD_HISTO(trace_over_trace_abs, histogram(0, 1.5, 100)), + ADD_HISTO(trace_over_bound, histogram(0, 1.5, 100)), + ADD_HISTO(trace_first_over_sec_term, histogram(0, 1.0, 100)), + ADD_HISTO(trace_first_term_trace, histogram(0, 1.0, 100)) {} +#undef ADD_HISTO + }; + std::unique_ptr histo; + + /// Stream insertion + friend std::ostream &operator<<(std::ostream &out, impurity_trace const & imp_trace); + }; +} // namespace triqs_cthyb From b8fff3bd20de283b58e8f1933d8f71b5a8242c6e Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:38:29 +0200 Subject: [PATCH 4/7] Add files via upload --- parameters.hpp | 272 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) create mode 100644 parameters.hpp diff --git a/parameters.hpp b/parameters.hpp new file mode 100644 index 00000000..59b50471 --- /dev/null +++ b/parameters.hpp @@ -0,0 +1,272 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2017, H. U.R. Strand, M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ + +#pragma once + +#include + +#include "./config.hpp" +#include "./types.hpp" +#include "./configuration.hpp" + +namespace triqs_cthyb { + + using namespace triqs::operators; + using indices_map_t = std::map; + + // All the arguments of the solver_core constructor + struct constr_parameters_t { + + /// Inverse temperature + double beta; + + ///block structure of the gf + gf_struct_t gf_struct; + + /// Number of Matsubara frequencies for gf + int n_iw = 1025; + + /// Number of tau points for gf + int n_tau = 10001; + + /// Number of Legendre polynomials for gf + int n_l = 50; + + /// Use Delta_tau and h_loc0 as input instead of G0_iw? + bool delta_interface = false; + + /// Number of tau points for delta + int n_tau_delta = -1; + + /// Write constr_parameters_t to hdf5 + friend void h5_write(h5::group h5group, std::string subgroup_name, constr_parameters_t const &sp); + + /// Read constr_parameters_t from hdf5 + friend void h5_read(h5::group h5group, std::string subgroup_name, constr_parameters_t &sp); + }; + + // All the arguments of the solve function + struct solve_parameters_t { + + /// Interacting part of the atomic Hamiltonian + /// type: Operator + many_body_op_t h_int; + + /// Number of QMC cycles + long n_cycles; + + /// Partition method + /// type: str + std::string partition_method = "autopartition"; + + /// Quantum numbers + /// type: list(Operator) + /// default: [] + std::vector quantum_numbers = {}; + + /// Restrict local Hilbert space to states with at least this number of particles + /// default: 0 + int loc_n_min = 0; + + /// Restrict local Hilbert space to states with at most this number of particles + /// default: INT_MAX + int loc_n_max = INT_MAX; + + /// Length of a single QMC cycle + /// default: 50 + long length_cycle = 50; + + /// Number of cycles for thermalization + /// default: 5000 + long n_warmup_cycles = 5000; + + /// Seed for random number generator + /// default: 34788 + 928374 * MPI.rank + long random_seed = 34788 + 928374 * mpi::communicator().rank(); + + /// Name of random number generator + /// type: str + std::string random_name = ""; + + /// Maximum runtime in seconds, use -1 to set infinite + /// default: -1 = infinite + long max_time = -1; + + /// Verbosity level + /// default: 3 on MPI rank 0, 0 otherwise. + int verbosity = ((mpi::communicator().rank() == 0) ? 3 : 0); // silence the slave nodes + + /// Add shifting an operator as a move? + bool move_shift = true; + + /// Add double insertions as a move? + bool move_double = true; + + /// Calculate the full trace or use an estimate? + bool use_trace_estimator = false; + + /// Measure G(tau)? :math:`G_{ij}(\tau)=G_{ji}^*(\tau)` is enforced for the resulting G(tau) + bool measure_G_tau = true; + + /// Measure G_l (Legendre)? Note, no hermiticity in G_l is ensured + bool measure_G_l = false; + + /// Measure O_tau by insertion + std::optional> measure_O_tau = {}; + + /// Minumum of operator insertions in: O_tau by insertion measure + int measure_O_tau_min_ins = 10; + + /// Measure G^4(tau,tau',tau'') with three fermionic times. + bool measure_G2_tau = false; + + /// Measure G^4(inu,inu',inu'') with three fermionic frequencies. + bool measure_G2_iw = false; + + /// Measure G^4(inu,inu',inu'') with three fermionic frequencies. + bool measure_G2_iw_nfft = false; + + /// Measure G^4(iomega,inu,inu') within the particle-particle channel. + bool measure_G2_iw_pp = false; + + /// Measure G^4(iomega,inu,inu') within the particle-particle channel. + bool measure_G2_iw_pp_nfft = false; + + /// Measure G^4(iomega,inu,inu') within the particle-hole channel. + bool measure_G2_iw_ph = false; + + /// Measure G^4(iomega,inu,inu') within the particle-hole channel. + bool measure_G2_iw_ph_nfft = false; + + /// Measure G^2(iomega,l,l') within the particle-particle channel. + bool measure_G2_iwll_pp = false; + + /// Measure G^2(iomega,l,l') within the particle-hole channel. + bool measure_G2_iwll_ph = false; + + /// Order of block indices in the definition of G^2. + block_order measure_G2_block_order = block_order::AABB; + + /// List of block index pairs of G^2 to measure. + /// default: measure all blocks + std::set> measure_G2_blocks = {}; + + /// Number of imaginary time slices for G^4 measurement. + int measure_G2_n_tau = 10; + + /// Number of bosonic Matsubara frequencies for G^4 measurement. + int measure_G2_n_bosonic = 30; + + /// Number of fermionic Matsubara frequencies for G^4 measurement. + int measure_G2_n_fermionic = 30; + + /// Number of Legendre coefficients for G^4(iomega,l,l') measurement. + int measure_G2_n_l = 20; + + /// NFFT buffer size for G^4(iomega,l,l') measurement. + int measure_G2_iwll_nfft_buf_size = 100; + + /// NFFT buffer sizes for different blocks + /// default: 100 for every block + std::map nfft_buf_sizes = {}; + + /// Measure perturbation order? + bool measure_pert_order = false; + + /// 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; + + /// Initial configuration + configuration::oplist_t initial_config = {}; + + /// Number of bins for the histograms + int nbins_histo = 100; + + /// Histogram for insert moves + std::vector> hist_insert; + + /// Histogram for remove moves + std::vector> hist_remove; + + /// Use the norm of the density matrix in the weight if true, otherwise use Trace + bool use_norm_as_weight = false; + + /// Analyse performance of trace computation with histograms (developers only)? + bool performance_analysis = false; + + /// Operator insertion/removal probabilities for different blocks + /// type: dict(str:float) + /// default: {} + std::map proposal_prob = {}; + + /// List of global moves (with their names). + /// Each move is specified with an index substitution dictionary. + /// type: dict(str : dict(indices : indices)) + /// default: {} + std::map move_global = {}; + + /// Overall probability of the global moves + double move_global_prob = 0.05; + + /// Threshold below which imaginary components of Delta and h_loc are set to zero + double imag_threshold = 1.e-13; + + /// The maximum size of the determinant matrix before a resize + int det_init_size = 100; + + /// Max number of ops before the test of deviation of the det, M^-1 is performed. + int det_n_operations_before_check = 100; + + /// Threshold for determinant precision warnings + double det_precision_warning = 1.e-8; + + /// Threshold for determinant precision error + double det_precision_error = 1.e-5; + + /// Bound for the determinant matrix being singular, abs(det) > singular_threshold. If <0, it is !isnormal(abs(det)) + double det_singular_threshold = -1; + + solve_parameters_t() {} + + solve_parameters_t(many_body_op_t h_int, long n_cycles) : h_int(h_int), n_cycles(n_cycles) {} + + /// Write solve_parameters_t to hdf5 + friend void h5_write(h5::group h5group, std::string subgroup_name, solve_parameters_t const &sp); + + /// Read solve_parameters_t from hdf5 + friend void h5_read(h5::group h5group, std::string subgroup_name, solve_parameters_t &sp); + + /// Threshold below which which off diagonal components of hloc are set to 0 + double off_diag_threshold = 0.0; + + /// Quadratic part of the local Hamiltonian. Must be provided if the Delta interface is used + std::optional h_loc0 = {}; + }; + + /// A struct combining both constr_params_t and solve_params_t + struct params_t : constr_parameters_t, solve_parameters_t { + params_t(constr_parameters_t constr_parameters_, solve_parameters_t solve_parameters_) + : constr_parameters_t(constr_parameters_), solve_parameters_t(solve_parameters_) {} + }; +} From 9a83248b1f3f7835c92472cb090cf060e1f426ca Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:46:49 +0200 Subject: [PATCH 5/7] Add files via upload --- qmc_data.hpp | 245 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 245 insertions(+) create mode 100644 qmc_data.hpp diff --git a/qmc_data.hpp b/qmc_data.hpp new file mode 100644 index 00000000..252848b5 --- /dev/null +++ b/qmc_data.hpp @@ -0,0 +1,245 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "impurity_trace.hpp" +#include +#include +#include +#include "container_set.hpp" + +namespace triqs_cthyb { + using namespace triqs::gfs; + using namespace triqs::mesh; + using namespace nda; + + /************************ + * The Monte Carlo data + ***********************/ + struct qmc_data { + + configuration config; // Configuration + time_segment tau_seg; + std::map, int> linindex; // Linear index constructed from block and inner indices + atom_diag const &h_diag; // Diagonalization of the atomic problem + mutable impurity_trace imp_trace; // Calculator of the trace + std::vector n_inner; + block_gf delta; // Hybridization function + bool updated; + double n_acc; + + /// This callable object adapts the Delta function for the call of the det. + struct delta_block_adaptor { + gf delta_block; // make a copy. Needed in the real case anyway. + + delta_block_adaptor(gf_const_view delta_block) : delta_block(std::move(delta_block)) {} + delta_block_adaptor(delta_block_adaptor const &) = default; + delta_block_adaptor(delta_block_adaptor &&) = default; + delta_block_adaptor &operator=(delta_block_adaptor const &) = delete; + delta_block_adaptor &operator=(delta_block_adaptor &&) = default; + + det_scalar_t operator()(std::pair const &x, std::pair const &y) const { + double dtau = double(x.first - y.first); + auto ind = delta_block.mesh().to_data_index(closest_mesh_pt(dtau)); + det_scalar_t res; + bool use_linear_interpolation = false; + if (use_linear_interpolation) { + double xa, xb, ya, yb; + xa = delta_block.mesh().to_value(ind); + ya = delta_block[ind](x.second, y.second); + if (dtau >= xa) { + xb = delta_block.mesh().to_value(ind+1); + yb = delta_block[ind+1](x.second, y.second); + } + else { + xb = xa; + xa = delta_block.mesh().to_value(ind-1); + yb = ya; + ya = delta_block[ind-1](x.second, y.second); + } + res = ya + (dtau - xa) * (yb - ya) / (xb - xa); + } + else + res = delta_block[ind](x.second, y.second); + return (x.first >= y.first ? res : -res); // x,y first are time_pt, wrapping is automatic in the - operation, but need to + // compute the sign + } + + friend void swap(delta_block_adaptor &dba1, delta_block_adaptor &dba2) noexcept { std::swap(dba1.delta_block, dba2.delta_block); } + }; + + std::vector> dets; // The determinants + int current_sign, old_sign; // Permutation prefactor + h_scalar_t atomic_weight; // The current value of the trace or norm + h_scalar_t atomic_reweighting; // The current value of the reweighting + + // Construction + qmc_data(double beta, solve_parameters_t const &p, atom_diag const &h_diag, std::map, int> linindex, + block_gf_const_view delta, std::vector n_inner, histo_map_t *histo_map) + : config(beta), + 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.time_invariance, p.performance_analysis), + n_inner(n_inner), + delta(map([](gf_const_view d) { return real(d); }, delta)), + updated(false), + n_acc(0.), + current_sign(1), + old_sign(1) { + std::vector>> X(delta.size()), Y(delta.size()); + for (auto const &o : p.initial_config) { + 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 config"; + + 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]), 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"; + if (p.verbosity >= 2) { + std::cerr << "WARNING: The Delta(tau) block number " << bl << " is not real in tau space\n"; + std::cerr << "WARNING: max(Im[Delta(tau)]) = " << max_element(abs(imag(delta[bl].data()))) << "\n"; + std::cerr << "WARNING: Dissregarding the imaginary component in the calculation.\n"; + } + } + 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 + qmc_data &operator=(qmc_data const &) = delete; + + void update_sign() { + + int s = 0; + size_t num_blocks = dets.size(); + std::vector n_op_with_a_equal_to(num_blocks, 0), n_ndag_op_with_a_equal_to(num_blocks, 0); + + // In this first part we compute the sign to bring the configuration to + // d^_1 d^_1 d^_1 ... d_1 d_1 d_1 d^_2 d^_2 ... d_2 d_2 ... d^_n .. d_n + + // loop over the operators "op" in the trace (right to left) + for (auto const &op : config) { + + // how many operators with an 'a' larger than "op" are there on the left of "op"? + for (int a = op.second.block_index + 1; a < num_blocks; ++a) s += n_op_with_a_equal_to[a]; + n_op_with_a_equal_to[op.second.block_index]++; + + // if "op" is not a dagger how many operators of the same a but with a dagger are there on his right? + if (op.second.dagger) + s += n_ndag_op_with_a_equal_to[op.second.block_index]; + else + n_ndag_op_with_a_equal_to[op.second.block_index]++; + } + + // Now we compute the sign to bring the configuration to + // d_1 d^_1 d_1 d^_1 ... d_1 d^_1 ... d_n d^_n ... d_n d^_n + for (int block_index = 0; block_index < num_blocks; block_index++) { + int n = dets[block_index].size(); + s += n * (n + 1) / 2; + } + + old_sign = current_sign; + current_sign = (s % 2 == 0 ? 1 : -1); + } + }; + + //--------- DEBUG --------- + + using det_type = det_manip::det_manip; + + // Print taus of operator sequence in dets + inline void print_det_sequence(qmc_data const &data) { + int i; + int block_index; + for (block_index = 0; block_index < data.dets.size(); ++block_index) { + auto det = data.dets[block_index]; + if (det.size() == 0) return; + std::cout << "BLOCK = " << block_index << std::endl; + for (i = 0; i < det.size(); ++i) { // c_dag + std::cout << " ic_dag = " << i << ": tau = " << det.get_x(i).first << std::endl; + } + for (i = 0; i < det.size(); ++i) { // c + std::cout << " ic = " << i << ": tau = " << det.get_y(i).first << std::endl; + } + } + } + + inline void print_det_sequence(det_type const &det) { + int i; + for (i = 0; i < det.size(); ++i) { // c_dag + std::cout << " ic_dag = " << i << ": tau = " << det.get_x(i).first << std::endl; + } + for (i = 0; i < det.size(); ++i) { // c + std::cout << " ic = " << i << ": tau = " << det.get_y(i).first << std::endl; + } + } + + // Check if dets are correctly ordered, otherwise complain + inline void check_det_sequence(det_type const &det, int config_id) { + if (det.size() == 0) return; + auto tau = det.get_x(0).first; + for (auto ic_dag = 0; ic_dag < det.size(); ++ic_dag) { // c_dag + if (tau < det.get_x(ic_dag).first) { + std::cout << "ERROR in det order in config " << config_id << std::endl; + print_det_sequence(det); + TRIQS_RUNTIME_ERROR << "Det ordering wrong: tau(ic_dag = " << ic_dag << ") = " << double(det.get_x(ic_dag).first); + } + tau = det.get_x(ic_dag).first; + } + tau = det.get_y(0).first; + for (auto ic = 0; ic < det.size(); ++ic) { // c + if (tau < det.get_y(ic).first) { + std::cout << "ERROR in det order in config " << config_id << std::endl; + print_det_sequence(det); + TRIQS_RUNTIME_ERROR << "Det ordering wrong: tau(ic = " << ic << ") = " << double(det.get_y(ic).first); + } + tau = det.get_y(ic).first; + } + } +} // namespace triqs_cthyb From 51c6bacf8d2195de42a59a16921807c037f043dd Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:50:08 +0200 Subject: [PATCH 6/7] Add files via upload --- solver_core.cpp | 490 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 490 insertions(+) create mode 100644 solver_core.cpp diff --git a/solver_core.cpp b/solver_core.cpp new file mode 100644 index 00000000..a84a6502 --- /dev/null +++ b/solver_core.cpp @@ -0,0 +1,490 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014, P. Seth, I. Krivenko, M. Ferrero and O. Parcollet + * Copyright (C) 2017, H. UR Strand, P. Seth, I. Krivenko, + * M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ +#include "./solver_core.hpp" +#include "./qmc_data.hpp" +#include "./configuration.hpp" + +#include +#include +#include +#include +#include +#include + +#include "./moves/insert.hpp" +#include "./moves/remove.hpp" +#include "./moves/double_insert.hpp" +#include "./moves/double_remove.hpp" +#include "./moves/shift.hpp" +#include "./moves/global.hpp" +#include "./measures/G_tau.hpp" +#include "./measures/G_l.hpp" +#include "./measures/O_tau_ins.hpp" +#include "./measures/perturbation_hist.hpp" +#include "./measures/density_matrix.hpp" +#include "./measures/average_sign.hpp" +#include "./measures/average_order.hpp" +#include "./measures/auto_corr_time.hpp" +#include "./measures/update_time.hpp" +#ifdef CTHYB_G2_NFFT +#include "./measures/G2_tau.hpp" +#include "./measures/G2_iw.hpp" +#include "./measures/G2_iw_nfft.hpp" +#include "./measures/G2_iwll.hpp" +#endif +#include "./measures/util.hpp" + +namespace triqs_cthyb { + + struct index_visitor { + std::vector indices; + void operator()(int i) { indices.push_back(std::to_string(i)); } + void operator()(std::string s) { indices.push_back(s); } + }; + + 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), n_tau_delta(p.n_tau_delta), delta_interface(p.delta_interface), constr_parameters(p) { + + if (n_tau_delta == -1) n_tau_delta = n_tau; + if (n_tau_delta < 2 * p.n_iw && ! delta_interface) + TRIQS_RUNTIME_ERROR + << "Must use as least twice as many tau points as Matsubara frequencies: n_iw = " << p.n_iw + //<< " but n_tau = " << p.n_tau << "."; + << " but n_tau_delta = " << p.n_tau_delta << "."; + + // Allocate single particle greens functions + if (not delta_interface) _G0_iw = block_gf({beta, Fermion, n_iw}, gf_struct); + _Delta_tau = block_gf({beta, Fermion, n_tau_delta}, gf_struct); + } + + /// ------------------------------------------------------------------------------------------- + + void solver_core::solve(solve_parameters_t const &solve_parameters_) { + + solve_parameters = solve_parameters_; + solve_parameters_t params(solve_parameters_); + + // Merge constr_params and solve_params + //params_t params(constr_parameters, solve_parameters); + + // http://patorjk.com/software/taag/#p=display&f=Calvin%20S&t=TRIQS%20cthyb + if (params.verbosity >= 2) + std::cout << "\n" + "╔╦╗╦═╗╦╔═╗ ╔═╗ ┌─┐┌┬┐┬ ┬┬ ┬┌┐ \n" + " ║ ╠╦╝║║═╬╗╚═╗ │ │ ├─┤└┬┘├┴┐\n" + " ╩ ╩╚═╩╚═╝╚╚═╝ └─┘ ┴ ┴ ┴ ┴ └─┘\n\n"; + + // determine basis of operators to use + fundamental_operator_set fops; + for (auto const &[bl, bl_size] : gf_struct) { + for (auto idx : range(bl_size)) { fops.insert(bl, idx); } + } + + // setup the linear index map + std::map, int> linindex; + int block_index = 0; + for (auto const &[bl, bl_size] : gf_struct) { + int inner_index = 0; + for (auto idx : range(bl_size)) { + linindex[std::make_pair(block_index, inner_index)] = fops[{bl, idx}]; + inner_index++; + } + block_index++; + } + + // Make list of block sizes + std::vector n_inner; + for (auto const &[bl, bl_size] : gf_struct) { n_inner.push_back(bl_size); } + + if (not delta_interface) { + + // ==== Assert that G0_iw fulfills the fundamental property G(iw)[i,j] = G(-iw)*[j,i] ==== + + if (not is_gf_hermitian(_G0_iw.value())) { + if (params.verbosity >= 2) + std::cout << "!-------------------------------------------------------------------------------------------!\n" + "! WARNING: S.G0_iw violates fundamental Green Function property G0(iw)[i,j] = G0(-iw)*[j,i] !\n" + "! Symmetrizing S.G0_iw ... !\n" + "!-------------------------------------------------------------------------------------------!\n\n"; + _G0_iw = make_hermitian(_G0_iw.value()); + } + + // ==== Compute Delta from G0_iw ==== + + // Initialize Delta with iw - inv(G0[iw]) + auto Delta_iw = inverse(_G0_iw.value()); + for (auto bl : range(gf_struct.size())) + for (auto iw : Delta_iw[bl].mesh()) Delta_iw[bl][iw] = iw - Delta_iw[bl][iw]; + + // Compute the constant part of Delta + Delta_infty_vec = map( + // Compute 0th moment of one block + [imag_threshold = params.imag_threshold](gf_const_view d) { + auto [tail, err] = fit_hermitian_tail(d); + if (err > 1e-5) std::cerr << "WARNING: Tail fit in extraction of delta(infty) has large error of: " << err << std::endl; + auto Delta_infty = matrix{tail(0, ellipsis())}; +#ifndef HYBRIDISATION_IS_COMPLEX + double max_imag = max_element(abs(imag(Delta_infty))); + if (max_imag > imag_threshold) + TRIQS_RUNTIME_ERROR << "Largest imaginary element of delta(infty) e.g. of the local part of G0: " << max_imag + << ", is larger than the set parameter imag_threshold " << imag_threshold; +#endif + return Delta_infty; + }, + Delta_iw); + + // Subtract constant part from Delta and perform Fourier transform + for (auto bl : range(gf_struct.size())) { + for (auto iw : Delta_iw[bl].mesh()) Delta_iw[bl][iw] = Delta_iw[bl][iw] - Delta_infty_vec.value()[bl]; + auto [Delta_tail_b, tail_err] = fit_hermitian_tail(Delta_iw[bl]); + _Delta_tau[bl]() = fourier(Delta_iw[bl], Delta_tail_b); + } + + // ==== Compute h_loc ==== + + _h_loc0 = {}; + + // Add non-interacting terms to h_loc + for (auto bl : range(gf_struct.size())) { + for (auto [n1, n2] : Delta_iw[bl].target_indices()) { +#ifdef LOCAL_HAMILTONIAN_IS_COMPLEX + dcomplex e_ij; + double max_imag = max_element(abs(imag(Delta_infty_vec.value()[bl]))); + if (max_imag > params.imag_threshold) + e_ij = Delta_infty_vec.value()[bl](n1, n2); + else + e_ij = Delta_infty_vec.value()[bl](n1, n2).real(); +#else + double e_ij = Delta_infty_vec.value()[bl](n1, n2).real(); +#endif + // Set off diagonal terms to 0 if they are below off_diag_threshold + if (n1 != n2 && abs(Delta_infty_vec.value()[bl](n1, n2)) < params.off_diag_threshold) e_ij = 0.0; + + auto bl_name = gf_struct[bl].first; + _h_loc0 = _h_loc0 + e_ij * c_dag(bl_name, n1) * c(bl_name, n2); + } + } + + } else { // Delta Interface + if (not is_gf_hermitian(_Delta_tau)) { + if (params.verbosity >= 2) + std::cout << "!---------------------------------------------------------------------------------------!\n" + "! WARNING: S.Delta_tau violates fundamental symmetry Delta(tau)[i,j] = Delta(tau)*[j,i] !\n" + "! Symmetrizing S.Delta_tau ... !\n" + "!---------------------------------------------------------------------------------------!\n\n"; + _Delta_tau = make_hermitian(_Delta_tau); + } + if (not params.h_loc0.has_value()) TRIQS_RUNTIME_ERROR << "h_loc0 must be provided when using the Delta interface"; + _h_loc0 = params.h_loc0.value(); + } + + _h_loc = params.h_int + _h_loc0; + +#ifndef HYBRIDISATION_IS_COMPLEX + // Check that diagonal components of Delta_tau are real + for (auto bl : range(gf_struct.size())) { + long bl_size = gf_struct[bl].second; + // Force all diagonal elements to be real + auto _ = range::all; + for (auto [i, j] : product_range(bl_size, bl_size)) { + auto Delta_tau_bl_ij = _Delta_tau[bl].data()(_, i, j); + double max_imag = max_element(abs(imag(Delta_tau_bl_ij))); + if (i == j && max_imag > 1e-10) { + std::cout << "Warning! Delta_tau diagonal term has max imaginary part: " << max_imag << "Disregarding imaginary part \n"; + Delta_tau_bl_ij = real(Delta_tau_bl_ij); + } else if (max_imag < params.imag_threshold) { + Delta_tau_bl_ij = real(Delta_tau_bl_ij); + } + } + } +#endif + + // Report what h_loc we are using + if (params.verbosity >= 2) + std::cout << "The local Hamiltonian of the problem:" << std::endl << _h_loc << std::endl; + // Reset the histograms + _performance_analysis.clear(); + histo_map_t *histo_map = params.performance_analysis ? &_performance_analysis : nullptr; + + // Determine block structure + if (params.partition_method == "autopartition") { + if (params.verbosity >= 2) + std::cout << "Using autopartition algorithm to partition the local Hilbert space" + << std::endl; + if (params.loc_n_min == 0 && params.loc_n_max == INT_MAX) + h_diag = {_h_loc, fops}; + else { + if (params.verbosity >= 2) + std::cout << "Restricting the local Hilbert space to states with [" << params.loc_n_min + << ";" << params.loc_n_max << "] particles" << std::endl; + h_diag = {_h_loc, fops, params.loc_n_min, params.loc_n_max}; + } + } else if (params.partition_method == "quantum_numbers") { + if (params.quantum_numbers.empty()) TRIQS_RUNTIME_ERROR << "No quantum numbers provided."; + if (params.verbosity >= 2) + std::cout << "Using quantum numbers to partition the local Hilbert space" << std::endl; + h_diag = {_h_loc, fops, params.quantum_numbers}; + } else if (params.partition_method == "none") { // give empty quantum numbers list + std::cout << "Not partitioning the local Hilbert space" << std::endl; + h_diag = {_h_loc, fops, std::vector{}}; + } else + TRIQS_RUNTIME_ERROR << "Partition method " << params.partition_method << " not recognised."; + + // FIXME save h_loc to be able to rebuild h_diag in an analysis program. + //if (_comm.rank() ==0) h5_write(h5::file("h_loc.h5",'w'), "h_loc", _h_loc, fops); + + if (params.verbosity >= 2) + std::cout << "Found " << h_diag.n_subspaces() << " subspaces." << std::endl; + + //if (params.performance_analysis) std::ofstream("impurity_blocks.dat") << h_diag; + + // If one is interested only in the atomic problem + if (params.n_warmup_cycles == 0 && params.n_cycles == 0) { + if (params.measure_density_matrix) _density_matrix = atomic_density_matrix(h_diag, beta); + return; + } + + // Initialise Monte Carlo quantities + qmc_data data(beta, params, h_diag, linindex, _Delta_tau, n_inner, histo_map); + auto qmc = + mc_tools::mc_generic(params.random_name, params.random_seed, params.verbosity); + + // -------------------------------------------------------------------------- + // Moves + // -------------------------------------------------------------------------- + + using move_set_type = mc_tools::move_set; + move_set_type inserts(qmc.get_rng()); + move_set_type removes(qmc.get_rng()); + move_set_type double_inserts(qmc.get_rng()); + move_set_type double_removes(qmc.get_rng()); + + auto &delta_names = _Delta_tau.block_names(); + auto get_prob_prop = [¶ms](std::string const &block_name) { + auto f = params.proposal_prob.find(block_name); + return (f != params.proposal_prob.end() ? f->second : 1.0); + }; + + bool use_improved_sampling = true; + std::vector taus_bin; + if (params.hist_insert.empty() || params.hist_remove.empty()) use_improved_sampling = false; + if (use_improved_sampling) { + if (params.hist_insert.size() != params.hist_remove.size() || params.hist_insert.size() != _Delta_tau.size()) + TRIQS_RUNTIME_ERROR << "Inconsistency in hist_insert and hist_remove: the first dimension should be equal " << + "to the number of blocks in gf_struct"; + for (size_t block = 0; block < _Delta_tau.size(); ++block) { + if (params.hist_insert[block].size() != params.hist_remove[block].size() || + params.hist_insert[block].size() != params.nbins_histo) + TRIQS_RUNTIME_ERROR << "Inconsistency in hist_insert and hist_remove: for each block, you need to provide an array " << + "of size nbins_histo"; + } + taus_bin = std::vector(params.nbins_histo + 1); // vector list of the tau corresponding to the beginning and end of each bin + // converts bin step into uint64 + uint64_t step = time_pt::Nmax / (params.nbins_histo - 1); // this is floored, so this ensures that the sum of all bin steps + // does not overshoot Nmax + time_pt step_t = time_pt(step, beta); + taus_bin[0] = time_pt(0, beta); + for (int i = 1; i < params.nbins_histo; ++i) { + if (i==1) + taus_bin[i] = time_pt(step / 2, beta); // first bin is half-sized + else + taus_bin[i] = taus_bin[i-1] + step_t; + } + taus_bin[params.nbins_histo] = time_pt(time_pt::Nmax, beta); // make sure the last tau point exactly reaches beta + } + for (size_t block = 0; block < _Delta_tau.size(); ++block) { + int block_size = _Delta_tau[block].data().shape()[1]; + auto const &block_name = delta_names[block]; + double prop_prob = get_prob_prop(block_name); + inserts.add(move_insert_c_cdag(block, block_size, block_name, data, qmc.get_rng(), histo_map, params.nbins_histo, + params.hist_insert[block], params.hist_remove[block], taus_bin, use_improved_sampling), "Insert Delta_" + block_name, prop_prob); + removes.add(move_remove_c_cdag(block, block_size, block_name, data, qmc.get_rng(), histo_map, params.nbins_histo, + params.hist_insert[block], params.hist_remove[block], taus_bin, use_improved_sampling), "Remove Delta_" + block_name, prop_prob); + if (params.move_double) { + for (size_t block2 = 0; block2 < _Delta_tau.size(); ++block2) { + int block_size2 = _Delta_tau[block2].data().shape()[1]; + auto const &block_name2 = delta_names[block2]; + double prop_prob2 = get_prob_prop(block_name2); + double_inserts.add( + move_insert_c_c_cdag_cdag(block, block2, block_size, block_size2, block_name, + block_name2, data, qmc.get_rng(), histo_map), + "Insert Delta_" + block_name + "_" + block_name2, prop_prob * prop_prob2); + double_removes.add( + move_remove_c_c_cdag_cdag(block, block2, block_size, block_size2, block_name, + block_name2, data, qmc.get_rng(), histo_map), + "Remove Delta_" + block_name + "_" + block_name2, prop_prob * prop_prob2); + } + } + } + + qmc.add_move(std::move(inserts), "Insert two operators", 1.0); + qmc.add_move(std::move(removes), "Remove two operators", 1.0); + if (params.move_double) { + qmc.add_move(std::move(double_inserts), "Insert four operators", 1.0); + qmc.add_move(std::move(double_removes), "Remove four operators", 1.0); + } + + if (params.move_shift) + qmc.add_move(move_shift_operator(data, qmc.get_rng(), histo_map), "Shift one operator", 1.0); + + if (params.move_global.size()) { + move_set_type global(qmc.get_rng()); + for (auto const &mv : params.move_global) { + auto const &name = mv.first; + auto const &substitutions = mv.second; + global.add(move_global(name, substitutions, data, qmc.get_rng()), name, 1.0); + } + qmc.add_move(std::move(global), "Global moves", params.move_global_prob); + } + + // -------------------------------------------------------------------------- + // Measurements + // -------------------------------------------------------------------------- + + // -------------------------------------------------------------------------- + // Two-particle correlators + + G2_measures_t G2_measures(_Delta_tau, gf_struct, params); + +#ifdef CTHYB_G2_NFFT + // Imaginary-time binning + if (params.measure_G2_tau) + qmc.add_measure(measure_G2_tau{G2_tau, data, G2_measures}, + "G2_tau imaginary-time measurement"); + + // NFFT Matsubara frequency measures + + if (params.measure_G2_iw_nfft) + qmc.add_measure(measure_G2_iw_nfft{G2_iw_nfft, data, G2_measures}, + "G2_iw nfft fermionic measurement"); + if (params.measure_G2_iw_pp_nfft) + qmc.add_measure(measure_G2_iw_nfft{G2_iw_pp_nfft, data, G2_measures}, + "G2_iw_pp nfft particle-particle measurement"); + if (params.measure_G2_iw_ph_nfft) + qmc.add_measure(measure_G2_iw_nfft{G2_iw_ph_nfft, data, G2_measures}, + "G2_iw_ph nfft particle-hole measurement"); + + // Direct Matsubara frequency measurement + + if (params.measure_G2_iw) + qmc.add_measure(measure_G2_iw{G2_iw, data, G2_measures}, + "G2_iw fermionic measurement"); + if (params.measure_G2_iw_pp) + qmc.add_measure(measure_G2_iw{G2_iw_pp, data, G2_measures}, + "G2_iw_pp particle-particle measurement"); + if (params.measure_G2_iw_ph) + qmc.add_measure(measure_G2_iw{G2_iw_ph, data, G2_measures}, + "G2_iw_ph particle-hole measurement"); + + // Legendre mixed basis measurements + if (params.measure_G2_iwll_pp) + qmc.add_measure(measure_G2_iwll{G2_iwll_pp, data, G2_measures}, + "G2_iwll_pp Legendre particle-particle measurement"); + if (params.measure_G2_iwll_ph) + qmc.add_measure(measure_G2_iwll{G2_iwll_ph, data, G2_measures}, + "G2_iwll_ph Legendre particle-hole measurement"); +#endif + + // -------------------------------------------------------------------------- + // Single-particle correlators + + if (params.measure_O_tau) { + + const auto &[O1, O2] = *params.measure_O_tau; + auto comm_0 = O1 * O2 - O2 * O1; + auto comm_1 = O1 * _h_loc - _h_loc * O1; + auto comm_2 = O2 * _h_loc - _h_loc * O2; + + if (!comm_0.is_zero() || !comm_1.is_zero() || !comm_2.is_zero()) { + if (params.verbosity >= 2) { + TRIQS_RUNTIME_ERROR << "Error: measure_O_tau, supplied operators does not commute with " + "the local Hamiltonian.\n" + << "[O1, O2] = " << comm_0 << "\n" + << "[O1, H_loc] = " << comm_1 << "\n" + << "[O2, H_loc] = " << comm_2 << "\n"; + } + } + qmc.add_measure( + measure_O_tau_ins{O_tau, data, n_tau, O1, O2, params.measure_O_tau_min_ins, qmc.get_rng()}, + "O_tau insertion measure"); + } + + if (params.measure_G_tau) { + G_tau = block_gf{{beta, Fermion, n_tau}, gf_struct}; + qmc.add_measure(measure_G_tau{data, n_tau, gf_struct, container_set()}, "G_tau measure"); + } + + if (params.measure_G_l) qmc.add_measure(measure_G_l{G_l, data, n_l, gf_struct}, "G_l measure"); + + // Other measurements + if (params.measure_pert_order) { + auto &g_names = _Delta_tau.block_names(); + perturbation_order = histo_map_t{}; + perturbation_order_total = histogram{}; + for (size_t block = 0; block < _Delta_tau.size(); ++block) { + auto const &block_name = g_names[block]; + qmc.add_measure(measure_perturbation_hist(block, data, (*perturbation_order)[block_name]), "Perturbation order (" + block_name + ")"); + } + 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"; + qmc.add_measure(measure_density_matrix{data, _density_matrix}, + "Density Matrix for local static observable"); + } + + qmc.add_measure(measure_average_sign{data, _average_sign}, "Average sign"); + qmc.add_measure(measure_average_order{data, _average_order}, "Average order"); + qmc.add_measure(measure_auto_corr_time{data, _auto_corr_time}, "Auto-correlation time"); + qmc.add_measure(measure_update_time{data, _update_time}, "Update time"); + // -------------------------------------------------------------------------- + + 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), sign); + qmc.collect_results(_comm); + + _final_config = data.config.get_oplist(); + + if (params.verbosity >= 2) { + std::cout << "Average sign: " << _average_sign << std::endl; + std::cout << "Average order: " << _average_order << std::endl; + std::cout << "Auto-correlation time: " << _auto_corr_time << std::endl; + std::cout << "Average update time: " << _update_time << std::endl; + } + + // Copy local (real or complex) G_tau back to complex G_tau + if (G_tau && G_tau_accum) *G_tau = *G_tau_accum; + } +} // namespace triqs_cthyb From 305dec3ce455aae2bc73726033ebe5939f222cb6 Mon Sep 17 00:00:00 2001 From: ec147 <134504305+ec147@users.noreply.github.com> Date: Mon, 8 Jul 2024 13:51:10 +0200 Subject: [PATCH 7/7] Add files via upload --- solver_core.hpp | 232 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 solver_core.hpp diff --git a/solver_core.hpp b/solver_core.hpp new file mode 100644 index 00000000..be08198a --- /dev/null +++ b/solver_core.hpp @@ -0,0 +1,232 @@ + +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2014-2017, H. U.R. Strand, P. Seth, I. Krivenko, + * M. Ferrero and O. Parcollet + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include +#include +#include +#include +#include +#include +#include + +#include "types.hpp" +#include "container_set.hpp" +#include "parameters.hpp" + +#include "configuration.hpp" + +namespace triqs_cthyb { + + /// Core class of the cthyb solver + class solver_core : public container_set_t { + + double beta; // inverse temperature + atom_diag h_diag; // diagonalization of the local problem + gf_struct_t gf_struct; // Block structure of the Green function + many_body_op_t _h_loc; // The local Hamiltonian = h_int + h0 + many_body_op_t _h_loc0; //noninteracting part of the local Hamiltonian + //int n_iw, n_tau, n_l; + int n_iw, n_tau, n_l, n_tau_delta; + bool delta_interface; + + std::vector _density_matrix; // density matrix, when used in Norm mode + mpi::communicator _comm; // define the communicator, here MPI_COMM_WORLD + histo_map_t _performance_analysis; // Histograms used for performance analysis + mc_weight_t _average_sign; // average sign of the QMC + double _average_order; // average perturbation order + double _auto_corr_time; // Auto-correlation time + double _update_time; + int _solve_status; // Status of the solve upon exit: 0 for clean termination, > 0 otherwise. + configuration::oplist_t _final_config; + + // Single-particle Green's function containers + std::optional _G0_iw; // Non-interacting Matsubara Green's function + G_tau_t _Delta_tau; // Imaginary-time Hybridization function + std::optional>> Delta_infty_vec; // Quadratic instantaneous part of G0_iw + + // Return reference to container_set + container_set_t &container_set() { return static_cast(*this); } + container_set_t const &container_set() const { return static_cast(*this); } + + public: + + // Struct containing the parameters relevant for the solver construction + constr_parameters_t constr_parameters; + + // Struct containing the parameters of the last call to the solve method + solve_parameters_t solve_parameters; + + /** + * Construct a CTHYB solver + * + * @param p Set of parameters specific to the CTHYB solver + */ + CPP2PY_ARG_AS_DICT + solver_core(constr_parameters_t const &p); + + // Delete assignement operator because of const members + solver_core(solver_core const &p) = default; + solver_core(solver_core &&p) = default; + solver_core &operator=(solver_core const &p) = delete; + solver_core &operator=(solver_core &&p) = default; + + /** + * Solve method that performs CTHYB calculation + * + * @param p Set of parameters for the CTHYB calculation + */ + CPP2PY_ARG_AS_DICT + void solve(solve_parameters_t const &p); + + /// The local Hamiltonian of the problem: :math:`H_{loc}` used in the last call to ``solve()``. + many_body_op_t const &h_loc() const { return _h_loc; } + + /// The noninteracting part of the local Hamiltonian. + many_body_op_t const &h_loc0() const {return _h_loc0; } + + /// Set of parameters used in the construction of the ``solver_core`` class. + constr_parameters_t last_constr_parameters() const { return constr_parameters; } + + /// Set of parameters used in the last call to ``solve()``. + solve_parameters_t last_solve_parameters() const { return solve_parameters; } + + /// :math:`G_0^{-1}(i\omega_n = \infty)` in Matsubara Frequency. + [[deprecated("Use h_loc0() instead.")]] + std::vector> Delta_infty() { + if (delta_interface) TRIQS_RUNTIME_ERROR << "Delta_infty cannot be accessed when using the Delta interface"; + return Delta_infty_vec.value(); + } + + /// Get a copy of the last container set. + // HACK TO GET CPP2PY TO WRAP THE container_set_t struct. + /* + CPP2PY_ARG_AS_DICT + void set_container_set(container_set_t &cs) { static_cast(*this) = cs; } + container_set_t last_container_set() { return static_cast(*this); } + */ + + /// :math:`\Delta(\tau)` in imaginary time. + block_gf_view Delta_tau() { return _Delta_tau; } + + /// :math:`G_0(i\omega)` in imaginary frequencies. + block_gf_view G0_iw() { + if (delta_interface) TRIQS_RUNTIME_ERROR << "G0_iw cannot be accessed when using the Delta interface"; + return _G0_iw.value(); + } + + /// Atomic :math:`G(\tau)` in imaginary time. + //block_gf_view atomic_gf() const { return ::triqs_cthyb::atomic_gf(h_diag, beta, gf_struct, _Delta_tau[0].mesh().size()); } + + /// Accumulated density matrix. + std::vector const &density_matrix() const { return _density_matrix; } + + /// Diagonalization of :math:`H_{loc}`. + atom_diag const &h_loc_diagonalization() const { return h_diag; } + + /// Histograms related to the performance analysis. + histo_map_t const &get_performance_analysis() const { return _performance_analysis; } + + /// Monte Carlo average sign. + mc_weight_t average_sign() const { return _average_sign; } + + /// Average perturbation order + double average_order() const { return _average_order; } + + /// Auto-correlation time + double auto_corr_time() const { return _auto_corr_time; } + + /// Average update time + double update_time() const { return _update_time; } + + /// Status of the ``solve()`` on exit. + int solve_status() const { return _solve_status; } + + configuration::oplist_t final_config() const { return _final_config; } + + /// is cthyb compiled with support for complex hybridization? + bool hybridisation_is_complex() const { +#ifdef HYBRIDISATION_IS_COMPLEX + return true; +#else + return false; +#endif + } + + /// is cthyb compiled with support for complex local Hamiltonian + bool local_hamiltonian_is_complex() const { +#ifdef LOCAL_HAMILTONIAN_IS_COMPLEX + return true; +#else + return false; +#endif + } + + CPP2PY_IGNORE + static std::string hdf5_format() { return "CTHYB_SolverCore"; } + + // Function that writes the solver_core to hdf5 file + friend void h5_write(h5::group h5group, std::string subgroup_name, solver_core const &s) { + h5::group grp = subgroup_name.empty() ? h5group : h5group.create_group(subgroup_name); + write_hdf5_format(grp, s); + h5_write_attribute(grp, "TRIQS_GIT_HASH", std::string(STRINGIZE(TRIQS_GIT_HASH))); + h5_write_attribute(grp, "CTHYB_GIT_HASH", std::string(STRINGIZE(CTHYB_GIT_HASH))); + h5_write(grp, "container_set", s.container_set()); + h5_write(grp, "constr_parameters", s.constr_parameters); + h5_write(grp, "solve_parameters", s.solve_parameters); + h5_write(grp, "G0_iw", s._G0_iw); + h5_write(grp, "Delta_tau", s._Delta_tau); + + h5_write(grp, "h_diag", s.h_diag); + h5_write(grp, "h_loc", s._h_loc); + h5_write(grp, "density_matrix", s._density_matrix); + h5_write(grp, "average_sign", s._average_sign); + h5_write(grp, "average_order", s._average_order); + h5_write(grp, "auto_corr_time", s._auto_corr_time); + h5_write(grp, "solve_status", s._solve_status); + h5_write(grp, "Delta_infty_vec", s.Delta_infty_vec); + } + + // Function that read all containers to hdf5 file + CPP2PY_IGNORE + static solver_core h5_read_construct(h5::group h5group, std::string subgroup_name) { + h5::group grp = subgroup_name.empty() ? h5group : h5group.open_group(subgroup_name); + auto constr_parameters = h5::h5_read(grp, "constr_parameters"); + auto s = solver_core{constr_parameters}; + h5_read(grp, "container_set", s.container_set()); + h5_read(grp, "solve_parameters", s.solve_parameters); + h5_read(grp, "G0_iw", s._G0_iw); + h5_read(grp, "Delta_tau", s._Delta_tau); + + h5::try_read(grp, "h_diag", s.h_diag); + h5::try_read(grp, "h_loc", s._h_loc); + h5::try_read(grp, "density_matrix", s._density_matrix); + h5::try_read(grp, "average_sign", s._average_sign); + h5::try_read(grp, "average_order", s._average_order); + h5::try_read(grp, "auto_corr_time", s._auto_corr_time); + h5::try_read(grp, "solve_status", s._solve_status); + h5::try_read(grp, "Delta_infty_vec", s.Delta_infty_vec); + + return s; + } + }; +} // namespace triqs_cthyb