From 1e608e25a611f9fc518275adddf8cf16040325f9 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Sun, 19 Apr 2026 20:51:15 +0800 Subject: [PATCH 1/6] Fix: switch to support the solid spherical harmonics that is consistent with other functionalities --- source/source_io/module_hs/cal_pLpR.cpp | 148 ++++++++++++++++++------ source/source_io/module_hs/cal_pLpR.h | 43 +++++-- 2 files changed, 146 insertions(+), 45 deletions(-) diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index e437335ec0..b9798ee369 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -18,75 +18,151 @@ * * FIXME: the following part will be transfered to TwoCenterIntegrator soon * + * Notation + * -------- + * ylm: complex spherical harmonics + * slm: solid (real) spherical harmonics + * + * Changelog + * --------- + * Switch to support the solid spherical harmonics to keep consistent with + * implementations of other parts. + * Formulation (in Chinese): + * https://my.feishu.cn/wiki/D0enwcUKfiJgtSkJ5scc9Dagntc */ -// L+|l, m> = sqrt((l-m)(l+m+1))|l, m+1>, return the sqrt((l-m)(l+m+1)) -double _lplus_on_ylm(const int l, const int m) +// L+ylm = sqrt((l-m)(l+m+1))ylm+1, return the sqrt((l-m)(l+m+1)) +double _lambda_plus(const int l, const int m) { - return std::sqrt((l - m) * (l + m + 1)); + return std::sqrt((l - m) * (l + m + 1)); // NOTE: complex spherical harmonics } -// L-|l, m> = sqrt((l+m)(l-m+1))|l, m-1>, return the sqrt((l+m)(l-m+1)) -double _lminus_on_ylm(const int l, const int m) +// L-ylm = sqrt((l+m)(l-m+1))ylm-1, return the sqrt((l+m)(l-m+1)) +double _lambda_minus(const int l, const int m) { - return std::sqrt((l + m) * (l - m + 1)); + return std::sqrt((l + m) * (l - m + 1)); // NOTE: complex spherical harmonics } +const std::complex i = {0., 1.}; +const double invsqrt2 = std::sqrt(2) * 0.5; + std::complex ModuleIO::cal_LzijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int mi, const int jt, const int ja, const int jl, const int jz, const int mj, const ModuleBase::Vector3& vR) { + if(mj == 0) { + return std::complex(0.); + } double val_ = 0; - calculator->calculate(it, il, iz, mi, jt, jl, jz, mj, vR, &val_); - return std::complex(mi) * val_; + calculator->calculate(it, il, iz, mi, jt, jl, jz, -mj, vR, &val_); + return std::sqrt(2) * mj * i * val_; } -std::complex ModuleIO::cal_LyijR( +std::complex ModuleIO::cal_LxijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int im, const int jt, const int ja, const int jl, const int jz, const int jm, const ModuleBase::Vector3& vR) { - // Ly = -i/2 * (L+ - L-) - const double plus_ = _lplus_on_ylm(jl, jm); - const double minus_ = _lminus_on_ylm(jl, jm); - double val_plus = 0, val_minus = 0; - if (plus_ != 0) - { - calculator->calculate(it, il, iz, im, jt, jl, jz, jm + 1, vR, &val_plus); - val_plus *= plus_; + const double lmbdp = _lambda_plus(jl, jm); + const double lmbdm = _lambda_minus(jl, jm); + // two-center-integral placeholders + double valp = 0.; + double valm = 0.; + if (jm > 1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm+1), vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm-1), vR, &valm); + } + return i * 0.5 * (lmbdp * valp + lmbdm * valm); } - if (minus_ != 0) - { - calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); - val_minus *= minus_; + if (jm == 1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -2, vR, &valp); + } + return i * 0.5 * lmbdp * valp; + } + if (jm == 0) { + const double lmbd = _lambda_plus(jl, 0); // std::sqrt(jl*(jl+1)) + if (std::fabs(lmbd) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -1, vR, &valp); + } + return i * invsqrt2 * lmbd * valp; + } + if (jm == -1) { + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valm); + } + return -i * 0.5 * lmbdm * valm; + } + if (jm < -1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm+1), vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm-1), vR, &valm); + } + return -i * 0.5 * (lmbdp * valp + lmbdm * valm); } - return std::complex(0, -0.5) * (val_plus - val_minus); + assert(false); // inaccessible } -std::complex ModuleIO::cal_LxijR( +std::complex ModuleIO::cal_LyijR( const std::unique_ptr& calculator, const int it, const int ia, const int il, const int iz, const int im, const int jt, const int ja, const int jl, const int jz, const int jm, const ModuleBase::Vector3& vR) { - // Lx = 1/2 * (L+ + L-) - const double plus_ = _lplus_on_ylm(jl, jm); - const double minus_ = _lminus_on_ylm(jl, jm); - double val_plus = 0, val_minus = 0; - if (plus_ != 0) - { - calculator->calculate(it, il, iz, im, jt, jl, jz, jm + 1, vR, &val_plus); - val_plus *= plus_; + const double lmbdp = _lambda_plus(jl, jm); + const double lmbdm = _lambda_minus(jl, jm); + // two-center-integral placeholders + double valp = 0.; + double valm = 0.; + if (jm > 1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm+1, vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm-1, vR, &valm); + } + return -i * 0.5 * (lmbdp * valp - lmbdm * valm); } - if (minus_ != 0) - { - calculator->calculate(it, il, iz, im, jt, jl, jz, jm - 1, vR, &val_minus); - val_minus *= minus_; + if (jm == 1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 0, vR, &valm); + } + return -i * 0.5 * (lmbdp * valp - lmbdm * valm); + } + if (jm == 0) { + const double lmbd = _lambda_plus(jl, 0); // std::sqrt(l*(l+1)) + if (std::fabs(lmbd) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 1, vR, &valp); + } + return -i * invsqrt2 * lmbd * valp; + } + if (jm == -1) { + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, -2, vR, &valm); + } + return -i * 0.5 * lmbdm * valm; + } + if (jm < -1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm+1, vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, jm-1, vR, &valm); + } + return i * 0.5 * (lmbdp * valp - lmbdm * valm); } - return std::complex(0.5) * (val_plus + val_minus); + assert(false); // inaccessible } ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator( diff --git a/source/source_io/module_hs/cal_pLpR.h b/source/source_io/module_hs/cal_pLpR.h index 978a89f466..9ecd15a612 100644 --- a/source/source_io/module_hs/cal_pLpR.h +++ b/source/source_io/module_hs/cal_pLpR.h @@ -1,7 +1,11 @@ /** - * calculate the matrix elements, in which the Lx/Ly/Lz - * are the angular momentum operators, |phi_i> and |phi_j> are the numerical - * atomic orbitals (NAOs). + * calculate the matrix elements with the ACA (atom-centered + * approximation), in which the Lx/Ly/Lz are the angular momentum operators, + * |phi_i> and |phi_j> are the numerical atomic orbitals (NAOs). + * + * Note: this module is the most suitable for the case of SOC (spin-orbit coupling) + * calculation. A strict implmentation requires the evaluation on the position operator + * and the momentum operator, is beyond the scope of present implementation. * * Formulation * ----------- @@ -13,8 +17,8 @@ * Lx = (L+ + L-) / 2 * Ly = (L+ - L-) / 2i * - * With L+, the spherical harmonic function Ylm (denoted as |l, m> in the following) - * can be raised: + * With L+, the COMPLEX spherical harmonic function Ylm (denoted as |l, m> in + * the following) can be raised: * * L+|l, m> = sqrt((l-m)(l+m+1))|l, m+1> * @@ -22,21 +26,42 @@ * * L-|l, m> = sqrt((l+m)(l-m+1))|l, m-1> * - * Therefore the Lx matrix element can be calculated as: + * Therefore, for the case where there is only one atom, the Lx matrix element + * can be calculated as: * * = sqrt((l-m)(l+m+1)) * delta(m, m'+1) / 2 * + sqrt((l+m)(l-m+1)) * delta(m, m'-1) / 2 * - * The Ly matrix element can be calculated as: + * Likewise the Ly matrix element can be calculated as: * * = sqrt((l-m)(l+m+1)) * delta(m, m'+1) / 2i * - sqrt((l+m)(l-m+1)) * delta(m, m'-1) / 2i * - * The Lz matrix element can be calculated as: + * and the Lz matrix element can be calculated as: * * = m * delta(m, m') * - * However, things will change when there are more than one centers. + * ABACUS employs the REAL spherical harmonics, the transformation between the + * REAL and COMPLEX spherical harmonics can be found from: + * https://abacus.deepmodeling.com/en/latest/advanced/pp_orb.html + * + * For the case where there are multiple atoms, present calculation is based on + * the Atom-Centered-Approximation (ACA), in which the global angular momentum + * operator is decomposed to the sum of the local opeartor that only considers + * one atom. Therefore, the matrix element can be calculated by combining the + * ladder operators with the two-center-integral technique: + * + * = ) + * + * the two-center-integral will be performed after the application of the L + * operator. For example after the application of Lz, the equation above yields + * + * = mj + * + * For Lx = 1/2 * (L+ + L-): + * + * = sqrt((lj-mj)(lj+mj+1)) * / 2 + * + sqrt((lj+mj)(lj-mj+1)) * / 2 * * Technical Details * ----------------- From 736f78afa90b8c95e744bfcf45c1a2954de9d5f5 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Sun, 19 Apr 2026 23:15:54 +0800 Subject: [PATCH 2/6] correct impl. and update tests --- source/source_io/module_hs/cal_pLpR.cpp | 19 +-------- source/source_io/test/cal_pLpR_test.cpp | 54 ++++++++++--------------- 2 files changed, 24 insertions(+), 49 deletions(-) diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index b9798ee369..89ab05688c 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -93,13 +93,7 @@ std::complex ModuleIO::cal_LxijR( } return i * invsqrt2 * lmbd * valp; } - if (jm == -1) { - if (std::fabs(lmbdm) > 1e-12) { - calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valm); - } - return -i * 0.5 * lmbdm * valm; - } - if (jm < -1) { + if (jm <= -1) { if (std::fabs(lmbdp) > 1e-12) { calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm+1), vR, &valp); } @@ -122,7 +116,7 @@ std::complex ModuleIO::cal_LyijR( // two-center-integral placeholders double valp = 0.; double valm = 0.; - if (jm > 1) { + if (jm >= 1) { if (std::fabs(lmbdp) > 1e-12) { calculator->calculate(it, il, iz, im, jt, jl, jz, jm+1, vR, &valp); } @@ -131,15 +125,6 @@ std::complex ModuleIO::cal_LyijR( } return -i * 0.5 * (lmbdp * valp - lmbdm * valm); } - if (jm == 1) { - if (std::fabs(lmbdp) > 1e-12) { - calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valp); - } - if (std::fabs(lmbdm) > 1e-12) { - calculator->calculate(it, il, iz, im, jt, jl, jz, 0, vR, &valm); - } - return -i * 0.5 * (lmbdp * valp - lmbdm * valm); - } if (jm == 0) { const double lmbd = _lambda_plus(jl, 0); // std::sqrt(l*(l+1)) if (std::fabs(lmbd) > 1e-12) { diff --git a/source/source_io/test/cal_pLpR_test.cpp b/source/source_io/test/cal_pLpR_test.cpp index 2dad5f663d..42d086273f 100644 --- a/source/source_io/test/cal_pLpR_test.cpp +++ b/source/source_io/test/cal_pLpR_test.cpp @@ -54,7 +54,7 @@ TEST_F(CalpLpRTest, CalLzijRTest) int it = 0, ia = 0, il = 0, iz = 0, mi = 0; int jt = 0, ja = 0, jl = 0, jz = 0, mj = 0; // self, the first s - // = 0: no magnetic moment + // = 0: Lz|s> = 0 out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); @@ -63,26 +63,26 @@ TEST_F(CalpLpRTest, CalLzijRTest) out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: orthogonal + // = 0: Lz|p(m=0)> = 0 il = 1; mi = -1; jl = 1; mj = 0; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 1: same + // = 0: sqrt(2)i = 0 il = 1; mi = 1; jl = 1; mj = 1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 1.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: orthogonal - il = 2; mi = -1; jl = 2; mj = 0; + // = -sqrt(2)i = -sqrt(2)i + il = 1; mi = 1; jl = 1; mj = -1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 1: same - il = 2; mi = 1; jl = 2; mj = 1; + EXPECT_NEAR(out.imag(), -sqrt(2.0), DOUBLETHRESHOLD); + // = sqrt(2)i = sqrt(2)i + il = 1; mi = -1; jl = 1; mj = 1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 1.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), sqrt(2.0), DOUBLETHRESHOLD); } TEST_F(CalpLpRTest, CalLxijRTest) @@ -102,21 +102,16 @@ TEST_F(CalpLpRTest, CalLxijRTest) out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0.5: Lx|p(m=0)> = 1/2 (|p(m=-1)> + |p(m=1)>) + // = i = i il = 1; im = -1; jl = 1; jm = 0; out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.5*sqrt(2.0), DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: expectation value is 0 - il = 1; im = 1; jl = 1; jm = 1; - out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0.5: Lx|d(m=0)> = 1/2 (|d(m=-1)> + |d(m=1)>) - il = 2; im = -1; jl = 2; jm = 0; + EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); + // = -i = -i + il = 1; im = 0; jl = 1; jm = -1; out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.5*sqrt(6.0), DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); } TEST_F(CalpLpRTest, CalLyijRTest) @@ -136,21 +131,16 @@ TEST_F(CalpLpRTest, CalLyijRTest) out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = -i/2: Ly|p(m=0)> = -i/2 (|p(m=1)> - |p(m=-1)>) + // = -i* = -i il = 1; im = -1; jl = 1; jm = 0; out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.5*sqrt(2.0), DOUBLETHRESHOLD); - // = 0: expectation value is 0 - il = 1; im = 1; jl = 1; jm = 1; - out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = -i/2: Ly|d(m=0)> = -i/2 (|d(m=1)> - |d(m=-1)>) - il = 2; im = -1; jl = 2; jm = 0; + EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); + // = i* = i + il = 1; im = 0; jl = 1; jm = 1; out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.5*sqrt(6.0), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); } int main(int argc, char **argv) From b04d621ae30945b627f5576bbfd2e1907b63c97b Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Mon, 20 Apr 2026 11:13:26 +0800 Subject: [PATCH 3/6] correct impl.: add the sqrt(2) coefficient for sl0 --- source/source_io/module_hs/cal_pLpR.cpp | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index 89ab05688c..ba9e6b18a8 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -93,7 +93,16 @@ std::complex ModuleIO::cal_LxijR( } return i * invsqrt2 * lmbd * valp; } - if (jm <= -1) { + if (jm == -1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 0, vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valm); + } + return -i * 0.5 * (std::sqrt(2) * lmbdp * valp + lmbdm * valm); + } + if (jm < -1) { if (std::fabs(lmbdp) > 1e-12) { calculator->calculate(it, il, iz, im, jt, jl, jz, -(jm+1), vR, &valp); } @@ -116,7 +125,7 @@ std::complex ModuleIO::cal_LyijR( // two-center-integral placeholders double valp = 0.; double valm = 0.; - if (jm >= 1) { + if (jm > 1) { if (std::fabs(lmbdp) > 1e-12) { calculator->calculate(it, il, iz, im, jt, jl, jz, jm+1, vR, &valp); } @@ -125,6 +134,15 @@ std::complex ModuleIO::cal_LyijR( } return -i * 0.5 * (lmbdp * valp - lmbdm * valm); } + if (jm == 1) { + if (std::fabs(lmbdp) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 2, vR, &valp); + } + if (std::fabs(lmbdm) > 1e-12) { + calculator->calculate(it, il, iz, im, jt, jl, jz, 0, vR, &valm); + } + return -i * 0.5 * (lmbdp * valp - std::sqrt(2) * lmbdm * valm); + } if (jm == 0) { const double lmbd = _lambda_plus(jl, 0); // std::sqrt(l*(l+1)) if (std::fabs(lmbd) > 1e-12) { From 3319eb054de9715a3abc45e3ed38e30bc610ea72 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Mon, 20 Apr 2026 11:20:34 +0800 Subject: [PATCH 4/6] correct the impl. of Lz --- source/source_io/module_hs/cal_pLpR.cpp | 2 +- source/source_io/test/cal_pLpR_test.cpp | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index ba9e6b18a8..41f1488a28 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -57,7 +57,7 @@ std::complex ModuleIO::cal_LzijR( } double val_ = 0; calculator->calculate(it, il, iz, mi, jt, jl, jz, -mj, vR, &val_); - return std::sqrt(2) * mj * i * val_; + return i * static_cast(mj) * val_; } std::complex ModuleIO::cal_LxijR( diff --git a/source/source_io/test/cal_pLpR_test.cpp b/source/source_io/test/cal_pLpR_test.cpp index 42d086273f..fedb48f6d9 100644 --- a/source/source_io/test/cal_pLpR_test.cpp +++ b/source/source_io/test/cal_pLpR_test.cpp @@ -68,21 +68,21 @@ TEST_F(CalpLpRTest, CalLzijRTest) out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: sqrt(2)i = 0 + // = 0: i = 0 il = 1; mi = 1; jl = 1; mj = 1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = -sqrt(2)i = -sqrt(2)i + // = -i = -i il = 1; mi = 1; jl = 1; mj = -1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), -sqrt(2.0), DOUBLETHRESHOLD); - // = sqrt(2)i = sqrt(2)i + EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); + // = i = i il = 1; mi = -1; jl = 1; mj = 1; out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), sqrt(2.0), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); } TEST_F(CalpLpRTest, CalLxijRTest) From d04f314771b2a8acca13a3637a66c55c3b52b2db Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Mon, 20 Apr 2026 11:47:51 +0800 Subject: [PATCH 5/6] update unittest --- source/source_io/test/cal_pLpR_test.cpp | 135 +++++++++++++----------- 1 file changed, 73 insertions(+), 62 deletions(-) diff --git a/source/source_io/test/cal_pLpR_test.cpp b/source/source_io/test/cal_pLpR_test.cpp index fedb48f6d9..19f53f5859 100644 --- a/source/source_io/test/cal_pLpR_test.cpp +++ b/source/source_io/test/cal_pLpR_test.cpp @@ -1,3 +1,21 @@ +/** + * Unit-test of cal_pLpR.cpp + * + * The representation matrices for Lx, Ly and Lz operators (SO(3) generators). + * On basis {s}: + * Lx: Ly: Lz: + * 0 0 0 + * + * On basis {px, py, pz}: + * Lx: Ly: Lz: + * px py pz px py pz px py pz + * px 0 px 0 i px 0 -i + * py 0 -i py 0 py i 0 + * pz i 0 pz -i 0 pz 0 + * + * will test if the calculated values are close to the above results. + */ + #include #include #include @@ -51,38 +69,29 @@ TEST_F(CalpLpRTest, CalLzijRTest) std::complex out; ModuleBase::Vector3 vR(0., 0., 0.); // home-cell - int it = 0, ia = 0, il = 0, iz = 0, mi = 0; - int jt = 0, ja = 0, jl = 0, jz = 0, mj = 0; // self, the first s + int it = 0, ia = 0, il = 0, iz = 0, im = 0; + int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; - // = 0: Lz|s> = 0 - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: orthogonal - jl = 1; - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: Lz|p(m=0)> = 0 - il = 1; mi = -1; jl = 1; mj = 0; - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: i = 0 - il = 1; mi = 1; jl = 1; mj = 1; - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); + // l=0 + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = -i = -i - il = 1; mi = 1; jl = 1; mj = -1; - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); - // = i = i - il = 1; mi = -1; jl = 1; mj = 1; - out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, mi, jt, ja, jl, jz, mj, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); + + // l=1 + il = 1; jl = 1; + std::vector> ans(9, 0.0); + ans[1] = {0.0, -1.0}; ans[3] = {0.0, 1.0}; + int idx = 0; + const std::vector m = {1, -1, 0}; // px, py, pz + for (auto im_: m) { + for (auto jm_: m) { + im = im_; jm = jm_; + out = ModuleIO::cal_LzijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), ans[idx].real(), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), ans[idx].imag(), DOUBLETHRESHOLD); + idx++; + } + } } TEST_F(CalpLpRTest, CalLxijRTest) @@ -91,27 +100,28 @@ TEST_F(CalpLpRTest, CalLxijRTest) ModuleBase::Vector3 vR(0., 0., 0.); // home-cell int it = 0, ia = 0, il = 0, iz = 0, im = 0; - int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; // self, the first s + int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; - // = 0: no anisotropy - out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: orthogonal - jl = 1; + // l=0 out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = i = i - il = 1; im = -1; jl = 1; jm = 0; - out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); - // = -i = -i - il = 1; im = 0; jl = 1; jm = -1; - out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); + + // l=1 + il = 1; jl = 1; + std::vector> ans(9, 0.0); + ans[5] = {0.0, -1.0}; ans[7] = {0.0, 1.0}; + int idx = 0; + const std::vector m = {1, -1, 0}; // px, py, pz + for (auto im_: m) { + for (auto jm_: m) { + im = im_; jm = jm_; + out = ModuleIO::cal_LxijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), ans[idx].real(), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), ans[idx].imag(), DOUBLETHRESHOLD); + idx++; + } + } } TEST_F(CalpLpRTest, CalLyijRTest) @@ -122,25 +132,26 @@ TEST_F(CalpLpRTest, CalLyijRTest) int it = 0, ia = 0, il = 0, iz = 0, im = 0; int jt = 0, ja = 0, jl = 0, jz = 0, jm = 0; // self, the first s - // = 0: no anisotropy + // l=0 out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = 0: orthogonal - jl = 1; - out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 0.0, DOUBLETHRESHOLD); - // = -i* = -i - il = 1; im = -1; jl = 1; jm = 0; - out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), -1.0, DOUBLETHRESHOLD); - // = i* = i - il = 1; im = 0; jl = 1; jm = 1; - out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); - EXPECT_NEAR(out.real(), 0.0, DOUBLETHRESHOLD); - EXPECT_NEAR(out.imag(), 1.0, DOUBLETHRESHOLD); + + // l=1 + il = 1; jl = 1; + std::vector> ans(9, 0.0); + ans[2] = {0.0, 1.0}; ans[6] = {0.0, -1.0}; + int idx = 0; + const std::vector m = {1, -1, 0}; // px, py, pz + for (auto im_: m) { + for (auto jm_: m) { + im = im_; jm = jm_; + out = ModuleIO::cal_LyijR(calculator_, it, ia, il, iz, im, jt, ja, jl, jz, jm, vR); + EXPECT_NEAR(out.real(), ans[idx].real(), DOUBLETHRESHOLD); + EXPECT_NEAR(out.imag(), ans[idx].imag(), DOUBLETHRESHOLD); + idx++; + } + } } int main(int argc, char **argv) From a6d6f14a109824e0a22bbaf4d0ec1f8554d5d6d7 Mon Sep 17 00:00:00 2001 From: kirk0830 Date: Mon, 20 Apr 2026 12:14:18 +0800 Subject: [PATCH 6/6] correct the unittest --- source/source_io/test/cal_pLpR_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/source_io/test/cal_pLpR_test.cpp b/source/source_io/test/cal_pLpR_test.cpp index 19f53f5859..99b394568e 100644 --- a/source/source_io/test/cal_pLpR_test.cpp +++ b/source/source_io/test/cal_pLpR_test.cpp @@ -110,7 +110,7 @@ TEST_F(CalpLpRTest, CalLxijRTest) // l=1 il = 1; jl = 1; std::vector> ans(9, 0.0); - ans[5] = {0.0, -1.0}; ans[7] = {0.0, 1.0}; + ans[5] = {0.0, 1.0}; ans[7] = {0.0, -1.0}; int idx = 0; const std::vector m = {1, -1, 0}; // px, py, pz for (auto im_: m) { @@ -140,7 +140,7 @@ TEST_F(CalpLpRTest, CalLyijRTest) // l=1 il = 1; jl = 1; std::vector> ans(9, 0.0); - ans[2] = {0.0, 1.0}; ans[6] = {0.0, -1.0}; + ans[2] = {0.0, -1.0}; ans[6] = {0.0, 1.0}; int idx = 0; const std::vector m = {1, -1, 0}; // px, py, pz for (auto im_: m) {