# Modeling the Kohn–Sham potential for molecular dissociation with orbital-independent functionals: A proof of principle

Sara Giarrusso, Federica Agostini

► To cite this version:

Sara Giarrusso, Federica Agostini. Modeling the Kohn–Sham potential for molecular dissociation with orbital-independent functionals: A proof of principle. The Journal of Chemical Physics, 2025, 162 (9), 10.1063/5.0243705. [https://arxiv.org/abs/2009.09486](https://arxiv.org/abs/2009.09486)

HAL Id: hal-04988564

https://hal.science/hal-04988564v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling the Kohn–Sham potential for molecular dissociation with orbital-independent functionals: A proof of principle

Sara Giarrusso ; Federica Agostini

<figure><img src="image_2.png" /></figure>

Check for updates

J. Chem. Phys. 162, 091103 (2025)

https://doi.org/10.1063/5.0243705

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

## Articles You May Be Interested In

KoopmanLab: Machine learning for solving complex physics equations

APL Mach. Learn. (September 2023)

Experimental realization of a quantum classification: Bell state measurement via machine learning

APL Mach. Learn. (September 2023)

<figure><img src="image_5.png" /></figure>

Nanotechnology &
Materials Science

<figure><img src="image_6.png" /></figure>

Optics &
Photonics

<figure><img src="image_7.png" /></figure>

Impedance
Analysis

<figure><img src="image_8.png" /></figure>

Scanning Probe
Microscopy

<figure><img src="image_9.png" /></figure>

Sensors

<figure><img src="image_10.png" /></figure>

Failure Analysis &
Semiconductors

<figure><img src="image_11.png" /></figure>

Unlock the Full Spectrum.
From DC to 8.5 GHz.

Your Application. Measured.

Find out more

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

Zurich
Instruments

<figure><img src="image_14.png" /></figure>

AIP
Publishing
# Modeling the Kohn-Sham potential for molecular dissociation with orbital-independent functionals: A proof of principle

Cite as: J. Chem. Phys. 162, 091103 (2025); doi: 10.1063/5.0243705

Submitted: 13 October 2024 • Accepted: 17 February 2025 •

Published Online: 6 March 2025

<figure><img src="image_15.png" /></figure>

View Online

<figure><img src="image_16.png" /></figure>

Export Citation

<figure><img src="image_17.png" /></figure>

CrossMark

Sara Giarrusso<sup>a)</sup> and Federica Agostini

## AFFILIATIONS

Université Paris-Saclay, CNRS, Institut de Chimie Physique UMR8000, 91405 Orsay, France

<sup>a)</sup>Author to whom correspondence should be addressed: sara.giarrusso@universite-paris-saclay.fr

## ABSTRACT

We model the Hartree-exchange-correlation potential of Kohn-Sham (KS) density-functional theory adopting a novel strategy inspired by the strictly correlated-electron limit and relying on the exact decomposition of the potential based on the exact factorization formalism. Starting with accurate density and conditional potential for a one-dimensional model of a stretched heteronuclear molecule, we provide a proof-of-principle example of an approximation that accurately reproduces the step of the exact KS potential without resorting to the KS orbitals, virtual or occupied. We also test our strategy using the approximate conditional potentials corresponding to the strictly correlated-electron and the exact-exchange functionals. The results are encouraging in that the initial approximations are modified in the qualitatively correct way: decreasing correlation in the former and increasing it in the latter.

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0243705

Kohn-Sham (KS) density functional theory (DFT) is an elegant and computationally efficient theory for electronic structures whose great success hinges on the mapping of a physical system of interacting electrons onto an auxiliary, formally non-interacting one, i.e., the KS system.<sup>1,2</sup> The KS external potential comprises the external potential of the original system and the Hartree-exchange-correlation (Hartree-XC) potential, whose XC part has to be approximated. Accurately capturing chemical and physical properties relies on the quality of the approximations of the XC potential. Nonetheless, fundamental processes in chemistry, such as bond breaking or charge localization, or properties, such as electric polarizabilities and bandgaps, can be captured correctly only if the XC potential manifests *non-intuitive* features, such as “peaks” and “steps.”<sup>3–25</sup>

Peaks and steps are necessary to enforce all the correlation effects of the original system in a non-interacting system. For instance, without the step structure around the more electronegative atom in the KS potential for the stretched LiH molecule, there would be a spurious accumulation of electron density on the hydrogen atom (more electronegative), resulting in fractionally charged fragments. Such accumulation of charge is impeded by electronic

repulsion in the interacting system. Conversely, the step structure around H in the exact KS potential *de facto* aligns the ionization potentials of the two fragments in the non-interacting system, preventing spurious electron density accumulation and, therefore, unphysical fractional charges.<sup>3,6</sup>

Mathematically, peaks and steps emerge as a consequence of the non-local dependence of the KS potential on the electronic density. It is logical, then, that the numerically efficient local or semilocal approximate XC functionals fail in capturing these features. On the other hand, many orbital-dependent functionals require virtual KS orbitals to incorporate non-locality,<sup>26–33</sup> which makes them expensive, thus limiting their applicability.

In this work, we propose a novel strategy aiming to model correctly the step structure of the exact potential in the molecular dissociation limit. We work within an orbital-independent theoretical framework, which allows us to approximate accurately the Hartree-XC potential based on a conceptually simple numerical procedure. The first key ingredient of this strategy is a decomposition of the Hartree-XC potential into local potentials, called kinetic, N − 1, and conditional potential.<sup>4</sup> Such decomposition is based on the theory of the conditional amplitude,<sup>34,35</sup> which, in the context
of the electronic wavefunction, has been more recently also referred to as exact electron factorization (EEF).<sup>17,23,24,36</sup> The great merit of this decomposition is its ability to isolate the contributions of the Hartree-XC potential responsible for peaks and steps, namely, the kinetic and N − 1 potentials, respectively.<sup>4–9,19,21,22</sup> Our idea is to approximate each component of the Hartree-XC potential individually, as it arises from the EEF formalism. Since this formalism is exact, accurately reproducing these components would yield a total Hartree-XC potential that behaves as the functional derivative of the Hartree-XC energy. This provides a rigorous foundation for directly modeling the potential. Another advantage of this framework is its ability to clearly distinguish between components that explicitly contribute to the energy (the conditional and kinetic potentials) and those enforcing appropriate functional behavior with respect to density changes (the N − 1 potential).<sup>4,9,21,22</sup> This distinction eliminates ambiguity when deriving the Hartree-XC energy from the approximate potential.

The second key ingredient of our approach is the conditional density. In fact, we model the individual potentials as explicit functionals of this quantity. To approximate the conditional density, we take inspiration from the mathematical structure of the strong-interaction or strictly correlated-electron (SCE) limit of DFT.<sup>37–40</sup> Throughout this paper, we adopt a model that is paradigmatic of a stretched heteronuclear molecule and for which accurate data are available.<sup>22</sup> Our strategy is tested at two different levels, both of which build an approximate conditional density: starting, first, from the exact conditional potential and, then, from approximate conditional potentials, coming from the SCE and exact-exchange (EXX) functionals.<sup>41,42</sup> Before illustrating our strategy, let us review the crucial definitions in the EEF framework.

The conditional amplitude, Φ( σ, 2 , . . . , N ; **r**), is

$$
\Phi(\sigma, 2, \dots, N; \mathbf{r}) := \frac{\Psi(1, \dots, N) \sqrt{N}}{\sqrt{n(\mathbf{r})}}, \quad (1)
$$

where Ψ and n are, respectively, the ground state (GS) wavefunc-
tion and GS density of the N-electron Hamiltonian under study and
i = σᵢ**rᵢ indicates the spatial and spin coordinates of the ith elec-
tron. The subscript 1 from the coordinates of electron 1 chosen as
reference is dropped, i.e., σ₁ = σ and **r₁ = r**.

According to the work of Buijse and collaborators,⁴ the
Hartree-XC potential, v<sub>Hxc</sub>(**r**), is exactly decomposed as

$$
v_{\text{Hxc}}(\mathbf{r}) = v_{\text{kin}}(\mathbf{r}) - v_{\text{s,kin}}(\mathbf{r}) + v_{\text{N-1}}(\mathbf{r}) - v_{\text{s,N-1}}(\mathbf{r}) + v_{\text{cond}}(\mathbf{r}) \quad (2)
$$

in terms of the kinetic potential v<sub>kin</sub>(**r**): = 1/2 ∫ |∇**r**| Φ|<sup>2</sup> dσd2···dN,
the N − 1 potential v<sub>N − 1</sub>(**r**): = ∫ Φ* ˆH<sup>N − 1</sup>Φ dσd2···dN − E<sub>0</sub><sup>N − 1</sup>,
and the conditional potential v<sub>cond</sub>(**r**): = ∫∑<sub>i=2</sub><sup>N</sup> v<sub>ee</sub>(**r** − **r**<sub>i</sub>)
|Φ|<sup>2</sup> dσd2···dN. We used the symbols ˆH<sup>N − 1</sup> and E<sub>0</sub><sup>N − 1</sup>, respectively,
for the Hamiltonian and the GS energy of the ionized system and
v<sub>ee</sub> for any interaction function. Finally, v<sub>s,kin</sub> and v<sub>s,N − 1</sub> for the KS
system are formally defined as for the interacting system, replacing
Φ with the KS conditional amplitude obtained from (1) when
substituting Ψ with the KS wavefunction. Note that, in the case of
two-electron singlets, the KS kinetic and N − 1 potentials are zero.<sup>22</sup>

In this work, we focus on the approximation of v<sub>kin</sub> and v<sub>N-1</sub>
in the case of two-electron singlets. Two-electron singlets are an
excellent proxy for studying the features of the KS potential in the

breaking of a chemical bond. This is essentially because the peaks
and the step structure arising specifically in dissociation processes
are located in very low-density regions, where the overall density
is dominated by the density of the sole highest-occupied molecular
orbital (HOMO), occupied by the two electrons involved in the bond
that is breaking. Therefore, we can regard these features as coming
from a two-electron singlet system where only v<sub>kin</sub> and v<sub>N-1</sub> play
a role. In turn, these are also the components of the total potential
that have proven the most delicate to model since they depend on
the interacting conditional amplitude, as opposed to v<sub>s,kin</sub> and v<sub>s,N-1</sub>
that can be written explicitly in terms of the occupied KS orbitals.

For two-electron singlets, the conditional density is simply
given by the modulus squared of the conditional amplitude,

$$
n_{\text{cond}}(\mathbf{r}, \mathbf{r}') = |\Phi(\mathbf{r}'; \mathbf{r})|^2, \quad (3)
$$

and, as a result, one can conveniently rewrite in terms of n<sub>cond</sub> the
kinetic potential as

$$
v_{\text{kin}}(\mathbf{r}) = \frac{1}{8} \int \frac{(\nabla_r n_{\text{cond}}(\mathbf{r}, \mathbf{r}'))^2}{n_{\text{cond}}(\mathbf{r}, \mathbf{r}')} d\mathbf{r}' \quad (4)
$$

and the N − 1 potential as

$$
v_{N-1}(\mathbf{r}) = \int \left( \frac{(\nabla_{\mathbf{r}'} n_{\text{cond}}(\mathbf{r}, \mathbf{r}'))^2}{8 n_{\text{cond}}(\mathbf{r}, \mathbf{r}')} + v(\mathbf{r}') n_{\text{cond}}(\mathbf{r}, \mathbf{r}') \right) d\mathbf{r}' - E_0^{N-1}, \quad (5)
$$

where v is the external potential of the problem under study. Impor-
tantly, the conditional potential can always be written in terms of the
conditional density,

$$
v_{\text{cond}}(\mathbf{r}) = \int n_{\text{cond}}(\mathbf{r}, \mathbf{r}') v_{\text{ee}}(\mathbf{r} - \mathbf{r}') d\mathbf{r}'. \quad (6)
$$

In conclusion, we can determine each individual contribution to the
total potential once we have a model for the conditional density.

To approximate this latter quantity, we propose the following
*ansatz*:

$$
n_{\text{cond}}^{f_i}(\mathbf{r}, \mathbf{r}') = n(\mathbf{r}') \sum_{i=1}^{N-1} G(\mathbf{r}' - f_i(\mathbf{r})), \quad (7)
$$

where fᵢ are special functions to be discussed and G(r) = e⁻β r² is a Gaussian function. Its parameter β, related to the width of the Gaussian, is to be determined from the condition ∫ n<sub>cond</sub>(r, r')dr' = 1∀r, which ensures that the probability of finding any of the N − 1 electrons over all space is one regardless of where the reference electron is located. Note that expression (7) is valid for a general number of electrons N, although in the examples we consider for testing, we limit ourselves to the case N = 2.

The *ansatz* of Eq. (7) is inspired by the structure of the SCE conditional density,<sup>40,43</sup>

$$
n_{\text{cond}}^{\text{SCE}}(\mathbf{r}, \mathbf{r}') = \sum_{i=2}^{N} \delta(\mathbf{r}' - f_{i,\text{SCE}}(\mathbf{r})), \quad (8)
$$

where f<sub>i,SCE</sub> are the infamous co-motion functions—or just co-
motions—of the SCE limit. They give the classical positions of the
N − 1 electrons as a function of the position of a reference elec-
tron. Two formal differences distinguish the SCE conditional density
[Eq. (8)] from our ansatz in Eq. (7), namely, (i) the Dirac delta has been fattened into a Gaussian function and (ii) the Gaussian is furthermore weighted by $n(r')$. Both these changes reduce the excessive correlation among electrons present in the SCE description. Correction (i) allows for quantum fluctuations around the parametric curve given by the co-motions. Such an approach has appeared in different studies that moved from the leading order in the strong-interaction limit to the next.39,44,45 Second, via correction (ii), the probability distribution of any of the $N-1$ electrons does not depend solely on the position $r$ of the reference electron but also on the probability of having any other electron in $r'$. Consequently, the term $\sum_{i=1}^{N-1} G(r' - f_i(r))$ is recognized as an approximate pair-correlation function.

In addition to (i) and (ii), the shape of the co-motions itself is another factor that greatly influences the structure of the approximate conditional density. The concept of the co-motions has been introduced in electronic-structure methods by the SCE framework and has remained relatively limited to this extreme scenario. Here, we explore the possibility of determining effective co-motions that can give a good description of the relative position of electrons at the physical regime. To distinguish them from the other effective models that will be proposed, we label the usual co-motions, associated with the SCE limit, as $f_{i,SCE}$. The symbol $f_i$ in Eq. (7) encompasses any function that plays the role of a co-motion: determining the position of any of the $N-1$ electrons as a function of the reference one.

Our approach is somewhat reminiscent of the construction of the multi-radii functional (MRF).46,47 There, one constructs an approximate spherically averaged conditional density from effective radial distances among electrons. Here, instead, we consider not only the radial distances but also effective co-motions where the information on the geometric orientation of electron pairs within the molecular structure is retained, aiming to model the full Hartree-XC potential, not (only) the energy density as in MRF.

To test our ansatz, we use input data from Ref. 22. The Hamiltonian considered there reads

$$ \hat{H} = \sum_{i=1,2} \left( -\frac{1}{2} \frac{\partial^2}{\partial x_i^2} + v(x_i) \right) + \frac{1}{\sqrt{c + (x_1 - x_2)^2}}, \quad (9) $$

with $v(x) = -\frac{1}{\sqrt{a+(x-\frac{R}{2})^2}} - \frac{1}{\sqrt{b+(x+\frac{R}{2})^2}}$. The repulsive interaction potential and the external potential are soft-Coulomb interaction functions where the Coulomb singularity in 1D is softened by a parameter at the denominator. The external potential, $v$, features two wells representing the (different) nuclei, separated by a distance $R$. Nucleus A is located at $x = \frac{R}{2}$, with potential value $1/\sqrt{a}$, and nucleus B is at $x = -\frac{R}{2}$, with potential value $1/\sqrt{b}$. The softening parameters have values $a = 0.7 a_0$, $b = 2.25 a_0$, and $c = 0.6 a_0$. With such parameters, the ionization potential of the isolated fragments, $\Delta I$, is $0.3 E_h$, reproducing $\Delta I$ of the real LiH molecule.

The exact conditional amplitude, $\Phi(x_2; x_1)$, and the corresponding kinetic, $N-1$, and conditional potentials have been calculated from the wavefunction [Eq. (1)] by solving the Schrödinger equation on a grid at various internuclear distances.22

The procedure to approximate the components of the Hartree-XC potential, including the peaks and the steps features, is as follows:

1. defining an effective distance, $D(x)$, between the two particles as a function of the position of the reference electron;

2. determining the co-motions in terms of $D(x)$;

3. constructing the conditional density according to Eq. (7) by using the co-motions obtained at point 2; and

4. using the conditional density to calculate the kinetic [Eq. (4)], $N-1$ [Eq. (5)], and conditional [Eq. (6)] potentials.

To initiate our procedure, thus for point 1, we take as a measure of the distance between the electron pair the inverse of the conditional potential, as done in the MRF approximation,

$$ D(x) = \frac{1}{v_{\text{cond}}(x)}, \quad (10) $$

essentially regarding the two particles as point charges. In Eq. (10), in principle, any $v_{\text{cond}}$ can be inserted since it is only used at this stage of the procedure to determine the co-motions and, in turn, the conditional density. We begin by using the exact conditional potential, and later, we shall consider approximate conditional potentials. Note that, because of the soft-Coulomb interaction in our model, the distance is more precisely $D(x) = \sqrt{\frac{1-c}{v_{\text{cond}}(x)^2}}$. However, this is virtually indistinguishable from Eq. (10) for large $R$ (we take $R = 11 a_0$ as the representative of large $R$).

<figure><img src="image_18.png" /><figcaption>FIG. 1. Upper panel: classical distance between particles $D(x)$ [Eq. (10)] as a function of the position of the reference electron $x$. The functions $|x|$ and $v_{\text{cond}}(x)$ are also plotted for comparison. This latter is enhanced by a factor 500, i.e., in units of $2 \times 10^{-3} E_h$. Lower panel: derivative of the distance between particles, $D'(x)$. The two transition regions where the derivative is not constant correspond to the two maxima of the conditional potential.</figcaption></figure>
In Fig. 1 (upper panel), we plot $D(x)$ of Eq. (10) using the exact conditional potential, together with $|x|$ (for reference) and the conditional potential. $D(x)$ is decreasing as $-x$ (modulo a constant) roughly up to the global maximum of the conditional potential, past which it starts increasing as $x$. The constant shift of $D(x)$ relative to $|x|$ is not artificial and simply follows from our definition of the distance [Eq. (10)]. In other words, this shift is “hidden” in the exact conditional potential. The switch from above to below $x$ takes place around $31 a_0$, where the repulsion between electrons is locally increasing. This overall behavior is fully appreciated in the derivative $D'(x)$, shown in the lower panel.

The local maximum (shoulder) in the conditional potential in addition to the global one around the bond midpoint was first identified in Ref. 22. Both maxima are related to positions of the reference electron for which the conditional amplitude changes from resembling the ground state of the ion to resembling its first excited state or vice versa. These positions also relate to the two peaks of the kinetic potential $v_{\text{kin}}$ and to the extent of the step structure in the $N-1$ potential $v_{N-1}$. The reader is referred to Ref. 22 for an extensive discussion of these features of the exact potentials.

We move to point 2, and we determine an effective co-motion to use in Eq. (7), for which we use the symbol $f_{\text{cond}}(x)$ to indicate that it is the one obtained by using the exact $v_{\text{cond}}$ in Eq. (10). Recalling that the co-motion gives the classical position of an electron as a function of the reference one in $x$, and thanks to the reduced dimensionality of the model, the distance between the electron pair is $D(x) = |x - f_{\text{cond}}(x)|$, or

$$f_{\text{cond}}(x) = \begin{cases} x + D(x), & x \le x_{\text{crit}}, \\ x - D(x), & x > x_{\text{crit}}, \end{cases} \quad (11)$$

where $x_{\text{crit}}$ is a critical value of $x$ to be assigned. Note first that the function $f_{\text{cond}}$ is inherently discontinuous at $x = x_{\text{crit}}$ since $D(x)$ is strictly positive. To ensure that it is as well-behaved as possible, we then select $x_{\text{crit}}$ such that its derivative, $f'_{\text{cond}}(x)$, is continuous. This corresponds to determining $x_{\text{crit}}$ from $D'(x)|_{x=x_{\text{crit}}} = 0$.

Figure 2 shows the co-motion $f_{\text{cond}}(x)$ obtained as described above. The physical interpretation of this co-motion is quite simple and reasonable, yielding the position of the remaining electron as a function of the position $x$ of the reference electron. In Fig. 2, we label $\Omega_A$ the region where the reference electron feels only the

<figure><img src="image_19.png" /><figcaption>FIG. 2. Effective co-motion, $f_{\text{cond}}(x)$, giving the position of the remaining electron as a function of the position of the reference one for internuclear distance $R = 11 a_0$.</figcaption></figure>

field of nucleus A and vice versa for $\Omega_B$.²² In $\Omega_A$ (red shaded area), when the reference electron is either sitting on nucleus A or on its right, up to $\approx 20 a_0$, the position of the remaining electron is constant and negative, $\approx -5.3 a_0$; thus, it is located around nucleus B. In turn, when the reference electron is on nucleus B or on its left, the position of the remaining electron is constant and positive, $\approx 5.4 a_0$; thus, it is located around nucleus B, as indicated in Fig. 2 in $\Omega_B$ (blue shaded areas). When the reference electron is being ionized and is located at large positive values of $x$, the remaining electron is close to nucleus A; thus, the reference electron effectively feels the field of nucleus B (therefore, the region on the right is labeled $\Omega_B$). The transition regions (in purple) correspond to regions where the reference electron feels both nuclei.

As indicated in point 3, we then construct an approximate conditional density via Eq. (7). The parameter of the Gaussian, $\beta^{f_{\text{cond}}}$, associated with the co-motion $f_{\text{cond}}$, is determined numerically from the normalization condition for the conditional density. The function $\beta^{f_{\text{cond}}}(x)$ is plotted in Fig. (SM1) of the supplementary material, together with the parameters $\beta^{\text{fsce}}(x)$ and $\beta^{\text{fex}}(x)$, which are associated with different kinds of co-motions (SCE and EXX) and will be illustrated later.

Figure 3 shows the first main result of this work: we show our approximate conditional density $n_{\text{cond}}^{f_{\text{cond}}}(x_1, x_2)$, obtained from the ansatz (7), and we compare it with the exact one. In the left column of the figure, where the conditional densities are plotted as 3D surfaces, the accuracy of the proposed ansatz stands out. The only region of appreciable deviation is around $x_2 \approx -\frac{R}{2}$ and $x_1 \in \Omega_A$, where there is an excess of probability density (as well as around $x_2 \approx \frac{R}{2}$ and $x_1 \in \Omega_B$, less visible in the figure). Plotting the natural logarithm of both conditional densities stresses other structural differences between the two, notably their different decay along the diagonal $x_1 = x_2$.

Finally, according to point 4, we can construct the local potentials of Eqs. (4)–(6) to assess how these deviations in the approximate conditional density impact the potentials. The resulting potentials are displayed individually in Fig. (SM2) of the supplementary material. In Fig. 4, we show their sum ($v_{\text{cond}} + v_{\text{kin}} + v_{N-1}$), which corresponds to an approximate Hartree-XC potential in two-electron singlets. In the uppermost panel, we report the results for $R = 11 a_0$. The accuracy of the model proposed is quite remarkable, particularly in the step structure, which has been modeled satisfactorily without any resort to virtual KS orbitals. To the best of our knowledge, previous models of this structure insisted on the necessity for using virtual KS orbital(s) as input.²⁷,³² Here, instead, we demonstrate that it is possible to reproduce this peculiar feature of the potential in the limit of a fully dissociated molecule within an orbital-independent framework. Concerning the peak around the bond midpoint, the lack of this feature in the approximate potential is to be attributed to the discontinuity of the co-motion. This problem seems specific to the one-dimensional setting and is analyzed in some detail in Sec. SM2 of the supplementary material.

To further assess the quality of our ansatz for the conditional density [Eq. (7)], we have applied the same strategy (points 1 to 4) on a geometry close to equilibrium ($R = 3 a_0$) and a midway geometry ($R = 7 a_0$). The resulting Hartree-XC potentials are shown in the lowest and middle panels of Fig. 4 and demonstrate that our strategy adapts seamlessly to the different geometries. The results
<figure><img src="image_20.png" /><figcaption>FIG. 3. Upper panels: exact conditional density $n_{\text{cond}}(x_1, x_2)$ (left) and contour plot of its logarithm (right). Lower panels: approximate conditional density $n_{\text{cond}}^{f_{\text{cond}}}(x_1, x_2)$ (left) and its logarithm (right). $R = 11$ $a_0$ everywhere.</figcaption></figure>

summarized in Fig. 4 represent the second main contribution of this work: the step structure at large interatomic distances can be modeled without invoking virtual orbitals using an approach that remains accurate across significantly varied molecular geometries and, therefore, across different correlation regimes.

Our results thus far serve as a proof of principle. However, the conditional potential used to initiate the procedure was the exact one, while in practice this ingredient is unavailable. To begin investigating the broader applicability of this strategy, we now test different approaches to constructing effective co-motion functions, starting from *approximate* forms of $v_{\text{cond}}$, and focusing on the dissociation limit.

First, let us use the SCE co-motion, $f_{\text{SCE}}$. To avoid the logarithmic divergence of $f_{\text{SCE}}$ at the point where the density integrates to one [which is $x_{\text{crit}}$ of Eq. (11) in the SCE co-motion], we adopt the model co-motion developed in Ref. 21 [see Sec. (SM3) of the *supplementary material* for further details]. The corresponding parameter, $\beta^{f_{\text{SCE}}}$, to construct our approximate conditional density, $n_{\text{cond}}^{f_{\text{SCE}}}$, is shown in blue in Fig. (SM1) of the *supplementary material*.

In Fig. 5 (upper left panel), the resulting conditional potential (blue) is shown in comparison with the exact (red) and the SCE (yellow) conditional potentials. Although it displays some incorrect

features and excessive repulsion, it appears qualitatively much closer to the shape of the exact potential, $v_{\text{cond}}$, than to the shape of $v_{\text{cond}}^{\text{SCE}}$. This result validates the strategy proposed here, aiming to modify the structure of the SCE conditional density [Eq. (8)] according to the *ansatz* (7) in order to capture more features of the electron-electron interaction in the physical regime. The total Hartree-XC potential coming from $n_{\text{cond}}^{f_{\text{SCE}}}$ (blue) is shown in Fig. 5 (lower left panel) in comparison with the exact (red) as well as the SCE (yellow) Hartree-XC potential.19,21 Also for the total potential, it is remarkable that the approximate potential from our *ansatz* more closely resembles the exact potential rather than the SCE potential, to which the specific co-motion function belongs.

Second, we apply our strategy to the exact-exchange approximation<sup>41,42</sup> (EXX).

As described in the *supplementary material*, the conditional potential and the total Hartree-XC potential in the EXX approximation for two-electron singlets coincide, i.e., $v_{\text{cond}}^{\text{EXX}} = v_{\text{Hxc}}^{\text{EXX}}$.

To construct an effective co-motion from $v_{\text{cond}}^{\text{EXX}}$, we first extract an effective distance as in Eq. (10). Then, we choose $x_{\text{crit}}$ such that the derivative of the resulting co-motion, $f_{\text{EXX}}$, is continuous, analogously to the construction of $f_{\text{cond}}$. The parameter, $\beta^{f_{\text{EXX}}}$, is shown
<figure><img src="image_21.png" /><figcaption>FIG. 4. Approximate Hartree-XC potential $v_{Hxc}^{\text{cond}}$ obtained from $\ell_{\text{cond}}$ and $n_{\text{cond}}^{\text{fexx}}$ (blue) in comparison with the exact one (red) for three different internuclear distances: $R = 11 a_0$ (uppermost panel), $7 a_0$ (middle), and $3 a_0$ (lowest).</figcaption></figure>

in yellow in Fig. (SM1) of the **supplementary material**. The resulting conditional potential (blue) is shown in Fig. 5 (upper right panel) together with the exact (red) and the EXX (yellow) conditional potentials. Although our approximation increases the repulsion around the positions of the nuclei if compared to the source potential ($v_{\text{cond}}^{\text{EXX}}$), it does introduce a qualitative correct feature: a small peak around the bond midpoint. This is a direct consequence of the structure of the proposed *ansatz*: the use of an effective co-motion introduces some correlation effects that are absent from the source potential. The same effect is visible in the total Hartree-XC potential, which now features a hint of step structure completely absent from the EXX potential. The EXX approximation alone cannot reproduce the step structure in two-electron singlet systems or regions dominated by the density of the HOMO because this feature arises purely from the correlation contribution in these regions.⁶ This limitation is related to the fact that the response component in the exchange potential for two-electron singlets is exactly zero.¹⁹

Unfortunately, neither the SCE-modified nor the EXX-modified potentials model the accurate extent of the step structure because their source potentials do not feature the local maximum of the exact conditional potential. Therefore, their effective co-motions do not model the transition around $31 a_0$ from one nucleus to the other, differently than $f_{\text{cond}}$ (Fig. 2). Furthermore, the EXX-modified potential (blue curve in the bottom right panel of Fig. 5) quantitatively worsens the already erroneous profile of the EXX potential, being way too repulsive on the nuclei. Therefore, we do not believe it is worthwhile to further explore the construction of an effective co-motion function for systems with a general number of electrons and geometries, starting from such an approximation. Nevertheless, it is instructive to note that even when starting from an approximation with exactly zero correlation, the mere adoption of a conditional density that structurally depends on a co-motion function can introduce some correlation.

In summary, in this work, we have presented a strategy to model non-local features of the exact Hartree-XC potential that

<figure><img src="image_22.png" /><figcaption>FIG. 5. Left column: Results obtained from $\ell_{\text{cond}}^{\text{fexx}}$ and $n_{\text{cond}}^{\text{fexx}}$. Upper panel: Approximate conditional potential coming from $n_{\text{cond}}^{\text{fexx}}$ (blue) in comparison with the exact one (red) and with the SCE conditional potential (yellow). Lower panel: Approximate Hartree-XC potential coming from $n_{\text{cond}}^{\text{fexx}}$ (blue) in comparison with the exact one (red) and with the SCE Hartree-XC potential (yellow). Right column: Results obtained from $\ell_{\text{EXX}}$ and $n_{\text{cond}}^{\text{fexx}}$. Upper panel: Approximate conditional potential coming from $n_{\text{cond}}^{\text{fexx}}$ (blue) in comparison with the exact one (red) and with the EXX conditional potential (yellow). Lower panel: Approximate Hartree-XC potential coming from $n_{\text{cond}}^{\text{fexx}}$ (blue) in comparison with the exact one (red) and with the EXX Hartree-XC potential (yellow). $R = 11 a_0$ everywhere.</figcaption></figure>
brings together the exact electron factorization and the strictly
correlated-electron limit of DFT. The features targeted are crucial
to describe the breaking of a chemical bond. We have provided a
proof-of-principle demonstration that it is possible to capture such
non-local features of the Hartree-XC potential, particularly its step
structure, without resorting to any KS virtual orbitals, which seems
to represent a shift compared to previous common understanding.
We have tested our strategy on a system that is paradigmatic of
the challenge faced by current KS-DFT approximations in mod-
eling molecular dissociation. The key ingredient of the proposed
procedure is the effective co-motion function used to construct an
approximate conditional density [Eq. (7)]. In this work, the effective
co-motion functions tested are extracted from the exact condi-
tional potential, along with the SCE and EXX conditional potentials.
We have obtained excellent accuracy using the exact potential as
input and important indications that our strategy is promising
even with approximate starting ingredients, correctly reducing the
perfect correlation in SCE and bringing in some correlation in
EXX (Fig. 5). Several aspects need to be further investigated and
refined, building upon the fundamental concepts proposed in this
work. For instance, the discontinuity obtained in the effective co-
motion functions appears to be an artifact of the one-dimensional
framework. Therefore, testing our strategy on simple spherical sys-
tems seems a natural next step. Additionally, to fully evaluate our
strategy, it will be important to test approximate functionals inspired
by our approach—such as the modified SCE potential—by assessing
their self-consistent performance. This would help verify whether
the imperfect but improved step structure resulting from our con-
struction can address the long-standing problem of molecular frag-
ments with fractional charges, which has plagued the application of
KS-DFT to molecular dissociation for decades. In general, models
of more widely applicable co-motion functions need to be devel-
oped for our strategy to become suitable for systems of chemical
interests. With widely applicable co-motion functions developed,
our strategy can be applied self-consistently: starting from the co-
motions and the density, an approximate conditional density can be
constructed, which in turn allows for the iterative construction of
a well-structured approximate Hartree-XC potential until conver-
gence is reached. The concept of intracule density or the Overhauser
model,⁴⁸ where the probability of having an electron pair separated
by a distance r₁₂ is the main variable, or the recently proposed
electron-centric approach⁴⁹ may be exploited to derive effective
co-motion functions in the future.

The supplementary material provides additional figures and
detailed discussions on specific aspects of this work, offering further
clarifications and context to complement the main text.

This project received funding from the European Union's
Horizon 2020 Research and Innovation program under Marie
Skłodowska-Curie Grant Agreement No. 101104947.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

**Sara Giarrusso:** Conceptualization (lead); Data curation (lead); Formal analysis (equal); Funding acquisition (lead); Writing – original draft (lead); Writing – review & editing (equal). **Federica Agostini:** Conceptualization (supporting); Data curation (supporting); Formal analysis (equal); Funding acquisition (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

DATA AVAILABILITY

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

REFERENCES

¹P. Hohenberg and W. Kohn, *Phys. Rev.* **136**(3B), 864 (1964).

²W. Kohn and L. J. Sham, *Phys. Rev.* **140**(4A), 1133 (1965).

³C. O. Almbladh and U. von Barth, in *Density Functional Methods in Physics*, edited by R. M. Dreizler and J. da Providência (Springer, 1985), pp. 209–231.

⁴M. A. Buijse, E. J. Baerends, and J. G. Snijders, *Phys. Rev. A* **40**, 4190 (1989).

⁵O. Gritsenko, R. van Leeuwen, and E. J. Baerends, *J. Chem. Phys.* **101**, 8955 (1994).

⁶O. V. Gritsenko and E. J. Baerends, *Phys. Rev. A* **54**, 1957 (1996).

⁷O. V. Gritsenko, R. van Leeuwen, and E. J. Baerends, *J. Chem. Phys.* **104**, 8535 (1996).

⁸O. V. Gritsenko and E. J. Baerends, *Theor. Chem. Acc.* **96**, 44 (1997).

⁹E. J. Baerends and O. V. Gritsenko, *J. Phys. Chem. A* **101**, 5383 (1997).

¹⁰S. van Gisbergen, P. Schipper, O. Gritsenko, E. Baerends, J. Snijders, B. Champagne, and B. Kirtman, *Phys. Rev. Lett.* **83**, 694 (1999).

¹¹D. G. Tempel, T. J. Martínez, and N. T. Maitra, *J. Chem. Theory Comput.* **5**, 770 (2009).

¹²N. Helbig, I. V. Tokatly, and A. Rubio, *J. Chem. Phys.* **131**, 224105 (2009).

¹³I. G. Ryabinkin, S. V. Kohut, and V. N. Staroverov, *Phys. Rev. Lett.* **115**, 083001 (2015).

¹⁴M. J. P. Hodgson, J. D. Ramsden, and R. W. Godby, *Phys. Rev. B* **93**, 155146 (2016).

¹⁵Z.-J. Ying, V. Brosco, G. M. Lopez, D. Varsano, P. Gori-Giorgi, and J. Lorenzana, *Phys. Rev. B* **94**, 075154 (2016).

¹⁶I. G. Ryabinkin, E. Ospadov, and V. N. Staroverov, *J. Chem. Phys.* **147**, 164117 (2017).

¹⁷M. J. P. Hodgson, E. Kraisler, A. Schild, and E. K. U. Gross, *J. Phys. Chem. Lett.* **8**, 5974 (2017).

¹⁸A. Kaiser and S. Kümmel, *Phys. Rev. A* **98**, 052505 (2018).

¹⁹S. Giarrusso, S. Vuckovic, and P. Gori-Giorgi, *J. Chem. Theory Comput.* **14**, 4151 (2018).

²⁰S. Giarrusso, P. Gori-Giorgi, and K. J. Giesbertz, *Eur. Phys. J. B* **91**, 186 (2018).

²¹S. Giarrusso and P. Gori-Giorgi, *J. Phys. Chem. A* **124**, 2473 (2020).

²²S. Giarrusso *et al.*, *J. Chem. Theory Comput.* **18**, 4762 (2022).

²³J. Kocák *et al.*, *J. Phys. Chem. Lett.* **12**, 3204 (2021).

²⁴J. Kocák, E. Kraisler, and A. Schild, *Phys. Rev. Res.* **5**, 013016 (2023).

²⁵D. Dar, L. Lacombe, and N. T. Maitra, *Chem. Phys. Rev.* **3**, 031307 (2022).

²⁶E. Baerends, *Phys. Rev. Lett.* **87**, 133004 (2001).

²⁷O. Gritsenko and E. J. Baerends, *Int. J. Quantum Chem.* **106**, 3167 (2006).

²⁸F. Furche, *Phys. Rev. B* **64**, 195120 (2001).

²⁹M. Fuchs and X. Gonze, *Phys. Rev. B* **65**, 235109 (2002).

³⁰M. Grüning, O. Gritsenko, and E. Baerends, *J. Chem. Phys.* **118**, 7183 (2003).

³¹J.-D. Chai, *J. Chem. Phys.* **136**, 154104 (2012).

³²M. Hellgren and T. Gould, *J. Chem. Theory Comput.* **15**, 4907 (2019).

³³J. Wang and E. J. Baerends, *Phys. Rev. Lett.* **128**, 013001 (2022).

³⁴G. Hunter, *Int. J. Quantum Chem.* **9**, 237 (1975).
<footnote n="35"><sup>35</sup>G. Hunter, *Int. J. Quantum Chem.* **9**, 311 (1975).</footnote>

<footnote n="36"><sup>36</sup>A. Schild and E. Gross, *Phys. Rev. Lett.* **118**, 163202 (2017).</footnote>

<footnote n="37"><sup>37</sup>M. Seidl, *Phys. Rev. A* **60**, 4387 (1999).</footnote>

<footnote n="38"><sup>38</sup>M. Seidl, P. Gori-Giorgi, and A. Savin, *Phys. Rev. A* **75**, 042511 (2007).</footnote>

<footnote n="39"><sup>39</sup>P. Gori-Giorgi, G. Vignale, and M. Seidl, *J. Chem. Theory Comput.* **5**, 743 (2009).</footnote>

<footnote n="40"><sup>40</sup>A. Mirtschink, M. Seidl, and P. Gori-Giorgi, *J. Chem. Theory Comput.* **8**, 3097 (2012).</footnote>

<footnote n="41"><sup>41</sup>R. Sharp and G. Horton, *Phys. Rev.* **90**, 317 (1953).</footnote>

<footnote n="42"><sup>42</sup>J. D. Talman and W. F. Shadwick, *Phys. Rev. A* **14**, 36 (1976).</footnote>

<footnote n="43"><sup>43</sup>R. Pederson, J. Chen, S. R. White, K. Burke *et al.*, *Phys. Rev. B* **105**, 245138 (2022).</footnote>

<footnote n="44"><sup>44</sup>M. Seidl, S. Giarrusso, S. Vuckovic, E. Fabiano, and P. Gori-Giorgi, *J. Chem. Phys.* **149**, 241101 (2018).</footnote>

<footnote n="45"><sup>45</sup>A. Gerolin, J. Grossi, and P. Gori-Giorgi, *J. Chem. Theory Comput.* **16**, 488 (2019).</footnote>

<footnote n="46"><sup>46</sup>S. Vuckovic and P. Gori-Giorgi, *J. Phys. Chem. Lett.* **8**, 2799 (2017).</footnote>

<footnote n="47"><sup>47</sup>S. Vuckovic and H. Bahmann, *J. Chem. Theory Comput.* **19**, 6172 (2023).</footnote>

<footnote n="48"><sup>48</sup>A. W. Overhauser, *Can. J. Phys.* **73**, 683 (1995).</footnote>

<footnote n="49"><sup>49</sup>P.-O. Roy, T. Henkes, and M. Ernzerhof, *J. Chem. Phys.* **159**, 151104 (2023).</footnote>