Benchmark Quasi-Variational Coupled Cluster Calculations of Multiple

Article pubs.acs.org/JCTC

Benchmark Quasi-Variational Coupled Cluster Calculations of Multiple Bond Breaking James B. Robinson and Peter J. Knowles* School of Chemistry, Cardiff University, Main Building, Park Place, Cardiff CF10 3AT, United Kingdom ABSTRACT: We present further evidence that closed-shell single-reference coupled cluster theory can be used as a reliable predictive tool even when multiple chemical bonds are broken, provided that a near-variational, rather than a projective, ansatz is used. Building on the Optimized-orbital Quasi-Variational Coupled Cluster Doubles (OQVCCD) method by adding the standard perturbative treatment of triple excitations, the OQVCCD(T) method provides outstanding accuracy for the dissociation of multiply bonded molecules and other problems involving strong nondynamic correlation of the electrons. We find that in the case of singly bonded molecules, OQVCCD and OQVCCD(T) perform similarly to the equivalent Brueckner Coupled Cluster Doubles approaches, BCCD and BCCD(T). However, when multiple bonds are broken, such as in acetylene and dicarbon, OQVCCD(T) is capable of predicting both qualitatively and quantitatively accurate potential energy curves, unlike the standard methods based on traditional coupled cluster theory, and for approximately the same computational cost.

1. INTRODUCTION Ab initio electronic structure calculations are those that attempt to find an approximate solution of the Born−Oppenheimer electronic Schrödinger equation. In single-reference calculations, the Hartree−Fock approximation is first made, and the optimum single-determinantal reference wavefunction is constructed. Since the true electronic wavefunction can be expanded exactly in the basis of the determinants that can be generated by the replacement of one or more of the orbitals occupied in the reference wavefunction with unoccupied orbitals, single-reference methods attempt to capture a representative subset of these terms by truncation to, for example, at most 2-fold excitations of the reference determinant. Single-reference methods thus have the advantage of being somewhat black-box. The widely used and extremely successful Traditional Coupled Cluster1−8 (TCC) method belongs to this category. Unfortunately, situations exist in which the Hartree−Fock approximation breaks down, so that the single-determinantal reference wavefunction provides a poor underlying description of the system. This is known as the regime of nondynamic (or static) correlation. The construction of a reference wavefunction consisting of multiple determinants has the potential to resolve this problem, and examples in this category are the multireference CI9,10 (MRCI) method and related formulations that seek to approximately correct for size-extensivity errors.11−13 However, the determinants to be included in the reference wavefunction must be chosen well, and these methods are thus difficult to use in a black-box fashion, especially on large molecules, are often highly expensive in terms of computational effort, and encounter problems due to the lack of rigorous extensivity of the energy. Other, more novel approaches to the treatment of static correlation have been proposed, such as the active-space CC methods of Head-Gordon et al.,14−16 and the spin-flip17,18 and double-ionization-potential19,20 EOM methods, but an allpurpose method has yet to emerge. Thus, there exist systems for which practical single-reference methods such as TCC © 2012 American Chemical Society

truncated to the singles and doubles level (CCSD) yield inadequate descriptions of the electronic structure, and for which multireference methods are impractical. A new method suitable for the treatment of problems in this niche would be highly desirable. Variational Coupled Cluster21 (VCC), a single-reference method, has been demonstrated22−29 to be, for a given truncation of the cluster operator, a significantly more robust electronic structure ansatz than the equivalent TCC method. This is because, unlike in VCC, the non-Hermitian nature of ̂ ̂ the similarity-transformed Hamiltonian, H̅ = e−TĤ eT, and the subsequent projective determination of the amplitude equations effectively eliminates the methodological property that calculated energies are upper bounds on the exact ground-state Schrödinger energy eigenvalue from the TCC approach. The widespread use of TCC, however, can be attributed to the Campbell−Baker−Hausdorff expansion30 of this effective Hamiltonian, which allows TCC calculations to be performed in polynomial time. In contrast, there is no such simplification for VCC, which, with factorial time complexity, is prohibitively expensive for all but extremely small systems. We have recently put forward a family of new singlereference quantum-chemical methods that function by constructing infinite-order approximations to VCC restricted to double excitations (VCCD).31−33 We have demonstrated that the most advanced of these methods, Quasi-Variational Coupled Cluster (QVCC), is capable of describing nondynamic electron correlation with a quality closer to VCC than TCC. However, to achieve quantitative chemical accuracy, it is wellknown that one must go beyond the double-excitation model by including the effects of triple excitations. For problems dominated by Hartree−Fock, most of the effect can be recovered by a simple perturbative treatment of the triples, and this is true whether one uses TCCSD or VCCSD, since in Received: April 11, 2012 Published: July 2, 2012 2653

dx.doi.org/10.1021/ct300416b | J. Chem. Theory Comput. 2012, 8, 2653−2660

Journal of Chemical Theory and Computation

Article

this regime they are nearly the same. However, when strong correlation is present, such as when multiple chemical bonds are broken, we have discovered some surprising evidence that the simple perturbative correction for triple excitations can remain excellent, provided that the leading-order doubleexcitation correlations are treated sufficiently well first.34 For example, the Optimized-orbital Quasi-Variational Coupled Cluster Doubles (OQVCCD) method, when combined with a perturbative treatment of triple excitations (OQVCCD(T)), is capable of achieving a physically correct and quantitatively accurate description of the potential energy curve of dinitrogen, N2,34 a system for which the current generation of practical single-reference methods catastrophically fails. In this article, we present the results of further benchmark calculations with the OQVCCD(T) method on systems exhibiting strong nondynamic correlation, giving more substantial evidence that the method is generally useful.

5. The series expansion should agree with a subset of the VCCD 6 (T4) contributions that can be computed with at most 6 (N2v4) effort. An energy functional, to be minimized variationally, that satisfies these aims can be written E0 = ⟨0|Ĥ |0⟩ † ̂ ∼ |̂ 0⟩ + ⟨0|T ̂ (Ĥ − E0)T |̂ 0⟩ E = E0 + 2⟨0|HT

(2)

∼ The T ̂ and T ̂ operators are transformations of the original cluster operator T̂ 2, with amplitudes defined by ij −1/2 c ij −1/2 i kj Tab ̅ = (1 − τab)(A )a Tcb + (1 − τij)(B )k Tab

−

1 −1/2 ij kl 1 (C )kl Tab − (1 − τij)(1 − τab)(D−1/2 )icak Tcbkj 2 2 (3)

2. REVIEW OF THEORY We here summarize the elements of the OQVCCD(T) theory, all of which have been presented and discussed in our earlier papers.31−34 The method is motivated by the desire to approximate VCC as well as possible, but without introducing a computational resource requirement that significantly exceeds that of TCC. At the level of a cluster operator that contains only double excitations, the projective coupled cluster (TCCD) energy agrees with VCCD up to and including terms of third order in the cluster operator. The leading-order difference between VCCD and TCCD35 is the two-electron part of † 1 ̂ 2̂ 2|0⟩L ⟨0|(T2̂ )2 HT 4

1 ij ac 1 ij ab TbcTij Cklij = δklij + Tab Tkl 2 2 1 ik ab Bji = δji + Tab T jk Dajib = δajib + TacikT jkbc 2

Aba = δba +

(4)

∼ where τpq permutes the labels p and q in what follows. T ̂ is defined similarly, but with the powers of the matrices A, B, C, and D changed from −1/2 to −1. This defines the QVCCD approximation. The rationale for the use of the four additive transformations, together with the choice of the various numerical factors, is discussed in detail in ref 33 and provides for the first four aims listed above. The binomial expansion of eq 3 also provides the oneelectron part of eq 1 exactly, as well as all of its two electron parts except for the terms where the Hamiltonian is connected to all four cluster operators, which are all absent. Although those missing terms (and many others at higher order) are individually nonzero for a two-electron (or two-virtual) subsystem, they sum to zero in that case. Thus, this approximation satisfies the fifth aim listed above and provides a pragmatic approximation to VCCD. The computational effort is similar to that needed for CCD, except for the calculation of the matrix powers, especially that of D, of dimensions Nv × Nv. The overall scaling is, however, still 6 (N3v3) plus the same 6 (N2v4) term that occurs in CCD. Another perspective on the role of the matrix transformations of the cluster amplitudes is to introduce partial local normalization. In the limit of two electrons, for example, these transformations correctly reintroduce division by the variational CID denominator into the above energy functional. Our functional may therefore be viewed as an extension of the coupled pair functional,36 and related approaches,37−39 with the matrix transformations replacing the explicit partial local denominators. It is additionally related to the Coupled Electron Pair Approximation (CEPA),35,40−44 modern developments of which have even taken a VCC viewpoint in order to facilitate their construction.45 One can show that the Linked Pair Functional (LPF)31 contains in full the two 6 (T3) diagrams for which the EPV contributions constitute the basis of the CEPA method.35 QVCCD goes beyond the LPF by including completely all four linked ⟨T̂ †Ĥ T̂ 2⟩ diagrams, as well as higher-order contributions. Our numerical results accordingly demonstrate superiority to CEPA, particularly as covalent bonds are stretched and the magnitude of T̂ increases.

(1)

where |0⟩ is the reference Slater determinant. Unfortunately, this contribution to the energy, which is expected to be important when VCC and TCC differ significantly, is difficult to compute in full. In particular, one of its parts contains the four-external orbital part of the Hamiltonian connected to each of the four T̂ operators, and its evaluation necessarily requires a computational effort that scales as 6 (v6), where v is the number of virtual orbitals. This is significantly more expensive than the effort needed for CCD, where, although the overall scaling is also sixth order in N, the number of electrons, it is at most 6 (N2v4), that is, fourth order in the virtual space. The QVCCD approximation is constructed with the following aims. 1. Exact results should be obtained for two electrons, or in the case where there are just two virtual orbitals. 2. A closed-form expression that exactly sums classes of important contributions would be preferable to a truncated series, since it is anticipated that the approximation will be deployed in cases where the correlation is strong, and therefore the series expansion of the VCCD energy will not necessarily be rapidly convergent. 3. The approximation should be a fully linked tensor contraction, in order to satisfy the requirement of extensivity and to give the desirable property of invariance with respect to rotations in the occupied and virtual orbital spaces. 4. The series expansion of the approximation should agree with that of VCCD through 6 (T3). 2654



Article

breaking of single bonds.64 We have additionally obtained CEPA results using Orca.65 To begin, we investigate two simple single bond breaking examples. We give calculated potential energy curves for BH with the cc-pVQZ66 basis in Figure 1 and for HF with the aug-

In order to account for the effects of single excitations, we minimize the functional also with respect to the orbitals,46 which defines the Optimized-orbital Quasi-Variational Coupled Cluster Doubles (OQVCCD) method, although it is also possible to make use of a Brueckner condition47−49 for this purpose, which defines the Brueckner Quasi-Variational Coupled Cluster Doubles (BQVCCD) method. The OQVCCD method is already correct to third order in Møller−Plesset perturbation theory and correctly constructs the terms involving single and double excitations at fourth order. However, connected triple excitations first enter at fourth order through the second-order wavefunction, and these terms are, of course, omitted by OQVCCD. This is the same position as CCSD. The [T] correction50 to CCSD is a minimal correction for this omission that is constructed from the converged CCSD singles and doubles cluster amplitudes, along with the one- and two-electron integrals, such that the CCSD[T] method is correct to fourth-order in Møller−Plesset theory. The CCSD(T)51 method also includes some additional terms at fifth-order and higher that have been justified in different ways,51−53 and further corrective terms have also been proposed.54 Mutual cancellation between the VCCD terms at fourth-order ensures that the standard [T] correction50 to CCSD is valid for OQVCCD, making it correct to fourth-order in Møller−Plesset theory also. Since, in OQVCCD, the effects of single excitations are introduced through orbital optimization, and the singles cluster amplitudes are thus formally zero, the [T] and (T) corrections are identical, since they differ only in terms that contain single excitation amplitudes. This defines the OQVCCD(T) method, which we investigate further in this contribution. We note that, since it is our intention to explore systems for which static correlation is strong, it might be preferable to design a correction for OQVCCD in the spirit of (2),55−57 rather than making use of the standard (T), in order to more effectively deal with the effects of connected triples (and quadruples) in the nondynamic regime by further decoupling the perturbative correction from the Hartree−Fock approximation. However, we defer this to later work and use only (T) in this initial exploration; it will be demonstrated shortly that the (T) correction itself can be greatly improved by a more robust description of static correlation phenomena at the doubles level, as was the case when we applied OQVCCD(T) to dinitrogen, N2.34 Here, we apply the OQVCCD(T) method to further benchmark problems, including several for which we have demonstrated OQVCCD to be a more faithful approximation to VCCSD than CCSD.33

Figure 1. Calculated potential energy curves for BH with the cc-pVQZ basis set.

Figure 2. Calculated potential energy curves for HF with the aug-ccpVQZ basis set.

cc-pVQZ67 basis in Figure 2. In both examples, the CCSD(T) method becomes poor as the bond is stretched. This is a wellknown problem associated with the (T) correction, since it becomes singular when the highest occupied and lowest unoccupied molecular orbitals become degenerate. However, while CCSD and OQVCCD perform similarly, the OQVCCD(T) method fares significantly better than CCSD(T). While this is promising, comparison with BCCD and BCCD(T) shows that this behavior may be a result of the different orbitals in use; CCSD(T) uses Hartree−Fock orbitals, whereas BCCD(T) uses Brueckner orbitals and OQVCCD(T) uses variationally optimal orbitals. Equivalently, these differences are a result of the different treatment of single excitations, due to the Thouless theorem,68 which states that any two singledeterminantal wavefunctions, |Φ⟩ and |Φ′⟩, may be related by ̂ |Φ′⟩ = eT1|Φ⟩. This behavior has been noted previously, for example by Nooijen and LeRoy, who found that the use of Brueckner orbitals improved the triples corrections substantially in HF, BeO, CN, and BN.45 In the calculations that

3. RESULTS AND DISCUSSION We have performed most calculations in this article with the Molpro59,60 quantum chemistry software package. For each molecule studied, we compare one-dimensional cuts of the potential energy surface obtained with various single-reference coupled-cluster methodologies with those obtained from internally contracted multireference configuration interaction9,10 (MRCI) calculations. These reference calculations used complete active space reference wavefunctions where the active space consists of the atomic valence orbitals, and the energy was corrected using the approximate extensivity correction of Davidson and Langhoff61 (MRCI+Q). Using the GAMESS62 package, we also examine the CR-CC(2,3) method,63 which has yielded some impressive results for the 2655



Article

follow, we use the BCCD and BCCD(T) methods to identify those systems for which the choice of orbitals affects the triples corrections more than the differences in the doubles-only theories. It should be noted, however, that the Brueckner orbitals are not always close to the variationally optimal orbitals,46 and we have checked, for additional clarity, that Brueckner Quasi-Variational Coupled Cluster Doubles (BQVCCD) performs in agreement with OQVCCD. We examine next the spectroscopic constants for a selection of diatomic molecules, given in Table 1. The CCSD and Table 1. Comparison of Equilibrium Bond Lengths and Spectroscopic Constants for Some Diatomic Moleculesa system

method

Re/Å

ωe/cm−1

ωexe/cm−1

HF

CCSD BCCD BQVCCD OQVCCD CCSD(T) BCCD(T) BQVCCD(T) OQVCCD(T) empirical58 CCSD BCCD BQVCCD OQVCCD CCSD(T) BCCD(T) BQVCCD(T) OQVCCD(T) empirical CCSD BCCD BQVCCD OQVCCD CCSD(T) BCCD(T) BQVCCD(T) OQVCCD(T) empirical

0.914 0.914 0.914 0.914 0.917 0.917 0.917 0.917 0.917 1.388 1.386 1.385 1.385 1.410 1.410 1.407 1.407 1.412 1.092 1.091 1.090 1.090 1.099 1.099 1.097 1.098 1.098

4198.7 4197.9 4200.9 4201.3 4148.4 4141.3 4146.5 4150.8 4138.3 1025.5 1030.3 1034.3 1034.6 929.2 933.3 942.1 942.9 916.6 2445.3 2456.2 2464.0 2461.0 2364.9 2370.5 2384.5 2382.0 2358.6

93.9 87.6 87.9 89.5 95.0 87.7 88.5 91.9 89.9 8.7 8.7 8.7 8.6 11.4 11.3 11.2 11.2 11.2 12.8 12.6 12.5 12.5 13.8 13.7 13.4 13.5 14.3

F2

N2

Figure 3. Calculated potential energy curves for D2h H2SiSiH2 as a function of the SiSi bond length with the cc-pV(D+d)Z basis set.

predicted by CCSD. The shape of the OQVCCD curve is, in fact, reminiscent of the shape of the VCCSD curve in examples benchmarked previously.33 Both CCSD(T) and BCCD(T) behave nonvariationally here as the double bond is stretched, whereas OQVCCD(T) remains above, and in good agreement with, the MRCI+Q curve throughout. Next, we examine the case of the symmetric stretching of the N−H bonds in ammonia, NH3. We give calculated potential energy curves for this system in Figure 4. This example involves

Basis set: cc-pV5Z, with correlation energy x−3-extrapolated using ccpVQZ and cc-pV5Z.

a

Figure 4. Calculated potential energy curves for the symmetric stretching of NH3 with the cc-pVTZ basis set.

OQVCCD results are of similar quality, as are the CCSD(T) and OQVCCD(T) results. This is to be expected due to the similarity of the potential energy curves around equilibria, evident in later figures. We point out that comparing the spectroscopic constants obtained from the BQVCCD, BCCD, BQVCCD(T), and BCCD(T) methods reveals deficiencies that arise in our current QVCCD approximation to VCCD, manifesting as incomplete recovery of dynamic correlation energy relative to CCD, leading to a slightly poorer description of the equilibrium region of potential energy surfaces. These defects are, however, quite small. Our first example that involves strong static correlation is the case of the symmetric stretching of a double bond, for which we take H2SiSiH2 constrained to a planar D2h geometry, and with the Si−H bond length and bond angle optimized at the CCSD level of theory for each Si−Si distance. In this system, the curve predicted by OQVCCD, shown with the cc-pV(D +d)Z basis in Figure 3, is qualitatively different from that

the simultaneous breaking of three single bonds, and the CCSD and BCCD methods struggle as these bonds are stretched. The BCCD method undergoes a nonvariational collapse to energies below the MRCI+Q curve, and the CCSD curve becomes unstable and begins to increase in energy too sharply from 2.2 Å. Both curves are clearly wrong in comparison to the MRCI +Q curve. However, the OQVCCD curve is smooth and continues fairly parallel to MRCI+Q throughout. The differences are much more obvious when the (T) correction is added to each method; both CCSD(T) and BCCD(T) now diverge unphysically, whereas OQVCCD(T) is in excellent agreement with MRCI+Q, representing a significant improvement over the methods based on TCC. While CR-CC(2,3) does not demonstrate a failure as pronounced as in CCSD(T), it exhibits instabilities similar to the CCSD method on which it is based. We have additionally examined the CEPA-2 method35,41 for this example, and our results show it to diverge from MRCI 2656



Article

+Q and become catastrophically unphysical at bond lengths even shorter than for which the breakdown of CCSD(T) occurs. The symmetric stretching of the C−H bonds in ethene, in which four bonds are now simultaneously broken, represents a similar, but more extreme, test case. In this example, the CCSD and BCCD methods both become unphysical and qualitatively incorrect in comparison to MRCI+Q from 2.2 Å, even before the corrections for triples are added. This is illustrated in Figure 5. However, OQVCCD possesses no unphysical maximum, and

Figure 6. Calculated cc-pVQZ energies of the H4 model system as a function of θ for R = 1.738 Å.

closely. We now investigate the effects of adding the (T) correction. In CCSD(T), the problematic behavior at θ = 90° is made worse: the energy gradient becomes steeper, and the curve is shifted further away from MRCI+Q. However, in OQVCCD(T), the shift to lower energies causes the peak to become closer to that of MRCI+Q, and, additionally, the slope at θ = 90° is significantly reduced relative to both CCSD(T) and OQVCCD, so that the (T) correction improves the situation. With the “wrong” reference, all single-reference curves are incorrect, and although the OQVCCD method remains an upper bound, this system is more sensibly treated by ensuring that the reference determinant in use is always the optimum determinant, which exploits the symmetry around θ = 90°. A plot of the forces, calculated by finite difference differentiation, is given in Figure 7 and shows the significant

Figure 5. Calculated potential energy curves for the symmetric stretching of the C−H bonds in ethene with the cc-pVDZ basis set.

the curve continues as one would expect from physical intuition. When the triples corrections are added, CCSD(T) and BCCD(T) diverge even more rapidly, and CR-CC(2,3) performs only slightly better. The contrast between these methods and OQVCCD(T), which is qualitatively correct and in good quantitative agreement with MRCI+Q throughout, is quite extraordinary in this example. Our next example is the previously studied H4 model system, in which four hydrogen atoms are arranged in a rectangle specified by the parameters R and θ.26 These control the distance of each H atom from the center of the rectangular arrangement and the angle subtended at the center of mass by radii to two neighboring H nuclei, respectively. The minimalenergy single-determinantal reference wavefunction for this system differs depending on whether θ belongs to the interval [0°,90°) or (90°,180°]. There is thus strong multireference character around θ = 90°, which corresponds to the square geometry. The symmetry of the system dictates that the exact potential energy surface is symmetrical about θ = 90°, whereas approximate methods based on a single-determinantal reference will not be, or, if for all geometries the lowest energy determinant is used as the reference, will show unphysical cusps at the square geometry. We have previously shown that when θ is varied with R fixed at 1.738 Å, the OQVCCD method faithfully reproduces the shape of the VCCSD potential energy curve33 that was first obtained by Van Voorhis and HeadGordon.26 Figure 6 shows the performance of various coupledcluster methods based on the locally optimum reference wavefunction. Their asymmetric performances with the single determinant optimum for small θ, continuing beyond θ = 90° with this “wrong” reference, are also indicated on this plot by the fainter points. While CCSD and OQVCCD both predict unphysical nonzero slope at the square geometry, the OQVCCD curve matches that of MRCI+Q much more

Figure 7. Calculated cc-pVQZ forces of the H4 model system as a function of θ for R = 1.738 Å.

mitigation of the discontinuity by OQVCCD(T) relative to either CCSD(T) or BCCD(T). This is especially true in the range θ = 80−85°, where the CCSD and CCSD(T) forces deviate strongly from MRCI+Q, but OQVCCD and OQVCCD(T) remain accurate. The acetylene molecule, C2H2, possesses an electronic structure that is analogous to the N2 molecule and thus represents a similarly extreme test of molecular electronic structure methods when the triple bond is stretched. The results of our calculations on this system are shown in Figure 8, for a fixed C−H bond length of 1.06 Å. From approximately 2657



Article

OQVCCD to predict a larger dissociation energy than CCSD. Again, this behavior of OQVCCD can be seen to reflect the underlying divergence between CCSD and VCCSD that occurs if the system is examined in a minimal basis.33 When the triples corrections are added, CCSD(T) and OQVCCD(T) are similar around equilibrium, but CCSD(T) becomes unphysical from approximately 2.3 Å. On the other hand, OQVCCD(T) remains physically correct throughout. It is noteworthy that, just as in our previous examples, the CCSD(T) and BCCD(T) methods turn over precisely when CCSD and BCCD begin to diverge from OQVCCD, from which we infer that it is the divergence of TCC from VCC at this point that primarily instigates the failure of the triples correction for these methods, and that the excellent performance of OQVCCD(T) here is a direct result of the more robust treatment of static correlation phenomena inherited by our methodology of approximating a parent theory, VCCD, which possesses true upper bound character. The curves presented have been plotted relative to the correct dissociation limits in order to demonstrate that there is no significant cancellation of errors between the asymptotic absolute energies and equivalent calculations on the atoms for either CCSD or OQVCCD, but there is for OQVCCD(T). To further support these results, we have included a plot of the forces, calculated by finite difference differentiation, in Figure 10, which shows the OQVCCD(T)

Figure 8. Calculated potential energy curves for the stretching of the carbon−carbon triple bond in acetylene with the aug-cc-pVTZ basis set.67

2.2 Å, none of the methods based on TCC is correct, each predicting an unphysical maximum in the potential energy curve followed by a nonvariational collapse to energies below MRCI+Q. The effect of the (T) correction on each of these methods is to push the energy even lower, causing the problem to become magnified. The CEPA-2 method also diverges to unphysically low energies on stretching, and the onset of this behavior occurs at even shorter bond lengths (from around 1.6 Å). In contrast, the OQVCCD method does not appear to degrade significantly in quality as the bond is stretched, predicting a potential energy curve with the characteristic VCCSD shape.33 The additional (T) perturbative correction of the energy results in a predicted curve that is in outstanding agreement with MRCI+Q. The C2 1Σ+g state is also a difficult test of single-reference correlation methods. We give calculated potential energy curves with the aug-cc-pVQZ basis in Figure 9. The energies for each

Figure 10. Calculated force curves for C2 with the aug-cc-pVQZ basis set.

curve to be almost flat, predicting a correct dissociation of the molecule, well beyond the region in which CCSD(T) and BCCD(T) become exceedingly inaccurate. Furthermore, we have included a plot of the square norms of the cluster amplitudes for the methods CCSD, BCCD, and OQVCCD in Figure 11. The OQVCCD square norm increases slowly with a decreasing gradient from around 2.4 Å, reaching a value of 2.40 at 1000 Å, suggesting that only small changes will occur at even longer bond lengths. In contrast, the CCSD and BCCD square norms approach their respective limits more slowly than OQVCCD. Asymptotically, the Hartree−Fock HOMO− LUMO gap decays inversely with distance, but this does not mean that at very long bond lengths the (T) correction will diverge linearly with distance. Single excitations from bonding orbitals to their corresponding antibonding orbitals have odd symmetry with respect to inversion in a center of symmetry, and therefore there are no symmetry-allowed triple excitations for which the Fock excitation energy tends to zero asymptotically. Our calculations demonstrate that where the covalent

Figure 9. Calculated energies for C2 relative to the atomic dissociation limits with the aug-cc-pVQZ basis set.

method are presented relative to twice the energy from equivalent atomic calculations; in the case of OQVCCD and OQVCCD(T), for technical reasons we are not yet able to perform open-shell calculations, and in these cases (and for CCSD and CCSD(T)), atomic spin-restricted CCSD and CCSD(T)69,70 values have been used. The OQVCCD method diverges from the CCSD method at long bond lengths; unlike CCSD, the curve does not plateau around 2.4 Å, leading 2658



Article

spin-contamination effects arising from the use of the UHF reference.72,73 Thus, although UHF-CCSD outperforms OQVCCD in the dissociation regime, as is to be expected, OQVCCD is more appropriate at intermediate radii. Furthermore, the UHF-VCCSD results suggest that the upper bound property potentially mitigates spin-contamination effects, and this has also been noted in previous work.28 A hypothetical UHF-OQVCCD theory may do the same.

4. CONCLUDING REMARKS We have provided further evidence that combining the OQVCCD approximation with the standard Møller−Plessetbased (T) perturbative introduction of connected triple excitations gives a powerful and versatile electronic structure ansatz. When a single-determinantal reference wavefunction is an adequate qualitative description of the molecular electronic structure, we find little difference between the BCCD(T) and OQVCCD(T) methods. Furthermore, from the examples of BH and HF, we have shown that some improvement in the (T) correction emerges when Brueckner or variationally optimal orbitals are used, although this has been noted previously by other authors.45 However, when the reference wavefunction becomes poor and static correlation effects become strong, not only is OQVCCD(T) significantly more robust than CCSD(T), but in extreme examples it is able to achieve a physically correct description of the energetics of bond breaking, and these improvements appear not to be related to the different orbitals in use. Thus, although the (T) correction can still be expected to break down when the HOMO−LUMO gap narrows sufficiently and the perturbative correction approaches singularity, unlike CCSD, the pseudovariational ansatz of OQVCCD theory does not fail entirely to capture the essential physics of nondynamic electron correlation, and the range of systems for which the remaining dynamic correlation effects can be legitimately included perturbatively appears greatly extended. In particular, it is remarkable that a quantitatively accurate description of the complete potential energy curve, from the repulsive domain through the equilibrium geometry to the dissociation limit, of molecules with electronic structures as complicated as acetylene and dicarbon can be achieved by a practical single-reference method.

Figure 11. Amplitude square norms for C2 with the aug-cc-pVQZ basis set.

bond has broken, OQVCCD(T) can give accurate results and tends to a constant energy close to the true dissociation limit. As a final test of the Quasi-Variational Coupled Cluster method, which appears particularly robust to the breakdown of the Restricted Hartree−Fock (RHF) approximation, we compare it to results that can be obtained by Coupled Cluster methods making use of an Unrestricted Hartree−Fock (UHF) reference wavefunction. We take the example of a model system consisting of six hydrogen atoms arranged to be equally spaced on the circumference of a circle. Errors relative to FCI energies are presented as a function of the radius of this circle in Figure 12. The UHF method and post-Hartree−Fock

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: KnowlesPJ@Cardiff.ac.uk. Figure 12. Errors relative to FCI for calculated energies for six H atoms equally spaced on the circumference of a circle as a function of the radius of the circle, with the STO-3G basis set.

Notes

methods based on it, such as UHF-CCSD, correctly describe the electron localization that occurs for large radii as the atoms become isolated, so that the UHF-CCSD method agrees with FCI at sufficiently long bond lengths. As is expected, although the OQVCCD method remains above FCI and physically correct throughout, it deteriorates in accuracy at longer bond lengths, overestimating the energy at dissociation due to the strong ionic contamination of the RHF reference.71 It still performs significantly better than either CCSD or BCCD, however, which are qualitatively wrong from 1.25 Å. Of particular interest, however, is the region between 1.0 and 2.0 Å, in which UHF-CCSD predicts a significantly larger error than the RHF-OQVCCD method. This can be attributed to

(1) Coester, F. Nucl. Phys. 1958, 7, 421. (2) Coester, F.; Kümmel, H. Nucl. Phys. 1960, 17, 477. (3) Č ížek, J. J. Chem. Phys. 1966, 45, 4256. (4) Č ížek, J. Adv. Chem. Phys. 1969, 14, 35. (5) Č ížek, J.; Paldus, J. Int. J. Quantum Chem. 1971, 5, 359. (6) Barlett, R. J.; Purvis, G. D., III. Int. J. Quantum Chem. 1978, 14, 561. (7) Bartlett, R. J. Annu. Rev. Phys. Chem. 1981, 32, 359. (8) Purvis, G. D., III; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910. (9) Werner, H.-J.; Knowles, P. J. J. Chem. Phys. 1988, 89, 5803−5814. (10) Knowles, P. J.; Werner, H.-J. Chem. Phys. Lett. 1988, 145, 514− 522. (11) Gdanitz, R. J.; Ahlrichs, R. Chem. Phys. Lett. 1988, 143, 413. (12) Szalay, P. G.; Bartlett, R. J. Chem. Phys. Lett. 1993, 214, 481.

The authors declare no competing financial interest.

■

2659

REFERENCES



Article

(13) Werner, H.-J.; Knowles, P. J. Theor. Chem. Acc. 1991, 78, 175− 187. (14) Krylov, A.; Sherrill, C.; Byrd, E.; Head-Gordon, M. J. Chem. Phys. 1998, 109, 10669. (15) Parkhill, J.; Head-Gordon, M. J. Chem. Phys. 2010, 133, 124102. (16) Parkhill, J.; Head-Gordon, M. J. Chem. Phys. 2010, 133, 024103. (17) Krylov, A. Chem. Phys. Lett. 2001, 338, 375. (18) Krylov, A. Chem. Phys. Lett. 2001, 350, 522. (19) Nooijen, M.; Bartlett, R. J. Chem. Phys. 1997, 107, 6812. (20) Sattelmeyer, K. W.; Schaefer, H. F., III; Stanton, J. F. Chem. Phys. Lett. 2003, 378, 42. (21) Bartlett, R. J.; Noga, J. Chem. Phys. Lett. 1988, 150, 29. (22) Fan, P.-D.; Kowalski, K.; Piecuch, P. Mol. Phys. 2005, 103, 2191−2213. (23) Fan, P.-D.; Piecuch, P. Adv. Quantum Chem. 2006, 51, 1−57. (24) Piecuch, P.; Pimienta, I. S. O.; Fan, P.-D.; Kowalski, K. ACS Symp. Ser. 2007, 958, 37−73. (25) Piecuch, P.; Kowalski, K.; Fan, P.-D.; Pimienta, I. S. O. Prog. Theor. Chem. Phys. 2010, 12, 119−206. (26) Van Voorhis, T.; Head-Gordon, M. J. Chem. Phys. 2000, 113, 8873. (27) Fan, P.-D.; Kowalski, K.; Piecuch, P. Mol. Phys. 2005, 103, 2191. (28) Cooper, B.; Knowles, P. J. J. Chem. Phys. 2010, 133, 234102. (29) Evangelista, F. A. J. Chem. Phys. 2011, 134, 224102. (30) Merzbacher, E. Quantum Mechanics, 2nd ed.; Wiley: New York, 1970. (31) Knowles, P. J.; Cooper, B. J. Chem. Phys. 2010, 133, 224106. (32) Robinson, J. B.; Knowles, P. J. J. Chem. Phys. 2011, 135, 044113. (33) Robinson, J. B.; Knowles, P. J. J. Chem. Phys. 2012, 136, 054114. (34) Robinson, J. B.; Knowles, P. J. Phys. Chem. Chem. Phys. 2012, 14, 6729. (35) Kutzelnigg, W. In Recent Progress in Coupled Cluster Methods; Č ársky, P., Paldus, J., Pittner, J., Eds.; Springer: New York, 2010; pp 1−61. (36) Ahlrichs, R.; Scharf, P.; Ehrhardt, C. J. Chem. Phys. 1985, 82, 890. (37) Kollmar, C.; Neese, F. J. Chem. Phys. 2011, 135, 084102. (38) Kollmar, C. J. Chem. Phys. 2006, 125, 084108. (39) Deprince, A. E.; Mazziotti, D. A. Phys. Rev. A 2007, 76, 042501. (40) Ahlrichs, R. J. Chem. Phys. 1968, 48, 1819−1832. (41) Meyer, W. J. Chem. Phys. 1973, 58, 1017. (42) Kollmar, C.; Neese, F. Mol. Phys. 2010, 108, 2449−2458. (43) Kutzelnigg, W.; von Herigonte, P. Electron Correlation at the Dawn of the 21st Century; Academic Press: New York, 1999; Vol. 36; Advances in Quantum Chemistry, pp 185−229. (44) Koch, S.; Kutzelnigg, W. Theor. Chim. Acta 1981, 59, 387. (45) Nooijen, M.; LeRoy, R. J. THEOCHEM 2006, 768, 25. (46) Kollmar, C.; Heßelmann, A. Theor. Chem. Acc. 2010, 127, 311. (47) Brueckner, K. A. Phys. Rev. 1954, 96, 508. (48) Nesbet, R. K. Phys. Rev. 1958, 109, 1632. (49) Handy, N. C.; Pople, J. A.; Head-Gordon, M.; Raghavachari, K.; Trucks, G. W. Chem. Phys. Lett. 1989, 164, 185. (50) Urban, M.; Noga, J.; Cole, S. J.; Bartlett, R. J. J. Chem. Phys. 1985, 83, 4041. (51) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (52) Pople, J. A.; Head-Gordon, M.; Raghavachari, K. J. Chem. Phys. 1987, 87, 5968. (53) Stanton, J. F. Chem. Phys. Lett. 1997, 281, 130. (54) Deegan, M. J. O.; Knowles, P. J. Chem. Phys. Lett. 1994, 227, 321. (55) Gwaltney, S. R.; Sherrill, C. D.; Head-Gordon, M. J. Chem. Phys. 2000, 113, 3548. (56) Gwaltney, S. R.; Head-Gordon, M. J. Chem. Phys. 2001, 115, 2014. (57) Gwaltney, S. R.; Head-Gordon, M. Chem. Phys. Lett. 2000, 323, 21−28. (58) Huber, K. P.; Herzberg, G. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979.

(59) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; et al. MOLPRO, version 2010.2, a package of ab initio programs; University College Cardiff Consultants Limited: Cardiff, U.K., 2011. See http://www.molpro.net (accessed July 2012). (60) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M. WIREs Comput. Mol. Sci. 2012, 2, 242−253. (61) Langhoff, S. R.; Davidson, E. R. Int. J. Quantum Chem. 1974, 8, 61−72. (62) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A., Jr. J. Comput. Chem. 1993, 14, 1347−1363. (63) Piecuch, P.; Włoch, M. J. Chem. Phys. 2005, 123, 224105. (64) Ge, Y.; Gordon, M.; Piecuch, P. J. Chem. Phys. 2007, 127, 174106. (65) Neese, F. ORCA, an ab initio, Density Functional and Semiempirical program package, version 2.6; Universität Bonn: Bonn, Germany, 2008. (66) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007−1023. (67) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796−6806. (68) Thouless, D. J. The Quantum Mechanics of Many-Body Systems; Academic: New York, 1961. (69) Knowles, P. J.; Hampel, C.; Werner, H.-J. J. Chem. Phys. 1993, 99, 5219. (70) Deegan, M. J. O.; Knowles, P. J. Chem. Phys. Lett. 1994, 227, 321−326. (71) Szabo, A.; Ostlund, N. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover Publications: Mineola, NY, 1996. (72) Nobes, R. H.; Pople, J. A.; Radom, L.; Handy, N. C.; Knowles, P. J. Chem. Phys. Lett. 1987, 138, 481−485. (73) Knowles, P. J.; Handy, N. C. J. Phys. Chem. 1988, 92, 3097− 3100.

2660


Benchmark Quasi-Variational Coupled Cluster Calculations of Multiple

Recommend Documents