From 06ab2be77f993b045a623948661663c86225c99b Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Sun, 1 Nov 2026 23:35:51 +0800 Subject: [PATCH 1/2] perf(lcao): reuse common terms across magnetic components --- source/source_lcao/center2orb_orb11.cpp | 216 ++++++++++++--- source/source_lcao/center2orb_orb11.h | 7 + source/source_lcao/center2orb_orb21.cpp | 58 ++++ source/source_lcao/center2orb_orb21.h | 11 + .../source_lcao/module_ri/matrix_orbs21.hpp | 53 ++-- source/source_lcao/test/CMakeLists.txt | 15 + .../test/test_center2orb_orb21.cpp | 261 ++++++++++++++++++ 7 files changed, 566 insertions(+), 55 deletions(-) create mode 100644 source/source_lcao/test/test_center2orb_orb21.cpp diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index 7241aaba141..bb5aa5635e1 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -119,9 +119,19 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, const int idx2 = this->MGT.get_lm_index(LB, mB); const double* Gaunt_Coefficients_ptr = &(this->MGT.Gaunt_Coefficients(idx1, idx2, 0)); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); for (const auto& tb_r: this->Table_r) { const int LAB = tb_r.first; + const double distance_power = std::pow(distance, LAB); + const int phase = (LA - LB - LAB) / 2; + const double i_exp = std::pow(-1.0, phase); + const double Interp_Tlm = (distance > tiny2) ? ModuleBase::PolyInt::Polynomial_Interpolation( + tb_r.second.data(), + rmesh, + dr_, + distance) + : tb_r.second.at(0); for (int mAB = 0; mAB != 2 * LAB + 1; ++mAB) // const int mAB = mA + mB; @@ -138,17 +148,7 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, { continue; } - const double ylm_real = (distance > tiny2) ? ylm_solid / pow(distance, LAB) : ylm_solid; - - const double i_exp = std::pow(-1.0, (LA - LB - LAB) / 2); - - const double Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) - : tb_r.second.at(0); + const double ylm_real = (distance > tiny2) ? ylm_solid / distance_power : ylm_solid; overlap += // pow(2*PI,1.5) @@ -197,9 +197,26 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 ModuleBase::Vector3 grad_overlap(0.0, 0.0, 0.0); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); for (const auto& tb_r: this->Table_r) { const int LAB = tb_r.first; + const double distance_power = std::pow(distance, LAB); + const int phase = (LA - LB - LAB) / 2; + const double i_exp = std::pow(-1.0, phase); + const double Interp_Tlm = (distance > tiny2) ? ModuleBase::PolyInt::Polynomial_Interpolation( + tb_r.second.data(), + rmesh, + dr_, + distance) + : tb_r.second.at(0); + const double grad_Interp_Tlm = (distance > tiny2) ? ModuleBase::PolyInt::Polynomial_Interpolation( + this->Table_dr.at(LAB).data(), + rmesh, + dr_, + distance) - + Interp_Tlm * LAB / distance + : 0.0; for (int mAB = 0; mAB != 2 * LAB + 1; ++mAB) // const int mAB = mA + mB; { @@ -211,30 +228,10 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 } const double ylm_solid = rly[idx3]; - const double ylm_real = (distance > tiny2) ? ylm_solid / pow(distance, LAB) : ylm_solid; + const double ylm_real = (distance > tiny2) ? ylm_solid / distance_power : ylm_solid; const ModuleBase::Vector3 gylm_solid = grly[idx3]; - const ModuleBase::Vector3 gylm_real - = (distance > tiny2) ? gylm_solid / pow(distance, LAB) : gylm_solid; - - const double i_exp = std::pow(-1.0, (LA - LB - LAB) / 2); - - const double Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(tb_r.second.data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) - : tb_r.second.at(0); - - const double grad_Interp_Tlm - = (distance > tiny2) - ? ModuleBase::PolyInt::Polynomial_Interpolation(this->Table_dr.at(LAB).data(), - Center2_Orb::get_rmesh(RcutA, RcutB, dr_), - dr_, - distance) // Interp(Table_dr) - - Interp_Tlm * LAB / distance - : 0.0; + const ModuleBase::Vector3 gylm_real = (distance > tiny2) ? gylm_solid / distance_power : gylm_solid; grad_overlap += i_exp // pow(2*PI,1.5) * Gaunt_real_A_B_AB @@ -244,3 +241,156 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 return grad_overlap; } + +template std::vector Center2_Orb::Orb11::contract_angular(const std::vector &components) const +{ + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const int count_A = 2 * LA + 1; + const int count_B = 2 * LB + 1; + const int count = count_A * count_B; + std::vector result(count); + for (int ma = 0; ma < count_A; ++ma) + { + const int idx_A = this->MGT.get_lm_index(LA, ma); + for (int mb = 0; mb < count_B; ++mb) + { + const int idx_B = this->MGT.get_lm_index(LB, mb); + const double *coefficients = &this->MGT.Gaunt_Coefficients(idx_A, idx_B, 0); + T &value = result[ma * count_B + mb]; + // Keep the scalar path's angular-momentum accumulation order. + for (const std::pair> &table : this->Table_r) + { + const int l = table.first; + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + if (coefficients[idx] != 0.0) + { + value += coefficients[idx] * components[idx]; + } + } + } + } + } + return result; +} + +std::vector Center2_Orb::Orb11::cal_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const +{ + const double tiny1 = 1e-12; + const double tiny2 = 1e-10; + const ModuleBase::Vector3 delta_R = RB - RA; + const double distance_true = delta_R.norm(); + const double distance = distance_true >= tiny1 ? distance_true : distance_true + tiny1; + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const double RcutA = this->nA.getRcut(); + const double RcutB = this->nB.getRcut(); + if (distance > RcutA + RcutB) + { + const int count = (2 * LA + 1) * (2 * LB + 1); + return std::vector(count); + } + + const int lmax = LA + LB; + std::vector rly; + ModuleBase::Ylm::rl_sph_harm( + lmax, + delta_R.x, + delta_R.y, + delta_R.z, + rly); + const int angular_count = (lmax + 1) * (lmax + 1); + std::vector components(angular_count); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_); + for (const std::pair> &table : this->Table_r) + { + const int l = table.first; + const double distance_power = std::pow(distance, l); + const int phase = (LA - LB - l) / 2; + const double i_exp = std::pow(-1.0, phase); + const double radial = distance > tiny2 ? ModuleBase::PolyInt::Polynomial_Interpolation( + table.second.data(), + rmesh, + dr_, + distance) + : table.second.at(0); + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + const double angular = distance > tiny2 ? rly[idx] / distance_power : rly[idx]; + components[idx] = i_exp * radial * angular; + } + } + return this->contract_angular(components); +} + +std::vector> Center2_Orb::Orb11::cal_grad_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const +{ + const double tiny1 = 1e-12; + const double tiny2 = 1e-10; + const ModuleBase::Vector3 delta_R = RB - RA; + const double distance_true = delta_R.norm(); + const double distance = distance_true >= tiny1 ? distance_true : distance_true + tiny1; + const int LA = this->nA.getL(); + const int LB = this->nB.getL(); + const double RcutA = this->nA.getRcut(); + const double RcutB = this->nB.getRcut(); + if (distance > RcutA + RcutB) + { + const int count = (2 * LA + 1) * (2 * LB + 1); + return std::vector>(count); + } + + const int lmax = LA + LB; + const int angular_count = (lmax + 1) * (lmax + 1); + const int gradient_count = 3 * angular_count; + std::vector rly(angular_count); + std::vector grly(gradient_count); + ModuleBase::Ylm::grad_rl_sph_harm( + lmax, + delta_R.x, + delta_R.y, + delta_R.z, + rly.data(), + grly.data()); + std::vector> components(angular_count); + const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_); + for (const std::pair> &table : this->Table_r) + { + const int l = table.first; + const double distance_power = std::pow(distance, l); + const int phase = (LA - LB - l) / 2; + const double i_exp = std::pow(-1.0, phase); + const double radial = distance > tiny2 ? ModuleBase::PolyInt::Polynomial_Interpolation( + table.second.data(), + rmesh, + dr_, + distance) + : table.second.at(0); + const double derivative = distance > tiny2 ? ModuleBase::PolyInt::Polynomial_Interpolation( + this->Table_dr.at(l).data(), + rmesh, + dr_, + distance) - + radial * l / distance + : 0.0; + for (int m = 0; m <= 2 * l; ++m) + { + const int idx = this->MGT.get_lm_index(l, m); + const double angular = distance > tiny2 ? rly[idx] / distance_power : rly[idx]; + const int offset = 3 * idx; + const double gx = grly[offset]; + const double gy = grly[offset + 1]; + const double gz = grly[offset + 2]; + const ModuleBase::Vector3 solid_gradient(gx, gy, gz); + const ModuleBase::Vector3 angular_gradient = + distance > tiny2 ? solid_gradient / distance_power : solid_gradient; + components[idx] = i_exp * (radial * angular_gradient + derivative * angular * delta_R / distance); + } + } + return this->contract_angular(components); +} diff --git a/source/source_lcao/center2orb_orb11.h b/source/source_lcao/center2orb_orb11.h index a7871eb1728..bbc64b7be83 100644 --- a/source/source_lcao/center2orb_orb11.h +++ b/source/source_lcao/center2orb_orb11.h @@ -38,7 +38,14 @@ class Center2_Orb::Orb11 const int& mA, const int& mB) const; + /// All magnetic components, ordered as [mA][mB]. + std::vector cal_overlap(const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const; + std::vector> cal_grad_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const; + private: + template std::vector contract_angular(const std::vector &components) const; + const Numerical_Orbital_Lm& nA; const Numerical_Orbital_Lm& nB; diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 261c7402526..8d09a300aae 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -169,3 +169,61 @@ ModuleBase::Vector3 Center2_Orb::Orb21::cal_grad_overlap(const ModuleBas return grad_overlap; } + +template +std::vector Center2_Orb::Orb21::cal_magnetic_components( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB, + std::vector (Orb11::*evaluate)(const ModuleBase::Vector3 &, const ModuleBase::Vector3 &) + const) const +{ + const int LA1 = this->nA1.getL(); + const int LA2 = this->nA2.getL(); + const int count_A1 = 2 * LA1 + 1; + const int count_A2 = 2 * LA2 + 1; + const int count_B = 2 * this->nB.getL() + 1; + const int count = count_A1 * count_A2 * count_B; + std::vector result(count); + for (const auto &orbital : this->orb11s) + { + const int LA = orbital.first; + const int count_A = 2 * LA + 1; + // These two-center overlaps are shared by all ways to couple the + // magnetic components of the two original orbitals on center A. + const std::vector overlaps = (orbital.second.*evaluate)(RA, RB); + for (int ma1 = 0; ma1 < count_A1; ++ma1) + { + const int idx1 = this->MGT.get_lm_index(LA1, ma1); + for (int ma2 = 0; ma2 < count_A2; ++ma2) + { + const int idx2 = this->MGT.get_lm_index(LA2, ma2); + const int offset = (ma1 * count_A2 + ma2) * count_B; + for (int ma = 0; ma < count_A; ++ma) + { + const int idx = this->MGT.get_lm_index(LA, ma); + const double coefficient = this->MGT.Gaunt_Coefficients(idx1, idx2, idx); + if (coefficient == 0.0) + { + continue; + } + for (int mb = 0; mb < count_B; ++mb) + { + result[offset + mb] += coefficient * overlaps[ma * count_B + mb]; + } + } + } + } + } + return result; +} + +std::vector Center2_Orb::Orb21::cal_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const +{ + return this->cal_magnetic_components(RA, RB, &Orb11::cal_overlap); +} + +std::vector> Center2_Orb::Orb21::cal_grad_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const +{ + return this->cal_magnetic_components>(RA, RB, &Orb11::cal_grad_overlap); +} diff --git a/source/source_lcao/center2orb_orb21.h b/source/source_lcao/center2orb_orb21.h index bc7cb5c7fa1..3170bfdf5f2 100644 --- a/source/source_lcao/center2orb_orb21.h +++ b/source/source_lcao/center2orb_orb21.h @@ -40,7 +40,18 @@ class Center2_Orb::Orb21 const int& mA2, const int& mB) const; + /// All magnetic components, ordered as [mA1][mA2][mB]. + std::vector cal_overlap(const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const; + std::vector> cal_grad_overlap( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB) const; + private: + template + std::vector cal_magnetic_components( + const ModuleBase::Vector3 &RA, const ModuleBase::Vector3 &RB, + std::vector (Orb11::*evaluate)(const ModuleBase::Vector3 &, const ModuleBase::Vector3 &) + const) const; + const Numerical_Orbital_Lm& nA1; const Numerical_Orbital_Lm& nA2; const Numerical_Orbital_Lm& nB; diff --git a/source/source_lcao/module_ri/matrix_orbs21.hpp b/source/source_lcao/module_ri/matrix_orbs21.hpp index b6f464d7440..1be8cc4deb3 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.hpp +++ b/source/source_lcao/module_ri/matrix_orbs21.hpp @@ -36,31 +36,35 @@ RI::Tensor Matrix_Orbs21::cal_overlap_matrix( default: throw std::invalid_argument(std::string(__FILE__)+" line "+std::to_string(__LINE__)); } + const ModuleBase::Vector3 position_A = tauA * lat0; + const ModuleBase::Vector3 position_B = tauB * lat0; for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; for( const auto &co4 : co3.second ) { const size_t NA1 = co4.first; - for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 ) + for (const auto &co5 : co4.second) { - for( const auto &co5 : co4.second ) + const int LA2 = co5.first; + for (const auto &co6 : co5.second) { - const int LA2 = co5.first; - for( const auto &co6 : co5.second ) + const size_t NA2 = co6.first; + for (const auto &co7 : co6.second) { - const size_t NA2 = co6.first; - for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 ) + const int LB = co7.first; + for (const auto &co8 : co7.second) { - for( const auto &co7 : co6.second ) + const size_t NB = co8.first; + const auto overlaps = co8.second.cal_overlap(position_A, position_B); + for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { - const int LB = co7.first; - for( const auto &co8 : co7.second ) + for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) { - const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const Tdata overlap = co8.second.cal_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB ); + const size_t offset = (MA1 * (2 * LA2 + 1) + MA2) * (2 * LB + 1) + MB; + const Tdata overlap = overlaps[offset]; const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; @@ -116,31 +120,36 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( } } + const ModuleBase::Vector3 position_A = tauA * lat0; + const ModuleBase::Vector3 position_B = tauB * lat0; for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; for( const auto &co4 : co3.second ) { const size_t NA1 = co4.first; - for( size_t MA1=0; MA1!=2*LA1+1; ++MA1 ) + for (const auto &co5 : co4.second) { - for( const auto &co5 : co4.second ) + const int LA2 = co5.first; + for (const auto &co6 : co5.second) { - const int LA2 = co5.first; - for( const auto &co6 : co5.second ) + const size_t NA2 = co6.first; + for (const auto &co7 : co6.second) { - const size_t NA2 = co6.first; - for( size_t MA2=0; MA2!=2*LA2+1; ++MA2 ) + const int LB = co7.first; + for (const auto &co8 : co7.second) { - for( const auto &co7 : co6.second ) + const size_t NB = co8.first; + const auto overlaps = co8.second.cal_grad_overlap(position_A, position_B); + for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { - const int LB = co7.first; - for( const auto &co8 : co7.second ) + for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) { - const size_t NB = co8.first; for( size_t MB=0; MB!=2*LB+1; ++MB ) { - const std::array grad_overlap = RI_Util::Vector3_to_array3(co8.second.cal_grad_overlap( tauA*lat0, tauB*lat0, MA1, MA2, MB )); + const size_t offset = (MA1 * (2 * LA2 + 1) + MA2) * (2 * LB + 1) + MB; + const std::array grad_overlap = + RI_Util::Vector3_to_array3(overlaps[offset]); const size_t iA1 = index_A1[TA][LA1][NA1][MA1]; const size_t iA2 = index_A2[TA][LA2][NA2][MA2]; const size_t iB = index_B[TB][LB][NB][MB]; diff --git a/source/source_lcao/test/CMakeLists.txt b/source/source_lcao/test/CMakeLists.txt index 53fde5e20ad..37ea3b65060 100644 --- a/source/source_lcao/test/CMakeLists.txt +++ b/source/source_lcao/test/CMakeLists.txt @@ -65,4 +65,19 @@ if(ENABLE_LCAO) device SOURCES test_init_chg_hr_error.cpp) + addtest( + TARGET + MODULE_LCAO_integral_reuse + LIBS + parameter + base + device + container + SOURCES + test_center2orb_orb21.cpp + ../center2orb.cpp + ../center2orb_orb11.cpp + ../center2orb_orb21.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_atomic_lm.cpp + ${ABACUS_SOURCE_DIR}/source_basis/module_ao/orb_gaunt_table.cpp) endif() diff --git a/source/source_lcao/test/test_center2orb_orb21.cpp b/source/source_lcao/test/test_center2orb_orb21.cpp new file mode 100644 index 00000000000..3bd581554fd --- /dev/null +++ b/source/source_lcao/test/test_center2orb_orb21.cpp @@ -0,0 +1,261 @@ +#include "source_lcao/center2orb_orb11.h" +#include "source_lcao/center2orb_orb21.h" + +#include "gtest/gtest.h" +#include +#include + +namespace +{ +Numerical_Orbital_Lm make_orbital( + const int l, + const int index, + const int nr, + const int nk) +{ + const double dr = 0.02; + const double cutoff = (nr - 1) * dr; + std::vector radial(nr); + std::vector rab(nr, dr); + std::vector values(nr); + for (int ir = 0; ir < nr; ++ir) + { + const double r = ir * dr; + radial[ir] = r; + values[ir] = std::pow(r, l) * std::pow(cutoff - r, 2) * std::exp(-(0.7 + 0.03 * index) * r); + } + Numerical_Orbital_Lm orbital; + const double *rab_ptr = rab.data(); + const double *radial_ptr = radial.data(); + const double *values_ptr = values.data(); + orbital.set_orbital_info( + "X", + 0, + l, + index, + nr, + rab_ptr, + radial_ptr, + Numerical_Orbital_Lm::Psi_Type::Psi, + values_ptr, + nk, + 0.1, + 0.01, + false, + true, + true); + return orbital; +} + +} // namespace + +TEST(Center2Orb11, AllMagneticComponentsMatchScalarEvaluation) +{ + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(6); + gaunt.init_Gaunt(6); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(10, 405, 65); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (int la : {0, + 1, + 3, + 6}) + { + for (int lb : {0, 2, 3}) + { + const auto left = make_orbital( + la, + 0, + 101, + 65); + const auto right = make_orbital( + lb, + 1, + 81, + 65); + Center2_Orb::Orb11 orbital( + left, + right, + &bessel, + gaunt); + orbital.init_radial_table(); + for (double distance : {0.0, + 1e-12, + 5e-11, + 2e-10, + 0.037, + 0.73, + 1.91, + 3.57, + 3.8}) + { + const double x = 0.8 * distance; + const double y = 0.36 * distance; + const double z = 0.48 * distance; + const ModuleBase::Vector3 position(x, y, z); + const auto values = orbital.cal_overlap(origin, position); + const auto gradients = orbital.cal_grad_overlap(origin, position); + const int count = (2 * la + 1) * (2 * lb + 1); + ASSERT_EQ(values.size(), count); + ASSERT_EQ(gradients.size(), count); + for (int ma = 0; ma <= 2 * la; ++ma) + { + for (int mb = 0; mb <= 2 * lb; ++mb) + { + const int offset = ma * (2 * lb + 1) + mb; + const double expected = orbital.cal_overlap( + origin, + position, + ma, + mb); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(values[offset], expected, tolerance); + const auto gradient = orbital.cal_grad_overlap( + origin, + position, + ma, + mb); + const double error = (gradients[offset] - gradient).norm(); + const double gradient_tolerance = 2e-10 * (1.0 + gradient.norm()); + EXPECT_LT(error, gradient_tolerance); + } + } + } + } + } +} + +TEST(Center2Orb21, AllMagneticComponentsMatchScalarEvaluation) +{ + // The l=6 auxiliary channel and f orbital produce channels through l=9. + const auto left = make_orbital( + 6, + 0, + 101, + 65); + const auto middle = make_orbital( + 3, + 1, + 101, + 65); + const auto right = make_orbital( + 3, + 2, + 81, + 65); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(9); + gaunt.init_Gaunt(9); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(13, 405, 65); + Center2_Orb::Orb21 orbital( + left, + middle, + right, + &bessel, + gaunt); + orbital.init_radial_table(); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + for (double distance : {0.0, + 1e-12, + 5e-11, + 2e-10, + 0.037, + 0.73, + 1.91, + 3.57, + 3.8}) + { + const double x = 0.8 * distance; + const double y = 0.36 * distance; + const double z = 0.48 * distance; + const ModuleBase::Vector3 position(x, y, z); + const auto values = orbital.cal_overlap(origin, position); + const auto gradients = orbital.cal_grad_overlap(origin, position); + ASSERT_EQ(values.size(), 13 * 7 * 7); + ASSERT_EQ(gradients.size(), values.size()); + for (int ma = 0; ma < 13; ++ma) + { + for (int mb = 0; mb < 7; ++mb) + { + for (int mc = 0; mc < 7; ++mc) + { + const int offset = (ma * 7 + mb) * 7 + mc; + const double expected = orbital.cal_overlap( + origin, + position, + ma, + mb, + mc); + const double tolerance = 2e-11 * (1.0 + std::abs(expected)); + EXPECT_NEAR(values[offset], expected, tolerance); + const auto gradient = orbital.cal_grad_overlap( + origin, + position, + ma, + mb, + mc); + const double error = (gradients[offset] - gradient).norm(); + const double gradient_tolerance = 2e-10 * (1.0 + gradient.norm()); + EXPECT_LT(error, gradient_tolerance); + } + } + } + } +} + +TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) +{ + const auto left = make_orbital( + 2, + 0, + 101, + 65); + const auto middle = make_orbital( + 1, + 1, + 101, + 65); + const auto right = make_orbital( + 3, + 2, + 101, + 65); + ORB_gaunt_table gaunt; + gaunt.init_Gaunt_CH(3); + gaunt.init_Gaunt(3); + ModuleBase::Sph_Bessel_Recursive::D2 bessel; + bessel.set_dx(0.001); + bessel.cal_jlx(7, 405, 65); + Center2_Orb::Orb21 orbital( + left, + middle, + right, + &bessel, + gaunt); + orbital.init_radial_table(); + const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); + const ModuleBase::Vector3 position(0.813, 0.321, 0.473); + const double step = 1e-5; + const auto gradients = orbital.cal_grad_overlap(origin, position); + for (int axis = 0; axis < 3; ++axis) + { + auto plus = position; + auto minus = position; + plus[axis] += step; + minus[axis] -= step; + const auto upper = orbital.cal_overlap(origin, plus); + const auto lower = orbital.cal_overlap(origin, minus); + ASSERT_EQ(upper.size(), gradients.size()); + ASSERT_EQ(lower.size(), gradients.size()); + for (std::size_t i = 0; i < gradients.size(); ++i) + { + const double finite_difference = (upper[i] - lower[i]) / (2 * step); + const double tolerance = 2e-7 * (1.0 + std::abs(finite_difference)); + EXPECT_NEAR(gradients[i][axis], finite_difference, tolerance); + } + } +} From 1550aaf820db082c1c9259f7e7d70d7d8616108e Mon Sep 17 00:00:00 2001 From: Jiacheng Xu <13862180016@163.com> Date: Mon, 2 Nov 2026 17:45:03 +0800 Subject: [PATCH 2/2] style(lcao): replace auto with explicit types --- source/source_lcao/center2orb_orb11.cpp | 4 +- source/source_lcao/center2orb_orb21.cpp | 6 +-- .../source_lcao/module_ri/matrix_orbs21.hpp | 36 +++++++++--------- .../test/test_center2orb_orb21.cpp | 38 +++++++++---------- 4 files changed, 42 insertions(+), 42 deletions(-) diff --git a/source/source_lcao/center2orb_orb11.cpp b/source/source_lcao/center2orb_orb11.cpp index bb5aa5635e1..47f6b335696 100644 --- a/source/source_lcao/center2orb_orb11.cpp +++ b/source/source_lcao/center2orb_orb11.cpp @@ -120,7 +120,7 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3& RA, const double* Gaunt_Coefficients_ptr = &(this->MGT.Gaunt_Coefficients(idx1, idx2, 0)); const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); - for (const auto& tb_r: this->Table_r) + for (const std::pair>& tb_r: this->Table_r) { const int LAB = tb_r.first; const double distance_power = std::pow(distance, LAB); @@ -198,7 +198,7 @@ ModuleBase::Vector3 Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2 ModuleBase::Vector3 grad_overlap(0.0, 0.0, 0.0); const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_); - for (const auto& tb_r: this->Table_r) + for (const std::pair>& tb_r: this->Table_r) { const int LAB = tb_r.first; const double distance_power = std::pow(distance, LAB); diff --git a/source/source_lcao/center2orb_orb21.cpp b/source/source_lcao/center2orb_orb21.cpp index 8d09a300aae..366e0041e2c 100644 --- a/source/source_lcao/center2orb_orb21.cpp +++ b/source/source_lcao/center2orb_orb21.cpp @@ -116,7 +116,7 @@ double Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3& RA, double overlap = 0.0; - for (const auto& orb11: this->orb11s) + for (const std::pair& orb11: this->orb11s) { const int LA = orb11.first; @@ -150,7 +150,7 @@ ModuleBase::Vector3 Center2_Orb::Orb21::cal_grad_overlap(const ModuleBas ModuleBase::Vector3 grad_overlap(0.0, 0.0, 0.0); - for (const auto& orb11: this->orb11s) + for (const std::pair& orb11: this->orb11s) { const int LA = orb11.first; @@ -183,7 +183,7 @@ std::vector Center2_Orb::Orb21::cal_magnetic_components( const int count_B = 2 * this->nB.getL() + 1; const int count = count_A1 * count_A2 * count_B; std::vector result(count); - for (const auto &orbital : this->orb11s) + for (const std::pair &orbital : this->orb11s) { const int LA = orbital.first; const int count_A = 2 * LA + 1; diff --git a/source/source_lcao/module_ri/matrix_orbs21.hpp b/source/source_lcao/module_ri/matrix_orbs21.hpp index 1be8cc4deb3..250d869996c 100644 --- a/source/source_lcao/module_ri/matrix_orbs21.hpp +++ b/source/source_lcao/module_ri/matrix_orbs21.hpp @@ -38,25 +38,25 @@ RI::Tensor Matrix_Orbs21::cal_overlap_matrix( const ModuleBase::Vector3 position_A = tauA * lat0; const ModuleBase::Vector3 position_B = tauB * lat0; - for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) + for( const std::pair>>>>> &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; - for( const auto &co4 : co3.second ) + for( const std::pair>>>> &co4 : co3.second ) { const size_t NA1 = co4.first; - for (const auto &co5 : co4.second) + for (const std::pair>>> &co5 : co4.second) { const int LA2 = co5.first; - for (const auto &co6 : co5.second) + for (const std::pair>> &co6 : co5.second) { const size_t NA2 = co6.first; - for (const auto &co7 : co6.second) + for (const std::pair> &co7 : co6.second) { const int LB = co7.first; - for (const auto &co8 : co7.second) + for (const std::pair &co8 : co7.second) { const size_t NB = co8.first; - const auto overlaps = co8.second.cal_overlap(position_A, position_B); + const std::vector overlaps = co8.second.cal_overlap(position_A, position_B); for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) @@ -122,25 +122,25 @@ std::array,3> Matrix_Orbs21::cal_grad_overlap_matrix( const ModuleBase::Vector3 position_A = tauA * lat0; const ModuleBase::Vector3 position_B = tauB * lat0; - for( const auto &co3 : center2_orb21_s.at(TA).at(TB) ) + for( const std::pair>>>>> &co3 : center2_orb21_s.at(TA).at(TB) ) { const int LA1 = co3.first; - for( const auto &co4 : co3.second ) + for( const std::pair>>>> &co4 : co3.second ) { const size_t NA1 = co4.first; - for (const auto &co5 : co4.second) + for (const std::pair>>> &co5 : co4.second) { const int LA2 = co5.first; - for (const auto &co6 : co5.second) + for (const std::pair>> &co6 : co5.second) { const size_t NA2 = co6.first; - for (const auto &co7 : co6.second) + for (const std::pair> &co7 : co6.second) { const int LB = co7.first; - for (const auto &co8 : co7.second) + for (const std::pair &co8 : co7.second) { const size_t NB = co8.first; - const auto overlaps = co8.second.cal_grad_overlap(position_A, position_B); + const std::vector> overlaps = co8.second.cal_grad_overlap(position_A, position_B); for (size_t MA1 = 0; MA1 != 2 * LA1 + 1; ++MA1) { for (size_t MA2 = 0; MA2 != 2 * LA2 + 1; ++MA2) @@ -189,14 +189,14 @@ std::map>>>>> matrixes; - for( const auto &co1 : center2_orb21_s ) + for( const std::pair>>>>>>> &co1 : center2_orb21_s ) { const size_t TA = co1.first; for( size_t IA=0; IA!=ucell.atoms[TA].na; ++IA ) { const ModuleBase::Vector3 &tauA( ucell.atoms[TA].tau[IA] ); - for( const auto &co2 : co1.second ) + for( const std::pair>>>>>> &co2 : co1.second ) { const size_t TB = co2.first; for( size_t IB=0; IB!=ucell.atoms[TB].na; ++IB ) @@ -214,10 +214,10 @@ std::map>>>>> m1 : matrixes) { const size_t T = m1.first; - for( auto m2 : m1.second ) + for( std::pair>>>> m2 : m1.second ) { const size_t I = m2.first; matrixes[T][I][T][I].resize(1); diff --git a/source/source_lcao/test/test_center2orb_orb21.cpp b/source/source_lcao/test/test_center2orb_orb21.cpp index 3bd581554fd..77d029251f2 100644 --- a/source/source_lcao/test/test_center2orb_orb21.cpp +++ b/source/source_lcao/test/test_center2orb_orb21.cpp @@ -65,12 +65,12 @@ TEST(Center2Orb11, AllMagneticComponentsMatchScalarEvaluation) { for (int lb : {0, 2, 3}) { - const auto left = make_orbital( + const Numerical_Orbital_Lm left = make_orbital( la, 0, 101, 65); - const auto right = make_orbital( + const Numerical_Orbital_Lm right = make_orbital( lb, 1, 81, @@ -95,8 +95,8 @@ TEST(Center2Orb11, AllMagneticComponentsMatchScalarEvaluation) const double y = 0.36 * distance; const double z = 0.48 * distance; const ModuleBase::Vector3 position(x, y, z); - const auto values = orbital.cal_overlap(origin, position); - const auto gradients = orbital.cal_grad_overlap(origin, position); + const std::vector values = orbital.cal_overlap(origin, position); + const std::vector> gradients = orbital.cal_grad_overlap(origin, position); const int count = (2 * la + 1) * (2 * lb + 1); ASSERT_EQ(values.size(), count); ASSERT_EQ(gradients.size(), count); @@ -112,7 +112,7 @@ TEST(Center2Orb11, AllMagneticComponentsMatchScalarEvaluation) mb); const double tolerance = 2e-11 * (1.0 + std::abs(expected)); EXPECT_NEAR(values[offset], expected, tolerance); - const auto gradient = orbital.cal_grad_overlap( + const ModuleBase::Vector3 gradient = orbital.cal_grad_overlap( origin, position, ma, @@ -130,17 +130,17 @@ TEST(Center2Orb11, AllMagneticComponentsMatchScalarEvaluation) TEST(Center2Orb21, AllMagneticComponentsMatchScalarEvaluation) { // The l=6 auxiliary channel and f orbital produce channels through l=9. - const auto left = make_orbital( + const Numerical_Orbital_Lm left = make_orbital( 6, 0, 101, 65); - const auto middle = make_orbital( + const Numerical_Orbital_Lm middle = make_orbital( 3, 1, 101, 65); - const auto right = make_orbital( + const Numerical_Orbital_Lm right = make_orbital( 3, 2, 81, @@ -173,8 +173,8 @@ TEST(Center2Orb21, AllMagneticComponentsMatchScalarEvaluation) const double y = 0.36 * distance; const double z = 0.48 * distance; const ModuleBase::Vector3 position(x, y, z); - const auto values = orbital.cal_overlap(origin, position); - const auto gradients = orbital.cal_grad_overlap(origin, position); + const std::vector values = orbital.cal_overlap(origin, position); + const std::vector> gradients = orbital.cal_grad_overlap(origin, position); ASSERT_EQ(values.size(), 13 * 7 * 7); ASSERT_EQ(gradients.size(), values.size()); for (int ma = 0; ma < 13; ++ma) @@ -192,7 +192,7 @@ TEST(Center2Orb21, AllMagneticComponentsMatchScalarEvaluation) mc); const double tolerance = 2e-11 * (1.0 + std::abs(expected)); EXPECT_NEAR(values[offset], expected, tolerance); - const auto gradient = orbital.cal_grad_overlap( + const ModuleBase::Vector3 gradient = orbital.cal_grad_overlap( origin, position, ma, @@ -209,17 +209,17 @@ TEST(Center2Orb21, AllMagneticComponentsMatchScalarEvaluation) TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) { - const auto left = make_orbital( + const Numerical_Orbital_Lm left = make_orbital( 2, 0, 101, 65); - const auto middle = make_orbital( + const Numerical_Orbital_Lm middle = make_orbital( 1, 1, 101, 65); - const auto right = make_orbital( + const Numerical_Orbital_Lm right = make_orbital( 3, 2, 101, @@ -240,15 +240,15 @@ TEST(Center2Orb21, ProductGradientMatchesFiniteDifference) const ModuleBase::Vector3 origin(0.0, 0.0, 0.0); const ModuleBase::Vector3 position(0.813, 0.321, 0.473); const double step = 1e-5; - const auto gradients = orbital.cal_grad_overlap(origin, position); + const std::vector> gradients = orbital.cal_grad_overlap(origin, position); for (int axis = 0; axis < 3; ++axis) { - auto plus = position; - auto minus = position; + ModuleBase::Vector3 plus = position; + ModuleBase::Vector3 minus = position; plus[axis] += step; minus[axis] -= step; - const auto upper = orbital.cal_overlap(origin, plus); - const auto lower = orbital.cal_overlap(origin, minus); + const std::vector upper = orbital.cal_overlap(origin, plus); + const std::vector lower = orbital.cal_overlap(origin, minus); ASSERT_EQ(upper.size(), gradients.size()); ASSERT_EQ(lower.size(), gradients.size()); for (std::size_t i = 0; i < gradients.size(); ++i)