Atomic cluster expansion for Pt-Rh catalysts: From ab initio to the simulation of nanoclusters in few steps

Insight into structural and thermodynamic properties of nanoparticles is crucial for designing optimal catalysts with enhanced activity and stability. We present a semi-automated workflow for parameterizing the atomic cluster expansion (ACE) from ab initio data. The main steps of the workflow are the generation of training data from accurate electronic structure calculations, an efficient fitting procedure supported by active learning and uncertainty indication, and a thorough validation. We apply the workflow to the simulation of binary Pt-Rh nanoparticles that are important for catalytic applications. We demonstrate that the Pt-Rh ACE is able to reproduce accurately a broad range of fundamental properties of the elemental metals as well as their compounds while retaining an outstanding computational efficiency. This enables a direct comparison of simulations to high resolution experiments.


Introduction
In the field of catalysis, noble metals play a key role.Among them, platinum (Pt) and rhodium (Rh) are used extensively as heterogeneous catalysts due to their superior activity and stability [1].For example, Pt-Rh based materials are widely used as effective vehicle exhaust catalysts to reduce emissions of harmful pollutants CO and NO x into less harmful CO 2 and N 2 [2,3].Furthermore, they have also been employed in fuel cells [4,5] to enhance the oxygen reduction and hydrogen oxidation reactions.
These catalytic materials are typically applied in the form of nanoparticles or nanoclusters to maximize the active surface area where the chemical reactions take place [6].An important aspect of the nanoparticles is their microstructure that can be tailored to optimize functional properties as well as stability.In multicomponent systems, the goal is often not to achieve a random distribution of chemical elements as in a solid solution, but to design nanoparticles with specific morphologies [7,8].For instance, despite the fact that Pt and Rh are miscible elements with face-centered cubic (fcc) structures their nanoparticles have usually been prepared in the form of concentric composites, typically with a spherical Rh core surrounded by a Pt-rich shell layer.Such core-shell morphologies can be prepared using a wet chemical approach, where the core is synthesized first and the surface is coated afterwards [7,8].They can exhibit an enhanced catalytic activity and structural stability compared to traditional alloyed catalysts [7] as a result of their non-equilibrium crystal structures [9] or shapes.However, the non-equilibrium nanoparticles may not be stable under operating conditions and undergo transformations that lead to deterioration of the catalytic activity.
In Fig. 1, an example of a Pt-Rh core-shell nanoparticle with a distorted octahedral shape is shown [10].The high angle annular dark field (HAADF) micrograph and energy dispersive X-ray spectroscopy (EDS) map acquired in a scanning transmission electron microscope (STEM) reveal a varying Pt shell thickness of 3-5 monolayers.STEM is widely used for characterizing catalyst nanoparticles, since it can provide information on the atomic arrangement and composition [11].However, investigations of dynamic processes at elevated temperatures requires to carry out challenging and time consuming in-situ experiments that are often limited by beam-induced artefacts [12].Therefore, it is desirable to complement the experimental characterization with atomistic modeling to obtain a deeper understanding of the thermodynamic and kinetic phenomena governing the structure of nanoparticles.
Atomistic simulations can provide valuable information about energetics of particular structural features, such as formation energies of point defects or various surface facets, but also directly follow the dynamical evolution of the entire system using molecular dynamics (MD) or explore the thermodynamic stability using Monte Carlo (MC) methods.Even though it is nowadays possible to study systems of millions of atoms, approaching realistic dimensions of nanoparticles, the essential prerequisite for reliable quantitative predictions is an accurate and efficient description of interatomic interactions.However, there exist, to our best knowledge, no interatomic potentials for the Pt-Rh binary system that could be applied in large-scale MD or MC studies of thermodynamic and kinetic phenomena.
In this work, we developed an atomic cluster expansion (ACE) parametrization for Pt-Rh using a semi-automatic workflow consisting of the generation of a training dataset based on accurate density functional theory (DFT) [13,14] calculations, efficient fitting and validation procedures [15], and uncertainty indication and active learning (AL) algorithms [16] that can be employed to improve the model if necessary.ACE has successfully been applied to model metallic, covalent and ionic materials and was shown to offer superior accuracy and computational efficiency [15][16][17].
The ACE methodology combines the advantages of machine-learned (ML) methods and physically based models of interatomic interactions.We summarize only the essentials of ACE and refer to original publications for more details [15,[18][19][20][21][22].
One of the key features of ACE is a complete and hierarchical set of basis functions B iv that span the space of local atomic environments.This enables to expand an atomic property ϕ (p) i such as the energy of atom i as with expansion coefficients c v where v is composed of several indices.The basis functions fulfill fundamental translation, rotation, inversion and permutation (TRIP) invariances for the representation of scalar variables, or equivariances for the expansion of vectorial or tensorial quantities.
In the simplest case, the energy is evaluated linearly as but most ACE parametrizations to date, including the presented Pt-Rh model, use two atomic properties with a Finnis-Sinclair square-root embedding, i .
We demonstrate that Pt-Rh ACE provides not only an accurate description of fundamental properties for both elemental metals as well as their compounds, but that it is well suited for large-scale atomistic simulations of nanoclusters.

Training of ACE
The parameterization of ACE models requires a consistent set of reference DFT data.For the Pt-Rh system, we generated about 30 000 atomic configurations and evaluated their energies and forces using the Vienna ab-initio simulation package (VASP) [23][24][25] and the PBE functional for approximating exchange and correlation [26].The considered structures span a wide portion of the configuration space for both elemental metals and their binary compounds.The configurations include perfect bulk structures, common defects, e.g., vacancies, surfaces and grain boundaries, and small random clusters.For all these configurations, we varied their volumes and randomly displaced atoms away from their equilibrium positions to sample vibrational degrees of freedom.A summary of energies per atom for all structures as a function of the nearest neighbour distance in each structure is shown in Fig. 2. The training data corresponds to 29800 structures with 437157 atoms, spanning nearest neighbour distances from 1.8 up to 6.3 Å and sampling the energy space from about −6.3 eV/atom to +2 eV/atom.The ACE parameterization was carried out using the pacemaker package [21].In the fitting process, a hierarchical basis extension scheme with a sequential addition of ACE basis functions was adopted.For the radial basis functions, spherical Bessel functions were used.Structures within an energy window of 1 eV/atom above the convex hull were assigned significantly larger weights in the loss function optimization, while structures with higher energy were given smaller weights.
After initial training, we employed AL to ensure a reliable description of Pt-Rh surfaces as these are crucial for simulations of core-shell nanoclusters (see below).This was achieved by computing the extrapolation grade γ based on the D-optimality criterion [16] for a number of (100), ( 110) and (111) fcc surface slabs containing up to 56 atoms with two compositions (Pt 0.3 Rh 0.7 and Pt 0.7 Rh 0.3 ).The slab configurations were obtained using a hybrid MD-MC method in the N V T ensemble at three different temperatures (300, 500 and 800 K).Structures with large values of the extrapolation grade (γ > 5) were selected using the MaxVol algorithm.These structures were then computed with DFT, added to the training set, and the ACE potential was retrained.See Section 4.4 for more details.
In Fig. 3, the energies and forces predicted by ACE are compared with the reference DFT data for the training set.The ACE model reproduces the DFT energies and forces with good accuracy.Obtained root mean square errors (RMSE) for the energy and force components are 13 meV/atom and 57 meV/ Å for the training set and 36 meV/atom and 140 meV/ Å for the test set, respectively.The relativley large errors in the energies are due to the very low weights that were assigned to high energy structures with atoms at close distance as these structures have no relevance for the simulation of the nanoparticles.In fact, for the lower energy structures within the 1 eV window above the convex hull that are most significant in the simulations, the RMSE are only 6 meV/atom and 30 meV/ Å for the training set and 6 meV/atom and 26 meV/ Å for the test set, respectively.

Validation of the Pt-Rh ACE
For an initial validation of the parameterized potential, various fundamental properties of the elemental metals were examined and compared to DFT and experimental data in Table 1.For both Pt and Rh, the equilibrium lattice parameters predicted by ACE are in excellent agreement with DFT, while being slightly larger than the experimental values.The elastic constants and bulk moduli also agree well with those of DFT, with somewhat larger discrepancies observed for C 11 and C 44 for Pt.The differences however remain comparable to those between DFT and experiment.The larger lattice constants and mostly smaller elastic moduli predicted by DFT and ACE in comparison with experiment are likely caused by the PBE-GGA exchange-correlation functional that leads to underbinding for nonmagnetic fcc metals [27].
The formation and migration energies of the mono-vacancy are in close agreement with the reference DFT values and comparable with the available experimental measurements for Rh.For Pt, the large discrepancy between the theoretical and experimental values is due to deficiencies in the description of the exchange-correlation contributions at metallic surfaces [28].A better agreement can be obtained after correcting for this intrinsic surface error, but the correction term depends on the electronic structure of the metal [28] and cannot be employed in classical atomistic models.
We also examined energies of surfaces with low-index orientations.As illustrated in Table 1, ACE captures correctly the ordering of Pt surface energies in comparison to DFT, as γ (100) > γ (110) > γ (111) .However, DFT strongly underestimate the absolute values of surface energies in comparison to experiment.The underlying reason for this discrepancy is identical as in the case of vacancies.For Rh, the stability of the lowest energy (111) surface is reproduced by ACE, but the energies of the (100) and (110) are almost degenerate according to DFT while ACE predicts the latter to have a slightly higher energy.Besides these common surface orientations, we also evaluated other surface configurations provided in the Crystalium database [29][30][31] and compared them to the available DFT data in Fig. 4. A close agreement is seen between ACE and DFT with a slight overestimation of Rh surface energies by ACE.When comparing both metals, it is evident that Pt has smaller surface energies than Rh, which is relevant for the structure of nanoclusters (see below).
The phase stabilities of various bulk phases over a broad range of volumes are examined for both metals in Fig. 5.It can be seen that the agreement between ACE and DFT is excellent up to the dissociation limit, which is often not the case for ML models [16,17].Both Pt and Rh are most stable in the  for Rh are all correctly captured and bcc and sc correspond to body-centered cubic and simple cubic phases, respectively.Phonon dispersions along high symmetry directions of the Brillouin zone and phonon densities of states for the fcc phases of Pt and Rh are shown in Fig. 6.ACE exhibits a close agreement with DFT for both elemental metals.
Apart from the properties of the elemental metals, we evaluated the equilibrium lattice parameters of random Pt-Rh fcc alloys as a function of temperature and compared them with available experimental data in Fig. 7(a).The linear decrease of the lattice parameter with increasing Rh concentration is well reproduced by ACE.As mentioned above, the overestimation of the values of the lattice parameters predicted by theory are due to the underbinding of the PBE exchange-correlation functional used for the DFT training data.
There exists experimental evidence for ordering tendencies in the Pt-Rh system [44] but no unequivocal confirmation of the existence of distinct ordered phases.Nevertheless, several theoretical studies [45][46][47][48][49] predicted a possible occurrence of ordered phases at low temperatures.The D0 22 phase and the socalled "phase 40" (NbP prototype, Pearson symbol tI8, space group I4 1 /amd) were identified as the lowest energy phases for the PtRh 3 and PtRh stoichiometries, respectively.More recently, Hart et al. [50] reported the Pd 2 Ti prototype (Pearson symbol oI6, space group Immm) to be a stable phase on the convex hull for the PtRh 2 stoichiometry.
We calculated the formation energies of these ordered structures using ACE and compared them with the reference DFT data in Fig. 7 (b).In addition, the plot contains formation energies of fcc random solid solutions.ACE predictions for the three ordered phases agree closely with the DFT results, reproducing correctly the convex hull.All random solid solutions lie above the convex hull, but their formation energies remain negative.

Simulations of Pt-Rh nanoclusters
As the Pt-Rh nanoclusters appear to be terminated predominantly with (111) and (100) planes (see Fig. 1), we started our simulations using truncated octahedrons that contain six (100) and eight (111) planes.The initial focus of the simulations was to examine the surface segregation that may play a role in stabilization of the core-shell geometry.We investigated two compositions of Pt 30 Rh 70 and Pt 45 Rh 55 , with a random initial distribution of both species in the nanoclusters (see supplementary Fig. S3).The simulations were done at T = 500 K using a hybrid MD/MC procedure for efficiently reaching thermodynamic equilibrium.
Figure 8 shows the final cluster configurations for both compositions.In both cases, we observed a strong tendency of the Pt atoms (blue spheres) to segregate to the surface and to form a continuous monolayer across the whole cluster surface.Based on the results for the Pt-poor composition Pt 30 Rh 70 , we surmise that the segregation of Pt occurs preferentially on the (111) facets as the coverage of the (100) facets remains incomplete (see the top left configuration in Fig. 8).For the Pt-rich composition Pt 45 Rh 55 , all surfaces of the  cluster are covered completely by Pt (see the top right configuration in Fig. 8).Interestingly, the excess Pt atoms in this case do not continue to segregate to adjacent subsurface layers but remain randomly dispersed in the cluster interior.While the formation of a Pt surface layer can be related to the low surface energies of Pt (cf.Fig. 4), experimental observations of the core-shell particles show a thicker coverage of 3-5 layers of Pt at the surface [10].However, according to our simulations Rh atoms have a slight preference for the occupation of the second subsurface layer (see the bottom right configuration in Fig. 8).
To analyze this discrepancy in more detail, we performed additional DFT calculations of Pt-Rh surface slabs with different orientations.All slabs consisted of ten atomic planes with two Pt and eight Rh layers.One of the Pt layers was located at the surface while the position of the second Pt layer was varied, from being the adjacent subsurface layer, i.e., forming two Pt layers on the surface to the opposite surface of the slab, i.e., a Rh slab with both surfaces covered by Pt layers.These calculations confirm that the single Pt layers on both slab surfaces corresponds to energetically most favorable configuration, see supplementary Fig. S4.In contrast, the formation of two adjacent Pt surface layers results mostly in configurations with higher energies than when the Pt layers are farther apart.The energetics of these interactions is described well by ACE, confirming the reliability of ACE predictions in the large-scale MD simulations.This analysis shows that in thermodynamic equilibrium a multilayer segregation of Pt atoms on the surface of the Pt-Rh nanoclusters is likely not favorable.This suggests that the core-shell morphologies with a Rh core and a thicker Pt shell observed in some STEM experiments are kinetically stabilized.This corroborates with the fact that the experimental synthesis procedure creates the core first followed by the coating of the surface element.
To obtain a complementary view on the stability of the core-shell clusters, we performed additional MD simulations at elevated temperatures.In contrast to the hybrid MD-MC simulations presented above, we generated the initial clusters as spheres with a Rh core and a Pt shell, to mimic real experimental geometries.The thickness of the shell was set to approximately three Pt layers.The clusters were then annealed at temperatures of 1000 and 1500 K.
In Fig. 9, the nanoclusters are visualized after 2 ns of annealing.At the lower temperature of 1000 K, the cluster undergoes small shape changes but its core-shell morphology remains intact.When the temperature is increased to 1500 K, we start observing diffusion events leading to a gradual intermixing of the Rh atoms from the core into the Pt shell.The duration of the MD simulation is too short to reach thermodynamic equilibrium and supports our

Conclusion
We presented the development and application of ACE for the Pt-Rh binary system.Our approach utilized a semi-automatic workflow combining highthroughput DFT calculations, efficient fitting and validation procedures, and active learning using uncertainty indication based on the D-optimality criterion.The resulting ACE potential accurately describes fundamental properties of the Pt-Rh system with ab initio accuracy and is suitable for large-scale atomistic simulations.We applied the ACE to investigate structural stability of various Pt-Rh nanoclusters using hybrid MC/MD simulations.Our results show a strong surface segregation of Pt and a preference to form a single monolayer coverage of the whole cluster.This is in contrast to experimental observations, where core-shell nanoclusters with thicker Pt shells were reported.By carrying out MD simulations at elevated temperatures, we showed that the core-shell cluster morphologies consisting of a Rh core with a thicker Pt shell observed in experiments are not thermodynamically favorable but rather kinetically stabilized.
The accuracy of ACE cannot exceed the accuracy of the reference DFT calculations.DFT predicts surface and defect energies in Pt poorly, which suggests that the simulations for Pt-Rh nanoclusters are limited primarily by the DFT reference data.This highlights that improved exchange correlation functionals that are suitable for high-ghroughput database generation are urgently required.

DFT calculations
All reference DFT calculations were done using the Vienna ab-initio simulation package (VASP) [23][24][25] with projector augmented waves [52].For the description of the exchange-correlation functional, we employed the Perdew-Burke-Ernzerhof (PBE) functional [26].For energy and force calculations parameters were chosen with a tight settings to obtain an converged accurate results.A plane-wave cutoff energy of E cut = 500 eV and Gaussian smearing of σ = 0.1 eV were used in all calculations.A dense Γ centered k-point mesh with a spacing of 0.125 Å−1 was used to sample the Brillouin zone.

ACE parametrization
The ACE model was parameterized using the pacemaker package [21].In the first stage, fitting was performed using a hierarchical basis extension, resulting in the ACE potential consisting of 8874 basis functions and 20 688 parameters with a maximum order ν = 4, corresponding to five-body interactions.Additional details of the model are summarized in Table 2. Relative weights of different structures in the loss function were distributed based on their distance to the convex hull.This method prioritized structures with lower energies, giving them a greater contribution to the loss function [21].In addition, the weights of non-periodic clusters were reduced by a factor of ten to provide regularization to the ACE model.The initial parametrization was improved using active learning as described in Section 4.4.

MD and MC simulations
The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) package [53] with the ML PACE [15] and MC [54] packages was used for the MD and MC atomistic simulations.We adopted an integration time step of ∆t = 1 fs for the MD simulations.All simulations were equilibrated at the target temperature/pressure for t eq = 50 ps before production runs.For MD simulations of lattice parameters of Pt-Rh alloys, we used a simulation cell containing 2048 atoms that corresponds to an 8 × 8 × 8 conventional fcc cell.
Simulations were conducted at two temperatures T = 300 K and T = 800 K in an isothermal-isobaric (N P T ) ensemble with a pressure of P = 0 bar for t = 500 ps.Nosé-Hoover thermostat and barostat were used to control temperature and pressure.Damping parameters were 0.1 ps for temperature and 1 ps for pressure to ensure the temperature and pressure fluctuations are appropriate.Two nanocluster geometries were created using the Atomistic Simulation Environment (ASE) [55]: a truncated octahedron for the surface segregation phenomena and a spherical nanocluster for annealing of nanoclusters.The truncated octahedron comprises six (100) and eight (111) surfaces and contained 2406 atoms.Different chemical compositions were initialized by randomly assigning Pt or Rh to the atomic positions.A spherical nanocluster was created with 3577 atoms consisting of a Rh core (1058 atoms) covered by three layers of Pt (2519 atoms).All simulations of nanoparticles were performed in the canonical (N V T ) ensemble using the Nosé-Hoover thermostat with a damping parameter of 0.1 ps.
A Monte Carlo scheme was used for the simulation of surface segregation using the MC [54] packages implemented in LAMMPS to allow for atomic swapping between different types of atoms.The simulations were performed for t = 200 ps for surface segregation and t = 2 ns for annealing of nanoparticles.The open visulization tool OVITO [56] was used for the visulization of nanoclusters.

Active learning
To capture atomic environments relevant for nanoclusters, we used small fcc surface structures with different crystallographic orientations and sizes, namely, slabs with the (100) surface orientation containing 32 atoms, the (110) surface of 56 atoms, and (111) surface orientation containing 36 atoms, respectively.To sample different slab configurations, we employed a hybrid MD-MC approach in the N V T ensemble, running each simulation for 100 ps at three different temperatures (300, 500 and 800 K).We calculated on-the-fly the extrapolation grade γ for each atom every 100th simulation step using the Doptimality criterion.Structures with γ values greater than 5 were selected as extrapolative structures, resulting in a few thousand candidate structures.All simulations were performed using the LAMMPS version, which calculates the extrapolation grades for ACE [16].To select the most representative atomic environments, we utilized the MaxVol algorithm to construct a new active set from the structures generated in the previous step.We gave priority to structures that had multiple atoms that entered the new active set [16] and selected the top hundred candidates for data efficiency.These selected structures were computed with DFT, added to the original training set, and used to upfit the ACE potential.As a result of the above approach, we observed a significant decrease in the extrapolation grades for representative nanocluster (see Supplementary Fig. S3), and a reduction in the error metrics for the energy and forces of the selected surface structures.The energy root mean square error (RMSE) decreased from 31 to 3 meV/atom and the force RMSE from 172 to 84 meV/ Å(see Supplementary Fig. S1).

STEM characterization
A probe corrected Titan Themis 60-300 operated at 300 kV was used for the STEM characterization.For imaging, a HAADF detector with a collection angle of 78-200 mrad was used.The EDS composition maps were acquired with a Brucker Super X detector.
Supplementary information.The supplementary information is provided in the following PDF file.

Fig. 1
Fig. 1 STEM-HAADF micrograph (a) and corresponding EDS composition map (b) of a Rh-Pt core shell nanoparticle.Pt is in blue and Rh in green.In (a) the crystallographic facets are indicated.

Fig. 2
Fig. 2 Energy per atom of structures in the Pt-Rh training set as a function of the nearest neighbour distance (NNB).

Fig. 3
Fig. 3 Comparison of energies and forces predicted by ACE with respect to the reference training DFT data.The "Low" subset includes data within a distance of 1 eV/atom above the convex hull.

Fig. 4
Fig. 4 Comparison of various surface energies obtained from ACE and DFT.The reference DFT data are taken from the Crystalium database[29-31].

Fig. 5 Fig. 6
Fig. 5 Energy as a function of the nearest neighbour distance (NNB) for various bulk phases of (a) Pt and (b) Rh.The lines represent ACE predictions while the dots correspond to the reference DFT data.The insets depict low-energy regions around the fcc ground state.

Fig. 7
Fig. 7 (a) Lattice constants of Pt-Rh alloys as a function of chemical composition at different temperatures.The square symbols represent the experimental values at T = 300 K [51].(b) Formation energies of most stable ordered phases from DFT (blue diamonds)[50] and their ACE predictions (red crosses).ACE predictions for fcc random solid solutions are marked by red circles.

(a) Pt 30 55 Fig. 8
Fig. 8 Visualization of the surface segregation in Pt-Rh octahedron nanoclusters with (a) Pt 30 Rh 70 and (b) Pt 45 Rh 55 composition.Pt atoms are colored in blue and Rh atoms in green.Top figures show the front view of truncated octahedron with four (111) and one (100) surface facets.Bottom figures show mid cross-sections of the clusters with four (111) and two (100) surface planes.

Fig. 9
Fig. 9 Visualization of core-shell nanocluster after t = 2 ns annealing at (a) 1000 K and (b) 1500 K. Rh atoms in green, Pt atoms in transparent blue.
Fig. S1 Energy (top figure) and forces (bottom figure) error metrics on selected surface configurations before and after upfitting of ACE potential.

Fig
Fig.S2A comparison showing the extrapolation grades for the surface atoms of the octahedron clusters.The top configuration was evaluated using the initial version of the ACE parametrization, the bottom configuration was simulated using the retrained ACE parametrization after the active learning procedure.

(a) Pt 30
Fig. S3 Visualization of the initial configurations of the octahedron nanoclusters used in the hybrid MD/MC simulations.Pt (blue) and Rh (green) atoms are randomly distributed in the nanoclusters.
Fig.S4Dependence of formation energies of surface slabs consisting of two Pt and eight Rh atomic layers.The first Pt layer is always positioned on the surface (index 1) while the position of the second Pt layer varies (i.e., index 2 for the 2nd Pt layer corresponds to two adjacent Pt layers at one surface of the slab while index 10 corresponds to a Rh slab covered on both surfaces with Pt).The formation energies are computed with respect to energies of the ground state fcc Pt and Rh.

Table 1
Comparison of various properties of elemental Pt and Rh obtained from ACE, DFT and experiment.For the vacancy, E f and E m represent the formation and migration energies, respectively.For surface energies, the experimental values are for average orientations.

Table 2
Details of the ACE parametrization: r c denotes the cutoff radius; RBF the type of radial basis function; "#func./element" represents the maximum number of functions per element for each body order ν.