- Email: [email protected]
First-order correlation orbitals for the MCSCF zeroth-order wave function
Volume 183, number 6
CHEMICAL PHYSICS LETTERS
13 September 1991
First-order correlation orbitals for the MCSCF zeroth-order wave function Andrzej LeS ’ and Ludwik Adamowicz * Deparfment of Chemistry, The University ofArizona, Tucson, AZ 85721, USA
Received 26 April 1991; in final form 5 June 199 1
The originally proposed first-order correlation orbital method (L. Adamowicz and R.J. Bartlett, J. Chem. Phys. 86 ( 1987) 63 14) has been extended in order to include the MCSCF zeroth-order wave function. The present method is based on the perturbation theory with the zeroth-order Hamiltonian suggested by Andersson et al. (J. Phys. Chem. 94 ( 1990) 5483). In the present work we restrict our consideration to the simplest case, i.e. the two-electron/two-orbital MCSCF reference wave function. The procedure is tested on the model system - lithium dimer -at a variety of interatomic distances. The results show that the calculated second-order correlation energy as well as the value of the second-order Hylleraas functional in the reduced space of the first-order correlation orbitals behave correctly upon dissociation of the system.
1. Introduction There are numerous chemical systems whose quantum mechanical description demands at least two Slater determinants in the molecular wave function. One can mention the following examples: ( 1) closed-shell system that dissociates into open-shell fragments, e.g. L&+2 Li, and (2) an open-shell (excited) singlet state, e.g. (2s) ’ (2~) I. These cases can be rather adequately studied by the multiconfigurational Hartree-Fock method (MCSCF) [ 11, In principle, this method may lead to the exact molecular wave function. In practice, however, due to the obvious numerical limitations, the MCSCF wave function involves rather few configurations (Slater determinants). Further improvement of the MCSCF wave function can be accomplished with the use of many body perturbation theory (MBPT) [2,3] or the multireference coupled cluster (CC) method based on the MCSCF reference (unperturbed) wave function. However, in order to perform calculations with those methods at an adequate level of accuracy one usually needs to expand the perturbation cor’ Permanent address: Department of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland. ’ Bitnet address: [email protected]
rection to the reference wave function in a rather extended basis of correlation orbitals. For larger systems this can present a prohibitive task. The purpose of the present work is to present a scheme for optimization of a compact set of correlation orbitals, which can be used in higher-order MBPT or CC calculations with the MCSCF reference wave function. The present method is similar to the previously proposed method for generating the first-order correlation orbitals (FOCOs) for a singledeterminantal reference wave function [4]. It employs the second-order Hylleraas functional and an exponential form of the orbital unitary transformation. The theoretical consideration presented here is restricted to the simplest case, i.e. the two-electron/ two-orbital MCSCF reference wave function. In the first part of this work we give an account of the theory, and next we present some numerical results on a model diatomic system.
2. Methods Our approach is based on many body perturbation theory. The zeroth-order N-electron Hamiltonian is taken in the form used by Andersson et al. [2] and Wolinski et al. [ 31:
0009-2614/91/$ 03.50 0 1991 Elsevier Science Publishers B.V. All rights reserved.
483
(1)
H( 1, .. .. N)=Ho+(H-Ho), H,, = PoFPo+ PKFPK
(2)
+ Pm FPSD+ P,,.. FP-r,.. , PO= 1Y(O)) ( Y(O)l 9
(3)
Y(O) = c A#&, K
(4) Y(‘)IFI
y/c’))
H-HOI
Yco)) ,
(6)
H-HOI
Y(O)) ,
(7) (8)
f(i)=h(i)+j(i)-k(i),
(9)
CA: K
(10)
+ C &&.Jk~, KL
1 kkk+ 1 A&kW. kcK
(11)
K,L
The C& are the reference Slater determinants contributing to the reference wave function Y(O). PK de-
notes the projection operator onto the orthogonal complement to Y(O)spanned by GK. PsD is the projection operator onto the space spanned by all determinants which can be achieved by all single, and double excitations from the reference determinants, and P,,., projects onto the triply, quadruply, ... etc. excited determinants. In the definition of the Coulomb and exchange operators, j and k, the first sum goes over the spin orbitals occupied in QK. The second sum runs over indices of those Slater reference determinants d& and aL which differ by only one spin orbital, e.g. I replacing k. The Coulomb and exchange operators are defined via their matrix elements in a usual way: (jkl)p@=
1k(lNlM2M2)
(Wv)
(kk/)pq=j-
+A2 I (core)n#>
2
(14)
This wave function provides an adequate description of the dissociation of a system with a single bond into open-shell fragments (atoms). From now onwards, the indexes c, m, n, p, q, r, s, t, u, a, /I, y, 6 will denote the molecular orbitals (not spin orbitals). By core we denote the set of doubly occupied orbitals that remain unchanged in all (normalized) reference Slater determinants, I (core&q>. In the present method we assumed that the Y(I) wave function contains only the doubly excited determinants with respect to those determinants that are present in the zeroth-order wave function I Y(l)) =
vinua1 1 yrsI (core)rS) .
(15)
rs
Such a limitation is exact for the single-determinant MBPT theory but is only approximately valid for the present perturbation theory based on the Ho Hamiltonian. In the present approach we also assume that the singly excited determinants with respect to Y(O) do not contribute significantly and may be neglected in the first approximation. Following Wolinski et al. [ 31 we define the matrix elements of the one-electron f operator in the following way: f,,=b,,+
care C U(cclpq)
c
+MWW~q) (12)
>
k(l )P( 1M2M2)
] .
UQ
-(wIW)
-OwIw>
1
dr
I
1
-0WW)
+A:()
dz
dr
2
+A,A,,~((mnlpq}-l(~~lnq))
r12
=(kalk).
I Y(O)) =A, I(core)mfi)
+AKWWpcO dz
r12
=
In the present work we show explicit formulas for the simplest case of the zeroth-order MCSCF wave function, namely,
(5)
F( 1, .... N)=
k(i)=
2.1. A particular case
+Aj -!- [ I (core)mA) + 1(core)nfi)
HoI Y(o)>=(
j(i) = E 4&hli K
13 September 199 1
CHEMICAL PHYSICS LETTERS
Volume 183, number 6
(13)
+AAJj:
((nmlm)-1(wImq)).
(16)
The zeroth-order energy can now be expressed using 484
Volume 183,number 6
CHEMICALPHYSICSLETTERS
the elements of the Fock matrix as follows:
J=
E(O)=E core+g I
virtual C yIsyru( (core)r$] Ho -E. I (core)tti>
1511(
(17) +2
where
13 September 1991
virtual C yrs( (core)rf]H-HoI
rS
Y(O))
(25)
can be written in a more compact form using the following matrix notation: f=IfXY}, r={rx,}, and v= { ( (core)@] H-HO I Y(O))}, with x, y belonging to the set of all virtual molecular orbitals, as
and
J=Tr[f(yyT+yTy)]
~=(2A:+A:)f,,,,+(2A~+A:)f,,
-,J?Tr(nT)+2
Tr(vy).
(26)
The linear configuration coefficients { Y,~}determining the first-order wave function Y(I) can be obtained via the minimization of the second-order Hylleraas functional J, which constitutes an upper bound to the second-order energy correction Et2):
In practice, the coefftcients {~,s}are determined by solving a set of linear equations, or if the virtual orbitals diagonalize the f matrix, the calculation of each yrSis reduced to the product of the ( QTvQ), matrix elements times the inverse of the expression e,t es-l?, where E,and Q denote the eigenvalues and eigenvectors off, respectively.
.I={ Yc’)]HO-E,,]YY(‘))
2.2. Reduction of the size of the virtual space
+ (A I +Az)A32$
fm, *
+2( Y”‘IH-H,,I
(19)
Yco))aE(*),
(20)
Let us now select an active space of orbitals which is a subset of all virtual orbitals. The active space orbitals will be called the first-order correlation orbitals (FOCOs ). FOCOs are obtained in the form of a rectangular unitary transformation, U, of all original virtual molecular orbitals { I r) },
where ( Yc’)l Ho-E, [ Y(l)} =
virtual C Y~J,, < (core)rfl Ho -&
I (coreNO
,
(21) ( (core)rSIHO-Eo]
(core)tti>
=f,.,+fss-E
ifr=tands=u,
=f,,
ifr#tands=u,
=t
ifr=tands#t,
( P’)l
( (core)rflH-HO
IO=
virtual
c I
U,vlr> ,
r’= 1,2, .,., NFoCo , NFoCo
(27)
We use an exponential representation of the U matrix in terms of an antisymmetric matrix R
(22)
H-Ho 1Yco))
virtual = 1 y,( (core)rfl H-H,, KS
(FOCOs method)
I Y(O)) ,
(23)
I Y(O))
=A, (rmIsm)+AZ(rn~sn) +A,l((rmIsn)+(r~lsm)). Jz
The second-order Hylleraas functional, J,
(24)
U=exp(R)
, IF= -R.
(28)
In the farther considerations the form of matrix R will be simplified by assuming that the FOCO-FOCO block as well as the block corresponding to mixing the non-active correlation orbitals (the orthogonal complement to the FOCO space) among themselves are equal to zero. This simplification results from the fact that neither active-active orbital rotations nor non-active-non-active orbital rotations can change the value of the second-order Hylleraas functional determined on the active space (spanned by FOCOs) [ 41. The Yy”) wave function is now expanded onto
485
Volume 183, number 6
the set of Slater determinants using the molecular or-
bitals from the active space FOCO
VI)=
I3 September 199 I
CHEMICAL PHYSICS LETTERS
P= P( UfUT)PT )
(34)
VR=P(UVUT)PT.
(35)
C yrJ3,1(core)r’.?) KV 3. Search for the minimum of J in the FOCO space
(29) where
Formally, this expression resembles the representation of the first-order wave function in the single-reference many body perturbation theory. In the FOCOs method [ 41, the coefficients g, are determined via the R and r”matrices. The Hylleraas functional J with FOCOs now adopts the following form:
In principle, a search for the minimal value of the J Foco(R, F) can be performed using any unconstrained minimization algorithm. The selection of an efficient optimization algorithm is an essential factor determining the practical applicability of the present method, because for a larger dimension of the FOCO space, when the number of orbital rotation parameters exceeds a few hundred, the minimization of JFoco becomes a time consuming task. In the present work we minimized JFoco using four different optimization algorithms. The function to be minimized can be defined as a function of the orbital rotation parameters only, def f(R)
JFoCo=Tr[ PUfUTPT(FpT+7TY)] -,??Tr(flT)+2Tr(PUvUTPTyl),
=minJFoCo(R,9) P
(36)
+2Tr[(UvUT)(PTyP)]
since optimal configuration coefficients, y, for a given set of FOCOs can be determined by solving a simple set of linear equations as it was discussed before. In the first algorithm (A 1) we performed a search for the optimal R following the conjugate gradient directions [5] with the gradient calculated through a numerical differentiation of p( R). The same method was used in the second algorithm (A2), but with the gradient calculated analytically [ 6 1. In the third algorithm (A3) we used the Newton-Raphson method with the gradient and Hessian matrices calculated according to the method suggested by Adamowicz and Bartlett [4]. Finally, in the fourth algorithm (A4) based also on the Newton-Raphson method the gradients and Hessian of Y were calculated using a different procedure. Below we describe the algorithm A4 in more details. We intend to solve the following non-linear equation:
=Tr(Y’r’p) +Tr(pTT)
kF( R)
(31)
where the operator P projects the entire virtual space onto the active space, i.e. r’= PUr. We would like to stress that the configuration coefficients, 7, correspond to reference determinants in which two occupied orbitals are replaced by FOCOs. Using the projector P one may write the g,, coefficients as g=(PU)TjqPU)=UT(PTpP)U
(32)
and rewrite the expression for the JFoco functional as JFoco= Tr [ ( PTTTP) (UfU’) ( PTyjP) ] +Tr[ (PTFP)(UfUT)(PTjjTP)] -ETr[
(PTjP)(PT$JTP)]
WY,
-,&?Tr(fi’) +2 Tr(vRp) =JFoco( R, y”), where 486
(33)
=o, RCR,i”
(37)
which constitutes the necessary condition for a minimum of 9;. The solution of eq. (37) is found in a series of two-step optimization cycles. In the first step
Volume 183, number 6
the elements of the R matrix are frozen (R set to R,) and the optimal 7 coefficients (i.e. jJ=j$) are obtained from (36). In the second step the linear coefficients 7 are frozen and we look for an optimal R matrix that fulfills the non-linear equation (37 ) (i.e. R=R,+ ,). The solution of eq. (37) is obtained by means of the Newton-Raphson method using the analytically calculated gradient and Hessian of 9 at R=R,,
E
y”=yn:
(46) (47)
a2R”
1
n! aR,, aR ys
=Tr(pT$jj)+Tr(jj$&pT)
(38) a26
” aR,,aR,
, ifn=2
(48)
n-l
a2R”
ii 3R,, aR,a = ,?I I(n-:)(?) ~(Aw,_,Cw,_,_,tw,_,Cw,_,_,A)
a2fR =Tr
= L (ACtCA) 2!
and 1
aR,,aR,
I3 September 1991
CHEMICAL PHYSICS LETTERS
’
a2fR aR,,aR,
‘=
(39)
ifn>3,
(40)
where A=aR/aR,, and C=aR/aR,. The first and second derivatives of zR are then assembled in the following way:
(41)
4 =p
where
a2R
and z=f, v. In order to evaluate expressions (38) and (39) we calculated a set of auxiliary matrices w,, I w,ck!R”,
k=O, l,..., h4,
IIRII”
a22R
a2eR
(49)
0 _aeR
se-R
aR,paR,& = p (- aR@ 2 e-R + aR,fl 2 aR,* aP
acR
+ aR,’
aR,,
’
(50)
(42)
where
FCC
deR
-2e-"teR2s)PT, aR,,
aR,,
(43)
and Ebeing a threshold. By means of We matrices we calculated (44)
(45)
where z=f, v. Formerly [4] (the algorithm A3), the solutions of eq. (2) have been accomplished using a procedure where the f and v matrices were redefined at each step of the Newton-Raphson loop, i.e. 2 replaced by eRfi2 ewRn, 2=f, v. In this procedure the evaluation of the next-step rotation matrix R=+, =R,+X is easy, because the first and second derivatives of (redefined) 2 are calculated at X=0. For example, the gradient and Hessian of redefined f can be expressed using the following commutators: aeX(tFf
ecR”)ecX ax,,
=[A,eRnfe-Rn],
(51)
x-o 487
Volume 183, number 6
-tic,
CHEMICAL PHYSICS LETTERS
[A,eRnfemRn]]),
(52)
where (53)
In the algorithm A4 the first and second derivatives of L are calculated directly according to expressions (49) and (50). These expressions become equivalent to expressions ( 5 1) and ( 52) at R,= 0 but they differ at R,#O. This is easily seen when comparing the power series:
a eX(eR”f
eeR”)eeX
axa,
+;
x-o
=[A,fl+[A, [Rmfll
[A, FL, [R,,flll+
1..
(54)
and aeAf e-R &,
R=Rn
=[A,fl+ $4
FLfll+;
+;(:W.,
[Rn,flll+;[~n,
+f[Rn, [R,, [A,flll)+
[Rn, [A,fll
... .
iA>[Rmflll (55)
The non-equivalence of formulas ( 54 ) and ( 55) might cause some differences during the minimization of 9. However, in the numerical example presented in the next section we found that both algorithms, A3 and A4, give nearly identical numerical results.
4. Results The calculations presented here have been done on the Li2 dimer at various internuclear separations. The zeroth-order wave functions have been obtained with the GAMESS program [ 71. First we calculated the second-order correlation correction to the MCSCF total energy, i.e. Ec2), using all virtual or488
13 September 1991
bitals. The results are presented in table 1 and in fig. 1. For comparison, in fig. 1 we plotted the SCF+E(2) energy calculated with the single-determinant reference wave function. The failure of the single-determinant approach is clear. The P2) correction does not help to improve the essentially wrong behaviour of the SCF energy. On the contrary, as one may expect, the MCSCF (two-determinantal) energy behaves correctly upon dissociation of the Li2 system. The addition of the second-order perturbation correction should farther improve this behaviour. In the next step we performed the calculations with FOCO orbitals. We reduced the dimension of the virtual space ( 14 orbitals) by 50% (7 orbitals in the FOCO space). The preliminary studies on the lithium dimer at R= 5.0a0 have shown that the computed minimal value of JFoco crucially depends on the applied optimization method and on the starting point. For example, starting with R=O the algorithms A3 and A4 converged to JFoco= -0.008473 au while the algorithms Al and A2 converged to JFoco= -0.008999 au. The latter value has been obtained with al1four algorithms, when some other preoptimized R was used. For the case of Li2 considered in the present paper the A 1 algorithm performs most efftciently and we used it in order to obtain FOCO orbitals for the description of the lithium bond stretching. For larger molecular systems, however, the Al algorithm may not be the best choice. The algorithm A3 is a good candidate. At least two reasons support its possible choice. Firstly, its viability was already demonstrated for the FOCO-SCF case [ 41. Secondly, both algorithms, A3 and A4, lead to numerically equivalent results. As we may see in table 1, such a dramatical reduction does not significantly alter the calculated values of the Hylleraas functional. The JFoCo values constitute more than 95Ohof the J values obtained with the full virtual space. It is then clear that the present FOCO-MCSCF method might become an efficient way to produce the correlation orbitals for multireference MBPT or coupled cluster calculations.
5. Discussion The present version of the FOCO-MCSCF method
Volume 183, number 6
CHEMICAL PHYSICS LETTERS
13 September 1991
Table 1 The MCSCF total energy ‘), configuration coefficients, second-order Hylleraas functional b, and its approximate value with the FOCOs method ‘) calculated for the lithium dimer MCSCF
A,
A2
JW
J’2’FocO
4.0 4.5 5.0 5.5 6.0
-
6.5
-
7.0 8.0 9.0
-
IO.0
-
0.982009 0.980953 0.911612 0.972053 0.963797 0.952531 0.938 104 0.901028 0.858682 0.818873
-0.188834 -0.194244 -0.210135 -0.234762 -0.266636 -0.304440 -0.346355 -0.433760 -0.512510 -0.513915
-0.011321 -0.010298 -0.009219 -0.008125 -0.007034 -0.005963 -0.004940 -0.003167 -0.001894 -0.001096
-0.011083 -0.010066 -0.008999 -0.007909 -0.006760 -0.005747 -0.004785 -0.003061 -0.001832 -0.001063
R
14.755312 14.770299 14.777348 14.779488 14.718852 14.776846 14.77437 1 14.769762 14.766650 14.764877
a) Distances in a,, energies in hartree, configuration coefficients in au. Li basis set (321/2 1) (contraction coefficients, exponents): stype (0.0696686, 36.8382), (0.381346, 5.48172), (0.681702, 1.11327), (-0.263127, 0.540205) (1.14339, 0.102255) (1.0, 0.0285645); andp-type (0.161546,0.540205), (0.915663,0.102255), (1.0,0.0285645). ‘r Virtual space is composed of 14 orbitals: 3 oU, 3 a,, 4 x., 4 xg ‘) Active space is composed of 7 orbitals: 2 cr., 1 o*,2 x,, 2 x*. The numerical values of all quantities calculated in the present work are available from the authors upon request.
The MCSCF
-14.‘0
I
-14.80
’ 0
+ E(2)
total
energy
for the lithium
dimer.
I 4
8 Internuclear
12 distance
16
20
a0
Fig. 1. The total energy E w+‘+~) of the Li, system calculated with the SCF and MCSCF wave functions.
for calculation of the correlation orbitals seems to be as viable and effective as the recently elaborated FOCO-SCF approach [ 4 ] +For a test case of a simple diatomic system the 5OWreduction of the size of the virtual space allows for recovery of more than 951 of the second-order correlation energy calculated with the non-reduced virtual space. We observed several numerical difficulties in the
optimization process. The probable reason for such difficulties is the existence of multiple local minima that may easily hinder the position of the global minimum. In the present work we tested four different algorithms in order to choose the most suitable way of optimizing the rotation matrix R. For the molecular system considered in the present paper we found that the algorithm Al based on a search following the conjugate directions [ 5 ] (without analytical gradient or Hessian) performs most efficiently. Such an algorithm may not be practical, however, for larger molecular systems, because its present implementation requires a large number of steps. An alternative algorithm A3 that employs analytical gradient and Hessian might be a better choice. This is the fastest way to minimize JFoco but it requires that the starting point is carefully chosen. One may raise a question whether some other choices of the active space could lead to similar results as those obtained with the FOCO space. A possible alternative can be active orbitals that constitute a subset of the eigenvectors of the first-order density matrix (approximate natural orbitals). Such a possibility has already been verified in the past with regard to the FOCO-SCF approach [ 81. The results have shown that these active orbitals are inferior with respect to the FOCO orbitals.
489
Volume 183, number 6
CHEMICAL PHYSICS LETTERS
I3 September 1991
6. Conclusions
References
The FOCO-MCSCF method allows to produce a compact set of correlation orbitals which, when used to perform higher-order multireference MBPT/CC calculations, should lead to a significant saving of the computer time. With the 50% reduction of the virtual space one may reproduce over 95O16 of the second-order correlation energy that has been calculated with the non-reduced space. The formalism of the FOCO-MCSCF method has been worked out for the simplest case, i.e. the twoelectron/two-orbital MCSCF wave function. The general case of an arbitrary CASSCF reference wave function is currently under development, and will be presented along with results of a number of test cases and timing information in the near future.
[ I] B.O. Roos, The multiconfigurational (MC) SCF method, in: Methods in computational physics, eds. G.H.F. Diercksen and S. Wilson (Reidel, Dordrecht, 1983). [ 21K. Andersson, P.-A. Malmqvist, B.O. Roos, A.J. Sadlej and K. Wolinski, J. Phys. Chem. 94 ( 1990) 5483. [ 31 K. Wolinski, H.L. Sellers and P. Pulay, Chem. Phys. Letters 140 (1987) 225. [4] L. Adamowicz and R.J. Bartlett, J. Chem. Phys. 86 (1987) 6314. [ 51R. Fletcher, TRUDGE, the Fortran subroutine based on the conjugate gradients algorithm. Practical methods of optimization (Wiley, New York, 1980). [ 61DUMING, Fortran subroutine, IMSL Inc., Houston (1987). [7] M. Dupuis, D. Spangler and J.J. Wendoloski, GAMESS, program QGOI, National Resource for Computations in Chemistry (University of California, Berkeley, 1980). [ 81 L. Adamowicz, J. Phys. Chem. 93 ( 1989) 1780.
490
CHEMICAL PHYSICS LETTERS
13 September 1991
First-order correlation orbitals for the MCSCF zeroth-order wave function Andrzej LeS ’ and Ludwik Adamowicz * Deparfment of Chemistry, The University ofArizona, Tucson, AZ 85721, USA
Received 26 April 1991; in final form 5 June 199 1
The originally proposed first-order correlation orbital method (L. Adamowicz and R.J. Bartlett, J. Chem. Phys. 86 ( 1987) 63 14) has been extended in order to include the MCSCF zeroth-order wave function. The present method is based on the perturbation theory with the zeroth-order Hamiltonian suggested by Andersson et al. (J. Phys. Chem. 94 ( 1990) 5483). In the present work we restrict our consideration to the simplest case, i.e. the two-electron/two-orbital MCSCF reference wave function. The procedure is tested on the model system - lithium dimer -at a variety of interatomic distances. The results show that the calculated second-order correlation energy as well as the value of the second-order Hylleraas functional in the reduced space of the first-order correlation orbitals behave correctly upon dissociation of the system.
1. Introduction There are numerous chemical systems whose quantum mechanical description demands at least two Slater determinants in the molecular wave function. One can mention the following examples: ( 1) closed-shell system that dissociates into open-shell fragments, e.g. L&+2 Li, and (2) an open-shell (excited) singlet state, e.g. (2s) ’ (2~) I. These cases can be rather adequately studied by the multiconfigurational Hartree-Fock method (MCSCF) [ 11, In principle, this method may lead to the exact molecular wave function. In practice, however, due to the obvious numerical limitations, the MCSCF wave function involves rather few configurations (Slater determinants). Further improvement of the MCSCF wave function can be accomplished with the use of many body perturbation theory (MBPT) [2,3] or the multireference coupled cluster (CC) method based on the MCSCF reference (unperturbed) wave function. However, in order to perform calculations with those methods at an adequate level of accuracy one usually needs to expand the perturbation cor’ Permanent address: Department of Chemistry, University of Warsaw, Pasteura 1, 02-093 Warsaw, Poland. ’ Bitnet address: [email protected]
rection to the reference wave function in a rather extended basis of correlation orbitals. For larger systems this can present a prohibitive task. The purpose of the present work is to present a scheme for optimization of a compact set of correlation orbitals, which can be used in higher-order MBPT or CC calculations with the MCSCF reference wave function. The present method is similar to the previously proposed method for generating the first-order correlation orbitals (FOCOs) for a singledeterminantal reference wave function [4]. It employs the second-order Hylleraas functional and an exponential form of the orbital unitary transformation. The theoretical consideration presented here is restricted to the simplest case, i.e. the two-electron/ two-orbital MCSCF reference wave function. In the first part of this work we give an account of the theory, and next we present some numerical results on a model diatomic system.
2. Methods Our approach is based on many body perturbation theory. The zeroth-order N-electron Hamiltonian is taken in the form used by Andersson et al. [2] and Wolinski et al. [ 31:
0009-2614/91/$ 03.50 0 1991 Elsevier Science Publishers B.V. All rights reserved.
483
(1)
H( 1, .. .. N)=Ho+(H-Ho), H,, = PoFPo+ PKFPK
(2)
+ Pm FPSD+ P,,.. FP-r,.. , PO= 1Y(O)) ( Y(O)l 9
(3)
Y(O) = c A#&, K
(4) Y(‘)IFI
y/c’))
H-HOI
Yco)) ,
(6)
H-HOI
Y(O)) ,
(7) (8)
f(i)=h(i)+j(i)-k(i),
(9)
CA: K
(10)
+ C &&.Jk~, KL
1 kkk+ 1 A&kW. kcK
(11)
K,L
The C& are the reference Slater determinants contributing to the reference wave function Y(O). PK de-
notes the projection operator onto the orthogonal complement to Y(O)spanned by GK. PsD is the projection operator onto the space spanned by all determinants which can be achieved by all single, and double excitations from the reference determinants, and P,,., projects onto the triply, quadruply, ... etc. excited determinants. In the definition of the Coulomb and exchange operators, j and k, the first sum goes over the spin orbitals occupied in QK. The second sum runs over indices of those Slater reference determinants d& and aL which differ by only one spin orbital, e.g. I replacing k. The Coulomb and exchange operators are defined via their matrix elements in a usual way: (jkl)p@=
1k(lNlM2M2)
(Wv)
(kk/)pq=j-
+A2 I (core)n#>
2
(14)
This wave function provides an adequate description of the dissociation of a system with a single bond into open-shell fragments (atoms). From now onwards, the indexes c, m, n, p, q, r, s, t, u, a, /I, y, 6 will denote the molecular orbitals (not spin orbitals). By core we denote the set of doubly occupied orbitals that remain unchanged in all (normalized) reference Slater determinants, I (core&q>. In the present method we assumed that the Y(I) wave function contains only the doubly excited determinants with respect to those determinants that are present in the zeroth-order wave function I Y(l)) =
vinua1 1 yrsI (core)rS) .
(15)
rs
Such a limitation is exact for the single-determinant MBPT theory but is only approximately valid for the present perturbation theory based on the Ho Hamiltonian. In the present approach we also assume that the singly excited determinants with respect to Y(O) do not contribute significantly and may be neglected in the first approximation. Following Wolinski et al. [ 31 we define the matrix elements of the one-electron f operator in the following way: f,,=b,,+
care C U(cclpq)
c
+MWW~q) (12)
>
k(l )P( 1M2M2)
] .
UQ
-(wIW)
-OwIw>
1
dr
I
1
-0WW)
+A:(
dz
dr
2
+A,A,,~((mnlpq}-l(~~lnq))
r12
=(kalk).
I Y(O)) =A, I(core)mfi)
+AKWWpcO dz
r12
=
In the present work we show explicit formulas for the simplest case of the zeroth-order MCSCF wave function, namely,
(5)
F( 1, .... N)=
k(i)=
2.1. A particular case
+Aj -!- [ I (core)mA) + 1(core)nfi)
HoI Y(o)>=(
j(i) = E 4&hli K
13 September 199 1
CHEMICAL PHYSICS LETTERS
Volume 183, number 6
(13)
+AAJj:
((nmlm)-1(wImq)).
(16)
The zeroth-order energy can now be expressed using 484
Volume 183,number 6
CHEMICALPHYSICSLETTERS
the elements of the Fock matrix as follows:
J=
E(O)=E core+g I
virtual C yIsyru( (core)r$] Ho -E. I (core)tti>
1511(
(17) +2
where
13 September 1991
virtual C yrs( (core)rf]H-HoI
rS
Y(O))
(25)
can be written in a more compact form using the following matrix notation: f=IfXY}, r={rx,}, and v= { ( (core)@] H-HO I Y(O))}, with x, y belonging to the set of all virtual molecular orbitals, as
and
J=Tr[f(yyT+yTy)]
~=(2A:+A:)f,,,,+(2A~+A:)f,,
-,J?Tr(nT)+2
Tr(vy).
(26)
The linear configuration coefficients { Y,~}determining the first-order wave function Y(I) can be obtained via the minimization of the second-order Hylleraas functional J, which constitutes an upper bound to the second-order energy correction Et2):
In practice, the coefftcients {~,s}are determined by solving a set of linear equations, or if the virtual orbitals diagonalize the f matrix, the calculation of each yrSis reduced to the product of the ( QTvQ), matrix elements times the inverse of the expression e,t es-l?, where E,and Q denote the eigenvalues and eigenvectors off, respectively.
.I={ Yc’)]HO-E,,]YY(‘))
2.2. Reduction of the size of the virtual space
+ (A I +Az)A32$
fm, *
+2( Y”‘IH-H,,I
(19)
Yco))aE(*),
(20)
Let us now select an active space of orbitals which is a subset of all virtual orbitals. The active space orbitals will be called the first-order correlation orbitals (FOCOs ). FOCOs are obtained in the form of a rectangular unitary transformation, U, of all original virtual molecular orbitals { I r) },
where ( Yc’)l Ho-E, [ Y(l)} =
virtual C Y~J,, < (core)rfl Ho -&
I (coreNO
,
(21) ( (core)rSIHO-Eo]
(core)tti>
=f,.,+fss-E
ifr=tands=u,
=f,,
ifr#tands=u,
=t
ifr=tands#t,
( P’)l
( (core)rflH-HO
IO=
virtual
c I
U,vlr> ,
r’= 1,2, .,., NFoCo , NFoCo
(27)
We use an exponential representation of the U matrix in terms of an antisymmetric matrix R
(22)
H-Ho 1Yco))
virtual = 1 y,( (core)rfl H-H,, KS
(FOCOs method)
I Y(O)) ,
(23)
I Y(O))
=A, (rmIsm)+AZ(rn~sn) +A,l((rmIsn)+(r~lsm)). Jz
The second-order Hylleraas functional, J,
(24)
U=exp(R)
, IF= -R.
(28)
In the farther considerations the form of matrix R will be simplified by assuming that the FOCO-FOCO block as well as the block corresponding to mixing the non-active correlation orbitals (the orthogonal complement to the FOCO space) among themselves are equal to zero. This simplification results from the fact that neither active-active orbital rotations nor non-active-non-active orbital rotations can change the value of the second-order Hylleraas functional determined on the active space (spanned by FOCOs) [ 41. The Yy”) wave function is now expanded onto
485
Volume 183, number 6
the set of Slater determinants using the molecular or-
bitals from the active space FOCO
VI)=
I3 September 199 I
CHEMICAL PHYSICS LETTERS
P= P( UfUT)PT )
(34)
VR=P(UVUT)PT.
(35)
C yrJ3,1(core)r’.?) KV 3. Search for the minimum of J in the FOCO space
(29) where
Formally, this expression resembles the representation of the first-order wave function in the single-reference many body perturbation theory. In the FOCOs method [ 41, the coefficients g, are determined via the R and r”matrices. The Hylleraas functional J with FOCOs now adopts the following form:
In principle, a search for the minimal value of the J Foco(R, F) can be performed using any unconstrained minimization algorithm. The selection of an efficient optimization algorithm is an essential factor determining the practical applicability of the present method, because for a larger dimension of the FOCO space, when the number of orbital rotation parameters exceeds a few hundred, the minimization of JFoco becomes a time consuming task. In the present work we minimized JFoco using four different optimization algorithms. The function to be minimized can be defined as a function of the orbital rotation parameters only, def f(R)
JFoCo=Tr[ PUfUTPT(FpT+7TY)] -,??Tr(flT)+2Tr(PUvUTPTyl),
=minJFoCo(R,9) P
(36)
+2Tr[(UvUT)(PTyP)]
since optimal configuration coefficients, y, for a given set of FOCOs can be determined by solving a simple set of linear equations as it was discussed before. In the first algorithm (A 1) we performed a search for the optimal R following the conjugate gradient directions [5] with the gradient calculated through a numerical differentiation of p( R). The same method was used in the second algorithm (A2), but with the gradient calculated analytically [ 6 1. In the third algorithm (A3) we used the Newton-Raphson method with the gradient and Hessian matrices calculated according to the method suggested by Adamowicz and Bartlett [4]. Finally, in the fourth algorithm (A4) based also on the Newton-Raphson method the gradients and Hessian of Y were calculated using a different procedure. Below we describe the algorithm A4 in more details. We intend to solve the following non-linear equation:
=Tr(Y’r’p) +Tr(pTT)
kF( R)
(31)
where the operator P projects the entire virtual space onto the active space, i.e. r’= PUr. We would like to stress that the configuration coefficients, 7, correspond to reference determinants in which two occupied orbitals are replaced by FOCOs. Using the projector P one may write the g,, coefficients as g=(PU)TjqPU)=UT(PTpP)U
(32)
and rewrite the expression for the JFoco functional as JFoco= Tr [ ( PTTTP) (UfU’) ( PTyjP) ] +Tr[ (PTFP)(UfUT)(PTjjTP)] -ETr[
(PTjP)(PT$JTP)]
WY,
-,&?Tr(fi’) +2 Tr(vRp) =JFoco( R, y”), where 486
(33)
=o, RCR,i”
(37)
which constitutes the necessary condition for a minimum of 9;. The solution of eq. (37) is found in a series of two-step optimization cycles. In the first step
Volume 183, number 6
the elements of the R matrix are frozen (R set to R,) and the optimal 7 coefficients (i.e. jJ=j$) are obtained from (36). In the second step the linear coefficients 7 are frozen and we look for an optimal R matrix that fulfills the non-linear equation (37 ) (i.e. R=R,+ ,). The solution of eq. (37) is obtained by means of the Newton-Raphson method using the analytically calculated gradient and Hessian of 9 at R=R,,
E
y”=yn:
(46) (47)
a2R”
1
n! aR,, aR ys
=Tr(pT$jj)+Tr(jj$&pT)
(38) a26
” aR,,aR,
, ifn=2
(48)
n-l
a2R”
ii 3R,, aR,a = ,?I I(n-:)(?) ~(Aw,_,Cw,_,_,tw,_,Cw,_,_,A)
a2fR =Tr
= L (ACtCA) 2!
and 1
aR,,aR,
I3 September 1991
CHEMICAL PHYSICS LETTERS
’
a2fR aR,,aR,
‘=
(39)
ifn>3,
(40)
where A=aR/aR,, and C=aR/aR,. The first and second derivatives of zR are then assembled in the following way:
(41)
4 =p
where
a2R
and z=f, v. In order to evaluate expressions (38) and (39) we calculated a set of auxiliary matrices w,, I w,ck!R”,
k=O, l,..., h4,
IIRII”
a22R
a2eR
(49)
0 _aeR
se-R
aR,paR,& = p (- aR@ 2 e-R + aR,fl 2 aR,* aP
acR
+ aR,’
aR,,
’
(50)
(42)
where
FCC
deR
-2e-"teR2s)PT, aR,,
aR,,
(43)
and Ebeing a threshold. By means of We matrices we calculated (44)
(45)
where z=f, v. Formerly [4] (the algorithm A3), the solutions of eq. (2) have been accomplished using a procedure where the f and v matrices were redefined at each step of the Newton-Raphson loop, i.e. 2 replaced by eRfi2 ewRn, 2=f, v. In this procedure the evaluation of the next-step rotation matrix R=+, =R,+X is easy, because the first and second derivatives of (redefined) 2 are calculated at X=0. For example, the gradient and Hessian of redefined f can be expressed using the following commutators: aeX(tFf
ecR”)ecX ax,,
=[A,eRnfe-Rn],
(51)
x-o 487
Volume 183, number 6
-tic,
CHEMICAL PHYSICS LETTERS
[A,eRnfemRn]]),
(52)
where (53)
In the algorithm A4 the first and second derivatives of L are calculated directly according to expressions (49) and (50). These expressions become equivalent to expressions ( 5 1) and ( 52) at R,= 0 but they differ at R,#O. This is easily seen when comparing the power series:
a eX(eR”f
eeR”)eeX
axa,
+;
x-o
=[A,fl+[A, [Rmfll
[A, FL, [R,,flll+
1..
(54)
and aeAf e-R &,
R=Rn
=[A,fl+ $4
FLfll+;
+;(:W.,
[Rn,flll+;[~n,
+f[Rn, [R,, [A,flll)+
[Rn, [A,fll
... .
iA>[Rmflll (55)
The non-equivalence of formulas ( 54 ) and ( 55) might cause some differences during the minimization of 9. However, in the numerical example presented in the next section we found that both algorithms, A3 and A4, give nearly identical numerical results.
4. Results The calculations presented here have been done on the Li2 dimer at various internuclear separations. The zeroth-order wave functions have been obtained with the GAMESS program [ 71. First we calculated the second-order correlation correction to the MCSCF total energy, i.e. Ec2), using all virtual or488
13 September 1991
bitals. The results are presented in table 1 and in fig. 1. For comparison, in fig. 1 we plotted the SCF+E(2) energy calculated with the single-determinant reference wave function. The failure of the single-determinant approach is clear. The P2) correction does not help to improve the essentially wrong behaviour of the SCF energy. On the contrary, as one may expect, the MCSCF (two-determinantal) energy behaves correctly upon dissociation of the Li2 system. The addition of the second-order perturbation correction should farther improve this behaviour. In the next step we performed the calculations with FOCO orbitals. We reduced the dimension of the virtual space ( 14 orbitals) by 50% (7 orbitals in the FOCO space). The preliminary studies on the lithium dimer at R= 5.0a0 have shown that the computed minimal value of JFoco crucially depends on the applied optimization method and on the starting point. For example, starting with R=O the algorithms A3 and A4 converged to JFoco= -0.008473 au while the algorithms Al and A2 converged to JFoco= -0.008999 au. The latter value has been obtained with al1four algorithms, when some other preoptimized R was used. For the case of Li2 considered in the present paper the A 1 algorithm performs most efftciently and we used it in order to obtain FOCO orbitals for the description of the lithium bond stretching. For larger molecular systems, however, the Al algorithm may not be the best choice. The algorithm A3 is a good candidate. At least two reasons support its possible choice. Firstly, its viability was already demonstrated for the FOCO-SCF case [ 41. Secondly, both algorithms, A3 and A4, lead to numerically equivalent results. As we may see in table 1, such a dramatical reduction does not significantly alter the calculated values of the Hylleraas functional. The JFoCo values constitute more than 95Ohof the J values obtained with the full virtual space. It is then clear that the present FOCO-MCSCF method might become an efficient way to produce the correlation orbitals for multireference MBPT or coupled cluster calculations.
5. Discussion The present version of the FOCO-MCSCF method
Volume 183, number 6
CHEMICAL PHYSICS LETTERS
13 September 1991
Table 1 The MCSCF total energy ‘), configuration coefficients, second-order Hylleraas functional b, and its approximate value with the FOCOs method ‘) calculated for the lithium dimer MCSCF
A,
A2
JW
J’2’FocO
4.0 4.5 5.0 5.5 6.0
-
6.5
-
7.0 8.0 9.0
-
IO.0
-
0.982009 0.980953 0.911612 0.972053 0.963797 0.952531 0.938 104 0.901028 0.858682 0.818873
-0.188834 -0.194244 -0.210135 -0.234762 -0.266636 -0.304440 -0.346355 -0.433760 -0.512510 -0.513915
-0.011321 -0.010298 -0.009219 -0.008125 -0.007034 -0.005963 -0.004940 -0.003167 -0.001894 -0.001096
-0.011083 -0.010066 -0.008999 -0.007909 -0.006760 -0.005747 -0.004785 -0.003061 -0.001832 -0.001063
R
14.755312 14.770299 14.777348 14.779488 14.718852 14.776846 14.77437 1 14.769762 14.766650 14.764877
a) Distances in a,, energies in hartree, configuration coefficients in au. Li basis set (321/2 1) (contraction coefficients, exponents): stype (0.0696686, 36.8382), (0.381346, 5.48172), (0.681702, 1.11327), (-0.263127, 0.540205) (1.14339, 0.102255) (1.0, 0.0285645); andp-type (0.161546,0.540205), (0.915663,0.102255), (1.0,0.0285645). ‘r Virtual space is composed of 14 orbitals: 3 oU, 3 a,, 4 x., 4 xg ‘) Active space is composed of 7 orbitals: 2 cr., 1 o*,2 x,, 2 x*. The numerical values of all quantities calculated in the present work are available from the authors upon request.
The MCSCF
-14.‘0
I
-14.80
’ 0
+ E(2)
total
energy
for the lithium
dimer.
I 4
8 Internuclear
12 distance
16
20
a0
Fig. 1. The total energy E w+‘+~) of the Li, system calculated with the SCF and MCSCF wave functions.
for calculation of the correlation orbitals seems to be as viable and effective as the recently elaborated FOCO-SCF approach [ 4 ] +For a test case of a simple diatomic system the 5OWreduction of the size of the virtual space allows for recovery of more than 951 of the second-order correlation energy calculated with the non-reduced virtual space. We observed several numerical difficulties in the
optimization process. The probable reason for such difficulties is the existence of multiple local minima that may easily hinder the position of the global minimum. In the present work we tested four different algorithms in order to choose the most suitable way of optimizing the rotation matrix R. For the molecular system considered in the present paper we found that the algorithm Al based on a search following the conjugate directions [ 5 ] (without analytical gradient or Hessian) performs most efficiently. Such an algorithm may not be practical, however, for larger molecular systems, because its present implementation requires a large number of steps. An alternative algorithm A3 that employs analytical gradient and Hessian might be a better choice. This is the fastest way to minimize JFoco but it requires that the starting point is carefully chosen. One may raise a question whether some other choices of the active space could lead to similar results as those obtained with the FOCO space. A possible alternative can be active orbitals that constitute a subset of the eigenvectors of the first-order density matrix (approximate natural orbitals). Such a possibility has already been verified in the past with regard to the FOCO-SCF approach [ 81. The results have shown that these active orbitals are inferior with respect to the FOCO orbitals.
489
Volume 183, number 6
CHEMICAL PHYSICS LETTERS
I3 September 1991
6. Conclusions
References
The FOCO-MCSCF method allows to produce a compact set of correlation orbitals which, when used to perform higher-order multireference MBPT/CC calculations, should lead to a significant saving of the computer time. With the 50% reduction of the virtual space one may reproduce over 95O16 of the second-order correlation energy that has been calculated with the non-reduced space. The formalism of the FOCO-MCSCF method has been worked out for the simplest case, i.e. the twoelectron/two-orbital MCSCF wave function. The general case of an arbitrary CASSCF reference wave function is currently under development, and will be presented along with results of a number of test cases and timing information in the near future.
[ I] B.O. Roos, The multiconfigurational (MC) SCF method, in: Methods in computational physics, eds. G.H.F. Diercksen and S. Wilson (Reidel, Dordrecht, 1983). [ 21K. Andersson, P.-A. Malmqvist, B.O. Roos, A.J. Sadlej and K. Wolinski, J. Phys. Chem. 94 ( 1990) 5483. [ 31 K. Wolinski, H.L. Sellers and P. Pulay, Chem. Phys. Letters 140 (1987) 225. [4] L. Adamowicz and R.J. Bartlett, J. Chem. Phys. 86 (1987) 6314. [ 51R. Fletcher, TRUDGE, the Fortran subroutine based on the conjugate gradients algorithm. Practical methods of optimization (Wiley, New York, 1980). [ 61DUMING, Fortran subroutine, IMSL Inc., Houston (1987). [7] M. Dupuis, D. Spangler and J.J. Wendoloski, GAMESS, program QGOI, National Resource for Computations in Chemistry (University of California, Berkeley, 1980). [ 81 L. Adamowicz, J. Phys. Chem. 93 ( 1989) 1780.
490