Excitation energies from particle-particle random phase approximation: Davidson algorithm and benchmark studies

To predict electronic excited states with high accuracy and efficiency has been a great desire for theoretical chemists for a long time. Many theoretical approaches^1–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 excitations⁶ 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 physics^16,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.⁸ 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(N⁶) to O(N⁴) within the lCCD formulation by using the tensor-hyper-contraction technique.²² 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.²³ 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 approach²⁴ 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.²⁴ 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.²⁵ 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,²⁶ which are prohibited by their formal first-order formulations.

Despite these efforts, we were still solving the pp-RPA equation using the direct diagonalization approach with a heavy computational cost O(N⁶), 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(N⁴) for a single eigenroot in Sec. II. Then in Sec. III we carry out benchmark tests with molecules in Thiel's^27,28 and Tozer's²⁹ 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,

in which

where

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

and the linear response equation with the integrated form is

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: δD^s = δD_int + δ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 K⁰ is simply

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

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

with

and

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

in which

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

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

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(N⁶) and the memory space cost is O(N⁴), with N being the maximum of N_vir and N_occ. 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 method³⁰ to solve the lowest few pair-addition eigenroots for the pp-RPA equation with an O(N⁴) time cost.

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

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,u_k), we approximate it by (⁠

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

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

Multiply V^T to the left, we obtain

with

The matrix

If the norm of r_k 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 u_k and

Then the original eigenvalue problem can be written as

Assume the approximated eigenvalue

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

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

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 (m ⩾ n). 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(N⁴) and memory cost is O(N²).

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(N²).

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

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

We implemented the above Davidson iterative method on the spin-separated and spin-adapted pp-RPA and pp-TDA equations¹⁸ in QM4D package.³⁴ Then we use it to benchmark excitation energies calculated with pp-RPA and compare with TDDFT results calculated with Gaussian 09.³⁵ We use molecules with the number of atoms larger than five in Thiel's^27,28 and Tozer's²⁹ 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.²³ 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.

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.

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.

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.³⁶ 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.²³ 

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.

We derived and implemented a Davidson algorithm for pp-RPA. The formal scaling is reduced to O(N⁴), 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.²³ 

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.