A comparison of the efficiency and accuracy of the quadratic configuration interaction (QCISD), coupled cluster (CCSD), and Brueckner coupled cluster (BCCD) methods

Volume 190, number I,2

CHEMICAL PHYSICS LETTERS

28 February 1992

A comparison of the efficiency and accuracy of the quadratic configuration interaction (QCISD ) , coupled cluster (CCSD ) , and Brueckner coupled cluster (BCCD ) methods Claudia Hampel, Kirk A. Peterson and Hans-Joachim Werner Fakultiit ftir Chemie, Universitiit Bielefeld, W-4800 Bielefeld, Germany Received 30 December 199 1; in final form 9 January 1992

The coupled-cluster method restricted to single and double excitations from a closed-shell reference function (CCSD) and the corresponding quadratic configuration interaction method (QCISD) are formulated in terms of quantities which can be computed directly from the two-electron integrals in A0 basis. A simple yet effective method to accelerate convergence in Brueckner coupled-cluster (BCCD) calculations is also described. Using this procedure BCCD calculations require no more effort than CCSD calculations. In order to compare the accuracy of all three methods potential energy functions and spectroscopic constants have been calculated for NS, CO, F2, and HF using large basis sets.

1. Introduction The coupled-cluster (CC) method originally proposed by Ciiek and Paldus [ l-3 ] is now one of the standard tools of quantum chemistry for computing correlated electronic wavefunctions. Purvis and Bartlett [4] were the first who implemented it including all single and double excitation operators (CCSD). Later work also included triple excitation operators [ 5,6]. The quadratic configuration interaction method (QCISD) recently proposed by Pople and co-workers [ 7 ] can be viewed as an approximation to the CCSD method. Another variant is the Brueckner CCD method [ 8-141, in which the f, amplitudes are eliminated by a suitable unitary transformation of the molecular orbitals. The main advantage of the CCSD and QCISD approaches as compared to the corresponding contiguration interaction (CISD) method is size-consistency. Since this becomes increasingly important with increasing number of electrons, CC methods should be particularly well suited for the description of electron correlation effects in molecules with a large number of electrons or with multiple bonds. The bottleneck in calculations for large systems is often

the storage of the two-electron integrals. Therefore, “integral-direct” methods have been developed [ 15 1, in which the two-electron integrals are recomputed whenever needed and never stored on disk. So far such direct computations have been limited to the Hartree-Fock self-consistent field (SCF) approximation [ 15 ] and Meller-Plesset perturbation theory [ 16- 18 1. Recently, we have also developed integralmulticontiguration direct self-consistent field (MCSCF) and multiconfiguration-reference configuration interaction (MRCI) methods [ 191. The common feature of these methods is that only partial integral transformation is needed. This integral transformation is done at the same time the twoelectron integrals are computed and the comparably small number of transformed integrals are stored on disk as in conventional methods. In the case of MCSCF calculations, the transformation has to be done once per iteration. In the case of MRCI calculations the transformation is needed only once, but a set of “external exchange operators” has to be computed from the A0 integrals in each iteration. These operators account for all contributions of two-electron integrals with three and four external orbitals. In the present work we show that a similar formu-

0009-2614/92/$ 05.00 0 1992 Elsevier Science Publishers B.V. All rights reserved.

1

Volume 190, number

1,2

CHEMICAL

PHYSICS

lation of the CCSD method is possible. As compared to the CI case, however, an additional generalized partial integral transformation is needed once per iteration. This transformation accounts for contributions of the disconnected triples in the T2 equation. Since these terms are not present in the QCISD approximation, integral-direct QCISD calculations will be easier and cheaper than CCSD calculations. As has been discussed in detail in two recent papers by Handy, Pople, and co-workers [ 12,13 1, the Brueckner CCD approximation offers several conceptual advantages over the CCSD or QCISD methods. For instance, the simplicity of the wavefunction facilitates the calculation of energy gradients. Moreover, it appears that Brueckner CCD wavefunctions are more stable with respect to spin contamination and other symmetry breaking problems than unrestricted Hat-tree-Fock wavefunctions. In the past, however, Brueckner calculations were considered too expensive to be competitive with other methods because an integral transformation had to be performed each time the orbitals are modified [ 10,111. In order to minimize the number of transformations, complete CCD (or CCDT) calculations were performed between each orbital transformation [ lo13 1. The computational effort was therefore at least 2-3 times larger than in a calculation with fixed orbitals. Another problem mentioned by Handy [ 20 ] is that slow convergence was often encountered. To circumvent this problem, a quite sophisticated NewtonRaphson optimization technique was proposed in ref. [ 131. This method, however, does not fully account for the coupling of the orbital changes with the T, amplitudes, and is therefore probably not very effective. In the present article we propose a very simple alternative method which is based on successive exponential transformations of the orbitals. This transformation is described by a unique and non-redundant set of parameters, which can be used in the usual DIIS (direct inversion of the iterative subspace) extrapolation technique [ 2 11. It will be demonstrated that with this technique the speed of convergence in BCCD calculations is similar to that in CCSD or QCISD calculations. Within the present A0 integral based formulation, which eliminates the expensive full transformation of the two-electron integrals, the computational expense of BCCD calcu2

LETTERS

28 February

1992

lations will be shown to be intermediate between the somewhat cheaper QCISD and the slightly more expensive CCSD calculations. In section 2 we will summarize the equations as implemented in the present work. In section 3 the optimization of Brueckner orbitals will be described. In section 4 the efficiency of the method will be demonstrated for some examples with a relatively large number of correlated electrons and basis functions. Finally, in section 5 benchmark calculations for Nz, CO, F2, and HF with very large basis sets will be presented. The results of these calculations will be compared to previous accurate MRCI calculations and experimental data.

2. The CCSD equations Our formulation of the closed-shell CCSD equations closely follows the self-consistent electron pair theory (SCEP) [22] given by Pulay, Saebo and Meyer [ 23 ] for the closed-shell CISD and CCD cases, but differs from this in the inclusion of the singles and the definitions of the intermediate quantities. Also the equations of Scuseria et al. [ 241 are closely related, even though they were not explicitly formulated in terms of matrix operations. The coupled-cluster wavefunction is defined as I U =w(O

IO> ,

(1)

where f is an excitation operator and 10) is the closed-shell reference wavefunction. The operator F can be written as f=?,

+f1+

...

(2)

with F, = c t$ai) ai

(3) (4)

Restricting F to F, + F yields the CCSD approximation. The & are the usual spin-coupled one-particle excitation operators, i.e. L+al=qt,?/i+rIfitT,*

(5)

Here and in the following i, j, k, 1 ... refer to the internal orbitals (occupied in the reference function),

CHEMICAL PHYSICS LETTERS

Volume 190, number 1,2

a, b, c, d, ... to the external orbitals, and r, s, t, u ...

to any orbitals. Equations which are sufficient to determine the amplitudes tf and Tibb=TJ& can be obtained by projecting the Schrijdinger equation from the left with the configurations I@P)=&lO)

>

I@$>=~~i~bjIO)

(i>j,alla,b).

[20$b+@]

.

I a$)

=2(akllb)-(allkb).

(16)

IfkIn= [Flak

(9)

and the normalization

=h&+ [kk”la=

1 [2(akIii)-(ail&)],

(17)

I

[Kk’],,=(aklZi)

,

(18)

[lk’~],=[Lk’],i=2(aklZi)-(allki),

(19)

where F is the closed-shell Fock matrix in the MO basis. We further define external exchange matrices K(D’i),=

(6$bI@$b)=1+&j&,.

ab

(7)

(8)

=S,,S,SikSj,+Sol,S,S,Sj~

(Lk’),b= [2Kk’-Kfk]

The integrals with one external orbital will be represented by the column vectors

This yields the orthogonality relation ( *zb

external orbitals. For convenience we also define the matrices

(6)

Since the doubly excited configurations are not orthogonal, it is convenient to use the orthogonal complement [23] of the space {0$}, 5+;

28 February 1992

1 D~u(rtlus). 111

(20)

(10)

In terms of these functions, the CCSD equations are defined as E=(oIri(l+~,++*+~~‘:)~o),

(11)

The quantity K( D”) (k) will denote the kth column of the matrix K( DO). Finally, we introduce the matrices (E’i),=G,t’,

(21)

v:, = (09 I (A-E) x

I+i,+f~+;f:+F,T*+;r; (

lO)=O, .

> (12)

vib=(d$b@-E)

(i>j, all a, b) .

l+f+,+f+z+;f+;+f,p,

(13)

The QCISD equations [ 7 ] are obtained by omitting all ?:, f’:, p’:termsandthe F,If;termineq. (13). For efficient use of modem computer hardware it is essential to compute the residual vectors vi and matrices Vu in terms of matrix operations. For this purpose we define internal Coulomb and exchange matrices (Jk’)ab = (k/lab) ,

(Kk’h=(akllb),

(14)

(15)

which represent the two-electron integrals with two

which in the MO basis have only one non-vanishing row. In analogy to eq. (20) these are used to define the generalized Coulomb and exchange matrices J(EU),,=

1 (ablic)tl,,

(22)

K(E”),,=

c (ailbc)t’,, C

(23)

C

which represent contributions of integrals with three external orbitals. The explicit expressions for the residual vectors Y’and matrices V” now take the simple form Yi=gi+ : [ (2Tik-Tki)rk-/&id‘]

,

V”= K’j+ K(W) + ; aij,,,Ck’+Gij+@t

(24) ,

(25)

where C”=T’i+titit

9

Dij=CU+ Ed+ pt 3

(26) (27)

and 3

Volume

G’JcSilit

190, number

_

+TQX+

CHEMICAL

1,2

$ [ K( D”)‘k’+k”k]i!=‘f-

PHYSICS

T /l&j

‘g [ (2T;k_Tk’)yki_fTkiZkj_

(Tkizki)t] (28)

(note that Tk’=Tzkt). The intermediate quantities used in eqs. (25) and (28) are scalars (II.. %kl=K$!+tr(G”K’k) + (fitpli+gtp)

(29)

Pki=Fkr+ 7 [tr(LkC~i)+r~tPki]+f~Ii:

(30)

vectors rk=f k+ 1 Lk’t’,

(31)

I

kl

+

$

[2K(Dik)-K(Dki)lck),

LETTERS

28 February

1992

left out. In particular, the operators J ( EU) and K( Eti) are not needed. As now used in various implementations of configuration interaction methods [ 22,23,25-271, the operators K’(D”) can be obtained directly from the A0 integrals by first transforming the matrices Do into the A0 basis, then computing the operators in this basis, and finally transforming them back into the MO basis. In the CISD and QCISD cases this avoids the necessity of a full integral transformation. As shown above, for the CCSD method one needs the additional operators J ( E”) and K(E”) . These can be obtained directly from the A0 integrals by a generalized partial integral transformation. In the following, quantities in the A0 basis will be denoted by a tilde, and the indices p, v, p, 0 ... will refer to AOs. Let R be the MO coefficient matrix with column vectors P and

(32) I..

E”

and matrices

PO

A= 1 Lk’Tlk,

(33)

-. = r” tJ PO.

(39)

The operators are first computed in the A0 basis

kl

X=F-A+

$ [2J(Ekk)-K(Ekk)]-

(40)

‘$ rktkt, (34)

(41) +f C Lkj(2CU_e/)_j I Zb=JkJ+J(EkJ)-

7

C I

p$~t+f Qjt,

[K’ke’+k’kjtlt]

,

(35)

.

(37)

As mentioned above, this formulation of the CCSD equations is equivalent to that given by Scuseria et al. [ 24 1. In fact, the T1 equations are easily shown to be identical, but ours and Scuseria’s T2 equations differ by terms which vanish when the T, equation is fulfilled (i.e. all v’=O). The QCISD equations are obtained from the above quantities by setting C”=T”, G”= tT”, and by replacing si by f i in the first term of eq. (28). Furthermore, in the second term of eq. (28) K( DO)(k) hastobeomittedandineqs. (29)-(31) and (33)(36) all contributions of the f, amplitudes must be 4

J(Eii)=Bt&i?J)fl,

(42)

K(E”)=WtR(E”)R.

(43)

(36)

with Gij= fTii+tiiit

and finally transformed into the MO basis

This procedure, which is similar to the computation of the operators Jk’ and K”, requires about ] N 4m + 4N ‘m * additions and multiplications (N basis functions, m correlated orbitals) rather than ;N’m* operations if the same quantities are computed from the fully transformed two-electron integrals (the full integral transformation scales with N5). The additional effort is, however, quite insignificant as compared to the tN4m2 operations (or 4N4m2 if the integrals are ordered in triples (PVlpo), (ppplau), (p~lpv) ) needed to evaluate the operators K (Dg ) and will therefore not introduce a bottleneck. We note that as in Scuseria’s treatment [24] the evaluation of the Gu, Yu, and ZkJterms requires only 4m3N3 operations, while 7m3N 3 operations were

needed

in the vectorized CCSD method of Lee and Rice [ 281. In order to achieve this reduction, it is not necessary to store an intermediate product of the size N4, as assumed in ref. [28]. Also, there is no out-of-core sort needed in each iteration. In our view, a significant advantage of avoiding the full integral transformation in CCSD calculations is the reduced disc space needed. Additionally, avoiding the full transformation is a prerequisite for the efficient implementation of the BCCD method and for integral-direct methods in which the storage of the two-electron integrals in A0 basis is entirely avoided. In the latter case QCISD will have a delinite advantage over CCSD since the operators J ( E”) and K( E”) are not needed. Such integral-direct methods are presently under development in our laboratory.

3. DIIS extrapolation in Brueckner calculations In the iterative solution of the CCSD or QCISD equations, the update of the amplitudes is obtained by first order perturbation theory, i.e. At6=-Vi/(ea_ei), AT~~=-V~~l(Ea+Eb-Ei-Ej)

(44) )

(45)

where t, are the eigenvalues of the Fock matrix. The simplest way of obtaining Brueckner orbitals is to absorb the Ar’ into the occupied orbitals in the beginning of each iteration h=@i+

C

a

tS@c2

28 February 1992

CHEMICAL PHYSICS LETTERS

Volume 190, number 1,2

(46)

with subsequent reorthornomalization. Afterwards, the f’ are set to zero. Unfortunately, the convergence of this procedure is often slow, which is due to a relatively strong coupling of the T ij amplitudes with the orbital changes. In the following we outline a simple extrapolation technique which considerably improves convergence. As demonstrated by Scuseria et al. [ 291, convergence in CCSD or QCISD calculations can be much accelerated by the DIIS extrapolation technique of Pulay [ 2 11. In the following we assume that all amplitudes TV and t’ form a vector T, and the residuals V’/ and vi form a corresponding vector Y. The DIIS method is based on the assumption that the residual

vector Y is approximately linear in small changes of the amplitudes T, i.e. that for withCI,=l n

T= C2,T(“)

n

(47)

the residual can be approximated by Y= Cl,?@‘. n

(48)

Here Y(“) are the residual obtained with the trial vectors T W) of the nth iteration. The coefficients ;1,, are determined by minimizing Y+V= 1 tr(V’j+V”)+ C (v’+y’) , I

iaj

(49)

or, alternatively, ATtAT=

C tr(AT”AT”)+ i>,J

C (At’+At’) I

(50)

subject to the auxiliary condition C,,n,,= 1. This yields a small set of linear equations in the unknown expansion coefficients A,. In the first case, the optimized residual vector Yis used to determine an update vector AT according to eqs. (44) and (45 ). A new expansion vector T("+ ‘) is then obtained by adding AT to the extrapolated T from eq. (47). In the second case, which corresponds only to a different weighting, the next expansion vector is found by adding the optimized AT to the extrapolated T. We found that both variants work about equally well. This method cannot straightforwardly be used for the optimization of Brueckner orbitals, since the f’ amplitudes go to zero and the coupling of the TV to the orbital changes would not be included. For the same reason DIIS extrapolation of the matrices To alone is not effective. The following procedure circumvents this problem. Instead of the simple absorption step of eq. (46) we employ unitary orbital transformations U=exp(A)=l

+A+tA’+

...

with flai=tL,

Ain=-&

Aij=flob=O.

(51)

In the nth iteration the accumulation transformation of the original orbitals is U(“)=exp(A(“))

=exp(h(“))

... exp(A(*)) exp(A(‘)) .

(52) 5

Volume 190, number

1,2

CHEMICAL

The transformation matrix Utn) is determined antisymmetric matrix A(“) ‘“‘_I)-f(U’“‘-I)2

A(“)=ln(IJ(“))=(IJ +f(U’“‘-1

by the

y+

...

from which a set of independent (Fa)(n)=(A(“)),i

(53) parameters (54)

can be extracted. The exponential and logarithmic expansions are necessary since the matrices A(“) do not commute and therefore A(“)# ,X;=, Acm). Using the vectors (?) (n) instead of (ii) (n) in the DIIS extrapolation yields improved vectors ?, which are employed to construct a new transformation matrix and improved orbitals. Convergence of this procedure is found to be equally fast as in the QCISD and CCSD methods with DIIS extrapolation as outlined above. The essential point is that always the original orbitals are transformed. The parameters ? describing this transformation converge to definite values rather than going to zero, and their coupling to the Tu amplitudes is appropriately taken into account. We note that in order to avoid numerical difficulties near convergence it is important to include a sufficiently large number of terms in the expansions of the exponential and logarithmic functions. In practice, we found that terms up to the sixth power are always sufficient. In the DIIS procedure the trial vectors 2”“’ and residual vectors Y(“) of each iteration must be stored on disk. Typically, we start DIIS extrapolation in the third iteration (not counting the initial one, which corresponds to MP2 ), and use at most six expansion vectors. Additional vectors then overwrite those with the smallest coefficients (Y, (which are usually the oldest ones).

4. Test calculations

In this section, we present some test calculations in order to analyse and compare the computational effort of the CISD, QCISD, CCSD, and BCCD methods. These methods have been implemented as 6

PHYSICS

LETTERS

28 February

1992

part of the MOLPRO suite of ab initio programs #I. This package uses symmetry adapted basis functions for Abelian point groups ( DZh and subgroups). Only the non-vanishing symmetry blocks of the matrices are stored and used in the matrix operations. Table 1 shows a timing analysis of some calculations for benzene and naphthalene. In both cases, all valence electrons were correlated and D2,, symmetry was used. The first CCSD calculation for benzene with a valence double zeta plus polarization (VDZ) basis set [ 301 takes only 11.5 s per iteration on a CRAY-YMP. About 28% of this time is needed for the operators K (Do) (eq. (20) ), 27% for the matrices G” (eq. (28)), and also 27% for the intermediate quantities Yb and Zb (eqs. (3 5 ) and ( 36 ) ). The partial integral transformation needed to compute the operators J( E”) and K( E”) (eqs. (38)(43) ) requires 15% of the time. In total, the CCSD calculation was about 10% more expensive than the BCCD calculation and 20% more expensive than the QCISD calculation. In the second calculation for benzene a valence triple zeta (VTZ) basis set [ 301 with 2d and If functions on carbon, and 2p functions on hydrogen has been employed (234 contractions). This yields about 16 times more integrals than the VDZ basis. Nevertheless, the computational effort for the K(D”) increases only by a factor of 10, which is due to the better vectorization obtained with the larger block sizes. As expected, the N4 steps become more dominant with increasing basis size (42% for K( De) ), but the effort for the operators J (E”) and K( E”) is still quite small ( 17%). These calculations demonstrate that with our method BCCD calculations are even less expensive than full CCSD calculations. This is due to the smaller effort for the partial integral transformation (there is only about half the number of operators). Moreover, the time per iteration for a CCSD calculation is only 1.6 times larger than for a simple CISD calculation. This ratio becomes somewhat less favourable if the increased number of iterations needed for CCSD is taken into account. In order to converge the energy to lo-’ hartree, 10 it-

MOLPRO is a package of ab initio programs written by H.-J. Werner and P.J. Knowles, with contributions of J. Almlaf, R. Amos, S. Elbert, W. Meyer, E.A. Reinsch, R. Pitzer and A. Stone.

28 February 1992

CHEMICAL PHYSICS LETTERS

Volume 190, number 1,2

Table 1 Timing analysis a) of CISD, QCISD, BCCD, CCSD calculations Quantity

CISD

BCCD

QCISD

CCSD

C6Hs, basis ‘) VDZ ( I 14 bf.), 30 el. correlated, D2,, symmetry ‘) K(D") G" y" 2" J’/

3.04 2.88

3.05 2.89 2.92

3.05 2.91 2.92 1.17

6.36

9.49

10.55

K”

J(E"),K(E")

time/iteration E corrd’

-0.668140

-0.821359

-0.818639

3.03 3.00 3.07 1.73 11.50 -0.820437

C6H6, basis b, VTZ (234 bf.), 30 el. correlated, Dzh symmetry ‘) K(D") Gtj yW zb J" KU

31.2 13.0

31.2 13.2 12.9

31.6 13.1 12.8 8.0

45.3

59.2

66.1

J(E"),K(E")

time/iteration E corrd’

-0.783210

-0.966997

-0.9643793

31.6 13.4 13.3 12.9 73.0 -0.965635

CIOHs,basis b, VDZ ( 180 bf.), 48 el. correlated, DZhsymmetry d, K(D") GU ykj ZkJ J’/

30.5 25.3

30.5 25.3 25.6

30.4 25.4 25.5 6.0

59.3

86.1

92.2

K’,

J(E"),K(E")

time/iteration E cond)

-0.997841

-1.347086

-1.343745

30.2 25.6 26.3 11.3 98.1 -1.345390

a1 All times for CRAY-YMP8/32 in seconds per iteration. b, Correlation consistent basis sets of Dunning [ 301. ‘) Rccz2.62 au, Ro=2.04 au. d, R,= 1.40 A, R,--= 1.08 A, all bond angles 120”.

erations were needed for QCISD, BCCD, and CCSD, but only 8 iterations in the CISD case (this can be reduced to 7 by using the Davidson method [ 3 1 ] instead of the DIIS method for convergence acceleration). The difference in the computational effort for QCI, BCCD, and CCSD becomes even less pronounced in the case of naphthalene. This is because the time required to compute the Gti, Yk’, and Zk’ increases with the third power of the number of electrons, while the dominant step of partial integral transformation scales only linearly. Unfortunately, due to limited disk space available, we were not able to repeat the calculation for naphthalene with the VTZ basis. Even without f functions and only lp on the hydrogen atoms, the integral file was larger than the maximum

of one Gbyte that was available for use. Otherwise it would certainly be possible to perform CCSD or QCISD calculations with more than 300 basis functions. This demonstrates the necessity of the development of integral-direct methods, in which this bottleneck is eliminated at the cost of more CPU time. Table 2 summarizes the CPU times needed for the complete calculations shown in table 1 and some other calculations with less symmetry. In the lowsymmetry cases the disk space limitations were even more stringent and only relatively small basis sets could be used. The integrals were computed using the vectorized program of Almlof and Taylor #*. The #’ J. Almliif and P.R. Taylor: MOLECULE, a vectorized Gaussian integral program.

7

Volume 190, number Table 2 Comparison

I,2

CHEMICAL

of computation

times b, for benzene, naphthalene,

PHYSICS

LETTERS

28 February

I992

aniline, and toluene calculations

Molecule

basis b,

C6H6

C6H6

C6H6

CIOH~

CeNH,

C7Hs

C,Hs

[3s2pld/

[ 4s3p2d/ 3~2~1

[4s3p2dlf/

[3s2pld/

[3s2pld/

[ 3s2pld/ 2SlPl

[3s2pld/

138 CsC’

138 C,

794 78 1914 2045 2288

985 183 5676 6062 6510

2SlPl

3~2~1

2SlPl

2SlPl

No. of functions

114

192

234

180

symmetry

DZh

D2h

DZh

D2h

133 C 2”

integrals d, SCF QCISD ‘) BCCD ‘) CCSD ‘)

125 1.5 98 109 118

298 6.5 360 399 433

1247 11.7 624 716 762

638 8.2 967 1036 1101

378 9.9 478 517 561

a) b, ‘) d, e,

All times in seconds for CRAY-YMP8/32. Generally contracted correlation consistent valence double zeta and triple zeta basis sets from ref. [ 301. The benzene ring lies in the symmetry plane. Including integral sort. Ten iterations were needed in all cases to converge energy to 10e6 hartree. All valence electrons were correlated.

integral times include input processing, computation of the one-electron integrals, and sorting of the twoelectron integrals. Table 2 demonstrates that for cases with high symmetry and no f functions the integration times are of the same order of magnitude as the times for the CCSD calculations. Hence, in these cases integral-direct methods would be quite expensive. This is even more so the case if f functions are included in the basis, since this strongly increases the integration time. On the other hand, in low symmetry cases with small basis sets the integration times are less dominant, and integral-direct methods would be quite efficient. The number of CSFs in the toluene calculations with C, and C, symmetry (depending on the orientation of the CH3 group) was 1082764 and 207 1630, respectively. In the latter case 50% of the CPU time was spent to compute the operators K( D”). This part could be reduced by a factor of two by reordering the two-electron integrals into triples (ptvjpa), (wplav), (polpv). In the current version of our program this order is only present for integrals over basis functions of different symmetry, and hence calculations with no symmetry are particularly slow. We note that the calculations in the last column of table 2 ran with an average speed of 245 Mflop/s on one processor of the CRAY-YMP 81 32. Finally we mention some timing comparisons with 8

2SlPl

other programs. Lee and Rice [ 28 ] have given timings for CCSD calculations on N2Fz. Their vectorized CCSD program needed 54 s CPU time per iteration on a CRAY-XMP. Our program needed 10.2 s per iteration on a CRAY-YMP, which would correspond to about 13 s on a CRAY-XMP. A more direct comparison was possible with some CCSD calculations provided by Scuseria [ 321. The first one was for HSOH with 50 basis functions in C, symmetry and all electrons correlated (for details see ref. [ 3 ] ). Using Scuseria’s program this calculation took 68 s per iteration, while our program needed 72 s per iteration (both on IBM 6000/550 workstations). The second calculation was for s-tetrazine with 106 basis functions in DZh symmetry (for details see ref. [ 33 ] ; no virtual orbitals were frozen). This calculation took 42 and 3 1 s per iteration with Scuseria’s and our programs, respectively. These results demonstrate that the advantage of avoiding the full integral transformation does not lead to a significant increase of CPU time, even though the operation count per iteration is somewhat higher.

5. Benchmark

calculations

for diatomic molecules

In order to investigate the accuracy of the CCSD, QCISD, and BCCD approximations, several bench-

Volume 190,number 1,2

CHEMICAL PHYSICS LETTERS

mark calculations for small molecules have been published in which the computed energies were compared to full CI (FCI) results [4,7-l 31. Since such

benchmark calculations are limited to very small basis sets, it is of interest to investigate how well their results can be extrapolated to large basis sets. The only possibility to judge the quality of different correlation treatments in large basis set calculations is by comparison to experimental values. Such comparisons are only meaningful, however, if the basis sets employed are sufftciently complete to render basis set errors small as compared to errors in the electron correlation treatment. In this section we report some calculations for the diatomic molecules NZ, FZ, CO, and HF with such large basis sets. These molecules were chosen because (i) they represent different types of molecular bonds; (ii) the computed spectroscopic properties can easily be compared to experimental data; and (iii) the results can be compared to accurate multi-reference CI calculations [ 34,351 carried out with the same basis sets. For each molecule we have used the largest van Duijneveldt [ 361 sp basis, 13s8p, with the innermost 8s, 4p contracted together according to the atomic SCF 1s orbital. Polarization and angular correlation effects are represented using Dunning’s [ 301 optimized 3d2flg primitive Gaussian set. This basis, which is the same as in ref. [ 34 1, will be denoted basis A. For CO, FZ, and HF, we have also used larger sp basis sets [ 371 (basis B) with additional diffuse d, f, g functions. The 1s core orbitals were not correlated. The potential energy functions were obtained by fitting polynomials of eighth degree to 9 or 10 energies. Varying the degree of the polynomials or adding further points had virtually no effect on the computed spectroscopic constants. The calculated spectroscopic constants are summarized in table 3. Generally, QCISD gives best agreement with experiment. The equilibrium distances are too small by 0.001 to 0.004 A, and the harmonic constants are too large by 30-60 cm-‘. An exception is the Fz molecule, in which case the error in r, is unusually large (0.02 A), and w, is about 100 cm-’ too large. The BCCD approximation gives the poorest agreement with experiment in all cases. The CCSD results are between QCISD and BCCD. This is in full accord with the comparisons to FCI data in ref. [ 12 1, which showed the best agreement between

28 February 1992

FCI and QCISD and largest deviations between FCI and BCCD. Also the fact that the o, values are too large is consistent with the findings of ref. [ 121, which demonstrated that all three methods yield energies above FCI and that the errors increase with increasing bond distance. Somewhat surprisingly, there are quite large variations in the differences between the three methods. While for F2 the o, values of BCCD and QCISD differ by only 10 cm-‘, they differ by almost 50 cm-’ for CO. In fact, the BCCD results show somewhat more uniform errors than the QCISD ones. For CO, Fz, and HF the potential energy functions have also been computed using the larger and more flexible basis set B. Only very small effects on the spectroscopic constants are found for FZ. However, for the strongly polar HF the o, values are reduced by about 18 cm-’ (QCISD: we=4175 cm-‘), and the equilibrium distances are slightly increased (QCISD: r, = 0.9 16 8, ). For CO the w, values are reduced by about 5 cm-’ (QCISD: w,=2198.4 cm-‘). We believe that the remaining basis sets errors are small as compared to the inaccuracies of the correlation treatments. Table 3 also includes a comparison with CEPA-1 calculations. This method [ 40,411, which can be viewed as an approximation to CCSD and needs the same computational effort as CISD, yields very similar results as the QCISD method, except in the case of FZ. In this case, CEPA- 1 gives considerably better results. The origin of this effect is not known. The main advantage of the CCSD or QCISD approximations appears to be their invariance with respect to unitary transformations of the internal orbitals, which is only approximately fulfilled in the CEPA- 1 approximation. The internally contracted CI calculations [ 27,341 carried out with the same basis set are more accurate than any of the single-reference calculations of the present work. In contrast to the single-reference methods, MRCI yielded somewhat too large equilibrium distances, but in all cases the o, values differed by less than 10 cm- ’ from the experimental data. The relatively large errors of the QCISD and CCSD results can be much reduced by taking into account connected triple excitations, which have been shown to significantly improve the agreement with FCI energies [ 5,6,12,13], in particular at extended bond 9

Volume 190, number Table 3 Spectroscopic

constants

I,2

CHEMICAL

‘) of diatomic

PHYSICS

28 February

Method

Energy b,

r,

N2

SCF BCCD CCSD QCISD CEPA- 1 MERCI ‘) exp. ‘) CCSD r’ CCSDT-1 r)

-

1.065 1.092 1.093 1.094 1.095 1.101 1.098 1.095 1.103

SCF BCCD CCSD QCISD CEPA- 1 MRCI ‘) exp. *) CCSD ‘) CCSDT- 1 r)

- 112.786327 -113.165262 - 113.166509 - 113.168399 - 113.171077 - 113.168938

F2

SCF BCCD CCSD QCISD CEPA- 1 MRCI g, exp. *)

HF

SCF BCCD CCSD QCISD CEPA- 1 MRCI =’ exp. *) CCSD f, CCSDT- 1 r)

B,

o,

W,

W&

2.122 2.020 2.016 2.012 2.010 1.987 1.998

0.0137 0.0158 0.0160 0.0162 0.0164 0.0171 0.0173

2732.1 2449.6 2438.6 2423.6 2416.5 2350.7 2358.6 2411 2323

10.8 12.7 12.9 13.2 13.4 13.9 14.3

1.102 1.123 1.124 1.127 1.127 1.132 1.128 1.127 1.136

2.025 1.950 1.946 1.938 1.937 1.919 1.931

0.0151 0.0160 0.0163 0.0172 0.0171 0.0173 0.0175

2429.1 2251.2 2236.5 2203.4 2204.9 2164.8 2169.8 2218 2112

11.3 11.8 12.1 13.4 13.2 13.1 13.3

- 198.764231 - 199.330455 - 199.331492 - 199.332601 -199.338189 - 199.323105

1.326 1.388 1.390 1.391 1.406 1.412 1.412

1.009 0.921 0.919 0.917 0.898 0.891 0.890

0.0080 0.0104 0.0104 0.0104 0.0122 0.0127 0.0138

1266.8 1023.8 1019.0 1014.5 941.6 921.2 916.6

6.5 9.2 9.2 9.2 10.8 11.6 11.2

-

0.897 0.914 0.914 0.915 0.915 0.918 0.917 0.917 0.919

0.747 0.776 0.778 0.780 0.785 0.788 0.798

4477.7 4207.6 4201.1 4192.5 4187.5 4138.2 4138.3 4187 4145

84.5 88.0 88.5 89.3 90.3 88.8 89.9

108.989919 109.379993 109.381034 109.382277 109.385644 109.386296

- 109.358131 - 109.376224

-113.141146 -113.158632

100.063913 100.362082 100.362466 100.362949 100.363168 100.370310

-100.312486 - 100.318362

21.90 21.10 21.08 21.06 21.06 20.92 20.96

‘) Basis A, see text. The constants have been computed from polynomial fits of 8th degree to 9 calculated b, Energies at calculated equilibrium distances. ‘) Data from ref. [ 341. ‘) Experimental values from ref. [ 381. ‘) CASSCF reference function with active space (2~So, lx, 2r). ‘) Results from ref. [ 391 for TZ+ 2P basis. *) CASSCF reference function with active space (3o,, 3o., lx,, In,, 2x.).

distances. It has been demonstrated that the accurate prediction of dipole moments and infrared intensities also requires the inclusion of triple excitations [ 42,43 1. Unfortunately, even in a non-iterative perturbational treatment, e.g., CCSD(T) [ 441, the triple cannot easily be computed directly from the A0 integrals. Therefore, an integral-direct CCSD (T) 10

1992

molecules

Molecule

co

LETTERS

energies.

method is not straightforward to implement. Scuseria and Schaefer [ 39,451 have carried out similar calculations for HF, N2, and CO in order to compare the reliability of predictions made with the CCSD, CCSDT, and other methods. In ref. [45] a standard DZP basis set was employed, and spectroscopic constants obtained with the CCSD method

Volume 190, number 1,2

CHEMICAL PHYSICS LETTERS

were found to be in much closer agreement with experiment than in the present work. Clearly, this is due to a cancellation of errors in the basis set and in the correlation treatment. Part of this error was removed in ref. [ 391, where a TZ+ 2P basis was used. The TZ+ 2P results for CCSD and CCSDT-1 of Scuseria and Schaefer are reproduced in table 3. It is seen that the o, values obtained with CCSD and the TZ+ 2P basis are still in error by 15-20 cm-’ as compared to our values obtained with basis B, which can be assumed to be close to the basis set limit. The corresponding r, values in ref. [ 391 are 0.002 A too long. The fact that these errors are of the same magnitude as the differences between the various methods supports our view that very large basis sets should be used when testing the reliability of correlation methods by comparison with experiment. The calculations in refs. [ 39,451 show that inclusion of connected triple excitations reduces the w, values, and in most cases the CCSD (T) method [ 441 would give excellent agreement with experiment when converged basis sets are used [ 39,45-471.

6. Conclusions It has been demonstrated that quadratic configuration interaction (QCISD), coupled cluster (CCSD), and Brueckner coupled-cluster (BCCD) calculations can be carried out without a full integral transformation. This is particularly important for an efficient implementation of the BCCD method and a prerequisite for integral-direct methods. A simple DIIS extrapolation technique is used to improve convergence in BCCD calculations. In contrast to previous implementations, the computational cost of our BCCD method is not higher than for CCSD calculations. It has been demonstrated that CCSD, QCISD, and BCCD calculations with about 250 basis functions and 50 correlated electrons are possible, the only bottleneck being the storage of the twoelectron integrals. The development of integral-direct methods, which will eliminate this problem, is in progress. The accuracy of spectroscopic properties computed with the CCSD, QCISD, and BCCD methods has been compared for some diatomic molecules using very large basis sets. It is found that the QCISD approximation yields the smallest errors, but

28 February 1992

satisfactory agreement with experiment is obtained with none of these methods without the inclusion of connected triple excitations.

Acknowledgement

We like to thank G. Scuseria for providing the test calculations mentioned in section 4. This work was supported by the Deutsche Forschungsgemeinschaft (SFB2 16 ) and the German Fonds der Chemischen Industrie. Computer time on a CRAY-YMP 8/32 was provided by the Hochstleistungsrechenzentrum (HLRZ) Jtilich.

References [ 1] J. Ciiek, J. Chem. Phys. 45 ( 1966) 4256. [2] J. Ciiek, Advan. Chem. Phys. 14 (1969) 35. [ 31 J. Ciiek and J. Paldus, Intern. J. Quantum Chem. 5 ( 1971) 359. , .. .. [4] G.D. Purvis and R.J. Bartlett, J. Chem. Phys. 76 (1982) 1910. [ 51 J. Noga and R.J. Bartlett, J. Chem. Whys. 86 ( 1987) 7041. [6] G.E. Scuseria and H.F. Schaefer III, Chem. Phys. Letters 152 (1988) 382. [ 7 ] M. Head-Gordon, J.A. Pople and K. Raghavachati, J. Chem. Phys. 87 (1987) 5968. [8]J.CiiekandJ.Paldus,Ph&aScripta21 (1980)251. [9] R.J. Bartlett and G.D. Purvis, Physica Scripta 21 (1980) 255. [ lo] R.A. Chiles and C.E. Dykstra, J. Chem. Phys. 74 (1981) 4544. [ 111 R.J. Bartlett, C.E. Dykstra and J. Paldus, NATO ASI, Vol. 133 (Reidel, Dordrecht, 1984) p. 127. 121 N.C. Handy, J.A. Pople, M. Head-Gordon, K. Raghavachari and G.W. Trucks, Chem. Phys. Letters 164 ( 1989) 185. 13 ] K. Raghavachari, J.A. Pople, E.S. Replogle, M. HeadGordon and N.C. Handy, Chem. Phys. Letters 167 ( 1990) 115. 141 G.E. Scuseria and H.F. Schaefer III, Chem. Phys. Letters 142 (1987) 354. 151 J. Almlbf, K. Frergi and K. Korsell, J. Comput. Chem. 3 (1982) 385. 161 M. Head-Gordon, J.A. Pople and M.J. Frisch, Chem. Phys. Letters 153 (1988) 503. 17) S. Sreo and J. Almliif, Chem. Phys. Letters 154 (1989) 83. 181 R. Ahlrichs, M. H&r, M. Bar, H. Horn and C. Kiilmel, Chem. Phys. Letters 162 ( 1989) 165. [ 191 H.-J. Werner, P.J. Knowles and J. Almliif, to be published. [20] N.C. Handy, work presented at the 27th Symposium flir Theoretische Chemie, Bielefeld, September 199 1. [21] P. Pulay, Chem. Phys. Letters 73 (1980) 393.

11

Volume 190, number I,2

CHEMICAL PHYSICS LETTERS

[22] W. Meyer, J. Chem. Phys. 64 (1976) 2901. [23] P. Pulay, S. Sz.berand W. Meyer, J. Chem. Phys. 81 ( 1984) 1901. [ 241 G.E. Scuseria, C.L. Janssen and H.F. Schaefer III, J. Chem. Phys. 89 (1988) 7382. [25] J.A. Pople, R. Seeger and R. Krishnan, Intern. J. Quantum Chem. Symp. 11 ( 1978) 149. [26] H.-J. Werner and E.A. Reinsch, J. Chem. Phys. 76 ( 1982) 3144. [27] H.-J. Werner and P.J. Knowles, J. Chem. Phys. 89 (1988) 5803. [ 28 ] T.J. Lee and J.E. Rice, Chem. Phys. Letters 150 ( 1988) 406. [ 291 G.E. Scuseria, T.J. Lee and H.F. Schaefer III, Chem. Phys. Letters 130 (1986) 236. [30] T.H. Dunning, J. Chem. Phys. 90 (1989) 1007. [ 3 1 ] E.R. Davidson, J. Comput. Phys. 17 ( 1975) 87. [ 321 G.E. Scuseria, private communication. [33] G.E. Scuseria, A.C. Scheiner, T.J. Lee, J.E. Rice and H.F. Schaefer III, J. Chem. Phys. 86 ( 1986) 288 1. [34] H.-J. Werner and P.J. Knowles, Theoret. Chim. Acta 78 (1990) 175.

12

28 February 1992

[ 351 H.-J. Werner and P.J. Knowles, J. Chem. Phys. 94 ( 199 1) 1264. [ 361 F.B. van Duijneveldt, IBM Research Report RJ 945 ( 1971). [37] P.J. Knowles, K. Stark and H.-J. Werner, Chem. Phys. Letters 185 (1991) 555. [38] K.P. Huber and G. Herzberg, Constants of diatomic molecules (Van Nostrand Reinhold, New York, 1979). [39] GE. Scuseria and H.F. Schaefer III, Chem. Phys. Letters 148 (1988) 205. [ 401 W. Meyer, Inter. J. Quantum Chem. Symp. 5 ( 197 1) 34 1. [41] W. Meyer, J. Chem. Phys. 58 (1973) 1017. [ 421 V. Kellij, J. Noga, G.H.F. Diercksen and A.J. Sadley, Chem. Phys. Letters 152 (1988) 387. [43] G.E. Scuseria, M.D. Miller, F. Jensen and J. Geertsen, J. Chem. Phys. 94 ( 1991) 6660. [44] K. Raghavachari, G.W. Trucks, J.A. Pople and M. HeadGordon, Chem. Phys. Letters 157 ( 1989) 479. [45] G.F. Scuseria and H.F. Schaefer III, Chem. Phys. Letters 152 (1988) 382. [46] T.J. Lee and G.E. Scuseria, J. Chem. Phys. 93 (1990) 489. [47] G.F. ScuseriaandT.J. Lee, J. Chem. Phys. 93 (1990) 5851.