A hybrid MC-SCF method: generalised valence bond (GVB) with complete active space SCF (CASSCF)

14 June 1996

ELSEVIER

CHEMICAL PHYSICS LETTERS Chemical Physics Letters 255 (1996) 320-326

A hybrid MC-SCF method: generalised valence bond (GVB) with complete active space SCF (CASSCF) Simon Clifford, Michael J. Bearpark, Michael A. Robb Department of Chemistry, King's College London Strand, London WC2R 2LS, UK

Received 14 February 1996; in final form 2 April 1996

Abstract A hybrid GVB and CASSCF method is proposed, GVBCAS, in which a full CAS CI is performed amongst the active orbitals (the CAS part), and GVB type excitations are allowed amongst a set of GVB orbitals. We demonstrate that the computational effort required for the energy and gradients is little more than that for a CASSCF calculation with the strongly occupied GVB orbitals allocated to the inactive space, and substantially less than that for a full CASSCF in the space of active + GVB orbitals. The efficacy of the method is illustrated with an example. 1. Introduction The multi-configuration SCF (MC-SCF) method is a useful computational approach for the study of chemical reactivity problems (especially in the excited state) that cannot be described with a single reference method such as SCF. MC-SCF is most commonly used in a complete active space (CAS) form where all possible configurations that can be generated from the arrangements of the active electrons in the active orbitals are used [1]. There are three computational bottlenecks in MCSCF theory: (a) the computation of the MO electron repulsion integrals (MO-ERI) needed in the Fock matrix construction through a partial four-index transformation; (b) the CI eigenvalue problem; and (c) the four-index transformation of the two-particle density matrix in the MO basis to the AO basis during the computation of the gradient. Several strategies for the Fock matrix construction using direct methods (i.e. the MO-ERI are constructed during the MC-SCF iterations) have been discussed

elsewhere [2], and one of these algorithms is suitable for large numbers of active orbitals. However, as the number of active orbitals increases, the number of configurations increases very rapidly and the solution of the CI eigenvalue problem (b) and the gradient construction (c) become the central bottlenecks. This fact restricts routine MC-SCF computations to around a few hundred thousand configurations (--~ 12 active orbitals). One may consider various ways of selecting C! configurations (the Restricted Active Space (RAS) method of Roos and co-workers is an example [1]) and so reduce the size of the eigenvalue problem. There are two strategies that might be used: (a) one can select configurations based upon perturbation theory or (b) one can use a physical model for configuration selection. However, any method based upon perturbation selection of configurations does not resolve the problem of the gradient computation, which involves the four-index transformation of the two-particle density matrix. In this Letter, we consider an example of the alternative strategy (b), in

0009-2614/96/$12.00 © 1996 Elsevier Science B.V. All rights reserved PII S0009- 261 4(96)003 65-X

S. Clifford et al./ Chemical Physics Letters 255 (1996) 320-326

which the active space is extended using generalised valence bond (GVB) [3] electron pairs. As we shall demonstrate, this approach has the advantage that there is a partial factorisation of the two-particle density matrix that removes the bottleneck associated with the transformation required in the gradient computation and facilitates the direct construction of the Fock matrix. Further, since this method is based upon a physical model, one can understand the basis of configuration selection and predict the range of applicability in advance. Perfect pairing generalised valence bond (GVBPP) [3] is a form of MC-SCF computation in which paired excitations are performed within distinct groups (pairs) of orbitals. Both the orbitals and the 2 X 2 GVB CI coefficients are optimised in the MC-SCF procedure. This scheme can describe any situation between a fully formed bond and a totally broken one; or it can be used to add dynamic correlation to an orbital pair. However, it cannot describe the 'resonance' between different pairing schemes. In this Letter we describe a hybrid approach in which a full CAS CI is still performed amongst a reduced set of active orbitals (the CAS part), and GVB type excitations are allowed amongst the remaining GVB orbitals. This scheme produces far fewer configurations than a CAS amongst all the valence orbitals (CAS and GVB) reducing the eigenvalue bottleneck. An example illustrates the motivation for this work. An 18-orbital, 18-electron, singlet state CASSCF has 449141836 configurations. While a single benchmark computation at such a level is near the limits of current computing power, routine geometry optimisation at such a level would be impossible. However, if six well defined pairs exist in the system and can be treated as GVB pairs (and the remaining six electrons in a space of six CAS orbitals), the total number of configurations becomes 11200 and routine geometry optimisations become possible.

2. Theory 2.1. GVBCAS wavefuncfion In GVBCAS we have GVB and CAS type orbitals. Each GVB type orbital is paired to another,

Core I GVBtype/

321

CAStype 4+ .

.

GVBtype Virtual .

.

ii Fig. 1. Schematic representation of the partition of the orbital space for a GVBCAS wavefunction with two GVB pairs and four active (CAS) orbitals. GVB orbital i is paired with l and orbital j is paired with k. The labels m, n, o, p refer to active CAS orbitals.

and each pair of orbitals is associated with two electrons. We shall use i, j, k, l to represent GVB orbitals, where i is paired with l and j with k; m, n, o, p stand for CAS orbitals, and r, s, t, u for any active orbital (CAS or GVB). Fig. 1 illustrates this for CAS(4,4) with two GVB pairs. GVB pair 1 comprises orbital i paired with I. We shall use the four active orbital plus two GVB pairs example to illustrate the important points of the theoretical development. There are two main ways to formulate the development: the first, and most general, we shall refer to as GVB* CAS. Here each GVB pair interacts with all the other active orbitais. One has a single CI Hamiltonian of dimension 2 P a i r s x Nconfs where Nco,fs is the number of configurations that would exist if we were only considering excitations within the active CAS orbitals. This scheme can be thought of as a truncated CASSCF. In a second method, which we shall call GVB + CAS, the GVB orbitais only 'see' the average effect of the CAS orbitals, and vice versa. Here separate CI Hamiltonians are used for the GVB and CAS parts, which have 2 pai~ and Nconfs configurations, respectively. The coupled eigenvalue problems are solved iteratively. We shall consider only the GVB* CAS method in detail since the GVB + CAS is less general and involves minor simplifications, and we shall show subsequently, that the savings in cpu time are minor.

2.2. GVBCAS eigenvalue computation In the GVB * CAS method we form configurations that have either one of the GVB orbitals in each pair doubly occupied with all possible arrangements of the remaining active electrons within the CAS orbitals. In the four active orbital plus two GVB pairs example there are 22 x 20 = 80 such configurations.

322

S. Clifford et a l . / Chemical Physics Letters 255 (1996) 3 2 0 - 3 2 6

The structure of the CI Hamiltonian corresponding to the GVB* CAS expansion has the simple general form shown in Fig. 2. Clearly the CI Hamiltonian is very sparse and has a simple block structure that can be utilised in the eigenvalue problem. Each block of the CI Hamiltonian can be labelled ( x h in Fig. 2) by the occupancy (0 or 1) of the GVB pairs. Thus the x = 01 block contains all the configurations where GVB pair 1 (i paired with l) is in the 'ground' or 0 state with configuration i21 °, and pair 2 ( j paired with k) is in the 'excited' 1 state with configuration

O0

GVBorbitals]

~i

O1

[ o !

10

'

11

: |o:

O0 i 2 j2... kOl 0 O1 i 2 jo... k210

j°k2.

The CI matrix elements H xL can be written in the usual form as

10 i° j2... k °l 2

(KIH]L)

11 i° j°... k 212

= HKL = y ' ( r [ h C [ s ) A ~ , t + -~l ~ (rslm)B~,t~, rs

rslu

(1) where h c is the core Hamiltonian and (rsltu) are the electron repulsion integrals. We shall use K' and E to index the gconfs CAS configurations. The configuration labels K and L are compound labels that involve both the GVB and CAS labels, so that K = xK' and L = AE. We use r, s, t, u for any active orbital (CAS or GVB). The A~ L a n d n mKnLo p matrix elements (m, n, o, p are CAS orbitais) in the diagonal blocks are identical to those for the same CAS calculation without the GVB contributions, that K'L' The remaining non-zero matrix eleis AK~ = Am,. ments that involve GVB type orbitals (with labels i, j, k, l) are presented in Table 1. Referring to Fig. 2, the diagonal blocks 00,00; 01,01; 10,10; and 11,11 contain the A~.z and B mKnLo p

Table 1 GVBCAS symbolic matrix elements A rKL B r Ks t uL xK' A rt` ...i 2. . . . BiiK.L . . . i2 . . . . . Birj£j . . . i 2 . . . j2 . . . . Bijg.Lj . i 2 . . . j2. . . . . . KL Biimn KL

Bi,n~n BjkK L#

.

i 2 • •. mXn y . . .

.

.

.

. . . . . .. .

.

.

.

.

Fig. 2. The structure of the GVBCAS Hamiltonian for an example with two GVB pairs and four active (CAS) orbitals.

matrix elements and differ only in the contribution KL KL KL KL KL KL f r o m A i i , n i i i i , n i i j j , n i j i j , n i i m n , and n i m i n . Note that the off-diagonal blocks contain matrix elements only along the diagonal of the block ( K ' = L'). Matrix elements between configurations that differ in the occupancy of more than one GVB pair are zero (e.g. the 01,10 block or the 00,11 block). Thus the GVBCAS CI Hamiltonian is very sparse and the extra eigenvector iteration time compared to that for the CAS part alone will be short.

i2... i 2... i 2 . . . j2 .. . i 2 . . . j2

Value 2 1 1 -- ±2

i 2 . . . mX+ lily- i . m x + lily - I

K'L ~ Amn 1

AE

.

mn... op

. i .2 . . . . m X n. y . .. . . . . i 2. . . . . . i 2 j.2 . . . k O.l 0. . . . . . . . i 2. j o . ,

k210

- ~ A m n K t Er -f1

K' =/-,', K = A K'=E,K=A K'=E, K=A K'=E,K=A K=A K=A K'=E,K~A

S. Clifford et aL / Chemical Physics Letters 255 (1996) 320-326

2.3. Density matrix factorisation and gradient computation

323

could be simplified via the factorisation of the twoparticle density matrix. Thus Fabcd = "r'Valabcd

a,b,c,d E activeorbitals,

(5)

The reduced density matrices are defined by Fabcd = ~/abVcd -- 1 ( ~/acTbd dr" Tad)/bc) KL %+ = ~_,CKCLAr+ ,

(2)

otherwise. (6)

KL

One can write Fabcd in general as: Frstu

KL = ErLC r CLarstu,

(3)

where c r is the coefficient of configuration K in the desired CI eigenvector, and the Ar~L and Brr,Z, are the same symbolic matrix elements as for the CI Hamiltonian. In Eqs. (2) and (3) the r, s, t, u include CAS and GVB orbitals, and K and L are compound configuration labels that involve both the GVB and CAS labels. From Table 1 there will clearly be some considerable simplification of this general form when the sums in Eqs. (2) and (3) are decomposed into contributions from GVB and CAS. This will simplify the construction of the Fock matrix and more importantly the computation of the gradient. We now discuss this latter point. The general expression for the MC-SCF energy gradient is [4]:

dE dq

dhab

d( ablcd)

ab

dq

abcd

d(alb) + EXab

ab

-

dq

dV,, c + - -

dq

(4)

Here a, b, c, d refer to any molecular orbital. The (d(alb))/(dq) are the MO overlap derivatives and Xab is a symmetric Lagrangian. The derivative integrals dh+b/d q and d( ablcd)/dq are contracted with the density matrices Tab and Fab~d as they are computed. However, the derivative integrals are computed in the atomic orbital basis, and so the density matrices must be transformed into this basis. Although this is trivial for the one-particle matrix, it will scale as n 5 for the two-particle matrix, where n is the number of atomic orbitals. Because the integral derivatives must usually be computed in batches and the fully transformed two-particle matrix is commonly too large to store in computer memory, this transformation must be performed many times. Schlegel and Robb [4] showed how this process

Fabc d ~ I" fae .{,. A "~bcd ~ '/~- 'b ~vald ,

(7)

FafaC I bed = ~ab]lcd -- ~ ( ~lacVbd + ~/ad~/bc),

(8)

for all a, b, c, d, and for a, b, c, d ~ Valence: A T " v a l = F ' v a l __ /-,fac

abcd

" abcd

(9)

+abcd"

Now only the ~'~ A ~r abcd vat needs to be transformed since Fafac bcd in the AO basis can be constructed from the AO one-particle density matrices. However, in general the transformation of just the valence orbitals will still be very slow for large active spaces. The GVBCAS wavefunction permits a further factorisation that we now discuss. For GVB* CAS this factorisation only holds if the CI eigenvectors are analytic, but the difference for numerically evaluated eigenfunctions is small. For GVB + CAS the factorisation is exact. The two-particle density matrix elements involving GVB orbitals from different pairs ( F , jj, F/jij, etc.), and both GVB and CAS orbitals (l~iimn, etc.) factorise as shown in Table 2. The only elements of the density matrix that do not factorise are the Fitit where i is the GVB pair of l (F,.m = 0). Thus one has factorisation within the space of the GVB and CAS orbitals that is similar to that discussed for the active and inactive orbitals. We now show how this can be used to simplify the gradient expression.

Table 2 Factorisation of GVB contributions to the two-particle density matrix F

Value

Fiiii ~iijj

"~ii I 5T.Tj~

Fiiran

I -~T.3%.

l~,.m ."

- y. ~.

S. Clifford et al. / Chemical Physics Letters 255 (1996) 320-326

324 9

The transformed two-particle density matrix can be written as a sum of a GVBCAS part DGVBCAS ",~t~,8 (analogous to Eq. (7)) and a core part P , ~ . The GVBCAS part is given as CAS

5 paCASGVB Fig. 3. A canonical VB structure for the S z state of azulene.

We use Greek subscripts for atomic orbitals for clarity. The fully transformed one-particle density matrix is defined as

876

~

2Ya[$ YV6

+ I_GVB_CAS

7P~

6r'CAS

ab

_GVB_GVB

-- PaT P06

_GVB_GVB

-- Pa~

[JI3V

I_CAS_GVB

P~,a + 7P.¢ Pva

GVB CAS

(10)

CAS CVB

CAS CVB

GVB _CAS

where the summation includes inactive, CAS and GVB orbitals and c~ is the coefficient of AO a in MO a. We now define partially transformed CAS and GVB density matrices in the AO basis: GVB pGVB

fly

CrnCn CoCplranop

.+_ I ,GVB ~GVB

MO

= E

a

E mnop

E

pCAS E c~cffy,,,. a,8 -~"

P~

GVB

"]- E Ct~CifJCTCf(I~iiii--'2 TiiTii "~ ")liiTii ) i +

C~Ci[3"Yii'

+

( 1 l) a b c d 1 + cici c, ct ( - ~TiiYu)-

i CAS mn

- p~

(12)

(13)

The index l in the last summation in Eq. (13) is just the label of the GVB orbital that is paired to i. Thus

Fig. 4. Electron density plots for the optimized four GVB orbitals (left) and six active orbitals (right) for the GVBCAS wavefunction for t h e S t state of azulene.

S. Clifford et aL / Chemical Physics Letters 255 (1996) 320-326

the first summation is a four-index transformation over the CAS orbitals and the second is a one-index transformation over the GVB orbitals arising from the fact that the two-particle density matrix does not factorise for contributions involving orbitals that are paired. The core contribution P~v8 core to the two-particle matrix is : core __

core

core

core

active

P~'t3~8 - Pa~ Pv8 + Pa~ Pv~ /

¼t Po,

core

core

core

core

core core

active core

+ P~

Pv~

active

active --

active active

+P~,8 Phv +p~8 Pt3r -1-p~

core

core

P~r ) ' (16)

where clearly cole

pcore ~--- E Cotab CB ~lab

(17)

ab

and active _ _GVB - -

Pa~

-t'a#

(18)

CAS

"l-pat3 •

3. An example We have implemented the method in a development version of the Glaussian program [5] and carried out several applications using this formalism. Our purpose here is to illustrate the type of problem where this method might be useful and to indicate the computational efficiency of the method. In other work [6] we have studied the conical intersection (accidental degeneracy) that occurs between the S t and S O state of azulene (Fig. 3). The S t state of azulene has the valence bond representation shown in Fig. 3 with two fully formed "rr bonds 8 - 9 and 5-6. In its full generality this is clearly a 10-orbital

325

10-electron CASSCF problem (CAS10). However, the fact that the fully formed ~ bonds 9 - 8 and 5 - 6 exist suggest that a 6-orbital CAS alone (CAS6) or a 6-orbital CAS with two GVB pairs should be a good representation. Thus in this section we discuss the energies and gradients computed at a single geometry with these three methods. The geometry chosen is the S J S 0 conical intersection as characterised by CAS10. In Fig. 4 we show plots of the two pairs of GVB orbitals and the six CAS orbitais. The optimized GVB orbitals associated with the fully formed 9 - 8 and 5 - 6 "rr bonds are localised as expected for GVB pairs. If such orbitals could not be found then the GVBCAS approach would not be applicable. Table 3 gives the energy of the S t states and the S t - S O energy difference. The differences in the rms gradients from CAS10 (taken as reference) give an indication of the effect of re-optimising the geometry at a lower level of theory. The average CPU times (IBM RS6000/590) for one iteration of the CI eigenvalue evaluation and for the gradient step are also presented. The CPU time taken for the CAS10 energy step cannot be compared with the CAS6 and GVBCAS energy steps because, in the CAS10 computation, the matrix elements are re-evaluated during each iteration of the eigenvalue solver. However, for the CAS6 and GVBCAS, it is possible to store the symbolic matrix elements on disk. The forces are computed with the same method for all three procedures so the timing data are comparable. It is clear that the GVBCAS reproduces the small S I - S O energy difference (a very sensitive test for the quality of the wavefunction) near the conical intersection quite well in comparison to CAS6. Further, the rms gradient is only changed slightly from

Table 3 Results and timing data a from CAS and G V B C A S e n e r g y and gradient calculations

of configurations

Type

Number

CAS6 GVBCAS GCB + CAS CAS10

175 700 175 + 4 19460

a 4-31G basis.

E n e r g y time (s)

Gradient time (s)

A E (S l - S o)

E S I ( E h)

rms gradient

82.98 141.3 145.9 -

194.0 399.9 410.4 1501.2

0.016 0.004 0.004 0.003

- 382.6805 - 382.7353 - 382.7353 -382.7865

0.0288 0.0029 0.0029 0.0

326

S. Clifford et a l . / Chemical Physics Letters 255 (1996) 320-326

CASI0 to GVBCAS. Thus geometry changes should be small in comparison to CAS6 errors. The timing data indicates that the GVBCAS computation can be carried out with only a modest increase in CPU time relative to the CAS6. Finally it can be seen that the timings and energies for GVB * CAS and GVB + CAS are almost identical. This is because the overhead in constructing the CI matrix elements and density matrices is similar. GVB + CAS will eventually become faster when there are many GVB orbitals, and the difference between 2 pairs )< Nconfs and 2 pairs -t-Nconfs becomes large.

4. Conclusions The GVBCAS method seems to be one route to very large active spaces. GVBCAS is based on a physical model: the 'resonance' effects must be described in the CAS space and the GVB pairs cannot resonate. Clearly, the method is only applicable if these constraints are valid physically. The purpose of this Letter has been to document the theoretical and technical details. In other work we have demonstrated the applicability of the method. In particular, in a recent example, we have performed geometry optimisations on a conical intersection that occurs in 18-annulene [7] using a CAS6 space with six GVB pairs. The CAS(6,6) energy computation required 658 s per iteration whilst CAS6 with six GVB pairs increased this time to 1544 s. A full CAS(18,18) geometry optimisation would be impossible. There is only one drawback with the GVBCAS algorithm proposed. The GVB orbitals will always be localized and thus the molecular symmetry cannot be fully exploited. There is an additional benefit in GVBCAS computations that we have found useful: the treatment of quasi-redundant orbital rotation parameters. In a CASSCF if the occupation of an orbital in the active space is nearly 2 or 0 then one has redundant orbital rotations between the doubly occupied orbital and the inactive space or the near empty orbital and the virtual space. These redundancies cause instabilities in convergence. Tightly bound lone pairs are an example. Putting such orbitals into the GVB space (e.g. a lone pair + virtual lone pair) eliminates such instabilities.

In this work, we have demonstrated one solution to problem of large active spaces in CASSCF and the use of RAS wavefunctions [1] yields another. However, the problem of the use of such wavefunctions in calculations of dynamic correlation effects remains. Several years ago we proposed a multi-reference MP2 method [8] that is based upon a perturb then diagonalize approach. Since this method deals with each reference state explicitly, it has been implemented with the GVBCAS algorithm used in this work.

Acknowledgements This work was supported in part by the EPSRC (UK) under grant No. GR/J25123 and by Gaussian Inc. through a studentship for S. Clifford. We are also grateful to Michael Frisch for suggesting the GVB + CAS approach to us.

References [I] B.O. Roos, Adv. Chem. Phys. 69 (1987) 399; J. Olsen, B.O. Roos, P. Jorgensen and HJ.A. Jensen, J. Chem. Phys. 89 (1988) 2185. [2] M. Frisch, I.N Ragazos, M.A. Robb and H.B. Schlegel Chem. Phys. Letters 189 (1992) 524. [3] F.W. Bobrowicz and W.A. Goddard III, in: Modern theoretical chemistry, Voi 3. Methods of electronic structure theory, ed. H. Schaefer (Plenum Press, New York, 1977). [4] H.B. Schlegel and M.A. Robb, Chem. Phys. Letters 93 (1982) 43. [5] M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B,C. Johnson, M.A. Robb, J.R. Cheeseman, T.A. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. AI-Laham, V.G. Zakrzewski, J.V. Ortiz, J.B. Foresman, J. Cioslowski, B.B. Stefanov, A. Nanayakkara, M. Challacombe, C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres, E. S. Replogle, R. Gomperts, R.L. Martin, D.J. Fox, J.S. Binkley, D.J. DeFrees, J. Baker, J.P. Stewart, M. Head-Gordon, C. Gonzalez and J. Pople, GAUSSIAN 94 (Gaussian, Inc., Pittsburgh, PA, 1995). [6] M.J. Beaxpark, F. Bernardi, S. Clifford, M. Olivucci, M.A. Robb, B.R. Smith and T. Vreven, J. Am. Chem. Soc. 118 (1996) 169. [7] M.J. Bearpark, F. Bemardi, S. Clifford, M. Olivucci, M.A. Robb and T. Vreven, Mol. Phys. accepted. [8] J.J. McDouall, K. Peasley and M.A. Robb, Chem. Phys. Letters 148 (1988) 183.