Spline-based neural network interatomic potentials: blending classical
  and machine learning models
While machine learning (ML) interatomic potentials (IPs) are able to achieve
accuracies nearing the level of noise inherent in the first-principles data to
which they are trained, it remains to be shown if their increased complexities
are strictly necessary for constructing high-quality IPs. In this work, we
introduce a new MLIP framework which blends the simplicity of spline-based MEAM
(s-MEAM) potentials with the flexibility of a neural network (NN) architecture.
The proposed framework, which we call the spline-based neural network potential
(s-NNP), is a simplified version of the traditional NNP that can be used to
describe complex datasets in a computationally efficient manner. We demonstrate
how this framework can be used to probe the boundary between classical and ML
IPs, highlighting the benefits of key architectural changes. Furthermore, we
show that using spline filters for encoding atomic environments results in a
readily interpreted embedding layer which can be coupled with modifications to
the NN to incorporate expected physical behaviors and improve overall
interpretability. Finally, we test the flexibility of the spline filters,
observing that they can be shared across multiple chemical systems in order to
provide a convenient reference point from which to begin performing
cross-system analyses.
Joshua A. Vita, Dallas R. Trinkle
# Spline-based neural network interatomic potentials: blending classical and machine learning models

###### Abstract

While machine learning (ML) interatomic potentials (IPs) are able to achieve accuracies nearing the level of noise inherent in the first-principles data to which they are trained, it remains to be shown if their increased complexities are strictly necessary for constructing high-quality IPs. In this work, we introduce a new MLIP framework which blends the simplicity of spline-based MEAM (s-MEAM) potentials with the flexibility of a neural network (NN) architecture. The proposed framework, which we call the spline-based neural network potential (s-NNP), is a simplified version of the traditional NNP that can be used to describe complex datasets in a computationally efficient manner. We demonstrate how this framework can be used to probe the boundary between classical and ML IPs, highlighting the benefits of key architectural changes. Furthermore, we show that using spline filters for encoding atomic environments results in a readily interpreted embedding layer which can be coupled with modifications to the NN to incorporate expected physical behaviors and improve overall interpretability. Finally, we test the flexibility of the spline filters, observing that they can be shared across multiple chemical systems in order to provide a convenient reference point from which to begin performing cross-system analyses.

interatomic potential splines machine learning interpretability

## 1 Introduction and Background

In the fields of computational materials science and chemistry, machine learning interatomic potentials (MLPs) are rapidly becoming an essential tool for running high-fidelity atomic-scale simulations. Notably, major breakthroughs in the field have typically been marked by the development of new methods for encoding the atomic environments into machine-readable descriptors, or using architectures from other fields of machine learning to improve regression from descriptors into energies and atomic forces. The seminal work using atom-centered symmetry functions (ACSFs) Behler and Parrinello (2007) as the encoder, followed by "smoooth overlap of atomic position" (SOAP) descriptors Bartok et al. (2010), and, most recently, equivariant message-passing networks Batzner et al. (2022) being milestones of particular importance in the field. In parallel with these improvements to the encoding functions of MLPs has been the development of new regression tools, favored either for the simplicity gained through the use of linear combinations of basis functions (e.g., Thompson et al. (2015), Shapeev (2016), Drautz (2019)) or the accuracy gained by increasing the effective cutoff radius of the model using message-passing neural networks (e.g., Gilmer et al. (2017), Schutt et al. (2018), Batatia et al. (2022)). While countless other models have been proposed using combinations or variations of different methods (an incomplete list: Mueller et al. (2020), Manzhos and Carrington (2021), Christensen et al. (2020), Musaelima et al. (2022), Gasteiger et al. (2021), Haghighatliar et al. (2021), Lubbers et al. (2018), Hu et al. (2021), Batatia et al. (2022)), the key insights remain the same: interatomic potentials can be greatly improved by leveraging architectures and optimization strategies drawn from other machine learning and deep learning applications.

Despite the success of these MLPs, there has been a persistent need in the community for the continued development of low-cost classical IPs, particularly for use with large-scale simulations Sosso et al. (2016), Ravelo et al. (2013),
Diemand et al. (2013); Phillips et al. (2020); Zepeda-Ruiz et al. (2017). Some of the major benefits that classical IPs Jones (1924); Daw and Baskes (1984); Buckingham (1938); Tersoff (1986); Brenner et al. (2002); Shan et al. (2010); van Duin et al. (2001); Baskes (1992) have over MLIPs are their low computational costs, strong foundations in known physics, and a history of scientific research analyzing their behaviors and theoretical limitations. Although improvements are being made to the speeds of MLIPs, the stark differences between classical and ML IPs in terms of size, design, and overall complexity make it difficult to leverage the well-established tools and knowledge from classical models in order to further improve their ML counterparts.

In this work, we develop a new IP model whose hyper-parameters can be tuned to transition smoothly between low-cost, low-complexity classical models and full MLIPs. By basing the proposed model off of a classical spline-based MEAM potential, then extending it using a basic ML architecture, we enable direct comparisons to well-established classical forms as well as modern ML models. These results help to bridge the gap between the two classes of models and highlight methods for improving speeds, interpretability, and transferability of MLIPs.

### s-MEAM

The model developed in this paper builds heavily upon the spline-based MEAM potential Lenosky et al. (2000) ("s-MEAM"), which we describe here in order to provide sufficient background to understand the new model architecture proposed in this work. The s-MEAM potential is a spline-based version of the popular analytical MEAM potential Baskes (1992) that was intended to provide additional flexibility to the model while maintaining the same overall functional form. In the s-MEAM formalism, the energy \(E_{i}\) of a given atom \(i\) is written as:

\[\begin{split} E&=\sum_{i}E_{i}=\sum_{i}\left[U_{c_{ i}}(n_{i})+\sum_{j<i}\phi_{c_{j}}(r_{ij})\right]\\ n_{i}&=\sum_{j\neq i}\rho_{c_{j}}(r_{ij})\\ &\qquad+\sum_{\begin{subarray}{c}j<k,\\ j,k\neq i\end{subarray}}f_{c_{j}}(r_{ij})f_{c_{k}}(r_{ik})g_{c_{j}c_{k}}[\cos \left(\theta_{jik}\right)]\,.\end{split} \tag{1}\]

In Eqn. 1 the energy of the \(i\)-th atom, \(E_{i}\), is composed of a pair term (\(\phi\)) and an embedding energy contribution (\(U\)) for a given "electron density" \(n_{i}\), where all five functions (\(\phi\), \(U\), \(\rho\), \(f\), and \(g\)) are represented using cubic Hermite splines. The pair term is calculated by summing over pair distances \(r_{ij}=|\vec{r_{j}}-\vec{r_{i}}|\) between each atom \(i\) and its neighbors \(j\) (with \(r_{ij}\) less than a chosen cutoff distance, \(r_{c}\)). The electron density \(n_{i}\) is further decomposed into 2-body (\(\rho\)) and 3-body (products of \(f\) and \(g\)) contributions. The 2-body term is similar to the summation over \(\phi\), while the 3-body term is computed by summing over the product of three spline functions that take \(r_{ij}\), \(r_{ik}\), and \(cos(\theta_{jik})\) as inputs, where \(\theta_{jik}\) is the bond angle formed by atoms \(i\), \(j\), and \(k\) with \(i\) at the center. The subscripts on the functions indicate that separate splines are used for evaluation depending on the chemistries \(c_{i}\), \(c_{j}\), and \(c_{k}\) of atoms \(i\), \(j\), and \(k\) (e.g., \(g_{AA}\) for A-A bonds, \(g_{AB}\) for A-B bonds, etc.).

In order to facilitate comparisons between s-MEAM and the model that will be proposed in this work, we will first define two new functions

\[\begin{split} G_{3,i}^{\alpha}(r_{ij},r_{ik},\cos\theta_{jik})= \\ \sum_{\begin{subarray}{c}j<k\\ j,k\neq i\end{subarray}}f_{c_{j}}^{\alpha}(r_{ij})f_{c_{k}}^{\alpha}(r_{ ik})g_{c_{jk}}^{\alpha}\left(\cos\left(\theta_{jik}\right)\right)\\ G_{2,i}^{\beta}(r_{ij})&=\sum_{j\neq i}\rho_{c_{j}}^{ \beta}(r_{ij}),\end{split} \tag{2}\]

then re-write Eqn. 1 as

\[\begin{split} E&=\sum_{i}\left[U_{c_{i}}(n_{i})+\frac {1}{2}G_{2,i}^{0}\right]\\ n_{i}&=\sum_{\beta=1}^{N_{2}}G_{2,i}^{\beta}+\sum_{ \alpha=1}^{N_{3}}G_{3,i}^{\alpha}\end{split}\,, \tag{3}\]

where \(N_{2}=1\) and \(N_{3}=1\). We will henceforth refer to \(G_{3,i}^{\alpha}\) and \(G_{2,i}^{\beta}\) as 3-body and 2-body "spline filters" respectively, in acknowledgement of the fact that they can be thought of as filters that characterize the local environment
around atom \(i\) in order to produce a scalar atomic environment descriptor \(n_{i}\). In Eqn. 3 we have introduced summations over the superscripts \(\alpha\) and \(\beta\) which currently only take on a single value of 1, but will be used later to denote different filters.

### Nnp

Shown to be universal function approximators Hornik et al. (1989), neural networks (NNs) are provably more flexible than classical IPs which use explicit analytical forms. Because of this, a sufficiently large NN would be expected to be able to accurately reproduce an arbitrary potential energy surface, assuming that it properly accounted for long-range interactions, was provided with enough fitting data, and did not suffer from limitations due to trainability. The original Behler-Parrinello NNP Behler and Parrinello (2007) was one of the first successful applications of NNs towards practical systems, where the atomic energy of atom \(i\) is written as

\[E=\sum_{i}N_{c_{i}}(\vec{D}_{i}) \tag{4}\] \[\vec{D_{i}}=\langle D_{3,i}^{1},\ldots,D_{3,i}^{N_{3}},D_{2,i}^{1 },\ldots,D_{2,i}^{N_{2}}\rangle.\]

In Eqn. 4, \(N_{c_{i}}\) is a neural network, \(c_{i}\) is the element type of atom \(i\), and \(\vec{D}_{i}\) is the atom-centered symmetry function (ACSF) descriptor Behler and Parrinello (2007) of atom \(i\). The ACSF local environment descriptor is comprised of radial symmetry functions

\[D_{2,i}^{\beta}(r_{ij})=\sum_{j\neq i}e^{-\eta^{\beta}(r_{ij}-R_{s}^{\beta})^{ 2}}v_{c}(r_{ij}), \tag{5}\]

parameterized by \(\eta^{\beta}\) for changing the width of the Gaussian distribution, and \(R_{s}^{\beta}\) to shift the distribution. A smooth cutoff function \(v_{c}(r_{ij})\) is used with the form:

\[v_{c}(r_{ij})=\left\{\begin{array}{ll}0.5\times\big{[}\cos\frac{\pi r_{ij}} {r_{c}}+1\big{]},&\mbox{ if }r_{ij}\leq r_{c}\\ 0,&\mbox{ if }r_{ij}>r_{c}.\end{array}\right. \tag{6}\]

Angular contributions are accounted for using the angular symmetry functions

\[D_{3,i}^{\alpha}(r_{ij},r_{ik},r_{jk},\theta_{jik})= \tag{7}\] \[2^{1-\zeta^{\alpha}}\sum_{j,k\neq i}(1+\lambda^{\alpha}\cos \theta_{jik})^{\zeta^{\alpha}}\] \[\times e^{-\eta^{\alpha}(r_{ij}^{2}+r_{ik}^{2}+r_{jk}^{2})}v_{c}( r_{ij})v_{c}(r_{ik})v_{c}(r_{jk}).\]

Multiple radial and angular symmetry functions are constructed by making \(N_{2}\) choices for \(\eta^{\beta}\) and \(R_{s}^{\beta}\), and \(N_{3}\) choices for \(\zeta^{\alpha}\), \(\lambda^{\alpha}\)\((=\pm 1)\), and \(\eta^{\alpha}\). The evaluations of all of these symmetry functions are then concatenated together into a single vector \(\vec{D}_{i}\) which is passed through a feed-forward neural network.

Obvious parallels can be drawn between the NNP form in Eqn. 4 and the re-written form of s-MEAM shown in Eqn. 3. What differentiates NNP from s-MEAM, however, are the details regarding the construction of the local descriptors, and the form of the embedding function. Where s-MEAM uses trainable spline filters for both the descriptor and the embedding function, NNP uses ACSFs and an NN respectively. Although an NN would have an increased fitting capacity over the \(U_{c_{i}}\) splines used in Eqn. 3, there are many similarities between the ACSF descriptors and the spline filters described in Eqn. 2. For example, the 2-body components of an ACSF descriptor, which are constructed by evaluating \(D_{2,i}^{\beta}\) with multiple radial shifts \(R_{s}^{\beta}\) for each neighboring atom \(j\), can be viewed as basis functions used for interpolating the desired range of atomic distances. This is conceptually related to how the basis functions of cubic Hermite splines allow \(G_{2,i}^{\beta}\) to interpolate over its domain as well. A similar argument can be made relating the angular components of ACSFs to the 3-body filters \(G_{3,i}^{\alpha}\), where Eqn. 2 and Eqn. 5 both multiply functions of pair distances (\(f_{c_{k}}^{\alpha}(r_{ij})\) and \(f_{c_{k}}^{\alpha}(r_{ik})\)) by a function of the triplet angle (\(g_{c_{jk}}^{\alpha}(\cos{(\theta_{jik})})\)). The ANI model Smith et al. (2017), which will be used in this work for comparison to the model which we developed, is nearly identical to the NNP form described above, with the modifications that only a single \(\eta^{\beta}\) is used, and Eqn. 7 is altered to introduce both radial (\(R_{s}^{\alpha}\)) and angular (\(\theta_{s}^{\alpha}\)) shifts:

\[D_{3,i}^{\alpha}(r_{ij},r_{ik},\theta_{jik})= \tag{8}\] \[2^{1-\zeta^{\alpha}}\sum_{j,k\neq i}[1+\cos(\theta_{jik}-\theta _{s}^{\alpha})]^{\zeta^{\alpha}}\] \[\times e^{-\eta^{\alpha}(\frac{r_{ij}+r_{ik}}{2}-R_{s}^{\alpha})} v_{c}(r_{ij})v_{c}(r_{ik}).\]
## 2 Methods

### s-NNP

The main contribution of this work is the development of a spline-based neural network potential, outlined in Fig. 1, which we call the "s-NNP" (short for "spline-NNP"). The s-NNP framework is intended to extend the fitting capacity of the s-MEAM class of potentials while maintaining the high interpretability and speed provided by the use of splines. s-NNP can be thought of as an s-MEAM potential with two critical changes: first, \(N_{3}\) and \(N_{2}\) in Eqn. 3, the numbers of \(G^{\alpha}_{3,i}\) and \(G^{\beta}_{2,i}\) spline filters, are taken to be hyper-parameters of the model; and second, the "embedding" spline \(U_{i}\) in an s-MEAM model is replaced by a fully-connected neural network. By including additional spline filters, the s-NNP is able to describe the local environment around an atom with increasingly fine resolution (since each spline can be thought of as a basis function for interpolating the atomic energies). The introduction of a neural network allows the model to achieve much more complex mappings into atomic energies than would be possible with the cubic splines \(U_{c_{i}}\). In the s-NNP formalism, the energy \(E_{i}\) of a given atom \(i\) is then written as

\[\begin{split} E&=\sum_{i}N(\vec{G}_{i}),\\ \vec{G}_{i}&=\langle G^{1}_{3,i},\dots,G^{N_{3}}_{3,i},G^{1}_{2,i},\dots,G^{N_{2}}_{2,i}\rangle.\end{split} \tag{9}\]

where \(N\) is a neural network, and \(\vec{G}_{i}\) is a vector of length \((N_{3}+N_{2})\). Notice that each component of \(\vec{G}_{i}\) is computed by evaluating a different 3- or 2-body spline filter for the local environment of atom \(i\). The parameters of an s-NNP that are trained during fitting are the positions of the knots of the \(f^{\alpha}\), \(g^{\alpha}\), and \(\rho^{\,\beta}\) splines from Eqn. 2, as well as the weights and biases in the neural network \(N\). The hyper-parameters of the model are \(N_{3}\), \(N_{2}\), the number of knots in each spline, the number of layers in \(N\), and the number of hidden nodes in each of those layers.
One benefit of the s-NNP framework is that it is closely related to both classical and machine learning interatomic potentials (see Section 4.2 for more discussion), making it possible to easily probe the performance gap between the two. For example, many classical potentials (LJ Jones [1924], EAM Daw and Baskes [1984], MEAM Baskes [1992], Stillinger-Weber Stillinger and Weber [1985], Tersoff Tersoff [1986], and Buckingham Buckingham [1938]) could be reformulated as s-NNP potentials with very few filters (\(N_{3}\in[0,1]\), \(N_{2}\in[1,2]\)) and custom embedding functions instead of a neural network (though given the universal approximation theorem, these embedding functions could be represented using an NN as well). Because of this, we can easily construct spline-based "classical" models by adjusting \(N_{3}\) and \(N_{2}\), but not including a neural network, then compare them directly to MLIPs by subsequently attaching networks with varying depths and widths. See Section 3 in the Results for details of such a study.

### Interpretability improvements

A central tenet of constructing interpretable models is designing the model architecture in a way that helps isolate the contributions of specific parameter sets to the final model predictions. With this in mind, we therefore propose the five modifications described in Fig. 2. Although the modifications proposed here are only discussed in the context of s-NNP, they can also be applied to many other existing MLIP frameworks.

When applying all of the modifications described in Fig. 2, the full form of s-NNP is written as:

\[\begin{split} E=\sum_{i}E_{i}&=\sum_{i}\left[E_{ \text{lin},i}+E_{\text{net},i}\right]\\ E_{\text{lin},i}&=E_{\text{se},i}+G_{2,i}^{\text{ short}}+b\\ E_{\text{sc},i}&=\sum_{n=1}^{N_{3}}G_{3,i}^{n}+ \sum_{n=1}^{N_{2}}G_{2,i}^{n}\\ E_{\text{net},i}&=\sigma\big{[}N_{\text{no-bias}}( \vec{G}_{i})\big{]},\end{split} \tag{10}\]

where \(E_{\text{sc},i}\) as a "skip connection" term (as used in other DL fields He et al. [2015]), \(G_{2,i}^{\text{short}}\) is a short-range 2-body filter, \(b\) is a trainable isolated atom energy, \(\sigma\) is the Softplus activation function \(\sigma(x)=\log{(1+e^{x})}\), and \(N_{\text{no-bias}}\) is a

Figure 2: Modifications to the s-NNP form described in Fig. 1 for improving interpretability. The four modifications are: 1) adding skip connections; 2) removing the internal network bias, and adding an external bias, \(b\); 3) wrapping the network outputs in a Softplus activation function; and 4) introducing a single short-range pair term using a 2-body filter. The short-range pair term is only allowed to be non-zero up to a cutoff of 2.5 Å.
neural network with no bias terms on any of its layers. Note that \(E_{\text{lin},i}\) encompasses all of the terms which are linearly dependent upon the spline filters, and \(E_{\text{net},i}\) captures the non-linear dependence. While a more in-depth discussion of the effects of each of these modifications can be found in Section 4.1, the practical result of Eqn. 10 is that the spline filter visualizations like those shown in Fig. 3 can be intuitively understood. For example, negative regions of the filters will usually correspond to negative \(E_{i}\), and positive regions of the filters will always correspond to positive \(E_{i}\). While Eqn. 10 still suffers from some drawbacks, predominantly associated with the non-linear effects of \(N_{\text{no-bias}}(\widetilde{G}_{i})\), we believe that it strikes a good balance between the accuracy achieved through a NN architecture, and the interpretability characteristic of classical potentials.

### Benchmarking dataset: AI

In order to test our framework, in Section 3.1 we fit s-NNP models to the aluminum dataset from Smith et al. [2021], which serves as a good benchmarking dataset due to its size and configurational complexity. This dataset was built using an active learning technique for generating interatomic potential training data Smith et al. [2018]. It has also been shown to contain extremely diverse atomic environments, such that the ANI model Smith et al. [2019] which was originally trained to the dataset was tested for use in shock simulations. After removing duplicate calculations (identical atomic positions and computed energies/forces, generated at different active learning steps), and one outlier configuration with particularly high forces, the final dataset consisted of 5,751 unique configurations (744,356 atoms total) with their corresponding DFT-computed energies and forces. A 90:10 train-test split was performed to partition the dataset into training/testing data. This data can be obtained from the original source ani.

### Additional datasets: Cu, Ge, Mo

As a test of the flexibility of the spline filters, in Section 3.3 we trained an s-NNP model simultaneously to the Al dataset described in Section 2.3 and the Cu, Ge, and Mo datasets from Zuo et al. [2020]. Each dataset from Zuo et al. [2020] was manually constructed and contains the ground state structure for the given element, strained supercells, slab structures, _ab initio_ molecular dynamics (AIMD) sampling of supercells at different temperatures (300 K and \(0.5\times\), \(0.9\times\), \(1.5\times\), and \(2.0\times\) the melting point), and AIMD sampling of supercells with single vacancies at 300 K and \(2.0\times\) the melting point. On average, each training dataset includes approximately 250 structures for a total of 27,416 atoms (Cu), 14,072 atoms (Ge), and 10,087 atoms (Mo). A detailed summary of the contents of the datasets can be found in Zuo et al. [2020]. This data can be obtained from the original source mle.

In order to attempt to balance the combined dataset, we took a random subset comprised of 20% of the Al dataset in addition to the full Cu, Ge, and Mo datasets, resulting in a total training set size of 1271 Al configurations, 262 Cu, 228 Ge, and 194 Mo. Another logical choice for constructing the multi-element training set would be to further downsample the Al dataset to better balance the relative concentrations of each element type. However, we observed that doing so resulted in worse property predictions for Al, presumably due to the large number of high-energy configurations in the original Al dataset making random sampling of low-energy configurations relevant to the properties of interest unlikely. In all cases, the atomic energies of each element type were shifted by the average energy of that type before training.

## 3 Results

### Benchmarking tests

Using the s-NNP framework, we fit a collection of models to the Al dataset to probe the effects of increasing model capacity in two distinct ways: first by increasing the number of spline filters, and second by increasing the size of the network used for mapping filter outputs to energies. Table 1 outlines the architectures and accuracies of the trained models, grouped by key architectural changes and sorted by total number of parameters. The trained s-NNP models can be conceptually broken down into three main groups, each highlighting the effects of specific architectural changes: 1) "\((N_{2},N_{3})\) linear" models, increasing the number of splines when no neural network is used; 2) "\((N_{2},N_{3})\)\(l\) =*" models, increasing the number of splines and network size; and 3) introducing the interpretability improvements described in Section 2.2.

The results in Table 1 show multiple avenues for systematically improving the performance of an s-NNP model, though each method appears to experience a saturation point beyond which the model suffers from diminishing returns as complexity increases. Using the "\((1,1)\) linear" model as a baseline, we see that increasing \(N_{3}\) from 1 to 8 can monotonically decrease both the energy and the force errors. Increasing \(N_{2}\) has no significant effect on the accuracy
of the linear models, which is consistent with the notion that a linear combination of cubic Hermite splines can be represented using a single spline.

The introduction of a neural network enables the model to better utilize the additional spline filters, leading to significant improvement of the "\((1,2)\)\(l=3\)" model over any of the linear models. Increasing the number of spline filters, and subsequently increasing the network width and depth to maintain a "funnel-like" structure (decreasing layer width with increasing depth) can also lead to steady improvements. However, with increasing model size we began to see a commensurate increase in training difficulty, often resulting in larger models with higher errors than what might be expected based on the performance of their smaller counterparts (e.g., "(2,4) \(l=6\)" compared to "(2,4) \(l=5\)", or "(1,8) \(l=5\), int. wide" compared to "(1,8) \(l=5\), int. all"). Notably, the interpretability improvements from Section 2.2 did not hurt model performance, demonstrating that accuracy and interpretability are not strictly opposing attributes of a model. Based on the results shown in Table 1, we will use the "(1,8) \(l=5\), int. all" for all experiments and analyses in the remainder of this paper.

### Model visualization

A major advantage of s-NNP over many other MLIPs, especially when coupled with the modifications from Section 2.2, is that the spline filters \(G^{\alpha}_{3,i}\) and \(G^{\beta}_{2,i}\) lend themselves to easy visualization. This can be valuable for helping model developers and users to better understand how their model is interacting with the data or influencing simulation results. The polar plots in Fig. 2(a), corresponding to the "(1,8) \(l=5\), int. all" model, represent the total activation of the 3-body filters induced by placing an atom at a given \((r,\theta)\). These visualizations make it easy to recognize

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline Model name & \((N_{2},N_{3})\) & Network & Parameters & \(E_{\text{RMSE}}\) & \(F_{\text{RMSE}}\) \\  & & & (splines, network) & (meV/atom) & & (eV/Å) \\  & & & Train & Test & Train & Test \\ \hline \((1,1)\) linear & \((1,1)\) & - & (66, 0) & 74 & 78 & 0.76 & 0.82 \\ \((1,2)\) linear & \((1,2)\) & - & (110, 0) & 56 & 56 & 0.54 & 0.58 \\ \((1,4)\) linear & \((1,4)\) & - & (198, 0) & 40 & 41 & 0.46 & 0.49 \\ \((2,4)\) linear & \((2,4)\) & - & (220, 0) & 38 & 36 & 0.46 & 0.50 \\ \((1,8)\) linear & \((1,8)\) & - & (374, 0) & 32 & 34 & 0.42 & 0.46 \\ \((1,2)\)\(l=3\) & \((1,2)\) & \((3,2,1)\) & (110, 23) & 44 & 44 & 0.33 & 0.35 \\ \((1,4)\)\(l=5\) & \((1,4)\) & \((5,4,3,2,1)\) & (198, 80) & 10 & 11 & 0.22 & 0.23 \\ \((2,4)\)\(l=4\) & \((2,4)\) & \((6,3,2,1)\) & (220, 74) & 22 & 21 & 0.21 & 0.21 \\ \((2,4)\)\(l=5\) & \((2,4)\) & \((6,4,3,2,1)\) & (220, 96) & 8.7 & 9.0 & 0.19 & 0.20 \\ \((2,4)\)\(l=6\) & \((2,4)\) & \((6,5,4,3,2,1)\) & (220, 127) & 10 & 12 & 0.24 & 0.25 \\ \((1,8)\)\(l=5\) & \((1,8)\) & \((9,8,4,2,1)\) & (374, 219) & 7.5 & 6.5 & 0.13 & 0.13 \\ \((1,8)\)\(l=5\), int. skip & \((1,8)\) & \((9,8,4,2,1)\) & (374, 219) & 5.5 & 4.6 & 0.15 & 0.15 \\ \((1,8)\)\(l=5\), int. all & \((1,8)\) & \((9,8,4,2,1)\) & (374, 219) & 5.5 & 5.9 & 0.12 & 0.12 \\ \((1,8)\)\(l=5\), int. wide & \((1,8)\) & \((128,64,32,16,1)\) & (374, 11,793) & 5.9 & 5.1 & 0.13 & 0.14 \\ ANI Smith et al. (2021) & - & \((96,96,64,1)\) & (-, 15,585) & - & 1.9 & - & 0.06 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Fitting results comparing linear s-NNP models, “\((N_{2},N_{3})\) linear”, s-NNP models with neural networks, “\((N_{2},N_{3})\)\(l=\)*”, s-NNP with interpretability improvements “\((1,8)\)\(l=5\), int. *”, and the ANI model from Smith et al. (2021), “ANI”. The \((N_{2},N_{3})\) notation for the spline layers indicates that \(N_{2}\) 2-body and \(N_{3}\) 3-body spline filters were used. Network architectures are denoted as a tuple of integers (“Network” column) specifying the number of nodes in each hidden layer of the model. Although the training/testing errors decrease significantly when adding additional splines to the linear models, their performance appears to begin to saturate with the “(1,8) linear” model. Model performance only begins to be competitive with the ANI results upon the inclusion of a neural network with sufficient depth, which allows for non-linear combinations of spline outputs. While the “(1,8) \(l=5\), int. skip” model only adds the use of skip connections, the “(1,8) \(l=5\), int. all” and “(1,8) \(l=5\), int. wide” use all of the interpretability improvements discussed in Section 2.2 and Section 4.1. Note that this means that the single 2-body filter specified using the “(1,8)” notation is the short-range pair term described in Fig. 2, and that there are therefore no 2-body filters being passed through the network. Testing errors for the “ANI” model were taken directly from Smith et al. (2021), which did not report training errors. Note that the ANI model uses ensemble-averaging over 8 networks, which they report yields energy and force errors that are “20% and 40% smaller, respectively, compared to a single ANI model”.
aspects of the local environments around an atom that are learned during training to have lower energy. For example, the averaged filter shown in Fig. 3b has an attractive behavior of the model for bond angles \(60^{\circ}\leq\theta_{jik}\leq 90^{\circ}\) in addition to repulsion for angles larger than approximately \(120^{\circ}\). Many of the individual filters in Fig. 3a also show characteristic features at various bond angles, especially for bond lengths near the first and second nearest neighbor distances in FCC aluminum. Fig. 3c shows the short-range pair term, which was learned to have a strongly repulsive contribution, as would be expected based on the Pauli exclusion principle.

It's worth mentioning that similar plots to the ones shown in Fig. 3 could also be generated for other NNP-like models, for example by visualizing each of the components of the vector output of the first hidden layer in an ANI model. However, most other neural network-based architectures would have significantly more filters to visualize given the size of their hidden layers. For example, the ANI model in Table 1 has 96 nodes in its first hidden layer, thus making it more difficult to interpret the results. Furthermore, since ANI (and most other models) does not use skip connections summing the hidden layer directly into the output, any resultant visualization will not necessarily be in units of energy, meaning it may undergo significant non-linear transformations as it passes through the network.

### Flexibility tests

The high interpretablility of the s-NNP spline filters, as discussed in Section 3.2, is particularly valuable when the filters can be applied to different chemical systems in order to enable cross-system comparisons. For example, if multiple s-NNP models using the same spline filters were trained to different chemical systems, then the networks of each model could be analyzed in order to understand how differences in chemistries result in different sensitivities to the local environment embeddings defined by the spline filters. The multi-element model trained to the datasets described in Section 2.4 uses the same spline filters for all of the data, but separate NNs for each element type. The network architecture maintained a funnel-like structure of \((12,8,4,2,1)\). It also used four additional 3-body filters in order to increase its fitting capacity, for a total of one short-range pair filter and 12 3-body filters, though without

Figure 3: Visualizations of spline filters of the “(1,8) final” model from Table 1 for atomic distances in the range [2.5Å, 7Å]. a) Plots of the individual \(G^{\alpha}_{3,i}\) filters, where integer labels above each plot indicate the index \(\alpha\). b) The average of all \(G^{\alpha}_{3,i}\) 3-body filters. c) The short-range pair term, as described in Section 2.2. Note that the radial splines \(f^{\alpha}_{c_{j}}\), \(f^{\alpha}_{c_{k}}\), and \(\rho^{\alpha}_{c_{j}}\) use linear extrapolation for distances outside of the domain defined by their knots. Dashed grey lines in the polar plots correspond to the first and second nearest-neighbor distances (2.86 Å and 4.05 Å respectively) for FCC Al at room temperature Nakashima [2020]. Each point in the polar plots is computed by fixing atom \(i\) at the origin, fixing atom \(j\) at the given \((r,\theta)\), fixing atom \(k\) along the \(\theta=0\) axis, then integrating \(G^{\alpha}_{3}\) over \(r_{k}\in[1.5,7.0]\). Note that there is a forced symmetry in the polar plots since the \(g^{\alpha}_{c_{jk}}\) splines in Eqn. 2 take \(\cos{(\theta)}\) as input. All polar plots use the same color scale, where values are clipped to fall within a chosen range to optimize for visibility while avoiding information loss. Values for \(r<1.5\) Å, which was the smallest distance sampled in the Al dataset as shown in Fig. B1, are omitted to ensure that high signals at small atomic distances would not wash out the rest of the colors in the plots. Though all plots in this figure are technically in units of eV, this would not be true if skip connections were not used.
the use of skip connections or the Softplus activation function. In the case of the multi-element model, the use of skip connections would imply a "background" energy that is consistent across all four element types and is augmented by the network contributions for each element. While this assumption sounds plausible in theory, we found that in practice a model using skip connections resulted in poorer property predictions than one which did not. This decreased performance when using skip connections can likely be attributed to the larger concentration of Al data dominating the fitting and causing the filters to learn a background energy that is influenced by the high-energy Al configurations and therefore not transferable to the lower energy Cu/Ge/Mo data. We believe that further research into constructing balanced training datasets could help address this issue, though it may also be possible that a larger number of splines is required for describing this more diverse dataset. Another alternative, which was not explored here, would be to only use skip connections for a subset of the spline filters in order to give the model the ability to learn a transferable background energy without enforcing the full constraint of skip connections for all filters.

The property predictions of the multi-element model plotted in Fig. 4 show that a model using shared filters for multiple datasets can learn to make reasonable property predictions for all four elements studied in this work. While the initial predictions for most of the elements were relatively good, the cubic elastic constants and bulk modulus for Mo were noticeably under-predicted, and the vacancy migration energy was far too large to be considered acceptable (dotted line in Fig. 4). Due to the relatively few number of spline filters used, and their high degree of interpretablity, we were able to fine-tune the model by zeroing out the network weights of hand-chosen filters in order to remove their contribution to the Mo energy predictions. For example, the contributions of each spline filter can be visualized

Figure 4: Property predictions of the lattice constant (\(a\)), cubic elastic constants (\(C_{ij}\)), bulk modulus (\(B\)), and vacancy formation and migration energies (\(E_{vf}\) and \(E_{m}\)) for s-NNP models compared to existing MLIPs. s-NNP results are from this work, ANI results from Smith et al. (2021), and all others from Zuo et al. (2020). All properties have been normalized with respect to the DFT-predicted values. Note that Smith et al. (2021) did not compute \(E_{m}\) for ANI. The “(1,8) \(l=5\), int. all” model (solid black line), which was trained only to Al, predicts all seven properties well, with the exception of \(C_{44}\), which has a relatively small magnitude compared to the other elastic constants. The performance of the multi-element s-NNP model (dashed black lines) is within the range of the other single-element MLIPs, though s-NNP’s predictions for Ge are somewhat distorted. The “fine-tuned” s-NNP model was manually adjusted after training was completed in order improve the property predictions, as described in Section 3.3.
individually in order to isolate the influence of each filter on the properties of interest (see Fig. C2). Following this approach, we were able to improve the Mo predictions to bring them within an acceptable range without altering the predictions for the other three elements (dashed lines in Fig. 4). However, we note that the original training errors of the model (before fine-tuning) as shown in Fig. C4 were comparable to those of the NNIPs from Zuo et al. (2020), suggesting that the property predictions of the multi-element s-NNP model may have been able to have been improved by re-balancing the dataset. Similar to what was described in Section 2.4, a possible explanation for this is that the Mo dataset did not fully constrain the portions of the spline filters relevant to computing the properties of interest, and their shapes were therefore dictated by the Al dataset, leading to corrupted Mo property predictions. Further research into methods for preventing this type of "cross-pollution" of information would be valuable for constructing more flexible and generalizable models.

In order to gain additional insights into the differences in energetics between the Al, Cu, Ge, and Mo systems, we compute the sensitivities of the energies predicted by the multi-element model with respect to each of the 3-body spline filters. The sensitivity of the total energy, \(E\), for all configurations \(N\) with respect to a given 3-body spline filter \(G_{3}^{\alpha}\) can be computed as:

\[s_{E}^{\alpha}=\sum_{i=1}^{N}|\partial E_{\text{s-NNP}}/\partial G_{3,i}^{ \alpha}|/E_{\text{DFT}}, \tag{11}\]

which can be easily computed via back-propagation.

Calculation and comparison of these sensitivities for each dataset, as shown in Fig. 5, highlight the learned similarities and differences between the four elemental datasets. Examination of the correlation between the sensitivities (right panel of Fig. 5) shows that the Cu and Mo datasets have similar filter sensitivities, while Ge is the only element to have a negative correlation with any of the other elements. The Al dataset appears to be uniformly similar to all three remaining elements, reflecting the fact that most of the geometric configurations of the Cu/Ge/Mo datasets are well-sampled by the Al dataset (see Fig. C3). We hypothesize that the negative correlation of the Ge sensitivities is due to an attempt by the model to encode the large cluster of outlying Ge points in the UMAP visualizations shown in Fig. C3.

## 4 Discussion

### Understanding model behavior

The main purpose behind the modifications proposed in Fig. 2 and Section 2.2 is to simplify the process of understanding the behavior of s-NNP models. This not only improves the usefulness of visualizations like those shown in Fig. 3 and Fig. C1, but can also aid in debugging model predictions in practice. In this section, we will discuss each modification from Section 2.2 in greater detail in order to highlight how they improve the interpretability of the model.

Figure 5: Analysis of sensitivities \(s_{E}^{\alpha}\) (Eqn. 11) of the predicted energies with respect to the 3-body spline filters, \(G_{3}^{\alpha}\), for each elemental model. We compute the sensitivity by calculating the derivative of the predicted energy with respect to the spline filter outputs, normalized by the DFT reference value, then summing over all atoms in the training set. Note that these sensitivities are computed with respect to the total energy \(E\), which includes the network contributions. The left panel plots the raw values, while the right panel shows the correlation between the sensitivities of each element. The empty bars for Mo correspond to the filters which were removed from the Mo energy predictions during fine-tuning as described in Section 3.3. The corresponding spline filters are visualized in Fig. C1.
The first of these modifications is the use of skip connections, which result in strong theoretical changes that can not only lead to improved trainability as seen in Table 1 and other deep learning applications Li et al. (2017), but also greatly increase the interpretability of the model. Modifying \(\vec{G}_{i}\) to have units of energy seemingly contrasts with the notion of the atomic "density", \(n_{i}\), from Eqn. 1 and the local environment descriptor, \(\vec{D}_{i}\), from Eqn. 4, both of which are considered to be intermediate representations which will only become energies once they have been transformed by a regression function (\(U_{c_{i}}\) or \(N_{c_{i}}\) respectively). However, we emphasize that these three quantities (\(n_{i}\), \(\vec{D}_{i}\), and \(\vec{G}_{i}\)) remain closely related even when \(\vec{G}_{i}\) is in units of energy while the others are not. Although local environment descriptors are usually seen as quantifying the geometry within a local neighborhood, there is no reason that the environment descriptor may not itself also be an energy. An energy-based descriptor could then be thought of as a type of energy partitioning scheme, where in the case of s-NNP with skip connections the atomic energies contribute both linearly to the total energy (through the skip connection) and non-linearly (through the network). In fact, in our previous work Vita and Trinkle (2021) we observed that \(U_{c_{i}}\) was often learned to be a nearly linear function, meaning that it essentially served the purpose of a simple unit conversion for \(n_{i}\). This linearity is essential to model interpretability, and is the main motivator for the use of skip connections in s-NNP, as it is a step towards simplifying the process of understanding how the spline filters contribute to the total energy.

Nevertheless, care must still be taken when interpreting the filter visualizations when the other modifications discussed in Section 2.2 are not also employed. For example, skip connections alone do _not_ guarantee that a negative filter value means \(E_{i}\) will also be negative, since the network contribution \(N(\vec{G}_{i})\) may outweigh the skip connection term. Similarly, adjusting the filter outputs (e.g., by tuning the knot positions) may have unexpected results due to the non-linear behavior of the network. These issues are also present with all non energy-based descriptors like \(n_{i}\) and \(\vec{D}_{i}\), but can be further improved upon using the remaining techniques discussed in this section.

In order to further facilitate straightforward interpretation of spline visualizations, we also choose to pin the knots of the radial splines to be 0 at the cutoff distance \(r_{c}\), and we remove bias terms from the layers of the NN. As a result of this, both \(N(\vec{G}_{i})\) and \(E_{\text{sc},i}\) will smoothly decay to 0 at the cutoff distance, thus incorporating the desired behavior which is enforced in the Behler-Parrinello NNP using the cutoff function \(v_{c}\) in Eqn. 6. In order to account for the fact that some datasets may have non-zero energy at the cutoff distance, an external bias term, \(b\) is added so that the atomic energy converges to \(b\) at \(r_{c}\). In this sense, \(b\) corresponds to the learned energy of the isolated atom. We note that this external bias term is preferable to using internal network bias or un-pinning the radial splines because it ensures that \(N(\vec{G}_{i})\) and \(E_{\text{sc},i}\) behave similarly as \(r_{ij}\) approaches \(r_{c}\).

Wrapping the network output \(N_{\text{no-bias}}(\vec{G}_{i})\) in the Softplus activation function guarantees that the NN contributions to the total energy are strictly non-negative. This ensures that any attractive behavior of the model arises solely from the spline filters through the skip connection term, \(E_{\text{sc},i}\), thus helping to isolate certain model behaviors to specific parameter sets. It is important to note, however, that the Softplus activation should only be used in conjunction with skip connections in order to ensure that the model can predict negative energies.

A common challenge for many MLPs is ensuring a repulsive behavior at small values of \(r_{ij}\). This difficulty is not reflected in classical potentials, however, as classical models often include explicit repulsive terms. While a data-driven solution to this problem is possible, by explicitly adding short-range dimers into the training dataset, many MLIP developers choose instead to adjust the form of their IP. One such method used in the literature is to augment the model by introducing an auxiliary potential designed to capture the repulsive behavior of a pair potential. For example, in Miladovich et al. (2023) a "composite potential" was constructed by first fitting a repulsive "auxiliary potential" to DFT data, then fitting a "main potential" to the residuals of the pair potential using a Gaussian Approximation Potential Bartok et al. (2010). We incorporate a similar idea by including a short-range 2-body filter, \(G_{2,i}\), which is summed directly into the total energy without ever passing through the NN. Note that this is different from the skip connections, which are summed directly into the energy _and_ passed through the NN. Having a spline filter which only contributes linearly to the total energy means that it is free of any of the interpretability issues described for the spline filters which are used as input to the NN. While this short-range pair term does not necessarily guarantee a repulsive behavior at small distances, we observe empirically that it is often learned during training to have a strongly repulsive shape. A similar technique is used by the SNAP interatomic potential Wood and Thompson (2018), where the repulsive ZBL pair potential is added to the potential form as a "reference potential".

### s-NNP's relation to other models

The s-NNP architecture shown in Fig. 1 can be further understood by drawing relationships between itself and other models from the literature, particularly the UF3 model Xie et al. (2021) and the original Behler-Parrinello NNP Behler and Parrinello (2007). s-NNP can be compared to UF3 and NNP by analyzing the differences between the two key
components of each model: the embedding function for encoding local atomic environments into a descriptor, and the regression function for mapping the descriptor into an energy. Although it can be difficult to clearly distinguish between the embedding/regression portions of most MLIPs due to the inability to definitively establish the roles of all parameters in a deep model, we will attempt to break each model down in intuitive ways to facilitate comparison. One can view the vector \(\vec{G}_{i}\) from Eqn. 9 as a descriptor generated by an embedding function defined by the spline filters \(G_{3,i}^{\alpha}\) and \(G_{2,i}^{\beta}\). This embedding technique is most similar to the UF3 method, which also decomposes the energy into two-body and three-body terms described by spline basis functions. Though there are some differences in the exact details of the UF3 and s-NNP spline functions, for example UF3's use of tricubic B-splines instead of the 1D cubic Hermite splines of s-NNP, the general principle is the same. While s-NNP's embedding function is most similar to that of UF3, its regression function is identical to that of the Behler-Parrinello NNP Behler and Parrinello (2007). Therefore, in order to help analyze the performance of s-NNP with respect to other models in the literature, it is suitable to think of s-NNP as a combination of a UF3-like embedding function with an NNP regressor. Or, equivalently, as an NNP using a spline-based embedding function instead of atom-centered symmetry functions Behler and Parrinello (2007). However, neither UF3 nor NNP utilize all of the interpretability improvements discussed in Section 2.2.

The recently-proposed EAM-R model Nitol et al. (2023) is the most closely related model to s-NNP in the literature, as it also incorporates components of both a classical model (EAM) and an MLIP (an NNP). EAM-R is a composite potential (using the terminology described in Section 4.1) where the auxiliary potential is an EAM model, and the main potential is an RANN (an NNP with descriptors inspired by the analytical MEAM equations Baskes (1992)). Similar to this work, the developers of EAM-R observed that combining a classical model with an MLIP resulted in both improved stability relative to an MLIP and improved accuracy to a classical potential alone. Despite s-NNP's similarity to EAM-R, s-NNP has some key differences, namely the use of spline filters (instead of analytical MEAM-inspired descriptors) and some of the interpretability improvements discussed in Section 2.2. In particular, the spline filters may be expected to be more flexible than the RANN descriptors (similar to how s-MEAM is more flexible than analytical MEAM) for a given computational cost, and benefit from the ability to enforce smoothness and convergence through curvature penalties and knot pinning. Furthermore, s-NNP's use of skip connections and the removal of the internal network bias greatly improve the interpretability of the model, as discussed in Section 2.2. Although EAM-R does not include these modifications, it is an excellent example of an MLIP which could easily incorporate these same interpretability improvements.

### Computational costs

The computational cost of s-NNP is dominated by the evaluation of \(\vec{G}_{i}\), and is particularly dependant upon the choice of \(N_{3}\). In fact, basic profiling tests revealed that the filter evaluations accounted for approximately 95% of the total CPU and GPU time. This behavior can be understood heuristically by the fact that Eqn. 2 involves approximately \(\mathcal{O}(N_{n}^{2})\) spline evaluations for each filter where \(N_{n}\) is the average number of neighbors within the cutoff distance (due to the summation over triplets of atoms), as opposed to the relatively few matrix multiplications associated with the evaluation of the neural network. In general, the computational cost of inference with an s-NNP potential will scale sub-linearly with \(N_{3}\) (some speedup can be achieved by performing batched spline evaluations for the filters in Eqn. 2). An important practical implication of this is that in order to improve the accuracy of a given s-NNP model (and other NNP-based MLIPs as well), it is much more computationally efficient to increase the size of the network rather than the size of the embedding function. On the other hand, increasing \(N_{3}\) may lead to larger increases in accuracy (up to a point) than what is achievable by only increasing the network size

Given the performances of the s-NNP models in Table 1 and the timing comparisons observed in our previous work between s-MEAM and a NNP Vita and Trinkle (2021), it may be expected that an s-NNP could be constructed that achieves identical errors to ANI while maintaining a higher speed. The "(1,8) \(l=5\), int. all" model is already nearing this threshold, especially taking into account that the values for ANI reported in Table 1 use ensemble-averages over 8 models, as reported by the original authors Smith et al. (2021), which they say can make the energy and force errors "20% and 40% smaller, respectively" and may account for the differences in performance as compared to "(1,8) \(l=5\), int. all".

## 5 Conclusion

In this work we developed a novel framework using spline-based filters coupled with a neural network regressor in order to blend the strengths of both classical and ML IPs. We use this framework to probe the gap between these two classes of models, observing performance limits of linear ("classical") models that can be overcome using even small neural networks. We then show that this improved performance can be maintained while incorporating architectural changes which improve the interpretability of the model, such as the use of skip connections, an external bias term,