We propose exchanging the energy functionals in ground-state density-functional theory with physically equivalent exact force expressions as a new promising route toward approximations to the exchange–correlation potential and energy. In analogy to the usual energy-based procedure, we split the force difference between the interacting and auxiliary Kohn–Sham system into a Hartree, an exchange, and a correlation force. The corresponding scalar potential is obtained by solving a Poisson equation, while an additional transverse part of the force yields a vector potential. These vector potentials obey an exact constraint between the exchange and correlation contribution and can further be related to the atomic shell structure. Numerically, the force-based local-exchange potential and the corresponding exchange energy compare well with the numerically more involved optimized effective potential method. Overall, the force-based method has several benefits when compared to the usual energy-based approach and opens a route toward numerically inexpensive nonlocal and (in the time-dependent case) nonadiabatic approximations.

## I. INTRODUCTION

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,^{1} 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.^{10} 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,^{11} 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.*^{12} 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,^{13} and then to assume functional differentiability^{14} 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.^{15} 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.^{16} While a generalized definition of functional differentiability (subdifferentiability) is enough to establish the mapping from *v*-representable densities to potentials,^{17} many of the commonly employed rules of differential calculus, such as linearity or the chain rule, might no longer hold in the same way.^{18} 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**).^{22} 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.^{31} This issue was also noted by Holas and March,^{23} 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^{32} 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.^{27} 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.^{33} 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.

## II. FORCE-BASED KOHN–SHAM SETTING

*N*-particle Hamiltonian [in Hartree atomic units $e=\u210f=me=(4\pi \u03f50)\u22121=1$], first in a time-dependent setting while we later switch to ground states. We have

*v*(

**r**

*σ*,

*t*) is the external, spin-resolved one-particle potential at time

*t*. Note that while the external potential can act separately on the spin components, we here do not take an external magnetic field nor spin–orbit coupling into account. We at the end comment how to extend the present force-based formalism to these cases as well. For antisymmetric wave functions Ψ(

**x**

_{1}, …,

**x**

_{N},

*t*), where

**x**

_{k}= (

**r**

_{k}

*σ*

_{k}), we define the spin-resolved

*p*th-order reduced density matrix as

*ρ*(

**x**,

*t*) =

*ρ*(

**r**

*σ*,

*t*) =

*ρ*

^{(1)}(

**r**

*σ*,

**r**

*σ*,

*t*) to express the (paramagnetic and spin-resolved) current density as

^{26,35}(also called “local force-balance equation”) as

**r**′ −

**r**|

^{−1}) indicates that the gradient only acts on the Coulomb interaction term. Those force terms can be linked directly to the quantum stress tensor

^{36}that includes information about the atomic shell structure.

^{37}Equation (4) has been the primary starting point for inquiries in time-dependent DFT (TDDFT). Among other things, it was used to provide a mapping from densities to potentials,

^{38}to analyze features of the time-dependent exchange–correlation potential,

^{39}to get exact constraints as well as formulations for nonadiabatic approximate functionals,

^{40,41}and to reformulate KS-TDDFT in terms of the second time derivative of the density.

^{42}While here we focus on the ground-state problem, some consequences for the time-dependent case will be discussed further in Sec. V.

**F**

_{W}[Ψ] and

**F**

_{T}[Ψ]. The auxiliary, noninteracting KS problem is controlled by the Hamiltonian $H\u0302s=T\u0302+V\u0302[vs]$, including a different external potential

*v*

_{s}(

**r**

*σ*,

*t*), and has a Slater determinant solution Φ. Analogous to Eq. (4), we then have for the auxiliary system

**j**

_{s}.

*ρ*(

**r**

*σ*) =

*ρ*

_{s}(

**r**

*σ*). In the ground state, it also holds that

*∂*

_{t}

**j**(

**r**

*σ*) =

*∂*

_{t}

**j**

_{s}(

**r**

*σ*) = 0 and we find with the definition of the Hartree exchange–correlation (Hxc) potential

*v*

_{Hxc}(

**r**

*σ*) =

*v*

_{s}(

**r**

*σ*) −

*v*(

**r**

*σ*) that

**F**

_{Hxc}for each spin channel. By virtue of the Hohenberg–Kohn theorem

^{12,43}and assuming non-degeneracy of the ground states for simplicity, the Slater determinant Φ as well as the interacting wave function Ψ are given solely and uniquely in terms of the density, which makes all the force densities determined by the density only. Equation (8) implies that

^{2}

*v*

_{Hxc}= −∇·

*f*_{Hxc}by applying the divergence and solve for

*v*

_{Hxc}using the corresponding Green’s function for the spatial domain $R3$,

**F**

_{W}[Φ] is the Hartree exchange (Hx) force density and

**F**

_{c}[Φ, Ψ] the correlation force density. If desirable, the correlation part can be split again into a kinetic-correlation contribution

**F**

_{T}[Ψ] −

**F**

_{T}[Φ] and an interaction-correlation contribution

**F**

_{W}[Ψ] −

**F**

_{W}[Φ]. The partition of Eq. (11) leads to the respective force-based potentials,

*v*

_{fHx}and

*v*

_{fc}, that add up to the exact Hxc potential,

*f*_{Hx}=

**F**

_{W}[Φ]/

*ρ*and

*f*_{c}=

**F**

_{c}[Φ, Ψ]/

*ρ*. Since the Hx force density is given in terms of the KS wave function only, we know this part explicitly and we can in principle calculate the exact force-based Hx potential for a given KS wave function.

*φ*

_{k}(

**r**

*σ*). We can then express

^{44}

**F**

_{H}(

**r**

*σ*) = −

*ρ*(

**r**

*σ*)∇

*v*

_{H}(

**r**), we read off the (spin-summed) Hartree potential

^{27}This relation provides an important link from forces, or approximations to them, back to the respective energies.

## III. FORCE-BASED EXACT CONSTRAINTS

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.^{45} A similar scheme could thus be developed on the basis of forces.

^{46}in the force-based formulation for the ground state take the simple form

**F**

_{W}[Φ] to the full

**F**

_{Hxc}[Φ, Ψ] available, we can tighten these constraints further. Equation (14) yields an antisymmetric integrand in

**r**

*σ*,

**r**′

*σ*′ in Eq. (17) that must be invariant under the exchange of

**r**

*σ*↔

**r**′

*σ*′. This means that Eq. (17) holds for

**F**

_{W}[Φ] independently and thus we also receive exact constraints for just the correlation force,

*et al.*

^{40}

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

## IV. DISCUSSION OF THE FORCE-BASED LOCAL-EXCHANGE POTENTIAL

*v*

_{fHx}=

*v*

_{H}+

*v*

_{fx}. Here, we used the usual definition of the exchange-hole density

^{44}(without a factor $12$ since it is spin-resolved) given by

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).

*E*

_{x}[

*ρ*] is defined as a

*density*functional by invoking the usual mapping

*ρ*↦ Φ. As was pointed out by van Leeuwen,

^{15}for an implicit density functional the (generalization of the) functional derivative is not straightforward and does not exist in general. On the other hand, if the functional derivative would exist, then

*by construction*it obeys a virial relation of the form (see Appendix B)

^{47,48}that needs to assume Fréchet (total functional) differentiability to allow for the application of the functional chain rule.

^{49,50}Yet, OEP exchange potentials, in accordance with non-differentiability of the exchange-energy functional, in general do not obey Eq. (22). Sometimes, this relation is additionally imposed, e.g., in Fritsche and Yuan

^{51}(also compare Table I); however, this will not restore differentiability. Consequently, the OEP procedure needs to be interpreted as a local-potential approximation but not to an actually existing local-exchange potential defined by a functional derivative. In the force-based approach that avoids any reference to functional differentiability, a different virial relation is derived. We start by applying the Helmholtz decomposition

^{52}to the exchange vector field. This yields a longitudinal (curl-free) and a transverse (divergence-free) vector field component as follows:

^{28}give the same relation, just without spin sum)

*α*_{fx}, and Table I). Hence,

*v*

_{fx}satisfies a virial relation of the form of Eq. (22) for closed-shell spherically symmetric systems, but in general we have the more involved Eq. (24) including a transverse component through the curl term. This is due to the fact that the exchange vector field

*f*_{x}is

*not purely longitudinal*, and hence while the exchange energy is directly linked to the exchange force density, the longitudinal part of the exchange vector field alone cannot yield the full exchange energy in general. Since OEP methods do not fulfill the virial relation of the form of Eq. (22) even in the spherically symmetric case (see Table I), this implies that the local-exchange potential from the force-balance approach is in general different from an exchange potential defined as a (generalized) exchange-energy derivative,

^{15}like those obtained by common OEP procedures. This was already pointed out in the work of Wang

*et al.*

^{53}when discussing the local-exchange potential of Harbola and Sahni,

^{27}which is equivalent to Eq. (19). To show this, they derived the second-order gradient expansion of both, the gradient of the energy-based exchange potential and

*f*_{x}=

**F**

_{x}[Φ]/

*ρ*from Eq. (14), and showed that the expressions do not match. However, note that in order to make this a strict statement about the potentials, we need to assume that

*f*_{x}is a gradient field, i.e.,

*α*_{fx}= 0, which does not hold in general.

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

Li | 245.3 | −0.049 | 0.647 | −1.464 |

Be | 415.2 | 0.037 | 32.87 | 17.97 |

Ne | 208.2 | −0.001 | 30.94 | 33.139 |

Na | 896.1 | 0.080 | −22.08 | −36.4 |

Mg | 1328.1 | 0.353 | 60.40 | −19.49 |

Ar | 221.95 | 0.000 | 7.61 | 8.21 |

Ca | 603.4 | −0.011 | 17.51 | −1.52 |

Zn | 6225.2 | 0.82 | 26.11 | −83.34 |

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

Li | 245.3 | −0.049 | 0.647 | −1.464 |

Be | 415.2 | 0.037 | 32.87 | 17.97 |

Ne | 208.2 | −0.001 | 30.94 | 33.139 |

Na | 896.1 | 0.080 | −22.08 | −36.4 |

Mg | 1328.1 | 0.353 | 60.40 | −19.49 |

Ar | 221.95 | 0.000 | 7.61 | 8.21 |

Ca | 603.4 | −0.011 | 17.51 | −1.52 |

Zn | 6225.2 | 0.82 | 26.11 | −83.34 |

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.^{54} 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.

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

Li | −29.34 | 1.496 | 1.234 | 0.787 |

Be | −39.33 | −2.255 | 0.040 | 0.909 |

Ne | −27.51 | −7.411 | −1.981 | 2.505 |

Na | −98.07 | 2.112 | 2.756 | 4.400 |

Mg | −118.2 | −2.106 | −0.467 | 4.099 |

Ar | −22.34 | 2.417 | 0.856 | 0.081 |

Ca | −91.12 | 0.113 | 1.903 | 1.508 |

Zn | −365.7 | −81.22 | 56.42 | 9.788 |

MARE (%) | 1.49 | 0.116 | 0.077 | 0.035 |

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

Li | −29.34 | 1.496 | 1.234 | 0.787 |

Be | −39.33 | −2.255 | 0.040 | 0.909 |

Ne | −27.51 | −7.411 | −1.981 | 2.505 |

Na | −98.07 | 2.112 | 2.756 | 4.400 |

Mg | −118.2 | −2.106 | −0.467 | 4.099 |

Ar | −22.34 | 2.417 | 0.856 | 0.081 |

Ca | −91.12 | 0.113 | 1.903 | 1.508 |

Zn | −365.7 | −81.22 | 56.42 | 9.788 |

MARE (%) | 1.49 | 0.116 | 0.077 | 0.035 |

*f*_{Hxc}is purely longitudinal and since the same holds by construction for the Hartree part, we must have a zero transverse contribution in

*f*_{x}+

*f*_{c}. If we now define

*α*_{fc}analogous to

*α*_{fx}, then this means that at each point in space and for every spin component it must hold that

Finally, let us comment on the homogeneous-density limit. In the work of Tchenkoue *et al.*,^{26} it was demonstrated how the usual Slater X*α*^{55} and local-density approximation (LDA)^{44} 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.*^{53} This directly connects the most fundamental functional approximations of DFT with the present formalism.

## V. THE FORCE-BASED APPROACH IN OTHER DFT VARIANTS

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.^{26} We even find that *v*_{fx} and *α*_{fx} can be chosen to be the local-exchange potential of CDFT and of time-dependent CDFT.^{26} This makes *v*_{fx} nicely compatible with this variant of DFT.

*∂*

_{t}

*ρ*(

**r**

*σ*,

*t*) = −∇·

**j**(

**r**

*σ*,

*t*) = −∇·

**j**

_{s}(

**r**

*σ*,

*t*) for both systems that share the density

*ρ*(

**r**

*σ*,

*t*) at all considered times. This is how we arrive at

^{5}

^{,}

*v*

_{fx}, yet the difference can be determined from

*α*_{fx}.

^{40}On the other hand, if

*f*_{Hxc}(

**r**

*σ*,

*t*) would be purely longitudinal, then Eq. (8) also holds in the time-dependent case and then Eq. (27) is a direct consequence of it by just multiplying with

*ρ*(

**r**

*σ*,

*t*) and taking the divergence. Conversely, it is only the transverse part of

*f*_{Hxc}(

**r**

*σ*,

*t*) that makes the difference when we compare

*instantaneously*the time-dependent case of Eq. (27) and the static case of Eq. (8). In other words, if we only consider the wave functions/forces at a given instant, it is only the nonzero phases/transverse forces that inform us whether we are considering a time-dependent situation. While we do not have access in this

*instantaneous*picture to all memory effects,

^{56,57}we nevertheless see that the transverse forces are important to generate memory over time. In an exchange-only approximation, this role is then taken over by

*α*_{fx}.

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.

## VI. NUMERICAL TESTS

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^{58} and ran simulations for a set of atoms in closed-shell configurations using norm-conserving pseudopotentials,^{59} 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^{60} (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.^{61} 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.

## VII. OUTLOOK AND CONCLUSIONS

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.^{26} 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,^{66} 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.^{40} 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,^{67} derive expressions for these terms for specific cases (e.g., the homogeneous limit) from wave function methods potentially augmented by modern machine-learning techniques,^{68} 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: NUMERICALLY CONVENIENT FORMS OF THE LOCAL-EXCHANGE POTENTIAL

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, $\rho ij*(r\u2032\sigma )/\rho (r\u2032\sigma )$ is uniformly equal to 1 and that then the second term vanishes.

### APPENDIX B: SCALING BEHAVIOR AND VIRIAL RELATION

*v*

_{fx,λ}(

**r**

*σ*) =

*λv*

_{fx}((

*λ*

**r**)

*σ*), where

*v*

_{fx,λ}is the expression from Eq. (19) with

*ρ*↦

*ρ*

_{λ}and $\rho (1)\u21a6\rho \lambda (1)$ replaced. Similarly, one finds

*E*

_{x}[

*ρ*

_{λ}] =

*λE*

_{x}[

*ρ*] directly from Eq. (21). This is the correct scaling behavior for the exchange energy.

^{44}Together with the assumption of functional differentiability of

*E*

_{x}[

*ρ*] as a

*density*functional and applicability of the usual chain rule, this suffices to derive the virial relation of Eq. (22). By virtue of the chain rule of functional calculus, we have

**r**= 3. It needs to be stressed that this form of the virial relation depends on the assumption of functional differentiability and that

*v*

_{x}(

**r**

*σ*) =

*δE*

_{x}[

*ρ*]/

*δρ*(

**r**

*σ*) defines a

*different*local-exchange potential than

*v*

_{fx}(

**r**

*σ*) given by a force-based approach in Eq. (19) or by OEPx.

*E*

_{x}[Φ] and

**F**

_{x}[Φ] by direct computation. Starting with the exchange-energy expression Eq. (21), we use the identity (

**r**−

**r**′)·∇|

**r**−

**r**′|

^{α}=

*α*|

**r**−

**r**′|

^{α}(which is also central for deriving the usual

*virial theorem*) with

*α*= −1 and the symmetry of the whole expression in

**r**↔

**r**′.

*v*

_{fx}, in Eq. (19) we first switch ∇′ over to the term 1/(4

*π*|

**r**−

**r**′|) by partial integration and then switch ∇′ → −∇ by symmetry. We have

**A**) = Δ

**A**+ ∇ × (∇ ×

**A**) and ∇ ×

*f*(

**r**)

**C**= (∇

*f*(

**r**)) ×

**C**were used. Now, the first part gives exactly

*E*

_{x}according to Eq. (B3) while the second line appears as an additional term in a virial relation between

*E*

_{x}and

*v*

_{fx}. However, since it appears as the curl of a vector expression it cannot be equal to the gradient of a scalar potential, so the difference comes from the transverse part of

*f*_{x}while

*v*

_{fx}corresponds only to the longitudinal part of

*f*_{x}. The nice thing is that this gives an explicit form for the transverse part of

*f*_{x}, while the longitudinal part is already given by −∇

*v*

_{fx}. We thus find the following Helmholtz decomposition:

*E*

_{x}[Φ] and

*v*

_{fx}holds in the form of Eq. (B2). We show that for spherically symmetric densities

*ρ*(

**r**

*σ*) =

*R*

_{σ}(|

**r**|), this is indeed the case. For this, we take the last integral of Eq. (B9) and perform integration by parts with the curl and vanishing boundary terms to get

**r**= 0 and $(\u2207\rho (r\sigma ))\xd7r=(r\xd7r)R\sigma \u2032(|r|)/|r|=0$, so the above expression evaluates as zero.

### APPENDIX C: NUMERICAL RESULTS FOR THE FBEx FUNCTIONAL

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.

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . | HF . |
---|---|---|---|---|---|

Li | 0.101 | 0.086 | 0.082 | 0.082 | 0.082 |

Be | 0.325 | 0.311 | 0.307 | 0.307 | 0.307 |

Ne | 0.900 | 0.859 | 0.843 | 0.845 | 0.844 |

Na | 0.118 | 0.083 | 0.074 | 0.074 | 0.074 |

Mg | 0.285 | 0.260 | 0.253 | 0.253 | 0.253 |

Ar | 0.622 | 0.585 | 0.590 | 0.590 | 0.590 |

Ca | 0.224 | 0.201 | 0.195 | 0.195 | 0.195 |

Zn | 0.368 | 0.332 | 0.300 | 0.300 | 0.300 |

Atom . | Slater . | FBEx . | OEPx-KLI . | OEPx . | HF . |
---|---|---|---|---|---|

Li | 0.101 | 0.086 | 0.082 | 0.082 | 0.082 |

Be | 0.325 | 0.311 | 0.307 | 0.307 | 0.307 |

Ne | 0.900 | 0.859 | 0.843 | 0.845 | 0.844 |

Na | 0.118 | 0.083 | 0.074 | 0.074 | 0.074 |

Mg | 0.285 | 0.260 | 0.253 | 0.253 | 0.253 |

Ar | 0.622 | 0.585 | 0.590 | 0.590 | 0.590 |

Ca | 0.224 | 0.201 | 0.195 | 0.195 | 0.195 |

Zn | 0.368 | 0.332 | 0.300 | 0.300 | 0.300 |

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_{2}, we employed a N–N distance of $1.09769A\u030a$. For CO_{2}, we considered a C–O bond length of $1.16A\u030a$. For CH_{4}, we considered a C–H bond length of $1.087A\u030a$. 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.

Molecule . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

N_{2} | 93.34 | −134.7 | −276.4 | −235.4 |

CO_{2} | 13.64 | −520.9 | −1157 | −660.1 |

CH_{4} | 72.36 | −12.517 | −38.66 | −19.15 |

Molecule . | Slater . | FBEx . | OEPx-KLI . | OEPx . |
---|---|---|---|---|

N_{2} | 93.34 | −134.7 | −276.4 | −235.4 |

CO_{2} | 13.64 | −520.9 | −1157 | −660.1 |

CH_{4} | 72.36 | −12.517 | −38.66 | −19.15 |

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.

Molecule . | Slater . | FBEx . | OEPx-KLI . | OEPx . | HF . |
---|---|---|---|---|---|

N_{2} | 0.635 | 0.607 | 0.629 | 0.630 | 0.617 |

CO_{2} | 0.619 | 0.586 | 0.545 | 0.544 | 0.546 |

CH_{4} | 0.566 | 0.540 | 0.543 | 0.545 | 0.546 |

Molecule . | Slater . | FBEx . | OEPx-KLI . | OEPx . | HF . |
---|---|---|---|---|---|

N_{2} | 0.635 | 0.607 | 0.629 | 0.630 | 0.617 |

CO_{2} | 0.619 | 0.586 | 0.545 | 0.544 | 0.546 |

CH_{4} | 0.566 | 0.540 | 0.543 | 0.545 | 0.546 |

### APPENDIX D: COMPARISON WITH THE EXACT EXCHANGE–CORRELATION POTENTIAL

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,^{69} or from CI densities.^{70} 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.

## REFERENCES

*Density Functional Theory: An Approach to the Quantum Many-Body Problem*

*Density Functional Theory: An Advanced Course*

*Time-Dependent Density-Functional Theory: Concepts and Applications*

*Fundamentals of Time-Dependent Density Functional Theory*

*Mathematical Methods in Physics: Distributions, Hilbert Space Operators, Variational Methods, and Applications in Quantum Physics*

*Nonequilibrium Many-Body Theory of Quantum Systems: A Modern Introduction*

*Density-Functional Theory of Atoms and Molecules*

*Mathematical Methods for Physicists*

_{3}Bi