diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
index 3128d47347..ee501561f7 100644
--- a/docs/advanced/input_files/input-main.md
+++ b/docs/advanced/input_files/input-main.md
@@ -3030,7 +3030,7 @@ These variables are relevant when using hybrid functionals with *[basis_type](#b
 ### exx_ccp_rmesh_times
 
 - **Type**: Real
-- **Description**: This parameter determines how many times larger the radial mesh required for calculating Columb potential is to that of atomic orbitals. The value should be at least 1. Reducing this value can effectively increase the speed of self-consistent calculations using hybrid functionals.
+- **Description**: This parameter determines how many times larger the radial mesh required for calculating Columb potential is to that of atomic orbitals. The value should be larger than 0. Reducing this value can effectively increase the speed of self-consistent calculations using hybrid functionals.
 - **Default**:
   - 5: if *[dft_functional](#dft_functional)==hf/pbe0/scan0/muller/power/wp22*
   - 1.5: if *[dft_functional](#dft_functional)==hse/cwp22*
diff --git a/source/source_base/sph_bessel_recursive-d2.cpp b/source/source_base/sph_bessel_recursive-d2.cpp
index 4418905d97..f23074e9c3 100644
--- a/source/source_base/sph_bessel_recursive-d2.cpp
+++ b/source/source_base/sph_bessel_recursive-d2.cpp
@@ -4,12 +4,12 @@
 //==========================================================
 
 #include "sph_bessel_recursive.h"
+#include "constants.h"
+#include "source_base/memory.h"
 
 #include<cmath>
 #include<stdexcept>
 
-#include "constants.h"
-
 namespace ModuleBase
 {
 
@@ -33,6 +33,7 @@ const std::vector<std::vector<std::vector<double>>> & Sph_Bessel_Recursive::D2::
 	cal_jlx_0( lmax+1, ix1_size, ix2_size );
 	cal_jlx_smallx( lmax+1, ix1_size, ix2_size );
 	cal_jlx_recursive( lmax+1, ix1_size, ix2_size );
+	ModuleBase::Memory::record("ORB::Jl(x)", sizeof(double) * (lmax+1) * ix1_size * ix2_size);
 	return jlx;
 }
 
@@ -41,7 +42,7 @@ void Sph_Bessel_Recursive::D2::cal_jlx_0( const int l_size, const size_t ix1_siz
 	if(jlx.size() < static_cast<size_t>(l_size))
 		jlx.resize(l_size);
 	
-	for( int l=0; l!=l_size; ++l )
+	for( int l=0; l<l_size; ++l )
 	{
 		const double jlx0 = (0==l) ? 1.0 : 0.0;
 
diff --git a/source/source_base/test/CMakeLists.txt b/source/source_base/test/CMakeLists.txt
index a39ba8436a..3481cba780 100644
--- a/source/source_base/test/CMakeLists.txt
+++ b/source/source_base/test/CMakeLists.txt
@@ -77,7 +77,7 @@ AddTest(
 AddTest(
   TARGET MODULE_BASE_sph_bessel_recursive
   LIBS parameter 
-  SOURCES sph_bessel_recursive_test.cpp ../sph_bessel_recursive-d1.cpp ../sph_bessel_recursive-d2.cpp
+  SOURCES sph_bessel_recursive_test.cpp ../sph_bessel_recursive-d1.cpp ../sph_bessel_recursive-d2.cpp ../memory.cpp ../global_variable.cpp
 )
 AddTest(
   TARGET MODULE_BASE_ylmreal
diff --git a/source/source_io/cal_r_overlap_R.cpp b/source/source_io/cal_r_overlap_R.cpp
index d947d792c8..f7512699e0 100644
--- a/source/source_io/cal_r_overlap_R.cpp
+++ b/source/source_io/cal_r_overlap_R.cpp
@@ -30,8 +30,8 @@ void cal_r_overlap_R::initialize_orb_table(const UnitCell& ucell,
     int Rmesh = static_cast<int>(orb.get_Rmax() / dr) + 4;
     Rmesh += 1 - Rmesh % 2;
 
-    int Lmax, Lmax_used;
-    std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_3(lmax_orb);
+    const int Lmax = lmax_orb + 1;
+    const int Lmax_used = 2 * lmax_orb + 1;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
                                              dk,
diff --git a/source/source_io/read_input_item_exx_dftu.cpp b/source/source_io/read_input_item_exx_dftu.cpp
index d235710a36..9889fd1156 100644
--- a/source/source_io/read_input_item_exx_dftu.cpp
+++ b/source/source_io/read_input_item_exx_dftu.cpp
@@ -332,9 +332,9 @@ void ReadInput::item_exx()
             }
         };
         item.check_value = [](const Input_Item& item, const Parameter& para) {
-            if (std::stod(para.input.exx_ccp_rmesh_times) < 1)
+            if (std::stod(para.input.exx_ccp_rmesh_times) <=0)
             {
-                ModuleBase::WARNING_QUIT("ReadInput", "exx_ccp_rmesh_times must >= 1");
+                ModuleBase::WARNING_QUIT("ReadInput", "exx_ccp_rmesh_times must > 0");
             }
         };
         this->add_item(item);
diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp
index b6898b57e4..60d44587cf 100644
--- a/source/source_io/test_serial/read_input_item_test.cpp
+++ b/source/source_io/test_serial/read_input_item_test.cpp
@@ -1397,7 +1397,7 @@ TEST_F(InputTest, Item_test2)
         it->second.reset_value(it->second, param);
         EXPECT_EQ(param.input.exx_ccp_rmesh_times, "1");
 
-        param.input.exx_ccp_rmesh_times = "0";
+        param.input.exx_ccp_rmesh_times = "-1";
         testing::internal::CaptureStdout();
         EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
         output = testing::internal::GetCapturedStdout();
diff --git a/source/source_io/to_wannier90_lcao.cpp b/source/source_io/to_wannier90_lcao.cpp
index 45414d15b9..a6bc0887b5 100644
--- a/source/source_io/to_wannier90_lcao.cpp
+++ b/source/source_io/to_wannier90_lcao.cpp
@@ -278,8 +278,8 @@ void toWannier90_LCAO::initialize_orb_table(const UnitCell& ucell)
     int Rmesh = static_cast<int>(orb_.get_Rmax() / dr) + 4;
     Rmesh += 1 - Rmesh % 2;
 
-    int Lmax, Lmax_used;
-    std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_3(lmax_orb);
+    const int Lmax = lmax_orb + 1;
+    const int Lmax_used = 2 * lmax_orb + 1;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
                                              dk,
diff --git a/source/source_io/unk_overlap_lcao.cpp b/source/source_io/unk_overlap_lcao.cpp
index 4d79d78a13..e9bf421bc4 100644
--- a/source/source_io/unk_overlap_lcao.cpp
+++ b/source/source_io/unk_overlap_lcao.cpp
@@ -40,8 +40,8 @@ void unkOverlap_lcao::init(const UnitCell& ucell,
     int Rmesh = static_cast<int>(orb.get_Rmax() / dr) + 4;
     Rmesh += 1 - Rmesh % 2;
 
-    int Lmax, Lmax_used;
-    std::tie(Lmax_used, Lmax) = Center2_Orb::init_Lmax_2_3(lmax_orb);
+    const int Lmax = lmax_orb + 1;
+    const int Lmax_used = 2 * lmax_orb + 1;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
                                              dk,
diff --git a/source/source_lcao/center2_orb-orb11.cpp b/source/source_lcao/center2_orb-orb11.cpp
index 9c4612d6fa..7e4304f01b 100644
--- a/source/source_lcao/center2_orb-orb11.cpp
+++ b/source/source_lcao/center2_orb-orb11.cpp
@@ -24,6 +24,8 @@ Center2_Orb::Orb11::Orb11(const Numerical_Orbital_Lm& nA_in,
 
 void Center2_Orb::Orb11::init_radial_table()
 {
+    const int rmesh = Center2_Orb::get_rmesh(this->nA.getRcut(), this->nB.getRcut(), dr_);
+
     const int LA = this->nA.getL();
     const int LB = this->nB.getL();
     for (int LAB = std::abs(LA - LB); LAB <= LA + LB; ++LAB)
@@ -33,8 +35,6 @@ void Center2_Orb::Orb11::init_radial_table()
             continue;
         }
 
-        const int rmesh = Center2_Orb::get_rmesh(this->nA.getRcut(), this->nB.getRcut(), dr_);
-
         this->Table_r[LAB].resize(rmesh, 0);
         this->Table_dr[LAB].resize(rmesh, 0);
 
@@ -43,8 +43,8 @@ void Center2_Orb::Orb11::init_radial_table()
                                     this->nA,
                                     this->nB,
                                     rmesh,
-                                    this->Table_r[LAB].data(),
-                                    this->Table_dr[LAB].data(),
+                                    this->Table_r[LAB],
+                                    this->Table_dr[LAB],
                                     psb_);
     }
     return;
@@ -81,8 +81,8 @@ void Center2_Orb::Orb11::init_radial_table(const std::set<size_t>& radials)
                                     this->nA,
                                     this->nB,
                                     radials_used,
-                                    this->Table_r[LAB].data(),
-                                    this->Table_dr[LAB].data(),
+                                    this->Table_r[LAB],
+                                    this->Table_dr[LAB],
                                     psb_);
     }
 }
diff --git a/source/source_lcao/center2_orb.cpp b/source/source_lcao/center2_orb.cpp
index b5ece1d3e3..b222716ce1 100644
--- a/source/source_lcao/center2_orb.cpp
+++ b/source/source_lcao/center2_orb.cpp
@@ -3,7 +3,6 @@
 #include "source_base/constants.h"
 #include "source_base/math_integral.h"
 #include "source_base/mathzone_add1.h"
-#include "source_base/memory.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_base/tool_title.h"
@@ -28,53 +27,6 @@ int Center2_Orb::get_rmesh(const double& R1, const double& R2, const double dr)
     return rmesh;
 }
 
-// used in <Phi|Phi> or <Beta|Phi>
-std::pair<int,int> Center2_Orb::init_Lmax_2_1(const int lmax_orb, const int lmax_beta)
-{
-    const int Lmax = std::max({-1, lmax_orb, lmax_beta});
-    const int Lmax_used = 2 * Lmax + 1;
-    assert(Lmax_used >= 1);
-    return {Lmax_used, Lmax};
-}
-
-// used in <jY|jY> or <Abfs|Abfs>
-std::pair<int,int> Center2_Orb::init_Lmax_2_2(const int& lmax_exx)
-{
-    const int Lmax = std::max(-1, lmax_exx);
-    const int Lmax_used = 2 * Lmax + 1;
-    assert(Lmax_used >= 1);
-    return {Lmax_used, Lmax};
-}
-
-// used in berryphase by jingan
-std::pair<int,int> Center2_Orb::init_Lmax_2_3(const int lmax_orb)
-{
-    const int Lmax = std::max(-1, lmax_orb) + 1;
-    const int Lmax_used = 2 * Lmax + 1;
-    assert(Lmax_used >= 1);
-    return {Lmax_used, Lmax};
-}
-
-// used in <jY|PhiPhi> or <Abfs|PhiPhi>
-std::pair<int,int> Center2_Orb::init_Lmax_3_1(const int& lmax_exx, const int lmax_orb)
-{
-    int Lmax = std::max(-1, lmax_orb);
-    int Lmax_used = 2 * Lmax + 1;
-    Lmax = std::max(Lmax, lmax_exx);
-    Lmax_used += lmax_exx;
-    assert(Lmax_used >= 1);
-    return {Lmax_used, Lmax};
-}
-
-// used in <PhiPhi|PhiPhi>
-std::pair<int,int> Center2_Orb::init_Lmax_4_1(const int lmax_orb)
-{
-    const int Lmax = std::max(-1, lmax_orb);
-    const int Lmax_used = 2 * (2 * Lmax + 1);
-    assert(Lmax_used >= 1);
-    return {Lmax_used, Lmax};
-}
-
 // Peize Lin update 2016-01-26
 void Center2_Orb::init_Table_Spherical_Bessel(const int Lmax_used,
                                               const double dr,
@@ -101,9 +53,7 @@ void Center2_Orb::init_Table_Spherical_Bessel(const int Lmax_used,
     }
 
     psb->set_dx(dr * dk);
-    psb->cal_jlx(Lmax_used, Rmesh, kmesh);
-
-    ModuleBase::Memory::record("ORB::Jl(x)", sizeof(double) * (Lmax_used + 1) * kmesh * Rmesh);
+    psb->cal_jlx(Lmax_used+1, Rmesh, kmesh);                // +1 for drs needs psb.jlx[l+1]. Peize Lin update 2025-12-27
 }
 
 // Peize Lin accelerate 2017-10-02
@@ -112,12 +62,15 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
                                  const Numerical_Orbital_Lm& n1,
                                  const Numerical_Orbital_Lm& n2,
                                  const int& rmesh,
-                                 double* rs,
-                                 double* drs,
+                                 std::vector<double> &rs,
+                                 std::vector<double> &drs,
                                  const ModuleBase::Sph_Bessel_Recursive::D2* psb)
 {
     ModuleBase::timer::tick("Center2_Orb", "cal_ST_Phi12_R");
 
+    assert(rmesh <= rs.size());
+    assert(rmesh <= drs.size());
+
     const int kmesh = n1.getNk();
     const double* kpoint = n1.getKpoint();
     const double dk = n1.getDk();
@@ -252,8 +205,8 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
                                  const Numerical_Orbital_Lm& n1,
                                  const Numerical_Orbital_Lm& n2,
                                  const std::set<size_t>& radials,
-                                 double* rs,
-                                 double* drs,
+                                 std::vector<double> &rs,
+                                 std::vector<double> &drs,
                                  const ModuleBase::Sph_Bessel_Recursive::D2* psb)
 {
     //	ModuleBase::TITLE("Center2_Orb","cal_ST_Phi12_R");
@@ -317,7 +270,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         // if(rs[ir])  => rs[ir]  has been calculated
         // if(drs[ir]) => drs[ir] has been calculated
         // Actually, if(ir[ir]||dr[ir]) is enough. Double insurance for the sake of avoiding numerical errors
-        if (rs[ir] && drs[ir]) {
+        if (rs.at(ir) && drs.at(ir)) {
             continue;
         }
 
@@ -330,7 +283,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         double temp = 0.0;
 
         ModuleBase::Integral::Simpson_Integral(kmesh, ModuleBase::GlobalFunc::VECTOR_TO_PTR(integrated_func), dk, temp);
-        rs[ir] = temp * ModuleBase::FOUR_PI;
+        rs.at(ir) = temp * ModuleBase::FOUR_PI;
 
         const std::vector<double>& jlm1_r = jlm1.at(ir);
         const std::vector<double>& jlp1_r = jlp1.at(ir);
@@ -353,7 +306,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
         }
 
         ModuleBase::Integral::Simpson_Integral(kmesh, ModuleBase::GlobalFunc::VECTOR_TO_PTR(integrated_func), dk, temp);
-        drs[ir] = -ModuleBase::FOUR_PI * (l + 1) / (2.0 * l + 1) * temp;
+        drs.at(ir) = -ModuleBase::FOUR_PI * (l + 1) / (2.0 * l + 1) * temp;
     }
 
     // cal rs[0] special
@@ -374,7 +327,7 @@ void Center2_Orb::cal_ST_Phi12_R(const int& job,
 
             // PLEASE try to make dualfac function as input parameters
             // mohan note 2021-03-23
-            rs[0] = ModuleBase::FOUR_PI / ModuleBase::Mathzone_Add1::dualfac(2 * l + 1) * temp;
+            rs.at(0) = ModuleBase::FOUR_PI / ModuleBase::Mathzone_Add1::dualfac(2 * l + 1) * temp;
         }
     }
 
diff --git a/source/source_lcao/center2_orb.h b/source/source_lcao/center2_orb.h
index 92288e97a8..f9a7688532 100644
--- a/source/source_lcao/center2_orb.h
+++ b/source/source_lcao/center2_orb.h
@@ -25,17 +25,6 @@ class Center2_Orb
 
     static int get_rmesh(const double& R1, const double& R2, const double dr);
 
-    // used in <Phi|Phi> or <Beta|Phi>
-    static std::pair<int,int> init_Lmax_2_1(const int lmax_orb, const int lmax_beta);
-    // used in <jY|jY> or <Abfs|Abfs>
-    static std::pair<int,int> init_Lmax_2_2(const int& lmax_exx);
-    // used in berryphase by jingan
-    static std::pair<int,int> init_Lmax_2_3(const int lmax_orb);
-    // used in <jY|PhiPhi> or <Abfs|PhiPhi>
-    static std::pair<int,int> init_Lmax_3_1(const int& lmax_exx, const int lmax_orb);
-    // used in <PhiPhi|PhiPhi>
-    static std::pair<int,int> init_Lmax_4_1(const int lmax_orb);
-
     static void init_Table_Spherical_Bessel(const int Lmax_used,
                                             const double dr,
                                             const double dk,
@@ -48,8 +37,8 @@ class Center2_Orb
                                const Numerical_Orbital_Lm& n1,
                                const Numerical_Orbital_Lm& n2,
                                const int& rmesh,
-                               double* rs,
-                               double* drs,
+                               std::vector<double> &rs,
+                               std::vector<double> &drs,
                                const ModuleBase::Sph_Bessel_Recursive::D2* psb);
 
     // Peize Lin add 2017-10-13
@@ -58,8 +47,8 @@ class Center2_Orb
                                const Numerical_Orbital_Lm& n1,
                                const Numerical_Orbital_Lm& n2,
                                const std::set<size_t>& radials, // only calculate ir in radials
-                               double* rs,
-                               double* drs,
+                               std::vector<double> &rs,
+                               std::vector<double> &drs,
                                const ModuleBase::Sph_Bessel_Recursive::D2* psb);
 };
 
diff --git a/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp b/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp
index f3efdfbd7e..51019df1b0 100644
--- a/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp
+++ b/source/source_lcao/module_ri/ABFs_Construct-PCA.cpp
@@ -173,7 +173,7 @@ RI::Tensor<double> get_column_mean0_matrix(const RI::Tensor<double>& m)
         = ModuleBase::Element_Basis_Index::construct_index(range_abfs);
 
 	Matrix_Orbs21 m_abfslcaos_lcaos;
-    m_abfslcaos_lcaos.init(abfs, lcaos, lcaos, ucell, orb, kmesh_times, orb.get_Rmax());
+    m_abfslcaos_lcaos.init(abfs, lcaos, lcaos, ucell, orb, kmesh_times);
 
     std::map<std::size_t, std::map<std::size_t, std::set<double>>> delta_R;
     for (std::size_t it = 0; it != abfs.size(); ++it)
diff --git a/source/source_lcao/module_ri/Exx_LRI_interface.hpp b/source/source_lcao/module_ri/Exx_LRI_interface.hpp
index 0494dce96b..06e0845da8 100644
--- a/source/source_lcao/module_ri/Exx_LRI_interface.hpp
+++ b/source/source_lcao/module_ri/Exx_LRI_interface.hpp
@@ -111,7 +111,7 @@ void Exx_LRI_Interface<T, Tdata>::exx_before_all_runners(
         this->symrot_.find_irreducible_sector(
             ucell.symm, ucell.atoms, ucell.st,
             RI_Util::get_Born_von_Karmen_cells(period), period, ucell.lat);
-        this->symrot_.set_abfs_Lmax(Exx_Abfs::get_Lmax(this->exx_ptr->abfs));
+        this->symrot_.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(this->exx_ptr->abfs));
         this->symrot_.cal_Ms(kv, ucell, pv);
     }
 }
diff --git a/source/source_lcao/module_ri/LRI_CV.hpp b/source/source_lcao/module_ri/LRI_CV.hpp
index 60670d5088..b343d4d20d 100644
--- a/source/source_lcao/module_ri/LRI_CV.hpp
+++ b/source/source_lcao/module_ri/LRI_CV.hpp
@@ -56,9 +56,6 @@ void LRI_CV<Tdata>::set_orbitals(
 
 	this->lcaos_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->lcaos);
     this->abfs_ccp_rcut = Exx_Abfs::Construct_Orbs::get_Rcut(this->abfs_ccp);
-    const double lcaos_rmax = Exx_Abfs::Construct_Orbs::get_Rmax(this->lcaos);
-    const double abfs_ccp_rmax
-        = Exx_Abfs::Construct_Orbs::get_Rmax(this->abfs_ccp);
 
 	const ModuleBase::Element_Basis_Index::Range
 		range_lcaos = ModuleBase::Element_Basis_Index::construct_range( lcaos );
@@ -71,11 +68,11 @@ void LRI_CV<Tdata>::set_orbitals(
 	this->m_abfs_abfs.MGT = this->m_abfslcaos_lcaos.MGT = MGT;
     this->m_abfs_abfs.init(
 		this->abfs_ccp, this->abfs,
-		ucell, orb, kmesh_times, lcaos_rmax + abfs_ccp_rmax);
+		ucell, orb, kmesh_times);
     if (init_C)
         this->m_abfslcaos_lcaos.init(
 			this->abfs_ccp, this->lcaos, this->lcaos,
-			ucell, orb, kmesh_times, lcaos_rmax);
+			ucell, orb, kmesh_times);
 
     this->m_abfs_abfs.init_radial_table();
     if (init_C) {
diff --git a/source/source_lcao/module_ri/Matrix_Orbs11.cpp b/source/source_lcao/module_ri/Matrix_Orbs11.cpp
index 942f5a4cc6..31d7a87dce 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs11.cpp
+++ b/source/source_lcao/module_ri/Matrix_Orbs11.cpp
@@ -5,6 +5,7 @@
 
 #include "Matrix_Orbs11.h"
 
+#include "exx_abfs-construct_orbs.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
 #include "source_pw/module_pwdft/global.h"
@@ -14,16 +15,15 @@ void Matrix_Orbs11::init(
     const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B, 
     const UnitCell& ucell,
     const LCAO_Orbitals& orb, 
-    const double kmesh_times, 
-    const double rmax)
+    const double kmesh_times)
 {
     ModuleBase::TITLE("Matrix_Orbs11", "init");
     ModuleBase::timer::tick("Matrix_Orbs11", "init");
 
     this->lat0 = &ucell.lat0;
 
-    const int Lmax = std::max({ Exx_Abfs::get_Lmax(orb_A), Exx_Abfs::get_Lmax(orb_B) });
-    const int Lmax_used = Exx_Abfs::get_Lmax(orb_A) + Exx_Abfs::get_Lmax(orb_B);
+    const int Lmax = std::max({ Exx_Abfs::Construct_Orbs::get_Lmax(orb_A), Exx_Abfs::Construct_Orbs::get_Lmax(orb_B) });
+    const int Lmax_used = Exx_Abfs::Construct_Orbs::get_Lmax(orb_A) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_B);
 
     //=========================================
     // (3) make Gaunt coefficients table
@@ -38,7 +38,10 @@ void Matrix_Orbs11::init(
     const double dr = orb.get_dR();
     const double dk = orb.get_dk();
     const int kmesh = orb.get_kmesh() * kmesh_times + 1;
-    int Rmesh = static_cast<int>(rmax / dr) + 4;
+    const double rmax
+        = Exx_Abfs::Construct_Orbs::get_Rmax(orb_A)
+        + Exx_Abfs::Construct_Orbs::get_Rmax(orb_B);
+    int Rmesh = static_cast<int>(rmax / dr) + 4;                            // extend Rcut, keep dR
     Rmesh += 1 - Rmesh % 2;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
diff --git a/source/source_lcao/module_ri/Matrix_Orbs11.h b/source/source_lcao/module_ri/Matrix_Orbs11.h
index 069127c355..7b520fe0f6 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs11.h
+++ b/source/source_lcao/module_ri/Matrix_Orbs11.h
@@ -26,8 +26,7 @@ class Matrix_Orbs11
         const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
         const UnitCell& ucell,
         const LCAO_Orbitals& orb,
-        const double kmesh_times,       // extend Kcut, keep dK
-        const double rmax);             // extend Rcut, keep dR
+        const double kmesh_times);       // extend Kcut, keep dK
 
     void init_radial_table();
     void init_radial_table(const std::map<size_t, std::map<size_t, std::set<double>>>& Rs); // unit: ucell.lat0
diff --git a/source/source_lcao/module_ri/Matrix_Orbs21.cpp b/source/source_lcao/module_ri/Matrix_Orbs21.cpp
index a11a1333b5..9317e9dbd8 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs21.cpp
+++ b/source/source_lcao/module_ri/Matrix_Orbs21.cpp
@@ -5,6 +5,7 @@
 
 #include "Matrix_Orbs21.h"
 
+#include "exx_abfs-construct_orbs.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
 #include "source_pw/module_pwdft/global.h"
@@ -15,17 +16,16 @@ void Matrix_Orbs21::init(
     const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
     const UnitCell& ucell,
     const LCAO_Orbitals& orb, 
-    const double kmesh_times, 
-    const double rmax)
+    const double kmesh_times)
 {
     ModuleBase::TITLE("Matrix_Orbs21", "init");
     ModuleBase::timer::tick("Matrix_Orbs21", "init");
     this->lat0 = &ucell.lat0;
 
     const int Lmax = std::max({
-        Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2),
-        Exx_Abfs::get_Lmax(orb_B) });
-    const int Lmax_used = Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2) + Exx_Abfs::get_Lmax(orb_B) + 1;
+        Exx_Abfs::Construct_Orbs::get_Lmax(orb_A1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_A2),
+        Exx_Abfs::Construct_Orbs::get_Lmax(orb_B) });
+    const int Lmax_used = Exx_Abfs::Construct_Orbs::get_Lmax(orb_A1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_A2) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_B);
 
     //=========================================
     // (3) make Gaunt coefficients table
@@ -36,11 +36,14 @@ void Matrix_Orbs21::init(
         { this->MGT->init_Gaunt_CH(Lmax); }
     if(this->MGT->get_Lmax_Gaunt_Coefficients() < Lmax)
         { this->MGT->init_Gaunt(Lmax); }
-
+    
     const double dr = orb.get_dR();
     const double dk = orb.get_dk();
     const int kmesh = orb.get_kmesh() * kmesh_times + 1;
-    int Rmesh = static_cast<int>(rmax / dr) + 4;
+    const double rmax
+        = std::min({Exx_Abfs::Construct_Orbs::get_Rmax(orb_A1), Exx_Abfs::Construct_Orbs::get_Rmax(orb_A2)})
+        + Exx_Abfs::Construct_Orbs::get_Rmax(orb_B);
+   int Rmesh = static_cast<int>(rmax / dr) + 4;                            // extend Rcut, keep dR
     Rmesh += 1 - Rmesh % 2;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
diff --git a/source/source_lcao/module_ri/Matrix_Orbs21.h b/source/source_lcao/module_ri/Matrix_Orbs21.h
index d59ef9672e..efd561a00f 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs21.h
+++ b/source/source_lcao/module_ri/Matrix_Orbs21.h
@@ -26,8 +26,7 @@ class Matrix_Orbs21
         const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B,
         const UnitCell& ucell,
         const LCAO_Orbitals& orb, 
-        const double kmesh_times,       // extend Kcut, keep dK
-        const double rmax);             // extend Rcut, keep dR
+        const double kmesh_times);       // extend Kcut, keep dK
 
     void init_radial_table();
     void init_radial_table(const std::map<size_t, std::map<size_t, std::set<double>>>& Rs); // unit: ucell.lat0
diff --git a/source/source_lcao/module_ri/Matrix_Orbs22.cpp b/source/source_lcao/module_ri/Matrix_Orbs22.cpp
index 69ea735fbc..1dca7a351b 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs22.cpp
+++ b/source/source_lcao/module_ri/Matrix_Orbs22.cpp
@@ -5,6 +5,7 @@
 
 #include "Matrix_Orbs22.h"
 
+#include "exx_abfs-construct_orbs.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
 #include "source_pw/module_pwdft/global.h"
@@ -16,8 +17,7 @@ void Matrix_Orbs22::init(
     const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B2,
     const UnitCell& ucell,
     const LCAO_Orbitals& orb, 
-    const double kmesh_times, 
-    const double rmax)
+    const double kmesh_times)
 {
     ModuleBase::TITLE("Matrix_Orbs22", "init");
     ModuleBase::timer::tick("Matrix_Orbs22", "init");
@@ -25,9 +25,9 @@ void Matrix_Orbs22::init(
     this->lat0 = &ucell.lat0;
 
     const int Lmax = std::max({
-        Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2),
-        Exx_Abfs::get_Lmax(orb_B1) + Exx_Abfs::get_Lmax(orb_B2) });
-    const int Lmax_used = Exx_Abfs::get_Lmax(orb_A1) + Exx_Abfs::get_Lmax(orb_A2) + Exx_Abfs::get_Lmax(orb_B1) + Exx_Abfs::get_Lmax(orb_B2) + 2;
+        Exx_Abfs::Construct_Orbs::get_Lmax(orb_A1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_A2),
+        Exx_Abfs::Construct_Orbs::get_Lmax(orb_B1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_B2) });
+    const int Lmax_used = Exx_Abfs::Construct_Orbs::get_Lmax(orb_A1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_A2) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_B1) + Exx_Abfs::Construct_Orbs::get_Lmax(orb_B2);
 
     //=========================================
     // (3) make Gaunt coefficients table
@@ -42,7 +42,10 @@ void Matrix_Orbs22::init(
     const double dr = orb.get_dR();
     const double dk = orb.get_dk();
     const int kmesh = orb.get_kmesh() * kmesh_times + 1;
-    int Rmesh = static_cast<int>(rmax / dr) + 4;
+    const double rmax
+        = std::min({Exx_Abfs::Construct_Orbs::get_Rmax(orb_A1), Exx_Abfs::Construct_Orbs::get_Rmax(orb_A2)})
+        + std::min({Exx_Abfs::Construct_Orbs::get_Rmax(orb_B1), Exx_Abfs::Construct_Orbs::get_Rmax(orb_B2)});
+     int Rmesh = static_cast<int>(rmax / dr) + 4;                            // extend Rcut, keep dR
     Rmesh += 1 - Rmesh % 2;
     Center2_Orb::init_Table_Spherical_Bessel(Lmax_used,
                                              dr,
diff --git a/source/source_lcao/module_ri/Matrix_Orbs22.h b/source/source_lcao/module_ri/Matrix_Orbs22.h
index e0c3790289..7451f80ed6 100644
--- a/source/source_lcao/module_ri/Matrix_Orbs22.h
+++ b/source/source_lcao/module_ri/Matrix_Orbs22.h
@@ -28,8 +28,7 @@ class Matrix_Orbs22
         const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb_B2,
         const UnitCell& ucell,
         const LCAO_Orbitals& orb,
-        const double kmesh_times,  // extend Kcut, keep dK
-        const double rmax);        // extend Rcut, keep dR
+        const double kmesh_times);  // extend Kcut, keep dK
 
     void init_radial_table();
     void init_radial_table(const std::map<size_t, std::map<size_t, std::set<double>>>& Rs); // unit: ucell.lat0
diff --git a/source/source_lcao/module_ri/RPA_LRI.hpp b/source/source_lcao/module_ri/RPA_LRI.hpp
index 641e07a7f4..8b039122ee 100644
--- a/source/source_lcao/module_ri/RPA_LRI.hpp
+++ b/source/source_lcao/module_ri/RPA_LRI.hpp
@@ -95,7 +95,7 @@ void RPA_LRI<T, Tdata>::cal_postSCF_exx(const elecstate::DensityMatrix<T, Tdata>
         symrot.find_irreducible_sector(ucell.symm, ucell.atoms, ucell.st,
             Rs, period, ucell.lat);
         // set Lmax of the rotation matrices to max(l_ao, l_abf), to support rotation under ABF
-        symrot.set_abfs_Lmax(Exx_Abfs::get_Lmax(exx_lri_rpa.abfs));
+        symrot.set_abfs_Lmax(Exx_Abfs::Construct_Orbs::get_Lmax(exx_lri_rpa.abfs));
         symrot.cal_Ms(kv, ucell, *dm.get_paraV_pointer());
         // output Ts (symrot_R.txt) and Ms (symrot_k.txt)
         ModuleSymmetry::print_symrot_info_R(symrot, ucell.symm, ucell.lmax, Rs);
diff --git a/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp b/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp
index 1556335572..0a0e22b21a 100644
--- a/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp
+++ b/source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp
@@ -86,6 +86,8 @@ std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> Exx_Abfs::Construct_
 	const double times_threshold )
 {
 	ModuleBase::TITLE("Exx_Abfs::Construct_Orbs::abfs_same_atom");
+	if(times_threshold>1)
+		{ return std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>(orb.get_ntype()); }
 
 	const std::vector<std::vector<std::vector<std::vector<double>>>>
 		abfs_same_atom_psi = psi_mult_psi( orbs );
diff --git a/source/source_lcao/module_ri/exx_abfs-construct_orbs.h b/source/source_lcao/module_ri/exx_abfs-construct_orbs.h
index ae698b33bd..01b2dc36a9 100644
--- a/source/source_lcao/module_ri/exx_abfs-construct_orbs.h
+++ b/source/source_lcao/module_ri/exx_abfs-construct_orbs.h
@@ -22,7 +22,7 @@ class Exx_Abfs::Construct_Orbs
 
 	static std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs_same_atom( 
 		const UnitCell &ucell,
-        const LCAO_Orbitals& orb,
+		const LCAO_Orbitals& orb,
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
 		const double kmesh_times_mot,
 		const double times_threshold=0);
@@ -32,30 +32,36 @@ class Exx_Abfs::Construct_Orbs
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs,
 		std::ostream &os);		
 
-    // get the max number of orbitals among all elements
-    // static int get_nmax_total(const
-    // std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb_in); get
-    // number of orbitals for each element static std::map<int, int>
-    // get_nw(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
-    // &orb_in);
+	// get the max number of orbitals among all elements
+	// static int get_nmax_total(const
+	// std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb_in); get
+	// number of orbitals for each element static std::map<int, int>
+	// get_nw(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>
+	// &orb_in);
 
-    // get multipole of orbitals for each element and angular moment
-    static std::vector<std::vector<std::vector<double>>> get_multipole(
-        const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>&
-            orb_in);
+	// get multipole of orbitals for each element and angular moment
+	static std::vector<std::vector<std::vector<double>>> get_multipole(
+		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb_in);
 
-    static std::vector<double> get_Rcut(
-        const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>&
-            orb_in);
-    static inline double get_Rmax(const std::vector<double>& rcut) {
-        return *std::max_element(rcut.begin(), rcut.end());
-    }
-    static inline double get_Rmax(
-        const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>&
-            orb_in) {
-        std::vector<double> rcut = get_Rcut(orb_in);
-        return get_Rmax(rcut);
-    }
+	static std::vector<double> get_Rcut(
+		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb_in);
+	static inline double get_Rmax(const std::vector<double>& rcut)
+	{
+		return *std::max_element(rcut.begin(), rcut.end());
+	}
+	static inline double get_Rmax(
+		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orb_in)
+	{
+		std::vector<double> rcut = get_Rcut(orb_in);
+		return get_Rmax(rcut);
+	}
+	template<typename T>
+	static int get_Lmax(const std::vector<std::vector<T>> &orb)
+	{
+		return	max_element(orb.begin(), orb.end(),
+				            [](const std::vector<T> &orb_A, const std::vector<T> &orb_B){ return orb_A.size() < orb_B.size(); })
+				->size() - 1;
+	}
 
 	static void filter_empty_orbs(
 		std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs);
@@ -74,7 +80,7 @@ class Exx_Abfs::Construct_Orbs
 
 	static std::vector<std::vector<std::vector<std::vector<double>>>> pca(
 		const UnitCell &ucell,
-        const LCAO_Orbitals& orb,
+		const LCAO_Orbitals& orb,
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs,
 		const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs,
 		const double kmesh_times_mot,
diff --git a/source/source_lcao/module_ri/exx_abfs-io.cpp b/source/source_lcao/module_ri/exx_abfs-io.cpp
index 3a8de0f1b3..2d78dcb4d8 100644
--- a/source/source_lcao/module_ri/exx_abfs-io.cpp
+++ b/source/source_lcao/module_ri/exx_abfs-io.cpp
@@ -15,6 +15,7 @@ std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> Exx_Abfs::IO::constr
 	const std::vector<std::string> &files_abfs,
 	const double kmesh_times )
 {
+	ModuleBase::TITLE("Exx_Abfs::IO::construct_abfs");
 	std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> abfs( files_abfs.size() );
 	for( size_t T=0; T!=files_abfs.size(); ++T )
 		abfs[T] = construct_abfs_T( 
diff --git a/source/source_lcao/module_ri/exx_abfs.h b/source/source_lcao/module_ri/exx_abfs.h
index 42cc41a01d..8a78307b36 100644
--- a/source/source_lcao/module_ri/exx_abfs.h
+++ b/source/source_lcao/module_ri/exx_abfs.h
@@ -22,14 +22,6 @@ class Exx_Abfs
 	
 	int rmesh_times = 5;				// Peize Lin test
 	int kmesh_times = 1;				// Peize Lin test
-
-	static int get_Lmax(const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>& orb)
-	{
-		int Lmax = -1;
-		for( const auto &orb_T : orb )
-			{ Lmax = std::max( Lmax, static_cast<int>(orb_T.size())-1 ); }
-		return Lmax;
-	}
 };
 
 #endif
diff --git a/source/source_lcao/module_ri/exx_opt_orb.cpp b/source/source_lcao/module_ri/exx_opt_orb.cpp
index 83ea26c18d..a00af11834 100644
--- a/source/source_lcao/module_ri/exx_opt_orb.cpp
+++ b/source/source_lcao/module_ri/exx_opt_orb.cpp
@@ -66,7 +66,7 @@ void Exx_Opt_Orb::generate_matrix(
 	{
 		if(judge_orbs_empty(lcaos))	{ return {}; }
 		Matrix_Orbs22 m_lcaoslcaos_lcaoslcaos;
-		m_lcaoslcaos_lcaoslcaos.init( lcaos, lcaos, lcaos, lcaos, ucell,orb, info.kmesh_times, orb.get_Rmax() );
+		m_lcaoslcaos_lcaoslcaos.init( lcaos, lcaos, lcaos, lcaos, ucell,orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_lcaoslcaos_lcaoslcaos.init_radial_table(radial_R);
 		#else
@@ -81,7 +81,7 @@ void Exx_Opt_Orb::generate_matrix(
 		if(judge_orbs_empty(lcaos))	{ return {}; }
 		if(judge_orbs_empty(jle))	{ return {}; }
 		Matrix_Orbs21 m_jyslcaos_lcaos;
-		m_jyslcaos_lcaos.init( jle, lcaos, lcaos, ucell , orb, info.kmesh_times, orb.get_Rmax() );
+		m_jyslcaos_lcaos.init( jle, lcaos, lcaos, ucell , orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_jyslcaos_lcaos.init_radial_table( radial_R);
 		#else
@@ -95,7 +95,7 @@ void Exx_Opt_Orb::generate_matrix(
 	{
 		if(judge_orbs_empty(jle))	{ return {}; }
 		Matrix_Orbs11 m_jys_jys;
-		m_jys_jys.init( jle, jle, ucell,orb, info.kmesh_times, orb.get_Rmax() );
+		m_jys_jys.init( jle, jle, ucell,orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_jys_jys.init_radial_table(radial_R);
 		#else
@@ -109,7 +109,7 @@ void Exx_Opt_Orb::generate_matrix(
 	{
 		if(judge_orbs_empty(abfs))	{ return {}; }
 		Matrix_Orbs11 m_abfs_abfs;
-		m_abfs_abfs.init( abfs, abfs, ucell, orb, info.kmesh_times, orb.get_Rmax() );
+		m_abfs_abfs.init( abfs, abfs, ucell, orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_abfs_abfs.init_radial_table(radial_R);
 		#else
@@ -124,7 +124,7 @@ void Exx_Opt_Orb::generate_matrix(
 		if(judge_orbs_empty(lcaos))	{ return {}; }
 		if(judge_orbs_empty(abfs))	{ return {}; }
 		Matrix_Orbs21 m_abfslcaos_lcaos;
-		m_abfslcaos_lcaos.init( abfs, lcaos, lcaos, ucell , orb, info.kmesh_times, orb.get_Rmax() );
+		m_abfslcaos_lcaos.init( abfs, lcaos, lcaos, ucell , orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_abfslcaos_lcaos.init_radial_table(radial_R);
 		#else
@@ -139,7 +139,7 @@ void Exx_Opt_Orb::generate_matrix(
 		if(judge_orbs_empty(jle))	{ return {}; }
 		if(judge_orbs_empty(abfs))	{ return {}; }
 		Matrix_Orbs11 m_jys_abfs;
-		m_jys_abfs.init( jle, abfs, ucell,orb, info.kmesh_times, orb.get_Rmax() );
+		m_jys_abfs.init( jle, abfs, ucell,orb, info.kmesh_times );
 		#if TEST_EXX_RADIAL>=1
 		m_jys_abfs.init_radial_table(radial_R);
 		#else