Short history: From DFT to machine learning

Heterogeneous catalysis is renowned for its great complexity in catalyst structure, and thus revealing how reactions occur on catalyst surfaces is an outstanding challenge. A major difficulty stems from the intimate coupling between catalyst surfaces and molecules during the catalytic conversion1,2. The catalyst surface may well reconstruct caused by the molecular adsorption, and molecules can in turn choose the best surface sites to achieve the highest reaction kinetics. Therefore, new techniques that can be operated under reaction conditions and have high spatial-temporal resolution were continuously pursued to characterize the in situ catalyst structure and to reveal the reaction mechanisms. However, to date, most experimental techniques remain frustrating to work properly under even the most common catalytic conditions, such as high-pressure conditions3,4 and solid-liquid reaction conditions5,6,7. On the other hand, atomic simulations, particularly those based on quantum mechanics (QM) calculations, have attracted much attention in the past two decades: they are not only an indispensable complement to experimental characterization techniques, but also could provide insightful predictions to guide experimental catalyst search8,9. It is the purpose of this review to introduce the advance of the latest machine learning atomic simulations in heterogeneous catalysis.

Theoretical calculation emerges as a key player in catalysis research as early as the 1980s. With the advent of density functional theory calculations for periodic systems in the 1990s10,11,12, the key obstacle in computing extended solid surfaces was removed, and the first principles atomic simulations soon became a popular and trustful tool in the community. Compared to empirical force field calculations and the finite system DFT calculations, the plane-wave DFT calculations for periodic systems can reach much better accuracy in computing the potential energy surface (PES) of surfaces and the interaction between molecules and surfaces13,14,15. This offers the possibility to systematically compare the adsorption and reaction of different molecules on different surfaces. Because of the relatively low scaling (O(NlnN)) of DFT calculations, the systems under 100 atoms (in periodicity) can now be routinely performed on modern workstation computers. To date, tremendous progress thus has been achieved in catalysis via atomic simulation, where perhaps all known important catalytic systems have been explored by using the single-crystal bulk structures and well-defined low Miller index surfaces as the model16,17,18,19,20,21,22. Even though, DFT calculations appear still not a game changer for at least two reasons. First, the catalytic systems are generally much more complex than single-crystal surfaces, which involve many more atoms and more complex surface geometries23,24,25. Second, the computation of chemical reactions that need to locate the reaction transition state (TS) is typically one or two orders of magnitude more expensive than that of the adsorption state (initial and final states)26,27,28,29,30,31. Therefore, the predictive power of DFT calculations for catalytic reactions is limited to even smaller systems, not to mention the concerned DFT intrinsic errors in computing the reaction kinetics.

To allow for efficient PES exploration as required, for example, to identify the best catalytic site, many elegant theoretical methods have been proposed. While the real catalyst structure evolution must be in contact with the reaction environment (e.g., gas molecules), most theoretical simulations have to work with fixed numbers of particles in the system (constant number of atoms). For example, molecular dynamics (MD) can produce a structure evolution trajectory in real-time on a finite-temperature free energy surface32,33, and enhanced MD methods can further speed up the exploration of designed chemical reaction events34,35. Unfortunately, due to the speed limit in PES evaluation, the time scale of MD simulations is generally below nano-seconds, which is far not enough for finding the most stable structure at a given composition, a global optimization problem. To this end, a number of global optimization methods were proposed since 1990 to break the time scale of MD simulation, which abandons the strict detailed-balance rules and the finite-temperature free energy effects in kinetics and focuses on finding the global minima (GM) by randomized structure perturbation. These methods include basin hopping36,37, genetic algorithm (GA)38,39,40, particle swarm optimization (PSO)41,42, and the stochastic surface walking (SSW) method43,44,45. As a complement, variable-atom global optimization methods have also been developed, which can find the most stable structure at different compositions46. Even with these advances, these current global optimization methods still cannot replace the grand canonical Monte-Carlo (GCMC) method47,48,49,50 that realizes the dynamic atom exchange under the fixed chemical potential, despite the fact that the GCMC method with the atom-wise operation in structure evolution is generally too expensive to apply to catalytic systems.

Another attractive feature of atomic simulations is the feasibility to locate the lowest energy reaction pathway, which is often hard to conclude from the experiment due to the transient nature of intermediates. The main stream of atomic simulation methods for finding reaction pathways is based on the Markov chain theory and transition state theory51,52, where a chemical process is considered to be composed by many elementary steps and the rate of each elementary step is controlled by the barrier height between its TS and the initial state (IS). For multiple-step chemical process, the mean-field microkinetics53,54 or kinetic Monte-Carlo (kMC) simulations55,56,57 are required to transfer the knowledge of the microscopic barrier of elementary reactions to macroscopic kinetics of catalysis. While the barrier is typically determined by TS location method as the energy barrier at zero K26,27,28,29,30,31, the free energy barrier of an elementary step could be computed by enhanced MD-based approaches once the reaction coordinate is known, such as the umbrella sampling (US) method that integrates out the potential of mean force along a designed reaction coordinate58,59,60. Nevertheless, all of these methods require a priori knowledge or pre-guess on the reaction mechanism to obtain the reaction coordinate for each elementary step. Compared to structure exploration, the finding of reaction pathways is thus not only much more demanding in computation but also needs a priori knowledge as input from experienced users.

In recent years, machine learning (ML) based atomic simulations are evolving rapidly, achieving huge progresses on both the methodology for PES evaluations to the algorithms for structure and reaction pathway sampling61,62. In particular, the latest neural network (NN) potential calculations can be more than 104 faster than DFT calculations without significant loss of accuracy63,64. These ML atomic simulations bypass the heavy QM calculations and utilize ML models to link the atomic coordinates with the total energies by learning the PES data from QM calculations64. These new techniques are bound to reshaping profoundly the research in heterogeneous catalysis61,62,65,66.

As a representative of ML potentials, the high-dimensional neural network (NN) proposed by Behler and Parrinello is one of the most widely utilized ML framework64,67,68. They propose that the total energy of a system is the sum of each atomic energy, which allows to learn the total energy, an extensive quantity, by using atom-wise NNs. The input layer of atomic NN is then a set of permutation-, translation- and rotation-invariant numerical functions, the so-called atom-centered symmetry functions, that can sensitively reflect the local chemical environment of atom. To date, there are many other flavors ML models, although most are still based on the same atom-wise NN framework68. The structural descriptors and thus the input layer of NN are generally the most different parts among different ML potentials, for example, the convolutional neural network (CNN) and graphic neural network (GNN) are utilized to distinguish the chemical environment of atom69,70,71,72. In our group, we proposed a set of power-type structural descriptors (PTSDs) as the structural descriptors for fitting the global PES data, where the spherical functions and the four-body terms are introduced73. The complex PTSD descriptors allow the ML atomic simulation in complex material systems, such as boron cluster73, zeolite66,74,75, phase interface76, and heterogeneous reactions77,78,79.

In this review, we will highlight the recent progress on ML-based methodologies for solving complex PES problems in heterogeneous catalysis, both on the structure determination and reaction pathway finding. We finally discuss the future directions in theoretical catalysis, focusing on large-scale atomic simulations.

Atomistic simulations for exploring structures

Since catalysts are dynamic, reconstructing and interacting constantly with coming molecules under catalytic conditions, atomic simulations to explore unknown catalyst structures ideally need to be held under the grand canonical (GC) ensemble, which means particles in a system can exchange with the environment as driven by the chemical potential80,81. Although the grand canonical Monte-Carlo (GCMC) simulations were developed for such a purpose47,48,49,50, it is difficult to combine GCMC with QM calculations due to the slow self-consistent-field cycles in solving Schrodinger equations. Additionally, the MC method as utilized in GCMC evolves structure via random perturbation, which is often not efficient enough to find stable configurations. Here we will introduce the progress of ML-powered global optimization and grand canonical structure search for catalyst structures, which are benefited significantly by the high speed and high accuracy of ML potentials.

Machine learning global optimization

The global optimization method was introduced to search for unknown structures from 1990s82. Compared to the typical MD and MC methods, the global optimization methods can better surmount high barriers between minima and avoid trapping at high-entropy regions of PES. Due to the dramatic structure change between configurations and the high computational costs, the global optimization methods have long been limited to a few toy models and small clusters83.

Unlike most global optimization methods, the stochastic surface walking (SSW) method developed in our group utilizes small structure perturbations as driven by bias potentials to move a structure from one minimum to another43,44,45. The transition region between minima can be properly visited by SSW which allows the reaction pathway searching during the global optimization. The early applications of SSW include both the finding of the global minimum of clusters43, and the phase transition pathways between bulk crystal materials45. The SSW method provides a powerful and convenient solution for PES exploration and thus becomes an ideal tool for PES data generation. Indeed, in 2017 a global-to-global scheme as introduced by our group to combine SSW with NN endows the global optimization ability to ML potentials84. The so-called global NN potential (G-NN) have good transferability and can be utilized to explore unknown structures with arbitrary compositions starting from random structures. The SSW-NN methods have been applied to a wide range of PES problems, in particularly those related to complex surface and interface structures in heterogeneous catalysis85.

Figure 1a gives an example of how SSW-NN determines the famous silver surface oxide structure, Ag12O6 phase on p(4×4) periodicity of Ag(111)86,87,88. First, a Ag-O global NN potential is generated by SSW-NN method, which is iteratively trained by learning a wide range of different Ag-O structures, including the bulk and surfaces for metal Ag and AgOx oxides. In searching for the surface oxide structure, the initial structure (‘IS’ in the figure) can randomly generated, for example, from a 1 ns of high-temperature (2000 K) NN-based MD simulations that contains O2 species, surface oxygen, subsurface oxygen, and silver vacancies. Next, the SSW-NN is conducted to locate the global minimum (GM) of this structure. A 2000-step SSW-NN trajectory at a given Ag12O6 composition is shown in Fig. 1a, where the energy drops rapidly in the early stage, for example, the structure in step-372 is only 0.2 eV higher than the GM. The SSW trajectory also reveals the stable structure patterns of surface oxide, namely (i) the subsurface oxygen is not stable (str 1 → 2, red part), (ii) molecular O2 species, if present, will dissociate on the surface (str 2 → 3, green part), and (iii) the optimal Ag-O arrangement (str 3 → 4, blue part). Finally, the GM structure, identical to that reported previously87,88, is identified after 1979 SSW steps, which is verified by continuing SSW search for another 10,000 steps and no more stable structures is found.

Fig. 1: The identification of silver surface oxide.
figure 1

a The 2000-step SSW-NN trajectory for identifying the structure of silver surface oxide (Ag12O6 phase on p(4 × 4) supercell of Ag(111)). The initial structure (IS) and four representative structures under different stage of optimization during SSW-NN are shown. b The building of Ag12O6 optimal surface phase with the given input of p(4 × 4) supercell and Ag12O6 composition by reinforcement learning based atomistic structure learning algorithm, in which the depth of blue in the figure represents the Q-values (the preference for the location of the next atom) predicted by the CNN, reproduced with permission from ref. 89. Copyright American Institute of Physics, 2019.

Full size image

For the same Ag-O surface oxide system, Hammer group has developed the atomistic structure learning algorithm (ASLA) to identify the global minimum, as shown in Fig. 1b89,90. The method builds 2D structures using an atom-by-atom strategy based on a reinforcement learning framework. In this method, CNN is used to acquire the knowledge of a surface structure and output the Q-values used for locating the next atom in ε-greedy policy. DFT calculations are performed to evaluate the stability of the predicted structures. It takes more than 3990 episodes (one trial for building Ag12O6 phase is one episode) to resolve the Ag12O6 global minimum phase. Apparently, if DFT can be replaced by ML potentials, the searching speed should be dramatically dropped. This demonstrates, on the other hand, that ML techniques can be combined versatility with PES evaluation methods to solve complex PES problems.

ASOP algorithm

The fast global optimization ability of SSW-NN opens the door for searching structures with variable compositions in a GC ensemble. We have developed the automated search for optimal surface phases (ASOP)91, which can explore the composition space under a predefined chemical potential condition. As shown in Fig. 2a, the ASOP simulation is based on a multi-grid framework to scan all likely compositions. It takes only simple inputs, including the bulk crystal structure, the surface Miller index, and the chemical potentials of exchanging particles, and outputs a phase diagram, including a list of stable phases.

Fig. 2: The ASOP algorithm and its application.
figure 2

a Flowchart of the ASOP algorithm that contains three steps. b The PES contour map and energy spectrum for silver surface oxides on Ag(100) from ASOP simulation. c The geometries of two phases, as labeled in both PES map and energy spectrum. Reproduced with permission from ref. 91. Copyright American Institute of Physics, 2022.

Full size image

In one ASOP simulation, the whole composition space is discretized into a series of grids from coarse to fine, in which each grid represents a unique surface periodicity ordered by the surface area from small to large. Then, several cycles of SSW-NN structure exploration and Monte-Carlo selection are conducted for the compositions in each grid. The simulation starts from the coarsest grid and progressively explores the larger grids, where the knowledge of the stability of compositions is inherited from grid to grid. The SSW-NN sweeps over each composition with only a few steps (e.g., smaller than 400) since the energetically favorable compositions will be biased and visited again in the subsequent cycles. The Monte-Carlo scheme is utilized to select the focus zone containing energetically favorable compositions.

Taking Ag(100) oxidation under typical ethene epoxidation conditions (500 K, 1 bar of O2) as the example, the ASOP algorithm can identify the top stable silver surface oxides within 90 hours on 80 CPU cores. Figure 2b shows the PES map and energy spectrum of surface oxides, in which the optimal zones locate at the Ag and O coverages being 0.6~1.0 and 0.5~0.9 ML, respectively (see the blue zone in the figure). Figure 2c plots the geometries of phase-1, the most stable structure from ASOP, which differs from an experimentally observed surface oxide at low O2 pressure, phase-292,93. The phase-1 possesses a Ag7O5 stoichiometry with O coverage of 0.625 ML, which contains the planar coordinated [AgO4] motif with high O-density. By contrast, the phase-2 has a lower O coverage (0.50 ML) and is less sable (0.025 J/m2, see energy spectrum of Fig. 2b), which features a missing-row reconstruction induced by atomic O adsorption. These results indicate that Ag(100) surface under typical ethene epoxidation conditions are in fact at a deeper oxidized state comparing to the surface science experiment. The Ag catalyst surface is highly dynamic under different O2 pressures.

Reaction activity prediction

The determination of the catalyst structures is generally not enough for predicting reactivity in heterogeneous catalysis. The molecules can adsorb, diffuse and react at a vast number of local sites, which may only have trivial differences in geometry but can provide very different reactivity. Traditionally, the reaction pathway at each site needs to be determined by locating all likely transition states (TS) in order to predict correctly the overall activity51,52. The TS search can be accelerated by the NN potential, but the pre-guess of reaction coordinate remains to be a key obstacle in unknown reactions. The Brønsted-Evans-Polanyi (BEP) relationship94,95,96, that correlates the reaction energy (adsorption energy difference between the initial and the final states) with the reaction barrier is often utilized to speed up the reactivity prediction, although the accuracy of BEP is generally poor and may only be used as a coarse screening tool.

Recently, Sun and Sautet demonstrated further the complexity of reaction pathways on nanocatalysts using a concept of catalyst fluxionality97,98,99. They utilized the HDNN potential together with a modified genetic algorithm to explore extensively the low energy metastable ensemble (LEME) of Pt13Hx clusters and then evaluates the hydrogen evolution reaction (HER) and methane activation activities on these clusters97. Fig. 3a shows the top three stable Pt13H26 clusters in its LEME, in which the Pt atoms exhibit various coordination numbers, such as four-coordinated Pt (PtH4, see GM0, GM1, and GM2), three-coordinated Pt (PtH3, see GM0 and GM1), and two-coordinated Pt (PtH2, see only GM1). The kinetic evaluation of all these clusters in Fig. 3b indicate the PtH2 site is the most active site with the lowest activation energy, which locates on energetically less favored GM1 compared to GM0. Notably, even with the consideration of presence probability, the weight reaction rate of GM1 is still 30 times higher than GM0, indicating the fluxionality of active site from the most stable structure to metastable isomer, which can hardly be noticed without the machine learning accelerated large-scale structure sampling.

Fig. 3: Pt clusters for methane activation.
figure 3

a The geometries of top three stable Pt13H26 clusters (denoted as GM0, GM1, and GM2) as explored by the HDNN potential together with modified genetic algorithm, in which the note ‘PtHx’ represents the coordination number (x) of Pt. b The activation energies (Ea) and relative contributions to the reaction rate of methane activation under 673 K. Reproduced with permission from ref. 97. Copyright American Chemical Society, 2018.

Full size image

The ultimate goal of reaction activity prediction in the heterogeneous reaction is to discover reaction channels in an automated way, which should take into account the intimate coupling between catalyst structure and molecular reactions. Since the complexity of the reaction network grows exponentially when many elementary reactions and likely reaction sites are present, an intelligent on-the-fly pathway explorer is essential in order to identify the kinetically relevant reactions by the algorithm itself. For this purpose, new methods have developed by exploiting the SSW global optimization for PES search and the G-NN potential for PES evaluation.

SSW-RS method

The SSW reaction sampling (SSW-RS) method modifies the SSW global optimization for pathway collection during the SSW PES search. Specifically, in the SSW-RS simulation, the structure characteristics, such as the bond matrix and chirality of the initial structure, are remembered and used to judge whether a new structure from SSW is identical to the initial structure. If it is false, SSW will output a reaction pair for the subsequent pathway search and the current structure is back to the initial structure to continue the SSW search. Finally, the TS location for all reaction pairs is performed by the double-ended surface walking (DESW)31, which can identify iteratively all TSs along the pathway connecting a reaction pair. The SSW-RS method can be combined with DFT and G-NN calculations for studying different type of reactions, from gas phase reactions100,101,102 to solid phase transitions103,104,105.

Figure 4(a) illustrates an SSW-RS example for finding the reaction pathways for the oxometallacycle (OMC) intermediate on Ag(100), which is known to be a critical intermediate in ethene epoxidation106,107,108. It was generally believed that OMC can further oxidize to AA and EO with similar barriers on metal Ag sites and thus the EO selectivity on Ag catalyst is about 50 %106,109,110.

Fig. 4: The SSW-RS method and its application.
figure 4

a The barrier-energy map of the sampled products of oxometallacycle (OMC) intermediate on p(4 × 4) supercell of Ag(100) after 10000-step SSW-RS. Abbreviations: OxoE, 2-oxoethyl; AA, acetaldehyde; VA, vinyl alcohol; EO, ethene oxide. b The diagram of the general SSW-RS based pathway sampling algorithm, reproduced with permission from ref. 79. Copyright American Chemical Society, 2021.

Full size image

By using SSW-RS, a total 864 qualified reaction pairs are collected from 10000 SSW-RS steps, and after DESW pathway search, 29 distinct products are identified. In these products, only five of them are kinetically favored with a barrier smaller than 1 eV (see Fig. 4a), i.e., OxoE (2-oxoethyl), AA (acetaldehyde), VA (vinyl alcohol), ethene, EO (ethene oxide). The SSW-RS identifies a new intermediate OxoE and much lower formation barrier of AA than EO, in which the route from OMC to AA should be indirect via OMC → OxoE → AA79. This finding rules out the metal Ag sites as the active site for EO production, and the active site for EO production must be Ag surface oxides.

For multiple-step reactions, SSW-RS can be iteratively performed and the whole reaction network can be resolved finally, as shown in Fig. 4b79. Starting from a given reactant A, the first SSW-RS is performed to sample the possible reaction pairs connecting to the reactant A and locate all the TSs, then the reaction network is updated based on the sampled new reaction channels. Next, according to the order of the calculated barriers, the kinetically favorable intermediates are chosen to start new SSW-RS tasks. By performing SSW-RS to scan all likely intermediates, the lowest reaction pathways from reactant A to any products can in principle be obtained in an automated manner. This simple algorithm meets difficulties when the reaction network is too complex with too many elementary reactions. The AI-Cat and MMLPS algorithms are thus developed to cope with the challenges relating to complex reaction networks.

AI-Cat algorithm

Kang et al. proposed an end-to-end artificial intelligence framework for the activity prediction of heterogeneous catalytic systems (AI-Cat)111. The AI-Cat method is a general framework for predicting the kinetics of heterogeneous catalytic processes where the reaction network is too complex to resolve by using SSW-RS. The method first needs a general reaction database containing information on common elementary reactions, including the structures and energetics of initial state (IS), final state (FS), and their TS. The database can be obtained by the SSW-RS simulations introduced above. The structures in the reaction database are then encoded by Surface-sensitive atom-centered Extended-Connectivity FingerPrint (s-ECFP). To learn the reaction information, two neural networks, namely the reaction pattern (R-Pat) unit and the kinetics information (K-Info) unit, are built for predicting the most likely reaction patterns (RPs) (i.e., the coding of IS/FS pair) from a given IS and predicting the kinetics data (reaction barrier and energy) of this IS/FS pair, respectively. The reaction pathways between any given reactant and product are determined using the Monte-Carlo tree search method, which inquires the R-Pat unit and the K-Info unit in making decisions.

Figure 5a shows a MC tree search step, namely the expansion step, from a father node to child nodes. The input reactant, the father node, firstly go through R-Pat unit to generate all the possible RPs. Next, the top RPs are selected and go through K-Info unit to predict the associated reaction energies and barriers. Finally, the kinetically favorable products, the child nodes are determined.

Fig. 5: The AI-Cat method and its application.
figure 5

a The expansion from father nodes (reactants) to child nodes (products) in AI-Cat method. b The Gibbs free energy profile for four low energy pathways of glycerol hydrogenolysis to 1,2-PDO and 1,3-PDO under typical experimental conditions (473 K, total pressure of 1 atm with glycerol: H2: H2O = 1: 140: 12). Reproduced with permission from ref. 111. Copyright Royal Society of Chemistry, 2022.

Full size image

By using SSW-RS to collect the reaction data for common C1-C4 organic compounds on Cu(111), Cu(100), and Cu(211), and training R-Pat unit and K-Info unit based on the database, the AI-Cat method is able to explore the glycerol hydrolysis pathways on Cu(111), which is significant in biomass-derived polyol utilization112,113. Fig. 5b shows the Gibbs free energy profile for the AI-Cat predicted lowest-energy pathways from glycerol to 1,2-propanediol (1,2-PDO, blue solid line) and 1,3-propanediol (1,3-PDO, red solid line). In the rate-determining step (TS1 and TS2), the overall barrier of 1,2-PDO path is 0.22 eV lower than 1,3-PDO path (TS2 – (TS1 + state 2)), indicating the hydroxyl on the terminal C favors dissociation more than central C. Therefore, the 1,2-PDO product after TS1 is more selective than 1,3-PDO after TS2, thus explaining the long-standing high selectivity puzzle of 1,2-PDO on Cu surfaces114. It is worth noting that the popular glyceraldehyde mechanism in literature that involves a dehydration process (TS3 and TS4)115,116, is neither kinetically favored nor able to explain the high 1,2-PDO selectivity (see TS1, TS2 < TS3, TS4 and TS4 > TS3 in the figure).

MMLPS algorithm

Shi et al. proposed a microkinetics-guided machine learning pathway search method (MMLPS) method to speed up the buildup of the reaction database by SSW-RS117. The MMLPS aims to fast build a reaction database for a target reaction and identify the kinetically favorable pathway. Compared to AI-Cat method which is more general and can predict unknown reactions, the MMLPS is the tool to get the accurate kinetics for a target reaction. Specifically, the MMLPS simulation divides a reaction network into several parts with different molecules and surface coverages, then each SSW-RS branch samples independently different parts of the reaction PES as guided by a fast microkinetics solver. A reaction dataset is established by merging reactions from all branches, and the microkinetics simulation is performed to identify the lowest barrier reaction pathway.

Figure 6(a) illustrates the complete 2D reaction map of CO and CO2 hydrogenation on Cu and Zn-alloyed Cu surface, as plotted from 14958 reaction pairs sampled by MMLPS. These reaction pairs are collected from three batches, i.e., CO2 + H2, HCOOH + H2, and HCHO + H2, which depicts the whole reaction channel from CO/CO2 to methanol. On all surfaces, CO2 hydrogenates via the formate pathway (CO2 − HCOO* − HCOOH* − H2COOH* − HCHO* − CH3O* − CH3OH* − CH3OH), and CO hydrogenates pathway via the formyl pathway (CO − CO* − CHO* − HCHO* − CH3O* − CH3OH* − CH3OH), as shown by the free energy profile in Fig. 6b, c. On Cu(211), only 1.40 eV is required for CO2 hydrogenation, which is 0.05 eV lower than that of CO, indicating that CO2 is the main carbon source in methanol products. The microkinetics based on MMLPS data reveals that Zn alloying has no obvious kinetics promotion to CO2/CO hydrogenation reaction and a high coverage of Zn would even poison the catalyst (yellow and blue lines in the figure).

Fig. 6: The MMLPS method and its application.
figure 6

a Contour plot for 14958 reaction pair (IS, TS, FS) obtained by MMLPS on Cu(211). b CO2 and c CO hydrogenation Gibbs energy profile on Cu and Zn alloyed Cu surface. Reproduced with permission from ref. 117. Copyright American Chemical Society, 2022.

Full size image

Conclusions

This review outlines recent advances in ML potential-based atomic simulations for heterogeneous catalysis. The high accuracy and high speed of ML potentials accelerate greatly the PES exploration and thus allow the development of new algorithms to solve long-standing challenges in heterogeneous catalysis.

In particular, we illustrate a few key advances in realizing ML potential atomic simulations, including the atomic-based ML model to fit the total energy that is assumed to be a sum of single atom energies, the structure descriptors to distinguish the chemical environment of atoms, and the global-to-global scheme to generate global PES dataset and train ML potential. As the representative, the SSW-NN method and the methods developed upon SSW-NN, such as ASOP, AI-Cat, and MMLPS, can now tackle complex catalysis problems, such as the surface phase diagram under the grand canonical ensemble, the end-to-end reaction prediction based on the known reaction database, and the automatic determination of the reaction kinetics in a complex reaction network.

Nevertheless, the current algorithms, even with ML potential calculations, generally only consider the weak coupling between catalyst structure and reaction by either focusing on the thermodynamics in structure evolution or resolving the reaction kinetics on well-defined surfaces. This is apparently due to the too-large degrees of freedom if both surface flexibility and the reaction varieties are treated at the same time. In some catalytic systems, this approximation may lead to misleading or even wrong conclusions, for example, the reactions on nanoclusters as shown in reactions on supported Pt13 cluster. Ag-catalyzed ethene epoxidation is another example: despite the metal Ag is now ruled out as the active site for epoxidation, the true active sites of Ag surface oxides remain largely unclear, which is strongly influenced by the reaction atmosphere (ethene and O2 pressures).

One possible and feasible solution to treat strong-coupling systems could be an automatic focus on transient structures during PES search, including the meta-stable phases and the TS structures of molecular reactions. This is likely achieved by the on-the-fly constrained PES searching procedure98,100, which identifies the critical reaction coordinate and the lowest energy reaction channel by constraining PES search at the target reaction coordinate (e.g., important chemical bonds). During such exploration of reaction pathways, the surface phases should be allowed to change by considering the grand canonical reaction conditions. Given the recent progresses, one can be sure that the methodology advance and better applications of ML techniques in heterogeneous catalysis are coming in the future.