Systematic optimization of long-range corrected hybrid density functionals

Systematic optimization of long-range corrected hybrid density functionals

Jeng-Da Chai^a兲 and Martin Head-Gordon^b兲

Department of Chemistry, University of California and Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

共Received 26 November 2007; accepted 3 January 2008; published online 27 February 2008兲

A general scheme for systematically modeling long-range corrected共LC兲 hybrid density functionals is proposed. Our resulting two LC hybrid functionals are shown to be accurate in thermochemistry, kinetics, and noncovalent interactions, when compared with common hybrid density functionals.

The qualitative failures of the commonly used hybrid density functionals in some “difficult problems,” such as dissociation of symmetric radical cations and long-range charge-transfer excitations, are significantly reduced by the present LC hybrid density functionals. © 2008 American Institute of Physics.关DOI:10.1063/1.2834918兴

I. INTRODUCTION

In the last two decades, density functional theory¹共DFT兲 based on the Kohn–Sham共KS兲 approach^2,3has been attract- ing considerable attention.^4,5 Due to its favorable scaling with system size and reasonable accuracy in many applica- tions, KS-DFT has been regarded as one of the most power- ful theoretical tools for studying both electronic and dynamic properties of medium to large ground-state systems. Re- cently, the development of time-dependent density functional theory 共TDDFT兲 for treating excited-state systems has also been making considerable progress.^6,7

In KS-DFT, the exact exchange-correlation energy func- tional E_xc关␳兴, however, remains unknown, and needs to be approximated. Functionals based on the local spin density approximation共LSDA兲 have been successful for nearly-free- electron systems.^4,5However, for molecular systems, where electron densities are highly nonuniform, the severe overbinding tendency of LSDA means it is not sufficiently accurate for most quantum chemical applications.

Functionals based on the semilocal generalized gradient approximations 共GGAs兲 have considerably reduced the er- rors associated with the LSDA and have shown reasonable accuracy for atomization energies of many strongly bound systems.⁸ For some weakly bound systems, such as hydrogen-bonded systems, GGAs are still reasonable for the energetics and geometries. However, GGAs can completely fail for van der Waals systems. For such systems, GGAs give insufficient binding or even unbound results. Moreover, GGAs tend to give predicted barrier heights of chemical re- actions that are usually seriously underestimated.

Both of the LSDA and GGAs 共commonly denoted as DFAs for density functional approximations兲 are based on the localized model exchange-correlation holes. The exact exchange-correlation hole is, however, fully nonlocal. There- fore, the success of DFAs is commonly believed to be due to a cancelation of errors between the DFAs for exchange and

correlation.^4,5In situations where the cancelation of errors is not complete, the DFAs can produce erroneous results. No- ticeably, some of these situations occur in the asymptotic regions of molecular systems, where the electron densities decay exponentially. In such regions, due to the severe self- interaction error 共SIEs兲 of DFAs, the DFA exchange- correlation potential exhibits an exponential decay, instead of the correct −1/r decay. This leads to many qualitative fail- ures for problems such as dissociation of cations with odd number of electrons or even the alkali halides.^9,10 In time- dependent DFT, SIE causes dramatic failures for long-range charge-transfer excitations of two well-separated molecules.^11–13 The spatially localized nature of DFAs also leads to the absence of London forces, which are a long- range correlation effect. Therefore, to circumvent the above difficulties, it seems necessary to incorporate part or all of the nonlocality of the exchange-correlation hole into the DFAs.

Hybrid DFT methods, which combine KS-DFT with wave function theory 共WFT兲, are promising as a cost- effective way to incorporate nonlocality of the exchange- correlation hole into the DFAs. They can provide reasonable accuracy for treating large-scale systems. In fact, the most widely used density functionals in quantum chemistry are all hybrid functionals! This happy marriage of KS-DFT and WFT was first proposed by Becke,¹⁴who argued that mixing a small fraction of the exact Hartree–Fock 共HF兲 exchange 共associated with the Kohn–Sham reference wave function兲 with DFAs will provide the desired nonlocality and thereby generally improve the DFA results. The general form of a hybrid density functional can be written as

E_xc= cxE_x^HF+ E_xc^DFA, 共1兲 where cxis a small fractional number, typically ranging from 0.2 to 0.25 for thermochemistry,¹⁴ and from 0.4 to 0.6 for kinetics.¹⁵

Indeed, a remarkable accuracy has been achieved by hy- brid density functionals. For example, one of the most widely used hybrid density functionals, B3LYP,^14,16 a com- bination of the B88 exchange functional¹⁷and the LYP cor- relation functional,¹⁸ has achieved a better accuracy for

a兲Electronic mail: [email protected].

b兲Author to whom correspondence should be addressed. Electronic mail:

[email protected].

0021-9606/2008/128共8兲/084106/15/$23.00 128, 084106-1 © 2008 American Institute of Physics

many strongly bound systems than the second-order Moller- Plesset perturbation theory 共MP2兲.¹⁹ Since then, there have been considerable efforts to improve E_xc关␳兴 relative to B3LYP. The development of hybrid density functionals has facilitated the transition of KS-DFT from solid-state physics into the realm of quantum chemistry as well.

In 1997, another significant advance in KS-DFT was also made by Becke, who proposed to model exchange- correlation functionals by a systematic procedure.²⁰ Similar to expanding molecular orbitals by linear combinations of atomic orbitals, he proposed to expand E_xc关␳兴 using power series expansions involving only the local spin density and its first derivative, in addition to a small fraction of the HF exchange. The linear coefficients in the expansions are opti- mized from a systematic fitting procedure to a set of reliable experimental data共the so-called training set兲. Due to the high flexibility of his functional forms, his resulting B97 func- tional has achieved impressive accuracy for thermochemis- try. Since then, additional attempts have been made to devise good basis functionals, which have led to many quite accu- rate E_xc关␳兴, such as VSXC,²¹ B97-1,²² B97-2,²³ B97-3,²⁴ BMK,¹⁵the HCTH family,^22,25,26 and M05-2X.²⁷

However, a few serious problems still remain in these global hybrid density functionals. As can been seen from Eq.

共1兲, the exchange-correlation potential decays as −c_X/r, not the correct −1/r decay. This still leads to qualitatively incor- rect results for charge-transfer 共CT兲 excited states of molecules.^11–13 Although the use of full HF exchange may remedy these difficulties, a DFA for correlation is, however, incompatible with the fully nonlocal HF exchange due to the absence of good cancelation of errors between them, al- though efforts to develop entirely new post-Hartree–Fock correlation functionals show promise.²⁸A similar difficulty applies to optimized effective potential approaches at present.²⁹

To make progress, the long-range corrected共LC兲 hybrid density functionals have been receiving increasing attention.^30–44 LC hybrids retain full HF exchange only for long-range electron-electron interactions共i.e., the asymptotic regime兲, and thereby resolve a significant part of the self- interaction problems associated with global hybrid function- als. However, the currently used LC hybrid functionals are still not as accurate as the best global hybrid functionals, especially for thermochemistry.

Aiming to improve on this situation, here we propose suitable basis functionals for constructing LC hybrid density functionals. Our two resulting LC hybrid density functionals are shown to be accurate in many applications, such as ther- mochemistry, kinetics, and noncovalent interactions, when compared with the widely used global hybrid functionals.

The rest of this paper is organized as follows. In Sec. II, we briefly describe the rationale for the LC hybrid approach. In Sec. III, we propose the suitable basis functionals for system- atically generating accurate LC hybrid functionals. The per- formance of the two new LC hybrid functionals is compared with that of other functionals in Sec. IV共on the training set兲, and in Sec. V共on some test sets兲. In Sec. VI, we give our conclusions.

II. RATIONALE FOR THE LC HYBRID SCHEME

For the LC hybrid scheme, one first defines long-range 共LR兲 and short-range 共SR兲 operators to partition the Cou- lomb operator. In the first LC scheme, proposed by Savin and co-workers, the LR part is treated by WFT 关such as configuration interaction共CI兲兴, and the SR part is treated by DFT.^30–36Its advantage is to reduce the cost of CI calculation in a finite set of one-electron basis functions, as the LR op- erator is chosen to be nonsingular at electron-electron coa- lescence共and hence the basis does not have to represent the cusp兲. Since DFAs perform well for the short-range interac- tion, this type of approach has been gaining some attention.

However, the need for high-level WFT for the LR interaction and the need to develop a generally accurate SR exchange- correlation functionals still hinder its progress.

Using the LSDA expression for the SR exchange from Refs. 31 and 45, a simplified LC hybrid scheme was first proposed by Iikura et al.³⁷In this scheme, the LR exchange is treated exactly by HF theory, while the SR exchange is approximated by DFAs, and the correlation functional re- mains the same as that of the full Coulomb interaction,

E_xc^LC-DFA= E_x^LR-HF+ E_x^SR-DFA+ E_c^DFA. 共2兲 This greatly reduces the computational cost of the LC hybrid scheme, as the cost now is almost the same as the existing global hybrid scheme! In this work, we therefore focus on this type of LC hybrid scheme. The remaining problems are the choice of the SR and LR operators, the development of an accurate SR exchange density functional, and the devel- opment of a correlation functional that is compatible with it.

The most popular type of splitting operator used in the LC hybrid scheme is the standard error function共erf兲,

1

r₁₂=erf共␻^r12兲

r₁₂ +erfc共␻^r12兲

r₁₂ , 共3兲

where r₁₂⬅兩r12兩=兩r1− r₂兩 共atomic units are used throughout this paper兲. On the right hand side of Eq.共3兲, the first term is long ranged, while the second term is short ranged. The pa- rameter␻defines the range of these operators. In principle, different types of operators can also be used in the LC hybrid approaches. In this work, we employ the erf/erfc partition, as it is particularly straightforward to implement efficiently.⁴⁶

For the simplest LC hybrid functional, the local spin density approximation is used for the DFAs. The LR HF exchange E_x^LR−HFis computed by the occupied spin orbitals

␺i␴共r兲 with the LR operator,

E_x^LR-HF= −¹₂

兺

␴

兺

i,j

occ.

冕冕

⫻erf共␻^r12兲

r₁₂ ␺i␴共r2兲␺j␴共r2兲dr1dr₂, 共4兲 while the analytical form of the SR LSDA exchange func- tional E_x^SR-LSDA can be obtained by the integration of the square of the LSDA density matrix with the SR operator,⁴⁵

E_x^SR-LSDA=

兺

␴

冕

Here, e_x^SR-LSDA_␴ 共␳␴兲 is the SR LSDA exchange energy density for ␴^-spin,

e_x^SR-LSDA_␴ 共␳_␴兲 = −3

2

冉

冊

where k_F␴⬅共6␲²␳_␴共r兲兲^1/3 is the local Fermi wave vector, and a_␴⬅␻/共2kF␴兲 is a dimensionless parameter controlling the values of the attenuation function F共a_␴兲,

F共a_␴兲 = 1 −8

3a_␴

冋 ^冑

^冉

冊

+共2a_␴− 4a_␴³兲exp

冉

冊 册

Retaining the LSDA correlation functional E_c^LSDA, one then has the simplest range-separated extended LDA 共RSHX- LDA兲 hybrid functional,³⁹

E_xc^RSHXLDA= E_x^LR-HF+ E_x^SR-LSDA+ E_c^LSDA. 共8兲 The optimal ␻ values for RSHXLDA were found to be 0.5 bohr⁻¹ 共Refs.39兲 for molecular systems, and 0.4 bohr⁻¹ for solid-state systems.⁴⁰However, due to its insufficient ac- curacy for thermochemistry, the development of gradient- corrected LC hybrid functionals continues attracting much attention, and will be our focus.

A relatively narrow range of ␻ ^values 共from 0.2 to 0.5 bohr⁻¹兲. 共Refs. 31,32,34, and 37–44兲 have been found for existing LC hybrid functionals by optimization of properties of interest. As can been seen in Eq.共3兲, the smaller the␻value is, the longer ranged the SR operator will be. As a result, the use of a small␻value in a LC hybrid functional implies that its SR exchange, which is actually not so short ranged, is approximated by spatially localized DFAs. Since the DFA exchange hole is semilocal and it strictly follows its reference electron, for relatively small ␻ values, the nonlo- cality of the exchange hole for this not-very-short-ranged electron-electron exchange interaction should still be impor- tant, and may not be adequately captured by SR DFA ex- change alone.

To remedy this, we argue that mixing with a small amount of the SR HF exchange should be helpful. Similar to Becke’s adiabatic-connection argument¹⁴ for mixing a frac- tion of HF exchange with DFT, mixing a small fraction of the SR HF exchange with the SR DFA exchange should also improve thermochemistry and provide the desired nonlocal correction to the SR exchange. Furthermore, this does not affect the already correct LR behavior of the LC hybrid func- tionals. Similar arguments for the importance of the SR HF exchange were also made in Ref.47.

Hence, we propose the following expression for the LC hybrid functionals:

E_xc^LC-DFA= E_x^LR-HF+ c_xE_x^SR-HF+ E_x^SR-DFA+ E_c^DFA, 共9兲 where E_x^SR-HFis the SR HF exchange,

E_x^SR-HF= −1 2

兺

兺

i,j

occ.

冕冕

⫻erfc共␻^r12兲

r₁₂ ␺i␴共r2兲␺j␴共r2兲dr1dr₂, 共10兲 and cxis a fractional number to be determined.

III. SYSTEMATIC OPTIMIZATION

From the above arguments, the key ingredient for a suc- cessful LC hybrid functional is to construct a generally ac- curate E_x^SR-DFAthat is compatible with the E_c^DFA, the fraction of E_x^SR-HF, and the full E_x^LR-HF. Since the optimal␻ ^{for LC} hybrid scheme is expected to be small, the optimal form of E_x^SR-DFAshould be close to that of E_x^DFA. Therefore, a minor modification to E_x^DFAmay provide a good starting point for developing accurate E_x^SR-DFA.

Since the uniform electron gas 共UEG兲 limit of SR ex- change is believed to be the leading contribution to the SR DFA exchange, and cannot be satisfied by any E_x^DFA共unless

␻^{= 0}兲, we remedy this by replacing the LSDA exchange en- ergy density e_x^LSDA_␴ 共␳_␴兲 with the SR-LSDA exchange energy density e_x^SR-LSDA_␴ 共␳␴兲 关in Eq.共6兲兴, while retaining its enhance- ment factor 共gradient-corrected terms兲. In general, the en- hancement factor of the SR-DFA exchange should be

␻-dependent, as the second-order gradient expansion of SR exchange depends on␻^.33 For a sufficiently small␻ ^value, however, our proposed functional form should be a good approximation.

To achieve a flexible functional form to represent the SR DFA exchange, we modify the B97 exchange functional²⁰by replacing e_x^LSDA_␴ 共␳_␴兲 with exSR-LSDA␴ 共␳_␴兲 关in Eq. 共6兲兴, and de- note this functional as SR-B97 共short-range B97兲 exchange, as it reduces to the B97 exchange functional at␻^{= 0.}

E_x^SR-B97=

兺

冕

g_x␴共s_␴²兲 =

兺

i=0 m

c_x␴,iu_xⁱ_␴, 共12兲

where gx␴共s_␴²兲 is a dimensionless inhomogeneity correction factor depending on the dimensionless reduced spin density gradient s_␴=兩ⵜ␳_␴兩/␳_␴^4/3, and the expansion function ux␴,

ux␴=␥x␴s_␴²/共1 +␥x␴s_␴²兲, 共13兲

␥x␴= 0.004. 共14兲

We use the same form for the correlation functional as the B97 correlation functional, which can be decomposed into same-spin E_c^B97_␴␴ and opposite-spin E_c^B97_␣␤components,

E_c^B97=

兺

␴

E_c^B97_␴␴+ E_c^B97_␣␤. 共15兲

For the same-spin terms,

E_c^B97_␴␴=

冕

gc␴␴共s_␴²兲 =

兺

i=0 m

cc␴␴,iu_cⁱ_␴␴, 共17兲

u_c␴␴=␥c␴␴s_␴²/共1 +␥c␴␴s_␴²兲, 共18兲

␥c␴␴= 0.2, 共19兲

and for the opposite-spin terms,

E_c^B97_␣␤=

冕

gc␣␤共s_av²兲 =

兺

i=0 m

cc␣␤,iu_cⁱ_␣␤, 共21兲

u_c␣␤=␥c␣␤s_av²/共1 +␥c␣␤s_av²兲, 共22兲

␥c␣␤= 0.006, 共23兲

s_av² =¹₂共s_␣²+ s_␤²兲. 共24兲 The correlation energy densities e_c^LSDA_␴␴ and e_c^LSDA_␣␤ are derived from Perdew–Wang⁴⁸ parametrization of the LSDA correla- tion energy, using the approach of Stoll et al.,⁴⁹

e_c^LSDA_␴␴ 共␳_␴兲 = ecLSDA共␳_␴^,0兲, 共25兲 e_c^LSDA_␣␤ 共␳_␣^,␳_␤兲 = ecLSDA共␳_␣^,␳_␤兲 − ecLSDA共␳_␣^,0兲

− e_c^LSDA共0,␳␤兲. 共26兲 Based on the above functional expansions, we propose two new LC hybrid functionals, ␻^{B97 and} ␻^B97X. ␻^B97 has no SR HF exchange共like most of the LC hybrid func- tionals兲,

E_xc^␻B97= E_x^LR-HF+ E_x^SR-B97+ E_c^B97. 共27兲 By contrast,␻B97X contains a small fraction of the SR HF exchange 共the “X” stands for the use of the SR HF ex- change兲,

E_xc^␻B97X= E_x^LR-HF+ cxE_x^SR-HF+ E_x^SR-B97+ E_c^B97. 共28兲 We determined the optimal ␻ values, the linear expan- sion coefficients, the expansion order m of E_xc^␻B97and E_xc^␻B97X by least-squares fittings to 412 accurate experimental and accurate theoretical results共the training set兲, including the 18 atomic energies from the H atom to the Ar atom,⁵⁰ the at- omization energies of the G3/99 set^51–53共223 molecules兲, the ionization potentials 共IPs兲 of the G2-1 set⁵⁴ 关40 molecules, excluding SH₂共²A₁兲 and N2共²⌸兲 cations due to the known convergence problems for pure density functionals⁵²兴, the electron affinities共EAs兲 of the G2-1 set 共25 molecules兲, the proton affinities共PAs兲 of the G2-1 set 共8 molecules兲, the 76 barrier heights of the NHTBH38/04 and HTBH38/04 sets,^55,56and the 22 noncovalent interactions of the S22 set.⁵⁷ All data are equally weighted in the least-squares fitting. By

choosing a diverse range of training data, our goal is to achieve optimized functionals whose performance is well- balanced across typical applications.

We enforce the exact UEG limit for the ␻^{B97 and}

␻B97X functionals by imposing the following constraints:

cc␴␴,0= 1, 共29兲

cc␣␤,0= 1, 共30兲

cx␴,0= 1, for E_xc^␻B97, 共31兲

and

cx␴,0+ cx= 1, for E_xc^␻B97X. 共32兲 Searching for the optimal parameters for ␻^{B97 and}

␻B97X naively seems impractical due to the use of large number of empirical parameters. Following Van Voorhis and Scuseria,²¹ our parameters are obtained by an iterative pro- cedure. First, we focus on a limited range of possible ␻ values between 0.0 and 0.5 bohr⁻¹based on those studied in previous LC hybrid functionals. For each␻ ^value共0.0, 0.1, 0.2, 0.3, 0.4, or 0.5 bohr⁻¹兲, the corresponding RSHXLDA orbitals³⁹ are used as the initial guess orbitals for least- squares fitting. We then obtain, for each␻value, a new set of linear expansion coefficients. With this new set of linear ex- pansion coefficients, the corresponding self-consistent orbit- als can be obtained and then used for another least-squares fitting. This procedure is repeated, for each␻value, until the energies and the linear expansion coefficients are sufficiently close to the previous ones.

Interestingly, we have found that the statistical errors obtained in the first cycle共using the RSHXLDA orbitals兲 are not very different from those obtained self-consistently, even though the linear expansion coefficients in different iterative cycles can be different. This indicates that one could have a good estimate of the performance of proposed functionals, for each␻value, even in its first iterative cycle共this is some-

TABLE I. Optimized parameters for the␻B97关in Eq.共27兲兴, and␻B97X关in Eq.共28兲兴.

␻^B97 ␻^B97X

␻ ^{0.4 bohr−1 0.3 bohr−1}

c_x␴,0 1.00000E + 00 8.42294E − 01

cc␴␴,0 1.00000E + 00 1.00000E + 00

cc␣␤,0 1.00000E + 00 1.00000E + 00

cx␴,1 1.13116E + 00 7.26479E − 01

cc␴␴,1 −2.55352E + 00 −4.33879E + 00

cc␣␤,1 3.99051E + 00 2.37031E + 00

cx␴,2 −2.74915E + 00 1.04760E + 00

cc␴␴,2 1.18926E + 01 1.82308E + 01

cc␣␤,2 −1.70066E + 01 −1.13995E + 01

cx␴,3 1.20900E + 01 −5.70635E + 00

cc␴␴,3 −2.69452E + 01 −3.17430E + 01

c_c␣␤,3 1.07292E + 00 6.58405E + 00

cx␴,4 −5.71642E + 00 1.32794E + 01

c_c␴␴,4 1.70927E + 01 1.72901E + 01

cc␣␤,4 8.88211E + 00 −3.78132E + 00

c_x 1.57706E − 01

thing we shall employ in assessing the usefulness of the new functional forms兲. All of the self-consistent optimizations are well converged within four iterative cycles.

During the optimization procedure, we found that the statistical errors are not significantly improved for m⬎4.

Thus, the functional expansions employed in ␻^{B97 and}

␻B97X are truncated at m = 4. For␻B97, the range separator

␻^{= 0.4 bohr−1}is found to be optimal, which is the same re- sult found by Vydrov et al.^41,42and in agreement with recent arguments made by Fromager et al.⁵⁸ However, a slightly smaller optimal value, ␻^{= 0.3 bohr−1}, is found for ␻^B97X.

As might be anticipated, the presence of a small fraction of the 共nonlocal兲 SR HF exchange allows the SR part to be longer ranged. The optimized parameters of the ␻^{B97 and}

␻B97X functionals are given in TableI.

The limiting cases where␻^{= 0 for}␻^{B97X and}␻^{B97 are} especially interesting, as these reduce to the same functional forms as the existing B97 共Ref. 20兲 and HCTH 共Ref. 22兲

functionals, respectively. Therefore, it is important to know how well B97 and HCTH perform here, when they are both optimized on the same training set. We thus reoptimize B97 and HCTH functionals on the same training set, truncate their functional expansions at the same order m = 4, and im- pose the same UEG limit. Their optimizations are done in a post-LSDA manner 共LSDA orbitals are used for the linear least-squares fittings, instead of using their self-consistent orbitals兲. As mentioned above, these post-LSDA results are believed to be quite close to the fully optimized self- consistent results. For comparisons within the training set, we denote these two reoptimized functionals as B97^* and HCTH^*, respectively.

IV. RESULTS FOR THE TRAINING SET

All calculations are performed with a development ver- sion ofQ-CHEM 3.0.⁵⁹Spin-restricted theory is used for singlet

TABLE II. Statistical errors共in kcal/mol兲 of the training set, including atomization energies 共AEs兲 of the G3/99 set共233 molecules兲 共Refs.51–53and61兲, ionization potentials 共IPs兲 of the G2-1 set 关40 molecules, except for SH₂共²A₁兲 and N2共²⌸兲 cations兴, electron affinities 共EAs兲 of the G2-1 set 共25 molecules兲, proton affinities 共PAs兲 of the G2-1 set 共8 molecules兲 共Ref.54兲, nonhydrogen transfer barrier heights of the NHTBH38/04 set 共38 barrier heights兲, hydrogen transfer barrier heights of the HTBH38/04 set 共38 barrier heights兲 共Refs.55and56兲, and the S22 set共22 molecules兲 for noncovalent interactions 共Ref.57兲. The B97^*and HCTH*functionals are defined in the text. For all cases, single-point calculations are performed using the 6-311+ + G共3df ,3pd兲 basis set. For the AE, IP, EA, and PA, the geometries and zero-point energies were obtained at the B3LYP/6-31G共2df ,p兲 level using a frequency scale factor of 0.9854 for zero-point energies. For the AE, the scaled共0.9854兲 thermal correction at the B3LYP/6-31G共2df ,p兲 level and experimental spin-orbital corrections for the atoms are also used for reversely converting experimental enthalpies of formation to atomization energies. For the S22 set, counterpoise corrections are used to reduce basis set superposition errors, and monomer deformations are not included in the interaction energies.

System Error ␻^B97X ␻^{B97 B97* HCTH* B97-1 B3LYP BLYP}

MSE −0.09 −0.20 0.54 1.74 −1.58 −4.30 −4.60

G3/99共223兲 MAE 2.09 2.56 2.99 4.80 4.85 5.46 9.77

rms 2.86 3.51 4.19 6.27 6.32 7.35 12.97

MSE −0.15 −0.48 2.33 0.37 −0.29 2.16 −1.51

IP共40兲 MAE 2.69 2.65 3.45 3.70 2.60 3.68 4.42

rms 3.59 3.58 4.59 4.49 3.22 4.80 5.27

MSE −0.43 −1.45 1.03 1.71 −0.90 1.73 0.39

EA共25兲 MAE 2.05 2.67 2.46 2.60 1.95 2.39 2.58

rms 2.59 3.10 3.25 3.73 2.42 3.31 3.20

MSE 0.60 0.68 −0.69 1.31 0.62 −0.75 −1.45

PA共8兲 MAE 1.22 1.45 1.23 1.77 0.99 1.14 1.57

rms 1.72 2.17 1.32 2.27 1.52 1.35 2.10

MSE 0.56 1.32 −2.21 −6.35 −3.14 −4.57 −8.68

NHTBH共38兲 MAE 1.75 2.31 2.67 6.70 3.52 4.69 8.72

rms 2.08 2.82 3.41 8.18 4.26 5.71 10.26

MSE −1.51 −0.34 −2.73 −6.25 −4.76 −4.48 −7.84

HTBH共38兲 MAE 2.24 2.24 2.89 6.34 4.76 4.56 7.84

rms 2.58 2.62 3.18 7.16 5.39 5.10 8.66

MSE 0.53 0.16 2.64 4.96 2.55 3.94 5.04

S22共22兲 MAE 0.87 0.60 2.69 4.96 2.55 3.94 5.04

rms 1.30 0.80 4.06 6.03 3.60 5.16 6.29

MSE −0.15 −0.14 0.27 0.22 −1.58 −2.77 −4.07

All共394兲 MAE 2.05 2.39 2.91 4.83 4.10 4.75 8.05

rms 2.75 3.23 3.98 6.22 5.42 6.39 10.88

states and spin-unrestricted theory for others, unless noted otherwise. Results for the training set are computed using the 6-311+ + G共3df ,3pd兲 basis set with the SG-1 grid⁶⁰ for nu- merically integrating the exchange-correlation contributions.

As is usual in hybrid density functional approaches, the elec- tronic energy is minimized with respect to the orbitals. The overall performance of the two new LC hybrid functionals is compared with B97^*, HCTH^*, B97-1,²² B3LYP,^14,16 and BLYP共Refs.17and18兲 in Table II. The first two compari- sons are particularly significant because they indicate how much improvement is possible with the addition of a single extra parameter corresponding to making long-range ex- change exact and thus self-interaction free. The comparisons with other functionals that are not optimized on the same training set are interesting but not as significant as the com- parisons on independent test sets discussed later.

The error for each entry is defined as 共error

= theoretical value− reference value兲. The notation used for characterizing statistical errors is as follows: Mean signed errors 共MSEs兲, mean absolute errors 共MAEs兲, root-mean- square共rms兲 errors, maximum negative errors 共Max共⫺兲兲, and maximum positive errors共Max共⫹兲兲.

A. Thermochemistry

Satisfactory accuracy for thermochemical calculations is one major criterion to judge the performance of density func- tionals. The G3/99 set,^51–53 compiled by Curtiss et al., is probably the leading standard test set for this purpose. Fol- lowing the procedure of G3X theory,⁶¹ the optimized B3LYP/6-31G共2df ,p兲 geometries and zero-point energies are used for all species. A frequency scale factor of 0.9854 is used for zero-point energies and thermal corrections,⁶¹ and atomic spin-orbital effects 共corrected by experimental results⁵¹兲 are included. Enthalpies of formation of free atoms are taken from experiment.⁵¹ The experimental atomization energies are then obtained from experimental standard en- thalpies of formation共at 298 K兲 with a reverse application of G3X theory.⁶¹

The IPs, EAs, and PAs of the G2-1 set⁵⁴are determined at zero temperature and atomic spin-orbit effects are not con- sidered. As can be seen in TableII,␻^B97X,␻B97, and B97^* provide very accurate thermochemical results, especially at- omization energies. Their performances on the IP, EA, and PA training sets are comparable. The important role of long- range exchange in obtaining good results can be clearly seen by comparing results for ␻B97 with HCTH^*: These func- tionals differ by just one parameter共16 versus 15兲 but results are qualitatively improved. The addition of short-range ex- change controlled by one more mixing parameter in␻^B97X leads to further improvements in the quality of results, which shows that this too is a physically important enhancement to the functional form.

In Table III, we compare the performance of various functionals on atomization energies of the G3/99 set. As can been seen,␻B97X performs best, followed by␻^{B97. Since} the G3/99 set is part of the training set for ␻^{B97X and}

␻B97, one cannot attach much significance to this result—it

is a necessary but not sufficient indication of their potential usefulness. MCY1, MCY2,⁶² and BMK 共Ref. 15兲 perform reasonably well on this test set, while M05-2X,²⁷ LC-␻^{PBE,41,42 and B97-1} 共Ref. 22兲 provide only slight im- provement over the most popular hybrid functional, B3LYP.^14,16 The last two columns in Table III tell whether the functionals obey the exact UEG limit, and whether they are LC hybrid functionals.

B. Kinetics

We evaluate the performance of functionals for barrier heights of chemical reactions in the NHTBH38/04 and HTBH38/04 sets.^55,56 The NHTBH38/04 set contains both forward and reverse barrier heights for 19 non-hydrogen- transfer reactions, and the HTBH38/04 set contains both for- ward and reverse barrier heights for 19 hydrogen-transfer reactions. The optimized geometries and the reference ener- gies are taken from Refs.55and56. As can be seen in Table II,␻^B97X,␻B97, and B97^*provide accurate kinetics, com- pared to HCTH^*, B97-1, B3LYP, and BLYP. It is noticeable that the pure density functional, HCTH^*, severely underesti- mates the barrier heights due to self-interaction errors, de- spite these data being part of its training set. Furthermore, the full inclusion of long-range exact exchange in ␻^B97X and␻B97 improves results relative to B97^*, showing that it is important for improving reaction barriers. Detailed infor- mations for the performance of functionals on these two sets are also given in TablesIVandV.

C. Noncovalent interactions

For noncovalent complexes in the S22 set,⁵⁷we perform calculations with the usual counterpoise corrections⁶³for re- ducing the basis set superposition error 共BSSEs兲. Monomer deformation energies are not included. In Table II, we ob- serve that all the functionals predict underbinding results, except for␻^{B97X and}␻B97. Clearly, B97^*and HCTH^*fail for noncovalent interactions, even though these data are in- cluded in their training set. This indicates that there is limited scope to simultaneously improve thermochemistry, kinetics, and noncovalent interactions by reoptimizing the parameters for B97 and HCTH. By contrast, inclusion of exact long-

TABLE III. Summary of performance of various functionals for the atomi- zation energies共in kcal/mol兲 of the G3/99 set. The last two columns indicate whether the functionals are exact in the uniform electron gas共UEG兲 limit and are long-range corrected共LC兲 hybrid functionals. The results for the MCY1 and MCY2 are taken from Ref.62, and the results for the BMK, M05-2X, and LC-␻PBE PBE are taken from Ref.42.

Functional MSE MAE UEG LC

␻^{B97X −0.09 2.09 Yes Yes}

␻^{B97 −0.20 2.56 Yes Yes}

MCY1 3.16 No Yes

MCY2 3.37 No Yes

BMK 2.68 3.69 No No

M05-2X 2.84 4.16 Yes No

LC-␻^{PBE 0.93 4.25 Yes Yes}

B97-1 −1.58 4.85 No No

B3LYP −4.30 5.46 No No

range exchange in ␻^{B97X and} ␻B97 leads to significant improvement. Detailed information on the S22 set can been seen in TableVI. This shows that B97-1, B3LYP, and BLYP are all underbinding, and completely fail for dispersion- bound complexes, while␻^{B97X and}␻B97 are still doing a reasonable共though not highly accurate兲 job.

V. RESULTS FOR THE TEST SETS

To test the transferability of the performance of the two new LC hybrid functionals, we evaluate their performance outside their training sets and compare the results with other three widely used functionals, B97-1,²² B3LYP,^14,16 and BLYP.^17,18 Calculations are performed on various test sets

involving 48 atomization energies, 30 reaction energies, 29 noncovalent interactions, 166 optimized geometry properties, four dissociation curves of symmetric radical cations, and one long-range charge-transfer excitation curve of two well- separated molecules. There are a total of 278 pieces of data in the test sets.

A. Atomization energies

The additional 48 atomization energies in the G3/05 test set⁶⁴ 共other than the 223 atomization energies in the G3/99 test set^51–53兲 are computed by various density functionals.

This test set can be regarded as one of the most stringent test sets, as it contains third-row elements共none is in our training

TABLE IV. Nonhydrogen transfer barrier heights共in kcal/mol兲 of the NHTBH38/04 set 共Ref.56兲.

Reactions ⌬Eref ␻B97X ␻B97 B97-1 B3LYP BLYP

Heavy-atom transfer reactions

H + N₂O→OH+N2 V^f 18.14 19.22 20.67 15.89 11.36 8.53

V^r 83.22 80.57 81.93 72.49 72.81 61.66

H + FH→HF+H V^f 42.18 43.10 44.78 37.93 31.01 26.03

V^r 42.18 43.10 44.78 37.93 31.01 26.03

H + ClH→HCl+H V^f 18.00 20.73 23.17 16.23 12.42 9.81

V^r 18.00 20.73 23.17 16.23 12.42 9.81

H + FCH₃→HF+CH3 V^f 30.38 32.14 33.46 27.55 21.78 16.11

V^r 57.02 55.41 55.83 49.63 48.63 42.27

H + F₂→HF+F V^f 2.27 0.86 1.96 −2.36 −7.54 −11.67

V^r 106.18 104.27 103.66 98.12 96.17 82.16

CH₃+ FCl→CH3F + Cl V^f 7.43 3.93 4.62 −2.15 −1.56 −6.95

V^r 60.17 58.52 59.96 51.22 51.08 41.90

Nucleophilic substitution reactions

F⁻+ CH₃F→FCH3+ F⁻ V^f −0.34 −2.27 −2.60 −3.74 −3.93 −7.90

V^r −0.34 −2.27 −2.60 −3.74 −3.93 −7.90

F⁻¯CH3F→FCH3¯F⁻ V^f 13.38 13.28 13.32 10.82 10.21 6.47

V^r 13.38 13.28 13.32 10.82 10.21 6.47

Cl⁻+ CH₃Cl→ClCH3+ Cl⁻ V^f 3.10 4.71 6.21 −1.03 −0.57 −3.96

V^r 3.10 4.71 6.21 −1.03 −0.57 −3.96

Cl⁻¯CH3Cl→ClCH3¯Cl⁻ V^f 13.61 16.09 17.74 9.90 9.30 5.81

V^r 13.61 16.09 17.74 9.90 9.30 5.81

F⁻+ CH₃Cl→FCH3+ Cl⁻ V^f −12.54 −13.11 −11.72 −17.00 −16.57 −19.37

V^r 20.11 20.83 20.15 17.99 18.25 13.02

F⁻¯CH3Cl→FCH3¯Cl⁻ V^f 2.89 4.23 5.39 0.68 0.29 −1.68

V^r 29.62 31.19 30.95 27.20 26.68 21.09

OH⁻+ CH₃F→HOCH3+ F⁻ V^f −2.78 −3.70 −4.05 −5.94 −5.83 −9.78

V^r 17.33 17.64 17.86 14.36 14.14 9.08

OH⁻¯CH3F→HOCH3¯F⁻ V^f 10.96 11.47 11.52 8.11 7.69 3.73

V^r 47.20 49.33 49.13 46.44 45.46 40.23

Unimolecular and association reactions

H + N₂→HN2 V^f 14.69 13.99 15.47 11.24 7.47 5.25

V^r 10.72 14.32 15.06 12.39 10.87 8.50

H + CO→HCO V^f 3.17 4.55 5.65 3.08 −0.60 −1.96

V^r 22.68 26.72 27.07 25.40 24.60 23.34

H + C₂H₄→CH3CH₂ V^f 1.72 4.07 4.94 3.23 −0.22 −0.74

V^r 41.75 47.07 48.49 43.14 41.71 38.08

CH₃+ C₂H₄→CH3CH₂CH₂ V^f 6.85 5.04 4.81 3.35 5.97 4.70

V^r 32.97 35.21 36.59 31.00 29.41 24.85

HCN→HCN V^f 48.16 46.29 45.89 46.07 47.39 46.77

V^r 33.11 33.12 32.80 32.71 33.33 31.69

MSE 0.56 1.32 −3.14 −4.57 −8.68

MAE 1.75 2.31 3.52 4.69 8.72

rms 2.08 2.82 4.26 5.71 10.26

Max共⫺兲 −3.50 −2.81 −10.73 −11.17 −24.02

Max共⫹兲 5.32 6.74 2.72 1.92 0.66