fxc contraction: GKS (noncollinear) for LDA

[susilehtola]

Adds the GKS overload

eval_fxc_contraction( Ps, Pz, Py, Px, tPs, tPz, tPy, tPx )
  -> ( FXCs, FXCz, FXCy, FXCx )

in the argument order of the GKS eval_exc_vxc, on the host replicated
integrator. Other integrators throw through defaulted virtuals, so no
device stubs are needed.

Method

The energy is the locally collinear one the GKS potential already uses,
n_+- = (rho_s +- |m|)/2 fed to the spin-polarized functional. The kernel
applied to a trial density, per noncollinear field slot, is generated
(host/gauxc_nc_kernel.hpp, xckernel's ncwriter: the mechanical second
derivative of that map) and assembled exactly as the potential is, so
FXC_X = d/dh VXC_X(P + h tP) for X = s, z, y, x.

Below gks_dtol the map is not twice differentiable (the transverse
kernel carries 1/|m|). There the generated collinear limit is used: each
magnetization component responds like the spin channel of a collinear
perturbation about the spin-symmetric reference, and charge and spin do
not couple (that coupling vanishes by spin symmetry). For LDA this is
the unique limit of the full kernel. Trial densities are symmetrized
(cf. #225).

Validation

H3 / cc-pVDZ (h3_blyp_cc-pvdz_ssf_gks.bin density), LDA_X+LDA_C_PW,
UltraFine, against central differences of the GKS potential. Relative max
error per block as the step halves:

reference	s	z	y	x
noncollinear	1.6e-7 -> 3.9e-8	2.9e-7 -> 7.2e-8	3.2e-7 -> 7.9e-8	9.9e-8 -> 2.5e-8
m = 0 (all points in the limit branch)	1.3e-7 -> 3.1e-8	2.8e-7 -> 6.9e-8	3.2e-7 -> 8.1e-8	9.9e-8 -> 2.5e-8

sum_X Tr(X F_X[Y]) = sum_X Tr(Y F_X[X]) to 1e-14.

Not included

• GGA throws. The Scalmani-Frisch gamma map used by the GKS GGA
  potential gives a kernel term (v_sigma++ - v_sigma--)/(2S) with
  S = |(grad rho_s . grad m_J)_J|, which diverges wherever S -> 0 at
  finite |m| -- including every density critical point (grad rho_s = 0).
  It needs a regularization of its own. In addition, the |m| -> 0 limit
  of the GGA kernel depends on how m vanishes (through grad|m|/|m|).
• meta-GGA: as for the GKS potential, there is no noncollinear
  kinetic-energy-density map yet.

Related, not changed here: below gks_dtol the GKS GGA potential uses
H = sign/3, whereas the gamma it evaluates implies H = s/|s|.

🤖 Generated with Claude Code

https://claude.ai/code/session_0135evJ9zgNL1y8U6T9bQ3UT