Linear response: CPHF/CPKS

A surprising number of seemingly unrelated questions — “is this SCF solution actually a minimum, or a saddle point?”, “how does the energy curve near a stationary point second-order (the Hessian)?”, “what Newton step converges SCF fastest?”, “at what photon energy does this molecule absorb light?” — all reduce to the same linear-algebra object: the action of the electronic (orbital) Hessian on a trial orbital rotation, \(\sigma=H\kappa\). qc-rs builds this action once, as a single injectable FockResponseKernel, and every consumer above is a thin wrapper around one of a small number of generic solvers applied to it. This chapter derives that unifying object and its four consumers, grounded in .design/51scf-linear-response.md.

The unifying object: \(A\) and \(B\)

Parametrize a trial orbital rotation the same way the SCF convergence chapter did, \(\mathbf C(\kappa)=\mathbf C\exp(\kappa)\) restricted to the non-redundant occupied-virtual block. The second-order expansion of the energy in \(\kappa\) defined the orbital Hessian \(H_{ai,bj}\) there; linear response theory just gives that same object (and its close relative) standard names and shows it governs far more than SCF convergence. For closed-shell (RHF/RKS) singlet rotations, in real canonical MOs with occupied \(i,j\) and virtual \(a,b\):

\[ A_{ia,jb} = \delta_{ij}\delta_{ab}(\varepsilon_a-\varepsilon_i) + 2(ia|jb) - (ij|ab)\ \bigl[+\,2f^{xc}_{ia,jb}\bigr], \qquad B_{ia,jb} = 2(ia|bj) - (ib|aj)\ \bigl[+\,2f^{xc}_{ia,bj}\bigr], \]

where the bracketed \(f^{xc}\) term (the exchange-correlation kernel, the second functional derivative of \(E_{\text{xc}}\)) is present for KS references and absent for pure HF. Two combinations of \(A,B\) recur everywhere: \(A+B\) is exactly the real-orbital-rotation Hessian from the SCF convergence chapter — the same object that drives the augmented-Hessian/QC-SCF and TRAH methods there — while \(A-B\) governs a different physical question (triplet/spin-flip response, dropping the Coulomb-like term and keeping only exchange). Which combination a given problem needs, and whether it needs an eigenvalue or a linear solve, is summarized below:

Consumer	Operator needed	Problem solved
Internal stability	\(A+B\)	lowest eigenvalue \(<0\) \(\Rightarrow\) unstable (real-orbital-rotation direction)
External/triplet stability	\(A-B\) (or the triplet kernel)	lowest eigenvalue \(<0\) \(\Rightarrow\) spin-flip/symmetry-breaking instability
QC-SCF / augmented Hessian / TRAH	\(A+B\)	lowest eigenpair of the bordered (augmented) matrix — already derived
CPHF/CPKS (nuclear/field response)	\(A+B\)	linear solve \((A+B)U=-b\) for a perturbation-specific right-hand side \(b\)
TDA (Tamm-Dancoff excited states)	\(A\) alone	lowest few eigenpairs, \(AX=\omega X\)
Full TDDFT/RPA	\(A,B\) jointly	\((A-B)(A+B)Z=\omega^2Z\) (Hermitian form)

Every row needs the same underlying \(A\)/\(B\) action — only the combination and the solver differ. This is the entire architectural point: build \(\sigma=(A{+}B)\kappa\) (or \(A\kappa\), or \((A{-}B)\kappa\)) once, as a single reusable primitive, and every consumer becomes a thin wrapper.

Building \(\sigma\) without ever forming \(A\) or \(B\)

Explicitly constructing \(A\) or \(B\) as a dense \((n_{\text{occ}}n_{\text{vir}})\times (n_{\text{occ}}n_{\text{vir}})\) matrix is exactly the kind of \(O(N^4)\)-scaling object the qc-hpc Workspace rules forbid materializing. Instead, \(\sigma\)-build evaluates the action of \(A+B\) on a trial vector \(X\) as a single response Fock build — a Fock-like contraction on a one-index transition density, reusing the exact same \(J/K\) (and, for KS, the same XC grid) machinery the ordinary SCF Fock build already uses:

Algorithm 4 (\(\sigma=(A{+}B)X\) via one response Fock build)

Input: trial vector \(X\) (shape \(n_{\text{vir}}\times n_{\text{occ}}\)), occupied/virtual MO coefficients \(C_o,C_v\), orbital energies \(\varepsilon\). Output: \(\sigma=(A+B)X\), same shape as \(X\).

1. Form the AO transition density \(P = C_v\,X\,C_o^{\mathsf T}\) (need not be symmetric in general; for the \(A+B\) combination specifically, use the symmetrized \(P+P^{\mathsf T}\) — this is what lets \(A+B\) collapse to a single symmetric response-Fock build rather than two).

2. Build the response Fock \(G[P]\) via the injected FockResponseKernel: \(2J[P]-a_xK[P]\) for HF/hybrid, with \(+2f^{xc}[P]\) added for KS (evaluated on the same DFT quadrature grid, just with libxc’s second functional derivative instead of the first).

3. Transform back to the occupied-virtual block: \(\sigma_{\text{2e}} = C_v^{\mathsf T}G[P]\,C_o\).

4. Add the diagonal (orbital-energy-difference) piece: \(\sigma = \sigma_{\text{2e}} + (\varepsilon_a-\varepsilon_i)X_{ia}\).

The cost of one \(\sigma\)-build is exactly one Fock build — the same \(J/K\) engine (incore/RI/direct, whichever the SCF itself used), the same Schwarz screening, the same threaded XC grid loop for KS. This is the single most important architectural decision in qc-rs’s response infrastructure: any future speedup to \(J/K\) or the XC grid loop benefits stability analysis, QC-SCF, CPHF, and TDA/TDDFT automatically, with no separate optimization work, because they all route through the identical Fock-build machinery — qc-rs deliberately does not give the response engine its own integral path. Multiple simultaneous trial vectors (Davidson’s block iterations, CPHF’s multiple right-hand sides — one per perturbation) batch naturally: stack several transition densities and do one batched Fock build instead of many separate ones, amortizing integral/grid access across all of them.

Reference coverage: one abstraction, three rotation structures

The response formalism must cover six references (RHF/UHF/ROHF × HF/KS), but it factors cleanly into two independent axes — spin structure (which determines the non-redundant rotation space and how \(\sigma\) is assembled) and kernel family (whether \(f^{xc}\) is added) — so qc-rs implements three rotation structures crossed with two kernel families rather than six separate paths:

Reference	Non-redundant rotation space	Spin channels	Kernel
RHF/RKS	occupied–virtual (spatial orbitals)	1 (singlet/triplet separated by kernel choice)	RestrictedSinglet/Triplet (+\(f^{xc}\) for RKS)
UHF/UKS	(\(\alpha\) occ–vir) \(\oplus\) (\(\beta\) occ–vir)	2 (\(\alpha\)/\(\beta\) coupled through \(J\) and \(f^{xc}_{\sigma\sigma'}\))	Unrestricted (+\(f^{xc}_{\sigma\sigma'}\) for UKS)
ROHF/ROKS	(closed–open) \(\oplus\) (closed–virtual) \(\oplus\) (open–virtual)	1 set, 3 spaces	a constrained spin-orbital Hessian + the ROHF coupling operator

Kohn-Sham-ification is structurally identical across all three: each HF-reference kernel simply gains an \(f^{xc}\) term (same-spin and cross-spin blocks for UKS), with the underlying rotation structure completely unchanged — RKS is RHF plus \(f^{xc}\), UKS is UHF plus \(f^{xc}_{\sigma\sigma'}\), ROKS is ROHF plus \(f^{xc}\). ROHF/ROKS is the genuinely hardest case (three coupled rotation spaces instead of one or two decoupled ones) and needs particular care in how half-occupied orbitals enter \(f^{xc}\), but it does not need a fourth architectural pattern — it needs a more careful accounting within the same \(\sigma\)-build abstraction.

External potentials (ECP, PCM, point charges/solvent) fold in as an organizing principle rather than a special case: \(\sigma\)-build is a sum of “the response of each density-dependent term,” so a density-independent external potential (a fixed ECP or point charge) contributes nothing extra to \(\sigma\) at all — its energy contribution is already baked into the orbital energies and converged density that the ordinary SCF Fock build used, and only genuinely density-dependent environment terms (like the PCM reaction field, which responds to any density including a transition density) need to add their own response contribution to the \(\sigma\)-build.

Four generic solvers, one shared engine

Once \(\sigma\)-build exists as an opaque closure, every consumer reduces to calling one of a small, reusable set of solvers on it — none of these solvers know or care whether \(\sigma\) came from RHF, UKS, or anything else:

• davidson_lowest — the lowest \(k\) eigenpairs of a large sparse-action operator via block Davidson with diagonal preconditioning. Used for stability analysis (\(k=1\), on \(A+B\) or \(A-B\)), for TDA (\(k=\) the number of excited states wanted, on \(A\) alone), and internally by the augmented-Hessian solver.

• augmented_hessian_step — the augmented-Hessian/QC-SCF eigenproblem (the bordered matrix’s lowest eigenpair, giving the automatic level shift), built on top of davidson_lowest.

• linear_solve — preconditioned conjugate gradient for the symmetric positive-definite system \((A+B)X=\text{rhs}\). This is the CPHF/CPKS solver, used identically whether the perturbation is a nuclear displacement (feeding the analytic Hessian), an external electric field (polarizability), or any other one-index perturbation.

• rpa_solve — the Hermitian reduction of the coupled \((A-B)(A+B)Z=\omega^2Z\) eigenproblem for full TDDFT/RPA (as opposed to TDA’s simpler \(A\)-only eigenproblem).

Internal stability analysis is davidson_lowest(apply_aplusb, diag_precond, k=1): a negative lowest eigenvalue of \(A+B\) means the current SCF solution is a saddle point in orbital-rotation space, not a true minimum, and its eigenvector gives the direction to rotate the orbitals toward the lower solution (the same Cayley-rotation machinery SOSCF uses). External/triplet stability is the same call on \(A-B\) or the dedicated triplet kernel instead.

Verified example — water/STO-3G stability analysis, both internal (real-orbital) and external (spin-flip) channels:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m = qc.chk.new(atom=water, ao="sto-3g", unit="angstrom").scf(ref="r", stability=True).run()
m.scf.energy, m.scf.converged
# (-74.96294665653868, True)
m.scf.stability
# {'internal_stable': True, 'internal_eigenvalue': 0.5236790238088829,
#  'external_stable': True, 'external_eigenvalue': 0.36925930146991304, 'stable': True}

Both the internal (\(A+B\)) and external (\(A-B\)) lowest eigenvalues are positive here — this restricted closed-shell solution is a genuine local minimum in both the real-orbital-rotation sense and the spin-symmetry-breaking sense, exactly the two questions this section’s two eigenvalue rows answer.

CPHF/CPKS: what makes it different from stability/QC-SCF

Stability analysis and QC-SCF only ever need \(A+B\)’s eigenvalues — they ask “in which direction, if any, does the energy decrease?” CPHF/CPKS asks a different question: “given a specific external perturbation (a nuclear displacement, an electric field), how do the orbitals respond, to first order?” This is a linear solve, not an eigenvalue problem, because the response \(U\) to a perturbation is uniquely determined (assuming \(A+B\) has no zero eigenvalue, i.e. the reference is stable) rather than extremal:

\[ (A+B)_{ai,bj}\,U_{bj} = -b_{ai}, \]

where \(b\) is a perturbation-specific right-hand side built from derivative integrals of whatever is being perturbed. For a nuclear-coordinate perturbation this \(b\) is built from the same first-derivative (gradient-level) integrals the analytic gradient chapter introduced — this is precisely the “genuinely needs orbital response” case that chapter flagged as out of its own scope, and the next chapter, on the analytic Hessian, derives that right-hand side and the resulting fold-back in full. For an external electric field, \(b\) is instead built from the AO dipole integrals — same solver, same \(A+B\) operator, different perturbation-specific right-hand side; this is what a static dipole polarizability calculation ultimately reduces to. Either way, linear_solve is called on the identical apply_aplusb closure stability analysis and QC-SCF already use — this is the sense in which “CPHF reuses the SCF convergence machinery” is not a loose analogy but a literal code-sharing fact.

Where excited-state response stands today

TDA and full TDDFT/RPA are architecturally accounted for in this design (davidson_lowest(apply_a, ...) and rpa_solve, respectively) — the same \(\sigma\)-build abstraction this chapter derives is exactly what a real implementation would consume. However, as of this writing qc.td(...) is a mock step in qc-rs: it accepts the same pending-step syntax and records a placeholder converged=True result, with no real eigenvalue solve behind it yet. This mirrors the post-HF correlation chapter’s honest accounting of cc2/caspt2/nevpt2 — the architectural seam is real and load-bearing (internal stability analysis already exercises the identical apply_a/apply_aplusb machinery TDA/TDDFT would need), but the excited-state solver itself has not yet been wired up to it.

Exercise 11

1. The verified stability example reports two eigenvalues, one for \(A+B\) and one for \(A-B\). Explain in one sentence why a restricted closed-shell reference needs to check both, rather than just one, to be confident it is a genuine minimum.

2. CPHF’s right-hand side \(b\) differs completely between a nuclear-displacement perturbation and an electric-field perturbation, yet both are solved with the exact same linear_solve(apply_aplusb, ...) call. What is it about the left-hand side operator \(A+B\) that makes it perturbation-independent, and why does that make code reuse possible here?

3. Why can a fixed (density-independent) ECP or point-charge external potential be completely ignored when building the response kernel \(\sigma\), even though it is very much not ignored when building the ordinary SCF Fock matrix or energy?

Solution to Exercise 11

1. \(A+B\) governs stability against real orbital rotations that preserve the closed-shell (spin-restricted) character of the reference — a negative eigenvalue there means a lower-energy restricted solution exists. \(A-B\) (or the dedicated triplet kernel) governs stability against symmetry-breaking rotations that would lower the energy only by leaving the restricted-singlet space (e.g. into an unrestricted, spin-polarized solution) — a negative eigenvalue there means the true lowest-energy solution is not even restricted-closed-shell at all. These are genuinely different failure modes (one within the ansatz, one escaping it), so a solution can pass one check and fail the other; checking only one leaves an entire class of instability undetected.

2. \(A+B\) is built entirely from the converged, unperturbed reference — its canonical orbital energies, its converged density, its (possibly hybrid/range-separated) exchange-correlation kernel — none of which depend on what kind of perturbation is being applied. The perturbation enters only through the right-hand side \(b\), which is built from derivative integrals specific to that perturbation (nuclear first-derivative integrals for a geometry perturbation, AO dipole integrals for a field). Because the expensive, iterative part of the solve (repeated \(\sigma\)-builds against the same operator) is perturbation-independent, the identical linear_solve/apply_aplusb machinery serves any perturbation just by swapping out the cheap right-hand-side construction.

3. The response kernel \(\sigma\)-build is, by construction, a sum of “the response of each density-dependent term” to a change in the transition density \(P\) — it differentiates each term in the Fock build with respect to the density. A fixed ECP operator or a fixed point charge contributes a term to the Fock matrix that does not depend on the density at all (it is added once, from the converged geometry, and stays constant across SCF iterations), so its derivative with respect to any density — including a transition density — is exactly zero. Its effect on the energy and the converged orbitals is already fully accounted for by the time you reach the response step; only genuinely density-dependent environment terms (like PCM’s density-responsive reaction field) need their own contribution to \(\sigma\).

The next chapter picks this machinery up directly: the analytic Hessian is, structurally, nothing more than “\(3N_{\text{atom}}\) CPHF solves plus a fold-back contraction” — every formula in this chapter reappears there, specialized to the nuclear-coordinate perturbation.