Interpretable Machine Learning for
Materials Design

Crystal structures for the perovskite systems were obtained from the “Stable Inorganic Perovskites” dataset published by Körbel, Marques, and Botti ⁶³, as hosted by NOMAD ¹⁹. This dataset contains a total of 144 DFT-relaxed inorganic perovskites identified via a high-throughput screening strategy. Using this dataset, we develop a model of perovskite volume. As we rely on the use of compositional descriptors for these systems, we have scaled the volume of the perovskite unit cell by the number of formula units, such that the volume has units of ų / formula unit.

Structures for 2D materials were obtained from the 2DMatPedia ¹¹, a large database containing a mixture of 6,351 real and hypothetical 2D systems. This database was generated via a DFT-based high-throughput screening approach, which investigated bulk structures hosted by the Materials Project database ¹⁸ to find systems which may plausibly form 2D structures. Among other things, the 2DMatPedia provides DFT-calculated exfoliation energies and bandgaps, along with a DFT-optimized structure for each material. We use this dataset to develop models for the bandgap and exfoliation energy of 2D materials. Although the dataset reports exfoliation energies in units of eV, we have converted these into units of J / m². Bandgaps are reported in units of eV.

To facilitate the development of ML algorithms capable of rapidly predicting material properties, we focus primarily on features that do not require further (computationally-intensive) DFT calculations. In the case of the 2D material bandgap, we include the DFT-calculated bandgap of the respective bulk material; we note that these values are tabulated on the Materials Project and can be looked up, thus circumventing the need for further DFT work. A variety of chemical featurization libraries were used to generate compositional and structural descriptors for the systems we investigated, and are listed in sections 2.2.1 and 2.2.2 respectively. Features with values of NaN (which occurred when a feature could not be calculated) were assigned a value of 0.

The choice of data filtering methodology was chosen based on the problem at-hand. The perovskite volume prediction problem did not utilize any data filtering. In the case of the 2D material bandgap and exfoliation energy prediction problems, the data obtained from the 2DMatPedia were required to satisfy all of the following criteria:

1. 1.

   No elements from the f-block, larger than U, or noble gases were allowed.

2. 2.

   Decomposition energy must be below 0.5 eV/atom.

3. 3.

   Exfoliation energy must be strictly positive.

Additionally, in the case of the 2D material bandgap, data were required to have a parent material defined on the Materials Project. This was done because we use the Materials Project’s tabulated DFT bandgap of the bulk system as a descriptor for the bandgap of the corresponding 2D system.

When training models, 10% of randomly selected data was held out as a testing set. The same train/test split was used for all 4 models considered (XGBoost, TPOT, Roost, and SISSO). To facilitate a transparent comparison between models, in all cases we report the Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), Maximum Error, and R² score of the test set.

XGBoost, TPOT, and SISSO were applied to investigate the volume of perovskites as a function of the compositional features described in section 2.2.1. Additionally, we trained a Roost model on the chemical formula of the perovskites to predict the volume. The train/test split resulted in a total of 129 entries in the training set, and 15 in the test set. We find generally good performance on the perovskite volume problem across all 4 models, although the TPOT and SISSO model display the best performance by all metrics investigated (see Table 2), including respective test-set R² of 0.979 and 0.990. The Roost model also performs well with a test-set R² of 0.935, but it also has a non-normal error, as can be seen in Figure 1. Finally, we find that while XGBoost is the worst performing method, it still has a relatively good test-set R² of 0.866.

The performance of all 4 models is summarized in Figure 1. Visually, we find a very tight fit by the TPOT model in both the training and test sets, with good correlation from the XGBoost and SISSO models. We also find a systematic under-prediction of perovskite volumes in the roost model in both the training and test set, with the under-prediction beginning at approximately 75 ų / formula unit, achieving a maximum deviation at approximately 130 ų / formula unit, and returning to parity at approximately 200ų / formula unit.

The good performance of the TPOT model results from a generated pipeline with four stages. The first two stages are MinMaxScaler units; although the second one is redundant, as the data is already scaled to be between 0 and 1 by the first scaler. The third stage is a OneHotEncoder unit, which leverages one-hot encoding for categorical features (the TPOT implementation defines a categorical feature as one with fewer than 10 unique values). Finally, the perovskite volume is predicted using support vector regression.

In the case of the XGBoost model, we can extract feature importances. Although various different feature importance metrics can be derived from XGBoost, in this case we use the “gain” metric, which describes how the model’s loss function improves when a feature is chosen for a split while constructing the trees. A large number of features (290) were input into this model, so we display only the 10 most-important features identified by XGBoost in Supporting Information Figure S1. Here, we find that the average Rahm atomic radii ^{79, 80} have the highest importance score, followed by the average Van der Waals radius used by the Universal Force Field (UFF) ⁸¹. The remaining 288 features fall off as a long tail of low importance scores, indicating that they did little to improve the model’s performance in predicting the perovskite volume.

For SISSO, we used reduced the feature space as outlined in Section 2.4.4, with the pre-screened features listed in the Supporting Information along with the assumption we made about the units of the descriptor when fed into SISSO. Generally, we find that the main descriptors selected by the procedure are related to volume and atomic radius. Some other descriptors with less interpretability are found, such as the C6 dispersion coefficients, polarizability, melting points, and Herfindahl-Hirschman Index (HHI)⁸² production and reserve values. Although typically used to help indicate the size of a company within a particular sector of the economy, the XenonPy definition of HHI appears to come from the work of Gaultois et al ⁸². In the referenced work, the HHI production value refers to the geographic distribution of elemental production (e.g. answering the question of concentrated the industry that produces pure Fluorine is), and HHI reserve value describes the geographic distribution of known deposits of these materials (e.g. whether they are spread out over a wide area, or concentrated in a small area).

We report the best descriptor found in Equation 1. In this equation, the variables $c_{0},a_{0},a_{1}$ are the regression coefficients determined by SISSO.

where $c_{0}=-10.547$, $a_{0}=4.556$, $c_{1}=3.050$, $Z^{ave}$ is the average atomic number, $C^{ave}$ is the average mass specific heat capacity of the elemental solid, $r^{ave}_{Slater}$ is the average atomic covalent radius predicted by Slater, $r^{ave}_{pyykko,triple}$ is the average triple bond covalent radius predicted by Pyyko, $r^{ave}_{pyykko}$ is the average single bond covalent radius predicted by Pyyko, and $V^{ave}_{gs}$ and $V^{min}_{gs}$ are the average and minimum ground state volume per atom as calculated by DFT. Unsurprisingly the ground state atomic volumes and covalent radii play an important role in determining the final volume of the perovskite structures. Interestingly, both the atomic number and specific heat capacity of the material appear in the final descriptor; however, these likely only act as minor corrections to the final results.

The bandgap predictions leveraged a data filtering strategy (described in Section 2.3. As a result of our data filtering approach, the 6,351 entries in the dataset were reduced to 1,412 entries. The train/test split divided the data into a training set of 1,270 rows, and a test set containing 142 entries. The performance metrics of the XGBoost, TPOT, Roost, and SISSO models of 2D Material Bandgap can be found in Table 3. Performance is generally worse on this problem when compared to the perovskite volume predictions. As a result, in addition to the compositional features of XenonPy (Section 2.2.1) we also used several structural features (section 2.2.2). We also leveraged the bulk bandgap of the parent-3D material for each of the 2D materials, as we observed the performance of the TPOT, SISSO and XGBoost models increased when this value was included.

Although test-set model performance was worse compared to the perovskite problem, the TPOT and SISSO models retained their status as having the best performance metrics for the test-set R², MAE, and RMSE. The XGBoost model additionally performed well, with the best performance in terms of maximum error, and nearly as good performance on the other metrics when compared to TPOT and SISSO. We find the Roost model overfit the data to some extent on the data, as the test-set error metrics are considerably worse than their training-set counterparts.

A parity plot summarizing these results can be found in Figure 2. In all cases, we can see a spike of misprediction for systems with a DFT bandgap of 0. We note here that a large portion of these entries had DFT bandgaps of 0: of the 382 of the 1,412 entries in the dataset, a total of 27% of all training data.

The pipeline generated by TPOT is less complex than that of the perovskite volume problem. The first stage of the pipeline is a MinMaxScaler unit, scaling each feature such that it is between 0 and 1. The second stage is then an ElasticNetCV unit, which uses 5-fold cross-validation to optimize the alpha and L1/L2 ratio of the Elastic Net model. The converged alpha value was $1.011*10^{-3}$, and the converged L1/L2 ratio was 0.95, which strongly leans towards the L1 ( Least Absolute Shrinkage and Selection Operator (LASSO)) regularization penalty.

We can also extract feature importances from the XGBoost model, and we report the 10 highest-ranked features in Supporting Information Figure S2. In contrast with the perovskite results, the features are generally ranked similarly; although there is a clear ranking, we do not see 1-2 features dominating followed by a long tail of unimportant features. The selected features are also less interpretable; we find the minimum Mendeleev number as the most important feature, an alternative numbering of the periodic table in which numbers are ordered by their group rather than period. Specifically, we use the variant proposed by Villars et al ⁸³. For example, in the referenced system of numbering, Li, Na, K, and Rb are elements 1, 2, 3, and 4. This numbering is generally consistent, with the exception is H, which is assigned a value of 92 to place it above F, which has value 93. By taking the minimum value of the composition, this metric can be intuited as identifying the element in the earliest group and earliest period of the periodic table. Other descriptors we find are the average atomic weight, the variance of the number of unfilled s-orbitals in the system, and DFT-calculated ground-state atomic volumes.

The results of the prescreening procedure for SISSO are presented in the Supporting Information Table S3. These features are similar to the previous results, with the addition of the parent-3D material bandgap and electronegativity information playing an important role as well.

The selected SISSO model is

where $c_{0}=0.143$, $a_{0}=8.054\times 10^{-3}$, $r^{min}_{vdw}$ is the minimum Van der Waals radius of the atoms in the material, $E_{Bandgap}^{3D,parent}$ is the bandgap of the 3D-parent material, and $Period^{min}$ is the minimum period of the elements in the material. This descriptor represents a simple rescaling of the bandgap of the 3D-parent material, which makes sense as both bandgaps are highly correlated to each other. Furthermore this results is consistent with the TPOT model which is also primarily controlled by the bandgap of the 3D-parent material.

In the case of the 2D material exfoliation energy problem, the training and test-set statistics for the XGBoost, TPOT, Roost, and SISSO models can be found in Table 4. In this case, our feature selection methodology down-selected the 6,351 rows of our dataset into 3,388 rows. The train/test split further divided this into a training-set of 3,049 entries, and a test set of 339 entries. Generally, we see the worst performance of the models in this problem, compared to the perovskite volume and 2D material bandgap problems.

A set of parity plots for all four models is presented in 2. To facilitate easier comparison at experimentally-relevant energy ranges, we have zoomed the plot in such that the highest exfoliation energy is 2 eV. Plots showing the entire energy range explored can be found in Supporting Information Figure S6 Here, we find that all models perform generally poorly, with the largest errors occurring at higher exfoliation energies in the case of XGBoost, TPOT, and SISSO. (see Figure 3). The best test-set R² and RMSE this time is only TPOT, although they are still relatively poor, with a test-set R² of only 0.543. Roost displays the best test-set MAE, although the model seems to have overfit, as it displays drastically better performance on the training set than it does on the test set. The XGBoost model performs slightly worse than either TPOT or Roost, and the SISSO approach did not perform well for this problem

The TPOT algorithm again results in a relatively simple model pipeline. The first stage is a SelectFwe unit, which down-selects the features according to the Family-wise Error (FWE) ²³. An alpha value of $0.047$ is selected for this purpose, with results in a down-selection of the 299 input features into 125 features. This is then fed into an ExtraTreesRegressor unit, which is an implementation of the Extremely Randomized Trees method proposed by Geurts, Ernst, and Wehenkel⁸⁴.

We again extract features from the XGBoost model (Supporting Information Figure S3), and find the Mendeleev Number again appears as an important feature, albeit as the maximum instead of the minimum. Additionally, we see descriptors related to bond strengths in the corresponding elemental systems: average melting points, and average heats of evaporation.

The list of preselected features can be found in the Supporting Information Table S4. Overall, we see that this set of features is the least similar out of all three problems with the bulk modulus, thermal conductivity, and decomposition energy being the most prevalent. Additionally, the bandgap of the material is also selected, suggesting the possible hierarchical learning. For this calculation additional features such as the decomposition energy, ratio between the perimeter and area of the surface, and the electronic bandgap of the material are also included.

The best SISSO model found for this problem is

where $c_{0}=0.327$, $a_{0}=1.24\times 10^{-4}$, $a_{1}=9.26\times 10^{8}$, $\frac{P}{A}$ is the perimeter to area ratio of the surface, $q_{evaporation}^{min}$ is the minimum atomic evaporation heat of each element in the material, $E_{decomposition}$ is the decomposition energy, $\kappa_{i}$ is the thermal conductivity of the elements at 298 K, $V_{ICSD}^{ave}$ is the average atomic volume in the ICSD database, and $r_{vdw,uff}^{min}$ is the minimum atomic Van der Waals radius from UFF for each element in the material.

In the case of the volume per formula unit of ABX₃ perovskites, we observe all four model types perform well. Overall, we find that the volume per formula unit for ABX₃ perovskites can be predicted using only compositional descriptors (i.e. with no structural descriptors). The likely reason all four models perform well despite having no structural information is the general similarity in crystal structure between these systems — they are all perovskites, and therefore all possess very similar crystal structures. Supporting this is that the Roost model, which only leverages the chemical formula as an input, and which we did not optimize the hyperparameters or architecture for, performed just as well on this problem — albeit with some systematic deviation from parity at intermediate volumes. Although interpretability is reduced by virtue of being a neural network, we can still achieve an important insight from this model — just by knowing the chemical formula of the system, we can achieve accurate predictions of perovskite volumes, which further justifies our use of compositional descriptors (see Section 2.2.1) on this problem as we move to the SISSO, XGBoost, and TPOT models. Additionally, we note this performance was achieved with a training set containing only 129 entries in the dataset - compared to the original Roost paper³¹ that leveraged approximately 275,000 entries from the OQMD datset⁸⁵.

Like the Roost model, we have difficulty in interpreting the pipeline generated by TPOT. The TPOT model delivers the best performance — which is clearly visible from the parity plot in Figure 1. This performance came at a price, however, and the rather complex pipeline containing scaling, one-hot encoding, and linear support vector regression does not lend itself well to interpretation.

Entering into the realm of interpretability, although the XGBoost model does not produce a direct formula for perovskite volumes, we can still gain some insight using it. It is still, however, relatively accurate — and allows us access to a feature importance metric (see Supporting Information Figure S2). In this case, we see the five most-important features are the average Rahm^{79, 80} atomic radius, average UFF⁸¹ atomic radius, sum of elemental velocities of sound in the material, average Ghosh⁸⁶ electronegativity, and the sum of the Pyykko⁸⁷ triple-bond covalent radii. Overall, we see a strong reliance on descriptors of atomic radius — which as we noted in the TPOT discussion makes intuitive sense.

Finally, the SISSO model (Equation 1) offers the most direct interpretation, as it is simply an equation. Immediately, we can see that ground state DFT atomic volumes are important. This result is highly intuitive and is not surprising when we consider that i) we are predicting volume, therefore using an average atomic volume makes sense and ii) the TPOT and XGBoost models extensively leveraged atomic radius descriptors that are also related to the volume of the atoms.

The overall good performance of SISSO for this application is promising, as it is one of the most accurate models, while being by far the most interpretable. This represents a key advantage to symbolic regression, as if you can find an accurate model, then it will be easy to understand and analyze the results. Moreover, we note that are not alone in the literature when it comes to leveraging SISSO to generate models of perovskite properties — the last several years have seen success in the creation of models of perovskite properties with this tool. The work of Xie et al. ⁶⁰ achieved good success in predicting the octahedral tilt in ABO₃ perovskites, the work of Bartel et al. ⁵⁹ resulted in the creation of a new tolerance factor for ABX₃ perovskite formation, and Ihalage and Hao ⁵⁸ leveraged descriptors generated by SISSO to predict the formation of quaternary perovskites with formula $(\text{A}_{1-x}\text{A}^{\prime}_{x})\text{BO}_{3}$ and $\text{A}(\text{B}_{1-x}\text{B}^{\prime}_{x})\text{O}_{3}$.

The 2D material bandgap models did not achieve the same performance as for the perovskite systems (see Figure 2). Even in the case of TPOT and SISSO, which still had the best test-set performance by most metrics, there were a few outliers. Specifically, we find that the test-set MAE for the models ranged between 0.29 eV (TPOT and SISSO) and 0.65 eV (Roost) relative to the PBE DFT calculations reported by the 2DMatPedia. Putting this number in perspective, we note the recent work of Tran et al⁸⁸, which benchmarked the bandgap predictions of several popular DFT functions for many of the systems in the C2DB; the work identified that the PBE functional exhibited a MAE of 1.50 eV relative to the G₀W₀ method. Other investigators have studied the prediction of 2D material bandgaps: Rajan et al.⁸⁹ also achieved a test-set MAE of 0.11 eV on a dataset of 23,870 MXene systems (which, as far as we are aware, has not been made publicly available) using a Gaussian Process regression approach, with DFT-calculated properties including the average M-X bond length, volume per atom, MXene phase, fand heat of formation, and compositional properties including the mean Van der Waals radius, standard deviation of periodic table group number, standard deviation of the ionization energy, and standard deviation of the meting temperature. Zhang et al.⁹⁰ improved on this error slightly, achieving a test-set MAE of 0.10 eV on the C2DB dataset (around 4000 entries)¹⁴ with both Support-Vector Regression and Random Forest approaches, albeit using descriptors such as the fermi-energy density of states and total energy of the system (requiring further DFT work for additional prediction). In contrast to both approaches, which used DFT-calculated values that would need to be obtained for new systems to be predicted, the only DFT-calculated value we leverage in our feature set is a bulk bandgap tabulated on the Materials Project¹⁸. Thus, although our TPOT model had a slightly higher MAE of 0.29 eV, we note that this would not require further DFT work to generate new predictions. As 2D systems are still relatively new, we note that much more work has been performed in the 3D materials space, particularly in the leveraging of neural networks to predict bandgaps. The recent Atomistic LIne Graph Neural Network (ALIGNN)⁹¹ reported a test-set MAE of 0.218 eV for the prediction of bulk materials hosted by Materials Project¹⁸ (which as of October 2021 has over 144,000 inorganic systems). The Materials Graph Network (MEGNet) architecture⁹² achieved a test-set MAE of 0.32 eV on the bulk systems of the Materials Project. Although these neural network models are on 3D systems, we note that they do not leverage DFT properties (which we re-iterate would cause any resulting model to require a DFT calculation for future prediction) and had access to much larger datasets than the training set we obtained after filtering the 2DMatPedia entries (see Section 2.3). Overall, although the systems we investigate are not 3D bulk systems, we believe this puts the TPOT MAE of 0.29 eV for the bandgap of 2D systems in perspective.

In all 4 models we trained, many of the incorrect predictions occur where the DFT bandgap is 0 eV (which represented 27% of the training set values). Because of this, we tried simplifying the bandgap problem, by training an XGBoost model to predict whether the system was a metal (see Supporting Information section S.6.3), and showed that we could achieve good results — although we incorporated a purely structural fingerprint, the Sine Matrix Eigenspectrum (see Supporting Information Section S.6.1). As this descriptor resulted in some rather large vectors (of length 40, the maximum number of atoms in any system) with little direct physical intuition, we do not directly include it for the purposes of this section. Instead, what we can derive from it is the knowledge that structural features can provide information to predict the bandgap.

If we take a closer look at the Roost model, which takes a purely compositional approach to materials property prediction, we can see a poor generalization to the test set (see Table 3). This indicates that we have likely caused it to over-fit (which could have been improved for example through the use of early stopping). Overall, we may also use this result to further show the need for structural descriptors when predicting bandgaps of these systems.

The TPOT, SISSO, and XGBoost models for 2D material bandgap achieved similar performance to one-another (see Table 2). Moreover, we again have the opportunity to extract meaning from the TPOT model here, as it leveraged Elastic Net to perform its prediction — a model which is achieved by mixing together Ridge Regression and the LASSO. Additionally, as the L1/L2 ratio is 0.95, this is nearly entirely L1 regularization (see section 3.2). As a result of this, and that the data were all scaled such that they ranged between 0 and 1, we can view the the feature coefficients approximately (but not entirely) through the lens of LASSO, which does perform feature selection. Doing this, we see that the corresponding bulk 3D-parent material’s bandgap (as found on Materials Project) is the most strongly-weighted feature by the elastic net model. This is extremely intuitive as well — it is reasonable for the bulk bandgap to be correlated with the 2D bandgap.

The XGBoost model, however, suggests an entirely different set of descriptors as being “important.” Instead, we find that the five most-important features are, in order from greatest to least, the minimum Mendeleev number, average atomic weight, variance of the number of unfilled s-orbitals, average ground-state DFT volume, and maximum Cordero⁹³ covalent radius. None of these descriptors have a particularly intuitive relationship with the bandgap either. This is a similar result to that found by Rajan et al. ⁸⁹, who leveraged several models to predict the GW bandgap of MXene systems; in this work, they found that several of the most-important features exhibited only a statistical (i.e. non-intuitive) relationship with the bandgap.

Finally, the SISSO model (Equation 2) shows an equivalent performance to TPOT and XGBoost, while using only three primary features. Similar to the TPOT model the bandgap of the bulk 3D-parent material is the most important feature, with only minor non-linear contributions from the minimum atomic Van der Waals radius and period. Interestingly, the SISSO model not only misses some of the metallic materials, but also incorrectly predicts some materials to metallic in the training set. This suggests the increase simplicity of the SISSO models may slightly reduce their reliability.

Future work on this problem may achieve better performance on the bandgap problem by investigating the bond strengths of the different elements in the system. We also note the very good performance that recent neural network approaches have had on the 3D bandgap problem^{91, 92}, likely due to their representation of the structure of the 3D systems. Similar to how the Roost model achieves good success when compositional descriptors are appropriate, we may find good success in leveraging neural network approaches when structural features are required. We note here that Deng et al.²⁸ achieved good results on a variety of molecule properties by incorporating various graph representations from different neural network architectures. Hence, future work in this domain may benefit from the incorporation of the information-dense structural fingerprints that may be obtained from neural network-based approaches.

We observed some of the worst model performance (across all models) in the case of the 2D material exfoliation energy. Despite being a larger dataset than either the perovskite (144 total, 129 in the training set) or bandgap (1,412 total, 1,270 in the training set), the 3,049 entries in the training set (out of 3,388 total) for the exfoliation energy proved insufficient to achieve good results for any of the models. Moreover, the compositional and structural features were not sufficient to adequately describe the system.

Some interpretation can at least be derived from the XGBoost and TPOT models. The five most-important features according to XGBoost (see Supporting Information Figure S3) are, in order, the maximum Mendeleev number, average melting point, average evaporation heat, maximum number of unfilled orbitals, and the sum of the melting points. In the case of the TPOT model, we arrived at an extremely randomized tree approach, which also has a feature importance metric. Here, we find that the average Van der Waals radius, maximum dipole polarizability, minimum atomic weight, maximum atomic weight, and maximum elemental (and tabulated) DFT bandgap are weighted as important. Between the two models, we see much difference in which features are deemed important. In the XGBoost model, the average melting point, sum of melting points, and average evaporation heat of the elemental systems can be rationalized if we realize that these are all driven by the strength of the interaction between atoms in the material, thus these descriptors may provide information relating to the forces that must be overcome when exfoliation is performed. Many of the other features in the two models correlate with size: the maximum Mendeleev number, average Van der Waals radius, and minimum / maximum atomic weights all provide a description of the size of the atoms in the system.

Finally, the SISSO model (Equation 3) echos these findings. Although it is less performant, it provides intuitive descriptors that tell a similar story. One term takes the ratio of the decomposition energy and approximations to the atomic volume. This captures a description of atomic size as well as the strength of atomic bonds. Additionally, the second term in the model incorporates information about the surface, thermal conductivity and the heat or evaporation. Although the second term is less descriptive, it still captures terms relating to the size of the atoms and bond strengths of that atoms involved in the 2D material. The relatively poor performance of the SISSO models, also highlights the need for better input features to describe the exfoliation energy of a material. In cases where only a loose correlation between a target property and the inputs exist, the relative simplicity of symbolic regression models can be detrimental. While TPOT and XGBoost can utilize information from all features in their final predictions, symbolic regression in general and SISSO in particular only acts on the order of tens of features. Because of this, it is likely that more structural information is needed to get accurate models with SISSO.

In effect, the models can be interpreted as collectively describing the exfoliation energy as a function of i) the size of the atoms in the 2D material and ii) the strength of the intermolecular forces between those atoms. When we predict exfoliation energies, we’re predicting the interaction between layers in an exfoliable material. Overall, finding better methods of cheaply approximating these weak interactions may provide better results in the prediction of exfoliation. Additionally, as the number of datasets which contain exfoliation energies increases (such as the 2DMatPedia¹¹, C2DB^{14, 15} and JARVIS¹⁷), further insight into this problem will be possible, and more-complex (albeit less-interpretible) models will become feasible.

Additionally, in order to obtain more-accurate predictions of exfoliation energy, data generated via a more thorough computational treatment may be required. We illustrate this by examining an outlier in the training set at 9.9 J / m², which all four models heavily under-predicted (see Supporting Information figure S6) (7-8 eV in the case of XGBoost and TPOT, and over 9 eV in the case of Roost and SISSO). Upon closer examination of this system, we find that it is actually a pair of layers containing N atoms (Figure 4 A). The 2DMatPedia ¹¹ reports that this system (2dm-id 5985) was not directly sourced by a simulated exfoliation from a bulk structure, but instead was obtained by substituting the atoms in a hypothetical 2D Sb structure (Figure 4 B). The Sb structure (2dm-id 4275) was obtained by a simulated exfoliation from a structure obtained from materials project (Figure 4 C). The parent bulk material (mp-567409) is reported by the Materials Project¹⁸ to be a monoclinic crystal which undergoes a favorable decomposition (energy above hull is reported as 0.121 eV / atom) to a triclinic system. That being said, as this is a hypothetical 2D system, comparison with the hypothetical 3D bulk system was necessary for the calculation of exfoliation energy. As the prediction of crystal structure is a very challenging field with few easy approximations⁹⁴, this may have contributed further to the extreme value of the exfoliation energy. Indeed, as Zhou et al report¹¹, the decomposition energy lends itself better to assessing whether a material is truly stable. Indeed, despite the extremely high exfoliation energy of this hypothetical 2D N system, it is reported by the 2DMatPedia to have a decomposition energy of 0 eV/atom. This too seems somehwat high, as systems containing N-N bonds tend to be high-energy materials, typically undergoing strongly exothermic decomposition to inert, gaseous N₂ ⁹⁵. With this in conjunction with the observation that our models all predict exfoliation energies significantly lower than the tabulated values, we have reason to believe that this system would be far easier to exfoliate than the  10 eV exfoliation energy implies. Moreover, this system may have a strong energetic preference to decompose further into N₂, which additional DFT work could reveal. Overall, this underscores the importance of obtaining high-quality data, and filtering that high-quality data, for the training of interpretable models.

As ML is further integrated into materials discovery workflows, we anticipate that the numerous successes neural networks have presented^{32, 33, 29, 30, 31, 34, 35, 36, 37, 38, 39, 40, 2, 41, 1, 3} will continue to propel them onto the cutting edge of chemical property prediction. This comes with the challenge of honing our techniques for their interpretation, an area which has seen much interest in recent years, and where there is still plenty of opportunity for further development^{96, 97}. We also expect AutoML techniques such as TPOT will continue gaining traction in materials discovery, due to the amount of success and attention they have recently had ^{42, 43, 44, 45, 46, 47}. This too presents the challenge of interpretability if highly complex pipelines are generated (see Sections 3.1 and 4.1). We note here that part of the value that AutoML techniques bring is the ability to make advanced techniques accessible to a wider audience of researchers by lowering the barrier of entry. Hence, we expect that the problem of interpretation may be compounded for AutoML (and especially NAS) systems: the ability to automatically extract some level of interpretation from the generated pipelines is important for automation to make ML truly accessible to non-experts. Overall, we expect that as neural network models and AutoML algorithms continue to grow in capability and complexity, work in developing the tools and techniques needed to interpret them will see a greater attention.

In contrast with challenge of interpreting neural networks or the pipelines found by AutoML systems, symbolic regression tools like Eureqa and SISSO yield an exact equation describing the model, and are thus easier to interpret. This makes it easier to achieve key insights with physical interpretations — such as the very intuitive way in which SISSO is able to describe the systems. Overall, despite its reduced ability to predict the exfoliation energy of a material when compared to the models of TPOT, XGBoost, and Roost, we note the mathematical equations returned by SISSO provide a direct relationship between the target properties and model predictions. Additionally, in the case of the exfoliation energy, we believe that we may see further improvements by including richer structural information. We base this on the observation that the Roost model performed poorly on both of this problems – recalling that Roost is only provided the chemical formula of the system, this could indicate that compositional descriptors alone are insufficient to describe these properties. Indeed, it is well-known that structure and energy are intimately related (the fundamental assumption of geometry optimization techniques is that energy is a function of atomic position), hence it can be inferred that exfoliation energy and structure are similarly related. In the case of bandgaps, we note that there is also a strong dependence on structure; Chaves et al⁹⁸ notes that the number of layers in a 2D material can strongly influence the band gap, reporting differences of up to several eV can occur between the bulk and monolayer form of a material.

Interoperability is still a challenge in the materials discovery ecosystem. Although it is possible to easily convert between different chemical file formats (e.g. by OpenBabel⁹⁹), and packages such as Pymatgen ⁶⁷, Atomic Simulation Environment (ASE) ¹⁰⁰, and RDKit ¹⁰¹ can easily convert to each others’ format, we note that there is a challenge of calculating features using a variety of different packages. Some tools expect Pymatgen objects (e.g. XenonPy), others expect ASE objects, whereas others require RDKit objects (e.g. all of the descriptors in the RDKit library) to perform a calculation of features, thus creating some standard for the interoperability of these packages would be beneficial. Additionally, further efforts should be made to report the sources of data used by featurization packages. We note that MatMiner ⁷⁰ is exemplary in this regard: each of the featurization classes it defines has a “citation” method returning the appropriate source to credit. Mendeleev ⁶⁶ is another good example of this; within its documentation, a table lists citations for many (though not all) of the elemental properties it can return. Overall, by placing a stronger focus on i) interoperability and ii) data provenance, the Python materials modeling ecosystem can be made stronger — and therefore help accelerate materials discovery.

All of the models we have investigated in this work required sufficient training data to avoid over-fitting. Although techniques such as cross-validation, early-stopping (in the case of neural networks and XGBoost), and train/test splitting can help guard against (and detect) over-fitting, having a sufficiently-large dataset is of the utmost importance to achieve truly generalizable models. As a result, there is a critical need for data management approaches that satsify the set of FAIR principles. This crucial need for effective data management has led to the incorporation of data storage tooling in popular chemistry packages including Pymatgen⁶⁷, ASE¹⁰⁰, and RDKit¹⁰¹. Moreover, advances in both computational capacity and techniques has given rise to studies performing the high-throughput screening of chemical systems^{102, 103, 104}. This has resulted in the development of tools focusing on the provenance of data, such as the Automated Interactive Infrastructure and Database for Computational Science (AiiDA) system ^{105, 106}.

Overall, we have identified a series of key issues should see more attention as the digital ecosystem surrounding materials modeling continues to develop. First, interpretability of models allows us to derive physical understanding from the available data. This is a key benefit of symbolic regression tools like SISSO, which result in the creation of human-readable equations describing the model. Additionally, increasing the accessibility of ML techniques through automation (such as in the field of AutoML) will allow a wider range of researchers the ability to benefit from advances in modeling techniques. Data management and data provenance is another major issue, which allows us to better understand which datasets can be combined (e.g. when combining DFT datasets, the methodologies should be consistent between them), and to help us understand if something intrinsic to the training data is affecting model performance. These data management goals are core focus of platforms such as Exabyte ¹⁰⁷, which provides an all-in-one solution for i) storing material data and metadata, ii) storing the methodology required to derive a property from a material, and iii) providing the means to automatically perform calculations, and iv) automatically extracting calculation results and storing them for the user. This focus on providing a tool that manages materials, workflows, and calculations has allowed Exabyte to be a highly successful platform, which has led to studies involving automated phonon calculations¹⁰⁸, high-throughput screening of materials for their band-structure ^{109, 110}. Future capabilities of the platform are slated to include a categorization scheme for computational models to provide even more metadata to track the provenance of calculated material properties¹¹¹.