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/5] 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 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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();