Atom-centered neural network (ANN) potentials have shown promise in computational simulations and are recognized as both efficient and sufficiently accurate to describe systems involving bond formation and breaking. A key step in the development of ANN potentials is to represent atomic coordinates as suitable inputs for a neural network, commonly described as fingerprints. The accuracy and efficiency of the ANN potentials depend strongly on the selection of these fingerprints. Here, we propose an optimization strategy of atomic fingerprints to improve the performance of ANN potentials. Specifically, a set of fingerprints is optimized to fit a set of pre-selected template functions in the *f*^{*}*g* space, where *f* and *g* are the fingerprint and the pair distribution function for each type of interatomic interaction (e.g., a pair or 3-body). With such an optimization strategy, we have developed an ANN potential for the Pd_{13}H_{2} nanoparticle system that exhibits a significant improvement to the one based upon standard template functions. We further demonstrate that the ANN potential can be used with the adaptive kinetic Monte Carlo method, which has strict requirements for the smoothness of the potential. The algorithm proposed here facilitates the development of better ANN potentials, which can broaden their application in computational simulations.

## I. INTRODUCTION

Atomic simulation approaches based on quantum mechanics, empirical interatomic potentials, or classical force fields have demonstrated a remarkable power in gaining chemical and physical insights and providing guidance for the design of novel materials. Quantum mechanical methods, primarily based on density functional theory (DFT), offer sufficient accuracy for the simulation of chemical process and prediction of material properties. However, the computational cost of such methods remains high, which limits their application in large systems containing hundreds of atoms and where extensive sampling is required. Relatively inexpensive methods based on classical force fields are capable of modeling systems with thousands or even millions of atoms. However, classical force fields are carefully parameterized to a specific chemical environment and suffer from a lack of transferability. In recent decades, machine learning (ML) algorithms have been introduced to describe interatomic interactions by fitting the potential energy surface with reference data obtained from quantum mechanical calculations.^{1–7} ML potentials have been demonstrated to have DFT-level accuracy and orders of magnitude lower computational cost and show applicability to chemical and material systems involving bond formation and breaking.^{8–10}

Until now, various ML techniques have been introduced for atomic simulation, including artificial neural networks (NNs),^{11–16} kernel-ridge regression,^{1,17–19} and Gaussian process regression.^{20–22} A key step in these techniques is to represent the atomic geometries with descriptors that are invariant to the operation of rotation and translation so as to preserve the intrinsic symmetries of the system. Such descriptors, commonly described as atomic fingerprints, are designed to reflect the unique features of each geometric structure. Several different fingerprint algorithms have been developed, i.e., atom-centered symmetry functions (ACSFs),^{11} the smooth overlap of atomic positions (SOAP),^{23} bispectrum fingerprints,^{20} and Zernike fingerprints.^{24} The parameters in these fingerprints are boundless and able to form infinite dimensional vectors. Theoretically, the more complete the fingerprints, the higher accuracy can be reached, but at the expense of a higher computational cost. Therefore, an optimal selection of fingerprint parameters is essential to the performance of ML techniques.

In this work, we present a simple but effective algorithm for selecting fingerprint functions without relying on pre-training of the ML potential. The fingerprints we use here are based upon the ACSFs proposed by Behler and Parrinello (denoted as BP fingerprints).^{11} The BP fingerprints have been widely used for the construction of atom-centered neural network (ANN) potentials that have a broad applicability across the periodic table and are able to provide accuracy comparable with the methods used to construct the reference data. Here, we first evaluate the performance of the ANN technique by fitting the potential energy surface of a 7-atom Lennard-Jones (LJ_{7}) particle and determine how the accuracy of the model varies with the set of BP fingerprints. Inspired by the LJ_{7} particle case, we developed an algorithm to optimize the BP fingerprints to improve the performance of the ANN model.

To test our algorithm, we have fit forces and energies from DFT calculations of a Pd_{13}H_{2} nanoparticle (NP) with the ANN technique. The PdH system is of interest in various areas, including hydrogen storage,^{25,26} hydrogen separation membranes,^{27} hydrocarbon reforming,^{28–31} etc. Many scientific issues, i.e., the reactivity and dynamics of PdH nanoparticles, are important for these relevant applications. Here, we choose to model the dynamics of a PdH NP over long time scales using the adaptive kinetic Monte Carlo (AKMC) method.^{32,33} In our tests, we have not been able to find a force field with even qualitative agreement with forces and energies from DFT for these PdH NPs. Thus, here, we developed an ANN model for the Pd_{13}H_{2} NP system and used AKMC to understand the time scales of the dynamics of the PdH NPs.

## II. COMPUTATIONAL DETAILS

### A. Atom-centered neural-network model

We used the TensorFlow version of the ANN algorithm implemented in the Atomistic Machine-Learning Package (AMP) developed by Khorshidi and Peterson for ML potential fitting.^{34} As shown in Scheme 1, the ANN algorithm first interprets the structural information into unique fingerprints by representing the atomic configurations as single-atom chemical environments with generalized numerical representations. Here, we applied modified BP fingerprints of radial ($fiI$) and angular ($fiII$) forms, with centers at *R*_{s} and *θ*_{s}, respectively,

where

We set the cutoff radius (*R*_{c}) as 3.0 Å and 6.0 Å for the LJ_{7} and Pd_{13}H_{2} particle system, respectively. The parameters *η* and *ζ* are used to control the radial and angular fingerprint widths, respectively. For example, a larger *η* (*ζ*) value yields a narrower distribution of the radial (angular) symmetry function. A set of symmetry functions with different *η* and *ζ* values is required to describe the chemical environment of each atomic center. The numerical representations generated were weighted through an (*n*, *m*) NN model. Here, the length of the tuple (*n*, *m*) represents the number of hidden layers in the NN model with *n* and *m* corresponding to the number of neurons in each hidden layer. In this study, a (5) and a (50) NN model was used for the LJ_{7} and Pd_{13}H_{2} particle system, respectively. This corresponds to a NN model with one single hidden layer consisting of 5 or 50 neurons. Each neuron uses an activation function (the sigmoid function is selected here) to transform its input values to the output value. The output atomic energies and forces are evaluated with a loss function as defined in Scheme 1. The weights connecting each layer were initialized from a standard normal distribution and then optimized to minimize the loss function with the Limited-Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm. Finally, the resulting ANN model was used as a LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) force field through the OpenKIM module (refer to Part II of the supplementary material for more details).^{35}

### B. Data collection

For the Pd_{13}H_{2} particle, high-temperature molecular dynamics (high-T MD) simulations followed by geometry optimization (GEO-OPT) and saddle-point searching based on the dimer method were performed to obtain reference data.^{36–38} In the high-T MD/GEO-OPT simulations, we first ran three independent MD simulations with an NVT ensemble and the temperature fixed at 1000 K to produce high-energy structures. Each MD simulation ran for 20 ps with a time step of 1 fs. In this step, we used Gaussian and plane-wave basis functions with the PBE functional as implemented in the CP2K package.^{39} A small single-ζ Gaussian basis set and the Goedecker–Teter–Hutter norm-conserved pseudopotentials were used to treat valence and the core electrons, respectively.^{40,41} The energy cutoffs for the finest grid level and Gaussian wavefunctions were set as 250 Ry and 30 Ry, respectively. Next, the fast inertial relaxation engine (FIRE) optimizer, as coded in the atomic simulation environment (ASE), was used to optimize the high-energy structures; 40 structures were collected along each optimization.^{42} Additional structures were also included to describe the chemical environments that were not found in either the dimer saddle search or high-T MD/GEO-OPT simulations. Specifically, geometries for the dissociation of the H_{2} molecule and H atoms embedding inside the Pd_{13} NP were also included. The data collected from each method was randomly split into training and testing sets in the ratio of 4:1 (Table I).

. | . | . | . | Embedded . |
---|---|---|---|---|

Categories . | MD/GEO-OPT . | Dimer . | H_{2} dissociation
. | H atom . |

Training | 5414 | 6312 | 1213 | 6863 |

Testing | 1354 | 1578 | 304 | 1716 |

. | . | . | . | Embedded . |
---|---|---|---|---|

Categories . | MD/GEO-OPT . | Dimer . | H_{2} dissociation
. | H atom . |

Training | 5414 | 6312 | 1213 | 6863 |

Testing | 1354 | 1578 | 304 | 1716 |

The reference energies and forces for the Pd_{13}H_{2} NPs were evaluated with DFT as implemented in the Vienna *ab initio* simulation package (VASP 5.4.4). Projected augmented wave based pseudopotentials and a plane-wave basis set with an energy cutoff of 300 eV were employed to treat the core and the valence electrons, respectively.^{43–45} The generalized gradient approximation with the PBE functional was selected to calculate the exchange–correlation interaction.^{46} A 15 × 15 × 15 Å^{3} cell was used to avoid interaction from neighboring particles. Note that the force field developed here was fit to the hydride system and is not expected to be transferrable even to similar systems such as Pd_{13} with molecular H_{2}. In fact, tests show that the ANN is qualitatively different from DFT for recombinative H_{2} desorption from Pd_{13} (gray areas shown in Fig. S1). This discrepancy is due to the lack of training data in those regions. The objective here, however, is not to develop a general PdH potential; rather, it is to develop a strategy to optimize fingerprints and demonstrate the feasibility of combining an ANN model with the AKMC algorithm.

## III. RESULTS AND DISCUSSION

### A. PES fitting of the LJ_{7} particle

First, we selected a simple system, the 7-atom LJ particle, to test the sensitivity of the ANN model with respect to variations of the fingerprint functions. We performed a basin-hopping calculation to locate all local minima of the LJ_{7} particle and then constructed a disconnectivity graph by finding transition states that connect the local minima using the Pele package.^{47} Reduced units are used for the Lennard-Jones potential, which is equivalent to setting the pair equilibrium well depth (*ε*) and separation (*σ*) to 1. The reason to focus on a disconnectivity graph as a measure of accuracy of the PES is that we are interested in modeling dynamics where the energy of minima and saddles are of particular importance. A total of 2609 structures were selected from the basin-hopping simulation, with 133 being selected from those with a minimum distance between atoms *R*_{min} < 1.0, 165 from those with *R*_{min} > 2.5, and 2307 from the rest plus the four local minima. This set of structures was used to fit the LJ_{7} PES. The performance of the ANN model was evaluated by comparing the calculated ANN disconnectivity graph with that obtained from the reference LJ potential.

Figure 1(a) shows the disconnectivity graph of the LJ_{7} particle obtained from the LJ potential. The LJ_{7} particle has four local minima (black dots) that are connected by three unique transition states (the vertices that connect each pair of local minima). The red counterpart shows the disconnectivity graph obtained from the ANN model; note that the red and black plots are so similar that they appear superimposed on one another. One hidden layer with five neurons was used in this ANN model. Five radial symmetry functions centered at zero with η = 0.05, 4, 20, 40, 80, and 120 were used as fingerprints. Note that this set of fingerprints was determined after numerous trials. We first followed the recommended fingerprints in Ref. 34 and then varied the interval between *η* values to find fingerprints that reproduced the disconnectivity graph of the LJ_{7} particle. The near perfect match between the black and red lines indicates a near perfect agreement of the disconnectivity graph obtained from the ANN model with that from the LJ potential and that we have used a suitable set of fingerprints.

To understand the importance of each fingerprint function, we retrained the ANN model with one fingerprint removed and calculated the difference of the disconnectivity graph between the LJ potential and the ANN model as the root mean square error

where $EiANN$ and $EiLJ$ are the energies of the local minima and saddle points in the disconnectivity graph from the ANN model and the LJ potential, respectively, and *N* is the total number of local minima and saddles. If the ANN finds more or fewer critical points than the LJ potential, the energies of the additional stationary points found in one model are evaluated in the other without relaxation, to provide a measure of the difference between the two landscapes in the RMSE.

The variation of the RMSE with one fingerprint removed is presented in Fig. 1(b) (top panel). Different sensitivities are observed for different fingerprints. For example, the removal of the *η* = 20 fingerprint, corresponding to a mean distance of 1.1, produces the largest error, suggesting that the ANN model is most sensitive to this fingerprint. In contrast, the ANN model with the *η* = 0.05 fingerprint removed, corresponding to a mean distance of 1.5, performs as accurately as with all fingerprints applied. The mean distance mentioned here corresponds to the center of the *f*(*r*)*g*(*r*) plot in the bottom panel of Fig. 1(b). This analysis indicates that the *η* = 0.05 fingerprint is redundant for the ANN model. The removal of other fingerprints had a similar influence on the accuracy of the ANN model, indicating that only the *η* = 20 fingerprint is important for accuracy, which is consistent with modeling a simple pair potential in which the local energy of an atom is based upon the number of neighbors.

To understand the different sensitivities of the fingerprints, we computed the product of each fingerprint, *f*(*r*), with the pair distribution function, *g*(*r*), of the training data [denoted as *f*(*r*)*g*(*r*)]. To see this more clearly, an enlarged plot for *η* = 4 and 20 is shown in Fig. S2. This product describes how much a specific *r* value is weighted in the training data by this fingerprint. Such a measure helps describe which fingerprint distances are present in the training data. As shown in the bottom panel in Fig. 1(b), the fingerprints selected here emphasize the description of various radial distributions. For example, only one sharp peak centered at the equilibrium distance (*r* = 1.1) is found for the *η* = 20 fingerprint. This indicates that this fingerprint describes the two-body interaction around the equilibrium distance of the system and explains why the ANN model is most sensitive to this fingerprint. In contrast to the *η* = 20 fingerprint, the plot of the *f*(*r*)*g*(*r*) with *η* = 0.05 has almost the same shape as that of the *g*(*r*). Besides the peak at *R*_{ij} = 1.1, contributions to this fingerprint from other regions beyond 1.3 are also considerable due to a significant population of those distances in the training data. Nevertheless, their contributions to the LJ potential are negligible compared to those from the region near the equilibrium distance. In other words, the *η* = 0.05 fingerprint is unable, on its own, to correctly describe the interaction near the equilibrium distance due to the inclusion of contribution from other regions. Therefore, the removal of the *η* = 20 fingerprint leads to a large deviation from the real LJ potential whereas a small error is induced upon removal of the *η* = 0.05 fingerprint. These results suggest that to reach a reasonable accuracy for the ANN model, each fingerprint should probe a specific region of the atom-centered chemical environment. The product of *f*(*r*) and *g*(*r*) is potentially a good criterion for the determination of fingerprints that are distinguishable from each other in this regard.

### B. Fingerprints optimization

A principle for designing an ANN algorithm is to have fingerprints distinguishable from each other such that each one probes a specific region of the atom-centered chemical environment. Following this rule, Smith *et al* developed extensible neural network potentials for multiple common organic molecules by evenly distributing the value of *R*_{s} and *θ*_{s} with fixed *η* and *ζ* values.^{16} An even distribution of the fingerprint center guarantees inclusion of all possible chemical environments to be explored. However, a biased population of a specific type of chemical environments in the training data will overweight or underweight the corresponding fingerprint, leading to a poor performance of the ANN model. Following the sensitivity test of the LJ_{7} particle, the product of each fingerprint with the pair distribution function *g*(*r*) can be used to determine how a specific fingerprint is biased in training data. In other words, we need to find a set of fingerprints that are distinguishable from each other upon taking the product with the pair distribution function. This means finding distinguishable fingerprints in the *f*(*r*)*g*(*r*) (simplified as *f*^{*}*g*) space instead of the *f*(*r*) space.

Here, we propose an optimization algorithm to find distinguishable fingerprints in the *f*^{*}*g* space and apply the algorithm to locate fingerprints for the Pd_{13}H_{2} nanoparticle system. Figure 3 presents the pair distribution function and radial fingerprints near the equilibrium distance for the Pd–H pair before and after optimization. We first defined a set of template fingerprints $fiTp(r)i=1,Nf$ that describe separate regions in the pair distribution space to be described. Then, an objective function (** τ**) defined by the following equation is minimized to find fingerprints $fiOptri=1,Nf$ that reproduce the template in the

*f*

^{*}

*g*space,

where *α* is the weight parameter set to ensure that the error in the tail region is properly weighted. The tail region is defined as that region where the value of the template fingerprint is smaller than 1/100 of its maximum value. The *α* value can be determined as the minimum value that eliminates the error in the tail region. An example is shown in Fig. 2, where the tail corresponds to the region with *r* > 2.0 Å. A small *α* leads to mismatch between the optimized and the temperate fingerprint in this tail region, i.e., the black line (the optimized fingerprint with *α* = 1) as compared to the gray-circled line (the template fingerprint). Increasing the value of *α* reduces the mismatch in the tail region and leads to a right shift of the fingerprint center (the *R*_{s} value) as shown in the right panel of Fig. 2. The *R*_{s} value reaches a maximum value with *α* = 100—that is, the minimum *α* value to be used to eliminate the mismatch in the tail region. As such, we set *α* as follows:

Three fingerprints with *R*_{s} = 1.6, 1.75, and 1.9 Å and *η* = 1500, as shown by dashed lines in the middle panel of Fig. 3, are selected as templates to describe the interactions centered at the respective *R*_{s} value in the *f* space. The *g*(*r*) of the H–Pd pair distribution is most populated around the first-nearest distance (1.77 Å) as shown in the top panel of Fig. 3. The strong population bias at the first-neighbor distance leads to a reduction of the fingerprints in the *f*^{*}*g* space toward the peak position of the *g*(*r*). The center of these populated fingerprints (dashed lines in the bottom panel of Fig. 3) ranges from 1.72 Å to 1.82 Å, much narrower than our target (1.6–1.9 Å). This indicates that the combination of template fingerprints with our training data tends to emphasize the interaction in a region close to the first-nearest neighbor distance (±0.05 Å).

By optimizing the objective function (** τ**), we shifted the centers of the populated fingerprints to our target range. The objective function is continuous and smooth with respect to

*η*and

*R*

_{s}. An example contour map of

**for the template fingerprint centered at 1.6 Å is given in Fig. S3. There exists only one critical point, which is the global minimum. The global minimum of**

*τ***can be located in the fast inertial relaxation engine (FIRE) as shown in the inset image of the bottom panel in Fig. 3. The fingerprints obtained after optimization are presented as solid lines in the middle panel of Fig. 3. Compared to the template fingerprints, the two outermost fingerprints exhibit a similar shape but are shifted away from the first-nearest distance in order to lower the weight in that region. The coverage of the middle fingerprint is expanded with a slightly shifted center. As a result, the optimized fingerprints cover radial regions as defined by the template in the**

*τ**f*

^{*}

*g*space (solid lines in the bottom panel of Fig. 3). To fully describe the Pd–H interaction, another six fingerprints are selected and optimized following the same procedure (Fig. S4). These fingerprints are expected to improve the description of interactions centered at the target

*R*

_{s}values for the training data.

### C. An improved ANN model with optimized fingerprints for the Pd_{13}H_{2} system

Next, we located fingerprints for other interactions following the proposed optimization procedure, including the H–H, and Pd–Pd radial, and the Pd–Pd–Pd angular terms. Only the Pd–Pd–Pd angular term is included given that the angular term with hydrogen has little contribution to the system. We followed the exact same procedure used in the radial term to optimize the angular term. The template and optimized fingerprints are presented in Figs. S5–S7 and the parameters are included in the attached potential files. To evaluate the performance of the optimized fingerprints, we trained ANN models to fit the DFT-based PES of the Pd_{13}H_{2} NP with both the template and optimized fingerprints. Ten independent models are trained for each set of fingerprints and the one with the best performance discussed here. In the models from both sets of fingerprints, the training and the test sets show comparable accuracy without obvious overfitting (Figs. S8 and S9).

Figure 4(a) shows the predicted energy (left panel) and forces (right panel) from the template (blue) and the optimized (red) models against the reference DFT values. In both panels, the optimized model exhibits a tighter distribution along the *y* = *x* line as compared to the template model. This indicates that the ANN model with optimized fingerprints reaches a higher accuracy compared to those with the template fingerprints. Specifically, the template model shows a lower performance on force prediction as indicated by the sparse distribution of blue dots in the right panel of Fig. 4(a). These results demonstrate that the ANN model is improved for the training and test data upon fingerprint optimization.

To further validate our models, we artificially pulled one Pd atom and one H atom away from the particle and evaluated the energy of structures along this trajectory as indicated in the inset in Fig. 4(b). The structures used in this validation were not included in the training data. The resulting energy profiles from the template (blue) and optimized (red) models are compared to the reference DFT values in Fig. 4(b). The energy profile from the optimized model shows better agreement with the reference DFT results than that from the template model. Specifically, the template model underestimates the energy at *r* < 2.7 Å but overestimates at *r* > 2.7 Å in the Pd–Pd case and predicts a flat basin from 2.5 Å to 2.7 Å, inconsistent with the reference DFT data. The optimized model largely minimizes this discrepancy. In the H–Pd case, the optimized model is able to reproduce DFT energies of structures with *r* < 2.0 Å and increases the prediction accuracy for structures beyond that region. More importantly, the optimized model yields a smooth energy profile along the trajectory, whereas a turning point is observed in the template model. These validation results further indicate that our fingerprint optimization method improves the global performance of the ANN model.

### D. AKMC simulations with the ANN model

The AKMC method is a coarse-grain dynamic simulation method with time evolution of the system accumulating on the time scale of the state-to-state transitions.^{32,33} In our previous study, we showed that the AKMC method is able to mitigate the time scale problem in molecular dynamics simulations and extend the time scale to that of experiments.^{48} Nevertheless, a significant challenge for the AMKC method is the availability of reliable and inexpensive force fields that are accurate and smooth enough to facilitate the location of transition states.

Here, we performed AKMC simulations to compute the dynamical evolution of a Pd_{13}H_{2} NP at 300 K with energies and forces evaluated by our optimized ANN model. In the AKMC simulation, a minimum distance (*R*_{min}) between each atomic type was set to avoid the collapse of the ANN model beyond the range of the training data. Any distance below *R*_{min} was set to *R*_{min}. The AKMC method as implemented in the EON package was used to perform the simulation.^{49} The dimer method was used to find all possible low-energy transition events, and the kinetic Monte Carlo algorithm was used to determine which event the system proceeds with. At each new state, such dimer/KMC cycles were repeated to simulate the state-to-state evolution of the system. More details about the AKMC method can be found in Ref. 48.

The blue line in Fig. 5(a) shows the time evolution of the total energy of the system from the AKMC simulation based on the ANN model. Roughly, five energy plateaus associated with five types of structures, as grouped in the disconnectivity graph [Fig. 5(b)], are observed in the simulation. Each drop in the energy evolution plot corresponds to a structural transition, except for the transition from 4 to 5, which leads to a slight increase in energy. For comparison, we also calculated the energy of the local-minima structures extracted from the AKMC simulation with DFT [red line in Fig. 5(a)]. The time evolution of the total energy for the AKMC simulation based on the ANN model is consistent with the DFT reference. Using DFT, we further relaxed the most stable structures, corresponding to the green dots in the disconnectivity graph. The structures obtained with the corresponding energy (*E*_{diff}) and structure (*S*_{diff}) difference are shown in Fig. 5(c). The gold (light blue) and green (white) colors represent Pd and H atoms in the DFT-relaxed (ANN-relaxed) structures, respectively. The structures obtained from both methods are in agreement with *S*_{diff} values smaller than 0.05 Å/atom. The energy differences range from 0.05 eV to 0.20 eV. These results demonstrate the reliability and accuracy of our ANN.

As discussed above, the NP undergoes a structural transition accompanied by diffusion of H atoms on the NP surface. The overall energy barrier for the structural transition from type *k* to *k* + 1 can be extracted from the disconnectivity graph, taken as the energy difference between the minimum-energy structure (green dots in Fig. 5(b)] of type *k* and the vertex that connects to a *k* + 1 state. For example, the 4 → 5 transition has an overall barrier ∼0.33 eV as indicated by the gray dotted lines in Fig. 5(b). In this way, we are able to characterize overall energy barriers for the *k* → *k* + 1 transitions ranging from 0.15 eV to 0.33 eV. Each structural type entails a specific Pd skeletal structure with H atoms that can almost freely diffuse on the surface of the Pd structure with very low barriers (<0.07 eV). These results suggest that the Pd_{13}H_{2} particle stays in a dynamic state with facile movement of H. Throughout the AKMC simulation, we did not observe embedding of H atoms into the nanoparticle. Likely, a higher barrier is required for such embedding events, whereas the diffusion of surface H atoms with low barriers limits the timescale of the AKMC simulation.

## IV. CONCLUSIONS

In this study, we first tested the sensitivity of the ANN model with respect to the ACSF by training an ANN model for the LJ_{7} particle and found that the product of the fingerprint functions with the corresponding pair distribution function are responsible for the performance variation of the ANN model. Inspired by this sensitivity test, we defined an objective function that can be minimized to find optimal fingerprints, thus improving the performance of the ANN model. In this optimization process, we first define a set of template fingerprints based on the radial and angular regions to be explored and then minimize the objective function to locate fingerprints that can reproduce the template fingerprints in the *f*^{*}*g* space. Accordingly, we successfully developed an ANN model for the Pd_{13}H_{2} NP system that is superior to the model trained with the template fingerprints. Furthermore, we applied the developed ANN model to AKMC simulations of dynamic behavior of the Pd_{13}H_{2} NP, proving the reliability and accuracy of the model.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional details about the structure of the optimized fingerprints and their efficacy.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

Funding for this project was provided by the National Science Foundation (Grant No. CHE-1764230) and Exxon-Mobil (Grant No. UTA17-000317 EM10480.13). Computational resources were provided by the Texas Advanced Computing Center and the XSEDE program (Grant No. CHE190010). G.H. acknowledges the ongoing support from the Welch Foundation (Grant No. F-1841).

## REFERENCES

_{147}nanocluster

_{2}O

_{3}catalyst