OF-DFT SCF calculations for an H₂ molecule using a nonlocal kinetic energy functional defined on energy coordinate

Abstract

A nonlocal kinetic energy density functional in the orbital-free density functional theory (OF-DFT) is utilized to perform the SCF calculations for an H₂ molecule. The functional is formulated within the framework of the novel DFT [Takahashi, J. Phys. B: At., Mol. Opt. Phys., 2018, 51, 055102] that uses the electron distribution n^e(ε) on the energy coordinate ε as a fundamental variable of DFT. The response function χ^e₀ in the integral kernel of the nonlocal functional is also represented on the coordinate ε. As a notable feature in our approach, χ^e₀ of the whole molecule is given by the composite of the response functions of the isolated fragments constituting the molecule. The electron density and the potential energy curve of the H₂ molecule are computed by the OF-DFT SCF calculations, showing reasonable agreements with those obtained by the Kohn–Sham DFT calculations.

---

1 Introduction

Density-functional theory (DFT) for electrons^1–3 is an indispensable theoretical tool to study the electronic properties of molecules and materials. A number of studies are published annually for the applications of DFT and for the theoretical developments. The theoretical framework of DFT was first provided by the Hohenberg–Kohn (HK) theorem⁴ which proves that the external potential υ_ext applied to an N-electron system uniquely corresponds to the resultant ground state electron density n of the system. The HK theorem offers a theoretical basis to develop the density functional to describe the energy components of the total electronic energy E_tot of the system of interest. Actually, E_tot[n] as a functional of the density n can be given by the sum of the energy contributions, thus,where r represents the coordinate of electrons. The term E_kin[n] on the right-hand side of eqn (1) is the kinetic energy functional, and E_H[n] stands for the Hartree potential of the density n(r). The functional E_xc[n], referred to as the exchange–correlation energy, represents the non-classical interaction between the electrons. The last term provides the Coulomb interaction between the electron density and the external potential. The functional E_tot[n] provides the ground state energy E₀ when the υ-representable density¹n coincides with the density n[υ_ext] that corresponds to the υ_ext of the system of interest.

We note that the existence of the functionals in eqn (1) is guaranteed solely by the HK theorem. Therefore, the range of definition of the functionals must be the set of υ-representative electron densities.¹ However, it can be extended to the set of N-rep. electron densities by virtue of the DFT formalism based on Levy's constraint search.⁵ Thus, the ground state energy E₀ for a given external potential υ_ext can be obtained at least formally through minimization of the energy functional of eqn (1) with respect to N-rep. electron densities on the basis of the variation principle. These discussions guarantee the existence of the route to find the energy E₀ and the corresponding density n using only the electron density as a fundamental variable. Such an approach is, of course, quite advantageous as compared with the wave-function theory because the number of variables to be optimized during the variational calculation can be drastically reduced. However, unfortunately, the accurate kinetic energy functional E_kin in eqn (1) is not currently available, although substantial efforts have been devoted to developing the functional.^6,7 To resolve this problem, a compromise was made by Kohn and Sham,⁸ who introduced a ‘non-interacting reference system’ represented with a set of wave functions {φ_i} for N electrons. Then, the kinetic energy of the non-interacting system with the density n is defined as the sum of the kinetic energies of N electrons, thus,

Accordingly, the energy functional of eqn (1) is rewritten aswhere Ē_xc[n] is supposed to include the difference E_kin–E^non-int_kin in the kinetic energy. In the approach of Kohn–Sham DFT (KS-DFT), the minimization of E_tot[n] with respect to n is replaced by the search in the orbital space {φ_i} as shown in eqn (2). The minimization under the orthonormalization conditions 〈φ_i|φ_j〉 = δ_ij between the orbital pairs (i, j) leads to the set of KS equations,for N electrons. The effective potential [small upsilon, Greek, macron]_eff in eqn (4) is given by the functional differentiation of the last three terms in eqn (3), thusThe exchange–correlation potential [small upsilon, Greek, macron]_xc in eqn (5) plays a decisive role in determining the accuracy of the KS-DFT calculation.

The reliability and accuracy of the KS-DFT have been well established through a lot of applications. However, as noted above, the introduction of the orbitals necessitates calculations to ensure the orthonormalization conditions between all the orbital pairs. This causes a drastic increase in the computational cost with the increase in the number of electrons in the system of interest, which prevents the application of KS-DFT to larger systems. It is, thus, desirable to develop a kinetic energy functional E_kin[n] which enables one to perform DFT calculations without orbitals, that is, orbital-free DFT (OF-DFT) calculations.^6,7,9 In the following paragraphs, we provide a brief review of the development of OF-DFT.

Thomas¹⁰ and Fermi¹¹ provided a formula for the exact kinetic energy for a homogeneous electron gas (HEG) with density n. By using the formula, the Thomas–Fermi (TF) functional E^TF_kin[n] is explicitly written as

with C_TF being . σ in eqn (6) represents the index for spin (σ = α, β). The first OF-DFT functional, referred to as the TF model, is then given byThe TF model can be improved by adding the LDA (local-density approximation) exchange functional E^LDA_x[n] to eqn (7), which gives the Thomas–Fermi–Dirac (TFD) model¹ written asThe LDA functional E^LDA_x[n] is explicitly given bywith the coefficient . Unfortunately, it was proved that neither the TF nor the TFD model realizes the formation of chemical bonds in principle.^12,13 The first primitive approach to resolve this issue was to add the von Weizsäcker (vW) correction¹⁴ to eqn (6). The vW correction term E_vW[n] is given byIt can be readily recognized that E_vW[n] itself gives the exact kinetic energy of a one-electron system or a two-electron system with a closed shell. The TF kinetic energy with the vW correction is simply given by

E^TFW_kin[n] = E^TF_kin[n] + λE_vW[n] (11)

where λ is a weight parameter to mix the vW correction. The analysis of the response function of eqn (11) in momentum space⁹ shows that E^TFW_kin[n] with λ = 1 reproduces the correct behavior at large k limit of the response function of the HEG with density n. On the other hand, λ = 1/9 realizes the correct behavior of the response function at small k. It was shown in ref. 15 that an intermediate value λ = 1/5 gives the best result for atoms with large nuclear charges. It is naturally anticipated that the inclusion of higher order gradient terms leads to better results. It is known, however, that it gives only a slight improvement on the kinetic energy.¹⁶

An effort has also been devoted to develop a series of GGA (generalized gradient approximation) kinetic energy functionals,^17–19 where the inhomogeneity parameter s_σ of the electron density plays a role. The definition of s_σ is common to those in the GGA corrections to E_xc functionals and is given by

It was demonstrated in ref. 19 that several GGA kinetic functionals applied to a set of 77 molecules show rather good agreements with the kinetic energies given by Hartree–Fock (HF) wave functions. Note, however, that the densities used as the inputs to the kinetic functionals were those obtained through HF-SCF calculations at their equilibrium structures. We note that optimizing the density self-consistently using a kinetic energy functional has its own difficulty. The gradient term in eqn (12) tends to flatten the density when the GGA correction term is applied to the density during the SCF procedure. As a consequence, the shell structures inherent in the electron density will be destroyed.

So far, we have reviewed the local or semi-local kinetic energy functionals, which refer only to the density n(r) or gradients of the density. It seems that inclusion of the nonlocal term into the kinetic functional is needed to realize the detailed structure in the electron density of the system of interest. To our knowledge, the nonlocal kinetic energy functional was first developed by Chacón et al. who utilized the weighted average of the density to construct the functional.²⁰ The weight function was determined so that the resultant kinetic functional realizes the response function of the HEG. In a development by Wang and Teter (WT),²¹ they also referred the response function of the HEG to construct the nonlocal term of the kinetic energy functional. The nonlocal term E^nloc_WT[n] in the WT functional is explicitly written as

where ω₀ is the kernel of the integral and only depends on the distance |r − r′|. k_F is the Fermi wave number for the HEG with the density n₀. Explicitly, k_F is determined as k_F = (3π²n₀)^1/3. For the bulk system, the average electron density of the unit cell can be taken as the density n₀. The parameter α in eqn (13) was set at 5/6 in the original study. The kernel ω₀ is directly related to the response function of the HEG with density n₀. Importantly, it was demonstrated in their study²¹ that the shell structure in the density profile of the Ne atom was found to be maintained after the SCF calculation employing the nonlocal term of eqn (13). A lot of developments have been conducted to improve the WT functional.^22–27 In the work by Huang and Carter,²⁴ they introduced a integration kernel ω₀ into eqn (13) that depends on the parameter s_σ in eqn (12) for applications with semiconductors. Interestingly, the functional was also applied to covalent bonds in homonuclear diatomic molecules.²⁸ It is quite impressive that their method is capable of describing the chemical bonds adequately. However, the method for determining the average electron density n₀ in a covalent bond seems somewhat ambiguous. We refer the reader to the literature^6,7 for a more detailed review of the recent advances in the OF-DFT.

In recent works, H. T. developed a series of kinetic energy functionals.^29–32 Importantly, some of the functionals^29,32 are based on the novel framework of the density-functional theory (DFT) for electrons,^33,34 where the electron distribution n^e(ε) defined on the energy coordinate ε plays a fundamental role in the DFT. Hereafter, we refer to the distribution n^e(ε) as the energy electron density. The explicit definitions of n^e(ε) and of the energy coordinate ε will be given in the next section. The external potential υ_ext can also be defined on the energy coordinate ε, that is, υ^e_ext(ε). As proved in the literature,³³ there exists one-to-one correspondence between a subset {n^e(ε)} of the energy electron densities and a subset {υ^e_ext(ε)} of the external potentials completely in parallel to the HK theorem.⁴ Thus, it is possible to formulate a kinetic energy functional E^e_kin[n^e] specific to each subset. The author developed a nonlocal kinetic energy functional²⁹ using the response function χ^e₀(ε, ε′) defined on the energy coordinate for describing covalent bonds in molecules. An important feature in the construction of the response function of a molecule is that it is given as a ‘composite’ of the response functions of its constituent atoms. The functional was applied to the calculation of the potential energy curve of the H₂ molecule. However, we note that the energy values were those obtained solely using a post SCF density by the KS-DFT.

The main purpose of the present work is to devise a method to optimize the electron density self-consistently using the nonlocal kinetic energy functional given in the literature.²⁹ A detailed analyses of the property of the composite response function χ^e,cmp₀(ε, ε′) is also an important issue. It will be shown that the composite response function χ^e,cmp₀(ε, ε′) of H₂ strictly reproduces the full response function χ^e,full₀(ε, ε′) obtained for the whole molecular system. The mechanism underlying the agreement of χ^e,cmp₀(ε, ε′) with χ^e,full₀(ε, ε′) will be discussed. The nonlocal kinetic functional using χ^e,cmp₀(ε, ε′) is applied in the OF-DFT SCF calculation to calculate the potential energy curve of the H₂ molecule. The result will be compared with that obtained by the corresponding KS-DFT calculation. We consider in the future the present method to be applied mainly to the condensed molecular systems such as proteins embedded in water.

The organization of the remainder of this article is as follows. In the next section, we first provide a brief review of the DFT that uses the energy electron density as the fundamental variable. Then, it will be followed by an outline of the construction of the nonlocal kinetic energy functional. The method for updating the electron density during the SCF of the OF-DFT is also described. The subsequent section will be devoted to providing the numerical details to implement the OF-DFT using the real-space grid formalism.^35–40 In the section ‘Results and discussion’ we first compare the composite response function χ^e,cmp₀(ε, ε′) with χ^e,full₀(ε, ε′) for H₂ molecules to show their equivalence. The potential energy curve of H₂ obtained by the SCF calculation of the OF-DFT will be compared with that given by the KS-DFT. We also discuss the advantages and deficiencies of our development in OF-DFT. In the last section, we summarize the present work, and provide a perspective for future developments.

2 Theory and method

This section is devoted to describing the theory and the method for the construction of the nonlocal kinetic energy functional. In the first subsection, the energy electron density n^e(ε), which serves as the fundamental variable in the novel DFT,³³ is introduced. Then, in the next subsection, the nonlocal kinetic energy functional (KEF) is formulated in terms of the functional of n^e(ε). In the third subsection, details of the methodology used to construct the composite response function are addressed. A method to update the electron density in the SCF of the OF-DFT is described in the last subsection.

2.1 Energy electron density

In this subsection, we briefly review the DFT for electrons^33,34 based on the energy electron density n^e(ε). To do this, we first introduce the definition of the function n^e(ε). Suppose that n_σ(r) is the electron density with spin σ of a molecule under consideration. Then, n^e_σ(ε) is given by the projection of n_σ(r) onto the energy coordinate ε,where υ_ext(r) is the external potential felt by the N electrons included in the system of interest. Eqn (14) indicates that the external potential υ_ext is used (though not necessarily) as a function to define the energy coordinate ε for a given spatial coordinate r. Hereafter, in this subsection, let us denote this special potential υ_ext as υ_def to avoid confusion. Provided that contour surfaces of an external potential are completely parallel to those of the potential υ_def(r), then it is possible to identify the potential using the potential defined on the one-dimensional coordinate ε, thus,Specifically, for υ_def(r) itself, the corresponding potential υ^e_def(ε) satisfies the relationship υ^e_def(ε) = ε by definition. Thus, a defining potential υ_def identifies a subset I of external potentials {υ_ext(r)}_I. Using the potential υ_def(r) as the defining potential for ε, all potentials in the subset can be represented by the potential on the coordinate ε without loss of information, that is, {υ^e_ext(ε)}_I. For each subset I of the external potentials, we obtain the corresponding set of ‘ground state’ electron densities {n(r)}_I. Projecting the set of the densities {n(r)}_I using eqn (14) leads to the corresponding set of energy electron densities {n^e(ε)}_I. Quite importantly, it is possible to prove that there exists one-to-one correspondence between the set of potentials {υ^e_ext(ε)}_I and the set of energy electron densities {n^e(ε)}_I.³³ The method of proof is completely parallel to the HK theorem.⁴ Thus, it is possible to formulate an energy functional such as E^e_xc[n^e] with n^e being an argument for each subset I.^33,34 In the literature,^29,32 the author developed kinetic energy functionals E^e_kin[n^e] including nonlocal terms.²⁹ In the following subsection, the outline of the method to construct the nonlocal kinetic energy functional for the energy electron density will be provided.

2.2 Nonlocal kinetic energy functional

An important feature of our development of the nonlocal kinetic energy functional (KEF) is that it refers to the response function of the constituent atoms in the molecule of interest, rather than that of the homogeneous electron gas (HEG). In the following, we first provide the general formulation for the KEF described with the inverse of the response function of the system of interest. To do this, we minimize the energy E_tot[n] in eqn (1) with respect to the density n under constraint , leading to the condition at the stationary state with the ground state density n₀,where μ is a chemical potential, and υ_eff[n](r) is defined as

The differentiation of eqn (16) with respect to n(r) at n₀ gives

where χ₀(r, r′) is the response function of the system with the density n₀ and defined as

Apart from the system of interest, we now introduce a reference system of which the ground state density is given as n₀. Then, for an arbitrary density n close to n₀, its kinetic energy E_kin[n] can be approximated by the second-order Taylor expansion with respect to the deviation δn = n − n₀, thus,

Using the relationship of eqn (18), we getThis equation constitutes the basis for the developments in the following.

The integration kernel of the nonlocal term of eqn (21) is a 6-dimensional function. Therefore, direct evaluation of eqn (21) requires a huge computational cost. However, when referring to the HEG as the reference system, the corresponding response function χ₀ can be reduced to a function of |r − r′| as in the kernel ω₀ in eqn (13). In our approach, we take a different approach based on the DFT that uses n^e(ε) as a fundamental variable. We consider the response function between the energy coordinates ε and ε′, not between the spatial coordinates r and r′. To this end, χ₀(r, r′) is projected onto the energy coordinate to yield χ^e₀(ε, ε′),

It is supposed that the external potential υ_ext(r) of the system of interest is used as the defining potential υ_def(r) in eqn (22). Due to the projection of eqn (22), large amounts of the information contained by the response function χ₀(r, r′) will be lost. However, note that when the contour surfaces of the polarization density δn(r) = n(r) − n₀(r) are parallel to those of the external potential υ_ext(r) used as υ_def(r), the projection does not degrade the energy of the nonlocal term in eqn (21). This can be readily verified as follows.To derive the third equality, we used the assumption that δn(r) is constant over the iso-energy surface of the potential υ_def. The quantity n̄^e(ε) is the value of δn(r) at r specified by ε = υ_def(r). n̄^e(ε) can also be obtained by dividing n^e(ε) by the spatial volume Ω(ε) with the coordinate ε, that is,

n̄^e(ε) = n^e(ε)·Ω(ε)⁻¹ (24)

where Ω(ε) is defined asEqn (23) is strictly valid, for example, when a spherically symmetric electron density of an atom is to be polarized due to a change in its nuclear charge.²⁹ Of course, in general, eqn (23) holds only approximately, and its validity also depends on the choice of the reference density n₀. This issue will be examined in the ‘Results and discussion’ section for H₂ molecules. Anyway, we now reduce the 6-dimensional integration kernel in eqn (21) to the 2-dimensional one; thus,With this kinetic energy, the total energy E^OF-DFT_tot[n] of the OF-DFT is given byFor later references, we also provide the potential υ_tot[n](r) given as the functional derivative of the total energy, thus,

2.3 Composite response function

The next issue is how to construct the response function χ^e₀ in eqn (26) for a reference system. A naive and time-consuming method to obtain χ₀ or χ^e₀ is to perform the inverse KS-DFT method^41–47 for a given reference density n₀. χ₀ is then constructed with the KS wave functions by exploiting the second-order perturbation theory. This approach is robust, but useless because the computational cost is much higher than that for the KS-DFT calculation itself for the whole system. Our strategy to solve this problem is to construct χ₀ as a ‘composition’ of the response functions of the subsystems comprising the entire system.²⁹ In what follows, we illustrate how to construct the composite response function. Suppose that a molecule of interest consists of two subsystems A and B. Then we solve the KS equation separately for each isolated subsystem. The reference electron density n₀(r) of the entire system is then given by the sum of the ground state densities n_A(r) and n_B(r) for the subsystems A and B, respectively, that is,

n₀(r) = n_A(r) + n_B(r) (29)

This method is completely the same as the Harris approach,⁴⁸ where the density of the coupled fragments is approximated by the sum of the isolated fragment densities. There might be another elaborate method to construct the reference density n₀. This issue will be discussed later. For the density n₀ of eqn (29), the corresponding reference system is defined as the system that has n₀ as its ground state density. Assuming that n₀ is non-interacting υ-representable, the KS equation for the reference system is explicitly written asIt is quite important to note in the equation that the external potential is given as a functional υ_ext[n₀] of the density n₀. Hence, υ_ext[n₀] is not equal to the external potential υ_ext of interest in general. For later discussion, we introduce the full response function for the system identified by eqn (30). χ^full₀(r, r′) is then projected onto the energy coordinate to yield χ^e,full₀(ε, ε′) given byNote that the potential υ_ext in eqn (31) is the external potential for the system of interest and not for the reference system. On the other hand, the composite response function χ^e,cmp₀(ε, ε′) defined on the energy coordinate is given by

χ^e,cmp₀(ε, ε′) = χ^e_A(ε, ε′) + χ^e_B(ε, ε′) (32)

In parallel to eqn (31), the individual response function χ^e_X(ε, ε′) for subsystem X = A or B is also given bywhere the response function χ_X(r, r′) is defined as with υ^X_eff being the KS effective potential of the isolated subsystem X.

Our primary concern to the composite response function χ^e,cmp₀(ε, ε′) is whether it works in eqn (26) as an approximation to χ^e,full₀. Note that a response property is nonlocal and, therefore, the correlation between the two subsystems is being incorporated into the ‘real’ response function. From this consideration, seemingly, the construction of the response function shown in eqn (32) does not work correctly. However, surprisingly, as will be demonstrated in the Results and discussion section, the full response function χ^e,full₀ can be safely replaced by χ^e,cmp₀ in describing the chemical bond in H₂ molecules. The mechanism underlying the equivalence of χ^e,cmp₀ to χ^e,full₀ for H₂ will also be elucidated. As shown in eqn (26), the calculation of the energy E^e,(2)_kin[n] requires the inversion of the response function χ^e,cmp₀ or χ^e,full₀. Since the response function has an eigenvector with a null eigenvalue, its inversion requires some devises.²⁹ Details of the numerical procedure for the inversion will be provided in the Computational details section.

2.4 OF-DFT SCF

For a post-SCF electron density given by e.g. KS-DFT calculations, the TFD model (eqn (8)) or a kinetic GGA functional provides rather appropriate kinetic energies. However, the density determined self-consistently with the TFD or a kinetic GGA functional tends to deviate from the KS density. It is well known that the shell structures in the density profiles of atoms disappear in the SCF calculations using the TFD model. The GGA correction to the model does not improve the situation because the gradient term in the functional tends to flatten the shell structures inherent in the density profile. In our preliminary SCF calculations using eqn (26), we find that the full optimization of the electron density n(r) in the real space gives rise to erroneous behaviors in the electron density. More precisely, a sharp peak was found to arise in a confined spatial region during the SCF procedure. We note that the nonlocal term in eqn (26) is always positive with respect to the deviation δn^e. Therefore, the nonlocal term effectively suppresses the deviation of the density n^e from the reference density n^e₀. However, the integration kernel χ^e,−1₀ applies to the projected density difference δn^e(ε). It is thus possible that n(r) deviates seriously from n₀(r) within a confined region of r without changing the amount of δn^e(ε). This constitutes the major drawback of the method using n^e(ε) in the SCF calculations.

To avoid such an unfavorable situation, we impose a constraint on the polarization density δn(r). To be explicit, δn(r) is subject to the constraint

δn(r) = n(r) − n₀(r) = δn^e(ε)|_ε=υdef(r) (34)

Thus, the amount of the density polarization at r is solely specified by its energy coordinate ε = υ_def(r). The update of the density is then performed through the equation as followswhere η is supposed to a small positive real number, and the initial guess of the density is supposed to be the reference density n₀ in eqn (28). υ^e_tot[n](ε) in eqn (35) is given by the projection of the potential υ_tot[n](r) shown in eqn (28). Explicitly, υ^e_tot[n](ε) is given byIt is easy to see that υ^e_tot[n](ε) = const. is achieved when the density converges, i.e., n_i+1(r) ≃ n_i(r). Remember that the initial guess of the density for the SCF is supposed to be n₀ in eqn (29). Hence, the density at the SCF convergence varies depending on the choice of the reference density n₀ taken as the initial guess.

It is helpful to see in more detail the integrand in eqn (35). First, we note that the following equation can be obtained directly from eqn (18), that is,

We also note that υ^e_tot can be expressed as υ^e_tot = υ^e_kin + υ^e_eff. Then the integrand can be formulated asThe notation δυ^e_kin[n] in the above equation means the potential change due to the variation of density from n₀ to n. In deriving the third equality in eqn (38), we used the fact that υ^e_kin[n₀](ε) + υ^e_eff[n₀](ε) = μ₀. The last term in the right hand side of the last equality becomes zero after integrating over ε′. The last equality in eqn (38) shows that the density change δn due to the variation δυ^e_kin of the kinetic potential compensates that given by the variation δυ^e_eff of the effective potential when the SCF convergence is achieved.

3 Computational details

This section provides the numerical details for the implementation of the OF-DFT calculations. The first subsection describes the numerical techniques used in the SCF of the OF-DFT. In our implementation of the OF-DFT, the real-space grid approach is employed to represent the potentials as well as density, where nonlocal pseudopotentials are usually used to describe the electron-nuclei potential. To perform the OF-DFT calculations for the H₂ molecule, a local potential is introduced for the hydrogen atom in the second subsection. In the third subsection, the numerical details for the implementation of the real-space grid method are provided. The fourth and fifth subsections describe the construction of the response function and the method of its inversion, respectively. The last subsection provides the method to construct the energy electron densities and the response functions.

3.1 OF-DFT SCF

In the OF-DFT SCF calculation, the electron density n(r) will be updated using eqn (35), where the non-negativity of the density must be guaranteed. To this end, we introduce a variable ϕ(r), satisfying n(r) = ϕ(r)². Then, eqn (35) is revised to update ϕ(r), thus²⁹,where [small chi, Greek, tilde]^e₀ is defined as²⁹The acceleration factor η in eqn (39) is set in the range 0.03 ≤ η ≤ 0.08 depending on the convergence behavior of the SCF. The criterion for the SCF convergence is specified by the condition −5.0 × 10⁻⁶ ≤ ΔE = E_i+1 − E_i ≤ 0.0 for the energy E^OF-DFT_tot[n] in eqn (27) given in the atomic units.

3.2 Local potential for H₂

As will be described in the following subsection, the present OF-DFT calculations are performed using the real-space grids.^35–39 Usually, in the KS-DFT calculations, the Coulomb potential formed by the nuclei is represented with nonlocal pseudopotentials to avoid steep behaviors in the resultant wave functions. However, to perform the OF-DFT calculations, the external potential υ_ext must be local. Thus, we introduce a local potential for H atoms to represent a chemical bond of the H₂ molecule. As a consequence, the resultant chemical bond is not for a real hydrogen molecule, but for a ‘pseudo’ H₂ molecule. Note, however, that it makes sense because the result given by the OF-DFT can be fully compared with that by the corresponding KS-DFT calculation using the same potential.

The local potential for the pseudo H atom is constructed by replacing the charge of the proton with a Gaussian charge distribution. Explicitly, the delta function δ(r − R) for the charge of a proton placed at R is replaced by a normalized Gaussian, thus,

The exponent α of the charge distribution was determined using the exponent α′ presented in Table I in the literature⁴⁹ which conducted the non-Born–Oppenheimer calculation for the H₂ molecule. The exponent α′ = 21.95 used in the literature for the initial guess of the proton wave function was translated to the exponent α of the density through α = 2 × 21.95. Then, the potential that corresponds to the distribution is given byThis local potential is used to describe the interaction between the electrons and the pseudo hydrogen atom H_A placed at R_A. The same is true for the potential of pseudo hydrogen atom H_B at R_B. The local potential of a pseudo hydrogen molecule H₂ is thus constructed.

3.3 Real-space grid method

The OF-DFT calculations are performed by conducting a Fortran module²⁹ combined with the ‘Vmol’ package^{33,37–40,50–52} based on the real-space grid formalism.^35,36 The grid width h of the real-space cell is set at h = 0.2867869 a.u., which corresponds to the cutoff energy E_c = 60 a.u. of the calculation using the plane waves. Then, h is multiplied by the number N_g = 64 of grids along each axis of the cell, leading to the cell size L = 18.35436 a.u. The H₂ molecule is enclosed in the cubic cell of the size L, where the mass center of the molecule is set at the center of the box and the molecular axis is aligned parallel to the x axis of the real-space cell. Each axis of the cell is uniformly discretized by 64 grids, leading to the grid width h = 0.2867869 a.u. The Hartree energy E_H[n] in eqn (27) is computed using the method in the literature.⁵³ The exchange–correlation energy E_xc[n] in eqn (27) is evaluated using the BLYP functional.^54,55 The KS-DFT calculations are also performed to yield the wave functions which are used to construct the response functions of the reference systems. The Laplacian in the KS equation shown in eqn (4) is represented with the 4th-order finite difference method.^35,36

3.4 Response function

As shown in eqn (31), the full response function χ^e,full₀(ε, ε′) is obtained by projecting χ^full₀(r, r′) onto the energy coordinate. χ^full₀(r, r′) can be built from the wave functions {φ_i0} given as the eigenvectors of the KS equation in eqn (30). Using the second-order perturbation theory (PT2), χ^full₀(r, r′) is given byThe effective potential υ_eff[n₀](r) in eqn (30) can be obtained by the inverse KS-DFT method^44–47 for example. However, in the case of the H₂ molecule, υ_eff[n₀](r) is directly obtained from the relationship;In applying eqn (44) to eqn (30), μ₀ is simply set to 0 since the value of μ₀ does not affect the construction of the response function. Eqn (44) is also used to determine the kinetic potential υ_kin[n₀](r). And the kinetic energy E_kin[n₀] in eqn (26) can be evaluated using n₀(r)^1/2 as the wave function of the H₂ molecule with the density n₀. Note that these methods can only be used for the one-electron system or two-electron system with the closed shell structure. However, for a ‘frozen’ electron density, its kinetic energy and the potential might be obtained adequately utilizing a series of kinetic GGA functionals.^19,32

The composite response function χ^e,cmp₀(ε, ε′) in eqn (32) can also be obtained from the solutions of the individual KS equations for the isolated subsystems H_A and H_B. The number of orbitals employed to construct the response function of H₂ is set to 10 for both χ^e,full₀ and χ^e,cmp₀.

3.5 Inversion of response function

When computing the energy E^e,(2)_kin[n] in eqn (26), it is necessary to invert the response function χ^e₀(ε, ε′). As proved in Appendix A of the literature,²⁹ the projected response function χ^e₀(ε, ε′) is positive semi-definite. Thus, it can be inverted by excluding the eigenvector with a null eigenvalue from the matrix. Furthermore, the eigenvectors with effectively null eigenvalues are also excluded to avoid unphysical oscillations that arise in the resultant kinetic potential. This approach is also known as the truncated singular value decomposition (TSVD)⁵⁶ in the context of the optimized effective potential (OEP) method.^2,57,58 In our present application, it is found that only a single eigenvector is dominant in the spectral decomposition of the response function. Indeed, for the hydrogen molecule with a bond length of 1.4 a.u., the ratio of the second largest eigenvalue to the largest eigenvalue in χ^e,full₀(ε, ε′), defined in eqn (31), is found to be only 3 × 10⁻⁴. Thus, when constructing the inversion of χ^e₀(ε, ε′), only the eigenvector with the largest eigenvalue is used as in the case of the literature.²⁹ The effect of this treatment was fully discussed in Appendix A of the literature.²⁹

3.6 Energy electron density

The energy electron density n^e_σ(ε) defined in eqn (14) is constructed on a natural logarithmic grid. The coordinate ε_k of the kth grid for the energy in the range from ε_min to ε_max is specified bywhere K is the number of grids on the energy axis. For the pseudo H₂ molecule described by the potential shown in eqn (42), we employ the following settings; ε_min = 0.12 a.u., ε_max = 8.3 a.u., and K = 20. The energy coordinate ε of the real-space grid specified with the index M = (l, m, n) is evaluated using the external potential with its opposite sign. Suppose that the coordinate of a grid is given by r_M = (x_l, y_m, z_n), then its energy coordinate ε(r_M) is given by using eqn (42), thus,where R_A and R_B are the coordinates of the hydrogen atoms. To increase the number of sampling points for constructing the distribution n^e_σ(ε), we introduce a double grid. Then, the electron densities on the double grids are evaluated using 4th-order Lagrange interpolations of the coarse grids. The width of the double grid is set at h/5 with h being the width of the coarse grid. However, for the projection of the response function χ₀(r, r′) onto ε and ε′, the double grid is not employed.

4 Results and discussion

In the first subsection of this section, we elucidate the mechanism underlying the equivalence between the full response function χ^e,full₀ and the composite response function χ^e,full₀ for the H₂ molecule. The density profiles of the H₂ molecule obtained by the OF-DFT SCF calculation are shown in the second subsection, where the densities are compared with those given by the corresponding KS-DFT calculations. The potential energy curves of the OF-DFT for the hydrogen molecule are presented in the third subsection and compared with the KS-DFT curve.

4.1 Response function

Before discussing the potential energy curve of H₂ given by the OF-DFT SCF calculation, we examine the composite response function χ^e,cmp₀(ε, ε′) of eqn (32) by comparing it with the full response function χ^e,full₀(ε, ε′) of eqn (31). Seemingly, the correlation between the two subsystems H_A and H_B is not incorporated in the composite response function χ^e,cmp₀(ε, ε′), since it is merely given by the sum of the individual response functions of the isolated subsystems as shown in eqn (32). In Fig. 1(a), the full response function χ^e,full₀(ε, ε′) is presented, and it is compared with the composite response function χ^e,cmp₀(ε, ε′) in Fig. 1(b). Surprisingly enough, the two response functions are shown to be hardly discernible. This means that χ^e,full₀ is essentially the same as χ^e,cmp₀ at least for the H₂ molecule.

Fig. 1 Response functions of the H₂ molecule with bond distance R(H–H) = 1.4 a.u. (a) The full response function χ^e,full₀(ε, ε′) defined in eqn (31). (b) The composite response function χ^e,cmp₀(ε, ε′) in eqn (32).

In fact, the projection of the response function onto the energy coordinate plays an essential role in the fact that χ^e,cmp₀ coincides with χ^e,cmp₀. For an intuitive understanding, we consider a density polarization due to the HOMO–LUMO transition in the full response function χ^full₀(r, r′) of eqn (43). Fig. 2 shows a schematic illustration of the profile of the product φ_HOMO(r) × φ_LUMO(r) in eqn (43) of the H₂ molecule. It can be seen in the figure that the transition from φ_HOMO to φ_LUMO induces a density polarization along the molecular axis. Thus, the correlation between the subsystems H_A and H_B is realized in the response function given by eqn (43). Such a correlation cannot be described with the composite response function. However, interestingly, the product φ_HOMO × φ_LUMO vanishes completely when projected onto the energy coordinate ε, because the contributions from H_A and H_B with opposite signs cancel with each other as illustrated in Fig. 2. One might think that the cancellation caused by the projection is nothing but a sort of degradation. However, it should be noted that the polarization from n₀(r) of eqn (29) to the ground state density n(r) does not include polarization that induces a dipole moment along the molecular axis. Thus, the cancellation of φ_HOMO × φ_LUMO due to the projection does not degrade the response function as long as the polarization from n₀ to n is concerned. On the other hand, polarization due to the transition from φ_HOMO to the orbital given by the linear combination of the two 2s orbitals on the hydrogen atoms may contribute to the variation from n₀ to n. Importantly, such a polarization in the response function does not vanish when it is projected onto the energy coordinate. Furthermore, the polarization can also be described with the composite response function.

Fig. 2 Profile of the product of the HOMO (highest occupied molecular orbital) and the LUMO (lowest unoccupied molecular orbital) of the H₂ molecule. The bold horizontal lines indicate the values of the product on the same energy coordinate ε defined by υ_def(r).

Now we provide mathematical proof of the fact that χ^e,cmp₀(ε, ε′) coincides with χ^e,full₀(ε, ε′). It is important to notice that n₀ is not independent from n_A. n₀ is solely determined from n_A due to the symmetry of the H₂ molecule. Then, χ^e,full₀ can be reformulated as

By projecting eqn (29) onto the energy coordinate, we directly obtain the following relationship,

n^e₀(ε) = n^e_A(ε) + n^e_B(ε) (48)

Since the potential υ_eff(r) specified by υ^e_eff(ε) has the same symmetry as the molecular symmetry of H₂, n^e_A(ε) in eqn (48) is completely identical to n^e_B(ε). Thus, the relationship

n^e₀(ε) = 2n^e_A(ε) (49)

holds. It should also be noted thatholds for the subsystem response functions in H₂ molecules. Then, eqn (47) can be rewritten aswhere we use the relationship of eqn (32). Thus, it is proved that the full response function χ^e,full₀ built from the wave functions of the whole molecule can be replaced by χ^e,cmp₀ given as the composite of the individual response functions χ^e_A and χ^e_B of the isolated subsystems. Of course, eqn (51) holds only when the system under consideration has high symmetry. However, it can be expected that the major part of the correlation between the subsystems can be incorporated into the composite response function even for systems with lower symmetry. Since the present study serves only as a preliminary test calculation using the simplest model system, more detailed analyses especially for hetero-nuclear diatomic molecules must be conducted in the future work. It should also be noted that the computational cost of constructing χ^e,cmp₀ scales only linearly with system size in contrast to the cost of χ^e,full₀.

Fig. 3(a) shows the contour plots of the local potential of the pseudo H₂ molecule introduced in Subsection 3.2. It is compared with the polarization density δn(r) = δn_KS(r) − δn₀(r) shown in Fig. 3(b). It should be noted that the contour lines of δn(r) are almost parallel to those of the external potential of the pseudo H₂ molecule, although the profile of δn along the molecular axis is rather different from that of the potential. Indeed, Fig. 3(b) shows that the polarization can be characterized by an increase in the electron population in the region of large energy coordinate ε and a decrease in the region of small ε. Therefore, it is reasonable to utilize a response function represented with the energy coordinate instead of the real-space coordinate. In fact, the response function χ^e,full₀, or equivalently, χ^e,cmp₀ presented in Fig. 1 shows bipolar behavior along ε in agreement with the profile of δn.

Fig. 3 (a) Contour plot of the electron-nuclei potential for the pseudo H₂ molecule. The potential of each hydrogen atom is given in eqn (42). The value of is plotted. (b) Polarization density δn(r) defined by the KS-DFT density n_KS(r) subtracted by the reference density n₀(r) of eqn (29). Unit of the density is a.u.⁻³. The white dots represent the positions of the hydrogen atoms.

4.2 OF-DFT SCF density

Here we present the electron density n_OF optimized through the OF-DFT SCF procedure. Fig. 4 shows the profile of n_OF along the molecular axis of the pseudo H₂ molecule with a bond length of 1.4 a.u. The reference electron density n₀(r) defined in eqn (29) as well as the density optimized by KS-DFT are also provided to make comparisons. As shown in the figure, after the SCF calculation of OF-DFT starting from the reference density n₀, the electron population on the two hydrogen atoms increases significantly as compared with n₀. The density profile can be directly compared with the density n_KS obtained by the KS-DFT calculation that uses the same exchange–correlation functional E_xc as in the OF-DFT calculation. It is seen in the figure that n_KS also increases on the hydrogen atoms as compared to n₀ in agreement with n_OF, although the peak height of n_KS is lower than that of the density n_OF. However, as a major difference between n_OF and n_KS, the density n_OF exhibits a depression at the center of the hydrogen atoms, which shows a clear contrast to n_KS. The source of the deficiency of the density n_OF can be attributed to the fact that the density polarization is constrained in the direction of the energy coordinate ε. As illustrated in eqn (35) or in eqn (39), in our approach, the amount of the density polarization δn is constrained to be constant on the surface with the same energy coordinate. However, in the spatial region between the two hydrogen atoms, δn remains almost constant as shown in Fig. 3(b), while the energy coordinate varies significantly within this region as in Fig. 3(a). A possible way to avoid the undesirable behavior in the polarization δn(r) is to construct n₀, for example, by means of the frozen density embedding (FDE) method.^6,59,60 In the FDE method, the electron population in the region between the hydrogen atoms will be increased, because the density of each subsystem is to be optimized under the influence of the surrounding environment. Thus, the polarization density δn̄ = δn_KS − δn_FDE will decrease in the region between the hydrogen atoms in contrast to δn = δn_KS − δn₀ plotted in Fig. 3(b). As a result, the description of the polarization δn̄ using the response function χ^e₀(ε, ε′) might be justified.

Fig. 4 Profile of the spin densities n_σ(r) along the molecular axis of the pseudo H₂ molecule with its bond distance being 1.4 a.u. The positions of the two hydrogen atoms are indicated with the arrows. The OF-DFT density is determined self-consistently by using the full response function χ^e,full₀ in eqn (39). The KS-DFT density is also shown for comparisons. The density n₀ is given by eqn (29). The BLYP functional^54,55 is used for the construction of these densities.

We also provide in Fig. 5 the density profile of the pseudo H₂ molecule along the y axis perpendicular to the molecular axis. It is seen in the figure that the density n_OF shows rather good agreement with the density n_KS, although the densities of n_OF at y = ∼±0.6 a.u. are observed to be smaller than those of n_KS by 0.02 a.u.⁻³.

Fig. 5 Profile of the spin densities n_σ(r) on the line parallel to the y axis, provided for the same H₂ molecule as in Fig. 4. The line intersects the x axis at the grid of x = −0.574 a.u. where the density exhibits its maximum value along the x axis. The other notations are the same as in Fig. 4.

4.3 Potential energy curve

The potential energy curves (PECs) of the pseudo H₂ molecule are presented in Fig. 6. The OF-DFT SCF calculations using the full response function χ^e,full₀ of eqn (31) and the composite response function χ^e,cmp₀ of eqn (32) are performed for the production of these PECs. To make comparisons, the PEC by the KS-DFT calculation using the same E_xc functional is also provided. Note that the von Weizsäcker kinetic energy functional E_vW[n] of eqn (10) gives the same potential energy curve as the KS-DFT, since E_vW[n] is equivalent to the KS-DFT in the case of the H₂ molecule. To see the effect of the SCF procedure on the energies of the OF-DFT, the PEC is also computed using eqn (27) with the reference density n₀ as an argument.

Fig. 6 Potential energy curves (PECs) of the pseudo H₂ molecule. ‘OF-DFT(full)’ stands for the OF-DFT SCF calculation with the full response function χ^e,full₀ defined in eqn (31). ‘OF-DFT(cmp)’ means the calculation with the composite response function χ^e,cmp₀ defined in eqn (32). The KS-DFT calculation is performed using the same E_xc functional (BLYP)^54,55 as in the OF-DFT calculations. For the PEC of ‘n₀’, the density n₀ of eqn (29) is used as the argument of the functional E_tot[n] defined in eqn (27). All the PECs are shifted so that their minimums become zero.

First, it should be noted that the PEC given by the OF-DFT with χ^e,full₀ shows an excellent agreement with that of the OF-DFT with χ^e,cmp₀. This is a direct consequence of the fact that χ^e,full₀(ε, ε′) is equivalent to χ^e,cmp₀(ε, ε′) as shown in Fig. 1. As shown in the PEC by the KS-DFT, the pseudo H₂ molecule has the potential minimum at R(H–H) = ∼1.5 a.u., which differs from the fact that the real H₂ has the energy minimum at the bond length of ∼1.4 a.u. It is found that the equilibrium bond length is elongated by ∼0.1 a.u. in the OF-DFT SCF calculations. Thus, it is revealed that the OF-DFT calculation using the response function represented on the energy coordinate ε degrades the shape of the PEC. However, the overall behavior of the PEC is reasonably described with the OF-DFT method. It is possible that some correction to the kinetic energy functional will improve the behavior of the PEC. A possible method for the correction is to improve the reference density n₀ given in the form of eqn (29). Constructing the density of each subsystem under the influence of its surroundings will deliver a better reference density for the initial guess of the interacting system. As noted above, the frozen-density embedding (FDE) approach^59,60 is a promising candidate to build the reference density. Another approach to improve the OF-DFT might be to make a correction to the kinetic energy functional itself. For example, it is reasonable to tune the kinetic energy density at a coordinate r according to the gradient of the external potential υ_ext(r) used as the defining potential for the energy coordinate. These issues will be addressed in future works. In the comparisons with the PEC of the KS-DFT, the behavior of the PEC given by eqn (27) using the non SCF density n₀ as an argument is found to be apparently worse than the PECs given by the self-consistent density. Actually, the equilibrium bond length is found to be further elongated by ∼0.1 a.u. as compared to the result of the OF-DFT SCF calculation. It is thus found that the present OF-DFT SCF procedure works successfully in describing the σ bond in the hydrogen molecule.

It will be useful to analyze the kinetic energy density d_kin(r) to see which spatial region in the molecule gives the larger kinetic energy error. For the H₂ molecule, we define d^vW_kin(r) of the von Weizsäcker functional or equivalently of the KS-DFT, thus,

Note that the integration of the kinetic energy density d^vW_kin(r) leads to the accurate total kinetic energy E^vW_kin[n] of the hydrogen molecule with the density n. For the kinetic energy functional E^e,(2)_kin[n] of eqn (26), given as the present work, the corresponding kinetic energy density d^PW_kin(r) is defined aswhere δn is defined as δn = n − n₀. Then, we define the difference Δd_kin(r) as

Δd_kin(r) = d^PW_kin[n_PW](r) − d^vW_kin[n_vW](r) (54)

where n_PW is the density optimized through the OF-DFT SCF calculation shown in subsection 2.4 and n_vW is that given by the KS-DFT calculation. Δd_kin(r) on the molecular plane is plotted in Fig. 7. In accordance with the trend that the present approach overestimates the electron densities at the atomic core regions as shown in Fig. 4, Δd_kin(r) exhibits sharp positive peaks at the core regions. Note that the larger curvatures of the density n_PW than n_vW at the atomic cores may also contribute to the growth of the peaks. Meanwhile, Δd_kin(r) shows a significant decrease in the region between the atomic cores. This is the direct consequence of the fact that the density n_PW becomes smaller than n_vM in this region. However, interestingly, it is shown that the depression of Δd_kin(r) is alleviated near the center of the chemical bond. This can be attributed to the increase in curvature of the density n_PW at the center of the bond as compared to that of n_vW. Anyway, it is revealed that the cancellation of Δd_kin(r) takes place between the overestimations of Δd_kin at the atomic cores and the underestimation in the region of the chemical bond.

Fig. 7 The difference Δd_kin(r) of the kinetic energy difference defined in eqn (54) for the H₂ molecule placed on the xy plane with R(H–H) = 1.5 a.u. The white dots on the contour plot represent the positions of the hydrogen atoms.

In closing this subsection, we make a remark on the accuracy of the present approach for the atomic kinetic energy. In Fig. 3 of ref. 29, kinetic energies of a pseudo hydrogen atom with respect to the variation of its nuclear charge Z_v were compared with those given by the KS-DFT. It was shown in the figure that the graph given by the present approach agrees well with that by the KS-DFT over the range of the core charge Z_v investigated, though the deviation becomes larger in the small value of Z_v. It is of substantial interest to compare the OF-DFT kinetic energies of atoms with the accurate reference values obtained e.g. by sophisticated molecular orbitals theories. However, unfortunately, the present calculations utilized the real-space grids to represent the density and the potentials, which necessitates the use of the pseudopotentials to ensure the smooth behaviors of the resultant densities on the potentials. This prevents the direct comparisons of the present method with some reference benchmark data. It should also be noted that the kinetic functional of eqn (26) will violate some physical constraint to be satisfied by the rigorous functionals. As an example of the constraint, it is known that the von Weizsäcker kinetic energy functional E_vW[n] serves as the lower bound of the non-interacting kinetic energy functional E_kin[n],⁷ that is, E_vW[n] ≤ E_kin[n]. Actually, for the hydrogen molecule with R(H–H) = 1.4 a.u., the kinetic energy of eqn (26) is evaluated as 0.9181 a.u., while E_vW[n] is computed as 1.0856 a.u. for the density n optimized by the present OF-DFT SCF calculation. Thus, it was revealed that the present kinetic functional violates the lower bound constraint. It is also possible in general that the kinetic potential becomes negative since the functional was not designed so that it fulfills these desired physical constraints. These unfavorable properties should be eliminated from the functional in the future developments.

5 Conclusions

In this article, a nonlocal kinetic energy functional E^nloc_kin was applied to the OF-DFT SCF calculation to describe the covalent bond in the H₂ molecule. An important feature of the functional E^nloc_kin is that it is based on a novel framework of DFT where the electron distribution on the energy coordinate serves as a fundamental variable of the DFT.^29,33,34 Thus, the response function χ^e₀(ε, ε′) as well as the polarization density δn^e(ε) in the nonlocal functional are described on the energy coordinate ε.

Another notable feature of the present approach is that the response function χ^e₀(ε, ε′) of the H₂ molecule is given as the sum of the response functions of the isolated subsystems. Surprisingly, it was demonstrated for the H₂ molecule that the response function χ^e,full₀ of the fully correlated system with a reference density n₀ can be faithfully reproduced by the composite response function χ^e,cmp₀ built from the response functions of the subsystems. The underlying mechanism for the equivalence between χ^e,full₀ and χ^e,cmp₀ was discussed, and a mathematical proof for the fact was provided with respect to the H₂ molecule.

A method to update the electron density for the SCF of the OF-DFT was also provided. The key of the method is to impose a constraint on the polarization density δn(r). That is, δn(r) is constrained to be constant on the surface with the same energy coordinate ε. In other words, the density polarization is considered only along the energy coordinate ε. The SCF density of the OF-DFT showed a reasonable agreement with the density given by the KS-DFT for the H₂ molecule although a significant density depression was observed in the region between the two hydrogen atoms. The potential energy curves (PECs) of the hydrogen molecule were also computed using the SCF densities of the OF-DFT. Although the equilibrium bond length was found to be elongated by ∼0.1 a.u. compared with the PEC by the KS-DFT, it was demonstrated that the SCF density can adequately describe the covalent bond. It was also found that the PEC can be significantly improved by the optimization of the density through the SCF that employed the reference density n₀ as an initial guess.

It is possible that the results given by the present approach can be improved e.g. by modifying the method to construct the reference density n₀. It might also be effective to modulate the kinetic energy density at a coordinate r according to the gradient of the external potential υ_ext(r). These issues will be examined in future works.

Finally, we also remark on the computational cost of the OF-DFT. Unfortunately, at the present stage of development, the computational time for the OF-DFT is found to be much larger than that for the KS-DFT. The major computational cost in our approach of the OF-DFT arises from the construction of the energy electron density and the 2-dimensional response function χ^e₀(ε, ε′). Note, however, that the cost associated with these calculations increases only linearly with the size of the system of interest. And we also expect that these calculations are suitable for the application of GPGPU (General-Purpose computing on Graphics Processing Units) and thus, they can be accelerated effectively.

Author contributions

Hideaki Takahashi: conceptualization; formal analysis; funding acquisition; investigation; methodology; software; visualization; writing – original draft; writing – review and editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information provides the raw data used to produce the graphs in Fig. 1(a) and (b). See DOI: https://doi.org/10.1039/d5cp04133c.

Acknowledgements

This paper was supported by the Grant-in-Aid for Scientific Research(C) (no. 17K05138, no. 22K12055, and no. 25K08570) from the Japan Society for the Promotion of Science (JSPS) and the Grant-in-Aid for Challenging Exploratory Research (no. 25620004) from the Japan Society for the Promotion of Science (JSPS).

Notes and references

1. R. G. Parr and W. Yang, Density-functional theory of atoms and molecules, Oxford university press, New York, 1989.
2. R. M. Martin, Electronic Structure, Basic Theory and Practical Methods, Cambridge University Press, Cambridge, 2004.
3. F. Jensen, Introduction to Computational Chemistry, Wiley, 2009.
4. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
5. M. Levy, Proc. Natl. Acad. Sci. U. S. A., 1979, 76, 6062–6065.
6. in Recent Progress in Orbital Free Density Functional Theory (Recent Advances in Computational Chemistry), ed. T. A. Wesolowski and Y. A. Wang, World Scientific, Singapore, 2013.
7. W. Mi, K. Luo, S. B. Trickey and M. Pavanello, Chem. Rev., 2023, 123, 12039–12104.
8. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
9. Y. A. Wang and E. A. Carter, Theoretical Methods in Condensed Phase Chemistry (Progress in Theoretical Chemistry and Physics), Kluwer, Dordrecht, 2000, ch. 5, vol. 5.
10. L. H. Thomas, Proc. Cambridge Philos. Soc., 1927, 23, 541–548.
11. E. Fermi, Z. Phys., 1928, 48, 73–79.
12. E. Teller, Rev. Mod. Phys., 1962, 34, 627–631.
13. N. L. Balàzs, Phys. Rev., 1967, 156, 42–47.
14. C. F. von Weizsacker, Z. Phys., 1935, 96, 431–458.
15. K. Yonei and Y. Tomishima, J. Phys. Soc. Jpn., 1965, 20, 1051–1057.
16. J. P. Perdew, M. Levy, G. S. Painter, S. Wei and J. B. Lagowski, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 838–843.
17. H. Lee, C. Lee and R. G. Parr, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 44, 768–771.
18. F. Tran and T. A. Wesołowski, Int. J. Quantum Chem., 2002, 89, 441–446.
19. A. J. Thakkar, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 46, 6920–6924.
20. E. Chacón, J. E. Alvarellos and P. Tarazona, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 32, 7868–7877.
21. L.-W. Wang and M. P. Teter, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13196–13220.
22. Y. A. Wang, N. Govind and E. A. Carter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 13465–13471.
23. Y. A. Wang, N. Govind and E. A. Carter, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 60, 16350–16358.
24. C. Huang and E. A. Carter, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 1–15.
25. L. A. Constantin, E. Fabiano and F. Della Sala, Phys. Rev. B, 2018, 97, 1–10.
26. W. Mi, A. Genova and M. Pavanello, J. Chem. Phys., 2018, 148, 184107.
27. W. Mi and M. Pavanello, Phys. Rev. B, 2019, 100, 1–6.
28. J. Xia, C. Huang, I. Shin and E. A. Carter, J. Chem. Phys., 2012, 136, 084102.
29. H. Takahashi, Int. J. Quantum Chem., 2022, 122, e26969.
30. H. Takahashi, J. Chem. Phys., 2023, 158, 014102.
31. H. Takahashi, J. Chem. Phys., 2023, 159, 124118.
32. H. Takahashi, Electron. Struct., 2025, 7, 025003.
33. H. Takahashi, J. Phys. B: At., Mol. Opt. Phys., 2018, 51, 055102.
34. H. Takahashi, J. Phys. B: At., Mol. Opt. Phys., 2020, 53, 245101.
35. J. R. Chelikowsky, N. Troullier and Y. Saad, Phys. Rev. Lett., 1994, 72, 1240–1243.
36. J. R. Chelikowsky, N. Troullier, K. Wu and Y. Saad, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 11355–11364.
37. H. Takahashi, T. Hori, T. Wakabayashi and T. Nitta, Chem. Lett., 2000, 222–223.
38. H. Takahashi, T. Hori, T. Wakabayashi and T. Nitta, J. Phys. Chem. A, 2001, 105, 4351.
39. H. Takahashi, T. Hori, H. Hashimoto and T. Nitta, J. Comput. Chem., 2001, 22, 1252–1261.
40. H. Takahashi, N. Matubayasi, M. Nakahara and T. Nitta, J. Chem. Phys., 2004, 121, 3989–3999.
41. W. M. C. Foulkes and R. Haydock, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 12520–12536.
42. R. van Leeuwen and E. J. Baerends, Phys. Rev. A: At., Mol., Opt. Phys., 1994, 49, 2421–2431.
43. Q. Wu and W. Yang, J. Chem. Phys., 2003, 118, 2498–2509.
44. E. S. Kadantsev and M. J. Stott, Phys. Rev. A: At., Mol., Opt. Phys., 2004, 69, 012502.
45. B. Kanungo, P. M. Zimmerman and V. Gavini, Nat. Commun., 2019, 10, 4497.
46. Y. Shi and A. Wasserman, J. Phys. Chem. Lett., 2021, 12, 5308–5318.
47. H. Takahashi, J. Chem. Phys., 2024, 161, 104108.
48. J. Harris, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 1770–1779.
49. Y. Shigeta, H. Takahashi, S. Yamanaka, M. Mitani, H. Nagao and K. Yamaguchi, Int. J. Quantum Chem., 1998, 70, 659–669.
50. H. Takahashi, H. Hashimoto and T. Nitta, J. Chem. Phys., 2003, 119, 7964–7971.
51. T. Hori, H. Takahashi and T. Nitta, J. Chem. Phys., 2003, 119, 8492–8499.
52. H. Takahashi, S. Sakuraba and A. Morita, J. Chem. Inf. Model., 2020, 60, 1376–1389.
53. R. N. Barnett and U. Landman, J. Chem. Phys., 1993, 48, 2081–2097.
54. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
55. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
56. A. J. C. Felipe, A. Bulat, T. Heaton-Burgess and W. Yang, J. Chem. Phys., 2007, 127, 174101.
57. R. T. Sharp and G. K. Horton, Phys. Rev., 1953, 90, 317.
58. J. D. Talman and W. F. Shadwick, Phys. Rev. A: At., Mol., Opt. Phys., 1976, 14, 36–40.
59. T. A. Wesolowsk and A. Warshel, J. Phys. Chem., 1993, 97, 8050–8053.
60. T. A. Wesolowski, S. Shedge and X. Zhou, Chem. Rev., 2015, 115, 5891–5928.

---