Local reduced-density-matrix-functional theory: Incorporating static correlation effects in Kohn-Sham equations

We propose a novel scheme to bring reduced density matrix functional theory (RDMFT) into the realm of density functional theory (DFT) that preserves the accurate density functional description at equilibrium, while incorporating accurately static and left-right correlation effects in molecules and keeping the good computational performance of DFT-based schemes. The key ingredient is to relax the requirement that the local potential is the functional derivative of the energy with respect to the density. Instead, we propose to restrict the search for the approximate natural orbitals within a domain where these orbitals are eigenfunctions of a single-particle hamiltonian with a local effective potential. In this way, fractional natural occupation numbers are accommodated into Kohn-Sham equations allowing for the description of molecular dissociation without breaking spin symmetry. Additionally, our scheme provides a natural way to connect an energy eigenvalue spectrum to the approximate natural orbitals and this spectrum is found to represent accurately the ionization potentials of atoms and small molecules.


I. INTRODUCTION
Computational simulations have opened a new avenue for the exploration and prediction of "à la carte" molecular complexes and materials, i.e., with tailored properties and functionality, due to the development of powerful algorithms and an increase in computational power. Generally, an increase in complexity or system size goes hand-in-hand with a loss of accuracy or an increase in numerical cost.
Highly accurate wave function methods, which have recently become available also for solids [1], are limited in the complexity of systems they can handle. This hinders the application of these methods to complex molecular structures, e.g. nanostructures and biomolecules. However, now, more than ever, there is a need for methods that are able to handle large-scale systems with high precision in order to assess the challenges in material, bio-, and nanosciences.
Density functional theory (DFT) [2,3], on the other hand, comes with a low computational cost, which allows its application to rather large systems. The lack of a systematic way to improve functionals is hampering the progress towards a more accurate and efficient DFT first-principles scheme. Currently, several properties are still difficult to predict within standard DFT, with the band gap and molecular dissociation being prominent examples. The latter is described correctly with a standard local density approximation (LDA) or generalized gradient approximation functional only if the spin symmetry is artificially broken. The recently introduced strictly correlated functional [4] rectifies this problem, however, at the expense of a wrong description of the equilibrium properties. Although the exact exchange-correlation (xc) functional of DFT is universal, one could of course use different approximations depending on the physical situation or quantity of interest. To this end, large databases collecting information from different functionals are compiled [5]. The lack of a proper description of static correlation is at the heart of many of the failures of present xc approximations of DFT for describing strongly correlated systems and molecular bond breaking and formation. This deficiency can be traced back to the fact that within Kohn-Sham (KS) DFT the density is reproduced via a single Slater determinant. When the true wave function has multireference character the KS kinetic energy is a poor approximation to the true kinetic energy. The difference has to be compensated by the xc energy which for many approximations is not done very successfully, leading to a wrong total energy in the case of molecular dissociation.
In electronic structure theory, a desirable feature of independent particle models, like the KS scheme in DFT, is the direct prediction of single-electron properties, like ionization potentials (IPs), from the eigenvalues of the single-particle Hamiltonians. Although the question about the physical content of the KS orbital energies raised a scientific debate [6], theoretical justification for this result was given by Baerends and coworkers [7,8] and by Bartlett and co-workers [9,10]. Unfortunately, orbital energies from approximate func-tionals tend to underestimate substantially the IPs of molecular systems [11,12].
Reduced-density-matrix-functional theory (RDMFT) [13] is an alternative to DFT to approximate the manyelectron problem. It is based on approximating the total energy of an electronic system in terms of the one-body reduced density matrix (1RDM). A main difference from DFT is the introduction of fractional occupation numbers which allows the exact treatment of the kinetic energy and potentially leads to improved accuracy whenever the ground-state many-body wave function is far from a single Slater determinant. Most approximations in RDMFT are explicit functionals of -and are minimized in terms of -the natural orbitals (NOs), φ j (r), and their occupation numbers, n j . So far, various different approximations for the total energy functional have become available [14][15][16][17][18][19][20][21][22][23][24][25] which have proven to describe correctly such diverse properties as molecular dissociation [16][17][18][19][20] or band gaps [25][26][27][28]. However, applications have been restricted to small molecules due to the computational cost to determine the orbitals. Although the optimization of the occupation numbers is a relatively inexpensive task, orbital minimization is complicated: it does not reduce to an iterative eigenvalue problem (as in DFT) and requires numerically expensive procedures. Significant effort has been devoted to devising effective Hamiltonians [29][30][31] to improve the efficiency, but with limited success so far.
If we were able to incorporate all the advantages of RDMFT functionals into DFT while keeping the cost of standard DFT functionals we would make a big step forward in constructing an efficient and accurate scheme able to cope with the challenge of describing structural and dynamical properties of many-electron systems including bond breaking and bond formation. In this paper, we propose such a framework that combines the best of both DFT and RDMFT. One can regard this approach as either an extension of DFT, where fractional occupations for the orbitals are introduced using an approximation for the xc energy functional borrowed from RDMFT, or, alternatively, as a constrained RDMFT calculation. In either case, we incorporate the proper nonidempotent nature of the density matrix in the calculation of the kinetic energy that is fundamental to the success of RDMFT approximations.
The central idea in our proposed framework, which we call local RDMFT, is to restrict the minimization with respect to the NOs to a domain where these orbitals are eigenfunctions of a single-particle Hamiltonian with a local potential. The best possible Hamiltonian is the one whose eigenorbitals minimize the total energy. The resulting equations are similar to the optimized effective potential (OEP) equations [32][33][34][35]. The OEP improves the accuracy in DFT and significant expertise has been developed in the past two decades [34] for its efficient implementation. Hence, our method can be implemented directly in existing DFT codes with only small modifications to address fractional occupation numbers. Fractional occupations, as in standard RDMFT, are provided by the minimization of the energy functional under the appropriate conditions. Our approach has some similarity with the idea explored by Grüning et al. [36], where the common energy denominator approximation is used together with the Müller functional. In that approach, occupancies are obtained in an empirical way and not through optimization.
The local RDMFT framework provides an energy eigenvalue spectrum connected to the NOs and as we show, single electron properties, like the IPs of small molecules, are well reproduced by the energies of the highest occupied molecular orbitals (HOMOs).
In Sec. II, we describe in detail the formalism of local RDMFT. Then in Sec. III, we show that the restriction to a local effective potential has little effect on the dissociation of dimers; hence, the accuracy of RDMFT is retained by the proposed method. In Sec. III we also illustrate the quality of the energy spectrum provided by the effective Hamiltonian by comparing the obtained IPs with experiment. In the Appendix we show that pure density xc functionals are not adequate for the present scheme since they cannot lead to fractional occupations in a minimization procedure and we need to employ functionals of the 1RDM as we do in the present work.

II. LOCAL RDMFT
Clearly, the integration of DFT with RDMFT is desirable as it could combine the best of both worlds. With fractional occupation numbers, static correlation would become accessible, while use of a common local potential to yield the NOs would improve dramatically the efficiency of the method. However, this target is not straightforward. A natural way to incorporate a local potential for the orbitals, is to consider the approximate RDMFT xc energy to be a density functional, so that its functional derivative would give the local potential. Unfortunately, such an approach leads in general to an idempotent density matrix (see Appendix) and we are back to square one.
Hence, we abandon the requirement that the xc energy must be a functional of the density and that the potential must arise as a functional derivative with respect to the density. Instead, we consider the total energy as a functional of the one-body-reduced-density matrix (1RDM), γ(r, r ), where v(r) is the external local potential and T [γ(r, r )] is the interacting kinetic energy which is a functional of the 1RDM. The electron-electron interaction energy can be cast into the last two terms in Eq. (1) where can be considered a functional of the occupation numbers n j and the natural orbitals φ j , i.e., the eigenvalues and the eigenvectors of γ: The central idea in our proposed local RDMFT scheme is to restrict the search for the optimal φ j 's within a domain where they are also eigenfunctions of a singleparticle Hamiltonian with a local potential, v rep (r), where v rep is the effective repulsive potential acting on any electron in the system, caused by the effective repulsion of the remaining N − 1 electrons (atomic units are used throughout the paper). Fractional occupations, n j , as in standard RDMFT, are provided by the minimization of the energy functional of Eq. (1) under the appropriate N -representability conditions. The natural orbitals in the exact theory, as well as the minimizing orbitals in RDMFT approximations, are typically satisfying a Schrödinger equation with a nonlocal effective potential. The local potential constraint leads to approximate natural orbitals (ANOs), φ j , which cannot become equal to the true natural orbitals.
By enforcing Eq. (3), the total energy becomes a functional of the local effective potential and of the occupation numbers, In the same way as in the OEP method [32][33][34][35], the optimal local potential is obtained by solving the integral equation where χ(r, r ), a generalized density-density response function, and b(r) are given by Hxc defined by δE Hxc δφ * j (r) Here E Hxc is the approximation for the electron-electron interaction energy, i.e., the last two terms in Eq. (1).
In the summations of Eqs. (5) and (6), we have excluded terms over pairs of orbitals differing in occupation by less than a cutoff value ∆n c . This choice excludes pairs of weakly occupied orbitals whose energies j are not accurate for finite localized basis sets; e.g., occasionally they violate the Aufbau principle and the negative definiteness of χ. By excluding these terms, we also partly alleviate a common inaccuracy of most RDMFT functionals to show a spurious excess of total occupation in weakly occupied orbitals [38]. In that way, we have examined the dependence of our solution on ∆n c . Typically, for very small values of ∆n c we run into convergence issues and IPs vary substantially as a function of ∆n c . As ∆n c increases to a typical value of 0.1 convergence improves dramatically and IPs remain unchanged as a function of ∆n c for a broad range of values. Even for values of ∆n c larger than the HOMO-lowest occupied molecular orbital (LUMO) occupation difference, we do not find any noticeable change in the results. Further increase deteriorates the accuracy since fewer terms are included in the summations. A choice for ∆n c ∼ 0.1−0.3 usually leads to a converged solution. Accidental exclusions of strongly-weakly occupied pairs have a negligible effect on the results.
To ensure a physical asymptotic decay of the effective repulsive potential, we do not solve Eq. (4) directly. Instead, we follow the methodology in Ref. [39] and express v rep (r) as the electrostatic potential of an effective repulsive density, ρ rep (r): The requirement [39] that the effective repulsive density corresponds to a (fictitious) system of N − 1 electrons repelling the electron at r yields the following two constraints: The first condition is necessary for the asymptotic decay of the effective repulsive potential as (N − 1)/r, which is a property of the exact Hxc potential [40]. The two conditions together become sufficient (although probably not necessary anymore) to guarantee the correct asymptotic behavior and a well-posed mathematical problem. Minimization of the total energy leads to an integral equation for the effective repulsive density: The two constraints can be incorporated with a Lagrange multiplier (9) and a penalty term (10) that introduces an energy cost for every point where ρ rep becomes negative.
We expand the ANOs in a basis set (orbital basis) and the effective repulsive density (rather than the potential) in another (auxiliary) basis and we obtain, similarly to Ref. [39], a linear system of equations for the expansion coefficients of the repulsive density. The inversion of the linear equations is complicated by the fact that the matrix of the response functionχ becomes singular when the auxiliary basis is large compared to the orbital basis [41][42][43][44]. As a result, the effective repulsive density becomes indeterminate in the null space of the (finiteorbital-basis) response function. This indeterminacy is substantially suppressed by the two constraints, Eqs. (9) and (10), that reduce drastically the form of the admissible effective repulsive densities and effective repulsive potentials. We have found that with a singular value decomposition, we obtain systematically smooth and physical densities and potentials (see Ref. 39).
A consequence of the local-potential approximation is that the asymptotic decay of the ANOs depends on the energy eigenvalue j and hence, differs from the (necessarily uniform) asymptotic decay of exact NOs with fractional occupancy [45,46]. As a result, the asymptotic exponential decay of the density is related to the highest energy eigenvalue with nonzero occupation and not to the IP. The effect in the total energy from the different asymptotics -compared with the behavior of the exact NOs -is negligible since the energy contribution of the asymptotic region is insignificant.
Similarly to standard RDMFT, nominal scaling of local RDMFT is N 4 × M , where M is the total number of different generalized Fock matrices [defined in Eq. (7)] that need to be evaluated and then transformed to obtain the matrix elements in Eq. (6). For most approximations equal occupations result in equal generalized Fock matrices. Hence, M is typically equal to the number of orbitals with fractional occupations plus an additional one common for all the fully occupied orbitals. Since M can be fixed, the nominal scaling of local RDMFT reduces to that of Hartree-Fock (HF) or hybrid DFT calculations. The number of different F (j) Hxc (r, r ) that need to be evaluated, and the additional cycle of convergence of the occupations apart from the ANOs add to the overall computational cost of local-RDMFT as compared to standard DFT and Hartree-Fock calculations.
The present scheme would benefit from the recent developments in reducing the scaling of wave-functionbased schemes [47]. Significant numerical cost saving is also expected by improving the enforcement of the positivity condition of Eq. (10). However, the search for the optimal repulsive density through an iterative eigenvalue equation is still much more efficient than full minimization for the whole set of natural orbitals. This computational efficiency represents a big advantage of local RDMFT over standard RDMFT allowing for the first time relatively large systems to be in the capability of RDMFT. The application of local RDMFT to larger systems gives very promising results [48]. Finally, one advantage of this method is that it can be easily implemented in standard electronic structure codes due to its similarity with the OEP method. Finally, with regard to the efficiency of our scheme, our target is not to replace well-established methods in routine calculations where these methods are known to be accurate but to complement them and improve over their results where they deviate from experiment. Such cases are for instance bond breaking and highly correlated systems where static correlations are important. At the same time, in our scheme, the orbital energies from the effective Hamiltonian offer improved spectral properties.

III. APPLICATIONS
In Fig. 1, we show the effective local potential for a Ne atom employing two RDMFT approximations, the Müller [14,16] and BBC3 [17] functionals. As we see, the optimal potentials are similar to the exact-exchange OEP (x-OEP), especially for BBC3 while local Müller is closer to the exact KS [37] potential. The comparison with the local potentials of different theories and the exact KS scheme are useful since a reasonable approximate potential resembling the exact KS potential will hopefully lead to a reasonable single-particle spectrum.
An important advantage of many RDMFT functionals is the qualitatively correct dissociation of small molecules, in contrast to available xc density functionals. One example is the H 2 molecule and, as we show in Fig. 2(top), this property holds for local RDMFT as well: for all three functionals, the Müller, BBC3, and the power functionals, local RDMFT reproduces the correct dissociation. The description at the equilibrium distance also agrees very well with the configuration interaction (CI) results [50,51] both for the position of the minimum and the curvature as can be seen in Table I. In the inset of Fig. 2(top), we show the difference between the energy eigenvalues corresponding to the HOMO and the LUMO, as they are defined by the number of electrons and the assumption of a single-electron picture. We find the same behavior as the one found for various DFT functionals in Ref. [52], i.e., the energy difference is approaching zero at small distances between the two hydrogen atoms. In Fig. 2(bottom), the dissociation curves for the triple bond of the N 2 molecule are shown for BBC3 and power functionals in comparison to HF and a DFT calculation [50,51] using the B3LYP functional [53]. Unfortunately, the Müller functional fails to describe the dissociation of N 2 . The binding energy as well as the vibrational frequency obtained with BBC3 are closer to the experimental values than those given by the power functional. Interestingly local RDMFT results are qualitatively much better than many density functional approximations including B3LYP, as seen in Fig. 2(bottom). It is worth mentioning that the correct dissociation of N 2 is not reproduced even at the level of accurate quantum chemistry methods like MP2 and single-reference coupled cluster.
We now focus on the single-electron spectrum of the local potential Hamiltonian. For this purpose, we obtained IPs as the corresponding eigenvalues of the local potential Hamiltonian for a test set of atomic and molecular systems and basis sets. This set comprises small atomic and molecular systems using the cc-pVTZ and uncontracted cc-pVTZ basis sets for the orbital and the auxiliary basis respectively, and we obtained IPs up to the third one. Our numerical results for several RDMFT functionals are shown in Table II where on the bottom we also show the average, absolute, percentage error of IPs. The same errors for the IPs obtained as the energy eigenvalues of Hartree-Fock (Koopmans' theorem) are also included for comparison. We find a remarkable agreement between the energy eigenvalues and the experimental IPs for the functionals we tested. All errors are below or around 5%, with the Marques-Lathiotakis (ML) functional [21] going as low as ∼3%. The agreement with experiment is even better for the first IPs with the Müller functional being the most accurate with an error of only ∼2%. Overall the agreement with experimental values is very good and substantially better than the HF Koopmans' theorem.
To estimate the effect of the local potential approximation we include in Table II the ratio E c of the correlation energies (defined as the energy differences from Hartree-Fock) with local RDMFT over those provided by the full RDMFT minimization. The average E c over the whole set of systems, included in the last row of Table II, is in the range 0.6-0.8% for all functionals considered. However, as in DFT calculations, the comparison of the obtained correlation energy to the exact one is not the only decisive factor to assess the accuracy of an approximation in reproducing many properties. For the Müller functional, the local RDMFT recovers on average 71% of the full minimization correlation energy. The Müller functional generally overestimates the correlation energy substantially in the full minimization and the constraint of the local potential offers an improvement. For more accurate functionals, however, this is not always the case.

IV. CONCLUSION
In conclusion, we presented a method on how to incorporate static correlation into KS-like equations by employing xc functionals from RDMFT. We have shown that when the xc energy is a density functional then the total energy minimization leads to an idempotent solution (Appendix). Consequently, we relaxed the requirement that the potential must be the functional derivative of the energy with respect to the density and decided to minimize the total energy with respect to the occupation numbers and the ANOs generated by a local effective  TABLE II. Ionization potentials (in eV) for atoms and small molecules compared to experiment, and the ratio, Ec, of correlation energy captured by the local approximation over that of a full RDMFT minimization for several functionals. ∆ all , ∆1st, are the average absolute errors for all IPs and for the first IP, respectively, defined as ∆ = 100 and Ec is the average of Ec over the whole set of systems. Vertical experimental IPs in the last column are obtained from the NIST Chemistry WebBook [54] and references in Ref. [55].
potential. In this way, we manage to describe the dissociation of diatomic molecules accurately. In addition, our approach allows us to connect a single-particle energy spectrum to the ANOs. This spectrum is in good agreement with experimental IPs and photoelectron spectra for molecules.
The proposed method provides a powerful tool which opens a new avenue: physically motivated approximations in density-matrix based schemes, able to cope with strongly correlated systems [25] and static correlation, can now be brought to the realm of DFT. The resulting KS-like approach is able to simultaneously describe ground-state properties, bond breaking, and single-electron spectral properties.
The scaling of our method can be easily reduced to that of hybrid DFT methods using standard techniques. The improved computational efficiency compared to full RDMFT minimization allows for the application to large systems opening the road for an improved description of electronic correlations in technologically important molecular systems.  In KS DFT the (exact or approximate) ground-state density ρ v and the ground-state energy E v of an interacting N -electron system in an external potential v(r) are obtained by minimizing, over N -electron densities ρ, the where T s [ρ(r)] is the noninteracting kinetic energy. For our analysis, it does not matter if density-functional of the the exchange and correlation energy in Eq. (A.1) is exact or approximate, e.g., given by LDA, E LDA xc [ρ]. Strictly, only the infimum of the KS total energy functional is defined. Nevertheless, we assume routinely that a minimizing density exists within the space of Nelectron, noninteracting v-representable densities, inside where we search for the minimum. This assumption (i.e., the existence of the minimum) amounts to restricting the domain of interacting v-representable densities under study to include only those that are also noninteracting v-representable.
On the other hand in RDMFT, the ground-state 1RDM, γ(r, r ), of the same interacting system, and its ground-state energy E v , are obtained by minimizing, over all N -representable 1RDMs γ(r, r ), the functional of Eq. (1). Again, it makes no difference for our analysis if E xc [γ(r, r )] is exact or an approximation.
It is desirable to bring the two theories together (KS DFT and RDMFT) and to obtain (approximately) the interacting, nonidempotent 1RDM γ v and the total energy E v by solving appropriate single-particle equations with a local effective potential v(r). It appears this aim can be achieved easily: (a) We can replace the noninteracting kinetic energy in Eq. (A.1) by the interacting kinetic energy, written in terms of the 1RDM. At the same time we keep the dependence of the xc energy on the density ρ(r), after correcting if necessary the xc energy functional for the kinetic part of the correlation energy. (b) Equivalently, we can regard the xc energy in Eq. (1) as a functional of the density γ(r, r), rather than the full 1RDM. Both ways amount to attempting to obtain approximately the interacting ground-state 1RDM γ v and the total energy E v from the minimization of a total energy expression having the following form: Compared with the KS DFT scheme, Eq. (A.1), where the 1RDM is constrained to be idempotent, the search for the 1RDM in this scheme, Eq. (A.2), has greater variational freedom. Below, we show the following corollary from Janak's theorem [56].
Excluding degeneracy of the noninteracting ground state, the minimizing N -representable γ s (r, r ) of This result is unexpected, since despite the greater variational freedom compared to Eq. (A.1), the solution of Eq. (A.2) is still idempotent. Also, from the corollary, it appears that employing a local potential to generate approximate NOs for the minimization of E RDM v [γ(r, r )] in Eq. (1), would always lead to an idempotent 1RDM (unless the ground state is degenerate). Although this is not true, it indicates the difficulty to combine a nonidempotent 1RDM with a local potential. The way out is to give up the requirement that the local potential must be the functional derivative of the xc energy with respect to the density. Proof: We must minimize E v [γ] under the usual Nrepresentability constraints for the 1RDM, i.e., that the eigenfunctions φ i (r) of γ(r, r ) are normalized and that their occupation numbers satisfy 0 ≤ n i ≤ 1 and i n i = N . Hence, we must minimize: where λ i , µ are Lagrange multipliers to satisfy the constraints. By varying Eq. (A.3) with respect to φ * i , we obtain that the minimizing φ i (approximate NOs) satisfy where i = λ i /n i . The optimal orbitals satisfy single-particle equations with a local potential, as intended. Next, from Janak's theorem (see also the book by Dreizler and Gross, [57] p. 54), or directly, by varying Eq. (A.3) with respect to the occupation number n i , we obtain i − µ = 0 .
(A.5) Equation (A.5) may not hold for any i. It must hold when the corresponding occupation number n i is not equal to zero or one, i.e., for any i such that 0 < n i < 1.
Using reductio ad absurdum, let us assume that the γ s (r, r ) which minimizes E v [γ(r, r )] is the 1RDM of an interacting wave function, exhibiting a large number of n i with fractional values, 0 < n i < 1 (an infinite number of fractional n i for a complete orbital basis). We are not interested in 1RDMs with just a few occupation numbers different from zero or one, which may arise when the many-body ground state is degenerate, since these 1RDMs do not correspond to interacting wave functions. It follows that Eq. (A.5) holds for all i for which 0 < n i < 1 and the single-particle Hamiltonian in Eq. (A.4) is degenerate for all these i with the same eigenvalue. However, it is unphysical (absurd) that the KS single-particle potential in Eq. (A.4), or indeed any local potential, can show such large degeneracy. We conclude that γ s (r, r ) cannot be the 1RDM of an interacting wave function. QED As with the KS total energy functional (A.1), strictly, only the infimum of the RDMFT total energy functional in Eq. (A.2) is defined. For the optimization of the ANOs leading to Eq. (A.4), where the occupation numbers are held fixed, we assume once more that, for the given set of occupations {n i }, a minimizing set of orbitals exists within the space of sets of orbitals where we search for the minimum. In this space, each set of orbitals originates from a common local potential v and the orbitals satisfy the Aufbau principle for the given set of occupations {n i }. For the special case where the fixed occupation numbers are zero or one, this assumption (as in the KS case) reduces to the familiar assumption that the density is non-interacting v-representable.