From 6c6c8711894f6379d47448f71fd718170995d2fe Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Mon, 9 Nov 2026 10:34:23 -0500
Subject: [PATCH 1/5] Separate GPU symmetry filtering from CPU loop conventions

---
 include/RI/physics/Exx.hpp                    |  8 +-
 include/RI/physics/GW.hpp                     |  8 +-
 include/RI/physics/RPA.hpp                    |  8 +-
 .../physics/symmetry/Filter_Atom_Symmetry.h   | 21 ------
 .../symmetry/Filter_Atom_Symmetry_GPU.h       | 74 +++++++++++++++++++
 include/RI/ri/LRI-cal_loop3_CPU.hpp           |  7 ++
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   |  7 ++
 include/RI/ri/LRI-cal_loop3_GPU.hpp           | 59 ++++++++-------
 include/RI/ri/LRI.h                           |  3 +
 unittests/gpu_symmetry/README.md              | 40 ++++++++++
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    | 13 +++-
 11 files changed, 179 insertions(+), 69 deletions(-)
 create mode 100644 include/RI/physics/symmetry/Filter_Atom_Symmetry_GPU.h

diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp
index 43b505c..42ff97d 100644
--- a/include/RI/physics/Exx.hpp
+++ b/include/RI/physics/Exx.hpp
@@ -9,7 +9,6 @@
 #include "../ri/Cell_Nearest.h"
 #include "../ri/Label.h"
 #include "../global/Map_Operator.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
 
 #include <cassert>
 
@@ -41,11 +40,8 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
-	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->lri.flag_symmetry = flag_symmetry;
+	this->lri.irreducible_sector = irreducible_sector;
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/GW.hpp b/include/RI/physics/GW.hpp
index 815d5f7..8733ee7 100644
--- a/include/RI/physics/GW.hpp
+++ b/include/RI/physics/GW.hpp
@@ -6,7 +6,6 @@
 
 #include "GW.h"
 #include "../ri/Label.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
 
 namespace RI
 {
@@ -36,11 +35,8 @@ void GW<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
-	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->lri.flag_symmetry = flag_symmetry;
+	this->lri.irreducible_sector = irreducible_sector;
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/RPA.hpp b/include/RI/physics/RPA.hpp
index 4f280f5..454e5cf 100644
--- a/include/RI/physics/RPA.hpp
+++ b/include/RI/physics/RPA.hpp
@@ -7,7 +7,6 @@
 
 #include "RPA.h"
 #include "../ri/Label.h"
-#include "./symmetry/Filter_Atom_Symmetry.h"
 
 #include <cassert>
 
@@ -37,11 +36,8 @@ void RPA<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, irreducible_sector);
-	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->lri.flag_symmetry = flag_symmetry;
+	this->lri.irreducible_sector = irreducible_sector;
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
index 77a0f1f..12698ff 100644
--- a/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry.h
@@ -45,14 +45,6 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b2_a2b0:
 				case Label::ab_ab::a1b2_a2b1:
 					return !this->symmetry.is_Aa_in_irreducible_sector(A1);
-			// GPU visits Aa01 before the output Aa2 is known in these cases.
-				case Label::ab_ab::a0b0_a1b1:
-				case Label::ab_ab::a0b1_a1b0:
-				case Label::ab_ab::a0b0_a1b2:
-				case Label::ab_ab::a0b1_a1b2:
-				case Label::ab_ab::a0b2_a1b0:
-				case Label::ab_ab::a0b2_a1b1:
-					return false;
 			default:
 				throw std::invalid_argument("label "+Label_Tools::get_name(label)+" error in "+std::string(__FILE__)+" line "+std::to_string(__LINE__));
 		}
@@ -70,11 +62,6 @@ 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__));
@@ -113,14 +100,6 @@ class Filter_Atom_Symmetry: public Filter_Atom<TA, std::pair<TA, TC>>
 				case Label::ab_ab::a1b0_a2b2:
 				case Label::ab_ab::a1b1_a2b2:
 					return !this->symmetry.in_irreducible_sector(A1, A3);
-			// GPU output is (Aa2, Ab2) or (Aa2, Ab01); Aa01 is intermediate.
-				case Label::ab_ab::a0b0_a1b1:
-				case Label::ab_ab::a0b1_a1b0:
-				case Label::ab_ab::a0b0_a1b2:
-				case Label::ab_ab::a0b1_a1b2:
-				case Label::ab_ab::a0b2_a1b0:
-				case Label::ab_ab::a0b2_a1b1:
-					return !this->symmetry.in_irreducible_sector(A2, A3);
 			// a01b2_a2b01:                                        Aa01,          Ab01,          Aa2
 			default:
 					return false;
diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry_GPU.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry_GPU.h
new file mode 100644
index 0000000..2275216
--- /dev/null
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry_GPU.h
@@ -0,0 +1,74 @@
+#pragma once
+
+#include "../../ri/Filter_Atom.h"
+#include "../../ri/Label_Tools.h"
+#include "Symmetry_Filter.h"
+
+namespace RI
+{
+
+// Arguments follow the GPU loop nesting, from outermost to innermost.
+// Only output atoms may be screened; intermediate atoms are summed over.
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_GPU: public Filter_Atom<TA, std::pair<TA, TC>>
+{
+  public:
+    using TAC = std::pair<TA, TC>;
+    using Group = Label::Aab_Aab;
+
+    Filter_Atom_Symmetry_GPU(
+        const TC &period,
+        const std::map<std::pair<TA,TA>, std::set<TC>> &irsec)
+        : symmetry(period, irsec) {}
+
+    bool filter_for1(const Label::ab_ab &label, const TA &A1) const override
+    {
+        switch (Label_Tools::to_Aab_Aab(label))
+        {
+            case Group::a01b01_a01b01: // Aa01 is intermediate.
+            case Group::a01b01_a01b2:
+                return false;
+            default: // Aa01 is the left output atom.
+                return !symmetry.is_Aa_in_irreducible_sector(A1);
+        }
+    }
+
+    bool filter_for1(const Label::ab_ab &, const TAC &A1) const override
+    {
+        // The first loop with a cell-bearing atom visits output Ab01.
+        return !symmetry.is_Ab_in_irreducible_sector(A1.first);
+    }
+
+    bool filter_for2(const Label::ab_ab &label, const TA &A1, const TAC &A2) const override
+    {
+        // Only the final stage of a01b2_a2b01 has both outputs at level 2.
+        if (Label_Tools::to_Aab_Aab(label) == Group::a01b2_a2b01)
+            return !symmetry.in_irreducible_sector(A1, A2);
+        return false;
+    }
+
+    bool filter_for32(const Label::ab_ab &label,
+                      const TA &A1, const TAC &A2, const TAC &A3) const override
+    {
+        switch (Label_Tools::to_Aab_Aab(label))
+        {
+            case Group::a01b01_a01b01: // Aa01, Aa2, Ab2
+                return !symmetry.in_irreducible_sector(A2, A3);
+            case Group::a01b01_a01b2: // Aa01, Ab01, Aa2
+                return !symmetry.in_irreducible_sector(A3, A2);
+            case Group::a01b01_a2b01: // Aa01, Ab01, Ab2
+            case Group::a01b2_a2b01: // Aa01, Aa2, Ab01
+                return !symmetry.in_irreducible_sector(A1, A3);
+            case Group::a01b01_a2b2: // Aa01, Ab01, Ab2
+                return !symmetry.in_irreducible_sector(A1, A2);
+            default:
+                throw std::invalid_argument("Unknown GPU symmetry label");
+        }
+    }
+
+    // filter_for31 builds intermediates and deliberately never screens pairs.
+    // All other callbacks retain the base class's no-op behavior.
+    Symmetry_Filter<TA,TC,Tdata> symmetry;
+};
+
+}
diff --git a/include/RI/ri/LRI-cal_loop3_CPU.hpp b/include/RI/ri/LRI-cal_loop3_CPU.hpp
index e702ebb..f618e25 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -6,6 +6,7 @@
 #pragma once
 
 #include "LRI.h"
+#include "../physics/symmetry/Filter_Atom_Symmetry.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Tensor_Multiply.h"
@@ -24,6 +25,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
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 ac08b09..308fa87 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
@@ -6,6 +6,7 @@
 #pragma once
 
 #include "LRI.h"
+#include "../physics/symmetry/Filter_Atom_Symmetry.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Tensor_Multiply.h"
@@ -24,6 +25,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 5ca1304..24090c2 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -8,6 +8,7 @@
 #ifdef __GPU_RI
 
 #include "LRI.h"
+#include "../physics/symmetry/Filter_Atom_Symmetry_GPU.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Map_Operator.h"
@@ -35,6 +36,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+
 	using namespace Array_Operator;
 
 	constexpr std::size_t memory_limit = std::size_t(5)*1024*1024*1024;
@@ -378,7 +385,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -431,7 +438,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -518,7 +525,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -571,7 +578,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -658,7 +665,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -711,7 +718,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -798,7 +805,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -851,7 +858,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -940,7 +947,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -1080,7 +1087,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -1220,7 +1227,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -1360,7 +1367,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -1553,7 +1560,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
 									// The output block is (Aa01, Ab01); Ab2 is only an intermediate atom.
-									if (this->filter_atom->filter_for32(label, Aa01, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1692,7 +1699,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, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1831,7 +1838,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, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1970,7 +1977,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, Ab2, Ab01))	continue;
+									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -2060,7 +2067,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -2109,8 +2116,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									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;
+									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -2225,7 +2231,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -2274,8 +2280,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									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;
+									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -2390,7 +2395,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -2439,8 +2444,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									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;
+									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -2555,7 +2559,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -2604,8 +2608,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									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;
+									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index c3db89d..e1c154d 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -96,6 +96,9 @@ class LRI
 		parallel = std::make_shared<Parallel_LRI_Equally<TA,Tcell,Ndim,Tdata>>();
 	std::unordered_map< Label::ab, RI_Tools::T_filter_func<Tdata> >
 		filter_funcs;
+	// Construct the symmetry filter for the selected backend at contraction time.
+	bool flag_symmetry = false;
+	std::map<std::pair<TA,TA>, std::set<TC>> irreducible_sector;
 	std::shared_ptr<Filter_Atom<TA,TAC>> filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
 public:		// private:
diff --git a/unittests/gpu_symmetry/README.md b/unittests/gpu_symmetry/README.md
index 167f271..6d0e8fe 100644
--- a/unittests/gpu_symmetry/README.md
+++ b/unittests/gpu_symmetry/README.md
@@ -33,3 +33,43 @@ 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.
+
+## GPU symmetry callback contract
+
+Arguments follow the enclosing loops from outermost to innermost. `filter_for31`
+always returns false (intermediate construction); `filter_for32` rejects only
+output pairs outside the sector. **Bold argument tuples changed in this refactor.**
+
+| Label | Family | `filter_for31` arguments | `filter_for32` arguments | Output pair screened by `filter_for32` |
+| --- | --- | --- | --- | --- |
+| `a0b0_a1b1` | `a01b01_a01b01` | `(Aa01, Ab01, Ab2)` | `(Aa01, Aa2, Ab2)` | `(Aa2, Ab2)` |
+| `a0b1_a1b0` | `a01b01_a01b01` | `(Aa01, Ab01, Ab2)` | `(Aa01, Aa2, Ab2)` | `(Aa2, Ab2)` |
+| `a0b0_a1b2` | `a01b01_a01b2` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Ab01, Aa2)`** | `(Aa2, Ab01)` |
+| `a0b1_a1b2` | `a01b01_a01b2` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Ab01, Aa2)`** | `(Aa2, Ab01)` |
+| `a0b2_a1b0` | `a01b01_a01b2` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Ab01, Aa2)`** | `(Aa2, Ab01)` |
+| `a0b2_a1b1` | `a01b01_a01b2` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Ab01, Aa2)`** | `(Aa2, Ab01)` |
+| `a0b0_a2b1` | `a01b01_a2b01` | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01, Ab2)` | `(Aa01, Ab2)` |
+| `a0b1_a2b0` | `a01b01_a2b01` | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01, Ab2)` | `(Aa01, Ab2)` |
+| `a1b0_a2b1` | `a01b01_a2b01` | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01, Ab2)` | `(Aa01, Ab2)` |
+| `a1b1_a2b0` | `a01b01_a2b01` | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01, Ab2)` | `(Aa01, Ab2)` |
+| `a0b0_a2b2` | `a01b01_a2b2` | `(Aa01, Aa2, Ab2)` | **`(Aa01, Ab01, Ab2)`** | `(Aa01, Ab01)` |
+| `a0b1_a2b2` | `a01b01_a2b2` | `(Aa01, Aa2, Ab2)` | **`(Aa01, Ab01, Ab2)`** | `(Aa01, Ab01)` |
+| `a1b0_a2b2` | `a01b01_a2b2` | `(Aa01, Aa2, Ab2)` | **`(Aa01, Ab01, Ab2)`** | `(Aa01, Ab01)` |
+| `a1b1_a2b2` | `a01b01_a2b2` | `(Aa01, Aa2, Ab2)` | **`(Aa01, Ab01, Ab2)`** | `(Aa01, Ab01)` |
+| `a1b2_a2b1` | `a01b2_a2b01` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01)` |
+| `a0b2_a2b0` | `a01b2_a2b01` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01)` |
+| `a0b2_a2b1` | `a01b2_a2b01` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01)` |
+| `a1b2_a2b0` | `a01b2_a2b01` | **`(Ab01, Ab2, Aa01)`** | **`(Aa01, Aa2, Ab01)`** | `(Aa01, Ab01)` |
+
+`filter_for1(Aa01)` is a no-op for the first two families (Aa01 is intermediate),
+and checks left-output membership for the other three. `filter_for1(Ab01)` checks
+right-output membership wherever it occurs. All `filter_for2` calls are no-ops
+except the final `(Aa01, Ab01)` stage of `a01b2_a2b01`, which checks that output
+pair. In that family's second stage, the former inner
+`filter_for2(Aa01, Ab01)` check is absorbed into
+**`filter_for32(Aa01, Aa2, Ab01)`**; `Aa2` remains an unrestricted summation atom.
+
+For two cell-bearing output atoms the relative cell is `(Rb - Ra) % period`;
+with `Aa01` at the origin it is `Rb % period`. An empty sector retains the existing
+unrestricted behavior. The numerical regression also disables symmetry after each
+filtered run to verify replacement of the previous backend's filter.
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index 233bc34..c4541c8 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -146,7 +146,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
             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.flag_symmetry = false;
         lri.cal_mode = RI::LRI_Cal_Mode::CPU;
         Tensors<T> full;
         lri.cal_loop3(labels, full, 1.0);
@@ -169,7 +169,8 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                     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);
+            lri.flag_symmetry = true;
+            lri.irreducible_sector = sector;
             for (const auto mode : modes)
             {
                 lri.cal_mode = mode;
@@ -204,6 +205,14 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                         }
                     compare(actual, accumulated, context + " factor=" + std::to_string(factor));
                 }
+
+                // Disable symmetry after a restricted run on this same backend.
+                // This must replace the previous backend-specific filter.
+                lri.flag_symmetry = false;
+                Tensors<T> unrestricted;
+                lri.cal_loop3(labels, unrestricted, 1.0);
+                compare(unrestricted, full, context + " symmetry disabled");
+                lri.flag_symmetry = true;
             }
         }
         std::cout << "PASS " << scalar << " label=" << name << std::endl;

From f6273c96598e4eb32c7d1b4f16d3a32c4b9e8135 Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Tue, 10 Nov 2026 01:21:54 -0500
Subject: [PATCH 2/5] Reuse backend symmetry filters bound to the current LRI

---
 include/RI/ri/LRI-cal_loop3_CPU.hpp           |  6 +----
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   |  6 +----
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |  6 +----
 include/RI/ri/LRI.h                           |  4 ++++
 include/RI/ri/LRI.hpp                         | 23 +++++++++++++++++++
 unittests/gpu_symmetry/README.md              |  6 +++++
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    | 23 +++++++++++++++++++
 7 files changed, 59 insertions(+), 15 deletions(-)

diff --git a/include/RI/ri/LRI-cal_loop3_CPU.hpp b/include/RI/ri/LRI-cal_loop3_CPU.hpp
index f618e25..6fed27d 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -25,11 +25,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->template prepare_filter_atom<Filter_Atom_Symmetry<TA,TC,Tdata>>();
 
 	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 308fa87..46e80d7 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
@@ -25,11 +25,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->template prepare_filter_atom<Filter_Atom_Symmetry<TA,TC,Tdata>>();
 
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 24090c2..0d7433f 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -36,11 +36,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	this->template prepare_filter_atom<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>();
 
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index e1c154d..df07fe9 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -101,6 +101,10 @@ class LRI
 	std::map<std::pair<TA,TA>, std::set<TC>> irreducible_sector;
 	std::shared_ptr<Filter_Atom<TA,TAC>> filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
+private:
+	template<typename Tfilter>
+	void prepare_filter_atom();
+
 public:		// private:
 	TC period;
 	MPI_Comm mpi_comm;
diff --git a/include/RI/ri/LRI.hpp b/include/RI/ri/LRI.hpp
index 127e656..eddfbf7 100644
--- a/include/RI/ri/LRI.hpp
+++ b/include/RI/ri/LRI.hpp
@@ -8,6 +8,7 @@
 #include "LRI.h"
 #include "../ri/Label.h"
 #include <limits>
+#include <typeinfo>
 
 namespace RI
 {
@@ -28,6 +29,28 @@ LRI<TA,Tcell,Ndim,Tdata>::LRI()
 		this->period[i] = std::numeric_limits<Tcell>::max()/4;		// /4 for not out of range when Array_Operator::operator%
 }
 
+template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
+template<typename Tfilter>
+void LRI<TA,Tcell,Ndim,Tdata>::prepare_filter_atom()
+{
+	if (!this->flag_symmetry)
+	{
+		if (!this->filter_atom || typeid(*this->filter_atom) != typeid(Filter_Atom<TA,TAC>))
+			this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+		return;
+	}
+
+	if (this->filter_atom && typeid(*this->filter_atom) == typeid(Tfilter))
+	{
+		const auto &filter = static_cast<const Tfilter &>(*this->filter_atom);
+		// A copied LRI must bind the filter to its own configuration.
+		if (&filter.symmetry.period == &this->period
+			&& &filter.symmetry.irreducible_sector_ == &this->irreducible_sector)
+			return;
+	}
+	this->filter_atom = std::make_shared<Tfilter>(this->period, this->irreducible_sector);
+}
+
 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,
diff --git a/unittests/gpu_symmetry/README.md b/unittests/gpu_symmetry/README.md
index 6d0e8fe..36f556d 100644
--- a/unittests/gpu_symmetry/README.md
+++ b/unittests/gpu_symmetry/README.md
@@ -73,3 +73,9 @@ For two cell-bearing output atoms the relative cell is `(Rb - Ra) % period`;
 with `Aa01` at the origin it is `Rb % period`. An empty sector retains the existing
 unrestricted behavior. The numerical regression also disables symmetry after each
 filtered run to verify replacement of the previous backend's filter.
+
+Backend entry points reuse an existing filter when its exact type matches the
+symmetry flag/backend and its configuration references belong to the current LRI.
+Both CPU modes share the CPU filter type. Sector/period content changes are read
+through those references; a copied LRI rebinds its filter on first use. The
+regression verifies pointer reuse, in-place sector updates and copy isolation.
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index c4541c8..148a978 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -147,6 +147,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
             name = RI::Label_Tools::get_name(all_labels[index]);
         }
         lri.flag_symmetry = false;
+        lri.filter_atom.reset(); // A missing filter must be initialized lazily.
         lri.cal_mode = RI::LRI_Cal_Mode::CPU;
         Tensors<T> full;
         lri.cal_loop3(labels, full, 1.0);
@@ -181,6 +182,20 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 Tensors<T> actual;
                 lri.cal_loop3(labels, actual, 1.0);
                 compare(actual, expected, context);
+                const auto active_filter = lri.filter_atom;
+                lri.irreducible_sector.clear();
+                Tensors<T> changed_sector;
+                lri.cal_loop3(labels, changed_sector, 1.0);
+                compare(changed_sector, full, context + " sector changed in place");
+                if (lri.filter_atom != active_filter)
+                    throw std::runtime_error(context + ": sector update rebuilt filter");
+                lri.irreducible_sector = sector;
+                auto copied = lri;
+                Tensors<T> copied_result;
+                copied.cal_loop3(labels, copied_result, 1.0);
+                compare(copied_result, expected, context + " copied LRI");
+                if (copied.filter_atom == active_filter)
+                    throw std::runtime_error(context + ": copied LRI retained source filter references");
 
                 // Seed every block, including those outside the selected sector.
                 // Filtering must preserve them while signed/zero increments work.
@@ -194,6 +209,8 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 for (const double factor : {-0.5, 0.0, 1.25})
                 {
                     lri.cal_loop3(labels, actual, factor);
+                    if (lri.filter_atom != active_filter)
+                        throw std::runtime_error(context + ": repeated contraction rebuilt filter");
                     sum += factor;
                     Tensors<T> accumulated = copy_tensors(seed);
                     for (const auto& left : expected)
@@ -212,6 +229,12 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 Tensors<T> unrestricted;
                 lri.cal_loop3(labels, unrestricted, 1.0);
                 compare(unrestricted, full, context + " symmetry disabled");
+                const auto disabled_filter = lri.filter_atom;
+                unrestricted.clear();
+                lri.cal_loop3(labels, unrestricted, 1.0);
+                compare(unrestricted, full, context + " symmetry disabled twice");
+                if (lri.filter_atom != disabled_filter)
+                    throw std::runtime_error(context + ": disabled filter was rebuilt");
                 lri.flag_symmetry = true;
             }
         }

From 4dfade3b8ce8d1d9134599faadfc0f01cc700f4b Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Tue, 10 Nov 2026 04:02:31 -0500
Subject: [PATCH 3/5] Revert "Reuse backend symmetry filters bound to the
 current LRI"

This reverts commit f6273c96598e4eb32c7d1b4f16d3a32c4b9e8135.
---
 include/RI/ri/LRI-cal_loop3_CPU.hpp           |  6 ++++-
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   |  6 ++++-
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |  6 ++++-
 include/RI/ri/LRI.h                           |  4 ----
 include/RI/ri/LRI.hpp                         | 23 -------------------
 unittests/gpu_symmetry/README.md              |  6 -----
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    | 23 -------------------
 7 files changed, 15 insertions(+), 59 deletions(-)

diff --git a/include/RI/ri/LRI-cal_loop3_CPU.hpp b/include/RI/ri/LRI-cal_loop3_CPU.hpp
index 6fed27d..f618e25 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -25,7 +25,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	this->template prepare_filter_atom<Filter_Atom_Symmetry<TA,TC,Tdata>>();
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
 	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 46e80d7..308fa87 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
@@ -25,7 +25,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	this->template prepare_filter_atom<Filter_Atom_Symmetry<TA,TC,Tdata>>();
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 0d7433f..24090c2 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -36,7 +36,11 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	this->template prepare_filter_atom<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>();
+	if (this->flag_symmetry)
+		this->filter_atom = std::make_shared<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>(
+			this->period, this->irreducible_sector);
+	else
+		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
 	using namespace Array_Operator;
 
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index df07fe9..e1c154d 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -101,10 +101,6 @@ class LRI
 	std::map<std::pair<TA,TA>, std::set<TC>> irreducible_sector;
 	std::shared_ptr<Filter_Atom<TA,TAC>> filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
-private:
-	template<typename Tfilter>
-	void prepare_filter_atom();
-
 public:		// private:
 	TC period;
 	MPI_Comm mpi_comm;
diff --git a/include/RI/ri/LRI.hpp b/include/RI/ri/LRI.hpp
index eddfbf7..127e656 100644
--- a/include/RI/ri/LRI.hpp
+++ b/include/RI/ri/LRI.hpp
@@ -8,7 +8,6 @@
 #include "LRI.h"
 #include "../ri/Label.h"
 #include <limits>
-#include <typeinfo>
 
 namespace RI
 {
@@ -29,28 +28,6 @@ LRI<TA,Tcell,Ndim,Tdata>::LRI()
 		this->period[i] = std::numeric_limits<Tcell>::max()/4;		// /4 for not out of range when Array_Operator::operator%
 }
 
-template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
-template<typename Tfilter>
-void LRI<TA,Tcell,Ndim,Tdata>::prepare_filter_atom()
-{
-	if (!this->flag_symmetry)
-	{
-		if (!this->filter_atom || typeid(*this->filter_atom) != typeid(Filter_Atom<TA,TAC>))
-			this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
-		return;
-	}
-
-	if (this->filter_atom && typeid(*this->filter_atom) == typeid(Tfilter))
-	{
-		const auto &filter = static_cast<const Tfilter &>(*this->filter_atom);
-		// A copied LRI must bind the filter to its own configuration.
-		if (&filter.symmetry.period == &this->period
-			&& &filter.symmetry.irreducible_sector_ == &this->irreducible_sector)
-			return;
-	}
-	this->filter_atom = std::make_shared<Tfilter>(this->period, this->irreducible_sector);
-}
-
 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,
diff --git a/unittests/gpu_symmetry/README.md b/unittests/gpu_symmetry/README.md
index 36f556d..6d0e8fe 100644
--- a/unittests/gpu_symmetry/README.md
+++ b/unittests/gpu_symmetry/README.md
@@ -73,9 +73,3 @@ For two cell-bearing output atoms the relative cell is `(Rb - Ra) % period`;
 with `Aa01` at the origin it is `Rb % period`. An empty sector retains the existing
 unrestricted behavior. The numerical regression also disables symmetry after each
 filtered run to verify replacement of the previous backend's filter.
-
-Backend entry points reuse an existing filter when its exact type matches the
-symmetry flag/backend and its configuration references belong to the current LRI.
-Both CPU modes share the CPU filter type. Sector/period content changes are read
-through those references; a copied LRI rebinds its filter on first use. The
-regression verifies pointer reuse, in-place sector updates and copy isolation.
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index 148a978..c4541c8 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -147,7 +147,6 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
             name = RI::Label_Tools::get_name(all_labels[index]);
         }
         lri.flag_symmetry = false;
-        lri.filter_atom.reset(); // A missing filter must be initialized lazily.
         lri.cal_mode = RI::LRI_Cal_Mode::CPU;
         Tensors<T> full;
         lri.cal_loop3(labels, full, 1.0);
@@ -182,20 +181,6 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 Tensors<T> actual;
                 lri.cal_loop3(labels, actual, 1.0);
                 compare(actual, expected, context);
-                const auto active_filter = lri.filter_atom;
-                lri.irreducible_sector.clear();
-                Tensors<T> changed_sector;
-                lri.cal_loop3(labels, changed_sector, 1.0);
-                compare(changed_sector, full, context + " sector changed in place");
-                if (lri.filter_atom != active_filter)
-                    throw std::runtime_error(context + ": sector update rebuilt filter");
-                lri.irreducible_sector = sector;
-                auto copied = lri;
-                Tensors<T> copied_result;
-                copied.cal_loop3(labels, copied_result, 1.0);
-                compare(copied_result, expected, context + " copied LRI");
-                if (copied.filter_atom == active_filter)
-                    throw std::runtime_error(context + ": copied LRI retained source filter references");
 
                 // Seed every block, including those outside the selected sector.
                 // Filtering must preserve them while signed/zero increments work.
@@ -209,8 +194,6 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 for (const double factor : {-0.5, 0.0, 1.25})
                 {
                     lri.cal_loop3(labels, actual, factor);
-                    if (lri.filter_atom != active_filter)
-                        throw std::runtime_error(context + ": repeated contraction rebuilt filter");
                     sum += factor;
                     Tensors<T> accumulated = copy_tensors(seed);
                     for (const auto& left : expected)
@@ -229,12 +212,6 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                 Tensors<T> unrestricted;
                 lri.cal_loop3(labels, unrestricted, 1.0);
                 compare(unrestricted, full, context + " symmetry disabled");
-                const auto disabled_filter = lri.filter_atom;
-                unrestricted.clear();
-                lri.cal_loop3(labels, unrestricted, 1.0);
-                compare(unrestricted, full, context + " symmetry disabled twice");
-                if (lri.filter_atom != disabled_filter)
-                    throw std::runtime_error(context + ": disabled filter was rebuilt");
                 lri.flag_symmetry = true;
             }
         }

From d8346c63f05e19176dde3f2b2ad188cfa1685ba9 Mon Sep 17 00:00:00 2001
From: linpz <linpz@mail.ustc.edu.cn>
Date: Tue, 10 Nov 2026 14:37:27 +0800
Subject: [PATCH 4/5] 1. add Filter_Atom_Symmetry_Wrapper

---
 include/RI/physics/Exx.hpp                    |   8 +-
 include/RI/physics/GW.hpp                     |   8 +-
 include/RI/physics/RPA.hpp                    |   8 +-
 .../symmetry/Filter_Atom_Symmetry_Wrapper.h   | 101 ++++++++++++++++++
 include/RI/ri/LRI-cal_loop3_CPU.hpp           |   7 --
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   |   7 --
 include/RI/ri/LRI-cal_loop3_GPU.hpp           |   7 --
 include/RI/ri/LRI.h                           |   3 -
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    |  13 +--
 9 files changed, 121 insertions(+), 41 deletions(-)
 create mode 100644 include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h

diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp
index 42ff97d..6c55afd 100644
--- a/include/RI/physics/Exx.hpp
+++ b/include/RI/physics/Exx.hpp
@@ -9,6 +9,7 @@
 #include "../ri/Cell_Nearest.h"
 #include "../ri/Label.h"
 #include "../global/Map_Operator.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 #include <cassert>
 
@@ -40,8 +41,11 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	this->lri.flag_symmetry = flag_symmetry;
-	this->lri.irreducible_sector = irreducible_sector;
+	if(flag_symmetry)
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
+	else
+		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/GW.hpp b/include/RI/physics/GW.hpp
index 8733ee7..9746ecf 100644
--- a/include/RI/physics/GW.hpp
+++ b/include/RI/physics/GW.hpp
@@ -6,6 +6,7 @@
 
 #include "GW.h"
 #include "../ri/Label.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 namespace RI
 {
@@ -35,8 +36,11 @@ void GW<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	this->lri.flag_symmetry = flag_symmetry;
-	this->lri.irreducible_sector = irreducible_sector;
+	if(flag_symmetry)
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
+	else
+		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/RPA.hpp b/include/RI/physics/RPA.hpp
index 454e5cf..51a8aed 100644
--- a/include/RI/physics/RPA.hpp
+++ b/include/RI/physics/RPA.hpp
@@ -7,6 +7,7 @@
 
 #include "RPA.h"
 #include "../ri/Label.h"
+#include "./symmetry/Filter_Atom_Symmetry_Wrapper.h"
 
 #include <cassert>
 
@@ -36,8 +37,11 @@ void RPA<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const bool flag_symmetry,
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
-	this->lri.flag_symmetry = flag_symmetry;
-	this->lri.irreducible_sector = irreducible_sector;
+	if(flag_symmetry)
+		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
+			this->period, irreducible_sector, this->lri.cal_mode);
+	else
+		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h
new file mode 100644
index 0000000..098176a
--- /dev/null
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h
@@ -0,0 +1,101 @@
+// ===================
+//  Author: Peize Lin
+//  date: 2026.11.10
+// ===================
+
+#pragma once
+
+#ifndef __GPU_RI
+
+#include "Filter_Atom_Symmetry.h"
+
+namespace RI
+{
+
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_Wrapper: public Filter_Atom_Symmetry<TA,TC,Tdata>
+{
+  public:
+	Filter_Atom_Symmetry_Wrapper(
+		const TC &period,
+		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
+		const LRI_Cal_Mode &cal_mode_in)
+		: Filter_Atom_Symmetry<TA,TC,Tdata>(period, irsec){}
+};
+
+}
+
+#else
+
+#include "../../ri/Filter_Atom.h"
+#include "Filter_Atom_Symmetry.h"
+#include "Filter_Atom_Symmetry_GPU.h"
+#include "../../ri/Label_Tools.h"
+
+#include <stdexcept>
+
+namespace RI
+{
+
+template<typename TA, typename TC, typename Tdata>
+class Filter_Atom_Symmetry_Wrapper: public Filter_Atom<TA, std::pair<TA, TC>>
+{
+  public:
+	using TAC = std::pair<TA, TC>;
+
+	Filter_Atom_Symmetry_Wrapper(
+		const TC &period,
+		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
+		const LRI_Cal_Mode &cal_mode_in)
+		:cal_mode(cal_mode_in),
+		 filter_CPU(period, irsec),
+		 filter_GPU(period, irsec){}
+
+	virtual bool filter_for1(const Label::ab_ab &label, const TA &A1) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for1(label, A1)
+			: this->filter_CPU.filter_for1(label, A1);
+	}
+	virtual bool filter_for1(const Label::ab_ab &label, const TAC &A1) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for1(label, A1)
+			: this->filter_CPU.filter_for1(label, A1);
+	}
+
+	virtual bool filter_for2(const Label::ab_ab &label, const TA &A1, const TAC &A2) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for2(label, A1, A2)
+			: this->filter_CPU.filter_for2(label, A1, A2);
+	}
+
+	virtual bool filter_for32(const Label::ab_ab &label, const TA &A1, const TAC &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? this->filter_GPU.filter_for32(label, A1, A2, A3)
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TA &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? false
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TAC &A2, const TAC &A3) const override
+	{
+		return this->cal_mode==LRI_Cal_Mode::GPU
+			? false
+			: this->filter_CPU.filter_for32(label, A1, A2, A3);
+	}
+
+  public:	// private
+	const LRI_Cal_Mode &cal_mode;
+	Filter_Atom_Symmetry<TA,TC,Tdata> filter_CPU;
+	Filter_Atom_Symmetry_GPU<TA,TC,Tdata> filter_GPU;
+};
+
+}
+
+#endif
\ 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 f618e25..e702ebb 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -6,7 +6,6 @@
 #pragma once
 
 #include "LRI.h"
-#include "../physics/symmetry/Filter_Atom_Symmetry.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Tensor_Multiply.h"
@@ -25,12 +24,6 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
-
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
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 308fa87..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
@@ -6,7 +6,6 @@
 #pragma once
 
 #include "LRI.h"
-#include "../physics/symmetry/Filter_Atom_Symmetry.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Tensor_Multiply.h"
@@ -25,12 +24,6 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
-
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 24090c2..0c602d7 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -8,7 +8,6 @@
 #ifdef __GPU_RI
 
 #include "LRI.h"
-#include "../physics/symmetry/Filter_Atom_Symmetry_GPU.h"
 #include "LRI_Cal_Aux.h"
 #include "../global/Array_Operator.h"
 #include "../global/Map_Operator.h"
@@ -36,12 +35,6 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
-	if (this->flag_symmetry)
-		this->filter_atom = std::make_shared<Filter_Atom_Symmetry_GPU<TA,TC,Tdata>>(
-			this->period, this->irreducible_sector);
-	else
-		this->filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
-
 	using namespace Array_Operator;
 
 	constexpr std::size_t memory_limit = std::size_t(5)*1024*1024*1024;
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index e1c154d..c3db89d 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -96,9 +96,6 @@ class LRI
 		parallel = std::make_shared<Parallel_LRI_Equally<TA,Tcell,Ndim,Tdata>>();
 	std::unordered_map< Label::ab, RI_Tools::T_filter_func<Tdata> >
 		filter_funcs;
-	// Construct the symmetry filter for the selected backend at contraction time.
-	bool flag_symmetry = false;
-	std::map<std::pair<TA,TA>, std::set<TC>> irreducible_sector;
 	std::shared_ptr<Filter_Atom<TA,TAC>> filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
 
 public:		// private:
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index c4541c8..233bc34 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -146,7 +146,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
             labels.push_back(all_labels[index]);
             name = RI::Label_Tools::get_name(all_labels[index]);
         }
-        lri.flag_symmetry = false;
+        lri.filter_atom = std::make_shared<RI::Filter_Atom<int, AtomCell>>();
         lri.cal_mode = RI::LRI_Cal_Mode::CPU;
         Tensors<T> full;
         lri.cal_loop3(labels, full, 1.0);
@@ -169,8 +169,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                     if (keep)
                         expected[left.first][right.first] = right.second.copy();
                 }
-            lri.flag_symmetry = true;
-            lri.irreducible_sector = sector;
+            lri.filter_atom = std::make_shared<RI::Filter_Atom_Symmetry<int, Cell, T>>(period, sector);
             for (const auto mode : modes)
             {
                 lri.cal_mode = mode;
@@ -205,14 +204,6 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                         }
                     compare(actual, accumulated, context + " factor=" + std::to_string(factor));
                 }
-
-                // Disable symmetry after a restricted run on this same backend.
-                // This must replace the previous backend-specific filter.
-                lri.flag_symmetry = false;
-                Tensors<T> unrestricted;
-                lri.cal_loop3(labels, unrestricted, 1.0);
-                compare(unrestricted, full, context + " symmetry disabled");
-                lri.flag_symmetry = true;
             }
         }
         std::cout << "PASS " << scalar << " label=" << name << std::endl;

From ab77f0bb222fd9420674f73d5d0fc02ba9c3ab1f Mon Sep 17 00:00:00 2001
From: maki49 <1579492865@qq.com>
Date: Tue, 10 Nov 2026 04:08:21 -0500
Subject: [PATCH 5/5] Select filter policy once before contraction loops

---
 include/RI/physics/Exx.hpp                    |   6 +-
 include/RI/physics/GW.hpp                     |   6 +-
 include/RI/physics/RPA.hpp                    |   6 +-
 .../symmetry/Filter_Atom_Symmetry_Wrapper.h   | 127 ++++------
 include/RI/ri/LRI-cal_loop3_CPU.hpp           | 148 ++++++------
 .../LRI-cal_loop3_CPU_fine_grained_lock.hpp   | 148 ++++++------
 include/RI/ri/LRI-cal_loop3_GPU.hpp           | 228 +++++++++---------
 include/RI/ri/LRI.h                           |   3 +-
 unittests/gpu_symmetry/README.md              |  23 +-
 .../gpu_symmetry/test_lri_gpu_symmetry.cpp    |  33 ++-
 10 files changed, 375 insertions(+), 353 deletions(-)

diff --git a/include/RI/physics/Exx.hpp b/include/RI/physics/Exx.hpp
index 6c55afd..ec837aa 100644
--- a/include/RI/physics/Exx.hpp
+++ b/include/RI/physics/Exx.hpp
@@ -42,10 +42,10 @@ void Exx<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
-			this->period, irreducible_sector, this->lri.cal_mode);
+		this->lri.filter_atom = Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>(
+			this->period, irreducible_sector);
 	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+		this->lri.filter_atom = {};
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/GW.hpp b/include/RI/physics/GW.hpp
index 9746ecf..ad982fb 100644
--- a/include/RI/physics/GW.hpp
+++ b/include/RI/physics/GW.hpp
@@ -37,10 +37,10 @@ void GW<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
-			this->period, irreducible_sector, this->lri.cal_mode);
+		this->lri.filter_atom = Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>(
+			this->period, irreducible_sector);
 	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+		this->lri.filter_atom = {};
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/RPA.hpp b/include/RI/physics/RPA.hpp
index 51a8aed..5d93190 100644
--- a/include/RI/physics/RPA.hpp
+++ b/include/RI/physics/RPA.hpp
@@ -38,10 +38,10 @@ void RPA<TA,Tcell,Ndim,Tdata>::set_symmetry(
 	const std::map<std::pair<TA,TA>, std::set<TC>> &irreducible_sector)
 {
 	if(flag_symmetry)
-		this->lri.filter_atom = std::make_shared<Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>>(
-			this->period, irreducible_sector, this->lri.cal_mode);
+		this->lri.filter_atom = Filter_Atom_Symmetry_Wrapper<TA,TC,Tdata>(
+			this->period, irreducible_sector);
 	else
-		this->lri.filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+		this->lri.filter_atom = {};
 }
 
 template<typename TA, typename Tcell, std::size_t Ndim, typename Tdata>
diff --git a/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h b/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h
index 098176a..23a538d 100644
--- a/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h
+++ b/include/RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h
@@ -5,97 +5,60 @@
 
 #pragma once
 
-#ifndef __GPU_RI
-
-#include "Filter_Atom_Symmetry.h"
-
-namespace RI
-{
-
-template<typename TA, typename TC, typename Tdata>
-class Filter_Atom_Symmetry_Wrapper: public Filter_Atom_Symmetry<TA,TC,Tdata>
-{
-  public:
-	Filter_Atom_Symmetry_Wrapper(
-		const TC &period,
-		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
-		const LRI_Cal_Mode &cal_mode_in)
-		: Filter_Atom_Symmetry<TA,TC,Tdata>(period, irsec){}
-};
-
-}
-
-#else
-
-#include "../../ri/Filter_Atom.h"
 #include "Filter_Atom_Symmetry.h"
+#include "../../ri/LRI.h"
+#ifdef __GPU_RI
 #include "Filter_Atom_Symmetry_GPU.h"
-#include "../../ri/Label_Tools.h"
-
-#include <stdexcept>
+#endif
+#include <memory>
 
 namespace RI
 {
 
+// Own both policies and select once before entering the contraction loops.
+// This callable is not an atom filter and does not forward filter callbacks.
 template<typename TA, typename TC, typename Tdata>
-class Filter_Atom_Symmetry_Wrapper: public Filter_Atom<TA, std::pair<TA, TC>>
+class Filter_Atom_Symmetry_Wrapper
 {
-  public:
-	using TAC = std::pair<TA, TC>;
-
-	Filter_Atom_Symmetry_Wrapper(
-		const TC &period,
-		const std::map<std::pair<TA,TA>, std::set<TC>> &irsec,
-		const LRI_Cal_Mode &cal_mode_in)
-		:cal_mode(cal_mode_in),
-		 filter_CPU(period, irsec),
-		 filter_GPU(period, irsec){}
+    using TAC = std::pair<TA, TC>;
+    struct Configuration
+    {
+        TC period;
+        std::map<std::pair<TA,TA>, std::set<TC>> sector;
+    };
+    // Shared immutable configuration keeps filter references valid after copies.
+    std::shared_ptr<const Configuration> configuration;
+    Filter_Atom_Symmetry<TA,TC,Tdata> filter_CPU;
+#ifdef __GPU_RI
+    Filter_Atom_Symmetry_GPU<TA,TC,Tdata> filter_GPU;
+#endif
 
-	virtual bool filter_for1(const Label::ab_ab &label, const TA &A1) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? this->filter_GPU.filter_for1(label, A1)
-			: this->filter_CPU.filter_for1(label, A1);
-	}
-	virtual bool filter_for1(const Label::ab_ab &label, const TAC &A1) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? this->filter_GPU.filter_for1(label, A1)
-			: this->filter_CPU.filter_for1(label, A1);
-	}
-
-	virtual bool filter_for2(const Label::ab_ab &label, const TA &A1, const TAC &A2) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? this->filter_GPU.filter_for2(label, A1, A2)
-			: this->filter_CPU.filter_for2(label, A1, A2);
-	}
-
-	virtual bool filter_for32(const Label::ab_ab &label, const TA &A1, const TAC &A2, const TAC &A3) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? this->filter_GPU.filter_for32(label, A1, A2, A3)
-			: this->filter_CPU.filter_for32(label, A1, A2, A3);
-	}
-	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TA &A2, const TAC &A3) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? false
-			: this->filter_CPU.filter_for32(label, A1, A2, A3);
-	}
-	virtual bool filter_for32(const Label::ab_ab &label, const TAC &A1, const TAC &A2, const TAC &A3) const override
-	{
-		return this->cal_mode==LRI_Cal_Mode::GPU
-			? false
-			: this->filter_CPU.filter_for32(label, A1, A2, A3);
-	}
-
-  public:	// private
-	const LRI_Cal_Mode &cal_mode;
-	Filter_Atom_Symmetry<TA,TC,Tdata> filter_CPU;
-	Filter_Atom_Symmetry_GPU<TA,TC,Tdata> filter_GPU;
+  public:
+    Filter_Atom_Symmetry_Wrapper(
+        const TC &period,
+        const std::map<std::pair<TA,TA>, std::set<TC>> &irsec)
+        : configuration(std::make_shared<Configuration>(Configuration{period, irsec})),
+          filter_CPU(configuration->period, configuration->sector)
+#ifdef __GPU_RI
+        , filter_GPU(configuration->period, configuration->sector)
+#endif
+    {}
+
+    const Filter_Atom<TA,TAC>* operator()(LRI_Cal_Mode mode) const
+    {
+        switch (mode)
+        {
+            case LRI_Cal_Mode::CPU:
+            case LRI_Cal_Mode::CPU_fine_grained_lock:
+                return &filter_CPU;
+#ifdef __GPU_RI
+            case LRI_Cal_Mode::GPU:
+                return &filter_GPU;
+#endif
+            default:
+                throw std::invalid_argument("Unknown contraction mode for symmetry filter");
+        }
+    }
 };
 
 }
-
-#endif
\ 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 e702ebb..ba60ef2 100644
--- a/include/RI/ri/LRI-cal_loop3_CPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_CPU.hpp
@@ -24,6 +24,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	const Filter_Atom<TA,TAC> unfiltered;
+	const Filter_Atom<TA,TAC> *filter_atom = this->filter_atom
+		? this->filter_atom(LRI_Cal_Mode::CPU) : &unfiltered;
+
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
@@ -71,19 +75,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Aa2 : list_Aa2)
 					{
-						if(this->filter_atom->filter_for1(label,Aa2))	continue;
+						if(filter_atom->filter_for1(label,Aa2))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa2,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa2,Ab01))	continue;
 							// D_mul = D_a * D_a0b0 * D_a1b1
 							Tensor<Tdata> D_mul;
 							for(const TA &Aa01 : list_Aa01)
 							{
-								if(this->filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
+								if(filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
@@ -102,7 +106,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								// a2b2 = a2b0b1 * b0b1b2
@@ -135,19 +139,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Aa2 : list_Aa2)
 					{
-						if(this->filter_atom->filter_for1(label,Aa2))	continue;
+						if(filter_atom->filter_for1(label,Aa2))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa2,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa2,Ab01))	continue;
 							// D_mul = D_a * D_a0b1 * D_a1b0
 							Tensor<Tdata> D_mul;
 							for(const TA &Aa01 : list_Aa01)
 							{
-								if(this->filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
+								if(filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
@@ -166,7 +170,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								// a2b2 = a2b0b1 * b0b1b2
@@ -201,21 +205,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 							if(D_a0b0.empty())	continue;
 							// D_mul = D_b * D_a1b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -230,7 +234,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a * D_a0b0
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// b1a1a0 = b0b1a1 * a0b0
@@ -265,21 +269,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 							if(D_a0b1.empty())	continue;
 							// D_mul = D_b * D_a1b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -294,7 +298,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a * D_a0b1
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// a0a1b0 = a0b1 * a1b0b1
@@ -329,21 +333,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 							if(D_a1b0.empty())	continue;
 							// D_mul = D_b * D_a0b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -358,7 +362,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a * D_a1b0
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// b1a0a1 = b0b1a0 * a1b0
@@ -393,21 +397,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 							if(D_a1b1.empty())	continue;
 							// D_mul = D_b * D_a0b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -422,7 +426,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a * D_a1b1
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// a1a0b0 = a1b1 * a0b0b1
@@ -459,21 +463,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 							if(D_a0b0.empty())	continue;
 							// D_mul = D_a * D_a2b1
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -488,7 +492,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -524,21 +528,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 							if(D_a0b1.empty())	continue;
 							// D_mul = D_a * D_a2b0
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -553,7 +557,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a0b1 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -589,21 +593,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 							if(D_a1b0.empty())	continue;
 							// D_mul = D_a * D_a2b1
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -618,7 +622,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a1b0 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -654,21 +658,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 							if(D_a1b1.empty())	continue;
 							// D_mul = D_a * D_a2b0
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -683,7 +687,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a1b1 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -721,19 +725,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -748,7 +752,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
@@ -786,19 +790,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -813,7 +817,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a0b1 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
@@ -851,19 +855,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -878,7 +882,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a1b0 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
@@ -916,19 +920,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -943,7 +947,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							// D_result = D_mul * D_a1b1 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
@@ -983,17 +987,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a1b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -1009,7 +1013,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -1048,17 +1052,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a0b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -1074,7 +1078,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -1113,17 +1117,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a0b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -1139,7 +1143,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -1178,17 +1182,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a1b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -1204,7 +1208,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
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 ac08b09..8407fec 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
@@ -24,6 +24,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	const Filter_Atom<TA,TAC> unfiltered;
+	const Filter_Atom<TA,TAC> *filter_atom = this->filter_atom
+		? this->filter_atom(LRI_Cal_Mode::CPU_fine_grained_lock) : &unfiltered;
+
 	using namespace Array_Operator;
 
 	const Data_Pack_Wrapper<TA,TC,Tdata> data_wrapper(this->data_pool, this->data_ab_name);
@@ -71,19 +75,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Aa2 : list_Aa2)
 					{
-						if(this->filter_atom->filter_for1(label,Aa2))	continue;
+						if(filter_atom->filter_for1(label,Aa2))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa2,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa2,Ab01))	continue;
 							// D_mul = D_a * D_a0b0 * D_a1b1
 							Tensor<Tdata> D_mul;
 							for(const TA &Aa01 : list_Aa01)
 							{
-								if(this->filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
+								if(filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
@@ -102,7 +106,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								// a2b2 = a2b0b1 * b0b1b2
@@ -135,19 +139,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Aa2 : list_Aa2)
 					{
-						if(this->filter_atom->filter_for1(label,Aa2))	continue;
+						if(filter_atom->filter_for1(label,Aa2))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa2,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa2,Ab01))	continue;
 							// D_mul = D_a * D_a0b1 * D_a1b0
 							Tensor<Tdata> D_mul;
 							for(const TA &Aa01 : list_Aa01)
 							{
-								if(this->filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
+								if(filter_atom->filter_for31(label,Aa2,Ab01,Aa01))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
@@ -166,7 +170,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa2,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								// a2b2 = a2b0b1 * b0b1b2
@@ -201,21 +205,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 							if(D_a0b0.empty())	continue;
 							// D_mul = D_b * D_a1b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -230,7 +234,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a * D_a0b0
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// b1a1a0 = b0b1a1 * a0b0
@@ -265,21 +269,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 							if(D_a0b1.empty())	continue;
 							// D_mul = D_b * D_a1b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -294,7 +298,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a * D_a0b1
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// a0a1b0 = a0b1 * a1b0b1
@@ -329,21 +333,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 							if(D_a1b0.empty())	continue;
 							// D_mul = D_b * D_a0b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -358,7 +362,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a * D_a1b0
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// b1a0a1 = b0b1a0 * a1b0
@@ -393,21 +397,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TAC &Ab01 : list_Ab01)
 					{
-						if(this->filter_atom->filter_for1(label,Ab01))	continue;
+						if(filter_atom->filter_for1(label,Ab01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for2(label,Ab01,Aa01))	continue;
+							if(filter_atom->filter_for2(label,Ab01,Aa01))	continue;
 							const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 							if(D_a1b1.empty())	continue;
 							// D_mul = D_b * D_a0b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Ab01,Aa01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -422,7 +426,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a * D_a1b1
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Ab01,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								// a1a0b0 = a1b1 * a0b0b1
@@ -459,21 +463,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 							if(D_a0b0.empty())	continue;
 							// D_mul = D_a * D_a2b1
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -488,7 +492,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -525,21 +529,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 							if(D_a0b1.empty())	continue;
 							// D_mul = D_a * D_a2b0
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -554,7 +558,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a0b1 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -591,21 +595,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 							if(D_a1b0.empty())	continue;
 							// D_mul = D_a * D_a2b1
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -620,7 +624,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a1b0 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -657,21 +661,21 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 							if(D_a1b1.empty())	continue;
 							// D_mul = D_a * D_a2b0
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -686,7 +690,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a1b1 * D_b
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 
@@ -725,19 +729,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -752,7 +756,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a0b0 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
@@ -791,19 +795,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -818,7 +822,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a0b1 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
@@ -857,19 +861,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -884,7 +888,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a1b0 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
@@ -923,19 +927,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						std::map<TAC,Tensor<Tdata>> Ds_result_fixed;
 
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 						{
 							const TAC &Ab2 = list_Ab2[ib2];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab2))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab2))	continue;
 							// D_mul = D_a * D_a2b2
 							Tensor<Tdata> D_mul;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab2,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
@@ -950,7 +954,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							// D_result = D_mul * D_a1b1 * D_b
 							for(const TAC &Ab01 : list_Ab01)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab2,Ab01))	continue;
+								if(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 Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
@@ -991,17 +995,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a1b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -1017,7 +1021,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -1055,17 +1059,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a0b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -1081,7 +1085,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
@@ -1119,17 +1123,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a0b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
@@ -1145,7 +1149,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
@@ -1183,17 +1187,17 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 
 					for(const TA &Aa01 : list_Aa01)
 					{
-						if(this->filter_atom->filter_for1(label,Aa01))	continue;
+						if(filter_atom->filter_for1(label,Aa01))	continue;
 						#pragma omp for schedule(dynamic) nowait
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for2(label,Aa01,Ab01))	continue;
+							if(filter_atom->filter_for2(label,Aa01,Ab01))	continue;
 							// D_mul1 = D_b * D_a1b2
 							Tensor<Tdata> D_mul1;
 							for(const TAC &Ab2 : list_Ab2)
 							{
-								if(this->filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for31(label,Aa01,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
@@ -1209,7 +1213,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_CPU_fine_grained_lock(
 							Tensor<Tdata> D_mul2;
 							for(const TAC &Aa2 : list_Aa2)
 							{
-								if(this->filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
+								if(filter_atom->filter_for32(label,Aa01,Ab01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
diff --git a/include/RI/ri/LRI-cal_loop3_GPU.hpp b/include/RI/ri/LRI-cal_loop3_GPU.hpp
index 0c602d7..7734160 100644
--- a/include/RI/ri/LRI-cal_loop3_GPU.hpp
+++ b/include/RI/ri/LRI-cal_loop3_GPU.hpp
@@ -35,6 +35,10 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 	std::map<TA, std::map<TAC, Tensor<Tdata>>> &Ds_result,
 	const Tdata fac_add_Ds)
 {
+	const Filter_Atom<TA,TAC> unfiltered;
+	const Filter_Atom<TA,TAC> *filter_atom = this->filter_atom
+		? this->filter_atom(LRI_Cal_Mode::GPU) : &unfiltered;
+
 	using namespace Array_Operator;
 
 	constexpr std::size_t memory_limit = std::size_t(5)*1024*1024*1024;
@@ -92,12 +96,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 								if(D_a0b0.empty())	continue;
 								const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
@@ -105,7 +109,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_for31(label, Aa01, Ab01, Ab2))	continue;
+									if (filter_atom->filter_for31(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -156,18 +160,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if (this->filter_atom->filter_for2(label, Aa01, Aa2))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if (D_a.empty())	continue;
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2) // G
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Aa2, Ab2))	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
 									if (!D_mul.exist)	continue;
 
@@ -228,12 +232,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 								if(D_a0b1.empty())	continue;
 								const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
@@ -241,7 +245,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_for31(label, Aa01, Ab01, Ab2))	continue;
+									if (filter_atom->filter_for31(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -292,18 +296,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if (this->filter_atom->filter_for2(label, Aa01, Aa2))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if (D_a.empty())	continue;
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2) // G
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if (this->filter_atom->filter_for32(label, Aa01, Aa2, Ab2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Aa2, Ab2))	continue;
 									const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab2);
 									if (!D_mul.exist)	continue;
 
@@ -367,18 +371,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -418,12 +422,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 								if(D_a0b0.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -431,7 +435,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -507,18 +511,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -558,12 +562,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 								if(D_a0b1.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -571,7 +575,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -647,18 +651,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -698,12 +702,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 								if(D_a1b0.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -711,7 +715,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -787,18 +791,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for (std::size_t ia01= 0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -838,12 +842,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 								if(D_a1b1.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -851,7 +855,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2) // F
 								{
 									const TAC &Aa2 = list_Aa2[ia2];
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Ab01, Aa2))	continue;
 									const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 									if (D_a.empty())	continue;
 
@@ -929,18 +933,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -980,12 +984,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 								if(D_a0b0.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -993,7 +997,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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -1069,18 +1073,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -1120,12 +1124,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 								if(D_a0b1.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -1133,7 +1137,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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -1209,18 +1213,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -1260,12 +1264,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 								if(D_a1b0.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -1273,7 +1277,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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -1349,18 +1353,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for (std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -1400,12 +1404,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01) // A
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01) // G
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 								if(D_a1b1.empty())	continue;
 								const GPU_Data::Pack &D_mul = rDs_mul.find_2nd(Aa01, Ab01);
@@ -1413,7 +1417,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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
 									const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 									if(D_b.empty())	continue;
 
@@ -1491,18 +1495,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
 									const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
 									if(D_a2b2.empty())	continue;
 
@@ -1541,19 +1545,19 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b0 = tools.get_Ds_ab(Label::ab::a0b0, Aa01, Ab01);
 								if(D_a0b0.empty())	continue;
 								for (std::size_t ib2 = 0; ib2 < list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
 									// The output block is (Aa01, Ab01); Ab2 is only an intermediate atom.
-									if (this->filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	continue;
+									if (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1631,18 +1635,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = Global_Func::find(Ds_a_transpose, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
 									const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
 									if(D_a2b2.empty())	continue;
 
@@ -1681,18 +1685,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a0b1 = tools.get_Ds_ab(Label::ab::a0b1, Aa01, Ab01);
 								if(D_a0b1.empty())	continue;
 								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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1770,18 +1774,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
 									const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
 									if(D_a2b2.empty())	continue;
 
@@ -1820,18 +1824,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b0 = tools.get_Ds_ab(Label::ab::a1b0, Aa01, Ab01);
 								if(D_a1b0.empty())	continue;
 								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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -1909,18 +1913,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ia2 = 0; ia2 < list_Aa2.size(); ++ia2)
 							{
 								const TAC &Aa2 = list_Aa2[ia2];
-								if(this->filter_atom->filter_for2(label,Aa01,Aa2))	continue;
+								if(filter_atom->filter_for2(label,Aa01,Aa2))	continue;
 								const Tensor<Tdata> &D_a = tools.get_Ds_ab(Label::ab::a, Aa01, Aa2);
 								if(D_a.empty())	continue;
 								for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 								{
 									const TAC &Ab2 = list_Ab2[ib2];
-									if(this->filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
+									if(filter_atom->filter_for31(label,Aa01,Aa2,Ab2))	continue;
 									const Tensor<Tdata> &D_a2b2 = tools.get_Ds_ab(Label::ab::a2b2, Aa2, Ab2);
 									if(D_a2b2.empty())	continue;
 
@@ -1959,18 +1963,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const Tensor<Tdata> &D_a1b1 = tools.get_Ds_ab(Label::ab::a1b1, Aa01, Ab01);
 								if(D_a1b1.empty())	continue;
 								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 (filter_atom->filter_for32(label, Aa01, Ab01, Ab2))	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);
@@ -2049,18 +2053,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -2098,7 +2102,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
@@ -2109,7 +2113,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
 									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if(filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -2147,12 +2151,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const GPU_Data::Pack &D_mul_1 = rDs_mul_1.find_2nd(Aa01, Ab01);
 								if (!D_mul_1.exist)	continue;
 								const GPU_Data::Pack &D_mul_2 = rDs_mul_2.find_2nd(Aa01, Ab01);
@@ -2213,18 +2217,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -2262,7 +2266,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
@@ -2273,7 +2277,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
 									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if(filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -2311,12 +2315,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const GPU_Data::Pack &D_mul_1 = rDs_mul_1.find_2nd(Aa01, Ab01);
 								if (!D_mul_1.exist)	continue;
 								const GPU_Data::Pack &D_mul_2 = rDs_mul_2.find_2nd(Aa01, Ab01);
@@ -2377,18 +2381,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a0b2 = tools.get_Ds_ab(Label::ab::a0b2, Aa01, Ab2);
 									if(D_a0b2.empty())	continue;
 
@@ -2426,7 +2430,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
@@ -2437,7 +2441,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
 									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if(filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b1 = tools.get_Ds_ab(Label::ab::a2b1, Aa2, Ab01);
 									if(D_a2b1.empty())	continue;
 
@@ -2475,12 +2479,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const GPU_Data::Pack &D_mul_1 = rDs_mul_1.find_2nd(Aa01, Ab01);
 								if (!D_mul_1.exist)	continue;
 								const GPU_Data::Pack &D_mul_2 = rDs_mul_2.find_2nd(Aa01, Ab01);
@@ -2541,18 +2545,18 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ib01 = 0; ib01 < list_Ab01.size(); ++ib01)
 						{
 							const TAC &Ab01 = list_Ab01[ib01];
-							if(this->filter_atom->filter_for1(label,Ab01))	continue;
+							if(filter_atom->filter_for1(label,Ab01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib2=0; ib2<list_Ab2.size(); ++ib2)
 							{
 								const TAC &Ab2 = list_Ab2[ib2];
-								if(this->filter_atom->filter_for2(label,Ab01,Ab2))	continue;
+								if(filter_atom->filter_for2(label,Ab01,Ab2))	continue;
 								const Tensor<Tdata> &D_b = tools.get_Ds_ab(Label::ab::b, Ab01, Ab2);
 								if(D_b.empty())	continue;
 								for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 								{
 									const TA &Aa01 = list_Aa01[ia01];
-									if(this->filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
+									if(filter_atom->filter_for31(label,Ab01,Ab2,Aa01))	continue;
 									const Tensor<Tdata> &D_a1b2 = tools.get_Ds_ab(Label::ab::a1b2, Aa01, Ab2);
 									if(D_a1b2.empty())	continue;
 
@@ -2590,7 +2594,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for (std::size_t ia01 = 0; ia01 < list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if(this->filter_atom->filter_for1(label,Aa01))	continue;
+							if(filter_atom->filter_for1(label,Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ia2=0; ia2<list_Aa2.size(); ++ia2)
 							{
@@ -2601,7 +2605,7 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 								{
 									const TAC &Ab01 = list_Ab01[ib01];
 									// Screen the output pair; Aa2 is summed over, not an output atom.
-									if(this->filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
+									if(filter_atom->filter_for32(label, Aa01, Aa2, Ab01))	continue;
 									const Tensor<Tdata> &D_a2b0 = tools.get_Ds_ab(Label::ab::a2b0, Aa2, Ab01);
 									if(D_a2b0.empty())	continue;
 
@@ -2639,12 +2643,12 @@ void LRI<TA,Tcell,Ndim,Tdata>::cal_loop3_GPU(
 						for(std::size_t ia01=0; ia01<list_Aa01.size(); ++ia01)
 						{
 							const TA &Aa01 = list_Aa01[ia01];
-							if (this->filter_atom->filter_for1(label, Aa01))	continue;
+							if (filter_atom->filter_for1(label, Aa01))	continue;
 							#pragma omp for schedule(dynamic) nowait
 							for(std::size_t ib01=0; ib01<list_Ab01.size(); ++ib01)
 							{
 								const TAC &Ab01 = list_Ab01[ib01];
-								if (this->filter_atom->filter_for2(label, Aa01, Ab01))	continue;
+								if (filter_atom->filter_for2(label, Aa01, Ab01))	continue;
 								const GPU_Data::Pack &D_mul_1 = rDs_mul_1.find_2nd(Aa01, Ab01);
 								if (!D_mul_1.exist)	continue;
 								const GPU_Data::Pack &D_mul_2 = rDs_mul_2.find_2nd(Aa01, Ab01);
diff --git a/include/RI/ri/LRI.h b/include/RI/ri/LRI.h
index c3db89d..62a3703 100644
--- a/include/RI/ri/LRI.h
+++ b/include/RI/ri/LRI.h
@@ -96,7 +96,8 @@ class LRI
 		parallel = std::make_shared<Parallel_LRI_Equally<TA,Tcell,Ndim,Tdata>>();
 	std::unordered_map< Label::ab, RI_Tools::T_filter_func<Tdata> >
 		filter_funcs;
-	std::shared_ptr<Filter_Atom<TA,TAC>> filter_atom = std::make_shared<Filter_Atom<TA,TAC>>();
+	// Select once per contraction; an empty selector means no filtering.
+	std::function<const Filter_Atom<TA,TAC>*(LRI_Cal_Mode)> filter_atom;
 
 public:		// private:
 	TC period;
diff --git a/unittests/gpu_symmetry/README.md b/unittests/gpu_symmetry/README.md
index 6d0e8fe..e93d384 100644
--- a/unittests/gpu_symmetry/README.md
+++ b/unittests/gpu_symmetry/README.md
@@ -71,5 +71,24 @@ pair. In that family's second stage, the former inner
 
 For two cell-bearing output atoms the relative cell is `(Rb - Ra) % period`;
 with `Aa01` at the origin it is `Rb % period`. An empty sector retains the existing
-unrestricted behavior. The numerical regression also disables symmetry after each
-filtered run to verify replacement of the previous backend's filter.
+unrestricted behavior. The numerical regression also disables symmetry on a copied selector after each
+filtered run to verify independent configuration of copied LRI objects.
+
+## Select filters before contraction loops
+
+`LRI::filter_atom` is an optional callable taking `LRI_Cal_Mode` and returning
+a `const Filter_Atom*`. Each backend invokes it once before its OpenMP/atom loops;
+all inner callbacks go directly to that pointer. An empty callable disables
+filtering. The returned pointer must remain valid for the contraction. Custom
+policies can use the same interface without changes to LRI.
+
+`Filter_Atom_Symmetry_Wrapper` is a callable rather than a filter subclass. It
+constructs CPU and (in GPU builds) GPU policies together and selects the existing
+policy on each contraction. Its shared immutable configuration owns copies of
+period and sector, keeping filter references valid after copies or temporary
+caller data. Call `set_symmetry` again to change that configuration.
+
+The regression asserts one selector invocation per contraction and stable filter
+identity across repeated calls, checks backend switching, copied selectors,
+caller configuration lifetime, and disabling symmetry. All 18 label expectations
+continue to use the independent output mask described above.
diff --git a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
index 233bc34..cb38bfd 100644
--- a/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
+++ b/unittests/gpu_symmetry/test_lri_gpu_symmetry.cpp
@@ -10,7 +10,7 @@
 #include <string>
 #include <vector>
 
-#include <RI/physics/symmetry/Filter_Atom_Symmetry.h>
+#include <RI/physics/symmetry/Filter_Atom_Symmetry_Wrapper.h>
 #include <RI/ri/LRI.h>
 
 namespace
@@ -146,7 +146,7 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
             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.filter_atom = {};
         lri.cal_mode = RI::LRI_Cal_Mode::CPU;
         Tensors<T> full;
         lri.cal_loop3(labels, full, 1.0);
@@ -169,7 +169,19 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                     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);
+            // The selector must own configuration even when caller temporaries change.
+            auto caller_period = period;
+            auto caller_sector = sector;
+            lri.filter_atom = RI::Filter_Atom_Symmetry_Wrapper<int, Cell, T>(caller_period, caller_sector);
+            caller_period[0] = 1;
+            caller_sector.clear();
+            const auto policy = lri.filter_atom;
+            auto selections = std::make_shared<std::size_t>(0);
+            lri.filter_atom = [policy, selections](RI::LRI_Cal_Mode mode) -> const RI::Filter_Atom<int, AtomCell>*
+            {
+                ++*selections;
+                return policy(mode);
+            };
             for (const auto mode : modes)
             {
                 lri.cal_mode = mode;
@@ -177,8 +189,12 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                                             + std::to_string(sector_index) + " mode="
                                             + std::to_string(static_cast<int>(mode));
                 std::cout << "CHECK " << context << std::endl;
+                const auto* selected_filter = lri.filter_atom(mode);
                 Tensors<T> actual;
+                const auto before = *selections;
                 lri.cal_loop3(labels, actual, 1.0);
+                if (*selections != before + 1)
+                    throw std::runtime_error(context + ": selector must run once per contraction");
                 compare(actual, expected, context);
 
                 // Seed every block, including those outside the selected sector.
@@ -192,7 +208,10 @@ 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})
                 {
+                    const auto before_factor = *selections;
                     lri.cal_loop3(labels, actual, factor);
+                    if (*selections != before_factor + 1 || lri.filter_atom(mode) != selected_filter)
+                        throw std::runtime_error(context + ": selection count or filter identity changed");
                     sum += factor;
                     Tensors<T> accumulated = copy_tensors(seed);
                     for (const auto& left : expected)
@@ -204,6 +223,14 @@ void check(const std::string& scalar, T imaginary, int selected, const std::stri
                         }
                     compare(actual, accumulated, context + " factor=" + std::to_string(factor));
                 }
+                auto copied = lri;
+                Tensors<T> copied_result;
+                copied.cal_loop3(labels, copied_result, 1.0);
+                compare(copied_result, expected, context + " copied selector");
+                copied.filter_atom = {};
+                copied_result.clear();
+                copied.cal_loop3(labels, copied_result, 1.0);
+                compare(copied_result, full, context + " symmetry disabled");
             }
         }
         std::cout << "PASS " << scalar << " label=" << name << std::endl;