Computational insights into hydrogen adsorption energies on medium-entropy oxides

Abstract

High entropy oxides (HEOs) have emerged as promising catalysts for several important chemical transformations including alkane activation. Hydrogen adsorption energy (HAE) has been used as a key descriptor for many reactions including methane C–H activation and hydrogen evolution reactions. Hence, understanding the relationship between HAEs and the surface chemistry of HEO surfaces could lay the foundation for meaningful correlations among methane C–H activation, HAE, and the complex, local environment of HEO surfaces. Here, we used a medium-entropy oxide as a prototypical system – Mg_0.25Ni_0.25Cu_0.25Zn_0.25O with a rock-salt structure – to interrogate these relationships. We sampled 2000 different surfaces of its (100) plane and calculated the HAEs at randomly chosen surface O sites using density functional theory (DFT). Our analysis of the 2000 data points reveals that the HAEs at the surface O sites are significantly influenced by the local environment around the adsorption sites, particularly the nature of the metal atom directly below the surface O site where H adsorbs. After comparing several popular graph-neural-network-based machine learning models, we found that the DimeNet++ model performed best achieving satisfactory accuracy in predicting HAEs for both Mg_0.25Ni_0.25Cu_0.25Zn_0.25O and slightly varied compositions. Our work underscores the promise of such models and the need for further refinement to address the complexity of HEOs.

---

1. Introduction

High-entropy oxides (HEOs),^1–3 a class of complex oxides composed of five or more principal metal cations arranged in a single-phase crystal structure, have garnered increasing attention due to their compositional diversity and disorder that enhance thermal stability and open up new opportunities in catalysis.^4–6 They are being actively studied as catalysts for reactions such as oxygen evolution^7–10 and methane combustion¹¹ due to their tunable compositions and properties. Moreover, HEOs exhibit intriguing exsolution behavior under reducing conditions such as in dry reforming of methane (DRM).^12,13 In addition, they can be employed as supports to stabilize single-atom catalysts.^14,15

Despite the configurational complexities of HEOs, computational and machine-learning approaches have been used to provide useful insights into their structure–function relationships. For example, Rossmeisl and coworkers have proposed a theoretical framework¹⁶ where they fit density functional theory (DFT) calculations performed for a limited number of sites to a model that predicts the reaction energies for all possible local atomic environments and applied it to HEOs for oxygen reduction reactions (ORRs).^9,17 This model worked well for ORRs on high-entropy alloys and HEOs partly because the key intermediates of ORRs such as *OH adsorb on the surface metal sites and the metal-element identity is the main factor dictating the adsorption energy. However, whether such an approach can apply to methane activation is unclear because intermediates such as *H will be adsorbed on O sites with more complex local environments due to the different combinations of metal cations around it.

A composition–activity relationship for HEOs for methane activation would be highly desirable for advancing more efficient HEO-based catalytic routes for methane and light alkane chemistries such as methane combustion¹¹ and photoelectrocatalytic coupling of methane to acetic acid.¹⁸ Although hydrogen adsorption energy (HAE) has been found to be a reliable descriptor for C–H activation energy on many oxides and other materials,^19,20 it is unclear whether it can be applied to HEOs. As a first step toward a composition–activity relationship for HEOs in methane activation, we will use a medium entropy oxide (MEO) containing only four principal metal cations to investigate if one can accurately predict HAEs with a machine-learning surrogate model using DFT data to train it.

2. Methods

2.1. Structure models and density functional theory computation

The overall workflow is shown in Fig. 1a. We chose the rock-salt Mg_0.25Ni_0.25Cu_0.25Zn_0.25O medium-entropy oxide (MEO) because of its bulk and surface structure simplicity; in addition, it was previously synthesized experimentally.²¹ To construct the initial structure of rock-salt Mg_0.25Ni_0.25Cu_0.25Zn_0.25O MEO, we created a periodic slab with a (2 × 2) cubic MgO supercell first of the (100) facet. We chose the (100) surface because it is the most stable, non-polar, and commonly exposed facet in rocksalt oxides such as MgO, NiO, and ZnO under typical experimental conditions.^22–24 The slab cell consists of five layers with a total of 80 atoms and a vacuum thickness of 13 Å and has the lattice parameters of a = b = 8.44 Å and c = 22.22 Å. Here, we built our surface model using a bulk lattice parameter of 4.22 Å, based on the experimental X-ray diffraction data of rock-salt HEOs containing Mg, Co, Ni, Cu, and Zn.²¹ The Brillouin zone was sampled with the 3 × 3 × 1 k-mesh. The atoms in the bottom two layers were frozen to their bulk position and only the top three layers were allowed to relax. To get the composition of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O, we randomly replaced 75% of the Mg atoms in the slab with equal numbers of Ni, Cu, and Zn atoms. All the randomization processes were performed using the “random” module in Python.²⁵ We created 2000 unique (100) surface structures by randomizing the positions of metal elements within the MEO structures while ensuring there were no duplicate structures across the generated dataset.²⁶ Similarly, to ensure the diversity of the local environment, the site of H adsorption was also completely randomly selected (only one surface O site was selected for each surface).

Fig. 1 Modeling hydrogen adsorption on the (100) surfaces of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O: (a) workflow; (b) a random (100) structure (the unit cell has 80 atoms in total) and H adsorption on a randomly chosen O site on the (100) surface (H coverage: 1.4 atom per nm²); (c) local environment of the hydrogen adsorption site: positions 1–4 denote the nearest metal neighbors of the O site on the surface, while position 5 refers to the subsurface metal atom located directly beneath the O site.

Spin-polarized DFT calculations were performed using the Vienna ab initio Simulation Package (VASP)^27,28 version 6.3.2. The Perdew–Burke–Ernzerhof (PBE)²⁹ functional of the generalized-gradient approximation (GGA) was used for electron exchange and correlation. The electron–core interaction was described using the projector-augmented wave method (PAW).^30,31 The kinetic energy cutoff was set to 500 eV for the plane wave basis set. Grimme's DFT+D3 method was used to account for van der Waals interactions.³² To better describe transition metal oxides that have strongly localized 3d electrons, DFT+U is often used. To assess the impact without including the U term in our systems, we conducted benchmark calculations on a representative set of configurations and found that applying DFT+U (with typical U values for Ni and Cu) alters the hydrogen adsorption energies by less than 0.15 eV compared to GGA-PBE. Given that our study primarily focuses on uncovering relative trends and structure–property correlations across 2000 surface configurations, we consider this level of deviation acceptable. Moreover, a recent study showed that in correlating catalytic activity with an electronic descriptor for perovskite oxides, the PBE functional performs better than DFT+U and the hybrid functionals.³³

2.2. Machine learning methods

Various ML models based on graph neural networks (GNNs) were employed and compared to predict HAEs using the MEO surface structure as input, including the Directional Message Passing Neural Network (DimeNet and DimeNet++),^34,35 MatErials Graph Network (MEGNet),³⁶ SchNet,³⁷ GemNet,^38,39 and Atomistic line graph neural network (ALIGNN).⁴⁰ These models were implemented in PyTorch⁴¹ and some of them were integrated in the FAIR-Chem framework.^42,43 After initial testing, DimeNet++ was selected for its accuracy and computational efficiency. The model architecture comprises four interaction blocks with a hidden dimension of 256. Geometric information was encoded using spherical harmonics up to order 7 and six radial basis functions with a neighbor cutoff radius of 8 Å. Graph construction was performed on-the-fly during training. The model was trained to minimize the mean absolute error (MAE) between predicted and reference adsorption energies. An initial learning rate of 1 × 10⁻⁴ was used and adjusted adaptively. To prevent overfitting, early stopping was employed with a patience of 20 epochs. The model was trained on the (100) surface configurations of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O, with the dataset randomly split into 80% for training, 10% for validation, and 10% for testing. It was evaluated on the designated test set and further tested on additional (100) surfaces with slightly different compositions (Mg_0.3Ni_0.2Cu_0.25Zn_0.25O and Mg_0.3Ni_0.2Cu_0.2Zn_0.3O).

3. Results and discussion

The key challenges in establishing structure–activity relationships for HEOs is how to statistically consider both the bulk configurations and the surface configurations. Below, we first present the DFT energies and their distributions to understand the data from the 2000 surface structure samples with a fixed bulk composition and then apply the GNN ML approach to the data and test it on sligthly varied bulk compositions.

3.1. DFT energetics and distribution

The surface structure, H adsorption site, and its local environment for a randomly generated (100) surface of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O MEO (shortened as MEO) of the rock-salt structure are shown in Fig. 1b and c. The H atom prefers to be adsorbed atop a surface O site, which is coordinated by four metal cations on the surface (labelled 1 to 4) and another metal cation directly beneath (labelled 5). We also tested other adsorption sites for H including the metal-top and metal–oxygen bridge sites on 10 surfaces and found that they are less favorable by at least ∼0.2 eV than the O-top site for H adsorption. To get a general sense of the variation and distribution of the reactivity of surface O sites toward H adsorption, we analyzed the total energy distributions for the clean MEO (100) surfaces (Fig. 2a) and the surface after hydrogen adsorption (Fig. 2b) relative to their respective average values, as well as the hydrogen adsorption energies (Fig. 2c). One can see that they all follow a normal distribution, with the total energies having a range of about 4 eV and the HAEs having a range of about 2 eV.

Fig. 2 Distributions of DFT-computed energies of the 2000 Mg_0.25Ni_0.25Cu_0.25Zn_0.25O (100) surfaces: (a) energy of the clean surface relative to the mean (eV per cell; each cell has 80 atoms); (b) energy of the H-adsorbed surface relative to the mean (eV per cell; each cell has 80 atoms plus 1 H atom); (c) hydrogen adsorption energy (eV, in reference to ½ H₂; H coverage: 1.4 nm⁻²).

3.2. Analysis of the clean surfaces and their energies

To understand the broad distribution of the clean surface energies shown in Fig. 2a, we have analyzed the 2000 random MEO (100) surfaces in detail by examining how their energies vary with the surface-layer concentrations of the four different types of metal elements. As shown in Fig. 3, the surfaces with fewer Ni and Mg atoms (Fig. 3a and b) and more Cu atoms (Fig. 3d) in the top layer are more stable, while the energies of the surfaces are less sensitive to the surface Zn coverage (Fig. 3c). This trend is further reflected in the plot of Pearson correlation coefficients between the energy of the surface and the concentration of a specific metal atom in the surface layer (Fig. 4a) and consistent with the computed (100) surface energies of the individual oxides (with the same rock-salt structure), which follow the order of NiO > MgO > ZnO > CuO (Fig. 4b). In other words, CuO has the lowest surface energy (Fig. 4b), so the Cu ions are more likely to be on the surface, resulting in an overall lower surface energy and more negative Pearson values for Cu, as shown in Fig. 4a. NiO has the highest surface energy (Fig. 4b), so the Ni ion tends to stay below the surface layer: the more Ni on the surface, the higher the surface energy, exhibiting more positive Pearson correlation for Ni, as shown in Fig. 4a.

Fig. 3 Distributions of the relative energies (eV per cell) of the 2000 clean (100) surfaces of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O for different surface concentrations of metal cations: (a) Ni; (b) Mg; (c) Zn; (d) Cu.

Fig. 4 (a) Pearson correlation coefficients between the energy of the surface and the surface-layer concentration of a metal atom for the 2000 Mg_0.25Ni_0.25Cu_0.25Zn_0.25O (100) surfaces; (b) computed (100) surface energies (J m⁻²) of different rock-salt metal oxides.

3.3. Analysis of the H-adsorbed surfaces and their energies

Next, we want to find out how the surface energies change after H adsorption. We first analyzed the correlation of the surface energies before and after H adsorption and found that there is a strong correlation: on average, the more stable clean surfaces retain their relative stability after H adsorption, as can be seen in Fig. 5. On one hand, it appears that the dominant metal cation on the surface has a strong influence on the energy of the H-adsorbed surface. On the other hand, the local environment of the surface O site where H adsorbs likely plays a significant role. In Fig. 6, we have plotted the relative energies of the 2000 H-adsorbed surfaces of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O in the form of a heat map, as a function of both the position-5 element (just beneath the O site) and the dominant surface element. It confirms the dominant surface metal type being the major factor, while it also shows that when Cu is in position-5, it further stabilizes the H-adsorbed surface no matter what the dominant surface-layer metal is. In other words, the most stable H-adsorbed surfaces have more Cu atoms on the surface as well as a Cu atom just beneath the O site where H adsorbs.

Fig. 5 Correlation between the relative energies of clean surfaces and H-adsorbed surfaces for the 2000 (100) surfaces of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O. R² = 0.85.

Fig. 6 Average relative energies (in eV per cell) of the 2000 H-adsorbed surfaces of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O as a function of the position-5 element (y-axis) and the dominant surface element (x-axis). Position-5 refers to the metal element directly beneath the O site where H adsorbs (see Fig. 1c).

3.4. Analysis of the DFT HAEs

To identify the factors influencing HAEs on the 2000 MEO (100) surfaces, we employ a divide-and-conquer approach by uniformly dividing the adsorption data into three subsets according to the adsorption strength based on the distribution in Fig. 2c: strong adsorption (HAE from −1.30 to −0.60 eV), intermediate adsorption (HAE from −0.60 to −0.39 eV), and weak adsorption (HAE from −0.39 to +0.60 eV). We found that the HAE has strong dependence on the positon-5 metal element. As shown in Fig. 7a, the structures with the strongest H adsorption tend to have Cu at position-5, while the structures with the weakest H adsorption tend to have Mg at position-5. In Fig. 7b, we sliced the data in a different way by plotting the HAE distribution for a given metal at position-5 which agrees with the trend in Fig. 7a: Cu at position-5 tends to have more negative HAEs and stronger H adsorption and Mg at position-5 tends to have less negative HAEs and weaker H adsorption, while Ni and Zn have similar intermediate H-adsorption values. On the other hand, we found that positions 1–4, corresponding to the metal sites surrounding the surface O site, have much less influence on the HAE values (Fig. 8).

Fig. 7 Relationship between the H adsorption strength and the metal element: (a) Percent distribution among the four metal elements at position 5 for different HAE ranges; (b) distribution of HAEs for a given metal element at position 5 which refers to the metal site directly beneath the O site where H adsorbs (see Fig. 1c).

Fig. 8 Percent distribution among the four metal elements at positions 1–4 for different HAE ranges. Positions 1–4 refer to the four nearest metal neighbors of the O site on the surface (see Fig. 1c).

We think that the sensitivity of the HAEs on the position-5 metal is closely related to the local environment of the O site on the (100) surface of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O. As shown in Fig. 1c in the side view, the nearest five metal neighbors of the surface O site form an inverted pyramid shape with the position-5 metal at the apex right beneath the O site. After hydrogen adsorption, the five metal atoms and the H atom on top form an octahedral coordination shell around the O site with the position-5 metal and the H atom axially coordinated to the O site. As a result, the adsorbed H atom is in a tug-of-war with the position-5 metal in bonding with the O site and the O–H strength is more sensitive to the push-and-pull with the position-5 metal, while the four equatorial metal atoms on the surface tend to balance out their interactions with the O site.

3.5. Machine learning models

The above analysis is based on DFT sampling of 2000 surface structures of a single composition for a medium-entropy oxide. To fully address the compositional variation and the resulting surface complexity of HEOs, we need an ML surrogate model. Here, we have evaluated five popular graph-neural-network-based ML models using the 2000 surface structures as input to predict HAEs: DimeNet++,^34,35 MEGNet,³⁶ SchNet,³⁷ GemNet,^38,39 and ALIGNN.⁴⁰ As shown in Fig. 9a, DimeNet++ demonstrated the best predictive performance for HAEs, followed closely by ALIGNN. DimeNet++ also exhibited superior predictive accuracy on an additional test set with compositional variations, demonstrating greater robustness and transferability. The good performance of DimeNet++ (Fig. 9b) can be attributed to its strong emphasis on local atomic environments, enabling effective capture of structure–property relationships.^34,35 The normal distribution of HAEs (Fig. 2c) places most of the data points in the middle and limits the model accuracy for very strong and very weak adsorption at the tails of the distribution (Fig. 9b), indicating the need for further improvement in HAE predictions.

Fig. 9 (a) Comparison of different machine-learning models trained on DFT-computed HAEs on 2000 Mg_0.25Ni_0.25Cu_0.25Zn_0.25O (100) surfaces: left axis, coefficient of determination; right axis, mean absolute error (MAE). (b) Parity plot for the DimeNet++ based on the 10% held-out test set.

To assess the predictive power of DimeNet++, we further tested it on two different compositions, Mg_0.3Ni_0.2Cu_0.25Zn_0.25O and Mg_0.3Ni_0.2Cu_0.2Zn_0.3O, from the composition of Mg_0.25Ni_0.25Cu_0.25Zn_0.25O used for training. As shown in Fig. 10, the model retained consistent predictive performance across these compositions, with an R² of 0.59 and a MAE of 0.12 eV. This consistency underscores DimeNet++'s generalization capability and potential as an ML surrogate model for predicting HAEs for HEOs. Looking ahead, we need to increase the complexity of the compositions to go beyond the four metal elements to include five or more metal elements, in order to get closer to the real experimental compositions for HEOs. This will be computationally challenging, but starting with five elements and leveraging the data efficiency of the DimeNet++ ML approach might prove fruitful.

Fig. 10 Performance of DimeNet++ for predicting the hydrogen adsorption energy (HAE) on the (100) surfaces of two different compositions, Mg_0.3Ni_0.2Cu_0.25Zn_0.25O and Mg_0.3Ni_0.2Cu_0.2Zn_0.3O, from Mg_0.25Ni_0.25Cu_0.25Zn_0.25O used for training.

The MAE of 0.12 eV from DimeNet++ is not bad, given that the intrinsic uncertainties in the DFT training data of adsorption energetics are often on the order of 0.1 to 0.2 eV (due to the choices of the kinetic energy cutoff, k-point sampling, pseudopotentials, convergence criteria, etc.). On the other hand, further improvement is needed in raising the R² value. Potential directions include improving the accuracy at the DFT level by using a higher kinetic energy cutoff, denser k-point sampling, and tighter geometry convergence criterion; sampling more configurations and diverse compositions at the DFT level; incorporating long-range electrostatic descriptors or electronic features (e.g., Bader charges and the local density of states) in the ML models.

4. Conclusions

We have investigated the hydrogen adsorption energies (HAEs) on the (100) surfaces of a prototypical medium-entropy oxide, Mg_0.25Ni_0.25Cu_0.25Zn_0.25O, through DFT calculations and machine learning (ML) methods. By sampling 2000 randomized surface configurations, we revealed that the HAE strongly depends on the local atomic environment, particularly the metal atom directly beneath the oxygen adsorption site, with copper significantly enhancing adsorption strength. In contrast, the identity of the surrounding surface metal sites had a relatively minor influence. Machine learning models, especially DimeNet++, effectively captured these relationships, achieving moderate predictive accuracy and demonstrating the ability to generalize across slightly varied compositions (±0.05). Despite these promising results, further refinement in modeling approaches and additional analysis of local atomic descriptors are necessary to increase prediction accuracy for wider HAE and composition ranges. Future studies will explore optimal compositional ratios for the targeted design of high-entropy oxide catalysts with tailored adsorption properties for reactions such as methane combustion.

Data availability

The data supporting this article have been deposited in the ioChem-BD Computational Chemistry repository (DOI: https://doi.org/10.19061/iochem-bd-6-522).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was sponsored by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, Catalysis Science Program. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under contract no. DE-AC02-05CH11231.

References

1. C. M. Rost, E. Sachet, T. Borman, A. Moballegh, E. C. Dickey, D. Hou, J. L. Jones, S. Curtarolo and J. P. Maria, Nat. Commun., 2015, 6, 8485.
2. N. Dragoe and D. Bérardan, Sci., 2019, 366, 573–574.
3. C. Oses, C. Toher and S. Curtarolo, Nat. Rev. Mater., 2020, 5, 295–309.
4. Y. Sun and S. Dai, Sci. Adv., 2021, 7, eabg1600.
5. Y. Pan, J.-X. Liu, T.-Z. Tu, W. Wang and G.-J. Zhang, Chem. Eng. J., 2023, 451, 138659.
6. Y. Wang, J. Mi and Z.-S. Wu, Chem Catal., 2022, 2, 1624–1656.
7. M. V. Kante, M. L. Weber, S. Ni, I. C. G. van den Bosch, E. van der Minne, L. Heymann, L. J. Falling, N. Gauquelin, M. Tsvetanova, D. M. Cunha, G. Koster, F. Gunkel, S. Nemšák, H. Hahn, L. Velasco Estrada and C. Baeumer, ACS Nano, 2023, 17, 5329–5339.
8. J. Baek, M. D. Hossain, P. Mukherjee, J. Lee, K. T. Winther, J. Leem, Y. Jiang, W. C. Chueh, M. Bajdich and X. Zheng, Nat. Commun., 2023, 14, 5936.
9. V. A. Mints, K. L. Svane, J. Rossmeisl and M. Arenz, ACS Catal., 2024, 14, 6936–6944.
10. K. Miao, W. Jiang, Z. Chen, Y. Luo, D. Xiang, C. Wang and X. Kang, Adv. Mater., 2024, 36, e2308490.
11. T. Li, Y. Yao, Z. Huang, P. Xie, Z. Liu, M. Yang, J. Gao, K. Zeng, A. H. Brozena, G. Pastel, M. Jiao, Q. Dong, J. Dai, S. Li, H. Zong, M. Chi, J. Luo, Y. Mo, G. Wang, C. Wang, R. Shahbazian-Yassar and L. Hu, Nat. Catal., 2021, 4, 62–70.
12. M. Zhang, J. Ye, Y. Gao, X. Duan, J. Zhao, S. Zhang, X. Lu, K. Luo, Q. Wang, Q. Niu, P. Zhang and S. Dai, ACS Nano, 2024, 18, 1449–1463.
13. Y. Shao, C. Wu, S. Xi, P. Tan, X. Wu, S. Saqline and W. Liu, Appl. Catal., B, 2024, 355, 124191.
14. M. Li, K. Michael Siniard, D. M. Driscoll, A. S. Ivanov, X. Lu, H. Chen, J. Zhang, F. Polo-Garzon, Z. Yang and S. Dai, J. Catal., 2024, 437, 115645.
15. H. Xu, Z. Zhang, J. Liu, C.-L. Do-Thanh, H. Chen, S. Xu, Q. Lin, Y. Jiao, J. Wang, Y. Wang, Y. Chen and S. Dai, Nat. Commun., 2020, 11, 3908.
16. T. A. A. Batchelor, J. K. Pedersen, S. H. Winther, I. E. Castelli, K. W. Jacobsen and J. Rossmeisl, Joule, 2019, 3, 834–845.
17. K. L. Svane and J. Rossmeisl, Angew. Chem., 2022, 61, e202201146.
18. S. Nie, L. Wu, Q. Zhang, Y. Huang, Q. Liu and X. Wang, Nat. Commun., 2024, 15, 6669.
19. A. A. Latimer, A. R. Kulkarni, H. Aljama, J. H. Montoya, J. S. Yoo, C. Tsai, F. Abild-Pedersen, F. Studt and J. K. Norskov, Nat. Mater., 2017, 16, 225–229.
20. V. Fung, G. Hu, Z. Wu and D. E. Jiang, J. Phys. Chem. Lett., 2020, 11, 7049–7057.
21. W. Mnasri, D. Bérardan, S. Tusseau-Nenez, T. Gacoin, I. Maurin and N. Dragoe, J. Mater. Chem. C, 2021, 9, 15121–15131.
22. R. S. Koster, C. M. Fang, M. Dijkstra, A. van Blaaderen and M. A. van Huis, J. Phys. Chem. C, 2015, 119, 5648–5656.
23. G. Pacchioni and H. Freund, Chem. Rev., 2013, 113, 4035–4072.
24. J. Goniakowski, F. Finocchi and C. Noguera, Rep. Prog. Phys., 2008, 71, 016501.
25. Python Software Foundation, random — Generate pseudo-random numbers, In Python 3.8.10 documentation, 2020, https://docs.python.org/3.8/library/random.html.
26. The dataset of computational results from this work is available in the ioChem-BD repository and can be accessed via  DOI:10.19061/iochem-bd-6-522.
27. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
28. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
29. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
30. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
31. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
32. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
33. R. Jacobs, J. Hwang, Y. Shao-Horn and D. Morgan, Chem. Mater., 2019, 31, 785–797.
34. J. Gasteiger, J. Groß and S. Günnemann, International Conference on Learning Representations (ICLR), 2020.
35. J. Gasteiger, S. Giri, J. T. Margraf and S. Günnemann, Machine Learning for Molecules Workshop, NeurIPS, 2020.
36. C. Chen, W. Ye, Y. Zuo, C. Zheng and S. P. Ong, Chem. Mater., 2019, 31, 3564–3572.
37. K. T. Schutt, F. Arbabzadah, S. Chmiela, K. R. Muller and A. Tkatchenko, Nat. Commun., 2017, 8, 13890.
38. J. Gasteiger, F. Becker and S. Günnemann, Conference on Neural Information Processing Systems (NeurIPS), 2021.
39. S. Stocker, J. Gasteiger, F. Becker, S. Günnemann and J. T. Margraf, Mach. Learn.: Sci. Technol., 2022, 3, 045010.
40. K. Choudhary and B. DeCost, npj Comput. Mater., 2021, 7, 185.
41. A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein and L. Antiga, Advances in neural information processing systems 32 (NeurIPS 2019), 2019.
42. L. Chanussot, A. Das, S. Goyal, T. Lavril, M. Shuaibi, M. Riviere, K. Tran, J. Heras-Domingo, C. Ho, W. Hu, A. Palizhati, A. Sriram, B. Wood, J. Yoon, D. Parikh, C. L. Zitnick and Z. Ulissi, ACS Catal., 2021, 11, 6059–6072.
43. R. Tran, J. Lan, M. Shuaibi, B. M. Wood, S. Goyal, A. Das, J. Heras-Domingo, A. Kolluru, A. Rizvi, N. Shoghi, A. Sriram, F. Therrien, J. Abed, O. Voznyy, E. H. Sargent, Z. Ulissi and C. L. Zitnick, ACS Catal., 2023, 13, 3066–3084.

---