diff --git a/source/source_io/module_hs/cal_pLpR.cpp b/source/source_io/module_hs/cal_pLpR.cpp index e437335ec0..41f1488a28 100644 --- a/source/source_io/module_hs/cal_pLpR.cpp +++ b/source/source_io/module_hs/cal_pLpR.cpp @@ -18,75 +18,154 @@ * * 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 i * static_cast(mj) * 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(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); } - return std::complex(0, -0.5) * (val_plus - val_minus); + 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); + } + 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 - 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) { + 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 * ----------------- diff --git a/source/source_io/test/cal_pLpR_test.cpp b/source/source_io/test/cal_pLpR_test.cpp index 2dad5f663d..99b394568e 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: no magnetic moment - 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: orthogonal - il = 1; mi = -1; jl = 1; mj = 0; - 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); - // = 1: same - 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); - // = 0: orthogonal - il = 2; mi = -1; jl = 2; 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 - il = 2; mi = 1; jl = 2; 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); + + // 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,32 +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 + // 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); - // = 0: orthogonal - jl = 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|p(m=0)> = 1/2 (|p(m=-1)> + |p(m=1)>) - 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; - 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); + + // 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) @@ -127,30 +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 - 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; + // 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); - // = -i/2: Ly|p(m=0)> = -i/2 (|p(m=1)> - |p(m=-1)>) - 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; - 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); + + // 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)