Efficient evaluation of three-center two-electron integrals over Gaussian functions

The RI (resolution of the identity) technique achieves significant increases in efficiency for various molecular electronic structure methods.

This results from the approximation of four-center two-electron integrals by corresponding three-center integrals.

It is shown that the three-center integrals required can be evaluated with a much simpler algorithm than for the general case.

This further increases the advantage of RI procedures.


The evaluation of electron repulsion integrals (ERI) over basis functions a,b,c,d(ab|cd) = ∫ a(r1)b(r1)r12−1c(r2)d(r2)dτis a basic task in molecular electronic structure methods.

The computation of ERIs is typically the dominant step in HF (Hartree–Fock) and DFT (density functional theory) treatments of molecules if the electronic energy is computed exactly within a basis set expansion of molecular orbitals.

It is thus desirable to develop approximations for eqn. (1) which combine efficiency and accuracy.

The RI (resolution of the identity) technique is a proven procedure for this purpose.

The present paper gives a brief summary of the RI method at first, then treats the evaluation of three-center integrals within the Obara–Saika1 (OS) procedure, and finally presents new vertical recursion relations for spherical harmonics as auxiliary functions.

An introduction of the RI technique2–7 conveniently starts from the definition〈f|g〉 = ∫ f(r1)r12−1g(r2)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 MPQ 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||δρ||2) 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(N4) effort to evaluate (ab|cd), N = number of basis functions, is thus reduced to O(N3).

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.

This has been demonstrated in a series of investigations considering J, the HF exchange K, and correlated treatments based on second order perturbation theory.7–11

Gains in efficiency resulting from the RI technique are most pronounced for the treatment of J.

Other procedures have been developed for the very same purpose, such as CFMM (continuous fast multipole moment)12 in combination with the J engine method13 and Fourier techniques.14

Multipole moment methodology can also be exploited within RI, thus combining the advantages of both procedures.15

Three-center integrals

We consider the evaluation of (ab|P) for the usual choice of atom-centered GTOs (Gaussian type orbitals)|la〉 = a(r) = (xAx)lx(yAy)ly(zAz)lzeα|rA|2l = (lx,ly,lz), l = lx + ly + lz.With the shorthand notation a(r) = |la〉 we drop the parameters A and α. b(r) = |lb〉 will be similarly specified by lb, β, B, and the auxiliary function P(r) = |L〉 by L, γ, C, and we write (ab|P) = (lalb|L).It is sufficient to consider (la0|L) since the general case is recovered by the horizontal recursion relation, e.g.

(la(lb + 1i|L) = ((la + 1i)lb|L) + (AiBi)(lalb|L), with 1i = (δ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 cL of (lalb|L) do not contain components of s type for a d set, or of p type for an f set, etc.

The choice of reduced auxiliary functions is not only more aesthetic; it also improves numerical stability.

We demonstrate simplifications of integral evaluations for the OS scheme since this has been considered and implemented.

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(ζγ)−1e−(αβ/ζ)|AB|2Fm(x)x = ρ|PQ|2The true ERI has index m = 0, the recursion then requires m > 0 for intermediate quantities, i.e.m = 0 to (la + lb + 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 last step typically dominates by far the evaluation of a complete integral batch.

The recursive increase of L can be simplified.

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

This is trivial and the corresponding term has only been included to show differences to the general case (ab|cd).

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

This implies the following asymptotic decay of integrals:(l0|L) ∝ |PC|L−1 for |PC| → ∞.(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)(0) ∝ |PC|−1The third term on the rhs of eqn. (12) maintains this asymptotic behavior, which in the final integral, m = 0, would lead to (l0|L) ∝ |PC|−1 in contradiction to eqn. (17).

This implies that the third term on the rhs of eqn. (12) can be neglected since it cancels after transformation.

The fourth term on the rhs of eqn. (12) has the same structure as the third and cancellation applies here as well.

That the third and fourth term cannot matter can also be seen directly.

For L = 2, i.e. a set of d functions, these terms give identical contributions ((1/2η)(l0|0)(m) − (ρ/2η2)(l0|0)(m+1)) to the integrals involving Cartesian functions dx2, dy2, dz2, and this cancels if one goes over to dx2y2 and d3z2r2.

The same reasoning applies for sets of f, g, etc., functions.

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)(0) 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).

The reduction in the index range of m following from eqn. (20) also leads to a reduced index range in the recursion eqn. (11), and eqn. (10) is required only for Lm ≤ (L + la + lb).

It should be mentioned that there is an alternative to eqn. (12),16 which is more efficient for four-center integrals, especially for large angular momentum functions.

This advantage is lost for the present case of three-center integrals.

Two features are relevant for an implementation of eqn. (20).

For the intermediate integrals, m > 0, one does not need all components of a shell: dx2, dy2, dz2, and dxy suffice to get a complete f shell, and six f components for a g set, etc.

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

It is more efficient to accumulate (l1i)0|L)(1) separately and to add the sum later on.

This offers the advantage that ((l1i)0|L)(1) has fewer components than (l0|(L + 1i)).


We summarize the simplification resulting from the RI technique in comparison to conventional treatments based on (ab|cd).

(i) The original five term recursion eqn. (12) is replaced by a two-term recursion eqn. (20).

(ii) Different from the general case in eqn. (12) each intermediate batch (l0|L) occurs only with a single m value; this reduces storage and computational requirements.

(iii) The simple form of eqn. (20) greatly facilitates the development of optimized hand-coded routines.

We note in passing that even larger simplifications can be achieved if only two-center integrals have to be computed17 and the present case lies in between the two-center and four-center cases.

The above described algorithm, including optimized hand-coded routines for la + lb + L ≤ 4, has been implemented in TURBOMOLE18 and will be included in the next release.

The efficiency of the procedure can be seen from timings of an example at hand, a treatment of Cu40Sb21(PH3)28, symmetry T, with 173 atoms, 1905 contracted GTOs in an SV(P) basis,19 and 5184 (contracted) auxiliary functions.7

The evaluation of the interelectronic Coulomb interaction J required (per iteration): 64 s in direct mode (all integrals computed on the fly); 46 s if 660 MB was used to store three-center integrals in memory; and 20 s if the far-field integrals were obtained by multipole moment techniques,15 all timings for a 2.4 GHz Xeon.