diff --git a/c++/triqs_ctseg/configuration.cpp b/c++/triqs_ctseg/configuration.cpp index 8473674..76f121b 100644 --- a/c++/triqs_ctseg/configuration.cpp +++ b/c++/triqs_ctseg/configuration.cpp @@ -222,6 +222,30 @@ double K_overlap(std::vector const &seglist, tau_t const &tau, bool i return is_c ? result : -result; } +// --------------------------- + +// List of operators containing all colors. +// Time are ordered in decreasing order, in agreement with the whole physic literature. +std::vector colored_ordered_ops(std::vector> const &seglists) { + int c = 0; // index of color + std::vector ops_list; // list of all the operators + for (auto const &seglist : seglists) { + for (auto const &s : seglist) { + ops_list.push_back({s.tau_c, c, false}); + ops_list.push_back({s.tau_cdag, c, true}); + if (is_cyclic(s)) { + ops_list.push_back({tau_t::beta(), c, false}); + ops_list.push_back({tau_t::zero(), c, true}); + } + } + ++c; + } + std::sort(ops_list.begin(), ops_list.end(), [](const colored_ops_t &a, const colored_ops_t &b) { + return b.tau < a.tau; // Note the order of b and a for descending sort + }); + return ops_list; +} + // =================== PRINTING ======================== std::ostream &operator<<(std::ostream &out, std::vector const &sl) { @@ -248,3 +272,11 @@ std::ostream &operator<<(std::ostream &out, configuration_t const &config) { return out; } +// --------------------------- + +std::ostream &operator<<(std::ostream &out, std::vector const &col) { + for (auto const &[i, co] : itertools::enumerate(col)) + out << "\n" << "i: " << i << ", tau: " << co.tau << ", color: " << co.color + << ", " << (co.is_cdag ? "cdag" : "c"); + return out; +} \ No newline at end of file diff --git a/c++/triqs_ctseg/configuration.hpp b/c++/triqs_ctseg/configuration.hpp index 73c0e16..1fec60f 100644 --- a/c++/triqs_ctseg/configuration.hpp +++ b/c++/triqs_ctseg/configuration.hpp @@ -41,6 +41,13 @@ struct segment_t { static segment_t full_line() { return {tau_t::beta(), tau_t::zero()}; } }; +// operator descriptor: time, color, is_cdag[creation/annihilation] +struct colored_ops_t { + tau_t tau; + int color; + bool is_cdag; +}; + // simple alias using vec_seg_iter_t = std::vector::const_iterator; @@ -141,9 +148,15 @@ double K_overlap(std::vector const &seglist, tau_t const &tau_c, tau_ double K_overlap(std::vector const &seglist, tau_t const &tau, bool is_c, gf const &K, int c1, int c2); +// List of operators containing all colors. +std::vector colored_ordered_ops(std::vector> const &seglists); + // =================== PRINTING ======================== std::ostream &operator<<(std::ostream &out, std::vector const &sl); std::ostream &operator<<(std::ostream &out, configuration_t const &config); template <> struct fmt::formatter : ostream_formatter {}; + +std::ostream &operator<<(std::ostream &out, std::vector const &col); +template <> struct fmt::formatter> : ostream_formatter {}; diff --git a/c++/triqs_ctseg/measures.hpp b/c++/triqs_ctseg/measures.hpp index e4f262c..4082f19 100644 --- a/c++/triqs_ctseg/measures.hpp +++ b/c++/triqs_ctseg/measures.hpp @@ -7,3 +7,4 @@ #include "./measures/density.hpp" #include "./measures/sign.hpp" #include "./measures/perturbation_order_histo.hpp" +#include "./measures/state_hist.hpp" diff --git a/c++/triqs_ctseg/measures/state_hist.cpp b/c++/triqs_ctseg/measures/state_hist.cpp new file mode 100644 index 0000000..722ca65 --- /dev/null +++ b/c++/triqs_ctseg/measures/state_hist.cpp @@ -0,0 +1,59 @@ +#include "./state_hist.hpp" +#include "../logs.hpp" + +namespace measures { + + state_hist::state_hist(params_t const &p, work_data_t const &wdata, configuration_t const &config, results_t &results) + : wdata{wdata}, config{config}, results{results} { + + beta = p.beta; + H = nda::zeros(ipow(2, config.n_color())); + } + + // ------------------------------------- + + void state_hist::accumulate(double s) { + + LOG("\n =================== MEASURE ================ \n"); + + /// Measure the state histogram + /** + *Measures the time the impurity spends in a certain atomic state: + * an atomic state is characterized by occupation number for each color, e.g. + * $|n_1 n_2 n_3 n_4\rangle = |0 1 1 0\rangle$ for a model with 4 colors + * + * - the index of a state in the histogram is given by $\sum_i n_i 2^i$ + * + * - the length of the histogram is 2^n_colors + */ + + Z += s; + + double tau_prev = beta; // time of prevous operator; start with beta + nda::vector state = nda::zeros(config.n_color()); + for (auto const &op : colored_ordered_ops(config.seglists)) { + int state_idx = 0; + for (auto c : range(config.n_color())) + if (state(c)) state_idx += ipow(2, c); // get the index of the impurity state + H(state_idx) += (tau_prev - op.tau); + tau_prev = (double)op.tau; + ALWAYS_EXPECTS((state(op.color) == op.is_cdag), "Operator error at color {}", op.color); + state(op.color) = !op.is_cdag; + } + + // get edge state contribution; tau_prev has time of last operator + ALWAYS_EXPECTS((state == nda::zeros(config.n_color())), "Operator error"); + H(0) += tau_prev; + } + // ------------------------------------- + + void state_hist::collect_results(mpi::communicator const &c) { + + Z = mpi::all_reduce(Z, c); + H = mpi::all_reduce(H, c); + H = H / (Z * beta); + + // store the result (not reused later, hence we can move it). + results.state_hist = std::move(H); + } +} // namespace measures \ No newline at end of file diff --git a/c++/triqs_ctseg/measures/state_hist.hpp b/c++/triqs_ctseg/measures/state_hist.hpp new file mode 100644 index 0000000..a55b715 --- /dev/null +++ b/c++/triqs_ctseg/measures/state_hist.hpp @@ -0,0 +1,25 @@ +#pragma once +#include "../configuration.hpp" +#include "../work_data.hpp" +#include "../results.hpp" +#include "../util.hpp" + +namespace measures { + + struct state_hist { + + work_data_t const &wdata; + configuration_t const &config; + results_t &results; + double beta; + + nda::vector H; + + double Z = 0; + + state_hist(params_t const ¶ms, work_data_t const &wdata, configuration_t const &config, results_t &results); + + void accumulate(double s); + void collect_results(mpi::communicator const &c); + }; +} // namespace measures \ No newline at end of file diff --git a/c++/triqs_ctseg/params.hpp b/c++/triqs_ctseg/params.hpp index 6c8873e..527456b 100644 --- a/c++/triqs_ctseg/params.hpp +++ b/c++/triqs_ctseg/params.hpp @@ -108,6 +108,9 @@ struct solve_params_t { /// Whether to measure langle s_x(tau)s_x(0)rangle (see [[measure_sperp_tau]]) bool measure_sperpt = false; + /// Whether to measure state histograms (see [[measure_statehist]]) + bool measure_statehist = false; + /// Hartree shift of the chem pot nda::vector hartree_shift = nda::vector{}; diff --git a/c++/triqs_ctseg/results.cpp b/c++/triqs_ctseg/results.cpp index 42dd620..f7e5120 100644 --- a/c++/triqs_ctseg/results.cpp +++ b/c++/triqs_ctseg/results.cpp @@ -15,6 +15,7 @@ void h5_write(h5::group h5group, std::string subgroup_name, results_t const &c) h5_write(grp, "densities", c.densities); h5_write(grp, "perturbation_order_histo_Delta", c.perturbation_order_histo_Delta); h5_write(grp, "perturbation_order_histo_Jperp", c.perturbation_order_histo_Jperp); + h5_write(grp, "state_hist", c.state_hist); } //------------------------------------ @@ -34,4 +35,5 @@ void h5_read(h5::group h5group, std::string subgroup_name, results_t &c) { h5_read(grp, "densities", c.densities); h5_read(grp, "perturbation_order_histo_Delta", c.perturbation_order_histo_Delta); h5_read(grp, "perturbation_order_histo_Jperp", c.perturbation_order_histo_Jperp); + h5_read(grp, "state_hist", c.state_hist); } diff --git a/c++/triqs_ctseg/results.hpp b/c++/triqs_ctseg/results.hpp index 58af0f6..8f921ef 100644 --- a/c++/triqs_ctseg/results.hpp +++ b/c++/triqs_ctseg/results.hpp @@ -38,6 +38,9 @@ struct results_t { /// Perturbation order histogram std::optional perturbation_order_histo_Jperp; + /// State histogram + std::optional> state_hist; + /// Average sign double sign; }; diff --git a/c++/triqs_ctseg/solver_core.cpp b/c++/triqs_ctseg/solver_core.cpp index b57832b..2a72c00 100644 --- a/c++/triqs_ctseg/solver_core.cpp +++ b/c++/triqs_ctseg/solver_core.cpp @@ -93,6 +93,7 @@ void solver_core::solve(solve_params_t const &solve_params) { if (p.measure_sperpt) CTQMC.add_measure(measures::sperp_tau{p, wdata, config, results}, "(tau)"); if (p.measure_perturbation_order_histograms) CTQMC.add_measure(measures::perturbation_order_histo{p, wdata, config, results}, "Perturbation orders"); + if (p.measure_statehist) CTQMC.add_measure(measures::state_hist{p, wdata, config, results}, "State histograms"); // Run and collect results CTQMC.warmup_and_accumulate(p.n_warmup_cycles, p.n_cycles, p.length_cycle, diff --git a/c++/triqs_ctseg/util.hpp b/c++/triqs_ctseg/util.hpp index 4374682..5f72313 100644 --- a/c++/triqs_ctseg/util.hpp +++ b/c++/triqs_ctseg/util.hpp @@ -26,3 +26,8 @@ long det_lower_bound_x(auto const &d, auto const &x) { long det_lower_bound_y(auto const &d, auto const &y) { return lower_bound([&d](long i) { return d.get_y(i).first; }, d.size(), y); } + +// Integer power +constexpr unsigned int ipow(unsigned int n, unsigned int m) { + return m == 0 ? 1 : m == 1 ? n : n * ipow(n, m - 1); +} diff --git a/python/triqs_ctseg/parameters_solve_params_t.rst b/python/triqs_ctseg/parameters_solve_params_t.rst index 065fee7..d3fd282 100644 --- a/python/triqs_ctseg/parameters_solve_params_t.rst +++ b/python/triqs_ctseg/parameters_solve_params_t.rst @@ -53,6 +53,8 @@ +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | measure_sperpt | bool | false | Whether to measure langle s_x(tau)s_x(0)rangle (see [[measure_sperp_tau]]) | +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ +| measure_statehist | bool | false | Whether to measure state histograms (see [[measure_statehist]]) | ++---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | hartree_shift | nda::vector | nda::vector{} | Hartree shift of the chem pot | +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | det_init_size | int | 100 | The maximum size of the determinant matrix before a resize | diff --git a/python/triqs_ctseg/solver_core_desc.py b/python/triqs_ctseg/solver_core_desc.py index 17077da..a740700 100644 --- a/python/triqs_ctseg/solver_core_desc.py +++ b/python/triqs_ctseg/solver_core_desc.py @@ -84,6 +84,11 @@ read_only= True, doc = r"""Perturbation order histogram""") +c.add_member(c_name = "state_hist", + c_type = "std::optional>", + read_only= True, + doc = r"""State histogram""") + c.add_member(c_name = "sign", c_type = "double", read_only= True, @@ -181,6 +186,8 @@ +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | measure_sperpt | bool | false | Whether to measure langle s_x(tau)s_x(0)rangle (see [[measure_sperp_tau]]) | +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ +| measure_statehist | bool | false | Whether to measure state histograms (see [[measure_statehist]]) | ++---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | hartree_shift | nda::vector | nda::vector{} | Hartree shift of the chem pot | +---------------------------------------+--------------------------------------+-----------------------------------------+-------------------------------------------------------------------------------------------------------------------+ | det_init_size | int | 100 | The maximum size of the determinant matrix before a resize | @@ -349,6 +356,11 @@ initializer = """ false """, doc = r"""Whether to measure langle s_x(tau)s_x(0)rangle (see [[measure_sperp_tau]])""") +c.add_member(c_name = "measure_statehist", + c_type = "bool", + initializer = """ false """, + doc = r"""Whether to measure state histograms (see [[measure_statehist]])""") + c.add_member(c_name = "hartree_shift", c_type = "nda::vector", initializer = """ nda::vector{} """,