Three-center expansion of electron repulsion integrals with linear combination of atomic electron distributions

7 July 1995

CHEMICAL PHYSICS LETTERS

ELSEVIER

Chemical Physics Letters 240 (1995) 578-584

Three-center expansion of electron repulsion integrals with linear combination of atomic electron distributions Seiichiro Ten-no, Suehiro Iwata Instihtte for Molecular Science, Myodaiji, Okazaki, Aichi 444, Japan

Received 9 January 1995; in final form 2 April 1995

Abstract A new systematic integrals is proposed density is expanded

way of constructing auxiliary basis functions for approximating the evaluation of electron repulsion and applied to SCF and MCSCF wavefunction calculations. In the approximation, the one-electron

in terms of a linear combination of atomic electron distributions (LCAD), and the four-center two-electron repulsion integrals are reduced to the three- and two-center quantities. This results in a high-accuracy approximation as well as a large reduction in disk storage and input/output requirement, proportional to N3 rather than N4, N being the number of basis functions. Numerical results indicate that the error from the present approximation decreases as the size of molecular basis functions increases and that the LCAD version of MCSCF calculations requires only a fractional amount of the CPU time required in the conventional

procedure without loss of accuracy.

1. Introduction Since the first use of Cartesian Gaussian-type functions as basis functions in ab initio methods, much effort has been devoted to efficient algorithms for evaluating two-electron repulsion integrals (ERIs) over primitive basis functions as well as over contracted functions [l-4]. The size of molecule we can deal with, however, is still strongly restricted by the N4 scaling of the total number of ERIs. To reduce the number of integrals to be calculated and transformed, several approximate methods have been proposed preserving the traditional linear combination of atomic orbital (LCAO) picture. One of the most important contributions is the pseudospectral method developed by Friesner and co-workers [5]; the method scales as N3 both in integral evaluation and Fock matrix assembly. Although the method has proven 0009-2614/95/$09.50

capable of achieving the high accuracy required in quantum chemical calculations, one is forced to handle additional numerical grids and dealiasing basis functions which may alter the numerical accuracy. Furthermore, the approximation is, in principle, valid only under the condition that the molecular basis set is nearly complete. In quantum chemical methods, it is frequently important to estimate energy differences with accuracy rather than the absolute quantities of reactant and product with a given set of basis functions. Another line of such approaches is on the basis of employing a set of new auxiliary functions to approximate the ERIs. A number of methods have been proposed to improve the quality of approximations [6,7]. Recently, a three-center expansion of the fourcenter ERIs in ab initio methods, the so-called V approximation, was developed by Vahtras, Almlof

0 1995 Elsevier Science B.V. All rights reserved

SD1 0009-2614(95)00564-l

S. Ten-no, S. Iwata / Chemical Physics Letters 240 (I 995) 578-584

and Feyereisen [8] and extended to coupled cluster theory by Rendell and Lee [9]. In their approximations, the orbital products are expanded in a set of auxiliary functions, and the coefficients are determined such that each self-Coulombic energy is reproduced. This is analogous to the idea used only for the direct Coulombic part in the LCAO local density approximation [lo] and has been accepted as an accurate fit. The asymptotic behavior of the numerical length of such expansions required to ensure accuracy would become proportional to the number of basis functions N rather than N 2 for large N, as shown by Beebe and Linderberg [ll] in the Cholesky decomposition of ERIs. The expansion, therefore, leads to an N3 behavior in the number of ERIs needed and an N4 behavior in the integral transformation. The expansion technique would enable us to remove the main bottlenecks concerning ERIs in the electronic structure calculation of large molecules, if one could obtain appropriate functions to mimic the one-electron densities in a systematic way. But, so far, an optimal choice for expansion functions corresponding to a various kind of molecular basis functions has not been well-addressed. In this Letter, we propose a method in which the auxiliary functions are simply the atomic electron distributions (ADS) consisting of the basis functions on each nucleus. The one-electron densities and the orbital products are approximated by a linear combination of atomic electron distributions. Then, the four-center ERIs, which include no AD, are approximately reduced to three- and two-center ERIs. The total number of corresponding ERIs scales again MN2, where M is the number of ADS which increases in proportion to the size of polyatomic molecules, and so to the number of the basis functions. The main feature of the present approximation is twofold. (1) One can construct a set of expansion functions in a systematic way for any molecular basis functions maintaining the N3 scaling of integral evaluation and disk storage, and correspondingly the N4 scaling of MO integral transformations. (2) The present approximation does not cause any errors in the electron density near the atomic nuclei, where the most accuracy is required to reproduce the total energy. In order to demonstrate the efficiency of the proposi-

579

tion, we briefly illustrate the approximation and present several numerical and timing results applied to SCF and MCSCF calculations in the following sections.

2. Method Within a Roothaan type expansion, one-electron densities of molecular orbitals are expressed by a linear combination of orbital products { x,( r)xb( t)} (a=1 ,..., N; b=l,..., a) for a given set of basis functions {x,(t)} (a = 1,. . . , N) which are in most cases located on each nucleus. Beebe and Linderberg termed these products distributions [ll]. In the three-center approximation for ERIs, the distributions are usually expanded in a new set of auxiliary functions. As suggested by Vahtras et al. [8], to reproduce the total energy to good accuracy, it is necessary to reproduce the behavior of the exact expansion near the atomic nuclei. On the other hand, the products between atomic orbitals belonging to two distinct nuclei are likely to have small amplitudes and/or to be similar to the other products; for instance, the interatomic distribution between s and s primitives would induce an electrostatic field similar to those induced by the other distributions such as between s and p primitives of the order of F,,(T) [l]. As a result, the effective rank of the distributions in the bonding region would be much smaller than that of the entire distributions. In this work, we therefore approximate the oneelectron density and the entire electron distributions by a linear combination of atomic electron distributions (LCAD),

where x,(r) and x,,(r) in the summation are couples of basis functions which share the coincident function centers. From the definition of ADS, the ERIs are divided into three categories, namely, four-, threeand two-center ERIs, the numbers of which are proportional to N4, N3 and N2 respectively, in polyatomic molecules. In the LCb approximation, we do not calculate the entire set of ERIs, but evaluate the four-center ERIs from the data of three-

S. Ten-no, S. Iwata / Chemical Physics Letters 240 (1995) 578-584

580

and two-center quantities. Using indices K, ADS, the ERIs are approximated [8,9] by (UblCff)

=

z(UblK)K;‘(hlCd),

A,.

..

for

(2)

Kh

or, in the orthogonalized (ablcd)

=

C(~bl~‘][~‘ld),

form [9], (3)

d

where V is the two-center ERIs in matrix notation and the square bracket denotes densities transformed by a vector decomposing V- ‘. One can consider that the accuracy of this expression would be improved with the quality of molecular basis function, since the convergence of distributions is doubly faster than that of orbitals. Therefore, the inner-projection of the above expressions in Eqs. (2) and (3) can be expected to become accurate for a large set of basis functions. We apply the present three-center approximation for ERIs with LCAD, which henceforward we simply call the LCAD approximation, to SCF and MCSCF wavefunction calculations. For the LCAD version of SCF calculations, one yields three-center expressions of Coulomb and exchange operators using Eq. (31, &?;= K$=

c

~&,(ablK][Klcd),

~c~(.il~][~ibi), K

(4) (5)

i

where P is the density matrix and i denotes occupied orbitals. These approximate operators satisfy the inequalities connected with the corresponding true operators, J3”
enter the CI part. Then, from the data of the first and the second order density matrix, the Lagrangian and Hessian are generated. In this step, a part of the Coulombic contributions can be calculated without preparing the four-index integrals. The accuracy of 10-3-10-4 for the approximate Hessian is found to be sufficient to ensure the same convergence. Further computational details for the LCAD version of MCSCF calculations will be presented in a forthcoming publication [12].

3. Result In order to investigate the behavior of the LCAD approximation, a couple of basis sets, namely STO3G minimal set [13] and Huzinaga-Dunning double zeta set [14], are used. We further examine the addition of polarization functions to the double zeta set for the first row atoms and for all of the atoms in the molecules [14]. The accuracy of the calculation depends both on a threshold to remove nearly-linear dependent vectors among the eigenvectors of the positive-definite two-center ERI matrix V, and on a cut-off criterion on the integrals. Generally speaking, the integrals should be computed to an accuracy of at least lo-* to ensure total energies correct to 10V6 hartree. In our calculation, the eigenvectors of V whose eigenvalues are smaller than a threshold of 10d9 are removed. The same value of the cut-off criterion is used for all types of integrals appearing in the present approximation. As an example, we show the distribution of eigenvalues of the ERIs for the neon atom in Fig. 1; they are ADS in the present case. Apparently, there exist 3 vectors for [3s2p] and 25 vectors for [3s2pld] of smaller eigenvalues than 10-9. These redundant linear dependences which do not alter the result of the total energy are classified in two categories. The first called the angular redundancy (AR) is due to the duplicate orbital products having the same total angular momentum, and the second called the Coulombic redundancy (CR) is due to the overcompleteness of ADS with respect to ERIs. In reality, there is no CR among ADS of these basis sets with this criteria in Fig. 1 and all of the linear dependent components are assigned to be AR. However, a large amount of CR could appear in the

S. Ten-no, S. Iwara / Chemical Physics Letters 240 (1995) 578-584

0

[2SlPl

fa

[3~21-4

1

3456189

-Log

Y

Fig. 1. Distributions of eigenvalues v of the V matrix for the neon atom.

581

use of a larger basis set such as triple zeta plus polarization function. In our code, the AR can be easily removed in generating the ADS, and it improves the timings for diagonalizing V and generating transformed integrals (ab 1 K] to remove CR. Table 1 shows the results of formaldehyde for four basis sets. One can see that the errors of both the total and orbital energies in the LCAD approximation decrease as the quality of basis functions is improved. This trend can be roughly elucidated by the ratios of the numbers of ADS to AOs, 2.67 for [2slp/ls], 4.09 for [3s2p/2s], 5.76 for [3s2pld/2s] and 5.50 for [3s2pld/2slp]. The increasing errors in 3C + 2C, where four-center ERIs are completely neglected, are due to the iterative nature of the calculation [9] along with the increasing variational

Table 1 Total and one-electron energies a of formaldehyde b

E (AE) E1 82 63

&4 ES &6

&7 &8

[2slp/lsl Full d

12/32 ’ LCAD

3c + 2c c

[3Qp/2sl Full

22/90 LCAD

3C+2C

- 112.352322

- 112.348336 (0.003986)

- 113.527268 (- 1.174946)

- 113.829150

- 113.829348 ( - 0.000198)

- 116.018733 (- 2.189583)

- 20.3177 - 11.1303 - 1.3573 - 0.8074 - 0.6464 - 0.5522 - 0.4590 - 0.3536

cLf

E1 E2 83 E4

Es &6

El &8

CL

1.489

- 20.3150 - 11.1306 - 1.3559 - 0.8073 - 0.6451 - 0.5507 - 0.4576 -0.3521 1.515

- 20.3437 - 11.1427 - 1.7621 - 1.1567 - 1.0111 - 0.6688 - 0.6086 -0.2758 1.645

- 20.5876 - 11.3550 - 1.4449 - 0.8742 - 0.7144 - 0.6510 - 0.5457 - 0.4423 3.006

- 20.5868 - 11.3548 - 1.4443 - 0.8745 - 0.7137 - 0.6504 - 0.5451 - 0.4427 3.032

- 20.8747 - 11.3964 - 2.1707 - 1.4519 - 1.2750 - 1.1288 - 0.9634 - 0.5162 1.864

[3s2pld/2s] Full

34/196 LCAD

3C+2C

[3s2pld/2slp] Full

40/220 LCAD

3c + 2c

- 113.891983

- 113.891933 (0.000051)

- 117.076117 ( - 3.184134)

- 113.895313

- 113.895300 (0.000013)

- 117.273377 ( - 3.378064)

- 20.5771 11.3451 - 1.4175 - 0.8703 - 0.7006 - 0.6558 - 0.5427 - 0.4400

-

2.732

- 20.5770 - 11.3451 - 1.4175 - 0.8703 - 0.7005 - 0.6557 - 0.5427 - 0.4400 2.734

- 21.0510 - 11.3410 - 2.3622 - 1.5823 - 1.3699 - 1.3138 - 1.1235 - 0.8029 0.201

- 20.5769 - 11.3462 - 1.4173 - 0.8699 - 0.6998 - 0.6552 - 0.5423 - 0.4394 2.738

a In E,. b A geometry of type I in Ref. [15] is used. ’ Total number of AOs/effective d Conventional. ’ Without four-center ERIs. f Dipole moment in debye units.

- 20.5769 - 11.3462 - 1.4173 - 0.8699 - 0.6998 - 0.6552 - 0.5422 - 0.4394 2.738

- 21.0655 - 11.2452 - 2.3764 - 1.6153 - 1.4078 - 1.3287 - 1.1321 - 0.8009 0.351

number of ADS without redundancies.

S. Ten-no, S. Iwata / Chemical Physics Letters 240 (1995) 578-584

582

lE-02

1E-03 f “,

lE-04

i IE-OS

lE-06

-

blplls]

. . ..0 ........

[3s2p,~~

----O----

[3s2pld2s]

----A---

[3s2pldZrlp]

Fig. 2. Absolute errors of SCF energies using the LCAD approximation.

degrees of freedom. Using the [3s2pld/2slp] set, the LCAD approximation recovers an amount of 99.99997% error caused from the absence of fourcenter ERIs. In the augmentation of d functions to the first-row atoms, the number of crude ADS increases by 150, but one third of them are discarded as AR in correspondence with the case of neon. Fig. Table 2 SCF energies a for the dissociation

of cis-glyoxal

’ using the LCAD approximation

Basis set

Integrals

cis-glyoxal

cis-glyoxal

i2slp/lsl

full LCAD

- 223.574509 - 223.565378 (0.009131)

(AE) [3s2p/2sl

full LCAD (BE)

[3s2pld/2s]

full LCAD

(AE) [3s2pld/2slp]

full LCAD

(AE)

2 shows the absolute errors in the LCAD approximation for several molecules; the optimized geometries at 3-21G [~~I/scF are used for all molecules. Similar tendencies with respect to the basis functions can be seen as in the case of formaldehyde. It should be noted that, in all of the examples calculated here, the LCAD errors with [3s2pld/2slp] set are less than 0.1 kcal/mol, which can be considered sufficiently accurate for ordinary chemical interests. This indicates that ADS containing polarization functions guarantee practical accuracy without any auxiliary functions in the bonding region. It would be required to add several auxiliary functions to reproduce the total energies even more qualitatively in the minimal and in the double zeta basis functions. We then show the results for the dissociation reaction of cis-glyoxal into 2C0 and H, in Table 2. The activation barriers are reproduced almost in quantitative agreement with those in the conventional scheme. In the SCF energies with the [3s2pld/2s] set, the absolute error increases as the chemical reaction proceeds. This originates from the internuclear distance between two accessing hydrogen atoms in the reactions. The discrepancy is also improved by the addition of p functions to the hydrogen atoms. Theoretical values for the dissociation energy are also reproduced except with the minimal basis set.

e

2C0 + Hz

E, d

Ed

- 223.385613 - 223.371476 (0.014138)

- 223.563057 - 223.560007 (0.003050)

118.449 121.588

7.186 3.371

- 226.506079 - 226.505668 (0.000411)

- 226.359253 - 226.355819 (0.003433)

- 226.494733 - 226.492978 (0.001755)

92.069 93.964

7.119 7.963

- 226.632369 - 226.632295 (0.000074)

- 226.499535 - 226.499132 (0.000404)

- 226.645425 - 226.646348 ( - 0.000923)

83.295 83.501

8.193 8.819

- 226.635885 - 226.635826 (0.000058)

- 226.504904 - 226.504858 (0.000046)

- 226.649971 - 226.649949 (0.000022)

82.133 82.125

8.833 8.856

a In E,. b Theoretical geometries at double zeta plus polarization/SCF ’ In a geometry of planar stationary point. d Activation barriers in kcal/mol. ’ Dissociation energies in kcal/moI.

(TS) ’

in Ref. [17] are used.

S. Ten-no, S. Iwata / Chemical Physics Letters 240 (1995) 578-584 Table 3 Characteristic

CPU times a, energies and the total number of ERIs used for the MCSCF calculation LCAD-MCSCF

583

of allene

Conventional

A0 integrals recorder integrals diagonalize V matrix orthogonalize integrals

13.3 2.3 12.4 30.1

MO integral transformation GUGA sort CI and density matrix Lagrangian and Hessian Newton-Raphson

14.1 1.2 0.2 22.3 d 20.2

105.2 2.5 0.2 12.3 20.2

one MCSCF iteration

51.9

140.4

total energy ’ (energy difference) A0 integrals MO integrals

- 115.784396 (0.000051)

- 115.784447

2.0 Mwd 2.3 Mwd

0.6 Mwd 0.7 Mwd

a Time in CPU seconds on a HP9000/735 computer. b In E,. ’ GAMESS quantum chemistry code [18] is used. d Three-index integrals whose eigenvalues are larger than 10m3 are assembled

Finally, Table 3 shows the results of an MCSCF calculation of allene in C, symmetry, using the [3s2pld/2slp] set totalling 65 basis functions and 345 ADS. There are 9 core orbitals, and the active space consists of 4 electrons in 4 orbitals, which generates 20 configuration state functions. In a conventional MCSCF calculation, 75% of the total CPU time is spent in the MO integrals transformation step. The LCAD version of the MCSCF calculation requires only 13% of the time required for the same step in the conventional calculations. This leads to an approximately 60% speedup for one iteration and accordingly for the total time, maintaining an accuracy within 0.1 kcal/mol of the true MCSCF energy. With larger basis sets, the efficiency of the present approximation turns out to rapidly increase both in CPU time and accuracy [12]. Further details of such calculations will be presented in a forthcoming article. The number of active ADS required in the present LCAD scheme is typically 5-6 times larger than the corresponding AOs with d functions. This is more efficient than the factor of 10 observed in the paper by Vahtras et al. [8], in which they used uncontracted

’

30.0 13.9 -

in the Hessian.

auxiliary functions. However, we have no meaningful timing data for the SCF calculation, since there still remain two steps whose computational works scale as N4. One is the integral transformation step to generate (ab ( K] and the other is the Fock matrix assembly for the exchange part in Eq. (5). Most of the computation time was spent for these steps. Although additional investigations are required in these areas, an efficiency of the present approximation for MCSCF wavefunctions is revealed. The LCAD approximation would also be efficient for other correlated wavefunctions such as MP2 and coupled cluster levels, as shown by previous investigators using three-center expansion [7,8].

Acknowledgement The authors are indebted to Dr. R.J. Harrison for his useful suggestions. A part of the work was supported by Grants-in-Aid for Scientific Research and for Priority Area (No. (No. 04640458) 04243102).

584

S. Ten-no, S. Iwata / Chemical Physics Letters 240 (1995) 578-584

References [l] S. Obara and A. Saika, J. Chem. Phys. 84 (1986) 3963. [2] R. Lindh, U. Ryu and B. Liu, J. Chem. Phys. 95 (1991) 5889. [3] P.M.W. Gill, Advan. Quantum Chem. 25 (1994) 141. [4] S. Ten-no, Chem. Phys. Letters 211 (1993) 259. [S] B.H. Greeley, T.V. Russo, D.T. Mainz, R.A. Friesner, JeanMarc Langlois, W.A. Goddard III, R.E. Donnelly Jr. and M.N. Rim&da, J. Chem. Phys. 101 (1994) 4028, and references therein. [6] C. Van AIsenoy, J. Comput. Chem. 9 (1988) 620, and references therein. [7] M. Feyereisen, G. Fitzgerald and A. Komornicki, Chem. Phys. Letters 208 (1993) 359. [8] 0. Vahtras, J. AImlijf and M.W. Feyereisen, Chem. Phys. Letters 213 (1993) 514; J. AImlof, in: Lecture notes in quantum chemistry II. VII, ed. B.O. Roos (Springer, Berlin, 1994). [9] A.P. Rendell and T.J. Lee, J. Chem. Phys. 101 (1994) 400. [lo] B.I. Dunlap, J.W.D. Connolly and J.R. Sabin, J. Chem. Phys. 71 (1979) 3396.

[ll] [12] [13] [14]

[15] 161 171

181

N.H.F. Beebe and J. Linderberg, Intern. J. Quantum. Chem. 12 (1977) 683. S. Ten-no and S. Iwata, in preparation. W.J. Hehre, R.F. Stewart and J.A. Pople, J. Chem. Phys. 51 (1969) 2657. T.H. Dunning Jr. and P.J. Hay, in: Methods of electronic structure theory, ed. H.F. Schaefer III (Plenum Press, New York). S. Ten-no, F. Hirata and S. Kato, J. Chem. Phys. 100 (1994) 7443 J.S. BinkIey, J.A. Pople and W.J. Hehre, J. Am. Chem. Sot. 102 (1980) 939 Y. Yamaguchi, Y. Osamura, J.D. Goddard and H.F. Schaefer III, A new dimension to quantum chemistry: analytic derivative methods in ab initio molecular electronic structure theory (Oxford Univ. Press, Oxford, 1994). M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S.J. Su, T.L. Windus, M. Dupuis and J.A. Montgomery, J. Comput. Chem 14 (1988) 158.