XC nuclear Hessian: basis, moving-grid and partition-weight classes

[susilehtola]

GauXC has nuclear gradients of the XC energy but no Hessian. This adds
eval_exc_hess(P, settings): the full 3N x 3N second derivative of
E_xc at fixed density matrix, RKS, for LDA, GGA and (tau-dependent)
meta-GGA functionals, on the host replicated integrator.

What is included

IntegratorSettingsEXC_HESS mirrors IntegratorSettingsEXC_GRAD:

• include_weight_derivatives = true (default): all three ways E_xc
  depends on the nuclei -- the basis functions ride their atoms, the grid
  points ride their parent atoms, and the partition weights change.
  Each task is accumulated in coordinates excluding its parent atom C
  (the task energy depends on R_A - R_C only), and C's rows and columns
  are restored by translational invariance. The weight class adds
  w'' e + w' e'^T + e' w'^T per point, with Becke or SSF weight
  derivatives in log form, which stay finite as the parent partition
  function goes to zero.
• include_weight_derivatives = false: the basis class alone (fixed
  grid, fixed weights), the analogue of the gradient's Hellmann-Feynman
  option.

The Pulay term is contracted with GEMMs (it is exactly bilinear in the
displaced basis rows), so its cost is BLAS-3 rather than
O(9 nbe^2 npts) scalar kernel calls.

Generated code

The per-point kernels (host/gauxc_hess_kernel.hpp) and their call
sites (host/gauxc_hess_call_*.inc) are generated by xckernel from the
symbolic nuclear Hessian, and should be regenerated rather than edited:

python -m xckernel.emitters.gauxcwriter --emit-dir src/xc_integrator/local_work_driver/host

Call-site arguments are bound by name from each kernel's own signature,
because GauXC's names (vgamma, ...) sort differently from Libxc's
(vsigma, ...), and a hand-written meta-GGA call site had silently
mis-bound five of seven arguments. The generator also proves, symbolically,
that the assembly recipe used here (per-atom rows contracted as an outer
product, plus the Pulay and same-atom terms) equals the unfactorised
Hessian.

Validation

Benzene / cc-pVDZ; LDA (SVWN), PBE, TPSS, LAK.

• Full Hessian vs central differences (Richardson) of
  eval_exc_grad with weight derivatives, displacing atom + basis and
  rebuilding the grid: diagonal and off-diagonal blocks agree to ~1e-12
  with Becke weights and ~1e-9 with SSF (PySCF0 grid). On UltraFine,
  Becke agrees to 8e-12. SSF agrees to 2e-7 there, which is the noise
  floor of the reference itself: its finite-difference Hessian is
  asymmetric by the same amount.
• H is symmetric to 3e-16 and satisfies the translational sum rule
  sum_B H[(A,x),(B,y)] = 0 to 3e-15.
• Basis class vs finite differences displacing only the basis
  function origins: agreement at the FD truncation level for all four
  functionals.

Limitations

• RKS only; host only (other integrators throw via a defaulted virtual,
  so no device stubs are needed).
• Laplacian-dependent meta-GGAs are rejected: they need fourth
  collocation derivatives.
• LKO weights throw when weight derivatives are requested.
• The Becke weight second derivative costs O(N^3) per grid point (every
  atom pair contributes); SSF skips inactive pairs and is much cheaper.

🤖 Generated with Claude Code

https://claude.ai/code/session_0135evJ9zgNL1y8U6T9bQ3UT