Quasi-particle energy spectra in local reduced density matrix functional theory

Recently, we introduced (e-print arXiv:1407.7128) {\em local reduced density matrix functional theory} (local RDMFT), a theoretical scheme capable of incorporating static correlation effects in Kohn-Sham equations. Here, we apply local RDMFT to molecular systems of relatively large size, as a demonstration of its computational efficiency and its accuracy in predicting single-electron properties from the eigenvalue spectrum of the single-particle Hamiltonian with a local effective potential. We present encouraging results on the photoelectron spectrum of molecular systems and the relative stability of C$_{20}$ isotopes. In addition, we propose a modelling of the fractional occupancies as functions of the orbital energies that further improves the efficiency of the method useful in applications to large systems and solids.


I. INTRODUCTION
In electronic structure theory, a desirable and elegant feature of independent particle models, such as the Hartree-Fock equations or the Kohn-Sham scheme, is the direct prediction of single-electron properties, like ionization potentials (IPs), from the eigenvalue spectrum of their corresponding effective single-particle Hamiltonians. For example, in Hartree-Fock (HF) theory, Koopmans showed 1 that the eigenenergies, i , of the occupied molecular orbitals are equal to minus the corresponding ionization potentials, I i = − i , within the approximation that the other occupied orbitals remain frozen. Also, in (exact) Kohn-Sham (KS) density functional theory (DFT) the energy of the highest occupied molecular orbital (HOMO) equals the first ionization potential of the system. 2 Further, by inverting accurate ground state densities to obtain a good approximation of the exact KS potential, it is found that occupied KS orbital energies approximate the experimental IPs of molecules much closer (∼0.1 eV difference) than those of HF (∼1 eV difference). [3][4][5] Although the question about the physical content of the KS orbitals and the KS orbital energies raised a scientific debate, 6 theoretical justification for this result was given by Baerends and co-workers 3,7 and by Bartlett and co-workers 8,9 who proved a generalization of Koopmans' theorem in KS-DFT. The KS molecular orbitals are routinely employed for chemical applications. [10][11][12] Orbital energies from local or semilocal density functional approximations (DFAs) underestimate substantially the IPs of molecules. 13,14 Nevertheless, they are still useful and the agreement with experimental IPs can be improved by applying a uniform shift 14,15 or linear scaling. 10 The wrong asymptotic behavior of the KS potential, a major deficiency of local DFAs like the local density approximation (LDA) or the generalized gradient approximation (GGA) and the manifestation of self-interactions 16 (SIs), is responsible for the large deviations of orbital energies from experimental IPs. For example, at large distances, the LDA exchange and correlation (xc) potential vanishes exponentially fast, rather than correctly as −1/r. Consequently, the electron-electron (e-e) part of the approximate KS potential decays as N/r and an electron of the system at infinity, feels the repulsion of all N electrons, itself included. The Perdew-Zunger self-interaction correction (SIC) method 17 in DFT offers a correction to this problem and was found to yield orbital energies closer to the experimental IPs 18 improving several other properties as well. [18][19][20] The GW method [21][22][23] was initially introduced to improve the obtained quasiparticle spectrum of solids but in the last decade, GW at various levels of approximations was also applied to finite systems [24][25][26][27][28][29][30] improving significantly the quasiparticle excitation energies with respect to standard DFT-approximations. Those calculations suffer from a strong initial state dependence and the good agreement found could be just fortuitous. In this context, self-consistent GW 27-29 was found to systematically improve ionization energies and total energies of closed shell systems. The single electron spectral properties are in very close agreement with experiment avoiding the starting-point dependence.
Recently, Gidopoulos-Lathiotakis 31 proposed to deal with the problem of SIs in DFAs by replacing the approximate Hartree exchange and correlation potential in the KS equations, with a different effective potential. The latter is obtained from the optimization of the same DFA energy, but is further constrained to satisfy conditions that enforce on it the asymptotic behavior of the exact KS potential. The resulting optimal potential was found to improve dramatically the agreement of orbital energies with the experimental IPs.
Reduced-density-matrix-functional theory (RDMFT) was introduced 32 as an alternative framework to DFT. In RDMFT, the one-body reduced density matrix (1-RDM) is the fundamental variable, in place of the electron density. Basic quantities associated with the 1-RDM are the occupation numbers and the natural orbitals, i.e., its eigenvalues and eigenfunctions. Several approximations for the total energy as a functional of the 1-RDM -or usually in terms of the occupation numbers and the natural orbitals -have become available. [33][34][35][36][37][38][39][40][41][42][43][44] They have proven to describe correctly many diverse properties such as molecular dissociation [35][36][37][38][39] or bandgaps. [44][45][46][47] However, so far, the computational cost has restricted the applications of RDMFT to prototype systems. Most of the computational expense is due to the determination of the orbitals which are not obtained from an eigenvalue equation but through a numerically expensive minimization. In contrast to DFT, in RDMFT there is no KS noninteracting system with the same 1-RDM as the interacting system. Thus different approaches have to be considered in order to define effective single-particle Hamiltonians. [48][49][50] Recently, we proposed local-RDMFT, 51 a theoretical framework that incorporates static correlation effects in the single-particle, Kohn-Sham equations. It is based on the adoption of RDMFT approximate functionals (optimized with fractional occupations) for the exchange and correlation energy, together with a search for an effective local potential, whose eigen-orbitals minimize the total energy. The search of the effective potential is performed as in Ref. 31, where, apart from correcting possible SIs, it was also found to avoid mathematical pathologies of finite-basis optimized effective potential (OEP). 52 Local-RDMFT can be viewed within the framework of the OEP method in DFT, where the correlated exchange and correlation (xc) functionals from RDMFT allow us to go beyond the level of an exchange only OEP (x-OEP) calculation. 53 Equally, local-RDMFT can be regarded an approximation in RDMFT, employing an effective single particle scheme to generate the approximate natural orbitals (ANOs). Thus, local RDMFT provides an energy eigenvalue spectrum directly connected to the ANOs and as we find in Ref. 51, this energy spectrum reproduces the IPs of small molecules in closer agreement with experiment than HF Koopmans'. In addition, it allows us to calculate accurately total energies at any geometry, from equilibrium all the way to the dissociation limit, which is well described without the need to break any spin symmetry.
In the present work, we demonstrate the efficiency of local-RDMFT by applying it to larger molecules. More specifically, we calculate the IPs of systems of ∼20 atoms and compare them with experiment. For some aromatic molecules, we compare the calculated orbital energies with the peaks of the corresponding photoelectron spectra (PES). We also study the relative stability of C 20 isomers and show that systems of this size are within the reach of our method. Finally, we propose that, in local-RDMFT, the optimization of the fractional occupations can be simplified by modelling them in terms of the orbital energies. We expect that such ideas will be very useful in the application to larger systems and solids.
In Sec. II, we summarize the basics of local-RDMFT and in Sec. III we discuss our results on the application to the C 20 isomers, the IPs of molecular systems and the comparison of the calculated orbital energies with the PES of aromatic molecules. Finally, in Sec. IV, we demonstrate that the optimization of fractional occupation numbers can be simplified by modeling their dependence on single particle energies.

II. LOCAL-RDMFT
Local-RDMFT combines two main features: (i) the nonidempotency of the optimal 1-RDM where the fractional occupation numbers are provided by minimizing the total energy functional under Coleman's N-representability conditions and (ii) the incorporation of a single particle effective Hamiltonian with a local potential. As we showed in Ref. 51, one has to depart from xc functionals that are explicit functionals of the electronic density alone, since they lead to idempotent solutions. Thus, we have to adopt either explicit functionals of the 1-RDM, or functionals of the orbitals and the occupation numbers.
The central assumption in local-RDMFT 51 is that the search for the set of optimal ANOs is restricted in the domain of orbitals that satisfy single-particle equations (KS equations) with a local potential. The search for the e-e repulsive part V rep (r) of the effective local potential (the analogue of the Hartree-exchange and correlation potential in the KS equations) is effected indirectly, by a search for the effective repulsive density (ERD) ρ rep (r) whose electrostatic potential is V rep (r), i.e., Additionally, following Ref. 31, two constraints are imposed in the minimization with respect to ρ rep (r): The first condition is a property of the exact KS potential and the x-OEP potential, while the second is a condition that gives physical content to the ERD as a density of N − 1 electrons. It is unknown if (3) is a property of the exact KS potential or of x-OEP but without it, the search for ERD is mathematically ill posed for finite basis sets. The two conditions, (2), (3), together lead to physical solutions. The optimal ERD and the effective local potential can be obtained, similar to the OEP method, by solving the integral equation 51 The response function χ (r, r ) and b(r) are given by with F (j ) Hxc defined by δE Hxc δφ * j (r) E Hxc is the approximation for the e-e interaction energy, φ j are the ANOs and n j , j their corresponding occupation numbers and orbital energies (eigenvalues of the effective hamiltonian). The two constraints can be incorporated with a Lagrange multiplier (2) and a penalty term (3) that introduces an energy cost for every point r where ρ rep (r) becomes negative. Terms over pairs of orbitals with almost equal occupations cause numerical instabilities in the sums of Eqs. (7) and (8) and we have decided to exclude them by introducing a small cutoff n c . The reader is referred to the discussion in Ref. 51. Our choice affects mostly pairs of weakly occupied orbitals, whose energies are in any case inaccurate for finite localized orbital basis sets, as occasionally they violate the aufbau principle and the negative definiteness of χ . For very small cutoff n c we observe convergence problems, mainly while attempting to enforce the positivity constraint (3). When n c is large enough to exclude erroneous terms involving weakly occupied orbitals (past a typical value ∼0.1), convergence issues improve dramatically and IPs remain unchanged for a broad range of values of n c . We have found that a choice for n c ∼ 0.1 − 0.3 leads to stable solutions where the IPs are insensitive to a change of n c .
The ANOs are expanded in a basis set (orbital basis), while the ERD in a separate (auxiliary) basis and Eq. (4) transforms into a linear system of equations. This linear system typically becomes singular and we use a singular value decomposition (SVD) to obtain smooth and physical densities and potentials (see Refs. 31 and 51). We note that a SVD for the matrix of the density-density response function may introduce singular behavior in the effective potential. 52 However, the two constraints (2), (3) reduce the variational freedom in the space of effective local potentials V rep (r), to such a degree that a discontinuity correction 52 in the null space of the (finite-orbital-basis) response function is no longer necessary.
We stress that there is no functional derivative relation linking our local effective potential with the total energy functional. As a result, in our formulation, we avoid the collapse of all eigenvalues corresponding to fractional occupations to the chemical potential. In addition, the effective local potential in local-RDMFT differs in general from the exact KS potential. However, a comparison with the exact KS potential is still meaningful and establishes the physical significance of the approximation.

III. APPLICATIONS
We applied local-RDMFT to the molecular systems shown in Table I, employing several approximate RDMFT functionals: Müller, 33, 35 the third Buijse-Baerends corrected (BBC3) approximation, 36 the Power functional 37,44 and the Marques-Lathiotakis (ML) approximation. 40 Values obtained with HF Koopmans' are also included. The cc-pVDZ and the uncontracted cc-pVDZ basis sets were employed for the orbital and the auxiliary basis sets, respectively, for all calculations. Our target is to examine the usefulness of the quasiparticle energy spectrum, i.e., how close the orbital energies are (in absolute value) to the corresponding IP. Only the first IPs are shown in Table I Fig. 1, the present results for larger systems are compared to those in Ref. 51 by plotting the absolute, percentage error of the IPs for the four different RDMFT energy functionals and HF. We find a remarkable agreement between the energy eigenvalues and the experimental IPs for the functionals we tested. For the small systems, all errors are below or around 5%, with the ML functional as low as ∼3%. We should emphasize that this agreement is found only if the positivity condition of Eq. (3) is enforced, otherwise the error increases to 20% for some functionals. The agreement with experiment is even better for the first IPs for small systems with the Müller functional the most accurate in this case with an error of only ∼2%. The trend is similar for the larger systems of Table I where the average error is relatively larger for all functionals except for the Müller functional, for which it remains below 3%. Overall the agreement with experimental values is very good and substantially better than HF Koopmans'. The good performance of the Müller functional is rather uniform. In contrast, the results of BBC3 deviate from experiment mainly for the IPs of alcanes, raising the error to ∼6.5%.
As an additional test for the physical interpretation of the obtained quasiparticle energy spectrum, we plot (see Fig. 2) the orbital energies of local-RDMFT together with the experimental PES for three aromatic molecules, benzene, naphthalene, and anthracene. For all systems, we show the eigenvalues obtained with the Müller functional as this functional is found to yield better results. For benzene, we also include the BBC3 and power functional eigenvalues. The agreement with the experimental spectrum is fair for all three functionals, however, the Müller orbital energies are more accurate. Especially for naphthalene, the Müller eigenvalues are in excellent agreement with experiment.
For comparison, for the LDA and GGA approximations, the errors of the KS eigenvalues compared with experimental IPs can be of the order of 30%-40% due to SIs. 30,31 Theinclusion of a percentage of HF exchange in hybrid functionals reduces SIs which, however, remain large. As an example, we applied the Becke 3 parameter exchange-correlation functional 55 (B3LYP) to the systems in Table I (using the same geometries and basis set) and the absolute, percentage error of the IPs is 26%. GW method applied on top of plain DFT and HF calculations is found to improve substantially the quasiparticle spectrum. To summarize a few applications: Blase et al. 24 obtained the IPs of photovoltaic-relevant molecules with a mean average error of 3.8%. Van Setten et al. 30 Table I compared with experimental results. the performance of a hierarchy of GW approximations for benzene, pyridine, and the diazines and compared the quasiparticle spectrum with PES. Caruso and co-workers 27-29 developed an all electron implementation of self-consistent GW (sc-GW ) with localized basis functions and showed that it is more accurate than other approximations lower in GW hierarchy. They applied sc-GW on five molecules relevant for organic photovoltaics 28 obtaining an average error of 0.4 eV (maximum error 1.2 eV).
The GW calculations come with a high computational cost despite employing routinely efficiency improving techniques like the resolution of identity. 56 Local-RDMFT on the other hand is more efficient method and the efficiency can further improve by adopting techniques like the RI. More importantly, it provides IPs of similar quality as GW approaches. Although we cannot make a quantitative comparison since we used different set of systems and basis sets, as we see in Table I, the root-mean-square deviation from experiment for the local Müller functional is as low as 0.32. This quantity has values similar to those reported for the GW approaches for all the functionals we employed.
As a demonstration of the efficiency of local-RDMFT, we calculated the three most stable isomers of C 20 namely the cage, the bowl, and the monocyclic ring structures using several functionals. These three isomers are energetically very close and the predicted most stable isomer differs from method to method. For example, the ring is found the most stable by HF 60 and DFT-GGA 61 and B3LYP, 62 the cage by DFT-LDA 63 and CCSD 64 and the bowl by a more recent CCSD 65 and quantum Monte Carlo calculations. 66 CCSD and QMC are the most accurate schemes, and there is a consensus that the bowl is the most stable structure with the cage being almost isoenergetic. Our local-RDMFT results for the total energies of the three isomers are shown in Table II. We employed cc-pVDZ and uncontracted cc-pVDZ as orbital and auxiliary basis sets, respectively, and the optimal geometries that were obtained at the MP2 level of theory. Apart from the approximations considered above in this application we also employed the functionals of Goedecker-Umrigar (GU), 34 the The results presented in this section demonstrate that the local-RDMFT formalism preserves the advantages of RDMFT in calculating correlation energies and, as we showed in Ref. 51, also in describing molecular dissociation. In addition, it provides energy eigenvalues which are in very good agreement with experimental IPs. Compared with standard RDMFT, the significant reduction in computational cost allows for applications to larger systems previously inaccessible to this theory.

IV. MODELLING THE FRACTIONAL OCCUPANCIES
Fractional occupation numbers are usually employed in DFT calculations in an ad hoc way to introduce temperature effects and to help the convergence of the self-consistent KS-equations loop in small-gap or metallic systems. In the case of local-RDMFT, fractional occupations are introduced naturally through an optimization procedure. The existence of fractional occupations and at the same time of a corresponding single electron energy spectrum allows for the modelling of the occupation numbers as functions of the energy eigenvalues. This modelling is no longer arbitrary and is implemented through the functional optimization. For instance, one can assume a smooth parametric form for the function n j ( j ) connecting the occupation numbers to singleparticle energies and optimize the model parameters such that the energy functional is minimized. The advantage is that occupation numbers are obtained in a simpler minimization procedure of a few variables only. As a demonstration we consider The parameters μ, β s , and β w are optimized by minimizing the energy functional with respect to them for a given set of orbitals. A Fermi distribution modelling of the occupation numbers as functions of the eigenvalues of a Hamiltonian with a local potential was also introduced by Grüning et al. 69 However, in that work, the model parameters were not optimized iteratively with the orbitals for each calculation but were chosen universally such that the obtained dissociation of H 2 molecule is as close to the exact as possible. In Fig. 3, we show the energy eigenvalues obtained for naphthalene compared with the PES and the eigenvalues of the standard local-RDMFT. As we see the obtained spectrum of such a model is reasonable. We believe that such a procedure will be useful in the application to periodic systems, especially metals, simplifying the optimization of the occupation numbers, offering a natural way to introduce occupation smearing quasiparticle renormalization factors and accounting for quasiparticle renormalization effects in the homogeneous electron gas case. 70

V. CONCLUSION AND PERSPECTIVES
Local-RDMFT, a novel scheme that incorporates static correlation in the KS equations and allows the accurate description of molecular dissociation, was applied to molecular systems of size up to 20 atoms.
The new approach associates a quasiparticle energy spectrum to the ANOs. This spectrum is in good agreement with experimental IPs and PES for molecules. The reduction in computational cost, permitted, for the first time, the calculation of larger molecules with the improved accuracy of RDMFT functionals. To demonstrate the efficiency of the new scheme, we applied it to the three most stable C 20 isomers although the tiny energy differences of these systems are probably beyond the accuracy of current RDMFT approximations.
The new method provides a powerful tool which opens a new avenue for bringing the advantages of RDMFT into DFT. Due to the similarity of the local-RDMFT and the OEP equations, the systematic and physically motivated approximations in density-matrix based schemes to cope with strongly correlated systems 44 and static correlation can now easily be brought to the realm of DFT. For the first time, a method is able to simultaneously describe ground-state properties, bond-breaking and photoelectron spectra.
Compared with orbital-dependent functionals in DFT, the additional cost in local-RDMFT comes from the iterative optimization of the occupation numbers and the ANOs. This extra cost can be reduced by connecting the occupation numbers directly to the energy eigenvalues through physically motivated models, see Eq. (10).
In the future, the method can be extended to the timedependent regime with the aim to provide more accurate energy spectra and description of electronic excitations. The development of a linear-response formalism will in addition give access to a large number of experimentally measurable properties.