Tests: add nspin=4 sigma_y regression case (scf_angle_spin4_y); #7759 is already fixed

[Critsium-xy]

Reminder

• I have read AGENTS.md and docs/developers_guide/agent_governance.md.
• I have linked an issue or explained why this PR does not need one.
• I have added adequate unit tests and/or case tests, or explained why not.
• I have listed the exact verification commands run and their results.
• I have described user-visible behavior changes, including INPUT parameter changes.
• I have explained core-module impact for ESolver, HSolver, ElecState, Hamilt, Operator, Psi, or other source/ changes.
• I have requested any needed governance exception below.

Linked Issue

Refs #7759 (and #7831). This PR contains no source change — please read the first section before merging.

First: #7759 is already fixed on develop

I set out to fix #7759 and found it had been fixed in the meantime, by #7832 (merged 2026-09-01), which is exactly "Option A" from the issue report. The history of the m_y convention:

PR	clx_j in merge_hr_part_to_hR	func_xyz_to_updown
before #7513	{0, 1, -1, 0}, no conj	tmp[1].imag() - tmp[2].imag()
#7513	{0, -1, 1, 0} + conj (textbook Pauli)	unchanged — inconsistent, this is the state #7759 reported
#7664	unchanged	-tmp[1].imag() + tmp[2].imag()
#7832	unchanged	tmp[1].imag() - tmp[2].imag() — consistent again

Verified empirically on develop 3937023d6 with the two-Fe non-collinear case, serial LAPACK solver:

case	E_KohnSham (eV)	AMAG	iters
scf_angle_spin4 (±x)	−6267.4651888506	6.70362	12
±y twin	−6267.4651888506	6.70362	12

Identical to 2.7e-9 eV, identical iteration count, identical DRHO trajectory, forces identically zero, stress agreeing to 4e-10 kbar. The σy loop on develop is self-consistent; #7759 can be closed. (Unrelated to #7913, which is about the DeltaSpin operator and the genelpa lower-triangle read, both still open.)

What's changed?

What is not fixed is that nothing in CI would have noticed any of this — so this PR closes that gap instead.

scf_angle_spin4 initialises the two Fe moments along ±x. m_x and m_z are read off the real part of the spinor density-matrix block and are invariant under the conjugation convention of that block; only m_y reads the imaginary part. The existing case is therefore structurally incapable of seeing a sign error in the σy channel — which is how this convention managed to change three times without a single test reacting.

tests/03_NAO_multik/scf_angle_spin4_y is the same case with the moments globally rotated by 90° about z (±x → ±y). Without spin–orbit coupling the spin frame is decoupled from the lattice, so every reference value must equal the scf_angle_spin4 one exactly. result.ref is that file byte for byte — the assertion is the invariant, rather than a freshly regenerated number. The anisotropic 1×2×1 k-mesh does not weaken this: at lspinorb 0 a global spin rotation cannot couple to the k-mesh.

Unit Tests and/or Case Tests for my changes

• Build: Linux, GCC, OpenMPI, cmake -B build -G Ninja -DBUILD_TESTING=ON -DENABLE_LCAO=ON -DENABLE_MPI=ON -DENABLE_OPENMP=ON; cmake --build build -j15 → 3318/3318, 0 errors.

• Ran both cases from their own directory with OMP_NUM_THREADS=1 ./abacus_basic_para → the table above.

• Teeth check. I reintroduced the Fix: correct Pauli-to-spinor Hamiltonian conversion for nspin=4  #7664 m_y sign in func_xyz_to_updown, rebuilt, and reran both cases:

  case	E_KohnSham (eV)	AMAG	iters
  scf_angle_spin4 (±x)	−6267.4651888506	6.70362	12 — unchanged, detects nothing
  scf_angle_spin4_y (±y)	−6264.2726978309	3.9e-05	50, never converged

  3.19 eV off, moment gone, no convergence. The new case catches the regression that the old one is blind to. The source change was then reverted; this PR ships only the test.

• Checks not run, with reason: I could not run the case under the CI solver. The box available to me links a gfortran-built ScaLAPACK against MKL, and pzpotrf → zdotc segfaults there for any rank count (an ABI mismatch in that environment, not an ABACUS defect), so both runs above use ks_solver lapack in serial. This is why result.ref is copied from the CI-generated scf_angle_spin4 reference rather than regenerated locally: the invariant is exactness against that file, and my serial numbers sit 7.7e-7 eV from it, well inside the case's threshold 0.00001. CI running this case for the first time is the real check.

Governance Notes

• INPUT/docs changes: none — no source/ change at all, only a new integration case plus its registration in CASES_CPU.txt. The governance check reports no findings.
• Core module impact: none.
• Coverage gap I am deliberately leaving open: scf_angle_spin4 is also listed in CASES_GPU.txt. I did not register the y twin there, because I had no GPU to verify the nspin=4 LCAO GPU path against. If a maintainer with GPU access confirms it, adding one line to CASES_GPU.txt extends the same guard to that path, and is worth doing.
• Exceptions requested: none.

🤖 Generated with Claude Code

[Critsium-xy]

Update after the first CI run of this case, which failed — correctly.

scf_angle_spin4_y: etotref        cal=-6269.25160172 ref=-6267.46518962 dev= 1.78641210
scf_angle_spin4_y: totalstressref cal= 4089.31939500 ref= 3912.92043700 dev=-176.39895800

I had verified the invariance only with the serial LAPACK solver (my box could not link ScaLAPACK against MKL). Working around that with LD_PRELOAD of OpenBLAS made the ScaLAPACK and ELPA paths usable, and the solver matrix is unambiguous:

case	solver	np	E_KohnSham (eV)	AMAG	iters
±x	scalapack_gvx	2	−6267.4651888505	6.70362	12
±y	scalapack_gvx	2	−6267.4651888505	6.70362	12
±x	scalapack_gvx	4	−6267.4651888505	6.70362	12
±y	scalapack_gvx	4	−6267.4651888505	6.70362	12
±x	genelpa	2	−6267.4651898827	6.70362	12
±y	genelpa	2	−6269.6774946269	6.80171	50, no conv.
±x	genelpa	4	−6267.4651902210	6.70362	12
±y	genelpa	4	−6267.7588409455	6.76183	50, no conv.

With scalapack_gvx the rotated pair is bitwise identical at both rank counts, so the σy convention this case exists to pin is correct and #7759 remains closed-out. With genelpa the ±y case does not converge and settles below the variational minimum by a rank-dependent amount — that is a solver defect, and it matches defect 2 of #7913 (the complex generalized transform reading the stale lower triangle of a matrix whose lower half LCAO never fills). Details and reasoning posted on #7913.

So the latest commit pins ks_solver scalapack_gvx, as scf_deltaspin4 and nine other cases in this directory already do, and the README records the measurements and the condition for removing the pin. Re-verified with the real harness rather than by hand this time:

cd tests/03_NAO_multik
../integrate/Autotest.sh -a <bin> -n 4 -r scf_angle_spin4
-> [ PASSED ] 8 test cases passed.     # 4 checks x 2 cases

Net effect: the case is green, it still pins the σy convention (reintroducing the #7664 sign still moves it by 3.19 eV), and once #7913 lands the pin can be dropped so it guards the genelpa path too.

[AsTonyshment]

Sry but... just a friendly heads-up: these convention issues are quite interconnected, so I would be cautious about treating this test as proof that the broader problem is resolved.

I've explored this on my local experimental branch. A more systematic cleanup may involve making DMK store the physical density matrix $P_{\mu\nu}(k) = \sum_n f_{nk} c_{\mu nk} c_{\nu nk}^{*}$ directly (currently it stores $\sum_n f_{nk} c_{\mu nk}^{*} c_{\nu nk}$), together with consistent changes to matrix indexing, the DMK-to-DMR Fourier phase (should be $\mathrm{e}^{-\mathrm{i}k\cdot R}$, currently it's $\mathrm{e}^{\mathrm{i}k\cdot R}$), Gint's lattice-vector convention, and the nspin=4 Pauli decomposition. For example, for a physical spin-density block $\Gamma$, we have $m_y = \mathrm{i}(\Gamma_{\uparrow\downarrow} - \Gamma_{\downarrow\uparrow})$; extracting the same physical component from $\Gamma^*$ requires the opposite sign. So an individual sign cannot be judged independently of what the input matrix actually represents.

These issues may need a thorough refactor by experts (like @dyzheng, he is very experienced😄), something like #7833, though it was closed due to out of sync :(

This matters beyond the SCF test, the conventions also enter Mulliken/DeltaSpin, DFT+U, forces/stress, and RT-TDDFT current contractions, where conjugation and index ordering can affect even the sign of the RT-TDDFT current. That's why I care about consistency across the code.

I'm not implying to block this PR, just some friendly comments. If it is merged, I would expect its reference values may need updating after independently validated fixes in the future. The no-SOC rotation invariant should still hold, but the current numbers aren't necessarily the final correct reference.