Heterogeneous catalysis from structure to activity via SSW-NN method

Modern research on heterogeneous catalysis calls for new techniques and methods to resolve the active site structure and reaction intermediates at the atomic scale. Here, we overview our recent progress on large-scale atomistic simulation via potential energy surface (PES) global optimization based on neural network (NN) potential, focusing on methodology details and recent applications on catalysis. The combination of stochastic surface walking (SSW) global optimization and the NN method provides a convenient and automated way to generate the transferable and robust NN potential for global PES, which can be utilized to reveal new chemistry from the unknown region of PES with an affordable computational cost. The predictive power of SSW-NN is demonstrated in several examples, where the method is applied to explore the material crystal phases, to follow the surface structure evolution under high pressure hydrogen and to determine the ternary oxide phase diagram. The limitations and future directions to develop the SSW-NN method are also discussed. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5113673., s


I. INTRODUCTION
Heterogeneous catalysts are renowned for their great complexity of material composition and surface structures. There has been a long history for both experiment and theory in the struggle to resolve the active site of catalysts that are responsible for activity and selectivity. [1][2][3][4] The latest progress in experiments are represented by high spatial resolution techniques, e.g., the spherical aberrationcorrected high-resolution transmission electron microscope, and synchrotron-based measurements, e.g., extended X-ray absorption fine structure. On the other hand, as a direct tool to correlate the atomic structure with its energy, theoretical simulations have been widely practiced in modeling catalyst structures and even predicting the activity, especially with the advent of density functional theory (DFT) calculations. 5 However, DFT simulation is generally limited to hundreds of atoms and thus fails to explore exhaustively the phase space of complex catalysts, such as amorphous structures, multicomponent oxide/alloys. [6][7][8] To identify the active site of catalysts from theory, it is essential to explore the potential energy surface (PES) of exposed surfaces.
This requires not only a fast and reliable approach to evaluate the energetics of structures but also an efficient method to explore the structural phase space. As for PES calculation methods, most are based on quantum mechanics as represented by DFT calculations, while the empirical force field calculations and more recently an artificial neural network (NN) 9 potential method ( Fig. 1), despite their limitation in transferability, are also often utilized in material applications. Since catalytic conversion involves the chemical bond making and breaking, DFT calculations have been the most popular way to provide an accurate description of reactions with reasonable computational cost. The computation speed of DFT is sensitive to the complexity of the employed density-functional (e.g., PBE 10 and HSE06 11 ) and has a poor scaling [at least O(NlnN)]. It is very difficult to explore the reaction network in large catalytic systems (e.g., >100 atoms). As a promising alternative, the NN potential method developed in the past decade demonstrates its power in treating complex PES problems, from the gas phase reactions to material dynamics. [12][13][14][15] Recently, it was also utilized for solving the structures of heterogeneous catalysts, for example, Pd(O), 16 Pt(H), 17 CuAu, 18,19 CuCeO, 20 and CuZnO. 21 Unlike the traditional force field FIG. 1. Scheme of the SSW-NN. The total energy can be obtained from the NN potential, which replaces the high cost procedure to solve the Schrödinger equation. The total energy of a system in the NN potential is a sum over all atomic energy. The atomic energy involves power-type structural descriptor (PTSD) computation (see Sec. II C 2) and a NN evaluation. The SSW global optimization is utilized to explore the PES and produce new structures. potentials, the NN potential is able to describe chemical reactions with high accuracy as long as the training dataset contains the reactive data, such as the transition state (TS).
In parallel to the progress in PES evaluation methods, many PES exploration methods have been developed through the past decades. As catalytic reactions and the structural reconstruction occur generally above ambient temperatures involving high barrier processes, traditional molecular dynamic (MD) simulation is often not appropriate for active site identification. Instead, the global optimization methods, which can overcome the high barrier on PES, are desirable. The common global optimization methods include simulated annealing, 22 genetic algorithm, 23 basin hopping, 24 and stochastic surface walking (SSW). 25,26 This perspective serves to outline our recent contributions in methodology development toward catalyst PES exploration and active site identification. We will show that the combination of the global NN potential with the SSW method, SSW-NN, provides a powerful platform to resolve the catalyst structure which can finally lead to the prediction of catalyst activity from first principles. The SSW-NN method together with other common atomic simulation techniques is now implemented in LASP software, Large-scale Atomic Simulation with neural network Potential, which allows one, within a user-friendly platform, to perform first principles calculations, NN potential generation, and atomic simulation using the NN potential. 27

A. SSW method
The SSW algorithm 25 is an unbiased global optimization method that can explore both minima and saddle points on PES.
SSW implements an automated climbing mechanism to manipulate a structural configuration moving smoothly from a local minimum to a high-energy configuration along one random mode direction. The method was initially developed for aperiodic systems, 26 such as molecules and clusters, and has been extended to periodic crystals. 28 The SSW method inherits the idea of bias-potential driven constrained-Broyden-dimer (BP-CBD) method for TS location. In one particular SSW step, labeled as t, a modified PES V mod , as shown in Eq. (1) and Fig. 1, is utilized for moving from the current minimum, R 0 t , to a high energy configuration, R H t (the climbing), in which a series of bias Gaussian potential vn (n is the index of the bias potential, n = 1, 2, . . ., H) is added one by one consecutively along the direction N n t (Fig. 1), where Rt is the coordination vector of the structure and V real represents the unmodified PES; R n t is the n-th local minima along the movement trajectory on the modified PES that is created after adding n Gaussian functions. The Gaussian function is controlled by its height w and its width ds and is always added along one particular walking direction as defined by N n . Once the R H t is reached, all bias potential are removed and local optimization is performed to quench the structure to a new minimum. Different from the BP-CBD method, each SSW step (from one minimum to another) will choose a random direction to perturb the structure after the direction is refined (softened) using the biased-CBD method. At the end of the SSW step, a structure selection module, e.g., in the Metropolis Monte Carlo scheme, is applied to accept/refuse the new minimum.

B. High-dimensional NN architecture
The artificial NN method was first developed to understand the signal processing in the brain. 29 In the following decades, NN evolves into a class of powerful algorithms applied to a variety of fields from numerical prediction, pattern recognition to data classification. It is most renowned for the powerful ability to establish the functional relationship between independent variables and target (dependent) values via nonlinear "black box" data processing.
In 2007, Behler and Parrinello implemented the highdimensional NN (HDNN) architecture for material simulation, which trains the first principles PES dataset to generate NN potential. [30][31][32] In the approach, the total energy Etot of system is written as the summation over all atoms [Eq. (2)] (see Fig. 1). Each atom is represented by a standard feed-forward (FFNN), where the input layer is a series of structural descriptors to describe the atom bonding environment and the output layer yields the atomic energy Ei (i indexes atoms), The structural descriptors can, in principle, be any functional forms of atomic coordinates. [33][34][35] In their original work, Behler et al. suggests a series of rotation-invariant symmetry functions as the The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcp structural descriptors, which are Gaussian type functions of intrinsic coordinates (pair distance, angles). 30 The NN potential can be trained by minimizing the cost function that measures the deviation of the NN output with respect to the training set properties, e.g., total energy, force, and stress. 36,37 The training set is a large structure dataset on PES with accurate energetics and forces, most often computed from first principles calculations. More details on HDNN can be found in previous reviews. 35,38,39 C. SSW-NN

Self-learning procedure
Unlike classical force field potentials, the NN potential is generally found to have a limited predictability beyond the training dataset, presumably due to the numerical function fitting by a large number of parameters. The key to improve the quality of the NN potential therefore relies much on the representativity of the PES dataset. To overcome this deficiency of the NN potential, we propose in 2017 a global-to-global scheme to generate the global NN (G-NN) potential for material simulation. 37 This scheme combines the SSW global optimization for data generation, the HDNN for PES description, and a self-learning procedure 17 to expand the dataset and upgrade the NN potential, the so-called SSW-NN method. The self-learning procedure of SSW-NN is described briefly as follows (see Fig. 2).
The first stage constructs an initial dataset by performing in parallel short-time SSW sampling based on first-principles DFT calculations. These DFT calculations are often restricted to small systems (typically below 20 atoms) and with a low accuracy calculation setup to speed up the global PES sampling. After the PES data are obtained from SSW, a small dataset is randomly selected and computed using DFT with a high accuracy calculation setup. This stage produces a dataset with the most common atomic environment for the target PES.
The second stage generates a NN potential using the current first principles dataset. The key features of our NN potential will be detailed in Subsection II C 2.

FIG. 2.
Self-learning procedure of the global NN potential. The global dataset is first generated using high accuracy DFT calculations, which is then trained to obtain the global NN potential (G-NN). Then, an additional dataset is generated by SSW sampling based on the G-NN potential. This additional dataset is then fed back into the global dataset and a new cycle self-learning starts.
The third stage expands the dataset by carrying out long-time SSW global PES sampling using the current NN potential. These SSW-NN simulations start from a variety of initial structures with different morphology, including bulk, surface, and clusters, different chemical compositions, and different number of atoms per unit cell. A small additional dataset is thus obtained from the SSW sampling trajectories, containing the structures on PES either randomly selected or exhibiting new atomic environment (e.g., out-of-bounds in structural descriptor, unrealistic energy/force/curvature). After calculating these additional data by DFT, they are added into the training dataset and the whole self-learning procedure returns back to the previous stage.
Typically, after ∼100 iterations, a robust and accurate NN potential can be obtained with a compact training set that contains the most representative structures. We emphasize that it is better to adopt consistent and high accuracy calculation setups in constructing the training dataset using first principles calculations, which can benefit greatly the data transferability and compatibility between systems and also help us to reduce the NN fitting error.

G-NN potential
Our G-NN potential utilizes our recently proposed powertype structural descriptors (PTSDs) 37 to distinguish different structures on PES. The general forms of the PTSD are described in the following equations: 37 where rij is the internuclear distance between atom i and j, θ ijk is the angle centered at i atom with j and k being neighbors (i, j, and k are atom indices). The key ingredients in PTSD are the cut-off function f c that decays to zero beyond the rc [Eq. (3)], power-type radial function, trigonometric angular functions, and spherical harmonic function.
In PTSD, the S 1 and S 2 are two-body functions, the S 3 , S 4 , and S 5 are three-body functions, and the S 6 is a four-body function. The replacement of the Gaussian-type structural descriptor which The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcp was proposed by Behler and Parrinello by the PTSD has several advantages: (i) the computational cost in numerical calculations is reduced; (ii) the adjustable parameters are reduced to one (n) that simplifies the search for optimal parameters for two-body function; (iii) the power function when combining with the decaying cut-off function can create radial distributions with flexible peak and shape, which fulfills the similar purpose of Gaussian function (see Fig. 3); (iv) the introduction of different powers (n, m, and p) in the threebody function can couple conveniently different radial distributions.
(v) The introduction of a spherical function greatly improves the description of the angular environment of atoms. Our G-NN potential typically has two to three hidden layers and 50-120 neurons for each layer, which leads to the number of network parameters (weights and bias in NN) in the range of 10 4 -10 6 . The optimization of such a large parameter space is achieved via the quasi-Newton Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) method by simultaneously fitting the total energy, force, and stress from first principles that constitutes the cost function Jtot as follows: where ρ = 1-10 and τ = 0.05-0.1. The training procedure is controlled by a number of hyperparameters, such as the relative weight among energy, force, and stress [1:ρ:τ in Eq. (11)], network weight initialization method, type of neuron activation function, and training epochs. 40,41 A prototype of these hyperparameters frequently utilized is summarized in Table I. Due to the structural variety of the global PES, the typical accuracy for G-NN is 5-10 meV/atom for

A. PES exploration
The thermodynamics and kinetics of materials are determined by the underlying PES. SSW-NN simulation is such a convenient and efficient tool to establish the global PES for materials. In the past several years, we have mapped out the PESs for a number of systems, e.g., single element crystal (boron), 37 molecular crystal (ice), 42 and metal oxide (TiO 2 , ZnCr 2 O 4 ). 43,44 The PESs of TiO 2 and ZnCr 2 O 4 are projected onto a two-dimensional contour map, as shown in Fig. 4, where the structural order parameter (OP) and the relative energy are the xand y-axis, respectively. The density of the map represented by different colors indicates the density of states (DOS) on PES, showing the energy degeneracy of structural configurations at the same structural OP.
The OP is defined by the following equation: 45 where Y lm is the spherical harmonic function of degree l and order m; n is the normalized direction between all bonded atoms; i and j are atoms in the lattice, rij is the distance between atom i and j, and rc is set at 60% of the typical single bond length for i and j atoms. N bonds is the number of bonds. By choosing a suitable degree l, the order parameter can measure the short-and medium-range ordering of atoms in the lattice and thus distinguish important crystal structures and amorphous structures. For example, it is often straightforward to tell the coordination number from the OP value, e.g., six-coordinated Ti atom with OP2 = 0.3-0.5 and five-coordinated Ti atom with OP2 = 0.6-0.8 in TiO 2 PES. As shown in Fig. 4  has only a well-defined global minimum, a spineltype structure, and all the other structures are much higher in energy. It is not surprising that more than 10 different TiO 2 crystal structures have been synthesized in experiment to date [e.g., anatase, rutile, (R), and (H)], while only a spinel phase is known for ZnCr 2 O 4 . In addition to the global minimum, the global PES can also reveal unknown metastable structures, which may exhibit attractive physicochemical properties. In the case of TiO 2 , it is known in general chemistry that [TiO 6 ] octahedron is the common building block for TiO 2 crystals. However, a new class of unprecedented microporous TiO 2 crystalline phases, named TiO 2 (TB) in Fig. 4, have been identified, which features a special [TiO5] trigonal bipyramid building block, a large pore size (5-7 Å), and high thermal stability. 43 These microporous materials are predicted to be the candidate of anode materials for Li-ion and Na-ion batteries. 46

B. Origin of amorphous TiOxHy for hydrogen evolution reaction
Black titania (TiO 2 ) with strong absorption in the entire visiblelight spectrum has been found to exhibit the hydrogen evolution reaction (HER) activity several orders higher than the conventional TiO 2 material. [47][48][49][50][51][52] Synthesized by hydrogenating pristine TiO 2 , the as-synthesized black TiO 2 commonly exhibits a core-shell structure with the amorphous shell, a few nanometers thick, coated on the anatase crystals. 48,[53][54][55] The amorphous shell is believed to provide the catalytic active site and responsible for the enhanced HER activity. However, the structure determination for the amorphous shell is challenging to both experiment and theory, not even mentioning to understand the high HER activity.
To resolve the HER active sites on the amorphous TiOxHy shell, we have performed SSW-NN simulation to obtain the TiOxHy global dataset with 143 786 structures and built a robust and accurate TiOxHy G-NN potential. The structures in the database cover Using the TiOxHy G-NN potential, the thermodynamics phase diagram of TiO 2 bulk and surface in contact with H 2 at different temperatures and pressures can thus be determined quantitatively. We found that among common anatase surfaces, only the ridged anatase (112) surface can reconstruct significantly by surface H atoms and a local high H coverage, 0.69 ML, can gradually be built up during the surface amorphization. The amorphous surface exposes different Ti cations, 25% four-coordinated Ti 4c , 50% fivecoordinated Ti5c, and 25% six-coordinated Ti 6c atoms [ Fig. 5(a)]. Consistently, the Ti-O bond length has a wide distribution, from 1.8 to 2.2 Å, as compared to 1.9-2.1 Å on a perfect TiO 2 surface. This high H coverage not only renders the black color of the amorphous TiO 2 but also provides a low energy reaction channel for HER: a transient Ti-H hydride becomes likely to form on the exposed Ti atoms of the amorphous surface. The nascent TiH hydride can react facilely with the neighboring OH to produce H 2 , where the barrier is more than 1 eV lower than the traditional H coupling channel via two surface OH groups (barrier > 1.9 eV), as shown in Fig. 5(b). 44,56 We note that after the theoretical work is published, a recent experimental work reveals the presence of TiH hydride on amor- phous TiO 2 , where the characteristic H chemical shift δ value in 1 H nuclear magnetic resonance spectroscopy has a negative peak (−0.6 ppm) and this peak of δ grows with the increase in the hydrogenation time. 57

C. Syngas conversion on ZnCrO catalysts with varied Zn:Cr ratio
Zinc-chromium oxide (ZnCrO) catalyst is the first generation industry catalysts for syngas-to-methanol. Experiment shows that the catalytic activity is significantly affected by Zn:Cr ratios: 58-62 the best activity and selectivity are achieved at Zn:Cr = ∼1:1, while the Zn:Cr = 1:2 catalysts yield rather poor activity. 62 The atomic structure of ZnCrO catalysts with Zn:Cr = 1:1 remains, however, uncertain as no clear evidence of ZnO formation was found, although the X-ray diffraction (XRD) patterns exhibit peak broadening and small peak shifting with respect to the 1:2 phase, ZnCr 2 O 4 spinel phase. [62][63][64] The structural uncertainty in the ZnCrO catalyst is, in fact, quite typical for multicomponent oxide catalysts, where the atomic structures of the active sites are often unknown.
To resolve where and how syngas conversion occurs on the ZnCrO catalyst, we have performed the SSW-NN simulation to obtain the ZnCrO global dataset and to build a robust and accurate ZnCrO G-NN potential. Our simulations contain structures from 10 to 84 atoms per cell and cover different Zn:Cr:O ratios, i.e., ZnO, CrOx, ZnCrxOy, with different morphology forms, e.g., bulk, layers, and clusters. The final ZnCrO global dataset contains 38 285 structures. The ZnCrO G-NN potential contains 324 PTSDs for each element, i.e., 132 two-body PTSDs, 170 three-body PTSDs, and 22 fourbody PTSDs, and compatibly, the network involves three-hidden layers (324-80-60-60-1 net), equivalent to 103 743 network parameters in total. The final RMSEs of energy and force are 4.3 meV/atom and 0.128 eV/Å, respectively. 65 Based on the ZnCrO G-NN potential, the global PES of ZnCrO is explored and the thermodynamics phase diagram of Zn-Cr-O is constructed. It reveals a small stable composition island, i.e., Zn:Cr:O = 6:6:16-3:8:16, where the oxide alloy crystallizes into a spinel phase [ Fig. 6(a)]. At Zn:Cr = 1:1, a Zn 3 Cr 3 O 8 metastable crystal phase is present, also with the spinel crystal structure, but contains the highest concentration of unusual [ZnO 6 ] octahedra (O h ) in the bulk compared to the other spinel crystals. This subtle structural difference turns out to be critical to affect the syngas conversion   Fig. 6(b)]. The microkinetics simulation based on DFT reaction energetics (methanol yield: ∼670 g kgcat −1 h −1 on Zn 3 Cr 3 O 8 ; methane yield: ∼0.08 g kgcat −1 h −1 on ZnCr 2 O 4 ) further rationalizes the sharp difference in activity and selectivity observed in experiment: the ZnCrO catalysts with Zn:Cr = 1:1 shows high activity and high methanol selectivity (methanol yield 80-600 g * kgcat −1 * h −1 and selectivity 80% at 573 K), but the ZnCrO catalysts with Zn:Cr = 1:2 shows low activity and low methanol selectivity [<5 g * kgcat −1 * h −1 with 14% selectivity to methanol at 573 K and ∼82% to hydrocarbon (91% CH 4 ) at 673 K]. 62,66 The findings on ZnCrO catalysts not only strengthen our understandings, from bulk to surfaces and to active sites, on the high-temperature syngas conversion on oxides but also suggest other complex oxide systems hotly debated in experiment, e.g., ZnZrO, CuZnO, ZnMnO, and ZnFeO, [66][67][68] should now be possible to explore with the advent of the SSW-NN method.

IV. REMARKS AND PROSPECTS
This perspective overviews the SSW-NN method and its recent applications in heterogeneous catalysis. While SSW-NN holds great promise in solving some challenging tasks in heterogeneous catalysis, ranging from catalyst active site structure to catalytic activity as demonstrated in the above examples, one must bear in mind that a successful application of the SSW-NN method relies much on the sampled dataset to the target problem. The SSW sampling becomes frustrated for some molecular reactions and for systems with too large degrees of freedom, as elaborated below.
(i) Reaction sampling for molecular reactions have a high computational cost even with the SSW method. It is especially problematic for reactions with a sharp TS region (i.e., characterized by a large negative curvature), where the exact TS position is difficult to capture by random sampling. As introduced above, the current NN potential generally focuses on materials such as TiO 2 and ZnCrO, and a giant step forward would couple molecular reactions with the solid materials to establish reactive NN potentials for heterogeneous catalysis. (ii) The iterative self-learning in building the global NN potential may become highly computational when the system size is large and the bonding complexity increases. In particular, the large configurational space of bonding environment in multielement systems is a key problem in constructing G-NN. For example, there are many possible bonding patterns for carbon with other elements (C, H, O, and N) in organic chemistry, which leads to great difficulty in building a general-purpose reactive NN potential for organic chemistry.
To enhance the reaction sampling and to reduce the system degree of freedom effectively, we expect that the combination of other techniques with SSW-NN, such as pattern recognition artificial intelligence methods, coarse graining, and rigid body methods could be the future direction to power up SSW-NN for complex catalysis problems. Our ongoing implementation of LASP software 27 could help us to fulfill this goal by integrating different techniques in one platform, which aims to accelerate the material simulation in the future.