Skip to content

perf(lcao): reuse common terms across magnetic components #8256

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mohanchen merged 3 commits into from
Nov 2, 2026
+573 −67
Merged

perf(lcao): reuse common terms across magnetic components #8256

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
06ab2be
perf(lcao): reuse common terms across magnetic components
Stardust0831 Nov 1, 2026
1550aaf
style(lcao): replace auto with explicit types
Stardust0831 Nov 2, 2026
e613e1a
Merge upstream develop into perf/ri-integral-reuse
Stardust0831 Nov 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
220 changes: 185 additions & 35 deletions source/source_lcao/center2orb_orb11.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,9 +119,19 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3<double>& RA,
const int idx2 = this->MGT.get_lm_index(LB, mB);
const double* Gaunt_Coefficients_ptr = &(this->MGT.Gaunt_Coefficients(idx1, idx2, 0));

for (const auto& tb_r: this->Table_r)
const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_);
for (const std::pair<const int, std::vector<double>>& 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;
Expand All @@ -138,17 +148,7 @@ double Center2_Orb::Orb11::cal_overlap(const ModuleBase::Vector3<double>& 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)
Expand Down Expand Up @@ -197,9 +197,26 @@ ModuleBase::Vector3<double> Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2

ModuleBase::Vector3<double> grad_overlap(0.0, 0.0, 0.0);

for (const auto& tb_r: this->Table_r)
const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, dr_);
for (const std::pair<const int, std::vector<double>>& 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;
{
Expand All @@ -211,30 +228,10 @@ ModuleBase::Vector3<double> 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<double> gylm_solid = grly[idx3];
const ModuleBase::Vector3<double> 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<double> gylm_real = (distance > tiny2) ? gylm_solid / distance_power : gylm_solid;

grad_overlap += i_exp // pow(2*PI,1.5)
* Gaunt_real_A_B_AB
Expand All @@ -244,3 +241,156 @@ ModuleBase::Vector3<double> Center2_Orb::Orb11::cal_grad_overlap( // caoyu add 2

return grad_overlap;
}

template <typename T> std::vector<T> Center2_Orb::Orb11::contract_angular(const std::vector<T> &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<T> 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<const int, std::vector<double>> &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<double> Center2_Orb::Orb11::cal_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const
{
const double tiny1 = 1e-12;
const double tiny2 = 1e-10;
const ModuleBase::Vector3<double> 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<double>(count);
}

const int lmax = LA + LB;
std::vector<double> 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<double> components(angular_count);
const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_);
for (const std::pair<const int, std::vector<double>> &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<ModuleBase::Vector3<double>> Center2_Orb::Orb11::cal_grad_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const
{
const double tiny1 = 1e-12;
const double tiny2 = 1e-10;
const ModuleBase::Vector3<double> 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<ModuleBase::Vector3<double>>(count);
}

const int lmax = LA + LB;
const int angular_count = (lmax + 1) * (lmax + 1);
const int gradient_count = 3 * angular_count;
std::vector<double> rly(angular_count);
std::vector<double> 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<ModuleBase::Vector3<double>> components(angular_count);
const int rmesh = Center2_Orb::get_rmesh(RcutA, RcutB, this->dr_);
for (const std::pair<const int, std::vector<double>> &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<double> solid_gradient(gx, gy, gz);
const ModuleBase::Vector3<double> 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);
}
7 changes: 7 additions & 0 deletions source/source_lcao/center2orb_orb11.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,14 @@ class Center2_Orb::Orb11
const int& mA,
const int& mB) const;

/// All magnetic components, ordered as [mA][mB].
std::vector<double> cal_overlap(const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const;
std::vector<ModuleBase::Vector3<double>> cal_grad_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const;

private:
template <typename T> std::vector<T> contract_angular(const std::vector<T> &components) const;

const Numerical_Orbital_Lm& nA;
const Numerical_Orbital_Lm& nB;

Expand Down
62 changes: 60 additions & 2 deletions source/source_lcao/center2orb_orb21.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ double Center2_Orb::Orb21::cal_overlap(const ModuleBase::Vector3<double>& RA,

double overlap = 0.0;

for (const auto& orb11: this->orb11s)
for (const std::pair<const int, Center2_Orb::Orb11>& orb11: this->orb11s)
{
const int LA = orb11.first;

Expand Down Expand Up @@ -150,7 +150,7 @@ ModuleBase::Vector3<double> Center2_Orb::Orb21::cal_grad_overlap(const ModuleBas

ModuleBase::Vector3<double> grad_overlap(0.0, 0.0, 0.0);

for (const auto& orb11: this->orb11s)
for (const std::pair<const int, Center2_Orb::Orb11>& orb11: this->orb11s)
{
const int LA = orb11.first;

Expand All @@ -169,3 +169,61 @@ ModuleBase::Vector3<double> Center2_Orb::Orb21::cal_grad_overlap(const ModuleBas

return grad_overlap;
}

template <typename T>
std::vector<T> Center2_Orb::Orb21::cal_magnetic_components(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB,
std::vector<T> (Orb11::*evaluate)(const ModuleBase::Vector3<double> &, const ModuleBase::Vector3<double> &)
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<T> result(count);
for (const std::pair<const int, Center2_Orb::Orb11> &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<T> 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<double> Center2_Orb::Orb21::cal_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const
{
return this->cal_magnetic_components<double>(RA, RB, &Orb11::cal_overlap);
}

std::vector<ModuleBase::Vector3<double>> Center2_Orb::Orb21::cal_grad_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const
{
return this->cal_magnetic_components<ModuleBase::Vector3<double>>(RA, RB, &Orb11::cal_grad_overlap);
}
11 changes: 11 additions & 0 deletions source/source_lcao/center2orb_orb21.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> cal_overlap(const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const;
std::vector<ModuleBase::Vector3<double>> cal_grad_overlap(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB) const;

private:
template <typename T>
std::vector<T> cal_magnetic_components(
const ModuleBase::Vector3<double> &RA, const ModuleBase::Vector3<double> &RB,
std::vector<T> (Orb11::*evaluate)(const ModuleBase::Vector3<double> &, const ModuleBase::Vector3<double> &)
const) const;

const Numerical_Orbital_Lm& nA1;
const Numerical_Orbital_Lm& nA2;
const Numerical_Orbital_Lm& nB;
Expand Down
Loading
Loading