Short- and long-range corrected hybrid density functionals with the D3 dispersion corrections

(1)

corrections

Chih-Wei Wang, Kerwin Hui, and Jeng-Da Chai

Citation: The Journal of Chemical Physics 145, 204101 (2016); doi: 10.1063/1.4967814 View online: http://dx.doi.org/10.1063/1.4967814

View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/145/20?ver=pdfcov Published by the AIP Publishing

Articles you may be interested in

Implementation and benchmark of a long-range corrected functional in the density functional based tight- binding method

J. Chem. Phys. 143, 184107 (2015); 10.1063/1.4935095 Wave function methods for fractional electrons

J. Chem. Phys. 139, 074107 (2013); 10.1063/1.4817849

Long-range corrected hybrid meta-generalized-gradient approximations with dispersion corrections J. Chem. Phys. 136, 154109 (2012); 10.1063/1.4704370

Simultaneous benchmarking of ground- and excited-state properties with long-range-corrected density functional theory

J. Chem. Phys. 129, 034107 (2008); 10.1063/1.2954017

Asymptotic correction of the exchange–correlation kernel of time-dependent density functional theory for long- range charge-transfer excitations

J. Chem. Phys. 121, 655 (2004); 10.1063/1.1759320

(2)

Short- and long-range corrected hybrid density functionals with the D3 dispersion corrections

Chih-Wei Wang,¹Kerwin Hui,¹and Jeng-Da Chai^1,2,a)

1Department of Physics, National Taiwan University, Taipei 10617, Taiwan

2Center for Theoretical Sciences and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan

(Received 30 August 2016; accepted 2 November 2016; published online 22 November 2016)

We propose a short- and long-range corrected (SLC) hybrid scheme employing 100% Hartree-Fock exchange at both zero and infinite interelectronic distances, wherein three SLC hybrid density functionals with the D3 dispersion corrections (SLC-LDA-D3, SLC-PBE-D3, and SLC-B97-D3) are developed. SLC-PBE-D3 and SLC-B97-D3 are shown to be accurate for a very diverse range of applications, such as core ionization and excitation energies, thermochemistry, kinetics, noncovalent interactions, dissociation of symmetric radical cations, vertical ionization potentials, vertical electron affinities, fundamental gaps, and valence, Rydberg, and long-range charge-transfer excitation energies. Relative to !B97X-D, SLC-B97-D3 provides significant improvement for core ionization and excitation energies and noticeable improvement for the self-interaction, asymptote, energy-gap, and charge-transfer problems, while performing similarly for thermochemistry, kinetics, and noncovalent interactions. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4967814]

I. INTRODUCTION

Due to its decent balance between cost and performance, Kohn-Sham density functional theory (KS-DFT)^1,2 has been a very popular electronic structure method for studying the ground-state properties of large systems.^3–6Recently, one of its most important extensions, time-dependent density functional theory (TDDFT),⁷has also been actively developed for studying the excited-state and time-dependent properties of large systems.^8–14Nonetheless, the exact exchange-correlation (XC) energy functional Exc[⇢], which is the essential ingre- dient of both KS-DFT and adiabatic TDDFT, has not been found, and hence, density functional approximations (DFAs) for Exc[⇢] have been successively developed to improve the accuracy of KS-DFT and TDDFT for general applications.

Functionals based on the conventional DFAs, such as the local density approximation (LDA),^15,16generalized gradient approximations (GGAs),¹⁷ and meta-GGAs (MGGAs),^18,19 are semilocal density functionals.²⁰They are reasonably accurate for the properties governed by short-range XC effects and are computationally favorable for very large systems. Nev- ertheless, owing to the inadequate treatment of nonlocal XC effects,^3–6,21 semilocal density functionals can perform very poorly for the problems related to the self-interaction error (SIE),²²noncovalent interaction error (NCIE),^23–25and static correlation error (SCE).^26–32

In particular, some of these situations happen in the asymptotic regions (r ! 1) of molecules, where the elec- tron densities decay exponentially. In these regions, owing to the pronounced SIEs associated with semilocal density functionals, the functional derivatives of most semilocal

a)Author to whom correspondence should be addressed. Electronic mail:

[email protected]

density functionals (i.e., the semilocal XC potentials) do not exhibit the correct ( 1/r) decay. Consequently, most semilo- cal density functionals can yield erroneous results for the highest occupied molecular orbital (HOMO) energies^33–38 and high-lying Rydberg excitation energies.^11,39–41 Even if the asymptote problems can be properly resolved by the recently developed semilocal density functionals with correct asymptotic behavior^42–46 and asymptotically corrected model XC potentials,^47–51the SIE problems may remain unre- solved.⁵²Besides, semilocal density functionals are inaccurate for charge-transfer (CT) excitation energies,41,42,52–61due to the lack of a space- and frequency-dependent discontinuity in the adiabatic XC kernel adopted in TDDFT.⁶²

In 1993, on the basis of the adiabatic-connection formal- ism, Becke proposed global hybrid density functionals,^63,64 combining semilocal density functionals with a small fraction (typically ranging from 0.2 to 0.25 for thermochemistry and from 0.4 to 0.6 for kinetics) of Hartree-Fock (HF) exchange.^63–72However, in certain situations, especially in the asymptotic regions of molecular systems, a very large fraction (even 100%) of HF exchange is needed. Widely used global hybrid density functionals, such as B3LYP,^64,65PBE0,^68,69and M06-2X,⁷¹do not qualitatively resolve the SIE, asymptote, and CT problems.^52,73

With the aim of resolving these problems, long-range corrected (LC) hybrid density functionals^74–87 have recently received considerable attention. A commonly used LC hybrid density functional (e.g., LC-!PBE⁷⁷ and !B97⁷⁹) employs 100% HF exchange for the long-range (LR) part of the interelectronic repulsion operator erf(!r12)/r12, a semilocal exchange for the complementary short-range (SR) opera- tor erfc(!r12)/r12, and a semilocal correlation for the entire Coulomb operator 1/r12, with the parameter ! (typically ranging from 0.2 to 0.5 bohr ¹) specifying the partitioning

0021-9606/2016/145(20)/204101/15/$30.00 145, 204101-1 Published by AIP Publishing.

(3)

of the interelectronic distance r12 = |r1 r2| (atomic units are used throughout this paper). Here, erf is the standard error function and erfc is the complementary error function. Besides, the inclusion of a small fraction of HF exchange at short range has been shown to improve the overall accuracy of conventional LC hybrid scheme (e.g., !B97X⁷⁹). Over the years, LC hybrid functionals have been shown to qualitatively resolve the SIE, asymptote, and CT problems, offering a cost-effective way to incorporate nonlocal exchange effects.

To properly account for noncovalent interactions, an accurate description of middle- and long-range dynamical correlation effects is essential. Accordingly, LC hybrid functionals can be combined with the DFT-D (KS-DFT with empiri- cal dispersion corrections) schemes^24,88–92(e.g., !B97X-D,⁸⁰

!M05-D,⁸⁴!M06-D3,⁸⁵and !B97X-D3⁸⁵) and the double- hybrid (adding a small fraction of second-order Møller- Plesset correlation) schemes^72,93–105 (e.g., !B97X-2⁸²).

Alternatively, LC hybrid functionals can also be incorporated with a fully nonlocal correlation density functional for van der Waals interactions (vdW-DF)^106–108(e.g., !B97X-V⁸⁶and

!B97M-V⁸⁷). Recently, we have shown that the !B97 series of functionals (!B97, !B97X, !B97X-D, etc.) has yielded impressive accuracy for various applications,^52,73,109 such as thermochemistry, kinetics, noncovalent interactions, dissociation of symmetric radical cations, frontier orbital energies, fundamental gaps, and valence, Rydberg, and long-range CT excitation energies.

In spite of its general applicability, there are some situations, however, where the !B97 series can fail qualitatively.

Very recently, Maier et al.¹¹⁰ have shown that popular LC hybrid functionals, such as LC-!PBE and !B97X-D, perform very poorly for core excitation energies. They have also shown that global hybrid functionals with a large fraction (about 50%) of HF exchange perform reasonably well for core excitation energies, showing consistency with the previous findings of Nakai and co-workers.^111,112However, global hybrid functionals with 50% HF exchange may not consistently perform well for thermochemistry and many other properties that do not require a large fraction of HF exchange. Within the framework of LC hybrid scheme, Hirao and co-workers have shown that the fraction of HF exchange at short range should be responsible for an accurate description of core excitation energies.¹¹³ Similarly, the short-range corrected hybrid density function- als proposed by Besley et al. have been shown to accurately describe core excitation energies.¹¹⁴

On the other hand, Chai and Head-Gordon have shown that the fraction of HF exchange in the middle-range (MR) region (0.5 bohr . r12 . 1.5 bohr) is important for a good balanced description of thermochemistry and kinetics.¹¹⁵ Besides, they have argued that the fraction of HF exchange in the LR region (r12& 1.5 bohr) should be crucial for the properties sensitive to the tail contributions (e.g., the SIE, asymptote, and CT problems), and the fraction of HF exchange in the SR region (r12 . 0.5 bohr) should be responsible for the proper- ties involving changes in the core contributions to Exc[⇢], such as core excitation energies. However, the SR region of the HF exchange operators adopted in the !B97 series has not been fully explored. Note that the fraction of HF exchange at zero interelectronic distance r12= 0 is only 0.00, 0.16, 0.22, and

0.20 for !B97, !B97X, !B97X-D, and !B97X-D3, respectively. Nonetheless, as the electron densities in the core region are rather high (i.e., close to the high-density limit, where HF exchange should dominate correlation), we argue that a very large fraction of HF exchange in the SR region should be adopted for an accurate description of the properties sensitive to the core contributions (e.g., core ionization and excitation energies).

In this work, we intend to improve the performance of the widely used LC hybrid functionals, LC-!PBE and the !B97 series, for core ionization and excitation energies, while retain- ing similar accuracy for many other applications. Specifically, we propose a new LC hybrid scheme employing 100% HF exchange at r12= 0 (i.e., the LC hybrid scheme is also short- range corrected), which is in strong contrast to the popular LC hybrid scheme (i.e., with the erf operator) and other LC hybrid schemes (e.g., with the erfgau113,116–118and terf^119,120opera- tors) employing vanishing HF exchange at r12= 0. The rest of this paper is organized as follows. We describe the short- and long-range corrected (SLC) hybrid scheme in SectionIIand develop three SLC hybrid density functionals with the D3 dispersion corrections in SectionIII. The performance of our new functionals is compared with other functionals in SectionIV (on the training set) and in SectionV(on various test sets).

Our conclusions are given in SectionVI.

II. SHORT- AND LONG-RANGE CORRECTED (SLC) HYBRID SCHEME

In the SLC hybrid scheme, we first define the short- and long-range (SLR) operator fSLR(r12)/r12, which is an operator that approaches 1/r12 at both the SR (r12= 0) and LR (r12! 1) limits, and the complementary MR opera- tor fMR(r12)/r12= (1 f_SLR(r12))/r12to partition the Coulomb operator,

1

r12 =f_SLR(r12)

r12 +f_MR(r12)

r12 . (1)

In this work, we adopt

f_SLR(r12) = erfc(!SRr₁₂) + erf(!LRr₁₂) (2) as a simple sum of the SR function erfc(!SRr₁₂) and LR function erf(!LRr₁₂). Here, !SRand !LRare parameters con- trolling the SR and LR behavior, respectively, of fSLR(r12).

Accordingly, we have

f_MR(r12) = 1 (erfc(!SRr₁₂) + erf(!LRr₁₂))

= erfc(!LRr12) erfc(!SRr12). (3) After the SLR/MR partition, a SLC hybrid density functional is defined as

E_xc^SLC= E_x^SLR-HF+ E_x^MR-DFA+ E_c^DFA. (4) Here, E_c^DFA is the DFA correlation energy of the Coulomb operator 1/r12 and E_x^SLR-HF is the HF exchange energy of the SLR operator fSLR(r12)/r12= erfc(!_SRr₁₂)/r12+ erf (!LRr₁₂)/r12 (computed by the occupied Kohn-Sham (KS) orbitals { i (r)}),

(4)

E_x^SLR-HF= 1 2

X↵, X^occ.

i, j

"

⇤i (r1) _j^⇤ (r2)

⇥f_SLR(r12)

r₁₂ ^j (r1) i (r2)dr1dr₂

= E_x^SR-HF(!SR) + E_x^LR-HF(!LR), (5) where

E_x^SR-HF(!SR) = 1 2

X↵, X^occ.

i, j

"

⇤i (r1) _j^⇤ (r2)

⇥erfc(!SRr12)

r₁₂ ^j (r1) i (r2)dr1dr2 (6) is the HF exchange energy of the SR operator erfc(!SRr12)/r12

and

E_x^LR-HF(!LR) = 1 2

X↵, ^occ.X

i, j

"

i⇤ (r1) _j^⇤ (r2)

⇥erf(!LRr₁₂)

r12 j (r1) i (r2)dr1dr2 (7) is the HF exchange energy of the LR operator erf(!LRr₁₂)/r12. In addition, E_x^MR-DFA, the DFA exchange energy of the MR operator fMR(r12)/r12= erfc(!_LRr₁₂)/r12 erfc(!SRr₁₂)/r12, can be expressed as

E_x^MR-DFA= E_x^SR-DFA(!LR) E_x^SR-DFA(!SR), (8) where E_x^SR-DFA(!) is the DFA exchange energy of the SR operator erfc(!r12)/r12.

On the basis of Eq.(5), fSLR(r12) can be regarded as the fraction of HF exchange at r12for the SLC hybrid density functional. Therefore, we impose the constraint 0  !LR  !^SR

<1 to ensure that 0  f^MR(r12)  1 and hence, 0  fSLR(r12)

 1 can be satisfied at each r12. For !SR= !_LR, we have f_MR(r12) = 0 and fSLR(r12) = 1, employing the full HF exchange and a DFA correlation. Note that fSLR(r12) (given by Eq.(2)) provides a smooth transition between the following two limits:

fSLR(r12= 0) = 1, lim

r12!1fSLR(r12) = 1, (9) employing 100% HF exchange at both the SR (r12= 0) and LR (r12 ! 1) limits. Note also that the SLC hybrid scheme reduces to the popular LC hybrid scheme (i.e., with the erf operator) as !SR ! 1, while it reduces to pure KS-DFT as

!SR! 1 and !LR= 0.

III. SLC HYBRID FUNCTIONALS WITH DISPERSION CORRECTIONS

On the basis of Eq.(4), here we introduce three SLC hybrid density functionals with the D3 dispersion corrections. As the simplest DFA is the LDA, we define the SLC-LDA functional as

E_xc^SLC-LDA= E_x^SLR-HF+ E^MR-LDA_x + E_c^LDA, (10) where E_c^LDAis the LDA correlation functional,¹⁶ E_x^SLR-HFis the SLR-HF exchange energy (given by Eq.(5)), and

E_x^MR-LDA= E_x^SR-LDA(!LR) E_x^SR-LDA(!SR)

= X ⌅↵,

e^MR-LDA_x dr (11)

is the MR-LDA exchange functional, which is known due to the analytical form of E_x^SR-LDA(!), the LDA exchange func- tional of the SR operator erfc(!r12)/r12.^74,121Here, e^MR-LDA_x is the MR-LDA exchange energy density for -spin,

e^MR-LDA_x = 3 2

3 4⇡

!1/3

⇢^4/3(r) ⇥F(aLR, ) F(aSR, )⇤ , (12) where aLR, ⌘ !^LR/(2(6⇡²⇢ (r))^1/3) and aSR, ⌘ !^SR/ (2(6⇡²⇢ (r))^1/3) are dimensionless parameters controlling the values of the attenuation function F(a),

F(a) = 1 8 3a

"p⇡erf 1 2a

!

3a + 4a³

+ (2a 4a³)exp 1 4a²

!#

. (13)

To go beyond the simplest SLC-LDA, we define the SLC- PBE functional as

E_xc^SLC-PBE= E_x^SLR-HF+ E_x^MR-PBE+ E_c^PBE, (14) where E_c^PBEis the PBE correlation functional,¹⁷E_x^SLR-HFis the SLR-HF exchange energy (given by Eq.(5)), and

E_x^MR-PBE= E_x^SR-PBE(!LR) E_x^SR-PBE(!SR) (15) is the MR-PBE exchange functional, with E_x^SR-PBE(!) being the PBE exchange functional of the SR operator erfc(!r12)/r12.¹²²

To further improve upon SLC-PBE, we adopt flexible functional forms in Eq.(4). Similar to the B97 ansatz,⁶⁷ we define the SLC-B97 functional as

E_xc^SLC-B97= E_x^SLR-HF+ E^MR-B97_x + E_c^B97. (16) Here, E_c^B97has the same functional form as the B97 correlation functional,⁶⁷which can be decomposed into same-spin E_c^B97 and opposite-spin E_c↵^B97components,

E_c^B97= X↵,

E_c^B97 + E_c↵^B97. (17) Here,

E_c^B97 =

⌅ e^LDA_c

Xm i=0

c_c ,i c s² 1 + c s²

!ⁱ

dr, (18)

E_c↵^B97=

⌅ e^LDA_c↵

Xm i=0

c_c↵ _,i* ,

c↵ s²_av 1 + c↵ s²_av+

-

i

dr, (19)

where c = 0.2, _c↵ = 0.006, s²_av = ¹₂(s²_↵+ s²), and s = |r⇢ (r)|/⇢^4/3(r). The correlation energy densities e^LDA_c = e^LDA_c (⇢ , 0) and e^LDA_c↵ = e^LDA_c (⇢↵, ⇢ ) e^LDA_c (⇢↵, 0) e^LDA_c (0, ⇢ ) are derived from the PW92 parametrization of the LDA correlation energy density e^LDA_c (⇢↵, ⇢ ),¹⁶using the approach of Stoll et al.¹²³In addition, E_x^SLR-HFis the SLR-HF

(5)

exchange energy (given by Eq.(5)) and E_x^MR-B97=

X ⌅↵,

e^MR-LDA_x Xm

i=0

c_x _,i ^x s² 1 + x s²

!i

dr (20) is the MR-B97 exchange functional, where x = 0.004 and e^MR-LDA_x is given by Eq. (12). Note that E^MR-B97_x has the same functional form as the SR-B97 exchange functional (see Eq. (11) of Ref.79) when !SR! 1 and has the same functional form as the B97 exchange functional⁶⁷when !SR! 1 and !LR= 0.

Following the DFT-D3 scheme,⁹⁰our total energy is given by

E_DFT-D3= E_KS-DFT+ E_disp(D3), (21) where EKS-DFTis the total energy in KS-DFT and

E_disp(D3) = X

n=6,8

X

A>B

C_n^AB Rⁿ_ABf

1 + 6(sr,nR₀ÂB/R_AB)ⁿ⁺⁸g (22) is the D3 dispersion correction (the unscaled version is adopted, and the three-body term is not included). Here, the second sum is over all atom pairs in the system, and RABis the interatomic distance of atom pair AB, while the cutoff radius RÂB₀ and the dispersion coefficients (C₆ÂB and C₈ÂB) for atom pair AB are provided in the DFT-D3 scheme.⁹⁰Therefore, sr,6

and sr,8, which control the strength of dispersion correction, are the parameters to be determined.

In this work, the SLC-LDA (Eq. (10)), SLC-PBE (Eq.(14)), and SLC-B97 (Eq. (16)) functionals with the D3 dispersion corrections (Eq.(22)) are denoted as SLC-LDA- D3, SLC-PBE-D3, and SLC-B97-D3, respectively. Note that SLC-LDA-D3 and SLC-PBE-D3 satisfy the exact uniform electron gas (UEG) limit by construction, while the exact UEG limit for SLC-B97-D3 is enforced by imposing the following constraints: cx ,0 = c_c _,0= c_c↵ _,0= 1.

The four parameters (!SR, !LR, sr,6, and sr,8) of SLC- LDA-D3 and SLC-PBE-D3 are determined by least-squares fittings to the accurate experimental and theoretical data in the training set, involving

• the 223 atomization energies (AEs) of the G3/99 set,¹²⁴

• the 40 ionization potentials (IPs), 25 electron affinities (EAs), and 8 proton affinities (PAs) of the G2-1 set,¹²⁵

• the 76 barrier heights of the NHTBH38/04 and HTBH38/04 sets,¹²⁶

• the 22 noncovalent interactions of the S22 set.^127,128 For the S22 set, an updated version of reference values from S22B¹²⁸is adopted. For the parameter optimization, we focus on a range of possible !SR(0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, and 4.0 bohr ¹) and !LR(0.30, 0.35, 0.40, 0.45, and 0.50 bohr ¹) values, and optimize the corresponding sr,6and sr,8in steps of 0.001, for 0 < sr,6<2 and 0 < sr,8 <2, respectively. The S22 data are weighted 10 times more than the others. As is usual in hybrid density functional approaches, the electronic energy is minimized with respect to the orbitals. Detailed information about the training set can be found in Refs.79,84and85.

The optimized parameters of SLC-LDA-D3 and SLC- PBE-D3 are summarized in Table I, and the HF exchange operators adopted in SLC-LDA-D3, SLC-PBE-D3, and the

!B97 series are plotted in Figure1. Note that the HF exchange

TABLE I. Optimized parameters for SLC-LDA-D3, SLC-PBE-D3, and SLC-B97-D3. Here, !SRand !LRare defined in Equations(5),(11),(15), and(20), sr,6and sr,8are defined in Equation(22), and the others are defined in Equations(18)–(20).

SLC-B97-D3 SLC-PBE-D3 SLC-LDA-D3

!SR(bohr¹) 2.0 2.0 1.5

!LR(bohr ¹) 0.40 0.40 0.45

s_{r ,6} 1.298 1.179 1.129

sr ,8 1.277 1.123 1.131

cx ,0 1.000000

c_x _,1 1.469313

cx ,2 6.185202

c_x _,3 23.053635

cx ,4 16.353923

c_c ,0 1.000000

cc ,1 2.154721

c_c ,2 10.271378

cc ,3 23.966521

cc ,4 15.345722

c_c↵ _,0 1.000000

cc↵ ,1 4.460711

c_c↵ _,2 25.043202

cc↵ ,3 22.506558

c_c↵ ,4 4.114590

operators adopted in LC-!PBE and !B97 are the same.

As can be seen, the fractions of HF exchange adopted in SLC-PBE-D3 and the !B97 series are similar in the MR region, showing consistency with the previous findings of Chai and Head-Gordon¹¹⁵ that the fine details of the MR region of the HF exchange operators adopted are important for good balanced performance in thermochemistry and kinetics.

Besides, as the LR-HF exchange contributions (see Eq.(7)) in SLC-PBE-D3, LC-!PBE, and !B97 are the same (with !LR

= 0.40 bohr ¹), SLC-PBE-D3, LC-!PBE, and !B97 should have similar performance for the properties sensitive to the tail contributions. In addition, the SR-HF exchange contribution (see Eq.(6)) in SLC-PBE-D3 is significant only in the region of r12 . 1/!SR = 0.5 bohr (i.e., the same as the SR region identified by Chai and Head-Gordon¹¹⁵) and hence, should be responsible only for the properties sensitive to the core contributions. By contrast, for SLC-LDA-D3, a larger fraction of HF exchange is needed to reduce the severe error associated with the underlying LDA. Interestingly, the HF exchange operators adopted in SLC-LDA-D3 and SLC-PBE-D3 look upside down, when compared with those adopted in the MR hybrid functionals developed by Henderson et al. for different purposes.¹²⁹

As the HF exchange operator adopted in SLC-PBE-D3 has been optimized, the same HF exchange operator is adopted in SLC-B97-D3 without further optimization. However, the remaining D3 parameters (sr,6and sr,8) and B97 linear expan- sion coefficients (cx ,i, cc ,i, and cc↵ ,i) of SLC-B97-D3 are determined self-consistently by a least-squares fitting proce- dure described in Ref.79(using the same training set), with the SLC-PBE-D3 orbitals being the initial guess orbitals. Dur- ing the parameter optimization, as the statistical errors of the training set for SLC-B97-D3 are not significantly improved for m >4, the functional expansions adopted in SLC-B97-D3 are

(6)

FIG. 1. Fraction of HF exchange f (r12) as a function of the interelectronic distance r12, for SLC-LDA-D3, SLC-PBE-D3, SLC-B97-D3, and the !B97 series.

truncated at m = 4. We summarize the optimized parameters of SLC-B97-D3 in TableI.

In Secs. IV–V, the overall performance of SLC-LDA- D3, SLC-PBE-D3, and SLC-B97-D3 will be compared with a popular semilocal functional:

• PBE,¹⁷

and several widely used LC hybrid functionals:

• LC-!PBE,⁷⁷

• !B97,⁷⁹

• !B97X,⁷⁹

• !B97X-D,⁸⁰

• !B97X-D3⁸⁵

TABLE II. Statistical errors (in kcal/mol) of the training set. PBE, LC-!PBE, and LC-!PBE-D3 (statistical errors are given in parentheses) were not particularly parametrized using this training set.

System Error PBE LC-!PBE(-D3) !B97 !B97X !B97X-D !B97X-D3 SLC-LDA-D3 SLC-PBE-D3 SLC-B97-D3

G3/99 MSE 20.90 3.12 (4.79) 0.29 0.20 0.24 0.14 2.63 0.57 0.32

(223) MAE 21.51 5.86 (7.17) 2.63 2.13 1.93 2.06 8.84 4.49 2.63

rms 26.30 7.43 (9.02) 3.58 2.88 2.77 2.81 11.09 5.91 3.49

IP MSE 0.04 2.86 (2.85) 0.50 0.14 0.20 0.07 11.60 1.85 0.22

(40) MAE 3.44 4.29 (4.29) 2.68 2.69 2.75 2.66 11.60 3.74 2.54

rms 4.35 5.39 (5.39) 3.60 3.59 3.62 3.53 12.40 4.70 3.45

EA MSE 1.72 0.18 (0.18) 1.52 0.47 0.07 0.37 8.96 0.54 1.66

(25) MAE 2.42 3.00 (3.01) 2.72 2.04 1.91 1.93 8.96 3.05 2.66

rms 3.06 3.50 (3.51) 3.11 2.57 2.38 2.41 9.62 3.52 3.06

PA MSE 0.83 0.86 (0.94) 0.67 0.56 1.42 1.10 1.91 0.84 0.80

(8) MAE 1.60 1.41 (1.45) 1.48 1.21 1.50 1.29 2.31 1.36 1.44

rms 1.91 2.04 (2.08) 2.18 1.70 2.05 1.92 2.54 2.01 2.19

NHTBH MSE 8.52 1.39 (1.01) 1.32 0.55 0.45 0.04 1.99 1.29 1.38

(38) MAE 8.62 2.47 (2.28) 2.32 1.75 1.51 1.53 3.32 2.38 2.13

rms 10.61 3.07 (2.83) 2.82 2.08 2.00 1.89 3.77 2.86 2.55

HTBH MSE 9.67 0.77 ( 1.23) 0.66 1.55 2.57 2.08 0.27 1.03 0.96

(38) MAE 9.67 1.39 (1.59) 2.11 2.27 2.70 2.40 1.99 1.41 2.04

rms 10.37 1.90 (2.07) 2.47 2.60 3.10 2.75 2.59 1.77 2.33

S22 MSE 2.71 2.82 ( 0.08) 0.10 0.47 0.14 0.07 0.34 0.11 0.20

(22) MAE 2.71 2.82 (0.26) 0.53 0.79 0.19 0.18 0.45 0.30 0.23

rms 3.73 3.58 (0.35) 0.63 1.11 0.25 0.25 0.61 0.39 0.33

Total MSE 10.32 2.30 (3.01) 0.23 0.21 0.38 0.28 3.38 0.12 0.26

(394) MAE 14.63 4.50 (5.10) 2.42 2.06 1.94 1.97 7.33 3.52 2.36

rms 20.40 6.09 (7.15) 3.27 2.76 2.73 2.69 9.65 4.90 3.15

(7)

TABLE III. Statistical errors (in eV) of the 23 core ionization energies of 14 molecules taken from Ref.134. The relativistic corrections are not considered.

System Error PBE LC-!PBE !B97 !B97X !B97X-D !B97X-D3 SLC-LDA-D3 SLC-PBE-D3 SLC-B97-D3

Core MSE 26.25 20.15 19.39 14.88 13.74 14.10 4.25 2.36 1.53

Ionization MAE 26.25 20.15 19.39 14.88 13.74 14.10 4.27 2.77 2.53

(23) rms 26.48 20.47 19.70 15.10 13.91 14.29 5.27 3.36 2.91

on the training set and various test sets (see thesupplementary material).

IV. RESULTS FOR THE TRAINING SET

All calculations are performed with a development version of Q-Chem 4.3.¹³⁰ Spin-restricted theory is used for singlet states and spin-unrestricted theory for others, unless noted otherwise. For the interaction energies of the weakly bound systems, the counterpoise correction¹³¹ is employed to reduce the basis set superposition error (BSSE).

Results for the training set are computed using the 6-311++G(3df,3pd) basis set with the fine grid EML(75,302), consisting of 75 Euler-Maclaurin radial grid points¹³²and 302 Lebedev angular grid points.¹³³ The error for each entry is defined as error = theoretical value reference value. The nota- tion adopted for characterizing statistical errors is as follows:

mean signed errors (MSEs), mean absolute errors (MAEs), and root-mean-square (rms) errors.

As shown in TableII, SLC-PBE-D3 consistently outperforms PBE and LC-!PBE for the AEs of the G3/99 set and noncovalent interactions of the S22 set, reflecting the effect of the improved HF exchange operator and dispersion correction, respectively. To provide the fairest comparison to SLC-PBE- D3, the performance of LC-!PBE-D3 (i.e., LC-!PBE with the D3 dispersion correction)⁹¹is also examined here. While LC-!PBE-D3 performs similarly to SLC-PBE-D3 for the S22 set due to the inclusion of dispersion correction, LC-

!PBE-D3 performs considerably worse than SLC-PBE-D3 for the G3/99 set.

Owing to its flexible functional forms, SLC-B97-D3 generally outperforms SLC-PBE-D3 and significantly outperforms SLC-LDA-D3 on the training set. Besides, as the fractions of HF exchange adopted in SLC-B97-D3 and !B97 are similar in the MR region, SLC-B97-D3 performs similarly to !B97 for thermochemistry and kinetics, implying that the SR-HF exchange contribution in SLC-B97-D3 does not

degrade its performance for normal chemistry. However, as mentioned previously, the HF exchange operator of SLC-PBE- D3 is adopted in SLC-B97-D3 (i.e., without further optimization), though the D3 parameters and B97 linear expansion coefficients of SLC-B97-D3 are optimized on the training set. Therefore, SLC-B97-D3 performs slightly worse than

!B97X-D3 (where the HF exchange operator, D3 parameters, and B97 linear expansion coefficients were fully optimized on the same training set). All the dispersion-corrected functionals perform reasonably well for the noncovalent interactions of the S22 set.

V. RESULTS FOR THE TEST SETS

To examine how SLC-LDA-D3, SLC-PBE-D3, and SLC- B97-D3 perform outside the training set, we also assess their performance on a wide variety of test sets, including

• the 23 core ionization energies of 14 molecules,¹³⁴

• the 38 core excitation energies of 13 molecules,¹¹⁴

• the 66 noncovalent interactions of the S66 set,¹³⁵

• four dissociation energy curves of symmetric radical cations,²²

• the 113 AEs of the AE113 database,^52,84

• the 131 vertical IPs of the IP131 database,⁸⁴

• the 131 vertical EAs of the EA131 database,^52,84

• the 131 fundamental gaps of the FG131 database,^52,84

• the 19 valence and 23 Rydberg excitation energies of five molecules,¹³⁶

• one long-range CT excitation energy curve of two well- separated molecules.^55,137

As will be discussed later, each vertical IP can be computed in two different ways, each vertical EA can be computed in three different ways, and each fundamental gap can be computed in three different ways. Consequently, there are in total 1335 pieces of data in the test sets, which are larger and more diverse than the training set.

TABLE IV. Statistical errors (in eV) of the 38 core excitation energies of 13 molecules taken from Ref.114. The relativistic corrections are not considered.

State Error PBE LC-!PBE !B97 !B97X !B97X-D !B97X-D3 SLC-LDA-D3 SLC-PBE-D3 SLC-B97-D3

Core ! MSE 42.32 41.30 40.31 31.95 28.74 30.05 4.81 1.46 0.38

Valence MAE 42.32 41.30 40.31 31.95 28.74 30.05 5.12 2.22 2.53

(15) rms 50.91 49.93 48.96 39.29 35.56 37.08 6.25 2.91 2.84

Core ! MSE 32.26 29.35 28.38 22.26 20.35 21.12 3.36 2.93 1.93

Rydberg MAE 32.26 29.35 28.38 22.26 20.35 21.12 3.50 3.22 2.94

(23) rms 39.91 37.64 36.78 29.43 26.90 27.93 4.94 3.81 3.26

(8)

TABLE V. Statistical errors (in kcal/mol) of the S66 set.¹³⁵

S66 MSE 2.22 2.46 0.15 0.16 0.30 0.23 0.04 0.06 0.35

(66) MAE 2.23 2.46 0.37 0.49 0.35 0.26 0.21 0.27 0.37

rms 2.75 2.80 0.47 0.65 0.51 0.35 0.30 0.35 0.46

A. Core ionization energies

To assess the accuracy of the density functionals on core ionization energies, the 23 core ionization energies of 14 molecules are collected from Ref. 134, where the atoms at which the 1s electrons are ionized are all first-row elements.

As discussed by Baerends and co-workers,¹³⁸ the ionization energies for all the occupied orbitals can be well approximated by the minus orbital energies, when the exact (or highly accurate) XC potential is adopted. Therefore, in this work, the core ionization energy of a molecule is calculated as the minus core orbital energy of the molecule, using the 6-311++G(3df,3pd) basis set and EML(75,302) grid.

As shown in TableIII, PBE performs worst for the core ionization energies, while LC-!PBE and !B97 only have minor improvement due to the vanishingly small fraction of HF exchange at small interelectronic distances. Besides,

!B97X, !B97X-D, and !B97X-D3, which include a small fraction of SR-HF exchange, perform slightly better than LC-

!PBE and !B97. Among the functionals examined on the core ionization energies, SLC-B97-D3 ranks first, while SLC-PBE- D3 and SLC-LDA-D3 rank second and third, respectively.

Overall, the SLC hybrid functionals are comparable in performance and are much more accurate than PBE, LC-!PBE, the !B97 series, and possibly, other LC hybrid functionals employing a small fraction of HF exchange in the SR region, reflecting that a very large fraction of HF exchange in the

SR region is indeed essential for an accurate description of core ionization energies. While the relativistic corrections are not considered here, our comments remain the same for the core ionization energies with the relativistic corrections (see thesupplementary material).

B. Core excitation energies

To examine if our SLC hybrid functionals also improve upon the other functionals for core excitation energies, we take the 38 core excitation energies of 13 molecules from Ref.114, containing a total of 15 core!valence and 23 core!Rydberg excitation energies for the first- and second-row nuclei (from the 1s core orbitals). In conventional TDDFT, the calcula- tions of core excited states can be prohibitively expensive, owing to the large number of roots required to obtain the high energy core excited states. Following Besley et al.,¹¹⁴ we perform TDDFT calculations using the Tamm-Dancoff approximation (TDA)¹³⁶ within a reduced single excitation space (which includes only excitations from the core orbitals of interest),¹³⁹ to reduce the computational costs of core excitation energies. The calculations are performed with the 6-311(2+, 2+)G** basis set and EML(100,302) grid.

For the core excitation energies (see Table IV), PBE, LC-!PBE, and !B97 perform very poorly, while !B97X,

!B97X-D, and !B97X-D3 only have minor improvement,

FIG. 2. Dissociation energy curve of H⁺₂. Zero level is set to E(H) + E(H⁺) for each method.

(9)

FIG. 3. Dissociation energy curve of He⁺₂. Zero level is set to E(He) + E(He⁺) for each method.

due to the small fraction of HF exchange in the SR region.

By contrast, SLC-PBE-D3 and SLC-B97-D3 perform compa- rably, slightly improve upon SLC-LDA-D3, and significantly outperform PBE, LC-!PBE, the !B97 series, and perhaps, other LC hybrid functionals adopting a small fraction of HF exchange in the SR region. For the core excitation energies, the statistical errors associated with SLC-PBE-D3 and SLC- B97-D3 are about one order of magnitude smaller than those associated with PBE, LC-!PBE, and the !B97 series! There- fore, the inclusion of a very large fraction of HF exchange at

small interelectronic distances is also important for accurately describing core excitation energies. While we do not consider the relativistic corrections here, our comments remain similar for the core excitation energies with the relativistic corrections (see thesupplementary material).

C. Noncovalent interactions

For the noncovalent interactions of the S66 set,¹³⁵ the performance of the functionals is evaluated using the

FIG. 4. Dissociation energy curve of Ne⁺₂. Zero level is set to E(Ne) + E(Ne⁺) for each method.

(10)

FIG. 5. Dissociation energy curve of Ar⁺₂. Zero level is set to E(Ar) + E(Ar⁺) for each method.

6-311++G(3df,3pd) basis set and EML(99,590) grid, and the counterpoise correction¹³¹is adopted to reduce the BSSE. As shown in TableV, PBE and LC-!PBE perform very poorly for the noncovalent interactions of the S66 set, due to the lack of a proper description of middle- and long-range dynamical correlation effects, while all the dispersion-corrected functionals perform reasonably well.

D. Dissociation of symmetric radical cations

Due to the pronounced SIEs associated with semilocal density functionals, unphysical fractional charge dissociation can happen, especially for symmetric charged radicals.²² Here, the dissociation energy curves of H⁺₂, He⁺₂, Ne⁺₂, and Ar⁺₂ are calculated using the 6-311++G(3df,3pd) basis set and EML(75,302) grid to examine the performance of the functionals upon the SIE problems. The results are compared with the H⁺₂ curve calculated using the HF theory (exact for any one-electron system) and the He⁺₂, Ne⁺₂, and Ar⁺₂ curves calculated using the highly accurate CCSD(T) theory (coupled-cluster theory with iterative sin- gles and doubles and perturbative treatment of triple substi- tutions).¹⁴⁰

As shown in Figures 2–5, unphysical barriers indeed appear in the PBE dissociation curves, owing to the significant SIEs of PBE. By contrast, the LC and SLC hybrid

functionals greatly reduce (or even remove) the unphysical barriers of the dissociation curves, due to the inclusion of 100% LR-HF exchange. SLC-LDA-D3, adopting the largest

!LR(0.45 bohr ¹), performs best, followed by SLC-PBE-D3, SLC-B97-D3, LC-!PBE, and !B97, adopting the second largest !LR(0.40 bohr ¹).

E. Atomization energies

Recently, we have developed the IP131, EA131, and FG131 databases,^52,84consisting of accurate reference values for the 131 vertical IPs, 131 vertical EAs, and 131 fundamental gaps, respectively, of 18 atoms and 113 molecules at their experimental geometries. In addition, we have developed the AE113 database,⁵² which contains accurate reference values for the atomization energies of 113 molecules in the IP131 database. Here, we examine the performance of the functionals on the AE113, IP131, EA131, and FG131 databases, using the 6-311++G(3df,3pd) basis set and EML(75,302) grid.

As shown in TableVI, owing to their flexible functional forms, SLC-B97-D3 and the !B97 series are comparable in performance, more accurate than SLC-PBE-D3 and much more accurate than PBE, LC-!PBE, and SLC-LDA-D3. Inter- estingly, SLC-PBE-D3 performs better than LC-!PBE, possibly due to the noticeable deviation of their HF exchange operators in the region of 0.5 bohr . r12. 0.8 bohr (where the

TABLE VI. Statistical errors (in eV) of the AE113 database.^52,84

AE113 MSE 0.83 0.10 0.05 0.05 0.04 0.05 0.04 0.03 0.04

(113) MAE 0.88 0.27 0.11 0.10 0.10 0.10 0.28 0.17 0.11

rms 1.06 0.41 0.15 0.13 0.14 0.13 0.34 0.23 0.14

(11)

TABLE VII. Statistical errors (in eV) of the IP131 database.⁸⁴

IP(1) = E_total(N 1) E_total(N)

IP131 MSE 0.26 0.10 0.00 0.00 0.03 0.02 0.57 0.09 0.02

(131) MAE 0.36 0.28 0.19 0.18 0.19 0.18 0.58 0.20 0.18

rms 0.52 0.46 0.26 0.26 0.27 0.26 0.64 0.28 0.26

IP(2) = ✏HOMO(N)

IP131 MSE 4.40 0.15 0.24 0.48 1.01 0.71 0.61 0.09 0.18

(131) MAE 4.40 0.42 0.40 0.51 1.01 0.72 0.70 0.36 0.37

rms 4.50 0.68 0.63 0.75 1.18 0.93 0.77 0.56 0.59

fractions of HF exchange adopted in SLC-PBE-D3,

!B97X-D, and !B97X-D3 are very similar!).

F. Vertical ionization potentials

The vertical IP of a molecule (containing N electrons) is defined as

IP(1) = Etotal(N 1) Etotal(N), (23) where Etotal(N) is the total energy of the N-electron system.

For the exact KS-DFT, the vertical IP of a molecule is the same as the minus HOMO energy of the molecule,^33–38

IP(2) = ✏HOMO(N). (24)

However, for an approximate XC density functional in KS- DFT, the computed IP(1) and IP(2) values may be different, showing the accuracy of the predicted total energies and HOMO energies, respectively.

Here, we examine the accuracy of the functionals on the IP131 database⁸⁴ and summarize our results in Table VII.

For IP(1), the !B97 series, SLC-PBE-D3, and SLC-B97- D3 are comparable in performance, outperforming the other functionals. For IP(2), LC-!PBE, !B97, SLC-PBE-D3, and SLC-B97-D3, which adopt !LR= 0.40 bohr ¹, perform com- parably and outperform the other functionals. By contrast, PBE severely underestimates IP(2), due to the incorrect XC potential asymptote. For the IP(1) and IP(2) values, SLC-PBE- D3 and SLC-B97-D3 achieve the best performance, followed closely by !B97.

G. Vertical electron affinities

The vertical EA of a molecule is defined as

EA(1) = Etotal(N) Etotal(N + 1). (25) By comparing Eq. (23)with Eq. (25), the vertical EA of a molecule is identical to the vertical IP of the corresponding anion, which is, for the exact KS-DFT, the minus HOMO energy of the anion,

EA(2) = ✏HOMO(N + 1). (26) In addition, the vertical EA of a molecule is traditionally approximated by the minus lowest unoccupied molecular orbital (LUMO) energy of the molecule,

EA(3) = ✏LUMO(N). (27)

Nonetheless, even for the exact KS-DFT, there is a fundamental difference between EA(3) and EA(2), owing to the deriva- tive discontinuity xc34,38,50,141–144of Exc[⇢]: EA(3) EA(2)

= ✏_HOMO(N + 1) ✏LUMO(N) = xc. Hybrid density functionals, which belong to the generalized Kohn-Sham (GKS) method¹⁴⁵ (not pure KS-DFT), effectively capture a fraction of xc of Exc[⇢] in KS-DFT. A recent study has found that the difference between ✏HOMO(N + 1) and ✏LUMO(N) is small for LC hybrid functionals.¹⁴⁶Therefore, EA(3) is expected to be close to EA(2) (i.e., the true vertical EA) for LC hybrid functionals.

TABLE VIII. Statistical errors (in eV) of the EA131 database.^52,84

EA(1) = E_total(N) E_total(N + 1)

EA131 MSE 0.10 0.11 0.24 0.19 0.14 0.16 0.23 0.08 0.22

(131) MAE 0.21 0.36 0.34 0.31 0.26 0.29 0.33 0.27 0.32

rms 0.34 0.54 0.44 0.41 0.36 0.39 0.45 0.35 0.42

EA(2) = ✏HOMO(N + 1)

EA131 MSE 2.03 0.00 0.13 0.17 0.32 0.23 0.40 0.08 0.10

(131) MAE 2.03 0.34 0.34 0.32 0.39 0.33 0.47 0.30 0.33

rms 2.30 0.43 0.43 0.41 0.51 0.44 0.62 0.38 0.41

EA(3) = ✏LUMO(N)

EA131 MSE 2.43 0.19 0.38 0.31 0.01 0.15 0.10 0.25 0.34

(131) MAE 2.45 0.42 0.49 0.52 0.54 0.49 0.30 0.39 0.44

rms 2.72 0.51 0.59 0.62 0.63 0.57 0.40 0.47 0.53