† Assessmentofdensityfunctionalmethodswithcorrectasymptoticbehavior PAPER PCCP

(1)

Cite this: Phys. Chem. Chem. Phys.,2013, 15, 8352

Assessment of density functional methods with correct asymptotic behavior†

Chen-Wei Tsai,âYu-Chuan Su,âGuan-De Liâand Jeng-Da Chai*âb

Long-range corrected (LC) hybrid functionals and asymptotically corrected (AC) model potentials are two distinct density functional methods with correct asymptotic behavior. They are known to be accurate for properties that are sensitive to the asymptote of the exchange-correlation potential, such as the highest occupied molecular orbital energies and Rydberg excitation energies of molecules. To provide a comprehensive comparison, we investigate the performance of the two schemes and others on a very wide range of applications, including asymptote problems, self-interaction-error problems, energy-gap problems, charge-transfer problems and many others. The LC hybrid scheme is shown to consistently outperform the AC model potential scheme. In addition, to be consistent with the molecules collected in the IP131 database [Y.-S. Lin, C.-W. Tsai, G.-D. Li and J.-D. Chai, J. Chem. Phys., 2012, 136, 154109], we expand the EA115 and FG115 databases to include, respectively, the vertical electron affinities and fundamental gaps of the additional 16 molecules and develop a new database, AE113 (113 atomization energies), consisting of accurate reference values for the atomization energies of the 113 molecules in IP131. These databases will be useful for assessing the accuracy of density functional methods.

I. Introduction

Over the past two decades, Kohn–Sham density functional theory (KS-DFT)^1,2 has been one of the most powerful approaches for the study of large ground-state systems due to its low computational cost and reasonable accuracy^3–5. Its extension to the real-time domain, time-dependent density functional theory (TDDFT),⁶ has also been actively developed for studying the time-dependent and excited-state properties of large systems.^7–10

In KS-DFT, the exact exchange-correlation (XC) energy functional, Exc[r], remains unknown and needs to be approximated.

Accurate approximations of Exc[r] have been successively developed to extend the applicability of KS-DFT to a wide range of systems. Among them, semilocal density functionals¹¹perform satisfactorily for some applications but can produce erroneous results in situations where the accurate treatment of the non- locality of the XC hole is important.^12–14

One of the most important and challenging topics in KS-DFT is the asymptotic behavior of the XC potential, nxc(r) = dExc[r]/dr(r). For finite systems, the exact nxc(r) should exhibit Coulombic (!1/r) decay as r - N.^15–18However, owing to the severe self-interaction error (SIE),¹⁹ the XC potential of semilocal functionals decays much faster than (!1/r) in the asymptotic region, yielding qualitatively incorrect predictions for properties sensitive to the asymptote of the XC potential, such as the highest occupied molecular orbital (HOMO) energies and high-lying (Rydberg) excitation energies of atoms and molecules.^20,21

To date, perhaps the most successful density functional methods in practice which include a correct XC potential in the asymptotic region of molecules are provided by the long- range corrected (LC) hybrid scheme^22–32 and asymptotically corrected (AC) model potential scheme.^33–42For the LC hybrid scheme, the nonlocal Hartree–Fock (HF) exchange for long- range electron–electron interactions is added to a semilocal functional, which thereby resolves a significant part of the SIE problems and provides an AC XC potential. In contrast, for the AC model potential scheme, an AC XC potential is directly modeled and added to a semilocal functional. In principle, a model XC potential should be a functional derivative of some Exc[r]. However, several existing AC model potentials are not functional derivatives. Consequently, when these AC model

aDepartment of Physics, National Taiwan University, Taipei 10617, Taiwan.

E-mail: [email protected]

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

† Electronic supplementary information (ESI) available. See DOI: 10.1039/

c3cp50441g

Received 30th January 2013, Accepted 27th March 2013 DOI: 10.1039/c3cp50441g www.rsc.org/pccp

PAPER

Downloaded by National Taiwan University on 08/05/2013 19:09:56. Published on 28 March 2013 on http://pubs.rsc.org | doi:10.1039/C3CP50441G

View Article Online

View Journal | View Issue

(2)

potentials are employed, several necessary conditions for a functional derivative, such as (i) minimization of the ground- state energy by the self-consistent KS orbitals, (ii) path inde- pendence of the van Leeuwen–Baerends line integral⁴³and (iii) net zero force and zero torque on the ground-state density, can be violated.⁴⁴ Besides, as these model potentials are not variationally stable, the associated XC energies are not well- defined and the predicted properties need to be carefully inter- preted. The proposals made by Staroverov and co-workers^45,46 are useful for transforming a model GGA exchange potential (which is not a functional derivative) into an exchange potential that has a parent GGA functional. However, for the LB94 potential (an AC model potential with no parent GGA functional),³³ the transformed potential no longer has the correct (!1/r) asymptote.⁴⁶ Further improvements on their scheme are needed to resolve this. Alternatively, one may adopt a density functional whose functional derivative has the correct (!1/r) asymptote,⁴⁷ which should be more desirable than AC model potentials. However, in this work, we only focus on the general performance of the AC model potential scheme.

The rest of this paper is organized as follows. In Section II, we evaluate the performance of the LC hybrid scheme, AC model potential scheme and other density functional methods on a very wide range of applications, including the asymptote problems, SIE problems, energy-gap problems, charge-transfer problems and many others. In Section III, we give our conclusions on the usefulness of these methods.

II. Results and discussion

For a comprehensive comparison of various density functional methods, we examine the performance of the local spin density approximation (LSDA),⁴⁸generalized gradient approximations

(GGAs) (PBE⁴⁹and BLYP⁵⁰), meta-generalized gradient approximations (MGGAs) (M06L⁵¹ and VS98⁵²), global hybrid functionals (B3LYP⁵³and M06-2X⁵⁴), LC hybrid functionals (oB97,²⁷ oB97X²⁷and oB97X-D²⁸) and AC model potentials (LB94³³and LBa³⁵) on various test sets 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 (BHs) of the NHTBH38/04 and HTBH38/04 sets,⁵⁷the 22 noncovalent interactions of the S22 set,⁵⁸the additional 48 AEs of the G3/05 set,⁵⁹the 113 AEs of the AE113 database,⁶⁰ the 131 vertical IPs of the IP131 database,³¹ the 131 vertical EAs of the EA131 database,^31,60 the 131 fundamental gaps (FGs) of the FG131 database,^31,60two dissociation curves of symmetric radical cations, 19 valence excitation energies, 23 Rydberg excitation energies and one long-range charge-transfer excitation curve of two well-separated molecules. As discussed in ref. 31, each vertical EA can be evaluated in two different ways and each FG can be evaluated in three different ways, so there are in total 1386 pieces of data in our test sets, which are very large and diverse. Unspecified detailed information of the test sets, basis sets and numerical grids used is given in ref. 31. The LB94 (or LBa) potential can be expressed as a linear combination of the LSDA exchange potential, the LSDA correlation potential and a gradient- dependent exchange potential (e.g., see eqn (55) of ref. 33).

Due to the inclusion of the gradient-dependent exchange potential, the LB94 (or LBa) potential is not a functional derivative.^44,46 In this work, the exchange energy from the LB94 (or LBa) exchange potential is calculated by the popular Levy–Perdew virial relation:⁶¹

E_x½r# ¼ Z

½3rðrÞ þ r ( rrðrÞ#nxðrÞdr; (1)

Table 1 Statistical errors (in kcal mol^!1) of the oB97 training set²⁷

System Error LSDA PBE BLYP M06L VS98 B3LYP M06-2X oB97 oB97X oB97X-D LB94 LBa

G3/99 MSE 120.59 20.90 !4.59 2.33 !1.11 !4.32 0.36 !0.29 !0.20 !0.24 !482.58 !56.44

(223) MAE 120.60 21.51 9.76 5.37 3.46 5.49 2.30 2.63 2.13 1.93 484.91 103.84

rms 142.50 26.30 12.97 7.07 4.50 7.36 3.39 3.58 2.88 2.77 609.20 171.61

IP MSE 3.42 0.03 !1.51 !1.88 0.52 2.17 0.16 !0.51 !0.15 0.19 97.67 104.03

(40) MAE 5.54 3.46 4.42 3.97 2.94 3.68 2.45 2.67 2.69 2.74 97.92 104.03

rms 6.66 4.35 5.27 4.95 3.66 4.79 3.58 3.60 3.59 3.61 116.03 115.94

EA MSE 6.45 1.72 0.36 !2.89 0.13 1.71 !1.37 !1.52 !0.47 0.07 127.66 116.53

(25) MAE 6.45 2.42 2.57 3.83 2.61 2.38 2.56 2.72 2.04 1.91 127.66 116.53

rms 7.29 3.06 3.17 4.47 3.60 3.27 3.07 3.11 2.57 2.38 136.16 120.87

PA MSE !5.91 !0.83 !1.47 0.97 1.00 !0.77 !1.21 0.67 0.56 1.42 !122.75 !95.66

(8) MAE 5.91 1.60 1.58 1.94 1.60 1.16 2.02 1.48 1.21 1.50 122.75 95.66

rms 6.40 1.91 2.10 2.42 2.18 1.36 2.21 2.18 1.70 2.05 126.54 98.45

NHTBH MSE !12.41 !8.52 !8.69 !3.22 !4.80 !4.57 !0.00 1.32 0.55 !0.45 !58.93 !40.86

(38) MAE 12.62 8.62 8.72 3.73 5.16 4.69 1.41 2.32 1.75 1.51 93.94 68.95

rms 16.13 10.61 10.27 5.26 6.49 5.71 1.91 2.82 2.08 2.00 127.21 90.80

HTBH MSE !17.90 !9.67 !7.84 !4.54 !5.25 !4.48 !0.80 !0.66 !1.55 !2.57 38.65 12.68

(38) MAE 17.90 9.67 7.84 4.54 5.26 4.56 1.32 2.11 2.27 2.70 44.31 22.80

rms 18.92 10.37 8.66 5.22 5.92 5.10 1.56 2.47 2.60 3.10 57.37 29.90

S22 MSE !1.97 2.77 5.05 0.91 !7.18 3.95 0.27 0.15 0.52 !0.08 51.70 36.16

(22) MAE 2.08 2.77 5.05 0.92 7.63 3.95 0.39 0.59 0.86 0.21 51.70 36.16

rms 3.18 3.89 6.31 1.06 11.80 5.17 0.55 0.77 1.28 0.26 61.66 42.55

Total MSE 65.86 10.32 !4.07 0.27 !1.91 !2.79 0.05 !0.23 !0.21 !0.38 !256.11 !16.49

(394) MAE 72.41 14.63 8.05 4.57 3.88 4.77 2.04 2.42 2.07 1.94 310.76 89.67

rms 107.53 20.40 10.88 6.12 5.38 6.39 3.02 3.27 2.77 2.73 463.27 141.88

(3)

while the correlation energy from the LB94 (or LBa) correlation potential is directly calculated by the LSDA correlation energy functional (i.e., its parent functional).

All the calculations are performed with a development version of Q-Chem 3.2.⁶² Spin-restricted theory is used for singlet states and spin-unrestricted theory is used for others unless noted otherwise. For the binding energies of the weakly bound systems, the counterpoise correction⁶³ is employed to reduce basis set superposition error (BSSE). 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) and root-mean-square (rms) errors.

A. xB97 training set

The oB97 training set²⁷consists of several popular databases, such as the 223 AEs of the G3/99 set,⁵⁵the 40 IPs, 25 EAs and 8 PAs of the G2-1 set,⁵⁶ the 76 BHs of the NHTBH38/04 and HTBH38/04 sets⁵⁷ and the 22 noncovalent interactions of the S22 set.⁵⁸ Table 1 summarizes the statistical errors of several methods for the oB97 training set. The results 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.⁶⁵

As shown, oB97, oB97X and oB97X-D perform the best as they all have very flexible functional forms that were parametrized using this training set. Although M06-2X was not particularly parametrized using this training set, it also performs very well.

B3LYP performs similarly to M06L and VS98 and significantly outperforms LSDA, PBE, BLYP, LB94 and LBa. Due to the lack of Exc[r], LB94 and LBa perform the worst (even worse than LSDA).

Therefore, one should avoid using the AC model potential scheme for the calculation of total energies and related properties.

B. Atomization energies

1. G3/05 set. Here, we examine the performance of these functionals on the additional 48 AEs of the G3/05 set⁵⁹using the G3LargeXP basis set⁶⁶and the extra fine grid EML(99 590).

As can be seen in Table 2, oB97, oB97X and oB97X-D remain the best in performance, even if they are not particularly parametrized using this set. On the other hand, due to the lack of Exc[r], LB94 and LBa remain the worst in performance.

2. AE113 database. Recently, we developed the IP131, EA115 and FG115 databases, consisting of accurate reference values for the 131 vertical IPs, 115 vertical EAs and 115 FGs, respectively, of 18 atoms and various molecules in their experimental geometries.³¹ To be consistent with the molecules collected in IP131, EA115 is expanded to EA131 (131 vertical electron affinities)^31,60to include the vertical EAs of the additional

16 molecules and FG115 is expanded to FG131 (131 fundamental gaps)^31,60to include the FGs of the additional 16 molecules. In addition, we develop a new database, AE113 (113 atomization energies),⁶⁰for the AEs of the 113 molecules in IP131.

As described in ref. 31, the reference values in EA131, FG131 and AE113 are obtained via the very accurate CCSD(T) theory (coupled-cluster theory with iterative singles and doubles and perturbative treatment of triple substitutions),⁶⁷ where the CCSD(T) correlation energies in the basis-set limit are extra- polated from calculations using the aug-cc-pVTZ and aug-cc- pVQZ basis sets:⁶⁸

E¹_XY¼E^corr_X X³! EY^corrY³

X³! Y³ ; (2)

with X = 3 and Y = 4 for the aug-cc-pVTZ and aug-cc-pVQZ basis sets, respectively. Note that the reference atomization energies for the AE113 database are calculated without the zero-point energy correction. For the completeness of this work, the reference values of IP131, EA131, FG131 and AE113 are all given in the supplementary material.⁶⁰

The performance of the functionals is examined on the AE113 database using the 6-311++G(3df,3pd) basis set and the fine grid EML(75 302). As shown in Table 3, oB97, oB97X and oB97X-D have the best performance, while LB94 and LBa have the worst performance, which is consistent with their performance on the atomization energies of the G3/99 and G3/05 sets.

C. Frontier orbital energies

For a system of N electrons (N is an integer), the vertical IP is defined as IP(N) ) E(N ! 1) ! E(N) and the vertical EA is defined as EA(N) ) E(N) ! E(N + 1), with E(N) being the ground-state energy of the N-electron system. For the exact KS-DFT, IP(N) is identical to the minus HOMO energy of the N-electron system,^69–74

IP(N) = !eN(N), (3)

and EA(N), which is IP(N + 1) by its definition, is identical to the minus HOMO energy of the (N + 1)-electron system,

EA(N) = !eN+1(N + 1), (4) where ei(N) is the i-th KS orbital energy of the N-electron system.

Conventionally, EA(N) is approximated by the minus LUMO (lowest unoccupied molecular orbital) energy of the N-electron system,

EA(N) E !eN+1(N). (5)

However, it has been found that a difference exists between eN+1(N + 1) and eN+1(N),

Dxc= eN+1(N + 1) ! eN+1(N), (6)

Table 2 Statistical errors (in kcal mol^!1) of the additional 48 atomization energies of the G3/05 set⁵⁹

G3/05 MSE 78.69 16.80 1.01 5.23 1.48 !3.57 !3.54 1.28 0.76 0.24 !427.33 !129.60

(48) MAE 78.83 18.07 9.89 7.50 4.62 6.29 5.14 4.25 3.60 3.01 439.45 155.28

rms 115.99 30.22 15.47 10.05 6.62 8.90 6.66 5.41 4.52 3.95 863.35 395.46

(4)

arising from the derivative discontinuity (DD) of Exc[r],41,70,74–83

Dxc¼ lim

Z!0^þ

dExc½r#

drðrÞ

!!

!!_N

þZ

!dExc½r#

drðrÞ

!!

!!_N

!Z

( )

: (7)

In contrast, as a hybrid functional, containing a fraction of the nonlocal HF exchange, belongs to the generalized KS (GKS) method, the corresponding GKS orbital energies incorporate part of the DD.^84–88A recent study shows that Dxc(see eqn (6)) is close to zero for LC hybrid functionals, so the !eN+1(N) calculated by LC hybrid functionals should be close to EA(N).⁸⁶

Here, we evaluate the performance of the functionals on frontier orbital energies for the IP131 and EA131 databases using the 6-311++G(3df,3pd) basis set with the fine grid EML(75 302) and summarize the results in Table 4.

For the IP131 database, IP(N) is calculated by !eN(N). As shown, the LSDA, GGAs and MGGAs severely underestimate the vertical IPs due to the incorrect asymptotic behavior of the XC potential. In contrast, the LC hybrid functionals and AC model potentials yield very accurate vertical IPs due to the correct asymptote of the XC potential.

For the EA131 database, EA(N) is calculated by both

!e^N+1(N + 1) and !eN+1(N). For the calculations using !eN+1(N + 1) (the minus HOMO energy of the (N + 1)-electron system), the LSDA, GGAs and MGGAs seriously underestimate the vertical EAs due to the incorrect asymptote of the XC potential, while the LC hybrid functionals and AC model potentials perform very well due to the correct asymptote of the XC potential.

On the other hand, for the calculations using !eN+1(N) (the minus LUMO energy of the N-electron system), the situation becomes very different. Although the AC model potentials yield accurate eN+1(N) due to the correct asymptote, the vertical EA remains severely overestimated due to the lack of the DD.

In spite of the lack of the DD, there is a fortuitous cancellation of errors in the vertical EA predicted by the LSDA, GGAs and MGGAs, as the eN+1(N) is (incorrectly) upshifted due to the incorrect asymptote of the XC potential, effectively incorporating part of the DD. As the DD is close to zero for LC hybrid functionals,⁸⁶ i.e., eN+1(N + 1) E eN+1(N), the LC hybrid functionals significantly outperform the others using this estimate.

Table 3 Statistical errors (in eV) of the atomization energies for the AE113 database

AE113 MSE 3.75 0.84 0.13 0.17 0.07 –0.04 0.07 0.05 0.05 0.04 !15.24 !2.36

(113) MAE 3.75 0.88 0.36 0.27 0.14 0.16 0.12 0.12 0.10 0.10 15.53 4.81

rms 4.22 1.06 0.48 0.36 0.20 0.23 0.16 0.15 0.13 0.14 20.02 7.44

Table 4 Statistical errors (in eV) of the frontier orbital energies for the IP131³¹and EA131 databases

!eN(N) ! IPvertical

Atoms MSE !4.88 !4.87 !5.00 !4.66 !4.80 !3.71 !2.22 !0.99 !1.17 !1.63 !0.21 !0.72

(18) MAE 4.88 4.87 5.00 4.66 4.80 3.71 2.22 0.99 1.17 1.63 0.48 0.73

rms 5.23 5.22 5.31 4.94 5.12 3.94 2.32 1.37 1.52 1.97 0.63 0.95

Molecules MSE !4.18 !4.32 !4.42 !4.09 !4.21 !3.06 !1.36 !0.12 !0.37 !0.91 1.08 0.53

(113) MAE 4.18 4.32 4.42 4.09 4.21 3.06 1.36 0.31 0.41 0.91 1.09 0.58

rms 4.24 4.38 4.47 4.14 4.26 3.10 1.40 0.41 0.53 1.00 1.18 0.68

Total MSE !4.28 !4.40 !4.50 !4.17 !4.29 !3.15 !1.48 !0.24 !0.48 !1.01 0.90 0.35

(131) MAE 4.28 4.40 4.50 4.17 4.29 3.15 1.48 0.40 0.51 1.01 1.00 0.60

rms 4.39 4.50 4.60 4.26 4.39 3.22 1.56 0.63 0.75 1.18 1.12 0.72

!eN+1(N + 1) ! EAvertical

Atoms MSE !2.84 !2.89 !3.02 !2.93 !3.01 !2.24 !1.32 !0.17 !0.25 !0.53 0.08 !0.32

(18) MAE 2.84 2.89 3.02 2.93 3.01 2.24 1.32 0.39 0.39 0.56 0.40 0.35

rms 3.05 3.10 3.22 3.14 3.21 2.40 1.42 0.65 0.66 0.83 0.59 0.56

Molecules MSE !1.91 !1.91 !2.03 !2.18 !1.89 !1.49 !0.88 !0.17 !0.20 !0.29 0.85 0.48

(113) MAE 1.91 1.92 2.03 2.18 1.89 1.49 0.90 0.36 0.34 0.39 0.90 0.58

rms 2.15 2.17 2.25 2.37 2.17 1.64 0.95 0.42 0.43 0.49 1.07 0.77

Total MSE !2.04 !2.05 !2.17 !2.28 !2.04 !1.59 !0.94 !0.17 !0.20 !0.32 0.75 0.37

(131) MAE 2.04 2.05 2.17 2.28 2.04 1.60 0.96 0.36 0.35 0.41 0.83 0.55

rms 2.29 2.32 2.41 2.49 2.34 1.76 1.03 0.46 0.47 0.55 1.02 0.75

!eN+1(N) ! EAvertical

Atoms MSE 3.11 2.62 2.77 2.02 2.47 2.08 0.77 !0.53 !0.36 !0.02 8.00 7.36

(18) MAE 3.11 2.76 2.77 2.14 2.62 2.08 1.00 0.64 0.64 0.75 8.00 7.36

rms 3.54 3.16 3.19 2.54 3.02 2.36 1.14 0.85 0.85 0.89 8.58 7.85

Molecules MSE 2.59 2.40 2.39 2.07 2.16 1.82 1.02 !0.36 !0.30 0.01 7.76 7.16

(113) MAE 2.59 2.40 2.39 2.07 2.16 1.82 1.02 0.47 0.50 0.51 7.76 7.16

rms 2.87 2.65 2.62 2.34 2.46 2.02 1.10 0.54 0.58 0.58 7.98 7.37

Total MSE 2.66 2.43 2.44 2.06 2.20 1.85 0.98 !0.39 !0.31 0.00 7.79 7.19

(131) MAE 2.66 2.45 2.44 2.08 2.22 1.85 1.02 0.49 0.52 0.54 7.79 7.19

rms 2.97 2.72 2.70 2.37 2.55 2.07 1.11 0.59 0.62 0.63 8.06 7.43

(5)

D. Fundamental gaps

The FG of an N-electron system is defined as Eg) IP(N) ! EA(N).

As described in detail in ref. 31, since IP(N) and EA(N) can both be calculated by their definitions, or by their relations to the frontier orbital energies, Egcan be evaluated by eN+1(N) ! eN(N), eN+1(N + 1) ! eN(N) and E(N ! 1) ! 2E(N) + E(N + 1). Here, we examine the performance of the functionals on FGs for the FG131 database using the 6-311++G(3df,3pd) basis set with the EML(75 302) grid and summarize the results in Table 5.

For the calculations using eN+1(N) ! eN(N) (the HOMO–

LUMO gaps), the LSDA, GGAs, MGGAs and AC model potentials perform very poorly due to the lack of the DD. As mentioned previously, the hybrid functionals incorporate part of the DD in the GKS orbital energies^84–88and therefore provide significant improvements, especially for the LC hybrid functionals.

For the calculations using eN+1(N + 1) ! eN(N) (the energy difference between the HOMOs of the N- and (N + 1)-electron systems), the LC hybrid functionals and AC model potentials perform very well due to the correct asymptote of the XC potential. In spite of the incorrect asymptote, there is a fortuitous cancellation of errors in the FGs calculated by the LSDA, GGAs and MGGAs, as the HOMO energies of the N- and (N + 1)-electron systems are both (incorrectly) upshifted.

For the calculations using E(N ! 1) ! 2E(N) + E(N + 1) (the IP(N) ! EA(N) values), the AC model potentials perform very

poorly, as the ground-state energies of the N- and (N * 1)- electron systems are involved in this estimate. All the other functionals perform very well, as the ground-state energies are insensitive to the asymptote.

E. Dissociation of symmetric radical cations

Due to the severe SIEs of semilocal functionals, spurious fractional charge dissociation can occur,⁸⁹ especially for symmetric charged radicals, X2+, such as H2+ 90 and He2+.⁹¹ To examine the performance of the density functional methods upon the SIE problems, we perform spin-unrestricted calculations using the aug-cc-pVQZ basis set and a high-quality EML(250 590) grid. The DFT results are compared with the results from HF theory and the highly accurate CCSD(T) theory.

The dissociation curves of H₂⁺and He₂⁺are shown in Fig. 1 and 2, respectively. As can be seen, the LC hybrid scheme can remove the unphysical barriers of the dissociation curves, while the AC model potential scheme cannot. We emphasize that a SIE-corrected functional must be AC, while an AC functional (or model potential) may not be SIE-corrected.

F. Excitation energies

1. Valence and Rydberg excitations. To assess the performance of the density functionals on valence and Rydberg excitation energies, we perform TDDFT calculations on five small molecules,⁹²including nitrogen gas (N2), carbon monoxide (CO),

Table 5 Statistical errors (in eV) of the fundamental gaps for the FG131 database

HOMO–LUMO gaps

Atoms MSE !7.93 !7.43 !7.70 !6.61 !7.20 !5.73 !2.93 !0.39 !0.74 !1.54 !8.15 !8.02

(18) MAE 7.93 7.43 7.70 6.61 7.20 5.73 2.93 0.77 1.04 1.78 8.15 8.02

rms 8.40 7.95 8.17 7.05 7.67 6.03 3.16 1.04 1.29 2.05 8.62 8.42

Molecules MSE !6.91 !6.86 !6.95 !6.29 !6.51 !5.02 !2.52 0.11 !0.20 !1.05 !6.81 !6.77

(113) MAE 6.91 6.86 6.95 6.29 6.51 5.02 2.52 0.45 0.49 1.06 6.81 6.77

rms 7.08 7.01 7.09 6.44 6.67 5.12 2.57 0.57 0.65 1.28 7.06 6.99

Total MSE !7.05 !6.94 !7.05 !6.34 !6.60 !5.11 !2.57 0.04 !0.27 !1.12 !6.99 !6.94

(131) MAE 7.05 6.94 7.05 6.34 6.60 5.11 2.57 0.49 0.57 1.16 6.99 6.94

rms 7.27 7.15 7.25 6.53 6.81 5.26 2.66 0.65 0.77 1.41 7.29 7.20

eN+1(N + 1) ! eN(N)

Atoms MSE !1.98 !1.92 !1.91 !1.66 !1.73 !1.40 !0.83 !0.75 !0.85 !1.04 !0.22 !0.33

(18) MAE 1.98 1.92 1.91 1.66 1.73 1.40 0.87 0.80 0.90 1.08 0.53 0.51

rms 2.36 2.33 2.28 1.98 2.10 1.66 1.00 1.05 1.12 1.30 0.73 0.67

Molecules MSE !2.41 !2.54 !2.53 !2.05 !2.46 !1.71 !0.62 !0.08 !0.31 !0.76 0.09 !0.09

(113) MAE 2.41 2.54 2.53 2.05 2.46 1.71 0.63 0.39 0.44 0.79 0.49 0.42

rms 2.60 2.73 2.69 2.24 2.67 1.84 0.75 0.49 0.56 0.90 0.68 0.59

Total MSE !2.35 !2.46 !2.44 !2.00 !2.36 !1.67 !0.65 !0.17 !0.38 !0.80 0.05 !0.12

(131) MAE 2.35 2.46 2.44 2.00 2.36 1.67 0.66 0.45 0.50 0.83 0.50 0.43

rms 2.57 2.68 2.64 2.21 2.60 1.82 0.79 0.60 0.66 0.96 0.69 0.61

IP(N) ! EA(N) values

Atoms MSE 0.15 0.22 0.24 0.32 0.35 0.31 0.29 0.37 0.32 0.29 0.45 0.76

(18) MAE 0.35 0.31 0.40 0.39 0.37 0.41 0.34 0.41 0.38 0.37 1.90 1.43

rms 0.58 0.51 0.60 0.62 0.64 0.60 0.51 0.66 0.61 0.59 2.67 2.16

Molecules MSE !0.42 !0.56 !0.57 !0.11 !0.54 !0.31 0.08 0.15 0.11 !0.01 !2.62 !1.60

(113) MAE 0.52 0.62 0.62 0.38 0.61 0.40 0.29 0.33 0.31 0.29 2.86 1.90

rms 0.70 0.80 0.80 0.48 0.76 0.53 0.36 0.40 0.41 0.41 3.35 2.35

Total MSE !0.34 !0.45 !0.46 !0.05 !0.42 !0.22 0.11 0.18 0.14 0.03 !2.20 !1.27

(131) MAE 0.50 0.58 0.59 0.38 0.57 0.40 0.29 0.34 0.32 0.30 2.73 1.83

rms 0.69 0.77 0.78 0.50 0.75 0.54 0.39 0.44 0.45 0.44 3.26 2.32

(6)

water (H2O), ethylene (C2H4) and formaldehyde (CH2O), using the 6-311(2+,2+)G** basis set and EML(99 590) grid. The molecular geometries and experimental values of the excitation energy are taken from ref. 92. For the AC model potential scheme, the adiabatic LSDA kernel is adopted in our TDDFT calculations.³³

As shown in Table 6, all the functionals perform reasonably well for the valence excitations. In contrast, for the Rydberg excitations, the LSDA and GGAs severely underestimate the excitation energies and the MGGAs only provide very minor improvements. On the other hand, the LC hybrid functionals and AC model potentials yield very accurate excitation energies for the Rydberg states owing to the correct asymptote of the XC potential.

2. Long-range charge-transfer excitations. Semilocal functionals qualitatively fail to describe long-range charge-transfer (CT) excitations between a donor and an acceptor.^93–98 The correct CT excitation energy, oCT(R), should have the following asymptote:⁹³

oCT(R - N) E IPD! EAA! 1/R (8) where IPDis the IP of the donor, EAAis the EA of the acceptor and R is their separation. As the (!1/R) dependency, which is

the consequence of Coulomb interactions, can be accurately described by the nonlocal HF exchange,⁹³ LC hybrid functionals, containing the full long-range HF exchange, can yield accurate CT excitation energies.27,28,31,32

On the other hand, within the framework of TDDFT, Hellgren and Gross showed that the discontinuity of the XC kernel is essential for the accurate description of CT excitations.⁹⁸They showed that the divergency of the discontinuity as r - N can produce a divergent kernel in the dissociation limit, which is essential for the asymptote of CT excitation energy (see eqn (8)).^95,98

Following Dreuw et al., we perform TDDFT calculations for the lowest CT excitations between ethylene and tetrafluoro- ethylene, when separated by a distance, R. For the AC model potentials, we adopt the adiabatic LSDA kernel in the TDDFT calculations.³³ Calculations are performed with the 6-31G*

basis and EML(99 590) grid. For comparison, high-level results by the symmetry-adapted-cluster configuration-interaction (SAC-CI) method are taken from Tawada et al.⁹⁹

Fig. 3 shows the values of the CT excitation energies. For the LSDA, GGAs, MGGAs and AC model potentials, the predicted CT excitation energies are severely underestimated due to the lack of the discontinuity of the XC kernel. By incorporating a fraction of nonlocal HF exchange, the hybrid functionals show some improvements, especially for the LC hybrid functionals. However, since the LC hybrid functionals still have short-range SIEs, a more flexible operator for HF exchange may be needed in the LC hybrid functional to further reduce such errors.²⁹

Fig. 4 shows the relative values of the CT excitation energies, where the excitation energy at 5 Å is set to zero for each method.

Therefore, the (!1/R) dependence dominates the relative CT excitation energies. The results of the LSDA, GGAs and MGGAs are overlapped, showing that none of them can capture the long-range Coulomb interaction in the CT excitations. The AC model potentials fail to improve semilocal functionals here.

While the global hybrid functionals provide some improvements, the LC hybrid functionals predict the relative values of the CT excitation curves, which are in excellent agreement with the high-level SAC-CI results.

III. Conclusions

In conclusion, we have examined the performance of the LC hybrid scheme and AC model potential scheme on a very wide range of applications. Due to the correct XC potential in the asymptotic region, both the LC hybrid scheme and AC model potential scheme are reliably accurate for properties related to the asymptote, such as the HOMO energies and Rydberg excitation energies of molecules.

Relative to the AC model potential scheme, the LC hybrid scheme has several advantages. First, the LC hybrid scheme shows significant improvements in the SIE problems, such as dissociation of symmetric radical cations, while the AC model potential scheme fails to resolve the SIE problems. Secondly, as the LC hybrid scheme belongs to the GKS method,

Fig. 1 Dissociation curve of H2+. Zero level is set to E(H) + E(H⁺) for each method.

Fig. 2 Dissociation curve of He₂⁺. Zero level is set to E(He) + E(He⁺) for each method.

(7)

the fundamental gap can be accurately predicted by the HOMO–LUMO gap. For the AC model potential scheme, the HOMO–LUMO gap is a poor approximation of the fundamental gap due to the lack of the DD. Thirdly, the long-range HF exchange kernel in the LC hybrid scheme is essential for the accurate description of the long-range CT excitations between two well-separated molecules, which cannot be described by the AC model potential scheme due to the lack of the discontinuity of the XC kernel. Finally, the LC hybrid scheme yields accurate predictions for atomization energies, barrier heights and noncovalent interactions. However, due to the lack of Exc[r], the AC model potential scheme performs very poorly for total energies and related properties.

Despite their similarity in correct asymptotic behavior, the AC model potential scheme may be accurate only for the asymptote problems, while the LC hybrid scheme, which has remedied several qualitative failures of semilocal functionals, could be reliably accurate for a very wide range of applications. Although only three LC hybrid functionals (oB97, oB97X, and oB97X-D) are examined in this work, we expect that other LC hybrid functionals should also perform well for properties sensitive to the long-range HF exchange, such as asymptote problems, SIE problems, energy- gap problems and charge-transfer problems. For properties insensitive to the long-range HF exchange (e.g., atomization energies, noncovalent interactions), further comparisons should be made for other LC hybrid functionals (e.g., see ref. 31 and 32).

Table 6 The 19 valence and 23 Rydberg excitation energies (in eV) of N2, CO, water, ethylene and formaldehyde, calculated by various methods. The molecular geometries and experimental reference values are taken from ref. 92

Molecule State LSDA PBE BLYP M06L VS98 B3LYP M06-2X oB97 oB97X oB97X-D LB94 LBa Expt.

N2 V¹Pg 9.09 9.12 9.10 9.44 9.42 9.27 9.05 9.49 9.44 9.38 8.72 9.01 9.31

V¹Su! 9.70 9.70 9.60 9.85 10.03 9.35 8.43 9.36 9.30 9.31 9.46 9.53 9.97

V¹Du 10.27 10.12 9.90 10.29 10.16 9.75 10.01 9.85 9.82 9.82 10.06 10.13 10.27

V³Su+ 7.93 7.56 7.50 7.32 7.24 7.11 7.50 7.15 7.15 7.17 7.58 7.65 7.75

V³Pg 7.58 7.42 7.46 7.84 7.74 7.60 7.76 7.98 7.89 7.82 7.22 7.47 8.04

V³Du 8.88 8.37 8.29 8.00 8.44 8.03 8.43 8.30 8.23 8.23 8.60 8.67 8.88

V³Su! 9.70 9.70 9.60 10.29 10.03 9.35 8.90 9.36 9.30 9.31 9.46 9.53 9.67

V³Pu 10.36 10.40 10.32 11.07 10.80 10.64 11.36 11.14 11.05 10.98 10.13 10.26 11.19

C4O V¹P 8.19 8.25 8.24 8.57 8.55 8.40 8.23 8.57 8.52 8.47 7.98 8.16 8.51

V¹S^! 9.89 9.86 9.78 10.26 10.19 9.72 9.12 9.88 9.80 9.78 9.92 9.94 9.88

V³P 5.96 5.74 5.82 6.10 6.03 5.86 6.29 6.15 6.11 6.07 5.60 5.75 6.32

V³S⁺ 8.42 8.11 8.09 8.06 7.93 7.94 8.16 8.00 7.98 8.00 8.41 8.44 8.51

V³D 9.21 8.77 8.72 8.90 8.89 8.67 9.12 8.97 8.89 8.88 9.22 9.25 9.36

V³S^! 9.89 9.86 9.78 10.51 10.19 9.72 9.31 9.88 9.80 9.78 9.92 9.94 9.88

H2O R¹B1 6.49 6.33 6.18 6.90 6.73 6.86 7.33 7.53 7.45 7.23 7.81 7.49 7.4

R¹A2 7.62 7.44 7.26 7.70 7.78 8.23 8.81 9.19 9.09 8.63 9.66 9.28 9.1

R¹A1 8.39 8.18 8.04 8.59 8.53 8.88 9.47 9.63 9.55 9.20 9.81 9.53 9.7

R¹B1 7.99 7.84 7.66 8.19 8.09 8.75 9.58 9.81 9.69 9.17 11.14 10.53 10.0

R¹A1 8.58 8.48 8.31 9.12 8.91 9.12 9.84 10.03 9.91 9.49 11.32 10.71 10.17

R³B1 6.24 6.00 5.89 6.60 6.38 6.49 6.98 7.12 7.09 6.89 7.33 7.06 7.2

C2H4 R¹B3u 6.61 6.40 6.16 6.63 6.63 6.58 6.88 7.57 7.41 7.02 7.77 7.38 7.11

V¹B1u 7.41 7.29 7.08 7.46 7.45 7.35 7.51 7.63 7.60 7.52 7.70 7.66 7.60

R¹B1g 7.08 6.87 6.61 7.01 7.06 7.10 7.39 8.19 8.05 7.59 7.22 7.53 7.80

R¹B2g 7.05 6.83 6.55 7.00 7.03 7.09 7.47 8.33 8.16 7.66 8.46 8.03 8.01

R¹Ag 7.40 7.18 6.92 7.23 7.37 7.44 7.86 8.54 8.38 7.87 9.68 9.13 8.29

R¹B3u 7.59 7.42 7.19 7.68 7.48 7.78 8.27 9.02 8.84 8.36 9.73 9.25 8.62

V³B1u 4.69 4.26 4.31 4.28 4.13 4.07 4.54 3.94 4.01 4.12 4.45 4.49 4.36

R³B3u 6.56 6.31 6.10 6.41 6.52 6.51 6.83 7.44 7.31 6.92 7.61 7.24 6.98

R³B1g 7.00 6.83 6.59 6.88 7.01 7.07 7.38 7.75 7.69 7.50 6.74 7.05 7.79

R³B2g 7.03 6.78 6.52 6.84 6.97 7.04 7.41 8.19 8.04 7.56 8.24 7.93 7.79

R³Ag 7.36 7.06 6.86 6.74 7.19 7.36 7.66 8.23 8.07 7.63 9.44 8.93 8.15

CH2O V¹A2 3.62 3.74 3.76 4.22 3.99 3.85 3.58 3.91 3.88 3.88 3.49 3.65 4.07

R¹B2 5.89 5.77 5.63 6.25 6.10 6.47 7.15 7.43 7.35 6.96 7.43 7.13 7.11

R¹B2 6.67 6.54 6.40 6.85 6.90 7.25 7.85 8.15 8.08 7.66 8.61 8.37 7.97

R¹A1 7.22 7.11 6.94 7.42 7.37 7.98 9.01 9.30 9.21 8.74 9.58 9.51 8.14

R¹A2 6.81 6.69 6.54 6.96 7.01 7.47 8.13 8.41 8.34 7.84 9.03 9.07 8.37

R¹B2 6.90 6.81 6.63 7.20 7.00 7.71 8.91 9.15 8.99 8.52 10.00 9.67 8.88

V³A2 3.01 3.00 3.06 3.56 3.29 3.13 3.05 3.26 3.22 3.21 2.86 3.03 3.50

V³A1 6.03 5.57 5.59 5.49 5.39 5.22 5.49 5.27 5.25 5.29 5.87 5.89 5.86

R³B2 5.84 5.63 5.55 5.96 5.92 6.36 7.04 7.24 7.19 6.81 7.14 6.91 6.83

R³B2 6.65 6.47 6.36 6.66 6.82 7.17 7.69 7.91 7.87 7.50 8.43 8.17 7.79

R³A1 6.53 6.35 6.22 6.47 6.60 7.16 7.86 8.07 8.00 7.56 8.47 8.17 7.96

MAE Valence (19) 0.24 0.32 0.36 0.29 0.29 0.42 0.41 0.28 0.31 0.32 0.36 0.27

Rydberg (23) 1.12 1.30 1.48 1.04 1.03 0.75 0.29 0.26 0.20 0.35 0.73 0.42

(8)

Acknowledgements

This work was supported by the National Science Council of Taiwan (Grant No. NSC101-2112-M-002-017-MY3), National Taiwan Univer- sity (Grant Nos. 99R70304, 101R891401 and 101R891403) and the National Center for Theoretical Sciences of Taiwan. Y.C.S. was partially supported by the NSC Undergraduate Fellowship. We are grateful for the partial support of the computer resources from the groups of Dr J.-L. Kuo (Academia Sinica) and Dr Y.-C. Cheng (NTU).

References

1 P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864.

2 W. Kohn and L. J. Sham, Phys. Rev. A, 1965, 140, A1133.

3 R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989.

4 S. Ku¨mmel and L. Kronik, Rev. Mod. Phys., 2008, 80, 3.

5 J. P. Perdew, A. Ruzsinszky, L. A. Constantin, J. Sun and G. I. Csonka, J. Chem. Theory Comput., 2009, 5, 902.

6 E. Runge and E. K. U. Gross, Phys. Rev. Lett., 1984, 52, 997.

7 M. Petersilka, U. J. Gossmann and E. K. U. Gross, Phys. Rev.

Lett., 1996, 76, 1212.

8 C. A. Ullrich, Time-Dependent Density-Functional Theory:

Concepts and Applications, Oxford University Press, New York, 2012.

9 M. E. Casida, Recent Advances in Density Functional Methods, Part I, World Scientific, Singapore, 1995.

10 E. K. U. Gross, J. F. Dobson and M. Petersilka, Density Functional Theory II, Springer, Heidelberg, 1996.

11 J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin and J. Sun, Phys. Rev. Lett., 2009, 103, 026403.

12 A. J. Cohen, P. Mori-Sa´nchez and W. Yang, Science, 2008, 321, 792.

13 A. J. Cohen, P. Mori-Sa´nchez and W. Yang, Chem. Rev., 2011, 112, 289.

14 J.-D. Chai, J. Chem. Phys., 2012, 136, 154104.

15 M. Levy, J. P. Perdew and V. Sahni, Phys. Rev. A: At., Mol., Opt. Phys., 1984, 30, 2745.

16 C.-O. Almbladh and U. von Barth, Phys. Rev. B, 1985, 31, 3231.

17 A. Go¨rling, Phys. Rev. Lett., 1999, 83, 5459.

18 P. W. Ayers and M. Levy, J. Chem. Phys., 2001, 115, 4438.

19 J. P. Perdew and A. Zunger, Phys. Rev. B, 1981, 23, 5048.

20 M. E. Casida, C. Jamorski, K. C. Casida and D. R. Salahub, J. Chem. Phys., 1998, 108, 4439.

21 D. J. Tozer and N. C. Handy, J. Chem. Phys., 1998, 109, 10180.

22 A. Savin, in Recent Developments and Applications of Modern Density Functional Theory, ed. J. M. Seminario, Elsevier, Amsterdam, 1996, pp. 327–357.

23 H. Iikura, T. Tsuneda, T. Yanai and K. Hirao, J. Chem. Phys., 2001, 115, 3540.

24 T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51.

25 O. A. Vydrov, J. Heyd, A. V. Krukau and G. E. Scuseria, J. Chem. Phys., 2006, 125, 074106.

26 E. Livshits and R. Baer, Phys. Chem. Chem. Phys., 2007, 9, 2932.

27 J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2008, 128, 084106.

28 J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615.

29 J.-D. Chai and M. Head-Gordon, Chem. Phys. Lett., 2008, 467, 176.

30 J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2009, 131, 174105.

31 Y.-S. Lin, C.-W. Tsai, G.-D. Li and J.-D. Chai, J. Chem. Phys., 2012, 136, 154109.

32 Y.-S. Lin, G.-D. Li, S.-P. Mao and J.-D. Chai, J. Chem. Theory Comput., 2013, 9, 263.

33 R. van Leeuwen and E. J. Baerends, Phys. Rev. A: At., Mol., Opt. Phys., 1994, 49, 2421.

34 D. J. Tozer, J. Chem. Phys., 2000, 112, 3507.

35 P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen and E. J. Baerends, J. Chem. Phys., 2000, 112, 1344.

Fig. 3 The lowest CT excitation energy of a C2H4 ( ( ( C2F4 dimer along the intermolecular distance R (in Å).

Fig. 4 Relative excitation energy for the lowest CT excitation of a C2H4( ( ( C2F4

dimer along the intermolecular distance R (in Å). The excitation energy at 5 Å is set to zero for each method.