Cation intercalation into metal oxide hosts is critical to many technologies, including energy storage [1], optoelectronics [2], and catalysis [3]. Numerous experimental and computational techniques have been used to understand this phenomenon, including analytical electrochemistry [4],

x-ray techniques [5], infrared measurements [6], and density functional theory (DFT) calculations [7]. These methods are regularly applied in concert [8] and provide multifaceted views of how a material engages in intercalation. However, these approaches can be costly and/or time-consuming, inhibiting the development of application-specific intercalation materials. This limitation is not restricted to cation intercalation and is generally relevant to processes involving the formation or breakage of chemical bonds.

Structure-based descriptors [9] have been employed to accelerate the characterization of bonding interactions in materials. These descriptors have been shown to correlate with difficult-to-obtain performance indicators, such as reaction rates and adsorption energies [10]. The simplest example is the local coordination number (LCN), which is an integer giving the number of bonds that connect a specific atom to its neighbors [11]. Notable extensions of the LCN are the generalized coordination number (for sites on metal nanoparticles) [12] and the adjusted coordination number (for sites on metal oxides) [13]. A further adaptation of the LCN is the orbitalwise coordination number [14], which incorporates additional information about electronic structure. Although each are powerful descriptors of chemical interactions, calculating these LCNs requires material-specific parametrization or the use of costly quantum chemical methods, limiting their ability to drive exploratory materials design.

In pursuit of increased generality and rapid prediction of bonding thermodynamics, more modern descriptors of structure venture beyond scalar LCNs to higher dimensional constructs that encode a given material as a graph [15]. While LCNs describe the local environment in which a chemical reaction occurs, graph descriptors detail the full global structure of a material. The state of the art in machine learning applied to describing chemical bonding [16,17] uses graph representations to accurately predict energetic behavior, suggesting that global connectivitynot just local coordination-could be used to describe and predict the thermodynamics of chemical bonding. It is therefore plausible that multiple scales of connectivity contribute to observable chemical behavior.

This Letter demonstrates how global and local connectivities can be used together to quantitatively predict the equilibrium energetics of hydrogen intercalation in metal oxides. We have previously investigated hydrogen intercalation in tungsten oxide using acid-base properties as descriptors [18]; however, these properties must be determined via computationally demanding quantum chemical calculations. This work expands the set of materials under investigation and uses inexpensive, structure-based descriptors to determine the energetics of hydrogen intercalation in oxides. Here, we introduce a mathematical definition of global connectivity, Ω. Importantly, Ω is a single number (not a full graph or matrix) that characterizes the overall connectivity of a structure. It is therefore as interpretable as a LCN and can be used to differentiate between materials based on their distinct geometries. Equations ( 1)-( 4) summarize the definition of Ω for a general atomic structure, and Fig. 1 illustrates the calculation of Ω for two contrasting examples of fictitious metal oxides.

Equation (1) gives the symmetric normalized Laplacian matrix [19] L, which catalogs a structure's atom-to-atom connections. L is an NxN matrix, with N being the total number of atoms contained in the molecule, nanoparticle, unit cell, etc. For the oxides in this work, N includes all metal, oxygen, and hydrogen atoms present in the unit cell. Therefore, N increases with higher levels of hydrogen insertion. In Eq. (1), u and v are atom indices and d i is the degree of atom i. Here, the degree is the LCN of an atom, i.e., the number of atoms which are directly bonded to it. Further discussion on L is included in Sec. S1 of the Supplemental Material [20]. Notably, L can be constructed without the use of bond distances and angles. It can be determined for an arbitrary system without any geometry optimization or specific knowledge of spatial arrangement. This is demonstrated in Fig. 1, which shows a more connected oxide in the top row and a less connected oxide in the bottom row. More connected systems are those that are more coordinatively saturated and densely bonded. The center column of Fig. 1 shows portions of the unique L representations of each oxide, computed according to Eq. ( 1).

The spectrum of eigenvalues for a given system, λ, is determined from L with Eq. (2). In this eigenvalue decomposition, λ is a vector containing N entries and X is the corresponding set of N eigenfunctions. λ is invariant to the ordering of rows and columns in L, i.e., the sequence of atom labeling does not impact the resulting spectrum. When ordering from smallest to largest, the first entry, λ 0 , is always zero. All remaining entries are strictly positive and lie between zero and 2 because of the normalization of L [19].

Ω is calculated from three quantities given by λ: the algebraic connectivity λ 1 [36], the spectral gap λ ϵ [37], and the Laplacian-energy-like invariant, LEL [38]. The invariance of λ asserts that these three values are also insensitive to atom numbering. As a result, they are known as graph invariants. Each holds unique information regarding the structure of the atomic system described by a particular L matrix.

λ 1 is the second entry of λ when the eigenvalues are ordered from smallest to largest. It has been shown to quantify the connectivity of networks represented by L, where smaller values are less connected and larger values are more connected [36]. For example, the more connected oxide in Fig. 1 has a λ 1 of 0.134, while the less connected oxide has a λ 1 of 0. A weakness of λ 1 is that it cannot distinguish different systems if they individually contain nonbonded sections: they will each have λ 1 values of 0, even if their geometries are dramatically dissimilar. Including λ ϵ in the definition of Ω helps alleviate this limitation. λ ϵ is the smallest nonzero eigenvalue of L [37]. It quantifies how easy it is to traverse from one point in a graph to another. Smaller values indicate that a structure contains multiple isolated parts which are connected by a small number of bonds. Because λ ϵ is defined to be nonzero, it can uniquely characterize layered systems that cannot be distinguished via λ 1 .

Finally, the LEL invariant is used to scale λ 1 and λ ϵ to make fair comparisons between systems containing different numbers of atoms. It is computed as the sum of the square roots of all entries in λ, according to Eq. (3). The LEL invariant is used as a normalization factor to account for effects on global connectivities computed for systems of different sizes.

Equation ( 4) defines the global connectivity, Ω, of a chemical system using each of the three graph invariants FIG. 1. Scheme for calculating Ω for two example oxide systems. The top system is more connected, and the bottom system is less connected; atoms in the top system are, on average, bonded to more neighboring atoms than those in the bottom system. Portions of the L matrices for each system are shown in the center column. Each entry is computed according to Eq. ( 1). Performing eigenvalue decompositions on L c and L l yields the full eigenvalue spectra, λ c and λ l . Ω c and Ω l are computed from these spectra according to Eq. ( 4). More connected systems have larger Ω values. discussed above. Ω ranges from 0 to 1 and can be interpreted as the fraction of λ associated with a system's connectivity. This is illustrated in Fig. 1, with the more connected oxide exhibiting a larger Ω than the less connected oxide. With this definition of Ω, we introduce the concept of global connectivity to explore the thermodynamics of hydrogen intercalation in metal oxides. First, we generated a dataset of hydrogen binding energies at various lattice oxygen sites in five bulk metal oxides. These binding energies were computed using DFT. Details are given in Sec. S2 of the Supplemental Material [20]. The oxides we considered include V 2 O 5 , MoO 3 , WO 3 , VO 2 , and TiO 2 , as described in Sec. S3 of the Supplemental Material [20]. Figure 2 shows how the DFT binding energies of hydrogen atoms in each oxide correlate with the inverse of Ω. As can be seen, oxides with smaller Ω values (less connected systems) exhibit more exothermic H-binding energies. These systems contain a large number of undercoordinated sites which are amenable to bonding with a guest hydrogen atom. Systems with higher Ω values have sites that are more coordinatively saturated, resulting in more endothermic (or less exothermic) binding.

Although Fig. 2 suggests a relationship between hydrogen binding energies and Ω, in some cases, single values of Ω returned multiple values of hydrogen binding energies. This is a consequence of the global nature of Ω: a given oxide has a unique Ω value, so that independent site types within its structure cannot be differentiated. These multiple site types also present many possible intercalation configurations [18]. The colored symbols in Fig. 2 show the energetics of the most stable binding configurations for each investigated oxide stoichiometry, and the black symbols show the unstable configurations. A line of best fit was drawn considering the most stable intercalated systems, since these are the ones with the highest probability of contributing to observed intercalation behavior.

Overall, Ω captured the trend of increasing H binding strength across oxides (TiO 2 < VO 2 < WO 3 < MoO 3 < V 2 O 5 ), illustrating its ability to describe meaningful structural traits of each metal oxide. However, using Ω alone as a descriptor for hydrogen binding energy gives an uncertainty > 0.3 eV, which is too large to quantitatively predict the equilibrium energetics of hydrogen intercalation. The invariability of Ω with atomic configuration also makes it impossible to obtain precise site-specific hydrogen binding energies from Ω alone.

Hydrogen (proton-coupled electron) intercalation into metal oxides can be viewed as two events: (1) binding of a proton to a lattice oxygen and (2) the addition of an electron in the diffuse oxide conduction bands. The global nature of Ω could capture information relevant to event 2. However, proton affinity (relevant to event 1) is dependent on local site coordination [39]. Therefore, we hypothesized that augmenting Ω with a local coordination descriptor could provide superior predictions of hydrogen insertion behavior in metal oxides. To construct an intercalation model capable of explaining site-specific energetics, we included local information through a scaled coordination number, CN s , defined by Eq. ( 5):

The numerator is the number of bonds within two nearest neighbors of the site of interest, demonstrated by Fig. S2 of the Supplemental Material [20]. The denominator is the total number of bonds in the unit cell (according to the relevant oxide crystal structure). Hence, CN s is the fraction of bonds in the unit cell which are associated with a given site. This description of local coordination differentiates each site type and assists Ω in describing hydrogen intercalation energies. Taken together, Ω and CN s quantify structure at global and local scales.

Figure 3 summarizes the performance of a kernel ridge regression (KRR) model that was used to calculate hydrogen intercalation energies using only Ω and CN s . This model was trained on DFT-calculated hydrogen binding energies in the investigated metal oxides. Full details on model construction are included in Sec. S5 of the Supplemental Material [20]. Intercalation chemistry is typically analyzed through the use of convex energy hulls, such as those shown in Fig. 3(a). The x axis shows the amount of hydrogen contained in the oxide lattice, and the y axis is the formation energy of a given intercalated oxide. As can be seen, the KRR model reproduced the DFT convex energy hulls with a root mean square error (RMSE) of 0.129 eV and a mean absolute error (MAE) of 0.089 eV FIG. 2. Correlation between hydrogen binding energies in extended metal oxides with the inverse of the global connectivity, Ω. Less connected oxides containing a larger number of undercoordinated sites bind hydrogen atoms more strongly, and more connected systems with fully saturated atoms bind hydrogen more weakly. The larger colored symbols represent stable configurations at each stoichiometry, and the black symbols correspond to unstable configurations. A line of best fit is drawn through the set of stable configurations.

(in terms of H-binding energy on a site-specific basis). Accordingly, Fig. 3(a) suggests that Ω and CN s could be useful descriptors of hydrogen binding in these metal oxides.

While the convex hulls identify stable phases of intercalated oxides, a more practical thermodynamic quantity is the hydrogen intercalation potential, U H . These potentials can be calculated from the slopes of the convex energy hulls (see Supplemental Material [20] Sec. S2). Figure 3(b) compares experimental U H values to those calculated from the KRR hulls. Each horizontal grouping of data points is a collection of experimentally measured U H values for a single H-intercalation equilibrium potential. For example, the red circles are the TiO 2 =H 0.125 TiO 2 redox couple. The vertical lines on each data point are the error associated with the model prediction, taken as the MAE of the KRR model (AE0.089 eV). The filled markers are data extracted from published literature reports [18,[40][41][42][43][44][45][46][47][48], and the open markers are from measurements we performed. Methods used to obtain the experimental data included in Fig. 3(b) are described in Secs. S6 and S7 of the Supplemental Material [20].

In general, the KRR model reproduced experimental potentials with an accuracy that would be expected from DFT calculated potentials, giving a RMSE of 0.116 eV. Furthermore, for oxides with data on more than one hydrogen intercalation equilibrium, the structure-based model correctly captured the phenomena of higher stoichiometries having more negative intercalation potentials. The performance of Ω and CN s versus other metrics of connectivity is provided in Sec. S8 of the Supplemental Material [20]. We reiterate that determining Ω and CN s did not require any quantum chemical calculations or other expensive computations other than generating data for developing the model. These features therefore provide a rapid route for determining hydrogen intercalation potentials in metal oxides that is based entirely on the structure of the parent oxides. In principle, the same statistical models could be built using experimental data from one set of metal oxides to predict insertion potentials of another set.

To interrogate the physical relevance of Ω to intercalation, we applied a convolutional neural network (CNN) to the same regression task addressed with the KRR model. The major difference in this approach was that the inputs were taken directly as the L matrices. During training, this allowed the CNN to autoembed representations of L which were significant to hydrogen intercalation. If the CNN embeddings correlate with Ω, that would suggest Ω is a meaningful metric of connectivity. The results of this analysis are shown in Fig. 4, and full implementation details of the CNN are described in Sec. S9 of the Supplemental Material [20].

The axes in Fig. 4 are the first and second principal components of the CNN's final hidden layer. Each data point is an intercalated oxide, and the color represents each oxide's global connectivity. There is a clear gradient in color moving from left to right in Fig. 4, suggesting that Ω is strongly correlated with the feature space that was generated by the network. The hidden layer node values, and therefore their principal components, are interpretable as features that have been self-generated by the network during training. Hence, without explicitly considering Eqs. ( 1)-( 4), the CNN arrived at a hidden layer feature space that closely resembled Ω. This interpretation of our model gives further evidence of system-scale, global connectivity as a physically meaningful descriptor of hydrogen intercalation thermodynamics in metal oxides.

In conclusion, this Letter introduced global connectivity Ω as a reactivity descriptor in redox-active metal oxides. We first derived an expression for Ω and illustrated its roots in spectral graph theory. Using DFT to calculate hydrogen intercalation energies in a set of five oxides, we showed that Ω captured the general thermodynamic trends of intercalation. Then, we combined Ω with a scaled local coordination, CN s , to build a regression model that depended only on structural information that is readily obtained from experiment. This structure-based model leveraged both global and local connectivities to accurately capture the trends and values of hydrogen intercalation potentials across all investigated oxides when compared against electrochemical experiments performed in the lab. To further support Ω as a meaningful descriptor of oxide redox reactivity, we applied a convolutional neural network to convert Laplacian matrices directly into hydrogen intercalation energies. We observed a strong correlation between the principal component space of the final hidden layer of the network and Ω, suggesting that the CNN embedded the linear algebra that links L and Ω. Future work can further develop the concept of global connectivity by incorporating more node-and edge-level information into L. These information-rich Laplacians could then be compared to more sophisticated network models, such as crystal graph convolutional neural networks [15]. Overall, this Letter introduces Ω as a useful metric for global connectivity and demonstrates how simplified representations of global and local structure can be used together to predict hydrogen binding in metal oxides.