Thermally-assisted-occupation density functional theory with generalized-gradient approximations

(1)

J. Chem. Phys. 140, 18A521 (2014); https://doi.org/10.1063/1.4867532 140, 18A521

Thermally-assisted-occupation density

functional theory with generalized-gradient approximations

Cite as: J. Chem. Phys. 140, 18A521 (2014); https://doi.org/10.1063/1.4867532

Submitted: 31 December 2013 . Accepted: 21 February 2014 . Published Online: 13 March 2014 Jeng-Da Chai

ARTICLES YOU MAY BE INTERESTED IN

Density functional theory with fractional orbital occupations

The Journal of Chemical Physics 136, 154104 (2012); https://doi.org/10.1063/1.3703894

Role of exact exchange in thermally-assisted-occupation density functional theory: A proposal of new hybrid schemes

A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu

(2)

THE JOURNAL OF CHEMICAL PHYSICS 140, 18A521 (2014)

Thermally-assisted-occupation density functional theory with generalized-gradient approximations

Jeng-Da Chai^a)

Department of Physics, Center for Theoretical Sciences, and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan

(Received 31 December 2013; accepted 21 February 2014; published online 13 March 2014)

We extend the recently proposed thermally-assisted-occupation density functional theory (TAO- DFT) [J.-D. Chai, J. Chem. Phys. 136, 154104 (2012)] to generalized-gradient approximation (GGA) exchange-correlation density functionals. Relative to our previous TAO-LDA (i.e., the local density approximation to TAO-DFT), the resulting TAO-GGAs are significantly superior for a wide range of applications, such as thermochemistry, kinetics, and reaction energies. For noncovalent interactions, TAO-GGAs with empirical dispersion corrections are shown to yield excellent performance. Due to their computational efficiency for systems with strong static correlation effects, TAO-LDA and TAO-GGAs are applied to study the electronic properties (e.g., the singlet-triplet energy gaps, vertical ionization potentials, vertical electron affinities, fundamental gaps, and symmetrized von Neumann entropy) of acenes with different number of linearly fused benzene rings (up to 100), which is very challenging for conventional electronic structure methods. The ground states of acenes are shown to be singlets for all the chain lengths studied here. With the increase of acene length, the singlet-triplet energy gaps, vertical ionization potentials, and fundamental gaps decrease monotonically, while the vertical electron affinities and symmetrized von Neumann en- tropy (i.e., a measure of polyradical character) increase monotonically. © 2014 AIP Publishing LLC.

[http://dx.doi.org/10.1063/1.4867532]

I. INTRODUCTION

Due to its favorable balance between cost and performance, Kohn-Sham density functional theory (KS-DFT)^1,2 has been a very popular method for the study of large ground- state systems.^3–5 However, its essential ingredient, the exact exchange-correlation (XC) density functional Exc[ρ], remains unknown and needs to be approximated. Functionals based on the conventional density functional approximations (DFAs), such as the local density approximation (LDA) and generalized gradient approximations (GGAs), are reliably accurate for properties governed by short-range XC effects, and are computationally efficient for very large systems. However, KS-DFAs can produce erroneous results in situations where the nonlocal XC effects are pronounced. Over the past two decades, numerous attempts have been made to resolve the qualitative failures of KS-DFAs.^4–8

Recently, we have shown that long-range corrected (LC) hybrid functionals,^9–20 incorporating the long-range Hartree- Fock (HF) exchange into the KS-DFAs, could be reliably accurate for a very wide range of applications,²¹ especially for those sensitive to the long-range HF exchange, such as the asymptote problems,^22–29self-interaction-error problems,^30,31 energy-gap problems,^32–45and charge-transfer problems.^46–51 To reduce the qualitative failures of KS-DFAs for noncovalent interactions,⁵²the DFT-D (KS-DFT with empirical dispersion corrections) schemes,^15,53–55adding empirical atom- atom dispersion potentials into the KS-DFAs, have shown

a)Electronic mail: [email protected]

generally satisfactory performance on a large set of noncovalent systems.^56,57 Alternatively, the double-hybrid (DH) methods,^17,58–65 mixing some of the HF exchange and some of the nonlocal orbital correlation energy from the second- order Møller-Plesset⁶⁶perturbation (MP2) theory into the KS- DFAs, can also be adopted to take into account nonlocal dy- namical correlation effects. DH functionals have shown an overall satisfactory accuracy for thermochemistry, kinetics, noncovalent interactions, and self-interaction-error problems.

Despite their computational efficiency, KS-DFAs, global hybrid functionals,⁶⁷ LC hybrid functionals, and DH functionals perform very poorly for strongly correlated (SC) systems (i.e., multi-reference systems), which are systems with strong static correlation effects, including bond- breaking reactions, conjugated polymers, and transition-metal compounds.^6–8,68 Within the framework of KS-DFT, fully nonlocal XC functionals, such as those based on random phase approximation (RPA), can be essential for the accurate treatment of SC systems. However, RPA-type functionals are computationally very demanding for large systems.^4,5,69,70

Aiming to reduce the static-correlation-error problems with minimum computational complexity, we have recently developed thermally-assisted-occupation DFT (TAO-DFT),⁸ a DFT with fractional orbital occupations given by the Fermi- Dirac distribution (controlled by a fictitious temperature θ ), for the study of large systems with strong static correlation effects. In contrast to finite-temperature DFT,⁷¹ TAO-DFT is developed for ground-state systems at zero temperature. TAO- DFT has similar computational cost as KS-DFT, and is re- duced to KS-DFT in the absence of strong static correlation

(3)

effects. Even at the simplest LDA level, the resulting TAO- LDA has been shown to consistently improve upon KS-LDA for multi-reference systems. However, TAO-LDA performs similarly to KS-LDA for single-reference systems, due to the absence of strong static correlation.

To improve the performance of TAO-LDA for single- reference systems, here we propose TAO-GGAs for the improved description of short-range XC effects. Relative to TAO-LDA, the resulting TAO-GGAs are significantly superior for a wide range of applications, such as thermochemistry, kinetics, and reaction energies. For noncovalent interactions, TAO-GGAs with empirical dispersion corrections are shown to yield excellent performance. The rest of this paper is or- ganized as follows. In Sec.II, we briefly describe the formulation of TAO-DFT and the DFAs to TAO-DFT. The performance of TAO-LDA and TAO-GGAs is compared with that of KS-LDA and KS-GGAs in Sec.III. In Sec.IV, we apply TAO-LDA and TAO-GGAs to study the electronic properties of linear acenes. Our conclusions are given in Sec.V.

II. TAO-DFT

For a system with Nα up-spin electrons and Nβ down- spin electrons in an external potential vext(r) at zero temper- ature, two noninteracting auxiliary systems at the same ficti- tious temperature θ (measured in energy units) are adopted in spin-polarized (spin-unrestricted) TAO-DFT:⁸ one described by the spin function α and the other by function β, with the re- spective thermal equilibrium density distributions ρ_{s, α}(r) and ρ_{s, β}(r) exactly equal to ρα(r) and ρβ(r), respectively, in the original interacting system at zero temperature. The resulting self-consistent equations for σ -spin electrons (σ = α or β) are given by (i runs for the orbital index)

− ¯²

2me∇² + vs,σ(r)

ψ_i,σ(r)= i,σψ_i,σ(r), (1) where the effective potential is

vs,σ(r)= vext(r)+ e²

ρ(r)

|r − r|dr+δExc[ρα, ρβ] δρ_σ(r) +δE_θ[ρα, ρ_β]

δρσ(r) . (2)

Here Exc[ρα, ρβ] is the XC energy defined in spin-polarized KS-DFT,^72,73 and Eθ[ρα, ρ_β]≡ A^θs⁼⁰[ρα, ρ_β]− A^θs[ρα, ρ_β] is the difference between the noninteracting kinetic free energy at zero temperature and that at the fictitious temperature θ. The σ -spin density is given by

ρσ(r)=^∞

i=1

fi,σ|ψi,σ(r)|², (3) where the occupation number fi, σis the Fermi-Dirac function f_i,σ = {1 + exp[(i,σ − μσ)/θ ]}⁻¹, (4) and the chemical potential μσ is chosen to conserve the num- ber of σ -spin electrons Nσ,

∞ i=1

{1 + exp[(i,σ − μσ)/θ ]}⁻¹= Nσ. (5)

The ground-state density ρ(r) is computed as the sum of the up-spin density ρα(r) and down-spin density ρβ(r):

ρ(r)= ρα(r)+ ρβ(r)=

σ=α,β

ρσ(r). (6)

The formulation of spin-polarized TAO-DFT has yielded two sets (one for each spin function) of self-consistent equations, Eqs.(1)–(5), for ρα(r) and ρβ(r), respectively, which are cou- pled with ρ(r) by Eq.(6).

The self-consistent procedure described in Ref.8may be adopted to obtain the spin densities and ground-state density.

After self-consistency is achieved, the noninteracting kinetic free energy A^θ_s is given by

A^θ_s[{fi,α, ψi,α}, {fi,β, ψi,β}] = Ts^θ[{fi,α, ψi,α}, {fi,β, ψi,β}]

+ES^θ[{fi,α}, {fi,β}], (7) which is the sum of the kinetic energy

T_s^θ[{fi,α, ψ_i,α}, {fi,β, ψ_i,β}]

= − ¯² 2me

σ=α,β

∞ i=1

f_i,σ

ψ_i,σ^∗ (r)∇²ψ_i,σ(r)dr

=

σ=α,β

∞

i=1

fi,σi,σ−

ρσ(r)vs,σ(r)dr

(8)

and entropy contribution E^θ_S[{fi,α}, {fi,β}] = θ

σ=α,β

∞ i=1

{fi,σ ln(fi,σ)+ (1 − fi,σ)

× ln(1 − fi,σ)} (9)

of noninteracting electrons at the fictitious temperature θ . The total ground-state energy E[ρα, ρβ] in spin-polarized TAO- DFT is evaluated by

E[ρα, ρβ]= A^θs[{fi,α, ψi,α}, {fi,β, ψi,β}]

+

ρ(r)vext(r)dr+ EH[ρ]+ Exc[ρα, ρβ]

+Eθ[ρα, ρ_β], (10)

where EH[ρ]≡^e₂² _ρ(r)ρ(r)

|r−r| drdr is the Hartree energy.

Spin-unpolarized (spin-restricted) TAO-DFT can be formu- lated by imposing the constraints of ψi, α(r)= ψi, β(r) and f_{i, α}= fi, βto spin-polarized TAO-DFT.

In spin-polarized TAO-DFT, as the exact Exc[ρα, ρβ] and Eθ[ρα, ρβ], in terms of the spin densities ρα(r) and ρβ(r), have not been known, DFAs for both of them (denoted as TAO-DFAs) are needed for practical applications. Accord- ingly, the performance of TAO-DFAs depends on the accu- racy of DFAs and the choice of the fictitious temperature θ . While E_xc^DFA[ρα, ρβ] can be readily obtained from those in spin-polarized KS-DFT, E_θ^DFA[ρα, ρβ] can be obtained with the knowledge of A^DFA,θ_s [ρα, ρβ] as follows:

E^DFA_θ [ρα, ρβ]≡ A^DFA,θs ⁼⁰[ρα, ρβ]− A^DFA,θs [ρα, ρβ]. (11) Note that E_θ^DFA₌₀[ρα, ρ_β]= 0 (i.e., an exact property of Eθ[ρα, ρ_β]) is ensured by Eq.(11). From the spin-scaling relation of

(4)

18A521-3 Jeng-Da Chai J. Chem. Phys. 140, 18A521 (2014)

A^θ_s[ρα, ρβ],⁷⁴Eq.(11)can be expressed in terms of A^DFA,θ_s [ρ]

(in its spin-unpolarized form):

E_θ^DFA[ρα, ρβ]= 1 2

A^DFA,θ_s ⁼⁰[2ρα]+ A^DFA,θs ⁼⁰[2ρβ]

−1 2

A^DFA,θ_s [2ρα]+ A^DFA,θs [2ρβ] . (12)

In our previous work,⁸ Perrot’s parametrization of A^LDA,θ_s [ρ], which is the LDA for A^θ_s[ρ] (see Appendix A of Ref.75), was adopted to obtain E_θ^LDA[ρα, ρ_β]. To go beyond the simple LDA, A^GEA,θ_s [ρ], which is the gradient expansion approximation (GEA) for A^θ_s[ρ] (see Appendices A and B of Ref.75), can be adopted to obtain E_θ^GEA[ρα, ρβ]. For the nearly uniform electron gas, E_θ^GEA[ρα, ρβ] is expected to im- prove upon E_θ^LDA[ρα, ρβ].

As discussed in Ref.8, TAO-DFT offers an explicit description of strong static correlation via the entropy contri- bution E_S^θ[{fi,α}, {fi,β}] (see Eq. (9)). Even at the simplest LDA level, the resulting TAO-LDA has been shown to perform reasonably well for multi-reference systems (due to the appropriate treatment of static correlation), when the distri- bution of orbital occupation numbers {f_{i, σ}} (related to the chosen θ ) is close to that of the natural orbital occupation numbers (NOONs).⁷⁶ This implies that a system-dependent θ (related to the distributions of NOONs) should be needed to capture the essential physics of strong static correlation effects. However, for the sake of simplicity, an optimal value of θ = 7 mhartree was previously defined for TAO-LDA, based on physical arguments and numerical investigations. Interest- ingly, TAO-LDA (with θ = 7 mhartree) was shown to consistently improve upon KS-LDA for multi-reference systems, while performing similarly to KS-LDA for single-reference systems. As TAO-GGAs should improve upon TAO-LDA mainly for properties governed by short-range XC effects, the optimal values of θ for TAO-GGAs should be similar to that for TAO-LDA, when the same physical arguments and nu- merical investigations are employed to define the optimal θ . Therefore, we adopt an optimal value of θ = 7 mhartree for all the TAO-LDA and TAO-GGAs calculations in this work, unless noted otherwise. The limiting case where θ = 0 for TAO-DFA is especially interesting, as this reduces to KS- DFA. Therefore, it is important to know how well KS-DFA performs here to assess the significance of the extra parame- ter θ for TAO-DFA.

III. RESULTS AND DISCUSSION FOR THE TEST SETS For a comprehensive comparison, we examine the performance of various XC functionals: LDA⁷⁷ and three popular GGAs (Perdew-Burke-Ernzerhof (PBE),⁷⁸ Becke-Lee- Yang-Parr (BLYP),⁷⁹ and BLYP-D⁵⁴) in both KS-DFT and TAO-DFT, 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 reaction energies of the 30 chemical reactions in the NHTBH38/04 and HTBH38/04

sets,⁸²the 166 optimized geometry properties of EXTS,⁸⁴and two dissociation curves for H2and N2. There are in total 592 pieces of data in our test sets, which are very large and di- verse. Detailed information about the test sets may be found in Ref.14.

For all the TAO-DFT calculations, θ = 7 mhartree is adopted, unless noted otherwise. Our preliminary TAO-DFT results show that the difference between using E_θ^LDA and E_θ^GEAis much smaller than the difference between using two different XC functionals. Unsurprisingly, as E_θ^LDA₌₀ = E^GEA_θ₌₀

= 0, the difference between Eθ^LDA and E_θ^GEA should remain small for a sufficiently small θ (i.e., 7 mhartree). Therefore, for brevity, we only adopt E^LDA_θ for the TAO-DFT calculations in this work.

All calculations are performed with a development ver- sion of Q-Chem 4.0.⁸⁵ Spin-restricted theory is used for singlet states and spin-unrestricted theory for others, unless noted otherwise. Results for the test sets 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.⁸⁷For the interaction energies of the weakly bound systems, the counterpoise correction⁸⁸is employed to reduce the 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), root-mean-square (rms) errors, maximum negative errors (Max(−)), and maximum positive errors (Max(+)).

A. ωB97 training set

The ωB97 training set¹⁴ contains several well-known 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 I summarizes the statistical errors of various functionals in both KS-DFT and TAO-DFT for the ωB97 training set. As shown, TAO- DFAs perform comparably to the corresponding KS-DFAs.

Unsurprisingly, as these systems do not have much static correlation, the exact NOONs should be close to either 0 or 1, which can be well simulated by the orbital occupation num- bers of TAO-DFAs (with a sufficiently small θ = 7 mhartree).

Due to the improved treatment of short-range XC effects, TAO-GGAs significantly outperform TAO-LDA for the AEs of the G3/99 set, the PAs of the G2-1 set, and the BHs of the NHTBH38/04 and HTBH38/04 sets.⁸⁹For the IPs and EAs of the G2-1 set, TAO-GGAs perform slightly better than TAO- LDA. For the noncovalent interactions of the S22 set, the dispersion corrected functionals (i.e., KS-BLYP-D and TAO- BLYP-D) are found to be very accurate, while all the other functionals perform poorly.⁸⁹This suggests that for noncovalent interactions, the successful DFT-D schemes for KS-DFT remain very accurate for TAO-DFT.

B. Reaction energies

The reaction energies of the 30 chemical reactions (those with different barrier heights for the forward and backward

(5)

TABLE I. Statistical errors (in kcal/mol) of the ωB97 training set.¹⁴

KS-DFT TAO-DFT

System Error LDA PBE BLYP BLYP-D LDA PBE BLYP BLYP-D

G3/99 MSE 120.60 20.90 − 4.59 − 0.83 95.02 7.91 − 16.24 − 12.27

(223) MAE 120.60 21.51 9.76 7.03 95.04 11.41 19.01 15.33

rms 142.51 26.30 12.96 9.17 114.19 15.07 24.24 19.35

IP MSE 3.42 0.03 − 1.50 − 1.50 1.79 − 1.08 − 2.61 − 2.61

(40) MAE 5.54 3.46 4.43 4.44 6.18 4.86 6.10 6.10

rms 6.66 4.35 5.28 5.29 7.63 6.00 7.40 7.40

EA MSE 6.45 1.72 0.36 0.36 4.20 0.22 − 1.08 − 1.07

(25) MAE 6.45 2.42 2.57 2.57 5.49 2.88 4.38 4.40

rms 7.29 3.06 3.17 3.17 6.45 3.44 5.44 5.47

PA MSE − 5.91 − 0.83 − 1.47 − 1.09 − 5.66 − 0.58 − 1.22 − 0.84

(8) MAE 5.91 1.60 1.58 1.55 5.66 1.47 1.50 1.55

rms 6.40 1.91 2.10 1.98 6.16 1.80 1.94 1.86

NHTBH MSE − 12.41 − 8.52 − 8.69 − 9.32 − 11.93 − 8.38 − 8.52 − 9.15

(38) MAE 12.62 8.62 8.72 9.35 12.15 8.49 8.56 9.19

rms 16.13 10.61 10.27 10.83 15.09 10.28 9.90 10.46

HTBH MSE − 17.90 − 9.67 − 7.84 − 8.89 − 16.34 − 9.20 − 7.25 − 8.33

(38) MAE 17.90 9.67 7.84 8.89 16.34 9.20 7.29 8.34

rms 18.92 10.37 8.66 9.52 17.06 9.87 8.24 9.17

S22 MSE − 1.97 2.77 5.05 0.23 − 2.30 2.44 4.70 − 0.12

(22) MAE 2.08 2.77 5.05 0.33 2.33 2.44 4.70 0.28

rms 3.18 3.89 6.31 0.45 3.40 3.57 5.95 0.37

Total MSE 65.86 10.33 − 4.07 − 2.37 51.26 2.81 − 10.81 − 8.99

(394) MAE 72.41 14.63 8.05 6.40 57.76 9.01 13.48 11.31

rms 107.53 20.40 10.88 8.44 86.26 12.38 18.92 15.43

directions) taken from the NHTBH38/04 and HTBH38/04 sets, are used to examine the performance of KS-DFAs and TAO-DFAs. As shown in TableII, TAO-DFAs have very similar performance to the corresponding KS-DFAs. TAO-GGAs are shown to significantly improve upon TAO-LDA for this test set.⁸⁹

C. Equilibrium geometries

Satisfactory predictions of molecular geometries are essential for practical applications. Geometry optimizations for various functionals in both KS-DFT and TAO-DFT are performed on the equilibrium experimental test set (EXTS),⁸⁴ consisting of 166 symmetry unique experimental bond lengths for small to medium size molecules. As the ground states of these molecules at their equilibrium geometries can be well described by single-reference wave func-

tions, TAO-DFAs are found to perform very similarly (see TableIII) to the corresponding KS-DFAs.⁸⁹

D. Dissociation of H2and N2

Due to the presence of strong static correlation effects, the dissociation of H2and N2remains an important and challenging subject in KS-DFT. Based on the symmetry con- straint, the difference between the spin-restricted and spin- unrestricted dissociation limits calculated by an approximate method can be adopted as a quantitative measure of the static correlation error (SCE) of the method.^6,8,68 Spin-restricted KS-DFAs, global hybrid functionals, LC hybrid functionals, and DH functionals have been shown to perform very poorly for both H₂and N₂dissociation curves, leading to very large SCEs. By contrast, as discussed in Ref. 8, spin-restricted TAO-LDA (with a θ between 30 and 50 mhartree) can prop- erly dissociate H2 and N2 (yielding vanishingly small SCEs)

TABLE II. Statistical errors (in kcal/mol) of the reaction energies of the 30 chemical reactions in the NHTBH38/04 and HTBH38/04 sets.⁸²

KS-DFT TAO-DFT

LDA PBE BLYP BLYP-D LDA PBE BLYP BLYP-D

MSE − 0.41 1.08 0.80 0.74 − 1.32 0.23 − 0.12 − 0.20

MAE 8.51 4.39 3.23 3.02 7.09 3.97 3.80 3.67

rms 11.10 6.24 4.37 4.20 9.38 5.97 4.95 4.89

Max(−) − 18.31 − 7.86 − 7.24 − 7.28 − 15.92 − 8.89 − 11.24 − 11.71

Max(+) 35.68 22.59 11.96 12.03 30.50 21.60 10.65 10.73

(6)

TABLE III. Statistical errors (in Å) of EXTS.⁸⁴

KS-DFT TAO-DFT

LDA PBE BLYP BLYP-D LDA PBE BLYP BLYP-D

MSE 0.004 0.014 0.018 0.019 0.005 0.014 0.019 0.020

MAE 0.013 0.015 0.019 0.020 0.013 0.015 0.020 0.021

rms 0.017 0.019 0.024 0.025 0.017 0.020 0.025 0.026

Max(−) − 0.091 − 0.069 − 0.064 − 0.063 − 0.091 − 0.069 − 0.064 − 0.064

Max(+) 0.081 0.095 0.103 0.103 0.080 0.095 0.102 0.102

to the respective spin-unrestricted dissociation limits, which is closely related to that the distribution of orbital occupation numbers {f_{i, σ}} (related to the chosen θ ) matches reasonably well with that of the NOONs.

To examine the performance of the present method, the potential energy curves (in relative energy) for the ground state of H2, calculated by spin-restricted TAO-DFAs (with various θ ), are shown in Figs. 1–4, where the zeros of energy are set at the respective spin-unrestricted dissociation limits. Near the equilibrium geometry of H2, where the single-reference character is dominant, TAO-DFAs (with θ = 7 mhartree) perform very similarly to the correspond- ing KS-DFAs (the θ = 0 cases). At the dissociation limit, where the multi-reference character becomes pronounced, TAO-DFAs (with θ = 40 mhartree) have vanishingly small SCEs, while TAO-DFAs (with θ= 7 mhartree) and KS-DFAs have noticeable SCEs. Overall, the SCEs of TAO-DFAs are not sensitive to the choice of DFAs (mainly responsible for short-range XC effects), but sensitive to the value of θ (closely related to the distribution of the NOONs). Similar results are also found for N2dissociation, as shown in Figs.5–8.

IV. ELECTRONIC PROPERTIES OF LINEAR ACENES Linear n-acenes (C4n+ 2H2n+ 4), consisting of n linearly fused benzene rings (see Fig.9), have recently attracted sig- nificant attention from many experimental and theoretical re-

−120

−100

−80

−60

−40

−20 0 20 40 60

1 2 3 4 5 6 7 8 9 10

Relative Energy (kcal/mol)

Internuclear Distance (Å )

( in mhartree)

=0

=7

=40

FIG. 1. Potential energy curves (in relative energy) for the ground state of H2, calculated by spin-restricted TAO-LDA (with various θ ). The zeros of energy are set at the respective spin-unrestricted dissociation limits. The θ

= 0 case corresponds to spin-restricted KS-LDA.

−120

−100

−80

−60

−40

−20 0 20 40 60

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 2. Same as Fig.1, but for spin-restricted TAO-PBE (with various θ ).

The θ= 0 case corresponds to spin-restricted KS-PBE.

searchers due to their fascinating electronic properties, open- ing up tremendous possibilities to realize electronic and spin- tronic nanodevices.^8,90–103The electronic properties of acenes largely depend on the chain lengths. Despite increasing interest in acenes, it remains very challenging to explore the properties of long-chain acenes from both experimental and theoretical perspectives. Experimentally, the difficulties in synthesizing long-chain acenes and their instability follow- ing isolation have been attributed to their radical character.

Accordingly, the experimental singlet-triplet energy gaps (ST

−120

−100

−80

−60

−40

−20 0 20 40 60

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 3. Same as Fig.1, but for spin-restricted TAO-BLYP (with various θ ).

The θ= 0 case corresponds to spin-restricted KS-BLYP.

(7)

−120

−100

−80

−60

−40

−20 0 20 40 60

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 4. Same as Fig.1, but for spin-restricted TAO-BLYP-D (with various θ). The θ= 0 case corresponds to spin-restricted KS-BLYP-D.

−300

−250

−200

−150

−100

−50 0 50 100 150

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 5. Potential energy curves (in relative energy) for the ground state of N2, calculated by spin-restricted TAO-LDA (with various θ ). The zeros of energy are set at the respective spin-unrestricted dissociation limits. The θ= 0 case corresponds to spin-restricted KS-LDA.

−300

−250

−200

−150

−100

−50 0 50 100 150

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 6. Same as Fig.5, but for spin-restricted TAO-PBE (with various θ ).

The θ= 0 case corresponds to spin-restricted KS-PBE.

−300

−250

−200

−150

−100

−50 0 50 100 150

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 7. Same as Fig.5, but for spin-restricted TAO-BLYP (with various θ ).

The θ= 0 case corresponds to spin-restricted KS-BLYP.

gaps) of n-acenes are only available up to pentacene.^90–93On the theoretical side, acenes, which belong to conjugated π - orbital systems, typically require high-level ab initio multi- reference methods, such as the density matrix renormalization group (DMRG) algorithm,⁹⁶ the variational two-electron re- duced density matrix (2-RDM) method,¹⁰²or other high-level methods,^95,99–101to capture the essential strong static correlation effects. However, these methods are prohibitively expen- sive for the study of large acenes.

As mentioned previously, for systems with pronounced strong static correlation effects, such as large acenes, the predictions from KS-DFAs can be problematic.^6–8 To examine how TAO-DFAs improve upon KS-DFAs here, spin- unrestricted KS-DFT and TAO-DFT calculations, employing various XC functionals: LDA, PBE, BLYP, and BLYP-D, are performed, using the 6-31G(d) basis set (up to 100-acene), for the lowest singlet and triplet energies on the respective geometries that were fully optimized at the same level. The ST gap of n-acene is calculated as (ET− ES), the energy difference between the lowest triplet (T) and singlet (S) states of n-acene.

−300

−250

−200

−150

−100

−50 0 50 100 150

1 2 3 4 5 6 7 8 9 10

( in mhartree)

=0

=7

=40

FIG. 8. Same as Fig.5, but for spin-restricted TAO-BLYP-D (with various θ). The θ= 0 case corresponds to spin-restricted KS-BLYP-D.

(8)

FIG. 9. Hexacene, consisting of 6 linearly fused benzene rings, is designated as 6-acene.

As shown in Figs.10and11, in contrast to the DMRG results,⁹⁶ the ST gaps calculated by spin-unrestricted KS- DFT, unexpectedly increase beyond 10-acene, due to unphysical symmetry-breaking effects.⁸⁹ By contrast, the ST gaps calculated by spin-unrestricted TAO-DFT, which are in good agreement with the existing experimental and high- level ab initio data,^96,101 are shown to decrease monotonically with the increase of chain length. This shows that the ground states of acenes are singlets for all the chain lengths studied. The ST gap of the largest acene studied here (100-acene) is 0.38 kcal/mol for TAO-PBE and TAO-BLYP, 0.39 kcal/mol for TAO-BLYP-D, and 0.49 kcal/mol for TAO- LDA.⁸⁹ To examine the possible symmetry-breaking effects, spin-restricted TAO-DFT calculations are also performed for the lowest singlet energies on the respective geometries that were fully optimized at the same level. The spin-unrestricted and spin-restricted TAO-DFT calculations are found to essen- tially yield the same energy value for the lowest singlet state of n-acene (i.e., no unphysical symmetry-breaking effects).

At the optimized geometry of the lowest singlet state (i.e., the ground state) of n-acene, containing N electrons, the ver- tical ionization potential IPv and the vertical electron affinity EAvcan be calculated by

IPv= EN−1− EN (13)

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70

2 4 6 8 10 12

Singlet-Triplet Energy Gap (kcal/mol)

n-acene

Experiment DMRG/cc-pVDZ DMRG/STO-3G CCSD(T)/CBS KS-LDA KS-PBE KS-BLYP KS-BLYP-D TAO-LDA TAO-PBE TAO-BLYP TAO-BLYP-D

FIG. 10. Singlet-triplet energy gap as a function of the acene length, calculated by various functionals in spin-unrestricted KS-DFT and TAO-DFT. The experimental data (uncorrected for zero-point vibrations, thermal vibrations, etc.) are taken from Refs.90–93, the DMRG data are taken from Ref.96, and the CCSD(T)/CBS data, calculated by the coupled-cluster theory with itera- tive singles and doubles and perturbative treatment of triple substitutions at the complete basis set limit, are taken from Ref.101.

0 1 2 3 4 5 6 7 8

10 20 30 40 50 60 70 80 90 100

Singlet-Triplet Energy Gap (kcal/mol)

n-acene

KS-LDA KS-PBE KS-BLYP KS-BLYP-D TAO-LDA TAO-PBE TAO-BLYP TAO-BLYP-D

FIG. 11. Same as Fig.10, but for the larger acenes.

and

EAv= EN− EN+1, (14) respectively, with ENbeing the total energy of the N-electron system. Accordingly, the fundamental gap Eg can be calculated by

Eg= IPv− EAv= EN+1+ EN−1− 2EN. (15) As the size of the acene increases, IPv (see Fig.12) monotonically decreases and EAv (see Fig.13) monotonically in- creases, yielding a monotonically deceasing Eg(see Fig.14).

In contrast to the calculated IPv and EAv, Eg seems to be rather insensitive to the choice of the XC functionals in TAO- DFT. For 100-acene, Egis 0.54 eV for TAO-PBE, TAO-BLYP, and TAO-BLYP-D, and 0.55 eV for TAO-LDA. Note that the calculated Egis within the most interesting range (1 to 3 eV) for n-acene (n: 44 to 10), giving promise for applications of acenes in nanoelectronics.⁸⁹

The orbital occupation numbers in TAO-DFT provide information useful in assessing the possible polyradical char- acter of n-acene. As they are closely related to the NOONs,⁷⁶

3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5

10 20 30 40 50 60 70 80 90 100

Vertical Ionization Potential (eV)

n-acene

Experiment CCSD(T)/CBS TAO-LDA TAO-PBE TAO-BLYP TAO-BLYP-D

FIG. 12. Vertical ionization potential for the lowest singlet state of n-acene as a function of the acene length, calculated by various functionals in spin- unrestricted TAO-DFT. The experimental data are taken from the compilation in Ref.97, and the CCSD(T)/CBS data are taken from Ref.95.

(9)

-1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4

10 20 30 40 50 60 70 80 90 100

Vertical Electron Affinity (eV)

n-acene

FIG. 13. Vertical electron affinity for the lowest singlet state of n-acene as a function of the acene length, calculated by various functionals in spin- unrestricted TAO-DFT. The experimental data are taken from the compilation in Ref.97, and the CCSD(T)/CBS data are taken from Ref.99.

we calculate the symmetrized von Neumann entropy (e.g., see Eq. (9) of Ref.103)

S_vN= −1 2

σ=α,β

∞ i=1

{fi,σ ln(fi,σ)+ (1 − fi,σ) ln(1− fi,σ)}, (16) for the lowest singlet state of n-acene as a function of the acene length, by spin-restricted TAO-DFT.⁸⁹Here SvNessen- tially provides no contributions for a single-reference system ({fi, σ} are close to either 0 or 1), and quickly increases with the number of active orbitals ({f_{i, σ}} are fractional for active orbitals, and are close to either 0 or 1 for others). Note that S_vN, which is simply (−_2θ¹E_S^θ[{fi,α}, {fi,β}]), can be readily obtained in TAO-DFT. As shown in Fig. 15, S_vN, which is closely related to the polyradical character of n-acene, in- creases monotonically with the chain length, supporting the previous finding that large acenes should exhibit polyradical character.^96,98

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9

10 20 30 40 50 60 70 80 90 100

Fundamental Gap (eV)

n-acene

FIG. 14. Fundamental gap for the lowest singlet state of n-acene as a func- tion of the acene length, calculated by various functionals in spin-unrestricted TAO-DFT. The experimental data are taken from the compilation in Ref.97, and the CCSD(T)/CBS data are taken from Refs.95and99.

0 2 4 6 8 10 12 14 16 18 20

10 20 30 40 50 60 70 80 90 100

Symmetrized Von Neumann Entropy

n-acene

TAO-LDA TAO-PBE TAO-BLYP TAO-BLYP-D

FIG. 15. Symmetrized von Neumann entropy for the lowest singlet state of n-acene as a function of the acene length, calculated by various functionals in spin-restricted TAO-DFT.

V. CONCLUSIONS

We have proposed generalized-gradient approximations to TAO-DFT. The resulting TAO-GGAs have been shown to significantly outperform our previous TAO-LDA for a wide range of applications, such as thermochemistry, kinetics, and reaction energies. For noncovalent interactions, TAO-GGAs with empirical dispersion corrections have been shown to yield excellent performance. Due to their computational efficiency, TAO-LDA and TAO-GGAs have been applied to study the electronic properties of acenes, including the ST gaps, vertical ionization potentials, vertical electron affinities, fundamental gaps, and symmetrized von Neumann entropy (i.e., a measure of polyradical character). The ground states of acenes have been shown to be singlets for all the chain lengths studied here. With the increase of acene length, the ST gaps, vertical ionization potentials, and fundamental gaps decrease monotonically, while the vertical electron affinities and polyradical character increase monotonically.

Although only three GGAs (PBE, BLYP, and BLYP-D) are examined in this work, the good properties of various GGAs (e.g., those recently developed for broad applicabil- ity or for specific properties^104–110) in KS-DFT are expected to be preserved in TAO-DFT (with a sufficiently small θ

= 7 mhartree). However, for some multi-reference systems (e.g., H₂and N₂dissociation curves), TAO-DFAs (with θ= 7 mhartree) may not provide a sufficient amount of static corre- lation energy. As a system-dependent θ (related to the distri- butions of NOONs) is expected to enhance the performance of TAO-DFAs for a wide range of single- and multi-reference systems, we plan to pursue this in the future.

ACKNOWLEDGMENTS

This work was supported by the National Science Council of Taiwan (Grant No. NSC101-2112-M-002-017- MY3), National Taiwan University (Grant No. NTU-CDP- 102R7855), the Center for Quantum Science and Engineering at NTU (Subproject Nos.: 102R891401 and 102R891403), and the National Center for Theoretical Sciences of Taiwan.