From e3ccc4c2db87bc35f4f17514714b0b8874dca6f6 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 17 Feb 2026 19:04:53 -0500 Subject: [PATCH 01/22] bread crumbs toward design + addition to dca_data --- include/dca/phys/dca_data/dca_data.hpp | 27 +++++++++++-------- include/dca/phys/dca_loop/dca_loop.hpp | 7 +++++ .../ctaux/walker/tools/g_tools.hpp | 10 ++++++- .../shared_tools/interpolation/shrink_G0.hpp | 7 ++++- .../interpolation/g0_interpolation.cpp | 2 ++ 5 files changed, 40 insertions(+), 13 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index d1fafc22b..3e7c5f53c 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -198,6 +198,9 @@ class DcaData { func::function, func::dmn_variadic> G0_r_w_cluster_excluded; func::function> G0_r_t_cluster_excluded; + func::function> + disordered_G0_r_t_cluster_excluded; + // Vector NuDmn * RDmn in size that is the disorder configuration func::function orbital_occupancy; @@ -365,19 +368,19 @@ void DcaData::read(const std::string& filename) { if (parameters_.isAccumulatingG4()) { concurrency_.broadcast_object(G_k_w); #ifndef NDEBUG - if (concurrency_.id() == concurrency_.first()) { - std::cout << "broadcasted G_k_w \n"; - } + if (concurrency_.id() == concurrency_.first()) { + std::cout << "broadcasted G_k_w \n"; + } #endif - for (auto& G4_channel : G4_) { + for (auto& G4_channel : G4_) { concurrency_.broadcast_object(G4_channel); #ifndef NDEBUG - if (concurrency_.id() == concurrency_.first()) { - std::cout << "broadcasted G4_channel \n"; - } + if (concurrency_.id() == concurrency_.first()) { + std::cout << "broadcasted G4_channel \n"; + } #endif - } + } } } @@ -555,7 +558,8 @@ void DcaData::initializeH0_and_H_i() { for (int nu2 = 0; nu2 < NuDmn::dmn_size(); ++nu2) for (int nu1 = 0; nu1 < NuDmn::dmn_size(); ++nu1) { if (std::abs(H_interactions(nu1, nu2, r) - H_interactions(nu2, nu1, minus_r)) > 1e-8) { - std::cout << r << " , " << minus_r << " , " << H_interactions(nu1, nu2, r) << " , " << H_interactions(nu2, nu1, minus_r) << "\n"; + std::cout << r << " , " << minus_r << " , " << H_interactions(nu1, nu2, r) << " , " + << H_interactions(nu2, nu1, minus_r) << "\n"; throw(std::logic_error("Double counting is not consistent.")); } } @@ -589,6 +593,7 @@ void DcaData::initialize_G0() { // Compute G0_r_w. math::transform::FunctionTransform::execute(G0_k_w, G0_r_w); Symmetrize::execute(G0_r_w, H_symmetry, true); + // At this point we can apply the disorder configuration // Compute G0_r_t. math::transform::FunctionTransform::execute(G0_k_t, G0_r_t); @@ -659,13 +664,13 @@ void DcaData::initializeSigma(const std::string& filename) { } } concurrency_.broadcast(parameters_.get_chemical_potential()); - #ifndef NDEBUG +#ifndef NDEBUG if (concurrency_.id() == concurrency_.first()) { std::cout << "broadcasted chemical potential: " << parameters_.get_chemical_potential(); } #endif concurrency_.broadcast(Sigma); - #ifndef NDEBUG +#ifndef NDEBUG if (concurrency_.id() == concurrency_.first()) { std::cout << "broadcasted Sigma \n"; } diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index bef3c0658..3aff5f15f 100644 --- a/include/dca/phys/dca_loop/dca_loop.hpp +++ b/include/dca/phys/dca_loop/dca_loop.hpp @@ -240,6 +240,13 @@ void DcaLoop::execute() { perform_cluster_exclusion_step(); + // I guess here is where we can put in the calculation of the + // disordered_g0_r_t + // The question is if all the previous intialization is immutable + // or whether some if will have to be done again or refactored + // If there is disorder do the disorder application here. + // + double L2_Sigma_difference = solve_cluster_problem(dca_iteration_); // returned from cluster_solver::finalize diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g_tools.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g_tools.hpp index e979c1c5d..2ce928d5d 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g_tools.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g_tools.hpp @@ -44,6 +44,13 @@ class G_TOOLS : public G_MATRIX_TOOLS { double get_Gflop(); + /** build the G_matrix for this configuration + * @param[in] full_configuration the full vertex configuration + * @param[in] N + * @param[in] G0 + * @param[inout] G the question is whether this is in out or just + * out. + */ template void build_G_matrix(configuration_type& full_configuration, const dca::linalg::Matrix& N, @@ -174,7 +181,8 @@ void G_TOOLS::build_G_matrix(configuration_type& full_conf G.ptr(0, 0), LD_G, thread_id, stream_id); if constexpr (dca::util::IsComplex_t::value) { - GFLOP += (8.0 * G.nrCols() * G.nrRows() * N.nrCols() + 12.0 * (G.nrCols() * G.nrRows())) * 1.0e-9; + GFLOP += + (8.0 * G.nrCols() * G.nrRows() * N.nrCols() + 12.0 * (G.nrCols() * G.nrRows())) * 1.0e-9; } else { GFLOP += 2. * G.nrCols() * G.nrRows() * N.nrCols() * 1.e-9; diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp index 0190d01e1..19036dc13 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp @@ -43,11 +43,16 @@ template auto shrinkG0(const SpGreensFunction& G0) { func::function, TDmn>> g0_trimmed; const int s = 0; + const spin_diff_limit = 10e-8; for (int b1 = 0; b1 < BDmn::dmn_size(); b1++) for (int b2 = 0; b2 < BDmn::dmn_size(); b2++) for (int r = 0; r < RDmn::dmn_size(); r++) - for (int t = 0; t < TDmn::dmn_size(); t++) + for (int t = 0; t < TDmn::dmn_size(); t++) { + assert(std::abs(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t)) > spin_diff_limit && + ("G0 0,0 and 1,1 spinsectors do not actually have a small difference (" + + std::to_string(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t))) c_str()); g0_trimmed(b1, b2, r, t) = G0(b1, s, b2, s, r, t); + } return g0_trimmed; } diff --git a/src/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.cpp b/src/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.cpp index 2127857c8..fb7557ae3 100644 --- a/src/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.cpp +++ b/src/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.cpp @@ -18,6 +18,8 @@ namespace solver { template void G0Interpolation::initialize(const FunctionProxy& G0_pars_t) { + // As the G0 comes in it is (b1, b2, r)->p, t + beta_ = PositiveTimeDomain::get_elements().back(); n_div_beta_ = Real(PositiveTimeDomain::get_size() - 1) / beta_; From 19135b87af07c87b6fcd4306031158333e583687 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Fri, 20 Feb 2026 16:23:34 -0500 Subject: [PATCH 02/22] disorder parameters added --- .../shared_tools/interpolation/shrink_G0.hpp | 5 +- .../phys/parameters/disorder_parameters.hpp | 108 ++++++++++++++++++ include/dca/phys/parameters/parameters.hpp | 21 ++-- 3 files changed, 124 insertions(+), 10 deletions(-) create mode 100644 include/dca/phys/parameters/disorder_parameters.hpp diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp index 19036dc13..9546dfa9d 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp @@ -43,14 +43,15 @@ template auto shrinkG0(const SpGreensFunction& G0) { func::function, TDmn>> g0_trimmed; const int s = 0; - const spin_diff_limit = 10e-8; + const double spin_diff_limit = 10e-8; for (int b1 = 0; b1 < BDmn::dmn_size(); b1++) for (int b2 = 0; b2 < BDmn::dmn_size(); b2++) for (int r = 0; r < RDmn::dmn_size(); r++) for (int t = 0; t < TDmn::dmn_size(); t++) { assert(std::abs(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t)) > spin_diff_limit && ("G0 0,0 and 1,1 spinsectors do not actually have a small difference (" + - std::to_string(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t))) c_str()); + std::to_string(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t))) + .c_str()); g0_trimmed(b1, b2, r, t) = G0(b1, s, b2, s, r, t); } return g0_trimmed; diff --git a/include/dca/phys/parameters/disorder_parameters.hpp b/include/dca/phys/parameters/disorder_parameters.hpp new file mode 100644 index 000000000..c1d65718b --- /dev/null +++ b/include/dca/phys/parameters/disorder_parameters.hpp @@ -0,0 +1,108 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak (doakpw@ornl.gov) + +#ifndef DCA_PHYS_PARAMETERS_DISORDER_PARAMETERS_HPP +#define DCA_PHYS_PARAMETERS_DISORDER_PARAMETERS_HPP + +#include + +namespace dca { +namespace phys { +namespace params { +// dca::phys::params:: + +/** This class handles the disorder parameters + */ +class DisorderParameters { +public: + DisorderParameters() {} + + template + int getBufferSize(const Concurrency& concurrency) const; + template + void pack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position) const; + template + void unpack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position); + + template + void readWrite(ReaderOrWriter& reader_or_writer); + + int get_num_disorder_configrations() const { + return num_disorder_configurations_; + } + int get_max_disorder_sites() const { + return max_disorder_sites_; + } + double get_disorder_density() const { + return disorder_density_; + } + +private: + int num_disorder_configurations_{0}; + int max_disorder_sites_{1}; + double disorder_density_{0.0}; +}; + +template +int DisorderParameters::getBufferSize(const Concurrency& concurrency) const { + int buffer_size = 0; + + buffer_size += concurrency.get_buffer_size(num_disorder_configurations_); + buffer_size += concurrency.get_buffer_size(max_disorder_sites_); + buffer_size += concurrency.get_buffer_size(disorder_density_); + + return buffer_size; +} + +template +void DisorderParameters::pack(const Concurrency& concurrency, char* buffer, int buffer_size, + int& position) const { + concurrency.pack(buffer, buffer_size, position, num_disorder_configurations_); + concurrency.pack(buffer, buffer_size, position, max_disorder_sites_); + concurrency.pack(buffer, buffer_size, position, disorder_density_); +} + +template +void DisorderParameters::unpack(const Concurrency& concurrency, char* buffer, int buffer_size, + int& position) { + concurrency.unpack(buffer, buffer_size, position, num_disorder_configurations_); + concurrency.unpack(buffer, buffer_size, position, max_disorder_sites_); + concurrency.unpack(buffer, buffer_size, position, disorder_density_); +} + +template +void DisorderParameters::readWrite(ReaderOrWriter& reader_or_writer) { + try { + reader_or_writer.open_group("disorder"); + + try { + reader_or_writer.execute("num-disorder-configurations", num_disorder_configurations_); + } + catch (const std::exception& r_e) { + } + try { + reader_or_writer.execute("max-disorder-sites", max_disorder_sites_); + } + catch (const std::exception& r_e) { + } + try { + reader_or_writer.execute("disorder-density", disorder_density_); + } + catch (const std::exception& r_e) { + } + reader_or_writer.close_group(); + } + catch (const std::exception& r_e) { + } +} + +} // namespace params +} // namespace phys +} // namespace dca + +#endif // DCA_PHYS_PARAMETERS_DISORDER_PARAMETERS_HPP diff --git a/include/dca/phys/parameters/parameters.hpp b/include/dca/phys/parameters/parameters.hpp index 56c4be6ff..f58af57b3 100644 --- a/include/dca/phys/parameters/parameters.hpp +++ b/include/dca/phys/parameters/parameters.hpp @@ -27,6 +27,7 @@ #include "dca/phys/parameters/analysis_parameters.hpp" #include "dca/phys/domains/cluster/cluster_domain_aliases.hpp" #include "dca/phys/parameters/dca_parameters.hpp" +#include "dca/phys/parameters/disorder_parameters.hpp" #include "dca/phys/parameters/domains_parameters.hpp" #include "dca/phys/parameters/double_counting_parameters.hpp" #include "dca/phys/parameters/ed_solver_parameters.hpp" @@ -65,6 +66,7 @@ template class Parameters : public AnalysisParameters, public DcaParameters, + public DisorderParameters, public DomainsParameters, public DoubleCountingParameters, public EdSolverParameters, @@ -165,10 +167,12 @@ template struct CheckParametersNumericTypes : public std::false_type {}; template -struct CheckParametersNumericTypes ::type>::value, bool>> - : public std::true_type {}; +struct CheckParametersNumericTypes< + Parameters, + std::enable_if_t::type>::value, + bool>> : public std::true_type {}; template @@ -202,8 +206,10 @@ Parameters::type>); + static_assert( + std::is_same_v< + typename Parameters::Scalar, + typename dca::util::ScalarSelect::type>); } template void Parameters::broadcast() -{ + NUMTRAITS>::broadcast() { concurrency_.broadcast_object(*this); } From 0aea7913907f8a19d2b31d7cf37b776c5046db89 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 24 Feb 2026 17:50:14 -0500 Subject: [PATCH 03/22] Sometimes asserts could use output --- include/dca/util/message_assert.hpp | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) create mode 100644 include/dca/util/message_assert.hpp diff --git a/include/dca/util/message_assert.hpp b/include/dca/util/message_assert.hpp new file mode 100644 index 000000000..248557c19 --- /dev/null +++ b/include/dca/util/message_assert.hpp @@ -0,0 +1,24 @@ +// Copyright (C) 2023 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak (doakpw@ornl.gov) +// + +#ifndef DCA_UTIL_MESSAGE_ASSERT_HPP +#define DCA_UTIL_MESSAGE_ASSERT_HPP + +#include +#include + +#define MESSAGE_ASSERT(condition, message) \ + do { \ + if (!(condition)) { \ + std::cerr << "Assertion failed: " << #condition << "\nMessage: " << (message) << std::endl; \ + assert(false); \ + } \ + } while (0) + +#endif From a85ae6ed496604ec15d43a92db03ab8e01d0a884 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 24 Feb 2026 17:50:59 -0500 Subject: [PATCH 04/22] assert guard for loss of spin sector symmetry --- .../shared_tools/interpolation/shrink_G0.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp index 9546dfa9d..79a3ec175 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/interpolation/shrink_G0.hpp @@ -19,6 +19,7 @@ #include "dca/phys/domains/quantum/electron_spin_domain.hpp" #include "dca/phys/domains/time_and_frequency/frequency_domain.hpp" #include "dca/phys/domains/time_and_frequency/time_domain.hpp" +#include "dca/util/message_assert.hpp" namespace dca { namespace phys { @@ -48,10 +49,10 @@ auto shrinkG0(const SpGreensFunction& G0) { for (int b2 = 0; b2 < BDmn::dmn_size(); b2++) for (int r = 0; r < RDmn::dmn_size(); r++) for (int t = 0; t < TDmn::dmn_size(); t++) { - assert(std::abs(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t)) > spin_diff_limit && - ("G0 0,0 and 1,1 spinsectors do not actually have a small difference (" + - std::to_string(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t))) - .c_str()); + MESSAGE_ASSERT(std::abs(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t)) < spin_diff_limit, + ("G0 0,0 and 1,1 spinsectors do not actually have a small difference (" + + std::to_string(G0(b1, 0, b2, 0, r, t) - G0(b1, 1, b2, 1, r, t)))); + g0_trimmed(b1, b2, r, t) = G0(b1, s, b2, s, r, t); } return g0_trimmed; From 1798e8ffa6f9271b8eccf52a3cbf2d44c95e8d5b Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 24 Feb 2026 17:51:23 -0500 Subject: [PATCH 05/22] data members for disorder --- include/dca/phys/dca_data/dca_data.hpp | 4 +++- include/dca/phys/dca_loop/dca_loop_data.hpp | 11 ++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index 3e7c5f53c..7e3f1d647 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -1,5 +1,5 @@ // Copyright (C) 2018 ETH Zurich -// Copyright (C) 2018 UT-Battelle, LLC +// Copyright (C) 2026 UT-Battelle, LLC // All rights reserved. // // See LICENSE for terms of usage. @@ -28,6 +28,7 @@ #include "dca/distribution/dist_types.hpp" #include "dca/function/domains.hpp" +#include "dca/function/domains/dmn_variadic.hpp" #include "dca/function/function.hpp" #include "dca/function/util/real_complex_conversion.hpp" #include "dca/io/reader.hpp" @@ -201,6 +202,7 @@ class DcaData { func::function> disordered_G0_r_t_cluster_excluded; // Vector NuDmn * RDmn in size that is the disorder configuration + func::function> disorder_configuration; func::function orbital_occupancy; diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp index 76e128f23..952fd0848 100644 --- a/include/dca/phys/dca_loop/dca_loop_data.hpp +++ b/include/dca/phys/dca_loop/dca_loop_data.hpp @@ -18,6 +18,7 @@ #include "dca/function/function.hpp" #include "dca/io/filesystem.hpp" #include "dca/io/reader.hpp" +#include "dca/phys/domains/cluster/cluster_definitions.hpp" #include "dca/phys/domains/cluster/cluster_domain.hpp" #include "dca/phys/domains/quantum/dca_iteration_domain.hpp" #include "dca/phys/domains/quantum/electron_band_domain.hpp" @@ -47,6 +48,10 @@ class DcaLoopData { using k_DCA = func::dmn_0>; + using r_DCA = + func::dmn_0>; + DcaLoopData(); template @@ -85,6 +90,9 @@ class DcaLoopData { func::function chemical_potential; func::function average_expansion_order; + func::function average_defect_density; + func::function defect_configurations; + int last_completed_iteration = -1; }; @@ -113,7 +121,8 @@ DcaLoopData::DcaLoopData() density("density"), chemical_potential("chemical-potential"), - average_expansion_order("expansion_order") {} + average_expansion_order("expansion_order"), + average_defect_density("defect_density") {} template template From 1d7bda1c8a55175b5a35546cd9259a4f1b3c2783 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Fri, 27 Feb 2026 12:26:28 -0500 Subject: [PATCH 06/22] sufficient disorder params + unit test --- .../phys/parameters/disorder_parameters.hpp | 51 +++++++++++++------ test/unit/phys/parameters/CMakeLists.txt | 1 + 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/include/dca/phys/parameters/disorder_parameters.hpp b/include/dca/phys/parameters/disorder_parameters.hpp index c1d65718b..39758f1e5 100644 --- a/include/dca/phys/parameters/disorder_parameters.hpp +++ b/include/dca/phys/parameters/disorder_parameters.hpp @@ -17,6 +17,13 @@ namespace params { // dca::phys::params:: /** This class handles the disorder parameters + * The default values result in no disorder configurations + * I'd rather have the code know whether the section was present and + * Therefore whether to enable any of the disorder flow at all. + * Even though that is a change to the semantics of the parameters + * I think it may be better design than defaults that result in no + * disorder being used do to zero length loops or branching on + * num_disorder_configurations. */ class DisorderParameters { public: @@ -32,29 +39,34 @@ class DisorderParameters { template void readWrite(ReaderOrWriter& reader_or_writer); - int get_num_disorder_configrations() const { - return num_disorder_configurations_; - } - int get_max_disorder_sites() const { - return max_disorder_sites_; + double get_disorder_potential() const { + return disorder_potential_; } double get_disorder_density() const { return disorder_density_; } + int get_disorder_num_configurations() const { + return disorder_num_configurations_; + } + int get_disorder_max_sites() const { + return disorder_max_sites_; + } private: - int num_disorder_configurations_{0}; - int max_disorder_sites_{1}; + double disorder_potential_{0.0}; double disorder_density_{0.0}; + int disorder_num_configurations_{0}; + int disorder_max_sites_{1}; }; template int DisorderParameters::getBufferSize(const Concurrency& concurrency) const { int buffer_size = 0; - buffer_size += concurrency.get_buffer_size(num_disorder_configurations_); - buffer_size += concurrency.get_buffer_size(max_disorder_sites_); + buffer_size += concurrency.get_buffer_size(disorder_potential_); buffer_size += concurrency.get_buffer_size(disorder_density_); + buffer_size += concurrency.get_buffer_size(disorder_num_configurations_); + buffer_size += concurrency.get_buffer_size(disorder_max_sites_); return buffer_size; } @@ -62,17 +74,19 @@ int DisorderParameters::getBufferSize(const Concurrency& concurrency) const { template void DisorderParameters::pack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position) const { - concurrency.pack(buffer, buffer_size, position, num_disorder_configurations_); - concurrency.pack(buffer, buffer_size, position, max_disorder_sites_); + concurrency.pack(buffer, buffer_size, position, disorder_potential_); concurrency.pack(buffer, buffer_size, position, disorder_density_); + concurrency.pack(buffer, buffer_size, position, disorder_num_configurations_); + concurrency.pack(buffer, buffer_size, position, disorder_max_sites_); } template void DisorderParameters::unpack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position) { - concurrency.unpack(buffer, buffer_size, position, num_disorder_configurations_); - concurrency.unpack(buffer, buffer_size, position, max_disorder_sites_); + concurrency.unpack(buffer, buffer_size, position, disorder_potential_); concurrency.unpack(buffer, buffer_size, position, disorder_density_); + concurrency.unpack(buffer, buffer_size, position, disorder_num_configurations_); + concurrency.unpack(buffer, buffer_size, position, disorder_max_sites_); } template @@ -81,17 +95,22 @@ void DisorderParameters::readWrite(ReaderOrWriter& reader_or_writer) { reader_or_writer.open_group("disorder"); try { - reader_or_writer.execute("num-disorder-configurations", num_disorder_configurations_); + reader_or_writer.execute("potential", disorder_potential_); + } + catch (const std::exception& r_e) { + } + try { + reader_or_writer.execute("density", disorder_density_); } catch (const std::exception& r_e) { } try { - reader_or_writer.execute("max-disorder-sites", max_disorder_sites_); + reader_or_writer.execute("num-configurations", disorder_num_configurations_); } catch (const std::exception& r_e) { } try { - reader_or_writer.execute("disorder-density", disorder_density_); + reader_or_writer.execute("max-sites", disorder_max_sites_); } catch (const std::exception& r_e) { } diff --git a/test/unit/phys/parameters/CMakeLists.txt b/test/unit/phys/parameters/CMakeLists.txt index 066c09b83..6e31c786a 100644 --- a/test/unit/phys/parameters/CMakeLists.txt +++ b/test/unit/phys/parameters/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(analysis_parameters) add_subdirectory(dca_parameters) +add_subdirectory(disorder_parameters) add_subdirectory(domains_parameters) add_subdirectory(double_counting_parameters) add_subdirectory(ed_solver_parameters) From 913795131db82ce270698951fe36150294aac38c Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Fri, 27 Feb 2026 13:14:38 -0500 Subject: [PATCH 07/22] finish adding disorder_params to parameters.hpp add data member to dca_loop data, should still be more --- include/dca/phys/dca_loop/dca_loop_data.hpp | 5 ++++- include/dca/phys/parameters/parameters.hpp | 5 +++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp index 952fd0848..69035b6ff 100644 --- a/include/dca/phys/dca_loop/dca_loop_data.hpp +++ b/include/dca/phys/dca_loop/dca_loop_data.hpp @@ -91,7 +91,10 @@ class DcaLoopData { func::function average_expansion_order; func::function average_defect_density; - func::function defect_configurations; + func::function num_defect_configurations; + // What should the data structure for the defect configurations look + // I think the same as the r_DCA domain with just a +/- 1 likely a + // Real to make it easy multiply with the disorder potential int last_completed_iteration = -1; }; diff --git a/include/dca/phys/parameters/parameters.hpp b/include/dca/phys/parameters/parameters.hpp index f58af57b3..89e5c6dcc 100644 --- a/include/dca/phys/parameters/parameters.hpp +++ b/include/dca/phys/parameters/parameters.hpp @@ -180,6 +180,7 @@ Parameters::Parameters(const std::string& version_stamp, concurrency_type& concurrency) : AnalysisParameters(Model::DIMENSION), DcaParameters(Model::BANDS), + DisorderParameters(), DomainsParameters(Model::DIMENSION), DoubleCountingParameters(), EdSolverParameters(), @@ -456,6 +457,10 @@ void Parameters::readWrite(reader_or_writer); From f8cc0aa29c6ee4cb476b5314559f6612e6b471a8 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Fri, 27 Feb 2026 18:30:49 -0500 Subject: [PATCH 08/22] wanting to get better output on electron bands dmn Unfortunate rippling --- applications/dca/CMakeLists.txt | 3 +- applications/dca/main_dca.cpp | 8 ++-- include/dca/linalg/util/gpu_event.hpp | 8 ++-- .../dca_loop/disorder/.#apply_disorder.cpp | 1 + .../dca_loop/disorder/.#apply_disorder.hpp | 1 + .../phys/dca_loop/disorder/apply_disorder.cpp | 11 ++++++ .../phys/dca_loop/disorder/apply_disorder.hpp | 17 +++++++++ .../domains/cluster/cluster_operations.hpp | 9 +++-- .../domains/quantum/electron_band_domain.hpp | 34 ++++++++++++----- .../dca/phys/parameters/main_parameters.hpp | 12 ++++++ src/phys/CMakeLists.txt | 2 +- src/phys/domains/quantum/CMakeLists.txt | 1 + .../domains/quantum/electron_band_domain.cpp | 24 ++++++++++++ src/phys/parameters/CMakeLists.txt | 4 ++ src/phys/parameters/main_parameters.cpp | 9 +++++ test/unit/math/nfft/CMakeLists.txt | 2 +- .../accumulation/sp/CMakeLists.txt | 8 ++-- .../accumulation/tp/ndft/CMakeLists.txt | 2 +- .../disorder_parameters/CMakeLists.txt | 6 +++ .../disorder_parameters_test.cpp | 38 +++++++++++++++++++ .../disorder_parameters/input_read_all.json | 8 ++++ 21 files changed, 179 insertions(+), 29 deletions(-) create mode 120000 include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp create mode 120000 include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp create mode 100644 include/dca/phys/dca_loop/disorder/apply_disorder.cpp create mode 100644 include/dca/phys/dca_loop/disorder/apply_disorder.hpp create mode 100644 include/dca/phys/parameters/main_parameters.hpp create mode 100644 src/phys/domains/quantum/electron_band_domain.cpp create mode 100644 src/phys/parameters/CMakeLists.txt create mode 100644 src/phys/parameters/main_parameters.cpp create mode 100644 test/unit/phys/parameters/disorder_parameters/CMakeLists.txt create mode 100644 test/unit/phys/parameters/disorder_parameters/disorder_parameters_test.cpp create mode 100644 test/unit/phys/parameters/disorder_parameters/input_read_all.json diff --git a/applications/dca/CMakeLists.txt b/applications/dca/CMakeLists.txt index f7a9c3762..94bbb9d39 100644 --- a/applications/dca/CMakeLists.txt +++ b/applications/dca/CMakeLists.txt @@ -7,8 +7,9 @@ if (DCA_BUILD_DCA) if (DCA_HAVE_GPU) target_link_libraries(main_dca PRIVATE ${DCA_KERNEL_LIBS}) endif() - + target_link_libraries(main_dca PUBLIC FFTW::Double signals ${DCA_LIBS} dca_io) install(TARGETS main_dca RUNTIME DESTINATION bin) + endif() diff --git a/applications/dca/main_dca.cpp b/applications/dca/main_dca.cpp index 50fcf77af..cba0d97b8 100644 --- a/applications/dca/main_dca.cpp +++ b/applications/dca/main_dca.cpp @@ -15,6 +15,7 @@ #include #include "dca/config/dca.hpp" +#include "dca/phys/parameters/main_parameters.hpp" #include "dca/application/dca_loop_dispatch.hpp" #include "dca/config/cmake_options.hpp" #include "dca/config/haves_defines.hpp" @@ -63,17 +64,16 @@ int dca_main(int argc, char** argv) { << std::endl; } - // Create the parameters object from the input file. ParametersType parameters(dca::util::GitVersion::string(), concurrency); parameters.read_input_and_broadcast(input_file); - if(concurrency.id() == concurrency.first()) + if (concurrency.id() == concurrency.first()) std::cout << "Input read and broadcast.\n"; parameters.update_model(); - if(concurrency.id() == concurrency.first()) + if (concurrency.id() == concurrency.first()) std::cout << "Model updated.\n"; parameters.update_domains(); - if(concurrency.id() == concurrency.first()) + if (concurrency.id() == concurrency.first()) std::cout << "Domains updated.\n"; dca::DistType distribution = parameters.get_g4_distribution(); diff --git a/include/dca/linalg/util/gpu_event.hpp b/include/dca/linalg/util/gpu_event.hpp index bb7614616..891d59199 100644 --- a/include/dca/linalg/util/gpu_event.hpp +++ b/include/dca/linalg/util/gpu_event.hpp @@ -64,7 +64,7 @@ class GpuEvent { }; // Returns the elapsed time in seconds between two recorded events. Blocks host. -float elapsedTime(cudaEvent_t stop, cudaEvent_t start) { +inline float elapsedTime(cudaEvent_t stop, cudaEvent_t start) { checkRC(cudaEventSynchronize(stop)); float msec(0); checkRC(cudaEventElapsedTime(&msec, start, stop)); @@ -91,8 +91,8 @@ class GpuEvent { #endif // DCA_HAVE_GPU -} // util -} // linalg -} // dca +} // namespace util +} // namespace linalg +} // namespace dca #endif // DCA_LINALG_UTIL_GPU_EVENT_HPP diff --git a/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp b/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp new file mode 120000 index 000000000..33c12799f --- /dev/null +++ b/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp @@ -0,0 +1 @@ +epd@a30four.3346297:1765815611 \ No newline at end of file diff --git a/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp new file mode 120000 index 000000000..33c12799f --- /dev/null +++ b/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp @@ -0,0 +1 @@ +epd@a30four.3346297:1765815611 \ No newline at end of file diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.cpp b/include/dca/phys/dca_loop/disorder/apply_disorder.cpp new file mode 100644 index 000000000..9e7411f81 --- /dev/null +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.cpp @@ -0,0 +1,11 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak, Oak Ridge National Lab, (doakpw@ornl.gov) +// +// Here we apply a disorder configuration and disorder potential to a g0 + +#include "apply_disorder.hpp" diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp new file mode 100644 index 000000000..ccb1ec8d8 --- /dev/null +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp @@ -0,0 +1,17 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak, Oak Ridge National Lab, (doakpw@ornl.gov) +// +// Here we apply a disorder configuration and disorder potential to a g0 + +#ifndef DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP +#define DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP +#include + +auto makeDisorderedG0(G0_r_t_clusted_excluded, disorder_configuration, disorder_potential); + +#endif diff --git a/include/dca/phys/domains/cluster/cluster_operations.hpp b/include/dca/phys/domains/cluster/cluster_operations.hpp index f62e84a0f..f4afe0b9f 100644 --- a/include/dca/phys/domains/cluster/cluster_operations.hpp +++ b/include/dca/phys/domains/cluster/cluster_operations.hpp @@ -96,7 +96,8 @@ int cluster_operations::index(const std::vector& element, assert(index > -1 and index < elements.size()); if (math::util::distance2(element, elements[index]) > 1.e-6) { - std::cout << "\n\t " << "cluster_operations::index" << "element mismatch " << "\t" << index << "\n"; + std::cout << "\n\t " << "cluster_operations::index" << "element mismatch " << "\t" << index + << "\n"; math::util::print(element); std::cout << "\n"; math::util::print(elements[index]); @@ -379,8 +380,8 @@ std::vector cluster_operations::find_closest_cluster_vector( return result_vec; } -} // domains -} // phys -} // dca +} // namespace domains +} // namespace phys +} // namespace dca #endif // DCA_PHYS_DOMAINS_CLUSTER_CLUSTER_OPERATIONS_HPP diff --git a/include/dca/phys/domains/quantum/electron_band_domain.hpp b/include/dca/phys/domains/quantum/electron_band_domain.hpp index ce2e37c1a..aa86e2462 100644 --- a/include/dca/phys/domains/quantum/electron_band_domain.hpp +++ b/include/dca/phys/domains/quantum/electron_band_domain.hpp @@ -15,6 +15,8 @@ #include #include +#include +#include "dca/util/to_string.hpp" namespace dca { namespace phys { @@ -45,24 +47,32 @@ class electron_band_domain { return elements_; } + static const auto& get_flavors() { + return flavors_; + } + template static void write(Writer& writer); template static void initialize(const Parameters& parameters); + template + static void print(const electron_band_domain& ebd, ss_type& ss); + // For testing purposes only. - static void setAVectors(const std::vector>& vecs){ - if(vecs.size() != get_size()){ - throw(std::logic_error(__PRETTY_FUNCTION__)); - } - for(int b = 0; b < get_size(); ++b){ - elements_[b].a_vec = vecs[b]; - } + static void setAVectors(const std::vector>& vecs) { + if (vecs.size() != get_size()) { + throw(std::logic_error(__PRETTY_FUNCTION__)); + } + for (int b = 0; b < get_size(); ++b) { + elements_[b].a_vec = vecs[b]; + } } private: static inline std::vector elements_; + static std::vector flavors_; }; template @@ -77,16 +87,22 @@ void electron_band_domain::initialize(const Parameters& /*parameters*/) { elements_.resize(Parameters::bands); using Lattice = typename Parameters::lattice_type; - auto flavours = Lattice::flavors(); + flavors_ = Lattice::flavors(); auto a_vecs = Lattice::aVectors(); for (size_t i = 0; i < a_vecs.size(); ++i) { elements_.at(i).number = i; - elements_.at(i).flavor = flavours.at(i); + elements_.at(i).flavor = flavors_.at(i); elements_.at(i).a_vec = a_vecs.at(i); } } +template +void electron_band_domain::print(const electron_band_domain& ebd, ss_type& ss) { + ss << "\t electron orbitals/bands: " << ebd.get_elements().size() << " of " + << vectorToString(ebd.get_flavors()) << '\n'; +} + } // namespace domains } // namespace phys } // namespace dca diff --git a/include/dca/phys/parameters/main_parameters.hpp b/include/dca/phys/parameters/main_parameters.hpp new file mode 100644 index 000000000..3d676d4d1 --- /dev/null +++ b/include/dca/phys/parameters/main_parameters.hpp @@ -0,0 +1,12 @@ +#ifndef DCA_PHYS_PARAMETERS_MAIN_PARAMETERS_HPP +#define DCA_PHYS_PARAMETERS_MAIN_PARAMETERS_HPP + +#include "dca/config/dca.hpp" + +extern template class dca::phys::params::Parameters< + Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solver_name, + dca::NumericalTraits::type>>; + +#endif diff --git a/src/phys/CMakeLists.txt b/src/phys/CMakeLists.txt index ce58bf35c..b93daefd6 100644 --- a/src/phys/CMakeLists.txt +++ b/src/phys/CMakeLists.txt @@ -3,4 +3,4 @@ add_subdirectory(dca_analysis) add_subdirectory(dca_step) add_subdirectory(domains) add_subdirectory(models) - +add_subdirectory(parameters) diff --git a/src/phys/domains/quantum/CMakeLists.txt b/src/phys/domains/quantum/CMakeLists.txt index 2a1c2bc2d..3be4f7f54 100644 --- a/src/phys/domains/quantum/CMakeLists.txt +++ b/src/phys/domains/quantum/CMakeLists.txt @@ -4,6 +4,7 @@ add_library(quantum_domains STATIC brillouin_zone_path_domain.cpp dca_iteration_domain.cpp electron_spin_domain.cpp + electron_band_domain.cpp numerical_error_domain.cpp point_group_symmetry_element.cpp symmetry_group_level.cpp) diff --git a/src/phys/domains/quantum/electron_band_domain.cpp b/src/phys/domains/quantum/electron_band_domain.cpp new file mode 100644 index 000000000..4b7055d5f --- /dev/null +++ b/src/phys/domains/quantum/electron_band_domain.cpp @@ -0,0 +1,24 @@ +// Copyright (C) 2018 ETH Zurich +// Copyright (C) 2018 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter Staar (taa@zurich.ibm.com) +// +// This file implements electron_spin_domain.hpp. + +#include "dca/phys/domains/quantum/electron_band_domain.hpp" +#include + +namespace dca { +namespace phys { +namespace domains { +// dca::phys::domains:: + +std::vector electron_band_domain::flavors_ = {}; + +} // namespace domains +} // namespace phys +} // namespace dca diff --git a/src/phys/parameters/CMakeLists.txt b/src/phys/parameters/CMakeLists.txt new file mode 100644 index 000000000..bd982f4dc --- /dev/null +++ b/src/phys/parameters/CMakeLists.txt @@ -0,0 +1,4 @@ +add_library(main_parameters STATIC main_parameters.cpp) +target_include_directories(main_parameters PRIVATE ${DCA_INCLUDE_DIRS} ${SIMPLEX_GM_RULE_INCLUDE_DIR}) +# somehow Parameters is dependent on locating fftw3.h! +target_link_libraries(main_parameters PRIVATE FFTW::Double dca_io quantum_domains) diff --git a/src/phys/parameters/main_parameters.cpp b/src/phys/parameters/main_parameters.cpp new file mode 100644 index 000000000..06fb75708 --- /dev/null +++ b/src/phys/parameters/main_parameters.cpp @@ -0,0 +1,9 @@ + + +#include "dca/phys/parameters/main_parameters.hpp" + +template class dca::phys::params::Parameters< + Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solver_name, + dca::NumericalTraits::type>>; diff --git a/test/unit/math/nfft/CMakeLists.txt b/test/unit/math/nfft/CMakeLists.txt index 327923ba6..9ceee3f8f 100644 --- a/test/unit/math/nfft/CMakeLists.txt +++ b/test/unit/math/nfft/CMakeLists.txt @@ -12,4 +12,4 @@ dca_add_gtest(dnfft_1d_gpu_test GTEST_MAIN CUDA INCLUDE_DIRS ${PROJECT_SOURCE_DIR} - LIBS FFTW::Double time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS}) + LIBS FFTW::Double quantum_domains time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS}) diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt index 9e121d50c..1282a6b47 100644 --- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt +++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt @@ -4,22 +4,22 @@ dca_add_gtest(sp_accumulator_gpu_test CUDA GTEST_MAIN INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR} - LIBS FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS} + LIBS FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS} ) dca_add_gtest(sp_accumulator_complex_gpu_test CUDA FAST INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR} - LIBS FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS} + LIBS FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS} ) - + dca_add_gtest(sp_accumulator_single_meas_G_test CUDA GTEST_MAIN INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR} - LIBS FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS} + LIBS FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS} ) #;mc_kernels diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt index df410d593..94d2ac38f 100644 --- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt +++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt @@ -18,4 +18,4 @@ dca_add_gtest(cached_ndft_gpu_test CUDA FAST INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR} - LIBS function_transform ${DCA_KERNEL_LIBS} hdf5::hdf5_cpp magma::magma magma::sparse ${THIS_TEST_LIBS}) + LIBS function_transform quantum_domains ${DCA_KERNEL_LIBS} hdf5::hdf5_cpp magma::magma magma::sparse ${THIS_TEST_LIBS}) diff --git a/test/unit/phys/parameters/disorder_parameters/CMakeLists.txt b/test/unit/phys/parameters/disorder_parameters/CMakeLists.txt new file mode 100644 index 000000000..9f46373dc --- /dev/null +++ b/test/unit/phys/parameters/disorder_parameters/CMakeLists.txt @@ -0,0 +1,6 @@ +# Domains parameters' unit tests + +dca_add_gtest(disorder_parameters_test + FAST + GTEST_MAIN + LIBS json) diff --git a/test/unit/phys/parameters/disorder_parameters/disorder_parameters_test.cpp b/test/unit/phys/parameters/disorder_parameters/disorder_parameters_test.cpp new file mode 100644 index 000000000..d67e5e4c3 --- /dev/null +++ b/test/unit/phys/parameters/disorder_parameters/disorder_parameters_test.cpp @@ -0,0 +1,38 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak, Oak Ridge National Lab, (doakpw@ornl.gov) +// +// This file tests disorder_parameters.hpp +// + +#include "dca/phys/parameters/disorder_parameters.hpp" +#include "gtest/gtest.h" +#include "dca/io/json/json_reader.hpp" + +TEST(DisorderParametersTest, DefaultValues) { + dca::phys::params::DisorderParameters pars; + + EXPECT_EQ(0.0, pars.get_disorder_potential()); + EXPECT_EQ(0.0, pars.get_disorder_density()); + EXPECT_EQ(0, pars.get_disorder_num_configurations()); + EXPECT_EQ(1, pars.get_disorder_max_sites()); +} + +TEST(DomainsParametersTest, ReadAll) { + dca::io::JSONReader reader; + dca::phys::params::DisorderParameters pars; + + reader.open_file(DCA_SOURCE_DIR + "/test/unit/phys/parameters/disorder_parameters/input_read_all.json"); + pars.readWrite(reader); + reader.close_file(); + + EXPECT_EQ(1.0, pars.get_disorder_potential()); + EXPECT_EQ(0.25, pars.get_disorder_density()); + EXPECT_EQ(10, pars.get_disorder_num_configurations()); + EXPECT_EQ(1, pars.get_disorder_max_sites()); +} diff --git a/test/unit/phys/parameters/disorder_parameters/input_read_all.json b/test/unit/phys/parameters/disorder_parameters/input_read_all.json new file mode 100644 index 000000000..0ea68c725 --- /dev/null +++ b/test/unit/phys/parameters/disorder_parameters/input_read_all.json @@ -0,0 +1,8 @@ +{ + "disorder": { + "potential": 1.0, + "density": 0.25, + "num-configurations": 10, + "max-sites": 1 + } +} From 87e10d2f7c8f5a763704ba1cd809b98ba46d306e Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Thu, 5 Mar 2026 15:34:00 -0500 Subject: [PATCH 09/22] initial API for makeDisorder --- include/dca/phys/dca_loop/disorder/apply_disorder.hpp | 4 +++- .../dca => src}/phys/dca_loop/disorder/apply_disorder.cpp | 3 +++ 2 files changed, 6 insertions(+), 1 deletion(-) rename {include/dca => src}/phys/dca_loop/disorder/apply_disorder.cpp (72%) diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp index ccb1ec8d8..b651ccf4e 100644 --- a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp @@ -10,8 +10,10 @@ #ifndef DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP #define DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP +#include "dca/phys/dca_data/dca_data.hpp" #include -auto makeDisorderedG0(G0_r_t_clusted_excluded, disorder_configuration, disorder_potential); +template +auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl); #endif diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.cpp b/src/phys/dca_loop/disorder/apply_disorder.cpp similarity index 72% rename from include/dca/phys/dca_loop/disorder/apply_disorder.cpp rename to src/phys/dca_loop/disorder/apply_disorder.cpp index 9e7411f81..c771c470a 100644 --- a/include/dca/phys/dca_loop/disorder/apply_disorder.cpp +++ b/src/phys/dca_loop/disorder/apply_disorder.cpp @@ -9,3 +9,6 @@ // Here we apply a disorder configuration and disorder potential to a g0 #include "apply_disorder.hpp" + +template +auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl) {} From ebf732175ce69c20c62993b670fd5cdc7a1dcd0b Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Thu, 5 Mar 2026 15:34:23 -0500 Subject: [PATCH 10/22] delete junk files --- include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp | 1 - include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp | 1 - 2 files changed, 2 deletions(-) delete mode 120000 include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp delete mode 120000 include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp diff --git a/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp b/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp deleted file mode 120000 index 33c12799f..000000000 --- a/include/dca/phys/dca_loop/disorder/.#apply_disorder.cpp +++ /dev/null @@ -1 +0,0 @@ -epd@a30four.3346297:1765815611 \ No newline at end of file diff --git a/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp deleted file mode 120000 index 33c12799f..000000000 --- a/include/dca/phys/dca_loop/disorder/.#apply_disorder.hpp +++ /dev/null @@ -1 +0,0 @@ -epd@a30four.3346297:1765815611 \ No newline at end of file From be6a33b874f61d04c2f16102657f1ef8fd7975e8 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Fri, 13 Mar 2026 14:30:19 -0400 Subject: [PATCH 11/22] starting work on inversion --- include/dca/phys/dca_data/dca_data.hpp | 2 ++ include/dca/phys/dca_loop/disorder/apply_disorder.hpp | 9 ++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index 7e3f1d647..b1280f306 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -198,6 +198,8 @@ class DcaData { G0_k_t_cluster_excluded; func::function, func::dmn_variadic> G0_r_w_cluster_excluded; + // Why these are the only G0 tensors that still can be reall or + // complex I'm not clear. func::function> G0_r_t_cluster_excluded; func::function> disordered_G0_r_t_cluster_excluded; diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp index b651ccf4e..d690f348d 100644 --- a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp @@ -14,6 +14,13 @@ #include template -auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl); +auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl) { + int matrix_dim = b::dmn_size() * s::dmn_size(); + + for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) + for (int it = 0; it < TDmn::dmn_size(); ++it) { + dca::linalg::matrixop::copyArrayToMatrix() + } +} #endif From fbbb6816af9945d6dbf1cfb6512ec603227a8b7a Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Mon, 23 Mar 2026 14:58:34 -0400 Subject: [PATCH 12/22] sketching in the loops stage 1 --- include/dca/phys/dca_data/dca_data.hpp | 13 +++++++-- include/dca/phys/dca_loop/dca_loop.hpp | 29 ++++++++++++++++--- .../phys/dca_loop/disorder/apply_disorder.hpp | 24 +++++++++++++-- 3 files changed, 56 insertions(+), 10 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index b1280f306..b40c1d461 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -57,6 +57,7 @@ #include "dca/util/timer.hpp" #include "dca/util/to_string.hpp" #include "dca/distribution/dist_types.hpp" +#include "dca/phys/dca_loop/disorder/apply_disorder.hpp" #ifdef DCA_WITH_ADIOS2 #include "dca/io/adios2/adios2_writer.hpp" #endif @@ -112,6 +113,9 @@ class DcaData { KClusterDmn, WVertexDmn, KExchangeDmn, WExchangeDmn>, DT>; + // Vector NuDmn * RDmn in size that is the disorder configuration + using DisorderConfiguration = func::function>; + DcaData(Parameters& parameters_ref); /** These reads are used by analysis programs only for now. @@ -203,8 +207,6 @@ class DcaData { func::function> G0_r_t_cluster_excluded; func::function> disordered_G0_r_t_cluster_excluded; - // Vector NuDmn * RDmn in size that is the disorder configuration - func::function> disorder_configuration; func::function orbital_occupancy; @@ -597,7 +599,6 @@ void DcaData::initialize_G0() { // Compute G0_r_w. math::transform::FunctionTransform::execute(G0_k_w, G0_r_w); Symmetrize::execute(G0_r_w, H_symmetry, true); - // At this point we can apply the disorder configuration // Compute G0_r_t. math::transform::FunctionTransform::execute(G0_k_t, G0_r_t); @@ -706,6 +707,12 @@ void DcaData::readSigmaFile(io::Reader& reader) { reader.close_group(); } +template +void DcaData::makeDisordedG0(const DisorderConfiguration& disorder_configuration) { + makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded, + disordered_G0_r_t_cluster_excluded); +} + template void DcaData::compute_single_particle_properties() { { diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index 3aff5f15f..f119df0e4 100644 --- a/include/dca/phys/dca_loop/dca_loop.hpp +++ b/include/dca/phys/dca_loop/dca_loop.hpp @@ -247,10 +247,7 @@ void DcaLoop::execute() { // If there is disorder do the disorder application here. // - double L2_Sigma_difference = - solve_cluster_problem(dca_iteration_); // returned from cluster_solver::finalize - - adjust_impurity_self_energy(); // double-counting-correction + workClusters(); perform_lattice_mapping(); @@ -288,6 +285,30 @@ void DcaLoop::execute() { } } +template +void DcaLoop::workClusters() { + if (parameters_.get_disorder_num_configurations() > 0) { + // additional things for summation and getting the post solve G + // to sum will need to be done here + for (int id = 0; id < parameters_.get_disorder_num_configurations()) { + MOMS.makeDisorderedG0(disorder_configurations[id]); + solvingTheCluster(); + } + } + else { + // here we solve just the single ordered cluster + solvingTheCluster(); + } +} + +template +void DcaLoop::solvingTheCluster() { + double L2_Sigma_difference = solve_cluster_problem(dca_iteration_); + // returned from cluster_solver::finalize + + adjust_impurity_self_energy(); // double-counting-correction +} + template void DcaLoop::finalize() { perform_cluster_mapping_self_energy(); diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp index d690f348d..83f9c92b6 100644 --- a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp @@ -10,16 +10,34 @@ #ifndef DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP #define DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP +#include "dca/linalg/matrixop.hpp" #include "dca/phys/dca_data/dca_data.hpp" #include +#include "dca/linalg/linalg.hpp" template -auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl) { - int matrix_dim = b::dmn_size() * s::dmn_size(); +auto makeDisorderedG0( + const decltype(dca::phys::DcaData::disorder_configuration)& disorder_configuration, + const decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl, + decltype(dca::phys::DcaData::disordered_G0_r_t_cluster_excluded)& disordered_g0_r_t_cl_exl) { + int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size(); + dca::linalg::Matrix> g0_rtcex_inverse(matrix_dim); + dca::linalg::Vector ipiv; + dca::linalg::Vector, dca::linalg::CPU> work; for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) for (int it = 0; it < TDmn::dmn_size(); ++it) { - dca::linalg::matrixop::copyArrayToMatrix() + dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, g0_r_t_cl_exl(0, 0, ir, it), + matrix_dim, g0_rtcex_inverse); + dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); + for (int imd = 0; imd < matrix_dim; ++imd) { + g0_rctex_inverse(imd, imd) += disorder_configuration(imd, ir); + } + // Then apply disorder potential to diagonal according to the + // configuration + dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); + dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse, + disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim); } } From eee0b49533d351e1df63ac4cb5ea00c06f3e1793 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Mon, 23 Mar 2026 16:43:04 -0400 Subject: [PATCH 13/22] next stage of getting the loops hooked through discovered much wider access to dca_data g0_r_t than expected --- include/dca/phys/dca_data/dca_data.hpp | 56 ++++++++++++++++--- include/dca/phys/dca_loop/dca_loop.hpp | 27 ++++++--- include/dca/phys/dca_loop/dca_loop_data.hpp | 8 +-- .../phys/dca_loop/disorder/apply_disorder.hpp | 12 ++-- .../tp/tp_equal_time_accumulator.hpp | 8 +-- .../ctaux/ctaux_cluster_solver.hpp | 7 ++- .../g0_interpolation_base.hpp | 16 +++--- .../ctint/walker/ctint_walker_base.hpp | 4 +- 8 files changed, 96 insertions(+), 42 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index b40c1d461..75fda805d 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -57,7 +57,7 @@ #include "dca/util/timer.hpp" #include "dca/util/to_string.hpp" #include "dca/distribution/dist_types.hpp" -#include "dca/phys/dca_loop/disorder/apply_disorder.hpp" +#include "dca/phys/types/dca_shared_types.hpp" #ifdef DCA_WITH_ADIOS2 #include "dca/io/adios2/adios2_writer.hpp" #endif @@ -113,8 +113,11 @@ class DcaData { KClusterDmn, WVertexDmn, KExchangeDmn, WExchangeDmn>, DT>; + using Type_G0_r_t = func::function>; + // Vector NuDmn * RDmn in size that is the disorder configuration - using DisorderConfiguration = func::function>; + using DST = DcaSharedTypes; + using DisorderConfiguration = typename DST::DisorderConfiguration; DcaData(Parameters& parameters_ref); @@ -191,6 +194,11 @@ class DcaData { func::function, func::dmn_variadic> G_r_w; func::function> G_r_t; + func::function, func::dmn_variadic> accumulated_G_k_w; + func::function, func::dmn_variadic> accumulated_G_k_t; + func::function, func::dmn_variadic> accumulated_G_r_w; + func::function> accumulated_G_r_t; + func::function, func::dmn_variadic> G0_k_w; func::function, func::dmn_variadic> G0_k_t; func::function, func::dmn_variadic> G0_r_w; @@ -204,9 +212,9 @@ class DcaData { G0_r_w_cluster_excluded; // Why these are the only G0 tensors that still can be reall or // complex I'm not clear. - func::function> G0_r_t_cluster_excluded; - func::function> - disordered_G0_r_t_cluster_excluded; + /// This G0 could be modified by the diagonal disorder scheme, but + /// normally will just be a copy of the G0_r_t_cluster_excluded + Type_G0_r_t mutable_G0_r_t_cluster_excluded; func::function orbital_occupancy; @@ -270,7 +278,15 @@ class DcaData { return (bool)non_density_interactions_; } + void makeDisordedG0(const DisorderConfiguration& disorder_configuration); + private: // Optional members. + void makeDisordedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& g0_r_t_cl_exl, Type_G0_r_t& disordered_G0_r_t_cl_exl); + + /// Due to the new ability to modify the G0 for disorder this is the + /// immutable g0 from the last iteration. + Type_G0_r_t G0_r_t_cluster_excluded; std::unique_ptr G_k_w_err_; std::unique_ptr G_r_w_err_; std::unique_ptr Sigma_err_; @@ -610,6 +626,8 @@ void DcaData::initialize_G0() { G0_k_t_cluster_excluded = G0_k_t; G0_r_w_cluster_excluded = G0_r_w; G0_r_t_cluster_excluded = G0_r_t; + // When there is no disorder in G0 mutable just gets the unmodified G0_r_t + mutable_G0_r_t_cluster_excluded = G0_r_t; } catch (const std::exception& exc) { std::throw_with_nested(std::runtime_error("Failure in initialization of G0!")); @@ -709,8 +727,7 @@ void DcaData::readSigmaFile(io::Reader& reader) { template void DcaData::makeDisordedG0(const DisorderConfiguration& disorder_configuration) { - makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded, - disordered_G0_r_t_cluster_excluded); + makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded, mutable_G0_r_t_cluster_excluded); } template @@ -900,6 +917,31 @@ void DcaData::print_Sigma_QMC_versus_Sigma_cg() { } } +template +void DcaData::makeDisordedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& g0_r_t_cl_exl, + Type_G0_r_t& disordered_g0_r_t_cl_exl) { + int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size(); + dca::linalg::Matrix, dca::linalg::CPU> g0_rtcex_inverse(matrix_dim); + dca::linalg::Vector ipiv; + dca::linalg::Vector, dca::linalg::CPU> work; + + for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) + for (int it = 0; it < TDmn::dmn_size(); ++it) { + dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, g0_r_t_cl_exl(0, 0, ir, it), + matrix_dim, g0_rtcex_inverse); + dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); + for (int imd = 0; imd < matrix_dim; ++imd) { + g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir); + } + // Then apply disorder potential to diagonal according to the + // configuration + dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); + dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse, + disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim); + } +} + } // namespace phys } // namespace dca diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index f119df0e4..c6cfdd672 100644 --- a/include/dca/phys/dca_loop/dca_loop.hpp +++ b/include/dca/phys/dca_loop/dca_loop.hpp @@ -94,19 +94,26 @@ class DcaLoop { void perform_cluster_exclusion_step(); - double solve_cluster_problem(int DCA_iteration); - void perform_lattice_mapping(); void update_DCA_loop_data_functions(int DCA_iteration); void logSelfEnergy(int i); + /// Temporary name, required to basically curry the disordered G0 + /// loop around the cluster solver. + double workTheClusters(); + /// The shared steps done once normally or side the accumulation + /// loop for disorder + double solvingTheCluster(); + ParametersType& parameters; DcaDataType& MOMS; concurrency_type& concurrency; private: + double solve_cluster_problem(int DCA_iteration); + DcaLoopData DCA_info_struct; cluster_exclusion_type cluster_exclusion_obj; @@ -247,7 +254,7 @@ void DcaLoop::execute() { // If there is disorder do the disorder application here. // - workClusters(); + auto L2_Sigma_difference = workTheClusters(); perform_lattice_mapping(); @@ -286,27 +293,29 @@ void DcaLoop::execute() { } template -void DcaLoop::workClusters() { - if (parameters_.get_disorder_num_configurations() > 0) { +double DcaLoop::workTheClusters() { + if (parameters.get_disorder_num_configurations() > 0) { // additional things for summation and getting the post solve G // to sum will need to be done here - for (int id = 0; id < parameters_.get_disorder_num_configurations()) { - MOMS.makeDisorderedG0(disorder_configurations[id]); + auto num_configurations = parameters.get_disorder_num_configurations(); + for (int id = 0; id < num_configurations; ++id) { + MOMS.makeDisorderedG0(DCA_info_struct.disorder_configurations[id]); solvingTheCluster(); } } else { // here we solve just the single ordered cluster - solvingTheCluster(); + return solvingTheCluster(); } } template -void DcaLoop::solvingTheCluster() { +double DcaLoop::solvingTheCluster() { double L2_Sigma_difference = solve_cluster_problem(dca_iteration_); // returned from cluster_solver::finalize adjust_impurity_self_energy(); // double-counting-correction + return L2_Sigma_difference; } template diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp index 69035b6ff..e23267e1d 100644 --- a/include/dca/phys/dca_loop/dca_loop_data.hpp +++ b/include/dca/phys/dca_loop/dca_loop_data.hpp @@ -23,6 +23,7 @@ #include "dca/phys/domains/quantum/dca_iteration_domain.hpp" #include "dca/phys/domains/quantum/electron_band_domain.hpp" #include "dca/phys/domains/quantum/electron_spin_domain.hpp" +#include "dca/phys/types/dca_shared_types.hpp" #ifdef DCA_HAVE_ADIOS2 #include "dca/io/adios2/adios2_writer.hpp" #endif @@ -52,6 +53,8 @@ class DcaLoopData { func::dmn_0>; + using DST = DcaSharedTypes; + using DisorderConfiguration = typename DST::DisorderConfiguration; DcaLoopData(); template @@ -92,10 +95,7 @@ class DcaLoopData { func::function average_defect_density; func::function num_defect_configurations; - // What should the data structure for the defect configurations look - // I think the same as the r_DCA domain with just a +/- 1 likely a - // Real to make it easy multiply with the disorder potential - + std::vector disorder_configurations; int last_completed_iteration = -1; }; diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp index 83f9c92b6..959849e31 100644 --- a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp +++ b/include/dca/phys/dca_loop/disorder/apply_disorder.hpp @@ -12,14 +12,14 @@ #define DCA_PHYS_DCA_LOOP_DISORDER_APPLY_DISORDER_HPP #include "dca/linalg/matrixop.hpp" #include "dca/phys/dca_data/dca_data.hpp" -#include #include "dca/linalg/linalg.hpp" -template -auto makeDisorderedG0( - const decltype(dca::phys::DcaData::disorder_configuration)& disorder_configuration, - const decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl, - decltype(dca::phys::DcaData::disordered_G0_r_t_cluster_excluded)& disordered_g0_r_t_cl_exl) { +template +void DcaData::makeDisorderedG0( + const DisorderConfiguration disorder_configuration, + const decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl, + decltype(dca::phys::DcaData::disordered_G0_r_t_cluster_excluded)& + disordered_g0_r_t_cl_exl) { int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size(); dca::linalg::Matrix> g0_rtcex_inverse(matrix_dim); dca::linalg::Vector ipiv; diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp index 535582b29..c9102e584 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp @@ -41,7 +41,7 @@ namespace solver { namespace ctaux { // dca::phys::solver::ctaux:: using dca::util::SignType; - + template class TpEqualTimeAccumulator { public: @@ -253,7 +253,7 @@ TpEqualTimeAccumulator::TpEqualTimeAccumulator(const Parameter G_r_t("G_r_t"), G_r_t_stddev("G_r_t_stddev"), - + G_r_t_accumulated("G_r_t_accumulated"), G_r_t_accumulated_squared("G_r_t_accumulated_squared"), @@ -339,7 +339,7 @@ void TpEqualTimeAccumulator::initialize_akima_coefficients() { for (int nu1_ind = 0; nu1_ind < nu::dmn_size(); nu1_ind++) { for (int nu0_ind = 0; nu0_ind < nu::dmn_size(); nu0_ind++) { for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++) - y[t_ind] = MOMS.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + y[t_ind] = MOMS.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); ai_obj.initialize(x, y); @@ -357,7 +357,7 @@ void TpEqualTimeAccumulator::initialize_akima_coefficients() { for (int nu0_ind = 0; nu0_ind < nu::dmn_size(); nu0_ind++) { for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size(); t_ind++) y[t_ind - t::dmn_size() / 2] = - MOMS.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + MOMS.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); ai_obj.initialize(x, y); diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp index b6edbe850..d089248db 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp @@ -117,7 +117,9 @@ class CtauxClusterSolver { return g0_; }; - typename Walker::Resource& getResource() { return dummy_walker_resource_; }; + typename Walker::Resource& getResource() { + return dummy_walker_resource_; + }; protected: void warmUp(Walker& walker); @@ -164,6 +166,7 @@ class CtauxClusterSolver { G0Interpolation g0_; typename Walker::Resource dummy_walker_resource_; + private: Rng rng_; @@ -230,7 +233,7 @@ template ::initialize(int dca_iteration) { dca_iteration_ = dca_iteration; - g0_.initializeShrinked(data_.G0_r_t_cluster_excluded); + g0_.initializeShrinked(data_.mutable_G0_r_t_cluster_excluded); Sigma_old_ = data_.Sigma; diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp index 6d6ee0d8a..6b3d990d1 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp @@ -127,10 +127,10 @@ void G0InterpolationBase::initialize_linear_coefficients(Data& data) for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) { for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) { G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind) = - data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); grad_G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind) = - (data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) - - data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind)); + (data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) - + data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind)); } } } @@ -141,10 +141,10 @@ void G0InterpolationBase::initialize_linear_coefficients(Data& data) for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) { for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) { G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind - 1) = - data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); grad_G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind - 1) = - (data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) - - data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind)); + (data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) - + data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind)); } } } @@ -169,7 +169,7 @@ void G0InterpolationBase::initialize_akima_coefficients(Data& data) for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) { for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) { for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++) - y[t_ind] = data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + y[t_ind] = data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); ai_obj.initialize(x, y); @@ -187,7 +187,7 @@ void G0InterpolationBase::initialize_akima_coefficients(Data& data) for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) { for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size(); t_ind++) y[t_ind - t::dmn_size() / 2] = - data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); + data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind); ai_obj.initialize(x, y); diff --git a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp index 469785075..295bfe138 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp @@ -41,7 +41,7 @@ #include "dca/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu.hpp" #endif -//#define DEBUG_SUBMATRIX +// #define DEBUG_SUBMATRIX namespace dca { namespace phys { @@ -330,7 +330,7 @@ void CtintWalkerBase::setInteractionVertices(const Data& data, vertices_.initialize(parameters.getDoubleUpdateProbability(), parameters.getAllSitesPartnership()); vertices_.initializeFromHamiltonian(data.H_interactions); if (data.has_non_density_interactions()) { - vertices_.checkForInterbandPropagators(data.G0_r_t_cluster_excluded); + vertices_.checkForInterbandPropagators(data.mutable_G0_r_t_cluster_excluded); vertices_.initializeFromNonDensityHamiltonian(data.get_non_density_interactions()); } } From 41a59a960a632fdcbe2f02843d3eb4661e1e8089 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 31 Mar 2026 11:27:26 -0400 Subject: [PATCH 14/22] moving makeDisorderedG0 defintion --- .../disorder/apply_disorder.hpp => dca_data/make_disorder.hpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename include/dca/phys/{dca_loop/disorder/apply_disorder.hpp => dca_data/make_disorder.hpp} (100%) diff --git a/include/dca/phys/dca_loop/disorder/apply_disorder.hpp b/include/dca/phys/dca_data/make_disorder.hpp similarity index 100% rename from include/dca/phys/dca_loop/disorder/apply_disorder.hpp rename to include/dca/phys/dca_data/make_disorder.hpp From 49fb54275c2a6f051758fff0a73ffe531763a17e Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 31 Mar 2026 13:51:13 -0400 Subject: [PATCH 15/22] actually compiling for square --- .../dca/linalg/util/atomic_add_cuda.cu.hpp | 39 +++------- include/dca/phys/dca_data/dca_data.hpp | 73 +++++++++++-------- 2 files changed, 51 insertions(+), 61 deletions(-) diff --git a/include/dca/linalg/util/atomic_add_cuda.cu.hpp b/include/dca/linalg/util/atomic_add_cuda.cu.hpp index a0afb4b61..d6875f636 100644 --- a/include/dca/linalg/util/atomic_add_cuda.cu.hpp +++ b/include/dca/linalg/util/atomic_add_cuda.cu.hpp @@ -21,28 +21,7 @@ namespace dca { namespace linalg { // dca::linalg:: - -#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600 -// Older devices do not have an hardware atomicAdd for double. -// See -// https://stackoverflow.com/questions/12626096/why-has-atomicadd-not-been-implemented-for-doubles -__device__ double inline atomicAddImpl(double* address, const double val) { - unsigned long long int* address_as_ull = (unsigned long long int*)address; - unsigned long long int old = *address_as_ull, assumed; - do { - assumed = old; - old = atomicCAS(address_as_ull, assumed, - __double_as_longlong(val + __longlong_as_double(assumed))); - // Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN) } - } while (assumed != old); - return __longlong_as_double(old); -} - -__device__ void inline atomicAdd(double* address, const double val) { - atomicAddImpl(address, val); -} - -#elif defined(DCA_HAVE_HIP) +#if defined(DCA_HAVE_HIP) // HIP seems to have some horrible problem with concurrent atomic operations. __device__ double inline atomicAddImpl(double* address, const double val) { unsigned long long int* address_as_ull = (unsigned long long int*)address; @@ -61,13 +40,12 @@ __device__ double inline atomicAddImpl(float* address, const float val) { unsigned long int old = *address_as_int, assumed; do { assumed = old; - old = atomicCAS(address_as_int, assumed, - __float_as_int(val + __int_as_float(assumed))); + old = atomicCAS(address_as_int, assumed, __float_as_int(val + __int_as_float(assumed))); // Note: uses integer comparison to avoid hang in case of NaN (since NaN != NaN) } } while (assumed != old); return __int_as_float(old); } - + __device__ void inline atomicAdd(float* address, const float val) { atomicAddImpl(address, val); } @@ -82,7 +60,7 @@ __device__ void inline atomicAdd(cuDoubleComplex* address, cuDoubleComplex val) atomicAddImpl(a_d + 1, val.y); } - __device__ void inline atomicAdd(magmaFloatComplex* const address, magmaFloatComplex val) { +__device__ void inline atomicAdd(magmaFloatComplex* const address, magmaFloatComplex val) { double* a_d = reinterpret_cast(address); atomicAddImpl(a_d, val.x); atomicAddImpl(a_d + 1, val.y); @@ -105,12 +83,13 @@ __device__ void inline atomicAdd(float* address, float val) { __device__ void inline atomicAdd(cuDoubleComplex* address, cuDoubleComplex val) { double* a_d = reinterpret_cast(address); - atomicAdd(a_d, val.x); - atomicAdd(a_d + 1, val.y); + ::atomicAdd(a_d, val.x); + ::atomicAdd(a_d + 1, val.y); } + #endif // atomic operation help -} // linalg -} // dca +} // namespace linalg +} // namespace dca #endif // DCA_LINALG_UTIL_ATOMIC_ADD_CUDA_CU_HPP diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index 75fda805d..a1eb7c50d 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -212,6 +213,8 @@ class DcaData { G0_r_w_cluster_excluded; // Why these are the only G0 tensors that still can be reall or // complex I'm not clear. + Type_G0_r_t G0_r_t_cluster_excluded; + /// This G0 could be modified by the diagonal disorder scheme, but /// normally will just be a copy of the G0_r_t_cluster_excluded Type_G0_r_t mutable_G0_r_t_cluster_excluded; @@ -278,15 +281,14 @@ class DcaData { return (bool)non_density_interactions_; } - void makeDisordedG0(const DisorderConfiguration& disorder_configuration); + void makeDisorderedG0(const DisorderConfiguration& disorder_configuration); private: // Optional members. - void makeDisordedG0(const DisorderConfiguration& disorder_configuration, - const Type_G0_r_t& g0_r_t_cl_exl, Type_G0_r_t& disordered_G0_r_t_cl_exl); + void makeDisorderedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& g0_r_t_cl_exl, Type_G0_r_t& disordered_G0_r_t_cl_exl); /// Due to the new ability to modify the G0 for disorder this is the /// immutable g0 from the last iteration. - Type_G0_r_t G0_r_t_cluster_excluded; std::unique_ptr G_k_w_err_; std::unique_ptr G_r_w_err_; std::unique_ptr Sigma_err_; @@ -592,7 +594,6 @@ void DcaData::initializeH0_and_H_i() { } Parameters::model_type::initialize_H_symmetries(H_symmetry); - compute_band_structure::execute(parameters_, band_structure); } @@ -726,7 +727,42 @@ void DcaData::readSigmaFile(io::Reader& reader) { } template -void DcaData::makeDisordedG0(const DisorderConfiguration& disorder_configuration) { +void DcaData::makeDisorderedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& g0_r_t_cl_exl, + Type_G0_r_t& disordered_g0_r_t_cl_exl) { + int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size(); + dca::linalg::Matrix> + g0_rtcex_inverse(matrix_dim); + // dca::linalg::Vector, dca::linalg::CPU>> + // ipiv; + // dca::linalg::Vector, dca::linalg::CPU, + // dca::linalg::util::DefaultAllocator, dca::linalg::CPU>> + // work; + + for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) + for (int it = 0; it < TDmn::dmn_size(); ++it) { + dca::linalg::Vector> + rt_block(std::size_t(NuDmn::dmn_size() * NuDmn::dmn_size())); + g0_r_t_cl_exl.slice(0, 1, {0, 0, ir, it}, rt_block.data()); + dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, rt_block.data(), matrix_dim, + g0_rtcex_inverse); + dca::linalg::matrixop::inverse(g0_rtcex_inverse); //, ipiv, work); + for (int imd = 0; imd < matrix_dim; ++imd) { + g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir); + } + // Then apply disorder potential to diagonal according to the + // configuration + dca::linalg::matrixop::inverse(g0_rtcex_inverse); + dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse, + &disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim); + } +} + +template +void DcaData::makeDisorderedG0(const DisorderConfiguration& disorder_configuration) { makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded, mutable_G0_r_t_cluster_excluded); } @@ -917,31 +953,6 @@ void DcaData::print_Sigma_QMC_versus_Sigma_cg() { } } -template -void DcaData::makeDisordedG0(const DisorderConfiguration& disorder_configuration, - const Type_G0_r_t& g0_r_t_cl_exl, - Type_G0_r_t& disordered_g0_r_t_cl_exl) { - int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size(); - dca::linalg::Matrix, dca::linalg::CPU> g0_rtcex_inverse(matrix_dim); - dca::linalg::Vector ipiv; - dca::linalg::Vector, dca::linalg::CPU> work; - - for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) - for (int it = 0; it < TDmn::dmn_size(); ++it) { - dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, g0_r_t_cl_exl(0, 0, ir, it), - matrix_dim, g0_rtcex_inverse); - dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); - for (int imd = 0; imd < matrix_dim; ++imd) { - g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir); - } - // Then apply disorder potential to diagonal according to the - // configuration - dca::linalg::matrixop::inverse(g0_rtcex_inverse, ipiv, work); - dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse, - disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim); - } -} - } // namespace phys } // namespace dca From ac0f742434538c8876141cd13a56ec8098c44edc Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Wed, 1 Apr 2026 18:05:18 -0400 Subject: [PATCH 16/22] I think this is a complete initial disorder impl --- include/dca/phys/dca_data/dca_data.hpp | 25 +++++-- include/dca/phys/dca_loop/dca_loop.hpp | 68 +++++++++++------- include/dca/phys/dca_loop/dca_loop_data.hpp | 3 +- .../ctaux/ctaux_cluster_solver.hpp | 70 +++++++++++++++++-- .../stdthread_qmci_cluster_solver.hpp | 27 +++++-- 5 files changed, 150 insertions(+), 43 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index a1eb7c50d..7744f51fd 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -114,6 +114,10 @@ class DcaData { KClusterDmn, WVertexDmn, KExchangeDmn, WExchangeDmn>, DT>; + using Type_G0_k_w = + func::function, func::dmn_variadic>; + using Type_G0_k_t = + func::function, func::dmn_variadic>; using Type_G0_r_t = func::function>; // Vector NuDmn * RDmn in size that is the disorder configuration @@ -205,10 +209,8 @@ class DcaData { func::function, func::dmn_variadic> G0_r_w; func::function> G0_r_t; - func::function, func::dmn_variadic> - G0_k_w_cluster_excluded; - func::function, func::dmn_variadic> - G0_k_t_cluster_excluded; + Type_G0_k_w G0_k_w_cluster_excluded; + Type_G0_k_t G0_k_t_cluster_excluded; func::function, func::dmn_variadic> G0_r_w_cluster_excluded; // Why these are the only G0 tensors that still can be reall or @@ -218,6 +220,8 @@ class DcaData { /// This G0 could be modified by the diagonal disorder scheme, but /// normally will just be a copy of the G0_r_t_cluster_excluded Type_G0_r_t mutable_G0_r_t_cluster_excluded; + Type_G0_k_w mutable_G0_k_w_cluster_excluded; + Type_G0_k_t mutable_G0_k_t_cluster_excluded; func::function orbital_occupancy; @@ -629,6 +633,7 @@ void DcaData::initialize_G0() { G0_r_t_cluster_excluded = G0_r_t; // When there is no disorder in G0 mutable just gets the unmodified G0_r_t mutable_G0_r_t_cluster_excluded = G0_r_t; + mutable_G0_k_w_cluster_excluded = G0_k_w; } catch (const std::exception& exc) { std::throw_with_nested(std::runtime_error("Failure in initialization of G0!")); @@ -750,15 +755,23 @@ void DcaData::makeDisorderedG0(const DisorderConfiguration& diso dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, rt_block.data(), matrix_dim, g0_rtcex_inverse); dca::linalg::matrixop::inverse(g0_rtcex_inverse); //, ipiv, work); + + // Then apply disorder potential to diagonal according to the + // configuration for (int imd = 0; imd < matrix_dim; ++imd) { g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir); } - // Then apply disorder potential to diagonal according to the - // configuration + dca::linalg::matrixop::inverse(g0_rtcex_inverse); dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse, &disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim); } + + // Giving up at this point in side effect free method. + math::transform::FunctionTransform::execute( + disordered_g0_r_t_cl_exl, mutable_G0_k_t_cluster_excluded); + // math::transform::FunctionTransform::execute(mutable_G0_k_t_cluster_excluded, + // mutable_G0_k_w_cluster_excluded); } template diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index c6cfdd672..f68cdba59 100644 --- a/include/dca/phys/dca_loop/dca_loop.hpp +++ b/include/dca/phys/dca_loop/dca_loop.hpp @@ -41,6 +41,8 @@ #include "dca/util/print_time.hpp" #include "dca/util/signal_handler.hpp" +#include "dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp" + namespace dca { namespace phys { // dca::phys:: @@ -103,16 +105,16 @@ class DcaLoop { /// Temporary name, required to basically curry the disordered G0 /// loop around the cluster solver. double workTheClusters(); - /// The shared steps done once normally or side the accumulation - /// loop for disorder - double solvingTheCluster(); ParametersType& parameters; DcaDataType& MOMS; concurrency_type& concurrency; private: - double solve_cluster_problem(int DCA_iteration); + void solve_cluster_problem(int DCA_iteration); + double finalize_cluster_problem(); + void accumulateGkw(); + void averageGkw(); DcaLoopData DCA_info_struct; @@ -256,6 +258,8 @@ void DcaLoop::execute() { auto L2_Sigma_difference = workTheClusters(); + adjust_impurity_self_energy(); // double-counting-correction + perform_lattice_mapping(); update_DCA_loop_data_functions(dca_iteration_); // Really just updates @@ -297,25 +301,22 @@ double DcaLoop::workTheClus if (parameters.get_disorder_num_configurations() > 0) { // additional things for summation and getting the post solve G // to sum will need to be done here - auto num_configurations = parameters.get_disorder_num_configurations(); + MakeDisorderConfigurations make_disorder_configurations; + make_disorder_configurations(parameters, DCA_info_struct.disorder_configurations); + int num_configurations = DCA_info_struct.disorder_configurations.size(); for (int id = 0; id < num_configurations; ++id) { MOMS.makeDisorderedG0(DCA_info_struct.disorder_configurations[id]); - solvingTheCluster(); + solve_cluster_problem(dca_iteration_); } + averageGkw(); } else { // here we solve just the single ordered cluster - return solvingTheCluster(); + solve_cluster_problem(dca_iteration_); + monte_carlo_integrator_.collectSingle(); } -} - -template -double DcaLoop::solvingTheCluster() { - double L2_Sigma_difference = solve_cluster_problem(dca_iteration_); - // returned from cluster_solver::finalize - - adjust_impurity_self_energy(); // double-counting-correction - return L2_Sigma_difference; + auto L2_Sigma_Difference = finalize_cluster_problem(); + return L2_Sigma_Difference; } template @@ -398,7 +399,7 @@ void DcaLoop::perform_clust } template -double DcaLoop::solve_cluster_problem( +void DcaLoop::solve_cluster_problem( int DCA_iteration) { // static_assert(std::is_same>::value); // static_assert(std::is_same>::value); @@ -414,18 +415,33 @@ double DcaLoop::solve_clust if (output_file_ && output_file_->isADIOS2()) output_file_->flush(); +} - { - if (concurrency.id() == concurrency.first()) - std::cout << "start Monte Carlo integration finalize.\n"; +template +void DcaLoop::accumulateGkw() { + monte_carlo_integrator_.accumulateGkw(); +} - profiler_type profiler("finalize cluster-solver", "DCA", __LINE__); - double L2_Sigma_difference = monte_carlo_integrator_.finalize(DCA_info_struct); - if (concurrency.id() == concurrency.first()) - std::cout << "Monte Carlo integration finalized.\n"; +template +void DcaLoop::averageGkw() { + MOMS.G_k_w = MOMS.accumulated_G_k_w; + MOMS.G_k_w /= parameters.get_disorder_num_configurations(); +} - return L2_Sigma_difference; - } +template +double DcaLoop::finalize_cluster_problem() { + if (concurrency.id() == concurrency.first()) + std::cout << "start Monte Carlo integration finalize.\n"; + + profiler_type profiler("finalize cluster-solver", "DCA", __LINE__); + + // So what we do here varies if disorder is activated. + + double L2_Sigma_difference = monte_carlo_integrator_.finalize(DCA_info_struct); + if (concurrency.id() == concurrency.first()) + std::cout << "Monte Carlo integration finalized.\n"; + + return L2_Sigma_difference; } template diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp index e23267e1d..297bee17f 100644 --- a/include/dca/phys/dca_loop/dca_loop_data.hpp +++ b/include/dca/phys/dca_loop/dca_loop_data.hpp @@ -1,11 +1,12 @@ // Copyright (C) 2018 ETH Zurich -// Copyright (C) 2018 UT-Battelle, LLC +// Copyright (C) 2026 UT-Battelle, LLC // All rights reserved. // // See LICENSE for terms of usage. // See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. // // Author: Peter Staar (taa@zurich.ibm.com) +// Peter W. Doak (doakpw@ornl.gov) // // This class contains physical and computational data of the DCA(+) loop. diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp index d089248db..ed695153b 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp @@ -29,6 +29,7 @@ #include "dca/function/domains.hpp" #include "dca/function/function.hpp" #include "dca/linalg/linalg.hpp" +#include "dca/linalg/util/allocators/allocators.hpp" #include "dca/math/function_transform/function_transform.hpp" #include "dca/math/statistics/util.hpp" #include "dca/parallel/util/get_workload.hpp" @@ -38,6 +39,7 @@ #include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/time_correlator.hpp" #include "dca/phys/dca_step/symmetrization/symmetrize.hpp" #include "dca/phys/domains/cluster/cluster_domain.hpp" +#include "dca/phys/domains/domain_aliases.hpp" #include "dca/phys/domains/quantum/electron_band_domain.hpp" #include "dca/phys/domains/quantum/electron_spin_domain.hpp" #include "dca/phys/domains/time_and_frequency/frequency_domain.hpp" @@ -75,6 +77,9 @@ class CtauxClusterSolver { static constexpr linalg::DeviceType device = device_t; protected: + // The following are deprecated because they don't follow the + // naming rules nor are they consistent with other classes aliases + // for these domains. using w = func::dmn_0; using b = func::dmn_0; using s = func::dmn_0; @@ -83,7 +88,9 @@ class CtauxClusterSolver { using CDA = ClusterDomainAliases; using RDmn = typename CDA::RClusterDmn; using KDmn = typename CDA::KClusterDmn; - + using DDA = DcaDomainAliases; + using WDmn = typename DDA::WDmn; + using NuDmn = typename DDA::NuDmn; using NuNuKClusterWDmn = func::dmn_variadic; using NuNuRClusterWDmn = func::dmn_variadic; @@ -121,6 +128,9 @@ class CtauxClusterSolver { return dummy_walker_resource_; }; + void accumulateGkw(); + void collectSingle(); + protected: void warmUp(Walker& walker); @@ -135,6 +145,8 @@ class CtauxClusterSolver { void collect_measurements(); void compute_G_k_w_from_M_r_w(); + void accumulateGkwFromMrw( + func::function, func::dmn_variadic>& G_k_w); double compute_S_k_w_from_G_k_w(); @@ -292,15 +304,24 @@ void CtauxClusterSolver::integrate() { } template -template -double CtauxClusterSolver::finalize( - dca_info_struct_t& dca_info_struct) { +void CtauxClusterSolver::accumulateGkw() { + collect_measurements(); + accumulate_G_k_w_from_M_r_w(data_.G_k_w); +} + +template +void CtauxClusterSolver::collectSingle() { collect_measurements(); symmetrize_measurements(); // Compute new Sigma. compute_G_k_w_from_M_r_w(); +} +template +template +double CtauxClusterSolver::finalize( + dca_info_struct_t& dca_info_struct) { // FT::execute(data_.G_k_w, data_.G_r_w); math::transform::FunctionTransform::execute(data_.G_k_w, data_.G_r_w); @@ -639,6 +660,47 @@ void CtauxClusterSolver::compute_G_k_w_from_M_ Symmetrize::execute(data_.G_k_w, data_.H_symmetry); } +template +void CtauxClusterSolver::accumulateGkwFromMrw( + func::function, func::dmn_variadic>& G_k_w) { + func::function, NuNuKClusterWDmn> M_k_w; + math::transform::FunctionTransform::execute(M_r_w_, M_k_w); + + const std::size_t matrix_size = b::dmn_size() * s::dmn_size(); + linalg::Matrix, dca::linalg::CPU> G0_times_M_matrix(matrix_size); + using MatrixView = + linalg::MatrixView, dca::linalg::CPU, + linalg::util::DefaultAllocator, dca::linalg::CPU>>; + + // G = G0 - G0*M*G0/beta + for (int k_ind = 0; k_ind < KDmn::dmn_size(); k_ind++) { + for (int w_ind = 0; w_ind < w::dmn_size(); w_ind++) { + // These views make strong assumptions about the function layouts! + const MatrixView G0_matrix(&data_.mutable_G0_k_w_cluster_excluded(0, 0, 0, 0, k_ind, w_ind), + matrix_size); + const MatrixView M_matrix(&M_k_w(0, 0, 0, 0, k_ind, w_ind), matrix_size); + + // G0 * M --> G0_times_M_matrix + linalg::matrixop::gemm(G0_matrix, M_matrix, G0_times_M_matrix); + + MatrixView G_matrix(&G_k_w(0, 0, 0, 0, k_ind, w_ind), matrix_size); + + // G0_times_M_matrix * G0 --> G_matrix + linalg::matrixop::gemm(G0_times_M_matrix, G0_matrix, G_matrix); + + // -G_matrix / beta + G0_cluster_excluded_matrix --> G_matrix + for (int j = 0; j < matrix_size; ++j) + for (int i = 0; i < matrix_size; ++i) + G_matrix(i, j) = -G_matrix(i, j) / parameters_.get_beta() + G0_matrix(i, j); + } + } + data_.accumulated_G_k_w += G_k_w; + // This seems pretty dicey to me I think with the disorder some + // symmetry is only guaranteed if enough disorder configurations are + // summed over. + // Symmetrize::execute(data_.G_k_w, data_.H_symmetry); +} + template double CtauxClusterSolver::compute_S_k_w_from_G_k_w() { static double alpha = parameters_.get_self_energy_mixing_factor(); diff --git a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp index 676a72559..c87463926 100644 --- a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp @@ -116,6 +116,11 @@ class StdThreadQmciClusterSolver : public QmciSolver { void logSingleMeasurement(StdThreadAccumulatorType& accumulator, int stamping_period, bool log_MFunction, bool log_MFunctionTime) const; + /// This collects measurements when dealing with disordered G0 + void accumulateGkw(); + /// This collects measurements when dealing with single G0 + void collectSingle(); + private: void startWalker(int id); void startAccumulator(int id, const Parameters& parameters); @@ -182,7 +187,7 @@ StdThreadQmciClusterSolver::StdThreadQmciClusterSolver( accumulators_queue_(), config_dump_(nr_walkers_) - //autocorrelation_data_(parameters_, 0, BaseClass::g0_) +// autocorrelation_data_(parameters_, 0, BaseClass::g0_) { if (nr_walkers_ < 1 || nr_accumulators_ < 1) { throw std::logic_error( @@ -367,8 +372,8 @@ void StdThreadQmciClusterSolver::startWalker(int id) { const int walker_index = thread_task_handler_.walkerIDToRngIndex(id); auto walker_log = last_iteration_ ? BaseClass::writer_ : nullptr; - Walker walker(parameters_, data_, rng_vector_[walker_index], BaseClass::getResource(), concurrency_.get_id(), id, - walker_log, BaseClass::g0_); + Walker walker(parameters_, data_, rng_vector_[walker_index], BaseClass::getResource(), + concurrency_.get_id(), id, walker_log, BaseClass::g0_); std::unique_ptr exception_ptr; @@ -499,6 +504,16 @@ auto StdThreadQmciClusterSolver::computeSingleMeasurement_G_k_w( return G_k_w; } +template +void StdThreadQmciClusterSolver::accumulateGkw() { + QmciSolver::accumulateGkw(); +} + +template +void StdThreadQmciClusterSolver::collectSingle() { + QmciSolver::collectSingle(); +} + template void StdThreadQmciClusterSolver::logSingleMeasurement( StdThreadAccumulatorType& accumulator_obj, int stamping_period, bool log_MFunction, @@ -589,8 +604,8 @@ void StdThreadQmciClusterSolver::startWalkerAndAccumulator(int id, // Create and warm a walker. auto walker_log = BaseClass::writer_; - Walker walker(parameters_, data_, rng_vector_[id], BaseClass::getResource(), concurrency_.get_id(), id, walker_log, - BaseClass::g0_); + Walker walker(parameters_, data_, rng_vector_[id], BaseClass::getResource(), + concurrency_.get_id(), id, walker_log, BaseClass::g0_); initializeAndWarmUp(walker, id, id); if (id == 0) { @@ -641,7 +656,7 @@ void StdThreadQmciClusterSolver::startWalkerAndAccumulator(int id, catch (...) { throw std::runtime_error("something mysterious went wrong in walker thread!"); } - + ++walk_finished_; if (BaseClass::writer_ && BaseClass::writer_->isADIOS2()) BaseClass::writer_->flush(); From ed84dab0eb63c1bbba2611572c8d6fdb1bb025cc Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Mon, 6 Apr 2026 19:13:52 -0400 Subject: [PATCH 17/22] added configuration weighting and manual t -> w ft --- include/dca/phys/dca_data/dca_data.hpp | 24 +++++++++++++++++-- include/dca/phys/dca_loop/dca_loop.hpp | 16 +++++++++---- include/dca/phys/dca_loop/dca_loop_data.hpp | 1 + .../ctaux/ctaux_cluster_solver.hpp | 14 ++++++----- .../stdthread_qmci_cluster_solver.hpp | 6 ++--- 5 files changed, 45 insertions(+), 16 deletions(-) diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index 7744f51fd..064147263 100644 --- a/include/dca/phys/dca_data/dca_data.hpp +++ b/include/dca/phys/dca_data/dca_data.hpp @@ -59,6 +59,7 @@ #include "dca/util/to_string.hpp" #include "dca/distribution/dist_types.hpp" #include "dca/phys/types/dca_shared_types.hpp" +#include "dca/util/type_help.hpp" #ifdef DCA_WITH_ADIOS2 #include "dca/io/adios2/adios2_writer.hpp" #endif @@ -770,8 +771,27 @@ void DcaData::makeDisorderedG0(const DisorderConfiguration& diso // Giving up at this point in side effect free method. math::transform::FunctionTransform::execute( disordered_g0_r_t_cl_exl, mutable_G0_k_t_cluster_excluded); - // math::transform::FunctionTransform::execute(mutable_G0_k_t_cluster_excluded, - // mutable_G0_k_w_cluster_excluded); + + // I can't figure out how to get this to work with the + // math::transform::FunctionTransform framwork so do this by hand + // here. + using Complex = std::complex>; + for (int k = 0; k < KClusterDmn::dmn_size(); ++k) { + const auto& k_val = KClusterDmn::get_elements()[k]; + for (int inu1 = 0; inu1 < NuDmn::dmn_size(); ++inu1) + for (int inu2 = 0; inu2 < NuDmn::dmn_size(); ++inu2) { + for (int w = 0; w < WDmn::dmn_size(); ++w) { + const auto& w_val = WDmn::get_elements()[w]; + Complex G_k_omega(0); + for (int t = 0; t < TDmn::dmn_size(); ++t) { + const auto& t_val = TDmn::get_elements()[t]; + G_k_omega += mutable_G0_k_t_cluster_excluded(inu1, inu2, k, t) * + std::exp(Complex(0, w_val * t_val)); + } + mutable_G0_k_w_cluster_excluded(inu1, inu2, k, w) = G_k_omega; + } + } + } } template diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index f68cdba59..ec0c1e26a 100644 --- a/include/dca/phys/dca_loop/dca_loop.hpp +++ b/include/dca/phys/dca_loop/dca_loop.hpp @@ -113,7 +113,7 @@ class DcaLoop { private: void solve_cluster_problem(int DCA_iteration); double finalize_cluster_problem(); - void accumulateGkw(); + void accumulateGkw(double weight); void averageGkw(); DcaLoopData DCA_info_struct; @@ -302,10 +302,12 @@ double DcaLoop::workTheClus // additional things for summation and getting the post solve G // to sum will need to be done here MakeDisorderConfigurations make_disorder_configurations; - make_disorder_configurations(parameters, DCA_info_struct.disorder_configurations); + make_disorder_configurations(parameters, DCA_info_struct.disorder_configurations, + DCA_info_struct.disorder_weights); int num_configurations = DCA_info_struct.disorder_configurations.size(); for (int id = 0; id < num_configurations; ++id) { MOMS.makeDisorderedG0(DCA_info_struct.disorder_configurations[id]); + std::cout << "Solving disorder configuration " << id << '\n'; solve_cluster_problem(dca_iteration_); } averageGkw(); @@ -418,14 +420,18 @@ void DcaLoop::solve_cluster } template -void DcaLoop::accumulateGkw() { - monte_carlo_integrator_.accumulateGkw(); +void DcaLoop::accumulateGkw(double weight) { + monte_carlo_integrator_.accumulateGkw(weight); } template void DcaLoop::averageGkw() { MOMS.G_k_w = MOMS.accumulated_G_k_w; - MOMS.G_k_w /= parameters.get_disorder_num_configurations(); + auto& disorder_weights = DCA_info_struct.disorder_weights; + auto total_disorder_weight = std::accumulate(disorder_weights.begin(), disorder_weights.end(), 0.0); + MOMS.G_k_w /= total_disorder_weight; + std::cout << " Averaged over " << disorder_weights.size() + << " disorder configurations with weight " << total_disorder_weight << '\n'; } template diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp index 297bee17f..113a598dd 100644 --- a/include/dca/phys/dca_loop/dca_loop_data.hpp +++ b/include/dca/phys/dca_loop/dca_loop_data.hpp @@ -97,6 +97,7 @@ class DcaLoopData { func::function average_defect_density; func::function num_defect_configurations; std::vector disorder_configurations; + std::vector disorder_weights; int last_completed_iteration = -1; }; diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp index ed695153b..73b3e2e6f 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp @@ -128,7 +128,7 @@ class CtauxClusterSolver { return dummy_walker_resource_; }; - void accumulateGkw(); + void accumulateGkw(double weight); void collectSingle(); protected: @@ -146,7 +146,8 @@ class CtauxClusterSolver { void compute_G_k_w_from_M_r_w(); void accumulateGkwFromMrw( - func::function, func::dmn_variadic>& G_k_w); + func::function, func::dmn_variadic>& G_k_w, + double weight); double compute_S_k_w_from_G_k_w(); @@ -304,9 +305,9 @@ void CtauxClusterSolver::integrate() { } template -void CtauxClusterSolver::accumulateGkw() { +void CtauxClusterSolver::accumulateGkw(double weight) { collect_measurements(); - accumulate_G_k_w_from_M_r_w(data_.G_k_w); + accumulateGkwfromMrw(data_.G_k_w); } template @@ -662,7 +663,8 @@ void CtauxClusterSolver::compute_G_k_w_from_M_ template void CtauxClusterSolver::accumulateGkwFromMrw( - func::function, func::dmn_variadic>& G_k_w) { + func::function, func::dmn_variadic>& G_k_w, + double weight) { func::function, NuNuKClusterWDmn> M_k_w; math::transform::FunctionTransform::execute(M_r_w_, M_k_w); @@ -694,7 +696,7 @@ void CtauxClusterSolver::accumulateGkwFromMrw( G_matrix(i, j) = -G_matrix(i, j) / parameters_.get_beta() + G0_matrix(i, j); } } - data_.accumulated_G_k_w += G_k_w; + data_.accumulated_G_k_w += G_k_w * weight; // This seems pretty dicey to me I think with the disorder some // symmetry is only guaranteed if enough disorder configurations are // summed over. diff --git a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp index c87463926..2c4a40343 100644 --- a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp @@ -117,7 +117,7 @@ class StdThreadQmciClusterSolver : public QmciSolver { bool log_MFunction, bool log_MFunctionTime) const; /// This collects measurements when dealing with disordered G0 - void accumulateGkw(); + void accumulateGkw(double weight); /// This collects measurements when dealing with single G0 void collectSingle(); @@ -505,8 +505,8 @@ auto StdThreadQmciClusterSolver::computeSingleMeasurement_G_k_w( } template -void StdThreadQmciClusterSolver::accumulateGkw() { - QmciSolver::accumulateGkw(); +void StdThreadQmciClusterSolver::accumulateGkw(double weight) { + QmciSolver::accumulateGkw(weight); } template From 64df292cf6364e9550a40741a6d807206c698f3d Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 7 Apr 2026 15:23:24 -0400 Subject: [PATCH 18/22] fix issue with writing "sample" conf with disorder --- .../stdthread_qmci_cluster_solver.hpp | 34 +++++++++++++------ 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp index 2c4a40343..0ef835a1c 100644 --- a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp +++ b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp @@ -128,8 +128,16 @@ class StdThreadQmciClusterSolver : public QmciSolver { void initializeAndWarmUp(Walker& walker, int id, int walker_id); - void readConfigurations(); - void writeConfigurations() const; + /** readConfigurations needs to be aware that there are different + * vertex configurations for difference disorder configurations + * but this is perhaps not very useful now since we generate new + * disorder configurations for each iteration. + */ + void readConfigurations(int disorder_configuiration = -1); + /** writeConfigurations needs to be aware that there are different + * vertex configurations for difference disorder configurations + */ + void writeConfigurations(int disorder_configuration = -1) const; void iterateOverLocalMeasurements(int walker_id, std::function&& f); @@ -266,7 +274,6 @@ void StdThreadQmciClusterSolver::integrate() { dca::profiling::Duration duration(end_time, start_time); total_time_ = duration.sec + 1.e-6 * duration.usec; - printIntegrationMetadata(); }; @@ -281,7 +288,7 @@ void StdThreadQmciClusterSolver::integrate() { print_metadata(); - if (parameters_.store_configuration()) { + if (parameters_.store_configuration() && parameters_.get_disorder_num_configurations() <= 0) { if (BaseClass::writer_) { if (BaseClass::writer_->isADIOS2()) { BaseClass::writer_->open_group("Configurations"); @@ -352,8 +359,8 @@ double StdThreadQmciClusterSolver::finalize(dca_info_struct_t& dca_i } // Write and reset autocorrelation. - std::cout << "Writing autocorrelation data\n"; - std::cout << "Autocorrelation incompatible with complex G0 and GPU"; + // std::cout << "Writing autocorrelation data\n"; + // std::cout << "Autocorrelation incompatible with complex G0 and GPU"; // autocorrelation_data_.write(*BaseClass::writer_, dca_iteration_); } // autocorrelation_data_.reset(); @@ -693,13 +700,16 @@ void StdThreadQmciClusterSolver::finalizeWalker(Walker& walker, int } template -void StdThreadQmciClusterSolver::writeConfigurations() const { +void StdThreadQmciClusterSolver::writeConfigurations(int disorder_configuration) const { if (parameters_.get_directory_config_write() == "") return; try { + std::string config_label; + if (disorder_configuration >= 0) + config_label = "_disorder_config_" + std::to_string(disorder_configuration); const std::string out_name = parameters_.get_directory_config_write() + "/process_" + - std::to_string(concurrency_.id()) + ".hdf5"; + std::to_string(concurrency_.id()) + config_label + ".hdf5"; io::HDF5Writer writer(false); writer.open_file(out_name); for (int id = 0; id < config_dump_.size(); ++id) @@ -711,15 +721,19 @@ void StdThreadQmciClusterSolver::writeConfigurations() const { } template -void StdThreadQmciClusterSolver::readConfigurations() { +void StdThreadQmciClusterSolver::readConfigurations(int disorder_configuration) { if (parameters_.get_directory_config_read() == "") return; Profiler profiler(__FUNCTION__, "stdthread-MC", __LINE__); try { + std::string config_label; + if (disorder_configuration >= 0) + config_label = "_disorder_config_" + std::to_string(disorder_configuration); + const std::string inp_name = parameters_.get_directory_config_read() + "/process_" + - std::to_string(concurrency_.id()) + ".hdf5"; + std::to_string(concurrency_.id()) + config_label + ".hdf5"; io::HDF5Reader reader(false); reader.open_file(inp_name); for (int id = 0; id < config_dump_.size(); ++id) From e3663f8f84be428a61085187d2fe02d5d8ffc60e Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Tue, 26 May 2026 17:00:45 -0400 Subject: [PATCH 19/22] lost main_parameters lib in rebase --- applications/dca/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/dca/CMakeLists.txt b/applications/dca/CMakeLists.txt index 94bbb9d39..b6cb96558 100644 --- a/applications/dca/CMakeLists.txt +++ b/applications/dca/CMakeLists.txt @@ -8,7 +8,7 @@ if (DCA_BUILD_DCA) target_link_libraries(main_dca PRIVATE ${DCA_KERNEL_LIBS}) endif() - target_link_libraries(main_dca PUBLIC FFTW::Double signals ${DCA_LIBS} dca_io) + target_link_libraries(main_dca PUBLIC FFTW::Double signals ${DCA_LIBS} dca_io main_parameters) install(TARGETS main_dca RUNTIME DESTINATION bin) From f8e5c7b1537e2d23a7f7064ef1a4f8dd3394ec38 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Wed, 27 May 2026 13:07:14 -0400 Subject: [PATCH 20/22] missing file --- include/dca/phys/types/dca_shared_types.hpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 include/dca/phys/types/dca_shared_types.hpp diff --git a/include/dca/phys/types/dca_shared_types.hpp b/include/dca/phys/types/dca_shared_types.hpp new file mode 100644 index 000000000..c8d4e182d --- /dev/null +++ b/include/dca/phys/types/dca_shared_types.hpp @@ -0,0 +1,18 @@ +#ifndef DCA_PHYS_TYPES_SHARED_TYPES_HPP +#define DCA_PHYS_TYPES_SHARED_TYPES_HPP + +#include "dca/phys/domains/cluster/cluster_domain_aliases.hpp" +#include "dca/phys/domains/domain_aliases.hpp" +#include "dca/function/function.hpp" +#include "dca/function/domains/dmn_variadic.hpp" + +namespace dca::phys { +template +struct DcaSharedTypes { + using DDA = DcaDomainAliases; + using Real = typename DDA::Real; + using DisorderConfiguration = + func::function>; +}; +} // namespace dca::phys +#endif From 0115061b2905f542f395613804d75aeb427527cd Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Wed, 27 May 2026 13:12:07 -0400 Subject: [PATCH 21/22] another missing file --- include/dca/phys/domains/domain_aliases.hpp | 34 +++++++++++++++++++++ 1 file changed, 34 insertions(+) create mode 100644 include/dca/phys/domains/domain_aliases.hpp diff --git a/include/dca/phys/domains/domain_aliases.hpp b/include/dca/phys/domains/domain_aliases.hpp new file mode 100644 index 000000000..7d4097ab1 --- /dev/null +++ b/include/dca/phys/domains/domain_aliases.hpp @@ -0,0 +1,34 @@ +#ifndef DCA_PHYS_DOMAINS_DOMAIN_ALIASES_HPP +#define DCA_PHYS_DOMAINS_DOMAIN_ALIASES_HPP + +#include "dca/phys/domains/cluster/cluster_domain_aliases.hpp" + +namespace dca::phys { +template +struct DcaDomainAliases { + using Real = typename Parameters::Real; + using Scalar = typename Parameters::Scalar; + using TpAccumulatorPrec = typename Parameters::TPAccumPrec; + using TpComplex = std::complex; + using TDmn = func::dmn_0; + using WDmn = func::dmn_0; + using WVertexDmn = func::dmn_0>; + using WExchangeDmn = func::dmn_0; + + using BDmn = func::dmn_0; + using SDmn = func::dmn_0; + using NuDmn = func::dmn_variadic; // orbital-spin index + using NuNuDmn = func::dmn_variadic; + static constexpr int Dimension = Parameters::lattice_type::DIMENSION; + using CDA = ClusterDomainAliases; + using RClusterDmn = typename CDA::RClusterDmn; + using KClusterType = typename CDA::KClusterType; + using KClusterDmn = typename CDA::KClusterDmn; + using RHostDmn = typename CDA::RSpHostDmn; + using KHostDmn = typename CDA::KSpHostDmn; + using KExchangeDmn = func::dmn_0; +}; + +} // namespace dca::phys + +#endif From 7c8e8b905126b45c0e749350c527c48393458e64 Mon Sep 17 00:00:00 2001 From: Peter Doak Date: Wed, 27 May 2026 13:13:21 -0400 Subject: [PATCH 22/22] functor to make disorder configurations --- .../disorder/makeDisorderConfigurations.hpp | 138 ++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp diff --git a/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp b/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp new file mode 100644 index 000000000..826f43f6a --- /dev/null +++ b/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp @@ -0,0 +1,138 @@ +// Copyright (C) 2026 UT-Battelle, LLC +// All rights reserved. +// +// See LICENSE for terms of usage. +// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications. +// +// Author: Peter W. Doak (doakpw@ornl.gov) +// +// This function creates a set of disorder configurations based on the +// sites in the cluster and requested density + +#ifndef DCA_PHYS_DCA_LOOP_MAKE_DISORDER_HPP +#define DCA_PHYS_DCA_LOOP_MAKE_DISORDER_HPP + +#include +#include +#include +#include "dca/phys/types/dca_shared_types.hpp" + +namespace dca::phys { + +namespace detail { +struct VectorHash { + size_t operator()(const std::vector& v) const { + size_t seed = 0; + for (int x : v) { + // XOR current hash with hash of element + constant + bit rotations + // This helps avoid collisions for similar vectors (e.g., + // [1,-1] vs [-1,1]) + // Based on the boost hash_combine function. + constexpr size_t golden_ratio = sizeof(int) == 8 ? 0x9e3779b97f4a7c15 : 0x9e3779b9; + seed ^= std::hash{}(x) + golden_ratio + (seed << 6) + (seed >> 2); + } + return seed; + } +}; + +} // namespace detail + +template +class MakeDisorderConfigurations { + using DST = DcaSharedTypes; + using DisorderConfiguration = typename DST::DisorderConfiguration; + using DDA = typename DST::DDA; + using RDmn = typename DDA::RClusterDmn; + using BDmn = typename DDA::BDmn; + +public: + void operator()(Parameters& parameters, std::vector& disorder_configurations, + std::vector& disorder_weights) { + auto m_configs = parameters.get_disorder_num_configurations(); + assert(m_configs > 0); + disorder_configurations.resize(m_configs); + disorder_weights.resize(m_configs); + + auto n_bands = BDmn::dmn_size(); + auto n_rsites = RDmn::dmn_size(); + auto n_sites = disorder_configurations[0].size(); + + // right now we always make spin synmmetric configurations + n_sites /= 2; + + // Check feasibility: only 2^n distinct ±1 vectors of length n exist + if (m_configs > (1 << n_sites)) { + throw std::runtime_error("More unique vectors requested than possible."); + } + + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_int_distribution<> dist(0, 1); + + // Set to track uniqueness using our custom hash + std::unordered_set, detail::VectorHash> seen; + + // here we put the disorder potential in while preserving spin + // symmetry for the potential + auto insert_config = [n_rsites, n_bands](const auto& this_config, auto& config) { + int i_disorder = 0; + for (int isite = 0; isite < n_rsites; ++isite) { + for (int iband = 0; iband < n_bands; ++iband, ++i_disorder) { + config(iband, 0, isite) = this_config[i_disorder]; + config(iband, 1, isite) = this_config[i_disorder]; + } + } + }; + + int max_ones = parameters.get_disorder_max_sites(); + int configs_generated = 0; + int attempts = 0; + int max_attempts = 1000; + // Generate vectors until we have m unique ones + while (configs_generated < m_configs && attempts < max_attempts) { + auto& config = disorder_configurations[configs_generated]; + auto& weight = disorder_weights[configs_generated]; + std::vector this_config(n_sites); + int num_dis = 0; + for (auto& site : this_config) { + // Map 0 → -1 and 1 → 1 + site = dist(gen) ? 1 : -1; + if (site == 1) + ++num_dis; + if (num_dis > max_ones) + break; + } + + if (num_dis <= max_ones) { + // Insert returns true only if this vector is new (unique) + if (seen.insert(this_config).second) { + ++configs_generated; + attempts = 0; + insert_config(this_config, config); + auto conc = parameters.get_disorder_density(); + weight = std::pow(conc, num_dis) * std::pow((1 - conc), (n_sites - num_dis)); + + // If disorder didn't have spin symmetry + // std::copy(this_config.begin(), this_config.end(), config.begin()); + } + else { + ++attempts; + } + } + else { + ++attempts; + } + } + auto disorder_half_pot = parameters.get_disorder_potential() * 0.5; + for (int ic = 0; ic < disorder_configurations.size(); ++ic) { + auto& config = disorder_configurations[ic]; + auto weight = disorder_weights[ic]; + config *= disorder_half_pot; + config.print_elements(std::cout); + std::cout << "weight: " << weight << '\n'; + } + } +}; +} // namespace dca::phys + +#endif