J. Phys. Chem. 1990,94, 5411-5482 It is clear that with the successful modeling of the copper atoms as one-electron systems it is not too surprising to find that the results for the copper clusters are often qualitatively similar to results for alkali-metal clusters. This fact has been noted several times in both experimental and theoretical papers. The most direct comparisons with the present copper clusters can be made with the corresponding anionic lithium clusters, which have been studied theoretically by Boustani and K o ~ t e c k p . ' Apart ~ from the similarities it can be noted that the smaller Li, clusters up to n = 5 were all found to be linear. This could be a real difference between lithium and copper clusters but could also be an effect of different treatments of core correlation. In our study we have noted that nonlinear structures were strongly stabilized in comparison to linear structures after core correlation effects were added and such effects were not included for Li; in ref 15. In this study we have shown that the present simplified treatment of transition-metal atoms can give useful and reliable information about cluster properties. In terms of reactivities copper clusters are not the most interesting metal clusters. For example, no copper cluster is reactive toward Hl, while for nickel clusters the opposite is true and all clusters are reactive. Iron and cobalt clusters are more interesting in that the reactivities vary consid-
5477
erably from cluster to cluster. In a recent paper' we have suggested that the oscillations in reactivity for these clusters are due to a tendency for iron and cobalt to mix in d"s2 atoms among the normal d*'s atoms and thereby cause some of the clusters to have a closed valence shell character. The closed-shell character will lead to a drastically reduced reactivity as compared to the corresponding nickel clusters. It seems reasonable to expect that the formation of closed valence shells for iron and cobalt clusters should be seen in some spectral properties, such as very low electron affinities for these clusters. The same type of experimental information as has been obtained for the copper clusters in terms of EDE's and HOMO-LUMO gaps would thus be even more interesting if it could be obtained for iron and cobalt clusters. It is likely that such information together with theoretical calculations could definitely establish the origin of the oscillations in reactivity that were observed already several years ago and for which conflicting explanations have been proposed. R d ~ t r yNO. C U ~107865-35-0; , Cuj-, 108658-50-0; C U ~108658, 51-1; C U ~108658-52-2; , C U ~108658-53-3; , CUT,108658-54-4; C U ~ , 108658-55-5; Cb-, 108658-56-6; CU~O-, 108658-57-7; CUI,12190-70-4; CU,, 66771-03-7; Cud, 65357-62-2; CU,, 66771-05-9; C U ~681 , 12-39-0; CUT,81 543-63-7; C U ~62863-28-9; , Cup, 81 543-64-8; CUI&127398-94-1.
The Restrlcted Active Space Self-Consistent-Field Method, Implemented with a Split Graph Unitary Group Approach Per-brke Malmqvist,* Alistair Rendell2 and Bjorn 0. Roos Department of Theoretical Chemistry, Chemical Centre, P.O.B. 124, S-221 00 Lund, Sweden (Received: November 1 , 1989; In Final Form: February 20, 1990)
An MCSCF method based on a restricted active space (RAS) type wave function has been implemented. The RAS concept is an extension of the complete active space (CAS) formalism, where the active orbitals are partitioned into three subspaces: RASI, which contains up to a given maximum number of holes; RAS2, where all possible distributionsof electrons are allowed; and RAS3, which contains up to a given maximum number of electrons. A typical example of a RAS wave function is all single, double, etc. excitations with a CAS reference space. Spin-adapted configurations are used as the basis for the MC expansion. Vector coupling coefficients are generated by using a split graph variant of the unitary group approach, with the result that very large MC expansions (up to around 1@) can be treated, still keeping the number of coefficients explicitly computed small enough to be kept in central storage. The largest MCSCF calculation that has so far been performed with the new program is for a CAS wave function containing 716418 CSF's. The MCSCF optimization is performed using an approximate form of the super-CI method. A quasi-Newton update method is used as a convergence accelerator and results in much improved convergence in most applications.
Introduction The complete active space (CAS) SCF method'+ has during the past 10 years been used in a large number of applications covering many areas of applied quantum chemistry.$ The pop ularity of the CASSCF approach is mostly due to its simple structure, where the only effort from the user is to select an adequate set of active orbitals in order to defme the wave function. In most applications this is a straightforward procedure, even if there exist cases where a deeper insight into the nature of the wave function is necessary before a proper choice can be made. The major drawback of the CASSCF approach is that the number of configuration state functions (CSFs) increases steeply with the number of active orbitals. For example, distributing 12 electrons among 12 active orbitals yields 226 5 12 singlet CSF's, if no reduction due to spatial symmetry is assumed. Adding one more active orbital increases the number to 7 3 6 164. This strong dependence on the size of the MC expansion on the number of active orbitals makes it difficult to perform CASSCF calculations Present address: NASA Am- Research Center, Computational Chemistry, Moffett Field, CA 94035.
with more than around 12 active orbitals, at least in cases where the number of active electrons is about the same as the number of active orbitals. In most applications this limit does not pose a problem, but there are situations where it would be desirable to be able to increase the number of correlating orbitals in the MCSCF wave function. Normally dynamic correlation effects are included in the calculation in a subsequent step, where a multireference (MR) CI calculation is performed using selected CSF's from the CASSCF wave function as reference configurations. The CASSCF calculation is then restricted to include only static correlation effects, and is used to determine the orbitals. (1) Roos, B. 0.;Taylor, P. R.;Siegbahn, P. E. M. Chem. Phys. 1980,48, 151. ( 2 ) Siegbahn, P. E. M.; Almlbf, J.; Heiberg, A.; Roos, B. 0. J . Chem. Phys. 1981, 74, 2384. (3) R.oos, B. 0. Inr. J. Quantum Chem. Symp. 1980, 14, 175. (4) Siegbahn,P. E. M.; Heiberg. A.; Roos, B. 0.;Levy,B. Phys. Scr. 1980, 21, 323.
(5) Roos, B. 0. In Ab Initio Methods in Quantum Chemistry, Part I&
Adv. Chem. Phys.. LXIX;Lawley, K.P., Ed.;Wiley: Chichester, U.K.,1987.
0 1990 American Chemical Society
Malmqvist et al.
5478 The Journal of Physical Chemistry, Vol. 94, No. 14, 1990
However, cases exist where the shapes of the strongly occupied orbitals depend on dynamical correlation effects. A typical example of this phenomenon is found for negative ions, which in many cases are only bound when a large fraction of the dynamical correlation is taken into account. Orbital optimization at a correlated level is then the only consistent approach. The polarity in a chemical bond is often sensitive to dynamic correlation effects. The calculation of, e.g., a dipole moment may in such cases pose difficulties if the orbitals have been obtained at a too low level of approximation. A third example relates to the calculation of wave functions for excited states of molecules. The shape of orbitals occupied in the excited state may in some cases depend strongly on dynamic correlation effects. A drastic example is given by the V state in the ethene molecule, which has the dominant electronic configuration ( A ~ A ~The ) ~ 7rB . orbital is very diffuse at the H F level of approximation. It contracts slowly when more and more correlation is included into the wave function.6 An extension of the CAS concept to include the types of wave functions needed in the examples given above has recently been suggested.' A restricted active space (RAS) is used, where the active orbitals are partitioned into three subsets: RASl, RAS2, and RAS3. Restrictions are introduced into the RASl and RAS3 subspaces such that in RASl a maximum number of holes are allowed and in RAS3, a maximum number of electrons. The remaining active electrons are distributed in all possible ways among the RAS2 orbitals, which in this respect corresponds to the active orbitals in CAS type wave functions. A core of inactive orbitals is also used as in the CAS wave function. The RAS CI method has previously been implemented using a determinant-based CI e~pansion.~ The present implementation is in terms of spin-adapted CSF's. The CI module has been designed to be able to handle large Full CI (FCI) problems efficiently. While it has been shown that determinant expansions offer the most readily vectorizable implementation, the increase in expansion length, and the fact that the wave function is not inherently a spin eigenstate, remain a disadvantage. The alternative is to look for improvements in the graphical unitary group approach (GUGA)," which allow the CI expansion to increase without increasing the list of coupling coefficients beyond feasible limits. A method has been found that achieves this goal. The basic idea is to split the graph into an upper and a lower part and compute explicit coupling coefficients for the two halves separately (split-GUGA). A major problem with a spin-adapted basis, in the context of FCI, is the generation, storage, and use of a large number of coupling coefficients. In particular, since the number of twoelectron coupling coefficients increases very rapidly, it soon became obvious that the explicit storage of them would seriously limit the size of such calculations. To circumvent this limitation the factorization of the two-electron coupling coefficients into products of one-electron coupling coefficients, which has been inherent in their generation? was incorporated into the FCI algorithm. This can be very efficiently implemented on a vector computerg but still placed a realistic limit of around 50000 CSF's to the size of the FCI due to storage requirements. A further reduction in the size of the external storage required can be achieved by a facorization of the oneelectron coupling coefficients. One possible way of doing this is the ghost orbital method, which has recently been used in the internally contracted CI code of Werner and Knowles.Io The other is a straightforward dissection of the segment product form of the coupling coeffcients, which is already at the heart of the GUGA method. The idea is to keep stored not the entire set of finished products, but a much smaller set of half-finished products to be combined when needed. This approach has already been used for several years in our transition density programs. Combined with a new storage scheme for the C1 ~
~
coefficients, these half-finished products, or new set of coupling coefficients, makes it possible to implement a directly vectorized u vector formation in terms of row or column DAXPY operations on a matrix. The orbital optimization in the RASSCF program has been implemented using an improved version of the approximate super-CI method, which had been used earlier for CASSCF calc u l a t i o n ~ . 'This ~ ~ method has in many applications shown good convergence properties, although the asymptotic convergence is only first order. The virtue of the method lies in its computational efficiency and the fact that the computational effort is almost independent of the number of inactive orbitals. Only a limited set of transformed two-electron integrals is needed, namely those with three of the four indexes corresponding to active orbitals. This method is therefore especially suited to treat large systems with many electrons and large basis sets. The first-order convergence can be improved by the use of quasi-Newton update methods.Il Such a scheme has been implemented here and results in greatly improved convergence properties in most applications. The update is made on the inverted approximate Hessian. The most important effect of the update is to reintroduce the CI coupling into the orbital optimization step. The RASSCF program is presently being used in a number of applications. Calculations on excited states of aromatic systems are done where the RASl and RAS3 spaces are used to introduce dynamic polarization terms into a A-electron CAS space.'2 The size of the A orbital in the V state of ethene is investigated as a function o! dynamic correlation.I2 Other applications of the new code involve studies of the ground and excited states of the CCCN radical, the HC1- ion, the dipole moment function of CrO, and a study of the stability of the CH,-. The update procedures described below have also been implemented into a CASSCF super-CI program and are used in transition-metal-cluster calculation~.~~ The Split Graph Unitary Group Approach The u routine is the main engine of a modern CI program, where most of the CPU time is spent. It is designed to perform a single operation as efficiently as possible: the application of some Hamiltonian operator to a wave function to generate the so-called 0 vector
either in a formally exact form (full CI or CAS), or approximated in some truncated wave function space. By adding some of the two-electron integrals to the one-electron integrals, the remaining two-electron part can be formulated algorithmically as loop over r and s urs :=
i,\k
End of r,s loop
(2)
Again u should be approximated by projection onto some wave function space. The algorithm assumes that the expansion space needed for the temporary wave functions u, is not much larger than that for the wave function and the final u function. This is obviously the case for CAS wave functions with no symmetry. To take symmetry into account merely requires that different symmetry blocks of the same CSF space is used, which is quite practical. It has been shown' that this is just as easy to do in the much more general RAS case, since the key property of the RAS space is closure under deexcitation: whenever the orbital index
R.; Roos, E. 0. Int. J . Quantum Chem. 1989,35, 813. (7) Ofsen, J.; Roos, E.0.; Jsrgensen, P.; Jensen, H. J. Aa. J. Chem. Phys. ( 6 ) Lindh,
1988,89,2185. (8) Shavitt. I. Inr. J. Quantum Chem. Symp. 1977,II,133; 1978,12, 5 ( 9 ) Siegbahn, P. E. M.Chem. Phys. Lett. 1984, 109,417. (IO) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1988,89, 5803
( I I ) Olsen, J.; Yeager, D. L.; Jolrgensen, P. Adv. Chem. Phys. 1983,54, I.
(12) Fuelscher, M.;Malmqvist, P.-A; Roos, B. 0. Work in progress. ( 1 3) Wahlgren, U.Private communication,
Restricted Active Space S C F Method
The Journal of Physical Chemistry, Vol. 94, No. 14, 1990 5419
order is such that may have compnents outside the expansion space, one merely commutes the E, and E,, factors, and the commutator is added to the one-electron part of the problem. An obvious improvement is to use the fact that (pqlrs) = (pqlsr) for the real-valued integrals, by a) skipping loops when r < s, and (b) using En 13, instead of when r > s. In a RAS space, this stratagem cannot be used for those pairs of operators, where the order needs to be interchanged. A basic operation in the CI is thus the formation of l?,$q.We now consider this operation with reference to the further factorization of the coupling coefficients themselves. Assume that the GUGA graph is split into two halves, one upper and one lower. Any spin-coupled configuration state function, or CSF, in the graph, of any symmetry, is normally represented by one specific walk from the top vertex to the bottom vertex. In the new scheme, it is represented by a pair: One upper walk from top vertex to any of a number of midvertices on the midleuel, which is where the graph has been split, and one lower walk from that midvertex down to the bottom vertex. Obviously all CSF‘s can then be represented by positions in a collection of two-dimensional arrays: There is one array for each midvertex, and each upper walk corresponds to one row, while each lower walk corresponds to one column. Point group symmetry is implemented by symmetryblocking the lists of upper and lower walks; then, only specific submatrices can contain nonzero coefficients, and only these submatrices need to be stored. Now cpnsider the calculation of one specific such matrix block of un = En*, for the case that r > s > m, where m is the midlevel. A particular matrix element has the following graph:
Am
+
9
LL
\
LR
1
The left-hand upper walk is denoted here by ULand the right-hand upper walk by UR,and similarly LLand LRdenote the lower walks. Here, LL = LR. According to the usual GUGA rules, this particular loop gives a contribution to the u vector which is where C stands for the expansion coefficients of \k and M denotes the specific matrix block belonging to this midvertex. Only the row index is changed by application of the operator. The column index remains the same, as indicated by the Kronecker symbol. Thus, all the lower walks can be treated at once in a DAXPY operation, using the UL-th row of this matrix block of C as input vector. The coefficient is simply a conventional one-electron coupling coefficient. Similarly, the m + 1 > r > s case is treated as a column DAXPY. The graph for the r > m > s - 1 case is shown below:
LL
y
‘R
where we have reverted to algebraic notation and suppressed the block superscript M. The point is that where previously we used ordinary coupling coefficients, but within a limited range of orbital indexes, we have now split the coupling coefficient into two precalculated factors B, which are simply the Shavitt segment products evaluated with the upper and lower halves of the graph, respectively. At first, this may appear to require twice as many operations. However, remember that this calculation is just one in a series, with both r and s varying. The intermediate result d obtained in an outer loop overs can be reused many times with varying r. The conjugate cases, Le., r < s, are obtained by simply interchanging the ket and bra side. Any formulas involving weight operators r = s are simple to treat specially, using bit-packed occupation numbers for the upper and lower walks. A number of strategems were used to improve the efficiency of the method sketched above. The r,s loop of eq 2 was enclosed by an outermost loop over the midvertices and symmetry labels, which define the specific submatrices containing intermediate results. This reduces the memory requirement since only the vectors C and u are stored in full. Bit flags were introduced to prevent redundant operations on rows or columns known to contain only zeroes. Occasional transpositions were made to replace row by column operations (stride 1 is faster on most machines). Furthermore, many of the DAXPY calls were intercepted by an “operation collecting” subroutine, which performed a large number of DAXPY operations all at once, when its buffer needed to be emptied. This last point enables us to perform many DAXPY operations into a vector register, before storing the register to memory, whenever possible. The final u subroutine, with its dedicated routines, is then about 1900 lines of F O R T U N - 7 7 code. The number of floating-point operations is roughly the same as with a good determinant code, so the effect of the large number of coupling coefficients is compensated for by the much lower number of CI coefficients. The u vector generation reaches, on our IBM3090/VF170 machine, an effective speed (floating-point operations divided with total CPU time) of about half of the speed of the DAXPY operations alone. This indicates a somewhat large overhead of integer operations with logic, which may be trimmed down. It also means that the speed will continue to increase with increasing CI expansion length. The largest CI calculation so far, with this program, used 7 16418 configurations, corresponding to about 3.2 million Slater determinants. CI Diagonalization Even if the u routine is the time-consuming part of the CI section, the other routines also deserve some attention, since it is possible to cut down on the number of u vectors that must be calculated in order to bring a calculation to convergence. First of all, the CI eigenvalue problem is solved by a generalization of the Davidson method.14 As in our previous programs, several CI vectors can be handled simultaneously. This modification, described by Liu,ls gives faster and safer convergence to excited states, in particular when there is a near-degeneracy problem. A new modification is the use of an explicit Hamiltonian matrix evaluated within a subspace, typically comprising the most important 50 or 100 CSF’s. As long as the orbitals are far from convergence, the eigenvectors of this Hamiltonian are used to start the Davidson procedure; later on, the converged CI vectors for the previous set of orbitals are used as start vectors. Furthermore, the spectral resolution of the explicit Hamiltonian matrix is used to obtain improved update vectors in the Davidson procedure. A similar scheme was used by Werner and Knowles in their MCSCF program.I6 When a Davidson-like scheme is used to find a particular eigenvector of a Hamiltonian matrix, the residue vectors (14) Davidson, E. R. J . Compur. Phys. 1975, 17, 87. (15) Liu, 9. In Numerical Algorithms in Chemistry: Algebraic Methods; Moler, C., Shavitt. I.. Eds.; LBL-8158, Lawrence Berkeley Laboratory: Berkeley, CA, 1978. (16) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1985,82, 5053.
Malmqvist et al.
5480 The Journal of Physical Chemistry, Vol. 94, No. 14, 1990
ri = CHijCj- EC, = ui - EC, I
(5)
= hpq + CDrsIWrs) - ‘/z(prIqs)l
(13)
rq
are used to obtain the update vectors by a so-called preconditioning formula, such as
This formula applies directly to the standard Davidson procedure, and E is then the current estimated eigenvalue. Various modifications are most easily understood as the solution of a perturbation problem (Ho- E l ) q = r (7) where & should approximate the full H matrix. Thus, in the standard Davidson, H is approximated simply by its diagonal part D. In the new scheme, it is replaced by
Ho = PHP + QDQ
(8)
where P projects out the selected subspace, and Q = 1 - P is its complement. The improved start vectors and the improved preconditioning are both important and cut down significantly on the number of CI iterations needed. When the eigenvector of the preceding RAS iteration is used as the trial vector in the diagonalization, Hois modified such that it also has this trial vector as an eigenvector.” The Super-Ci Method An approximate form of the super-CI method’* was used in the first implementation of the CASSCF method.’.’ It has since then been used in a large number of application^.^ In most cases good convergence behavior is obtained with this method. Problems have occurred in specific cases, especially in calculations on closely spaced excited states, where the solution has been to optimize a set of average orbitals for the nearly degenerate states. In the super-CI method one finds the orbital rotation parameters from the solution to the super-CI secular problem, defined by a wave function comprising the MCSCF reference state 10)and the singly excited states (the Brillouin states) Ipr) = NpAkpr - k r p ) l O )
(9)
where Npris a normalization constant. These states are not orthogonal. The overlap integrals have the following form (qslpr) = 4(Pprsq - Pprqs) + W
fpq
p q
+ s#rs
- 6spDqr - a q 3 s p (10)
where D and P are the first- and second-order reduced density ma trices. The super-CI wave function is given as ISCI) = 10) + E t p @ )
(11)
P’r
In the original CASSCF formulation we obtained the new orbitals as the natural orbitals for this wave function. This procedure cannot be used for the more general case of a RASSCF wave function, since the RAS space is not invariant to a transformation into natural orbitals. Instead new orbitals are obtained by applying the unitary transformation exp(t), where t is an anti-Hermitian matrix obtained from the expansion coefficients in ( 1 1). The super-CI secular problem is not solved exactly. Instead an effective one-electron Hamiltonian is used in the calculation of the interaction elements between the Brillouin states
where the Fock matrix f is defined as (17) The authors are grateful to Jeppe Olsen for suggesting this improvement . (18) Grein. F.; Chang, T. C. Chcm. Phys. Lcrr. 1971, 12, 44. Grein, F.; Banerjee, A. In?.J . Quonrum Chem. Symp. 1975, 9, 147.
Here h and (pqlrs) are the one- and two-electron integrals in the mo&ular orbital basis. It is easy to show that for the inactive orbitals i, fn corresponds to (spin averaged) orbital energies in the Koopmans sense. In the same way f, corresponds to orbital energies for virtual orbitals a. The operator (12) thus resembles the zeroth-order Hamiltonian used in Mdler-Plesset perturbation theory. The matrix element of (12) between two Brillouin states is given in terms of a few simple matrices (qsIfiei,rr
- E O I P ~ ) = Nqs-NpA6pqGrs - 6rsHpq + DrJpq + DP$rsI
where Grs
= C(2PaBrs - Da&’rS)fap a$
H, = GP, + 4 Ipq
4 + I4P
= CDpaL,
..
(14)
where Eois the RASSCF energy and P the second-order reduced density matrix, Equation 14 has been obtained under the assumption that ErpIO) = 0. This is strictly true only for CAS type wave functions, where the only pr rotations used are r inactive and p active or external, or r active and p external. For RASSCF wave functions we also have to introduce rotations between the three different RAS spaces (but not rotations within a given subspace). The assumption made above is then approximately fulfilled if the RASl space contains only orbitals with occupation numbers c l p e to two, and RAS3 orbitals with small occupation numbers. Erpis a deexcitation operator and ErplO)will therefore contain only configurations within the RAS space. It can be excluded from (9) without loss of generality. An interesting alternative to (9) is then to define Ipr) as
kr) = ( I - P ) ~ , , J o )
(15)
where B is a projector onto the RAS CI space. This approach has, however, not been implemented in the present version of the RASSCF program. The construction of the u vector in the Davidson procedure used for solving the super-CI secular problem is performed as a series of matrix operations using the one-particle matrices defined in (14), together with the energy gradient vector, which is of course computed exactly. Both the construction of these matrices and the subsequent solution to the secular problem are very fast processes, which can also be effectively vectorized. Note also that the whole formalism is set up to need only first order transformed integrals, since it is based entirely on the Fock matrix (1 3) and the MCSCF Fock matrix used to compute the gradient vector. As a result the computational effort of the super-CI section is small, and there is no problem to use large basis sets and many inactive electrons. Actually the bottleneck in these calculations is the storage of the A 0 two-electron integrals, which limits the present implementation to about 250 A 0 basis functions. A Quasi-Newton Update Method The convergence pattern of the super-CI method is in most CASSCF calculations satisfactory. However, when the extension to RASSCF is made, one can expect more problems with convergence, since the internal rotations between the different RAS subspaces are more strongly coupled to the CI rotations. As always, when active-active rotations must be treated in an MCSCF optimization, the variation problem will contain near linear dependencies. Such near (rather than exact) linear dependencies aggravate the CI coupling problem and can also enhance the effects of nonlinearity. As a result, poor convergence is obtained. Our present implementation of the super-CI method does not explicitly include the CI coupling, and it is instead essential to look for methods that can be used to restore the coupling. The super-CI method can be regarded as an approximate form of the full second-order Newton-Raphson approach. It is therefore
Restricted Active Space SCF Method
The Journal of Physical Chemistry, Vol. 94, No. 14, 1990 5481
TABLE I: A CASSCF Chleulltion 011 the C d State of Acetykae, witb Four Electroas Distributed in Four Active Orbitals. ComprrisOa of Convergence witb and witbout a Quai-Newtoa Update (E.erpks in Atomic Units) without QUNE' with QUNEb no. of CI energy max energy max CPU
iteration 1 2 3 4 5 6 7 8 9
IO I1 12 13 14 15
iterations 2 3 3 2 3 3 4 2 2 2 2 2 2 2 2 2
lowering
derivative
-8.9E-03 -5.4E-02 -3.7E-03 -9.3 E-04 -2.4E-04 -6.5E-05 -1.8E-05 -5.28-06 -1.6E-06 -4.9E-07 -1.6E-07 -5.5E-08 -1.9E-08 -7.0E-09 -2.6E-09
6.8E-04d 1.1E-02 2.3E-02 7.5E-03 4.78-03 2.6E-03 1.5E-03 8.48-04 4.7E-04 2.7E-04 1.6E-04 1.1E-04 6.78-05 4.3E-05 2.88-05 1.8E-05
lowering
derivative
-8.98-03 -5.4E-02 -3.7E-03 -1.1E-03 -1.8E-04 -2.1E-05 -2.6E-06 -5.8E-07 -5.3E-08 -9.2E-09 -1.4E-09 -2.1 E-1 0 -1.9E-11 -2.9E-12 -1.6E-12
6.8E-04 1.1E-02 2.3E-02 7.5E-03 5.7E-03 2.3E-03 6.0E-04 2.9E-04 1.2E-04 3.3E-05 1.3E-05 6.3E-06 1.8E-06 7.8E-07 3.1E-07 3.1E-08
methodC
time, s
sx sx sx
2 2 2 2 2 2 2
LS QN QN QN QN QN QN QN QN QN QN QN
2 2 2 2 2 2 2 2
final 'Calculation without a quasi-Newton (QUNE) update procedure. bCalculation with the (QUNE) update procedure. 'SX, step taken without using Q U N E LS, a line search is made using the direction given by the last iteration; QN, a QUNE step is taken. d6.8E-04 6.8 X ol,-' natural to look for update methods that can be used to improve the approximate Hessian. As a result, the super-CI method has been augmented with a convergence accelerator, similar to that used in quasi-Newton update methods. Provided that the CI calculation is sufficiently converged for any given orbital set, the Brillouin-Levy-Berthier (BLB) vector, which contains essentially the derivatives of the energy with respect to the orbital rotation parameters, will be a function of the orbital parameters. If this function could be inverted, it would provide a correctionformula, which for any given BLB vector gives the exact orbital rotation needed. When the C I calculation is performed using the new orbitals, the exact MCSCF solution would be obtained. This is of course the foundation of all iterative methods, which may be regarded as approximating the correction formula in one way or another. In particular, in the Newton-Raphson method, one uses the Hessian matrix, containing second derivatives with respect to both rotation parameters and to C I coefficients, to obtain a correction formula that gives quadratic asymptotic convergence. In a quasi-Newton method, an approximate Hessian (or inverse Hessian) is improved in each iteration, based on the actual gradients obtained. This is done by storing a growing set of vector pairs, which provide improvements to the (inverse) Hessian." A similar scheme has been implemented in the present program, where the proposed rotation parameter matrix returned from the super43 routine is regarded as an approximate correction function of the BLB matrix. This function is successively improved by collecting a set of vector pairs, just as in a quasi-Newton method. The only difference is that this correction function is nonlinear, rather than being an approximate Hessian matrix. A very appealing feature of this method is that, although it is based on orbital rotation parameters and derivatives only, the update implicitly takes account of the coupling between rotation parameters and CI coefficients. This feature has proved to be essential when the coupling is large, as it is in many near-degenerate cases, for many excited states, and for negative ions. Table I illustrates the convergence improvement that can be obtained by using these methods. The step taken, i.e., the change in orbital parameters, in each iteration is indicated by the letters SX,QN, or LS. SX is the step proposed by the super-CI method. If the quasi-Newton step is accepted, it is denoted QN. When we are not in a quadratic minimum region, or if it is clear that the update correction is not very good, a line search is performed in the same direction as the last step, and the new step is denoted LS.The table also includes the energy lowering, and maximum gradient element, that is obtained at that iteration step if the QUNE procedure is omitted. These two columns thus refer to a completely separate calculation and have been included for comparison. A number of test calculations has been performed using the update procedure. In almost all cases the convergence is much improved, compared to
TABLE II: Results from a CASSCF Calculation on Phenol witb Eight Electrons Distributed in Seven Active Orbitals (Energies in atomic units)' CASSCF no. of CI energy max CPU
iteration
iterations
lowering
derivative method time, s
1 2 3 4 5 6 7 8 9 10 11 final
3 3 4 4 4 4 4 5 3 3 3 2
-4.3E-03 -1 SE-02 -1.3E-02 -6.7E-03 -1.3E-03 -3.7E-04 -5.2E-05 -7.8E-06 -1.2E-06 -9.2E-08 -1.3E-08
1.2E-02 7.4E-03 1.7E-02 1.7E-02 7.68-03 4.1 E-03 1.3E-03 5.4E-04 2.6E-04 8.3E-05 2.5E-05 8.2E-06
a
sx sx sx LS
QN QN QN QN QN QN QN
34 35 40 38 42 46 42 44 50 48 44
For definitions see Table I.
the super-CI method. Theoretically, the convergence is superlinear; i.e., the asymptotic convergence is superior to any first-order method, but not to a second order method, such as NewtonRaphson.
Some Test Applications The RASSCF program is now at production stage and is being employed in a number of application problems. For CASSCF calculations the new program is about 3 times as fast as the old program, which was based on the factorized two-electron coupling coefficient method.9 The upper limit for CASSCF calculations is about lo6 CSFs, depending somewhat on what computer is being used. For RASSCF calculations the limit is somewhat lower due to the larger number of M O s involved. The program stores the second-order density matrix in core, which sets the upper limit for the active M O s to about 50 for the present implementation. This number and also the CI limit can certainly be increased on systems with larger central storage (the present version uses 2 megawords (REAL*g)). Below we give some test examples from calculations with the new program. They should be considered as examples, where the aim is merely to show the program at work and not to give a full description of the problem under study. As a first example, consider a straightforward CASSCF calculation on the acetylene ground state. The orbital basis comprised 136 functions, and 4 active electrons (the *-electrons) were distributed in 4 active orbitals, giving 8 CSFs. This is a very small calculation, which could just as well have been performed by using the old program. However, a very strict convergence criterion had to be used, since the calculation was one in a series aimed at determining the quadrupole moment derivative by a difference
5482 The Journal of Physical Chemistry, Vol. 94, No. 14. 1990 TABLE III: CASSCF a d RASSCF CakuIatiOg am the V State of Ethene, Illwtrating tbe Spatial Extent of the r* Orbital as a Function of the Number of Correlating Orbitals calculation’ no. of CSF‘s eneriv.. au (zZP sec/itc 22.8 1 5d 28 848 -77.899 306 CAS (108.2~) 23.4 60 92406 -17.923856 CAS 110u.3~) 23.0 227 268800 -77.934 101 CAS iioa&j CAS ( 1 0 u , 5 ~ ) 716523 -77.943303 21.9 878 RAS (5,2,22) 4 704 -78.048 389 20.8 68
-.
’All calculations correlate the 12 valence electrons. The active orbital space is given for each CAS calculation. The RAS calculation has 5 u orbitals in RASI, 2 T orbitals in RAS2, and 16 u and 6 T orbitals in the RAS3 space. The basis set is of A N 0 type: C(14~,9p,4d,3f/6s,Sp,3d,2fj;H(7s,3p/3s,2p). bThe expectation value (E$) is a convenient measure of the amount of dynamical correlation included. The true value is believed to be- about I5 au. CSecondsper CI iteration on an IBM 3090-1708 computer. dThe same calculation with the old CASSCF program takes 40 s per CI iteration.
appro~imation.’~As Table I shows such highly converged solutions are easily obtained with the new program. The next example has a much larger number of inactive orbitals. Table I1 shows the progress of a calculation on the ground state of phenol, using 138 A N 0 (atomic natural orbital) basis functions. This is again a CASSCF calculation, with 8 active electrons distributed in 7 A orbitals, yielding 490 CSFs. The calculation time is in this case completely dominated by the orbital optimization. Even if the CI times in these two examples are irrelevant, it should be noted that only a few CI iterations are needed in each macro iteration. The number of CI iterations grows in the beginning due to an increasing demand on accuracy. Toward convergence, the number decreases again, since each CI calculation starts with the converged CI vector of the previous iteration. The time used in a super-CI macro iteration is then comparable to a micro iteration in the second-order Newton-Raphson method. The example in Table 111 illustrates the use of the RASSCF method in a study of the V state of the ethene molecule. This is a singlet PA* excited state and is dominated by two resonance VB structures with one positive and one negative methylene ion. The V state of ethene has been much studied as a prototype for excited states of high ionicity. The main difficulty here is to include enough dynamic correlation in order to obtain a realistic description of the shape of, in particular, the A* orbitals6 As is illustrated in Table 111, it is very costly to include a larger portion of the dynamic correlation by simply extending the active space in a conventional CASSCF calculation. By contrast, the RASSCF requires only a modest number of configurations. The RASSCF calculation was performed with at most two holes in RASl and at most two electrons in RAS3. The RASl space comprised all occupied valence u orbitals, RASZ the two A orbitals, and RAS3 a number of u and ?r orbitals. A RASSCF calculation of comparable accuracy (measured by the size of the ?r* orbital) to CASSCF requires less than 1% of the number of CSF‘s used in the CASSCF calculation. The time required does, however, not decrease as much, since the RASSCF calculation uses a much larger number of correlated orbitals. As a rule of thumb, the number of operations grows proportionately to the number of CSF’s, the square of the number of active electrons, and the square of the number of active orbitals. The final example, in Table IV shows an attempt to investigate whether the extra electron in the methylene anion, in a planar configuration, is bound. A large A N 0 basis set was used and RASSCF calculations were performed with an increasing number of orbitals in RAS3. The RASl space contained the CH bonding orbitals, while the RASZ space contained two A orbitals. Three different types of wave functions were used: SD with up to two holes in RASl and up to two electrons in RAS3; SDTQ with up to four holes and electrons, respectively; and SDTQ‘ which includes (19) Lindh, R.;Liu, B. Manuscript in preparation.
Malmqvist et al. TABLE IV: RASSCF C.krl.tioaS 011 the Methylene Anion, CH3; in a Planar Configuration calcuhb CIc energy no. of CSF’s ( z ~ ) ~ 1 SD -39.623 112 127 144 SDTQ -39.624633 437 146 2 SD 363 -39.654 392 SDTQ 57 978 -39.654 498 -39.657615 56 4210 SDTQ 49 3 SD -39.689 543 801 -39.962 469 48 478 1 SDTQ 1361 -39.706 120 43 4 SD 1498 SDTQ -39.707 008 38 5 SD 1695 -39.712 198 38 883 1 SDTQ -39.713 337 33 2473 6 SD -39.7 18 962 35 6’ SD 2473 -39.720 51 2 83 ‘Basis set: ANO, (C/14s,l lp,4d,3f/6~,7~,3d,2f;(H/7~,3p/3s,Zp); planar geometry; R(CH) = 1.08 A; D,, symmetry. bCarbon 1s inactive. The RASl space is 2al’, le’, RAS2: I a r . 2 a r . The RAS3 space is calculation 1: 3al’, 2e‘; 2: 1 + 3e‘,3ar; 3: 2 + 4a1’,4a/,le”; 4: 3 + 5aI’,4e’; 5 : 4 + 2e”; 6: 5+ 6a1’,5e’; 6‘: as 6 but with more diffuse A O s on carbon: 2p (0.0112, 0.0048); 3d (0.03). ‘SD: one or two holes in RASl and one or two electrons in RAS3. SDTQ: up to four holes in RASl and up to four electrons in RAS3. SDTQ: one or two holes in RASl and up to four electrons in RAS3. dThe expectation value (z2}, where z is the axis perpendicular to the molecular plane.
only quadruple excitations involving the two A electrons (up to two holes in RASl but up to four electrons in RAS3). These calculations are presented here to illustrate the flexibility in the choice of the RAS wave function. The outcome of the calculation is negative: in no case was the extra electron found to be bound, as is shown by the expectation value ( zz) given in the last column in Table IV. Adding more diffuse basis functions, as in calculation 6’, increases the value of (z2) even further. By contrast the value decreases to 33.1 au, when the molecule is bent such that the angle between a CH bond an the symmetry axis is 7 5 O .
Conclusions We have presented a new MCSCF program, which is based on RAS type wave functions. The super43 method has been used for the orbital optimization in conjunction with a quasi-Newton update procedure introduced to accelerate convergence. A split graph unitary group approach has been implemented to compute a small set of coupling coefficients, which is used to calculate the CI Hamiltonian matrix elements. The new approach allows the treatment of large CI expansions, in a spin-adapted CSF basis, without having to store a large number of coupling coefficients on disk. The limit in the present implementation is around lo6 CSFs for CASSCF wave functions. For restricted RAS spaces it is somewhat smaller and here depends on the size of the three RAS spaces. A number of examples have been given which illustrate the types of problem that can be treated with the new program. The RAS concept allows the optimization of the molecular orbitals for a wave function which includes dynamical correlation, which in many cases is necessary in order to obtain a correct description of the electronic structure. The usefulness of the RAS concept has recently been illustrated in calculations of energies and dynamic polarizabilities of CH+ and N2.20s21 Acknowledgment. We are grateful to Jeppe Olsen for many fruitful discussions. We also thank Markus Fuelscher for help with the calculations on ethene and phenol. The research reported in this paper has been supported by a grant from the Swedish Natural Science Research Council (NFR) and by IBM Sweden under a joint study contract. (20) Olsen, J.; Sinchez de Meras, A. M.;Jensen, H. J. Aa.; Jergenscn, P.Chem. Phys. Lett. 1989, 154, 380. (21) Jenscn, H. J. Aa.; Jergensen, P.;Hclgaker, T.; Olsen, J. Chem. Phys. Lett. 1989, 162, 355.
