On the performance of natural orbital functional approximations in the Hubbard model

Strongly correlated materials are now under intense development, and natural orbital functional (NOF) methods seem to be able to capture the physics of these systems. We present a benchmark based on the Hubbard model for a class of commonly used NOF approximations (also known as reduced density matrix functional approximations). Our findings highlight the importance of imposing ensemble N-representability conditions in order to obtain consistent results in systems with either weak or strong electronic correlation, such as the Hubbard system with a varying two-particle interaction parameter. Based on the accuracy of the results obtained using PNOF7, which retrieves a large amount of the total strong nondynamic correlation, the Hubbard model points out that N-representability gives solid foundations for NOF development.

the PNOF7 natural orbital occupancies (figure 6) are not always adapted to the symmetry of D 4h point group. PNOF7 yields two configurations of orbitals depending on the correlation regime. In the region U/t 2.7, natural orbitals remain adapted to the symmetry, however, for U/t > 2.7 PNOF7 leads to localization in the natural orbital representation. Orbitals adapted to the symmetry can be retrieved by a unitary transformation that leads to the canonical representation as we have shown for six site hexagon Hubbard model in the article.
For the six site hexagon, no significant changes appear for PNOF7 double occupancy (figure 5) and natural orbital occupancies ( figure 6).
Finally, for the ten site Hubbard model including an Aubry-André on-site potential, the MBB and CGA energies have been corrected. In contrast to figure 10 of the article, MBB and CGA are able to reproduce the asymptotic behavior as U/t gets larger, although, energies are always below the FCI ones (figure 10). We thank Mr Robert Schade for his kind comments regarding the MBB energies.
The theorem on the existence of the functional [1,2,4] of the 1RDM does not provide the means for its construction. Exact properties of the hole function [18,19] have been used to develop NOFAs mantaining the general form of the Hartree-Fock exchange term. These approximations include Coulomb and exchange type integrals, so they are denoted as J K-only NOFAs. Alternatively, Piris has proposed [11] employing the known ensemble N-representability conditions of the two-electron reduced density matrix (2RDM) [20] and progressively enforcing them to generate improved functionals. The latter includes exchange and time-inversion integrals [21], so they must be referred to as J KL-only NOFAs.
Recently the interest related to the Hubbard model has grown in electronic structure theory [22][23][24][25][26]. This system can be viewed as the simplest possible model of correlated fermions, since the compromise between kinetic energyincluding via particle hopping-and Coulomb repulsion is included in the Hamiltonian. Thus, the Hubbard model exhibits magnetic ordering, metal-insulator transition, superconductivity, and Tomonaga-Luttinger liquid in 1D, among others, so making it possible to study many properties in strongly correlated materials. In other words, this minimal model captures the basic nature of electron correlation, and offers full tunability in order to explore different correlation regimes. As such, the Hubbard model is an ideal candidate for calibration and benchmarking of approximate electronic structure methods.
Strongly correlated systems are challenging for theoretical methods, since the independent electron picture fails to describe them. Variational wavefunction approaches such as Gutzwiller [27,28] and Baeriswyl [29][30][31] have proven to be able to describe fermionic lattice models, and together with DMRG or quantum Monte Carlo methods, they are commonly employed to study the fundamental properties of strongly correlated lattice systems. An attempt to use another approach based on RDMFT is done here-in fact, apart from the good performance shown by variational 2RDM calculations [32,33], NOFAs seem also to be suitable to give an accurate description of such systems [34,35]. Therefore, the objective of the present paper is to use the Hubbard model as a stringent validation tool for commonly used NOFAs. In particular, the one-dimensional Hubbard model with and without external potential is studied, focusing on energies, double occupancy, and natural ONs. The latter are an indicator of electron correlation [25], so together with the energy results they show how accurately the Hubbard model is being described. The paper is organized as follows. In section 2.1 we describe the main features of NOF theory and we present the approximations employed in this work. Then we give a brief description of the Hubbard model, and we describe the properties that will be studied throughout the paper. Results obtained using a set of popular NOFAs are shown in section 3, together with the exact results obtained from FCI calculations.

Natural orbital functional theory
The Hamiltonian of a many-electron system is the sum of oneand two-particle operators, hence the energy is determined exactly by the 1 and 2 RDM, denoted hereafter by Γ and D, respectively. The exact functional can be used with an approximate 2RDM built from 1RDM to obtain an energy functional E [Γ]. The major advantage of this formulation with respect to the density functional theory (DFT) is that the kinetic energy does not require the construction of a functional. The unknown functional only needs to incorporate electron correlation. This is, however, a formidable task-since the functional N-representability [36], which refers to the conditions that guarantee the one-to-one correspondence between E [Ψ] and E [Γ], is a problem related to the N-representability of the 2RDM [20].
In NOF theory [37], the spectral decomposition of the Here, H ii is the one-particle part of the Hamiltonian involving the kinetic energy and the external potential operators, whereas kl|ij are the matrix elements of the two-particle represents the reconstructed 2RDM from the ONs. We neglect any explicit dependence of D on the NOs themselves because the energy functional (1) already has a strong dependence on the NOs via the one-and two-particle integrals.
It is noteworthy that the resulting approximate functional alone can be implicitly dependent on Γ [3], since the Gilbert's theorem [1] on the existence of the explicit functional E [Γ] is valid only for the exact ground-state energy. In this vein, NOs are the orbitals that diagonalize the one-matrix corresponding to an approximate energy that still depends on the 2RDM. A detailed account of the state of the art of the NOF theory can be found elsewhere [38,39].
We consider S z eigenstates, so only RDM blocks that conserve the number of each spin type are non-vanishing. Specifically, the 1RDM has two nonzero blocks Γ αα and Γ ββ , whereas the 2RDM has three independent nonzero blocks: D αα , D αβ , and D ββ . The parallel-spin components of the 2RDM must be antisymmetric, but D αβ presents no special symmetry [37].
We address only singlet states in this work, so the parallel spin blocks of the RDMs are equal. A spin-restricted formalism is adopted; accordingly, the NOs have identical spatial parts and equal ONs f i = f α i = f β i , with 0 f i 1 and 2 i f i = N. The latter represents the necessary and sufficient conditions for ensemble N-representability of the 1RDM [40].
The Hartree-Fock (HF)-like approximation for the 2RDM is defined as Some approximations to the 2RDM are produced by replacing only the σσ-elements of the HF like approximation, viz.
The F( f i , f j ) functions are collected in table 1. These approximations lead to J K-only NOFAs, since the electronic energy involves only the Coulomb (J ij = ij|ij ) and exchange integrals (K ij = ij| ji ). Note that equation (3) violates the antisymmetric requirement unless F = f i f j ; consequently, none of these functionals affords an N-representable 2RDM. Extensive N-representability violations for these NOFAs have been reported [41,42].
The use of the (2,2)-positivity N-representability conditions, also known as P, Q, and G conditions, for generating an N-representable 2RDM was proposed by Piris [11]. This Table 1. F( f i , f j ) functions used for 2RDM reconstruction in equation (3).
particular reconstruction is based on the introduction of two auxiliary matrices and Π, expressed in terms of the ONs, to reconstruct the cumulant part of 2RDM [43], viz.
Both auxiliary matrices are constrained to certain bounds due to the enforced positivity conditions, symmetric properties and sum rules satisfied by the 2RDM. The conservation of the total spin makes it possible to derive the diagonal elements ∆ ii = f 2 i and Π ii = f i [44]. Appropriate approximations for off-diagonal elements lead to different implementations of the NOF known in the literature as PNOFi (i = 1-7). The performance of these functionals is comparable to those of the best wavefunction-based methods, and has been recently reviewed in [45]. In the present study, the latest functional of this series is used.
PNOF7 [17] corresponds to an interacting electron-pair model. Let us divide the spatial orbital space Ω into N/2 mutually disjoint subspaces Ω g , so each orbital belongs only to one subspace. Consider each subspace contains one orbital g below the Fermi level ( N/2), and N g orbitals above it, which is reflected in additional sum rules for the ONs: Taking into account the spin, each subspace contains solely an electron pair, hence the normalization condition for the 1RDM is automatically fulfilled. It is important to note that orbitals satisfying the pairing conditions (5) are not required to remain fixed throughout the orbital optimization process [46]. Given the functional form of the auxiliary matrices and Π [17], the energy (1) of the PNOF7 can be conveniently written as where From equation (6), we can observe that the energy for singlet states contains not only J and K integrals, but also the exchange and time-inversion integral L ij = ii| jj [21], so PNOF7 is a J KL-only approximation. The first term of the equation (6) is the sum of the pair energies described accurately by the two-electron functional E g . In the second term, E pq correlates the motion of the electrons in different pairs with parallel and opposite spins. It is clear that the main weakness of the approach (6) is the absence of the inter-pair dynamic electron correlation, since Π inter ij has significant values only when the ONs differ substantially from 1 and 0. Consequently, PNOF7 is expected to be able to recover the complete intrapair, but only the nondynamic inter-pair correlation.
On the other hand, when the Π inter ij elements are neglected, equation (6) yields the independent pair model PNOF5 [14,15]. Interestingly, an antisymmetrized product of strongly orthogonal geminals wavefunction, with the expansion coefficients explicitly expressed by the ONs, also leads to PNOF5 [15]; therefore, it is strictly N-representable. Indeed, 1RDM functional theory can also be connected to the geminal functional theory [47,48].
Minimization of the energy functional equation (1) is performed under the orthonormality requirement for the NOs, and the particle number fixed to N, by using the Lagrange multipliers λ and μ respectively. Occupancies are expressed by the free variable γ with f i = cos 2 γ i to assure 1RDM ensemble N-representability conditions [40]. Thus, we can define an auxiliary functional An iterative diagonalization procedure [46] has been employed to make L stationary with respect to variations in the ONs and the NOs separately, which is based on the hermiticity of the matrix of the Lagrange multipliers λ at the extremum. Basically the 1RDM and λ − λ † , where † denotes the conjugate transpose, can be brought simultaneously to a diagonal form at the solution, hence the set of NOs is given by those orbitals that make λ Hermitian. Occupation optimization is done by sequential quadratic programming (SQP) for all J K-only NOFAs. Since both PNOF5 [14] and PNOF7 [17] are electron-pairing approaches, the particle number is automatically fixed to N, so a Conjugate Gradient (CG) algorithm is used to solve the unconstrained occupation optimization in the case of these NOFAs, which yields substantial savings of computational time.

Hubbard model
The Hubbard Hamiltonian is possibly the simplest prototype for modeling strongly correlated systems. This model has been largely used to benchmark electronic structure methods [22][23][24][25][26], but also to describe the electronic properties of many materials, e.g. metal-insulator transitions, charge-and spindensity waves in superlattices, etc. In one dimension and standard notation the Hubbard model Hamiltonian reads where greek indexes μ and υ denote sites, µ, υ indicates only near-neighbors hopping, t > 0 is the hopping parameter, fermionic creation(annihilation) operator, v µ,σ is the on-site energy and U is the site interaction parameter. Hereafter, we use the term homogeneous Hubbard whenever v µ,σ = 0; ∀ {µ, σ}, whereas inhomogeneous Hubbard is used whenever nonzero external potential is present. Thus, for the Hubbard dimer it will be set v SA,σ = −v SB,σ in section 3.3, so that it gives rise to a nonzero potential ∆v. An even more interesting external potential will be used to study inhomogeneous systems for the ten site model, namely the one-dimensional Aubry-André model [49] including electron-electron repulsion where V is the modulation amplitude of the on-site potential, α determines the periodicity, and δ fixes the modulation phase. This model, which is intimately related with the Harper model, has been used to explore topological properties in 1D systems, among others. The electron-electron repulsion is extremely local in the Hubbard model, and can be tuned by the parameter U. Although additional complexities can be included by setting variable parameters U µ and t µ [22], let us fix t and vary only the particle-particle interaction U in order to cover different correlation regimes, thus U/t will be used as a dimensionless measure for the relative contribution of both terms. In the limit U/t → 0, the system can be described by mean-field theories-hence the Hartree-Fock approximation recovers the exact wavefunction for the symmetric Hubbard dimer if U = 0; however, as long as U/t gets larger the system is nevermore weakly correlated, and methods including electron correlation are needed to give an accurate description of the model.
A general solution for the Hubbard model requires a numerical treatment. We will restrict attention to half-filling cases throughout this paper, so there is on average one particle per site in the model, such that corresponding distribution among the sites depends on the correlation regime and the number of sites of the Hubbard model. Overall, if the hopping parameter is larger than the on-site interaction the electrons tend to occupy doubly the sites, while at the U/t 1 limit, also known as the strong correlation limit, electrons try to keep away from each other by half-filling the sites, which corresponds to the Mott-Hubbard regime [13,24]. This transition is quantified by the natural orbital occupancies of the system, f i , which actually are an indicator of correlation [25]. Our analysis will focus on the latter, but also on E/t values and on the double occupancy of the sites dE/dU, which is defined as

Results and discussions
The Hubbard model is here exploited in order to clarify the performance of commonly used NOFAs in many correlation regimes (specially in the strong correlation regime). The present paper focuses on ground state E/t values, double occupancy (equation (11)), and natural occupancies in comparison with the exact result; which are obtained from FCI calcul ations for a range of U/t values in order to cover all correlation regimes in each case. To this purpose we have used a modified version of the code developed by Knowles and Handy [50,51]. 1RDM and 2RDM have been calculated from the FCI expansion coefficients using a homemade code, the DMN code, developed by Matito and Feixas [52]. The results presented here for all NOFAs have been computed using the DoNOF code developed by Piris and co-workers. Note that double occupancy of sites has been numerically evaluated by using the five-point formula where the step size is set to h = 0.001 (the error within this approximation is of the order of h 4 ).

Exact results for homogeneous Hubbard model
In this section we discuss several properties of the homogeneous two site, four site square, and six site hexagon Hubbard systems, which have been previously employed to study many electronic structure methods [22][23][24][25][26], [53,54]. Energy values for the homogeneous 10 sites, 12 sites, and 14 sites are also included. Exact FCI E/t values for a range of U/t values corre sponding to these systems are shown in figure 1. There is a similar trend for the three curves corresponding to the graph on the left-all of them show negative E/t values at the zero correlation point (U = 0), and converge asymptotically to E/t = 0 at the strong-correlation limit (U/t → ∞). As expected, the absolute energies are larger as the number of sites is increased. Interestingly, relative differencies are smaller in the case of 10, 12, and 14 site systems, and the asymptotic limit is located at larger U/t values when more sites are added to the model.
In order to get a more reliable indication of the energy, the derivative with respect to parameter U/t is also studied for the two, four, and six site systems (figure 2), since physical interpretations can be obtained due to their conceptual simplicity. In fact, this is a measure of double occupancy of the sites according to equation (11). According to figure 2 the double occupancy is maximum in the weak correlation region, since the exact ground-state wave function is recovered by an independent particle model when there are no two-particle interactions. The population of the sites spreads out as the correlation increases, so for large U values the double occupancy tends to zero due to particle-particle repulsion. Neither the energy nor dE/dU show significant differences when the number of sites varies in the Hubbard model, but the study of the ONs displays another situation.
The half-filled two site Hubbard model is the simplest case that we have studied. It is known that when there is no correlation (U = 0) the exact wave function for the ground-state may be written as a single Slater determinant |Ψ = |σ 2 g built using a σ g orbital; labeled attending to D ∞h symmetry of the system. The σ g orbital is defined as and is also known as the bonding orbital in quantum chemistry. Note that sites are labeled by S X , where X is replaced by a letter in alphabetic order. Nevertheless, when correlation increases the single Slater determinant picture no longer holds and the wave function may be written as |Ψ = C g |σ 2 g + C u |σ 2 u , where the expansion coefficients C g and C u are determined variationally, and the second Slater determinant is built using the σ u orbital This orbital is also known as the antibonding orbital in quantum chemistry. The σ g and σ u orbitals form a basis which is adapted to the symmetry of the D ∞h point group. These orbitals are also the NOs for the system, since the 1RDM obtained from the wavefunction is already diagonal in this basis for the homogeneous two site Hubbard model. Corresponding ONs are shown in the top of figure 3. Note that coefficients C g and C u are equalized in the strong-correlation limit (U/t → ∞), so the occupancies of the NOs become equal to one (spins summed). In other words, the 1RDM is diagonal and equal to the identity matrix (this limit resembles the H 2 molecule in the dissociation limit for a minimal basis). For the four site square Hubbard model the orbitals adapted to the D 4h symmetry also correspond to the NOs; these are the non-degenerate the two-fold degenerate e g orbitals and the non-degenerate According to the natural occupancies shown in the middle of figure 3, both e g orbitals are half-occupied independently of the interaction strength U, whereas orbitals a 1g and b 1g play the role of the σ g and σ u orbitals in the two site model. Regarding the orbitals adapted to the symmetry of the system in the case of the six site hexagon Hubbard model, which also correspond to the exact NOs, these orbitals read as the two-fold degenerate e 2g orbitals the two-fold degenerate e 1u orbitals and the non-degenerate The plot in the bottom of figure 3 shows the exact NO occupancies. The behavior is similar to that obtained for the two site model; at the weak correlation limit the double occupancy of NOs is maximum, while all orbitals become half-occupied in the limit U/t → ∞.

NOF results for homogeneous Hubbard model
In this section, the Hubbard model is used as a stringent validation tool for commonly used NOFAs, which have been described in section 2.1. Note that a spin-restricted formalism has been employed for all calculations, in order to hold S 2 = 0 and S z = 0. Here, we will use the adaptedto-symmetry basis set described in the previous section for each model, which is conceptually similar to a molecular-like basis set. Interestingly, the latter correspond to the NOs for the homogeneous two site and four site square systems, independently of the interaction strength U/t and the NOFA. Figure 4 shows the differences obtained for E/t values with respect to FCI (∆ [E/t] = E NOFA /t − E FCI /t ) for a wide range of U/t values. As we observe from the left-top of figure 4, only MBB, PNOF5 and PNOF7 functionals reproduce the exact FCI energy for the Hubbard dimer (the result for MBB has been obtained previously [26]). Apart from MBB, PNOF5 and PNOF7; CA and CGA functionals show a good asymptotic behavior in the U/t → ∞ limit. Conversely, neither restricted Hartree-Fock (RHF) nor the Power functional recover the exact energy or attain a good asymptotic behavior, though RHF shows a good performance at the weak correlation limit. (Results corresponding to some other NOFAs like the Goedecker and Umrigar [7], Marques and Lathiotakis [12], and Gritsenko, Pernal and Baerends [10] have not been included, due to showing catastrophic performance even for the two site case.) MBB and Power approximations are known to violate the N-representability P-, Q-, and G-conditions [20], as some of us have recently shown [42], even for the twoelectron case. E/t values corresponding to the Power functional differ from those reported by Kamil et al (see figure 1 in [26]). The difference between Kamil's et al energies and ours is due to the fact that we stick to a spin restricted formalism for our optimizations of the orbitals, so the σ g and σ u basis is not altered. Therefore, the orbitals for α and β electrons were the same. Of course, the usage of an unrestricted formalism improves the E/t values, as Kamil et al show, but the price that we have to pay for using an unrestricted formalism is that S 2 = 0, since a mixture between the singlet and triplet solutions may arise (also known as the loss of the exact nonmagnetic character). In order to conserve S 2 = 0, we have to use a restricted formalism and expect that the accuracy of the NOFA produces the correct asymptotic behavior in the U/t → ∞ limit (which is not the case of the Power functional, as illustrated in the left-top of figure 4).
The robustness of well-behaved functionals has been tested beyond the dimer using a system with more degrees of freedom, such as the half-filled four site square Hubbard model. A similar problem, the H 4 molecule, has recently been studied in a NOFT context [55]. Contrary to what happened for the two site homogeneous Hubbard model (where intra-pair nondynamic correlation effects were dominant), there is correlation between pairs in the four site square Hubbard model. Therefore, a good functional for the four site square homogeneous Hubbard model should have a good description of nondynamic intra-pair correlation effects, but also a reasonable description of inter-pair correlation. Similar to the results obtained for the two site Hubbard model, E NOFA /t − E FCI /t values displayed in the left-middle of figure 4 obtained with RHF and Power present bad asymptotic curves, so they do not correspond to a solution of the problem. Conversely, all other NOFAs show proper behavior for small and large U/t values. PNOF7 is in outstanding agreement with FCI for all correlation regimes, whereas CGA differs from the exact curve in the interval 0 < U/t < 10. PNOF5, MBB, and CA yield qualitatively correct curves, but their corresponding energies are not as accurate as those obtained with CGA and PNOF7. The importance of inter-pair correlation must be emphasized in order to describe this system properly, because the only difference between PNOF5 and PNOF7 is exclusively due to the addition of a term to account for inter-pair electron correlation.
The half-filled six site hexagon Hubbard model presents an added complexity with respect to the two and four site systems. This model has been previously used as a benchmark, e.g. for testing new approaches based on coupled cluster methods like CCD0 in order to retrieve nondynamic correlation effects [23]. PNOF7 provides the most accurate total energies as illustrated in the left-bottom of figure 4, followed closely by CGA. Both functionals yield energies slightly below the FCI energy beyond a certain U/t value, so even PNOF7, which is built imposing the analytic necessary N-representability conditions of the 2RDM [41,20], is not stricly N-representable in constrast to PNOF5. CA, MBB, and PNOF5 show correct behavior in the weak and strong correlation limits; however, they yield poor results in the intermediate region.
Within the limitations imposed by the exact diagonalization calculation, results for larger systems have been obtainednamely, the 10, 12, and 14 site homogeneous Hubbard systems. Thus, results obtained for PNOF7, PNOF5, CGA, and MBB approximations are shown in the right-hand plots of figure 4, whereas other NOFAs have been omitted due to present bad performance for tiny systems. Contrary to the errors shown for smaller systems, CGA yields larger errors going from 10 to 14 sites, as MBB, and both of them show energies below the exact one throughout all energy curves, so the variational principle is strongly violated. PNOF approximations do not breakdown for large systems, and in spite of providing an error larger than CGA in absolute terms near the U/t ≈ 50, they approach closely to the exact result from above, for any correlation regime.
Let us focus now on the results obtained for the double occupancy in the case of the two, four and six site models. Since PNOF7 and CGA yield best energies, dE NOFA /dU − dE FCI /dU differences obtained employing these approximations are plotted in figure 5. Results obtained by using RHF and PNOF5 are also included for comparison. Note that RHF gives a constant value for each model independently of the site interaction strength U. PNOF7 is the only functional able to go parallel to the exact dE/dU according to figure 5, whereas CGA shows large errors and discontinuities for the two and four site systems at small U/t values, which is the region where it goes noticeably below the exact energies and, as shown below, CGA gives wrong pinned natural ONs. The independent pair model PNOF5 fails significantly for systems beyond two particles, specially in the four site case, due to the missing inter-pair correlation. This functional shows a behavior similar to RHF for small U values, which can be associated with the result obtained for the natural ONs (see below).
NO occupancies obtained for FCI, CGA, PNOF5, and PNOF7 in the four site square and six site hexagon Hubbard systems are plotted in figure 6. Results corresponding to the Hubbard dimer have been omitted since MBB, CGA, PNOF5 and PNOF7 reproduce the exact occupancies. Regarding the result obtained for the four site model, only PNOF7 is able to provide precise ONs for any interaction strength, in contrast to the rest of NOFAs, which show good asymptotic behavior in the U/t → ∞ limit but tend to stick to f a1g = 2.0 (spins summed) occupancies in the weak correlation region related to the bad performance obtained for dE/dU in the same correlation region. CGA exhibits the largest occupancies pinned to 2.0 in the a 1g orbital but the changes of b 1g occupancy, so it shows a wrong break of the unitary occupation of the twofold degenerate e g orbitals that remain fixed in the exact solution. Pinned natural ONs are intimately related with the wrong description of the double occupancy. According to equation (11), a bad description of the latter is due to an incorrect 2RDM needed to evaluate Ψ|n µ,↑ n µ,↓ |Ψ in the sites basis, so the double occupancy of the sites computed from first derivative of the energy also reflects the pitfalls in the 2RDM built from pinned to 2.0 ONs.
The situation is completely different for the six site model, for which NOs corresponding to the ground state for both PNOF5 and PNOF7 depend on each approximation, whereas the NOs are still those adapted to the symmetry for the J K-only NOFAs. Thus, PNOF NOs present three-fold degeneracy and thereby are not related to D 6h point group symmetry; conversely, there is a double degeneracy of the occupations associated with the e 2g and e 1u orbitals for the rest of approximations. Similarly to the results obtained for the four site Hubbard model, CGA exhibits pinned occupancies at the weak correlation regime, and therefore is not able to give an accurate description of the model. All of the NOFAs shown tend to unitary ONs in the strong correlation limit. The three-fold degeneracy can be explained by considering that PNOFi {i = 5, 7} functionals lead generally to localization of the molecular orbitals in the NO representation. Nevertheless, there is an equivalent canonical representation that can afford delocalized molecular orbitals adapted to the symmetry of the system upon diagonalization of the matrix of Lagrange multipliers (λ in equation (8)) obtained from optimized NOs [56]. Thus, the NOs obtained by using the J KL-only approximations, which are plotted on the right side of figure 7, transform into the symmetry-adapted orbitals shown on the left side of figure 7, so PNOFi {i = 5, 7} functionals are able to retrieve the orbitals adapted to the symmetry of the system in the canonical orbital representation [56]. The three-fold degeneracy is only a matter of the nature of PNOFi {i = 5, 7} functionals, but does not introduce any artifacts in the description of properties.

Inhomogeneous Hubbard model
MBB, PNOF5 and PNOF7 reproduce the exact results for the energy, double occupancy of sites, and NO occupancies in the homogeneous two site Hubbard model. Nevertheless, the MBB functional does not really recover the exact expression of Löwdin-Shull [39,[57][58][59] for any two-electron system (apart from some phases [58]), but PNOFi (i = 5 and 7) functionals do. In order to prove the deviation of the MBB functional, and that PNOFi (i = 5 and 7) actually recover the exact two-electron functional, we include two additional tests for these approximations. Results obtained by using the CGA approximation are also included, to test whether the accuracy of the method still holds. The tests proposed use an inhomogeneous two site Hubbard model, so a nonzero on-site energy is set such that it gives rise to a potential difference ∆v = v SB − v SA = 2v. Since the system is nevermore symmetric, the orbitals adapted to the symmetry are no longer the NOs, so the latter may arise from optimization of the energy for each NOFA. This model is also used in benchmarking [24], e.g. for the widely used Bethe ansatz local density approximation (BALDA) developed by Capelle and collaborators [60][61][62][63][64] (which fails to recover the exact results). Then, we have compared the exact results with the results obtained by using these NOFAs for the total energy and for the difference between site occupancies (i.e. |∆n| = |n SB − n SA | where n SX is the occupancy of site X).
In figure 8 we have plotted for a fixed value of t = 0.5 and four values of U the exact energies and the energies obtained with MBB, CGA, PNOF5 and PNOF7. In general, the ∆v = 0 limit recovers the homogeneous two site Hubbard model where all functionals produce exact energies for any U, as we observe in figure 8. For U = 0 the exact energy is given by E = − (2t) 2 + (∆v) 2 , where the RHF energy is exact since there are no electron-electron interactions for any ∆v; hence MBB, CGA, PNOF5 and PNOF7 reproduce the exact Hartree-Fock energy. Nevertheless, once the on-site interaction is turned on, MBB and CGA no longer produce exact results for the two site inhomogeneous Hubbard model, and only PNOF5 and PNOF7 yield exact energies.
Despite that the energy obtained with MBB and CGA functionals for U = 1 (top right plot in figure 8) seems still to be exact, a property such as the |∆n| reveals that MBB and CGA results are no longer exact (see figure 9). In figure 9 we have plotted |∆n| for t = 0.5 and the same four values of U that we used in figure 8. We observe from the |∆n| test that MBB and CGA only yield exact results in the non-interaction limit (U = 0); even in the weak correlation region (U = 1) MBB and CGA |∆n| values differ from the exact ones. A remarkable result is that PNOF5 and PNOF7 have proven to retrieve the exact functional of Löwdin-Shull. Thus, energies and properties, like for example |∆n| versus ∆v or the charge-transfer to Mott-insulator transitions (which happens around ∆v ≈ U), obtained with PNOF5 and PNOF7 in any correlation regime for the two site inhomogeneous Hubbard model will be exact.
The on-site potential can be more interestingly modulated if a system with more sites is employed. Hence, we have used an oscillatory potential such as that introduced     in the Aubry-André model [49] to carry out calculations on the ten site inhomogeneous Hubbard model. Note that the Hamiltonian of the system is now given by equation (10) with α = 1/10, V = 2.0 and δ = −2π/10. According to the results given in figure 10, and similarly to the results obtained for the inhomogeneous Hubbard dimer, CGA and MBB approximations fail dramatically when studying inhomogeneous systems, except in the weak correlation limit. PNOF5 and PNOF7 approx imations not only show good asymptotic behavior, but also lie close to the exact curve for any correlation regime. Therefore, systems with spatial inhomogeneities are also well-described by PNOF approximations.

Conclusion
The resurgence of interest in reduced density-matrix functionals has brought a new class of theoretical methods that have appeared recently in the literature [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. These approximations are given in terms of NOs and ONs, so they should be referred to as NOFAs. A distinction must be drawn between NOFAs developed with simple assumptions about the hole function, e.g. MBB, Power, CA, or CGA approximations; and NOFAs that are built imposing necessary known N-representability conditions of the 2RDM, such as PNOF5 or PNOF7. In this paper the robustness of a set of NOFAs is put into test by using the homogeneous Hubbard model with 2, 4, 6, 10, 12, and 14 sites, and the inhomogeneous model which includes nonzero on-site potential for the Hubbard model with two and ten sites. A simple comparison between the exact FCI calculations and approximated results for varying interaction strengths U/t is carried out.
Both PNOF5 and PNOF7 reproduce the exact Löwdin-Shull functional for the two site system, while the other NOFAs fail for one or another test. For all test cases studied in the homogeneous Hubbard model, PNOF7 results are consistently in outstanding agreement with FCI. Among J K-only NOFAs only CGA is able to provide energies close to the exact ones, but fails significantly for ONs of the NOs, and yield a discontinuous curve for the double occupancy of the sites. The enhanced accuracy of CGA is in fact due to the inclusion of a proper particle-hole symmetry, which is not the case for some other J K-only NOFAs. Nevertheless, CGA violates many fundamental properties, such as N-representability conditions or antisymmetry, of the 2RDM [42]. Thus, CGA fails dramatically for the inhomogeneous Hubbard model, i.e when the onsite potential varies. PNOF7 presents particle-hole symmetry at the same time that the fundamental properties of the 2RDM are conserved, and thereby the inhomogeneous Hubbard is also well described.
In view of the results obtained for the Hubbard model with varying interaction strength, PNOF7 can describe not only weakly correlated systems, but also problems where strong correlation effects arise. This opens a new avenue where PNOF7 becomes a tool for studying many interesting applications which include confined fermions, disorder and critical behavior in optical lattices, effects of spatial inhomogeneity in strongly correlated systems, and various critical phenomena in 1D chains, among others. This is where-apart from the wellestablished variational wavefunction methods of Gutzwiller and Baeriswyl, different perturbative expansions, or quantum Monte Carlo methods-the widely used BALDA approximation has so far led, despite the pitfalls obtained with it even in the simple two site nonsymmetric Hubbard model [24].
To conclude, the present study proves the importance of developing functionals that satisfy at least the analytically necessary N-representability conditions of the 2RDM, in order to obtain consistent results in systems with either weak or strong electronic correlation. Besides, concerning the improvement of PNOF7 over the independent pair model PNOF5, it shows that a well-balanced inter-pair correlation is crucial to account properly for electron correlation effects. Our calculations shine a light on the insight, successes, and limits of current NOFT approaches, and may be useful for the development of new natural orbital functional approximations.