Fix and extend DFT+U for PW basis (nspin=1/2/4), occupation matrix mixing, and onsite kernel corrections

[dyzheng]

Describe the bug

Background

DFT+U in ABACUS was previously only available for the plane-wave (PW) basis with nspin=4. While porting DFT+U to PW nspin=1/2 and reviewing the implementation, several correctness bugs in the onsite force/stress kernels and MPI parallelization were found. These fixes are independent of the new feature work and should land first.

Bugs to fix

1. Wrong Pauli matrix contraction order in nonlocal onsite force/stress for nspin=4/SOC.
   The contraction ps[1]*dbb2 + ps[2]*dbb1 vs ps[1]*dbb1 + ps[2]*dbb2 must match the storage convention of the coefficient array (deeq_nc uses index 1 = σ_↓↑, index 2 = σ_↑↓; the DFT+U vu array uses the opposite convention; the DeltaSpin lambda_coeff array follows the deeq_nc convention). The formula must be chosen per coefficient-array convention, not by a universal rule. Affects CPU (force_op.cpp, stress_op.cpp), CUDA (force_op.cu, stress_op.cu), and ROCm (force_op.hip.cu, stress_op.hip.cu) kernels.

2. nspin=2 crash in onsite_op. The npol==1/2 branches were missing, causing crashes for nspin=2 DFT+U/DeltaSpin calculations. Restore the npol branching in all kernel paths (CPU/CUDA/ROCm).

3. DFT+U locale double-counting when kpar > 1. The occupation matrix (locale) is accumulated multiple times across k-point parallel pools.

4. Incorrect spin-channel selection in cal_occ_pw. Use isk[ik] instead of ik >= nk/2 to select the spin channel.

5. GPU path leading dimension. Use vkbnc instead of vkb.nc as the leading dimension in the GPU kernel launch.

Feature work bundled with the fix

• Relax DFT+U nspin validation to accept nspin=1/2/4 for the PW basis.
• Encapsulate Plus_U with typed accessors (get/set_locale, get_orbital_corr, get_hubbard_u, is_locale_initialized, mark_locale_dirty, enable_mixing) and replace direct member accesses in LCAO operator code and DeePKS.
• Rewrite cal_occ_pw for nspin=1/2/4 with correct becp indexing.
• Restructure eff_pot_pw layout: nspin=2 uses split [spin_up|spin_down]; add get_eff_pot_pw_spin(isk) for nspin-aware access.
• Add DFT+U occupation matrix mixing via the Broyden method in Charge_Mixing (uom_mdata, allocate_mixing_uom(), mix_uom(), conserve_setting()), allocated at the first SCF iteration for PW.
• Propagate ld_psi for correct GEMM strides when ngk[ik] < npwx.
• Extract setup_pw_dftu_indices() from cal_ps_dftu.

Scope

• source/source_lcao/module_dftu/ (dftu.cpp, dftu_io.cpp, dftu_occup.cpp, dftu_pw.cpp)
• source/source_estate/module_charge/charge_mixing.{h,cpp}
• source/source_pw/module_pwdft/ (onsite_proj*, forces_onsite.cpp, stress_onsite.cpp, kernels/force_op., kernels/stress_op.)
• source/source_lcao/module_operator_lcao/dftu_lcao.cpp
• source/source_lcao/module_deepks/ (accessor migration only)
• Unit tests: source/source_lcao/module_dftu/test/dftu_pw_test.cpp
• Integration tests: DFTU-only cases under tests/17_DS_DFTU/
• Docs: construct_H.md, dft_plus_u entry in input-main.md

Testing

• PW DFT+U unit tests for nspin=1/2/4.
• Integration tests: DFT+U-only cases (PW/LCAO, nspin=1/2/4), FeO atom-order invariance cases (50/51), SOC force/stress regression (035_PW_15_SO, 099_PW_DJ_SO).
• Verify kpar>1 runs produce identical occupation matrices.

Related

Split out from PR #7693.

Expected behavior

No response

To Reproduce

No response

Environment

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

• Verify the issue is not a duplicate.
• Describe the bug.
• Steps to reproduce.
• Expected behavior.
• Error message.
• Environment details.
• Additional context.
• Assign a priority level (low, medium, high, urgent).
• Assign the issue to a team member.
• Label the issue with relevant tags.
• Identify possible related issues.
• Create a unit test or automated test to reproduce the bug (if applicable).
• Fix the bug.
• Test the fix.
• Update documentation (if necessary).
• Close the issue and inform the reporter (if applicable).