The relationship between activated H2 bond length and adsorption
  distance on MXenes identified with graph neural network and resonating
  valence bond theory
Motivated by the recent experimental study on hydrogen storage in MXene
multilayers [Nature Nanotechnol. 2021, 16, 331], for the first time we propose
a workflow to computationally screen 23,857 compounds of MXene to explore the
general relation between the activated H2 bond length and adsorption distance.
By using density functional theory (DFT), we generate a dataset to investigate
the adsorption geometries of hydrogen on MXenes, based on which we train
physics-informed atomistic line graph neural networks (ALIGNNs) to predict
adsorption parameters. To fit the results, we further derived a formula that
quantitatively reproduces the dependence of H2 bond length on the adsorption
distance from MXenes within the framework of Pauling's resonating valence bond
(RVB) theory, revealing the impact of transition metal's ligancy and valence on
activating dihydrogen in H2 storage.
Jiewei Cheng, Tingwei Li, Yongyi Wang, Ahmed H. Ati, Qiang Sun
The relationship between activated H2 bond length and adsorption distance on MXenes identified with graph neural network and resonating valence bond theory

###### Abstract

Motivated by the recent experimental study on hydrogen storage in MXene multilayers [_Nature Nanotechnol._ 2021, 16, 331], for the first time we propose a workflow to computationally screen 23,857 compounds of MXene to explore the general relation between the activated H2 bond length and adsorption distance. By using density functional theory (DFT), we generate a dataset to investigate the adsorption geometries of hydrogen on MXenes, based on which we train physics-informed atomistic line graph neural networks (ALIGNNs) to predict adsorption parameters. To fit the results, we further derived a formula that quantitatively reproduces the dependence of H2 bond length on the adsorption distance from MXenes within the framework of Pauling's resonating valence bond (RVB) theory, revealing the impact of transition metal's legacy and valence on activating dihydrogen in H2 storage.

## Introduction

One of the most promising solutions to achieving carbon neutrality is to use hydrogen energy, which has renewability and zero CO2 emission[1, 2, 3]. However, the most difficult challenge is to find materials that can store hydrogen with large gravimetric and volumetric density and operate under ambient thermodynamic conditions[4, 5, 6]. While the bonding of hydrogen existed in nature is either too strong as in metal hydrides, or too weak as in MOFs[7]. To balance the thermodynamics and kinetics, hydrogen binding needs to be between physisorption and chemisorption, namely, in quasi-molecular form, where H2 is activated with elongated H-H bond length.

The activation can be achieved either by electron transferring like in Kubas effect[8] or by charge polarization unveiled by Jena and co-workers[9]. The former takes advantage of the unfilled \(d\) orbitals of transition-metal atoms where H2 molecules donate electrons to the unfilled \(d\) orbitals and the transition metals back donate the electrons to the H2 molecules, and the latter relies on the charge polarization induced by the local electrical field. In both cases, H2 retains its molecular bond but with stretched H-H bond length.

Figure 1: Workflow for the high-throughput screening of MXene for hydrogen storage via graph neural network and multiscale simulation.
The H-H bond length is like a barometer that directly measures the interaction strength of H\({}_{2}\) molecule with materials, while the interaction directly depends on the distance of H\({}_{2}\) from the substrates. Therefore, it is of vital significance to explore the relation of H-H bond length changing with the adsorption distance, which motivates the studies in this direction. For example, Grundemann _et al._[10] developed an equation based on Pauling's bond order concept and experimental evidences to explain the H-H bond length and metal-H distances' correlations in transition metal dihydrides and dihydrogen complexes. Zhou _et al._[11] computationally studied external electric field's impact on stretching dihydrogen bond and adsorption distance on boron nitride sheet. Friederich _et al._[12] revealed the correlation between H-H bond length and activation process of dihydrogen in Vaska's complexes via _ab-initio_ simulations and machine learning. These studies clearly showed that the pursued relationship depends on the substrate materials.

Recently, MXenes have been found to be promising candidates both as hydrogen storage materials[13, 14, 15] and as catalysts for enhancing metal hydrides' performance[16, 17, 18] due to its superior ability in activating H\({}_{2}\) bonds. However, current studies have only covered a very small part of the MXene compounds and the underlying hydrogen activation mechanisms on MXenes are still unclear. Intriguingly, the transition metal (TM) sites are normally coordinated in MXenes, while one functional group is bonded to three TM sites. Therefore, the dihydrogen bond activation process is neither on open metal sites nor on normally saturated TM sites. Two-dimensional (2D) MXene monolayers are building blocks of multilayer MXenes. It is of importance to explore the activation of dihydrogen on 2D MXene layers, especially, the general relationship between H-H bond length and the adsorption distance from MXenes.

In this work, we design a workflow by combining high-throughput simulation and machine learning to screen the MXenes' material space for hydrogen activation, as shown in Figure 1, where the state-of-art ALIGNN model[19] and physics-informed machine learning[20] are applied.

**Results and discussion**

**MXene's hydrogen adsorption geometry datasets**

According to experiment[14] and aNANt database[21], the MXenes' H\({}_{2}\) adsorption geometry dataset is constructed. MXenes' H\({}_{2}\) adsorption geometry refers to H\({}_{2}\)'s bond length and its adsorption distance from MXene. In a unit cell of MXene (see Figure 1), five joint layers are arranged in an array of FG\({}_{1}\), TM\({}_{1}\), C/N, TM\({}_{2}\) and FG\({}_{2}\) (FG\({}_{1/2}\): functional group, TM\({}_{1/2}\): early transition metal, C/N: carbon or nitrogen). According to aNANt database[22], there are 11 and 13 choices for TM sites and FG layers respectively, which eventually generates 23,857 distinct computationally designed 2D MXenes monolayers. To avoid jeopardizing hydrogen storage performance, MXenes with fifth period transition metals (i.e. Hf, Ta and W) are excluded due to their heavy atomic weight. The remaining 12,647 structures are randomly divided into 1,305 and 11,342 subsets (Figure S2) as the training and prediction subsets. We assign one hydrogen molecule on each TM site and conducted high-throughput calculation to obtain hydrogen adsorption geometry for each MXene monolayer in the first subset. More details of our simulation can be found in Supplementary Information.

We then extract four labels to identify the adsorption geometry of H\({}_{2}\) on MXenes: the bond lengths of two hydrogen molecules ( \(d_{\text{H1}}\), \(d_{\text{H2}}\)), the distances between TM sites and geometry centers of hydrogen molecules (d\({}_{\text{MH1}}\), d\({}_{\text{MH2}}\)). The H\({}_{2}\) bond length and adsorption distance can faithfully reflect the process of H\({}_{2}\)'s adsorption, as is illustrated in Figure 2. A potential energy surface scan is performed to H\({}_{2}\) adsorbed on a \(3\times 3\times 1\) supercell of Ti\({}_{2}\)CH\({}_{2}\) (Figure 2(a-b)), where H\({}_{2}\)'s z coordinate and sorbent atoms are fixed and H\({}_{2}\)'s bond length is able to relax. Figure 2(c) shows that when H\({}_{2}\) approaches the sorbent, the \(d_{\text{H}}\) scarcely stretches and adsorption energy E\({}_{ads}\) increases slightly due to charge polarization. When \(d_{\text{MH}}\) is smaller than 4 A, H-H bond and E\({}_{ads}\) observably increases with the decrease of adsorption distance

Figure 3: Distribution of distances between (a) two hydrogen atoms and (b) TM and center of H\({}_{2}\). PI-ALIGNN predictions versus DFT calculation results of (c) \(d_{\text{H}}\) and (d) \(d_{\text{MH}}\) in the test set. The data related to up/bottom hydrogen molecule are plotted with aqua/light-blue color.

Figure 2: (a) side view and (b) top view of H\({}_{2}\) adsorbed on a \(3\times 3\times 1\) supercell of Ti\({}_{2}\)CH\({}_{2}\); (c) adsorption energy of H\({}_{2}\) changing with \(d_{\text{MH}}\) and \(d_{\text{H}}\) on Ti\({}_{2}\)CH\({}_{2}\).
primarily due to Kubas effect. Therefore, the four adsorption geometry labels are sufficient to give an insight into Hz's asorption on MXenes.

The distributions of these four labels in our dataset are plotted in Figure 3 (a-b) and the relation between \((\,d_{\rm H1}\,,d_{\rm H2})\) and \((\,d_{\rm MH1}\,,d_{\rm MH2})\) is plotted in Figure S1. Figure 3 (a) shows that the molecular or quasi-molecular adsorption of H\({}_{2}\) happens with \(d_{\rm H}\) close to the equilibrium distance and atomic H binding also happens.

**Two-step ALIGNN model**

The first ALIGNN model as H\({}_{2}\) quasi-molecule classifier is designed to predict whether the hydrogen maintains a molecule or quasi-molecule form, since the adsorbed H\({}_{2}\) can be dissociated in many ways (Figure S3), which will reduce the accuracy of H\({}_{2}\) adsorption geometry predictor due to limited amount of data. The model is composed of two identical ALIGNNs in the classification form to predict the state of one hydrogen molecule. If the distance between hydrogen atoms in the adsorbed molecule is smaller than 0.9 A (1.2 times of the equilibrium distance), the MXene structure will be labeled quasi-molecular configuration. The precision of our two models on test set reaches 0.938 and 0.957, respectively (Table S2).

Next, 1,021 structures in our dataset with two quasi-molecule labels are picked out to train the second model, a H\({}_{2}\) adsorption geometry predictor. This model is a multi-output regression ALIGNN trained to predict four hydrogen adsorption geometry labels, namely two H\({}_{2}\) bond lengths and two H\({}_{2}\) adsorption distances. To achieve a better performance, we inform our model of the fact that when H\({}_{2}\) is far away from TM, it tends to maintain its equilibrium bond length by adding a bond-order term in loss function, making it a physics-informed (PI) version. Compared to normal ALIGNN, our PI-ALIGNN's mean MAE of 5 folds is reduced by 0.00013 A, which is about 4% of total MAE (Table S3-4). The mean MAE of 0.00295 A in 5-fold cross-validation is in the same order of magnitude as the error range in DFT calculations for geometry (around 0.002 A), making our model reliable in predicting \(d_{\rm H}\). As shown in Figure 3(c-d), the adsorption geometries of the H\({}_{2}\) are well predicted with a MAE of 0.00312 A for \(d_{\rm H}\) and 0.166 A for \(d_{\rm MH}\) in test set. More details of ALIGNN model's implementation are given in Supporting Information.

Then, this two-step model is applied to the remaining database. The first H\({}_{2}\) quasi-molecule classifier ALIGNN picked out 9,883 structures with two quasi-molecule labels from 12,647 pieces of data. The H\({}_{2}\) adsorption geometry predictor PI-ALIGNN predicted four hydrogen adsorption geometry labels of these 9,883 MXenes.

**Analysis to MXene's hydrogen storage mechanism**

To understand the mechanism of hydrogen adsorption on MXenes both qualitatively and quantitively, we choose the five mono-component FGs, i.e. H, O, F, Cl and Br, which have relatively simple FG-H\({}_{2}\) interaction, to analyze the dependence of H\({}_{2}\) bond length on adsorption distance (Figure 4). The rest data containing binary and ternary component FGs are presented in Figure S4. Note that we don't discriminate up and bottom H\({}_{2}\) in the following discussions.

In the existing research literatures, it has been well accepted that the \(d_{\rm H}\) increases as \(d_{\rm MH}\) reduces due to the electron donation and back-donation between the transition metal's \(d\) orbital and H\({}_{2}\)'s \(\sigma^{*}\) orbital. Such picture gives rise to a natural assumption that the total valence bond order of hydrogen atom is unity (since each H has one electron) and can be composed of individual valence bond orders, namely, we have the following equation[10],

\[p_{\rm HH}+p_{\rm MH}=1 \tag{1}\]

where \(p_{\rm HH}\) and \(p_{\rm MH}\), the bond order between H-H and TM-H, is given by[22],

\[p_{\rm KH}=exp\left(-\frac{r_{\rm MH}-r_{\rm MH}^{0}}{p_{\rm KH}}\right) \tag{2}\]

where X is H/M,\(r_{\rm XH}\) is the distance between X and H atom, \(r_{\rm MH}^{0}\) represents the equilibrium X-H distance of free donor XH (when \(p_{\rm KH}=1\)) and \(b_{\rm KH}\) is a decay parameter. A quantitative relation reads[10],

\[r_{\rm HH}=r_{\rm HH}^{0}-b_{\rm HH}\,ln\left[1-exp\left(-\frac{r_{\rm MH}-r_{ \rm MH}^{0}}{b_{\rm MH}}\right)\right] \tag{3}\]

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \begin{tabular}{c} element \\ (ligancy 6) \\ \end{tabular} & Sc & Ti & V & Cr & Y & Zr & Nb & Mo \\ \hline 
\begin{tabular}{c} valence \\ number \(\nu\) \\ \end{tabular} & 3 & 4 & 5 & 6 & 3 & 4 & 5 & 6 \\ \hline \(\nu_{L}^{2}\) & 6.25 & 5.48 & 3.62 & 1.86 & 6.25 & 5.48 & 3.62 & 1.86 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Different values of \(\nu_{L}^{\nu}\) for TMs in MXenes

Figure 4: Predicted along with calculated hydrogen adsorption geometries of the five mono-component FGs separated by different TMs with (a) Sc, (b) Ti, (c) V, (d) Cr, (e) Y, (f) Zr, (g) Nb and (h) Mo. Dashed lines are from Eq. (3) and solid lines from Eq. (6), our formula for H\({}_{2}\) adsorption geometry.
where we set \(r_{\text{MH}}^{0}=1.6\) A and \(b_{\text{H}}=b_{\text{MH}}=0.404\) A according to the literature [10] and equilibrium H\({}_{2}\) bond length \(r_{\text{HH}}^{0}=0.75\) A in line with our computational result. Given that

\[r_{\text{HH}}=d_{\text{H}},\ r_{\text{MH}}=\sqrt{d_{\text{MH}}^{2}+\left(\frac{ d_{\text{H}}}{2}\right)^{2}}, \tag{4}\]

The dashed lines in Figure 4 are from Eq. (3), suggesting that the Eq. (3) is not suitable for the functionalized MXenes. It was originally derived for open TM sites in complexes and TM hydrides, which overestimates \(d_{\text{H}}\) in our systems with normally coordinated metal sites (NCMSs). TM sites in MXenes have a ligancy 6 which lowers their ability to activate dihydrogen bonds. A revised equation is needed for MXenes and other materials by taking ligancy and valence into consideration.

For this purpose, we adapt a new formula for the bond order \(p_{\text{MH}}\) between TM and H atom, deduced from Pauling's RVB theory through a statistic treatment of valence electrons [24, 23],

\[p_{\text{MH}}^{*}=exp\left(-\frac{r_{\text{MH}}-r_{\text{MH}}^{0}}{b_{\text{ MH}}}\right)\times\frac{1}{1+A\left(v_{\text{L}}^{0}-1\right)} \tag{5}\]

where A is a coupling coefficient of resonating between different structures and \(v_{\text{L}}^{0}\) is the number of unsynchronized resonance structures per atom, determined with valence \(v\) and ligancy \(L\), given by Eq. (54) (see Supporting Information for detailed deductions). Eq. (5) gives a basic picture of the resonating process of valence bonds. Strength of resonance is represented with coupling constant A and possible resonating structures \(v_{\text{L}}^{0}\). The valence electrons are delocalized due to resonance so that the bond strength, or bond order, decreases under the same bond distance \(r_{\text{MH}}\) The bond order of dihydrogen bond is still described with Eq. (2) since H\({}_{2}\) has no ligand.

The values of \(v_{\text{L}}^{0}\) for TMs in this work is shown in Table 1, and the revised formula for H\({}_{2}\) adsorption geometry on NCMS reads

\[r_{\text{HH}}=r_{\text{HH}}^{0}-b_{\text{HH}}\,n\left[1-exp\left(-\frac{r_{ \text{MH}}-r_{\text{HH}}^{0}}{b_{\text{MH}}}\right)\times\frac{1}{1+A\left(v _{\text{L}}^{0}-1\right)}\right]. \tag{6}\]

We choose A = 0.3 to achieve a good agreement with our data (solid lines in Figure 4). Through \(v_{\text{L}}^{0}\), the activation of H\({}_{2}\) bond on MXenes is properly described by Eq. (6) as the solid lines plotted in Fig.3. What's more, TM with more valence electrons tends to have longer H-H bond length at the same adsorption distance, implying a better activation of H-H bond. From the perspective of Pauling's RVB theory, the increase of \(v_{\text{L}}^{0}\) for the unsynchronized resonance structures per atom results from the lack of valence electron compared with atom's ligancy. The resonance between different electron configurations delocalizes the electrons, thus TM's ability to activate H\({}_{2}\) bond is weakened. Therefore, TM with more valence electrons tends to perform as better dihydrogen bond activator in MXenes, whose TM sites all have same ligancy of 6.

Eq. (6) gives us a quantitative correlation between dihydrogen bond length and H\({}_{2}\)'s adsorption distance from NCMSs on MXenes with mono-component FGs, which can be readily extended to other crystal material surfaces with similar local structures.

## Conclusion

In summary, to better understand the recent experiment on hydrogen storage in multilayer MXenes, we focus on the activation of H\({}_{2}\) on the functionalized MXenes, which can be measured by the H-H bond length changing with the adsorption distance. Based on high-throughput simulation, machine learning and Pauling's RVB theory, we derived a general formula for describing the relationship between activated H\({}_{2}\) bond length and adsorption distance, which could be extended to other 2D materials and materials surfaces to get insight on H\({}_{2}\)'s activation process.

## Associated Content

The implementation and detailed analysis of high-throughput simulation, ALIGNN model and related deductions in Pauling's RVB theory are available in the Supporting Information. This material is available free of charge via the Internet at [http://pubs.acs.org](http://pubs.acs.org).

## Author Information

* School of Materials Science and Engineering, Peking University, Beijing 100871, China; Center for Applied Physics and Technology, Peking University, Beijing 100871, China; Email: sunqiang@pku.edu.cn

## Acknowledgment

This work was partially supported by grants from the National Key Research and Development Program of China (2021YFB4000601), and from the China Scholarship Council (CSC). The calculations were supported by the High-performance Computing Platform of Peking University.
