Artificial intelligence pathway search to resolve catalytic glycerol hydrogenolysis selectivity

The complex interaction between molecules and catalyst surfaces leads to great difficulties in understanding and predicting the activity and selectivity in heterogeneous catalysis. Here we develop an end-to-end artificial intelligence framework for the activity prediction of heterogeneous catalytic systems (AI-Cat method), which takes simple inputs from names of molecules and metal catalysts and outputs the reaction energy profile from the input molecule to low energy pathway products. The AI-Cat method combines two neural network models, one for predicting reaction patterns and the other for providing the reaction barrier and energy, with a Monte Carlo tree search to resolve the low energy pathways in a reaction network. We then apply AI-Cat to resolve the reaction network of glycerol hydrogenolysis on Cu surfaces, which is a typical selective C–O bond activation system and of key significance for biomass-derived polyol utilization. We show that glycerol hydrogenolysis features a huge reaction network of relevant candidates, containing 420 reaction intermediates and 2467 elementary reactions. Among them, the surface-mediated enol–keto tautomeric resonance is a key step to facilitate the primary C–OH bond breaking and thus selects 1,2-propanediol as the major product on Cu catalysts. 1,3-Propanediol can only be produced under strong acidic conditions and high surface H coverage by following a hydrogenation–dehydration pathway. AI-Cat further discovers six low-energy reaction patterns for C–O bond activation on metals that is of general significance to polyol catalysis. Our results demonstrate that the reaction prediction for complex heterogeneous catalysis is now feasible with AI-based atomic simulation and a Monte Carlo tree search.


Introduction
Glycerol is an important biomass-derived platform chemical [1][2][3][4] and its selective conversion to other chemicals, 5-10 particularly, 1,2propanediol (1,2-PDO) 11-16 and 1,3-propanediol (1,3-PDO) [17][18][19][20] are of great signicance in industry. Among the many catalysts tested, Cu-based catalysts are regarded to be the best to selectively convert glycerol to 1,2-PDO with low/zero selectivity to 1,3-PDO and other deep hydrogenolysis products (e.g. propanol). Considering that 1,3-PDO is a more valuable chemical than 1,2-PDO, 7 there has been considerable interest to quantify the reaction kinetics and to improve the 1,3-PDO selectivity. However, due to the complex nature of the polyol interaction with catalyst surfaces, the reaction mechanism remains hotly debated: 3,7 the high selectivity to 1,2-PDO for glycerol hydrogenolysis on Cu-based catalysts has been puzzling for decades. The complexity originates mainly from the huge reaction network with many likely reaction intermediates due to the presence of all 13 chemical bonds, i.e. ve C-H, three O-H, three C-O and two C-C, in glycerol that are apparently possible to break on catalysts. 3 Fig . 1 shows the proposed reaction mechanisms for glycerol hydrogenolysis in the literature dating back to the 1980s: six different reaction routes to 1,2-PDO and 1,3-PDO 5,12,18,21-23 were proposed, each in a different color. Since 1,2-PDO is observed to be the major product on Cu-based catalysts (see the ESI, Table S1 † for collected experimental data), the pathways to 1,2-PDO assume the preference of breaking one of the two primary C-O bonds but with varied patterns, including 1,2 b-elimination dehydration (blue arrows, 12,22 Fig. 1), direct C-O bond breaking (red, 23 green 23 ), and retro-Michael addition reaction (black 5 and purple 21 ). These pathways are, however, in contrast with the general chemistry knowledge that the secondary C is chemically more active than the primary C, and more importantly, cannot explain satisfactorily the high selectivity to 1,2-PDO but low/zero selectivity to 1,3-PDO (brown 18 ) and the deep hydrogenolysis products.
The glycerol hydrogenolysis reaction thus illustrates well the dilemma of selectivity control in heterogeneous catalysis where the current knowledge on elementary reaction steps, being difficult to glean from experiment, is much short of one's expectation towards the rational design of the reaction. The participation of catalysts (e.g. metal surfaces) in reactions not only greatly increases the likely intermediates (e.g. adsorbed atoms, radicals), but also creates catalyst-specic preference for specic reaction patterns. The theoretical methods available to date generally fail to provide fast and accurate guidance to catalyst and reaction design. In particular, the valence bond theory well recognized in chemistry oen cannot predict even qualitatively the activity and selectivity in heterogeneous catalysis.
Instead, quantum mechanics based atomistic simulations are now popularly resorted to for elucidating the reaction kinetics of heterogeneous catalysis, which are however limited largely to small reaction systems (e.g. less than 100 atoms) on well-dened catalyst surfaces due to the high cost and poor scaling in computation. In fact, theoretical attempts were made to expedite the surface reaction exploration as early as the 1990s. These methods involve rst the generation of the likely reaction patterns of molecules according to the octet rule or the bond matrix, such as Netgen (Broadbelt et al.), 24,25 RMG (Green group), 26,27 RING (Rangarajan et al.) 28 and pReSt (Valchos group). 29 The kinetics data of reactions between enumerated reaction pairs, i.e. the reaction energy (energy difference between the product and the reactant) and reaction barrier (energy difference between the transition state and the reactant), while being difficult to obtain, can then be estimated by empirical/approximate approaches, such as the Brønsted-Evans-Polanyi (BEP) relationship, 30,31 scaling relationship, [32][33][34] and UBI-QET method. 35,36 In these methods, the reaction energy may be fetched either from experimental data, or estimated by the group additivity approach, the gas phase bond energies, and the atomic binding energy, or directly calculated by on-the-y density functional theory (DFT) calculations, whilst the reaction barrier is estimated from the reaction energy by the BEP relationship or can be explicitly determined by locating the transition state (TS). [37][38][39][40][41][42][43] While the scheme appears logically sound, there are still many methodological difficulties associated with the accurate reaction prediction. For example, the BEP relationship is a crude approach, i.e. the error bar of reaction barrier estimation is generally above 0.2 eV; the TS location requires the pre-knowledge of reaction coordinates and has a high computational cost even with DFT.
Recent years have seen great success in developing articial intelligence (AI) methods for organic retrosynthesis, 44-56 which relies on two key merits that appear to be unique to organic reactions: rst, a large reaction database from organic experiments is available on the market (e.g. Reaxys, and SciFinder databases); second, the general validity of the octet rule in describing elementary organic transformation. It is thus possible to relate an organic reaction to the reactant and product using numerical descriptors, such as Simplied Molecular Input Line Entry Specication 57 (SMILES) and Extended-Connectivity ngerprints 58 (ECFP) where the bond order matrix can be used to detect the structural change. On the other hand, the situation in heterogeneous catalysis is rather different. The reaction steps occurring on surfaces are numerous and generally not separable until the nal products are obtained. As exemplied in Fig. 1, acetol is the only one among all proposed intermediates that can be conrmed via experiment. 12,23 Furthermore, an accurate account of heterogeneous reactions is virtually not possible due to the structural complexity of solid catalysts (various surface structures and compositions). It is no wonder that the reaction database for heterogeneous catalysis is not established from experiment, not even mentioning the AI-based reaction prediction based on the heterogeneous reaction database.
Herein we propose a general AI framework for predicting heterogeneous catalysis activity, namely the AI-Cat method. Our method provides a one-go solution, from the building of a heterogeneous catalytic reaction database, to the end-to-end activity prediction using machine learning. The reaction database is collected by using our recently developed global neural network (G-NN) potential 59-61 that, in combination with stochastic surface walking (SSW) global optimization, 62 can fast explore reactions on surfaces. The end-to-end activity predictor combines two neural network (NN) models, one for predicting possible reaction patterns and the other for providing the reaction barrier and energy, with the Monte Carlo (MC) tree search to resolve the low energy pathways in a reaction network. We illustrate our method in resolving the complex reaction network of glycerol hydrogenolysis on Cu and rationalizing the high selectivity of glycerol to 1,2-PDO.

Theoretical methods
Owing to the lack of experimental kinetics data on elementary reactions in heterogeneous catalysis, any successful AI model for predicting heterogeneous reaction activity must be able to build a reaction database by its own, from which the reaction pattern and kinetics information can be abstracted. The subsequent machine learning of these reaction data would then require sensitive reaction descriptors that can distinguish molecules/fragments on surfaces. Before proceeding to present our AI-Cat method, we rst introduce the molecule descriptor for surface reactions developed in this work.
It is well established that an organic molecule can be uniquely identied using the atom-centered extended-connectivity ngerprints, i.e. ECFP-x where x represents the largest possible fragment that has a width of x bonds (diameter in bonds) with respect to the center atom. ECFP can abstract the topological information of the molecule from the bond matrix, which can be readily established for heavy atoms (non-H atoms) if the octet-rule is satised. The ECFP is thus not quite compatible with heterogeneous catalytic systems for two reasons: (i) unsaturated reaction intermediates (adsorbates on surfaces) are common; and (ii) H atoms, either as adsorbates or attached to molecular fragments, are of key signicance to distinguish different molecular states.
For describing surface reactions uniquely, we modied ECFP descriptors in three places: (i) the bond matrix generation allows the adsorbed atoms to not obey the octet rules and takes into account the periodicity, if necessary; (ii) the information of the catalyst, e.g. the type of metal, is embodied into the initial identier for adsorbed atoms; (iii) H atoms are treated equally as the heavy atoms. The thus-modied ECFP-x is named surface sensitive ECFP-x (s-ECFP-x), which can produce the unique identier for adsorbed molecules/fragments. For example, for an adsorbed HCOO that has an O adsorbed on the surface (HCOO*, * indicates the adsorbed atom), s-ECFP-x will consider the fragment to be formate HC(]O)O* by maximally letting non-adsorbed atoms (HCO) to obey the octet rule. By this way, the adsorbed carboxyl (C*(]O)OH) and formate (HC(]O)O*) can be readily distinguished.
An example of generating the s-ECFP-4 digital code for a coadsorbed formate fragment and H atom is given in Fig. 2 (also see ESI Section 2 † for details). The new features in s-ECFP as compared to the traditional ECFP are highlighted by the red explanatory text. An iterative procedure is applied to every atom to create a series of identiers by gradually increasing the radius around the central atom, as indicated by the red circles around 1 O. Every iteration takes into account the topological information related to the circled atoms and the bonds connecting them. Aer all atoms are treated, a set of s-ECFP-4 ngerprints are generated to form a vector as the unique descriptor of the structure. We would like to comment that the s-ECFP-4 descriptor can distinguish different adsorbed molecules, including single molecules and coadsorbed molecules (e.g. HCOO and H in Fig. 2), but does not record their relative physical distance. This treatment is largely valid for metal catalysis since the adsorbate diffusion is generally facile and thus not important for evaluating the catalytic activity.
We now turn to the AI-Cat method, which relies on the s-ECFP-4 descriptor to distinguish different reactions. There are two parts in the AI-Cat method, i.e. (i) reaction database construction; (ii) AI model for activity prediction. The AI model can be further decomposed into two NN prediction models and an MC tree search module. They are elaborated in the following section.

Algorithm for constructing the reaction database and the reaction database for catalytic reactions on Cu surfaces
With the advent of G-NN potential based global potential energy surface (PES) sampling as implemented in the SSW-NN method, it is now feasible to predict unknown chemical systems in an automated way. In the past few years our group has managed to generate a series of four-element G-NN potentials M-C-H-O (M ¼ Fe, Co, Ni, Cu, Ru, Rh, Pd Ag, Ir, Pt and Au) covering the most common catalytic metals and these G-NN potentials are accessible from the LASP 61 project website (https:// www.lasphub.com). These G-NN potentials (see ESI Section 2 † for details) allow for fast PES evaluation in the SSW reaction sampling (SSW-RS) 39,63 of common C1-C3 reactions on these metal surfaces. 64 This work focuses on the building of the reaction database for catalytic reactions on Cu surfaces.
In initiating the automatic reaction database generation, we require a minimum amount of information as input, including the metal species (e.g. Cu), the type of surface (e.g. (111), (100) surfaces), and the C : H : O composition (e.g. the type and the number of molecules, such as CO, H 2 , H 2 O, methanol, and glycerol). The program then utilizes the SSW-RS method to search for all the likely distinct reactions at the given C : H : O composition on the surface, which is an iterative procedure to explore all the likely reactions until no new low energy minima are identied.
The scheme of the SSW-RS method is shown in Fig. 3a with each small circle representing a minimum and the dotted large circle representing the range of SSW-RS explorations centering Fig. 2 Procedure to generate identifiers for molecules on surfaces using Surface-sensitive atom-centered Extended-Connectivity FingerPrint (s-ECFP). The iterative procedure for generating s-ECFP is illustrated using the adsorbed O end ( 1 O) of HCOO as the example. Each iteration creates an identifier for the substructure in a circle around the central atom. The identifier at iteration 1 only represents the information of the 1 O atom, its attached bonds and the surface in bonding. By enlarging the circles, the next iterations include more and more information associated with the substructure centering on the 1 O. Duplicate identifiers will be removed, for example, the substructure of 1 O at the iteration 3 is the same as the one of atom 2 C at the iteration 2. The explanatory text in red indicates the added features for s-ECFP compared to the ECFP. this minimum. By automatically switching the starting minima in different iterations as indicated by the cycle n (n ¼ 1, 2, ., see Fig. 3a), the global PES will be mapped out and all the likely reaction pairs, initial state (IS) and nal state (FS), together with the identied transition state (TS) can be obtained. The SSW-RS method has been described in detail previously, for example, to explore the reaction space of gas phase organic reactions. 65 In this work, the s-ECFP descriptor is utilized to judge whether a newly found reaction is distinct from others. Fig. 3b and c illustrate how a reaction database is established from the raw SSW-RS data. We use the following two strategies for purifying the raw data to maintain the key kinetics information. First, as there may well be many different pathways for the same elementary reaction (Fig. 3b), we only abstract the lowest energy minima for stationary points (IS and FS), and the lowest energy for the saddle points (TS) as long as the TS connects correctly to the IS and the FS as indicated by s-ECFP. This treatment ignores the kinetics for the transformation between molecular congurations on the surface (e.g. diffusion) since they are in general low-barrier fast reactions, particularly for C1-C3 small molecules. This stores the correct reaction energy and the lowest energy barrier between the reaction pair. Second, too rare reactions, e.g. those identied only once from SSW-RS, are not considered since kinetically relevant reactions should occur frequently.
Based on the elementary reaction pairs, the reaction patterns can then be extracted and also stored in the reaction database with a unique reaction pattern identier (indexed by integers from 1 to n). Following the shell-based scheme to describe a reaction, proposed by Christ et al. 66 and Saller and coworkers, 67 we dene the reaction patterns of surface reactions by the bond matrix elements for the reaction center atoms and their rst neighbors, and the reaction center atoms' identiers (s-ECFP) (also see ESI Section 3 † for details). As illustrated in In addition, we also record the s-ECFP identier of reaction center atoms of reactants for the purpose of reaction center recognition. The reaction pattern identier will be the output for the AI model for reaction prediction, and utilized for predicting the kinetics data of a reaction, as presented below.

AI-Cat model for activity prediction
The AI model for activity prediction utilizes the MC tree search algorithm to identify the low energy pathways connecting the Each small circle represents a distinguishable state that has many connections (pathways) to other small circles. Each SSW reaction sampling is represented by the dotted large circle that collects the pathways linking to the centering state. By automatically selecting the state for sampling, the reaction network for a complex reaction system can be established (also see text). (B). Illustration of multiple pathways connecting to the same reaction pair, where only the lowest energy IS, TS and FS, i.e. the reaction pair information, are recorded in the reaction database to provide the reaction barrier and the reaction energy. (C). Illustration of extracting the reaction pattern associated with the reaction pair using the bond matrix where the reaction center atom is uniquely identified by s-ECFP. Finally, the reaction database includes both the reaction pattern and the reaction pair information.
input reactant to the target product. During the reaction network search, the MC tree search algorithm needs to acquire two NNs persistently, namely a reaction pattern (R-Pat) unit to enumerate the most likely reaction patterns of a given reactant, and a kinetics information (K-Info) unit to provide the essential kinetics data.
2.2.1 MC tree search. The MC tree search is the central unit governing the reaction network exploration and realizing the activity prediction. The method was commonly utilized for the sequential decision-making problems, such as games, automated theorem proving and more recently retrosynthesis. 45,68 We here briey introduce the algorithm with emphasis on its connection with the R-Pat and K-Info units.
We consider a whole reaction network as a tree structure with each node representing an intermediate state (containing all species, including both the reacting species and spectator species) and the edges being the state-to-state connections governed by kinetics. Each node is associated with three values that can be updated during tree search, i.e. the action value Q, the visit count N and the reaction probability P. For a newly identied node, the initial value for Q and N is 0, and that for P is set as the probability from the father node to the current node computed from eqn (1), which is a weighted rate constant, where i denotes the current node and E a(j) is the reaction barrier from the father node to its jth child node in eqn (1). Terminal nodes need to be dened for stopping a search along one tree branch, which belong to either one of the three situations, (i) the target; (ii) the dead-end states with no further connections or connecting only to those with no further connections, and (iii) the states with too high overall barriers to access. The overall barrier of a multiple-step process corresponds to the energy difference calculated from the TOF-determining transition state (TDTS) and the TOF-determining intermediate (TDI) in the energetic span model 69 (the maximum energy height needed to overcome along a reaction pathway).
As shown in Fig. 4A, MC tree search has four phases that are executed iteratively until the lowest energy pathway is identied.
Phase 1 selection: starts from the root node, a selection procedure is performed iteratively to identify the most promising next node according to the priority score V i .
In eqn (2), the subscript i and fa indicate the current node and its father node. The selection traverses the tree by selecting the node with maximum V until reaching an unvisited node (N ¼ 0). The search will then move to phase 2 if the selected node is unvisited, or to phase 4 if the node needs to be expanded (N ¼ 1). Phase 2 rollout: starts from the selected, non-visited node, a rollout procedure is carried out by enquiring the R-Pat unit to obtain the target product. The rollout stops if a target node is met or the maximum depth D max of the search is reached. Every enquiry of the R-Pat unit gives a number of possible reaction patterns, C RP , leading to the new product nodes. Due to the fast increase of nodes, we limit C RP to 3, i.e. the top 3 likely reaction patterns, and D max to 10.
Phase 3 reward: Q and N are updated immediately aer phase 2 for the nodes traversed in the branch from the selected node to the root node. Q will be rewarded by adding 1 if the target product is reached, or be punished by deducting 1. The visit count N increases by 1 for the whole branch nodes. The current node thus changes to a visited node. Phase 4 expansion: starts from the selected, visited node (N ¼ 1), the expansion of the tree is performed by enquiring both the R-Pat and the K-Info units to generate the child nodes for the current node. Fig. 4B shows the owchart of the expansion phase. As shown, the reactant is rstly converted to an s-ECFP-4 code and fed into the R-Pat unit, which returns a probability distribution over all possible RPs (RP 1 to RP n ). Then, for the k most probable RPs, by applying each one to the reactant for generating the product (the child node) and thus the reaction R, the K-Info unit is utilized to obtain the reaction barrier and reaction energy associated with each child node. Finally, the child nodes are ranked by the reaction barrier info. To avoid the explosion of the tree size, we limit the maximum number of child nodes to 15, where the child node should not appear previously in the same branch and the overall barrier to the child node needs to be less than a preset maximum barrier (e.g. 2.5 eV or 0.3 eV above the current lowest energy pathway, the reference).
By iteratively executing these four phases, a tree of reaction networks is established and the pathways from the reactant to the target product are recorded and ranked according to the predicted overall reaction barrier. It is worth mentioning that the key feature of our AI-Cat model is the accessibility to the kinetics information of predicted elementary reactions. By setting the current lowest energy pathway as the reference, we can efficiently rule out the node with a too high overall barrier (e.g. 0.3 eV higher than the reference with the consideration of the likely error of K-Info unit prediction) in the phase 1 selection (by setting as terminal nodes). This leads to a highly efficient removal of irrelevant and redundant branches and thus the tree search converges quickly to the low energy pathways. It is worth mentioning that the complexity of MCTS depends on the input chemical composition: the molecules dened by the input and the available reaction patterns learned from the database.
2.2.2 R-Pat unit. The R-Pat unit takes the input from the reactant identier and yields the likely reaction pattern identi-ers as the output. The input is given by a vector with the value indicating the presence/absence of a particular s-ECFP-4 descriptor in the reactant. In this work, we utilize feed-forward NN to train the R-Pat unit by learning the reaction database, where the network has two hidden layers with 512 and 512 neurons of Relu nonlinearities, and the last layer is a Somax layer giving the probability distribution for all possible reaction patterns appearing in the reaction database. By applying these reaction patterns to the input structure, we can recover the associated reactions dened by the reactant and product.

K-Info unit.
The K-Info unit takes the input from the reaction information and outputs the reaction barrier and reaction energy. The input is a vector concatenated by three vectors, two from the s-ECFP-4 of the reactant and product, and one from the general reaction information that has nine dimensions, containing the highest, lowest and average reaction barriers and reaction energy, and the surface type. Similar to that in the R-Pat unit, a feed-forward NN is utilized to train the K-Info unit by learning the reaction database, where the network has two hidden layers with 512 and 256 neurons of Sigmoid nonlinearities, and the last layer is a linear layer that gives give the reaction barrier and energy.

G-NN potential and DFT calculations
The quaternary Cu-C-H-O G-NN potential used for SSW-RS reaction data collection was trained by using LASP soware, 61 which iteratively learns the global PES data generated by SSW global optimization (see ESI Section 2 † for the detailed methodology). The Cu-C-H-O G-NN potential has a double-net architecture implemented in LASP 3.0, 70 a ve-layer (410-80-80-80-1) feed-forward NN for all elements plus an auxiliary four-layer (258-80-80-1) feedforward NN for C, H and O elements, equivalent to 266 823 parameters in total. The potential was trained using a dataset of 90 726 distinct structures that are computed using plane wave DFT calculations 71,72 (DFT calculation setup can be found in the next paragraph). The nal RMS accuracy of the G-NN potential is 3.77 meV per atom for energy and 0.10 eVÅ À1 for atomic force (the G-NN PES accuracy is also addressed in ESI Fig. S5 †).
DFT calculations have been utilized to generate the dataset for G-NN training and also for validating the surface reactions reported in this work. All the total energy calculations were carried out using the VASP soware, 71,72 where the basis set is expanded in plane waves and the electron-ion interaction is represented by the projector augmented wave (PAW) pseudopotential. For the G-NN training dataset, the electron exchange and correlation effects are described by the GGA-PBE functional. For reaction pathway validations, the van der Waals correction using the Grimme D3 approach (DFT-D3 (ref. 73)) has been added to the GGA-PBE total energy. In all calculations, the kinetic energy cutoff utilized is 450 eV and the rst Brillouin zone k-point sampling utilizes the Monkhorst-Pack scheme with an automated mesh determined by 25 times the reciprocal lattice vectors. The nite difference method with the displacement of 0.01Å for an atom is utilized to calculate the vibrational modes and thus to determine the zero-point energy (ZPE).
To obtain the free energy prole of glycerol hydrogenolysis (see Section 3.1), the ab initio thermodynamics approach 74 is utilized to correct the Gibbs free energy under the reaction conditions (473 K). For all of the adsorbed states on the surface, the vibrational entropy is taken into account based on the harmonic frequency calculation using DFT-D3. The free energy of gas phase molecules (glycerol, acetol, 1,2-PDO, and 1,3-PDO) is computed under the typical experimental conditions (473 K, total pressure of 1 atm, glycerol : H 2 : H 2 O ¼ 1 : 140 : 12) used in ref. 14, according to the standard thermodynamic equations.

Results
In this work we have explored extensively the common organic reactions up to four C atoms (C1-C4) on three Cu surfaces, namely, close-packed Cu(111), open Cu(100) and stepped Cu(211), by using automated SSW-RS sampling with the quaternary Cu-C-H-O G-NN potential. Our simulation yields more than 200 000 reaction pairs and 24 130 different intermediates, and aer removing the duplicate ones, the nal reaction database contains 30 508 distinct reactions, 7776 reaction intermediates composed of 1188 distinct fragments/ molecules, and 5136 reaction patterns.
In our AI-Cat model, the R-Pat unit learns 7776 intermediates and 5136 reaction patterns associated with them, 90% of them being set as the training set and the remainder as the testing set. The probability of occurrence of a reaction pattern is determined by the weighted rate constant in eqn (1), which is used to supervise the training. Aer training the model, the top 3 accuracy for predicting the most probable reaction pattern is perfect (99%) for the training set and 57% for the testing set; and the top 15 accuracy for the testing set reaches 81%.
The K-Info unit learns the reaction dataset of 30 508 reactions on three different Cu surfaces. By splitting the dataset into 90% training set and 10% testing set, the NN nally achieves the accuracy of prediction with MAE of 0.061 and 0.054 eV for the reaction barrier and reaction energy on the training dataset, and MAE of 0.157 and 0.149 eV for the barrier and reaction energy on the testing dataset, respectively. For the purpose of the benchmarking of the AI-Cat model, we have rst analyzed the AI-Cat accuracy on the reaction dataset and conrmed the good performance of AI-Cat for different types of reactions, which is much better than the simple BEP model (see ESI Section 5 † for details). We then tested AI-Cat prediction ability for 21 different compositions by mixing 7 common small molecules, including H 2 , CO, CO 2 , H 2 O, CH 3 OH, CH 4 and CH 2 ]CH 2 , and AI-Cat gives the correct activity prediction for the known reactions (see ESI Section 6 † for details). We thus can move on to treat the concerned system with a complex reaction network, the glycerol hydrogenolysis on Cu.
By using the AI-Cat method, we are able to predict the likely hydrogenolysis pathways starting from glycerol and H 2 reactants on Cu(111). Aer 150 000 MCTS, we obtained 420 distinct intermediates and 2467 distinct elementary reactions, from which the lowest energy pathways are collected.

Reaction pathways and microkinetics of glycerol hydrolysis on Cu(111)
Our AI-Cat model predicts that 1,2-PDO is preferred over 1,3-PDO with a lower overall barrier among total eight products that have relatively low formation barriers (<1.6 eV). The key pathways are summarized in Fig. 5 and all pathway data from AI-Cat are detailed in the ESI. † Here we summarize the lowest barrier products, which are dihydroxyacetone, acetol, 1,2-PDO, prop-2ene-1,2-diol, 2-oxopropanal, glyceraldehyde, 3-hydroxypropanal, and 1,3-PDO, where the predicted energy barrier from AI-Cat differs by less than 0.3 eV. Among these products, dihydroxyacetone and glyceraldehyde, while being facile to form, are both endothermic intermediates and can further convert to 1,2-PDO and 1,3-PDO. Acetol and 3-hydroxypropanal are intermediates leading to 1,2-PDO and 1,3-PDO, respectively.
The remaining dehydration products, i.e. prop-2-ene-1,2-diol and 2-oxopropanal, are less exothermic to form than 1,2-PDO, but share the similar rate-limiting steps with 1,2-PDO. Therefore, from the AI-Cat energetics, only 1,2-PDO and 1,3-PDO are likely products for glycerol hydrogenolysis on Cu surfaces, where both molecules are thermodynamically highly favored products (exothermic by more than 0.7 eV). . The solid lines are the dehydrogenation (deH) pathways initiated by blue arrows in (A) and the dotted lines are the dehydration (deH 2 O) pathways initiated by hydrogenation (green arrow) and dehydration (red arrow) in (A). The energy zero refers to glycerol in the gas phase at 473 K and 0.0065 atm; and its dehydrogenated product, state 2 (see (A)), is the most stable species (À0.08 eV) before C-O bond breaking. The snapshots and free energy of key reaction transition states, TS1 to TS4, are also shown.
In order to determine the selectivity accurately and perform the microkinetics simulation, we have further utilized DFT calculations with van der Waals correction (DFT-D3 (ref. 73)) to rene the reaction prole of the four lowest energy pathways, where all minimum structures are fully relaxed and the TSs are re-searched (the detailed comparison between G-NN and DFT energetics is detailed in the ESI, Fig. S6 †). All Gibbs free energies reported hereaer are from DFT based energetics.
Both 1,2-PDO and 1,3-PDO have two low energy pathways as shown in Fig. 5, the dehydration pathway (Fig. 5a middle-to-top/ bottom vertical lines) and the dehydrogenation pathway (Fig. 5a le/right circles), as distinguished by their initial steps. In addition to these four low energy pathways, 1,2-PDO also has a retro-Michael addition pathway, bifurcating from the dehydrogenation pathway to 1,3-PDO (see ESI Fig. S7 †). Fig. 5b highlights the lowest energy reaction proles for glycerol hydrogenolysis to 1,2-PDO and to 1,3-PDO, which all belong to the dehydrogenation pathway (the circles in Fig. 5a), together with the two higher barrier pathways via initial dehydration.
By analyzing the low energy pathways from AI-Cat, we can glean key information on this catalytic system. (i) The lowest barrier elementary reactions are the dehydrogenation (blue arrow in Fig. 5a) and the hydrogenation (green arrow) reactions, which are followed by the C-O bond breaking (OH or H 2 O removal) reactions (red arrow). As a result, in between every C-O breaking, there are a series of dehydrogenation/hydrogenation fast equilibrium steps, e.g. glycerol with dihydroxyacetone, 1,2-PDO with acetol; (ii) the secondary C is always the rst reaction site for both the dehydrogenation and the direct C-O bond breaking, which is in line with the general knowledge in organic chemistry; (iii) the 1,2 b-elimination dehydration (Fig. 1) and the C-C bond breaking are kinetically prohibited on Cu surfaces. These are consistent with the general tendency to retain C-O and C-C bonds for common C1-C2 reactions on Cu catalysts as shown in ESI Table S3. † In the lowest pathway to 1,2-PDO, the secondary OH (linking with the secondary C) of glycerol is rst dehydrogenated on Cu(111) (state 2 Fig. 5a) followed by the secondary C-H bond breaking to form dihydroxyacetone, which is also the key step determining the selectivity. Next, dihydroxyacetone breaks one of its terminal OH groups via a Cu-mediated keto-enol tautomeric resonance (surface Cu replacing H in the enol form) and then undergoes a hydrogenation reaction to acetol, which is the rate-limiting step with 1.25 eV overall free energy barrier with respect to the adsorbed state 2, the most stable state along the pathway. Aer the C-O breaking, acetol readily converts to 1,2-PDO with three low barrier steps, involving hydrogenation and water generation. The further hydrogenolysis from acetol to 2propanol is however kinetically prohibited (free energy barrier 1.45 eV) because of the high stability of acetol (see ESI Fig. S8 † for the reaction prole of deep hydrogenolysis).
The lowest energy pathway to 1,3-PDO follows exactly the same pattern as that to 1,2-PDO except that the reaction is initiated from the primary C and the overall barrier increases by 0.2 eV. Obviously, the higher barrier in forming 1,3-PDO can be attributed to the fact that the aldehyde (glyceraldehyde) is less stable than the acetone form (dihydroxyacetone) and so does the TS for C-O bond dissociation. It should be mentioned that, as shown in ESI Fig. S7, † even if glyceraldehyde is produced, the presence of a low-barrier dehydration channel via the retro-Michael addition reaction will still favor the selectivity to 1,2-PDO.
Compared to the dehydrogenation pathways, the dehydration pathways ( Fig. 5b dotted lines) generally have a higher free energy barrier, 1.71 eV to 1,3-PDO and 1.77 eV to 1,2-PDO and are thus essentially non-selective to these two products. The higher barrier is mainly due to the poor ability of the protonated glycerol precursor (states 1 and 7 in Fig. 5a) on Cu(111). From our DFT calculations, the adsorbed state 2 in the dehydrogenation pathway is 1.07 eV more stable than the protonated glycerol, in which the vibrational zero-point-energy and van der Waals interaction together account for 0.35 eV energy difference.
Microkinetics simulations based on the DFT free energy prole are then applied to compute the steady-state rate, where the species coverages are determined explicitly to compute the rate. The computation is done by solving the ordinary differential equation (ODE) of all reactions using the Backward Differentiation Formula (BDF) method. 75 The simulation is under the typical reaction conditions, i.e. 473 K, p H 2 ¼ 1 bar, P others ¼ 0.0065 bar (other saturated glycerol derivatives). All the kinetics data are obtained from the above results and listed in ESI Table S6. † Our simulation results show that the major pathways of the 1,2-PDO or 1,3-PDO generation are dehydrogenation pathways, which contribute 2.1 Â 10 À4 s À1 and 2.3 Â 10 À7 s À1 for 1,2-PDO and 1,3-PDO, respectively, and the dehydration pathways are kinetically forbidden, which only contribute 1.6 Â 10 À9 s À1 and 6.9 Â 10 À9 s À1 for 1,2-PDO and 1,3-PDO, respectively. The nal TOF is 2.1 Â 10 À4 s À1 for 1,2-PDO and 2.4 Â 10 À7 s À1 for 1,3-PDO, which leads to a high selectivity of 99.88% to 1,2-PDO, in agreement with the experiment selectivity of 96.1%. 14, 16 We would like to point out that the previous proposed pathways to 1,2-PDO (see Fig. 1) are not favored kinetically: they do not show up as the lowest energy pathways from our AI-Cat search and are also conrmed as high barrier pathways by manually checking with DFT calculations. For example, the glyceraldehyde formation followed by the retro-Michael reaction as suggested by Montassier et al. 5 (black arrow in Fig. 1) has been a popular mechanism for 1,2-PDO in the literature. 13,21 However, our results show that glyceraldehyde formed by the dehydrogenation at the primary C is kinetically less favored than that at the secondary C in forming dihydroxyacetone by 0.22 eV (see Fig. 5b), although the subsequent retro-Michael addition reaction is kinetically feasible. The other proposed pathways include the formation of acetol (observed intermediate in experiments) as an intermediate, but those involving either 1,2 b-elimination dehydration or direct C-O bond breaking are non-site-selective, which are readily ruled out from our AI-Cat search due to their high barriers.

Discussions on glycerol hydrolysis
Our new reaction mechanism for glycerol hydrogenolysis on Cu catalysts can rationalize the experimental observations on Cubased catalysts reported so far, which are summarized as follows.
First, acetol is a stable product together with 1,2-PDO. Ample evidence from experiments 12,23 suggests that acetol is indeed generated from glycerol on Cu surfaces, instead of being a dehydrogenated product from 1,2-PDO. In particular, the selectivity could be dominated by acetol (90 mol%) when the reaction is operated in a N 2 atmosphere 23 (250 C and Cu/Al 2 O 3 catalyst). Indeed, from our results, acetol can be thermodynamically as stable as 1,2-PDO depending on H 2 pressure and thus can readily desorb from Cu surfaces. The high stability of acetol also leads to further C-O bond breaking (to 2-propanol) becoming kinetically impossible and the hydrogenation to 1,2-PDO is le as the only viable route (see ESI Fig. S8 †).
Second, dihydroxyacetone is a key intermediate in forming acetol and 1,2-PDO. By comparing glycerol conversions on Cu-ZnO (1 : 1) and Cu in a N 2 atmosphere, 21 Liu et al. found that both acetol and its hydrogenated form 1,2-PDO are produced in the absence of H 2 . This indicates that dehydrogenation should be involved in the catalysis process, which would otherwise lead only to the formation of the dehydration product (acetol) due to the lack of H supply. Our results conrm that the dehydrogenation steps are kinetically most favored on Cu surfaces and dihydroxyacetone is the precursor to acetol and 1,2-PDO.  19 From our results, the presence of strong acidity as introduced by H 4 SiW 12 O 40 /WO x helps to stabilize the protonated glycerol (states 1 and 7) and thus benet the nonselective direct dehydration pathways over the selective dehydrogenation pathways. In line with this, the high H 2 pressure that depresses the dehydrogenation channel should also be effective to increase 1,3-PDO selectivity. Indeed, we note that under acidic conditions, by replacing Cu with Pt, the 1,3-PDO selectivity can be much improved, 76,77 which is in line with the fact that the Pt surface can adsorb H stronger and achieve the saturated H coverage. This effectively blocks the 1,2-PDO selectivity that needs rst dehydrogenation. For example, the Qiao group reported that tungsten-doped siliceous mesocellular foam supported Pt catalysts (Pt/W-MCFs) can reach a high 1,3-PDO selectivity of 65% at 423 K and H 2 pressure of 4.0 MPa. 77 We note that the mechanism identied here for glycerol hydrogenolysis can be readily applied to understand the dehydration of other polyols. For example, 1, 2-butanediol is reported to be selectively converted to 1-hydroxy-2-butanone (>95%) over the CuZnCrAlO (CuO : ZnO : ZrO 2 : Al 2 O 3 ¼ 12 : 1 : 2 : 2) catalyst at 210 C under a N 2 atmosphere. 78 The selective dehydrogenation to ketone products is in line with our nding that the dehydrogenation which is kinetically favored selectively occurs at the secondary C and further dehydration is kinetically unfavored because of the high stability (similar to acetol). Similarly, 1,3-butanediol is found to be mainly converted to 4-hydroxy-2-butanone ($42%) and 2-butanone ($35%) over 25.0 wt% Cu-SiO 2 catalyst at 250 C under a N 2 atmosphere. 79 From our mechanisms, the main products should be produced from the selective secondary C-OH dehydrogenation (4-hydroxy-2-butanone) which is then followed by dehydration via the retro-Michael addition (2-butanone).

Low-energy C-O bond cleavage channels in general
Our AI-Cat conrms that despite the huge number of possible reaction intermediates from glycerol hydrogenolysis, the C-O bond cleavage must follow several low energy reaction patterns. This leads us to wonder how many low energy reaction channels are available for C-O bond cleavage on Cu surfaces, which should be of general signicance for understanding polyol catalysis on metals in general.
By analyzing our reaction database which contains 1948 reactions and 592 reaction patterns that are related to C-O bond breaking, we identied six low energy reaction patterns, which start from saturated or close-to-saturated molecules/fragments where the unsaturated bond is terminated by the surface. These six types of reactions are summarized in the le panel of Fig. 6 and ordered according to their reaction barriers on Cu surfaces. They are further elaborated in the following.
Types I and II: the most feasible C-O bond breaking (barrier < 0.5 eV) occurs on the C atom with dual C-O bonds. This includes COOH / CO + OH (type I) and R 2 C(OH)OR / R 2 C]O + HOR (type II).
Types III, IV and V: these reaction types feature a non-sp 3 C (a-C) neighbouring the C-O bond, which generally have the barrier of 0.5-1 eV. This unsaturated C can either be sp 2 C atom with a double-bond or a C radical site, which can adsorb on metal surfaces and is thus stabilized as a reaction intermediate. These a-C are in analogy to the allylic site or the a-site of carbonyl in organic chemistry, and the presence of catalysts stabilizes the site and makes these patterns widely present in the reaction database. Among them, both III and IV are typical examples with asite sp 2 C. Type III, RC(]O)CH(R)CH 2 (OR) / RC(]O)C(R)] CH 2 + *H + *OR, is the retro-Michael addition, also found in the glycerol hydrogenolysis reaction network (glyceraldehyde to 2hydroxyacrylaldehyde); type IV, RC(]O)CH 2 (OR) / RC(-O*)] CH 2 + *OR, is the keto-enol tautomerization facilitated C-O bond breaking, which is the key reaction in glycerol hydrogenolysis (dihydroxyacetone / 6); type V reaction, R 2 C-*CH 2 (OR) / R 2 C]CH 2 + *OR, can be considered as an a-C radical pattern appearing in the base catalyzed E2 elimination.
Type VI: the direct C-O bond breaking can occur by a rst hydrogenation followed by dehydration with the barrier being