The particle-particle random phase approximation (pp-RPA) has been used to investigate excitation problems in our recent paper [Y. Yang, H. van Aggelen, and W. Yang, J. Chem. Phys.139, 224105 (2013)]. It has been shown to be capable of describing double, Rydberg, and charge transfer excitations, which are challenging for conventional time-dependent density functional theory (TDDFT). However, its performance on larger molecules is unknown as a result of its expensive O(N6) scaling. In this article, we derive and implement a Davidson iterative algorithm for the pp-RPA to calculate the lowest few excitations for large systems. The formal scaling is reduced to O(N4), which is comparable with the commonly used configuration interaction singles (CIS) and TDDFT methods. With this iterative algorithm, we carried out benchmark tests on molecules that are significantly larger than the molecules in our previous paper with a reasonably large basis set. Despite some self-consistent field convergence problems with ground state calculations of (N − 2)-electron systems, we are able to accurately capture lowest few excitations for systems with converged calculations. Compared to CIS and TDDFT, there is no systematic bias for the pp-RPA with the mean signed error close to zero. The mean absolute error of pp-RPA with B3LYP or PBE references is similar to that of TDDFT, which suggests that the pp-RPA is a comparable method to TDDFT for large molecules. Moreover, excitations with relatively large non-HOMO excitation contributions are also well described in terms of excitation energies, as long as there is also a relatively large HOMO excitation contribution. These findings, in conjunction with the capability of pp-RPA for describing challenging excitations shown earlier, further demonstrate the potential of pp-RPA as a reliable and general method to describe excitations, and to be a good alternative to TDDFT methods.

To predict electronic excited states with high accuracy and efficiency has been a great desire for theoretical chemists for a long time. Many theoretical approaches1–6 have been developed over the past few decades. Among these approaches, time-dependent density functional theory (TDDFT)4,7 has been widely used in large systems because of its low computational cost and satisfying accuracy. However, its inability to well describe double, Rydberg, charge transfer, and extended π-systems excitations6 within its formal linear-response formulation has seen its limitation in many practical applications.

Recently, the particle-particle (pp-) random phase approximation (RPA),8–10 a counterpart of the more well-known particle-hole (ph-) RPA,11–15 has been adopted from nuclear physics16,17 for applications in electronic systems. Furthermore, the pp-RPA correlation energy can be viewed as a simplest approximation within the exact formulation of electronic correlation energy in terms of the adiabatic-connection and the pairing matrix fluctuation.8 The discussion of pp-RPA has been extended to many more aspects, including: (1) the pp-RPA’s equivalence to ladder-coupled-cluster doubles (lCCD),9,10 (2) extensive benchmark tests on enthalpies of formation, reaction energies, reaction barriers, and van der Waals interactions,18,19 (3) the ladder channel's performance on homogeneous electron gas,20,21 (4) the cost reduction from O(N6) to O(N4) within the lCCD formulation by using the tensor-hyper-contraction technique.22 These above discussions are all within the field of correlation energy calculations. Meanwhile, we developed an approach to calculate excitation energies based on the pp-RPA.23 It has been tested on a small number of atoms and molecules to illustrate its ability to describe double, Rydberg, and charge transfer excitations with a proper choice of DFT/Hartree-Fock (HF) references. We also carried out calculations with particle-particle Tamm-Dancoff approximation (pp-TDA) and the results are very similar to pp-RPA for those small systems. Moreover, like the linear response TDDFT, we formulate our approach24 by investigating the linear response of systems to pairing fields, which we name it as the TDDFT-P theory. We conclude that our approach of calculating excitation energies with pp-RPA is a special case for the TDDFT-P theory with a bare Coulomb kernel, without considering the contributions of the pairing matrix to the exchange correlation energy. The TDDFT-P formulation requires the use of the (generalized) Kohn-Sham DFT reference in the formulation and therefore naturally justifies the use of approximate DFT references in practically applications of pp-RPA.24 The relation between pp-RPA and pp-TDA and the approaches to derive them can be summarized in Fig. 1. In addition, Zhang et al.25 have illustrated the calculations of gradient for excited states as well as ground states, which makes the geometry optimization for ground states and excited states applicable through pp-RPA. We also have extended our attempts to the second-ph-RPA/second-pp-RPA to include double/non-HOMO excitations,26 which are prohibited by their formal first-order formulations.

FIG. 1.

Schematic sketch of the relation between pp-RPA and pp-TDA and the ways to derive them.

FIG. 1.

Schematic sketch of the relation between pp-RPA and pp-TDA and the ways to derive them.

Close modal

Despite these efforts, we were still solving the pp-RPA equation using the direct diagonalization approach with a heavy computational cost O(N6), and therefore we are limited to tests on a small number of systems with few electrons and small basis sets. To further make this method applicable to larger systems, a computational algorithm with less computational cost is highly demanded. More importantly, since the performance of pp-RPA on large systems excitations is still not known, an extensive benchmark test is needed.

In this article, we first review the method of calculating excitations with pp-RPA followed by a Davidson-like approach to lower the computational cost to O(N4) for a single eigenroot in Sec. II. Then in Sec. III we carry out benchmark tests with molecules in Thiel's27,28 and Tozer's29 test sets and compare our results with the well-known low-cost methods including configuration interaction singles (CIS), time-dependent Hartree-Fock (TD-HF), and TDDFT with B3LYP(TD-B3LYP) or PBE(TD-PBE) references. Finally, in Sec. IV we give our concluding remarks.

The formal derivation for the TDDFT-P theory has been presented in Ref. 24. Here, we review some of the most important parts to give our readers a general impression of this fundamental theory.

We now consider a Hamiltonian within a pairing field,

\begin{equation}\hat{H}=\hat{T}+\hat{V}+\hat{W}+\hat{D},\end{equation}
Ĥ=T̂+V̂+Ŵ+D̂,
(1)

in which

$\hat{T}$
T̂⁠,
$\hat{V}$
V̂
, and
$\hat{W}$
Ŵ
represent kinetic energy, external potential, and two-electron interactions, respectively, and
$\hat{D}$
D̂
is the external pairing field,

\begin{equation}\hat{D}=\frac{1}{2}\int d\mathbf {x}d\mathbf {x^{\prime }}[D^*\hat{\psi }(\mathbf {x^{\prime }})\hat{\psi }(\mathbf {x})+\mathit {h.c.}], \end{equation}
D̂=12dxdx[D*ψ̂(x)ψ̂(x)+h.c.],
(2)

where

$\hat{\psi }(\mathbf {x^{\prime }})\hat{\psi }(\mathbf {x})$
ψ̂(x)ψ̂(x) stands for the pair removal part, while
$\mathit {h.c.}$
h.c.
is short for Hermitian conjugate and it stands for the pair addition part. In presence of the above pairing field, the pairing matrix

\begin{equation}\kappa (\mathbf {x},\mathbf {x^{\prime }})=\langle \Psi |\psi ^\dagger (\mathbf {x^{\prime }})\psi ^\dagger (\mathbf {x})|\Psi \rangle \end{equation}
κ(x,x)=Ψ|ψ(x)ψ(x)|Ψ
(3)

is not zero. In fact, a perturbative pairing field δD(y, y′; τ) at time τ results in a tiny change of the pairing matrix δκ(x, x′; t) at time t. If we ignore the higher order terms and only investigate the first order change, then through the linear response theory, we obtain a linear pp-response function

\begin{eqnarray}&&K(\mathbf {x},\mathbf {x^{\prime }};\mathbf {y},\mathbf {y^{\prime }};t)\nonumber\\&&\quad =-i\theta (t)\langle \Psi |[\psi (\mathbf {x^{\prime }};t)\psi (\mathbf {x};t),\psi ^\dagger (\mathbf {y};0)\psi ^\dagger (\mathbf {y^{\prime }};0)]|\Psi \rangle\qquad\end{eqnarray}
K(x,x;y,y;t)=iθ(t)Ψ|[ψ(x;t)ψ(x;t),ψ(y;0)ψ(y;0)]|Ψ
(4)

and the linear response equation with the integrated form is

\begin{equation}\delta \kappa (\mathbf {x},\mathbf {x^{\prime }};t)= \int d\tau d\mathbf {y}d\mathbf {y^{\prime }} K(\mathbf {x},\mathbf {x^{\prime }};\mathbf {y},\mathbf {y^{\prime }};t-\tau ) \delta D(\mathbf {y},\mathbf {y^{\prime }};\tau ). \end{equation}
δκ(x,x;t)=dτdydyK(x,x;y,y;tτ)δD(y,y;τ).
(5)

Note here, the creation and annihilation operators in Eq. (4) are all in Heisenberg picture.

A non-interacting system |Φs⟩ is now assumed to hold the same electron density and pairing matrix as the interacting real system at every time. Instead of complicated two-electron interactions, this non-interacting system only has effective one-body normal and pairing potentials, and the exchange-correlation part exists in both potentials. The total pairing field includes both the internal pairing potential and the external pairing field: δDs = δDint + δD. This non-interacting system is a mapping image from the real many-body system and it has a pp-response function that is much easier to calculate. By changing the coordinate basis to orbital basis and performing Fourier transform, its pp-response function K0 is simply

\begin{eqnarray}K^{0}_{pq,rs}(\omega )&=&(\delta _{pr}\delta _{qs}-\delta _{qr}\delta _{ps})\nonumber\\&&\times\,\frac{\theta (p-F)\theta (q-F)-\theta (F-p)\theta (F-q)}{\omega -(\epsilon _p+\epsilon _q)+i\eta},\nonumber\\[-4pt]\end{eqnarray}
Kpq,rs0(ω)=(δprδqsδqrδps)×θ(pF)θ(qF)θ(Fp)θ(Fq)ω(εp+εq)+iη,
(6)

with F standing for the Fermi level. Therefore, the linear response equation for the non-interacting system becomes two separate ones,

\begin{align}\delta \kappa ^s_{ij}(\omega )&=-\frac{\delta D^s_{ij}(\omega )}{\omega -(\epsilon _i+\epsilon _j)+i\eta },\end{align}
δκijs(ω)=δDijs(ω)ω(εi+εj)+iη,
(7a)
\begin{align}\delta \kappa ^s_{ab}(\omega )&=\frac{\delta D^s_{ab}(\omega )}{\omega -(\epsilon _a+\epsilon _b)+i\eta },\end{align}
δκabs(ω)=δDabs(ω)ω(εa+εb)+iη,
(7b)
where i, j stand for occupied orbitals while a, b for virtual orbitals.

Apart from above response equations that build up the dependence of δκ on the total pairing field δDs, the changes of pairing matrix also in turn affect the total pairing field, or more specifically, the internal mean-field pairing potential δDint. We use a response kernel L to represent the dependence of δDint on δκ and it can be derived that

\begin{eqnarray}L_{pq,rs}&=&\langle pq||rs\rangle +2\int d\mathbf {x}_1d\mathbf {x}_2d\mathbf {x}^{\prime }_1d\mathbf {x}^{\prime }_2\phi ^*_p(\mathbf {x}_1)\phi ^*_q(\mathbf {x}^{\prime }_1) \nonumber\\&&\times\,\left(\frac{\delta ^2E_{xc}[\rho ,\kappa ]}{\delta \kappa ^*(\mathbf {x}_1,\mathbf {x}^{\prime }_1)\delta \kappa (\mathbf {x}_2,\mathbf {x}^{\prime }_2)}\right)_{\rho } \phi _r(\mathbf {x}_2)\phi _s(\mathbf {x}^{\prime }_2),\end{eqnarray}
Lpq,rs=pqrs+2dx1dx2dx1dx2ϕp*(x1)ϕq*(x1)×δ2Exc[ρ,κ]δκ*(x1,x1)δκ(x2,x2)ρϕr(x2)ϕs(x2),
(8)

with

\begin{equation}\langle pq||rs\rangle =\langle pq|rs\rangle -\langle pq|sr\rangle \end{equation}
pqrs=pq|rspq|sr
(9)

and

\begin{equation}\langle pq|rs\rangle =\int d\mathbf {x}d\mathbf {x}^{\prime }\frac{\phi ^*_p(\mathbf {x})\phi ^*_q(\mathbf {x}^{\prime })\phi _r(\mathbf {x})\phi _s(\mathbf {x}^{\prime })}{|\mathbf {r}-\mathbf {r}^{\prime }|}. \end{equation}
pq|rs=dxdxϕp*(x)ϕq*(x)ϕr(x)ϕs(x)|rr|.
(10)

Plug in the response equation and eliminate the internal pairing potential, after some rearrangement, we come to the TDDFT-P equation

\begin{equation}{\left[\begin{array}{c@{\quad}c}\bf A & \bf B\\\bf B^\dagger & \bf C\\\end{array}\right]} {\left[\begin{array}{c}\mathbf {X}\\\mathbf {Y} \end{array}\right]} - \omega {\left[\begin{array}{c@{\quad}c}\bf I & \bf 0\\\bf 0 & \bf -I\\\end{array}\right]} {\left[\begin{array}{c}\mathbf {X}\\\mathbf {Y} \end{array}\right]} = - {\left[\begin{array}{c}\delta \mathbf {D}^{pp}\\\delta \mathbf {D}^{hh} \end{array}\right]}, \end{equation}
ABBCXYωI00IXY=δDppδDhh,
(11)

in which

\begin{equation}A_{ab,cd}=\delta _{ac}\delta _{bd}(\epsilon _a+\epsilon _b)+L_{ab,cd},\end{equation}
Aab,cd=δacδbd(εa+εb)+Lab,cd,
(12a)
\begin{equation}\hspace*{-30pt}B_{ab,ij}=L_{ab,ij},\hspace*{43pt}\end{equation}
Bab,ij=Lab,ij,
(12b)
\begin{equation}\hspace*{2pt}C_{ij,kl}=-\delta _{ik}\delta _{jl}(\epsilon _i+\epsilon _j)+L_{ij,kl},\hspace*{-2pt}\end{equation}
Cij,kl=δikδjl(εi+εj)+Lij,kl,
(12c)
\begin{equation}\hspace*{-30pt}X_{ab}=\delta \kappa _{ab}(\omega ),\hspace*{27pt}\end{equation}
Xab=δκab(ω),
(12d)
\begin{equation}\hspace*{-30pt}Y_{ij}=\delta \kappa _{ij}(\omega ),\hspace*{27pt}\end{equation}
Yij=δκij(ω),
(12e)
\begin{equation}\hspace*{-40pt}[\delta \mathbf {D}^{pp}]_{ab}=\delta D_{ab}(\omega ),\hspace*{35pt}\end{equation}
[δDpp]ab=δDab(ω),
(12f)
\begin{equation}\hspace*{-38pt}[\delta \mathbf {D}^{hh}]_{ij}=\delta D_{ij}(\omega ).\hspace*{36pt}\end{equation}
[δDhh]ij=δDij(ω).
(12g)
Note here, we have further simplified the equation by restricting a > b, c > d, i > j, and k > l. In real atomic or molecular systems, we take the limit that the external pairing field D goes to 0. Therefore, setting the right-hand side of Eq. (11) to be zero, we obtain a generalized eigenvalue equation

\begin{equation}{\left[\begin{array}{c@{\quad}c}\bf A & \bf B\\\bf B^\dagger & \bf C\\\end{array}\right]} {\left[\begin{array}{c}\mathbf {X}\\\mathbf {Y} \end{array}\right]} = \omega {\left[\begin{array}{c@{\quad}c}\bf I & \bf 0\\\bf 0 & \bf -I\\\end{array}\right]} {\left[\begin{array}{c}\mathbf {X}\\\mathbf {Y} \end{array}\right]}. \end{equation}
ABBCXY=ωI00IXY.
(13)

Since the response kernel L is still not well known, if we simply ignore the exchange-correlation part in L in Eq. (8), we arrive at the pp-RPA equation, with matrix elements being

\begin{align}A_{ab,cd}&=\delta _{ac}\delta _{bd}(\epsilon _a+\epsilon _b)+\langle ab||cd\rangle ,\end{align}
Aab,cd=δacδbd(εa+εb)+abcd,
(14a)
\begin{equation}\hspace*{-34pt}B_{ab,ij}=\langle ab||ij\rangle ,\hspace*{42pt}\end{equation}
Bab,ij=abij,
(14b)
\begin{equation}C_{ij,kl}=-\delta _{ik}\delta _{jl}(\epsilon _i+\epsilon _j)+\langle ij||kl\rangle .\end{equation}
Cij,kl=δikδjl(εi+εj)+ijkl.
(14c)

The dimension of the square matrices A and C in Eq. (13) are

$\frac{N_{vir}(N_{vir}-1)}{2}$
Nvir(Nvir1)2 and
$\frac{N_{occ}(N_{occ}-1)}{2}$
Nocc(Nocc1)2
, respectively, and the dimension of matrix B is
$\frac{N_{vir}(N_{vir}-1)}{2}\times \frac{N_{occ}(N_{occ}-1)}{2}$
Nvir(Nvir1)2×Nocc(Nocc1)2
. There are totally
$\frac{N_{vir}(N_{vir}-1)+ N_{occ}(N_{occ}-1)}{2}$
Nvir(Nvir1)+Nocc(Nocc1)2
real eigenroots. Among them,
$\frac{N_{vir}(N_{vir}-1)}{2}$
Nvir(Nvir1)2
eigenvectors normalize to 1, while the rest
$\frac{N_{occ}(N_{occ}-1)}{2}$
Nocc(Nocc1)2
normalize to −1

\begin{equation}X^T_IX_I-Y^T_IY_I=\pm 1. \end{equation}
XITXIYITYI=±1.
(15)

If a real eigenvector normalizes to 1, it is for a pair addition process and if normalizes to −1, it is for a pair removal process. When we calculate excitations with pp-RPA, we usually start with a (N − 2)-electron system and then add two electrons back. Therefore, the eigenvalues of our interests are those lowest ones with normalization to 1.

The time cost for a direct diagonalization approach to solve the pp-RPA equation (Eq. (13)) is O(N6) and the memory space cost is O(N4), with N being the maximum of Nvir and Nocc. For problems with small number of electrons and also small basis sets, this direct diagonalization approach is applicable. However, for larger systems, the computational cost significantly increases with more electron numbers and basis functions.

We now use the basic idea behind the Davidson method30 to solve the lowest few pair-addition eigenroots for the pp-RPA equation with an O(N4) time cost.

We first simplify the notation and write the pp-RPA equation (13) as

\begin{equation}\mathbf {Mu=\omega W u}, \end{equation}
Mu=ωWu,
(16)

where M is the pp-RPA matrix [A, B; B†, C] and W is the diagonal matrix [I, 0; 0, −I]. Suppose for an exact eigenpair (ωk,uk), we approximate it by (

$\tilde{\omega }_k$
ω̃k⁠,
$\tilde{\mathbf {u}}_k$
ũk
) such that
$\tilde{\mathbf {u}}_k$
ũk
is a linear combination of v1, v2, …, vn

\begin{equation}\mathbf {u}_k\approx \mathbf {\tilde{u}}_k= {\left[\begin{array}{c}\mathbf {v}_1 \mathbf {v}_2 \cdots \mathbf {v}_n \end{array}\right]} {\left[\begin{array}{c}c_1\\c_2\\\vdots \\c_n \end{array}\right]} =\mathbf {V}\mathbf {c}, \end{equation}
ukũk=v1v2vnc1c2cn=Vc,
(17)

where v's are basis vectors and they are orthonormalized with respect to W,

\begin{equation}{[\mathbf {V}^T\mathbf {WV}]}_{ij}=\delta _{ij}\operatorname{sgn}\big(\mathbf {v}_i^T\mathbf {W}\mathbf {v}_i\big), \end{equation}
[VTWV]ij=δijsgnviTWvi,
(18)

and c's are linear combination coefficients. Therefore, the pp-RPA eigenvalue equation can be approximated by

\begin{equation}\mathbf {MVc}=\tilde{\omega }_{k} \mathbf {W Vc}. \end{equation}
MVc=ω̃kWVc.
(19)

Multiply VT to the left, we obtain

\begin{equation}\tilde{\mathbf {M}}\mathbf {c}=\tilde{\omega }_{k}(\mathbf {V}^T\mathbf {WV})\mathbf {c}, \end{equation}
M̃c=ω̃k(VTWV)c,
(20)

with

\begin{equation}\tilde{\mathbf {M}}=\mathbf {V}^T\mathbf {MV}. \end{equation}
M̃=VTMV.
(21)

The matrix

$\tilde{\mathbf {M}}$
M̃ in Eq. (21) has dimension n, which is the number of basis vectors and it is usually much smaller than the original matrix M. By solving the eigenvalue problem (Eq. (20)), we are able to get the approximated eigenvalue
$\tilde{\omega }_k$
ω̃k
and the coefficients c, and therefore the approximated eigenvector
$\tilde{\mathbf {u}}_k$
ũk
through Eq. (27). To test whether
$\tilde{\mathbf {u}}_k$
ũk
is a good approximation, we calculate the residual rk,

\begin{equation}\mathbf {r}_k=\mathbf {M}\mathbf {\tilde{u}}_k-\tilde{\omega }_k \mathbf {W} \mathbf {\tilde{u}}_k. \end{equation}
rk=Mũkω̃kWũk.
(22)

If the norm of rk is within a given threshold, we consider the result to be converged, otherwise, we need to expand the basis vector space to help obtain a better approximation.

Suppose the difference between uk and

$\tilde{\mathbf {u}}_k$
ũk is

\begin{equation}\mathbf {t}=\mathbf {u}_k-\tilde{\mathbf {u}}_k. \end{equation}
t=ukũk.
(23)

Then the original eigenvalue problem can be written as

\begin{eqnarray}\mathbf {M}(\mathbf {\tilde{u}}_k+\mathbf {t})&=&\omega _k \mathbf {W}(\mathbf {\tilde{u}}_k+\mathbf {t}),\nonumber\\[-6pt]\\[-6pt](\mathbf {M}-\omega _k \mathbf {W})\mathbf {t}&=&-(\mathbf {M}\mathbf {\tilde{u}}_k-\omega _k \mathbf {W} \mathbf {\tilde{u}}_k).\nonumber \end{eqnarray}
M(ũk+t)=ωkW(ũk+t),(MωkW)t=(MũkωkWũk).
(24)

Assume the approximated eigenvalue

$\tilde{\omega }_k$
ω̃k is already a good approximation to the real one ωk, we would then need to get t by solving (recall the definition of rk in Eq. (22))

\begin{equation}(\mathbf {M}-\tilde{\omega }_k \mathbf {W})\mathbf {t}=-\mathbf {r}_k.\end{equation}
(Mω̃kW)t=rk.
(25)

To strictly solve the difference vector t, we need to calculate the inverse of

$\mathbf {M}-\tilde{\omega }_k \mathbf {W}$
Mω̃kW⁠, but it is expensive to do so. Fortunately, following Davidson's suggestion, a good preconditioner P can be constructed using the diagonal part of M and the diagonal metric matrix W, since usually, the orbital energy parts in A and C are much larger than the other matrix elements,

\begin{equation}\mathbf {P}_{ij}=\delta _{ij}(\mathbf {M}_{ii}-\tilde{\omega }_k \mathbf {W}_{ii}). \end{equation}
Pij=δij(Miiω̃kWii).
(26)

The inverse for P is easy to calculate and thus the ith element for the approximated difference vector t is

\begin{eqnarray}\mathbf {t}_i\approx \Delta \mathbf {u}_i&=&-[\mathbf {P}^{-1}\mathbf {r}_k]_i\nonumber\\&=& -\frac{[\mathbf {r}_k]_i}{\mathbf {M}_{ii}-\tilde{\omega }_k \mathbf {W}_{ii}}.\end{eqnarray}
tiΔui=[P1rk]i=[rk]iMiiω̃kWii.
(27)

As a good approximation to t, the newly calculated vector Δu can be added to the existing basis vector space V to help further obtain a better approximation to the exact eigenroot pair until finally converged.

It is worth noting that the way to augment such a basis vector space is not limited to the Davidson approach only. Other methods to expand the subspace, such as the Jacobi-Davidson approach,31,32 are also applicable. However, because our goal in this work is simply to lower the computational cost to run benchmark tests rather than to carry out a method comparison, we only implemented the original Davidson flavor of this subspace expansion.

Now let us present our algorithm with a detailed work flow.

  1. Perform a self-consistent field (SCF) calculation for the (N − 2)-electron system with HF or a chosen DFT functional.

  2. Generate an initial guess for the basis vector set. Because we aim for the lowest n pair addition eigenroots, a good initial guess can be generated by sorting the sum of any two virtual orbital energies ɛab = ɛa + ɛb and getting the lowest m ones (mn). Suppose the lth lowest value is ɛcd, then the lth initial basis vector elements are
    \begin{equation}{[\mathbf {v}_{l}]}_{pq}=\delta _{cp}\delta _{dq}. \end{equation}
    [vl]pq=δcpδdq.
    (28)
  3. Calculate matrix-vector product MV using the approach in Ref. 33. This is one of the most expensive steps in the whole calculation: for each basis vector, the time cost is O(N4) and memory cost is O(N2).

  4. Calculate vector-vector product

    $\mathbf {\tilde{M}}=\mathbf {V}^T(\mathbf {MV})$
    M̃=VT(MV)⁠. For this step, the time and memory costs are both O(N2).

  5. Solve the reduced generalized eigenvalue problem (Eq. (20)) to obtain approximated eigenvalues

    $\tilde{\omega }$
    ω̃ and coefficients c. Because the number of basis vectors is much smaller than N, the cost in this step is negligible.

  6. Sort the approximated eigenpairs (

    $\tilde{\omega }$
    ω̃⁠,Vc) and pick up the lowest n pair addition ones according to the normalization constraint Eq. (15).

  7. Calculate residual vectors using Eq. (22). Note here, instead of calculating

    $\mathbf {M}\mathbf {\tilde{u}}_k=\mathbf {M}(\mathbf {Vc})$
    Mũk=M(Vc)⁠, we calculate (MV)c, in which MV is already calculated in step 3. Therefore, the cost is also negligible.

  8. Calculate the norm of the residual vectors, if all converged, exit. Otherwise, for non-converged roots, calculate the approximated difference vectors Δu using Eq. (27).

  9. Orthogonalize Δu with respect to V using the generalized Gram-Schmidt method. If it is not numerically zero, calculate the self-product with metric W. Normalize it and add the to basis space V.

  10. Go back to step 3 and continue the loop.

The pp-TDA is simply solving the equation

\begin{equation}\mathbf {A}\mathbf {u}=\mathbf {\omega }\mathbf {u}. \end{equation}
Au=ωu.
(29)

The matrix A is a diagonally dominated Hermitian matrix. Therefore, the pp-TDA equation can be solved using the canonical Davidson algorithm.30 

We implemented the above Davidson iterative method on the spin-separated and spin-adapted pp-RPA and pp-TDA equations18 in QM4D package.34 Then we use it to benchmark excitation energies calculated with pp-RPA and compare with TDDFT results calculated with Gaussian 09.35 We use molecules with the number of atoms larger than five in Thiel's27,28 and Tozer's29 test sets. If a molecule exists in both test sets, we use Thiel's geometry and reference values. Note that most of these test molecules do not contain particularly challenging cases for TDDFT calculations of charge transfer, double and Rydberg excitations, which we have shown pp-RPA can describe well.23 In other words, these tests are for the general cases where TDDFT performs well normally.

Limited by the convergence difficulty in the ground state calculations (especially calculations with the PBE functional) for some (N − 2)-electron systems, we present a relatively complete comparison and discussion with different functionals for 15 molecules with the aug-cc-pVDZ basis set. There are 24 singlet excitations and 19 triplet excitations in this set. For each excitation, we perform pp-RPA and pp-TDA calculations with HF, B3LYP, and PBE references. We also compare these results with the well-known computationally efficient methods including CIS, TD-HF, TD-B3LYP, and TD-PBE in order to further assess the performance for our approaches. For the rest molecules, because the PBE functional cannot converge their (N − 2)-electron systems and methods with HF references (CIS, TD-HF, pp-RPA-HF) always give unsatisfying results, we only present converged pp-RPA-B3LYP results and compare with TD-B3LYP calculations. Note here, all these excitations chosen are all HOMO excitations or excitations with HOMO excitation characters. HOMO excitations are usually the most important low-lying excitations that are of interest, therefore, although we will miss those pure non-HOMO excitations using the canonical pp-RPA treatment, in general, we only miss a few important low-lying excitations in the systems we tested.

Basis set convergence is tested along the cc-pVXZ series, X = D,T,Q, as well as the aug-cc-pVXZ series, X = D,T. The QM4D program uses Cartesian atomic orbitals and removes basis functions with angular momentum higher than “f.” We choose butadiene and furan as test molecules. For each molecule, two lowest singlet and two lowest triplet excitations are investigated.

The results are shown in Fig. 2. Because of a convergence failure, aug-cc-pVTZ results for HF reference are missing. It can be seen that for the two DFT references, the excitation energies decrease from cc-pVDZ to cc-pVQZ. The energies further decrease when it comes to aug-cc-pVDZ. The difference between the results of aug-cc-pVTZ and aug-cc-pVDZ are very small (≈0.02 eV). Therefore, we can consider that the excitation energy is already converged for the aug-cc-pVDZ basis set. Note here, even though cc-pVTZ and cc-pVQZ have more contracted Gaussian-type orbitals (CGTOs) than aug-cc-pVDZ, they do not reach the basis set convergence for excitation energies. The reason is that excited states are more diffuse than the ground state and in order to describe a balanced ground state and excited states well, adding diffuse functions is more crucial and more efficient than adding angular-momentum functions. For calculations with HF references, we can observe the similar convergence trend. However, it converges much slower than DFT references. With the aug-cc-pVTZ result missing, we cannot guarantee the aug-cc-pVDZ result is converged. Judging from the figure, it is very likely that more basis functions are needed considering its slow convergence behavior. Since the aug-cc-pVDZ basis is already converged for DFT references and considering the computational cost, we decide to use aug-cc-pVDZ basis in all the rest calculations, even though we cannot guarantee a convergence with the HF reference.

FIG. 2.

Basis set convergence test for pp-RPA. All calculations show the lowest two singlet and lowest two triplet excitations. Calculations with B3LYP and PBE references converge fast, with aug-cc-pVDZ showing converged results. Calculations with HF references converge much slower. The importance of diffuse function is also observed in the convergence behavior, with aug-cc-pVDZ basis providing lower excitation energies than cc-pVQZ. (a) Butadiene by pp-RPA-HF, (b) butadiene by pp-RPA-B3LYP, (c) butadiene by pp-RPA-PBE, (d) furan by pp-RPA-HF, (e) furan by pp-RPA-B3LYP, (f) furan by pp-RPA-PBE.

FIG. 2.

Basis set convergence test for pp-RPA. All calculations show the lowest two singlet and lowest two triplet excitations. Calculations with B3LYP and PBE references converge fast, with aug-cc-pVDZ showing converged results. Calculations with HF references converge much slower. The importance of diffuse function is also observed in the convergence behavior, with aug-cc-pVDZ basis providing lower excitation energies than cc-pVQZ. (a) Butadiene by pp-RPA-HF, (b) butadiene by pp-RPA-B3LYP, (c) butadiene by pp-RPA-PBE, (d) furan by pp-RPA-HF, (e) furan by pp-RPA-B3LYP, (f) furan by pp-RPA-PBE.

Close modal

The calculation results for the 15 molecules are shown in Table I. The pp-RPA with HF references (pp-RPA-HF) accidentally has a 0 eV mean signed error, while the error for (pp-RPA-B3LYP and pp-RPA-PBE) are 0.04 eV and −0.14 eV, respectively. As to the mean absolute error, the errors are 0.92 eV, 0.40 eV, and 0.38 eV, respectively. These numbers mean that for pp-RPA, the calculated results distribute around the benchmark values with very little systematic bias. However, when we look at the mean absolute error, pp-RPA-HF has such a large one that even doubles those of pp-RPA-B3LYP or pp-RPA-PBE. Therefore, pp-RPA-HF is not as accurate as pp-RPA-DFT to treat lowest few single excitations. Moreover, when we analyse data from pp-RPA-HF, it happens in many cases that there lie many excitations with incorrect symmetry below the targeted excitations. This also happens sometimes in TD-HF or CIS calculations, which might result from the problem with HF references. For pp-RPA-B3LYP and pp-RPA-PBE, the mean absolute errors are much smaller which means that they are much more accurate. The excitation spectrum is also clean without inserted wrong symmetry states and makes it easier to analyze for these DFT references.

Table I.

Vertical excitation energies (in eV) from pp-RPA, pp-TDA, CIS, TD-HF, and TDDFT. Numbers in the parenthesis indicate the error. TD-HF instabilities are characterized by imaginary excitation energy, but are denoted with negative numbers and thus having large negative errors.

MoleculeExciRefRPA-HFRPA-B3LYPRPA-PBETDA-HFTDA-B3LYPTDA-PBECISTD-HFTD-B3LYPTD-PBE
Ethene 3B1u 4.5 3.92(−0.58) 3.62(−0.88) 3.48(−1.02) 3.90(−0.60) 3.58(−0.92) 3.44(−1.06) 3.59(−0.91) 0.74(−3.76) 4.05(−0.45) 4.24(−0.26) 
Ethene 1B1u 7.8 6.26(−1.54) 8.45(0.65) 8.85(1.05) 6.24(−1.56) 8.45(0.65) 8.86(1.06) 7.71(−0.09) 7.36(−0.44) 7.38(−0.42) 7.39(−0.41) 
Butadiene 3Bu 3.2 3.22(0.02) 2.53(−0.67) 2.28(−0.92) 3.11(−0.09) 2.31(−0.89) 2.02(−1.18) 2.63(−0.57) −2.18(−5.38) 2.79(−0.41) 2.95(−0.25) 
Butadiene 3Ag 5.08 5.60(0.52) 6.07(0.99) 5.85(0.77) 5.49(0.41) 5.85(0.77) 5.59(0.51) 4.33(−0.75) 2.95(−2.13) 4.85(−0.23) 4.99(−0.09) 
Butadiene 1Bu 6.18 5.49(−0.69) 6.57(0.39) 6.51(0.33) 5.38(−0.80) 6.38(0.20) 6.28(0.10) 6.19(0.01) 5.90(−0.28) 5.56(−0.62) 5.44(−0.74) 
Butadiene 1Ag 6.55 5.92(−0.63) 6.44(−0.11) 6.11(−0.44) 5.82(−0.73) 6.30(−0.25) 5.95(−0.60) 7.39(0.84) 7.26(0.71) 6.49(−0.06) 6.11(−0.44) 
Hexatriene 3Bu 2.4 2.59(0.19) 1.88(−0.52) 1.61(−0.79) 2.49(0.09) 1.68(−0.72) 1.36(−1.04) 2.11(−0.29) −2.58(−4.98) 2.12(−0.28) 2.28(−0.12) 
Hexatriene 3Ag 4.15 5.20(1.05) 4.86(0.71) 4.47(0.32) 5.10(0.95) 4.64(0.49) 4.19(0.04) 3.57(−0.58) 1.56(−2.60) 3.94(−0.21) 4.04(−0.11) 
Hexatriene 1Ag 5.09 5.42(0.33) 4.99(−0.10) 4.49(−0.60) 5.34(0.25) 4.86(−0.23) 4.32(−0.77) 6.88(1.79) 6.71(1.62) 5.50(0.41) 5.02(−0.07) 
Hexatriene 1Bu 5.1 5.03(−0.07) 5.29(0.19) 5.13(0.03) 4.94(−0.16) 5.09(−0.01) 4.88(−0.22) 5.33(0.23) 5.05(−0.05) 4.60(−0.50) 4.44(−0.66) 
Octetraene 3Bu 2.2 2.15(−0.05) 1.49(−0.71) 1.21(−0.99) 2.07(−0.13) 1.31(−0.89) 0.99(−1.21) 1.79(−0.41) −2.74(−4.94) 1.71(−0.49) 1.87(−0.33) 
Octetraene 3Ag 3.55 4.75(1.20) 4.00(0.45) 3.57(0.02) 4.66(1.11) 3.80(0.25) 3.30(−0.25) 3.01(−0.54) −1.12(−4.67) 3.26(−0.29) 3.36(−0.19) 
Octetraene 1Ag 4.47 5.02(0.55) 4.08(−0.39) 3.53(−0.94) 4.95(0.48) 3.96(−0.51) 3.37(−1.10) 6.39(1.92) 6.32(1.85) 4.80(0.33) 4.17(−0.30) 
Octetraene 1Bu 4.66 4.58(−0.08) 4.50(−0.16) 4.29(−0.37) 4.50(−0.16) 4.31(−0.35) 4.05(−0.61) 4.73(0.07) 4.47(−0.19) 3.96(−0.70) 3.78(−0.88) 
Decapentaene 1Bu 4.27 4.69(0.42) 4.60(0.33) 4.40(0.13) 4.61(0.34) 4.40(0.13) 4.13(−0.14) 4.88(0.61) 4.63(0.36) 4.07(−0.20) 3.86(−0.41) 
Cyclopropene 3B2 4.34 4.22(−0.12) 4.08(−0.26) 3.91(−0.43) 4.14(−0.20) 3.89(−0.45) 3.66(−0.68) 3.46(−0.88) 0.67(−3.67) 3.70(−0.64) 3.79(−0.55) 
Cyclopropene 1B2 7.06 5.87(−1.19) 7.29(0.23) 7.31(0.25) 5.80(−1.26) 7.12(0.06) 8.15(1.09) 6.67(−0.39) 6.40(−0.66) 6.09(−0.97) 5.91(−1.15) 
Cyclopentadiene 3B2 3.25 3.12(−0.13) 2.66(−0.59) 2.53(−0.72) 3.03(−0.22) 2.48(−0.77) 2.31(−0.94) 2.49(−0.76) −2.06(−5.31) 2.74(−0.51) 2.91(−0.34) 
Cyclopentadiene |$\underline{{}^3\!A_1}$|A13̲ 5.09 5.45(0.36) 5.33(0.24) 5.07(−0.02) 5.36(0.27) 5.15(0.06) 4.86(−0.23) 4.29(−0.80) 2.94(−2.15) 4.75(−0.34) 4.87(−0.22) 
Cyclopentadiene 1B2 5.55 5.16(−0.39) 5.46(−0.09) 5.45(−0.10) 5.07(−0.48) 5.31(−0.24) 5.26(−0.29) 5.43(−0.12) 5.12(−0.43) 4.95(−0.60) 4.88(−0.67) 
Cyclopentadiene |$\underline{{}^1\!A_1}$|A11̲ 6.31 6.01(−0.30) 6.42(0.11) 6.16(−0.15) 5.94(−0.37) 6.33(0.02) 6.05(−0.26) 7.87(1.56) 7.78(1.47) 6.40(0.09) 5.78(−0.53) 
Norbornadiene 3A2 3.72 3.77(0.05) 3.68(−0.04) 3.58(−0.14) 3.57(−0.15) 3.36(−0.36) 3.20(−0.52) 2.87(−0.85) −1.48(−5.20) 3.10(−0.62) 3.16(−0.56) 
Norbornadiene |$\underline{{}^3\!B_2}$|B23̲ 4.16 4.49(0.33) 4.76(0.60) 4.51(0.35) 4.28(0.12) 4.41(0.25) 4.11(−0.05) 3.27(−0.89) −1.32(−5.48) 3.63(−0.53) 3.78(−0.38) 
Norbornadiene 1A2 5.34 3.95(−1.39) 5.36(0.02) 5.34(0.00) 3.75(−1.59) 5.05(−0.29) 4.99(−0.35) 5.53(0.19) 5.30(−0.04) 4.70(−0.64) 4.40(−0.94) 
Norbornadiene |$\underline{{}^1\!B_2}$|B21̲ 6.11 4.58(−1.53) 6.84(0.73) 6.99(0.88) 4.38(−1.73) 6.51(0.40) 6.64(0.53) 7.02(0.91) 6.81(0.70) 5.28(−0.83) 4.94(−1.17) 
Furan 3B2 4.17 3.82(−0.35) 3.49(−0.68) 3.37(−0.80) 3.74(−0.43) 3.33(−0.84) 3.17(−1.00) 3.28(−0.89) −1.14(−5.31) 3.70(−0.47) 3.89(−0.28) 
Furan |$\underline{{}^3\!A_1}$|A13̲ 5.48 5.79(0.31) 5.37(−0.11) 5.12(−0.36) 5.71(0.23) 5.22(−0.26) 4.93(−0.55) 4.87(−0.61) 4.07(−1.41) 5.18(−0.30) 5.25(−0.23) 
Furan 1B2 6.32 5.08(−1.24) 6.57(0.25) 6.58(0.26) 5.00(−1.32) 6.42(0.10) 6.40(0.08) 6.28(−0.04) 5.97(−0.35) 5.93(−0.39) 5.88(−0.44) 
Furan |$\underline{{}^1\!A_1}$|A11̲ 6.57 6.61(0.04) 6.88(0.31) 6.65(0.08) 6.54(−0.03) 6.80(0.23) 6.56(−0.01) 7.89(1.32) 7.77(1.20) 6.58(0.01) 6.28(−0.29) 
Pridazine 1B1 3.78 1.51(−2.27) 2.68(−1.10) 2.91(−0.87) 1.27(−2.51) 2.33(−1.45) 2.50(−1.28) 4.93(1.15) 4.72(0.94) 3.60(−0.18) 3.15(−0.63) 
Pridazine 1A2 4.32 2.91(−1.41) 3.41(−0.91) 3.48(−0.84) 2.60(−1.72) 2.90(−1.42) 2.91(−1.41) 6.13(1.81) 5.95(1.63) 4.20(−0.12) 3.53(−0.79) 
s-tetrazine 3B3u 1.89 3.28(1.39) 2.19(0.30) 1.86(−0.03) 2.94(1.05) 1.64(−0.25) 1.24(−0.65) 2.40(0.51) 1.89(0.00) 1.47(−0.42) 1.15(−0.74) 
s-tetrazine 3Au 3.52 5.38(1.86) 3.67(0.15) 3.16(−0.36) 4.91(1.39) 2.91(−0.61) 2.30(−1.22) 4.67(1.15) 4.38(0.86) 3.15(−0.37) 2.54(−0.98) 
s-tetrazine 1B3u 2.24 3.78(1.54) 2.73(0.49) 2.41(0.17) 3.53(1.29) 2.33(0.09) 1.96(−0.28) 3.54(1.30) 3.32(1.08) 2.27(0.03) 1.85(−0.39) 
s-tetrazine 1Au 3.48 5.59(2.11) 3.91(0.43) 3.40(−0.08) 5.13(1.65) 3.18(−0.30) 2.58(−0.90) 5.71(2.23) 5.55(2.07) 3.54(0.06) 2.87(−0.61) 
Formaldehyde 3A2 3.5 1.65(−1.85) 3.15(−0.35) 3.40(−0.10) 1.53(−1.97) 2.85(−0.65) 3.04(−0.46) 3.67(0.17) 3.34(−0.16) 3.10(−0.40) 2.97(−0.53) 
Formaldehyde 1A2 3.88 2.00(−1.88) 3.68(−0.20) 3.97(0.09) 1.88(−2.00) 3.43(−0.45) 3.68(−0.20) 4.49(0.61) 4.31(0.43) 3.83(−0.05) 3.71(−0.17) 
Acetone 3A2 4.05 3.10(−0.95) 4.13(0.08) 4.24(0.19) 2.92(−1.13) 3.78(−0.27) 3.83(−0.22) 4.42(0.37) 4.13(0.08) 3.68(−0.37) 3.52(−0.53) 
Acetone 1A2 4.4 3.38(−1.02) 4.56(0.16) 4.66(0.26) 3.20(−1.20) 4.25(−0.15) 4.31(−0.09) 5.15(0.75) 4.98(0.58) 4.31(−0.09) 4.14(−0.26) 
Benzoquinone 3B1g 2.51 4.74(2.23) 2.93(0.42) 2.44(−0.07) 4.07(1.56) 2.05(−0.46) 1.53(−0.98) 3.31(0.80) 3.01(0.50) 1.93(−0.58) 1.42(−1.09) 
Benzoquinone 1B1g 2.78 4.98(2.20) 3.14(0.36) 2.65(−0.13) 4.36(1.58) 2.36(−0.42) 1.85(−0.93) 4.00(1.22) 3.84(1.06) 2.43(−0.35) 1.87(−0.91) 
DMABN 1A1 4.56 6.14(1.58) 4.96(0.40) 4.56(0.00) 6.04(1.48) 4.70(0.14) 3.90(−0.66) 5.45(0.89) 5.16(0.60) 4.52(−0.04) 4.19(−0.37) 
DMABN 1B1 4.25 5.68(1.43) 4.68(0.43) 4.24(−0.01) 5.58(1.33) 4.42(0.17) 4.21(−0.04) 5.54(1.29) 5.35(1.10) 4.39(0.14) 3.97(−0.28) 
Total   MSE 0.00 0.04 −0.14 −0.16 −0.24 −0.44 0.31 −0.95 −0.32 −0.49 
Total   MAE 0.92 0.40 0.38 0.89 0.43 0.60 0.79 1.82 0.37 0.49 
Singlets   MSE −0.23 0.10 −0.04 −0.38 −0.16 −0.30 0.84 0.62 −0.24 −0.56 
Singlets   MAE 1.08 0.36 0.34 1.08 0.34 0.54 0.89 0.83 0.33 0.56 
Triplets   MSE 0.29 −0.05 −0.27 0.12 −0.34 −0.62 −0.36 −2.93 −0.42 −0.41 
Triplets   MAE 0.71 0.46 0.44 0.64 0.54 0.67 0.67 3.08 0.42 0.41 
MoleculeExciRefRPA-HFRPA-B3LYPRPA-PBETDA-HFTDA-B3LYPTDA-PBECISTD-HFTD-B3LYPTD-PBE
Ethene 3B1u 4.5 3.92(−0.58) 3.62(−0.88) 3.48(−1.02) 3.90(−0.60) 3.58(−0.92) 3.44(−1.06) 3.59(−0.91) 0.74(−3.76) 4.05(−0.45) 4.24(−0.26) 
Ethene 1B1u 7.8 6.26(−1.54) 8.45(0.65) 8.85(1.05) 6.24(−1.56) 8.45(0.65) 8.86(1.06) 7.71(−0.09) 7.36(−0.44) 7.38(−0.42) 7.39(−0.41) 
Butadiene 3Bu 3.2 3.22(0.02) 2.53(−0.67) 2.28(−0.92) 3.11(−0.09) 2.31(−0.89) 2.02(−1.18) 2.63(−0.57) −2.18(−5.38) 2.79(−0.41) 2.95(−0.25) 
Butadiene 3Ag 5.08 5.60(0.52) 6.07(0.99) 5.85(0.77) 5.49(0.41) 5.85(0.77) 5.59(0.51) 4.33(−0.75) 2.95(−2.13) 4.85(−0.23) 4.99(−0.09) 
Butadiene 1Bu 6.18 5.49(−0.69) 6.57(0.39) 6.51(0.33) 5.38(−0.80) 6.38(0.20) 6.28(0.10) 6.19(0.01) 5.90(−0.28) 5.56(−0.62) 5.44(−0.74) 
Butadiene 1Ag 6.55 5.92(−0.63) 6.44(−0.11) 6.11(−0.44) 5.82(−0.73) 6.30(−0.25) 5.95(−0.60) 7.39(0.84) 7.26(0.71) 6.49(−0.06) 6.11(−0.44) 
Hexatriene 3Bu 2.4 2.59(0.19) 1.88(−0.52) 1.61(−0.79) 2.49(0.09) 1.68(−0.72) 1.36(−1.04) 2.11(−0.29) −2.58(−4.98) 2.12(−0.28) 2.28(−0.12) 
Hexatriene 3Ag 4.15 5.20(1.05) 4.86(0.71) 4.47(0.32) 5.10(0.95) 4.64(0.49) 4.19(0.04) 3.57(−0.58) 1.56(−2.60) 3.94(−0.21) 4.04(−0.11) 
Hexatriene 1Ag 5.09 5.42(0.33) 4.99(−0.10) 4.49(−0.60) 5.34(0.25) 4.86(−0.23) 4.32(−0.77) 6.88(1.79) 6.71(1.62) 5.50(0.41) 5.02(−0.07) 
Hexatriene 1Bu 5.1 5.03(−0.07) 5.29(0.19) 5.13(0.03) 4.94(−0.16) 5.09(−0.01) 4.88(−0.22) 5.33(0.23) 5.05(−0.05) 4.60(−0.50) 4.44(−0.66) 
Octetraene 3Bu 2.2 2.15(−0.05) 1.49(−0.71) 1.21(−0.99) 2.07(−0.13) 1.31(−0.89) 0.99(−1.21) 1.79(−0.41) −2.74(−4.94) 1.71(−0.49) 1.87(−0.33) 
Octetraene 3Ag 3.55 4.75(1.20) 4.00(0.45) 3.57(0.02) 4.66(1.11) 3.80(0.25) 3.30(−0.25) 3.01(−0.54) −1.12(−4.67) 3.26(−0.29) 3.36(−0.19) 
Octetraene 1Ag 4.47 5.02(0.55) 4.08(−0.39) 3.53(−0.94) 4.95(0.48) 3.96(−0.51) 3.37(−1.10) 6.39(1.92) 6.32(1.85) 4.80(0.33) 4.17(−0.30) 
Octetraene 1Bu 4.66 4.58(−0.08) 4.50(−0.16) 4.29(−0.37) 4.50(−0.16) 4.31(−0.35) 4.05(−0.61) 4.73(0.07) 4.47(−0.19) 3.96(−0.70) 3.78(−0.88) 
Decapentaene 1Bu 4.27 4.69(0.42) 4.60(0.33) 4.40(0.13) 4.61(0.34) 4.40(0.13) 4.13(−0.14) 4.88(0.61) 4.63(0.36) 4.07(−0.20) 3.86(−0.41) 
Cyclopropene 3B2 4.34 4.22(−0.12) 4.08(−0.26) 3.91(−0.43) 4.14(−0.20) 3.89(−0.45) 3.66(−0.68) 3.46(−0.88) 0.67(−3.67) 3.70(−0.64) 3.79(−0.55) 
Cyclopropene 1B2 7.06 5.87(−1.19) 7.29(0.23) 7.31(0.25) 5.80(−1.26) 7.12(0.06) 8.15(1.09) 6.67(−0.39) 6.40(−0.66) 6.09(−0.97) 5.91(−1.15) 
Cyclopentadiene 3B2 3.25 3.12(−0.13) 2.66(−0.59) 2.53(−0.72) 3.03(−0.22) 2.48(−0.77) 2.31(−0.94) 2.49(−0.76) −2.06(−5.31) 2.74(−0.51) 2.91(−0.34) 
Cyclopentadiene |$\underline{{}^3\!A_1}$|A13̲ 5.09 5.45(0.36) 5.33(0.24) 5.07(−0.02) 5.36(0.27) 5.15(0.06) 4.86(−0.23) 4.29(−0.80) 2.94(−2.15) 4.75(−0.34) 4.87(−0.22) 
Cyclopentadiene 1B2 5.55 5.16(−0.39) 5.46(−0.09) 5.45(−0.10) 5.07(−0.48) 5.31(−0.24) 5.26(−0.29) 5.43(−0.12) 5.12(−0.43) 4.95(−0.60) 4.88(−0.67) 
Cyclopentadiene |$\underline{{}^1\!A_1}$|A11̲ 6.31 6.01(−0.30) 6.42(0.11) 6.16(−0.15) 5.94(−0.37) 6.33(0.02) 6.05(−0.26) 7.87(1.56) 7.78(1.47) 6.40(0.09) 5.78(−0.53) 
Norbornadiene 3A2 3.72 3.77(0.05) 3.68(−0.04) 3.58(−0.14) 3.57(−0.15) 3.36(−0.36) 3.20(−0.52) 2.87(−0.85) −1.48(−5.20) 3.10(−0.62) 3.16(−0.56) 
Norbornadiene |$\underline{{}^3\!B_2}$|B23̲ 4.16 4.49(0.33) 4.76(0.60) 4.51(0.35) 4.28(0.12) 4.41(0.25) 4.11(−0.05) 3.27(−0.89) −1.32(−5.48) 3.63(−0.53) 3.78(−0.38) 
Norbornadiene 1A2 5.34 3.95(−1.39) 5.36(0.02) 5.34(0.00) 3.75(−1.59) 5.05(−0.29) 4.99(−0.35) 5.53(0.19) 5.30(−0.04) 4.70(−0.64) 4.40(−0.94) 
Norbornadiene |$\underline{{}^1\!B_2}$|B21̲ 6.11 4.58(−1.53) 6.84(0.73) 6.99(0.88) 4.38(−1.73) 6.51(0.40) 6.64(0.53) 7.02(0.91) 6.81(0.70) 5.28(−0.83) 4.94(−1.17) 
Furan 3B2 4.17 3.82(−0.35) 3.49(−0.68) 3.37(−0.80) 3.74(−0.43) 3.33(−0.84) 3.17(−1.00) 3.28(−0.89) −1.14(−5.31) 3.70(−0.47) 3.89(−0.28) 
Furan |$\underline{{}^3\!A_1}$|A13̲ 5.48 5.79(0.31) 5.37(−0.11) 5.12(−0.36) 5.71(0.23) 5.22(−0.26) 4.93(−0.55) 4.87(−0.61) 4.07(−1.41) 5.18(−0.30) 5.25(−0.23) 
Furan 1B2 6.32 5.08(−1.24) 6.57(0.25) 6.58(0.26) 5.00(−1.32) 6.42(0.10) 6.40(0.08) 6.28(−0.04) 5.97(−0.35) 5.93(−0.39) 5.88(−0.44) 
Furan |$\underline{{}^1\!A_1}$|A11̲ 6.57 6.61(0.04) 6.88(0.31) 6.65(0.08) 6.54(−0.03) 6.80(0.23) 6.56(−0.01) 7.89(1.32) 7.77(1.20) 6.58(0.01) 6.28(−0.29) 
Pridazine 1B1 3.78 1.51(−2.27) 2.68(−1.10) 2.91(−0.87) 1.27(−2.51) 2.33(−1.45) 2.50(−1.28) 4.93(1.15) 4.72(0.94) 3.60(−0.18) 3.15(−0.63) 
Pridazine 1A2 4.32 2.91(−1.41) 3.41(−0.91) 3.48(−0.84) 2.60(−1.72) 2.90(−1.42) 2.91(−1.41) 6.13(1.81) 5.95(1.63) 4.20(−0.12) 3.53(−0.79) 
s-tetrazine 3B3u 1.89 3.28(1.39) 2.19(0.30) 1.86(−0.03) 2.94(1.05) 1.64(−0.25) 1.24(−0.65) 2.40(0.51) 1.89(0.00) 1.47(−0.42) 1.15(−0.74) 
s-tetrazine 3Au 3.52 5.38(1.86) 3.67(0.15) 3.16(−0.36) 4.91(1.39) 2.91(−0.61) 2.30(−1.22) 4.67(1.15) 4.38(0.86) 3.15(−0.37) 2.54(−0.98) 
s-tetrazine 1B3u 2.24 3.78(1.54) 2.73(0.49) 2.41(0.17) 3.53(1.29) 2.33(0.09) 1.96(−0.28) 3.54(1.30) 3.32(1.08) 2.27(0.03) 1.85(−0.39) 
s-tetrazine 1Au 3.48 5.59(2.11) 3.91(0.43) 3.40(−0.08) 5.13(1.65) 3.18(−0.30) 2.58(−0.90) 5.71(2.23) 5.55(2.07) 3.54(0.06) 2.87(−0.61) 
Formaldehyde 3A2 3.5 1.65(−1.85) 3.15(−0.35) 3.40(−0.10) 1.53(−1.97) 2.85(−0.65) 3.04(−0.46) 3.67(0.17) 3.34(−0.16) 3.10(−0.40) 2.97(−0.53) 
Formaldehyde 1A2 3.88 2.00(−1.88) 3.68(−0.20) 3.97(0.09) 1.88(−2.00) 3.43(−0.45) 3.68(−0.20) 4.49(0.61) 4.31(0.43) 3.83(−0.05) 3.71(−0.17) 
Acetone 3A2 4.05 3.10(−0.95) 4.13(0.08) 4.24(0.19) 2.92(−1.13) 3.78(−0.27) 3.83(−0.22) 4.42(0.37) 4.13(0.08) 3.68(−0.37) 3.52(−0.53) 
Acetone 1A2 4.4 3.38(−1.02) 4.56(0.16) 4.66(0.26) 3.20(−1.20) 4.25(−0.15) 4.31(−0.09) 5.15(0.75) 4.98(0.58) 4.31(−0.09) 4.14(−0.26) 
Benzoquinone 3B1g 2.51 4.74(2.23) 2.93(0.42) 2.44(−0.07) 4.07(1.56) 2.05(−0.46) 1.53(−0.98) 3.31(0.80) 3.01(0.50) 1.93(−0.58) 1.42(−1.09) 
Benzoquinone 1B1g 2.78 4.98(2.20) 3.14(0.36) 2.65(−0.13) 4.36(1.58) 2.36(−0.42) 1.85(−0.93) 4.00(1.22) 3.84(1.06) 2.43(−0.35) 1.87(−0.91) 
DMABN 1A1 4.56 6.14(1.58) 4.96(0.40) 4.56(0.00) 6.04(1.48) 4.70(0.14) 3.90(−0.66) 5.45(0.89) 5.16(0.60) 4.52(−0.04) 4.19(−0.37) 
DMABN 1B1 4.25 5.68(1.43) 4.68(0.43) 4.24(−0.01) 5.58(1.33) 4.42(0.17) 4.21(−0.04) 5.54(1.29) 5.35(1.10) 4.39(0.14) 3.97(−0.28) 
Total   MSE 0.00 0.04 −0.14 −0.16 −0.24 −0.44 0.31 −0.95 −0.32 −0.49 
Total   MAE 0.92 0.40 0.38 0.89 0.43 0.60 0.79 1.82 0.37 0.49 
Singlets   MSE −0.23 0.10 −0.04 −0.38 −0.16 −0.30 0.84 0.62 −0.24 −0.56 
Singlets   MAE 1.08 0.36 0.34 1.08 0.34 0.54 0.89 0.83 0.33 0.56 
Triplets   MSE 0.29 −0.05 −0.27 0.12 −0.34 −0.62 −0.36 −2.93 −0.42 −0.41 
Triplets   MAE 0.71 0.46 0.44 0.64 0.54 0.67 0.67 3.08 0.42 0.41 

Unlike the results we observed for small molecules, in large systems, there is relatively a larger difference between pp-RPA and pp-TDA. The excitation energies calculated from pp-TDA is always lower than its corresponding pp-RPA. Therefore, the mean signed error shifted from close to zero to −0.16 eV ∼ −0.44 eV depending on the reference that is chosen. Even with such a shift, the mean absolute error does not change much for HF reference or B3LYP reference. However, for pp-TDA-PBE, both the mean signed error and mean absolute error are much larger than pp-RPA-PBE. These suggest that pp-TDA gives excitation energies slightly lower than pp-RPA and its performance is similar or slightly worse than pp-RPA. This observation is very different from TDDFT, where TDA results are often better than RPA results.

We presented CIS, TD-HF, TD-B3LYP, TD-PBE results in the last four columns of Table I. It can be seen that CIS overestimates excitation energies. Among all the tested methods, TD-HF is the only one that has stability issues: instead of real ones, the excitation energies are imaginary (shown with negative numbers in the table). When we calculate the error for TD-HF, we use the negative numbers as a punishment for instability issue. Therefore, for TD-HF there is a large negative signed error and large mean absolute error. For TD-B3LYP and TD-PBE, they mostly underestimate excitation energies with negative mean signed errors. Compare CIS and pp-RPA-HF, TD-B3LYP and pp-RPA-B3LYP, TD-PBE and pp-RPA-PBE, we can see the mean absolute errors are very similar in each comparison group, which suggests that our pp-RPA methods are comparable to the convention CIS and TDDFT methods for large molecules.

We investigated the error for the singlet excitation group and the triplet excitation group. For pp-RPA-HF and pp-TDA-HF, the singlet excitations have negative errors while triplet excitations have positive errors, and triplets excitations have slightly smaller mean absolute error. For the two DFT references, the results are opposite: triplet excitations have more negative errors and larger mean absolute errors. Therefore, singlet excitations are better described by pp-RPA with DFT references. A B3LYP example is shown in Fig. 3.

FIG. 3.

Error distributions for pp-RPA-B3LYP, pp-TDA-B3LYP, and TD-B3LYP for (a) singlet and (b) triplet excitations. It can be observed that (1) pp-RPA has nearly even distribution around the zero line with no systematic error, while TD-B3LYP tends to underestimate excitations; (2) pp-TDA excitations are lower than pp-RPA with slightly larger error; (3) for the pp- methods, most singlet excitations' errors are within −0.5 eV to 0.5 eV, they are much better described than triplet excitations with error distributed between −1 eV and 1 eV.

FIG. 3.

Error distributions for pp-RPA-B3LYP, pp-TDA-B3LYP, and TD-B3LYP for (a) singlet and (b) triplet excitations. It can be observed that (1) pp-RPA has nearly even distribution around the zero line with no systematic error, while TD-B3LYP tends to underestimate excitations; (2) pp-TDA excitations are lower than pp-RPA with slightly larger error; (3) for the pp- methods, most singlet excitations' errors are within −0.5 eV to 0.5 eV, they are much better described than triplet excitations with error distributed between −1 eV and 1 eV.

Close modal

It had been our concern for the pp-RPA and pp-TDA methods because within our regular treatment, non-HOMO excitations cannot be captured. When the molecule gets larger, the non-HOMO excitations might play a more important role. However, according to the large systems we investigated, the lowest few excitations very commonly have some HOMO excitation characters. To our surprise, even if the HOMO excitation contribution is as small as half in a TDDFT calculation, our approach still can capture that state well with a reasonably good excitation energy. These excitations with significant non-HOMO excitation characters are marked with underlines in Table I. This phenomenon is very similar to double excitations for TDDFT: although TDDFT cannot capture pure double excitations, it can accurately predict the excitation energies for many excitations with both double and single excitation characters.36 Therefore, even though our methods with regular treatments cannot capture those excitations with almost pure non-HOMO excitation characters, they are reliable in predicting most of the low-lying excitations with some HOMO-excitations characters. Note there is another way to capture non-HOMO excitations with pp-RPA: we can use the non-ground state of (N − 2)-electron systems as the reference, somewhat like the use of HF* in our previous work.23 

The starting points for pp-RPA calculations are converged (N − 2)-electron systems. Most systems of interest are closed-shell systems, therefore, the corresponding (N − 2)-electron systems sometimes show diradical nature if HOMO and HOMO-1 orbitals are very close. This diradical nature makes it difficult to converge (N − 2)-electron systems.

One extreme case is molecules with high symmetry. If HOMO and HOMO-1 happen to be degenerate, then the (N − 2) SCF calculation will be difficult. Even if in some cases the convergence can be reached, the orbital degenerate symmetry is not conserved. This category includes molecules such as benzene and triazine. For them, we still cannot solve the problem gracefully with the current approach. However, fortunately, these high degenerated cases are rare, especially in large molecules.

There are also some molecules without high symmetry but with (N − 2) convergence difficulties. This often happens with the PBE functional. However, when some Hartree-Fock exchange is present, which expands the HOMO-LUMO gap, SCF convergence is easier to be reached. In Table II, we show these molecules that cannot be converged with PBE functional but can be converged with B3LYP functional. It can be seen that in this set, pp-RPA-B3LYP still has small mean signed error and TD-B3LYP has a negative error. As to the absolute error, pp-RPA-B3LYP is 0.03 eV larger than the set in Table I while TD-B3LYP is 0.03 eV smaller which means pp-RPA does not perform as well as TD-B3LYP in this set. This relatively large error for pp-RPA might be related to the (N − 2)-electron system: PBE cannot converge, therefore, we anticipate the converged B3LYP results might also be not good enough. However, the 0.43 eV mean absolute error for pp-RPA-B3LYP is acceptable although it is larger than 0.34 eV for TD-B3LYP.

Table II.

Vertical excitation energies (in eV) from pp-RPA-B3LYP and TD-B3LYP. Numbers in the parenthesis indicate the error. Only B3LYP results are reported for the following molecules because PBE cannot converge their corresponding (N − 2)-electron systems.

MoleculeExciRefRPA-B3LYPTD-B3LYP
Pyrrole 3B2 4.48 3.79(−0.69) 4.07(−0.41) 
Pyrrole 3A1 5.51 5.28(−0.23) 5.22(−0.29) 
Pyrrole 1A1 6.37 6.69(0.32) 6.33(−0.04) 
Pyrrole 1B2 6.57 6.80(0.23) 5.98(−0.59) 
Imidazole 3A 4.69 3.94(−0.75) 4.24(−0.45) 
Imidazole 3A 5.79 5.14(−0.65) 5.35(−0.44) 
Imidazole 1A 6.19 6.54(0.35) 6.16(−0.04) 
Imidazole 1A 6.93 7.27(0.34) 6.35(−0.58) 
Pyridine 3A1 4.06 4.01(−0.05) 3.91(−0.15) 
Pyridine 3B2 4.64 3.95(−0.69) 4.46(−0.18) 
Pyridine 1B2 4.85 6.00(1.15) 5.45(0.60) 
Pyridine 1A1 6.26 6.97(0.71) 6.19(−0.07) 
Pyrazine 1B3u 3.95 4.41(0.46) 3.93(−0.02) 
Pyrazine 1Au 4.81 5.15(0.34) 4.69(−0.12) 
Pyrimidine 1B1 4.55 4.32(−0.23) 4.25(−0.30) 
Pyrimidine 1A2 4.91 4.84(−0.07) 4.60(−0.31) 
Formamide 3A 5.74 4.39(−1.35) 5.10(−0.64) 
Formamide 1A 7.39 8.21(0.82) 7.53(0.14) 
Acetamide 3A 5.88 5.49(−0.39) 5.25(−0.63) 
Acetamide 1A 7.27 6.68(−0.59) 7.10(−0.17) 
Propanamide 3A 5.9 6.47(0.57) 5.27(−0.63) 
Propanamide 1A 7.2 7.33(0.13) 7.00(−0.20) 
Cytosine 1A 4.66 4.74(0.08) 4.59(−0.07) 
Thymine 1A 5.2 5.57(0.37) 4.89(−0.31) 
Uracil 1A 5.35 5.90(0.55) 5.08(−0.27) 
Adenine 1A 5.25 5.16(−0.09) 4.91(−0.34) 
N-phenylpyrrole 1B2 4.85 4.83(−0.02) 4.51(−0.34) 
N-phenylpyrrole 1A1 5.13 5.23(0.10) 4.57(−0.56) 
N-phenylpyrrole 1A1 5.94 6.19(0.25) 4.87(−1.07) 
MSE     0.03 −0.29 
MAE     0.43 0.34 
MoleculeExciRefRPA-B3LYPTD-B3LYP
Pyrrole 3B2 4.48 3.79(−0.69) 4.07(−0.41) 
Pyrrole 3A1 5.51 5.28(−0.23) 5.22(−0.29) 
Pyrrole 1A1 6.37 6.69(0.32) 6.33(−0.04) 
Pyrrole 1B2 6.57 6.80(0.23) 5.98(−0.59) 
Imidazole 3A 4.69 3.94(−0.75) 4.24(−0.45) 
Imidazole 3A 5.79 5.14(−0.65) 5.35(−0.44) 
Imidazole 1A 6.19 6.54(0.35) 6.16(−0.04) 
Imidazole 1A 6.93 7.27(0.34) 6.35(−0.58) 
Pyridine 3A1 4.06 4.01(−0.05) 3.91(−0.15) 
Pyridine 3B2 4.64 3.95(−0.69) 4.46(−0.18) 
Pyridine 1B2 4.85 6.00(1.15) 5.45(0.60) 
Pyridine 1A1 6.26 6.97(0.71) 6.19(−0.07) 
Pyrazine 1B3u 3.95 4.41(0.46) 3.93(−0.02) 
Pyrazine 1Au 4.81 5.15(0.34) 4.69(−0.12) 
Pyrimidine 1B1 4.55 4.32(−0.23) 4.25(−0.30) 
Pyrimidine 1A2 4.91 4.84(−0.07) 4.60(−0.31) 
Formamide 3A 5.74 4.39(−1.35) 5.10(−0.64) 
Formamide 1A 7.39 8.21(0.82) 7.53(0.14) 
Acetamide 3A 5.88 5.49(−0.39) 5.25(−0.63) 
Acetamide 1A 7.27 6.68(−0.59) 7.10(−0.17) 
Propanamide 3A 5.9 6.47(0.57) 5.27(−0.63) 
Propanamide 1A 7.2 7.33(0.13) 7.00(−0.20) 
Cytosine 1A 4.66 4.74(0.08) 4.59(−0.07) 
Thymine 1A 5.2 5.57(0.37) 4.89(−0.31) 
Uracil 1A 5.35 5.90(0.55) 5.08(−0.27) 
Adenine 1A 5.25 5.16(−0.09) 4.91(−0.34) 
N-phenylpyrrole 1B2 4.85 4.83(−0.02) 4.51(−0.34) 
N-phenylpyrrole 1A1 5.13 5.23(0.10) 4.57(−0.56) 
N-phenylpyrrole 1A1 5.94 6.19(0.25) 4.87(−1.07) 
MSE     0.03 −0.29 
MAE     0.43 0.34 

We derived and implemented a Davidson algorithm for pp-RPA. The formal scaling is reduced to O(N4), which makes it possible to be applied to larger systems. We used this algorithm and carried out benchmark tests on pp-RPA and pp-TDA on HF, B3LYP, and PBE references. Despite some convergence problems with (N − 2)-electron systems, we are able to accurately capture lowest few excitations for systems with converged SCF calculations. Among them, DFT references have better performance with both clean ordered spectrum and accurate excitation energies. Excitations calculated from pp-TDA are lower than pp-RPA. Compared to TDDFT, there is no systematic bias for pp-RPA with the mean signed error close to zero. The mean absolute error with B3LYP or PBE references is similar to that of TDDFT, which suggests the pp-RPA is a comparable method to TDDFT for the test molecules. Moreover, despite some concerns for non-HOMO excitations, in many cases, excitations with relatively large non-HOMO excitation contributions are also well described in terms of the excitation energy, as long as there is also a considerably amount of HOMO excitation contribution. Therefore, the pp-RPA is also a reliable method to solve the lowest few excitations problems in large systems and it can potentially be a good alternative to TDDFT methods when TDDFT faces significant challenges, as in the case of charge transfer, double and Rydberg excitations.23 

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

1.
A. D.
Mclachlan
and
M. A.
Ball
,
Rev. Mod. Phys.
36
,
844
(
1964
).
3.
T.
Ziegler
,
A.
Rauk
, and
E.
Baerends
,
Theor. Chim. Acta
43
,
261
(
1977
).
4.
M. E.
Casida
, “
Time-dependent density functional response theory of molecular systems: Theory, computational methods, and functionals
,” in
Recent Developments and Applications of Modern Density Functional Theory
,
Theoretical and Computational Chemistry
Vol.
4
, edited by
J.
Seminario
(
Elsevier
,
1996
), pp.
391
439
.
5.
A. I.
Krylov
,
Acc. Chem. Res.
39
,
83
(
2006
).
6.
A.
Dreuw
and
M.
Head-Gordon
,
Chem. Rev.
105
,
4009
(
2005
).
7.
E.
Runge
and
E. K. U.
Gross
,
Phys. Rev. Lett.
52
,
997
(
1984
).
8.
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
Phys. Rev. A
88
,
030501
(R) (
2013
).
9.
D.
Peng
,
S. N.
Steinmann
,
H.
van Aggelen
, and
W.
Yang
,
J. Chem. Phys.
139
,
104112
(
2013
).
10.
G. E.
Scuseria
,
T. M.
Henderson
, and
I. W.
Bulik
,
J. Chem. Phys.
139
,
104113
(
2013
).
11.
D.
Bohm
and
D.
Pines
,
Phys. Rev.
82
,
625
(
1951
).
12.
D.
Pines
and
D.
Bohm
,
Phys. Rev.
85
,
338
(
1952
).
13.
A.
Heßelmann
and
A.
Gorling
,
Mol. Phys.
109
,
2473
(
2011
).
14.
H.
Eshuis
,
J.
Bates
, and
F.
Furche
,
Theor. Chem. Acc.
131
,
1084
(
2012
).
15.
X.
Ren
,
P.
Rinke
,
C.
Joas
, and
M.
Scheffler
,
J. Mater. Sci.
47
,
7447
(
2012
).
16.
P.
Ring
and
P.
Schuck
,
The Nuclear Many-Body Problem
(
Springer
,
2004
).
17.
J.-P.
Blaizot
,
Quantum Theory of Finite Systems
(
MIT Press
,
Cambridge, MA
,
1986
).
18.
Y.
Yang
,
H.
van Aggelen
,
S. N.
Steinmann
,
D.
Peng
, and
W.
Yang
,
J. Chem. Phys.
139
,
174110
(
2013
).
19.
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
J. Chem. Phys.
140
,
18A511
(
2014
).
20.
J. J.
Shepherd
,
T. M.
Henderson
, and
G. E.
Scuseria
,
J. Chem. Phys.
140
,
124102
(
2014
).
21.
J. J.
Shepherd
,
T. M.
Henderson
, and
G. E.
Scuseria
,
Phys. Rev. Lett.
112
,
133002
(
2014
).
22.
N.
Shenvi
,
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
J. Chem. Phys.
141
,
024119
(
2014
).
23.
Y.
Yang
,
H.
van Aggelen
, and
W.
Yang
,
J. Chem. Phys.
139
,
224105
(
2013
).
24.
D.
Peng
,
H.
van Aggelen
,
Y.
Yang
, and
W.
Yang
,
J. Chem. Phys.
140
,
18A522
(
2014
).
25.
D.
Zhang
,
D.
Peng
,
P.
Zhang
, and
W.
Yang
, “
Analytic gradient, geometric optimization and excited state potential energy surface from the particle-particle random phase approximation
Phys. Chem. Chem. Phys.
(unpublished).
26.
D.
Peng
,
Y.
Yang
,
P.
Zhang
, and
W.
Yang
, “
Restricted second random phase approximations and Tamm-Dancoff approximations for electronic excitation energy calculations
,”
J. Chem. Phys.
(submitted).
27.
M.
Schreiber
,
M. R.
Silva-Junior
,
S. P. A.
Sauer
, and
W.
Thiel
,
J. Chem. Phys.
128
,
134110
(
2008
).
28.
M. R.
Silva-Junior
,
M.
Schreiber
,
S. P. A.
Sauer
, and
W.
Thiel
,
J. Chem. Phys.
129
,
104103
(
2008
).
29.
M. J. G.
Peach
,
P.
Benfield
,
T.
Helgaker
, and
D. J.
Tozer
,
J. Chem. Phys.
128
,
044118
(
2008
).
30.
E. R.
Davidson
,
J. Comput. Phys.
17
,
87
(
1975
).
31.
G. G.
Sleijpen
and
H.
Van der Vorst
,
SIAM J. Matrix Anal. Appl.
17
,
401
(
1996
).
32.
G.
Sleijpen
and
H.
Van der Vorst
,
SIAM Rev.
42
,
267
(
2000
).
33.
H.
Weiss
,
R.
Ahlrichs
, and
M.
Hser
,
J. Chem. Phys.
99
,
1262
(
1993
).
34.
An in-house program for qm/mm simulations (http://www.qm4d.info).
35.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
 et al, Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford, CT,
2009
.
36.
S.
Hirata
and
M.
Head-Gordon
,
Chem. Phys. Lett.
302
,
375
(
1999
).