diff --git a/source/module_hamilt_general/module_xc/exx_info.h b/source/module_hamilt_general/module_xc/exx_info.h
index 24acbdb9b2..9aca0a4088 100644
--- a/source/module_hamilt_general/module_xc/exx_info.h
+++ b/source/module_hamilt_general/module_xc/exx_info.h
@@ -15,7 +15,10 @@ struct Exx_Info
     {
         bool cal_exx = false;
 
-        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, 
+            std::pair<bool, 
+                std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, 
+                    std::vector<std::map<std::string,std::string>>>>> coulomb_settings;
 		// Fock:
 		//		"alpha":		"0"
 		//		"Rcut_type":	"limits" / "spencer"
@@ -52,7 +55,10 @@ struct Exx_Info
 
     struct Exx_Info_RI
     {
-        const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> &coulomb_param;
+        const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, 
+            std::pair<bool, 
+                std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, 
+                    std::vector<std::map<std::string,std::string>>>>> &coulomb_settings;
 
         bool real_number = false;
 
@@ -74,7 +80,7 @@ struct Exx_Info
         int abfs_Lmax = 0; // tmp
 
         Exx_Info_RI(const Exx_Info::Exx_Info_Global& info_global)
-            : coulomb_param(info_global.coulomb_param)
+            : coulomb_settings(info_global.coulomb_settings)
         {
         }
     };
diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
index 32c4b30885..bc7fcca691 100644
--- a/source/module_io/input_conv.cpp
+++ b/source/module_io/input_conv.cpp
@@ -328,10 +328,12 @@ void Input_Conv::Convert()
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Hf;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
             {"alpha", "1"},
             {"Rcut_type", "spencer"},
             {"lambda", std::to_string(PARAM.inp.exx_lambda)} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 	// use the error function erf(w|r-r'|), exx just has the short-range part
     else if (dft_functional_lower == "hse"
@@ -339,24 +341,28 @@ void Input_Conv::Convert()
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erfc;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
             {"alpha", "1"},
 			{"omega", std::to_string(PARAM.inp.exx_hse_omega)},
             {"Rcut_type", "limits"} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 	// use the error function erf(w|r-r'|), exx just has the long-range part
     else if ( dft_functional_lower == "wp22" )
     {
         GlobalC::exx_info.info_global.cal_exx = true;
         GlobalC::exx_info.info_global.ccp_type = Conv_Coulomb_Pot_K::Ccp_Type::Erf;
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
+        std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>> coulomb_param;
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Fock] = {{
             {"alpha", "1"},
             {"Rcut_type", "spencer"},
             {"lambda", std::to_string(PARAM.inp.exx_lambda)} }};
-        GlobalC::exx_info.info_global.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
+        coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{
             {"alpha", "-1"},
 			{"omega", std::to_string(PARAM.inp.exx_hse_omega)},
             {"Rcut_type", "limits"} }};
+        GlobalC::exx_info.info_global.coulomb_settings[Conv_Coulomb_Pot_K::Coulomb_Method::Center2] = std::make_pair(true, coulomb_param);
     }
 #ifdef __EXX
     else if (dft_functional_lower == "opt_orb")
diff --git a/source/module_ri/Exx_LRI.h b/source/module_ri/Exx_LRI.h
index ae53393ee6..079a5b0ffe 100644
--- a/source/module_ri/Exx_LRI.h
+++ b/source/module_ri/Exx_LRI.h
@@ -15,6 +15,7 @@
 #include <vector>
 #include <array>
 #include <map>
+#include <unordered_map>
 #include <deque>
 #include <mpi.h>
 
@@ -37,6 +38,15 @@
 		class OperatorLREXX;
 	}
 
+template<typename Tdata>
+class Exx_Obj
+{
+	// match with Conv_Coulomb_Pot_K::Coulomb_Method
+	public:
+		LRI_CV<Tdata> cv;
+		std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
+};
+
 template<typename Tdata>
 class Exx_LRI
 {
@@ -85,9 +95,9 @@ class Exx_LRI
 
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
-	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
-
-	LRI_CV<Tdata> cv;
+	//std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
+	std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Method, Exx_Obj<Tdata>> exx_objs;
+	//LRI_CV<Tdata> cv;
 	RI::Exx<TA,Tcell,Ndim,Tdata> exx_lri;
 
 	void post_process_Hexx( std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> &Hexxs_io ) const;
diff --git a/source/module_ri/Exx_LRI.hpp b/source/module_ri/Exx_LRI.hpp
index 46663f0d75..8c370e71cb 100644
--- a/source/module_ri/Exx_LRI.hpp
+++ b/source/module_ri/Exx_LRI.hpp
@@ -48,19 +48,18 @@ void Exx_LRI<Tdata>::init(const MPI_Comm &mpi_comm_in,
 		{ this->abfs = Exx_Abfs::IO::construct_abfs( abfs_same_atom, orb, this->info.files_abfs, this->info.kmesh_times ); 	}
 	Exx_Abfs::Construct_Orbs::print_orbs_size(ucell, this->abfs, GlobalV::ofs_running);
 
-	const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>>
-		coulomb_param_updated = RI_Util::update_coulomb_param(this->info.coulomb_param, ucell.omega, this->p_kv->get_nkstot_full());
-	this->abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, coulomb_param_updated, this->info.ccp_rmesh_times);
-
 	for( size_t T=0; T!=this->abfs.size(); ++T )
 		{ GlobalC::exx_info.info_ri.abfs_Lmax = std::max( GlobalC::exx_info.info_ri.abfs_Lmax, static_cast<int>(this->abfs[T].size())-1 ); }
 
-	this->cv.set_orbitals(
-		ucell,
-		orb,
-		this->lcaos, this->abfs, this->abfs_ccp,
-		this->info.kmesh_times, this->info.ccp_rmesh_times );
-
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		const std::unordered_map<Conv_Coulomb_Pot_K::Coulomb_Type, std::vector<std::map<std::string,std::string>>>
+			coulomb_param_updated = RI_Util::update_coulomb_param(settings_list.second.second, ucell.omega, this->p_kv->get_nkstot_full());
+		this->exx_objs[settings_list.first].abfs_ccp = Conv_Coulomb_Pot_K::cal_orbs_ccp(this->abfs, coulomb_param_updated, this->info.ccp_rmesh_times);
+		this->exx_objs[settings_list.first].cv.set_orbitals(ucell, orb,
+															this->lcaos, this->abfs, this->exx_objs[settings_list.first].abfs_ccp,
+															this->info.kmesh_times, this->info.ccp_rmesh_times );
+	}
 	ModuleBase::timer::tick("Exx_LRI", "init");
 }
 
@@ -95,26 +94,40 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
 		list_As_Vs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Vs, 2, false);
 
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
-		Vs = this->cv.cal_Vs(ucell,
-			list_As_Vs.first, list_As_Vs.second[0],
-			{{"writable_Vws",true}});
-	this->cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs);
+	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Vs;
+	std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>> dVs;
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>
+			Vs_temp = this->exx_objs[settings_list.first].cv.cal_Vs(ucell,
+				list_As_Vs.first, list_As_Vs.second[0],
+				{{"writable_Vws",true}});
+		this->exx_objs[settings_list.first].cv.Vws = LRI_CV_Tools::get_CVws(ucell,Vs_temp);
+		Vs = Vs.empty() ? Vs_temp : LRI_CV_Tools::add(Vs, Vs_temp);
+
+		if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
+		{
+			std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, Ndim>>>
+				dVs_temp = this->exx_objs[settings_list.first].cv.cal_dVs(ucell,
+					list_As_Vs.first, list_As_Vs.second[0],
+					{{"writable_dVws",true}});
+			this->exx_objs[settings_list.first].cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs_temp);
+			dVs = dVs.empty() ? dVs_temp : LRI_CV_Tools::add(dVs, dVs_temp);
+		}
+		
+	}
 	if (write_cv && GlobalV::MY_RANK == 0)
 		{ LRI_CV_Tools::write_Vs_abf(Vs, PARAM.globalv.global_out_dir + "Vs"); }
 	this->exx_lri.set_Vs(std::move(Vs), this->info.V_threshold);
 
 	if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
 	{
-		std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
-			dVs = this->cv.cal_dVs(ucell,
-				list_As_Vs.first, list_As_Vs.second[0],
-				{{"writable_dVws",true}});
-		this->cv.dVws = LRI_CV_Tools::get_dCVws(ucell,dVs);
-		this->exx_lri.set_dVs(std::move(dVs), this->info.V_grad_threshold);
+		std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
+			dVs_order = LRI_CV_Tools::change_order(std::move(dVs));
+		this->exx_lri.set_dVs(std::move(dVs_order), this->info.V_grad_threshold);
 		if(PARAM.inp.cal_stress)
 		{
-			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs);
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dVRs = LRI_CV_Tools::cal_dMRs(ucell,dVs_order);
 			this->exx_lri.set_dVRs(std::move(dVRs), this->info.V_grad_R_threshold);
 		}
 	}
@@ -123,29 +136,47 @@ void Exx_LRI<Tdata>::cal_exx_ions(const UnitCell& ucell,
 	const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA,std::array<Tcell,Ndim>>>>>
 		list_As_Cs = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
 
-	std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>, std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
-		Cs_dCs = this->cv.cal_Cs_dCs(
-			ucell,
-			list_As_Cs.first, list_As_Cs.second[0],
-			{{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
-			 {"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}});
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs = std::get<0>(Cs_dCs);
-	this->cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs);
+	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Cs;
+	std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;
+	for(const auto &settings_list : this->info.coulomb_settings)
+	{
+		if(settings_list.second.first)
+		{
+			std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+				std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
+					Cs_dCs = this->exx_objs[settings_list.first].cv.cal_Cs_dCs(
+						ucell,
+						list_As_Cs.first, list_As_Cs.second[0],
+						{{"cal_dC",PARAM.inp.cal_force||PARAM.inp.cal_stress},
+						{"writable_Cws",true}, {"writable_dCws",true}, {"writable_Vws",false}, {"writable_dVws",false}});
+			std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> &Cs_temp = std::get<0>(Cs_dCs);
+			this->exx_objs[settings_list.first].cv.Cws = LRI_CV_Tools::get_CVws(ucell,Cs_temp);
+			Cs = Cs.empty() ? Cs_temp : LRI_CV_Tools::add(Cs, Cs_temp);
+		
+			if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
+			{
+				std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> &dCs_temp = std::get<1>(Cs_dCs);
+				this->exx_objs[settings_list.first].cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs_temp);
+				dCs = dCs.empty() ? dCs_temp : LRI_CV_Tools::add(dCs, dCs_temp);
+			}
+		}
+	}
 	if (write_cv && GlobalV::MY_RANK == 0)
 		{ LRI_CV_Tools::write_Cs_ao(Cs, PARAM.globalv.global_out_dir + "Cs"); }
 	this->exx_lri.set_Cs(std::move(Cs), this->info.C_threshold);
 
 	if(PARAM.inp.cal_force || PARAM.inp.cal_stress)
 	{
-		std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3> &dCs = std::get<1>(Cs_dCs);
-		this->cv.dCws = LRI_CV_Tools::get_dCVws(ucell,dCs);
-		this->exx_lri.set_dCs(std::move(dCs), this->info.C_grad_threshold);
+		std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, Ndim>
+			dCs_order = LRI_CV_Tools::change_order(std::move(dCs));
+		this->exx_lri.set_dCs(std::move(dCs_order), this->info.C_grad_threshold);
 		if(PARAM.inp.cal_stress)
 		{
-			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs);
+			std::array<std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>,3> dCRs = LRI_CV_Tools::cal_dMRs(ucell,dCs_order);
 			this->exx_lri.set_dCRs(std::move(dCRs), this->info.C_grad_R_threshold);
 		}
 	}
+	
 	ModuleBase::timer::tick("Exx_LRI", "cal_exx_ions");
 }
 
diff --git a/source/module_ri/LRI_CV.h b/source/module_ri/LRI_CV.h
index 99b52476ea..e05f27ed5f 100644
--- a/source/module_ri/LRI_CV.h
+++ b/source/module_ri/LRI_CV.h
@@ -47,14 +47,14 @@ class LRI_CV
 		const std::vector<TA> &list_A0,
 		const std::vector<TAC> &list_A1,
 		const std::map<std::string,bool> &flags);						// "writable_Vws"
-	inline std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
+	inline std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
 	cal_dVs(
 		const UnitCell &ucell,
 		const std::vector<TA> &list_A0,
 		const std::vector<TAC> &list_A1,
 		const std::map<std::string,bool> &flags);						// "writable_dVws"
-	std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,
-	          std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
+    std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+              std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
 	cal_Cs_dCs(
 		const UnitCell &ucell,
 		const std::vector<TA> &list_A0,
diff --git a/source/module_ri/LRI_CV.hpp b/source/module_ri/LRI_CV.hpp
index 2e4dc71e64..2147b75c35 100644
--- a/source/module_ri/LRI_CV.hpp
+++ b/source/module_ri/LRI_CV.hpp
@@ -146,15 +146,14 @@ auto LRI_CV<Tdata>::cal_dVs(
 	const std::vector<TA> &list_A0,
 	const std::vector<TAC> &list_A1,
 	const std::map<std::string,bool> &flags)					// + "writable_dVws"
--> std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>
+-> std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>
 {
 	ModuleBase::TITLE("LRI_CV","cal_dVs");
 	const T_func_DPcal_data<std::array<RI::Tensor<Tdata>,3>>
 		func_DPcal_dV = std::bind(
 			&LRI_CV<Tdata>::DPcal_dV, this,
 			std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
-	return LRI_CV_Tools::change_order(
-		this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV));
+	return this->cal_datas(ucell,list_A0, list_A1, flags, this->ccp_rmesh_times, func_DPcal_dV);
 }
 
 template<typename Tdata>
@@ -163,7 +162,9 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
 	const std::vector<TA> &list_A0,
 	const std::vector<TAC> &list_A1,
 	const std::map<std::string,bool> &flags)					// "cal_dC" + "writable_Cws", "writable_dCws", "writable_Vws", "writable_dVws"
--> std::pair<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>, std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3>>
+-> std::pair<
+        std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
+        std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
 {
 	ModuleBase::TITLE("LRI_CV","cal_Cs_dCs");
 	const T_func_DPcal_data<std::pair<RI::Tensor<Tdata>, std::array<RI::Tensor<Tdata>,3>>>
@@ -173,16 +174,16 @@ auto LRI_CV<Tdata>::cal_Cs_dCs(
 	std::map<TA,std::map<TAC, std::pair<RI::Tensor<Tdata>, std::array<RI::Tensor<Tdata>,3>>>>
 		Cs_dCs_tmp = this->cal_datas(ucell,list_A0, list_A1, flags, std::min(1.0,this->ccp_rmesh_times), func_DPcal_C_dC);
 
-	std::map<TA,std::map<TAC,RI::Tensor<Tdata>>> Cs;
-	std::array<std::map<TA,std::map<TAC,RI::Tensor<Tdata>>>,3> dCs;
-	for(auto &Cs_dCs_A : Cs_dCs_tmp)
-		for(auto &Cs_dCs_B : Cs_dCs_A.second)
-		{
-			Cs[Cs_dCs_A.first][Cs_dCs_B.first] = std::move(std::get<0>(Cs_dCs_B.second));
-			if(flags.at("cal_dC"))
-				for(int ix=0; ix<3; ++ix)
-					dCs[ix][Cs_dCs_A.first][Cs_dCs_B.first] = std::move(std::get<1>(Cs_dCs_B.second)[ix]);
-		}
+	std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Cs;
+    std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>> dCs;
+    for (auto& Cs_dCs_A: Cs_dCs_tmp)
+        for (auto& Cs_dCs_B: Cs_dCs_A.second) {
+            Cs[Cs_dCs_A.first][Cs_dCs_B.first]
+                = std::move(std::get<0>(Cs_dCs_B.second));
+            if (flags.at("cal_dC"))
+                dCs[Cs_dCs_A.first][Cs_dCs_B.first]
+                    = std::move(std::get<1>(Cs_dCs_B.second));
+        }
 	return std::make_pair(Cs, dCs);
 }
 
diff --git a/source/module_ri/LRI_CV_Tools.h b/source/module_ri/LRI_CV_Tools.h
index c97f502516..060552fcf4 100644
--- a/source/module_ri/LRI_CV_Tools.h
+++ b/source/module_ri/LRI_CV_Tools.h
@@ -6,114 +6,264 @@
 #ifndef LRI_CV_TOOLS_H
 #define LRI_CV_TOOLS_H
 
+#include "abfs.h"
 #include "module_base/abfs-vector3_order.h"
+#include "module_ri/abfs.h"
 
+#include <RI/global/Array_Operator.h>
 #include <RI/global/Tensor.h>
-
-#include <cstddef>
 #include <array>
-#include <vector>
+#include <cstddef>
 #include <map>
-
-#include "module_ri/abfs.h"
+#include <vector>
 
 namespace LRI_CV_Tools
 {
-	template<typename Tdata> extern RI::Tensor<Tdata>                           cal_I( const RI::Tensor<Tdata>                           &m  );
-	template<typename Tdata> extern std::vector<std::vector<RI::Tensor<Tdata>>> cal_I( const std::vector<std::vector<RI::Tensor<Tdata>>> &ms );
-
-	template<typename Tdata> inline RI::Tensor<Tdata>               transform_Rm(const RI::Tensor<Tdata>               &V );
-	template<typename Tdata> inline std::array<RI::Tensor<Tdata>,3> transform_Rm(const std::array<RI::Tensor<Tdata>,3> &dV);
-
-	//template<typename T> inline bool exist(const T &V);
-
-	//template<typename T1, typename T2, typename Treturn>
-	//extern Treturn mul1(const T1 &t1, const T2 &t2);
-	//template<typename T1, typename T2, typename Treturn>
-	//extern Treturn mul2(const T1 &mat, const T2 &vec);
-
-    template<typename Tdata> inline bool exist(const RI::Tensor<Tdata> &V);
-    template<typename T, std::size_t N> inline bool exist(const std::array<T,N> &dV);
-
-    template<typename Tdata>
-    extern RI::Tensor<Tdata> mul1(const RI::Tensor<Tdata> &t1, const RI::Tensor<Tdata> &t2);
-    template<typename T>
-    extern std::array<T,3> mul1(const std::array<T,3> &t1, const T &t2);
-
-    template<typename Tdata>
-    extern std::vector<RI::Tensor<Tdata>> mul2(const std::vector<std::vector<RI::Tensor<Tdata>>> &mat,
-                                                             const std::vector<RI::Tensor<Tdata>> &vec);
-    template<typename T1, typename T2>
-    extern std::array<T2,3> mul2(const T1 &t1, const std::array<T2,3> &t2);
-
-	//template<typename T, std::size_t N>
-	//std::array<T,N> operator-(const std::array<T,N> &v1, const std::array<T,N> &v2);
-	//template<typename T>
-	//std::vector<T> operator-(const std::vector<T> &v1, const std::vector<T> &v2);
-	template<typename T, std::size_t N>
-	extern std::vector<std::array<T,N>> minus(
-		const std::vector<std::array<T,N>> &v1,
-		const std::vector<std::array<T,N>> &v2);
-
-	template<typename T, std::size_t N>
-	extern std::array<T,N> negative(const std::array<T,N> &v_in);
-
-	//template<typename T> T transpose12(const T &c_in);
-    template<typename Tdata> RI::Tensor<Tdata> transpose12(const RI::Tensor<Tdata> &c_in);
-    template<typename T, std::size_t N> std::array<T,N> transpose12(const std::array<T,N> &c_in);
-
-	template<typename T, std::size_t N>
-	extern std::array<std::vector<T>,N>
-	change_order(
-		std::vector<std::array<T,N>> &&ds_in);
-	template<typename T, std::size_t N>
-	std::vector<std::array<T,N>>
-	change_order(
-		std::array<std::vector<T>,N> &&ds_in);
-	template<typename T, std::size_t N>
-	extern std::array<std::vector<std::vector<T>>,N>
-	change_order(
-		std::vector<std::vector<std::array<T,N>>> &&ds_in);
-	template<typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
-	extern std::array<std::map<TkeyA,std::map<TkeyB,Tvalue>>,N>
-	change_order(
-		std::map<TkeyA,std::map<TkeyB,std::array<Tvalue,N>>> && ds_in);
-
-	template<typename Tcell>
-	extern std::array<Tcell,3> cal_latvec_range(const double &rcut_times, 
-												const UnitCell &ucell,
-												const std::vector<double>& orb_cutoff);
-
-	template<typename TA, typename Tcell, typename Tdata>
-	extern std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,RI::Tensor<Tdata>>>>
-	get_CVws(
-		const UnitCell &ucell,
-		const std::map<TA,std::map<std::pair<TA,std::array<Tcell,3>>,RI::Tensor<Tdata>>> &CVs);
-	template<typename TA, typename Tcell, typename Tdata>
-	extern std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,std::array<RI::Tensor<Tdata>,3>>>>
-	get_dCVws(
-		const UnitCell &ucell,
-		const std::array<std::map<TA,std::map<std::pair<TA,std::array<Tcell,3>>,RI::Tensor<Tdata>>>,3> &dCVs);
-
-	template<typename TA, typename TC, typename Tdata>
-	extern std::array<std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3>,3>
-	cal_dMRs(
-		const UnitCell &ucell,
-		const std::array<std::map<TA,std::map<std::pair<TA,TC>,RI::Tensor<Tdata>>>,3> &dMs);
-
-    using TC = std::array<int, 3>;
-    using TAC = std::pair<int, TC>;
-    template <typename T>
-    using TLRI = std::map<int, std::map<TAC, RI::Tensor<T>>>;
-    template <typename T>
-    TLRI<T> read_Cs_ao(const std::string& file_path, const double& threshold = 1e-10);
-    template <typename T>
-    void write_Cs_ao(const TLRI<T>& Vs, const std::string& file_path);
-    template <typename T>
-    TLRI<T> read_Vs_abf(const std::string& file_path, const double& threshold = 1e-10);
-    template <typename T>
-    void write_Vs_abf(const TLRI<T>& Vs, const std::string& file_path);
-}
+template <typename Tdata>
+extern RI::Tensor<Tdata> cal_I(const RI::Tensor<Tdata>& m);
+template <typename Tdata>
+extern std::vector<std::vector<RI::Tensor<Tdata>>> cal_I(const std::vector<std::vector<RI::Tensor<Tdata>>>& ms);
+
+template <typename Tdata>
+inline RI::Tensor<Tdata> transform_Rm(const RI::Tensor<Tdata>& V);
+template <typename Tdata>
+inline std::array<RI::Tensor<Tdata>, 3> transform_Rm(const std::array<RI::Tensor<Tdata>, 3>& dV);
+
+// template<typename T> inline bool exist(const T &V);
+
+// template<typename T1, typename T2, typename Treturn>
+// extern Treturn mul1(const T1 &t1, const T2 &t2);
+// template<typename T1, typename T2, typename Treturn>
+// extern Treturn mul2(const T1 &mat, const T2 &vec);
+
+template <typename Tdata>
+inline bool exist(const RI::Tensor<Tdata>& V);
+template <typename T, std::size_t N>
+inline bool exist(const std::array<T, N>& dV);
+
+template <typename Tdata>
+extern RI::Tensor<Tdata> mul1(const RI::Tensor<Tdata>& t1, const RI::Tensor<Tdata>& t2);
+template <typename T>
+extern std::array<T, 3> mul1(const std::array<T, 3>& t1, const T& t2);
+
+template <typename Tdata>
+extern std::vector<RI::Tensor<Tdata>> mul2(const std::vector<std::vector<RI::Tensor<Tdata>>>& mat,
+                                           const std::vector<RI::Tensor<Tdata>>& vec);
+template <typename T1, typename T2>
+extern std::array<T2, 3> mul2(const T1& t1, const std::array<T2, 3>& t2);
+template <typename T>
+extern RI::Tensor<T> mul2(const T& t1, const RI::Tensor<T>& t2);
+template <typename T, typename TkeyA, typename TkeyB, typename Tvalue>
+extern std::map<TkeyA, std::map<TkeyB, Tvalue>> mul2(const T& t1, const std::map<TkeyA, std::map<TkeyB, Tvalue>>& t2);
+
+// template<typename T, std::size_t N>
+// std::array<T,N> operator-(const std::array<T,N> &v1, const std::array<T,N>
+// &v2); template<typename T> std::vector<T> operator-(const std::vector<T> &v1,
+// const std::vector<T> &v2);
+template <typename T, std::size_t N>
+extern std::vector<std::array<T, N>> minus(const std::vector<std::array<T, N>>& v1,
+                                           const std::vector<std::array<T, N>>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+extern std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> minus(
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v1,
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+inline std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> minus(
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v1,
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue>
+extern std::map<TkeyA, std::map<TkeyB, Tvalue>> minus(std::map<TkeyA, std::map<TkeyB, Tvalue>>& v1,
+                                                      std::map<TkeyA, std::map<TkeyB, Tvalue>>& v2);
+
+template <typename T, std::size_t N>
+extern std::vector<std::array<T, N>> add(const std::vector<std::array<T, N>>& v1,
+                                         const std::vector<std::array<T, N>>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+extern std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> add(
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v1,
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+inline std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> add(
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v1,
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v2);
+template <typename TkeyA, typename TkeyB, typename Tvalue>
+extern std::map<TkeyA, std::map<TkeyB, Tvalue>> add(std::map<TkeyA, std::map<TkeyB, Tvalue>>& v1,
+                                                    std::map<TkeyA, std::map<TkeyB, Tvalue>>& v2);
+
+template <typename T, std::size_t N>
+extern std::array<T, N> negative(const std::array<T, N>& v_in);
+
+// template<typename T> T transpose12(const T &c_in);
+template <typename Tdata>
+RI::Tensor<Tdata> transpose12(const RI::Tensor<Tdata>& c_in);
+template <typename T, std::size_t N>
+std::array<T, N> transpose12(const std::array<T, N>& c_in);
+
+template <typename T, std::size_t N>
+extern std::array<std::vector<T>, N> change_order(std::vector<std::array<T, N>>&& ds_in);
+template <typename T, std::size_t N>
+std::vector<std::array<T, N>> change_order(std::array<std::vector<T>, N>&& ds_in);
+template <typename T, std::size_t N>
+extern std::array<std::vector<std::vector<T>>, N> change_order(std::vector<std::vector<std::array<T, N>>>&& ds_in);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+extern std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> change_order(
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>&& ds_in);
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+extern std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> change_order(
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>&& ds_in);
+
+template <typename Tcell>
+extern std::array<Tcell, 3> cal_latvec_range(const double& rcut_times,
+                                             const UnitCell& ucell,
+                                             const std::vector<double>& orb_cutoff);
+
+template <typename TA, typename Tcell, typename Tdata>
+extern std::map<int, std::map<int, std::map<Abfs::Vector3_Order<double>, RI::Tensor<Tdata>>>> get_CVws(
+    const UnitCell& ucell,
+    const std::map<TA, std::map<std::pair<TA, std::array<Tcell, 3>>, RI::Tensor<Tdata>>>& CVs);
+template <typename TA, typename Tcell, typename Tdata>
+extern std::map<int, std::map<int, std::map<Abfs::Vector3_Order<double>, std::array<RI::Tensor<Tdata>, 3>>>> get_dCVws(
+    const UnitCell& ucell,
+    const std::map<TA, std::map<std::pair<TA, std::array<Tcell, 3>>, std::array<RI::Tensor<Tdata>, 3>>>& dCVs);
+template <typename TA, typename TC, typename Tdata>
+extern std::array<std::array<std::map<TA, std::map<std::pair<TA, TC>, RI::Tensor<Tdata>>>, 3>, 3> cal_dMRs(
+    const UnitCell& ucell,
+    const std::array<std::map<TA, std::map<std::pair<TA, TC>, RI::Tensor<Tdata>>>, 3>& dMs);
+
+using TC = std::array<int, 3>;
+using TAC = std::pair<int, TC>;
+template <typename T>
+using TLRI = std::map<int, std::map<TAC, RI::Tensor<T>>>;
+template <typename T>
+TLRI<T> read_Cs_ao(const std::string& file_path, const double& threshold = 1e-10);
+template <typename T>
+void write_Cs_ao(const TLRI<T>& Vs, const std::string& file_path);
+template <typename T>
+TLRI<T> read_Vs_abf(const std::string& file_path, const double& threshold = 1e-10);
+template <typename T>
+void write_Vs_abf(const TLRI<T>& Vs, const std::string& file_path);
+
+template <typename T>
+struct is_std_array : std::false_type
+{
+};
+template <typename T, std::size_t N>
+struct is_std_array<std::array<T, N>> : std::true_type
+{
+};
+template <typename T>
+struct is_tensor : std::false_type
+{
+};
+template <typename T>
+struct is_tensor<RI::Tensor<T>> : std::true_type
+{
+};
+
+template <typename Tout>
+struct TinType;
+
+template <typename T>
+struct TinType<RI::Tensor<T>>
+{
+    using type = T;
+};
+
+template <typename T, std::size_t N>
+struct TinType<std::array<RI::Tensor<T>, N>>
+{
+    using type = T;
+};
+
+template <typename Tdata, typename = std::enable_if_t<!is_std_array<Tdata>::value>>
+inline void init_elem(Tdata& data, const size_t ndim0, const size_t ndim1)
+{
+    data = Tdata({ndim0, ndim1});
+};
+template <typename T, std::size_t N>
+extern void init_elem(std::array<RI::Tensor<T>, N>& data, const size_t ndim0, const size_t ndim1);
+
+template <typename Tdata, typename = std::enable_if_t<!is_std_array<Tdata>::value && !is_tensor<Tdata>::value>>
+inline void add_elem(Tdata& data, const Tdata& val, const Tdata& frac)
+{
+    data += frac * val;
+};
+template <typename T, std::size_t N>
+extern void add_elem(std::array<T, N>& data, const T& val, const T& frac);
+template <typename Tdata, typename = std::enable_if_t<is_tensor<Tdata>::value>>
+inline void add_elem(const Tdata& data,
+                     const int lmp,
+                     const int lmq,
+                     const typename TinType<Tdata>::type& val,
+                     const typename TinType<Tdata>::type& frac)
+{
+    data(lmp, lmq) += frac * val;
+};
+template <typename T, std::size_t N>
+extern void add_elem(std::array<RI::Tensor<T>, N>& data,
+                     const int lmp,
+                     const int lmq,
+                     const std::array<T, N>& val,
+                     const T& frac);
+template <typename Tdata, typename = std::enable_if_t<is_tensor<Tdata>::value>>
+inline void add_elem(Tdata& data,
+                     const int lmp0,
+                     const int lmq0,
+                     const Tdata& val,
+                     const int lmp1,
+                     const int lmq1,
+                     const typename TinType<Tdata>::type& frac)
+{
+    data(lmp0, lmq0) += frac * val(lmp1, lmq1);
+};
+template <typename T, std::size_t N>
+extern void add_elem(std::array<RI::Tensor<T>, N>& data,
+                     const int lmp0,
+                     const int lmq0,
+                     const std::array<RI::Tensor<T>, N>& val,
+                     const int lmp1,
+                     const int lmq1,
+                     const T& frac);
+
+template <typename Tout, typename Tin>
+inline RI::Tensor<Tout> convert(RI::Tensor<Tin>&& data);
+template <typename Tout, typename Tin, std::size_t N>
+extern std::array<RI::Tensor<Tout>, N> convert(std::array<RI::Tensor<Tin>, N>&& data);
+
+// template <typename T>
+// typename std::enable_if<!RI::Global_Func::is_complex<T>::value, T>::type
+// inline check_zero(T value) {
+//     return (std::abs(value) < 1e-8) ? static_cast<T>(0) : value;
+// }
+
+// template <typename T>
+// typename std::enable_if<RI::Global_Func::is_complex<T>::value, T>::type
+// inline check_zero(const T& value) {
+//     using RealType = typename T::value_type;
+//     RealType real_part = std::real(value);
+//     RealType imag_part = std::imag(value);
+
+//     real_part = (std::abs(real_part) < 1e-8) ? 0 : real_part;
+//     imag_part = (std::abs(imag_part) < 1e-8) ? 0 : imag_part;
+
+//     return std::complex<RealType>(real_part, imag_part);
+// }
+
+// template <typename T>
+// extern RI::Tensor<T> check_zero(RI::Tensor<T>&& data);
+// template <typename T, std::size_t N>
+// extern std::array<RI::Tensor<T>, N> check_zero(std::array<RI::Tensor<T>, N>&& data);
+
+template <typename T>
+struct plus
+{
+    T operator()(const T& lhs, const T& rhs) const
+    {
+        using namespace RI::Array_Operator;
+        return lhs + rhs;
+    }
+};
+} // namespace LRI_CV_Tools
 
 #include "LRI_CV_Tools.hpp"
 #include "write_ri_cv.hpp"
diff --git a/source/module_ri/LRI_CV_Tools.hpp b/source/module_ri/LRI_CV_Tools.hpp
index 39087e6641..24d3493e04 100644
--- a/source/module_ri/LRI_CV_Tools.hpp
+++ b/source/module_ri/LRI_CV_Tools.hpp
@@ -6,251 +6,366 @@
 #ifndef LRI_CV_TOOLS_HPP
 #define LRI_CV_TOOLS_HPP
 
-#include "LRI_CV_Tools.h"
-#include "Inverse_Matrix.h"
 #include "../module_base/mathzone.h"
+#include "Inverse_Matrix.h"
+#include "LRI_CV_Tools.h"
+#include "RI_Util.h"
+
+#include <RI/global/Global_Func-1.h>
+#include <RI/global/Map_Operator.h>
 #include "../module_hamilt_pw/hamilt_pwdft/global.h"
 
-template<typename Tdata>
-RI::Tensor<Tdata>
-LRI_CV_Tools::cal_I( const RI::Tensor<Tdata> &m )
-{
-	Inverse_Matrix<Tdata> I;
-	I.input(m);
-	I.cal_inverse( Inverse_Matrix<Tdata>::Method::potrf );
-	return I.output();
+template <typename Tdata>
+RI::Tensor<Tdata> LRI_CV_Tools::cal_I(const RI::Tensor<Tdata>& m) {
+    Inverse_Matrix<Tdata> I;
+    I.input(m);
+    I.cal_inverse(Inverse_Matrix<Tdata>::Method::potrf);
+    return I.output();
 }
 
-template<typename Tdata>
+template <typename Tdata>
 std::vector<std::vector<RI::Tensor<Tdata>>>
-LRI_CV_Tools::cal_I( const std::vector<std::vector<RI::Tensor<Tdata>>> &ms )
-{
-	Inverse_Matrix<Tdata> I;
-	I.input(ms);
-	I.cal_inverse( Inverse_Matrix<Tdata>::Method::potrf );
-	return I.output({ms[0][0].shape[0], ms[1][0].shape[0]}, {ms[0][0].shape[1], ms[0][1].shape[1]});
+    LRI_CV_Tools::cal_I(const std::vector<std::vector<RI::Tensor<Tdata>>>& ms) {
+    Inverse_Matrix<Tdata> I;
+    I.input(ms);
+    I.cal_inverse(Inverse_Matrix<Tdata>::Method::potrf);
+    return I.output({ms[0][0].shape[0], ms[1][0].shape[0]},
+                    {ms[0][0].shape[1], ms[0][1].shape[1]});
 }
 
-
-
-template<typename Tdata>
-RI::Tensor<Tdata> LRI_CV_Tools::transform_Rm(const RI::Tensor<Tdata> &V)
-{
-	return V.transpose();
+template <typename Tdata>
+RI::Tensor<Tdata> LRI_CV_Tools::transform_Rm(const RI::Tensor<Tdata>& V) {
+    return V.transpose();
 }
 
-template<typename Tdata>
-std::array<RI::Tensor<Tdata>,3> LRI_CV_Tools::transform_Rm(const std::array<RI::Tensor<Tdata>,3> &dV)
-{
-	return std::array<RI::Tensor<Tdata>,3>{-dV[0].transpose(), -dV[1].transpose(), -dV[2].transpose()};
+template <typename Tdata>
+std::array<RI::Tensor<Tdata>, 3>
+    LRI_CV_Tools::transform_Rm(const std::array<RI::Tensor<Tdata>, 3>& dV) {
+    return std::array<RI::Tensor<Tdata>, 3>{-dV[0].transpose(),
+                                            -dV[1].transpose(),
+                                            -dV[2].transpose()};
 }
 
-template<typename Tdata>
-bool LRI_CV_Tools::exist(const RI::Tensor<Tdata> &V)
-{
-	return !V.empty();
+template <typename Tdata>
+bool LRI_CV_Tools::exist(const RI::Tensor<Tdata>& V) {
+    return !V.empty();
 }
 
-template<typename T, std::size_t N>
-bool LRI_CV_Tools::exist(const std::array<T,N> &dV)
-{
-	for(size_t i=0; i<3; ++i)
-	{
-		if(!dV[i].empty())
-			return true;
-	}
-	return false;
+template <typename T, std::size_t N>
+bool LRI_CV_Tools::exist(const std::array<T, N>& dV) {
+    for (size_t i = 0; i < 3; ++i)
+        if (!dV[i].empty())
+            return true;
+    return false;
 }
 
-
-template<typename Tdata>
-RI::Tensor<Tdata> LRI_CV_Tools::mul1(
-	const RI::Tensor<Tdata> &t1,
-	const RI::Tensor<Tdata> &t2)
-{
-	const size_t sa0=t1.shape[0], sa1=t2.shape[0], sl0=t2.shape[1], sl1=t2.shape[2];
-	return (t1 * t2.reshape({sa1,sl0*sl1})).reshape({sa0,sl0,sl1});
+template <typename Tdata>
+RI::Tensor<Tdata> LRI_CV_Tools::mul1(const RI::Tensor<Tdata>& t1,
+                                     const RI::Tensor<Tdata>& t2) {
+    const size_t sa0 = t1.shape[0], sa1 = t2.shape[0], sl0 = t2.shape[1],
+                 sl1 = t2.shape[2];
+    return (t1 * t2.reshape({sa1, sl0 * sl1})).reshape({sa0, sl0, sl1});
 }
-template<typename T>
-std::array<T,3> LRI_CV_Tools::mul1(
-	const std::array<T,3> &t1,
-	const T &t2)
-{
-	return std::array<T,3>{
-		mul1(t1[0],t2), mul1(t1[1],t2), mul1(t1[2],t2) };
+template <typename T>
+std::array<T, 3> LRI_CV_Tools::mul1(const std::array<T, 3>& t1, const T& t2) {
+    return std::array<T, 3>{mul1(t1[0], t2), mul1(t1[1], t2), mul1(t1[2], t2)};
 }
 /*
 template<typename T>
 std::array<T,3> LRI_CV_Tools::mul1(
-	const T &t1,
-	const std::array<T,3> &t2)
+    const T &t1,
+    const std::array<T,3> &t2)
 {
-	return std::array<T,3>{
-		mul1(t1,t2[0]), mul1(t1,t2[1]), mul1(t1,t2[2]) };
+    return std::array<T,3>{
+        mul1(t1,t2[0]), mul1(t1,t2[1]), mul1(t1,t2[2]) };
 }
 */
 
-template<typename Tdata>
-std::vector<RI::Tensor<Tdata>> LRI_CV_Tools::mul2(
-	const std::vector<std::vector<RI::Tensor<Tdata>>> &mat,
-	const std::vector<RI::Tensor<Tdata>> &vec)
-{
-	const size_t sa0=vec[0].shape[0], sa1=vec[1].shape[0], sl0=vec[0].shape[1], sl1=vec[0].shape[2];
-	const RI::Tensor<Tdata> vec0=vec[0].reshape({sa0,sl0*sl1}), vec1=vec[1].reshape({sa1,sl0*sl1});
-	return std::vector<RI::Tensor<Tdata>>
-		{( mat[0][0]*vec0 + mat[0][1]*vec1 ).reshape({sa0,sl0,sl1}),
-		 ( mat[1][0]*vec0 + mat[1][1]*vec1 ).reshape({sa1,sl0,sl1})};
+template <typename Tdata>
+std::vector<RI::Tensor<Tdata>>
+    LRI_CV_Tools::mul2(const std::vector<std::vector<RI::Tensor<Tdata>>>& mat,
+                       const std::vector<RI::Tensor<Tdata>>& vec) {
+    const size_t sa0 = vec[0].shape[0], sa1 = vec[1].shape[0],
+                 sl0 = vec[0].shape[1], sl1 = vec[0].shape[2];
+    const RI::Tensor<Tdata> vec0 = vec[0].reshape({sa0, sl0 * sl1}),
+                            vec1 = vec[1].reshape({sa1, sl0 * sl1});
+    return std::vector<RI::Tensor<Tdata>>{
+        (mat[0][0] * vec0 + mat[0][1] * vec1).reshape({sa0, sl0, sl1}),
+        (mat[1][0] * vec0 + mat[1][1] * vec1).reshape({sa1, sl0, sl1})};
 }
 /*
 template<typename T1, typename T2>
 std::array<T2,3> LRI_CV_Tools::mul2(
-	const std::array<T1,3> &t1,
-	const T2 &t2)
+    const std::array<T1,3> &t1,
+    const T2 &t2)
 {
-	return std::array<T2,3>{
-		mul2(t1[0],t2), mul2(t1[1],t2), mul2(t1[2],t2) };
+    return std::array<T2,3>{
+        mul2(t1[0],t2), mul2(t1[1],t2), mul2(t1[2],t2) };
 }
 */
-template<typename T1, typename T2>
-std::array<T2,3> LRI_CV_Tools::mul2(
-	const T1 &t1,
-	const std::array<T2,3> &t2)
-{
-	return std::array<T2,3>{
-		mul2(t1,t2[0]), mul2(t1,t2[1]), mul2(t1,t2[2]) };
+template <typename T1, typename T2>
+std::array<T2, 3> LRI_CV_Tools::mul2(const T1& t1,
+                                     const std::array<T2, 3>& t2) {
+    return std::array<T2, 3>{mul2(t1, t2[0]), mul2(t1, t2[1]), mul2(t1, t2[2])};
+}
+
+template <typename T>
+RI::Tensor<T> LRI_CV_Tools::mul2(const T& t1, const RI::Tensor<T>& t2) {
+    return t1 * t2;
+}
+
+template <typename T, typename TkeyA, typename TkeyB, typename Tvalue>
+std::map<TkeyA, std::map<TkeyB, Tvalue>>
+    LRI_CV_Tools::mul2(const T& t1,
+                       const std::map<TkeyA, std::map<TkeyB, Tvalue>>& t2) {
+    std::map<TkeyA, std::map<TkeyB, Tvalue>> res;
+    for (const auto& outerPair: t2) {
+        const TkeyA keyA = outerPair.first;
+        const std::map<TkeyB, Tvalue>& innerMap = outerPair.second;
+        std::map<TkeyB, Tvalue> newInnerMap;
+
+        for (const auto& innerPair: innerMap) {
+            const TkeyB keyB = innerPair.first;
+            const Tvalue value = innerPair.second;
+            newInnerMap[keyB] = mul2(t1, value);
+        }
+
+        res[keyA] = newInnerMap;
+    }
+
+    return res;
 }
 
 /*
 template<typename T, std::size_t N>
-std::array<T,N> LRI_CV_Tools::operator-(const std::array<T,N> &v1, const std::array<T,N> &v2)
+std::array<T,N> LRI_CV_Tools::operator-(const std::array<T,N> &v1, const
+std::array<T,N> &v2)
 {
-	std::array<T,N> v;
-	for(std::size_t i=0; i<N; ++i)
-		v[i] = v1[i] - v2[i];
-	return v;
+    std::array<T,N> v;
+    for(std::size_t i=0; i<N; ++i)
+        v[i] = v1[i] - v2[i];
+    return v;
 }
 template<typename T>
-std::vector<T> LRI_CV_Tools::operator-(const std::vector<T> &v1, const std::vector<T> &v2)
+std::vector<T> LRI_CV_Tools::operator-(const std::vector<T> &v1, const
+std::vector<T> &v2)
 {
-	assert(v1.size()==v2.size());
-	std::vector<T> v(v1.size());
-	for(std::size_t i=0; i<v.size(); ++i)
-		v[i] = v1[i] - v2[i];
-	return v;
+    assert(v1.size()==v2.size());
+    std::vector<T> v(v1.size());
+    for(std::size_t i=0; i<v.size(); ++i)
+        v[i] = v1[i] - v2[i];
+    return v;
 }
 */
-template<typename T, std::size_t N>
-std::vector<std::array<T,N>> LRI_CV_Tools::minus(
-	const std::vector<std::array<T,N>> &v1,
-	const std::vector<std::array<T,N>> &v2)
-{
-	assert(v1.size()==v2.size());
-	std::vector<std::array<T,N>> v(v1.size());
-	for(std::size_t i=0; i<v.size(); ++i)
-		for(std::size_t j=0; j<N; ++j)
-			v[i][j] = v1[i][j] - v2[i][j];
-	return v;
+template <typename T, std::size_t N>
+std::vector<std::array<T, N>>
+    LRI_CV_Tools::minus(const std::vector<std::array<T, N>>& v1,
+                        const std::vector<std::array<T, N>>& v2) {
+    assert(v1.size() == v2.size());
+    std::vector<std::array<T, N>> v(v1.size());
+    for (std::size_t i = 0; i < v.size(); ++i)
+        for (std::size_t j = 0; j < N; ++j)
+            v[i][j] = v1[i][j] - v2[i][j];
+    return v;
 }
 
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> LRI_CV_Tools::minus(
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v1,
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v2) {
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> v1_order
+        = change_order(std::move(v1));
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> v2_order
+        = change_order(std::move(v2));
+    auto dv = minus(v1_order, v2_order);
+    return change_order(std::move(dv));
+}
 
-template<typename T, std::size_t N>
-std::array<T,N> LRI_CV_Tools::negative(const std::array<T,N> &v_in)
-{
-	std::array<T,N> v_out;
-	for(std::size_t i=0; i<N; ++i)
-		v_out[i] = -v_in[i];
-	return v_out;
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> LRI_CV_Tools::minus(
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v1,
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v2) {
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> dv;
+    for (size_t i = 0; i != N; ++i)
+        dv[i] = minus(v1[i], v2[i]);
+    return dv;
 }
 
+template <typename TkeyA, typename TkeyB, typename Tvalue>
+std::map<TkeyA, std::map<TkeyB, Tvalue>>
+    LRI_CV_Tools::minus(std::map<TkeyA, std::map<TkeyB, Tvalue>>& v1,
+                        std::map<TkeyA, std::map<TkeyB, Tvalue>>& v2) {
+    assert(v1.size() == v2.size());
+    using namespace RI::Map_Operator;
+    using namespace RI::Array_Operator;
 
-template<typename Tdata>
-RI::Tensor<Tdata> LRI_CV_Tools::transpose12(const RI::Tensor<Tdata> &c_in)
-{
-	RI::Tensor<Tdata> c_out({c_in.shape[0], c_in.shape[2], c_in.shape[1]});
-	for(size_t i0=0; i0<c_in.shape[0]; ++i0)
-		for(size_t i1=0; i1<c_in.shape[1]; ++i1)
-			for(size_t i2=0; i2<c_in.shape[2]; ++i2)
-				c_out(i0,i2,i1) = c_in(i0,i1,i2);
-	return c_out;
+    std::map<TkeyA, std::map<TkeyB, Tvalue>> dv;
+    auto it1 = v1.begin();
+    auto it2 = v2.begin();
+    while (it1 != v1.end() && it2 != v2.end()) {
+        assert(it1->first == it2->first);
+        const TkeyA& keyA = it1->first;
+        const std::map<TkeyB, Tvalue>& map1 = it1->second;
+        const std::map<TkeyB, Tvalue>& map2 = it2->second;
+        dv[keyA] = map1 - map2;
+        ++it1;
+        ++it2;
+    }
+    return dv;
 }
 
-template<typename T, std::size_t N>
-std::array<T,N> LRI_CV_Tools::transpose12(const std::array<T,N> &c_in)
-{
-	std::array<T,N> c_out;
-	for(size_t i=0; i<N; ++i)
-		c_out[i] = transpose12(c_in[i]);
-	return c_out;
+template <typename T, std::size_t N>
+std::vector<std::array<T, N>>
+    LRI_CV_Tools::add(const std::vector<std::array<T, N>>& v1,
+                        const std::vector<std::array<T, N>>& v2) {
+    assert(v1.size() == v2.size());
+    std::vector<std::array<T, N>> v(v1.size());
+    for (std::size_t i = 0; i < v.size(); ++i)
+        for (std::size_t j = 0; j < N; ++j)
+            v[i][j] = v1[i][j] + v2[i][j];
+    return v;
 }
 
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> LRI_CV_Tools::add(
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v1,
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>& v2) {
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> v1_order
+        = change_order(std::move(v1));
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> v2_order
+        = change_order(std::move(v2));
+    auto dv = add(v1_order, v2_order);
+    return change_order(std::move(dv));
+}
 
-template<typename T, std::size_t N>
-std::array<std::vector<T>,N>
-LRI_CV_Tools::change_order(std::vector<std::array<T,N>> &&ds_in)
-{
-	std::array<std::vector<T>,N> ds;
-	for(int ix=0; ix<N; ++ix)
-	{
-		ds[ix].resize(ds_in.size());
-		for(int iv=0; iv<ds_in.size(); ++iv)
-			ds[ix][iv] = std::move(ds_in[iv][ix]);
-	}
-	return ds;
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> LRI_CV_Tools::add(
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v1,
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>& v2) {
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> dv;
+    for (size_t i = 0; i != N; ++i)
+        dv[i] = add(v1[i], v2[i]);
+    return dv;
 }
 
-template<typename T, std::size_t N>
-std::vector<std::array<T,N>>
-LRI_CV_Tools::change_order(std::array<std::vector<T>,N> &&ds_in)
-{
-	std::vector<std::array<T,N>> ds(ds_in[0].size());
-	for(int ix=0; ix<N; ++ix)
-	{
-		assert(ds.size() == ds_in[ix].size());
-		for(int iv=0; iv<ds.size(); ++iv)
-			ds[iv][ix] = std::move(ds_in[ix][iv]);
-	}
-	return ds;
+template <typename TkeyA, typename TkeyB, typename Tvalue>
+std::map<TkeyA, std::map<TkeyB, Tvalue>>
+    LRI_CV_Tools::add(std::map<TkeyA, std::map<TkeyB, Tvalue>>& v1,
+                        std::map<TkeyA, std::map<TkeyB, Tvalue>>& v2) {
+    assert(v1.size() == v2.size());
+    using namespace RI::Map_Operator;
+    using namespace RI::Array_Operator;
+
+    std::map<TkeyA, std::map<TkeyB, Tvalue>> dv;
+    auto it1 = v1.begin();
+    auto it2 = v2.begin();
+    while (it1 != v1.end() && it2 != v2.end()) {
+        assert(it1->first == it2->first);
+        const TkeyA& keyA = it1->first;
+        const std::map<TkeyB, Tvalue>& map1 = it1->second;
+        const std::map<TkeyB, Tvalue>& map2 = it2->second;
+        dv[keyA] = map1 + map2;
+        ++it1;
+        ++it2;
+    }
+    return dv;
 }
 
-template<typename T, std::size_t N>
-std::array<std::vector<std::vector<T>>,N>
-LRI_CV_Tools::change_order(std::vector<std::vector<std::array<T,N>>> &&ds_in)
-{
-	std::array<std::vector<std::vector<T>>,N> ds;
-	for(int ix=0; ix<N; ++ix)
-	{
-		ds[ix].resize(ds_in.size());
-		for(int i0=0; i0<ds_in.size(); ++i0)
-		{
-			ds[ix][i0].resize(ds_in[i0].size());
-			for(int i1=0; i1<ds_in[i0].size(); ++i1)
-				ds[ix][i0][i1] = std::move(ds_in[i0][i1][ix]);
-		}
-	}
-	return ds;
+template <typename T, std::size_t N>
+std::array<T, N> LRI_CV_Tools::negative(const std::array<T, N>& v_in) {
+    std::array<T, N> v_out;
+    for (std::size_t i = 0; i < N; ++i)
+        v_out[i] = -v_in[i];
+    return v_out;
 }
 
-template<typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
-std::array<std::map<TkeyA,std::map<TkeyB,Tvalue>>,N>
-LRI_CV_Tools::change_order(std::map<TkeyA,std::map<TkeyB,std::array<Tvalue,N>>> && ds_in)
-{
-	std::array<std::map<TkeyA,std::map<TkeyB,Tvalue>>,N> ds;
-	for(auto &ds_A : ds_in)
-		for(auto &ds_B : ds_A.second)
-			for(int ix=0; ix<N; ++ix)
-				ds[ix][ds_A.first][ds_B.first] = std::move(ds_B.second[ix]);
-	return ds;
+template <typename Tdata>
+RI::Tensor<Tdata> LRI_CV_Tools::transpose12(const RI::Tensor<Tdata>& c_in) {
+    RI::Tensor<Tdata> c_out({c_in.shape[0], c_in.shape[2], c_in.shape[1]});
+    for (size_t i0 = 0; i0 < c_in.shape[0]; ++i0)
+        for (size_t i1 = 0; i1 < c_in.shape[1]; ++i1)
+            for (size_t i2 = 0; i2 < c_in.shape[2]; ++i2)
+                c_out(i0, i2, i1) = c_in(i0, i1, i2);
+    return c_out;
+}
+
+template <typename T, std::size_t N>
+std::array<T, N> LRI_CV_Tools::transpose12(const std::array<T, N>& c_in) {
+    std::array<T, N> c_out;
+    for (size_t i = 0; i < N; ++i)
+        c_out[i] = transpose12(c_in[i]);
+    return c_out;
+}
+
+template <typename T, std::size_t N>
+std::array<std::vector<T>, N>
+    LRI_CV_Tools::change_order(std::vector<std::array<T, N>>&& ds_in) {
+    std::array<std::vector<T>, N> ds;
+    for (int ix = 0; ix < N; ++ix) {
+        ds[ix].resize(ds_in.size());
+        for (int iv = 0; iv < ds_in.size(); ++iv)
+            ds[ix][iv] = std::move(ds_in[iv][ix]);
+    }
+    return ds;
 }
 
+template <typename T, std::size_t N>
+std::vector<std::array<T, N>>
+    LRI_CV_Tools::change_order(std::array<std::vector<T>, N>&& ds_in) {
+    std::vector<std::array<T, N>> ds(ds_in[0].size());
+    for (int ix = 0; ix < N; ++ix) {
+        assert(ds.size() == ds_in[ix].size());
+        for (int iv = 0; iv < ds.size(); ++iv)
+            ds[iv][ix] = std::move(ds_in[ix][iv]);
+    }
+    return ds;
+}
+
+template <typename T, std::size_t N>
+std::array<std::vector<std::vector<T>>, N> LRI_CV_Tools::change_order(
+    std::vector<std::vector<std::array<T, N>>>&& ds_in) {
+    std::array<std::vector<std::vector<T>>, N> ds;
+    for (int ix = 0; ix < N; ++ix) {
+        ds[ix].resize(ds_in.size());
+        for (int i0 = 0; i0 < ds_in.size(); ++i0) {
+            ds[ix][i0].resize(ds_in[i0].size());
+            for (int i1 = 0; i1 < ds_in[i0].size(); ++i1)
+                ds[ix][i0][i1] = std::move(ds_in[i0][i1][ix]);
+        }
+    }
+    return ds;
+}
+
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>
+    LRI_CV_Tools::change_order(
+        std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>&& ds_in) {
+    std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N> ds;
+    for (auto& ds_A: ds_in)
+        for (auto& ds_B: ds_A.second)
+            for (int ix = 0; ix < N; ++ix)
+                ds[ix][ds_A.first][ds_B.first] = std::move(ds_B.second[ix]);
+    return ds;
+}
+
+template <typename TkeyA, typename TkeyB, typename Tvalue, std::size_t N>
+std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>>
+    LRI_CV_Tools::change_order(
+        std::array<std::map<TkeyA, std::map<TkeyB, Tvalue>>, N>&& ds_in) {
+    std::map<TkeyA, std::map<TkeyB, std::array<Tvalue, N>>> ds;
+    for (int ix = 0; ix < N; ++ix)
+        for (auto& ds_A: ds_in[ix])
+            for (auto& ds_B: ds_A.second)
+                ds[ds_A.first][ds_B.first][ix] = std::move(ds_B.second);
+    return ds;
+}
 
-template<typename Tcell>
-std::array<Tcell,3>
-LRI_CV_Tools::cal_latvec_range(const double &rcut_times, 
+template <typename Tcell>
+std::array<Tcell, 3> LRI_CV_Tools::cal_latvec_range(const double& rcut_times,
 							   const UnitCell &ucell,
-							   const std::vector<double>& orb_cutoff)
-{
-	double Rcut_max = 0;
-	for(int T=0; T<ucell.ntype; ++T)
+							   const std::vector<double>& orb_cutoff) {
+    double Rcut_max = 0;
+    for(int T=0; T<ucell.ntype; ++T)
 		Rcut_max = std::max(Rcut_max, orb_cutoff[T]);
 	const ModuleBase::Vector3<double> proj = ModuleBase::Mathzone::latvec_projection(
 		std::array<ModuleBase::Vector3<double>,3>{ucell.a1, ucell.a2, ucell.a3});
@@ -289,36 +404,115 @@ LRI_CV_Tools::get_CVws(
 	return CVws;
 }
 
-template<typename TA, typename Tcell, typename Tdata>
-std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,std::array<RI::Tensor<Tdata>,3>>>>
-LRI_CV_Tools::get_dCVws(
-	const UnitCell &ucell,
-	const std::array<std::map<TA,std::map<std::pair<TA,std::array<Tcell,3>>,RI::Tensor<Tdata>>>,3> &dCVs)
+template <typename TA, typename Tcell, typename Tdata>
+std::map<int, std::map<int, std::map<Abfs::Vector3_Order<double>, std::array<RI::Tensor<Tdata>, 3>>>> LRI_CV_Tools::
+    get_dCVws(const UnitCell& ucell,
+              const std::map<TA, std::map<std::pair<TA, std::array<Tcell, 3>>, std::array<RI::Tensor<Tdata>, 3>>>& dCVs)
 {
-	std::map<int,std::map<int,std::map<Abfs::Vector3_Order<double>,std::array<RI::Tensor<Tdata>,3>>>> dCVws;
-	for(int ix=0; ix<3; ++ix)
-	{
-		for(const auto &dCVs_A : dCVs[ix])
-		{
-			const TA iat0 = dCVs_A.first;
-			const int it0 = ucell.iat2it[iat0];
-			const int ia0 = ucell.iat2ia[iat0];
-			const ModuleBase::Vector3<double> tau0 = ucell.atoms[it0].tau[ia0];
-			for(const auto &dCVs_B : dCVs_A.second)
-			{
-				const TA iat1 = dCVs_B.first.first;
-				const int it1 = ucell.iat2it[iat1];
-				const int ia1 = ucell.iat2ia[iat1];
-				const std::array<int,3> &cell1 = dCVs_B.first.second;
-				const ModuleBase::Vector3<double> tau1 = ucell.atoms[it1].tau[ia1];
-				const Abfs::Vector3_Order<double> R_delta = -tau0+tau1+(RI_Util::array3_to_Vector3(cell1)*ucell.latvec);
-				dCVws[it0][it1][R_delta][ix] = dCVs_B.second;
-			}
-		}
-	}
-	return dCVws;
+    std::map<int, std::map<int, std::map<Abfs::Vector3_Order<double>, std::array<RI::Tensor<Tdata>, 3>>>> dCVws;
+    for (const auto& dCVs_A: dCVs)
+    {
+        const TA iat0 = dCVs_A.first;
+        const int it0 = ucell.iat2it[iat0];
+        const int ia0 = ucell.iat2ia[iat0];
+        const ModuleBase::Vector3<double> tau0 = ucell.atoms[it0].tau[ia0];
+        for (const auto& dCVs_B: dCVs_A.second)
+        {
+            const TA iat1 = dCVs_B.first.first;
+            const int it1 = ucell.iat2it[iat1];
+            const int ia1 = ucell.iat2ia[iat1];
+            const std::array<int, 3>& cell1 = dCVs_B.first.second;
+            const ModuleBase::Vector3<double> tau1 = ucell.atoms[it1].tau[ia1];
+            const Abfs::Vector3_Order<double> R_delta
+                = -tau0 + tau1 + (RI_Util::array3_to_Vector3(cell1) * ucell.latvec);
+            dCVws[it0][it1][R_delta] = dCVs_B.second;
+        }
+    }
+    return dCVws;
+}
+
+template <typename T, std::size_t N>
+void LRI_CV_Tools::init_elem(std::array<RI::Tensor<T>, N>& data,
+                             const size_t ndim0,
+                             const size_t ndim1) {
+    for (size_t i = 0; i < N; ++i) {
+        data[i] = RI::Tensor<T>({ndim0, ndim1});
+    }
+}
+
+template <typename T, std::size_t N>
+void LRI_CV_Tools::add_elem(std::array<T, N>& data,
+                            const T& val,
+                            const T& frac) {
+    for (size_t i = 0; i < N; ++i)
+        data[i] += frac * val;
+}
+
+template <typename T, std::size_t N>
+void LRI_CV_Tools::add_elem(std::array<RI::Tensor<T>, N>& data,
+                            const int lmp,
+                            const int lmq,
+                            const std::array<T, N>& val,
+                            const T& frac) {
+    for (size_t i = 0; i < N; ++i) {
+        data[i](lmp, lmq) += frac * val[i];
+    }
+}
+
+template <typename T, std::size_t N>
+void LRI_CV_Tools::add_elem(std::array<RI::Tensor<T>, N>& data,
+                            const int lmp0,
+                            const int lmq0,
+                            const std::array<RI::Tensor<T>, N>& val,
+                            const int lmp1,
+                            const int lmq1,
+                            const T& frac) {
+    for (size_t i = 0; i < N; ++i) {
+        data[i](lmp0, lmq0) += frac * val[i](lmp1, lmq1);
+    }
 }
 
+template <typename Tout, typename Tin>
+RI::Tensor<Tout> LRI_CV_Tools::convert(RI::Tensor<Tin>&& data) {
+    return RI::Global_Func::convert<Tout>(data);
+}
+
+template <typename Tout, typename Tin, std::size_t N>
+std::array<RI::Tensor<Tout>, N>
+    LRI_CV_Tools::convert(std::array<RI::Tensor<Tin>, N>&& data) {
+    std::array<RI::Tensor<Tout>, N> out;
+    for (size_t i = 0; i != N; ++i)
+        out[i] = RI::Global_Func::convert<Tout>(data[i]);
+    return out;
+}
+
+// template <typename T>
+// RI::Tensor<T> LRI_CV_Tools::check_zero(RI::Tensor<T>&& data) {
+//     RI::Tensor<T> result(data.shape);
+
+//     const std::size_t rows = data.shape[0];
+//     const std::size_t cols = data.shape[1];
+
+//     for (std::size_t i = 0; i < rows; ++i) {
+//         for (std::size_t j = 0; j < cols; ++j) {
+//             result(i, j) = LRI_CV_Tools::check_zero(data(i, j));
+//         }
+//     }
+
+//     return result;
+// }
+
+// template <typename T, std::size_t N>
+// std::array<RI::Tensor<T>, N>
+//     LRI_CV_Tools::check_zero(std::array<RI::Tensor<T>, N>&& data) {
+//     std::array<RI::Tensor<T>, N> result;
+
+//     for (size_t i = 0; i != N; ++i)
+//         result[i] = LRI_CV_Tools::check_zero(std::move(data[i]));
+
+//     return result;
+// }
+
 
 // dMRs[ipos0][ipos1] = \nabla_{ipos0} M R_{ipos1}
 template<typename TA, typename TC, typename Tdata>
@@ -361,4 +555,5 @@ LRI_CV_Tools::cal_dMRs(
 	}
 	return dMRs;
 }
+
 #endif
diff --git a/source/module_ri/RPA_LRI.h b/source/module_ri/RPA_LRI.h
index e194a8ea8d..ff6c05f0f8 100644
--- a/source/module_ri/RPA_LRI.h
+++ b/source/module_ri/RPA_LRI.h
@@ -69,7 +69,6 @@ template <typename T, typename Tdata> class RPA_LRI
 
     std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> lcaos;
     std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs;
-    std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_ccp;
 
     // Exx_LRI<double> exx_postSCF_double(info);
     // LRI_CV<Tdata> cv;
diff --git a/source/module_ri/RPA_LRI.hpp b/source/module_ri/RPA_LRI.hpp
index 79c851fa3a..567bd9d78e 100644
--- a/source/module_ri/RPA_LRI.hpp
+++ b/source/module_ri/RPA_LRI.hpp
@@ -23,7 +23,6 @@ void RPA_LRI<T, Tdata>::init(const MPI_Comm& mpi_comm_in, const K_Vectors& kv_in
     this->orb_cutoff_ = orb_cutoff;
     this->lcaos = exx_lri_rpa.lcaos;
     this->abfs = exx_lri_rpa.abfs;
-    this->abfs_ccp = exx_lri_rpa.abfs_ccp;
     this->p_kv = &kv_in;
 
     //	this->cv = std::move(exx_lri_rpa.cv);
@@ -44,7 +43,7 @@ void RPA_LRI<T, Tdata>::cal_rpa_cv(const UnitCell& ucell)
     const std::pair<std::vector<TA>, std::vector<std::vector<std::pair<TA, std::array<Tcell, Ndim>>>>> list_As_Vs
         = RI::Distribute_Equally::distribute_atoms(this->mpi_comm, atoms, period_Vs, 2, false);
 
-    std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs = exx_lri_rpa.cv.cal_Vs(ucell,
+    std::map<TA, std::map<TAC, RI::Tensor<Tdata>>> Vs = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.cal_Vs(ucell,
                                                                               list_As_Vs.first,
                                                                               list_As_Vs.second[0],
                                                                               {
@@ -57,8 +56,8 @@ void RPA_LRI<T, Tdata>::cal_rpa_cv(const UnitCell& ucell)
         = RI::Distribute_Equally::distribute_atoms_periods(this->mpi_comm, atoms, period_Cs, 2, false);
 
     std::pair<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>,
-              std::array<std::map<TA, std::map<TAC, RI::Tensor<Tdata>>>, 3>>
-        Cs_dCs = exx_lri_rpa.cv.cal_Cs_dCs(ucell,
+              std::map<TA, std::map<TAC, std::array<RI::Tensor<Tdata>, 3>>>>
+        Cs_dCs = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.cal_Cs_dCs(ucell,
                                            list_As_Cs.first,
                                            list_As_Cs.second[0],
                                            {
@@ -331,7 +330,7 @@ void RPA_LRI<T, Tdata>::out_coulomb_k(const UnitCell &ucell)
     for (int I = 0; I != ucell.nat; I++)
     {
         mu_shift[I] = all_mu;
-        all_mu += exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[I]);
+        all_mu += exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.get_index_abfs_size(ucell.iat2it[I]);
     }
     const int nks_tot = PARAM.inp.nspin == 2 ? (int)p_kv->get_nks() / 2 : p_kv->get_nks();
     std::stringstream ss;
@@ -344,7 +343,7 @@ void RPA_LRI<T, Tdata>::out_coulomb_k(const UnitCell &ucell)
     for (auto& Ip: this->Vs_period)
     {
         auto I = Ip.first;
-        size_t mu_num = exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[I]);
+        size_t mu_num = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.get_index_abfs_size(ucell.iat2it[I]);
 
         for (int ik = 0; ik != nks_tot; ik++)
         {
@@ -373,7 +372,7 @@ void RPA_LRI<T, Tdata>::out_coulomb_k(const UnitCell &ucell)
             {
                 auto iJ = vq_Jp.first;
                 auto& vq_J = vq_Jp.second;
-                size_t nu_num = exx_lri_rpa.cv.get_index_abfs_size(ucell.iat2it[iJ]);
+                size_t nu_num = exx_lri_rpa.exx_objs[Conv_Coulomb_Pot_K::Coulomb_Method::Center2].cv.get_index_abfs_size(ucell.iat2it[iJ]);
                 ofs << all_mu << "   " << mu_shift[I] + 1 << "   " << mu_shift[I] + mu_num << "  " << mu_shift[iJ] + 1
                     << "   " << mu_shift[iJ] + nu_num << std::endl;
                 ofs << ik + 1 << "  " << p_kv->wk[ik] / 2.0 * PARAM.inp.nspin << std::endl;
diff --git a/source/module_ri/conv_coulomb_pot_k.h b/source/module_ri/conv_coulomb_pot_k.h
index 7daca7e1d7..c9870e71c4 100644
--- a/source/module_ri/conv_coulomb_pot_k.h
+++ b/source/module_ri/conv_coulomb_pot_k.h
@@ -14,6 +14,7 @@ namespace Conv_Coulomb_Pot_K
 		Hf,						//		"hf_Rcut"
 		Erfc,					//		"hse_omega"
 		Erf};					//		"hse_omega", "hf_Rcut"
+	enum class Coulomb_Method{Center2, Ewald}; // Different methods for constructing the Coulomb matrix.
 
 	template<typename T> extern T cal_orbs_ccp(
 		const T &orbs,