Handover: Phase 17e — polarised GGA f_xc + UKS-GGA second-order SCF

Filed: 2026-05-20, molecular-methods chat Branch: claude/exciting-goldberg-029ff6 (worktree), targeting main Status: ✅ COMPLETE — the 4-commit stack landed on origin/main (pushed 3d13fcf5..403067d7). See § “Status snapshot” below.

TL;DR

The only thing blocking UKS-GGA / hybrid-GGA second-order SCF (Newton + TRAH) was the missing spin-polarised GGA exchange-correlation kernel (f_xc = second derivatives of the XC energy density). Every other second-order SCF path already ships: RHF/UHF Newton/SOSCF/TRAH, RKS LDA+GGA Newton/SOSCF/TRAH, UKS-LDA Newton/SOSCF/TRAH. This work closes the last gap — “Phase 17e” in the roadmap / xc_kernel.hpp / uks.hpp comments.

Three code milestones + one docs milestone:

  1. M1 — Functional::eval_polarised_gga_fxc: libxc xc_gga_fxc wrapper for spin = 2 (the 15 second-derivative pieces).

  2. M2 — PolarisedGGAXCKernelBuilder + factories: the AO-basis open-shell GGA W^XC matvec.

  3. M3 — wire the dispatcher into the UKS Newton / TRAH drivers.

  4. M4 — roadmap / CHANGELOG / user-guide / this handover.

Why this is the right next item

docs/roadmap.md and the in-tree comments (cpp/include/vibeqc/uks.hpp, cpp/include/vibeqc/xc_kernel.hpp, cpp/src/xc.cpp::eval_polarised_lda_fxc) all explicitly name “Phase 17e” as the gate. Before this work:

  • make_polarised_lda_xc_kernel_builder raised on any non-LDA functional.

  • uks.cpp Newton + TRAH called the LDA-only builder, so UKS-GGA Newton/TRAH raised "only LDA functionals are supported".

  • Functional had no polarised GGA second-derivative method at all (eval_polarised_lda_fxc was LDA-only; the GGA signature was “reserved for Phase 17e”).

Design

libxc layout (M1)

xc_gga_fxc(p, np, rho, sigma, v2rho2, v2rhosigma, v2sigma2) for spin = 2 returns, per grid point:

  • v2rho2 — 3: (aa, ab, bb)

  • v2rhosigma — 6: (a·aa, a·ab, a·bb, b·aa, b·ab, b·bb) (full 2×3)

  • v2sigma2 — 6: (aa·aa, aa·ab, aa·bb, ab·ab, ab·bb, bb·bb) (upper triangle of the symmetric 3×3)

= 15 unique pieces, bundled into struct PolarisedGGAFxc (declared in xc.hpp). LDA components inside a composite (e.g. VWN5 inside B3LYP) contribute only to the v2rho2_* block. Meta-GGA raises.

Open-shell GGA W^XC matvec (M2)

Static UKS GGA potential, spin σ (mirrors uks.cpp::build_xc):

V_σ,μν = Σ_g w_g [ v_ρσ χ_μχ_ν + U_σ·∇(χ_μχ_ν) ]
U_α = 2 v_σαα ∇ρ_α + v_σαβ ∇ρ_β       (σαβ cross channel, factor 1)
U_β = 2 v_σββ ∇ρ_β + v_σαβ ∇ρ_α

Linearising about (D_α^used, D_β^used) gives the perturbations δρ_σ, δ∇ρ_σ, δσ_αα, δσ_αβ, δσ_ββ; the five first-derivative perturbations δv_* are the symmetric 5×5 f_xc Hessian contracted against that 5-vector. δV_σ then splits into the same three contributions as the closed-shell GGA builder (UnpolarisedGGAXCKernelBuilder):

  • T1 — δv_ρσ · χ_μχ_ν

  • T2 — δv_σ in the flow prefactor (∇ρ pinned at the reference)

  • T3 — δ∇ρ in the flow vector (v_σ pinned at the reference)

The αβ coupling enters through v2rho2_ab, the v2rhosigma_* / v2sigma2_* cross pieces, and the σ_αβ cross term. spin_block() assembles one spin given its “self” (σσ, factor 2) and “cross” (αβ, factor 1) channels.

Dispatch (M3)

New unified factory make_polarised_xc_kernel_builder dispatches on func.kind(): LDA → existing make_polarised_lda_xc_kernel_builder, GGA / hybrid-GGA → new make_polarised_gga_xc_kernel_builder, meta-GGA → raises. uks.cpp Newton + TRAH now call the dispatcher. make_polarised_lda_* kept for back-compat (still bound, still used by existing tests).

Correctness witnesses

  • M1 — eval_polarised_gga_fxc vs central-difference of eval_polarised’s first derivatives (the 5×5 Hessian), all 15 pieces, on PBE / BLYP / B3LYP / PBE0. tests/test_xc.py::test_polarised_gga_fxc_*.

  • M2 — W^XC_σ[δD_α,δD_β] vs central-difference of the static per-spin UKS GGA potential V_xc,σ(D_α,D_β). tests/test_xc_kernel.py::test_polarised_gga_kernel_*.

  • M3 — UKS Newton + TRAH converge LDA / PBE / B3LYP on the OH· radical to the physical broken-symmetry doublet (⟨S²⟩ ≈ 0.75); Newton/TRAH agree, and where plain DIIS converges (LDA, PBE) Newton matches it. Energy-parity asserts use abs=1e-6 — different converged second-order methods sit at slightly different depths within conv_tol_grad = 1e-6 on the radical’s flat direction (a wrong minimum would differ by tens of mHa). For B3LYP plain DIIS collapses to the ⟨S²⟩ = 0 closed-shell solution — second-order SCF recovers the doublet. tests/test_uks_second_order.py.

Status snapshot

All four milestones are code-complete, built, green, and landed on origin/main (the 4-commit stack pushed as 3d13fcf5..403067d7). Commits are referenced by subject below (SHAs shift on rebase):

Milestone

State

Lands in commit

M1 libxc wrapper

done, 9 test cases

xc: polarised GGA fxc kernel …

M2 kernel builder

done, 13 test cases

xc: polarised GGA fxc kernel …

M3 wire UKS Newton/TRAH

done, 8 test cases

uks: wire UKS-GGA/hybrid-GGA …

M4 docs / roadmap

done

docs: Phase 17e …

Commit stack (oldest → newest):

fix(bindings): escape nested quotes in CCSD docstrings
xc: polarised GGA fxc kernel + open-shell GGA W^XC matvec (Phase 17e)
uks: wire UKS-GGA/hybrid-GGA second-order SCF (Newton + TRAH)
docs: Phase 17e — UKS-GGA second-order SCF (roadmap, changelog, guide)

origin/main was broken on arrival. cpp/src/bindings.cpp had unescaped nested double-quotes in the CCSD aux_basis / triples_memory_mode docstrings — clang rejects it, so the whole extension fails to build. The fix(bindings): … commit fixes it and is ordered first so every commit in the stack ships green. This is a pre-existing upstream defect (a CCSD-bindings commit in the 61-commit window since 2e8970b) — flag to the maintainer; CI evidently did not gate it.

Commit grouping — M1 + M2 ship as one commit (xc: polarised GGA fxc kernel …) because cpp/src/bindings.cpp is touched by both and cannot be split into two individually-compiling commits.

Verified green on the rebased base (origin/main @ 755acef4 + the 4-commit stack, rebuilt clean):

  • test_xc.py 44 pass, test_xc_kernel.py 13 pass (+1 skip), test_uks_second_order.py 8 pass.

  • Regression: test_newton.py + test_trah.py + test_rks_second_order.py + test_uhf_newton.py + test_hessian_analytic_uks.py — 46 pass, 1 skip, 1 xfail (the xfail is the pre-existing 17e-2 nuclear-Hessian GGA case, not in scope here; the skip is the meta-GGA placeholder).

One test tolerance was adjusted after the rebase: upstream landed feat(scf): TRAH proper (755acef4) in the 61-commit window, which changed TRAH’s stopping depth. Newton and TRAH now converge to the same UKS minimum but to slightly different depths within conv_tol_grad = 1e-6 (~2.6e-7 Ha apart for B3LYP on OH·), so the Newton-vs-TRAH energy-parity assert was relaxed from abs=1e-8 to abs=1e-6 (“agrees to the SCF gradient convergence tolerance”).

Files touched

  • cpp/include/vibeqc/xc.hpp — struct PolarisedGGAFxc + eval_polarised_gga_fxc declaration.

  • cpp/src/xc.cpp — eval_polarised_gga_fxc implementation.

  • cpp/include/vibeqc/xc_kernel.hpp — make_polarised_gga_* + make_polarised_xc_kernel_builder declarations.

  • cpp/src/xc_kernel.cpp — PolarisedGGAXCKernelBuilder + factories.

  • cpp/src/bindings.cpp — bindings for the new Functional method + the two new factory functions.

  • cpp/src/uks.cpp — Newton + TRAH call sites use the dispatcher.

  • cpp/include/vibeqc/uks.hpp — Phase 17e comment refresh.

  • python/vibeqc/__init__.py — re-export the two new factories.

  • tests/test_xc.py, tests/test_xc_kernel.py, tests/test_uks_second_order.py — new coverage.

Next steps / priorities

Phase 17e (this work) is done and on main. Nothing else is required for it.

Natural follow-ups, in rough priority order:

  1. Close the 17e-2 analytic-Hessian xfail. eval_polarised_gga_fxc is the exact primitive the analytic nuclear UKS Hessian needs for hybrid-GGA (B3LYP-UKS) — see tests/test_hessian_analytic_uks.py’s xfail. Wiring it into the CPKS response there closes the last FD-fallback in the analytic-Hessian path.

  2. Polarised meta-GGA f_xc (τ-dependent) — unblocks UKS-MGGA second-order SCF (TPSSh, M06-2X, r²SCAN UKS). Needs libxc xc_mgga_fxc (spin = 2) plus the τ second-derivative block in the kernel builder. Today UKS-MGGA Newton/TRAH raise a pointer.

  3. Periodic-Γ UKS second-order — the Γ-only periodic UKS driver could reuse make_polarised_xc_kernel_builder once a periodic grid + AO table are threaded through, mirroring how molecular UKS does it now.

Out of scope / follow-ups

  • Phase 17e-2 — the analytic nuclear UKS Hessian (vibrational frequencies) for GGA also needs σ-coupled polarised second derivatives. eval_polarised_gga_fxc is exactly the primitive that unblocks it; tests/test_hessian_analytic_uks.py has an xfail for the GGA case. NOT done here — wiring CPKS for nuclear displacements is a separate, larger deliverable.

  • Polarised meta-GGA f_xc (τ-dependent) — UKS-MGGA second-order SCF still routes to DIIS / quadratic-fallback; the factory raises.

Assumptions / decisions

  • libxc v2rhosigma ordering taken as rho_index*3 + sigma_index (the standard libxc spin-polarised layout); validated by the M1 finite-difference test against eval_polarised — any index slip shows up as an O(1) relative error.

  • make_polarised_lda_xc_kernel_builder kept rather than removed: it is bound to Python and used by existing tests; the dispatcher delegates to it for LDA.

  • Worktree build: third_party/ is shared from the parent checkout (<vibe-qc-checkout>/third_party) since it is gitignored; a worktree-private .venv was created per the standing worktree-venv rule.