Capturing the electron-electron cusp with the coupling-constant averaged exchange–correlation hole: a case study for Hooke’s atoms

.


I. INTRODUCTION
Due to its relatively low computational scaling combined with high accuracy in the study of electronic structure of many-body systems, density functional theory (DFT) has become the most widely used electronic structure method with an increasing range of applications in condensed-matter physics, quantum chemistry, and materials science.In principle, DFT is an exact method with which the ground-state energy and electron density can be computed, from which many important physical and chemical properties can be extracted 1 .In practice, approximations must be introduced to DFT to make it computationally useful; in the Kohn-Sham formulation of DFT (KS-DFT), 2 the exchange-correlation (XC) component of the energy which carries the many-electron effects must be approximated.Therefore it is the quality of the XC approximation that determines the quality of a DFT calculation in predicting the total energy and other ground-state properties of interest.
An exact expression for the XC energy can be obtained in terms of the electron density n(r) and the coupling-constant (λ ) averaged XC hole density nxc (r, r ) via their Coulomb interaction 3 as where nxc (r, r ) is the probability depletion of finding an electron at r , given an electron located at r. nxc (r, r ) is entirely attributed to quantum effects, which include the selfinteraction correction, the Pauli exclusion principle (arising from the exchange symmetry of indistinguishable electrons), and the electron-electron correlation resulting from the Coulombic repulsion. 4The first two effects give rise to the exchange hole density n x (r, r ), which is completely negative and independent of the coupling constant.The remaining quantum effects produce the correlation hole, which is defined by subtracting the exchange hole density from nxc (r, r ) as nc (r, r ) = nxc (r, r ) − n x (r, r ), yielding the λ -averaged correlation hole.
Eq. ( 1) guarantees an accurate evaluation of XC energy if an accurate XC hole model is provided.Thus the quality of XC hole models underpin the XC energy approximation and play a fundamental role in understanding the quality and assessing the performance of a diverse range of density functional approximations (DFAs) when applied to different systems and properties.However, practical DFT calculations only require approximations of the XC energy, leading to a tendency to neglect the importance of XC holes in favor of directly modeling the XC energy.This trend has led to there being relatively few XC holes studies.Notably, early successful DFAs such as the PW91 approximation of Perdew and Wang 5 were based on modeling the XC hole, and the construction of the strongly constrained and appropriately normed (SCAN) density functional was also grounded in the understanding of XC holes. 6ecently, there have been new DFA developments based on XC holes. 7lthough being formally defined in Eq. ( 1), XC holes are challenging to evaluate accurately, contributing to the scarcity of the XC hole studies.There are two significant challenges associated with this: i) the XC hole has to be calculated for each coupling constant λ to evaluate the coupling-constant integrated XC hole; ii) high-level electronic structure methods are required to obtain accurate ground-state wave functions for each λ .These methods typically have high-rank polynomial scaling with system size and become computationally in-tractable for large systems.The Lieb optimization approach 8 can address challenge i) by transforming the problem of finding the ground-state electron density of a λ -interacting system into maximizing the Lieb functional of the λ -dependent external potential, 9 while keeping the electron density fixed.In combination with the coupled-cluster singles and doubles excitation method (CCSD), the Lieb optimization method has been applied to two-electron systems, such as the Helium isoelectronic series, with a focus on the XC energy. 10,11The CCSD method is exact in the complete basis set limit, equivalent to the full configuration interaction (FCI) approach for two-electron systems.
However, the λ -averaged XC hole has not been studied using the Lieb optimization with a CCSD reference wave function, even for the simple two-electron systems.Therefore, it is currently not known how the basis set influences the quality of the calculated XC hole and the associated electronelectron cusp condition 12,13 of the correlation hole 14 when the coupling-constant averaged quantities are considered.The electron-electron cusp condition describes the behavior of a many-electron wave function when two anti-parallel electrons come infinitesimally close to each other, arising due to the singularity of the Coulomb repulsion at the coalescence point.This dynamical correlation effect at zero separation introduces non-smoothness into the many-body wave function, which cannot be effectively represented by orbital product expansion wave functions 15 .Increasing the basis set size can help reduce the cusp error, but this approach is limited by the unfavorable computational scaling of high-level electronic structure methods.
In this study, we examine the electron-electron cusp condition and basis set effects on the XC hole through the calculation of the Lieb functional at the CCSD level for a simple model system, namely the Hooke atom (Hookium).By introducing the harmonic-oscillator potential as the external potential in the Hamiltonian of a two-electron system, given in atomic units as the resulting problem is one of the few examples of a twoelectron system for which a series of exact solutions exist, in this case an infinite set of solutions corresponding to different harmonic confinement constants, k. 16,17 The Hookium atom is therefore a useful reference for evaluating XC hole models since the exact XC holes can be computed.We commence in Section II by providing an overview of the theoretical framework for computing the λ -dependent XC hole, the Lieb optimization method, and the electron-electron cusp condition in Coulombic systems, and the solvable Hookium model.Computational details are then discussed in Section III.In Section IV we examine the basis set effects and cusp condition effects on the XC hole calculated at the CCSD level at λ = 1 (the physical system), for which the exact wave function solution is known.We then compare and benchmark the local density approximation (LDA) and Perdew-Burke-Ernzerhof (PBE) XC hole models with the coupling-constant averaged XC hole from Lieb optimizations at the CCSD level.
System-and angle-averaged XC holes are calculated to enable direct comparison between the benchmark data and these simple desnsity-functional models.Finally, we conclude our work with a brief summary in Section V.

II. THEORY AND METHODOLOGY
A. The exchange-correlation hole and the coupling constant In KS-DFT, the ground-state energy of a many-electron system in an external potential v ext (r) is obtained by mapping the interacting system of electrons to an auxiliary noninteracting system of electrons with the same electron density.The Schrödinger equation for this auxiliary system can then be solved in a basis of one-electron orbitals. 2The ground-state energy is thus expressed as a functional of the electron density n(r), which can be resolved into the sum of several contributions as where T s is the non-interacting kinetic energy, which is evaluated exactly using the KS orbitals, and E H the classical electrostatic Hartree energy, which is evaluated exactly in terms of n(r).The only term in Eq. ( 3) which must be approximated is the XC energy E xc [n], which describes all of the many-electron effects in the system.The KS non-interacting system may be linked to the physically-interacting system by continuously varying the strength of the electron-electron interaction between the noninteracting and physically-interacting limits by scaling the two-electron operator Vee by a coupling-constant λ between zero and one.The electronic state evolves through a family of solutions to the λ -interacting Hamiltonian, where T is the kinetic energy operator and v λ a modified external potential, thus estabilishing an adiabatic connection between the non-interacting and physically-interacting systems. 18The modified external potential v λ is determined for each interaction strength such that the density remains constant at the physical (λ = 1) density for all λ .Clearly, v λ reduces to the local KS potential v s when λ = 0 and is equal to the physical external potential v ext when λ = 1.Supposing Ψ λ is the normalized ground-state manyelectron wave function of the λ -interacting system with N e electrons, the second-order reduced density matrix is expressed as 13,19 This two-particle density may be used to evaluate the expectation value of two-body operators 4 , but it cannot be diagonalized by a unitary transformation of one-electron basis functions 13 .The XC hole density at each coupling strength λ is defined as, where the second term removes the classical Hartree contribution to the two-particle density n(r)n(r ), with the remaining n λ xc (r, r ) accounting for only the XC effects.The λ -averaged XC hole density is given by coupling constant integration over this quantity, nxc (r, r ) = from which an exact expression for E xc can be obtained, shown in Eq. ( 1).At λ = 0, the XC hole is reduces to the exchange hole, where ψ iσ (r) are the KS spin-orbitals.Therefore the λaveraged correlation hole can be defined by Furthermore, since the Coulomb operator has spherical symmetry, the XC energy may be computed exactly from the spherically-averaged XC hole.As a result the systemand spherically-averaged XC hole density nxc (u) is a useful quantity that can be modelled in order to construct XC energy functionals.This may be written in terms of the distance vector u = r − r as where Ω u is the solid angle around direction u and integration is carried out to average over this angle and the spatial coordinates of the entire system.It is a remarkable result that the XC energy may then be expressed precisely as one-dimensional integral over u = |r − r| for any system, where we identify ε xc (u) = 4πu nxc (u) .The exact systemand spherically-averaged exchange and correlation holes satisfy the following sum rules respectively,

B. The Lieb optimization
Given a Hamiltonian Ĥλ [v λ ], the ground-state energy E λ [v λ ] for an N e -electron system is given by the Rayleigh-Ritz variation principle as where W N e is the set of all L 2 -normalized, antisymmetric N eelectron wave functions with a finite kinetic energy.The ground-state energy in Eq. ( 15) is well-defined for all potentials v λ ∈ χ * with χ * = L 3 2 + L ∞ , a vector space containing all Coulomb potentials.For a variationally-determined solution to Eq. ( 15), E λ [v λ ] is concave and continuous in v λ .
Following the convex-conjugate formulation of DFT by Lieb, 9 the universal density functional F λ [n] may be defined as the Legendre-Fenchel transform to the ground-state energy of Eq. ( 15) as which is convex in n by construction and thus may be defined for arbitrary as defined in Eq. ( 16) yields an expression for the Hohenberg-Kohn variation principle in which the biconjugate functional ] and χ = L 3 ∩ L 1 is the dual vector space to χ * and which encompasses all N erepresentable densities.The conjugate functionals Eq. ( 16) and Eq. ( 17) are related by Fenchel's inequality as which becomes an equality by maximization of the right-hand side with respect to v λ which is the same, for non-degenerate solutions, as satisfying the stationary condition By definition, E * λ [v λ ] is concave in v λ hence has no more than one stationary point; if a solution to Eq. ( 19) exists, it is therefore unique.This can also be expressed by re-arrangement of Eq. ( 18) to the form , which becomes an equality by minimization of the right-hand side with respect to n(r) thus satisfying the stationary condition where v λ is the optimizing potential.In the Lieb optimization method, the universal density functional F λ is maximized with respect to the potential v λ (r) for a given electronic structure method with energy functional E λ and yielding density n(r).To construct the density-fixed adiabatic connection, the optimizing potential v λ (r) is that for which E λ yields the physically-interacting λ = 1 density for all values of λ ∈ [0, 1]. 10,11he universal density functional F λ may be written as a sum of terms according to the Kohn-Sham decomposition as 2 in which T s is the non-interacting kinetic energy, E H is the classical Coulomb energy, E x is the exchange energy and E c,λ is the λ -interacting correlation energy.Substituting Eq. ( 21) into Eq.( 20) yields an expression for the optimizing potential in terms of its individual contributions, Identifying that v λ =1 = v ext , the external potential due to the electrostatic potential of the nuclei, and v λ =0 = v s , the Kohn-Sham potential may be eliminated from Eq. ( 22) to yield the expression for the optimizing potential at interaction strength λ as In order to optimize F λ with respect to the potential, it is expanded in a Gaussian basis as proposed by Wu and Yang as 20,21 in which v H is the Coulomb potential evaluated with an input λ = 1 density n in , v ref is a reference exchange potential also evaluated on this density to ensure that v λ has the correct asymptotic behaviour and g t are a set of Gaussian functions with expansion coefficients b t .The form of the reference potential employed in this work is that of the a localized Hartree-Fock potential, 22 corrected at long-range by an approximate Fukui potential. 23The details of the construction of the reference potential are given in Appendix A.
With the parameterization of the potential in Eq. ( 24) the Lieb functional can be defined as an optimization of the objective function with respect to variations in the potential basis coefficients b; the gradient of Eq. ( 25) with respect to the potential basis coefficients is given by whilst the second derivative of the objective function with respect to the potential basis coefficients is given by It can be seen from Eq. ( 26) that the stationary condition of Eq. ( 19) will be satisfied where the iterating density n λ ,b becomes equal to the input density n in .In this work, the objective function is optimized by an approximate Newton approach implemented in the QUEST code; this is a secondorder optimization algorithm in which the Hessian is approximated by the non-interacting Hessian, given by Eq. ( 27) at λ = 0. 24 In this process, the potential basis coefficients are updated at each iteration using a backtracking line-search and the wave function E λ evaluated with the corresponding potential v λ ,b , yielding the energy and iterating density n λ ,b from which the objective function Eq. ( 25), gradient Eq. ( 26) and approximate Hessian are constructed.At the point of convergence, for which Eq. ( 26) becomes zero, the optimizing potential may be used to obtain the λ -interacting one-and twoparticle reduced density matrices required for the construction of the λ -interacting XC hole as described in Subsection II A.
With the above calculations completed for each λ , a series of λ -dependent and then λ -averaged quantities such as the XC holes and XC energies given in Eqs.(5−11) can be readily obtained.

C. The electron-electron cusp condition
For a Coulombic system, the electron-electron cusp condition describes the behavior the electrons in exact eigenfunctions of the Schrödinger equation, which exhibit a cusp at the points of electron coalescence due to singularities in the Coulomb potential at such points. 25Specifically, the first derivative of the wave function is discontinuous at these points.The electron-electron cusp condition may be expressed using the pair-correlation function, defined as the ratio of the two-particle density to the product of the one-particle densities 14 Given the spherically-averaged pair-correlation functional defined analogously to the spherically-averaged XC hole as g(r, u) = dΩ u 4π g(r, r + u), the electron-electron cusp condition is written as 12,13 ∂ g(r, u) This may be written in terms of the system and sphericallyaveraged XC hole defined in Subsection.II A, using the relation between g(r, r ) and n 2 (r, r ), as (30) Due to the Pauli principle, the cusp condition only arises between electrons with anti-parallel spin and is thus exclusively a correlation effect.The electronic cusp condition can therefore be written in terms of the system and spherically-averaged correlation hole as

D. Hookium atoms
A Hookium atom is a model system comprising two electrons confined by a harmonic potential rather than a Coulomb potential, 16 with electronic Hamiltonian given in Eq. ( 2).Introducing the center of mass coordinate R = (r 1 + r 2 )/2 and the relative separation vector u = r 1 − r 2 , the Hookium atom Hamiltonian may be resolved into a center of mass and relative separation term as The second term Ĥ(u) is of particular interest as it describes the relative motion between the two interacting electrons bound by the harmonic potential and is effectively a one-body problem with Schrödinger equation Ĥ(u)ϕ(u) = εϕ(u).Using a separation of variables to write ϕ(u) in terms of the product of radial and angular components where Y lm is the spherical harmonic function describing the angular wave function, a second-order differential equation for T (u) can be obtained.By inserting the regular solution into the differential equation, a recurrence relation 16 can be found for the coefficients {a i , i ≥ 2} where we only consider the ground state with the angular momentum l = 0.A series of exact solutions can be determined by imposing the condition a N = a N+1 = 0 at i = N, leading to a i = 0 for all i ≥ N. Consequently, N represents the polynomial order of T (u) in the radial wave function ϕ(u).
For the ground state with l = 0, N is roughly proportional to k −7.9 as observed by fitting the values of N against k 16 .Since k is the harmonic constant which determines the strength with which electrons are confined, an increase in N implies less confinement and a more radially-diffuse electron density.
However, it is obvious that there doesn't exist an analytical wave function solution for the Hookium atom with the electron-electron interaction scaled by an arbitrary λ = 1.Therefore, the coupling-constant-averaged correlation hole for the Hookium atom has seldom been studied, and only the correlation hole at λ = 1 has been comprehensively studied 26 and used to benchmark correlation hole models [27][28][29] .
For example, the exchange hole and the correlation hole of the Hookium atom with k = 1/4 (corresponding to N = 2) for the λ = 1 case have been carefully studied in Ref. 26, which is also used to benchmark the system-and angle-averaged XC hole models of different meta-GGAs 27 .Using only the λ = 1 results, the validity of the electronic cusp condition in the ground state of the Hookium atom for arbitrary harmonic confinement k has been demonstrated 28 and it has been demonstrated that the LDA hole model can precisely capture the cusp condition of the Hookium atom 29 .
In this work, we employ the exact solutions of the Hookium atom at λ = 1 to benchmark those calculated from the CCSD wave function.We then use the Lieb optimization with a CCSD wave function to calculate λ -averaged XC holes for the Hookium atom, which can serve as a benchmark for XC hole models.

III. COMPUTATIONAL DETAILS
In this work all calculations are carried out using the QUEST code with a the spin-restricted CCSD wave function as the reference method.The convergence of self-consistent field calculations was accelerated using the C1-DIIS method, with a convergence threshold of 10 −12 a.u. on the norm of the DIIS error vector.For the CCSD calculations the convergence threshold for both the excitation amplitudes and deexcitation amplitudes, required for evaluation of the CCSD one-and two-particle densities, was 10 −10 a.u.for the norm of the difference of the amplitudes between iterations.
Lieb optimizations were carried out with at the CCSD level for a range of interaction strengths λ ∈ [0, 1] using the approximate Netwon method described in Subsection II B. In each case the CCSD λ = 1 density was used as input to the Lieb functional, in order to fix the density along the adibatic connection at its physical value.The potential was regularised using the smoothing norm method with a regularization parameter of 10 −5 a.u.Convergence of the Lieb optimization was based on the norm of the gradient with respect to potential expansion coefficients, with a convergence threshold of 10 −8 a.u.used throughout.1][32] In each case, the basis sets were uncontracted spherical Gaussians with exponents for the Helium atom used throughout.To evaluate the system-and spherically-averaged XC holes, a nested numerical quadrature was employed.The spherically-averaged XC hole n λ xc was constructed by angular integration using an order-41 Lebedev quadrature grid at each reference point 33,34 , leading to, with quadrature weights w Ω i and nodes r Ω i associated to the angular quadrature nodes (ϕ i , θ i ) by The system-averaging was then carried-out by numerical integration of the reference point using a full quadrature grid, with angular component again given by the order-41 Lebedev quadrature and radial component constructed using the scheme of Lindh, Malmqvist and Gagliardi 35 with a relative error threshold of 10 −10 a.u., where w r i are the weights of the radial quadrature and r i j the product of radial quadrature nodes r r i and angular quadrature nodes r Ω j .Exact analytical results for the Hookium atom at λ = 1, as described in Subsection II D, were also calculated with Mathematica, allowing us to carefully assess the accuracy of the finite basis CCSD calculations.

A. Accuracy of finite-basis CCSD Hookium solutions
As described in Subsection.II B, n(r) given by CCSD is used as the reference electron density of the physical interacting system for the Lieb optimization in Eq. (25).Therefore, we first assess the quality of CCSD calculated total energies and densities with a range of orbital basis set sizes for Hookium by comparing them with the exact analytical results.

Total energies
Table I lists the percentage errors (PEs) for the finite-basis CCSD total energies with respect to the exact results for different solutions to the Hookium atom labelled by N, as described in Section II D, computed with different basis sets.All PEs are positive, as expected since CCSD is equivalent to FCI for these two-electron systems, and so the energy approaches the complete basis FCI energy from above.In general, the accuracy of the energies can be improved systematically by using basis sets with a higher cardinal number X or higher augmentation with diffuse functions Y .This leads to a reduction in the PEs to be in the range 0.2% -2.5%.It is clear that for the solutions with N < 5 PEs below 1% can be achieved with triply augmented basis sets with cardinal numbers of 4 or above.Indeed, adding extra diffuse functions does not further improve the accuracy of the results for these solutions.However, for larger values of N is it is essential to include many more diffuse functions to obtain reasonable accuracy.For 5 ≤ N ≤ 8 pentuple augmentation is required to achieve PEs below 2% and for N > 8 hextuple augmentation is required.The dependence on cardinal number X is less significant, once sufficient diffuse functions are included for a given value of N, there appears to be little benefit in using basis sets with X > 4.

Electron densities
In Figure 1, we plot the CCSD electron densities of Hookium atom solutions with 2 ≤ N ≤ 5 radially from the atomic nucleus.For comparison the densities of the corresponding exact solutions are also shown.The convergence of CCSD electron densities at each Hookium solution N is examined by gradually increasing the size of the t-aug-cc-pVXZ basis set by changing the cardinal number X from 2 to 6.
Figure 1 shows that, as the order of the Hookium solution N increases from 2 to 5, the corresponding electron density becomes increasingly spatially diffuse.Interestingly, the convergence behavior of the density with respect to basis set appears to be dependent on whether the value of N is even or odd.Specifically, for the N = 2 and N = 4 solutions, the CCSD densities converge to the corresponding exact densities relatively quickly with increasing basis set size, and there is no discernible difference in the results obtained with a basis sets with X ≥ 3.However, convergence of the density with respect to the basis set is considerably slower for the N = 3 and N = 5 solutions.It should be noted that, for N = 5, the CCSD energies and densities in the largest basis set t-aug-cc-pV6Z both have a greater error than those from the smallest basis set taug-cc-pVDZ considered here, as can be seen in Table I.This indicates that for N = 5 the triply augmented basis sets are not adequately diffuse and errors could be reduced by further augmentation of the basis set.Indeed, from the analysis of the electron density it is clear that a sufficiently diffuse basis set would be required to to represent the electron density accurately as N increases, consistent with the analysis of the CCSD total energies in Section IV A 1.

Orbital basis N
When comparing the CCSD λ = 1 XC holes with those of the exact solutions, the errors are dominated by incompleteness of the finite basis set in which the orbitals are expanded.However, when comparing the CCSD λ = 0 exchange and λ = 1 correlation holes with those of the exact solutions, an additional source of error is introduced; the incompleteness of the basis set in which the potential is expanded, shown in Eq. ( 24), and the associated numerical errors in the convergence of the Lieb optimization at λ = 0.

The exchange hole
For closed-shell two-electron systems, the exchange energy is related to the Hartree energy as and this dominates the XC energy.In addition, the exchange hole is related to the electron density for the closed-shell two-electron system as n x (r, r ) = −n(r ).As a result, the convergence of the exchange hole with respect to basis set size is the same as that observed for the electron density.This can be seen in the upper panels of Figure 2, Figure 1 and Table II.Furthermore, different convergence patterns are observed in Figure 2  Table III presents the PEs of exchange energy E x obtained using Lieb optimization at the CCSD level for solutions to the Hookium atom of different order N with increasing basis set size.An initial observation that can be made is that different error characteristics are again exhibited for solutions with odd and even values N respectively.Specifically, for N = 2, 4 the error is relatively small for all basis-sets with X > 2 while for solutions with N = 3 and N = 5, the PEs are generally larger in magnitude by comparison.Secondly, for a solution with any given N, once the number of diffuse basis functions is sufficient, increasing the cardinal number of the basis set will not increase the accuracy of the energy.For solutions with N = 2, 3, 4, 5, the accuracy does not significantly improve beyond Y = d, t, q, q respectively.For N ≥ 6, the improvements in accuracy with increasing cardinal number ceases at Y =p.

The correlation hole and the description of the cusp
We now consider the system-and spherically-averaged correlation holes for Hookium atoms with N = 2, 3, 4 in Figure 3 and compare them with those of the corresponding exact solutions.Figure 3 illustrates that, as the order of the Hookium solution N increases from 2 to 4, the exact correlation holes become increasingly shallow.This trend is consistent with the behavior observed in the electron densities plotted in Figure 1 and the exchange holes in Figure 2. In addition, the cusp at u = 0 becomes shallower as the order of the solution N increases, indicating that the cusp effect is less significant for more diffuse electron densities.
Figure 3 displays the effect of basis set size on the correlation holes obtained via Lieb optimization at the CCSD level.With the exception of the t-aug-cc-pV6Z basis set for the most diffuse solution with N = 4, enlarging the basis set results in an overall improved representation of the correlation holes with respect to those of the analytical solutions.In the lower panels of Figure 3, the error in the system-and sphericallyaveraged correlation holes are plotted radially from the atomic nuclei.Compared with higher-order solutions of larger N, the maximum error for N = 2 arises at u = 0, indicating that the cusp condition is more significant the more localized the electrons are, consistent with Eq. (31).For the Hookium solution with N = 4, the error is more uniformly distributed radially than for solutions with N = 2 and N = 3.
To quantitatively estimate the effect of the electron-electron cusp, we define a cusp-effect driven error δ E PE c in the correlation energy as (39) where a characteristic distance u c defining the electronelectron cusp region is determined by the solution to

N=4
FIG. 2. System-and spherically-averaged exchange holes n x (u) calculated with the Lieb functional at the CCSD level for the Hookium solutions with N = 2, 3, 4 (upper panels) and the corresponding deviations with respect to exact results (lower panels).The orbital basis sets t-aug-cc-pVXZ with X =D, T, Q, 5, 6 are employed for each N, with corresponding potential basis sets as described in Section III.Table IV presents the cusp-effect driven errors δ E PE c computed for Hookium solutions of N = 2, 3, 4 with increasing basis set sizes.The trends of PEs with respect to N are in agreement with the observations in Figure 3; solutions with higher N values exhibit smaller errors resulting from the cusp effect.This is also consistent with the cusp condition for the correlation hole expressed in Eq. (31).Since solutions of increasing N have an increasingly diffuse electron density, as shown in Figure 1, the integral dr n 2 (r)/(2N e ) decreases, and n c (0) becomes smaller 36 .This, in turn, leads to a reduction in n c (0) , resulting in a flatter n c (u) approaching u = 0 for solutions of larger N, as demonstrated in Figure 3.It follows therefore that the cusp-driven error becomes much less significant when N is large or the electron density is diffuse.
Table IV also shows that, as the cardinal number X of the basis set is increased from 3 to 6, the cusp errors become consistently smaller for solutions with N = 2, 3.However, the sit-TABLE III.Percentage errors of the exchange energy E x calculated by Lieb optimization at the CCSD level for the Hookium solutions with N = 2 − 11.Orbital basis sets of Y -aug-cc-pVXZ (X = D, T, Q, 5, 6; Y = d, t, q, p, s) and potential basis sets of (Y − 1)-aug-cc-pVXZ (X = D, T, Q, 5, 6) are employed, respectively.
Figure 4 plots the spherically-averaged correlation energy density with . It shows that the cuspeffect driven error that arises at short interelectronic separations is significantly attenuated at larger values of the interelectronic distance u.Conversely, the errors in the correlation hole become more pronounced at larger values of u in ε λ =1 c (u).Overall, calculations with the basis set with X = 6 yield a markedly improved description of ε λ =1 c (u) compared to the results obtained with smaller basis sizes for the less diffuse N = 2, 3 Hookium solutions.These observations validate the application of Lieb optimization at the CCSD level with the appropriate Gaussian basis sets to calculate XC holes accurately.
Table V collects the PEs of the correlation energy E λ =1 c calculated via Lieb optimization at the CCSD level, by subtracting the Lieb λ = 0 energy from the CCSD energy.Similar trends are observed for the correlation energy with respect to Hookium atom solution N and basis set size, compared with E x shown in Table III.While the greatest accuracy is obtained for Hookium atom solutions with N ≤ 4, the PEs for the correlation energy are usually 3-7 times larger than those for E x in the same calculation.For N = 2, 3, 4, the Y -aug-cc-pV6Z basis sets (Y =t, q, p, s) consistently yield accurate correlation energies with PEs of 1-2%.However, as N increases beyond 4, the improvement in accuracy achieved by increasing the basis set size (either via larger cardinal numbers X or increased augmentation Y ) is not as significant as was observed for E x .

The exchange and correlation hole
Figure 5 shows the XC holes n λ =1 xc (u) for Hookium atom solutions with N = 2, 3, 4. Convergence of the XC holes with respect to the basis set is smoothly achieved for all three cases.In the case of N = 2, 3, the correlation component of the XC holes have significant errors due to the electron-electron cusp, and increasing the cardinal number of the basis set leads to improved accuracy.It is worth noting that for the N = 2 solution, the correlation hole around u = 0 from Lieb optimization at the CCSD level are too shallow, which is compensated by the X hole being too deep, with the resulting error cancellation yielding a better accuracy for the XC hole than for either component individually.Overall, Lieb optimization at the CCSD level in the largest basis set t-aug-cc-pV6Z provides a satisfactory description of the XC hole for both the N = 2 and N = 3 Hookium atom solutions, for which the electron densities are relatively localized and not too diffuse.

C. Coupling-constant averaged XC holes and hole models
To assess the quality of the λ -averaged DFT XC holes, we first analyze how closely they align with the Lieb optimization results at the CCSD level calculated in the same basis set.This is important because CCSD-based Lieb optimizations have a much higher computational cost than DFT calculations and hence have a much greater limitation in terms of system and basis set size.
In Figure 6, we have plotted the accuracy of LDA hole densities evaluated by n LDA .The cardinal number X of the basis set was increased continuously from 2 to 6 to examine the convergence of this error with respect to basis set size.It is important to note that the errors corresponding to the X = 2 basis set show a different pattern of errors to those of the larger basis sets, indicating this small basis is typically insufficient to accurately evaluate the λ -averaged XC holes.This is consistent with the previous discussion concerning the λ = 1 case.
Figure 6 shows that the largest change in n LDA x (u) − n CCSD x (u) with respect to the basis set size occurs at the value of u with the largest error, particularly in the case of the N = 3 Hookium atom solution.Nevertheless, overall the LDA hole model calculations of the exchange holes exhibit a rapid convergence in their errors with respect to basis set size, suggesting that the accuracy of the DFT exchange holes is relatively insensitive to basis set size.
However, for the λ -averaged correlation holes, the situation is somewhat different.Changes in the error of the LDA correlation holes with respect to increasing basis set size are most visible around the first peak in the radial plots of the correlation holes, extending to subsequent peaks further from the nucleus for the N = 3, 4 Hookium atom solutions.As the basis set grows, the position of the first peak tends to shift toward u = 0.This could be attributed to the fact that the cusp-driven deviations become less significant with larger basis sets with CCSD-based calculations, whereas LDA is designed to satisfy the cusp condition and converges much more rapidly with basis set size.
Although the exchange hole is the dominant component of the XC hole, the errors of the LDA exchange holes and correlation holes are comparable in size.In Figure 6, the errors at short-range for the LDA XC holes are dominated by the correlation hole contribution, while errors at mid-range mainly arise from LDA exchange hole or both.Overall, accurate λaveraged XC hole (or correlation hole) calculations require the use of basis-sets with X ≥ 4. It is worth mentioning that similar trends are observed for the PBE hole model, the results TABLE V. Absolute percentage errors of the correlation energy E λ =1 c using Lieb optimization at the CCSD level for Hookium solutions with N = 2 − 11.Orbital basis sets of Y -aug-cc-pVXZ (X =D, T, Q, 5, 6; Y =d, t, q, p, s) and potential basis sets of (Y − 1)-aug-cc-pVXZ (X =D, T, Q, 5, 6) are employed.of which are presented in the Supporting Information.

Orbital basis N
The XC holes obtained with different basis sets are used to calculate the corresponding LDA XC energies, and Table VI presents the PEs of the LDA XC energies relative to Lieb optimization values at the CCSD level for basis sets of increasing size.Table VI shows that the PE variations of the LDA exchange energy are relatively small, within 0.2%, while the convergence behavior of correlation energy for the Hookium solutions with N = 3, 4 exhibits no clear trend.However, for the N = 2 solution, in which the cusp-effect driven error is the most significant in the CCSD calculations, the PE changes of the LDA correlation energy with increasing basis set size are slightly larger.The changes in the PE of E xc with increasing basis-set size are similar to that of the exchange energy, with only a change of 0.3% for the N = 2 solution from the smallest to the largest basis set; this is because E x represents the vast majority of E xc .
In Figure 7, the exchange holes, λ -averaged correlation holes, and XC holes from the LDA hole model, PBE hole model, and Lieb optimizations at the CCSD level are presented with the t-aug-cc-pV6Z basis set. Figure 7 shows that both the LDA and PBE hole models, in particular the LDA one, tend to localize the exchange hole, regardless of the order of the solution N. For the correlation holes, although both LDA and PBE hole models capture the cusp condition, they exhibit an almost linear behaviour before reaching their maximum value, resulting in an overly shallow correlation hole density in the small u region but an overly deep correlation hole at the intermediate u region.Figure 7 also shows that, for both the exchange and correlation holes, the PBE hole model is superior to the corresponding LDA hole model.However, the LDA and PBE model XC holes appear much more similar due to error cancellation between their respective X and C holes.

V. CONCLUSION
In this study, we have employed the Lieb optimization approach with CCSD used as the reference wave function method to obtain accurate representations of the XC hole of the Hookium atom -a model system for which exact solutions    can be obtained.Our investigation focuses on the difficulty in representing the electron-electron cusp condition within a finite Gaussian basis set, the manifestation of this in the correlation hole and effect on the cusp-related error of increasing the basis set size.We have found that the error resulting from the cusp effect can be effectively and sufficiently reduced by using a larger basis set, and that the cusp condition in the correlation hole becomes less significant for larger N Hookium atom solutions with diffuse electron densities.For smaller N Hookium solutions with electron densities that are more localized, the coupling-constant-averaged XC holes were calcu-lated using the Lieb optimization with CCSD reference wave function and used as a reference to benchmark DFT XC hole models.We confirmed the presence of significant error cancellation between the exchange hole and correlation hole for both PBE and LDA hole models and this results in their XC holes having a greater accuracy than either the exchange or correlation holes alone.
for solutions of even and odd values of N. Plots of the deviation of the finite-basis exchange holes from the exact solutions in Figure 2 indicate that basis-set convergence is generally reached with the t-aug-cc-pV6Z basis set for solutions of N = 2, 3, 4.

FIG. 3 .
FIG.3.System-and spherically-averaged correlation hole densities n λ =1 c (u) calculated by Lieb optimization at the CCSD level (upper panels) for Hookium solutions with N = 2, 3, 4 and corresponding errors with respect to exact results (lower panels).

FIG. 7 .
FIG.7.System-and spherical-averaged n x (u) , and λ -averaged nc (u) , nxc (u) calculated by LDA and PBE hole models as compared with the Lieb optimization results at the CCSD level for Hookium solutions with N = 2, 3, 4 in the t-aug-cc-pV6Z orbital basis set and d-aug-cc-pV6Z potential basis set.

TABLE I .
Percentage errors of the CCSD total energy E CCSD tot
B. Accuracy of finite-basis CCSD XC holes at λ = 1

TABLE VI .
Percentage errors of energy components E LDA CCSD+Lieb results under the same orbital basis set for Hookium solutions of N = 2, 3, 4.