“Know how to solve every problem that has been solved.” “What I cannot create, I do not understand.” — Richard Feynman

McMurchie-Davidson: change basis to where Coulomb is diagonal

Quantum Chemistry

⛓ What you need to know first 5 concepts, 4 layers

The requisite-knowledge inventory for this page, bottom-up: the primitives at the base, combined upward until you reach what this page assumes. Skim the layers you already own; start wherever the ground gets unfamiliar.

base
- Quantum mechanics (states & operators)
L1
- Particle in a box
- Two-electron integrals (ERIs)
L2
- Slater orbitals
L3
- Gaussian product theorem
↳you are here

McMurchie-Davidson is the alternative recurrence scheme for Gaussian-basis molecular integrals. The structure: instead of climbing the Cartesian-Gaussian angular-momentum ladder directly (as Obara-Saika does), change basis to Hermite Gaussians on the product center, where the Coulomb integral is essentially diagonal. The work splits into two cleanly separated pieces: a Cartesian-to-Hermite expansion (cheap, geometric) and a Hermite-integral table (Boys-function ladder). Contract the two to recover the original integral.

The conceptual reason this is worth doing: the Hermite expansion is structurally tidy — once you've built the $E_{t}^{ab}$ coefficients for a shell pair, they're reused across the whole contraction set. And the $R_{t uv}^{n}$ auxiliaries depend only on the geometry of the product centers $P, Q$ and the combined exponent $η = pq / (p + q)$ , not on angular momentum, so a single $R$ -table per shell-block feeds every Cartesian recombination on that block. Total work decouples into two independent recurrences instead of one intertwined one.

What's a Hermite Gaussian

A Hermite Gaussian centered at $P$ with exponent $p$ and Hermite indices $(t, u, v)$ is defined by differentiating a centered Gaussian with respect to the center:

Λ_{t uv} (r; p, P) = (\frac{\partial}{\partial P _{x}})^{t} (\frac{\partial}{\partial P _{y}})^{u} (\frac{\partial}{\partial P _{z}})^{v} e^{- p ∣ r - P ∣^{2}}

These are conveniently called "Hermite Gaussians" because the partial derivatives generate Hermite polynomials in $(r - P)$ times the original Gaussian — specifically $Λ_{t uv} (r; p, P) = (2 p)^{t + u + v} H_{t uv} (p (r - P)) e^{- p ∣ r - P ∣^{2}}$ in one convention. They form a basis for Gaussian-modulated polynomials just like Cartesian Gaussians do, but they're better-suited to multi-center integration because differentiation under the integral sign is clean.

The MD strategy in three steps

The full algorithm:

Cartesian → Hermite expansion. Each pair of Cartesian Gaussians on different centers, $χ_{a} (r - A) χ_{b} (r - B)$ , is written as a finite sum of Hermite Gaussians on the product center $P$ :

χ_{a} (r - A) χ_{b} (r - B) = K_{A B} t = 0 \sum a_{x} + b_{x} u = 0 \sum a_{y} + b_{y} v = 0 \sum a_{z} + b_{z} E_{t uv}^{ab} Λ_{t uv} (r; p, P)

The coefficients $E_{t uv}^{ab}$ are independent $(x, y, z)$ products, $E_{t uv}^{ab} = E_{t}^{a_{x} b_{x}} E_{u}^{a_{y} b_{y}} E_{v}^{a_{z} b_{z}}$ , each computed by a small 1D recursion (below). They're geometric — pure functions of the centers, exponents, and angular-momentum indices, not of any Coulomb operator. Two-electron integrals get two such expansions, one for the $(ab)$ pair on $P$ and one for the $(cd)$ pair on $Q$ .

Hermite-integral auxiliaries. The integral of Hermite Gaussians weighted by $1/ r_{12}$ has a closed form in terms of derivatives of the Boys function. Define

R_{t uv}^{n} (η, P - Q) = (\frac{\partial}{\partial P _{x}})^{t} (\frac{\partial}{\partial P _{y}})^{u} (\frac{\partial}{\partial P _{z}})^{v} (- 2 η)^{n} F_{n} (T)

with base case $R_{000}^{n} = (- 2 η)^{n} F_{n} (T)$ and a three-term recurrence

R_{t + 1, u, v}^{n} = t R_{t - 1, u, v}^{n + 1} + (P_{x} - Q_{x}) R_{t, u, v}^{n + 1}

(analogous on $u$ and $v$ ). This recurrence descends in $(t, u, v)$ indices while ascending in level $n$ , exactly as the Obara-Saika one does for Cartesian indices. The base case is the Boys function at increasing levels; everything is geometric (depends only on centers and exponents), not on the original $a, b, c, d$ angular-momentum vectors. So one $R_{t uv}^{n}$ table works for every shell on a given ( $P, Q$ ) pair.

Contraction. The final integral is

(ab ∣ cd) = \frac{2 π ^{5/2}}{pq p + q} K_{A B} K_{C D} t uv \sum τυν \sum E_{t uv}^{ab} E_{τυν}^{c d} (- 1)^{τ + υ + ν} R_{t + τ, u + υ, v + ν}^{0}

A finite sum of products: each $E_{t uv}^{ab}$ coefficient times each $E_{τυν}^{c d}$ coefficient times the $R^{0}$ value at the index sum. The $(- 1)^{τ + υ + ν}$ factor arises because the Hermite derivatives on electron 2 are with respect to $Q$ , which flips sign relative to the $P$ -derivatives that build the $R$ table.

The E-coefficient recurrence (1D)

The Hermite expansion factor in the $x$ direction satisfies

E_{t}^{i + 1, j} = X_{P A} E_{t}^{i, j} + \frac{1}{2 p} E_{t - 1}^{i, j} + (t + 1) E_{t + 1}^{i, j}

E_{t}^{i, j + 1} = X_{PB} E_{t}^{i, j} + \frac{1}{2 p} E_{t - 1}^{i, j} + (t + 1) E_{t + 1}^{i, j}

where $X_{P A} = P_{x} - A_{x}$ and $X_{PB} = P_{x} - B_{x}$ , starting from $E_{0}^{0, 0} = 1$ and $E_{t}^{0, 0} = 0$ for $t \neq = 0$ . The recurrence raises one of the Cartesian angular-momentum indices $(i, j)$ while building up the Hermite index $t$ . After climbing to the target $(a_{x}, b_{x})$ , the $E_{t}^{a_{x} b_{x}}$ values for $t = 0, 1, \dots, a_{x} + b_{x}$ are what you need. Analogous recurrences run independently on $(y, z)$ directions.

Worked example: $(p_{x} s ∣ ss)$

Compute the same $(p_{x} s ∣ ss)$ integral as the Obara-Saika page, with identical parameters ( $p_{x}$ on $A = (0, 0, 0)$ with $α = 0.5$ , $s$ on $B = (1, 0, 0)$ with $β = 0.3$ , etc.).

# Same integral as the Obara-Saika page: (p_x s | s s) on four distinct centers.
import math, numpy as np
from scipy.special import hyp1f1

def boys(n, T):
    return (1.0 / (2 * n + 1)) * hyp1f1(n + 0.5, n + 1.5, -T)

alpha, A = 0.5, np.array([0.0, 0.0, 0.0])     # p_x on A
beta_, B = 0.3, np.array([1.0, 0.0, 0.0])     # s   on B
gamma_, C = 0.4, np.array([0.0, 1.0, 0.0])    # s   on C
delta, D = 0.6, np.array([0.0, 0.0, 1.0])     # s   on D

p = alpha + beta_                              # = 0.8
q = gamma_ + delta                             # = 1.0
P = (alpha * A + beta_ * B) / p                # = (0.375, 0, 0)
Q = (gamma_ * C + delta * D) / q               # = (0, 0.4, 0.6)
eta = p * q / (p + q)                          # = 0.4444
K_AB = math.exp(-alpha * beta_ / p * np.sum((A - B)**2))   # = 0.8290
K_CD = math.exp(-gamma_ * delta / q * np.sum((C - D)**2))  # = 0.6188
T = eta * np.sum((P - Q)**2)                               # = 0.2936
prefactor = 2 * math.pi**2.5 / (p * q * math.sqrt(p + q)) * K_AB * K_CD

The MD pipeline for this case is short because most coefficients are zero. For the $(p_{x}, s)$ pair on $(A, B)$ :

$E_{0}^{0, 0} = 1$ , all others zero (the $(s, s)$ base case).
Apply the $i \to i + 1$ recurrence once with $X_{P A} = P_{x} - A_{x} = 0.375$ : $E_{0}^{1, 0} = X_{P A} \cdot 1 = 0.375$ , $E_{1}^{1, 0} = (1/ (2 p)) \cdot 1 = 0.625$ , others zero.

For the $(s, s)$ pair on $(C, D)$ , only $E_{0}^{0, 0} = 1$ is non-zero. So the contraction reduces to

(p_{x} s ∣ ss) = \frac{2 π ^{5/2}}{pq p + q} K_{A B} K_{C D} (E_{0}^{1, 0} R_{000}^{0} + E_{1}^{1, 0} R_{100}^{0})

Two terms. Build the $R_{t uv}^{n}$ auxiliaries needed: $R_{000}^{0} = F_{0} (T) = 0.9102$ , $R_{000}^{1} = - 2 η F_{1} (T) = - 0.1027$ , and the recurrence gives $R_{100}^{0} = (P_{x} - Q_{x}) R_{000}^{1} = 0.375 \cdot (- 0.1027) = - 0.0385$ . Plug in:

# ─── McMurchie-Davidson: Hermite expansion + R_tuv recurrence ──────────────
# Step 1: Cartesian-to-Hermite expansion coefficients for (p_x, s) on (A, B).
#         The s on B contributes nothing extra; (s, s) on (C, D) gives E^{0,0}_0 = 1.
#         Recurrence for E^{i+1, j}_t (in the i direction):
#             E^{i+1, j}_t = X_PA · E^{i, j}_t + (1/(2p))·E^{i, j}_{t-1} + (t+1)·E^{i, j}_{t+1}
#         At i = 0, j = 0:   E^{0,0}_0 = 1, all others zero.
#         Raise i to 1 (with j = 0, single x-step):
E_10_0 = P[0] - A[0]      # = 0.375  (the X_PA piece, since E^{0,0}_{-1}, E^{0,0}_1 are zero)
E_10_1 = 1.0 / (2 * p)    # = 0.625  (the 1/(2p) piece picks up E^{0,0}_0 = 1)

# Step 2: Hermite-integral auxiliaries R^n_{tuv}(eta, P-Q).
#         Base case:  R^n_{0,0,0} = (-2 eta)^n · F_n(T)
#         Recurrence: R^n_{t+1, u, v} = t · R^{n+1}_{t-1, u, v} + (P_x - Q_x) · R^{n+1}_{t, u, v}
R0_000 = boys(0, T)                               # = 0.9101805
R1_000 = (-2 * eta) * boys(1, T)                  # = -0.1027020
R0_100 = (P[0] - Q[0]) * R1_000                    # = -0.03851324

# Step 3: contract. Only t=0,1 and τυν=000 contribute (everything else multiplied by zero).
md_val = prefactor * (E_10_0 * R0_000 + E_10_1 * R0_100)
print(f"(p_x s | s s)  via MD  =  {md_val:.10e}")

(p_x s | s s)  via MD  =  4.7309112359e+00

Same $4.7309112359...$ as the Obara-Saika page, computed via completely different intermediates: OS used auxiliary integrals $[ss ∣ ss]^{(0)} = 15.22$ and $[ss ∣ ss]^{(1)} = 4.69$ ; MD used Hermite-expansion coefficients $E_{0}^{1, 0} = 0.375, E_{1}^{1, 0} = 0.625$ and Hermite-integral values $R_{000}^{0} = 0.9102, R_{100}^{0} = - 0.0385$ . Two different sets of intermediate quantities, two different recurrences, same number to machine precision.

What scales when you go to higher angular momentum

For a generic $(ab ∣ cd)$ integral, MD needs:

1D $E_{t}^{ij}$ tables in each Cartesian direction, dimension up to $(a_{x} + b_{x} + 1) \times (a_{x} + b_{x} + 1)$ in the $x$ direction etc. Computed by the small 3-term recurrence above. One pair of tables per shell pair on the $A, B$ side, one pair on the $C, D$ side.
3D $R_{t uv}^{n}$ table on the $(P, Q)$ geometry, dimensions up to $L_{m a x} + 1$ in each Hermite index where $L_{m a x} = L_{a} + L_{b} + L_{c} + L_{d}$ . Climbs from $R_{000}^{n} = (- 2 η)^{n} F_{n} (T)$ via the 3-term recurrence on $(t, u, v)$ while descending in $n$ .
Final contraction: a sum of products of $E$ coefficients with $R$ values. For an $L_{a} + L_{b}$ shell pair times an $L_{c} + L_{d}$ shell pair, this is $O ((L_{a} + L_{b} + 1)^{3} \cdot (L_{c} + L_{d} + 1)^{3})$ term sums.

The amortization advantage of MD is in step 2: the $R$ table is computed once per shell-block geometry and reused for every Cartesian recombination on that block. With many contracted primitives per basis function (typical for cc-pVDZ and beyond), this amortization can dominate the total cost.

Why the Hermite basis makes the Coulomb integral easy

The structural fact: the Coulomb kernel $1/ r_{12}$ has the Laplace-transform identity

\frac{1}{r _{12}} = \frac{2}{π} \int_{0}^{\infty} e^{- s^{2} r_{12}^{2}} d s

so the 2-electron Coulomb integral becomes a 1D $s$ integral of a product of two Gaussian integrals — one for each electron, both evaluated against a third Gaussian $e^{- s^{2} r_{12}^{2}}$ . After completing the squares (which invokes the Gaussian product theorem in 6D), the integrand reduces to a function of $(P - Q)$ and the combined exponent. The Hermite-Gaussian definition above — $(\partial / \partial P)^{(t uv)}$ acting on a centered Gaussian — converts the $(P - Q)$ dependence into derivative operations on the Boys function. That decouples the angular momentum from the radial integral cleanly, which is what makes the $R_{t uv}^{n}$ recurrence so much simpler than the OS one.

The OS vertical recurrence does essentially the same thing in Cartesian-Gaussian language, but it has to carry the polynomial pieces around explicitly through every step, which is why the eight-term form is so dense. MD pre-organizes the work into Hermite intermediates, runs the Boys ladder once, and only re-introduces the Cartesian structure at the very end via the $E$ coefficients.

When to use MD vs OS

Use McMurchie-Davidson when:

Angular momentum is high ( $L \geq 3$ , $f, g$ shells).
You need lots of derivative integrals (nuclear gradient or Hessian, geometric optimizations).
Many contractions share the same $(P, Q)$ geometry and you want the $R_{t uv}^{n}$ table amortized.
You're implementing relativistic or two-component methods where the integrand has extra derivative structure that MD's basis change handles naturally.

For routine low-angular-momentum work, Obara-Saika is usually faster (less storage, fewer auxiliary tables, simpler bookkeeping). Modern production codes typically use OS for the SCF inner loop and switch to MD for derivative integrals and high- $L$ shells.

Related on this site

Recursion strategies overview — the shared parent page framing OS and MD as two ladders for the same problem.
Obara-Saika — the direct-Cartesian alternative algorithm.
Gaussian Product Theorem — the identity that lets the $(s, s)$ Cartesian product collapse into a single Gaussian on $P$ in the first place.
Two-electron integrals (from scratch) — full pipeline for the H₂ molecule where everything is $s$ .

References

McMurchie, L. E. & Davidson, E. R. (1978). One- and two-electron integrals over Cartesian Gaussian functions. J. Comput. Phys., 26, 218.
Saunders, V. R. (1983). Molecular integrals for Gaussian-type functions. In Methods in Computational Molecular Physics, NATO ASI Series. (Detailed exposition of the MD scheme.)
Helgaker, T., Jørgensen, P., Olsen, J. (2000). Molecular Electronic-Structure Theory, ch. 9. Wiley.