diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
index 12698ff..77a0f1f 100644
--- a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
@@ -45,6 +45,14 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b2_a2b0:
 				case Label::ab_ab::a1b2_a2b1:
 					return !this->symmetry.is_Aa_in_irreducible_sector(A1);
+			// GPU visits Aa01 before the output Aa2 is known in these cases.
+				case Label::ab_ab::a0b0_a1b1:
+				case Label::ab_ab::a0b1_a1b0:
+				case Label::ab_ab::a0b0_a1b2:
+				case Label::ab_ab::a0b1_a1b2:
+				case Label::ab_ab::a0b2_a1b0:
+				case Label::ab_ab::a0b2_a1b1:
+					return false;
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
 		}
@@ -62,6 +70,11 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a0b1_a1b2:
 				case Label::ab_ab::a0b2_a1b0:
 				case Label::ab_ab::a0b2_a1b1:
+			// a01b2_a2b01: Ab01 is the output right atom on GPU.
+				case Label::ab_ab::a0b2_a2b0:
+				case Label::ab_ab::a0b2_a2b1:
+				case Label::ab_ab::a1b2_a2b0:
+				case Label::ab_ab::a1b2_a2b1:
 					return !this->symmetry.is_Ab_in_irreducible_sector(A1.first);
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
@@ -100,6 +113,14 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b0_a2b2:
 				case Label::ab_ab::a1b1_a2b2:
 					return !this->symmetry.in_irreducible_sector(A1, A3);
+			// GPU output is (Aa2, Ab2) or (Aa2, Ab01); Aa01 is intermediate.
+				case Label::ab_ab::a0b0_a1b1:
+				case Label::ab_ab::a0b1_a1b0:
+				case Label::ab_ab::a0b0_a1b2:
+				case Label::ab_ab::a0b1_a1b2:
+				case Label::ab_ab::a0b2_a1b0:
+				case Label::ab_ab::a0b2_a1b1:
+					return !this->symmetry.in_irreducible_sector(A2, A3);
 			// a01b2_a2b01:                                        Aa01,          Ab01,          Aa2
 			default:
 					return false;
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index c05794f..d10bee9 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -1554,7 +1554,8 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
+									// The output block is (Aa01, Ab01); Ab2 is only an intermediate atom.
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
 									const Tensor<Tdata> &D_b = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1693,7 +1694,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1832,7 +1833,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
 									const Tensor<Tdata> &D_b = Global_Func::find(Ds_b_transpose, Ab01.first, TAC{Ab2.first, (Ab2.second-Ab01.second)%this->period});
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -1971,7 +1972,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
@@ -2106,10 +2107,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
+									// Screen the output pair; Aa2 is summed over, not an output atom.
+									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
@@ -2270,10 +2272,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
+									// Screen the output pair; Aa2 is summed over, not an output atom.
+									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
@@ -2434,10 +2437,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
+									// Screen the output pair; Aa2 is summed over, not an output atom.
+									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
@@ -2598,10 +2602,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								const TAC &Aa2 = list_Aa2[ia2];
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
+									// Screen the output pair; Aa2 is summed over, not an output atom.
+									if(this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 									if(this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
diff --git a/include/RI/ri/LRI_Cal_Aux.h b/include/RI/ri/LRI_Cal_Aux.h
index 1eae641..c5985c8 100644
--- a/include/RI/ri/LRI_Cal_Aux.h
+++ b/include/RI/ri/LRI_Cal_Aux.h
@@ -95,7 +95,7 @@ namespace LRI_Cal_Aux
 		if (D_result.empty())
 		{
 			if(Tdata(1.0)==fac)
-				D_result = D_add;
+				D_result = D_add.copy();
 			else if(Tdata(-1.0)==fac)
 				D_result = -D_add;
 			else
diff --git a/include/RI/ri/LRI_k-cal_cvc_mo.hpp b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
index d540770..594bd18 100644
--- a/include/RI/ri/LRI_k-cal_cvc_mo.hpp
+++ b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
@@ -467,7 +467,7 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_hartree_onthefly(
 		{
 			for (int k1_index=0; k1_index<k1_indices.size(); ++k1_index)
 			{
-				const TA& mu = list_J[mu_index];
+				const TA& mu = list_I[mu_index];
 				const int k1 = k1_indices[k1_index];
 				// 2.1 calculate C^mu (i,a^*)[k1,k1] on-the-fly
 				const Tensor<Tdata> C_mu_ia = Cs_ao_mo_to_Cs_mo(Cs_ao_mo, map_psi, k1, k1,
diff --git a/include/RI/ri/LRI_k-cal_hartree.hpp b/include/RI/ri/LRI_k-cal_hartree.hpp
index 9705106..7f6063a 100644
--- a/include/RI/ri/LRI_k-cal_hartree.hpp
+++ b/include/RI/ri/LRI_k-cal_hartree.hpp
@@ -121,7 +121,7 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 #pragma omp for schedule(dynamic) collapse(2)
 		for (int v_index=0; v_index<list_I.size(); v_index++)
 		{
-			for (int u_index=0; u_index<list_I.size(); u_index++)
+			for (int u_index=0; u_index<list_J.size(); u_index++)
 			{
 				const TA v = list_I[v_index];
 				const TA u = list_J[u_index];
diff --git a/unittests/gpu_symmetry/CMakeLists.txt b/unittests/gpu_symmetry/CMakeLists.txt
new file mode 100644
index 0000000..f385c36
--- /dev/null
+++ b/unittests/gpu_symmetry/CMakeLists.txt
@@ -0,0 +1,36 @@
+cmake_minimum_required(VERSION 3.18)
+project(libri_gpu_symmetry_tests LANGUAGES CXX)
+
+option(ENABLE_GPU "Run the DDLA/CUDA contraction checks" OFF)
+set(LIBRI_ROOT "${CMAKE_CURRENT_SOURCE_DIR}/../.." CACHE PATH "LibRI source tree")
+set(LIBCOMM_ROOT "" CACHE PATH "LibComm installation")
+set(CEREAL_ROOT "" CACHE PATH "cereal installation")
+set(DDLA_ROOT "" CACHE PATH "DDLA installation")
+
+find_package(MPI REQUIRED COMPONENTS CXX)
+find_package(OpenMP REQUIRED COMPONENTS CXX)
+find_package(BLAS REQUIRED)
+find_package(LAPACK REQUIRED)
+
+add_executable(test_lri_gpu_symmetry test_lri_gpu_symmetry.cpp)
+target_compile_features(test_lri_gpu_symmetry PRIVATE cxx_std_14)
+target_include_directories(test_lri_gpu_symmetry PRIVATE
+  "${LIBRI_ROOT}/include" "${LIBCOMM_ROOT}/include" "${CEREAL_ROOT}/include")
+target_link_libraries(test_lri_gpu_symmetry PRIVATE
+  MPI::MPI_CXX OpenMP::OpenMP_CXX ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
+
+if(ENABLE_GPU)
+  find_package(CUDAToolkit REQUIRED)
+  find_library(DDLA_LIBRARY NAMES ddla PATHS "${DDLA_ROOT}/lib" REQUIRED NO_DEFAULT_PATH)
+  target_compile_features(test_lri_gpu_symmetry PRIVATE cxx_std_17)
+  target_compile_definitions(test_lri_gpu_symmetry PRIVATE __GPU_RI __DDLA_RI DDLA_USE_CUDA)
+  target_include_directories(test_lri_gpu_symmetry PRIVATE "${DDLA_ROOT}/include")
+  target_link_libraries(test_lri_gpu_symmetry PRIVATE
+    ${DDLA_LIBRARY} CUDA::cudart CUDA::cublas CUDA::curand CUDA::cusolver)
+endif()
+
+enable_testing()
+add_test(NAME lri_symmetry COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 1
+  ${MPIEXEC_PREFLAGS} $<TARGET_FILE:test_lri_gpu_symmetry> ${MPIEXEC_POSTFLAGS})
+set_tests_properties(lri_symmetry PROPERTIES ENVIRONMENT
+  "OMP_NUM_THREADS=1;OPENBLAS_NUM_THREADS=1;MKL_NUM_THREADS=1")
diff --git a/unittests/gpu_symmetry/README.md b/unittests/gpu_symmetry/README.md
new file mode 100644
index 0000000..167f271
--- /dev/null
+++ b/unittests/gpu_symmetry/README.md
@@ -0,0 +1,35 @@
+# GPU symmetry regression
+
+Compare all 18 contractions, individually and together, against full CPU
+contractions with an independently applied output-sector mask. The fixture uses
+two unequal atoms and three periodic cells, real and complex tensors, both CPU
+modes, and optionally DDLA/CUDA. It checks unrestricted, restricted and empty
+output, as well as seeded signed, zero and fractional accumulation.
+
+```sh
+cmake -S unittests/gpu_symmetry -B build-symmetry \
+  -DLIBCOMM_ROOT=/path/to/LibComm -DCEREAL_ROOT=/path/to/cereal \
+  -DENABLE_GPU=ON -DDDLA_ROOT=/path/to/DDLA
+cmake --build build-symmetry --parallel 2
+ctest --test-dir build-symmetry --output-on-failure
+```
+
+The GPU test requires one GPU, one MPI rank, and DDLA built for CUDA. Omit
+`ENABLE_GPU` and `DDLA_ROOT` to check a CPU-only build. The library headers can
+be selected independently with `LIBRI_ROOT`, allowing the same test to reproduce
+failures in an earlier revision.
+
+For OpenMP coverage or an isolated failing contraction:
+
+```sh
+OMP_NUM_THREADS=4 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 \
+  mpirun -np 1 build-symmetry/test_lri_gpu_symmetry
+OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 \
+  mpirun -np 1 build-symmetry/test_lri_gpu_symmetry 1 gpu real
+```
+
+Arguments select the contraction index (`0`–`17`, `18` for the combined sum,
+or `-1` for all), backend (`cpu`, `gpu`, `all`) and scalar (`real`, `complex`,
+`all`). Index 1 is `a1b2_a2b1`: filtering `(Aa01,Aa2)` while summing over `Aa2`
+incorrectly removes contributions to the accepted output pair `(Aa01,Ab01)`.
+The off-diagonal sector with a nonzero translation detects this error.
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
new file mode 100644
index 0000000..233bc34
--- /dev/null
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -0,0 +1,242 @@
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <map>
+#include <memory>
+#include <set>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include <RI/physics/symmetry/Filter_Atom_Symmetry.h>
+#include <RI/ri/LRI.h>
+
+namespace
+{
+using Cell = std::array<int, 1>;
+using AtomCell = std::pair<int, Cell>;
+using Sector = std::map<std::pair<int, int>, std::set<Cell>>;
+using Label = RI::Label::ab_ab;
+
+template <typename T>
+using Tensors = std::map<int, std::map<AtomCell, RI::Tensor<T>>>;
+
+template <typename T>
+RI::Tensor<T> make_tensor(const RI::Shape_Vector& shape, int left, int right,
+                          int cell, int label, T imaginary)
+{
+    RI::Tensor<T> tensor(shape);
+    for (std::size_t i = 0; i < tensor.data->size(); ++i)
+    {
+        // Different atoms and cells must not accidentally produce equal sums.
+        const double real = 0.2 + 0.07 * left + 0.03 * right + 0.02 * cell
+                            + 0.005 * label + 0.013 * i;
+        const double imag = 0.04 + 0.02 * left - 0.01 * right + 0.003 * i;
+        (*tensor.data)[i] = T(real) + imaginary * T(imag);
+    }
+    return tensor;
+}
+
+template <typename T>
+Tensors<T> copy_tensors(const Tensors<T>& source)
+{
+    Tensors<T> result;
+    for (const auto& left : source)
+        for (const auto& right : left.second)
+            result[left.first][right.first] = right.second.copy();
+    return result;
+}
+
+template <typename T>
+void compare(const Tensors<T>& actual, const Tensors<T>& expected,
+             const std::string& context)
+{
+    if (actual.size() != expected.size())
+        throw std::runtime_error(context + ": output atom count differs");
+    for (const auto& left : expected)
+    {
+        const auto found_left = actual.find(left.first);
+        if (found_left == actual.end() || found_left->second.size() != left.second.size())
+            throw std::runtime_error(context + ": output block count differs");
+        for (const auto& right : left.second)
+        {
+            const auto found = found_left->second.find(right.first);
+            if (found == found_left->second.end())
+                throw std::runtime_error(context + ": missing output block");
+            const auto& value = found->second;
+            const auto& reference = right.second;
+            if (value.shape.size() != reference.shape.size()
+                || !std::equal(value.shape.begin(), value.shape.end(), reference.shape.begin()))
+                throw std::runtime_error(context + ": output shape differs");
+            for (std::size_t i = 0; i < reference.data->size(); ++i)
+            {
+                const double error = std::abs((*value.data)[i] - (*reference.data)[i]);
+                if (!std::isfinite(error) || error > 1e-10 * (1.0 + std::abs((*reference.data)[i])))
+                {
+                    std::cerr << context << " atoms=" << left.first << ',' << right.first.first
+                              << " cell=" << right.first.second[0] << " element=" << i
+                              << " error=" << error << '\n';
+                    throw std::runtime_error(context + ": contraction differs from masked CPU reference");
+                }
+            }
+        }
+    }
+}
+
+template <typename T>
+void check(const std::string& scalar, T imaginary, int selected, const std::string& backend)
+{
+    RI::LRI<int, int, 1, T> lri;
+    const std::vector<Label> all_labels(RI::Label::array_ab_ab.begin(), RI::Label::array_ab_ab.end());
+    const std::map<int, std::array<double, 1>> atoms{{0, {0.0}}, {1, {0.4}}};
+    const std::array<std::array<double, 1>, 1> lattice{{{{1.0}}}};
+    const Cell period{{3}};
+    lri.set_parallel(MPI_COMM_WORLD, atoms, lattice, period, all_labels);
+    int label_index = 0;
+    for (const auto label : RI::Label::array_ab)
+    {
+        const bool coefficient = label == RI::Label::ab::a || label == RI::Label::ab::b;
+        const RI::Shape_Vector shape = coefficient ? RI::Shape_Vector{2, 2, 2} : RI::Shape_Vector{2, 2};
+        Tensors<T> tensors;
+        for (int left = 0; left < 2; ++left)
+            for (int right = 0; right < 2; ++right)
+                for (int cell = 0; cell < period[0]; ++cell)
+                    tensors[left][{right, {cell}}] = make_tensor(shape, left, right, cell, label_index, imaginary);
+        lri.set_tensors_map2(tensors, {label});
+        ++label_index;
+    }
+
+    // Empty means no reduction in the existing LibRI API. The other sectors
+    // exercise unequal atom identities, relative translations, and empty output.
+    const std::vector<Sector> sectors{
+        {},
+        {{{1, 0}, {Cell{{1}}}}},
+        {{{0, 1}, {Cell{{0}}, Cell{{-1}}}}, {{1, 1}, {Cell{{1}}}}},
+        {{{7, 9}, {Cell{{1}}}}}
+    };
+    std::vector<RI::LRI_Cal_Mode> modes;
+    if (backend != "gpu")
+    {
+        modes.push_back(RI::LRI_Cal_Mode::CPU);
+        modes.push_back(RI::LRI_Cal_Mode::CPU_fine_grained_lock);
+    }
+#ifdef __GPU_RI
+    if (backend != "cpu")
+        modes.push_back(RI::LRI_Cal_Mode::GPU);
+#endif
+    if (modes.empty())
+        throw std::runtime_error("Requested backend is not compiled in");
+
+    const int count = static_cast<int>(all_labels.size());
+    for (int index = 0; index <= count; ++index)
+    {
+        if (selected >= 0 && index != selected)
+            continue;
+        std::vector<Label> labels;
+        std::string name;
+        if (index == count)
+        {
+            labels = all_labels;
+            name = "all";
+        }
+        else
+        {
+            labels.push_back(all_labels[index]);
+            name = RI::Label_Tools::get_name(all_labels[index]);
+        }
+        lri.filter_atom = std::make_shared<RI::Filter_Atom<int, AtomCell>>();
+        lri.cal_mode = RI::LRI_Cal_Mode::CPU;
+        Tensors<T> full;
+        lri.cal_loop3(labels, full, 1.0);
+        if (full.size() != 2)
+            throw std::runtime_error("CPU reference must contain both atoms");
+        for (const auto& left : full)
+            if (left.second.size() != 6)
+                throw std::runtime_error("CPU reference must contain all atom/cell blocks");
+
+        for (std::size_t sector_index = 0; sector_index < sectors.size(); ++sector_index)
+        {
+            const auto& sector = sectors[sector_index];
+            Tensors<T> expected;
+            for (const auto& left : full)
+                for (const auto& right : left.second)
+                {
+                    const auto pair = sector.find({left.first, right.first.first});
+                    const bool keep = sector.empty()
+                                      || (pair != sector.end() && pair->second.count(right.first.second));
+                    if (keep)
+                        expected[left.first][right.first] = right.second.copy();
+                }
+            lri.filter_atom = std::make_shared<RI::Filter_Atom_Symmetry<int, Cell, T>>(period, sector);
+            for (const auto mode : modes)
+            {
+                lri.cal_mode = mode;
+                const std::string context = scalar + " label=" + name + " sector="
+                                            + std::to_string(sector_index) + " mode="
+                                            + std::to_string(static_cast<int>(mode));
+                std::cout << "CHECK " << context << std::endl;
+                Tensors<T> actual;
+                lri.cal_loop3(labels, actual, 1.0);
+                compare(actual, expected, context);
+
+                // Seed every block, including those outside the selected sector.
+                // Filtering must preserve them while signed/zero increments work.
+                Tensors<T> seed;
+                for (const auto& left : full)
+                    for (const auto& right : left.second)
+                        seed[left.first][right.first] = make_tensor(right.second.shape, left.first,
+                            right.first.first, right.first.second[0], 0, imaginary);
+                actual = copy_tensors(seed);
+                double sum = 0;
+                for (const double factor : {-0.5, 0.0, 1.25})
+                {
+                    lri.cal_loop3(labels, actual, factor);
+                    sum += factor;
+                    Tensors<T> accumulated = copy_tensors(seed);
+                    for (const auto& left : expected)
+                        for (const auto& right : left.second)
+                        {
+                            auto& target = accumulated.at(left.first).at(right.first);
+                            for (std::size_t i = 0; i < target.data->size(); ++i)
+                                (*target.data)[i] += T(sum) * (*right.second.data)[i];
+                        }
+                    compare(actual, accumulated, context + " factor=" + std::to_string(factor));
+                }
+            }
+        }
+        std::cout << "PASS " << scalar << " label=" << name << std::endl;
+    }
+}
+} // namespace
+
+int main(int argc, char** argv)
+{
+    int provided = 0;
+    if (MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided) != MPI_SUCCESS)
+        return 1;
+    int status = 0;
+    try
+    {
+        if (RI::MPI_Wrapper::mpi_get_size(MPI_COMM_WORLD) != 1 || provided < MPI_THREAD_MULTIPLE)
+            throw std::runtime_error("This test requires one MPI rank and MPI_THREAD_MULTIPLE");
+        const int selected = argc > 1 ? std::stoi(argv[1]) : -1;
+        const std::string backend = argc > 2 ? argv[2] : "all";
+        const std::string scalar = argc > 3 ? argv[3] : "all";
+        if (selected < -1 || selected > 18 || (backend != "cpu" && backend != "gpu" && backend != "all")
+            || (scalar != "real" && scalar != "complex" && scalar != "all"))
+            throw std::runtime_error("Usage: test_lri_gpu_symmetry [-1..18] [cpu|gpu|all] [real|complex|all]");
+        if (scalar != "complex")
+            check<double>("real", 0.0, selected, backend);
+        if (scalar != "real")
+            check<std::complex<double>>("complex", {0.0, 1.0}, selected, backend);
+    }
+    catch (const std::exception& error)
+    {
+        std::cerr << error.what() << '\n';
+        status = 1;
+    }
+    MPI_Finalize();
+    return status;
+}
diff --git a/unittests/physics/test_exx_gpu_symmetry.cpp b/unittests/physics/test_exx_gpu_symmetry.cpp
new file mode 100644
index 0000000..309336a
--- /dev/null
+++ b/unittests/physics/test_exx_gpu_symmetry.cpp
@@ -0,0 +1,99 @@
+// Standalone callback regression: g++ -std=c++11 -Iinclude unittests/physics/test_exx_gpu_symmetry.cpp -o /tmp/test_exx_gpu_symmetry
+#include <cmath>
+#include "RI/physics/symmetry/Filter_Atom_Symmetry.h"
+
+#include <array>
+#include <cassert>
+#include <map>
+#include <set>
+#include <utility>
+
+int main()
+{
+    using Cell = std::array<int, 1>;
+    using AtomCell = std::pair<int, Cell>;
+    using Label = RI::Label::ab_ab;
+    const Cell zero{{0}};
+    const Cell period{{2}};
+    const Cell one{{1}};
+    const std::map<std::pair<int, int>, std::set<Cell>> sector{
+        {{1, 0}, {zero}}
+    };
+    RI::Filter_Atom_Symmetry<int, Cell, double> filter(period, sector);
+    const AtomCell left{1, zero};
+    const AtomCell right{0, zero};
+    const AtomCell intermediate{2, zero};
+    const AtomCell shifted_right{0, one};
+    std::set<Label> checked;
+
+    // GPU output: (Aa2, Ab2); Aa01 is only an intermediate atom.
+    const std::array<Label, 2> output_a2_b2{{Label::a0b0_a1b1, Label::a0b1_a1b0}};
+    for (const Label label : output_a2_b2)
+    {
+        checked.insert(label);
+        assert(!filter.filter_for1(label, 2));
+        assert(!filter.filter_for32(label, 2, left, right));
+        assert(filter.filter_for32(label, 2, right, left));
+        assert(filter.filter_for32(label, 2, left, shifted_right));
+    }
+
+    // GPU output: (Aa2, Ab01); Ab01 can already be screened on its own.
+    const std::array<Label, 4> output_a2_b01{{Label::a0b0_a1b2, Label::a0b1_a1b2,
+                                              Label::a0b2_a1b0, Label::a0b2_a1b1}};
+    for (const Label label : output_a2_b01)
+    {
+        checked.insert(label);
+        assert(!filter.filter_for1(label, 2));
+        assert(!filter.filter_for1(label, right));
+        assert(filter.filter_for1(label, left));
+        assert(!filter.filter_for32(label, 2, left, right));
+        assert(filter.filter_for32(label, 2, right, left));
+        assert(filter.filter_for32(label, 2, left, shifted_right));
+    }
+
+    // GPU output: (Aa01, Ab2); Ab01 is only an intermediate atom.
+    const std::array<Label, 4> output_a01_b2{{Label::a0b0_a2b1, Label::a0b1_a2b0,
+                                              Label::a1b0_a2b1, Label::a1b1_a2b0}};
+    for (const Label label : output_a01_b2)
+    {
+        checked.insert(label);
+        assert(!filter.filter_for1(label, left.first));
+        assert(filter.filter_for1(label, right.first));
+        assert(!filter.filter_for32(label, left.first, intermediate, right));
+        assert(filter.filter_for32(label, left.first, intermediate, shifted_right));
+        assert(filter.filter_for32(label, right.first, intermediate, left));
+    }
+
+    // GPU output: (Aa01, Ab01); Ab2 is only an intermediate atom.
+    const std::array<Label, 4> output_a01_b01_a2b2{{Label::a0b0_a2b2, Label::a0b1_a2b2,
+                                                     Label::a1b0_a2b2, Label::a1b1_a2b2}};
+    for (const Label label : output_a01_b01_a2b2)
+    {
+        checked.insert(label);
+        assert(!filter.filter_for1(label, left.first));
+        assert(filter.filter_for1(label, right.first));
+        assert(!filter.filter_for32(label, left.first, intermediate, right));
+        assert(filter.filter_for32(label, left.first, intermediate, shifted_right));
+        assert(filter.filter_for32(label, right.first, intermediate, left));
+    }
+
+    // GPU output: (Aa01, Ab01); this group screens the pair in filter_for2.
+    const std::array<Label, 4> output_a01_b01_a01b2{{Label::a0b2_a2b0, Label::a0b2_a2b1,
+                                                      Label::a1b2_a2b0, Label::a1b2_a2b1}};
+    for (const Label label : output_a01_b01_a01b2)
+    {
+        checked.insert(label);
+        assert(!filter.filter_for1(label, right));
+        assert(filter.filter_for1(label, left));
+        assert(!filter.filter_for1(label, left.first));
+        assert(!filter.filter_for2(label, left.first, right));
+        assert(filter.filter_for2(label, left.first, shifted_right));
+        assert(filter.filter_for2(label, right.first, left));
+    }
+
+    assert(checked.size() == RI::Label::array_ab_ab.size());
+    for (const Label label : RI::Label::array_ab_ab)
+    {
+        assert(checked.count(label) == 1);
+    }
+}
diff --git a/unittests/ri/LRI-k-atom-lists-test.hpp b/unittests/ri/LRI-k-atom-lists-test.hpp
new file mode 100644
index 0000000..e43ee86
--- /dev/null
+++ b/unittests/ri/LRI-k-atom-lists-test.hpp
@@ -0,0 +1,104 @@
+#pragma once
+
+#include "RI/ri/LRI_k.h"
+
+#include <cmath>
+#include <complex>
+#include <iostream>
+#include <stdexcept>
+
+namespace LRI_K_Atom_Lists_Test
+{
+using Scalar = std::complex<double>;
+using LRI = RI::LRI_k<int, int, 3, Scalar>;
+using Tensors = std::map<int, std::map<LRI::TAC, RI::Tensor<Scalar>>>;
+
+inline RI::Tensor<Scalar> tensor(const RI::Shape_Vector &shape, double value)
+{
+    RI::Tensor<Scalar> result(shape);
+    for (auto &element : *result.data)
+        element = value;
+    return result;
+}
+
+inline void setup(LRI &lri)
+{
+    lri.set_parallel(MPI_COMM_SELF, {{0, {0, 0, 0}}, {1, {0.3, 0, 0}}, {2, {0.7, 0, 0}}},
+                     {{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}}, {1, 1, 1},
+                     RI::Global_Func::to_vector(RI::Label::array_ab_ab));
+    Tensors cs, vs;
+    for (int a = 0; a < 3; ++a)
+        for (int b = 0; b < 3; ++b)
+        {
+            cs[a][{b, {0, 0, 0}}] = tensor({1, 1, 1}, 1 + a + 2*b);
+            vs[a][{b, {0, 0, 0}}] = tensor({1, 1}, 2 + a + 3*b);
+        }
+    lri.set_tensors_map2(cs, {RI::Label::ab::a}, {}, "Cs");
+    lri.set_tensors_map2(vs, {RI::Label::ab::a0b0}, {}, "Vs");
+}
+
+inline void check(const RI::Tensor<Scalar> &actual, double expected)
+{
+    if (actual.empty() || actual.data->size() != 1)
+        throw std::runtime_error("Missing atom-list contraction output");
+    const double error = std::abs((*actual.data)[0] - expected);
+    if (!std::isfinite(error) || error > 1e-11 * (1 + std::abs(expected)))
+        throw std::runtime_error("Atom-list contraction differs from scalar reference");
+}
+
+inline void test_hartree()
+{
+    LRI lri;
+    setup(lri);
+    const std::vector<int> left = {0, 1, 2}, right = {0, 1};
+    std::map<int, std::map<int, std::map<int, RI::Tensor<Scalar>>>> density;
+    double m[3] = {}, n[3] = {};
+    for (int v : left)
+        for (int u : right)
+        {
+            const double d = 0.1 * (1 + 2*v + u);
+            density[v][u][0] = tensor({1, 1}, d);
+            m[u] += (1 + u + 2*v) * d;
+            m[v] += (1 + v + 2*u) * d;
+        }
+    for (int mu : left)
+        for (int nu : left)
+            n[mu] += (2 + mu + 3*nu) * m[nu];
+
+    const auto result = lri.cal_cvcd_k_hartree(
+        density, {{0, 0, 0}}, {0}, left, right, left, "Cs", "Vs");
+    for (int s : left)
+        for (int t : right)
+            check(result.at(s).at(t).at(0),
+                  (1 + s + 2*t) * n[s] + (1 + t + 2*s) * n[t]);
+    std::cout << "PASS Hartree with unequal left/right atom lists" << std::endl;
+}
+
+inline void test_cvc()
+{
+    LRI lri;
+    setup(lri);
+    std::map<int, std::map<int, RI::Tensor<Scalar>>> cs_mo, psi;
+    for (int atom = 0; atom < 3; ++atom)
+    {
+        cs_mo[0][atom] = tensor({1, 1, 2}, atom + 1);
+        psi[0][atom] = tensor({2, 1}, 1);
+    }
+    const std::vector<std::pair<std::vector<int>, std::vector<int>>> partitions = {
+        {{0}, {1}}, {{0}, {1, 2}}, {{0, 1}, {2}}};
+    for (const auto &partition : partitions)
+        for (bool is_A : {false, true})
+        {
+            double expected = 0;
+            for (int mu : partition.first)
+                for (int nu : partition.second)
+                    // Each MO coefficient contains its two conjugate contributions.
+                    expected += 4 * (mu + 1) * (2 + mu + 3*nu) * (nu + 1);
+            const auto result = lri.cal_cvc_mo_k_hartree_onthefly(
+                cs_mo, psi, {0}, {0}, partition.first, partition.second,
+                {"O", "V", "O", "V"}, 1, 1, "Vs", is_A);
+            check(result.at(0).at(0), expected);
+        }
+    std::cout << "PASS CVC with distinct and unequal left/right atom lists" << std::endl;
+}
+}
diff --git a/unittests/ri/LRI-loop3-test.hpp b/unittests/ri/LRI-loop3-test.hpp
index 44dd0ed..0c871e0 100644
--- a/unittests/ri/LRI-loop3-test.hpp
+++ b/unittests/ri/LRI-loop3-test.hpp
@@ -13,6 +13,11 @@
 #include"unittests/print_stl.h"
 #include"unittests/global/Tensor-test.h"
 
+#include<cmath>
+#include<complex>
+#include<limits>
+#include<stdexcept>
+#include<vector>
 #include<map>
 #include<unordered_map>
 #include<iostream>
@@ -396,12 +401,124 @@ namespace LRI_Loop3_Test
 	}
 	
 
+	template<typename Tdata>
+	void test_const_tensor_accumulation()
+	{
+		using Treal = RI::Global_Func::To_Real_t<Tdata>;
+		RI::Tensor<Tdata> input({2});
+		for(std::size_t i=0; i<2; ++i)
+			(*input.data)[i] = RI::Global_Func::convert<Tdata>(std::complex<Treal>(i+1, 0.5));
+		const RI::Tensor<Tdata> &source = input;
+		RI::Tensor<Tdata> first, second;
+		RI::LRI_Cal_Aux::add_Ds(source, first, Tdata(1));
+		RI::LRI_Cal_Aux::add_Ds(source, second, Tdata(1));
+		RI::LRI_Cal_Aux::add_Ds(source, first, Tdata(2));
+		for(std::size_t i=0; i<2; ++i)
+		{
+			const Tdata expected = RI::Global_Func::convert<Tdata>(std::complex<Treal>(i+1, 0.5));
+			if((*source.data)[i]!=expected || (*second.data)[i]!=expected
+				|| (*first.data)[i]!=Tdata(3)*expected)
+				throw std::runtime_error("Accumulating a const tensor changed its input or another result");
+		}
+		std::cout<<"PASS const tensor accumulation preserves input and independent results"<<std::endl;
+	}
+
+
+	template<typename Tdata>
+	void test_loop3_accumulation(
+		const RI::LRI_Cal_Mode cal_mode,
+		const std::vector<Tdata> &factors = {Tdata(0), Tdata(-1), Tdata(1), Tdata(1),
+			Tdata(-1), Tdata(0.5), Tdata(0), Tdata(-0.5), Tdata(0.3)})
+	{
+		using LRI = RI::LRI<int,int,1,Tdata>;
+		using Treal = typename LRI::Tdata_real;
+		using T_Ds = std::map<int, std::map<typename LRI::TAC, RI::Tensor<Tdata>>>;
+		const typename LRI::TAC key = {0,{0}};
+		const auto labels = RI::Global_Func::to_vector(RI::Label::array_ab_ab);
+		const auto make_tensor = [](const RI::Shape_Vector &shape)
+		{
+			RI::Tensor<Tdata> D(shape);
+			for(std::size_t i=0; i<D.data->size(); ++i)
+				(*D.data)[i] = RI::Global_Func::convert<Tdata>(
+					std::complex<Treal>(0.5+0.07*i, 0.1+0.03*i));
+			return D;
+		};
+
+		LRI lri;
+		lri.set_parallel(MPI_COMM_SELF, {{0,{0.0}}}, {{{1.0}}}, {1}, labels);
+		for(const RI::Label::ab label : RI::Label::array_ab)
+		{
+			const bool coefficient = label==RI::Label::ab::a || label==RI::Label::ab::b;
+			T_Ds Ds;
+			Ds[0][key] = make_tensor(coefficient ? RI::Shape_Vector{2,2,2} : RI::Shape_Vector{2,2});
+			lri.set_tensors_map2(Ds, {label});
+		}
+
+		const RI::Tensor<Tdata> seed = make_tensor({2,2});
+		const Treal tolerance = 128*std::numeric_limits<Treal>::epsilon();
+		std::vector<std::vector<RI::Label::ab_ab>> label_groups;
+		for(const RI::Label::ab_ab label : labels)
+			label_groups.push_back({label});
+		label_groups.push_back(labels);
+
+		for(const auto &cal_labels : label_groups)
+		{
+			T_Ds reference;
+			lri.cal_mode = RI::LRI_Cal_Mode::CPU;
+			lri.cal_loop3(cal_labels, reference);
+			const RI::Tensor<Tdata> contribution = reference.at(0).at(key).copy();
+			const Treal norm = contribution.norm(2);
+			if(!std::isfinite(norm) || norm==0)
+				throw std::runtime_error("Invalid CPU reference contraction");
+
+			lri.cal_mode = cal_mode;
+			for(const bool seeded : {false, true})
+			{
+				T_Ds result;
+				// Tensor copies share storage unless copy() is used explicitly.
+				if(seeded)
+					result[0][key] = seed.copy();
+				Tdata accumulated_factor = 0;
+				Treal accumulated_magnitude = 0;
+				for(const Tdata factor : factors)
+				{
+					lri.cal_loop3(cal_labels, result, factor);
+					accumulated_factor += factor;
+					accumulated_magnitude += std::abs(factor);
+					const RI::Tensor<Tdata> &actual = result.at(0).at(key);
+					if(result.size()!=1 || result.at(0).size()!=1 || actual.shape.size()!=2
+						|| actual.shape[0]!=2 || actual.shape[1]!=2)
+						throw std::runtime_error("Unexpected contraction output shape");
+					for(std::size_t i=0; i<seed.data->size(); ++i)
+					{
+						const Tdata initial = seeded ? (*seed.data)[i] : Tdata(0);
+						const Tdata expected = initial + accumulated_factor*(*contribution.data)[i];
+						const Treal error = std::abs((*actual.data)[i] - expected);
+						const Treal scale = 1 + std::abs(initial)
+							+ accumulated_magnitude*std::abs((*contribution.data)[i]);
+						if(!std::isfinite(error) || error>tolerance*scale)
+						{
+							std::cerr<<RI::Label_Tools::get_name(cal_labels)
+								<<" mode="<<static_cast<int>(cal_mode)<<" seeded="<<seeded
+								<<" factor="<<factor<<" element="<<i<<" error="<<error<<std::endl;
+							throw std::runtime_error("Accumulation differs from the CPU reference");
+						}
+					}
+				}
+			}
+		}
+		std::cout<<"Accumulation: all 18 contractions and their sum passed, mode="<<static_cast<int>(cal_mode)<<std::endl;
+	}
+
+
 	template<typename Tdata>
 	void main(int argc, char *argv[])
 	{
 		int mpi_init_provide;
 		MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &mpi_init_provide);
 
+		test_const_tensor_accumulation<Tdata>();
+
 		const std::vector<RI::LRI_Cal_Mode> cal_modes = {
 			RI::LRI_Cal_Mode::CPU,
 			RI::LRI_Cal_Mode::CPU_fine_grained_lock,
@@ -410,7 +527,16 @@ namespace LRI_Loop3_Test
 			#endif
 		};
 		for(const RI::LRI_Cal_Mode cal_mode: cal_modes)
+		{
 			test_loop3<Tdata>(cal_mode);
+			test_loop3_accumulation<Tdata>(cal_mode);
+			// Exercise imaginary coefficients, including the first write to an empty result.
+			using Treal = RI::Global_Func::To_Real_t<Tdata>;
+			test_loop3_accumulation<Tdata>(cal_mode, {
+				RI::Global_Func::convert<Tdata>(std::complex<Treal>(0.25, 0.5)),
+				RI::Global_Func::convert<Tdata>(std::complex<Treal>(0, -1)),
+				RI::Global_Func::convert<Tdata>(std::complex<Treal>(-0.25, 0.5))});
+		}
 
 		MPI_Finalize();
 	}