Obara-Saika: climb the angular-momentum ladder directly

Obara-Saika is the most direct recurrence scheme for computing Gaussian-basis molecular integrals. The structure: start from $(ss|ss{)}^{(m)}$ (the Boys-function auxiliary integral — closed form) and climb up the Cartesian angular-momentum ladder one increment at a time. Each step adds one unit of angular momentum on one center in one Cartesian direction, and expresses the new integral as a linear combination of integrals you've already computed. Eight terms in the most general form, fewer when angular momenta are zero on the unaffected centers.

The conceptual difficulty in reading Obara-Saika in the original paper or in Helgaker chapter 9 is that the recurrences are stated in maximum generality (eight terms, four centers, all Cartesian components) without first pinning down what's going on. This page works through the simplest non-trivial case — a $(p_{x}s|ss)$ integral on four distinct centers, which needs one recurrence step — and shows why the eight-term general form is the eight-term general form.

Start here: the recurrence on the simplest integral

Before the four-center two-electron integral, meet the recurrence on the simplest integral it applies to: the overlap of two 1D Gaussians carrying angular momentum, $S(i,j)=\int (x-A{)}^i{e}^{-\alpha (x-A{)}^2}(x-B{)}^j{e}^{-\beta (x-B{)}^2}dx$. Same machine, three terms instead of eight. The bottom rung $S(0,0)$ is just two s-type Gaussians — that is the Gaussian product theorem, a closed form. Every rung up adds one unit of angular momentum as a shift relative to the product center $P$ plus lower-rung corrections:

# The 3-term baby version: the OVERLAP of two 1D Gaussians WITH angular momentum.
# g_i(x) = (x - A)^i exp(-alpha (x - A)^2). Same recurrence as the real integral, fewer terms.
import numpy as np
from scipy.integrate import quad

A, B, alpha, beta = -0.4, 0.7, 1.3, 0.9
p  = alpha + beta
P  = (alpha*A + beta*B) / p             # product-Gaussian center (Gaussian product theorem)
mu = alpha*beta / p
S00 = np.sqrt(np.pi/p) * np.exp(-mu*(A-B)**2)   # BASE RUNG: two s-Gaussians, closed form

memo = {(0, 0): S00}
def S(i, j):                            # Obara-Saika overlap recurrence
    if i < 0 or j < 0: return 0.0
    if (i, j) in memo: return memo[(i, j)]
    if i > 0:        # raise the first index, reach back down by one
        v = (P - A)*S(i-1, j) + ((i-1)*S(i-2, j) + j*S(i-1, j-1)) / (2*p)
    else:            # i is maxed; raise the second index instead
        v = (P - B)*S(i, j-1) + (i*S(i-1, j-1) + (j-1)*S(i, j-2)) / (2*p)
    memo[(i, j)] = v
    return v

def brute(i, j):                        # ground truth: just integrate it numerically
    f = lambda x: (x-A)**i*np.exp(-alpha*(x-A)**2) * (x-B)**j*np.exp(-beta*(x-B)**2)
    return quad(f, -25, 25)[0]

for (i, j) in [(0,0), (1,0), (2,0), (2,2), (4,4)]:
    print(f"S({i},{j}): OS ={S(i,j): .9f}   brute ={brute(i,j): .9f}")

S(0,0): OS = 0.627906671   brute = 0.627906671
S(1,0): OS = 0.282558002   brute = 0.282558002
S(2,0): OS = 0.269857162   brute = 0.269857162
S(2,2): OS = 0.073246124   brute = 0.073246124
S(4,4): OS = 0.102759009   brute = 0.102759009

Every value matches direct numerical integration to machine precision, up to $S(4,4)$ — the recurrence is exact, not an approximation. The $(p_{x}s|ss)$ integral below is this same ladder with more moving parts: two electrons (so indices on both), four centers, and the $1/r_{12}$ operator, which replaces the Gaussian base rung with the Boys function and adds an auxiliary index. That is the only reason the general form grows from three terms to eight.

The shape of the recurrence

A primitive Cartesian Gaussian on center $A$ with angular-momentum quantum numbers $\mathbf{a}=(a_x,a_y,a_z)$ is

Total angular momentum $L_a=a_x+a_y+a_z$: 0 for s, 1 for p, 2 for d, etc. A 2-electron repulsion integral on four primitives is

For $\mathbf{a}=\mathbf{b}=\mathbf{c}=\mathbf{d}=(0,0,0)$ (all s) this is closed-form via the Boys function. For higher angular momentum, the polynomial prefactors don't go away under the Gaussian product theorem — they survive, get multiplied together, and the resulting integrand is (high-degree polynomial) × Gaussian × $1/r_{12}$. The Obara-Saika recurrence is a way to compute these by induction on the angular momentum.

The auxiliary integral and the level $m$

The recurrence requires a one-parameter family of integrals, not just one integral. Define the auxiliary integral at level $m$ as

At $m=0$ this is the actual physical integral $(\mathbf{ab}|\mathbf{cd})$. The higher-$m$ auxiliaries are mathematical bookkeeping — they don't correspond to any physical quantity by themselves, but they appear as intermediate terms when you do the angular-momentum recurrence. Each angular-momentum step shuffles between $m$ and $m+1$, and at the end of the climb you read off the $m=0$ result.

The base case is

where $p=\alpha +\beta$, $q=\gamma +\delta$, $K_{AB}=e^{-\alpha \beta |\mathbf{A}-\mathbf{B}{|}^2/p}$, same for $K_{CD}$, $T=\eta |\mathbf{P}-\mathbf{Q}{|}^2$, $\eta =pq/(p+q)$, and $\mathbf{P},\mathbf{Q}$ are the Gaussian-product centers on each electron pair. $F_m(T)$ is the Boys function at level $m$.

The vertical recurrence (raising on one center)

The vertical recurrence relation (VRR) raises angular momentum on center $A$ by one unit in the Cartesian direction $i\in \{x,y,z\}$, with all other angular momenta fixed at zero on $B$ and $D$:

where $N_i(\mathbf{a})$ is the $i$-th angular-momentum component of $\mathbf{a}$, $\mathbf{W}=(p\mathbf{P}+q\mathbf{Q})/(p+q)$ is the four-center "Coulomb-weighted" center, and $1_i$ is the unit vector in the $i$ direction. The structure: the new integral on the LHS equals a geometric piece ($P_i-A_i$ and $W_i-P_i$ factors with lower-angular-momentum integrals at level $m$ and $m+1$) plus correction terms involving integrals at angular momentum $\mathbf{a}-1_i$ (one less than the input) — which is why this is an induction-on-angular-momentum scheme.

The horizontal recurrence (transferring between centers)

The vertical recurrence only raises on $A$ and $C$ — it never touches the $B$ or $D$ angular momenta. To handle integrals where $\mathbf{b}$ or $\mathbf{d}$ is non-zero, OS uses the horizontal recurrence relation (HRR), which transfers angular momentum between centers on the same electron:

No $(m)$ superscript: the HRR is geometric, depending only on the centers $\mathbf{A},\mathbf{B}$, not on the Coulomb operator. So the strategy is: VRR up to $[\mathbf{a}0|\mathbf{c}0]$ at all needed levels $m$, then apply HRR (level-free) to redistribute angular momentum onto $B$ and $D$.

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

Compute a $(p_{x}s|ss)$ integral on four distinct centers with the parameters below. The four basis primitives are: $p_x$ on $A=(0,0,0)$ with exponent $\alpha =0.5$; $s$ on $B=(1,0,0)$ with $\beta =0.3$; $s$ on $C=(0,1,0)$ with $\gamma =0.4$; $s$ on $D=(0,0,1)$ with $\delta =0.6$.

# Concrete integral: (p_x s | s s) on four distinct centers.
import math, numpy as np
from scipy.special import hyp1f1

def boys(n, T):
    """F_n(T) = (1/(2n+1)) · M(n+1/2; n+3/2; -T)   via confluent hypergeometric."""
    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

# Gaussian product centers and combined exponents
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)                          # = 8/18 = 0.4444
W = (p * P + q * Q) / (p + q)                  # = (1/6, 2/9, 1/3)
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

Since $\mathbf{b}=\mathbf{d}=0$ already, no HRR is needed. We need one VRR step: raise $A$-angular-momentum from $(0,0,0)$ to $(1,0,0)$, in the $x$ direction. With $\mathbf{a}=0$ on the RHS, the $N_i(\mathbf{a})=N_x(0)=0$ terms vanish. With $\mathbf{c}=0$, the $N_i(\mathbf{c})=0$ term vanishes too. What remains is the bare geometric piece:

Two terms, each a number times a previously-computed base-case auxiliary. Plug in the numerics — $P_x-A_x=0.375$, $W_x-P_x=-0.2083$, $[ss|ss{]}^{(0)}=15.220$, $[ss|ss{]}^{(1)}=4.688$ — and we get the integral.

# ─── Obara-Saika vertical recurrence ────────────────────────────────────────
# Base case: [ss|ss]^(m) = prefactor · F_m(T)
ss_ss_m0 = prefactor * boys(0, T)              # = 1.521999e+01
ss_ss_m1 = prefactor * boys(1, T)              # = 4.687611e+00

# Vertical step in the x-direction, raising angular momentum on center A
# from a = (0,0,0) to a = (1,0,0):
#   [(a+1ᵢ)0|c0]^(m) = (Pᵢ-Aᵢ)·[a0|c0]^(m) + (Wᵢ-Pᵢ)·[a0|c0]^(m+1) + (Nᵢ terms = 0)
# i = x = 0
px_s_ss = (P[0] - A[0]) * ss_ss_m0 + (W[0] - P[0]) * ss_ss_m1
print(f"(p_x s | s s)  via OS  =  {px_s_ss:.10e}")

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

That's it. One recurrence step, two terms, the answer comes out to machine precision. The McMurchie-Davidson page computes the same integral through a completely different set of intermediates and lands on the same $\mathrm{4.7309112359...}$ number — a strong cross-check on both implementations.

What scales when you go to higher angular momentum

For a generic $(\mathbf{ab}|\mathbf{cd})$ integral, OS needs:

1. Boys-function evaluations at levels $m=0$ through $L_a+L_b+L_c+L_d$. One call to the confluent hypergeometric routine per level, or one downward recurrence from $F_{m_{\max }}$ using $F_{m-1}(T)=(2T{F}_m(T)+e^{-T})/(2m-1)$.
2. VRR climbs up the $(\mathbf{a},\mathbf{c})$ angular-momentum lattice from $(0,0)$ to $(L_a+L_b,L_c+L_d)$. The intermediate integrals $[\mathbf{a}0|\mathbf{c}0{]}^{(m)}$ are stored in a 4-index table (3 Cartesian indices per angular vector + 1 level index per center pair).
3. HRR redistributions from $(\mathbf{a}+\mathbf{b},0,\mathbf{c}+\mathbf{d},0)$ down to $(\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d})$. Pure linear-algebra recombinations, no level dependence.

The work per integral grows roughly as $(L+1{)}^4$ for a target angular momentum $L$, but contraction effects (multiple primitives per contracted basis function) dominate the real-world cost at typical basis-set sizes. The Boys-function table is the same for every shell-block; the OS algebra is the per-shell-block work.

Why the vertical recurrence works

The structural reason OS works: the polynomial prefactor $(x-A_x{)}^{a_x}$ can be raised by one unit using a differentiation identity on the Gaussian itself:

So differentiating with respect to a Gaussian's center coordinate multiplies the function by $2\alpha (x-A_x)$ — effectively raising angular momentum by one. Combined with the chain rule (the Boys function inside the integral depends on $T=\eta |\mathbf{P}-\mathbf{Q}{|}^2$ and hence on $\mathbf{A}$ via $\mathbf{P}$), differentiating a complete integral $[\mathbf{a}0|\mathbf{c}0{]}^{(m)}$ with respect to $A_i$ yields exactly the right-hand side of the vertical recurrence above, with the $[\cdots {]}^{(m+1)}$ terms coming from differentiating the Boys function (which raises its level by one). The eight terms are the eight contributions from a product-rule expansion of $\partial /\partial A_i$ acting on a product structure in the integrand.

Helgaker derives this from generating-function manipulations on the Gaussian product — equivalent and cleaner-looking but less intuitive. The bare differentiation-and-product-rule derivation is in Obara & Saika 1986 §II.

When to use OS vs MD

Use Obara-Saika when:

• Angular momentum is moderate ($L\le 3$, i.e. up to $f$ orbitals).
• You want minimal intermediate storage.
• The integrand is the bare Coulomb operator (no derivative operators, no relativistic corrections).

McMurchie-Davidson becomes preferable when you have high $L$ with lots of derivatives (geometric forces, nuclear gradients) or when you want the recurrence ladder to be shared across a whole shell-block at once — the $R_{tuv}^n$ auxiliaries are $\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d}$-independent so they amortize beautifully when you have many primitives in a contraction.

Related on this site

• Recursion strategies overview — the shared parent page framing OS and MD as two ladders for the same problem.
• McMurchie-Davidson — the alternative algorithm.
• Gaussian Product Theorem — the identity that gives the $[ss|ss{]}^{(m)}$ base case in closed form.
• Lanczos and the three-term recurrence — the same skeleton: this angular-momentum ladder is the Hermite recurrence, which is the tridiagonal "multiply-by-$x$" operator behind orthogonal polynomials and Lanczos.
• Two-electron integrals (from scratch) — full pipeline for the H₂ molecule where everything is $s$ and the OS machinery is invisible.

References

• Obara, S. & Saika, A. (1986). Efficient recursive computation of molecular integrals over Cartesian Gaussian functions. J. Chem. Phys., 84, 3963.
• Head-Gordon, M. & Pople, J. A. (1988). A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations. J. Chem. Phys., 89, 5777.
• Helgaker, T., Jørgensen, P., Olsen, J. (2000). Molecular Electronic-Structure Theory, ch. 9. Wiley.