Fix: correct Pauli-to-spinor Hamiltonian conversion for nspin=4

source/source_estate/module_dm/density_matrix.cpp

@@ -655,7 +655,7 @@ void DensityMatrix_Tools::func_xyz_to_updown<double>(const std::complex<double>
 {
     target_DMR_mat[icol + step_trace[0]] = tmp[0].real() + tmp[3].real();  // rho_0 = (rho_upup + rho_downdown).real()
     target_DMR_mat[icol + step_trace[1]] = tmp[1].real() + tmp[2].real();  // rho_x = (rho_updown + rho_downup).real()
-    target_DMR_mat[icol + step_trace[2]] = tmp[1].imag() - tmp[2].imag(); // rho_y = Im(rho_updown - rho_downup)
+    target_DMR_mat[icol + step_trace[2]] = -tmp[1].imag() + tmp[2].imag(); // rho_y = -Im(rho_updown) + Im(rho_downup)
     target_DMR_mat[icol + step_trace[3]] = tmp[0].real() - tmp[3].real();  // rho_z = (rho_upup - rho_downdown).real()
 }
 
@@ -664,7 +664,7 @@ void DensityMatrix_Tools::func_xyz_to_updown<std::complex<double>>(const std::co
 {
     target_DMR_mat[icol + step_trace[0]] = tmp[0] + tmp[3];                                         // rho_0 = (rho_upup + rho_downdown)
     target_DMR_mat[icol + step_trace[1]] = tmp[1] + tmp[2];                                         // rho_x = (rho_updown + rho_downup)
-    target_DMR_mat[icol + step_trace[2]] = -ModuleBase::IMAG_UNIT * (tmp[1] - tmp[2]); // rho_y = -i*(rho_updown - rho_downup)
+    target_DMR_mat[icol + step_trace[2]] = ModuleBase::IMAG_UNIT * (tmp[1] - tmp[2]); // rho_y = i*(rho_updown - rho_downup)
     target_DMR_mat[icol + step_trace[3]] = tmp[0] - tmp[3];                                         // rho_z = (rho_upup - rho_downdown)
 }
 

source/source_lcao/module_deltaspin/spin_constrain.cpp

@@ -142,7 +142,7 @@ int SpinConstrain<TK>::get_spin_sign(int ik) const
  * where P_at = sum_{l,m} |alpha_{l,m}><alpha_{l,m}| is the atomic projector.
  *
  * The magnetic moment is extracted via Pauli matrix traces:
- *   Mx = Re(occ[1] + occ[2]), My = Im(occ[1] - occ[2]), Mz = Re(occ[0] - occ[3])
+ *   Mx = Re(occ[1] + occ[2]), My = -Im(occ[1] - occ[2]), Mz = Re(occ[0] - occ[3])
  *
  * @par Algorithm (npol=1, collinear):
  * Only the z-component (spin projection) is computed:

source/source_lcao/module_deltaspin/spin_constrain.h

@@ -63,11 +63,11 @@ namespace spinconstrain
  * @brief Convert spinor occupation matrix to magnetic moment vector using Pauli matrices.
  *
  * @details For a two-component spinor wavefunction, the spin density matrix is:
- *   rho = |a|^2    a*b  |   = | (1+Mz)/2    (Mx+iMy)/2 |
- *         |b*a    |b|^2  |     | (Mx-iMy)/2   (1-Mz)/2  |
+ *   rho = |a|^2    a*b  |   = | (1+Mz)/2    (Mx-iMy)/2 |
+ *         |b*a    |b|^2  |     | (Mx+iMy)/2   (1-Mz)/2  |
  * The magnetic moment components are extracted via Pauli matrix traces:
  *   Mx = Tr(rho * sigma_x) = occ[1] + occ[2]           (real part)
- *   My = Tr(rho * sigma_y) = Im(occ[1] - occ[2])        (imaginary part)
+ *   My = Tr(rho * sigma_y) = -Im(occ[1] - occ[2])      (from sigma_y = [[0,-i],[i,0]])
  *   Mz = Tr(rho * sigma_z) = occ[0] - occ[3]            (real part)
  * where occ = {|a|^2, a*b, b*a, |b|^2} from becp coefficients.
  *
@@ -79,7 +79,7 @@ inline ModuleBase::Vector3<double> pauli_to_moment(const std::complex<double> oc
 {
     return ModuleBase::Vector3<double>(
         weight * (occ[1] + occ[2]).real(),
-        weight * (occ[1] - occ[2]).imag(),
+        -weight * (occ[1] - occ[2]).imag(),
         weight * (occ[0] - occ[3]).real()
     );
 }

source/source_lcao/module_deltaspin/test/deltaspin_core_test.cpp

@@ -23,15 +23,15 @@ struct Vec3i { int x, y, z; };
 // 1. pauli_to_moment: spinor -> magnetic moment
 //
 // Mx = w * (occ[1] + occ[2]).real()
-// My = w * (occ[1] - occ[2]).imag()
+// My = -w * (occ[1] - occ[2]).imag()  (from sigma_y = [[0,-i],[i,0]])
 // Mz = w * (occ[0] - occ[3]).real()
 // =====================================================================
 
 static Vec3 pauli_to_moment(const std::complex<double> occ[4], double weight)
 {
     return {
         weight * (occ[1] + occ[2]).real(),
-        weight * (occ[1] - occ[2]).imag(),
+        -weight * (occ[1] - occ[2]).imag(),
         weight * (occ[0] - occ[3]).real()
     };
 }
@@ -82,10 +82,10 @@ TEST_F(PauliToMomentTest, GeneralCase_AllComponents)
     occ[3] = {0.4, 0.0};
     auto M = pauli_to_moment(occ, 1.0);
     // Mx = (0.1+0.2i + 0.1-0.2i).real = 0.2
-    // My = (0.1+0.2i - (0.1-0.2i)).imag = (0+0.4i).imag = 0.4
+    // My = -(0.1+0.2i - (0.1-0.2i)).imag = -(0+0.4i).imag = -0.4
     // Mz = (0.6 - 0.4) = 0.2
     EXPECT_NEAR(M.x, 0.2, 1e-15);
-    EXPECT_NEAR(M.y, 0.4, 1e-15);
+    EXPECT_NEAR(M.y, -0.4, 1e-15);
     EXPECT_NEAR(M.z, 0.2, 1e-15);
 }
 

source/source_lcao/module_dftu/dftu_pw.cpp

@@ -309,8 +309,8 @@ void Plus_U::cal_occ_pw(const int iter,
                     }
                     vu_iat[index[0]] = 0.5 * (vu_tmp[0] + vu_tmp[3]);
                     vu_iat[index[3]] = 0.5 * (vu_tmp[0] - vu_tmp[3]);
-                    vu_iat[index[1]] = 0.5 * (vu_tmp[1] + std::complex<double>(0.0, 1.0) * vu_tmp[2]);
-                    vu_iat[index[2]] = 0.5 * (vu_tmp[1] - std::complex<double>(0.0, 1.0) * vu_tmp[2]);
+                    vu_iat[index[1]] = 0.5 * (vu_tmp[1] - std::complex<double>(0.0, 1.0) * vu_tmp[2]);
+                    vu_iat[index[2]] = 0.5 * (vu_tmp[1] + std::complex<double>(0.0, 1.0) * vu_tmp[2]);
                 }
             }
         }

source/source_lcao/module_operator_lcao/dspin_lcao.cpp

@@ -56,10 +56,10 @@ inline void cal_coeff_lambda(const std::vector<double>& current_lambda, std::vec
     coefficients[1] = -current_lambda[0];
 }
 inline void cal_coeff_lambda(const std::vector<double>& current_lambda, std::vector<std::complex<double>>& coefficients)
-{// {\lambda^{I,3}, \lambda^{I,1}+i\lambda^{I,2}, \lambda^{I,1}-i\lambda^{I,2}, -\lambda^{I,3}}
+{// {\lambda^{I,3}, \lambda^{I,1}-i\lambda^{I,2}, \lambda^{I,1}+i\lambda^{I,2}, -\lambda^{I,3}}
     coefficients[0] = std::complex<double>(current_lambda[2], 0.0);
-    coefficients[1] = std::complex<double>(current_lambda[0] , current_lambda[1]);
-    coefficients[2] = std::complex<double>(current_lambda[0] , -1 * current_lambda[1]);
+    coefficients[1] = std::complex<double>(current_lambda[0] , -current_lambda[1]);
+    coefficients[2] = std::complex<double>(current_lambda[0] , current_lambda[1]);
     coefficients[3] = std::complex<double>(-1 * current_lambda[2], 0.0);
 }
 

tests/01_PW/099_PW_DJ_SO/result.ref

@@ -1,5 +1,5 @@
-etotref -5662.3908859903258417
-etotperatomref -2831.1954429952
-totalforceref 17.965510
-totalstressref 100582.607209
+etotref -5662.3881388456420609
+etotperatomref -2831.1940694228
+totalforceref 15.774740
+totalstressref 100840.559090
 totaltimeref 1.26

tests/03_NAO_multik/scf_angle_spin4/result.ref

@@ -1,5 +1,5 @@
-etotref -6267.4651896196382950
-etotperatomref -3133.7325948098
-totalforceref 0.000000
-totalstressref 3912.920437
+etotref -6267.4651944939805617
+etotperatomref -3133.7325972470
+totalforceref 0.000008
+totalstressref 3912.920542
 totaltimeref 15.08

tests/03_NAO_multik/scf_out_dos_spin4/result.ref

@@ -1,6 +1,6 @@
-etotref -1964.0663947982770878
+etotref -1964.0663947982982336
 etotperatomref -982.0331973991
-totalforceref 0.162158
-totalstressref 1877.059089
+totalforceref 0.162298
+totalstressref 1877.059021
 totaldosref 38
 totaltimeref 16.23

tests/03_NAO_multik/scf_u_spin4/result.ref

@@ -1,5 +1,5 @@
-etotref -6789.2816406266510967
-etotperatomref -3394.6408203133
-totalforceref 11.331534
-totalstressref 4697.832232
+etotref -6789.1423886377124290
+etotperatomref -3394.5711943189
+totalforceref 14.359774
+totalstressref 4333.997
 totaltimeref 9.71