MP2 and double hybrids¶

vibe-qc ships second-order Møller-Plesset perturbation theory (MP2) as a post-SCF correlation method for molecules, plus the spin-component-scaled variants (SCS-MP2, SOS-MP2) and the first double-hybrid density functional (B2PLYP). All variants compose with density fitting (RI-MP2) and with the unrestricted (UHF) reference.

At a glance¶

Variant	Entry point	Reference	Recommended default
Canonical MP2	`run_mp2(mol, basis, rhf)`	Szabo-Ostlund eq. 6.71	Parity validation only
RI-MP2	`run_mp2(..., MP2Options(density_fit=True, aux_basis="cc-pvtz-ri"))`	Vahtras-Almlöf-Feyereisen 1993	Production
SCS-MP2	`run_scs_mp2(mol, basis, rhf)`	Grimme 2003 (c_os=6/5, c_ss=1/3)	Thermochemistry, reactions
SOS-MP2	`run_sos_mp2(mol, basis, rhf)`	Jung-Lochan-Dutoi-Head-Gordon 2004 (c_os=1.3, c_ss=0)	Even cheaper than SCS
Canonical / RI / SCS / SOS UMP2	`run_ump2`, `run_scs_ump2`, `run_sos_ump2`	UHF reference	Open-shell
B2PLYP	`run_b2plyp(mol, basis)`	Grimme 2006	First double hybrid

The four run_* convenience wrappers (run_scs_mp2, run_sos_mp2, run_scs_ump2, run_sos_ump2) default to RI-MP2 with the per-zeta RIfit aux auto-resolved from the orbital basis name. Canonical MP2 is reachable only via the explicit density_fit=False flag — the brief’s recommended setup throughout, since canonical MP2 scales as O(N⁵) in the AO→MO ERI transform while RI-MP2 scales as O(N⁴) with sub-µHa drift on the standard bases.

Canonical MP2 + RI-MP2¶

from vibeqc import Atom, BasisSet, Molecule, RHFOptions, MP2Options
from vibeqc import run_rhf, run_mp2, default_aux_basis_for

mol = Molecule.from_xyz("h2o.xyz")
basis = BasisSet(mol, "cc-pvtz")

# RHF reference. MP2 requires a converged RHF result.
rhf_opts = RHFOptions(); rhf_opts.conv_tol_energy = 1e-10
hf = run_rhf(mol, basis, rhf_opts)

# Canonical MP2 (O(N⁵)) — useful for parity validation.
mp2 = run_mp2(mol, basis, hf)
print(mp2.e_total, mp2.e_correlation)

# RI-MP2 (O(N⁴), recommended for production).
mp2_opts = MP2Options()
mp2_opts.density_fit = True
mp2_opts.aux_basis = default_aux_basis_for(basis.name, kind="ri")
mp2_df = run_mp2(mol, basis, hf, mp2_opts)

MP2Result carries the per-component Grimme decomposition (see next section). Sum invariant: e_os + e_ss == e_correlation at canonical c_os = c_ss = 1.

Spin-component-scaled MP2¶

The MP2 correlation energy splits into an opposite-spin (αβ singlet pair) piece and a same-spin (αα + ββ triplet pair) piece. The Grimme decomposition (Grimme, J. Chem. Phys. 118, 9095 (2003)) sets:

e_os = Σ (ia|jb)² / Δ                          (singlet pair)
e_ss = Σ [(ia|jb)² − (ia|jb)(ib|ja)] / Δ       (triplet pair)
e_correlation(c_os, c_ss) = c_os · e_os + c_ss · e_ss

with Δ = ε_i + ε_j − ε_a − ε_b. Setting different (c_os, c_ss) gives the named recipes:

Recipe	c_os	c_ss	Origin
Canonical MP2	1	1	Szabo-Ostlund eq. 6.71
SCS-MP2	6/5	1/3	Grimme 2003 (semi-empirical fit)
SOS-MP2	1.3	0	Jung et al. 2004 (drop same-spin)
B2PLYP correction	0.27	0.27	Grimme 2006 (double hybrid)
DSD double hybrids	varies	varies	Kozuch-Martin 2011

vibe-qc’s MP2Options exposes c_os and c_ss (defaults 1.0, 1.0). Both compose with density_fit=True without code duplication — the same kernel computes (e_os, e_ss); the scales apply at the final sum.

from vibeqc import run_scs_mp2, run_sos_mp2

scs = run_scs_mp2(mol, basis, hf)            # RI by default, c_os=6/5, c_ss=1/3
sos = run_sos_mp2(mol, basis, hf)            # RI by default, c_os=1.3, c_ss=0

# Custom scaling factors:
mine = run_scs_mp2(mol, basis, hf, c_os=1.10, c_ss=0.50)

The unrestricted variants (run_scs_ump2, run_sos_ump2) apply the same coefficients to UMP2’s channel decomposition: c_os scales the αβ channel (UMP2Result.e_ab), c_ss scales the αα + ββ sum (UMP2Result.e_aa + UMP2Result.e_bb).

When to use which¶

SCS-MP2 is a well-validated improvement over MP2 across thermochemistry, reaction barriers, conformer ordering. Default if you don’t have a specific reason to pick something else.
SOS-MP2 drops the more expensive same-spin term entirely; with Laplace-transformed denominators it is O(N⁴) without RI. vibe-qc’s vanilla implementation is O(N⁴) (via RI) for both; SOS still outperforms SCS on cost-per-accuracy in some regimes (Jung et al.).
Canonical MP2 is for parity validation against reference codes (PySCF.MP2, ORCA ! MP2) — every published MP2 number used the unscaled formula.

Verifying the RI fit: `report_ri_residual`¶

RI-MP2 replaces the canonical four-index (ia|jb) with the Coulomb-metric Dunlap fit

\[(ia|jb)_{\mathrm{RI}} = \sum_{PQ} (ia|P)\, [V^{-1}]_{PQ}\, (Q|jb).\]

For any finite auxiliary basis this introduces a per-bucket residual. The residual is structurally biased: \(E_{\mathrm{os}}\) is over-estimated (less negative than canonical) and \(E_{\mathrm{ss}}\) is under-estimated (more negative than canonical), and the two errors partially cancel in the unscaled sum but survive SCS / SOS scaling. Concretely, on H₂O / cc-pVTZ:

Aux basis	\(n_{\mathrm{aux}}\)	\(\Delta E_{\mathrm{os}}\)	\(\Delta E_{\mathrm{ss}}\)	\(\Delta E_{\mathrm{SCS}}\)
cc-pvtz-ri	141	+5.05e-5	−2.43e-5	+5.25e-5
cc-pvqz-ri	242	+2.09e-5	−4.56e-6	+2.35e-5
def2-qzvpp-rifit	253	+1.56e-5	−5.42e-6	+1.69e-5
cc-pvtz-jkfit (mis-pair)	139	+3.33e-4	−4.21e-4	+2.60e-4

(Bucket deltas in Hartree; SCS uses Grimme’s \(c_{\mathrm{os}}=6/5\), \(c_{\mathrm{ss}}=1/3\).)

MP2Options.report_ri_residual = True turns on an opt-in diagnostic that, alongside the RI build, also computes the canonical (ia|jb) and returns the per-bucket residuals on MP2Result:

opts = MP2Options()
opts.density_fit = True
opts.aux_basis = "cc-pvtz-ri"
opts.report_ri_residual = True              # opt-in

mp2 = run_mp2(mol, basis, hf, opts)

assert mp2.ri_residual_reported             # True iff the flag fired
print(mp2.e_os_ri_residual)                 # +5.048e-05 on H2O/cc-pvtz
print(mp2.e_ss_ri_residual)                 # -2.430e-05

The diagnostic doubles the cost of the RI-MP2 call (it forces the O(N⁵) canonical AO→MO transform). It is intended for verification of aux-basis adequacy in benchmarks and unit tests — not production. The flag is silently no-op when density_fit = False (ri_residual_reported stays False).

The four SCS/SOS convenience wrappers accept the flag directly, so you don’t have to hand-build an options struct:

scs = run_scs_mp2(mol, basis, hf, report_ri_residual=True)
print(scs.e_os_ri_residual, scs.e_ss_ri_residual)

run_sos_mp2, run_scs_ump2, and run_sos_ump2 take the same report_ri_residual=False keyword.

When the diagnostic earns its keep¶

Aux-basis sizing for a new system class. Run once on a small representative geometry; if |e_os_ri_residual| exceeds your tolerance (rule of thumb: ≤10 µHa for SCS work), bump to the next zeta of the matching RIfit family.
Catching a wrong aux family. Pairing a JKfit aux with MP2 is technically valid but over-fit — the bucket residuals are typically 10× larger than with the matching RIfit family. The diagnostic surfaces this immediately.
Cross-code parity setup. Use it to confirm both codes are running RI (not canonical) before chasing sub-µHa OS/SS differences — if both codes report similar OS+/SS− residuals against their own canonicals, the comparison is well-framed.

What the diagnostic does not tell you¶

The residuals are vs. canonical at the same RHF reference, not vs. the basis-set limit. A canonical-MP2 number itself carries basis-set-incompleteness error of ~mHa on cc-pVTZ; the RI residual is strictly the fit-on-top-of-that. For total-energy accuracy budgets, combine the RI residual with a basis-set-extrapolation estimate.

Open-shell: UMP2 and its scaled variants¶

run_ump2(mol, basis, uhf, UMP2Options(...)) takes a converged UHF reference. UMP2Result.e_aa / e_bb / e_ab are the three spin-channel energies (already in the natural Grimme convention).

from vibeqc import run_uhf, run_scs_ump2, UHFOptions

mol = Molecule([Atom(8,[0,0,0]),Atom(1,[0,0,1.83])], multiplicity=2)
basis = BasisSet(mol, "cc-pvdz")
uhf_opts = UHFOptions(); uhf_opts.conv_tol_energy = 1e-10
uhf = run_uhf(mol, basis, uhf_opts)

scs_uhf = run_scs_ump2(mol, basis, uhf)
# scs_uhf.e_correlation == 6/5 · e_ab + 1/3 · (e_aa + e_bb)

The RI fit-residual diagnostic works identically for UMP2: UMP2Options.report_ri_residual = True (alongside density_fit = True) populates per-channel residuals on the result — UMP2Result.e_aa_ri_residual, e_bb_ri_residual, e_ab_ri_residual, plus the ri_residual_reported flag. The αβ channel carries the opposite-spin-like bias and the αα / ββ channels the same-spin-like bias, mirroring the closed-shell signature documented above. Same opt-in cost caveat: the diagnostic doubles the runtime of the RI-UMP2 call.

opts = UMP2Options()
opts.density_fit = True
opts.aux_basis = "def2-svp-rifit"
opts.report_ri_residual = True

ump2 = run_ump2(mol, basis, uhf, opts)
assert ump2.ri_residual_reported
print(ump2.e_aa_ri_residual, ump2.e_bb_ri_residual, ump2.e_ab_ri_residual)

Double hybrids¶

vibe-qc ships two double hybrids out of the box:

Functional	SCF mix	MP2 mix (c_os, c_ss)	Reference
B2PLYP	0.53·HF + 0.47·B88, 0.73·LYP	0.27, 0.27	Grimme 2006
DSD-PBEP86	0.69·HF + 0.31·PBE, 0.44·P86	0.55, 0.09	Kozuch-Martin 2011

Both are dispatched via name-specific wrappers (run_b2plyp, run_dsd_pbep86) or the generic run_double_hybrid(mol, basis, "name", ...). All three return a :class:DoubleHybridResult with the same shape.

B2PLYP¶

A double hybrid combines a hybrid-DFT SCF step with a scaled MP2 correlation correction on the converged KS orbitals:

E_DH = E_RKS[xc_scf] + c_os · E_os + c_ss · E_ss

B2PLYP (Grimme, J. Chem. Phys. 124, 034108 (2006)) is the first double hybrid wired up in vibe-qc:

E_B2PLYP = E_RKS[0.53·HF + 0.47·B88, 0.73·LYP] + 0.27 · E_MP2_corr

The orchestrator does both steps in one call:

from vibeqc import run_b2plyp

result = run_b2plyp(mol, basis)
print(result)
# DoubleHybridResult(functional='b2plyp', e_total=..., e_rks=..., e_mp2_corr=...)

print(result.e_total)         # the published double-hybrid energy
print(result.rks.energy)      # the hybrid SCF step alone
print(result.mp2.e_correlation)  # 0.27 · (e_os + e_ss)

run_b2plyp defaults to DF on both steps: RI-J + RI-K for the hybrid SCF, RI-MP2 for the correction. Aux bases are auto-resolved from the orbital basis name (JKfit for SCF, RIfit for MP2). Pass density_fit_mp2=False for an explicit canonical-MP2 parity validation run.

Anatomy¶

`Functional("b2plyp")` accessor	Value
`hf_exchange_fraction`	0.53
`is_hybrid`	`True`
`is_double_hybrid`	`True`
`mp2_c_os`	0.27
`mp2_c_ss`	0.27
`kind`	`XCKind.GGA`

Non-double-hybrid functionals (LDA, GGA, hybrid-GGA: PBE0, B3LYP, PW1PW, …) carry is_double_hybrid = False and mp2_c_os = mp2_c_ss = 0. The same Functional accessors generalise the dispatcher pattern to DSD-PBEP86 / PWPB95 / ωB97M(2) once each is registered.

DSD-PBEP86¶

The second double hybrid, Kozuch-Martin 2011 (PCCP 13, 20104):

E_DSD-PBEP86 = E_RKS[0.69·HF + 0.31·PBE, 0.44·P86]
             + 0.55 · E_os + 0.09 · E_ss

from vibeqc import run_dsd_pbep86, Functional

result = run_dsd_pbep86(mol, basis)
print(result.e_total)

fn = Functional("dsd-pbep86")  # or "dsdpbep86"
print(fn.hf_exchange_fraction)   # 0.69
print(fn.mp2_c_os, fn.mp2_c_ss)  # (0.55, 0.09)  — asymmetric scaling

DSD-PBEP86 is the prototypical asymmetric-c_os/c_ss double hybrid — c_os ≠ c_ss is what distinguishes the DSD family from B2PLYP-style symmetric mixing. vibe-qc’s MP2Options.c_os / c_ss cover this case natively; the Functional resolver carries the published coefficients.

The published DSD-PBEP86 method usually pairs with the D3(BJ) dispersion correction; the dispatcher returns the un-dispersed XC + MP2 total. Add D3(BJ) post-hoc via vibeqc.compute_d3bj when you need the published “with-D” energy. A dedicated D4 parameter set for DSD-PBEP86 is a separate upcoming item.

Dispersion: D3(BJ) and D4¶

Double hybrids inherit the dispersion question from regular hybrid DFT: the published B2PLYP / DSD-PBEP86 totals in the literature almost always include a D3(BJ) or D4 dispersion correction on top of the XC + MP2 piece. The dispatcher folds it in when you ask for it:

# Un-dispersed B2PLYP (the "B2PLYP / no-D" total):
no_d = run_b2plyp(mol, basis)

# Full B2PLYP-D4 (the published total — Caldeweyher-Bannwarth-Grimme 2019):
d4   = run_b2plyp(mol, basis, dispersion="d4")
print(d4.e_total)         # XC + MP2 + D4
print(d4.dispersion.energy)  # the D4 piece alone (Hartree)

dispersion="d4" reaches the reference Caldeweyher-Bannwarth-Grimme 2019 implementation via the optional dftd4 Python package (install with pip install -e '.[dispersion]'). When dftd4 is not installed the dispatcher raises with an actionable hint; dispersion=None (the default) returns the un-dispersed total and never touches the dispersion path.

For arbitrary (non-double-hybrid) functionals, call vibeqc.compute_d4(mol, functional) directly and add the energy to your SCF total. Any functional name dftd4 catalogs is accepted (B3LYP, PBE0, PW1PW, r²SCAN, …).

D3(BJ) is also available (vibe-qc’s own framework + optional dftd3 backend) — see :func:vibeqc.compute_d3bj and the Density fitting page for the dispersion section. Direct dispersion="d3bj" in the double-hybrid dispatcher is a future addition; for now call compute_d3bj after run_b2plyp(...) and add the energies.

Adding more double hybrids¶

The same Functional + dispatcher pattern generalises to other double-hybrid families. To register a new one:

Add an alias entry to cpp/src/xc.cpp::resolve_alias with the SCF mix (HF fraction + weighted libxc components) and the MP2 c_os / c_ss pair.
Add a thin Python wrapper run_<name> in python/vibeqc/__init__.py that delegates to run_double_hybrid(mol, basis, "<name>", ...).
Add a parity test against the reference code’s same-recipe composition.

Partially unblocked: PWPB95 (Goerigk-Grimme 2011) is a meta-GGA double hybrid — vibe-qc’s molecular meta-GGA τ-density path is now live (tpssh, m06-2x, …; see functionals.md § Meta-GGA functionals), so a pwpb95 registration in resolve_alias() should now compile and run. Still missing: the dispersion piece (PWPB95-D3(BJ) parameters in dispersion_params.cpp) and a parity test against the published reference. ωB97M(2) (Mardirossian-Head-Gordon 2018) is a range-separated meta-GGA double hybrid — the τ side is in, but it still needs the range-separated Coulomb attenuation machinery (orthogonal workstream).

Cross-code parity¶

tests/test_b2plyp.py and tests/test_dsd_pbep86.py pin vibe-qc’s dispatchers against PySCF’s hand-composed equivalents (matching dft.RKS(xc='...') for the SCF + mp.MP2(mf).density_fit().kernel() correlation correction scaled by the published coefficients):

Method	H₂O/cc-pVDZ gap	Tolerance
B2PLYP	1.93e-7 Ha	1e-5 Ha (grid-coupling slack)
DSD-PBEP86	1.07e-5 Ha	5e-5 Ha (P86 is more grid-sensitive)

The gap is set by the DFT-grid difference between vibe-qc and PySCF; total energies at our own grid are reproducible to machine precision.

For a full multi-system, multi-variant validation harness against ORCA, see examples/molecular/mp2_benchmarks/ — S22 (closed shell) + open-shell complement, same-converger discipline, automated comparator.

Citations¶

For published work that uses these methods:

Canonical MP2 — Møller, Plesset, Phys. Rev. 46, 618 (1934); operational form in Szabo-Ostlund, Modern Quantum Chemistry (eq. 6.71).
RI-MP2 — Vahtras, Almlöf, Feyereisen, Chem. Phys. Lett. 213, 514 (1993); Weigend, Häser, Patzelt, Ahlrichs, Chem. Phys. Lett. 294, 143 (1998) (the RIfit auxiliary basis family vibe-qc uses).
SCS-MP2 — Grimme, J. Chem. Phys. 118, 9095 (2003).
SOS-MP2 — Jung, Lochan, Dutoi, Head-Gordon, J. Chem. Phys. 121, 9793 (2004).
B2PLYP — Grimme, J. Chem. Phys. 124, 034108 (2006).
DSD-PBEP86 — Kozuch, Martin, Phys. Chem. Chem. Phys. 13, 20104 (2011).
D4 dispersion — Caldeweyher, Bannwarth, Grimme, J. Chem. Phys. 150, 154122 (2019); reference implementation: Grimme group dftd4 Python package.