Extracting an Empirical Intermetallic Hydride Design Principle from Limited Data via Interpretable Machine Learning

An open question in the metal hydride community is whether there are simple, physics-based design rules that dictate the thermodynamic properties of these materials across the variety of structures and chemistry they can exhibit. While black box machine learning-based algorithms can predict these properties with some success, they do not directly provide the basis on which these predictions are made, therefore complicating the a priori design of novel materials exhibiting a desired property value. In this work we demonstrate how feature importance, as identiﬁed by a gradient boosting tree regressor, uncovers the strong dependence of the metal hydride equilibrium H 2 pressure on a volume-based descriptor that can be computed from just the elemental composition of the intermetallic alloy. Elucidation of this simple structure-property relationship is valid across a range of compositions, metal substitutions, and structural classes exhibited by intermetallic hydrides. This permits rational targeting of novel intermetallics for high-pressure hydrogen storage (low-stability hydrides) by their descriptor values, and we predict a known intermetallic to form a low-stability hydride (as conﬁrmed by density functional theory calculations) that has not yet been experimentally investigated.

Development of renewable energy technologies is more critical now than ever to avoid some of the catastrophic consequences of climate change. 1 Hydrogen is a clean energy carrier poised to make an impact throughout the renewable energy space, but storage and transportation of hydrogen gas (H 2 ) remains a significant challenge. 2,3 Decades of research have been devoted to storing hydrogen more economically/efficiently, and metal hydrides [4][5][6][7][8][9][10] are one of the most extensively studied materials for applications in H 2 storage for transportation, H 2 compressors, thermal energy storage, and H 2 getters. [11][12][13][14] Their practical applicability varies widely as a function of their thermodynamic properties which, when combined with other factors such as sustainability, cost, kinetics, capacity, has lead to thousands of metal hydrides being investigated experimentally. Thus, an open question is whether there exist simple materials design rules that dictate their thermodynamic properties across their varying chemical and structural space. Correlations have been elucidated from various experimental results [15][16][17][18][19][20] and empirical design rules derived, [21][22][23][24][25] such as the pressure dependence on interstitial volumes for a given intermetallic series. Computational screenings have also been performed, [26][27][28] but this problem has received comparatively little attention from a "big-data" perspective. In other energy related applications such as hydrogen storage or xenon/krypton separations in porous materials, big data approaches have been able to identify relatively simple materials descriptors and models that predict thermodynamic performance across a wide swathe of material space. 29,30 Can similar results be achieved for metal hydrides?
Statistical or machine learning (ML) techniques have the potential to answer this question and, despite sometimes lingering skepticism over their utility, are now continually employed to solve problems in the physical sciences. 31 Some prominent examples include generative models for drug design, 32 prediction of conductive metal organic frameworks, 33 classification of stable perovskites, 34 among others. [35][36][37][38][39][40] One natural way to segregate these studies is by those that use "black-box" vs. "explainable" ML techniques. 41 Black-box ML techniques are well-suited to make accurate predictions of materials properties but provide little visibility into how the algorithm utilizes feature space to make a prediction, a potential limitation that explainable ML techniques seek to address. For example, "explainable" insight can be derived by simply extracting a given feature's importance when training a model 42 or by using ML methods whose mapping from features to prediction is directly interpretable by design. 43,44 A few studies have applied ML techniques to make black-box predictions on the thermodynamic properties of metal hydrides. Rahnama et al. trained a model to use measured properties of metal hydrides to predict other measured properties. 45,46 However, this is not particularly predictive since one would have to perform an experiment or simulation to measure the materials' properties to use the model, and at that point the properties would be known, negating the need for a model. One of their main conclusions from this approach is that "composition formula was found to be an insignificant variable", 45 which is surprising given the large body of literature on doping and destabilization of metal hydrides.
Hattrick-Simpers et al. trained a model on the Department of Energy's experimental metal hydride (HydPARK) database to predict hydriding enthalpies solely from the composition of the intermetallic phase, which was then used as a surrogate model to quickly evaluate the performance of novel intermetallic compositions for use in hydrogen compressors. 47 These studies did not exploit insights from explainable ML to determine what properties of intermetallic compounds dictate their thermodynamic performance.
In this work we also train an ML model on the HydPARK database using features derived solely from the intermetallic composition (no structural or hydride information); however, our major contribution is to use feature importance from gradient boosting trees to gain "explainable" insight into simple structure-property relationships that govern the thermodynamics of hydride formation. While our ML model can accurately predict the room temperature equilibrium H 2 pressure of intermetallic hydrides, its interpretability allows us to generalize the pressure dependence on the lattice volume in the LaNi 5 substitution series (a historically known design correlation [18][19][20] ) over a surprisingly wide range of intermetallic chemistries and structures. This unifies disparate experimental results onto a single structure-property relationship. Its elucidation provides thermodynamic insight into the underpinnings of the ML model predictions, which we further corroborate with density functional theory (DFT). We then utilize this to predict a known intermetallic for high-pressure H 2 storage applications whose hydriding properties have not yet been experimentally tested. impractical to compute for hundreds or thousands of structures. In contrast, the HydPARK database contains thermodynamic data that is not easily calculated but readily measurable, such as the equilibrium pressure of H 2 , P eq , at a given temperature, T . Therefore, ∆S can be computed using experimental data and the van't Hoff relationship, Often we are interested in predicting the equilibrium H 2 pressure, ln P o eq (@25 o C), as it indicates how much H 2 a material can deliver at room temperature and provides a standardized metric for comparing metal hydrides that accounts for both entropic and enthalpic hydriding effects. 51,52 Figure 1 shows the strong enthalpy-entropy correlation (which is not unique to this application 53  tains many complex stoichiometries for which the exact crystal structure is not known and therefore cannot be included in a computational database of specific crystal structures. We therefore proceed with the HydPARK database for developing our ML model.
Next we clean and prepare the HydPARK database before training an ML model with additional details included in the Supplementary Information section S1. Briefly, we remove compositions with missing or unusable data (for which ln P o eq cannot be calculated), thereby reducing the size of the dataset from 2732 entries to 570 entries. Further investigation reveals significant spread in the reported experimental data for duplicate compositions (e.g. 6 different CaNi 5 entries) as well as incorrectly collected data in HydPARK, something we expect will challenge the development of a highly accurate ML model. We therefore remove duplicates, yielding 409 unique compositions, while minimizing the bias introduced by this process (more details in S1); however, these literature references need to be revisited individually and experiments repeated when necessary. Additionally, there is a large imbalance in both the distribution of thermodynamic properties (Table 1) and the population sizes of different metal hydride classes. For example, less than 3% of complete and unique entries are complex hydrides 55,56 which we discard since these ∼10 sample points are insufficient for an ML model to learn from. S12 contains our final version of the "ML ready" HydPARK database. Table 1: The distribution of ln P o eq for complete and unique compositions in Hyd-PARK. Even within the ln P o eq > −10 subset, the data is non-uniformly distributed as shown in Figure 2d.

%
Machine learning with feature importance.-The Magpie code 57 was used to generate a set of 145 features derived solely from the intermetallic chemical composition for each HydPARK material. Therefore, no structurally specific features were included for training the ML model (other than what is implicitly encoded by the properties of the material's constituent elements), an approach which has shown great success in a variety of materials science applications. 57,58 In other words, we try to discover whether the thermodynamic properties of intermetallic hydrides can be a priori predicted from the intermetallic composition without any information about the hydride composition or structure. Next a gradient boosting tree regressor (the best performer in comparison to other regression techniques as shown in S2) was trained using scikit-learn 59 to predict ln P o eq . A 10-fold validation was performed, and the combined test and train sets for each of the 10 models is shown in Fig where f i is the composition fraction of element i, and ν i is the volume occupied per atom in the ground state elemental solid of species i. In other words, it is Magpie's estimation of the specific volume per atom for a given composition. Notably, its average relative importance across all 10-fold validation models is close to 100 (not exactly 100 since one model yields most_GSvolume_pa as the most important), indicating that it is the most important descriptor regardless of the test/train split; it even remains so when removing data from the training set (see S4). We discuss secondary features in more detail later as they constitute important features when training individual models to predict ln P o eq 's constituent components, ∆S and ∆H. Ultimately, the interpretable model suggests that a simple volume-based descriptor may be the single most useful feature, a surprising result given the wide ranging chemistries and structure types displayed by these materials.
An intermetallic hydride design principle.-The high importance of the ν M agpie pa descriptor warrants further investigation into a structurally specific volume descriptor. We cross-reference the compositions in the training set with the MP database to identify ∼80 overlapping structures for which we extract the DFT relaxed crystal structure corresponding to the lowest formation energy per atom. We then derive an analogous descriptor to ν M agpie pa based on the cell volume, V cell , and the number of atoms in the cell, n atoms , from the MP crystal structure Hence ν pa refers to the volume per atom in a crystal which can be either estimated (ν M agpie   Empirical modeling of binary alloy formation enthalpies 21,63,64 has utilized cellular models incorporating properties like electronegativity differences and the difference in electron density at the boundary between dissimilar atoms. Extension to ternary hydrogen-containing alloys often relied on knowledge of structurally specific features. 16,22,25 ∆H has also been rationalized in terms of a qualitative rule of reversed stability: substituting La or Ni to stabilize the binary intermetallic results in a less stable hydride phase. 22 Interestingly, the same insights from these different modeling efforts can be qualitatively rationalized across the HydPARK materials via our data-driven approach. Specifically, ν pa is a "synthetic feature" that encodes many of the other chemically specific features that affect ∆H; therefore, this feature can even be removed when training the ML model without significant loss in accuracy (S7). As a simple illustration in Figure 4d-e, materials with larger ν pa tend to have larger average pairwise electronegativity differences between elements (Magpie's MeanIonic-Char feature), as well as lower mean melting temperatures. We can rationalize such trends as indicators for increasing hydrogen absorption strength or decreased energy penalty for lattice deformation, which we investigate further with DFT.
These structure-property relationships are also invaluable for outlier identification. For example, note that the thermodynamic properties of materials with ν pa > ∼17Å 3 /atom deviate significantly from the trend exhibited by materials with ν pa < 17Å 3 /atom. These exactly correspond to the materials for which the model generalizes poorly (Figure 2b) due to a lack of materials in this descriptor regime, an insight only derived because the simple structure property relationships elucidated by the ML model's feature importance. This breakdown of the log-linear (ln P o eq ) or linear (∆H) correlation at this critical threshold suggests that these materials require a different physical understanding than the ν pa structure-property relationship, but the lack of data in this regime must be addressed before ML models have the chance to provide further data-driven insight. Other secondary benefits of the model's interpretability are discussed in S7, S8, and S10.
Thermodynamic insights from DFT.-We can further corroborate our insights into the ν pa structure-property relationship by examining A site substitutions in the LaNi 5 series with DFT. First we define E f as the formation energy of the intermetallic alloy with respect to the elemental crystals, ∆E def as the energy penalty required to deform the intermetallic lattice to accommodate H absorption, and ∆E H as the binding energy between hydrogen and metal atoms in the hydride lattice (see S11 for the definition). V is the volume of the hydrided lattice, and V 0 is the volume of the intermetallic lattice. Specific details on these calculations are provided in S11, [65][66][67][68][69][70][71] and we summarize DFT computed ∆H, E f , ∆E def and ∆E H in Table 2. Note that in these calculations the final state is the hydride and the initial state is the intermetallic, i.e. the ∆'s correspond to the hydriding reaction, not the dehydriding reaction. We also note that we considered a hydride composition of ANi 5 H 7 (A = U, Ce or La), which corresponds to a hydrogen/metal ratio of 1.17 and is close to the maximum hydrogen uptake reported for LaNi 5 . 71 This demonstrates the experimental trend of increased hydride stability (∆H) with ν pa while the individual energy terms yield additional insight. First, there is no apparent correlation with E f . Rather, increasing ν pa more importantly indicates a propensity for a lower energy penalty of deformation ∆E def , which can be rationalized by the reduced volume expansion required to form the hydride phase. It should be pointed out that in addition to the magnitude of volume expansion V /V 0 , the elastic modulus of the intermetallic will also affect the energy penalty of deformation ∆E def : assuming the same volume expansion, a stiffer intermetallic would require a bigger energy penalty to deform the lattice in comparison with a less stiff intermetallic. Second, the binding energy between hydrogen and metal atoms in the hydride lattice increases with ν pa , which can be rationalized by the lower electronegativity of La in comparison with U (1.1/La vs. 1.38/U on Pauling scale), i.e. the binding between hydrogen and LaNi 5 is expected to be more ionic than that between hydrogen and UNi 5 .
Both of these effects promote greater hydride stability, and it is these trends which underpin the ν pa structure-property relationship. The general trend of ∆H as a function of ν pa may also be inferred from simple "chemical intuition", i.e. intermetallic alloys with larger ν pa usually consist of elements with larger atomic radii, and elements with larger atomic radii tend to have lower electronegativities (see S9), because of reduced attractions to electrons in the outer valence shell. The result is that an intermetallic alloy with a larger ν pa value tends to form more ionic bonds with hydrogen in the hydride phase, and therefore the enthalpy of formation of the hydride, ∆H, is usually bigger than that of an intermetallic alloy with a smaller ν pa value. We note that for intermetallic alloys with similar ν pa values, the energy penalties of deforming the alloy lattice to accommodate the chemically absorbed hydrogen can be very different, which may result in a large scattering of the ∆H values.
ML informed targeting of a novel hydride phase.-We propose UNi 5 in Table 2 for two reasons. First, based on our ML informed results, we predict this A site substitution to LaNi 5 to reduce the stability of the metal hydride phase since it significantly reduces ν pa (U has a smaller atomic radius than La). Second, UNi 5 is an experimentally synthesized intermetallic 72 in the Inorganic Crystal Structure Database (ICSD), but its hydrided form has not yet been reported in the ICSD nor is it contained in the HydPARK database (potentially due to the large H 2 pressures that may be necessary to synthesize it near room temperature). As confirmed by our DFT calculations, UNi 5 H 7 has a very small reaction enthalpy of -0.6 kJ/(mol·H 2 ) and therefore should be a low stability hydride. As seen from Figure 1, UNi 5 H 7 , should it be synthesized in the future, would be one of the least stable hydrides in the entire HydPARK database. We note the hydriding reaction enthalpy of a metal alloy may differ depending on the amount of hydrogen that is absorbed. Nevertheless, the qualitative knowledge generated by our interpretable ML provides a path to rationally target novel hydride phases with a desired thermodynamic property, i.e. very low stability for high-pressure H 2 or hydrogen isotope storage applications. 73,74 In conclusion, utilizing the HydPARK experimental metal hydride database, we have trained an ML model to predict the equilibrium plateau pressure, one of the most relevant thermodynamic quantities for practical applications which is also unique to this database (i.e. not contained in any computationally derived databases due to its dependence on ∆S).
Exploiting the explainability of gradient boosting trees with our data-driven approach enables several key understandings. First, basic thermodynamic insight into intermetallic metal hydride formation can be derived from features generated only from the elemental composition of the intermetallic phase (a particularly useful capability if the exact crystal structure of an experimental compound is not known). Past experimental studies have elucidated the dependence of equilibrium H 2 pressure on cell volume or structurally specific interstitial volumes, and the identification of our ν pa structure-property relationship encompasses these observations across a range of intermetallic chemistries and structures. The thermodynamic basis for this correlation is attributed to the underlying structure-property relationship between ∆H and ν pa ; furthermore, materials not described by this simple structure-property relationship can now be investigated to determine the chemistry behind their outlying behavior. All of these insights are predicated on the physical interpretability of an ML model, which, when corroborated with DFT calculations, is ultimately used to propose a novel hydride of a known intermetallic with significant potential as a high-pressure hydrogen storage material.
Furthermore, we utilized a noisy, imbalanced database which required multiple heuristic steps to process and clean. This simply highlights that statistical learning techniques still have the power to help extract useful information in materials science applications, even when approximations must be made to prepare the training data. Future efforts in this space will benefit greatly from a concerted effort of the metal hydride community to centralize the reporting of experimental measurements such as ∆H, ∆S, P eq , T , V cell (if possible), etc. 75 There should also be a standardized method for reporting more complex phenomena such as hysteresis, existence of multiple hydride phases, sloping plateaus, etc. This could be incorporated into the framework of one of the many existing materials databases (MP, OQMD, AFLOWLIB) which would better position data-driven/ML based approaches to impact the discovery and understanding of metal hydrides. Less than 20% of the HydPARK database was used due to missing and/or duplicate information. Several errors were found in the dataset from a cursory manual investigation of the literature references therein, and standardized reporting in a central repository could help avoid such inconsistencies. Moreover, if more complete material entries existed in the HydPARK database with larger volumes per atom, our explainable ML approach might be able to elucidate the structure-property relationship(s) that differentiate them from the ν pa discussed in this work. However, this will be an unlikely accomplishment from an ML perspective until more/better data is acquired in this regime. Having gained explainable insights into the thermodynamics of hydride formation, future efforts can now be directed towards explainable ML models that discern whether a given composition (out of the essentially infinite number that may exist) will form a hydride and, if so, what its hydrogen content may be. We propose that combining all of these efforts will result in the data-driven discovery of novel, high-performing hydrides.