# Multireference methods: an outlook

Every method this manual has derived so far — Hartree-Fock, Kohn-Sham DFT, RI-MP2 — starts from a single
Slater determinant (or single Kohn-Sham determinant) and treats everything beyond it as a perturbative or
density-functional correction. That single-reference picture fails qualitatively, not just
quantitatively, whenever more than one electron configuration contributes comparably to the true
wavefunction — stretched bonds approaching dissociation, biradicals, many transition-metal complexes, and
excited states near a conical intersection. This closing chapter of Part II derives what a genuine
multireference treatment (CASSCF, and the CASPT2/NEVPT2 perturbative corrections built on top of it) would
require, and states plainly what qc-rs has actually built toward that goal versus what remains a
placeholder — grounded in `.design/55casscf-fci-spine.md` and the current, honestly-mock `execute_casscf`/
`execute_fci` in `qc-workflow`.

## Why single-reference methods break down

Hartree-Fock and Kohn-Sham DFT both assume the true ground state is *dominated* by a single determinant —
formally, that the CI expansion coefficient of the reference determinant is close to 1 and every other
determinant contributes only a small perturbative correction. This assumption fails whenever two or more
determinants become comparably important. The textbook example is stretching an H$_2$ bond toward
dissociation: at equilibrium, the bonding-orbital-doubly-occupied determinant dominates overwhelmingly, but
as the bond stretches toward two separated hydrogen atoms, the bonding and antibonding orbitals become
degenerate, and the true wavefunction needs *both* the bonding-doubly-occupied and
antibonding-doubly-occupied determinants with comparable weight to describe two independent, singly-occupied
atoms correctly. A single RHF determinant cannot represent this — it either stays artificially
bonding-dominated (wrong dissociation limit, a qualitatively wrong potential energy surface) or a UHF
solution breaks spin symmetry to compensate (a documented, structural failure mode [the SCF convergence
chapter](scf-convergence-theory.md) and [linear-response chapter](linear-response-cphf.md) both flag via
the external/triplet stability eigenvalue). The same qualitative breakdown recurs for biradicals (two
nearly-degenerate frontier orbitals, each singly occupied in the true ground state), many open-shell
transition-metal complexes (near-degenerate d-orbitals), and any excited state near a conical intersection
with the ground state.

## CASSCF: the active-space Hamiltonian and the CI expansion

The **Complete Active Space Self-Consistent Field (CASSCF)** method addresses this by partitioning the
orbital space into three groups — **core** (always doubly occupied, never correlated explicitly), **active**
(a small, chemically chosen set spanning the near-degenerate orbitals responsible for the multireference
character), and **virtual** (always empty) — and solving *simultaneously* for the best orbitals *and* the
best CI expansion coefficients of every determinant that distributes the active electrons among the active
orbitals in every possible way. Freezing the core orbitals as doubly occupied lets the full Hamiltonian be
projected onto a much smaller **active-space Hamiltonian**,

$$
\hat H_{\text{cas}} = E_{\text{core}} + \sum_{tu}h^{\text{eff}}_{tu}\hat E_{tu} +
\tfrac12\sum_{tuvw}(tu|vw)\,\hat e_{tuvw}, \qquad
h^{\text{eff}}_{tu} = h_{tu} + \sum_i^{\text{core}}\bigl[2(tu|ii)-(ti|iu)\bigr],
$$

where $t,u,v,w$ index only the active orbitals, $\hat E_{tu}$ is the spin-free singlet excitation operator,
and $E_{\text{core}}$ folds the nuclear repulsion and every core-orbital contribution into a single
constant. This active-space Hamiltonian is genuinely small (its dimension scales only with the *active*
orbital count, not the full AO/MO count), which is exactly what makes an essentially-exact diagonalization
— **Full CI within the active space** — computationally tractable even though Full CI over the whole
orbital space would be impossible for any but the smallest molecules.

The CI expansion itself is characterized entirely by two numbers a user chooses chemically: $n_{\text{cas}}$
(the number of active orbitals) and the active electron count, and its result is a full configuration
expansion $|\Psi\rangle=\sum_IC_I|\Phi_I\rangle$ over every determinant $\Phi_I$ distributing those
electrons among those orbitals — capturing, by construction, any near-degeneracy among the active orbitals
exactly, rather than perturbatively. From this CI solution follow the **one- and two-particle reduced
density matrices** $D_{tu}=\langle\Psi|\hat E_{tu}|\Psi\rangle$ and
$d_{tuvw}=\langle\Psi|\hat e_{tuvw}|\Psi\rangle$, which are what actually drives everything downstream —
the CASSCF energy itself, $E=E_{\text{core}}+\sum_{tu}h^{\text{eff}}_{tu}D_{tu}+\tfrac12\sum_{tuvw}
(tu|vw)d_{tuvw}$, and, as the next section shows, the orbital optimization step.

## Orbital optimization reuses exactly the machinery already derived

CASSCF is a *simultaneous* orbital-and-CI optimization, iterated as coupled macro/micro cycles: at each
macro iteration, solve the CI problem for the current orbitals to get updated RDMs, then take a
step in orbital-rotation space using those RDMs, then re-solve the CI problem at the new orbitals, and so
on until both the orbital gradient and the CI residual converge. The orbital-rotation step is the genuinely
interesting part for this manual, because **it reuses, unchanged, exactly the same augmented-Hessian
machinery [the SCF convergence theory chapter](scf-convergence-theory.md) and [the linear-response
chapter](linear-response-cphf.md) already derived** for ordinary SCF stability analysis and CPHF — this is
not a coincidental resemblance, it is a deliberate architectural decision (`.design/55` explicitly notes:
"the AH/Davidson/rotation machinery is not reimplemented"). The non-redundant orbital-rotation blocks now
include core-active and active-virtual pairs in addition to the ordinary occupied-virtual ones (active-active
rotations are redundant — they only relabel determinants within the CI expansion, not genuinely new
orbitals), the orbital gradient and generalized Fock matrix are built from the RDMs rather than a simple
density matrix, and the Hessian-vector product needed by the augmented-Hessian solver is obtained by a
**one-index transformation** — contracting the trial rotation $\kappa$ into one integral index at a time,
re-evaluating the generalized Fock at frozen RDMs, and reading off $(\mathbf H\kappa)_{pq}=2(F^\kappa_{pq}-
F^\kappa_{qp})$ — never forming the orbital Hessian explicitly, exactly the same matrix-free discipline
`apply_aplusb` follows for ordinary SCF. The lowest eigenpair of the same bordered (augmented) matrix from
[the SCF convergence chapter](scf-convergence-theory.md) gives the orbital step, with the same built-in,
automatic level shift.

## CASPT2/NEVPT2: perturbation theory on top of a multireference wavefunction

A converged CASSCF wavefunction already captures the *qualitative* multireference physics (the correct
near-degenerate mixing), but it is missing the same kind of dynamic electron correlation that [ordinary
RI-MP2](post-hf-correlation.md) recovers on top of a single-reference HF wavefunction — correlation among
electrons that are *not* strongly coupled by near-degeneracy, but still contribute a real energetic
correction. **CASPT2** and **NEVPT2** are both second-order perturbation theories built *on top of* a
converged CASSCF reference, generalizing Møller-Plesset perturbation theory to a multi-determinantal
zeroth-order wavefunction — the same basic idea as MP2 (partition the Hamiltonian, treat the reference as
zeroth order, compute a second-order energy correction), but with a genuinely harder zeroth-order problem:
where MP2's zeroth order is one determinant with trivially known excitation amplitudes, a multireference
perturbation theory's zeroth order is an entire CI expansion, and the perturbation must be applied
consistently across every determinant in it. CASPT2 and NEVPT2 differ mainly in how they partition the
Hamiltonian into zeroth-order and perturbation pieces (CASPT2's Fock-operator-based partitioning versus
NEVPT2's genuinely two-electron Dyall Hamiltonian), each with known trade-offs in intruder-state behavior
and computational cost that are outside this outlook's scope to derive in full.

What *is* directly relevant here is a concrete complexity fact the design work already surfaced: these
perturbative corrections need **higher-order reduced density matrices** than CASSCF's own energy and
orbital-optimization steps do — the CI-solver interface's design already reserves `make_rdm3` (the
three-particle RDM) for CASPT2 and `make_rdm123_4` (up to the four-particle RDM) for NEVPT2, exposed only
"lazily," on demand, precisely because these higher-order RDMs are substantially more expensive to build
and store than the 1-/2-RDMs CASSCF itself needs. This is the concrete computational reason CASPT2/NEVPT2
are harder engineering problems than CASSCF alone, not merely "more perturbation theory bolted onto the
same CASSCF output" — a real implementation needs new density-matrix machinery beyond what CASSCF's own
orbital optimization already requires.

## What qc-rs has actually built, honestly

**Today, `qc.casscf(...)`, `qc.fci(...)`, and `qc.lct(method="caspt2"/"nevpt2")` are all mock steps.** They
accept the full pending-step syntax this manual describes elsewhere, execute without error, and record a
deterministic placeholder energy and `converged=True` — but that energy is not a real CASSCF/FCI/CASPT2/
NEVPT2 calculation, and changing the active-space size changes nothing about it:

```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
m0 = qc.chk.new(atom=water, ao="sto-3g", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").run()

m_44 = m0.casscf(ncas=4, nelecas=4).run()
m_44.casscf.energy   # -48.076100000000004

m_66 = m0.casscf(ncas=6, nelecas=6).run()
m_66.casscf.energy   # -48.076100000000004 -- identical to the (4,4) active space above

m_44_again = m0.casscf(ncas=4, nelecas=4).run()
m_44_again.casscf.energy  # -48.076100000000004 -- same again
```

The energy is completely insensitive to the active-space size — a physically impossible result for a real
CASSCF calculation, and exactly the tell that this is a deterministic mock rather than genuine physics.
This is the same honest accounting [the post-HF correlation chapter](post-hf-correlation.md) and [the
linear-response chapter](linear-response-cphf.md) already gave for `cc2` and `qc.td(...)` respectively —
qc-rs is explicit and consistent about which pending-step syntax is backed by real physics (RHF/UHF/ROHF/
KS-DFT SCF, analytic gradients/Hessians where implemented, RI-MP2/SCS-MP2/SOS-MP2) and which is scaffolding
reserved for a feature that is designed but not yet built.

What genuinely *is* built, and is worth taking seriously as forward progress rather than dismissing
alongside the mock energies: the **architectural seam** this chapter derived — the `ActiveHamiltonian`
projection, the `CiSolver` trait (structurally mirroring PySCF's `fcisolver` protocol, deliberately not
ORZ's disk-and-string-dispatch design), and above all the decision to reuse [the linear-response
chapter](linear-response-cphf.md)'s exact augmented-Hessian/Davidson machinery for CASSCF's orbital
optimization rather than reimplementing it. When a real `qc-fci`/`qc-casscf` implementation lands, it slots
into infrastructure this manual has already derived in full, rather than needing a parallel, from-scratch
optimization engine.

:::{exercise}
:label: ex-multireference-outlook

1. The verified example shows `casscf(ncas=4, ...)` and `casscf(ncas=6, ...)` giving the exact same energy.
   Explain in one sentence why this specific observation — insensitivity to active-space size — is a
   stronger tell of a mock implementation than, say, an energy that merely looked "a bit too round" or
   suspiciously simple.
2. CASSCF's orbital-rotation step explicitly excludes active-active rotations from its non-redundant
   parametrization, while including core-active and active-virtual rotations. Explain in one sentence why
   rotating two active orbitals into each other is redundant specifically for CASSCF (as opposed to
   ordinary SCF, where there is no "active" space at all).
3. CASPT2 and NEVPT2 both need substantially higher-order reduced density matrices (3-particle, up to
   4-particle) than CASSCF's own energy/orbital-optimization steps require (1-/2-particle). Why does adding
   a *perturbative correction on top of* a CASSCF reference require higher RDM orders than the CASSCF
   reference calculation needed in the first place?
:::

:::{solution} ex-multireference-outlook
:class: dropdown

1. A "too round" or suspiciously simple-looking number could, in principle, be a genuine coincidence of a
   particular molecule/basis/active-space combination — it is weak evidence at best. Complete insensitivity
   to the active-space size, however, is not merely suspicious, it is a **direct logical contradiction**
   with what CASSCF actually computes: the active-space Hamiltonian, the CI expansion, and hence the
   energy are all explicitly functions of $n_{\text{cas}}$ and the active electron count by construction —
   a genuine CASSCF energy *must* change (generically, though not universally, decrease in magnitude) as
   the active space grows to include more correlation, so an energy that provably cannot respond to that
   input is definitive proof the calculation behind it is not the one being requested, not just suspicious
   in style.
2. In CASSCF, all determinants generated by distributing the active electrons among the active orbitals are
   already included in the CI expansion with their own optimized coefficients — rotating one active orbital
   into a linear combination with another active orbital does not access any electronic configuration that
   was not already reachable by adjusting those CI coefficients directly. The active-active block is
   therefore a redundant parametrization of the same wavefunction, not new variational freedom, which is
   exactly why it is excluded (in ordinary SCF, by contrast, there is no CI expansion doing this job at
   all — every occupied-virtual rotation genuinely changes the one-determinant wavefunction, so no rotation
   block is redundant there).
3. A second-order perturbative energy correction is, by the general structure of second-order perturbation
   theory, a matrix element of the perturbation operator between the zeroth-order reference and *doubly
   excited* (or more) intermediate states, contracted twice against reference quantities — one extra particle-
   rank of density matrix is needed for each additional power of the two-electron perturbation operator
   entering the energy expression relative to what the zeroth-order energy itself required. Since CASPT2/
   NEVPT2 apply a genuine two-electron perturbation operator on top of an already-two-electron zeroth-order
   Hamiltonian expectation value, the resulting energy expression needs density-matrix contractions of a
   correspondingly higher particle rank (3-particle for CASPT2, up to 4-particle for NEVPT2's more general
   Dyall-Hamiltonian partitioning) than the plain CASSCF energy, which only ever needs the reference
   wavefunction's own 1- and 2-RDMs.
:::

This closes Part II's derivation of every method this manual documents in depth. [Part III's guide
chapters](../20-guide/scf.md) return to the compact usage-focused presentation of these same ideas; every
digest there links back to the full derivation a chapter in this Part provides.