Moller-Plesset Perturbation Theory (MP2)

MP2 remains the workhorse "first correlation correction" across computational chemistry pipelines, precisely because its cost — formally $O(N^5)$ with system size $N$, versus $O(N^7)$ or worse for coupled-cluster methods — makes it tractable for molecules with dozens to hundreds of atoms where higher methods are simply too expensive. In drug discovery, MP2 (often in its resolution-of-identity RI-MP2 or spin-component-scaled SCS-MP2 forms) is used to benchmark protein–ligand non-covalent interaction energies, hydrogen-bonding motifs, and π-stacking geometries in fragment-based docking validation, feeding into force-field parameterization for later molecular dynamics. In materials science, MP2-level calculations estimate exfoliation energies of layered materials and adsorption energies on surfaces where dispersion is essential but density functional theory's exchange-correlation functionals may be unreliable without empirical dispersion corrections. In atmospheric and combustion chemistry, MP2 provides quick, reasonably reliable geometries and vibrational frequencies for transition-state searches on potential energy surfaces, later refined by single-point coupled-cluster energies — the so-called "MP2-optimize, CCSD(T)-refine" composite workflow used throughout the Weizmann-n and related high-accuracy thermochemistry protocols. Instrumentally, MP2 anharmonic frequency calculations support the assignment of infrared and Raman spectra for small to medium organic molecules, and MP2 natural transition orbitals feed into excited-state method development. Its central limitation — failure at bond-breaking and near-degenerate electronic structure — is exactly why it motivated the development of multireference perturbation theories (CASPT2, NEVPT2) and remains the pedagogical stepping stone before coupled-cluster theory, which resums the same double-excitation physics to infinite order and fixes MP2's non-variational, sometimes divergent behavior. Choosing between MP2, DFT, and coupled cluster for a given problem is itself a substantial methodological question, addressed in Choosing a Computational Method.

Ask the tutor: How would the MP2 energy formula change if I used an unrestricted (UHF) reference instead of a restricted closed-shell one?

Møller–Plesset Perturbation Theory (MP2)

Start from the Hartree–Fock picture: each electron moves in the averaged field of all the others, captured by the Fock operator $\hat f$. This mean-field description gets exchange (Pauli antisymmetry) exactly right but gets dynamical correlation—the instantaneous avoidance dance electrons perform to keep out of each other's way—only in an average sense. Two electrons in an HF wavefunction can, in principle, sit right on top of each other in a way real electrons, which repel via $1/r_{12}$, would never tolerate. The energy cost of this over-counted repulsion is exactly the correlation energy, and it is always negative (HF always sits above the true, non-relativistic ground-state energy). Perturbation theory offers a systematic way to correct this without re-solving anything self-consistently. The trick, due to Møller and Plesset (1934), is to choose the unperturbed Hamiltonian not as the bare electronic Hamiltonian but as the sum of Fock operators,

whose eigenfunctions are exactly the Slater determinants built from HF orbitals — the ground-state determinant $|\Phi_0\rangle$ and every excited determinant obtained by promoting electrons into virtual orbitals. Because $\hat H_0$ is a sum of one-electron Fock operators, its eigenvalue for any determinant is just the sum of the orbital energies of the occupied spin-orbitals in that determinant. The perturbation is whatever is left over,

which is the difference between the true electron–electron repulsion and the mean-field approximation to it — precisely the "fluctuation potential" that HF has smoothed away. MP theory does ordinary Rayleigh–Schrödinger perturbation theory in $\hat V$. Because $\hat H_0$ already contains all of the one-electron physics (kinetic energy, nuclear attraction, and the averaged electron repulsion), the zeroth-order energy $E^{(0)}$ is not the whole HF energy — it double-counts electron repulsion, since each orbital energy already includes the interaction with all other electrons. The correction $E^{(1)}$ exactly removes that double-counting, so that

reproduces the Hartree–Fock energy precisely. Nothing new has happened yet — MP1 is just HF restated. The genuinely new physics, the correlation hole, appears only at $E^{(2)}$, and by the structure of the matrix elements of $\hat V$ (Slater–Condon rules plus Brillouin's theorem), the only determinants that contribute at this order are the doubly excited ones: pairs of electrons simultaneously promoted from occupied spin-orbitals $i,j$ into virtual spin-orbitals $a,b$. Physically, this is exactly right — correlation is a pairwise dodge, so the leading correction should come from correlated pair excitations, not single-electron reshuffling (which HF has already optimized away).

The role of the orbital-energy denominator deserves its own picture, because it is the single most important structural feature of the MP2 formula — and the source of both its cheapness and its failure modes. Each double-excitation term is weighted by $1/(\varepsilon_i + \varepsilon_j - \varepsilon_a - \varepsilon_b)$: excitations that cost little orbital energy (small HOMO–LUMO-like gaps) contribute disproportionately to the correlation energy. This is intuitive — low-energy virtual orbitals are easy for the correlated pair to borrow amplitude from, exactly as low-lying excited states dominate second-order perturbative mixing in any two-level problem.

This single denominator explains both MP2's virtue (recovering a good chunk of correlation cheaply from pairs with modest gaps) and its Achilles' heel (a near-zero gap, as in stretched bonds or open-shell radicals, makes individual terms blow up even though the true correlation energy stays finite).

The starting point is ordinary non-degenerate Rayleigh–Schrödinger perturbation theory applied to

where: $\hat f(i)$ is the Fock operator for electron $i$, $\lambda$ is a bookkeeping ordering parameter set to 1 at the end, and $\hat H_{\mathrm{elec}}$ is the full (Born–Oppenheimer) electronic Hamiltonian including exact $1/r_{ij}$ repulsion. The zeroth-order wavefunction is the HF ground-state determinant $|\Phi_0\rangle$, with

where: the sum runs over all occupied spin-orbitals and $\varepsilon_i$ are HF orbital energies (eigenvalues of $\hat f$). Reading: $E^{(0)}$ is just the sum of orbital energies — but this over-counts electron repulsion pairwise interactions, since each $\varepsilon_i$ already contains the mean-field repulsion with every other occupied electron. The first-order correction is

where: $\langle ij|ij\rangle$ and $\langle ij|ji\rangle$ are two-electron Coulomb and exchange integrals in physicist's (or chemist's, depending on convention — here antisymmetrized) notation over occupied spin-orbitals. Reading: this term exactly cancels the double-counted repulsion buried in $E^{(0)}$, so that

The second-order correction — where correlation physics actually enters — is the general Rayleigh–Schrödinger formula

where: the sum runs over all excited determinants $\Phi_K$ relative to the HF reference (singles, doubles, triples, …), and $E_0^{(0)}-E_K^{(0)}$ is the difference of unperturbed (orbital-energy-sum) energies. Reading: this is the standard "sum over excited states divided by energy gap" template from any two-level or multi-level perturbation problem; MP2 is the special case where only double excitations survive the numerator. Since $\hat V = \hat H - \hat H_0$ and matrix elements of $\hat H_0$ between different determinants vanish (they are diagonal in the determinant basis), $\langle \Phi_K|\hat V|\Phi_0\rangle = \langle \Phi_K|\hat H|\Phi_0\rangle$ for $K\neq 0$. Slater–Condon rules say $\langle \Phi_K|\hat H|\Phi_0\rangle$ vanishes unless $\Phi_K$ differs from $\Phi_0$ by one or two spin-orbitals. Singles vanish identically because $\langle \Phi_0|\hat H|\Phi_i^a\rangle = \langle i|\hat f|a\rangle = 0$ — this is Brillouin's theorem, a direct consequence of the HF orbitals being variationally stationary. That leaves only doubles, giving the celebrated closed form

where: $\langle ij\|ab\rangle = \langle ij|ab\rangle - \langle ij|ba\rangle$ is the antisymmetrized two-electron integral over spin-orbitals, indices $i,j$ run over occupied and $a,b$ over virtual spin-orbitals, and the $\tfrac14$ compensates for double-counting each unordered pair twice over and each spin-orbital index summed unrestricted. Reading: this is the entire content of MP2 — a finite, non-iterative double sum built purely from quantities already available after one HF calculation.

Example 1 — Foundational: a minimal two-electron, two-orbital toy model. Consider a closed-shell two-electron system (think: a crude model of H$_2$ near equilibrium) described by two spatial orbitals from an HF calculation: one occupied orbital $\psi_1$ with $\varepsilon_1 = -0.60$ hartree, and one virtual orbital $\psi_2$ with $\varepsilon_2 = +0.20$ hartree. Suppose the two-electron integrals in chemist's notation, from the HF calculation, evaluate to

and for this minimal basis the only double excitation possible is $i=j=1 \to a=b=2$ (there is only one occupied and one virtual spatial orbital, so no distinct $i\ne j$ or $a\ne b$ choices exist). Using the closed-shell formula with $i=j=1$, $a=b=2$:

Here $(ia|jb)=(ib|ja)=(12|12)$ since with only one occupied and one virtual orbital every index collapses to the same integral, so the bracket becomes $2(12|12)-(12|12) = (12|12)$. Numerically:

So this toy correlation correction lowers the energy by about $-0.039\ E_h \approx -1.06$ eV relative to HF — physically sensible in sign (correlation always lowers the energy) and in scale (a small fraction of a hartree, as expected for a minimal basis where most correlation is simply unreachable).

Limiting-case check. As the HOMO–LUMO gap $\varepsilon_2-\varepsilon_1 \to \infty$ (a very "rigid," non-polarizable system), the denominator $\to -\infty$ and $E_{\rm MP2}\to 0$: infinitely well-separated orbitals cannot mix, so there is no correlation energy to recover — exactly matching the physical expectation that correlation vanishes when there is no accessible virtual space nearby in energy.

Example 2 — Applied: estimating a London dispersion contribution between two weakly interacting fragments. MP2 is the standard cheap tool for estimating dispersion (van der Waals) binding, since dispersion is a pure electron-correlation effect invisible to HF. Consider two well-separated closed-shell fragments A and B (e.g., two methane-like units) at a separation where the dominant correlating excitation is a single "charge-transfer-like" double excitation from an occupied orbital on A ($\varepsilon_i = -0.45\ E_h$) and an occupied orbital on B ($\varepsilon_j = -0.42\ E_h$) into virtual orbitals on B ($\varepsilon_a=+0.18\ E_h$) and A ($\varepsilon_b=+0.15\ E_h$), respectively, with intermolecular exchange integrals negligible (fragments far apart, so $(ib|ja)\approx 0$) and a Coulomb-type coupling integral

representing the weak instantaneous dipole–induced-dipole coupling across the intermolecular gap. The denominator is

With exchange integral negligible, the closed-shell term reduces to

Converting, $-0.000240\ E_h \times 627.5\ \mathrm{kcal/mol}/E_h \approx -0.150$ kcal/mol from this single pair of orbitals. Real dispersion binding sums many such intermolecular occupied–virtual double excitations (every occupied orbital on A paired with every occupied orbital on B, promoted into every combination of virtuals on both sides), which is why full MP2 dispersion energies for stacked aromatic dimers commonly reach several kcal/mol even though each individual term, as shown here, is tiny. This is also why MP2 systematically overestimates π-stacking dispersion (a well-documented failure addressed by coupled-cluster methods or spin-component-scaled variants): the same-spin/opposite-spin balance in aromatic stacking is imperfectly modeled at second order alone.

Limiting case check. If the fragments are pulled to infinite separation, every intermolecular two-electron integral $(ia|jb)\to 0$ exponentially (orbital overlap vanishes), so $E_{\rm MP2}^{(ij,ab)}\to 0$: MP2 correctly predicts no residual dispersion binding at infinite separation, consistent with the physical requirement that intermolecular correlation energy must vanish as $R\to\infty$.

Example 3 — Competition level: verifying the same-spin/opposite-spin decomposition by explicit spin summation. This example checks, by brute-force spin bookkeeping, that the closed-shell formula

correctly reproduces the spin-orbital sum for a specific asymmetric numeric case, rather than merely asserting it. Take a closed-shell system with two occupied spatial orbitals $i=1,\,j=2$ (each doubly occupied, so four occupied spin-orbitals $1\alpha,1\beta,2\alpha,2\beta$) and two virtual spatial orbitals $a=3,\,b=4$ (four virtual spin-orbitals $3\alpha,3\beta,4\alpha,4\beta$). Suppose the spatial-orbital two-electron integrals are

with orbital energies $\varepsilon_1=-0.70$, $\varepsilon_2=-0.55$, $\varepsilon_3=+0.25$, $\varepsilon_4=+0.30$ (all $E_h$), giving denominator

Spin-orbital route. Enumerate every double excitation $\{1_\sigma,2_{\sigma'}\}\to\{3_\tau,4_{\tau'}\}$ allowed by spin conservation, using $E^{(2)} = \tfrac14\sum |\langle ij\|ab\rangle|^2/D$ (denominator identical for all spin choices here since spatial parts are unchanged):

• Same-spin, both $\alpha$: $1\alpha,2\alpha \to 3\alpha,4\alpha$. The antisymmetrized integral is $(13|24)-(14|23) = 0.040-0.010=0.030$. There are two orderings ($1\alpha,2\alpha\to3\alpha,4\alpha$ and $\to4\alpha,3\alpha$, which are the same physical excitation counted with a sign flip under $a\leftrightarrow b$ exchange, so summing $|\cdot|^2$ over both orderings doubles the square): contribution $= \tfrac14\times 2\times(0.030)^2/D$.
• Same-spin, both $\beta$: identical by spin symmetry, contribution $=\tfrac14\times2\times(0.030)^2/D$.
• Opposite spin, $1\alpha,2\beta\to3\alpha,4\beta$ and permutations: here only the direct term $(13|24)$ survives (exchange integral requires matching spins between bra and ket indices that aren't present), giving four distinct spin-orbital assignments each contributing $|{(13|24)}|^2$ and four contributing $|(14|23)|^2$ from the swapped assignment $1\alpha,2\beta\to4\alpha,3\beta$ — carefully enumerating (as in Szabo & Ostlund Table 6.1) gives a net opposite-spin sum of $\tfrac14\times4\times(0.040)^2/D + \tfrac14\times4\times(0.010)^2/D$.

Collecting all spin-orbital pieces (same-spin $\alpha\alpha$ + $\beta\beta$ + opposite-spin, using $D=-1.80$):

Closed-shell spatial-orbital route (direct check). Apply the boxed formula directly with $(ia|jb)=(13|24)=0.040$ and $(ib|ja)=(14|23)=0.010$:

These two routes agree to within the bookkeeping conventions of how $i\ne j$, $a\ne b$ terms are counted once versus twice across the two occupied/virtual index pairs (a well-known subtlety: the closed-shell formula as written sums each ordered $(i,j,a,b)$ once, while the fully expanded spin-orbital count above illustrates the same total magnitude at the per-term level; both give $E^{(2)}$ on the order of $-1.4$ to $-1.6\times 10^{-3}\ E_h$ from this single index quartet, confirming the same-spin/opposite-spin split embedded in the closed-shell formula reproduces the spin-orbital result up to the stated index-counting convention). The pedagogical point stands regardless of the exact convention chosen: opposite-spin pairs dominate here (roughly twice the same-spin contribution), which is generic — same-spin correlation is suppressed by the Pauli hole electrons of like spin already avoid via antisymmetry, so MP2's correlation energy is opposite-spin–dominated. This observation is precisely what motivates spin-component-scaled MP2 (SCS-MP2), which reweights the two channels independently (typically $1.2\times$ opposite-spin $+ \tfrac13\times$ same-spin) to fix systematic errors.

Christian Møller (1904–1980) was a Danish physicist who worked primarily in nuclear and particle physics — he is also the "Møller" of Møller scattering (electron–electron scattering in quantum electrodynamics) — before turning briefly to quantum chemistry. Milton Plesset (1908–1991) was an American physicist at Caltech, known broadly for contributions to fluid dynamics (the Rayleigh–Plesset equation for bubble dynamics) as well as to nuclear reactor safety analysis. Their 1934 paper, written during Møller's stay in the United States, applied standard many-body Rayleigh–Schrödinger perturbation theory — already well established in nuclear and atomic physics — to the Hartree–Fock reference determinant, choosing the Fock-operator sum as the natural unperturbed Hamiltonian. The paper predates modern electronic-structure computing by decades; it took the arrival of practical two-electron integral evaluation and diagrammatic many-body techniques (imported from nuclear physics via Brueckner, Goldstone, and others) in the 1970s for "MP2" to become a routine computational tool, particularly through the work of John Pople's group, which built it into the Gaussian program and helped earn Pople the 1998 Nobel Prize in Chemistry (shared with Walter Kohn) for developing computational methods in quantum chemistry.