From bc79bb2ca52af2dbd9ab612c38e248f6c9dd96f1 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Fri, 23 Oct 2026 09:54:06 +0800
Subject: [PATCH 1/9] Fix scaled accumulation in GPU contractions

Apply fac_add_Ds while merging downloaded GPU contributions through the
existing tensor accumulation helper. Remove the unused GPU result map
and locks; the result merge runs outside the OpenMP packing regions.

Add regression coverage for signed, zero and fractional coefficients,
seeded outputs, repeated calls and combined contraction labels.
---
 include/RI/ri/LRI-cal_loop3_GPU.hpp | 45 ++++++--------
 include/RI/ri/gpu/GPU_Data_Output.h |  9 ++-
 unittests/ri/LRI-loop3-test.hpp     | 93 ++++++++++++++++++++++++++++-
 3 files changed, 115 insertions(+), 32 deletions(-)

diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index d2b220f..30d7323 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -53,11 +53,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	GPU_Backend::Context gpu_context(this->mpi_comm);
 	GPU_Backend::Queue &queue = gpu_context.queue();
 
-	std::map<TA, omp_lock_t> lock_Ds_result_add_map = LRI_Cal_Aux::init_lock_result(labels, this->parallel->list_A, Ds_result);
-
 	{
-		std::map<TA, std::map<TAC, Tensor<Tdata>>> Ds_result_thread;
-
 		for(const Label::ab_ab &label : labels)
 		{
 			const std::vector<TA>  list_Aa01_Da = LRI_Cal_Aux::filter_list_map( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].a01, data_wrapper(Label::ab::a).Ds_ab );
@@ -202,7 +198,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b0_a1b1
 
 				case Label::ab_ab::a0b1_a1b0:
@@ -338,7 +334,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b1_a1b0
 
 			  // Aab_Aab::a01b01_a01b2
@@ -480,7 +476,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b0_a1b2
 
 				case Label::ab_ab::a0b1_a1b2:
@@ -620,7 +616,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b1_a1b2
 
 				case Label::ab_ab::a0b2_a1b0:
@@ -760,7 +756,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b2_a1b0
 
 				case Label::ab_ab::a0b2_a1b1:
@@ -900,7 +896,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b2_a1b1
 
 			  // Aab_Aab::a01b01_a2b01
@@ -1042,7 +1038,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b0_a2b1
 
 				case Label::ab_ab::a0b1_a2b0:
@@ -1182,7 +1178,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b1_a2b0
 
 				case Label::ab_ab::a1b0_a2b1:
@@ -1322,7 +1318,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b0_a2b1
 
 				case Label::ab_ab::a1b1_a2b0:
@@ -1462,7 +1458,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b1_a2b0
 
 			  // Aab_Aab::a01b01_a2b2
@@ -1603,7 +1599,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b0_a2b2
 
 				case Label::ab_ab::a0b1_a2b2:
@@ -1742,7 +1738,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b1_a2b2
 
 				case Label::ab_ab::a1b0_a2b2:
@@ -1881,7 +1877,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b0_a2b2
 
 				case Label::ab_ab::a1b1_a2b2:
@@ -2020,7 +2016,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), rDs_tmp_segments_size, queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b1_a2b2
 
 			  // Aab_Aab::a01b2_a2b01
@@ -2186,7 +2182,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b2_a2b1
 
 				case Label::ab_ab::a0b2_a2b0:
@@ -2350,7 +2346,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b2_a2b0
 
 				case Label::ab_ab::a0b2_a2b1:
@@ -2514,7 +2510,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a0b2_a2b1
 
 				case Label::ab_ab::a1b2_a2b0:
@@ -2678,7 +2674,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						rDs_output.h_array.size(), queue);
 					GPU_Backend::sync(queue);
 
-					rDs_output.download(Ds_result, queue);
+					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b2_a2b0
 
 				default:
@@ -2686,11 +2682,8 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 			} // end switch(label)
 		} // end for label
 
-		// 最后结束前必须等待直至加成功
-		LRI_Cal_Aux::add_Ds_omp_wait_map(Ds_result_thread, Ds_result, lock_Ds_result_add_map, fac_add_Ds);
-	} // end #pragma omp parallel
+	}
 
-	LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map);
 	LRI_Cal_Aux::destroy_Ds_result(Ds_result);
 	#ifdef __MKL_RI
 	mkl_set_num_threads(mkl_threads);
diff --git a/include/RI/ri/gpu/GPU_Data_Output.h b/include/RI/ri/gpu/GPU_Data_Output.h
index c47001f..4c9056c 100644
--- a/include/RI/ri/gpu/GPU_Data_Output.h
+++ b/include/RI/ri/gpu/GPU_Data_Output.h
@@ -9,6 +9,7 @@
 
 #include "../../global/gpu/GPU_Backend.h"
 #include "GPU_Data_Pack.h"
+#include "../LRI_Cal_Aux.h"
 
 #include <cstring>
 #include <type_traits>
@@ -59,7 +60,8 @@ class Output
 
 	void download(
 		std::map<TA, std::map<TAC, Tensor<Tdata_CPU>>> &Ds_result,
-		GPU_Backend::Queue queue) const
+		GPU_Backend::Queue queue,
+		const double fac_add_Ds = 1.0) const
 	{
 		std::vector<Tdata_GPU> h_data(this->totalSize);
 		GPU_Backend::download(
@@ -74,10 +76,7 @@ class Output
 				download_data_impl(D_output, h_data, pack, Same_data_layout());
 
 				Tensor<Tdata_CPU> &D_result = Ds_result[ptrList_A.first][ptrList_B.first];
-				if(D_result.empty())
-					D_result = std::move(D_output);
-				else
-					D_result += D_output;
+				LRI_Cal_Aux::add_Ds(std::move(D_output), D_result, fac_add_Ds);
 			}
 		}
 	}
diff --git a/unittests/ri/LRI-loop3-test.hpp b/unittests/ri/LRI-loop3-test.hpp
index 3aae9f1..faf38e2 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,6 +401,89 @@ namespace LRI_Loop3_Test
 	}
 	
 
+	template<typename Tdata>
+	void test_loop3_accumulation(const RI::LRI_Cal_Mode cal_mode)
+	{
+		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();
+				Treal accumulated_factor = 0, accumulated_magnitude = 0;
+				for(const double factor : {0.0, -1.0, 1.0, 1.0, -1.0, 0.5, 0.0, -0.5, 0.3})
+				{
+					lri.cal_loop3(cal_labels, result, factor);
+					accumulated_factor += Treal(factor);
+					accumulated_magnitude += std::abs(Treal(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[])
 	{
@@ -406,11 +494,14 @@ namespace LRI_Loop3_Test
 			RI::LRI_Cal_Mode::CPU,
 			RI::LRI_Cal_Mode::CPU_fine_grained_lock,
 			#ifdef __GPU_RI
-			, RI::LRI_Cal_Mode::GPU
+			RI::LRI_Cal_Mode::GPU,
 			#endif
 		};
 		for(const RI::LRI_Cal_Mode cal_mode: cal_modes)
+		{
 			test_loop3<Tdata>(cal_mode);
+			test_loop3_accumulation<Tdata>(cal_mode);
+		}
 
 		MPI_Finalize();
 	}

From cdeab5750d748381e661f356027ade1c45d712b0 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Tue, 20 Oct 2026 12:57:33 -0400
Subject: [PATCH 2/9] Fix GPU EXX real-space symmetry filtering

(cherry picked from commit 2f77684417313e78e2da67aeee2ebf9b44f1e9fe)
(cherry picked from commit 87f521fa2a8c159bbe955dd729f2ec4c862ae688)
---
 .../physics/symmetry/Filter_Atom_Symmetry.h   |  8 ++++
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |  3 +-
 unittests/physics/test_exx_gpu_symmetry.cpp   | 46 +++++++++++++++++++
 3 files changed, 56 insertions(+), 1 deletion(-)
 create mode 100644 unittests/physics/test_exx_gpu_symmetry.cpp

diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
index 12698ff..c1f0b9a 100644
--- a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
@@ -45,6 +45,10 @@ 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 computes these contractions before their output atoms are known.
+				case Label::ab_ab::a0b0_a1b1:
+				case Label::ab_ab::a0b0_a1b2:
+					return false;
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
 		}
@@ -100,6 +104,10 @@ 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), respectively.
+				case Label::ab_ab::a0b0_a1b1:
+				case Label::ab_ab::a0b0_a1b2:
+					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 30d7323..c991bdf 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -1553,7 +1553,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);
diff --git a/unittests/physics/test_exx_gpu_symmetry.cpp b/unittests/physics/test_exx_gpu_symmetry.cpp
new file mode 100644
index 0000000..0da1347
--- /dev/null
+++ b/unittests/physics/test_exx_gpu_symmetry.cpp
@@ -0,0 +1,46 @@
+// 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 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 Cell one{{1}};
+    const AtomCell shifted_right{0, one};
+
+    assert(!filter.filter_for1(Label::a0b0_a1b1, 2));
+    assert(!filter.filter_for32(Label::a0b0_a1b1, 2, left, right));
+    assert(filter.filter_for32(Label::a0b0_a1b1, 2, right, left));
+    assert(filter.filter_for32(Label::a0b0_a1b1, 2, left, shifted_right));
+
+    assert(!filter.filter_for1(Label::a0b0_a1b2, 2));
+    assert(!filter.filter_for1(Label::a0b0_a1b2, right));
+    assert(filter.filter_for1(Label::a0b0_a1b2, left));
+    assert(!filter.filter_for32(Label::a0b0_a1b2, 2, left, right));
+    assert(filter.filter_for32(Label::a0b0_a1b2, 2, right, left));
+
+    assert(!filter.filter_for1(Label::a0b0_a2b1, 1));
+    assert(!filter.filter_for32(Label::a0b0_a2b1, 1, intermediate, right));
+    assert(filter.filter_for32(Label::a0b0_a2b1, 1, intermediate, left));
+
+    assert(!filter.filter_for1(Label::a0b0_a2b2, 1));
+    assert(!filter.filter_for32(Label::a0b0_a2b2, 1, intermediate, right));
+    assert(filter.filter_for32(Label::a0b0_a2b2, 1, intermediate, left));
+}

From 4a41b91335013e36c8e623370389cc7243c1dc26 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Wed, 21 Oct 2026 02:04:43 -0400
Subject: [PATCH 3/9] Fix GPU symmetry filtering for all LibRI contraction
 labels

(cherry picked from commit a0db48f56ffadf448c8057788df57562b42cd3c8)
(cherry picked from commit dc1084707f51615ca34e479bb15299e4b1da4540)
---
 .../physics/symmetry/Filter_Atom_Symmetry.h   | 17 +++-
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |  6 +-
 unittests/physics/test_exx_gpu_symmetry.cpp   | 85 +++++++++++++++----
 3 files changed, 87 insertions(+), 21 deletions(-)

diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
index c1f0b9a..77a0f1f 100644
--- a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
@@ -45,9 +45,13 @@ 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 computes these contractions before their output atoms are known.
+			// 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__));
@@ -66,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__));
@@ -104,9 +113,13 @@ 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), respectively.
+			// 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:
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index c991bdf..6fd4949 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -1693,7 +1693,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 +1832,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 +1971,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);
diff --git a/unittests/physics/test_exx_gpu_symmetry.cpp b/unittests/physics/test_exx_gpu_symmetry.cpp
index 0da1347..309336a 100644
--- a/unittests/physics/test_exx_gpu_symmetry.cpp
+++ b/unittests/physics/test_exx_gpu_symmetry.cpp
@@ -15,6 +15,7 @@ int main()
     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}}
     };
@@ -22,25 +23,77 @@ int main()
     const AtomCell left{1, zero};
     const AtomCell right{0, zero};
     const AtomCell intermediate{2, zero};
-    const Cell one{{1}};
     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));
+    }
 
-    assert(!filter.filter_for1(Label::a0b0_a1b1, 2));
-    assert(!filter.filter_for32(Label::a0b0_a1b1, 2, left, right));
-    assert(filter.filter_for32(Label::a0b0_a1b1, 2, right, left));
-    assert(filter.filter_for32(Label::a0b0_a1b1, 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));
+    }
 
-    assert(!filter.filter_for1(Label::a0b0_a1b2, 2));
-    assert(!filter.filter_for1(Label::a0b0_a1b2, right));
-    assert(filter.filter_for1(Label::a0b0_a1b2, left));
-    assert(!filter.filter_for32(Label::a0b0_a1b2, 2, left, right));
-    assert(filter.filter_for32(Label::a0b0_a1b2, 2, right, 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));
+    }
 
-    assert(!filter.filter_for1(Label::a0b0_a2b1, 1));
-    assert(!filter.filter_for32(Label::a0b0_a2b1, 1, intermediate, right));
-    assert(filter.filter_for32(Label::a0b0_a2b1, 1, 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(!filter.filter_for1(Label::a0b0_a2b2, 1));
-    assert(!filter.filter_for32(Label::a0b0_a2b2, 1, intermediate, right));
-    assert(filter.filter_for32(Label::a0b0_a2b2, 1, intermediate, left));
+    assert(checked.size() == RI::Label::array_ab_ab.size());
+    for (const Label label : RI::Label::array_ab_ab)
+    {
+        assert(checked.count(label) == 1);
+    }
 }

From 11de9ee6458ee8acc4ba21a3b961277198aa0d00 Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Fri, 23 Oct 2026 22:46:58 +0800
Subject: [PATCH 4/9] fix(gpu): filter output pairs without truncating
 intermediate sums

In the four a01b2_a2b01 contractions, Aa2 is an intermediate atom
summed into the (Aa01, Ab01) output. Applying filter_for2 to
(Aa01, Aa2) drops valid terms. Apply it once Ab01 is available.

Add standalone CPU/DDLA regression coverage against unfiltered CPU
contractions masked independently by output atom pairs and cells.
Cover all 18 labels and their sum, real/complex tensors, empty and
seeded outputs, and signed/zero/fractional increments.

Validated CPU-only and DDLA v0.0.4/CUDA builds on V100 with 1 and 4
OpenMP threads. The two preceding symmetry fixes are cherry-picked
from maki49; this regression additionally catches the four remaining
incorrect contraction paths.

(cherry picked from commit 2781215e73813de9d58fc45a6087b7b857659f45)
---
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |  12 +-
 unittests/gpu_symmetry/CMakeLists.txt         |  36 +++
 unittests/gpu_symmetry/README.md              |  35 +++
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    | 242 ++++++++++++++++++
 4 files changed, 321 insertions(+), 4 deletions(-)
 create mode 100644 unittests/gpu_symmetry/CMakeLists.txt
 create mode 100644 unittests/gpu_symmetry/README.md
 create mode 100644 unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp

diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 6fd4949..18e0273 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -2106,10 +2106,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 +2271,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 +2436,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 +2601,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/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..6db0f99
--- /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, T(1));
+        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, T(1));
+                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, T(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;
+}

From 27951e74f428decceca7805297838db6130cf8da Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Sat, 24 Oct 2026 00:12:40 +0800
Subject: [PATCH 5/9] test: use real scaling coefficients in symmetry
 regression

---
 unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index 6db0f99..233bc34 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -149,7 +149,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
         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, T(1));
+        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)
@@ -178,7 +178,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                                             + std::to_string(static_cast<int>(mode));
                 std::cout << "CHECK " << context << std::endl;
                 Tensors<T> actual;
-                lri.cal_loop3(labels, actual, T(1));
+                lri.cal_loop3(labels, actual, 1.0);
                 compare(actual, expected, context);
 
                 // Seed every block, including those outside the selected sector.
@@ -192,7 +192,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 double sum = 0;
                 for (const double factor : {-0.5, 0.0, 1.25})
                 {
-                    lri.cal_loop3(labels, actual, T(factor));
+                    lri.cal_loop3(labels, actual, factor);
                     sum += factor;
                     Tensors<T> accumulated = copy_tensors(seed);
                     for (const auto& left : expected)

From 1c2a8905a95a34dc80556e1b82965a50a773f186 Mon Sep 17 00:00:00 2001
From: linpz <linpz@mail.ustc.edu.cn>
Date: Fri, 23 Oct 2026 00:57:44 +0800
Subject: [PATCH 6/9] 1. fix fac_add_Ds bug in LRI-cal_loop3_GPU

(cherry picked from commit fe27bb14d9abbccd9d05296a0833278acc7f34f4)
---
 include/RI/ri/LRI-cal_loop3_GPU.hpp | 34 ++++++++++++-----------------
 include/RI/ri/LRI.h                 |  2 +-
 include/RI/ri/LRI_k-cal_cvc_mo.hpp  | 24 +++++++++++++-------
 include/RI/ri/LRI_k-cal_hartree.hpp | 18 ++++++++++-----
 include/RI/ri/gpu/GPU_Data_Output.h | 22 ++++++++++++++++---
 5 files changed, 62 insertions(+), 38 deletions(-)

diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 30d7323..c05794f 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -33,7 +33,7 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	const std::vector<Label::ab_ab> &labels,
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-	const double fac_add_Ds)
+	const Tdata fac_add_Ds)
 {
 	using namespace Array_Operator;
 
@@ -53,15 +53,16 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	GPU_Backend::Context gpu_context(this->mpi_comm);
 	GPU_Backend::Queue &queue = gpu_context.queue();
 
+	std::map<TA, omp_lock_t> lock_Ds_result_add_map = LRI_Cal_Aux::init_lock_result(labels, this->parallel->list_A, Ds_result);
+
+	for(const Label::ab_ab &label : labels)
 	{
-		for(const Label::ab_ab &label : labels)
+		const std::vector<TA>  list_Aa01_Da = LRI_Cal_Aux::filter_list_map( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].a01, data_wrapper(Label::ab::a).Ds_ab );
+		const std::vector<TAC> list_Ab01_Db = LRI_Cal_Aux::filter_list_map( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].b01, data_wrapper(Label::ab::b).Ds_ab );
+		const std::vector<TAC> list_Aa2_Da  = LRI_Cal_Aux::filter_list_set( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].a2,  data_wrapper(Label::ab::a).index_Ds_ab[0] );
+		const std::vector<TAC> list_Ab2_Db  = LRI_Cal_Aux::filter_list_set( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].b2,  data_wrapper(Label::ab::b).index_Ds_ab[0] );
+		switch(label)
 		{
-			const std::vector<TA>  list_Aa01_Da = LRI_Cal_Aux::filter_list_map( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].a01, data_wrapper(Label::ab::a).Ds_ab );
-			const std::vector<TAC> list_Ab01_Db = LRI_Cal_Aux::filter_list_map( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].b01, data_wrapper(Label::ab::b).Ds_ab );
-			const std::vector<TAC> list_Aa2_Da  = LRI_Cal_Aux::filter_list_set( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].a2,  data_wrapper(Label::ab::a).index_Ds_ab[0] );
-			const std::vector<TAC> list_Ab2_Db  = LRI_Cal_Aux::filter_list_set( this->parallel->list_A[Label_Tools::to_Aab_Aab(label)].b2,  data_wrapper(Label::ab::b).index_Ds_ab[0] );
-			switch(label)
-			{
 
 			  // Aab_Aab::a01b01_a01b01
 
@@ -2677,23 +2678,16 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 					rDs_output.download(Ds_result, queue, fac_add_Ds);
 				} break; // end case a1b2_a2b0
 
-				default:
-					throw std::invalid_argument(std::string(__FILE__)+std::to_string(__LINE__));
-			} // end switch(label)
-		} // end for label
-
-	}
+			default:
+				throw std::invalid_argument(std::string(__FILE__)+std::to_string(__LINE__));
+		} // end switch(label)
+	} // end for label
 
+	LRI_Cal_Aux::destroy_lock(lock_Ds_result_add_map);
 	LRI_Cal_Aux::destroy_Ds_result(Ds_result);
 	#ifdef __MKL_RI
 	mkl_set_num_threads(mkl_threads);
 	#endif
-	//for (const auto &Ds_A : Ds_result)
-	//	for (const auto &Ds_B : Ds_A.second)
-	//	{
-	//		std::ofstream ofs("H_" + std::to_string(Ds_A.first) + "_" + std::to_string(Ds_B.first.first) + ".txt");
-	//		ofs << Ds_B.second << std::endl;
-	//	}
 
 }	// end LRI::cal_loop3()
 
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index d9f7a07..272c75f 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -88,7 +88,7 @@ class LRI
 	void cal_loop3_GPU(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const double fac_add_Ds = 1.0);
+		const Tdata fac_add_Ds = 1.0);
 	#endif
 
 public:
diff --git a/include/RI/ri/LRI_k-cal_cvc_mo.hpp b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
index f76463e..f92123d 100644
--- a/include/RI/ri/LRI_k-cal_cvc_mo.hpp
+++ b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
@@ -131,10 +131,12 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_Csk_ao_mo(
 #endif
 
 #pragma omp parallel for schedule(static) collapse(2)
-    for (const int ik : k_indices)
+    for (int ik_index=0; ik_index<k_indices.size(); ++ik_index)
     {
-        for (const int iat1 : list_IJ)
+		const int ik = k_indices[ik_index];
+        for (int iat1_index=0; iat1_index<list_IJ.size(); ++iat1_index)
         {
+			const int iat1 = list_IJ[iat1_index];
             Tk k = kindex_map.at(ik);
             auto& tensor_ao_mo = Csk_ao_mo.at(ik).at(iat1); // C'^\mu (s,m)[k]
             auto& psi_k = map_psi.at(ik);
@@ -303,12 +305,15 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_onthefly(
 		// 2. calculate CVC_mo_k
 		std::map<int, std::map<int, Tensor<Tdata>>> cvc_mo_k_thread;
 #pragma omp for schedule(dynamic, 64) collapse(3)
-		for (const Tk q: q_list)
+		for (int q_index=0; q_index<q_list.size(); ++q_index)
 		{
-			for (const TA mu: list_I)
+			const Tk q=q_list[q_index];
+			for (int mu_index=0; mu_index<list_I.size(); ++mu_index)
 			{
-				for (const TA nu: list_J)
+				const TA mu = list_I[mu_index];
+				for (int nu_index=0; nu_index<list_J.size(); ++nu_index)
 				{	// 2.1 calculate V(q)_{mu,nu} on-the-fly
+					const TA nu = list_J[nu_index];
 					Tensor<Tdata> Vq_mu_nu;
 					const auto& V_mu = Vs.at(mu);
 					const std::vector<TAC> list_nuR = Divide_Atoms::traversal_atom_period(std::vector<TA>{nu}, this->period);
@@ -432,8 +437,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_hartree_onthefly(
 			const auto& V_mu = Vs.at(mu);
 			auto& Vq_mu_thread = Vq_thread[mu];
 	#pragma omp for schedule(dynamic)
-			for (const TAC& nu_R : list_JR)
+			for (int nu_R_index=0; nu_R_index < list_JR.size(); ++nu_R_index)
 			{
+				const TAC& nu_R = list_JR[nu_R_index];
 				const Tensor<Tdata>& V_mu_nu_R = Global_Func::find(V_mu, nu_R);
 				if (V_mu_nu_R.empty()) continue;
 				const TA nu = nu_R.first;
@@ -457,10 +463,12 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_hartree_onthefly(
 		// 2 calculate CVC_mo_k
 		std::map<int, std::map<int, Tensor<Tdata>>> cvc_mo_k_thread;
 #pragma omp for schedule(dynamic, 64) collapse(2)
-		for (const TA mu : list_I)
+		for (int mu_index=0; mu_index < list_I.size(); ++mu_index)
 		{
-			for (const int k1: k1_indices)
+			const TA& mu = list_J[mu_index];
+			for (int k1_index=0; k1_index<k1_indices.size(); ++k1_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,
 												mu, psi_type[0], psi_type[1], nocc, nvirt, false);
diff --git a/include/RI/ri/LRI_k-cal_hartree.hpp b/include/RI/ri/LRI_k-cal_hartree.hpp
index 4934efe..1ac14cf 100644
--- a/include/RI/ri/LRI_k-cal_hartree.hpp
+++ b/include/RI/ri/LRI_k-cal_hartree.hpp
@@ -63,8 +63,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 			const auto& V_mu = Vs.at(mu);
 			auto& Mu_Vq_nu_thread = Vq_thread[mu];
 	#pragma omp for schedule(dynamic) nowait
-			for (const TAC& nu_R : list_IJR)
+			for (int nu_R_index=0; nu_R_index<list_IJR.size(); nu_R_index++)
 			{
+				const TAC& nu_R = list_IJR[nu_R_index];
 				const Tensor<Tdata>& V_mu_nu_R = Global_Func::find(V_mu, nu_R);
 				if (V_mu_nu_R.empty()) continue;
 				const TA nu = nu_R.first;
@@ -83,8 +84,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 			const auto& Cs_mu = Cs.at(mu);
 			auto& Mu_C_nu_k_thread = Csk_thread[mu];
 	#pragma omp for schedule(dynamic) nowait
-			for (const TAC& nu_R : list_IJR)
+			for (int nu_R_index=0; nu_R_index<list_IJR.size(); nu_R_index++)
 			{
+				const TAC& nu_R = list_IJR[nu_R_index];
 				const Tensor<Tdata>& Cs_mu_nu_R = Global_Func::find(Cs_mu, nu_R);
 				if (Cs_mu_nu_R.empty()) continue;
 				const TA nu = nu_R.first;
@@ -117,10 +119,12 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 		// 3. calculate M_nu =\sum_{uvk} (C^nu_u_v[k] + C^nu*_v_u[k]) D_v_u[k]		
         std::map<TA, Tensor<Tdata>> M_nu_thread;
 #pragma omp for schedule(dynamic) collapse(2)
-		for (const TA v : list_I)
+		for (int v_index=0; v_index<list_I.size(); v_index++)
 		{
-			for (const TA u : list_J)
+			const TA v = list_I[v_index];
+			for (int u_index=0; u_index<list_I.size(); u_index++)
 			{
+				const TA u = list_J[u_index];
 				const std::map<int, Tensor<Tdata>>& D_v_u = Global_Func::find(Ds, v, u);
 				if (D_v_u.empty()) continue;
 				const std::size_t nwt1 = D_v_u.begin()->second.shape[0];
@@ -199,10 +203,12 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 #endif
 
 #pragma omp parallel for schedule(dynamic) collapse(2)
-	for (const TA s : list_I)
+	for (int s_index=0; s_index<list_I.size(); s_index++)
 	{
-		for (const TA t : list_J)
+		const TA s = list_I[s_index];
+		for (int t_index=0; t_index<list_J.size(); t_index++)
 		{
+			const TA t = list_J[t_index];
 			auto& hartree_st = hartree_k.at(s).at(t);
 			const std::map<int, Tensor<Tdata>>& Csk_s_t = Global_Func::find(Csk, s, t);
 			const std::map<int, Tensor<Tdata>>& Csk_t_s = Global_Func::find(Csk, t, s);
diff --git a/include/RI/ri/gpu/GPU_Data_Output.h b/include/RI/ri/gpu/GPU_Data_Output.h
index 4c9056c..c2b7214 100644
--- a/include/RI/ri/gpu/GPU_Data_Output.h
+++ b/include/RI/ri/gpu/GPU_Data_Output.h
@@ -9,7 +9,6 @@
 
 #include "../../global/gpu/GPU_Backend.h"
 #include "GPU_Data_Pack.h"
-#include "../LRI_Cal_Aux.h"
 
 #include <cstring>
 #include <type_traits>
@@ -61,7 +60,7 @@ class Output
 	void download(
 		std::map<TA, std::map<TAC, Tensor<Tdata_CPU>>> &Ds_result,
 		GPU_Backend::Queue queue,
-		const double fac_add_Ds = 1.0) const
+		const Tdata_CPU fac) const
 	{
 		std::vector<Tdata_GPU> h_data(this->totalSize);
 		GPU_Backend::download(
@@ -76,7 +75,24 @@ class Output
 				download_data_impl(D_output, h_data, pack, Same_data_layout());
 
 				Tensor<Tdata_CPU> &D_result = Ds_result[ptrList_A.first][ptrList_B.first];
-				LRI_Cal_Aux::add_Ds(std::move(D_output), D_result, fac_add_Ds);
+				if(D_result.empty())
+				{
+					if(Tdata_CPU(1.0) == fac)
+						D_result = std::move(D_output);
+					else if(Tdata_CPU(-1.0)==fac)
+						D_result = -D_output;
+					else
+						D_result = fac * D_output;
+				}
+				else
+				{
+					if(Tdata_CPU(1.0)==fac)
+						D_result += D_output;
+					else if(Tdata_CPU(-1.0)==fac)
+						D_result -= D_output;
+					else
+						D_result += fac * D_output;
+				}
 			}
 		}
 	}

From c39328132b9b2579ecf665c2ce914ce9ce92b622 Mon Sep 17 00:00:00 2001
From: linpz <linpz@mail.ustc.edu.cn>
Date: Fri, 23 Oct 2026 19:28:00 +0800
Subject: [PATCH 7/9] 1. change fac_add_Ds to Tdata

(cherry picked from commit 6534722c2c4925d521793b3d279052ef3c7be886)
---
 include/RI/global/Global_Func-2-Tensor.h      | 17 +++++++
 include/RI/ri/LRI-cal_loop3_CPU.hpp           |  2 +-
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   |  2 +-
 include/RI/ri/LRI.h                           |  8 +--
 include/RI/ri/LRI.hpp                         |  2 +-
 include/RI/ri/LRI_Cal_Aux.h                   | 49 ++++++++++++-------
 include/RI/ri/LRI_k-cal_cvc_mo.hpp            |  8 +--
 include/RI/ri/LRI_k-cal_hartree.hpp           |  4 +-
 include/RI/ri/gpu/GPU_Data_Output.h           | 20 +-------
 9 files changed, 62 insertions(+), 50 deletions(-)

diff --git a/include/RI/global/Global_Func-2-Tensor.h b/include/RI/global/Global_Func-2-Tensor.h
index 68ba301..586d704 100644
--- a/include/RI/global/Global_Func-2-Tensor.h
+++ b/include/RI/global/Global_Func-2-Tensor.h
@@ -39,6 +39,23 @@ namespace Global_Func
 			ptr_out[i] = Global_Func::convert<Tout>(ptr_in[i]);
 		return t_out;
 	}
+
+	template<typename T>
+	struct Tensor_Tdata;
+
+	template<typename Tdata>
+	struct Tensor_Tdata<Tensor<Tdata>>
+	{
+		using type = Tdata;
+	};
+
+	template<typename Tkey, typename Tvalue, typename TComparator, typename TAllocator>
+	struct Tensor_Tdata<std::map<Tkey, Tvalue, TComparator, TAllocator>>
+		: Tensor_Tdata<Tvalue>
+	{};
+
+	template<typename T>
+	using Tensor_Tdata_t = typename Tensor_Tdata<T>::type;
 }
 
 }
\ No newline at end of file
diff --git a/include/RI/ri/LRI-cal_loop3_CPU.hpp b/include/RI/ri/LRI-cal_loop3_CPU.hpp
index 7027469..e702ebb 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -22,7 +22,7 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	const std::vector<Label::ab_ab> &labels,
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-	const double fac_add_Ds)
+	const Tdata fac_add_Ds)
 {
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI-cal_loop3_CPU_fine_grained_lock.hpp b/include/RI/ri/LRI-cal_loop3_CPU_fine_grained_lock.hpp
index b915c4c..ac08b09 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU_fine_grained_lock.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU_fine_grained_lock.hpp
@@ -22,7 +22,7 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	const std::vector<Label::ab_ab> &labels,
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-	const double fac_add_Ds)
+	const Tdata fac_add_Ds)
 {
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index 272c75f..c3db89d 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -72,23 +72,23 @@ class LRI
 	void cal_loop3(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const double fac_add_Ds = 1.0);
+		const Tdata fac_add_Ds = Tdata(1.0));
 
 	void cal_loop3_CPU(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const double fac_add_Ds = 1.0);
+		const Tdata fac_add_Ds = Tdata(1.0));
 
 	void cal_loop3_CPU_fine_grained_lock(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const double fac_add_Ds = 1.0);
+		const Tdata fac_add_Ds = Tdata(1.0));
 	
 	#ifdef __GPU_RI
 	void cal_loop3_GPU(
 		const std::vector<Label::ab_ab> &labels,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-		const Tdata fac_add_Ds = 1.0);
+		const Tdata fac_add_Ds = Tdata(1.0));
 	#endif
 
 public:
diff --git a/include/RI/ri/LRI.hpp b/include/RI/ri/LRI.hpp
index 1c75440..127e656 100644
--- a/include/RI/ri/LRI.hpp
+++ b/include/RI/ri/LRI.hpp
@@ -32,7 +32,7 @@ template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
 void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3(
 	const std::vector<Label::ab_ab> &labels,
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
-	const double fac_add_Ds)
+	const Tdata fac_add_Ds)
 {
 	if(this->cal_mode == LRI_Cal_Mode::CPU)
 		this->cal_loop3_CPU(labels, Ds_result, fac_add_Ds);
diff --git a/include/RI/ri/LRI_Cal_Aux.h b/include/RI/ri/LRI_Cal_Aux.h
index 81bd2ed..1eae641 100644
--- a/include/RI/ri/LRI_Cal_Aux.h
+++ b/include/RI/ri/LRI_Cal_Aux.h
@@ -9,6 +9,7 @@
 #include "Label.h"
 #include "../global/Map_Operator.h"
 #include "../global/Array_Operator.h"
+#include "../global/Global_Func-2-Tensor.h"
 #include "../parallel/Parallel_LRI.h"
 
 #include <map>
@@ -61,43 +62,53 @@ namespace LRI_Cal_Aux
 	inline void add_Ds(
 		Tensor<Tdata> &&D_add,
 		Tensor<Tdata> &D_result,
-		const double fac = 1.0)
+		const Tdata fac = Tdata(1.0))
 	{
 		if(D_result.empty())
 		{
-			if(1.0==fac)
+			if(Tdata(1.0)==fac)
 				D_result = std::move(D_add);
-			else if(-1.0==fac)
+			else if(Tdata(-1.0)==fac)
 				D_result = -D_add;
 			else
-				D_result = Tdata(fac) * D_add;
+				D_result = fac * D_add;
 		}
 		else
 		{
-			if(1.0==fac)
+			if(Tdata(1.0)==fac)
 				D_result += D_add;
-			else if(-1.0==fac)
+			else if(Tdata(-1.0)==fac)
 				D_result -= D_add;
 			else
-				D_result += Tdata(fac) * D_add;
+				Blas_Interface::axpy(fac, D_add, D_result);
 		}
 	}
 
-	// D_result = D_add * scale  (init)
-	// or D_result += D_add * scale (accumulate via BLAS axpy)
+	// D_result = D_add * fac  (init)
+	// or D_result += D_add * fac (accumulate via BLAS axpy)
 	template<typename Tdata>
 	inline void add_Ds(
 		const Tensor<Tdata> &D_add,
 		Tensor<Tdata> &D_result,
-		const Tdata scale)
+		const Tdata fac = Tdata(1.0))
 	{
 		if (D_result.empty())
 		{
-			D_result = D_add * scale;
+			if(Tdata(1.0)==fac)
+				D_result = D_add;
+			else if(Tdata(-1.0)==fac)
+				D_result = -D_add;
+			else
+				D_result = fac * D_add;
 		}
 		else
 		{
-			Blas_Interface::axpy(scale, D_add, D_result);
+			if(Tdata(1.0)==fac)
+				D_result += D_add;
+			else if(Tdata(-1.0)==fac)
+				D_result -= D_add;
+			else
+				Blas_Interface::axpy(fac, D_add, D_result);
 		}
 	}
 
@@ -105,11 +116,11 @@ namespace LRI_Cal_Aux
 	void add_Ds(
 		std::map<Tkey, Tvalue> &&Ds_add,
 		std::map<Tkey, Tvalue> &Ds_result,
-		const double fac = 1.0)
+		const Global_Func::Tensor_Tdata_t<Tvalue> fac = Global_Func::Tensor_Tdata_t<Tvalue>(1.0))
 	{
 		if(Ds_add.empty())
 			return;
-		if(Ds_result.empty() && 1.0==fac)
+		if(Ds_result.empty() && Global_Func::Tensor_Tdata_t<Tvalue>(1.0)==fac)
 			Ds_result = std::move(Ds_add);
 		else
 		{
@@ -152,7 +163,7 @@ namespace LRI_Cal_Aux
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result_thread,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 		omp_lock_t &lock_Ds_result_add,
-		const double &fac)
+		const Tdata fac = Tdata(1.0))
 	{
 		if( !Ds_result_thread.empty() && omp_test_lock(&lock_Ds_result_add) )
 		{
@@ -167,7 +178,7 @@ namespace LRI_Cal_Aux
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &&Ds_result_thread,
 		std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 		omp_lock_t &lock_Ds_result_add,
-		const double &fac)
+		const Tdata fac = Tdata(1.0))
 	{
 		if(!Ds_result_thread.empty())
 		{
@@ -185,7 +196,7 @@ namespace LRI_Cal_Aux
 		std::map<Tkey, Tvalue, TComparator>& Ds_result_thread,
 		std::map<Tkey, Tvalue, TComparator>& Ds_result,
 		std::map<Tkey, omp_lock_t> &lock_Ds_result_add_map,
-		const double &fac)
+		const Global_Func::Tensor_Tdata_t<Tvalue> fac = Global_Func::Tensor_Tdata_t<Tvalue>(1.0))
 	{
 		for(auto ptr=Ds_result_thread.begin(); ptr!=Ds_result_thread.end(); )
 		{
@@ -208,7 +219,7 @@ namespace LRI_Cal_Aux
 		std::map<Tkey0, std::map<Tkey1, Tvalue, TComparator1>, TComparator0> &Ds_result_thread,
 		std::map<Tkey0, std::map<Tkey1, Tvalue, TComparator1>, TComparator0> &Ds_result,
 		std::map<Tkey0, std::map<Tkey1, omp_lock_t>> &lock_Ds_result_add_map,
-		const double &fac)
+		const Global_Func::Tensor_Tdata_t<Tvalue> fac = Global_Func::Tensor_Tdata_t<Tvalue>(1.0))
 	{
 		for(auto ptr=Ds_result_thread.begin(); ptr!=Ds_result_thread.end(); )
 		{
@@ -228,7 +239,7 @@ namespace LRI_Cal_Aux
 		std::map<Tkey, Tvalue, TComparator>& Ds_result_thread,
 		std::map<Tkey, Tvalue, TComparator>& Ds_result,
 		std::map<Tkey, Tlock> &lock_Ds_result_add_map,
-		const double &fac)
+		const Global_Func::Tensor_Tdata_t<Tvalue> fac = Global_Func::Tensor_Tdata_t<Tvalue>(1.0))
 	{
 		if(Ds_result_thread.empty())
 			return;
diff --git a/include/RI/ri/LRI_k-cal_cvc_mo.hpp b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
index f92123d..d540770 100644
--- a/include/RI/ri/LRI_k-cal_cvc_mo.hpp
+++ b/include/RI/ri/LRI_k-cal_cvc_mo.hpp
@@ -133,9 +133,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_Csk_ao_mo(
 #pragma omp parallel for schedule(static) collapse(2)
     for (int ik_index=0; ik_index<k_indices.size(); ++ik_index)
     {
-		const int ik = k_indices[ik_index];
         for (int iat1_index=0; iat1_index<list_IJ.size(); ++iat1_index)
         {
+			const int ik = k_indices[ik_index];
 			const int iat1 = list_IJ[iat1_index];
             Tk k = kindex_map.at(ik);
             auto& tensor_ao_mo = Csk_ao_mo.at(ik).at(iat1); // C'^\mu (s,m)[k]
@@ -307,12 +307,12 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_onthefly(
 #pragma omp for schedule(dynamic, 64) collapse(3)
 		for (int q_index=0; q_index<q_list.size(); ++q_index)
 		{
-			const Tk q=q_list[q_index];
 			for (int mu_index=0; mu_index<list_I.size(); ++mu_index)
 			{
-				const TA mu = list_I[mu_index];
 				for (int nu_index=0; nu_index<list_J.size(); ++nu_index)
 				{	// 2.1 calculate V(q)_{mu,nu} on-the-fly
+					const Tk q=q_list[q_index];
+					const TA mu = list_I[mu_index];
 					const TA nu = list_J[nu_index];
 					Tensor<Tdata> Vq_mu_nu;
 					const auto& V_mu = Vs.at(mu);
@@ -465,9 +465,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvc_mo_k_hartree_onthefly(
 #pragma omp for schedule(dynamic, 64) collapse(2)
 		for (int mu_index=0; mu_index < list_I.size(); ++mu_index)
 		{
-			const TA& mu = list_J[mu_index];
 			for (int k1_index=0; k1_index<k1_indices.size(); ++k1_index)
 			{
+				const TA& mu = list_J[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 1ac14cf..9705106 100644
--- a/include/RI/ri/LRI_k-cal_hartree.hpp
+++ b/include/RI/ri/LRI_k-cal_hartree.hpp
@@ -121,9 +121,9 @@ 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++)
 		{
-			const TA v = list_I[v_index];
 			for (int u_index=0; u_index<list_I.size(); u_index++)
 			{
+				const TA v = list_I[v_index];
 				const TA u = list_J[u_index];
 				const std::map<int, Tensor<Tdata>>& D_v_u = Global_Func::find(Ds, v, u);
 				if (D_v_u.empty()) continue;
@@ -205,9 +205,9 @@ LRI_k<TA, Tcell, Ndim, Tdata>::cal_cvcd_k_hartree(
 #pragma omp parallel for schedule(dynamic) collapse(2)
 	for (int s_index=0; s_index<list_I.size(); s_index++)
 	{
-		const TA s = list_I[s_index];
 		for (int t_index=0; t_index<list_J.size(); t_index++)
 		{
+			const TA s = list_I[s_index];
 			const TA t = list_J[t_index];
 			auto& hartree_st = hartree_k.at(s).at(t);
 			const std::map<int, Tensor<Tdata>>& Csk_s_t = Global_Func::find(Csk, s, t);
diff --git a/include/RI/ri/gpu/GPU_Data_Output.h b/include/RI/ri/gpu/GPU_Data_Output.h
index c2b7214..ec0bcea 100644
--- a/include/RI/ri/gpu/GPU_Data_Output.h
+++ b/include/RI/ri/gpu/GPU_Data_Output.h
@@ -9,6 +9,7 @@
 
 #include "../../global/gpu/GPU_Backend.h"
 #include "GPU_Data_Pack.h"
+#include "../LRI_Cal_Aux.h"
 
 #include <cstring>
 #include <type_traits>
@@ -75,24 +76,7 @@ class Output
 				download_data_impl(D_output, h_data, pack, Same_data_layout());
 
 				Tensor<Tdata_CPU> &D_result = Ds_result[ptrList_A.first][ptrList_B.first];
-				if(D_result.empty())
-				{
-					if(Tdata_CPU(1.0) == fac)
-						D_result = std::move(D_output);
-					else if(Tdata_CPU(-1.0)==fac)
-						D_result = -D_output;
-					else
-						D_result = fac * D_output;
-				}
-				else
-				{
-					if(Tdata_CPU(1.0)==fac)
-						D_result += D_output;
-					else if(Tdata_CPU(-1.0)==fac)
-						D_result -= D_output;
-					else
-						D_result += fac * D_output;
-				}
+				LRI_Cal_Aux::add_Ds(std::move(D_output), D_result, fac);
 			}
 		}
 	}

From a02f626d7955d3e636fba814401b51e34c21c38b Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Sat, 24 Oct 2026 14:44:37 +0800
Subject: [PATCH 8/9] Fix atom-list indexing in indexed Hartree loops

Use the right atom list to bound the density loop and the left atom list to select mu in the CVC loop. Add scalar-reference regressions for distinct and unequal atom partitions.
---
 include/RI/ri/LRI_k-cal_cvc_mo.hpp     |   2 +-
 include/RI/ri/LRI_k-cal_hartree.hpp    |   2 +-
 unittests/ri/LRI-k-atom-lists-test.hpp | 104 +++++++++++++++++++++++++
 3 files changed, 106 insertions(+), 2 deletions(-)
 create mode 100644 unittests/ri/LRI-k-atom-lists-test.hpp

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/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;
+}
+}

From 2cf829054e2bbbc3a5ebd2bb9e2bb6d38e6060fa Mon Sep 17 00:00:00 2001
From: Jiacheng Xu <13862180016@163.com>
Date: Sat, 24 Oct 2026 14:44:43 +0800
Subject: [PATCH 9/9] Keep const tensor inputs independent of accumulators

Copy the input when initializing an empty result with unit coefficient. Tensor assignment shares storage, so subsequent in-place additions otherwise change the input and other accumulators. Cover independent outputs and complex accumulation coefficients.
---
 include/RI/ri/LRI_Cal_Aux.h     |  2 +-
 unittests/ri/LRI-loop3-test.hpp | 45 +++++++++++++++++++++++++++++----
 2 files changed, 41 insertions(+), 6 deletions(-)

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/unittests/ri/LRI-loop3-test.hpp b/unittests/ri/LRI-loop3-test.hpp
index faf38e2..0c871e0 100644
--- a/unittests/ri/LRI-loop3-test.hpp
+++ b/unittests/ri/LRI-loop3-test.hpp
@@ -402,7 +402,33 @@ namespace LRI_Loop3_Test
 	
 
 	template<typename Tdata>
-	void test_loop3_accumulation(const RI::LRI_Cal_Mode cal_mode)
+	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;
@@ -452,12 +478,13 @@ namespace LRI_Loop3_Test
 				// Tensor copies share storage unless copy() is used explicitly.
 				if(seeded)
 					result[0][key] = seed.copy();
-				Treal accumulated_factor = 0, accumulated_magnitude = 0;
-				for(const double factor : {0.0, -1.0, 1.0, 1.0, -1.0, 0.5, 0.0, -0.5, 0.3})
+				Tdata accumulated_factor = 0;
+				Treal accumulated_magnitude = 0;
+				for(const Tdata factor : factors)
 				{
 					lri.cal_loop3(cal_labels, result, factor);
-					accumulated_factor += Treal(factor);
-					accumulated_magnitude += std::abs(Treal(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)
@@ -490,6 +517,8 @@ namespace LRI_Loop3_Test
 		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,
@@ -501,6 +530,12 @@ namespace LRI_Loop3_Test
 		{
 			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();