diff --git a/applications/dca/CMakeLists.txt b/applications/dca/CMakeLists.txt index f7a9c3762..b6cb96558 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) + + target_link_libraries(main_dca PUBLIC FFTW::Double signals ${DCA_LIBS} dca_io main_parameters) 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/function/function.hpp b/include/dca/function/function.hpp index 359402f1e..5b40279c9 100644 --- a/include/dca/function/function.hpp +++ b/include/dca/function/function.hpp @@ -408,6 +408,21 @@ class function { // These are the linear start and end indexes with respect to the complete function. std::size_t start_; std::size_t end_; + +public: +#ifdef DEBUG + // Variadic debug accessor: forwards all args to operator() + template + const scalartype& value_at_debug(Ts&&... indices) const { + return (*this)(std::forward(indices)...); + } + + // Non-const version (if needed for writing) + template + scalartype& value_at_debug(Ts&&... indices) { + return (*this)(std::forward(indices)...); + } +#endif }; template 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/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/linalg/util/util_lapack.hpp b/include/dca/linalg/util/util_lapack.hpp index 3a29639be..80df3679b 100644 --- a/include/dca/linalg/util/util_lapack.hpp +++ b/include/dca/linalg/util/util_lapack.hpp @@ -60,7 +60,7 @@ inline void checkLapackInfoInternal(int info, std::string function_name, std::st #define warnLapackInfo(info) \ dca::linalg::lapack::util::warnLapackInfoInternal(info, __FUNCTION__, __FILE__, __LINE__) inline void warnLapackInfoInternal(int info, std::string function_name, std::string file_name, - int line) { + int line) { if (info < 0) { std::stringstream s; s << "Error in function: " << function_name << " (" << file_name << ":" << line << ")" @@ -69,15 +69,16 @@ inline void warnLapackInfoInternal(int info, std::string function_name, std::str throw LapackException(s.str(), info); } - else if (info > 0) { - std::cout << "warning lapack info = " << info << " at " << file_name << ":" << line << '\n'; - } +#ifndef NDEBUG + // else if (info > 0) { + // std::cout << "warning lapack info = " << info << " at " << file_name << ":" << line << '\n'; + // } +#endif } - -} // util -} // lapack -} // linalg -} // dca +} // namespace util +} // namespace lapack +} // namespace linalg +} // namespace dca #endif // DCA_LINALG_UTIL_UTIL_LAPACK_HPP diff --git a/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp b/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp index ef86b651c..052a56ae3 100644 --- a/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp +++ b/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp @@ -22,8 +22,7 @@ namespace transform { // dca::math::transform:: // Empty template declaration. -template +template struct hermite_cubic_spline {}; // Template specialization for the equidistant periodic case. @@ -117,8 +116,98 @@ typename hermite_cubic_spline +class hermite_cubic_spline { +private: + typedef typename lh_dmn_type::dmn_specifications_type lh_spec_dmn_type; + typedef typename rh_dmn_type::dmn_specifications_type rh_spec_dmn_type; + + typedef typename lh_spec_dmn_type::scalar_type lh_scalar_type; + typedef typename rh_spec_dmn_type::scalar_type rh_scalar_type; + + typedef typename lh_spec_dmn_type::element_type lh_element_type; + typedef typename rh_spec_dmn_type::element_type rh_element_type; + + typedef lh_scalar_type f_scalar_type; + +public: + static f_scalar_type execute(int i, int j); +}; + +template +typename hermite_cubic_spline::f_scalar_type hermite_cubic_spline< + lh_dmn_type, rh_dmn_type, EQUIDISTANT, INTERVAL, DIMENSION>::execute(int i, int j) { + const static rh_scalar_type a = -0.5; + + lh_element_type x = lh_dmn_type::get_elements()[i]; + rh_element_type y = rh_dmn_type::get_elements()[j]; + + rh_scalar_type* super_basis = rh_dmn_type::get_super_basis(); + + rh_scalar_type* inv_basis = rh_dmn_type::get_inverse_basis(); + rh_scalar_type* inv_super_basis = rh_dmn_type::get_inverse_super_basis(); + + rh_element_type delta(DIMENSION, 0.); + rh_element_type delta_affine(DIMENSION, 0.); + + f_scalar_type result = 0; + + { + for (int li = 0; li < DIMENSION; li++) + delta[li] = (y[li] - x[li]); + + { + for (int li = 0; li < DIMENSION; li++) + delta_affine[li] = 0.; + + for (int li = 0; li < DIMENSION; li++) + for (int lj = 0; lj < DIMENSION; lj++) + delta_affine[li] += inv_super_basis[li + lj * DIMENSION] * delta[lj]; + } + + for (int li = 0; li < DIMENSION; li++) { + while (delta_affine[li] > 0.5 - 1.e-6) + delta_affine[li] -= 1.; + + while (delta_affine[li] < -0.5 - 1.e-6) + delta_affine[li] += 1.; + } + + { + for (int li = 0; li < DIMENSION; li++) + delta[li] = 0.; + + for (int li = 0; li < DIMENSION; li++) + for (int lj = 0; lj < DIMENSION; lj++) + delta[li] += super_basis[li + lj * DIMENSION] * delta_affine[lj]; + } + + { + for (int li = 0; li < DIMENSION; li++) + delta_affine[li] = 0.; + + for (int li = 0; li < DIMENSION; li++) + for (int lj = 0; lj < DIMENSION; lj++) + delta_affine[li] += inv_basis[li + lj * DIMENSION] * delta[lj]; + } + + for (int li = 0; li < DIMENSION; li++) + delta[li] = (delta_affine[li] > -2. and delta_affine[li] < 2.) + ? hermite_spline::cubic(0., delta_affine[li], 1., a) + : 0.; + + f_scalar_type t_result = 1; + for (int li = 0; li < DIMENSION; li++) + t_result *= delta[li]; + + result += t_result; + } + + return result; +} + +} // namespace transform +} // namespace math +} // namespace dca #endif // DCA_MATH_FUNCTION_TRANSFORM_HERMITE_SPLINES_HERMITE_CUBIC_SPLINE_HPP diff --git a/include/dca/math/util/vector_operations.hpp b/include/dca/math/util/vector_operations.hpp index cd8bf7923..638622b3f 100644 --- a/include/dca/math/util/vector_operations.hpp +++ b/include/dca/math/util/vector_operations.hpp @@ -104,16 +104,49 @@ std::complex innerProduct(const std::vector>& x, return res; } +template +std::complex innerProduct(const std::vector& x, const std::vector>& y) { + assert(x.size() == y.size()); + + std::complex res(0.); + for (std::size_t i = 0; i < x.size(); ++i) + res += x[i] * std::conj(y[i]); + + return res; +} + +template +std::complex innerProduct(const std::vector>& x, const std::vector& y) { + assert(x.size() == y.size()); + + std::complex res(0.); + for (std::size_t i = 0; i < x.size(); ++i) + res += x[i] * y[i]; + + return res; +} + // Treats scalars as vectors of size 1. template T innerProduct(const T x, const T y) { return x * y; } + template std::complex innerProduct(const std::complex x, const std::complex y) { return x * std::conj(y); } +// template +// std::complex innerProduct(const T x, const std::complex y) { +// return x * std::conj(y); +// } + +// template +// std::complex innerProduct(const std::complex x, const T y) { +// return x * y; +// } + // Computes the square of the L^2 norm of the vector x. template auto l2Norm2(const std::vector& x) { diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp index d1fafc22b..efc7a2fba 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. @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -28,6 +29,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" @@ -56,6 +58,8 @@ #include "dca/util/timer.hpp" #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 @@ -102,8 +106,12 @@ class DcaData { using SpGreensFunction = func::function, func::dmn_variadic>; + using SpDisorderedGreensFunction = + func::function, func::dmn_variadic>; using SpRGreensFunction = func::function, func::dmn_variadic>; + using SpRDisorderedGreensFunction = + func::function, func::dmn_variadic>; using TpGreensFunction = func::function, 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>; + + using Disordered_G0_k_k_w = + func::function, func::dmn_variadic>; + using Disordered_G0_k_k_t = + func::function, func::dmn_variadic>; + using Disordered_G0_r_r_t = + func::function>; + + // Vector NuDmn * RDmn in size that is the disorder configuration + using DST = DcaSharedTypes; + using DisorderConfiguration = typename DST::DisorderConfiguration; + DcaData(Parameters& parameters_ref); /** These reads are used by analysis programs only for now. @@ -125,6 +150,9 @@ class DcaData { template void write(Writer& writer); + template + void writeDisorderConfiguration(Writer& writer); + #ifdef DCA_WITH_ADIOS2 void writeDistributedG4Adios(io::ADIOS2Writer& writer); #endif @@ -186,18 +214,28 @@ 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; 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; - func::function> G0_r_t_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; + + /// These G0 are only used when diagonal disorder is applied + Disordered_G0_r_r_t disordered_G0_r_r_t_cl_exl; + Disordered_G0_k_k_w disordered_G0_k_k_w_cl_exl; + Disordered_G0_k_k_t disordered_G0_k_k_t_cl_exl; func::function orbital_occupancy; @@ -261,7 +299,14 @@ class DcaData { return (bool)non_density_interactions_; } + void makeDisorderedG0(const DisorderConfiguration& disorder_configuration); + private: // Optional members. + void makeDisorderedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& 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. std::unique_ptr G_k_w_err_; std::unique_ptr G_r_w_err_; std::unique_ptr Sigma_err_; @@ -365,19 +410,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 - } + } } } @@ -533,6 +578,14 @@ void DcaData::write(Writer& writer) { writer.close_group(); } +template +template +void DcaData::writeDisorderConfiguration(Writer& writer) { + if (parameters_.dump_cluster_Greens_functions()) { + writer.execute(G_k_w); + } +} + template void DcaData::initialize() { initializeH0_and_H_i(); @@ -555,7 +608,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.")); } } @@ -566,7 +620,6 @@ void DcaData::initializeH0_and_H_i() { } Parameters::model_type::initialize_H_symmetries(H_symmetry); - compute_band_structure::execute(parameters_, band_structure); } @@ -659,13 +712,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"; } @@ -697,6 +750,84 @@ void DcaData::readSigmaFile(io::Reader& reader) { reader.close_group(); } +template +void DcaData::makeDisorderedG0(const DisorderConfiguration& disorder_configuration, + const Type_G0_r_t& g0_r_t_cl_exl) { + int nu_matrix_dim = NuDmn::dmn_size(); + dca::linalg::Matrix> + g0_rtcex_inverse(nu_matrix_dim); + + int nu_r_matrix_dim = NuDmn::dmn_size() * RClusterDmn::dmn_size(); + dca::linalg::Matrix> + g0_dis_tcex_inverse(nu_matrix_dim); + // dca::linalg::Vector, dca::linalg::CPU>> + // ipiv; + // dca::linalg::Vector, dca::linalg::CPU, + // dca::linalg::util::DefaultAllocator, dca::linalg::CPU>> + // work; + + dca::linalg::Vector> + rt_block(std::size_t(NuDmn::dmn_size() * NuDmn::dmn_size())); + + for (int it = 0; it < TDmn::dmn_size(); ++it) { + for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir) { + g0_r_t_cl_exl.slice(0, 1, {0, 0, ir, it}, rt_block.data()); + dca::linalg::matrixop::copyArrayToMatrix(nu_matrix_dim, nu_matrix_dim, rt_block.data(), + nu_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 < nu_matrix_dim; ++imd) { + g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir); + } + dca::linalg::matrixop::inverse(g0_rtcex_inverse); + for (int imd = 0; imd < nu_matrix_dim; ++imd) { + disordered_G0_r_r_t_cl_exl(imd, ir, imd, ir, it) = g0_rtcex_inverse(imd, imd); + } + } + } + + math::transform::FunctionTransform::execute(disordered_G0_r_r_t_cl_exl, + disordered_G0_k_k_t_cl_exl); + // At this point I beleive the matrix is nondiagonal in k + + // disordered_G0_k_k_t_cl_exl.print_elements(std::cout); + + // 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 k1 = 0; k1 < KClusterDmn::dmn_size(); ++k1) { + const auto& k1_val = KClusterDmn::get_elements()[k1]; + for (int k2 = 0; k2 < KClusterDmn::dmn_size(); ++k2) { + const auto& k2_val = KClusterDmn::get_elements()[k2]; + 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 += disordered_G0_k_k_t_cl_exl(inu1, k1, inu2, k2, t) * + std::exp(Complex(0, w_val * t_val)); + } + disordered_G0_k_k_w_cl_exl(inu1, k1, inu2, k2, w) = G_k_omega; + } + } + } + } +} + +template +void DcaData::makeDisorderedG0(const DisorderConfiguration& disorder_configuration) { + makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded); +} + template void DcaData::compute_single_particle_properties() { { diff --git a/include/dca/phys/dca_data/make_disorder.hpp b/include/dca/phys/dca_data/make_disorder.hpp new file mode 100644 index 000000000..959849e31 --- /dev/null +++ b/include/dca/phys/dca_data/make_disorder.hpp @@ -0,0 +1,44 @@ +// 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 "dca/linalg/matrixop.hpp" +#include "dca/phys/dca_data/dca_data.hpp" +#include "dca/linalg/linalg.hpp" + +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; + 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_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); + } +} + +#endif diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp index bef3c0658..ec0c1e26a 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:: @@ -94,19 +96,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(); + ParametersType& parameters; DcaDataType& MOMS; concurrency_type& concurrency; private: + void solve_cluster_problem(int DCA_iteration); + double finalize_cluster_problem(); + void accumulateGkw(double weight); + void averageGkw(); + DcaLoopData DCA_info_struct; cluster_exclusion_type cluster_exclusion_obj; @@ -240,8 +249,14 @@ void DcaLoop::execute() { perform_cluster_exclusion_step(); - double L2_Sigma_difference = - solve_cluster_problem(dca_iteration_); // returned from cluster_solver::finalize + // 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. + // + + auto L2_Sigma_difference = workTheClusters(); adjust_impurity_self_energy(); // double-counting-correction @@ -281,6 +296,31 @@ void DcaLoop::execute() { } } +template +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 + MakeDisorderConfigurations make_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(); + } + else { + // here we solve just the single ordered cluster + solve_cluster_problem(dca_iteration_); + monte_carlo_integrator_.collectSingle(); + } + auto L2_Sigma_Difference = finalize_cluster_problem(); + return L2_Sigma_Difference; +} + template void DcaLoop::finalize() { perform_cluster_mapping_self_energy(); @@ -361,7 +401,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); @@ -377,18 +417,37 @@ 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(double weight) { + monte_carlo_integrator_.accumulateGkw(weight); +} - 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; + 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'; +} - 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 76e128f23..113a598dd 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. @@ -18,10 +19,12 @@ #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" #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 @@ -47,6 +50,12 @@ class DcaLoopData { using k_DCA = func::dmn_0>; + using r_DCA = + func::dmn_0>; + + using DST = DcaSharedTypes; + using DisorderConfiguration = typename DST::DisorderConfiguration; DcaLoopData(); template @@ -85,6 +94,10 @@ class DcaLoopData { func::function chemical_potential; func::function average_expansion_order; + func::function average_defect_density; + func::function num_defect_configurations; + std::vector disorder_configurations; + std::vector disorder_weights; int last_completed_iteration = -1; }; @@ -113,7 +126,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 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/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 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..4b1eb9296 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"), diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp index 343b3f37e..906af5d38 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp @@ -70,7 +70,7 @@ class CtauxAccumulator : public MC_accumulator_data using ParametersType = Parameters; using DataType = Data; using BaseClass = MC_accumulator_data; - + typedef vertex_pair vertex_pair_type; typedef vertex_singleton vertex_singleton_type; @@ -255,8 +255,8 @@ class CtauxAccumulator : public MC_accumulator_data * in the initializers. */ template -CtauxAccumulator::CtauxAccumulator( - const Parameters& parameters_ref, Data& data_ref, int id) +CtauxAccumulator::CtauxAccumulator(const Parameters& parameters_ref, + Data& data_ref, int id) : MC_accumulator_data(), parameters_(parameters_ref), @@ -307,9 +307,8 @@ void CtauxAccumulator::initialize(int dca_iter parameters_.dump_every_iteration()); if (perform_equal_time_accumulation_) { - equal_time_accumulator_ptr_ = - std::make_unique>(parameters_, data_, - thread_id); + equal_time_accumulator_ptr_ = std::make_unique>( + parameters_, data_, thread_id); equal_time_accumulator_ptr_->resetAccumulation(); } } @@ -406,8 +405,8 @@ void CtauxAccumulator::updateFrom(walker_type& walker.get_error_distribution() = 0; #endif // DCA_WITH_QMC_BIT - //single_particle_accumulator_obj.syncStreams(*event); - //two_particle_accumulator_.syncStreams(*event); + // single_particle_accumulator_obj.syncStreams(*event); + // two_particle_accumulator_.syncStreams(*event); } template @@ -519,7 +518,8 @@ void CtauxAccumulator::accumulate_equal_time_q template void CtauxAccumulator::accumulate_two_particle_quantities() { profiler_type profiler("tp-accumulation", "CT-AUX accumulator", __LINE__, thread_id); - gflop_ += 1e-9 * two_particle_accumulator_.accumulate(M_, hs_configuration_, current_phase_.getSign()); + gflop_ += + 1e-9 * two_particle_accumulator_.accumulate(M_, hs_configuration_, current_phase_.getSign()); } template 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..11eb77237 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" @@ -71,10 +73,15 @@ class CtauxClusterSolver { using Walker = ctaux::CtauxWalker; using Accumulator = ctaux::CtauxAccumulator; using SpGreensFunction = typename Data::SpGreensFunction; + using SpDisorderedGreensFunction = typename Data::SpDisorderedGreensFunction; + using SpRDisorderedGreensFunction = typename Data::SpRDisorderedGreensFunction; 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 +90,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; @@ -95,6 +104,7 @@ class CtauxClusterSolver { void write(Writer& writer); void initialize(int dca_iteration); + void initialize(int dca_iteration, SpRDisorderedGreensFunction& g0_disordered); void integrate(); @@ -117,7 +127,12 @@ class CtauxClusterSolver { return g0_; }; - typename Walker::Resource& getResource() { return dummy_walker_resource_; }; + typename Walker::Resource& getResource() { + return dummy_walker_resource_; + }; + + void accumulateGkw(double weight); + void collectSingle(); protected: void warmUp(Walker& walker); @@ -133,6 +148,10 @@ class CtauxClusterSolver { void collect_measurements(); void compute_G_k_w_from_M_r_w(); + // Now this needs to have k1, k2 + void accumulateGkwFromMrw( + func::function, func::dmn_variadic>& G_k_w, + double weight); double compute_S_k_w_from_G_k_w(); @@ -164,6 +183,7 @@ class CtauxClusterSolver { G0Interpolation g0_; typename Walker::Resource dummy_walker_resource_; + private: Rng rng_; @@ -226,6 +246,27 @@ void CtauxClusterSolver::write(Writer& writer) writer.close_group(); } +template +void CtauxClusterSolver::initialize( + int dca_iteration, SpRDisorderedGreensFunction& g0_disordered) { + dca_iteration_ = dca_iteration; + + g0_.initialize(g0_disordered); + + Sigma_old_ = data_.Sigma; + + accumulator_.initialize(dca_iteration_); + + averaged_ = false; + compute_jack_knife_ = + (dca_iteration == parameters_.get_dca_iterations() - 1) && + (parameters_.get_error_computation_type() == ErrorComputationType::JACK_KNIFE); + + if (concurrency_.id() == concurrency_.first()) + std::cout << "\n\n\t CT-AUX Integrator has initialized (DCA-iteration : " << dca_iteration + << ")\n\n"; +} + template void CtauxClusterSolver::initialize(int dca_iteration) { dca_iteration_ = dca_iteration; @@ -289,15 +330,24 @@ void CtauxClusterSolver::integrate() { } template -template -double CtauxClusterSolver::finalize( - dca_info_struct_t& dca_info_struct) { +void CtauxClusterSolver::accumulateGkw(double weight) { + collect_measurements(); + accumulateGkwfromMrw(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); @@ -636,6 +686,48 @@ 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, + double weight) { + 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 * weight; + // 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/ctaux/ctaux_walker.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp index 54c117ddc..11c31581d 100644 --- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp +++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp @@ -352,7 +352,7 @@ class CtauxWalker : public WalkerBIT, public CtauxWalkerData, 2> M_; + // std::array, 2> M_c; std::array, 2> exp_v_minus_one_; std::array, 2> exp_v_minus_one_dev_; std::array m_computed_events_; @@ -454,7 +454,6 @@ CtauxWalker::CtauxWalker(Parameters& parameters_ref, } } - template CtauxWalker::CtauxWalker(Parameters& parameters_ref, Data& MOMS_ref, Rng& rng_ref, int id) @@ -518,7 +517,6 @@ CtauxWalker::CtauxWalker(Parameters& parameters_ref, } } - template void CtauxWalker::printSummary() const { // std::defaultfloat is only supported by GCC 5 or later. @@ -923,13 +921,13 @@ int CtauxWalker::generateDelayedSpinsAbortAtBennett( } } -// #ifndef NDEBUG -// if (delayed_spins.size() > max_num_delayed_spins) { -// std::cout << "Delayed spins = " << delayed_spins.size() << " max_num_delayed_spins = " << max_num_delayed_spins << '\n'; -// if (delayed_spins.size() >= 256) -// throw std::runtime_error("delayed spins hit 256 or greater bailing out!"); -// } -// #endif + // #ifndef NDEBUG + // if (delayed_spins.size() > max_num_delayed_spins) { + // std::cout << "Delayed spins = " << delayed_spins.size() << " max_num_delayed_spins = " << + // max_num_delayed_spins << '\n'; if (delayed_spins.size() >= 256) + // throw std::runtime_error("delayed spins hit 256 or greater bailing out!"); + // } + // #endif // We need to unmark all "virtual" interacting spins, that we have temporarily marked as // annihilatable in CT_AUX_HS_configuration::get_random_noninteracting_vertex(). // TODO: Eliminate the need to mark and unmark these spins. @@ -940,11 +938,12 @@ int CtauxWalker::generateDelayedSpinsAbortAtBennett( assert(single_spin_updates_proposed == currently_proposed_creations_ + currently_proposed_annihilations_ + num_statics); -// #ifndef NDEBUG -// if (single_spin_updates_proposed >= max_num_delayed_spins) { -// std::cout << "single_spin_updates_proposed = " << single_spin_updates_proposed << " max_num_delayed_spins = " << max_num_delayed_spins << '\n'; -// } -// #endif + // #ifndef NDEBUG + // if (single_spin_updates_proposed >= max_num_delayed_spins) { + // std::cout << "single_spin_updates_proposed = " << single_spin_updates_proposed << " + // max_num_delayed_spins = " << max_num_delayed_spins << '\n'; + // } + // #endif return single_spin_updates_proposed; } @@ -1124,8 +1123,9 @@ void CtauxWalker::read_Gamma_matrices(e_spin_states case e_DN: linalg::util::syncStream(thread_id, stream_id); #ifndef NDEBUG - if(Gamma_dn.capacity().first < vertex_indixes.size()) - std::cout << "gamma_dn.capacity().first == " << Gamma_dn.capacity().first << " vertex_indixes.size() == " << vertex_indixes.size() << '\n'; + if (Gamma_dn.capacity().first < vertex_indixes.size()) + std::cout << "gamma_dn.capacity().first == " << Gamma_dn.capacity().first + << " vertex_indixes.size() == " << vertex_indixes.size() << '\n'; #endif CT_AUX_WALKER_TOOLS::compute_Gamma( Gamma_dn, N_dn, G_dn, vertex_indixes, exp_V, exp_delta_V, thread_id, stream_id); @@ -1340,21 +1340,21 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde #ifndef NDEBUG // The espin_up corresponds the Gamma matrix and diags used. - auto isDebugTestMinMax = [this](bool HS_field_up, bool espin_up, auto& ratio_HS_field, auto& Gamma_index_HS_field, auto& Gamma_CPU, auto& Gamma_diag_max, auto& Gamma_diag_min, int trace_depth) { + auto isDebugTestMinMax = [this](bool HS_field_up, bool espin_up, auto& ratio_HS_field, + auto& Gamma_index_HS_field, auto& Gamma_CPU, auto& Gamma_diag_max, + auto& Gamma_diag_min, int trace_depth) { if (std::abs(ratio_HS_field) <= std::abs(Scalar(1.e-16))) - this->pushOnTrace("gamma matrix was ill conditioned!"); - if (Gamma_index_HS_field > 0 && std::abs(ratio_HS_field) > 1.e-16 && - !ctaux_tools.test_max_min(Gamma_index_HS_field, Gamma_CPU, Gamma_diag_max, - Gamma_diag_min)) { - std::cerr - << "ratiod: " << ratio_HS_field - << "After e_spin_HS_field_ " << (HS_field_up ? "UP" : "DN") - << " == e_" << (espin_up ? "UP" : "DN") << ", solve_Gamma test_max_min on Gamma_LU failed!\n"; - this->dumpTrace(trace_, parameters_, thread_id, trace_depth); - this->last_failed_min_max_ = true; - } - else - this->last_failed_min_max_ = false; + this->pushOnTrace("gamma matrix was ill conditioned!"); + if (Gamma_index_HS_field > 0 && std::abs(ratio_HS_field) > 1.e-16 && + !ctaux_tools.test_max_min(Gamma_index_HS_field, Gamma_CPU, Gamma_diag_max, Gamma_diag_min)) { + std::cerr << "ratiod: " << ratio_HS_field << "After e_spin_HS_field_ " + << (HS_field_up ? "UP" : "DN") << " == e_" << (espin_up ? "UP" : "DN") + << ", solve_Gamma test_max_min on Gamma_LU failed!\n"; + this->dumpTrace(trace_, parameters_, thread_id, trace_depth); + this->last_failed_min_max_ = true; + } + else + this->last_failed_min_max_ = false; }; #endif @@ -1390,7 +1390,8 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde #ifndef NDEBUG if (Gamma_index_HS_field_DN != Gamma_up_size) { - std::cerr << "Gamma_index_HS_field_DN = " << Gamma_index_HS_field_DN << " Gamma_up_size = " << Gamma_up_size << '\n'; + std::cerr << "Gamma_index_HS_field_DN = " << Gamma_index_HS_field_DN + << " Gamma_up_size = " << Gamma_up_size << '\n'; } #endif if constexpr (Lattice::spin_symmetric) @@ -1401,9 +1402,10 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde std::cout << "warning Gamma_up delayed spin for spin_symmetric=false model"; #ifndef NDEBUG - isDebugTestMinMax(false, true, ratio_HS_field_DN, Gamma_index_HS_field_DN, Gamma_up_CPU, Gamma_up_diag_max, Gamma_up_diag_min, trace_depth); + isDebugTestMinMax(false, true, ratio_HS_field_DN, Gamma_index_HS_field_DN, Gamma_up_CPU, + Gamma_up_diag_max, Gamma_up_diag_min, trace_depth); #endif - + Gamma_up_size += 1; } else { @@ -1412,7 +1414,8 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde #ifndef NDEBUG if (Gamma_index_HS_field_DN != Gamma_dn_size) { - std::cerr << "Gamma_index_HS_field_DN = " << Gamma_index_HS_field_DN << " Gamma_dn_size = " << Gamma_dn_size << '\n'; + std::cerr << "Gamma_index_HS_field_DN = " << Gamma_index_HS_field_DN + << " Gamma_dn_size = " << Gamma_dn_size << '\n'; } #endif @@ -1420,9 +1423,9 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde exp_delta_V_HS_field_DN, Gamma_dn_diag_max, Gamma_dn_diag_min); - #ifndef NDEBUG - isDebugTestMinMax(false, false, ratio_HS_field_DN, Gamma_index_HS_field_DN, Gamma_dn_CPU, Gamma_dn_diag_max, Gamma_dn_diag_min, trace_depth); + isDebugTestMinMax(false, false, ratio_HS_field_DN, Gamma_index_HS_field_DN, Gamma_dn_CPU, + Gamma_dn_diag_max, Gamma_dn_diag_min, trace_depth); #endif Gamma_dn_size += 1; @@ -1434,7 +1437,8 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde #ifndef NDEBUG if (Gamma_index_HS_field_UP != Gamma_up_size) { - std::cerr << "Gamma_index_HS_field_UP = " << Gamma_index_HS_field_UP << " Gamma_up_size = " << Gamma_up_size << '\n'; + std::cerr << "Gamma_index_HS_field_UP = " << Gamma_index_HS_field_UP + << " Gamma_up_size = " << Gamma_up_size << '\n'; } #endif @@ -1446,7 +1450,8 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde std::cout << "warning Gamma_up delayed spin for spin_symmetric=false model"; #ifndef NDEBUG - isDebugTestMinMax(true, true, ratio_HS_field_UP, Gamma_index_HS_field_UP, Gamma_up_CPU, Gamma_up_diag_max, Gamma_up_diag_min, trace_depth); + isDebugTestMinMax(true, true, ratio_HS_field_UP, Gamma_index_HS_field_UP, Gamma_up_CPU, + Gamma_up_diag_max, Gamma_up_diag_min, trace_depth); #endif Gamma_up_size += 1; } @@ -1456,16 +1461,18 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde #ifndef NDEBUG if (Gamma_index_HS_field_UP != Gamma_dn_size) { - std::cerr << "Gamma_index_HS_field_UP = " << Gamma_index_HS_field_UP << " Gamma_dn_size = " << Gamma_dn_size << '\n'; + std::cerr << "Gamma_index_HS_field_UP = " << Gamma_index_HS_field_UP + << " Gamma_dn_size = " << Gamma_dn_size << '\n'; } #endif - + ratio_HS_field_UP = ctaux_tools.solve_Gamma_blocked(Gamma_index_HS_field_UP, Gamma_dn_CPU, exp_delta_V_HS_field_UP, Gamma_dn_diag_max, Gamma_dn_diag_min); #ifndef NDEBUG - isDebugTestMinMax(true, false, ratio_HS_field_UP, Gamma_index_HS_field_UP, Gamma_dn_CPU, Gamma_dn_diag_max, Gamma_dn_diag_min, trace_depth); + isDebugTestMinMax(true, false, ratio_HS_field_UP, Gamma_index_HS_field_UP, Gamma_dn_CPU, + Gamma_dn_diag_max, Gamma_dn_diag_min, trace_depth); #endif Gamma_dn_size += 1; } @@ -1523,34 +1530,34 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde std::cout << "Rejected but last vertex insert failed min_max test!\n"; #endif - /* - Gamma_up_diag_max = tmp_up_diag_max<1. ? 1. : tmp_up_diag_max; - Gamma_dn_diag_max = tmp_dn_diag_max<1. ? 1. : tmp_dn_diag_max; - Gamma_up_diag_min = tmp_up_diag_min>1. ? 1. : tmp_up_diag_min; - Gamma_dn_diag_min = tmp_dn_diag_min>1. ? 1. : tmp_dn_diag_min; + /* + Gamma_up_diag_max = tmp_up_diag_max<1. ? 1. : tmp_up_diag_max; + Gamma_dn_diag_max = tmp_dn_diag_max<1. ? 1. : tmp_dn_diag_max; + Gamma_up_diag_min = tmp_up_diag_min>1. ? 1. : tmp_up_diag_min; + Gamma_dn_diag_min = tmp_dn_diag_min>1. ? 1. : tmp_dn_diag_min; - delayed_spins[delayed_index].is_accepted_move = false; + delayed_spins[delayed_index].is_accepted_move = false; - if(delayed_spins[delayed_index].e_spin_HS_field_DN == e_UP) - { - CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_up_CPU, Gamma_up_size-1); - } - else - { - CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_dn_CPU, Gamma_dn_size-1); - } + if(delayed_spins[delayed_index].e_spin_HS_field_DN == e_UP) + { + CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_up_CPU, Gamma_up_size-1); + } + else + { + CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_dn_CPU, Gamma_dn_size-1); + } - if(delayed_spins[delayed_index].e_spin_HS_field_UP == e_UP) - { - CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_up_CPU, Gamma_up_size-1); - } - else - { - CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_dn_CPU, Gamma_dn_size-1); - } - */ + if(delayed_spins[delayed_index].e_spin_HS_field_UP == e_UP) + { + CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_up_CPU, Gamma_up_size-1); + } + else + { + CT_AUX_WALKER_TOOLS::set_to_identity(Gamma_dn_CPU, Gamma_dn_size-1); + } + */ - // delayed_spins[delayed_index].is_accepted_move = false; + // delayed_spins[delayed_index].is_accepted_move = false; #ifndef NDEBUG auto getNewDiagMin = [](auto& Gamma, auto gamma_size, auto& gamma_min) -> Scalar { @@ -1616,7 +1623,7 @@ void CtauxWalker::add_delayed_spin(int& delayed_inde pushOnTrace("double delayed_spins[", std::to_string(delayed_index), "].e_spin_HS_field_[DN,UP] == [e_UP,e_DN] G_dn", std::to_string(Gamma_dn_size - 1), ", G_up:", std::to_string(Gamma_up_size - 1)); -#endif //NDEBUG +#endif // NDEBUG } if (delayed_spins[delayed_index].e_spin_HS_field_DN == e_DN and 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/ctint/walker/ctint_walker_base.hpp b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp index 469785075..405ab627d 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 { diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp index 66b31d2f7..fd8e865d0 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp @@ -54,11 +54,15 @@ class SpAccumulator { using RDmn = typename Parameters::RClusterDmn; using NuDmn = func::dmn_variadic; // orbital-spin index - using PDmn = func::dmn_variadic; +#ifdef DISORDERED_G0 + using PDmn = func::dmn_variadic; using MFunction = - func::function, func::dmn_variadic>; - + func::function, func::dmn_variadic>; +#else + using PDmn = func::dmn_variadic; + using MFunction = func::function, func::dmn_variadic>; +#endif constexpr static int oversampling = 8; using NfftType = math::nfft::Dnfft1D; using MFunctionTime = NfftType; @@ -187,8 +191,11 @@ void SpAccumulator::accumulate( const Real t_i = config[i].get_tau(); const int delta_r = RDmn::parameter_type::subtract(r_j, r_i); const Real scaled_tau = (t_i - t_j) * one_div_two_beta; // + (i == j) * epsilon; - +#ifdef DISORDER_G0 + const int index = bbr_dmn(b_i, r_i, b_j, r_j); +#else const int index = bbr_dmn(b_i, b_j, delta_r); +#endif const Scalar f_val = Ms[s](i, j); (*cached_nfft_obj_)[s].accumulate(index, scaled_tau, factor * f_val); @@ -204,7 +211,7 @@ void SpAccumulator::accumulate( template void SpAccumulator::finalizeFunction(MFunctionTimePair& ft_objs, - MFunction& function) { + MFunction& function) { func::function, func::dmn_variadic> tmp("tmp"); const Real normalization = 1. / RDmn::dmn_size(); 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..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 { @@ -43,11 +44,17 @@ template auto shrinkG0(const SpGreensFunction& G0) { func::function, TDmn>> g0_trimmed; const int s = 0; + 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++) + for (int t = 0; t < TDmn::dmn_size(); t++) { + 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; } 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..0fed9f13f 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 @@ -59,6 +59,8 @@ class StdThreadQmciClusterSolver : public QmciSolver { using typename BaseClass::Accumulator; using Walker = stdthreadqmci::StdThreadQmciWalker; using SpGreensFunction = typename BaseClass::SpGreensFunction; + using SpDisorderedGreensFunction = typename BaseClass::SpDisorderedGreensFunction; + using SpRDisorderedGreensFunction = typename BaseClass::SpRDisorderedGreensFunction; using StdThreadAccumulatorType = stdthreadqmci::StdThreadQmciAccumulator; using MFunction = typename StdThreadAccumulatorType::MFunction; @@ -75,6 +77,7 @@ class StdThreadQmciClusterSolver : public QmciSolver { const std::shared_ptr>& writer); void initialize(int dca_iteration); + void initialize(int dca_iteration, SpRDisorderedGreensFunction& g0_disordered); void integrate(); @@ -116,6 +119,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(double weight); + /// This collects measurements when dealing with single G0 + void collectSingle(); + private: void startWalker(int id); void startAccumulator(int id, const Parameters& parameters); @@ -123,8 +131,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); @@ -182,7 +198,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( @@ -227,6 +243,21 @@ void StdThreadQmciClusterSolver::initialize(int dca_iteration) { measurements_done_ = 0; } +template +void StdThreadQmciClusterSolver::initialize(int dca_iteration, + SpRDisorderedGreensFunction& g0_disordered) { + Profiler profiler(__FUNCTION__, "stdthread-MC-Integration", __LINE__); + + last_iteration_ = dca_iteration == parameters_.get_dca_iterations() - 1; + + measurements_ = parameters_.get_measurements()[dca_iteration]; + + BaseClass::initialize(dca_iteration, g0_disordered); + + walk_finished_ = 0; + measurements_done_ = 0; +} + template void StdThreadQmciClusterSolver::integrate() { Profiler profiler(__FUNCTION__, "stdthread-MC-Integration", __LINE__); @@ -261,7 +292,6 @@ void StdThreadQmciClusterSolver::integrate() { dca::profiling::Duration duration(end_time, start_time); total_time_ = duration.sec + 1.e-6 * duration.usec; - printIntegrationMetadata(); }; @@ -276,7 +306,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"); @@ -347,8 +377,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(); @@ -367,8 +397,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 +529,16 @@ auto StdThreadQmciClusterSolver::computeSingleMeasurement_G_k_w( return G_k_w; } +template +void StdThreadQmciClusterSolver::accumulateGkw(double weight) { + QmciSolver::accumulateGkw(weight); +} + +template +void StdThreadQmciClusterSolver::collectSingle() { + QmciSolver::collectSingle(); +} + template void StdThreadQmciClusterSolver::logSingleMeasurement( StdThreadAccumulatorType& accumulator_obj, int stamping_period, bool log_MFunction, @@ -589,8 +629,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 +681,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(); @@ -678,13 +718,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) @@ -696,15 +739,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) 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/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 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/domains/time_and_frequency/time_domain.hpp b/include/dca/phys/domains/time_and_frequency/time_domain.hpp index f7bdb02e5..e6b17655a 100644 --- a/include/dca/phys/domains/time_and_frequency/time_domain.hpp +++ b/include/dca/phys/domains/time_and_frequency/time_domain.hpp @@ -32,6 +32,7 @@ class time_domain { using element_type = double; using scalar_type = double; + static constexpr int DIMENSION = 1; // Needed in function transform. using dmn_specifications_type = math::transform::interval_dmn_1D_type; @@ -114,8 +115,8 @@ void time_domain::to_JSON(stream_type& ss) { ss << "]\n"; } -} // domains -} // phys -} // dca +} // namespace domains +} // namespace phys +} // namespace dca #endif // DCA_PHYS_DOMAINS_TIME_AND_FREQUENCY_TIME_DOMAIN_HPP diff --git a/include/dca/phys/parameters/disorder_parameters.hpp b/include/dca/phys/parameters/disorder_parameters.hpp new file mode 100644 index 000000000..39758f1e5 --- /dev/null +++ b/include/dca/phys/parameters/disorder_parameters.hpp @@ -0,0 +1,127 @@ +// 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 + * 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: + 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); + + 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: + 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(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; +} + +template +void DisorderParameters::pack(const Concurrency& concurrency, char* buffer, int buffer_size, + int& position) const { + 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, 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 +void DisorderParameters::readWrite(ReaderOrWriter& reader_or_writer) { + try { + reader_or_writer.open_group("disorder"); + + try { + 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("num-configurations", disorder_num_configurations_); + } + catch (const std::exception& r_e) { + } + try { + reader_or_writer.execute("max-sites", disorder_max_sites_); + } + 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/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/include/dca/phys/parameters/mci_parameters.hpp b/include/dca/phys/parameters/mci_parameters.hpp index 8624733bf..ddcf5e1b6 100644 --- a/include/dca/phys/parameters/mci_parameters.hpp +++ b/include/dca/phys/parameters/mci_parameters.hpp @@ -160,6 +160,8 @@ class MciParameters { bool fix_meas_per_walker_ = true; bool adjust_self_energy_for_double_counting_; ErrorComputationType error_computation_type_; + /// actually a runtime variable that should be moved from + /// parameters. bool store_configuration_; DistType g4_distribution_; int stamping_period_ = 0; @@ -322,41 +324,43 @@ void MciParameters::readWrite(ReaderOrWriter& reader_or_writer) { reader_or_writer.close_group(); if constexpr (ReaderOrWriter::is_reader) { - // The input file can contain an integral seed or the seeding option "random". - // Check parameters consistency. - if (g4_distribution_ == DistType::BLOCKED) { + // The input file can contain an integral seed or the seeding option "random". + // Check parameters consistency. + if (g4_distribution_ == DistType::BLOCKED) { #ifdef DCA_HAVE_MPI - // Check for number of accumulators and walkers consistency. - if (!shared_walk_and_accumulation_thread_ || walkers_ != accumulators_) { - throw std::logic_error( - "\n With distributed g4 enabled, 1) walker and accumulator should share " - "thread, " - "2) #walker == #accumulator\n"); - } + // Check for number of accumulators and walkers consistency. + if (!shared_walk_and_accumulation_thread_ || walkers_ != accumulators_) { + throw std::logic_error( + "\n With distributed g4 enabled, 1) walker and accumulator should share " + "thread, " + "2) #walker == #accumulator\n"); + } - // Check for number of ranks and g4 measurements consistency. - // This is potentially wrong if fancy rank ganging or the like is used. - int mpi_size; - MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); - int local_meas = measurements_.back() / mpi_size; - if (measurements_.back() % mpi_size != 0 && local_meas % accumulators_ != 0) { - throw std::logic_error( - "\n With distributed g4 enabled, 1) local measurements should be same across " - "ranks, " - "2) each accumulator should have same measurements\n"); - } + // Check for number of ranks and g4 measurements consistency. + // This is potentially wrong if fancy rank ganging or the like is used. + int mpi_size; + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + int local_meas = measurements_.back() / mpi_size; + if (measurements_.back() % mpi_size != 0 && local_meas % accumulators_ != 0) { + throw std::logic_error( + "\n With distributed g4 enabled, 1) local measurements should be same across " + "ranks, " + "2) each accumulator should have same measurements\n"); + } #else - throw(std::logic_error("MPI distribution requested with no MPI available.")); + throw(std::logic_error("MPI distribution requested with no MPI available.")); #endif // DCA_HAVE_MPI - if (stamping_period_ != 0) { - if (!(shared_walk_and_accumulation_thread_ && walkers_ == accumulators_)) - throw std::runtime_error("Individual measurement stamping not available unless shared-walk-and-accumulation-thread = true and walkers == acceptors!"); + if (stamping_period_ != 0) { + if (!(shared_walk_and_accumulation_thread_ && walkers_ == accumulators_)) + throw std::runtime_error( + "Individual measurement stamping not available unless " + "shared-walk-and-accumulation-thread = true and walkers == acceptors!"); + } + // Solve conflicts } - // Solve conflicts - } } if (!time_correlation_window_) - compute_G_correlation_ = false; + compute_G_correlation_ = false; } void MciParameters::solveDcaIterationConflict(int iterations) { diff --git a/include/dca/phys/parameters/parameters.hpp b/include/dca/phys/parameters/parameters.hpp index 56c4be6ff..89e5c6dcc 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 @@ -176,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(), @@ -202,8 +207,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); } @@ -451,6 +457,10 @@ void Parameters::readWrite(reader_or_writer); 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 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 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/dca_loop/disorder/apply_disorder.cpp b/src/phys/dca_loop/disorder/apply_disorder.cpp new file mode 100644 index 000000000..c771c470a --- /dev/null +++ b/src/phys/dca_loop/disorder/apply_disorder.cpp @@ -0,0 +1,14 @@ +// 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" + +template +auto makeDisorderedG0(decltype(dca::phys::DcaData::G0_r_t_cluster_excluded)& g0_r_t_cl_exl) {} 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_; diff --git a/src/phys/domains/cluster/CMakeLists.txt b/src/phys/domains/cluster/CMakeLists.txt index efa0df526..557719833 100644 --- a/src/phys/domains/cluster/CMakeLists.txt +++ b/src/phys/domains/cluster/CMakeLists.txt @@ -4,5 +4,6 @@ add_library(cluster_domains STATIC cluster_definitions.cpp momentum_exchange_domain.cpp symmetries/symmetry_operations/group_action.cpp) +# value_at_debug.cpp) dca_gpu_runtime_link(cluster_domains) 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/function_transform/function_transform_test.cpp b/test/unit/math/function_transform/function_transform_test.cpp index 6a1f68dba..5bab7ae53 100644 --- a/test/unit/math/function_transform/function_transform_test.cpp +++ b/test/unit/math/function_transform/function_transform_test.cpp @@ -31,15 +31,16 @@ using McOptions = MockMcOptions; #include "dca/parallel/no_concurrency/no_concurrency.hpp" #include "dca/parallel/no_threading/no_threading.hpp" #include "dca/profiling/null_profiler.hpp" +#include "dca/phys/domains/domain_aliases.hpp" using Model = dca::phys::models::TightBindingModel>; using Concurrency = dca::parallel::NoConcurrency; -using Parameters = - dca::phys::params::Parameters::type>>; +using Parameters = dca::phys::params::Parameters< + Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler, Model, void, + dca::ClusterSolverId::CT_AUX, + dca::NumericalTraits< + double, typename dca::util::ScalarSelect::type>>; const std::string input_dir = DCA_SOURCE_DIR "/test/unit/math/function_transform/"; @@ -49,6 +50,8 @@ using NuDmn = dca::func::dmn_variadic; using WDmn = dca::func::dmn_0; using KDmn = Parameters::KClusterDmn; using RDmn = Parameters::RClusterDmn; +using DDA = dca::phys::DcaDomainAliases; +using TDmn = typename DDA::TDmn; const std::vector> a_vecs{std::vector{0, 0}, std::vector{0.25, 0.25}}; @@ -141,6 +144,80 @@ void spTestImplementation(const bool direct) { } } +template +void timeFreqTestImplementation(const bool direct) { + using namespace dca::func; + using Real = double; + using Complex = std::complex; + function> f_b_b_in; + function> f_in; + + // Initialize the input function. + for (int k = 0; k < KDmn::dmn_size(); ++k) + for (int i = 0; i < InpDmn::dmn_size(); ++i) { + const Complex val(i * i + 0.5 * i, k); + f_in(i, k) = val; + for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2) + for (int b1 = 0; b1 < BDmn::dmn_size(); ++b1) + f_b_b_in(b1, 0, b2, 0, k, i) = f_b_b_in(b1, 1, b2, 1, k, i) = val; + } + + function> f_b_b_out; + function> f_out; + + // Regular transform. + dca::math::transform::FunctionTransform::execute(f_in, f_out); + // Transform with phase factor. + dca::math::transform::FunctionTransform::execute(f_b_b_in, f_b_b_out); + + const int exp_sign = direct ? 1 : -1; + const double norm = direct ? 1 : 1. / InpDmn::dmn_size(); + + constexpr double tolerance = 1e-10; + + for (int k = 0; k < KDmn::dmn_size(); ++k) { + // Test regular transform. + for (int o = 0; o < OutDmn::dmn_size(); ++o) { + const auto& o_dmn_val = OutDmn::get_elements()[o]; + Complex val(0); + + for (int i = 0; i < InpDmn::dmn_size(); ++i) { + const auto& i_dmn_val = InpDmn::get_elements()[i]; + val += f_in(i, k) * std::exp(Complex(0., exp_sign * dca::math::util::innerProduct( + static_cast(i_dmn_val), + static_cast(o_dmn_val)))); + } + val *= norm; + + EXPECT_LE(std::abs(val - f_out(o, k)), tolerance); + + // // Test transform with phase factors. + // for (int b1 = 0; b1 < BDmn::dmn_size(); ++b1) { + // for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2) { + // const auto a_diff = dca::math::util::subtract(a_vecs[b2], a_vecs[b1]); + // for (int o = 0; o < OutDmn::dmn_size(); ++o) { + // const auto& o_dmn_val = OutDmn::get_elements()[o]; + // Complex val(0); + + // for (int i = 0; i < InpDmn::dmn_size(); ++i) { + // const auto& i_dmn_val = InpDmn::get_elements()[i]; + // const auto& k_val = direct ? o_dmn_val : i_dmn_val; + // const Complex phase_factor = std::exp(Complex( + // 0, exp_sign * dca::math::util::innerProduct(static_cast(k_val), a_diff))); + // val += f_b_b_in(b1, 0, b2, 0, k, i) * phase_factor * + // std::exp(Complex( + // 0., exp_sign * dca::math::util::innerProduct(i_dmn_val, o_dmn_val))); + // } + // val *= norm; + + // EXPECT_LE(std::abs(val - f_b_b_out(b1, 0, b2, 0, k, o)), tolerance); + // } + // } + // } + } + } +} + TEST(FunctionTransformTest, SpaceToMomentumCmplx) { initialize(); spTestImplementation(true); @@ -150,3 +227,8 @@ TEST(FunctionTransformTest, MomentumToSpaceCmplx) { initialize(); spTestImplementation(false); } + +// TEST(FunctionTransformTest, TimeToFrequency) { +// initialize(); +// timeFreqTestImplementation(true); +// } 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/CMakeLists.txt b/test/unit/phys/CMakeLists.txt index 9d64b1f5e..a6ab19113 100644 --- a/test/unit/phys/CMakeLists.txt +++ b/test/unit/phys/CMakeLists.txt @@ -6,3 +6,4 @@ add_subdirectory(dca_step) add_subdirectory(domains) add_subdirectory(models) add_subdirectory(parameters) +add_subdirectory(dca_loop) 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/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) 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 + } +}