{'fraction_non_alphanumeric': 0.031151211801896735, 'fraction_numerical': 0.010603266596417282, 'mean_word_length': 4.9922388845040775, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Machine learning approaches for feature engineering of the crystal structure: Application to the prediction of the formation energy of cubic compounds

In this study, we present a novel approach along with the needed computational strategies for efficient and scalable feature engineering of the crystal structure in compounds of different chemical compositions. This approach utilizes a versatile and extensible framework for the quantification of a three-dimensional (3-D) voxelized crystal structure in the form of 2-point spatial correlations of multiple atomic attributes and performs principal component analysis to extract the low-dimensional features that could be used to build surrogate models for material properties of interest. An application of the proposed feature engineering framework is demonstrated on a case study involving the prediction of the formation energies of crystalline compounds using two vastly different surrogate model building strategies - local Gaussian process regression and neural networks. Specifically, it is shown that the top 25 features (i.e., principal component scores) identified by the proposed framework serve as good regressors for the formation energy of the crystalline substance for both model building strategies.


Introduction
Although physics-based modeling approaches such as density-functional theory (DFT) [1,2] offer the preferred avenue for estimating the physical and chemical properties of crystal structures, they are not ideally suited for materials discovery and innovation efforts. Such efforts demand inverse solutions such as finding the molecular chemistry and structure that meets a targeted combination of physical and chemical properties [3]. One of the most practical strategies for addressing the inverse solutions of materials design is to first produce high-fidelity low-computational cost surrogate models trained to the available data (e.g., collections of DFT computations), and then use the surrogate model for addressing the inverse problems of interest. Such a strategy can prove to be highly beneficial in rapid screening of vast design spaces [3][4][5]. In this context, the emerging toolsets and paradigms in Machine Learning (ML) offer new opportunities for building the surrogate models needed to accelerate the discovery of new crystal structures. In particular, it is now possible to train these models using collections of publicly accessible DFT data in repositories such as the Quantum Machine (QM) [6][7][8][9][10][11][12][13][14][15][16].
ML models have been used in prior work to predict a broad range of material properties based on DFT datasets. The properties modeled have included optical and electronic bandgaps [17,18], formation energy [6,8,[19][20][21][22][23][24], atomization energies [25][26][27][28][29] and polarizability of crystalline compounds [30]. ML models have also been used in prior work for the prediction and classification of crystal compositions and structures [31][32][33][34][35] and in the development of many-body interatomic potentials for atomistic simulations [36][37][38]. One of the foundational elements of ML common to most model building approaches is feature engineering, which aims to identify a small list of transformed input variables (called features) from the original large list of input variables that potentially could influence the predictions of the output variables (called targets). As one would expect, feature engineering governs the accuracy and utility of the surrogate model. In some of the ML approaches (e.g., convolutional neural networks (CNNs)), feature engineering occurs implicitly in the model training phase. In general, the more implicit one makes the feature engineering task, the more training data would be needed. This is because the model needs to learn the salient features first before learning their quantitative relationships with the targets. In much of the prior work [12,21,25,33], large feature lists have been manually cultivated by researchers based largely on their intuition of the atomic physics. Consequently, these feature lists have included various simplified attributes of the chemical compositions (e.g., atomic fractions) as well as the atomic structure (e.g., bond lengths and bond angles). These approaches face significant challenges. First, the creation of these lists has been pursued largely in an ad-hoc manner, which inevitably reflects the bias of the individual researchers. As such, these approaches may fail to account for certain potentially salient features controlling the target property. Second, the specification of chemical composition using distinct labels for the different atomic species hinders learning across crystal structures with different elemental compositions. One strategy to address this challenge has been to supplement the feature lists with certain atomically weighted physical properties [8,39,40]. Additionally, there have also been various efforts at enhancing the specification of the atomic structure in the feature lists by retaining much of the actual three-dimensional (3-D) atomic structure information. For instance, Faber et al. [20] utilized an extension of the Coulomb matrix representation to account for the periodicity of the crystal structure. In a different approach, Ward et al. [23] represented the crystal structure with a Voronoi decomposition and employed local atomic species-level descriptors between neighboring species in addition to global structure descriptors such as maximum packing efficiency. In similar efforts, Schutt et al. and Honrao et al. used a feature set containing values of the partial Radial Distribution Function (RDF) between pairs of constituent species in the crystal structure [41,42]. The partial RDF is a measure of the frequency of a pair of species separated by a certain distance within the material volume. However, the partial RDF does not contain directional information about the distribution of the atoms. One approach to overcome this limitation has been to supplement the partial RDF with an Angular Distribution Function (ADF) that contains information on the distribution of inter-atomic angles between species [43]. Nevertheless, this approach still utilizes distinct labels for constituent species, and is thus typically useful only for description of crystal structures with few species (such as Al-Ni and Cd-Te systems [43]). Another interesting approach proposed recently relied on the development of graph embeddings that encode species-specific features (such as the periodic table group number, period number, electronegativity) and a limited set of structural features (such as the interatomic distance) as nodes and edges of a multigraph, respectively [24,44,45]. This representation is amenable to usage in graph convolutional neural networks [46,47] that sequentially build up localized representations for each node by iteratively including information from neighboring nodes.
Using the feature engineering schemes described above, current efforts have employed various regression approaches such as kernel ridge regression (KRR) [20,42], support vector machines [42,43], deep neural networks [21], gradient boosted decision trees [8], and graph convolutional neural networks [24,44,45] to build the desired reduced-order models. While the regression-based learning methods are computationally very efficient, they are often prone to over-fitting because they invariably employ a large number of implicit model parameters that need to be trained on a large collection of ground-truth data. The overfitting may be mitigated to a certain extent with the adoption of well-known techniques such as early stopping [48], loss function regularization and dropout [49,50]. With the adoption of such techniques, the high computational efficiency of neural networks (NNs) makes them an attractive option for building reduced-order models. However, a significant challenge arises from the availability of relatively small training datasets due to the very high cost of DFT computations. The role of the feature engineering procedure is especially critical in such smaller sized datasets because of the need to restrict the number of model fit parameters that can be used in these cases.
ML approaches based on Bayesian inference offer an alternate option that is likely to prove beneficial for building surrogate models from DFT computations. More specifically, Gaussian process regression (GPR) [51] offers a powerful non-parametric Bayesian approach to building surrogate models from small training datasets. There are many potential benefits to the utilization of Gaussian Processes in the context of regression. First, they allow a formal treatment of uncertainty in the model predictions. In other words, they do not just estimate the expected values for the model outputs, but also their distributions. Second, the formal treatment of model uncertainty allows the design and implementation of strategies that could potentially reduce the effort spent in the generation of the training data. More specifically, it is possible to formulate and maximize the expected information gain to the surrogate model with the addition of each specific new training data point. This is particularly important to building surrogate models with limited data from the computationally expensive DFT datasets. Furthermore, GPR models are typically more interpretable than NNs.
In this paper, we present a systematic and comprehensive approach to feature engineering of crystal structure that directly addresses the challenges described earlier. In this new approach, we systematically and efficiently assemble an extremely large number of spatial correlations (specifically, 2-point correlations) [5,52] that implicitly account for atomic attributes (e.g., Pauling electronegativity, ionization energy, heat of fusion). After processing the large feature list for the complete ensemble of crystal structures of interest, we employ principal component analysis (PCA) to obtain a suitable low-dimensional representation of the feature list. This strategy offers many advantages compared to the approaches used in current literature. First, the approach presented in this work has the potential to systematically generate a very large set of physics-inspired features for the rigorous and comprehensive quantification of the crystal structure. Second, the proposed protocols utilize digital (i.e., voxelized) representations of the 3-D crystal structure [53,54] along with compact Fourier representations and associated computational algorithms (e.g., Fast Fourier transform (FFT)) for highly efficient computation of the features (i.e., spatial correlations). Third, the encoding of individual chemical species through their individual physical properties enables effective learning across crystal structures of different elemental compositions. Fourth, the use of PCA provides an objective (i.e., data-driven) path for establishing low-dimensional representations of the compound crystal structure that can then be used with a broad variety of model-building approaches. Fifth, the features identified by the proposed framework are independent of the target property predicted by the model.
The primary goal of this paper is to develop and demonstrate a versatile feature engineering methodology for materials problems involving different chemical compositions of compounds and their crystal structures. This methodology is aimed at extracting reliable models from small data sets using physics-inspired feature engineering approaches. We demonstrate the utility of this feature engineering scheme by building predictive models for the crystal formation energy using two drastically different model-building approaches: (i) a localized variant of GPR [55], and (ii) a feedforward neural network (NN). A second goal of this work is to critically compare the relative merits of both model building techniques for addressing materials problems.

Background
The central goal of the feature engineering step is to establish a compact set of salient inputs (i.e., features) that serve as suitable predictors for the selected targets. In the context of our problem, this would entail establishing a set of salient crystal structure descriptors for capturing high-fidelity reduced-order Structure-Property (SP) relationships [5,56]. In this effort, the formalism of npoint spatial correlations (also called n-point statistics) [5,[57][58][59][60][61] offers a systematic approach to statistical quantification of the underlying morphological patterns in the heterogeneous material internal structure. Several existing statistical atomic physics models such as the Ising model (for predicting ferromagnetism) [62] and the Potts model (a general model of interacting atomic spins) [63] predict the bulk properties of materials as a sum of contributions arising from local structure features that can be easily interpreted as components of the n-point spatial correlations. Spatial correlations are indeed the features dictated by the governing physics in studies of heterogeneous material structure at all length scales, spanning from the atomistic to the mesoscale. However, the comprehensive set of n-point spatial correlations is typically too large and unwieldy to serve directly as inputs for the reduced-order models. As such, there is a critical need for a compact representation of the material structure features that can serve effectively as inputs to produce the desired highfidelity reduced-order models. In this section, we present a brief overview of the concepts of spatial correlations, salient feature extraction (using PCA), and the strategies for building reduced-order models (i.e., GPR and NNs) used in this work.

Spatial Correlations
Typically, efficient computations of spatial correlations for the quantification of the heterogeneity of the material structure utilize digital representations (i.e., voxelized representations with suitable assignments of material states to each voxel) and FFT algorithms for computing the convolution operations involved in these computations [5]. Mathematically, these voxelized representations can be expressed as a high-dimensional array m p s , whose elements capture the value of a local feature (reflecting the local material state) indexed by p in a voxel indexed by s (defined as vector of integers {s 1 , s 2 , s 3 } for 3-D material volumes). This array of size (P × n(S)), with P denoting the number of distinct material states assigned to each voxel in the material volume and n(S) denoting the cardinality of the set of spatial voxels S, defines one instantiation of a material structure.
Our interest lies in extracting microstructure statistics from each instantiation (i.e., m p s array) that can serve as features in the formulation of surrogate models connecting material structure to the properties/performance characteristics of interest. The most systematic set of such microstructure statistics are provided by the formalism of n-point spatial correlations, which capture the relevant statistics related to the morphological details (i.e., size and shape distributions) of the distinct local states in the material structure. The most basic set of these statistics are the 2-point correlations denoted by f pq r , which track how often distinct local states p and q are separated by a discretized vector indexed by r. Mathematically, 2-point correlations may be computed as [5,61,64] One way to interpret the above definition is to recognize that the numerator denotes the number of successes in locating the local states p and q separated by the vector index r, and that the denominator represents the total number of valid trials (for periodic microstructures, this is equal to the number of voxels). It should be noted that the RDF (used commonly in the quantification of the molecular structures) [59,65] is a particular variant of the 2-point correlations described in Eq.
(1) in which one does not consider the orientation of the vector, but only its magnitude. Noting that the central operation in Eq. (1) is essentially a discrete cross-correlation over the spatial domain of voxels S, the 2-point correlations are efficiently computed by taking advantage of FFT algorithms [5,52]. The computations naturally exploit the periodicity of the crystal structures. There also exist several redundancies that can be leveraged for obtaining a more compact set of spatial correlations.
Specifically, it can be seen that Eq. (1) produces P 2 sets of spatial correlations. It has been shown that only P of these are adequate to compute the rest of the spatial correlations [61,66]. As a result, it is often adequate to compute the autocorrelations (which correspond to n = p in Eq. (1)) for the dominant local state in the collected ensemble of material structures and its cross-correlations with the rest of the local states.

Extraction of salient features
In prior work, our research group has established a generalized framework referred as Materials The benefits of the MKS pipeline have been expounded in prior work [5,56,[69][70][71][72][73]. First, this approach yields an objective low-dimensional representation of the material structure that is not influenced by either the materials manufacturing process or property information, thereby providing a consistent representation of the material structure in process-structure-property (PSP) linkages [5]. Second, since PCA is essentially a Fourier representation, it allows for approximate but optimal reconstructions of the material structure statistics (controlled by the truncation levels applied in retaining the PC scores), which in turn can then be used with sophisticated algorithms [74,75] for the statistical reconstruction of the material structure. Third, the concepts based on spatial correlations and PCA are broadly applicable to virtually all different classes of material internal structures found at multiple hierarchical length scales, spanning from the atomic to the macroscale.
Fourth, and perhaps the most importantly, it has been seen that only a handful of PC scores (≈ 5 to 10) are often adequate in producing high fidelity PSP linkages needed to drive materials innovation efforts [5,70,76]. It should be noted that the use of the PC representations of the spatial correlations as the features of the material structure has thus far been explored mostly at the meso-length scales [5,70,73,77]. This concept is only now beginning to be explored at the atomic structure length scales.

Reduced-order model building strategies
In the final step of the MKS framework, the reduced-order features of interest are correlated with the target property using a variety of model-building strategies. NNs provide a nonlinear modeling approach to accomplish this task and are known to be sufficient to capture any arbitrary mapping between the inputs and the output. The most basic feedforward NNs consist of multiple fully connected layers of multiple neurons, with each neuron capturing a linear transformation followed by a nonlinear activation (e.g., ReLu, sigmoid function) [78]. The architecture of a typical feedforward NN is shown in Fig. 1. The weights and the bias (i.e., model fitting parameters) associated with each neuron are calibrated with backpropagation. This is accomplished through minimization of a user-specified loss function between the predicted output and the corresponding observed values in the training data. NNs derive their scalability from the computationally efficient algorithms used for the calibration of the model parameters. These are readily accessible in many software packages (e.g., PyTorch [79], TensorFlow [80]), most of which are amenable to computation on Graphics Processing Units (GPUs). NNs also typically iterate through multiple epochs (i.e., passes over the entire training data set) in order to optimize the model parameters. The main limitation of the NNs is that they need a substantially large training dataset, without which they are most likely to produce a model overfit (i.e., significantly larger errors for test data points compared to the training data points). Overfits to training data often occur due to an over-parametrization of the NN (i.e., more trainable weights and biases to be calibrated than the size of the dataset). In this work, we utilize feedforward NNs with two hidden layers to predict the formation energies of the compounds in our dataset. For an NN with two hidden layers consisting of h 1 and h 2 neurons, respectively, the number of trainable model parameters to be calibrated (denoted by C) for a single output from an N-dimensional input is given by An alternate modeling technique that is better suited for small datasets is GPR [51,81,82], which offers a nonparametric Bayesian approach to building surrogate models. GPR (and its variants) have successfully been applied to model structure-property linkages in the mesoscale from relatively small datasets [70,83,84]. GPR models the target as a Gaussian Process (GP) that is fully defined by the specification of a mean and a covariance. The GP is then tuned using conditional distributions defined on available training data. Let X, X * and y denote the N × D matrix of training data points, N * × D matrix of test data points and N × 1 vector of target values in the training set, respectively. Additionally, let K(X, X ), K(X * , X * ) and K(X, X * ) denote the N ×N , N * ×N * and N × N * covariance matrices assembled (using a kernel function described later) using the inputs in the training dataset, test dataset, and between the training and test datasets, respectively. In GPR, the predictive mean and variance for test points is expressed as where denotes the transpose of a matrix. As stated earlier, one typically uses a suitable kernel function to compute the various covariance matrices in Eqs. 3 and 4. One of the most commonly used kernel functions is the ARD-SE (Automatic Relevance Determination -Squared Exponential) [51,85] function expressed as where x and x denote any two input vectors selected from the data matrices for which the covariance matrix is being computed, subscript d corresponds to the d th feature in the input vector, σ s is a scaling parameter that controls the scaling of the output variance, θ d is the correlation-length hyperparameter corresponding to the d th feature in the input vector, σ n denotes the noise in the target, and δ xx denotes the Kronecker delta. The treatment of the noise using a single parameter added to the diagonal of the covariance indicates independence of noise with the input data (i.e., homoscedasticity) [51]. The central benefit of the ARD-SE kernel is that it allows us to tune independently the correlation-length hyperparameter for each of the input features. This allows for better interpretability of the model, because the salient features with relatively smaller length scale hyperparameters exhibit a higher sensitivity to the predicted target. One should also note that extremely small values of the hyperparameters tend to make the predictions very noisy. Therefore, it is important to tune the values of the hyperparameters in order to produce robust and reliable predictions. This is typically accomplished with Maximum Likelihood Estimation (MLE), which lacks a closed-form solution; consequently, one must resort to iterative schemes such as quasi-Newton algorithms to optimize the hyperparameters [86][87][88].
One of the central challenges in the implementation of GPR comes from the computational cost, especially with large training datasets. The main computational bottleneck arises from the computation of the inverse of K(X, X ) in Eq. 3, which typically scales as O(N 3 ) for dense matrices. This makes traditional GPR intractable even for moderately sized datasets (N > 1, 000).
This difficulty is compounded by the need to optimize the large number of hyperparameters (D + 2) in the ARD-SE kernel. Several methods have been suggested in literature to address these challenges.
These include the utilization of low-rank sparse approximations for the covariance matrix [51,89,90], tree-based partitioning to develop smaller datasets and fit individual GP predictors [91], and Localized GPR (L-GPR) proposed by Gramacy and Apley (discussed next and used in this study) [92].
L-GPR involves identifying a local subset of training data points to construct a separate GP for each predictive point. In addition to allowing the application of GPR to larger datasets, a key advantage of L-GPR is the ability to accommodate non-stationarity in the model by allowing optimization of the interpolation hyper-parameters in the kernel functions depending on the location of the test point in the high-dimensional input space. Different criteria may be considered while selecting the local training points; an intuitive but suboptimal method would be using a certain number of nearest neighbors as the local neighborhood [93]. However, a more informative design criterion would be to sequentially build the local neighborhood by evaluating a tradeoff between expected information gain by adding a training point and the increase in the predictive variance of the GPR. This method of local training point selection promotes an efficient exploration of the sample space and provides a natural guidance to avoid overfitting of the models. This approach is referred as Active Learning Cohn (ALC) [55,94], and is utilized in this work.

3. Novel framework for quantifying crystal structures
In this section, we describe a new framework for effective low-dimensional representation (i.e., fingerprints) of the crystal structure of the compounds. This is accomplished by first establishing computationally efficient protocols to arrive at voxelized representations of the atomic structure of compounds, and then suitably extending and applying the current MKS framework (described in Sections 2.1 and 2.2) on such voxelized representations. In prior work [74,95], very simple descriptors were utilized for the microstructure array, m n s , where it was assigned a value of one for voxels within the atomic volume and a value of zero otherwise. A key challenge encountered in this simple representation is that it does not allow an efficient interpolation between material structures (i.e., compounds) of varying chemical compositions. In this work, we expand the previous framework to allow for an enhanced representation of multiple atomic attributes of interest in each voxel. This is accomplished by representing the local material state as a vector set of attributes that take continuous values. Such continuous local states have been utilized previously in the treatment of mesoscale polycrystalline microstructures and mesoscale fields of chemical compositions [5,96]. In this work, we will extend these representations to the quantification of crystal structures of compounds.
Choudhary et al. [8] have studied the order of importance of various atomic attributes and their role in the formation energy of a compound. Specifically, they have identified ionization energy, Pauling electronegativity and heat of fusion as the most important attributes. Therefore, in this work, we define the local material state in the voxel indexed by s using a vector of attributes a 0s , a 1s , a 2s , v s , where a 0s , a 1s and a 2s denote the heat of fusion, ionization energy, and Pauling electronegativity, respectively, and v is a binary attribute that takes a value of zero for voxels within any of the atomic volumes (defined by a sphere with radius equal to the atomic radius of the species) and a value one otherwise. Our model building efforts (described later) found it beneficial to further enhance the attribute list using various monomials of the atomic attributes. In other words, the full list of local states included in the atomic state vector m s may be defined as 0s , a 1s , a 2s , a 0s a 1s , a 1s a 2s , a 0s a 2s , a 0s a 1s a 2s , a 2 0s , a 2 0s a 1s , . . . , a l0 0s a l1 1s a l2 2s , v s where l 0 , l 1 and l 2 represent the maximum degree considered for each of the three atomic attributes.
A total of (l 0 + 1) × (l 1 + 1) × (l 2 + 1) local states were used in the definition of the atomic state vector. Through multiple trials, it was found that the combination of l 0 , l 1 and l 2 set to 5, 5 and 2, respectively, yielded the optimal performance in terms of the tradeoff between computational cost (which scales with the number of local states), and the accuracy of our models. For these selections, the total number of local states in m s is 108. Elements of m s are then used in the computation of the 2-point correlations using Eq. (1). As already noted, the convolution operation involved is best accomplished using the FFT algorithm. Mathematically, this is expressed as where M p k and F pq k denote the discrete Fourier transforms of m p s and f pq r , respectively, and the superscript * denotes the complex conjugate. It is important to note that Eq. (7) implicitly assumes the fields involved are spatially periodic. Although all the compounds considered in this work exhibit periodicity (i.e., all spatial fields defined based on the atomic structure of the compound are inherently spatially periodic), the sizes of their periodic unit cells vary substantially. For example, the lattice parameter (reflecting the length scale of the periodicity) for the compounds included in this study varies over the range of 1.8Å to 12.4Å. Therefore, if one were computing the spatial correlations using voxelization schemes based on the respective unit cell sizes, Eq. (7) can be applied directly without any problems. However, a constant voxel size is essential for obtaining the desired low-dimensional representations using PCA in the later steps. Since the unit cell size L of the different compounds considered in this study is unlikely to be an integer multiplier of a common voxel size, t, a suitable strategy is needed to efficiently compute the spatial correlations using a standardized voxel size. Given the size of the dataset, these computations are only practical if we can continue to exploit the cost savings provided by the FFT algorithm in evaluating the spatial correlations. A suitable computational protocol has been devised to address this challenge, which is described next.
This new protocol builds on prior efforts involving mesoscale material microstructures [5].
The computational scheme presented in this work is designed to use a single voxel size and a common set of corresponding vectors indexed by r, for the computation of sets of 2-point correlations (i.e., f pq r ) using Eqs. (1) and (7) for all the compounds considered in this study. For 3-D volumes, it is convenient to represent r itself as a vector of integers r 1 , r 2 , r 3 . The 2-point correlations used in this study are computed only for the set of positive vectors up to a selected maximum length for each of the three components, i.e., {r|0 ≤ r 1 ≤ r max , 0 ≤ r 2 ≤ r max , 0 ≤ r 3 ≤ r max }.
Therefore, the total number of discrete vectors for which we seek to compute f pq r is r 3 max . Our approach to efficiently compute these correlations considers the numerator and denominator in Eq.
(1) separately. As mentioned earlier, for each selected r, the numerator represents the number of successes in locating the local states p and q separated by the selected vector. Our strategy to accurately (and efficiently) compute this quantity is to express the numerator as a cross-correlation between two different material structure spatial fields,m p s and m q s , both of which are derived from the original material structure. The construction of these structure fields is illustrated in Fig. 2 for a simple crystalline compound comprising two chemical species A and B. Note that all structure fields shown in Fig. 2 are 2-D sections of the 3-D material volume of this compound. Fig. 2a shows the unit cell of the original compound structure within a volume that is closely approximated by an integer number of voxels of the selected size, t. Therefore, the voxelated structure shown in Fig. 2a is no longer exactly periodic. In fact, the size of this structure is L t + 0.5 t, where · indicates the greatest integer function (also referred to as the floor function). Clearly, the size of the voxelized unit cell shown in Fig. 2a is not equal to L. This discrepancy introduces a small error in the computation of the spatial correlations that scales with the value of t. In other words, one would have to select the value of t carefully, as a compromise between minimizing the voxelization error and keeping the computational cost reasonable.
The structure fieldsm p s and m q s are defined over larger spatial domains as shown in Figs. 2b and 2c. The size of these fields is specified as L t + 0.5 t + r max t, while the corresponding set of voxels is labelled as S * . Note that the fields are extended so that when any of the vectors of interest (identified  Fig. 3b) reflects the void volume fraction in the crystal structure. Since the void state in the present application is set to either zero or one in each voxel, the autocorrelation map in Fig. 3a actually provides statistical information for a large set of vectors, whose heads can be selected anywhere in the depicted vector space, but tails fixed at 0, 0, 0 . As a specific example, the pink-colored vector in Fig. 3b corresponds to 3.8, 0, 0 Å , with a void autocorrelation of 0.58. This autocorrelation value simply reflects the probability of finding two void (green) voxels in the crystal structure shown in Fig. 3a separated by the selected vector. Fig. 3a shows two example placements of the selected vector -a successful placement (pink vector with a black outline) connecting two void (green) voxels and an unsuccessful placement (pink vector without a black outline). Due to the periodicity of the crystal structure, it is seen that the autocorrelation value for the pink-colored vector is the same as autocorrelation for the zero vector (i.e., the void volume fraction). The crosscorrelation map shown in Fig. 3c similarly captures the spatial correlations between states v s and a 2 s present in the structure. Note that since we do not allow the void state and any material state to co-exist in a single voxel, the cross-correlation for the zero vector exhibits a zero value. A cross-correlation peak is observed for 0, 0, 1.8 Å in Fig. 3c, shown using a blue-colored vector. As before, two example placements of this vector are shown in Fig. 3a, with one reflecting a successful sampling of the desired spatial correlation.

Dataset
A dataset of compound crystal structures and formation energies was obtained from the publicly accessible JARVIS-DFT repository [7,[9][10][11][12]. The extracted data included the crystal structure information (i.e., Bravais lattice type and the coordinates of the atoms present in the unit cell) and their DFT-computed chemical properties. Specifically, the formation energies of ∼ 1740 cubic crystalline compounds that were computed with the OptB88vdW functional [97] were extracted from this repository, and were used to train our models and evaluate the utility of our feature  The feature engineering framework developed in Section 3 is applied to the selected dataset.
As described in section 2.1, the complete set of P 2 (or 108 2 ) spatial statistics contains several redundancies. Building on prior work [refs], the following sets of spatial correlations were included in this work: (i) autocorrelations of v s , a 0s and a 1s (this were the most dominant signals in the dataset), and (ii) the cross-correlations of these three local states with the rest of the local states, while eliminating the trivially related sets (the f pq r are trivially related to f qp r ). This resulted in a total of 321 sets of spatial correlations for each crystal compound. The dimensionality reduction procedure described in section 2.2 was then applied on the complete set of spatial statistics to obtain 25 PC scores for each compound. Since each PC basis in our application carries (weighted) information from a total of ∼ 3 million 2-point statistics, its precise interpretation is currently impractical.

Results and discussion
As mentioned earlier, two different modeling strategies have been explored in this study to evaluate critically the efficacy our feature engineering framework. These modeling strategies were specifically chosen to represent distinctly different approaches to building surrogate models. Therefore, our main goal here is to evaluate the performance of the proposed feature engineering framework (involving PCA on a large feature vector of spatial correlations that account rigorously for the atomic details in the crystal structure) on distinctly different model building strategies. The selected model building strategies for this study include L-GPR (implemented using the laGP package in the R programming language [55]) and NN (implemented using PyTorch [79]). The inputs and the output to both models were the set of the top 25 PC scores (scaled to unit variance) generated using the feature engineering framework presented in this work and the DFT-computed formation energy, respectively.
Both modeling strategies benefitted significantly from the scaling of inputs (to unit variance), due to their usage of gradient-based optimization in the regression procedure (i.e., hyperparameter tuning in L-GPR, and loss function minimization in NNs).
In this study, Mean Absolute Error (MAE) and Median Absolute Error (MedAE) have been chosen as error metrics for quantifying the predictive accuracy of the surrogate models. For a set of J predictive test samples, these error measures are expressed as M edAE = M edian(|y actual |, |y (2) pred − y actual |, . . . , |y in which y These are described next.

L-GPR modeling approach
In this study, five different models were built using the L-GPR strategy. Recall that in this strategy, a separate GP is formulated for each point in the dataset using training points selected from the local neighborhood based on the ALC criterion described in Section 2.3. Since predictions were made at each point independently, the M AE and M edAE were computed directly with Eqs. (9) and (10) Table 1 summarizes the accuracies of these models, as evaluated by the error metrics described previously. In this table, the value of χ is also shown as a percentage of the total number of points in the dataset. The main purpose of this analysis was to understand the effectiveness of the ALC criterion in the L-GPR scheme for selecting optimal training points for the desired model.
Since GPR is essentially an interpolation strategy, L-GPR offers significant computational savings compared to the conventional GPR strategy (which utilizes all points in the dataset). This is because the computational cost of the GPR scales as O(N 3 ), as mentioned previously. Also, the use of ALC in the L-GPR to select the training points offers an organic strategy for regularization and avoiding over-fit. More specifically, since the ALC criterion selects neighborhood points based on the goal of maximizing information gain, one would expect that after the training dataset has included most of the "important" points, there would only be an incremental improvement in performance with a further increase in χ. This tendency was indeed observed in our results. Table 1    in the range of (-2 to 0) eV/atom. This range is also the region with the highest concentration of training data points (see Fig. 4b). As expected, the accuracy of the L-GPR models built is quite sensitive to the amount of available training data close to test data points. For example, the average error for the samples with DFT-computed formation energies greater than 2 eV/atom is quite high, as there are only 20 data points with the corresponding formation energies. The predictive uncertainty of the L-GPR model also followed a similar distribution over these ranges of formation energies, with the highest uncertainty associated with the crystals with DFT-computed formation energies greater than 2 eV/atom. This indicates that prediction confidence in regions with more data is significantly higher than in the regions with sparse data.

NN modeling approach
Five different feedforward NN models (all consisting of two hidden layers) were evaluated in this study. For each model, the loss function (to be minimized) was chosen as L1 loss [98], which reflects the mean absolute error in the prediction of the target for the training dataset. This minimization was performed using the Adam gradient-based optimizer [99] with a learning rate set to 5e − 5.    Figure  6).  Table 2).

Comparison of the L-GPR and NN surrogate model building strategies
While the accuracy of both modeling strategies explored in this study are reported using the same metrics (MAE and MedAE), the performance of these models are not easily compared directly be- This correlated with the use of ∼ 500 nodes for almost all estimators built, indicating that a large number of parameters were employed in the final model. While the test MAE reported in this study is lower, the feature engineering procedure used in this study is significantly more computationally intensive and difficult to scale to larger datasets.
It is emphasized that in comparison with the prior studies, the strategy developed in this work relies on unsupervised (i.e., independent of the selected output variable) feature engineering to identify a small number of salient features (PC scores). These features are then used as inputs to models with relatively fewer fit parameters. This confirms the value and utility of the feature engineering protocols presented in this work. Both models produced in this study are likely to find uses in potentially different applications. The main strength of the L-GPR strategy is that it can provide objective guidance on what new training data points should be generated to improve the model fidelity. This is because the GPR models not only predict the expected values of the output for test data points, but also their variance. This attribute can be used in a suitable design of experiments strategy [51,100] to identify which specific inputs exhibit the highest potential for model improvement. This is particularly useful in the initial stages of generating the training data sets, especially in situations where the cost of data generation is high (e.g., data generated using DFT models). In contrast, NNs could become the preferred approach after a substantial amount of data has been collected. This is mainly because of their ability to cover a much richer space of model functions, and their superior computational efficiency in handling large data sets. The data set used here is small enough that the L-GPR strategy is likely to be the better option for its size.
The overall predictive accuracy obtained by the models in this work is not sufficient for the accurate computation of the crystal formation energy. We believe that the model accuracy can be improved by (i) implementing the feature engineering scheme on a larger dataset (by lifting the constraint of cubic crystal symmetry), and (ii) using the obtained spatial correlations directly as inputs to other model architectures like CNNs. Larger training datasets would allow for models with a larger number of model fit parameters, and thereby improve the model accuracy.

Conclusion
In this study, we have introduced a novel systematic and computationally efficient feature engineering framework based on 2-point spatial correlations for the quantification of cubic crystal structures with varied chemistries. The approach presented in this work employed a voxelized representation of the crystal structure and computed the spatial correlations on suitably selected atomic attributes (here taken as heat of fusion, ionization energy, and Pauling electronegativity). The proposed approach offers a novel avenue for feature engineering that will allow interpolations across different chemistries of the compounds. This work demonstrates that the proposed feature engineering approach combined with PCA for dimensionality reduction is capable of generating a compact set of salient features (i.e., PC scores) representing the many details of the crystal structure. These PC scores can be used effectively in building surrogate models needed to screen for materials exhibiting potential for improved properties. In this work, this was demonstrated by using the same features as inputs to two very different model-building strategies (i.e., L-GPR and NN) for the predictions of crystal formation energy. The results of this work indicate that the utility of the generated features is independent of the model-building strategy. Moreover, the distribution of errors in predictions (and their associated uncertainties in case of Bayesian approaches) offer objective guidance on where additional training data should be targeted to further improve the performance of the surrogate models. This work also found that the L-GPR modeling strategy produces more robust predictions when dealing with relatively smaller datasets, as the one utilized in this study.

Declaration of Competing Interest
The authors declare that no known competing financial interests have influenced the work reported in this paper.