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, C_{4}H_{2}, C_{4}N_{2}, and C_{3}O_{2}. 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.

## I. INTRODUCTION

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 *Q*_{1} + *Q*_{3} 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, H_{8} and C_{4}N_{2}, 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.

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 (*Q*_{1} + *Q*_{3}) proposed in our previous work (Ref. 15) is able to capture the important features of kinetic energy density, but the R^{2} 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.

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, C_{4}H_{2}, C_{4}N_{2}, and C_{3}O_{2}. 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.

## II. KINETIC ENERGY FUNCTIONAL MODELS

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,

Here, *c*_{j} and *d*_{i,j} are the coefficients of the model, *β*_{i,j} are the widths of the Gaussian functions, and *N*_{G} is the number of Gaussian basis functions used to model the nonlocal kernel. The three terms present in the local kinetic energy expression, *t*_{local}, 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:

where **t**_{KSDFT} is the vector that contains information about kinetic energy densities at different grid points, **c** = {*c*_{1}, *c*_{2}, …, *d*_{1,1}, *d*_{1,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 **t**_{KSDFT}, 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 A–II 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.

### A. A nonlocal model inspired by a model proposed by Chacón *et al.*

Chacón *et al.*^{26} introduced the notion of a locally averaged density, $\rho \u0304r$, to incorporate spatial correlations into the kinetic energy functional,

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

Here, $TvW\rho $ is the von Weizsäcker term, $TLDA\rho $ is the Thomas–Fermi term of the kinetic energy functional, and *t* is the function of density such that $TLDA\rho =\u222b\rho rt\rho (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

where *d*_{p} 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 $\rho \u03042x,\rho \u03043x$, 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, $\rho x$ and $\u2207\rho x$, and consider the following nonlocal model:

and

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 *Q*_{4} by using a set of basis functions,

Interestingly, the model is linear with respect to all the coefficients; however, there is a coupling between the coefficients *d*_{k} 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 *d*_{k} while keeping the parameters *α*_{1}, *α*_{2} constant. In the second step, we fix *d*_{k} and solve another linear system to obtain **c** and *α*_{1}, *α*_{2}. This two-step scheme is repeated until convergence.

### B. A quadratic model

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\rho x,c$,

Here, *Q*_{5} is the nonlocal kernel that is modeled using a set of Gaussian basis functions as discussed in Sec. II A, and $tlocal\rho 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

We note that the functional form of both $tnl5x$ and $tnl5,scalingx$ is quadratic $c\u0304$, and we adopt a two-step optimization scheme to obtain an optimal set of coefficients: In the first step, we initialize $c\u0304$ 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\u0304$ by using a trust-region algorithm. These two steps are repeated until convergence.

### C. A nonlocal model based on a kinetic energy functional proposed by Liu and Parr

Liu and Parr^{55} argued that the Kohn–Sham kinetic energy can be expressed as a functional of degree one in electron density,

where

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

where *d*_{k} are the coefficients and *β*_{k} are the widths of the Gaussian functions used to model the function $Px,x\u2032$.

It is easy to see that after the derivatives of $Px,x\u2032$ 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:

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.

## III. RESULTS AND DISCUSSIONS

### A. DFT data and model training

We use six different one-dimensional systems, namely, H_{8}, LiH, LiF, C_{4}N_{2}, C_{4}H_{2}, and C_{3}O_{2}, 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 *Q*_{1}. Similarly, the model containing the local terms and the third (fourth) term in Eq. (1) is denoted by *Q*_{2} (*Q*_{3}).

### B. Kinetic energy density prediction

The performance of the eight models proposed in Ref. 15 is measured in terms of the coefficient of determination (i.e., R^{2} 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 R^{2} 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 *Q*_{1} + *Q*_{3} nonlocal kernels perform better than the other models for almost all systems. For some systems, *Q*_{1} + *Q*_{2} + *Q*_{3} performs better than *Q*_{1} + *Q*_{3}, but it also contains a significantly larger number of coefficients compared to *Q*_{1} + *Q*_{3} 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 paper^{15} 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.

Systems . | Parameters . | $tnl4$ . | $tnl5$ (Quadratic) . | $tnl6$ (Liu–Parr) . | $tnl5,scalingx$ (Coordinate scaling) . |
---|---|---|---|---|---|

H_{8} | R^{2} | 0.9941 | 0.9996 | 0.9898 | 0.9953 |

RMSE | 0.0029 | 0.0015 | 0.0038 | 0.0025 | |

T_{KSDFT} (Ha) | 1.1147 | 1.1145 | 1.1145 | 1.1147 | |

T_{model} (Ha) | 1.1136 | 1.1152 | 1.1202 | 1.1651 | |

T_{error} (Ha) | 0.0011 | 0.0007 | 0.0057 | 0.0504 | |

LiH | R^{2} | 0.9984 | 0.9987 | 0.9982 | 0.9981 |

RMSE | 0.0038 | 0.0035 | 0.0040 | 0.0041 | |

T_{KSDFT} (Ha) | 0.7583 | 0.7428 | 0.8068 | 0.7338 | |

T_{model} (Ha) | 0.7651 | 0.7642 | 0.8566 | 0.7592 | |

T_{error} (Ha) | 0.0075 | 0.0061 | 0.0498 | 0.0253 | |

LiF | R^{2} | 0.9999 | 0.9999 | 0.9994 | 0.9998 |

RMSE | 0.0043 | 0.0031 | 0.0087 | 0.0065 | |

T_{KSDFT} (Ha) | 6.7109 | 7.0165 | 6.7263 | 6.9105 | |

T_{model} (Ha) | 6.7355 | 7.0304 | 6.7752 | 7.0209 | |

T_{error} (Ha) | 0.0246 | 0.0138 | 0.0489 | 0.1103 | |

C_{4}H_{2} | R^{2} | 0.9996 | 0.9996 | 0.9986 | 0.9973 |

RMSE | 0.0069 | 0.0075 | 0.0100 | 0.0185 | |

T_{KSDFT} (Ha) | 13.1728 | 12.4828 | 13.0915 | 13.4613 | |

T_{model} (Ha) | 13.2473 | 12.4428 | 13.1402 | 13.5490 | |

T_{error} (Ha) | 0.0745 | 0.0399 | 0.0995 | 0.0878 | |

C_{4}N_{2} | R^{2} | 0.9999 | 0.9999 | 0.9993 | 0.9997 |

RMSE | 0.0031 | 0.0039 | 0.0089 | 0.0065 | |

T_{KSDFT} (Ha) | 20.7398 | 20.2343 | 20.7399 | 20.6913 | |

T_{model} (Ha) | 20.6410 | 20.2450 | 20.8654 | 20.7679 | |

T_{error} (Ha) | 0.0987 | 0.0107 | 0.1834 | 0.0765 | |

C_{3}O_{2} | R^{2} | 0.9999 | 0.9998 | 0.9995 | 0.9998 |

RMSE | 0.0030 | 0.0039 | 0.0066 | 0.0061 | |

T_{KSDFT} (Ha) | 19.3773 | 19.3775 | 19.2867 | 19.3338 | |

T_{model} (Ha) | 19.3316 | 19.4445 | 19.1988 | 19.4088 | |

T_{error} (Ha) | 0.0456 | 0.0670 | 0.0878 | 0.0750 |

Systems . | Parameters . | $tnl4$ . | $tnl5$ (Quadratic) . | $tnl6$ (Liu–Parr) . | $tnl5,scalingx$ (Coordinate scaling) . |
---|---|---|---|---|---|

H_{8} | R^{2} | 0.9941 | 0.9996 | 0.9898 | 0.9953 |

RMSE | 0.0029 | 0.0015 | 0.0038 | 0.0025 | |

T_{KSDFT} (Ha) | 1.1147 | 1.1145 | 1.1145 | 1.1147 | |

T_{model} (Ha) | 1.1136 | 1.1152 | 1.1202 | 1.1651 | |

T_{error} (Ha) | 0.0011 | 0.0007 | 0.0057 | 0.0504 | |

LiH | R^{2} | 0.9984 | 0.9987 | 0.9982 | 0.9981 |

RMSE | 0.0038 | 0.0035 | 0.0040 | 0.0041 | |

T_{KSDFT} (Ha) | 0.7583 | 0.7428 | 0.8068 | 0.7338 | |

T_{model} (Ha) | 0.7651 | 0.7642 | 0.8566 | 0.7592 | |

T_{error} (Ha) | 0.0075 | 0.0061 | 0.0498 | 0.0253 | |

LiF | R^{2} | 0.9999 | 0.9999 | 0.9994 | 0.9998 |

RMSE | 0.0043 | 0.0031 | 0.0087 | 0.0065 | |

T_{KSDFT} (Ha) | 6.7109 | 7.0165 | 6.7263 | 6.9105 | |

T_{model} (Ha) | 6.7355 | 7.0304 | 6.7752 | 7.0209 | |

T_{error} (Ha) | 0.0246 | 0.0138 | 0.0489 | 0.1103 | |

C_{4}H_{2} | R^{2} | 0.9996 | 0.9996 | 0.9986 | 0.9973 |

RMSE | 0.0069 | 0.0075 | 0.0100 | 0.0185 | |

T_{KSDFT} (Ha) | 13.1728 | 12.4828 | 13.0915 | 13.4613 | |

T_{model} (Ha) | 13.2473 | 12.4428 | 13.1402 | 13.5490 | |

T_{error} (Ha) | 0.0745 | 0.0399 | 0.0995 | 0.0878 | |

C_{4}N_{2} | R^{2} | 0.9999 | 0.9999 | 0.9993 | 0.9997 |

RMSE | 0.0031 | 0.0039 | 0.0089 | 0.0065 | |

T_{KSDFT} (Ha) | 20.7398 | 20.2343 | 20.7399 | 20.6913 | |

T_{model} (Ha) | 20.6410 | 20.2450 | 20.8654 | 20.7679 | |

T_{error} (Ha) | 0.0987 | 0.0107 | 0.1834 | 0.0765 | |

C_{3}O_{2} | R^{2} | 0.9999 | 0.9998 | 0.9995 | 0.9998 |

RMSE | 0.0030 | 0.0039 | 0.0066 | 0.0061 | |

T_{KSDFT} (Ha) | 19.3773 | 19.3775 | 19.2867 | 19.3338 | |

T_{model} (Ha) | 19.3316 | 19.4445 | 19.1988 | 19.4088 | |

T_{error} (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 Manzhos^{43,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.

### C. Behavior outside the training region: Performance for charged systems

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, C_{4}N_{2} and H_{8}, obtained from KSDFT calculations using the exact-exchange functional and from the three models (proposed in Sec. II) trained on three neutral structures of C_{4}N_{2} and H_{8}. 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.

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 C_{4}N_{2}. 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 H_{8}, but are about 5% in the case of C_{4}N_{2}, 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.

### D. A model for all systems

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 R^{2} 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 R^{2} (more than 0.99) except for H_{8} and LiH.

. | Systems/metrics . | C_{4}N_{2}
. | C_{4}H_{2}
. | C_{3}O_{2}
. | H_{8}
. | LiH . | LiF . |
---|---|---|---|---|---|---|---|

Universal model 1 | R^{2} | 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 | |

T_{KSDFT} (Ha) | 20.4839 | 19.1199 | 13.1749 | 0.4397 | 6.7529 | 0.8068 | |

T_{model} (Ha) | 20.4454 | 19.1496 | 13.2281 | 0.2873 | 6.9459 | 0.8486 | |

T_{error} (Ha) | 0.0385 | 0.0297 | 0.0532 | 0.1524 | 0.1934 | 0.0418 | |

Universal model 2 | R^{2} | 0.9999 | 0.9999 | 0.9998 | ⋯ | ⋯ | ⋯ |

RMSE (Ha/bohr) | 0.0028 | 0.0022 | 0.0034 | ⋯ | ⋯ | ⋯ | |

T_{KSDFT} (Ha) | 20.4745 | 19.1213 | 13.1876 | ⋯ | ⋯ | ⋯ | |

T_{model} (Ha) | 20.4842 | 19.1289 | 13.1703 | ⋯ | ⋯ | ⋯ | |

T_{error} (Ha) | 0.0098 | 0.0076 | 0.0172 | ⋯ | ⋯ | ⋯ | |

Universal model 3 | R^{2} | ⋯ | ⋯ | ⋯ | 0.9976 | 0.9999 | 0.9999 |

RMSE (Ha/bohr) | ⋯ | ⋯ | ⋯ | 0.0089 | 0.0056 | 0.0049 | |

T_{KSDFT} (Ha) | ⋯ | ⋯ | ⋯ | 0.4410 | 6.8143 | 0.7986 | |

T_{model} (Ha) | ⋯ | ⋯ | ⋯ | 0.4841 | 6.8634 | 0.8162 | |

T_{error} (Ha) | ⋯ | ⋯ | ⋯ | 0.0431 | 0.0490 | 0.0176 |

. | Systems/metrics . | C_{4}N_{2}
. | C_{4}H_{2}
. | C_{3}O_{2}
. | H_{8}
. | LiH . | LiF . |
---|---|---|---|---|---|---|---|

Universal model 1 | R^{2} | 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 | |

T_{KSDFT} (Ha) | 20.4839 | 19.1199 | 13.1749 | 0.4397 | 6.7529 | 0.8068 | |

T_{model} (Ha) | 20.4454 | 19.1496 | 13.2281 | 0.2873 | 6.9459 | 0.8486 | |

T_{error} (Ha) | 0.0385 | 0.0297 | 0.0532 | 0.1524 | 0.1934 | 0.0418 | |

Universal model 2 | R^{2} | 0.9999 | 0.9999 | 0.9998 | ⋯ | ⋯ | ⋯ |

RMSE (Ha/bohr) | 0.0028 | 0.0022 | 0.0034 | ⋯ | ⋯ | ⋯ | |

T_{KSDFT} (Ha) | 20.4745 | 19.1213 | 13.1876 | ⋯ | ⋯ | ⋯ | |

T_{model} (Ha) | 20.4842 | 19.1289 | 13.1703 | ⋯ | ⋯ | ⋯ | |

T_{error} (Ha) | 0.0098 | 0.0076 | 0.0172 | ⋯ | ⋯ | ⋯ | |

Universal model 3 | R^{2} | ⋯ | ⋯ | ⋯ | 0.9976 | 0.9999 | 0.9999 |

RMSE (Ha/bohr) | ⋯ | ⋯ | ⋯ | 0.0089 | 0.0056 | 0.0049 | |

T_{KSDFT} (Ha) | ⋯ | ⋯ | ⋯ | 0.4410 | 6.8143 | 0.7986 | |

T_{model} (Ha) | ⋯ | ⋯ | ⋯ | 0.4841 | 6.8634 | 0.8162 | |

T_{error} (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 C_{4}N_{2}, C_{4}H_{2}, and C_{3}O_{2}. However, deviations from one-dimensional model KSDFT results are clearly evident in H_{8} 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 H_{8} 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 H_{8} 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.

Interestingly, we also observe large errors in the total kinetic energy for LiH (0.19 Ha) and H_{8} (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 H_{8} and universal model 3 with structures from C_{4}H_{2}, C_{4}N_{2}, and C_{3}O_{2}. 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 *Q*_{1} + *Q*_{3} model proposed in Ref. 15. However, for the 162 structures in the test set, the RMSE and R^{2} 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.

### E. Model for modified kinetic energy density

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\u0304KE\rho =Ts[\rho ]+EX[\rho ({\psi k})]\u2212ELDA[\rho ]$. 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

where $T\u0304KE$ is the modified kinetic energy for which we want to develop a density-dependent model, *T*_{s} is the Kohn–Sham kinetic energy, *E*_{X} is the exact-exchange energy, and *E*_{LDA} 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 *R*^{2} 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.

Systems . | Parameters . | Individual model . | Universal model 1 . | Universal model 2 . | Systems . | Parameters . | Individual model . | Universal model 1 . | Universal model 3 . |
---|---|---|---|---|---|---|---|---|---|

H_{8} | R^{2} | 0.9914 | 0.8817 | 0.9876 | C_{4}H_{2} | R^{2} | 0.9989 | 0.9984 | 0.9992 |

RMSE | 0.0039 | 0.0140 | 0.0054 | RMSE | 0.0074 | 0.0086 | 0.0081 | ||

$T\u0304KSDFT$ (Ha) | −1.7549 | −1.1145 | −1.243 | $T\u0304KSDFT$ (Ha) | 3.1179 | 3.2522 | 3.1183 | ||

$T\u0304model$ (Ha) | −1.6669 | −1.8698 | −1.3449 | $T\u0304model$ (Ha) | 3.0108 | 3.3805 | 3.2315 | ||

$T\u0304error$ (Ha) | 0.0879 | 0.2696 | 0.1029 | $T\u0304error$ (Ha) | 0.1071 | 0.1282 | 0.1132 | ||

LiH | R^{2} | 0.9969 | 0.9711 | 0.9903 | C_{4}N_{2} | R^{2} | 0.9997 | 0.9992 | 0.9993 |

RMSE | 0.0030 | 0.0093 | 0.0065 | RMSE | 0.0044 | 0.0082 | 0.0051 | ||

$T\u0304KSDFT$ (Ha) | −0.6192 | −0.6193 | −0.6192 | $T\u0304KSDFT$ (Ha) | 7.1003 | 6.8673 | 6.8673 | ||

$T\u0304model$ (Ha) | −0.6063 | −0.5541 | −0.5905 | $T\u0304model$ (Ha) | 7.1551 | 6.7514 | 6.9454 | ||

$T\u0304error$ (Ha) | 0.0129 | 0.0652 | 0.0287 | $T\u0304error$ (Ha) | 0.0548 | 0.1158 | 0.0781 | ||

LiF | R^{2} | 0.9999 | 0.9988 | 0.9994 | C_{3}O_{2} | R^{2} | 0.9996 | 0.9995 | 0.9992 |

RMSE | 0.0035 | 0.0112 | 0.0045 | RMSE | 0.0064 | 0.0069 | 0.0073 | ||

$T\u0304KSDFT$ (Ha) | 2.4865 | 2.4040 | 2.4576 | $T\u0304KSDFT$ (Ha) | 6.7275 | 7.1814 | 6.8234 | ||

$T\u0304model$ (Ha) | 2.5129 | 2.6726 | 2.4963 | $T\u0304model$ (Ha) | 6.6423 | 7.0795 | 6.8357 | ||

$T\u0304error$ (Ha) | 0.0265 | 0.2686 | 0.0387 | $T\u0304error$ (Ha) | 0.0851 | 0.1018 | 0.0123 |

Systems . | Parameters . | Individual model . | Universal model 1 . | Universal model 2 . | Systems . | Parameters . | Individual model . | Universal model 1 . | Universal model 3 . |
---|---|---|---|---|---|---|---|---|---|

H_{8} | R^{2} | 0.9914 | 0.8817 | 0.9876 | C_{4}H_{2} | R^{2} | 0.9989 | 0.9984 | 0.9992 |

RMSE | 0.0039 | 0.0140 | 0.0054 | RMSE | 0.0074 | 0.0086 | 0.0081 | ||

$T\u0304KSDFT$ (Ha) | −1.7549 | −1.1145 | −1.243 | $T\u0304KSDFT$ (Ha) | 3.1179 | 3.2522 | 3.1183 | ||

$T\u0304model$ (Ha) | −1.6669 | −1.8698 | −1.3449 | $T\u0304model$ (Ha) | 3.0108 | 3.3805 | 3.2315 | ||

$T\u0304error$ (Ha) | 0.0879 | 0.2696 | 0.1029 | $T\u0304error$ (Ha) | 0.1071 | 0.1282 | 0.1132 | ||

LiH | R^{2} | 0.9969 | 0.9711 | 0.9903 | C_{4}N_{2} | R^{2} | 0.9997 | 0.9992 | 0.9993 |

RMSE | 0.0030 | 0.0093 | 0.0065 | RMSE | 0.0044 | 0.0082 | 0.0051 | ||

$T\u0304KSDFT$ (Ha) | −0.6192 | −0.6193 | −0.6192 | $T\u0304KSDFT$ (Ha) | 7.1003 | 6.8673 | 6.8673 | ||

$T\u0304model$ (Ha) | −0.6063 | −0.5541 | −0.5905 | $T\u0304model$ (Ha) | 7.1551 | 6.7514 | 6.9454 | ||

$T\u0304error$ (Ha) | 0.0129 | 0.0652 | 0.0287 | $T\u0304error$ (Ha) | 0.0548 | 0.1158 | 0.0781 | ||

LiF | R^{2} | 0.9999 | 0.9988 | 0.9994 | C_{3}O_{2} | R^{2} | 0.9996 | 0.9995 | 0.9992 |

RMSE | 0.0035 | 0.0112 | 0.0045 | RMSE | 0.0064 | 0.0069 | 0.0073 | ||

$T\u0304KSDFT$ (Ha) | 2.4865 | 2.4040 | 2.4576 | $T\u0304KSDFT$ (Ha) | 6.7275 | 7.1814 | 6.8234 | ||

$T\u0304model$ (Ha) | 2.5129 | 2.6726 | 2.4963 | $T\u0304model$ (Ha) | 6.6423 | 7.0795 | 6.8357 | ||

$T\u0304error$ (Ha) | 0.0265 | 0.2686 | 0.0387 | $T\u0304error$ (Ha) | 0.0851 | 0.1018 | 0.0123 |

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 *R*^{2} scores for universal model 1 are similar to the models that were trained for the individual systems except for LiH and H_{8}. 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 R^{2} scores are very high. These results suggest that quadratic models can be used to describe kinetic energy and exact-exchange energy with good accuracy.

## IV. CONCLUSION

In this study, using data from a few model one-dimensional systems, such as H_{8}, C_{4}H_{2}, C_{4}N_{2}, LiH, LiF, and C_{3}O_{2}, 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 R^{2} 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 R^{2} 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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to disclose..

## DATA AVAILABILITY

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

### APPENDIX A: FORMULATION

In this section, we briefly outline formulation and implementation of a one-dimensional model analog of Kohn–Sham DFT^{1,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:

where *ψ*_{i} and *ɛ*_{i} are the model Kohn–Sham eigenfunctions and eigenvalues, respectively, *v*_{ext}(*x*) is the external potential, *v*_{H}[*ρ*](*x*) is the Hartree potential, $VX{\psi j}(x,x\u2032)$ is the exact-exchange operator, *N*_{s} is the total number of occupied orbitals, and *f*_{i} is the occupation number of the *i*th orbital. In our model system, we replace the Coulomb form of Hartree potential and external nuclear potential by a soft form,^{57}

Here, *Z*_{I} and *X*_{I} denote the nuclear charge and position of the *I*th atom, respectively. The exact-exchange operator *V*_{X} is a nonlocal operator that depends on the occupied orbitals ${\psi j}j=1Ns$ and has the following form:^{58,59}

The action of *V*_{X} on a given orbital $\psi x$ is given by^{58,59}

Therefore, the Hamiltonian in the eigenvalue problem outlined in Eq. (A1) depends both on the electron density, $\rho x$, and the orbitals, ${\psi 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 mixing^{60–62} while keeping the exchange operator constant.

### APPENDIX B: IMPLEMENTATION

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}

where *g*^{(i)} denotes the value of the function at the *i*th grid point and *N*_{FD} is the order of finite-difference approximation. The weights, *W*_{k}, associated with the function values are

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

where $\u2207h2$ is the finite-difference discrete second derivative matrix and $\psi ji$ is the value of the *j*th wavefunction at the *i*th 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:

Here, *T*_{s} = ∫*t*_{s}(*x*)*dx* and *t*_{s}(*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.

### APPENDIX C: ANALYTIC EXPRESSION FOR THE DERIVATIVES OF *P*(r, r′)

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,

Here, *d*_{k} are the coefficients and *β*_{k} and *N*_{G} are the widths and the number of the Gaussian functions used to model the function *P*(**r**, **r**′).