- Email: [email protected]
A simple MC SCF perturbation theory: Orthogonal valence bond Møller-Plesset 2 (OVB MP2)
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
A SIMPLE MC SCF PERTURBATION THEORY: ORTHOGONAL VALENCE BOND MBLLER-PLESSET
8 July 1988
2 (OVB MP2)
Joseph J.W. McDOUALL, Kathryn PEASLEY and Michael A. ROBB Department of Chemistry, King’s College London, Strand, London WCZR ZLS, UK
Received 7 April 1988
A multi-reference Meller-Plesset perturbation theory - orthogonal valence bond Moller-Plesset 2 (OVB MP2) - is derived using an effective Hamiltonian formalism. The method is characterized by the use of localized CAS SCF orbitals in order to achieve a quasi-degenerate (in terms of the one-electron reference Hamiltonian) reference space. Thus the reference CI expansion corresponds to orthogonalvalence bond. The effkiency of the method is illustrated for two molecules (CH2 3B, ‘A, :H20 r,, 1.5r,, 2r,) and the Be atom, where full CI results are available.
1. Introduction
The use of low-order perturbation theory in a “sum over orbitals” formulation (i.e. many-body perturbation theory, MBPT [ 1] or Meller-Plesset [ 21 MP as it is often currently refered to) as opposed to a “sum over states” formulation (see for example the CIPSI algorithm [ 31) is now widely accepted as a reliable method for the computation of correlation energy where errors of the order of a few percent are acceptable. This acceptability is due to the experience gained using widely available standard computer programs such as the GAUSSIAN [4] series. The efficiency of these MP methods arise from the fact that the time-consuming evaluation of matrix elements over arbitrary configuration state functions is completely avoided and the summations are performed directly over a small number of four-index transformed integrals. At the moment, these techniques are only available for SCF or UHF wavefunctions. However, with the availability of full CI benchmarks, the limitations of SCF MP methods have become apparent in computations where one has partly broken bonds and an MC SCF reference is essential (see for example the review paper of Handy [ 51). Clearly, a simple MP2 method is needed that can be used with MC SCF wavefunctions. General multi-reference MBPT was first formulated by Brandow [ 61 more than 20 years ago. How-
ever, progress in applications has been very slow (see for example ref. [ 71). Recently, Wolinski et al. [ 81 have developed a very promising MC SCF MP theory. In this work we shall use an approach that lies between this later type and the original quasi-degenerate MBPT formulation of Brandow [ 61. We shall show that the essentialstep in turning the quasi-degenerate MBPT of Brandow into a practical procedure involves the use of an orthogonal valence bond OVB expansion of the reference space. Recently we have shown in some detail the way in which an MC SCF wavefunction can be transformed to VB space. On the one hand we have demonstrated that the localization of the active orbitals of a CAS SCF wavefunction [ 91 yields an orthogonal valence bond wavefunction OVB and on the other we have been able to show how the transformation to a non-orthogonal (spin-coupled) VB expansion [ lo] can be accomplished.
2. Theory In general there are two variants of multi-reference perturbation theory: ( 1) a “perturb then diagonalize” approach involving the construction of an effective Hamiltonian (for a review see ref. [ 111): the quasi-degenerate MBPT of Brandow [ 61 is of this type or (2 ) a “diagonalize then perturb” approach
0 009-2614/88/$ 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division )
183
Volume148,number2,3
CHEMICAL PHYSICSLETTERS
(see for example the CIPSI method ref. [ 3 ] ): the recent algorithm of Wolinski et al. [ 8 ] or multi-reference coupled cluster (see for example ref. [ 12 ] and references therein) approaches are of this type. Of course, the “diagonalize then perturb” approach or the “perturb then diagonalize” approach will yield the same result at infinite order but the results at low order may differ significantly. We begin with a brief summary of the equations of multi-reference perturbation theory in a simple matrix formulation, We divide the full CI space into a reference space P spanned by the MC SCF (CAS SCF) configurations and a secondary space Q containing the remaining configurations. We seek the solution of the Bloch equations (see ref. [ 111 for a full discussion)
where H, means the projection of H onto the P sub space ( H, = PHP), etc. The diagonalization of Herr gives a subset of the eigenvalues of H exactly. The MC SCF CI energy is obtained from the solution of the zerotharder problem HppAA= EjA* .
(2)
The transformation U is given from the solution of H,U,+H,U,-U,H,,=O
(3a)
and H,U,
-t H, U, = U,, Herr .
( ;
-;=
> .
(4)
We then obtain a set of coupled equations for each reference function of the form
184
where the subscripts R and S run over the reference configurations and {H,}R denotes the Rth column of the sub-block H, of H. {Herr}RS is given by (from eq. (3b)) {He&={HppJRS+
({H,}R)~-G.
(6)
If we ignore the coupling in the third term in eq. (5) we have (H,}R+(Hqq-(Heff}RR1)*CR=O.
(7)
Eq. (7) has exactly the same form as the infinite-order perturbation equations for single reference perturbation theory (see for example the derivation of MBPT in ref. [ 131). The neglected coupling terms (I:,+,{ HeffJR.$S) give rise to the so-called “folded diagrams” in Brandow’s theory [ 6 ] and do not enter until third order. We can now proceed to some simple working equations. First, we observe that the total energy is given (approximately) from Ra = {AA)~H,K{AA 1,
(8)
where the A, are obtained from the zero&order (MC SCF) problem eq. (2), as opposed to diagonalizing HGrr.The use of eq. (8) rather than the diagonalization of Heff has the advantage that the single excitation contribution remains nil. Thus we are using a method that is mid-way between “diagonalize then perturb” and “perturb then diagonalize”. Secondly, one expects that the off-diagonal corrections to H,* will be small so that eq. (6 ) is replaced by
(3b)
If H, is first diagonalized then we have a “diagonalize then perturb” method otherwise we have a “perturb then diagonalize” method. We can choose intermediate normalization so that U,= 1 and U then has the form UE
8July 1988
{He&S= (Hpp}RS+~~RS({Hqp)R)T.CS.
(9)
Now, eq. (8) takes the simple form (lOa)
-& =G! + 5 IA&& with A& = ({H&)=G
,
(lob)
where the CR are obtained from the solution of eq. (7). Note that eqs. (10a) and (lob) differ from a “diagonalize then perturb” type of perturbation theory (see for example ref. [ 31) only in the way in which the CR are computed.
CHEMICAL PHYSICS LETTERS
Volume 148, number 2,3
A simple MP perturbation theory now arises by ( 1) expanding the inverse occurring in the solution of eq. (7) CR=-(H*~-{H~~~}RR~)-‘IH~~]R
(11)
in a power series of the form (A-B)-‘=A-‘+A-‘BA-‘+...
,
where A is the matrix of the diagonal elements of H,- %rLr~l, (2) replacing the “sum over states” with a “sum over orbitals” and (3) replacing the diagonal elements of H,,{HeffJRR1 with orbital energy differences obtained from the diagonal elements of some reference oneelectron Hamiltonian hR. Note thatwehavethefreedom to use a difirent reference hR for each referencestateR in a “‘perturbthen diagonalize”formalism(i.e. wecan use a differentA/ B partitionfor each reference function in eq. (11)). It has been the choice of hR that has caused the difficulty with possible implementations of multireference MP theory. (See the discussion of the choice of hR in ref. [ 8 ] ). The use of a different reference hRfor each referencestateR and the use of an orthogonalvalencebond reference space is criticalto thesuccessof the present method.When the “sum over states” is replaced by a “sum over orbitals” one must distinguish three types of orbitals: inactive (or core), active (or valence) and virtual orbitals. For a given reference state R certain active orbitals are occupied while the others are unoccupied. In some other reference state S the roles of these same occupied/unoccupied orbitals will be completely reversed. Now we suppose that the occupied active orbitals in R have negative orbital energies while the unoccupied orbitals have positive energies. If the same reference h is used for reference state S then obviously the occupied active orbitals in S will have positive orbital energies while the unoccupied active orbitals have negative energies and the perturbation expansion may be hopelessly divergent. The preceding simple argument is enough to convince one that a low-order muiti-reference MP theory will only be successful if (a) the orbitals are quasi-degenerate and (b) a different hR is used for each reference state R. The use of a different hR for each reference state R is simple
8 July 1988
to implement. Requiring that the orbitals are quasidegenerate implies the use of a valence bond reference space. In previous work [9] we have shown that localization of the active orbitals of a CAS SCF wavefunction leads to an orthogonal VB expansion. Thus the CI expansion becomes a linear combination of atomic or fragment electronic structures (neutral and ionic). In this context, the obvious form of hR is the UHF form,
and
where 1.;1”and lffl include the inactive t “occupied” active in reference state R. For canonical MC SCF orbitals hg and h; will be approximately diagonal. In MP methods one uses the diagonal elements of hi and hj as orbital energies. For canonical orbitals we expect a substantial gap between the occupied and unoccupied orbital energies. However, for localized CAS SCF orbitals the orbital energies for a particular reference state R should be quasidegenerate (i.e. similar to those in the atom or fragment rather than in the molecule). The AER of cq. ( lob) become similar (identical at infinite interfragment distance) to sums of the correlation energies of neutral, anionic or cationic states of the isolated atomic or molecular fragments. However, they must of course exclude the correlation of two active “holes” with two active “particles” to avoid overcounting since this effect is included in the reference space. This idea has been followed by Ruedenberg et al. [ 141 in a semi-empirical manner for diatomic mOlcCUleSby extracting &??R from atomic computations. Of course, since hE and h$ are not diagonal in a basis of localized orbitals (except at infinite interfragment distance), there must be terms at third order that arise from this effect. We can now summarise an algorithm for multi-reference MP2 according to the scheme outlined above: ( 1) A CAS SCF computation in a basis of determinants (rather than configuration state functions). (2) Localization of CAS SCF orbitals [ 9 ] to give 185
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
an orthogonal valence bond expansion OVB. (3 ) If I {A~}R I 2 “threshold” then compute hz and h$ for reference configuration R. (4) If ({A,}, 12 “threshold” then compute AE, for reference configuration R using the usual MP2 formula excluding the “reference states” from the summation (i.e. the correlation of two active “holes” by two active “particles” is excluded since this state is in the reference space). (5) Compute the total energy accordig to eq. (10a).
3. Results and discussion As has become the accepted practise, we now discuss the results of OVB MP2 in comparison with some full CI computations [ 51. We have chosen the full CI computations of Bauschlicher and Taylor [ 15] on the singlet-triplet separation in CHr (doublezeta + polarization ) , the computations of Handy and co-workers [ 5,161 and Brown et al. [ 171 on Hz0 at (at the double-zeta level) r,, 1.5r,, and 2r, and the Be results of Heully and Daudey [ 17] since they are classic problems where SCF plus single reference MP2 fails. The CH2 problem is now well understood. Aside from basis set effects, the error in any single reference treatment (SCF, MP2 or SDCI) of the ‘A, state arises from the fact this state is a combination of two configurations - (a, ) * and a small amount of (b, ) *. On the other hand the ‘B, state has the configuration (a, ) ’(b, ) ’and is well described at the SCF level. The results are collected in table 1 along with those of Bauschlicher and Taylor [ 151. It should be noted that for the 3B, state, MC SCF is identical to restricted open shell SCF (ROSCF), and the MC SCF SDCI and OVB MP2 are equivalent to a restricted open shell. Similarity, the UHF result for the ‘A1 state is jsut ordinary closed shell. Of course the MC SCF and MC SCF SDCI results are in very good agreement with the full CI with respect to the ‘A1-3B1energy difference. While UHF MP2 is an improvement over the UHF result, it requires OVB MP2 to reduce the error to about the sort of level that one might expect for a single reference MP2 in an appropriate problem. The results on Hz0 are collected in table 2, along with the full CI results of Handy [ 161 and the MC 186
8 July 1988
Table 1 Comparison of OVB MP2 results with CI results for CH2 Method
‘A,
3B,
A()B,-‘A,) (kJ mol-‘)
UHF MC SCF UHF MP2 OVB MP2 MC SCF SDCI =’ exact (full CI) c,
-38.8863 =’ -38.9076 -39.0099 -39.0154
-38.9326 - 38.9276 b, - 39.0382 - 39.0366 b,
121.5 52.7 74.1 55.6
- 39.0222
-39.0416 b,
51.0
- 39.0272
- 39.0463
so. 1
a) Equivalent to RHF. b, The triplet 2 determinant (S,=O) MC SCF is identical to a restricted open shell SCF (ROHF). “Ref. [lS].
SCF SDCI results of Brown et al. [ 171. We have included the results at r,, 1.5r, and 2r,. Further we have added the results of a computation run at r,= 10 8, for comparison with the results at infinite r, of Brown et al. [ 171. In fact the reference space chosen in some of the computations is not balanced at infinite re so we include results at r,= 10 8, merely to indicate that our method behaves sensibly as r gets large. We have included results for two different (CAS SCF) reference spaces denoted in table 2 as “n in m” where n is the number of active electrons and m is the number of active orbitals. In each case the valence orbitals were localized using the procedure described in ref. [ 91 subject to the constraint that the 0 atom orbitals retained Czv symmetry (i.e. rather than a hybridized form). The “4 in 4” reference space included the 0 2p(b2), the 0 2p(a,) and the two H 1s orbitals with the 0 2p (b, ) lone pair and the 0 2s in the inactive space. This is the smallest space that can describe the dissociation into atoms. The corresponding MC SCF SDCI computations of Brown et al. [ 171 had no inactive space. In the “6 in 6” reference space the 0 2p( b, ) lone pair and a corresponding 0 3p (b, ) virtual orbital was added to the active space. The localized orbitals deserve some comment for the “6 in 6” reference space. In this case the two 0 p(b, ) orbitals localized into a “tight” orbital and a “diffuse” orbital where the diagonal elements of the density matrix were both approximately 1 (as opposed to the canonical delocalized situation where the occupancies were approximately 2 and 0).
Volume 148, number 2,3
CHEMICAL PHYSICS LEI-IERS
8 July 1988
Table 2 Comparison of MC SCF MP2 results with CI results for Hz0 Method
r.
I.%,
2.or,
1oA
SCF MP2 exact (full CI) 9)
-76.0098 -76.1493 -76.1579
-75.8035 - 75.9946 -76.0145
-15.5952 -75.8525 -75.9052
- 75:*648 b’
4in4”) MC SCF MC SCF SDCI d, OVB MP2
-76.0596 -76.1559 -76.1259
- 75.9244 - 76.0125 - 75.9843
-75.8273 -75.9033 - 75.8904
-15.7958 -75.8631 bi -75.8535
6 in 6 ” MC SCF MC SCF SCDI f’ OVB MP2
- 76.0970 -76.1575 -76.1442
- 75.9527 -76.0140 -76.0028
-75.8442 -75.9047 -75.8942
-75.8057 -15.8642 b’ -75.8573
a1 Refs. [5,16]. b, Ref. [ 171 result obtained at 100 A. =) Active space consisted of four electrons in the four O-H (H 2s and 0 2p) bonding orbitals. *) Ref. [ 171. In these calculations all the inactive orbitals (0 Is, 2s, 2p) were included in the reference space as well as the orbitals in c). e, To the active space in c) was added the additional 0 2p (lone pair) and a corresponding 0 3p orbital. ‘) Ref. [ 171. In these calculations all the inactive orbitals (0 ls, 2s) were included in the reference space as well as the orbitals in e).
This situation is essential in order to active a quasidegenerate reference. If the two 0 p (b , ) orbitals are not localized the 0 3p (b, ) virtual orbital has a high orbital energy and the quasi-degeneracy requirement cannot be satisfied. From table 2 it can be seen that while the OVB MP2 results for the “4 in 4” reference space are quite encouraging for r > 2r, the results for r < 2r, are rather poor. However, the addition of the two 0 p (b, ) orbitals in the “6 in 6” reference space improves the agreement of OVB MP2 with MC SCF SCDI and full CI considerably. In this case, while the absolute error at each value of r is still quite large ( x 10 rn& ), the relative error is similar at each point ( 10 m& and 2re and 14 n$, at re). Notice that the error in single reference MP2 at r, is 9 mEi,. Thus as r approaches r = co, where the OVB wavefunction is dominated by one term, we should expect an error of the same order of magnitude. We have thus accomplished our objective since the error in OVB MP2 is about 10 rnE,, at all values of r. The addition of the two 0 p(b, ) orbitals to the reference space at r, has a very large effect because at this geometry the OVB CI expansion is dominated by the 02-2H+ configuration so the additional atomic intra-atomic correlation is very important. The extension of the computations beyond this
point is not easy from a technical point of view. On the one hand, the addition of any further weakly occupied orbitals to the active space of the MC SCF is not without problems. As the occupation number of a weakly occupied orbital begins to approach 0.0, the MC SCF orbital rotation Hessian becomes almost singular. Thus the added weakly occupied orbital is not well determined and the localization step is not well determined either. Further, since such weakly occupied orbitals cannot be localized, they will have high orbital energies and the quasi-degenerate requirement becomes difficult to satisfy. Of course one could go to higher order (a) by computing the offdiagonal terms in H,, and (b) including third-order effects. It is of interest to comment on the semi-empirical scheme of Ruedenberg et al. [ 141 in relation to the present approach. In the approach of Ruedenberg et al., the AER (eq. ( lob) ) are replaced by correlation energies of neutral and ionic states of the constituent atoms. On the one hand this semi-empirical estimate does not allow for changes in the correlation energy of the atom when placed in the molecule. On the other hand, the correlation of two active “holes” by two active “particles” is overcounted in Ruedenberg’s scheme since this term is included in AERas well as in the reference space. In table 3 we have shown the 187
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
8 July I988
Table 3
Contribution of some important AER(eq. ( 1Ob) ) for Hz0 (4 in 4 CAS SCF) Configuration
&IIG%I
R r,
1.5r,
2.Or,
IOA
O('P)HH
-0.057/0.27
-O.OS5/0.38
-0.058/0.50
-0.058/1.0
0-H+H 02-H+H+
-0.064/0.27 -0.084/0.33
-0.068/0.21 -0.089/0.20
-0.084/0.12 -0.128/0.06
values of AER and 1{A,}, 1 for the most important OVB configurations for H20. As expected for an OVB
expansion, the 02-H+H+ configuration is very important at r,. Thus in order to describe the difference E(2r,) -E(r,) it is necessary to get AE(O*-2H+) correct. Thus the “6 in 6” reference space which includes the dominant 0 p (b , ) * electron correlation is essential. Note that while AE(03P 2H’S) is almost independent of r,, the value of AE(02-2H+) is strongly dependent on r. This fact seems to preclude the use of semi-empirical estimates. Finally, we give a brief discussion on the comparison of the results obtained with the method used in this work with the exhaustive study of Heully and Daudey [ 18 ] on the Be atom. Ordinary MP methods converge very slowly for Be because of the near degeneracy of the 2s and 2p orbitals. Further as demonstrated by Heully and Daudey [ 18 ] the multi-reference MP series diverges because the (2s) l (ns)’ state lies lower in energy than the (2p)* state. Heully and Daudey avoid this problem in an elegant way using an “intermediate” effective Hamiltonian. Our results are collected in table 4 along with the results of Heully and Daudey [ 181. We have used the same 6-3 1 1 + C ( 1d) basis as Heully and Daudey. The results of standard MP2-MP4 calculations indicate the expected slow convergence. The corresponding effective Hamiltonian results of Heully and Daudy [ 18 ] are also given. Their results are much improved over the usual MP2 at second order. It can be seen that the OVB MP2 result (CAS SCF “2 in 4” reference space) has a very small error indeed. The localized orbitals in this case were an exactly degenerate sp3-type combination of the CAS SCF active orbitals. Heully and Daudey [ 181 show that the multi-reference perturbation expansion is truly divergent in this example while the intermediate Hamiltonian (which gives only the Be ‘S state energy 188
/o.o IO.0
Table 4 Comparison of OVB MP2 results for Be( ‘S) with the MBPT results of Heully and Daudey [ 181 Method
E
“exact” (full CI ref. [ 181)
- 14.6445
0.0
this work SCF MP2 MP3 MP4 MC SCF OVB MP2
-
14.5720 14.6186 14.6298 14.6343 14.6156 14.6435
0.0725 0.0259 0.0147 0.0102 0.0289 0.0010
- 14.6251 - 14.6486
0.0194 -0.0041
Error
ref. 1181 Hen SC MBPT-2 Herr MC MBPT-2
exactly) is convergent. Thus, in order to extend the present computations beyond second order, one would be forced to use similar methods.
4. Conclusions A practical implementation of multi-reference MP2 depends upon the quasi-degeneracy of the active orbitals (which implies the use of an OVB expansion of the MC SCF wavefunction) and the use of a different reference one-electron Hamiltonian for each reference configuration. The results on Hz0 show the importance of the correlation energy of ionic fragments. Since we have implemented the method in the most obvious way, a detailed discussion of computational efficiency is impossible. However, at its simplest level computation times for OVB MP2 will be simply Nxsingle reference MP2 where N is the number of reference configurations. The MC SCF method is often criticized because
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
the choice of reference space (i.e. active orbitals) may seem rather arbitrary. The transformation to an OVB expansion via the localization of the active orbitals [ 91, yields a wavefunction which is easily interpreted and thus demonstratably correct for a given problem. While the orthogonality of the orbitals in an OVB expansion yields a wavefunction with large ionic contributions, we have demonstrated how these contributions can be eliminated [ lo]. However, in the present method, these ionic terms cause no difficulties. The use of a non-orthogonal valence bond reference space would in fact destroy the very simple interpretation of the correlation energy at r, in H20 (see the previous discussion of the contribution AE(O’-2H+)).
Acknowledgement The authors are ‘grateful for the support of the SERC (UK) via grant No. GR/E 4499.4. The localized orbitals were computed using the algorithm described in ref. [ 91; however, the starting orbitals were in each case generated using the Boys method and the authors are grateful to Dr. A. Bottoni for his programs. Finally, we have used the GAUSSIAN [ 41 suite of programs in all our computations.
References [ 1 ] J. Goldstone, Proc. Roy. Sot. A 239 ( 1959) 267. [ 2 ] C. Meller and MS. Plesset, Phys. Rev. 46 ( 1934) 6 18.
8 July 1988
[ 31 B. Huron, P. Rancurrel and J.P. Malrieu, J. Chem. Phys. 73 (1973) 5745. [4] J.S. Binkley, R.A. Whiteside, R. Krishnan, R. Seeger, D.J. DeFrees, H.B. Schlegel, S. Topiol, L.R. Kahn and J.A. Pople, GAUSSIAN 80, QCPE 13 (1980); J, Binkley, M, Frisch, K. Raghavachari, D.J. DeFrees, H.B. Schlegel, R. Whiteside, E. Fluder, R. Seeger and J.A. Pople, GAUSSIAN 82, Carnegie-Mellon University. [ 51 N.C. Handy, Faraday Symp. Chem. Sot. 19 (1984) 17. [6] B.H. Brandow, Rev. Mod. Phys. 39 (1967) 771. [7] R.J. Bartlet, Ann. Rev. Phys. Chem. 32 (1981) 359; I. Lindgren, Physica Scripta 32 (1985); D. Mukhergec, Intern. J. Quantum Chem. Symp. 20 (1986) 409. [ 81 K. Wolinski, H. Sellers and P. Pulay, Chem. Phys. Letters 140 (1987) 225. [9] J.J.W. McDouall and M.A. Robb, Chem. Phys. Letters 132 (1987) 319. [lo] J.J.W. McDouall and M.A. Robb, Chem. Phys. Letters 142 (1987) 131. [ 111 Ph. Durand and J.P. Malrieu, Advaa. Chem. Phys. 67 (1987). [ 121M.R. Hoffman and J. Simons, J. Chem. Phys. 88 (1988) 993. [ 131 M.A. Robb, in: Computational techniques in quantum chemistry and molecular physics, eds. G.H.F. Diercksen, B.T. Sutcliffe and A. Veillard (Reidel, Dordtecht, 1975) pp. 435 t-f. [ 141M.M. Schmidt, M.T.B. Lam, S.T. Elbert and K. Ruedenberg, Theoret. Chim. Acta 68 ( 1985) 69. [ 151 C.W. Bauschlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 5510. [ 161N.C. Handy, Chem. Phys. Letters 74 (1980) 280; R.J. Harrison and N.C. Handy, Chem. Phys. Letters 95 (1983) 386; P. Saxe, H.F. Schaefer III and NC. Handy, Chem. Phys. Letters 79 ( 1981) 202. [ 171F.B. Brown, I. Shavitt and R. Shepard, Chem. Phys. Letters 105 (1984) 363. [ 181J. Heully and J. Daudey, J. Chem. Phys. 88 (1988) 1046.
189
CHEMICAL PHYSICS LETTERS
A SIMPLE MC SCF PERTURBATION THEORY: ORTHOGONAL VALENCE BOND MBLLER-PLESSET
8 July 1988
2 (OVB MP2)
Joseph J.W. McDOUALL, Kathryn PEASLEY and Michael A. ROBB Department of Chemistry, King’s College London, Strand, London WCZR ZLS, UK
Received 7 April 1988
A multi-reference Meller-Plesset perturbation theory - orthogonal valence bond Moller-Plesset 2 (OVB MP2) - is derived using an effective Hamiltonian formalism. The method is characterized by the use of localized CAS SCF orbitals in order to achieve a quasi-degenerate (in terms of the one-electron reference Hamiltonian) reference space. Thus the reference CI expansion corresponds to orthogonalvalence bond. The effkiency of the method is illustrated for two molecules (CH2 3B, ‘A, :H20 r,, 1.5r,, 2r,) and the Be atom, where full CI results are available.
1. Introduction
The use of low-order perturbation theory in a “sum over orbitals” formulation (i.e. many-body perturbation theory, MBPT [ 1] or Meller-Plesset [ 21 MP as it is often currently refered to) as opposed to a “sum over states” formulation (see for example the CIPSI algorithm [ 31) is now widely accepted as a reliable method for the computation of correlation energy where errors of the order of a few percent are acceptable. This acceptability is due to the experience gained using widely available standard computer programs such as the GAUSSIAN [4] series. The efficiency of these MP methods arise from the fact that the time-consuming evaluation of matrix elements over arbitrary configuration state functions is completely avoided and the summations are performed directly over a small number of four-index transformed integrals. At the moment, these techniques are only available for SCF or UHF wavefunctions. However, with the availability of full CI benchmarks, the limitations of SCF MP methods have become apparent in computations where one has partly broken bonds and an MC SCF reference is essential (see for example the review paper of Handy [ 51). Clearly, a simple MP2 method is needed that can be used with MC SCF wavefunctions. General multi-reference MBPT was first formulated by Brandow [ 61 more than 20 years ago. How-
ever, progress in applications has been very slow (see for example ref. [ 71). Recently, Wolinski et al. [ 81 have developed a very promising MC SCF MP theory. In this work we shall use an approach that lies between this later type and the original quasi-degenerate MBPT formulation of Brandow [ 61. We shall show that the essentialstep in turning the quasi-degenerate MBPT of Brandow into a practical procedure involves the use of an orthogonal valence bond OVB expansion of the reference space. Recently we have shown in some detail the way in which an MC SCF wavefunction can be transformed to VB space. On the one hand we have demonstrated that the localization of the active orbitals of a CAS SCF wavefunction [ 91 yields an orthogonal valence bond wavefunction OVB and on the other we have been able to show how the transformation to a non-orthogonal (spin-coupled) VB expansion [ lo] can be accomplished.
2. Theory In general there are two variants of multi-reference perturbation theory: ( 1) a “perturb then diagonalize” approach involving the construction of an effective Hamiltonian (for a review see ref. [ 111): the quasi-degenerate MBPT of Brandow [ 61 is of this type or (2 ) a “diagonalize then perturb” approach
0 009-2614/88/$ 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division )
183
Volume148,number2,3
CHEMICAL PHYSICSLETTERS
(see for example the CIPSI method ref. [ 3 ] ): the recent algorithm of Wolinski et al. [ 8 ] or multi-reference coupled cluster (see for example ref. [ 12 ] and references therein) approaches are of this type. Of course, the “diagonalize then perturb” approach or the “perturb then diagonalize” approach will yield the same result at infinite order but the results at low order may differ significantly. We begin with a brief summary of the equations of multi-reference perturbation theory in a simple matrix formulation, We divide the full CI space into a reference space P spanned by the MC SCF (CAS SCF) configurations and a secondary space Q containing the remaining configurations. We seek the solution of the Bloch equations (see ref. [ 111 for a full discussion)
where H, means the projection of H onto the P sub space ( H, = PHP), etc. The diagonalization of Herr gives a subset of the eigenvalues of H exactly. The MC SCF CI energy is obtained from the solution of the zerotharder problem HppAA= EjA* .
(2)
The transformation U is given from the solution of H,U,+H,U,-U,H,,=O
(3a)
and H,U,
-t H, U, = U,, Herr .
( ;
-;=
> .
(4)
We then obtain a set of coupled equations for each reference function of the form
184
where the subscripts R and S run over the reference configurations and {H,}R denotes the Rth column of the sub-block H, of H. {Herr}RS is given by (from eq. (3b)) {He&={HppJRS+
({H,}R)~-G.
(6)
If we ignore the coupling in the third term in eq. (5) we have (H,}R+(Hqq-(Heff}RR1)*CR=O.
(7)
Eq. (7) has exactly the same form as the infinite-order perturbation equations for single reference perturbation theory (see for example the derivation of MBPT in ref. [ 131). The neglected coupling terms (I:,+,{ HeffJR.$S) give rise to the so-called “folded diagrams” in Brandow’s theory [ 6 ] and do not enter until third order. We can now proceed to some simple working equations. First, we observe that the total energy is given (approximately) from Ra = {AA)~H,K{AA 1,
(8)
where the A, are obtained from the zero&order (MC SCF) problem eq. (2), as opposed to diagonalizing HGrr.The use of eq. (8) rather than the diagonalization of Heff has the advantage that the single excitation contribution remains nil. Thus we are using a method that is mid-way between “diagonalize then perturb” and “perturb then diagonalize”. Secondly, one expects that the off-diagonal corrections to H,* will be small so that eq. (6 ) is replaced by
(3b)
If H, is first diagonalized then we have a “diagonalize then perturb” method otherwise we have a “perturb then diagonalize” method. We can choose intermediate normalization so that U,= 1 and U then has the form UE
8July 1988
{He&S= (Hpp}RS+~~RS({Hqp)R)T.CS.
(9)
Now, eq. (8) takes the simple form (lOa)
-& =G! + 5 IA&& with A& = ({H&)=G
,
(lob)
where the CR are obtained from the solution of eq. (7). Note that eqs. (10a) and (lob) differ from a “diagonalize then perturb” type of perturbation theory (see for example ref. [ 31) only in the way in which the CR are computed.
CHEMICAL PHYSICS LETTERS
Volume 148, number 2,3
A simple MP perturbation theory now arises by ( 1) expanding the inverse occurring in the solution of eq. (7) CR=-(H*~-{H~~~}RR~)-‘IH~~]R
(11)
in a power series of the form (A-B)-‘=A-‘+A-‘BA-‘+...
,
where A is the matrix of the diagonal elements of H,- %rLr~l, (2) replacing the “sum over states” with a “sum over orbitals” and (3) replacing the diagonal elements of H,,{HeffJRR1 with orbital energy differences obtained from the diagonal elements of some reference oneelectron Hamiltonian hR. Note thatwehavethefreedom to use a difirent reference hR for each referencestateR in a “‘perturbthen diagonalize”formalism(i.e. wecan use a differentA/ B partitionfor each reference function in eq. (11)). It has been the choice of hR that has caused the difficulty with possible implementations of multireference MP theory. (See the discussion of the choice of hR in ref. [ 8 ] ). The use of a different reference hRfor each referencestateR and the use of an orthogonalvalencebond reference space is criticalto thesuccessof the present method.When the “sum over states” is replaced by a “sum over orbitals” one must distinguish three types of orbitals: inactive (or core), active (or valence) and virtual orbitals. For a given reference state R certain active orbitals are occupied while the others are unoccupied. In some other reference state S the roles of these same occupied/unoccupied orbitals will be completely reversed. Now we suppose that the occupied active orbitals in R have negative orbital energies while the unoccupied orbitals have positive energies. If the same reference h is used for reference state S then obviously the occupied active orbitals in S will have positive orbital energies while the unoccupied active orbitals have negative energies and the perturbation expansion may be hopelessly divergent. The preceding simple argument is enough to convince one that a low-order muiti-reference MP theory will only be successful if (a) the orbitals are quasi-degenerate and (b) a different hR is used for each reference state R. The use of a different hR for each reference state R is simple
8 July 1988
to implement. Requiring that the orbitals are quasidegenerate implies the use of a valence bond reference space. In previous work [9] we have shown that localization of the active orbitals of a CAS SCF wavefunction leads to an orthogonal VB expansion. Thus the CI expansion becomes a linear combination of atomic or fragment electronic structures (neutral and ionic). In this context, the obvious form of hR is the UHF form,
and
where 1.;1”and lffl include the inactive t “occupied” active in reference state R. For canonical MC SCF orbitals hg and h; will be approximately diagonal. In MP methods one uses the diagonal elements of hi and hj as orbital energies. For canonical orbitals we expect a substantial gap between the occupied and unoccupied orbital energies. However, for localized CAS SCF orbitals the orbital energies for a particular reference state R should be quasidegenerate (i.e. similar to those in the atom or fragment rather than in the molecule). The AER of cq. ( lob) become similar (identical at infinite interfragment distance) to sums of the correlation energies of neutral, anionic or cationic states of the isolated atomic or molecular fragments. However, they must of course exclude the correlation of two active “holes” with two active “particles” to avoid overcounting since this effect is included in the reference space. This idea has been followed by Ruedenberg et al. [ 141 in a semi-empirical manner for diatomic mOlcCUleSby extracting &??R from atomic computations. Of course, since hE and h$ are not diagonal in a basis of localized orbitals (except at infinite interfragment distance), there must be terms at third order that arise from this effect. We can now summarise an algorithm for multi-reference MP2 according to the scheme outlined above: ( 1) A CAS SCF computation in a basis of determinants (rather than configuration state functions). (2) Localization of CAS SCF orbitals [ 9 ] to give 185
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
an orthogonal valence bond expansion OVB. (3 ) If I {A~}R I 2 “threshold” then compute hz and h$ for reference configuration R. (4) If ({A,}, 12 “threshold” then compute AE, for reference configuration R using the usual MP2 formula excluding the “reference states” from the summation (i.e. the correlation of two active “holes” by two active “particles” is excluded since this state is in the reference space). (5) Compute the total energy accordig to eq. (10a).
3. Results and discussion As has become the accepted practise, we now discuss the results of OVB MP2 in comparison with some full CI computations [ 51. We have chosen the full CI computations of Bauschlicher and Taylor [ 15] on the singlet-triplet separation in CHr (doublezeta + polarization ) , the computations of Handy and co-workers [ 5,161 and Brown et al. [ 171 on Hz0 at (at the double-zeta level) r,, 1.5r,, and 2r, and the Be results of Heully and Daudey [ 17] since they are classic problems where SCF plus single reference MP2 fails. The CH2 problem is now well understood. Aside from basis set effects, the error in any single reference treatment (SCF, MP2 or SDCI) of the ‘A, state arises from the fact this state is a combination of two configurations - (a, ) * and a small amount of (b, ) *. On the other hand the ‘B, state has the configuration (a, ) ’(b, ) ’and is well described at the SCF level. The results are collected in table 1 along with those of Bauschlicher and Taylor [ 151. It should be noted that for the 3B, state, MC SCF is identical to restricted open shell SCF (ROSCF), and the MC SCF SDCI and OVB MP2 are equivalent to a restricted open shell. Similarity, the UHF result for the ‘A1 state is jsut ordinary closed shell. Of course the MC SCF and MC SCF SDCI results are in very good agreement with the full CI with respect to the ‘A1-3B1energy difference. While UHF MP2 is an improvement over the UHF result, it requires OVB MP2 to reduce the error to about the sort of level that one might expect for a single reference MP2 in an appropriate problem. The results on Hz0 are collected in table 2, along with the full CI results of Handy [ 161 and the MC 186
8 July 1988
Table 1 Comparison of OVB MP2 results with CI results for CH2 Method
‘A,
3B,
A()B,-‘A,) (kJ mol-‘)
UHF MC SCF UHF MP2 OVB MP2 MC SCF SDCI =’ exact (full CI) c,
-38.8863 =’ -38.9076 -39.0099 -39.0154
-38.9326 - 38.9276 b, - 39.0382 - 39.0366 b,
121.5 52.7 74.1 55.6
- 39.0222
-39.0416 b,
51.0
- 39.0272
- 39.0463
so. 1
a) Equivalent to RHF. b, The triplet 2 determinant (S,=O) MC SCF is identical to a restricted open shell SCF (ROHF). “Ref. [lS].
SCF SDCI results of Brown et al. [ 171. We have included the results at r,, 1.5r, and 2r,. Further we have added the results of a computation run at r,= 10 8, for comparison with the results at infinite r, of Brown et al. [ 171. In fact the reference space chosen in some of the computations is not balanced at infinite re so we include results at r,= 10 8, merely to indicate that our method behaves sensibly as r gets large. We have included results for two different (CAS SCF) reference spaces denoted in table 2 as “n in m” where n is the number of active electrons and m is the number of active orbitals. In each case the valence orbitals were localized using the procedure described in ref. [ 91 subject to the constraint that the 0 atom orbitals retained Czv symmetry (i.e. rather than a hybridized form). The “4 in 4” reference space included the 0 2p(b2), the 0 2p(a,) and the two H 1s orbitals with the 0 2p (b, ) lone pair and the 0 2s in the inactive space. This is the smallest space that can describe the dissociation into atoms. The corresponding MC SCF SDCI computations of Brown et al. [ 171 had no inactive space. In the “6 in 6” reference space the 0 2p( b, ) lone pair and a corresponding 0 3p (b, ) virtual orbital was added to the active space. The localized orbitals deserve some comment for the “6 in 6” reference space. In this case the two 0 p(b, ) orbitals localized into a “tight” orbital and a “diffuse” orbital where the diagonal elements of the density matrix were both approximately 1 (as opposed to the canonical delocalized situation where the occupancies were approximately 2 and 0).
Volume 148, number 2,3
CHEMICAL PHYSICS LEI-IERS
8 July 1988
Table 2 Comparison of MC SCF MP2 results with CI results for Hz0 Method
r.
I.%,
2.or,
1oA
SCF MP2 exact (full CI) 9)
-76.0098 -76.1493 -76.1579
-75.8035 - 75.9946 -76.0145
-15.5952 -75.8525 -75.9052
- 75:*648 b’
4in4”) MC SCF MC SCF SDCI d, OVB MP2
-76.0596 -76.1559 -76.1259
- 75.9244 - 76.0125 - 75.9843
-75.8273 -75.9033 - 75.8904
-15.7958 -75.8631 bi -75.8535
6 in 6 ” MC SCF MC SCF SCDI f’ OVB MP2
- 76.0970 -76.1575 -76.1442
- 75.9527 -76.0140 -76.0028
-75.8442 -75.9047 -75.8942
-75.8057 -15.8642 b’ -75.8573
a1 Refs. [5,16]. b, Ref. [ 171 result obtained at 100 A. =) Active space consisted of four electrons in the four O-H (H 2s and 0 2p) bonding orbitals. *) Ref. [ 171. In these calculations all the inactive orbitals (0 Is, 2s, 2p) were included in the reference space as well as the orbitals in c). e, To the active space in c) was added the additional 0 2p (lone pair) and a corresponding 0 3p orbital. ‘) Ref. [ 171. In these calculations all the inactive orbitals (0 ls, 2s) were included in the reference space as well as the orbitals in e).
This situation is essential in order to active a quasidegenerate reference. If the two 0 p (b , ) orbitals are not localized the 0 3p (b, ) virtual orbital has a high orbital energy and the quasi-degeneracy requirement cannot be satisfied. From table 2 it can be seen that while the OVB MP2 results for the “4 in 4” reference space are quite encouraging for r > 2r, the results for r < 2r, are rather poor. However, the addition of the two 0 p (b, ) orbitals in the “6 in 6” reference space improves the agreement of OVB MP2 with MC SCF SCDI and full CI considerably. In this case, while the absolute error at each value of r is still quite large ( x 10 rn& ), the relative error is similar at each point ( 10 m& and 2re and 14 n$, at re). Notice that the error in single reference MP2 at r, is 9 mEi,. Thus as r approaches r = co, where the OVB wavefunction is dominated by one term, we should expect an error of the same order of magnitude. We have thus accomplished our objective since the error in OVB MP2 is about 10 rnE,, at all values of r. The addition of the two 0 p(b, ) orbitals to the reference space at r, has a very large effect because at this geometry the OVB CI expansion is dominated by the 02-2H+ configuration so the additional atomic intra-atomic correlation is very important. The extension of the computations beyond this
point is not easy from a technical point of view. On the one hand, the addition of any further weakly occupied orbitals to the active space of the MC SCF is not without problems. As the occupation number of a weakly occupied orbital begins to approach 0.0, the MC SCF orbital rotation Hessian becomes almost singular. Thus the added weakly occupied orbital is not well determined and the localization step is not well determined either. Further, since such weakly occupied orbitals cannot be localized, they will have high orbital energies and the quasi-degenerate requirement becomes difficult to satisfy. Of course one could go to higher order (a) by computing the offdiagonal terms in H,, and (b) including third-order effects. It is of interest to comment on the semi-empirical scheme of Ruedenberg et al. [ 141 in relation to the present approach. In the approach of Ruedenberg et al., the AER (eq. ( lob) ) are replaced by correlation energies of neutral and ionic states of the constituent atoms. On the one hand this semi-empirical estimate does not allow for changes in the correlation energy of the atom when placed in the molecule. On the other hand, the correlation of two active “holes” by two active “particles” is overcounted in Ruedenberg’s scheme since this term is included in AERas well as in the reference space. In table 3 we have shown the 187
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
8 July I988
Table 3
Contribution of some important AER(eq. ( 1Ob) ) for Hz0 (4 in 4 CAS SCF) Configuration
&IIG%I
R r,
1.5r,
2.Or,
IOA
O('P)HH
-0.057/0.27
-O.OS5/0.38
-0.058/0.50
-0.058/1.0
0-H+H 02-H+H+
-0.064/0.27 -0.084/0.33
-0.068/0.21 -0.089/0.20
-0.084/0.12 -0.128/0.06
values of AER and 1{A,}, 1 for the most important OVB configurations for H20. As expected for an OVB
expansion, the 02-H+H+ configuration is very important at r,. Thus in order to describe the difference E(2r,) -E(r,) it is necessary to get AE(O*-2H+) correct. Thus the “6 in 6” reference space which includes the dominant 0 p (b , ) * electron correlation is essential. Note that while AE(03P 2H’S) is almost independent of r,, the value of AE(02-2H+) is strongly dependent on r. This fact seems to preclude the use of semi-empirical estimates. Finally, we give a brief discussion on the comparison of the results obtained with the method used in this work with the exhaustive study of Heully and Daudey [ 18 ] on the Be atom. Ordinary MP methods converge very slowly for Be because of the near degeneracy of the 2s and 2p orbitals. Further as demonstrated by Heully and Daudey [ 18 ] the multi-reference MP series diverges because the (2s) l (ns)’ state lies lower in energy than the (2p)* state. Heully and Daudey avoid this problem in an elegant way using an “intermediate” effective Hamiltonian. Our results are collected in table 4 along with the results of Heully and Daudey [ 181. We have used the same 6-3 1 1 + C ( 1d) basis as Heully and Daudey. The results of standard MP2-MP4 calculations indicate the expected slow convergence. The corresponding effective Hamiltonian results of Heully and Daudy [ 18 ] are also given. Their results are much improved over the usual MP2 at second order. It can be seen that the OVB MP2 result (CAS SCF “2 in 4” reference space) has a very small error indeed. The localized orbitals in this case were an exactly degenerate sp3-type combination of the CAS SCF active orbitals. Heully and Daudey [ 181 show that the multi-reference perturbation expansion is truly divergent in this example while the intermediate Hamiltonian (which gives only the Be ‘S state energy 188
/o.o IO.0
Table 4 Comparison of OVB MP2 results for Be( ‘S) with the MBPT results of Heully and Daudey [ 181 Method
E
“exact” (full CI ref. [ 181)
- 14.6445
0.0
this work SCF MP2 MP3 MP4 MC SCF OVB MP2
-
14.5720 14.6186 14.6298 14.6343 14.6156 14.6435
0.0725 0.0259 0.0147 0.0102 0.0289 0.0010
- 14.6251 - 14.6486
0.0194 -0.0041
Error
ref. 1181 Hen SC MBPT-2 Herr MC MBPT-2
exactly) is convergent. Thus, in order to extend the present computations beyond second order, one would be forced to use similar methods.
4. Conclusions A practical implementation of multi-reference MP2 depends upon the quasi-degeneracy of the active orbitals (which implies the use of an OVB expansion of the MC SCF wavefunction) and the use of a different reference one-electron Hamiltonian for each reference configuration. The results on Hz0 show the importance of the correlation energy of ionic fragments. Since we have implemented the method in the most obvious way, a detailed discussion of computational efficiency is impossible. However, at its simplest level computation times for OVB MP2 will be simply Nxsingle reference MP2 where N is the number of reference configurations. The MC SCF method is often criticized because
Volume 148, number 2,3
CHEMICAL PHYSICS LETTERS
the choice of reference space (i.e. active orbitals) may seem rather arbitrary. The transformation to an OVB expansion via the localization of the active orbitals [ 91, yields a wavefunction which is easily interpreted and thus demonstratably correct for a given problem. While the orthogonality of the orbitals in an OVB expansion yields a wavefunction with large ionic contributions, we have demonstrated how these contributions can be eliminated [ lo]. However, in the present method, these ionic terms cause no difficulties. The use of a non-orthogonal valence bond reference space would in fact destroy the very simple interpretation of the correlation energy at r, in H20 (see the previous discussion of the contribution AE(O’-2H+)).
Acknowledgement The authors are ‘grateful for the support of the SERC (UK) via grant No. GR/E 4499.4. The localized orbitals were computed using the algorithm described in ref. [ 91; however, the starting orbitals were in each case generated using the Boys method and the authors are grateful to Dr. A. Bottoni for his programs. Finally, we have used the GAUSSIAN [ 41 suite of programs in all our computations.
References [ 1 ] J. Goldstone, Proc. Roy. Sot. A 239 ( 1959) 267. [ 2 ] C. Meller and MS. Plesset, Phys. Rev. 46 ( 1934) 6 18.
8 July 1988
[ 31 B. Huron, P. Rancurrel and J.P. Malrieu, J. Chem. Phys. 73 (1973) 5745. [4] J.S. Binkley, R.A. Whiteside, R. Krishnan, R. Seeger, D.J. DeFrees, H.B. Schlegel, S. Topiol, L.R. Kahn and J.A. Pople, GAUSSIAN 80, QCPE 13 (1980); J, Binkley, M, Frisch, K. Raghavachari, D.J. DeFrees, H.B. Schlegel, R. Whiteside, E. Fluder, R. Seeger and J.A. Pople, GAUSSIAN 82, Carnegie-Mellon University. [ 51 N.C. Handy, Faraday Symp. Chem. Sot. 19 (1984) 17. [6] B.H. Brandow, Rev. Mod. Phys. 39 (1967) 771. [7] R.J. Bartlet, Ann. Rev. Phys. Chem. 32 (1981) 359; I. Lindgren, Physica Scripta 32 (1985); D. Mukhergec, Intern. J. Quantum Chem. Symp. 20 (1986) 409. [ 81 K. Wolinski, H. Sellers and P. Pulay, Chem. Phys. Letters 140 (1987) 225. [9] J.J.W. McDouall and M.A. Robb, Chem. Phys. Letters 132 (1987) 319. [lo] J.J.W. McDouall and M.A. Robb, Chem. Phys. Letters 142 (1987) 131. [ 111 Ph. Durand and J.P. Malrieu, Advaa. Chem. Phys. 67 (1987). [ 121M.R. Hoffman and J. Simons, J. Chem. Phys. 88 (1988) 993. [ 131 M.A. Robb, in: Computational techniques in quantum chemistry and molecular physics, eds. G.H.F. Diercksen, B.T. Sutcliffe and A. Veillard (Reidel, Dordtecht, 1975) pp. 435 t-f. [ 141M.M. Schmidt, M.T.B. Lam, S.T. Elbert and K. Ruedenberg, Theoret. Chim. Acta 68 ( 1985) 69. [ 151 C.W. Bauschlicher Jr. and P.R. Taylor, J. Chem. Phys. 85 (1986) 5510. [ 161N.C. Handy, Chem. Phys. Letters 74 (1980) 280; R.J. Harrison and N.C. Handy, Chem. Phys. Letters 95 (1983) 386; P. Saxe, H.F. Schaefer III and NC. Handy, Chem. Phys. Letters 79 ( 1981) 202. [ 171F.B. Brown, I. Shavitt and R. Shepard, Chem. Phys. Letters 105 (1984) 363. [ 181J. Heully and J. Daudey, J. Chem. Phys. 88 (1988) 1046.
189