Accurate real-time evolution of electron densities and ground-state properties from generalized Kohn-Sham theory

The exact time-dependent Kohn-Sham (KS) exchange-correlation (xc) potential is not easily approximated due to its nonphysical properties; the burden of capturing xc effects solely within a multiplicative potential gives rise to pathological features, such as spatial steps. The generalization of the KS approach allows one to lessen this burden via an additional, more physical non-local potential. We present unrestricted Hartree-Fock (HF) theory -- in which each electron occupies a distinct orbital regardless of its spin -- as well as restricted HF, each with a corresponding additional multiplicative correlation potential which in principle ensures the single-particle electron density is exactly equal to the many-body density. The exact form of the local correlation potential when all electrons are permitted to occupy distinct orbitals is largely free from the pathological features which are present in the exact KS-xc potential. Hence we find that an (adiabatic) local density approximation to the correlation potential yields accurate ground-state properties and real-time dynamic densities for one-dimensional few-electron test systems -- we compare our results to the exact many-body quantities.

Models which reliably describe excited many-body systems for a low computational cost have remained elusive within solid state physics and quantum chemistry, despite the growing need for them. Density functional theory [1] (DFT), within the Kohn-Sham (KS) approach [2], is an extremely popular method for ground-state electronic structure calculations owing to its computationally efficiency and accuracy for most solids [3]. However, for modeling systems with strong electron localization, such as molecules, the spatially local exact exchange-correlation (xc) potential of KS theory has been shown to exhibit important features that depend non-locally on the density [4][5][6][7][8][9], which common approximations fail to capture [10][11][12][13][14][15][16], e.g., sudden changes in the level of the potential termed 'steps'. These missing features lead to an inaccurate description of the system, e.g., when the atoms of a diatomic molecule are dissociated [12], and within time-dependent DFT [17] (TDDFT) for dynamic systems in which strong currents flow [18][19][20][21][22].
Accurate dynamic densities are crucial for predicting currents and electronic properties of molecules which are perturbed by an external field. In particular, when acting as a molecular junction, the system is well beyond linear-response and steady-state regimes. Yet commonly used approximations within DFT predict current-voltage characteristics within molecular electronic systems which are incorrect by orders of magnitude [23]. An accurate description of excited electronic systems requires a reliable model of exchange and correlation. The spatially non-local potential of Hartree-Fock (HF) theory offers an exact description of exchange to which correlation effects may be introduced, for example via the commonly used GW approximation [24] within many-body perturbation theory (MBPT) or hybrid density functionals [25,26].
However, these methods do not yield reliable dynamic densities: the GW approximation is widely used to calculate photoemission properties of molecules [24,[27][28][29][30][31][32]. These calculations can be computationally expensive and are inaccurate when electron correlation is strong [33][34][35][36][37][38]. For time-dependent systems, the computational challenges within MBPT prevents its use for practical calculations in many cases. Hybrid density functionals have been successfully used to calculate quasiparticle energies and ground-state densities [26], yet do not offer real-time densities at present.
In principle KS theory within TDDFT yields an exact time-dependent density; however, in practice the exact time-dependent xc potential is difficult to approximate [18][19][20][21][22]. TDDFT has been very successfully employed within the linear response regime, however this is limited to weakly perturbed systems, e.g., photo-absorption, and the approximations employed are also known to be unreliable for describing important phenomena, such as charge transfer [39][40][41][42][43][44].
Generalized Kohn-Sham (GKS) theory [45] offers a computationally cheap way of incorporating correlation effects into HF theory, e.g., hybrid functionals. GKS establishes that for a given non-local potential, such as Fock exchange, there exists a unique spatially local potential which ensures the single-particle density equals the many-body density exactly [45]. The form of this local potential thus depends on the choice of non-local potential. Hence with an appropriate choice, the corresponding local potential can be more amenable to approximation. Recently GKS has been extended to systems undergoing excitations [43,46].
In this paper we consider two choices for the non-local potential. The first is the Fock operator of restricted Hartree-Fock (RHF), which leads to a set of RHF-Kohn-Sham (RHFKS) equations, derived in Ref. 45. For the second we employ unrestricted Hartree-Fock (UHF) [47], which leads to a set of UHF-Kohn-Sham (UHFKS) equations, which we derive below.
Within GKS the electrons must be described by a single Slater determinant (SD), Φ. We employ an 'unrestricted SD' in which electrons with different spins occupy different single-particle orbitals, i.e., ψ k (x, σ) = φ k (x)γ k (σ) where γ · = α for up-spin electrons and γ · = β for down-spin electrons. We construct the SD from the orbitals {ψ k (x, σ)}. This is in contrast to the SD of RHF in which two electrons with opposite spins occupy the same orbital. We then define the functional whereT is the kinetic energy operator andÛ is the electron-electron interaction operator. From Eq. (1) a unique density functional can be defined via the con-strained search formalism of DFT [48]: where the minimization searches over all SDs that yield the electron density n. This then allows us to define our correlation energy functional as the difference between Q S [n] and the Hohenberg-Kohn functional (which captures all xc effects exactly) [1]: where Ψ[n] is the ground-state many-body wavefunction which yields the electron density n. Because UHF captures exchange and static correlation effects [49], the correlation energy defined by Eq. (3) corresponds to dynamic correlation.
The total ground-state many-body energy can be written in terms of these functionals, as such where v ext is the external potential of the many-body system. Finally the set of single-particle UHFKS equations can be derived via a minimization of this energy employ-ing Lagrangian multipliers to ensure orthogonality of the single-particle orbitals, {φ k } (see Supplemental Material for more details): (atomic units are used throughout this paper; = e = m e = 4πε 0 = 1) where and All electrons experience the same Hartree potential (v H ) and spatially local correlation potential, v c . In principle the single-particle density n s (x) = i,γ |φ γ i (x)| 2 exactly equals the many-body density, n(x); however, in practice the local correlation potential must be approximated, to which we now turn our attention.
In this paper we model one-dimensional systems consisting of a two opposite-spin electrons in their groundstate and for time-dependent scenarios. We focus on these small yet challenging model systems as all of the exact quantities can be numerically determined, including the fully-correlated many-body wavefunction, by solving the (time-dependent) many-body Schrödinger equation for any external potential, v ext (x, t), with the appropriately softened 1D Coulomb interaction (|x − x | + 1) −1 [50] by employing our iDEA code [22]. We also model our systems with (time-dependent) R-and UHF theory in order to assess their accuracy for various scenarios. Finally, we use our exact solution to the Schrödinger equation to develop a local density approximation (LDA) to the correlation potential of R-and UHFKS. We term these approximations R-and ULDA+ to distinguish them from the usual LDA employed within standard KS theory. We assess their accuracy for various ground-state and time-dependent systems.
First, we calculate the exact v c (x) of RHFKS theory and that of UHFKS theory for ground-state one-dimensional model of the hydrogen molecule (H 2 ). As we have access to the exact many-body density for H 2 we can 'reverse-engineer' the R-and UHFKS equations in order to find the corresponding exact correlation potential in each case.   Fig. 3). The two potentials are indistinguishable, indicating the lack of static correlation at this atomic separation. As the atoms are dissociated these two correlation potentials diverge starting at the point where UHF breaks the spin symmetry (Coulson-Fischer point) [51] . The form of the exact v c (x) of RHFKS is non-local with respect to the density and is important to correctly distribute the electronic charge, in a similar fashion to KS theory (for an asymmetric system the peak in the center would be a step [9]; see Supplemental Material); see Fig. 1(b). These types of features are notoriously difficult to capture in approximate methods. However, the exact correlation potential of UHFKS does not contain these features, and hence is more amenable to a local approxi-mation. We attribute this to the 'nearsightedness' [52] of the exchange potential felt by each electron within UHF, which originates from the use of different orbitals for different spins [53] -because each electron is free to independently experience the Hartree potential of the other electron, there is no need for the spatial peak (or step) in the local potential in order to correctly distribute the change density. This is demonstrated in Ref. 52. As a result, we find that an LDA to the correlation energy, and thus the spatially local correlation potential, yields accurate electron densities and energies for our model systems; see Figs. 3 and 4 below. We also find that the adiabatic LDA (ALDA) yields accurate real-time densities within time-dependent UHFKS theory; see below. In order to construct an LDA for the correlation energy of U-and RHFKS we use a set of 'slab systems' [54] which are homogeneous in the center of the system and tend to zero towards the edge -as in Ref. 55. We vary the height of the slabs (n) and calculate the correlation energy (E exact − E HF ) for each system [56]; Fig. 2. Once an LDA for each correlation energy is constructed we apply it to its own training systems and find small errors in E c (∆E c ) owing to the inhomogeneity of the slab systems. We use these errors to calculate a refined LDA+: ε c (n) → ε c (n) − ∆E c (n)/2. Figure 2(b) shows the correlation energy per electron for the range of slab systems (ε c ) as a function of the density, n. The form of the 'restricted' correlation energy, i.e., that which corresponds to RHFKS, is similar to the 'unrestricted' correlation energy for high density regions because for these systems static correlation is negligible. However, for the lower density regions the two energies differ to a large degree. Note that the restricted correlation energy is much larger than the unrestricted correlation energy because it has the burden of capturing static correlation. R-and UHF are compared against the many-body exact. As expected, RHF is incorrect for all separations and UHF is inaccurate at and around the bonding length owing to an absence of dynamical correlation effects. Our RLDA+ only slightly improves the energy. Our ULDA+ gives a large improvement upon UHF at and around the bonding length and slightly worsens the performance of UHF in the dissociation limit.
We now calculate the molecular energy of H 2 as the atoms are separated employing R-, UHF, R-and ULDA+ and compare them to the exact case. Figure 3 shows, as expected, that the RHF energy is inaccurate owing to a complete absence of correlation. UHF correctly gives the dissociation energy of the molecule by capturing static correlation [49] but yields the incorrect energy at and around the bonding length (∼ 1.4 a.u.) owing to the absence of dynamic correlation in the approximation. At the bonding length R-and UHF yield exactly the same energy; the error in the total electron energy is ∼ 1.3% [57]. The error in the ionization potential (IP) predicted by U-and RHF is 2.1%. At an atomic separation of 7.2 a.u. the error in the IP of RHF is 34.3% and 1.7% for UHF. ULDA+ yields an accurate molecular energy at the bonding length with an error of ∼ 0.5% and introduces an error as the atoms are separated (∼ 2%) because of the use of an LDA to v c ; as the electrons localize each to an H atom, the approximate correlation potential introduces a 'self-correlation error'. In principle this error could be removed via a more sophisticated approximate to v c [38]. Whereas RLDA+ gives a relatively poor total energy at the bonding length 1.6% and an only slightly improved total energy when the molecule is stretched 9.0% (compared to the RHF error of 15.2%). Because the LDA for the restricted correlation energy aims to add both static and dynamic correlation, and there is negligible static correlation in the systems with relatively small atom separations, RLDA+ yields relatively inaccurate energies by introducing spurious static correlation. This issue is not present for ULDA+ as it only attempts to introduce dynamic correlation which is present in these systems. (Note that R-and ULDA+ employ our LDA to each correlation energy density functional in order to calculate the total electron energy, as is appropriate within DFT.) Next we turn to the ground-state density. We compare the density from R-and ULDA+ to R-, UHF and the exact at the bonding length (see Fig. 4(a)) and when the atoms are dissociated; see Fig. 4(b). When the atoms are relatively close together, R-and UHF yield the same density (like for the energy). R-and ULDA+ also yield very similar densities, with ULDA+ performing slightly better. On the other hand, when the atoms are separated the RHF and RLDA+ densities are inaccurate. The fact that RLDA+ only slightly improves the electron density compared to RHF shows how ineffective an LDA to the restricted correlation potential is. This is not surprising as the exact correlation potential depends non-locally on the density (shown above). UHF is extremely accurate when the atoms are dissociated, as already discussed, and ULDA+ is only slightly worse for the density demonstrating that the self-correlation error introduced by our LDA to the correlation potential does not have a detrimental effect.
Let us now assess the accuracy of employing the ALDA within R-and ULDA+ to calculate the time-dependent density and current when the H 2 molecule is perturbed by an electric field (−0.03x). The time-dependent version of Eq. (5) isĥ whereĥ γ is the single-particle Hamiltonian given by Eq. (5) but with the time-dependent potentials. (We do not derive this equation as the existence of a unique local potential cannot be ensured in general [46], but is instead an ansatz.) Figure 5(a) shows the electron density at t = 30 a.u. The density sloshes back and forth within the molecule and generates a current; see Fig. 5(b). The current is the clearest indicator of how accurate each approximation is: 4. (a) The H2 molecule at the bonding length (atomic separation ≈ 1.4 a.u.). U-and RHF densities are indistinguishable: both are inaccurate in the central region when compared to the many-body exact. The decay of the approximate densities is also incorrect. R-and ULDA+ yield accurate densities. (b) The dissociated H2 molecule (separation = 8.0 a.u.). UHF density is extremely accurate. The RHF density is poor when compared to the exact owing to the absence of static correlation. RLDA+ does not improve the density much over RHF, showing the in inability to capture static correlation with a local approximation. The ULDA+ density is only slightly worse than UHF owing to the self-correlation error (see text) introduced by the local approximation to the correlation energy. RLDA+ is the worst. R-and UHF perform identically and yield a reasonably good account of the electron current. The best method is ULDA+; although the current is by no means exact, it is quantitatively correct throughout the whole simulation (t = 30 a.u.). The relative errors of the time-dependent charge densities as the system evolves is shown in Fig. 7(a).
When the stretched molecule is perturbed, the story is similar: RLDA+ performs the worst, followed by RHF. However, for this system UHF is the most accurate followed by ULDA+. This is not a surprising result; it follows from the advantageous nearsightedness of UHF's non-local potential which leads to near exact densities in the ground state. Like for the ground state, the LDA potential has not significantly worsened the ULDA+ result. Again, in principle this error can be reduced by more advanced approximations to the correlation potential. the system evolves is shown in Fig. 7(b). (Videos of all time-dependent results can be found in the Supplemental Material.) In conclusion, we have shown that the local potential which ensures an exact electron density, which is characteristic of Kohn-Sham (KS) theory, can be more easily approximated when a physically justified non-local potential is employed in conjunction within generalized KS theory. We find that the exact local correlation potential which corresponds to unrestricted Hartree-Fock (UHF) theory within generalized KS theory is missing the pathological features which are present in the exact exchangecorrelation potential of standard KS theory. This is not the case for the local correlation potential which corresponds to restricted HF theory. As a result, we find that our LDA to the correlation potential of UHFKS theory yields accurate ground-state densities and energies and can even yield relatively accurate currents when the system is perturbed by an electric field. We demonstrate this for one-dimensional systems which consist of two electrons, thus allowing us to compare our calculations to the exact many-body solution.
Our LDA to the correlation potential of UHFKS introduces an intrinsic spurious self-correlation error. We find that this leads to inaccurate energies and densities The same is the case for the current. Now UHF yields the most accurate current but ULDA+ yields a quantitatively good description of the current throughout the simulation. This current could be improved by a more advanced approximation to the correlation potential.
when the atoms of H 2 are dissociated, although these errors are significantly small relative to other approximate methods. This error can in principle be removed by a more advanced approximation to the correlation energy. Yet overall, our 'ULDA+' yields accurate dynamic densities and currents as well as ground-state properties even when correlation is strong.
We thank the University of York for computational resources.