Verifying an RI-MP2 auxiliary basis¶

When you reach for RI-MP2 in production (a production route, label vocabulary in Methods: theory and status), you replace the four-index AO integrals by three-index integrals over an auxiliary basis. Both the canonical and the RI method scale formally as O(N⁵), but the fitted integrals cut the prefactor and the storage by an order of magnitude, and you inherit the fit’s residual on top of the basis-set-incompleteness error you’d have had anyway. For most workflows the per-zeta RIfit family (cc-pvNz-ri, def2-Nzv-rifit) is tight enough that the fit residual stays well below ~µHa per atom. But when you move to a new system class, switch zeta, or accidentally pair an orbital basis with the wrong aux family, the residual can balloon by an order of magnitude, and the silent failure mode is a wrong SCS / SOS / double-hybrid number, not a crash.

This tutorial walks through MP2Options.report_ri_residual, the opt-in diagnostic added in #mp2-ri-residual, on water with two aux choices. It complements the double-hybrid B2PLYP tutorial (which is where you’d actually consume this diagnostic in anger, the post-SCF MP2 step in a B2PLYP run inherits the same RI fit error).

Theory¶

The Coulomb-metric fit. MP2 needs only the integrals \((ia|jb)\) over the pair densities \(\rho_{ia}(\mathbf r) = \phi_i(\mathbf r)\phi_a(\mathbf r)\) of an occupied and a virtual orbital. RI-MP2 (Feyereisen, Martin, Fitzgerald, George, and Komornicki, Andrew (1993). *Use of approximate integrals in ab initio theory. An application in MP2 energy calculations.* Chemical Physics Letters, **208**(5), 359–363. doi:[10.1016/0009-2614(93)87156-W](https://doi.org/10.1016/0009-2614(93)87156-W)) expands each pair density in an auxiliary basis \(\{\chi_P\}\), \(\tilde\rho_{ia} = \sum_P d^{ia}_P \chi_P\), with coefficients that minimise the Coulomb self-repulsion of the residual \(\delta_{ia} = \rho_{ia} - \tilde\rho_{ia}\). That gives the familiar three-index form

\[ (ia|jb)_{\mathrm{RI}} = \sum_{PQ} (ia|P)\, [V^{-1}]_{PQ}\, (Q|jb), \quad V_{PQ} = (P|Q), \qquad (ia|jb) = (ia|jb)_{\mathrm{RI}} + (\delta_{ia}|\delta_{jb}) , \]

where the second identity is exact for the Coulomb metric: the residual is Coulomb-orthogonal to the fit space, so the cross terms vanish and the error of every integral is a Gram matrix of residuals, positive semidefinite as a matrix over \((ia),(jb)\). The fit error is therefore second order in the size of the residual, which is why a moderately sized per-zeta auxiliary family (Weigend, Florian, Köhn, Andreas, and Hättig, Christof (2002). *Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations.* Journal of Chemical Physics, **116**(8), 3175–3183. doi:[10.1063/1.1445115](https://doi.org/10.1063/1.1445115): cc-pVXZ-RI, def2-XZVP-RIFIT) is enough for microhartree-per-atom accuracy, and why the residual shrinks when the auxiliary basis grows. The MP2 energy in the Grimme decomposition is

\[ E_{\mathrm{os}} = -\sum_{ijab}\frac{(ia|jb)^2}{\varepsilon_a+\varepsilon_b-\varepsilon_i-\varepsilon_j}, \qquad E_{\mathrm{ss}} = -\sum_{ijab}\frac{(ia|jb)\,[(ia|jb)-(ib|ja)]}{\varepsilon_a+\varepsilon_b-\varepsilon_i-\varepsilon_j}, \]

with \(\varepsilon\) the orbital energies, and the scaled variants use \(c_\mathrm{os}E_\mathrm{os} + c_\mathrm{ss}E_\mathrm{ss}\). The diagnostic in this tutorial is the difference of each bucket against the canonical value, \(\Delta E_x = E_x^{\mathrm{RI}} - E_x^{\mathrm{canonical}}\), computed by running the canonical AO to MO transformation alongside the RI build.

The structural fingerprint. The residual has a definite sign pattern in the Grimme decomposition: \(E_{\mathrm{os}}\) is over-estimated (less negative than canonical), and \(E_{\mathrm{ss}}\) is under-estimated (more negative than canonical). The two errors partially cancel in the unscaled \(E_{\mathrm{corr}} = E_{\mathrm{os}} + E_{\mathrm{ss}}\) but survive SCS / SOS scaling, because SCS weights the OS bucket (c_os = 6/5) more than three times as heavily as the SS bucket (c_ss = 1/3).

A wrong-direction signature in the diagnostic (ΔE_os < 0, ΔE_ss > 0, or both with the same sign) almost certainly indicates a numerical bug in the fit pipeline, not a fit-basis issue. The diagnostic is exactly the right tool to catch that quickly.

Where it helps and where it fails. RI-MP2 is the right default for anything beyond the smallest systems, and the residual is a property of the orbital-basis and auxiliary-basis pairing, not of the molecule, so one check per (basis family, element set) is enough. It fails, silently, when the auxiliary basis does not match the orbital family or lacks coverage for an element (the wrong number is then a plausible-looking energy, not a crash), and it inherits every limitation of MP2 itself (near-degenerate and stretched systems). The diagnostic detects the first kind of failure; it cannot detect the second.

Running the diagnostic¶

The example runs RHF on water/cc-pVTZ, then an SCS-MP2 with cc-pvtz-ri and report_ri_residual=True, which reports the per-bucket fit error alongside the OS and SS correlation energies:

import vibeqc as vq

mol = vq.Molecule(
    [
        vq.Atom(8, [0.0,  0.00,  0.00]),
        vq.Atom(1, [0.0,  1.43, -0.98]),
        vq.Atom(1, [0.0, -1.43, -0.98]),
    ]
)
basis = vq.BasisSet(mol, "cc-pvtz")

rhf_opts = vq.RHFOptions(); rhf_opts.conv_tol_energy = 1e-10
hf = vq.run_rhf(mol, basis, rhf_opts)

opts = vq.MP2Options()
opts.density_fit = True
opts.aux_basis = "cc-pvtz-ri"
opts.c_os = 6.0 / 5.0          # SCS-MP2 (Grimme 2003)
opts.c_ss = 1.0 / 3.0
opts.n_frozen_core = 0          # all electrons; the default freezes the O 1s core
opts.report_ri_residual = True  # ← opt-in flag

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

assert mp2.ri_residual_reported
print(f"E_os = {mp2.e_os:+.8f} Ha   "
      f"ΔE_os = {mp2.e_os_ri_residual:+.2e}")
print(f"E_ss = {mp2.e_ss:+.8f} Ha   "
      f"ΔE_ss = {mp2.e_ss_ri_residual:+.2e}")

Output:

E_os = -0.20613056 Ha   ΔE_os = +5.20e-05
E_ss = -0.06591474 Ha   ΔE_ss = -2.43e-05

These are all-electron numbers from a vibe-qc 0.17.10 development build (each MP2 call takes about 3 s). With the default frozen core (drop the n_frozen_core line) the energies become E_os = -0.19535593 and E_ss = -0.06281455 Ha and the residuals +5.08e-05 and -2.38e-05 Ha, essentially unchanged. ΔE_os > 0 and ΔE_ss < 0 is the expected Coulomb-metric Dunlap fingerprint. With c_os = 6/5 and c_ss = 1/3:

\[ \Delta E_{\mathrm{SCS}} \approx 1.2 \cdot (+5.20 \times 10^{-5}) + 0.333 \cdot (-2.43 \times 10^{-5}) = +5.4 \times 10^{-5} \mathrm{\,Ha} \approx 0.034 \mathrm{\,kcal/mol}. \]

For chemical accuracy (1 kcal/mol) this is well inside tolerance; for chasing sub-µHa parity against a reference calculation it is not.

Aux-basis sweep, confirming a genuine fit residual¶

A fit residual should shrink monotonically as you enlarge the aux basis. Re-running on cc-pvqz-ri (242 functions vs 141 for cc-pvtz-ri):

Aux basis	\(n_{\mathrm{aux}}\)	\(\Delta E_{\mathrm{os}}\)	\(\Delta E_{\mathrm{ss}}\)	\(\Delta E_{\mathrm{SCS}}\)
cc-pvtz-ri	141	+5.2e-5	−2.4e-5	+5.4e-5
cc-pvqz-ri	242	+2.0e-5	−4.4e-6	+2.2e-5

Both rows show the same sign pattern; the magnitudes drop roughly 2×, the right signature for a fit-residual sweep.

If you run the same sweep and the residual stays flat or grows with aux size, that is not a fit-residual signature, it’s evidence of a code bug (Cholesky / SVD propagation, mis-routed weights), an ill-conditioned \(V_{PQ}\), or an aux family that is misaligned with the orbital basis. Open an issue with the diagnostic output and the \(V_{PQ}\) condition number (the latter is on vibeqc.density_fitting.DensityFitting(...).info.metric_condition_estimate).

When to pair this diagnostic with a code¶

Before publishing a new aux-basis recommendation. Run the diagnostic on a representative set; the per-system bucket residual is your “this aux is fit for this purpose” evidence.
Before cross-code parity work. If two codes’ RI-MP2 numbers disagree by more than the per-code RI residuals against canonical, one of them is doing something other than the standard Coulomb-metric fit (e.g. one might be running canonical silently, always confirm both codes are running RI as advertised).
As a unit-test guard. Pin the residual magnitude alongside the energy itself in your test suite, a regression in the RI pipeline shows up as a residual shift, not just an energy shift.

Cost note¶

The diagnostic doubles the runtime of an RI-MP2 call, it forces the O(N⁵) canonical AO→MO transform alongside the O(n_aux · …) RI build. Keep it off in production loops; only switch it on when you need the answer it gives.

opts.report_ri_residual = False   # ← default, production

References¶

RI-MP2 (the algorithm). F. Weigend, M. Häser, “RI-MP2: first derivatives and global consistency,” Theor. Chem. Acc. 97, 331 (1997). doi:10.1007/s002140050269. The original derivation of the Coulomb-metric-fitted MP2 amplitudes vibe-qc implements.
Dunlap-fit theory. B. I. Dunlap, “Robust and variational fitting: Removing the four-center integrals from center stage in quantum chemistry,” J. Mol. Struct. THEOCHEM 529, 37 (2000). doi:10.1016/S0166-1280(00)00528-5. The “robust” fitting framework that pins the sign pattern of the residual.
Aux-basis families. F. Weigend, A. Köhn, C. Hättig, “Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations,” J. Chem. Phys. 116, 3175 (2002). doi:10.1063/1.1445115 (the cc-pVNZ-RI family); F. Weigend, “Accurate Coulomb-fitting basis sets for H to Rn,” Phys. Chem. Chem. Phys. 8, 1057 (2006). doi:10.1039/B515623H (the def2-*-rifit family).
SCS-MP2. S. Grimme, J. Chem. Phys. 118, 9095 (2003). doi:10.1063/1.1569242, the scaling weights c_os = 6/5, c_ss = 1/3 used in this tutorial.

All four fire automatically in the bundled citation database routes for density_fit=True MP2 / c_os / c_ss jobs.