Fix: correct Pauli-to-spinor Hamiltonian conversion for nspin=4 (#7664)

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

Fix two bugs in LCAO non-collinear Hamiltonian construction:

1. Wrong sign in off-diagonal elements: H_{up,down} = B_x + i*B_y (wrong)
   should be B_x - i*B_y (correct), and vice versa for H_{down,up}.
   Fixed by correcting clx_j coefficients in merge_hr_part_to_hR().

2. Missing complex conjugate in lower triangle fill: H(-R) used transpose
   instead of conjugate transpose, breaking Hermiticity for complex matrices.
   Fixed by using std::conj() when filling lower triangle.

These errors caused the non-collinear Hamiltonian to be the complex conjugate
of the correct result, leading to incorrect spin textures in nspin=4 calculations.
The PW code path was not affected.

Add test case and verification script to validate:
- H(R=0) Hermiticity: max|H - H^dagger| < 1e-10
- Off-diagonal phase: Im(H_{up,down}) < 0 for m||+y direction

See tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md for details.

* fix: correct Pauli-to-spinor conversion in DFT+U and DeltaSpin for nspin=4

Fix three critical bugs in non-collinear (nspin=4) LCAO calculations:

1. DFT+U transfer_vu (dftu_lcao.cpp): Fix sign error in Pauli-to-spinor
   conversion. The off-diagonal elements had wrong imaginary part sign:
   - Before: V_{up,down} = 0.5*(V_x + i*V_y)  (wrong)
   - After:  V_{up,down} = 0.5*(V_x - i*V_y)  (correct, from sigma_y)

2. DFT+U force/stress (dftu_force_stress.hpp): Convert VU from Pauli basis
   to spinor basis before force calculation. The old code incorrectly mixed
   Pauli-basis VU with spinor-basis DM.

3. DeltaSpin force/stress (dspin_force_stress.hpp): Convert lambda from
   Pauli basis to spinor basis. The constraint force F = lambda·dM/dR
   requires proper Pauli-to-spinor conversion:
   - lambda_spinor = (lambda_z, lambda_x, lambda_x, -lambda_z)
   for (uu, ud, du, dd) components.

These fixes ensure consistent Pauli-to-spinor conversion across all modules:
- H construction (gint_common.cpp): already fixed
- DFT+U Hamiltonian: fixed in this commit
- DFT+U force/stress: fixed in this commit
- DeltaSpin force/stress: fixed in this commit

Verified by scf_u_spin4 test (nspin=4 + DFT+U): SCF converges correctly.

* chore: remove .py and .md test files from PR

* fix: correct DFT+U force for nspin=4 - DM is stored in Pauli basis, not spinor basis

Two bugs fixed in dftu_force_stress.hpp:
1. Removed incorrect VU Pauli-to-spinor conversion: DM for nspin=4 is stored in
   Pauli basis (rho_0, rho_x, rho_y, rho_z) per func_xyz_to_updown(), so VU must
   also stay in Pauli basis for the force trace formula F = -Tr(VU * dDM/dR).
2. Removed force *= 2.0 for nspin=4: Pauli basis already includes all spin
   channels, unlike nspin=1 where the factor of 2 accounts for spin degeneracy.
Updated scf_u_spin4 result.ref accordingly.

* fix: remove force*=2.0 for nspin=4 in DeltaSpin - Pauli basis already covers all spin channels

* fix: add missing blacs_context to ELPA Constructor 1 for nspin=4 support

Constructor 1 of ELPA_Solver was missing elpa_set_integer("blacs_context", ...)
while Constructor 2 (otherParameter) already had it. Without blacs_context,
ELPA's internal MPI operations (e.g. MPI_Bcast in complex Cholesky and
invert_triangular) can fail with INVALID DATATYPE when using complex eigensolves.
Also update scf_angle_spin4 result.ref with corrected reference energy.

* fix: correct rho_y sign in spinor-to-Pauli DM conversion (func_xyz_to_updown)

For Pauli decomposition: rho = rho_0*I + rho_x*sigma_x + rho_y*sigma_y + rho_z*sigma_z
sigma_y = [[0,-i],[i,0]], so rho_updown = rho_x + i*rho_y, rho_downup = rho_x - i*rho_y
Thus rho_y = Im(rho_updown - rho_downup) = tmp[1].imag() - tmp[2].imag].

Previously the real version had -tmp[1].imag()+tmp[2].imag() = -2*rho_y (wrong sign),
and the complex version had i*(tmp[1].imag()-tmp[2].imag()) = 2i*rho_y (wrong formula).
This broke rotational invariance: mag along y gave wrong energy (~4 eV deviation vs x/z).

* test: update scf_angle_spin4 and scf_u_spin4 result.ref after DM rho_y fix

* chore: revert density_matrix.cpp rho_y fix (wrong branch) and remove verify_hamiltonian_convention test dir

- Revert density_matrix.cpp func_xyz_to_updown rho_y sign fix from
  commit 52ee608 (belongs on a separate DM-fix branch)
- Remove tests/03_NAO_multik/verify_hamiltonian_convention/ (debug helper)
- Update result.ref for scf_angle_spin4 and scf_u_spin4 to match
  current code (Pauli-to-spinor + ELPA fixes only)

* fix: restore density_matrix.cpp rho_y sign fix (paired with gint_common clx_j fix)

The gint_common.cpp fix corrects Pauli→spinor (H construction) and the
density_matrix.cpp fix corrects spinor→Pauli (DM Fourier transform).
Both must use the same σ_y convention for self-consistency.
Also update result.ref files for both test cases.

* fix: correct DFT+U force/stress reference values and clean up empty nspin=4 block

Both VU (from cal_v_of_u) and DMR are stored in Pauli basis for nspin=4,
so the Pauli-to-spinor conversion in force/stress calculation is NOT needed.
The previous result.ref for scf_u_spin4 (totalforceref=6.562) was incorrect
because it was generated with code that mixed Pauli-basis VU with incorrectly
converted values. The correct force is 11.33, consistent with the physical
Pauli-basis trace Tr(VU * dDM/dR).

Changes:
- Remove empty if(nspin==4) block in dftu_force_stress.hpp (no conversion needed)
- Update scf_u_spin4/result.ref: totalforceref 6.562 -> 11.332 (correct value)
- Update scf_angle_spin4/result.ref: energy/stress to match computed values
- Add scf_angle_spin4/threshold: relax energy threshold to 1e-5 eV for
  non-collinear calculation numerical reproducibility
- Update scf_out_dos_spin4/result.ref: force/stress to match computed values

* fix: correct sigma_y sign convention in Pauli-spinor conversions for nspin=4

Fix inconsistent Pauli-to-spinor and spinor-to-Pauli conversion signs
across multiple modules, which broke rotational invariance in
non-collinear (nspin=4) LCAO calculations.

The standard sigma_y = [[0,-i],[i,0]] convention requires:
  H_{up,down} = B_x - i*B_y
  rho_y = -Im(rho_updown - rho_downup)

Changes:
- density_matrix.cpp: fix rho_y sign in func_xyz_to_updown
  (real: restore correct -Im(updown)+Im(downup);
   complex: use i*(updown-downup) with full complex values)
- dftu_pw.cpp: fix Pauli-to-spinor sign in DFT+U transfer_vu (PW path)
- spin_constrain.h: fix pauli_to_moment My sign, update comments
- dspin_lcao.cpp: fix cal_coeff_lambda Pauli-to-spinor sign
- deltaspin_core_test.cpp: update unit test for corrected formula
- scf_u_spin4/result.ref: update reference values after correction

* docs: fix My sign in spin_constrain.cpp comment

Correct the comment to match the sigma_y = [[0,-i],[i,0]] convention:
My = -Im(occ[1] - occ[2]) instead of Im(occ[1] - occ[2])

* test: update 099_PW_DJ_SO ref after sigma_y sign fix in dftu_pw

* test: update scf_out_dos_spin4 ref after sigma_y sign fix

* test: update scf_angle_spin4 ref after sigma_y sign fix

Author: dyzheng

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