Electronic structure calculations based on Kohn–Sham density functional theory (KSDFT) that incorporate exact-exchange or hybrid functionals are associated with a large computational expense, a consequence of the inherent cubic scaling bottleneck and large associated prefactor, which limits the length and time scales that can be accessed. Although orbital-free density functional theory (OFDFT) calculations scale linearly with system size and are associated with a significantly smaller prefactor, they are limited by the absence of accurate density-dependent kinetic energy functionals. Therefore, the development of accurate density-dependent kinetic energy functionals is important for OFDFT calculations of large realistic systems. To this end, we propose a method to train kinetic energy functional models at the exact-exchange level of theory by using a dictionary of physically relevant terms that have been proposed in the literature in conjunction with linear or nonlinear regression methods to obtain the fitting coefficients. For our dictionary, we use a gradient expansion of the kinetic energy nonlocal models proposed in the literature and their nonlinear combinations, such as a model that incorporates spatial correlations between higher order derivatives of electron density at two points. The predictive capabilities of these models are assessed by using a variety of model one-dimensional (1D) systems that exhibit diverse bonding characteristics, such as a chain of eight hydrogens, LiF, LiH, C4H2, C4N2, and C3O2. We show that by using the data from model 1D KSDFT calculations performed using the exact-exchange functional for only a few neutral structures, it is possible to generate models with high accuracy for charged systems and electron and kinetic energy densities during self-consistent field iterations. In addition, we show that it is possible to learn both the orbital dependent terms, i.e., the kinetic energy and the exact-exchange energy, and models that incorporate additional nonlinearities in spatial correlations, such as a quadratic model, are needed to capture subtle features of the kinetic energy density that are present in exact-exchange-based KSDFT calculations.

Kohn–Sham density functional theory (KSDFT)1,2 is one of the most popular and powerful techniques to perform electronic structure calculations,3–6 with recent advances allowing for simulations of systems containing thousands of atoms.7–9 The main limitations of KSDFT are its cubic scaling and a large prefactor associated with the underlying electronic-eigenvalue problem due to which it cannot be applied to many large systems that are encountered in various materials science and physics applications. In contrast, orbital-free density functional theory (OFDFT) scales linearly with the number of atoms and has a significantly smaller prefactor, which opens the possibility for its application to systems containing as many as a million atoms,10,11 thereby enabling us to handle realistic systems that are currently not accessible in KSDFT simulations.12–14 However, as discussed in the literature (e.g., Refs. 15–19), one of the main challenges in OFDFT calculations is that the existing kinetic energy functionals are not accurate enough to capture subtle features of the kinetic energy density obtained from a KSDFT calculation.

Early efforts to derive such a functional focused primarily on using gradient expansions of the kinetic energy,20–25 but such expansions are typically good only under limiting conditions (for example, homogeneous systems). Therefore, for many years, researchers have focused on approximating the kinetic energy functional by using a nonlocal kernel to capture spatial correlations in the density, along with the terms present in gradient expansion based models.26–31,33–36 However, these nonlocal kernels are also valid under limiting conditions (for example, homogeneous periodic systems), meaning that models containing these nonlocal correlations may not be accurate for general systems. Many researchers have used a different route to design kinetic energy functionals by using the generalized gradient approximation (GGA).37,38 For example, Constantin et al. showed that a properly designed Laplacian-level kinetic energy functional can have similar accuracy to that of a nonlocal model.40,41 Systematic comparisons of many such models have been performed in Refs. 32, 38, and 39, and these results suggest that the true kinetic energy functional is perhaps very complex and must include GGA-type functional forms, nonlocal terms, and their higher order generalizations.27,42

In the last few years, many researchers have used machine learning techniques to learn kinetic energy functionals directly from KSDFT data.43–50 These efforts have resulted in good approximations, but one of the fundamental limitations of data-driven methods is the transferability to systems outside the training domain. A hybrid model that combines the data-driven learning and physics-based model parameters can be a promising path forward toward developing an accurate and transferable kinetic energy functional model. In line with this philosophy, in Ref. 15, we developed a simple methodology to learn a kinetic energy functional by using a dictionary of physically motivated terms, and the coefficients were obtained by using linear regression. Among the multiple models explored in that study, the model containing spatial correlations between electron densities and between electron densities and their derivatives (using two different kernels and, hereafter, denoted as the Q1 + Q3 model) was optimal both in terms of number of parameters and in terms of performance.

In this contribution, we expand the scope of the work proposed in Ref. 15 and assess the performance of the models of kinetic energy functionals generated by using data from one-dimensional model KSDFT calculations performed with the exact-exchange functional. Figure 1 shows electron densities and kinetic energy densities of two representative one-dimensional systems, H8 and C4N2, obtained from one-dimensional model KSDFT calculations performed by using the local density approximation (LDA) exchange and exact-exchange functionals. It is clear that the electron densities and kinetic energy densities obtained from exact-exchange-based calculations exhibit features that are absent in LDA exchange-based calculations. Here, we note that these subtle features are related to the fact that the exchange operator is full-rank because it incorporates nonlocal correlations between different wavefunctions of the Hamiltonian. Thus, it is natural to ask if the models proposed in Ref. 15 can reliably capture subtle features present in the kinetic energy density profiles. In addition, a density functional model for kinetic energy alone is not sufficient to run OFDFT calculations with the exact-exchange level of accuracy. This is because the exact-exchange energy itself depends on the wavefunctions. Hence, to run an OFDFT calculation with the exact-exchange level of accuracy, we need a density functional of both the kinetic energy and the exact-exchange energy. In that regard, the second objective of this study is to explore the effectiveness of the proposed models to describe a modified kinetic energy functional, which is obtained by the sum of non-interacting kinetic energy and the exact-exchange energy minus the LDA exchange energy.

FIG. 1.

Plots showing the difference in ground state electron densities and kinetic energy densities of two representative one-dimensional (1D) systems, H8 [(a) and (b)] and C4N2 [(c) and (d)], obtained from model KSDFT calculations using local density approximation (LDA) exchange and exact-exchange functionals.

FIG. 1.

Plots showing the difference in ground state electron densities and kinetic energy densities of two representative one-dimensional (1D) systems, H8 [(a) and (b)] and C4N2 [(c) and (d)], obtained from model KSDFT calculations using local density approximation (LDA) exchange and exact-exchange functionals.

Close modal

Another goal of this study is to explore the possibility of using these empirical models to train generic or universal kinetic energy models for a large number of heterogeneous systems. Such a generic model will allow us to perform large-scale OFDFT calculations. For example, Fig. 2 shows that a generic model trained by using structures from a variety of systems (further details are given in Sec. III A) using the model (Q1 + Q3) proposed in our previous work (Ref. 15) is able to capture the important features of kinetic energy density, but the R2 score (i.e., the coefficient of determination) is only 0.985, and deviations from the ground truth are clearly evident. To this end, it is instructive to explore different ways to enhance the predictive capabilities of such generic models.

FIG. 2.

Kinetic energy densities predicted for a representative one-dimensional system, C4N2, using a model that incorporates spatial correlations of the electron density and its derivative (denoted by Q1 + Q3 in Ref. 15) are compared with exact-exchange-based calculations. The insets show the deviations of the model from direct calculations. Predictions made by using a quadratic model (see details in Sec. II) proposed in this work are in better agreement with the KSDFT results.

FIG. 2.

Kinetic energy densities predicted for a representative one-dimensional system, C4N2, using a model that incorporates spatial correlations of the electron density and its derivative (denoted by Q1 + Q3 in Ref. 15) are compared with exact-exchange-based calculations. The insets show the deviations of the model from direct calculations. Predictions made by using a quadratic model (see details in Sec. II) proposed in this work are in better agreement with the KSDFT results.

Close modal

With these objectives in mind, we propose physics-inspired empirical models of the kinetic energy functional by including spatial correlations between higher order derivatives of the electron density. As discussed in Sec. II, we start with a dictionary of functional forms that have been proposed in the literature to capture the local and nonlocal contributions and use a basis function expansion to model the nonlocal kernel. To illustrate the effectiveness of our models, in Sec. III, we use a variety of heterogeneous one-dimensional systems, such as LiF, LiH, a one-dimensional chain of eight hydrogens, C4H2, C4N2, and C3O2. These systems exhibit a variety of bonding characteristics, such as the presence of single, double, and triple bonds, and form a good test set to explore the capabilities of the models. We show that reliable density-dependent models for the kinetic energy can be trained by using data from only a few structures.

To model the kinetic energy obtained from solving the model one-dimensional Kohn–Sham equations with the exact-exchange functional, we consider the nonlocal models proposed in Ref. 15 and also generalize these models to incorporate spatial correlations and higher order nonlinearities. Details of the model one-dimensional Kohn–Sham equations for exact-exchange are given in  Appendix A.

Equation (1) combines all the models proposed in Ref. 15. It contains a variety of local terms and a few nonlocal terms. The local terms in the model are the first three terms in the gradient expansion of the kinetic energy functional for a one-dimensional system derived by Šamaj and Percus.20 To model the nonlocal contributions, three different nonlocal kernels are used, and these nonlocal kernels are modeled by using a set of basis functions,

tx=tlocalx,c+j=1NGd1,j[ρ(x)]2eβ1,jxx2u(x,x)2[ρ(x)]2dx+j=1NGd2,j[ρ(x)]2eβ2,jxx2u(x,x)2ρ(x)ρ(x)2dx+j=1NGd3,jxρ(x)eβ3,jxx2u(x,x)2xρ(x)dx,
tlocalx,c=c1ρ3(x)+c2[ρ(x)]2ρ(x)+c32ρ(x),
(1)
u(x,x)=ρ(x)+ρ(x).

Here, cj and di,j are the coefficients of the model, βi,j are the widths of the Gaussian functions, and NG is the number of Gaussian basis functions used to model the nonlocal kernel. The three terms present in the local kinetic energy expression, tlocal, are based on the gradient expansion and satisfy the uniform coordinate scaling constraint.51–54 The kernels in the nonlocal terms are dimensionless quantities, and the exponents of density and gradient of density have been selected to satisfy coordinate scaling requirements. The optimal set of coefficients is obtained using linear regression by fitting the model to kinetic energy density data obtained by solving the eigenvalue problem of the Kohn–Sham Hamiltonian using the exact-exchange functional. A detailed description of the numerical implementation of model KSDFT with the exact-exchange functional is given in  Appendix B. To calculate the coefficients, we optimize the following objective function:

Ω=tKSDFTAc2,
(2)

where tKSDFT is the vector that contains information about kinetic energy densities at different grid points, c = {c1, c2, …, d1,1, d1,2, …} is the vector that contains the unknown coefficients, and A is the design matrix such that the number of rows is equal to the number of grid points used for fitting (i.e., length of tKSDFT, and the number of columns is equal to the number of terms used to model the kinetic energy functional. Therefore, each row of A contains values of the terms shown in Eq. (1) at a grid point, and for a given set of values of the widths of the Gaussian functions, the vector c is obtained by solving a linear system. The widths of the Gaussian kernels are optimized by using the Nelder–Mead method.

In Secs. II AII C, we propose a few other nonlocal models to approximate the kinetic energy density. These models are different from the models discussed above and are also inspired by physics-based models proposed in the literature. However, parameters present in these new nonlocal models cannot be obtained by using a linear regression framework, and a nonlinear minimization method has to be employed to learn the coefficients.

Chacón et al.26 introduced the notion of a locally averaged density, ρ̄r, to incorporate spatial correlations into the kinetic energy functional,

ρ̄(r)=ρ(r)Wrr,ρrdrandWr,ρrdr=1.
(3)

Using this averaged electron density, the authors proposed the following form to model the kinetic energy:

Tρ=TvWρ13TLDAρ+43ρrtρ̄rdr.
(4)

Here, TvWρ is the von Weizsäcker term, TLDAρ is the Thomas–Fermi term of the kinetic energy functional, and t is the function of density such that TLDAρ=ρrtρ(r)dr. Intuitively, this model can be generalized by using models proposed in Ref. 15 for one-dimensional systems. For example, if we focus on the local terms [Eq. (1)], then a probable model for the kinetic energy is

Tρ=tlocalρx,c+t̄localρ̄x,ρx,c̄dx,t̄localρ̄x,ρx,c̄=c̄2ρ̄3x+c̄2ρ̄xρxρx+c̄32ρ̄x,ρ̄(x)=ρ(x)p=1Ngdpexx2ρ(x)+ρ(x)2dxp=1Ngdpexx2ρ(x)+ρ(x)2dx,
(5)

where dp for p = 1, 2, 3, … are the coefficients used to model the nonlocal kernel. Our preliminary analysis suggests that this model is, indeed, very powerful, but due to the nonlinearities embedded in the higher order polynomial terms, such as ρ̄2x,ρ̄3x, obtaining the best set of fitting parameters becomes challenging. To this end, we simplify this model and consider local averaging of only a few quantities, ρx and ρx, and consider the following nonlocal model:

tnl4x=tlocalρx,c+α1[ρx]2Q4xx,ux,x[ρx]2dx+α2xρxQ4xx,ux,xxρxdx,

and

ux,x=ρx+ρx.
(6)

Here, α1 and α2 are the two parameters introduced to make the two terms dimensionally consistent. We note that the exponent of 2 in the first nonlocal term is required to satisfy the coordinate scaling criterion. The model described above is similar to the model containing two nonlocal kernels (with spatial correlations in the electron density and the derivative of the electron density) as discussed in Ref. 15, except that in Eq. (6), the same kernel is used to capture spatial correlations between the electron density and its derivative. We model the nonlocal kernel Q4 by using a set of basis functions,

Q4xx,ux,x=k=1NGdkeβkxx2ux,x2.
(7)

Interestingly, the model is linear with respect to all the coefficients; however, there is a coupling between the coefficients dk and the parameters α1, α2. Thus, we adopt a two-step iterative optimization scheme to obtain the optimal set of coefficients. In the first step, we use linear regression to solve for c and dk while keeping the parameters α1, α2 constant. In the second step, we fix dk and solve another linear system to obtain c and α1, α2. This two-step scheme is repeated until convergence.

The model proposed in Sec. II A includes spatial correlations between electron densities and between derivatives of the electron densities at two points. To further understand the relevance of spatial correlations between higher order derivatives of the electron density, we consider a model that is quadratic in terms of tlocalρx,c,

tnl5x=tlocalρx,c+tlocalρx,c̄×Q5xx,u(x,x)tlocalρx,c̄dx,tlocalρ(x),c=c1ρ3(x)+c2[ρ(x)]2ρ(x)+c32ρ(x).
(8)

Here, Q5 is the nonlocal kernel that is modeled using a set of Gaussian basis functions as discussed in Sec. II A, and tlocalρx,c contains terms that depend on the density and its derivatives at x. We also note that the nonlocal part of tnl5x does not satisfy the coordinate scaling law; however, the local part of the model satisfies it. Therefore, to design a quadratic model that satisfies the coordinate scaling condition, we introduce the dimensionless variables p and q that are commonly used in GGA models of the kinetic energy functional.37,38 The expression of the modified quadratic model in 1D is

tnl5,scalingx=tlocalρx,c+t̄ρx,c̄×Q5scalingxx,u(x,x)t̄ρx,c̄dx,tlocalρ(x),c=c1ρ3(x)+c2[ρ(x)]2ρ(x)+c32ρ(x),t̄ρ(x),c̄=1+c̄1p(x)+c̄2q(x)ρ2(x),p(x)=|ρ(x)|2ρ4(x),q(x)=2ρ(x)ρ3(x).
(9)

We note that the functional form of both tnl5x and tnl5,scalingx is quadratic c̄, and we adopt a two-step optimization scheme to obtain an optimal set of coefficients: In the first step, we initialize c̄ with random numbers and solve the linear system in c and d. In the second step, we optimize the objective function with respect to c̄ by using a trust-region algorithm. These two steps are repeated until convergence.

Liu and Parr55 argued that the Kohn–Sham kinetic energy can be expressed as a functional of degree one in electron density,

Tρ=18r2γ1r,rdr=18rρr2ρrdr+12ρrrrPr,rr=r+rρrrPr,rr=rdr,

where

P2r,r=γ1r,r2ρrρr
(10)

and γ1r,r is the first order density matrix associated with the ground state electron density. The authors further suggested that Pr,r is a function of degree zero in ρr, but to construct nonlocal models, we relax this constraint and model the function Pr,r by using a set of Gaussian functions. Therefore, we use (in 1D)

Px,x=k=1NGdkeβkxx2ux,x2,ux,x=ρx+ρx,
(11)

where dk are the coefficients and βk are the widths of the Gaussian functions used to model the function Px,x.

It is easy to see that after the derivatives of Px,x are evaluated and substituted in Eq. (10), imposing the condition x = x′ reduces the kinetic energy to a model that depends on only electron density and its derivatives at x. Therefore, to construct a nonlocal model, we do not impose this constraint and instead integrate over x′ to incorporate spatial correlations in the electron density and its derivatives. This leads to the following nonlocal model:

tnl6x=tlocalρx,c+ρxxxPx,x+xρxxPx,xPx,xdxdx.
(12)

An additional term in the denominator is introduced to preserve the overall dimensionality of the terms on the right-hand side of Eq. (10), and this normalization introduces additional nonlinearity into the model. The local terms are the same as those in the quadratic model. The analytic expressions for the first and second order derivatives of P(x, x′) are given in  Appendix C. To evaluate the coefficients and the widths of the Gaussian functions, we use a version of the Nelder–Mead algorithm implemented in MATLAB.

We use six different one-dimensional systems, namely, H8, LiH, LiF, C4N2, C4H2, and C3O2, to assess the performance of our models. Model one-dimensional Kohn–Sham calculations (using exact-exchange) are performed for all of these six systems with open boundary conditions, and kinetic energy densities and electron densities for 30 different structures for each of these six systems are generated for training and testing purposes. For each system, three structures are randomly chosen to train models, and the remaining 27 structures are used for testing the trained model. We use six Gaussian functions to approximate the nonlocal kernels in tnl4x, tnl5x, tnl5,scalingx, and tnl6x. Combining local and nonlocal terms, there are nine parameters in the model tnl4, 12 parameters in the models tnl5 and tnl5,scalingx, and nine parameters in the model tnl6. To assess the performance of the models proposed in Ref. 15, for exact-exchange calculations, we use the same nomenclature, i.e., the model containing the local terms and the first nonlocal term in Eq. (1) (the second term on the right-hand side) is denoted by Q1. Similarly, the model containing the local terms and the third (fourth) term in Eq. (1) is denoted by Q2 (Q3).

The performance of the eight models proposed in Ref. 15 is measured in terms of the coefficient of determination (i.e., R2 score) and root mean square error (RMSE), and the maximum absolute error in the predicted total kinetic energy is shown in the supplementary material. Interestingly, even for data obtained from KSDFT calculations with the exact-exchange functional, kinetic energy models developed for each system have R2 scores larger than 0.99, and maximum absolute errors in the predicted kinetic energies are less than 1% of the KSDFT values. For all systems considered in this study, it is clear that inclusion of nonlocal terms leads to improved performance. For example, models containing one nonlocal term have lower RMSE scores than the local models. Similarly, models containing two nonlocal terms have lower RMSE scores than models containing one nonlocal term. However, we do not see an improvement in the predictive capability when all the three nonlocal terms are present in the models. Finally, we find that the models containing Q1 + Q3 nonlocal kernels perform better than the other models for almost all systems. For some systems, Q1 + Q2 + Q3 performs better than Q1 + Q3, but it also contains a significantly larger number of coefficients compared to Q1 + Q3 and is prone to over-fitting.

In Table I, the performance of the three nonlocal models introduced in this work is summarized. We see that for most of the systems listed in Table I, kinetic energies predicted using the quadratic models have the lowest RMSE scores and lowest maximum absolute errors. For many systems, RMSE scores of the models generated by using tnl4 are similar to those of the quadratic models. However, models generated by using tnl6 are unable to properly learn the kinetic energy density. We also observe that the quadratic model, tnl5, performs best in terms of the error in the total kinetic energy among the three models. This is interesting because tnl5 does not satisfy the coordinate scaling condition, and the quadratic model that satisfies the scaling condition has much higher RMSE values. Therefore, for the subsequent analyses, we do not consider tnl5,scalingx. Here, we note that the quadratic models have 12 fitting parameters compared to nine in the case of tnl6, which can potentially make tnl6 a bit inflexible. In terms of the computational cost, training the quadratic models is expensive, while training tnl4 incurs the lowest computational cost. We also note that the performance of the quadratic model, tnl5, is comparable to the best models in our previous paper15 and shown in the supplementary material in terms of RMSE and significantly better in terms of the maximum absolute error in the predicted total kinetic energy. These results indicate the effectiveness of including higher order terms in the gradient expansion of kinetic energy into the nonlocal models. However, this also suggests that by using a relatively small dataset (which may or may not cover the whole phase space of electron and kinetic densities), it is possible to generate models that violate traditional constraints but have lower RMSE values than models that satisfy those constraints.

TABLE I.

A summary of kinetic energy densities and total kinetic energies obtained from nonlocal models and model KSDFT calculations for one-dimensional structures in test datasets for the three nonlocal models proposed in Sec. II. For each system, three structures are used for training and 27 structures are used for testing. Here, R2 and root mean square error (RMSE) are the coefficient of determination and the root mean square error of kinetic energy density (Ha/bohr), respectively, TKSDFT and TModel are the total kinetic energies (Ha) obtained from model one-dimensional KSDFT calculations and the various models, respectively, and Terror is the error in total kinetic energy (Ha) between the one-dimensional model KSDFT value and the model. For each model type, TKE, Tmodel, and Terror only for the structure with the maximum absolute error (i.e., the worst case scenario) are shown, meaning that for different models, the worst case scenarios can correspond to different structures.

SystemsParameterstnl4tnl5 (Quadratic)tnl6 (Liu–Parr)tnl5,scalingx (Coordinate scaling)
H8 R2 0.9941 0.9996 0.9898 0.9953 
RMSE 0.0029 0.0015 0.0038 0.0025 
TKSDFT (Ha) 1.1147 1.1145 1.1145 1.1147 
Tmodel (Ha) 1.1136 1.1152 1.1202 1.1651 
Terror (Ha) 0.0011 0.0007 0.0057 0.0504 
LiH R2 0.9984 0.9987 0.9982 0.9981 
RMSE 0.0038 0.0035 0.0040 0.0041 
TKSDFT (Ha) 0.7583 0.7428 0.8068 0.7338 
Tmodel (Ha) 0.7651 0.7642 0.8566 0.7592 
Terror (Ha) 0.0075 0.0061 0.0498 0.0253 
LiF R2 0.9999 0.9999 0.9994 0.9998 
RMSE 0.0043 0.0031 0.0087 0.0065 
TKSDFT (Ha) 6.7109 7.0165 6.7263 6.9105 
Tmodel (Ha) 6.7355 7.0304 6.7752 7.0209 
Terror (Ha) 0.0246 0.0138 0.0489 0.1103 
C4H2 R2 0.9996 0.9996 0.9986 0.9973 
RMSE 0.0069 0.0075 0.0100 0.0185 
TKSDFT (Ha) 13.1728 12.4828 13.0915 13.4613 
Tmodel (Ha) 13.2473 12.4428 13.1402 13.5490 
Terror (Ha) 0.0745 0.0399 0.0995 0.0878 
C4N2 R2 0.9999 0.9999 0.9993 0.9997 
RMSE 0.0031 0.0039 0.0089 0.0065 
TKSDFT (Ha) 20.7398 20.2343 20.7399 20.6913 
Tmodel (Ha) 20.6410 20.2450 20.8654 20.7679 
Terror (Ha) 0.0987 0.0107 0.1834 0.0765 
C3O2 R2 0.9999 0.9998 0.9995 0.9998 
RMSE 0.0030 0.0039 0.0066 0.0061 
TKSDFT (Ha) 19.3773 19.3775 19.2867 19.3338 
Tmodel (Ha) 19.3316 19.4445 19.1988 19.4088 
Terror (Ha) 0.0456 0.0670 0.0878 0.0750 
SystemsParameterstnl4tnl5 (Quadratic)tnl6 (Liu–Parr)tnl5,scalingx (Coordinate scaling)
H8 R2 0.9941 0.9996 0.9898 0.9953 
RMSE 0.0029 0.0015 0.0038 0.0025 
TKSDFT (Ha) 1.1147 1.1145 1.1145 1.1147 
Tmodel (Ha) 1.1136 1.1152 1.1202 1.1651 
Terror (Ha) 0.0011 0.0007 0.0057 0.0504 
LiH R2 0.9984 0.9987 0.9982 0.9981 
RMSE 0.0038 0.0035 0.0040 0.0041 
TKSDFT (Ha) 0.7583 0.7428 0.8068 0.7338 
Tmodel (Ha) 0.7651 0.7642 0.8566 0.7592 
Terror (Ha) 0.0075 0.0061 0.0498 0.0253 
LiF R2 0.9999 0.9999 0.9994 0.9998 
RMSE 0.0043 0.0031 0.0087 0.0065 
TKSDFT (Ha) 6.7109 7.0165 6.7263 6.9105 
Tmodel (Ha) 6.7355 7.0304 6.7752 7.0209 
Terror (Ha) 0.0246 0.0138 0.0489 0.1103 
C4H2 R2 0.9996 0.9996 0.9986 0.9973 
RMSE 0.0069 0.0075 0.0100 0.0185 
TKSDFT (Ha) 13.1728 12.4828 13.0915 13.4613 
Tmodel (Ha) 13.2473 12.4428 13.1402 13.5490 
Terror (Ha) 0.0745 0.0399 0.0995 0.0878 
C4N2 R2 0.9999 0.9999 0.9993 0.9997 
RMSE 0.0031 0.0039 0.0089 0.0065 
TKSDFT (Ha) 20.7398 20.2343 20.7399 20.6913 
Tmodel (Ha) 20.6410 20.2450 20.8654 20.7679 
Terror (Ha) 0.0987 0.0107 0.1834 0.0765 
C3O2 R2 0.9999 0.9998 0.9995 0.9998 
RMSE 0.0030 0.0039 0.0066 0.0061 
TKSDFT (Ha) 19.3773 19.3775 19.2867 19.3338 
Tmodel (Ha) 19.3316 19.4445 19.1988 19.4088 
Terror (Ha) 0.0456 0.0670 0.0878 0.0750 

Now, we compare the performance of our models with that of the other data-driven models proposed in the literature. Recently, Golub and Manzhos43,45 and Seino et al.56 proposed models to learn the kinetic energy functional of three-dimensional (3D) systems with pseudopotentials using deep neural networks and Gaussian process regression. The performance of their models in terms of RMSE scores and errors in kinetic energy is similar to that of the models proposed in this paper. However, we have only focused on all electron one-dimensional systems. One important difference is that our model has much fewer parameters to fit compared to the neural network based models. Hence, our model requires significantly lesser data points to train. Snyder et al. also proposed a kernel ridge regression based model to fit the total kinetic energy (instead of kinetic energy density) for model one-dimensional systems. The performance of our model in terms of total error in kinetic energy is similar. However, their model requires a much larger dataset to train as they have only a scalar quantity (i.e., total kinetic energy) to fit for a given structure.

Next, to assess the transferability of models outside the training domain, we analyze the performance of the three new models (proposed in Sec. II) in charged systems obtained by inserting or removing electrons. Figure 3 shows the kinetic energy densities of two prototypical systems, C4N2 and H8, obtained from KSDFT calculations using the exact-exchange functional and from the three models (proposed in Sec. II) trained on three neutral structures of C4N2 and H8. Electron densities obtained from KSDFT calculations for these charged systems are used as model inputs. We observe that the quadratic model, tnl5, performs better than the two other models; the insets in Fig. 3 show that kinetic energy densities predicted from the two other models show significant oscillations near local minima and inflection points (of the kinetic energy density). While nonlocal correlations between electron density and its derivatives are present in all three models, the oscillations are absent only in the results obtained from the quadratic model, tnl5. Therefore, we conjecture that this satisfactory performance of quadratic models is related to the cancellation of errors caused by spatial correlations (Fig. 5) between the different derivative terms.

FIG. 3.

Kinetic energy densities obtained from model KSDFT calculations using the exact-exchange functional are compared with model predictions for negatively charged one-dimensional systems, i.e., H82 (a) and C4N22 (b). Electron densities obtained from model one-dimensional KSDFT calculations (using exact-exchange functional) for charged systems are taken as inputs for model predictions. The insets highlight the behavior of the models close to the local minima of the kinetic energy density profiles.

FIG. 3.

Kinetic energy densities obtained from model KSDFT calculations using the exact-exchange functional are compared with model predictions for negatively charged one-dimensional systems, i.e., H82 (a) and C4N22 (b). Electron densities obtained from model one-dimensional KSDFT calculations (using exact-exchange functional) for charged systems are taken as inputs for model predictions. The insets highlight the behavior of the models close to the local minima of the kinetic energy density profiles.

Close modal

To further understand the behavior of these models, variations in total kinetic energies with respect to the total charge on these systems are shown in Fig. 4. For the one-dimensional chain of hydrogens, the total kinetic energies lie in the range of 0.4–1 Ha, which is much smaller than the magnitudes of total kinetic energies of C4N2. The errors in predicted total kinetic energies in the case of systems with ±1 charge are in the range of 10%–15% in the case of H8, but are about 5% in the case of C4N2, suggesting that all the models perform satisfactorily (in terms of errors in kinetic energy) outside the training domain. However, if we consider both the slope of kinetic energy vs net charge on the system and kinetic energy density profiles, the quadratic model performs better than the other two models. These results suggest that a generic model for kinetic energy with good transferability (far away from the training domain) can be designed by using only a few structures for training.

FIG. 4.

Predicted total kinetic energies for charged H8 (a) and C4N2 (b) systems from tnl4, tnl5 (quadratic model), and tnl6 (based on Liu–Parr formulation) models are compared to assess their performance outside the training domain.

FIG. 4.

Predicted total kinetic energies for charged H8 (a) and C4N2 (b) systems from tnl4, tnl5 (quadratic model), and tnl6 (based on Liu–Parr formulation) models are compared to assess their performance outside the training domain.

Close modal

Inspired by the performance of the proposed models for charged systems, we train a universal model by using data from all the one-dimensional systems considered in this study. We use the quadratic model, tnl5, for this purpose because of its superior performance for test systems both inside and outside the training domain. We train a model (universal model 1) by using three structures from each of the six systems, i.e., 18 structures in total. The RMSE and R2 values for the predicted kinetic energy densities of structures not used in training (i.e., 162 structures) for this universal model 1 are 0.0035 and 0.999 Ha/bohr, respectively, and the results are summarized in Table II. The good predictive capability of the quadratic model, tnl5, is clear from the high R2 (more than 0.99) except for H8 and LiH.

TABLE II.

The performance of the three universal quadratic models (described in tnl5) is assessed using different metrics, root mean square errors and R2 scores corresponding to the predicted kinetic energy densities for each system and the maximum absolute errors in predicted total kinetic energies. The universal model 2 corresponds to a model trained with three systems, C4N2, C4H2, and C3O2, and the universal model 3 corresponds to a model trained with three systems, LiH, LiF, and H8. The training set contains three structures from each of the six one-dimensional model systems, and the remaining structures are set aside for testing purposes.

Systems/metricsC4N2C4H2C3O2H8LiHLiF
Universal model 1 R2 0.9998 0.9999 0.9997 0.9068 0.9997 0.9933 
RMSE (Ha/bohr) 0.0060 0.0046 0.0068 0.0113 0.0088 0.0080 
TKSDFT (Ha) 20.4839 19.1199 13.1749 0.4397 6.7529 0.8068 
Tmodel (Ha) 20.4454 19.1496 13.2281 0.2873 6.9459 0.8486 
Terror (Ha) 0.0385 0.0297 0.0532 0.1524 0.1934 0.0418 
Universal model 2 R2 0.9999 0.9999 0.9998 ⋯ ⋯ ⋯ 
RMSE (Ha/bohr) 0.0028 0.0022 0.0034 ⋯ ⋯ ⋯ 
TKSDFT (Ha) 20.4745 19.1213 13.1876 ⋯ ⋯ ⋯ 
Tmodel (Ha) 20.4842 19.1289 13.1703 ⋯ ⋯ ⋯ 
Terror (Ha) 0.0098 0.0076 0.0172 ⋯ ⋯ ⋯ 
Universal model 3 R2 ⋯ ⋯ ⋯ 0.9976 0.9999 0.9999 
RMSE (Ha/bohr) ⋯ ⋯ ⋯ 0.0089 0.0056 0.0049 
TKSDFT (Ha) ⋯ ⋯ ⋯ 0.4410 6.8143 0.7986 
Tmodel (Ha) ⋯ ⋯ ⋯ 0.4841 6.8634 0.8162 
Terror (Ha) ⋯ ⋯ ⋯ 0.0431 0.0490 0.0176 
Systems/metricsC4N2C4H2C3O2H8LiHLiF
Universal model 1 R2 0.9998 0.9999 0.9997 0.9068 0.9997 0.9933 
RMSE (Ha/bohr) 0.0060 0.0046 0.0068 0.0113 0.0088 0.0080 
TKSDFT (Ha) 20.4839 19.1199 13.1749 0.4397 6.7529 0.8068 
Tmodel (Ha) 20.4454 19.1496 13.2281 0.2873 6.9459 0.8486 
Terror (Ha) 0.0385 0.0297 0.0532 0.1524 0.1934 0.0418 
Universal model 2 R2 0.9999 0.9999 0.9998 ⋯ ⋯ ⋯ 
RMSE (Ha/bohr) 0.0028 0.0022 0.0034 ⋯ ⋯ ⋯ 
TKSDFT (Ha) 20.4745 19.1213 13.1876 ⋯ ⋯ ⋯ 
Tmodel (Ha) 20.4842 19.1289 13.1703 ⋯ ⋯ ⋯ 
Terror (Ha) 0.0098 0.0076 0.0172 ⋯ ⋯ ⋯ 
Universal model 3 R2 ⋯ ⋯ ⋯ 0.9976 0.9999 0.9999 
RMSE (Ha/bohr) ⋯ ⋯ ⋯ 0.0089 0.0056 0.0049 
TKSDFT (Ha) ⋯ ⋯ ⋯ 0.4410 6.8143 0.7986 
Tmodel (Ha) ⋯ ⋯ ⋯ 0.4841 6.8634 0.8162 
Terror (Ha) ⋯ ⋯ ⋯ 0.0431 0.0490 0.0176 

To further investigate this, the kinetic energy density obtained from model one-dimensional KSDFT calculations using exact-exchange functional and the universal model are shown in Fig. 5. It is clear that the model is able to capture subtle features of the kinetic energy density for C4N2, C4H2, and C3O2. However, deviations from one-dimensional model KSDFT results are clearly evident in H8 and LiH. We conjecture that this discrepancy is related to the magnitude of kinetic energy densities of these systems; since the kinetic energy densities of H8 and LiH are about an order of magnitude smaller than other systems, the universal model 1 is biased toward larger values of kinetic energy densities. The kernels for quadratic models obtained for the individual systems and for the universal model 1 (highlighted in black) are shown in Fig. 6. It is clear that the kernels for LiH and H8 are different from the other systems, and the kernel in the universal model is similar to the systems with larger kinetic energy densities. One possible remedy is to explore different objective functions to obtain the model parameters. Another option is to use feature selection to design a better quadratic model with one or more kernels, and we plan to explore this in the future.

FIG. 5.

Kinetic energy densities (for each of the six one-dimensional systems explored in this study) predicted using the universal quadratic model 1 for structures in the test set: C4N2 (a), C4H2 (b), C3O2 (c), H8 (d), LiF (e), and LiH (f). The maximum absolute errors in predicted total kinetic energies, RMSEs, and high R2 scores for individual systems are also shown. The insets show that model predictions are able to capture subtle features of the KSDFT kinetic energy densities. The training set used to train the universal model contains three structures from each of the six systems, and the remaining structures are set aside for testing.

FIG. 5.

Kinetic energy densities (for each of the six one-dimensional systems explored in this study) predicted using the universal quadratic model 1 for structures in the test set: C4N2 (a), C4H2 (b), C3O2 (c), H8 (d), LiF (e), and LiH (f). The maximum absolute errors in predicted total kinetic energies, RMSEs, and high R2 scores for individual systems are also shown. The insets show that model predictions are able to capture subtle features of the KSDFT kinetic energy densities. The training set used to train the universal model contains three structures from each of the six systems, and the remaining structures are set aside for testing.

Close modal
FIG. 6.

Kernels obtained from quadratic models trained for the each of the six one-dimensional model systems explored in this study and the kernel obtained from the universal quadratic model described in Sec. III D.

FIG. 6.

Kernels obtained from quadratic models trained for the each of the six one-dimensional model systems explored in this study and the kernel obtained from the universal quadratic model described in Sec. III D.

Close modal

Interestingly, we also observe large errors in the total kinetic energy for LiH (0.19 Ha) and H8 (0.15 Ha), suggesting that there are no cancellation of errors. Hence, we train two different quadratic models, universal model 2 with structures from LiH, LiF, and H8 and universal model 3 with structures from C4H2, C4N2, and C3O2. We use three structures from each material, hence a total of nine structures for each of the two models. We observe that errors in kinetic energies and kinetic energy densities (summarized in Table II) for all systems have reduced significantly for universal model 2 and universal model 3 compared to universal model 1. This validates our previous point that universal models work better when similar types of systems are used for training, but this behavior (i.e., failure of universal model 1 to handle a diverse set of systems) suggests that our model is still not complete, and the true kinetic energy density functional perhaps is more complex and incorporates higher order spatial correlations. For comparison, we also train a universal model using the Q1 + Q3 model proposed in Ref. 15. However, for the 162 structures in the test set, the RMSE and R2 scores for the predicted kinetic energy densities are 0.061 and 0.983 Ha/bohr, respectively. Therefore, the quadratic model has better transferability and can be used to develop generic models.

Inspired by the results of the universal model based on the quadratic model, we use the quadratic model to fit a density functional that includes the kinetic energy and the exact-exchange energy. The goal is to obtain a functional which in conjunction with the external potential and Hartree energy functionals can be used to perform an OFDFT calculation with the exact-exchange level of accuracy. Thus, we define a modified kinetic energy by adding the two orbital dependent terms (i.e., kinetic energy and exact-exchange energy densities). However, for consistency, we also subtract a local or semi-local exchange energy, i.e., we define T̄KEρ=Ts[ρ]+EX[ρ({ψk})]ELDA[ρ]. Here, we use the LDA exchange energy as we have already fitted the Kohn–Sham kinetic energy for the LDA in our previous work. Therefore, the expression for the modified kinetic energy functional is given by

T̄KE=Ts+EXELDA,EX=12fifjiNsjNsψi(x)ψj(x)ψi(x)ψj(x)×11+|xx|2dxdx,ELDA[ρ]=vLDA(ρ(x))ρ(x)dx,vLDA=Aπtan1πρ2κ,
(13)

where T̄KE is the modified kinetic energy for which we want to develop a density-dependent model, Ts is the Kohn–Sham kinetic energy, EX is the exact-exchange energy, and ELDA is the local density approximation (LDA) exchange energy. The modified kinetic energy essentially captures the kinetic energy of interacting electrons, hence justifying the use of the kinetic energy density models developed in this paper.

We fit the quadratic model for modified kinetic energy [in Eq. (13)] for the six systems. For training the model, we randomly select three structures and set aside the remaining 27 structures for testing. The parameters of the fitting, such as the number of Gaussian functions and degree of polynomial terms in the local term, are the same as the ones used in training the Kohn–Sham kinetic energy density in Secs. III B and III D. We summarize the performance of the quadratic model for modified kinetic energy in the third column of Table III. The plots showing modified kinetic energy densities obtained from KSDFT calculations and the model predictions are also shown in Fig. 7. We observe that the R2 scores for all systems are greater than 0.99, indicating a very good performance of the quadratic model. From Fig. 7, we observe that the quadratic models are able to capture subtle features present in the modified kinetic energy density profiles. The maximum error in the fit is observed near the minima points in the kinetic energy density, a pattern that has also been observed in the models developed by Golub and Manzhos.43 One of the reasons for this behavior may be the lack of data points in these regions, as pointed out in Ref. 43.

TABLE III.

A summary of modified kinetic energy densities (T̄KE=Ts+EXELDA) and total kinetic energies obtained from the quadratic model and one-dimensional model KSDFT calculations for structures in test datasets for the three nonlocal models proposed in Sec. II. The individual model is trained separately for each system with three structures. Universal model 1 is a single global model trained using a collection of 18 structures with three structures taken from each system, and universal model 2 is trained with H8, LiH, and LiF with three structures from each system, and universal model 3 is trained with C4H2, C4N2, and C3O2 with three structures from each system. Here, R2 and RMSE are the coefficients of determination and the root mean square error of kinetic energy density (Ha/bohr), respectively, TKSDFT and TModel are the total kinetic energies (Ha) obtained from one-dimensional model KSDFT calculations and the various models, respectively, and Terror is the error in total kinetic energy (Ha) between the 1D model KSDFT value and the model. For each model type, TKSDFT, Tmodel, and Terror only for the structure with the maximum absolute error are shown.

SystemsParametersIndividual modelUniversal model 1Universal model 2SystemsParametersIndividual modelUniversal model 1Universal model 3
H8 R2 0.9914 0.8817 0.9876 C4H2 R2 0.9989 0.9984 0.9992 
RMSE 0.0039 0.0140 0.0054 RMSE 0.0074 0.0086 0.0081 
T̄KSDFT (Ha) −1.7549 −1.1145 −1.243 T̄KSDFT (Ha) 3.1179 3.2522 3.1183 
T̄model (Ha) −1.6669 −1.8698 −1.3449 T̄model (Ha) 3.0108 3.3805 3.2315 
T̄error (Ha) 0.0879 0.2696 0.1029 T̄error (Ha) 0.1071 0.1282 0.1132 
LiH R2 0.9969 0.9711 0.9903 C4N2 R2 0.9997 0.9992 0.9993 
RMSE 0.0030 0.0093 0.0065 RMSE 0.0044 0.0082 0.0051 
T̄KSDFT (Ha) −0.6192 −0.6193 −0.6192 T̄KSDFT (Ha) 7.1003 6.8673 6.8673 
T̄model (Ha) −0.6063 −0.5541 −0.5905 T̄model (Ha) 7.1551 6.7514 6.9454 
T̄error (Ha) 0.0129 0.0652 0.0287 T̄error (Ha) 0.0548 0.1158 0.0781 
LiF R2 0.9999 0.9988 0.9994 C3O2 R2 0.9996 0.9995 0.9992 
RMSE 0.0035 0.0112 0.0045 RMSE 0.0064 0.0069 0.0073 
T̄KSDFT (Ha) 2.4865 2.4040 2.4576 T̄KSDFT (Ha) 6.7275 7.1814 6.8234 
T̄model (Ha) 2.5129 2.6726 2.4963 T̄model (Ha) 6.6423 7.0795 6.8357 
T̄error (Ha) 0.0265 0.2686 0.0387 T̄error (Ha) 0.0851 0.1018 0.0123 
SystemsParametersIndividual modelUniversal model 1Universal model 2SystemsParametersIndividual modelUniversal model 1Universal model 3
H8 R2 0.9914 0.8817 0.9876 C4H2 R2 0.9989 0.9984 0.9992 
RMSE 0.0039 0.0140 0.0054 RMSE 0.0074 0.0086 0.0081 
T̄KSDFT (Ha) −1.7549 −1.1145 −1.243 T̄KSDFT (Ha) 3.1179 3.2522 3.1183 
T̄model (Ha) −1.6669 −1.8698 −1.3449 T̄model (Ha) 3.0108 3.3805 3.2315 
T̄error (Ha) 0.0879 0.2696 0.1029 T̄error (Ha) 0.1071 0.1282 0.1132 
LiH R2 0.9969 0.9711 0.9903 C4N2 R2 0.9997 0.9992 0.9993 
RMSE 0.0030 0.0093 0.0065 RMSE 0.0044 0.0082 0.0051 
T̄KSDFT (Ha) −0.6192 −0.6193 −0.6192 T̄KSDFT (Ha) 7.1003 6.8673 6.8673 
T̄model (Ha) −0.6063 −0.5541 −0.5905 T̄model (Ha) 7.1551 6.7514 6.9454 
T̄error (Ha) 0.0129 0.0652 0.0287 T̄error (Ha) 0.0548 0.1158 0.0781 
LiF R2 0.9999 0.9988 0.9994 C3O2 R2 0.9996 0.9995 0.9992 
RMSE 0.0035 0.0112 0.0045 RMSE 0.0064 0.0069 0.0073 
T̄KSDFT (Ha) 2.4865 2.4040 2.4576 T̄KSDFT (Ha) 6.7275 7.1814 6.8234 
T̄model (Ha) 2.5129 2.6726 2.4963 T̄model (Ha) 6.6423 7.0795 6.8357 
T̄error (Ha) 0.0265 0.2686 0.0387 T̄error (Ha) 0.0851 0.1018 0.0123 
FIG. 7.

Shown here are the modified kinetic energy densities (for each of the six one-dimensional systems explored in this study) predicted using the quadratic model for structures in the test set: C4N2 (a), C4H2 (b), C3O2 (c), H8 (d), LiF (e), and LiH (f). The insets show that model predictions are able to capture subtle features of the model one-dimensional KSDFT kinetic energy densities. The training set used to train the universal model contains three structures from a given system, and the remaining structures are set aside for testing.

FIG. 7.

Shown here are the modified kinetic energy densities (for each of the six one-dimensional systems explored in this study) predicted using the quadratic model for structures in the test set: C4N2 (a), C4H2 (b), C3O2 (c), H8 (d), LiF (e), and LiH (f). The insets show that model predictions are able to capture subtle features of the model one-dimensional KSDFT kinetic energy densities. The training set used to train the universal model contains three structures from a given system, and the remaining structures are set aside for testing.

Close modal

To further check the validity of the quadratic model for modified kinetic energy density, we train three universal models with data from all the systems using the modified kinetic energy density in Eq. (13). We use three structures from each system for training the model. The results summarized in the fourth column of Table III show that R2 scores for universal model 1 are similar to the models that were trained for the individual systems except for LiH and H8. The discrepancy for these two materials is related to the significantly smaller kinetic energies of these two systems compared to the others as we noted in the case of universal model for the Kohn–Sham kinetic energy. However, for the universal models 2 and 3, the model predictions agree very well with KSDFT results, and the R2 scores are very high. These results suggest that quadratic models can be used to describe kinetic energy and exact-exchange energy with good accuracy.

In this study, using data from a few model one-dimensional systems, such as H8, C4H2, C4N2, LiH, LiF, and C3O2, we explored a variety of empirical models of the kinetic energy functional. To capture the subtle features of the kinetic energy density obtained from one-dimensional model DFT calculations with the exact-exchange functional, we used a variety of models proposed in our previous work and three nonlocal models developed in this work. These nonlocal models are inspired by physics-based models that have been proposed in the literature and incorporate spatial correlations between electron density and its higher order derivatives. The nonlocal kernels used to incorporate spatial correlations in the electron density and its derivatives are modeled by using a set of Gaussian basis functions. The predictive capabilities of the proposed models are systematically tested for structures both similar to those present in training sets and for charged systems (not present in training sets). The models are able to predict kinetic energy densities to a very high level of accuracy, i.e., with R2 scores >0.99, and total kinetic energies with maximum errors of less than 0.001 Ha/atom and a mean absolute error of about 0.004 Ha/atom. We do observe large errors (maximum absolute error of 0.19 Ha) in total kinetic energy for universal models trained with all the systems; however, those errors are greatly reduced (maximum absolute error of 0.04 Ha and mean absolute errors are less than 0.01 Ha) if the universal models are trained on structures with similar magnitudes of kinetic energy density values. In addition to kinetic energy density, we show that the models are able to predict the exact-exchange energy density to a high level of accuracy as well.

If we compare the R2 scores, RMSE values for predicted kinetic energy densities, and the maximum absolute errors in the total kinetic energies, the quadratic model performs much better than the two other nonlocal models proposed in this study. This implies that incorporating spatial correlations between the electron density and its higher order derivatives (up to the fourth order) helps in capturing salient features of kinetic energies obtained from one-dimensional model KSDFT with the exact-exchange functional. This implies that it is, indeed, possible to generate reliable models within the framework of a data-informed physics-based learning paradigm by using only a few structures. Furthermore, we show that it is possible to generate generic kinetic energy models for a broad class of heterogeneous systems. Even in this case, the models have good predictive capability, and the quadratic model performs better than the other models analyzed in this study.

Although the models are validated for model one-dimensional systems, they can be generalized to three-dimensional systems. We plan to tackle this as well as the issues related to learning the nonlocal part of pseudopotential in a future publication. There are still a few challenges associated with these models, namely, effective optimization scheme to obtain the best set of parameters for the quadratic model, data selection, and avoiding bias toward systems with large kinetic energies, and we plan to tackle them in the future. Overall, the models considered in this study are a step forward in the direction of generating accurate, reliable, and computationally efficient data-driven orbital-free density functional theory models.

See the supplementary material for the table of results for models proposed in Ref. 15 and the Matlab codes to fit the models in this paper.

This work was performed with LDRD funding (Grant No. 21-ERD-005) under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Computing support for this work came from the Lawrence Livermore National Laboratory (LLNL) institutional computing facility. The authors acknowledge helpful discussions with Xin Jing about the implementation of the exact-exchange functional.

The authors have no conflicts of interest to disclose..

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

In this section, we briefly outline formulation and implementation of a one-dimensional model analog of Kohn–Sham DFT1,2 with the exact-exchange potential with open boundary conditions. In what follows, x denotes the coordinate in 1D space. The ground state can be determined by solving the following nonlinear eigenvalue problem self-consistently:

122x2+vextx+vHρx+VX{ψj}(x)ψix=εiψix,ψixψjxdx=δij,ρx=j=1Nsfjψj2x,fj{1,2},
(A1)

where ψi and ɛi are the model Kohn–Sham eigenfunctions and eigenvalues, respectively, vext(x) is the external potential, vH[ρ](x) is the Hartree potential, VX{ψj}(x,x) is the exact-exchange operator, Ns is the total number of occupied orbitals, and fi is the occupation number of the ith orbital. In our model system, we replace the Coulomb form of Hartree potential and external nuclear potential by a soft form,57 

vHρ(x)=Rρxdx1+xx2,vextx=IZI1+xXI2.
(A2)

Here, ZI and XI denote the nuclear charge and position of the Ith atom, respectively. The exact-exchange operator VX is a nonlocal operator that depends on the occupied orbitals {ψj}j=1Ns and has the following form:58,59

VX{ψj}(x,x)=j=1Nsfjψj(x)ψj(x)1+xx2.
(A3)

The action of VX on a given orbital ψx is given by58,59

VX{ψj}ψ(x)=j=1Nsfjψjxψxψjx1+xx2dx.
(A4)

Therefore, the Hamiltonian in the eigenvalue problem outlined in Eq. (A1) depends both on the electron density, ρx, and the orbitals, {ψj}j=1Ns. Hence, we solve this eigenvalue problem by using two sets of self-consistent field iterations;58 in the outer loop, the exchange operator is updated, and in the inner loop, the electron density is updated using Pulay mixing60–62 while keeping the exchange operator constant.

We solve the partial differential equation in Eq. (A1) numerically by using a finite-difference discretization method in which a supercell with edge length L is discretized using an uniform finite-difference grid of spacing h. Grid points are indexed from 1 to N, where N = L/h. The second derivative in Eq. (A1) is approximated using the central finite-difference approximation,63–65 

2gx2(i)=k=0NFDWkg(i+k)+g(ik)+O(h2NFD),
(B1)

where g(i) denotes the value of the function at the ith grid point and NFD is the order of finite-difference approximation. The weights, Wk, associated with the function values are

W0=1h2k=1NFD1k2,Wk=2(1)k+1h2k2(NFD!)2(NFDk)!(NFD+k)!,k=1,2,,NFD.
(B2)

With the above discretization scheme, the eigenvalue equation in Eq. (A1) can be written in terms of matrix eigenvalue problem,

Hψ=εψ,HRN×N,ψRN×1,Hψ=12h2ψ+vH[ρ]+vextψ+VX{ψj}ψ,ρ(i)=j=1Nsfjψj(i)ψj(i),
(B3)

where h2 is the finite-difference discrete second derivative matrix and ψji is the value of the jth wavefunction at the ith grid point.

The eigenvalue problem outlined above in Eq. (B3) is solved in Matlab by using the inbuilt eigensolver. With a given initial guess of electron density and the exchange operator, these equations are solved by using two sets of nested self-consistent iterations until convergence to a desired accuracy is reached. The convergence for both the inner and outer loops is checked with respect to relative error in electron density, and the change in the energy of the system and the self-consistent iterations are considered to have converged when the change in energy per atom is less than 10−6 Ha. The kinetic energy and kinetic energy density are obtained by using the converged orbitals using the following expressions:

Ts=i=1Nj=1Nshψj(i)h2ψj(i),ts(i)=j=1Nsψj(i)h2ψj(i).
(B4)

Here, Ts = ∫ts(x)dx and ts(x) are the total kinetic energy and kinetic energy density, respectively. Information about the mesh sizes used to discretize the domain is given in Table IV below.

TABLE IV.

Parameters used for solving the one-dimensional model Kohn–Sham equation using the exact-exchange functional.

SystemsH8, LiH, LiFC4N2C4H2, C3O2
Grid size (Bohr) 0.05 0.06 0.06 
Edge length (Bohrs) 50 100 100 
SystemsH8, LiH, LiFC4N2C4H2, C3O2
Grid size (Bohr) 0.05 0.06 0.06 
Edge length (Bohrs) 50 100 100 

We obtain the analytical expressions of the first and second order derivative terms of P(r, r′) appearing in Eq. (12) for the nonlocal model based on the density matrix proposed by Liu and Parr. These expressions can be substituted in Eq. (12) to obtain the expression of kinetic energy density tnl(6)(r) as a function of the fitting coefficients, which are later obtained by minimizing the error from the kinetic energy density obtained from one-dimensional model KSDFT simulations,

P(r,r)=k=1NGdkeβk|rr|2ρ(r)+ρ(r)2,rP(r,r)=k=1NGdkeβkα2βk|rr|2[ρ(r)+ρ(r)]rρ(r)2βk(rr)[ρ(r)+ρ(r)]2,rP(r,r)=k=1NGdkeβkα2βk|rr|2[ρ(r)+ρ(r)]rρ(r)2βk(rr)[ρ(r)+ρ(r)]2,rrP(r,r)=k=1NGdkeβkα[4βk(rr)[ρ(r)+ρ(r)]]rρ(r)rρ(r)2βk|rr|2rρ(r)rρ(r)+2βk[ρ(r)+ρ(r)]2+k=1NGdkeβkα2βk|rr|2[ρ(r)+ρ(r)]rρ(r)2βk(rr)[ρ(r)+ρ(r)]22βk|rr|2[ρ(r)+ρ(r)]rρ(r)+2βk(rr)[ρ(r)+ρ(r)]2,α=|rr|2ρ(r)+ρ(r)2.
(C1)

Here, dk are the coefficients and βk and NG are the widths and the number of the Gaussian functions used to model the function P(r, r′).

1.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
2.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
3.
F.
Brockherde
,
L.
Vogt
,
L.
Li
,
M. E.
Tuckerman
,
K.
Burke
, and
K.-R.
Müller
,
Nat. Commun.
8
,
872
(
2017
).
4.
M.
Brack
,
Rev. Mod. Phys.
65
,
677
(
1993
).
5.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
47
,
558
(
1993
).
6.
G.
Kresse
and
J.
Hafner
,
Phys. Rev. B
48
,
13115
(
1993
).
7.
S.
Ghosh
and
P.
Suryanarayana
,
Comput. Phys. Commun.
212
,
189
(
2017
).
8.
Q.
Xu
,
A.
Sharma
,
B.
Comer
,
H.
Huang
,
E.
Chow
,
A. J.
Medford
,
J. E.
Pask
, and
P.
Suryanarayana
,
SoftwareX
15
,
100709
(
2021
).
9.
M. M. G.
Alemany
,
X.
Huang
,
M. L.
Tiago
,
L. J.
Gallego
, and
J. R.
Chelikowsky
,
Solid State Commun.
146
,
245
(
2008
), special section: Negative Refraction and Metamaterials for Optical Science and Engineering.
10.
W.
Mi
and
M.
Pavanello
,
Phys. Rev. B
100
,
041105
(
2019
).
11.
X.
Shao
,
K.
Jiang
,
W.
Mi
,
A.
Genova
, and
M.
Pavanello
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
11
,
e1482
(
2021
).
12.
S.
Carr
,
S.
Fang
, and
E.
Kaxiras
,
Nat. Rev. Mater.
5
,
748
(
2020
).
13.
H.-S.
Do
and
B.-J.
Lee
,
Sci. Rep.
8
,
4851
(
2018
).
14.
A. C.
Pan
,
D.
Jacobson
,
K.
Yatsenko
,
D.
Sritharan
,
T. M.
Weinreich
, and
D. E.
Shaw
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
4244
(
2019
).
15.
S.
Kumar
,
E.
Landinez
,
B.
Sadigh
,
S.
Hamel
,
B.
Gallagher
,
V.
Bulatov
,
J.
Klepeis
, and
A.
Samanta
, “
Accurate parametrization of the kinetic energy functional
,”
J. Chem Phys.
(to be published).
16.
W. C.
Witt
,
B. G.
del Rio
,
J. M.
Dieterich
, and
E. A.
Carter
,
J. Mater. Res.
33
,
777
(
2018
).
17.
V.
Bach
and
L.
Delle Site
,
Many-Electron Approaches in Physics, Chemistry and Mathematics
(
Springer
,
2014
).
18.
Y. A.
Wang
and
E. A.
Carter
, “
Orbital-free kinetic energy density functional theory
,” in
TheoreticalMethods in Condensed Phase Chemistry
, Progress in Theoretical Chemistry and Physics, edited by
S. D.
Schwartz
(
Kluwer, Dordrecht
,
2000
), pp.
117
184
.
19.
E. A.
Carter
, “
Challenges in modeling materials properties without experimentalInput
,”
Science
321
(
5890
),
800
803
(
2008
).
20.
L.
Šamaj
and
J.
Percus
,
J. Chem. Phys.
111
,
1809
(
1999
).
21.
D. R.
Murphy
,
Phys. Rev. A
24
,
1682
(
1981
).
22.
W.-P.
Wang
,
R. G.
Parr
,
D. R.
Murphy
, and
G. A.
Henderson
,
Chem. Phys. Lett.
43
,
409
(
1976
).
23.
C. Q.
Ma
and
V.
Sahni
,
Phys. Rev. B
16
,
4249
(
1977
).
24.
C.
Lee
and
S. K.
Ghosh
,
Phys. Rev. A
33
,
3506
(
1986
).
25.
R. G.
Parr
,
Horizons of Quantum Chemistry
(
Springer
,
1980
), pp.
5
15
.
26.
E.
Chacón
,
J. E.
Alvarellos
, and
P.
Tarazona
,
Phys. Rev. B
32
,
7868
(
1985
).
27.
L.-W.
Wang
and
M. P.
Teter
,
Phys. Rev. B
45
,
13196
(
1992
).
28.
W. C.
Witt
and
E. A.
Carter
,
Phys. Rev. B
100
,
125107
(
2019
).
29.
V. L.
Lignères
and
E. A.
Carter
,
Handbook of Materials Modeling
(
Springer
,
2005
), pp.
137
148
.
30.
Y. A.
Wang
,
N.
Govind
, and
E. A.
Carter
,
Phys. Rev. B
60
,
16350
(
1999
).
31.
D.
Garcia-Aldea
and
J. E.
Alvarellos
,
Phys. Rev. A
76
,
052504
(
2007
).
32.
D.
Garcia-Aldea
and
J. E.
Alvarellos
,
J. Chem. Phys.
127
,
144109
(
2007
).
33.
W.
Mi
,
A.
Genova
, and
M.
Pavanello
,
J. Chem. Phys.
148
,
184107
(
2018
).
34.
Q.
Xu
,
J.
Lv
,
Y.
Wang
, and
Y.
Ma
,
Phys. Rev. B
101
,
045110
(
2020
).
35.
D.
García-Aldea
and
J. E.
Alvarellos
,
Phys. Chem. Chem. Phys.
14
,
1756
(
2012
).
36.
C.
Huang
and
E. A.
Carter
,
Phys. Rev. B
81
,
045206
(
2010
).
37.
J. P.
Perdew
and
L. A.
Constantin
,
Phys. Rev. B
75
,
155109
(
2007
).
38.
S.
Laricchia
,
L. A.
Constantin
,
E.
Fabiano
, and
F.
Della Sala
,
J. Chem. Theory Comput.
10
,
164
(
2014
).
39.
S.
Liu
,
D.
Zhao
,
C.
Rong
,
T.
Lu
, and
S.
Liu
,
J. Chem. Phys.
150
,
204106
(
2019
).
40.
L. A.
Constantin
,
E.
Fabiano
, and
F.
Della Sala
,
J. Chem. Theory Comput.
15
,
3044
(
2019
).
41.
L. A.
Constantin
,
E.
Fabiano
, and
F.
Della Sala
,
J. Phys. Chem. Lett.
9
,
4385
(
2018
).
42.
A. A.
Kugler
,
Phys. Rev. A
41
,
3489
(
1990
).
43.
P.
Golub
and
S.
Manzhos
,
Phys. Chem. Chem. Phys.
21
,
378
(
2019
).
44.
J. C.
Snyder
,
M.
Rupp
,
K.
Hansen
,
K.-R.
Müller
, and
K.
Burke
,
Phys. Rev. Lett.
108
,
253002
(
2012
).
45.
S.
Manzhos
and
P.
Golub
,
J. Chem. Phys.
153
,
074104
(
2020
).
46.
S.
Manzhos
,
Mach. Learn.: Sci. Technol.
1
,
013002
(
2020
).
47.
M. M.
Denner
,
M. H.
Fischer
, and
T.
Neupert
,
Phys. Rev. Res.
2
,
033388
(
2020
).
48.
F.
Imoto
,
M.
Imada
, and
A.
Oshiyama
,
Phys. Rev. Res.
3
,
033198
(
2021
).
49.
S. A.
Ghasemi
and
T. D.
Kühne
,
J. Chem. Phys.
154
,
074107
(
2021
).
50.
K.
Yao
and
J.
Parkhill
,
J. Chem. Theory Comput.
12
,
1139
(
2016
).
51.
E.
Fabiano
and
L. A.
Constantin
,
Phys. Rev. A
87
,
012511
(
2013
).
52.
A.
Borgoo
and
D. J.
Tozer
,
J. Chem. Theory Comput.
9
,
2250
(
2013
).
53.
A.
Borgoo
,
A. M.
Teale
, and
D. J.
Tozer
,
Phys. Chem. Chem. Phys.
16
,
14578
(
2014
).
54.
M.
Levy
and
J. P.
Perdew
,
Phys. Rev. A
32
,
2010
(
1985
).
55.
S.
Liu
and
R. G.
Parr
,
Chem. Phys. Lett.
278
,
341
(
1997
).
56.
J.
Seino
,
R.
Kageyama
,
M.
Fujinami
,
Y.
Ikabata
, and
H.
Nakai
,
J. Chem. Phys.
148
,
241705
(
2018
).
57.
L. O.
Wagner
,
E. M.
Stoudenmire
,
K.
Burke
, and
S. R.
White
,
Phys. Chem. Chem. Phys.
14
,
8581
(
2012
).
58.
L.
Lin
,
J. Chem. Theory Comput.
12
,
2242
(
2016
).
59.
X.
Wu
,
A.
Selloni
, and
R.
Car
,
Phys. Rev. B
79
,
085102
(
2009
).
60.
P.
Pulay
,
Chem. Phys. Lett.
73
,
393
(
1980
).
61.
P. P.
Pratapa
and
P.
Suryanarayana
,
Chem. Phys. Lett.
635
,
69
(
2015
).
62.
A. S.
Banerjee
,
P.
Suryanarayana
, and
J. E.
Pask
,
Chem. Phys. Lett.
647
,
31
(
2016
).
63.
P.
Suryanarayana
and
D.
Phanish
,
J. Comput. Phys.
275
,
524
(
2014
).
64.
B.
Fornberg
,
Math. Comput.
51
,
699
(
1988
).
65.
S.
Ghosh
and
P.
Suryanarayana
,
J. Comput. Phys.
307
,
634
(
2016
).

Supplementary Material