当前位置：搜档网 › 09 Relative Energies and Thermochemistry

09 Relative Energies and Thermochemistry

A Chemist’s Guide to Density Functional Theory. Second Edition

Wolfram Koch, Max C. Holthausen

ISBNs: 3-527-30372-3 (Softcover); 3-527-60004-3 (Electronic) 9Relative Energies and Thermochemistry

Information about the energetic properties of molecules is at the heart of every quantum chemical investigation. In this chapter we will take a close look at the performance of current approximate density functionals when it comes to the determination of such impor-tant energetic properties as atomization energies, ionization energies (IE) and electron af-finities (EA). In particular the first quantity is of great value since it provides an idea of the ballpark figure we are talking about when we are interested in computing the thermochem-istry of chemical reactions. The energies needed to remove (IE) or add an electron (EA) to an atom or molecule are obviously of significant interest in their own right. For example, the interpretation of photoelectron spectroscopy experiments is greatly facilitated if accu-rate ionization energies are available. Furthermore, since the computation of all these dif-ferent quantities poses severe and often different demands on the method chosen, the accu-racy with which a functional delivers such energies is a probe for its versatility. Somewhat counterintuitively at first sight, the reliable calculation of energetic information for atoms, in particular for transition-metals, is especially difficult and not without ambiguity. We therefore devote a complete section to this problem. Finally, we take up the discussion on excitation energies from Section 5.3.7 and give an overview of the current state of the art in the determination of electronically excited states and the corresponding transition energies using density functional theory. We conclude with a few remarks about the ability of ap-proximate density functional theory to reproduce singlet-triplet gaps in carbenes and re-lated species.

9.1Atomization Energies

Chemical reactions involve the cleavage and formation of bonds within molecules. The calculation and prediction of thermochemical data has long been a vivid field for quantum chemistry (see, e. g., Irikura and Frurip, 1998). For example, whenever experimental data is not available and empirical estimates fail, some type of quantum chemistry usually be-comes involved to obtain the missing information, but ‘computational thermochemistry’ is of great relevance also in many other areas. In practice, there is always the need to reach some compromise between accuracy and computational effort. Hartree-Fock theory pro-vides an exact treatment of exchange and scales well with the molecular size, but it suffers from severe deficiencies in describing chemical bonding due to the neglect of correlation energy contributions. Except for isodesmic (or related) reactions (see, e. g., Hehre et al., 1986) it cannot be used for thermochemical predictions. The introduction of dynamic and static correlation effects by means of post-HF wave function based methods improves the situation to a desired accuracy, but severely suffers from the notoriously unfavorable scal-ing with molecular size. Notwithstanding the rather positive appraisal of the local density approximation for the evaluation of molecular geometries and vibrational frequencies in the preceding chapter, binding energies obtained with this method are generally very inac-

137

9 Relative Energies and Thermochemistry

Table 9-1. Deviations between computed atomization energies and experiment for the JGP test set employing the 6-31G(d) basis set [kcal/mol]. Taken from Johnson, Gill and Pople, 1993.

HF MP2QCISD SVWN SLYP BVWN BLYP mean abs. dev. 86 22 2936 (40)a38 4 (4)a6

mean dev.–86–22–2936 (40)a380 (4)a1

a Basis set free results taken from Becke, 1992.

curate. The literature is full of results giving testimony to the inability of this method to deliver even qualitatively meaningful answers for problems related to chemical energetics. This problem was recognized more than two decades ago and was actually the main stimu-lus for the development of gradient corrections and, later, hybrid functionals.

Computed atomization energies, i. e., of the (hypothetical) reactions in which a mol-ecule is broken up into its constituent ground state atoms, are often very error-prone since their evaluation necessitates the breaking of each bond in a molecule. While we usually have a closed-shell molecule on the left hand side of the reaction, the ground state atoms defining the right hand side are open-shell with varying multiplicity. Hence large differen-tial correlation effects are typical for these reactions. Therefore, the calculation of atomiza-tion energies is a stringent test for any computational strategy and the deviations from experiments seen in such studies can probably be considered as upper bounds. As an in-structive example, Table 9-1 shows some error statistics for atomization energies obtained with different methods in combination with the rather small 6-31G(d) basis set for the JGP set of 32 first and second-row species.

It is apparent that the Hartree-Fock level is characterized by an enormous average devia-tion from experiment, but standard post-HF methods for including correlation effects such as MP2 and QCISD also err to an extent that renders their results completely useless for this kind of thermochemistry. We should not, however, be overly disturbed by these errors since the use of small basis sets such as 6-31G(d) is a definite ‘no-no’ for correlated wave function based quantum chemical methods if problems like atomization energies are to be addressed. It suffices to point out the general trend that these methods systematically un-derestimate the atomization energies due to an incomplete recovery of correlation effects, a reliable assessment of which requires sufficiently large and flexibly polarized basis sets.30 The errors are systematic because correlation effects are always stronger in molecular sys-tems than in their fragments (most correlation effects are roughly proportional to the number of spin-paired electrons). An insufficient recovery of electron correlation leads to a lack of stabilization of the parent molecular systems, which causes the underbinding tendency. The LDA, in turn, shows a notorious overbinding for every single molecule in the test set except Li2, and deviations from experimental atomization energies as large as 220 % occur 30The G2 extrapolation scheme – which is a prescription to extrapolate the quality of QCISD(T)/ 6-311++G(3df,2p) calculations – actually reaches chemical accuracy for the G2 test set of species, but only with an empirical correction, depending on the number of electron pairs in a molecule in order to better account for these effects.

138

9.1 Atomization Energies Table 9-2. Signed deviations [kcal/mol] between computed atomization energies (employing a 6-31G(d) basis) and experiment. Taken from Johnson, Gill, and Pople, 1993.

Molecule SVWN SLYP BVWN BLYP Molecule SVWN SLYP BVWN BLYP CH7630F247551118

CH2(3B1)21192–2O257681019

CH33128 3–2N23239–1 6

CH444404–3CO3746–51

C2H25055–9–6CN3745 3 9

C2H47071–3–4NO4453 513

C2H68685–1–6CO28299–311

for particular species like F2 (see below). Although better than the HF approximation, the LDA is certainly not a useful thermochemical tool with a mean absolute deviation of 36 kcal/ mol and it has largely been abandoned for this kind of studies. In an attempt to amend the situation by inclusion of a gradient-corrected correlation functional (SLYP) one ends up with even larger errors, which is irksome at first sight. A spectacular improvement, though, results from the application of gradient corrections to exchange (but not to correlation): the BVWN functional affords a reduced mean absolute deviation of 4 kcal/mol. Inspection of the last two columns in Table 9-1 again shows that inclusion of gradient corrections to correlation (to yield the BLYP functional) slightly decreases the overall accuracy. The com-parison with basis set free results for the same set of species taken from a paper published by Becke, 1992a, reveals no significant influence of the basis set size on the overall per-formance of the SVWN or the BVWN functional.31 Becke has shown in subsequent work that neither the choice of the PW91 parameterization of the uniform electron gas (instead of VWN) nor the addition of PW91 gradient corrections for correlation significantly changes the overall picture (Becke, 1992b). Also in this latter study, a slightly larger mean absolute deviation occurred for the G2 set of molecules upon inclusion of gradient corrections to correlation as compared to exchange-only corrections.

A closer look at the original data published by Johnson, Gill, and Pople, 1993, reveals that use of the gradient-corrected LYP correlation functional instead of VWN increases atomization energies for non-hydride species quite significantly, while those containing hydrogen atoms are reduced. This rather systematic trend is portrayed in Table 9-2 for a few example cases.

Use of the LYP correlation functional apparently reduces the overbinding proportional to the number of hydrogen atoms by about 1 kcal/mol per H from CH to CH4. On the contrary, for the nonhydride diatomics listed, the atomization energies increase by 4–5 kcal/ mol per atom upon substitution of VWN by the LYP functional. For species like C2H2, 31The largest deviations occur for different species, however. The BVWN results published by Becke were obtained in a post-LDA manner at LDA optimized geometries as opposed to the data published by Johnson, Gill, and Pople, 1993, which were computed selfconsistently at geometries corresponding to the respective level of theory. Hence, it is difficult to unambiguously pin down the origin of these differences.

139

9 Relative Energies and Thermochemistry

C2H4, and C2H6, in which both types of bonds are broken, the particular stoichiometry determines which functional performs better. The very same trends are observed for other correlation functionals (see Becke, 1992b). Such behavior is clearly unsatisfying and an indication that an only incomplete error cancellation is operative. Consequently, quite large errors can occur for unfortunate cases (ranging from –12 kcal/mol for H2O to +19 kcal/mol in the case of NO for the BLYP functional). While BP86 compares favorably with other GGA functionals for the evaluation of molecular structures and vibrational frequencies, it is defeated by BLYP when it comes to determining atomization energies. In combination with the 6-31G(d) basis set, the BLYP functional yields atomization energies for the G2 set with a mean absolute deviation 2 kcal/mol smaller than BP86 and a significantly smaller maximum deviation (Bauschlicher, 1995a). Using a large 6-311+G(3df,2p) basis, the situ-ation becomes even worse for BP86: the mean unsigned error is almost 5 kcal/mol smaller at the BLYP level and the difference in maximum errors is nearly 10 kcal/mol in favor of the latter.32 The newer PW91 correlation functional performs equally well or marginally better than LYP when used in combination with the Becke exchange term. Both show very similar mean absolute deviations for over 100 atomization energies evaluated in the com-prehensive study of Scheiner, Baker, and Andzelm, 1997 (cf. Table 9-3).

Although observed maximum deviations are sometimes substantial and evidently far from chemical accuracy, the fact that the overall errors for gradient-corrected functionals are more than five times smaller than those of the traditional wave function based methods shows nonetheless their general suitability for thermochemical studies at a modest level of computational effort. The BLYP functional in combination with small basis sets would lend itself particularly well to thermochemical studies on extended systems, where the computational demands of larger bases or correlated post-HF methods are prohibitive. However, for medium sized hydrocarbons a large underestimation of atomization energies has been observed for the BLYP functional. The BP86 functional overestimated the same atomization energies twice as much as BLYP underestimated them, so here BPW91 seems to be the GGA functional of choice with only very moderate deviations (Curtiss et al., 1997). All in all, the introduction of gradient corrections to exchange is the key to improved thermochemical data, whereas the influence of corrections to the correlation term is rather modest (inclusion of the latter, however, has important consequences for the accuracy of computed ionization energies, see below). The importance of gradient corrections to ex-change is not completely unexpected considering the fact that exchange, which is the domi-nant component of the exchange-correlation energy in Kohn-Sham theory, is largely in error in the LDA. Gunnarson and Jones, 1985, have argued convincingly that much of the exchange error inherent to the LDA stems from an improper incorporation of angular char-acteristics and nodal structures of the Kohn-Sham orbitals. Differential exchange effects are generally overestimated and errors in atomization energies are largest for molecular systems, in which substantial changes in the orbital nodal structure occur – e. g., upon bond formation from atoms resulting in occupied antibonding orbitals as in O2. Similar argu-32Remarkably, the BLYP functional approaches or sometimes surpasses the accuracy of hybrid functionals if small basis sets are used (see Table 9-5).

140

141

ments can be used to rationalize, for example, the huge overbinding of F 2 by the LDA (borrowing from an illuminating paper by Ernzerhof, Perdew, and Burke, 1997). The ex-perimental dissociation energy of the fluorine molecule into two ground state atoms (F 2

)(g 1+Σ→ 2F (2P u )) amounts to 37 kcal/mol. Hartree-Fock and the local density approxima-tion both give ridiculously wrong results: with the 6-31G(d) basis set the former yields ?34 kcal/mol (the fluorine molecule is unbound at the HF level with respect to the two constituting F atoms) while the latter gives a binding energy of 84 kcal/mol, overshooting the experimental D 0 by more than a factor of two. Since F 2 formally contains a single bond,these large errors are somewhat irritating. It is, however, well known that the lone pairs are strongly interacting in molecular F 2, which in fact is the origin of the problems for both methods. HF fails because it does not account for correlation energy and the electron pairs repel each other too strongly. If correlation is included through second or fourth order M?ller-Plesset perturbation theory, very realistic binding energies of 35 and 30 kcal/mol,respectively, result. To LDA, on the contrary, overlapping lone pairs are nothing but a higher electron density. This method overestimates the exchange stabilization brought about by these orbital interactions, which leads to an overestimation of the molecular binding en-ergy. Incorporation of explicit density gradient dependencies into the exchange terms re-pairs the shortcomings to a large extent (BVWN gives 47 kcal/mol), but still, GGA functionals do not quite reach chemical accuracy. It is clear that the gradient-corrected functionals represent a major improvement over the local density approximation and de-liver average errors which sometimes get close to our target accuracy of 2 kcal/mol. Like-wise, the data indicate that even though the various GGA functionals differ significantly in their mathematical appearance, they all perform quite similarly. However, we also note that the maximum deviations are significant and that we are still a long way from a density functional approach that is able to generally provide chemical accuracy.

In view of this situation, Becke has taken the next logical step and improved the GGA performance by admixture of exact exchange as we have already discussed to some extent in Chapter 6. His first approach, the half-and-half scheme (Becke, 1993a), did not include gradient corrections and was not much of an improvement over GGA functionals in terms of thermochemical accuracy. However, a subsequently suggested parameterized version,which included gradient corrections to exchange and correlation, gave impressively re-duced mean errors for atomization energies of the G2 set (Becke 1993b). This was the forerunner of the now widely used B3LYP hybrid (Stephens et al., 1994) which today is the most popular density functional and is implemented into most major computer codes. For details on these functionals, the reader should leaf back to Section 6.6.

Bauschlicher and Partridge, 1995, tested the B3LYP functional in combination with different basis sets on the G2 set of molecules. In combination with the 6-31G(d) basis, it yields an accuracy comparable to that of the pure BLYP functional (5.2 kcal/mol average error). This only mediocre performance improves significantly to a mean absolute error of only 2.2 kcal/mol, if the larger 6-311+G(3df,2p) basis is used, regardless of whether the geometries were obtained at this level or with the much more affordable 6-31G(d) basis set.Also, use of the aug-cc-pVTZ basis gave an improved average error of 2.3 kcal/mol, almost reaching the desired goal of chemical accuracy. However, in this study the atomization

9.1 Atomization Energies

9 Relative Energies and Thermochemistry

Table 9-3. Mean absolute deviations (MAD) from experiment [kcal/mol] for 44 atomization energies and number of results that deviate by less than 5, between 5 and 10, and over 10 kcal/mol from experiment. Taken from Martell, Goddard, and Eriksson, 1997.

Basis set MAD< 5[5..10]> 10MAD< 5[5..10]> 10

BLYP B3LYP

6-31G(d,p) 7.6211211 5.62615 3

6-311G(d,p) 6.8231110 6.82314 7 cc-pVDZ 7.32013118.5132011 cc-pVTZ 7.2201212 3.136 5 3

BP86B3P86

6-31G(d,p)12.7 9 926 6.7211211

6-311G(d,p)10.415 722 6.22412 8 cc-pVDZ 9.9131021 5.42515 4 cc-pVTZ14.2 5 831 6.9211112

BPW91B3PW91

6-31G(d,p) 7.0221111 5.62417 3

6-311G(d,p) 6.12411 9 6.92117 6 cc-pVDZ 7.22012128.4131911 cc-pVTZ 7.8191114 3.836 6 2 energy for SO2 proved to be very problematic and the results for this molecule were found to be extremely sensitive to the basis set quality. Martell, Goddard and Eriksson, 1997, studied the performance of the three commonly used GGA functionals, namely BP86, BLYP, and BPW91 together with the corresponding hybrid functionals B3P86, B3LYP, and B3PW91 on a set of 44 small first and second-row molecules. They used four different basis sets (6-31G(d,p), 6-311G(d,p), cc-pVDZ, and cc-pVTZ) in order to assess the reli-ability of predicted atomization energies with these moderately sized bases.33 Comparing the results for different basis sets and functionals it is important to firstly realize that no particular method provides results superior to all others. The authors noted a general under-estimation of atomization energies for the two hybrid functionals B3LYP and B3PW91, which contrasts with the overestimation found for the three pure GGA protocols and, to a smaller extent, for the B3P86 hybrid. Looking at the mean absolute errors for all method/ basis set combinations documented in Table 9-3, the B3LYP functional gives the highest accuracy, closely followed by B3PW91. Larger errors occur for the B3P86 hybrid, which only marginally surpasses the BPW91 and BLYP gradient-corrected functionals in terms of accuracy. The worst performance is found for the BP86 functional, for which deviations from experiment below 5 kcal/mol are the exception. The error pattern compiled in Table

33Molecular geometries, which were not reported in this study, have been obtained using the 6-31G(d,p) and cc-pVDZ basis sets. Thus, the 6-311G(d,p) and cc-pVTZ results refer to single point energy calculations only. 142

9.1 Atomization Energies Table 9-4. Average shifts in atomization energies upon basis set enlargement [kcal/mol]. Based on data taken from Martell, Goddard, and Eriksson, 1997.

BLYP B3LYP BP86B3P86BPW91B3PW91 6-31G(d,p)→6-311G(d,p)–2.6–2.6–2.6–2.0–2.5–2.4

cc-pVDZ→cc-pVTZ 4.8 5.5 5.2 5.9 5.0 5.9

9-3 reveals that the B3LYP and B3PW91 hybrid functionals give quite reliable atomization energies in combination with the cc-pVTZ basis set. The most problematic systems in this study were SO2, ClO2 and CCl, which also pose severe difficulties for traditional quantum chemical methods.

A disturbing trend in the basis set dependence is seen from the mean unsigned errors listed in Table 9-3. Reduced errors occur for the pure GGA functionals and B3P86 when improving the quality of the Pople-type basis from 6-31G(d,p) to 6-311G(d,p) – exactly what one would expect for any well-behaved quantum chemical method. Yet the opposite trend emerges for the correlation consistent basis sets when going from the polarized dou-ble-zeta cc-pVDZ to the polarized triple-zeta cc-pVTZ basis set. Better results are obtained with the smaller basis set. Only the B3LYP and B3PW91 results show the expected behavior, these two functionals deliver the smallest mean absolute errors of all methods if combined with the large cc-pVTZ basis set. These baffling findings can be rationalized by inspection of the data summarized in Table 9-4.

In spite of large differences in mean errors obtained with the various methods tested, substitution of the 6-31G(d,p) basis by the larger 6-311G(d,p) set yields a systematic shift to reduced atomization energies on average by 2.5 kcal/mol, irrespective of the method used. Conversely, use of the larger cc-pVTZ instead of the cc-pVDZ basis set brings about an increase in atomization energies, on average the order of +5.5 kcal/mol. Apparently, the 6-311G(d,p) basis yields a better description of isolated atoms, whereas the improved cor-relation consistent basis stabilizes molecular systems quite significantly with respect to the atoms. These shifts explain why GGA functionals, which usually overestimate atomization energies, perform better with the larger Pople-type 6-311G(d,p) and the smaller cc-pVDZ correlation consistent basis set. The clearly visible reason is once again error cancellation.

Another study by Scheiner, Baker, and Andzelm, 1997, has quite extensively addressed the evaluation of atomization energies with respect to different functionals and basis sets of varying quality. In this work it has been observed that polarized double- and triple-zeta basis sets which have been explicitly optimized at the LDA level (denoted DZVP LDA and TZVP LDA) are better suited for LDA and GGA calculations than for hybrid functionals. Use of a standard TZV2P basis for BLYP and BPW91 led to ca. 1 kcal/mol larger errors compared to results obtained with the TZVP LDA basis set, whereas just the opposite has been observed for the B3LYP hybrid functional. Furthermore, the latter functional still showed a remarkable drop by 2.5 kcal/mol in mean absolute deviations if a large uncontracted aug-cc-pVTZ basis set was used – in contrast, only marginal improvements (below 0.3 kcal/ mol) were seen for BLYP and BPW91. Apparently, the basis set requirements for con-

143

9 Relative Energies and Thermochemistry

verged results are higher for hybrid methods than for GGA functionals. Furthermore, it seems that hybrid methods get along much better with standard basis sets taken from the wave function ab initio world than GGA functionals do.

Redfern, Blaudeau, and Curtiss, 1997, have tested the BLYP, B3LYP, BPW91, and B3PW91 functionals with respect to the accuracy of atomization energies computed for a set of 19 molecules containing third-row, non-transition-metal elements. They used a rather large 6-311+G(3df,2p) basis set for single point energy calculations on top of MP2/6-31G(d) geometries (which might not give the highest accuracy possible for the density functional treatment). Among the functionals tested, the B3PW91 approach afforded the lowest aver-age unsigned and maximum error (2.1 and 5.7 kcal/mol, respectively) which compared nicely to much more costly G2 calculations (1.2 and 5.2 kcal/mol). The B3LYP functional gave slightly worse energetics (3.3 and 6.2 kcal/mol), whereas the pure GGA functionals led to larger mean errors and substantial maximum deviations (BLYP: 5.3 and 24.2 kcal/ mol; BPW91: 4.5 and 24.7 kcal/mol).

Before we end this discussion let us take up the thread from Section 6-9 and present some results pertaining to the large number of new functionals that have emerged in the past few years. The literature contains a variety of attempts to further improve the accuracy of density functional methods, which essentially follow two distinct lines, namely (a) the fitting of adjustable functional parameters to some kind of experimental data and (b) the fulfillment of theoretically derived and physically meaningful conditions. Examples, for instance, belonging to the first category are the CAM(A) and CAM(B) exchange functionals reported by Laming, Termath, and Handy, 1993. Two different fitting procedures to experi-mental data have been applied and for a G2 subset the resulting functionals showed a non-uniform performance. If combined with the LYP correlation functional, the CAM(A)-LYP functional gave significantly improved geometries as compared to BLYP but at the same time much worse mean errors for atomization energies. In contrast, the CAM(B)-LYP func-tional showed the reverse behavior with worse geometric parameters than CAM(A)-LYP but smaller errors for atomization energies (CAM(A)-LYP: 21.9 kcal/mol, CAM(B)-LYP: 6.5 kcal/mol, BLYP: 9.5 kcal/mol). Such a situation is of course by no means satisfying. On the one hand, the results show that there is definitely room for improvement within the particular formulations of GGA functionals by means of fitting procedures and CAM(B)-LYP might indeed appear as a useful improvement over the BLYP functional. On the other hand, one could expect from theoretical reasoning that an improved description of molecu-lar binding also leads naturally to a better performance in structure prediction, which obvi-ously is not the case. It is hence apparent that reparameterization does not necessarily im-prove the fundamental physics but rather exerts some shift on the outcome of error com-pensation effects. Therefore, one can rightly argue that this might lead to a non-systematic performance for molecules or properties not included in the fitting set, which would render the quality of such corrections very difficult to judge a priori. Obviously, only extensive testing can identify the particular advantages and potential caveats of such functionals.

The idea of modifying existing functionals by fitting particular terms to accurate experi-mental data has been tempting to others, too. Stewart and Gill, 1995, have reparameterized a simplified LYP correlation functional formalism, and tested its performance for atomiza-144

9.1 Atomization Energies tion energies in combination with Becke’s 1988 (B) exchange term. This new simple func-tional, referred to as Becke-Wigner (BW), has been tested with the 6-31G(d) and 6-311+G(3df,2p) basis sets to evaluate the atomization energies of the G2 set. The result-ing mean absolute errors were slightly in favor of BLYP (6-31G(d): 5.5 and 4.5, 6-311+G(3df,2p): 4.9 and 4.7 kcal/mol, for BW and BLYP, respectively). However, the computed data for individual molecules were quite different. If combined with the rather modestly sized 6-31+G(d) basis set, the empirical density functional EDF1 by Adamson, Gill, and Pople, 1998, yields reasonably good results with mean absolute errors for atomi-zation energies of 3.2 kcal/mol compared to 4.4 and 5.9 kcal/mol for BLYP and B3LYP in this basis, respectively. No improvement has been found upon exact exchange admixture. It will be interesting to see a further assessment of the accuracy of this functional in future applications. Neumann and Handy, 1995, implemented the Becke-Roussel exchange func-tional (BR), which was fitted to model the shape of the Hartree-Fock exchange hole in a two-term Taylor expansion without any reference to the electron gas model. This func-tional (which depends on the density, its gradient and Laplacian as well as on the kinetic energy density) was tested in combination with the P86 correlation functional, employing

a polarized triple-zeta TZ2P basis set, and atomization energies were obtained for a set of

27 diatomic first and second-row molecules with a mean error of 5.5 kcal/mol (compared to 6.0 kcal/mol for BP86 applied to the same set). When used in a refitted three parameter hybrid framework with some 20 % exact exchange admixture (BR3P86), the resulting at-omization energies were improved with respect to the pure GGA, but slightly larger overall and maximum deviations occurred as compared with B3P86 results. As briefly mentioned in Chapter 6, Becke, 1997, proposed another functional, which contains exact exchange admixture and was derived from a systematic fitting to thermochemical data of the G2 set by adjusting 10 parameters, i. e., B97. Atomization energies were obtained with a mean absolute deviation of 1.8 kcal/mol and an absolute maximum error of 5.5 kcal/mol. This accuracy closely approaches that of the G2 extrapolation scheme, for which 1.2 and 5.1 kcal/ mol result for mean and maximum absolute deviations, respectively. Hence, this fitting scheme created a functional which definitely surpasses the quality of hybrid functionals like B3LYP or B3PW91. However, this new method is awaiting the extensive testing which the latter two hybrids have seen in the recent past, and it remains to be verified whether its performance endures as favorably as found for the test set for which it has been parameterized. B97-1, the self-consistent reparameterization of the B97 functional leads to a slightly larger mean absolute error (2.2 kcal/mol) for atomization energies of a G2 subset of species if combined with a TZ2P basis set (Hamprecht et al., 1998). Rather impressive thermochemical results as documented in Table 9-5 have been reported by van V oorhis and Scuseria, 1998, for their VSXC functional. This functional was the outcome of a fitting procedure adjusting no less than 21 different parameters. In addition, it goes beyond the standard GGA functionals by depending also on the non-interacting kinetic energy density. Further developments along similar lines have been reported in the recent literature and are discussed in Chap-ter6.

Neumann and Handy, 1996, implemented the recent B95 correlation functional and tested it on a G2 subset. This functional was originally proposed by Becke and obeys

145

9 Relative Energies and Thermochemistry

some physically motivated minimal requirements, thus representing our first example of path (b) among the lines of modern functional development. In Becke’s original work (Becke, 1996a) the new method has been applied in a post-LDA manner (i. e., the func-tional was applied on KS orbitals and the corresponding density obtained from a con-verged SVWN calculation – as usually done by this author) whereas Neumann and Handy tested a fully selfconsistent implementation. Benchmarked for atomization energies of the G2 set (or a subset thereof by the latter authors), the new functional led to large overbinding for non-hydrogen species, exaggerating the above mentioned observations for the inclu-sion of gradient corrections to correlation even more. Species containing hydrogen atoms on the other hand, were described with a better accuracy. The error of the pure GGA (i. e., non-hybrid) form (BB95, 8.8 kcal/mol) was found inferior even to BP86 (6.0 kcal/mol) by Neumann and Handy, confirming the disappointing results of Becke’s original investi-gation. Better results, however, were obtained for a fitted single parameter hybrid imple-mentation, blending B with exact exchange (dubbed B1B95). Becke found this functional superior to his three parameter fit (mean unsigned error 2.0 vs. 2.4 kcal/mol). Correspond-ingly, a smaller error was also reported for B1B95 (2.6 kcal/mol) than for B3P86 (3.2 kcal/ mol) in the study by Neumann and Handy. Another functional, which does not rely on empirical adjustments, is the PBE functional introduced by Perdew, Burke, and Ernzerhof, 1996. When applied to atomization energies for the G2 set of species (in combination with the rather flexible 6-311+G(3df,2p) basis set on top of MP2 optimized structures), this functional performs much better than SVWN, but does not reach the accuracy of BLYP (mean absolute errors are SVWN: 36.4, PBE: 8.6, BLYP: 4.7 kcal/mol). Admixing of 25 % exact exchange does ameliorate the performance but, as reported by Ernzerhof and Scuseria, 1999a, the resulting PBE1PBE hybrid functional still falls short of the B3LYP hybrid. For the extended G2 set (148 molecules) the errors amount to 4.8 and 3.1 kcal/mol for the PBE1PBE and B3LYP functionals, respectively. For the original G2 set consisting of 55 molecules the errors are reduced to 3.5 kcal/mol (PBE1PBE) and 2.4 kcal/mol (B3LYP). In a related study, Adamo and Barone, 1999, report an absolute mean error for the PBE1PBE functional combined with the 6-311++G(3df,3pd) basis set on the 32 mol-ecule JGP subset of the G2 database of 2.6 kcal/mol. Finally we mention the recent contri-bution by Rabuck and Scuseria, 1999, who applied the B3LYP, VSXC, PBE1PBE and PBE functionals to the determination of enthalpies of formation for molecules which are not included in the typical density functional training sets and which are known to be problematic. As expected, the average errors are significantly larger than for the G2 or related references. The best performance is achieved with the VSXC functional (8.8 kcal/ mol absolute average deviation, 24.3 kcal/mol maximum deviation) followed by B3LYP (10.6 and –37.2 kcal/mol, respectively) and PBE1PBE (11.5 and 39.5 kcal/mol, respec-tively). The pure GGA functional PBE works significantly worse and shows an average error of 38.2 kcal/mol and a maximum error exceeding 100 kcal/mol, rendering it fairly useless in this context. Remember that VSXC achieves its good performance without con-taining any Hartree-Fock exchange. Rather, it differs from regular GGA functionals by the fact that it depends not only on the density gradient but also on the kinetic energy density.

146

9.1 Atomization Energies

In conclusion, according to the results of a variety of systematic studies, the introduction of hybrid functionals can be considered a successful step towards the ultimate goal of chemi-cal accuracy for the evaluation of atomization energies of main group species, provided that basis sets of polarized triple-zeta quality or better are used. Although functionals like B3LYP and B3PW91 do not quite reach the target accuracy of below 2 kcal/mol, they provide a pragmatic means to predict atomization energies with a pleasing accuracy. As such, they constitute highly efficient alternatives to far more demanding post-HF methods, which show comparable mean and maximum deviations in a variety of cases. These two hybrid methods in particular are available in several major computer codes and provide a significant improvement over results for pure GGA functionals with only few exceptions. From the data given above, the rough hierarchy of functionals given in Section 6-9, i. e., LDA < GGA < hybrid functionals, is confirmed. If we go one step further and also ask the question, which of the widely available functionals is to be recommended with respect to the quality of the computational prediction of atomization energies, we arrive at the follow-ing conclusion (with the expected accuracy increasing from left to right):

SVWN << BP86 < BLYP ≈ BPW91 < B3P86 < B3LYP ≈ B3PW91.

In terms of basis sets, there is compelling evidence that sets smaller than polarized tri-ple-zeta quality significantly reduce the accuracy that can be obtained with modern hybrid functionals and cannot be recommended if quantitative energetic results are the prime tar-get.

In Table 9-5 we summarize the performance of various functionals discussed above as collected from many sources, which highlights the conclusions of the above discussion. Table 9-5. Compilation of mean absolute and maximum absolute deviations (in parentheses) for atomization energies [kcal/mol] of small main group molecules from different sources.

32 1st row species, 6-31G(d) basis set, Johnson, Gill, and Pople, 1993

HF85.9SVWN35.7

MP222.4BVWN 4.4

QCISD28.8BLYP 5.6

33 1st and 2nd row diatomic molecules, TZ2P basis, Laming, Termath, and Handy, 1993

LDA43.6(18.3)CAM(A)-LYP21.9(14.5) BLYP9.5 (9.3)CAM(B)-LYP 6.5(12.0) G2 set, B1 = 6-31G(d), B2 = 6-311+G(3df,2p), Bauschlicher, 1995

HF/B180.5(184.3)HF/B274.5(170.0) MP2/B116.0 (40.3)MP2/B2 7.3 (25.4) BLYP/B1 5.3 (18.8)BLYP/B2 5.0 (15.8) BP86/B1 7.2 (24.0)BP86/B210.3 (25.4) B3LYP/B1 5.2 (31.5)BP86/B2 2.2 (8.4) B3P86/B1 5.9 (22.6)BP86/B2 7.8 (22.7)

147

9 Relative Energies and Thermochemistry

Table 9-5, continued.

G2 set, Bauschlicher and Partridge, 1995

6-31G(d) aug-cc-pVTZ6-311+G(3df,2p)

B3LYP 5.2(31.5) 2.6(18.2) 2.2(8.1)

44 1st and 2nd row species, Martell, Goddard, and Eriksson, 1997

BLYP BP86BPW91B3LYP B3P86B3PW91 6-31G(d,p)7.612.77.0 5.6 6.7 5.6

6-311G(d,p) 6.810.4 6.1 6.8 6.2 6.9

cc-pVDZ7.3 9.97.28.5 5.48.4

cc-pVTZ7.214.27.8 3.1 6.9 3.8

G2 set, Stewart and Gill, 1995

BW/6-31G(d) 5.5(25.5)BLYP/6-31G(d) 4.5(16.3) BW/6-311+G(3df,2p) 4.9(15.3)BLYP/6-311+G(3df,2p) 4.7(15.3)

27 1st and 2nd row diatomic molecules, TZ2P basis, Neumann and Handy, 1995, 1996

BP86 6.0(18.3)B3P86 3.2 (9.3) BRP86 5.5(14.5)BR3P86 3.1(12.0)

B1B958.8(24.1)B1B95 2.6 (9.4)

19 species incl. 3rd row atoms, 6-311+G(3df,2p) basis, Redfern, Blaudeau, and Curtiss, 1997

G2 1.2 (5.2)

BLYP 5.3(24.2)B3LYP 3.3(6.2) BPW91 4.5(24.7)B3PW91 2.1(5.7) 108 1st and 2nd row species, Scheiner, Baker, and Andzelm, 1997

DZVP LDA TZVP LDA DZP6-31G(d)TZ2P UCC SVWN47.652.147.052.250.156.4 BLYP 7.4 6.910.2 7.0 7.4 7.1 BPW91 6.4 6.2 9.7 7.4 7.3 7.0

B3LYP 8.8 7.810.1 6.8 6.5 4.1

G2 set, 6-31+G(d) basis set, Adamson, Gill, and Pople, 1998

EDF1 3.2(15.3)BLYP 4.4(16.3)

B3LYP 5.9(35.9)

G2 (first 2 cols.) and ext. G2 set, 6-311+G(3df,2p) basis, Ernzerhof and Scuseria, 1999a SVWN36.4(84)83.7(216) PBE 8.6(26)17.1 (52) BLYP 4.7(15) 7.1 (28)

B3LYP 2.4(10) 3.1 (20) PBE1PBE 3.5(10) 4.8 (24) VSXC 2.5(10) 2.7 (8) 148

149

9.2Atomic Energies

Now that we have considered the performance of DFT for the prediction of atomization energies for main group species in some detail, we focus a little closer on the right hand side of such reactions: the atoms. Atoms are not only the smallest subunits in chemistry,they are – seemingly paradoxically – also among the most difficult systems to describe for approximate density functional theory. The only exceptions which are completely unproblematic include closed-shell atoms, such as the ground states of the rare gases, but these are not the subject of this section. Although the errors seen in the preceding section stem at least to some extent from problems describing main group atoms in general and atomic states in particular, the main thrust of the following discussion will be geared to-wards transition-metal atoms and ions. The multifaceted chemistry of transition-metals is largely determined by their variable occupation of nd, (n+1)s, and (n+1)p valence orbitals which poses severe challenges for a theoretical treatment. From a physical point of view,subtle differential correlation and exchange effects of the various nd p (n+1)s q occupations are realized in the different atomic states. A method which aims at an accurate description of atomic states must be capable of providing a balanced and unbiased representation of the many possible electronic situations. This is anything but an easy task for any current quan-tum chemical strategy, including sophisticated approaches such as configuration interac-tion or coupled cluster methods. One should therefore not be surprised that problems arise with Kohn-Sham methods based on approximate density functionals.

A second major reason why atoms are so difficult, in particular for methods rooted in approximate density functional theory, has been touched upon already in Chapter 5. In Kohn-Sham theory, by definition, we do not have access to the correct many-electron wave function and its symmetry requirements. It is therefore not clear how to deal with atomic

terms whose wave functions are eigenfunctions of the 22??L

,S and related operators (see Section 5.3.7). The usual way out is to select a single-determinantal non-interacting Kohn-Sham reference system for defining the values of the conserved quantum numbers. This leads to ambiguities and possible inconsistencies in the description of these states. Con-sider the high spherical symmetry of atomic species and recall from Chapter 5 the inability of current approximate functionals to properly account for the related degeneracy effects which occur in open-shell situations. A comprehensive computational study to investigate these problems has been reported by Baerends, Branchadell, and Sodupe 1997. They dem-onstrate that such difficulties already show up for main group atoms with partially occu-pied p-orbitals. Let us consider the example of a 2P ground state for a boron atom with its [He] (2s)2 (2p)1 electron configuration. The energy differences between the spherical den-sity with 1/3 of an electron in each of the three real p-orbitals and a non-spherical density derived from occupying the real p z orbital 34 amounts to some non-negligible 0.2 eV if the BP86 protocol is used. If instead one of the complex p-orbitals is occupied (e. g., p x + i p y which corresponds to M L = 1), the resulting energy is roughly in-between the previous two

This corresponds to the component of 2P with M L = 0. Occupying the real p x or p y orbital results in the same

energy – but note that the real p x and p y orbitals are no longer eigenfunctions of the 2?L

o perator.9.2 Atomic Energies

150

9 Relative Energies and Thermochemistry

results. The exact energy density functional would be invariant over the set of charge den-sities belonging to a degenerate ground state and would produce precisely the same energy for all these possible representations. However, none of the currently known approximate density functionals is able to meet this requirement. This type of problem is particularly prominent when it comes to describing transition-metal atoms or ions with partially filled d-shells. Here, the energy even depends on which of the real d-orbitals are selected for specifying the configuration, because the shape of the d z 2 orbital differs from the shape of the other four orbitals of this set. We already illustrated this problem in Chapter 5 for the d 1configuration of the ten-fold degenerate 2D ground state of the scandium dication. The important take home message here is that due to the deficiencies of the currently used density functionals, there is no unambiguous reference energy for atoms in approximate density functional theory .

What is obviously needed is a generally accepted recipe for how atomic states should be dealt with in approximate density functional theory and, indeed, a few empirical rules have been established in the past. Most importantly, due to the many ways atomic energies can be obtained, one should always explicitly specify how the calculations were performed to ensure reproducibility. From a technical point of view (after considerable discussions in the past among physicists) there is now a general consensus that open-shell atomic calcu-lations should employ spin polarized densities, i. e. densities where not necessarily )r (21

)r (2

1)r (H H H βαρ+ρ=ρ. Note that this does not mean that the unrestricted Kohn-Sham

formalism has to be used, restricted open-shell variants are in principle equally eligible (but recall the discussion in Section 5.3.5). All this condition states is that the α-density does not have to be equal to the β-density. Actually, this rule must seem trivial and enforces itself almost automatically, because spin unpolarized open-shell calculations are – if pos-sible at all – usually difficult to perform with most current program packages. By the same token, densities that are allowed to be non-spherical should be used. The corresponding atomic orbitals should be occupied either by one or two electrons rather than distributing the N electrons equally over the n degenerate orbitals.35 This rule is also – in principle –automatically obeyed in most calculations done with standard codes, since it represents the default way of performing such calculations. However, even in cases where a calculation is started with an integer occupation of d-orbitals, unphysical mixings between d-orbitals and the (n+1)s-orbital may occur, depending on the symmetry imposed. The solutions resulting from such scrambling of the original occupation pattern cannot usually be related to any physical state anymore, as outlined further below. When spin-polarized, non-spheri-cal densities are allowed, the additional variational freedom leads to solutions which are usually significantly lower in energy than if these restrictions are enforced.

The physical reasoning for why these densities were frequently employed in the earlier days of density func-tional theory was that in this way the degeneracy of the partially filled d-orbitals could be retained. A technical reason why these densities still have to be employed in some recent investigations is that calculations with integral orbital occupations simply do not converge in the self consistent field procedure (see, e. g., Blanchet,Duarte, and Salahub, 1997). Such densities correspond to a representation of a particular state 2S+1L with M s = S and a spherical averaging over M L .

151

Following the discussion in Section 5.3.7, among the possible occupations of the real atomic orbitals connected to a formal configuration one should select only those which correspond to a single-determinantal representation of the desired atomic term. As stated very early on by Ziegler, Rauk, and Baerends, 1977, only such single-determinantal states are valid for a description using the current Kohn-Sham technology and the corresponding approximate density functionals. This also means that there are states of atoms or mol-ecules which cannot be computed directly owing to their inherent multi-determinantal char-acter. In these cases alternative routes such as the sum method introduced in Section 5.3.7must be used. For the transition-metal atoms and their positive ions, only the lowest multiplet components of a particular configuration are needed for the ground and first excited states and no such complications occur. Rather, all these terms can be represented by single deter-minants. However, with the exception of m l = 0 all the 2l +1 components of atomic orbitals for a given l are complex and therefore not directly accessible for a representation using real orbitals. Instead, linear combinations of the complex determinants sharing the same ±m l need to be formed such that they lead to real representations. The adequate occupations of real d-orbitals, which correspond to single Slater determinants with the correct angular momentum and spin symmetry, have been summarized in an appendix of a frequently quoted paper by Hay, 1977, and are reproduced in Table 9-6.

It has been noted in Hay’s paper that the occupations for the d 1, d 4, d 6, and d 9 states are in principle arbitrary. This does not strictly hold true for density functional applications because of the above-mentioned dependence of the energy on the shape of the occupied orbitals. The density generated from occupying the d z 2 differs from the one obtained from placing the electron in, e. g., the d xy orbital. Feeding an approximate density functional with these two unequal densities may lead to non-identical energies (cf. Figure 5-2). In most practical applications, however, the errors introduced in this way should be much smaller than those caused by other limitations of the functional or basis set employed.

Table 9-6. Open-shell d-configurations after Hay, 1977.Configuration Ground State Term

Occupation

d 12D

z )d (2d 23F 1y x 1z )d ()d (222?d 34F 1

yz 1xz 1xy )d ()d ()d (d 45D

1yz 1xz 1xy 1y x )d ()d ()d ()d (22?d 56S 1yz 1xz 1xy 1z 1y x )d ()d ()d ()d ()d (222?d 65D

1yz 1xz 1xy 1y x 2z )d ()d ()d ()d ()d (222?d 74F 1yz 1xz 1xy 2y x 2z )d ()d ()d ()d ()d (222?d 83F 2yz 2xz 2xy 1y x 1z )d ()d ()d ()d ()d (222?d 9

yz 2xz 2xy 2y x 1z )d ()d ()d ()d ()d (222?9.2 Atomic Energies

9 Relative Energies and Thermochemistry

The generation of a clean occupation pattern for atomic d-orbitals is greatly facilitated if point-group symmetries can be exploited which prevent unphysical mixing. This has been discussed in some detail in a paper on density functional atomic calculations (using a re-stricted open-shell rather than the usual unrestricted strategy) of 3d transition-metal atoms by Russo, Martin and Hay, 1994. What one needs are symmetry constraints which allow the occupation of the nd- and (n+1)s-orbitals in such a way that any unwanted mixing leading to unclear, non-integer occupations is prevented by symmetry. Let us elucidate this by using the 5F term of a d3s1 occupation, as realized in the first excited state of the V+ ion as an example. If we are fortunate enough to have a program which supports non-Abelian point-group symmetries at our disposal, octahedral symmetry with inversion (point-group O h) should be employed. The three singly occupied d-orbitals are – following Table 9-6 –chosen as the d xy, d xz, and d yz orbitals which span the three-dimensional irreducible repre-sentation t2g. Hence, these three orbitals are equivalent and each is occupied by one elec-tron of the same spin. Since the Pauli principle excludes the occupation of each spin orbital by more than one electron and because the other d-orbitals belong to a different irreducible representation of O h (namely e g) the electrons deposited in the t2g orbitals cannot move to any another d-orbital. The s-orbital belongs to the a1g irreducible representation. Also since none of the d-orbitals transforms as a1g and in particular the occupied d xy, d xz, and d yz orbitals are separated from the s-orbital by symmetry, also s/d mixing is impossible. As a consequence the originally defined orbital assignment of electrons is frozen by symmetry and will not change during the course of the calculation. If instead we have only Abelian point-groups at our disposal and select C2v symmetry (the typical choice) it is still possible to unambiguously assign the three d-electrons to the d xy, d xz, and d yz orbitals. In this case, these three d-orbitals are all in an irreducible representation of their own, namely a2, b1 and b2, respectively. Thus, no mixing between these or other d-orbitals can occur. However, the formally singly occupied 4s orbital belongs to the totally symmetric a1 representation which happens to be the same irreducible representation in which the d z2 and d x2-y2 orbitals can be found. Hence, these three orbitals can now mix and there is no guarantee that the electron which was initially assigned to the 4s orbital stays there and does not partially move into the d z2 or the d x2-y2 orbital. Unfortunately, the use of high point-group symmetries and other symmetry arguments is not a panacea. First, many programs, such as for example the cur-rent versions of Gaussian simply do not support non-Abelian point-groups and O h is thus out of reach. Second, there is also a number of cases where symmetry alone, even if symmetries such as O h were accessible, does not help. A case in point is provided by the 3d14s1 occupation of the 3D ground term of the scandium cation. No point-group is avail-able which could both exclude mixing of the s- and d-orbitals and still prevent unpaired d-electrons from moving between degenerate symmetry-equivalent d-orbitals. Neither the high O h symmetry nor the Abelian C2v point-group (nor any other point-group) assures that this electron distribution will persist throughout the calculation. Let us be specific: in O h we can prevent the 4s-orbital from mixing with any of the 3d-orbitals because they are separated by symmetry. But because the 3d-orbitals transform as two- (e g) or three-dimen-sional (t2g) irreducible representations of O h, the smearing of the one d-electron between symmetry equivalent orbitals cannot be ruled out. For example, if we initially place the 152

9.2 Atomic Energies electron in the d z2 orbital, there is no symmetry related reason why this electron should not partially occupy the other d-orbital transforming as e g, i. e.,d x2-y2. Fractional occupations of these orbitals would be the result. On the other hand, while a unique assignment of the 3d electron is possible in C2v because there are three equivalent one-dimensional irreduc-ible representations of C2v in which this single electron could be uniquely accommodated, the 4s/3d mixing cannot be prevented in that point-group for the same reason as explained above. In such cases, one usually prefers C2v, but the user has to observe carefully the progress of the calculation in order to ensure that the initially adjusted occupation pattern persists until convergence of the self consistent field, at least as much as possible. These guidelines can be summarized in the following rule of thumb for the calculation of atomic energies: accept only solutions with the correct occupation pattern and with integer d-orbital occupations (see Ricca and Bauschlicher, 1995a). Finally, we should note that by lifting all restrictions with respect to symmetry and occupation pattern, a further lowering of the atomic energy can sometimes be achieved, see Baerends, Branchadell, and Sodupe 1997. The physical meaning of such solutions, which are often characterized by fractional occupations of d- and s-orbitals is, however, questionable.

In addition to these more technical problems, there are other inconsistencies which re-strict the quality of atomic energies. The most prominent issue in this context is the bias of current approximate density functionals towards preferentially occupying d- rather than s-orbitals (for detailed discussions see Gunnarsson and Jones, 1985, Ziegler and Li, 1994, Holthausen et al., 1995). This is just the opposite of what is generally seen for traditional wave function based approaches, which favor states with fewer d-electrons. We can ration-alize the contrasting shortcomings of these two different schools of theory by the following considerations exemplified for the atomic excitation energies of 3d transition-metal cati-ons. The separations between the d n versus d n–1s1 states are determined by Coulomb repul-sion, exchange energy, and electron correlation. As a consequence of the different sizes of d- and s-orbitals36 the average interelectronic distances are larger for d n–1s1 configurations, resulting in a reduced Coulomb repulsion. From the traditional wave function based view-point, this means that the correlation problem is less severe for this occupation pattern. Hartree-Fock theory, which treats Coulomb and exchange interactions exactly but neglects correlation of electrons with antiparallel spin, indeed shows a pronounced bias towards d n?1s1 configurations. This is particularly so for the late transition-metals for which this change of the electronic configuration is accompanied by a spin flip and HF favors the high-spin states. The very same situation causes a different problem if density functional theory is invoked. Due to the more compact electron arrangements in the valence d-shells, the average exchange stabilization per d-electron pair, K dd, is stronger than that between an s- and a d-electron, K sd. As to the current situation, it has been shown that the LDA overes-timates the absolute exchange terms K dd, K sd, as well as those between s- and d-electrons and core-orbitals. As discussed in a key paper by Gunnarsson and Jones, 1985, this exag-geration of exchange stabilizations is more pronounced for d n than for d n–1s1 situations 36d-orbitals are generally more compact than s-orbitals of the following main quantum number by a factor of about 1.5 to 3.4. The / ratio is particularly large for the 3d-series of elements, see Bauschlicher, 1998.

153

9 Relative Energies and Thermochemistry

which leads to the bias mentioned above. Fortunately, these effects are compensated for in

part by account of correlation energies, which operate in the opposite direction, so that

even the LDA usually performs better than the HF approach. Absolute exchange energies,

however, are much larger in transition-metal atoms than absolute correlation energies. For

example, a SVWN treatment of the Cu+ ion yields E x in the order of –60 E h, whereas E c is

about –3 E h, which illustrates the predominant influence of shortcomings in the exchange

part of functionals on observed errors.

We can already anticipate that gradient corrections to exchange or an admixture of HF

exchange in the hybrid functional scheme can help to improve agreement with experiment.

This has indeed been observed, but a tendency to artificially stabilize d n over d n–1s1 con-

figurations remains. The inclusion of atomic state splittings of transition-metals into the

databases used to construct new functionals appears to be a logical consequence, but as far

as we are aware, this has not been done as of yet. As a consequence, even using the most

advanced functionals, there are a number of cases where either wrong atomic ground states

are predicted by density functional theory, or where s/d mixing results in intermediate

occupations, which cannot be connected to physically reasonable configurations. For ex-

ample, the atomic ground state term of cobalt is 4F, dominated by a 3d74s2 occupation. The

second 4F term with a 3d84s1 occupation is – after correcting for the differential relativistic

effects – 0.17 eV higher in energy. Density functional calculations of various flavor give

the reverse result with the 3d84s1 occupation being lower in energy than the 3d74s2 one.

Finally we need to mention that heavier elements exhibit strong relativistic effects, which

also have a significant influence on the physical properties of the d-block elements. While

common wisdom has it that relativistic contributions have even qualitative consequences

for bonding or electronic state splittings of 5d transition-metals,37 their influence is not that

dramatic, but still non-negligible, for the 4d elements. The relativistic contributions to 3d

elements are often ignored, but yet, for the later elements of this row, they are larger than

inexperienced newcomers to the field might anticipate. In particular for Cu there is quite a

deal of influence on the stability of atomic states. The experimental value for the 3d10(1S)→ 3d94s1(3D) state splitting in the copper cation is 2.81 eV. However, a sophisticated post-HF treatment of these two states gives a value of 3.11 eV.38 Such a deviation is outside the

expected error range for a high quality level of theory, and indeed, the inclusion of scalar

relativistic effects by means of perturbation theory (mass-velocity and Darwin corrections)

gives a value of 2.85 eV in good agreement with experiment. Hence, the 3d94s1 configura-

tion is stabilized by 0.26 eV with respect to 3d10(1S) due to relativistic contributions. This

is certainly a non-negligible source of error for non-relativistic calculations if triplet copper

cations are involved. Hence, while non-relativistic calculations on 3d elements are very

well acceptable for many purposes, one has to take into account relativity if a higher accu-

racy is aimed at, at least for the later elements in the row. The most convenient way to

37In fact, for 5d transition-metals relativistic contributions, and in particular spin-orbit coupling, can be of the same order of magnitude as chemical bonding.

38At the highly correlated CASSCF-AQCC level of conventional ab initio theory using a very large [7s6p4d3f2g]

ANO basis set.

154

9.2 Atomic Energies include scalar relativity to some extent into the calculations is the use of relativistic effec-tive core potentials or explicit one-component schemes, but the neglect of spin-orbit cou-pling effects can be problematic for the heavier elements.

A typical set of results for excitation energies of 3d transition-metal atoms and their cations based on the results of Koch and Hertwig, 1998, is summarized in Tables 9-7 and 9-8. The SVWN, BLYP and B3LYP functionals were combined with a sufficiently flexible contracted GTO basis set of 8s7p4d2f quality to expand the Kohn-Sham orbitals. In these calculations the rules outlined above were followed. For clarity the occupations used in the respective point-group symmetries are also included in the Tables. Since the computationally predicted excitation energies have been obtained in a completely non-relativistic scheme, they are compared to experimental energies that have been empirically corrected for differ-ential scalar relativistic effects taken from Raghavachari and Trucks, 1989a and 1989b. These corrections are based on approximate calculations and neglect the influence of rela-tivity on the correlation energy (and vice versa), but have been shown to provide a good approximation. The preference of current approximate density functionals for d-rich occu-pations – in particular with the LDA – can easily be inferred from these results. Use of the gradient-corrected BLYP protocol or the hybrid B3LYP approach does lead to a significant reduction of the deviations but also to a less systematic behavior. Nevertheless, as com-pared to other strategies the overall accuracy of these results, in particular for the B3LYP functional, is satisfying. For example, the mean absolute deviations for the neutral excita-tion energies (i. e.,d n s2 versus d n+1s1 configurations) as determined with the Hartree-Fock model, second order M?ller-Plesset perturbation theory or the QCISD(T) model amount to 0.86, 0.55 eV, and 0.14 eV, respectively (Raghavachari and Trucks, 1989a). The corre-sponding density functional results are 0.74eV for the SVWN, 0.55 eV for the BLYP functional, and 0.33 eV if the B3LYP scheme is used. For the cations (i. e.,d n s1 versus d n+1) the density functional approaches perform significantly better with mean unsigned errors of 0.32, 0.18 and only 0.16 eV for LDA, BLYP, and B3LYP, respectively. For comparison, Raghavachari and Trucks, 1989b report mean deviations of the HF, MP2 and QCISD(T) schemes of 1.32, 0.35, and 0.23 eV, respectively. It should be noted that due to the above mentioned limitations in the symmetry treatment for some of the difficult cases, the con-verged wave functions were not absolutely clean. A case in point is provided by the Fe 5D ground term. The corresponding d6s2 occupation cannot be treated in a unique, point-group symmetry determined way. As indicated in Table 9-7, one has to resort to C2v symmetry and face the problem of s/d mixing. In fact, in all final Kohn-Sham wave functions some of the s-electron population had moved to the symmetry related d-orbitals, creating a slightly fuzzy picture and adding to the inherent uncertainty of the results. Overall, this uncertainty in atomic Kohn-Sham calculations is at least in the order of some tenths of an eV, but larger deviations may also occur using even state-of-the-art functionals.

155

156

9 Relative Energies and Thermochemistry

T a b l e 9-7. E x c i t a t i o n E n e r g i e s [e V ] o f 3d -T r a n s i t i o n M e t a l s .

A t o m n

d n s 2

P o i n t g r o u p

O c c u p a t i o n

d n +1s 1P o i n t g r o u p

O c c u p a t i o n L S D B L Y P B 3L Y P E x p .a

S c

C 2v

a 11a 124

O h

e g 2a 1g 10.560.670.791.33

T i

O h

e g 2a 1g 2

O h

t 2g 3a 1g 1

–0.260.060.240.69

O h

t 2g 3a 1g 2

C 2v

a 11

b 11b 21a 21a 11

–0.94–0.39–0.19

0.11

C r

C 2v

a 11

b 11b 21a 21a 12

O h

e g 2 t 2g 3a 1g 1

–1.15–0.81–0.79

–1.17

M n

O h

e g 2 t 2g 3a 1g 2

C 2v

a 12a 11

b 11b 21a 21a 11

1.001.001.42

1.97

F e

C 2v

a 12a 11

b 11b 21a 21a 12

O h

e g 4t 2g 3a 1g 1

0.220.210.28

0.65

C o

O h

e g 4t 2g 3a 1g 2

O h

e g 2t 2g 6a 1g 1

–0.72–0.47

–0.10

0.17

N i 8

O h

e g 2t 2g 6a 1g 2

C 2v

a 12a 11

b 12b 22a 22a 11

–1.26–0.80

–0.41

–0.33

C u 9

C 2v

a 12a 11

b 12b 22a 22a 122

O h e g 4t 2g 6a 1g 1

–2.57–2.14

–1.82

–1.85

C o r r e c t e d f o r r e l a t i v i s t i c e f f e c t s , t a k e n f r o m R a g h a v a c h a r i a n d T r u c k s , 1989a .

T a b l e 9-8. E x c i t a t i o n E n e r g i e s [e V ] o f 3d -T r a n s i t i o n M e t a l C a t i o n s .

C a t i o n n

d n +1 s 1P o i n t g r o u p O c c u p a t i o n

d n +2

P o i n t g r o u p O c c u p a t i o n L S D

B L Y P

B 3L Y P

E x p .a

S c +

C 2v

a 11a 113

O h e g 2

0.36

0.33

0.30

0.44

T i +

O h

e g 2a 1g 14

F O h

t 2g 3

–0.39

–0.26

–0.21

–0.07

V +

O h

t 2g 3a 1g 15D C 2v

a 11

b 11b 21a 21

–0.94

–0.65

–0.55

C r +

C 2v

a 11

b 11b 21a 21a 116S

O h

e g 2 t 2g 3

–1.99

–1.70

–1.67

–1.78

M n

O h

e g 2 t 2g 3a 1g 1

C 2v

a 12a 11

b 11b 21a 21

1.21

0.95

1.22

1.54

F e

C 2v

a 12a 11

b 11b 21a 21a 11

O h

e g 4t 2g 3

–0.28

–0.41

–0.18

–0.07

C o

O h

e g 4t 2g 3a 1g 1

O h

e g 2t 2g 6

–1.15

–0.96

–0.73

–0.80

N i +

O h

e g 2t 2g 6a 1g 1

C 2v

a 12a 11

b 12b 22a 22

–0.80

–1.48

–1.22

–1.48

C u +

C 2v

a 12a 11

b 12b 22a 22a 11

O h e g 4t 2g 6

–3.63–3.17–2.94–3.26

C o r r e c t e d f o r r e l a t i v i s t i c e f f e c t s , t a k e n f r o m R a g h a v a c h a r i a n d T r u c k s , 1989b .

09 Relative Energies and Thermochemistry

相关文档

最新文档