Charge transfer excitations from particle-particle random phase approximation—Opportunities and challenges arising from two-electron deficient systems

Charge transfer (CT) excitations are frequently encountered in materials and biological sciences.^1–4 They play an important role in electron transfer and energy conversion. Therefore, it is highly desirable to understand and describe CT excitations well. However, the accurate theoretical description of CT excitations remains a great challenge. Although the time-dependent density functional theory (TDDFT)^5,6 has been widely applied in many electronically excited systems thanks to its good balance of theoretical accuracy and computational efficiency, under commonly used approximate density functional kernels, it frequently fails for CT excitations.^7,8 In particular, TDDFT tends to greatly underestimate CT excitation energies as a result of the self-interaction error.^7,8 In order to reduce this error, one has to introduce the Hartree-Fock (HF) exchange in the kernel, and therefore, a common remedy is the range-separated (RS) hybrid kernels. The correct asymptotic trend for CT excitations is −1/R (R is the distance between the donor and the acceptor), and, in principle, the full HF exchange in the long range is required to reproduce this trend. Apart from the qualitative asymptotic behavior, the quantitative prediction of CT excitation energies is also challenging. Within the framework of TDDFT, Stein et al.⁹ have developed a system-dependent RS functional with a tuned parameter to overcome this challenge. The tuning procedure is performed to minimize the difference between the ionization potential (IP) and the orbital energy of the highest occupied molecular orbital (HOMO), based on the exact theoretical condition that the IP should equal the HOMO energy in the Kohn-Sham (KS) theory. Although this tuning is system-dependent and, in principle, a small change in the distance between the donor and the acceptor will lead to a new system and therefore a new tuning parameter, this tuned functional is shown to perform extremely well in quantitatively predicting CT excitations. The M06-HF functional by Zhao and Truhlar¹⁰ is also capable of accurately predicting CT energies using a global full HF exchange.

Other methods that are outside the framework of TDDFT have also been developed to describe CT excitations. Constrained variational density functional theory (CV-DFT), especially a relaxed theory to the fourth order (R-CV(4)-DFT), can achieve similar results to those using the scheme by Stein et al. even at the local density approximation (LDA) level.¹¹ Delta self-consistent field (ΔSCF)^11,12 and perturbative ΔSCF¹³ can also yield good results, even though these methods have a longer history than the linear-response TDDFT. The constrained DFT by Wu and Van Voorhis^14,15 and a scheme by Theophilou et al.^16,17 based on a minimization on the HF functional within a constrained space are also able to give good CT results.

The particle-particle random phase approximation (pp-RPA) has been a textbook method for treating correlation in nuclear physics for a long time.^18,19 In recent years, it was introduced by van Aggelen et al. to atomic and molecular systems for treating electron correlation and was successfully combined with DFT as well as the conventional HF theory.^20,21 Recently, the pp-RPA was developed into a tool for studying atomic and molecular excitations using a two-electron deficient (N-2) system as the initial reference.²² It is able to describe double excitations,^22,23 Rydberg excitations,^22,24 CT excitations,^22,23 diradical systems,^23,25 and conical intersections,²⁶ all of which are challenging for TDDFT. The pp-RPA has a similar computational cost to TDDFT.²⁷ Therefore, it could be a promising method that complements TDDFT.

In previous preliminary investigations on CT excitations,²² we found that the pp-RPA is capable of correctly describing the asymptotic −1/R trend with a specially designed initial reference named HF*. Nonetheless, some important questions still remain. First, the reason of the correct asymptotic trend has not been fully discussed. Second, the qualitative performance has not been tested when the donor and acceptor are spatially close. Third, our previous tests only used the specially designed HF* reference,²² and a more comprehensive functional performance test is in need. In this paper, we will focus on these three remaining issues and discuss the potential opportunities and challenges.

The pp-RPA can be derived from multiple independent approaches, including the adiabatic connection and pairing matrix fluctuations,^20,21 the equations of motion,^18,28,29 and time-dependent density functional theory with a pairing field (TDDFT-P).³⁰ The final pp-RPA equation is a generalized eigenvalue problem,

with

where a, b, c, d are virtual orbital indices and i, j, k, l are occupied orbital indices with the restrictions that $a>b$⁠, $c>d$⁠, $i>j$⁠, and $k>l$⁠. The brackets are defined as

and $⟨pq||rs⟩≡⟨pq|rs⟩−⟨pq|sr⟩$⁠.

The pp-RPA equation (Eq. (1)) describes the two-electron addition and two-electron removal processes. If we switch to an N-2 reference and study the two-electron addition processes, we will obtain a series of neutral N-electron states, including the ground state as well as some electronically excited states.²² The excitation energies can be predicted by taking the differences of two two-electron addition energies,

where N and N − 2 denote the number of electrons, and 0 and n denote the ground and excited states, respectively.

In practice, in order to improve the computational efficiency, we often use the spin-separated or spin-adapted forms of the pp-RPA equation.²⁹ The spin separation can always be performed whereas the singlet-triplet spin adaptation is only possible for restricted closed-shell N-2 reference. Because we will use both spin-adapted and spin-separated (⁠$αβ,αβ$⁠) equations in later analysis process, here we briefly review their working forms.

The spin-separated (⁠$αβ,αβ$⁠) equation describes the process of adding or removing one α-β electron pair. Its form is

where $spin=(αβ,αβ)$⁠, and the matrix elements are

When the initial SCF N-2 calculation is restricted closed-shell, the spin separated (⁠$αβ,αβ$⁠) equation can be further resolved into one equation for singlet pairs and another equation for triplet pairs. This is the singlet-triplet spin adaptation. The spin-adapted singlet and triplet equations describe adding or removing a singlet or a triplet pair of electrons, respectively. Their forms are

where the multiplicity mult is either singlet (s) or triplet (t). The elements in the triplet matrix are

with the restriction that $a>b$⁠, $c>d$⁠, $i>j$⁠, and $k>l$⁠, while the elements in the singlet matrix are

with the restriction that $a≥b$⁠, $c≥d$⁠, $i≥j$⁠, and $k≥l$⁠.

When a system of considerable dimensions is the target of pp-RPA calculations, the above formalism may have serious limitations concerning computational resources. One important bottleneck in the regular treatment of pp-RPA happens during the construction of the pp-RPA matrix elements as the two-electron integrals involved (Eq. (3)) have to be computed on the fly. Diminishing the related computational load or completely avoiding the expensive evaluation of similar integrals has been a long-pursued goal in quantum chemistry. Different approaches have been developed in molecular orbital methods to make the study of extended systems feasible. The Mulliken approximation for the evaluation of multicenter integrals has proven to be satisfactory in many scenarios while considerably reducing the number of two-electron integrals to be handled. In this approach, the general four-center integrals are expressed as a linear combination of integrals involving only two atoms.

Alongside the Mulliken approach, other techniques have been applied in methods like the density functional tight binding (DFTB),^31,32 which have demonstrated their effectiveness for the description of a wide range of problems in physics, chemistry, and biology.^33–37 Those techniques include the monopolar truncation of the multipole expansion of the integrals and their parametrization. Furthermore, in cases where the Mulliken approach fails, the so-called on-site correction has been shown to be a successful fix.^38,39 The particle-particle DFTB (pp-DFTB) employs such techniques for the evaluation of the integrals appearing in the pp-RPA equations and uses ground state DFTB energies and orbitals as reference. Note that DFTB is even roughly three orders of magnitude faster than those relatively inexpensive DFT approaches using local or semilocal exchange-correlation (XC) functionals. Therefore, pp-DFTB is a highly efficient scheme, not only for the solution of the pp-RPA equations but also for the calculation of the reference KS orbitals and energies.

A detailed formulation and benchmark of the pp-DFTB formalism can be found elsewhere (paper in preparation). When applying the Mulliken and monopole approximations, the integrals $⟨pq|rs⟩$ can be reduced to

where the sums run over every atom in the system and $qApr$ are the Mulliken transition charges for atom A,

In Eq. (11), $cμp$ are the expansion coefficients of the molecular orbital p in a Slater-type atomic basis, and $Sμν$ are the elements of the orbital overlap matrix in such basis. (A|B) represents pairwise integrals with a Coulomb kernel, which are calculated via an interpolation formula (A|B) = U_AB(U_A,U_B,R_AB) depending on the Hubbard-like parameters U_A = (A|A) and the inter-atomic distance R_AB. The parameters U_A are computed for every atom type in advance by numerical integration and stored in a file that is read during the calculation. Thus, the evaluation of every integral is practically instantaneous. If the on-site correction is applied, additional one-center terms appear, which are also calculated in advance for every element and read on the fly. In this manner, the on-site-corrected scheme does not introduce any significant further computational cost.

The delocalized orbital problem for CT systems, which will be discussed later, can be overcome with the use of range-separated hybrid functionals in DFTB^40,41 (this method will be hereafter referred to as LC-DFTB) which has been recently implemented and benchmarked. This formalism employs a Baer-Neuhauser-Livshits (BNL)-like XC functional^42,43 with Yukawa range separation^44–46 of the Coulomb potential into short and long range contributions,

Here ω is the range separation parameter, which is set to $ω=0.3$ in this work.

LC-DFTB is based on the second-order expansion of the local part of the XC functional in the density matrix fluctuations $Δ𝐏$ around a reference density matrix. The reference is chosen to be the superposition of atomic density matrices, which are obtained from LC-DFT calculations. As basis set we use compressed atomic orbitals, obtained from pseudo-atomic DFT calculations with an additional parabolic potential (r/r₀)² for confinement purposes. The confinement radius, r₀, is chosen for every element to yield an optimal transferability. An independent confinement radius is employed for the compression of the input atomic density. There are thus two independent confinement radii per element in LC-DFTB. In a recent benchmark of the method,⁴¹ these parameters were taken from the well established parametrization set mio and no systematic optimization was conducted. Here, the set of compression radii were fitted to reproduce the CT excitation energies of the donor-acceptor systems, Benzene-TCNE and Toluene-TCNE. We then found a good transferability to the Stein set of CT complexes.⁹ In Table I, we summarize the compression radii employed in this work.

By applying the variational principle to the LC-DFTB total energy, one obtains the following Kohn-Sham equations:

with

whose solution for the N − 2 electron system provides the reference orbital energies $ϵi$ and orbital coefficients $cνi$ for the pp-DFTB equation. In Eq. (14), the zeroth-order Hamiltonian, H⁰, is evaluated at the reference density and treated in the two-center approximation. This quantity along with the overlap matrix elements $Sμν$ is tabulated as a function of the inter-atomic distance. $γlr$ and $γfr$ are, respectively, long-range and full-range two-center integrals, which are computed using analytical formulas depending on the distance between the two centers and atomic parameters calculated at the LC-DFT level of theory. Eq. (13) is solved self-consistently with an iterative update of the density matrix.

We will first carry out a qualitative analysis based on model systems. Then we will present qualitative results for the CT set by Stein et al.⁹ All pp-RPA calculations were performed with the QM4D package.⁴⁷ We performed a comprehensive study on the effect of different functionals including PBE,⁴⁸ B3LYP,^49,50 CAMB3LYP,⁵¹ RCAMB3LYP,⁵² and HF. All calculations used Cartesian cc-pVDZ basis set.⁵³ The pp-DFTB method was implemented in a development version of the DFTB+ code.⁵⁴ 

Since the pp-RPA calculation starts with an N-2 system, a good description of the electronic structure for this N-2 system is very important. This is not a challenge for most non-CT systems. However, for CT complexes, which contain two differentiated subsystems — a donor and an acceptor, their corresponding N-2 systems are more complicated and need further discussion.

Here we restrict ourselves to systems with singlet or triplet ground states and limit the bond order between the donor and the acceptor to be no larger than 1. As is shown in Figure 1, we divide the CT systems into two major categories. The first one includes singly bonded systems such as $H2$⁠, LiH, and NaCl. At a short donor-acceptor distance, their neutral species can be either covalent or ionic, whereas at a long distance, they separate into two radicals, forming a strongly correlated singlet system. The second category includes non-bonded systems, such as the famous CT models $He2$⁠, $C2H4$$⋯$$C2F4$⁠, and acene$⋯$tetracyanoethylene (TCNE). Their subsystems interact weakly through van der Waals interactions.

For the systems in the first category, although the neutral system may be strongly correlated, their N-2 picture is relatively simple. They can mostly be thought of as species that lose the two electrons that form the single bond and left with two cations. Many of these cations are closed-shell species, and therefore when the distance between the donor and the acceptor is large, the cations are much easier to describe with DFT than the neutral strongly correlated diradical system. We call them Type I CT systems. For the systems in the second category, further discussion is needed. In their N-2 picture, the two electrons can be both removed from the donor subsystem, a case that usually happens when the second ionization potential of the donor is even smaller than the first ionization potential of the acceptor. This includes simple atomic cases such as HeBe and HeMg, as well as much larger molecular cases such as anthracene⋯TCNE. We call them Type II CT systems. It is also possible that one electron is removed from the donor and the other from the acceptor. We call them Type III systems, and this type also includes many model CT systems such as $He2$ and $C2H4$$⋯$$C2F4$⁠. In these Type III systems, the singlet N-2 species become strongly correlated when the two subsystems are far away from each other.

It is worth noting that we can also do categorization based on where the electrons are removed from. Then there are two categories of electron removal: one with both electrons removed from the same subsystem and forming a neutral-dication pair (Type II), and one with electrons removed from different subsystems and forming two cations, which can be either strongly correlated (Type III) or not (Type I).

Also note that when the two electrons are removed from different subsystems, $Li2$⁠, for example, the disassociation profile for its N-2 system will be very different from that of the normal neutral systems. As a result of the Coulomb repulsion between cations, the ground state energy of this reference N-2 system follows a decreasing 1/R trend in the long distance, indicating the instability of this N-2 system.

A direct consequence of the two-electron removal is the shift of orbital energies. In a neutral system, when the donor and the acceptor are far apart with barely any orbital overlap, they can be thought of as isolated systems without any influence on each other. However, this is no longer true for the N-2 systems. For Type II systems, if we ignore the higher-order polarization effect, the neutral subsystem has little impact on the orbital energies or the orbital shapes of the dication. By contrast, the dication subsystem behaves like an external +2 point charge to the neutral subsystem. This point charge simply lowers all the orbital energies of the neutral subsystem by 2/R, where R is the separation distance.

Similarly, for Type I and Type III systems, if all orbitals are localized, then both charged subsystems act like a +1 point charge to each other, systematically lowering the orbital energies of all the orbitals in the entire system by 1/R.

Here we will adopt the single-pole approximation to analyze the reason why the pp-RPA is capable of describing the asymptotic −1/R trend for CT excitations.

We first look at the Type I N-2 systems. Their neutral ground state can be recovered by adding one electron each to the empty lowest unoccupied molecular orbitals (LUMOs) of the donor and the acceptor, while one simple CT state can be reached by adding them both to the LUMO of the acceptor. We denote the LUMOs for the donor and the acceptor as p and q, respectively. Then using the spin-adapted pp-RPA equation and considering only the diagonal elements concerning the ground state and CT excited states, we obtain the transition energies to the ground states (both singlet and triplet),

and the transition energy to the CT state

Then the CT excitation energies are

Both $ϵp$ and $ϵq$ display a −1/R trend. The first integral

is a regular electron integral within the acceptor subsystem, and therefore it does not depend on R. The second integral

is the Coulomb repulsion between two orbitals that are localized on the donor and the acceptor, respectively and therefore has a 1/R trend. The third integral

decreases to zero quickly as a result of the vanishing orbital overlap between p and q. Taking all these terms into consideration, the net trend for the CT excitation is

3211

We then look at Type II N-2 systems. Their neutral ground state can be recovered by adding two electrons back to the LUMO of the donor subsystem, and the simple CT state can be achieved by adding one electron to the LUMO of the donor and another one to the LUMO of the acceptor. We also denote these two orbitals for donor and acceptor as p and q, respectively. The single-pole approximation on the spin-adapted pp-RPA equation shows that the transition energy to the ground state is

and the transition energies to the CT singlet and triplet excited states are

Then the CT excitation energies are

Orbital q belongs to the neutral acceptor subsystem and therefore $ϵq$ has a −2/R trend due to the influence of the dication donor subsystem, whereas the orbital energy $ϵp$ almost keeps constant when R increases. Similar to the earlier discussion, the integral $⟨pq|pq⟩$ has a 1/R trend. $⟨pq|qp⟩$ vanishes quickly and $⟨pp|pp⟩$ keeps constant. Then the CT excitation energies of Type II systems have a trend of

For Type III systems, the situation becomes much more complicated. When removing two electrons, the singlet ground state becomes a strongly correlated system, just like the stretched $H2$⁠. The single reference canonical molecular orbitals often delocalize onto two subsystems. This difficult delocalized picture will be discussed later. Now we simply assume a localized picture, for example, from a broken-symmetry calculation with one β electron localizing on the donor and the other α electron on the acceptor. We denote these two singly occupied molecular orbitals on donor and acceptor as p and q, respectively. Then the ground state can be recovered by adding one α electron to p and one β electron to q. A simple S_z conserving CT state is reached by adding the β electron to q and the α electron to the local orbital of the acceptor whose energy is slightly higher than q, which we denote as r. Therefore this is a CT excitation from p to r orbital. With the spin-separated $αβ,αβ$ equation, the transition energy to ground state and CT excited states are

and the excitation energy is

Orbitals q and r are in the same acceptor subsystem, while p is in the donor subsystems. Therefore, the integral $⟨qr|qr⟩$ keeps constant while $⟨qp|qp⟩$ has a 1/R trend. The −1/R trend in orbital energies $ϵr$ and $ϵp$ cancels each other and therefore for this type of CT systems, the CT excitation still has a −1/R trend.

All in all, if we neglect the possible strong correlation and orbital delocalization problem in Type III systems, using a single-pole analysis, we can conclude that the description of CT excitations by pp-RPA is always asymptotically correct.

It should be noted that in a more realistic multiple-pole picture, as long as the configurations involved lead to the same CT type, the asymptotic behavior remains the same. This is because the configuration interactions among those CT configurations keep constant at long distances and have no asymptotic R dependence. Therefore only the −1/R trend from each single pole picture between the ground state configuration and every CT configuration remains. However, we need to remind the readers that the current implementation of pp-RPA is only limited to HOMO excitations. Therefore, those CT excitation configurations excited from below the HOMO of the donor are not included in the multiple-pole picture.

It is also worth noting that the discussions above are independent of any functional that is employed in the N-2 SCF procedure. Therefore, the correct −1/R asymptotic behavior always holds for any functional reference.

Now let us come to the orbital delocalization problem for Type III systems. This problem is by far the greatest challenge for pp-RPA in describing CT excitations. The delocalized orbital problem is closely related to the challenging strongly correlated systems. For example, the simplest molecule of this type is $He2$⁠, and after removing two electrons, it resembles the famous challenging $H2$ molecule.⁵⁵ In the KS molecular orbital picture, in which all molecular orbitals diagonalize the one-body effective Hamiltonian, both the HOMO and the LUMO delocalize onto two atomic centers with a spin-restricted calculation. Only in the dissociation limit can they be transformed into localized orbitals and still be eigenvectors of the Hamiltonian. Note that this delocalized picture does not frequently occur for orbitals below the HOMO or above the LUMO. For example, Figure 2 shows the delocalized frontier orbitals and localized non-frontier orbitals for N-2 o-xylene⋯TCNE. Even in challenging symmetric cases such as $He2$⁠, a small perturbative external field can help localize these non-frontier orbitals.

In this scenario, because the LUMO is delocalized, all previous discussions on the asymptotic trend based on localized molecular orbitals are no longer valid. Although it is still true that adding two electrons to the delocalized LUMO recovers the ground state, yet adding one to the LUMO and another one to a localized higher virtual orbital gives a state that has only about half electron transferred. Besides the problem with the number of electrons transferred, for these N-2 systems, it is often difficult to converge the DFT SCF iterations, which is a prerequisite for all pp-RPA calculations. Because of these two obstacles, in our previous work,²² we tried a different reference named HF*, in which we first perform a SCF calculation for neutral systems to obtain localized orbitals, and then we manually remove two electrons by changing the occupation number. This approach essentially changes the challenging Type III systems to an easier Type II system, and it gives the correct −1/R trend for $He2$ and $C2H4$$⋯$$C2F4$ CT model systems.²² Furthermore, it is also shown to be equivalent to configuration interaction singles and doubles (CISD) for $H2$⁠. However, it currently can only be derived using the equation-of-motion method and only combines with HF references. Another way of overcoming these problems is to adopt the high-spin triplet (S_z = 1) or broken-symmetry singlet as initial references, but they are temporarily not pursued due to potential spin symmetry breaking and spin incompleteness. The most desirable way is to obtain localized molecular orbitals with a spin-restricted calculation. This might be possible if we can combine the pp-RPA with constrained DFT,¹⁴ thus changing Type III systems to Type II.

The delocalization problem is also especially challenging for pp-DFTB calculations. In this scheme, the reference orbitals inherit this problem as DFTB is based on semi-local DFT functionals. For the Stein set of CT complexes, the HOMO and LUMO are also delocalized within DFTB. As a result, both the donor and acceptor subsystems have nearly +1 charge, thus leading to a complicated Type III system. As mentioned in Section II B, this issue can be overcome by adopting reference orbitals calculated with the LC-DFTB method. LC-DFTB calculations for the N-2 electron CT systems produce more localized HOMO and LUMO compared to DFTB results. This leads to a simple Type II CT system, where the donor is +2 charged and the acceptor is neutral.

Another challenge is rooted in the error of DFT frontier orbital energies. In principle, within the Kohn-Sham or the generalized Kohn-Sham formulations, the HOMO energy from the exact functional should equal the ionization potential, while the LUMO energy equals the electron affinity.^56,57 However, for currently used density functionals, this exact relation is violated. Especially, the LUMO energy is always greatly underestimated. This quantitative error can bring about qualitative problems in some CT systems. The Be⋯He system is one example. In principle, the two electrons should be both removed from the Be atom because the second ionization potential of Be is smaller than the first ionization potential of He. Therefore, it should form a Type II N-2 system without any convergence issue or orbital delocalization problem. However, in practice, most density functionals cannot converge the SCF calculation because they predict the LUMO of Be²⁺, which should correspond to the second ionization potential of Be, to be lower in energy than the HOMO of He, which corresponds to the first ionization potential of He.

Figure 3 shows the orbital energies as a function of the donor-acceptor distance for the N-2 Be⋯He system. It can be seen that the HOMO of the whole system, which is also the HOMO of He, is affected by the +2 charge on Be and its energy has a −2/R trend. In principle, its asymptotic limit should equal the first ionization potential of He. The LUMO of the whole system is also the LUMO of Be²⁺ and its energy does not change with R. This energy should equal the second ionization potential of Be. The B3LYP functional underestimates the LUMO energy but overestimates HOMO energy of the whole complex, leading to an incorrect energy crossing when the distance is about 7.5 Å. Before the crossing, the two electrons are both removed from Be, while after that it will in principle predict an unphysical Type III system and make both orbitals singly occupied. However, the change of occupation numbers leads to orbital re-optimization and a new set of orbital and orbital energies, finally making the SCF very difficult to converge. Even if with some convergence trick, we can converge the SCF calculation here and obtain Type II systems, the final solution is non-Aufbau with the HOMO energy higher than the LUMO energy. This is not quite likely to be physical for a ground state. In contrast to B3LYP, HF performs much better in this system with an accurate prediction of HOMO and LUMO energies both at short and long donor-acceptor distances. However, it should be noted that this good accuracy for HF is because the system is small and simple. In more general and more complicated systems, it is already known that HF tends to overestimate HOMO-LUMO gaps.

Therefore, the underestimated HOMO-LUMO gaps by commonly used density functionals may qualitatively change an easy Type II CT system into a challenging Type III system or into an unphysical Type II system with a non-Aufbau occupation.

Regarding pp-DFTB, the challenge associated with the correct description of frontier orbital energies can be addressed by employing the LC-DFTB formalism as the reference. LC-DFTB has been shown to deliver fair accuracy for the calculation of HOMO and LUMO energies, with results comparable to first-principle approaches using the BNL functional and outperforming DFT calculations based on local and hybrid functionals.⁴¹ For example, the mean unsigned deviation (MUD) of the negative of the LC-DFTB HOMO energies of a set of organic molecules from their experimental ionization potential was reported to be 0.50 eV.⁴¹ This is comparable to the BNL/cc-pVTZ value of 0.30 eV and substantially smaller than the B3LYP/cc-pVTZ and PBE/cc-pVTZ values of 2.04 and 2.87 eV, respectively. When using the confinement radii summarized in Table I, the MUD slightly increases with a value of 0.62 eV, but remains much smaller than the results for the ab initio semilocal approach. HOMO-LUMO gaps within LC-DFTB are also satisfactory in general. For the same set of molecules, we obtained a MUD with respect to LC-DFT (BNL) of 1.69 eV, considerably lower than for PBE (5.21 eV) or B3LYP (3.83 eV). For the +2 charged CT complexes under investigation, HOMO-LUMO gaps are also accurately predicted within LC-DFTB for both short and long donor-acceptor distances.

We use the CT set by Stein et al.⁹ for our testing purpose. All the systems in this set have a strong electron withdrawing molecule — TCNE, and an electron donating molecule ranging from derivatives of benzene to derivatives of anthracene.

We find that for the complexes containing benzene or its derivative, the N-2 systems suffer severely from delocalized frontier orbital problems with DFT references and therefore become Type III systems. The N-2 benzene⋯TCNE is difficult even for SCF convergence. For those systems that can be converged, the degree of delocalization varies with the choice of the XC functional. In general, the less HF exchange, the more delocalized the frontier orbitals are. As is discussed earlier, the delocalized orbitals lead to a picture with less than one electron transferred.

By contrast, the systems including molecules from the anthracene family do not have the orbital delocalization problem, probably due to the relatively larger electron donating capacity of anthracene. Therefore, the two electrons are both removed from the anthracene or its derivatives. The systems containing naphthalene is in the middle — some density functionals with large HF exchange, for example, CAMB3LYP and RCAMB3LYP, predict localized orbitals while those with small HF exchange, for example, B3LYP, predict delocalized orbitals.

We perform our pp-RPA calculations with PBE, B3LYP, CAMB3LYP, RCAMB3LYP, and HF references, which roughly have an increasing portion of HF exchange (RCAMB3LYP has a small portion of HF exchange in the short distance but more than 100% HF exchange correction in the long range). We also compute the CT excitation energies using pp-DFTB, except for the two complexes that contain the Cl element, because the mio set of parameters does not include the Cl element. We first present the results for Type II naphthalene and the anthracene family. Afterwards we will discuss the results for the challenging benzene family. The available experimental reference values were obtained either in the gas phase or in the liquid phase. In order to get a complete set of references for both gas phase and liquid phase, as is suggested by Stein et al.,⁹ we add 0.32 eV to liquid phase results to obtain approximate gas phase results or subtract 0.32 eV from gas phase results to obtain approximate liquid phase results.

The pp-RPA and its approximation pp-DFTB results are listed in Table II and then plotted in Figure 4. It can be seen that with the increasing portion of HF exchange, there is a systematic trend for the results obtained with pp-RPA. The SCF procedure for the N-2 systems frequently fails to converge when the PBE functional is employed, but when convergence is reached, the pp-RPA calculations give the lowest excitation energies. The pp-RPA calculations with B3LYP reference give larger excitation energies than those with PBE, but still greatly underestimate the excitation energies compared with both gas phase and liquid phase references. Compared to liquid phase results, the pp-RPA-CAMB3LYP only slightly underestimates the results by 0.14 eV, while pp-RPA with RCAMB3LYP and HF references overestimate the results by 0.12 eV and 0.15 eV, respectively. However, if we compare them with the gas phase reference, which is assumed to be 0.32 eV higher than the liquid phase reference, then all the pp-RPA variants underestimate the results and give larger errors. The best performance is achieved with the RCAMB3LYP and HF references, for which the excitation energies are underestimated by 0.20 eV and 0.19 eV, respectively. As can be seen from Table II and Figure 4, the overall performance of pp-DFTB with LC-DFTB is comparable to that of pp-RPA with the RCAMB3LYP and HF references, that is, those with the greatest percentage of exact exchange. Similarly, the CT energies are slightly overestimated compared to the liquid phase reference whereas they are in average underestimated with respect to the gas phase reference. Importantly, the set of parameters employed is shown to be transferable across the table, which covers different donor compounds and ranges of CT energies. pp-DFTB has the advantage of being very efficient, which is particularly useful when tackling systems with problematic convergence. In such cases, for which a large number of SCF steps are required, the calculation takes only few seconds to converge.

We choose the RCAMB3LYP results as one of the best pp-RPA-DFT results, as well as the pp-DFTB results to compare with other methods based on TDDFT. The results are plotted in Figure 5. It can be seen that the TD tuned-BNL gives excellent results that closely match the gas phase reference. The pp-RPA and pp-DFTB perform slightly worse and give results between the gas phase reference and the liquid phase reference. Results from TDDFT with B3LYP and CAMB3LYP functionals are even lower.

It should also be noted that if we look at the relative values between different systems, TDDFT with conventional functionals performs badly. The results with TDB3LYP is especially strange — when the reference value is relatively large, it predicts a small number and when the reference value is relatively small, it predicts a relatively large number. In Figure 5, if we compare its trend line with that of the liquid phase reference, there seems to be a mirror lying at around 1.5 eV. This interesting phenomenon is surprising to us and worth further investigation. However, since this paper is mainly about the pp-RPA, we will temporarily ignore this strange phenomenon by TDB3LYP and only focus on the results by pp-RPA.

In order to understand better the underestimated CT values by pp-RPA, we designed and performed a detailed analysis based on the structure of the pp-RPA working equation. The procedure is sketched in Figure 6.

The most direct approximation (step 1) is to simply look at the orbital energy contribution from the orbitals involved. The ground state is recovered by adding two electrons to the LUMO of the N-2 system, while the lowest CT state is obtained by adding one to the LUMO and another one to the lowest virtual orbitals localized on TCNE, which is often LUMO+1 or LUMO+2. For convenience, we now denote the LUMO as L₀, while the localized virtual orbital on TCNE as L_n. Therefore, the orbital energy contribution to the CT excitation energy is

It is nothing but the difference of two virtual orbital energies. Then we can move one step further (step 2) and perform the single-pole approximation, in which we also include the kernel contribution but exclude the coupling between the two configurations. As we analyzed before, the approximate CT energy is

The third step is to consider the coupling, but only within the $2×2$ single-pole approximation space. The coupling element is

Note here, the coupling can be zero if the symmetries of L₀ and L_n are different. The fourth step is to move to pp-TDA, which considers all the coupling among those configurations involving two-electron addition in unoccupied orbitals, and therefore it potentially becomes a multi-pole picture. Finally, the last step is the full pp-RPA, which also includes some contribution from two-electron addition in occupied orbitals.

Figure 7 shows the results for the molecules in benzene ((a) and (b)), naphthalene (c), and anthracene ((d)–(f)) families. From the first step results, we can see that for all molecules, it always holds that the more HF exchange, the larger the L_n − L₀ gap. This trend is the same as the well studied HOMO-LUMO gap.

When we move one step further and carry out the single-pole approximation, the CT energies decrease dramatically. Nonetheless, for molecules in the anthracene family, this decrease is almost the same for all functionals, and therefore, the lines are almost parallel. By contrast, for those molecules in the benzene family, the amount of decrease varies depending on the functional. Interestingly, the larger the original L_n − L₀ gap, the larger decrease in this second step, finally almost reversing the order of the results for different functionals. The cause of this behavior lies in the shape of frontier orbitals. As we discussed earlier, all systems containing anthracenes are Type II systems with localized frontier orbitals that are almost the same for different functionals. However, the systems containing benzenes are Type III systems, and their frontier orbitals differ from each other due to the varied degree of orbital delocalization for different functionals. Therefore, the kernel correction also varies greatly for different functionals. The results for naphthalene (Figure 7(c)) lie between those of the benzene family and the anthracene family. The pp-RPA-B3LYP predicts delocalized frontier orbitals while the pp-RPA with CAMB3LYP, RCAMB3LYP, and HF references predicts localized orbitals. Therefore, for naphthalene, the behavior of pp-RPA-B3LYP in this step is very different from the rest.

Next, when we go to step 3 and take the coupling effect into consideration, we can see that in the benzene and naphthalene families, the correction is almost zero, indicating the coupling element is almost zero due to a symmetry mismatch. The correction for anthracenes is much larger, but still considerably smaller than the kernel correction. This correction in this step is also very similar for different functional references.

For the fourth and fifth steps, the CT excitation energies become further larger. However, even if we add up the effects from the last three steps, they are usually far from canceling the correction effect of the second step. For the last two steps, the corrections for different functionals are also similar.

If we look at the five steps together, for anthracene, the lines for different functionals almost go parallel with each other, indicating that the differences among different functionals practically originate from the initial orbital energy difference. Therefore, the orbital energy difference between L_n and L₀ directly affects the final results for Type II systems. By contrast, for Type III systems such as benzenes, because of the varying frontier orbital, the kernel also makes a difference.

Could the errors in these orbital energies be the source of the underestimation in CT excitations? Unfortunately, although the L_n − L₀ gap plays the most important role, current DFT studies cannot give a clear interpretation for the virtual orbital energies beyond LUMO. Therefore, we still do not have a way of judging whether the L_n − L₀ gap is overestimated or underestimated. This is another theoretical challenge arising from N-2 systems. Nonetheless, from the systematically underestimated data in this study, it is very likely that currently used density functionals also underestimate the L_n − L₀ gaps. It is possible that even HF might underestimate them as well. However, it should also be specially noted that even an exact L_n − L₀ gap does not necessarily lead to the most accurate pp-RPA results, because the pp-RPA kernel is already an approximation within the TDDFT-P theory.³⁰ At the current stage, a good excitation result will rely on some error cancellation between the approximate orbital energies and the pp-RPA kernel. This error cancellation may as well be the reason why the pp-RPA with a reference with a small portion of HF exchange is even better than those with a large HF exchange for the benzene family.

In summary, as long as the orbitals are localized, the pp-RPA can always give the asymptotically correct −1/R trend for CT excitations because of a concerted contribution from the orbital energies and the kernel. This asymptotic trend is independent of any functional references. In practical CT systems, the pp-RPA performs better than TDDFT with traditional hybrid or range-separated functionals, but it also underestimates the CT excitation energies. It still cannot compete with the TD tuned-BNL approach. However, our pp-RPA approach is universal in that it uses the same functional for all systems. The pp-DFTB can greatly reduce the computational cost and yields CT results comparable with the best pp-RPA results. There remain three major challenges in the description of CT excitations for pp-RPA. First, the delocalized orbitals in combination with strongly correlated systems often account for the SCF convergence failure and a picture involving the transfer of roughly half an electron. Second, the currently used density functionals often underestimate the LUMO energy of the dication species, turning the easy Type II systems into difficult Type III systems with delocalized orbitals. These two challenges can be overcome with a suitable range-separated hybrid functional within pp-DFTB. Third, the performance of pp-RPA greatly depends on the energy difference between the L_n and L₀ orbitals, but the meaning of this energy difference is still unclear, and therefore it calls for further theoretical investigations. It is worth noting that these three challenges are also relevant to TDDFT, which also sometimes suffers from SCF convergence failure, underestimated HOMO-LUMO gaps, and dependence of higher virtual orbitals, although not in the context of charge transfer excitations we discuss here. Nevertheless, a fix to these three issues will benefit both pp-RPA and the conventional TDDFT. In summary, although facing some challenges, the pp-RPA and pp-DFTB are promising general approaches for describing charge transfer excitations.

Y.Y. and W.Y. appreciate the support from National Science Foundation (Grant No. CHE-1362927). D.Z. appreciates the support as part of 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.