feat(xc): add gga_grad 1/2 with variationally consistent mode 2

[MagTheoryLab]

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 #7730. This PR supersedes the modes 1/2 implementation in #7758, with the original author's agreement to withdraw that PR and have this work submitted here. It builds on the original proposal by @dyzheng. This PR carries forward the gga_grad interface, the local-magnitude option, the projected-gradient construction, and the associated builtin/LibXC integration from that work. These are prior contributions, not features independently originated by this PR's submitter.

This submission covers modes 1 and 2 only. Mode 3 from the earlier proposal is outside this PR's scope; the handoff does not imply that its functionality is included here. In particular, it does not claim to implement or disprove the complete Scalmani–Frisch functional. The additional work here is the C2 radial map, the exact reverse of the specified discrete energy, the density/sigma sanitizer derivatives, consistent XC stress/bookkeeping, and the accompanying verification.

Dependency attribution: the current branch also includes TaoXia's genelpa and Pauli/y corrections from #7914, with original authorship and cherry-pick trailers preserved, plus Chen Chengbing's native ELPA/cuSolver upper-triangle fixes. These dependency changes must be separated from the GGA contribution during review. Once their upstream equivalents merge, the feature branch should be rebased to avoid submitting them twice. No separate issue is needed for the implementation and verification of the modes already proposed in #7758.

What's changed?

For magnetic nspin=4 calculations, add an explicit choice of local spin mapping:

• gga_grad=0 remains the default and selects the existing algorithm.
• gga_grad=1 carries forward Feature: add gga_grad parameter for noncollinear spin GGA functional #7758's local-magnitude option for the builtin GGA path, independent of the initial global quantization axis. For LibXC, modes 0 and 1 retain the same path.
• gga_grad=2 evaluates LDA/GGA using a C2-regularized local spin map and differentiates the same discrete energy used in the forward evaluation. The potential, gradient stress, and vtxc use that common definition.
• gga_grad=3 is rejected. Other spin configurations retain their existing behavior.

Problem in #7758 addressed by this PR

The reviewed #7758 snapshot is e641dba50e4ea86a03d126191dd271695fda9276. Its mode-2 implementation constructs projected spin gradients from the FFT gradients of the four density channels, then computes the magnetic divergence contribution by applying grad_dot to the spin-difference flux and multiplying the result by the local magnetization direction. See the builtin implementation at that snapshot and the commit describing the projected/full divergence distinction.

The resulting mode-2 potential is generally not the exact derivative of its own energy on a fixed FFT grid. In the smooth, unsaturated unregularized case, write $p_\mu=m_\mu/|m|$, $Q_{\nu\mu}=\partial p_\nu/\partial m_\mu$, $A=(a_+-a_-)/2$, and $H=(h_+-h_-)/2$, with the energy derivatives $a_s,h_s$ defined below. The distinction is

$$ v_\mu^{\mathrm{7758,proj}}=p_\mu A-p_\mu D_hH, $$

$$ v_\mu^{\mathrm{exact}}=p_\mu A +\sum_\nu Q_{\nu\mu},H\cdot G_hm_\nu-D_h(p_\mu H). $$

Equivalently, the missing fixed-grid correction is

$$ R_\mu=v_\mu^{\mathrm{exact}}-v_\mu^{\mathrm{7758,proj}} =\sum_\nu Q_{\nu\mu},H\cdot G_hm_\nu-D_h(p_\mu H)+p_\mu D_hH. $$

$R_\mu$ cancels in the continuum by the chain and product rules, but does not vanish identically for sampled fields and a finite FFT derivative. The correction here is to reverse the actual discrete energy graph, rather than assume that cancellation. Sections 1 and the five-point counterexample below provide the derivation and a reproducible fixed-grid failure of the shortcut.

The C2 map is an additional, explicitly documented change to the functional near zero magnetization. It is not presented as an algebraically equivalent rewrite of #7758. Its first and second derivatives, and the active density/sigma branches, must then be propagated consistently as described below. Likewise, this PR does not carry forward the assertion that retaining a full divergence alone makes mode 3 more accurate; it implements and validates mode 2 only.

The numerical evidence is separated accordingly: the FFT example establishes the discrete formula mismatch, the ABACUS real-PW tests validate the new implementation, and the Fe2 table compares SCF convergence against develop and common-fix controls. The Fe2 results are not a direct before/after run of #7758, and do not independently prove that the formula mismatch caused its SCF behavior. The default develop algorithm remains available as mode 0.

1. Differentiate the energy that is actually evaluated

Let $x=(n+n_{\mathrm{core}},m_x,m_y,m_z)$, let $N_s(x)$ be the local spin density for $s\in{+,-}$, and define

$$ J_{sA}=\frac{\partial N_s}{\partial x_A},\qquad K_{sAB}=\frac{\partial J_{sA}}{\partial x_B}. $$

For the discrete FFT gradient $G_h$, mode 2 uses the projected gradient and grid energy

$$ g_s=\sum_A J_{sA}G_hx_A,\qquad E_h=w\sum_i f\bigl(N_+(x_i),N_-(x_i),g_+(i),g_-(i)\bigr), \quad w=\Omega/N_{\mathrm{grid}}. $$

Write $a_s=\partial f/\partial N_s$ and $h_s=\partial f/\partial g_s$. On a fixed differentiable branch,

$$ \delta g_s=\sum_{A,B}K_{sAB},\delta x_B,G_hx_A +\sum_A J_{sA}G_h\delta x_A. $$

Using the grid inner-product adjoint $G_h^\dagger=-D_h$, substitution into $\delta E_h=w\sum_{i,B}v_B(i)\delta x_B(i)$ gives

$$ \boxed{ v_B=\sum_s J_{sB}a_s +\sum_{s,A}K_{sAB},h_s!\cdot G_hx_A -D_h!\left(\sum_sJ_{sB}h_s\right). } $$

The middle term is the response of the local-map Jacobian itself. Dropping it, or moving a spatially varying Jacobian outside the discrete divergence without the corresponding discrete identity, generally changes the derivative while leaving the forward energy unchanged. Also, $G_hN_s(x)$ and $\sum_AJ_{sA}G_hx_A$ need not be equal on a finite grid; replacing one by the other changes the discrete functional.

To see why the continuum formula alone is insufficient, consider the unsaturated positive-density map $N_\pm=(n\pm|m|)/2$. Put $p_\mu=m_\mu/|m|$, $A=(a_+-a_-)/2$, and $H=(h_+-h_-)/2$. In the continuum,

$$ v_\mu=p_\mu A+H\cdot\nabla p_\mu-\nabla\cdot(p_\mu H) =p_\mu(A-\nabla\cdot H). $$

Thus the projected expression is consistent with this continuum energy after the two cross terms cancel. Keeping only the full divergence is not, by itself, a more complete derivative of the same energy. On a finite FFT grid, however, the continuum product/chain rules cannot simply be assumed. The implementation therefore evaluates the boxed discrete reverse, including the Hessian contribution. This is the precise scope of the formula correction relative to the projected/full-divergence distinction described in #7758.

Code: xc_functional_ncgga_sf.cpp implements the builtin reverse; reverse_ncl_sf_discrete in libxc_tools.cpp implements its LibXC counterpart. Real-PW tests exercise the actual gradient/divergence operators rather than a pointwise FFT mock.

A small, reproducible discrete counterexample

This is a dimensionless toy functional illustrating the algebra, not a PBE/ABACUS benchmark or a measured failure of the #7758 executable. On five periodic points $t_j=2\pi j/5$, take the real spectral derivative

$$ (G_h)_{ij}=-\frac15\sum_{k=-2}^{2}k\sin[k(t_i-t_j)], $$

$$ m_j=(1+0.2\cos t_j)(\cos(0.7\sin t_j),\sin(0.7\sin t_j),0), \qquad \delta m_j=(\cos2t_j,0.3\sin t_j+0.2,0). $$

Let $q=p\cdot G_hm$, $p=m/|m|$, and $E_h=\tfrac12\sum_jq_j^2$. This is obtained from the projected-spin-gradient energy with $f=\tfrac12|g_+-g_-|^2$ on an unsaturated branch. Here $H=q$ and $A=0$. The exact discrete directional derivative is 0.006394957429, whereas using only $-pG_hq$ gives 0.149387398889.

Perturbation step	Central energy FD	Error of exact discrete reverse	Error of continuum-projected shortcut
1e-3	0.006393862364	1.0951e-6	0.1429935
5e-4	0.006394683663	2.7377e-7	0.1429927
2.5e-4	0.006394888988	6.8441e-8	0.1429925
1.25e-4	0.006394940318	1.7111e-8	0.1429925

The exact reverse approaches the energy difference at second order. The shortcut's error remains finite as the perturbation decreases at fixed grid size. This establishes the fixed-grid derivative mismatch; it does not assert failure of a continuum limit under grid refinement.

Reproduce with Python 3 (standard library only)

"""Five-point periodic spectral counterexample; Python standard library only.
Toy projected-gradient energy, not PBE or an ABACUS before/after benchmark.
"""
import math,json
N=5;x=[2*math.pi*j/N for j in range(N)]
D=[[-sum(k*math.sin(k*(xi-xj)) for k in range(-2,3))/N for xj in x] for xi in x]
def grad(v):return [sum(D[i][j]*v[j] for j in range(N)) for i in range(N)]
m=[[(1+.2*math.cos(t))*math.cos(.7*math.sin(t)),(1+.2*math.cos(t))*math.sin(.7*math.sin(t))] for t in x]
u=[[math.cos(2*t),.3*math.sin(t)+.2] for t in x]
def evaluate(m):
 r=[math.hypot(*a) for a in m];p=[[a/r[i] for a in m[i]] for i in range(N)];gm=[grad([a[b] for a in m]) for b in range(2)]
 q=[sum(p[i][b]*gm[b][i] for b in range(2)) for i in range(N)]
 E=.5*sum(a*a for a in q);dq=grad(q)
 projected=[[-p[i][b]*dq[i] for b in range(2)] for i in range(N)]
 flux=[grad([p[i][b]*q[i] for i in range(N)]) for b in range(2)]
 exact=[[-flux[b][i]+sum(((1. if a==b else 0.)-p[i][a]*p[i][b])/r[i]*q[i]*gm[a][i] for a in range(2)) for b in range(2)] for i in range(N)]
 return E,sum(exact[i][b]*u[i][b] for i in range(N) for b in range(2)),sum(projected[i][b]*u[i][b] for i in range(N) for b in range(2))
E,a,b=evaluate(m);rows=[]
for eps in [1e-3,5e-4,2.5e-4,1.25e-4]:
 plus=[[m[i][j]+eps*u[i][j] for j in range(2)] for i in range(N)];minus=[[m[i][j]-eps*u[i][j] for j in range(2)] for i in range(N)]
 fd=(evaluate(plus)[0]-evaluate(minus)[0])/(2*eps)
 rows.append(dict(eps=eps,fd=fd,exact_error=abs(fd-a),projected_error=abs(fd-b)))
result=dict(energy=E,exact_directional_derivative=a,continuum_projected_directional_derivative=b,steps=rows)
print(json.dumps(result,indent=2))
assert rows[-1]['exact_error']<1e-7 and rows[-1]['projected_error']>1e-3
assert all(rows[i+1]['exact_error']<.3*rows[i]['exact_error'] for i in range(3))

The script was run successfully; the decreasing-step and residual-error assertions passed.

2. C2 regularization must propagate through both derivatives

With $r=|m|$, $t=r/\eta$, and fixed $\eta=10^{-3}$ in atomic density units, define

$$ S_\eta(r)=\begin{cases} \eta(3t^5-8t^4+6t^3),&0\le r<\eta,\\ r,&r\ge\eta. \end{cases} $$

Inside the regularized interval,

$$ S_\eta'=15t^4-32t^3+18t^2,\qquad S_\eta''=(60t^3-96t^2+36t)/\eta. $$

At $r=\eta$, the value/first/second derivatives are $(\eta,1,0)$, matching the outer branch. As $r\to0$, both $S_\eta'/r$ and $S_\eta''$ vanish. Consequently the Cartesian gradient and Hessian have finite, continuous limits:

$$ p_\mu=S_\eta'\widehat m_\mu,\qquad Q_{\mu\nu}=\frac{S_\eta'}r\delta_{\mu\nu} +\left(S_\eta''-\frac{S_\eta'}r\right)\widehat m_\mu\widehat m_\nu. $$

The code uses this $Q$ in the middle term of the discrete reverse. Replacing only the forward $r$ by $S_\eta(r)$ while retaining an unregularized magnetic derivative would define a different potential.

The actual spin map additionally uses

$$ q=|n+n_{\mathrm{core}}|,\qquad c=\min(S_\eta(r),q),\qquad N_\pm=(q\pm c)/2. $$

Its Jacobian follows the selected abs/saturation branch. C2 regularity of the radial function does not make abs/min or all LibXC clipping boundaries globally smooth; tests distinguish branch interiors and radial joins, and the implementation specifies branch choices at the boundaries.

3. LibXC energy weights and clipping belong to the derivative

For one active density branch, let $y_s=\max(T,\rho_s)$, $Y=\sum_sy_s$, and $M=\sum_sb_s\rho_s$, where $b_s$ is the fixed mask on that branch. ABACUS accumulates $M\epsilon$, whereas LibXC's vrho differentiates $Y\epsilon$. Therefore

$$ \frac{\partial(M\epsilon)}{\partial\rho_s} =b_s\epsilon+\frac MY,\mathbf1_{\rho_s>T}(v_s^{\mathrm{LibXC}}-\epsilon). $$

For sanitized sigma variables $\widetilde\sigma=C(\sigma)$,

$$ \frac{\partial(M\epsilon)}{\partial\sigma} =\frac MY(DC)^T v_\sigma^{\mathrm{LibXC}}. $$

For example, when the cross component is clipped at $\pm(\widetilde\sigma_{++}+\widetilde\sigma_{--})/2$, its derivative transfers to the diagonal components with weights $\pm1/2$; it is not sufficient to zero only the cross derivative. make_libxc_weighted_derivatives accounts for these density and sigma branches before the reverse above. The inactive total-density branch has zero energy and derivative under the implemented threshold convention.

Finally, bookkeeping uses the completed returned four-component potential:

$$ v_{\mathrm{txc}}=w\sum_i\left[v_n(i)n(i)+\mathbf v_m(i)\cdot\mathbf m(i)\right]. $$

The frozen core participates in XC evaluation but is not added to the valence-density contraction in this ledger.

4. Stress uses the same gradient graph

For a homogeneous infinitesimal strain $\varepsilon$ with grid field values held fixed, $\delta g_{s,i}=-\varepsilon_{ji}g_{s,j}$. Hence

$$ \delta E_h^{\mathrm{metric}} =-w\sum_{k,s,i,j}h_{s,i}(k)g_{s,j}(k),\delta\varepsilon_{ji}. $$

This gives the gradient-metric contribution used by the stress assembly, with the same map and sanitized derivatives as the energy. Volume, density-rescaling, local-XC and core terms must be included separately for the full stress. The tests below distinguish metric-only, frozen-field full-XC, and total self-consistent stress; those are different verification scopes.

Formula contracts of the imported dependencies (credited to their own changes)

For Hermitian $H,S$ with authoritative upper triangles, the lower entries must obey $A_{ji}=\overline{A_{ij}}$. The generalized eigenproblem $HC=SC\Lambda$ must not depend on arbitrary unused lower-buffer contents. Native ELPA reconstructs private Hermitian inputs and cuSolver selects the upper triangle. This belongs to the separate solver fix, not to the GGA derivation.

For the standard Pauli convention,

$$
\lambda\cdot\sigma=
\begin{pmatrix}
\lambda_z&\lambda_x-i\lambda_y\
\lambda_x+i\lambda_y&-\lambda_z
\end{pmatrix},\qquad
m_y=\operatorname{Tr}(\rho\sigma_y)=-2\operatorname{Im}\rho_{\uparrow\downarrow}.
$$

The operator and moment must use the same y convention. TaoXia's #7914 correction supplies this mapping; it is not a new correction claimed by this GGA PR.

Unit Tests and/or Case Tests for my changes

Derivative evidence: the real-PW suite passed 32/32 in serial and 32/32 with two MPI ranks. It covers the discrete adjoint, density-channel derivatives, radial/saturation branches, core response and translation, local rotation torque, XC bookkeeping, and stress. The density-direction checks include a decreasing-step error test, not just one favorable finite difference.

Frozen-field check	Maximum of the best scaled errors over tested directions
Density-direction derivative, 56 grouped directions	2.3252e-8
Gradient-metric stress, 54 grouped components	2.5108e-11
Full-XC diagonal stress, 15 grouped components	3.2159e-8

Here the scaled error is $|a-d|/\max(1,|a|,|d|)$ in the test's atomic-unit convention, and the best value is selected within each tested step sequence. The density-direction threshold is 3e-8 with additional step-decrease gates. These are current implementation results; this table is not a measured failing-then-passing A/B against the head of #7758. Such a #7758-head comparison has not been run in this validation campaign.

Self-consistent Fe2 comparison: use the same two-atom bcc Fe LCAO/SOC case, 100 Ry cutoff, kspacing=0.4, Gaussian smearing 0.01 Ry, Broyden mixing beta 0.2/ndim 30, scf_thr=scf_ene_thr=1e-6, and a 200-step limit. CPU/GPU use the same CUDA-enabled binary with device switched, native elpa, and OMP_NUM_THREADS=1. Environment: Tesla V100-SXM2-32GB, CUDA 12.9.1, OpenMPI 5.0.10, ELPA 2026.02.001, LibXC 7.0.0.

Source / mode	CPU iterations; final DRHO	GPU iterations; final DRHO
Official develop fb9ce1d	200, not converged; 1.1802e-5	200, not converged; 1.5394e-4
Develop + common solver/Pauli fixes	200, not converged; 4.5907e-5	200, not converged; 1.1845e-5
This branch, mode 0	200, not converged; 4.5907e-5	200, not converged; 9.5818e-6
This branch, mode 1	21, converged; 9.5448e-7	21, converged; 9.5448e-7
This branch, mode 2	17, converged; 7.8127e-7	17, converged; 7.8123e-7

The common-fix control separates GGA effects from eigensolver dependencies. The CPU mode-0 printed trajectory and final energy/forces/stress match the common-fix control exactly; the GPU mode-0 unconverged trajectories are not identical and are not claimed as a converged parity result. Convergence on one case does not prove a general convergence guarantee or determine the correctness of a different functional.

Converged comparison	Absolute energy difference (eV)	Maximum stress difference (kbar)
Mode 1 CPU/GPU	9.094947e-13	2.799993e-9
Mode 2 CPU/GPU	1.151420e-9	2.926000e-7
Mode 2 one/four CPU MPI ranks	1.237822e-9	2.799510e-5

All meet the declared 1e-5 eV / 1e-4 eV/Angstrom / 0.01 kbar energy/force/stress tolerances. Forces are zero at the reported precision in this symmetric case; their agreement is not a nonzero-force derivative test. The four-rank mode-2 run also converged in 17 steps.

Commands run: $CPU_BUILD and $GPU_BUILD denote the configured CPU and CUDA build directories; $SOURCE_DIR is the verified source snapshot. Expanded build/job scripts and input/output evidence are retained in the validation record.

export OMP_NUM_THREADS=1
cmake --build "$CPU_BUILD" -j8 --target abacus_basic_para MODULE_HAMILT_XCTest_VXC MODULE_HAMILT_XCTest_NCGGA_RADIAL MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD MODULE_HAMILT_XCTest_GRADCORR MODULE_IO_input_test_para deltaspin_core_test
ctest --test-dir "$CPU_BUILD" --output-on-failure -R 'MODULE_HAMILT_XCTest_(VXC|GRADCORR|NCGGA)|MODULE_IO_input_test_para$|deltaspin_core_test$'
cmake --build "$GPU_BUILD" -j8 --target abacus_basic_gpu MODULE_HAMILT_XCTest_NCGGA_RADIAL MODULE_HAMILT_XCTest_NCGGA_DISCRETE_FD MODULE_HSOLVER_LCAO MODULE_HSOLVER_LCAO_cusolver deltaspin_core_test
ctest --test-dir "$GPU_BUILD" --output-on-failure -R 'MODULE_HAMILT_XCTest_NCGGA|deltaspin_core_test$'
"$CPU_BUILD/abacus_basic_para" --version
"$CPU_BUILD/abacus_basic_para" -h gga_grad
"$CPU_BUILD/abacus_basic_para" --generate-parameters-yaml
# In each prepared valid case directory, with gga_grad=0/1/2/3:
"$CPU_BUILD/abacus_basic_para" --check-input
# In each prepared Fe2 directory, with the selected immutable executable:
mpirun -np 1 "$ABACUS"
mpirun -np 4 "$ABACUS"  # mode-2 CPU MPI comparison

git diff --check upstream/develop HEAD
python3 tools/03_code_analysis/agent_governance_check.py --base upstream/develop --head HEAD --format text

CPU CTest 7/7 passed; GPU-build CTest 4/4 passed. The CUDA-build real-PW derivative test still exercises host-side FFT/density code; actual GPU SCF coverage comes from the case comparison. CLI accepted 0/1/2 and rejected 3; version v3.11.0-beta9; the generated gga_grad documentation matches the tracked metadata. A separate ENABLE_LIBXC=OFF configuration compiled target xc_; a complete no-LibXC executable was not built in that check. The final native ELPA regressions passed 3/3 on one and four ranks, and the cuSolver/ScaLAPACK group passed 6/6.

Limitations and checks not accepted:

• Full-SCF force/stress finite differences are not validated. After displacing the second atom by +0.03 in fractional x and tightening both SCF thresholds to 1e-9, the central structures exhausted 250 steps: DRHO=2.053e-5 (mode 1) and 1.0015e-5 (mode 2). Dependent perturbation runs were stopped/skipped rather than interpreting unconverged energies as derivatives.
• The supplied case has no atom-wise spin-constraint mask. Its SCF results do not establish a nonzero y-constraint response; the imported Pauli change has focused unit coverage.
• The installed ELPA two-stage CPU path fails in multi-rank tests independently of ABACUS. That library was not changed; the SCF comparison uses native one-stage elpa. No blanket genelpa MPI pass is claimed.
• Full repository, ROCm, broad hybrid/meta-GGA, and OMP>1 validation were not performed. No mode-3 data or historical FeGe results are used as evidence for this modes-1/2 submission.

Governance Notes

• INPUT/docs: new gga_grad, default 0, accepted 0/1/2; fixed mode-2 eta=1e-3, not a new user parameter. C++ registration and generated docs/parameters.yaml / docs/advanced/input_files/input-main.md are updated.
• Core-module impact: XC energy, four-component potential, gradient stress and core-response call sites; ElecState potential assembly explicitly forwards the mode. HSolver and DeltaSpin/Operator changes belong to the credited dependency fixes. No new ESolver or Psi workflow switch is introduced. New XC sources/tests are wired through CMake.
• Global/header review: global references decrease by 28 in the total diff; remaining new references forward configuration at existing upper-layer boundaries. Complete matrix/array/vector types are required by the affected public declarations. Governance reports warnings but no errors.
• Exceptions requested: none.

[mohanchen]

I will handle this PR later.

[mohanchen]

I'd recommend not modifying the code this way. It would be better to pass the required variables directly as function input parameters, rather than using PARAMETERS/PARAM in any form.

[mohanchen]

I'd recommend not modifying the code this way. It would be better to pass the required variables directly as function input parameters, rather than using PARAMETERS/PARAM in any form.

[mohanchen]

I'd recommend not modifying the code this way. It would be better to pass the required variables directly as function input parameters, rather than using PARAMETERS/PARAM in any form.

[mohanchen]

I would suggest not using the number 0 directly, because other developers will not know what it stands for, I would suggest using a local variable (for example xxx=0 and then use xxx) here.