# Post-Hartree-Fock correlation theory

Hartree-Fock gives each electron the *average* field of the others — its systematic error, the
**correlation energy**, is defined as $E_{\text{corr}}=E_{\text{exact}}-E_{\text{HF}}$ and is typically
$\sim1\%$ of the total energy but $\gg1\ \text{kcal/mol}$, i.e. chemically decisive. This chapter derives
Møller-Plesset second-order perturbation theory (MP2) — the correlation method qc-rs actually implements —
in its **RI/density-fitted** form, since that is the *only* form qc-rs ships (per `AGENTS.md`: "these
correlation methods are RI-C only"), grounded in `.design/25qc.lct.mp2.md` and
`crates/qc-ints-engine/src/trafo/mp2.rs`.

## Møller-Plesset perturbation theory: the setup

MP2 treats the Hartree-Fock Hamiltonian as the zeroth-order problem and the difference between the true
electronic Hamiltonian and the Fock operator sum as the perturbation:

$$
\hat H = \hat H_0 + \hat V, \qquad \hat H_0 = \sum_i\hat f(i), \qquad
\hat V = \hat H_{\text{elec}} - \hat H_0,
$$

where $\hat f$ is the Fock operator (whose eigenfunctions are exactly the canonical HF orbitals — this is
why HF, not an arbitrary reference, is the natural zeroth order). The zeroth-order energy is the sum of
occupied orbital energies, and the first-order correction recovers exactly the HF energy —
$E^{(0)}+E^{(1)}=E_{\text{HF}}$ — so **MP2 is the first correlation-recovering order**, the second-order
correction $E^{(2)}$. Standard Rayleigh-Schrödinger perturbation theory gives

$$
E^{(2)} = \sum_{n\ne0}\frac{|\langle\Psi_0|\hat V|\Psi_n\rangle|^2}{E_0^{(0)}-E_n^{(0)}}.
$$

Because $\hat V$ is a two-electron operator, only *doubly excited* determinants $\Psi_{ij}^{ab}$ (replacing
occupied orbitals $i,j$ with virtuals $a,b$) have a nonzero matrix element with $\Psi_0$ — singles vanish by
Brillouin's theorem (the HF orbitals already make the singles-HF coupling zero) and triples-and-higher
vanish because $\hat V$ only connects determinants differing by at most two spin-orbitals. Working through
the two-electron matrix elements and energy denominators (each virtual/occupied pair contributing
$\varepsilon_a-\varepsilon_i$ to the energy gap) gives the **closed-shell spin-adapted MP2 energy**:

$$
E_{\text{MP2}} = \sum_{i\le j}\sum_{ab}\frac{(ia|jb)\bigl[2(ia|jb)-(ib|ja)\bigr]}
{\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b},
$$

with $(ia|jb)$ the two-electron integral in chemists' notation over occupied ($i,j$) and virtual ($a,b$)
canonical MOs. This is exactly qc-rs's working equation (`rimp2_energy`'s doc comment), summed over
occupied pairs $j\le i$ with an implicit factor 2 for the off-diagonal pairs (the summand is $i
\leftrightarrow j$ symmetric).

## Opposite-spin / same-spin decomposition

Splitting the numerator's two terms gives a natural physical decomposition used throughout qc-rs's API and
output:

$$
E_{\text{OS}} = \sum_{ijab}\frac{(ia|jb)^2}{\Delta_{ijab}}, \qquad
E_{\text{SS}} = \sum_{ijab}\frac{(ia|jb)\bigl[(ia|jb)-(ib|ja)\bigr]}{\Delta_{ijab}}, \qquad
E_{\text{MP2}} = E_{\text{OS}} + E_{\text{SS}},
$$

with $\Delta_{ijab}=\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b$ the denominator. $E_{\text{OS}}$
is the Coulomb-hole (opposite-spin electron pairs avoiding each other) contribution, and $E_{\text{SS}}$ is
the Fermi-hole (same-spin, already partly separated by the Pauli exclusion built into HF) contribution —
empirically, canonical MP2 systematically **overestimates** $E_{\text{SS}}$ relative to the exact
correlation energy and underestimates $E_{\text{OS}}$ less severely, which motivates the two spin-scaled
variants qc-rs implements alongside plain MP2:

$$
E_{\text{SCS-MP2}} = c_{\text{os}}E_{\text{OS}} + c_{\text{ss}}E_{\text{SS}}, \qquad
(c_{\text{os}},c_{\text{ss}}) = \Bigl(1.2,\ \tfrac13\Bigr)\ \text{(Grimme, 2003)},
$$

$$
E_{\text{SOS-MP2}} = c_{\text{os}}'E_{\text{OS}}, \qquad c_{\text{os}}' = 1.3\ \text{(Head-Gordon, 2004;
qc-rs default)},
$$

SOS-MP2 discarding the same-spin term entirely (rather than merely down-weighting it), which — as the next
section shows — is what makes an $O(N^4)$ evaluation possible instead of $O(N^5)$.

## Why RI-only: the density-fitted 3-index factorization

A direct 4-index $(ia|jb)$ evaluation costs $O(N^5)$ in the AO→MO integral transformation alone (four
successive $O(N^4)$-scaling half-transforms) and needs the 4-center ERI tensor materialized or built
on-the-fly per quartet. qc-rs instead factorizes every two-electron integral through the [RI/density-fitting
machinery](density-fitting-ri.md)'s whitened 3-index MO factor,

$$
(ia|jb) = \sum_P B^P_{ia}\,B^P_{jb}, \qquad B^P_{pq} = \sum_Q\bigl(V^{-1/2}\bigr)_{PQ}(Q|pq),
$$

so every occupied-virtual pair integral is a single inner product over the (much smaller) auxiliary index
$P$, and the whole MP2 energy becomes a sequence of `nfit`$\times$`nvir` GEMMs rather than a raw 4-index
transform. This is why `AGENTS.md` states plainly that qc-rs's correlation methods are **RI-C only** — no
4-center MP2 transformation path exists or is planned; you must supply a correlation-fitting auxiliary
basis (the `ric=` argument, distinct from the SCF's own `rijk=` auxiliary basis — the two serve different
accuracy/purpose trade-offs and are typically different published basis sets, e.g. `cc-pvdz-ri/mp2fit`
rather than `cc-pvdz-jkfit`).

:::{prf:algorithm} Canonical RI-MP2 pair loop
:label: alg-rimp2-pairloop

**Input:** whitened occ×vir RI factor $B^P_{ia}$ (built once via [the RI half-transform](density-fitting-ri.md)),
canonical orbital energies $\varepsilon_i,\varepsilon_a$.
**Output:** $E_{\text{OS}}, E_{\text{SS}}$ (hence $E_{\text{MP2}}$, $E_{\text{SCS-MP2}}$, $E_{\text{SOS-MP2}}$).

1. For each occupied pair $(i,j)$ with $j\le i$ (weight 2 if $j\ne i$, else 1):
   1. Form $I_{ab} = (ia|jb) = \sum_P B^P_{ia}B^P_{jb}$ — one GEMM $(n_{\text{vir}}\times n_{\text{fit}})
      \cdot(n_{\text{fit}}\times n_{\text{vir}})$.
   2. $(ib|ja)$ is simply $I_{ba}$ (the same matrix, transposed indices) — no second GEMM needed.
   3. Accumulate $E_{\text{OS}}\mathrel{+}= \text{weight}\sum_{ab}I_{ab}^2/\Delta_{ijab}$ and
      $E_{\text{SS}}\mathrel{+}= \text{weight}\sum_{ab}I_{ab}(I_{ab}-I_{ba})/\Delta_{ijab}$.
2. Return $(E_{\text{OS}}, E_{\text{SS}})$.
:::

The occupied-pair loop is embarrassingly parallel — each pair reads only the shared, immutable $B$ factor
and orbital energies, and writes only its own thread-local $(E_{\text{OS}},E_{\text{SS}})$ accumulator, so
qc-rs threads it directly over `qc_hpc::thread_count()` with one scratch buffer allocated per worker
**outside** the loop (the same no-hot-loop-allocation discipline as everywhere else in the codebase).

## The Laplace transform: SOS-MP2 at $O(N^4)$

SOS-MP2's decision to discard the same-spin term entirely enables an asymptotically cheaper algorithm.
Write the denominator's reciprocal as a **numerical Laplace transform**,

$$
\frac{1}{\Delta_{ijab}} \approx \sum_\alpha w_\alpha\, e^{-\Delta_{ijab}\,t_\alpha}, \qquad
\Delta_{ijab} = (\varepsilon_a-\varepsilon_i)+(\varepsilon_b-\varepsilon_j) > 0,
$$

using a small set of quadrature nodes/weights $(t_\alpha,w_\alpha)$ fit to reproduce $1/\Delta$ to
$<10^{-9}$ relative error over the physical range of orbital-energy gaps (qc-rs tries a $\sim8$–$12$-point
minimax-optimized fit first, falling back to a $\sim40$-point tanh-sinh/double-exponential grid if the
optimized fit can't reach that tolerance — either way accurate to microhartree in the final energy). The
exponential **separates** as a product over the four orbital indices, $e^{-\Delta_{ijab}t}=
e^{-(\varepsilon_a-\varepsilon_i)t/2}e^{-(\varepsilon_a-\varepsilon_i)t/2}e^{-(\varepsilon_b-\varepsilon_j)t/2}
e^{-(\varepsilon_b-\varepsilon_j)t/2}$, letting the sum over $(i,j,a,b)$ **factorize**: defining the
Laplace-scaled factor $\tilde B^P_{ia}(t)=B^P_{ia}\,e^{-(\varepsilon_a-\varepsilon_i)t/2}$,

$$
E_{\text{OS}} = -\sum_\alpha w_\alpha \bigl\|G(t_\alpha)\bigr\|_F^2, \qquad
G^{PQ}(t) = \sum_{ia}\tilde B^P_{ia}(t)\,\tilde B^Q_{ia}(t),
$$

where $G(t)$ is an $n_{\text{fit}}\times n_{\text{fit}}$ matrix built by a single GEMM contraction over the
combined occupied-virtual index — no explicit loop over occupied *pairs* $(i,j)$ at all. Building $G(t)$
costs $O(n_{\text{fit}}^2\,n_{\text{occ}}\,n_{\text{vir}})$ per quadrature point instead of the pair loop's
$O(n_{\text{occ}}^2\,n_{\text{vir}}^2\,n_{\text{fit}})$ — trading the $O(N^5)$-scaling occupied-pair sum for
an $O(N^4)$-scaling sum over a handful ($\sim40$) of quadrature points, each cheap. This is exactly the
`sos_mp2_laplace` implementation, and it is verified to match the direct pair-loop `rimp2_energy(...).os`
to quadrature accuracy — the two algorithms compute the *same* sum, only reorganized. In a distributed
run, only the small $n_{\text{fit}}\times n_{\text{fit}}$ matrix $G(t)$ needs an `Allreduce` per quadrature
point — far cheaper communication than any tensor scaling with $n_{\text{occ}}\cdot n_{\text{vir}}$.

## Open-shell references: UHF and the ROHF ambiguity

**UHF** has no subtlety: each spin channel has its own canonical orbitals and orbital energies, and the
same-spin sum splits into $\alpha\alpha$ and $\beta\beta$ contributions while the opposite-spin sum becomes
purely $\alpha\beta$ (no exchange term crosses spin channels) — `riump2_energy` is a direct generalization
of the closed-shell pair loop with independent $\alpha$/$\beta$ orbital spaces.

**ROHF is genuinely ambiguous** for MP2, and this is worth understanding rather than treating as an
implementation detail: ROHF diagonalizes one **effective Roothaan Fock operator** shared by both spins, so
the eigenvalues of that shared operator are *not* meaningful per-spin MP2 denominators — a well-defined
per-spin Fock matrix $F^\sigma$ exists (it built the ROHF energy), but it is not diagonal in the shared
ROHF MO basis. Two defensible choices exist, and qc-rs implements both (selected by the `lct.rohf_mp2` IOP
key):

| Variant | IOP value | What it does |
|---|---|---|
| **pyscf** (qc-rs default) | `"pyscf"` | Keep the raw ROHF orbitals; take the *diagonal* $\operatorname{diag}(C^{\mathsf T}F^\sigma C)$ in the shared ROHF MO basis as the per-spin orbital energies. Matches PySCF's `mp.MP2(ROHF)` (internally `UMP2` on the ROHF-converted orbitals). |
| **qcrs** | `"qcrs"` | **Semicanonicalize**: diagonalize $F^\sigma$ *within* the occupied block and *within* the virtual block of the shared ROHF orbitals separately (occupied and virtual never mix, so the ROHF density is unchanged); the resulting eigenvalues become the per-spin orbital energies and the orbitals themselves rotate within each block. |

Both are legitimate — they differ only in how the un-diagonalized off-diagonal $F^\sigma_{ij}$/$F^\sigma_{ab}$
blocks are handled (dropped in `pyscf`, folded into a rotation in `qcrs`) — and the resulting correlation
energies are close but not identical, which the verified example below demonstrates directly. Frozen-core
truncation (`lct.frozen_core`, an integer count of lowest occupied orbitals excluded from the correlation
sum — default `0`, i.e. all-electron) applies identically to both variants and to UHF/RHF, simply by
dropping the lowest `n_frozen` rows of the occupied block before building $B^P_{ia}$.

Verified examples — water/cc-pVDZ (with `cc-pvdz-ri/mp2fit` as the correlation-fitting auxiliary basis):

```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="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="mp2").run()
m.scf.energy, m.lct.energy, m.lct.e_corr, m.lct.e_os, m.lct.e_ss
# (-76.02679364497408, -76.23074681436064, -0.20395316938655841, -0.15237768872141289, -0.05157548066514553)

m_scs = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="scs-mp2").run()
m_scs.lct.e_corr   # -0.2000450533540773 == 1.2*(-0.15237768872141289) + (1/3)*(-0.05157548066514553)

m_sos = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="sos-mp2").run()
m_sos.lct.e_corr   # -0.19809099564379767 == 1.3*(-0.15237768872141289)

# frozen-core: freeze the O 1s core orbital
m_fc = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom",
                   iop={"lct.frozen_core": 1}).scf(ref="r").lct(method="mp2").run()
m_fc.lct.e_corr   # -0.20161465527331127 -- smaller in magnitude: the frozen core orbital
                  # can no longer correlate into virtuals

# ROHF-MP2: CH3 radical, both variants
ch3 = "C 0 0 0; H 0 1.079 0; H 0.934 -0.540 0; H -0.934 -0.540 0"
for variant in ("pyscf", "qcrs"):
    r = qc.chk.new(atom=ch3, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom", spin=2,
                    iop={"lct.rohf_mp2": variant}).scf(ref="ro").lct(method="mp2").run()
    print(variant, r.scf.energy, r.lct.e_corr)
# pyscf -39.55963721124465 -0.13088782816695832
# qcrs  -39.55963721124465 -0.13083446557503792
```

The `pyscf`/`qcrs` ROHF-MP2 correlation energies agree to about $5\times10^{-5}\ E_h$ here — close, as
expected for two defensible resolutions of the same ambiguity, but not identical, which is exactly why the
choice is an explicit, documented IOP key rather than a silently-picked internal detail.

## Where correlation stops today

`lct(method=...)` is qc-rs's umbrella for every SCF-reference and CASSCF-reference correlation method.
Today only the RI-MP2 family (`mp2`, `scs-mp2`, `sos-mp2`) is a real, energy-producing implementation;
`cc2` (an SCF-reference coupled-cluster-like method) and the CASSCF-reference methods `caspt2`/`nevpt2`
remain **mock** — they accept the same pending-step syntax and produce a deterministic placeholder energy,
but not a physically meaningful one yet. [The multireference outlook chapter](multireference-outlook.md)
covers what a real CASPT2/NEVPT2 implementation would need and why CASSCF's own mock status blocks it.

:::{exercise}
:label: ex-post-hf-theory

1. SOS-MP2 discards the same-spin term entirely rather than merely reweighting it like SCS-MP2. Using the
   Laplace-transform derivation, explain in one or two sentences why discarding $E_{\text{SS}}$ specifically
   (not $E_{\text{OS}}$) is what enables the exponential-factorization trick and the $O(N^5)\to O(N^4)$
   speedup.
2. The verified frozen-core example shows $|E_{\text{corr}}|$ *decreasing* in magnitude when the oxygen 1s
   orbital is frozen ($-0.20395$ vs $-0.20161\ E_h$). Explain why removing an occupied orbital from the
   correlation sum can only ever decrease $|E_{\text{corr}}|$ (for a fixed sign of each term), never
   increase it.
3. The `pyscf` and `qcrs` ROHF-MP2 variants agree on the *SCF* energy but differ slightly on the MP2
   correlation energy. Why does the ambiguity affect only the correlation step and not the mean-field step
   that precedes it?
:::

:::{solution} ex-post-hf-theory
:class: dropdown

1. The same-spin term $E_{\text{SS}}$ has an *exchange-like* structure with an antisymmetrized numerator,
   $(ia|jb)-(ib|ja)$ — the second piece couples $a$ and $b$ *across* the two occupied indices $i,j$ in a
   way that does not separate into a product of independent per-orbital exponential factors under the
   Laplace substitution. The opposite-spin term's numerator $(ia|jb)^2$, in contrast, is already a plain
   product of two identical-structure integrals, so the exponential denominator factorizes cleanly into
   $i,a$ and $j,b$ pieces that recombine into the $G(t)=\sum_{ia}\tilde B_{ia}\tilde B_{ia}$ contraction.
   Discarding $E_{\text{SS}}$ is precisely what removes the non-factorizable term from the sum.
2. Each term in the MP2 double sum, $(ia|jb)[2(ia|jb)-(ib|ja)]/\Delta_{ijab}$ for a fixed occupied-virtual
   quadruple, contributes independently and (for the physically normal case of a negative-definite MP2
   energy expression) each occupied index $i$ only ever *adds* excitation channels into virtuals. Freezing
   an occupied orbital removes every term in the sum that involves that orbital as $i$ or $j$ — it can only
   subtract already-non-negative-magnitude contributions from the total, never add new ones, so
   $|E_{\text{corr}}|$ can only shrink (or stay the same, if the frozen orbital happened to contribute
   exactly zero) when orbitals are removed from the active correlation space.
3. The ROHF self-consistent field is entirely determined by the **shared** Roothaan effective Fock
   operator and its variational minimization — there is only one well-defined ROHF energy and one set of
   ROHF orbitals/occupations, so both variants start from *literally the same converged mean-field state*
   (confirmed by the identical `scf.energy` above). The ambiguity is not in what ROHF *is*; it is in how to
   assign meaningful **per-spin single-particle energies** to a state whose defining operator is spin-shared
   — a question that only arises once you try to build a spin-resolved *post*-HF perturbative correction,
   which is exactly the MP2 step, not the SCF step.
:::

The RI factor $B^P_{ia}$ this chapter builds on is derived in full in [the density-fitting/RI
chapter](density-fitting-ri.md); [the multireference outlook](multireference-outlook.md) picks up where
this chapter's "mock" methods leave off.