J. Phys. Chem. A 2010, 114, 8687–8696
8687
Electronic Structure Analysis of the Ground-State Potential Energy Curve of Be2† Michael W. Schmidt,*,‡ Joseph Ivanic,§ and Klaus Ruedenberg‡ Department of Chemistry and Ames Laboratory (US-DOE), Iowa State UniVersity, Ames, Iowa 50011, and AdVanced Biomedical Computing Center, Information Systems Program, SAIC-Frederick Inc., NCI-Frederick, Frederick, Maryland 21702 ReceiVed: February 18, 2010; ReVised Manuscript ReceiVed: May 11, 2010
The recently measured ground-state potential energy curve of the diatomic beryllium molecule is reproduced to within an accuracy of 20 cm-1 by a full valence configuration interaction calculation based on augmented correlation-consistent double-, triple-, and quadruple-ζ basis sets, followed by a two-tier extrapolation to the complete basis set limit and complemented by a configuration interaction estimate of the core and core-valence correlations. The origin of binding in Be2 as well as the unusual shape of its potential energy curve is elucidated by an in-depth analysis of the contributions of the various components of this wave function to the bonding process. Beyond the bonding region, the 6/8 London dispersion interaction is recovered. I. Introduction A highly unusual bond is formed by two beryllium atoms inasmuch as the attraction between them results entirely from changes in the dynamic electron correlations. Although weak by average chemical measures, this bond is an order of magnitude stronger than those caused by dispersion forces. Over the last 30 years, quantum chemists1-14 and spectroscopists15-18 together have converged on the ground-state potential curve of this molecule. Pioneering ab initio calculations in 1980 by Liu et al.1,2 first revealed several unexpected features in the interaction potential of these two closed-shell atoms, viz., (i) the potential well is surprisingly deep for a correlation-caused attraction, (ii) correspondingly, its minimum occurs at a bond distance shorter than is typical for van der Waals-bound systems, and (iii) there is a substantial change in slope at about 3.2 Å that results in the distinction of an “inner well” from an “outer well” (although not in two minima). Liu’s early work1 predicted a binding energy (De) of 806 cm-1 and a bond distance (Re) of 2.49 Å. Presently accepted values18 for these quantities are De ) 929.7 ( 2.0 cm-1 and Re ) 2.454 Å, which stand in marked contrast to data19 for the other possible closed-shell second-row molecule, viz., Ne2: De ) 29.4 cm-1 and Re ) 3.091 Å. Experimental confirmation of this rather deep binding potential first came in the mid-1980s15-17 when the first four vibrational levels (of 11 total) of the ground state were measured. This early vibrational data covers, however, only the “inner well” of the potential. Additional quantum chemistry calculations were performed at regular intervals, with advancing accuracy, leading to steadily increasing estimates for De. In 2006, Spirko20 combined the best quantum calculations with the mid-1980s experimental data to produce a potential curve that turned out to be in essentially perfect agreement with the analytic potential extracted from the experiment in 2009.18 As might be expected, the seven additional recently measured vibrational levels18 are highly anharmonic, sampling the “outer well” up to 5 cm-1 below the dissociation limit. All of the vibrational states are plotted in a commentary21 †
Part of the “Klaus Ruedenberg Festschrift”. * To whom correspondence should be addressed. ‡ Iowa State University. § SAIC-Frederick Inc.
accompanying the new experimental measurements. Following the new experiments, Patkowski et al22 have suggested there may be one more vibrational level (V ) 11) just 0.44 cm-1 below the dissociation limit. The special electronic characteristic of Be and the other alkaline earth elements is the easy accessibility to an entirely empty valence p shell. This vacant p shell can influence bonding in two ways: (i) by s,p hybridization, i.e., mixing p-character into the doubly occupied s-orbitals of the restricted Hartree-Fock (RHF) function, and (ii) by a strong near-degeneracy correlation effect through the double excitation s2fp2. Although 2s,2p hybridization certainly occurs in Be2, closed-shell RHF calculations still yield a repulsive potential. In a qualitative comparison of the bonding in all second-row diatomics, Kutzelnigg23 has credited hybridization with keeping this repulsive character small, so that dynamic correlation can bind Be2. Turning to the role of the p orbitals in electron correlation, the beryllium atom has orbital occupancies of (2s)1.80, (2p)0.20 when a multireference wave function is used. Thus, Be makes slightly larger use of the valence p shell than do either Mg [(3s)1.85, (3p)0.15] or Ca [(4s)1.82, (4p)0.18]. It is worth mentioning that data14 for Mg2 (Re ) 3.89 Å, De ) 431 cm-1) and Ca2 (Re ) 4.3 Å, De ) 1102 cm-1) show them to be more similar to Be2 than to the respective rare gas dimers. This large multireference character in the atom, which is also present in the molecule, means both Be and Be2 have an unrestricted Hartree-Fock (UHF) instability.4 Perhaps surprisingly, use of the full-valence-space multiconfiguration selfconsistent-field (MCSCF) wave function (eight orbitals generated by the 2s, 2p shells, which are occupied by four valence electrons) still results in a repulsive potential (except for a shallow well at very long R), which is in fact nearly parallel to the RHF curve.2,6,9 Nonetheless, this wave function, which contains both s,p orbital hybridization and s2fp2 excitation, actually recovers the lion’s share of the correlation energy. However, the true potential shape emerges only when additional excitations outside the valence 2s and 2p orbitals are included, either in extended multireference configuration interaction (CI) calculations1,2,6-9,11 or by high-order many-body calculations.3,5,10,14 Thus, the binding in Be2 is contingent on the effects of
10.1021/jp101506t 2010 American Chemical Society Published on Web 05/27/2010
8688
J. Phys. Chem. A, Vol. 114, No. 33, 2010
dynamical electron correlation, which is uncommon and therefore of considerable interest. It should be remarked that, because of its zeroth-order multireference character, Be2 is much more complicated than Ne2, where one may simply utilize the coupled cluster singles, doubles, and perturbative triples (CCSD(T)) method with correlation-consistent basis sets to extrapolate to a very accurate van der Waals-type potential curve.24 The mentioned investigations have led up to a very accurate ab initio recovery of the experimental spectroscopic data20,22 and, by implication, to a correspondingly accurate quantitative dynamic correlation recovery. By contrast, the goal of the present study is to disentangle the essential correlation features, and their interconnections, that are responsible for this strong correlative bond and for the unusual shape of the dissociation curve, i.e., to provide insight rather than mere numerical accuracy. In order to be credible, insight must of course be based on accurate wave functions, but these must also be amenable to transparent interpretations. To this end, we have constructed consistent sets of such wave functions that are sufficiently systematic to allow the deduction of the essential correlation features. The exhibition of these features is achieved by comparing potential energy curves rather than data tables.7,14 The best wave functions that were generated in the pursuit of this endeavor yielded potential energy curves that are in fact very close to most of the previous high-accuracy results14 (The binding energy at the equilibrium distance is recovered within (20 cm-1 of the best theoretical and the experimental values), although we cannot claim the accuracy of Gdanitz’ r12-CI-based core/valence energies.11 Section II contrasts a number of Be2 potentials generated with standard methods, all using the same modern basis set and briefly exhibits the basis set dependence. Section III presents a very accurate potential computed from valence full configuration interaction (FCI) results, extrapolated to a complete basis set (CBS) limit and complemented by an accurate estimate of core/ valence correlation effects. Section IV, the heart of the study, then develops an in-depth analysis of the dynamic correlations embodied in the valence FCI wave function and leads to the desired elucidations. Section V shows that the calculated potential has, in fact, the dispersion interaction form in the longrange regime. II. Comparison of Methods for Calculating the Valence Correlation Much of the quantitative information in this study will be presented in the form of plots of potential energy curves, computed on a fairly dense adaptive grid.25 Tables of data are not provided, as these are standard calculations, which can be easily regenerated. All calculations were performed with the GAMESS program package26 except for the CCSDT computations, which were done with ACESII.27 II.1. Comparison of Correlation Recovery Methods. It is useful to exhibit how a number of methods, which are popular in current quantum chemistry packages, perform (or fail to do so) for Be2 in order to illustrate the difficulty in treating this system successfully. Many of these methods1-14 have been evaluated earlier (for thorough reviews of previous computations, see references 7 and 14), but some previously unreported potentials are also considered. To facilitate a consistent and reliable comparison, we use the same large basis set, and present the results as figures. Our basis set of choice is a revised version of the aug-cc-pVQZ set,28 containing 5s,4p,3d,2f,1g contracted functions, which are augmented by a set of s,p,d,f,g functions with small exponents.
Schmidt et al.
Figure 1. Comparison of methods, using aug-cc-pVQZ.
A number of single reference methods are explored, since Be2 does dissociate properly into two closed-shell atoms when an RHF molecular wave function is used. Additional potentials based on multireference methods have also been computed. Unless explicitly stated otherwise, the multiconfigurational reference in the present work is always a function in the space of the 112 determinants that are generated in D2h symmetry by the four valence electrons occupying eight molecular orbitals corresponding to the 2s, 2p valence shells of the two atoms. In this section, only correlations in the valence shell are taken into account, the two cores being optimized but correlation inactive. For the single-configuration SCF approximation, both RHF and UHF results are considered. The literature does not contain data for the use of many-body methods based on the use of a UHF reference instead of an RHF reference. Nor is any provided here. One density functional theory (DFT) flavor, namely revised-TPSS,29 is probed, although it is generally accepted that DFT wave functions do not comprise dispersion effects.30 Single reference methods based on closed-shell SCF (RHF) are presented (i) for second-order perturbation theory (MP231), (ii) for coupled cluster at several levels [viz., CCSD,32 CCSD(T),32 the new perturbative triples correction known as CR-CC(2,3),33 perturbative quadruples via CCSD(TQ)-b,34 and full CCSDT35], (iii) for CI36 at the doubles (CISD) and triples (CISDT) levels. Multireference methods include (i) second-order perturbation theory (MRMP37), (ii) multireference singles and doubles CI (MR-CISD38) from MCSCF, and (iii) the FCI from the valence space. The FCI is in fact identical to multireference CISDTQ36 in the four valence electron case of Be2. It is the “exact” result for a given basis set and essentially identical to RHF-based CISDTQ, since the inactive 1s core orbitals are nearly the same in the RHF or MCSCF references. The various methods are compared in Figure 1. The same “exact” FCI potential is shown in black in all panels to help
Analyzing the Ground State Potential Energy Curve of Be2 the eye orient itself to each scale. Figure 1A shows absolute energies, where the energy zero is set at the computed FCI value for two Be atoms. On this large scale, it is difficult to see that the FCI curve shows an “inner” and an “outer” well. On the other hand, it is immediately apparent that the RHF and MCSCF interaction potentials are nearly parallel, and that the neardegenerate 2s2 f 2p2 excitations recover the great majority of the total electron correlation. The MCSCF potential is, however, still repulsive, except for a very shallow minimum around R ) 5.5 Å, which implies that the attractive interaction potential of Be2 requires dynamical electron correlation effects. Figure 1A also illustrates the UHF instability due to the large neardegeneracy correlation. Note that the RHF and UHF curves do not dissociate to the same asymptote, since this instability originates in the atomic 2s,2p near-degeneracy. Finally, there are three SD-excited methods shown in this panel. The size consistency of CCSD means it dissociates to the correct asymptote, namely, that for FCI, whereas the lack of size extensivity of CISD causes it to dissociate to a higher energy. MP2 is also size extensive, dissociating correctly to twice the atomic MP2 value, but perturbative recovery of the correlation energy of the valence two-electron atom is not quantitative. On the remaining panels in Figure 1, the asymptotic energy for each method is shifted to the same dissociation limit. Figure 1B shows in greater detail how closely the MCSCF potential parallels the RHF potential. It also shows that FCI is a good approximation to the experimentally derived potential,18 with an inner well appearing inside R ) 3.2 Å. It is clear that the UHF and DFT potentials dramatically overbind, and UHF underpredicts Re. The particular DFT functional chosen, namely revTPSS, Perdew’s most recent nonempirical (also nonhybrid) meta-GGA,29 yields a potential that does not exhibit a slope change. Note that MP2, a formally less exact method than CCSD since it does not have converged amplitudes, manages to closely track the outer well, but does not possess the inner well. CCSD has an interesting “wiggle” with a shallow well at large distance. CISD might be considered the most successful doubles method investigated here inasmuch as it has a small inner well. However, none of the doubles methods are satisfactory. The CISDT method (not shown to avoid clutter) has both inner and outer wells, but slightly overbinds with a binding energy of 970 cm-1 (at R ) 2.5 Å) compared to 786 cm-1 for FCI. Figure 1C illustrates the performance of single reference coupled-cluster methods. Fully converged CCSDT does partially incorporate quadruples effects and closely follows the FCI curve. Approximate treatments of triple excitations, namely CCSD(T) and the recently developed CR-CC(2,3) method, are not as accurate. Perhaps surprisingly, since Be2 is a valence fourelectron system, the perturbative treatment of both triples and quadruples by CCSD(TQ)-b is not any better. Of course, in this case, fully iterated CCSDTQ would be equivalent to FCI. Figure 1B/1C shows that excitation levels higher than doubles with respect to RHF are necessary to reproduce the qualitative shape of the Be2 interaction potential. Figure 1D focuses on multireference methods. The interaction potential from MR-CISD is nearly indistinguishable from that of FCI, showing that, at most, two electrons need to be excited from the MCSCF reference of the valence 2s and 2p orbitals in order to obtain good accuracy. However, perturbative treatment of these external double excitations, via the MRMP approach, results in a disappointing interaction potential that contains a more pronounced wiggle than that of CCSD. To summarize, qualitatively correct curves, possessing both inner and outer wells, are obtained only if the total electron
J. Phys. Chem. A, Vol. 114, No. 33, 2010 8689
Figure 2. Illustration of the basis set requirements, using unaugmented basis sets.
excitation level with respect to RHF exceeds 2, although no more than two electrons need be excited into orbitals higher than those constructed from the valence 2s and 2p shell. Theories that are limited to double excitations with respect to RHF, or that fully excite only within the valence 2s and 2p orbitals (i.e., valence-MCSCF) are inadequate. The MR-CISD and CCSDT methods are the closest approximates to FCI, with MR-CISD being almost exact. II.2. Comparison of Basis Sets. It is well-known that large basis sets are needed to obtain accurate results for Be2,1,2,6,7,10-12,14 and these requirements are reviewed in Figure 2. Results are presented for the multireference valence FCI in each basis, so that the correlation method is “exact”. Since the differences in basis sets are more visible in the absence of augmentation by diffuse functions, Figure 2 shows only results for the systematic basis set sequence28 cc-pVDZ, cc-pVTZ, and cc-pVQZ, as solid curves. Also shown, as dashed curves, are two basis sets that are intermediate between the DZ and TZ sets. It is known6 that omitting all basis functions higher than p produces an FCI curve that is entirely repulsive. The DZ basis contains just s,p,d basis functions and exhibits an incorrect double well. Adding the f orbital of the TZ set to the DZ set produces the dashed green curve DZ(+f), while deleting this f orbital from the TZ set produces the dashed blue curve TZ(-f). In both cases, this f basis function has a 325 cm-1 impact around Re, i.e., on the binding energy, which is why Liu’s seminal work1,2 included an f STO. On the other hand, going from DZ to TZ adds one function each of s,p,d,f type, with the smallest s,p,d exponents becoming smaller. It is apparent that these s,p,d improvements in the TZ(-f) curve remove the double minimum artifact in DZ as well as most of the quantitative error around R ) 4 Å without help from the f orbital. Almlo¨f7 has noted the importance of saturating the small exponent s,p,d orbital space. The diffuse augmentation remains important even at the QZ size. The FCI curve utilizing the aug-cc-pVQZ basis in Figure 1D closely hugs the red experimental curve at long R, while that employing the unaugmented cc-pVQZ curve in Figure 2 does not. At R ) 2.45 Å, the diffuse augmentation to the QZ basis contributes 41 cm-1 to the binding energy of Be2. In summary, accurate potential curves require saturation of the low-angular-momentum diffuse exponent function space, particularly at long interatomic separations, as well as inclusion of f (and higher) angular momentum functions to recover dynamic correlation around the equilibrium distance. The unaugumented 4s,3p,2d,1f basis set cc-pVTZ is the smallest basis qualitatively capable of generating inner and outer wells.
8690
J. Phys. Chem. A, Vol. 114, No. 33, 2010
III. Accurate Ground-State Potential Energy Curve With the exception of the noble gases, where the CCSD(T) approach is very effective,24 the determination of near-quantitative potentials for diatomics with many electrons, such as C2,39 N2,39 O2,40 F2,41 requires considerable effort, such as the correlation energy extrapolation by intrinsic scaling (CEEIS42) to the FCI limit. In Be2, by contrast, with only four valence electrons, the execution of the full valence CI calculation is quite feasible so that quantum chemistry leads in a straightforward manner to a very accurate potential. The best of the potential curves discussed in the preceding section was that obtained by the full valence CI calculation using the aug-cc-pVQZ basis set. This valence FCI wave function will be used in section IV, to understand the bonding in Be2. In the present section, we demonstrate that the remaining errors in this potential can be understood, and that they do not affect the conclusions about bonding. The differences between this potential and the experimental curve, which are manifest in Figures 1C,D, must be due to the following shortcomings: (i) lack of a complete atomic orbital (AO) basis since aug-cc-pVQZ only includes up to g-type AOs, (ii) lack of electron correlation effects involving the four core electrons, (iii) lack of relativistic effects, from scalar interactions and spin-orbit coupling. III.1. Relativistic Effects. Consider first the scalar relativistic effects, viz., the mass-velocity and Darwin corrections. Calculating them via the third-order Douglas-Kroll-Hess transformation of the one-electron integrals,43 we find nearly the same corrections at all bond distances, e.g., 1057 cm-1 at R ) 2.5 Å and 1061 cm-1 at R ) 20.0 Å, for both the MR-CISD or FCI wave functions. Although a 4 cm-1 correction to De is within the noise level of the basis set and core correlation estimates, this correction is easily evaluated at all distances (at the MRCISD level) and is included. Regarding spin-orbit coupling, we note that both the orbital angular momentum as well as the spin angular momentum vanish for the 1Σg+ molecule and at the 1S atomic limit. Hence, this interaction is considered to be negligible. III.2. CBS Limit of Valence FCI. The basis set error must be treated by separate extrapolations of the MCSCF and correlation energy values obtained with the correlation consistent basis set series aug-cc-pVDZ, aug-cc-pVTZ, and aug-cc-pVQZ. The MCSCF energies are extrapolated using the exponential form commonly used for SCF and MCSCF energies:44
EMCSCF(∞) ) EMCSCF(X) + ae-bX The correlation energies (EMCSCF - EFCI) are extrapolated by inverse powers:45
Ecorr(∞) ) Ecorr(X) +
A B + 5 3 X X
Each formula requires three successive energies, taken from the augmented basis sets X ) 2,3,4, to determine the three constants. The CBS limit is then given by
ECBS ) EMCSCF(∞) - Ecorr(∞) At R ) 2.45 Å, the MCSCF energy extrapolation lowers the total energy of Be2 by 60.9 cm-1. However, it lowers the energy of one Be atom by 32.8 cm-1, thus reducing the binding energy
Schmidt et al. by 4.7 cm-1. Correlation extrapolation for Be2 increases the valence-only dynamic correlation by 131.5 cm-1, to a CBS limit of 4221 cm-1. Extrapolation for a Be atom increases the atomic dynamic correlation by 33.8 cm-1, to a CBS limit of 540 cm-1. Correlation extrapolation thus increases the binding energy by 63.9 cm-1. III.3. Core Correlation. Core correlation corrections are of similar magnitude as the valence CBS extrapolation. Previous work in this laboratory41 as well as that of other groups46 relied on calculations using the cc-pCVQZ basis set47 and an active space of (1s,2s,2p) followed by an all-electron MR-CISD with a Davidson +Q correction. The present work uses the new, recently developed cc-pwCVQZ basis set,28 which emphasizes the recovery of core/valence correlations over core/core correlations. Moreover, a different configurational scheme is used, which is faster and perhaps more accurate. The core correlation is taken as the energy increment from our standard MCSCF energy to the corresponding core-excited-only MR-CISD energy, which results from a computation allowing single and double excitations out of the two 1s core orbitals into (i) the active four-orbital (2s,2p) valence space and (ii) the external, dynamically correlating orbital space beyond 2s,2p. In the counting language of the occupation restricted multiple active space (ORMAS) CI program,36 the active spaces have the following size and occupation specifications: Space (i): the two core orbitals contain between two and four electrons; Space (ii): the eight valence orbitals possess between four and six electrons; Space (iii): the external correlating orbitals contain between zero and two electrons. Thus the valence electrons are only excited within the valence 2s,2p orbital space, but not into the external correlating space. Note that this core correction does not include simultaneous core/valence excitations into the external space. These core correlations are then added directly to the CBS estimate of the valence FCI to produce the estimated Be2 potential. The accuracy of this scheme was tested by a benchmark MRCISDTQ calculation (ORMAS occupancies: core ) 0-4, valence ) 0-8, external ) 0-4) in which every conceivable core-core, core-valence, and valence-valence correlation effect is recovered up to the level of quadruple excitations. The calculation is so large that it could only be done with the double-ζ basis cc-pwCVDZ (even then involving 63 million determinants). This benchmark calculation results in the energy difference [E(R ) 2.5) - E(R ) 20)] ) -221 cm-1, while the procedure described in the preceding paragraph, when used with the same cc-pwCVDZ basis, yields the value -225 cm-1 for this energy difference. These binding energies correspond to core correlation corrections of 69 cm-1 (MR-CISDTQ) and 73 cm-1 (new procedure). This near-identity implies that simultaneous core/valence excitations beyond the active orbitals recover only intra-atomic correlations (at great effort) and, hence, generate only a curve-nonparallelity of 4 cm-1 from the described approximate accounting for the core correlation contributions to the potential curve. (Since DZ bases are rather poor, the valence FCI binding energy of -152 cm-1 is, of course, rather far from the value of -786 obtained with the larger basis aug-cc-pVQZ.) We believe that this approach to core correlation corrections can be generalized to other systems. It relies, however, on the use of a full valence active space and on basis sets containing reliable core-correlating functions. The core correlation increments obtained with the ccpwCVQZ basis set are 18 526 at R ) 2.45 and 18 456 cm-1 at R ) 20 Å. These correlations of ∼84 millihartree (mh) are somewhat smaller than the corresponding values in F2 (∼122 mh)41 and in O2 (∼124 mh).40 They increase the binding energy
Analyzing the Ground State Potential Energy Curve of Be2
Figure 3. The CBS limit curve includes basis extrapolations to both MCSCF and correlation energies. Additional corrections consist of core correlation and scalar relativity (see text). The FCI curve, drawn in black in Figures 1 and 2, is shown in purple here.
in Be2 by 70 cm-1 (0.320 mh), which is similar in magnitude to the contributions in F2 (decreasing binding by 0.314 mh)41 and in O2 (increasing binding by 0.503 mh).40 Meyer9 estimated the core contribution to binding in Be2 to be 83 cm-1 from a core polarization potential. Almlo¨f7 estimated 89 cm-1 by a more direct calculation. More recently, Martin10 gave estimates in the range 71-82 cm-1 using various CI and coupled cluster methods. Gdanitz11 explicitly correlated the core with r12 methods and obtained a directly calculated De ) 898 ( 8 cm-1 without basis set extrapolation or additivity schemes. He did not provide a separate estimate of the core correlation correction. Szalewicz14 estimated core correlation corrections to be 85 cm-1. Thus, the present value of 70 cm-1 is smaller than most previous estimates, perhaps because our work uses the newly developed cc-pwCVQZ core/valence basis set for Be.28 III.4. Ground State Potential Energy Curve. After correction of the valence basis set deficiencies by CBS extrapolation and inclusion of the core electron correlation effects, the final potential shown in Figure 3 is a close match to experiment. At the distance R ) 2.45 Å, we obtain the following accounting of the various contributions to the binding energy ∆E ) [E(R ) 2.45) - E(R ) 20)], expressed in cm-1: RHF-SCF energy difference without p,d,f,g hybridization Hybridization in RHF-SCF with p,d,f,g from aug-cc-pVQZ basis
+8139.7 -5513.0
RHF-SCF energy difference ∆E with aug-cc-pVQZ basis near-degeneracy correlation within 2s-2p valence shell
+2626.7
full valence shell MCSCF value for ∆E with aug-cc-pVQZ basis dynamic valence shell correlation
+2290.3
- 336.4
-3076.6
full CI binding energy ∆E with aug-cc-pVQZ basis scalar relativistic correction CBS extrapolation of MCSCF CBS extrapolation of dynamic valence correlation core correlation energy contribution
- 786.3
total estimated binding energy ∆E at 2.45 Å
- 911.7
+ 3.8 + 4.7 - 63.9 - 70.0
The experimental estimate of the binding energy (at R ) 2.4536 Å) is -929.7 ( 2 cm-1. The deviation of theory from 18
J. Phys. Chem. A, Vol. 114, No. 33, 2010 8691
Figure 4. Be2 compared to He2. Solid lines use the full aug-cc-pVQZ basis, while dashed lines use only the s AOs from this set. The distance axis is scaled by each atom’s single bond radius, d ) R(X2)/R1(X).
experiment, which is about the width of the curves displayed in Figure 3, corresponds to the theoretical error bars that are associated with the CBS extrapolation, the core correlation approximation, and, possibly, the neglect of spin orbit coupling. Our interaction energy of about -912 cm-1 can be compared to the directly computed result -898 cm-1 of Gdanitz,11 and to the -939 cm-1 reported by Szalewicz,14 which is based on similar corrections involving CBS extrapolations and core/ valence correlation estimates. IV. Dynamic Correlation, Binding, and the Shape of the Potential IV.1. Hybridization and Near-Degeneracy Correlation. Although the RHF curve is still repulsive, 2s,2p orbital hybridization decreases this repulsion sufficiently so that addition of the much smaller dynamic correlation is able to produce a bound potential. Figure 4 illustrates the impact that hybridization has on the RHF potential of Be2, as was already anticipated by Kutzelnigg.23 This figure compares the RHF potential of Be2 to that of He2, which lacks an accessible p orbital. In order to facilitate the comparison, the x-axis is scaled by expressing the bond distance for Be2 and He2 as multiples of their respective single bond radii,48 namely, R1(Be) ) 1.02 Å and R1(He) ) 0.46 Å. The dashed curves are the potentials that result when p and higher angular momentum functions are deleted from the aug-cc-pVQZ basis of each system so that they cannot hybridize. The solid curves include hybridization (the solid Be2 curve is thus identical with the RHF curves in Figure 1). It is apparent that (i) the potential curve of Be2 parallels that of He2 when both systems are constrained to s-type functions only; (ii) the inclusion of basis functions other than s-type has only a very small effect on the potential of He2; (iii) by contrast, the potential of Be2 is shifted downward by many thousands of wavenumbers throughout its bonding region when hybridization with p basis functions is possible. The energy lowering at R ) 2.45 is 5513 cm-1, a value that was already included into the energy summary given in the previous section. Thus, as Kutzelnigg predicted,23 this lowering of the RHF potential by hybridization is an essential contributing factor in making it possible for electron correlation to overcome the RHF repulsion and bring about bonding. In addition to hybridization, the MCSCF potential also contains the effects of nondynamic correlations due to excitations from the SCF orbitals into the valence orbitals generated
8692
J. Phys. Chem. A, Vol. 114, No. 33, 2010
Schmidt et al.
TABLE 1: Atomic Be Natural Orbital Occupation Numbers per Angular Shell and Individual Orbital, from FCI/ aug-cc-pVQZ orbital
total shell
individual orbital
1s 2s 2p 3s 3d 4f 3p 4s 4d 5g 5f 4p 5p 6g 5s 5d 6f 6s 6d 6p
2.0 1.8194048850 0.1765742463 0.0032187492 0.0005653645 0.0000789838 0.0000630279 0.0000301106 0.0000255515 0.0000181440 0.0000075334 0.0000035778 0.0000034986 0.0000024129 0.0000017696 0.0000013140 0.0000008162 0.0000000103 0.0000000035 0.0000000000
2.0 1.8194048850 0.0588580821 0.0032187492 0.0001130729 0.0000112834 0.0000210093 0.0000301106 0.0000051103 0.0000020160 0.0000010762 0.0000011926 0.0000011662 0.0000002681 0.0000017696 0.0000002628 0.0000001166 0.0000000103 0.0000000007 0.0000000000
by the 2p shell. This near-degeneracy correlation energy lowering also lowers the energy by thousands of wavenumbers, as was displayed in Figure 1A. Surprisingly, however, as is apparent from Figure 1B, it hardly changes the shape of the potential, and it does not remove the repulsive character of the potential (except for a very shallow minimum at large R). IV.2. Dynamic Correlation. Details of dynamical correlation are normally considered difficult to analyze because many small contributions are involved, being energy increments stemming from many small CI coefficients. In Be2, it is possible to design a systematic sequence of calculations, incorporating increasing amounts of dynamical correlation, to understand the origin of its molecular binding. The key to this is Be’s atomic electron correlation. Table 1 lists the natural orbital occupation numbers for Be, computed using valence FCI (equivalent to CISD) with the aug-cc-pVQZ basis set. There are 1.996 electrons in the Valence shell of the atom, so the molecular valence space occupancy is 3.992 for two atoms at the dissociative limit. This may be compared to 3.976 in the valence natural orbitals of Be2’s FCI at R ) 2.45 Å. It is this slightly larger degree of excitation into the external orbitals that lowers the molecule’s dynamic correlation energy enough to bind the atoms. Table 1 suggests a strategy to consistently enlarge the orbital space beyond the valences 2s and 2p. In fact, Meyer’s inclusion of the two molecular orbitals arising from the 3s orbital of Be in a MCSCF reference9 must have been motivated by its clear predominance in the external space of the atom. Note that the remaining orbitals are not ordered as one might expect. The orbitals 3s, 3d, and 3p trail off in importance, by a factor of about 6 and then about 9, rather than being roughly similar. In fact, the 4f shell as a whole is more important than the 3p shell. Thus, whereas the constructions of modern basis set families, such as the approximate natural orbital or the correlation consistent bases (which are used in the present work), are derived from the notion that shells with the same principal quantum number make similar contributions to the correlation energy, this type of quasi-degeneracy concept does not hold up for the n ) 3 shell of the beryllium atom. Guided by Table 1, we shall add molecular orbitals to the Be2 MCSCF ansatz one atomic shell at a time, e.g., adding σg x σu orbitals to include both 3s orbitals of the atoms. Specifically, the order in which atomic shells are added is
Figure 5. Systematic increase in the orbital spaces that afford the dynamic correlation in the context of MCSCF calculations.
dictated by the anticipated energy impact of the orbitals, as the MCSCF process will change the shape of orbitals to minimize the energy. The energy impact of any orbital is then presumed to be predicted by its occupation number. Thus, our criterion is the occupancy of individual orbitals, not shells, since MCSCF optimization acts on individual orbitals. Therefore, the order in which we add new shells is taken from the final column of Table 1, which lists the shell’s occupancy divided by the shell’s degeneracy, namely, 3s, 3d, 4s, 3p, 4f, 4d, and 5g. It is straightforward to successively add these shells at long bond distances, in exactly that order, and obtain molecular orbitals that possess the desired atomic character. For example, increasing the active space by σg x σu x πu x πg x δg x δu orbitals converges on orbitals that exhibit 3d character at long R, provided this is done after the +3s shell has been added (otherwise the σg x σu part of the +3d extension will collapse to 3s shapes). The orbitals at long R with the desired atomic 3d character can then be used as starting orbitals for the same MCSCF calculation at successively shorter bond distances. Even though these orbitals will begin to hybridize as the distance R shortens, they do not lose their predominant “3d” character. All such shell by shell extensions of the active space yield continuous potentials. We shall use the notation {MCSCF(4,x)+nl} to refer to exciting the four valence electrons within x orbitals, indicating by “+nl” the last shell added. (The entire function space of aug-cc-pVQZ would therefore be reached for x ) 158.) The potential curves for such systematic shell-driven increases in the MCSCF active space are displayed in Figure 5, where the asymptotic energy of each curve is shifted to the same dissociation limit. Even though the 3s shell has the largest occupation number and generates a correspondingly substantial energy lowering, the +3s curve remains essentially parallel to the valence MCSCF curve. Manifestly, the orbitals that arise from the 3d shell are the ones that account for the presence of the bond in Be2, as the MCSCF(4,20)+3d curve is the first to exhibit binding as well as the inner and the outer sections of the potential curve. Beyond the +3d curve, the largest energy impact is provided by +4f, being some 148 cm-1. This is smaller than the ∼325 cm-1 impact of the first f basis function shown in Figure 2. Presumably, hybridization of orbital shapes at finite R incorporates some of the effect of the f basis function into earlier shell increments. Beyond +4f, there is a steady convergence toward the basis set’s FCI limit, as shells of gradually decreasing importance are added.
Analyzing the Ground State Potential Energy Curve of Be2 TABLE 2: Be2 Natural Orbital Occupation Numbers from MCSCF(4,20)+3d/aug-cc-pVQZa Calculations interatomic distance (Å) orbital
2.5
3.0
3.2
3.4
5.0
σg σu σg πu πg σu 2s,2p total σg σu +3s total πu σg σu πg δg δu +3d total grand total
1.8593 1.7107 0.1750 0.0649 0.0415 0.0210 3.9788 0.0059 0.0039 0.0098 0.0036 0.0019 0.0005 0.0004 0.0004 0.0002 0.0080 4.0002
1.8413 1.7716 0.1042 0.0651 0.0524 0.0324 3.9845 0.0043 0.0035 0.0078 0.0020 0.0018 0.0003 0.0003 0.0003 0.0002 0.0077 4.0000
1.8357 1.7846 0.0907 0.0646 0.0548 0.0368 3.9866 0.0040 0.0033 0.0073 0.0016 0.0016 0.0003 0.0003 0.0002 0.0002 0.0065 4.0004
1.8312 1.7932 0.0822 0.0639 0.0564 0.0408 3.9880 0.0038 0.0032 0.0070 0.0012 0.0013 0.0002 0.0002 0.0002 0.0002 0.0051 4.0001
1.8186 1.8121 0.0633 0.0606 0.0598 0.0572 3.9920 0.0035 0.0030 0.0065 0.0002 0.0003 0.0001 0.0001 0.0001 0.0001 0.0014 3.9999
a
Rounding to four digits leads to small errors in the sums. The actual calculations had of course exactly four electrons. The occupancies of doubly degenerate shells are given for one of their orbitals, with atomic shell totals including such occupations twice.
These inferences are confirmed by an examination of the occupation numbers for the MCSCF(4,20)+3d active orbitals, which are listed in Table 2 for representative bond distances. An examination of the contours of these orbitals confirms that they possess the expected atomic shell ordering 2s, 2p, 3s, 3d and, moreover, the +3s and +3d orbitals still possess s and d shapes at any bond distance. At large distances, the molecular occupation numbers closely resemble the atomic occupancies given in Table 1. As the bond distance is shortened from R ) 5 Å, the valence 2s, 2p shell shows the typical increase in occupation numbers for the bonding orbitals σg and πu coupled with the coresponding decreased use of the antibonding σu and πg. In addition, hybridization of 2s and 2p takes place within the σ orbitals, which also causes some redistribution of the σ occupation numbers. Finally, since dynamical correlation (excitation into orbitals outside the 2s,2p shell) increases as the bond forms, the total electron count in the valence shell decreases slowly. While these occupancy changes may seem small, they are nonetheless entirely responsible for the bond in Be2. The first shell in the external space is formed from 3s AOs. But, its occupancy increases only by 50% when the bond is formed. Since this is a small relative change, the shell does not contribute to the bond energy, as the +3s curve in Figure 5 shows. In contrast, the +3d shell orbitals increase their occupation number by a factor of 7, from 0.0011 at R ) ∞ to 0.0080 near Re. Notably, the bonding orbitals of πu-type (from dxz, dyz) and of σg-type (from dz2) change much more than the others. The growth in importance of the πu accelerates over the range R ) 3.4-3.0, where the potential bends from its outer slope into the inner well. Also, the contributions from diffuse dxz (or dyz) AOs to the πu orbital doubles in size between R ) 5.0 and R ) 3.2, before decreasing back nearly to its R ) 5 value at R ) 2.5 Å. The σg’s diffuse dz2 coefficient steadily decreases from R ) 5 to R ) 2.5, declining to about 1/3 of its original value. These large relative changes in the +3d shell occupations confirm the role of this shell in generating the binding energy.
J. Phys. Chem. A, Vol. 114, No. 33, 2010 8693 TABLE 3: Atomic Shell Occupation Numbersa orbital
He
Be
Ne
1s 2s 2p 3d 3p 3s 4f 4d 4p 4s 5g 5f 5d 5p 5s 6g 6f 6d 6p 6s
1.9838894300 0.0076085695 0.0076928880 0.0003297475 0.0002494758 0.0001148404 0.0000409654 0.0000299520 0.0000194988 0.0000106197
2.0 1.8194048850 0.1765742463 0.0005653645 0.0000630279 0.0032187492 0.0000789838 0.0000255515 0.0000035778 0.0000301106 0.0000181440 0.0000075334 0.0000013140 0.0000034986 0.0000017696 0.0000024129 0.0000008162 0.0000000035 0.0000000000 0.0000000103
2.0 1.9888943850 5.9330929176 0.0219996815 0.0387434352 0.0079628739 0.0027913816 0.0023990740 0.0020590002 0.0003334311 0.0006161175 0.0004145589 0.0002269810 0.0002751606 0.0000094395 0.0000797823 0.0000503153 0.0000377675 0.0000119475 0.0000017494
0.0000073395 0.0000051180 0.0000015294 0.0000000255
a Occupation numbers for each shell, using valence FCI for He and Be and CI-SDTQ for Ne, for aug-cc-pVQZ. Shells higher than n ) 2 have been sorted by angular momentum, in reverse order. Note that the highest n shell is formed using the “aug” part of the basis set, and so might not correspond to results that would be obtained with cc-pV5Z.
In the MCSCF(4,20)+3d function the 10 orbitals of the +3d shell are the highest active molecular orbitals and, as mentioned above, they look like diatomic orbitals constructed from 3d AOs. This is because, even at shorter distances, the s,p,f,g admixtures into them during MCSCF orbital optimization remain smaller than their d components. Is this weighting of the correlating orbitals in fact also the case for the corresponding natural orbitals of the actual FCI wave function? Checking the natural orbitals of the FCI wave function at R ) 2.45, one finds that the first eight of these arise from 2s and 2p orbitals, with the next two having predominantly 3s character. Of the following orbitals, the next eight have in fact d shapes as do the 11th and 12th NO. The only difference is that the 9th and the 10th NOs are 3p-like. The NO structure of the accurate FCI wave function therefore closely resembles that of its MCSCF(4,20)+3d approximant. IV.3. Why Do the 3d (and 4f) Angular Correlations Dominate the Bonding in Be2? It is apparent from the preceding discussion that the relative importance of the various correlating orbitals in the molecule is a consequence of their importance in the Be atom, as indicated by the intra-atomic NO occupation numbers. Table 3 contains orbital shell occupations for the closed-shell atoms He, Be, and Ne. For He,49 the orbitals naturally divide into shells of nearly equal importance, starting with a correlation shell containing 2s,2p, then 3s,3p,3d, and finally 4s,4p,4d,4f. Within any set for a given principal quantum number, the occupation numbers differ by no more than four, and the results are perfectly sorted by placing them in the reverse order of the angular momenta. For Ne, the same ordering scheme generally applies, although the ns orbital’s occupation is closer to the n+1 shell than to its own. Neon also occupies the np orbital more than nd, for n ) 3 and 5, because a p orbital can provide both angular correlation to 2s and radial correlation to 2p. Beryllium differs from He and Ne by having the 2p orbitals available for angular correlation in the same energetic shell as the occupied 2s orbital, and this 2p-near-degeneracy-type of
8694
J. Phys. Chem. A, Vol. 114, No. 33, 2010
Schmidt et al.
TABLE 4: Long-Range Be2 Potentialsa FCI R
RHF ACCQ
4.0 4.5 5.0 5.5 6.0 6.5 7.0 8.0 9.0 10.0 C6 C8
192.916 71.974 24.893 8.055 2.415 0.626 0.109 -0.054 -0.059 -0.039
b
MCSCF ACCQ
ACCD
ACCT
ACCQ
CBS limitc
estimated potentiald
analytic fit to experimente
149.709 35.477 1.310 -5.924 -5.748 -4.242 -2.920 -1.366 -0.693 -0.375
-145.454 -127.109 -91.168 -59.286 -37.575 -24.535 -16.850 -8.785 -4.912 -2.826
-196.136 -137.953 -89.296 -55.100 -33.591 -20.828 -13.381 -6.273 -3.308 -1.840
-196.742 -135.829 -86.586 -52.811 -31.806 -19.340 -12.031 -5.095 -2.435 -1.268 228 11400
-195.034 -134.977 -85.711 -52.014 -31.146 -18.784 -11.502 -4.474 -1.907 -0.905
-200.108 -136.931 -86.447 -52.297 -31.262 -18.839 -11.534 -4.491 -1.919 -0.913 146 22800
-199.625 -135.447 -85.306 -50.891 -29.401 -16.689 -9.384 -2.931 -0.910 -0.282
a Distances are in Å, energies are in cm-1, dispersion coefficients Cn are in Hartree-Bohrn. The latter are obtained by linear fits to the three points R ) 6.5, 7, and 8 Å. The theoretical potentials are obtained as V(R) - V(20 Å), while the analytic formula from experiment uses V(R) - V(∞). b The notation ACCX is shorthand for aug-cc-pVXZ, X ) D,T,Q. c Using the two-tier extrapolation of both MCSCF and correlation energies, but no other corrections. d The estimated potential includes CBS extrapolation, core/valence correlation, and relativistic corrections as described in section III. e See ref 18. This analytic fit has an exponential decay toward R ) ∞.
correlation downgrades the value of the correlation that the higher 3p and 4p orbitals can proffer. In addition, as discussed above, the near degeneracy of the 2p orbitals also leads to their substantive hybridization mixing into the Hartree-Fock orbital. This upgrades the value of the 3d correlation since it provides the most important angular correlation to the 2p admixture into the Hartree-Fock orbitals. This importance of the 3d orbitals relative to the 3s and 3p orbitals is much larger in the Be2 wave function than the role played by the 3d orbitals in Ne2, as we confirmed by a comparative examination of the two molecules.50 As noted above, Table 2 moreover reveals the following difference between the contributions of the two πu-type orbitals (coming from dxz, and dyz) and that of the σg-type orbital (coming from dz2). All of them increase as the atoms approach each other in the outer region of the potential well from 5 Å up to around 3.2 Å. From there to the equilibrium region at 2.5 Å, the two πu-type contributions further almost double, while the σg-type contribution hardly changes anymore. This increasing predominance of the πu-type therefore appears to be related to the marked change in slope of the potential curve at 3.2 Å. In the absence of a more detailed energy decomposition, one might speculate that the πu-type 3d may have less repulsion by the other atom at shorter distances due to its node directly along the bond axis. These relative changes in the comparatively small correlation contributions of a few 3d orbitals and some 4f orbitals would not, however, be able to determine the shape of the potential energy curve and the binding energy, if it were not for the extraordinary insensitiVity of the much larger contributions provided by the 2p-type angular correlation and the 3s-type radial correlation as the internuclear distance changes from equilibrium to dissociation. So far, we have not gained any insight into why these energy lowerings of close to 20 000 cm-1 change so little with the internuclear distance. V. Long-Range Potential In the very long-range region, the dynamic correlation is expected to turn into the London dispersion forces51 between the two Be atoms, which confer to the potential V(R) the distance dependence
V(R) ) -C6/R6 - C8/R8 - C10/R10...
Recent experiments18 have established the level for V ) 10, which is bound by 4.8 cm-1 and has a classical turning point near R ) 7.9 Å. It therefore clearly samples the long-range potential experimentally. Nonetheless, the experimentalists reported18 that, among the available analytical potentials, a generalized Morse potential with an exponential long-range decay yielded a better overall fit to the 10 experimental level differences than a generalized Morse potential hybridized with a built-in R6/R8 decay.52 Since they report an absolute mean fitting error of 0.5 cm-1, which is 10% of the binding energy of the highest level mentioned above, it would appear that the fitting optimization was dominated by the lower levels so that an experimental determination of the dispersion constants in Be2 is still outstanding. Since the publication of these experiments, a new study22 based on the potential of Szalewicz14 has appeared, suggesting the existence of a V ) 11 level, bound by just 0.44 cm-1. The dispersion-type asymptotic shape of the potential is crucial to the existence of such a level. In ref 24, it was shown that the coefficient C10 cannot be determined from a long-range potential energy curve unless the latter is known with an accuracy of better than 10-9 hartree and that, because of the mutual linear dependence of R-8 and R-10, the value obtained for C8 depends on the arbitrary choice of C10. In the present context (as in ref 24), one can therefore omit C10 as well as all higher coefficients without loss of fitting accuracy. The dispersion expansion can then be written as
-R6V(R) ) C6 + C8(R-2) and the dispersion character of our theoretical V(R) can be tested by plotting [-R6 V(R)] versus (R-2). If it turns into a straight line for (R-2) f 0, then the values of C6 and C8 are obtained as the intercept and slope of its linear least mean squares fit. Several long-range potentials are listed in Table 4. The RHF potential goes to zero by R ) 8 Å, but the MCSCF and CI curves extend several angstroms further out. Since the augmented double and triple-ζ potentials approach the augmented quadruple-ζ potential from below, the complete basis limit curve may be less reliable regarding the curVe shape at long distances than the directly computed FCI/aug-cc-pVQZ results. This may be a symptom of the basis set superposition error in the smaller, notably the double-ζ, basis sets. While the CBS, core/valence,
Analyzing the Ground State Potential Energy Curve of Be2
Figure 6. Demonstration of London dispersion form. The estimated potential (including CBS extrapolation and core/valence correlation) is an improvement over the directly computed FCI/aug-cc-pVQZ potential in the bonding region, but not at the large R limit. Linear fits to extract Cn data used the points at R ) 6.5, 7, and 8 Å, ignoring the last two points at R ) 9, 10 Å. Units on x axis ) bohr-2, units on y axis ) hartree-bohr.6
and relativity corrections manifestly improve the direct FCI at small R values, this is not assured to be so in the delicate region governed by dispersion forces. The core/valence correlation and relativistic corrections are, however, negligible at large R. Figure 6 therefore displays dispersion plots for the following three potentials from Table 4: FCI-ACCQ (black), estimated potential (green), and (in red) the potential inferred from the vibrational spectrum in reference.18 It is apparent that all theoretical curves display dispersion character from about R ) 6.5 Å out. Also shown (blue) is the total electron correlation energy potential Vcorr ) V(RHF) - V(FCI), which expectedly goes into the 6-8 form at somewhat shorter distances and then becomes identical with the FCI curve. The last two rows of Table 4 list the values of C6 and C8 that are obtained by leastmean-squares fitting. As stated above, we consider the values from the uncorrected FCI/aug-cc-pVQZ results as the most trustworthy. The most recent (2006) calculations of the dispersion coefficients53 of Be2 using atomic polarizabilities have given the values C6 ) 214 ( 3, C8 ) 10 230 ( 60, and C10 ) 504 300. These authors also cite older values, which are rather similar to their data. Sta¨rck and Meyer’s direct fit9 to their CI potential yielded C6 ) 212.8, C8 ) 10 200, C10 ) 487 400. It is gratifying that the C6 value from fitting our directly computed FCI data, 228 hartree-bohr,6 lies within 7% of the polarizability-derived value. The recovery of long-range dispersion interactions by our best FCI potential makes it again manifest that the correlation effects that determine binding near the equilibrium distance are of a very different nature. VI. Conclusions A very accurate potential for Be2 has been achieved by combining full CI on a four valence electron system in a large basis set with the usual error corrections: two-tier extrapolation using systematic basis set sequences to complete basis limits and estimation of core/valence correlation changes by a novel CI technique. Relativistic corrections are much smaller than basis set defects in the valence FCI or the core/valence correlation. While the potential is based on FCI, a MR-CISD (with the valence electrons fully exploiting the 2s-2p active
J. Phys. Chem. A, Vol. 114, No. 33, 2010 8695 space, but only two electrons excited into higher orbitals) would be equally accurate. Diffuse s,p,d basis functions are important at all bond distances, but particularly in the “outer” potential, whereas an f basis function is quantitatively important in the “inner well”. The final estimated potential is a very good match to recent experiments.18 At R ) 2.45 Å, the directly computed De ) 786 cm-1 is corrected (as described above) to 912 cm-1, compared to the experimental value 929 cm-1. The straightforward structure of the valence FCI calculation (unextrapolated and uncorrected, using an aug-cc-pVQZ basis) made an analysis possible that elucidated the origin and shape of the Be2 potential. It is found that bonding in Be2 comes about as follows. The largest contribution to the binding energy is due to hybridization of 2s and 2p orbitals in the RHF wave function,23 which lowers the RHF potential at Req by 5,500 cm-1, leaving a repulsion of only ∼2600 cm-1. While the near-degeneracy 2s2f2p2 electron correlation is large at all distances [E(MCSCF) lies 19 600 cm-1 below E(RHF) at Req], it is nearly constant with R,2,6,9 so that the MCSCF potential is almost exactly parallel to RHF and still repulsive by ∼2300 cm-1 at Req.. Net bonding is only realized by the additional dynamical correlation energy, which increases from ∼1,100 cm-1 at R ) ∞ to ∼4200 cm-1 at Req, thus yielding the ca. -800 cm-1 FCI bond energy. The crucial contributions of this dynamic correlation are those of the 10 active orbitals that are derived from the atomic 3d shells. They produce a qualitatively correct potential, namely, binding at about the correct distance as well as the peculiar slope change at R ) 3.2 Å. Tellingly, the occupation of this orbital group increases by a factor of 7 in going from R ) ∞ to Re, primarily in the πu and σg orbitals. In particular, the occupation of the 3d-πu orbitals increases very rapidly around R ) 3.2, exactly where the slope changes between “inner” and “outer” wells. Inclusion of higher dynamic correlation orbitals, following this essential 3d-type increment, converges the potential toward its quantitative shape, the largest of these contributions deriving from the atomic 4f shells. The orbitals deriving from the atomic 3s and 3p shells have rather little correlating effect, presumably because the usefulness of the correlation they can provide is diminished by preemption from the more easily accessible near-degeneracy correlation in the lower 2s-2p space. All aforementioned observations derive in some way from the fact that, in the Be atom, the 2p orbitals as well as the 2s orbital are energetically accessible to the two valence electrons. Only one part of the total picture remains difficult to understand, namely, the surprising constancy of the 2s2-2p2 near-degeneracy correlation over the entire relevant range of the internuclear distance, which allows the higher dynamic correlations, notably by the 3d orbitals, to determine the shape of the potential energy curve. Acknowledgment. The authors are grateful to David Feller for arranging access to the various new Correlation Consistent basis sets for Be, in advance of their publication.28 They thank Laimutis Bytautas for stimulating discussions and help in the literature search. The present work was supported by the DOE Chemical Physics program (M.W.S., K.R.), the NSF Petascale Applications grant (M.W.S.), and the NSF Cyberinfrastructure grant (M.W.S.). J.I. thanks Jack Collins for valuable discussions and the staff and administration of the Advanced Biomedical Computing Center for their help and support. M.W.S. (35 years) and J.I. (15 years) wish to thank K.R. for his long-standing willingness to share with them his deep insights into quantum
8696
J. Phys. Chem. A, Vol. 114, No. 33, 2010
chemistry, just one of which is that much of the explanation for the chemistry should be sought in the shapes of the orbitals. References and Notes (1) Liu, B.; McLean, A. D. J. Chem. Phys. 1980, 72, 3418–3419. (2) Lengsfield, B. H.; McLean, A. D.; Yoshimine, M.; Liu, B. J. Chem. Phys. 1983, 79, 1891–1895. (3) Petersson, G. A.; Shirley, W. A. Chem. Phys. Lett. 1989, 160, 494– 501. (4) Lepetit, M. B.; Malrieu, J. P. Chem. Phys. Lett. 1990, 169, 285– 291. (5) Noga, J.; Kutzelnigg, W.; Klopper, W. Chem. Phys. Lett. 1992, 199, 497–504. (6) Evangelisti, S.; Bendazzoli, G. L.; Gagliardi, L. Chem. Phys. 1994, 185, 47–56. (7) Røeggen, I.; Almlo¨f, J. Int. J. Quantum Chem. 1996, 60, 453–466. (8) Fusti-Molnar, L.; Szalay, P. G. Chem. Phys. Lett. 1996, 258, 400– 408. (9) Sta¨rck, J.; Meyer, W. Chem. Phys. Lett. 1996, 258, 421–426. (10) Martin, J. M. L. Chem. Phys. Lett. 1999, 303, 399–407. (11) Gdanitz, R. J. Chem. Phys. Lett. 1999, 312, 578–584. (12) Røeggen, I.; Veseth, L. Int. J. Quantum Chem. 2005, 101, 201– 210. (13) Harkless, J. A. W.; Irikura, K. K. Int. J. Quantum Chem. 2006, 106, 2372–2378. (14) Patkowski, K.; Podeszwa, R.; Szalewicz, K. J. Phys. Chem. A 2007, 111, 12822–12838. (15) Bondybey, V. E.; English, J. M. J. Chem. Phys. 1984, 80, 568– 570. (16) Bondybey, V. E. Chem. Phys. Lett. 1984, 109, 436–441. (17) Bondybey, V. E. Science 1985, 227, 125–131. (18) Merritt, J. M.; Bondybey, V. E.; Heavens, M. C. Science 2009, 324, 1548–1551. (19) Aziz, R. A.; Slaman, M. Chem. Phys. 1989, 130, 187–194. (20) Spirko, V. J. Mol. Spectrosc. 2006, 235, 268–270. (21) Bernath, P. F. Science 2009, 324, 1526–1527. (22) Patkowski, K.; Spirko, V.; Szalewicz, K. Science 2009, 326, 1382– 1384. (23) Kutzelnigg, W. Theoretical Models of Chemical Bonding; Maksic, Z. B., Ed.; Springer-Verlag: Berlin, 1990; Part 2, pp 1-43. (24) Bytautas, L.; Ruedenberg, K. J. Chem. Phys. 2008, 128, 21438/ 1–12. (25) The potential was sampled from R ) 2.0 to 3.0 in steps of 0.1, from 3.0 to 4.0 in steps of 0.2, from 4.0 to 7.0 in steps of 0.5, then 8, 9, 10, and 20 Å. A few calculations were performed at 2.45 Å, which is close to the FCI minimum and experiment. Since size-extensive methods at R ) 20 Å differ from the sum of two atoms by less than 1 microhartree, this distance was considered “separated atoms” for any method that is not size-extensive. (26) 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. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347–1363. (27) ACES II is a program product of the Quantum Theory Project, University of Florida. Authors: Stanton, J. F.; Gauss, J.; Watts, J. D.; Nooijen, M.; Oliphant, N.; Perera, S. A.; Szalay, P. G.; Lauderdale, W. J.; Kucharski, S. A.; Gwaltney, S. R.; Beck, S.; Balkova´, A.; Bernholdt, D. E.; Baeck, K. K.; Rozyczko, P.; Sekino, H.; Hober, C.; Bartlett, R. J. (28) Some previous workers have used unpublished basis sets for Be that were available on the famous PNNL basis set exchange web site. Every correlation consistent basis set used herein was taken from a forthcoming paper (Prascher, B. P.; Peterson, K. A.; Woon, D. E.; Dunning, T. H.; Wilson A. K. To be submitted for publication), which makes small adjustments to the polarization exponents and diffuse augmentations of the previously
Schmidt et al. posted Be basis sets, and defines the new “omega” basis sets for core/ valence correlation. (29) Perdew, J. P.; Ruzsinsky, A.; Csonka, G. I.; Constantin, L. A.; Sun, J. Phys. ReV. Lett. 2009, 103, 026403/1-4. (30) Two studies involving DFT-based methods for Be2 are of interest: (a) Lotrich, V. F.; Bartlett, R. J.; Grabowski, I. Chem. Phys. Lett. 2005, 405, 43–48. (b) Ruzsinszky, A.; Perdew, J. P.; Csonka, G. I. J. Phys. Chem. A 2005, 109, 11015–11021. (31) (a) Fletcher, G. D.; Schmidt, M. W.; Gordon, M. S. AdV. Chem. Phys. 1999, 110, 267–294. (b) Aikens, C. A.; Webb, S. P.; Bell, R. L.; Fletcher, G. D.; Schmidt, M. W.; Gordon, M. S. Theor. Chem. Acc. 2003, 110, 233–253. (32) Piecuch, P.; Kucharski, A. A.; Kowalski, K.; Musial, M. Comput. Phys. Commun. 2002, 149, 71–96. (33) Piecuch, P.; Wloch, M. J. Chem. Phys. 2005, 123, 214107/1–10. (34) The specific perturbative Q correction used is variant b (eq 26) in Kowalski, K.; Piecuch, P. J. Chem. Phys. 2000, 113, 5644–5652. (35) Noga, J.; Bartlett, R. J. J. Chem. Phys. 1987, 86, 7041–7050. (36) Ivanic, J. J. Chem. Phys. 2003, 119, 9364–9376; 9377-9385. (37) This particular MR-PT method is called MCQDPT when applied to multiple states, but MRMP when only a single state is treated, as here: (a) Hirao, K. Chem. Phys. Lett. 1992, 190, 374–380. (b) Nakano, H. J. Chem. Phys. 1993, 99, 7983–7992. (38) An early use of full valence MCSCF plus single and double excitations, which is sometimes referred to as “second-order CI”, is Lie, G. C.; Hinze, J.; Liu, B. J. Chem. Phys. 1973, 59, 1872–1886. (39) Bytautas, L.; Ruedenberg, K. J. Chem. Phys. 2005, 122, 154110/ 1–21. (40) (a) Bytautas, L.; Ruedenberg, K. J. Chem. Phys. 2010, 132, 074307-1/15. (b) Bytautas, L.; Matsunaga, N.; Ruedenberg, K. J. Chem. Phys. 2010, 132, 074109/1-10. (41) Bytautas, L.; Matsunaga, N.; Nagata, T.; Gordon, M. S.; Ruedenberg, K. J. Chem. Phys. 2007, 127, 204313/1–19. (42) Bytautas, L.; Ruedenberg, K. J. Chem. Phys. 2004, 121, 10905– 10918. (43) (a) Nakajima, T.; Hirao, K. J. Chem. Phys. 2000, 113, 7786–7789. (b) Wolf, A.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2002, 117, 9215– 9226. (44) Feller, D. J. Chem. Phys. 1992, 96, 6104–6114. (45) Wilson, A. K.; Dunning, T. H. J. Chem. Phys. 1997, 106, 8718– 8726. (46) Peterson, K. A.; Wilson, A. K.; Woon, D. E.; Dunning, T. H. Theor. Chem. Acc. 1997, 97, 251–259. (47) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1995, 103, 4572– 4585. (48) Pyykko, P.; Atsumi, M. Chem. Eur. J. 2009, 15, 186–197. (49) Cressy, N.; Miller, K. R.; Ruedenberg, K. Int. J. Quantum Chem. 1969, 3, 107–113. (50) Single reference CISD/aug-cc-pVQZ calculations for Ne2 at R ) 10 Å have 15.882981 valence electrons in the 2s,2p natural orbitals, scarcely different from 15.882925 at Re ) 3.091 Å. The occupation numbers of the correlating orbitals are also nearly invariant with bond distance: 3p ) 0.055107, 3d ) 0.035422, 3s ) 0.012069, and all others ) 0.014421 e- at R ) 10 Å; and 3p ) 0.055116, 3d ) 0.035441, 3s ) 0.012066, and all others ) 0.014452 e- at Re ) 3.091 Å. (51) (a) London, F. Trans. Faraday Soc. 1937, 33, 8–26. (b) Buckingham, A. D. AdV. Chem. Phys. 1967, 12, 107–142. (c) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular Forces, Their Origin and Determination; Clarendon Press: Oxford, 1981. (d) Stone, A. J. The Theory of Intermolecular Forces; Clarendon Press: Oxford, 1996. (52) LeRoy, R. J.; Huang, Y.; Jary, C. J. Chem. Phys. 2006, 125, 164310/ 1–12. (53) Porsev, S. G.; Derevianko, A. J. Exp. Theor. Phys. 2006, 102, 195–205.
JP101506T
