Bond-breaking excitations with diverging coupling matrix of response density functional theory from highest-level functionals

Bond-breaking excitations ωα are the problematic case of adiabatic time-dependent density functional theory (TDDFT). The calculated ωα erroneously vanishes with the bond elongation, since the Hartree-exchange-correlation kernel and the corresponding response coupling matrix K of standard approximations lack the characteristic divergence in the dissociation limit. In this paper an approximation for K is proposed constructed from the highest-level functionals, in which both occupied and virtual Kohn-Sham orbitals participate with the weights wp. The latter provide the correct divergence of K in the limit of dissociating two-electron bond. The present K brings a decisive contribution to the energy of the 1Σu+ in the prototype H2 molecule calculated for various H-H separations. At shorter separations it improves ωα compared to the zero-order TDDFT estimate, while at the largest separation it reproduces near-saturation of the reference excitation energy.


Introduction
The key element of response density functional theory (DFT) is the Hartree-exchange-correlation (Hxc) kernel f Hxc (x 1 , x 2 , ω), which was rigorously derived in the paper of E. K. U. Gross and collaborators [1]. It produces the response coupling matrix K, which couples orbital transitions in the matrix diagonalization problem of timedependent DFT (TDDFT) [2]. In the customary adiabatic approximation this problem assumes the following form Here, E is the diagonal matrix of the differences between the energies a and i of the virtual φ a and occupied φ i Kohn-Sham (KS) orbitals. In (1) F α is the solution vector, which contains the real parts of the responses δγ R s,ai (ω α ) of the KS independent-particle first-order reduced density matrix (1RDM).
Typically, the zero-order TDDFT (2) fairly approximates the energies ω α of single excitations in compact molecules and localized excitations in molecular complexes [3]. Then, the inclusion of K evaluated with standard approximate DFT functionals further improves the quality of ω α calculated with (1). On the other hand, there are also well-known problematic cases of adiabatic TDDFT, such as the simultaneous vanishing of both ω α and ( a − i ) for the excitation corresponding to the bonding-antibonding orbital transition in the case of bond breaking [4,5]. As was established in [4], in this case of strong left-right electron correlation, in order to provide the true finite ω α , the exact Hxc-kernel and the corresponding element K ai,ia diverge in the dissociation limit.
In this paper, to resolve the above mentioned problematic case of adiabatic TDDFT, it is proposed to construct the response coupling matrix K from an approximate effective electron-electron interaction energy functional W [{φ i }, {φ a }], which can be attributed to the highest level of the "Jacob's ladder" of DFT [6]. Unlike functionals of lower levels, the highest-level functionals depend not only on the occupied KS orbitals φ i or on the electron density ρ obtained from φ i , but also on the virtual KS orbitals φ a . Previously, functionals of this level were employed to properly bring strong left-right correlation to calculation of ground-state potential energy curves [7,8].
In Section 2 a hybrid response approach is proposed with K obtained from the functional W [{φ i }, {φ a }]({w}), in which both occupied and virtual KS orbitals φ p participate with the weights w p . With this, the functional naturally accounts for strong left-right correlation. In Section 3 the expression for K is derived with the orbital differentiation of W [{φ i }, {φ a }]({w}) applied previously in the orbital extension [9,10] of time-dependent density matrix functional theory (TDDMFT) [11][12][13]. The divergence of the derived K for dissociating two-electron bond is demonstrated. In Section 4 the problematic case of adiabatic TDDFT is illustrated with the comparative calculations of the 1 Σ + u excitation along the bond-breaking coordinate of the prototype H 2 molecule. In Section 5 the matrix K derived from the adapted to the KS orbitals Löwdin-Shull (LS) functional [14] is applied together with the model potential of van Leeuwen and Baerends (LB) [15] to calculation of the problematic excitation. The correction from K produces the proper qualitative effect within the LBKα scheme, keeping ω α from a rapid decline with R(H − H). In Section 6 the conclusions are drawn and further development of the present methodological strategy is envisaged.

Strong left-right correlation with highest-level functionals
In this paper we consider a hybrid approach to the eigenvalue problem (1) and (2), in which the matrix E is composed from the energies p of the conventional KS orbitals φ p . They are the solutions of the one-electron KS equations.
whereĥ(r) is the one-electron operator (with the external potential v ext ) and v Hxc is the KS Hxc potential.
On the other hand, the response coupling matrix K in (1) is obtained within the adiabatic approximation with double differentiation with respect to the orbitals of the following highest-level functional Here, both occupied and virtual KS orbitals from (3) participate with the weights w p in the first term, the one-electron energy (h pp is the diagonal matrix element of the operator (4) for the orbital φ p ), as well as as in the second term, the effective electron-electron interaction energy W [{φ i }, {φ a }]({w}). The weights w p are chosen to provide the minimum of the functional (5) for the KS orbitals φ p .
With its dependence on orbitals and their weights, the GKS functional (5) has the same form as functionals of density matrix functional theory (DMFT) [14,[16][17][18][19][20][21][22][23][24][25] depending on the natural orbitals (NOs) χ p and their occupations n p . Because of this, the former naturally describes the considered strong left-right electron correlation in a dissociating two-electron bond. In this particular case the relevant bonding occupied φ g and the antibonding unoccupied φ u KS orbitals tend to the corresponding NOs χ g and χ u , while the weights w g and w u together with the NO occupations n g and n u approach the dissociation limit 1/2. The left-right correlation can be captured already in the minimal two-orbital model, which involves the frontier orbitals φ g , φ u as well as their weights w g and w u . In this model, the adapted to the KS orbitals LS functional [14] of DMFT can be employed, the electron-electron interaction energy part of which reads where φ p φ q |φ r φ s is the two-electron integral in its physicists' notation In the limit of the dissociating two-electron bond all twoelectron integrals in (6) tend to the same value I, while w u approaches w g . As a result, the energy (6) vanishes in this limit, thus correctly describing two individual one-electron atoms. Because of the above mentioned formal similarity of the highest-level functional (5) and DMFT functionals, the orbital differentiation developed previously in DMFT and TDDMFT will be applied in Section 3 to obtain the response coupling matrix K.

Diverging response coupling matrix from highest-level functionals
In this section the expression for K is presented, which is derived within the orbital extension of the TDDMFT [9,10]. This extension and the KS approach of TDDFT both consider the orbital changes in response to the time-dependent external perturbation δv ext (rt). Here, U is the evolution matrix the antihermiticity of which is required to ensure orthogonality of the perturbed orbitals. The orbital changes (8) induce the change in the electron-electron interaction which, in a general case, is represented with the following matrix equation [10,26] δv ef f ai (t) = jb K ai,jb δU bj (t).
Here, v ef f ai is the matrix element of the effective potential, which in the case of TDDFT is the Hxc potential v Hxc ai , while in the case of TDDMFT this is the potential of the electron-electron interaction v ee ai . To evaluate from (10) the elements of the coupling matrix K ai,jb in the adiabatic approximation, the latter potential is used, which is obtained with the orbital differentiation of the functional (5). The orbital derivatives of (5) enter in v ee ai in the form of the matrix W and its adjoint W † ai , so v ee ai reads [10,26] Then, the adiabatic coupling matrix K is obtained with further orbital differentiation of (12), so K ai,jb are expressed through the projected orbital derivatives of v ee ai as follows [10,26] . (13) Inserting (12) in (13), we obtain the final general expression for K For the introduced two-orbital model the following explicit expression for the element K ug,gu is derived from (6), (11), (12), and (14) As was already mentioned in Section 2, in the case of strong left-right correlation in the limit of the dissociating two-electron bond all two-electron integrals tend to the same value I, while w u approaches w g . As a result, K has the following asymptotics so it correctly diverges due to its vanishing denominator. In Section 4 the coupling matrix (15) will be used to evaluate the energy of the problematic excitation along the bond-breaking coordinate.

Problem of the bond-breaking excitation
To illustrate the problem of adiabatic TDDFT with bondbreaking excitations, Table 1 presents the energies ω α of the lowest 1 Σ + u excitation in the prototype H 2 molecule calculated at four different R(H − H) separations with TDDFT with the xc functionals BLYP and B3LYP. The former is one of the standard functionals of generalized gradient approximation (GGA) with the exchange functional of Becke (B88) [27] and the correlation functional of Lee, Yang, and Parr (LYP) [28,29], while B3LYP is the related hybrid functional [30]. The BLYP ω BLY P and B3LYP ω B3LY P excitation energies are compared with the reference values ω F CI obtained with the full configuration interaction (FCI) method. All calculations are carried out in the augmented with polarization functions correlationconsistent triple-zeta aug-cc-pVTZ basis [31,32] with the GAMESS-US program package [33,34].
The characteristic feature of standard DFT approximations is a severe artificial upshift of their KS orbital energies [35]. In particular, the exact energy of the highest occupied KS orbital is equal to (minus) the first vertical ionization potential (VIP) I p [36][37][38][39][40][41]. At variance with this, the BLYP and B3LYP energies 1g of the occupied KS orbital of H 2 are much higher than VIPS −I F CI p calculated with FCI (see Tab. 1).
Apparently, since the energies 1u of the lowest occupied KS orbital are also severely upshifted, the corresponding difference ∆ , the zero-order TDDFT estimate (2), qualitatively correctly vanishes with R(H − H). Note, that due to the admixture of the Hartree-Fock (HF) functional in B3LYP, B3LY P 1g is consistently lower than BLY P 1g , while B3LY P 1u is higher than BLY P 1u , so that ∆ B3LY P is consistently larger than ∆ BLY P . However, while the contribution of the BLYP response coupling matrix always produces an appreciable increase of the resultant excitation ω BLY P compared to ∆ BLY P , the B3LYP coupling matrix produces a smaller increase at the longer R(H − H) = 3 and R(H − H) = 4 Bohr, while at shorter separations it, actually, decreases ω B3LY P compared to ∆ B3LY P . As a result, ω BLY P is rather close to ω B3LY P at all R(H − H).
The above mentioned positive contributions of the BLYP and B3LYP coupling matrices at longer R(H − H) cannot compensate vanishing ∆ , so ω BLY P and ω B3LY P experience artificial too fast decline at these separations.

LBKα scheme
Can the proposed in this paper response coupling matrix K prevent the demonstrated fast decline of the bondbreaking excitation calculated by adiabatic TDDFT? In order to answer this question, K is combined within the LBKα scheme with the approximate xc potential of van Leeuwen and Baerends (LB) [15]. The latter is employed in the present paper in the form where v LDA x (r) is the LDA exchange potential [42], v V W N c (r) is the LDA correlation potential in the parametrization of Vosko, Wilk, and Nusair (VWN) [43], and B x (β, γ; r) is the B88-type correction to the energy density of the LDA exchange functional B x (β, γ; r) = − βρ(r) 1/3 x(r) 2 1 + γβx(r)sinh −1 (x(r)) , with x(r) being the dimensionless gradient-dependent argument Because of the choice of the parameter γ = 3, v LBα xc has the correct long-range asymptotics With the present parameters α = 1.19 and β = 0.01, the potential (17) and (18) represents the outermost component of that with statistical average over (different) orbital potentials (SAOP) [44]. The latter was developed with the special purpose to fairly approximate the orbital energies of the accurate KS potential [45]. The solution of equation (3) for H 2 is obtained in the same aug-cc-pVTZ basis. One can see from = −0.130 Hartree. Yet, the difference ∆ LBα , the zero-order TDDFT LBα estimate of ω α , is rather close to ∆ BLY P , due to the compensation of the artificial upshifts of the orbital energies in the latter.
Within LBKα, the difference ∆ BLY P is combined with the element K ug,gu in the small matrix approximation to (1) [46] The element K ug,gu is composed according to (15) from the LBα orbitals φ 1g and φ 1u , while their weights w g and w u are obtained from the minimization of the energy (5) with the LS functional (6) with the additional requirement w g + w u = 1. One can see from Table 2, that the resultant weights are rather close to the occupations n 1g and n 1u of the corresponding FCI NOs. The correction from K has a major influence on the calculated excitation energy ω LBKα (see Tab. 1). This correction greatly improves the quality of ω LBKα compared to the zero-order estimate ∆ LBα . At the separations R(H − H) = 1.4−3 Bohr the resultant ω LBKα is close to the reference ω F CI . At R(H − H) = 4 Bohr the proposed coupling matrix produces the required qualitative effect, preventing the bond-breaking excitation calculated within adiabatic TDDFT from the fast decline and reproducing near-saturation of ω F CI .

Discussion and conclusions
In this paper it is proposed to evaluate the adiabatic response coupling matrix K from the GKS energy functional, in which both occupied and virtual KS orbitals are involved through the orbital weights. The matrix K is to be combined with the conventional KS orbital energy matrix E in the response matrix diagonalization problem and the concrete combination LBKα is considered.
The present K naturally incorporates the effect of strong left-right correlation and this is demonstrated in a two-fold way. First, a derivation is given of the correct divergence of K in the limit of the dissociating twoelectron bond. Second, it is shown that, within small matrix approximation, K brings a decisive contribution to These results demonstrate the principal applicability of the proposed K within the response matrix diagonalization problem (1). Further development will use reliable highest-level functionals (5) also for the general N -electron case. Then, the corresponding coupling matrix K of (14) will exhibit a similar divergence as in the considered two-electron case, if the bond breaking in an N -electron system is similarly reflected in near-degeneracy of the relevant orbitals of the bonding and antibonding character. This near-degeneracy means, that the orbitals possess nearly equal energies and weights and the presence of the difference between the latter quantities in the denominator of (14) will cause the divergence of K. Just as in the two-electron case, this divergence is required, in order to counter the effect of near-degenerate orbital energies and to provide a finite bond-breaking excitation energy.
Moreover, the present methodological strategy can also be applied to tackle another problem of adiabatic TDDFT, the lack of double excitations. Based on the recent unified description of TDDFT, time-dependent Hartree-Fock (TDHF) theory, and time-dependent density matrix functional theory (TDDMFT) [26], the consistent differentiation of the properly chosen GKS functional can produce the following response matrix diagonalization problem Here, the solution vector F α reads where γ R (ω α ) is the sub-vector of the real off-diagonal elements of the response of the GKS density matrix and δw(ω α ) is that of w p changes. In (22) and (23) the matrix E is an extension of the orbital energy matrix E of (2), while D is the following compound matrix Here, Ω ss is the block responsible for coupling of single excitations. It is an extension of the matrix (E + 2K) in (1) and it contains the derivatives of the generic functional with respect to the KS orbitals. The block Ω dd is responsible for coupling of double excitations and it contains the second derivatives of the generic functional with respect to the orbital weights w p . Finally, the block Ω sd is responsible for coupling of single and double excitations and it contains mixed second derivatives with respect to both KS orbitals and their weights. The corresponding extension of the present approach is in progress. Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.