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 <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();
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<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
diff --git a/include/dca/linalg/util/gpu_event.hpp b/include/dca/linalg/util/gpu_event.hpp
index bb7614616..891d59199 100644
--- a/include/dca/linalg/util/gpu_event.hpp
+++ b/include/dca/linalg/util/gpu_event.hpp
@@ -64,7 +64,7 @@ class GpuEvent {
 };
 
 // Returns the elapsed time in seconds between two recorded events. Blocks host.
-float elapsedTime(cudaEvent_t stop, cudaEvent_t start) {
+inline float elapsedTime(cudaEvent_t stop, cudaEvent_t start) {
   checkRC(cudaEventSynchronize(stop));
   float msec(0);
   checkRC(cudaEventElapsedTime(&msec, start, stop));
@@ -91,8 +91,8 @@ class GpuEvent {
 
 #endif  // DCA_HAVE_GPU
 
-}  // util
-}  // linalg
-}  // dca
+}  // namespace util
+}  // namespace linalg
+}  // namespace dca
 
 #endif  // DCA_LINALG_UTIL_GPU_EVENT_HPP
diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp
index d1fafc22b..064147263 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 <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() {
   {
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 <class Parameters, DistType DT>
+void DcaData<Parameters, DT>::makeDisorderedG0(
+    const DisorderConfiguration disorder_configuration,
+    const decltype(dca::phys::DcaData<Parameters, DIST>::G0_r_t_cluster_excluded)& g0_r_t_cl_exl,
+    decltype(dca::phys::DcaData<Parameters, DIST>::disordered_G0_r_t_cluster_excluded)&
+        disordered_g0_r_t_cl_exl) {
+  int matrix_dim = dca::phys::DcaData::NuDmn::dmn_size();
+  dca::linalg::Matrix<std::complex<double, dca::linalg::CPU>> g0_rtcex_inverse(matrix_dim);
+  dca::linalg::Vector<int, dca::linalg::CPU> ipiv;
+  dca::linalg::Vector<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::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<ParametersType> DCA_info_struct;
 
   cluster_exclusion_type cluster_exclusion_obj;
@@ -240,8 +249,14 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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<ParametersType, DcaDataType, MCIntegratorType, DIST>::execute() {
   }
 }
 
+template <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
+double DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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<ParametersType> 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 <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
 void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::finalize() {
   perform_cluster_mapping_self_energy();
@@ -361,7 +401,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::perform_clust
 }
 
 template <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
-double DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::solve_cluster_problem(
+void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::solve_cluster_problem(
     int DCA_iteration) {
   //  static_assert(std::is_same<DcaDataType, ::DcaDataType<DIST>>::value);
   // static_assert(std::is_same<MCIntegratorType, ::ClusterSolver<DIST>>::value);
@@ -377,18 +417,37 @@ double DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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 <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
+void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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 <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
+void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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 <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
+double DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::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 <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
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<domains::cluster_domain<double, Parameters::lattice_type::DIMENSION, domains::CLUSTER,
                                           domains::MOMENTUM_SPACE, domains::BRILLOUIN_ZONE>>;
+  using r_DCA =
+      func::dmn_0<domains::cluster_domain<double, Parameters::lattice_type::DIMENSION, domains::CLUSTER,
+                                          domains::REAL_SPACE, domains::BRILLOUIN_ZONE>>;
+
+  using DST = DcaSharedTypes<Parameters>;
+  using DisorderConfiguration = typename DST::DisorderConfiguration;
   DcaLoopData();
 
   template <typename WRITER>
@@ -85,6 +94,10 @@ class DcaLoopData {
   func::function<double, DCA_iteration_domain_type> chemical_potential;
   func::function<double, DCA_iteration_domain_type> average_expansion_order;
 
+  func::function<double, DCA_iteration_domain_type> average_defect_density;
+  func::function<int, DCA_iteration_domain_type> num_defect_configurations;
+  std::vector<DisorderConfiguration> disorder_configurations;
+  std::vector<double> disorder_weights;
   int last_completed_iteration = -1;
 };
 
@@ -113,7 +126,8 @@ DcaLoopData<Parameters>::DcaLoopData()
 
       density("density"),
       chemical_potential("chemical-potential"),
-      average_expansion_order("expansion_order") {}
+      average_expansion_order("expansion_order"),
+      average_defect_density("defect_density") {}
 
 template <typename Parameters>
 template <typename Writer>
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 <vector>
+#include <unordered_set>
+#include <random>
+#include "dca/phys/types/dca_shared_types.hpp"
+
+namespace dca::phys {
+
+namespace detail {
+struct VectorHash {
+  size_t operator()(const std::vector<int>& 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<int>{}(x) + golden_ratio + (seed << 6) + (seed >> 2);
+    }
+    return seed;
+  }
+};
+
+}  // namespace detail
+
+template <typename Parameters>
+class MakeDisorderConfigurations {
+  using DST = DcaSharedTypes<Parameters>;
+  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<DisorderConfiguration>& disorder_configurations,
+                  std::vector<double>& 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<std::vector<int>, 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<int> 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..c9102e584 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp
@@ -41,7 +41,7 @@ namespace solver {
 namespace ctaux {
 // dca::phys::solver::ctaux::
 using dca::util::SignType;
-  
+
 template <class Parameters, class Data>
 class TpEqualTimeAccumulator {
 public:
@@ -253,7 +253,7 @@ TpEqualTimeAccumulator<Parameters, Data>::TpEqualTimeAccumulator(const Parameter
 
       G_r_t("G_r_t"),
       G_r_t_stddev("G_r_t_stddev"),
-      
+
       G_r_t_accumulated("G_r_t_accumulated"),
       G_r_t_accumulated_squared("G_r_t_accumulated_squared"),
 
@@ -339,7 +339,7 @@ void TpEqualTimeAccumulator<Parameters, Data>::initialize_akima_coefficients() {
       for (int nu1_ind = 0; nu1_ind < nu::dmn_size(); nu1_ind++) {
         for (int nu0_ind = 0; nu0_ind < nu::dmn_size(); nu0_ind++) {
           for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++)
-            y[t_ind] = MOMS.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+            y[t_ind] = MOMS.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
 
           ai_obj.initialize(x, y);
 
@@ -357,7 +357,7 @@ void TpEqualTimeAccumulator<Parameters, Data>::initialize_akima_coefficients() {
         for (int nu0_ind = 0; nu0_ind < nu::dmn_size(); nu0_ind++) {
           for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size(); t_ind++)
             y[t_ind - t::dmn_size() / 2] =
-                MOMS.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+                MOMS.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
 
           ai_obj.initialize(x, y);
 
diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp
index b6edbe850..73b3e2e6f 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp
@@ -29,6 +29,7 @@
 #include "dca/function/domains.hpp"
 #include "dca/function/function.hpp"
 #include "dca/linalg/linalg.hpp"
+#include "dca/linalg/util/allocators/allocators.hpp"
 #include "dca/math/function_transform/function_transform.hpp"
 #include "dca/math/statistics/util.hpp"
 #include "dca/parallel/util/get_workload.hpp"
@@ -38,6 +39,7 @@
 #include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/time_correlator.hpp"
 #include "dca/phys/dca_step/symmetrization/symmetrize.hpp"
 #include "dca/phys/domains/cluster/cluster_domain.hpp"
+#include "dca/phys/domains/domain_aliases.hpp"
 #include "dca/phys/domains/quantum/electron_band_domain.hpp"
 #include "dca/phys/domains/quantum/electron_spin_domain.hpp"
 #include "dca/phys/domains/time_and_frequency/frequency_domain.hpp"
@@ -75,6 +77,9 @@ class CtauxClusterSolver {
   static constexpr linalg::DeviceType device = device_t;
 
 protected:
+  // The following are deprecated because they don't follow the
+  // naming rules nor are they consistent with other classes aliases
+  // for these domains.
   using w = func::dmn_0<domains::frequency_domain>;
   using b = func::dmn_0<domains::electron_band_domain>;
   using s = func::dmn_0<domains::electron_spin_domain>;
@@ -83,7 +88,9 @@ class CtauxClusterSolver {
   using CDA = ClusterDomainAliases<Parameters::lattice_type::DIMENSION>;
   using RDmn = typename CDA::RClusterDmn;
   using KDmn = typename CDA::KClusterDmn;
-
+  using DDA = DcaDomainAliases<Parameters>;
+  using WDmn = typename DDA::WDmn;
+  using NuDmn = typename DDA::NuDmn;
   using NuNuKClusterWDmn = func::dmn_variadic<nu, nu, KDmn, w>;
   using NuNuRClusterWDmn = func::dmn_variadic<nu, nu, RDmn, w>;
 
@@ -117,7 +124,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 +145,9 @@ class CtauxClusterSolver {
   void collect_measurements();
 
   void compute_G_k_w_from_M_r_w();
+  void accumulateGkwFromMrw(
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KDmn, WDmn>>& G_k_w,
+      double weight);
 
   double compute_S_k_w_from_G_k_w();
 
@@ -164,6 +179,7 @@ class CtauxClusterSolver {
   G0Interpolation<device, Scalar> g0_;
 
   typename Walker::Resource dummy_walker_resource_;
+
 private:
   Rng rng_;
 
@@ -230,7 +246,7 @@ template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistTy
 void CtauxClusterSolver<device_t, Parameters, Data, DIST>::initialize(int dca_iteration) {
   dca_iteration_ = dca_iteration;
 
-  g0_.initializeShrinked(data_.G0_r_t_cluster_excluded);
+  g0_.initializeShrinked(data_.mutable_G0_r_t_cluster_excluded);
 
   Sigma_old_ = data_.Sigma;
 
@@ -289,15 +305,24 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::integrate() {
 }
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
-template <typename dca_info_struct_t>
-double CtauxClusterSolver<device_t, Parameters, Data, DIST>::finalize(
-    dca_info_struct_t& dca_info_struct) {
+void CtauxClusterSolver<device_t, Parameters, Data, DIST>::accumulateGkw(double weight) {
+  collect_measurements();
+  accumulateGkwfromMrw(data_.G_k_w);
+}
+
+template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
+void CtauxClusterSolver<device_t, Parameters, Data, DIST>::collectSingle() {
   collect_measurements();
   symmetrize_measurements();
 
   // Compute new Sigma.
   compute_G_k_w_from_M_r_w();
+}
 
+template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
+template <typename dca_info_struct_t>
+double CtauxClusterSolver<device_t, Parameters, Data, DIST>::finalize(
+    dca_info_struct_t& dca_info_struct) {
   // FT<k_DCA,r_DCA>::execute(data_.G_k_w, data_.G_r_w);
   math::transform::FunctionTransform<KDmn, RDmn>::execute(data_.G_k_w, data_.G_r_w);
 
@@ -636,6 +661,48 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::compute_G_k_w_from_M_
   Symmetrize<Parameters>::execute(data_.G_k_w, data_.H_symmetry);
 }
 
+template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
+void CtauxClusterSolver<device_t, Parameters, Data, DIST>::accumulateGkwFromMrw(
+    func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KDmn, WDmn>>& G_k_w,
+    double weight) {
+  func::function<std::complex<double>, NuNuKClusterWDmn> M_k_w;
+  math::transform::FunctionTransform<RDmn, KDmn>::execute(M_r_w_, M_k_w);
+
+  const std::size_t matrix_size = b::dmn_size() * s::dmn_size();
+  linalg::Matrix<std::complex<Real>, dca::linalg::CPU> G0_times_M_matrix(matrix_size);
+  using MatrixView =
+      linalg::MatrixView<std::complex<Real>, dca::linalg::CPU,
+                         linalg::util::DefaultAllocator<std::complex<Real>, 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<Parameters>::execute(data_.G_k_w, data_.H_symmetry);
+}
+
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
 double CtauxClusterSolver<device_t, Parameters, Data, DIST>::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/walker/tools/g0_interpolation/g0_interpolation_base.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp
index 6d6ee0d8a..6b3d990d1 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_base.hpp
@@ -127,10 +127,10 @@ void G0InterpolationBase<Parameters>::initialize_linear_coefficients(Data& data)
       for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) {
         for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) {
           G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind) =
-              data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+              data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
           grad_G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind) =
-              (data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) -
-               data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind));
+              (data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) -
+               data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind));
         }
       }
     }
@@ -141,10 +141,10 @@ void G0InterpolationBase<Parameters>::initialize_linear_coefficients(Data& data)
       for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) {
         for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) {
           G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind - 1) =
-              data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+              data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
           grad_G0_r_t_shifted(nu0_ind, nu1_ind, r_ind, t_ind - 1) =
-              (data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) -
-               data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind));
+              (data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind + 1) -
+               data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind));
         }
       }
     }
@@ -169,7 +169,7 @@ void G0InterpolationBase<Parameters>::initialize_akima_coefficients(Data& data)
       for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) {
         for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) {
           for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++)
-            y[t_ind] = data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+            y[t_ind] = data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
 
           ai_obj.initialize(x, y);
 
@@ -187,7 +187,7 @@ void G0InterpolationBase<Parameters>::initialize_akima_coefficients(Data& data)
         for (int nu0_ind = 0; nu0_ind < b::dmn_size() * s::dmn_size(); nu0_ind++) {
           for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size(); t_ind++)
             y[t_ind - t::dmn_size() / 2] =
-                data.G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
+                data.mutable_G0_r_t_cluster_excluded(nu0_ind, nu1_ind, r_ind, t_ind);
 
           ai_obj.initialize(x, y);
 
diff --git a/include/dca/phys/dca_step/cluster_solver/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<device_t, Parameters> {
 
   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 <class configuration_type>
   void build_G_matrix(configuration_type& full_configuration,
                       const dca::linalg::Matrix<Scalar, device_t>& N,
@@ -174,7 +181,8 @@ void G_TOOLS<device_t, Parameters>::build_G_matrix(configuration_type& full_conf
                                                  G.ptr(0, 0), LD_G, thread_id, stream_id);
 
     if constexpr (dca::util::IsComplex_t<Scalar>::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..295bfe138 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp
@@ -41,7 +41,7 @@
 #include "dca/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu.hpp"
 #endif
 
-//#define DEBUG_SUBMATRIX
+// #define DEBUG_SUBMATRIX
 
 namespace dca {
 namespace phys {
@@ -330,7 +330,7 @@ void CtintWalkerBase<Parameters, DIST>::setInteractionVertices(const Data& data,
   vertices_.initialize(parameters.getDoubleUpdateProbability(), parameters.getAllSitesPartnership());
   vertices_.initializeFromHamiltonian(data.H_interactions);
   if (data.has_non_density_interactions()) {
-    vertices_.checkForInterbandPropagators(data.G0_r_t_cluster_excluded);
+    vertices_.checkForInterbandPropagators(data.mutable_G0_r_t_cluster_excluded);
     vertices_.initializeFromNonDensityHamiltonian(data.get_non_density_interactions());
   }
 }
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 <int dimension, typename Scalar>
 auto shrinkG0(const SpGreensFunction<dimension, Scalar>& G0) {
   func::function<Scalar, func::dmn_variadic<BDmn, BDmn, RDmn<dimension>, 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<dimension>::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..0ef835a1c 100644
--- a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp
@@ -116,6 +116,11 @@ class StdThreadQmciClusterSolver : public QmciSolver {
   void logSingleMeasurement(StdThreadAccumulatorType& accumulator, int stamping_period,
                             bool log_MFunction, bool log_MFunctionTime) const;
 
+  /// This collects measurements when dealing with disordered G0
+  void accumulateGkw(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 +128,16 @@ class StdThreadQmciClusterSolver : public QmciSolver {
 
   void initializeAndWarmUp(Walker& walker, int id, int walker_id);
 
-  void readConfigurations();
-  void writeConfigurations() const;
+  /** readConfigurations needs to be aware that there are different
+   *  vertex configurations for difference disorder configurations
+   *  but this is perhaps not very useful now since we generate new
+   *  disorder configurations for each iteration.
+   */
+  void readConfigurations(int disorder_configuiration = -1);
+  /** writeConfigurations needs to be aware that there are different
+   *  vertex configurations for difference disorder configurations
+   */
+  void writeConfigurations(int disorder_configuration = -1) const;
 
   void iterateOverLocalMeasurements(int walker_id, std::function<void(int, int, bool)>&& f);
 
@@ -182,7 +195,7 @@ StdThreadQmciClusterSolver<QmciSolver>::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(
@@ -261,7 +274,6 @@ void StdThreadQmciClusterSolver<QmciSolver>::integrate() {
 
     dca::profiling::Duration duration(end_time, start_time);
     total_time_ = duration.sec + 1.e-6 * duration.usec;
-
     printIntegrationMetadata();
   };
 
@@ -276,7 +288,7 @@ void StdThreadQmciClusterSolver<QmciSolver>::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 +359,8 @@ double StdThreadQmciClusterSolver<QmciSolver>::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 +379,8 @@ void StdThreadQmciClusterSolver<QmciSolver>::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<std::exception> exception_ptr;
 
@@ -499,6 +511,16 @@ auto StdThreadQmciClusterSolver<QmciSolver>::computeSingleMeasurement_G_k_w(
   return G_k_w;
 }
 
+template <class QmciSolver>
+void StdThreadQmciClusterSolver<QmciSolver>::accumulateGkw(double weight) {
+  QmciSolver::accumulateGkw(weight);
+}
+
+template <class QmciSolver>
+void StdThreadQmciClusterSolver<QmciSolver>::collectSingle() {
+  QmciSolver::collectSingle();
+}
+
 template <class QmciSolver>
 void StdThreadQmciClusterSolver<QmciSolver>::logSingleMeasurement(
     StdThreadAccumulatorType& accumulator_obj, int stamping_period, bool log_MFunction,
@@ -589,8 +611,8 @@ void StdThreadQmciClusterSolver<QmciSolver>::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 +663,7 @@ void StdThreadQmciClusterSolver<QmciSolver>::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 +700,16 @@ void StdThreadQmciClusterSolver<QmciSolver>::finalizeWalker(Walker& walker, int
 }
 
 template <class QmciSolver>
-void StdThreadQmciClusterSolver<QmciSolver>::writeConfigurations() const {
+void StdThreadQmciClusterSolver<QmciSolver>::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 +721,19 @@ void StdThreadQmciClusterSolver<QmciSolver>::writeConfigurations() const {
 }
 
 template <class QmciSolver>
-void StdThreadQmciClusterSolver<QmciSolver>::readConfigurations() {
+void StdThreadQmciClusterSolver<QmciSolver>::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<scalar_type>& 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<scalar_type> 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 <class Parameters>
+struct DcaDomainAliases {
+  using Real = typename Parameters::Real;
+  using Scalar = typename Parameters::Scalar;
+  using TpAccumulatorPrec = typename Parameters::TPAccumPrec;
+  using TpComplex = std::complex<TpAccumulatorPrec>;
+  using TDmn = func::dmn_0<domains::time_domain>;
+  using WDmn = func::dmn_0<domains::frequency_domain>;
+  using WVertexDmn = func::dmn_0<domains::vertex_frequency_domain<domains::COMPACT>>;
+  using WExchangeDmn = func::dmn_0<domains::FrequencyExchangeDomain>;
+
+  using BDmn = func::dmn_0<domains::electron_band_domain>;
+  using SDmn = func::dmn_0<domains::electron_spin_domain>;
+  using NuDmn = func::dmn_variadic<BDmn, SDmn>;  // orbital-spin index
+  using NuNuDmn = func::dmn_variadic<NuDmn, NuDmn>;
+  static constexpr int Dimension = Parameters::lattice_type::DIMENSION;
+  using CDA = ClusterDomainAliases<Dimension>;
+  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<domains::MomentumExchangeDomain>;
+};
+
+}  // 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 <string>
 #include <vector>
+#include <stdexcept>
+#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 <typename Writer>
   static void write(Writer& writer);
 
   template <typename Parameters>
   static void initialize(const Parameters& parameters);
 
+  template <typename ss_type>
+  static void print(const electron_band_domain& ebd, ss_type& ss);
+
   // For testing purposes only.
-  static void setAVectors(const std::vector<std::vector<double>>& 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<std::vector<double>>& 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<element_type> elements_;
+  static std::vector<int> flavors_;
 };
 
 template <typename Writer>
@@ -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 <typename ss_type>
+void electron_band_domain::print(const electron_band_domain& ebd, ss_type& ss) {
+  ss << "\t electron orbitals/bands: " << ebd.get_elements().size() << " of "
+     << vectorToString(ebd.get_flavors()) << '\n';
+}
+
 }  // namespace domains
 }  // namespace phys
 }  // namespace dca
diff --git a/include/dca/phys/parameters/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 <stdexcept>
+
+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 <typename Concurrency>
+  int getBufferSize(const Concurrency& concurrency) const;
+  template <typename Concurrency>
+  void pack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position) const;
+  template <typename Concurrency>
+  void unpack(const Concurrency& concurrency, char* buffer, int buffer_size, int& position);
+
+  template <typename ReaderOrWriter>
+  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 <typename Concurrency>
+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 <typename Concurrency>
+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 <typename Concurrency>
+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 <typename ReaderOrWriter>
+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<dca::config::McOptions::MC_REAL,
+                         typename dca::util::ScalarSelect<dca::config::McOptions::MC_REAL,
+                                                          Model::lattice_type::complex_g0>::type>>;
+
+#endif
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 <typename Concurrency, typename Threading, typename Profiler, typename
           typename RandomNumberGenerator, ClusterSolverId solver_name, class NUMTRAITS>
 class Parameters : public AnalysisParameters,
                    public DcaParameters,
+                   public DisorderParameters,
                    public DomainsParameters,
                    public DoubleCountingParameters,
                    public EdSolverParameters,
@@ -165,10 +167,12 @@ template <class Parameters, typename = bool>
 struct CheckParametersNumericTypes : public std::false_type {};
 
 template <class Parameters>
-struct CheckParametersNumericTypes <Parameters,
-				    std::enable_if_t<std::is_same<typename Parameters::Scalar,
-								  typename dca::util::ScalarSelect<typename Parameters::Real, Parameters::complex_g0>::type>::value, bool>>
-    : public std::true_type {};
+struct CheckParametersNumericTypes<
+    Parameters,
+    std::enable_if_t<std::is_same<typename Parameters::Scalar,
+                                  typename dca::util::ScalarSelect<typename Parameters::Real,
+                                                                   Parameters::complex_g0>::type>::value,
+                     bool>> : public std::true_type {};
 
 template <typename Concurrency, typename Threading, typename Profiler, typename Model,
           typename RandomNumberGenerator, ClusterSolverId solver_name, typename NUMTRAITS>
@@ -176,6 +180,7 @@ Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solve
            NUMTRAITS>::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<Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solve
 #endif
 
   // check consistency between the value of the Parameters::complex_g0 and Parameters NUMTRAITS
-  static_assert(std::is_same_v<typename Parameters::Scalar,
-		typename dca::util::ScalarSelect<typename Parameters::Real, Parameters::complex_g0>::type>);
+  static_assert(
+      std::is_same_v<
+          typename Parameters::Scalar,
+          typename dca::util::ScalarSelect<typename Parameters::Real, Parameters::complex_g0>::type>);
 }
 
 template <typename Concurrency, typename Threading, typename Profiler, typename Model,
@@ -266,8 +273,7 @@ void Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator,
 template <typename Concurrency, typename Threading, typename Profiler, typename Model,
           typename RandomNumberGenerator, ClusterSolverId solver_name, typename NUMTRAITS>
 void Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator, solver_name,
-                NUMTRAITS>::broadcast()
-{
+                NUMTRAITS>::broadcast() {
   concurrency_.broadcast_object(*this);
 }
 
@@ -451,6 +457,10 @@ void Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator,
   AnalysisParameters::readWrite(reader_or_writer);
   DcaParameters::readWrite(reader_or_writer);
   DomainsParameters::readWrite(reader_or_writer);
+  // It makes sense for this to come after because validity of input
+  // for disorder is dependent on particulars of the domains. That is when
+  // we have some input validation someday
+  DisorderParameters::readWrite(reader_or_writer);
   DoubleCountingParameters::readWrite(reader_or_writer);
   EdSolverParameters::readWrite(reader_or_writer);
   FourPointParameters<Model::DIMENSION>::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 <class Parameters>
+struct DcaSharedTypes {
+  using DDA = DcaDomainAliases<Parameters>;
+  using Real = typename DDA::Real;
+  using DisorderConfiguration =
+      func::function<Real, func::dmn_variadic<typename DDA::NuDmn, typename DDA::RClusterDmn>>;
+};
+}  // 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 <cassert>
+#include <iostream>
+
+#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 <class Parameters, DistType DIST>
+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 <typename Scalar>
 void G0Interpolation<linalg::CPU, Scalar>::initialize(const FunctionProxy<Scalar, PTdmn>& 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/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 <stdexcept>
+
+namespace dca {
+namespace phys {
+namespace domains {
+// dca::phys::domains::
+
+std::vector<int> 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<dca::config::McOptions::MC_REAL,
+                         typename dca::util::ScalarSelect<dca::config::McOptions::MC_REAL,
+                                                          Model::lattice_type::complex_g0>::type>>;
diff --git a/test/unit/math/nfft/CMakeLists.txt b/test/unit/math/nfft/CMakeLists.txt
index 327923ba6..9ceee3f8f 100644
--- a/test/unit/math/nfft/CMakeLists.txt
+++ b/test/unit/math/nfft/CMakeLists.txt
@@ -12,4 +12,4 @@ dca_add_gtest(dnfft_1d_gpu_test
   GTEST_MAIN
   CUDA
   INCLUDE_DIRS ${PROJECT_SOURCE_DIR}
-  LIBS FFTW::Double time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS})
+  LIBS FFTW::Double quantum_domains time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS})
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt
index 9e121d50c..1282a6b47 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/CMakeLists.txt
@@ -4,22 +4,22 @@ dca_add_gtest(sp_accumulator_gpu_test
   CUDA
   GTEST_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS     FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS}
+  LIBS     FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS}
   )
 
 dca_add_gtest(sp_accumulator_complex_gpu_test
   CUDA
   FAST
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS     FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS}
+  LIBS     FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS}
   )
 
-  
+
 dca_add_gtest(sp_accumulator_single_meas_G_test
   CUDA
   GTEST_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS     FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS}
+  LIBS     FFTW::Double quantum_domains ${DCA_LIBS} ${DCA_KERNEL_LIBS}
   )
 
 #;mc_kernels
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt
index df410d593..94d2ac38f 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/ndft/CMakeLists.txt
@@ -18,4 +18,4 @@ dca_add_gtest(cached_ndft_gpu_test
   CUDA
   FAST
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS function_transform ${DCA_KERNEL_LIBS} hdf5::hdf5_cpp magma::magma magma::sparse ${THIS_TEST_LIBS})
+  LIBS function_transform quantum_domains ${DCA_KERNEL_LIBS} hdf5::hdf5_cpp magma::magma magma::sparse ${THIS_TEST_LIBS})
diff --git a/test/unit/phys/parameters/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
+  }
+}