In this article, we develop systematically second random phase approximations (RPA) and Tamm-Dancoff approximations (TDA) of particle-hole and particle-particle channels for calculating molecular excitation energies. The second particle-hole RPA/TDA can capture double excitations missed by the particle-hole RPA/TDA and time-dependent density-functional theory (TDDFT), while the second particle-particle RPA/TDA recovers non-highest-occupied-molecular-orbital excitations missed by the particle-particle RPA/TDA. With proper orbital restrictions, these restricted second RPAs and TDAs have a formal scaling of only O(N4). The restricted versions of second RPAs and TDAs are tested with various small molecules to show some positive results. Data suggest that the restricted second particle-hole TDA (r2ph-TDA) has the best overall performance with a correlation coefficient similar to TDDFT, but with a larger negative bias. The negative bias of the r2ph-TDA may be induced by the unaccounted ground state correlation energy to be investigated further. Overall, the r2ph-TDA is recommended to study systems with both single and some low-lying double excitations with a moderate accuracy. Some expressions on excited state property evaluations, such as

$\langle \hat{S}^{2}\rangle$
Ŝ2 are also developed and tested.

The particle-hole random phase approximation (ph-RPA)1,2 has been a convenient method to study particle-hole excitations and correlation energies for nuclei,3–8 molecules9–12 and solid.13–17 The ph-RPA can be viewed as a response in time-dependent Hartree theory where the exchange correlation contribution is omitted in time-dependent density-functional theory (TDDFT),18 as the correlation energies of all ring diagrams,7 or alternatively as a ring approximation in coupled-cluster doubles.19 Viewed as a correlation energy functional from the adiabatic connection in density-functional theory (DFT),10 there is a renaissance of the ph-RPA in quantum chemistry community due to its good description of van der Waals interactions10 and the correct dissociation limit of H2.20 These features have incentivized developments of fast algorithms for ph-RPA correlation energies.12,17 Nonetheless, the ph-RPA has formidable fractional charge errors which prohibit the applications of the ph-RPA to many systems.21 

On the other hand, the pp-RPA,22–27 also known as Brueckner's theory,28–30 has been a textbook method in nuclear physics to study pairing vibrations.7,8 The pp-RPA can be interpreted as time-dependent Hartree-Fock-Bogoliubov approximation,7 as linear-response time-dependent density-functional theory with pairing field at the zero pairing field limit,31 as the adiabatic connection of the pairing matrix fluctuations,32,33 as sum of all ladder diagrams,7 or as a ladder approximation in coupled-cluster doubles.34–36 Recent applications of the pp-RPA in molecular systems reveal that the pp-RPA satisfies the flat-plane condition32,37 and has better thermochemistry behavior than the ph-RPA.38 The pp-RPA can also capture double excitations from an (N − 2)-electron reference, which is impossible for adiabatic linear-response TDDFT.39,40 However, due to the limitation of the (N − 2)-electron reference construction, single excitations from non-highest occupied molecular orbitals (non-HOMO) are absent in the pp-RPA as normally applied—although the pp-RPA with a non-ground state reference can in principle lead to such excitations. In this article, we develop the restricted second random-phase approximation and the restricted second particle-particle random phase approximation that can calculate full single excitation spectrum and some double excitations.

The second particle-hole random phase approximation (2ph-RPA)6,41–44 is a natural extension of the ph-RPA or time-dependent Hartree-Fock (TDHF) in the equation-of-motion (EOM)26 framework, where the excitation operators include both one-particle-one-hole (1p1h) and two-particle-two-hole (2p2h) excitations. The 2ph-RPA has been applied to nuclear physics42,43 and metal clusters44 to study double excitations, but not in chemistry in general. McKoy and coworkers’ higher RPA45,46 is similar albeit involves a very sophisticated ground state treatment. 2p2h excitations have also been used to improve 1-particle excitations.47–51 

Parallel to the 2ph-RPA, we devise here the second particle-particle random phase approximation (2pp-RPA) that supplements the pp-RPA excitation operators with three-particle-one-hole (3p1h) and one-particle-three-hole (1p3h) operators. From an (N − 2)-electron reference, the 2pp-RPA has all critical double excitations in the pp-RPA, plus all single excitations. The philosophy behind the 2pp-RPA is very similar to including 3p1h/1p3h operators and above in double-ionization-potential/double-electron-affinity equation-of-motion coupled-cluster (DIP/DEA-EOM-CC) methods.52,53 The computational scaling of second RPAs (the 2ph-RPA and the 2pp-RPA) can be reduced by placing some rational restrictions on the excitation operators, leading to a formal scaling of O(N4), the same as TDDFT and the pp-RPA. Some molecular tests show that these restricted second RPAs can capture the most important low-lying excitations while keeping the computational complexity manageable. Preliminary results54 show that ph series RPAs and TDAs are probably size extensive while pp series RPAs and TDAs are not. Although these second RPAs can be used to study correlation energies as well, the focus of this article is on excitation energies.

This article is organized as follows. The theories of the 2ph-RPA and the 2pp-RPA and their restrictions are described in Sec. II. Section III presents the implementation and calculation details for these methods. Results are shown in Sec. IV with discussions. Finally, Sec. VI concludes this article.

Second RPA theories are expressed in Rowe's EOM formalism.26,55 For an electronic Hamiltonian

$\hat{H}$
Ĥ⁠, we have its eigenvalues and eigenvectors

(1)

Note that

$\hat{H}$
Ĥ is expressed in second quantization and |M⟩'s do not have to have the same number of electrons. Suppose we have an initial state |0⟩ and want to study its excitation spectrum to some final state |F⟩ and F ≠ 0. Defining an excitation operator

(2)

we have

(3)

and

(4)

The EOM equation for the transition energy is

(5)

where ωF = EFE0, and the double commutator is

(6)

for Bosonic operators, or

(7)

for Fermionic operators.6,55

Equation (5) alone is of little use as neither |0⟩ nor

$\hat{O}_{F}^{\dagger }$
ÔF is an easy task for a general system. By approximating the initial state |0⟩ and the excitation operator
$\hat{O}_{F}^{\dagger }$
ÔF
, we can obtain an eigenvalue equation derived from the stationary condition without any constraint of Eq. (5). Expanding
$\hat{O}_{F}^{\dagger }$
ÔF
as a linear combination of some operators,

(8)

and approximating |0⟩ as a model state |MS⟩, ωF is then a function of of the vectors of XIF and YIF. Setting the variation of ωF zero leads to the generalized eigenvalue equation

(9)

where

(10)
(11)
(12)

and

(13)

The 2ph-RPA is the EOM with |0⟩ approximated by the Hartree-Fock (HF) wave function and

$\hat{Q}_{I}^{\dagger }$
Q̂I approximated by all single and double particle-hole excitation operators, i.e., operators of {ai} and {aibj}. In this article, a, b, c, ⋅⋅⋅ represent unoccupied spin orbitals, and i, j, k, ⋅⋅⋅ represent occupied spin orbitals, while p, q, r, s, ⋅⋅⋅ represent general spin orbitals. {⋅⋅⋅} indicates the operator is normal ordered with respect to the Fermi sea.56 The resulting eigenvalue equation is42–44,57

(14)

where subscripts S and D denote the single and double excitation block, respectively, and I is an identity matrix. The matrix elements are

(15)
(16)
(17)

and

(18)

where Fpq is the Fock matrix element, ⟨pq||rs⟩ is an antisymmetrized two-electron integral

(19)

with

(20)

and U(pq) is an operator that antisymmetrizes the term with respect to p and q,

(21)

Note that only the single-single block of B matrix is nonzero. For a HF reference, Fpq = δpqεp where εp is the molecular orbital eigenvalue.

It is interesting to compare Eq. (14) with the configuration interaction singles and doubles (CISD) equation58 

(22)

where

(23)

and ASS, ASD, and ADD are defined the same as in the 2ph-RPA. Therefore, the CISD matrix and the 2ph-RPA matrix contain exactly the same amount of information, i.e., we can build the CISD equation from the 2ph-RPA equations, and vice versa. However, CISD includes the double excitation configurations to approximate both the ground and excited states, while the 2ph-RPA only approximates the excited states with no direct reference to the ground state. The rearrangement of matrix elements makes CISD and the 2ph-RPA dramatically different.

It is also noted that the matrix in the TDHF equation

(24)

is a submatrix of the 2ph-RPA matrix. It is well-known that TDHF suffers from HF instability issues for many systems.3,59 Therefore, any instability in TDHF will lead to instability in the 2ph-RPA (as well as the restricted 2ph-RPA introduced later). The Tamm-Dancoff approximation (TDA) is the approximation to set the B matrix zero, which could ameliorate some instability issues and improve the result. The TDA for TDHF is configuration interaction singles (CIS). The TDA for the 2ph-TDA, dubbed 2ph-TDA here, is then

(25)

Note that the 2ph-TDA matrix is a submatrix of the CISD matrix.

Parallel to the 2ph-RPA, by complementing the 3p1h {abci} and 1p3h {aijk} operators in the pp-RPA derivation, we obtain here the 2pp-RPA equation

(26)

with

(27)
(28)
(29)
(30)
(31)
(32)

and

(33)

All the expressions are derived using the Wick's theorem contraction (see, for example, Ref. 56). Again, only the 2p-2h block of the B matrix is nonzero. Note U(def) is an antisymmetrized operator with three arguments,

Similarly, we can obtain the 2pp-TDA by setting B zero,

(34)

In the EOM formalism, since there is no explicit excited state wave function, properties of excited states such as density matrices and S2 expectation values are not directly available. The motivation for excited state property calculations originates from distinguishing spin states in the results of an unrestricted calculation. There were a few discussions on this topic in the literature. The original proposal by Rowe is84 

(35)

Alternatively, we can use the expression

(36)

with A(W) and B(W) defined as in Eqs. (10) and (11) but with

$\hat{H}$
Ĥ replaced by
$\hat{W}$
Ŵ
. Yeager and co-workers60–62 proposed some other formulas of
$\langle F|\hat{W}|F\rangle$
F|Ŵ|F
. However, due to the complexity of their expressions, Eq. (35) is still the most used expression in practice. Equation (35) does not always produce reasonably values. For example, triplet excited states from a closed-shell singlet reference in TDDFT or TDHF has an S2 expectation value of Ref. 63,

(37)

which is larger than 2 unless YF is zero. On the other hand, Casida explicitly assigned the excited state wave function |F⟩ in TDDFT,18 

(38)

where ZF is the orthonormal eigenvector of a transformed TDDFT equation,

(39)

In TDDFT, ZF is related to XF and YF,

(40)

Casida's explicit construction of the wave function enables the calculation of excited state expectation values. In this way, Casida's expectation value is equivalent to

(41)

with Δεia, jb = δijδaba − εi). Similarly, Ipatov et al.63 proposed another expectation value formula related to Casida's ansatz,

(42)

Equations (41) and (42) are preferable for delivering exact

$\langle F|\hat{S}^{2}|F\rangle$
F|Ŝ2|F for closed-shell triplet excitations. In this subsection, we will discuss some other possible expressions for excited state properties.

Equations (41) and (42) actually resemble the expectation value expressions in the TDA,

(43)

Therefore, we deem that B(W) and YF may not be useful in expectation evaluation, and propose to use Eq. (43) for general EOM solutions even though YF may not be zero. Note that for hole-hole excitations in pp-RPA/2pp-RPA solutions, Eq. (43) should be slightly modified,

(44)

Note also that in EOM the normalization requires that

Therefore the normalization factors in Eqs. (43) and (44) are required. Alternatively, Eq. (37) delivers an intuition that reversing the sign of the

$\mathbf {Y}_{F}^{\dagger }\mathbf {Y}_{F}$
YFYF term may be a good way to calculate expectation values. Thus we have a revised expression

(45)

where G(W) = A(W)* for particle-hole excitations and G(W) = C(W) for particle-particle or hole-hole excitations. Note that for

$\hat{W}=\hat{H}$
Ŵ=Ĥ⁠, Eq. (45) could produce excitation energies quite different from the eigenvalues from the original EOM equation. The performance of expressions above will be discussed in Sec. III.

Before we conclude this subsection, we have a brief description of obtaining the matrix elements of A(W), B(W), and C(W) for

$\hat{S}^{2}$
Ŝ2⁠. The
$\hat{S}^{2}$
Ŝ2
operator can be expressed as

(46)

Since both the reference and excited states can be

$\hat{S}_{z}$
Ŝz eigenvectors from construction, the operator we need to study is

(47)

Employing normal order techniques,

$\hat{\Sigma }$
Σ̂ can be partitioned into a constant, a one-body operator, and a two-body operator in normal order,

(48)

Therefore, the normal ordered matrix elements of

$\hat{\Sigma }$
Σ̂ are

(49)

and

(50)

All these matrix elements can be expressed via the spatial overlap matrix between α and β orbitals. Then the resulting matrix elements of A(Σ), B(Σ), and C(Σ) are the same as those in A, B, and C except that Fpq is replaced by

$(\Sigma _{N}^{1})_{pq}$
(ΣN1)pq and ⟨pq||sr⟩ is replaced by
$(\Sigma _{N}^{2})_{pq,sr}$
(ΣN2)pq,sr
. Refer to the supplementary material for detailed expressions.54 

Those eigenvalue problems in second RPAs and their TDAs stated above are still too time-consuming, with a scaling of O(N6) even with Davidson's algorithm.64 Our goal for utilizing these extended models is to capture all single excitations and the most important—low-lying—double excitations that remedies the lack of double excitations in the ph-RPA/ph-TDA or the missing of non-HOMO single excitations in the pp-RPA/pp-TDA. Therefore, for the 2ph-RPA/2ph-TDA, we will restrict the double excitation tensors to include only double excitations from HOMOs. The resulting models are thus named the restricted second particle-hole random phase approximation (r2ph-RPA) and the restricted second Tamm-Dancoff approximation (r2ph-TDA). Depending on the degeneracy of the frontier orbitals, HOMOs may include multiple orbitals. The resulting excitation operator for the r2ph-TDA is then

(51)

and similar for the r2ph-RPA.

On the other hand, for the 2pp-TDA, we limit the 3p1h excitation tensors to those that place two electrons back to HOMOs of the N-electron state (the LUMOs of the (N − 2)-electron state). Then the 2pp-TDA operator is

(52)

Note that the HOMOs in Eqs. (51) and (52) are both HOMOs for N-electron states. If the N-electron state has degenerate HOMOs (e.g., HOMOs of benzene have twofold spatial degenerate), the self-consistent (N − 2)-electron state is likely to break the spatial symmetry which could cause problems for excitation energy calculations. Therefore, HOMOs in the 2pp-TDA usually contains one spatial orbital (two spin orbitals). Due to the asymmetry of 3p1h operators and 1p3h operators, the restricting 3p1h operators while keeping 1p3h operators unchanged seems unbalanced. Therefore, only the restricted 2pp-TDA but not the restricted 2pp-RPA be discussed. Note that the restricted orbitals could potentially be non-HOMOs due to the state of interest; however, in this paper we only focus on restrictions on HOMOs as we want to study the low lying double excitations. This construction also implies that all the restricted methods are not orbital invariant under unitary transformation among occupied orbitals since we can only use canonical orbitals to determine HOMOs.

With these orbital restrictions, the scaling of the r2ph-RPA, the r2ph-TDA, and the r2pp-TDA can be reduced to O(N4) with Davidson's diagonalization.64 

All the methods described above are implemented in QM4D,65 except for the unbalanced restricted 2pp-RPA. For testing purpose, all excitations have ΔSz = 0, and only closed-shell systems are tested. At the current stage, we only implemented direct diagonalization which scales as O(N6) for the restricted second RPAs. Full CI (FCI), multireference configuration interaction (MRCI), and EOM-CCSD calculations are done in GAMESS (US),66MOLPRO,67,68 and Gaussian 09,69 respectively. By default, neutral excitation energies for electron number conserving EOMs are calculated at N-electron reference, while for pp-methods excitations are calculated at (N − 2)-electron reference, with the excitation energies expressed as the difference of the double electron affinity of an N-electron state and that of the lowest N-electron state. All geometries used in this study are adopted from G2/97 set of MP2 neutrals70 except that BH molecule uses a bond length of 1.232 Å adopted from Ref. 51. Excitation energies studied in this article are all vertical excitation energies. Cartesian basis sets are used when possible. Refer to the supplementary material for geometries and details of MRCI calculations.54 

By default, second RPAs use HF references except for TDDFT and pp-RPA/pp-TDA@B3LYP71,72 calculations. Note that @B3LYP indicates that the normal ordered Hamiltonian used in the EOM is

(53)

with

$\epsilon _{r}^{\text{B3LYP}}$
εrB3LYP the B3LYP molecular orbital eigenvalue, rather than the true Hamiltonian which uses the nondiagonal Fock matrix.39 B3LYP reference was tried for some second RPA and second TDA calculations with many unphysical low-lying double excitations, probably due to the ad hoc choice of the Hamiltonian of Eq. (53). Therefore, B3LYP reference is only used in the pp-RPA/pp-TDA@B3LYP but not in any second RPAs/TDAs. In contrast, the use of B3LYP or other DFT reference is well justified for pp-RPA calculations, because one can view it as a linear response calculation for time-dependent DFT in an external paring field.31 Another way to understand this is from the fact that the pairing perturbation to a molecule can be described by time-dependent Kohn-Sham Bogoliubov theory, and the pp-RPA/pp-TDA@B3LYP equation is just the linear response of the time-dependent Kohn-Sham Bogoliubov theory at the limit of zero pairing field.31 No such connection exists for the 3p1h or 1p3h perturbation.

For comparisons, Second-Order Polarizability Propagator Approximation (SOPPA) and Random Phase Approximation with Doubles corrections (RPA(D)) results from Dalton73,74 are also added, as these two methods contains corrections from double excitations.

We first study the behavior of various

$\langle \hat{S}^{2}\rangle$
Ŝ2 formulas for CO with 6-31G basis set as an example for the lowest triplet excited states a3Σ+ (π → π*) (see Table I). It is noted that results from Rowe's formula deviate from the exact value (2) very much, especially for TDHF and the 2ph-RPA. Both TDA and Rev formulas deliver the exact value for all methods. They are indistinguishable for all examples in this article. The
$\langle \hat{S}^{2}\rangle$
Ŝ2
assignments for the rest of the calculations in this article use the TDA formula in Eq. (43) by default.

Table I.

$\hat{S}^{2}$
Ŝ2 expectation values for different formulas. The state to study is the a3Σ+(π → π*) state of CO, with the basis set 6-31G. Rowe, TDA, and Rev refer to Eqs. (35), (43), and (45), respectively.

FormulaRoweTDARev
TDHF 2.133 2.000 2.000 
2ph-RPA 2.215 2.000 2.000 
pp-RPA 2.003 2.000 2.000 
2pp-RPA 2.004 2.000 2.000 
FormulaRoweTDARev
TDHF 2.133 2.000 2.000 
2ph-RPA 2.215 2.000 2.000 
pp-RPA 2.003 2.000 2.000 
2pp-RPA 2.004 2.000 2.000 

H2O is used as a test case for all methods with a relatively small basis set 6-31G. Since the basis set is small enough, all methods mentioned above can be tested. All data are listed in Table II, with FCI calculations in the same basis set as a reference. Numbers in Table II for HF-based methods are not supposed to compare to experimental values because of the insufficient size of the basis set. The molecule lies in the yz plane while the principal axis is the z axis. The ground state configuration of the frontier orbitals in the basis set 6-31G is 4(a1)25(b1)26(a1)07(b2)0. EOM-CCSD is very accurate for all the single excitations listed in the table, closely followed by the r2ph-RPA, the r2ph-TDA, SOPPA, and RPA(D). CIS, TDHF, and TDDFT with the B3LYP functional are most practiced methods to calculate single excitations, with the average errors 0.5–1.0 eV shown in Table II. As a general trend, the 2ph-RPA and the 2ph-TDA systematically pull down excitations from TDHF and CIS. The 2ph-RPA and the 2ph-TDA correct CIS and TDHF excitation energies in the right direction. This is probably caused by the unbalanced treatment between the excited and ground states. Surprisingly, with hole restriction, the r2ph-RPA and the r2ph-TDA produce much better results, with errors below 0.5 eV. The pp-RPA and the pp-TDA (with HF references by default) substantially underestimate the excitation energies by more than 4 eV. The pp-RPA@B3LYP and the pp-TDA@B3LYP, as DFT-based methods, are different from the FCI results for 1.5 eV, but because of the small basis set for FCI, this difference does not represent the real error. A more fair comparison will be presented in the following paragraph and Table III. At this time, the 2pp-RPA and the 2pp-TDA produce better results while r2pp-TDA overestimate most excitations. The defect of the pp-RPA and the pp-TDA that single excitations from a non-HOMO level are absent is conspicuous in 3A1 and 1A1 states, which is ameliorated in the 2pp-RPA, the 2pp-TDA and the r2pp-TDA. Note that for all states studied in Table II, no double excitations are involved.

Table II.

Vertical excitation energies of H2O with different methods using a small basis set 6-31G. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

   EOM     2ph-2ph-r2ph-r2ph-  pp-RPApp-TDA2pp-2pp-r2pp-
SymbolExcitationsFCI-CCSDCISTDHFTDDFTaSOPPARPA(D)bRPATDARPATDApp-RPApp-TDA@B3LYPc@B3LYPcRPATDATDA
3B1 5(b1) → 6(a17.63 7.53 8.33 8.20 6.92 7.27 NA 5.83 5.95 7.68 7.80 3.09 2.98 5.84 5.65 7.58 7.46 11.25 
1B1 5(b1) → 6(a18.37 8.27 9.30 9.24 7.72 8.06 8.13 6.61 6.67 8.75 8.80 3.60 3.49 6.47 6.28 8.34 8.22 11.79 
3A1 4(a1) → 6(a19.85 9.77 10.12 9.80 8.66 9.47 NA 7.89 8.15 9.80 10.11 NA NA NA NA 12.65 12.53 12.56 
3A2 5(b1) → 7(b210.07 10.02 10.55 10.42 9.24 9.80 NA 8.27 8.38 9.98 10.10 5.49 5.38 8.25 8.06 9.98 9.86 13.67 
1A2 5(b1) → 7(b210.58 10.52 11.20 11.13 9.83 10.33 10.43 8.79 8.86 10.70 10.77 5.76 5.65 8.61 8.43 10.49 10.37 13.95 
1A1 4(a1) → 6(a110.93 10.86 11.77 11.69 9.88 10.56 10.62 9.11 9.19 11.68 11.77 NA NA NA NA 13.65 13.53 13.53 
MSEd … −0.08 0.64 0.51 −0.87 −0.32 −0.23 −1.82 −1.71 0.19 0.32 −4.68 −4.79 −1.87 −2.06 0.87 0.75 3.22 
MAEe … 0.08 0.64 0.53 0.86f 0.32 0.23 1.82 1.71 0.24 0.32 4.68 4.79 1.87f 2.06f 0.96 1.00 3.22 
   EOM     2ph-2ph-r2ph-r2ph-  pp-RPApp-TDA2pp-2pp-r2pp-
SymbolExcitationsFCI-CCSDCISTDHFTDDFTaSOPPARPA(D)bRPATDARPATDApp-RPApp-TDA@B3LYPc@B3LYPcRPATDATDA
3B1 5(b1) → 6(a17.63 7.53 8.33 8.20 6.92 7.27 NA 5.83 5.95 7.68 7.80 3.09 2.98 5.84 5.65 7.58 7.46 11.25 
1B1 5(b1) → 6(a18.37 8.27 9.30 9.24 7.72 8.06 8.13 6.61 6.67 8.75 8.80 3.60 3.49 6.47 6.28 8.34 8.22 11.79 
3A1 4(a1) → 6(a19.85 9.77 10.12 9.80 8.66 9.47 NA 7.89 8.15 9.80 10.11 NA NA NA NA 12.65 12.53 12.56 
3A2 5(b1) → 7(b210.07 10.02 10.55 10.42 9.24 9.80 NA 8.27 8.38 9.98 10.10 5.49 5.38 8.25 8.06 9.98 9.86 13.67 
1A2 5(b1) → 7(b210.58 10.52 11.20 11.13 9.83 10.33 10.43 8.79 8.86 10.70 10.77 5.76 5.65 8.61 8.43 10.49 10.37 13.95 
1A1 4(a1) → 6(a110.93 10.86 11.77 11.69 9.88 10.56 10.62 9.11 9.19 11.68 11.77 NA NA NA NA 13.65 13.53 13.53 
MSEd … −0.08 0.64 0.51 −0.87 −0.32 −0.23 −1.82 −1.71 0.19 0.32 −4.68 −4.79 −1.87 −2.06 0.87 0.75 3.22 
MAEe … 0.08 0.64 0.53 0.86f 0.32 0.23 1.82 1.71 0.24 0.32 4.68 4.79 1.87f 2.06f 0.96 1.00 3.22 
a

TDDFT with B3LYP functional.

b

Due to the limitation of the program Dalton, only singlet data are listed.

c

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

d

Mean signed error. Errors are with respect to FCI values. NA values are excluded.

e

Mean absolute error. Errors are with respect to FCI values. NA values are excluded.

f

The errors for DFT-based methods are not fairly justified errors because of the small basis set and thus inaccurate FCI results. More fair comparison is presented in Table III with a larger basis set.

Table III.

H2O excitation spectra with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 except for MRCI in MOLPRO,67,68 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

    EOM-         pp-RPApp-TDA 
SymbolExcitationsExpt.aMRCICCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3B1 5(b1) → 6(a17.1 7.18 7.05 7.91 7.78 6.48 6.93 7.05 6.59 NA 3.08 3.01 6.34 6.19 10.27 
1B1 5(b1) → 6(a17.5 7.57 7.46 8.60 8.55 6.86 7.57 7.61 6.94 6.96 3.30 3.23 6.76 6.61 10.51 
3A2 5(b1) → 7(b29.0 9.18 9.05 9.91 9.78 8.15 8.94 9.05 8.49 NA 4.82 4.75 8.50 8.35 12.03 
1A2 5(b1) → 7(b29.1 9.36 9.22 10.27 10.22 8.28 9.30 9.34 8.63 8.75 4.88 4.81 8.71 8.57 12.09 
3A1 4(a1) → 6(a19.3 9.48 9.37 10.00 9.75 8.55 9.73 9.98 8.89 NA NA NA NA NA 12.07 
1A1 4(a1) → 6(a19.7 9.95 9.86 10.91 10.88 9.05 10.62 10.65 9.31 9.33 NA NA NA NA 12.43 
MSEe … … −0.12 0.82 0.71 −0.89 −0.06 0.16 −0.48 −0.42 −4.30 −4.37 −0.75 −0.89 2.78 
MAEf … … 0.12 0.82 0.71 0.89 0.25 0.26 0.48 0.42 4.30 4.37 0.75 0.89 2.78 
    EOM-         pp-RPApp-TDA 
SymbolExcitationsExpt.aMRCICCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3B1 5(b1) → 6(a17.1 7.18 7.05 7.91 7.78 6.48 6.93 7.05 6.59 NA 3.08 3.01 6.34 6.19 10.27 
1B1 5(b1) → 6(a17.5 7.57 7.46 8.60 8.55 6.86 7.57 7.61 6.94 6.96 3.30 3.23 6.76 6.61 10.51 
3A2 5(b1) → 7(b29.0 9.18 9.05 9.91 9.78 8.15 8.94 9.05 8.49 NA 4.82 4.75 8.50 8.35 12.03 
1A2 5(b1) → 7(b29.1 9.36 9.22 10.27 10.22 8.28 9.30 9.34 8.63 8.75 4.88 4.81 8.71 8.57 12.09 
3A1 4(a1) → 6(a19.3 9.48 9.37 10.00 9.75 8.55 9.73 9.98 8.89 NA NA NA NA NA 12.07 
1A1 4(a1) → 6(a19.7 9.95 9.86 10.91 10.88 9.05 10.62 10.65 9.31 9.33 NA NA NA NA 12.43 
MSEe … … −0.12 0.82 0.71 −0.89 −0.06 0.16 −0.48 −0.42 −4.30 −4.37 −0.75 −0.89 2.78 
MAEf … … 0.12 0.82 0.71 0.89 0.25 0.26 0.48 0.42 4.30 4.37 0.75 0.89 2.78 
a

Experimental data collected from Ref. 75.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Errors are with respect to the MRCI values. NA values are excluded.

f

Mean absolute error. Errors are with respect to the MRCI values. NA values are excluded.

We also test H2O with a larger basis set, aug-cc-pVTZspd, i.e., the aug-cc-pVTZ basis set with no f functions. The frontier orbital configuration of the ground state in this basis set remains the same: 4(a1)25(b1)26(a1)07(b2)0. Results are shown in Table III. Experimental values are imprecise for H2O75, therefore we use MRCI results as reference because MRCI and EOM-CCSD results are very close to each other and similar to experimental data. Again, CIS and TDHF generally overestimate excitations by about 0.8 eV, while TDDFT underestimates excitations by about 0.9 eV. The full 2ph-RPA and the full 2ph-TDA cannot be calculated because of memory limit. We see that both the r2ph-RPA and the r2ph-TDA perform reasonably well with the MAE (mean absolute error) of about 0.25 eV. The pp-RPA and the pp-TDA substantially underestimate excitations with the HF reference. In contrast, using B3LYP reference reduces the error to less than 1 eV. The r2pp-TDA is better than the pp-RPA and the pp-TDA, but the results are significantly worse than those of the r2ph-RPA and the r2ph-TDA. Again, note that the pp-RPA and the pp-TDA fail to capture non-HOMO excitations using the ground state (N − 2)-electron reference. Again, the accuracy of SOPPA and RPA(D) results are close to the those of r2ph-RPA and r2ph-TDA.

The systematic overestimation of r2pp-TDA results and be understood using configuration interaction (CI) reasoning, since these TDAs are actually special CI methods. Due to the two-body nature of the electronic Hamiltonian, only configurations differing no more than two electrons can directly interact with a specific configuration. In the r2pp-TDA, the N −electron ground state has correlations from all single excitations and many double excitations, while N −electron excited states has significantly less single (de)excitations that couple with the zeroth order configurations. Therefore, the ground state is more correlated than the excited states and the excitation energies are usually overestimated.

Be is a small atom with low-lying double excitations. Table IV shows excitation energies of various methods with the basis set aug-cc-pVTZspd. Note that the HF reference of the Be atom suffers from instability, and those imaginary eigenvalues in TDHF, 2ph-RPA and r2ph-RPA are marked with a suffix “i”. This presentation is different from some literature that treats imaginary eigenvalues as negative rather than imaginary values.51,63 In this case, pp-series methods based on the HF reference are very accurate, since there are only two valence electrons in the system. The pp-RPA@B3LYP and the pp-TDA@B3LYP perform quite well for low-lying excitations, but their errors grow rapidly for higher excitations. Such behavior for B3LYP reference has been observed in a previous study.39 CIS, TDHF, and TDDFT cannot capture double excitations as expected. Among methods based on N-electron reference, the r2ph-TDA is the best as it is free from instability issues and capable of describing double excitations while keeping the error relatively low. In fact, r2ph-TDA systematically underestimates the excitation energies, showing that there seems to be a missing ground state correlation energy. SOPPA and RPA(D) have difficulty capturing the double excitation dominated state 1D.

Table IV.

Be excitation spectra with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 FCI in GAMESS (US),66 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

    EOM-   2ph-2ph-r2ph-r2ph-    pp-RPApp-TDAr2pp-
SymbolConfigurationExpt.aFCICCSDCISTDHFbTDDFTcRPAbTDARPAbTDASOPPARPA(D)dpp-RPApp-TDA@B3LYPe@B3LYPeTDA
3P 2s12p1 2.73 2.72 2.73 1.70 −0.96i 2.10 −1.27i 1.45 −1.22i 1.49 1.81 NA 2.73 2.73 2.82 2.81 2.73 
1P 2s12p1 5.28 5.33 5.36 5.10 4.84 4.90 3.74 4.07 3.77 4.10 4.86 4.87 5.36 5.36 6.15 6.15 5.34 
3S 2s13s1 6.46 6.43 6.44 5.53 5.49 5.69 5.12 5.16 5.15 5.19 6.47 NA 6.44 6.44 8.40 8.40 6.44 
1S 2s13s1 6.78 6.77 6.77 6.17 6.16 5.97 5.45 5.47 5.47 5.49 6.01 6.48 6.77 6.77 8.68 8.68 6.77 
1D 2p2 7.05 7.16 7.18 NA NA NA 5.92 5.92 5.94 5.94 8.83 NA 7.18 7.18 7.98 7.97 7.18 
3P 2s13p1 7.30 7.42 7.42 6.56 6.55 6.08 6.15 6.16 6.17 6.18 7.02 NA 7.43 7.42 9.46 9.46 7.43 
3P 2p2 7.46 7.40 7.45 NA NA NA 6.19 6.19 6.21 6.21 8.64 NA 7.46 7.45 7.84 7.83 7.46 
MSEf … 0.04 0.06 −0.67 −0.65 −0.75 −1.26 −1.21 −1.24 −1.18 0.36 −0.38 0.07 0.06 1.00 1.00 0.06 
MAEg … 0.07 0.07 0.67 0.65 0.75 1.26 1.21 1.24 1.18 0.95 0.38 0.07 0.07 1.00 1.00 0.07 
    EOM-   2ph-2ph-r2ph-r2ph-    pp-RPApp-TDAr2pp-
SymbolConfigurationExpt.aFCICCSDCISTDHFbTDDFTcRPAbTDARPAbTDASOPPARPA(D)dpp-RPApp-TDA@B3LYPe@B3LYPeTDA
3P 2s12p1 2.73 2.72 2.73 1.70 −0.96i 2.10 −1.27i 1.45 −1.22i 1.49 1.81 NA 2.73 2.73 2.82 2.81 2.73 
1P 2s12p1 5.28 5.33 5.36 5.10 4.84 4.90 3.74 4.07 3.77 4.10 4.86 4.87 5.36 5.36 6.15 6.15 5.34 
3S 2s13s1 6.46 6.43 6.44 5.53 5.49 5.69 5.12 5.16 5.15 5.19 6.47 NA 6.44 6.44 8.40 8.40 6.44 
1S 2s13s1 6.78 6.77 6.77 6.17 6.16 5.97 5.45 5.47 5.47 5.49 6.01 6.48 6.77 6.77 8.68 8.68 6.77 
1D 2p2 7.05 7.16 7.18 NA NA NA 5.92 5.92 5.94 5.94 8.83 NA 7.18 7.18 7.98 7.97 7.18 
3P 2s13p1 7.30 7.42 7.42 6.56 6.55 6.08 6.15 6.16 6.17 6.18 7.02 NA 7.43 7.42 9.46 9.46 7.43 
3P 2p2 7.46 7.40 7.45 NA NA NA 6.19 6.19 6.21 6.21 8.64 NA 7.46 7.45 7.84 7.83 7.46 
MSEf … 0.04 0.06 −0.67 −0.65 −0.75 −1.26 −1.21 −1.24 −1.18 0.36 −0.38 0.07 0.06 1.00 1.00 0.06 
MAEg … 0.07 0.07 0.67 0.65 0.75 1.26 1.21 1.24 1.18 0.95 0.38 0.07 0.07 1.00 1.00 0.07 
a

Experimental values from NIST.76 

b

The HF solution of Be atom suffers from instability. Therefore, TDHF, the 2ph-RPA, and the r2ph-RPA have imaginary eigenvalues presented.

c

TDDFT with B3LYP functional.

d

Due to the limitation of the program Dalton, only singlet data are listed.

e

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

f

Mean signed error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA and imaginary values are excluded.

g

Mean absolute error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA and imaginary values are excluded.

The BH molecule has both low-lying double excitations and non-HOMO single excitations. Table V lists some lowest excitations. The ground state configuration for BH is 1σ222, where 1σ is the 1s orbital of the B atom. FCI excitation energies are used as the reference owing to the incompleteness of experimental data.77 Note that EOM-CCSD is very accurate for single excitations but not so at all double excitations. It is suggested that although the EOM-CCSD formalism contains double excitations, the quantitative description is not as good as that of single excitations.78 Expectantly, CIS, TDHF, and TDDFT miss quite a few double excitations from 3σ to π. Again, the instability of the HF and B3LYP reference leads to imaginary eigenvalues for TDHF, TDDFT, and r2ph-RPA. A systematic downshift of r2ph-TDA excitation energies is observed. The pp-RPA and the pp-TDA have similar results with an MAE of about 0.5 eV. The r2pp-TDA systematically pulls up pp-TDA results, thus overestimating most excitations. The pp-RPA and the pp-TDA with B3LYP reference still present unbalanced excitation energies of low and high excitations. Overall, the r2ph-TDA, the pp-RPA, the pp-TDA, and the r2pp-TDA perform reasonably well for states studied in Table V. SOPPA significantly overestimates the first excitation energy, probably due to the singlet-triplet instability of the HF reference. SOPPA also has a very large error on the double excitation dominated transition 3σ2 → π2. Overall, SOPPA and RPA(D) are fairly accurate for most single dominated excitations.

Table V.

Vertical excitation energies of BH with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 MRCI in MOLPRO,67,68 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

     EOM-   r2ph-r2ph-    pp-RPApp-TDAr2pp-
TermTransitionExpt.aFCIMRCICCSDCISTDHFbTDDFTb,cRPAbTDASOPPARPA(D)dpp-RPApp-TDA@B3LYPe@B3LYPeTDA
a3Π 3σ → π … 1.34 1.33 1.33 0.55 −1.63i −0.43i −1.63i 0.24 7.35 NA 1.64 1.61 1.38 1.32 2.40 
A1Π 3σ → π 2.87 2.90 2.87 2.96 2.86 2.66 2.69 1.81 2.08 2.46 2.46 3.19 3.16 3.19 3.13 4.21 
b3Σ 2 → π2 … 4.71 4.69 5.69 NA NA NA 4.42 4.42 9.33 NA 5.50 5.47 5.12 5.06 6.56 
C1Δ 2 → π2 5.72 5.78 5.73 6.74 NA NA NA 5.09 5.09 6.37 6.38 6.14 6.11 6.04 5.98 7.20 
c3Σ+ 3σ → 3s σ … 6.37 6.39 6.38 6.05 NA 5.58 5.34 5.39 6.18 NA 5.58 5.55 7.64 7.58 6.52 
B1Σ+ 3σ → 3s σ 6.49 6.49 6.50 6.52 6.42 NA 5.66 5.53 5.54 7.50 7.50 5.70 5.67 8.41 8.36 6.65 
C1Σ+ 2 → π2 6.86 6.92 6.89 7.47 NA 7.23 NA 5.93 5.93 8.78 8.78 6.83 6.80 7.09 7.14 7.72 
23Σ+ 3σ → 3p σ … 7.33 7.34 7.33 6.87 7.82 6.37 6.21 6.30 7.06 NA 6.61 6.58 9.02 8.96 7.47 
23Π 3σ → 3p π … 7.54 7.55 7.56 7.27 8.29 6.56 6.51 6.54 8.00 NA 6.88 6.85 9.90 9.84 7.76 
D1Π 3σ → 3p π … 7.61 7.60 7.62 7.44 10.89 6.61 6.59 6.60 7.45 7.46 6.97 6.94 10.04 9.98 7.88 
MSEf … … −0.01 0.24 −0.31 −0.28 −0.77 3.57 −0.90 1.33 0.36 −0.15 −0.18 1.08 1.03 0.78 
MAEg … … 0.02 0.25 0.31 0.28 0.77 3.63 0.90 1.55 0.66 0.51 0.52 1.08 1.03 0.78 
     EOM-   r2ph-r2ph-    pp-RPApp-TDAr2pp-
TermTransitionExpt.aFCIMRCICCSDCISTDHFbTDDFTb,cRPAbTDASOPPARPA(D)dpp-RPApp-TDA@B3LYPe@B3LYPeTDA
a3Π 3σ → π … 1.34 1.33 1.33 0.55 −1.63i −0.43i −1.63i 0.24 7.35 NA 1.64 1.61 1.38 1.32 2.40 
A1Π 3σ → π 2.87 2.90 2.87 2.96 2.86 2.66 2.69 1.81 2.08 2.46 2.46 3.19 3.16 3.19 3.13 4.21 
b3Σ 2 → π2 … 4.71 4.69 5.69 NA NA NA 4.42 4.42 9.33 NA 5.50 5.47 5.12 5.06 6.56 
C1Δ 2 → π2 5.72 5.78 5.73 6.74 NA NA NA 5.09 5.09 6.37 6.38 6.14 6.11 6.04 5.98 7.20 
c3Σ+ 3σ → 3s σ … 6.37 6.39 6.38 6.05 NA 5.58 5.34 5.39 6.18 NA 5.58 5.55 7.64 7.58 6.52 
B1Σ+ 3σ → 3s σ 6.49 6.49 6.50 6.52 6.42 NA 5.66 5.53 5.54 7.50 7.50 5.70 5.67 8.41 8.36 6.65 
C1Σ+ 2 → π2 6.86 6.92 6.89 7.47 NA 7.23 NA 5.93 5.93 8.78 8.78 6.83 6.80 7.09 7.14 7.72 
23Σ+ 3σ → 3p σ … 7.33 7.34 7.33 6.87 7.82 6.37 6.21 6.30 7.06 NA 6.61 6.58 9.02 8.96 7.47 
23Π 3σ → 3p π … 7.54 7.55 7.56 7.27 8.29 6.56 6.51 6.54 8.00 NA 6.88 6.85 9.90 9.84 7.76 
D1Π 3σ → 3p π … 7.61 7.60 7.62 7.44 10.89 6.61 6.59 6.60 7.45 7.46 6.97 6.94 10.04 9.98 7.88 
MSEf … … −0.01 0.24 −0.31 −0.28 −0.77 3.57 −0.90 1.33 0.36 −0.15 −0.18 1.08 1.03 0.78 
MAEg … … 0.02 0.25 0.31 0.28 0.77 3.63 0.90 1.55 0.66 0.51 0.52 1.08 1.03 0.78 
a

Experimental values from Ref. 77.

b

Both HF and B3LYP solutions of BH molecule suffer from instability. Therefore, TDHF, TDDFT, and the r2ph-RPA have imaginary eigenvalues presented.

c

TDDFT with B3LYP functional.

d

Due to the limitation of the program Dalton, only singlet data are listed.

e

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

f

Mean signed error. Errors are with respect to FCI values. Spatial multiplicity is accounted. NA and imaginary values are excluded.

g

Mean absolute error. Errors are with respect to FCI values. Spatial multiplicity is accounted. NA and imaginary values are excluded.

Various methods are tested on the CO molecule with results tabulated in Table VI. No double excitations are involved in the range of spectra of interest. The pp-RPA and the pp-TDA with ground state (N − 2)-electron reference do not capture single excitations from non-HOMO orbitals (the HOMO is σ), such as those π → π* transitions. The r2pp-TDA qualitatively captures these excitations, however, systematically overestimating all excitation energies. For this system, the HF reference is stable, and the r2ph-TDA has the same level of accuracy of CIS, TDHF, and TDDFT, with relatively low average errors. The accuracy of SOPPA and RPA(D) is close to the r2ph-TDA.

Table VI.

Vertical excitation energies of CO with different methods using the basis set aug-cc-pVTZspd with no f functions. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

             pp-RPApp-TDA 
StateTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
a3Π σ → π* 6.32 6.26 5.67 5.07 5.69 4.73 5.34 5.74 NA 5.57 5.55 5.82 5.67 9.12 
a3Σ+ π → π* 8.51 7.98 7.34 5.79 7.51 5.77 7.33 7.65 NA NA NA NA NA 12.85 
A1Π σ → π* 8.51 8.52 8.84 8.55 8.21 7.74 8.04 8.03 8.19 7.82 7.83 8.03 8.03 11.54 
d3Δ π → π* 9.36 8.92 8.28 7.37 8.23 7.35 8.27 8.49 NA NA NA NA NA 13.90 
e3Σ π → π* 9.88 9.49 9.27 8.89 9.30 8.77 9.15 9.15 NA NA NA NA NA 14.82 
I1Σ π → π* 9.88 9.70 9.27 8.89 9.30 8.88 9.26 9.19 9.56 NA NA NA NA 15.03 
D1Δ π → π* 10.23 9.82 9.69 9.49 9.63 9.38 9.57 9.42 9.73 NA NA NA NA 15.35 
b3Σ+ σ → 3s 10.40 10.65 11.20 11.10 9.89 10.59 10.68 10.72 NA 8.96 8.88 12.24 11.96 12.50 
B1Σ+ σ → 3s 10.78 11.13 12.18 12.14 10.37 11.19 11.23 11.14 11.33 9.41 9.34 12.83 12.57 12.88 
j3Σ+ σ → 3p σ 11.30 11.49 12.44 12.42 10.62 11.75 11.77 11.54 NA 9.76 9.68 13.72 13.44 13.28 
C1Σ+ σ → 3p σ 11.40 11.68 12.78 12.76 10.77 11.95 11.96 11.73 11.65 10.02 9.94 14.03 13.75 13.46 
c3Π σ → 3p π 11.55 11.77 12.57 12.39 10.87 11.88 11.93 11.35 NA 10.06 9.98 13.90 13.62 13.59 
E1Π σ → 3p π 11.53 11.93 12.89 12.88 10.98 12.06 12.08 11.93 11.97 10.19 10.12 13.95 13.67 13.76 
MSEe … −0.03 0.17 −0.19 −0.64 −0.63 −0.28 −0.32 −0.03 −1.19 −1.24 1.38 1.17 3.27 
MAEf … 0.28 0.90 1.14 0.64 0.98 0.66 0.54 0.40 1.19 1.24 1.71 1.54 3.27 
             pp-RPApp-TDA 
StateTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
a3Π σ → π* 6.32 6.26 5.67 5.07 5.69 4.73 5.34 5.74 NA 5.57 5.55 5.82 5.67 9.12 
a3Σ+ π → π* 8.51 7.98 7.34 5.79 7.51 5.77 7.33 7.65 NA NA NA NA NA 12.85 
A1Π σ → π* 8.51 8.52 8.84 8.55 8.21 7.74 8.04 8.03 8.19 7.82 7.83 8.03 8.03 11.54 
d3Δ π → π* 9.36 8.92 8.28 7.37 8.23 7.35 8.27 8.49 NA NA NA NA NA 13.90 
e3Σ π → π* 9.88 9.49 9.27 8.89 9.30 8.77 9.15 9.15 NA NA NA NA NA 14.82 
I1Σ π → π* 9.88 9.70 9.27 8.89 9.30 8.88 9.26 9.19 9.56 NA NA NA NA 15.03 
D1Δ π → π* 10.23 9.82 9.69 9.49 9.63 9.38 9.57 9.42 9.73 NA NA NA NA 15.35 
b3Σ+ σ → 3s 10.40 10.65 11.20 11.10 9.89 10.59 10.68 10.72 NA 8.96 8.88 12.24 11.96 12.50 
B1Σ+ σ → 3s 10.78 11.13 12.18 12.14 10.37 11.19 11.23 11.14 11.33 9.41 9.34 12.83 12.57 12.88 
j3Σ+ σ → 3p σ 11.30 11.49 12.44 12.42 10.62 11.75 11.77 11.54 NA 9.76 9.68 13.72 13.44 13.28 
C1Σ+ σ → 3p σ 11.40 11.68 12.78 12.76 10.77 11.95 11.96 11.73 11.65 10.02 9.94 14.03 13.75 13.46 
c3Π σ → 3p π 11.55 11.77 12.57 12.39 10.87 11.88 11.93 11.35 NA 10.06 9.98 13.90 13.62 13.59 
E1Π σ → 3p π 11.53 11.93 12.89 12.88 10.98 12.06 12.08 11.93 11.97 10.19 10.12 13.95 13.67 13.76 
MSEe … −0.03 0.17 −0.19 −0.64 −0.63 −0.28 −0.32 −0.03 −1.19 −1.24 1.38 1.17 3.27 
MAEf … 0.28 0.90 1.14 0.64 0.98 0.66 0.54 0.40 1.19 1.24 1.71 1.54 3.27 
a

Experimental values from Ref. 79.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA values are excluded.

f

Mean absolute error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA values are excluded.

We also test our methods on N2 with the basis set aug-cc-pVTZspd. The HF ground state configuration of a neutral N2 molecule is

$\sigma _{u}^{2}\sigma _{g}^{2}\pi _{u}^{4}$
σu2σg2πu4⁠, while its doubly charged cation has a ground state configuration of
$\sigma _{u}^{2}\pi _{u}^{4}$
σu2πu4
. In other words, the HOMO in the N-electron HF configuration is πu, but the molecule will lose two σg electrons to form the lowest doubly charged cation. We think this is considered an artifact of HF orbitals, as σg is generally considered the HOMO for N2,80 which is also confirmed by a B3LYP calculation. That said, the pp-RPA and the pp-TDA with HF or B3LYP reference are blind to all πu and σu excitations. The r2pp-TDA can capture but seriously overestimate all excitations listed in Table VII. CIS, TDHF, TDDFT, and the r2ph-TDA all tend to underestimate low-lying excitation energies but overestimate higher excitation energies. Note that the HF HOMOs are two degenerate πu orbitals, thus the double excitations are restricted to the two πu orbitals. N2 is the only example with degenerate HOMOs in this study. In this case, SOPPA and RPA(D) perform surprisingly well for a wide range of states.

Table VII.

Vertical excitation energies of N2 with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

            pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
$A{}^{3}\Sigma _{u}^{+}$
A3Σu+
 
πu → πg 7.75 7.04 5.48 1.66 6.38 4.88 6.88 NA NA NA NA NA 9.25 
B3Πg σg → πg 8.04 7.76 7.52 7.11 7.18 9.48 7.18 NA 7.82 7.54 7.16 7.16 11.17 
W3Δu πu → πg 8.88 8.40 6.57 4.86 7.32 5.84 8.09 NA NA NA NA NA 10.29 
a1Πg σg → πg 9.31 9.12 9.50 9.23 8.85 9.48 8.49 8.78 9.34 9.19 9.07 9.07 12.77 
$B^{\prime }{}^{3}\Sigma _{u}^{-}$
BΣu3
 
πu → πg 9.67 9.32 7.75 7.11 8.68 6.79 9.05 NA NA NA NA NA 11.45 
$a^{\prime }{}^{1}\Sigma _{u}^{-}$
aΣu1
 
πu → πg 9.92 9.57 7.75 7.11 8.68 6.88 9.10 9.69 NA NA NA NA 11.44 
w1Δu πu → πg 10.27 10.00 8.33 8.01 9.10 7.30 9.58 10.00 NA NA NA NA 11.99 
C3Πu σu → πg 11.19 11.16 11.56 11.09 10.49 11.39 10.59 NA NA NA NA NA 15.62 
$E{}^{3}\Sigma _{g}^{+}$
EΣg+3
 
σg → 3s σg 12.00 12.18 13.22 13.14 11.57 12.70 12.18 NA 11.33 10.94 13.80 13.80 14.38 
$a^{\prime \prime }{}^{1}\Sigma _{g}^{+}$
aΣg+1
 
σg → 3s σg 12.20 12.76 14.44 14.08 12.10 14.54 12.76 14.75 11.81 11.45 14.40 14.40 15.09 
c1Πu σg → 3p πu 12.90 13.25 14.91 15.45 12.64 14.94 13.05 14.09 12.35 11.95 15.25 15.25 15.84 
$c^{\prime }{}^{1}\Sigma _{u}^{+}$
cΣu+1
 
σg → 3p σu 12.98 14.52 14.24 14.90 12.22 13.35 12.87 13.89 12.04 11.64 15.24 15.24 15.33 
MSEe … −0.05 −0.33 −0.90 −0.83 −0.54 −0.49 0.45 −0.39 −0.70 0.97 0.97 2.59 
MAEf … 0.38 1.43 2.02 0.83 1.77 0.61 0.85 0.40 0.70 1.47 1.47 2.59 
            pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
$A{}^{3}\Sigma _{u}^{+}$
A3Σu+
 
πu → πg 7.75 7.04 5.48 1.66 6.38 4.88 6.88 NA NA NA NA NA 9.25 
B3Πg σg → πg 8.04 7.76 7.52 7.11 7.18 9.48 7.18 NA 7.82 7.54 7.16 7.16 11.17 
W3Δu πu → πg 8.88 8.40 6.57 4.86 7.32 5.84 8.09 NA NA NA NA NA 10.29 
a1Πg σg → πg 9.31 9.12 9.50 9.23 8.85 9.48 8.49 8.78 9.34 9.19 9.07 9.07 12.77 
$B^{\prime }{}^{3}\Sigma _{u}^{-}$
BΣu3
 
πu → πg 9.67 9.32 7.75 7.11 8.68 6.79 9.05 NA NA NA NA NA 11.45 
$a^{\prime }{}^{1}\Sigma _{u}^{-}$
aΣu1
 
πu → πg 9.92 9.57 7.75 7.11 8.68 6.88 9.10 9.69 NA NA NA NA 11.44 
w1Δu πu → πg 10.27 10.00 8.33 8.01 9.10 7.30 9.58 10.00 NA NA NA NA 11.99 
C3Πu σu → πg 11.19 11.16 11.56 11.09 10.49 11.39 10.59 NA NA NA NA NA 15.62 
$E{}^{3}\Sigma _{g}^{+}$
EΣg+3
 
σg → 3s σg 12.00 12.18 13.22 13.14 11.57 12.70 12.18 NA 11.33 10.94 13.80 13.80 14.38 
$a^{\prime \prime }{}^{1}\Sigma _{g}^{+}$
aΣg+1
 
σg → 3s σg 12.20 12.76 14.44 14.08 12.10 14.54 12.76 14.75 11.81 11.45 14.40 14.40 15.09 
c1Πu σg → 3p πu 12.90 13.25 14.91 15.45 12.64 14.94 13.05 14.09 12.35 11.95 15.25 15.25 15.84 
$c^{\prime }{}^{1}\Sigma _{u}^{+}$
cΣu+1
 
σg → 3p σu 12.98 14.52 14.24 14.90 12.22 13.35 12.87 13.89 12.04 11.64 15.24 15.24 15.33 
MSEe … −0.05 −0.33 −0.90 −0.83 −0.54 −0.49 0.45 −0.39 −0.70 0.97 0.97 2.59 
MAEf … 0.38 1.43 2.02 0.83 1.77 0.61 0.85 0.40 0.70 1.47 1.47 2.59 
a

Experimental values from Ref. 81.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA values are excluded.

f

Mean absolute error. Errors are with respect to experimental values. Spatial multiplicity is accounted. NA values are excluded.

Calculation results of CH2O are listed in Table VIII. The ground state configuration of CH2O in the basis set aug-cc-pVTZspd is σ2π2n2. The r2pp-TDA recovers the non-HOMO excitations missing in all pp-RPA(TDA) methods, but with a large average error of above 2 eV. CIS tends to overestimate most of the excitations shown in the table, while the systematic downshift in the r2ph-TDA greatly improves the accuracy. Again, SOPPA and RPA(D) have similar accuracy compared to that of the r2ph-TDA.

Table VIII.

Vertical excitation energies of CH2O with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

             pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3A2 n → π* 3.50 3.48 3.67 3.33 3.11 3.03 3.37 2.88 NA 1.64 1.51 3.17 2.87 6.69 
1A2 n → π* 3.94 3.94 4.50 4.32 3.84 3.98 4.16 3.36 3.58 1.98 1.86 3.70 3.45 6.97 
3A1 π → π* 5.53 5.77 4.68 1.38 5.21 1.38 4.68 5.45 NA NA NA NA NA 8.11 
3B2 n → 3sa1 6.83 7.00 8.27 8.18 6.36 7.69 7.76 6.2 NA 3.71 3.55 7.42 6.98 8.73 
1B2 n → 3sa1 7.09 7.15 8.62 8.61 6.49 8.00 8.01 6.32 6.51 3.82 3.66 7.92 7.47 8.80 
3A1 n → 3pb2 7.79 8.02 9.31 9.25 7.24 8.74 8.78 7.23 NA 4.75 4.59 8.81 8.36 9.67 
1A1 n → 3pb2 7.97 8.13 9.46 9.61 7.32 8.91 8.98 7.31 9.12 4.91 4.75 9.48 9.05 9.72 
3B2 n → 3pa1 7.96 7.83 9.05 8.99 7.33 8.48 8.53 7.22 NA 4.98 4.83 8.92 8.50 9.77 
1B2 n → 3pa1 8.12 8.03 9.41 9.40 7.47 8.81 8.82 7.35 7.42 5.08 4.92 9.15 8.72 9.89 
1B1 σ → π* 8.68 9.20 9.69 9.44 8.86 9.42 9.67 8.55 8.87 NA NA NA NA 11.96 
1A2 n → 3pb1 8.38 8.57 10.07 10.06 8.02 9.50 9.51 7.83 7.93 5.65 5.50 10.17 9.72 10.52 
MSEe … 0.12 0.99 0.62 −0.41 0.20 0.59 −0.55 −0.13 −2.79 −2.93 0.79 0.39 2.28 
MAEf … 0.17 1.15 1.40 0.44 1.04 0.77 0.55 0.57 2.79 2.93 0.92 0.64 2.28 
             pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3A2 n → π* 3.50 3.48 3.67 3.33 3.11 3.03 3.37 2.88 NA 1.64 1.51 3.17 2.87 6.69 
1A2 n → π* 3.94 3.94 4.50 4.32 3.84 3.98 4.16 3.36 3.58 1.98 1.86 3.70 3.45 6.97 
3A1 π → π* 5.53 5.77 4.68 1.38 5.21 1.38 4.68 5.45 NA NA NA NA NA 8.11 
3B2 n → 3sa1 6.83 7.00 8.27 8.18 6.36 7.69 7.76 6.2 NA 3.71 3.55 7.42 6.98 8.73 
1B2 n → 3sa1 7.09 7.15 8.62 8.61 6.49 8.00 8.01 6.32 6.51 3.82 3.66 7.92 7.47 8.80 
3A1 n → 3pb2 7.79 8.02 9.31 9.25 7.24 8.74 8.78 7.23 NA 4.75 4.59 8.81 8.36 9.67 
1A1 n → 3pb2 7.97 8.13 9.46 9.61 7.32 8.91 8.98 7.31 9.12 4.91 4.75 9.48 9.05 9.72 
3B2 n → 3pa1 7.96 7.83 9.05 8.99 7.33 8.48 8.53 7.22 NA 4.98 4.83 8.92 8.50 9.77 
1B2 n → 3pa1 8.12 8.03 9.41 9.40 7.47 8.81 8.82 7.35 7.42 5.08 4.92 9.15 8.72 9.89 
1B1 σ → π* 8.68 9.20 9.69 9.44 8.86 9.42 9.67 8.55 8.87 NA NA NA NA 11.96 
1A2 n → 3pb1 8.38 8.57 10.07 10.06 8.02 9.50 9.51 7.83 7.93 5.65 5.50 10.17 9.72 10.52 
MSEe … 0.12 0.99 0.62 −0.41 0.20 0.59 −0.55 −0.13 −2.79 −2.93 0.79 0.39 2.28 
MAEf … 0.17 1.15 1.40 0.44 1.04 0.77 0.55 0.57 2.79 2.93 0.92 0.64 2.28 
a

Experimental values from Ref. 79.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Errors are with respect to experimental values. NA values are excluded.

f

Mean absolute error. Errors are with respect to experimental values. NA values are excluded.

Table IX tabulates excitation energies from the HOMO π orbital of ethylene from various methods. The pp-RPA and the pp-TDA systematically underestimate all the excitations while the r2pp-TDA does the opposite. The pp-RPA and the pp-TDA with B3LYP reference produce better accuracy. CIS performs in the same level as EOM-CCSD in this special case. The first triplet excitation energy of TDHF is very low, probably because the system is close to the stability-instability transition. The r2ph-TDA this time underestimates most excitations by about 0.6 eV, similar to the error of TDDFT. For this molecular, surprisingly, SOPPA and RPA(D) have smaller errors than EOM-CCSD.

Table IX.

Vertical excitation energies of C2H4 with different methods using the basis set aug-cc-pVTZspd. All calculations are done in QM4D,65 except for FCI in GAMESS (US),66 and EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 The molecule lies in the yz plane with the C=C bond aligned at the z axis. All excitations are from the HOMO π orbital. Unit: eV.

             pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3B3u 3s = ag 6.98 7.24 6.91 6.87 6.54 6.51 7.05 NA 4.42 4.40 7.72 7.68 8.41   
1B3u 3s = ag 7.11 7.37 7.13 7.12 6.62 6.62 7.17 7.18 4.51 4.49 7.81 7.77 8.51   
3B1g 3p σ = b2u 7.79 7.98 7.63 7.60 7.13 7.23 7.75 NA 5.09 5.07 8.43 8.39 9.09   
1B1g 3p σ = b2u 7.80 8.03 7.73 7.71 7.17 7.30 7.82 7.83 5.09 5.09 8.50 8.46 9.11   
1B2g 3p σ = b1u 7.90 8.08 7.88 7.88 7.16 7.40 7.88 7.89 5.08 5.06 8.53 8.49 9.08   
1B1u π* = b2g 8.00 8.01 7.71 7.36 7.37 6.99 7.43 7.51 6.14 6.12 8.44 8.43 9.94   
3Ag 3p π = b3u 8.15 8.51 7.97 7.93 7.89 7.69 8.35 NA 5.76 5.73 9.47 9.43 9.74   
1Ag 3p π = b3u 8.28 8.80 8.53 8.49 8.15 7.86 8.61 8.62 6.28 6.26 9.90 9.87 9.98   
MSEe … 0.23 −0.14 −0.54 −0.48 −0.60 −0.04 −0.01 −2.23 −2.25 0.68 0.64 1.70   
MAEf … 0.23 0.20 0.58 0.48 0.60 0.19 0.19 2.23 2.25 0.84 0.81 1.70   
             pp-RPApp-TDA 
SymbolTransitionExpt.aEOM-CCSDCISTDHFTDDFTbr2ph-RPAr2ph-TDASOPPARPA(D)cpp-RPApp-TDA@B3LYPd@B3LYPdr2pp-TDA
3B3u 3s = ag 6.98 7.24 6.91 6.87 6.54 6.51 7.05 NA 4.42 4.40 7.72 7.68 8.41   
1B3u 3s = ag 7.11 7.37 7.13 7.12 6.62 6.62 7.17 7.18 4.51 4.49 7.81 7.77 8.51   
3B1g 3p σ = b2u 7.79 7.98 7.63 7.60 7.13 7.23 7.75 NA 5.09 5.07 8.43 8.39 9.09   
1B1g 3p σ = b2u 7.80 8.03 7.73 7.71 7.17 7.30 7.82 7.83 5.09 5.09 8.50 8.46 9.11   
1B2g 3p σ = b1u 7.90 8.08 7.88 7.88 7.16 7.40 7.88 7.89 5.08 5.06 8.53 8.49 9.08   
1B1u π* = b2g 8.00 8.01 7.71 7.36 7.37 6.99 7.43 7.51 6.14 6.12 8.44 8.43 9.94   
3Ag 3p π = b3u 8.15 8.51 7.97 7.93 7.89 7.69 8.35 NA 5.76 5.73 9.47 9.43 9.74   
1Ag 3p π = b3u 8.28 8.80 8.53 8.49 8.15 7.86 8.61 8.62 6.28 6.26 9.90 9.87 9.98   
MSEe … 0.23 −0.14 −0.54 −0.48 −0.60 −0.04 −0.01 −2.23 −2.25 0.68 0.64 1.70   
MAEf … 0.23 0.20 0.58 0.48 0.60 0.19 0.19 2.23 2.25 0.84 0.81 1.70   
a

Experimental values from Ref. 79.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Errors are with respect to experimental values.

f

Mean absolute error. Errors are with respect to experimental values.

We also test different methods on E-butadiene (Table X) with the basis set aug-cc-pVDZ. Due to the relatively small size of the basis set, we only look at the lowest two singlet and triplet excitations. Note that the HF solution of this system is unstable with an imaginary singlet-triplet transition energy. The pp-RPA and the pp-TDA with HF and B3LYP reference all have relatively low average errors of about 0.5 eV. The r2pp-TDA again systematically overestimates all excitations. CIS, TDDFT, and the r2ph-TDA are of similar accuracy for the four states studied. The accuracy of SOPPA and RPA(D) is not quite different from that of the r2ph-TDA and the r2ph-RPA.

Table X.

Vertical excitation energies of E-Butadiene with different methods using the basis set aug-cc-pVDZ. All calculations are done in QM4D,65 except for EOM-CCSD in Gaussian 09,69 and SOPPA and RPA(D) in Dalton.74 Unit: eV.

StateRef.aEOM-CCSDCISTDHFTDDFTbr2ph-TDASOPPARPA(D)cpp-RPApp-TDApp-RPA @B3LYPdpp-TDA @B3LYPdr2pp-TDA
3Bu(π → π*3.20 3.29 2.64 −2.17i 2.80 2.51 2.81 NA 3.22 3.12 2.54 2.32 5.55 
3Ag(π → π*5.08 5.16 4.34 2.97 4.86 4.29 4.69 NA 5.60 5.50 6.08 5.86 7.70 
1Bu(π → π*6.18 6.35 6.20 5.91 5.57 5.78 5.62 5.81 5.49 5.38 6.58 6.38 8.11 
1Ag(π → π*6.55 7.06 7.39 7.27 6.50 6.79 6.78 7.05 5.92 5.83 6.45 6.32 8.59 
MSEe … 0.21 −0.11 −0.55 −0.32 −0.41 −0.28 0.07 −0.19 −0.30 0.16 −0.03 2.23 
MAEf … 0.21 0.54 1.03 0.32 0.53 0.39 0.44 0.47 0.50 0.54 0.53 2.23 
StateRef.aEOM-CCSDCISTDHFTDDFTbr2ph-TDASOPPARPA(D)cpp-RPApp-TDApp-RPA @B3LYPdpp-TDA @B3LYPdr2pp-TDA
3Bu(π → π*3.20 3.29 2.64 −2.17i 2.80 2.51 2.81 NA 3.22 3.12 2.54 2.32 5.55 
3Ag(π → π*5.08 5.16 4.34 2.97 4.86 4.29 4.69 NA 5.60 5.50 6.08 5.86 7.70 
1Bu(π → π*6.18 6.35 6.20 5.91 5.57 5.78 5.62 5.81 5.49 5.38 6.58 6.38 8.11 
1Ag(π → π*6.55 7.06 7.39 7.27 6.50 6.79 6.78 7.05 5.92 5.83 6.45 6.32 8.59 
MSEe … 0.21 −0.11 −0.55 −0.32 −0.41 −0.28 0.07 −0.19 −0.30 0.16 −0.03 2.23 
MAEf … 0.21 0.54 1.03 0.32 0.53 0.39 0.44 0.47 0.50 0.54 0.53 2.23 
a

Best estimated theoretical values adopted from Ref. 82.

b

TDDFT with B3LYP functional.

c

Due to the limitation of the program Dalton, only singlet data are listed.

d

The pp-RPA and the pp-TDA using B3LYP reference as in Ref. 39.

e

Mean signed error. Imaginary values are excluded.

f

Mean absolute error. Imaginary values are excluded.

Results above show that the r2ph-TDA usually has similar errors compared to TDDFT, while capable of describing double excitations. We now study the detailed error distributions of the r2ph-TDA, compared to EOM-CCSD, TDDFT, CIS, SOPPA, and RPA(D) in Figure 1. These aggregated results are also fitted to linear regressions, with results shown in Table XI. Figure 1 shows that N2 results are the main outliers of r2ph-TDA results, thus another set regressions without N2 are also performed. Overall, EOM-CCSD has high correlation coefficients (R2) above 0.99 with very little bias (|b| < 0.2 eV). TDDFT also has high correlation coefficients (R2 > 0.98), yet with a bias of about −0.5 eV. The R2 value for the r2ph-TDA without N2 (0.973) is much better than that including N2 (0.944), albeit the change of the bias is not significant. CIS, on which the r2ph-TDA is a correction based, also undergoes improvement through excluding N2 from the data. N2 results have the largest errors probably because of the unphysical description of the HOMO in the HF reference. Interestingly, SOPPA and RPA(D) are very accurate for N2. The most significant outliers of SOPPA and RPA(D) is BH. The regression results show that these two methods have better overall accuracy and small biases. Although the r2ph-TDA is slightly less accurate than TDDFT in these test data, the r2ph-TDA also has advantage of capturing some double excitations and free from instability. Considering double excitations, the r2ph-TDA is still preferable to SOPPA and RPA(D). The negative bias of the r2ph-TDA data (∼−1.5 eV) is compatible with results of the 2ph-RPA in literature as a result of some missing ground state correlation.43,44,83

FIG. 1.

Error distributions of various methods. The reference values are either experimental values, or other theoretical values. See Tables III–X for details. Data of H2O are those using the basis set aug-cc-pVTZspd. NA and imaginary excitation energies are excluded. (a) EOM-CCSD, (b) TDDFT(B3LYP), (c) The r2ph-TDA, (d) CIS, (e) SOPPA, (f) RPA(D).

FIG. 1.

Error distributions of various methods. The reference values are either experimental values, or other theoretical values. See Tables III–X for details. Data of H2O are those using the basis set aug-cc-pVTZspd. NA and imaginary excitation energies are excluded. (a) EOM-CCSD, (b) TDDFT(B3LYP), (c) The r2ph-TDA, (d) CIS, (e) SOPPA, (f) RPA(D).

Close modal
Table XI.

Linear regression results of excitation energies. The equation to fit is y = ax + b. Spatial multiplicity is accounted. NA and imaginary values are excluded.

 AllWithout N2
MethodsR2ab (eV)R2ab (eV)
EOM-CCSD 0.992 1.003 0.048 0.994 0.991 0.160 
TDDFT(B3LYP) 0.987 0.986 −0.559 0.986 0.976 −0.456 
r2ph-TDA 0.944 1.143 −1.624 0.973 1.142 −1.526 
CIS 0.952 1.143 −1.133 0.971 1.158 −1.068 
SOPPA 0.886 0.864 1.073 0.849 0.842 1.263 
SOPPA without BH 0.962 1.025 −0.418 0.954 1.043 −0.465 
RPA(D) 0.972 1.102 −0.778 0.970 1.030 −0.256 
 AllWithout N2
MethodsR2ab (eV)R2ab (eV)
EOM-CCSD 0.992 1.003 0.048 0.994 0.991 0.160 
TDDFT(B3LYP) 0.987 0.986 −0.559 0.986 0.976 −0.456 
r2ph-TDA 0.944 1.143 −1.624 0.973 1.142 −1.526 
CIS 0.952 1.143 −1.133 0.971 1.158 −1.068 
SOPPA 0.886 0.864 1.073 0.849 0.842 1.263 
SOPPA without BH 0.962 1.025 −0.418 0.954 1.043 −0.465 
RPA(D) 0.972 1.102 −0.778 0.970 1.030 −0.256 

Second RPAs and second TDAs of ph- and pp-channels are introduced in this article to study molecular excitations. These extensions enable capturing double excitations in the RPA/TDA and non-HOMO excitations in the pp-RPA/TDA. Based on orbital restrictions, the r2ph-RPA, the r2ph-TDA, and the r2pp-TDA can describe all single excitations, and double excitations from HOMO, with a formal scaling of O(N4). Since the r2ph-RPA and the 2ph-RPA inherit all instability issues from TDHF, we suggest that r2ph-TDA is preferable than the r2ph-RPA. In theory, the r2pp-TDA and the r2ph-TDA can capture the same states according to the formalism. However, r2pp-TDA usually overestimates the excitation energies. Additionally, the possibility that the (N − 2)-electron reference may be degenerate and that the SCF result could deteriorate the symmetry. Moreover, pp series methods are almost impossible to be size extensive, while ph series methods are probably size extensive according our preliminary data.54 These restricted second RPAs and TDAs are tested on various systems. The r2ph-TDA have similar results compared to TDDFT and CIS, but with a larger negative bias, which indicates some ground state correlation energies are unaccounted in the r2ph-TDA. SOPPA and RPA(D) are very accurate overall, but for double dominated excitations, their errors are still fairly large. Considering EOM-CCSD is unbalanced for single and double excitations even with a high scaling of O(N6), and that SOPPA and RPA(D) do not treat double excitations well with a scaling of O(N5), the r2ph-TDA is recommended to study systems with both single and some low-lying double excitations with a moderate accuracy, when the excitation energies are the central concern (transition moments from TDA are less unsatisfactory). Beyond the excitation energy tests, we also develop expressions on excited state property evaluations that are at least suitable for

$\langle \hat{S}^{2}\rangle$
Ŝ2 calculations, which are very useful to distinguish different spin states of the excitations.

Funding support from the National Science Foundation (CHE-1362927) is greatly appreciated. Y.Y. has also been supported by the Paul Gross Fellowship, Duke University.

1.
D.
Bohm
and
D.
Pines
,
Phys. Rev.
82
,
625
(
1951
).
2.
D.
Pines
and
D.
Bohm
,
Phys. Rev.
85
,
338
(
1952
).
5.
D.
Thouless
and
J.
Valatin
,
Nucl. Phys.
31
,
211
(
1962
).
6.
D. J.
Rowe
,
Nuclear Collective Motion: Models and Theory
(
Methuen
,
London
,
1970
).
7.
J.
Blaizot
and
G.
Ripka
,
Quantum Theory of Finite Systems
(
The MIT Press
,
Cambridge, MA
,
1986
).
8.
P.
Ring
and
P.
Schuck
,
The Nuclear Many-Body Problem
(
Springer
,
2004
).
9.
S.
Kurth
and
J. P.
Perdew
,
Phys. Rev. B
59
,
10461
(
1999
).
10.
11.
F.
Furche
,
J. Chem. Phys.
129
,
114105
(
2008
).
12.
H.
Eshuis
,
J.
Yarkony
, and
F.
Furche
,
J. Chem. Phys.
132
,
234114
(
2010
).
13.
D. C.
Langreth
and
J. P.
Perdew
,
Solid State Commun.
17
,
1425
(
1975
).
14.
A. G.
Eguiluz
,
Phys. Rev. Lett.
51
,
1907
(
1983
).
15.
G. F.
Giuliani
and
G.
Vignale
,
Quantum Theory Of The Electron Liquid
(
Cambridge University Press
,
2005
).
16.
J.
Harl
,
L.
Schimka
, and
G.
Kresse
,
Phys. Rev. B
81
,
115126
(
2010
).
17.
X.
Ren
,
P.
Rinke
,
C.
Joas
, and
M.
Scheffler
,
J. Mater. Sci.
47
,
7447
(
2012
).
18.
M. E.
Casida
, “
Time-dependent density functional response theory for molecules
,” in
Recent Advances in Computational Chemistry
, edited by
D. P.
Chong
(
World Scientific
,
Singapore
,
1995
), Vol.
1
, p.
155
.
19.
G. E.
Scuseria
,
T. M.
Henderson
, and
D. C.
Sorensen
,
J. Chem. Phys.
129
,
231101
(
2008
).
20.
M.
Fuchs
,
Y.-M.
Niquet
,
X.
Gonze
, and
K.
Burke
,
J. Chem. Phys.
122
,
094116
(
2005
).
21.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
,
Phys. Rev. A
85
,
042507
(
2012
).
22.
N.
Fukuda
,
F.
Iwamoto
, and
K.
Sawada
,
Phys. Rev.
135
,
A932
(
1964
).
23.
J.
Toivanen
and
J.
Suhonen
,
Phys. Rev. Lett.
75
,
410
(
1995
).
24.
W. J.
Mulhall
,
R. J.
Liotta
,
J. A.
Evans
, and
R. P.
Perazzo
,
Nucl. Phys. A
93
,
261
(
1967
).
25.
26.
D. J.
Rowe
,
Rev. Mod. Phys.
40
,
153
(
1968
).
27.
G.
Ripka
and
R.
Padjen
,
Nucl. Phys. A
132
,
489
(
1969
).
28.
K. A.
Brueckner
and
C. A.
Levinson
,
Phys. Rev.
97
,
1344
(
1955
).
29.
R.
Eden
,
Proc. R. Soc. London, Ser. A
235
,
408
(
1956
).
31.
D.
Peng
,
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
J. Chem. Phys.
140
,
18A522
(
2014
).
32.
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
Phys. Rev. A
88
,
030501
(
2013
).
33.
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
J. Chem. Phys.
140
,
18A511
(
2014
).
34.
J.
Čížek
,
J. Chem. Phys.
45
,
4256
(
1966
).
35.
D.
Peng
,
S. N.
Steinmann
,
H.
van Aggelen
, and
W.
Yang
,
J. Chem. Phys.
139
,
104112
(
2013
).
36.
G. E.
Scuseria
,
T. M.
Henderson
, and
I. W.
Bulik
,
J. Chem. Phys.
139
,
104113
(
2013
).
37.
P.
Mori-Sánchez
,
A.
Cohen
, and
W.
Yang
,
Phys. Rev. Lett.
102
,
066403
(
2009
).
38.
Y.
Yang
,
H.
van Aggelen
,
S. N.
Steinmann
,
D.
Peng
, and
W.
Yang
,
J. Chem. Phys.
139
,
174110
(
2013
).
39.
Y.
Yang
,
H.
van Aggelen
, and
W.
Yang
,
J. Chem. Phys.
139
,
224105
(
2013
).
40.
Y.
Yang
,
D.
Peng
,
J.
Lu
, and
W.
Yang
,
J. Chem. Phys.
141
,
124104
(
2014
).
41.
T.
Tamura
and
T.
Udagawa
,
Nucl. Phys.
53
,
33
(
1964
).
42.
C.
Yannouleas
,
Phys. Rev. C
35
,
1159
(
1987
).
43.
G.
Lauritsch
and
P.-G.
Reinhard
,
Nucl. Phys. A
509
,
287
(
1990
).
44.
D.
Gambacurta
and
F.
Catara
,
J. Phys.: Conf. Ser.
168
,
012012
(
2009
).
45.
T.-I.
Shibuya
and
V.
McKoy
,
Phys. Rev. A
2
,
2208
(
1970
).
46.
J.
Rose
,
T.-i.
Shibuya
, and
V.
McKoy
,
J. Chem. Phys.
58
,
74
(
1973
).
47.
P.
Jørgensen
,
J.
Oddershede
, and
M. A.
Ratner
,
Chem. Phys. Lett.
32
,
111
(
1975
).
48.
J.
Oddershede
,
Adv. Quantum Chem.
11
,
275
(
1978
).
49.
O.
Christiansen
,
K.
Bak
,
H.
Koch
, and
S.
Sauer
,
Chem. Phys. Lett.
284
,
47
(
1998
).
50.
M. E.
Casida
,
J. Chem. Phys.
122
,
054111
(
2005
).
51.
D.
Zhang
,
S. N.
Steinmann
, and
W.
Yang
,
J. Chem. Phys.
139
,
154109
(
2013
).
52.
K. W.
Sattelmeyer
,
H. F.
Schaefer
 III
, and
J. F.
Stanton
,
Chem. Phys. Lett.
378
,
42
(
2003
).
53.
J.
Shen
and
P.
Piecuch
,
J. Chem. Phys.
138
,
194102
(
2013
).
54.
See supplementary material at http://dx.doi.org/10.1063/1.4901716 for detailed derivations, size-extensivity discussion, molecular geometries, and MRCI calculation details.
55.
R.
McWeeny
,
Methods of Molecular Quantum Mechanics
,
Theoretical Chemistry
, 2nd ed. (
Academic Press
,
London
,
1989
), Vol.
2
.
56.
I.
Shavitt
and
R. J.
Bartlett
,
Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory
(
Cambridge University Press
,
New York
,
2009
).
57.
C.
Yannouleas
,
M.
Dworzecka
, and
J.
Griffin
,
Nucl. Phys. A
397
,
239
(
1983
).
58.
A.
Szabo
and
N. S.
Ostlund
,
Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory
(
Dover Publications
,
New York
,
1996
).
59.
J.
Čížek
and
J.
Paldus
,
J. Chem. Phys.
47
,
3976
(
1967
).
60.
D.
Yeager
,
M.
Nascimento
, and
V.
McKoy
,
Phys. Rev. A
11
,
1168
(
1975
).
61.
D. L.
Yeager
and
V.
McKoy
,
J. Chem. Phys.
67
,
2473
(
1977
).
62.
D.
Lynch
,
M. F.
Herman
, and
D. L.
Yeager
,
Chem. Phys.
64
,
69
(
1982
).
63.
A.
Ipatov
,
F.
Cordova
,
L. J.
Doriol
, and
M. E.
Casida
,
J. Mol. Struct.: THEOCHEM
914
,
60
(
2009
).
64.
E. R.
Davidson
,
J. Comput. Phys.
17
,
87
(
1975
).
65.
An in-house program for QM/MM simulations (http://www.qm4d.info).
66.
M. W.
Schmidt
 et al.,
J. Comput. Chem.
14
,
1347
(
1993
).
67.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
,
M.
Schütz
 et al., molpro, version 2012.1, a package of ab initio programs,
2012
, see http://www.molpro.net.
68.
H.-J.
Werner
,
P. J.
Knowles
,
G.
Knizia
,
F. R.
Manby
, and
M.
Schütz
,
WIREs Comput. Mol. Sci.
2
,
242
(
2012
).
69.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
, et al., Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford, CT,
2009
.
70.
L. A.
Curtiss
,
K.
Raghavachari
,
P. C.
Redfern
, and
J. A.
Pople
,
J. Chem. Phys.
106
,
1063
(
1997
).
71.
C. T.
Lee
,
W. T.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
72.
A. D.
Becke
,
J. Chem. Phys.
98
,
5648
(
1993
).
73.
K.
Aidas
 et al.,
WIREs: Comput. Mol. Sci.
4
,
269
(
2014
).
74.
Dalton, a molecular electronic structure program, release dalton2013.4 (
2013
), see http://daltonprogram.org.
75.
Z.-L.
Cai
,
D. J.
Tozer
, and
J. R.
Reimers
,
J. Chem. Phys.
113
,
7084
(
2000
).
76.
A.
Kramida
,
Yu.
Ralchenko
,
J.
Reader
, and
NIST ASD Team
, NIST Atomic Spectra Database (ver. 5.1), [Online]. Available: http://physics.nist.gov/asd [2014, April 17]. National Institute of Standards and Technology, Gaithersburg, MD,
2013
.
77.
H.
Koch
,
O.
Christiansen
,
P.
Jørgensen
, and
J.
Olsen
,
Chem. Phys. Lett.
244
,
75
(
1995
).
79.
D. J.
Tozer
and
N. C.
Handy
,
J. Chem. Phys.
109
,
10180
(
1998
).
80.
I.
Levine
,
Quantum Chemistry
, 6th ed. (
Prentice-Hall
,
2008
).
81.
S. B.
Ben-Shlomo
and
U.
Kaldor
,
J. Chem. Phys.
92
,
3680
(
1990
).
82.
M.
Schreiber
,
M. R.
Silva-Junior
,
S. P. A.
Sauer
, and
W.
Thiel
,
J. Chem. Phys.
128
,
134110
(
2008
).
83.
K.
Takayanagi
,
K.
Shimizu
, and
A.
Arima
,
Nucl. Phys. A
477
,
205
(
1988
).
84.
Note that in Ref. 6 the normalization factor
$\langle 0|[{\hat{O}}_{F},{\hat{O}}_{F}^{\dagger }]0\rangle$
0|[ÔF,ÔF]0
is absent, probably because they study the excited state but not the deexcited state and the normalization is 1 by default.

Supplementary Material