Excitation energies from particle-particle random phase approximation with accurate optimized effective potentials

The accurate description of electronic excited states has long been an important and challenging topic in theoretical chemistry. Many methods that can deal with excited states have been developed in the past few decades. Among these methods, time-dependent density functional theory (TDDFT)^1–3 has been widely used, thanks to its good balance of theoretical accuracy and computational efficiency. However, TDDFT has inherent problems in some challenging cases, such as double excitations, Rydberg excitations, charge transfer excitations, and diradical systems.⁴ Recently our group has developed the particle-particle random phase approximation (pp-RPA)^5–11 into an effective approach to address these challenging excitation problems and we found that the pp-RPA can successfully describe many of them. Nevertheless, as a post self-consistent-field (SCF) method, the final pp-RPA results depend on the SCF density functional reference, or more specifically, the orbital energies and orbital shapes. It has been demonstrated that better results of TDDFT calculations can be achieved^12,13 with the exact KS density functional theory,¹⁴ utilizing the optimized effective potential (OEP).^15–18 In this paper, we implement the pp-RPA with the OEP, based on the Wu-Yang approach.^17,19 With an accurate input density, the Wu-Yang direct optimization method can produce highly accurate local Kohn-Sham (KS) exchange-correlation potential and therefore the corresponding KS orbitals and orbital energies. This makes it possible to investigate the theoretical limit of the pp-RPA with the exact KS theory.

We first briefly review the pp-RPA and OEP methods in Sec. II. Then in Sec. III, we apply the pp-RPA with the OEP reference to double excitations, Rydberg excitations, diradicals, and regular single excitations, in order to investigate the theoretical limit of the pp-RPA based on local KS references. We also compare the results with calculations based on approximate functionals from our previous work.^5,9,20 We give our final conclusion in Sec. IV.

The pp-RPA formula can be derived in a variety of ways, including the equation-of-motion,^21,22 adiabatic connection pairing matrix fluctuation,^5,6 and time-dependent density functional theory with a pairing field (TDDFT-P).⁷ Here, we review the TDDFT-P method as this method justifies the use of the DFT references.

The general Hamiltonian with a pairing potential takes the form

where $T^$⁠, $V^$⁠, and $Ŵ$ represent the kinetic energy operator, external regular potential operator, and the two electron interaction operator, respectively. $D^$ is the external pairing field operator,

where h.c. stands for the Hermitian conjugate. Note that the total number of electrons in the system can change under a finite $D*(x,x′)$⁠. The pairing matrix under the pairing field is defined as

Now we perturb the interacting system with a small pairing field δ$D(x,x′;t)$⁠. Then the change of the pairing matrix that responds to the perturbation is

where the particle-particle (pp) response function K is

Note that H denotes the Heisenberg picture. For a non-interacting system with the Fourier transform from the time domain to the energy domain and the variables from coordinate representation to orbital representation, the pp linear response function is

where F denotes the Fermi level. In the following discussion, we will use p,q,r,s for general orbitals, i,j,k,l for occupied orbitals, and a,b,c,d for virtual orbitals. Thus, the linear response equation of the pairing matrix to the pairing field can be expressed as

According to the linear response theory, when δκ changes, there arises an internal mean-field paring potential δD_int, which further influences δκ back. Here we make another approximation that the response of δD^s to δκ is adiabatic. The adiabatic pp response kernel L is thus,

where $⟨pq||rs⟩$ is the antisymmetrized two-electron integral. Combining two parts of perturbations together, we reach the TDDFT-P equation⁷ 

where

Note that to avoid redundancy, we restrict the index $a>b$ and $i>j$⁠. Furthermore, by taking the external perturbation δD to 0, the eigenvalue equation is thus

Since the kernel L remains unknown, we use the antisymmetrized bare Coulomb interaction as the kernel, which leads to the pp-RPA equation. And in the pp-RPA, the simplified matrix elements are known explicitly now,

with the eigenvalues being two-electron addition energies

and two-electron removal energies

If we compute from an two-electron deficient (N − 2) system, the excitation energy is then

In order to get the excitation energies, we only need the orbitals and orbital energies of the non-interacting (N − 2)-electron system. The orbitals of (N − 2)-electron systems are more contracted, but qualitatively similar to those of N-electron systems, making it possible to qualitatively identify the character of each double electron addition excitation and consequently the corresponding N-electron system configuration. In this work, we have used four ways to obtain this information: (1) approximate density functional like PBE,²³ (2) approximate functionals of density matrix, like HF and B3LYP^24,25 in the generalized KS calculations, (3) OEP from accurate electron densities, and (4) OEP from densities calculated with approximate functionals of the density matrix.

We use the accurate and efficient direct optimization method for the determination of the OEP and its associated orbitals and orbital energies from a given electron density, developed by Wu and Yang.^17,19 Consider an N-electron closed-shell system with a given input electron density ρ_in, the Levy constrained search²⁶ shows that the kinetic energy is

Therefore, the evaluation of the kinetic energy is turned into a constrained minimization problem, where the Lagrangian is

And the stationary condition will be

The determinant, $Ψdet$⁠, is thus an implicit functional of $v$(r),

Then we can rewrite the Lagrangian W_s as a functional only depending on $v$(r),

As a consequence, the OEP problem is simplified into the maximization of W_s with respect to the potential $v$(r). The most efficient way to carry out the OEP is by expanding the potential $v$(r),

where $v$_ext(r) is the external potential, $v$₀ is a fixed reference potential, and the last term is a linear combination with a set of finite basis functions {g_t(r)} and the coefficients {b_t}. We choose the Fermi-Amaldi potential as the reference potential,²⁷ namely,

which can provide the correct long-range behavior and a set of Gaussian functions for {g_t(r)} for calculation convenience. Note that ρ₀(r) is a fixed reference electron density. Now the OEP problem turns into the optimization problem of $Ws[Ψdet[v(r)],v(r)]$ with respect to the expansion coefficients {b_t}. One way to carry out this optimization is by applying the Newton method. One can obtain the first and second order derivatives (gradient and Hessian),

and then iteratively solve the equation

where the gradient vector is defined as $gt=∂Ws∂bt$⁠. We also apply the BFGS method²⁸ in some cases as an alternative optimization method to avoid the inversion of the Hessian.

We use the following procedure for excitation energy calculations:

1. Perform a self-consistent field (SCF) or post-SCF calculation for the (N − 2)-electron system to obtain the total electron density;

2. use this density as an input for the OEP calculation to obtain orbitals and eigenvalues;

3. calculate the excitation energies with the pp-RPA.

All electron densities are calculated with Gaussian 09.²⁹ Our OEP and spin-adapted pp-RPA²⁰ methods are implemented in the QM⁴D package.³⁰ All calculations are compared with Refs. 5, 8, and 9 and for the regular single excitation calculations, we use molecules from Ref. 31. In the (N − 2)-electron system calculations, aug-cc-pVDZ is used for the molecules in Secs. III B and III C and aug-cc-pVQZ is used for the molecules in Secs. III D and III E. In the OEP calculations, the potential basis set is expanded in s, p, and d-type even tempered functions. In the benchmark test of regular single excitations, we adopt the basis set in Wu’s work,³² where for non-hydrogen atoms, the basis has 19s/12p/3d functions. And for small molecules, the basis has 38s/12p/7d functions considering the computational cost.

In general, the basis set convergence needs to be verified with both atomic orbital basis and OEP potential basis. The basis set convergence test with correlation consistent basis set (cc-pVXD, with X = D, T, Q) and augmented correlation consistent basis set (aug-cc-pVXD, with X = D, T, Q) has been investigated previously⁸ and here the aug-cc-pVDZ basis is chosen for the convergence test of different potential basis sets. We used butadiene for the potential basis convergence test. The result is shown in Fig. 1. The convergence with all three input densities (HF, B3LYP, and CCSD) is fast. According to Table I, 1 and 2 only contain s-type Gaussian functions in the potential basis set, and 3, 4, and 5 contain s, p, and d−type Gaussian functions. We observed that excitation energies slightly changed with the existence of p and d−type functions. This is because p and d−type functions can capture the non-spherical effects. The potential basis set with 19s/12p/3d functions (x-label 4 in Fig. 1) is considered to be converged. Therefore, we took this potential basis set and aug-cc-pVDZ as the atomic basis set for regular single excitations. And for small molecules, 38s/12p/7d is used since the computational cost of small molecules is low.

We tested four small molecules, shown in Table II. Both single and double excitations are captured by the pp-RPA. The pp-RPA with OEP orbitals and orbital energies with an input coupled-cluster singles and doubles (CCSD) density (pp-RPA-OEP-CCSD) provides the best excitation energies, where the mean signed error is 0.00 eV and the mean absolute error is 0.10 eV. However, other density sources do not perform as well. For the pp-RPA-OEP-HF method, the mean absolute error is 0.85 eV, while that of the pp-RPA-HF is 0.60 eV, similarly for the pp-RPA-OEP-B3LYP and the pp-RPA-B3LYP (0.38 eV and 0.21 eV, respectively). This is reasonable because the CCSD density is an accurate density source so that the OEP with a CCSD density can provide accurate orbitals and orbital energies, while the quality of the electron densities calculated by the HF and B3LYP is not good enough for these small molecules. It is encouraging that the best KS reference here from O-CCSD (OEP with CCSD density) provides excellent accuracy, the best overall in this test set. However, the KS OEP results from other densities (O-HF and O-B3LYP) do not perform better than the corresponding generalized KS (GKS) results (HF or B3LYP).

In terms of the accuracy, The equation-of-motion coupled cluster method restricted to single and double excitations (EOM-CCSD)³³ performs slightly better than the OEP-pp-RPA. However, the main goal of this work is to investigate the limit of the pp-RPA, rather than promoting the computational method of the pp-RPA with the electron density from CCSD calculations.

Rydberg excitations of the above molecules are also computed, shown in Table III. The pp-RPA-OEP-CCSD overestimates the excitation energies, which behaves slightly worse than the pp-RPA with PBE or B3LYP. In general, the pp-RPA does not predict Rydberg excitations as accurately as it does for the low-lying excitations with both approximated density functionals or functionals of density matrices.

When two electrons are removed from a diradical system, the remaining (N − 2)-electron system is generally a closed-shell singlet system, which can be well described by some well-known density functionals like B3LYP or PBE. This is ideal for the pp-RPA calculation because it starts with an (N − 2)-electron system. Therefore, the pp-RPA is expected to describe properties of diradicals well, for example, singlet-triplet (ST) gaps. Two types of ST gaps are usually considered: vertical ST gaps, $Egv$⁠, are directly calculated from the pp-RPA equation, while adiabatic gaps, $Ega$⁠, need a correction term given by the geometry difference of the N-electron system. In general, two types of corrections are carried out in calculations⁹ 

In our cases, Eq. (35) is adopted as the correction equation.

We tested the singlet-triplet gaps for diradicals with five molecules, four of which are diatomic diradicals and the other one is a carbene-like diradical. The ST gaps for diatomics are in Table IV. In this case, the pp-RPA-OEP results are consistently better than the pp-RPA results calculated with approximate density functionals. Especially for the pp-RPA-OEP-HF, the mean absolute error decreases from 4.1 to 2.7 kcal/mol. This result shows that the pp-RPA calculated with the DFT references is better than the HF references. Here we observe the OEP HOMO-LUMO degeneracy problem with the CCSD densities of NH and OH⁺. This phenomenon was addressed by Rohr³⁵ and it also affects the ST gaps with the PBE densities. The pp-RPA-OEP-PBE differs from pp-RPA-PBE mainly in the long-range behavior of the local potential; however, for $NH$ and $OH+$⁠, their results differ by 1.9 and 4.2 kcal/mol, respectively.

The ST gaps of methylene are shown in Table V. Methylene has two non-degenerate frontier orbitals—one σ orbital and one π orbital with different symmetries. The ground state is $B13$⁠, where one electron occupies the σ orbital and one occupies the π orbital. Three excited states, $A11$⁠, $B11$⁠, and 2¹A₁, represent three singlet states, σ², σ¹π¹, and π², respectively. Note that $B11$ is an open-shell singlet. The result shows that for $A11$⁠, the pp-RPA-OEP-HF gives the best estimation, while for other two states, the pp-RPA-OEP results with other three densities are better than those with the HF density. In fact, for methylene, the first adiabatic gap had historically been a puzzle. The reason is that the singlet energy with $A11$ geometry is fairly close to or lower than the triplet energy.^38–40 In our calculation, the pp-RPA-OEP-CCSD underestimates this gap (0.0 kcal/mol) and pp-RPA-HF and pp-RPA-B3LYP also underestimate this gap.

The regular single excitation results of 12 molecules are shown in Table VI, which contains 19 triplet and singlet excitations. The total mean absolute error (MAE) shows that the pp-RPA-OEP methods perform consistently better than the pp-RPA calculated with approximate density functionals. The pp-RPA-OEP-CCSD method still gives the best excitation energies. For HF, the MAE of the pp-RPA-OEP-HF is 0.52 eV, which significantly improves the pp-RPA-HF result (0.92 eV). And the MAE of the pp-RPA-OEP-B3LYP is slightly better the pp-RPA-B3LYP (0.38 eV versus 0.40 eV). This reflects that excitation energies from the pp-RPA are better described with KS density functional references rather than GKS references from density matrix functional in our regular single excitations. In addition, from the total mean signed error (MSE), we notice that the pp-RPA-OEP underestimates excitation energies, which behaves similar with the pp-RPA-PBE, shown in Fig. 2. However, compared with the TDDFT results from B3LYP and PBE functionals, pp-RPA has better MSEs. Standard TDDFT calculations with B3LYP and PBE strongly underestimate both singlet and triplet excitations, as shown in Table VI.

Moreover, the singlet and triplet excitations behave differently with the pp-RPA-OEP methods. Our method improves the accuracy of predictions of singlet excitations. For HF, the improvement is 0.6 eV in terms of the MAE. However, the error in triplet excitations is nearly unchanged. And for B3LYP, triplet excitations with the pp-RPA-OEP-B3LYP have an MAE of 0.48 eV, which is slightly larger than the pp-RPA-B3LYP (0.46 eV). Therefore, singlet excitations are better described than triplet excitations by the pp-RPA for regular single excitations. MSEs are visually shown in Fig. 2, where the pp-RPA-OEP underestimates both singlet and triplet excitations. This shows that the pp-RPA with DFT references will underestimate excitation energies, which is similar to previous results.²⁰ 

We combined the OEP and the pp-RPA method together to calculate excitation energies. We have applied our method to double excitations, Rydberg excitations, diradicals, and regular single excitations. Among our SCF and post-SCF methods, the CCSD method provides the most accurate electron densities. And consequently, the pp-RPA-OEP-CCSD approach overall gives the best estimation of excitations. Other input densities are not as good in all systems. For example, the MAE of the pp-RPA-OEP-B3LYP (0.38 eV) is larger than that of the pp-RPA directly calculated with B3LYP (0.21 eV) in small molecules calculations; however, for regular single excitations, the pp-RPA-OEP performs slightly better. Moreover, the MSE shows that the pp-RPA-OEP method generally underestimates excitation energies. The pp-RPA-OEP-CCSD approaches the theoretical limit of the pp-RPA in excitation energy calculations based on KS references. However, the MAE in the regular single excitations of the pp-RPA-OEP-CCSD is still 0.35 eV. This MAE illustrates the intrinsic limitation of the pp-RPA method based on KS references. In order to achieve excitation energies with a higher accuracy, one needs to go beyond the pp-RPA approximation with the consideration of more general kernel within the TDDFT-P formulation.

Y.J. appreciates funding from GPNano fellowship from nanoscience at Duke University. We acknowledge support from the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0012575(Y.J. and D.Z.), and the National Science Foundation (Grant No. CHE-1362927) (Y.Y. and W.Y.). Y.J. also appreciates Mr. Xinchang Xie for helpful discussions.