6288
J. Phys. Chem. 1996, 100, 6288-6297
New Versions of Approximately Extensive Corrected Multireference Configuration Interaction Methods† La´ szlo´ Fu1 sti-Molna´ r and Pe´ ter G. Szalay* Department of Theoretical Chemistry, Eo¨ tVo¨ s Lora´ nd UniVersity, H-1518 Budapest, P.O. Box 32, Hungary ReceiVed: September 25, 1995; In Final Form: December 13, 1995X
In a previous paper (Szalay,P. G.; Bartlett, R. J. J. Chem. Phys. 1995, 103, 3600) we reviewed all approximately extensive corrected multireference configuration interaction (MR-CI) methods and also proposed a new method called the multireference averaged quadratic coupled-cluster (MR-AQCC) method. In this paper we continue the investigation of these methods and propose new variants by including the effect of (a) virtual orbitals in the EPV (exclusion principle violating) contributions, (b) redundancies due to the MR nature of the reference space. Results of test calculations on water, Be2, Be/H2, ozone and NO2/HF systems are presented which show that these new methods are able to describe potential energy surfaces very accurately even if small, but qualitatively correct, reference spaces are used.
1. Introduction methods1,2
have Though single reference coupled-cluster (CC) been proven to be the most accurate methods to describe the ground state of molecules,3 they fail if the single determinant wave function is not a good zeroth order description of the system. Generalization to multireference (MR) cases is not straightforward.4,5 On the other hand, truncated configuration interaction (CI) methods (where MR generalization is simple) suffer from not being extensive,6 meaning that energy and other properties do not scale properly with the size of the system. Since MR-CI is conceptually simple, much effort has been invested recently to correct CI for the extensivity error. In a previous paper7 we reviewed all these methods and also proposed a new method called MR-AQCC (multireference averaged quadratic CC)8 which is very similar to the successful MR-ACPF (multireference averaged coupled pair functional) method of Gdanitz and Ahlrichs9 but it is based on a different extensivity correction.10 MR-AQCC does not seem to fail in pathological cases when MR-ACPF, due to overcompensation of the extensivity error, does.7,8 The comparison in ref 7 was based on both theoretical analysis and numerical tests. The theoretical analysis showed that the essential differences among the methods can be classified as: (1) handling the EPV (exclusion principle violating)11 terms, (2) whether they consider or neglect redundancies due to the MR nature of the reference function, (3) choice of the zeroth order wave function (Φ0) and energy (E0), and (4) existence of an energy functional. The functional of the energy must be such that (a) it depends only on the wave function parameters [As discussed in the case of multireference CC methods,12 a functional can be always defined, if this requirement is neglected. In this case, however, additional stationary equations13 have to be solved. Therefore the cost of the calculations increases.] (b) Its stationary equations reproduce the equations characterizing the method exactly, and (c) at the stationary point it gives the energy. For more detailed discussion see ref 7. The importance of the functional lies in the calculation of energy derivatives:12,14,15 if such a functional exists, calculation of energy derivatives can † This paper is dedicated to Prof. Isaiah Shavitt, who had a great influence on this paper by his publications, personal discussions, and even by presenting encouraging results with our MR-AQCC method. X Abstract published in AdVance ACS Abstracts, March 1, 1996.
0022-3654/96/20100-6288$12.00/0
be performed via the generalized Hellmann-Feynman theorem.16 We believe that in practical applications this is very important, since it allows the easy and inexpensive calculation of energy derivatives. It turned out that the existence of the functional is determined by the way the other three factors listed above are handled.7 As for the third point, numerical tests showed that the choice of the zeroth order problem (using the same reference space) is much less important than the way the EPV terms are approximated: different choices of the former with equal handling of the EPV terms causes eventually different numerical (convergence) behavior, but the results are essentially the same. Handling the EPV terms represents the main difference among the methods, and the numerical tests showed that (a) those terms cannot be neglected and (b) crude approximations can cause a proliferation of the effect of EPV terms leading to unreliable results. There is only one method, the MR-CEPA (MR coupled electron pair approximation) by Ruttink et al.20 which takes the redundancies due to the MR nature of the reference function into account, but this method neglects the effect of the EPV terms, by using the so-called L-CPMET or CEPA(0) approximation. [This approximation was first introduced by C ˇ ´ızˇek under the name L-CPMET.17 It is equivalent to LCCD2,6 and D-MBPT(∞),18 but the name CEPA(0)19 is the most frequently used name in the literature.] In summary, we concluded in ref 7 that MR-ACPF and MRAQCC seem to be the best choice, since (a) in tests they perform very well and (b) a functional exists for them. Redundancies handled in MR-CEPA20 have a large effect, as well. This paper reports our newest results in search of the optimal approximation. Guided by the results of ref 7, we propose two modifications to MR-AQCC: (1) removing the electron pair concept from handling the EPV terms to allow distinguishing virtual orbitals and (2) including the effect of redundancies due to the MR nature of the reference function. Using test systems similar to those in ref 7, we will show that both effects are important and improvement is possible. In addition, we also focus on the choice of the reference space. Our goal is to find methods which work well with the smallest reference space possible. The smallest reference space must give a qualitatively correct description of the system at all geometries of interest. Such a space can be best constructed © 1996 American Chemical Society
Approximately Extensive MR-CI Methods
J. Phys. Chem., Vol. 100, No. 15, 1996 6289
by Pulay’s scheme21 first used for the UNO-CAS method.22 In this scheme the natural orbital occupations of the UHF wave function are calculated. Active orbitals are identified by occupations between 1.98 and 0.02, i.e., if the correlation from or to these orbitals is substantial. 2. The CEPA(0) Approximation In this section we are going to deal with the single reference case; generalization to multireference cases will be discussed later. The starting point of our derivations is the full-CI energy expression and equation for the single and double excited determinants:
restrictions are due to the fact that in the summations of the left-hand side occupied (or virtual) labels cannot coincide, since this would mean the excitation of more than one electron from (or to) the same orbital. These restrictions will be essential in the following discussion. a cd Using the cluster condition cacd ikl ≈ ci ckl in the single ab cd abcd excitation equation and cijkl ≈ cij ckl in the double excitation equation we arrive at
〈φia|HN|φ0〉 + ∑cck〈φia|HN|φck〉 + k,c “*i,a”
cklcd〈φia|HN|φklcd〉 + ∑ k>l,c>d
ciacklcd〈φ0|HN|φklcd〉 ) cia ∑ cklcd〈φ0|HN|φklcd〉 ∑ k>l,c>d k>l,c>d
〈φ0|HN|Ψfci〉 ) ∆E and
〈φia|HN|Ψfci〉
)
cia∆E
〈φijab|HN|Ψfci〉
)
cijab∆E
〈φijab|HN|φ0〉 + ∑cck〈φijab|HN|φck〉 + k,c
where φ0 is a Hartree-Fock wave function with energy E0. HN ) H - E0 is the operator for the correlation energy ∆E. The fullCI wave function is written as a combination of excited determinants (using intermediate normalization):
Ψfci ) φ0 + ∑ciaφia + i,a
∑
cijabφijab +
i>j,a>b
abc abc cijk φijk + ... ∑ i>j>k,a>b>c
Here indices i, j, k, ... (a, b, c, ...) represent occupied (virtual) orbitals in φ0. We expand the above expressions by inserting Ψfci:
∆E )
∑
cklcd〈φ0|HN|φklcd〉
(1)
k>l,c>d
〈φia|HN|φ0〉 + ∑cck〈φia|HN|φck〉 +
∑ k>l>m,c>d>e
k,c cde a cde cklm 〈φi |HN|φklm 〉
∑
cklcd〈φia|HN|φklcd〉 + k>l,c>d ) cia cklcd〈φ0|HN|φklcd〉 k>l,c>d
∑
∑ k>l>m,c>d>e
“*ij,ab”
cde ab cde cklm 〈φij |HN|φklm 〉+
“)i,a”
(2)
cde ab cde cklm 〈φij |HN|φklm 〉+
k>l>m,c>d>e
(3)
According to the Slater rules “*i,a”
cde a cde cda 〈φi |HN|φklm 〉 ) ∑ ckli 〈φ0|HN|φklcd〉 ∑ cklm k>l>m,c>d>e k>l,c>d
and
∑ k>l>m>n,c>d>e>f
“*ij,ab”
cdef cdef cklmn 〈φijab|HN|φklmn 〉
)
cklcd〈φ0|HN|φklcd〉
(4)
“)ij,ab”
cklcd〈φijab|HN|φklcd〉 + ∑ k>l,c>d
∑ cdef cdef cklmn 〈φijab|HN|φklmn 〉) ∑ k>l>m>n,c>d>e>f cijab ∑ cklcd〈φ0|HN|φklcd〉 k>l,c>d
∑
k>l,c>d
Rijab ) cijab k,c
cijabcklcd〈φ0|HN|φklcd〉 ) ∑ k>l,c>d cijab ∑ cklcd〈φ0|HN|φklcd〉 k>l,c>d
We can apply the cluster condition in the double-triple matrix element in the double excited equation as well and then, as usual,23 neglect this term, since it depends on the single excitation coefficient and is therefore higher order in perturbation theory. Now one can easily recognize that in both equations the righthand side cancels the last terms of the left-hand side completely, but some terms remain on the right-hand side because of the restricted summations on the left-hand side. These remaining terms consist of contributions where occupied and/or virtual indices coincide (denoted by “) i, a” and “) ij, ab”)
Ria ) cia
and
〈φijab|HN|φ0〉 + ∑cck〈φijab|HN|φck〉 +
cklcd〈φijab|HN|φklcd〉 + ∑ k>l,c>d
cdab cklij 〈φ0|HN|φklcd〉 ∑ k>l,c>d
We denote by “* i, a” and “* ij, ab” that the summation indices cannot coincide with i or a and i, j or a, b, respectively. These
cklcd〈φ0|HN|φklcd〉 ∑ k>l,c>d
(5)
These terms are often called EPV (exclusion principle violating)11 terms, but (as, for example, our derivation shows clearly) there is nothing unphysical here. The simplest approximation is that one neglects the EPV terms. This is the L-CPMET17 or LCCD2,6 or MBPT(∞)18 or CEPA(0)19 approximation. This approximation yields very similar equations to the CISD equations, but instead of the correlation energy, we have zero on the right-hand side. Therefore, the method is extensive and also includes the effect of higher excitations, although they do not show up explicitly. Computationally, it costs about the same as the corresponding CISD calculation. One way of generalizing the single reference CEPA(0) method to the MR case is that we use the result of the SR derivation in the MR case and we simply delete the correlation energy from the right-hand side of the MR-CI equations. The MR-LCCM (linearized coupled-cluster method) of Bartlett and co-workers24 was the first method using this approximation. In the actual derivation they inspected the linearized form of the MR-CC equations.25 Later, Gdanitz and Ahlrichs introduced the MR-CEPA(0) method,9 which differs from MR-LCCM only
6290 J. Phys. Chem., Vol. 100, No. 15, 1996
Fu¨sti-Molna´r and Szalay
in handling the orthogonal complement of the reference space. For detailed comparison of these and other related methods26-29 see ref 7. The numerical experience showed that neglecting the EPV terms (eqs 4 and 5) causes an overestimation of the effect of higher excitations. Therefore, it is necessary to calculate EPV terms exactly or to approximate them. 3. Beyond the CEPA(0) Approximation As Daudey et al.30 showed recently, it is possible to calculate EPV terms exactly in the single reference case (the so-called self-consistent size-consistent CI ((SC)2 CI) method). In the MR case the exact handling of the EPV terms is not straightforward, but Malrieu and co-workers have already showed two ways of handling this problem.30,31 Another approximation has been published recently by Fink and Staemmler under the acronym MC-CEPA.32 As we see from eqs 4 and 5, the EPV corrections are different for all excitations and depend on a subset of the expansion coefficients. Therefore, as discussed in ref 7, one cannot construct the energy functional in this case. For this reason, we prefer those methods where the EPV terms are approximated by averaged quantities and therefore ensure the existence of the energy functional. As in the previous section, we consider the SR case and then apply the results for the MR equation. This procedure is obviously an approximation, and we will go beyond it in the next section. First we define the usual pair energy as
As of the EPV terms of the single excitation equation (eq 4), MR-ACPF is using the same formula:
∆E Ria ) cia ne/2
although the above derivation is strictly valid for the double excitations only. Though the approximation of the EPV terms in this form is rather simple, MR-ACPF represents a substantial improvement over the CEPA(0) approximation. The MR-AQCC Method.7,8 The multireference averaged quadratic coupled-cluster (MR-AQCC) method uses a slightly different approximation and tries to account for the interaction of the electrons through another averaging procedure.10 The electrons are divided into (n2e) equivalent electron pairs:
k,l ≈ j )
Rijab ≈ cijab
(( ) ( )) (
Ria ) cia 1 -
k 2 or l1 + l2 > 2 the resulting operator produces higher than double excitation. For these EiEj products the last term of the left-hand side of eq 13 completely cancels the corresponding terms on the right-hand side. If k1 + k2 e 2 and l1 + l2 e 2, Ej will not show up on the left-hand side; therefore, the corresponding term remains on the right-hand side. One can characterize all possibilities by a 9×9 matrix with rows and columns corresponding to the excitation classes: valid combinations can be marked by 0, others by 1. By this, it is easy to find whether a combination is valid or not; we just have to know which class Ei and Ej belong to. The definition of classes and, consequently, the subset of Ej operators depend on the representation of the space spanned by φi’s. Ruttink et al.20 used a determinental basis, while we prefer configurational space for conceptual and programtechnical reasons. As noted by Ruttink et al.,20 this has only a small effect on the results (see also Table 1 and 5). We define the quantity
(k,l) )
∑ cj〈φ0|HN|φj〉
j∈(k,l)
which is simply the energy contribution from excitation class (k, l). Now we can write the right-hand side of the i-th equation (φi ∈ (k1, l1)):
Ri ) ci( ∑ P[(k1,l1),(k2,l2)] (k2,l2))
(14)
(k2,l2)
Here the P[(k1,l1),(k2,l2)] is the elements of the 9×9 matrix, and its elements can be one or zero according to whether the two classes give or do not give a contribution to Ri. This is the essence of the MR-CEPA method. It does not consider the EPV terms, just the redundant excitations as described above. Therefore, in the single reference case the method is identical to CEPA(0).
Now we are going to consider the EPV terms as well. In this procedure we get a modified version of the ACPF, AQCC, and AQCC-v methods, i.e., we will use a similar averaged counting to get new values for the P matrix as used in the original methods. The most logical generalization can be done for the MRAQCC-v method. Besides the redundancy, we want to consider the possible coincidence of the indices in classes (k1, l1) and (k2, l2). Take as the first example the classes (2, 2) and (2, 2) with a product belonging to quadruple excited functions which does not give a contribution to Ri in the MR-CEPA case. Since this setup is exactly the same as in the SR case, we recognize immediately that both the double occupied and the virtual orbitals can coincide and the corresponding P element will differ form zero:
(
)( ) ( )( )
No - 2 Nv - 2 2 2 P(2, 2; 2, 2) ) 1 No Nv 2 2
As another example take the classes (2, 2) and (2, 1). Here in the second class one excited electron goes into the active space; therefore, only one virtual index can be in coincidence:
(
)( ) ( )( )
No - 2 Nv - 2 2 1 P(2, 2; 2, 1) ) 1 No Nv 2 1
As a third example consider the classes (1, 1) and (1, 1). The product of these two excitations is a double excitation, i.e., the whole term survives on the right-hand side:
P(1, 1; 1, 1) ) 1 as in MR-CEPA. The general formula can be given as
P(a, b; c, d) ) No - a Nv - b na - f(b-a) (Na - na)f(a-b) c d f(d-c) f(c-d) 1na No Nv (Na - na) c d f(d-c) f(c-d)
(
)(
)( ( )( )(
)(
)(
)
)
where no is the number of active electrons and f(x) ) x if x g 0 and f(x) ) 0 otherwise. Note that the denominator is the number of all possible excitations in class (c, d) and the numerator is the number of those excitations in class (c, d), which are allowed together with an excitation in class (a, b). In the case of MR-AQCC, we do not consider the coincidence of the virtual labels. Therefore, for classes (2, b) and (2, d), where b and d can be 0, 1, or 2, we have
(
No - 2 2 P(2, b; 2, d) ) 1 No 2
( )
)
which is exactly the MR-AQCC factor in the SR case. For any pair of classes in form (2, b) and (1, d)
(
No - 2 1 P(2, b; 1, d) ) 1 No 1 The general formula can be given as
( )
)
Approximately Extensive MR-CI Methods
(
)( ( )(
J. Phys. Chem., Vol. 100, No. 15, 1996 6293
No - a na - f(b-a) c f(d-c) P(a, b; c, d) ) 1 na No c f(d-c)
)
)
We use the acronym MR-AQCC-vmc and MR-AQCC-mc for the two methods, respectively. “mc” refers to the fact that redundancies due to the multireference φ0 are considered. (In this nomenclautre, MR-CEPA of Ruttink et al.20 could be called MR-CEPA(0)-mc.) Generalization of MR-ACPF is not straightforward. The reason is that the “electron pair” picture can hardly be combined with the concept of excitation classes. To be able to test the effect of the redundant terms, we tried to find an acceptable generalization. We distinguish four cases. First, if the product EiEj is not in the R space, we have
P(a, b; c, d) ) 1 Here a + c e 2 and b + d e 2 must be fulfilled. If this is not the case, and at the same time any of a or c is zero (i.e. the electrons are not excited from the same space), EPV terms are not possible and we have
P(a, b; c, d) ) 0 for a ) 0 or c ) 0, but a * c. If Ei is a double excitation, we get the MR-ACPF formula:
ne - 2 2 ) P(a, b; c, d) ) 1 ne ne where a + c > 2 and/or b + d > 2, a ) 2 and/or b ) 2. Finally, if Ei is a single excitation, we have ne - 1 electrons to put into pairs and therefore we get:
P(a, b; c, d) ) 1 -
ne - 1 1 ) ne ne
where a + c > 2 and/or b + d > 2, a e 1 and c e 1. We refer to this method as MR-ACPF-mc. The discussion above showed how the redundant terms due to the MR nature of the reference wave function can be handled in the MR-AQCC, MR-AQCC-v, and MR-ACPF methods. Two main differences occur due to these changes. First, the energy functional does not exist, since the approximate form of the EPV contribution as given by eq 14 does depend on a subset of ci coefficients. As discussed above and also in ref 7, the functional is very important in practical applications. From this point of view we still prefer the original methods (MR-AQCC, MR-AQCC-v, and MR-ACPF), but the effect of the redundancy should be tested in practical applications. This is done in the following section. The second difference lies in handling the single excitations: the original methods and the corresponding one marked by “mc” are not equivalent if applied to a single reference case. This sounds a bit illogical, since in the generalization we considered additional redundancies due to the MR nature of the reference function, which is not present in the SR case. As mentioned earlier, the derivation of the original methods deals only with the double excitations and the approximate EPV contribution is applied for the single excitations, as well. This was necessary to save the functional form. In the mc methods, on the other hand, we can separately handle the single excitation, since we do not have the functional, anyway. Moreover, the logic of classifying excitations into classes automatically allows the derivation of the different
TABLE 1: Relative Energies for Water with the DZP Basis Using a 4×4 CAS Reference Spacea geometry (hartree)36
full-CI MR-CI9 MR-CEPA(0) MR-CEPAb MR-ACPF9 MR-ACPF-mc MR-AQCC MR-AQCC-mc MR-AQCC-v MR-AQCC-vmc
1Re
1.5Re
2Re
-76.256624 4.98 -2.06 -0.90 -0.12 0.13 1.41 1.49 1.78 1.76
-76.071405 4.52 -2.41 -0.76 -0.22 0.23 1.21 1.51 1.55 1.73
-75.952269 3.71 -2.59 -0.58 0.18 0.32 1.23 1.45 1.49 1.62
Energies relative to full-CI in millihartrees. Re )1 Å, ∠HOH ) 104°. b With the original MR-CEPA procedure using a determinental basis the relative energies are -0.97, -0.84, and -0.59 for the three geometries, respectively.20 a
factors of the single excitation equations. In this sense, too, the mc methods are less approximate. Again, numerical tests below will show how important these effects are. In this section we generalized MR-AQCC-v, MR-AQCC, and MR-ACPF by considering the redundancy due to the MR nature of the reference space. We would like to stress again that the approximation used in MR-AQCC-v is the most suitable for this generalization. 5. Numerical Tests Several approximations were used in the methods discussed in the previous section. Although we tried to introduce those approximations as carefully as possible, it is not entirely clear how they affect the results. The interaction of some of the approximations is especially unpredictable. Therefore, we decided to test the methods on examples where full-CI or other high-level calculations are available. As in refs 7 and 8, the tests include the stretching of both bonds in water, the potential curve for Be2, insertion of Be into H2, and the equilibrium geometry and harmonic frequencies of ozone. Finally we have tested the extensivity of the methods. Of course we are not going to repeat the conclusion of ref 7, we rather update the calculations by comparing to recent full-CI results and by focusing on the following questions: (1) what is the effect of the coincidence of virtual orbitals as included in MRAQCC-v, (2) what is the effect of the redundancy due to the MR nature of the reference function as included in MR-CEPA20 and in our mc methods, and (3) how large is the necessary reference space. All calculations have been performed by the COLUMBUS MR-CI program system.33 For ozone the analytical first derivative of the energy has been calculated according to the method described in ref 34. To obtain the UHF natural orbital occupations, the TX90 program of Pulay35 has been used. Bond Stretching in Water. Simultaneous stretching of both OH bonds in water requires a four-orbital/four-electron complete active space (CAS) (1a212a211b21(3a14a11b22b2)4) including the bonding and antibonding orbitals. The results of a full-CI calculation are available with a DZP basis in a study by Bauschlicher and Taylor,36 where the lowest occupied (1a1) orbital is frozen, and with the cc-PVDZ basis by Olsen et al.37 for all electrons. We performed calculations with both bases using the same geometries (see Tables 1 and 2) as in refs 36 and 37. The number of electrons was chosen to be ne ) 8 for the corrections in both calculations. Table 1 shows the relative energies with respect to full-CI at different distances for the DZP basis. As for the effect of redundancies due to the reference functions, one can clearly
6294 J. Phys. Chem., Vol. 100, No. 15, 1996
Fu¨sti-Molna´r and Szalay
TABLE 2: Relative Energies for Water with the cc-PVDZ Basis Using a 4×4 CAS Reference Spacea
a
geometry
1Re
1.5Re
2Re
2.5Re
3Re
full-CI (hartree)37 MR-CI MR-CEPA(0) MR-CEPA MR-ACPF MR-ACPF-mc MR-AQCC MR-AQCC-mc MR-AQCC-v MR-AQCC-vmc CCSDTQ37
-76.241860 4.96 -1.80 -0.79 0.05 0.22 1.52 1.56 1.89 1.83 0.02
-76.072348 4.72 -2.00 -0.57 0.07 0.43 1.47 1.69 1.82 1.92 0.13
-75.951665 3.72 -2.37 -0.54 0.23 0.34 1.28 1.45 1.54 1.64 -0.45
-75.917985 3.14 -0.92 -0.62 0.20 0.20 1.07 1.23 1.29 1.40
-75.911946 3.01 -0.82 -0.64 0.20 0.18 1.03 1.18 1.24 1.34
Energies relative to full-CI in millihartrees. Re ) 1.84345 bohr, ∠HOH ) 110.6°.
Figure 1. Water molecule (4×4 CAS reference, DZP basis): Shifted E-E(FCI) at different distances.
conclude that it increases the absolute energy in all cases. Since the methods are not variational (i.e., full-CI is not lower bound), this does not mean immediately that the description is worse. In contrast, we concluded in ref 7 that MR-CEPA(0) and even MR-ACPF overshoot the effect of the higher excitations; therefore, this correction by decreasing the energy acts in the right direction. Since AQCC and AQCC-v do not overshoot full-CI, AQCC-mc and AQCC-vmc give a larger absolute error than AQCC and AQCC-v. The absolute value of the error is not the most important test, however, since our goal is to get potential energy surfaces, which are parallel to the full-CI one. In Figure 1 we shifted the curves of relative energies obtained by the different methods to the same origin at Re, since the resulting graph is the best choice to present the error in parallelism. The figure shows that the error in parallelism (the sum of the absolute value of the largest positive and negative error) for MR-AQCC and MR-ACPF is about the same (a few tenths of a millihartree). The improvement due to the mc correction is substantial, and this error reduces to less than 0.05 mhartree for the MR-AQCC-mc and MR-AQCC-vmc methods. In the case of MR-ACPF-mc we have an improvement with respect to MR-ACPF, as well. Comparing MR-AQCC and MRAQCC-v results, the former is somewhat better both with and without the mc correction. Table 2 shows the results obtained with the cc-PVDZ basis. With this basis the absolute errors are somewhat larger for all methods at 1Re and 1.5Re but about the same at 2Re. This is a rather unusual basis set effect considering that the two basis sets are of the same size. We suspect that this difference is related to the unusually large HOH bond angle obtained with the cc-PVDZ basis. Interestingly, the increase of the absolute error corrects the error in parallelism for MR-ACPF, while making it worse for the AQCC methods. Nevertheless, this
error is smaller than 0.5 mhartree, or less than 2% of the absolute energy change between 1Re and 2Re. Olsen et al.37 also report the results for several single reference coupled-cluster methods using both RHF and UHF references. We included the most advanced RHF-CCSDTQ results in Table 2. The error of parallelism of the RHF-CCSDTQ method is smaller than that of the MR-CI method, but the MR-ACPF and MR-AQCC methods are already clearly better. Since CCSDTQ is very expensive, the performance of the corrected MR-CI methods is encouraging. Potential Energy Surface of Be2. In a previous paper7 we found good agreement between the MR-AQCC potential energy curve using a DZP basis and 2×2 CAS reference and the fullCI one at the available geometries (near equilibrium).38 At larger distances we could compare to CCSDT39 and the agreement was satisfactory, as well. In the meantime, a new full-CI curve with a more diffuse basis (the generally contracted basis of Widmark et al.40 using a 3s2p1d contraction, WMR) has been published by Nebot-Gil et al.41 We performed calculations with this basis at the geometries published in ref 41 (four core electrons frozen), and the results were rather disappointing for us: there was a good agreement of the MR-AQCC and MR-AQCC-v methods near equilibrium, but the depth of the second minimum41 [The second minimum is not physical but rather a basis set effect. Since the full-CI curve shows this second minimum, the approximate methods should, as well.] was highly overestimated by all methods. This effect is the most extreme for MR-ACPF, since absolutely no minimum is found.7 Therefore we performed the UHF natural orbital occupation analysis according to Pulay.21 and found that a qualitative description requires a 4×4 CAS including the highest occupied and lowest unoccupied σg and σu orbitals. Therefore, we repeated the calculations using a 4×4 CAS reference space. The potential energy curves are given in Figure 2. The relative energies with respect to full-CI are summarized in Table 3, and the corresponding error curves (shifted to the same origin at 5 au) are plotted in Figure 3. Redundancy due to the MR nature of the reference space now lowers the energy for all but the MR-CEPA method which overestimates full-CI considerably. Though it is the curve of the MR-ACPF method which is the closest to full-CI, both MR-ACPF and MR-ACPFmc describe the very small barrier and second minimum rather purely, since the error in parallelism of the MR-ACPF curve is large for the entire geometry range. On the other hand MRAQCC and MR-AQCC-v are extremely parallel to full-CI up to 5.5 au, while the mc methods are very accurate for the whole geometry range, and the error is smaller than 0.01 mhartree. We are aware that a 4×4 CAS for this system is large, since all electrons are already correlated in the reference space, but it seems to be conclusive that the AQCC methods are much
Approximately Extensive MR-CI Methods
Figure 2. Potential energy surface of a Be2 molecule (4×4 CAS reference, WMR basis).
more accurate for this system than the other methods. Therefore, we feel encouraged to investigate the Be2 potential energy surface using larger basis sets and focus on the physical problem rather than testing the methods. Those results will be published in a subsequent publication. The Be/H2 System. Insertion of Be into H2 is often used as a model to test the performance of both single reference and multireference methods.42 As in the full-CI calculation,42 we used a 3s1p and 2s basis for Be and H, respectively.43 Although no electron was frozen in the calculations, we used ne ) 4, since the correlation of the the core electrons is not important,9 and therefore they should not be included in the averaging procedure.7,9 In ref 7 we have found a large error for MR-AQCC using a 2×2 CAS reference (1a212a213a21 and 1a212a211b22 configurations), especially at the third geometry. We were interested to see whether inclusion of virtual orbitals and the effect of the MR nature of the reference space on the EPV correction improve the results. As Table 4 shows, MR-AQCC-v is slightly worse than MR-AQCC. As for Be2, the mc correction (except in MCCEPA) decreases the energy and now MR-AQCC-mc performs very well for the third geometry, as well. The results of the UHF natural orbital test21 are listed in Table 5. It shows the importance of the 1b2 and 3a1 orbitals in the correlation, since their occupation changes considerably with the geometry. According to Pulay’s suggestion,21 the 2a1 and 2b2 orbitals should also be included in the reference space, resulting a 4×4 CAS. The results with the 4×4 CAS reference are listed in Table 6. We find excellent agreement for all AQCC methods with full-CI. The results of ACPF are improved as well, but the error for the third geometry is larger. Ozone. It turned out that ozone is difficult to describe with small reference spaces.8 According to the UHF natural orbital test,21 the qualitative description requires two configurations:
Φ1 ) ...6a214b221b211a22 Φ2 ) ...6a214b221b212b21 Table 7 shows the equilibrium geometries and harmonic vibrational frequencies of ozone with the different methods. In the calculations we used the same DZP basis as in refs 8 and 44, and the three lowest occupied orbitals were kept frozen. Therefore, the number of correlated electrons was ne ) 18. As already concluded in ref 8, MR-ACPF overestimates the effect of higher excitations resulting in too long bonds and too low stretching frequencies. MR-AQCC improves the results for both
J. Phys. Chem., Vol. 100, No. 15, 1996 6295 geometry and vibrational frequencies considerably, while MRAQCC-v gives excellent agreement with CCSDT.44 We did not test the mc methods in this example, since analytic energy derivatives are not available for them and the calculation from single energy points would be tedious. This fact underlines again that the existence of the energy functional and therefore the cheap calculation of the energy derivatives are very important in practical applications. Extensivity. Because of the approximations, the multireference methods discussed above are not rigorously extensive, and only numerical tests can provide an idea about their extensivity error. As usual, we calculate the extensivity error of the noninteracting subsystems, although it must be understood that this error is present in interacting systems (molecules) as well. This latter fact can be well demonstrated by approaching the noninteracting systems toward each other: the interaction energy will grow continuously, and there is no reason why the extensivity error should disappear at infinitely small interaction energy. As in ref 7, we have chosen N2O plus HF system45 as a test and used two basis sets:the 6-31+Gdp basis used in ref 45 and a DZP basis.46 The reference space was SR on HF and a 4×4 CAS on N2O. The results are summarized in Table 8. Both the inclusion of virtual orbitals and the mc effect cause a decrease of ∆E ) Edimer - EHF - EH2O, i.e., the extensivity error. In case of MR-ACPF-mc this means that it does not overestimate the extensivity error, as MR-ACPF does. On the other hand, extensivity error of the mc methods becomes a bit larger, while MR-AQCC-v has an even larger error. In any event the extensivity error of all AQCC and ACPF methods is of the same magnitude, and its absolute value is 1.5-2 orders of magnitude smaller than that of MR-CI. Therefore, all methods can be classified as “approximately extensive”. 6. Conclusions In this paper we continued the investigation of approximately extensive corrected MR-CISD methods. In the theoretical section we proposed new versions of these methods by (a) including the effect of virtual orbitals in the EPV contribution (MR-AQCC-v method), (b) including the effect of redundancies due to the MR nature of the reference space similarly to the MR-CEPA method of Ruttink et al.20 does (MR-ACPF-mc, MRAQCC-mc, and MR-AQCC-vmc methods). Numerical test calculations on water, Be2, the Be/H2 system, ozone, and NO2/HF system lead us to the following conclusions: (1) inclusion of the EPV terms due to the coincidence of virtual orbitals in an averaged way, as presented in this paper, has a very small effect. Since this small change in some cases acts in the wrong direction, we prefer the original MR-AQCC over MR-AQCC-v, (2) the sign of the energy change due to the mc correction varies for the different systems. The parallelism of the approximate curves to the full-CI improves by use of the mc correction especially for larger bond lengths, (3) contrary to previous attempts (see ref 7), qualitative description of both Be2 and the Be/H2 system requires a 4×4 CAS reference space, and (4) as already stated in our previous paper,7 MR-AQCC is preferred to MR-ACPF in some difficult cases, but normally they give almost identical results. The same conclusion can be drawn from the results of this paper. It was only the water calculation using the cc-PVDZ basis, where MRACPF (and also MR-ACPF-mc) performed slightly better than the MR-AQCC methods. We think that the MR-AQCC methods presented here and in previous papers7,8 can provide a very good approximation to
6296 J. Phys. Chem., Vol. 100, No. 15, 1996
Fu¨sti-Molna´r and Szalay
TABLE 3: Relative Energies for Be2 with the WMR Basis Using a 4×4 CAS Reference Spacea rBeBe (au)
full-CI41 (hartree)
4.0 4.4 4.8 5.0 5.2 5.6 6.0 7.0 8.0 9.0 10.0 12.0 14.0 16.0
-29.229 824 -29.235 853 -29.237 368 -29.237 467 -29.237 414 -29.237 227 -29.237 140 -29.237 231 -29.237 255 -29.237 158 -29.237 044 -29.236 905 -29.236 853 -29.236 834
a
MR-CI MR-CEPA(0) MR-CEPA MR-ACPF 3.20 3.22 3.27 3.28 3.28 3.24 3.19 3.06 2.98 2.95 2.93 2.93 2.93 2.93
-5.85 -6.24 -6.50 -6.49 -6.38 -6.02 -5.66 -5.10 -4.87 -4.77 -4.73 -4.70 -4.69 -4.69
-3.76 -3.97 -4.31 -4.48 -4.62 -4.79 -4.81 -4.60 -4.42 -4.34 -4.30 -4.27 -4.27 -4.26
0.21 0.12 0.00 -0.05 -0.11 -0.20 -0.24 -0.24 -0.21 -0.19 -0.19 -0.18 -0.18 -0.18
MR-ACPFMR-AQCCMR-AQCCmc MR-AQCC mc MR-AQCC-v vmc 0.13 0.03 -0.08 -0.14 -0.20 -0.28 -0.32 -0.29 -0.26 -0.24 -0.23 -0.23 -0.22 -0.22
2.28 2.27 2.27 2.26 2.24 2.19 2.13 2.04 2.00 1.98 1.97 1.97 1.97 1.97
1.98 1.97 1.97 1.97 1.96 1.95 1.94 1.93 1.93 1.92 1.92 1.92 1.92 1.92
2.36 2.35 2.35 2.34 2.32 2.27 2.22 2.12 2.08 2.06 2.05 2.05 2.05 2.05
2.07 2.05 2.06 2.06 2.06 2.05 2.04 2.02 2.01 2.00 2.00 2.00 2.00 2.00
Energies relative to full-CI in millihartree.
TABLE 6: Relative Energies for the Be/H2 System Using a 4×4 CAS Referencea r(BeH2) r(H2)
2.50 2.78
2.75 2.55
3.00 2.32
full-CI (hartree)42 MR-CI MR-CEPA(0) MR-CEPA MR-ACPF MR-ACPF-mc MR-AQCC MR-AQCC-mc MR-AQCC-v MR-AQCC-vmc
-15.62288 0.04 -0.13 -0.04 -0.01 -0.01 0.02 0.01 0.03 0.02
-15.60292 0.08 -0.93 -0.35 -0.14 -0.14 0.01 -0.05 0.03 -0.01
-15.62496 0.18 -1.29 -1.11 -0.42 -0.43 -0.01 -0.05 0.04 0.03
a
Energies relative to full-CI in millihartrees.
Figure 3. Potential energy surface (shifted relative energies) of a Be2 molecule (4×4 CAS reference, WMR basis).
TABLE 7: Ground State Properties of Ozone Using a 2×2 CAS Reference Space and the DZP Basisa
TABLE 4: Relative Energies for the Be/H2 System Using a 2×2 CAS Referencea
method
r(BeH2) r(H2)
2.50 2.78
2.75 2.55
3.00 2.32
full-CI (hartree)42 MR-CI24 MR-CEPA(0)9 MR-CEPAc MR-ACPF9 MR-ACPF-mc MR-AQCC MR-AQCC-mc MR-AQCC-v MR-AQCC-vmc
-15.62288 0.84 -3.28 -1.54 -0.90 -0.93 0.29 -0.19 0.45 0.02
-15.60292 2.01 4.30b -1.94 -0.90 -0.76 1.11 0.57 1.37 0.87
-15.62496 3.08 -5.50 -4.81 -0.53 -2.20 1.98 0.36 2.29 1.11
a
rOO ∠OOO ω1 ω2 ω3 Eb
MR-CI
MRACPF
MRAQCC
MRAQCC-v
CCSDT44
1.261 116.5 1235 761 1338 -0.871296
1.315 115.8 930 640 898 -0.942051
1.292 116.1 1070 694 1070 -0.928506
1.287 116.2 1100 705 1110 -0.923854
1.286 116.7 1141 705 1077 -0.941197
a Bond lengths in angstroms, bond angles in degrees, harmonic frequencies in wavenumbers. b Total energy +224 hartrees.
TABLE 8: Extensivity Error of Different Methods on the N2O/HF System Using a 4×4 CAS Referencea basis
b
Energies relative to full-CI in millihartrees. Saddle point, see ref 9. c Using the original MR-CEPA procedure with a determinental basis the relative energies are -1.65, -2.55, and -5.88 for the three geometries, respectively.20
MR-CI MR-CEPA(0) MR-CEPA MR-ACPF MR-ACPF-mc MR-AQCC MR-AQCC-mc MR-AQCC-v MR-AQCC-vmc
TABLE 5: UHF Natural Orbital Occupations for the Be/H2 System r(BeH2) r(H2)
2.50 2.78
2.75 2.55
3.00 2.32
2a1 1b2 3a1 2b2
1.97 1.41 0.59 0.03
1.95 1.11 0.89 0.05
1.95 0.80 1.20 0.05
full-CI in terms of parallelism of the potential energy curves. The requirement is that the reference space is able to qualitatively describe the system at all geometries. To find this minimal reference space, the inspection of UHF natural orbital occupation can be used as proposed by Pulay.21,22 In this way
a
DZP
6-31+Gdp
0.031 518 -0.086 055 -0.000 097 -0.000 100 0.000 126 0.000 247 0.000 349 0.000 879 0.000 955
0.032 695 -0.044 842 0.000 161 -0.000 327 0.000 140 0.000 002 0.000 290 0.000 539 0.000 835
Edimer - EN2O - EHF in hartrees.
a “black box” type algorithm can be set up to construct the wave function for the MR methods. Though MR-AQCC-mc and MR-AQCC-vmc show better overall performance then MR-AQCC, the latter method performs very well near the equilibrium geometry. Since analytic energy derivatives are available for MR-AQCC and not for the mc methods, we prefer the use of MR-AQCC for properties near
Approximately Extensive MR-CI Methods equilibrium (geometry, vibrational frequencies, etc.). [Even if the formulas of analytical gradients were derived and programmed for the mc methods, they will cost more (see the introduction).] For potential energy surfaces, where single point energy calculations are sufficient, MR-AQCC-mc (and eventually MR-AQCC-vmc) can be used. Acknowledgment. The paper is using parts of the masters thesis of L.F.M.. The authors thank Prof. R. J. Bartlett for initiating this project and Prof. J. H. van Lenthe for providing outputs of their MR-CEPA calculations. This work was supported by the Hungarian Scientific Research Foundation (No. F-013962) and by a grant of the Department of Education (No. 333/94). References and Notes (1) Coester, F. Nucl. Phys. 1958, 1, 421. Coester, F.; Ku¨mmel, H. Nucl. Phys. 1960, 17, 477. C ˇ ´ızˇek, J. J. Chem. Phys. 1966, 45, 4256; AdV. Chem. Phys. 1969, 14, 35. Bartlett, R. J. J. Phys. Chem. 1989, 93, 1697; Annu. ReV. Phys. Chem. 1981, 32, 359. (2) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910. (3) Bartlett, R. J. J. Phys. Chem. 1989, 93, 1697 and references therein. (4) Paldus, J. In Methods of Computational Molecular Physics; Wilson, S., Dircksen, G. H. F., Eds.; NATO Advanced Study Institute Series; Plenum: New York, 1992; Vol. 4, p 99. (5) Paldus, J. In RelatiVistic and Electron Correlation Effects in Molecules and Solids; Malli, G. L., Ed.; NATO Advanced Study Institute Series; Plenum: New York, 1994; p 207. (6) Bartlett, R. J.; Purvis, G. D. Int. J. Quantum Chem. 1978, 14, 561. Bartlett, R. J.; Purvis, G. D. Phys. Scr. 1980, 21, 251. (7) Szalay, P. G.; Bartlett, R. J. J. Chem. Phys. 1995, 103, 3600. (8) Szalay, P. G.; Bartlett, R. J. Chem. Phys. Lett. 1993, 214, 481. (9) Gdanitz, R.; Ahlrichs, R. Chem. Phys. Lett. 1988, 143, 413. (10) Meissner, L. Chem. Phys. Lett. 1988, 146, 205. (11) Urban, M.; Hubacˇ, I.; Kello¨, V.; Noga, J. J. Chem. Phys. 1980, 72, 3378. (12) Szalay, P. G. Int. J. Quantum Chem. 1995, 55, 151. (13) Salter, E. A.; Trucks, G. W.; Bartlett, R. J. J. Chem. Phys. 1989, 90, 1752. (14) Szalay, P. G.; Nooijen, M.; Bartlett, R. J. J. Chem. Phys. 1995, 103, 281 and references therein. (15) Jørgensen, P.; Helgaker, T. J. Chem. Phys. 1988, 89, 1560. Helgaker, T.; Jørgensen, P. Theor. Chim. Acta 1989, 75, 111. (16) Lo¨wdin, P.-O. AdV. Chem. Phys. 1959, 2, 270. Helgaker, T. U.; Almlo¨f, J. Int. J. Quantum Chem. 1984, 26, 275. (17) C ˇ ´ızˇek, J. J. Chem. Phys. 1966, 45, 4256; AdV. Chem. Phys. 1969, 14, 35.
J. Phys. Chem., Vol. 100, No. 15, 1996 6297 (18) Bartlett, R. J.; Shavitt, I. Chem. Phys. Lett. 1977, 50, 190. (19) Koch, S.; Kutzelnigg, W. Theor. Chim. Acta 1981, 59, 387. (20) Ruttink, P. J. A.; van Lenthe, J. H.; Zwaans, R.; Groenenboom, G. C. J. Chem. Phys. 1991, 94, 7212. (21) Pulay, P.; Hamilton, T. P. J. Chem. Phys. 1988, 88, 4926. (22) Bofill, J. M.; Pulay, P. J. Chem. Phys. 1989, 90, 3637. (23) Kutzelnigg, W. In Modern theoretical chemistry; Schaefer, H. F., Ed.; Plenum Press: New York, 1977; p 129. (24) Laidig, W. D.; Bartlett, R. J. Chem. Phys. Lett. 1984, 104, 424. Laidig, W. D.; Saxe, P.; Bartlett, R. J. J. Chem. Phys. 1987, 86, 887. (25) Paldus, J. In New horizons in quantum chemistry; Lo¨vdin, P-O., Pullman, B., Eds.; Reidel: Dodrecht, 1983; p 31. (26) Hoffmann, M. R.; Simons, J. J. Chem. Phys. 1989, 90, 3671. (27) Cave, R. J.; Davidson, E. R. J. Chem. Phys. 1988, 88, 5770. (28) Cave, R. J.; Davidson, E. R. J. Chem. Phys. 1988, 89, 6798. (29) Tanaka, K.; Sakai, T.; Terashima, H. Theor. Chim. Acta 1989, 76, 213. Sakai, T.; Tanaka, K. Theor. Chim. Acta 1993, 83, 451. (30) Daudey, J-P.; Heully, J-L.; Malrieu, J-P. J. Chem. Phys. 1993, 99, 1240. (31) Malrieu, J-P.; Daudey, J-P.; Caballol, R. J. Chem. Phys. 1994, 101, 8908. (32) Fink, R.; Staemmler, V. Theor. Chim. Acta 1993, 87, 129. (33) Shepard, R.; Shavitt, I.; Pitzer, R. M.; Comeau, D. C.; Pepper, M.; Lischka, H.; Szalay, P. G.; Ahlrichs, R.; Brown, F. B.; Zhao, J.-G. Int. J. Quantum Chem. 1988, S22, 149. (34) Shepard, R.; Lischka, H.; Szalay, P. G.; Kovar, T.; Ernzerhof, M. J. Chem. Phys. 1992, 96, 2085 and references therein. (35) Pulay, P. TX90; Fayetteville, AR, 1990. Pulay, P. Theor. Chim. Acta 1979, 50, 299. (36) Bauschlicher, C. W.; Taylor, P. R. J. Chem. Phys. 1986, 85, 2779. (37) Olsen, J.; Jørgensen, P.; Koch, H.; Balkova, A.; Bartlett, R. J. Full Configuration-Interaction and State of the Art Correlation Calculations on Water in a Valence Double Zeta Basis With Polarization Functions. To be published. (38) Harrison, R. J.; Handy, N. C. Chem. Phys. Lett. 1983, 98, 97. (39) Sosa, C.; Noga, J.; Bartlett, R. J. J. Chem. Phys. 1988, 88, 5974. (40) Widmark, P. O.; Malquist, P. A.; Roos, B. O. Theor. Chim. Acta 1990, 77, 291. (41) Nebot-Gil, I.; Sa´nchez-Marı´n, J.; Malrieu, J-P.; Heully, J-L.; Maynau, D. J. Chem. Phys. 1995, 103, 2576. (42) Purvis, G. D.; Shepard, R.; Brown, F. B.; Bartlett, R. J. Int. J. Quantum Chem. 1983, 23, 835. (43) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1918. (44) Watts, J. D.; Stanton, J. F.; Bartlett, R. J. Chem. Phys. Lett. 1991, 178, 471. (45) Del Bene, J. E.; Stahlberg, E. A.; Shavitt, I. Int. J. Quantum Chem. 1990, S24, 455. (46) Dunning, T. H., Jr. J. Chem. Phys. 1970, 53, 2823.
JP952840J
