Use mapped-density FFT gradients for noncollinear gga_grad=2

[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

Related to #8128. Fermi-level convergence is handled separately in #8130; this PR contains no occupation/Fermi-solver changes and does not depend on that PR.

What's changed?

For magnetic nspin=4, gga_grad=2, evaluate FFT gradients after mapping the four density channels to two local spin densities, then differentiate that same discrete energy. Apply this definition consistently to built-in XC, LibXC, the returned four-component potential/vtxc, and the GGA metric-stress contribution.

This changes the finite-grid local-collinear functional. The upstream projected-gradient implementation already differentiates its own discrete energy; this PR is not a claim that its Hessian/divergence terms are simply missing. At finite FFT cutoff, the two forward energy definitions differ. In particular, the projected definition can respond to pointwise spin-direction changes that preserve the local mapped densities, and its branch-dependent gradient can jump at the density-saturation boundary. The mapped-density definition preserves the local-collinear invariance and continuity properties tested below.

gga_grad=0/1, default selection (0), radial regularization, clipping thresholds, the LDA local map, occupations, DeltaSpin, force/stress sign fixes and solvers retain the upstream implementation. No additional gga_grad mode, input switch, CUDA XC kernel, or local-branch GPU infrastructure is introduced. CUDA-enabled ABACUS uses the upstream XC dispatch; its SCF regression must not be confused with testing a new device-native XC kernel.

Discrete energy and potential

Let $x=(n+\rho_{\rm core},m_x,m_y,m_z)$, $N_s(x)$ be the existing regularized/clipped spin-density map, and $J_{sA}=\partial N_s/\partial x_A$. Let $G_h$ be the retained plane-wave FFT gradient, $D_h=-G_h^*$ its divergence, and $w=\Omega/N_{\rm grid}$.

The existing projected graph uses $g_s^{\rm proj}=\sum_A J_{sA}G_hx_A$. The proposed graph is:

$$E_{\rm XC}^{\rm map}=w\sum_i f\!\left(N_+(x_i),N_-(x_i),(G_hN_+)_i,(G_hN_-)_i\right).$$

With $h_s=\partial f/\partial g_s$, reversing FFT first, local map second gives:

$$v_A(i)=\sum_{s=\pm}J_{sA}(x_i)\left[\frac{\partial f}{\partial N_s}(i)-(D_hh_s)_i\right].$$

No pointwise discrete product/chain rule is assumed: in general $G_hN_s(x)\ne\sum_AJ_{sA}(x)G_hx_A$. The projected graph's Hessian term disappears because that graph has been replaced, not because a derivative term is being dropped. The built-in metric stress and LibXC gradient data now use the same $G_hN_s$ as the energy. Density/sigma branch masks and the existing radial map are retained; clipping kinks are not claimed to become differentiable.

Code anchors at the tested commit: built-in forward gradient, built-in reverse, LibXC reverse, and new regression tests.

Unit Tests and/or Case Tests for my changes

Fresh verification uses upstream 870ce1589445c41726af9f655c78df2f3cfa60c1 and candidate d18f769ddfd6a300c5717382e1faee3e6d7a48fa, built with the same toolchain. All three fresh SCF runs below passed their joint charge/energy/moment gates.

Fresh negative control: the upstream production XC files pass all 40 pre-existing derivative/covariance tests, but fail 14 of the 16 new property tests (six local-rotation cases and all eight clipping-boundary cases). The candidate passes 56/56. For the smooth-density fixture, rotating only the local spin direction changes the upstream XC energy by 0.0102146466 Ry (both built-in and LibXC); the mapped implementation satisfies the 1e-10 Ry energy and 1e-11 Ry component-potential tolerances. These are fixture-level quantities, not per-atom material errors. At the clipping boundary, the upstream built-in smooth fixture has energy jumps of 0.02125655 and 0.02385027 Ry as the relative offset halves from 5e-6 to 2.5e-6; the candidate passes the shrinking-jump criterion at all four offsets.

Build/test environment: NVIDIA V100, one node, GCC with CUDA 12/NVHPC 26.3, OpenMPI 5.0.10, LibXC 7.0.0, FFTW 3.3.10 and ELPA 2026.02.001. OMP_NUM_THREADS=1; CTest exercises one/two ranks; each SCF uses one MPI rank. The CPU case uses device cpu in the same CUDA-enabled executable.

Commands run (src/build denote the checkout and build directory; library roots are provided by the above modules):

export OMP_NUM_THREADS=1
cmake -S "$src" -B "$build" \
  -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_COMPILER=g++ -DMPI_CXX_COMPILER=mpicxx \
  -DUSE_CUDA=ON -DCMAKE_CUDA_ARCHITECTURES=70 -DENABLE_LCAO=ON \
  -DENABLE_MPI=ON -DENABLE_OPENMP=ON -DENABLE_ELPA=ON -DENABLE_LIBXC=ON \
  -DENABLE_LIBRI=OFF -DENABLE_JSON=OFF -DENABLE_CUSOLVERMP=OFF \
  -DBUILD_TESTING=ON -DCOMMIT_INFO=OFF -DGIT_SUBMODULE=OFF \
  -DFETCHCONTENT_SOURCE_DIR_GOOGLETEST="$gtest_source" \
  -DLAPACK_DIR="$OPENBLAS_ROOT/lib" -DSCALAPACK_DIR="$ScaLAPACK_ROOT/lib" \
  -DELPA_DIR="$ELPA_ROOT" -DFFTW3_DIR="$FFTW3_ROOT" -DLibxc_DIR="$LIBXC_ROOT"
cmake --build "$build" -j8 --target MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD \
  MODULE_HAMILT_XCTest_NCGGA_RADIAL MODULE_HAMILT_XCTest_VXC abacus_basic_gpu
ctest --test-dir "$build" --output-on-failure \
  -R '^MODULE_HAMILT_XCTest_(NCGGA_DISCRETE_FD|NCGGA_DISCRETE_FD_MPI2|NCGGA_RADIAL|VXC)$'
"$build/source/source_hamilt/module_xc/test/MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD" \
  --gtest_output=xml:candidate-tests.xml

Build passed; 4/4 CTest targets passed, including the 56-case production fixture on both rank counts. The first attempt at MPI2 failed before test execution because the batch allocation supplied one slot; rerunning with two allocated slots passed, with no source/tolerance change.

The production-PW fixture covers built-in PBE and LibXC PBE, smooth/negative-density/saturated/radially regularized branches, all four density derivatives plus the core-density derivative, the vtxc contraction, all six metric-stress components, and gradient/divergence adjointness, on one and two MPI ranks. Added regressions check pointwise spin rotation and energy continuity when crossing saturation. Existing derivative tolerances are retained.

The negative control rebuilds the three XC production files from the exact upstream base against the same new tests, then restores/rebuilds the candidate. FeRh SCF uses identical structure, constraints, pseudopotentials, orbitals and initialization; ecutwfc=100 Ry, ecutrho=400 Ry, 4×4×4 k points, noncollinear PBE, SOC off. No mixing retuning or Fermi patch is used in this fresh comparison.

CLI and generated documentation checks passed; SCF runtime status is listed separately below:

"$build/abacus_basic_gpu" --version
"$build/abacus_basic_gpu" -h gga_grad
"$build/abacus_basic_gpu" --generate-parameters-yaml > docs/parameters.yaml
python3 docs/generate_input_main.py docs/parameters.yaml --output docs/advanced/input_files/input-main.md
# In each of the three FeRh case directories:
"$build/abacus_basic_gpu" --check-input
mpirun -np 1 --map-by "$MAP_OPT" "$build/abacus_basic_gpu"
git diff --cached --check
python3 tools/03_code_analysis/agent_governance_check.py --staged --format text

The executable reports v3.11.0-beta10. Generated documentation changes only the gga_grad=2 paragraph; governance reports no findings. Source SHA-256 checks before testing and after baseline restoration match this patch. The SCF launch uses the site's one-rank GPU mapping; its executable is selected separately for the baseline and candidate runs.

Fresh same-base FeRh results

Two-atom distorted B2-derived FeRh, using the ten-structures/05_FeRh_B2 case in the public reproduction archive. The lattice constant is 3.06 Å and the second fractional position is (0.517, 0.491, 0.506). Fe/Rh vector targets are (0.888594398178, 0.513030214989, 2.81907786236) and (0.0855050358314, 0.484923155196, −0.0868240888335) μB.

Implementation / device	SCF steps	Final energy (eV/cell)	Final charge residual	Final magnetic RMS (μB)	Wall time (s)
develop / GPU	23	−6601.263315878381	9.1101e-8	7.54827e-7	126.39
This PR / GPU	23	−6601.263256223581	7.3484e-8	4.29083e-7	128.11
This PR / CPU	23	−6601.263256224059	7.3484e-8	4.29079e-7	688.17

All three runs pass charge residual ≤ 1e-7, absolute last-step energy change ≤ 1e-7 eV, and magnetic RMS ≤ 1e-6 μB, with the latter independently recomputed from final moment vectors. The XC-only change raises the converged GPU energy by 0.0298274 meV/atom. Candidate CPU/GPU final energies differ by only 4.78e-10 eV/cell, and both converge in 23 steps. This case converges with both definitions; it does not establish a universal convergence improvement. Timings are single-run observations, not a statistically established speedup/regression.

Existing weak-moment SCF regression: reference and solver settings

The initial PW integration run failed only 038_PW_NC_GGA_GRAD2. Energy passed; its force/stress fingerprints differed from the old reference by 0.000624 eV/Å / 0.007202 kbar. Rebuilding the exact upstream XC files with the same toolchain reproduced the old reference, while candidate CPU/GPU runs reproduced the CI values. This Fe2 fixture starts with nonparallel moments that decay towards zero: it tests the weak-moment limit, not a stable magnetic ground state.

At the original scf_thr=1e-10, the force/stress output depends on the finite-accuracy SCF terminal state even when energies agree within 1e-9 eV. The fixture now uses the existing BPCG solver with scf_thr=1e-14 and a regenerated CPU one-rank reference. Tight-threshold Davidson probes hit psi_norm <= 0 on CPU and GPU, so tightening Davidson alone was not retained. No comparison tolerances, eigensolver implementation, production defaults, or Fermi code are changed.

XC graph / solver / SCF threshold / device	Force fingerprint (eV/Å)	Stress fingerprint (kbar)
Upstream / Davidson / 1e-10 / CPU1	18.267150	73803.586958
Mapped / Davidson / 1e-10 / CPU1	18.266528	73803.579752
Upstream / BPCG / 1e-14 / CPU1	18.270388	73803.529922
Mapped / BPCG / 1e-14 / CPU1 (new reference)	18.270388	73803.529922
Mapped / BPCG / 1e-14 / CPU4	18.270388	73803.529922
Mapped / BPCG / 1e-14 / GPU1	18.270408	73803.529885

These are sums of absolute components, not vector norms or derivative-error bounds. The upstream and mapped CPU1 energy/force/stress fingerprints agree at the reported precision under the final settings. CPU4 and GPU1 pass against the new CPU1 reference with the unchanged default tolerances (energy 1e-7 eV, force fingerprint 1e-4 eV/Å, stress fingerprint 1e-3 kbar). This weak-moment limit does not imply that the two XC definitions give identical energies for finite magnetization.

To check the physical meaning rather than merely update a reference, both implementations were also run at atom-2 x displacements ±0.002/±0.001 Å and xx strains ±1e-4/±5e-5. All 16 finite-difference corners converged. At the final BPCG settings, both implementations give force errors 2.80e-5 / 4.24e-5 eV/Å and xx stress errors 0.0168 / 0.0171 kbar. For comparison, the original Davidson/1e-10 mapped force errors were 3.55e-4 / 3.40e-4 eV/Å (upstream: 3.14e-4 / 2.99e-4). Residual stress errors remain in this 20-Ry fixture, equally for the two implementations; this is not a claim of complete end-to-end stress closure and does not replace the finite-magnetization XC module derivative tests.

Focused reproduction, from tests/01_PW, with this case listed in CASES_CPU.txt:

OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a /path/to/abacus -n 1 -o 1 -r '^038_PW_NC_GGA_GRAD2$' -g
OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a /path/to/abacus -n 4 -o 1 -r '^038_PW_NC_GGA_GRAD2$'

The candidate was additionally checked against the committed reference using CPU1, CPU4, and GPU1 (device gpu for the last), with explicit checks of #SCF IS CONVERGED# in each application log. This final check ran on a second V100 node with the same toolchain and one OpenMP thread. All three passed; the largest force/stress fingerprint deviations from the committed reference were 1.6e-5 eV/Å / 4.72e-4 kbar, below the unchanged tolerances. The first CI revision passed all build configurations and its CUDA test job; the PW integration job is rerun for this fixture update.

LCAO references for the changed finite-grid functional

The next CPU integration run passed 01_PW, 02_NAO_Gamma, and all preceding unit-test groups, then exposed two stale references in 03_NAO_multik. A repository-wide INPUT audit found only three explicit gga_grad=2 integration cases: the PW case above and these two LCAO cases. Only the two LCAO references and their README explanations are updated here; INPUT, iteration budgets, comparison tolerances, and production source are unchanged.

To distinguish the intended functional change from a new defect, the exact upstream XC files were rebuilt in the same checkout/toolchain, with the candidate restored and rebuilt afterwards. Upstream one-/four-rank CPU runs reproduce the previous references (energy within 1e-10 eV, identical force/stress fingerprints). Mapped one-/four-rank runs reproduce the failing CI numbers, then pass all 8 energy/force/stress comparisons against the regenerated one-rank references at the unchanged default tolerances.

Case / XC graph	State being compared	Energy (eV/cell)	Force fingerprint (eV/Å)	Stress fingerprint (kbar)
scf_gga_grad2_spin4 / upstream	Converged, 15 Ry	-6267.500675032697	0.000000	3941.148224
scf_gga_grad2_spin4 / mapped	Converged, 15 Ry	-6267.393141439436	0.000000	3873.105155
scf_deltaspin4 / upstream	Step 10, not converged, 12 Ry	-6807.615155408896	99.992332	30048.363075
scf_deltaspin4 / mapped	Step 10, not converged, 12 Ry	-6807.797913831613	93.374790	30504.021979

The converged 15-Ry case shifts by +53.7668 meV/atom. Thus the small 100-Ry FeRh shift above is not a universal energy-error bound: coarse FFT grids can produce substantially different finite-grid functionals. The 12-Ry DeltaSpin case intentionally stops after ten SCF steps; its -0.1827584227 eV/cell trajectory shift includes density/constraint response and must not be described as a converged energy difference or a force/stress finite-difference result. The absolute-component sums in the table are regression fingerprints, not component-wise accuracy bounds.

Focused same-toolchain commands, from tests/03_NAO_multik, with both cases listed in CASES_CPU.txt:

OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a /path/to/abacus -n 1 -o 1 -r '^(scf_gga_grad2_spin4|scf_deltaspin4)$' -g
OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a /path/to/abacus -n 4 -o 1 -r '^(scf_gga_grad2_spin4|scf_deltaspin4)$'

Both commands completed for each XC implementation. Application logs independently confirm convergence for scf_gga_grad2_spin4 and nonconvergence at the intentional iteration limit for scf_deltaspin4. This follow-up does not claim new GPU LCAO or end-to-end derivative validation; the previous CI CUDA job passed at revision 67b0e1914. The complete CPU Integration Test and Unit Test run passed at 0e6337bf4, including 03_NAO_multik and all downstream test groups that were skipped after the previous failure.

Earlier ten-structure convergence screen

The earlier screen in #8128 covers 10 fixed structures × 4 archived implementations (40 completed runs). Mapped-density XC passed 10/10, versus 7/10 for the tested official develop, 8/10 for tmp, and 6/10 for the historical accel build.

Green cells give iterations after all convergence gates pass; red cells reached the 200-step limit without joint convergence. The mapped candidate converged the three cases that failed for the tested official baseline: fcc Fe (28 steps), CoO (87), and NiO (31). It is not uniformly faster: FeO took 178 steps versus 93 for official develop; FeRh took 23 in both. Thus the observed benefit is broader convergence in this fixed sample, not universal acceleration or a material-independent success probability.

Exact iteration counts for all ten structures

Structure	tmp	Historical accel	Tested develop	Mapped candidate
Fe bcc	20	37	40	26
Fe fcc	28	Not converged (200)	Not converged (200)	28
Co hcp	18	34	36	21
Ni fcc	20	26	29	22
FeRh B2	24	27	23	23
FeCo B2	24	38	34	25
FeNi B2	27	29	30	27
FeO rocksalt	Not converged (200)	Not converged (200)	93	178
CoO rocksalt	Not converged (200)	Not converged (200)	Not converged (200)	87
NiO rocksalt	26	Not converged (200)	Not converged (200)	31
Accepted / 10	8	6	7	10

Matched inputs and gates. Within each four-version group, structures, vector-moment targets, pseudopotentials, projector orbitals and explicit INPUT settings were identical. Noncollinear constrained PBE; SOC off; 100 Ry wavefunction / 400 Ry density cutoff; 4×4×4 k points; Gaussian 0.01 Ry; Broyden beta=0.2, beta_mag=0.8; Davidson subspace; 1 MPI rank / 1 V100; OMP=1; 200-step budget. Acceptance required normal exit, explicit SCF convergence, DRHO ≤1e-7, absolute final energy change ≤1e-7 eV, and full-PW magnetic RMS ≤1e-6 μB, independently checked against the final atomic moments. Four initial wall-time truncations were rerun unchanged with an extended safety limit: none of the final 40 results is time-censored. The oxides use globally tilted AFM-I metal targets, unconstrained oxygen, and PBE without U; these are numerical test structures rather than experimental-ground-state claims.

Version boundary. These are archived whole-build results, not a rerun of this exact PR on the current base. The official baseline was f8118a0de8c8818c7c35f32ed443c18ad70a486d; accel was 4cc4a9769965c025045c07a1fb074d3e86de8d14; the historical tmp build has five XC files byte-matched to tmp@19ade18599bfdf9ee78a289bf798a9b2583d52cd. For explicit gga_grad=2 and built-in PBE, the audited accel/develop XC cores use the same projected-gradient definition; their different whole-build iteration counts are not evidence of a fourth independent XC method. The candidate contains the mapped-density XC patch and Pauli-y force correction and predates the separate Fermi-solver fix. Its build also has local GPU infrastructure absent from this upstream-only PR. Therefore 10/10 must not be presented as proof that this XC-only PR alone caused every improvement. The fresh same-base FeRh A/B above supplies the isolated control; this larger screen supplies supporting coverage without repeating 40 costly runs. Inputs, fingerprints and numerical summaries are in the public evidence archive.

Across the seven jointly converged official/candidate cases, candidate-minus-official energies range from −0.012533 to +1.695695 meV/atom (FeO and bcc Fe respectively; FeRh: +0.029827 meV/atom). These are whole-build self-consistent energy differences, potentially including different constrained local minima, not XC-only errors. No energy difference is assigned to an unconverged baseline. Forces/stress were disabled in this screen; passing magnetic RMS is not a magnetic-force derivative test.

Earlier evidence and numerical limits

The wider materials and derivative study is documented in #8128. Those measurements used archived whole builds, so they are supporting evidence rather than a same-base XC-only acceptance test for this PR. Earlier grid benchmarks held wavefunction cutoff at 100 Ry:

Single-grid refinement at density cutoff 400 Ry	Mapped implementation: SCF steps	Mapped energy shift (meV/atom)	Archived develop energy shift (meV/atom)
bcc Fe, 36³ → 60³	26 → 27	+0.014708	+0.558346
B2 FeRh, 40³ → 60³	23 → 22	+0.001299	+0.002997

The mapped benchmark build included the separate Fermi repair, so these rows are explicitly not labeled as measurements of this PR's exact executable. They show no need for a denser sampling grid to obtain SCF convergence in those two cases, not a general grid-independence guarantee.

At fixed 60³ density grid, increasing the density cutoff 400 → 600 → 800 Ry changed the mapped FeRh energy by 0, −0.012492, −0.016771 meV/atom; iterations were 49/50/49. The nonlinear map has measurable high-frequency sensitivity even though its discrete derivative is consistent. These are finite-range sensitivity data, not an infinite-cutoff reference. The double-grid path also showed shared bcc Fe instability in archived implementations; three diagnostic trajectories were stopped early and must not be presented as 200-step failures or blamed on this XC change.

Checks not run: no repeat of the full ten-material campaign, the FeRh atomic/magnetic-force finite-difference campaign, or tiny-angle J/D series on this exact branch; no ROCm validation. The fresh short SCF cases disable force/stress output. The fresh six-component unit check validates the XC metric contribution, not a complete end-to-end strained-cell stress finite difference. The full DFT functional is not invariant under arbitrary local spin rotations; the new invariant test applies specifically to this local-collinear XC energy.

Governance Notes

• INPUT/docs changes: update the gga_grad=2 description in the C++ Input_Item and header comment; regenerate the corresponding YAML/Markdown documentation. No default, range, new parameter, or parsing change.
• Core module impact: replace the forward/reverse graph within XC and its GGA metric stress; remove unused gradient storage from NclSfDiscreteData. No global-dependency increase, new interface default, source file, or CMake target.
• Exceptions requested: none. Whole-file quality scores include historical debt: libxc_abacus.h 41 → 41 and libxc_tools.cpp 49 → 58; built-in XC 83 → 93 and the production-PW test remains 92. Existing low-scoring unrelated interfaces are not rewritten in this numerical change.