Response Potential in the Strong-Interaction Limit of Density Functional Theory: Analysis and Comparison with the Coupling-Constant Average.

Using the formalism of the conditional amplitude, we study the response part of the exchange-correlation potential in the strong-coupling limit of density functional theory, analyzing its peculiar features and comparing it with the response potential averaged over the coupling constant for small atoms and for the hydrogen molecule. We also use a simple one-dimensional model of a stretched heteronuclear molecule to derive exact properties of the response potential in the strong-coupling limit. The simplicity of the model allows us to unveil relevant features also of the exact Kohn-Sham potential and its different components, namely the appearance of a second peak in the correlation kinetic potential on the side of the most electronegative atom.

Introduction

Kohn–Sham (KS) density functional theory (DFT)[1] is the most used tool in quantum chemistry calculations thanks to its ability to predict properties of interest of a variety of physical, chemical, and biochemical systems at an acceptable computational cost with a reasonable accuracy. Nonetheless, there are still many relevant cases, typically when electron correlation plays a prominent role, in which current KS DFT methodologies are deficient, making the quest for new approximations to the unknown piece of information in DFT, the so-called exchange–correlation (XC) functional, an active research field.[2−4] This quest for a better (more versatile or accurate) and at the same time computationally affordable XC functional cannot proceed without a synchronized understanding of its exact properties and a constant search to find new ones that can act as constraints to build approximations.[5−8] In this context, a very important role is played by studies[9−27] focusing on the XC potential given by the functional derivative of the XC functional, , whose properties are crucial, for example, to accurately predict static electric polarizabilities and band gaps, and to correctly describe strongly correlated systems and bond breaking.

Pioneering work in this direction was pursued by Baerends and co-workers, who have analyzed the XC potential, deriving exact expressions in terms of wave functions and KS quantities.[9,11,13−15,28] Their work builds on the theory of conditional probability amplitudes first developed by Hunter,[29,30] which yields an exact differential equation for the square root of the density,[30] and was introduced in a DFT context by Levy, Perdew, and Sahni.[31] Baerends and co-workers have applied the same formalism to the KS Hamiltonian, deriving an insightful and exact decomposition of the XC potential, into so-called kinetic, response, and XC hole terms,[9,11,13−15,28] showing that each contribution has different properties and peculiarities that should be approximated with different standards.[28,32−34] For example, the response part builds a step structure in the KS potential of a stretched heteronuclear molecule, and the kinetic part builds a peak in the midbond region of a stretched bond.[15] Lately this subject has gained renovated interest for various reasons, spanning from the construction of KS potentials from wave functions in finite basis sets,[20−22,26] to the use of response potential approximations to compute band gaps[33] and correcting semilocal functionals,[34] to further investigations of the step structure for molecular dissociation[16,27] and of the kinetic peak for Mott insulators.[17,24]

At the same time, in recent years, the mathematical structure of the limit of infinite interaction strength in DFT, corresponding to the so-called strictly correlated electrons (SCE) functional, has been thoroughly investigated.[35−40] The SCE functional has a highly nonlocal density dependence, but its functional derivative can be computed via a physically transparent and rigorous auxiliary equation, which provides a powerful shortcut to access the corresponding XC potential.[41] This SCE XC potential has been used in the KS framework to compute properties of electrons confined at low density, close to the “Wigner molecular” regime.[41−44] Despite how extreme the SCE limit might sound, it has the advantage of unveiling explicitly how the density is transformed into an electron–electron interaction, in a well-defined asymptotic case (low density or strong interaction) for the exact XC functional.[39,40] Its peculiar mathematical structure has already inspired new approximations, in which instead of the traditional DFT ingredients (local density, density gradients, KS orbitals) certain integrals of the density play a crucial role.[45−47]

The SCE limit, however, has never been analyzed from the point of view of the conditional amplitude framework, and nothing is known about the behavior of the different components of the corresponding XC potential. It is the main purpose of this work to fill this gap. We start by generalizing the Schrödinger equation for the square root of the density to any coupling-strength λ value, analyzing its features in the λ → ∞ (or SCE) limit (section ). We derive and analyze the response part of the SCE exchange–correlation potential (section ), and we compare it with the one from the coupling-constant average formalism for small atoms and the H2 molecule (section ). Using a one-dimensional model[16,17,27] for the dissociation of a heteroatomic molecule, we analyze in this limit the SCE and exact exchange–correlation potentials, focusing on the step structure and further analyzing the kinetic potential for the physical coupling strength (section ).

Strong-Interaction Limit of the Effective Equation for the Square Root of the Density

Consider the λ-dependent Hohenberg–Kohn functional within the constrained-search definition[48]assuming that ρ is v-representable for all λ, one can write a series of λ-dependent Hamiltonian with fixed densitywhere V̂λ = ∑ vλ(r) andis the local external potential that delivers the prescribed density as the ground-state density of Hamiltonian eq at each λ, i.e., ρλ(r) = ρ1(r) ≡ ρ(r), and Ψλ(1, ..., N)—with 1, ..., N the spin-spatial coordinates of the N electrons—is the ground state wave function of Hamiltonian eq at each λ. Following refs (9 and 29−31), we partition the Hamiltonian eq aswhere Ĥλ* is the Hamiltonian eq deprived of one particle (which is in general different from the λ-dependent Hamiltonian of the physical (N – 1)-electron system), and factorize the wave function Ψλ(1, ..., N) aswith 1 = rσ being the spatial-spin coordinates of electron 1 taken as a reference. The function Φλ(σ, 2, ..., N|r) is called conditional amplitude and describes the behavior of the remaining N – 1 electrons as a parametric function of the position r of electron 1. Notice that the conditional amplitude, when integrated over its N – 1 variables, is normalized to 1 for all values of r.

By applying eq to eq , multiplying by Φλ*(σ, 2, ..., N|r) to the left, and integrating over all variables except r as in refs (9 and 29−31), one obtains an effective equation for the square root of the density for any λ-valuewhereand E0,λ* is the ground-state energy of the N – 1 system in the same effective potential as the N-particle one, i.e., of Ĥλ* of eq (thus E0,λ* = E0,λ only for λ = 1). The various components of the effective potential have each its own physical meaning and peculiar features, and have been carefully studied by many authors, at λ = 1 and λ = 0.[9,11,13−17,22,24,27,28] The term vλ,(r) is related to the response potential (see section below) and is given bywhere the subtraction of the quantity E0,λ* makes this potential go to zero when |r| → ∞, as in this case the conditional amplitude usually collapses to the ground state of the system deprived of one electron, if accessible (see refs (49 and 50) for an in-depth discussion and exceptions). The kinetic potential isand it also goes usually to zero when |r| → ∞, as the conditional amplitude in this case becomes insensitive to the position of the reference electron (again, see refs (49 and 50) for an in-depth discussion and exceptions). Finally, the conditional potential iswhere it should be noted that vλ,cond = vH + vxc,holeλ and it tends manifestly to zero when |r| → ∞. For any finite λ, the difference E0,λ – E0,λ* in eq equals minus the exact ionization potential Ip of the physical system, which dictates the asymptotic decay of the density[31,51]

Similarly, the sum vλ, eff(r) + vλ(r) is obviously λ-independent, as the density is the same for all coupling strengths λ. It is exactly by equating vλ, eff(r) + vλ(r) at λ = 0 and λ = 1 that Baerends and co-workers could derive their insightful decomposition of the KS potential,[9,11,13−15,28] as this gives an equation for vλ=0 (i.e., the KS potential) in terms of wave function and KS orbital quantities.

General Structure of the λ → ∞ Limit

When λ → ∞, the Hamiltonian of eq has the expansion[35,36,39,40]where V̂SCE = ∑ vSCE(r) is the one-body potential that minimizes the classical potential energy operator V̂ee + V̂SCE and delivers the prescribed ground-state density ρ(r) according to eq .[35,36,42] The modulus squared of the corresponding wave function usually collapses into a distribution that can be written as[35,36,52]where the co-motion functions f(r) describe the perfect correlation between the N electrons. They are nonlocal functionals of the density satisfying the equationwhich ensures that the probability of finding one electron at position r in the volume element dr be the same as finding electron i at position f(r) in the volume element df(r). They also satisfy cyclic group properties (for a recent review on the mathematical properties of the co-motion functions, see ref (53)):

The corresponding SCE functional, given by[35,54]yields the strong-coupling (or low-density) asymptotic value of the exact Hartree-exchange–correlation functional.[39,40] Despite the extreme nonlocality of VeeSCE[ρ], its functional derivative can be computed from the exact force equation[35,41]

According to eq , the one-body potential vSCE(r) of eq is exactly equal to minus vHxcSCE(r): in fact, the gradient of vHxcSCE(r) represents the net repulsion felt by an electron in r due to the other N – 1 electrons at positions f(r), while vSCE(r) appearing in the λ → ∞ Hamiltonian of eq exactly compensates this net force, in such a way that the classical potential energy operator V̂ee + V̂SCE is stationary (and minimum) on the manifold parametrized by the co-motion functions. Equation defines vHxcSCE(r) up to a constant, which is fixed by imposing that both vHxcSCE(r) and vSCE(r) = −vHxcSCE(r) go to zero when |r| → ∞.

The effective eq for in the SCE limit can be easily understood if we divide both sides by :

When λ → ∞, we see that the first term in the left-hand-side of eq goes to zero, as the density ρ(r) does not change with λ and it is well-behaved, with the exception of the values of r on top of the nuclear positions R, where the density has a cusp and yields back the Coulombic divergence. Nagy and Jánosfalvi[55] have carefully analyzed the λ → ∞ behavior at the nuclear cusps in , showing that for all λ values the kinetic divergence at a nucleus of charge Z at position R cancels exactly the external potential . We can then safely disregard both the kinetic divergence and the Coulombic divergence in the λ → ∞ limit. The other case, which we do not consider here, where this term may diverge is when the KS highest-occupied molecular orbital (HOMO) has a nodal plane that extends to infinity.[49,50,56]

All the remaining terms, except for vλ,kin(r), will tend to a finite, in general nonzero, limiting value, as they grow linearly with λ (for example vλ(r) → −λvHxcSCE(r) of eq . Notice that vλ,cond(r) has been already defined with the factor λ in front; see eqs and 10. The only delicate term is vλ,kin(r) of eq , which contains the gradient of a conditional amplitude that is collapsing into a distribution. Several results in the literature suggest[36,39,57] that this term grows with λ only as , thus still vanishing with respect to the other terms. However, we should keep in mind that no rigorous proof of this statement is available at present. Nonetheless, as shown below, the SCE limit provides a perfectly consistent treatment of the leading order of eq when λ → ∞, providing further evidence that the kinetic potential vλ,kin(r) should be subleading in eq .

Conditional Probability Amplitude and Ionization Potential at the SCE Limit

We can now use eq to find the conditional amplitude in the SCE limit and to partition the corresponding effective potential into its two components of eqs and 10 (as said, the kinetic part disappears in this limit). Notice that although in eq we have considered only one possible permutation of the N electrons (compare the expression, e.g., with eq 14 in ref (54)), this does not affect the derivations below, as explicitly shown in Appendix A. Integrating over s, we getand applying eq we findEquation shows that the conditional amplitude gets a very transparent meaning in the SCE limit, as it simply gives the position of the other N – 1 electrons as a function of the position r of the first electron.

In what follows we label with “SCE” the terms that survive when we take the limit λ → ∞ of eq . We then use eq to evaluate in this limit vSCE(r):

Now we use the fact that the ground-state energy of the N-particle system with density ρ(r) at the SCE limit is simply given by the value of the classical potential energy V̂ee + V̂SCE on the manifold parametrized by the co-motion functionswhich allows us to rewrite the first term of eq as

The last two terms in the right-hand side of eq vanish for |r| → ∞. On the other hand, by construction v(r) → 0 when |r| → ∞, and thus necessarilyand we obtain the final simple expression for vSCE(r):

Equation might look puzzling, but one could also expect it from the fact that, as said, in the SCE limit we obtain the quantities that survive in eq when we take the λ → ∞ limit. This means that the difference E0,λ – E0,λ* grows linearly with λ for large λ:

Then we see that the only way in which eq can be satisfied when λ goes to infinity is if eq holds. Indeed this result was already implicit in ref (35), where it was noticed that the configuration with one electron at infinity must belong to the degenerate minimum of the classical potential energy operator V̂ee + V̂SCE. Equation shows that also for the next leading order there should be no energy cost to remove one electron, a statement that is implicitly contained in ref (36).

Notice that the zero ionization energy of eq concerns the λ → ∞ Hamiltonian in the adiabatic connection of eq . A very different result is obtained if vHxcSCE(r) is used as an approximation for the Hartree-XC potential in the self-consistent KS equations, where the corresponding KS HOMO eigenvalue has been found to be very close to minus the exact ionization potential for low-density systems,[41,44] displaying the correct step structure when the number of electrons is changed in a continuous way.[43]

Different Types of Response Potentials: vresp(r), v̅resp(r), and vrespSCE(r)

In order to compare the SCE response potential with the physical one, we first review the different possible definitions that appear in the literature[9,11,13−15,28,58] for this term, and fully define the response potential in the SCE limit.

We start from the pair density P2λ(r,r′), associated with the Hamiltonian in eq according to the formulaand the corresponding exchange–correlation pair-correlation function gxcλ(r,r′) at a given coupling strength λ

We also define the coupling-constant averaged (CCA) pair-correlation function g̅xc(r,r′)

In what follows we use the subscript s when the quantity of interest refers to the KS or λ = 0 case and we omit the subscript λ when it refers to the physical system λ = 1.

Response Potential in Terms of Kinetic and Interaction Components

The XC functional of KS DFT can be written aswhere

If we now take the functional derivative of the XC energy with respect to the density, we can recognize four different contributions to the XC potential:[11,15]whereand

We can also group the potentials in eqs and 35 into one total response potential, vresp(r)

By inserting the KS Slater determinant and the λ = 1 wave function into eq , it has been shown[11,15] that

Response Potential from the Coupling-Constant Averaged XC Hole and Comparison between vresp(r) and v̅resp(r)

The XC energy can be also written in terms of the CCA g̅xc(r,r′)as the integration over λ allows recovering the kinetic contribution to Exc[ρ].[59−61] Taking the functional derivative of eq , we obtain two terms[34]whereand

Equation defines the quantity v̅resp(r), but looking at eq one can also determine it aswhich is how we have computed the response potential in section . Comparing eqs and 39, we have

Intuitively, one would expect that the sum of the response parts of the left-hand side of eq equals the response part in the right-hand side, and that so do the remainders on both sides. However, this is not true, and in general we have (see, e.g., eq in section and the CCA response potentials shown in Figures and 7 in comparison with their response potentials at physical λ shown in the works cited in the figures’ captions)

Figure 2

Comparison between v̅resp(r) and vrespSCE(r) for the He atom. In the top-right inset the CCA response potential of He is zoomed in to allow a closer comparison with its response potential at full coupling strength, vresp(r), shown in Figure 3c of ref (15).

Figure 7

Comparison between v̅resp(r) and vrespSCE(r) for the H2 molecule at the equilibrium distance, plotted along the internuclear axis, with the origin of the axes being at the bond midpoint. In the top-right inset the CCA response potential of H2 is zoomed in to allow a closer comparison with its response potential, vresp(r), shown in Figure 3a of ref (15).

Response Potential for the SCE Limit by Means of Energy Densities

The two response potentials defined in eq and eq can be both thought of as a measure that answers the question,[11,15,28] “How sensitive is the pair-correlation function on average to local changes in the density?” Therefore, it seems interesting to ask what happens to it when electrons are perfectly correlated to each other, i.e., in the SCE limit.

From the AC formalism of eq , the integrated form of eq iswhere is the (global) AC integrand, defined as

We can generalize eq to any XC energy along the adiabatic connection as

Using the expansion of the (global) AC integrand in the strongly interacting limit[35,36,39,40,52]to first order we obtain

Defining the SCE XC energy asand inserting eq into eq , we get the simple relation

In order to derive the response potential in the SCE limit, we can now start by taking the functional derivative of with respect to the density, similarly to what is done in eqs and 39. However, in recent years, focus has been brought to the importance of using the local counterpart of the global integrand , i.e. the so-called energy density, wλ[ρ](r). This different approach is especially important for DFAs (density functional approximations) in view of the fact that local models are generally more amenable to the construction of size-consistent and accurate methods than their global counterparts.[46,62,63] The local analogue of eq for the XC energy becomes

Whenever energy densities are used, it is crucial to define a “gauge” within which all the quantities taken into account are computed consistently at different λ-values, being the choice of wλ(r) not unique. A physically sound and commonly used gauge of the energy density is the one given in terms of the electrostatic potential of the XC hole, which corresponds toandwhere w1[ρ](r) in the literature is also labeled as w[ρ](r) or wxc[ρ](r), while for w̅[ρ](r) the symbol w̅xc[ρ](r) or ϵxc[ρ](r) is also commonly used.

The corresponding energy density at λ = 0, w0[ρ](r), is usually also labeled ϵ(r) or w[ρ](r). For λ → ∞ we have, in this gauge[54]where vH(r) is the Hartree potential. Let us write the AC integrand at λ → ∞, in analogy to eq , as

Taking the functional derivative of with respect to the density, we obtainwhereand

Finally, inserting the explicit expression for the energy density for λ → ∞ (eq ), we find the response potential at the SCE limitwhich is exactly equal to vSCE(r) of eq . Notice that the SCE response potential of eq scales linearly with respect to uniform scaling of the density:[5]where ργ(r) ≡ γ3ρ(γ r) is a scaled density.

SCE Response Potential for a Two-Electron Density

When the number of electrons equals two, we also have another expression for computing vrespSCE(r). In this case the SCE total energy E0,SCE of section is equal towhere the right-hand side is the value of the SCE potential energy on the manifold parametrized by the co-motion function. This value is a degenerate minimum, meaning that we can evaluate it at any point lying on the manifold, such as for |r| → ∞ (for a nice illustration of the degenerate minimum of the SCE potential energy, the interested reader is addressed to Figure 1 of ref (57)). When |r| → ∞, the potential vHxcSCE(r) is gauged to go to zero. At the same time, the co-motion function f(r) will tend to a well-defined position r0 well inside the density, i.e., f(r→∞) → r0. We thus have

Combining eqs and 64 we find

Examples of CCA and SCE Response Potentials

We have computed the SCE response potential, vrespSCE(r), for small atoms and for the hydrogen molecule at equilibrium distance; for this latter case and for the species H–, He, Be, and Ne also, accurate CCA response potentials v̅resp(r) have been obtained. Notice that, in previous works, several authors[9,11,13−15,20−22,26,28,64] have computed the response potential at physical coupling strength, vresp(r), of eqs and 37. To our knowledge, accurate CCA response potentials v̅resp(r) (eq ) are reported here for the first time. In Appendix B we also briefly discuss the extent of the error resulting from combining data coming from different methods, namely from the Lieb maximization procedure[63,65,66] and Hylleraas-type wave functions[67,68] or quantum Monte Carlo calculations[10,12,67] as explained in the next sections. In the figures all quantities are reported in atomic units.

Computational Details for Atomic Densities

For the sake of clarity, we treat in separate sections the computation of vrespSCE(r) (section ) and v̅resp(r) for atoms (section ).

SCE Response Potential

The calculation of vrespSCE(r) for spherical atoms is based on the ansatz for the radial part of the co-motion functions reported in ref (35). These co-motion functions are exact for N = 2,[37] and for N > 2 they give either the exact SCE solution or get very close to it.[53] Moreover, even when they are not truly optimal, the corresponding potential still satisfies eq .[53] This means that we are in any case using a perfectly correlated wave function to compute a meaningful response potential. The radial co-motion functions f(r) of ref (35) are givenwhere N is the number of electrons; Ne(r) is the cumulant function:Ne–1(y) is its inverse, defined for y ∈ [0,N), and a are the (radial) distances for which Ne(a) = i, with i integer. These radial co-motion functions give the distances from the nucleus of the remaining N – 1 electrons as a function of the distance r of the first one. The relative angles between the electrons are found by minimizing the total repulsion energy for each given r.[35,44] The SCE potential, vHxcSCE(r), is then obtained by integration of eq . Finally, we apply eq (or, equivalently for N = 2, eq ) to get the SCE response potential.

This procedure is very “robust”, meaning that we have obtained comparable SCE response potentials using densities of different levels of accuracy. The densities we have used were obtained from the following:

(A) CCSD calculations and aug-cc-pCVTZ basis set stored on a 0.01 bohr grid; see ref (63)

(B) Hylleraas-type wave functions; see refs (67 and 68) for the two-electron systems and quantum Monte Carlo calculations; see refs (10, 12, and 69) for the others

The cumulant function of eq was computed either with simple interpolations between the grid points of a given density or in some cases (for H–, He, and Li+) with explicitly fitted densities, constrained to satisfy the cusp condition and the correct asymptotic behavior.

Group A regards all the systems taken into account. Group B regards the species H–, He, Be, and Ne. The figures in section only show the SCE response potential coming each time from the most accurate available density.

Coupling-Constant Averaged Response Potential

The equation used in practice to compute v̅resp(r) iswhere w̅(r) is given in eq , and was calculated by averaging the energy densities wλ(r) at each r over the interval [0,1] with an increment Δλ = 10–1. The wλ(r) were obtained through the Lieb maximization procedure and taken from refs (62, 63, and 70). The XC potentials were taken instead from Hylleraas-type calculations[67] or quantum Monte Carlo results,[10,12,67] as they were overall more accurate. This choice is further validated in Appendix B.

Computational Details for the Hydrogen Molecule

For the hydrogen molecule a different approach–i.e. the “dual Kantorovich formulation” in the framework of optimal transport theory[37,71] – was used for the computation of the SCE potential and thus of the SCE response potential. The basic idea relies on finding the SCE potential as a result of a nested optimization on a parametrized expression which has the correct asymptotic behavior, the correct cylindrical symmetry and models the barrier region in the midbond. From the optimized potential one derives the co-motion function by inverting eq ; for details see ref (71).

For the CCA energy density, w̅(r), exactly the same procedure described for atoms has been used.

The XC potential for the physical system in this case was obtained within the Lieb Maximisation procedure itself as in ref (63), namely as the optimized effective potential that keeps the density fixed minus the Hartree potential and the potential due to the field of the nuclei (see also Appendix B for data validation).

Results and Discussion

We start by showing in Figure the CCA and SCE response potentials for the H– anion: we see that on average the SCE response potential is larger than the CCA one, but there is an intermediate region, in the range 1.7 ≲ r ≲ 5.2, where the CCA values are above the SCE ones. Since the SCE response potential does not contain any information on how the kinetic potential is affected by a change in the density, this could be a region where the contributions coming from the kinetic correlation response effects overcome the Coulomb correlation ones, even though we cannot exclude that already the mere Coulombic contribution to correlation is higher in the physical case. Indeed, it has been shown that the SCE pair density can be insensitive to changes in certain regions of the density.[72]

Figure 1

Comparison between v̅resp(r) and vrespSCE(r) for the H– anion.

In Figure we report a similar comparison for the He atom density. Since He is less correlated than H–, in this case the CCA potential v̅resp(r) differs even more from the SCE one. Comparing the two species H– and He between each other, one can further observe that the value of the distance at which the response potential of the species i has a maximum, r, is also shifted leftward (closer to the nucleus) when going from Z = 1 to Z = 2, reflecting the contraction of the density. This information is also mirrored in the SCE limit by the shift in the a1 values appearing in 66b for the computation of the co-motion functions for the two species. Indeed we find that .

As could be expected from eq , the response potential at the SCE limit shows an almost perfect scaling behavior along the He series when we increase the nuclear charge Z. This is shown in Figure , where we report the scaled potentials, , as a function of the scaled coordinate Zr. More diffuse densities, like He and H–, deviate from the linear-scaling trend, showing increasing correlation effects in their densities. Such correlation effects (curve lying below the uniformly scaled trend for small r and above for large r) are stronger closer to the nucleus. In the top-right inset of Figure , we show only the values of the maxima of the SCE response potential of each species divided by its nuclear charge, as a function of Z. In this inset also a hypothetical system with nuclear charge Zcrit = 0.911 028 9, the minimum nuclear charge that can still bind two electrons (see ref (67)), is included.

Figure 3

Scaled SCE response potentials, , as a function of the scaled coordinate Zr for the He series from H– up to Ne8+. In the inset, in which only the “slice” at r = 0 (i.e., the maximum values of the SCE response potentials) is plotted as a function of the nuclear charge Z; also the hypothetical system with Z = Zcrit (see text) is considered.

In the upper panel of Figure we show the SCE and the CCA response potentials for the Be atom together with the exchange contribution vresp,(r) (corresponding to λ = 0), and the correlation contributions obtained by subtracting vresp,(r) from v̅resp(r) and vrespSCE(r). As it was found in ref (28), the exchange-only response potential shows a clear step structure in the region of the shell boundary. The total CCA response potential also shows a step at the same position, while the SCE response potential has a kink. The kink can be understood by looking at the shape of the radial co-motion functions (see eq and the lower panel of Figure ), which determine the structure of the SCE response potential according to eq . The SCE reference system correlates two adjacent electron positions in such a way that the density between them exactly integrates to 1; therefore, the a appearing in eqs and 66b are simply the shells that contain always one electron each.[73] For the case of Be, the kink appears at the corresponding a2 value, which is very close to the shell boundary. In fact, when the reference electron is at distance r ≈ a2 from the nucleus, the second electron is found at this same distance (but on the opposite side with respect to the nucleus), while the third electron is very close to the nucleus and the fourth is almost at infinity. This situation results in an abrupt change of the pair density for small variations of the density, as particularly the position of the fourth electron changes very rapidly with small density variations. Another interesting feature we can observe from Figure is that the Coulomb correlation contribution to the CCA response potential, v̅resp,(r), appears to be negative inside the entire 1s shell region. Furthermore, while the total physical response potential is always below the SCE one, the exchange part appears to be higher in a region quite close to the shell boundary (0.6 ≲ r ≲ 1.0). This results in the Coulomb correlation contribution for the SCE-limit case, vresp,SCE(r), to be also negative in that region.

Figure 4

Total response potentials v̅resp(r) and vrespSCE(r), and their components vresp,(r), v̅resp,(r), and vresp,SCE(r) (upper panel) and radial co-motion functions (lower panel) for the Be atom.

In the upper panel of Figure we show the SCE response potential and its correlation part for the Ne atom. The SCE response potentials vrespSCE(r) and vresp,SCE(r) are numerically less accurate, due to the higher dimensional angular minimization. Nevertheless, the relation between its structure and the corresponding co-motion functions in the lower panel of Figure is clearly visible. We also show the CCA response potentials together with the separate exchange and correlation contributions. Differently from the Be atom, neither the total response potential nor any single correlation contribution (CCA or SCE) is anywhere negative. Still the structure is very similar, showing two steps in the vresp,(r), one very tiny at around 0.1 and another at around 0.4 distance from the nucleus, and two wells in the v̅resp,(r). In Figure we show only the CCA correlation contributions to the CCA response potential of the two species for closer comparison.

Figure 5

Total response potentials v̅resp(r) and vrespSCE(r), and their components vresp,(r), v̅resp,(r), and vresp,SCE(r) (upper panel) and radial co-motion functions (lower panel) for the Ne atom.

Figure 6

Correlation parts of the CCA response potential, v̅resp,, for the Be and Ne atoms.

In Figure the CCA response potential for the hydrogen molecule at equilibrium distance is shown, together with the SCE one. It is interesting to compare Figure with Figure 3a of ref (15), where the response potential vresp(r) of eq was reported, together with other components of the XC potential. The response potential at full coupling strength for the same system is also shown in Figure 4 of ref (64), albeit a minus sign and a constant shift. The overall structure is completely different: in the case shown here there is a local minimum of v̅resp(r) at approximately 1 bohr distance from the bond midpoint, while vresp(r) shown in refs (15 and 64) has a maximum located at the nuclei. This must necessarily be due to the coupling-constant average procedure, in which the responses of the kinetic and Coulombic contributions are taken into account in two different ways. It is then important to keep these different features in mind when one wants to model the response potential, depending on whether the target is v̅resp(r) or vresp(r).

Simple Model for a Stretched Heteronuclear Dimer

The purpose of this section is to analyze the response potential in the SCE limit for the very relevant case of a dissociating heteroatomic molecule, where the exact response potential is known to develop a characteristic step structure.[13,15−17,23,25,74] Although numerically stable KS potentials have been presented and discussed in the literature for small molecules,[20,26,27] an accurate calculation of the SCE potential for a stretched heterodimer is still not available. In fact, while with the dual Kantorovich procedure[71,75] it is possible to obtain accurate values of VeeSCE[ρ] for small molecules, the quality of the corresponding SCE potentials, particularly in regions of space where the density is very small, is not good enough to allow for any reliable analysis.

We then used a simplified one-dimensional (1D) model system, where only the two valence electrons involved in the stretched bond are treated explicitly. Several authors have used this kind of 1D model; the models have been proven to reproduce and to allow understanding of the most relevant features appearing in the exact KS potential of real molecules.[16,17,23,25] Here we approximate the density of the very stretched molecule as just the sum of the two “atomic” densitieswhere a and b mimic the different ionization potentials of the “atoms” (pseudopotentials or frozen cores) and the density is normalized to 2. We have chosen a > b; therefore, the more electronegative atom will be found to the right side of the origin (at a distance from it) and the less electronegative atom will be to the left.

In the last part of this section (section ), we inspect and reveal further features of the response potential also at physical coupling strength and put them closely in relation with the SCE scenario discussed in the first part (section ); this investigation is indeed made possible thanks to the simplicity of the model.

SCE Response Potential for the Model Stretched Heterodimer

In 1D, we have (see eq for comparison)and, as we have two electrons, there is only one of the “SCE shell” borders, a, appearing in eqs and 66b:

We have used the subscript “R ” because the distance a1 is a function of the separation between the centers of the exponentials in eq . Also, there is only one co-motion function that describes the position of one electron given the position x of the other, equal to[52,73,76]

We have stressed in section that the border of a shell that contains one electron coincides with the reference position at which one of the co-motion functions diverges. The same is true when x → a, except that in the one-dimensional case the electron that goes to infinity has to “reappear” on the other side, lim f(x) = ∓∞. Moreover, as we have only two electrons, we can use eq to compute vrespSCE(r)which further shows that

In Figure we show the SCE response potential compared to the “exact” v̅resp(x) for the model density of eq at internuclear separation R = 8, using a = 2 and b = 1. In Figure , we also show the local-density-approximation (LDA) CCA response potential v̅respLDA(x) computed, as in ref (34), via eq .

Figure 8

SCE response potential compared to the “exact” and the LDA v̅resp(x) for the model density in eq with a = 2, b = 1, and R = 8. The red dashed line highlights the position where x = a.

We stress that eq is the correct definition of v̅respLDA(x), since the energy density in LDA does not have any gauge ambiguity, being given exactly in terms of the electrostatic potential associated with the CCA exchange–correlation hole of the uniform electron gas.[77] For the one-dimensional ϵxcLDA, we have used the parametrization of Casula et al.,[78] in which the electron–electron Coulomb interaction is renormalized at the origin,[77] with thickness parameter b = 0.1. Notice that the SCE response potentials evaluated with the full Coulomb interaction 1/|x| or with the interaction renormalized at the origin[77] are indistinguishable on the scale of Figure , since in the SCE limit the electron–electron distance |x – f(x)| for a stretched two-electron “molecule” never explores the short-range part of the interaction.

The “exact” v̅resp(x) has been computed by inverting the KS equation for the doubly occupied ground-state orbital , disregarding the external potential given by attractive delta functions located at the “nuclei”, and assuming that, for the stretched molecule, the interaction between fragments is negligible (which is asymptotically true), while the contributions coming from the Hartree potential on each fragment (the self-interaction error) are exactly canceled by the XC hole. In other words, when R is large, we have vHxc(x) ≈ v̅resp(x) ≈ v(x) + vresp(x).

We see that, as is well-known, the LDA response potential completely misses the peak with the step structure of the “exact” v̅resp(x) being, instead, way too repulsive on the atoms,[34] and following essentially the density shape. The SCE response potential, instead, even though clearly not in agreement with the “exact” one, shows an interesting structure located at the peak of v̅resp(x), and also a sort of steplike feature.

In Figure we illustrate the behavior of the SCE response potential alone as the internuclear separation R grows, for the same values of a and b of Figure . We see that the SCE response potential, contrary to the exact one, does not saturate to a step height equal to the difference of the ionization potentials of the two fragments, ΔIp = |I – I|. On the contrary, vrespSCE(x) goes (although very slowly) to zero in the dissociation limit, similarly to what happens for the midbond peak in a homodimer, as explained in refs (24 and 73). This has to be expected, in view of the fact that, in the SCE limit, we are only taking into account the expectation of the Coulomb electron–electron interaction, which, when considering two distant one-electron fragments as in this case, is a vanishing contribution.[24] The fact that we still observe the SCE response structure for quite large R values is related to the nonlocality of the SCE potential and to the long-range nature of the Coulomb interaction. A kinetic contribution to SCE is clearly needed, something that is being currently investigated by looking at the next leading terms in the λ → ∞ expansion.[36,57]

Figure 9

SCE response potential for the model density in eq with a = 2, b = 1, and increasing internuclear distances, R.

The peak structure of the SCE response potential is located at a of eq , which is given by

If we compare this result with the one for the location of the step in the exact KS potential, given by eqs 27 and 29 of ref (16), we see that the two expressions differ by the termwhich becomes comparatively less important as the bond is stretched. In Figure we have reported the case a = 2, b = 1, and R = 8, for which eq gives , and the correction term for the actual position of the step,[16] which is also the position at which the kinetic peak has its maximum, xstep = xpeak, givesThe reason why, in spite of this significative correction, in Figure the peak of the “exact” v̅resp(x) visibly coincides with a will be clear in section .

Behavior of the Co-Motion Function for Increasing Internuclear Distance

The features of the SCE response potential can be understood by looking at how the co-motion function changes with increasing internuclear separation R. In the 1D two-electron case considered here, eq becomes

For R ≫ 0, when the reference electron (e1) is in the center of one the two “atomic” densities, e.g., at , the other electron (e2) is in the center of the other “atom”, . This is a simple consequence of the fact that the overall density is normalized to 2 and, if the overlap in the midbond region is negligible, for symmetry reasons, the area from to is exactly equivalent to the sum of the areas outside that range.

We see that after a critical internuclear distance, Rc, at which the overlap between the densities of the separated fragments becomes negligible, the slope of the co-motion function when e1 is in becomes equal toso that there is a region where , and, similarly, another region where , by interchanging e1 with e2. Notice that the extension of these regions is different for the two branches of eq and it is wider when the reference electron is around the less electronegative “atom” as can be seen in Figure , where we show the (numerically) exact

Figure 10

Derivative of the co-motion function for the model density in eq with a = 2, b = 1, and increasing internuclear distances, R.

There, the two regions clearly appear as left and right plateaus, with their extent increasing linearly with R. These plateaus are the signature of molecular dissociation: they are absent at equilibrium distance, and they start to appear as the overlap between the two densities is small. We see from eq that they encode information on the ratio between the ionization potentials of the two fragments.

Careful Inspection of the Exact Features of the KS Potential for the Dissociating AB Molecule

The model density ρ(x) of eq corresponds to an asymptotic simplification of different models that appeared in the literature to study the KS potential in the dimer dissociation limit.[16,17,23,25] Here we review in detail the properties of the KS potential and the two single contributions that can be extracted from this model, v(x) (eq ) and vresp(x) (eq or 37), also showing that a second peak in the kinetic potential appears on the side of the more electronegative “atom”, a feature that seemed to have been overlooked in previous studies. In order to study the dissociation regime, we use the Heitler–London wave function:where SAB = ∫ϕ(x) ϕ dx, and To compute the kinetic potential, in the dissociation limit, we can use eq and the conditional amplitude coming from the Heitler–London wave function considering SAB = 0, which yields the well-known expression[16,17]where we have used the fact that vkin(x) = v(x) as the kinetic KS potential is zero for a closed-shell two-electron system. Analogously, vresp(x) can be obtained from v of eq :whereand . Comparing these two contributions with the KS potential obtained from the density by inversion (subtracting the external potential due to the attractive delta peaks at the “nuclear” positions), we have in this limit, as already discussedsince vcond(x) goes to zero when the fragments are very far from each other. In Figure we show the potential obtained from the inversion of the KS equation with its two components v(x) and vresp(x).

Figure 11

Hartree-XC potential, vHxc(x), and its contributions v(x) and vresp(x), for a = 2, b = 1, and R = 8. The red dashed line highlights the position where x = a.

For this simple model, we have exact expressions regarding each component of the potential and their maxima, inflection points, and so forth. Some of these relevant analytic expressions are listed in Table .

Table 1

Some of the Relevant Analytic Features of the Analytic 1D Model Dimera

a ⩾ b
Part I
		84
	vc,kin(xeq) = vresp(xeq)
Part II
		85
		vresp(xstep(2)) = vresp(xstep(1))

The table has two parts: , and . In part I, xpeak(1) is the position at which the kinetic potential, v(x), has a maximum in between the two nuclear centers; xstep(1) is the (coinciding) position at which the response potential, vresp(x), has an inflection point. With the subscript “flex” we indicate the inflection point of both the total Hartree-XC potential and the kinetic potential; they are distinguished via an additional subscript, respectively “Hxc” and “k”. Finally, xeq is used to label the x-value at which v(x) and vresp(x) cross. In part II, the analogous quantities appearing, in this case, somewhere far from the midbond on the side of the more electronegative fragment are listed. For example, xpeak(2) is the second maximum of the kinetic potential, eq (top-right entry of part II), which also coincides with the second inflection point of the response potential as argued in the main text.

By looking at Table , one sees, for example, that the peak of the total Hartree-XC potential is not located where the peak of the kinetic correlation builds up. In particular the maximum of the Hartree-XC potential is found atwhich is exactly a (see eq and compare also Figure ). Thus, the Hartree-XC potential reaches its maximum when the density integrates to one electron (or the correct integer number of electrons in a general two-fragment case) because this is where the two fragments must be detached from one another. From a different perspective, this is a manifestation that the response and the kinetic correlation contributions in the dissociation limit are not independent and that their sum can be sometimes more meaningful than the separate contributions. Also, by playing around with the expressions in Table , one realizes that there can be misleading coincidental features. For example, the last entry of part I of Table , which is the analytic expression for the distance at which the kinetic correlation potential and the response potential equate, xeq, is such that the two contributions v(x) and vresp(x) cross exactly at a if a = 2 and b = 1 as in Figure , but this is not a general feature. Similarly, if we choose , then the height of the kinetic peak becomes equal to the height of the step and so on.

Note here that the features listed in Table are obtained for the zero-overlap case, SAB = 0, in eq . Nonetheless, they should become asymptotically exact in the dissociation limit.

Another feature that came to our attention and that—to the best of our knowledge—has not been discussed before, is the fact that the kinetic correlation potential has a second peak on the side of the more electronegative atom. This second maximum is located where the second inflection point of the response potential is; see Figure and eq in Table . To understand the appearance of the second peak, we can identify two regimes, A and B, by the leading exponential coefficient: for example, in our case, in the region starting from −∞ the density of the fragment with the smallest coefficient, ρ(x), is larger than the other, ρ(x); approaching the A center there is a point in which ρ(x) becomes larger than the other density. This transition between regimes determines both the kinetic peak and the response step. In particular the distance x(1) at which the orbitals ϕ (or the fragment densities, which are simply their square) equateis found to coincide with that of eq in Table , i.e., the maximum of the first kinetic peak as well as of the flex coming from the building up of the response potential step, x(1) = xpeak(1) = xstep(1). Note also that this distance is always somewhere in between the two centers of the fragments, .

Nonetheless, since ρ(x) is asymptotically dominating, by going further in the direction of +∞, the “B regime” is to be encountered again and the two fragment densities, though both very small in magnitude, will be equal again, at some point, x(2):

At this distance also another kinetic peak is appearing and another flex is coming from the exhaustion of the response potential step or, in short, x(2) = xpeak(2) = xstep(2). This is in agreement with the observation in the work of Baerends and co-workers that steps in the response potential and peaks in the kinetic correlation potential are always related.[9,11] Note that this secondary peak is not visible when looking instead at the CCA response potential, again showing the importance of keeping in mind which contribution of the XC potential one is targeting when designing approximations.

Conclusions

In the present work we have generalized the concept of effective and response potentials, as well as of conditional amplitude, for any λ-value, and derived the modulus squared of this latter in the λ → ∞ (SCE) limit. A consistent definition of the response potential in the SCE limit arises from our treatment. In the simple 1D model of a dissociating molecule (eq ), it is found that interesting similarities between dissociation features of the exchange–correlation potential and SCE features, such as the behavior of the co-motion function for increasing internuclear distance or the structure of the SCE response potential itself, can be established. For example, in the dissociation regime, the slope of the co-motion function is determined by the ratio between the ionization potentials of the fragments (compare Figure ), whereby the step height of the exchange–correlation potential is determined by their difference. In addition, the co-motion function confers to the SCE response potential an asymmetric structure which indicates on which side of the system the more electronegative fragment is located.

Further analyzing the different components of the exchange–correlation potential that are relevant in the dissociation limit, namely vresp and v, or v̅resp, we have identified the presence of a second peak of lower intensity in the kinetic correlation potential on the side of the more electronegative atom, and by comparison, we have observed that the peak of the coupling-constant averaged response potential asymptotically coincides with that of the SCE response potential itself. Our work, together with a very recent and promising study,[47] shows that the SCE framework encodes more than few pieces of information on the physical system, and that useful guidelines in the design of highly nonlocal density functional approximations (based on integrals of the density) can fruitfully be drawn from it. A step further in this direction will be to study exact properties of the kinetic potential that appears as the next leading term (∼λ–1/2) in the expansion of the adiabatic connection integrand in the λ → ∞ limit,[36] as well as spin effects that have been shown[57] to enter at orders .

We have also reported, for some small systems (He series, Be, Ne, and H2), the response potential coupling-constant averaged along the adiabatic connection; the study of this different response potential complements that of the response potential at full coupling strength and could provide other hints for the construction of approximate XC functionals, especially of a new generation of DFAs based on local quantities along the adiabatic connection.[46,62,63]