Diagonal Disorder

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()

applications/dca/main_dca.cpp

@@ -15,6 +15,7 @@
 #include <iostream>
 
 #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<dca::io::JSONReader>(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();

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<double*>(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<double*>(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

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

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 <algorithm>
 #include <cassert>
 #include <complex>
+#include <cstddef>
 #include <cstring>
 #include <iostream>
 #include <string>
@@ -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
@@ -111,6 +115,16 @@ class DcaData {
                                         KClusterDmn, WVertexDmn, KExchangeDmn, WExchangeDmn>,
                      DT>;
 
+  using Type_G0_k_w =
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, WDmn>>;
+  using Type_G0_k_t =
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, TDmn>>;
+  using Type_G0_r_t = func::function<Scalar, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, TDmn>>;
+
+  // Vector NuDmn * RDmn in size that is the disorder configuration
+  using DST = DcaSharedTypes<Parameters>;
+  using DisorderConfiguration = typename DST::DisorderConfiguration;
+
   DcaData(Parameters& parameters_ref);
 
   /** These reads are used by analysis programs only for now.
@@ -186,18 +200,29 @@ class DcaData {
   func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>> G_r_w;
   func::function<Scalar, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, TDmn>> G_r_t;
 
+  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, WDmn>> accumulated_G_k_w;
+  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, TDmn>> accumulated_G_k_t;
+  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>> accumulated_G_r_w;
+  func::function<Scalar, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, TDmn>> accumulated_G_r_t;
+
   func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, WDmn>> G0_k_w;
   func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, TDmn>> G0_k_t;
   func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>> G0_r_w;
   func::function<Scalar, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, TDmn>> G0_r_t;
 
-  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, WDmn>>
-      G0_k_w_cluster_excluded;
-  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, TDmn>>
-      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<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>>
       G0_r_w_cluster_excluded;
-  func::function<Scalar, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, TDmn>> 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;
+
+  /// This G0 could be modified by the diagonal disorder scheme, but
+  /// normally will just be a copy of the G0_r_t_cluster_excluded
+  Type_G0_r_t mutable_G0_r_t_cluster_excluded;
+  Type_G0_k_w mutable_G0_k_w_cluster_excluded;
+  Type_G0_k_t mutable_G0_k_t_cluster_excluded;
 
   func::function<Real, NuDmn> orbital_occupancy;
 
@@ -261,7 +286,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, Type_G0_r_t& disordered_G0_r_t_cl_exl);
+
+  /// Due to the new ability to modify the G0 for disorder this is the
+  /// immutable g0 from the last iteration.
   std::unique_ptr<SpGreensFunction> G_k_w_err_;
   std::unique_ptr<SpRGreensFunction> G_r_w_err_;
   std::unique_ptr<SpGreensFunction> Sigma_err_;
@@ -365,19 +397,19 @@ void DcaData<Parameters, DT>::read(const std::string& filename) {
   if (parameters_.isAccumulatingG4()) {
     concurrency_.broadcast_object(G_k_w);
 #ifndef NDEBUG
-  if (concurrency_.id() == concurrency_.first()) {
-    std::cout << "broadcasted G_k_w \n";
-  }
+    if (concurrency_.id() == concurrency_.first()) {
+      std::cout << "broadcasted G_k_w \n";
+    }
 #endif
 
-  for (auto& G4_channel : G4_) {
+    for (auto& G4_channel : G4_) {
       concurrency_.broadcast_object(G4_channel);
 #ifndef NDEBUG
-  if (concurrency_.id() == concurrency_.first()) {
-    std::cout << "broadcasted G4_channel \n";
-  }
+      if (concurrency_.id() == concurrency_.first()) {
+        std::cout << "broadcasted G4_channel \n";
+      }
 #endif
-  }
+    }
   }
 }
 
@@ -555,7 +587,8 @@ void DcaData<Parameters, DT>::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 +599,6 @@ void DcaData<Parameters, DT>::initializeH0_and_H_i() {
   }
 
   Parameters::model_type::initialize_H_symmetries(H_symmetry);
-
   compute_band_structure<Parameters>::execute(parameters_, band_structure);
 }
 
@@ -600,6 +632,9 @@ void DcaData<Parameters, DT>::initialize_G0() {
     G0_k_t_cluster_excluded = G0_k_t;
     G0_r_w_cluster_excluded = G0_r_w;
     G0_r_t_cluster_excluded = G0_r_t;
+    // When there is no disorder in G0 mutable just gets the unmodified G0_r_t
+    mutable_G0_r_t_cluster_excluded = G0_r_t;
+    mutable_G0_k_w_cluster_excluded = G0_k_w;
   }
   catch (const std::exception& exc) {
     std::throw_with_nested(std::runtime_error("Failure in initialization of G0!"));
@@ -659,13 +694,13 @@ void DcaData<Parameters, DT>::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 +732,73 @@ void DcaData<Parameters, DT>::readSigmaFile(io::Reader<Concurrency>& reader) {
   reader.close_group();
 }
 
+template <class Parameters, DistType DT>
+void DcaData<Parameters, DT>::makeDisorderedG0(const DisorderConfiguration& disorder_configuration,
+                                               const Type_G0_r_t& g0_r_t_cl_exl,
+                                               Type_G0_r_t& disordered_g0_r_t_cl_exl) {
+  int matrix_dim = dca::phys::DcaData<Parameters, DT>::NuDmn::dmn_size();
+  dca::linalg::Matrix<Scalar, dca::linalg::CPU,
+                      dca::linalg::util::DefaultAllocator<Scalar, dca::linalg::CPU>>
+      g0_rtcex_inverse(matrix_dim);
+  // dca::linalg::Vector<int, dca::linalg::CPU,
+  //                     dca::linalg::util::DefaultAllocator<std::complex<double>, dca::linalg::CPU>>
+  //     ipiv;
+  // dca::linalg::Vector<std::complex<double>, dca::linalg::CPU,
+  //                     dca::linalg::util::DefaultAllocator<std::complex<double>, dca::linalg::CPU>>
+  //     work;
+
+  for (int ir = 0; ir < RClusterDmn::dmn_size(); ++ir)
+    for (int it = 0; it < TDmn::dmn_size(); ++it) {
+      dca::linalg::Vector<Scalar, dca::linalg::CPU,
+                          dca::linalg::util::DefaultAllocator<Scalar, dca::linalg::CPU>>
+      rt_block(std::size_t(NuDmn::dmn_size() * NuDmn::dmn_size()));
+      g0_r_t_cl_exl.slice(0, 1, {0, 0, ir, it}, rt_block.data());
+      dca::linalg::matrixop::copyArrayToMatrix(matrix_dim, matrix_dim, rt_block.data(), matrix_dim,
+                                               g0_rtcex_inverse);
+      dca::linalg::matrixop::inverse(g0_rtcex_inverse);  //, ipiv, work);
+
+      // Then apply disorder potential to diagonal according to the
+      // configuration
+      for (int imd = 0; imd < matrix_dim; ++imd) {
+        g0_rtcex_inverse(imd, imd) += disorder_configuration(imd, ir);
+      }
+
+      dca::linalg::matrixop::inverse(g0_rtcex_inverse);
+      dca::linalg::matrixop::copyMatrixToArray(g0_rtcex_inverse,
+                                               &disordered_g0_r_t_cl_exl(0, 0, ir, it), matrix_dim);
+    }
+
+  // Giving up at this point in side effect free method.
+  math::transform::FunctionTransform<RClusterDmn, KClusterDmn>::execute(
+      disordered_g0_r_t_cl_exl, mutable_G0_k_t_cluster_excluded);
+
+  // I can't figure out how to get this to work with the
+  // math::transform::FunctionTransform framwork so do this by hand
+  // here.
+  using Complex = std::complex<util::RealAlias<Scalar>>;
+  for (int k = 0; k < KClusterDmn::dmn_size(); ++k) {
+    const auto& k_val = KClusterDmn::get_elements()[k];
+    for (int inu1 = 0; inu1 < NuDmn::dmn_size(); ++inu1)
+      for (int inu2 = 0; inu2 < NuDmn::dmn_size(); ++inu2) {
+        for (int w = 0; w < WDmn::dmn_size(); ++w) {
+          const auto& w_val = WDmn::get_elements()[w];
+          Complex G_k_omega(0);
+          for (int t = 0; t < TDmn::dmn_size(); ++t) {
+            const auto& t_val = TDmn::get_elements()[t];
+            G_k_omega += mutable_G0_k_t_cluster_excluded(inu1, inu2, k, t) *
+                         std::exp(Complex(0, w_val * t_val));
+          }
+          mutable_G0_k_w_cluster_excluded(inu1, inu2, k, w) = G_k_omega;
+        }
+      }
+  }
+}
+
+template <class Parameters, DistType DT>
+void DcaData<Parameters, DT>::makeDisorderedG0(const DisorderConfiguration& disorder_configuration) {
+  makeDisorderedG0(disorder_configuration, G0_r_t_cluster_excluded, mutable_G0_r_t_cluster_excluded);
+}
+
 template <class Parameters, DistType DT>
 void DcaData<Parameters, DT>::compute_single_particle_properties() {
   {