fix: restore product-rule factor in LCAO DeltaSpin forces

[aboys-cb]

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

No separate issue: this PR documents and fixes the force-normalization regression introduced by #7513 (698e1876211b9fdf4f334fa8ac3efa40b9e4801f, 2026-06-30). Base: upstream develop at 97697c26b4562faf7fbcfffd75ba970c635d945a; no downstream optimization commits are included.

Unit Tests and/or Case Tests for my changes

Update the existing tests/03_NAO_multik/scf_deltaspin4 case. The regular 03_NAO_multik CI job already selects this case through CASES_CPU.txt. There is no new case directory, orbital asset, Python test runner, shared comparison logic, or CTest entry.

The case uses the existing Fe.upf and Fe_gga_6au_100Ry_4s2p2d1f.orb from tests/PP_ORB. It is periodic Fe2 in a 2.86 Angstrom cubic cell, LCAO + SOC (nspin 4), 12 Ry, a 2x2x2 k mesh, no Hubbard U, and target moments (0,0,3.0) and (2.59807621135332,0,1.5) Bohr magnetons: both have magnitude 3.0, with a 60-degree relative tilt. Fe1 is at the origin; Fe2 is at (1.573,1.287,1.430) Angstrom, displaced by (+0.143,-0.143,0) from the body center. The geometric displacement produces nonzero atomic forces, while the magnetic magnitude and direction perturbation exercise the DeltaSpin constraint-force contribution.

The old case deliberately accepted an unconverged 100-step SCF with tolerances of 1 eV, 10 eV/Angstrom and 500 kbar. Those thresholds could hide the force error. The lightweight CI case uses ecutwfc 12 and scf_nmax 10 to bound its cost. SCF convergence and basis/grid convergence are not acceptance criteria for this output regression. Its regenerated result.ref checks reproducibility at those exact settings, with the force threshold kept at 1e-6 eV/Angstrom. Energy/stress use the standard 1e-7 eV / 1e-3 kbar defaults. The short, extensionless README has two lines.

The unchanged Autotest.sh reads reference values from the case's result.ref. Its totalforceref is the sum of absolute values of the final printed atomic-force components, not the net vector force. This is a standard output regression, not a finite-difference calculation inside CI. Independent derivative evidence is reported separately below.

Current lightweight regression: 12 Ry, at most 10 SCF steps

Per review, restore the original benchmark cutoff (ecutwfc 100 -> 12) and reduce scf_nmax 240 -> 10 to bound runtime. All other computational inputs, assets, comparison tolerances and shared test scripts are unchanged in this follow-up. The README still has two lines. The initial 12 Ry / 240-step probe was stopped after observing charge oscillation; retaining 240 steps would defeat the low-cost regression purpose.

Job 1730113, Sai 16v100n08, 16V100 / flood-1o2gpu, completed successfully. The CPU scalapack_gvx one-rank reference and independent four-rank corrected run both completed exactly 10 steps without charge convergence, as allowed for this test. OMP/OpenBLAS=1. The checked fixed/original executable SHA256 values are the same as those recorded below. The shared Autotest.sh commands below were rerun unchanged at the new input settings.

Assertion	Previous 100 Ry reference	New 12 Ry / 10-step reference	Unchanged tolerance
etotref (eV)	-6778.9302831721670373	-6807.4629648003083275	1e-7
etotperatomref (eV)	-3389.4651415861	-3403.7314824002	1e-7
totalforceref (eV/Angstrom)	7.902208	93.685868	1e-6
totalstressref (kbar)	861.017880	30117.121565	1e-3

These reference changes follow the deliberately different cutoff and finite SCF trajectory; they are not converged physical predictions or a force-only code change in SCF energy. totaltimeref is refreshed from 142.51 to 27.24 seconds and is not a numerical acceptance assertion. Observed runtimes were 27.24 seconds (one rank) and 29.70 seconds (four ranks), single observations on this node, not a controlled performance speedup claim.

The four-rank corrected run passed every standard assertion: energy differed from the one-rank reference by 1.34e-9 eV, aggregate force was identical at six decimals, and aggregate stress differed by 2e-6 kbar. The original buggy executable at four ranks reproduced the corrected energy and stress exactly but gave totalforceref=94.929118; its force error against the new reference is 1.243250 eV/Angstrom, and the standard runner returned 1 (expected force-only failure). Thus lowering the cutoff and limiting the iteration count preserve detection of the missing factor of two, without relaxing tolerances.

# Fresh case copies, ecutwfc=12 and scf_nmax=10:
bash ../integrate/Autotest.sh -a "$FIXED" -n 1 -o 1 -j 1 -r '^scf_deltaspin4$' -g
# Copy generated result.ref into independent fixed/original case directories.
bash ../integrate/Autotest.sh -a "$FIXED" -n 4 -o 1 -j 1 -r '^scf_deltaspin4$'     # PASS
bash ../integrate/Autotest.sh -a "$ORIGINAL" -n 4 -o 1 -j 1 -r '^scf_deltaspin4$'  # expected force-only FAIL

The earlier converged 100 Ry controls and independent derivative evidence below remain separate. No new 12 Ry GPU or finite-difference claim is made. git diff --check and the staged governance checker pass for this follow-up.

Earlier 100 Ry converged negative control (separate physical validation)

Slurm job 1726202, node 16v100n01: all SCFs passed the convergence audit. The one-rank reference converged in 37 steps; both four-rank runs converged in 35 steps.

CPU, 4 MPI ranks	Original condition	Fixed
Fe1 x final printed force (eV/Angstrom)	1.9455369007	1.9745000028
Fe1 x DeltaSpin force (eV/Angstrom)	0.0289631020	0.0579262041
Fe1 y DeltaSpin force (eV/Angstrom)	-0.0288718778	-0.0577437555
Fe1 z DeltaSpin force (eV/Angstrom)	-0.0000227294	-0.0000454589
totalforceref (eV/Angstrom)	7.786582	7.902208
Total energy (eV)	-6778.930283172707	-6778.930283172707
totalstressref (kbar)	861.017880	861.017880
SCF steps	35	35
Final DRHO	6.1093e-11	6.1093e-11
Final constrained-moment RMS	8.4338e-11	8.4338e-11
Final DeltaE_womix (eV)	-7.53286e-10	-7.53286e-10
Standard Autotest exit code	1 (expected FAIL)	0 (PASS)

The independently cold-started one-rank reference gives totalforceref=7.902208, identical to the four-rank fixed result at the collector's six-decimal precision. Its total energy differs from the four-rank result by 5.40e-10 eV. The original-condition run fails on force only, with an aggregate difference of 0.115626 eV/Angstrom (115,626 times the force threshold); energy and stress pass. Thus restoring the bug is demonstrably detected by the existing shared comparison flow.

The stronger constraint was chosen using a controlled magnetic perturbation: at the same geometry and 60-degree tilt, increasing both target magnitudes from 2.2 to 3.0 Bohr magnetons increases Fe1's DeltaSpin x force from 0.0084497251 to 0.0579262041 eV/Angstrom (6.86x). The force contribution itself is measured here, rather than assuming a larger total atomic force implies a stronger DeltaSpin signal.

Executed from the staged tests/03_NAO_multik directory, with the repository's unmodified integration scripts and the final case name:

bash ../integrate/Autotest.sh -a "$FIXED" -n 1 -o 1 -j 1 -r '^scf_deltaspin4$' -g
# Copy that generated result.ref into fresh fixed/original case directories.
bash ../integrate/Autotest.sh -a "$FIXED" -n 4 -o 1 -j 1 -r '^scf_deltaspin4$'
bash ../integrate/Autotest.sh -a "$ORIGINAL" -n 4 -o 1 -j 1 -r '^scf_deltaspin4$'

Reference generation uses the fixed executable at one rank; validation uses a fresh cold start at four ranks. The original executable is a matched negative control with upstream dspin_fs.cpp. Both use the unmodified upstream AO implementation and the same base objects and link settings. The check also audits the SCF-converged marker, DRHO < 1e-10, constrained-moment RMS < 1.01e-10 (printed rounding), and |DeltaE_womix| < 1e-9 eV. Those extra convergence audits are external validation, not new logic in the shared CI script.

Environment: Sai 16V100 / flood-1o2gpu, CPU scalapack_gvx, one or four MPI ranks, OMP/OpenBLAS=1; GCC 13.3.0, CUDA 12.9, OpenMPI 5.0.10 and Libxc 7.0.0. Reused the previously built/verified v3.11.0-beta10 executables because this follow-up changes only case data/layout. No performance comparison is asserted.

Fixed executable SHA256: 97f9ec738013190ff33e0d4012b56737fcce8f9099469704f5e349c52da6dbee.
Original executable SHA256: 194b6e94a96df0ee9fe51b3b42964bdfe51136ddac30c6bc5f073bb714c942a3.

Independent derivative at 100 Ry (separate physical validation)

Slurm job 1726205, node 16v100n01, one V100 GPU / one MPI rank, OMP/OpenBLAS=1, cusolver. This diagnostic uses the same geometry, orbital, k mesh and target moments, but retains the earlier 100 Ry cutoff and converged SCF. It is independent of the current 12 Ry / 10-step CI case; the device/solver and asset paths also differ. A cold-start SCF converged in 35 steps (DRHO=6.1242e-11, moment RMS=9.5650e-11, DeltaE_womix=-2.5522e-11 eV).

For Fe1 x, the DeltaSpin contribution is:

Method	Force (eV/Angstrom)
Explicit analytical derivative of both overlaps	0.057926204064935
Fixed production expression in the diagnostic	0.057926204064948
Frozen-trace central difference, h=1e-4 Bohr	0.057926205940877
Frozen-trace central difference, h=2e-4 Bohr	0.057926197248923
Original production expression, CPU A/B above	0.0289631020

The finite-difference errors relative to the fixed expression are 1.88e-9 and 6.82e-9 eV/Angstrom. The two explicit analytical overlap terms and the fixed production expression agree within 1.4e-14 eV/Angstrom. This independently supports the product-rule normalization rather than merely observing that the new output is twice the old output.

# Run from a fresh copy of the final case with device gpu / ks_solver cusolver.
# PROBE is the separately instrumented diagnostic, not the production executable.
OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 mpirun -np 1 "$PROBE" > stdout.log 2> stderr.log

The diagnostic executable was checked against SHA256 2ee90c3678159f9d676870b8e425bbd467990f37518efd5519ca100640f19ab2 before execution. It was built from the same fixed source with only a separately instrumented force translation unit.

The separate diagnostic holds the converged density matrix D and constraint field lambda fixed, recomputes the two projector overlaps at displaced coordinates, and central-differences their contracted trace. It also evaluates both analytical overlap derivatives explicitly. It is not part of the generic test runner and introduces no diagnostic code into this PR. The CI result.ref is an output regression reference at 12 Ry / 10 steps; its values are not finite-difference reference values. The converged 100 Ry derivative evidence remains separate.

Local layout and governance checks

# A small registration harness includes the actual tests/03_NAO_multik/CMakeLists.txt.
cmake -S "$CHECK/ctest-src" -B "$CHECK/ctest-build"
ctest --test-dir "$CHECK/ctest-build" --show-only=json-v1
git diff --cached --check
python3 tools/03_code_analysis/agent_governance_check.py --staged --event-path "$PR_EVENT" --format text
python3 tools/03_code_analysis/agent_governance_check.py \
  --base 97697c26b4562faf7fbcfffd75ba970c635d945a --head HEAD \
  --event-path "$PR_EVENT" --format text

Registration-only CMake/CTest inspection passed: the existing 03_NAO_multik test invokes Autotest.sh, and its unchanged CPU case list contains scf_deltaspin4. Byte comparisons confirm the shared integration scripts, case lists and category CMake files are unchanged from the PR base. The README has exactly two lines. For that earlier 100 Ry validation, INPUT/STRU/KPT/threshold matched the tested files byte-for-byte. The current 12 Ry / 10-step validation is recorded above. Whitespace, staged governance, and full base-to-head governance checks passed with no findings. This registration harness does not claim a full repository build or full category runtime test.

Verification limits and negative results

In earlier validation on the smaller-displacement geometry (not the updated regression above), full-SCF total-energy differences at h=0.0005 Angstrom did not pass a 1e-4 eV/Angstrom force tolerance on the fixed upstream-based executable. The mean net-force subtraction was undone when comparing an individual atom to its energy derivative:

Trial	Analytical force	Energy derivative	Absolute error (eV/Angstrom)
Fe2 x, CPU / 4 MPI, 8-Bohr orbital	-0.169318452488	-0.183493657460	0.014175205
Fe2 x, GPU / 1 MPI, 8-Bohr orbital	-0.169318452426	-0.183493579243	0.014175127
Relative x displacement, GPU, 8-Bohr orbital	0.158986329600	0.166849738889	0.007863409
Relative x displacement, GPU, existing 6-Bohr orbital	0.324536400400	0.340216838595	0.015680438

For the relative displacement, R1x += q/2, R2x -= q/2, and the analytical derivative is (F1x-F2x)/2. For Fe2-only displacement the raw force is F2x_printed + net_force_x/2.

All SCFs met the electronic/magnetic convergence gates. The CPU/GPU agreement rules out calling this a GPU-only residual. These full-force discrepancies are not claimed fixed by this PR. The total-force regression above detects changes in the final printed force, while the frozen-state derivative isolates the diagnosed DeltaSpin error. Neither establishes full total-energy/force consistency for production relaxation or MD. The existing loose thresholds are tightened, not relaxed.

Full repository CI, PW finite differences, a broader scan of noncollinear directions, stress finite differences and ROCm execution were not run. This is a correctness fix, with no performance claim.

What's changed?

For basis_type lcao, sc_mag_switch 1 and nspin 4, the DeltaSpin atomic-force contribution was half of its complete value. Total forces are affected whenever this contribution is nonzero; the entire total force is not halved.

At fixed D and lambda, the relevant trace contains products of orbital-projector overlaps:

E_lambda = sum_(mu,nu,a) lambda_a D^a_(mu,nu) B_mu B_nu

d(B_mu B_nu)/dR = (dB_mu/dR) B_nu + B_mu (dB_nu/dR)

cal_force_IJR accumulates only the first derivative. To derive the missing factor, write the contracted trace as E = Re sum C_mu,nu B_mu* B_nu, with C_nu,mu = C_mu,nu* (the spin/constraint indices are implicit). Its derivative is Re(S1 + S2), where S1 = sum C_mu,nu (dB_mu*) B_nu and S2 = sum C_mu,nu B_mu* (dB_nu). Exchanging mu and nu in S2 and using Hermiticity gives S2 = S1*. Hence dE/dR = 2 Re(S1), and F = -2 Re(S1). Both terms survive in the Pauli representation as well. The final factor of two accounts for the two overlap derivatives, not spin degeneracy. Transforming D to Pauli components does not remove either derivative.

#7513 restricted the final multiplication to nspin != 4, reasoning that the Pauli basis already covers the spin channels. This conflated spin counting with the product rule. This PR restores unconditional multiplication by two and explains why in a code comment. The other Pauli-to-spinor corrections from #7513 are retained.

Scope:

• The shared LCAO CPU/GPU DeltaSpin atomic-force implementation is corrected.
• nspin=2 retains its existing multiplication.
• Stress already explicitly differentiates both overlaps; its normalization is unchanged.
• This force-only patch does not modify fixed-geometry SCF, energies, moments or magnetic-force definitions.
• PW has a separate force implementation with an explicit factor two. No PW or DFT+U code is changed or newly certified here.

Governance Notes

• INPUT/docs changes: No documentation update required: no parameter/default/parser/interface changes; no update to docs/parameters.yaml or input-main.md is required. Existing parameters are used for a bounded-cost force regression; the current case uses 12 Ry and at most 10 SCF steps without requiring convergence. Case references changed because the geometry, SOC, cutoff and iteration settings changed; the before/after numerical evidence is above.
• Core module impact: Only the final normalization of the shared LCAO DeltaSpin Operator force changes. No ESolver/HSolver/ElecState/Hamilt/Psi interface, global dependency, header or MPI-ownership changes. nspin=2 and stress retain their existing formulas.
• Exceptions requested: None. The previously recorded whole-file quality score was 53 at the base and 54 after the fix; below-60 historical complexity/style debt is not expanded in this focused correction.

[mohanchen]

maybe you can set up a smaller value to accelerate the calculations？