diff --git a/.gitignore b/.gitignore
index 91e050f0e..3ed603dee 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,6 @@ bin*
 *#
 *.swp
 .idea
+
+#Install directory
+install*
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/cmake/dca_config.cmake b/cmake/dca_config.cmake
index 976f24253..94ad1b81c 100644
--- a/cmake/dca_config.cmake
+++ b/cmake/dca_config.cmake
@@ -450,6 +450,17 @@ if(DCA_WITH_CTAUX_TRACING)
   add_compile_definitions(CTAUX_DEBUG_TRACING)
 endif()
 
+option(DCA_WITH_DISORDER "Enable the site-local disorder potential in the CT-AUX cluster solver." OFF)
+if(DCA_WITH_DISORDER)
+  # GPU disorder path is not implemented yet (Step G); the two-site walker G0 breaks the GPU build.
+  if(DCA_WITH_CUDA OR DCA_WITH_HIP)
+    message(FATAL_ERROR
+      "DCA_WITH_DISORDER is not supported together with a GPU build (DCA_WITH_CUDA / "
+      "DCA_WITH_HIP) yet. Reconfigure with the GPU backend disabled.")
+  endif()
+  add_compile_definitions(DISORDERED_G0)
+endif()
+
 configure_file("${PROJECT_SOURCE_DIR}/include/dca/config/mc_options.hpp.in"
         "${CMAKE_BINARY_DIR}/include/dca/config/mc_options.hpp" @ONLY)
 
diff --git a/include/dca/function/function.hpp b/include/dca/function/function.hpp
index 359402f1e..5b40279c9 100644
--- a/include/dca/function/function.hpp
+++ b/include/dca/function/function.hpp
@@ -408,6 +408,21 @@ class function {
   // These are the linear start and end indexes with respect to the complete function.
   std::size_t start_;
   std::size_t end_;
+
+public:
+#ifdef DEBUG
+  // Variadic debug accessor: forwards all args to operator()
+  template <typename... Ts>
+  const scalartype& value_at_debug(Ts&&... indices) const {
+    return (*this)(std::forward<Ts>(indices)...);
+  }
+
+  // Non-const version (if needed for writing)
+  template <typename... Ts>
+  scalartype& value_at_debug(Ts&&... indices) {
+    return (*this)(std::forward<Ts>(indices)...);
+  }
+#endif
 };
 
 template <typename scalartype, class domain, DistType DT>
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/linalg/util/util_lapack.hpp b/include/dca/linalg/util/util_lapack.hpp
index 3a29639be..80df3679b 100644
--- a/include/dca/linalg/util/util_lapack.hpp
+++ b/include/dca/linalg/util/util_lapack.hpp
@@ -60,7 +60,7 @@ inline void checkLapackInfoInternal(int info, std::string function_name, std::st
 #define warnLapackInfo(info) \
   dca::linalg::lapack::util::warnLapackInfoInternal(info, __FUNCTION__, __FILE__, __LINE__)
 inline void warnLapackInfoInternal(int info, std::string function_name, std::string file_name,
-                                    int line) {
+                                   int line) {
   if (info < 0) {
     std::stringstream s;
     s << "Error in function: " << function_name << " (" << file_name << ":" << line << ")"
@@ -69,15 +69,16 @@ inline void warnLapackInfoInternal(int info, std::string function_name, std::str
 
     throw LapackException(s.str(), info);
   }
-  else if (info > 0) {
-    std::cout << "warning lapack info = " << info << " at " << file_name << ":" << line << '\n';
-  }
+#ifndef NDEBUG
+  // else if (info > 0) {
+  //   std::cout << "warning lapack info = " << info << " at " << file_name << ":" << line << '\n';
+  // }
+#endif
 }
 
-  
-}  // util
-}  // lapack
-}  // linalg
-}  // dca
+}  // namespace util
+}  // namespace lapack
+}  // namespace linalg
+}  // namespace dca
 
 #endif  // DCA_LINALG_UTIL_UTIL_LAPACK_HPP
diff --git a/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp b/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp
index ef86b651c..052a56ae3 100644
--- a/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp
+++ b/include/dca/math/function_transform/hermite_splines/hermite_cubic_spline.hpp
@@ -22,8 +22,7 @@ namespace transform {
 // dca::math::transform::
 
 // Empty template declaration.
-template <typename lh_dmn_type, typename rh_dmn_type, ELEMENT_SPACINGS ES, BOUNDARY_CONDITIONS BC,
-          int DIMENSION>
+template <typename lh_dmn_type, typename rh_dmn_type, ELEMENT_SPACINGS ES, BOUNDARY_CONDITIONS BC, int DIMENSION>
 struct hermite_cubic_spline {};
 
 // Template specialization for the equidistant periodic case.
@@ -117,8 +116,98 @@ typename hermite_cubic_spline<lh_dmn_type, rh_dmn_type, EQUIDISTANT, PERIODIC, D
   return result;
 }
 
-}  // transform
-}  // math
-}  // dca
+template <typename lh_dmn_type, typename rh_dmn_type, int DIMENSION>
+class hermite_cubic_spline<lh_dmn_type, rh_dmn_type, EQUIDISTANT, INTERVAL, DIMENSION> {
+private:
+  typedef typename lh_dmn_type::dmn_specifications_type lh_spec_dmn_type;
+  typedef typename rh_dmn_type::dmn_specifications_type rh_spec_dmn_type;
+
+  typedef typename lh_spec_dmn_type::scalar_type lh_scalar_type;
+  typedef typename rh_spec_dmn_type::scalar_type rh_scalar_type;
+
+  typedef typename lh_spec_dmn_type::element_type lh_element_type;
+  typedef typename rh_spec_dmn_type::element_type rh_element_type;
+
+  typedef lh_scalar_type f_scalar_type;
+
+public:
+  static f_scalar_type execute(int i, int j);
+};
+
+template <typename lh_dmn_type, typename rh_dmn_type, int DIMENSION>
+typename hermite_cubic_spline<lh_dmn_type, rh_dmn_type, EQUIDISTANT, INTERVAL, DIMENSION>::f_scalar_type hermite_cubic_spline<
+    lh_dmn_type, rh_dmn_type, EQUIDISTANT, INTERVAL, DIMENSION>::execute(int i, int j) {
+  const static rh_scalar_type a = -0.5;
+
+  lh_element_type x = lh_dmn_type::get_elements()[i];
+  rh_element_type y = rh_dmn_type::get_elements()[j];
+
+  rh_scalar_type* super_basis = rh_dmn_type::get_super_basis();
+
+  rh_scalar_type* inv_basis = rh_dmn_type::get_inverse_basis();
+  rh_scalar_type* inv_super_basis = rh_dmn_type::get_inverse_super_basis();
+
+  rh_element_type delta(DIMENSION, 0.);
+  rh_element_type delta_affine(DIMENSION, 0.);
+
+  f_scalar_type result = 0;
+
+  {
+    for (int li = 0; li < DIMENSION; li++)
+      delta[li] = (y[li] - x[li]);
+
+    {
+      for (int li = 0; li < DIMENSION; li++)
+        delta_affine[li] = 0.;
+
+      for (int li = 0; li < DIMENSION; li++)
+        for (int lj = 0; lj < DIMENSION; lj++)
+          delta_affine[li] += inv_super_basis[li + lj * DIMENSION] * delta[lj];
+    }
+
+    for (int li = 0; li < DIMENSION; li++) {
+      while (delta_affine[li] > 0.5 - 1.e-6)
+        delta_affine[li] -= 1.;
+
+      while (delta_affine[li] < -0.5 - 1.e-6)
+        delta_affine[li] += 1.;
+    }
+
+    {
+      for (int li = 0; li < DIMENSION; li++)
+        delta[li] = 0.;
+
+      for (int li = 0; li < DIMENSION; li++)
+        for (int lj = 0; lj < DIMENSION; lj++)
+          delta[li] += super_basis[li + lj * DIMENSION] * delta_affine[lj];
+    }
+
+    {
+      for (int li = 0; li < DIMENSION; li++)
+        delta_affine[li] = 0.;
+
+      for (int li = 0; li < DIMENSION; li++)
+        for (int lj = 0; lj < DIMENSION; lj++)
+          delta_affine[li] += inv_basis[li + lj * DIMENSION] * delta[lj];
+    }
+
+    for (int li = 0; li < DIMENSION; li++)
+      delta[li] = (delta_affine[li] > -2. and delta_affine[li] < 2.)
+                      ? hermite_spline::cubic(0., delta_affine[li], 1., a)
+                      : 0.;
+
+    f_scalar_type t_result = 1;
+    for (int li = 0; li < DIMENSION; li++)
+      t_result *= delta[li];
+
+    result += t_result;
+  }
+
+  return result;
+}
+
+}  // namespace transform
+}  // namespace math
+}  // namespace dca
 
 #endif  // DCA_MATH_FUNCTION_TRANSFORM_HERMITE_SPLINES_HERMITE_CUBIC_SPLINE_HPP
diff --git a/include/dca/math/util/vector_operations.hpp b/include/dca/math/util/vector_operations.hpp
index cd8bf7923..638622b3f 100644
--- a/include/dca/math/util/vector_operations.hpp
+++ b/include/dca/math/util/vector_operations.hpp
@@ -104,16 +104,49 @@ std::complex<T> innerProduct(const std::vector<std::complex<T>>& x,
   return res;
 }
 
+template <typename T>
+std::complex<T> innerProduct(const std::vector<T>& x, const std::vector<std::complex<T>>& y) {
+  assert(x.size() == y.size());
+
+  std::complex<T> res(0.);
+  for (std::size_t i = 0; i < x.size(); ++i)
+    res += x[i] * std::conj(y[i]);
+
+  return res;
+}
+
+template <typename T>
+std::complex<T> innerProduct(const std::vector<std::complex<T>>& x, const std::vector<T>& y) {
+  assert(x.size() == y.size());
+
+  std::complex<T> res(0.);
+  for (std::size_t i = 0; i < x.size(); ++i)
+    res += x[i] * y[i];
+
+  return res;
+}
+
 // Treats scalars as vectors of size 1.
 template <typename T>
 T innerProduct(const T x, const T y) {
   return x * y;
 }
+
 template <typename T>
 std::complex<T> innerProduct(const std::complex<T> x, const std::complex<T> y) {
   return x * std::conj(y);
 }
 
+// template <typename T>
+// std::complex<T> innerProduct(const T x, const std::complex<T> y) {
+//   return x * std::conj(y);
+// }
+
+// template <typename T>
+// std::complex<T> innerProduct(const std::complex<T> x, const T y) {
+//   return x * y;
+// }
+
 // Computes the square of the L^2 norm of the vector x.
 template <typename T>
 auto l2Norm2(const std::vector<T>& x) {
diff --git a/include/dca/parallel/no_concurrency/serial_processor_grouping.hpp b/include/dca/parallel/no_concurrency/serial_processor_grouping.hpp
index e15e0e0ab..16acf73ef 100644
--- a/include/dca/parallel/no_concurrency/serial_processor_grouping.hpp
+++ b/include/dca/parallel/no_concurrency/serial_processor_grouping.hpp
@@ -38,6 +38,7 @@ class SerialProcessorGrouping {
   int last() const {
     return nr_threads_ - 1;
   }
+  void barrier() const {}
 
   auto get() const {
 #if defined (CRAY_MPICH_VERSION) && defined (DCA_HAVE_MPI)
diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp
index d1fafc22b..1a88f3439 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"
@@ -42,6 +44,7 @@
 #include "dca/phys/domains/cluster/cluster_operations.hpp"
 #include "dca/phys/domains/cluster/interpolation/hspline_interpolation/hspline_interpolation.hpp"
 #include "dca/phys/domains/cluster/momentum_exchange_domain.hpp"
+#include "dca/phys/domains/cluster/disorder_configuration_domain.hpp"
 #include "dca/phys/domains/quantum/brillouin_zone_cut_domain.hpp"
 #include "dca/phys/domains/quantum/electron_band_domain.hpp"
 #include "dca/phys/domains/quantum/electron_spin_domain.hpp"
@@ -56,6 +59,8 @@
 #include "dca/util/timer.hpp"
 #include "dca/util/to_string.hpp"
 #include "dca/distribution/dist_types.hpp"
+#include "dca/phys/types/dca_shared_types.hpp"
+#include "dca/util/type_help.hpp"
 #ifdef DCA_WITH_ADIOS2
 #include "dca/io/adios2/adios2_writer.hpp"
 #endif
@@ -102,8 +107,12 @@ class DcaData {
 
   using SpGreensFunction =
       func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KClusterDmn, WDmn>>;
+  using SpDisorderedGreensFunction =
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, KClusterDmn, NuDmn, KClusterDmn, WDmn>>;
   using SpRGreensFunction =
       func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RClusterDmn, WDmn>>;
+  using SpRDisorderedGreensFunction =
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RClusterDmn, NuDmn, RClusterDmn, WDmn>>;
 
   using TpGreensFunction =
       func::function<TpComplex,
@@ -111,6 +120,23 @@ 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>>;
+
+  // Two real-index (R1,R2) disordered G0 types. The imaginary-time version follows the Scalar
+  // convention of the other real-space-time G0s, while the Matsubara-frequency version is complex.
+  using Disordered_G0_r_r_t =
+      func::function<Scalar, func::dmn_variadic<NuDmn, RClusterDmn, NuDmn, RClusterDmn, TDmn>>;
+  using Disordered_G0_r_r_w =
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RClusterDmn, NuDmn, RClusterDmn, WDmn>>;
+
+  // 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.
@@ -125,6 +151,9 @@ class DcaData {
   template <typename Writer>
   void write(Writer& writer);
 
+  template <typename Writer>
+  void writeDisAvgGreensFunction(Writer& writer);
+
 #ifdef DCA_WITH_ADIOS2
   void writeDistributedG4Adios(io::ADIOS2Writer<Concurrency>& writer);
 #endif
@@ -186,18 +215,42 @@ 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;
+
+  /// These G0 are only used when diagonal disorder is applied
+  Disordered_G0_r_r_t disordered_G0_r_r_t_cl_exl;
+  /// Disordered G0 in r-space Matsubara frequency, computed directly in makeDisorderedG0 by
+  /// per-site frequency inversion. Feeds the r-space Dyson equation in accumulateGrrwFromMrrw.
+  Disordered_G0_r_r_w disordered_G0_r_r_w_cl_exl;
+  /// Disorder-averaged interacting Green's function, accumulated over configs in
+  /// accumulateGrrwFromMrrw.
+  Disordered_G0_r_r_w disorder_G_r_r_w;
+  /// Rank-local (not MPI-collected) analog of disorder_G_r_r_w: each rank sums its own per-config
+  /// Dyson here so computeErrorBars can take the cross-rank stddev (QMC error of the disorder
+  /// average). Only populated when error bars are requested.
+  Disordered_G0_r_r_w disorder_G_r_r_w_local;
+
+  /// Disorder realizations of the current DCA step: each configuration is the per-site random
+  /// potential (box in [-W/2,W/2], binary +/-V/2), with uniform averaging weights. Generated in
+  /// the DCA loop; written to output when the "dump-disorder-configs" parameter is set.
+  std::vector<DisorderConfiguration> disorder_configurations;
+  std::vector<double> disorder_weights;
 
   func::function<Real, NuDmn> orbital_occupancy;
 
@@ -261,7 +314,14 @@ class DcaData {
     return (bool)non_density_interactions_;
   }
 
+  void makeDisorderedG0(const DisorderConfiguration& disorder_configuration);
+
 private:  // Optional members.
+  void makeDisorderedG0(const DisorderConfiguration& disorder_configuration,
+                        const Type_G0_r_t& g0_r_t_cl_exl);
+
+  /// Due to the new ability to modify the G0 for disorder this is the
+  /// immutable g0 from the last iteration.
   std::unique_ptr<SpGreensFunction> G_k_w_err_;
   std::unique_ptr<SpRGreensFunction> G_r_w_err_;
   std::unique_ptr<SpGreensFunction> Sigma_err_;
@@ -365,19 +425,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
-  }
+    }
   }
 }
 
@@ -530,9 +590,39 @@ void DcaData<Parameters, DT>::write(Writer& writer) {
     }
 #endif
   }
+
+  // Dump the disorder realizations of this DCA step (the per-site random potentials) as a
+  // single function over <nu, r, configuration index>. The disorder potential is site-local
+  // (one value per site), so a single R index is the faithful representation.
+  if (parameters_.dump_disorder_configs() && !disorder_configurations.empty()) {
+    domains::DisorderConfigurationDomain::initialize(static_cast<int>(disorder_configurations.size()));
+    using ConfigDmn = func::dmn_0<domains::DisorderConfigurationDomain>;
+    func::function<Real, func::dmn_variadic<NuDmn, RClusterDmn, ConfigDmn>> disorder_randoms(
+        "disorder_configurations");
+
+    // Each configuration is a contiguous nu*r block; copy it into the ic-th slot (the
+    // configuration index is the slowest-varying dimension of disorder_randoms).
+    const int block = disorder_configurations[0].size();
+    for (int ic = 0; ic < static_cast<int>(disorder_configurations.size()); ++ic) {
+      const Real* src = disorder_configurations[ic].values();
+      Real* dst = disorder_randoms.values() + ic * block;
+      for (int j = 0; j < block; ++j)
+        dst[j] = src[j];
+    }
+    writer.execute(disorder_randoms);
+  }
+
   writer.close_group();
 }
 
+template <class Parameters, DistType DT>
+template <typename Writer>
+void DcaData<Parameters, DT>::writeDisAvgGreensFunction(Writer& writer) {
+  if (parameters_.dump_cluster_Greens_functions()) {
+    writer.execute(G_k_w);
+  }
+}
+
 template <class Parameters, DistType DT>
 void DcaData<Parameters, DT>::initialize() {
   initializeH0_and_H_i();
@@ -555,7 +645,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 +657,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);
 }
 
@@ -659,13 +749,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 +787,75 @@ 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) {
+  const int nu_matrix_dim = NuDmn::dmn_size();
+  const int r_matrix_dim = RClusterDmn::dmn_size();
+  const int nu_r_matrix_dim = nu_matrix_dim * r_matrix_dim;
+
+  // The frequency clean G0 is kept in k-space; bring it to (r,w).
+  math::transform::FunctionTransform<KClusterDmn, RClusterDmn>::execute(G0_k_w_cluster_excluded,
+                                                                        G0_r_w_cluster_excluded);
+
+  // Unfold the single-displacement clean G0 to two indices: G0(R1,R2) = g(subtract(R2,R1)).
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int R1 = 0; R1 < r_matrix_dim; ++R1)
+      for (int R2 = 0; R2 < r_matrix_dim; ++R2) {
+        const int d = RClusterDmn::parameter_type::subtract(R2, R1);  // index of r_{R1} - r_{R2}
+        for (int inu1 = 0; inu1 < nu_matrix_dim; ++inu1)
+          for (int inu2 = 0; inu2 < nu_matrix_dim; ++inu2)
+            disordered_G0_r_r_w_cl_exl(inu1, R1, inu2, R2, w) =
+                G0_r_w_cluster_excluded(inu1, inu2, d, w);
+      }
+
+  // Snapshot the clean unfold as the w->t tail-subtraction reference (shares G0_dis's 1/(iw) moment).
+  Disordered_G0_r_r_w clean_G0_r_r_w(disordered_G0_r_r_w_cl_exl);
+
+  // Per-frequency disorder Dyson: invert, add static V(R) on the diagonal, invert back.
+  dca::linalg::Matrix<std::complex<Real>, dca::linalg::CPU> g0_dis_block(nu_r_matrix_dim);
+  for (int w = 0; w < WDmn::dmn_size(); ++w) {
+    dca::linalg::matrixop::copyArrayToMatrix(nu_r_matrix_dim, nu_r_matrix_dim,
+                                             &disordered_G0_r_r_w_cl_exl(0, 0, 0, 0, w),
+                                             nu_r_matrix_dim, g0_dis_block);
+    dca::linalg::matrixop::inverse(g0_dis_block);
+    for (int ir = 0; ir < r_matrix_dim; ++ir)
+      for (int imd = 0; imd < nu_matrix_dim; ++imd)
+        g0_dis_block(imd + nu_matrix_dim * ir, imd + nu_matrix_dim * ir) +=
+            disorder_configuration(imd, ir);
+    dca::linalg::matrixop::inverse(g0_dis_block);
+    dca::linalg::matrixop::copyMatrixToArray(g0_dis_block,
+                                             &disordered_G0_r_r_w_cl_exl(0, 0, 0, 0, w),
+                                             nu_r_matrix_dim);
+  }
+
+  // Disordered G0 in imaginary time: the Dyson is a convolution in tau, so transform w->t with the
+  // tail trick instead of inverting per slice. Subtract the clean reference (O(1/w^2) remainder).
+  Disordered_G0_r_r_w delta_G0_r_r_w("delta_disordered_G0_r_r_w");
+  for (int i = 0; i < delta_G0_r_r_w.size(); ++i)
+    delta_G0_r_r_w(i) = disordered_G0_r_r_w_cl_exl(i) - clean_G0_r_r_w(i);
+
+  // FT the smooth remainder w->t (narrows complex->Scalar; t is trailing, (nu,R,nu,R) spectate).
+  Disordered_G0_r_r_t delta_G0_r_r_t("delta_disordered_G0_r_r_t");
+  math::transform::FunctionTransform<WDmn, TDmn>::execute(delta_G0_r_r_w, delta_G0_r_r_t);
+
+  // Add back the exact clean tau reference (V=0 => remainder is zero => clean tau G0 exactly).
+  for (int it = 0; it < TDmn::dmn_size(); ++it)
+    for (int R1 = 0; R1 < r_matrix_dim; ++R1)
+      for (int R2 = 0; R2 < r_matrix_dim; ++R2) {
+        const int d = RClusterDmn::parameter_type::subtract(R2, R1);
+        for (int inu1 = 0; inu1 < nu_matrix_dim; ++inu1)
+          for (int inu2 = 0; inu2 < nu_matrix_dim; ++inu2)
+            disordered_G0_r_r_t_cl_exl(inu1, R1, inu2, R2, it) =
+                delta_G0_r_r_t(inu1, R1, inu2, R2, it) + g0_r_t_cl_exl(inu1, inu2, d, it);
+      }
+}
+
+template <class Parameters, DistType DT>
+void DcaData<Parameters, DT>::makeDisorderedG0(const DisorderConfiguration& disorder_configuration) {
+  makeDisorderedG0(disorder_configuration, 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_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp
index bef3c0658..ead7a6a95 100644
--- a/include/dca/phys/dca_loop/dca_loop.hpp
+++ b/include/dca/phys/dca_loop/dca_loop.hpp
@@ -25,6 +25,7 @@
 #include "dca/io/filesystem.hpp"
 #include "dca/io/writer.hpp"
 #include "dca/io/io_types.hpp"
+#include "dca/math/function_transform/function_transform.hpp"
 #include "dca/phys/dca_algorithms/compute_greens_function.hpp"
 #include "dca/phys/dca_data/dca_data.hpp"
 #include "dca/phys/dca_loop/dca_loop_data.hpp"
@@ -41,6 +42,9 @@
 #include "dca/util/print_time.hpp"
 #include "dca/util/signal_handler.hpp"
 
+#include "dca/phys/dca_loop/disorder/disorder_average.hpp"
+#include "dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp"
+
 namespace dca {
 namespace phys {
 // dca::phys::
@@ -94,19 +98,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 +251,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 +298,62 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::execute() {
   }
 }
 
+template <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
+double DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::workTheClusters() {
+#ifdef DISORDERED_G0
+  // The ordered single-cluster path is disabled under DISORDERED_G0, so with zero configurations
+  // fall back to one V=0 configuration, driving the disorder pipeline through its clean limit.
+  // Silent for a missing disorder section (deliberate clean run); warn if a section requested 0.
+  if (parameters.get_disorder_num_configurations() == 0) {
+    if (parameters.get_disorder_present() && concurrency.id() == concurrency.first())
+      std::cerr << "Warning: a disorder section was provided but disorder-num-configurations is 0; "
+                   "running a single configuration with no disorder.\n";
+    parameters.set_disorder_potential(0.0);
+    parameters.set_disorder_num_configurations(1);
+  }
+#endif
+  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, MOMS.disorder_configurations, MOMS.disorder_weights);
+    int num_configurations = MOMS.disorder_configurations.size();
+
+    if (concurrency.id() == concurrency.first()) {
+      if (parameters.get_disorder_unique_configs() &&
+          parameters.get_disorder_distribution() == "box")
+        std::cerr << "Warning: Uniqueness of configurations is checked only for binary disorder\n";
+      if (num_configurations < parameters.get_disorder_num_configurations())
+        std::cerr << "Warning: Cluster too small to generate the requested "
+                  << parameters.get_disorder_num_configurations()
+                  << " unique disorder configurations; generated only " << num_configurations
+                  << ".\n";
+    }
+#ifdef DISORDERED_G0
+    MOMS.disorder_G_r_r_w = 0;
+    MOMS.disorder_G_r_r_w_local = 0;
+    MOMS.G_k_w = 0;
+#else
+    MOMS.accumulated_G_k_w = 0;
+#endif
+    for (int id = 0; id < num_configurations; ++id) {
+      MOMS.makeDisorderedG0(MOMS.disorder_configurations[id]);
+      if (concurrency.id() == concurrency.first())
+        std::cout << "Solving disorder configuration " << id << '\n';
+      solve_cluster_problem(dca_iteration_);
+      monte_carlo_integrator_.accumulateGkw(MOMS.disorder_weights[id]);
+    }
+    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();
@@ -353,7 +426,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::adjust_impuri
 template <typename ParametersType, typename DcaDataType, typename MCIntegratorType, DistType DIST>
 void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::perform_cluster_exclusion_step() {
   if (concurrency.id() == concurrency.first())
-    std::cout << "\n\t\t cluster-exclusion-step " << dca::util::print_time();
+    std::cout << "\n\t\t cluster-exclusion-step " << dca::util::print_time() << std::endl;
 
   profiler_type profiler("cluster-exclusion-step", "DCA", __LINE__);
 
@@ -361,7 +434,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 +450,55 @@ 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() {
+  auto& disorder_weights = MOMS.disorder_weights;
+  const double total_disorder_weight =
+      std::accumulate(disorder_weights.begin(), disorder_weights.end(), 0.0);
+
+#ifdef DISORDERED_G0
+  using RClusterDmn = typename DcaDataType::RClusterDmn;
+  using KClusterDmn = typename DcaDataType::KClusterDmn;
+  using NuDmn = typename DcaDataType::NuDmn;
+  using WDmn = typename DcaDataType::WDmn;
+
+  // Disorder average (weighted mean) of the two-site interacting G.
+  MOMS.disorder_G_r_r_w /= total_disorder_weight;
+
+  // Translational average <nu1,R1,nu2,R2,w> -> <nu1,nu2,R,w> (restores translational symmetry).
+  disorder::translationalAverage<NuDmn, RClusterDmn, WDmn>(MOMS.disorder_G_r_r_w, MOMS.G_r_w);
+
+  // r -> k (cluster-momentum-diagonal), the same transform as the clean DCA path.
+  math::transform::FunctionTransform<RClusterDmn, KClusterDmn>::execute(MOMS.G_r_w, MOMS.G_k_w);
+#else
+  MOMS.G_k_w = MOMS.accumulated_G_k_w;
+  MOMS.G_k_w /= total_disorder_weight;
+#endif
 
-    return L2_Sigma_difference;
-  }
+  if (concurrency.id() == concurrency.first())
+    std::cout << " Averaged over " << disorder_weights.size()
+              << " disorder configurations with weight " << total_disorder_weight << '\n';
+}
+
+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__);
+
+  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..de1ed13cb 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,8 @@ 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;
   int last_completed_iteration = -1;
 };
 
@@ -113,7 +124,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/disorder_average.hpp b/include/dca/phys/dca_loop/disorder/disorder_average.hpp
new file mode 100644
index 000000000..dc86ccaf2
--- /dev/null
+++ b/include/dca/phys/dca_loop/disorder/disorder_average.hpp
@@ -0,0 +1,78 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Numeric kernels of the DISORDERED_G0 disorder-averaging path, factored out of the cluster
+// solver and the DCA loop so they can be unit tested in isolation:
+//
+//   accumulateDisorderDyson - per-frequency dense (nu*N_R) Dyson G = G0_dis - G0_dis*M*G0_dis/beta,
+//                             accumulated weighted into a two space-index function. Used by
+//                             CtauxClusterSolver::accumulateGrrwFromMrrw.
+//   translationalAverage    - project the disorder-averaged two-site G(<nu,R,nu,R,w>) onto the
+//                             single-displacement g(<nu,nu,R,w>) with the cluster subtract table.
+//                             Used by DcaLoop::averageGkw.
+
+#ifndef DCA_PHYS_DCA_LOOP_DISORDER_DISORDER_AVERAGE_HPP
+#define DCA_PHYS_DCA_LOOP_DISORDER_DISORDER_AVERAGE_HPP
+
+#include <complex>
+
+#include "dca/linalg/matrix.hpp"
+#include "dca/linalg/matrix_view.hpp"
+#include "dca/linalg/matrixop.hpp"
+
+namespace dca::phys::disorder {
+
+// Per-frequency disorder Dyson on the dense two-site basis. At each w the (nu,R) indices form a
+// (nu*N_R) square block (row index nu + nu_size*R, leading dimension nu*N_R; the contiguous layout
+// of <nu,R,nu,R,w> functions). Computes G = G0_dis - G0_dis*M*G0_dis/beta and accumulates
+// weight*G into G_acc (G_acc is NOT zeroed - the caller accumulates over configurations).
+//
+// G0_dis and m are inputs (taken by non-const ref so MatrixView can alias their storage); G_acc is
+// the running accumulator. FG0/FM/FAcc are the (possibly distinct) two space-index function types.
+template <class NuDmn, class RDmn, class WDmn, typename Real, class FG0, class FM, class FAcc>
+void accumulateDisorderDyson(FG0& G0_dis, FM& m, double beta, double weight, FAcc& G_acc) {
+  const std::size_t matrix_size = NuDmn::dmn_size() * RDmn::dmn_size();
+  dca::linalg::Matrix<std::complex<Real>, dca::linalg::CPU> G0_times_M(matrix_size);
+  dca::linalg::Matrix<std::complex<Real>, dca::linalg::CPU> G0_M_G0(matrix_size);
+  using MatrixView = dca::linalg::MatrixView<std::complex<Real>, dca::linalg::CPU>;
+
+  for (int w = 0; w < WDmn::dmn_size(); ++w) {
+    const MatrixView G0_matrix(&G0_dis(0, 0, 0, 0, w), matrix_size);
+    const MatrixView M_matrix(&m(0, 0, 0, 0, w), matrix_size);
+
+    dca::linalg::matrixop::gemm(G0_matrix, M_matrix, G0_times_M);  // G0 * M
+    dca::linalg::matrixop::gemm(G0_times_M, G0_matrix, G0_M_G0);   // (G0 * M) * G0
+
+    MatrixView acc(&G_acc(0, 0, 0, 0, w), matrix_size);
+    for (int j = 0; j < static_cast<int>(matrix_size); ++j)
+      for (int i = 0; i < static_cast<int>(matrix_size); ++i)
+        acc(i, j) += static_cast<Real>(weight) * (G0_matrix(i, j) - G0_M_G0(i, j) / beta);
+  }
+}
+
+// Translational average of the disorder-averaged two-site G: project <nu1,R1,nu2,R2,w> onto the
+// single-displacement g <nu1,nu2,R,w>. Disorder averaging restores translational symmetry, so g
+// depends only on the displacement R. Summing over R1 = R0 at fixed displacement R uses
+// R2 = subtract(R,R0) (subtract(i,j) = index of r_j - r_i), the exact inverse of the unfold
+// g(R1-R2) -> G(R1,R2) that uses subtract(R2,R1). g_r is overwritten (zeroed then filled).
+template <class NuDmn, class RDmn, class WDmn, class TwoRFunction, class OneRFunction>
+void translationalAverage(const TwoRFunction& g_rr, OneRFunction& g_r) {
+  const int N = RDmn::dmn_size();
+  g_r = 0;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int R = 0; R < N; ++R)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          for (int R0 = 0; R0 < N; ++R0)
+            g_r(n1, n2, R, w) += g_rr(n1, R0, n2, RDmn::parameter_type::subtract(R, R0), w);
+  g_r /= static_cast<double>(N);
+}
+
+}  // namespace dca::phys::disorder
+
+#endif  // DCA_PHYS_DCA_LOOP_DISORDER_DISORDER_AVERAGE_HPP
diff --git a/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp b/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp
new file mode 100644
index 000000000..53840ab0a
--- /dev/null
+++ b/include/dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp
@@ -0,0 +1,135 @@
+// 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 <algorithm>
+#include <random>
+#include <string>
+#include <unordered_set>
+#include <vector>
+#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) {
+    const int m_configs = parameters.get_disorder_num_configurations();
+    assert(m_configs > 0);
+    disorder_configurations.resize(m_configs);
+    disorder_weights.resize(m_configs);
+
+    const int n_bands = BDmn::dmn_size();
+    const int n_rsites = RDmn::dmn_size();
+
+    const std::string& distribution = parameters.get_disorder_distribution();
+    const double potential = parameters.get_disorder_potential();
+    const double density = parameters.get_disorder_density();
+    const bool enforce_unique = parameters.get_disorder_unique_configs();
+
+    const bool is_box = (distribution == "box");
+
+    // Seed from the Monte Carlo seed, shifted by a fixed offset so the disorder realizations
+    // use a stream distinct from the QMC walkers (which derive their seeds from get_seed()).
+    // The seed is read on the root rank and broadcast as part of Parameters, so it is identical
+    // on every rank (every rank generates the same disorder ensemble) and fixed for the whole
+    // run (the ensemble is the same in every DCA iteration, i.e. quenched disorder). No separate
+    // seed broadcast is needed here.
+    constexpr unsigned disorder_seed_offset = 0x9e3779b9u;  // golden ratio; arbitrary but fixed
+    std::mt19937 gen(static_cast<unsigned>(parameters.get_seed()) + disorder_seed_offset);
+    // box:    site potential ~ Uniform[-W/2, W/2], with W = potential.
+    std::uniform_real_distribution<double> box_dist(-0.5 * potential, 0.5 * potential);
+    std::uniform_real_distribution<double> unit_dist(0.0, 1.0);
+
+    // Draw one disorder value for a single (band, site):
+    //  - box:    uniform in [-W/2, W/2].
+    //  - binary: +V/2 with probability `density`, else -V/2 (V = potential).
+    auto draw_site_value = [&]() -> double {
+      if (is_box)
+        return box_dist(gen);
+      return (unit_dist(gen) < density ? 1.0 : -1.0) * 0.5 * potential;
+    };
+
+    // Uniqueness only makes sense for the discrete binary distribution; a continuous
+    // box draw is distinct with probability one.
+    const bool tracking_unique = enforce_unique && !is_box;
+    std::unordered_set<std::vector<int>, detail::VectorHash> seen;
+
+    int generated = 0;
+    int attempts = 0;
+    const int max_attempts = 1000;
+    while (generated < m_configs && attempts < max_attempts) {
+      auto& config = disorder_configurations[generated];
+
+      std::vector<int> pattern;  // sign pattern, only populated for the uniqueness check
+      if (tracking_unique)
+        pattern.reserve(static_cast<std::size_t>(n_bands) * n_rsites);
+
+      // Fill the per-(band, site) random potential, spin-symmetrically.
+      for (int isite = 0; isite < n_rsites; ++isite)
+        for (int iband = 0; iband < n_bands; ++iband) {
+          const double value = draw_site_value();
+          config(iband, 0, isite) = value;
+          config(iband, 1, isite) = value;
+          if (tracking_unique)
+            pattern.push_back(value > 0.0 ? 1 : -1);
+        }
+
+      if (tracking_unique && !seen.insert(pattern).second) {
+        ++attempts;  // duplicate binary configuration, redraw
+        continue;
+      }
+      ++generated;
+      attempts = 0;
+    }
+    // If uniqueness was requested but the cluster is too small to supply m_configs
+    // distinct configurations, keep only the ones actually generated.
+    if (generated < m_configs) {
+      disorder_configurations.resize(generated);
+      disorder_weights.resize(generated);
+    }
+
+    // Equal weight for every configuration in the disorder average.
+    const double weight = 1.0 / static_cast<double>(disorder_configurations.size());
+    std::fill(disorder_weights.begin(), disorder_weights.end(), weight);
+  }
+};
+}  // namespace dca::phys
+
+#endif
diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp
index 535582b29..4b1eb9296 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp
@@ -41,7 +41,7 @@ namespace solver {
 namespace ctaux {
 // dca::phys::solver::ctaux::
 using dca::util::SignType;
-  
+
 template <class 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"),
 
diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp
index 343b3f37e..906af5d38 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp
@@ -70,7 +70,7 @@ class CtauxAccumulator : public MC_accumulator_data<typename Parameters::Scalar>
   using ParametersType = Parameters;
   using DataType = Data;
   using BaseClass = MC_accumulator_data<Scalar>;
-  
+
   typedef vertex_pair<Parameters> vertex_pair_type;
   typedef vertex_singleton vertex_singleton_type;
 
@@ -255,8 +255,8 @@ class CtauxAccumulator : public MC_accumulator_data<typename Parameters::Scalar>
  *  in the initializers.
  */
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
-CtauxAccumulator<device_t, Parameters, Data, DIST>::CtauxAccumulator(
-    const Parameters& parameters_ref, Data& data_ref, int id)
+CtauxAccumulator<device_t, Parameters, Data, DIST>::CtauxAccumulator(const Parameters& parameters_ref,
+                                                                     Data& data_ref, int id)
     : MC_accumulator_data<Scalar>(),
 
       parameters_(parameters_ref),
@@ -307,9 +307,8 @@ void CtauxAccumulator<device_t, Parameters, Data, DIST>::initialize(int dca_iter
                                       parameters_.dump_every_iteration());
 
   if (perform_equal_time_accumulation_) {
-    equal_time_accumulator_ptr_ =
-        std::make_unique<ctaux::TpEqualTimeAccumulator<Parameters, Data>>(parameters_, data_,
-                                                                                thread_id);
+    equal_time_accumulator_ptr_ = std::make_unique<ctaux::TpEqualTimeAccumulator<Parameters, Data>>(
+        parameters_, data_, thread_id);
     equal_time_accumulator_ptr_->resetAccumulation();
   }
 }
@@ -406,8 +405,8 @@ void CtauxAccumulator<device_t, Parameters, Data, DIST>::updateFrom(walker_type&
   walker.get_error_distribution() = 0;
 #endif  // DCA_WITH_QMC_BIT
 
-  //single_particle_accumulator_obj.syncStreams(*event);
-  //two_particle_accumulator_.syncStreams(*event);
+  // single_particle_accumulator_obj.syncStreams(*event);
+  // two_particle_accumulator_.syncStreams(*event);
 }
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
@@ -519,7 +518,8 @@ void CtauxAccumulator<device_t, Parameters, Data, DIST>::accumulate_equal_time_q
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
 void CtauxAccumulator<device_t, Parameters, Data, DIST>::accumulate_two_particle_quantities() {
   profiler_type profiler("tp-accumulation", "CT-AUX accumulator", __LINE__, thread_id);
-  gflop_ += 1e-9 * two_particle_accumulator_.accumulate(M_, hs_configuration_, current_phase_.getSign());
+  gflop_ +=
+      1e-9 * two_particle_accumulator_.accumulate(M_, hs_configuration_, current_phase_.getSign());
 }
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
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..acb2ddb44 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,15 +29,18 @@
 #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"
+#include "dca/phys/dca_loop/disorder/disorder_average.hpp"
 #include "dca/phys/dca_step/cluster_solver/ctaux/ctaux_accumulator.hpp"
 #include "dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp"
 #include "dca/phys/dca_step/cluster_solver/shared_tools/interpolation/g0_interpolation.hpp"
 #include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/time_correlator.hpp"
 #include "dca/phys/dca_step/symmetrization/symmetrize.hpp"
 #include "dca/phys/domains/cluster/cluster_domain.hpp"
+#include "dca/phys/domains/domain_aliases.hpp"
 #include "dca/phys/domains/quantum/electron_band_domain.hpp"
 #include "dca/phys/domains/quantum/electron_spin_domain.hpp"
 #include "dca/phys/domains/time_and_frequency/frequency_domain.hpp"
@@ -71,10 +74,14 @@ class CtauxClusterSolver {
   using Walker = ctaux::CtauxWalker<device_t, Parameters, Data>;
   using Accumulator = ctaux::CtauxAccumulator<device_t, Parameters, Data, DIST>;
   using SpGreensFunction = typename Data::SpGreensFunction;
+  using SpDisorderedGreensFunction = typename Data::SpDisorderedGreensFunction;
 
   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,9 +90,12 @@ 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>;
+  using NuRNuRClusterWDmn = func::dmn_variadic<nu, RDmn, nu, RDmn, w>;
 
 public:
   CtauxClusterSolver(Parameters& parameters_ref, Data& MOMS_ref,
@@ -117,7 +127,12 @@ class CtauxClusterSolver {
     return g0_;
   };
 
-  typename Walker::Resource& getResource() { return dummy_walker_resource_; };
+  typename Walker::Resource& getResource() {
+    return dummy_walker_resource_;
+  };
+
+  void accumulateGkw(double weight);
+  void collectSingle();
 
 protected:
   void warmUp(Walker& walker);
@@ -133,6 +148,19 @@ class CtauxClusterSolver {
   void collect_measurements();
 
   void compute_G_k_w_from_M_r_w();
+  // Now this needs to have k1, k2
+  void accumulateGkwFromMrw(
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, KDmn, WDmn>>& G_k_w,
+      double weight);
+  // r-space disorder analog of accumulateGkwFromMrw: per-config Dyson into the two-r accumulator.
+  void accumulateGrrwFromMrrw(
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>>& G_r_r_w,
+      double weight);
+  // Rank-local (un-collected) analog of accumulateGrrwFromMrrw for QMC error bars: Dysons this
+  // rank's own per-config M (normalized by the per-rank phase) into G_r_r_w.
+  void accumulateLocalGrrwFromMrrw(
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>>& G_r_r_w,
+      double weight);
 
   double compute_S_k_w_from_G_k_w();
 
@@ -164,6 +192,7 @@ class CtauxClusterSolver {
   G0Interpolation<device, Scalar> g0_;
 
   typename Walker::Resource dummy_walker_resource_;
+
 private:
   Rng rng_;
 
@@ -176,6 +205,10 @@ class CtauxClusterSolver {
   FPScalar accumulated_sign_;
   func::function<std::complex<Real>, NuNuRClusterWDmn> M_r_w_;
   func::function<std::complex<Real>, NuNuRClusterWDmn> M_r_w_squared_;
+#ifdef DISORDERED_G0
+  func::function<std::complex<Real>, NuRNuRClusterWDmn> M_r_r_w_;
+  func::function<std::complex<Real>, NuRNuRClusterWDmn> M_r_r_w_squared_;
+#endif
 
   bool averaged_;
   bool compute_jack_knife_;
@@ -207,8 +240,19 @@ CtauxClusterSolver<device_t, Parameters, Data, DIST>::CtauxClusterSolver(
 
       M_r_w_("M_r_w"),
       M_r_w_squared_("M_r_w_squared"),
+#ifdef DISORDERED_G0
+      M_r_r_w_("M_r_r_w"),
+      M_r_r_w_squared_("M_r_r_w_squared"),
+#endif
       averaged_(false),
       writer_(writer) {
+#ifdef DISORDERED_G0
+  // The two-particle (G4) accumulator measures against the clean G0_k_w_cluster_excluded, which is
+  // inconsistent with the disordered walker; fail fast rather than produce wrong G4 under disorder.
+  if (parameters_.isAccumulatingG4())
+    throw std::logic_error(
+        "Two-particle (G4) accumulation is not supported with DISORDERED_G0.");
+#endif
   if (concurrency_.id() == concurrency_.first())
     std::cout << "\n\n\t CT-AUX Integrator is born \n" << std::endl;
 }
@@ -289,15 +333,31 @@ 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();
+#ifdef DISORDERED_G0
+  accumulateGrrwFromMrrw(data_.disorder_G_r_r_w, weight);
+  // Per-rank (un-collected) accumulation for the cross-rank QMC error bar (only when requested).
+  if (accumulator_.compute_std_deviation())
+    accumulateLocalGrrwFromMrrw(data_.disorder_G_r_r_w_local, weight);
+#else
+  accumulateGkwFromMrw(data_.accumulated_G_k_w, weight);
+#endif
+}
+
+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);
 
@@ -412,13 +472,31 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::computeErrorBars() {
   if (concurrency_.id() == concurrency_.first())
     std::cout << "\n\t\t compute-error-bars on Self-energy\t" << dca::util::print_time() << "\n\n";
 
+  accumulator_.finalize();
+
+#ifdef DISORDERED_G0
+  // QMC error bars on the disorder-averaged G_k_w/Sigma. disorder_G_r_r_w_local holds this rank's
+  // own (un-collected) disorder-averaged two-site G; run averageGkw's pipeline (normalize ->
+  // translational average -> FT r->k) per rank, then take the cross-rank mean and stddev.
+  const double total_disorder_weight =
+      std::accumulate(data_.disorder_weights.begin(), data_.disorder_weights.end(), 0.0);
+  data_.disorder_G_r_r_w_local /= total_disorder_weight;
+
+  func::function<std::complex<Real>, NuNuRClusterWDmn> G_r_w_new("G_r_w_new");
+  disorder::translationalAverage<nu, RDmn, w>(data_.disorder_G_r_r_w_local, G_r_w_new);
+
   func::function<std::complex<Real>, func::dmn_variadic<nu, nu, KDmn, w>> G_k_w_new("G_k_w_new");
+  math::transform::FunctionTransform<RDmn, KDmn>::execute(G_r_w_new, G_k_w_new);
 
-  func::function<std::complex<Real>, func::dmn_variadic<nu, nu, RDmn, w>> M_r_w_new("M_r_w_new");
-  func::function<std::complex<Real>, func::dmn_variadic<nu, nu, KDmn, w>> M_k_w_new("M_k_w_new");
+  compute_S_k_w_new(G_k_w_new, Sigma_new_);
 
-  accumulator_.finalize();
+  concurrency_.average_and_compute_stddev(Sigma_new_, data_.get_Sigma_stdv());
+  concurrency_.average_and_compute_stddev(G_k_w_new, data_.get_G_k_w_stdv());
+#else
+  func::function<std::complex<Real>, func::dmn_variadic<nu, nu, KDmn, w>> G_k_w_new("G_k_w_new");
+  func::function<std::complex<Real>, func::dmn_variadic<nu, nu, KDmn, w>> M_k_w_new("M_k_w_new");
 
+  func::function<std::complex<Real>, func::dmn_variadic<nu, nu, RDmn, w>> M_r_w_new("M_r_w_new");
   M_r_w_new = accumulator_.get_sign_times_M_r_w();
   M_r_w_new /= static_cast<typename decltype(M_r_w_new)::this_scalar_type>(
       accumulator_.get_accumulated_phase());
@@ -430,6 +508,7 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::computeErrorBars() {
 
   concurrency_.average_and_compute_stddev(Sigma_new_, data_.get_Sigma_stdv());
   concurrency_.average_and_compute_stddev(G_k_w_new, data_.get_G_k_w_stdv());
+#endif  // DISORDERED_G0
 
   // sum G4
   if (accumulator_.perform_tp_accumulation()) {
@@ -470,6 +549,18 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::collect_measurements(
     concurrency_.delayedSum(accumulator_.get_Gflop());
     accumulated_sign_ = accumulator_.get_accumulated_phase();
     collect_delayed(accumulated_sign_);
+#ifdef DISORDERED_G0
+    // get_sign_times_M_r_w() returns the two-site M(R1,R2,w) here; static_assert pins the domain.
+    static_assert(
+        std::is_same_v<decltype(M_r_r_w_), std::decay_t<decltype(accumulator_.get_sign_times_M_r_w())>>);
+    M_r_r_w_ = accumulator_.get_sign_times_M_r_w();
+    collect_delayed(M_r_r_w_);
+
+    if (accumulator_.compute_std_deviation()) {
+      M_r_r_w_squared_ = accumulator_.get_sign_times_M_r_w_sqr();
+      concurrency_.delayedSum(M_r_r_w_squared_);
+    }
+#else
     static_assert(
         std::is_same_v<decltype(M_r_w_), std::decay_t<decltype(accumulator_.get_sign_times_M_r_w())>>);
     M_r_w_ = accumulator_.get_sign_times_M_r_w();
@@ -479,6 +570,7 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::collect_measurements(
       M_r_w_squared_ = accumulator_.get_sign_times_M_r_w_sqr();
       concurrency_.delayedSum(M_r_w_squared_);
     }
+#endif
 
     if (accumulator_.perform_equal_time_accumulation()) {
       Profiler profiler("Additional time measurements.", "QMC-collectives", __LINE__);
@@ -518,9 +610,16 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::collect_measurements(
     concurrency_.resolveSums();
   }
 
+#ifdef DISORDERED_G0
+  M_r_r_w_ /= static_cast<typename decltype(M_r_r_w_)::this_scalar_type>(accumulated_sign_);
+  M_r_r_w_squared_ /=
+      static_cast<typename decltype(M_r_r_w_squared_)::this_scalar_type>(accumulated_sign_);
+  // The disorder Dyson consumes M_r_r_w_ directly; the single-R M_r_w_ is unused here.
+#else
   M_r_w_ /= static_cast<typename decltype(M_r_w_)::this_scalar_type>(accumulated_sign_);
   M_r_w_squared_ /=
       static_cast<typename decltype(M_r_w_squared_)::this_scalar_type>(accumulated_sign_);
+#endif
   if (accumulator_.perform_tp_accumulation()) {
     for (auto& G4 : data_.get_G4())
       G4 /= static_cast<typename std::remove_reference<decltype(G4)>::type::this_scalar_type>(
@@ -570,8 +669,13 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::symmetrize_measuremen
   if (concurrency_.id() == concurrency_.first())
     std::cout << "\n\t\t symmetrize measurements has started \t" << dca::util::print_time() << "\n";
 
+#ifdef DISORDERED_G0
+  // No Symmetrize overload for the two-site <nu,R,nu,R,w> domain, and crystal symmetry only holds
+  // after disorder averaging, so per-config symmetrization is skipped.
+#else
   Symmetrize<Parameters>::execute(M_r_w_, data_.H_symmetry);
   Symmetrize<Parameters>::execute(M_r_w_squared_, data_.H_symmetry);
+#endif
 }
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
@@ -603,6 +707,7 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::computeG_k_w(
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
 void CtauxClusterSolver<device_t, Parameters, Data, DIST>::compute_G_k_w_from_M_r_w() {
+#ifndef DISORDERED_G0
   func::function<std::complex<double>, NuNuKClusterWDmn> M_k_w;
   math::transform::FunctionTransform<RDmn, KDmn>::execute(M_r_w_, M_k_w);
 
@@ -634,6 +739,78 @@ void CtauxClusterSolver<device_t, Parameters, Data, DIST>::compute_G_k_w_from_M_
   }
 
   Symmetrize<Parameters>::execute(data_.G_k_w, data_.H_symmetry);
+#endif  // DISORDERED_G0
+}
+
+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) {
+#ifndef DISORDERED_G0
+  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_.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);
+    }
+  }
+  G_k_w *= weight;
+  data_.accumulated_G_k_w += G_k_w;
+  // No symmetrization here: crystal symmetry only holds after disorder averaging.
+#endif  // DISORDERED_G0
+}
+
+template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
+void CtauxClusterSolver<device_t, Parameters, Data, DIST>::accumulateGrrwFromMrrw(
+    func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>>& G_r_r_w,
+    double weight) {
+#ifdef DISORDERED_G0
+  // Per-frequency Dyson G = G0_dis - G0_dis*M*G0_dis/beta on the dense (nu*N_R) two-site basis,
+  // accumulated with the config weight into G_r_r_w.
+  disorder::accumulateDisorderDyson<nu, RDmn, w, Real>(
+      data_.disordered_G0_r_r_w_cl_exl, M_r_r_w_, parameters_.get_beta(), weight, G_r_r_w);
+#endif  // DISORDERED_G0
+}
+
+template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
+void CtauxClusterSolver<device_t, Parameters, Data, DIST>::accumulateLocalGrrwFromMrrw(
+    func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>>& G_r_r_w,
+    double weight) {
+#ifdef DISORDERED_G0
+  // Same Dyson as accumulateGrrwFromMrrw, but on THIS rank's own (un-collected) M instead of the
+  // MPI-combined M_r_r_w_. Each rank keeps its own per-config sum so computeErrorBars can take the
+  // cross-rank stddev. Normalized by the per-rank phase, matching the clean error-bar branch.
+  func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>> M_r_r_w_local(
+      accumulator_.get_sign_times_M_r_w(), "M_r_r_w_local");
+  M_r_r_w_local /= static_cast<typename decltype(M_r_r_w_local)::this_scalar_type>(
+      accumulator_.get_accumulated_phase());
+  disorder::accumulateDisorderDyson<nu, RDmn, w, Real>(
+      data_.disordered_G0_r_r_w_cl_exl, M_r_r_w_local, parameters_.get_beta(), weight, G_r_r_w);
+#endif  // DISORDERED_G0
 }
 
 template <dca::linalg::DeviceType device_t, class Parameters, class Data, DistType DIST>
@@ -899,15 +1076,36 @@ auto CtauxClusterSolver<device_t, Parameters, Data, DIST>::local_G_k_w() const {
     throw std::logic_error("The local data was already averaged.");
 
   func::function<std::complex<double>, func::dmn_variadic<nu, nu, KDmn, w>> G_k_w_new("G_k_w_new");
+#ifdef DISORDERED_G0
+  // Per-node local G: the production r-space pipeline (Dyson -> translational average -> FT r->k)
+  // on THIS config alone, with no cross-rank averaging. All intermediates are local scratch.
+
+  // Two-site M(R1,R2,w), normalized by the local (per-node) sign.
+  func::function<std::complex<Real>, NuRNuRClusterWDmn> M_r_r_w_new(
+      accumulator_.get_sign_times_M_r_w(), "M_r_r_w_new");
+  M_r_r_w_new /= static_cast<typename decltype(M_r_r_w_new)::this_scalar_type>(
+      accumulator_.get_accumulated_phase());
+
+  // r-space disorder Dyson, single config (weight 1); G_r_r_w_new must start zeroed to accumulate.
+  func::function<std::complex<Real>, NuRNuRClusterWDmn> G_r_r_w_new("G_r_r_w_new");
+  disorder::accumulateDisorderDyson<nu, RDmn, w, Real>(
+      data_.disordered_G0_r_r_w_cl_exl, M_r_r_w_new, parameters_.get_beta(), 1.0, G_r_r_w_new);
+
+  // Translational average then FT r->k.
+  func::function<std::complex<Real>, NuNuRClusterWDmn> G_r_w_new("G_r_w_new");
+  disorder::translationalAverage<nu, RDmn, w>(G_r_r_w_new, G_r_w_new);
+  math::transform::FunctionTransform<RDmn, KDmn>::execute(G_r_w_new, G_k_w_new);
+#else
   func::function<std::complex<double>, func::dmn_variadic<nu, nu, KDmn, w>> M_k_w_new("M_k_w_new");
   func::function<std::complex<double>, func::dmn_variadic<nu, nu, RDmn, w>> M_r_w_new(
       accumulator_.get_sign_times_M_r_w(), "M_r_w_new");
 
-  M_r_w_new /= accumulator_.get_accumulated_sign();
+  M_r_w_new /= static_cast<typename decltype(M_r_w_new)::this_scalar_type>(
+      accumulator_.get_accumulated_phase());
 
   math::transform::FunctionTransform<RDmn, KDmn>::execute(M_r_w_new, M_k_w_new);
-
   compute_G_k_w_new(M_k_w_new, G_k_w_new);
+#endif  // DISORDERED_G0
 
   return G_k_w_new;
 }
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..c56029564 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
@@ -52,10 +52,17 @@ class G0InterpolationBase {
   typedef typename Parameters::profiler_type profiler_t;
 
   typedef func::dmn_0<domains::time_domain_left_oriented> shifted_t;
-  typedef func::dmn_variadic<nu, nu, r_dmn_t, shifted_t> nu_nu_r_dmn_t_shifted_t;
-
   typedef func::dmn_0<func::dmn<4, int>> akima_dmn_t;
+
+#ifdef DISORDERED_G0
+  // Disorder breaks translational invariance: key G0 by the site pair (R0,R1), not a displacement.
+  typedef func::dmn_variadic<nu, r_dmn_t, nu, r_dmn_t, shifted_t> nu_nu_r_dmn_t_shifted_t;
+  typedef func::dmn_variadic<akima_dmn_t, nu, r_dmn_t, nu, r_dmn_t, shifted_t>
+      akima_nu_nu_r_dmn_t_shifted_t;
+#else
+  typedef func::dmn_variadic<nu, nu, r_dmn_t, shifted_t> nu_nu_r_dmn_t_shifted_t;
   typedef func::dmn_variadic<akima_dmn_t, nu, nu, r_dmn_t, shifted_t> akima_nu_nu_r_dmn_t_shifted_t;
+#endif
 
 public:
   G0InterpolationBase(int id, const Parameters& parameters);
@@ -122,6 +129,29 @@ void G0InterpolationBase<Parameters>::initialize(Data& data) {
 template <typename Parameters>
 template <class Data>
 void G0InterpolationBase<Parameters>::initialize_linear_coefficients(Data& data) {
+#ifdef DISORDERED_G0
+  // Two-site source G0_dis(nu0,R0,nu1,R1,t); same shifted_t packing, extra R0 loop.
+  const auto& g0 = data.disordered_G0_r_r_t_cl_exl;
+  for (int t_ind = 0; t_ind < t::dmn_size() / 2 - 1; t_ind++)
+    for (int r1 = 0; r1 < r_dmn_t::dmn_size(); r1++)
+      for (int r0 = 0; r0 < r_dmn_t::dmn_size(); r0++)
+        for (int nu1 = 0; nu1 < b::dmn_size() * s::dmn_size(); nu1++)
+          for (int nu0 = 0; nu0 < b::dmn_size() * s::dmn_size(); nu0++) {
+            G0_r_t_shifted(nu0, r0, nu1, r1, t_ind) = g0(nu0, r0, nu1, r1, t_ind);
+            grad_G0_r_t_shifted(nu0, r0, nu1, r1, t_ind) =
+                (g0(nu0, r0, nu1, r1, t_ind + 1) - g0(nu0, r0, nu1, r1, t_ind));
+          }
+
+  for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size() - 1; t_ind++)
+    for (int r1 = 0; r1 < r_dmn_t::dmn_size(); r1++)
+      for (int r0 = 0; r0 < r_dmn_t::dmn_size(); r0++)
+        for (int nu1 = 0; nu1 < b::dmn_size() * s::dmn_size(); nu1++)
+          for (int nu0 = 0; nu0 < b::dmn_size() * s::dmn_size(); nu0++) {
+            G0_r_t_shifted(nu0, r0, nu1, r1, t_ind - 1) = g0(nu0, r0, nu1, r1, t_ind);
+            grad_G0_r_t_shifted(nu0, r0, nu1, r1, t_ind - 1) =
+                (g0(nu0, r0, nu1, r1, t_ind + 1) - g0(nu0, r0, nu1, r1, t_ind));
+          }
+#else
   for (int t_ind = 0; t_ind < t::dmn_size() / 2 - 1; t_ind++) {
     for (int r_ind = 0; r_ind < r_dmn_t::dmn_size(); r_ind++) {
       for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) {
@@ -149,6 +179,7 @@ void G0InterpolationBase<Parameters>::initialize_linear_coefficients(Data& data)
       }
     }
   }
+#endif
 }
 
 template <typename Parameters>
@@ -164,6 +195,42 @@ void G0InterpolationBase<Parameters>::initialize_akima_coefficients(Data& data)
   for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++)
     x[t_ind] = t_ind;
 
+#ifdef DISORDERED_G0
+  // Two-site source G0_dis(nu0,R0,nu1,R1,t): one spline per (R0,R1,nu0,nu1), same shifted_t split.
+  const auto& g0 = data.disordered_G0_r_r_t_cl_exl;
+  {  // negative-tau half
+    for (int r1 = 0; r1 < r_dmn_t::dmn_size(); r1++)
+      for (int r0 = 0; r0 < r_dmn_t::dmn_size(); r0++)
+        for (int nu1 = 0; nu1 < b::dmn_size() * s::dmn_size(); nu1++)
+          for (int nu0 = 0; nu0 < b::dmn_size() * s::dmn_size(); nu0++) {
+            for (int t_ind = 0; t_ind < t::dmn_size() / 2; t_ind++)
+              y[t_ind] = g0(nu0, r0, nu1, r1, t_ind);
+
+            ai_obj.initialize(x, y);
+
+            for (int t_ind = 0; t_ind < t::dmn_size() / 2 - 1; t_ind++)
+              for (int l = 0; l < 4; l++)
+                akima_coefficients(l, nu0, r0, nu1, r1, t_ind) = ai_obj.get_alpha(l, t_ind);
+          }
+  }
+
+  {  // positive-tau half
+    for (int r1 = 0; r1 < r_dmn_t::dmn_size(); r1++)
+      for (int r0 = 0; r0 < r_dmn_t::dmn_size(); r0++)
+        for (int nu1 = 0; nu1 < b::dmn_size() * s::dmn_size(); nu1++)
+          for (int nu0 = 0; nu0 < b::dmn_size() * s::dmn_size(); nu0++) {
+            for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size(); t_ind++)
+              y[t_ind - t::dmn_size() / 2] = g0(nu0, r0, nu1, r1, t_ind);
+
+            ai_obj.initialize(x, y);
+
+            for (int t_ind = t::dmn_size() / 2; t_ind < t::dmn_size() - 1; t_ind++)
+              for (int l = 0; l < 4; l++)
+                akima_coefficients(l, nu0, r0, nu1, r1, t_ind - 1) =
+                    ai_obj.get_alpha(l, t_ind - t::dmn_size() / 2);
+          }
+  }
+#else
   {
     for (int r_ind = 0; r_ind < r_dmn_t::dmn_size(); r_ind++) {
       for (int nu1_ind = 0; nu1_ind < b::dmn_size() * s::dmn_size(); nu1_ind++) {
@@ -199,6 +266,7 @@ void G0InterpolationBase<Parameters>::initialize_akima_coefficients(Data& data)
       }
     }
   }
+#endif
 
   /*
     {
diff --git a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_cpu.hpp b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_cpu.hpp
index 1c1c5f8e6..05db69d85 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_cpu.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_cpu.hpp
@@ -66,11 +66,13 @@ class G0Interpolation<dca::linalg::CPU, Parameters> : public G0InterpolationBase
   /// For testing
   auto& getAkimaCoefficients() { return Base::akima_coefficients; }
 private:
-  auto interpolate(int nu_0, int nu_1, int delta_r, Real delta_time) -> Scalar const;
+  // These take both absolute sites (r_0, r_1) rather than a displacement; the clean build rebuilds
+  // the displacement r_1 - r_0 internally.
+  auto interpolate(int nu_0, int r_0, int nu_1, int r_1, Real delta_time) -> Scalar const;
 
-  auto interpolate_on_diagonal(int nu_i) -> Scalar const;
+  auto interpolate_on_diagonal(int nu_i, int r_i) -> Scalar const;
 
-  auto interpolate_akima(int nu_0, int nu_1, int delta_r, Real tau) -> Scalar const;
+  auto interpolate_akima(int nu_0, int r_0, int nu_1, int r_1, Real tau) -> Scalar const;
 
 private:
   using Base::parameters;
@@ -130,22 +132,22 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::build_G0_matrix(
       for (int i = 0; i < j; i++) {
         vertex_singleton_type& v_i = configuration_e_spin[i];
 
-        G0_e_spin(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                            r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+        G0_e_spin(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                            v_j.get_spin_orbital(), v_j.get_r_site(),
                                             (v_i.get_tau() - v_j.get_tau()));
       }
     }
 
     {  // i == j
-      G0_e_spin(j, j) = interpolate_on_diagonal(v_j.get_spin_orbital());
+      G0_e_spin(j, j) = interpolate_on_diagonal(v_j.get_spin_orbital(), v_j.get_r_site());
     }
 
     {  // j > i
       for (int i = j + 1; i < configuration_size; i++) {
         vertex_singleton_type& v_i = configuration_e_spin[i];
 
-        G0_e_spin(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                            r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+        G0_e_spin(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                            v_j.get_spin_orbital(), v_j.get_r_site(),
                                             (v_i.get_tau() - v_j.get_tau()));
       }
     }
@@ -175,22 +177,22 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::build_G0_matrix(
       for (int i = 0; i < j; i++) {
         vertex_singleton_type& v_i = configuration[i];
 
-        G0(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                     r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+        G0(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                     v_j.get_spin_orbital(), v_j.get_r_site(),
                                      (v_i.get_tau() - v_j.get_tau()));
       }
     }
 
     {  // i == j
-      G0(j, j) = interpolate_on_diagonal(v_j.get_spin_orbital());
+      G0(j, j) = interpolate_on_diagonal(v_j.get_spin_orbital(), v_j.get_r_site());
     }
 
     {  // j > i
       for (int i = j + 1; i < configuration_size; i++) {
         vertex_singleton_type& v_i = configuration[i];
 
-        G0(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                     r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+        G0(i, j) = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                     v_j.get_spin_orbital(), v_j.get_r_site(),
                                      (v_i.get_tau() - v_j.get_tau()));
       }
     }
@@ -227,8 +229,8 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::update_G0_matrix(
 
       vertex_singleton_type& v_i = configuration_e_spin[i];
 
-      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                    r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                    v_j.get_spin_orbital(), v_j.get_r_site(),
                                     (v_i.get_tau() - v_j.get_tau()));
     }
   }
@@ -244,13 +246,13 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::update_G0_matrix(
 
       vertex_singleton_type& v_i = configuration_e_spin[i];
 
-      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                    r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                    v_j.get_spin_orbital(), v_j.get_r_site(),
                                     (v_i.get_tau() - v_j.get_tau()));
     }
 
     {  // i == j
-      G0_ptr[j] = interpolate_on_diagonal(v_j.get_spin_orbital());
+      G0_ptr[j] = interpolate_on_diagonal(v_j.get_spin_orbital(), v_j.get_r_site());
     }
 
     // i>j
@@ -259,8 +261,8 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::update_G0_matrix(
 
       vertex_singleton_type& v_i = configuration_e_spin[i];
 
-      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_j.get_spin_orbital(),
-                                    r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+      G0_ptr[i] = interpolate_akima(v_i.get_spin_orbital(), v_i.get_r_site(),
+                                    v_j.get_spin_orbital(), v_j.get_r_site(),
                                     (v_i.get_tau() - v_j.get_tau()));
     }
   }
@@ -277,8 +279,8 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::update_G0_matrix(
 
     vertex_singleton_type& v_i = configuration_e_spin[i];
 
-    G0_ptr[i] = interpolate(v_i.get_spin_orbital(),v_j.get_spin_orbital(),
-    r1_minus_r0(v_j.get_r_site(), v_i.get_r_site()),
+    G0_ptr[i] = interpolate(v_i.get_spin_orbital(), v_i.get_r_site(),
+    v_j.get_spin_orbital(), v_j.get_r_site(),
     (v_i.get_tau()-v_j.get_tau()));
     }
     }
@@ -286,7 +288,7 @@ void G0Interpolation<dca::linalg::CPU, Parameters>::update_G0_matrix(
 }
 
 template <typename Parameters>
-auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate(int nu_0, int nu_1, int delta_r,
+auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate(int nu_0, int r_0, int nu_1, int r_1,
                                                                 Real tau) -> Scalar const {
   // make sure that new_tau is positive !!
   new_tau = tau + beta;
@@ -300,7 +302,11 @@ auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate(int nu_0, int nu
   delta_tau = scaled_tau - t_ind;
   assert(delta_tau > -1.e-16 && delta_tau <= 1 + 1.e-16);
 
-  linind = nu_nu_r_dmn_t_t_shifted_dmn(nu_0, nu_1, delta_r, t_ind);
+#ifdef DISORDERED_G0
+  linind = nu_nu_r_dmn_t_t_shifted_dmn(nu_0, r_0, nu_1, r_1, t_ind);
+#else
+  linind = nu_nu_r_dmn_t_t_shifted_dmn(nu_0, nu_1, r1_minus_r0(r_1, r_0), t_ind);
+#endif
 
   f_0 = G0_r_t_shifted(linind);
   grad = grad_G0_r_t_shifted(linind);
@@ -309,17 +315,23 @@ auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate(int nu_0, int nu
 }
 
 template <typename Parameters>
-auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate_on_diagonal(int nu_i)
+auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate_on_diagonal(int nu_i,
+                                                                            [[maybe_unused]] int r_i)
     -> Scalar const {
   const static int t_0_index = shifted_t::dmn_size() / 2;
+#ifdef DISORDERED_G0
+  // On-site G0_dis(R,R) varies site to site under disorder, so read the vertex's own site r_i
+  // (still the R-diagonal) rather than the clean origin.
+  return -G0_r_t_shifted(nu_i, r_i, nu_i, r_i, t_0_index);
+#else
   const static int r_0_index = r_cluster_type::origin_index();
-
   return -G0_r_t_shifted(nu_i, nu_i, r_0_index, t_0_index);
+#endif
 }
 
 template <typename Parameters>
-auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate_akima(int nu_0, int nu_1,
-                                                                      int delta_r, Real tau)
+auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate_akima(int nu_0, int r_0, int nu_1,
+                                                                      int r_1, Real tau)
     -> Scalar const {
   // make sure that new_tau is positive !!
   new_tau = tau + beta;
@@ -333,7 +345,11 @@ auto G0Interpolation<dca::linalg::CPU, Parameters>::interpolate_akima(int nu_0,
   delta_tau = scaled_tau - t_ind;
   assert(delta_tau > -1.e-16 && delta_tau <= 1 + 1.e-16);
 
-  linind = 4 * nu_nu_r_dmn_t_t_shifted_dmn(nu_0, nu_1, delta_r, t_ind);
+#ifdef DISORDERED_G0
+  linind = 4 * nu_nu_r_dmn_t_t_shifted_dmn(nu_0, r_0, nu_1, r_1, t_ind);
+#else
+  linind = 4 * nu_nu_r_dmn_t_t_shifted_dmn(nu_0, nu_1, r1_minus_r0(r_1, r_0), t_ind);
+#endif
 
   auto* a_ptr = &akima_coefficients(linind);
 
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..405ab627d 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_base.hpp
@@ -41,7 +41,7 @@
 #include "dca/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu.hpp"
 #endif
 
-//#define DEBUG_SUBMATRIX
+// #define DEBUG_SUBMATRIX
 
 namespace dca {
 namespace phys {
diff --git a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu_submatrix.hpp b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu_submatrix.hpp
index f183c1620..75f3db35c 100644
--- a/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu_submatrix.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/ctint/walker/ctint_walker_cpu_submatrix.hpp
@@ -76,15 +76,15 @@ class CtintWalkerSubmatrixCpu : public CtintWalkerSubmatrixBase<Parameters, DIST
 
   DMatrixBuilder<linalg::CPU, Scalar>& d_matrix_builder_;
 
-  BaseClass::MatrixPair getM();
+  typename BaseClass::MatrixPair getM();
 
   /** The following methods are really only here to get decent unit testing
       they shouldn't really be called outside of the base implementations
   */
   void computeMInit() override;
   void computeGInit() override;
-  BaseClass::MatrixPair getRawM();
-  BaseClass::MatrixPair getRawG();
+  typename BaseClass::MatrixPair getRawM();
+  typename BaseClass::MatrixPair getRawG();
 
 private:
   /** This does a bunch of things, this is the majority of a step
@@ -161,7 +161,7 @@ class CtintWalkerSubmatrixCpu : public CtintWalkerSubmatrixBase<Parameters, DIST
   using SubmatrixBase::n_init_;
   using SubmatrixBase::D_;
   using SubmatrixBase::G_;
-  using SubmatrixBase::M_;
+  
   using SubmatrixBase::s_;
   using SubmatrixBase::Gamma_indices_;
   using SubmatrixBase::Gamma_inv_;
@@ -616,7 +616,7 @@ void CtintWalkerSubmatrixCpu<Parameters, DIST>::computeMixedInsertionAndRemoval(
 }
 
 template <class Parameters, DistType DIST>
-CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<
+typename CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<
     Parameters, DIST>::getRawM() {
   typename BaseClass::MatrixPair M;
   M = M_;
@@ -626,7 +626,7 @@ CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubm
 }
 
 template <class Parameters, DistType DIST>
-CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<
+typename CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<
     Parameters, DIST>::getRawG() {
   typename BaseClass::MatrixPair G;
   G = G_;
@@ -636,7 +636,7 @@ CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubm
 }
 
 template <class Parameters, DistType DIST>
-CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<Parameters,
+typename CtintWalkerSubmatrixCpu<Parameters, DIST>::BaseClass::MatrixPair CtintWalkerSubmatrixCpu<Parameters,
                                                                                          DIST>::getM() {
   typename BaseClass::MatrixPair M;
   computeM(M);
diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp
index 66b31d2f7..0681f4562 100644
--- a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp
@@ -54,11 +54,15 @@ class SpAccumulator<Parameters, linalg::CPU> {
   using RDmn = typename Parameters::RClusterDmn;
 
   using NuDmn = func::dmn_variadic<BDmn, SDmn>;  // orbital-spin index
-  using PDmn = func::dmn_variadic<BDmn, BDmn, RDmn>;
 
+#ifdef DISORDERED_G0
+  using PDmn = func::dmn_variadic<BDmn, RDmn, BDmn, RDmn>;
   using MFunction =
-      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RDmn, WDmn>>;
-
+      func::function<std::complex<Real>, func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>>;
+#else
+  using PDmn = func::dmn_variadic<BDmn, BDmn, RDmn>;
+  using MFunction = func::function<std::complex<Real>, func::dmn_variadic<NuDmn, NuDmn, RDmn, WDmn>>;
+#endif
   constexpr static int oversampling = 8;
   using NfftType = math::nfft::Dnfft1D<Scalar, WDmn, PDmn, oversampling, math::nfft::CUBIC>;
   using MFunctionTime = NfftType;
@@ -187,8 +191,11 @@ void SpAccumulator<Parameters, linalg::CPU>::accumulate(
         const Real t_i = config[i].get_tau();
         const int delta_r = RDmn::parameter_type::subtract(r_j, r_i);
         const Real scaled_tau = (t_i - t_j) * one_div_two_beta;  // + (i == j) * epsilon;
-
+#ifdef DISORDERED_G0
+        const int index = bbr_dmn(b_i, r_i, b_j, r_j);
+#else
         const int index = bbr_dmn(b_i, b_j, delta_r);
+#endif
         const Scalar f_val = Ms[s](i, j);
 
         (*cached_nfft_obj_)[s].accumulate(index, scaled_tau, factor * f_val);
@@ -204,18 +211,32 @@ void SpAccumulator<Parameters, linalg::CPU>::accumulate(
 
 template <class Parameters>
 void SpAccumulator<Parameters, linalg::CPU>::finalizeFunction(MFunctionTimePair& ft_objs,
-                                                                    MFunction& function) {
+                                                              MFunction& function) {
   func::function<std::complex<Real>, func::dmn_variadic<WDmn, PDmn>> tmp("tmp");
-  const Real normalization = 1. / RDmn::dmn_size();
 
   for (int s = 0; s < 2; ++s) {
     ft_objs[s].finalize(tmp);
+#ifdef DISORDERED_G0
+    // M_single(dr) = (1/N_R) sum_R0 M_two(R0, R0+dr): the origin-average 1/N_R is deferred to averageGkw, so keep per-pair M here (norm=1).
+    const Real normalization = 1.;
+    // tmp leaves [W,B,R,B,R] (PDmn=<B,R,B,R>); function leaves [B,S,R,B,S,R,W].
+    for (int w_ind = 0; w_ind < WDmn::dmn_size(); w_ind++)
+      for (int r2_ind = 0; r2_ind < RDmn::dmn_size(); r2_ind++)
+        for (int b2_ind = 0; b2_ind < BDmn::dmn_size(); b2_ind++)
+          for (int r1_ind = 0; r1_ind < RDmn::dmn_size(); r1_ind++)
+            for (int b1_ind = 0; b1_ind < BDmn::dmn_size(); b1_ind++)
+              function(b1_ind, s, r1_ind, b2_ind, s, r2_ind, w_ind) +=
+                  tmp(w_ind, b1_ind, r1_ind, b2_ind, r2_ind) * normalization;
+#else
+    // 1/N_R averages the difference-binned M over the absolute origin (the translational average).
+    const Real normalization = 1. / RDmn::dmn_size();
     for (int w_ind = 0; w_ind < WDmn::dmn_size(); w_ind++)
       for (int r_ind = 0; r_ind < RDmn::dmn_size(); r_ind++)
         for (int b2_ind = 0; b2_ind < BDmn::dmn_size(); b2_ind++)
           for (int b1_ind = 0; b1_ind < BDmn::dmn_size(); b1_ind++)
             function(b1_ind, s, b2_ind, s, r_ind, w_ind) +=
                 tmp(w_ind, b1_ind, b2_ind, r_ind) * normalization;
+#endif
   }
 }
 
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 75fee3a43..a28a3e42d 100644
--- a/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp
+++ b/include/dca/phys/dca_step/cluster_solver/stdthread_qmci/stdthread_qmci_cluster_solver.hpp
@@ -59,6 +59,7 @@ class StdThreadQmciClusterSolver : public QmciSolver {
   using typename BaseClass::Accumulator;
   using Walker = stdthreadqmci::StdThreadQmciWalker<typename BaseClass::Walker, Data>;
   using SpGreensFunction = typename BaseClass::SpGreensFunction;
+  using SpDisorderedGreensFunction = typename BaseClass::SpDisorderedGreensFunction;
   using StdThreadAccumulatorType =
       stdthreadqmci::StdThreadQmciAccumulator<Accumulator, typename BaseClass::SpGreensFunction>;
   using MFunction = typename StdThreadAccumulatorType::MFunction;
@@ -116,6 +117,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 +129,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 +196,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 +275,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 +289,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 +360,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 +380,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 +512,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 +612,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) {
@@ -678,13 +701,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 +722,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/cluster/disorder_configuration_domain.hpp b/include/dca/phys/domains/cluster/disorder_configuration_domain.hpp
new file mode 100644
index 000000000..a93243ac4
--- /dev/null
+++ b/include/dca/phys/domains/cluster/disorder_configuration_domain.hpp
@@ -0,0 +1,72 @@
+// 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: Anirudha Mirmira
+//
+// Runtime-sized index domain enumerating the disorder configurations of a DCA step. It is
+// used to assemble the combined "disorder-randoms" function written when the disorder
+// "dump-randoms" parameter is set. Its size is established by initialize(num_configurations)
+// and is not part of the global domain initialization.
+
+#ifndef DCA_PHYS_DOMAINS_CLUSTER_DISORDER_CONFIGURATION_DOMAIN_HPP
+#define DCA_PHYS_DOMAINS_CLUSTER_DISORDER_CONFIGURATION_DOMAIN_HPP
+
+#include <numeric>
+#include <string>
+#include <vector>
+
+namespace dca {
+namespace phys {
+namespace domains {
+// dca::phys::domains::
+
+class DisorderConfigurationDomain {
+public:
+  using scalar_type = int;
+  using element_type = int;
+
+  // Set the number of disorder configurations. Call before the domain is used.
+  static void initialize(int num_configurations);
+
+  // Returns the number of disorder configurations.
+  static int get_size() {
+    return elements_.size();
+  }
+
+  // The configuration indices 0 .. get_size()-1.
+  static const std::vector<int>& get_elements() {
+    return elements_;
+  }
+
+  static const std::string& get_name() {
+    const static std::string name = "Disorder configuration domain.";
+    return name;
+  }
+
+  template <class Writer>
+  static void write(Writer& writer);
+
+private:
+  inline static std::vector<int> elements_{};
+};
+
+inline void DisorderConfigurationDomain::initialize(int num_configurations) {
+  elements_.resize(num_configurations);
+  std::iota(elements_.begin(), elements_.end(), 0);
+}
+
+template <class Writer>
+void DisorderConfigurationDomain::write(Writer& writer) {
+  writer.open_group(get_name());
+  writer.execute("element_indices_", elements_);
+  writer.close_group();
+}
+
+}  // namespace domains
+}  // namespace phys
+}  // namespace dca
+
+#endif  // DCA_PHYS_DOMAINS_CLUSTER_DISORDER_CONFIGURATION_DOMAIN_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/domains/time_and_frequency/time_domain.hpp b/include/dca/phys/domains/time_and_frequency/time_domain.hpp
index f7bdb02e5..e6b17655a 100644
--- a/include/dca/phys/domains/time_and_frequency/time_domain.hpp
+++ b/include/dca/phys/domains/time_and_frequency/time_domain.hpp
@@ -32,6 +32,7 @@ class time_domain {
   using element_type = double;
   using scalar_type = double;
 
+  static constexpr int DIMENSION = 1;
   // Needed in function transform.
   using dmn_specifications_type = math::transform::interval_dmn_1D_type;
 
@@ -114,8 +115,8 @@ void time_domain::to_JSON(stream_type& ss) {
   ss << "]\n";
 }
 
-}  // domains
-}  // phys
-}  // dca
+}  // namespace domains
+}  // namespace phys
+}  // namespace dca
 
 #endif  // DCA_PHYS_DOMAINS_TIME_AND_FREQUENCY_TIME_DOMAIN_HPP
diff --git a/include/dca/phys/parameters/disorder_parameters.hpp b/include/dca/phys/parameters/disorder_parameters.hpp
new file mode 100644
index 000000000..7c2389de7
--- /dev/null
+++ b/include/dca/phys/parameters/disorder_parameters.hpp
@@ -0,0 +1,178 @@
+// 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>
+#include <string>
+
+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()
+      : disorder_distribution_("binary"),
+        disorder_potential_(0.0),
+        disorder_density_(0.0),
+        disorder_num_configurations_(0),
+        // disorder_max_sites_(1),
+        disorder_unique_configs_(false),
+        disorder_present_(false) {}
+
+  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);
+
+  const std::string& get_disorder_distribution() const {
+    return disorder_distribution_;
+  }
+  double get_disorder_potential() const {
+    return disorder_potential_;
+  }
+  void set_disorder_potential(double potential) {
+    disorder_potential_ = potential;
+  }
+  double get_disorder_density() const {
+    return disorder_density_;
+  }
+  int get_disorder_num_configurations() const {
+    return disorder_num_configurations_;
+  }
+  void set_disorder_num_configurations(int num_configurations) {
+    disorder_num_configurations_ = num_configurations;
+  }
+  // int get_disorder_max_sites() const {
+  //   return disorder_max_sites_;
+  // }
+  bool get_disorder_unique_configs() const {
+    return disorder_unique_configs_;
+  }
+  // Whether a "disorder" section was present in the input. Lets the disorder flow distinguish a
+  // deliberate clean run (no section) from a misconfigured one (section present, zero configs).
+  bool get_disorder_present() const {
+    return disorder_present_;
+  }
+
+private:
+  std::string disorder_distribution_;
+  double disorder_potential_;
+  double disorder_density_;
+  int disorder_num_configurations_;
+  // int disorder_max_sites_;
+  bool disorder_unique_configs_;
+  bool disorder_present_;
+};
+
+template <typename Concurrency>
+int DisorderParameters::getBufferSize(const Concurrency& concurrency) const {
+  int buffer_size = 0;
+
+  buffer_size += concurrency.get_buffer_size(disorder_distribution_);
+  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_);
+  buffer_size += concurrency.get_buffer_size(disorder_unique_configs_);
+  buffer_size += concurrency.get_buffer_size(disorder_present_);
+
+  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_distribution_);
+  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_);
+  concurrency.pack(buffer, buffer_size, position, disorder_unique_configs_);
+  concurrency.pack(buffer, buffer_size, position, disorder_present_);
+}
+
+template <typename Concurrency>
+void DisorderParameters::unpack(const Concurrency& concurrency, char* buffer, int buffer_size,
+                                int& position) {
+  concurrency.unpack(buffer, buffer_size, position, disorder_distribution_);
+  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_);
+  concurrency.unpack(buffer, buffer_size, position, disorder_unique_configs_);
+  concurrency.unpack(buffer, buffer_size, position, disorder_present_);
+}
+
+template <typename ReaderOrWriter>
+void DisorderParameters::readWrite(ReaderOrWriter& reader_or_writer) {
+  try {
+    disorder_present_ = reader_or_writer.open_group("disorder");
+
+    try {
+      reader_or_writer.execute("distribution", disorder_distribution_);
+    }
+    catch (const std::exception& r_e) {
+    }
+    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) {
+    // }
+    try {
+      reader_or_writer.execute("unique-configs", disorder_unique_configs_);
+    }
+    catch (const std::exception& r_e) {
+    }
+    reader_or_writer.close_group();
+  }
+  catch (const std::exception& r_e) {
+  }
+
+  // Check values.
+  if (!(disorder_distribution_ == "binary" || disorder_distribution_ == "box"))
+    throw std::logic_error("Illegal value for disorder distribution (expected 'binary' or 'box').");
+}
+
+}  // 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/mci_parameters.hpp b/include/dca/phys/parameters/mci_parameters.hpp
index 8624733bf..ddcf5e1b6 100644
--- a/include/dca/phys/parameters/mci_parameters.hpp
+++ b/include/dca/phys/parameters/mci_parameters.hpp
@@ -160,6 +160,8 @@ class MciParameters {
   bool fix_meas_per_walker_ = true;
   bool adjust_self_energy_for_double_counting_;
   ErrorComputationType error_computation_type_;
+  /// actually a runtime variable that should be moved from
+  /// parameters.
   bool store_configuration_;
   DistType g4_distribution_;
   int stamping_period_ = 0;
@@ -322,41 +324,43 @@ void MciParameters::readWrite(ReaderOrWriter& reader_or_writer) {
   reader_or_writer.close_group();
 
   if constexpr (ReaderOrWriter::is_reader) {
-    // The input file can contain an integral seed or the seeding option "random".  
-  // Check parameters consistency.
-  if (g4_distribution_ == DistType::BLOCKED) {
+    // The input file can contain an integral seed or the seeding option "random".
+    // Check parameters consistency.
+    if (g4_distribution_ == DistType::BLOCKED) {
 #ifdef DCA_HAVE_MPI
-    // Check for number of accumulators and walkers consistency.
-    if (!shared_walk_and_accumulation_thread_ || walkers_ != accumulators_) {
-      throw std::logic_error(
-          "\n With distributed g4 enabled, 1) walker and accumulator should share "
-          "thread, "
-          "2) #walker == #accumulator\n");
-    }
+      // Check for number of accumulators and walkers consistency.
+      if (!shared_walk_and_accumulation_thread_ || walkers_ != accumulators_) {
+        throw std::logic_error(
+            "\n With distributed g4 enabled, 1) walker and accumulator should share "
+            "thread, "
+            "2) #walker == #accumulator\n");
+      }
 
-    // Check for number of ranks and g4 measurements consistency.
-    // This is potentially wrong if fancy rank ganging or the like is used.
-    int mpi_size;
-    MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
-    int local_meas = measurements_.back() / mpi_size;
-    if (measurements_.back() % mpi_size != 0 && local_meas % accumulators_ != 0) {
-      throw std::logic_error(
-          "\n With distributed g4 enabled, 1) local measurements should be same across "
-          "ranks, "
-          "2) each accumulator should have same measurements\n");
-    }
+      // Check for number of ranks and g4 measurements consistency.
+      // This is potentially wrong if fancy rank ganging or the like is used.
+      int mpi_size;
+      MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
+      int local_meas = measurements_.back() / mpi_size;
+      if (measurements_.back() % mpi_size != 0 && local_meas % accumulators_ != 0) {
+        throw std::logic_error(
+            "\n With distributed g4 enabled, 1) local measurements should be same across "
+            "ranks, "
+            "2) each accumulator should have same measurements\n");
+      }
 #else
-    throw(std::logic_error("MPI distribution requested with no MPI available."));
+      throw(std::logic_error("MPI distribution requested with no MPI available."));
 #endif  // DCA_HAVE_MPI
-    if (stamping_period_ != 0) {
-      if (!(shared_walk_and_accumulation_thread_ && walkers_ == accumulators_))
-	throw std::runtime_error("Individual measurement stamping not available unless shared-walk-and-accumulation-thread = true and walkers == acceptors!");
+      if (stamping_period_ != 0) {
+        if (!(shared_walk_and_accumulation_thread_ && walkers_ == accumulators_))
+          throw std::runtime_error(
+              "Individual measurement stamping not available unless "
+              "shared-walk-and-accumulation-thread = true and walkers == acceptors!");
+      }
+      // Solve conflicts
     }
-    // Solve conflicts
-  }
   }
   if (!time_correlation_window_)
-      compute_G_correlation_ = false;
+    compute_G_correlation_ = false;
 }
 
 void MciParameters::solveDcaIterationConflict(int iterations) {
diff --git a/include/dca/phys/parameters/output_parameters.hpp b/include/dca/phys/parameters/output_parameters.hpp
index c1b870074..ad2fb7c1f 100644
--- a/include/dca/phys/parameters/output_parameters.hpp
+++ b/include/dca/phys/parameters/output_parameters.hpp
@@ -52,7 +52,8 @@ class OutputParameters {
         dump_lattice_self_energy_(false),
         dump_cluster_Greens_functions_(false),
         dump_Gamma_lattice_(false),
-        dump_chi_0_lattice_(false) {
+        dump_chi_0_lattice_(false),
+        dump_disorder_configs_(false) {
   }
 
   template <typename Concurrency>
@@ -129,6 +130,9 @@ class OutputParameters {
   bool dump_chi_0_lattice() const {
     return dump_chi_0_lattice_;
   }
+  bool dump_disorder_configs() const {
+    return dump_disorder_configs_;
+  }
   bool dump_every_iteration() const {
     return dump_every_iteration_;
   }
@@ -151,6 +155,7 @@ class OutputParameters {
   bool dump_cluster_Greens_functions_;
   bool dump_Gamma_lattice_;
   bool dump_chi_0_lattice_;
+  bool dump_disorder_configs_;
   bool dump_every_iteration_ = false;
 };
 
@@ -175,6 +180,7 @@ int OutputParameters::getBufferSize(const Concurrency& concurrency) const {
   buffer_size += concurrency.get_buffer_size(dump_cluster_Greens_functions_);
   buffer_size += concurrency.get_buffer_size(dump_Gamma_lattice_);
   buffer_size += concurrency.get_buffer_size(dump_chi_0_lattice_);
+  buffer_size += concurrency.get_buffer_size(dump_disorder_configs_);
   buffer_size += concurrency.get_buffer_size(dump_every_iteration_);
 
   return buffer_size;
@@ -200,6 +206,7 @@ void OutputParameters::pack(const Concurrency& concurrency, char* buffer, int bu
   concurrency.pack(buffer, buffer_size, position, dump_cluster_Greens_functions_);
   concurrency.pack(buffer, buffer_size, position, dump_Gamma_lattice_);
   concurrency.pack(buffer, buffer_size, position, dump_chi_0_lattice_);
+  concurrency.pack(buffer, buffer_size, position, dump_disorder_configs_);
   concurrency.pack(buffer, buffer_size, position, dump_every_iteration_);
 }
 
@@ -223,6 +230,7 @@ void OutputParameters::unpack(const Concurrency& concurrency, char* buffer, int
   concurrency.unpack(buffer, buffer_size, position, dump_cluster_Greens_functions_);
   concurrency.unpack(buffer, buffer_size, position, dump_Gamma_lattice_);
   concurrency.unpack(buffer, buffer_size, position, dump_chi_0_lattice_);
+  concurrency.unpack(buffer, buffer_size, position, dump_disorder_configs_);
   concurrency.unpack(buffer, buffer_size, position, dump_every_iteration_);
 }
 
@@ -256,6 +264,7 @@ void OutputParameters::readWrite(ReaderOrWriter& reader_or_writer) {
     try_to_read_or_write("dump-cluster-Greens-functions", dump_cluster_Greens_functions_);
     try_to_read_or_write("dump-Gamma-lattice", dump_Gamma_lattice_);
     try_to_read_or_write("dump-chi-0-lattice", dump_chi_0_lattice_);
+    try_to_read_or_write("dump-disorder-configs", dump_disorder_configs_);
     try_to_read_or_write("dump-every-iteration", dump_every_iteration_);
 
     reader_or_writer.close_group();
diff --git a/include/dca/phys/parameters/parameters.hpp b/include/dca/phys/parameters/parameters.hpp
index b9bf52383..4dc500758 100644
--- a/include/dca/phys/parameters/parameters.hpp
+++ b/include/dca/phys/parameters/parameters.hpp
@@ -31,6 +31,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"
@@ -69,6 +70,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,
@@ -169,10 +171,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>
@@ -180,6 +184,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(),
@@ -206,8 +211,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,
@@ -270,8 +277,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);
 }
 
@@ -394,6 +400,7 @@ int Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator, s
 
   buffer_size += AnalysisParameters::getBufferSize(concurrency);
   buffer_size += DcaParameters::getBufferSize(concurrency);
+  buffer_size += DisorderParameters::getBufferSize(concurrency);
   buffer_size += DomainsParameters::getBufferSize(concurrency);
   buffer_size += DoubleCountingParameters::getBufferSize(concurrency);
   buffer_size += EdSolverParameters::getBufferSize(concurrency);
@@ -414,6 +421,7 @@ void Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator,
                                  int& position) const {
   AnalysisParameters::pack(concurrency, buffer, buffer_size, position);
   DcaParameters::pack(concurrency, buffer, buffer_size, position);
+  DisorderParameters::pack(concurrency, buffer, buffer_size, position);
   DomainsParameters::pack(concurrency, buffer, buffer_size, position);
   DoubleCountingParameters::pack(concurrency, buffer, buffer_size, position);
   EdSolverParameters::pack(concurrency, buffer, buffer_size, position);
@@ -432,6 +440,7 @@ void Parameters<Concurrency, Threading, Profiler, Model, RandomNumberGenerator,
                                    int& position) {
   AnalysisParameters::unpack(concurrency, buffer, buffer_size, position);
   DcaParameters::unpack(concurrency, buffer, buffer_size, position);
+  DisorderParameters::unpack(concurrency, buffer, buffer_size, position);
   DomainsParameters::unpack(concurrency, buffer, buffer_size, position);
   DoubleCountingParameters::unpack(concurrency, buffer, buffer_size, position);
   EdSolverParameters::unpack(concurrency, buffer, buffer_size, position);
@@ -462,6 +471,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_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/cluster/CMakeLists.txt b/src/phys/domains/cluster/CMakeLists.txt
index efa0df526..557719833 100644
--- a/src/phys/domains/cluster/CMakeLists.txt
+++ b/src/phys/domains/cluster/CMakeLists.txt
@@ -4,5 +4,6 @@ add_library(cluster_domains STATIC
   cluster_definitions.cpp
   momentum_exchange_domain.cpp
   symmetries/symmetry_operations/group_action.cpp)
+#  value_at_debug.cpp)
 
 dca_gpu_runtime_link(cluster_domains)
diff --git a/src/phys/domains/quantum/CMakeLists.txt b/src/phys/domains/quantum/CMakeLists.txt
index 2a1c2bc2d..3be4f7f54 100644
--- a/src/phys/domains/quantum/CMakeLists.txt
+++ b/src/phys/domains/quantum/CMakeLists.txt
@@ -4,6 +4,7 @@ add_library(quantum_domains STATIC
   brillouin_zone_path_domain.cpp
   dca_iteration_domain.cpp
   electron_spin_domain.cpp
+  electron_band_domain.cpp
   numerical_error_domain.cpp
   point_group_symmetry_element.cpp
   symmetry_group_level.cpp)
diff --git a/src/phys/domains/quantum/electron_band_domain.cpp b/src/phys/domains/quantum/electron_band_domain.cpp
new file mode 100644
index 000000000..4b7055d5f
--- /dev/null
+++ b/src/phys/domains/quantum/electron_band_domain.cpp
@@ -0,0 +1,24 @@
+// Copyright (C) 2018 ETH Zurich
+// Copyright (C) 2018 UT-Battelle, LLC
+// All rights reserved.
+//
+// See LICENSE for terms of usage.
+// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
+//
+// Author: Peter Staar (taa@zurich.ibm.com)
+//
+// This file implements electron_spin_domain.hpp.
+
+#include "dca/phys/domains/quantum/electron_band_domain.hpp"
+#include <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/integration/statistical_tests/square_lattice/CMakeLists.txt b/test/integration/statistical_tests/square_lattice/CMakeLists.txt
index 2321e4702..43c7ca18b 100644
--- a/test/integration/statistical_tests/square_lattice/CMakeLists.txt
+++ b/test/integration/statistical_tests/square_lattice/CMakeLists.txt
@@ -37,6 +37,27 @@ dca_add_gtest(ctint_square_lattice_validation_stattest
     LIBS     ${TEST_LIBS}
     )
 
+# Disorder clean-limit stat test (Step 1): drives the r-space disorder pipeline at V(R) = 0 and
+# checks it reproduces the clean ED reference. Only built with the disorder option, never for GPU.
+if(DCA_WITH_DISORDER AND NOT DCA_WITH_CUDA AND NOT DCA_WITH_HIP)
+  dca_add_gtest(disorder_ctaux_square_lattice_stattest
+      STOCHASTIC
+      # Use the machine's core count rather than a fixed number; more ranks give better statistics.
+      MPI MPI_NUMPROC ${CPUS}
+      INCLUDE_DIRS ${TEST_INCLUDES}
+      LIBS     ${TEST_LIBS}
+      )
+
+  # Disorder finite-V robustness test (Step 2): runs the binary + box disorder ensemble at V=0.1 and
+  # checks G_k_w stays finite/sane. A smoke test (no accuracy/baseline); grouped with the stochastic
+  # tests so it is not built unless DCA_WITH_TESTS_STOCHASTIC is on.
+  dca_add_gtest(disorder_finite_v_stattest
+      STOCHASTIC
+      MPI MPI_NUMPROC 1
+      INCLUDE_DIRS ${TEST_INCLUDES}
+      LIBS     ${TEST_LIBS}
+      )
+endif()
 
 configure_file(data.ed.hdf5 ${CMAKE_CURRENT_BINARY_DIR}/data.ed.hdf5 COPYONLY)
 configure_file(verification_covariance_input.hdf5 ${CMAKE_CURRENT_BINARY_DIR}/verification_covariance_input.hdf5 COPYONLY)
diff --git a/test/integration/statistical_tests/square_lattice/disorder_ctaux_square_lattice_stattest.cpp b/test/integration/statistical_tests/square_lattice/disorder_ctaux_square_lattice_stattest.cpp
new file mode 100644
index 000000000..9f4e9c6d5
--- /dev/null
+++ b/test/integration/statistical_tests/square_lattice/disorder_ctaux_square_lattice_stattest.cpp
@@ -0,0 +1,140 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Disorder stat test, Step 1 (clean limit): drives the CT-AUX r-space disorder pipeline with a
+// zero potential (V(R) = 0) and checks that the QMC-sampled local_G_k_w statistically agrees
+// with the clean ED reference. With V = 0 the per-configuration r-space Dyson, the translational
+// average, and the FT r->k must collapse to the clean answer, so any corruption introduced by the
+// disorder machinery (driven here by genuine M sampled at U = 4) shows up as a failed p-value.
+// This is the clean-limit regression of the full disorder code path.
+//
+// Only built when DCA_WITH_DISORDER is ON (see CMakeLists.txt), which is CPU-only, so there is no
+// GPU branch here; the clean path is already covered by ctaux_square_lattice_validation_stattest.
+
+#include <string>
+#include <iostream>
+#include <cmath>
+
+#include "gtest/gtest.h"
+
+#include "dca/math/statistical_testing/statistical_testing.hpp"
+#include "dca/config/profiler.hpp"
+
+// The shared setup header uses a global Scalar and includes ctint_cluster_solver.hpp, which needs
+// dca::config::McOptions. Provide both before pulling in the setup header (same pattern as the
+// rashba stat test). The square lattice Green's function is real, so Scalar is double.
+using Scalar = double;
+
+#include "test/mock_mcconfig.hpp"
+namespace dca {
+namespace config {
+using McOptions = MockMcOptions<Scalar>;
+}  // namespace config
+}  // namespace dca
+
+#include "test/integration/statistical_tests/square_lattice/square_lattice_setup.hpp"
+
+dca::testing::DcaMpiTestEnvironment* dca_test_env;
+
+TEST(DisorderCtauxSquareLatticeCleanLimitTest, GreensFunction) {
+  using namespace dca::testing;
+  const std::string ed_data_name = "data.ed.hdf5";
+
+  const int id = dca_test_env->concurrency.id();
+  const int number_of_samples = dca_test_env->concurrency.number_of_processors();
+
+  if (id == dca_test_env->concurrency.first()) {
+    dca::util::GitVersion::print();
+    dca::util::Modules::print();
+  }
+
+  ParametersType<ClusterSolverId::CT_AUX> parameters(dca::util::GitVersion::string(),
+                                                     dca_test_env->concurrency);
+  parameters.read_input_and_broadcast<dca::io::JSONReader>(dca_test_env->input_file_name);
+  parameters.update_model();
+  parameters.update_domains();
+
+  parameters.set_measurements(parameters.get_measurements().back() * number_of_samples / 50);
+
+  DcaData<ClusterSolverId::CT_AUX> data(parameters);
+  data.initialize();
+
+  // The stat test drives the solver directly (initialize -> integrate -> local_G_k_w) and never
+  // calls DcaLoop::workTheClusters, which is what normally populates the disordered bare
+  // propagator data.disordered_G0_r_r_{t,w}_cl_exl that the walker and the post-integration
+  // r-space Dyson read. Populate it here with a zero potential (V(R) = 0); the default-constructed
+  // DisorderConfiguration is identically zero.
+  typename DcaData<ClusterSolverId::CT_AUX>::DisorderConfiguration zero_config;
+  data.makeDisorderedG0(zero_config);
+
+  // Do one QMC iteration.
+  QuantumClusterSolver<ClusterSolverId::CT_AUX> qmc_solver(parameters, data, nullptr);
+  qmc_solver.initialize(0);
+  qmc_solver.integrate();
+
+  using dca::func::function;
+  function<double, SigmaCutDomain> G_k_w_sample =
+      cutFrequency(qmc_solver.local_G_k_w(), n_frequencies);
+  G_k_w_sample.set_name("G_k_w");
+
+  // Catastrophic-failure guard: the clean-limit disorder pipeline must produce finite numbers.
+  for (int i = 0; i < G_k_w_sample.size(); ++i)
+    ASSERT_TRUE(std::isfinite(G_k_w_sample(i)))
+        << "Non-finite value in disorder clean-limit local_G_k_w at linear index " << i;
+
+  // Read the expected (clean) ED result.
+  function<double, SigmaCutDomain> G_k_w_expected;
+  if (id == dca_test_env->concurrency.first()) {
+    function<std::complex<double>, SigmaDomain> G_k_w_full("cluster_greens_function_G_k_w");
+    dca::io::HDF5Reader reader;
+    reader.open_file(ed_data_name);
+    reader.open_group("functions");
+    reader.execute(G_k_w_full);
+    reader.close_group();
+    reader.close_file();
+    G_k_w_expected = cutFrequency(G_k_w_full, n_frequencies);
+  }
+  dca_test_env->concurrency.broadcast(G_k_w_expected);
+
+  // Compute covariance and average the QMC result across ranks.
+  function<double, CovarianceDomain> G_k_w_covariance("G_k_w_covariance");
+  dca_test_env->concurrency.computeCovarianceAndAvg(G_k_w_sample, G_k_w_covariance);
+
+  // Compute the p-value against the clean ED reference.
+  if (id == dca_test_env->concurrency.first()) {
+    dca::math::StatisticalTesting test(G_k_w_sample, G_k_w_expected, G_k_w_covariance, 1);
+    double p_value = test.computePValue(false, number_of_samples);
+    test.printInfo("disorder_ctaux_square_testinfo.out", true);
+    double p_value_default = 0.05;
+    std::cout << "\n***\nThe p-value is " << p_value << "\n***\n";
+    EXPECT_LT(p_value_default, p_value);
+  }
+}
+
+int main(int argc, char** argv) {
+  int result = 0;
+
+  ::testing::InitGoogleTest(&argc, argv);
+
+  dca::parallel::MPIConcurrency concurrency(argc, argv);
+  dca_test_env = new dca::testing::DcaMpiTestEnvironment(
+      concurrency, dca::testing::test_directory + "square_lattice_input.json");
+  ::testing::AddGlobalTestEnvironment(dca_test_env);
+
+  ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
+
+  if (dca_test_env->concurrency.id() != 0) {
+    delete listeners.Release(listeners.default_result_printer());
+    listeners.Append(new dca::testing::MinimalistPrinter);
+  }
+
+  result = RUN_ALL_TESTS();
+
+  return result;
+}
diff --git a/test/integration/statistical_tests/square_lattice/disorder_finite_v_stattest.cpp b/test/integration/statistical_tests/square_lattice/disorder_finite_v_stattest.cpp
new file mode 100644
index 000000000..a97607150
--- /dev/null
+++ b/test/integration/statistical_tests/square_lattice/disorder_finite_v_stattest.cpp
@@ -0,0 +1,155 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Disorder stat test, Step 2 (finite disorder robustness): runs the full CT-AUX disorder ensemble
+// average at a small finite potential (V = 0.1) for both the binary and box distributions, and
+// checks that the disorder-averaged G_k_w stays finite and sane (no NaN/Inf, nonzero, bounded).
+// The aim is NOT accuracy -- it is to confirm the disorder path does not blow up catastrophically
+// once V != 0. It replicates the production disorder loop (DcaLoop::workTheClusters + averageGkw):
+// generate the ensemble -> per config makeDisorderedG0 -> integrate -> accumulateGkw(weight) ->
+// normalize -> translational average -> FT r->k.
+//
+// Only built when DCA_WITH_DISORDER is ON (see CMakeLists.txt), which is CPU-only.
+
+#include <cmath>
+#include <numeric>
+#include <string>
+#include <vector>
+
+#include "gtest/gtest.h"
+
+#include "dca/config/profiler.hpp"
+#include "dca/io/json/json_reader.hpp"
+
+// The shared setup header uses a global Scalar and includes ctint_cluster_solver.hpp, which needs
+// dca::config::McOptions. Provide both before pulling in the setup header (same pattern as the
+// rashba stat test). The square lattice Green's function is real, so Scalar is double.
+using Scalar = double;
+
+#include "test/mock_mcconfig.hpp"
+namespace dca {
+namespace config {
+using McOptions = MockMcOptions<Scalar>;
+}  // namespace config
+}  // namespace dca
+
+#include "test/integration/statistical_tests/square_lattice/square_lattice_setup.hpp"
+
+// These pull in dca_shared_types.hpp / domain_aliases.hpp, so they must follow the setup header
+// (which defines the time/frequency/cluster domains they alias).
+#include "dca/math/function_transform/function_transform.hpp"
+#include "dca/phys/dca_loop/disorder/disorder_average.hpp"
+#include "dca/phys/dca_loop/disorder/makeDisorderConfigurations.hpp"
+
+dca::testing::DcaMpiTestEnvironment* dca_test_env;
+
+namespace {
+
+using dca::ClusterSolverId;
+using SolverParameters = dca::testing::ParametersType<ClusterSolverId::CT_AUX>;
+using SolverData = dca::testing::DcaData<ClusterSolverId::CT_AUX>;
+
+// Drive the full disorder ensemble average for the given input file and return the disorder-averaged
+// G_k_w. Mirrors DcaLoop::workTheClusters (the disorder loop) followed by DcaLoop::averageGkw.
+SolverData::SpGreensFunction runDisorderEnsemble(const std::string& input_file) {
+  SolverParameters parameters(dca::util::GitVersion::string(), dca_test_env->concurrency);
+  parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+  parameters.update_model();
+  parameters.update_domains();
+
+  SolverData data(parameters);
+  data.initialize();
+
+  // Generate the disorder ensemble (binary or box, per the input's distribution).
+  dca::phys::MakeDisorderConfigurations<SolverParameters> make_configs;
+  std::vector<typename SolverData::DisorderConfiguration> configs;
+  std::vector<double> weights;
+  make_configs(parameters, configs, weights);
+  EXPECT_GT(configs.size(), 0u) << "disorder generator produced no configurations";
+
+  dca::testing::QuantumClusterSolver<ClusterSolverId::CT_AUX> solver(parameters, data, nullptr);
+
+  // --- DcaLoop::workTheClusters disorder loop ---
+  data.disorder_G_r_r_w = 0;
+  for (std::size_t i = 0; i < configs.size(); ++i) {
+    data.makeDisorderedG0(configs[i]);
+    solver.initialize(0);
+    solver.integrate();
+    solver.accumulateGkw(weights[i]);
+  }
+
+  // Catastrophic-failure guard on the two-site accumulator, before the transforms, to localize a NaN.
+  for (int i = 0; i < data.disorder_G_r_r_w.size(); ++i)
+    EXPECT_TRUE(std::isfinite(data.disorder_G_r_r_w(i).real()) &&
+                std::isfinite(data.disorder_G_r_r_w(i).imag()))
+        << "non-finite disorder_G_r_r_w at linear index " << i;
+
+  // --- DcaLoop::averageGkw ---
+  const double total_weight = std::accumulate(weights.begin(), weights.end(), 0.0);
+  data.disorder_G_r_r_w /= total_weight;
+
+  using NuDmn = SolverData::NuDmn;
+  using RClusterDmn = SolverData::RClusterDmn;
+  using KClusterDmn = SolverData::KClusterDmn;
+  using WDmn = SolverData::WDmn;
+
+  dca::phys::disorder::translationalAverage<NuDmn, RClusterDmn, WDmn>(data.disorder_G_r_r_w,
+                                                                      data.G_r_w);
+  dca::math::transform::FunctionTransform<RClusterDmn, KClusterDmn>::execute(data.G_r_w, data.G_k_w);
+
+  return data.G_k_w;
+}
+
+// Hard robustness assertions: every element finite, the result is not identically zero, and its
+// magnitude is bounded (a blow-up would show up as a huge or non-finite value).
+void expectFiniteAndSane(const SolverData::SpGreensFunction& g_k_w) {
+  double max_abs = 0.;
+  for (int i = 0; i < g_k_w.size(); ++i) {
+    ASSERT_TRUE(std::isfinite(g_k_w(i).real()) && std::isfinite(g_k_w(i).imag()))
+        << "non-finite G_k_w at linear index " << i;
+    max_abs = std::max(max_abs, std::abs(g_k_w(i)));
+  }
+  EXPECT_GT(max_abs, 0.) << "disorder-averaged G_k_w is identically zero";
+  // |G(k,iw)| <= beta for a fermionic propagator; 1e3 is a very loose blow-up ceiling.
+  EXPECT_LT(max_abs, 1e3) << "disorder-averaged G_k_w magnitude is implausibly large (blow-up)";
+}
+
+}  // namespace
+
+TEST(DisorderFiniteVTest, BinaryDisorderIsFinite) {
+  expectFiniteAndSane(
+      runDisorderEnsemble(dca::testing::test_directory + "input_disorder_finite_v_binary.json"));
+}
+
+TEST(DisorderFiniteVTest, BoxDisorderIsFinite) {
+  expectFiniteAndSane(
+      runDisorderEnsemble(dca::testing::test_directory + "input_disorder_finite_v_box.json"));
+}
+
+int main(int argc, char** argv) {
+  int result = 0;
+
+  ::testing::InitGoogleTest(&argc, argv);
+
+  dca::parallel::MPIConcurrency concurrency(argc, argv);
+  dca_test_env = new dca::testing::DcaMpiTestEnvironment(
+      concurrency, dca::testing::test_directory + "input_disorder_finite_v_binary.json");
+  ::testing::AddGlobalTestEnvironment(dca_test_env);
+
+  ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
+
+  if (dca_test_env->concurrency.id() != 0) {
+    delete listeners.Release(listeners.default_result_printer());
+    listeners.Append(new dca::testing::MinimalistPrinter);
+  }
+
+  result = RUN_ALL_TESTS();
+
+  return result;
+}
diff --git a/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_binary.json b/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_binary.json
new file mode 100644
index 000000000..35c749faa
--- /dev/null
+++ b/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_binary.json
@@ -0,0 +1,51 @@
+{
+  "physics": {
+    "beta"                      : 1.0,
+    "chemical-potential"        : 0.0
+  },
+
+  "single-band-Hubbard-model": {
+    "t"       : 1.0,
+    "U"       : 4.0
+  },
+
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[2, 0],
+        [0, 2]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 128
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 128
+    }
+  },
+
+  "CT-AUX": {
+    "expansion-parameter-K": 1.0,
+    "initial-configuration-size": 10,
+    "max-submatrix-size": 16
+  },
+
+  "DCA": {
+    "iterations": 1,
+    "self-energy-mixing-factor": 0.0,
+    "interacting-orbitals": [0]
+  },
+
+  "Monte-Carlo-integration": {
+    "warm-up-sweeps": 200,
+    "sweeps-per-measurement": 1,
+    "measurements": 2000,
+    "seed": 42
+  },
+
+  "disorder": {
+    "distribution": "binary",
+    "potential": 0.1,
+    "density": 0.5,
+    "num-configurations": 4,
+    "unique-configs": false
+  }
+}
diff --git a/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_box.json b/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_box.json
new file mode 100644
index 000000000..5bf472ca1
--- /dev/null
+++ b/test/integration/statistical_tests/square_lattice/input_disorder_finite_v_box.json
@@ -0,0 +1,51 @@
+{
+  "physics": {
+    "beta"                      : 1.0,
+    "chemical-potential"        : 0.0
+  },
+
+  "single-band-Hubbard-model": {
+    "t"       : 1.0,
+    "U"       : 4.0
+  },
+
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[2, 0],
+        [0, 2]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 128
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 128
+    }
+  },
+
+  "CT-AUX": {
+    "expansion-parameter-K": 1.0,
+    "initial-configuration-size": 10,
+    "max-submatrix-size": 16
+  },
+
+  "DCA": {
+    "iterations": 1,
+    "self-energy-mixing-factor": 0.0,
+    "interacting-orbitals": [0]
+  },
+
+  "Monte-Carlo-integration": {
+    "warm-up-sweeps": 200,
+    "sweeps-per-measurement": 1,
+    "measurements": 2000,
+    "seed": 42
+  },
+
+  "disorder": {
+    "distribution": "box",
+    "potential": 0.1,
+    "density": 0.5,
+    "num-configurations": 4,
+    "unique-configs": false
+  }
+}
diff --git a/test/unit/math/function_transform/function_transform_test.cpp b/test/unit/math/function_transform/function_transform_test.cpp
index 6a1f68dba..5bab7ae53 100644
--- a/test/unit/math/function_transform/function_transform_test.cpp
+++ b/test/unit/math/function_transform/function_transform_test.cpp
@@ -31,15 +31,16 @@ using McOptions = MockMcOptions<double>;
 #include "dca/parallel/no_concurrency/no_concurrency.hpp"
 #include "dca/parallel/no_threading/no_threading.hpp"
 #include "dca/profiling/null_profiler.hpp"
+#include "dca/phys/domains/domain_aliases.hpp"
 
 using Model =
     dca::phys::models::TightBindingModel<dca::phys::models::RashbaHubbard<dca::phys::domains::D4>>;
 using Concurrency = dca::parallel::NoConcurrency;
-using Parameters =
-    dca::phys::params::Parameters<Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler,
-                                  Model, void, dca::ClusterSolverId::CT_AUX,               dca::NumericalTraits<double,
-                                   typename dca::util::ScalarSelect<double,
-                                                                    Model::lattice_type::complex_g0>::type>>;
+using Parameters = dca::phys::params::Parameters<
+    Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler, Model, void,
+    dca::ClusterSolverId::CT_AUX,
+    dca::NumericalTraits<
+        double, typename dca::util::ScalarSelect<double, Model::lattice_type::complex_g0>::type>>;
 
 const std::string input_dir = DCA_SOURCE_DIR "/test/unit/math/function_transform/";
 
@@ -49,6 +50,8 @@ using NuDmn = dca::func::dmn_variadic<BDmn, SDmn>;
 using WDmn = dca::func::dmn_0<dca::phys::domains::frequency_domain>;
 using KDmn = Parameters::KClusterDmn;
 using RDmn = Parameters::RClusterDmn;
+using DDA = dca::phys::DcaDomainAliases<Parameters>;
+using TDmn = typename DDA::TDmn;
 
 const std::vector<std::vector<double>> a_vecs{std::vector<double>{0, 0},
                                               std::vector<double>{0.25, 0.25}};
@@ -141,6 +144,80 @@ void spTestImplementation(const bool direct) {
   }
 }
 
+template <class InpDmn, class OutDmn>
+void timeFreqTestImplementation(const bool direct) {
+  using namespace dca::func;
+  using Real = double;
+  using Complex = std::complex<Real>;
+  function<Complex, dmn_variadic<NuDmn, NuDmn, KDmn, InpDmn>> f_b_b_in;
+  function<Complex, dmn_variadic<KDmn, InpDmn>> f_in;
+
+  // Initialize the input function.
+  for (int k = 0; k < KDmn::dmn_size(); ++k)
+    for (int i = 0; i < InpDmn::dmn_size(); ++i) {
+      const Complex val(i * i + 0.5 * i, k);
+      f_in(i, k) = val;
+      for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2)
+        for (int b1 = 0; b1 < BDmn::dmn_size(); ++b1)
+          f_b_b_in(b1, 0, b2, 0, k, i) = f_b_b_in(b1, 1, b2, 1, k, i) = val;
+    }
+
+  function<Complex, dmn_variadic<NuDmn, NuDmn, KDmn, OutDmn>> f_b_b_out;
+  function<Complex, dmn_variadic<KDmn, OutDmn>> f_out;
+
+  // Regular transform.
+  dca::math::transform::FunctionTransform<InpDmn, OutDmn>::execute(f_in, f_out);
+  // Transform with phase factor.
+  dca::math::transform::FunctionTransform<InpDmn, OutDmn>::execute(f_b_b_in, f_b_b_out);
+
+  const int exp_sign = direct ? 1 : -1;
+  const double norm = direct ? 1 : 1. / InpDmn::dmn_size();
+
+  constexpr double tolerance = 1e-10;
+
+  for (int k = 0; k < KDmn::dmn_size(); ++k) {
+    // Test regular transform.
+    for (int o = 0; o < OutDmn::dmn_size(); ++o) {
+      const auto& o_dmn_val = OutDmn::get_elements()[o];
+      Complex val(0);
+
+      for (int i = 0; i < InpDmn::dmn_size(); ++i) {
+        const auto& i_dmn_val = InpDmn::get_elements()[i];
+        val += f_in(i, k) * std::exp(Complex(0., exp_sign * dca::math::util::innerProduct(
+                                                                static_cast<double>(i_dmn_val),
+                                                                static_cast<double>(o_dmn_val))));
+      }
+      val *= norm;
+
+      EXPECT_LE(std::abs(val - f_out(o, k)), tolerance);
+
+      // // Test transform with phase factors.
+      // for (int b1 = 0; b1 < BDmn::dmn_size(); ++b1) {
+      //   for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2) {
+      //     const auto a_diff = dca::math::util::subtract(a_vecs[b2], a_vecs[b1]);
+      //     for (int o = 0; o < OutDmn::dmn_size(); ++o) {
+      //       const auto& o_dmn_val = OutDmn::get_elements()[o];
+      //       Complex val(0);
+
+      //       for (int i = 0; i < InpDmn::dmn_size(); ++i) {
+      //         const auto& i_dmn_val = InpDmn::get_elements()[i];
+      //         const auto& k_val = direct ? o_dmn_val : i_dmn_val;
+      //         const Complex phase_factor = std::exp(Complex(
+      //             0, exp_sign * dca::math::util::innerProduct(static_cast<double>(k_val), a_diff)));
+      //         val += f_b_b_in(b1, 0, b2, 0, k, i) * phase_factor *
+      //                std::exp(Complex(
+      //                    0., exp_sign * dca::math::util::innerProduct(i_dmn_val, o_dmn_val)));
+      //       }
+      //       val *= norm;
+
+      //       EXPECT_LE(std::abs(val - f_b_b_out(b1, 0, b2, 0, k, o)), tolerance);
+      //     }
+      //   }
+      // }
+    }
+  }
+}
+
 TEST(FunctionTransformTest, SpaceToMomentumCmplx) {
   initialize();
   spTestImplementation<RDmn, KDmn>(true);
@@ -150,3 +227,8 @@ TEST(FunctionTransformTest, MomentumToSpaceCmplx) {
   initialize();
   spTestImplementation<KDmn, RDmn>(false);
 }
+
+// TEST(FunctionTransformTest, TimeToFrequency) {
+//   initialize();
+//   timeFreqTestImplementation<TDmn, WDmn>(true);
+// }
diff --git a/test/unit/math/nfft/CMakeLists.txt b/test/unit/math/nfft/CMakeLists.txt
index 327923ba6..9ceee3f8f 100644
--- a/test/unit/math/nfft/CMakeLists.txt
+++ b/test/unit/math/nfft/CMakeLists.txt
@@ -12,4 +12,4 @@ dca_add_gtest(dnfft_1d_gpu_test
   GTEST_MAIN
   CUDA
   INCLUDE_DIRS ${PROJECT_SOURCE_DIR}
-  LIBS FFTW::Double time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS})
+  LIBS FFTW::Double quantum_domains time_and_frequency_domains random function mc_kernels gpu_utils nfft ${THIS_TEST_LIBS})
diff --git a/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp b/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
index 8f4e76978..9a019bd3b 100644
--- a/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
+++ b/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
@@ -177,8 +177,9 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorReal) {
   EXPECT_DOUBLE_EQ(f_copy(0), f(0));
   EXPECT_DOUBLE_EQ(f_copy(1), f(1));
 
-  EXPECT_DOUBLE_EQ(err_expected(0), err_no_overwriting(0));
-  EXPECT_DOUBLE_EQ(err_expected(1), err_no_overwriting(1));
+  // Jackknife error is a floating-point reduction: compare within epsilon_, not bit-exact.
+  EXPECT_NEAR(err_expected(0), err_no_overwriting(0), epsilon_);
+  EXPECT_NEAR(err_expected(1), err_no_overwriting(1), epsilon_);
 
   // Overwrite the jackknife estimates with their average.
   auto err_overwriting = sum_interface_.jackknifeError(f, true);
@@ -188,8 +189,8 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorReal) {
   EXPECT_DOUBLE_EQ(rank_avg, f(0));
   EXPECT_DOUBLE_EQ(d, f(1));
 
-  EXPECT_DOUBLE_EQ(err_expected(0), err_overwriting(0));
-  EXPECT_DOUBLE_EQ(err_expected(1), err_overwriting(1));
+  EXPECT_NEAR(err_expected(0), err_overwriting(0), epsilon_);
+  EXPECT_NEAR(err_expected(1), err_overwriting(1), epsilon_);
 }
 
 TEST_F(MPICollectiveSumTest, JackknifeErrorComplex) {
@@ -260,10 +261,11 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorComplex) {
   EXPECT_DOUBLE_EQ(f_copy(0).real(), f(0).real());
   EXPECT_DOUBLE_EQ(f_copy(1).real(), f(1).real());
 
-  EXPECT_DOUBLE_EQ(err_expected(0).real(), err_no_overwriting(0).real());
-  EXPECT_DOUBLE_EQ(err_expected(0).imag(), err_no_overwriting(0).imag());
-  EXPECT_DOUBLE_EQ(err_expected(1).real(), err_no_overwriting(1).real());
-  EXPECT_DOUBLE_EQ(err_expected(1).imag(), err_no_overwriting(1).imag());
+  // Jackknife error is a floating-point reduction: compare within epsilon_, not bit-exact.
+  EXPECT_NEAR(err_expected(0).real(), err_no_overwriting(0).real(), epsilon_);
+  EXPECT_NEAR(err_expected(0).imag(), err_no_overwriting(0).imag(), epsilon_);
+  EXPECT_NEAR(err_expected(1).real(), err_no_overwriting(1).real(), epsilon_);
+  EXPECT_NEAR(err_expected(1).imag(), err_no_overwriting(1).imag(), epsilon_);
 
   // Overwrite the jackknife estimates with their average.
   auto err_overwriting = sum_interface_.jackknifeError(f, true);
@@ -273,10 +275,10 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorComplex) {
   EXPECT_DOUBLE_EQ(std::complex<double>(rank_avg, rank_avg + r).real(), f(0).real());
   EXPECT_DOUBLE_EQ(std::complex<double>(rank_avg, d).real(), f(1).real());
 
-  EXPECT_DOUBLE_EQ(err_expected(0).real(), err_overwriting(0).real());
-  EXPECT_DOUBLE_EQ(err_expected(0).imag(), err_overwriting(0).imag());
-  EXPECT_DOUBLE_EQ(err_expected(1).real(), err_overwriting(1).real());
-  EXPECT_DOUBLE_EQ(err_expected(1).imag(), err_overwriting(1).imag());
+  EXPECT_NEAR(err_expected(0).real(), err_overwriting(0).real(), epsilon_);
+  EXPECT_NEAR(err_expected(0).imag(), err_overwriting(0).imag(), epsilon_);
+  EXPECT_NEAR(err_expected(1).real(), err_overwriting(1).real(), epsilon_);
+  EXPECT_NEAR(err_expected(1).imag(), err_overwriting(1).imag(), epsilon_);
 }
 
 TEST_F(MPICollectiveSumTest, ComputeCovarianceScalar) {
diff --git a/test/unit/phys/CMakeLists.txt b/test/unit/phys/CMakeLists.txt
index 9d64b1f5e..a6ab19113 100644
--- a/test/unit/phys/CMakeLists.txt
+++ b/test/unit/phys/CMakeLists.txt
@@ -6,3 +6,4 @@ add_subdirectory(dca_step)
 add_subdirectory(domains)
 add_subdirectory(models)
 add_subdirectory(parameters)
+add_subdirectory(dca_loop)
diff --git a/test/unit/phys/dca_data/CMakeLists.txt b/test/unit/phys/dca_data/CMakeLists.txt
index e27e5959c..70ae89e4d 100644
--- a/test/unit/phys/dca_data/CMakeLists.txt
+++ b/test/unit/phys/dca_data/CMakeLists.txt
@@ -4,3 +4,30 @@ dca_add_gtest(dca_data_test
   FAST
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS function ${DCA_LIBS} parallel_no_concurrency)
+
+# Disorder (DISORDERED_G0) tests are only built with the disorder option, never for a GPU build.
+if(DCA_WITH_DISORDER AND NOT DCA_WITH_CUDA AND NOT DCA_WITH_HIP)
+  # Tests the real-space (DISORDERED_G0) disordered-G0 construction: k->r FT of the clean cluster-
+  # excluded G0, the two space-index unfolding, and the per-frequency disorder Dyson inversions.
+  dca_add_gtest(make_disordered_g0_test
+    FAST
+    GTEST_MAIN
+    INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
+    LIBS function FFTW::Double ${DCA_LIBS} parallel_no_concurrency)
+
+  # Tests the DISORDERED_G0 walker Green's-function path: the dense disordered tau G0 (w->t Fourier
+  # with tail handling), the two-site walker akima storage, and absolute-site keying in build_G0_matrix.
+  dca_add_gtest(disordered_walker_g0_test
+    FAST
+    GTEST_MAIN
+    INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
+    LIBS function FFTW::Double ${DCA_LIBS} parallel_no_concurrency)
+
+  # Step F (the U != 0 teeth): drives the real CT-AUX walker at U != 0 and checks that the accumulated
+  # two-site M responds to the disorder potential, i.e. the walker samples the disordered G0.
+  dca_add_gtest(disordered_walker_u_test
+    FAST
+    GTEST_MAIN
+    INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
+    LIBS function FFTW::Double ${DCA_LIBS} ${DCA_KERNEL_LIBS} parallel_no_concurrency)
+endif()
diff --git a/test/unit/phys/dca_data/disordered_walker_g0_test.cpp b/test/unit/phys/dca_data/disordered_walker_g0_test.cpp
new file mode 100644
index 000000000..0c6a92fba
--- /dev/null
+++ b/test/unit/phys/dca_data/disordered_walker_g0_test.cpp
@@ -0,0 +1,217 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Unit tests for the DISORDERED_G0 walker Green's-function path, covering three chunks:
+//   A. DcaData::makeDisorderedG0 builds the dense disordered G0 in imaginary time via the w->t
+//      subtract-the-clean-reference Fourier transform (replaces the old per-slice site-diagonal
+//      loop). At zero potential the result must reproduce the clean tau G0 unfolded to two indices.
+//   B. The CT-AUX walker's G0Interpolation stores the now two-site (R0,R1) disordered G0 in its
+//      akima splines: the spline value at each tau node must equal the input G0.
+//   C. build_G0_matrix keys the lookups by the two ABSOLUTE sites, not a single displacement: two
+//      vertex pairs sharing a displacement but at different absolute sites must give different G0.
+
+#include <complex>
+#include <memory>
+#include <vector>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/io/json/json_reader.hpp"
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+#include "dca/linalg/matrix.hpp"
+#include "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp"
+#include "dca/phys/parameters/parameters.hpp"
+#include "dca/phys/models/analytic_hamiltonians/square_lattice.hpp"
+#include "dca/parallel/no_concurrency/no_concurrency.hpp"
+#include "dca/parallel/no_threading/no_threading.hpp"
+#include "dca/profiling/null_profiler.hpp"
+#include "dca/phys/dca_data/dca_data.hpp"
+// vertex_singleton must precede the g0 interpolation header, which uses it without including it.
+#include "dca/phys/dca_step/cluster_solver/ctaux/structs/vertex_singleton.hpp"
+#include "dca/phys/dca_step/cluster_solver/ctaux/walker/tools/g0_interpolation/g0_interpolation_cpu.hpp"
+#include "test/unit/phys/dca_step/cluster_solver/stub_rng.hpp"
+
+namespace {
+
+using Scalar = double;
+using Concurrency = dca::parallel::NoConcurrency;
+using Lattice = dca::phys::models::square_lattice<dca::phys::domains::D4>;
+using Model = dca::phys::models::TightBindingModel<Lattice>;
+using Parameters =
+    dca::phys::params::Parameters<Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler,
+                                  Model, dca::testing::StubRng, dca::ClusterSolverId::CT_AUX,
+                                  dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>;
+using Data = dca::phys::DcaData<Parameters>;
+
+using NuDmn = Data::NuDmn;
+using RDmn = Data::RClusterDmn;
+using KDmn = Data::KClusterDmn;
+using WDmn = Data::WDmn;
+using TDmn = Data::TDmn;
+using Complex = std::complex<double>;
+
+using G0Interp = dca::phys::solver::ctaux::G0Interpolation<dca::linalg::CPU, Parameters>;
+using Vertex = dca::phys::solver::ctaux::vertex_singleton;
+
+constexpr char input_file[] = DCA_SOURCE_DIR "/test/unit/phys/dca_data/input_make_disorder.json";
+
+// One set of globally initialized domains / parameters shared by all tests.
+Parameters& sharedParameters() {
+  static Concurrency concurrency(0, nullptr);
+  static Parameters parameters("", concurrency);
+  static bool initialized = false;
+  if (!initialized) {
+    parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+    parameters.update_model();
+    parameters.update_domains();
+    initialized = true;
+  }
+  return parameters;
+}
+
+std::unique_ptr<Data> makeData() {
+  return std::make_unique<Data>(sharedParameters());
+}
+
+// Clean, k-diagonal, cluster-excluded G0 in frequency: invertible (diagonal-dominant) nu-blocks,
+// distinct per k so the displacement structure is nontrivial. Drives makeDisorderedG0's freq path.
+void fillCleanG0kw(Data& data) {
+  auto& g0 = data.G0_k_w_cluster_excluded;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int k = 0; k < KDmn::dmn_size(); ++k)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          g0(n1, n2, k, w) = (n1 == n2)
+                                 ? Complex(3.0 + 0.1 * k + 0.01 * n1, 0.2 + 0.02 * w)
+                                 : Complex(0.3, n1 < n2 ? 0.1 : -0.1);
+}
+
+// Clean cluster-excluded G0 in (single-displacement) imaginary time: the tau reference that the
+// w->t step adds back. Structured in (nu, r, t) so the unfold/add-back has teeth.
+void fillCleanG0rt(Data& data) {
+  auto& g = data.G0_r_t_cluster_excluded;
+  for (int t = 0; t < TDmn::dmn_size(); ++t)
+    for (int r = 0; r < RDmn::dmn_size(); ++r)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          g(n1, n2, r, t) = 0.5 + 0.3 * r + 0.05 * t + 0.1 * (n1 == n2 ? 1.0 : -0.5);
+}
+
+// A synthetic two-site disordered G0 in tau that depends on BOTH absolute sites (R0 and R1), not
+// just their displacement -- the property the walker rewrite must respect.
+void fillSyntheticDisorderedG0rt(Data& data) {
+  auto& g = data.disordered_G0_r_r_t_cl_exl;
+  for (int t = 0; t < TDmn::dmn_size(); ++t)
+    for (int r1 = 0; r1 < RDmn::dmn_size(); ++r1)
+      for (int r0 = 0; r0 < RDmn::dmn_size(); ++r0)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          for (int n0 = 0; n0 < NuDmn::dmn_size(); ++n0)
+            g(n0, r0, n1, r1, t) =
+                (1.0 + 0.5 * r0 + 0.25 * r1 + 0.1 * (n0 + n1)) * (1.0 + 0.1 * t);
+}
+
+}  // namespace
+
+// Chunk A: at zero potential the dense disordered tau G0 must equal the clean tau G0 unfolded with
+// the same subtract(R2,R1) displacement (delta == 0, so only the added-back reference survives).
+TEST(DisorderedWalkerG0Test, ImaginaryTimeUnfoldAtZeroPotential) {
+  auto data = makeData();
+  fillCleanG0kw(*data);
+  fillCleanG0rt(*data);
+
+  typename Data::DisorderConfiguration config;  // zero potential
+  data->makeDisorderedG0(config);
+
+  const int n_r = RDmn::dmn_size();
+  constexpr double tol = 1e-8;
+  for (int t = 0; t < TDmn::dmn_size(); ++t)
+    for (int R1 = 0; R1 < n_r; ++R1)
+      for (int R2 = 0; R2 < n_r; ++R2) {
+        const int d = RDmn::parameter_type::subtract(R2, R1);
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2) {
+            const double got = data->disordered_G0_r_r_t_cl_exl(n1, R1, n2, R2, t);
+            const double ref = data->G0_r_t_cluster_excluded(n1, n2, d, t);
+            EXPECT_NEAR(got, ref, tol)
+                << "t=" << t << " R1=" << R1 << " R2=" << R2 << " n1=" << n1 << " n2=" << n2;
+          }
+      }
+}
+
+// Chunk B: the walker akima storage is now keyed by the two sites (R0,R1). The cubic's constant
+// term (value at the left tau node) must equal the input disordered G0 at that node, per site pair.
+// The two-site akima storage only exists under DISORDERED_G0; in the clean build the storage is
+// keyed by a single displacement, so this property is not defined.
+TEST(DisorderedWalkerG0Test, WalkerAkimaStoresTwoSiteG0) {
+  auto data = makeData();
+  fillSyntheticDisorderedG0rt(*data);
+
+  G0Interp g0(0, sharedParameters());
+  g0.initialize(*data);
+
+  const auto& ak = g0.getAkimaCoefficients();
+  const int n_r = RDmn::dmn_size();
+  const int half = TDmn::dmn_size() / 2;
+  constexpr double tol = 1e-9;
+
+  // Negative-tau half: shifted node s in [0, half-2] holds the spline for [node s, node s+1], and
+  // its constant coefficient a[0] equals the data value at full-tau index s.
+  for (int s = 0; s < half - 1; ++s)
+    for (int r1 = 0; r1 < n_r; ++r1)
+      for (int r0 = 0; r0 < n_r; ++r0)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          for (int n0 = 0; n0 < NuDmn::dmn_size(); ++n0) {
+            const double a0 = ak(0, n0, r0, n1, r1, s);
+            const double ref = data->disordered_G0_r_r_t_cl_exl(n0, r0, n1, r1, s);
+            EXPECT_NEAR(a0, ref, tol) << "s=" << s << " r0=" << r0 << " r1=" << r1;
+          }
+}
+
+// Chunk C: build_G0_matrix must key by absolute sites. Two vertex pairs with the SAME displacement
+// but at different absolute sites must yield different G0 (off-diagonal and site-resolved diagonal).
+// Site-resolved G0 only exists under DISORDERED_G0; the clean build collapses both pairs to their
+// shared displacement (giving identical G0), so this property is not defined there.
+TEST(DisorderedWalkerG0Test, BuildG0MatrixUsesAbsoluteSites) {
+  auto data = makeData();
+  fillSyntheticDisorderedG0rt(*data);
+
+  G0Interp g0(0, sharedParameters());
+  g0.initialize(*data);
+
+  const double beta = sharedParameters().get_beta();
+  const double tau_a = 0.3 * beta;
+  const double tau_b = 0.1 * beta;
+
+  // Two pairs sharing a displacement but at different absolute sites (4-site chain cluster).
+  const int rA0 = 0, rA1 = 1, rB0 = 1, rB1 = 2;
+  ASSERT_EQ(RDmn::parameter_type::subtract(rA1, rA0), RDmn::parameter_type::subtract(rB1, rB0))
+      << "test setup expects the two pairs to share a displacement";
+
+  auto vtx = [](int nu, int r, double tau) {
+    return Vertex(0, dca::phys::e_UP, nu, 0, r, 0, tau, dca::phys::solver::ctaux::HS_UP,
+                  dca::phys::solver::ctaux::HS_FIELD_UP, 0);
+  };
+  std::vector<Vertex> cfg1 = {vtx(0, rA0, tau_a), vtx(0, rA1, tau_b)};
+  std::vector<Vertex> cfg2 = {vtx(0, rB0, tau_a), vtx(0, rB1, tau_b)};
+
+  dca::linalg::Matrix<Scalar, dca::linalg::CPU> G0_1, G0_2;
+  g0.build_G0_matrix(cfg1, G0_1);
+  g0.build_G0_matrix(cfg2, G0_2);
+
+  // Off-diagonal: same displacement, different absolute sites -> different G0.
+  EXPECT_GT(std::abs(G0_1(0, 1) - G0_2(0, 1)), 1e-6)
+      << "build_G0_matrix collapsed to a displacement; disorder requires absolute sites";
+  // Diagonal is site-resolved (on-site G0_dis(R,R) varies site to site).
+  EXPECT_GT(std::abs(G0_1(0, 0) - G0_2(0, 0)), 1e-6);
+
+  // Sanity: rebuilding the same configuration is deterministic.
+  dca::linalg::Matrix<Scalar, dca::linalg::CPU> G0_1b;
+  g0.build_G0_matrix(cfg1, G0_1b);
+  EXPECT_DOUBLE_EQ(G0_1(0, 1), G0_1b(0, 1));
+}
diff --git a/test/unit/phys/dca_data/disordered_walker_u_test.cpp b/test/unit/phys/dca_data/disordered_walker_u_test.cpp
new file mode 100644
index 000000000..b9a4a1b5e
--- /dev/null
+++ b/test/unit/phys/dca_data/disordered_walker_u_test.cpp
@@ -0,0 +1,161 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Step F (the U != 0 teeth) for the DISORDERED_G0 walker path.
+//
+// All other disorder tests run at U = 0, where M = 0 and the reconstructed G = G0_dis regardless of
+// what the walker did -- so they cannot see whether the QMC walker actually samples in the disordered
+// medium. This test drives the real CT-AUX walker + accumulator at U != 0 and inspects the sampled
+// two-site M(R1,R2,w) directly:
+//
+//   * Determinism: same disorder config + same seed  ->  bitwise-identical M.
+//   * Teeth:       different disorder potential V     ->  M MUST change.
+//
+// The teeth isolates the walker. Comparing the reconstructed local_G_k_w instead would NOT work:
+// it differs with V through the Dyson's G0_dis even if the walker ignored V. M, by contrast, is
+// produced solely by the walker's sampling, so a V-dependent M proves the walker's bare line is the
+// disordered G0. Under the pre-fix bug (walker read the clean G0_r_t_cluster_excluded; V entered only
+// the post-QMC Dyson) M would be V-independent and the sensitivity assertion would fail.
+
+#include <algorithm>
+#include <complex>
+#include <memory>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/io/json/json_reader.hpp"
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+#include "dca/math/random/std_random_wrapper.hpp"
+#include "dca/parallel/no_concurrency/no_concurrency.hpp"
+#include "dca/parallel/no_threading/no_threading.hpp"
+#include "dca/profiling/null_profiler.hpp"
+#include "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp"
+#include "dca/phys/models/analytic_hamiltonians/square_lattice.hpp"
+#include "dca/phys/models/tight_binding_model.hpp"
+#include "dca/phys/parameters/parameters.hpp"
+
+using Scalar = double;
+
+#include "test/mock_mcconfig.hpp"
+namespace dca {
+namespace config {
+using McOptions = MockMcOptions<Scalar>;
+}  // namespace config
+}  // namespace dca
+
+#include "dca/phys/dca_data/dca_data.hpp"
+#include "dca/phys/dca_step/cluster_solver/ctaux/ctaux_cluster_solver.hpp"
+
+namespace {
+
+using Concurrency = dca::parallel::NoConcurrency;
+using Lattice = dca::phys::models::square_lattice<dca::phys::domains::D4>;
+using Model = dca::phys::models::TightBindingModel<Lattice>;
+using RngType = dca::math::random::StdRandomWrapper<std::mt19937_64>;
+using Parameters =
+    dca::phys::params::Parameters<Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler,
+                                  Model, RngType, dca::ClusterSolverId::CT_AUX,
+                                  dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>;
+using Data = dca::phys::DcaData<Parameters>;
+using Solver = dca::phys::solver::CtauxClusterSolver<dca::linalg::CPU, Parameters, Data>;
+using Walker = Solver::Walker;
+using Accumulator = Solver::Accumulator;
+
+using NuDmn = Data::NuDmn;
+using RDmn = Data::RClusterDmn;
+using WDmn = Data::WDmn;
+
+constexpr char input_file[] = DCA_SOURCE_DIR "/test/unit/phys/dca_data/input_disorder_walker_u.json";
+
+// Domain of the accumulated two-site M(nu0,R0,nu1,R1,w) (mirrors the solver's NuRNuRClusterWDmn).
+using MDmn = dca::func::dmn_variadic<NuDmn, RDmn, NuDmn, RDmn, WDmn>;
+using MFunction = dca::func::function<std::complex<double>, MDmn>;
+
+double maxAbs(const MFunction& f) {
+  double m = 0.;
+  for (int i = 0; i < f.size(); ++i)
+    m = std::max(m, std::abs(f(i)));
+  return m;
+}
+
+double maxAbsDiff(const MFunction& a, const MFunction& b) {
+  double m = 0.;
+  for (int i = 0; i < a.size(); ++i)
+    m = std::max(m, std::abs(a(i) - b(i)));
+  return m;
+}
+
+}  // namespace
+
+TEST(DisorderedWalkerUTest, WalkerMSamplesDisorderedMedium) {
+  Concurrency concurrency(0, nullptr);
+  Parameters parameters("", concurrency);
+  parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+  parameters.update_model();
+  parameters.update_domains();
+
+  Data data(parameters);
+  data.initialize();  // populates a valid bare cluster-excluded G0 from the model H_DCA.
+
+  // Drive the real CT-AUX walker + accumulator for one disorder config and return the accumulated,
+  // sign-weighted two-site M the walker sampled. Fresh walker/accumulator/rng each call (same seed),
+  // so two calls with the same config are bit-for-bit identical.
+  auto sampleM = [&](const typename Data::DisorderConfiguration& config) {
+    data.makeDisorderedG0(config);  // sets disordered_G0_r_r_{w,t}_cl_exl; the walker reads the tau one.
+
+    // StdRandomWrapper derives its engine seed from a static per-construction counter
+    // (engine_seed = hash(counter_ + seedArg)), so a fixed seed argument alone does NOT reproduce a
+    // stream. Reset the counter so every run draws the IDENTICAL random sequence; then the only thing
+    // that can move M between runs is the disorder G0 the walker samples in.
+    RngType::resetCounter();
+    RngType rng(0, 1, parameters.get_seed());
+    Walker walker(parameters, data, rng, 0);
+    Accumulator accumulator(parameters, data, 0);
+
+    walker.initialize(0);
+    accumulator.initialize(0);
+
+    for (int i = 0; i < parameters.get_warm_up_sweeps(); ++i)
+      walker.doSweep();
+    walker.markThermalized();
+
+    const int n_meas = parameters.get_measurements()[0];
+    for (int i = 0; i < n_meas; ++i) {
+      walker.doSweep();
+      accumulator.updateFrom(walker);
+      accumulator.measure();
+    }
+    accumulator.finalize();
+
+    return MFunction(accumulator.get_sign_times_M_r_w());
+  };
+
+  typename Data::DisorderConfiguration zero_config;  // V = 0 everywhere.
+
+  typename Data::DisorderConfiguration nonzero_config;  // staggered site potential, breaks symmetry.
+  for (int r = 0; r < RDmn::dmn_size(); ++r)
+    for (int nu = 0; nu < NuDmn::dmn_size(); ++nu)
+      nonzero_config(nu, r) = (r % 2 == 0) ? 0.6 : -0.4;
+
+  const MFunction M_zero_a = sampleM(zero_config);
+  const MFunction M_zero_b = sampleM(zero_config);
+  const MFunction M_nonzero = sampleM(nonzero_config);
+
+  // Sanity: U != 0, so the walker actually sampled a nontrivial M.
+  EXPECT_GT(maxAbs(M_zero_a), 1e-6) << "expected a nontrivial interacting M at U != 0";
+
+  // Determinism anchor: identical config + identical seed must reproduce M exactly.
+  EXPECT_LT(maxAbsDiff(M_zero_a, M_zero_b), 1e-12)
+      << "same disorder config and seed must give identical M";
+
+  // The teeth: the walker's sampled M must respond to the disorder potential, i.e. the walker
+  // propagates in the disordered G0. Under the old bug (walker in the clean G0) this difference is 0.
+  EXPECT_GT(maxAbsDiff(M_zero_a, M_nonzero), 1e-6)
+      << "walker M did not change with V: the walker is not sampling the disordered G0";
+}
diff --git a/test/unit/phys/dca_data/input_disorder_walker_u.json b/test/unit/phys/dca_data/input_disorder_walker_u.json
new file mode 100644
index 000000000..ecf7b2766
--- /dev/null
+++ b/test/unit/phys/dca_data/input_disorder_walker_u.json
@@ -0,0 +1,40 @@
+{
+  "physics": {
+    "beta": 2,
+    "chemical-potential": 0
+  },
+  "single-band-Hubbard-model": {
+    "t": 1,
+    "U": 2
+  },
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[4, 0],
+        [0, 1]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 64
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 16,
+      "four-point-fermionic-frequencies": 8
+    }
+  },
+  "four-point": {
+    "type": "NONE",
+    "momentum-transfer": [0., 0.],
+    "frequency-transfer": 0
+  },
+  "Monte-Carlo-integration": {
+    "seed": 42,
+    "warm-up-sweeps": 20,
+    "sweeps-per-measurement": 1,
+    "measurements": 150
+  },
+  "CT-AUX": {
+    "expansion-parameter-K": 1.0,
+    "initial-configuration-size": 10,
+    "initial-matrix-size": 64,
+    "max-submatrix-size": 1
+  }
+}
diff --git a/test/unit/phys/dca_data/input_make_disorder.json b/test/unit/phys/dca_data/input_make_disorder.json
new file mode 100644
index 000000000..02a81250f
--- /dev/null
+++ b/test/unit/phys/dca_data/input_make_disorder.json
@@ -0,0 +1,28 @@
+{
+  "physics": {
+    "beta": 2,
+    "chemical-potential": 0
+  },
+  "single-band-Hubbard-model": {
+    "t": 1,
+    "U": 2
+  },
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[4, 0],
+        [0, 1]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 64
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 16,
+      "four-point-fermionic-frequencies": 8
+    }
+  },
+  "four-point": {
+    "type": "NONE",
+    "momentum-transfer": [0., 0.],
+    "frequency-transfer": 0
+  }
+}
diff --git a/test/unit/phys/dca_data/make_disordered_g0_test.cpp b/test/unit/phys/dca_data/make_disordered_g0_test.cpp
new file mode 100644
index 000000000..8db07e509
--- /dev/null
+++ b/test/unit/phys/dca_data/make_disordered_g0_test.cpp
@@ -0,0 +1,222 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Unit tests for the real-space (R1,R2) disordered-G0 construction in DcaData::makeDisorderedG0
+// used by the DISORDERED_G0 path:
+//   1. the single k->r FunctionTransform of the cluster-excluded clean G0,
+//   2. unfolding the translationally invariant clean G0 into the full two space-index quantity
+//      disordered_G0_r_r_w_cl_exl(nu1,R1,nu2,R2,w) = g(R1 - R2) via the cluster subtract table,
+//   3. the per-frequency disorder Dyson in the dense (nu*N_R) basis (invert -> add a block ->
+//      invert).
+//
+// The first two tests drive the real makeDisorderedG0 (the disorder block it injects is currently
+// a zero dummy, so its output equals the clean unfolded G0 -- a clean inversion round trip). The
+// third test exercises the invert -> add nonzero block -> invert pipeline that the real disorder
+// potential will use, on a genuine clean block extracted from makeDisorderedG0's output.
+
+#include <complex>
+#include <memory>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/io/json/json_reader.hpp"
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+#include "dca/linalg/matrix.hpp"
+#include "dca/linalg/matrixop.hpp"
+#include "dca/math/util/vector_operations.hpp"
+#include "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp"
+
+#include "dca/phys/parameters/parameters.hpp"
+#include "dca/phys/models/analytic_hamiltonians/square_lattice.hpp"
+#include "dca/parallel/no_concurrency/no_concurrency.hpp"
+#include "dca/parallel/no_threading/no_threading.hpp"
+#include "dca/profiling/null_profiler.hpp"
+#include "dca/phys/dca_data/dca_data.hpp"
+#include "test/unit/phys/dca_step/cluster_solver/stub_rng.hpp"
+
+namespace {
+
+using Scalar = double;
+using Concurrency = dca::parallel::NoConcurrency;
+using Lattice = dca::phys::models::square_lattice<dca::phys::domains::D4>;
+using Model = dca::phys::models::TightBindingModel<Lattice>;
+using Parameters =
+    dca::phys::params::Parameters<Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler,
+                                  Model, dca::testing::StubRng, dca::ClusterSolverId::CT_AUX,
+                                  dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>;
+using Data = dca::phys::DcaData<Parameters>;
+
+using NuDmn = Data::NuDmn;
+using RDmn = Data::RClusterDmn;
+using KDmn = Data::KClusterDmn;
+using WDmn = Data::WDmn;
+using TDmn = Data::TDmn;
+using Complex = std::complex<double>;
+
+constexpr char input_file[] = DCA_SOURCE_DIR "/test/unit/phys/dca_data/input_make_disorder.json";
+
+// Builds a fresh DcaData sharing one set of (globally initialized) domains.
+std::unique_ptr<Data> makeData() {
+  static Concurrency concurrency(0, nullptr);
+  static Parameters parameters("", concurrency);
+  static bool initialized = false;
+  if (!initialized) {
+    parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+    parameters.update_model();
+    parameters.update_domains();
+    initialized = true;
+  }
+  return std::make_unique<Data>(parameters);
+}
+
+// Clean, k-diagonal, cluster-excluded G0 in frequency: diagonal-dominant (invertible) nu-blocks,
+// distinct per k so the displacement structure is nontrivial.
+void fillCleanG0(Data& data) {
+  auto& g0 = data.G0_k_w_cluster_excluded;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int k = 0; k < KDmn::dmn_size(); ++k)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          g0(n1, n2, k, w) = (n1 == n2)
+                                 ? Complex(3.0 + 0.1 * k + 0.01 * n1, 0.2 + 0.02 * w)
+                                 : Complex(0.3, n1 < n2 ? 0.1 : -0.1);
+}
+
+// The active (legacy) first loop in makeDisorderedG0 inverts nu-blocks of G0_r_t_cluster_excluded;
+// give it invertible (identity) blocks so it is well defined. Its output is unused here.
+void fillG0rtIdentity(Data& data) {
+  auto& g = data.G0_r_t_cluster_excluded;
+  g = Scalar(0);
+  for (int t = 0; t < TDmn::dmn_size(); ++t)
+    for (int r = 0; r < RDmn::dmn_size(); ++r)
+      for (int n = 0; n < NuDmn::dmn_size(); ++n)
+        g(n, n, r, t) = 1.0;
+}
+
+}  // namespace
+
+// makeDisorderedG0 should FunctionTransform the clean k-space G0 to (r,w) and unfold it to the full
+// two space-index quantity. With the (currently zero) disorder block, the invert/invert round trip
+// returns the clean values, so disordered_G0_r_r_w_cl_exl must equal the manual inverse-DFT of the
+// clean G0 unfolded via subtract(R2,R1).
+TEST(MakeDisorderedG0Test, KtoRFtAndUnfold) {
+  auto data = makeData();
+  fillCleanG0(*data);
+  fillG0rtIdentity(*data);
+
+  typename Data::DisorderConfiguration config;  // zero potential
+  data->makeDisorderedG0(config);
+
+  const int n_k = KDmn::dmn_size();
+  const int n_r = RDmn::dmn_size();
+  constexpr double tol = 1e-9;
+
+  // Manual clean g(d) = (1/N_k) sum_k G0_k exp(-i k.r_d) (FunctionTransform K->R convention).
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int R1 = 0; R1 < n_r; ++R1)
+      for (int R2 = 0; R2 < n_r; ++R2) {
+        const int d = RDmn::parameter_type::subtract(R2, R1);  // index of r_{R1} - r_{R2}
+        const auto& r_d = RDmn::get_elements()[d];
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2) {
+            Complex g_r(0);
+            for (int k = 0; k < n_k; ++k) {
+              const auto& k_val = KDmn::get_elements()[k];
+              g_r += data->G0_k_w_cluster_excluded(n1, n2, k, w) *
+                     std::exp(Complex(0, -dca::math::util::innerProduct(k_val, r_d)));
+            }
+            g_r /= static_cast<double>(n_k);
+
+            const Complex got = data->disordered_G0_r_r_w_cl_exl(n1, R1, n2, R2, w);
+            EXPECT_NEAR(got.real(), g_r.real(), tol)
+                << "mismatch at w=" << w << " (n1,R1,n2,R2)=(" << n1 << "," << R1 << "," << n2 << ","
+                << R2 << ")";
+            EXPECT_NEAR(got.imag(), g_r.imag(), tol);
+          }
+      }
+}
+
+// The unfolded clean G0 must be translationally invariant: every (R1,R2) pair sharing a
+// displacement subtract(R2,R1) carries the same nu-block.
+TEST(MakeDisorderedG0Test, UnfoldingIsTranslationallyInvariant) {
+  auto data = makeData();
+  fillCleanG0(*data);
+  fillG0rtIdentity(*data);
+
+  typename Data::DisorderConfiguration config;
+  data->makeDisorderedG0(config);
+
+  const int n_r = RDmn::dmn_size();
+  constexpr double tol = 1e-10;
+  const int w = 0;
+
+  for (int R1 = 0; R1 < n_r; ++R1)
+    for (int R2 = 0; R2 < n_r; ++R2) {
+      const int d = RDmn::parameter_type::subtract(R2, R1);
+      // Reference pair for this displacement: (R1'=d, R2'=0) also has subtract(0,d) == d.
+      for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+        for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2) {
+          const Complex a = data->disordered_G0_r_r_w_cl_exl(n1, R1, n2, R2, w);
+          const Complex b = data->disordered_G0_r_r_w_cl_exl(n1, d, n2, 0, w);
+          EXPECT_NEAR(a.real(), b.real(), tol);
+          EXPECT_NEAR(a.imag(), b.imag(), tol);
+        }
+    }
+}
+
+// Disorder Dyson on the dense (nu*N_R) basis: invert the clean block, add a nonzero disorder block
+// V, invert again. Verifies the algebraic identity inverse(G_dis) - inverse(G_clean) == V (so the
+// two inversions compose correctly through the contiguous block layout) and that V actually changes
+// the Green's function.
+TEST(MakeDisorderedG0Test, DisorderBlockInversionIdentity) {
+  auto data = makeData();
+  fillCleanG0(*data);
+  fillG0rtIdentity(*data);
+
+  typename Data::DisorderConfiguration config;
+  data->makeDisorderedG0(config);  // disordered_G0_r_r_w_cl_exl now holds the clean unfolded G0.
+
+  const int nu_r = NuDmn::dmn_size() * RDmn::dmn_size();
+  const int w = 0;
+
+  // Clean block for frequency w (same pointer idiom as makeDisorderedG0).
+  dca::linalg::Matrix<Complex, dca::linalg::CPU> clean(nu_r);
+  dca::linalg::matrixop::copyArrayToMatrix(
+      nu_r, nu_r, &data->disordered_G0_r_r_w_cl_exl(0, 0, 0, 0, w), nu_r, clean);
+
+  // Nonzero, modest disorder block (general, not site-diagonal).
+  dca::linalg::Matrix<Complex, dca::linalg::CPU> V(nu_r);
+  for (int j = 0; j < nu_r; ++j)
+    for (int i = 0; i < nu_r; ++i)
+      V(i, j) = Complex(0.05 * (i + 1) - 0.02 * j, 0.01 * (i - j));
+
+  dca::linalg::Matrix<Complex, dca::linalg::CPU> inv_clean(clean);
+  dca::linalg::matrixop::inverse(inv_clean);
+
+  dca::linalg::Matrix<Complex, dca::linalg::CPU> disordered(inv_clean);
+  for (int j = 0; j < nu_r; ++j)
+    for (int i = 0; i < nu_r; ++i)
+      disordered(i, j) += V(i, j);
+  dca::linalg::matrixop::inverse(disordered);
+
+  dca::linalg::Matrix<Complex, dca::linalg::CPU> inv_dis(disordered);
+  dca::linalg::matrixop::inverse(inv_dis);
+
+  constexpr double tol = 1e-8;
+  bool differs = false;
+  for (int j = 0; j < nu_r; ++j)
+    for (int i = 0; i < nu_r; ++i) {
+      const Complex recovered = inv_dis(i, j) - inv_clean(i, j);
+      EXPECT_NEAR(recovered.real(), V(i, j).real(), tol) << "(" << i << "," << j << ")";
+      EXPECT_NEAR(recovered.imag(), V(i, j).imag(), tol) << "(" << i << "," << j << ")";
+      if (std::abs(disordered(i, j) - clean(i, j)) > 1e-6)
+        differs = true;
+    }
+  EXPECT_TRUE(differs) << "the disorder block must change the Green's function";
+}
diff --git a/test/unit/phys/dca_loop/CMakeLists.txt b/test/unit/phys/dca_loop/CMakeLists.txt
new file mode 100644
index 000000000..f1e27fb4c
--- /dev/null
+++ b/test/unit/phys/dca_loop/CMakeLists.txt
@@ -0,0 +1,6 @@
+# test/unit/phys/dca_loop
+
+# Disorder unit tests are only built with the disorder option, and never for a GPU build.
+if(DCA_WITH_DISORDER AND NOT DCA_WITH_CUDA AND NOT DCA_WITH_HIP)
+  add_subdirectory(disorder)
+endif()
diff --git a/test/unit/phys/dca_loop/disorder/CMakeLists.txt b/test/unit/phys/dca_loop/disorder/CMakeLists.txt
new file mode 100644
index 000000000..cdd6af5fe
--- /dev/null
+++ b/test/unit/phys/dca_loop/disorder/CMakeLists.txt
@@ -0,0 +1,19 @@
+# test/unit/phys/dca_loop/disorder
+
+# DISORDERED_G0 numeric kernels: translationalAverage orientation round-trip (Nc=1, 2x2, 3x4)
+# and the per-frequency disorder Dyson (accumulateDisorderDyson). Uses mock domains, so it needs
+# no FFTW and no concurrency-initialized cluster domain.
+dca_add_gtest(disorder_average_test
+  FAST
+  GTEST_MAIN
+  INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
+  LIBS function ${DCA_LIBS} parallel_no_concurrency)
+
+# DISORDERED_G0 disorder-averaging pipeline (loop body of workTheClusters + averageGkw) on the real
+# DcaData with U=0 (M=0): clean-limit identity and the weighted disorder average. Uses the real
+# cluster (2x4) and the r<->k FunctionTransform, so it links FFTW.
+dca_add_gtest(disorder_loop_test
+  FAST
+  GTEST_MAIN
+  INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
+  LIBS function FFTW::Double ${DCA_LIBS} parallel_no_concurrency)
diff --git a/test/unit/phys/dca_loop/disorder/disorder_average_test.cpp b/test/unit/phys/dca_loop/disorder/disorder_average_test.cpp
new file mode 100644
index 000000000..60c18b1f7
--- /dev/null
+++ b/test/unit/phys/dca_loop/disorder/disorder_average_test.cpp
@@ -0,0 +1,263 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Unit tests for the DISORDERED_G0 numeric kernels in
+// include/dca/phys/dca_loop/disorder/disorder_average.hpp:
+//
+//   1. translationalAverage  - orientation / round-trip: unfolding a single-displacement g0 into
+//      the two space-index quantity (using subtract(R2,R1)) and then projecting it back (using
+//      subtract(R,R0)) must recover g0 exactly. Tested on three clusters: single site (Nc=1), a
+//      2x2 cluster (Nc=4), and an irregular 3x4 cluster (Nc=12). The first two are self-inverse
+//      (-r == r for every site) and CANNOT distinguish subtract(R,R0) from its negation; the 3x4
+//      cluster has sites with -r != r, so it is the one that actually catches an add-vs-subtract
+//      (orientation) flip in the projection.
+//
+//   3. accumulateDisorderDyson - the per-frequency dense (nu*N_R) Dyson
+//      G = G0_dis - G0_dis*M*G0_dis/beta, accumulated weighted. Checked against an independent
+//      naive matrix computation, and that repeated calls accumulate (+=) rather than overwrite.
+//
+// The tests use lightweight mock domains (a band/spin/frequency placeholder and a rectangular
+// cluster mock that carries its own subtract table) so all three cluster shapes coexist in one
+// executable - the real cluster_domain is a global singleton initialized once per binary.
+
+#include "dca/phys/dca_loop/disorder/disorder_average.hpp"
+
+#include <complex>
+#include <string>
+#include <vector>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+
+namespace {
+
+using dca::func::dmn_0;
+using dca::func::dmn_variadic;
+using Complex = std::complex<double>;
+
+// Minimal fixed-size leaf domain (band/spin, frequency placeholders). Satisfies the dmn_0
+// parameter_type interface: get_size / element_type + get_elements / get_name.
+template <int N>
+struct FixedDmn {
+  using element_type = int;
+  static int get_size() {
+    return N;
+  }
+  static int dmn_size() {
+    return N;
+  }
+  static std::vector<int>& get_elements() {
+    static std::vector<int> e = [] {
+      std::vector<int> v(N);
+      for (int i = 0; i < N; ++i)
+        v[i] = i;
+      return v;
+    }();
+    return e;
+  }
+  static std::string get_name() {
+    return "FixedDmn";
+  }
+};
+
+// Rectangular LX x LY cluster mock with periodic wrap. Site index = x + LX*y. Carries the cluster
+// subtract table with the production convention subtract(i,j) = index of r_j - r_i. For LX or LY
+// odd-or-not, sites with -r != r exist iff the cluster is not self-inverse (e.g. 3x4).
+template <int LX, int LY>
+struct RectCluster {
+  using element_type = int;
+  static int get_size() {
+    return LX * LY;
+  }
+  static int dmn_size() {
+    return LX * LY;
+  }
+  static std::vector<int>& get_elements() {
+    static std::vector<int> e = [] {
+      std::vector<int> v(LX * LY);
+      for (int i = 0; i < LX * LY; ++i)
+        v[i] = i;
+      return v;
+    }();
+    return e;
+  }
+  static std::string get_name() {
+    return "RectCluster";
+  }
+  // index of r_j - r_i, periodic in the LX x LY cell.
+  static int subtract(int i, int j) {
+    const int xi = i % LX, yi = i / LX;
+    const int xj = j % LX, yj = j / LX;
+    const int dx = ((xj - xi) % LX + LX) % LX;
+    const int dy = ((yj - yi) % LY + LY) % LY;
+    return dx + LX * dy;
+  }
+  // index of -r_i (used only by the test to reason about self-inverseness).
+  static int negate(int i) {
+    return subtract(i, 0);  // r_0 - r_i = -r_i
+  }
+};
+
+using NuMock = dmn_0<FixedDmn<2>>;  // 1 band x 2 spin placeholder
+using WMock = dmn_0<FixedDmn<3>>;   // 3 frequencies
+
+// A distinct, displacement-sensitive value for the reference single-displacement g0.
+Complex g0Value(int n1, int n2, int d, int w) {
+  return Complex(1.0 + d + 0.1 * n1 + 0.01 * n2 + 0.001 * w, 0.5 * d - 0.2 * w + 0.03 * n1);
+}
+
+// Builds g0, unfolds it to the two space-index quantity via subtract(R2,R1) (as makeDisorderedG0
+// does), runs translationalAverage, and asserts the projection recovers g0 exactly.
+template <int LX, int LY>
+void runTranslationalAverageRoundTrip() {
+  using RDmn = dmn_0<RectCluster<LX, LY>>;
+  using TwoR = dca::func::function<Complex, dmn_variadic<NuMock, RDmn, NuMock, RDmn, WMock>>;
+  using OneR = dca::func::function<Complex, dmn_variadic<NuMock, NuMock, RDmn, WMock>>;
+
+  const int N = RDmn::dmn_size();
+  const int nu = NuMock::dmn_size();
+  const int nw = WMock::dmn_size();
+
+  OneR g0("g0");
+  for (int w = 0; w < nw; ++w)
+    for (int d = 0; d < N; ++d)
+      for (int n2 = 0; n2 < nu; ++n2)
+        for (int n1 = 0; n1 < nu; ++n1)
+          g0(n1, n2, d, w) = g0Value(n1, n2, d, w);
+
+  // Unfold: GG(n1,R1,n2,R2,w) = g0(n1,n2, subtract(R2,R1), w).
+  TwoR gg("gg");
+  for (int w = 0; w < nw; ++w)
+    for (int R1 = 0; R1 < N; ++R1)
+      for (int R2 = 0; R2 < N; ++R2) {
+        const int d = RDmn::parameter_type::subtract(R2, R1);
+        for (int n2 = 0; n2 < nu; ++n2)
+          for (int n1 = 0; n1 < nu; ++n1)
+            gg(n1, R1, n2, R2, w) = g0(n1, n2, d, w);
+      }
+
+  OneR projected("projected");
+  dca::phys::disorder::translationalAverage<NuMock, RDmn, WMock>(gg, projected);
+
+  constexpr double tol = 1e-12;
+  for (int w = 0; w < nw; ++w)
+    for (int R = 0; R < N; ++R)
+      for (int n2 = 0; n2 < nu; ++n2)
+        for (int n1 = 0; n1 < nu; ++n1) {
+          const Complex got = projected(n1, n2, R, w);
+          const Complex ref = g0(n1, n2, R, w);
+          EXPECT_NEAR(got.real(), ref.real(), tol)
+              << "Nc=" << N << " (n1,n2,R,w)=(" << n1 << "," << n2 << "," << R << "," << w << ")";
+          EXPECT_NEAR(got.imag(), ref.imag(), tol);
+        }
+}
+
+}  // namespace
+
+TEST(DisorderAverageTest, TranslationalAverageRoundTripSingleSite) {
+  runTranslationalAverageRoundTrip<1, 1>();  // Nc=1
+}
+
+TEST(DisorderAverageTest, TranslationalAverageRoundTripFourSite) {
+  runTranslationalAverageRoundTrip<2, 2>();  // Nc=4, self-inverse
+}
+
+TEST(DisorderAverageTest, TranslationalAverageRoundTripIrregular3x4) {
+  // Nc=12: this cluster has sites with -r != r, so it (unlike the 1- and 4-site clusters) fails if
+  // the projection uses the wrong (negated) displacement orientation.
+  using RDmn = RectCluster<3, 4>;
+  bool has_non_self_inverse = false;
+  for (int i = 0; i < RDmn::get_size(); ++i)
+    if (RDmn::negate(i) != i)
+      has_non_self_inverse = true;
+  ASSERT_TRUE(has_non_self_inverse) << "the 3x4 cluster must not be self-inverse to be a teeth test";
+
+  runTranslationalAverageRoundTrip<3, 4>();
+}
+
+// accumulateDisorderDyson on a small (Nc=2) cluster: compare to an independent naive computation of
+// weight * (G0 - G0*M*G0/beta) on the dense (nu*N_R) block, and verify accumulation.
+TEST(DisorderAverageTest, DisorderDysonMatchesNaiveAndAccumulates) {
+  using RDmn = dmn_0<RectCluster<2, 1>>;  // Nc=2
+  using TwoR = dca::func::function<Complex, dmn_variadic<NuMock, RDmn, NuMock, RDmn, WMock>>;
+
+  const int N = RDmn::dmn_size();
+  const int nu = NuMock::dmn_size();
+  const int nw = WMock::dmn_size();
+  const int dim = nu * N;  // matrix dimension
+  const double beta = 2.5;
+  const double weight = 0.75;
+
+  // Controlled, invertible-ish G0_dis and a modest M.
+  TwoR g0("g0_dis"), m("m"), g_acc("g_acc");
+  auto block_index = [nu](int n, int R) { return n + nu * R; };  // row/col index in the dense block
+  for (int w = 0; w < nw; ++w)
+    for (int R2 = 0; R2 < N; ++R2)
+      for (int n2 = 0; n2 < nu; ++n2)
+        for (int R1 = 0; R1 < N; ++R1)
+          for (int n1 = 0; n1 < nu; ++n1) {
+            const int i = block_index(n1, R1), j = block_index(n2, R2);
+            g0(n1, R1, n2, R2, w) =
+                (i == j) ? Complex(2.0 + 0.1 * i + 0.05 * w, 0.3) : Complex(0.2 - 0.01 * (i + j), 0.1);
+            m(n1, R1, n2, R2, w) = Complex(0.4 - 0.02 * i + 0.03 * j, 0.05 * (i - j) + 0.01 * w);
+          }
+
+  g_acc = 0;
+  dca::phys::disorder::accumulateDisorderDyson<NuMock, RDmn, WMock, double>(g0, m, beta, weight, g_acc);
+
+  // Independent naive computation per frequency: result = weight*(G0 - G0*M*G0/beta).
+  constexpr double tol = 1e-11;
+  for (int w = 0; w < nw; ++w) {
+    std::vector<Complex> G0(dim * dim), M(dim * dim);
+    auto at = [dim](int i, int j) { return i + dim * j; };
+    for (int R2 = 0; R2 < N; ++R2)
+      for (int n2 = 0; n2 < nu; ++n2)
+        for (int R1 = 0; R1 < N; ++R1)
+          for (int n1 = 0; n1 < nu; ++n1) {
+            const int i = block_index(n1, R1), j = block_index(n2, R2);
+            G0[at(i, j)] = g0(n1, R1, n2, R2, w);
+            M[at(i, j)] = m(n1, R1, n2, R2, w);
+          }
+    std::vector<Complex> G0M(dim * dim, 0), G0MG0(dim * dim, 0);
+    for (int i = 0; i < dim; ++i)
+      for (int j = 0; j < dim; ++j) {
+        Complex s = 0;
+        for (int k = 0; k < dim; ++k)
+          s += G0[at(i, k)] * M[at(k, j)];
+        G0M[at(i, j)] = s;
+      }
+    for (int i = 0; i < dim; ++i)
+      for (int j = 0; j < dim; ++j) {
+        Complex s = 0;
+        for (int k = 0; k < dim; ++k)
+          s += G0M[at(i, k)] * G0[at(k, j)];
+        G0MG0[at(i, j)] = s;
+      }
+    for (int R2 = 0; R2 < N; ++R2)
+      for (int n2 = 0; n2 < nu; ++n2)
+        for (int R1 = 0; R1 < N; ++R1)
+          for (int n1 = 0; n1 < nu; ++n1) {
+            const int i = block_index(n1, R1), j = block_index(n2, R2);
+            const Complex ref = weight * (G0[at(i, j)] - G0MG0[at(i, j)] / beta);
+            const Complex got = g_acc(n1, R1, n2, R2, w);
+            EXPECT_NEAR(got.real(), ref.real(), tol) << "w=" << w << " (i,j)=(" << i << "," << j << ")";
+            EXPECT_NEAR(got.imag(), ref.imag(), tol);
+          }
+  }
+
+  // Accumulation: a second identical call must double the accumulator (+= not =).
+  TwoR after_one("after_one");
+  after_one = g_acc;
+  dca::phys::disorder::accumulateDisorderDyson<NuMock, RDmn, WMock, double>(g0, m, beta, weight, g_acc);
+  for (int idx = 0; idx < g_acc.size(); ++idx) {
+    EXPECT_NEAR(g_acc(idx).real(), 2.0 * after_one(idx).real(), tol);
+    EXPECT_NEAR(g_acc(idx).imag(), 2.0 * after_one(idx).imag(), tol);
+  }
+}
diff --git a/test/unit/phys/dca_loop/disorder/disorder_loop_test.cpp b/test/unit/phys/dca_loop/disorder/disorder_loop_test.cpp
new file mode 100644
index 000000000..90f39ac35
--- /dev/null
+++ b/test/unit/phys/dca_loop/disorder/disorder_loop_test.cpp
@@ -0,0 +1,200 @@
+// 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: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// Unit test for the DISORDERED_G0 disorder-averaging pipeline as wired in DcaLoop::workTheClusters
+// / averageGkw, exercised on the real DcaData (2x4 cluster, Nc=8) WITHOUT running the QMC.
+//
+// With U = 0 the interaction vertices vanish, so M = 0 and the per-configuration Dyson collapses to
+// G_c = G0_dis,c. We therefore drive the real pipeline kernels directly per configuration:
+//   makeDisorderedG0(config) -> accumulateDisorderDyson(G0_dis, M=0, ...) -> /total_weight
+//   -> translationalAverage -> FunctionTransform r->k,
+// exactly as the loop does, and check:
+//
+//   1. Clean limit (zero disorder potential): the whole machinery is the identity on G0 - the
+//      disorder-averaged G_k_w reproduces the clean cluster-excluded G0_k_w (the unfold, projection
+//      and r<->k transforms compose to identity).
+//   2. The weighted disorder average equals the explicit weighted mean of the per-configuration
+//      disordered G0 (validates the accumulate + normalize wiring), and the result is finite and
+//      actually shifted away from the clean G0 by the disorder.
+
+#include "dca/phys/dca_loop/disorder/disorder_average.hpp"
+
+#include <cmath>
+#include <complex>
+#include <memory>
+#include <numeric>
+#include <vector>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+#include "dca/io/json/json_reader.hpp"
+#include "dca/function/domains.hpp"
+#include "dca/function/function.hpp"
+#include "dca/math/function_transform/function_transform.hpp"
+#include "dca/phys/domains/cluster/symmetries/point_groups/2d/2d_square.hpp"
+
+#include "dca/phys/parameters/parameters.hpp"
+#include "dca/phys/models/analytic_hamiltonians/square_lattice.hpp"
+#include "dca/parallel/no_concurrency/no_concurrency.hpp"
+#include "dca/parallel/no_threading/no_threading.hpp"
+#include "dca/profiling/null_profiler.hpp"
+#include "dca/phys/dca_data/dca_data.hpp"
+#include "test/unit/phys/dca_step/cluster_solver/stub_rng.hpp"
+
+namespace {
+
+using Scalar = double;
+using Concurrency = dca::parallel::NoConcurrency;
+using Lattice = dca::phys::models::square_lattice<dca::phys::domains::D4>;
+using Model = dca::phys::models::TightBindingModel<Lattice>;
+using Parameters =
+    dca::phys::params::Parameters<Concurrency, dca::parallel::NoThreading, dca::profiling::NullProfiler,
+                                  Model, dca::testing::StubRng, dca::ClusterSolverId::CT_AUX,
+                                  dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>;
+using Data = dca::phys::DcaData<Parameters>;
+
+using NuDmn = Data::NuDmn;
+using RDmn = Data::RClusterDmn;  // 2x4 cluster, Nc=8
+using KDmn = Data::KClusterDmn;
+using WDmn = Data::WDmn;
+using Complex = std::complex<double>;
+
+constexpr char input_file[] =
+    DCA_SOURCE_DIR "/test/unit/phys/dca_loop/disorder/input_disorder_loop.json";
+
+std::unique_ptr<Data> makeData() {
+  static Concurrency concurrency(0, nullptr);
+  static Parameters parameters("", concurrency);
+  static bool initialized = false;
+  if (!initialized) {
+    parameters.read_input_and_broadcast<dca::io::JSONReader>(input_file);
+    parameters.update_model();
+    parameters.update_domains();
+    initialized = true;
+  }
+  return std::make_unique<Data>(parameters);
+}
+
+// Clean, cluster-excluded G0 in frequency: distinct, diagonal-dominant (invertible) nu-blocks per k.
+void fillCleanG0(Data& data) {
+  auto& g0 = data.G0_k_w_cluster_excluded;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int k = 0; k < KDmn::dmn_size(); ++k)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1)
+          g0(n1, n2, k, w) = (n1 == n2) ? Complex(3.0 + 0.1 * k + 0.01 * n1, 0.2 + 0.02 * w)
+                                        : Complex(0.3, n1 < n2 ? 0.1 : -0.1);
+}
+
+// makeDisorderedG0's legacy time-domain loop inverts nu-blocks of G0_r_t_cluster_excluded; give it
+// invertible (identity) blocks. Its output is unused by the DISORDERED_G0 path.
+void fillG0rtIdentity(Data& data) {
+  auto& g = data.G0_r_t_cluster_excluded;
+  g = Scalar(0);
+  for (int t = 0; t < Data::TDmn::dmn_size(); ++t)
+    for (int r = 0; r < RDmn::dmn_size(); ++r)
+      for (int n = 0; n < NuDmn::dmn_size(); ++n)
+        g(n, n, r, t) = 1.0;
+}
+
+// Runs the real disorder-averaging pipeline (loop body of workTheClusters + averageGkw) with M = 0.
+// beta is irrelevant here because the G0*M*G0 term vanishes for M = 0.
+void runDisorderPipeline(Data& data, const std::vector<Data::DisorderConfiguration>& configs,
+                         const std::vector<double>& weights) {
+  Data::Disordered_G0_r_r_w m_zero("m_zero");
+  m_zero = 0;
+  data.disorder_G_r_r_w = 0;
+  for (std::size_t c = 0; c < configs.size(); ++c) {
+    data.makeDisorderedG0(configs[c]);
+    dca::phys::disorder::accumulateDisorderDyson<NuDmn, RDmn, WDmn, Scalar>(
+        data.disordered_G0_r_r_w_cl_exl, m_zero, /*beta=*/1.0, weights[c], data.disorder_G_r_r_w);
+  }
+  const double total = std::accumulate(weights.begin(), weights.end(), 0.0);
+  data.disorder_G_r_r_w /= total;
+  dca::phys::disorder::translationalAverage<NuDmn, RDmn, WDmn>(data.disorder_G_r_r_w, data.G_r_w);
+  dca::math::transform::FunctionTransform<RDmn, KDmn>::execute(data.G_r_w, data.G_k_w);
+}
+
+}  // namespace
+
+// Clean limit: zero disorder potential => G0_dis = clean G0, and the full disorder machinery
+// (unfold -> dense identity Dyson -> disorder average -> translational average -> FT r->k) must
+// reproduce the clean cluster-excluded G0_k_w.
+TEST(DisorderLoopTest, CleanLimitReproducesClusterExcludedG0) {
+  auto data = makeData();
+  fillCleanG0(*data);
+  fillG0rtIdentity(*data);
+
+  // Three identical zero-potential configurations, equal weights.
+  std::vector<Data::DisorderConfiguration> configs(3);  // default-constructed = zero potential
+  std::vector<double> weights(3, 1.0 / 3.0);
+
+  runDisorderPipeline(*data, configs, weights);
+
+  constexpr double tol = 1e-9;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    for (int k = 0; k < KDmn::dmn_size(); ++k)
+      for (int n2 = 0; n2 < NuDmn::dmn_size(); ++n2)
+        for (int n1 = 0; n1 < NuDmn::dmn_size(); ++n1) {
+          const Complex got = data->G_k_w(n1, n2, k, w);
+          const Complex ref = data->G0_k_w_cluster_excluded(n1, n2, k, w);
+          EXPECT_NEAR(got.real(), ref.real(), tol)
+              << "(n1,n2,k,w)=(" << n1 << "," << n2 << "," << k << "," << w << ")";
+          EXPECT_NEAR(got.imag(), ref.imag(), tol);
+        }
+}
+
+// With distinct nonzero configurations and weights, the accumulated/normalized two-site G must equal
+// the explicit weighted mean of the per-config disordered G0 (M = 0 => G_c = G0_dis,c). Also checks
+// the averaged G_k_w is finite and actually shifted away from the clean G0 by the disorder.
+TEST(DisorderLoopTest, DisorderAverageIsWeightedMeanOfConfigs) {
+  auto data = makeData();
+  fillCleanG0(*data);
+  fillG0rtIdentity(*data);
+
+  const int nu = NuDmn::dmn_size();
+  const int nr = RDmn::dmn_size();
+
+  // Three distinct, nonzero, site-local disorder configurations.
+  std::vector<Data::DisorderConfiguration> configs(3);
+  for (int c = 0; c < 3; ++c)
+    for (int r = 0; r < nr; ++r)
+      for (int n = 0; n < nu; ++n)
+        configs[c](n, r) = 0.4 * (c + 1) * (1.0 + 0.1 * r) - 0.2 * n;
+  const std::vector<double> weights = {0.2, 0.5, 0.3};
+  const double total = std::accumulate(weights.begin(), weights.end(), 0.0);
+
+  // Independent expected weighted mean of the per-config disordered G0 (M = 0 => G_c = G0_dis,c).
+  Data::Disordered_G0_r_r_w expected("expected");
+  expected = 0;
+  for (int c = 0; c < 3; ++c) {
+    data->makeDisorderedG0(configs[c]);
+    for (int idx = 0; idx < expected.size(); ++idx)
+      expected(idx) += weights[c] * data->disordered_G0_r_r_w_cl_exl(idx);
+  }
+  for (int idx = 0; idx < expected.size(); ++idx)
+    expected(idx) /= total;
+
+  runDisorderPipeline(*data, configs, weights);
+
+  // The pipeline never modifies G0_k_w_cluster_excluded, so it still holds the clean G0.
+  constexpr double tol = 1e-9;
+  for (int idx = 0; idx < expected.size(); ++idx) {
+    EXPECT_NEAR(data->disorder_G_r_r_w(idx).real(), expected(idx).real(), tol) << "idx=" << idx;
+    EXPECT_NEAR(data->disorder_G_r_r_w(idx).imag(), expected(idx).imag(), tol) << "idx=" << idx;
+  }
+
+  // The averaged G_k_w is finite, and the disorder shifts it away from the clean G0.
+  bool differs = false;
+  for (int idx = 0; idx < data->G_k_w.size(); ++idx) {
+    ASSERT_TRUE(std::isfinite(data->G_k_w(idx).real()) && std::isfinite(data->G_k_w(idx).imag()));
+    if (std::abs(data->G_k_w(idx) - data->G0_k_w_cluster_excluded(idx)) > 1e-6)
+      differs = true;
+  }
+  EXPECT_TRUE(differs) << "nonzero disorder must shift the averaged G_k_w away from the clean G0";
+}
diff --git a/test/unit/phys/dca_loop/disorder/input_disorder_loop.json b/test/unit/phys/dca_loop/disorder/input_disorder_loop.json
new file mode 100644
index 000000000..42492e369
--- /dev/null
+++ b/test/unit/phys/dca_loop/disorder/input_disorder_loop.json
@@ -0,0 +1,28 @@
+{
+  "physics": {
+    "beta": 2,
+    "chemical-potential": 0
+  },
+  "single-band-Hubbard-model": {
+    "t": 1,
+    "U": 0
+  },
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[2, 0],
+        [0, 4]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 64
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 16,
+      "four-point-fermionic-frequencies": 8
+    }
+  },
+  "four-point": {
+    "type": "NONE",
+    "momentum-transfer": [0., 0.],
+    "frequency-transfer": 0
+  }
+}
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..27b70f35a 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,33 @@ 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}
   )
 
+# Test of the two-site (DISORDERED_G0) sp accumulator: correct two-site indexing and off-diagonal
+# (r1 != r2) population of the M function. Only built with the disorder option, never for a GPU build.
+if(DCA_WITH_DISORDER AND NOT DCA_WITH_CUDA AND NOT DCA_WITH_HIP)
+  dca_add_gtest(sp_accumulator_disorder_test
+    FAST
+    GTEST_MAIN
+    INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
+    LIBS     FFTW::Double quantum_domains ${DCA_LIBS}
+    )
+endif()
+
 #;mc_kernels
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator_disorder_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator_disorder_test.cpp
new file mode 100644
index 000000000..505f71338
--- /dev/null
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator_disorder_test.cpp
@@ -0,0 +1,170 @@
+// Copyright (C) 2018 ETH Zurich
+// Copyright (C) 2018 UT-Battelle, LLC
+// All rights reserved.
+//
+// See LICENSE.txt for terms of usage.
+// See CITATION.txt for citation guidelines if you use this code for scientific publications.
+//
+// Author: Anirudha Mirmira (anirudha.mirmira@gmail.com)
+//
+// CPU test for the disordered (DISORDERED_G0) single-particle accumulator. It verifies that the
+// two-site M function M(b1,s,r1,b2,s,r2,w) is filled with genuine two-site structure:
+//   * the value M_ij from a configuration pair (i,j) lands in bin (r_i, r_j) -- ALL entries,
+//     including the off-diagonal r1 != r2 ones, not only the diagonal;
+//   * bins with no contributing pair stay exactly zero (correct, leak-free indexing);
+//   * the two spin sectors are written into their own s slot.
+//
+// The test only exercises real code under DISORDERED_G0 (the macro that switches the accumulator's
+// PDmn to <B,R,B,R> and MFunction to <Nu,R,Nu,R,W>). When the macro is off it degenerates to a
+// trivial passing test so the target still builds in a non-disorder configuration.
+
+#include <array>
+#include <cmath>
+#include <complex>
+#include <vector>
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+
+using Scalar = double;
+
+#include "test/mock_mcconfig.hpp"
+namespace dca {
+namespace config {
+using McOptions = MockMcOptions<Scalar>;
+}  // namespace config
+}  // namespace dca
+
+#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp"
+#include "test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/accumulation_test.hpp"
+
+constexpr int n_bands = 1;
+constexpr int n_sites = 4;
+constexpr int n_frqs = 32;
+
+using SpAccumulatorDisorderTest =
+    dca::testing::AccumulationTest<Scalar, n_bands, n_sites, n_frqs, false>;
+
+using Parameters = SpAccumulatorDisorderTest::Parameters;
+using Configuration = SpAccumulatorDisorderTest::Configuration;  // std::array<vector<Vertex>, 2>
+using MatrixPair = SpAccumulatorDisorderTest::Sample;            // std::array<Matrix, 2>
+using Vertex = dca::testing::Vertex;
+
+using HostAccumulator = dca::phys::solver::accumulator::SpAccumulator<Parameters, dca::linalg::CPU>;
+using MFunction = typename HostAccumulator::MFunction;
+using RDmn = typename HostAccumulator::RDmn;
+using BDmn = typename HostAccumulator::BDmn;
+using WDmn = typename HostAccumulator::WDmn;
+
+namespace {
+
+// Largest |M(b=0, s, r1, b=0, s, r2, w)| over all Matsubara frequencies for a given bin/spin.
+double maxAbsOverW(const MFunction& f, const int s, const int r1, const int r2) {
+  double m = 0.;
+  for (int w = 0; w < WDmn::dmn_size(); ++w)
+    m = std::max(m, std::abs(f(0, s, r1, 0, s, r2, w)));
+  return m;
+}
+
+}  // namespace
+
+// All sites carry band 0 (n_bands == 1); each spin's vertices sit on distinct sites with distinct
+// imaginary times, so every config pair (i,j) maps to its own (r_i, r_j) bin with no cancellation.
+// We then assert exactly which bins are populated.
+TEST_F(SpAccumulatorDisorderTest, OffDiagonalIndexPlacement) {
+  // Spin up lives on sites {0,1,2}; spin down on {2,3}. Site 3 is never touched by spin up, site 0
+  // and 1 never by spin down -- this exposes both off-diagonal placement and leak-free indexing.
+  const std::array<std::vector<int>, 2> sites{std::vector<int>{0, 1, 2}, std::vector<int>{2, 3}};
+
+  Configuration config;
+  MatrixPair M;
+  for (int s = 0; s < 2; ++s) {
+    const int n = sites[s].size();
+    config[s].resize(n);
+    M[s].resize(n);
+    for (int k = 0; k < n; ++k)
+      config[s][k] = Vertex{0, sites[s][k], 0.1 + 0.3 * k + 0.05 * s};  // distinct, in [0, beta]
+    for (int j = 0; j < n; ++j)
+      for (int i = 0; i < n; ++i)
+        M[s](i, j) = 0.5 + i + 2. * j + 0.1 * s;  // never zero
+  }
+
+  HostAccumulator acc(parameters_);
+  acc.resetAccumulation();
+  acc.accumulate(M, config, 1.);  // factor / sign = 1
+  acc.finalize();
+  const MFunction& Mrw = acc.get_sign_times_M_r_w();
+
+  constexpr double zero_tol = 1e-14;
+  int off_diag_nonzero = 0;
+
+  for (int s = 0; s < 2; ++s) {
+    const auto& used = sites[s];
+    auto is_used = [&](int r) { return std::find(used.begin(), used.end(), r) != used.end(); };
+
+    for (int r1 = 0; r1 < RDmn::dmn_size(); ++r1)
+      for (int r2 = 0; r2 < RDmn::dmn_size(); ++r2) {
+        const double mag = maxAbsOverW(Mrw, s, r1, r2);
+        if (is_used(r1) && is_used(r2)) {
+          EXPECT_GT(mag, zero_tol) << "expected nonzero bin s=" << s << " (" << r1 << "," << r2 << ")";
+          if (r1 != r2)
+            ++off_diag_nonzero;
+        }
+        else {
+          EXPECT_EQ(0.0, mag) << "expected exactly-zero bin s=" << s << " (" << r1 << "," << r2 << ")";
+        }
+      }
+  }
+
+  // Spin up contributes 3x3 - 3 = 6 off-diagonal bins, spin down 2x2 - 2 = 2 -> 8 total.
+  EXPECT_EQ(8, off_diag_nonzero) << "off-diagonal (r1 != r2) entries must be populated, not skipped";
+}
+
+// Exact, NFFT-independent check that the VALUE M(i,j) is placed in bin (r_i, r_j) for the full
+// matrix (every off-diagonal included). With all imaginary times equal, every bin shares one
+// frequency profile g(w), so M_rw(bin_ij, w) * M(0,0) == M_rw(bin_00, w) * M(i,j) for all w; g(w)
+// (whatever the NFFT produces) cancels, leaving an exact proportionality of the placed values.
+TEST_F(SpAccumulatorDisorderTest, ValuePlacementFullMatrix) {
+  const std::vector<int> sites{0, 1, 2};  // distinct -> one pair per bin
+  const int n = sites.size();
+
+  Configuration config;
+  MatrixPair M;
+  for (int s = 0; s < 2; ++s) {
+    config[s].resize(n);
+    M[s].resize(n);
+    for (int k = 0; k < n; ++k)
+      config[s][k] = Vertex{0, sites[k], 0.0};  // identical tau -> shared frequency profile
+    for (int j = 0; j < n; ++j)
+      for (int i = 0; i < n; ++i)
+        M[s](i, j) = 1.0 + i + 3. * j;  // real, nonzero, distinct per (i,j)
+  }
+
+  HostAccumulator acc(parameters_);
+  acc.resetAccumulation();
+  acc.accumulate(M, config, 1.);
+  acc.finalize();
+  const MFunction& Mrw = acc.get_sign_times_M_r_w();
+
+  const int s = 0;
+  const int r0 = sites[0];
+  const double w00 = M[s](0, 0);  // reference weight, nonzero
+
+  // bin_00 must be nonzero so the proportionality reference is meaningful.
+  ASSERT_GT(maxAbsOverW(Mrw, s, r0, r0), 1e-14);
+
+  constexpr double tol = 1e-10;
+  for (int j = 0; j < n; ++j)
+    for (int i = 0; i < n; ++i) {
+      const int ri = sites[i];  // right vertex i -> first R index
+      const int rj = sites[j];  // left  vertex j -> second R index
+      const double wij = M[s](i, j);
+      for (int w = 0; w < WDmn::dmn_size(); ++w) {
+        const std::complex<double> lhs = Mrw(0, s, ri, 0, s, rj, w) * w00;
+        const std::complex<double> rhs = Mrw(0, s, r0, 0, s, r0, w) * wij;
+        EXPECT_NEAR(lhs.real(), rhs.real(), tol)
+            << "value misplaced for pair (" << i << "," << j << ") bin (" << ri << "," << rj
+            << ") at w=" << w;
+        EXPECT_NEAR(lhs.imag(), rhs.imag(), tol);
+      }
+    }
+}
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 1acf66cc9..e06ae3d55 100644
--- a/test/unit/phys/parameters/CMakeLists.txt
+++ b/test/unit/phys/parameters/CMakeLists.txt
@@ -1,5 +1,8 @@
 add_subdirectory(analysis_parameters)
 add_subdirectory(dca_parameters)
+if(DCA_WITH_DISORDER AND NOT DCA_WITH_CUDA AND NOT DCA_WITH_HIP)
+  add_subdirectory(disorder_parameters)
+endif()
 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..112db1c33
--- /dev/null
+++ b/test/unit/phys/parameters/disorder_parameters/disorder_parameters_test.cpp
@@ -0,0 +1,61 @@
+// 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("binary", pars.get_disorder_distribution());
+  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());
+  EXPECT_FALSE(pars.get_disorder_unique_configs());
+  EXPECT_FALSE(pars.get_disorder_present());
+}
+
+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("box", pars.get_disorder_distribution());
+  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());
+  EXPECT_TRUE(pars.get_disorder_unique_configs());
+  // A "disorder" section was present in the input.
+  EXPECT_TRUE(pars.get_disorder_present());
+}
+
+// A real input that has no "disorder" section: the parameters keep their (clean) defaults and
+// get_disorder_present() stays false, which is how the disorder flow tells a deliberate clean run
+// apart from a section that is present but misconfigured.
+TEST(DisorderParametersTest, ReadNoDisorderSection) {
+  dca::io::JSONReader reader;
+  dca::phys::params::DisorderParameters pars;
+
+  reader.open_file(DCA_SOURCE_DIR
+                   "/test/unit/phys/parameters/disorder_parameters/input_no_disorder.json");
+  pars.readWrite(reader);
+  reader.close_file();
+
+  EXPECT_FALSE(pars.get_disorder_present());
+  EXPECT_EQ(0, pars.get_disorder_num_configurations());
+}
diff --git a/test/unit/phys/parameters/disorder_parameters/input_no_disorder.json b/test/unit/phys/parameters/disorder_parameters/input_no_disorder.json
new file mode 100644
index 000000000..ecf7b2766
--- /dev/null
+++ b/test/unit/phys/parameters/disorder_parameters/input_no_disorder.json
@@ -0,0 +1,40 @@
+{
+  "physics": {
+    "beta": 2,
+    "chemical-potential": 0
+  },
+  "single-band-Hubbard-model": {
+    "t": 1,
+    "U": 2
+  },
+  "domains": {
+    "real-space-grids": {
+      "cluster": [[4, 0],
+        [0, 1]]
+    },
+    "imaginary-time": {
+      "sp-time-intervals": 64
+    },
+    "imaginary-frequency": {
+      "sp-fermionic-frequencies": 16,
+      "four-point-fermionic-frequencies": 8
+    }
+  },
+  "four-point": {
+    "type": "NONE",
+    "momentum-transfer": [0., 0.],
+    "frequency-transfer": 0
+  },
+  "Monte-Carlo-integration": {
+    "seed": 42,
+    "warm-up-sweeps": 20,
+    "sweeps-per-measurement": 1,
+    "measurements": 150
+  },
+  "CT-AUX": {
+    "expansion-parameter-K": 1.0,
+    "initial-configuration-size": 10,
+    "initial-matrix-size": 64,
+    "max-submatrix-size": 1
+  }
+}
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..bb64c665b
--- /dev/null
+++ b/test/unit/phys/parameters/disorder_parameters/input_read_all.json
@@ -0,0 +1,9 @@
+{
+  "disorder": {
+    "distribution": "box",
+    "potential": 1.0,
+    "density": 0.25,
+    "num-configurations": 10,
+    "unique-configs": true
+  }
+}
diff --git a/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp b/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
index 7a9ac7cfb..e99ea101a 100644
--- a/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
+++ b/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
@@ -38,6 +38,7 @@ TEST(OutputParametersTest, DefaultValues) {
   EXPECT_FALSE(pars.dump_cluster_Greens_functions());
   EXPECT_FALSE(pars.dump_Gamma_lattice());
   EXPECT_FALSE(pars.dump_chi_0_lattice());
+  EXPECT_FALSE(pars.dump_disorder_configs());
   EXPECT_FALSE(pars.dump_every_iteration());
 
   // // And nothing except required input should get updated by reading the minimal input