Volume 184, number 1,2,3
CHEMICAL PHYSICS LETTERS
20 September 1991
Two-electron integrals and integral derivatives revisited Itai Panas Department of Inorganic Chemistry, University ofGiiteborg, Vasaparken, S-41 I 24 Gothenburg, Sweden, and Chalmers University of Technology, S-41 2 96 Gothenburg, Sweden
Received 4 June I99 1
A new vertical recurrence relation (NVRR) for electron-repulsion integrals (ERIs) is introduced. NVRR connects the Cartesian Gaussian ERI with the auxiliary functions, introduced by McMurchie and Davidson. Aspects of NVRR were used by McMurchie and Davidson in their fundamental work of 1978,by Gill, Head-Gordon and Pople in their scaling method ( 1990) and by Hamilton and Schaefer in their algebra for moving angular momenta between electrons (1990). Here, NVRR is employed to obtain an efficient formulation for the ERI gradients. Also, a straightforward procedure for constructing two-electron integrak over general spherical tensor operators is outlined.
1. Introduction The interest in efficient algorithms for generating electron-repulsion integrals ( ERIs) has increased gradually since the time when a direct scheme was first introduced for solving the Hartree-Fock problem [ 1 ] and peaked with the introduction of direct second-order many-body perturbation theory MP2 [ 2,3]. One reason for this is that much larger calculations can be performed with the direct technique than with the conventional one, as the integral storage bottleneck is replaced by a CPU time bottleneck. An important incentive in this endeavour is that surprisingly good results for a wide variety of systems can be obtained with MP2 [ 41. The success of the contracted configuration-interaction externally method [ 51 for describing electron correlation can, to some extent, be attributed to this fact. In that method, the contraction coefficients for the contracted configurations are determined from the MP 1 wavefunction. A major accomplishment of the direct schemes is that ab initio calculations can formally be performed on large molecular systems. Crucial in these cases is to be able to benefit from asymptotic simplifications of the formalism. Generally, this is done by taking the asymptotic value of the incomplete gamma-function, which obviously still depends on all four prim-
86
itive basis-function indices. In ref. [ 6-8 ] we present an alternative route, which exploits the separability of the electronic coordinates at the asymptotic region. Thus, we are able to avoid the construction of primitive and contracted integrals in this case, to proceed directly to the construction of the Fock matrix. We believe this to be the most efficient technique to obtain the long-range contribution to the Fock matrix. The focus of the present study is set at investigating the short-range two-electron integrals, corresponding to the region where the electronic coordinates cannot be separated. A new vertical recurrence relation (NVRR) is deduced in section 2.1. NVRR is the simplest recursive connection between the ERI and the McMurchie and Davidson (MD) auxiliary functions [ 91. The NVRR is similar in spirit to the recurrence relation (RR) described by Obara and Saika (OS) [ lo], although their RR connects the ERI with the incomplete gamma-function. An RR, which is equivalent to the NVRR but expressed somewhat differently, is central in the Gill, Head-Gordon and Pople (GHP ) scaling method [ 11,121. This method is a synthesis of the MD algorithm and the Head-Gordon and Pople ( HGP ) horizontal recurrence relation (HRR ) scheme [ 131. By scaling the MD auxiliary functions with the exponents of the primitive Gaussians, GHP
0009-2614/91/$ 03.50 0 1991 Elsevier Science Publishers B.V. All rights reserved.
Volume184,numberI ,2,3
CHEMICAL PHYSICS LETTERS
are able to contract over the primitive basis functions before performing their RR. This is the most efficient way to construct contracted ERIs to date
20 September 1991
Cartesian Gaussian functions ( CGFs), A CGF is defined as
L1-21. A direct consequence of the understanding of the construction of ERIs and ERI derivatives in terms of the NVRR, is that it enables us to obtain readily an expression for Cartesian Gaussian-type two-electron integrals with a general spherical tensor operator. This is shown in section 2.2. Some special cases were presented in refs. [ 9,14 1 but the present work formally generalizes and clarifies their deductions. In section 3, a modification of the Hamilton and Schaefer (HS) approach [ 15 ] to the evaluation of integral derivatives is demonstrated. It is shown that although some aspects of the NVRR were used by HS, it was not exploited fully. Section 4 comments on results obtained and relationships between the NVRR and the GHP methods. Section 5 is a summary.
Xexp[ -rc(r-K)2]
,
(5)
where r= (x, y, z), k= (k,, ky, k,) and K= (K,, KY, K,) designate the position of the electron described, the angular momentum and the origin of the CGF, respectively. All positions are taken relative to a laboratory-fmed coordinate system. There are two conventional ways to proceed from ( 1). One is to expand the integrand in Rys poly nomials and perform a Gaussian quadrature [ 14,161, and the other is to expand the integral in incomplete gamma-functions and functions related to them (see, e.g., refs. [ 9,10,14] and references therein). In what follows, we proceed in the spirit of the latter. Section 2.1 introduces a new vertical recurrence relation (NVRR), while section 2.2 demonstrates how to benefit from the NVRR in the case of two-electron integrals with spherical tensor operators.
2. Two-electron integral evaluation 2. I. A new vertical recurrence relation The common denominator for all conventional ERI evaluation schemes using Gaussian basis functions is having to deal with an integral of the form
s
(1)
W=K(P---(1)2)
(2)
K&Y_
(3)
SA(s2)exp( - Ws’) cls,
p+4
where S,(s2) is a polynomial of degree ,l in s2, q=y+6 and (Y,8, y and 6 are the exponential factors belonging to a general quadruple of basis functions. In the case of the origin positions of the basis functions being A, B, C and D, respectively, P and Q become
p=a+b,
p,_ d,+BB, I-
and
a+B
whereiisx,yorz.
The standard notation to represent primitive and contracted integrals will be used. Thus, by [ ab 1cd] is meant a two-electron integral, consisting of a quadruple of primitive CGFs that belong to four contractions with angular momenta a, b, c and d. The contracted counterpart of [ ab 1cd] is (ub 1cd). Thus, it is common to exclude the primitive and contracted basis-function indices, since these are deemed irrelevant for the discussion, We proceed by differentiating the primitive integral to obtain
(n “.) aA,
+
as.
[UblCdl=2a[(u+l,)bIcd]
-ai[(a--li)blcd]+ZP[a(b+li)lcd]
-4[a(b-l,)lcdl,
Qi=y, (4)
In what follows, we will work with un-normalized
(6)
where by 1; we mean (6,, S;,,,~3,). Writing r,i-Bi= (r,i-Ai) + (Ai-Bi) and introducing it in the third term of (6) we have 87
Volume 184, number I ,2,3
CHEMICAL PHYSICS LETTERS
erators in (8) and ( 10) commute with p. We also have
(ih“,> + z
[abjcd] =2(cr+/3) [ (a+l;)alcd]
-a,[ (U-li)blCGf]
t2fl(Ai_Bj)
20 September 1991
[UblCd]
-b,[a(b--1,)WI
(7)
Or
(13)
[(a+li)blcdl This is obtained from the definition of P, since
+ai[ (a-l,)blcd]
tbJu(b-li)
-2/q/$-B;)
($ +&)m-Q)
[abled]
[cd] . 1
=(s+~)&A-Q,
(8)
=&fV-Q)
By writing
.
I
a
_a+2
an,=aA;
(9)
aJ3j’
we conclude that a general ERI [ablcd] can be expressed in terms of a linear combination of expressions,
The right-hand side of ( 13) is related to the MD auxiliary function by R(t
, u>
v)=
m-q)“2Pq 2rP2
X(&)&&&I]O]
(10)
(14)
3
(15a)
and as a generalization of the relations for two-electron field- and field-gradient integrals obtained by MD, it is easy to show that
where
x ev( -&d dTl dTz, and thus,
We introduce the definition
I
2n:5/2
[WI = tp+qj,,2pq 5o exC--~2) b.
(lib)
The integral in (1 lb) is the zeroth-order imcomplete gamma-function and p= exp - SP
(A-B)2
Xexp - -$ (
(C-0)2).
(
)
(12)
From (12), it is apparent that the differential op88
(15b)
UW
IO17
(15c)
in order to simplify the discussions to follow. We conclude that since a/an, commutes with each of the differential operators, with p and with the coefficients of (B), it also commutes with the superposition of operators corresponding to the original ERI. Hence, we can write the ERI as this superposition of operators acting on the expression
Z$IOOIOOI=[OOl~(L
dip
Jiz)
1001
3
(16)
1
and obtain
Thus, eq. (8) can be written as
[(a+li)bld
+bi[a(b-li)]cd]].
(18)
This is the new vertical recurrence relation (NVRR ). The NVRR at center B is obtained by interchanging A with B and LYwith 8. In order to obtain the NVRR for the angular momentum of center C, we proceed in the same manner as for A, noting that
(&+~)[~lc~l =-(&+&)WWl
>
I19)
obtained from the well-known translational invariance of ERIs. Thus, we obtain
[abl(c+li)dl
(C-li)d]-2d(C,-Di)
+di[UbjC(d-li)]}
.
VRR of OS as applied by HGP inside the innermost loop over primitives. GHP observed that scaling the MD functions with the exponents emerging in the coefficients of the MD analog to the NVRR, and contracting before performing the corresponding recursion relation, made this technique more efficient than previous ones. This is due to the fact that it is preferable in many cases to work with the representation of the integrals rather than with the integrals themselves [ 121, because this reduces redundancies in the data handling. 2.2. Integrals with spherical tensor operators
+Uj[ (U-l~)blCd]-2fi(Ai-Bi)[UblCd]
+ci[abl
20September1991
CHEMICAL PHYSICSLETTERS
Volume184,number1,2,3
A deeper understanding of the relationship between the Coulombic interaction operator and its representation in terms of ERIs emerges. It is apparent that any linear combination of differential operators of the form found in ( 10) commutes with the NVRR. This directly implies that this operator will penetrate the integral expression and differentiate the Coulombic interaction operator in the same manner as found in ( 15b). Thus, we conclude: the fact that the electron-repulsion integrals representation is a projection of the Coulombic interaction operator on the P’s and Q’s of the basis set, is also valid for the derivatives of this operator if the corresponding derivatives of the integral are defined as suggested above. In particular, we have that the result of the wellknown formula for irregular solid harmonics as given by Hobson [ 171,
[Ublcd] (20)
To obtain the NVRR on D we interchange C with D and y with 6. It is evident that the contracted integrals (uO]cO)+( (u+b)O] (c+d)O) needed for the HRR described in ref. [ 13 1, can be constructed by applying the NVRR to the MD auxiliary functions. We also observe that only three of the four terms in eq. ( 18) survive in case of the HRR. Furthermore, the NVRR allows for the constructions of angular momenta associated with electrons 1 and 2, independently. This should be compared to the live-term
=
Pj”‘(cose~2)exp(+iImIq,z) ,
r::’
(21)
where a~~~~+-~‘m’($
&i-$-)-I,
(22a)
when replacing rb’ on the left-hand side of (2 1) by [ab )cd] and defining a differential operator analogous to (22a),
89
Volume184,number 1,2,3
20 September 1991
CHEMICAL PHYSICS LETTERS
becomes
~abtai++)I4
=
X(+i)l”‘-‘[abID(e,
5;(I;‘)
-C!P(Ai-Bi) [abId] +bi[U(b-li)
The fact that the coefficients in (23b) are exponentindependent obviously implies that the sum can be taken at the contracted ERI level,
X(&i)‘“‘-‘(abID(e,
By allowing for the differential operators on the righthand side of (25a) to act and employing the same technique as was done between (6) and (7), we obtain
Jml-e,Z-]ml)]cd), (23~)
where each term can be solved by the HRR and the NVRR.
3. ERI gradients As observed by OS, the expense for evaluating an ERI derivative, & [ablcd]=2a[(a_tl,)b(cd] ,
+Ci[Ubt (c-l,)d]
tdi[UblC(d-li)]
-(28(Ai-B,)+2S(Ci_Di))[UbtCd]
-2 (Y+d)[&I (C+li)dl> .
Wb)
This is the relationship obtained by HS for shifting angular momentum between the electrons. In the following, we assume the overall strategy of HS but proceed from ( 18 ) in conjunction with (24) to obtain
tUi[ (U_l,)blCd]-2~(A;-Bi)[UblCd]
1
(24)
using a recurrence relation, is related to the evaluation of the first term in (24), since the second is an intermediary of the former. The strategy for ERI derivatives that emerges from this observation, assumed by Hamilton and Schaefer (HS) [ 151, is to try and minimize the number of integrals such as the first term in (24). In order to show the close relationship of the NVRR to the technique of HS, we begin by reproducing two of their results. Introducing ( 19) into (8), we obtain
tb~[a(b-l,)]cd]}-u,[(a-li)alcd]
,
(26a)
.
(26b)
and by analogy,
+Cl[Ubl
(C-li)d]-26(Cj_Di)
tdi[~IC(d-Ii)])-C,[~t
[UblCd]
(C-lt)d]
Eq. (26b) is equivalent to the corresponding derivative given by HS. From the definition of [ pb I D( Six, dip Siz) [cd], we readily obtain [ahi tCd=[ablD(&, sly,hz) 1~41
- [aibWl ,
90
s
lml-e,I-Im])lcd]. Wb)
+[(a--lj)blcd]
1Wa)
Id]
VW
Volume184, number 1,2,3
- [ablc,4
CHEMICAL PHYSICS LETTERS
,
Wb)
which can be taken at the contracted level,
(28a)
-(UbICid)
.
(28b)
Attractive features are: (i) All the integrals involved except for (ublD( S,,, &,, 6,,) Icd) contain lower than or the same angular momenta as (abl cd). (ii) The construction of (ub I D( 6ix, Siy, &) Icd) is parallel to (ublcd) in that the MD auxiliary functions building (ub ID( d,, ~5,~&) Icd) are shifted (Six, d,, d,,), as compared to (ub (cd). This suggests that it is cheaper to construct (ubl D(&., 8, S,,) Icd) than ((U+li)blCd). The latter is constructed in the HS scheme. (iii) The handling of constituents resulting from derivatives with respect to center B is avoided, as opposed to HS or any other conventional integral derivative approach.
20 September 1991
traction loops, since redundant work is reduced. A first step is to make the auxiliary functions angular dependent. This is done by turning to the MD auxiliary functions, which are angle-dependent constructions from the incomplete gamma-function. The second step is to scale these auxiliary functions with the exponents to obtain a VRR which can be applied at the contracted level. This can be done for each of the electrons, resulting in a choice whether to perform the scaling procedure on both of them (“early contraction”), on one (“late contraction”) or not at all. The NVRR (section 2.1) displays the features exploited by GHP and simplifies their concepts in that no reference to the transformation between Hermite and Cartesian Gaussian functions needs to be considered. The connection lies in the translational invariance of ERIs as used in ( 19 ). This is easily understood by virtue of eq. ( 14). Hence, we have implicitly exploited the translational invariance of the- integral [KIZ] in that
and for which the equivalent of the NVRR thus becomes
4. Comments A natural ambition, when trying to construct efficient algorithms for contracted ERIs, is to seek to reduce the amount of work inside the innermost loops over primitives. Thus, the philosophy of the HGP and GHP methods is to form the true contracted ERI from a few contracted auxiliary functions outside the loops over primitives_ In the HGP method, a vertical recurrence relation (VRR) is employed inside the contraction to obtain a class of contracted ERIs with angular momenta collected to two centers. The true contracted ERI is obtained by applying their horizontal recursion-relation technique to the contracted set. Being a modification of the OS method, the HGP approach constructs the contracted set directly from the incomplete gamma-functions. The success of the GHP technique lies in realizing that moving information from the VRR to the auxiliary functions reduces the amount of work inside the innermost con-
tkt[K-iiILlj
(30a)
and
+rJKlz-l;]j.
(job)
Thus, obviously, [KI Z] can be constructed from the MD auxiliary functions. But rather than forming these intermediary Cartesian Gaussian integrals, section 2.1 constructs the four-center quantities directly. As for ERI derivatives, we believe that the NVRR has brought new insights concerning the nature of the two-electron integral in general (section 2.2 ), and ideas as to how to improve on gradient evaluation algorithms in particular (section 3).
91
Volume 184, number 1,2,3
CHEMICAL PHYSICS LETTERS
5. Summary The present work has addressed the two-electron integral evaluation from a somewhat new angle: (i) A new vertical recurrence relation for the evaluation of electron-repulsion integrals was deduced. (ii) The relationship between the Coulombic interaction operator and its electron-repulsion integrals representation was discussed. (iii) A general expression for the two-electron integral with spherical tensor operators was presented. (iv) An improvement of the approach, assumed by Hamilton and Schaefer to simplify the evaluation of ERI derivatives, was obtained.
References [ 11J. Almliif, K. Faegri and K. Korsell, J. Comput. Chem. 3 (198,)
385.
[ 21 M. Head-Gordon, J.A. Pople and MJ. Frisch, Chem. Phys. Letters 153 ( 1988) 503.
[ 31 S. Sambaand J. Almliif, Chem. Phys. Letters 154 (1989) 83.
92
20 September 1991
[4] N.C. Handy, private communication. [5] P.E.M. Siegbahn, Intern. J. Quantum Chem. 23 (1983) 1869. [6] 1. Panas, J. Almhiif and M.W.F. Feyereisen, Intern. J. Quantum Chem., in press. 1711. Panas and J. AlmhGf, Intern. J. Quantum Chem., submitted for publication. [S] 1. Panas, subsequent paper. [P] L.E. McMurchie and E.R. Davidson, J. Comput. Phys. 26 (1978) 218. [lo] S. Obam and A. Saika, J. Chem. Phys. 84 ( 1986) 3963. [I 11P.M.W. Gill, M. Head-Gordon and J.A. Pople, Intern. J. Quantum Chem. Symp. 23 (1989) 503. [ 12 ] P.M.W. Gill, M. Head-Gordon and J.A. Pople, I. Phys. Chem. 94 (1990) 5564. [ I3 ] M. Head-Gordon and J.A. Pople, J. Chem. Phys. 89 ( 1988) 5777. [ 141V.R. Saunders, in: Methods in computational molecular physics, eds. G.H.F. Diercksen and S. Wilson (Reidel, Dordrecht, 1983). [ 15] T-P. Hamilton and H.F. Schaefer, Chem. Phys. 150 ( 1991) 163. [16]M.Dupuis,J.RysandH.F.King, J.Chem.Phys.65 (1976) 111. [ 171 E.W. Hobson, The theoty of spherical and ellipsoidal harmonics (Cambridge Univ. Press, Cam bridge, 193 1) .