Hole–hole Tamm–Dancoff-approximated density functional theory: A highly efficient electronic structure method incorporating dynamic and static correlation

Linear response time-dependent density functional theory (TDDFT) is the most commonly used electronic structure method for excited states due to its low computational cost and relative accuracy for absorption spectra involving valence excited states.¹ Unfortunately, due to the difficulties that presently available approximate exchange-correlation (XC) functionals encounter with near-degeneracies and static correlation, it is unsuitable for photochemical problems involving a conical intersection between the ground and first excited states.² While complete active space self-consistent field (CASSCF) methods can treat static electron correlation, they struggle to describe dynamic correlation.³ Adding a correction to recover dynamic correlation (as in multireference perturbation theory⁴ or multireference configuration interaction⁵) significantly increases the computational cost of the method and renders nonadiabatic dynamics simulations computationally intractable for many interesting medium to large sized molecules (although such dynamics is feasible for small molecules^6–12).

There has long been interest in the development of inexpensive excited-state methods that can simultaneously treat dynamic and static correlation. One approach to this problem is to augment traditional multiconfigurational methods with semiempirical parameters. These parameters can be as simple as a constant scaling parameter applied to the energy¹³ or as complex as a full semiempirical treatment of the Hamiltonian matrix elements.^14,15 While often successful, these methods may require cumbersome system-specific parameterization.^15–17 Alternatively, conventional multiconfigurational wavefunctions can be combined with a Kohn–Sham (KS) density functional theory (DFT)-based treatment of dynamic electron correlation.^18,19 Many of these methods use conventional multiconfigurational wavefunctions to capture static correlation and include a DFT-based treatment of dynamic correlation. Methods of this type include range-separated wavefunction/DFT methods,^20,21 multiconfigurational pair-density functional theory^22,23 (MC-PDFT), and combinations of DFT and configuration interaction,^24,25 such as in DFT/MRCI.^26,27 Alternative approaches attempt to incorporate multireference character directly in the DFT formalism by using an unconventional reference²⁸ or appealing to ensemble densities. Spin-flip TDDFT (SF-TDDFT)²⁹ is an example of a scheme based on an unconventional reference. Spin contamination problems have plagued DFT-based spin–flip methods;^30,31 however, this has recently been addressed by combining SF-TDDFT with ad hoc corrections from DFT/MRCI.³² Ensemble formulations of DFT include the spin-restricted ensemble-referenced Kohn–Sham (REKS) methods.³³ To date, REKS methods need to be formulated specifically for the chosen ensemble (defined by an active space, as in CASSCF), and a general formulation applicable to arbitrary active spaces is lacking. Current implementations including gradients and nonadiabatic couplings are only applicable to an active space of two electrons in two orbitals.^34,35

The recent development of particle–particle random phase approximation (pp-RPA) methods has opened new possibilities for inexpensive excited-state methods. The original aim of pp-RPA was to provide ground state correlation energies via an adiabatic connection fluctuation dissipation theorem (ACFDT) approach.^36–41 While the particle–hole RPA (ph-RPA) ACFDT approach recovers the ring channel of the correlation energy from the coupled cluster doubles (CCD) equations,⁴² the ladder channel of the CCD correlation energy is obtained from the pp-RPA ACFDT approach.^36,41 Yang and co-workers also highlighted pp-RPA and its Tamm–Dancoff approximated pp-TDA variant as effective methods to compute electronic ground and excited state energies.^43–52 Starting from a doubly cationic (N−2)-electron reference, the N-electron ground state and excited states generated by excitations from the highest occupied molecular orbital (HOMO) are recovered by performing two-electron attachments.^45,49,50 This allows the treatment of the N-electron ground and excited states on equal footing (derived as simultaneous eigenvalues of a common Hamiltonian) at a computational cost comparable to the simplest excited state methods, e.g., TDDFT/ph-RPA and configuration interaction singles (CIS). Because the ground and excited states are treated on equal footing, pp-TDA based on an (N−2)-electron reference is able to predict the correct topography around conical intersections, as has been shown explicitly for H₃ and NH₃.⁵³ This is a major advance over conventional ph-TDDFT and CIS methods, which cannot reproduce conical intersections involving the ground state.² In 2014, Peng et al. derived the pp-RPA equations from linear response theory by choosing a pairing field perturbation (termed TDDFT-P, see below in Sec. II A).⁵⁴ This formally allowed for the combination of pp-RPA and DFT references. Hence, the effect of dynamic correlation on the orbitals is incorporated through the exchange-correlation (XC) potential, while the pp-RPA and pp-TDA schemes ensure that ground and excited states are treated on equal footing and therefore can treat exact degeneracies correctly (see also Sec. II D).

In a procedure complementary to pp-TDA and pp-RPA based on an (N–2)-electron reference, the N-electron ground state and excited states can also be generated through double annihilations from a double anionic (N+2)-electron reference in which the lowest unoccupied molecular orbital (LUMO) is populated with two additional electrons. This (N+2)-electron pp-RPA scheme and its corresponding hole–hole (hh) Tamm–Dancoff approximation (hh-TDA) were first presented by Yang and co-workers.⁴⁹ Although the pp-RPA and pp-TDA methods based on an (N–2)-electron reference have now become quite established,^43–53 less attention has been paid to the hh-TDA method based on an (N+2)-electron reference. Yang and co-workers applied the hh-TDA method to oxygen and sulfur atoms⁴⁹ (noting “relatively large errors” in the results), but we have found no published reports of further developments or applications of hh-TDA to molecules. In this paper, we suggest that the hh-TDA method is worthy of renewed attention. Similar to the pp-TDA method, hh-TDA can effectively capture dynamic and static correlation in ground and low-lying excited states, including near- and exact degeneracies. Furthermore, the active orbital space in the (N+2)-electron-based hh-TDA appears to be suitable for the description of excited states in many organic molecules that are inaccessible to (N−2)-electron pp-TDA. Examples include molecules with low-lying nπ^* and ππ^* excited states that cannot be simultaneously described within pp-TDA (which, by construction, can only describe excited states where an electron is excited from the HOMO).

In contrast to previous work on pp-RPA, for both pp- and hh-TDA, we formulate the response kernel in a functional-specific way that resembles the kernels occurring in linear density matrix response theory. This new formulation of the response kernel is compared to the previously used functional-independent variant. We put particular emphasis on the utility of the hh-TDA method in the treatment of organic and biologically relevant systems involving both ππ^* and nπ^* transitions.

The particle–particle random phase approximation (pp-RPA) equations were derived by Peng and co-workers⁵⁴ by means of coupled time-dependent perturbation theory in analogy to well-established linear response theory by choosing a pairing field perturbation within the framework of Hartree–Fock (HF)/Kohn–Sham–Bogoliubov theory (termed TDDFT-P).⁵⁵ The ground state for this non-interacting particle system is defined by the zero-temperature grand potential as^54–57 

Here, T_s[γ, κ] is the independent particle kinetic energy and V_ext[γ] is the external potential energy, which contains the nuclear-electron attraction. D_ext[κ] denotes the external pairing potential, and the last term preserves the total electron number. E_JXC[γ, κ] is the mean-field potential energy due to the electrons and includes the particle–hole and particle–particle channels via the one-particle density matrix γ and the pairing matrix κ, respectively. The indices J, X, and C denote the Coulomb, exchange, and correlation components of this functional, respectively. The one-particle density matrix in the canonical molecular spin orbital basis is given by

while the pairing matrix (or anomalous density matrix) is defined as

Ψ refers to the single-determinant ground state wavefunction of the non-interacting Kohn–Sham (KS) system. Here and in the following, i, j, k, l denote occupied, a, b, c, d denote unoccupied, and p, q, r, s denote general molecular orbitals.

For the hypothetical true functional E_JXC[γ, κ], the real-space density and anomalous density of the non-interacting KS system are identical to the respective quantities of the true interacting system. In the work of Peng et al.,⁵⁴ a non-superconducting system with no external pairing field in the ground state Hamiltonian is used. In that case, D_ext[κ] = 0. As a result, the mean-field potential energy is free from indirect electron–electron interactions (e.g., phonon-mediated interactions)⁵⁵ and contains only direct electron–electron interaction contributions via the (anti-symmetrized) Coulomb operator and XC potential. Based on the ground state defined by Eq. (1), Peng et al. derived the pp-RPA method by perturbing the ground state with an external pairing field. The TDDFT-P equations, which describe the coupled response to the pairing field perturbation, are given as⁵⁴ 

Here, δκ_pq(ω) are the first order changes in the pairing matrix induced by the external frequency-dependent pairing field perturbation $δDpqext(ω)$⁠, ω is the frequency of the perturbation, ϵ_p is the orbital energy, and μ is the chemical potential or Fermi energy. The pairing matrix response L_pq,rs of the mean field potential is contained in $δDpqJXC(ω)$⁠,

In our formulation of pp-RPA, the choice of L_pq,rs differs from the one presented by Yang and co-workers,^44,45,54 as discussed in detail below.

For a non-superconducting system in the absence of D_ext[κ], the pairing matrix and its associated XC potential are zero in the electronic ground state. Consequently, only the γ-dependent terms in Eq. (1) survive, and minimizing the zero-temperature grand potential becomes equivalent to minimizing the standard KS or Hartree–Fock (HF) energy expression. It is only after taking the second derivatives with respect to κ that the particle–particle channels of the mean-field potential as in Eqs. (4) and (5) become non-zero. In the framework of ab initio wavefunction theory, the pp-RPA equations are based on a HF reference, and the mean-field potential takes the form of an anti-symmetrized Coulomb integral.⁴¹ In that case, the pp and particle–hole (ph) channels of the mean-field potential are formally equivalent. We note, however, that only the exchange-type integral survives the spin integration in the pp channel (see the supplementary material).

In the context of approximate KS DFT, there may exist some liberty in choosing the explicit form of this response kernel. Given that contemporary semi-local density functional approximations (DFA) are defined for κ = 0, the pairing matrix response of the density functional drops out,

Yang and co-workers have already employed this approximation for the XC density functional response. Furthermore, they chose to use the HF-type mean-field response, i.e., a bare anti-symmetrized Coulomb integral, as the response kernel.⁵⁴ It is therefore independent of the underlying density functional approximation (DFA). With this kernel, the chosen DFA only affects the orbitals and orbital energies that enter the pp-RPA calculation, but not the pp response expression itself.

In this work, we employ a different choice for the pp response kernel. Given the formally equivalent ph and pp response kernels when using a HF reference,⁴¹ which may be regarded as a special choice of the KS system, we draw analogy from the well-known ph linear response TDDFT. In spin-preserving^1,58 and spin–flip (SF)²⁹ formulations of TDDFT, the response kernel reflects the underlying ground state DFA in the way the non-local Fock exchange enters the response. Hence, we choose to use the same modification of the non-local Fock exchange, i.e., the global scaling and/or range-separation employed by the corresponding DFA,

The first term on the right-hand side corresponds to the Coulomb-type two-electron integral (Mulliken notation for spin MOs), whereas the remaining terms denote the modified exchange-type two-electron integrals. $aXFR$ is the functional-specific scaling factor for the global [full-range (FR)] exchange integral,

Here, ϕ refers to a molecular spin orbital and x₁/x₂ denote all spatial and spin coordinates of electrons 1/2. The short range-attenuated (SRA) part of the exchange integral employs the modified Coulomb operator and is given as

Like Yang and co-workers, we also assume the pairing matrix response of the XC functional to be zero [see Eq. (7)], allowing the response kernel to become frequency independent. Due to spin integration (note the definition of the pairing matrix), the Coulomb part in the pp response kernel is zero; hence,

Here, we have used the Mulliken notation for spatial MOs and overbarred orbital indices to highlight orbitals corresponding to β spin. Our choice for the response kernel [Eq. (11)] is identical to the one in collinear SF-TDDFT but with inverted sign [note that L_pq,rs is subtracted in Eq. (14) below]. This can be understood from a second quantization picture in which both the spin–flip and pp-RPA approaches apply creation/annihilation operators acting on different spin, i.e., $a↑†a↓$ in spin–flip and $a↑†a↓†$/$a↑a↓$ in the pp/hh channels. Therefore, only the exchange integral does not vanish upon spin integration, as in pp-RPA, and the flipped sign results from the commutation relation of the annihilation and creation operators (see the supplementary material for a comparison). SF-TDDFT can be derived from linear-response ph TD-DFT if changes in the density matrix that do not preserve the $Ŝz$ expectation value are allowed.²⁹ Due to its origin from an analogy to linear response SF-TDDFT, we refer to our choice of the response kernel, Eq. (11), as the linear response-type kernel in the following to distinguish it from the DFA-independent, HF-like kernel used by Yang and co-workers.⁴⁹ 

Using the response kernel in Eq. (11), we can rewrite Eqs. (4) and (5) in matrix notation as

Here, the superscripts pp and hh define the particle–particle (double creation in the virtual space) and hole–hole (double annihilation in the occupied space) blocks, respectively. A^pp is of dimension $nvirt2×nvirt2$⁠, and A^hh is of dimension $nocc2×nocc2$⁠. Their elements are given as

and

B^hp = [B^ph]^T is of dimension $nocc2×nvirt2$ and is given as

If we now let the perturbation [right-hand side in Eq. (12)] go to zero, we find that there are non-trivial solutions to this system of equations, which correspond to the poles of the particle–particle/hole–hole response function. This leads to the non-Hermitian eigenvalue problem

Employing the commonly used notation,⁵⁴ the matrices X and Y contain the eigenvectors, which correspond to the pp and hh changes in the pairing matrix κ. ω is a diagonal matrix, which contains the eigenvalues of the non-Hermitian eigenvalue problem. If we neglect the coupling matrices B^hp and B^ph between the pp and hh blocks, we obtain the so-called Tamm–Dancoff approximated (TDA) eigenvalue problems,

for the pp-TDA case and

for the hh-TDA problem.

Note that within TDA, the constant diagonal shift of $∓2μ$ in Eqs. (13) and (14) can be removed from the definition of the corresponding matrix elements. Both of these eigenvalue problems are guaranteed to have only real eigenvalues. This makes them more robust than the full pp-RPA case, which can have complex solutions that spoil the potential energy surfaces (PES) around conical intersections (as noted in the context of coupled cluster theory^59,60). While the pp-TDA method based on an (N–2)-electron reference has been extensively investigated as a method to describe low-lying excited states^45,49 and S₀/S₁ conical intersections,⁵³ the hh-TDA method has been less thoroughly explored for the molecular electronic structure.⁴⁹ In this work, we suggest that hh-TDA should be revisited as an efficient DFT-based method capable of computing low-lying excited states even in the presence of conical intersections.

Before benchmarking the hh-TDA method along with our choice for the response kernel [Eq. (11)], we will briefly discuss some practical aspects of the pp-TDA and hh-TDA methods in Sec. II D.

In both the pp-TDA and hh-TDA methods, the N-electron target system is described by first solving the ground state electronic structure for a system that differs by two electrons from the N-electron target system. This corresponds to a double cation (N–2 electrons) in the case of pp-TDA or a double anion (N+2 electrons) in the case of hh-TDA. The ground and excited states of the N-electron target system are then obtained by creation (pp-TDA) or annihilation (hh-TDA) of two electrons (see Fig. 1).

The inability to describe either doubly excited states or conical intersections involving the ground state are well-known shortcomings of standard single-reference methods such as linear-response (ph) TD-DFT and CIS. These are both rectified by pp- and hh-TDA. Since the ground and the excited states are solutions of the same eigenvalue problem, near- and exact degeneracies (such as conical intersections) can be properly described. Furthermore, some doubly excited states can be computed (see Fig. 2). At the same time, the preceding ground state calculation of the double ionic state is of single-reference complexity and incorporates dynamic correlation by virtue of the XC potential. Therefore, both pp- and hh-TDA schemes are, in principle, able to capture static and dynamic correlation, making them promising methods for use in nonadiabatic dynamics simulations. Furthermore, in contrast to the SF-TDDFT scheme, pure eigenfunctions of $Ŝ2$ are obtained trivially in both pp- and hh-TDA.

However, both methods also have some obvious shortcomings in common. First, the excited state expansion space is highly restricted. Taking the N-electron target system as reference, pp-TDA only includes excitations from the HOMO, while hh-TDA is restricted to excitations to the LUMO. Second, both methods take a detour by computing the orbitals for a reference that differs in its total charge from the target state, as SF-TDDFT does with a differing spin state. As a result, two-state degeneracies of the N-electron state can be treated at the cost of a single reference method. However, the orbitals are not optimized for the N-electron system. This may have a significant implication for dynamics simulations: a degeneracy on the (N+2) or (N−2) surface itself can occur and cause instabilities in the SCF procedure.

Zhang et al. noted that the choice of a (N+2) reference may be hampered by the existence of unbound orbitals.⁵⁰ Presumably, this is one of the reasons that the hh-TDA method has not received serious attention in molecular electronic structure so far. There are two distinct consequences of unbound orbitals that can be envisioned. The first consequence is quite practical, namely, that it may be difficult to converge the SCF procedure for the (N+2) reference. In the usual hh-TDA method, non-convergence of the SCF procedure for the (N+2) reference would be fatal. The second consequence is somewhat more formal. Even if SCF convergence can be achieved, orbitals with positive orbital energies do not correspond to bound electrons. They can thus be very sensitive to the employed basis set. Indeed, with a sufficiently flexible basis set, they would be expected to be highly delocalized continuum functions and poor representatives of a low-lying valence excited state.

Regarding the first consideration, we note that poor convergence of the SCF procedure will be exacerbated if the long-range potential is incorrect. DFAs with 100% long-range Fock exchange are known to capture the asymptotic potential correctly. Accordingly, using DFAs with 100% long-range exchange, we have not observed any instability of the (N+2)-electron SCF procedure for the molecules we tested. Even the green fluorescent protein chromophore (HBI) anion with a net charge of −3 for the reference state converges without difficulty when an asymptotically correct, range-separated DFA is used. However, the hh-TDA method is likely to encounter serious SCF convergence difficulties for DFAs that do not incorporate long-range exact exchange.

It has been previously observed that finite-basis set DFT is capable of providing reasonable electron affinities computed as energy differences EA = E_N − E_N+1, even when the anion has unbound occupied orbitals.⁶¹ This suggests that the second consideration is largely formal. However, there are also some practical concerns. The success of hh-TDA in describing excited states is predicated on the HOMO orbital of the (N+2)-electron determinant being a good approximation to the bound orbital that predominantly contributes to the low-lying N-electron excited states. This will certainly not hold if the HOMO is a good approximation to a continuum orbital, which is the expected outcome in a sufficiently large basis set. We take a practical approach here and recommend that diffuse basis sets should be avoided in the context of hh-TDA. For the finite orbital basis sets of double- or triple-zeta quality used in this study, there are very few cases where the hh-TDA scheme based on an (N+2)-electron reference fails, in spite of the presence of occupied orbitals with positive energies. This implies that the HOMO of the (N+2)-electron system is a good approximation to the LUMO of the N-electron system, albeit with shifted orbital energies. This shift is effectively removed during the hh-TDA step, leading to a reasonable description of the relative energies between the N-electron ground and low-lying excited states. We will further elaborate on this issue in Sec. IV E, but we emphasize again that diffuse basis sets should usually be avoided in combination with the (N+2)-electron reference calculation.

For pp-TDA, there may also exist a practical implication with respect to the underlying DFA. It was shown that pp-RPA based on a HF reference leads to less accurate singlet-triplet splittings,⁴⁸ slower basis set convergence,⁴⁵ and overall less accurate excitation energies for small molecules⁴⁵ than pp-RPA based on DFAs with little or no exact exchange. Similar results for excitation energies were reported for pp-TDA.⁴⁵ A convenient test case is the ππ^* state of ethylene, where both hh-TDA and pp-TDA provide a suitable active orbital space. Accurate methods such as coupled cluster and multireference perturbation theory agree that this should be the lowest valence excited state. For ethylene, we observe a stronger dependence on the amount of Fock exchange in the DFA for pp-TDA (Table S13 in the supplementary material). The ππ^* state is the fourth excited state for pp-TDA-HF, while it is the second excited state for hh-TDA-HF. It is well-known that adding Fock exchange to a DFA increases the differential treatment of virtual and occupied orbitals. In molecular systems,⁶² the virtual orbitals in HF are best suited to describe electron attachment, while virtual orbitals in (semi-)local DFT are more appropriate to describe electronic excitations.⁶³ This implies that DFAs with large amounts of exact exchange will raise the energy of the HOMO of the N-electron target system in pp-TDA. Therefore, when using practical DFAs, the optimal amount of exact exchange could be different for pp-TDA and hh-TDA, and this point deserves further study. In this work, we focus on DFAs for the linear-response-type kernel hh-TDA framework.

We will note (details in Table S13) that the ordering of the ππ^* state for ethylene in pp-TDA is improved by using our linear response-type kernel compared to the previously introduced functional-independent HF-type response kernel.^45,49,54 For clarity, it should be noted that purely (semi-)local DFAs have a vanishing response kernel in our linear-response-type formulation, i.e., A^pp and A^hh are diagonal in that case. In the HF-type response kernel, the bare anti-symmetrized Coulomb integral is used as the response kernel, regardless of the underlying DFA.

In hh-TDA, within the limitations of the given basis set (see above), the LUMO (with respect to the N-electron target system) is generated on equal footing with the occupied orbitals for HF or any DFA. Another difference between hh- and pp-TDA is the dimension of the A matrix to be diagonalized. In pp-TDA, this dimension increases strongly with the size of the atomic orbital (AO) basis. For hh-TDA, the dimensions are independent of the basis set size (unless effective core potentials are used).

We mention for the sake of completeness that pp-TDA-HF is equivalent to full configuration interaction (FCI) for any two-electron system, while hh-TDA-HF is equivalent to FCI for an N-electron system in a basis that provides exactly two virtual molecular spin orbitals. We have confirmed this via calculations on H₂ in an STO-3G basis set (see the supplementary material). Thus, both schemes have similarities to CASCI methods with a fixed active space (i.e., N electrons in N/2 + 1 orbitals for hh-TDA and two electrons in n_virt + 1 orbitals for pp-TDA). One can speculate that the ability of the pp-TDA and hh-TDA to describe static correlation arises from their similarity to CASCI methods. For the same reason, however, neither pp-TDA nor hh-TDA is size consistent. As long as no degeneracy of the LUMO/LUMO+1 (hh-TDA) or the HOMO/HOMO–1 (pp-TDA) is present, this is not expected to have practical implications for the description of relative energies between electronic states.

The hh-TDA and pp-TDA methods have been implemented in the electronic structure code TeraChem.^64,65 All hh-TDA and pp-TDA calculations in this work are performed with this development version of TeraChem. To compare with the pp-TDA data in the literature, we have implemented our choice for the response kernel [Eq. (11)] as well as the HF-type response kernel of Yang and co-workers.⁵⁴ We confirmed that our implementation is correct by comparison of pp-TDA excitation energies with values reported in the literature⁴⁵ as well as comparisons between CASCI and hh/pp-TDA using Hartree–Fock reference states.

In the functional assessment for vertical excitation energies with hh-TDA in Sec. IV A, we considered the global hybrid functionals B3LYP,^66–69 PBE0,^70,71 and BHLYP^66,67,72 as well as Hartree–Fock. Thus, we cover different amounts of non-local Fock exchange in the mean field Hamiltonian (20%, 25%, 50%, and 100%, respectively). Additionally, we test different range-separated hybrid functionals, namely, CAM-B3LYP,^66–69,73 ωPBEh,⁷⁴ and the B97⁷⁵ type functionals: ωB97,⁷⁶ ωB97X,⁷⁶ and ωB97X-D3.⁷⁷ For clarity, it is emphasized that the latter two are different in their functional parameters (see the supplementary material) and hence should give rise to different electronic structures. The DFA assessment is restricted to the linear response-type kernel variant of hh-TDA [Eq. (11)] For the assessment of different functionals, we employed the spherical split-valence atomic orbital def2-SV(P)^78,79 basis set by Ahlrichs and co-workers. All SCF calculations used the converged Hartree–Fock orbitals as guess orbitals. The systems considered in this work are given in the supplementary material.

When comparing pp-TDA and hh-TDA in Sec. IV B, we use the same spherical TZVP basis set⁸⁰ that has been used in the original benchmark papers by Thiel and co-workers.⁸¹ Both schemes are employed in a setting, in which the occupied and active orbitals have been generated on equal footing, i.e., where these orbitals experience a mean-field potential with the same number of electrons. The pp-TDA is hence combined with the generalized gradient approximation (GGA) functional PBE⁷⁰ (see also Sec. II D). hh-TDA is combined with the range-separated functional ωB97X,⁷⁶ which is shown to be a well-performing combination in Sec. IV A (and in the supplementary material). We consider both the linear response-type (this work) and the HF-type (Yang and co-workers) response kernels for pp-TDA and hh-TDA. Additionally, we have performed ph Tamm–Dancoff approximated TDDFT calculations with the ωB97X functional. For reference, we use the best estimates for excitation energies of Thiel and co-workers⁸¹ and complement these with equation-of-motion singles and doubles coupled cluster (EOM-CCSD) calculations (the same basis set) as implemented in Q-Chem v5.0.0.⁸² 

In Sec. IV C 1, we calculate the ethylene potential energy surface (PES) with hh-TDA-ωB97X (linear response kernel) using the Cartesian def2-SVP^78,79 basis set. The scans along the pyramidalization and torsion angles were performed without relaxation of the remaining degrees of freedom.

The excited state energies of thymine at three critical structures are considered in Sec. IV C 2. We use hh-TDA in combination with the ωPBE functional and the Cartesian 6-31G^** basis set.^83–85 A range separation parameter ω = 0.2 a.u. was selected after a coarse scan and comparison against high level reference values for this molecule. We calculated these reference values with EOM-CCSD/aug-cc-pVDZ with Q-Chem. While the Franck–Condon (FC), the S₁ minimum, and the S₁/S₂ minimum energy conical intersection (MECI) were optimized with hh-TDA, we used the respective coupled cluster (CC) geometries of Ref. 86 for single-points at the EOM-CCSD/aug-cc-pVDZ level of theory. The S₂ (ππ^*) saddle point geometry therein is claimed to be “in close proximity” to the conical intersection seam.⁸⁶ Due to the lack of a properly CC-optimized MECI structure, we use that structure as a substitute and estimate the MECI energy by averaging the S₁ and S₂ energies computed with EOM-CCSD/aug-cc-pVDZ.

In Sec. IV C 3, we apply hh-TDA with the ωB97⁷⁶ functional in combination with the Cartesian 6-31G^** basis set^83–85 and compare to results obtained with multireference configuration interaction methods including single and double excitations (MRSDCI).⁸⁷ 

The performance of density functionals for the calculation of vertical excitation energies in the linear-response ph time-dependent density functional theory framework is well known⁸⁸ and the role of Fock exchange is fairly well understood.¹ In contrast, for the hh-TDA method proposed in this work, the impact of Fock exchange and the use of an anionic reference on the excitation energies is not known. Therefore, we start by benchmarking several common DFAs along with Hartree–Fock (HF) in the calculation of vertical excitation energies with hh-TDA.

To gain more insight, we consider the lowest vertical excitation energies using molecules and reference data from previously published datasets where highly accurate excitation energies are available.^81,89,90 We classify the excitation type into different categories: intermolecular charge-transfer (CT) and predominantly local excitations. We also consider molecules that have push–pull type excitations, i.e., excitations with partial intramolecular CT character. We find that the functionals behave similarly for this set as for the local excitations. Thus, we have added them to the latter set and distinguish only intermolecular CT and intramolecular excitations (see the supplementary material for separate results of the local and the push–pull set). In the latter set, we have selected molecules for which hh-TDA, across different DFAs, produces the same character for the lowest vertical excitation as the reference method (SCS-CC2⁹¹).

We have not considered any purely semi-local GGA-type DFAs since these show severe convergence problems in the self-consistent field (SCF) procedure of the double anionic reference (see also Sec. II D). Although we also expect that global hybrids will not be optimal in that case, for completeness, we have included some prototypical global hybrid functionals that differ in the amount of Fock exchange. All functionals considered are listed in Sec. III, but we restrict the discussion in this paper to the globally constant Fock-exchange DFAs B3LYP, BHLYP, and HF (these functionals contain 20%, 50%, and 100% of Fock exchange, respectively). Furthermore, we find that all tested range-separated DFAs with 100% long-range exchange perform similarly; therefore, we restrict the current discussion to the ωB97X-D3 functional. The results for the other functionals that were tested can be found in the supplementary material.

Figure 3 shows the Gaussian error distributions based on the mean deviation (MD) and standard deviation (SD) obtained for the aforementioned reference datasets of excitation energies. While hh-TDA-HF yields CT excitation energies that are on average underestimated by almost 1 eV, the excitation energies for locally excited states show a systematic overestimation of about the same magnitude. Furthermore, the error is significantly more systematic for the intermolecular CT excitations. For DFAs with low amounts of Fock exchange, i.e., B3LYP, these trends are reversed, however, with a smaller magnitude for the underestimation of local excitation energies. It is noteworthy that the DFAs behave quite different than in a ph linear response TD-DFT framework. In the latter, DFAs with low amounts of Fock exchange significantly underestimate CT excitations.⁸⁹ The BHLYP functional, which is in between in terms of the percentage of Fock exchange, shows more consistent errors for both sets (MD ≈ 0.4 eV–0.5 eV), as well as more narrow error distributions. The range-separated DFA ωB97X-D3 behaves even slightly better than BHLYP. This, along with the fact that the 100% asymptotic Fock exchange makes these DFAs less prone to SCF convergence issues with the (N+2)-electron reference, suggests the use of asymptotically correct range-separated hybrid functionals in combination with the hh-TDA scheme.

In addition to benchmarking the DFAs for CT and local states, we investigate state splittings between ππ^* and locally excited states with different character (mostly nπ^*). For this purpose, we have selected a set of systems mostly comprised of molecules from Thiel’s 2008 benchmark set.^81,90 Problematic systems, which are impossible to treat with hh-TDA due to the restricted orbital space, were discarded from the benchmark set. These are systems in which either the population of the frontier orbitals in the (N+2)-electron system is not clear (degenerate LUMOs in the N-electron determinant) or the excited states cannot predominantly be described by a single orbital–orbital transition (benzene, naphthalene, pyrrole, and s-triazine). Furthermore, butadiene, hexatriene, and octatetraene are not considered here but in Sec. IV B. The results are visualized in Fig. 4, while the respective symmetry labels and results for other density functionals can be found in the supplementary material. First, we stress the relevance of high amounts of Fock exchange in the DFA, which we mentioned earlier in Sec. II D. We find that B3LYP and BHLYP, but also a few range-separated functionals (CAM-B3LYP or ωPBEh, see supplementary material), show incorrect orbital occupations for some carbonyl systems in the (N+2)-electron SCF calculation, i.e., the HOMO of the (N+2)-electron system does not correspond to the LUMO obtained from an SCF calculation on the N-electron system. Hence, the excited states for these systems cannot be described properly with these functionals.

We also find that the global hybrid functionals perform poorly for the state splittings compared to Hartree–Fock or range-separated functionals (e.g., ωB97X-D3). Except in two cases (nπ^*/ππ^* in aspirin and σπ^*/ππ^* in cyclopropane), the latter two perform similarly. Both describe the nπ^* splittings reasonably well, while systematically overestimated ππ^* excitations are found for aromatic systems, regardless of which DFA is used. This seems to be caused by the absence of accessible higher lying anti-bonding π orbitals in the expansion space, which would lower these excitations and improve the energetic splittings with respect to the nπ^* excitations. This is also observed for the splittings between the different ππ^* states in cyclopentadiene and norbornadiene, presumably for the same reasons. Interestingly, reduced short-range Fock exchange leads to increased σπ^* excitation energies in cyclopropene, which then lead to qualitatively incorrect energy splittings with respect to the ππ^* state. This is observed for all DFAs, and only HF ranks these states qualitatively correctly.

Overall, the results indicate that the hh-TDA scheme in combination with a range-separated hybrid functional appears to be a promising method for describing the ground and low-lying excited states of organic molecules on equal footing.

In this section, we compare the hh-TDA and pp-TDA schemes as well as the two different response kernel choices. In Table I, we have tabulated the statistical data for excitation energies computed for a subset of excitations from the Thiel benchmark set⁸¹ with the linear response (LR) and HF-type kernels for hh-TDA-ωB97X and pp-TDA-PBE. We note that our pp-TDA-PBE (HF) results are in agreement with the ones presented earlier by Yang and co-workers (note that we use a slightly different basis set here).⁴⁵ The statistical data are listed in Table I, while detailed values of the individual systems and states are given in the supplementary material. A significant difference between the hh-TDA and pp-TDA schemes becomes apparent when looking at systems that involve both ππ^* and other types of excitations (i.e., nπ^* or σπ^*). Here, we find that hh-TDA typically provides the better option since the configuration space in hh-TDA is more suited to describe both excitation types simultaneously. This manifests itself in a larger number of states on this set, which can actually be described by hh-TDA. The nucleobases and carbonyl systems are particularly good cases for hh-TDA (see the supplementary material). For pp-TDA-PBE, SCF convergence problems for these systems eliminate the possibility to include them in the analysis. This seems to be due to (nearly) degenerate frontier molecular orbitals in the (N−2) reference and was confirmed using alternative electronic structure software (Turbomole version 6.4).^74,92 The HF-type variant performs better for pp-TDA-PBE than the LR variant, the latter being equivalent to plain PBE orbital energy differences. For hh-TDA-ωB97X, this is reversed and the LR-type kernel performs better, as reflected by smaller mean (absolute) and standard deviations. Both pp-TDA-PBE schemes on average show underestimated excitation energies, whereas the opposite is true for hh-TDA-ωB97X. From all of the considered pp/hh schemes, hh-TDA-ωB97X (LR) performs the best on this set.

Compared to regular ph-TDA-ωB97X and EOM-CCSD, hh-TDA-ωB97X (LR) provides almost comparable accuracy, with a slightly less systematic error distribution as reflected by a SD, which is twice as large as for ph-TDA-ωB97X. In summary, we find that hh-TDA provides a reasonable method typically with a broader applicability to many organic systems compared to pp-TDA. In Sec. IV C, we will only consider the LR-type kernel variant of hh-TDA and investigate potential energy surface properties beyond the Franck–Condon point.

As addressed in Sec. II D, one of the main advantages of the hh-TDA method is the ability to accurately treat degeneracies involving the ground state and low-lying excited states. The topologically correct description of the potential energy surfaces in regions of degeneracy is especially important in nonadiabatic photochemical dynamics where population transfer between adiabatic potential surfaces is often mediated by these seams of conical intersection. A prototypical example of such an intersection can be accessed through torsion and pyramidalization degrees of freedom in ethylene; this particular case has previously been used as a model system for assessing the ability of electronic structure methods to describe statically correlated character and ground-excited state degeneracies.² 

Typical low-cost single reference excited state methods, such as CIS and TDDFT, suffer from instability related to HOMO–LUMO orbital degeneracies as an MECI involving the ground state and an excited state is approached. As a result, these methods cannot properly describe the twisted-pyramidalized MECI in ethylene. By contrast, hh-TDA is well-suited to describing conical intersection topologies, as demonstrated in combination with the ωB97X functional in Fig. 5(a), in which the double cone topology characterizes the vicinity of the S₀/S₁ MECI. Here, the potential energy surfaces are plotted along torsion and pyramidalization coordinates, which are the branching plane coordinates corresponding to the twisted-pyramidalized minimum energy conical intersection (MECI) between the S₀ ground state and the ππ^* excited state.

Moreover, CIS and TDDFT cannot describe double excitations, which then leads to a description of the S₁ minimum at a purely twisted geometry (90° torsion, no pyramidalization).² Methods that incorporate doubly excited states describe the S₁ minimum as simultaneously twisted (90° torsion) and pyramidalized. hh-TDA captures this feature of the S₁ potential energy surface (PES), as demonstrated by the appearance of the minimum at a twisted (90°) and pyramidalized (60°) geometry that lies 4.76 eV above the ground state minimum. This is detailed in Fig. 5(b). The S₁ minimum computed by the the-state-averaged, extended multi-state complete active space second-order perturbation method [SA3-XMS-CAS(2,2)PT2]⁹³ on the same set of ethylene coordinates is found at 90° torsion and 72° pyramidalization and lies 5.46 eV above the ground state minimum [additional information regarding the SA3-XMS-CAS(2,2)PT2 comparison can be found in the supplementary material; these potential energy surfaces can be compared to previous results obtained with quasidegenerate multi-reference perturbation theory²]. These results indicate that hh-TDA describes the potential energy surfaces reasonably well and thus is a promising candidate for applications in efficient nonadiabatic dynamics simulations. This finding complements earlier studies on H₃ and NH₃, which showed that pp-TDA is able to describe conical interactions.⁵³ 

The thymine molecule is a biologically relevant prototype for internal conversion dynamics between nπ^* and ππ^* states. At the Franck–Condon point, highly accurate coupled cluster (CC) computations predict that thymine has nearly degenerate S₁ and S₂ states of nπ^* and ππ^* character, respectively.⁸⁶ Furthermore, accurate CC-based methods predict a reaction coordinate from the Franck–Condon point to the MECI between the S₁ and S₂ states without barriers or the presence of minima on the S₂ state.⁸⁶ At lower levels of theory—particularly those that do not include dynamic electron correlation, such as CASSCF—it is common to find a larger (greater than 1 eV) splitting between the S₁ and S₂ states at the Franck–Condon point as well as a minimum on the S₂ state.^94,95 It remains unclear how these properties of the thymine PES affect its photodynamics; highly efficient electronic structure methods that agree more closely with the accurate CC-based methods would be required to answer these questions. We find hh-TDA to be a promising candidate for use in nonadiabatic dynamics simulations of thymine.

The relative excited state energies computed with hh-TDA with the ωPBE(ω = 0.2 au) functional and EOM-CCSD (see Sec. III) at three different thymine structures are shown in Fig. 6. The vertical excitation energies to the S₁ and S₂ states are underestimated by about 0.5 eV–0.8 eV by hh-TDA-ωPBE(ω = 0.2 a.u.) compared to EOM-CCSD. The latter themselves are, however, higher by about 0.2 eV–0.3 eV compared to the higher-level CC3 data.⁸⁶ The S₁/S₂ splittings are reproduced to within 0.17 eV–0.22 eV by hh-TDA. As a result, hh-TDA does a particularly good job describing the relative energetics of the S₁ and S₂ states at stationary points relevant to the photodynamics of thymine. The stabilization of the S₁ (nπ^*) minimum is correct to 0.04 eV compared to both the Franck–Condon point and S₁/S₂ MECI. Furthermore, no minimum is present on the S₂ PES of thymine between the Franck–Condon point and S₁/S₂ MECI. The relative energy of that intersection with respect to the S₂ energy at the Franck–Condon point is also in good agreement with EOM-CCSD. Remembering that the computational cost of hh-TDA scales formally as O(N⁴) but, in practice, as O(N²) with system size (see Sec. IV D), while EOM-CCSD formally scales as O(N⁶), the agreement between these methods is quite remarkable. What remains to be done in the future studies is to apply hh-TDA directly in a nonadiabatic dynamics simulation of thymine and assess the performance against experimental ultrafast spectroscopic observables.

Malonaldehyde is a molecule that is a prototype for excited-state proton transfer, nπ^*/ππ^* internal conversion dynamics, and photoisomerization. Here, we apply hh-TDA with the ωB97⁷⁶ functional and compare to the results obtained with multireference configuration interaction methods including single and double excitations (MRSDCI, see Fig. 7).⁸⁷ Again, we find good agreement in the absolute excitation energies and in the splitting between nπ^* and ππ^* states at the Franck–Condon point. However, the relative energetics of the S₁ and S₂ states at important stationary points show some potential problems. The minimum on the S₁ (nπ^*) state is a bit too stable relative to both the Franck–Condon point and the S₀/S₁ MECI. This indicates that the population may become spuriously trapped on the S₁ state leading to longer excited-state lifetimes and perhaps increased involvement of triplet states (in view of the El-Sayed rules and the fact that the electronic wavefunction has nπ^* character in the region). In addition, the S₁/S₂ MECI reached by proton transfer is energetically unfavorable. This may indicate that the proton transfer reaction will play a different role in the internal conversion dynamics than MRSDCI would predict. The issues encountered by hh-TDA for malonaldehyde seem to be related to the restriction on the orbital space to include only one π^* orbital. For this molecule, a second π^* orbital may be essential to an accurate description of the photodynamics. Although we expect that malonaldehyde is not an ideal application of the hh-TDA method, nonadiabatic dynamics simulations, in which the molecular symmetry is lifted, would provide a more definitive assessment of the performance of hh-TDA for this molecule.

In Fig. 8, we report the computational time required to solve for the ground and first excited state of the Nile blue chromophore in aqueous solution using the hh-TDA and TDA-TDDFT methods with the Cartesian def2-SVP basis set and the ωPBE functional with ω = 0.8 a.u. The principal conclusions that should be drawn from these timings are that the hh-TDA and TDA-TDDFT methods are of similar computational cost and exhibit similar scaling behavior with the system size, i.e., O(N²) in our AO direct implementation (see below). In TDDFT, the solution of the KS equations defines the ground state, whereas in hh-TDA, both the ground and excited states are obtained as eigenvectors of the TDA response matrix. We find that multiplication of trial vectors against the response matrix is less expensive in the hh-TDA method (because there are no contributions from the Coulomb-type integrals or derivatives of the exchange-correlation potential). In cases where only one or two excited states are required, it should be expected that the cost of TDA-TDDFT and hh-TDA will be quite similar, with the particulars of a given molecule determining which method will be faster. For example, in the present case, TDDFT required extra guess vectors in the Davidson diagonalization in order to avoid solving for higher-lying charge transfer excited states rather than the locally excited state of Nile blue.

Formally, the hh-TDA method appears to scale as O(o⁶), where o represents the number of occupied orbitals. However, since only a few electronic states are actually of interest, the full diagonalization of the response matrix can be avoided, and the scaling of the method would be O(o⁴)—the cost of the multiplication of a trial vector by the response matrix. However, this ignores the cost of forming the response matrix in the molecular orbital basis. In practice, our implementation of hh-TDA is fully AO integral direct and scales, formally, as O(N⁴), where N is the number of AO basis functions. The advantage of the AO formulation of hh-TDA is that sparsity can be exploited and the contractions between trial vectors and the response matrix scale more like O(N²). This is clearly demonstrated in Fig. 8, where the apparent quadratic scaling of the method allows us to apply it to systems with more than 2100 atoms in a polarized double-ζ basis set (roughly 18 000 basis functions). It should be noted that the underlying SCF procedure costs approximately the same as the solution of the eigenvalue problem (for two states) and the SCF procedure scales slightly worse with the system size.

Following the remarks above on the basis set requirements for hh-TDA due to the (N+2)-electron reference, all calculations in this paper used double- and triple-zeta basis sets without diffuse functions. In order to explore the basis set dependence more thoroughly, we collect the lowest excitation energies computed with hh-TDA-ωB97X and different basis sets STO-3G, def2-SVP, def2-SVPD, def2-TZVP, and def2-QZVP (the latter two without f and g functions) for a few molecular systems (acetamide, p-benzoquinone, butadiene, ethylene, formaldehyde, norbornadiene, uracil, and water) in Table S11 of the supplementary material. As could be expected, these results show that excitation energies obtained with the STO-3G basis set are widely disparate from larger basis sets. However, with the exception of acetamide, the excitation energies vary mostly by less than 0.2 eV over the other choices of basis sets without diffuse functions. This indicates that the presence of unbound orbitals need not be a significant problem (as long as diffuse basis sets are avoided, see results for def2-SVPD).

Acetamide was among the problematic systems mentioned in Sec. IV A, where DFAs with insufficient exact exchange populated an incorrect orbital in the (N+2)-electron reference. Similar occupation problems are observed for the def2 basis sets considered here, leading to strong basis set dependence of the excitation energy. However, this seems to be distinct from the unbound orbital problem and is instead related to the limited active space in hh-TDA.

Based on this assessment and the promising results obtained in this work, we therefore recommend the use of hh-TDA in combination with asymptotically correct range-separated hybrid DFAs and double- and triple-zeta basis sets without diffuse basis functions. The latter can be added if the excited states to be probed by hh-TDA are actually of Rydberg character (cf. results for water). In the future, the unbound orbital issue could be addressed by using an N-electron reference in the orbital generation. Such an approach was already applied in one of the early works on pp-RPA.⁴⁹ We plan to investigate this issue in order to facilitate the routine use of hh-TDA in geometry optimizations and nonadiabatic dynamics simulations.

We have shown that the hole–hole Tamm–Dancoff approximation (hh-TDA) to the particle–particle random-phase approximation (pp-RPA) represents an efficient DFT-based electronic structure scheme for the electronic ground and low-lying excited states. The method shares similarities with the previously presented particle–particle (pp-TDA) approach. We find some advantages of the hh-TDA scheme, particularly in its ability to describe different low-lying excitation types (nπ^*, σπ^*, and ππ^*). Starting from an (N+2)-electron reference determinant, the ground and the low-lying excited states of the N-electron target system are obtained by a double annihilation of two electrons. Since the ground and excited states are obtained from the same eigenvalue problem, static correlation cases can be handled. Dynamic correlation is included in the orbitals by virtue of the density functional. We have also introduced an alternative choice for the response kernel that differs from the functional-independent Hartree–Fock-type kernel employed in the previous works on pp-RPA and pp-TDA.⁵⁴ Our choice is functional-dependent and in line with the response kernels appearing in particle–hole linear response theory.

To avoid SCF convergence problems of the (N+2)-electron reference, the method should be combined with range-separated density functionals that have 100% asymptotic Fock exchange. We recommend avoidance of diffuse atom-centered basis sets since these can lead to continuum-like orbitals in the (N+2)-electron reference. We then do not observe any SCF convergence problems and find that the hh-TDA can describe low-lying excited states reasonably well at a computational complexity comparable to a configuration interaction singles (CIS) calculation. As with any density functional approximation-based method, case-specific functional assessment is recommended. The range-separation parameter and amount of short-range Fock exchange can be viewed as parameters to optimize the method performance for a specific molecule (as is often done by variation of the active space in CAS methods).

On the other hand, the main limitation of the hh-TDA is the restricted “virtual” orbital space, which formally consists exclusively of the lowest unoccupied orbital (with respect to the N-electron system). This precludes its application to systems with degenerate LUMOs, such as benzene. In spite of this shortcoming, the hh-TDA method is well-suited for the calculation of ground and low-lying excited states for many organic molecules. It offers an inexpensive alternative to the existing and widely used CASSCF methods, with the advantage that dynamic correlation effects are included at the orbital level by means of the used density functional.

The method is implemented in the GPU-accelerated electronic structure code TeraChem. Future directions will include the computation of nonadiabatic couplings and testing the method in nonadiabatic dynamics simulations. In this context, it will be important to address the problem of potentially occupying unbound orbitals in the (N+2)-electron SCF calculation. Generating orbitals for an N-electron reference with fractional orbital occupation and using them in an (N+2)-electron-type hh-TDA procedure could remove this issue entirely, while guaranteeing stability in nonadiabatic dynamics simulations. This will be explored in the future work.

See the supplementary material for a zip-file containing hh-TDA optimized geometries of thymine (Fig. 6) and malonaldehyde (Fig. 7) and a pdf file containing detailed results on the DFA assessment (Sec. IV A) and the pp/hh-TDA benchmarks (Sec. IV B) and a potential curve of ethylene (at 90° torsion).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was supported by the AMOS program of the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Chemical Sciences, and Biosciences Division and by the Department of Energy, Laboratory Directed Research and Development program at the SLAC National Accelerator Laboratory under Contract No. DE-AC02-76SF00515. T.J.M. is a co-founder of PetaChem, LLC. C.B. thanks the German National Academy of Sciences Leopoldina for support through the Leopoldina Fellowship Program (Project No. LPDS 2018-09). The authors thank Deniz Tuna for helpful discussions and for early testing of the method.