Excitation energies through Becke's exciton model within a Cartesian-grid KS DFT
aa r X i v : . [ phy s i c s . c h e m - ph ] F e b Excitation energies through Becke’s exciton model within aCartesian-grid KS DFT
Abhisek Ghosal, Tarun Gupta, Kishalay Mahato and Amlan K. Roy ∗ Department of Chemical SciencesIndian Institute of Science Education and Research (IISER) KolkataNadia, Mohanpur-741246, WB, India
Abstract
Photon-induced electronic excitations are ubiquitously observed in organic chromophore. Inthis context, we present a simple, alternative time-independent DFT procedure, for computationof single-particle excitation energies, in particular, the lower bound excited singlet states, whichare of primary interest in photochemistry. This takes inspiration from recently developed Becke’sexciton model, where a key step constitutes the accurate evaluation of correlated singlet-tripletsplitting energy. It introduces a non-empirical model, both from “adiabatic connection theorem”and “virial theorem” to analyze the role of 2e − integral in such calculations. The latter quantityis efficiently mapped onto a real grid and computed accurately using a purely numerical strategy.Illustrative calculations are performed on 10 π -electron organic chromophores within a Cartesian-grid implementation of pseudopotential Kohn-Sham (KS) DFT, developed in our laboratory, takingSBKJC-type basis functions within B3LYP approximation. The triplet and singlet excitationenergies corresponding to first singly excited configuration, are found to be in excellent agreementwith TD-B3LYP calculations. Further, we perform the same for a set of larger molecular systemsusing the asymptotically corrected LC-BLYP, in addition to B3LYP. A systematic comparisonwith theoretical best estimates demonstrates the viability and suitability of current approach indetermining optical gaps, combining predictive accuracy with moderate computational cost. Keywords:
Density functional theory, adiabatic connection theorem, virial theorem, Becke’sexciton model, range-separated hybrid functional, singlet-triplet splitting. ∗ Email: [email protected], [email protected]. . INTRODUCTION Over the past decades, Kohn-Sham density functional theory (KS-DFT) [1] has emergedas a very powerful and successful tool for ground-state calculations of many-electron sys-tems, such as atoms, molecules and solids [2]. Moreover, its time-dependent (TD) variant,generally entitled linear-response (LR)-TDDFT has been developed for computing excita-tion energies [3, 4], mostly for discrete spectrum. In comparison to the conventional wavefunction-based formalisms, such as configuration interaction (CI), multireference CI, com-plete active space self-consistent field, equation-of-motion coupled cluster etc., LR-TDDFThas gained popularity due to its reasonable trade-off between accuracy and efficiency. Thisis evident from its applicability to medium and large systems [5]. However, despite its hugesuccess, it has well-known difficulties regarding double excitation, charge transfer and Ry-dberg excitation. This is mainly due to the adverse effects of exchange-correlation (XC)potential within adiabatic approximation. A detailed analysis could be found in [6].Besides TDDFT, attempts were also made to obtain excited states within a time-independent DFT rubric–several elegant approaches are available in the literature [7–11].In particular, the ∆SCF method [12, 13]was developed within the standard self-consistentfield iteration, by adopting non-Aufbau occupations at each iterations, and it targets thenon-Aufbau solutions using the ground-state functional. It has favorable scaling like ground-state DFT, and hence smaller computational resource compared to LR-TDDFT methods.Although it already provides a good estimation of excitation energies for molecular systems[14], its tendency to variationally collapse to the ground state is a serious concern. Suchfeatures during the ∆SCF procedure has been characterized in recent times [15]; this isparticularly severe in systems with a dense energy spectrum near the Fermi energy levels.Several sophisticated schemes were put forth to alleviate this issue, which include a fewconstrained-DFT formalisms [16–18] as well as gentlest ascent dynamics and meta dynamicsrelated methods [19], etc. These are often somehow involved, but are quite successful whereTDDFT fails to perform well. Furthermore, a lot of recent developments have been placedin the field of unconstrained excited-state orbital relaxation that have the same complexityas normal ground-state or TDDFT, while using a simple ∆SCF-style ansatz [20–24]. On theother hand, orbital energies themselves have physical meaning as excitation energies [25],but are quite sensitive to the choice of density functionals and its implementations. In a2ecent article [26], KS-DFT has been used successfully to extract the photo-electron andelectronic excitation spectrum for molecules from orbital energies through range-separated(RS) XC functionals. Moreover, efforts are also known the development of multi-reference(MR)-DFT by means of configurations mixing such as, multiconfiguration pair DFT [27]and MR-DFT with generalized auxiliary systems [28].In a series of articles, Becke [29–32], introduced a simple novel model for estimation oflowest single-particle excitation energies via correlated singlet-triplet splitting (STS) energy.It offers an accurate, economical way to compute optical gaps in large molecules. Themain focus here lies in the lowest singlet excited state, which are of fundamental interestin photochemistry. While triplet states can be easily obtained, in principle , standard DFTcannot be used in a singlet excited state due to its multi-determinantal nature. The presentscheme advocates two separate single-determinant DFT calculations—one for the closed-shell ground state and another for open-shell lowest triplet excited state, followed by a simpletwo-electron integral (Coulomb self energy) evaluation corresponding to the HOMO-LUMOtransition. We note that there is no concern about standard density functionals appliedon a given triplet excited state as it is represented by a single Slater determinant and ischaracterized by Fermi hole. Thus the main ingredient is the correlated STS energy whichcan be approached by means of (i) adiabatic connection theorem [29] and (ii) virial theorem[30]. In a sense, this avoids the configuration mixing and is altogether non-empirical.In the present work, we have adopted the above approach for lowest single-particle excita-tion energy corresponding to singlet excited state in molecules. The required computationsare performed following a pseudopotential KS-DFT implemented in Cartesian coordinategrid (CCG), as developed in our laboratory [33–41]. The pertinent two-electron integral(from appropriate orbitals obtained after solving the lowest triplet excited state) is carriedout numerically using a recently designed algorithm [41]. This employs the Fourier convolu-tion theorem in conjunction with a RS Coulomb interaction kernel. The latter is efficientlymapped onto the real grid through a simple grid optimization prescription, giving rise tosome constraints in RS parameter. The “adiabatic” model [29] corresponding to correlatedSTS energy is also provided accurately, leading to an easy route to compute the singletexcitation energy. Following the arguments of Becke [30], an elegant “virial” theorem isthen engaged to analyze the role of two-electron integral in determining the correlated STSenergy. Both procedures are applied to calculate optical gaps arising from both π → π ⋆ and3 → π ⋆ transitions in a decent number of molecules adopting B3LYP [42] functional. For aproper comparison, we carry out parallel calculations corresponding to lowest triplet excitedstates, from GAMESS quantum chemistry package [43] using TD-B3LYP method, whichreveals a better accuracy. Further, to extend the scope and applicability of this approxi-mation, we execute the same using “virial” theorem for a set of large molecular systems,and employing an additional XC functional from RS hybrid [44] family. The effectivenessof this scheme is illustrated by the respective statistical analysis. The article is organizedas follows. The underlying theorem of single-particle excitation energy along with a briefsummary of our KS-DFT framework in CCG, is provided in Sec. II. Section III offers thenecessary computational and technical details. Finally, the feasibility, performance and ac-curacy of our results are critically assessed in Sec. IV. Some concluding remarks as well asthe future prospects are given in Sec. V. II. METHODOLOGY
Let us consider an excitation of a given system, corresponding to an electronic configu-ration ϕ i ϕ f from a closed-shell ground state. Assuming completely filled closed-shell core,this is spanned by four Slater determinants: | ϕ αi ϕ αf i , | ϕ αi ϕ βf i , | ϕ βi ϕ αf i and | ϕ βi ϕ βf i , where α and β denote up and down spin. Therefore, the coupled excited states are found by diag-onalizing the Hamiltonian matrix in the space of above four determinants. As such, thesinglet state is given by | ψ S i = √ {| ϕ αi ϕ βf i − | ϕ βi ϕ αf i} , whereas the three degenerate tripletstates are written as follows: ψ T = | ϕ αi ϕ αf i or √ {| ϕ αi ϕ βf i + | ϕ βi ϕ αf i} or | ϕ βi ϕ βf i . Thecorresponding singlet and triplet energies (identified by “S” and “T” subscripts) are, E S = E αβ + K if (1)and E T = E αα or E αβ − K if or E ββ , (2)where E σ σ is the energy of a given determinant of form | ϕ σ i ϕ σ f i , ( σ , σ ) ∈ { α, β } , and K if is the 2e − integral (or Coulomb self-energy of product of transition orbitals) defined as, K if = Z Z ϕ i ( r ) ϕ f ( r ) ϕ i ( r ) ϕ f ( r ) | r − r | d r d r . (3)Now combining Eqs. (1), (2), one can connect singlet and triplet excitation energy as follows, E = E + 2 K if , (4)4here E = E S − E , E = E T − E , while E signifies the ground-state energy of acertain closed-shell system. But Eq. (4) is highly inaccurate for singlet excitation energy, E . The problem lies in the determination of 2 K if term, also called the zeroth-order (oruncorrelated) STS energy. In order to tackle this issue, recently some novel proposals (semi-empirical [29] as well as non-empirical [29, 30]) have appeared in the literature to derive asimplified formula for correlated STS energy. A. The adiabatic connection theorem
One can make use of the well-known “adiabatic connection” theorem [45] to approachthe correlated STS energy in single-particle excitations, which we are discussed here. Fora specific excited configuration, all the orbitals involved in Eqs. (1), (2) are same. In thatcase, singlet and triplet states have same density and non-interacting kinetic energy. Thentheir energy difference can be written as (∆ E = 2 K if ),∆ E STS = ∆ E + ∆ E corrST , ∆ E STS = 2 K if + ∆ E corrST , (5)where ∆ E corrST represents the singlet-triplet correlation energy difference. Recently, a non-empirical formula [29] has been proposed for it, based on the inter-electronic cusp condition,and its effect on electron correlation. Accordingly, it can be expressed as,∆ E corrST = − . Z ϕ i ( r ) ϕ j ( r ) z (cid:20) − ln(1 + z C ) z C (cid:21) d r . (6)The only unknown quantity, the correlation length z C , measures the spatial extent of electroncorrelation in configuration, ϕ i ϕ f . Now if one allows a “strictly correlated electrons” limit, z C can be written in terms of two-electron integral, K if , as in the following,0 . z Z ϕ i ( r ) ϕ j ( r ) d r = 2 K if . (7)A detailed derivation could be found in [29]. B. The virial exciton model
From the above discussion, it is clear that, one may add a correlation correction term,∆ E corrST , to 2 K if to recover the desired correlated STS energy. As suggested in [30], one5ay also further proceed to approximate this by noting that it comprises the kinetic andpotential energy contributions, ∆ E corrSTS = ∆ T corrSTS + ∆ V corrSTS . (8)Now invoking the standard virial theorem, the above equation can be simplified as,∆ E corrSTS = 12 ∆ V corrSTS . (9)One may further presume that the non-local part of ∆ V corrSTS is dominated by pair-densityeffects. In that scenario, correlation would then lower the potential energy of singlet staterelative to triplet state, by ∆ V corrSTS = − K if . Hence, one can write,∆ E corrSTS = − K if . (10)which surprisingly leads to a very simple relation as below,∆ E STS = 2 K if − K if = K if . (11)It may be noted that, this model (as well as that in II(A)) is also a purely non-empiricalone. But obviously, this expression for correlated STS energy is much simpler involving onlythe two-electron integral. C. The semi-empirical approach
As implied, both the routes in Eqs. (5), (11) are purely non-empirical, suggesting ∆ E STS to be less than 2 K if . This prompts one to formally define a molecule-independent re-scalingparameter f such that, ∆ E STS = 2 f K if , < f < . (12)and if determined semi-empirically, it might offer good excitation energies overall, as longas de-localization error [46] is not a serious issue. This led to a semi-empirical technique[29] to obtain an optimum value of f . This was accomplished in [29] by fitting the resultswith best estimated theoretical excitation energy data set of [47], resulting in a value for f as 0.486. It is surprisingly close to 0.5, as also appears from a consideration of II(B). Thesame author has also reported that the obtained mean average error (using this f ) remainssimilar to that in non-empirical calculations, in II(A).6 II. COMPUTATIONAL DETAILS
Our desired quantity, the singlet excitation energy E can be obtained from two single-determinant calculations (both are spin polarized), namely, one for closed-shell ground stateand another for fully relaxed triplet state. These are computed using our in-house pseu-dopotential KS-DFT program in CCG [48]. An initial version was first developed by thecorresponding author in 2008, on the basis of works [33–37], which has been later extendedby his group [38–41]. Since the details have been published in the above references, here wegive only some essential aspects. Accordingly, the single-particle KS equation in presence ofpseudopotential is written as (atomic unit employed unless stated otherwise), (cid:20) − ∇ + v pion ( r ) + v ext ( r ) + v h [ ρ ( r )] + v xc [ ρ ( r )] (cid:21) ϕ σi ( r ) = ǫ i ϕ σi ( r ) , (13)where v pion denotes the ionic pseudopotential, written as below, v pion ( r ) = X R a v pion , a ( r − R a ) . (14)In the above equation, v pion , a signifies ion-core pseudopotential associated with atom A, sit-uated at R a . The classical Coulomb (Hartree) term, v h [ ρ ( r )] describes usual electrostaticinteraction amongst valence electrons, whereas v xc [ ρ ( r )] represents the non-classical XC partof latter, and { ϕ σi , σ = α or β } corresponds to a set of N occupied orthonormal spinmolecular orbitals (spin-MOs).Now various quantities like localized basis function, two-electron potential, MOs andelectron density are directly set up on a 3D cubic box, r i = r + ( i − h r , i = 1 , , , ...., N r , r = − N r h r , r ∈ { x, y, z } , (15)where h r , N r denote grid spacing and total number of points along three directions respec-tively. While one-electron contributions of KS-Fock matrix are evaluated through well-established recursion relations, the two-electron matrix elements are computed via directnumerical integration in the grid, F hxc µν = h χ µ ( r g ) | v hxc ( r g ) | χ ν ( r g ) i = h x h y h z X g χ µ ( r g ) v hxc ( r g ) χ ν ( r g ) . (16)where v hxc refers to the classical (Hartree) and non-classical (XC) potential combined. Theconstruction of HF exchange in KS-Fock matrix has been documented in our earlier work[41]; hence not repeated here. 7he actual implementation of the central quantity, K if , is quite different from the originalprescription in [29]. Here, we use the Fourier convolution theorem (FCT) using a range-separated technique in Coulomb interaction kernel, instead of the multi-center numericalintegral procedure put forth in [49]. According to [41], this two-electron integral K if can berewritten as, K if = Z Z ϕ i ( r ) ϕ f ( r ) ϕ i ( r ) ϕ f ( r ) | r − r | d r d r = Z ϕ i ( r ) ϕ f ( r ) v if ( r ) d r . (17)Now, the main concern lies in the evaluation of v if integral which is related to electro staticpotential integral evaluation in [41]. Again, this can be recast as, v if ( r ) = Z ϕ i ( r ) ϕ f ( r ) | r − r | d r = Z ϕ if ( r ) | r − r | d r = ϕ if ( r ) ⋆ v c ( r ) (18)The last expression is in terms of convolution integral, where ϕ if denotes simple multiplica-tion of ith and fth MO from lowest triplet excited state and v c ( r ) represents the Coulombinteraction kernel. Now one can invoke the FCT to further simplify this integral, v if ( r ) = F − { v c ( k ) ϕ if ( k ) } where ϕ if ( k ) = F { ϕ if ( r ) } (19)Here v c ( k ) and ϕ if ( k ) stand for Fourier integrals of Coulomb kernel and MOs respectively.The main concern lies in an accurate mapping of the former, which has singularity at r → ζ ).This gives rise to the following expression, v c ( r ) = erf( ζ r ) r + erfc( ζ r ) r v c ( r g ) = v clong ( r g ) + v cshort ( r g ) . (20)In the above equation, erf( x ) and erfc( x ) denote error function and its complement re-spectively, while the second expression is written in real-space CCG. The Fourier integralof Coulomb kernel can be separated out as follows: (i) FT of short-range part is treatedanalytically and (ii) long-range portion is computed directly from FFT of correspondingreal-space CCG values. Then the remaining problem lies in finding an optimum value of8arameter ζ opt for successful mapping of Coulomb kernel in CCG from first principles . Inthis regard, we proposed a simple procedure which is found to be sufficiently accurate forlowest triplet excited states (as exemplified in results that follows) as well. This prompts usto write, ζ opt ≡ opt ζ E tot = opt N x ,N y ,N z E tot , at fixed h r , (21)using a suitably defined constraint [41] ( ζ × L = 7), where L (= N r h r ; r ∈ { x, y, z } ) isthe smallest side of simulating box. In the same spirit, other necessary quantities suchas correlation length, z C and singlet-triplet correlation energy difference, ∆ E C are directlycomputed in real-space grid using pseudo KS orbitals ϕ i and ϕ f .The strategy as outlined for HF exchange component, is implemented in case of B3LYP,containing a fixed amount of former with conventional DFT XC functional. As such, theXC energy corresponding to this functional is defined as follows, E B3LYPxc = (1 − a ) E xLSDA + a E xHF + a x E xB88 + a c E cLYP + (1 − a c ) E cVWN . (22)Throughout our current presentation, we use the recommended values of a , a x , a c , as advo-cated in [50], i.e., 0 .
2, 0 .
72, 0 .
81 for B3LYP. This is computed using KS orbitals obtainedfrom the solution of Eq. (13) in real-space CCG. The relevant LDA- and GGA-type function-als in connection with B3LYP are: (i) Vosko-Wilk-Nusair (VWN)–with the homogeneouselectron gas correlation proposed in parametrization formula V [51] (ii) B88–incorporatingBecke [52] semi-local exchange (iii) Lee-Yang-Parr (LYP) [53] semi-local correlation. Allthese are adopted from density functional repository program [54] except LDA.The present work employs SBKJC [55] effective core potential basis sets for species con-taining Group II elements. These are imported from EMSL Basis Set Library [56]. Thenorm-conserving pseudopotential matrix elements in contracted basis are collected fromGAMESS [43]. The triplet calculations refer to restricted open-shell, so that orbitals ϕ i and ϕ f are well defined. The convergence criteria imposed in this communication are slightlytighter than our earlier works [33, 34, 37]; this is to generate accurate orbital energies,especially for triplet excited states. Changes in following quantities were followed duringself-consistent field process, viz., (i) orbital energy difference between two successive itera-tions and (ii) absolute deviation in density matrix elements. They both were required toremain below a prescribed threshold set to 10 − a.u.; this ensured that total energy main-tained a convergence of at least this much, in general. In order to perform discrete Fourier9ransform, standard FFTW3 package [57] was invoked. The resulting generalized matrix-eigenvalue problem is solved through standard LAPACK routine [58] easily. All molecularcalculations are performed with calculated geometry (B3LYP XC functional and cc-PVTZbasis set), taken from NIST database [59] (otherwise stated below). Other details includingscaling properties may be found in references [33–41]. IV. RESULT AND DISCUSSION
A set of 10 selective organic chromophores ( π containing molecules) from our previouswork [41] has been chosen, excluding those which have degenerate frontier orbitals, as themethod in its present form is not applicable for these. These results pertain to the pseu-dopotential approximation neglecting the effects of core electrons. This may, by and large,be justifiable on the ground that only valence excitations are taken into account. By doingso, one can make a balance between accuracy and cost. To establish the effectiveness ofthis approach, we have chosen “SBKJC” type of pseudopotential basis set which does notcontain any diffuse function, and this suffices in our present scenario. It may be emphasizedthat, in this exploratory study, only lowest singlet and triplet excited state correspondingto first single excitation, is pursued for each molecule.Each calculation begins by optimizing the non-uniform grid for a given molecule, bothin its closed-shell ground and lowest excited triplet state separately, employing a simpleprocedure as reported earlier [38, 39]. Individual singlet excitation energies (in eV) throughB3LYP XC functional, are displayed in Table I, along with mean absolute error (MAE)and mean error (ME) statistics. It contains both non-empirical and empirical excitationvalues as following: (i) PR represents Eq. (12) using semi-empirical re-scaling parameter f = 0 .
486 (ii) PR presents Eq. (5) using non-empirical model from “adiabatic connection”(iii) PR refers to results from Eq. (11) using non-empirical model from virial theorem. Toput things in perspective, column 6 reports the respective TD-B3LYP energies employingsame functional and basis set, computed from GAMESS[43]. The correlation lengths z C manifestly fall into two distinct groups with π → π ⋆ and n → π ⋆ transition–a feature thathas been observed in [29]. Furthermore, in consonance with their finding, we also witnessa similar kind of pattern in the calculated z C values, having smaller ϕ i ϕ f overlap for the n → π ⋆ transition for current set of molecules. Overall, the obtained singlet excitation10 ABLE I: Singlet excitation energies (in eV) using B3LYP functional. See text for details.
Molecule State PR PR PR TD-B3LYP [43] z C Ethylene B u ( π → π ⋆ ) 7.78 7.63 7.87 8.09 2.97Propene A ′ ( π → π ⋆ ) 7.18 7.05 7.26 7.81 2.991,3-Butadiene (E) B( π → π ⋆ ) 5.63 5.42 5.70 6.02 3.291,3-Pentadiene (E) A ′ ( π → π ⋆ ) 5.47 5.28 5.54 5.88 3.291,3,5-Hexatriene (E) B u ( π → π ⋆ ) 4.38 4.14 4.44 4.79 3.532,4-Hexadiene (E,E) B u ( π → π ⋆ ) 5.39 5.20 5.45 5.79 3.271,3-Cyclo-pentadiene A ′ ( π → π ⋆ ) 5.12 5.03 5.17 5.28 2.98Thiophene B ( π → π ⋆ ) 5.61 5.31 5.66 6.02 3.99Formaldehyde A ( n → π ⋆ ′′ ( n → π ⋆ ) 4.67 4.75 4.68 5.07 1.44MAE – 0.40 0.53 0.34 – –ME – 0.40 0.53 0.34 – – energies are within fraction of eV compared to the TD-B3LYP results. As expected from[30], the correlated virial theorem (PR ) confers remarkable improvement over the adiabaticconnection formula (PR ). Here, the MAE value is about two times smaller compared tothe latter one. On the other hand, the semi-empirical approach (PR ) performs very closeto PR results, as evidenced from the proximity in their corresponding MAE values. Theworst case candidate is propene, with an excitation energy 0 .
55 eV too low, in comparisonto TD-B3LYP result. We note that the MAE and ME remain same in magnitude and sign,indicating a systematic error in calculated values. Out of three methods, the virial PR routeprovides the best estimate in exciting energies; quite competitive with that of TD-B3LYP,with an MAE of 0.34 eV.Next we analyze the effects of triplet excitations and correlated STS term separatelyfor B3LYP functional in Table II, giving a side-by-side comparison with TD-B3LYP in bothcases. As singlet excitation energies are obtained from these two above quantities, a detailedassessment of them is worthwhile considering. Due to the superior performance, henceforthin future proceedings, we engage only the PR . The E error metrics are: MAE = 0.40eV and ME = − results of previoustable. On the other hand, STS measures read: MAE = 0.70 eV and ME = 0.69 eV, showingvery close agreement in magnitude and also importantly without changing sign. Howevermagnitude of these errors (MAE and ME) in STS are relatively larger than the E T errors.11 ABLE II: Triplet excitation energies and correlated STS energies (in eV) using B3LYP functional.
Molecule State E ∆ E STS
Ref. [48] TD-B3LYP PR TD-B3LYPEthylene B u ( π → π ⋆ ) 4.47 4.03 3.40 4.06Propene A ′ ( π → π ⋆ ) 4.44 4.03 2.82 3.781,3-Butadiene (E) B( π → π ⋆ ) 3.26 2.71 2.44 3.311,3-Pentadiene (E) A ′ ( π → π ⋆ ) 3.24 2.71 2.30 3.171,3,5-Hexatriene (E) B u ( π → π ⋆ ) 2.42 1.85 2.02 2.942,4-Hexadiene (E,E) B u ( π → π ⋆ ) 3.22 2.70 2.23 3.091,3-Cyclo-pentadiene A ′ ( π → π ⋆ ) 3.21 2.70 1.96 2.58Thiophene B ( π → π ⋆ ) 3.88 3.47 1.78 2.55Formaldehyde A ( n → π ⋆ ) 3.19 3.20 0.34 0.78Acetaldehyde A ′′ ( n → π ⋆ ) 4.39 4.44 0.29 0.24MAE – 0.40 – 0.70 –ME – − A careful look on respective ME in these two quantities reveals that, ME in STS is about75% larger compared to that in E T , but with an opposite sign. And that makes the singletexcitation energies more close to TD-B3LYP results, occurring through a systematic errorcancellation. This leads us to conclude that, STS term and not E , is the major source oferror. It is worthwhile mentioning that, the success of this approach relies on an accurateestimation of the two-electron integral, which in turn depends on the accuracy of tripletexcited state. Thus it is evident that our current prescription can be reliably used fordetermining optical gaps using the non-empirical time-independent methods of Sec. II, atleast for small molecular systems.At this stage, in order to assess the validity and usefulness, an additional set of largermolecular systems are taken up in Tables III and IV. We consider a representative set of16 organic chromophores and 5 linear acene molecules with 2–6 rings for which accuratetheoretical reference data is available [60]. All the molecular geometry of Table III are takenfrom [60] whereas for Table IV, these are optimized with B3LYP and cc-pVDZ basis set usingGAMESS. Since the present scheme has earlier been verified to deliver practically identicalresults as that of GAMESS, for brevity, we have employed the same for calculation of E and E T in the same basis. In addition to B3LYP, another XC functional from RS hybridfamily, has been invoked. Accordingly, the inter-electronic space is separated through the12 ABLE III: Excitation energies (in eV) of organic chromophores from “virial theorem”.
Molecule State E E (PR ) Lit. ‡ B3LYP LC-BLYP B3LYP LC-BLYPCyclopropene B ( π → π ⋆ ) 4.03 4.05 7.04 7.07 7.01Norbornadiene A ( π → π ⋆ ) 4.62 4.23 5.77 5.45 4.91Naphthalene B u ( π → π ⋆ ) 3.22 3.53 4.74 5.20 4.64Furan B ( π → π ⋆ ) 4.15 4.25 6.65 6.85 6.57Pyrrole B ( π → π ⋆ ) 4.45 4.52 6.88 7.07 6.85Pyridine B ( n → π ⋆ ) NC 4.22 NC 4.65 4.74Pyridazine B u ( n → π ⋆ ) 2.76 2.78 3.35 3.34 3.57s-tetrazine B ( π → π ⋆ ) 1.39 1.63 1.83 2.07 2.15Acetone A ( n → π ⋆ ) 3.71 3.80 3.98 4.08 4.11Formamide A ′′ ( n → π ⋆ ) NC 5.03 NC 5.25 5.42Acetamide A ′′ ( n → π ⋆ ) 5.20 5.15 5.45 5.37 5.46Propanamide A ′′ ( n → π ⋆ ) 5.64 5.19 6.07 5.40 5.49Thymine A ′ ( π → π ⋆ ) 3.51 3.59 5.50 5.73 5.48Uracil A ′′ ( n → π ⋆ ) 4.28 3.76 5.66 5.86 4.66Imidazole A ′ (Π → π ⋆ ) NC 4.75 NC 7.15 6.89E-Octatetraene B u ( π → π ⋆ ) 1.93 1.92 4.32 3.69 4.22MAE from Lit. results[30] . 0.27 0.23ME from Lit. results [30]. − − − − a MAE and ME values from RO-PBE0 calculations using 6-31G ⋆ are 0.37 and − b MAE and ME values from RO-LC ω PBE0 calculations using 6-31G ⋆ are 0.31 and − † NC denotes “not converged”. ‡ This corresponds to E , from [30]. use of a RS-operator g ( γ, r ) such that,1 r = ˜ g ( γ, r ) r + g ( γ, r ) r , (23)where ˜ g ( γ, r ) represents the complementary RS-operator. Here, γ is an important RS param-eter which has a pivotal role to adjust the contribution of HF exchange between short-range(SR) and long-range (LR) region, for a specific g ( γ, r ). In particular, we consider the long-range correction (LC) scheme presented in [61], given by, E LCxc = E dfa , sr ( γ ) + E x , lr ( γ ) + E c , dfa ,g ( γ, r ) = erf( γ r ) and ˜ g ( γ, r ) = erfc( γ, r ) . (24)It uses full HF exchange at LR region ( E x , lr ) in contrast to the conventional hybrid func-tionals, and successfully mitigates the lack of having correct asymptotic behaviour of exact13 ABLE IV: Excitation energies (in eV) of linear acenes from “virial theorem”.
Number of rings E E (PR ) Lit. † B3LYP LC-BLYP B3LYP LC-BLYP2 3.24 3.56 4.75 5.22 4.653 2.22 2.46 3.71 4.43 3.584 1.55 1.81 2.85 3.41 2.755 1.08 1.32 2.30 2.96 2.226 0.75 0.98 1.89 2.50 1.82MAE from Lit. results [30] 0.09 0.70ME from Lit. results [30] − − − − † This corresponds to E , from [30]. exchange at LR region present in hybrid family. But at the same time, it fails to achieve thestandard energy calculations at the level of hybrid ladder. To be consistent with B3LYP,we utilize the particular LC-BLYP (with γ = 0 . E dfa , sr ( γ ) and E c , dfa respectively.The calculated E , as computed above, are compared with the all-electron results fromBecke’s “virial theorem” method [30] along with respective MAE and ME values for twoconcerned functionals at the bottom. As in Table II, here also, due to the proximity of PR results to that of TD-B3LYP, it suffices to report only PR excitation energies. For sakeof completeness, this also accompanies the lowest triplet excitation energies. Note that, forthree molecules (pyridine, formamide, imidazole), the spin-restricted triplet calculations didnot lead to convergence in the default options, for B3LYP. Hence they are not included inMAE and ME evaluations. It is evident that the performance of B3LYP (excitation energiesare overestimated) is usually much more consistent than that of LC-BLYP. A systematicerror cancellation may be involved in these calculation. On the other hand, the MAE valueof LC-BLYP (0.23) is slightly lower than that of B3LYP (0.27), but did not outperformextensively as was expected from the foregoing discussion. Initially, one may draw theconclusion that the effect of full HF exchange at long-range does not have significant impacton excitation energies, despite the fact that it enriches the behavior of frontier orbitalsinvolved in K if calculations. Such discrepancy has occurred due to the fact that we assume γ
14s independent of system size. However, it has been shown that it is possible to attain a higherlevel of performance provided one treats γ as a system-dependent parameter (functional of ρ ) tuned from first principles [44]. It is our belief that, an optimally tuned γ (perhaps in thespirit of size dependency principle) will provide better results than the conventional hybridand RS hybrid functionals. Note that, the less dramatic performance of LC-BLYP in bothground and triplet excited state has negligible effect on singlet excitation energy, relative tothe STS term, which is in harmony with the observation made in [31]. Additionally, we havealso estimated MAE and ME values from TBE-2 results [47]. It is found that these in TableIII (with respect to TBE-2), do not alter significantly from their counterparts in Becke [30],which was of course, expected. So the basic conclusions will remain same. Furthermore,we have also quoted restricted open-shell KS (ROKS) error statistics (MAE and ME) fordifferent XC functionals in footnote of Table III from [21]. They are all-electron ROKS lowestvertical excitation results computed using 6-31G ⋆ basis set. Our performance (MAE andME from TBE-2) are in harmony with ROKS MAE and ME values. A similar qualitativetrend is also found for linear acenes, given in Table IV, where again both E as well as E ,are given. In this occasion, the better compatibility of MAE and ME values signify that theerror is systematic in nature in case of LC-BLYP. But these are always overestimated andthe extent is larger than B3LYP. This may arise due to a greater involvement of γ on systemsize ; hence possibly requires an optimally tuned γ in theses cases. This demands furtherelaborate investigation, and may be considered in future works. V. FUTURE AND OUTLOOK
We have examined the feasibility and practicability of a simple yet accurate time-independent DFT approach for realistic calculation of single excitation energies throughDFT in CCG. This was applied to a host of organic chromophores and linear acenes having π network. The excitation energies were estimated using two representative set of functionals–each from respective hybrid and RS hybrid within a pseudopotential approximation. Theobtained results for all these species, from virial theorem are in quite good agreement withthe reference results. It can be easily extended to the excitation energies in charge transfercomplexes [32, 63] as well. Its application with RS hybrid, hyper, as well as local hybridXC functionals would further enhance the scope in a wide range of systems. It may also be15esirable to examine its performance in various excited configurations other than the lowestexcited state. But this may pose some computational challenge regarding the state-specificSCF convergence. In such occasions, some solution, as suggested in [17, 22–24, 64], may offersome valuable guidelines. To conclude, the present work has demonstrated the usefulness ofa simple scheme in predicting optical gap in organic chromophores within a CCG-DFT. VI. ACKNOWLEDGEMENT
AG is grateful to UGC for a senior research fellowship. TG acknowledges INSPIREprogram for financial support. AKR thankfully acknowledges funding from DST SERB,New Delhi, India (sanction order: EMR/2014/000838). [1] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[2] A. D. Becke, J. Chem. Phys. , 18A301 (2014).[3] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997 (1984).[4] M. E. Casida, in Recent Advances in Computational Chemistry; D. P. Chong (Ed.) (WorldScientific, NJ, 1995), pp. 155–192.[5] M. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem. , 287 (2012).[6] N. T. Maitra, J. Chem. Phys. , 220901 (2016).[7] A. K. Theophilou, J. Phys. C , 5419 (1979).[8] E. K. U. Gross, L. N. Oliveira, and W. Kohn, Phys. Rev. A , 2809 (1988).[9] R. Singh and B. M. Deb, J. Chem. Phys. , 5892 (1996).[10] A. K. Roy, R. Singh, and B. M. Deb, J. Phys. B , 4763 (1997).[11] P. Ayers, M. Levy, and ´A. Nagy, J. Chem. Phys. , 191101 (2015).[12] T. Ziegler, A. Rauk, and E. J. Baerends, Theor. Chem. Acc. , 261 (1977).[13] T. Kowalczyk, S. R. Yost, and T. V. Voorhis, J. Chem. Phys. , 054128 (2011).[14] I. Seidu, M. Krykunov, and T. Ziegler, J. Phys. Chem. A , 5107 (2014).[15] I. Seidu, M. Krykunov, and T. Ziegler, Mol. Phys. , 661 (2014).[16] T. Ziegler, M. Krykunov, and J. Cullen, J. Chem. Phys. , 124107 (2012).[17] G. M. J. Barca, A. T. B. Gilbert, and P. M. W. Gill, J. Chem. Phys. , 111104 (2014).
18] P. Ramos and M. Pavanello, Phys. Chem. Chem. Phys. , 21172 (2016).[19] C. Li, J. Lu, and W. Yang, J. Chem. Phys. , 224110 (2015).[20] I. Frank, J. Hutter, D. Marx, and M. Parrinello, J. Chem. Phys. , 4060 (1998).[21] T. Kowalczyk, T. Tsuchimochi, P.-T. Chen, L. Top, and T. Van Voorhis, J. Chem. Phys. ,164101 (2013).[22] J. A. Shea and E. Neuscamman, J. Chem. Phys. , 081101 (2018).[23] D. Hait and M. Head-Gordon, J. Chem. Theory Comput. , 1699 (2020).[24] K. Carter-Fenk and J. M. Herbert, J. Chem. Theory Comput. (2020).[25] R. Van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Theory Comput. , 4432 (2014).[26] R. L. A. Haiduke and R. J. Bartlett, J. Chem. Phys. , 131101 (2018).[27] G. Li Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar, and L. Gagliardi,J. Chem. Theory Comput. , 3669 (2014).[28] Z. Chen, D. Zhang, Y. Jin, Y. Yang, N. Q. Su, and W. Yang, J. Phys. Chem. Lett. , 4479(2017).[29] A. D. Becke, J. Chem. Phys. , 194107 (2016).[30] A. D. Becke, J. Chem. Phys. , 044112 (2018).[31] A. D. Becke, J. Chem. Phys. , 081102 (2018).[32] A. D. Becke, S. G. Dale, and E. R. Johnson, J. Chem. Phys. , 211101 (2018).[33] A. K. Roy, Int. J. Quant. Chem. , 837 (2008).[34] A. K. Roy, Chem. Phys. Lett. , 142 (2008).[35] A. K. Roy, in Handbook of Computational Chemistry Research , edited by C. T. Collett andC. D. Robson (Nova publishers, New York, 2009).[36] A. K. Roy, Trends in Phys. Chem. , 27 (2010).[37] A. K. Roy, J. Math. Chem. , 1687 (2011).[38] A. Ghosal and A. K. Roy, In Specialist Periodical Reports: Chemical Modelling, Applicationsand Theory; M. Springborg and J.-O. Joswig (Eds.) Vol. 13 (Royal Society of Chemistry,London, 2016).[39] A. Ghosal, T. Mandal, and A. K. Roy, Int. J. Quant. Chem. , e25708 (2018).[40] T. Mandal, A. Ghosal, and A. K. Roy, Theor. Chem. Acc. , 10 (2019).[41] A. Ghosal, T. Mandal, and A. K. Roy, J. Chem. Phys. , 064104 (2019).[42] A. D. Becke, J. Chem. Phys. , 5648 (1993).
43] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Hensen,S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, et al., J. Comput. Chem. , 1347 (1993).[44] R. Baer, E. Livshits, and U. Salzner, Annu. Rev. Phys. Chem. , 85 (2010).[45] J. Harris and R. Jones, J. Phys. F , 1170 (1974).[46] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr., Phys. Rev. Lett. , 1691 (1982).[47] M. R. Silva-Junior, M. Schreiber, S. P. Sauer, and W. Thiel, J. Chem. Phys. , 174318(2010).[48] A. K. Roy, A. Ghosal and T. Mandal, InDFT: A DFT Program for Atoms and Molecules inCCG (Theoretical Chemistry Laboratory, IISER Kolkata, India, 2019), This is based on theextension of an initial version of the code, established by A. K. Roy, in 2008, whose resultswere published in refs. [26-30].[49] A. D. Becke and R. M. Dickson, J. Chem. Phys. , 2993 (1988).[50] P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Chem. Phys. , 11623(1994).[51] S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. , 1200 (1980).[52] A. D. Becke, Phys. Rev. A , 3098 (1988).[53] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B , 785 (1988).[54] Density Functional Repository, Quantum Chemistry Group (CCLRC Daresbury Laboratory,Daresbury, Cheshire, UK, 2001).[55] W. J. Stevens, H. Basch, and M. Krauss, J. Chem. Phys. , 6026 (1984).[56] D. Feller, J. Comp. Chem. , 1571 (1996).[57] M. Frigo and S. G. Johnson, Proceedings of the IEEE , 216 (2005).[58] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Dongarra, J. D. C. A. Greenbaum, S. Ham-marling, A. A. McKenney, and D. Sorensen, LAPACK Users’ Guide , vol. 9 (SIAM, 1999).[59] R. D. Johnson,
III (Ed.) NIST Computational Chemistry Comparisons and BenchmarkDatabase, NIST Standard Reference Database, Number, Release 18 (NIST, Gaithersburg,MD, 2016).[60] M. R. Silva-Junior, M. Schreiber, S. P. Sauer, and W. Thiel, J. Chem. Phys. , 174318(2010).[61] H. Iikura, T. Tsuneda, and K. Hirao, J. Chem. Phys. , 3540 (2001).[62] S. Grimme and M. Parac, Chem. Phys. Chem , 292 (2003).
63] X. Feng, A. D. Becke, and E. R. Johnson, J. Chem. Phys. , 231101 (2018).[64] A. T. B. Gilbert, N. A. Besley, and P. M. W. Gill, J. Phys. Chem. A , 13164 (2008)., 13164 (2008).