Predicting hydrogen storage in MOFs via machine learning.

The H2 capacities of a diverse set of 918,734 metal-organic frameworks (MOFs) sourced from 19 databases is predicted via machine learning (ML). Using only 7 structural features as input, ML identifies 8,282 MOFs with the potential to exceed the capacities of state-of-the-art materials. The identified MOFs are predominantly hypothetical compounds having low densities (<0.31 g cm-3) in combination with high surface areas (>5,300 m2 g-1), void fractions (∼0.90), and pore volumes (>3.3 cm3 g-1). The relative importance of the input features are characterized, and dependencies on the ML algorithm and training set size are quantified. The most important features for predicting H2 uptake are pore volume (for gravimetric capacity) and void fraction (for volumetric capacity). The ML models are available on the web, allowing for rapid and accurate predictions of the hydrogen capacities of MOFs from limited structural data; the simplest models require only a single crystallographic feature.

Introduction

Hydrogen (H2) is considered to be a future automotive fuel.1, 2, 3, 4, 5, 6 This potential reflects its high specific energy compared with competing fuels, such as natural gas and gasoline, and the ability of H2 to be produced renewably and consumed without CO2 emissions., Nevertheless, the adoption of hydrogen in mobile applications, such as fuel cell (FC) vehicles has been limited by its low volumetric energy density.,, Consequently, the design of low-cost H2 storage systems that overcome these volumetric limitations has been the focus of recent research.,8, 9, 10, 11, 12 At present, FC vehicles employ storage systems based on gaseous H2 compressed to pressures up to 700 bar. This approach is costly and can incur limitations in driving range.,,,

Storage based on adsorption in porous hosts is an alternative to high-pressure compression. Due to their high gravimetric densities, fast kinetics, and reversibility, metal-organic frameworks (MOFs) have emerged as one of the most promising classes of hydrogen sorbents., MOFs are crystalline materials formed by the self-assembly of inorganic metal clusters and organic linkers.16, 17, 18, 19, 20, 21, 22 By virtue of their building-block structure and the large number of potential components, the number of MOFs is potentially limitless.21, 22, 23, 24, 25 Further modifications to MOF chemistry can be achieved by introducing functional groups, substituting different metals, and by mixing metals and/or linkers.26, 27, 28

Despite these many possibilities, a relatively small fraction of MOFs have been synthesized., While the crystal structures of these “real” MOFs are available in the Cambridge Structural Database (CSD),, many exhibit disorder, missing atoms, or have negligible porosity; consequently, these materials are not immediately amenable to assessment via computational modeling.,31, 32, 33, 34, 35

One way to bypass these complications is through computational design. To date, nearly a million “hypothetical” MOFs have been reported,,36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46 and it is reasonable to expect that many more materials will be proposed.47, 48, 49, 50, 51 High-throughput screening using Grand Canonical Monte Carlo (GCMC)52, 53, 54, 55, 56 has been successful in identifying promising candidates with superior gas storage capacities on sub-sets of these catalogs.,,,,,57, 58, 59, 60 Nevertheless, given the large number of possibilities, a systematic search across all of these materials is challenging even with high-throughput techniques., Furthermore, differences in the implementation (i.e., use of different temperature/pressure conditions or interatomic potentials) can complicate comparisons between screening studies. Thus, more efficient and consistent screening approaches are desirable for predicting the gas storage properties of MOFs in existing and future databases.

Machine learning (ML) could provide a path forward.62, 63, 64, 65 For ML to be helpful, access to high-quality training data is essential. Unfortunately, training on experimental H2 storage data in MOFs is non-trivial,,,66, 67, 68: experimental uptake data are generally restricted to a relatively small number of MOFs, and can depend sensitively upon the experimental conditions and the purity of the sample.,, Employing a dataset based on a consistent set of computational predictions may be a better choice.,

Earlier work has demonstrated that accurate isotherms for H2 uptake in MOFs can be predicted using the pseudo-Feynman-Hibbs potential (to describe H2) combined with general interatomic potentials to describe the MOF.,,, This approach was used to screen a database of 5,309 real MOFs, from which IRMOF-20 was identified and experimentally demonstrated to have a favorable balance of high gravimetric and volumetric H2 density. In a follow-on study, a larger database of 495,305 MOFs was compiled from several publicly available databases (see Table S1 for details).,,,,36, 37, 38, 39, 40, Following a pre-screen based on crystallographic properties and empirical correlations, the H2 capacities of a subset of 43,777 MOFs were evaluated using GCMC. Three additional MOFs—SNU-70, UMCM-9, and PCN-610/NU-100—were identified and shown experimentally to out-perform the leading MOF candidate, IRMOF-20.

The database of MOF properties generated in these previous studies presents an opportunity to develop ML models that can predict H2 uptake across even larger MOF datasets., Table 1 summarizes previous ML studies of H2 storage in MOFs. (Reports employing ML for other adsorbates, such as CH4,, CO2,, and N2, are summarized in Table S2.) To the best of our knowledge, ML was first used to predict H2 uptake in compounds from the Nanoporous Materials Genome. A neural network (NN) was used to predict usable capacities on a test set of ∼1,000 compounds, including MOFs. In the same year, Borboudakis et al. predicated H2 capacities in 100 MOFs using 92 binary features related to a MOF's linker, metal cluster, and functional group(s). Ridge linear regression (RR)76, 77, 78 and support vector machine (SVM), algorithms were used to predict gravimetric capacity. Later, Bucior et al predicted the H2 capacities of 54,776 MOFs extracted from the CSD using multilinear regression (MLR). The models were trained using the energetics of H2-MOF interactions and the usable volumetric capacities predicted by GCMC. More recently, ML was used to predict H2 storage capacities in 105 hypothetical MOFs constructed from 17 different topologies, 4 distinct metal clusters, and 5 unique organic linkers. NN models employing 11 features were trained to predict total volumetric uptake at various temperatures and pressures.

Table 1

Summary of recent studies that use machine learning to predict H2 adsorption in MOFs

Study	ML features	ML method	Properties predicted	Accuracy
Anderson et al.43	epsilon, temperature, pressure, ρcrys, vf, vsa, mpd, lcd, alchemical catecholate site density, unit cell volume	neural network76	total volumetric H2 for pressures 0.1, 1, 5, 35, 65, and 100 bar at 77, 160, and 295 K	AUE = 0.75–2.93 g-H2 L−1
Bucior et al.80	energetics of MOF-guest interactions	multilinear regression with LASSO76	deliverable H2 storage capacity between 2 and 100 bar at 77 K	R2 = 0.96; AUE = 1.4–3.4 g-H2 L−1; RMSE = 3.1–4.4 g-H2 L−1
Borboudakis et al.63	92 binary features based on linker, metal cluster, and 12 functional groups	ridge linear regression and support vector machine with polynomial/Gaussian kernel76, 77, 78	total H2 storage capacity at 1 bar and 77 K	AUE = 0.47 (ridge regression), 0.50 (SVM) g-H2 g−1-MOF
Thornton et al.61	adsorption energy, ρcrys, vf, gsa, vsa, lcd	neural network76	net H2 capacity for pressure swing between 1 and 100 bar at 77 and 298 K	R2 = 0.88; RMSE = 3.6 g-H2 L−1

ρcrys, vf, vsa, mpd, lcd represent single-crystal density, void fraction, volumetric surface area, maximum pore diameter, and largest cavity diameter, respectively. R2, AUE, and RMSE represent the coefficient of determination, average unsigned error, and root-mean-square error, respectively.

Expanding upon these previous reports, this study applies ML to explore a large database of 918,734 known and proposed MOFs. The database was assembled from a diverse collection of publicly available MOF repositories,,,,,,36, 37, 38, 39, 40, 41, 42, 43, 44, 45,, and allows for a wide-ranging and consistent assessment of H2 uptake in MOFs.

Here, the extremely randomized trees (ERT), algorithm was identified as the most accurate ML model for predicting H2 uptake. A training set comprising 24,674 MOFs was sufficient to enable accurate predictions of usable capacities across 820,039 unseen compounds. These predictions were made using a small set of seven crystallographic features as input: single-crystal density, pore volume, gravimetric and volumetric surface area, void fraction, largest cavity diameter, and pore limiting diameter. Importantly, ML identified 8,282 MOFs—8,187 appropriate for pressure swing (PS) operation and 95 for temperature-PS (TPS) use—with the potential to exceed both the gravimetric and volumetric capacities of state-of-the-art materials. These compounds are comprised predominantly of hypothetical MOFs, and exhibit low densities (<0.31 g cm−3) in combination with high surface areas (>5,300 m2 g−1), void fractions (∼0.90), and pore volumes (>3.3 cm3 g−1). In addition to identifying high-capacity MOFs, the relative importance of the input features is quantified; dependencies on the ML algorithm and training set size and are also assessed. The most important features for predicting H2 uptake are pore volume (for gravimetric capacity) and void fraction (for volumetric capacity). A simplified model using only two input features is demonstrated to predict capacities with high accuracy—within 0.2 wt % and 1.4 g-H2 L−1 of more expensive Monte Carlo calculations. The ML models are available for use via the web, allowing for rapid and accurate predictions of hydrogen capacities with only a small amount of structural data required as input.

Methods

MOF database

A database of crystal structures for 918,734 MOFs was created by combining 19 existing databases.,,,,,36, 37, 38, 39, 40, 41, 42, 43, 44, 45,, Table 2 summarizes the source databases and the number of MOFs contained in each. Out of these 19 databases, only the UM, CSD,, and CoRE, databases contain data on MOFs that have been previously synthesized. (MOFs listed in these datasets are referred to as “real” MOFs.) The remaining databases contain data for proposed, or “hypothetical”, MOFs. The seven crystallographic properties for all MOFs in the database were calculated using the zeo++ code, with a probe radius of 1.86 Å. These data are available at the HyMARC data hub. Additional details can be found in our previous work. These properties include: single-crystal density (d), pore volume (pv), gravimetric surface area (gsa), volumetric surface area (vsa), void fraction (vf), largest cavity diameter (lcd), and pore limiting diameter (pld).

Table 2

MOF datasets employed in this study

Source	Databaseidentity	No. of MOFs
Goldsmith et al.,31 Chung et al.,33 Moghadam et al.,29 Groom et al.30	real MOFs:UM31+CoRE33+CSD29,30	15,235
Chung et al.34	CoRE 201934	14,142
Moghadam et al.,29 Groom et al.30	aCSD 2017 additional29,30	48,696
Martin et al.38	mail-order38	112
Bao et al.46	in silico deliverable46	2,816
Bao et al.39	in silico surface39	8,885
Witman et al.40	MOF-74 analogs40	61
Colón et al.59	ToBaCCo59	13,512
Gomez-Gualdron et al.45	Zr-MOFs45	204
Wilmer et al.36	Northwestern36	137,000
Aghaji et al.,37 Boyd et al.85,86	bUniv. of Ottawa37,85,86	317,462
Lan et al.81	BJT MOFs81	303,793
Chung et al.41,87	cR-WLLFHS41,87	51,163
Li et al.82	MTV82	11,555
Anderson et al.42	CSM-2018-I42	117
Anderson et al.43	CSM-2018-II43	32
Anderson et al.44	CSM-2019-I44	99
Ahmed et al.1	in-house1	18
	total	918,734

A subset of the CSD 2017 MOF dataset, whose crystallographic properties were found to exhibit extremely low values (e.g., GSA ~0) in a previous study.

A recent version of this database is available publicly;, however, this study employs an earlier version that was shared privately.

A curated subset of the Northwestern database.

A previous study examined a subset of the present database, wherein the hydrogen uptake in 495,305 MOFs was estimated using the Chahine rule.,, Subsequently, usable uptake in a portion of this subset comprising 43,777 MOFs predicted to be promising based on the Chahine rule was evaluated using GCMC. This GCMC-evaluated dataset contained a mix of real and hypothetical MOFs: 15,235 real MOFs were sourced from the UM, CoRE, and CSD,, and 28,542 hypothetical MOFs were extracted from the mail-order, in silico deliverable, in silico surface, MOF-74 analogs, ToBaCCo, Zr-MOFs, Northwestern, University of Ottawa,,, and in-house hypothetical MOF databases (see Ahmed et al. or Table S1 for details).,,,,36, 37, 38, 39, 40 Hydrogen uptake isotherms for two operating conditions were predicted: for an isothermal PS at T = 77 K between 5 and 100 bar, and for a combined TPS between 77 K/100 bar (filled state) and 160 K/5 bar (empty state). UG and UV capacities were then calculated based on the isotherm data.

In addition to the 43,777 MOFs examined in Ahmed et al., in this study GCMC isotherms were evaluated for an additional 54,918 MOFs (see Ahmed et al. and Table S1 for further details). These additional MOFs were selected at random from the 495,305-entry HyMARC database and therefore represent a more diverse sampling of the MOF property space. To this dataset, 423,429 additional compounds were added from 7 additional datasets: BJT (Beijing, Jiangsu, Tianjin) MOFs, R-WLLFHS,, MTV, CSM-2018-I, CSM-2018-II, and CSM-2019-I, and selected MOFs from the CSD 2017 dataset., Subsequently, the capacities of the MOFs from these additional datasets were predicted by the ML models without retraining (i.e., no MOFs from these datasets were used for training or testing, and none of their isotherms were evaluated in advance with GCMC). In total, the dataset employed in this study contains H2 uptake data for 98,695 MOFs and crystallographic property data for 918,734 MOFs.

The present dataset includes approximately 74,000 MOFs having open metal sites (OMS), comprising roughly 8% of the total dataset. As the interatomic potential used in our GCMC calculations is not tuned to capture the unique aspects of the H2-OMS interaction, it is possible that the calculated capacities for this class of MOFs will be less accurate. Figure S1 and Table S3 compare experiments and the present GCMC calculations of H2 capacities across a benchmark set of OMS MOFs discussed by García-Holley et al. and in our previous work. These data show that GCMC calculations using the pseudo-Feynman-Hibbs potential are in good agreement with experimental data for these OMS MOFs. The good agreement between theory and experiments is a consequence of the low temperature operating conditions used in our study, combined with the relatively low density of OMS in these MOFs.

ML models

The No Free Lunch Theorem implies that the optimal choice of ML algorithm is problem specific. The differing performance of the algorithms summarized in Tables 1 and S2 is consistent with this notion. Identifying the best algorithm for a given dataset requires comparing multiple ML methods, each with optimized hyperparameters. Unfortunately, few comparisons of ML methods for gas adsorption exist; although dozens of ML algorithms are available,76, 77, 78, 79,,90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104 only RR,76, 77, 78 MLR, SVM,, and NN have been examined for predicting H2 storage.,,,, This study casts a wider net by comparatively assessing 14 ML algorithms (Table 3).76, 77, 78, 79,,90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104

Table 3

Machine learning regression algorithms employed in this work

Machine learning algorithm	Abbreviation
Extremely randomized trees76,83,103,104	ERT
Boosted decision trees76,92,102, 103, 104	BDT
Bagging with decision trees76,9093,103,104	B/DT
Random forest76,90,94,103,104	RF
Bagging with random forest76,93,94,103,104	B/RF
Gradient boosting76,92,95,102, 103, 104	GB
Decision trees76,90,103104	DT
Nu-support vector machine with radial basis function (RBF) kernel76,79,90,96,98,103104	Nu-SVM/RBF-K
Support vector machine RBF kernel76,79,90,97,98,103104	SVM/RBF-K
Support vector machine with linear kernel76,79,96,99,103104	SVM/L-K
Linear regression76, 77, 78,99,100,103104	LR
Ridge regression76, 77, 78,99,100,103104	RR
K-nearest neighbors76,90,101,103104	K-NN
AdaBoost76,92,102, 103, 104	AB

The crystallographic properties of MOFs are known to correlate with H2 capacities.,,,105, 106, 107, 108 The ML models developed here exploit these correlations by adopting only crystallographic properties as input features. Moreover, the number of features was restricted to a small set comprising seven properties: d, pv, gsa, vsa, vf, lcd, and pld. These are the same properties employed in our previous work.,,, Figure S2 shows the distribution of crystallographic properties for the training, test, and unseen datasets. Also, Table S4 summarizes five descriptive (minimum, maximum, mean, median, and percent of 0's) and two distribution statistics (skew and kurtosis) of all crystallographic features for the training, test, and unseen datasets. (The details regarding these statistics and the definitions of skew and kurtosis can be found in Table S4.) The maxima and minima of the features in the training set establish the validity ranges of the ML models developed here.

The goal of the ML models is to predict four output properties: UG and UV for each of PS and TPS operating conditions. This was accomplished by developing separate ML models for each of the four targeted capacities. Figure S3 illustrates the overall work flow.

The existing dataset of 98,695 MOFs (for which both crystallographic and capacity data are available) was initially split into training and test sets of 74,201 and 24,674 MOFs, respectively, after shuffling the entire dataset. ML algorithms,,, (Table 3) were implemented using the Scikit-learn library. Both scaled and unscaled features were used in training ML models. Ten-fold cross-validation was used to optimize the hyperparameters of each model. The performance of the ML algorithms was assessed by comparing the predicted H2 capacities with the capacity predicted by GCMC for the MOFs in the test set. The metrics used for the performance assessment of ML models were the R2, AUE, RMSE, MAE, and . Additional details regarding these calculations can be found in supplemental note S2 of the supplemental information.

Dataset size

An obstacle to wider adoption of ML in materials science is the availability of sufficient quantities of high-quality training data., Unfortunately, it is not yet clear how much data are needed to construct a useful ML model for a given system. Fernandez et al. found that a reasonable balance between accuracy (R2 ∼ 0.85–0.93) and computational expense for predicting methane storage in MOFs was achieved for a training set containing data on 10,000 MOFs with 3 features. In contrast, Fanourgakis et al. showed that a much smaller training set of ∼1,000 MOFs was sufficient to predict methane uptake when using six crystallographic features and four fictitious features. The different training set sizes required in these previous studies arise from the differing numbers and types of features used.

This study explores this issue further by systematically examining the effect of training set size, and the training set to test set ratio, on ML accuracy. For each of the four targeted capacity outputs, 100 independent ML models were developed by varying the size of the training set between 100 and 74,000 MOFs (see Table S5 for a list of the training set sizes). The four best-performing ERT ML algorithms identified earlier were used with 10-fold cross-validation. The resulting models were assessed using a common test set of 24,674 MOFs.

Feature importance/selection

The well-known Chahine rule proposes a linear correlation between gravimetric surface area and excess gravimetric H2 capacity in adsorbents., Nevertheless, the Chahine rule overpredicts H2 capacities for MOFs with high surface areas, and has not been extended to predict usable capacities.,, Hence, a model for predicting H2 uptake that is more general than the Chahine rule, yet requires limited input data, would be very helpful. In principle, ML could be used to generate such a predictive model if the features that are the most important for predicting H2 uptake could be identified. Along these lines, Pardakhti et al. reported improved accuracy in predicting CH4 adsorption when using a combination of (7) crystallographic and (19) chemical features. Recently, Moosavi et al. explored feature importance in predicting the synthesis of MOFs.

This study determines the minimum number and optimal combination of crystallographic features necessary to achieve a specified accuracy in predicting H2 uptake. The relative importance of the input features was assessed for all possible univariate and multivariate feature combinations using ERT ML models. The number of multivariate feature combinations, M, is given by:, where n = 7 is the total number of available features, and 1 ≤ n ≤ 7 is the number of features used as input to a given ML model. A total of 127 feature combinations are possible. ML models were developed for each of these feature combinations for each of the 4 output capacities, resulting in a total of 508 distinct ML models. All models were trained using a dataset of 74,021 MOFs and tested on a common set of 24,674 MOFs. Ten-fold cross-validation was used for tuning and validating the models using only the training set. Univariate feature importance was further assessed using (1) Pearson's correlation coefficient (r),116, 117, 118 (2), Breiman and Friedman's tree-based algorithm as implemented in Scikit-learn,, and (3) the permutation importance method as implemented in rfpimp package. Additional details regarding these methods can be found in Figure S7.

Results

Evaluating ML algorithms

Tables S6–S9 illustrate the effect of several feature scaling methods on the performance of the ML algorithms examined here. Only the SVM family of models (SVM/L-K, SVM/RBF-K, and Nu-SVM/RBF-K),,,,, were impacted by the choice of scaling method.

Figure 1 compares the accuracy of the ML algorithms for predicting hydrogen uptake in MOFs. Coefficient of determination (R2) and average unsigned error (AUE) were used as performance metrics. SVM variants were trained using min-max feature scaling; unscaled features were used in training the remaining models. The performance of the algorithms as measured by four additional metrics—root-mean-square error (RMSE), explained variance (EV), median absolute error (MAE), and Kendall rank correlation coefficient ()—is reported in Tables S6–S9.

Figure 1

Comparison of ML algorithms for predicting hydrogen uptake in MOFs

(A and C) Left and (B and D) right panels report performance for PS and TPS conditions, respectively. (A and B) Top and (C and D) bottom panels report performance for usable gravimetric and volumetric capacities, respectively. The abbreviations for the ML methods are defined in Table 3.

Overall, these data indicate that the tree-based ensemble methods are superior to the other methods examined. In particular, the ERT,, algorithm exhibited the best performance overall. Boosted decision trees,,90, 91, 92,, random forest,,, and Bagging algorithm variants,,,, (with tree-based base estimators) are nearly as accurate. The R2 values for ERT predictions exceed 0.997 for gravimetric capacities, which are equivalent to errors of ∼0.14 wt %. Volumetrically, the accuracy of the ERT algorithm is slightly worse than its gravimetric performance: R2 = 0.967–0.984, equivalent to errors of ∼1.1 g-H2 L−1 on average. In general, the worst-performing algorithms were linear regression, ridge regression, and SVM with linear kernel. For these algorithms R2 varies between 0.913 and 0.992 depending on the conditions (i.e., gravimetric/volumetric and PS/TPS). As expected, the linear nature of these algorithms fails to fully capture the nonlinear dependence of output capacities on the multiple input features.

Figure 1 also shows that all the algorithms tested yield more accurate predictions of usable gravimetric (UG) capacities compared with those for usable volumetric (UV) capacities. Likewise, all algorithms more accurately predict usable capacities under PS conditions than under TPS conditions. This reflects the fact that the functional relationships between output capacities (UG/UV) and input features under PS and TPS conditions are likely different, as was observed in previously reported structure(feature)-property(capacity) relationships.,, Table 4 summarizes the performance of the ERT algorithm in further detail. A comparison of Tables 1 and 4 indicates that the accuracy of the present ML models surpass previously reported models for H2 uptake. Furthermore, the present models also appear to be an improvement over earlier models that aim to predict the adsorption capacities of MOFs for any gas species, Table S2. This improved performance can be attributed to the exploration and optimization of multiple ML algorithms, use of an appropriate feature set, and the relatively large size of the present training set.

Table 4

Performance of the extremely randomized trees ML algorithm in predicting UG and UV H2 capacities of MOFs under PS and TPS conditions

H2 capacity type	R2	AUE (capacity units)	RMSE (capacity units)	Kendall τ	MAE (capacity units)
UG at PS (wt %)	0.997	0.14	0.18	0.961	0.10
UV at PS (g-H2 L−1)	0.984	0.97	1.40	0.922	0.69
UG at TPS (wt %)	0.997	0.16	0.23	0.966	0.10
UV at TPS (g-H2 L−1)	0.967	1.32	1.92	0.819	0.91

R2, AUE, RSME, and MAE represent the coefficient of determination, average unsigned error, root-mean-squared error, and median absolute error, respectively.

Figure 2 illustrates the degree of agreement between ERT ML predictions and GCMC calculations of usable H2 capacities under PS conditions as a function of MOF source database (Figure S4 shows similar data for TPS conditions; see also Table 4). As mentioned above, the present ML models more accurately predict UG capacities than UV capacities. The largest differences between ML and GCMC capacities (Figures 2C, 2F, S4C, and S4F) primarily occur for the real MOF dataset. In principle, these differences may arise either from ML overfitting or from inaccurate GCMC predictions caused by non-ideal/incomplete MOF crystal structure data (i.e., missing atoms, disorder, etc.), as mentioned in previous studies.,,,123, 124, 125 ERT algorithms are fairly robust against overfitting. To examine the possibility for overfitting, test set errors were compared with training set errors, as shown in Figure S5 and Table 4. These data suggest that the outliers are not a consequence of over fitting; hence, inaccuracies in the crystal structure data are proposed as the most likely source of this disagreement.,,,123, 124, 125

Figure 2

Performance of the ERT algorithm with respect to GCMC calculations for predicting usable H2 capacities in MOFs

Data were collected at 77 K for a pressure swing (PS) between 100 and 5 bar on a test set of 24,674 MOFs. Different colors represent different categories of MOFs. (A–C) Top and (D–F) bottom panels illustrate performance for usable gravimetric and volumetric capacities, respectively. (A and D) Agreement between ML and GCMC predictions. (B and E) Difference between ML and GCMC as a function of GCMC capacity. (C and F) Distribution of differences in predictions between ML and GCMC.

Effect of training set size

Figure 3 illustrates the impact of training set size on the accuracy of the ERT ML models, as quantified using R2 and AUE (Table S5 summarizes the dataset sizes used in these plots). For training sets containing more than 5,000 MOFs, R2 and AUE vary slowly and in a monotonic fashion, with AUE decreasing and R2 increasing. The accuracy of the models is more sensitive to the size of the training set for smaller training sets containing roughly 5,000 or fewer MOFs. Figure S6 highlights the variation in performance for these smaller training sets.

Figure 3

ML performance versus training set size

Performance of ERT ML models for predicting usable (A) gravimetric and (B) volumetric H2 capacity as a function of training set size and the ratio of training to test set size. One hundred different training sets, ranging in size between 100 and 74,021 MOFs were examined. A common set of 24,674 MOFs was used for testing. Performance is quantified using R2 (left axis, black) and the AUE (right axis, blue and red for UG and UV, respectively). Lines represent a power law fit to the data.

The trends AUE as a function of training set size can be fit to a power law expression of the form AUE(m) = αmβ + γ, wherem represents the size of the training set and β is the power law exponent. Fitting this model to the data shown in Figure 3 reveals that the AUE for UG converges faster with training set size (β = −0.37 and −0.43) than it does for UV (β = −0.16 and −0.23). A full tabulation of the power law parameters is given in Table S10. Based on these power law expressions, one can determine the necessary size of the training set to achieve a desired level of accuracy. For example, assuming PS operation, to achieve an AUE of approximately 0.25 wt % and 1.5 g-H2 L−1 requires training set sizes (for UG and UV) of less than 300 MOFs randomly selected from the diverse datasets used here.

Univariate feature importance

Figure 4 illustrates the relative importance of the seven crystallographic features in predicting usable hydrogen uptake in MOFs. Feature importance was determined by developing ERT models for each single feature individually. Additional details for these models are provided in the supplemental information. Based on these models, it is evident that pore volume (pv) and void fraction (vf) are the dominant features in predicting H2 capacity; these two properties appear as the first- or second-most important single features regardless of operating condition or capacity type. The importance of these features can be rationalized by two factors. First, based on the empirical Chahine rule, the pore volume of an MOF correlates with its excess uptake. Second, pore volume and void fraction are related (since pv = vf d−1)—MOFs with larger pv have larger vf, and vice versa.

Figure 4

Univariate feature importance in predicting usable H2 capacities in MOFs

Feature importance was determined by developing distinct ERT models for each individual feature. The accuracy of the resulting models was assessed using R2 (left axis; black dataset) and AUE (right axis; red dataset). Models were trained on a dataset of 74,201 MOFs and tested on a set of 24,674 MOFs. pv, pore volume; d, density; vf, void fraction; gsa, gravimetric surface area; pld, pore limiting diameter; lcd, largest cavity diameter; vsa, volumetric surface area.

Conversely, the largest cavity diameter (lcd) and volumetric surface area (vsa) are the single features whose ML models yield the lowest accuracy. The relative importance of the individual features for predicting UG capacities is: pv > d > vf > gsa > pld > lcd > vsa. This ordering is the same for PS and TPS conditions. In contrast, the importance ordering for UV capacities differs based on the operating condition. Nevertheless, vf and pv remain the two most important single features for both UV conditions, in that order (Figure 4).

Despite their limited input, the single-feature ML models illustrated in Figure 4 achieve high accuracy. For example, any of the three independent models for UG-PS based only on pv, d, or vf can predict capacities with R2 > 0.95 and with AUE of less than 0.5 wt %. The accuracy and simplicity of the univariate ML models suggest that they can be used to quickly screen new MOFs for their utility in hydrogen storage. To that end, optimized single-feature ML models for the four categories of usable capacities considered here have been made available for use on the web with an interactive web form or with a python API. Furthermore, the ML models can be downloaded via figshare. These models take as input either pv (for UG predictions) or vf (for UV predictions) of a given MOF. These input data can be quickly calculated from a MOF's crystal structure using modern structure analysis codes.,,127, 128, 129, 130 As shown in Figure 4, these models can predict UG with an average error of less than 0.4 wt %, and UV with errors less than 2.2 g-H2 L−1.

Figure S7 compares the single-feature importance assessments based on ERT ML models (as reported in Figure 4) with three popular methods for determining feature importance: Pearson's correlation coefficient (r),116, 117, 118 Breiman and Friedman's tree-based algorithm as implemented in Scikit-learn,, and the permutation importance method as implemented in the rfpimp package. It is clear that the feature importance methods do not reproduce in detail the rank ordering of feature importance that is suggested by our ERT ML models. Nevertheless, good agreement is evident more broadly. For example, in the case of UG (Figures S7A and S7C), the three feature importance methods suggest that in aggregate pv is the most important feature, while vsa is the least, in agreement with the ERT models (Figures 4A and 4B). Similarly, for UV, the importance methods suggest that vf and lcd are among the most and least important features, respectively. This is the same trend found in the univariate ERT models (Figures 4C and 4D).

Multivariate feature importance

Figure 5 illustrates how the accuracy of the ML models varies with the number and combination of features. Assuming 7 features, 27 – 1 = 127 possible combinations exist. For a given number of features, Figure 5 plots the combination of features resulting in the highest accuracy model. (The supplementary file [Table S11] summarizes the performance for all 508 possible feature combinations and capacity/operating condition types.) As expected, Figure 5 shows that ML accuracy generally increases as the number of input features increases. As previously discussed, when limited to a single feature, vf yields the best accuracy for predicting UV, while pv is the best choice for UG. When the feature set is extended to 2 features, the combination of d and pv is the optimal choice among the possible pairs regardless of the capacity (UG versus UV) or operating condition (PS versus TPS). For larger numbers of features, the optimal feature combination depends upon the operating condition and the capacity type. Based on the AUE, whose value tends to plateau as more features are added, highly accurate ML models can be generated using only 5 input features (Table 5). These data lend further support to the notion that the accuracy of a given ML model depends on both the number and identity of the input features. As a slightly more accurate alternative to the univariate web models described above, a subset of the present multivariate ML models that use 4, 5, and 7 input features are also available on the web using an interactive web form and via a python API. The ML models can also be downloaded via figshare.

Figure 5

Multivariate feature importance in predicting usable H2 capacities in MOFs

The accuracy of ERT ML models, as determined by R2 and AUE, was determined as a function of the number and combination of input features. Each data point represents the most accurate feature combination for a given number of features. ERT models were trained on a dataset of 74,201 MOFs. R2 and AUE were calculated using a test of 24,674 MOFs. Feature abbreviations are defined in Figure 4.

Table 5

Optimal combinations of features for predicting UG and UV H2 storage capacities at PS and TPS conditions

Condition	Feature combination	No. of features	R2	AUE	RMSE	Kendall τ
UG at PS	gsa, vf, pv, lcd, pld	5	0.997	0.14 wt %	0.19 wt %	0.959
UG at TPS	d, vsa, pv, lcd, pld	5	0.996	0.18 wt %	0.25 wt %	0.959
UV at PS	vsa, vf, pv, lcd, pld	5	0.983	1.01 g-H2 L−1	1.45 g-H2 L−1	0.920
UV at TPS	vsa, vf, pv, lcd, pld	5	0.961	1.41 g-H2 L−1	2.10 g-H2 L−1	0.814

H2 uptake in unseen MOFs

Figure 6 illustrates the H2 storage capacities of 820,039 MOFs as predicted by the 7-feature ERT ML models developed here. (This dataset is publicly accessible via HyMARC data hub.) These MOFs are referred to as “unseen”, in that they have not been included in the training or test sets used to develop the models. Figures 6A and 6B show UV capacities as functions of UG capacities under PS and TPS conditions, respectively. Both plots exhibit a rapid increase in UV at low values of UG, and reach a maximum in UV at UG values of approximately 9 wt %. Beyond the maximum, UV decreases relatively slowly with increasing UG. These trends are consistent with our earlier findings derived from GCMC calculations on smaller datasets.,,

Figure 6

ML predictions of H2 capacities for 820,093 unseen MOFs

Predicted capacities for (A) PS and (B) temperature + PS operation. Colors indicate the originating database for a given MOF. (C and D) Validation of ML-predicted capacities for the highest-capacity MOFs identified by ML; shown in the rectangular regions in (C and D) using GCMC simulations. For comparison, the capacities of PCN-610/NU-100 (PS: 10.1 wt %, 35.5 g-H2 L−1) and MOF-5 (TPS: 7.8 wt %, 51.9 g -H2 L−1) are shown.

In the case of PS operation, the maximum UV across the MOFs in the dataset is 37.4 g-H2 L−1; for TPS operation the maximum UV is 48.5 g-H2 L−1. In the case of UG, the maximum value predicted is 39 wt % for PS operation and 42 wt % for TPS. These values can be placed in context by comparing against the Department of Environment hydrogen storage targets, which stipulate system-level hydrogen densities of 5.5 wt % and 40 g-H2 L−1 by 2025 and 6.5 wt %/50 g-H2 L−1 longer-term (“Ultimate target”). Given that the tank and balance-of-plant for the storage system have non-zero mass and volume, the MOFs examined here cannot meet the Ultimate target for UV, regardless of operating condition. More optimism exists, however, for meeting the gravimetric targets given the high UG exhibited by these systems on a MOF-only basis. Of course, an additional challenge is to identify MOFs that excel both gravimetrically and volumetrically.,,,,

It is also helpful to compare the performance predictions in Figures 6A and 6B with that of state-of-the-art materials. In the case of PS operation, our previous study demonstrated that PCN-610 (NU-100) exhibits a hydrogen capacity of 10.1 wt % and 35.5 g-H2 L−1, which, to our knowledge, is the best combination of gravimetric and volumetric capacities reported for any MOF under these conditions. The data in Figure 6A reveal that 16,345 MOFs can, in principle, exceed this capacity on both a UG and UV basis. In the case of TPS operation (Figure 6B), MOF-5 remains the benchmark, which a measured capacity of 7.8 wt % and 51.9 g-H2 L−1. Figure 6D shows that only 21 MOFs out-perform MOF-5 under these conditions.

Regarding the accuracy of the present ML predictions, Table 4 shows that the AUE of these models are on the order of 0.15 wt % and 1.3 g-H2 L−1. Although these errors are small, a more rigorous validation of the ML can be achieved with GCMC calculations. Thus, GCMC calculations were performed on a subset of MOFs that ML predicted to exhibit high UV and UG capacities. These MOFs fall within the rectangular regions shown in Figures 6A and 6B, and exhibit capacities that meet or exceed 36 g-H2 L−1 and 7.5 wt % for PS conditions and 48 g-H2 L−1 and 7.5 wt % under TPS conditions. In total, 21,700 compounds were re-examined with GCMC based on their ML-predicted PS capacities, and another 7,901 were re-examined for TPS.

Figure 6C compares ML and GCMC predictions for usable capacities for 21,700 high-capacity MOFs under PS conditions. The strong overlap in the two datasets further highlights the accuracy of the ML models. A total of 8,187 MOFs were predicted by GCMC to out-perform PCN-610/NU-100 under these conditions. A summary of the 10 highest-capacity MOFs, sorted based on their GCMC capacities, is provided in Table 6 (a more extensive listing is provided in Table S12). The highest-capacity MOFs are all hypothetical compounds: five originate from the ToBaCCo database, two are from the University of Ottawa database, and the remainder are from the Northwestern database. These MOFs all exhibit high surface areas (average = 5,746, range = 4,346–7,835 m2 g−1) and large void fractions of 0.89, on average. The range of these property values are consistent with those reported in an earlier study,,, and suggest that maximizing the surface area is an important design guideline for PS operation. The highest-capacity MOF, mof_7642,59 is predicted to exhibit capacities of 11.1 wt % and 40.5 g-H2 L−1, surpassing that of PCN-610/NU-100, the record-holder under PS conditions. The crystal structure of mof_7642 is shown in Figure 7A.

Table 6

Highest-capacity MOFs, as identified by ML and verified with GCMC, under pressure swing and temperature + pressure swing conditions

Name	Source	Density (g cm−3)	Grav. surface area (m2 g−1)	Vol. surface area (m2 cm−3)	Void fraction	Pore volume (cm3 g−1)	Largest cavity diameter (Å)	Pore limiting diameter (Å)	Usable grav. capacity (wt %)		Usable vol. capacity (g-H2 L−1)
									GCMC	ML	GCMC	ML
Pressure swing
mof_7642	ToBaCCo	0.30	5,561	1,695	0.89	2.93	12.8	11.8	11.1	10.3	40.5	37.4
mof_7690	ToBaCCo	0.30	5,715	1,706	0.89	2.98	12.8	12.0	11.3	10.4	40.3	37.3
mof_7594	ToBaCCo	0.40	5,070	2,031	0.86	2.15	11.2	9.7	8.6	7.9	39.9	37.0
mof_7210	ToBaCCo	0.29	5,936	1,730	0.89	3.04	13.4	11.7	11.4	10.5	39.8	37.1
mof_7738	ToBaCCo	0.25	6,054	1,502	0.90	3.64	14.5	13.5	13.0	12.0	39.7	37.0
hypotheticalMOF_5045702_i_1_j_24_k_20_m_2	NW	0.31	5,926	1,820	0.88	2.87	16.0	11.0	10.9	10.1	39.7	37.2
str_m3_o19_o19_f0_nbo.sym.1.out	UO	0.31	5,073	1,583	0.90	2.88	17.7	12.9	10.8	10.1	39.7	37.1
hypotheticalMOF_5037315_i_1_j_20_k_12_m_1	NW	0.31	5,818	1,787	0.88	2.86	16.0	11.0	10.9	10.0	39.7	37.0
hypotheticalMOF_5037467_i_1_j_20_k_12_m_8	NW	0.31	5,860	1,800	0.88	2.85	16.0	11.0	10.9	10.0	39.7	37.0
str_m3_o5_o20_f0_nbo.sym.1.out	UO	0.39	4,772	1,882	0.87	2.22	14.1	9.6	8.7	8.1	39.7	37.2
Temperature + pressure swing
str_m1_o1_o11_f0_pcu.sym.102.out	UO	0.45	4,352	1,974	0.84	1.84	12.9	10.1	10.4	9.7	53.1	48.1
str_m1_o1_o11_f0_pcu.sym.117.out	UO	0.47	4,162	1,977	0.83	1.74	12.8	9.9	9.9	9.0	52.8	48.0
str_m1_o1_o11_f0_pcu.sym.121.out	UO	0.47	4,263	2,006	0.83	1.76	12.1	10.2	10.0	9.4	52.7	48.1
str_m1_o1_o11_f0_pcu.sym.13.out	UO	0.46	4,326	2,005	0.83	1.79	12.7	9.9	10.1	9.3	52.6	48.0
str_m1_o1_o11_f0_pcu.sym.159.out	UO	0.58	3,703	2,138	0.80	1.38	10.4	8.6	8.3	7.6	52.6	48.5
str_m1_o1_o11_f0_pcu.sym.200.out	UO	0.45	4,359	1,978	0.84	1.84	12.9	10.1	10.3	9.6	52.6	48.1
str_m1_o1_o11_f0_pcu.sym.212.out	UO	0.60	3,417	2,035	0.83	1.39	12.0	10.1	8.1	7.5	52.5	48.1
str_m1_o1_o11_f0_pcu.sym.51.out	UO	0.46	4,330	2,007	0.83	1.79	11.9	9.9	10.1	9.3	52.5	48.1
str_m1_o1_o11_f0_pcu.sym.71.out	UO	0.45	4,436	1,980	0.84	1.87	13.0	10.9	10.4	9.7	52.5	48.1
str_m1_o1_o11_f0_pcu.sym.89.out	UO	0.58	3,507	2,043	0.83	1.42	12.4	9.8	8.2	7.7	52.5	48.1

Here, NW and UO refer to the Northwestern and University of Ottawa databases. Grav., gravimetric; Vol., volumetric.

Figure 7

Crystal structures of high-capacity MOFs

Highest-capacity MOFs under (A) PS and (B) temperature + PS conditions. These MOFs originate from the ToBaCCo and University of Ottawa databases, respectively.

A search in the CCDC was performed to identify MOFs that have been synthesized that are similar to the high-capacity compounds identified here. The existence of similar MOFs may suggest synthetic procedures that could be adapted to the present systems. The top 5 MOFs under PS conditions contain relatively long tritopic linkers. In the case of mof_7642, this search identified the interpenetrated MOF RANCEQ as having a similar index of 0.82. Interpenetration is fairly common in MOFs (such as mof_7642) with longer linkers, and is generally undesirable for achieving high uptake. Nevertheless, several examples of successful synthesis of MOFs with long, multi-topic linkers that do not undergo interpenetration, have been reported. These include MOF-180 and MOF-200, the PCN-6X series, and NOTT-112. The next four PS candidates in Table 6 exhibit pillared Zn paddlewheel clusters with long ditopic linkers. Karagiaridi et al. demonstrated the feasibility of synthesizing pillared paddlewheel MOFs with long linkers; the SALEM-X series are examples. Finally, str_m3_o5_o20_f0_nbo.sym.1.out is based on a Zn paddlewheel cluster and a ditopic linker. HOFSUS (CSD Refcode) is an example of such a MOF.

Figure 6D provides a similar comparison between ML predictions and GCMC calculations for MOFs expected to exhibit high capacities under TPS conditions. Under these conditions, only 95 MOFs were predicted by GCMC to out-perform MOF-5. A summary of the 10 highest-capacity MOFs, sorted by their GCMC capacities, is provided in Table 6 (see Table S13 for a more extensive tabulation). As found for PS operation, all of the top performing candidates are hypothetical compounds. One difference with the PS case is that all of these MOFs originate from the University of Ottawa database. Furthermore, none of the highest-capacity MOFs identified for PS operation appear as top candidates for TPS. Comparing the highest-capacity MOFs for both operating conditions, it can be seen that the high-capacity TPS MOFs systematically exhibit lower surface areas (average = 4,073 m2 g−1), smaller void fractions (average = 0.83), and higher densities. Hence, the categories of MOFs that maximize uptake under PS and TPS conditions exhibit distinct properties. These differences suggest that maximizing the surface area—which, as discussed above, is desirable for maximizing PS capacity—is not advantageous for TPS operation. This behavior can be explained by trends in total capacities, which the TPS capacities reported here approximate. More specifically, it is known that total volumetric capacities are maximized for intermediate values of the surface area; for larger surface areas the volumetric capacity decreases.

Returning to the list of promising MOFs for TPS operation, Table 6 reports that the highest-capacity MOF, str_m1_o1_o11_f0_pcu.sym.102.out, has a GCMC-predicted capacity of 10.4 wt % and 53.1 g-H2 L−1. This capacity surpasses that of MOF-5, which, to our knowledge, holds the capacity record under these conditions. The crystal structure of this MOF is shown in Figure 7B.

The top 10 MOFs under TPS conditions contain the same Zn metal cluster and terephthalic acid linkers, where the linkers have been modified with varying functional groups. The slight differences in the capacities of these MOFs can be traced to differences in the functional groups. A similarity search based on str_m1_o1_o11_f0_pcu.sym.117.out identified 40 similar MOFs. Approximately 30 of these (for example, HIFTOG, MIBQAR, UNIGEE, VUSJUP, and ZELROZ) contain Zn metal clusters and linkers based on variants of terephthalic acid.

Figures S8 and S9 and Table S14 quantify the differences between ML and GCMC predictions on the subset of high-capacity MOFs shown in Figures 6C and 6D. For PS operation, the AUE of ML relative to GCMC is 0.24 wt % and 0.66 g-H2 L−1, while for TPS the AUE is 0.24 wt % and 1.28 g-H2 L−1. Both sets of errors are comparable with the errors reported in Table 4 for the original test set of MOFs. Figures S8C and S8F and S9(c,f) plot the frequency distribution of the differences between GCMC and ML. These distribution plots suggest that the largest differences occur for predictions involving real MOFs and for hypothetical MOFs extracted from databases other than those from Northwestern, University of Ottawa, and BJT. (These MOFs are referred to as “other hypothetical MOFs” in Figure 6). These MOFs, along with the real compounds, exhibit higher structural diversity than those contained in the other databases. For example, the diversity of the topologies used in the ToBaCCo and Zr-MOFs databases and in the linkers used in MTV-MOF database are larger than what is found in the databases from Northwestern, University of Ottawa, and BJT.

Discussion

Limitations of this study

As described previously, some of the high-capacity MOFs identified here may prove difficult to synthesize. Although this limitation applies primarily to the hypothetical MOFs, in some cases real MOFs are also known to undergo framework collapse during activation, which would reduce capacity., Nevertheless, future improvements to synthesis techniques may overcome these limitations—what is difficult to make today may be possible in the future. Secondly, our models do not distinguish between realistic MOFs having non-defective crystal structures and those for which the structures are defective/unrealistic. Unrealistic structures can result from incomplete or imperfect virtual solvent removal and the presence of partial occupancies or symmetry disorder in the crystal structure. Consequently, a defective/unrealistic MOF could be erroneously predicted to be a promising candidate. Follow-up calculations using GCMC and visual inspection of the crystal structure are recommended for all promising candidates identified by ML. Finally, the ML models developed here are non-interpretable, “black-box” models. Although these models are demonstrated to be highly accurate, additional effort is required to assess the relative importance of their input data. (The approach demonstrated here for evaluating feature importance involved the development of multiple models with varying numbers and combinations of features.) Alternatively, interpretable white-box ML models could be developed to provide more insight into feature importance. However, our experience suggests that white-box models generate less accurate predictions.

Concluding remarks

The H2 storage capacities of nearly a million MOFs have been predicted via ML. The predictions span a diverse collection of MOFs sourced from 19 databases and reveal performance under two operating conditions: PS and temperature + PS. More than a dozen ML algorithms were benchmarked, with the ERT method found to be the most accurate. The resulting ML models are accessible on the web at the HyMARC data hub. These models allow for accurate, rapid screening of the hydrogen storage properties of new MOFs using minimal structural data as input; only a single feature is needed for the simplest models.

The accuracy of the ML models was characterized as a function of training set size and the number/combination of input features. Regarding the dependence on the training set, the accuracy of the models can be well described using a simple power law function of the training set size. The dependence on the number and combination of input features was determined by evaluating 508 independent ML models generated from all possible combinations of the seven features. The most important features for predicting H2 uptake are pore volume (for gravimetric capacity) and void fraction (for volumetric capacity).

Using these models, 8,282 MOFs are identified that have the potential to exceed the capacities of state-of-the-art materials under usable conditions. The identified MOFs are predominantly hypothetical compounds, which (for PS operation) exhibit low densities (<0.31 g cm−3) in combination with high surface areas (>5,300 m2 g−1), void fractions (∼0.90), and pore volumes (>3.3 cm3 g−1). These MOFs are suggested as targets for experimental synthesis.

Experimental procedures

Resource availability

Lead contact

Prof. Donald Siegel, djsiege@umich.edu.

Materials availability

This study did not generate new reagents.

Data and code availability

Original data have been deposited to HyMARC Data Hub: https://datahub.hymarc.org/dataset/computational-prediction-of-hydrogen-storage-capacities-in-mofs

Interactive ML models: https://sorbent-ml.hymarc.org/

Python API: https://sorbent-ml.hymarc.org/

Downloadable ML models and instructions: https://doi.org/10.6084/m9.figshare.14173520.v1