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 and linear-response chapter 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 and the linear-response chapter 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 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 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:

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 and the linear-response chapter 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’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 19

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 to Exercise 19

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 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.