Density functional theory with fractional orbital occupations Jeng-Da Chai

(1)

Jeng-Da Chai

Citation: J. Chem. Phys. 136, 154104 (2012); doi: 10.1063/1.3703894 View online: http://dx.doi.org/10.1063/1.3703894

View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v136/i15 Published by the American Institute of Physics.

Additional information on J. Chem. Phys.

Journal Homepage: http://jcp.aip.org/

Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors

(2)

Density functional theory with fractional orbital occupations

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 23 January 2012; accepted 29 March 2012; published online 17 April 2012)

In contrast to the original Kohn-Sham (KS) formalism, we propose a density functional theory (DFT) with fractional orbital occupations for the study of ground states of many-electron systems, wherein strong static correlation is shown to be described. Even at the simplest level represented by the local density approximation (LDA), our resulting DFT-LDA is shown to improve upon KS-LDA for multi-reference systems, such as dissociation of H2and N2, and twisted ethylene, while performing similar to KS-LDA for single-reference systems, such as reaction energies and equilibrium geometries. Because of its computational efficiency (similar to KS-LDA), this DFT-LDA is applied to the study of the singlet-triplet energy gaps (ST gaps) of acenes, which are “challenging problems” for conventional electronic structure methods due to the presence of strong static correlation effects. Our calculated ST gaps are in good agreement with the existing experimental and high-level ab initio data.

The ST gaps are shown to decrease monotonically with the increase of chain length, and become vanishingly small (within 0.1 kcal/mol) in the limit of an infinitely large polyacene. In addition, based on our calculated active orbital occupation numbers, the ground states for large acenes are shown to be polyradical singlets. © 2012 American Institute of Physics. [http://dx.doi.org/10.1063/1.3703894]

I. INTRODUCTION

As the problem of solving the N-electron Schrödinger equation quickly becomes intractable as the size of a system increases, the development of efficient and accurate electronic structure methods for large systems, continues being the subject of intense current interest. Over the past two decades, Kohn-Sham density functional theory (KS-DFT) (Refs.1and 2) has become one of the most popular theoretical approaches for calculations of electronic properties of large ground-state systems (up to a few thousand electrons), due to its favorable balance between cost and performance.^3–7 Recently, its time-dependent extension, time-dependent density functional theory (TDDFT), has also been actively developed for treat- ing electron dynamics and excited states of large systems with considerable success.^8,9

Although the exact exchange-correlation (XC) functional Exc[ρ] in KS-DFT has not been known, functionals based on the standard density functional approximations (DFAs), such as the local density approximation (LDA) and generalized gradient approximations (GGAs), can accurately describe short-range XC effects (due to the accurate treatment of on-top hole density), and are computationally favorable for large systems.^3–7 Although KS-DFAs have been successful in many applications, due to the lack of accurate treatment of nonlocality of XC hole, KS-DFAs can exhibit the following three types of qualitative failures: (i) self-interaction error (SIE), (ii) noncovalent interaction error (NCIE), and (iii) static correlation error (SCE). In situations where these failures occur, KS-DFAs can produce erroneous results.¹⁰There- fore, resolving the qualitative failures of KS-DFAs at a rea-

a)Electronic mail: [email protected].

sonable computational cost seems to be the first step toward finding an efficient and accurate electronic structure method for large systems.

The SIEs of KS-DFAs lead to drastic failures for problems such as barrier heights of chemical reactions, band gaps of solids, and dissociation of symmetric radical cations.¹¹ In TDDFT, SIE causes failures for problems such as Rydberg excitations in molecules and long-range charge transfer excitations between two well-separated molecules.¹²The SIEs of KS-DFAs can be greatly reduced by hybrid DFT methods,¹³ incorporating some of the exact Hartree-Fock (HF) exchange into the KS-DFAs. Over the past 20 years, several hybrid functionals, such as global hybrid functionals^14,15 and long- range corrected (LC) hybrid functionals,^16–20 have been de- veloped to improve the accuracy of Exc[ρ], extending the ap- plicability of KS-DFT to a wide range of systems.

The proper treatment of noncovalent interactions requires the accurate description of dynamical correlation effects at medium and long ranges, which cannot be properly captured by KS-DFAs.²¹ In particular, for dispersion-dominated (van der Waals (vdW)) interactions, KS-LDA tends to overesti- mate the binding energies, while KS-GGAs tend to give insufficient binding or even unbound results. The NCIEs of KS-DFAs can be efficiently reduced by the DFT-D (KS-DFT with empirical dispersion corrections) schemes,^22–25 which have shown generally satisfactory performance on a large set of noncovalent systems.^26,27 The dispersion corrections can also be computed less empirically by the exchange-hole dipole moment method²⁸ or by the local response dispersion method.²⁹ Alternatively, a fully nonlocal density functional for vdW interactions (vdW-DF) (Ref. 30) can also be adopted to reduce the NCIEs of KS-DFAs. Currently, perhaps the most successful approach to taking into account nonlocal

(3)

dynamical correlation effects is provided by the double hybrid (DH) methods,^31–33 which mix some of the HF exchange and some of the nonlocal orbital correlation energy from the second-order Møller-Plesset perturbation theory³⁴ into the KS-DFAs. DH functionals have shown an over- all satisfactory accuracy for thermochemistry, kinetics, and noncovalent interactions. In addition, the sharp increase in HF exchange in typical DH functionals has also greatly reduced the SIEs relative to KS-DFAs and conventional hybrid functionals.

Systems with strong static (nondynamical) correlation effects, such as bond-breaking reactions, diradicals, conjugated polymers, magnetic materials, and transition metal com- pounds, belong to the class of strongly correlated (SC) systems (multi-reference systems). Such a system is usually char- acterized by a small (or vanishing) energy gap between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO), the HOMO-LUMO gap. Despite their computational efficiency, the accurate treatment of SC systems poses a great challenge to KS-DFAs and hybrid DFT methods.^10,35The DH methods may also be inad- equate for SC systems, as the second-order perturbation energy components diverge to minus infinity for systems with vanishing HOMO-LUMO gaps. Within the framework of KS- DFT, fully nonlocal XC functionals, such as those based on random phase approximation (RPA), are believed to be essential for the accurate treatment of SC systems.^6,7However, compared with KS-DFAs, hybrid functionals, and DH functionals, RPA-type functionals are computationally very de- manding for large systems, and their applications to SC systems are too scarce to make a firm judgment on their accuracy.

Recently, the SIEs of RPA-type functionals have been shown to be severe, even for a simple one-electron system such as H⁺₂.³⁶

Aiming to reduce the SCEs of KS-DFAs without extra computational cost, we propose a DFT with fractional orbital occupations, rather than a fully nonlocal XC functional in KS- DFT. The rest of this paper is organized as follows. In Sec.II, we briefly describe the rationale for fractional orbital occupations. In Sec.III, we describe the formulation of this DFT and explain how strong static correlation is described by this DFT. At the simple LDA level, the performance of our resulting DFT-LDA is compared with KS-LDA for several single- and multi-reference systems in Sec.IV. Based on physical arguments and numerical investigations, the optimal parameter for this DFT-LDA is defined in Sec. V. Our conclusions are given in Sec.VI.

II. RATIONALE FOR FRACTIONAL ORBITAL OCCUPATIONS

For the exact DFT, the exact ground-state energy can be obtained, only when the exact ground-state density is available to insert into the exact ground-state energy functional. Therefore, the development of a generally accurate DFT method should involve not only an accurate approximation for the ground-state energy functional but also an appropriate representation (possibly in terms of orbitals and occupation numbers) of the ground-state density. However, much

less attention has been paid to the latter (representation of ground-state density) than to the former (ground-state energy functional). Due to the search over a restricted domain of densities, the exact ground-state density of interest may not be obtained within the framework of KS-DFT, in which case even the exact KS-DFT will fail (i.e., the exact XC functional may not be differentiable at the exact ground-state density).⁷ Noticeably, some of these situations occur for systems with vanishingly small HOMO-LUMO gaps (SC systems), indi- cating the close relationship between strong static correlation effects and representations of the ground-state density. There- fore, to accurately describe SC systems, it seems intuitive to focus on devising an appropriate representation for the exact ground-state density and a DFT associated with such a representation.

A ground-state density ρ(r) is said to be interacting v- representable, if it can be obtained from a ground-state wave function of an interacting N-electron Hamiltonian for some external potential v(r). The exact ρ(r) can be obtained by the full configuration interaction (FCI) method at the com- plete basis set limit (i.e., interacting v-representable),³⁷ and can be compactly expressed in terms of the natural or- bitals (NOs) {χi(r)} and natural orbital occupation numbers (NOONs) {ni},³⁸

ρ_FCI(r)=

∞ i=1

ni|χi(r)|², (1) where {ni}, obeying the following two conditions:

∞ i=1

ni= N, 0 ≤ ni≤ 1, (2) are related to the variationally determined coefficients of the FCI expansion. As shown in Eq. (1), the exact ρ(r) can be represented by orbitals and occupation numbers, highlighting the importance of ensemble-representable densities.

By contrast, in KS-DFT,² ρ(r) is assumed to be nonin- teracting pure-state (NI-PS) vs-representable, as it belongs to a one-determinantal ground-state wave function of a nonin- teracting N-electron Hamiltonian (KS Hamiltonian) for some local potential vs(r) (KS potential).^39–41Correspondingly, the KS orbital occupation numbers should be either 0 or 1. As the Aufbau principle (i.e., filling the KS orbitals in order of in- creasing energy) should be obeyed, ρ(r) can be expressed as the sum of the densities of the lowest N occupied KS orbitals {φi(r)},³

ρ_KS-DFT(r)=

N i=1

|φi(r)|². (3)

As ground-state densities of most nondegenerate atomic and molecular systems (e.g., closed-shell systems with sizable HOMO-LUMO gaps) are likely to be NI-PS vs-representable, the commonly used XC functionals in KS-DFT are reliably accurate for these systems (assuming that the SIEs and NCIEs of these functionals are not severe). However, if a one-determinantal ground-state wave function is insufficient to represent ρ(r), these XC functionals may not be reliably accurate, as they are all developed based on general properties

(4)

of systems where the KS wave function is a one-determinantal wave function.⁷

As shown by several researchers,^40–44 there are some reasonable ground-state densities that are not NI-PS vs- representable. Clearly, such densities cannot be obtained within the framework of KS-DFT. To remedy this situa- tion, KS-DFT has been extended to ensemble DFT (E-DFT), wherein ρ(r) is assumed to be noninteracting ensemble (NI- E) vs-representable, as it belongs to an ensemble of pure determinantal states of the noninteracting KS system.^45,46Cor- respondingly, ρ(r) can be expressed as

ρ_E-DFT(r)=

∞ i=1

g_i|φi(r)|², (4)

where the occupation numbers {gi}, obeying the following two conditions:

∞ i=1

g_i= N, 0 ≤ gi ≤ 1, (5)

are given by

g_i =

⎧⎪

⎨

⎪⎩

1, for i < _F x_i,for i = F

0, for i > F.

(6)

Here, i is the orbital energy of the i-th KS orbital φi(r), F

is the Fermi energy, and xi is a fractional number (between 0 and 1). As can be seen, fractional orbital occupations can occur only for the orbitals at the Fermi level.

In 1998, the close relationship between strong static cor- relation effects and non-NI-PS vs-representable densities was observed by Baerends and co-workers,⁴²who studied the¹_g⁺ ground state of the C2molecule (a system where the ground- state wave function is nondegenerate but has a strong multi- reference character), and investigated the possibility that the ground-state density ρ(r) may be NI-E vs-representable (as assumed in E-DFT), rather than NI-PS vs-representable (as assumed in KS-DFT). In their study, ρ(r) was first obtained from the highly accurate ab initio CI wave function method, wherein ρ(r) can be represented by the NOs {χi(r)} and NOONs {ni} in Eq. (1). An iterative method for the con- struction of the KS orbitals and the KS potential from the CI density was then adopted,⁴⁷combined with the constraint of integer occupations of the KS orbitals. The ρ(r) of C2

was found to be NI-PS vs-representable at a bond distance shorter than the equilibrium distance, while being non-NI- PS vs-representable at the longer bond distances, leading to non-Aufbau solutions with unoccupied KS orbitals having energies lower than those of the highest occupied KS orbitals (i.e., holes below the Fermi level). On the other hand, when the ρ(r) of C2 was represented by the ensemble solution during the above iterative procedure, no holes below the Fermi level were found, and the corresponding KS orbitals were close to the NOs. Based on their results, Baerends and co-workers⁴²argued that the KS orbitals (generated from Aufbau solutions) for a NI-PS vs-representable ground-state density are comparable to the NOs, while the KS orbitals (generated from non-Aufbau solutions) for a non-NI-PS vs-

representable ground-state density can be distinctly different from the NOs in order to incorporate the effect of the configuration mixing on the ground-state density. They concluded that the ground-state density of a system with strong static correlation effects is likely to be non-NI-PS vs-representable, in which case an ensemble representation (via fractional orbital occupations) of the ground-state density is crucial. Ar- guments in support of this are also available in Refs. 43 and48.

The idea of simulating strong static correlation effects by fractional occupation numbers (FONs) in DFT has spurred the development of the DFT-FON method,^48–52the spin-restricted ensemble-referenced KS method,^53,54and the fractional-spin DFT method,^10,35 with great success for some SC systems.

The practical implementation of E-DFT and related methods has, however, been impeded by a few factors as follows. First, a double-counting of correlation effects may occur, when the conventional approximate XC functional is evaluated with an ensemble density in Eq. (4), rather than a NI-PS vs-representable density in Eq. (3). By taking into account the double-counting effects, the XC functional in E- DFT may need to be re-derived. Second, the computational cost of E-DFT is more expensive than that of KS-DFT, which makes E-DFT less practical for the study of large ground-state systems. Third, the computational cost of analytical nuclear gradients (if available) for E-DFT is more expensive than that for KS-DFT, which makes it a formidable computational task to perform geometry optimization of large molecules for E-DFT.

In view of the above difficulties, we focus on the representation of the ground-state density from the exact theory in Eq. (1). Although the exact orbital occupation numbers for interacting electrons are intractable for large systems due to the exponential complexity, the approximate ones can, however, be properly simulated. Based on the statistical properties of strongly correlated eigenstates, Flambaum et al.

argued that the distribution of occupation numbers (the mi- crocanonical averaging of the occupation numbers) for a finite number of interacting Fermi particles (with two-body interaction) practically does not depend on a particular many- body system and has a universal form that can be ap- proximately described by the Fermi-Dirac distribution with renormalized parameters (i.e., orbital energies, chemical potential, and temperature).^55–57 Statistical effects of the interaction have been shown to be absorbed by introduction of the effective temperature.

In view of the close relationship between strong static correlation effects and representations of the ground-state density as well as the close relationship between the distribution of occupation numbers for interacting electrons (in the sense of statistical average) and the Fermi-Dirac distribution for noninteracting electrons, in this work, the ground-state density ρ(r) of a system of N interacting electrons (at zero temperature) in the presence of an external potential vext(r), is assumed to be noninteracting thermal ensemble (NI-TE) vs- representable, as it is represented by the thermal equilibrium density of an auxiliary system of N noninteracting electrons at a fictitious temperature θ ≡ kBTel(where kBis the Boltzmann constant, Tel is the temperature measured in absolute

(5)

temperature, and θ is the temperature measured in energy units) in the presence of a local potential vs(r). Correspond- ingly, ρ(r) can be expressed as

ρ(r)=^∞

i=1

f_i|ψi(r)|², (7) where the occupation number fiis the Fermi-Dirac function

f_i = {1 + exp[(i− μ)/θ]}⁻¹, (8) which obeys the following two conditions:

∞ i=1

f_i= N, 0 ≤ fi≤ 1, (9)

iis the orbital energy of the i-th orbital ψi(r), and μ is the chemical potential chosen to conserve the number of electrons N.

In Sec.III, we demonstrate how a DFT associated with the NI-TE vs-representable ρ(r) in Eq. (7), can be formu- lated. In other words, for a given fictitious temperature θ , the remaining “renormalized parameters” ({i} and μ) and the {ψi(r)}, can be self-consistently determined to represent the ρ(r) in Eq. (7), which then determines the ground-state energy of the system. Strong static correlation is shown to be described by a term related to the θ and {fi} in this DFT.

To avoid any possible confusion with KS-DFT, E-DFT, and finite-temperature DFT (FT-DFT),⁵⁸we refer to this DFT as thermally-assisted-occupation DFT (TAO-DFT). We wish to develop TAO-DFT with the following characteristics.

r

It is developed for ground-state systems at zero temperature.

r

It represents the ground-state density with orbitals and occupation numbers.

r

It may be used together with existing XC functionals in KS-DFT.

r

It reduces to KS-DFT in the absence of strong static correlation effects.

r

It treats single- and multi-reference systems in a more balanced way than KS-DFT.

r

It has similar computational cost as KS-DFT (e.g., energy and analytical nuclear gradients).

III. TAO-DFT

A. Self-consistent equations

Consider a system of N interacting electrons moving in an external potential vext(r) at zero temperature. Based on the HK theorems,¹the ground-state energy E[ρ], a functional of the ground-state density ρ(r), can be written as

E[ρ]= F [ρ] +

ρ(r)vext(r)dr, (10) where the universal functional

F[ρ]= T [ρ] + Vee[ρ] (11) is the sum of the interacting kinetic energy T[ρ] and the electron-electron repulsion energy Vee[ρ].

In KS-DFT,^1,2F[ρ] is usefully partitioned as

F[ρ]= Ts[ρ]+ EH[ρ]+ (T [ρ] + Vee[ρ]− Ts[ρ]− EH[ρ])

= Ts[ρ]+ EH[ρ]+ Exc[ρ]. (12) Here, Ts[ρ] is the noninteracting kinetic energy at zero temperature,

E_H[ρ]≡e² 2

ρ(r)ρ(r)

|r − r| drdr (13) is the Hartree energy, and Exc[ρ]≡ (T [ρ] + Vee[ρ]

− Ts[ρ]− EH[ρ]) is the XC energy defined in KS-DFT. In KS-DFT, Ts[ρ], the big unknown in terms of the density, is exactly treated by the use of KS orbitals. However, as previously discussed, the basic ansatz of KS-DFT (i.e., NI-PS v_s-representability of the given ρ(r)) can be violated for SC systems, in which case even the exact KS-DFT will fail to convey reliably accurate results.⁷

To make progress, a different representation of the ground-state density is adopted in TAO-DFT, wherein ρ(r) is represented by the thermal equilibrium density of an auxiliary system of N noninteracting electrons at a fictitious tempera- ture θ in the presence of some local potential vs(r). Aiming to achieve this representation, in contrast to the original KS partition, F[ρ] is partitioned into the following set of terms:

F[ρ]= A^θs[ρ]+EH[ρ]+ (T [ρ] + Vee[ρ]−A^θs[ρ]−EH[ρ])

= A^θs[ρ]+ EH[ρ]+ (T [ρ] + Vee[ρ]

− Ts[ρ]− EH[ρ])+ (Ts[ρ]− A^θs[ρ])

= A^θs[ρ]+ EH[ρ]+ Exc[ρ]+ Eθ[ρ]. (14) Here, Ts[ρ], EH[ρ], and Exc[ρ] are the same as those defined in KS-DFT, A^θ_s[ρ] is the noninteracting kinetic free energy at temperature θ , and Eθ[ρ]≡ Ts[ρ]− A^θs[ρ]= A^θs⁼⁰[ρ]

− A^θs[ρ] is the difference between the noninteracting kinetic free energy at zero temperature and that at temperature θ .

Substituting Eq. (14) into Eq.(10) and minimizing the E[ρ] with respect to ρ(r) (subject to the constraint that the number of electrons be N), yields the following Euler equation for the ground-state density ρ(r),

μ= δA^θ_s[ρ]

δρ(r) + vext(r)+ e²

ρ(r)

|r − r|dr +δE_xc[ρ]

δρ(r) +δE_θ[ρ]

δρ(r) , (15)

where μ is the chemical potential of the system.

To bypass the exact functional form of A^θ_s[ρ] (the big unknown in terms of the density) needed in Eq.(15), consider an auxiliary system of N noninteracting electrons moving in a local potential vs(r) at temperature θ . Based on Mermin’s theorems,⁵⁸the grand-canonical potential ^θ_s of this reference system, a functional of the thermal equilibrium density ρs(r), can be written as

^θ_s[ρs]= A^θs[ρs]+

ρs(r)[vs(r)− μs]dr, (16) where μsis the chemical potential of the reference system.

(6)

Minimization of the ^θ_s[ρs] with respect to ρs(r), gives the following Euler equation for the thermal equilibrium den- sity ρs(r):

μs = δA^θ_s[ρs]

δρ_s(r) + vs(r). (17) Comparing Eq.(17)with Eq.(15), shows that both minimiza- tions have the same solution ρs(r)= ρ(r), if we choose vs(r) (up to a constant) as

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

ρ(r)

|r − r|dr+δE_xc[ρ]

δρ(r) +δE_θ[ρ]

δρ(r) . (18) Alternatively, as A^θ_s[ρ] can be expressed exactly in terms of orbitals and occupation numbers (see below), Eq.(17)can also be handled, in an exact manner, by solving the one- electron Schrödinger equations for the potential vs(r), given

by

− ¯² 2me

∇² + vs(r)

ψi(r)= iψi(r), (19) and construct

ρ(r)= ρs(r)=

∞ i=1

fi|ψi(r)|², (20) where the occupation number fiis the Fermi-Dirac function

fi = {1 + exp[(i− μ)/θ]}⁻¹, (21) and the chemical potential μ is chosen to conserve the number of electrons N,

∞ i=1

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

The formulation of TAO-DFT leads to a set of self-consistent equations, Eqs.(18)–(22).

To obtain a self-consistent ground-state density in TAO- DFT: (i) Choose a trial ρ(r) to construct vs(r) by Eq.(18); (ii) solve Eq.(19), which gives {i, ψi(r)}; (iii) find μ by solving Eq. (22); (iv) determine {fi} by Eq. (21) and new ρ(r) by Eq.(20). This process is coupled with Eq.(18)to achieve self- consistency. When converged, the entropy functional reads

S_s^θ[{fi}] = −kB

∞ i=1

[filn(fi)+ (1 − fi) ln(1− fi)].

(23) The exact noninteracting kinetic free energy A^θ_s at the ficti- tious temperature θ , can be expressed in terms of {fi, ψi},

A^θ_s[{fi, ψ_i}] = Ts^θ[{fi, ψ_i}] − θ kB

S^θ_s[{fi}], (24) which is the sum of the kinetic energy

T_s^θ[{fi, ψi}] = − ¯² 2me

∞ i=1

fi

ψ_i^∗(r)∇²ψi(r)dr

=

∞ i=1

f_i_i−

ρ(r)vs(r)dr (25)

and entropy contribution

− θ

k_BS^θ_s[{fi}] = θ

∞ i=1

[filn(fi)+ (1 − fi) ln(1− fi)]

(26) of noninteracting electrons at temperature θ . Based on Eqs.(10)and(14), the total ground-state energy E[ρ] can be evaluated by

E[ρ]= A^θs[{fi, ψ_i}] +

ρ(r)vext(r)dr

+ EH[ρ]+ Exc[ρ]+ Eθ[ρ]. (27) To sum up, in TAO-DFT, the partition of F[ρ] in Eq. (14) and the exact treatment of A^θ_s[ρ] in Eq. (24), are shown to yield a set of self-consistent equations for the NI- TE vs-representable ρ(r) in Eq.(20). Note that these equations resemble the finite-temperature KS equations,^2,58so the implementation of TAO-DFT can be easily achieved using existing FT-DFT codes with a slight modification (i.e., replacing the XC free energy in FT-DFT to the sum of Exc[ρ] and Eθ[ρ]

in TAO-DFT). Hence, the computational cost of TAO-DFT is similar to that of KS-DFT or FT-DFT. Similar to FT-DFT,^2,58 due to the explicit inclusion of Fermi-Dirac occupation function in TAO-DFT, the entropy contribution (−_k^θ_BS_s^θ[{fi}]) in Eq. (26) is essential to make the total ground-state en- ergy functional E[ρ] variational⁵⁹ (e.g., making the nu- clear gradients of E[ρ] equal to the Hellmann-Feynman forces⁶⁰).

B. Spin-polarized formalism

For a system with Nαup-spin electrons and Nβdown-spin electrons, the standard computational approach is the spin- polarized (spin-unrestricted) formalism, wherein the funda- mental variables are the up-spin density ρα(r) and down-spin density ρβ(r) of the ground-state density

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

σ=α,β

ρ_σ(r). (28)

In analogy to the two-Fermi-level picture of spin-polarized KS-DFT,^48,61 spin-polarized TAO-DFT can also be formu- lated with the two-chemical-potential picture, wherein two noninteracting auxiliary systems at the same fictitious tem- perature θ are adopted: one described by the spin function α and the other by function β, with the corresponding ther- mal equilibrium density distributions ρs, α(r) and ρs, β(r) ex- actly equal to ρα(r) and ρβ(r), respectively, in the original interacting system at zero temperature. Similar to the previ- ous derivations (but using the spin-polarized extensions of the HK theorems^48,61 and the Mermin theorems^62,63 for the physical and auxiliary systems, respectively), one-electron Schrödinger equations for electrons with σ -spin (σ = α or β), can be obtained as follows (i runs for the orbital index):

− ¯² 2me

∇² + vs,σ(r)

ψi,σ(r)= i,σψi,σ(r), (29)

(7)

with a local potential

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

ρ(r)

|r − r|dr +δE_xc[ρα, ρ_β]

δρ_σ(r) +δE_θ[ρα, ρ_β]

δρ_σ(r) . (30) Here, Exc[ρα, ρβ] is the same as the XC energy defined in spin-polarized KS-DFT,^48,61 and Eθ[ρα, ρ_β]

≡ Ts[ρα, ρβ]− A^θs[ρα, ρβ]= A^θs⁼⁰[ρα, ρβ]− A^θs[ρα, ρβ] is the spin-polarized version of Eθ[ρ]. The σ -spin density can be constructed by

ρ_σ(r)=

∞ i=1

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

∞ i=1

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

The formulation of spin-polarized TAO-DFT yields two sets (one for each spin function) of self-consistent equations, Eqs. (29)–(33), for ρα(r) and ρβ(r), respectively, which are coupled with ρ(r) by Eq.(28).

To obtain self-consistent spin densities (and the ground- state density) in spin-polarized TAO-DFT: (i) Choose trial spin densities ρα(r) and ρβ(r) to compute the ground-state density ρ(r) by Eq.(28); (ii) for the σ -spin (σ= α or β) elec- trons, construct vs,σ(r) by Eq.(30); (iii) solve Eq.(29), which gives {i, σ, ψi, σ(r)}; (iv) find μσby solving Eq.(33); (v) de- termine {f_{i, σ}} by Eq.(32)and new ρσ(r) by Eq.(31). This process is coupled with Eq.(28)to achieve self-consistency.

When converged, A^θ_s,σ, the sum of the kinetic energy and en- tropy contribution of noninteracting σ -spin electrons at the fictitious temperature θ , is given by

A^θ_s,σ[{fi,σ, ψi,σ}]

= − ¯² 2me

∞ i=1

fi,σ

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

+ θ

∞ i=1

[fi,σ ln(fi,σ)+ (1 − fi,σ) ln(1− fi,σ)]

=^∞

i=1

{fi,σ_i,σ + θ [fi,σ ln(fi,σ)+ (1 − fi,σ)

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

ρσ(r)vs,σ(r)dr, (34) and the total ground-state energy E[ρα, ρβ] in spin-polarized TAO-DFT is evaluated by

E[ρα, ρβ]=

σ=α,β

A^θ_s,σ[{fi,σ, ψi,σ}] +

ρ(r)vext(r)dr + EH[ρ]+ Exc[ρα, ρ_β]+ Eθ[ρα, ρ_β]. (35)

Spin-unpolarized (spin-restricted) TAO-DFT can be for- mulated by imposing the constraints of ψi, α(r)= ψi, β(r) and fi, α= fi, βto spin-polarized (spin-unrestricted) TAO-DFT.

C. Analytical nuclear gradients

The analytical computation of nuclear gradients is crucial for the efficient optimization of molecular geometries. In light of the similarity of TAO-DFT and FT-DFT,^2,58analytical nuclear gradients for TAO-DFT can be easily obtained from those for FT-DFT with a slight modification (mentioned previously). Therefore, the computational cost of the analytical nuclear gradients for TAO-DFT is similar to that for KS-DFT or FT-DFT. For nonorthogonal atomic orbital representations (e.g., Gaussian-type orbitals), the generalized Pulay force for a noninteracting thermal ensemble⁶⁴ should be included to add the effect of the basis-set dependent response to the analytical nuclear gradients for TAO-DFT.

D. Local density approximation

In TAO-DFT, as the exact Exc[ρ] and Eθ[ρ], in terms of the ground-state density ρ(r), remain unknown, DFAs for both of them (denoted as TAO-DFAs) are needed for practical applications. The performance of TAO-DFAs depends on the accuracy of DFAs and the choice of the fictitious tem- perature θ . In this work, we adopt the LDA (the simplest DFA) for both the Exc[ρ] and Eθ[ρ] in TAO-DFT (denoted as TAO-LDA). As TAO-LDA is exact for a uniform electron gas, it provides a good starting point for more accurate and so- phisticated TAO-DFAs. Besides, TAO-LDA is readily avail- able, as E_xc^LDA[ρ] can be directly obtained from that of KS- LDA,^65,66 and E_θ^LDA[ρ] can be obtained with the knowledge of A^LDA,θ_s [ρ] as follows:

E^LDA_θ [ρ]≡ A^LDA,θs ⁼⁰[ρ]− A^LDA,θs [ρ]. (36) Here, Perrot’s parametrization of A^LDA,θ_s [ρ] (in its spin- unpolarized form) (Ref.67) is adopted to obtain E_θ^LDA[ρ] (in its spin-unpolarized form) by Eq.(36). For completeness of this work, A^LDA,θ_s [ρ] (after correcting some typos in Ref.67) is explicitly shown here (in atomic units),

A^LDA,θ_s [ρ]=

a_s^LDA,θ(r)dr, (37) where a^LDA,θ_s (r)≡ θρ(r)f (y) and y ≡ (π²/√

2)θ^−3/2ρ(r).

The function f(y) was parametrized separately for the two re- gions y≤ y0and y≥ y0(y0≡ ₄^3π^√₂),⁶⁷

f(y)= lny − 0.8791880215 + 0.1989718742y + 0.1068697043 × 10⁻²y²

− 0.8812685726 × 10⁻²y³+ 0.1272183027

× 10⁻¹y⁴− 0.9772758583 × 10⁻²y⁵ + 0.3820630477 × 10⁻²y⁶− 0.5971217041

× 10⁻³y⁷, for y≤ y0, (38)

f(y)= 0.7862224183u − 0.1882979454

× 10¹u⁻¹+ 0.5321952681u⁻³

(8)

+ 0.2304457955 × 10¹u⁻⁵− 0.1614280772

×10²u⁻⁷+ 0.5228431386 × 10²u⁻⁹

− 0.9592645619 × 10²u⁻¹¹ + 0.9462230172 × 10²u⁻¹³

− 0.3893753937 × 10²u⁻¹⁵,

for y≥ y0with u≡ y^2/3. (39) The θ = 0 case, A^LDA,θs ⁼⁰[ρ], is the same as the Thomas- Fermi kinetic energy density functional,^3,4

A^LDA,θ=0_s [ρ]= CF

ρ^5/3(r)dr, (40) where CF = ₁₀³(3π²)^2/3.

For spin-polarized (spin-unrestricted) TAO-LDA, the corresponding spin-polarized forms, E_xc^LDA[ρα, ρ_β] (available from that of spin-polarized KS-LDA; Refs. 65 and 66) and E_θ^LDA[ρα, ρ_β], should be adopted. From the spin- scaling relation of A^θ_s[ρα, ρ_β] (same as that of Ts[ρα, ρβ]⁶⁸), E_θ^LDA[ρα, ρ_β] can be conveniently expressed by its spin- unpolarized form E_θ^LDA[ρ],

E_θ^LDA[ρα, ρ_β]≡ A^LDA,θs ⁼⁰[ρα, ρ_β]− A^LDA,θs [ρα, ρ_β]

= 1

2(A^LDA,θ=0_s [2ρα]+ A^LDA,θ=0s [2ρβ])

−1

2(A^LDA,θ_s [2ρα]+ A^LDA,θs [2ρβ])

= 1

2{(A^LDA,θs ⁼⁰[2ρα]− A^LDA,θs [2ρα]) + (A^LDA,θs ⁼⁰[2ρβ]− A^LDA,θs [2ρβ])}

= 1

2{Eθ^LDA[2ρα]+ E^LDAθ [2ρβ]}. (41) E. Strong static correlation from TAO-LDA

The ground-state density ρ(r) of a strongly correlated system containing a sufficiently large number of electrons, can be represented by Eq. (1) (with the exact NOs {χi(r)}

and NOONs {ni}). Assume that the ρ(r) can also be represented by Eq.(20)(with the orbitals {ψi(r)} and their occu- pation numbers {fi} from the exact TAO-DFT). Note that for such a NI-TE vs-representable ρ(r), its “internal variables”

{fi} and {ψi(r)} in Eq.(20), can still be tuned by changing the fictitious temperature θ . If a θ is chosen so that {ni}≈ {fi} (in the sense of statistical average, as mentioned previously), we have {χi(r)}≈ {ψi(r)} (as both Eqs.(1)and(20)repre- sent the same ρ(r)). In fact, the NI-TE vs-representability of ρ(r) is likely to be fulfilled for this θ , due to the similarity of Eqs.(1)and(20).

Consequently, the exact kinetic energy of interacting electrons T[ρ] can be properly simulated by T_s^θ[{fi, ψi}]

(as appeared in A^θ_s[{fi, ψi}] = Ts^θ[{fi, ψi}] − _k^θ_BS_s^θ[{fi}]), namely,

T[ρ]= T [{ni, χi}] ≈ Ts^θ[{fi, ψi}], (42) due to their similar expressions,⁵ while the electron-electron repulsion energy Vee[ρ]= F [ρ] − T [ρ] (see Eqs. (11)

and(14)) is given by

V_ee[ρ]≈ EH[ρ]+ Exc[ρ]+ Eθ[ρ]− θ kB

S_s^θ[{fi}]. (43) On the right-hand side of Eq. (43), the first term is the Hartree energy, and the sum of the remaining terms (Exc[ρ]

+ Eθ[ρ]−_k^θ_BS_s^θ[{fi}]) should properly describe the XC energy defined in the exact wave function theory.

Here, we explain how strong static correlation is described by TAO-LDA, with arguments similar to the above.

Suppose that the exact ρ(r) in Eq.(1)can be reasonably represented by Eq. (20)(with the orbitals {ψi(r)} and their oc- cupation numbers {fi} from TAO-LDA). When applying the above arguments for TAO-LDA (i.e., choosing a θ so that {ni}

≈ {fi}, which gives {χi(r)}≈ {ψi(r)}), T[ρ] can still be prop- erly simulated by T_s^θ[{fi, ψi}] (see Eq.(42)), while Vee[ρ] is only approximated by

V_ee[ρ]≈ EH[ρ]+ Exc^LDA[ρ]+ E^LDAθ [ρ]− θ kB

S_s^θ[{fi}].

(44) On the right-hand side of Eq.(44), the first two terms (EH[ρ]

and E_xc^LDA[ρ]) are the same as those defined in KS-LDA, the third term E_θ^LDA[ρ] locally accounts for the difference be- tween the exact Ts and A^θ_s (at the LDA level), and the last term is the entropy contribution (−_k^θ_BS_s^θ[{fi}] ≈ −_k^θ_BS^θ_s[{ni}]

= θ_∞

i=1[niln(ni)+ (1 − ni) ln(1− ni)]) (see Eq.(26)).

Due to their local treatment, E_xc^LDA[ρ] and E_θ^LDA[ρ]

are not expected to properly describe nonlocal XC effects (e.g., long-range dynamical correlation and strong static correlation). However, as the entropy contribution is a fully nonlocal density functional ({fi} are implicit density functionals), it may describe nonlocal correlation effects. There is cer- tainly a close relationship between the entropy (defined by the NOONs {ni}) and correlation energy of a system. A famous example is given by the Collins conjecture⁶⁹ that the correlation energy of a system is proportional to the Jaynes (in- formation) entropy S_Jaynes[{ni}] = −_∞

i=1n_i ln(ni).⁷⁰ Inter- estingly, the entropy contribution in Eq.(44)is proportional to the Gibbs (thermodynamic) entropy (S^θ_s[{fi}] ≈ S^θs[{ni}]

= −kB

_∞

i=1[ni ln(ni)+ (1 − ni) ln(1− ni)]), with the constant of proportionality being explicitly given by (−_k^θ_B). Note that the similarity of information entropy and thermodynamic entropy in a many-body quantum system (with strong interactions) has been shown in Ref. 71, based on statistical arguments.

As the entropy contribution (−_k^θ_BS_s^θ[{fi}]) in Eq.(44)es- sentially provides no contributions for a single-reference sys- tem ({fi}≈ {ni} are close to either 0 or 1), and significantly lowers the total energy of a multi-reference system ({fi}

≈ {ni} are fractional for active orbitals, and are close to either 0 or 1 for others), we expect that this term (absent in KS-LDA) plays a crucial role in simulating strong static correlation (rather than dynamical correlation).

IV. NUMERICAL INVESTIGATIONS OF AN OPTIMALθ VALUE

The fictitious (reference) temperature θ for TAO-DFT, controlling the orbital occupation numbers {fi}, is closely

(9)

TABLE I. Statistical errors (in kcal/mol) of the reaction energies of 30 chemical reactions, Ref.20, calculated by TAO-LDA (with various θ (in mHartree)). The θ= 0 case corresponds to KS-LDA.

θ 0 1 3 5 7 10 15 20 30 40

MSE − 0.41 − 0.72 − 0.94 − 1.13 − 1.32 − 1.59 − 1.96 − 2.25 − 2.73 − 3.04

MAE 8.51 8.27 7.75 7.36 7.09 6.95 7.53 8.92 12.28 15.21

rms 11.10 10.89 10.31 9.76 9.38 9.16 9.75 11.25 15.20 19.03

Max(−) − 18.31 − 17.43 − 15.65 − 15.73 − 15.92 − 16.55 − 18.63 − 22.61 − 30.65 − 39.72

Max(+) 35.68 35.59 33.88 32.18 30.50 28.01 23.95 20.08 17.09 25.45

related to the strength of static correlation. At the LDA level, an immediate question is how the θ for the resulting TAO-LDA should be chosen. As previously argued, the entropy contribution, −_k^θ_BS_s^θ[{fi}] = θ_∞

i=1[fi ln(fi)+ (1 − f_i) ln(1− fi)], can be responsible for strong static correla- tion effects, especially when the {fi} (tunable by the θ ) prop- erly simulate the exact NOONs {ni}. To numerically inves- tigate this conjecture, the performance of TAO-LDA (with θ

= 0, 1, 3, 5, 7, 10, 15, 20, 30, and 40 mHartree) is exam- ined for both single-reference systems (reaction energies and equilibrium geometries) and multi-reference systems (dissociation of H2and N2, twisted ethylene, and singlet-triplet (ST) energy gaps of linear acenes). The limiting case where θ = 0 for TAO-LDA is especially interesting, as this reduces to KS- LDA. Therefore, it is important to know how well KS-LDA performs here to assess the significance of the extra parameter θfor TAO-LDA.

All calculations are performed with a development version of Q-Chem 3.2.⁷² The error for each entry is defined as (error = theoretical value − reference value). The nota- tion 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(+)). Results are computed using the 6-311++G(3df,3pd) basis set, unless noted otherwise.

A. Single-reference systems 1. Reaction energies

The accurate prediction of reaction energies is usually one of the major criteria in the assessment of the performance of electronic structure methods. The reaction energies of 30 chemical reactions (a test set described in Ref.20) are used to examine the performance of TAO-LDA. As shown in Table I, TAO-LDA (with a θ smaller than 10 mHartree)

has similar performance to KS-LDA.⁷³ This is unsurpris- ing, 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 numbers of TAO-LDA (with a sufficiently small θ ). Consequently, T_s^θ[{fi, ψ_i}] (see Eq.(42)) is close to T_s^θ⁼⁰[{ψi}] (KS kinetic energy), and E_θ^LDA[ρ] (see Eq.(36)) and the entropy contribution (−_k^θ_BS^θ_s[{fi}] = θ_∞

i=1[filn(fi)+ (1 − fi) ln(1− fi)]) have insignificant contributions to the total energy, relative to E_xc^LDA[ρ] (see Eq.(44)).

2. Equilibrium geometries

Geometry optimizations for TAO-LDA 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 near their equilibrium geometries can be well described by single-reference wave functions, TAO-LDA (with a θ smaller than 10 mHartree) is also found to perform simi- larly to KS-LDA,⁷³as shown in TableII.

B. Multi-reference systems 1. Dissociation of H₂and N₂

H₂dissociation, a single-bond breaking system, is partic- ularly challenging for KS-DFT. Figure1shows the potential energy curves (in total energy) for the ground state of H₂, calculated by both the spin-restricted and spin-unrestricted formalisms of the HF theory and KS-DFT (with LDA (Refs. 65 and66) and B3LYP (Refs. 13 and 14) functionals), where the exact potential energy curve is calculated by the coupled-cluster theory with iterative singles and dou- bles (CCSD).⁷⁵ Due to the symmetry constraint, the spin- restricted and spin-unrestricted potential energy curves, calculated by the exact theory, should be the same. Therefore,

TABLE II. Statistical errors (in Å) of EXTS,⁷⁴calculated by TAO-LDA (with various θ (in mHartree)). The θ

= 0 case corresponds to KS-LDA.

θ 0 1 3 5 7 10 15 20 30 40

MSE 0.004 0.004 0.004 0.005 0.005 0.005 0.006 0.008 0.015 0.030

MAE 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.014 0.021 0.036

rms 0.017 0.017 0.017 0.017 0.017 0.018 0.019 0.021 0.037 0.080

Max(−) − 0.091 − 0.091 − 0.091 − 0.091 − 0.091 − 0.092 − 0.095 − 0.101 − 0.110 − 0.110

Max(+) 0.081 0.081 0.080 0.080 0.080 0.080 0.078 0.083 0.222 0.581