# 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](scf-convergence-theory.md)
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](scf-convergence-theory.md) |
| 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](../30-hpc/threads-and-blas.md) 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:

:::{prf:algorithm} $\sigma=(A{+}B)X$ via one response Fock build
:label: alg-sigma-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](dft-xc-quadrature.md), 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](scf-convergence-theory.md)
  (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](analytic-hessian-thermo.md)), 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](scf-convergence-theory.md) 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:

```python
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](analytic-derivatives.md) 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](analytic-hessian-thermo.md), 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](post-hf-correlation.md)'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}
:label: ex-cphf-theory

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} ex-cphf-theory
:class: dropdown

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](analytic-hessian-thermo.md) 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.