Post-SCF methods: RI-MP2

Post-SCF methods: RI-MP2#

Hartree–Fock is a mean-field method: each electron feels only the average field of the others, so it misses the way electrons instantaneously avoid one another. The energy of that missing effect is the correlation energy, and recovering it is the job of post-SCF (post-Hartree–Fock) methods. This chapter covers the family qc-rs implements today: RI-MP2 and its spin-scaled variants.

Theory: correlation and second-order perturbation theory#

Recall the definition from Part II:

\[ E_{\text{corr}} = E_{\text{exact}} - E_{\text{HF}} \;<\; 0 . \]

It is a small fraction of the total energy but chemically decisive — bond energies, reaction barriers, and dispersion all live there. The cheapest systematic way to estimate it is second-order Møller–Plesset perturbation theory (MP2): treat the difference between the true electron–electron repulsion and the HF mean field as a perturbation, and take the second-order energy correction. For a closed shell it is

\[ E^{(2)} = \sum_{ij}^{\text{occ}} \sum_{ab}^{\text{virt}} \frac{ (ia\,|\,jb)\,\big[\,2(ia\,|\,jb) - (ib\,|\,ja)\,\big] } { \varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b } , \]

a sum over occupied (\(i,j\)) and virtual (\(a,b\)) orbitals of the two-electron integrals \((ia|jb)\) divided by orbital-energy denominators. The two terms in the bracket are the opposite-spin and same-spin contributions — qc-rs reports them separately (e_os, e_ss), which is what makes the spin-scaled variants below possible.

The RI (“density-fitting”) approximation#

The bottleneck is the four-index integral \((ia|jb)\). The resolution-of-identity (RI), or density-fitting, approximation fits the AO pair density into an auxiliary basis \(\{\chi_P\}\) under the Coulomb metric,

\[ (\mu\nu\,|\,\lambda\sigma) \;\approx\; \sum_{PQ} (\mu\nu\,|\,P)\,\big(V^{-1}\big)_{PQ}\,(Q\,|\,\lambda\sigma), \qquad V_{PQ} = (P|Q), \]

collapsing the four-index object into a product of three-index ones through a single whitened factor

\[ \boxed{\,B_{P,\mu\nu} = \sum_Q \big(V^{-1/2}\big)_{PQ}\,(Q\,|\,\mu\nu)\,} \qquad\Longrightarrow\qquad (\mu\nu\,|\,\lambda\sigma) \approx \sum_P B_{P,\mu\nu}\,B_{P,\lambda\sigma}. \]

\(V^{-1/2}\) comes from eigendecomposing the symmetric metric \(V=U\Lambda U^\top\) and dropping near-linearly- dependent modes below a cutoff (ri.metric_cutoff, default \(10^{-10}\)) — the only approximation RI makes, controlled entirely by the auxiliary basis and that cutoff, identically for every qc-rs backend. Crucially, \(B\) depends only on integrals, never on the density, so it is built once and reused; MP2 needs it in the molecular-orbital form \((ia|jb)\approx\sum_P B_{P,ia}B_{P,jb}\) (the same \(B\) that RI-JK uses for the SCF, half-transformed to occupied×virtual indices), which is why RI is a single unified machinery (full derivation) shared by the SCF’s eri="ri-*" and MP2’s ric=. This is the only MP2 path in qc-rs — fast and accurate, the fitting error well below chemical accuracy — which is why you supply a correlation-fitting auxiliary basis with ric= on qc.chk.new(...).

Usage#

MP2 is an LCT method (lct, for the correlation methods that sit on top of an SCF reference). Run the SCF first, then add the lct step:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

hf = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom",
                ric="cc-pvdz-ri/mp2fit").scf(ref="r").run()   # note the ric= aux basis
m  = hf.lct(method="mp2").run()

m.lct.energy    # -76.230747   total MP2 energy (E_HF + E_corr)
m.lct.e_corr    # -0.203953    the correlation energy
m.lct.e_os      # opposite-spin part
m.lct.e_ss      # same-spin part

The ric= auxiliary basis (here cc-pvdz-ri/mp2fit, the correlation-fitting set matched to cc-pvdz) is required — RI-MP2 has no non-RI fallback. Choose the aux to match the orbital basis (cc-pvtz → cc-pvtz-ri/mp2fit, def2 sets → their def2-*/mp2fit).

Tip

How good is RI-MP2? The RI fitting error is tiny. For this water/cc-pVDZ case, qc’s RI-MP2 correlation energy (−0.203953) agrees with a conventional (non-RI) MP2 to 1.5 × 10⁻⁵ Ha — far below chemical accuracy (~1.6 mHa).

Spin-component-scaled variants#

Because MP2 slightly overcounts correlation, two empirical rescalings of the opposite-/same-spin parts often improve accuracy at no extra cost — just change method:

scs = hf.lct(method="scs-mp2").run().lct.energy    # -76.226839   SCS-MP2 (Grimme)
sos = hf.lct(method="sos-mp2").run().lct.energy    # -76.224885   SOS-MP2 (same-spin dropped, O(N⁴))

  • scs-mp2 scales opposite- and same-spin by different factors (⅚ and ⅓).

  • sos-mp2 keeps only the opposite-spin term (scaled 1.3), which enables an O(N⁴) Laplace algorithm.

Open shells and frozen core#

RI-MP2 works on RHF, UHF, and ROHF references (an open-shell radical uses its UHF/ROHF orbitals). The frozen-core approximation (excluding chemically inert core orbitals from the correlation sum, the usual default in the field) and other tuning are iop keys under lct.* — e.g. iop={"lct.frozen_core": True} — documented in the IOP reference.

Note

What is implemented today The RI-MP2 family (mp2 / scs-mp2 / sos-mp2) computes real energies. The other lct methods (cc2) and the multireference methods (casscf, caspt2, nevpt2, fci, dmrg) and excited-state td are recognized but mock (they set up the workflow without a real correlation energy yet). This chapter covers the methods that produce real numbers.

Worked example: the correlation energy#

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

hf = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom",
                ric="cc-pvdz-ri/mp2fit").scf(ref="r").run()
print(f"E(RHF)     = {hf.scf.energy:.6f}")            # -76.026794

for method in ("mp2", "scs-mp2", "sos-mp2"):
    e = hf.lct(method=method).run().lct.energy
    print(f"E({method:8}) = {e:.6f}")
# E(mp2     ) = -76.230747   (E_corr = -0.203953)
# E(scs-mp2 ) = -76.226839
# E(sos-mp2 ) = -76.224885

The RHF energy drops by ~0.20 Ha once correlation is added — small next to the total, but this is the part that gets bond energies and reaction thermochemistry right.

Exercise 5

  1. You call hf.lct(method="mp2").run() but hf was built without ric=. What goes wrong, and why does MP2 need it when the SCF did not?

  2. From m.lct.e_os and m.lct.e_ss, how would you reconstruct the plain MP2 correlation energy? And what does SOS-MP2 do to e_ss?

RI-MP2 gave us a correlated energy at one geometry. To find equilibrium structures and forces, the next chapter turns to analytic gradients and geometry optimization.