A Paradigm for Density Functional Theory Using Electron Distribution on the Energy Coordinate
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug A Paradigm for Density Functional Theory UsingElectron Distribution on the Energy Coordinate
Hideaki Takahashi ∗ Department of Chemistry, Graduate School of Science,Tohoku University, Sendai, Miyagi 980-8578, Japan (Dated: August 10, 2020)Static correlation error(SCE) inevitably emerges when a dissociation of a covalent bond is de-scribed with a conventional denstiy-functional theory (DFT) for electrons. SCE gives rise to aserious overshoot in the potential energy at the dissociation limit even in the simplest molecules.The error is attributed to the basic framework of the approximate functional for the exchangecorrelation energy E xc which refers only to local properties at coordinate r , namely, the electrondensity n ( r ) and its derivatives. To solve the problem we developed a functional E exc which usesthe energy electron distribution n e ( ǫ ) as a fundamental variable in DFT. n e ( ǫ ) is obtained by theprojection of the density n ( r ) onto an energy coordinate ǫ defined with the external potential ofinterest. The functional was applied to the dissociations of single, double, and triple bonds in smallmolecules showing reasonable agreements with the results given by a high level molecular orbitalstheory. We also applied the functional to the computation of the energy change associated with spindepolarization and symmetrization in Carbon atom, which made an improvement over the conven-tional functional. This work opens the way for development of tougher functional that necessitatesnon-local properties of electrons such as kinetic energy functional. I. THE STATIC CORRELATION ERROR INDENSITY FUNCTIONAL THEORY
The density-functional theory (DFT) offers a ver-satile framework to study properties of interactingmany-particles systems such as liquids or electronsin molecules or bulk materials[1, 2]. DFT owes itssuccess in physics to the simple treatments of thecomplex correlations among the particles in terms oftheir distribution functions. For instance free energychange associated with an insertion of a solute into apure solvent can be formulated by a functional of thedistribution ρ of the solvent around the solute. In theelectronic structure calculations the electron density n plays an essential role in describing the exchange andcorrelation energies E xc of electrons. The key approachcommon to these theories is to project the effect of theparticles correlations onto the 1-body potential givenas a functional of the distribution ρ or of the electrondensity n . The construction of the functionals is, ofcourse, based on the statistical mechanics in the theoryof solutions and on the quantum mechanics in the DFTfor electrons.The DFT for electrons, Kohn-Sham(KS) DFT[3] inparticular, has been extensively utilized in the field ofphysics and chemistry. A vast amount of work has beendevoted to the foundation of the theoretical basis of theDFT and the developments of the functionals as well astheir applications. Actually, KS-DFT has been estab-lished as a reliable and efficient method in the electronicstructure calculation. However, the functionals forKS-DFT still suffer from two major problems[4]. First,in the approximate DFT functional, the self-interaction ∗ [email protected] of an electron due to Coulomb interaction does not fullycancel the exchange energy of itself, that is often referredto as self-interaction error (SIE). Second, static correla-tion error (SCE) gives rise to erroneous destabilizationin energy in describing the dissociation of a chemicalbond when spin-symmetry adapted wave functions areused. Indeed, these errors seriously affect the energeticseven in the simplest molecules. In the dissociation of theH molecule for instance, a GGA(generalized gradientapproximation)[2] level E xc functional overshoots thepotential energy at the dissociation by ∼
50 kcal/mol[5].The SCE arises more or less in stretching every covalentbond. The major objective of the present Letter is todevelop a simple and efficacious functional to solve theSCE problem on the basis of a new DFT paradigmwhere the electron distribution on the energy coordinate plays a role.The mechanism underlying the emergence of SCEcan be elucidated in terms of the exchange hole modelused in the local density approximation (LDA)[2] whichconstitutes the basis of the developments of all the E xc functionals. The exchange hole in the homogeneouselectron gas (HEG) is in particular the most prevailinghole model. Here, we consider the bond dissociationof an H molecule for simplicity. When the bond isinfinitely stretched with the spin symmetry maintained,the population of the electron with a spin σ on eachatom becomes half an electron. Hence, the hole depthmodeled by the LDA approach is also just half of theexact exchange hole since the depth of the hole is equalto the density n σ ( r ), which causes serious underestima-tion in the exchange energy. One way to compensatethe error is to incorporate some fraction of the exchangehole of the opposite spin into the total hole[5] althoughthere is no theoretical justification for the treatment. FIG. 1. Electron densities n σ ( σ = α, β ) of an H molecule atthe dissociation limit. The overlap of the electron densities ofindependent H atoms (upper) and the density of the entan-gled H atoms (lower). The shaded regions show the electrondistributions for spin α . The real numbers are fractional pop-ulations with α spin. The coupling parameter λ representsthe gradual variation of the density from n to n . II. FORMULATION OF A FUNCTIONAL WITHA DISTRIBUTION ON THE ENERGYCOORDINATE
To clarify the properties of the functional to be de-veloped in this Letter we present in Fig. 1 the gradualvariation of the electron density with spin σ from a refer-ence density n to n of the interacting system through acoupling parameter λ (0 ≤ λ ≤
1) for a hydrogen moleculewhen the bond distance R is infinitely extended( R = ∞ ). n is merely the overlap of the electron densities with op-posite spins for the independent H atoms placed at sitesH A and H B . n is the spin restricted density of the ‘en-tangled’ H atoms. In the following we drop the spinindex for the sake of notational simplicity unless other-wise stated. We also assume that the reference density n is given from the outset. It is apparent that the spatialbehavior of the electron density n ( r ) completely differsfrom that of n ( r ). The target exchange-correlation en-ergy E xc [ n ] can be expressed in the form of couplingparameter integration, thus, E xc [ n ] = E xc [ n ] + Z dλ d e E xc [ n λ ] dλ = E xc [ n ] + Z dλ Z d r dn λ ( r ) dλ d e E xc [ n λ ] dn λ ( r ) (1)where E xc is a certain conventional functional and e E xc represents a functional specific to the evaluation of thedifference in the exchange-correlation energy from λ = 0to 1. Of course, the integration in Eq. (1) should beevaluated as zero when the exact exchange-correlationfunctional E ex xc is employed( e E xc = E ex xc ). It is also impor-tant to note that any intermediate density expressed as n λ = (1 − λ ) n + λn should also give the same result, i.e. E xc [ n λ ] = E xc [ n ] = E xc [ n ] (0 < λ < e E xc .Here, we introduce the energy electron density n e ( ǫ )[6]by projecting the density n ( r ) n e ( ǫ ) = Z dǫ n ( r ) δ ( ǫ − υ ( r )) (2)where υ ( r ) is the external potential of interest used todefine the energy coordinate ǫ . Explicitly, υ ( r ) is givenby υ ( r ) = P A Z A / | r − R A | , where Z A and R A arecharge and position of nuclear A in the molecule. As wasproved in the previous work[6], there exists one-to-onecorrespondence between a subset of the external poten-tials and a subset of the energy electron distributions.Further, a rigorous framework of the density functionaltheory can be established almost in parallel to the con-ventional DFT for electrons[6]. The advantage of using FIG. 2. The spatial electron density n λ ( r ) of an H molecule(upper) and its projection n eλ ( ǫ ) onto the energy co-ordinate ǫ (lower left). At the dissociation limit n eλ ( ǫ ) staysconstant with respect to the variation of λ (lower right). the distribution n e ( ǫ ) is that n eλ ( ǫ ) stays constant for anarbitrary ǫ with respect to the variation of λ at R = ∞ .Hereafter, we attach the superscript e on the functionsand functionals that are relevant to the energy coordinateto avoid confusions. We replace the distribution n λ ( r ) inthe integrand by the corresponding energy electron den-sity n eλ ( ǫ ), which reads E xc [ n ] = E xc [ n ] + Z dλ Z dǫ dn eλ ( ǫ ) dλ d e E exc [ n eλ ] dn eλ ( ǫ ) (3)Assuming the linear dependence of n eλ ( ǫ ) on λ , we have E xc [ n ] = E xc [ n ]+ Z dǫ ( n e ( ǫ ) − n e ( ǫ )) Z dλ e υ exc [ n λ ]( ǫ ) (4)where e υ exc [ n eλ ]( ǫ ) is the exchange-correlation potentialdefined on the energy coordinate. Obviously, the inte-gral in Eq. (4) becomes zero at the dissociation limit forany choice of the approximate functional e E exc [ n e ] sincethe relation n e ( ǫ ) − n e ( ǫ ) = 0 is fully satisfied. Even ina weakly interacting system it is anticipated that n e isvery close to n e and, hence, the integral is evaluated tobe a finite but small value. The energy electron distribu-tion n e ( ǫ ) also provides a benefit in the evaluation of thecoupling parameter integration for e υ exc [ n eλ ] . The weakdependence of n eλ on λ justifies the application of the lin-ear response scheme in constructing e υ exc [ n e ]( ǫ ) from n e , υ xc [ n e ] ( ǫ ) = υ xc [ n e ] ( ǫ )+ Z dǫ ′ δ e υ exc [ n e ] ( ǫ ) δn e ( ǫ ′ ) ( n e ( ǫ ′ ) − n e ( ǫ ′ )) (5)Thus, the exchange-correlation potential υ xc [ n e ] can beconstructed only from the knowledge of the reference den-sity n e and the response function at n e . Equation (5)allows us to determine υ xc [ n e ] as well as n e in a self-consistent manner[6].Now, we formulate a functional E xc [ n ] provided thatthe exact spin-restricted density n is given. Without thecoupling parameter, E xc [ n ] in Eq. (3) can be merely ex-pressed as E xc [ n ] = E xc [ n ] + e E exc [ n e ] − e E exc [ n e ] (6)We define in Eq. (6) the static correlation E sc as E sc = E xc [ n ] − e E exc [ n e ] since e E exc [ n e ] itself describesthe E xc energy of the system with the density n e . Asdiscussed above Eq. (6) realizes E xc [ n ] = E xc [ n ] when R = ∞ since n e = n e is ensured. In other words, Eq.(6) works exactly when the corresponding exchange holeis equally fragmented onto the two atomic sites. Im-portantly, at the other limit where the two atomic sitescoincide( R = 0) E sc = 0 is attained under an appropriatechoice of e E exc . However, in a system with an intermedi-ately dissociated bond full inclusion of E sc gives rise toan artificial lowering of E xc as demonstrated in Ref. 6.A reasonable way to solve the problem is to attenu-ate the energy E sc in the intermediate region accordingto the hole delocalization. To this end we first rewriteEq. (6) equivalently in terms of the exchange-correlationenergy density U xc , E xc [ n ] = Z dǫ U exc [ n e ] ( ǫ ) n e ( ǫ )+ Z dǫ (cid:16) e U xc [ n ] ( ǫ ) − e U exc [ n e ] ( ǫ ) (cid:17) n e ( ǫ ) (7)where U exc ( ǫ ) is the exchange-correlation energy densityat ǫ and e U xc [ n ] ( ǫ ) is specifically the weighted averageof U xc [ n ]( r ) over the hypersurface with the energy co-ordinate ǫ and defined as e U xc [ n ] ( ǫ ) = 1 n e ( ǫ ) Z U xc [ n ] ( r ) n ( r ) δ ( ǫ − υ ( r )) d r (8)It is easy to see that R dǫ e U xc [ n ] ( ǫ ) n e ( ǫ ) = E xc [ n ].Then, we introduce a factor D e ( ǫ ) as a function of theenergy coordinate ǫ . D e ( ǫ ) measures the degree of the fragmentation of the exchange hole over multiple sites.Using the factor Eq. (7) is to be revised as E xc [ n ] = Z dǫ U exc [ n e ] ( ǫ ) n e ( ǫ )+ Z dǫF ( D e ( ǫ )) (cid:16) e U xc [ n ] ( ǫ ) − e U exc [ n e ] ( ǫ ) (cid:17) n e ( ǫ ) (9)where F is a function of D e ( ǫ ). Equation (9) shows thatthe static correlation E sc is incorporated into E xc ac-cording to the weight F ( D e ( ǫ )). F is supposed to bea monotonically increasing function ranging from 0 to 1and can be designed to realize E sc appropriately in theintermediate region. III. NUMERICAL IMPLEMENTATION OF THEFUNCTIONAL
We implement a numerical functional E xc with theform of Eq. (9) on the basis of the conventional DFTfunctional. To this end, we adapt BLYP[7, 8] knownas a representative GGA functional to evaluate e E exc interms of the energy electron density n e ( ǫ ). The vari-ables G [ n ( r )] needed as arguments for the BLYP func-tional, namely, n ( r ) itself and its derivatives ( |∇ n ( r ) | and ∇ n ( r )) are transformed to the corresponding quan-tities G e [ n e ( ǫ )] defined on the energy coordinate. Explic-itly, G e [ n e ( ǫ )] are obtained as the average of G [ n ( r )] overthe hypersurface with energy coordinate ǫ , thus, G e [ n e ( ǫ )] = Ω( ǫ ) − Z d r G [ n ( r )] δ ( ǫ − υ ( r )) (10)where Ω( ǫ ) is the volume of the region with coordinate ǫ and defined by Ω( ǫ ) = Z d r δ ( ǫ − υ ( r )) (11)It is readily recognized that each quantity G e [ n e ( ǫ )] hasthe same dimension with the original variable G [ n ( r )].Therefore, it is possible to substitute G e [ n e ( ǫ )] for G [ n ( r )] in adopting the conventional GGA. Of course,this is not a unique way of developing the functionalfor n e . However, it will be an instant but reasonableapproach to realize the functional based on the energycoordinate. It should be stressed that e E exc with the vari-ables on the energy coordinate (Eqs. (10) and (11)) pro-vides exactly the same results with the original GGAfunctional in principle when it is applied to hydrogenicatoms, i.e. spherically symmetric 1-electron systems, be-cause G [ n ( r )] are constant over a contour surface withan energy coordinate.Evaluation of the delocalization factor D e ( ǫ ) and thedetermination of the correlation factor F in Eq. (9)are the keys to the success of the application of thepresent theory. In this work we utilize the generalizedBecke-Roussel(GBR)[5] exchange hole for the evaluationof D e ( ǫ ). In the GBR approach the exchange hole ata reference point r is modeled by the hole for a hydro-genic atom. Explicitly, the hole function n GBR x ( A, ζ, r ) isexpressed by a Slater function, n GBR x ( A, ζ, r ) = A exp( − ζr ) (12)where r is the distance between the reference point andthe center of the exchange hole. ζ is a positive real num-ber and represents the width of the hole. In the originalBR approach[9] the factor A in Eq. (12) is determined asa function of ζ to ensure that the hole population is nor-malized to one. In contrast, in the GBR method the fac-tor A is left as a flexible parameter which can accommo-date the fragmentation of the exchange hole associatedwith a bond dissociation. Explicitly, A is determined sothat the Coulomb potential of the hole given by Eq. (12)realizes the exact exchange energy density U exact x ( r ) ofthe real system U exact x ( r ) = 12 Z d r ′ | r − r ′ | n x ( r , r ′ ) (13)where n x ( r , r ′ ) stands for the exchange hole associatedwith r . In addition, we have two additional conditionsto obtain the other parameters ( ζ, r ), for which we referthe readers to Ref. 5. It is known that the GBR param-eters ( A, ζ, r ) can be uniquely determined for any givenreference point through a certain numerical approach. Itis, then, possible to define the projected hole localized ata reference point r . The population P GBR ( r ) of the pro-jected hole is evaluated as 8 πA/ζ . Obviously, P GBR ( r )is relevant to the degree of the fragmentation of the realexchange hole. Actually, for an H molecule, it is antici-pated that P GBR ( r ) becomes 0 . R → ∞ . Hence, anatural choice of the delocalization factor D ( r ) in termsof the GBR exchange hole model is to take D ( r ) = N s (1 − P GBR ( r )) (14)where N s is the number of sites onto which the exchangehole is to be fragmented and N s = 2 is supposed forthe dissociation limit of a diatomic molecule. In moregeneral approach that can be applied to the multi-sitefragmentation it would be appropriate to take N s as theinverse of weighted average of P GBR ( r ), thus, N s = (cid:18) N Z d r P GBR ( r ) n ( r ) (cid:19) − (15)where N is the number of electrons with the spin σ . Thedelocalization factor D e ( ǫ ) dependent on the energy co-ordinate ǫ is merely obtained by the weighted averageover the isosurface with coordinate ǫ , D e ( ǫ ) = n e ( ǫ ) − Z d r D ( r ) n ( r ) δ ( ǫ − υ ( r )) (16)which can be justified when D ( r ) is almost constant overthe isosurfaces of the external potential υ . In Fig. 3 D ( r ) FIG. 3. The delocalization factor D ( r ) mapped onto the iso-surfaces of the external potential ǫ = υ ( r ) for H moleculewith various bond distances R . Depicted isosurfaces are for ǫ = 0 . , .
8, and 1 . ǫ on the molecularplane are also drawn with the gray thin lines ( ǫ = 0.4, 0.5,0.7, 1.0, 1.2, 2.0, and 4.0 au). is mapped with colors on the isosurfaces of the potential υ for the various bond distances R of H . It is shownin the figure that D ( r ) increases with increasing bonddistance R as expected. At each bond distance it is alsofound that D ( r ) tends to have relatively larger valuesnear the nuclei, that is, D ( r ) is large in the region wherethe coordinate ǫ has larger values. It is also interestingto note that D ( r ) is relatively smaller in the spatial re-gion between the nuclei. The projection of Eq. (16) canbe justified when D ( r ) is almost constant over the iso-surfaces of the external potential υ . When R is large, itseems that D ( r ) is nearly constant on an isosurface. Onthe other hand, for smaller R the constancy of D ( r ) isnot fulfilled on the isosurfaces with small υ in particu-lar. However, D ( r ) itself stays small in the region whereboth R and υ are small. Futhermore, the energy electrondensity n e ( ǫ ) is also small on the smaller coordinate ǫ .Therefore, the error due to Eq. (16) is supposed to besmall by consulting Eq. (9). In Fig. 3 number of sites N s given by Eq. (15) are also presented, which shows N s increases with increasing distance R as expected. Inprinciple, N s becomes 1 . . R = 0 and ∞ ,respectively.For the correlation function F ( D e ( ǫ )) in Eq. (9), wefirst tested the simplest model that F ( D e ( ǫ )) = D e ( ǫ ).However, we found it underestimates E sc in the interme-diate region of R as will be shown later. Thus, we intro-duce upward-convex polynomial functions F k ( D e ( ǫ )) F k ( D e ( ǫ )) = 1 − (1 − D e ( ǫ )) k ( k = 0 , , · · · ) (17)Specifically, F is identically zero and F k asymptoticallygoes to 1 when k → ∞ . Thus, the series of the poly-nomial functions F k covers the full range of the corre-lation strength. It is also important to note that therelations F k (0) = 0 and F k (1) = 1 is guaranteed forany choice of k ≥
1. In this work, the order k solelyserves as an adjustable parameter in the numerical im-plementation of Eq. (9) besides those intrinsic to theBLYP functional. Although k can be a real number, weonly consider a positive integer because a precise opti-mization of the functional is not the major objective ofthis work. 0 ≤ D e ( ǫ ) ≤ D e ( ǫ ) possibly becomeslarger than 1. Hence, F k ( x ) = 1 . < x ) has to be im-posed on F k ( x )( k = 1 , , · · · ).For the numerical calculations we utilized ‘Vmol’program[10–12] revised for the present development. Theprogram is based on the real-space grid approach[13, 14]where the interaction between valence electrons and nu-clei is represented with non-local pseudopotentials in theKleinman and Bylander separable form[15]. The chargesof the core electrons are summed to the nuclear chargeto evaluate the energy coordinate for the valence elec-trons. The energy distribution functions n e ( ǫ ) is con-structed for the densities on the discrete grid points.To increase the sampling points we introduced a dou-ble grid technique[16] where the density and its deriva-tives on a dense grid are obtained through the 4th-orderpolynomial interpolations of the values on the coarsegrids. The number of the double grids between the coarsegrids is set at N den = 5. As was done in the previ-ous paper[6], the axis for the energy coordinate ǫ wasuniformly descretized with N e points in the lower en-ergy region( ǫ min ≤ ǫ ≤ ǫ core ), while N e points were uni-formly placed in the log-scaled coordinate in the higherenergy region( ǫ core ≤ ǫ ≤ ǫ max ). The explicit values for( N e , N e ) and ( ǫ min , ǫ core , ǫ max ) as well as the method todiscretize the energy coordinate are presented in ‘Supple-mental Material’. The reference electron density n wasbuild with the LDA functional and the BLYP functionalwas utilized in the functional of Eq. (9). IV. RESULTS AND DISCUSSIONA. Bond Dissociations
We first applied Eq. (9) to the dissociation of an H molecule. Here, we leave the sole adjustable parameter k in the correlation factor(Eq. (17)) undetermined. Weperformed a series of calculations with various k valuesincluding pure BLYP functional based on the conven-tional DFT. As a reference we also carried out CCSD(T)calculation by conducting ‘Gaussian 09’ program pack-age, where aug-cc-pVQZ basis set was employed. Theresults are shown in Fig. 4 where the graphs for F k aresuperposed. It is shown that the BLYP energy with thestatic correlation (SC) using factor F does not fully re-produce the dissociation limit though it provides a sig-nificant improvement on the pure BLYP functional. Thelarger correlation factors F and F almost successfullyrealize the dissociation behavior though factor F rather overestimates SC in the intermediate H − H distance. Ac-tually, dissociation energies E d are, respectively, com-puted as 107.8 and 107.1 kcal/mol with these factors,showing good agreements with the CCSD(T) calculation(109.2 kcal/mol). Thus, we determined the parameter as k = 3 in the following applications.Next, we examined the dissociation of a double bond k = 1 k = 3 k = 5 k = 0 F k ( D ) D BLYPCCSD(T)BLYP+SC with F BLYP+SC with F BLYP+SC with F P o t e n t i a l E n e r gy / a . u . ! ! ! ! ! ! ! ! R (H-H) / a.u. FIG. 4. Potential energy curves for dissociating H molecule.The graphs for the correlation factors F k of Eq. (17) as func-tions of the delocalization factor D e ( ǫ )(Eq. (16)) are super-posed. in C H . To this end the bond length is graduallyincreased with the geometry of the CH groups beingfixed at that of the equilibrium structure in ethene. Instretching the bond we notice that an anti-bonding σ ∗ orbital comes in the set of the low-lying N val orbitalsin place of the π orbital, which causes numerical insta-bilities in the potential energy(PE). To avoid such anunfavorable situation the wave function obtained for aC − C distance is used as the initial guess for the SCFof the next ethene configuration created by a small in-crement of the bond length. The results are presentedin Fig. 5 where the delocalization factor D ( r ) mappedonto the isosurfaces of υ ( r ) are also shown in the inset.It is noticeable that the use of the correlation factor F gives the dissociation energy of 180 . . F fails to attain the dissociationlimit of F although the behavior of PE in the intermedi-ate region is more favorable. This is simply because D ( r )is evaluated to be rather small even at R (C − C) = 7.5 a.u.as shown in the inset. The population of the projectedhole by GBR approach is anticipated to be 0.5 on eachsite at the dissociation limit. Such a situation is almostrealized in H as shown in Fig. 3. However, in othersystems the hole does not perform ideally as suggestedin Refs. 17 and 18. That is, the hole population tendsto be overestimated, which gives rise to the underestima-tion of D ( r ) according to Eq. (14). As a consequencethe number of site N s of Eq. (15) will also be underes-timated. Actually, N s is estimated to be only 1 .
124 at R (C − C) = 9.5 a.u. Thus, the major role of the correla-tion factor F k is to compensate the drawback. In otherwords, another method to determine D ( r ) will possiblyimprove the functional.Another severe test was performed by studying the BLYPUCCSD(T)BLYP+SC with F BLYP+SC with F R (C-C) = 7.5 a.u. P o t e n t i a l E n e r gy / a . u . ! ! ! ! ! ! R (C-C) / a.u. FIG. 5. Potential energy curves for the dissociating doublebond in C H . The delocalization factor D ( r ) is depicted inthe inset. The color assignment for D ( r ) is common to thatfor Fig. 3. The contour surfaces of υ ( r ) are for 3.0, 4.0, 5.0,and 7.0 a.u. and the contour lines are for 2.0, 2.5 a.u. All thepotential energy curves are shifted so that the energies of theleft ends of the graphs match that of ‘BLYP+SC with F ’. dissociation of the triple bond in an N molecule. The re-sults are shown in Fig. 6. As a reference the result givenby UCCSD(T)/cc-pVQZ calculation is also presented.From the UCCSD(T) potential energy curve, dissociationenergy E d is estimated to be 222.4 kcal/mol. In contrast,the present calculations with factors F and F provide E d = 272 . . E d of the triple bond. It is also notable that our methodunderestimates PE of the intermediate N − N distances.As shown in the inset the values of D ( r ) are again esti-mated to be rather small although they are larger thanthose obtained for C H at R (C − C) = 7.5 a.u. This maycause the excess rise in the PE of the dissociating region.However, it is quite curious that even the factor F couldnot fully compensate the error as it performed nicely forC H . Anyway, it was demonstrated that the present ap-proach with F reduces E d of the spin-restricted BLYPby ∼
136 kcal/mol, showing a better performance thanB13[17] and KP[18] functionals.As a overall remark the functional of Eq. (9) is almostsuccessful in realizing the correct E d of the various chem-ical bonds. We found, however, the static correlationsfor the intermediately dissociated bonds tend to be over-estimated, which leads artificial lowering of the PE. Theerror might be attributed to the underestimation of D ( r )even in the molecules with sufficiently stretched bonds.Thus, a strong correlation factor F ( ǫ ) in Eq. (9) is neededto compensate the drawback, which unexpectedly lowersthe exchange energies in the intermediate bond distances.Sophistication of the method to evaluate D ( r ) would im- BLYPUCCSD(T)BLYP+SC with F BLYP+SC with F R (N-N) = 7.0 a.u. P o t e n t i a l E n e r gy / a . u . ! ! ! ! ! ! R (N-N) / a.u. FIG. 6. Potential energy curves for the dissociating triplebond in N . The delocalization factor D ( r ) is depicted in theinset. The color assignment for D ( r ) is common to that forFig. 3. The contour surfaces of υ ( r ) are for 3.0, 4.0, and 5.0a.u. and the contour lines are for 2.0, 2.5, and 8.0 a.u. All thepotential energy curves are shifted so that the energies of theleft ends of the graphs match that of ‘BLYP+SC with F ’. prove the overall behavior of PE. B. Spin-depolarization and Symmetrization
It is proved in Ref. 19 that the exact functional E exact for the total energy satisfies the following relation, E exact [ n mix ] = E υ ( N ) (18)where E υ ( N ) is the ground state energy of the N elec-trons subject to the external potential υ . n mix is thelinear combination of the g -degenerate ground state den-sities { n i } ( i = 1 , , · · · , g ) for υ , n mix = g X i =1 C i n i (19)with the normalization condition P gi C i = 1. It shouldbe noted that n mix is not υ -representable in general[20]and assessing the realization of Eq. (18) in an ap-proximate functional is known as a stringent test forstrong correlation[17]. To examine the performanceof the functional of Eq. (9) we construct n mix bymixing the three degenerate densities( g = 3) of groundstate Carbon atom. We suppose that each density n i is spin depolarized , that is, 2p-shell occupancies of(2 p x α ) . (2 p y α ) . (2 p z α ) (2 p x β ) . (2 p y β ) . (2 p z β ) areassumed in a component density. Then, a mixing param-eter λ (0 ≤ λ ≤
1) is introduced so that the three densitiesgive the equivalent contributions at λ = 1 /
2. Explicitly,we adopt C = λ q , C = (1 − λ ) q , C = 1 − C − C with q = ln 3ln 2 . Accordingly, the occupancy of eachp orbital is given as (1 − C λi ). It is readily recog-nized that in the densities n λ =0 , two of the three porbitals are individually occupied with an electron (0.5electron with spin α and 0.5 electron with spin β oneach orbital), while the three orbitals have equivalentoccupancies at λ = 1 /
2. In Fig. 7 we provide theenergy difference ∆ E ( λ ) between the total energy ofCarbon atom with the spin polarized density, namely(2 p x α ) (2 p y α ) (2 p z α ) (2 p x β ) (2 p y β ) (2 p z β ) and theenergy E [ n λ mix ] computed with Eq. (9). The exactfunctional realizes ∆ E ( λ ) = 0 irrespective of λ . Inapplying Eq. (9) we adopt the energy density n e ofthe supermolecule as the reference, which consists ofinfinitely separated Carbon atoms with the spin polar-ized densities. The result by the conventional BLYPfunctional is also presented to make comparisons. TheCoulomb interaction J [ n λ mix ] = R d r d r ′ n λ mix ( r ) n λ mix ( r ′ ) | r − r ′ | is also shown where J [ n λ =0mix ] = J [ n λ =1mix ] is taken as thestandard.The specific energy differences ∆ E (0) or ∆ E (1) showthe energy ∆ E depol due to the spin depolarization, while∆ E (1 /
2) subtracted by ∆ E depol is the energy ∆ E sym dueto symmetrization, i.e. ∆ E sym = ∆ E (1 / − ∆ E depol .It is notable in the figure that ∆ E depol is estimated tobe almost 0 by BLYP+SC(Eq. (9)) with the correlationfactor F although use of the factor F slightly degradesthe energy. It was demonstrated that the functionalworks nicely in computing the energy change ∆ E depol associated with the spin-depolarization. On the otherhand, ∆ E depol given by the conventional GGA func-tional(BLYP) becomes positively much larger ( ∼ J [ n λ mix ] decreaseswith the increase in λ from 0 to 1 / n λ mix fully delocalizes at λ = 1 / E ( λ ) of Eq. (9) withthe factor F unfavorably decreases along with J [ n λ mix ]suggesting E exc [ n λ mix ] itself stays almost constant. Theenergy ∆ E ( λ ) given by the factor F shows almost thesame tendency as F . Although the values ∆ E (1 / F as well as F are much closer to 0 than thatcomputed with the BLYP functional, ∆ E sym itself byEq. (9) is larger than that with the BLYP functional.The constancy of ∆ E ( λ ) with BLYP implies that thedecrease in J [ n λ mix ] is perfectly compensated with theincrease in the E xc [ n λ mix ]. The problem in ∆ E sym hasa relevance with the other error in DFT referred to asself-interaction error(SIE) and should be considered inthe forthcoming issues. V. CONCLUSION
We developed an electron density functional to com-pute strong electron correlation referred to as static cor-relation(SC) that emerges when a covalent bond is split.In contrast to the conventional DFT the energy electrondistribution n e ( ǫ ) instead of n ( r ) plays a fundamental Px Py Pz
BLYPBLYP+SC with F BLYP+SC with F Coulomb ! E ( ! ) / a . u . " " ! FIG. 7. Energy differences ∆ E between total energies of thespin restricted and unrestricted densities of Carbon atom areshown as functions of the mixing parameter λ . The 2p orbitalsin Carbon atom are depicted as the set of three boxes, inwhich α spins in the spin restricted density are representedwith upward bold arrows. The value on each box shows thespin population dependent on the mixing parameter. It issupposed that β spin (not shown) has the same population as α on each orbital. role in the functional development on the basis of the rig-orous framework of DFT as discussed in Ref. 6. The keyfeature of the new distribution is that n eλ ( ǫ ) stays con-stant with respect to the variation of parameter λ whichcouples the spin polarized density with that for the en-tangled state when the fragments are infinitely separatedin a homolytic dissociation.A remarkable progress was made in the functional byintroducing the correlation factor as a function of the de-gree D e ( ǫ ) of delocalization of the exchange hole on thecoordinate ǫ . In the present work, D e ( ǫ ) is obtained byutilizing the generalized Becke-Roussel(GBR) machinery.The functional was applied to the computation of thesplits of single, double and triple bonds. Namely, the dis-sociations of the bonds in H , C H , and N were stud-ied. It was demonstrated that the dissociations of thesingle and the double bonds are faithfully realized withthe correlation factor F showing reasonable agreementswith those obtained by sophisticated molecular orbitaltheories. However, the dissociation energy was found tobe overestimated in N although substantially improvedas compared with the conventional GGA functional. Thesource of the discrepancy might be attributed to the un-derestimation in D e ( ǫ ) by means of the GBR method.The present functional was also applied to the computa-tion of the energy change due to the spin depolarizationand symmetrization in Carbon atom, which showed a re-markable improvement over the conventional one.The present work suggests the possibility of a paradigmfor the electronic density functional based on the energycoordinate ǫ . In the conventional DFT functional theelectronic correlation of an electron at r is solely deter-mined by the local properties, namely, the electron den-sity and its derivatives at r . This condition will impose asevere limitation on the applicability of functional. Thepresent DFT approach in terms of the energy electrondensity offers a simple framework to overcome the prob-lem due to the spatial locality intrinsic to the LDA-basedDFT. The development of the kinetic energy functionalis undoubtedly the toughest challenge in DFT, where in-clusion of a non-local term is considered to be essentialas discussed in Ref. 21. The use of the energy coordinateinstead of r in DFT will serve a potential approach to akinetic energy functional. ACKNOWLEDGMENTS
This paper was supported by the Grant-in-Aid for Sci-entific Research on Innovative Areas (No. 23118701)from the Ministry of Education, Culture, Sports, Science,and Technology (MEXT); the Grant-in-Aid for Challeng-ing Exploratory Research (No. 25620004) from the JapanSociety for the Promotion of Science (JSPS); and theGrant-in-Aid for Scientific Research(C) (No. 17K05138)from the Japan Society for the Promotion of Science(JSPS). This research also used computational resourcesof the HPCI system provided by Kyoto, Nagoya, andOsaka university through the HPCI System ResearchProject (Project IDs: hp170046, hp180030, hp180032,hp190011, and hp200016). [1] J.-P. Hansen and I. R. Mcdonald,
Theory of Simple Liq-uids (Academic Press, 2013).[2] R. G. Parr and W. Yang,
Density-functional theory ofatoms and molecules (Oxford university press, New York,1989).[3] W. Kohn and L. J. Sham, Self-consistent equations in-cluding exchange and correlation effects, Phys. Rev. ,A1133 (1965).[4] A. J. Cohen, P. Mori-S´anchez, and W. Yang, Insightsinto Current Limitations of Density Functional Theory,Science , 792 (2008).[5] A. D. Becke, A real-space model of nondynamical corre-lation, J. Chem. Phys. , 2972 (2003).[6] H. Takahashi, Density-functional theory based on theelectron distribution on the energy coordinate, J.Phys. B: Atomic, Molecular and Optical Physics ,055102(11pp) (2018).[7] A. D. Becke, Density-functional exchange-energy approx-imation with correct asymptotic behavior, Phys. Rev. A , 3098 (1988).[8] C. Lee, W. Yang, and R. G. Parr, Development ofthe Colle-Salvetti correlation-energy formula into a func-tional of the electron density, Phys. Rev. B , 785(1988).[9] A. D. Becke and M. R. Roussel, Exchange holes in in-homogeneous system: a coordinate-space model, Phys.Rev. A , 3761 (1989).[10] H. Takahashi, T. Hori, T. Wakabayashi, and T. Nitta,A density functional study for hydrogen bond energy byemploying real space grids, Chem. Lett. , 222 (2000).[11] H. Takahashi, T. Hori, T. Wakabayashi, and T. Nitta,Real Space Ab Initio Molecular Dynamics Simulations forthe Reactions of OH Radical/OH Anion with Formalde- hyde, J. Phys. Chem. A , 4351 (2001).[12] H. Takahashi, T. Hori, H. Hashimoto, and T. Nitta, Ahybrid qm/mm method employing real space grids forqm water in the tip4p water solvents, J. Comp. Chem. , 1252 (2001).[13] J. R. Chelikowsky, N. Troullier, K. Wu, and Y. Saad,High-order finite-difference pseudopotential method: anapplication to diatomic molecules, Phys. Rev. B ,11355 (1994).[14] J. R. Chelikowsky, N. Troullier, and Y. Saad, Finite-difference pseudopotential method: electronic structurecalculations without a basis, Phys. Rev. Lett. , 1240(1994).[15] L. Kleinman and D. M. Bylander, Efficacious formfor model pseudopotentials, Phys. Rev. Lett. , 1425(1982).[16] T. Ono and K. Hirose, Timesaving double-grid methodfor real-space electronic-structure calculations, Phys.Rev. Lett. , 5016 (1999).[17] A. D. Becke, Density functionals for static, dynamical,and strong correlation, J. Chem. Phys. , 074109(10)(2013).[18] J. Kong and E. Proynov, Density functional model fornondynamic and strong correlation, J. Chem. TheoryComput. , 133 (2016).[19] W. Yang, Y. Zhang, and P. W. Ayers, Degenerate groundstates and a fractional number of electrons in densityand reduced density matrix functional theory, Phys. Rev.Lett. , 5172 (2000).[20] M. Levy, Electron densities in search of hamiltonians,Phys. Rev. A , 1200 (1982).[21] L.-W. Wang and M. P. Teter, Kinetic-energy functionalof the electron density, Phys. Rev. B45