Simplify Lmax and rmesh in Matrix_Orbs and Center2_Orb

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*

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;
 

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

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,

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);

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();

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,

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,

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_);
     }
 }

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

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);
 };
 

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)

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);
     }
 }

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) {