Fix: switch to support the solid spherical harmonics that is consistent with other functionalities

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<double> i = {0., 1.};
+const double invsqrt2 = std::sqrt(2) * 0.5;
+
 std::complex<double> ModuleIO::cal_LzijR(
     const std::unique_ptr<TwoCenterIntegrator>& 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<double>& vR)
 {
+    if(mj == 0) {
+        return std::complex<double>(0.);
+    }
     double val_ = 0;
-    calculator->calculate(it, il, iz, mi, jt, jl, jz, mj, vR, &val_);
-    return std::complex<double>(mi) * val_;
+    calculator->calculate(it, il, iz, mi, jt, jl, jz, -mj, vR, &val_);
+    return i * static_cast<double>(mj) * val_;
 }
 
-std::complex<double> ModuleIO::cal_LyijR(
+std::complex<double> ModuleIO::cal_LxijR(
     const std::unique_ptr<TwoCenterIntegrator>& 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<double>& 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<double>(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<double> ModuleIO::cal_LxijR(
+std::complex<double> ModuleIO::cal_LyijR(
     const std::unique_ptr<TwoCenterIntegrator>& 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<double>& 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<double>(0.5) * (val_plus + val_minus);
+    assert(false); // inaccessible
 }
 
 ModuleIO::AngularMomentumCalculator::AngularMomentumCalculator(

source/source_io/module_hs/cal_pLpR.h

@@ -1,7 +1,11 @@
 /**
- * calculate the <phi_i|Lx/Ly/Lz|phi_j> 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 <phi_i|Lx/Ly/Lz|phi_j> 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,30 +17,51 @@
  *                                  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>
  * 
  * Likely, with L-, the spherical harmonic function Ylm can be lowered:
  * 
  *                          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:
  * 
  *                <l, m|Lx|l, m'> =   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:
  * 
  *                <l, m|Ly|l, m'> =   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:
  * 
  *                          <l, m|Lz|l, m'> = 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:
+ * 
+ *               <phi_i|Lx/Ly/Lz|phi_j> = <phi_i| (Lx/Ly/Lz |phi_j>)
+ * 
+ * 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
+ * 
+ *               <phi_i|Lz|phi_j> = mj <phi_i|phi_j>
+ * 
+ * For Lx = 1/2 * (L+ + L-):
+ * 
+ *               <phi_i|Lx|phi_j> = sqrt((lj-mj)(lj+mj+1)) * <phi_i|phi_j' > / 2
+ *                                + sqrt((lj+mj)(lj-mj+1)) * <phi_i|phi_j''> / 2
  * 
  * Technical Details
  * -----------------