The attraction density functional theory (DFT) has for electronic structure theory is that it is easier to do computationally than *ab initio*, correlated wavefunction methods, due to its effective one-particle structure. On the contrary, *ab initio* theorists insist on the ability to converge to the right answer in appropriate limits, but this requires a treatment of the reduced two-particle density matrix. DFT avoids that by appealing to an “existence” theorem (not a constructive one) that all its effects are subsummed into a DFT functional of the one-particle density. However, the existence of thousands of DFT functionals emphasizes that there is no satisfactory way to systematically improve the Kohn-Sham (KS) version as most changes in parameterization or formulation seldom lead to a new functional that is genuinely better than others. Some researchers in the DFT community try to address this issue by imposing conditions rigorously derived from exact DFT considerations, but to date, no one has shown how this route will ever lead to converged results even for the ground state, much less for all the other electronic states obtained from time-dependent DFT that are critically important for chemistry. On the contrary, coupled-cluster (CC) theory and its equation-of-motion extensions provide rigorous results for both that KS-DFT methods are attempting to emulate. How to use them and their exact formal properties to tie CC theory to an effective one-particle form is the target of this perspective. This route addresses the devil’s triangle of KS-DFT problems: the one-particle spectrum, self-interaction, and the integer discontinuity.

## INTRODUCTION

Coupled-cluster (CC) theory and its applications will offer the underlying framework for what this perspective will say about density functional theory (DFT), so before embarking on that task, a few words about *ab initio* correlated theory and coupled-cluster theory (CC) will be useful. One might say that our foray into DFT and Kohn-Sham (KS-DFT) is geared toward reducing CC theory to an exact one-particle form, where the correlation described by a CC wavefunction might be replaced by adding a “correlation potential” to define an effective one-particle operator whose eigenfunctions are “correlated molecular orbitals” and whose eigenvalues are principal Ips for occupied levels and electron affinities for at least the lumo. KS-DFT is one vehicle for trying to do that.

The most fundamental objective of *ab initio* correlated theory is that the results should converge to the right answer in appropriate limits (unlike virtually all DFT methods). Thus, the theory can, and should be, “predictive,” providing results in the absence of experiment or resolving experimental discrepancies. The goal of any *ab initio*, correlated calculation should be the full configuration interaction (FCI) that has to provide the exact answer for the nonrelativistic Schrödinger equation in a basis set. Thus, for the small systems where it can be done, a comparison provides a litmus test for the accuracy of a correlated method.^{1} I do not think there is much argument that today, the CC route toward correlated results—even with a single reference (SR)—and its operationally single reference equation-of-motion (EOM-CC) extensions^{2} offer the best answers for the largest number of problems encountered. This is often the case even when conventional thinking would say that the problem being addressed requires a multireference (MR) (determinantal) starting point, as a correct zeroth-order approximation. The reason, of course, is that the single reference (SR-CC) gets to the FCI very quickly, and the FCI does not care what single determinant reference is used in its generation, as long as it overlaps the exact wavefunction. Thus, SR-CC answers at some reasonable truncation of the cluster operator will frequently transcend those obtained from obviously multireference (MR) situations, either MR-CC or configuration interaction (MR-CI). One does have to use enough cluster operators in SR-CC and its EOM-CC extensions to give the method the freedom it needs to provide a correct description, but unlike MR methods, where the choice of reference space introduces another potential uncertainly that has to be addressed to ensure a converged calculation, there are virtually no decisions for the user in CC and EOM-CC applications beyond the level of correlation and basis chosen.

Of course, CI single or MR does not scale appropriately with the number of particles in the system, i.e., it is not (size-) extensive^{3} as it contains inappropriate unlinked diagrams in it. They are also responsible for its variational bound, so such a bound is inconsistent with extensivity, until the FCI is reached. The fact that a FCI for a formally infinite system such as a polymer or crystal is impossible also means that only “many-body” methods can be used to describe electron-correlation in such systems. For electron correlation, *any* CI approximation for an infinite-system has zero overlap with the correct wavefunction.

The SR-CC equations are now well known,^{1}

where H is the Hamiltonian for the system and for T itself, one has projections of the Schrödinger equation onto the space of single, double, triple, … excitations, $Q\u2282\Phi ia,\Phi ijab,\Phi ijkabc,\u2026$,

to derive the nonlinear CC equations that can have products of up to four Ts in them. $H\xaf$ naturally deals with the one- and two-particle reduced density matrices, required to describe the total energy and other two-particle properties.

The extension to EOM-CC is also well-known. On top of $H\xaf$ which describes the ground state (or suitable reference state), all other states are available as target states. For example, for electronic excited (EE) states, the right- and left-hand operators, R (L), consist of particle-hole excitations {*a*^{†}*i*, *a*^{†}*ib*^{†}*j*, …}, while if the target state is for an ionized state (IP-EOM-CC), we create a “hole” in the target state by having one more hole annihilation operator, *i*, than creation operators, *a*^{†}. If we want electron attached states (EA-EOM-CC), our operator creates one more electron than in the reference state, {*a*^{†}, *a*^{†}*b*^{†} *j*, *a*^{†}*b*^{†}*c*^{†}*ij*, … }, and so forth for double IP and double EA (DIP/DEA-EOM). Then, one simply diagonalizes the (non-Hermitian) matrix,

to get the spectrum of states of interest in an operationally single reference way.

Note that unlike most quantum chemical methods, there is also a seamless connection between electronic and ionized or electron attached states. This transition is important to our time-dependent DFT (TDDFT) eigenvalue theorem used in the Quantum Theory Project (QTP) functionals discussed later.

Certain obvious multideterminantal situations such as open-shell singlets in excited states (EE), or two electrons in two orbitals in all possible ways (DIP, DEA), can be readily addressed by such methods.^{100} Although much progress has been made in genuinely multideterminantal reference CC methods,^{4} for this perspective, the focus is on providing SR-CC/EOM-CC quality results from some kind of effective one-particle KS-DFT theory. This objective seems to be the next best option to making KS-DFT predictive.

Now, we see that we have a Hamiltonian in CC, $H\xaf$, that has any single determinant as its eigenstate, with eigenvalue E_{CC}. Furthermore, important to this perspective is that the IP/EA-EOM-CC equations can also be put into the form of *an effective one-particle equation,* for an h^{eff} operator consisting of the Fock operator, f(1), that contains the kinetic energy operator, the electron-nuclear potential, the Fock effective one-particle operator, J-K, and a CC correlation potential, Σ_{CC}(1),

with

The quantity, $\u2211CC(1)$, is a *frequency independent* form of the Dyson self-energy,^{101} i.e., a one-particle correlation potential that can be derived from manipulations involving the IP/EA-EOM-CC equations to assume an effective, one-particle form^{5,6} to augment the Fock operator, *f*. Its right (and left) eigenfunctions are correlated molecular orbitals defining a correlated orbital theory (COT) whose common occupied eigenvalues, *ϵ*_{p}*,* correspond to “principal Ips.” This means the energy required to remove an electron from the system for *any* of the occupied orbitals with no “shake-up” effects. The actual eigenvalues will be the same as those computed from IP-EOM-CC, whose FCI limit is well-defined.

At least some of the unoccupied ones will have an eigenvalue that corresponds to an electron affinity. But since most EAs lie in the continuum (within a square-integrable basis set, a spectrum of EAs could be obtained from EA-EOM-CC as a discretization of the continuum), but without due regard to the real and imaginary parts of the spectrum, the real nature of the unoccupied orbitals and their EAs would be unclear. So only the lowest EA states are considered in this perspective. But like the IP-EOM-CC result, the EA-EOM-CC has to also provide the full CI result in the limit. Thus, this COT one-particle theory is a model for KS-DFT, with particular relevance to how one can potentially provide CC accuracy within an effective one-particle theory: the avowed goal of our work in KS-DFT. In DFT the exact density is given by the square of n-orbitals, but not in COT.

As we explore DFT and its KS realization, the above discussion of CC/EOM-CC will color all developments. In particular, the rigorous one-particle spectrum obtained from IP/EA-EOM-CC is fundamental to electronic structure theory providing a COT or *quantitative* molecular orbital theory (QMOT). Yet KS-DFT approximations prior to those discussed here, seldom correctly describe it.^{7} There are also two-particle effects (ultimately arising from the reduced two-particle density matrix) that besides determining the total ground state energy pertain to electronic excitation energies and to operators such as spin (S^{2}). Yet, any direct consideration of the two-matrix is bypassed in KS-DFT, a goal of any one-particle, effective correlated method.

The following will focus on having the correct one-particle system and the correct density from just n-orbitals, with the two-particle effects arising from the associated functionals.

## ADVENTURE I. HF-DFT

Our first foray into DFT was to explore whether a typical DFT functional, such as Becke, Lee-Yang-Parr (BLYP), would give us accurate correlation energy corrections when added on top of a Hartree-Fock (HF) density. From a molecular viewpoint, doing a HF calculation is easy, so if BLYP or another DFT functional would provide a correlation correction once the HF density is dropped into it—that is what is meant by HF-DFT—then this is an extremely easy way to get a correlation energy, compared to CCSD(T), for example.

To answer this question, Oliphant added a numerical integration procedure to ACES II. He did a number of tests on how best to do the numerical integrations required by DFT functionals.^{8} He also did side-by-side comparisons between HF-DFT and CCSD(T) in large (TZ2P) basis sets. Others had partly explored the same idea starting with Clementi.^{9–14} We attempted to be thorough, treating energies, geometries, vibrational frequencies, and dipole moments and used large enough bases, TZ2P, to provide reliable conclusions. Our results were surprising. For the set of molecules we considered, the correlation corrections from HF-DFT were in excellent agreement with CCSD(T). See Table I.

. | SCF . | HF-DFT . | CCSD(T) . |
---|---|---|---|

Error in bond length | |||

Mean | 0.022 | 0.005 | 0.005 |

Maximum | 0.083 | 0.022 | 0.016 |

Error in bond angle | |||

Mean | 2.7 | 1.7 | 1.9 |

Maximum | 9.3 | 5.7 | 8.9 |

Error in harmonic frequency | |||

Mean | 144 | 40 | 31 |

Maximum | 394 | 166 | 136 |

Error in atomization energy | |||

Mean | 84.5 | 3.7 | 11.5 |

Maximum | 161.0 | 16.3 | 20.7 |

Error in dipole moment | |||

Mean | 0.29 | 0.14 | 0.10 |

Maximum | 1.18 | 0.55 | 0.29 |

% Error in dipole moment | |||

Mean (%) | 52 | 17 | 9 |

Maximum (%) | 241 | 65 | 20 |

. | SCF . | HF-DFT . | CCSD(T) . |
---|---|---|---|

Error in bond length | |||

Mean | 0.022 | 0.005 | 0.005 |

Maximum | 0.083 | 0.022 | 0.016 |

Error in bond angle | |||

Mean | 2.7 | 1.7 | 1.9 |

Maximum | 9.3 | 5.7 | 8.9 |

Error in harmonic frequency | |||

Mean | 144 | 40 | 31 |

Maximum | 394 | 166 | 136 |

Error in atomization energy | |||

Mean | 84.5 | 3.7 | 11.5 |

Maximum | 161.0 | 16.3 | 20.7 |

Error in dipole moment | |||

Mean | 0.29 | 0.14 | 0.10 |

Maximum | 1.18 | 0.55 | 0.29 |

% Error in dipole moment | |||

Mean (%) | 52 | 17 | 9 |

Maximum (%) | 241 | 65 | 20 |

However, at first sight why should this be? Becke’s local exchange in BLYP^{15} was determined by fitting the nonlocal HF exchange energy for the inert gas atoms to a local form. Clearly, the intent of that procedure was to represent exchange, so “B” should have nothing to do with correlation. “LYP,” on the other hand, takes the correlation expression of Colle and Salvetti (CS)^{16} and represents it as a density functional.^{17} As the CS expression is meant to be the correlation correction on top of HF, logically, one would drop the HF density into that to get correlation corrections. But doing so is a disaster. Hence, regardless of its origins, “BLYP” itself is a “correlation functional” as defined by giving the right correlation corrections on top of a HF density. We also observed this for the Perdew, Burke, Ernzerhof (PBE) functional, so this is not peculiar to BLYP.

Now, it appears many think that a local exchange has correlation effects of left-right type,^{18,19} perhaps especially important for bond breaking. The idea of a double hybrid (DH),^{20} where one uses a weighted combination of local and nonlocal exchange (NLX), and local and nonlocal correlation like that given by second-order many-body perturbation theory (MBPT2), would seem to pertain to this HF-DFT observation.

But there is more to the HF-DFT story. A follow-up paper by Sekino *et al.*^{21} addressed the nature of the HF-DFT density matrix. In other words, this very simple correlated method provides all the components of an honest correlated method such as CC or MBPT theory beyond the energy correction itself, like the correlated density matrix and its natural orbitals. The mathematical structure should be especially interesting as HF-DFT has undergone something of a renaissance.^{22}

When we were interested in providing a hierarchy of correlated approximations for rate constants from the simplest possible, HF-DFT, to high-level CC theory, Verma took on the task of using HF-DFT to compute transition states (TSs) and activation barriers.^{23} One would think that is a demanding task for such a simple, nonvariational method. Yet, once again, it worked amazingly well. Verma, Perera, and I developed the nonvariational equations to analytically obtain transition states and applied them to the set of reactions provided by Truhlar.^{24,25} In that paper, we also redid the TSs using our ΛCCSD(T) method^{26} [an improvement over the perturbative CCSD(T) by using the left-hand eigenvector, Λ of the CC functional], whose results at consistently determined TS geometries were accurate to ∼1 kcal/mol. The mean absolute deviation (MAD) results obtained from HF-DFT for H-transfer reactions were accurate to 2.5 kcal/mol, compared to 6.3 for those using BLYP in a self-consistent manner, with HF at 12.3 and ΛCCSD(T) at 1.4. For the nonhydrogen transfer reactions at fixed quadratic CI (QCISD) geometries, the results for HF-DFT with BLYP were 2.67 and 3.17 with PBE. Those results were competitive with most KS-DFT studies done at that time and remain impressive for such a simple scheme.

As discussed below, in the intervening years between our study of HF-DFT density matrices and its application to activation barriers, we had explored a number of other issues in DFT and its KS form. These caused us to focus on the KS potential more than the functional as from our viewpoint, it is the *weak link* in KS-DFT applications.

This leads to the message made in the above HF-DFT paper. In our opinion, the reason HF-DFT works well compared to normal, self-consistent KS-DFT, is that

Its potential, the Fock operator, does not suffer from the KS potential’s self-consistency problem, as the latter, Vxc[ρ] = δExc[ρ]/δρ(1), is only correct to first-order in the density, ρ(1). That is, while a normal (convex) variational functional Exc[ρ] is correct to second-order in ρ(1), its functional derivative is only correct to first order (note that the full meaning of “order” depends on the relevant perturbation theory, discussed later). This is the KS-DFT analog of variational wavefunction theory as the energy is correct to second order, but the wavefunction itself, only to first order. Thus, the self-consistent treatment of the KS-DFT equations that requires that Vxc[ρ] be updated each iteration to give new KS orbitals will suffer from this first-order error unless precautions are taken. Consequently, the converged density, ρ(1) = Σ|φ

_{i}|^{2}, will be flawed, often being much worse than the non-self-consistent route taken in HF-DFT. Thus, a poor potential can cause the so-called “density error” in KS-DFT.^{22}Whereas the self-consistent potential Vxc[ρ] suffers from the self-interaction error (SIE) that permeates virtually all KS-DFT calculations, in HF-DFT, the SIE is removed by the simple expedient of using the Fock potential to generate the HF density.

Third, no Vxc[ρ] in virtually any widely used KS-DFT approximation provides a meaningful one-particle spectrum from the obtained KS orbitals.

^{7}On the contrary, the Fock operator does, as is well-known from Koopmans’ theorem.

Hence, we believe that progress in KS-DFT requires (1) eliminating self-interaction; (2) that the equations provide a correct zeroth-order one-particle spectrum for the ground state KS calculation; and (3) that the potential be vetted by the ionization eigenvalue theorem. The latter is the precaution required to avoid suffering from its first-order error in self-consistent KS calculations. This will be discussed in depth after our second adventure.

## ADVENTURE II. EXACT LOCAL EXCHANGE, OPTIMIZED EFFECTIVE POTENTIAL (OEPx)

The optimized effective potential (OEP) approach has been long known in KS-DFT^{27–31} as it allows one to use the exchange functional,

where the spin orbitals i and j are both occupied in the KS determinant to define an exchange potential. EXX is the definition of “exact” exchange in KS-DFT as it is the HF exchange energy written in terms of KS orbitals. The next step is to turn this EXX functional into a potential, V_{x}^{OEP} [ρ(1)], to enable self-consistent KS calculations. This uses the OEP procedure which means define the *local* potential in a basis set. The most important contribution of our work led by Ivanov^{32} demonstrates how to do OEPx in a Gaussian basis set, the S algorithm. In the adjacent paper, a slightly different approach was offered by Görling,^{33} the V algorithm. Prior work with OEPx used either numerical grids or plane waves. Yet if this method is to be widely used in quantum chemistry, one has to be able to do it in a Gaussian basis set.

The relevant equations can be obtained from the functional derivative chain rule,

The first term is straightforward. The second can be handled with coupled-perturbed KS theory. But the last depends on the inverse of the density response function, X(1, 2)^{−1}. The density response is

Note that X(1, 2) depends only on the orbital products of occupied, φ_{i}(1), and unoccupied orbitals, φ_{a}(1),^{102} contrary to the more general form used by Krieger-Lee and Iafrate (KLI).^{31}

The V_{x}^{OEP} is a local, multiplicative form of the exchange potential unlike the usual exchange operator (K) of HF theory. The latter is usually called nonlocal exchange (NLX). Because V_{x}^{OEP} is local, one can plot it and analyze how it behaves.

Examples taken from our OEPx work^{34–36} are shown in Figs. 1 and 2. One should note that the OEPx potential has the correct long-range (−1/r) behavior. It also shows some shell structure lacking in the Becke exchange and those from LDA and Slater. But what happens when used in calculations?

The total energy from OEPx, meaning the expectation value of the KS determinant, ⟨KS|H|KS⟩, is quite close to HF as one would expect but is always above HF because this OEPx calculation has an extra *local potential* constraint compared to that of HF. For Ar, e.g., HF is −526.818 a.u. compared to 526.812, for OEPx.^{32}

We know that the OEPx homo energy should be exactly that from HF.^{37} This “shift” is included as a boundary condition in the OEPx potentials shown in Figs. 1 and 2. Nonetheless, we can observe the results. For N_{2}, where the Koopmans’ HF approximation has the wrong order for the σ and π orbitals, OEPx gets it right, giving 17.16 eV compared to experiment, 15.58 eV. For H_{2}O and CO, the values are 13.71 and 13.76 (HF) (exp 12.62) and 15.03 vs 15.07 (HF) (exp 14.01). If an LDA correlation potential is added, the orbital energies are far worse compared to known IPs.^{34} (When OEP2sc correlation is added after our third adventure, the results agree well with experiment attesting to the essential role of correlation.)

Another aspect of this work is that there is no need for any numerical integration. All manipulations are done directly in the Gaussian basis just as in quantum chemistry.

Finally, importantly, there is no (one-electron) self-interaction problem in OEPx since the functional is the correct exchange that cancels the Coulomb term, E_{J} = ∫∫ρ(1)ρ(2)/|**r**_{1} – **r**_{2}|, when some orbital, {φ_{i}}, is the same. Thus, the $\Sigma iiiii$ (the occupied orbitals) defines the SIE, and this number amounts to hartrees even for small molecules, so its elimination is a necessity if KS-DFT is to be developed to give the “right answer for the right reason” (this phrase is associated with Ernest Davidson). OEPx is a critical step toward doing KS-DFT in as *ab initio* a way as possible.

A follow-up paper by Hirata,^{35} with some help from the mathematician Talman of OEP fame, enabled us to prove an important theorem that the OEPx potential is unique in the exact case and although ill-posed for a finite number of functions can be made to be unique by projecting the equation and the exchange potential upon the accessible function space of the kernel.

The next step was to use OEPx in TDDFT. This required that we evaluate the OEPx kernel. Again led by Hirata, my group formulated this subject to the adiabatic approximation that all terms including the kernel are frequency independent and implemented it to make the first TDDFT study using OEPx.^{36}

The numerical results are interesting. Basically, the TDOEP results are a little better than those from TDHF, using experiment and STEOM-CCSD^{38} as the reference. The latter provides MAD from experiment as 0.22 eV for Be, 0.13 for N_{2}, and 0.09 for water. TDOEP provides 0.46, 0.92, and 0.83, respectively, while TDHF gives 0.54, 1.50, and 1.12. Clearly, there is a difference for valence excited states, as in N_{2}, compared to the Rydberg excited states in H_{2}O. But perhaps the most interesting numerical result is that the unoccupied orbital energies, {*ε*_{a}}, determined from the *local* OEPx potential, provide an amazingly good, Hückel approximation for excited states as $\omega ia\u2248\epsilon iOEPx\u2212\epsilon aOEPx$. Needless to say, when using *ε*_{a} from a HF calculation that approximates the electron affinity, this would be a terrible approximation except for charge-transfer examples.

Further work by Ivanov *et al.*^{34} addressed how best to solve the OEPx equations in a basis set. Hirata *et al.* further used TDOEPx to describe dynamic polarizabilities,^{39} demonstrating that unlike the incorrect behavior encountered with normal KS-DFT approximations,^{40} OEPx correctly describes them. Bokhan also applied OEPx to static hyperpolarizabilities.^{41}

## ADVENTURE III. *AB INITIO* DFT. OEP CORRELATION POTENTIALS (OEP2sc)

Following the above, one could seriously focus on the rigorous inclusion of electron correlation in KS-DFT. At the time, there was said to be no correlation potential that could properly tie to the OEPx potential. Clearly, none of the standard ones like LYP or Vosko, Wilks, Nussair (VWN) used in LDA would work. There had been some mostly failed efforts^{42,43,103} to use a second-order MBPT orbital-dependent functional as a source for a correlation potential, but this turns out to be very difficult. It was first realized numerically by Grabowski and my group.^{44} We simply took the functional derivative of MBPT2 analogous to how we did so for EXX, and got some results for what has been called Goerling-Levy perturbation theory.^{45} We also introduced the terminology “*ab initio DFT*,” or perhaps better, “*ab initio KS-DFT*,” which meant that this approach was, in principle, capable of converging to the exact, FCI result, unlike other KS-DFT efforts. That is, if one successfully takes the functional derivative (or its equivalent in our density based method, discussed later) of the FCI orbital-dependent energy expression to define a Vc^{(FCI)}, then along with Vx^{(OEP)}coupled through the KS iterations (note the FCI is invariant to orbital choice), one should be able to obtain the FCI result.

The term “*ab initio*” also emphasizes that we are using the orbital dependent EXX and MBPT2 functionals known from wavefunction theory and since there is no need for numerical integration, the correct reference result would be the FCI in the chosen basis as in our case, the same basis was used for the orbitals and the potential instead of some auxiliary basis for the latter.

But in practice there was a fatal flaw. All efforts including our first had assumed that the KS problem should provide the logical zeroth-order approximation for a converging perturbation series. After all, this defines the “adiabatic path” from the noninteracting to the interacting Hamiltonian. But when that interacting Hamiltonian is analyzed,^{46–49} it is apparent that its perturbation on top of the KS unperturbed problem has diagonal and other terms in it, which means it is likely to be too large to offer an adequate low-order MBPT2 approximation. As such, one is led to divergence in trying to get the OEP correlation potential. Many excuses for this inability to obtain a meaningful Vc^{(OEP)}were made,^{103,104} but the simple fact is that the obvious PT based on the noninteracting KS Hamiltonian as the “unperturbed” problem and MBPT2 as the correlation functional is very difficult to converge because of the too-large perturbation in the latter. These failures arise from undue reliance on the adiabatic path. Our path differs. We first look at the interacting Hamiltonian and then tie its density to that for the KS determinant by insisting that the correlated density corrections vanish to some order of perturbation theory. This procedure uniquely defines the correct V_{XC} to that order (up to a constant) that makes the solution of the KS equations provide the correlated density as the KS one.

The need for realistic correlation potentials is apparent from Fig. 3, for Ar,^{50} showing a comparison between PBE and OEP2sc, defined below.

Writing the interacting Hamiltonian in second-quantization,

With the two-electron integrals ⟨pq||rs⟩ in the antisymmetrized, 1, 2, 12 form, and expanding the Fock matrix elements in terms of occupied orbital {i, j, k, l} and unoccupied orbital {a, b, c, d}, while {p, q, r, s} denotes any orbital,

Then, introducing the fact that f = h + J − K = h^{s} − K − Vxc, and expressing H_{N} in terms of the h^{s} (KS) eigenfunctions, where ⟨p|h^{s}|q⟩ = δ_{pq}ϵ_{p}, the Hamiltonian becomes

W indicates the two-particle term.

When evaluating this interacting Hamiltonian in MBPT, there are several choices that could be made for the separation into *H* = *H*_{0} + *V*. If one chooses to define the unperturbed Hamiltonian by

the Kohn-Sham choice, then the perturbation for that choice, *V*^{KS}, has to retain *all other* terms including the diagonal elements, ⟨*p*|*K* + *Vxc*|*p*⟩.

Having such diagonal terms in *V*^{KS} plus the remainder of the terms makes the perturbation unnecessarily large, frequently leading to divergence when trying to solve the OEP equations with a MBPT2 orbital dependent functional.^{103,104,49} We call this correlation potential OEP2ks. It is the appropriate one for Goerling-Levy perturbation theory so that approach will generally fail.

From many years of MBPT results, it is now well-established that even though some other choices of *H*_{0} like that due to Epstein-Nesbet appear to potentially have some attractive features, the Moeller-Plesset (MP) choice is properly orbital invariant (occ-occ, virt-virt) unlike the Epstein-Nesbet choice and tends to offer the best balance between the numerator and denominator in perturbation theory to generally ensure good convergence with orders. Pathological cases can diverge, of course, but even those are usually easily resumed with Pade’ approximants^{51} along with more general resummation techniques such as those in CC theory.

The MP choice in generalized MBPT (GMBPT),^{52} meaning for any choice of orbitals such as the KS ones used here, is to employ the full Fock operator involving occupied-occupied and unoccupied-unoccupied terms. This means the first three terms generated by *f* in Eq. (13). This defines H_{0}^{GMBPT}. The occupied-virtual term, *f*_{ia}, has a different role, so it, along with W, defines the “correlation” perturbation, V^{GMBPT}. Besides observed numerical performance, this H_{0}^{GMBPT} choice is invariant to rotations within the occupied and virtual spaces and is self-interaction free as the Coulomb and exchange terms are treated consistently. We consider such requirements to be mandatory. They also facilitate computations. Its invariance offers a convenient way to handle this H_{0}^{GMBPT} in practice by making a semicanonical (SC) orbital transformation to put it into the form

in terms of the rotated orbitals, $p\u0303$. Except for the new orbitals, the *V*^{SC} perturbation remains the same as V^{GMBPT}, ensuring that the “correlation” perturbation remains manageable in MBPT2 and higher approximations.^{105} Henceforth, we will assume semicanonical orbitals and suppress the ∼ in any equations.

With this choice of “correlation” perturbation, we derive the OEP equations from the condition that the correct *correlated* density should be given by the KS determinant, $\Phi KS,$ and ρ(1) = Σ|ϕ_{i}|^{2}. Consequently, we want

where *δ*(*r* − *r*′) is the Dirac delta function and the usual resolvent is *R*_{0} = (*E*_{0} − *H*_{0})^{−1}*Q* If we insist that this density vanish, then we have the OEPx condition,

With the help of the density-density response function, we can write this as

to determine *V*_{X}._{.}

The MBPT2 orbital dependent functional makes it possible to also obtain the $VCOEP2sc$ potential with care from the OEP equations^{46–48} properly using singular value decomposition for the generalized inverse of X(1, 2). For correlation, the key equation is to make the second-order in V^{SC} density vanish,

The latter is turned into a real space equation for *L*(1) from which the exchange and correlation potential is defined by the solution to the response equation,^{47,48}

then

The correlation potential arises from MBPT2, so it is second-order, while the exchange potential arises in first-order. Of course, once both potentials are inserted into self-consistent KS cycles, one simply goes to convergence losing any order identity to the final exchange-correlation potential. The KS orbital product is then a set of n-occupied, correlated orbitals that constitute the KS determinant and represent the corrected density as $\rho 1=\u2211i|\phi i|2$.

Note that these equations come from the “density condition,” Eqs. (17) and (20), not any functional derivative. The relationship between the functional derivative and our much-preferred density based approach is detailed in the Appendix of Ref. 53. Whereas any functional derivative re-introduces the adiabatic path, the density condition adds appropriate infinite summations that enable the theory to overcome the divergence encountered with the KS choice of H_{0}. Thus, the density condition naturally avoids the adiabatic path as it introduces Eq. (16) for H_{0} that leads to a resummation of terms that avoids the divergence. This is the proper second-order perturbation theory, not GL perturbation theory!

This correlation potential answered the question of what potential properly fits the $VXOEP$ exchange potential. These correlation potentials were the first self-consistent potentials successfully obtained from an orbital dependent correlation functional such as MBPT2, linearized coupled-cluster singles and doubles (LCCSD)^{49} (also called CCPT2), CCD, and RPA.^{7} See Fig. 4 for RPA and CCD.

Using the RPA ground state energy^{7} as a correlated energy functional is recommended from the fluctuation-dissipation theorem of DFT,^{54,55} but it is usually used as an add-on rather than as a self-consistent potential. But as we showed, our OEP approach provides it self-consistently as in Fig. 4, as an approximation to CCD that sums only the direct-ring diagrams. See the work of Heβelmann^{56} for further studies of RPA. In principle, these OEP methods could be used with any other correlation methods including those from MCSCF, to MR-CI, to MR-CC, etc.

Grabowski, Teale, Smiga, and I made a study of the convergence of exchange and correlation potentials with the level of theory from MBPT2 to LCCSD, CCSD, and CCSD(T).^{57} The idea was to demonstrate that just as the results from *ab initio* correlated methods converge systematically (but not variationally) to the right answer, the potentials generated do, too. Hence, it is valid to talk about converged exchange and correlation potentials. Examples for Ar are shown in Figs. 5 and 6.

These potentials are obtained in two complementary ways. The *ab initio DFT* ones directly use an orbital dependent correlation functional such as MBPT2 and LCCSD (CCPT2) in the OEP procedure that defines the potential and then uses them in KS iterations until converged. If this were done at the CCSD or beyond level, the well-known insensitivity of CCSD to the orbital choice would make it difficult to gain anything further from the KS iterations. Hence, the CCSD and CCSD(T) potentials are obtained from an inversion procedure, asking what one-electron potential produces the computed CCSD or CCSD(T) densities. *There is no energy functional.* However, when done both ways, as in the MP2 (MBPT2) potentials in the figures, the results are entirely consistent. When quantum Monte Carlo (QMC) results are also available,^{58} there is agreement with those, too. The fact that these potentials often have no correspondence with those used in DFT for correlation such as LDA and Slater-Vosko-Wilks-Nussair (SVWN) (see Figs. 5 and 6) raises an important issue about whether any of the correlation potentials actually used in the vast majority of KS-DFT applications really have “genuine” correlation in them. Instead, the two potentials together, *pseudo*-exchange and *pseudo*-correlation, might show some changes in results that can be in the direction of correlation corrections, but until one is controlled to focus on the other, it is very hard to say that there is “real correlation” included in the calculation. The only viable control at this time seems to be to treat the exchange at the OEPx level or 100% HF (NLX) and then ask what happens to the correlation.

In the correlated OEP case (OEP2sc), correlation is immediately apparent in the orbital energies and other correlated properties, including the density.^{46,47} But the uncertainty in those used in the other functionals camouflages the real role of the correlation potential in any application. Hence, if one chooses a KS-DFT calculation because one thinks it is “correlated,” that argument needs to be reconsidered. Also, considering that pseudo-exchange and pseudo-correlation have the same sign but scale differently, the oft quoted argument that the sum of the two is more accurate than the separate parts is hard to justify. If (T_{coor} = −E_{coor}) could be included in the analysis, there might be some potential cancellation.

An effort has been made to fit the observed exchange and correlation potentials from Fig. 5, as a function of the density, i.e., V_{XC}[ρ(r)], instead of Vxc(r) shown in the plots.^{59} By making the potential density dependent, one adheres to the DFT structure and would hope that Vxc could then have a larger domain of application. But this is quite difficult to do, as this paper discusses. However, one might say that is an essential step if a truly seamless connection between WFT and DFT is to be achieved. The attempt to use the adiabatic connection between weak and strong correlation as a route to develop better correlation potentials is commendable.^{60}

One final example pertains to weak interactions where correlation is essential. *Ab initio DFT* provides relatively correct weak interactions even though the potentials are local and consequently do not overlap at its van-der Waals minimum for the simple reason that the functional, in this case MBPT2, is nonlocal. So despite the potentials used to generate the orbitals, as Lotrich *et al.* showed,^{61} OEP-2sc will still show a proper minimum for weak interactions that is usually a little better than CCSD. This should be contrasted with the “tack-on” approaches used to try to fix KS-DFT for weak interactions. A satisfactory theory should be able to describe such interactions as a standard part of the theory.

In more recent work from the Grabowski group, instead of using MBPT2 as the orbital dependent, correlation functional, they chose to use the so-called spin-scaled version (SCS-OEP) of second-order that weights the different spin components differently than in MBPT2 itself. This resulted in some successes for improving correlation potentials that compare favorably to the original OEP2ks choice^{62} and for the KS eigenvalues that according to our theorem^{47} (discussed below) should be good approximations to principal Ips.^{63}

In the case of MBPT2 and LCCSD, the failures of the OEP2ks potentials are apparent compared to either the wavefunction generated potentials and radial distribution functions, and those from *ab initio DFT.* It is particularly apparent when applied to Be, where the quasidegeneracy of the 2s and 2p orbitals causes slower convergence, and OEP2ks cannot be converged at all. In Fig. 7 is shown the first iteration of an ultimately diverging attempt to solve for the OEP2ks potential. For other systems, the lack of convergence was shown by Bokhan^{64} for some 12 atom and small molecule examples, where 7 diverged. Schweigert’s work^{65} showed that 19 out of 35 molecules in the G2 set did not converge. Yet, OEP2-sc does okay for Be,^{49} while its infinite order CCPT2 (LCCD) [in Fig. 7 D-MBPT(∞)] does quite well in overcoming the effects of the quasidegeneracy. Both have no problem with divergence for all the other examples investigated.

Armed with the newly derived OEPx and OEP2sc local potentials, Schweigert and I asked some interesting questions about mixing nonlocal HF (NLX) into them.^{53,65} For example, is there anything more to gain from making the transition to what is called KS-DFT’s generalized form (GKS)? In practice, nearly all KS-DFT calculations today have some admixture of NLX, such as B3LYP and PBE0, all the range separated ones such as CAM-B3LYP, where the degree of NLX varies with distance, etc. One nice thing about OEPx and NLX is that any combination of the two still avoids any self-interaction error as the EXX functional is the same for both potentials, so one can look at V_{XC} = (1 − κ)V_{X}^{OEPx} + κK^{NLX} + V_{C}^{OEPsc}, as a function of κ.

The first interesting observation is that as little as 6% NLX is sufficient to overcome the divergence issue encountered with the KS unperturbed Hamiltonian. In fact, adding 37% NLX to the Be example results in the potentials generated in Fig. 8. But otherwise, the addition of NLX to the OEPx result shows no useful improvement for most properties and the optimum amount is also system specific. Furthermore, we concluded “that the semicanonical modification of the KS-based MBPT renders the presence of the nonlocal exchange redundant, which permits retaining the formal consistency of purely local operators in KS DFT. However, if the KS eigen-spectrum is concerned, by adjusting the amount of the nonlocal exchange in the hybrid Hamiltonian, in this case by a fixed 50%, the correlated orbital energy spectrum provides an accurate approximation to the principal IPs, including those for core electrons.” The results are shown in Table II.

. | I_{k}
. | −ε_{k} − I_{k}
. | ΔI_{k}
. | ||
---|---|---|---|---|---|

. | Expt. . | HF + PT2H(100) . | EXX + PT2SC . | $12$(HF + EXX) + PT2SC . | IP-EOM CCSD . |

H_{2}O | |||||

1b_{1} | 12.62 | −0.44 | −0.13 | −0.18 | 0.02 |

1a_{1} | 14.74 | −0.40 | −0.20 | −0.20 | 0.08 |

1b_{2} | 18.51 | −0.18 | −0.11 | −0.05 | 0.57 |

2a_{1} | 32.61 | 2.54 | −1.93 | 0.40 | 0.23 |

1a_{1} | 539.7 | 20.9 | −22.6 | −0.9 | 1.2 |

CO | |||||

5σ | 14.01 | 0.41 | 0.38 | 0.42 | 0.21 |

1π | 16.91 | −0.57 | 0.16 | −0.13 | 0.21 |

4σ | 19.72 | 0.89 | −0.32 | 0.36 | 0.11 |

3σ | 38.30 | 1.34 | −3.68 | −1.07 | −1.29 |

2σ | 296.2 | 13.8 | −16.9 | −1.6 | 0.8 |

1σ | 542.6 | 21.5 | −22.4 | −0.5 | 1.1 |

. | I_{k}
. | −ε_{k} − I_{k}
. | ΔI_{k}
. | ||
---|---|---|---|---|---|

. | Expt. . | HF + PT2H(100) . | EXX + PT2SC . | $12$(HF + EXX) + PT2SC . | IP-EOM CCSD . |

H_{2}O | |||||

1b_{1} | 12.62 | −0.44 | −0.13 | −0.18 | 0.02 |

1a_{1} | 14.74 | −0.40 | −0.20 | −0.20 | 0.08 |

1b_{2} | 18.51 | −0.18 | −0.11 | −0.05 | 0.57 |

2a_{1} | 32.61 | 2.54 | −1.93 | 0.40 | 0.23 |

1a_{1} | 539.7 | 20.9 | −22.6 | −0.9 | 1.2 |

CO | |||||

5σ | 14.01 | 0.41 | 0.38 | 0.42 | 0.21 |

1π | 16.91 | −0.57 | 0.16 | −0.13 | 0.21 |

4σ | 19.72 | 0.89 | −0.32 | 0.36 | 0.11 |

3σ | 38.30 | 1.34 | −3.68 | −1.07 | −1.29 |

2σ | 296.2 | 13.8 | −16.9 | −1.6 | 0.8 |

1σ | 542.6 | 21.5 | −22.4 | −0.5 | 1.1 |

I trust that another pertinent observation during our third adventure was to eliminate the dogma surrounding the role of the KS orbitals and their energies.^{106} On the contrary, for some time, Baerends has called attention to the fact that the KS orbitals are better than the HF ones for various purposes, including having a degree of orbital relaxation compared to the normal Koopmans’ HF orbitals. In particular, inserting them into the Dyson equation provided evidence that their eigenvalues are properly associated with ionization potentials.^{66}

We took this a step further. By the simple expedient of looking at the solutions of the adiabatic TDDFT equations when an electron is excited into the continuum, that is, “ionized,” it has to follow that the ground state KS solution has to have eigenvalues attached to each canonical orbital that correspond to the principal Ips.^{67} This pertains to “all orbitals” not just the homo where the usual exact DFT condition deduced from the long-range behavior of the density forces its eigenvalue to be the first IP. The proof of this IP-eigenvalue theorem is that by simply adding a zero Gaussian function, or z ∼ exp(−10^{70})r^{2},^{47} into the basis set used in the adiabatic TDDFT calculation of excited states, the IPs will fall out as eigenvalues that correspond to excitations into the continuum. If the KS eigenvalues do not have this property—and they do not in the vast majority of KS-DFT applications^{7}—then the KS-DFT method is not “consistent.”^{68} Consistency means four things: (1) The functional provides reasonably well described total energy properties. (2) The Vxc potential is vetted to be accurate by the eigenvalue-IP condition. (3) The orbital eigenvalues properly tie to the TDDFT excitation energies, making them correct, principal ionization potentials with a seamless connection to the bound electronic states of a Rydberg series. (4) The potential is the functional derivative of the energy functional or its density based equivalent as we use, not just some arbitrary potential used to give a set of orbitals.

The second of these is particularly important. The DFT tautology means start with an energy functional that is convex in ρ, meaning that its variational KS solution is correct to second-order in ρ. By taking the functional derivative to define the Sham potential, $\delta E\rho \delta \rho =vS$, *v*_{S} is formally only accurate to first order in ρ. As discussed above, this is the KS-DFT analog to wavefunction theory. To ensure the potential is accurate requires that it be constrained. We propose to do that by insisting that for any occupied orbital, ϕ_{m},

up to a small threshold. We believe this is a critical condition to provide improved solutions from KS-DFT. First, and foremost, one now gets the qualitatively right one-particle spectrum that other inconsistent KS-DFT methods will not provide. Second, the Ip condition is tantamount to ensuring that the method is largely self-interaction free. In terms of our *ab initio* COT Eq. (7), as long as

self-interaction is mitigated as there is none on the right side of this equation. Third, those two properties along with the integer discontinuity constitute the “Devil’s Triangle” of what is wrong with KS-DFT. The devil is chipping away at “Jacob’s ladder.”^{69} But, nonetheless, I believe all three corners of the triangle can be accommodated mostly by exploiting this eigenvalue IP condition.

We will see some of the consequences in Adventure IV.

## ADVENTURE IV. THE QTP CONSISTENT SET OF FUNCTIONALS

The basic idea behind all the QTP functionals, or better, reparameterization of existing functionals, is the ionization eigenvalue theorem proven from adiabatic TDDFT. To satisfy this, the orbital energies in a ground state KS-DFT are expected to be the principal IPs to within a small threshold. This is a kind of “Koopmans” theorem for KS-DFT, but I avoid using that term since in DFT, we do not have “wavefunctions,” and as Nesbet points out,^{70} the real point of Koopmans’ theorem in HF theory is that after ionizing or attaching an electron, the wavefunction is “optimum” (although not orbitally relaxed), with the orbital energy difference simply being a consequence. This interpretation permits the same logic to be used in deriving the extended Koopmans’ theorem (EKT).^{71}

On the contrary, the TDDFT proof does not depend on anything more than the Runge-Gross theorem based exclusively on the density.^{72} In unpublished work,^{73} we further derive the Ip-eigenvalue theorem for KS-DFT from the Janak theorem. Hence, there is no doubt that this IP-eigenvalue theorem is a valid constraint and it should be imposed on new functionals and used to update old functionals to enable KS-DFT to significantly “increase its capability.”

The Frontier article we wrote on the QTP functionals summarizes much of this,^{74} so for this perspective, I will try to avoid duplication while limiting myself to some essentials.

Let us look at a few examples. We took the Coulomb-attenuated, CAM-B3LYP functional,^{75} as an example of a relatively standard, extensively parameterized, range-separated approximation and minimally reparameterized it to define CAM-QTP(00).^{76} It has the form

with range separated NLX, B88 local exchange, and LYP and VWN correlation potentials. The *μ* parameter controls the degree of range separation, while α and β weight the local and nonlocal exchange and γ, the mix of the two “correlation” potentials. When α + β = 1 as in QTP(01) and (02), we eliminate one parameter. In QTP(02), we also dispense with γ choosing just LYP as the *pseudo-*correlation term.

To test the idea of our IP-eigenvalue theorem, our first attempt, QTP(00), simply adjusted these four parameters to reproduce the five principal IPs of water. We were able to do that to a threshold of ∼0.8 eV with a set of parameters that are far from unique. From the homo down, the errors compared to experiment are 0.08, 0.18, 0.34, 0.77, and 0.69. The surprise is that when applied to 42 other arbitrary, closed shell molecules, their one-particle, vertical IP spectrum agreed to the same ∼0.8 threshold! As the one-particle spectrum for all 42 molecules makes sense, the results show that the choice of the water molecule in the parameterization is not an important factor. In other words, there is some “universality” to this IP-eigenvalue theorem.

The second surprise is that when applied to a *total energy property* such as transition states and activation barriers, this very simple, minimally parameterized method using B3LYP geometries for the TS gave MAD errors of 2.3 kcal/mol for H transfer reactions and 3.02 for non-H transfer reactions, competitive with most KS-DFT attempts.

The third surprise was that this method also gave quite reliable results for the hard to describe core IPs of light-atom molecules and even their TDDFT excitation spectrum.^{77} The errors for both are ∼1 eV. So somehow, this simple parameterization accounts for the very large, core orbital relaxation on the *orbital eigenvalues* that KS-DFT using CAM-QTP(00) produces. There is no delta KS calculation involved. Hence, we conclude that enforcing the Ip-eigenvalue theorem has a significant role on improving the performance of KS-DFT calculations at no additional cost. In other words, unlike *ab initio DFT* that is one model built upon the solution of the OEP equations, the QTP functionals take no more time nor trouble to apply than any standard, range separated KS-DFT calculation. In any program that does CAM-B3LYP, one can simply change from its parameters to ours to obtain the QTP family of results, for better or worse.

The QTP functionals fit the parameters in Eq. (28),

CAMQTP00

Fitted to 5 IPs of water.

CAMQTP01

Fitted to 4 valence IPs of water, 34 excitation states of 4 molecules (ethylene, formaldehyde, acetaldehyde, and acetone), and total atomization energies of G2-1 test set.

CAMQTP02

Fit 2 parameters to homo-lumo condition that makes IP^{*}(OH^{−}) = EA(OH), IP(H_{2}O) = EA(H_{2}O^{+}) as close as possible.

The actual parameter values are shown in Table III, compared to the original CAM-B3LYP values.

. | Exchange . | Correlation . | ||
---|---|---|---|---|

Name . | α
. | β
. | μ
. | ν
. |

CAMB3LYP | 0.19 | 0.46 | 0.33 | 0.81 |

CAMQTP00 | 0.54 | 0.37 | 0.29 | 0.80 |

CAMQTP01 | 0.23 | 0.77 | 0.31 | 0.80 |

CAMQTP02 | 0.28 | 0.72 | 0.335 | 1.00 |

. | Exchange . | Correlation . | ||
---|---|---|---|---|

Name . | α
. | β
. | μ
. | ν
. |

CAMB3LYP | 0.19 | 0.46 | 0.33 | 0.81 |

CAMQTP00 | 0.54 | 0.37 | 0.29 | 0.80 |

CAMQTP01 | 0.23 | 0.77 | 0.31 | 0.80 |

CAMQTP02 | 0.28 | 0.72 | 0.335 | 1.00 |

We also saw that CAM-QTP(00) did not do everything well. It is inferior for geometry predictions and heats of atomization. It also can be improved for excitation energies. All three are much improved by Jin’s work to define CAM-QTP(01),^{78} e.g., still retaining the IP condition for the four valence IPs of water but sacrificing the core IP, while adding some other info. QTP(01) is also a bit better for IPs (see Fig. 9 and Table IV). Table IV includes results from the extended Koopmans’ theorem (EKT) and QTP(17),^{79} the latter unlike all others of the QTP family is not range-separated but is instead a global hybrid, based on BLYP, but targeted to core ionization effects. Nonetheless, if it is true that imposing the IP-eigenvalue theorem will help any KS-DFT approximation, this is an initial illustration. However, it is notably poorer for the homo IP.

. | Valence orbitals . | HOMO . | ||||
---|---|---|---|---|---|---|

. | No. . | MUD . | MAD . | No. . | MUD . | MAD . |

LDA—Expt. | 401 | −4.66 | 4.66 | 63 | −4.33 | 4.33 |

B3LYP—Expt. | 401 | −3.26 | 3.26 | 63 | −3.23 | 3.23 |

CAM-B3LYP—Expt. | 401 | −1.37 | 1.37 | 63 | −1.50 | 1.50 |

CAM-QTP00—Expt. | 401 | 0.76 | 0.81 | 63 | 0.16 | 0.32 |

CAM-QTP01—Expt. | 401 | 0.13 | 0.40 | 63 | −0.15 | 0.30 |

CAM-QTP02—Expt. | 401 | 0.38 | 0.49 | 63 | 0.00 | 0.25 |

QTP17—Expt. | 401 | −0.18 | 0.49 | 63 | −0.78 | 0.78 |

HF—Expt. | 401 | 1.84 | 1.86 | 63 | 0.63 | 0.74 |

EKT(CCSD)—Expt. | 401 | 0.97 | 0.97 | 63 | 0.34 | 0.35 |

LDA—CCSD | 401 | −4.84 | 4.84 | 63 | −4.34 | 4.34 |

B3LYP—CCSD | 401 | −3.44 | 3.44 | 63 | −3.23 | 3.23 |

CAM-B3LYP—CCSD | 401 | −1.55 | 1.55 | 63 | −1.51 | 1.51 |

CAM-QTP00—CCSD | 401 | 0.58 | 0.65 | 63 | 0.16 | 0.32 |

CAM-QTP01—CCSD | 401 | −0.05 | 0.31 | 63 | −0.16 | 0.29 |

CAM-QTP02—CCSD | 401 | 0.20 | 0.34 | 63 | 0.00 | 0.24 |

QTP17—CCSD | 401 | −0.36 | 0.47 | 63 | −0.79 | 0.79 |

HF—CCSD | 401 | 1.66 | 1.69 | 63 | 0.63 | 0.76 |

EKT(CCSD)—CCSD | 401 | 0.79 | 0.80 | 63 | 0.34 | 0.34 |

. | Valence orbitals . | HOMO . | ||||
---|---|---|---|---|---|---|

. | No. . | MUD . | MAD . | No. . | MUD . | MAD . |

LDA—Expt. | 401 | −4.66 | 4.66 | 63 | −4.33 | 4.33 |

B3LYP—Expt. | 401 | −3.26 | 3.26 | 63 | −3.23 | 3.23 |

CAM-B3LYP—Expt. | 401 | −1.37 | 1.37 | 63 | −1.50 | 1.50 |

CAM-QTP00—Expt. | 401 | 0.76 | 0.81 | 63 | 0.16 | 0.32 |

CAM-QTP01—Expt. | 401 | 0.13 | 0.40 | 63 | −0.15 | 0.30 |

CAM-QTP02—Expt. | 401 | 0.38 | 0.49 | 63 | 0.00 | 0.25 |

QTP17—Expt. | 401 | −0.18 | 0.49 | 63 | −0.78 | 0.78 |

HF—Expt. | 401 | 1.84 | 1.86 | 63 | 0.63 | 0.74 |

EKT(CCSD)—Expt. | 401 | 0.97 | 0.97 | 63 | 0.34 | 0.35 |

LDA—CCSD | 401 | −4.84 | 4.84 | 63 | −4.34 | 4.34 |

B3LYP—CCSD | 401 | −3.44 | 3.44 | 63 | −3.23 | 3.23 |

CAM-B3LYP—CCSD | 401 | −1.55 | 1.55 | 63 | −1.51 | 1.51 |

CAM-QTP00—CCSD | 401 | 0.58 | 0.65 | 63 | 0.16 | 0.32 |

CAM-QTP01—CCSD | 401 | −0.05 | 0.31 | 63 | −0.16 | 0.29 |

CAM-QTP02—CCSD | 401 | 0.20 | 0.34 | 63 | 0.00 | 0.24 |

QTP17—CCSD | 401 | −0.36 | 0.47 | 63 | −0.79 | 0.79 |

HF—CCSD | 401 | 1.66 | 1.69 | 63 | 0.63 | 0.76 |

EKT(CCSD)—CCSD | 401 | 0.79 | 0.80 | 63 | 0.34 | 0.34 |

From Fig. 9 and Table IV, it is clear that the QTP functionals offer results that are competitive with IP-EOM-CCSD, one goal of our attempt to reduce CC to a one-particle form. The paper offers further comparisons to IP-EOM-CCSDT and CCSDTQ.^{67}

For other properties, QTP(01) provides excitation energies compared to the Caricato data set^{80} accurate to 0.29 eV compared to 0.60 for QTP(00), or 0.27 eV for EE-EOM-CCSD. Anharmonically corrected vibrational frequencies are accurate to 42 cm^{−1} compared to 97 for QTP(00), and atomization energies are accurate to 3.6 kcal/mol compared to 10.7.

The avowed intent of KS-DFT calculations is to provide an accurate approximation to the electronic density itself. Our approach to document this is to evaluate all the moments of the charge distribution from r^{−2} to r^{3}. Then, one can ask how well do these QTP parameterizations do compared to MBPT and CC theory. This is shown in Fig. 10 for light atoms. CAM-QTP(00) provides rather good agreement. Reference 81 also addresses the density for transition metal atoms where CAM-QTP(01) does better.

A different example of the benefit due to imposing the IP-eigenvalue theorem can be seen from Margraf ’s study of double hybrid (DH) KS-DFT.^{82} In this case, the mixed wavefunction, DFT functional

is composed of local, multiplicative GGA components and nonlocal HF exchange and MBPT2. The local exchange is B88, and the local correlation is LYP. The original functional of Grimme termed B2-PLYP uses an inconsistent potential

that ignores the MBPT2 part of the functional. For a number of molecules, this potential results in a mean absolute error of 1.53 eV for the homo, and 1.50 eV for the valence.

Substantial improvement is achieved by making the potential consistent. To do so, one has to add to it the OEP generated expressions for *V*^{HF} and for $VCMBPT2$. This is done in two forms: an (a) form that is the purely local OEP2sc potential and a (b) form that takes the liberty of adding any amount of NLX and local OEPx exchange together (50-50), without compromising the one-electron self-interaction,

Using the former, purely local B2-PLYP + OEPsc(a), values are 0.89 for the homo and 0.80 for the valence. OEP2sc(b) itself gives 0.53 and 0.42 and (except for H2S) provides a 0.9 eV error for the core.

But setting aside consistency between functional and potential, if one takes the inconsistent potential of Eq. (30) and fixes the two parameters *a* and *c* by invoking the IP-eigenvalue theorem for a few molecules, a = 0.60, with c = 0.35, primarily from some heats of atomization, and then compares this form of the functional, B2IP-PLYP to the original, B2-PLYP, where a = 0.53 and b = 0.27 and an extensively reparameterized version of Martin, B2GP-PLYP,^{83} who used a large set of W4 atomization energies and barrier heights that gives a = 0.65 and c = 0.36.

As shown in Fig. 11, for 13 of the data sets used in KS-DFT,^{107} the IP form is slightly better. In an independent assessment of excitation energies, B2IP-PLYP is deemed the best.^{84}

The other major problem that affects KS-DFT is the integer discontinuity caused by adding or removing electrons in KS-DFT calculations.^{85} What happens in the real world is that a straight-line should connect the results for N electrons, and N + 1, and N − 1. See Fig. 12 ^{86} for the C atom. This figure shows how the QTP parameterizations fare compared to many other approximations.

The exact result is shown by the dashed black line, while the one concave curve is HF, with the rest showing the convex “curvature” shown by most KS-DFT calculations. The red curve is CCD, now generalized to treat fractional occupation numbers,^{108,87} while brown is MBPT2. The orange, blue, and gray curves are from the QTP choices. All are closer to linear than any of the other functionals except MBPT2 and maybe, CCD. These types of plots are important to understanding much about the correct KS-DFT behavior, albeit with finite Gaussian basis sets. See Ref. 88 for a thorough study of this general behavior for electron addition and removal of fractional electrons with particular regard for the effect of finite basis sets in such curves.

If an approximate KS-DFT calculation offers the right IP and EA, then, at least the integer end-points of this figure would make the curves coincide. This offers part of the inspiration for the completely *nonempirical* CAM-QTP(02).

The IP eigenvalue theorem applies to any bound system, neutral, cation, anion, etc. Just as it applies to F^{−}, it has to also apply to F. Thus, the homo eigenvalue for F^{−} has to be the ionization potential, but the electron affinity of F arising from its lumo eigenvalue has to be the same number. This is called the homo-lumo condition. This obvious exact condition is seldom satisfied in any KS-DFT method, and even IP-EOM-CC and EA-EOM-CC will not typically provide the exact equivalence, partly due to one being a closed shell and the other an open shell causing different rates of convergence (see Table V).^{108}

Species . | EA-EOM-CCSD . | Species . | IP-EOM-CCSD . |
---|---|---|---|

B^{+} | 8.23 | B | 8.31 |

CN^{+} | 14.17 | CN | 14.11 |

CH^{+} | 10.57 | CH | 10.67 |

NO^{+} | 9.7 | NO | 9.74 |

LiH | 0.29 | LiH^{−} | 0.3 |

LiF | 0.33 | LiF^{−} | 0.33 |

F_{2} | 0.14 | F_{2}^{−} | 0.29 |

F | 3.22 | F^{−} | 3.41 |

OH | 1.61 | OH^{−} | 1.84 |

H_{2}O^{+} | 12.70 | H_{2}O | 12.68 |

Species . | EA-EOM-CCSD . | Species . | IP-EOM-CCSD . |
---|---|---|---|

B^{+} | 8.23 | B | 8.31 |

CN^{+} | 14.17 | CN | 14.11 |

CH^{+} | 10.57 | CH | 10.67 |

NO^{+} | 9.7 | NO | 9.74 |

LiH | 0.29 | LiH^{−} | 0.3 |

LiF | 0.33 | LiF^{−} | 0.33 |

F_{2} | 0.14 | F_{2}^{−} | 0.29 |

F | 3.22 | F^{−} | 3.41 |

OH | 1.61 | OH^{−} | 1.84 |

H_{2}O^{+} | 12.70 | H_{2}O | 12.68 |

^{a}

Only the valence electrons are included in the active space of these calculations that used unrestricted wavefunctions as a reference for open-shell cases. The molecular calculations were performed at the experimental geometries of the neutral systems.

But, nonetheless, in appropriate correlation and basis set limits, this should be the case. This is the requirement Haiduke and I impose to define the two parameters in QTP(02).^{89} CAM-BLYP is the underlying functional restricting the correlation to just LYP, eliminating γ. We also insist that α + β = 1, eliminating another parameter. Clearly, this condition should help with the integer discontinuity/curvature error (note that there is still a difference between the left side and right side associated with the discontinuity). Also, this procedure indirectly addresses EAs as a special case of IPs avoiding any “fitting” to any EA information. The systems we actually use are the OH^{−}. OH^{•} and the H_{2}O, H_{2}O^{+} pairs. Figure 13 shows the behavior as a function of fractional occupation number compared to straight lines.

As we were gradually moving to transition metal systems, we report some results for them, too, but primarily, we redid all the other tests for activation barriers, vibrational frequencies, and excitation energies, providing uniform results for QTP(02). They are usually close to QTP(00) and QTP(01) but often a little better. See Table VI.

. | . | CAM- B3LYP . | CAM- QTP-00 . | CAM- QTP-01 . | CAM- QTP-02 . | EOM-CCSD . |
---|---|---|---|---|---|---|

Vertical excitation energies | All | 0.33 | 0.60 | 0.29 | 0.36 | 0.27 |

Valence excitations | 0.31 | 0.71 | 0.46 | 0.49 | 0.47 | |

Rydberg excitations | 0.35 | 0.52 | 0.17 | 0.26 | 0.11 | |

Charge-transfer excitations | Aromatic- | 0.70 | 0.14 | 0.08 | 0.21 | |

tetracyanoethylene | ||||||

complexes | ||||||

Electron affinities | Atoms | 1.02 | 0.19 | 0.27 | 0.16 | |

Molecules | 0.82 | 0.19 | 0.18 | 0.10 | ||

Transition metals | 0.89 | 0.15 | 0.15 | 0.09 | ||

Barrier heights | Hydrogen transfer | 2.8 | 2.3 | 2.7 | 2.2 | |

reactions | ||||||

Nonhydrogen | 2.38 | 3.02 | 2.22 | 2.15 | ||

transfer reactions | ||||||

Geometry | Bond length | 0.011 | 0.022 | 0.015 | 0.017 | |

Angle | 2.0 | 2.4 | 2.2 | 2.3 | ||

Harmonic frequency | 48 | 145 | 75 | 91 | ||

Dipole moments | 0.08 | 0.14 | 0.11 | 0.13 |

. | . | CAM- B3LYP . | CAM- QTP-00 . | CAM- QTP-01 . | CAM- QTP-02 . | EOM-CCSD . |
---|---|---|---|---|---|---|

Vertical excitation energies | All | 0.33 | 0.60 | 0.29 | 0.36 | 0.27 |

Valence excitations | 0.31 | 0.71 | 0.46 | 0.49 | 0.47 | |

Rydberg excitations | 0.35 | 0.52 | 0.17 | 0.26 | 0.11 | |

Charge-transfer excitations | Aromatic- | 0.70 | 0.14 | 0.08 | 0.21 | |

tetracyanoethylene | ||||||

complexes | ||||||

Electron affinities | Atoms | 1.02 | 0.19 | 0.27 | 0.16 | |

Molecules | 0.82 | 0.19 | 0.18 | 0.10 | ||

Transition metals | 0.89 | 0.15 | 0.15 | 0.09 | ||

Barrier heights | Hydrogen transfer | 2.8 | 2.3 | 2.7 | 2.2 | |

reactions | ||||||

Nonhydrogen | 2.38 | 3.02 | 2.22 | 2.15 | ||

transfer reactions | ||||||

Geometry | Bond length | 0.011 | 0.022 | 0.015 | 0.017 | |

Angle | 2.0 | 2.4 | 2.2 | 2.3 | ||

Harmonic frequency | 48 | 145 | 75 | 91 | ||

Dipole moments | 0.08 | 0.14 | 0.11 | 0.13 |

But perhaps the most surprising result is the excellent description of EAs obtained solely from the lumo eigenvalue. For all the main group atoms through Kr and 21 small closed shell molecules, the MAD was 0.16 eV for the atoms and less than 0.1 eV for the molecules. Even for molecules containing transition metals, the MAD is 0.09 eV. This is to be compared to the original CAM-B3LYP whose MAD is 1 eV for the atoms, 0.8 eV for the molecules, and 0.9 eV for those with transition metal atoms.

Note that even the computed sign of an EA is a very sensitive property. According to experiment when available and high-level predictions from the energy difference between the anion and neutral using ΔCCSDT-3,^{90} otherwise, QTP(02) almost always gets the correct sign. This is usually the case for QTP(00) and (01), too, but their MAD is slightly higher, as shown in Table VI.

It should be apparent that if the IPs and EAs are correct, then the charge transfer (CT) energy (E(A^{+}) + E(B^{−}) + 1/R_{AB}) in the limit of complete separation has to be ϵ_{A}(homo) + ϵ_{B}(lumo). Hence, the QTP functionals also solve the CT problem that has been a long-term problem for KS-DFT.^{74} See Table VII.

. | CAM-B3LYP . | CAM-QTP(00) . | CAM-QTP(01) . | CAM-QTP(02) . | Expt. . |
---|---|---|---|---|---|

Benzene | 2.94 | 3.79 | 3.69 | 3.85 | 3.59 |

Toluene | 2.63 | 3.46 | 3.35 | 3.50 | 3.36 |

o-xylene | 2.38 | 3.21 | 3.10 | 3.25 | 3.15 |

Naphthalene | 1.85 | 2.71 | 2.65 | 2.80 | 2.60 |

2.62 | 3.48 | 3.42 | 3.57 | 3.23 | |

MAD | 0.70 | 0.14 | 0.08 | 0.21 | |

MAX | 0.77 | 0.25 | 0.19 | 0.34 |

. | CAM-B3LYP . | CAM-QTP(00) . | CAM-QTP(01) . | CAM-QTP(02) . | Expt. . |
---|---|---|---|---|---|

Benzene | 2.94 | 3.79 | 3.69 | 3.85 | 3.59 |

Toluene | 2.63 | 3.46 | 3.35 | 3.50 | 3.36 |

o-xylene | 2.38 | 3.21 | 3.10 | 3.25 | 3.15 |

Naphthalene | 1.85 | 2.71 | 2.65 | 2.80 | 2.60 |

2.62 | 3.48 | 3.42 | 3.57 | 3.23 | |

MAD | 0.70 | 0.14 | 0.08 | 0.21 | |

MAX | 0.77 | 0.25 | 0.19 | 0.34 |

One more property that is essential to solid state theory is the bandgap. The bandgap is defined to be the E(IP) − E(EA) for an infinite system such as a crystalline solid. It is meaningless to talk about any kind of Δ or functional, energy difference such as E(N) − E(N − 1) for an IP and E(N + 1) − E(N) for an EA for an infinite system. But once the theory is such that these differences are computed directly as eigenvalues, then the bandgap becomes the difference between the homo and lumo energy bands, gap = ϵ_{homo} − ϵ_{lumo}.

Using the eigenvalues produced from IP-EOM-CC and EA-EOM-CC, we can study the bandgap as a function of the size of molecular oligomers with conventional square-integrable functions. Ranasinghe^{91} did this using ACES II and its massively parallel ACES 3 version for polyacetylene and polyacenes. See Figs. 14 and 15.

The *ab initio* CC calculations stop at 9 and 6 units, respectively, while the rest of the results are obtained from various KS-DFT functionals. The important fact is that the QTP(00) and QTP(01) results exactly follow the *ab initio* values, permitting extrapolating beyond that. Of course, one could even take the KS-DFT results to the infinite polymer limit. The other functionals considered do not provide the right CC answer.

### The Devil’s triangle

As shown in Fig. 16, the “devil’s triangle” of (1) self-interaction, (2) correct one-particle spectrum, and (3) the integer discontinuity and associated curvature errors threaten the ability of KS-DFT to give the right answers for problems. Indeed, they even potentially destroy the foundation for KS-DFT calculations. All three issues are different manifestations of the same underlying problem. When Jónsson corrects for self-interaction in KS-DFT by generalizing the Perdew-Zunger^{92} method to use complex orbitals,^{93} or Der-you Kao *et al.*^{94} and Pederson^{95,96} accomplish much the same thing by exploiting orbital localization, the first quantity used to demonstrate success is the resultant KS orbital eigenvalues for the homos. If SI is mitigated, the eigenvalues are vastly more accurate. We depend on the converse. If the eigenvalues provide correct IPs, SI is, necessarily, ameliorated.

Furthermore, none of our work is limited to the homo as in our opinion, all occupied state eigenvalues should be similarly well-described as proven by our TDDFT argument. When Gidopoulos and Lathiotakis^{97,98} correct for self-interaction by insisting that the resultant density obtained from a KS-DFT calculation integrates to N − 1 electrons instead of the usual situation where the density would correspond to N, Gidopoulos finds that he gets the eigenvalues for all the valence occupied orbitals to be good approximations to their respective principal IPs. This even includes the core once adding 50% NLX,^{99} as we did in OEP2sc. The above three routes are “active” as they change elements of the theory to eliminate self-interaction. In the case of the QTP functionals, the role is “passive” since we expect self-interaction to be mitigated by the simple expedient that once the method provides accurate IPs, particularly for all occupied orbitals, then SIE is attenuated.

If one needed more proof, then the fact that one also gets excellent results for delicate anions further attests to the fact that there is little SIE. The existence of anions in KS-DFT when providing the calculation the freedom for the extra electron to float away (see also Ref. 88), by potentially having some plane waves in the basis set or making the basis more and more diffuse, is a perfect example of a failure due to the SIE error. This should be the first condition checked to assess SIE.

We also see that QTP(02) directly uses the condition that a molecule and its cation, or a molecule and its anion, should be simply related by straight lines, when plots of energies as a function of fractional occupation numbers are obtained in KS-DFT. This is an exact, i.e., nonempirical, condition imposed on KS-DFT to fix two parameters in CAM-BLYP to define CAM-QTP(02). Yet even without this explicit requirement, we still get much straighter lines from CAM-QTP(00) and CAM-QTP(01) than other widely used functionals reflecting that the QTP family fulfill our “consistency conditions.” If the IPs and EAs are right, this helps a great deal in avoiding any error due to the integer discontinuity, not to mention charge-transfer.

The bottom line is all three problems of the devil’s triangle will be reduced when KS-DFT approximations adhere to the IP eigenvalue theorem. The IP-eigenvalue issue should be under control before considering other empirical modifications. Much the same benefit would appear to pertain to solid state applications for crystals, subject to due regard for the differences between finite and infinite systems and the limitations caused by 100% NLX at large r in the range-separated models discussed here.

The second level consequences of the devil’s triangle pertain to excitation energies to Rydberg states, for charge transfer situations, and for transition states and activation barriers. As our results show, all of these are also relatively well described by the QTP functionals. If one also ensures size-extensivity and intensivity, which would be guaranteed by the exponential wavefunction of CC theory and EOM-CC and is implicit in the COT equations, most of the exact conditions required in electronic structure theory would be incorporated. Putting these elements together provides a “quantitative molecular orbital framework” to underlie the conceptual benefits of qualitative MO theory but now augmented with a complementary numerical, quantitative approach to the theory that can be used hand-in-hand with its qualitative arguments. Together, one would hope this would explain and, now, quantify a lot of chemistry.

My image for an effective one-particle theory to provide CC accuracy for most properties, even energy bands and excitation spectra for solids, would seem to be achievable. Future work in this direction should include introducing new concepts for functionals that contain genuine correlation and exchange effects. This will benefit from further study of the *ab initio* correlation potential defined by COT, including adding another step to the procedure that besides getting the IPs and EAs for a given set of orbitals, as shown in Eq. (7), will also permit mixing the occupied and virtual orbitals to converge onto an improved set of MOs for a problem. Future work should also address TDDFT beyond the adiabatic approximation to provide better approaches to obtaining EOM-CC quality results for all sectors of Fock space, including the principal IPs. Our adiabatic argument is the first step, but achieving a “predictive” level for principal IPs should make everything better. Furthermore, it is artificial and limiting to act as if KS-DFT is only a ground state theory, when the problems in chemistry demand a knowledge of many states. Further synergy between KS-DFT and TDDFT will help to define much improved methodology, as they do for the QTP family.

One final point: All the conditions we consider in trying to make KS-DFT provide CC quality answers are exact ones, starting with the COT orbital eigenvalues. Augmenting them with further orbital optimization should provide a framework for a purely *ab initio* MO theory that reduces the correlation problem to that for a correlated orbital potential (COP). Second, there should never be a self-interaction error. The third condition is to ensure there is no “curvature error” when adding or removing an electron. Thus, all three components of the devil’s triangle are addressed to some degree by the QTP functionals and their numerical successes demonstrate the critical role these conditions should have in future attempts to define improved KS-DFT methods. To this point, what we have reported is only the first, simple step, but it seems to be a relevant one to build upon to attain a truly seamless connection between wavefunction theory and KS-DFT to provide new, reliable, and easily applied electronic structure methods to the wealth of problems in chemistry, condensed matter physics, and materials science.

## ACKNOWLEDGMENTS

This work has been supported by the U. S. Air Force Office of Scientific Research over many years (Grant Nos. FA9550-11-1-0065, FA9550-14-1-0281, FA9550-04-0119, and FA9550-04-01-01) under the kind oversight of Dr. Mike Berman. I thank all the exceptional co-workers that contributed their talents to the work reported. I particularly appreciate Professor Nikitas Gidopoulos, Durham University, UK, for daily discussions and helpful suggestions for this manuscript. This perspective was written while enjoying the hospitality of the Institute for Advanced Study, Durham University, UK, during a three month stay. Dr. Young Choon Park was instrumental in updating many of the results shown, assisted by Ms. Moneesha Ravi.

## REFERENCES

For a multideterminantal problem such as two electrons in two orbitals, the DIP and DEA-EOM-CC can be used since the target state relative to a core with ±2 electrons has the four determinants in it, {$\Phi 0,\Phi hl,\Phi h\u21c0,l\u21c0\Phi hh\u21c0ll\u21c0$}, for a homo, h, and a lumo, l, e.g., Φ_{0} has the homo doubly occupied. The weights of these determinants will be determined by the EOM-CC-eigenvector step.^{2} This means that they will have whatever weight required, avoiding the limitation of a single reference method where the weights of quasidegenerate determinants in the Q space could be constrained from growing to appropriate values by truncation of the excitation manifold in the CCSDT … method.

The frequency dependence of the Dyson self-energy is removed by replacing the operator by a matrix representation composed of the usual single, double, … excitations used in IP/EA-EOM-CC while first isolating its principal IPs and EAs. These two steps provide the correlated one-particle operator expression for the MOs. In addition to the principal IPs obtained as the eigenvalues for each of the orbitals, the remaining “shake-up” IPs composed of ionizations followed by electronic excitation would also be described by the additional IP-EOM-CC solutions obtained from $H\xaf$ but will typically require triple or higher excitations for an accurate description, so those effects are conveniently separated from the COT eigenvalue problem.

There is no need as in the KLI approximation to allow these orbitals to be arbitrary since the sum over i, j cancels identically. The so-called CEDA—O. V. Gritsenko and E. J. Baerends, Phys. Rev. A **64**, 042506 (2001), M. Grüning *et al.*, J. Chem. Phys. **116**, 6435 (2002) simply applies an Unsöld (constant energy denominator approximation) to this correct form^{32} instead of the arbitrary orbital form used in KLI.^{31}

“We have shown in a variety of ways that the correlation potential corresponding to MP2 diverges for finite systems. Consequently, a direct self-consistent application of MP2 is not possible.”—A. F. Bonetti *et al.*, Phys. Rev. Lett. **86**, 2241 (2001).

In addition to Ref. 103, a few other examples include, P. Mori-Sanchez, Q. Wu, and W. Yang, J. Chem. Phys. **123**, 062204 (2005). Who say “This failure basically excludes PT based orbital functionals for practical applications in DFT—Jiang and Engel show that exact degeneracy can be encountered in doing their OEP calculations—H. Jiang and E. Engel, J. Chem. Phys. **123**, 224102 (2005). D. Rohr, O. Gritsenko, and E. J. Baerends, Chem. Phys. Lett. **432**, 336 (2005) construct an example where the homo and lumo are degenerate, to argue that OEP exchange or second-order correlation cannot work—Mori-Sanchez *et al.* further erroneously blame the failure of using doubles only MBPT2 as a correlation potential on the neglect of singles. All of these problems are manifestations of using the wrong PT as our results show.^{49}

When the orbitals must be of a certain form like localized, then one cannot use an orbital rotation to account for the off-diagonal f_{ij} and f_{ab} terms, so the iterative equations would need to be solved instead. Similarly, if one wants to keep the KS orbitals as they are, this can be accomplished also via the iterative solution.

Quoting R. G. Parr and Y. Weitao, in *Density-Functional Theory of Atoms and Molecules* (Oxford University Press, 1994), “Given the auxiliary nature of the KS orbitals-just N orbitals the sum of squares of which add up to the true total electron density-one should expect no simple physical meaning for the Kohn-Sham orbital energies. There is none.” And from W. Kohn, A. D. Becke, and R. G. Parr, J. Phys. Chem. **100**, 12974, (1996), “The individual eigenfunctions and eigenvalues, ε_{j} and φ_{j}, of the KS equations have no strict physical significance….”

(1) Basic properties, which include main group atomization energies (MGAE109), electron affinities (EA13), ionization potentials (IP13^{*}), and proton affinities. (2) Reaction barrier heights, which include hydrogen and nonhydrogen atom transfer barrier heights (HTBH38 and NHTBH38, respectively). (3) Reaction energies, which include alkyl bond dissociation energies (ABDE12), hydrocarbon reactions (HC7), π-system thermochemistry (pTC13), noncovalent complexation energies (NCCE31), multireference bond energies (MRBE5^{*}), isomerizations of large molecules (ISOL6), and “difficult cases” (DC9).

Note that for the self-interaction free OEPx method, also done in an effectively complete basis set, the curvature is not the convex curve usually obtained as shown but is essentially flat or even slightly concave, demonstrating the effect of removing SI—T. W. Hollins, “Local exchange potentials in density functional theory local exchange potentials in density functional theory,” Ph.D. thesis, University of Durham, 2014, and any ambiguity due to finite bases. This behavior offers a hint about correlation that harkens back to HF-DFT using BLYP. A SI free HF calculation still shows a concave curvature, only removed by MBPT2 and CCD in this example. Thus, it appears that a correct correlation potential, at least relative to HF, needs to include the “correlation” effect achieved by the local OEPx.