Design of lightweight BCC multi-principal element alloys with enhanced hydrogen storage using a machine learning-driven genetic algorithm

Abstract

Body-centered cubic (BCC) based multi-principal element alloy (MPEA) hydrides have demonstrated significant potential for compact and efficient hydrogen storage. In this work, we first leverage machine learning (ML) models to predict the hydrogen affinity, storage capacity and phase stability of BCC MPEAs, creating a unique hydrogen-to-metal (H/M) predictor for materials with unprecedented performance. We developed a metaheuristic optimizer high-throughput framework by interfacing ML models with a genetic algorithm for the accelerated search of {Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, Mo} based lightweight BCC MPEAs with improved hydrogen storage characteristics. We report five new MPEAs with a predicted gravimetric hydrogen storage capacity of around 3.5 wt% or more, including Cr_0.09Mg_0.73Ti_0.18 (4.25 wt% H) and Cr_0.21Nb_0.11Ti_0.35V_0.33 (3.5 wt% H). The electronic structure of the top-performing composition, Cr_0.09Mg_0.73Ti_0.18, was analyzed using density functional theory (DFT) to understand the reasons for its improved hydrogen storage properties compared to TiFe (1.90 wt% H), LaNi₅ (1.37 wt% H) or BCC MPEAs like TiVNbCr (3.70 wt% H). Temperature-dependent molecular dynamics (MD) studies were further performed on optimized BCC MPEAs to qualitatively study hydrogen mobility and analyze the effect of different elemental composition on bulk hydrogen diffusion. Our findings demonstrate how a ML assisted genetic algorithm framework can be used for efficient search of stable, lightweight and cost-effective MPEAs while minimizing the need for expensive ab initio calculations.

---

1. Introduction

Body-centered cubic (BCC) high-entropy alloys (HEAs), a subset of multi-principal element alloys (MPEAs), have shown strong potential for effective interstitial hydrogen storage.^1–8 HEAs are a class of alloys defined by their roughly equiatomic composition (5–35% of each element) of five or more elements. This mixture of elements of different atomic radii into a single lattice creates a significant distortion within the lattice which can potentially create an extremely high volumetric density of accessible interstitials for hydrogen occupation.^9,10 Furthermore, the high mixing entropy of several elements promotes the formation of predictable solid solution phases rather than more complex intermetallic phases.^11,12 The large number of elements and possible proportions within MPEAs also allows increased flexibility in designing materials to fit several different objective properties, which is especially important for designing effective, practical, and cost-efficient hydrogen storage materials.^13–15 However, the high dimensional space of MPEAs limits the search efficiency of conventional methods in finding an optimal composition for hydrogen storage.¹⁶ With several orders of magnitude more possible combinations of elements than traditional alloys, it is practically impossible to test every possible MPEA composition with simulation or experiment.

Machine learning (ML) methods have shown great promise in the search for better hydrogen storage materials.^17–21 With the ability to run high-throughput analysis of materials for trivial computational costs compared to ab initio methods, several ML models for predicting various features related to hydrogen storage such as weight percent or hydrogen solubility have already been developed.²² However, even with a functioning ML model capable of predicting hydrogen storage properties of hypothetical materials, discovering promising compositions for further testing is a nontrivial task. The ideal composition identified by ML then must be optimized experimentally for wt% H, hydrogen solubility, and phase stability. This results in a slow discovery rate of MPEAs for advanced energy storage materials; therefore, it remains a critical challenge in the transition to renewable energy.

Notably, a lower dimensional space problem can easily be solved using a basic grid search, but for higher dimensions like MPEAs it becomes increasingly impractical to search through every possible variation in every possible combination of elements. Other techniques, such as Bayesian optimization, are similarly unsuited for optimization with such high dimensional space.^23–25 Therefore, the development of efficient computational tools that enable the accelerated search of multi-dimensional space is essential for advancing energy storage materials.^26–29 Recently, Jennings et al. have shown that genetic algorithms (GAs), a type of metaheuristic optimizer based on Darwinian evolution, are well suited for this optimization of machine learning-based functions, being extremely robust and able to handle spaces with large numbers of dimensions.^30–32

In this work, we demonstrate the utility of combining ML models with a metaheuristic optimizer, i.e., the genetic algorithm (GA), to perform high-throughput compositional search of MPEAs for enhanced hydrogen storage. The search spans a multi-dimensional MPEA space, constrained to compositions that are thermodynamically stable in the BCC phase. The ML models, trained on literature datasets,^33–36 predict phase stability (BCC, face-centered cubic (FCC), hexagonal closed-packed (HCP), or intermetallic), bulk modulus, and hydrogen-to-metal (H/M) ratios including hydrogen solubility^16–18,37 with quantitative accuracy, achieving an average prediction error below 5% across validation datasets. We performed a systematic optimization of hydrogen storage properties on {Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, Mo} based binary, ternary, quaternary, and quinary MPEAs. Compositions optimized via the ML-GA framework are further refined using density functional theory (DFT) to analyze electronic structure contributions to hydrogen storage properties, revealing a direct correlation between d-band center shifts in transition metals and hydrogen adsorption energy. Additionally, molecular dynamics simulations quantify hydrogen diffusion coefficients at varying temperatures, demonstrating a composition-dependent activation energy barrier that influences storage kinetics.

2. Methods

2.1. Data-selection and preparation

Our framework utilizes four datasets for training four different ML models: (i) our H/M dataset consists of 431 unique compositions compiled from the literature and compositions sent to us from Sandia National Laboratories.³⁸ A subset of the H/M dataset for metal hydrides was used from the publicly available HydPark dataset.³³ All compositions with reported non-solid solution phases such as C14 Laves structures were removed. (ii) The HydPark dataset was used only for hydrogenation enthalpy prediction with 558 data points where enthalpy is distributed among 427 unique compositions of varying material classes including intermetallics, solid solutions, and Mg-alloys. (iii) A DFT-generated dataset by Mei et al. was used to predict the bulk modulus,³⁹ and an experimental dataset compiled by Huang et al. of MPEA phases was used.⁴⁰ The bulk modulus training dataset consisted mostly of 2486 refractory metals with some non-refractory element representation (except for Mg). (iv) The training dataset for phase prediction contained 1207 unique MPEAs with the vast majority of compositions being BCC or FCC. In our model, we simplified this phase classification to four categories: BCC, FCC, BCC + FCC, and IM, with all compositions with a minor intermetallic phase being classified as IM.

2.2. Featurization

We generated 44 features using the statistical values of the properties of pure elements within a composition.^41,42 For each of the 11 properties of pure elements (see Table S1), 4 features were generated: the minimum and maximum values of that property among elements within the composition, the variance of that property weighted by the proportions, and the weighted mean. The atomic packing factor, i.e., the fraction of an alloy taken up by atoms, is calculated as where is the density of the material (calculated with the rule of mixture in the given crystal phase in the training dataset), r_i and c_i are the atomic radii and proportion of element i respectively, M_i is the molar mass of element i, and N_a is Avogadro's number. These features were employed to train the bulk modulus, phase classification and hydride enthalpy models. The final three features, the probability of the composition being BCC, the probability of the composition being FCC, and the bulk modulus, were calculated using ML models and used as features only within the final H/M predictor. The bulk modulus ML model was used purely for feature creation within the H/M prediction model; the physical significance of the bulk modulus in representing the ability of the lattice to stretch to accommodate hydrogen atoms is valuable in building a more well-rounded ML model for H/M prediction.

2.3. Machine learning training

Given the small datasets (400–2000 points), we opted not to use neural networks due to their inability to predict accurately from small datasets without significant adjustment. Instead, three models were tested for each ML model: a random forest regressor (RFR), a gradient boosting regressor (GBR), and an Extreme Gradient Boosting Regressor (XGBoost). RFRs were chosen for use in H/M and bulk modulus prediction models for their high coefficient of determination, and a GBR was chosen for (de)hydrogenation enthalpy prediction for a similar reason. A gradient boosting classifier was chosen for use in the phase classification model given its similar accuracy performance to the original model used by Huang et al.⁴⁰ The equilibrium hydrogen pressures P_eq and temperature are determined using the thermodynamic parameters of the metal hydride (de)hydrogenation reaction (enthalpy change (ΔH°) and entropy change (ΔS°)), as described by the van't Hoff equation ln(P_eq/P₀) = −ΔH°/RT + ΔS°/R, where P₀ and R represent the standard atmospheric pressure and the universal gas constant, respectively. Here, enthalpy change (ΔH°) refers to the hydride reaction enthalpy and is usually closer to the desorption enthalpy than to the absorption enthalpy. In SI Fig. S1 we present the training and test loss function and other ML performance parameters (MAE, RMSE, and MAPE) with respect to the number of trained and tested datapoints and the results show that the loss function and corresponding error bars in prediction go down eventually, which suggests that the models perform well and they did not overfit. Also, hyperparameter tuning is critical in machine learning models as it controls the balance between bias and variance, and in our case, systematic tuning improved model accuracy and generalization. Hyperparameter optimization for GBR, XGBoost and RFR was carried out using grid search with 5-fold cross-validation to systematically explore combinations of model parameters. For GBR, parameters such as the learning rate, number of estimators, maximum tree depth, and subsampling ratio were tuned to balance bias and variance. In XGBoost, additional regularization parameters (gamma and alpha) and column subsampling were optimized to improve generalization and reduce overfitting. For the RFR, tuning focused on the number of trees, maximum features, and minimum samples per split/leaf, which improved model stability and predictive accuracy. Overall, careful hyperparameter selection enhanced model performance by reducing error metrics and improving robustness across training and test sets.

2.4. Uncertainty analysis

Uncertainty analysis in the RF model is rooted in its ensemble structure, where predictions are generated by aggregating outputs from N decision trees. The RF model is an ensemble of N decision trees, and for a given feature vector x, each tree T_i produces an estimate y_i(x). The RF prediction is taken to be the mean of these individual tree outputs, i.e., . To capture the RF model's uncertainty, we computed the variance of tree predictions: , while standard deviation of uncertainty in the RF model was computed using formula . This standard deviation is visualized as vertical error bars in the prediction plots, allowing a direct assessment of confidence in the model outputs across training and test datasets. Uncertainty analysis in the GBR model is performed by examining the variability of predictions across consecutive boosting stages. For each input, the final prediction is taken as the average of all stage-wise outputs, while the spread (standard deviation) of these outputs is used as a measure of predictive uncertainty. This standard deviation, interpreted as a confidence interval, was visualized as ±1σ error bars in the parity plots, providing a direct indication of the model's reliability.

2.5. ML accelerated GA

The ML-GA framework incorporates a two-step evaluation process: first, ML functions estimate a predicted fitness and second, an energy calculator (phase stability) determines the actual fitness. A nested GA is employed to explore the surrogate model, which is constructed using ML predictions. This surrogate ML model interfaced with GA serves as a high-throughput screening mechanism within the overarching “master” GA, allowing rapid exploration based solely on predicted fitness. It operates by iterating over the current population, where evaluation and selection rely on the ML-predicted wt% H, enthalpy change, and bulk modulus. To ensure robustness of the GA search and avoid ad hoc parameter selection, the fitness function—including the “k” penalty term—was systematically tuned via sensitivity analysis (SI Fig S2 and S3). Multiple values of k were tested to quantify its influence on selection pressure and convergence behavior. The analysis confirmed that the chosen value balances exploration and exploitation, providing reliable identification of promising candidates. Once the nested GA completes its search, the final population of unrelaxed candidates is returned to the master GA for further refinement.^30–32 Although reinforcement learning (RL) could theoretically provide an adaptive control mechanism for compositional optimization, its applicability is limited in this study due to the discontinuous and data-sparse nature of the objective function. The GA-ML hybrid framework presented here offers a more interpretable and computationally efficient alternative for exploring non-differentiable design spaces within limited datapoints.

2.6. Methods applied for molecular dynamics simulation of hydrogen diffusion

The interatomic potentials for the molecular simulations in this study can be represented as U = U_M–M + U_H–H + U_M–H, where U_M–M denotes the metal–metal interaction and U_H–H and U_M–H represents the hydrogen–hydrogen and metal atom–hydrogen interactions respectively. Lennard-Jones (LJ) 12–6 potential, presented in eqn (1) (ref. 43 and 44), was applied with a cut-off distance of 10 Å to model the van der Waals interactions between each of these metal elements with gaseous hydrogen atoms (U_M–H and U_H–H).where σ_ij and ∈_ij is the distance where potential energy becomes zero and the potential well depth, respectively, and r_ij is the distance between one LJ site and another. The Lorentz–Berthelot mixing rule was applied to model the interaction parameters (σ and ∈) for each of the U_H–H and U_M–H cross-interactions.⁴³ Inter- and intra-atomic LJ interaction parameters for different BCC MPEA combinations with atomic hydrogen were taken from reported literature.⁴⁵ The LJ potential was employed as a computationally efficient first-order approximation to capture qualitative hydrogen diffusion behavior in MPEAs, particularly under conditions dominated by weak, non-reactive physisorption. Although it cannot explicitly model bond formation, breaking, or electronic effects inherent to metal–hydrogen interactions, it effectively reflects how lattice geometry and atomic size mismatch influence relative trends of gaseous hydrogen diffusion across alloys over a large temperature range (−50 °C to +90 °C).²⁰ This simplified approach enables rapid screening and comparison of compositional effects before applying more accurate, but computationally demanding, reactive (embedded atom method (EAM), modified embedded atom method (MEAM) and ReaxFF) potentials in future studies.

Each simulation supercell initially consisted of 250 atoms depending on the different elemental compositions and the concentration selected, and we have added 4 hydrogen atoms at four different octahedral positions within the crystal structure to assess the hydrogen diffusion inside bulk alloys. First, the energy minimization of each structure was performed using the conjugate gradient method,⁴⁶ followed by isothermal-isobaric (NPT) equilibration at 300 K for 100 picoseconds (ps) and constant pressure heating to 600 K/900 K at 4 K ps⁻¹ and a subsequent 100 ps NPT and another 100 ps NVT equilibration at that final temperature to relax the structures. Finally, the hydrogen diffusion simulations were carried out on the relaxed structures at the final temperature for 1 nanosecond (ns). The mean squared average displacements (MSAD) and subsequent radial distribution functions (RDFs) were generated to qualitatively analyze hydrogen diffusion in different ML-GA optimized BCC MPEAs. All the simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software¹⁶ and a 1 femtosecond (fs) timestep size was used with the velocity Verlet algorithm^47,48 for performing the time integration of the equation of motions and the system temperature was controlled using the Nosé–Hoover thermostat.^45,49,50

2.7. Density-functional theory (DFT) calculations

The electronic-structure, local-lattice distortion, and hydrogen diffusion analysis was done by employing first-principles DFT as implemented in the Vienna Ab initio Simulation Package (VASP).^51,52 The generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) was employed in all calculations⁵³ with a plane-wave cut-off energy of 520 eV. The choice of PBE over LDA or meta-GGA^54,55 functionals is based on the work of Söderling et al.⁵⁶ and Giese et al.⁵⁷ that establishes the effectiveness of GGA functionals. Large Supercell Random Approximates (SCRAPs), i.e., 60 atoms per cell, with the optimized disorder (zero-correlation) were generated (a single, optimized configurational representation) for DFT calculations.⁵⁸ Energy and force convergence criteria of 10⁻⁸ eV and 10⁻⁶ eV Å⁻¹, respectively, were used for full (volume and atomic) relaxation of SCRAPs. The Monkhorst–Pack k-mesh was used for Brillouin zone integration during structural optimization and charge self-consistency calculations.⁵⁹

3. Results and discussion

Fig. 1a shows an ML assisted GA flowchart or framework for the accelerated design of MPEAs with improved hydrogen capacity. To form the starting population, a sample of 400 randomly generated equiatomic compositions is created using elements sampled from a starting list of Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, and Mo, typically used in MPEAs for H storage.⁶⁰ This sample is then ranked according to a fitness function, which is as follows: each composition with a reaction enthalpy change below the minimum threshold is scored the lowest, with higher reaction enthalpies being scored higher. Each composition with a valid enthalpy changes but an invalid phase probability is scored one tier above, with compositions having higher probability of being in the correct phase being scored higher. Finally, compositions with a valid phase and enthalpy change are scored the highest, with compositions with a higher hydrogen weight percent (calculated as the H/M predicted by ML divided by the atomic weight of the composition) scored higher. In this way, the algorithm will first optimize for valid reaction enthalpy, then for a valid phase and finally optimize for high weight percent. For each optimization, the genetic algorithm evaluated the ML models from 100–400 times each round across 10–15 rounds, for a total of 4000–6000 evaluations per optimization. This is significantly less than the total number of evaluations required for a grid search of the entire space with the same accuracy, which requires calculations in the order of billions).⁶¹ In this order, with more optimization, it may even be possible for future studies to interface ab initio calculations with genetic algorithms rather than only ML models.

Fig. 1 (a) Schematic flowchart of the ML accelerated genetic algorithm (GA) framework for the metaheuristic optimization and generation of MPEA compositions with optimized stability and improved H/M. The optimization stops when ML-assisted GA fails to produce improved hydrogen storage MPEA compositions. (b) Convergence test of the ML-GA optimizer framework on {Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, Mo} based MPEAs showing ultrafast convergence to compositions with improved hydrogen capacity.

Crossovers were not implemented in this algorithm, as they apply poorly to compositions which must maintain a total compositional sum of 100%. Rather, we created a custom mutation algorithm to mutate compositions where each parent composition has a 2/3rd probability of undergoing mutation, where 1–8% of an element's proportion is redistributed to another element, ensuring that the total remains 100%. The redistribution percentage is randomly selected from a uniform distribution between 1% and 8%. For example, if an element is at 30% and the mutation results in a 5% redistribution, it could shift to 25%, balancing another element's proportion accordingly. There is also a one-third probability that an element will be swapped by a randomly selected element from a predefined pool of elements, provided it is not already present in the composition. The selection mechanism retains the top 50 compositions based on fitness function scores.

The fitness score of a sample s, under a set of constraints c_i (e.g., BCC probability and hydrogen solubility energy) for which we require the values of those constraints of the sample v_i(s) > c_i for all i, is calculated as f(s) = v_j(s) – k × j, where j is the smallest integer such that for all i > j, v_i(s) > c_i and k is an arbitrarily large constant such that k > > max(v_i(s)) for all s. Here, we set c₀ such that c₀ > > max(v₀(s)), and v₀(s) should be set as the performance metric being sought (e.g. hydrogen weight percent). This ensures that the GA first optimizes materials to satisfy the constraints before the performance metric, reducing the need for searching in regions which fail the constraints. However, this method also risks overlooking compositions in regions that show potential for high performance but initially appear to poorly meet the constraints. This can at least in part be mitigated using a diverse initial sample. Finally, the algorithm terminates once the top 5 compositions maintain unchanged fitness scores for five consecutive generations, indicating convergence to an optimal solution. This approach effectively enhances compositional diversity while optimizing for hydrogen storage performance using quantitative metrics.

Our ML-GA optimizer in Fig. 1b demonstrates rapid convergence toward MPEA compositions with improved hydrogen storage capacity. Initially, both the average and maximum composition hydrogen storage capacities increase sharply within the first five optimization steps, indicating an efficient search of the composition space. The maximum composition (orange curve) reaches a plateau slightly earlier than the average composition (blue curve), suggesting that the optimizer quickly identifies a promising candidate with maximized hydrogen capacity. After approximately 15 optimization steps, both curves plateau, showing that further iterations yield minimal improvements. This rapid convergence highlights the GA's effectiveness in efficiently navigating the compositional space to identify MPEAs with enhanced hydrogen storage characteristics. Sensitivity analysis of the genetic algorithm to the number of children allowed to survive each round, the k parameter, was performed. We found k = 30 as an optimal value (see SI Fig. S2), minimizing unnecessary evaluations and guaranteeing that almost all trials reached the optimal value. The optimal value was found in an average of approximately 1100 evaluations (SI Fig. S3), far outperforming the number of evaluations required in a random or grid search. Our ML-GA pipeline leverages fast RF/XGBoost surrogates to rapidly screen very large surrogate-level composition sets and combines this broad search capability with targeted DFT validation, enabling diverse candidate discovery alongside focused high-fidelity evaluation. In contrast, Bayesian optimization (BO) has successfully delivered Pareto-optimal MPEA candidates with only a small fraction of the composition space explored.⁶²

Notably, Ti, V, and Mg are among the best hydrogen absorbing elements, i.e., hydrogen absorption can be maximized while maintaining favorable thermodynamics for storage {Ti, Mg} rich alloys. To visually show the optimization process, we choose three representative MPEAs, i.e., {Mg–Ti–V}, {Ti–V–Cr}, and {Ti–V–Nb}, from {Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, Mo}. We show ternary map for hydrogen storage for the three alloy classes including the starting optimization path as shown in Fig. 2. The starting point for each case was set at the center of the triangle, i.e., equiatomic MPEAs. Notably, the ML-GA optimizer always shifts away from equiatomic points toward Ti/V-rich or Mg-rich formulations. The Nb–Ti–V system exhibits a peak hydrogen storage capacity of 3.6 wt% in the Ti-rich region (Fig. 2a), aligning with theoretical capacities of 4.0 wt% H for TiH₂ and 3.8 wt% H for VH₂. The lower capacity of NbH (≈2.6 wt% H) accounts for the Nb diminished role, while in the Mg–Ti–V system in Fig. 2b, the highest storage reaches 3.9 wt% near Mg ∼60–70 at%, reflecting MgH₂'s high capacity (7.6 wt% H). Although the model does not explicitly account for desorption and kinetics, it necessitates the inclusion of Ti and V into the final alloy composition, which are known to improve rates of hydrogen uptake and release. While the Cr–Ti–V system in Fig. 2c shows a peak storage of 3.8 wt% H, Cr's limited hydrogen solubility (∼0.4 wt% H) leads the MLGA optimizer to favor Ti–V compositions similar to those in the Nb–Ti–V system. Comparatively, Mg–Ti–V offers the highest storage capacity but requires destabilization for practical use, while Nb–Ti–V and Cr–Ti–V achieve slightly lower gravimetric capacities but offer better kinetic properties.⁶³ Notably, across all three systems, Ti and V play critical roles in optimizing hydrogen uptake, reinforcing their importance in hydrogen storage alloy design. This also effectively justifies the reason the ML-GA optimizer picks the elements favorable for hydrogen absorption.

Fig. 2 ML-GA driven hydrogen storage optimization in (a) Nb–Ti–V, (b) Mg–Ti–V, and (c) Cr–Ti–V systems. Contours show hydrogen storage capacity (wt% H), with red regions indicating higher storage. White stars mark optimal paths of compositions, favoring Ti-rich and Mg-rich regions. The Mg–Ti–V system achieves the highest capacity (∼3.9 wt% H), while Nb–Ti–V and Cr–Ti–V peak around 3.6–3.8 wt% H, highlighting Ti and V as key enhancers.

Notably, the strength of the metaheuristic ML-GA optimizer used for searching new alloying compositions over multi-dimensional space as shown in Fig. 2 is based on the foundations of best-performing ML models that are used as drivers for providing seeds for searching new MPEAs with hydrogen solubility, phase stability, and high storage. Fig. 3 shows three ML models, each interfaced with a metaheuristic optimizer to identify key descriptors. The H/M ratio model in Fig. 3b shows strong accuracy (R² = 0.828, MAE = 0.179, and RMSE = 0.253), with Mulliken electronegativity, melting point, thermal expansion, and bulk modulus as dominant features (Fig. 3c). The associated error bars indicate moderate predictive uncertainty, which increases at higher H/M values, reflecting reduced model confidence in extrapolated regimes. The observed variation in the error bars arises because they represent ensemble-predicted uncertainty, reflecting the model's internal confidence rather than the true deviation from actual values. Regions in the feature space with dense data points yield smaller uncertainty estimates, as the predictions from ensemble models (RF and XGBoost) are more confident there, whereas sparsely populated regions naturally exhibit larger uncertainties. Thus, low uncertainty does not necessarily imply perfect prediction accuracy, and vice versa.

Fig. 3 Training and test results of the best-performing machine-learning models, along with feature importance plots, for (a–c) H/M, (d–f) reaction enthalpy, and (g–i) bulk modulus.

The reaction enthalpy model in Fig. 3e shows moderate correlation (R² = 0.608, MAE = 9.781 kJ mol⁻¹, and RMSE = 17.411 kJ mol⁻¹), with electronegativity, electron affinity, and covalent and atomic radii playing key roles (Fig. 3f). However, this model exhibits noticeably larger uncertainty bands, particularly for higher enthalpy values, suggesting limited generalization beyond the training domain. The bulk modulus model in Fig. 3h demonstrates good agreement (R² = 0.800, MAE = 8.397 GPa, and RMSE = 11.325 GPa), primarily influenced by electronegativity, covalent radius, thermal expansion, and mean density (Fig. 3i). Compared to H/M and reaction enthalpy, the bulk modulus predictions show narrower error bars and more stable uncertainties, indicating greater robustness. Across all models, electronegativity and atomic size emerge as critical descriptors of property trends. Overall, the models performed adequately in terms of 5-fold R² cross-validation (CV), and the uncertainty analysis highlights property-dependent generalization behavior: bulk modulus predictions are the most reliable and H/M shows moderate robustness, while reaction enthalpy exhibits higher variance and weaker generalization. Cross-validation was carefully carried out to account for potential data imbalance, particularly for materials with very low or high ΔH values.⁶⁴ The bulk modulus and phase classification models performed similarly to what was reported by Mei et al.³⁷ and Huang et al.,⁴⁰ respectively, further validating the suitability of the chosen feature set and model architecture. We also tested our phase predictor and reaction enthalpy predictor models within the updated framework incorporating additional featurization and selection strategies (SI Fig. S4), which led to a slight improvement in model performance (see Table S3), particularly with LightGBM⁶⁵ over random forest. However, this enhancement was marginal and did not qualitatively alter the predictions or the outcomes of our GA-based framework.

We have also performed SHAP analysis for an interpretable or explainable understanding of the most important features shown in Fig. 3 for all three models (SI Fig. S5). SHAP analysis agrees quite well with results in Fig. 3, and it also shows that bulk modulus predictions in hydrogen storage alloys are primarily controlled by electronegativity descriptors, atomic radii, and density. Higher electronegativity and density increase stiffness, while larger atomic radii and higher thermal expansion reduce it. On the other hand, higher values of Mulliken electronegativity (blue points on the right of SI Fig. S5a) reduce predicted H/M, whereas lower values (pink/red points) increase H/M. For the reaction enthalpy predictive model, SHAP analysis reveals that the predictions are dominated by the covalent radius, electronegativity descriptors, and density-related features. Larger covalent radii and lower electronegativity values decrease predicted stability, while higher electronegativity and structural variability (e.g., chemistry, lattice or volume, crystal structure, etc.) increase the enthalpy. These findings highlight the combined role of electronic bonding and lattice heterogeneity in governing hydrogen reaction thermodynamics (SI Fig. S5).

Notably, the new H/M ML model demonstrates a 5-fold CV test score of 0.79. Although dedicated models for H/M prediction are largely absent in the literature, several studies^66,67 have already estimated hydrogen weight percent, which can be directly converted to H/M ratios using a straightforward mathematical transformation. This allows our ML model to effectively bridge the gap by leveraging existing hydrogen storage data to predict H/M with high accuracy. Our model achieves a testing R² of 0.89 for wt% H, which is a notable improvement from existing literature for weight percent prediction in metal hydrides.²² However, it is also important to note that these values do not necessarily accurately represent the true capabilities of the model. The training and testing set for this model is limited in scope with a low diversity of data points, with some compositions within the data set being almost identical apart from a small deviation in elemental composition. Many of the high-importance features in the models in Fig. 3c also align with key physical parameters. Hydrogenation enthalpy was dominated by the atomic and electronic parameters like the covalent radius, electronegativity and electron affinity, while phase classification (Fig. 4) had a strong dependence on the packing fraction (which is known to be lower in BCC structures). The cause behind the importance of the mean melting point in H/M prediction is unclear, but it may be possible that materials with lower mean melting points have a greater ability to expand their lattices to accommodate more hydrogen.

Fig. 4 (a) The confusion matrix for phase classification, highlighting the model's predictive accuracy across different phases, and (b) the top 10 most important features influencing the classification decision. The RF classifier is shown in SI Fig. S6.

In Fig. 4, we present an optimized phase classifier for analyzing the confusion matrix that ensures reliable predictions of phase stability, aiding in the computational screening of alloys with targeted properties.⁶⁴ Each row in Fig. 4a shows the true phase label, and each column corresponds to the predicted phase label. The diagonal elements indicate the correct predictions, while the off-diagonal elements reflect misclassifications. The model performs well in predicting FCC (85% accuracy) and intermediate (IM) phases (74%), whereas the BCC phase has a moderate accuracy of 79%, with some misclassification into IM (13%) and FCC + BCC (5.7%). The FCC + BCC phase has the lowest accuracy (59%) due to increased overlap with other phases, particularly BCC and IM. Feature importance ranking derived from the trained model is shown in Fig. 4b. The most critical predictor is the minimum covalent radius, followed by the variance in the covalent radius and variance in the atomic number. These features play a significant role in determining phase stability due to their influence on atomic packing and electronic interactions. The packing fraction, melting point, and variance of the lattice constant are also important, reinforcing the role of structural and thermodynamic properties in phase formation. Entropic descriptors like ΔS_mix contribute less but still provide valuable information on phase stability in MPEAs. The results emphasize the dominant role of atomic-scale features in phase classification and provide guidance for material design strategies.

3.1. Model performance analysis

The performance of different ML models for predicting H/M ratios, phase classification (confusion matrix), reaction enthalpy, and bulk modulus is shown in Table 1. Interestingly, the phase classification model achieves an accuracy of 0.76 on testing and 0.96 on training, with a 5-fold cross-validation score of 0.77. While training performance is very high, the testing score suggests slight overfitting, but the model still performs well in classifying alloy phases. Overall, all other models exhibit strong predictive performance and remain robust across all validation tests.

Table 1 A summary of model performance across four models, i.e., H/M, phase classification (confusion matrix), reaction enthalpy, and bulk moduli

Model	Testing RMSE	Training RMSE	Testing R^2	Training R^2	5-Fold CV score
a Accuracy score.
H/M	0.25	0.13	0.83	0.95	0.79
Phase classification	N/A	N/A	0.76^a	0.96^a	0.77^a
Reaction enthalpy	17.4 kJ mol^−1	4.6 kJ mol^−1	0.61	0.97	0.71
Bulk modulus	11.3 GPa	10.3 GPa	0.80	0.83	0.81

---

The performance test of H/M model in Table 1 demonstrates excellent generalization with a low testing RMSE of 0.25 and a training RMSE of 0.13, suggesting that it effectively captures the underlying trends in the data. The high R² score (0.95 for training and 0.83 for testing) indicates strong predictive power, while the 5-fold cross-validation score (0.79) confirms model robustness. For reaction enthalpy predictions, the model achieves a testing RMSE of 17.4 kJ mol⁻¹ and a training RMSE of 4.6 kJ mol⁻¹. The testing R² score (0.61) indicates a reasonably strong correlation between predictions and actual values, while the training R² (0.97) is higher, suggesting that the model may slightly overfit. Predictions may critically depend on other parameters (mostly thermodynamic parameters like temperature and pressure) that were not considered features to remove bias, and the data points for reaction enthalpy >80 kJ mol⁻¹ H are quite sparse for better predictions (Fig. 3e). The 5-fold cross-validation score (0.73) confirms its reliability. The bulk modulus model achieves a testing RMSE of 11.3 GPa and a training RMSE of 10.3 GPa, with testing and training R² values of 0.80 and 0.83, respectively. The 5-fold CV score of 0.81 indicates consistency and reliability across different data splits. All these suggest well-balanced models with good predictive ability.

3.2. Validation of the ML-GA framework

Due to the lack of experimental data on wt% H for specified ternary or beyond MPEA space in the literature, the ML-GA was tested on binary Mg–Ni systems to test the framework's ability to accurately locate weight percent maxima within a composition space. Fig. 5a shows the relationship between the Mg fraction and hydrogen storage capacity (weight percent), highlighting both experimental and predicted results.⁶⁸ The x-axis represents the Mg fraction in the alloy, while the y-axis denotes the hydrogen storage capacity in weight percent. Different markers distinguish between experimental non-BCC (blue circles), experimental BCC (orange circles), predicted maximum capacity (green triangle), and the predicted maximum for BCC structures with a confidence level greater than 0.95 (red triangle). This means that our ML-GA optimizer is able to accurately predict the hydrogen storage capacity (in terms of wt% H) in binary Mg–Ni alloys with respect to the experimentally available data from the literature.⁶⁹ We found that Mg-rich alloys show higher hydrogen storage capacity compared to the one observed in Mg-poor regions, reaching its peak around 80–90% Mg. The experimental results show that BCC structures tend to have higher hydrogen weight percent⁶⁹ compared to non-BCC structures, which agrees with ML-GA predictions. The green and red triangles indicate ML-predicted optimal compositions, suggesting that specific Mg–Ni ratios can achieve even higher hydrogen storage capacities. Notably, the genetic algorithm was able to locate the maximum weight percent for both BCC and non-BCC alloys to within a 5% accuracy in Mg fractions. Since large deviations from equiatomic compositions often stabilize ordered or FCC/intermetallic phases that alter hydrogen-storage mechanisms and preclude direct comparison with the BCC experimental data, restricting sampling to near-equiatomic ranges (±10 at% around equiatomic) provides the most relevant, reproducible, and experimentally comparable domain for ML-GA optimization and validation.

Fig. 5 (a) The validation test of hydrogen storage capacity was done on binary Mg–Ni systems, comparing experimental and ML-predicted maxima for BCC and non-BCC phases. Optimizer predicted wt% H for the Mg–Ni system shows excellent agreement with experiments.^47,66 (b) Pareto-front analysis of wt% H vs. hydrogenation enthalpy for Cr–Nb–Ti–V MPEAs (color represents the time evolution of the GA search).

Fig. 5b illustrates the Pareto-front relationship between two objectives, i.e., hydrogen weight percent (wt% H) and reaction enthalpy, for Cr–Nb–Ti–V MPEAs that defines the trade-off between phase stability (minimizing) and hydrogen storage capacity (maximizing). For the given set of feasible solutions, i.e., best MPEA compositions, each composition satisfies the above two objectives. The colors indicate the time evolution of the optimization process with red color representing later stages of compositional evolution, while blue points correspond to the initial stage of the optimization process. Notably, there exists a strong positive correlation between hydrogen capacity and desorption enthalpy, where increasing hydrogen weight percent is generally accompanied by an increase in reaction enthalpy, while lower reaction enthalpy often corresponds to a decrease in total hydrogen capacity. Over time, the evolutionary algorithm tends to favor compositions with higher hydrogen storage capacity, leading to clustering of points at higher wt% H and reaction enthalpy. The dense cluster on the top right in Fig. 5b highlights this trade-off, where increased hydrogen capacity is often associated with higher hydrogenation enthalpy, suggesting a need for further refinement to balance storage performance and thermodynamic stability. Ultimately, the plot confirms that the ML-GA framework effectively guides the search toward promising hydrogen storage materials while indicating potential areas for improvement in optimization constraints.

3.3. Metaheuristic optimizer predicted MPEAs with improved hydrogen storage capacity

The ML-GA optimized search results for the five best {Mg–Ti–V–Cr–Fe–Nb–Mo} MPEAs are shown in Table 2 (the extended table with the top eight MPEA compositions is given in supplementary Table S2). Ni and Co were deprioritized because in many Mg- or Ti-rich BCC MPEA compositions they promote Ni-rich intermetallics/ordered phases or FCC phases that change hydrogen uptake pathways, while Al was excluded for causing phase separation or brittle intermetallics in pre-screening. In this work, the alloy search scheme is designed for hydrogen storage, focusing on maximizing wt% H, while ensuring stability in the BCC phase as shown by Pareto-front analysis in Fig. 5b. Each of the five best compositions in Table 2 consists of multiple elements with varying atomic fractions, and the table includes estimates for hydrogen storage capacity, thermodynamic stability, and phase probability under specific constraints. The estimated hydrogen weight percent across the composition range is from 3.01% to 4.22%, with the highest value observed for Cr_0.09Mg_0.73Ti_0.18, a ternary alloy with a high magnesium content (see Table S2). Most other compositions exhibit hydrogen storage capacities between 3.01% and 3.45%, indicating a relatively narrow range of optimal solutions identified by the genetic algorithm. The estimated ΔH spans from −61.6 kJ mol⁻¹ to −38.3 kJ mol⁻¹, with Cr_0.11Fe_0.05Mg_0.32Ti_0.35V_0.17 exhibiting the most negative value at −61.6 kJ mol⁻¹, indicating the strongest hydrogen–metal interactions. In contrast, alloys with less negative enthalpy values, such as Cr_0.21Nb_0.11Ti_0.35V_0.33 at −38.3 kJ mol⁻¹ H, may offer improved or similar hydrogen desorption compared to LaNi₅ (−40 ± 7 kJ mol⁻¹ of H)⁷⁰ or BCC MPEAs with unfavorable desorption (50–70 kJ mol⁻¹ H).^71–73 The BCC phase probability (pBCC) varies between 85% and 98%, with the composition Cr_0.11Mg_0.15Mo_0.15Ti_0.29V_0.30 achieving the highest probability at 98% (Table S2). Most compositions exceed 90% BCC phase probability, reinforcing the preference for this structure in hydrogen storage applications.

Table 2 Rank ordering of five best Mg–Ti–V–Cr–Fe–Nb–Mo MPEA compositions (see Table S2) for hydrogen storage predicted by the ML-GA optimizer, displaying estimated hydrogen weight percent, formation enthalpy, BCC phase probability, and applied constraints. MD calculated mean squared average displacement (Ų) of hydrogen diffusion is provided for comparison. The results highlight BCC MPEAs with high H storage capacity

Rank	Composition	Estimated wt% H	ΔH [meV per datom]	% BCC phase formation (pBCC)	Constraints	H diffusion, MSAD [Å^2] [300 K/600 K]
MPEA1	Cr0.09Mg0.73Ti0.18	4.22	−60.7	91	|ΔH| > 40 pBCC > 0.8	0.28/0.81
					cmin = 1
					cmax = 0.0
					3 elements
MPEA2	Cr0.21Nb0.11Ti0.35V0.33	3.45	−38.3	91	|ΔH| > 30 pBCC > 0.7	0.02/0.04
					cmin = 0.05
					cmax = 0.35
					4 elements
MPEA3	Cr0.16Mg0.09Nb0.08Ti0.35V0.32	3.37	−54.5	85	|ΔH| > 40 pBCC > 0.8	0.04/0.63
					cmin = 0.05
					cmax = 0.35
					5 elements
MPEA4	Cr0.25Mg0.35Ti0.35V0.05	3.24	−59.7	95	|ΔH| > 40 pBCC > 0.8	0.09/0.12
					cmin = 0.05
					cmax = 0.35
					4 elements
MPEA5	Cr0.11Fe0.05Mg0.32Ti0.35V0.17	3.18	−61.6	96	|ΔH| > 60 pBCC > 0.8	0.12/0.23
					cmin = 0.05
					cmax = 0.35
					5 elements

---

The ML-GA optimizer predicted compositions in Table 2 are all >75% either Mg, Ti, or V. All three of these elements are both light weight and have a relatively high affinity for hydrogen, suggesting that the model is capable of effectively selecting elements with strong hydrogen storage potential. The model also incorporated several elements with high BCC stability, such as Mo, Cr, and Fe in order to fit the minimum required BCC probability. For the composition with the highest probabilistic score of BCC phase formation, [Cr_0.11Mg_0.15Mo_0.15Ti_0.29V_0.30 (Table S2)] the model incorporated molybdenum despite its low hydrogen affinity in order to increase the probability of forming a BCC alloy. This behavior has been experimentally verified in BCC titanium alloys, where it is known that molybdenum can make BCC phases more thermodynamically stable.⁷⁴ All the compositions in Table 2 had very high (>3) weight percentages for hydrogen storage, with Cr_0.09Mg_0.73Ti_0.18 being an outlier at an unprecedented 4.22 weight percent. However, given that it is not a high entropy alloy, it is likely that our phase prediction model may not be able to accurately predict its phase. The DFT results show that this composition is more energetically stable in an HCP phase (due to the high percent of magnesium) than a BCC phase, although not by a significant margin. It is conceivable that a similar composition could easily be designed to have a more favorable BCC phase while still maintaining a very high weight percentage. The 0 K DFT result showing a slight HCP preference for Cr_0.09Mg_0.73Ti_0.18 can be offset at finite temperature by configurational entropy (19.5 meV per atom at 300 K and 65 meV per atom at 1000 K) and vibrational free-energy contributions (10–50 meV per atom), both favoring BCC. Under hydrogenation, interstitial entropy and strain relief further stabilize BCC, making its appearance in experiments and ML predictions reasonable despite the small 0 K enthalpy gap.

Notably, the vanadium-based BCC alloys with a higher concentration (≥30 at%) have demonstrated excellent hydrogen absorption capacities, reaching approximately 3.5 wt% H.^75–78 Thus, to develop compositionally versatile hydrogen storage materials, multicomponent BCC alloys incorporating elements such as Mg, Ti, V, Cr, Fe, Nb, and Mo have been explored.³³ These elements contribute to structural stability while optimizing hydrogen absorption and desorption properties. To stabilize the BCC phase at room temperature, elements such as V and Cr are commonly used as the phase stabilizers.^78–80 Cr and Fe enhance mechanical strength, Nb and Mo improve phase stability, and Mg facilitates hydrogen diffusion, making such alloys promising candidates for hydrogen storage applications. The ML-GA predicted compositions in Table 2 show good trends when compared to existing work on Cr–Ti–V based alloys. For example, Okada et al.⁷⁵ also reported higher hydrogen weight percent (3 wt% H) for a Ti_0.27V_0.40Cr_0.27Fe_0.06 MPEA. Additionally, a classical molecular dynamics (MD) study evaluated hydrogen mobility in the five MPEAs, with results in Table 2 showing MPEA1 exhibiting the highest mean squared average displacement (MSAD) at room temperature (300 K), aligning with ML-GA predictions.

3.4. A quantitative analysis establishing the relationship between hydrogen diffusion rates and alloy composition

To further substantiate, we expand our analysis to present a quantitative composition–diffusivity relationship among the five MPEAs in Table 2. The temperature dependent diffusion rates obtained from MD simulations are shown in Fig. 6a and b, which suggest a clear, consistent pattern, i.e., hydrogen is the most mobile species in every alloy and its mobility grows far more with temperature than the metallic constituents. From the plotted bars at 300 K hydrogen MSADs are visually in the order of a few 10⁻¹ Ų (for example MPEA1 ≈0.27 Ų versus typical metal species ≈0.12–0.17 Ų), whereas MPEA2 exhibits very low MSADs for both metals and hydrogen (≈0.01–0.04 Ų), indicating a much more constrained environment. At 600 K the hydrogen MSADs increase dramatically (MPEA1 near ∼0.8 Ų and MPEA3 ∼0.6 Ų by visual estimate) while most metal MSADs increase by smaller factors; hydrogen therefore often shows a roughly threefold or larger temperature-driven increase compared with the roughly twofold increases seen for many metal species. Element-specific differences are also evident: certain atomic species in MPEA1 and MPEA3 show disproportionately large increases at 600 K, suggesting composition- or short-range-order-dependent activation of diffusion for particular elements. Physically, the data are consistent with hydrogen behaving as a fast interstitial diffuser with relatively low barriers in MPEA1 and MPEA3 and with deep traps or steric constraints in MPEA2 that suppress both lattice and hydrogen motion. This interpretation has direct implications for materials selection: alloys like MPEA1 and MPEA3 would permit rapid hydrogen transport (advantageous for permeation or fast uptake) but may be more susceptible to hydrogen-induced damage unless traps are engineered, while MPEA2-like chemistries could be preferable when immobilizing hydrogen is desired. To quantify these implications one should convert the reported MSADs to diffusion coefficients using the Einstein relation (MSAD = 6Dt for 3D diffusion) with the simulation time window and if possible obtain D(T) across multiple temperatures to extract activation energies from an Arrhenius fit; complementary analyses such as van Hove functions, residence-time statistics, and site-specific binding energies will then distinguish true long-range diffusion from local hopping or trapping and complete the mechanistic picture. CALPHAD-based thermodynamic phase stability analyses (see SI Fig. S7) and atomistic simulations confirm that no structural or phase transitions occur within the analyzed temperature range, ensuring that the observed differences in hydrogen MSAD arise from compositional variation and diffusion phenomena rather than phase instabilities. Correspondingly, Arrhenius fitting of temperature-dependent diffusivities yields lower activation energies for MPEA1 and MPEA3 (0.18 to 0.25 eV) relative to MPEA2 (0.35 to 0.40 eV), making them low-barrier hydrogen conductors. Thus, the stronger temperature dependence observed in MPEA2 stems from its deeper trapping sites and higher activation barrier, while MPEA1 and MPEA3 maintain intrinsically higher diffusivity, consistent with rapid interstitial diffusion facilitated by local short-range-order variations.

Fig. 6 MD calculated mean squared average displacements (MSAD in Ų) for each element in the presence of hydrogen for each MPEA at (a) 300 K, and (b) 600 K. The corresponding radial distribution function is shown in SI Fig. S8.

3.5. Understanding electronic-structure origin of high gravimetric capacity

We employed DFT to analyze the role of electronic structure on hydrogen storage properties of top performing MPEAs (see method description in the SI for more details). Fig. 7a illustrates the relationship between DFT-calculated volume (V) and ML-predicted wt% H, which shows that larger specific volumes tend to accommodate more hydrogen, leading to higher weight percentages. Our predictions show Cr_0.09Mg_0.73Ti_0.18 as the composition with the highest hydrogen storage capacity, while Cr_0.11Mg_0.15Mo_0.15Ti_0.29V_0.30 is positioned at a lower wt% H despite its relatively small volume. In Fig. 7b, we show correlation between phase stability, the ratio of the bulk modulus (B₀) and E_form in the presence of hydrogen, i.e., |B₀/E_{form, H2}|, which indicates that hydrogen absorption will be easier for alloys with moderate stability and higher |B₀/E_{form, H2}|. Notably, a lower |B₀/E_{form, H2}| suggests very high stability, making hydrogen absorption difficult. A strong positive trend in Fig. 7b also supports this hypothesis. Cr_0.09Mg_0.73Ti_0.18 is positioned at the upper end of the trend, making it a promising candidate for hydrogen storage. In contrast, Cr_0.11Mg_0.15Mo_0.15Ti_0.29V_0.30 is situated at a relatively lower formation energy, suggesting reduced stability.

Fig. 7 (a) Correlation between DFT calculated volume vs. ML-GA predicted wt% H for in Mg–Ti–V–Cr–Fe–Nb–Mo MPEA. (b) The formation energy (E_form,ML) as a function of the bulk modulus-to-formation energy ratio (|B₀/E_{form, H2}|), highlighting the stability and mechanical properties of the alloys. (c) The mean squared average displacement (Ų) shows 2.5x (600 K) and 5x (900 K) increases in hydrogen diffusion in Cr_0.09Mg_0.73Ti_0.18 compared to room temperature (300 K). (d) Total DOS of BCC Cr_0.09Mg_0.73Ti_0.18 pristine (shaded), with H in the octahedral site (black) and with H in the tetrahedral site (orange), and (e) zoomed TDOS near the Fermi-level (E_Fermi). The vertical dashed line marks E_Fermi.

In Fig. 7c, we show MD calculated MSAD for hydrogen diffusion for the Cr_0.09Mg_0.73Ti_0.18 MPEA that shows the highest hydrogen diffusivity at room temperature compared to the other alloys (Table 2). The higher ML-GA predicted hydrogen storage correlated with higher MSAD (see MSAD of individual elements in SI Fig. S9). In Fig. 6a and b, the MD calculated MSAD (Ų) values for all 5 MPEAs in Table 2 are given for 300 K and 600 K showing the temperature dependency for hydrogen diffusion in MPEAs. While Cr_0.09Mg_0.73Ti_0.18 (MPEA1) exhibits the highest hydrogen diffusivity at room temperature, Cr_0.16Mg_0.09Nb_0.08Ti_0.35V_0.32 (MPEA3) shows an approximately 15x increase in hydrogen diffusivity at 600 K compared to room temperature (also shown in Table 2). Notably, the pair distribution function (see SI Fig. S8a and b) for all ranked MPEAs in Table 2 shows structural integrity in the BCC phase (good agreement with DFT (Fig. 7b) and ML predictions (Fig. 4a)) despite prolonged hydrogen exposure, i.e., >one nanosecond (on the MD scale). Interestingly, Cr_0.09Mg_0.73Ti_0.18 shows nearly a 3x increase in hydrogen diffusivity at 600 K compared to room temperature (300 K). This strong temperature dependence highlights the influence of alloy composition on the activation of diffusion pathways. The overall temperature-dependent diffusion simulations reveal that the mean squared displacement of hydrogen atoms inside the alloy crystals increases dramatically, rising by approximately 200% at 600 K and further surging to 400% at 900 K. These results underscore the crucial role of temperature in modulating hydrogen mobility within MPEAs, with significant implications for optimizing hydrogen storage materials under varying operational conditions.

Finally, the electronic structure of BCC Cr_0.09Mg_0.73Ti_0.18 was analyzed through the total density of states (DOS) shown in Fig. 7d and e for three configurations, i.e., pristine, hydrogen at octahedral sites, and hydrogen at tetrahedral sites. The pristine alloy displays a broad, feature-rich DOS with finite states at the Fermi level (E_Fermi), reflecting metallic stability and compositional disorder.^20,81 When hydrogen occupies the octahedral site, a distinct H-1s feature appears deep below the E_Fermi with only weak overlap with metal d-states, indicating localized bonding and minimal perturbation to the metallic framework, conditions favorable for reversible hydrogen uptake. In contrast, tetrahedral hydrogen introduces sharper DOS peaks closer to E_Fermi, evidencing stronger H-1s-d hybridization and enhanced binding strength but potentially slower desorption kinetics. The hydrogen-projected DOS in Fig. 7e confirms this site dependence (see SI Fig. S10), with tetrahedral H contributing more near E_Fermi. Retention of metallicity in the octahedral configuration suggests reduced phase instability and an optimal balance between absorption and desorption. Thus, the tunable electronic structure and controlled hybridization in Cr_0.09Mg_0.73Ti_0.18 underpins its promise for efficient, reversible hydrogen storage.

4. Conclusion

In summary, this work successfully demonstrates the potential of applying an ML-GA framework to predict and optimize the hydrogen storage properties of BCC MPEA hydrides. By creating a novel H/M predictor, we were able to accelerate the identification of promising hydrogen storage alloys based on {Mg, Al, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Nb, Mo} BCC MPEAs showing high (>3 wt%) gravimetric hydrogen storage capacities. The integration of ML with a genetic algorithm-based optimizer further enhanced the search process, resulting in the prediction of several MPEAs with high hydrogen capacity and tunable thermodynamics. The ML-GA framework successfully identifies compositions that optimize both hydrogen storage capacity and phase stability, highlighting the trade-off between formation enthalpy and wt% H as a critical factor in material selection. Additionally, we also discuss the role of the electronic structure using DFT considering both octahedral and tetrahedral site occupation by hydrogen and highlighting the origin of higher storage capacity in Cr_0.09Mg_0.73Ti_0.18. MS simulations show that Cr_0.09Mg_0.73Ti_0.18 exhibits the highest hydrogen diffusivity at room temperature among the optimized BCC MPEAs, with H diffusion increasing nearly threefold at 600 K. Temperature-dependent simulations reveal a dramatic increase in hydrogen mobility, emphasizing the critical role of alloy composition and temperature in hydrogen storage properties. This work highlights the powerful synergy that arises from combining ML methodologies and evolutionary algorithms with atomistic simulations, allowing efficient exploration and identification of stable MPEAs with high hydrogen capacity, reducing reliance on costly ab initio calculations and paving the way for more efficient hydrogen storage materials. This framework could also be adapted to accelerate the discovery of materials for other applications, such as predicting hydrogen diffusion and permeability in metal alloys to enable facile H separation and purification.⁸²

Conflicts of interest

There are not conflicts to declare.

Data availability

Data are available upon request from the authors on reasonable request. The codes used in this work are hosted on GitHub (linked): (1) https://github.com/kevjidji/Genetic-Material-Framework, (2) https://github.com/Tanumoybanerjee88/ML-GA-prediction-for-BCC-alloy-for-H-storage.

Supplementary information: details of the density functional theory methods, feature engineering protocols, AI/ML validation procedures, an extended table of predicted hydrogen-storage MPEAs, supporting tests of the ML-GA framework, CALPHAD-based thermodynamic results, and molecular dynamics data on pair distribution functions. See DOI: https://doi.org/10.1039/d5ta06903c.

Acknowledgements

This work was supported by the Laboratory Directed Research and Development (LDRD) program at Ames National Laboratory. Kevin Ji is grateful for the research opportunity at Ames National Laboratory supported by the U.S. Department of Energy (DOE), Office of Science, Science Undergraduate Laboratory Internships (SULI) program. The Ames National Laboratory is supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. The Ames National Laboratory is operated for the U.S. DOE by Iowa State University under contract DE-AC02-07CH11358. We thank Dr Kavita Joshi of CSIR, NCL Pune (India), for sharing literature data on hydrogen storage materials. The authors gratefully acknowledge research support from the U.S. Department of Energy (DOE), Office of Energy Efficiency and Renewable Energy, Fuel Cell Technologies Office, through the Hydrogen Storage Materials Advanced Research Consortium (HyMARC). Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration (DOE/NNSA) under contract DE-NA0003525. This written work was authored by an employee of NTESS. The employee, not NTESS, owns the right, title and interest in and to the written work and is responsible for its contents. Any subjective views or opinions that might be expressed in the written work do not necessarily represent the views of the U.S. government. The publisher acknowledges that the U.S. government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this written work or allow others to do so, for U.S. government purposes. The DOE will provide public access to results of federally sponsored research in accordance with the DOE Public Access Plan.

References

1. Y. Jiang and W. Jiang, High Entropy Alloys: Emerging Materials for Advanced Hydrogen Storage, Energy Technol., 2024, 12, 2401061.
2. L. Kong, et al., A review on BCC-structured high-entropy alloys for hydrogen storage, Front. Mater., 2023, 10, 1135864.
3. G. Ek, et al., Elucidating the Effects of the Composition on Hydrogen Sorption in TiVZrNbHf-Based High-Entropy Alloys, Inorg. Chem., 2021, 60, 1124–1132.
4. R. B. Strozi, D. R. Leiva, J. Huot, W. J. Botta and G. Zepon, An approach to design single BCC Mg-containing high entropy alloys for hydrogen storage applications, Int. J. Hydrogen Energy, 2021, 46, 25555–25561.
5. A. Bouzidi, L. Laversenne, V. Nassif, E. Elkaim and C. Zlotea, Hydrogen Storage Properties of a New Ti-V-Cr-Zr-Nb High Entropy Alloy, Hydrogen, 2022, 3, 270–284.
6. A. H. Silva, C. Zlotea, Y. Champion, W. J. Botta and G. Zepon, Design of TiVNb-(Cr, Ni or Co) multicomponent alloys with the same valence electron concentration for hydrogen storage, J. Alloys Compd., 2021, 865, 158767.
7. M. M. Nygård, et al., The average and local structure of TiVCrNbDx (x=0,2.2,8) from total scattering and neutron spectroscopy, Acta Mater., 2021, 205, 116496.
8. F. Marques, et al., Mg-containing multi-principal element alloys for hydrogen storage: A study of the MgTiNbCr0.5Mn0.5N0.5 and Mg0.68TiNbNi0.55 compositions, Int. J. Hydrogen Energy, 2020, 45, 19539–19552.
9. E. Halpren, X. Yao, Z. W. Chen and C. V. Singh, Machine learning assisted design of BCC high entropy alloys for room temperature hydrogen storage, Acta Mater., 2024, 270, 119841.
10. P. Singh, B. Vella, G. Ouyang, N. Argibay, G. Cui, R. Arroyave and D. D. Johnson, A ductility metric for refractory-based multi-principal element alloys, Acta Mater., 2023, 257, 119104.
11. K. T. Møller, T. R. Jensen, E. Akiba and H.-W. Li, Hydrogen - A sustainable energy carrier, Prog. Nat. Sci. Mater. Int., 2017, 27, 34–40.
12. Q. Hou, X. Yang and J. Zhang, Review on Hydrogen Storage Performance of MgH2: Development and Trends, ChemistrySelect, 2021, 6, 1589–1606.
13. D. B. Miracle and O. N. Senkov, A critical review of high entropy alloys and related concepts, Acta Mater., 2017, 122, 448–511.
14. P. Singh, A. Sharma, A. V. Smirnov, M. S. Diallo, P. K. Ray, B. Balasubramanian and D. D. Johnson, Design of high-strength refractory complex solid-solution alloys, npj Comput. Mater., 2018, 4(1), 16.
15. E. P. George, D. Raabe and R. O. Ritchie, High-entropy alloys, Nat. Rev. Mater., 2019, 4, 515–534.
16. M. Witman, et al., Data-driven discovery and synthesis of high entropy alloy hydrides with targeted thermodynamic stability, Chem. Mater., 2021, 33(11), 4067–4076.
17. M. D. Witman, et al., Extracting an empirical intermetallic hydride design principle from limited data via interpretable machine learning, J. Phys. Chem. Lett., 2019, 11(1), 40–47.
18. M. D. Witman, et al., Towards Pareto optimal high entropy hydrides via data-driven materials discovery, J. Mater. Chem. A, 2023, 11(29), 15878–15888.
19. M. D. Allendorf, et al., Challenges to developing materials for the transport and storage of hydrogen, Nat. Chem., 2022, 14(11), 1214–122.
20. T. Banerjee, K. Ji, W. Xia, G. Ouyang, T. Del Rose, I. Z. Hlova, B. Ueland, D. D. Johnson, C.-Z. Wang, G. Balasubramanian and P. Singh, Machine-learning and first-principles investigation of lightweight medium-entropy alloys for hydrogen-storage applications, Int. J. Hydrogen Energy, 2025, 154, 149916.
21. N. Wilson, et al., HyStor: An experimental database of hydrogen storage properties for various metal alloy classes, Int. J. Hydrogen Energy, 2024, 90, 460–469.
22. A. Rahnama, G. Zepon and S. Sridhar, Machine learning based prediction of metal hydrides for hydrogen storage, part I: Prediction of hydrogen weight percent, Int. J. Hydrogen Energy, 2019, 44, 7337–7344.
23. D. Khatamsaz, et al., Multi-Objective Materials Bayesian Optimization with Active Learning of Design Constraints: Design of Ductile Refractory Multi-Principal-Element Alloys, Acta Mater., 2022, 236, 118133.
24. D. Khatamsaz, et al., Bayesian optimization with active learning of design constraints using an entropy-based approach, npj Comput. Mater., 2023, 9(1), 49.
25. R. Arroyave, et al., A perspective on Bayesian methods applied to materials discovery and design, MRS Commun., 2022, 12(6), 1037–1049.
26. S. Nations, T. Nandi, A. Ramazani, S. N. Wang and Y. H. Duan, Metal hydride composition-derived parameters as machine learning features for material design and H2 storage, J. Energy Storage, 2023, 70, 107980.
27. J. M. Kim, T. Ha, J. Lee, Y.-S. Lee and J.-H. Shim, Prediction of pressure-composition-temperature curves of AB2-type hydrogen storage alloys by machine learning, Met. Mater. Int., 2022, 29, 861–869.
28. H. V. Thanh, et al., Hydrogen storage on porous carbon adsorbents: rediscovery by nature-derived algorithms in random forest machine learning model, Energies, 2023, 16, 2348.
29. S. Shekhar and C. Chowdhury, Prediction of hydrogen storage in metal-organic frameworks: a neural network based approach, Results Surf. Interfaces, 2024, 14, 100166.
30. P. C. Jennings, et al., Genetic algorithms for computational materials discovery accelerated by machine learning, npj Comput. Mater., 2019, 5, 46.
31. J. H. Holland, Adaptation in Natural and Artificial Systems, The University of Michigan Press, Ann Arbor, MI, 1975, p. 211.
32. D. E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Boston, MA, 1989, p. 412.
33. Data were obtained from the Hydrogen Materials Advanced Research Consortium (HyMARC) Data Hub at datahub.hymarc.org.
34. S. Suwarno, et al., Machine learning analysis of alloying element effects on hydrogen storage properties of AB2 metal hydrides, Int. J. Hydrogen Energy, 2022, 47(23), 11938–11947.
35. T. R. Somo, et al., Hydrogen storage behaviours of high entropy alloys: A review, J. Energy Storage, 2023, 73, 108969.
36. H. Huo, et al., Machine-Learning Rationalization and Prediction of Solid-State Synthesis Conditions, Chem. Mater., 2022, 34, 7323–7336.
37. X.-Y. Zhou, et al., Machine learning assisted design of FeCoNiCrMn high-entropy alloys with ultra-low hydrogen diffusion coefficients, Acta Mater., 2022, 224, 117535.
38. M. Witman, M. Allendorf and V. Stavila, Database for machine learning of hydrogen storage materials properties, 2022, https://zenodo.org/records/7324809.
39. W. Mei, G. Zhang and K. Yu, Predicting elastic properties of refractory high-entropy alloys via machine-learning approach, Comput. Mater. Sci., 2023, 226, 112249.
40. W. Huang, P. Martin and H. L. Zhuang, Machine-learning phase prediction of high-entropy alloys, Acta Mater., 2019, 169, 225–236.
41. L. Ward, A. Agrawal and A. Choudhary, et al., A general-purpose machine learning framework for predicting properties of inorganic materials, npj Comput. Mater., 2016, 2, 16028.
42. L. Ward, et al., Matminer: An open source toolkit for materials data mining, Comput. Mater. Sci., 2018, 152, 60–69.
43. J. Delhommelle, et al., Inadequacy of the Lorentz-Berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation, Mol. Phys., 2001, 99, 619–625.
44. V. P. Filippova, et al., Calculation of the parameters of the Lennard-Jones potential for pairs of identical atoms based on the properties of solid substances, Inorg. Mater. Appl. Res., 2015, 6, 1–4.
45. D. W. Jacobson, et al., Revisting Lennard Jones, Morse, and N-M potentials for metals, Comput. Mater. Sci., 2022, 205, 111206.
46. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, Sandia Report, 1993,  DOI:10.2172/10176421.
47. W. B. Paul, Molecular Dynamics Simulation, Elementary Methods. J. M. Haile, Wiley, Chichester, 1992, vol. 489, pp. 223–224.
48. W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters, J. Chem. Phys., 1982, 76, 637–649.
49. S. Nose, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81, 511–519.
50. W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A, 1985, 31, 1695.
51. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
52. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
53. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1997, 77, 3865.
54. P. Singh, M. K. Harbola, B. Sanyal and A. Mookerjee, Accurate determination of band gaps within density functional formalism, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 87(23), 235110.
55. P. Singh, M. K. Harbola, M. Hemanadhan, A. Mookerjee and D. D. Johnson, Better band gaps with asymptotically corrected local exchange potentials, Phys. Rev. B, 2016, 93(8), 085204.
56. P. Söderlind, P. E. A. Turchi PEA, A. Landa and V. Lordi, Ground-state properties of rare-earth metals: an evaluation of density-functional theory, J. Phys.: Condens. Matter, 2014, 26, 416001.
57. T. J. Giese and D. M. York, A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method, J. Chem. Theory Comput., 2018, 14, 1564–1582.
58. R. Singh, et al., Accelerating computational modeling and design of high-entropy alloys, Nat. Comput Sci, 2021, 1, 54–61.
59. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
60. F. Marques, M. Balcerzak, F. Winkelmann, G. Zepon and M. Felderhoff, Review and outlook on high-entropy alloys for hydrogen storage, Energy Environ. Sci., 2021, 14(10), 5191–5227.
61. S. Heiles and R. L. Johnston, Global optimization of clusters using electronic structure methods, Int. J. Quantum Chem., 2013, 113, 2091–2109.
62. D. Khatamsaz, B. Vela and P. Singh, et al., Bayesian optimization with active learning of design constraints using an entropy-based approach, npj Comput. Mater., 2023, 9, 49.
63. B. Cheng, Y. Li and X. Li, et al., Solid-State Hydrogen Storage Properties of Ti–V–Nb–Cr High-Entropy Alloys and the Associated Effects of Transitional Metals (M = Mn, Fe, Ni), Acta Metall. Sin. (Engl. Lett.), 2023, 36, 1113–1122.
64. M. D. Witman, et al., Machine learning-guided materials and system co-design for high-pressure hydrogen compression, ChemRxiv, 2025,  DOI:10.26434/chemrxiv-2025-352v.
65. Y. Zhang, et al., A New Phase Classifier with an Optimized Feature Set in ML-Based Phase Prediction of High-Entropy Alloys, Appl. Sci., 2023, 13(20), 11327.
66. A. Rahnama, et al., Machine learning based prediction of metal hydrides for hydrogen storage, part I: Prediction of hydrogen weight percent, Int. J. Hydrogen Energy, 2019, 44, 7337–7344.
67. K. S. Chandana and R. Kamesh, Machine learning insights into prediction of H2 gravimetric capacity in Mg-based pure metal alloys, Int. J. Hydrogen Energy, 2024, 77, 695–711.
68. Z. Ding, Y. Li and H. Jiang, et al., The integral role of high-entropy alloys in advancing solid-state hydrogen storage, Interdiscip. Mater., 2025, 4, 75–108.
69. H. Shao, et al., Preparation and hydrogen storage properties of nanostructured Mg–Ni BCC alloys, J. Alloys Compd., 2009, 477, 301–306.
70. A. Resnik, M. Stioui, A. Grayevsky and D. Shaltiel, Reaction kinetic Thermal desorption spectra of hydrogen from LaNi₅, J. Less-Common Met., 1987, 131, 117–124.
71. A. Agafonov, et al., Destabilizing high-capacity high entropy hydrides via earth abundant substitutions: From predictions to experimental validation, Acta Mater., 2024, 276, 120086.
72. M. M. Nygard, et al., Local order in high-entropy alloys and associated deuterides – a total scattering and Reverse Monte Carlo study, Acta Mater., 2020, 199, 504–513.
73. N. Pineda-Romero, et al., The effect of 10 at.% Al addition on the hydrogen storage properties of the Ti0.33V0.33Nb0.33 multi-principal element alloy, intermetallic, 2022, 146, 107590.
74. P. Mohan, et al., Influence of β-phase stability in elemental blended Ti-Mo and Ti-Mo-Zr alloys, Micron, 2021, 142, 102992.
75. M. Okada, T. Kuriiwa, T. Tamura, H. Takamura and A. Kamegawa, Ti–V–Cr b.c.c. alloys with high protium content, J. Alloys Compd., 2002, 330–332, 511–516.
76. H. Z. Hu, C. M. Ma and Q. J. Chen, Mechanism and microstructural evolution of TiCrVFe hydrogen storage alloys upon de-/hydrogenation, J. Alloys Compd., 2021, 877, 160315.
77. X. Y. Chen, B. Liu, S. B. Zhang, X. Ding and R. R. Chen, Effect of heat treatment on microstructure and thermal stability of Ti₁₉Hf₄V₄₀Mn₃₅Cr₂ hydrogen storage alloy, J. Alloys Compd., 2022, 917, 165355.
78. H. C. Lin, K. M. Lin, K. C. Wu, H. H. Hsiung and H. K. Tsai, Cyclic hydrogen absorption–desorption characteristics of TiCrV and alloys, Int. J. Hydrogen Energy, 2007, 32, 4966–4972.
79. H. Itoh, H. Arashima, K. Kubo, T. Kabutomori and K. Ohnishi, Improvement of cyclic durability of BCC structured Ti–Cr–V alloys, J. Alloys Compd., 2005, 404–406, 417–420.
80. P. Ruz, S. Banerjee, R. Halder, A. Kumar and V. Sudarsan, Thermodynamics, kinetics and microstructural evolution of Ti_0.43Zr_0.07Cr_0.25V_0.25 alloy upon hydrogenation, Int. J. Hydrogen Energy, 2017, 42, 11482–11492.
81. P. Singh, A. V. Smirnov and D. D. Johnson, Ta-Nb-Mo-W refractory high-entropy alloys: Anomalous ordering behavior and its intriguing electronic origin, Phys. Rev. Mater., 2018, 2, 055004.
82. G. M. Lu, et al., Explainable machine learning for hydrogen diffusion in metals and random binary alloys, Phys. Rev. Mater., 2023, 7(10), 105402.

---