Exchange energies with forces in density-functional theory

It is with great pleasure that we provide this contribution to the special issue of the Journal of Chemical Physics honoring John Perdew and his work in quantum chemistry. John Perdew has made groundbreaking advancements in developing exchange–correlation functionals within density-functional theory (DFT), which are crucial for an accurate description of the interactions between electrons. DFT,¹ with its many variants,^2–9 is nowadays the workhorse of first-principle simulations in quantum chemistry, solid state physics, and materials science, and John Perdew greatly contributed to this success story. His work has focused on improving the accuracy and efficiency of density-functional calculations by proposing more precise and robust functionals, allowing researchers to study a wide range of chemical and physical properties of materials. Specifically, he has explored fundamental (exact) constraints that an exchange–correlation functional must satisfy to accurately describe the electronic interactions in a system. Enforcing such exact constraints can greatly improve the reliability of approximate functionals.¹⁰ This work follows this general idea by taking up previous suggestions on how to rephrase the exchange–correlation potential in terms of forces. We give known and novel exact force-based constraints and show how to translate these ideas into an efficient numerical scheme.

Most DFT simulations are performed using the Kohn–Sham (KS) scheme,¹¹ where the density of the interacting system is predicted by solving an auxiliary noninteracting system. It is precisely the mentioned exchange–correlation potential that relates the interacting and the noninteracting system through the underlying density-potential mapping v(r) ↔ ρ(r). For a recent review on this mapping in the context of DFT, we point to the work of Penz et al.¹² It is common practice to derive approximations for the (in general unknown) exchange–correlation potential by re-expressing the universal density functional as a sum of noninteracting kinetic, Hartree, and exchange-energy functionals, as well as the unknown correlation energy functional,¹³ and then to assume functional differentiability¹⁴ with respect to the density. While for approximate functionals that are given explicitly in terms of the density, potentials can be determined this way by direct differentiation, for implicit functionals this is no longer possible in general.¹⁵ To make matters worse, it has been shown that the universal density functional of DFT is not functionally differentiable with respect to the usual function spaces.¹⁶ While a generalized definition of functional differentiability (subdifferentiability) is enough to establish the mapping from v-representable densities to potentials,¹⁷ many of the commonly employed rules of differential calculus, such as linearity or the chain rule, might no longer hold in the same way.¹⁸ This fact therefore questions this common way to infer exchange–correlation potentials from exchange–correlation energy functionals. We note that the theoretical setting of an exact regularization procedure is available that renders the involved functionals differentiable^19–21 and that this surprisingly links to the Zhao–Morrison–Parr method for mapping ρ(r) ↦ v(r).²² Importantly, in this work we highlight that also the exchange-only energy is non-differentiable with respect to densities, thus allowing local-exchange potentials only in the form of generalized constructions such as the optimized effective potential (OEP), or, alternatively, leading to an additional vector potential for exchange effects. This vector potential naturally appears in a force-based approach and acts semi-locally on the wave function. This is in contrast to the exchange term in Hartree–Fock that acts fully nonlocally on the wave functions.

From a physical point of view, one can always exchange the description of a many-particle quantum system in its ground state in terms of energies by a description based on forces. Both views have been viable routes toward getting the desired potential. Indeed, the exact exchange–correlation potential of DFT can be expressed directly in terms of the difference in force densities between the interacting and the auxiliary noninteracting system,^23–26 thus bypassing functional differentiation and related issues. A method for deriving DFT potentials from the electric field due to the Fermi–Coulomb hole charge distribution was pioneered by Harbola and co-workers.^27–30 However, this approach misses the kinetic-correlation contribution and thus does not retrieve the full exchange–correlation potential.³¹ This issue was also noted by Holas and March,²³ who first used a force-based approach to give an expression for the exchange–correlation potential of DFT in the form of a (path-independent) line integral. Building upon this important work, Sahni³² was able to extend the method of Harbola.

In this work, we show that a force-based approach is not only conceptually very appealing but also practically relevant. In doing so, we stick to a fully spin-resolved, collinear formulation. Specifically, we show that besides the usual Hartree potential, we can also derive the simple explicit form of the local-exchange potential previously suggested by Harbola and Sahni.²⁷ This potential we show to be directly linked to the exchange force density and it enters a generalized exchange virial relation. A different form of a generalized exchange virial relation is actually discussed in another paper of this special edition.³³ We further find a relation between the exchange and the correlation force densities that takes the form of a novel exact constraint. As we demonstrate, the formulation of the force-based local-exchange potential is consistent with current-density-functional theory (CDFT)^6,34 and we discuss its connection to the time-dependent case. In the context of ground-state DFT, we then show that the explicit force-based local-exchange potential performs similarly to the numerically much more involved optimized effective potential (OEP) approach in exchange approximation. We show that the difference between OEP and force-based local-exchange potential can be connected to the abovementioned exact constraint that exchange and correlation force densities need to fulfill. We finally comment on practical ways on how to treat the remaining correlation force densities. In this, we highlight how the force-based approach provides a route toward numerically inexpensive nonlocal (in how it depends on the density) and nonadiabatic functionals that also act semi-locally on the wave function if they contain a vector-potential contribution.

Let us now comment on some exact constraints for the force densities. If, for the sake of consistency, just a single particle with wave function φ(rσ) is considered, then it directly follows F_W[φ] = 0 (no self-force) and naturally F_c[φ, φ] = 0. Note that F_W[φ] = 0 can also be deduced from Eq. (14). Moreover, in the one-particle case, because of Ψ = Φ = φ the kinetic-correlation and the interaction-correlation force densities must vanish independently. We further remark that these self-interaction properties are directly related to the corresponding expressions for the energy, which serve as a basis for the construction of self-interaction corrections, as pioneered by Perdew and Zunger.⁴⁵ A similar scheme could thus be developed on the basis of forces.

A further exact constraint that holds locally for the exchange and correlation vector fields is derived at the end of Sec. IV.

The potential v_fx is therefore the exchange potential that originates from only the longitudinal part of the exchange vector field f_x = F_x[Φ]/ρ. We will come back to this point and its implications for density-functional approximations below. To complete the picture, the missing correlation potential is given uniquely in terms of the (unknown) force density difference F_c[Φ, Ψ] from Eq. (11) and a simple Coulomb integral [see Eq. (12)]. For the force-based local-exchange potential given by Eq. (19), a numerically more convenient form in terms of the Slater exchange potential plus correction terms can be derived (see  Appendix A). It also obeys the usual coordinate scaling relations (see  Appendix B).

On the practical side, if one is only interested in finding a local potential that minimizes the exchange energy, then the common OEP approaches will usually perform better than the force-based local-exchange potential (see Table II in  Appendix C for comparison to Hartree–Fock results). This is by design, since the exchange-only OEP procedure is precisely such that it seeks the local potential v_OEPx that minimizes the energy with an uncorrelated state.⁵⁴ Herein, the state is always a Slater determinant constructed from the orbitals of a one-particle Hamiltonian with the chosen potential v_OEPx. Note, however, that due to the restriction of the OEP to local potentials, the obtained energies will be higher than the Hartree–Fock results that allow for nonlocal potentials. Furthermore, as it is clear from the previous discussion, the common OEP procedures applied to the exchange energy do not give the correct exchange force density. Instead they will generate a purely longitudinal vector field that is not related to the exchange force density in a direct manner. We usually lose control over the connection between the energy terms and the corresponding force densities (which leads, among others, to a violation of the virial relation). An important exception is the exchange-only local-density approximation, where the connection still holds as can be shown directly. The same holds for correlation approximation. Any approximate correlation energy can always only lead to a longitudinal vector field via the corresponding (generalized) energy derivative, while an approximation based on forces will usually include a transverse component. This means that there is no strict connection between the energy-based and force-based approximations. To put it differently, if we want to build approximations in DFT based on Hartree and exchange terms beyond the local-density approximation, we have to decide whether we use the exchange energy or the exchange force density. Both strategies will only agree when we use the exact exchange and correlation terms together.

Finally, let us comment on the homogeneous-density limit. In the work of Tchenkoue et al.,²⁶ it was demonstrated how the usual Slater Xα⁵⁵ and local-density approximation (LDA)⁴⁴ formulas for the local-exchange potential can be derived directly from the exchange force expression of Eq. (14). Since f_x(rσ) is purely longitudinal for a homogeneous density, the exact same derivation can also be started immediately from the local-exchange potential expression of Eq. (19). A related derivation of the same fact based on the second-order gradient expansion of the exchange-hole density was already given in the work of Wang et al.⁵³ This directly connects the most fundamental functional approximations of DFT with the present formalism.

Another advantage of the force-based approach is the inherent compatibility to CDFT and time-dependent DFT (TDDFT). The generalized exchange virial relation of Eq. (24) highlights the connection of the force-based approach to CDFT. If besides the density ρ we also intend to control the current density j, then we would need a transverse exchange–correlation vector potential as well, where α_fx contributes to the exchange vector potential.²⁶ We even find that v_fx and α_fx can be chosen to be the local-exchange potential of CDFT and of time-dependent CDFT.²⁶ This makes v_fx nicely compatible with this variant of DFT.

Finally, let us shift attention back to the time-independent setting again. Therein, besides Eq. (8), the exact ground-state exchange–correlation potential and force density still also obey Eq. (27). This gives rise to a different version of the local-exchange potential. Here, we will not investigate this alternative force-based formulation further but will compare these different definitions in a forthcoming publication. It however highlights a route to more, possibly useful conditions: Higher-order equations of motions bring with them new exact constraints.

Finally, we consider the differences between the force-based approach and the energy-based approach in practice, with a focus on the effects of the transverse part of the exchange vector field f_x expressed through the vector potential α_fx. For this, we solve the KS equation in exchange approximation (FBEx), i.e., we take v_fx from Eq. (19) and assume v_fc = 0 for the total v_Hxc in Eq. (12) in every KS iteration step, and check how this performs in comparison to common exchange approximations. In this investigation, we do not yet employ the transverse part of f_x somehow beneficially. Yet, involving only the longitudinal component of f_x in the calculation is equivalent to considering the full f_x plus the transverse part from f_c since Eq. (25) holds as an exact constraint. A force-based approximation focusing purely on exchange effects thus would need to also consider the transverse contribution from the exchange vector field. To summarize, there are two possible viewpoints on this approximation that are equally justified: When considering only the v_fx exchange potential then exchange effects from α_fx are missing, or alternatively, that this procedure additionally includes correlation effects from α_fc.

We have implemented the force-based local-exchange potential in the real-space code Octopus⁵⁸ and ran simulations for a set of atoms in closed-shell configurations using norm-conserving pseudopotentials,⁵⁹ a grid spacing of 0.15 bohr, and a radius of 10 bohrs for Be and Ne, a radius of 12 bohrs for Mg, Ar, and Zn, and a radius of 14 bohrs for Ca. We found that the FBEx potential performs similar to the much more involved OEP in exchange approximation (OEPx) or its further approximation OEPx-KLI⁶⁰ (see Fig. 1). While the pure Slater, FBEx, and OEPx-KLI potentials all share the same computational scaling as Hartree–Fock, the OEPx method only works as an iterative procedure and is more costly. Furthermore, we demonstrate that the local-exchange potential adheres to the virial relation of the form of Eq. (22) up to numerical inaccuracies (see Table I) because of spherical symmetry, while the OEPx and the OEPx-KLI violate this relation. This numerically confirms that the exchange functional is not functionally differentiable. Further numerical tests and comparisons, also for small molecules, can be found in  Appendix C.

In fact, Fig. 1 shows that the FBEx and the OEPx potentials are almost identical, apart from the small “bumps” that are an indication of the shell structure of the atoms. The Slater potential also does not capture them and a suitable correction for it based on the kinetic-energy density is available.⁶¹ Due to the use of pseudopotentials in our simulations displayed in Fig. 1, we see here either one or no bump. To check that the differences between these two potentials are indeed only present at the shells of the atoms, we also performed all-electron calculations for Ne and Ar using a grid spacing of 0.05 bohr and a radius of 14 bohrs and we obtain that the potentials differ only at the location of the bumps, see the top panels of Fig. 2. The lower panels of Fig. 2 show the difference between OEPx-KLI and FBEx together with ‖α_fx‖, i.e., the transverse part of f_x. The same comparison is conducted with OEPx for those atoms where a bump is visible despite using pseudopotentials, see Fig. 3. In each case, we find that the FBEx force has a nonvanishing transverse part only at the position of the bumps and that the norm of α_fx follows a pattern similar to the difference between the FBEx and OEPx potentials, clearly showing that there is a connection between the transverse part of the force and the bumps of the OEPx potential. In fact, we interpret our result in the following way: The bumps appear in the OEPx potential as the procedure tries to impose a longitudinal vector field at places where the exchange vector field actually has a transverse component. Thus, even if we do not employ the transverse part of the forces explicitly, they contain physical information (related to the shell structure of atoms) that can be potentially used for more advanced approximations. For instance, this feature is related to the correlation forces f_c, since α_fx needs to precisely compensate α_fc. In this manner, we get local information about the correlation vector field that could provide quite stringent constraints on future approximations.

A further comparison of the FBEx and OEPx-KLI methods with inversion procedures that yield the full exchange–correlation potential is performed in  Appendix D.

Considering all the different insights obtained by this investigation, we want to highlight two specific results that we deem important for the future of force-based approximations. On the one hand, we have seen that the transverse part of the exchange vector field contains important physical information. It stands to reason that the standard OEP procedure tries to turn these transverse parts into longitudinal contributions of the corresponding OEP potential, which are responsible for the appearance of the “bumps.” An obvious way of including these contributions is to employ an auxiliary system that also contains a vector potential instead of the usual KS system with only a scalar potential. Using the beneficial connection of the force-based approach to CDFT, the corresponding exchange vector field is given via a nonlinear partial-differential equation.²⁶ This paves the way to obtain a semi-locally acting vector potential in the context of electronic ground-state DFT.

On the other hand, for time-dependent DFT, we have seen that the appearance of the transverse vector field implies nonadiabaticity. That is, if we solve the corresponding Sturm–Liouville Eq. (27) instead of the Poisson Eq. (12), we automatically get a nonadiabatic functional based on force densities. These two aspects make the force-based approach quite promising to find more accurate yet numerically inexpensive approximations within density-functional theories. It is even relatively easy to extend the present approach to other variants of DFT, for instance to forms that include noncollinear magnetism and spin–orbit coupling.^62–65 Using the corresponding equations of motions for the current density,⁶⁶ one can apply the same Hartree exchange and correlation force density decomposition and hence is able to derive the corresponding potentials also for this case. Furthermore, in order to address the still unknown correlation force density we highlight that the transverse part of the exchange vector field provides us with local constraints on approximate correlation force densities. In the correlation force density, the interaction part F_W[Ψ] − F_W[Φ] can be expressed by the correlation hole, while the kinetic part F_T[Ψ] − F_T[Φ] can be expressed as the difference between the interacting and the noninteracting one-body reduced density matrix close to the diagonal.⁴⁰ Approximations can then be tested by comparing to the transverse exchange vector field. From this perspective, the success of LDA-based approximations can be explained by the fact that already on the exchange level no transverse forces appear, such that the virial relation is fulfilled, and hence adding purely longitudinal correlations obeys the zero transverse vector field constraint of Eq. (9). Alternatively, one can start from approximated correlated reduced density matrices and derive the corresponding forces. One can therefore either try to build approximate models based on physical intuition,⁶⁷ derive expressions for these terms for specific cases (e.g., the homogeneous limit) from wave function methods potentially augmented by modern machine-learning techniques,⁶⁸ or devise perturbative expansions on top of the KS Slater determinants. Even though the force densities are three-dimensional vector fields and thus more involved than energy expressions, the previously successful application of the aforementioned approaches to construct correlation energy functionals makes it plausible that similar methods are well applicable to the force-based approach to KS-DFT.

In conclusion, we have shown that defining the Hx potential and energy of KS-DFT by forces not only is conceptually beneficial but also has certain advantages in practice over the common energy-based approach. It is numerically straightforward to construct the corresponding potential from a given force density, the method allows to avoid various problems of the energy-based approach such as determining implicit functional derivatives, and it further provides an explicit form for the local-exchange potential and exchange energy from the exchange force density. This force-based local (in the sense on how it acts on the wave function) exchange approximation depends nonlocally on all other points and all occupied orbitals and is numerically as cheap as the Slater potential. The non-explicit correlation potential is defined uniquely by the correlation force density and in contrast to the energy-based approach, the role of correlations in compensating the transverse part of the exchange vector field is transparent. It is seen that the exchange vector field provides local information about the properties of the correlation vector field. We also have a straightforward connection to the current density variant of DFT and to the time-dependent case. Furthermore, the approach can be seamlessly applied to atomic, molecular, and solid state systems. We showed numerically that the well-known bumps of the OEPx potential are connected to the transverse exchange vector field and with this also to the correlation vector field due to the exact constraint that the transverse exchange vector field is exactly compensated by the corresponding correlation effects. We think, following the ideas of John Perdew and others, that such local exact constraints are a good starting point to help in devising correlation force density approximations, in DFT and its variants.

This work was supported by the European Research Council (Grant No. ERC-2015-AdG694097), by the Cluster of Excellence “CUI: Advanced Imaging of Matter” of the Deutsche Forschungsgemeinschaft (DFG) – EXC 2056–Project ID 390715994, and the Grupos Consolidados (Grant No. IT1249-19). M.P., M.A.C., and A.L. have received funding from the ERC-2021-STG under Grant Agreement No. 101041487 REGAL. M.A.C. and A.L. were also supported by the Research Council of Norway through funding of the CoE Hylleraas Centre for Quantum Molecular Sciences Grant No. 262695 and CCerror Grant No. 287906.

The authors have no conflicts to disclose.

Nicolas Tancogne-Dejean: Conceptualization (equal); Formal analysis (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (equal). Markus Penz: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Andre Laestadius: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Mihály A. Csirik: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Michael Ruggenthaler: Conceptualization (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Angel Rubio: Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This form has a few advantages. First, we only need to solve one Poisson equation and compute one gradient per pair of indices i, j. Therefore, the numerical cost only increases by one gradient per pair of indices i, j compared to the Slater potential. Furthermore, it provides an analytical expression for beyond-Slater approximations and might serve as the starting point for the development of novel functionals. Finally, from this expression it is also clear that in the single orbital case, $ρij*(r′σ)/ρ(r′σ)$ is uniformly equal to 1 and that then the second term vanishes.

Here, we show further numerical comparisons of the force-based local-exchange potential to well-established exchange potentials in DFT. First, we investigate how the force-based local-exchange potential compares to the Hartree–Fock energies. Since the force-based local-exchange potential is not derived directly from the exchange-energy expression of Eq. (21), it is not designed to approximate the nonlocal Hartree–Fock exchange-energy expression. Still, the resulting energies of the force-based local-exchange potential determined from Eq. (21) together with the respective orbitals (see Table II) are in good agreement with the Hartree–Fock exchange energies. Note that due to the nonlinear core correction from the pseudopotential and the larger number of valence electrons, the results for Zn show a larger discrepancy with Hartree–Fock.

In Table III, we further report the eigenvalue of the highest occupied orbital for different exchange functionals. While OEPx-KLI and OEPx are yielding similar ionization energies as Hartree–Fock, within a meV precision, the force-based local-exchange potential leads to only slightly different results. The Slater potential shows a stronger deviation from the Hartree–Fock values.

Table IV lists the difference in exchange energy computed from the orbitals and the virial relation of Eq. (22) for small molecules. This shows that for non-spherically symmetric systems, this virial relation is not respected by the force-based local-exchange potential either. For N₂, we employed a N–N distance of $1.09769Å$⁠. For CO₂, we considered a C–O bond length of $1.16Å$⁠. For CH₄, we considered a C–H bond length of $1.087Å$⁠. In all cases, we employed a grid spacing of 0.15 bohr and a simulation box made of atom-centered spheres of radii 12 bohrs.

The corresponding ionization energies for these molecules are given in Table V. Similar to the atomic case, we find that the force-based local-exchange potential performs much better than the Slater potential and yields ionization energies close to the ones obtained from Hartree–Fock or OEPx.

The fact that we only include the longitudinal part of the exchange force density can be viewed as implicitly including a correlation force density that imposes Eq. (25). It is therefore interesting to compare not only to OEPx results but also to the exact exchange–correlation potential. For the atoms considered in the present work, this was done by the mean of the Kohn–Sham inversion procedure, for instance based on Green’s function densities,⁶⁹ or from CI densities.⁷⁰ The comparisons are shown in Fig. 4. From these results, it is clear that the implicitly included correlation part does not seem to agree with the exact potential, as the bumps representing the atomic shells are still a dominant feature in this potential.