Post-Hartree-Fock correlation theory

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

Algorithm 3 (Canonical RI-MP2 pair loop)

Input: whitened occ×vir RI factor \(B^P_{ia}\) (built once via the RI half-transform), 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):

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 covers what a real CASPT2/NEVPT2 implementation would need and why CASSCF’s own mock status blocks it.

Exercise 8

  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?

The RI factor \(B^P_{ia}\) this chapter builds on is derived in full in the density-fitting/RI chapter; the multireference outlook picks up where this chapter’s “mock” methods leave off.