Calculating excitation energies by extrapolation along adiabatic connections

In this paper, an alternative method to range-separated linear-response time-dependent density-functional theory and perturbation theory is proposed to improve the estimation of the energies of a physical system from the energies of a partially interacting system. Starting from the analysis of the Taylor expansion of the energies of the partially interacting system around the physical system, we use an extrapolation scheme to improve the estimation of the energies of the physical system at an intermediate point of the range-separated or linear adiabatic connection where either the electron--electron interaction is scaled or only the long-range part of the Coulomb interaction is included. The extrapolation scheme is first applied to the range-separated energies of the helium and beryllium atoms and of the hydrogen molecule at its equilibrium and stretched geometries. It improves significantly the convergence rate of the energies toward their exact limit with respect to the range-separation parameter. The range-separated extrapolation scheme is compared with a similar approach for the linear adiabatic connection, highlighting the relative strengths and weaknesses of each approach.


I. INTRODUCTION
The calculation of excitation energies in densityfunctional theory (DFT) is nowadays mostly done by means of linear response theory in the timedependent framework.
Linear-response timedependent DFT (TDDFT) [1] exhibits an excellent price-performance ratio and is, within the usual adiabatic semi-local approximations, very successful at describing excitations to low-lying valence states. However, these approximations introduce several limitations, especially for the treatment of static correlation [2], Rydberg and charge-transfer excitations [3,4], and double or multiple excitations [5].
Range-separated DFT constitutes an alternative to standard KS DFT [6,26] where the physical electronic Hamiltonian is replaced not by an effective non-interacting Hamiltonian but by a partially interacting Hamiltonian that incorporates the long-range part only of the electronelectron interaction [27][28][29][30][31]. This partially interacting Hamiltonian corresponds to an intermediate point along a range-separated adiabatic connection, which links the KS Hamiltonian to the physical Hamiltonian by progressively switching on the long-range part of the two-electron interaction, whilst simultaneously modifying the oneelectron potential so as to maintain a constant ground-state density.
In range-separated time-dependent DFT, the excitation energies of the long-range interacting Hamiltonian act as starting approximations to the excitation energies of the physical system and are corrected using a short-range density-functional kernel, in the same manner as the KS excitation energies act as starting approximations in linearresponse TDDFT. Several such range-separated linear-response schemes have been developed, in which the short-range part is described by an approximate adiabatic semi-local density-functional kernel and the long-range linear-response part is treated at the Hartree-Fock [32][33][34][35], multiconfiguration self-consistent field (MCSCF) [34,35], second-order polarization-propagator approximation (SOPPA) [35], or density-matrix functional theory (DMFT) [36] level.
Within the time-independent framework, a standard method for improving upon the excitation energies of the partially interacting Hamiltonian would be to use perturbation theory. However, given that perturbation theory in its standard Rayleigh-Schrödinger based formulation does not keep the ground-state density constant at each order in the perturbation, it has not led to a systematic improvement [37].
In this work, we propose a time-independent alternative method for correcting the excitation energies of the partially-interacting system, based on extrapolation along the range-separated adiabatic connection. Given that the long-range part of the interaction is included in the partially interacting system, its excitation energies constitute better approximations to the energies of the physical system than do the excitation energies of the KS system. The analysis of the Taylor expansion of the energies in the range-separation parameter µ about the physical system (µ → +∞) presented in Ref. 38 shows that the energies of the partially interacting system converge towards their physical limits as µ −2 . Using this information, it is possible to develop a scheme for extrapolating the energies of the physical system from the energies of the partially interacting system by following the ideas of Refs. 39, 40. This extrapolation scheme involves low-order derivatives of the energies with respect to µ and constitutes an alternative to perturbation theory and to range-separated TDDFT [32,41,42].
The extrapolation scheme is also applied to the linear adiabatic connection, where the interaction is scaled by a parameter λ going from 0 to 1, and where the analysis of the excitation energies around λ = 1 provides the required information to improve the estimation of the energies of the physical system from an intermediate point of the connection.
The expression for the energies of the partially interacting system and for their extrapolations are given in Section II for the range-separated and linear adiabatic connections. The extrapolation is subsequently applied to the range-separated energies of the helium and beryllium atoms and of the hydrogen molecule at its equilibrium and stretched geometries; for helium, we also use the linear adiabatic connection. The computational details are given in Section III and the results are discussed in Section IV.

A. Range-separated adiabatic connection
Range-separated DFT uses a partially interacting system, where the long-range part of the Coulomb interaction is included instead of the more traditional non-interacting KS system-see, for example, Ref. 31. In terms of the long-range (lr) electron-electron interaction operator W lr,µ ee = 1 2 w lr,µ ee (r 12 )n 2 (r 1 , r 2 )dr 1 dr 2 , (1) wheren 2 (r 1 , r 2 ) is the pair-density operator and w lr,µ ee (r 12 ) is the error-function interaction the Hamiltonian of the partially interacting system is given bŷ H lr,µ =T +V ne +Ŵ lr,µ ee +V sr,µ Hxc .
The parameter µ controls the range of the separation, with 1/µ acting as a smooth cut-off radius. This Hamiltonian also contains the shortrange Hartree-exchange-correlation potential op-eratorV sr,µ Hxc , whose role is to ensure that the ground-state density of the partially-interacting system is equal to the ground-state density of the physical system for all µ. Here Ψ µ 0 is the ground-state wave function of the partially interacting Hamiltonian andn(r) is the density operator. The remaining terms in the Hamiltonian of Eq. (3) are the usual kinetic-energy operatorT and nuclear-electron interaction operatorV ne = v ne (r)n(r)dr.
The eigenvectors and eigenvalues ofĤ lr,µ are the ground-and excited-state wave functions |Ψ µ k and energies E µ k of the partially interacting system These excited-state wave functions and energies provide natural first approximations to the excited-state wave functions and energies of the physical system. For µ = 0, they reduce to the single-determinant eigenstates and associated energies of the non-interacting KS Hamiltonian, while, for µ → ∞, they reduce to the excitedstate wave functions and energies of the physical HamiltonianĤ In Ref. 38, it was shown that the asymptotic expansion of the total energy of state k around the physical system is are the corrections entering at the second and third powers of 1/µ, respectively. Following the scheme proposed in Refs. 39, 40, it is possible to estimate the energy of the physical system E k from the energy of the partially interacting system E k and its first-and second-order derivatives with respect to µ.
¿From the Taylor expansion of the energies when µ → ∞, the first-order derivatives of the energies with respect to µ behave as around the real system. Inserting this into Eq. (8), the exact energies E k can be written as a function of the energies along the adiabatic connection and of their first-order derivative as This scheme gives extrapolated energies that are correct up to and including the second power of 1/µ relative to the energies of the physical system. The correction given by the extrapolation scheme vanishes at µ = 0 by construction, but should improve the description of the energies as soon as the interaction is switched on. One should note that the absence of a correction at µ = 0 is only due to the choice of 1/µ k as the basis for the expansion. Other basis functions such as µ 2 /(a + µ 5 ) would lead to a correction at µ = 0 but are not considered in this work. A more elaborate scheme can be developed by using also the correction E (−3) k and the secondorder derivative. In this case, the first-and secondorder derivatives are given by and, after eliminating E Higher-order derivatives should further reduce errors. Additionally, several points along the adiabatic connection could be used to perform the extrapolation to increase the accuracy of the extrapolated energies. However, only first-and secondorder corrections at a single point of the adiabatic connection are considered hereinafter.

B. Linear adiabatic connection
If the linear adiabatic connection is performed, then the partially interacting Hamiltonian is defined asĤ λ =T + λŴ ee +V λ whereV λ is adjusted to keep the ground-state density constant. This potential can be expressed in terms of the connecting parameter λ aŝ whereV λ c enters at second order in λ and is equal toV c at λ = 1. The energies of the partially interacting system can then be expanded around the physical system as where E (1) k and E (2) k are the contributions entering at the first and second power of (λ − 1), respectively. As in the range-separated case, by differentiation with respect to λ, it is then possible to extrapolate the energies of the physical system at first order by considering only the correction E When λ = 0, this extrapolation is equivalent to the first-order correction of Görling-Levy perturbation theory [23,25]. A second-order correction can be obtained by using also the correction E (2) k . The first-and second-order derivatives are and the extrapolated energies become 3

III. COMPUTATIONAL DETAILS
Calculations were performed for the He and Be atoms and the H 2 molecule with a development version of the DALTON program [43], see Refs. [44][45][46]. Following the procedure of Ref. 38, a full CI (FCI) calculation was first carried out to get the exact ground-state density within the basis set considered: uncontracted t-aug-cc-pV5Z for He, uncontracted d-aug-cc-pVDZ for Be, and uncontracted d-aug-cc-pVTZ for H 2 . A Lieb optimization of the short-range potential v sr,µ (r) was then performed to reproduce the FCI density with the long-range electron-electron interaction w lr,µ ee (r 12 ). Finally, an FCI calculation was carried out with the partially-interacting Hamiltonian constructed from w lr,µ ee (r 12 ) and v sr,µ (r) to obtain the zerothorder energies and wave functions.
Starting from the analytical form of the fit given in the supplementary material of Ref. 38, it is then straightforward to calculate the analytical derivatives of the energies with respect to µ. In the linear case, a cubic fit of the energies was performed. The extrapolated energies were calculated using Eqs. (11), (14), (17), and (20).
All the unextrapolated curves shown hereinafter correspond to the curves of Ref. 38.

IV. RESULTS AND DISCUSSION
A. Range-separated adiabatic connection of the helium atom

Ground-state energy
The results of the first-and second-order extrapolation schemes on the ground-state total energy of the helium atom are shown in Figure 1 (top). By construction, the extrapolation has no effect at µ = 0 and the ground-state energy of the KS system is, therefore, unaffected by the extrapolation. However, for µ > 0, the extrapolated energies show a systematic improvement with respect to the unextrapolated ground-state energy, that is, the ground-state energy of the partially interacting Hamiltonian without any correction.
Without extrapolation, a range-separation parameter of about 6 bohr −1 is needed to give an error smaller than 10 mhartree relative to the energy of the physical system. With the first-and second-order corrections added, the same accuracy is achieved with a range-separation parameter of only 2.8 and 1.5 bohr −1 , respectively. µ (bohr −1 ) Excitation energies (hartree) Figure 1: (Color online) Helium ground-state energy (top) E µ 0 , and first S (middle) and P (bottom) excitation energies ∆E µ k = E µ k − E µ 0 calculated without extrapolation (full lines), with first-order extrapolation (dashed) and second-order extrapolation (dot-dashed) as a function of µ. The dotted horizontal lines are the physical energies. The colored regions represent errors of ±10 and ±1 mhartree for the ground-state and excitation energies, respectively. Figure 1 also shows the effects of extrapolation on the lowest Rydberg S (middle) and P (bottom) excitation energies of helium. The convergence of the excitation energies towards their physical limit is overall improved by the first-and secondorder corrections with respect to the unextrapolated curves. In fact, the range-separation parameter required to achieve an accuracy of 1 mhartree is divided by approximately a factor 2 or a factor 3 by the first-and second-order schemes, respectively. For the excitation energies considered here, a range-separation value of 2 and 1 bohr −1 suffices to reduce the error to less than 1 mhartree with the first-and second-order schemes, respectively. The 1 S and 3 S excitation energies change monotonically with increasing µ. Accordingly, extrapolation provides a systematic improvement, the sign of the derivative pulling the excitation energies towards their physical limits at both first-and second-order levels.

Rydberg excitation energies
The 3 P excitation energy also changes monotonically with µ and the first-order extrapolation provides therefore a systematic improvement. However, the first-order extrapolated energy does not converge monotonically towards its physical limit (not visible), leading to a slight degradation of the excitation energies around µ = 1.5 bohr −1 at second order.
Finally, the 1 P excitation energy shows a nonmonotonic behavior even before extrapolation, exhibiting a "bump" for small µ that is probably a basis-set effect [38]. In fact, for such small µ, only the very long-range part of the interaction is modified, which is poorly described by Gaussian basis functions. The higher a given state is in energy, the more sensitive it becomes to this basis-set defect.
As a consequence, the 1 P excitation energy approaches its physical limit from above, its firstorder derivative changing sign around 0.7 bohr −1 . In this region, the extrapolated energies become less accurate than the unextrapolated energy. However, this behavior is observed only in a small region. As soon as the excitation energy recovers a monotonic convergence towards its physical limit (for µ larger than 0.7 bohr −1 ), the energy is improved by the extrapolation and converges faster to its physical limit.

B. Range-separated adiabatic connection for the valence excitation of the beryllium atom
The ground-state energy of the beryllium atom is shown in Figure 2 (top). Since beryllium has a core orbital, the convergence of its total energies µ (bohr −1 ) Excitation energies (hartree) Figure 2: (Color online) Ground-state energy E µ 0 (top) and excitation energies ∆E µ k = E µ k − E µ 0 (bottom) of beryllium as a function of µ. The unextrapolated energies are shown as full lines, the first-order extrapolated energies are plotted in dashed lines and the secondorder ones in dot-dashed lines. The energies of the physical system are given as horizontal dotted lines. An error of ±50 mhartree is colored around the physical ground-state energy and an error of ±2 mhartree is colored around the physical excitation energies.
is slower than for helium as the density is more contracted and a larger range-separation parameter is needed to describe correctly the core region. However, this affects all valence states in a similar fashion (not shown here).
Extrapolation systematically improves the convergence of the ground-state energy along the adiabatic connection. First-order extrapolation reduces the error to less than 50 mhartree with µ ≈ 5 bohr −1 , an order of magnitude smaller than the error without the extrapolation correction but still large. Second-order extrapolation gives the same error reduction already with µ ≈ 3 bohr −1 .
The effect of the extrapolation on the valence excitation energies of beryllium is shown in Figure 2 (bottom). As the errors associated with the core largely cancel in the excitation energies, the unextrapolated excitation energies already converge faster than do the total energies. With first-order extrapolation, an error smaller than 2 mhartree is reached with µ ≈ 0.5 bohr −1 , to be compared with a much larger error of 4 hartree in the total energies with the same µ value. The second-order extrapolation allows one to reach the same accuracy with a range-separation parameter as small as 0.3 bohr −1 . However, once again, some bumps are observed in the extrapolated energies probably due to the limited size of the basis set. This fast convergence of the excitation energies with respect to the range-separation parameter is due to the fact that in beryllium, static correlation is important and the multi-configurational character of the wave function is quickly established when the interaction is switched on; see Ref. 47.

C. Range-separated adiabatic connection for the hydrogen molecule
Finally, we consider extrapolation of the lowest excitation energies of the hydrogen molecule along the range-separated adiabatic connection, at the equilibrium geometry and at a stretched geometry. The results of the first-and second-order extrapolations on the singlet and triplet Σ + g → Σ + u excitation energies at the equilibrium geometry are shown in Figure 3 (top). First-and second-order extrapolations provide a systematic improvement in the excitation energies, µ ≈ 2 bohr −1 for first order and µ ≈ 1 bohr −1 for second order being sufficient to reproduce the physical energies to within 1 mhartree.
Having stretched the hydrogen molecule to three times the equilibrium distance, we apply extrapolation to the singlet and triplet excitations to the 1Σ + u state and to the double excitation to the 2Σ + g state-see the bottom part of Figure 3. Again, the improvement is systematic. The triplet extrapolated energy shows a monotonic behavior with respect to µ, whereas the singlet energy shows a slight bump at 0.8 bohr −1 . However, all extrapolated excitation energies converge faster than their unextrapolated counterparts. Extrapolation works remarkably well, reducing errors to less than 5 mhartree with µ ≈ 0.6 bohr −1 , compared with 2 bohr −1 without extrapolation. In particular, extrapolation allows us to describe double and single excitation energies equally well. In this case, one should note that the second-order scheme does not improve significantly the convergence of the µ (bohr −1 ) Excitation energies (hartree) 3Req Figure 3: (Color online) Unextrapolated (full lines), first-order extrapolated (dashed lines) and secondorder extrapolated (dot-dashed lines) excitation energies ∆E µ k = E µ k − E µ 0 of the H2 molecule at the equilibrium internuclear distance Req (top) and three times the equilibrium distance (bottom) as a function of µ. The excitation energies of the physical system ∆E k = ∆E µ→∞ k are plotted as horizontal dotted lines. An error of ±1 mhartree is colored around the physical excitation energies at equilibrium and an error of ±5 mhartree at stretched geometry. 1 1 Σ + g → 2 3 Σ + g (σ + u ) 2 and 1 1 Σ + g → 1 1 Σ + u excitation energies because of their nonmonotonicity probably due to the limited basis set.

Total energies
The total energies of the helium atom along the linear adiabatic connection are plotted in Figure 4 (top). When λ = 0, no interaction is included, so the KS energies are recovered as for µ = 0. When λ = 1, the full interaction is present and the energies of the physical system are recovered, which corresponds to the limit µ → ∞. The two limiting cases are, therefore, identical for the two adiabatic connections but the way they are connected differs.
The evolution of the total energies with respect to λ is almost linear. Although this behavior is easier to predict and should provide an efficient framework for extrapolations, the value of λ required to have an error of 10 mhartree is very close to 1; in the range-separation case, an intermediate value of µ is sufficient. In general, we note that, in the linear case, the calculation of the wave function is by and large equally expensive at all points along the connection since the electron-electron cusp must always be described (except at λ = 0); in the range-separated case, by contrast, the calculation becomes less expensive with decreasing µ since the description of the cusp is then avoided.

Excitation energies and comparison with the range-separated case
The lowest (Rydberg) excitation energies of helium along the linear adiabatic connection are also plotted in Figure 4 (middle and bottom). As for the total energies, the end points are the same as in the range-separated case but the behavior of the energies along the connection is more linear. As in the range-separated case, the 1 P excitation energy does not evolve monotonically with λ, probably because of basis-set limitations.
When the first-and second-order extrapolation corrections are added, a systematic improvement is observed for the triplet excitation energies. With the second-order correction, the physical energies are already reproduced to 1 mhartree at λ = 0. The singlet excitation energies are less affected by the correction but are still overall improved, the amount of interaction required to reproduce the physical limit within an accuracy of 1 mhartree dropping to 50%. Moreover, unlike in the range-separated case, the KS excitation energy also benefits from this correction, which no longer vanishes in this limit. Indeed, at λ = 0, the extrapolated excitation energy matches the results obtained in Ref. 25, using first-order Görling-Levy λ Error in the excitation energies (hartree) Figure 5: Error in the 3 S excitation energy along the range-separated and linear adiabatic connections for the helium atom as functions of µ and λ. The unextrapolated energies are given in full lines and the extrapolated ones in dashed lines. An error smaller than 1 mhartree around the physical limit is given by the colored region.
To compare the range-separated and linear cases, the effects of first-order extrapolation on the lowest excitation energy of helium are shown in Figure 5 in both cases. Without extrapolation, the scaling parameter λ must be greater than 0.95 to reproduce energies within 1 mhartree. By contrast, a small change in µ near zero gives a large change in the energies; thus, a range-separation parameter of 2 bohr −1 is sufficient to ensure the same accuracy. Clearly, the range-separated connection includes the most significant region of the interaction first, whereas the linear connection treats all ranges equally, independently of their importance for the excitation energies. For the KS system, it is obviously better to use the correction obtained from the linear connection as the corresponding range-separated correction vanishes. For a partially interacting system, the comparison is more difficult.

V. CONCLUSION
In this work, we have exploited the asymptotic behavior of the energies of a partially interacting system along the range-separated adiabatic connection to design an energy correction that allows us to extrapolate to the physical energies of the system from its partially interacting energies. The simplest possible extrapolations were obtained by using either only the first-order derivative of the energies with respect to the range-separation parameter at a given point or by using the first-and second-order derivatives.
This extrapolation scheme was tested at the FCI level of theory on the helium and beryllium atoms and on the hydrogen molecule (at equilibrium and at stretched geometry), where it significantly improves the convergence of energies and excitation energies towards their physical limits. Moreover, the improvements are systematic, except at µ = 0 (where the correction is zero by construction) and in a few cases where the partially interacting energies present a bump for small µ. In all cases, with respect to the unextrapolated case, the extrapolation schemes reduce the smallest value of the range-separation parameter required to reproduce the physical energies of the system with a given accuracy by approximately a factor of 2 with the first-order scheme and by a factor of 3 with the second-order scheme. This is of particular relevance for truncated wave functions, as the smaller the range-separation parameter is, the fewer Slater determinants are needed to describe the wave functions with an equivalent accuracy.
Finally, the extrapolation scheme was applied along the linear adiabatic connection, where it also improves significantly the description of the excitation energies along the connection.
All results discussed here were obtained without the use of approximate functionals. The proposed extrapolation scheme should now be tested in a more pragmatic case, where the potential is not obtained by Lieb optimization but from different approximations such as the (semi)local approximations or more interestingly approximations where the long-range behavior of the potential is correct, such as the optimized-effective potential (OEP) [48,49]. The effects of the inclusion of higher-order derivatives and of multiple points on this extrapolation should also be explored.