The rys quadrature revisited: a novel formulation for the efficient computation of electron repulsion integrals over gaussian functions

The Rys quadrature revisited: A novel formulation for the efficient computation of electron repulsion integrals over Gaussian functions Michel Dupuis Pacific Northwest National Laboratory, Environmental Molecular Sciences Laboratory, K8-91, P.O. Box 999, Battelle Boulevard, Richland, Washington 99352 Antonio Marquez Department of Chemistry, University of Sevilla, Sevilla, Spain 共Received 7 August 2000; accepted 7 November 2000兲 A novel formulation of the Rys quadrature algorithm for the calculation of the electron repulsion integrals over Gaussian basis functions is presented. The new algorithm is specifically designed for high contractions. As for the original Rys quadrature algorithm, the new algorithm is very efficient for high angular momentum functions. In addition it is also equally efficient for low angular momentum functions. The new algorithm takes unique advantage of 共1兲the numerical Rys quadrature methodology in 共2兲dealing with charge distributions a la McMurchie–Davidson and in 共3兲scaling integral blocks as a means of transferring angular momentum a la Gill–Head–Gordon– Pople. An analysis of the algorithm suggests very favorable floating-point operation counts. ©2001 American Institute of Physics. 关DOI: 10.1063/1.1336541兴 I. INTRODUCTION The evaluation of the electron repulsion integrals, 共 ␹ i ␹ j 储 ␹ k ␹ l兲⫽ 冕冕 ␹ i共1兲 ␹ j共1兲1 r12 ␹ k共2兲 ␹ l共2兲dr1dr2,共1兲 where ␹ i, ␹ j, ␹ k, and ␹ lare one-electron Gaussian basis functions, and of the derivatives of these integrals with respect to their centers, continues to be one of the major computational tasks in ab initio molecular orbital theory 共MO兲 and density functional theory 共DFT兲methods using atomcentered Gaussian bases. In traditional algorithms where the electron repulsion integrals are computed and stored on disk, their computation and handling is most often the timelimiting step of the calculation. This is even more so with the ‘‘direct’’ algorithms where the integrals are computed several times, whenever they are needed. The new generation of basis sets, with their large degree of contraction of primitive Gaussian functions and very high angular momentum functions, bring also a heavy computational demand. Thus advances in integral computation technology will continue to help push the limits of molecular calculations on large molecules and on molecules with many heavy elements. In atom-centered bases the necessary angular dependence of a basis function is achieved explicitly through the use of spherical harmonics or equivalently through the use of integer powers of the Cartesian coordinates. When the latter is combined with a Gaussian radial factor, as suggested by Boys,1a primitive basis function takes the form ␩ i共1,ai,I兲⫽共x1⫺xi兲nxi共y1⫺yi兲nyi共z1⫺zi兲nzi ⫻exp共⫺air1i 2兲,共2兲 where ‘‘1’’ denotes the electron with coordinates (x1,y1,z1), (xi,yi,zi) represent the coordinates of the basis function center I,r1iis the distance between electron ‘‘1’’ and the center I, and aiis the Gaussian exponent. The sum of the powers, ␭i⫽nxi⫹nyi⫹nzi 共3兲 is closely related to the total angular momentum quantum number. In practice the basis functions ␹ ’s are linear combinations of primitive basis functions, ␹ i共1,I兲⫽兺 ai K cai ␩ i共1,ai,I兲,共4兲 where K共not to be confused with center Kof ␹ kand ␩ k兲is called the degree or level of contraction. The attractive feature of Gaussian functions resides in the fact that there exist practical closed-form expressions for the electron repulsion integrals.1,2 A major breakthrough in electron repulsion integral technology was the method introduced by Pople and Hehre 共PH兲,3a method that proved to be very efficient for lowangular momentum highly contracted basis functions. The successive rotations involving local axes, albeit extremely efficacious for sand pfunctions, proved however difficult to extend to higher angular momentum functions. McMurchie and Davidson 共MD兲4proposed a method highly efficient for contracted functions and readily extendable to high angular momentum functions. Dupuis, Rys, and King 共DRK兲5devel- oped a radically different method based on an exact numerical quadrature using roots and weights generated from a set of orthogonal polynomials, the Rys polynomials. The implementation of the DRK method proved efficient for basis functions with a low degree of contraction. In the DRK approach the six-dimensional primitive integrals take the form of sums of products of three two-dimensional integrals, JOURNAL OF CHEMICAL PHYSICS VOLUME 114, NUMBER 5 1 FEBRUARY 2001 20670021-9606/2001/114(5)/2067/12/$18.00 © 2001 American Institute of Physics Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 each associated with a Cartesian direction. Rys et al.6devel- oped also several recursion formulas for the elementary twodimensional integrals. Over 10 years after the work of MD and DRK, Obara and Saika 共OS兲7proposed some new recursion formulas for the six-dimensional integrals of Eq. 共1兲. These recursion formulas turned out to be closely related to those of DRK, except that they apply to the six-dimensional primitive integrals directly. OS’s work generated a new wave of research in electron repulsion evaluation methodologies with significant contributions from Head-Gordon and Pople 共HGP兲,8 Gill, Head-Gordon, and Pople 共GHGP兲,9Lindh, Ryu, and Liu 共LRL兲,10 all algorithms that have received much attention, and the lesser-known algorithms by Ten-no,11 and by Ishida 共ISH兲.12 HGP improved OS’s method by applying the transfer relations, already noted in DRK’s work, outside the quadruple sum over contractions, so that only primitive integrals where ␩ jand ␩ lare of stype need be computed inside the contraction loops. This idea was taken further in GHGP’s method where additionally these authors borrowed the proposal of MD and limited the mandatory primitive integrals to those where only ␩ iis not an sorbital. Generation of other intermediate integrals is carefully accomplished via a socalled ‘‘scaling’’ of integrals, so that much intermediate work is moved stepwise outside the contraction loops. The floating-point operation count 共FLOP兲shows that GHGP’s method is indeed superior for contracted functions. Note that the FLOP counts reported in Refs. 8, 9, 11 for the DRK method are all erroneous. Only the count given in Ref. 10 is approximately correct, in accord with the efficiency observed in practice with the DRK algorithm.13,14 Finally, Ishida’s algorithm12 has actually an even better FLOP count than GHGP’s although it looks to be very complex to implement in practice. The driving force for the present research was the belief that applying recursion formulas to two-dimensional integrals and then combining these integrals to form sixdimensional integrals ought to be more efficient than applying somewhat more complex recursion formulas directly to six-dimensional integrals. Comparison of the LRL method10 with the GHGP method8supports this belief. It clearly shows that for low contraction level and high angular momentum functions, the Rys quadrature technique as implemented by LRL and in some cases as implemented by DRK remains the method of choice. The superiority of the GHGP algorithm9 over the LRL algorithm for even moderately contracted low angular momentum functions is however quite clear. In the present work we aimed at devising an algorithm based on the Rys quadrature that would be especially efficient at handling highly contracted basis functions for all levels of basis function angular momentum. The algorithm that we outline below succeeds in doing this. In what follows we give a detailed description of how the new algorithm works for contracted sand pfunctions, because for such basis functions the presentation remains tractable and helps understand how and why the algorithm works. Modern basis sets do emphasize high angular momentum basis functions. The new algorithm will work in these cases as well. We note however that the higher the angular momentum, the lower the effective contraction level for any group of four shells. Thus it appears that an optimal strategy might involve two 共or more兲different algorithms, one for moderately to highly contracted low angular momentum functions 共s,p, most likely d兲, combined with one for lightly contracted high angular momentum basis functions 共possibly dand higher angular momentum兲. The ‘‘revisited’’ Rys quadrature algorithm described below combines the following features: 共1兲the numerical quadrature of DRK; 共2兲the ‘‘scaling’’ technique of GHGP; and 共3兲the transfer equations of HGP. The deterministic numerical quadrature step replaces the complex tree algorithm of GHGP, and relies on simple recursion relations for the two-dimensional integrals. In addition the numerical quadrature makes it possible to evaluate the needed elementary integrals on centers Iand Kdirectly rather than on the exponent-weighted charge-distribution centers. In this manner the ‘‘scaling’’ procedure is simplified. The resulting algorithm is anticipated to have a very attractive FLOP count. In Sec. II we introduce the notations and give the basic equations of the numerical quadrature for completeness. The ‘‘scaling’’ procedure is described in Sec. III along with a detailed outline of the algorithm. Conclusions and outlook are given in Sec. IV. In what follows we adopt the notation of GHGP.9Parentheses are used for integrals over contracted functions, square-brackets are used for primitive integrals, a mixed notation is used to indicate that the contraction summation has been carried out for only one of the two charge distributions, either ␹ i ␹ jor ␹ k ␹ l. II. THE RYS QUADRATURE A. The numerical quadrature It has been known since the early work of Boys1that electron repulsion integrals over primitive Gaussian functions can be expressed in the form, 关 ␩ i ␩ j 储 ␩ k ␩ l兴⫽兺 m⫽0 L CoefmFm共X兲,共5兲 where Fm共X兲⫽ 冕 0 1t2mexp共⫺Xt2兲dt,共6兲 and L⫽␭i⫹␭j⫹␭k⫹␭l.共7兲 The value of Xdepends upon the exponential parameters ai, aj,ak, and al, and the positions of the centers of the four Gaussians, but it is independent of the angular momentum indices, the nx’s, ny’s, and nz’s. It follows immediately that Eq. 共5兲can be written, 关 ␩ i ␩ j 储 ␩ k ␩ l兴⫽ 冕 0 1PL共t2兲exp共⫺Xt2兲dt,共8兲 where PL(t2) is a polynomial of degree Lin t2with coefficients Coefm. The integral in Eq. 共8兲can be evaluated exactly by an N-point numerical quadrature formula, 2068 J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 M. Dupuis and A. Marquez Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 关 ␩ i ␩ j 储 ␩ k ␩ l兴⫽兺 ␣ ⫽1 N PL共t ␣ 2兲W ␣ ,共9兲 where Nis any integer satisfying the condition, N⬎L 2,共10兲 t ␣ is a positive root of the Nth polynomial RN(t,X), RN共t ␣ ,X兲⫽0, 共11兲 and W ␣ is a weight factor, which depends upon the value of X.RN(t,X) is the member of degree 2Nof a set of polynomials 共the Rys polynomials兲of even order in the variable t, orthonormal in the interval 关0,1兴with respect to the weight factor exp(⫺Xt2), 冕 0 1RN共t,X兲RM共t,X兲exp共⫺Xt2兲dt⫽ ␦ NM .共12兲 In the early work of King and Dupuis,15 the roots and weights were evaluated by polynomial approximations valid over finite intervals of Xalthough interpolation approaches should prove more efficient. B. The polynomial P L „ t 2… We express the Coulomb operator 1/r12 as a Gaussian transform, 1/r12⫽2 冑 ␲ 冕 0 ⬁exp共⫺u2r12 2兲du.共13兲 A change in order of integration leads to 关 ␩ i ␩ j 储 ␩ k ␩ l兴⫽2 冑 ␲ 冕 0 ⬁ 关 ␩ i ␩ j 兩 exp共⫺u2r12 2兲 兩 ␩ k ␩ l兴du, 共14兲 and the integrand in Eq. 共14兲factors into the product of three two-dimensional integrals associated with the three axes of a coordinate system, b ␩ i ␩ j 兩 exp共⫺u2r12 2兲 兩 ␩ k ␩ lc⫽Ix ⬘Iy ⬘Iz ⬘.共15兲 Of the twelve nx, ny, nz indices, only four nx values enter into the Ix ⬘term, Ix ⬘共nxi,nxj,nxk,nxl兲 ⫽ 冕冕 共x1⫺xi兲nxi共x1⫺xj兲nxj共x2⫺xk兲nxk共x2⫺x1兲nxl ⫻exp共⫺Qx兲dx1dx2,共16兲 where Qx⫽ai共x1⫺xi兲2⫹aj共x1⫺xj兲2⫹ak共x2⫺xk兲2 ⫹al共x2⫺xl兲2⫹u2共x1⫺x2兲2.共17兲 Here xiis the xcoordinate of the center Iof primitive ␩ iand nxi is the corresponding nx angular momentum index, and similar notation is used for ␩ j, ␩ k, and ␩ l. We define some new quantities, xA⫽共aixi⫹ajxj兲/共ai⫹aj兲,共18兲 xB⫽共akxk⫹alxl兲/共ak⫹al兲,共19兲 A⫽ai⫹aj,共20兲 B⫽ak⫹al,共21兲 ␳ ⫽AB/共A⫹B兲,共22兲 Dx⫽ ␳ 共xA⫺xB兲2,共23兲 Gx⫽共aiaj/共ai⫹aj兲兲共xi⫺xj兲2⫹共akal/共ak⫹al兲兲共xk⫺xl兲2. 共24兲 In terms of these, Qxbecomes Qx⫽Gx⫹A共xl⫺xA兲2⫹B共x2⫺xB兲2⫹u2共x1⫺x2兲2.共25兲 We then make a change of variable from uto t, u2⫽ ␳ t2/共1⫺t2兲,共26兲 t2⫽u2/共 ␳ ⫹u2兲,共27兲 dt⫽ ␳ 共 ␳ ⫹u2兲⫺3/2 du.共28兲 As uvaries from zero to infinity, tvaries from zero to unity. Finally we define a modified form of the two-dimensional integrals, Ix⫽exp共Dxt2兲共1⫺t2兲⫺1/2Ix ⬘,共29兲 and similarly Iyrelated to Iy ⬘and Izto Iz ⬘. Substitution of Eqs. 共15兲,共26兲–共29兲into Eq. 共14兲yields Eq. 共8兲, where PL共t2兲⫽2共 ␳ / ␲ 兲1/2IxIyIz,共30兲 and using yand zanalogs of Eq. 共23兲, it follows that X⫽Dx⫹Dy⫹Dz.共31兲 Substitution of Eq. 共30兲into Eq. 共9兲yields the working formula for the electron repulsion integral, 关 ␩ i ␩ j 储 ␩ k ␩ l兴⫽兺 ␣ ⫽1 N Ix共u ␣ 兲Iy共u ␣ 兲Iz #共u ␣ 兲,共32兲 where Iz #共u ␣ 兲⫽2共 ␳ / ␲ 兲1/2Iz共u ␣ 兲W ␣ .共33兲 C. Recursion formulas for the two-dimensional integrals We define gn,m x共x1,x2兲⬅共x1⫺xi兲n共x2⫺xk兲mexp关⫺A共x1⫺xA兲2 ⫺B共x2⫺xB兲2⫺u2共x1⫺x2兲2兴,共34兲 and the two-dimensional integral, Gn,m x⫽ 冕冕 gn,m x共x1,x2兲dx1dx2.共35兲 Gn,m xis related to Ix ⬘of Eq. 共16兲by placing the polynomial factors of the ␩ jand ␩ lprimitives on center Iand K. Thus, Gn,m xexp共⫺Gx兲⫽Ix ⬘共nxi⫹nxj,0,nxk⫹nxl,0,u兲.共36兲 Rys et al.6showed that 2069J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 Rys quadrature Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 m 2Gn,m⫺1 x⫹ 兵 B共xB⫺xk兲⫹u2共xi⫺xk兲 其 Gn,m x⫹u2Gn⫹1,m x ⫺共B⫹u2兲Gn,m⫹1 x⫽0, 共37兲 n 2Gn⫺1,m x⫹ 兵 A共xA⫺xi兲⫹u2共xk⫺xi兲 其 Gn,m x ⫺共A⫹u2兲Gn⫹1,m x⫹u2Gn,m⫹1 x⫽0. 共38兲 We make use of Eqs. 共18兲–共21兲and sum Eqs. 共37兲and 共38兲 to obtain a new recursion formula, BGn,m⫹1 x⫽⫺AGn⫹1,m x⫹n 2Gn⫺1,m x⫹m 2Gn,m⫺1 x ⫹ 兵 aj共xj⫺xi兲⫹al共xl⫺xk兲 其 Gn,m x.共39兲 Equation 共39兲is the key equation that underlies the new formulation. It provides a means to transfer angular momentum from electron 1 to electron 2. Indeed, Eq. 共39兲says that, in the 共n,m兲plane, point (n,m⫹1) which lies east of 共n,m兲 can be obtained from 共n,m兲and its north (n⫹1,m), west (n,m⫺1), and south (n⫺1,m) neighbors. Analogous equations hold for yand z. Equation 共39兲is depicted in Fig. 1. In their earlier work DRK made use of other recursion relations, also derived from Eqs. 共27兲and 共28兲. Of those we repeat only one that will be useful to us, Gn⫹1,0 x⫽nB10Gn⫺1,0 x⫹C00 xGn,0 x,共40兲 where C00 xand B10 are given by C00 x⫽共xA⫺xi兲⫹B共xB⫺xA兲 A⫹Bt2,共41兲 B10⫽1 2A⫺B 2A共A⫹B兲t2.共42兲 Finally we note the transfer equations that transfer angular momentum from center Ito center J, and from center Kto center L, Ix ⬘共nxi,nxj,m,0兲⫽Ix ⬘共nxi⫹1,nxj⫺1,m,0兲 ⫹共xi⫺xj兲Ix ⬘共nxi,nxj⫺1,m,0兲,共43兲 Ix ⬘共nxi,nxj,nxk,nxl兲 ⫽Ix ⬘共nxi,nxj,nxk⫹1,nxl⫺1兲 ⫹共xk⫺xl兲Ix ⬘共nxi,nxj,nxk,nxl⫺1兲.共44兲 Similar equations hold for Ix,Iy, and Iz. D. Recursion formulas for primitive integrals The recursion formulas given above are valid for the two-dimensional integrals that are used to build the primitive electron repulsion integrals as sums of products of the twodimensional integrals for Nroots, according to Eq. 共32兲. These recursion relations can be turned into recursion formulas applicable directly to the six-dimensional primitive integrals. Multiplication of Eq. 共39兲by exp(⫺Gx), then by exp(Dxt ␣ 2)(1⫺t ␣ 2)⫺1/2, and by Iyand Iz, followed by summation over roots and weights ␣ , leads to the following equation: B关n,0 储 m⫹1,0兴 ⫽⫺A关n⫹1,0 储 m,0兴⫹n 2关n⫺1,0 储 m,0兴⫹m 2关n,0 储 m⫺1,0兴 ⫹ 兵 aj共xj⫺xi兲⫹al共xl⫺xk兲 其 关n,0 储 m,0兴.共45兲 Recall that the bracket notation refers to six-dimensional integrals, in this case integrals over primitive Gaussian 共uncontracted兲functions. There is one more remark to make about Eq. 共45兲. In contrast to Eq. 共39兲that is valid for the twodimensional integral along the xdirection, and that offers no possibility of confusion, in Eq. 共45兲the direction for the increment or decrement 1of angular momentum is dictated by the coefficient of the fourth term in the sum. We really ought to use the notation 1x. However, to avoid complexity, we will keep the notation 1without subscript, and we will have to remember the implication due to the coefficient 兵 aj(xj⫺xi)⫹al(xl⫺xk) 其 . Equation 共45兲is closely related to Eq. 共30兲of LRL.10 The difference lies in our working with centers I,J,K, and Ldirectly, in contrast with these authors working with centers Aand B. Just as before, from Eqs. 共43兲and 共44兲substituted into Eq. 共29兲and then Eq. 共32兲it follows that 关ni,nj 储 m,0兴⫽关ni⫹1,nj⫺1 储 m,0兴⫹共xi⫺xj兲 ⫻关ni,nj⫺1 储 m,0兴共46兲 and 关ni,nj 储 nk,nl兴⫽关ni,nj 储 nk⫹1,nl⫺1兴⫹共xk⫺xl兲 ⫻关ni,nj 储 nk,nl⫺1兴.共47兲 Equations 共46兲and 共47兲are the ‘‘horizontal recursion relations’’ of Head-Gordon and Pople8applied in the xdirection. In these equations the variables ni,nj,nk, and nl have a collective meaning, for example, ni refers to the triplet 共nxi, nyi,nzi兲and similarly for nj,nk, and nl. The direction to which the transfer applied is implied by the factors (xi ⫺xj) and (xk⫺xl). Note that these relations are independent of the Gaussian exponents, something that is not the case for FIG. 1. Diagramatic representation of the recursion formula Eq. 共45兲.u stands for unscaled integral block, and sstands for scaled integral block. 2070 J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 M. Dupuis and A. Marquez Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 Eq. 共45兲. Thus Eqs. 共46兲and 共47兲may be applied outside the sums over primitive functions, while Eq. 共45兲must be applied inside the summation loop. Thus the computational strategy takes shape as follows: inside the loop over primitives, calculate the primitive integrals of the form 关n⫹m,0 储 0,0兴, use Eq. 共45兲to generate 关n,0 储 m,0兴, and sum up into the integrals over contracted functions (n,0 储 m,0). Outside the loop over primitives, apply the transfer equations 共46兲and 共47兲to form the final integrals (ni,nj 储 nk,nl). This algorithm is equivalent to the second LRL algorithm in Ref. 10. III. THE RYS QUADRATURE REVISITED Electron repulsion integrals over contracted basis functions can be written as a double sum over the ␩ i ␩ jand ␩ k ␩ l products of primitive Gaussians with the primitive integrals multiplied by the charge distribution coefficients CA ⫽caicaj and CB⫽cakcal 关the c’s are the contraction coefficients of Eq. 共4兲兴. Thus we write 共n,0 储 m,0兲⫽兺 B兺 A关n,0 储 m,0兴CBCA.共48兲 The notation is adapted from GHGP, the parentheses are the result of contraction of the square brackets for the bra 关n,0兩as well as for the ket 兩m,0兴. Let us assume for the moment that we have at our disposal the primitive integrals 关n,0 储 m,0兴,关n⫹1,0 储 m,0兴, 关n⫺1,0 储 m,0兴, and 关n,0 储 m⫺1,0兴, and that we would like to calculate not only the contracted integrals (n,0 储 m,0), (n ⫹1,0 储 m,0), (n⫺1,0 储 m,0), and (n,0 储 m⫺1,0), but also the integrals (n,0 储 m⫹1,0). For the first four groups of integrals we apply Eq. 共48兲directly. For the last group, (n,0 储 m ⫹1,0), we follow the procedure described below. This step is at the heart of the method. It corresponds to transferring angular momentum from electron 1 to electron 2. Substitution of Eq. 共45兲into Eq. 共48兲yields 共n,0 储 m⫹1,0兲 ⫽兺 B 1 BCB兺 A ⫺A关n⫹1,0 储 m,0兴CA ⫹n 2兺 B 1 BCB兺 A关n⫺1,0 储 m,0兴CA ⫹m 2兺 B 1 BCB兺 A关n,0 储 m⫺1,0兴CA ⫹共xj⫺xi兲兺 B 1 BCB兺 Aaj关n,0 储 m,0兴CA ⫹共xl⫺xk兲兺 B al BCB兺 A关n,0 储 m,0兴CA.共49兲 The strategy to be followed is explicitly depicted in Eq. 共49兲. Contraction of the primitive integrals 关n⫺1,0 储 m,0兴, 关n,0 储 m⫺1,0兴, and 关n,0 储 m,0兴over ␩ i ␩ j关terms 2, 3, and 5 in Eq. 共49兲兴 can be carried out inside the loop over A. The contraction of the primitive integrals 关n⫹1,0 储 m,0兴and 关n,0 储 m,0兴关terms 1 and 4 in Eq. 共49兲兴 can be carried out also inside the loop over A, but only after they have been scaled by ⫺Aand aj, respectively. Inside the loop over B共loop over ␩ k ␩ l兲contraction over Bis performed after scaling the first four partial sums by 1/Band the fifth sum by al/B共in the present case the scaling factor 1/Bcould be absorbed into CB, this is not true in general兲. The final contracted integrals (n,0 储 m⫹1,0) can be generated outside the loop over primitive functions by adding the five terms as in Eq. 共49兲after multiplication by the appropriate factors (xj⫺xi) and (xl ⫺xk), or their yand zanalogs, factors that are independent of the primitive functions. The situation is however a little more complicated than that. Indeed the various integral groups in Eq. 共45兲are themselves obtained through application of Eq. 共45兲for 0,1,...,n, and 0,1,...,m. Repeated application of Eq. 共45兲provides a means to compute any integral group 关i,0 储 k,0兴from 关n ⫹m,0 储 0,0兴with the condition i⫹k⭐n⫹m. This is truly the heart of the method. In what follows we derive the expressions for the 关i,0 储 k,0兴group directly in terms of the 关i,0 储 0,0兴 groups. These expressions are the key equations of the new algorithm and show precisely which scaled integrals need be calculated. Since there is no possibility of confusion we will use the notation 关n 储 m兴to represent the 关n,0 储 m,0兴integral group, and 关spdfg... 储 m兴to represent the ‘‘row vector’’ of integral groups 关s 储 m兴,关p 储 m兴,.... The derivation that follows initially deals with the ‘‘diagonal’’ integrals where the angular momentum exponents are accumulated on x,ory, or zonly. The generalization is not very much more complicated, it requires that one pays attention to the exact combination of angular momentum exponents that are being created. This point is illustrated towards the end of this section.The recursion relation Eq. 共45兲substituted into Eq. 共48兲 can be conveniently represented in a matrix notation. Let us denote Lthe total angular momentum for a group of four shells, we define two matrices of order L⫹1 denoted Iand R(L,1x). Istands for the identity matrix of order L⫹1 here 共not the basis function center兲. We define the quantity labeled cx, cx⫽ 兵 aj共xj⫺xi兲⫹a1共x1⫺xk兲 其 .共50兲 The R(L,1x) matrix is a tridiagonal matrix with ⫺Aon the lower diagonal, with i/2 on the upper diagonal where i ⫽1¯L⫹1 is the row index, and cxon the diagonal. The notation reflects that R(L,1x) is a matrix function of Land of one angular momentum direction, xin the present case. This is because Eq. 共45兲raises the power of the xangular monomial by 1, and equivalent equations can be written for yand z. The angular momentum direction 1xis linked to the diagonal term cx. For the time being we may ignore the subscript in cxfor ease of notation. At some point however, we will not be able to ignore it any longer, as discussed later. Thus for L⫽4 we have 2071J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 Rys quadrature Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 R共L,1兲⫽ 冤 c1 2000 ⫺Ac 2 200 0⫺Ac 3 20 00⫺Ac4 2 000⫺Ac 冥 .共51兲 Using this notation Eq. 共49兲takes the form 关spd¯L 储 m⫹1兴⫽1/B 兵 关spd¯L 储 m兴R共L,1兲 ⫹关spd¯L 储 m⫺1兴 其 m/2. 共52兲 Successive applications of Eq. 共52兲yield the expression relating 关spd¯L 储 m兴to 关spd¯L 储 s兴for any m,m ⫽1¯L/2. We have 关spd¯L 储 p兴⫽关spd¯L 储 s兴Rp, 关spd¯L 储 d兴⫽关spd¯L 储 s兴Rd, 共53兲 ] 关spd¯L 储 L/2兴⫽关spd¯L 储 s兴RL/2 . The notation in Eq. 共53兲indicates that each integral group 关s 储 d兴,关p 储 d兴,...,关L 储 d兴is, in a generic sense, equal to a linear combination of the integral groups 关s 储 s兴,关p 储 s兴,...,关L 储 s兴, with the coefficients being the matrix elements of the matrix Rm,m⫽1¯L/2, Rp⫽1 BR, Rd⫽1 B2 冉 R2⫹B 2I 冊 , 共54兲 Rf⫽1 B3 冉 R3⫹3B 2R 冊 , Rg⫽1 B4 冉 R4⫹3BR2⫹3B2 4I 冊 . ¯ In Eq. 共54兲we used the shortened notation Rfor the R(L,1) matrix of Eq. 共51兲. These equations provide a general recipe to express any group 关sp¯L 储 m兴in terms of 关sp¯L 储 s兴for any mand L. Symbolic algebra software can be conveniently used to obtain the matrices Rd,Rf,Rg,.... For an integral group 关ni,nj 储 nk,nl兴the target integral groups are 关max(ni,nj)¯ni⫹nj 储 max(nk,nl)兴up to 关max(ni,nj)¯ni ⫹nj 储 nk⫹nl兴. As an example a 关pp 储 pp兴group requires the four groups 关p 储 p兴,关p 储 d兴,关d 储 p兴,关d 储 d兴. The formula for 关p 储 p兴, in terms of 关s 储 s兴,关p 储 s兴, and 关d 储 s兴is found as the p,p elements of Rp, for 关d 储 p兴it is the element d,p of Rp, for 关p 储 d兴it is the element p,d of the matrix Rd, and for 关d 储 d兴it is the element d,d of Rdin accord with Eq. 共54兲. We have 关p 储 p兴⫽1 B 冉 1 2关s 储 s兴⫹c关p 储 s兴⫺A关d 储 s兴 冊 , 关d 储 p兴⫽1 B共关p 储 s兴⫹c关d 储 s兴⫺A关f 储 s兴兲, 共55兲 关p 储 d兴⫽1 B2 冉 c关s 储 s兴⫹ 再 c2⫺3 2A⫹B 2 冎 关p 储 s兴⫺2Ac关d 储 s兴 ⫹A2关f 储 s兴 冊 , 关d 储 d兴⫽1 B2 冉 1 2关s 储 s兴⫹2c关p 储 s兴⫹ 再 c2⫺5 2A⫹B 2 冎 关d 储 s兴 ⫺2Ac关f 储 s兴⫹A2关g 储 s兴 冊 . In this example, Eq. 共55兲says that in order to get the contracted (p 储 p) block we need the 关s 储 s兴block, the 关p 储 s兴block scaled with aj关recall the definition of cgiven earlier Eq. 共50兲兴 and the 关d 储 s兴block scaled with ⫺A. All these are computed inside the loop over the pairs of ␩ i ␩ jprimitives. Inside the loop over the ␩ k ␩ lprimitives, all the ␩ i ␩ j-contracted blocks get scaled with alor B-dependent factors. While these formulas may appear cumbersome, it must be remembered that the algorithm makes extensive reuse of the calculated basic integral blocks 关s 储 s兴,关p 储 s兴,...,关g 储 s兴and that scaling is relatively inexpensive, some of the scaling occurring outside the inner loops as seen in Eq. 共49兲. The target integrals blocks for the low angular momentum integrals up to (dd 储 dd) are given in the Appendix. The Appendix provides a global view of block scaling, i.e., which block must be computed, so that it is possible to generate all needed individual integrals. We do find there the terms in Eq. 共55兲for the (d 储 d) integrals. The complexity of the formulas grows quickly with the target integral blocks, i.e., with the angular momentum. As the effective degree of contraction Kdecreases with increasing angular momentum we may want to compromise between ease of programming and computational gains. In practice we can decide to write a computer code that covers all cases with basis functions up to d type. The recursion formula is represented diagrammatically as in Fig. 1. Before returning to the discussion of the efficiency of the proposed algorithm, we must recall that equations such as Eq. 共55兲were derived for the case when the angular momentum grows along one direction only. For example, Eq. 共55兲is correct when both p’s are identical, let us say the (s 储 s) and (px 储 s) integrals lead to the (s 储 px) integrals and the (s 储 s), (px 储 s), and (dxx 储 s) integrals lead to the (s 储 dxx) integral. To obtain cross integrals such as (s 储 dxy) we must raise the angular momentum in the xdirection and in the ydirection. We accomplish this by combining the recursion formulas in two directions, 2072 J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 M. Dupuis and A. Marquez Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 关spd¯L 储 mx⫹1兴⫽1/B 兵 关spd¯L 储 mx兴R共L,1x兲 ⫹关spd¯L 储 mx⫺1兴 其 共mx/2兲, 共56兲 关spd¯L 储 my⫹1兴⫽1/B 兵 关spd¯L 储 my兴R共L,1y兲 ⫹关spd¯L 储 my⫺1兴 其 共my/2兲. 共57兲 Recall that R(L,1x) and R(L,1y) differ through the diagonal terms, cxand cy, respectively. Application of Eq. 共56兲leads from (py 储 s) and (dxy 储 s)to(py 储 px), and from (s 储 s) and (px 储 s)to(s 储 px). Application of Eq. 共57兲then leads from (py 储 px) and (s 储 px) to the desired integral (s 储 dxy). In place of the second equation in Eq. 共53兲we can write bspd¯L 储 pp⬘c⫽关spd¯L 储 s兴Rpp⬘,共58兲 with Rpp⬘⫽1 B2RR⬘⬅1 B2R共L,1x兲R共L,1y兲.共59兲 Equation 共59兲is closely related to the second equation in Eq. 共54兲giving Rdfor like-angular momentum raising. For higher angular momenta we can define other needed transformation matrices Rdp ,Rpp⬘p⬙,Rdd⬘,Rdpp⬘, and so on. The list of scaled blocks for ‘‘diagonal’’ and ‘‘nondiagonal’’ integrals are identical, albeit the individual integral expressions are different. The approach is completely general and the formulas can be obtained with a symbolic algebra software. IV. DISCUSSION The key factor that makes the proposed algorithm computationally efficient lies in the very small number of terms in the ‘‘row vector’’ of integral groups 关spd¯L 储 s兴that must be computed in the inner loop. This number is equal to (L⫹1)(L⫹2)(L⫹3)/6. In the (pp 储 pp) the ‘‘row vector’’ of blocks 关spdfg 储 s兴contains 35 integrals, in the (dd 储 dd) case it contains 165 integrals. In contrast the number of individual target integrals grows as 兵 (l⫹1)关7(l⫹1)2 ⫺1兴/6 其 2where l⫽L/4. This number is 92⫽81 for (pp 储 pp) and 312⫽961 for (dd 储 dd). These latter numbers correspond to the work done inside all four loops over primitives in the LRL algorithm. To the extent that the additional work due to scaling and building of the target integrals from the scaled integrals stays at a reasonable level, this algorithm will be superior to LRL’s for moderately and highly contracted functions. Indeed this should the case. We note as well that the computational dependence for the 关spd¯L 储 s兴of (L ⫹1)(L⫹2)(L⫹3)/6 is smaller than the one reported by Ishida.12 Comparison of the present algorithm with the GHGP algorithm shows that both require computation of the same number of basic integrals inside the inner loop over primitives, mainly 关spd¯L 储 s兴integrals just as the MD method. For MD and GHGP these primitive integrals are centered on the exponent-weighted charge distribution center of the primitives, labeled A and B in Eqs. 共18兲and 共19兲. A step required for MD and GHGP involves the transferring of the target primitive integrals from the exponent-weighted primitive centers to the actual primitive centers. Although creating the target integrals 关spd¯ 储 spd¯兴from the basic integrals 关spd¯L 储 s兴is trivial in their approaches, the transfer from exponent-weighted centers to primitive centers requires additional work. In our approach this step is completely bypassed as the numerical quadrature yields primitive integrals directly centered on the primitive centers. An added benefit of the numerical quadrature is the extreme simplicity of computing the basic integrals, in contrast to the complexity of the tree search in the GHGP algorithm. Of course once we have calculated (n0 储 m0), then we can use the best possible strategy to generate the (ij 储 kl) integrals.16 It has become usual to express the performance of algorithms for the calculation of electron repulsion integrals by expressing the number of floating-point operations 共FLOP兲 needed to compute (pp 储 pp),(dd 储 dd),... integrals as a function of the degree of contraction K, assuming that all four shells p,d,... have the same degree of contraction. The expression is FLOP⫽xK4⫹yK2⫹z.共60兲 Equation 共60兲reflects the loop structure of the algorithm, x corresponds to the work done inside all four loops over primitives, ycorresponds to the work done outside the two inner loops 共 ␩ i ␩ j-contraction兲but inside the two outer loops 共 ␩ k ␩ l-contraction兲, and zcorresponds to the work done outside all four contraction loops. Table I gives x,y, and zfor the PH, GHGP, LRL, and DRK algorithms, along with estimates for the present algorithm. These estimates are based on the analysis that follows. The analysis becomes quickly complex and higher angular momentum cases are not considered here. A. Estimates for „ pp 储 pp …and „ dd 储 dd … To gain insight in the potential of the present approach, we consider the cases of (pp 储 pp) and (dd 储 dd) integrals. 1. Evaluation of the roots and weights of the Rys polynomials The evaluation of the Rys roots and weights is done inside the K4loop since the argument X关Eq. 共31兲兴 depends on the exponents of the primitive functions. For the (pp 储 pp) case we need 3 roots and 3 weights, and for the (dd 储 dd) case 5 roots and 5 weights. We can adopt an efficient TABLE I. FLOP counts for various shell quartets using different algorithms. PH Ref. 3 GHGP Ref. 9 LRL Ref. 10 DRK Ref. 5 DM This worka (pp 储 pp)x220 300 691 904 440 y2300 610 30 30 140 z4000 680 330 0 650 (dd 储 dd)x2450 10 255 14 660 2575 y25 800 30 30 2550 z28 900 11 300 0 18 800 (sp sp 储 sp sp)x220 450 753 1056 y2300 1300 30 30 z4000 1700 800 0 aEstimates. See text for explanations. 2073J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 Rys quadrature Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 approach17 involving lookup tables and interpolations. A 3-point interpolation, with 3 multiplications and 2 additions, can yield roots and weights of sufficient accuracy for a total of 6*5⫽30 FLOPs for (pp 储 pp) and 10*5⫽50FLOPs for (dd 储 dd). No additional work is necessary for sp shells 共Table II兲. 2. Evaluation of the two-dimensional integrals I x (nx ,0, t ␣ ) , I y (ny ,0, t ␣ ) , and I z (nz ,0, t ␣ ) This step is done inside the K4loop. Gx,Gy, and Gzare calculated by means of the 2-point recursion formula of Eq. 共40兲for each root t ␣ . There are 3 roots for the (pp 储 pp) and (spsp 储 spsp) cases, and 5 roots for the (dd 储 dd) case. For (pp 储 pp) integrals the indices nx,ny, and nz run up to 4. The cost of this step is 3*(7 multiplicationsplus6additions) ⫽39FLOPs. For the (dd 储 dd) case we have 5*(15 multiplications plus 14 additions)⫽145FLOPs. About 16 FLOPs per root are needed to create the coefficients in Eqs. 共41兲and 共42兲, about 48 FLOPs for the (pp 储 pp) case and 80 FLOPs for the (dd 储 dd) case. 3. Evaluation of the primitive integrals [n 储 s] For (pp 储 pp) we need 关nx ny nz 储 0兴for nx,ny,nz⫽1to 4, i.e., 关g 储 s兴integrals. There are 35 such integrals. This K4 step is carried out according to Eq. 共32兲in the most general case. In the case at hand, we can take advantage of the fact that whenever nx,orny,ornz is equal to zero, the corresponding Gx,orGy,orGzfactor is equal to 1.0. We do not have to multiply by that factor. The cost per integral is essentially one multiplication to form the product Gx*Gy*Gz for each root plus (nroots-1) additions for the summation over roots 共Table II兲. For (pp 储 pp) this amounts to (35-5)*3⫹35*2⫽160 FLOPs. For (dd 储 dd) it amounts to (165-9)*5⫹165*4⫽1440 FLOPs. 4. Evaluation of the scaled integrals Tables III and IV illustrate the block scaling, at the K4 level and at the K2level, respectively. Scaling has to be performed to generate the complete set of terms needed for the (pp 储 pp) integrals. The tables are based on the formulas given in the Appendix. The work amounts to 79 共multiply⫹add兲inside the K4loop and 53 共multiply⫹add兲 inside the K2loop. For sp shells these numbers are 90 and 59, respectively. For (dd 储 dd) integrals, a similar analysis shows that the K4scaling can be realized with 856 共multiply⫹add兲, while the K2scaling amounts to 1258 共multiply⫹add兲. 5. K 0 work and horizontal recursion formulas Once the scaled integrals have been calculated, the final integrals (n,0 储 m,0) have to be assembled before the horizontal recursion formulas can be applied. Inspection of the expressions given in the Appendix and of the actual expres- TABLE II. Algorithm for constructing the 关g 储 s兴primitive integrals for the (pp 储 pp) integrals. 1兲s21兲xxxx ⫽1)*G(xxxx) 2兲x⫽1)*G(x)22兲xxx y ⫽3)*G(xxx) 3兲y⫽1)*G(y)23兲xxx z ⫽4)*G(xxx) 4兲z⫽G(z)24兲xx yy ⫽5)*G(yy) 5兲xx ⫽1)*G(xx)25兲xx y z ⫽9)*G(x) 6兲xy ⫽2)*G(y)26兲xx zz ⫽10)*G(xx) 7兲xz ⫽4)*G(x)27兲x yyy ⫽2)*G(yyy) 8兲yy ⫽1)*G(yy)28兲xyyz⫽18)*G(x) 9兲yz⫽4)*G(y)29兲xyzz ⫽19)*G(x) 10兲zz ⫽G(zz)30兲x zzz ⫽20*G(x) 11兲xxx ⫽1)*G(xxx)31兲yyyy ⫽1)*G(yyyy) 12兲xx y ⫽3)*G(xx)32兲yyy z ⫽4)*G(yyy) 13兲xx z ⫽4)*G(xx)33兲yy zz ⫽10)*G(yy) 14兲xyy ⫽2)*G(yy)34兲y zzz ⫽20)*G(y) 15兲xyz⫽9)*G(x)35兲zzzz ⫽G(zzzz) 16兲xzz ⫽10)*G(x) 17兲yyy ⫽1)*G(yyy) 18兲yy z ⫽4)*G(yy) 19兲yzz ⫽10)*G(y) 20兲zzz ⫽G(zzz) TABLE III. K4scaling for the (pp 储 pp) and (dd 储 dd) integrals.a (pp 储 pp)(dd 储 dd) K4K2K4K2 ajspd spdfg Aspdf spdfgh aj 2spd spdfg Aajpdf pdfgh A2dfg dfghi aj 3pdfg Aaj 2dfgh A2ajfghi A3ghij aj 4dfg Aaj 3fgh A2aj 2ghi A3ajhij A4ijk aEntry nmeans that a primitive block 关n 储 0兴gets scaled by the coefficient in column 1. For example in the (pp 储 pp) case, the first entry in the table indicates that the 关s 储 0兴,关p 储 0兴, and 关d 储 0兴blocks get scaled with the aj coefficients to form three scaled blocks. Entries in italics are needed for sp shells. 2074 J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 M. Dupuis and A. Marquez Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17 sions in the (pp 储 pp) case suggest that in average each of the (pd 储 pd) integrals needed for a (pp 储 pp) block is the sum of about 4 unscaled or scaled contributions, and each of the (dfg 储 dfg) integrals needed for a (dd 储 dd) block is the sum of about 8 contributions, for a total of about 320 FLOPs and 7500 FLOPs, respectively. Finally the horizontal recursion transformation from (n,0 储 m,0) to (i,j 储 kl) can be performed with the most efficient algorithm of LRL.10 Together these two steps yield the K0coefficient in Table I. B. Discussion It is clear from Table I that DRK remains the method of choice for uncontracted functions. This is not really unexpected. The other methods are more efficient for contracted functions. They were designed to take advantage of the concept of charge distribution to varying extents. This can be seen in Table I. MD’s method falls in that category as well. Our goal was to devise an algorithm based on the Rys quadrature that would deal efficiently with contracted functions. The critical factor for success resides in the reduction in the number of fundamental integrals that need to be computed inside the K4loop. In the spirit of minimizing this number, it is clear that PH’s method remains the most successful 共many are equal to zero due to using local axis frames兲, albeit at the cost of much higher K2and K0coeffi- cients. It is clear also from Table I that the present algorithm is superior to LRL’s method for all the cases. The reason lies in this minimum number of integrals computed in the K4 loop. The trick about transferring angular momentum between electrons, a technique that we borrowed from MD and GHGP, is responsible for the gain, since both the LRL method and the present method use the same numerical quadrature. Table I shows that the K4coefficient in the present method is close to the same coefficient in the GHGP’s method, while the K2and K0coefficients are much smaller than those in GHGP. The similar magnitude of the K4coef- ficients between the two methods comes in part from the fact that both methods compute the same number of fundamental integrals. The much-reduced K2and K0coefficients in the present method is a direct reflection of the advantage of using a numerical quadrature and working directly with the basis function centers, rather than the charge distribution centers. The ‘‘scaling’’ step is much simpler because of that. That the new method does not do quite as well as GHGP for (pp 储 pp) integrals shows that there is some inherent overhead in the numerical quadrature when dealing with low angular momentum functions. We do not seem to have succeeded in overcoming this point. However the present algorithm appears to be largely superior to GHGP’s for (dd 储 dd) integrals. We expect this to hold true for higher angular momentum functions. We expect the present method to be competitive for (spsp 储 spsp) integrals, although we have not fully analyzed the effect of the different contraction coefficients for the sfunction and the pfunctions of an sp shell. V. SUMMARY We have outlined a novel approach to the calculation of electron repulsion integrals over Gaussian basis functions. It takes advantage of the best features of several approaches: the angular momentum transfer of McMurchie–Davidson, the scaling technique of Gill–Head–Gordon–Pople, the horizontal transfer equations proposed by Head–Gordon–Pople, and the Rys numerical integration of Dupuis–Rys–King. The resulting algorithm is anticipated to be quite competitive for contracted Gaussian functions. The derivation and implementation are simplified through the use of symbolic algebraic software. The method generalizes to high angular momentum integrals. Implementation is in progress and performance data will be presented in future reports. The extension of the method to integral 共first- and second-兲de- rivatives with respect to nuclear coordinates should also benefit from the advances devised in the present work. ACKNOWLEDGMENTS We dedicate this work to the memory of John Rys. We thank Dr. A. Komornicki for stimulating discussions and his careful reading of the paper. We thank also Dr. Graham Fletcher for useful discussions. This work was supported in part by the U.S. Department of Energy’s Office of Biological and Environmental Research, and by the Office of Basic Energy Sciences, Chemical Physics Program. The Pacific Northwest National Laboratory is a multiprogram national laboratory operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract No. DE-AC06- 76RLO-1830. TABLE IV. K2scaling for the (pp 储 pp) and (dd 储 dd) integrals.a (pp 储 pp)(dd 储 dd) K4K2K4K2 alspd spdf (aj)alspd spdfg (A)alpdf pdfgh (aj 2)alpdfg (Aaj)aldfgh (A2)alfghi (aj 3)aldfg (Aaj 2)alfgh (A2aj)alghi (A3)alhij Ba lpdfg (aj)Ba ldfg (A)Ba lfgh al 2spd spdfg (aj)al 2pdfg (A)al 2dfgh (aj 2)al 2dfg (Aaj)al 2fgh (A2)al 2ghi Ba l 2dfg al 3pdfg (aj)al 3dfg (A)al 3fgh Bspd spdfg (aj)Bpdfg (A)Bdfgh (aj 2)Bdfg (A2)Bghi B2dfg aSee explanations in Table III. 2075J. Chem. Phys., Vol. 114, No. 5, 1 February 2001 Rys quadrature Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 150.214.230.47 On: Mon, 14 Nov 2016 10:01:17