- Email: [email protected]
Convergence of a matrix-oriented coupled cluster (doubles) method
Volume 145, number 6
CHEMICAL PHYSICS LETTERS
22 April 1988
CONVERGENCE OF A MATRIX-ORIENTED COUPLED CLUSTER (DOUBLES) METHOD Clifford E. DYKSTRA and Joseph D. AUGSPURGER Department of Chemistry, University ofIllinois, Urbana, IL 61801, USA
Received 29 January 1988
The convergence of an iterative matrix-oriented coupled cluster approach is investigated and compared with a more recent implementation. Where the reference wavefunction dominates, the special features of the matrix-oriented implementation result in rapid convergence.
The recent paper by Scuseria, Lee and Schaefer (SLS ) [ 11 discussed their implementation of a sin-
gles and doubles coupled cluster (CCSD) treatment for electron correlation. It showed one way, the direct inversion of the iterative subspace (DIIS) method of Pulay [ 2 1, in which rapid convergence could be achieved within the coupled cluster approach. Trucks, Noga and Bartlett [ 31 report that this implementation is equivalent to the reduced linear equation (RLE) method of Purvis and Bartlett [ 41 introduced several years ago. There are, though, other effective means, including developments with matrix-oriented methods [ 51. SLS offered no comparison with these, yet their effectiveness reinforces the idea that coupled cluster methods are computationally workable, with several fast convergence options. The method that SLS implemented for CC wavefunctions is designed for extrapolation in the iterative solution of a system of coupled equations. It requires storing correction vectors and SLS adopted eight as the number of CCSD expansion coefficient (amplitude) vectors to keep. They found that starting extrapolation between the fifth and eighth CCSD iteration and continuing it for two or three iterations produced satisfactory convergence. Our earlier implementation of the coupled cluster approach [5 ] not only achieves good convergence in most cases, but also can take advantage of other computational features. This simple procedure is an iterative first-order correction to the correlated wavefunction. Computed on each interaction is a matrix element be-
tween each configuration being included in the wavefunction expansion and the current guess of the wavefunction, Y. If t, stands for the coupled cluster amplitude of the vf substituted configuration, then the iterative cycle for finding the wavefunction is, @+I)
=t’“’ I
_
(Y,;IH-E’“‘I
!P)
(1)
This requires saving only one set of amplitudes, {tl}, at a time, which is helpful for large expansions. It is most rapidly convergent when Y/en)dominates. Initially, ylco)= t&F. This correction cycle can be used with any chosen set of substitutions, doubles, singles and doubles, singles, doubles and triples, etc. The value of ( yI 1H( u/l) in the denominator in eq. ( 1) may be approximated since the numerator goes to zero upon convergence. For this reason, a freely chosen parameter, d, is inserted to aid convergence. The best choice is a value similar in size to the final total correlation energy since at the outset, the correlation energy (in E(O)) is zero. In practice, ( w,IH( yr> is conveniently obtained just from modified orbital energies, and A is chosen to be a value slightly larger than the expected correlation energy. In most cases, there is little difference in convergence upon changing A by the lesser of 0.05 au or + 20%. Another useful part of our iterative treatment is that virtual orbitals sets used to construct substituted configurations are chosen for each electron pair so as to minimize the convergence-retarding coupling between doubly substituted configurations. The use of
0 009-2614/88/$ 03.50 0 Elsevier Science Publishers B.V. ( North-Holland Physics Publishing Division )
545
Volume145,number6
CHEMICAL PHYSICSLETTERS
Table1
effective. Our iteration cycles were at the CCD level and those of SLS were at the CCSD level (with singles in the cluster operator), but this should make little difference in convergence rate. As was demonstrated in the course of a CIDS study by Schaefer and co-workers [ 8 1, there is very little difference when relaxing the double substitution coefficients because of the inclusion of the single substitutions versus not relaxing them (i.e. not including singles). CI and CC approaches should be similar on this point, and the work of Purvis and Bartlett has amply demonstrated the negligible impact of the single substitutions on the convergence rate for a coupled cluster wavefunction [ 41. Trucks et al. [ 31 also make this point in their report. Our procedure is not faster than that of SLS or Bartlett’s approach [ 3 1, for the H,Be system, which is basically at a transition state structure where WseFdoes not strongly dominate. For such cases, extrapolation via the RLE procedure is preferred. However, since some number of initial iterations are required to extrapolate, perhaps the overall best process is one that starts as ours and then extrapolates. As an aside, we point out that the treatment of single substitutions has been done in two ways in our CC implementation. The most complete of these is a full CCSD treatment but with the orbitals fully relaxed so that the pi (singles) amplitudes or expansion coefficients converge to values that are identically zero [ 5 1. This is a step beyond what is normally designated CCDS. The other way is to retain the original occupied orbitals and optimally incorporate the single substitutions onto the CCD
CCD convergence example: Hz0 with a DZ basis following Scuseria et al. [ 1 ] Iteration a)
Total energy (au)
Energy change (au)
SCF 1 2 3 4 5 6 7 8 CCD+S
-76.14358708 -76.15451921 -76.15523150 -76.15526708 -76.15527341 -76.15527285 -76.15527297 -76.15527290 -76.156 18155
-75.00983770 -0.13374739 -0.010932 13 -0.00071228 -0.000035 59 -0.00000633 -0.00000056 -0.00000012 -0.00000007 -0.00090865
22April1988
a) 4~0.24 (eq. (1)).
pair-selected virtual orbitals is unique to approaches that represent the wavefunction in terms of atomic basis functions, and Chiles and Dykstra’s treatment [ 5 ] was the first CCD method of that sort. This convergence feature and the iterative process of eq. ( 1) were originally used for the CID method called selfconsistent electron pair (SCEP) theory, which was developed by Meyer [6]. Dykstra, Schaefer and Meyer [ 71 showed in 1976 how these features make for good convergence for SCEP (CID) wavefunctions. Table 1 illustrates convergence for a specific case using matrix-oriented CCD: Table 2 shows convergence data using most of the examples reported by SLS [ 11. Our procedure converges more rapidly than the procedure of SLS in several cases, though it may be fair to regard the two as comparable and both very
Table 2 Convergence of iterative (eq. ( 1) ) matrix-oriented CCD+S versus convergence of CCSD using DIIS ‘) Molecule
Basis
Number of basis functions
HB H,Be Hz0
DZ DZ DZ
12 14 14
H20
DZP
26
db’
0.08 0.12 0.24 0.21 0.24
Iterations required to convergence to 10e6 au
lo-’ au
lo-* au
6 14 6 6 5
7 20 8 8 8
8 28 9 9 10
(SLS: 11) (SLS: 14) (SLS: 11) (SLS: 10)
‘) Basis sets and geometries are the same as used by SLS [ 11. However, their specification of the DZP basis for water leads to slightly different SCF and correlation energies than they reported. b’A (eq. ( 1) ) is chosen as a rough guess of the correlation energy. The calculations on Hz0 with two different choices of d are representative of the sensitivity to this value.
546
Volume 145, number 6
CHEMICAL PHYSICS LETTERS
wavefunction. What this means is that upon convergence of a CCD calculation, the small CI matrix formed from the usual set of singly substituted configurations plus the CCD wavefunction (as if it were one configuration) is diagonalized. The doubly substituted configuration expansion coefficients are only indirectly relaxed in this procedure (through renormalization), and the singles correlation energy is usually quite close to that obtained with CCSD. This process, which we may designate as CCD + S is quite useful because the computational cost of the “+S” step is entirely unimportant relative to the CCD process, whereas adding the single substitutions into the CC exponential operator does add a noticeable number of computational steps. In the example in table 1, the singles correlation energy obtained with CCD+S is 0.000909 au, while SLS obtained 0.000803 au with CCSD. A final point about computational efficiency within the matrix-oriented coupled cluster approach is the important savings that can be achieved by tint carrying out the full coupled cluster treatment on each iteration. Dykstra has shown that CCD calculations can converge quickly by a simple scaling of CID expansion coefficients [ 9 1. Alternatively, the approximate form of CCD, ACCD [ lo], which involves many less computational steps, can be used for the
22 April 1988
first several iterations. Tests that we have carried out show that this does not retard convergence, but does yield a net computational savings. Either way, significant computational savings can be realized by reserving the full CC treatment for only the last several iterations. This work was supported, in part, by the Physical Chemistry Program of the National Science Foundation (Grant CHE-19496). References [ 1] G.E. Scuseria, T.J. Lee and H.F. Schaefer III, Chem. Phys. Letters 130 (1986) 236. [2] P. Pulay, J. Comput. Chem. 3 (1982) 556. [ 31 G.W. Trucks, J. Noga and R.J. Bartlett, Chem. Phys. Letters 145 (1988) 548. [4] G.D. Purvis and R.J. Bartlett, J. Chem. Phys. 75 (1981) 1284. [ 51 R.A. Chiles and C.E. Dykstra, J. Chem. Phys. 74 (1981) 4544. [6] W. Meyer, J. Chem. Phys. 64 (1976) 2901. [7] C.E. Dykstra, H.F. Schaefer III and W. Meyer, J. Chem. Phys. 65 (1976) 2740. [ 81 C.E. Dykstra, M. Hereld, R.R. Lucchese, H.F. Schaefer III and W. Meyer, J. Chem. Phys. 67 (1977) 4071. [9] C.E. Dykstra, Chem. Phys. Letters 88 (1982) 202. [ lo] SM. Bachrach, R.A. Chiles and C.E. Dykstra, J. Chem. Phys. 75 (1981) 2270.
547
CHEMICAL PHYSICS LETTERS
22 April 1988
CONVERGENCE OF A MATRIX-ORIENTED COUPLED CLUSTER (DOUBLES) METHOD Clifford E. DYKSTRA and Joseph D. AUGSPURGER Department of Chemistry, University ofIllinois, Urbana, IL 61801, USA
Received 29 January 1988
The convergence of an iterative matrix-oriented coupled cluster approach is investigated and compared with a more recent implementation. Where the reference wavefunction dominates, the special features of the matrix-oriented implementation result in rapid convergence.
The recent paper by Scuseria, Lee and Schaefer (SLS ) [ 11 discussed their implementation of a sin-
gles and doubles coupled cluster (CCSD) treatment for electron correlation. It showed one way, the direct inversion of the iterative subspace (DIIS) method of Pulay [ 2 1, in which rapid convergence could be achieved within the coupled cluster approach. Trucks, Noga and Bartlett [ 31 report that this implementation is equivalent to the reduced linear equation (RLE) method of Purvis and Bartlett [ 41 introduced several years ago. There are, though, other effective means, including developments with matrix-oriented methods [ 51. SLS offered no comparison with these, yet their effectiveness reinforces the idea that coupled cluster methods are computationally workable, with several fast convergence options. The method that SLS implemented for CC wavefunctions is designed for extrapolation in the iterative solution of a system of coupled equations. It requires storing correction vectors and SLS adopted eight as the number of CCSD expansion coefficient (amplitude) vectors to keep. They found that starting extrapolation between the fifth and eighth CCSD iteration and continuing it for two or three iterations produced satisfactory convergence. Our earlier implementation of the coupled cluster approach [5 ] not only achieves good convergence in most cases, but also can take advantage of other computational features. This simple procedure is an iterative first-order correction to the correlated wavefunction. Computed on each interaction is a matrix element be-
tween each configuration being included in the wavefunction expansion and the current guess of the wavefunction, Y. If t, stands for the coupled cluster amplitude of the vf substituted configuration, then the iterative cycle for finding the wavefunction is, @+I)
=t’“’ I
_
(Y,;IH-E’“‘I
!P)
(1)
This requires saving only one set of amplitudes, {tl}, at a time, which is helpful for large expansions. It is most rapidly convergent when Y/en)dominates. Initially, ylco)= t&F. This correction cycle can be used with any chosen set of substitutions, doubles, singles and doubles, singles, doubles and triples, etc. The value of ( yI 1H( u/l) in the denominator in eq. ( 1) may be approximated since the numerator goes to zero upon convergence. For this reason, a freely chosen parameter, d, is inserted to aid convergence. The best choice is a value similar in size to the final total correlation energy since at the outset, the correlation energy (in E(O)) is zero. In practice, ( w,IH( yr> is conveniently obtained just from modified orbital energies, and A is chosen to be a value slightly larger than the expected correlation energy. In most cases, there is little difference in convergence upon changing A by the lesser of 0.05 au or + 20%. Another useful part of our iterative treatment is that virtual orbitals sets used to construct substituted configurations are chosen for each electron pair so as to minimize the convergence-retarding coupling between doubly substituted configurations. The use of
0 009-2614/88/$ 03.50 0 Elsevier Science Publishers B.V. ( North-Holland Physics Publishing Division )
545
Volume145,number6
CHEMICAL PHYSICSLETTERS
Table1
effective. Our iteration cycles were at the CCD level and those of SLS were at the CCSD level (with singles in the cluster operator), but this should make little difference in convergence rate. As was demonstrated in the course of a CIDS study by Schaefer and co-workers [ 8 1, there is very little difference when relaxing the double substitution coefficients because of the inclusion of the single substitutions versus not relaxing them (i.e. not including singles). CI and CC approaches should be similar on this point, and the work of Purvis and Bartlett has amply demonstrated the negligible impact of the single substitutions on the convergence rate for a coupled cluster wavefunction [ 41. Trucks et al. [ 31 also make this point in their report. Our procedure is not faster than that of SLS or Bartlett’s approach [ 3 1, for the H,Be system, which is basically at a transition state structure where WseFdoes not strongly dominate. For such cases, extrapolation via the RLE procedure is preferred. However, since some number of initial iterations are required to extrapolate, perhaps the overall best process is one that starts as ours and then extrapolates. As an aside, we point out that the treatment of single substitutions has been done in two ways in our CC implementation. The most complete of these is a full CCSD treatment but with the orbitals fully relaxed so that the pi (singles) amplitudes or expansion coefficients converge to values that are identically zero [ 5 1. This is a step beyond what is normally designated CCDS. The other way is to retain the original occupied orbitals and optimally incorporate the single substitutions onto the CCD
CCD convergence example: Hz0 with a DZ basis following Scuseria et al. [ 1 ] Iteration a)
Total energy (au)
Energy change (au)
SCF 1 2 3 4 5 6 7 8 CCD+S
-76.14358708 -76.15451921 -76.15523150 -76.15526708 -76.15527341 -76.15527285 -76.15527297 -76.15527290 -76.156 18155
-75.00983770 -0.13374739 -0.010932 13 -0.00071228 -0.000035 59 -0.00000633 -0.00000056 -0.00000012 -0.00000007 -0.00090865
22April1988
a) 4~0.24 (eq. (1)).
pair-selected virtual orbitals is unique to approaches that represent the wavefunction in terms of atomic basis functions, and Chiles and Dykstra’s treatment [ 5 ] was the first CCD method of that sort. This convergence feature and the iterative process of eq. ( 1) were originally used for the CID method called selfconsistent electron pair (SCEP) theory, which was developed by Meyer [6]. Dykstra, Schaefer and Meyer [ 71 showed in 1976 how these features make for good convergence for SCEP (CID) wavefunctions. Table 1 illustrates convergence for a specific case using matrix-oriented CCD: Table 2 shows convergence data using most of the examples reported by SLS [ 11. Our procedure converges more rapidly than the procedure of SLS in several cases, though it may be fair to regard the two as comparable and both very
Table 2 Convergence of iterative (eq. ( 1) ) matrix-oriented CCD+S versus convergence of CCSD using DIIS ‘) Molecule
Basis
Number of basis functions
HB H,Be Hz0
DZ DZ DZ
12 14 14
H20
DZP
26
db’
0.08 0.12 0.24 0.21 0.24
Iterations required to convergence to 10e6 au
lo-’ au
lo-* au
6 14 6 6 5
7 20 8 8 8
8 28 9 9 10
(SLS: 11) (SLS: 14) (SLS: 11) (SLS: 10)
‘) Basis sets and geometries are the same as used by SLS [ 11. However, their specification of the DZP basis for water leads to slightly different SCF and correlation energies than they reported. b’A (eq. ( 1) ) is chosen as a rough guess of the correlation energy. The calculations on Hz0 with two different choices of d are representative of the sensitivity to this value.
546
Volume 145, number 6
CHEMICAL PHYSICS LETTERS
wavefunction. What this means is that upon convergence of a CCD calculation, the small CI matrix formed from the usual set of singly substituted configurations plus the CCD wavefunction (as if it were one configuration) is diagonalized. The doubly substituted configuration expansion coefficients are only indirectly relaxed in this procedure (through renormalization), and the singles correlation energy is usually quite close to that obtained with CCSD. This process, which we may designate as CCD + S is quite useful because the computational cost of the “+S” step is entirely unimportant relative to the CCD process, whereas adding the single substitutions into the CC exponential operator does add a noticeable number of computational steps. In the example in table 1, the singles correlation energy obtained with CCD+S is 0.000909 au, while SLS obtained 0.000803 au with CCSD. A final point about computational efficiency within the matrix-oriented coupled cluster approach is the important savings that can be achieved by tint carrying out the full coupled cluster treatment on each iteration. Dykstra has shown that CCD calculations can converge quickly by a simple scaling of CID expansion coefficients [ 9 1. Alternatively, the approximate form of CCD, ACCD [ lo], which involves many less computational steps, can be used for the
22 April 1988
first several iterations. Tests that we have carried out show that this does not retard convergence, but does yield a net computational savings. Either way, significant computational savings can be realized by reserving the full CC treatment for only the last several iterations. This work was supported, in part, by the Physical Chemistry Program of the National Science Foundation (Grant CHE-19496). References [ 1] G.E. Scuseria, T.J. Lee and H.F. Schaefer III, Chem. Phys. Letters 130 (1986) 236. [2] P. Pulay, J. Comput. Chem. 3 (1982) 556. [ 31 G.W. Trucks, J. Noga and R.J. Bartlett, Chem. Phys. Letters 145 (1988) 548. [4] G.D. Purvis and R.J. Bartlett, J. Chem. Phys. 75 (1981) 1284. [ 51 R.A. Chiles and C.E. Dykstra, J. Chem. Phys. 74 (1981) 4544. [6] W. Meyer, J. Chem. Phys. 64 (1976) 2901. [7] C.E. Dykstra, H.F. Schaefer III and W. Meyer, J. Chem. Phys. 65 (1976) 2740. [ 81 C.E. Dykstra, M. Hereld, R.R. Lucchese, H.F. Schaefer III and W. Meyer, J. Chem. Phys. 67 (1977) 4071. [9] C.E. Dykstra, Chem. Phys. Letters 88 (1982) 202. [ lo] SM. Bachrach, R.A. Chiles and C.E. Dykstra, J. Chem. Phys. 75 (1981) 2270.
547