An introduction of the RI technique^2–7 conveniently starts from the definition〈f|g〉 = ∫ f(r₁)r₁₂⁻¹g(r₂)dτwhich fulfils all requirements of a scalar product: it is linear and positive definite, since 〈f|f〉 = 0 if and only if f = 0.

ERIs are then simply written as (ab|cd) = 〈ab|cd〉.

Let us now introduce a set of functions labelled P, Q and the projection operator onto this space where M_PQ denotes matrix elements of the inverse of 〈P|Q〉. is optimal in the sense is usually called ‘resolution of the identity’; its insertion into eqn. (1) yields the RI approximation for ERIs (using parentheses for two-electron integrals since this is common usage for the charge density notation) We similarly get a concise expression for the Coulomb energy J of a charge distribution ρ(r)which has an error (O||δρ||²) in the sense of eqn. (4).

The functions P, Q are usually termed auxiliary or fitting functions, since eqn. (6) implies that ρ is approximated as ρ ≈  = ρ.

The approximation in eqn. (5) decomposes four-index ERIs into two- and three-index-terms which is its most important feature.

The formal O(N⁴) effort to evaluate (ab|cd), N = number of basis functions, is thus reduced to O(N³).

If the auxiliaries P, Q are chosen as atom-centered functions one further has to deal with three-centre integrals only, which leads to additional simplifications as will be demonstrated below.

The errors introduced by the approximation in eqn. (5) are of no concern if optimized auxiliary functions are employed.

We consider the evaluation of (ab|P) for the usual choice of atom-centered GTOs (Gaussian type orbitals)|l_a〉 = a(r) = (x − A_x)^l_x(y − A_y)^l_y(z − A_z)^l_ze^{−α|r − A|2}l = (l_x,l_y,l_z), l = l_x + l_y + l_z.With the shorthand notation a(r) = |l_a〉 we drop the parameters A and α. b(r) = |l_b〉 will be similarly specified by l_b, β, B, and the auxiliary function P(r) = |L〉 by L, γ, C, and we write (ab|P) = (l_al_b|L).It is sufficient to consider (l_a0|L) since the general case is recovered by the horizontal recursion relation, e.g.

(l_a(l_b + 1_i|L) = ((l_a + 1_i)l_b|L) + (A_i − B_i)(l_al_b|L), with 1_i = (δ_ix, δ_iy, δ_iz,) for i = x,y,z.

Although we consider Cartesian Gaussians, eqn. (7), we will assume that |L〉 is always transformed later on to reduced (real) spherical harmonics comprising 5, 7, 9 components for d, f, g sets.

This can be done explicitly by a transformation step, which is advantageous if integrals have to be transformed into an MO representation as in correlated treatments.

For HF or DFT it is easier to ensure simply that contraction coefficients c_L of (l_al_b|L) do not contain components of s type for a d set, or of p type for an f set, etc.

The relevant equations and definitions are, if we stay close to the nomenclature of OS (keeping their P, since this can hardly be confused with the auxiliary functions):(00|0)^(m) = 2π^5/2(ζ + γ)^1/2(ζγ)⁻¹e^{−(αβ/ζ)|A−B|2}F_m(x)x = ρ|P − Q|²The true ERI has index m = 0, the recursion then requires m > 0 for intermediate quantities, i.e.m = 0 to (l_a + l_b + L) for the start, eqn. (10).

In general the easiest way is to first get the necessary (l0|0)^(m), eqn. (11), and then to increase L by means of eqn. (12).

The recursive increase of L can be simplified.

The first term on the rhs of eqn. (12) clearly vanishes since Q = C.

Next we exploit that components of |L〉 are (transformed to) spherical harmonics.

This implies the following asymptotic decay of integrals:(l0|L) ∝ |P − C|^−L−1 for |P − C| → ∞.(l0|L) may vanish even faster if l0 does not include a partial wave of s character but this is of no concern for the present considerations.

Withone can identify all contributions to the final integral that vanish too slowly, and these terms can be neglected since they cancel in the transformation to spherical harmonics.

From eqns. (10) and (11) one has (l0|0)⁽⁰⁾ ∝ |P − C|⁻¹The third term on the rhs of eqn. (12) maintains this asymptotic behavior, which in the final integral, m = 0, would lead to (l0|L) ∝ |P − C|⁻¹ in contradiction to eqn. (17).

For L = 2, i.e. a set of d functions, these terms give identical contributions ((1/2η)(l0|0)^(m) − (ρ/2η²)(l0|0)^(m+1)) to the integrals involving Cartesian functions d_x², d_y², d_z², and this cancels if one goes over to d_x²−y² and d_3z²−r².

We thus can replace the five-term recursion in eqn. (12) by a two-term recursion Since the index m on the lhs is connected only with (m + 1) on the rhs, one starts the recursion with (l0|0)^(L), then gets (l0|1)^(L−1), etc., until the final integral (l0|L)⁽⁰⁾ is reached.

For each L one has only a single m value and this index can simply be implied with corresponding savings in memory and overhead necessary to implement eqn. (20).

For the final batch, (l0|(L + 1_i))⁽⁰⁾ , one will not use eqn. (20) directly if integrals over contracted GTOs are computed.

It is more efficient to accumulate (l − 1_i)0|L)⁽¹⁾ separately and to add the sum later on.

This offers the advantage that ((l − 1_i)0|L)⁽¹⁾ has fewer components than (l0|(L + 1_i)).