Fix Pauli-to-Spinor Conversion in LCAO Non-Collinear Calculations

[dyzheng]

Summary

This PR fixes critical bugs in the Pauli-to-spinor basis conversion for LCAO non-collinear spin (nspin=4) calculations. The fixes ensure correct Hamiltonian construction, DFT+U potential handling, and force/stress calculations for systems with non-collinear magnetization and spin-orbit coupling (SOC).

Background

In ABACUS LCAO calculations with nspin=4 (non-collinear spin), the effective potential and other operators are initially computed in the Pauli basis (V₀, V_x, V_y, V_z) and must be converted to the spinor basis for the 2×2 spinor Hamiltonian matrix. The conversion follows:

H = V₀·I + V_x·σ_x + V_y·σ_y + V_z·σ_z

where the Pauli matrices are:

σ_x = [[0, 1], [1, 0]]
σ_y = [[0, -i], [i, 0]]
σ_z = [[1, 0], [0, -1]]

This yields the spinor basis representation:

H_{↑,↑}   = V₀ + V_z                  (real)
H_{↑,↓}   = V_x - i·V_y               (complex, imaginary part negative)
H_{↓,↑}   = V_x + i·V_y               (complex, imaginary part positive)
H_{↓,↓}   = V₀ - V_z                  (real)

Issues Fixed

Issue 1: Incorrect Sign in Off-Diagonal Hamiltonian Elements

File: source/source_lcao/module_gint/gint_common.cpp

Problem: The merge_hr_part_to_hR() function used incorrect signs for the imaginary components when converting from Pauli to spinor basis.

Original Code:

std::vector<int> clx_j = {0, 1, -1, 0};
// is=1 (up,down): +i  → H_{up,down} = V_x + i·V_y  ❌ WRONG
// is=2 (down,up): -i  → H_{down,up} = V_x - i·V_y  ❌ WRONG

Fixed Code:

std::vector<int> clx_j = {0, -1, 1, 0};
// is=1 (up,down): -i  → H_{up,down} = V_x - i·V_y  ✅ CORRECT
// is=2 (down,up): +i  → H_{down,up} = V_x + i·V_y  ✅ CORRECT

Additional Fix: Added complex conjugate transpose for lower triangle filling to ensure Hermiticity:

// Before: only transpose
lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol);

// After: conjugate transpose
lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol));

This ensures H(-R) = H†(R), which is required for the Hamiltonian to be Hermitian in k-space.

Issue 2: Incorrect Fourier Transform Sign in DM(R) Calculation

Files:

• source/source_estate/module_dm/density_matrix.cpp
• source/source_lcao/module_lr/dm_trans/dmr_complex.cpp

Problem: The inverse Fourier transform from k-space to real-space used e^{+ik·R} instead of the correct e^{-ik·R}.

Mathematical Background:
The Fourier transform pair should be:

• Forward (R→k): H(k) = Σ_R H(R) · e^{+ik·R}
• Inverse (k→R): DM(R) = (1/N_k) Σ_k DM(k) · e^{-ik·R}

Fix: Changed the sign in the phase factor:

// Before:
kphase_vec[ik][iR] = TK(cosp, sinp);   // e^{+ik·R} ❌

// After:
kphase_vec[ik][iR] = TK(cosp, -sinp);  // e^{-ik·R} ✅

Issue 3: Incorrect Pauli-to-Spinor Conversion in DFT+U

File: source/source_lcao/module_operator_lcao/dftu_lcao.cpp

Problem: The transfer_vu() function had the same sign error as Issue 1 when converting the DFT+U potential from Pauli to spinor basis.

Fix:

// Before:
vu[index[1]] = 0.5 * (vu_tmp[index[1]] + i * vu_tmp[index[2]]);  // V_x + i·V_y ❌
vu[index[2]] = 0.5 * (vu_tmp[index[1]] - i * vu_tmp[index[2]]);  // V_x - i·V_y ❌

// After:
vu[index[1]] = 0.5 * (vu_tmp[index[1]] - i * vu_tmp[index[2]]);  // V_x - i·V_y ✅
vu[index[2]] = 0.5 * (vu_tmp[index[1]] + i * vu_tmp[index[2]]);  // V_x + i·V_y ✅

Issue 4: Inconsistent Basis in DFT+U Force/Stress Calculation

File: source/source_lcao/module_operator_lcao/dftu_force_stress.hpp

Problem: The force/stress calculation mixed Pauli-basis VU with spinor-basis DM, computing an incorrect trace.

Mathematical Background:
The force should be computed as:

F = -Tr(VU · dDM/dR)

In the spinor basis, this expands to:

F = V_{uu}·dDM_{uu} + V_{ud}·dDM_{du} + V_{du}·dDM_{ud} + V_{dd}·dDM_{dd}

Where:

V_{uu} = 0.5·(V₀ + V_z)
V_{dd} = 0.5·(V₀ - V_z)
V_{ud} = 0.5·(V_x - i·V_y)
V_{du} = 0.5·(V_x + i·V_y)

For real DM (dDM_{ud} = dDM_{du}), this simplifies to:

F = 0.5·(V₀+V_z)·dDM_{uu} + V_x·dDM_{ud} + 0.5·(V₀-V_z)·dDM_{dd}

Fix: Convert VU from Pauli to spinor basis before force calculation:

if(this->nspin == 4) 
{
    const int m_size2 = tlp1 * tlp1;
    std::vector<double> VU_pauli = VU;
    for (int m = 0; m < m_size2; m++)
    {
        VU[m]                  = 0.5 * (VU_pauli[m] + VU_pauli[m + 3 * m_size2]); // V_uu
        VU[m + m_size2]        = 0.5 * VU_pauli[m + m_size2];                     // Re(V_ud)
        VU[m + 2 * m_size2]    = 0.5 * VU_pauli[m + m_size2];                     // Re(V_du)
        VU[m + 3 * m_size2]    = 0.5 * (VU_pauli[m] - VU_pauli[m + 3 * m_size2]); // V_dd
    }
}

Issue 5: Inconsistent Basis in DeltaSpin Force/Stress Calculation

File: source/source_lcao/module_operator_lcao/dspin_force_stress.hpp

Problem: Similar to Issue 4, the DeltaSpin constraint force calculation used lambda (Pauli basis) directly with DM (spinor basis).

Fix: Convert lambda from Pauli to spinor basis:

if (nspin == 4)
{
    // lambda in spinor basis: (lambda_uu, lambda_ud, lambda_du, lambda_dd)
    const double lambda_spinor[4] = {lambda[2], lambda[0], lambda[0], -lambda[2]};
    for (int is = 0; is < 4; is++)
    {
        const double lambda_tmp = lambda_spinor[is];
        // ... force calculation using lambda_tmp and DM[step_trace[is]]
    }
}

Verification

Code Review

Comprehensive code review was performed on all modules handling nspin=4:

Module	File	Status	Notes
H construction	gint_common.cpp	✅ Fixed	Pauli→spinor conversion
DM Fourier	density_matrix.cpp, dmr_complex.cpp	✅ Fixed	k→R transform sign
DFT+U transfer_vu	dftu_lcao.cpp	✅ Fixed	Pauli→spinor conversion
DFT+U force/stress	dftu_force_stress.hpp	✅ Fixed	Basis consistency
DeltaSpin force/stress	dspin_force_stress.hpp	✅ Fixed	Basis consistency
SOC matrix elements	soc.cpp	✅ Correct	Uses Clebsch-Gordan coefficients
Nonlocal force/stress	nonlocal_force_stress.hpp	✅ Correct	Spinor basis consistent
DM construction	cal_dm_psi.cpp	✅ Correct	Standard DM = C†·f·C

Mathematical Verification

Hermiticity Check: H(R=0) satisfies H = H† to machine precision (max deviation: 1.00e-12)

Phase Convention: For magnetization along +y direction:

• Before fix: Im(H_{↑,↓}) > 0 ❌
• After fix: Im(H_{↑,↓}) < 0 ✅

This matches the correct physics: H_{↑,↓} = V_x - i·V_y has negative imaginary part when V_y > 0.

Test Results

Test Suite: tests/03_NAO_multik/ (nspin=4 LCAO tests)

Test Case	Description	Status	Notes
scf_out_dos_spin4	DOS output with SOC	✅ PASS	All 5 checks passed
scf_out_mul_spin4	Mulliken analysis	✅ PASS	All 6 checks passed
scf_u_spin4	DFT+U with SOC	✅ PASS	Reference updated
scf_angle_spin4	Angle-constrained magnetization	⚠️ MPI error	Pre-existing bug
scf_out_hsr_spin4	H(R)/S(R) output	⚠️ MPI error	Pre-existing bug

Note: Two tests (scf_angle_spin4 and scf_out_hsr_spin4) fail with MPI_Bcast errors. Investigation shows this is a pre-existing bug in the develop branch, unrelated to the Pauli matrix fixes. These tests also fail on the unmodified develop branch.

Reference Value Updates

The scf_u_spin4 test reference values were updated to reflect the corrected physics:

Quantity	Old Value	New Value	Change
Total energy (eV)	-6789.2818	-6789.1424	+0.1394 eV
Force	11.335	14.446	+3.111
Stress (kbar)	4892.275	4327.911	-564.364

These changes are expected and correct, as the old reference values were computed with buggy code.

Impact Assessment

Affected Calculations

• ✅ nspin=4 LCAO calculations (non-collinear spin, SOC)
• ✅ DFT+U with nspin=4
• ✅ DeltaSpin constraints with nspin=4
• ✅ Force/stress calculations with nspin=4

Unaffected Calculations

• ✅ nspin=1,2 calculations (collinear spin) - no Pauli matrix conversion needed
• ✅ PW basis calculations - different code path
• ✅ Gamma-only calculations - no k-point Fourier transform

Backward Compatibility

Breaking Change: Results for nspin=4 LCAO calculations will change. This is intentional and correct - the old results were wrong due to the bugs fixed in this PR.

Users should:

1. Re-run nspin=4 calculations with this fix
2. Expect different (correct) energies, forces, and stresses
3. Update any reference data or workflows that depended on the old (incorrect) results

Files Modified

Source Code

1. source/source_lcao/module_gint/gint_common.cpp - H construction fix
2. source/source_estate/module_dm/density_matrix.cpp - DM Fourier fix
3. source/source_lcao/module_lr/dm_trans/dmr_complex.cpp - DM Fourier fix
4. source/source_lcao/module_operator_lcao/dftu_lcao.cpp - DFT+U transfer_vu fix
5. source/source_lcao/module_operator_lcao/dftu_force_stress.hpp - DFT+U force fix
6. source/source_lcao/module_operator_lcao/dspin_force_stress.hpp - DeltaSpin force fix

Tests

7. tests/03_NAO_multik/scf_u_spin4/result.ref - Updated reference values

Testing Instructions

To verify the fixes:

# Build ABACUS
cd build
cmake ..
make -j

# Run nspin=4 tests
cd tests/03_NAO_multik
export OMP_NUM_THREADS=2
bash ../integrate/Autotest.sh -a /path/to/abacus -n 4 -r "scf_out_dos_spin4|scf_out_mul_spin4|scf_u_spin4"

Expected output:

[PASSED] 15 test cases passed.

Related Issues

This PR addresses the Pauli-to-spinor conversion issues identified in the code review. The MPI errors in scf_angle_spin4 and scf_out_hsr_spin4 are separate pre-existing issues that require dedicated investigation.

Checklist

• Code changes implement correct Pauli-to-spinor conversion
• All mathematical formulas verified
• Comprehensive code review of all nspin=4 modules
• Test suite passes (excluding pre-existing MPI bugs)
• Reference values updated for affected tests
• Documentation updated with detailed explanations
• No impact on nspin=1,2 or PW calculations
• Changes pushed to feature branch