Multi-scale modeling and experimental investigation of oxidation behavior in platinum nanoparticles

Abstract

Understanding the impact of oxidation on the reactivity and performance of Pt nanoparticles (NPs) is crucial for developing durable and efficient catalysts. In this study, we investigate the oxidation process of a realistic Pt NP using a multistep approach combining computational methods (ReaxFF, MACE-MP-0, and DFT) with experimental techniques (XRD, TEM, and EXAFS). Our workflow aims to measure oxidation extent, compare different computational models, analyze electronic structure changes, and provide guidance for selecting appropriate computational models in catalytic atomistic studies. We perform hybrid MD-MC simulations using ReaxFF which reveal significant oxidation with oxygen penetrating deep into the core at high oxygen partial pressure, with the formation of detached small cluster oxide Pt₆O₈ species. We investigate the plausibility of these configurations and possible degradation mechanism by carrying out XRD, TEM, and EXAFS measurements on samples of various average particle sizes. Experimental measurements show partial agreement with our simulations in terms of coordination numbers, bond distances, oxygen fractional occupancy and onset/place-exchange potentials. Despite these agreements, we find poor matches between the binding energies calculated by ReaxFF and DFT, casting doubt on the predictive power of ReaxFF and the existence of Pt₆O₈ species. In contrast, the universal MACE-MP-0 model shows significant improvement in the prediction of energetics. Comparing these force fields with DFT calculations on oxidized and non-idealized systems is essential for understanding the limitations of such models in predicting catalytically relevant properties at high potentials and was previously unexplored in the literature. Our study provides a foundation for understanding the complex interplay between nanoparticle structure, oxidation state, and catalytic performance, aiming to guide the rational design of advanced catalytic materials through atomistic modeling.

---

1 Introduction

Platinum nanoparticles (Pt NPs) have emerged as important catalysts in various industrial applications, including fuel cells, catalytic converters, and chemical synthesis. However, the long-term performance of Pt NPs can be compromised when exposed to reactive conditions such as the acidic environment in proton-exchange membrane hydrogen fuel cells, potentially leading to a loss of catalytic activity.^1–4 Oxide species (O*, OH*) are present at high potentials due to electrochemical oxidation reactions occurring close to 0.8 V.^5,6 These species are thought to be important site-blockers for the oxygen reduction reaction (ORR), and oxidized catalysts generally show hindered activity, mainly due to site blocking and altered kinetics.⁷ For the ORR oxide saturation is also known to lead to hydrolysis i.e. parasitic, non-electrochemical reactions which might contribute to the high onset potential on platinum-based catalysts.⁸ This decline in performance presents challenges for the development of durable and efficient catalysts, emphasizing the need for a better understanding of the complex relationship between reactivity conditions and catalyst performance.^9,10

Experimental studies on the oxidation of platinum nanoparticles have faced challenges, especially when trying to directly link the oxidation state to the resulting chemical reactivity. Early studies, primarily focused on low-index platinum surfaces in ultra-high vacuum (UHV) conditions^11,12 agreed that atomic oxygen binds on Pt(111) forming a (2 × 2)-O structure at low temperature,^13–15 with some persistence at higher temperature.¹⁶ These studies also indicated that molecular oxygen is unlikely to be present on Pt(111) at room temperature.^12,16 Platinum catalysts can undergo oxidation even at low pressure, inducing modification of electronic properties and structural changes, such as surface buckling.¹⁷ Recent in situ studies have advanced our understanding by investigating platinum oxidation under near-ambient pressure conditions for longer exposure times. Fantauzzi et al.¹⁸ reported the oxidation of Pt(111) up to 3 monolayers (ML) through a combination of experimental and theoretical approaches. Furthermore, a study on stepped platinum surfaces revealed the formation of platinum-oxide nanoclusters at a pressure of 1 Torr.¹⁹ Finally, an in situ study by Ellinger et al.²⁰ at a partial pressure of 500 mbar and high temperature showed the growth of two atomic layers of bulk-like platinum oxide.

Modern platinum catalysts often consist of deposited platinum nanoparticles, which possibly exhibit complex structures that differ from their perfect counterparts.²¹ These structural differences can significantly influence the oxidation behavior and catalytic properties.²² The oxidation of these platinum nanoparticles has been extensively studied in the literature, with the formation of platinum-oxide phases reported in numerous studies.^23–29 Saveleva et al.²⁵ reported the growth of mixed PtO and PtO₂ layers when studying the electro-oxidation of Pt nanoparticles under applied voltage, at 450 K. Additionally Imai et al.²⁶ characterized the initial platinum-oxide phase to have an α-PtO₂ local structure which transitioned to a β-PtO₂ structure as the oxidation time increased under sustained potential. The presence of an ordered oxide phase was confirmed by other studies under various experimental conditions, with a high dependence on the average particle size, where smaller particles were found to be more prone to oxidation.^27–29

These studies have provided valuable insights into the oxidation behavior of platinum nanoparticles, but using experimental techniques the underlying mechanisms and implications for catalytic activity are often out of reach. As a result, computational studies have emerged as a valuable tool to complement experimental approaches, providing insights into the oxidation process, and helping to overcome the limitations of traditional techniques. Many theoretical studies have supported the claims made by experimental investigations, elucidating the electronic and structural changes induced by oxygen adsorption on platinum catalysts.^30–32 However, the computational modeling of realistic nanoparticles remains largely unexplored, with most studies focusing on idealized surfaces or perfect nanoparticles,^33,34 leaving room for further computational investigation into the complexity of these systems. Here, we investigate the oxidation behavior of an experimentally reconstructed platinum nanoparticle shown in Fig. 1 which was obtained from scanning transmission electron microscopy (STEM) experimental work.³⁵ These experimentally reconstructed nanoparticles have been shown to exhibit unique structural features not present in perfect nanoparticles.³⁶

Fig. 1 Experimentally reconstructed platinum nanoparticle used in this study. The nanoparticle consists of 353 atoms and was obtained from STEM images.

We use a multistep approach that combines Monte Carlo simulations with molecular dynamics using a fast ReaxFF³⁷ forcefield, to provide sufficient sampling. We then validate our results by comparing the coordination numbers and bond lengths with experimental data obtained from X-ray diffraction (XRD), transmission electron microscopy (TEM) and extended X-ray absorption fine structure (EXAFS) measurements. We further discuss results obtained with the ReaxFF forcefield by progressively climbing the accuracy ladder, e.g. using the recent MACE-MP-0 model³⁸ and linear scaling DFT calculations. Finally, we present our results which will aim to highlight multiple key aspects such as the electronic structure modifications of the nanoparticle under oxidation conditions and the difference in thermodynamic and geometric predictions between both forcefields.

2 Computational methods

Our workflow is conducted on a platinum nanoparticle composed of 353 atoms which was obtained through a three-dimensional reconstruction process using STEM images from a study of Aarons et al.³⁵ Our workflow is divided into 5 steps as shown in Fig. 2 and will be detailed below.

Fig. 2 The computational workflow followed in this study. Monte Carlo simulations provide initial configurations, molecular dynamics simulations are then performed to relax the geometries and further explore the phase space. Finally, ONETEP linear-scaling DFT calculations are carried out to obtain the electronic structure and refined energetics of the oxidized nanoparticles.

(1) We perform hybrid Monte Carlo-molecular dynamics (MC-MD) simulations in the grand canonical ensemble using the LAMMPS software.^39,40 We use a Pt/O/H/C/Ni ReaxFF forcefield developed by Gai et al.,⁴¹ which is originally based on a forcefield from Fantauzzi et al.⁴²

This forcefield has been previously used in studies focusing on the oxidation of platinum surfaces and nanoparticles using various Monte Carlo minimization algorithms.^18,31,32,41 In this study we deviate from these approaches by incorporating dynamics for the platinum atoms. Grand-canonical moves are then performed while the platinum nanoparticle is allowed to move, sampling the canonical ensemble. The main parameter of these simulations is the oxygen chemical potential of the gas reservoir which is calculated in eqn (1) from the chemical potential of O₂, assuming ideal gas and adding the energy of an oxygen molecule calculated with the ReaxFF forcefield for reference consistency.

is the standard Gibbs free energy at the given temperature for an oxygen molecule, calculated from and , the standard enthalpy and entropy respectively, both taken from the NIST-JANAF thermodynamics tables.⁴³P is the pressure, P° is the reference pressure. Throughout the study, we use the pressure as a reference, since it is the only variable term in eqn (1). We study a broad range of pressures from 10⁻²⁵ to 1.0 atm at the same temperature of 350 K.

(2) After our MC-MD simulations, we randomly extract 32 geometries throughout the highest pressure run at evenly spaced oxygen coverages, which we can then study systematically.

(3) The geometries are then subjected to an additional 100 ns of NVT simulations at 350 K. This is done to amplify potential geometrical features such as small platinum-oxide clusters that were previously observed in past studies using similar forcefields. The grand canonical Monte Carlo moves are stopped, and the system is allowed to relax in the canonical ensemble.

(4) We then retrieve relaxed geometries after the NVT simulations by selecting the most stable configurations in the last 1 ns of the simulation for each trajectory. These geometries are then optimized using the ReaxFF forcefield and the recent MACE-MP-0 foundation model separately. The MACE-MP-0 model³⁸ is based on a message-passing version of the atomic cluster expansion (ACE)⁴⁴ and is trained on the material project database. MACE was recently benchmarked to be one of the most accurate forcefield to predict the energetics of materials across the periodic table.⁴⁵

(5) DFT calculations are then performed on the relaxed geometries by both forcefields to compare them and obtain the electronic structure of the oxidized nanoparticles. Computationally elucidating the electronic structure of such large platinum nanoparticles remains a challenging task, often hindered by the computational expense associated with such calculations. We overcome these limitations by using the ONETEP linear-scaling DFT software,⁴⁶ which exhibits a linear scaling behavior with respect to the number of atoms. This enables us to readily access electronic structure calculations for a wide range of nanoparticle sizes and configurations.⁴⁷ We make use of the PBE functional⁴⁸ without using any dispersion corrections. Additional details about ONETEP and the parameters used in this study can be found in the ESI† along with the parameters used for the ReaxFF and MACE-MP-0 forcefields.

3 Results

3.1 Monte Carlo simulations

Monte Carlo simulations are run until accepted additions and deletions exchange moves stabilize to the same count over the last 5 × 10⁷ exchange moves. This criterion is not applied for higher pressures (1.0, 10⁻², and 10⁻⁵ atm), where the simulations are stopped before reaching this criterion due to low deletion acceptance rate. Throughout the study we will refer to the oxidized nanoparticles using the concept of oxygen ratio N_O/N_Pt instead of coverage. This choice is motivated by the difficulty of defining a proper surface area for the nanoparticle and the presence of subsurface oxygen atoms.

Atomic oxygen tends to adsorb on low-index platinum surfaces on hollow sites (3-fold or 4-fold),⁴⁹ this is highlighted in Fig. 3 where we report that for most pressures, more than 70% of the oxygen atoms are in a 3-fold coordination. We also report a low number of atomic oxygen binding on bridge sites (2-fold coordination) which does not seem to be impacted by the pressure, indicating that these are specific sites, where oxygen opportunistically binds. After visual inspection, these sites are often found on dangling platinum atoms or on the edge of the nanoparticle. These sites were previously reported to be the most stable sites for atomic oxygen binding by DFT studies on smaller nanoparticles.^49,50 Finally, as the pressure increases 4-fold subsurface sites become more populated, which occurs between 10⁻¹⁰ and 10⁻² atm when the outer coverage is between 0.7 and 0.9 ML. As suggested by previous DFT studies^51–53 subsurface sites require sufficient lattice disruption to be populated, which is the case here due to oxygen adsorption at the outer shell.

Fig. 3 Pressure dependence of oxygen coordination number population. The percentage is calculated with respect to the total number of oxygen atoms for each pressure. Lines are added to guide the eye.

We report the dependence of oxygen ratio versus pressure in Table 1 which further highlights that the strong adsorption energy of atomic oxygen on Pt catalysts leads to sustained coverage even at extremely low pressures (10⁻²⁵, 10⁻¹⁵ atm). We compare our results to similar studies using similar methods and systems, and also include an approximation of the outer coverage which is done by computing the alpha shape of oxygen atoms and platinum atoms separately. The coverage is then calculated as the ratio of these two atomic counts. At low pressures, our results compare relatively well with the study of Gai et al.⁴¹ which was done at a lower temperature on a Pt(111) surface. Discrepancies between our results and Pt(111) surface widen when considering the experimental study of Getman et al.⁵⁴ where a coverage of 0.48 was reported at a pressure of 10⁻⁴ atm and a temperature of 353 K. A restriction on coverage is identified by the authors to be due to kinetic limitations of O₂, which becomes a less effective oxidant as the oxygen coverage increases. When switching to NO₂, a kinetically favorable oxidant, the authors found a coverage of 0.75 ML, closer to our results. We note that our results are not completely comparable to clean surface studies since our nanoparticle has a different surface area and binding sites, generally displaying slightly higher coverages due to stronger binding on a nanoparticle of this size.⁵⁵

Table 1 Summary of oxygen fractional occupancy N_O/N_Pt at different pressures. Our values are reported as oxygen ratio along with an approximation of the outer coverage (calculated as described in the text) in parentheses. Our values are compared to similar studies using similar methods and systems, all reported as coverage, unless specified otherwise

Pressure (atm)	1.0	10^−2	10^−3	10^−5	10^−10	10^−15	10^−25
a Gai et al.^41 Pt(111) (300 K). b Octahedral NP (500 K). c Cuboctahedral NP (500 K). d Cubic NP (500 K). e Kirchhoff et al.^31 3 nm cuboctahedral NP (400 K). f Reported as oxygen ratio.
This work	1.20 (0.90)	1.14 (0.90)	—	0.73 (0.91)	0.47 (0.75)	0.37 (0.59)	0.23 (0.36)
Ref^a	0.8	0.72	—	0.7	0.65	0.5	—
						0.40^b
						0.36^c
						0.26^d
Ref^e	—	—	1.38^f (1 mbar)	—	—	—	—

---

Finally, our simulations systematically display a lower oxygen ratio than the studies of Kirchhoff et al.³¹ at similar pressures and temperatures despite using a similar ReaxFF forcefield. We attribute it to our methodology which differs from the Monte Carlo minimization scheme used in these studies. This minimization scheme includes a full geometry relaxation carried before the acceptance criteria, ultimately leading to faster convergence and higher oxygen ratio, but does not correctly sample the canonical ensemble and resembles more a basin-hopping algorithm. Our current hybrid MC-MD methodology was not driven by the need to reach specific coverages, but to also include the dynamics of the platinum atoms, leading to correct canonical sampling, despite lower acceptance rate, and higher computational cost. In addition, including the dynamics of the Pt skeleton was motivated by previous report of small platinum-oxide cluster forming, potentially leading to a degradation mechanism.

To understand the structure of the oxidized nanoparticle we plot the oxygen distribution profiles in Fig. 4. The plots represent the spherical density or concentration of oxygen atoms as a function of distance from the center of mass. As the pressure increases, the amount of oxygen penetrating inside the nanoparticle increases, first in subsurface sites at 10⁻⁵ atm, and then in the core region for pressures above 10⁻² atm, in agreement with the previous coordination number analysis shown in Fig. 3. Additionally, the same plots reveal a significant amount of lattice strain caused by subsurface and core oxygen atoms at higher pressures. At a pressure of 10⁻¹⁵ atm, the outer shell boundary is clearly visible and located approximately 11 Å away from the center of mass. As the pressure increases, the boundary becomes less defined and oxygen atoms are found at distances up to 20 Å away from the center of mass.

Fig. 4 Distribution of oxygen with respect to the center of mass of the oxidized nanoparticle. Subsurface and core oxygen atoms are observed as the pressure increases from 10⁻¹⁵ to 10⁻² atm. No normalization has been applied. Oxygen coordination averaged over platinum atoms located around 1 Å of each radial subshell is shown as color gradient. The figure has been created by averaging over the 100 ns of NVT simulations in step 3, but is shown here for clarity.

Typically, oxidation of metallic nanoparticles is thought to proceed from the surface inwards, with the formation of an oxide shell surrounding a metallic core. The presence of oxygen atoms near the center suggests that under sufficiently oxidizing conditions, platinum nanoparticles can become completely oxidized throughout their entire volume. Various experimental studies have shown that this is most likely not the case for larger nanoparticles (>2.0 nm).^29,56 For smaller nanoparticles with diameters close or lower than 2.0 nm, previous studies suggest that core oxidation is possible under specific conditions.²⁷ To further validate these claims and our results, three carbon black-supported Pt catalysts were prepared and characterized by XRD, TEM, and XAS (EXAFS) with the aim of determining crystalline domains, average Pt–Pt/O distances, and coordination numbers as a function of particle size distribution. Measurements were made on prepared catalysts which had been chemically reduced during the synthesis process and subsequently passivated by exposure to air. Our results are summarized in Table 2, the size distributions as measured by TEM gradually increase as a function of Pt loading, XRD measurements showed that the 10 and 20 wt% Pt samples did not show any discrete crystalline Pt domains suggesting that the Pt NPs are both disordered and contain a significant amount of oxide. The 50 wt% Pt sample did show a crystalline Pt domain profile, but the average crystallite size was estimated to be less than 2 nm. The EXAFS CNs support these results with the 10 wt% Pt sample showing Pt–O but no significant Pt–Pt character, the 20 wt% Pt sample showing a mix of Pt–Pt and Pt–O CNs suggesting very small metallic Pt domains which are significantly oxidized. The 50 wt% Pt sample shows significant Pt–Pt character with a relatively small Pt–O component, indicating the presence of metallic Pt NPs with a defined passivated oxide layer at the surface. The Pt–Pt distances for the 20 and 50 wt% Pt samples are 2.76 ± 0.01 Å consistent with that shown by bulk Pt–Pt (2.77 Å), and the Pt–O distances vary between 1.97 and 2.03 Å. Our Monte Carlo simulations are done on a nanoparticle consisting of 353 atoms, whose size is calculated to be 2.03 nm, falling in between the 10 and 20 wt% Pt samples. In Table 2 our simulation done at 1.0 atm shows reasonable agreement with the experimental results, the Pt–Pt coordination number is low, due to the significant oxidation predicted by the Monte Carlo simulations. Similarly, the Pt–O CN is close to the expected range for smaller particle distributions (between 3.4 and 3.6). However, the theoretical Pt–Pt average bond distance is higher than the experimentally measured value, which might be due to higher deformation and lack of a defined oxide phase in the theoretical model. We note that our theoretical coordination numbers were calculated using a distance cutoff determined by the first minima of the corresponding radial distribution function: 3.04 Å for Pt–Pt and 2.3 Å for Pt–O. Due to the particularly high sensitivity of the CNs to the cutoff the comparison with experimental estimate should be taken with caution. The experimental bond distances suggest the presence of Pt(IV) rather than Pt(II) (2.01–2.03 vs. 2.07 Å).⁵⁷ Additionally, the mixed oxidation state Pt oxide, Pt₃O₄ shows Pt–O distances of 1.97–2.00 Å.^58,59

Table 2 Experimental and theoretical coordination numbers and bond distances for Pt–Pt and Pt–O. Theoretical values are averaged over all atoms, which is closer to the experimental definition. The Pt₃₅₃ diameter was estimated by calculating the averaged maximum distance between every two atom pairs of the outer platinum shell

	Cryst. size (nm)	Particle size (nm)	Coord. number		Bond distances (Å)
			Pt–Pt	Pt–O	Pt–Pt	Pt–O
10 wt%	—	1.7 ± 0.5	—	3.6 ± 0.2	—	2.03 ± 0.01
20 wt%	—	2.4 ± 0.7	2.4 ± 0.7	3.4 ± 0.3	2.76 ± 0.01	2.02 ± 0.01
50 wt%	<2	2.7 ± 1.3	7.8 ± 0.6	1.9 ± 0.4	2.77 ± 0.01	1.97 ± 0.01
Pt353 MC (1.0 atm)	—	2.03	2.68	3.63	2.94	2.05
Pt353 MC (10^−25 atm)	—	2.03	8.57	0.67	2.82	1.97

---

In contrast with our results, various studies on similar nanoparticle sizes have reported oxidation to occur mainly at the outer platinum shell. Takagi et al.²⁸ performed in situ hard X-ray photoelectron spectroscopy (HAXPES) measurements at various potentials and reported that for a sample with a size distribution centered around 2.6 nm, only the outer platinum shell was oxidized. The presence of oxygen atoms near the center of the nanoparticle in our simulation is justified by the nature of the Monte Carlo simulations which disregard kinetic aspects. Hence, atomic rearrangement required for oxygen to penetrate deep into the core is bypassed i.e. no activation energy of oxygen diffusion into the platinum lattice, which otherwise, renders this process kinetically unfavorable.⁶⁰ Our experimental evidence along with previous experimental work, suggests that the size range between 1.5 nm and 2.5 nm represents a critical transition region where the oxidation behavior changes dramatically.²⁷ Within this narrow size window, nanoparticles appear to shift abruptly from partial surface oxidation to complete oxidation, indicating a sharp threshold effect. The precise point of this shift likely depends on specific experimental conditions such as temperature, type of oxidant and partial pressure.

3.2 Molecular dynamics

We further progress in our workflow by performing molecular dynamics simulations at 350 K for 100 ns. Our interest in extended dynamics can be explained by considering previous literature using similar forcefields where small platinum-oxide clusters have been reported to form.⁴¹ As illustrated in Fig. 5, we also report similar dissociation-like events taking place during the dynamics and leading to a number of partially detached Pt₆O₈ species on the nanoparticle. The occurrence of these events strongly depends on the oxygen fraction occupancy, with higher occupancies leading to a greater number of detachments. Each of these events is accompanied by an energy reduction as schematically shown in Fig. 5a. These species were first reported by Kirchhoff et al.^31,32 where the authors investigated the possibility of a stability region for the small cluster in the aqueous phase with specific charge and spin states.

Fig. 5 An example of a dissociation event taking place during the dynamics at high oxygen fractional occupancy (N_O/N_Pt > 0.6). These events lead to an energy reduction (a) and geometrical deformation shown by an increase in mean square deviation (msd) of nuclear coordinates (b). The final state of the nanoparticle after 100 ns is shown with an oxygen ratio of 1.07 (c). A Pt₆O₈ species is partially detached from the nanoparticle (d).

The resulting geometry after 100 ns of NVT simulations for a high oxygen ratio is shown in Fig. 5c where the nanoparticle geometry exhibits significant structural deformations, the morphology appears highly distorted, with pronounced surface irregularities and a loss of the expected quasi-spherical or faceted shape. The structure displays extensive porosity, with the formation of voids and channels permeating throughout the entire nanoparticle. This gives rise to fragmented, loosely bound domains interconnected by these more stable oxidized platinum clusters. The Pt₆O₈ species emerges as a distinctive motif, standing out in contrast to the overall amorphous structure. Importantly, no other such well-defined cluster species can be clearly distinguished within the nanoparticle structure. These structural perturbations collectively contribute to an increased surface-to-volume ratio which should elevate the surface free energy, with the Pt₆O₈ units serving as structural anchors. At high oxygen ratios, 100 ns of extended dynamics is not sufficient to fully equilibrate the system, and as shown in Fig. 5a the energy of the system is still decreasing due to the dissociation events produced by the ReaxFF forcefield.

If, as described by the forcefield, these events are taking place frequently within 100 ns of simulation, it would be expected to find experimental evidence of this species even in UHV conditions for pressures above 1 Pa. Despite this, we were not able to find experimental evidence of the existence of the Pt₆O₈ species. Instead, platinum dissolution is believed to occur at the interface at specific potentials, and to be an atomic process that might be influenced by oxide formation, but not via Pt₆O₈ species.^61,62 The presence of this species will be further investigated in the next sections, by performing geometry optimizations with MACE-MP-0 and subsequent DFT calculations.

3.3 Geometry optimizations

We selected the last snapshots from the previous molecular dynamics runs for all our oxygen ratios and performed geometry optimizations. For comparison, we also performed geometry optimizations using the recent MACE-MP-0 model using the same geometries. We quickly noticed a global trend where the MACE-MP-0 forcefield attempts to put back the Pt₆O₈ species together with the rest of the nanoparticle, in contrast with the ReaxFF forcefield which keeps the Pt₆O₈ species in a semi-detached state.

An example of the process is shown in Fig. 6 where the Pt₆O₈ species is packed back together with the rest of the nanoparticle. Although geometry optimizing will only partially change the initial geometry from ReaxFF, these results indicate that MACE-MP-0 would give geometries that do not include Pt₆O₈ species. In order to systematically compare geometries between the two forcefields we define the concept of sphericity score which quantifies the deviation from a perfect spherical shape by calculating the variance in the distances between the centroid and the platinum atoms at the nanoparticle surface.

S = Var[R([scr S, script letter S]_outer)] (2)

Fig. 6 Illustration of the commonly observed phenomena during a MACE-MP-0 geometry optimization, the Pt₆O₈ species is packed back into the nanoparticle.

[scr S, script letter S]_outer is the set of atoms belonging to the outer shell, obtained by computing the alpha shape of the nanoparticle, and R([scr S, script letter S]_outer) represents the distances from each atom in [scr S, script letter S]_outer to the centroid. The obtained scores are then normalized between 0 and 1 within the dataset, with 0 being the least spherical and 1 being the most spherical nanoparticle. Fig. 7 shows the sphericity scores for the entire dataset. Generally, ReaxFF produces geometries with lower sphericity scores, while the geometries optimized using MACE-MP-0 are consistently more spherical. In addition, our sphericity concept is a good measure of the geometrical deformation of the nanoparticle against oxygen ratio. This tendency towards more spherical geometries suggests that MACE-MP-0 favors more compact and densely packed structures, compared to the more irregular structures predicted by ReaxFF.

Fig. 7 Sphericity scores calculated for geometries predicted by both forcefields against oxygen ratio, the dataset has been scaled between 0 (least spherical) and 1 (most spherical). As the oxygen ratio increases, particles are less spherical. MACE-MP-0 geometry optimizations consistently make the geometries more spherical.

With both of our forcefields having different behaviors in terms of geometrical deformation, we expect that thermodynamic properties will also be different. We further explore these disparities by focusing on thermodynamic aspects, using eqn (3). The equation makes use of the computational hydrogen electrode (CHE) scheme which enables extrapolation of the free energy at different potentials as a posteriori correction.⁶³ Additionally, this approach also avoids calculations involving atomic or molecular oxygen, which are known to be problematic for DFT.⁶⁴

where G_nO* is the Gibbs free energy of the nanoparticle with n oxygen atoms, G_* is the Gibbs free energy of the nanoparticle with no oxygen atoms. G_H2O and G_H2 are the Gibbs free energies of water and dihydrogen in their standard state, respectively. U_SHE is the potential of the standard hydrogen electrode. The Gibbs free energies in eqn (3) are often reduced to G = E_DFT + E_ZPE − TS_vib for adsorbates.⁶⁵ Considering the large size of our systems we decided to ignore zero-point and thermal contributions for the oxidized nanoparticles and reduce the Gibbs free energy to the DFT energy, G ≈ E_DFT. For reference species, the total energy is computed with DFT on the isolated molecule, zero point energies (ZPEs) are available elsewhere⁶³ (0.27 eV for H₂, 0.56 eV for H₂O), and entropy contributions are taken from the NIST-JANAF thermodynamic tables.⁴³ In practice eqn (3) is related to the total (ΔG_O) or per atom (ΔḠ_O) binding energy of oxygen atoms on the nanoparticle, using water and hydrogen as reference states.

We report binding energies per oxygen atom (ΔḠ_O) in Fig. 8 for geometry optimized structures using each forcefield. As expected, the binding strength lowers as the coverage increases as the interaction between adsorbates intensifies.¹⁶ From the same figure, two regimes can be distinguished, starting with an initial linear progression for both forcefields. ReaxFF then predicts a near-constant regime for oxygen ratios above 0.7, while MACE-MP-0 is more moderate but still displays a plateauing behavior. This phenomenon was also observed in previous studies of oxidized platinum nanoparticles including experimental evidence.⁶⁶

Fig. 8 Oxygen binding energy per oxygen atom plotted vs. oxygen ratio calculated using eqn (3). Two regimes can be observed, a linear, monotonic regime for ratios below 0.7, which then plateau at higher ratios. Compared to MACE-MP-0, ReaxFF predicts stronger binding at higher coverages.

Finally, it is possible to calculate the thermodynamic expected oxidation state at a given potential. To do so, we select evenly spaced points in the ΔG_Ovs. N_O/N_Pt fitting. For each of these points, we calculate the potential dependence using the CHE on a fine grid of U for each oxidized configuration. This process yields a set of linear curves representing the potential dependence for the selected systems. By finding the lowest curve in this set for multiple potential subranges, we determine the thermodynamically favored state at each potential. As a result, we obtain an oxygen ratio vs. potential curve which represents the expected oxidation state across a range of potentials. This process is carried out for both forcefield and shown in Fig. 9. Starting at 0.7 V, oxidation steadily grows as the potential increases, which is the expected behavior and in reasonable agreement with previous studies.⁶⁷ Our value of onset potential for oxidation is also in reasonable agreement with the one of 0.8 V reported in the literature.^67–70 At high oxygen coverages, we observe a notable difference between the predictions of ReaxFF and MACE-MP-0. When the potential approaches 1.05 V ReaxFF predicts a sharp and almost vertical progression. We link this phenomenon to the experimentally observed place-exchange mechanism, where the oxidation starts to complete, and the nanoparticle evolves towards an oxide phase. This phenomenon was experimentally reported occurring at ≈1.05 V on platinum nanoparticles and surfaces.^71–74 In our framework the event is transcribed via thermodynamic arguments and is due to the plateauing regime of the Gibbs free energies vs. oxygen ratio. In contrast, MACE-MP-0 predicts a more gradual increase in oxidation, most likely because it does not recognize the oxide phase formed by the ReaxFF forcefield. The potential reported by ReaxFF for the place-exchange mechanism (≈1.05 V) is in very good agreement with the reported experimental values. Results before the molecular dynamics (step 3) are also included to highlight the effects that the additional platinum-oxide clusters have on the predicted stability diagram. For ReaxFF we see that the curve is shifted towards lower potentials with the difference between the two curves becoming larger as the potential increases. MACE-MP-0 displays the inverted trend and smaller differences, as the forcefield does not predict the species.

Fig. 9 Oxygen ratio vs. potential curves obtained by following the procedure described in the text. Results before the molecular dynamics (step 3) and their MACE-MP-0 geometry optimized version are shown for comparison.

Being able to produce an accurate version of this curve is of paramount importance since oxygen population is suspected to be the cause of the high onset potential of the ORR on Pt for both thermodynamic and kinetic reasons.⁷⁰ A catalyst that shifts this curve by 0.1–0.2 V would be expected to outperform compared to current catalysts for the ORR. This shift would require the binding energy of the oxygen atoms to be reduced by the same amount across the whole potential range, which was previously done by tuning the nanoparticle size, composition or shape.⁴ Performing similar simulations using a universal potential such as MACE-MP-0 would allow for the exploration of these effects, along with external factors such as a support or a solvent, given that these more advanced forcefields retain reasonable accuracy. Additionally, going beyond the CHE approximation will be required to produce an accurate version of Fig. 9. For example, our simulations do not account for electrochemical phenomena i.e. electric fields that arise due to the potential drop at the interface. Similarly, in the CHE framework, electrochemical adsorption is represented by a coupled electron-proton transfer,⁷⁵ most likely misrepresenting any electrosorption valency phenomena normally involved in these adsorption processes.⁷⁶

3.4 ONETEP calculations

In an attempt to explain the underlying physics for the slope changes of the Gibbs free energy in Fig. 8, we compute electronic properties by performing single-point energy evaluations at the DFT level. Due to the cost of these calculations, we only selected half of our previous dataset of 32 oxygen ratios previously optimized by ReaxFF and MACE-MP-0. This also allows for a direct comparison of the binding energies predicted by the two forcefields. Fig. 10 reveals that ReaxFF tends to overestimate the binding energy ΔG_O, while MACE-MP-0 is in better agreement with DFT. We note that for all our tested oxygen ratios, a geometry optimization with the MACE-MP-0 forcefield also leads to an energy reduction at the DFT level. With structure containing more than 700 atoms at high oxygen ratios we report important total binding energy difference in Fig. 10. When normalized per atom, the observed differences for MACE-MP-0 fall within the reasonable range of 0.25 eV per atom. Considering how far these structures are from the convex hull of Pt and O these values are probably current state-of-the art for any non-fine tuned forcefield.⁴⁵ This is not the case for the ReaxFF forcefield, showing values of up to 0.7 eV per atom at high oxygen ratios. This suggests inability of the forcefield to predict the thermodynamic properties of such oxidized nanoparticles. Furthermore, this reinforces the hypothesis that the Pt₆O₈ species are most probably an artifact due to the parametrization of the ReaxFF forcefield on specific bulk platinum oxides.

Fig. 10 Total binding energies ΔG_O predicted by both forcefields plotted against DFT. Colors are used to represent the binding energy difference per atom between the forcefields and DFT ΔΔḠ^DFT/FF_O. The dashed line serves as a reference for the perfect agreement between the forcefields and DFT.

It is interesting that despite the large differences in binding energies with DFT, and apparent inability of ReaxFF to predict the thermodynamic properties of the oxidized nanoparticles, various predicted properties, such as the oxygen ratio vs. potential curves are in good agreement with the expected behavior. This aspect will be discussed in-depth in the last section.

In an attempt to explain the saturation regime, i.e., the slope change of ΔG_O at high oxygen ratios predicted by both forcefield, we computed the projected density of states (PDOS) of the Pt 5d orbitals (Fig. 11a) and the O 2p orbitals (Fig. 11b) for a range of oxygen ratios. As the oxygen ratio increases, the PDOS of the Pt 5d orbital displays a significant peak increase just below the Fermi energy centered around −2.0 eV. In parallel, we observe an important decrease in intensity for states at lower energies as the oxygen ratio increases. With oxygen addition, the Pt 5d PDOS displays a decreased metallic character as the peak becomes narrower and more localized, indicating a change of the electronic structure towards a platinum-oxide-like cluster. We report minimal change in the d-band center, despite the significant changes in the PDOS, and the apparent large changes in oxygen binding energy. This suggests that the d-band center is not a good descriptor for our systems, most likely due to the complete oxidation of the nanoparticle. Such behavior indicate altered Pt–O character due to a change in the nature of the Pt–O bond and do not follow the conventional model brought forward by explaining bond strength in terms of d-band center position.

Fig. 11 Pt 5d (a) and O 2p (b) PDOS for a range of oxygen ratios as indicated by the color gradient. The Fermi energy, set to 0 eV, is indicated by the vertical dashed line.

In contrast, the O 2p PDOS does not display substantial variations, with only a slight decrease in the peak located at −6.0 eV and the appearance of states above the Fermi level. The peak centered around −6.0 eV shows a noticeable broadening and shifts towards higher energies with increased oxidation. This shift suggests a weakening of the initial Pt/O 2p interaction, possibly due to increased competition among oxygen adsorbates for available Pt states.

The appearance of the peak on the Pt 5d PDOS possibly signifies a more molecular-like interaction with O 2p states and the formation of a distinct bonding–antibonding pair with these orbitals. The emergence of O 2p states just above the Fermi level with increasing coverage further supports this interpretation, which would likely represent the antibonding molecular-like orbital in this picture. In Fig. 11b, the states lying on the Fermi energy, increasing with coverage, are harder to interpret. Considering the average destabilization of the Pt–O bond with increasing oxidation, we suggest that these states might be non-bonding or weakly antibonding states. In opposition with the first argument, the slight increase in bandwidth for the Pt 5d PDOS indicates that the Pt–O interaction should retain metallic character. This suggests that the Pt–O interaction is not purely ionic, but rather a metallic bond with a significant ionic character.

For the O 2p PDOS the discussed variations mainly occur at low oxygen ratios (0.25 to 0.6) but vary less at higher ratios. In contrast, changes in Pt 5d PDOS are more gradual, suggesting that the Pt structure continuously accommodate to additional oxide species. This indicates oxide species in fixed electronic states for higher oxygen ratios, which can also explain the initial linear regime at low oxygen ratios followed by the plateauing trend at higher oxygen ratios in Fig. 8.

Interestingly, our Pt 5d PDOS trends cannot be linked to both previous theoretical literature and experimental evidence, as both seem to predict inverted trends, i.e. a decrease of intensity close to the Fermi energy with a diffuse increase in intensity for lower states. A DFT study of Verga et al.³⁰ investigated various oxygen coverages using DFT on perfect nanoparticles, and showed that the density of states clearly shift to lower energies as the oxygen ratio increases, without the presence of a narrow peak below the Fermi energy. This trend is also reported for other metals, such as Pd, Rh, and Cu.^77,78 Similarly, Takagi et al.²⁸ reported the same behavior experimentally on similarly sized Pt nanoparticles using HAXPES. However, in these studies, nanoparticles are only partially oxidized, e.g. oxygen atoms are only located on the outer shell. This strengthens our conclusion that the features observed in our PDOS are due to a progression towards complete oxidation of the nanoparticle. An additional reason for this discrepancy could be the peculiar geometries given by the ReaxFF forcefield, not entirely addressed even after a MACE-MP-0 geometry optimization.

4 Discussions

We investigated the oxidation behavior of a realistic platinum nanoparticle using a multistep approach combining ReaxFF and MACE-MP-0 forcefields, complemented by DFT calculations. Our hybrid MD-MC simulations revealed that the nanoparticle becomes highly oxidized at pressures above 10⁻⁵ atm, with coverages that are in reasonable agreement with existing literature. We performed multiple XRD, TEM, and EXAFS measurements on samples of different sizes to confirm the plausibility of such oxidation. For our lowest platinum loading sample (10 wt%, 1.7 nm average diameter), we observed predominantly oxidized species with no metallic character from our XRD measurements. Similarly, for the larger sample (20 wt%, 2.4 nm average diameter) XRD did not show a significant crystallite domain. However, our EXAFS analysis showed a low Pt–Pt coordination number from which we conclude that the particles are in a mixed state of metallic and oxidized platinum. These findings, along with the existing literature suggest that the 1.5 to 2.5 nm size range represents a transitional behavior between larger, partially oxidized nanoparticles and smaller, fully oxidized nanoparticles. We suggest that such Monte Carlo simulations will consistently predict complete oxidation regardless of particle sizes due to the lack of kinetic aspects. This aspect is crucial in experimental conditions, where the kinetic associated with oxygen diffusion into the platinum lattice hinders the oxidation process. It is then desirable to develop computational methods that can balance between kinetic and thermodynamic aspects to provide more relevant results for catalytic conditions.

Our simulations revealed the presence of Pt₆O₈ clusters previously reported by theoretical studies using similar ReaxFF forcefields. While these clusters emerged as potential catalytic degradation products and prompted further dynamics, subsequent MACE-MP-0 forcefield optimizations showed systematic reincorporation of these species into the nanoparticle. Comparison with DFT calculations supported MACE-MP-0 reliability, and, along with the lack of experimental evidence for such species, suggested that the Pt₆O₈ clusters were likely artifacts from the ReaxFF forcefield.

We further explored the thermodynamic properties of the oxidized nanoparticles by calculating the oxygen binding energy using the CHE scheme. By computing most stable configurations in potential intervals, we were able to retrieve the expected place-exchange mechanism observed experimentally at ≈1.05 V. In our simulations this phenomenon is directly caused by the plateauing behavior of the binding energy per atom vs. oxygen ratio curve predicted by ReaxFF. Interestingly, despite the inaccuracies of the ReaxFF forcefield, we found many points of agreement with the experimental and theoretical literature, such as the onset potential for oxidation and the place-exchange mechanism. This suggests that such properties do not require accurate capture of the thermodynamics and correct identification of the oxide phase forming at high potential. This hints at an underlying phenomenon, which would be common to various oxide-phases, most likely linked to the electronic structure of the catalyst. We attempted to explore this possibility by computing Pt 5d and O 2p PDOS. We found modifications of the platinum electronic structure, leading to a localized peak close to the Fermi energy interpreted as an increase of the ionic character of the Pt–O bonds. In contrast, we did not observe the expected decrease in bandwidth of the d-band, suggesting that some metallic character is retained, resulting in a hybrid, metallic-oxide-like structure. The former ionic character can possibly be linked to the high onset potential of the ORR on Pt, as it was previously shown that the reaction thermodynamic is altered on oxidized catalyst.

Our study is a first step toward modelling of catalysts under realistic conditions and offer solid foundation to pave the way for more focused studies on the oxidized nanoparticles and their underlying physics. Future work should focus on providing a more accurate picture of the potential energy surface by using fine-tuned foundation models. Additionally, the kinetic aspects of oxidation should be investigated, as the thermodynamic stability of highly oxidized states may not necessarily translate to their practical relevance under catalytic conditions. By combining accurate forcefields with DFT calculations and experimental validation, it may be possible to design more efficient catalysts for the oxygen reduction reaction and other important applications.

Data availability

Data supporting this article, including initial and final snapshots, LAMMPS and ONETEP input/output files are available at the following DOI: https://doi.org/10.5281/zenodo.14630262.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work. We are grateful for access to the ARCHER2 national supercomputer which was obtained via the United Kingdom Car–Parrinello Consortium (UKCP) funded by EPSRC (Grant No. EP/X035956/1). We thank H. Marchbank, A. Phillips, G. Goodlet, and T. Tran for their assistance with experimental measurements. We thank E. Gerouville for proofreading.

Notes and references

1. M. Shao, Q. Chang, J.-P. Dodelet and R. Chenitz, Chem. Rev., 2016, 116, 3594–3657.
2. S. Sui, X. Wang, X. Zhou, Y. Su, S. Riffat and C.-J. Liu, J. Mater. Chem. A, 2017, 5, 1808–1825.
3. J. Stacy, Y. N. Regmi, B. Leonard and M. Fan, Renewable Sustainable Energy Rev., 2017, 69, 401–414.
4. N. Jung, D. Y. Chung, J. Ryu, S. J. Yoo and Y.-E. Sung, Nano Today, 2014, 9, 433–456.
5. K. C. Neyerlin, W. Gu, J. Jorne and H. A. Gasteiger, J. Electrochem. Soc., 2006, 153, A1955.
6. R. Jinnouchi, K. K. T. Suzuki and Y. Morimoto, Catal. Today, 2016, 262, 100–109.
7. Z. Wu, G. Yang, Q. Zhang, Z. Liu and F. Peng, Energy Environ. Sci., 2024, 17, 3338–3346.
8. Y. Li and Z.-F. Liu, Phys. Chem. Chem. Phys., 2020, 22, 22226–22235.
9. A. Kulkarni, S. Siahrostami, A. Patel and J. K. Nørskov, Chem. Rev., 2018, 118, 2302–2312.
10. M. M. Montemore, M. A. van Spronsen, R. J. Madix and C. M. Friend, Chem. Rev., 2018, 118, 2816–2862.
11. J. L. Gland, Surf. Sci., 1980, 93, 487–514.
12. J. L. Gland, B. A. Sexton and G. B. Fisher, Surf. Sci., 1980, 95, 587–602.
13. G. N. Derry and P. N. Ross, Surf. Sci., 1984, 140, 165–180.
14. C. T. Campbell, G. Ertl, H. Kuipers and J. Segner, Surf. Sci., 1981, 107, 220–236.
15. M. Peuckert and H. P. Bonzel, Surf. Sci., 1984, 145, 239–259.
16. D. H. Parker, M. E. Bartram and B. E. Koel, Surf. Sci., 1989, 217, 489–510.
17. S. P. Devarajan, J. A. Hinojosa and J. F. Weaver, Surf. Sci., 2008, 602, 3116–3124.
18. D. Fantauzzi, S. Krick Calderón, J. E. Mueller, M. Grabau, C. Papp, H.-P. Steinrück, T. P. Senftle, A. C. T. van Duin and T. Jacob, Angew. Chem., Int. Ed., 2017, 56, 2594–2598.
19. Z. Zhu, F. F. Tao, F. Zheng, R. Chang, Y. Li, L. Heinke, Z. Liu, M. Salmeron and G. A. Somorjai, Nano Lett., 2012, 12, 1491–1497.
20. C. Ellinger, A. Stierle, I. K. Robinson, A. Nefedov and H. Dosch, J. Phys.: Condens. Matter, 2008, 20, 184013.
21. J. Chen, B. Lim, E. P. Lee and Y. Xia, Nano Today, 2009, 4, 81–95.
22. N. Seriani and F. Mittendorfer, J. Phys.: Condens. Matter, 2008, 20, 184023.
23. H. Wang, Y. Wang, Z. Zhu, A. Sapi, K. An, G. Kennedy, W. D. Michalak and G. A. Somorjai, Nano Lett., 2013, 13, 2976–2979.
24. F. Şen and G. Gökağaç, J. Phys. Chem. C, 2007, 111, 5715–5720.
25. V. A. Saveleva, V. Papaefthimiou, M. K. Daletou, W. H. Doh, C. Ulhaq-Bouillet, M. Diebold, S. Zafeiratos and E. R. Savinova, J. Phys. Chem. C, 2016, 120, 15930–15940.
26. H. Imai, K. Izumi, M. Matsumoto, Y. Kubo, K. Kato and Y. Imai, J. Am. Chem. Soc., 2009, 131, 6293–6300.
27. R. Banerjee, D. A. Chen, S. Karakalos, M.-L. C. Piedboeuf, N. Job and J. R. Regalbuto, ACS Appl. Nano Mater., 2018, 1, 5876–5884.
28. Y. Takagi, H. Wang, Y. Uemura, T. Nakamura, L. Yu, O. Sekizawa, T. Uruga, M. Tada, G. Samjeské, Y. Iwasawa and T. Yokoyama, Phys. Chem. Chem. Phys., 2017, 19, 6013–6021.
29. H. Yoshida, H. Omote and S. Takeda, Nanoscale, 2014, 6, 13113–13118.
30. L. G. Verga, J. Aarons, M. Sarwar, D. Thompsett, A. E. Russell and C.-K. Skylaris, Faraday Discuss., 2018, 208, 497–522.
31. B. Kirchhoff, L. Braunwarth, C. Jung, H. Jónsson, D. Fantauzzi and T. Jacob, Small, 2020, 16, 1905159.
32. B. Kirchhoff, C. Jung, H. Jónsson, D. Fantauzzi and T. Jacob, J. Phys. Chem. C, 2022, 126, 6773–6781.
33. R. Slapikas, I. Dabo and S. B. Sinnott, Comput. Mater. Sci., 2022, 209, 111364.
34. D.-H. Lim and J. Wilcox, J. Phys. Chem. C, 2011, 115, 22742–22747.
35. J. Aarons, L. Jones, A. Varambhia, K. E. MacArthur, D. Ozkaya, M. Sarwar, C.-K. Skylaris and P. D. Nellist, Nano Lett., 2017, 17, 4003–4012.
36. T. Ellaby, J. Aarons, A. Varambhia, L. Jones, P. Nellist, D. Ozkaya, M. Sarwar, D. Thompsett and C.-K. Skylaris, J. Phys.: Condens. Matter, 2018, 30, 155301.
37. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
38. I. Batatia, P. Benner, Y. Chiang, A. M. Elena, D. P. Kovács, J. Riebesell, X. R. Advincula, M. Asta, W. J. Baldwin, N. Bernstein, A. Bhowmik, S. M. Blau, V. Cărare, J. P. Darby, S. De, F. Della Pia, V. L. Deringer, R. Elijošius, Z. El-Machachi, E. Fako, A. C. Ferrari, A. Genreith-Schriever, J. George, R. E. A. Goodall, C. P. Grey, S. Han, W. Handley, H. H. Heenen, K. Hermansson, C. Holm, J. Jaafar, S. Hofmann, K. S. Jakob, H. Jung, V. Kapil, A. D. Kaplan, N. Karimitari, N. Kroupa, J. Kullgren, M. C. Kuner, D. Kuryla, G. Liepuoniute, J. T. Margraf, I.-B. Magdău, A. Michaelides, J. H. Moore, A. A. Naik, S. P. Niblett, S. W. Norwood, N. O'Neill, C. Ortner, K. A. Persson, K. Reuter, A. S. Rosen, L. L. Schaaf, C. Schran, E. Sivonxay, T. K. Stenczel, V. Svahn, C. Sutton, C. van der Oord, E. Varga-Umbrich, T. Vegge, M. Vondrák, Y. Wang, W. C. Witt, F. Zills and G. Csányi, A Foundation Model for Atomistic Materials Chemistry, 2023.
39. H. M. Aktulga, J. C. Fogarty, S. A. Pandit and A. Y. Grama, Parallel Comput., 2012, 38, 245–259.
40. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
41. L. Gai, Y. K. Shin, M. Raju, A. C. T. Van Duin and S. Raman, J. Phys. Chem. C, 2016, 120, 9780–9793.
42. D. Fantauzzi, J. Bandlow, L. Sabo, J. E. Mueller, A. C. T. van Duin and T. Jacob, Phys. Chem. Chem. Phys., 2014, 16, 23118–23133.
43. T. C. Allison, NIST-JANAF Thermochemical Tables – SRD 13, 2013.
44. I. Batatia, D. P. Kovács, G. N. C. Simm, C. Ortner and G. Csányi, MACE: Higher Order Equivariant Message Passing Neural Networks for Fast and Accurate Force Fields, 2023.
45. J. Riebesell, R. E. A. Goodall, P. Benner, Y. Chiang, B. Deng, A. A. Lee, A. Jain and K. A. Persson, Matbench Discovery – A Framework to Evaluate Machine Learning Crystal Stability Predictions, 2024.
46. J. C. A. Prentice, J. Aarons, J. C. Womack, A. E. A. Allen, L. Andrinopoulos, L. Anton, R. A. Bell, A. Bhandari, G. A. Bramley, R. J. Charlton, R. J. Clements, D. J. Cole, G. Constantinescu, F. Corsetti, S. M.-M. Dubois, K. K. B. Duff, J. M. Escartín, A. Greco, Q. Hill, L. P. Lee, E. Linscott, D. D. O'Regan, M. J. S. Phipps, L. E. Ratcliff, Á. R. Serrano, E. W. Tait, G. Teobaldi, V. Vitale, N. Yeung, T. J. Zuehlsdorff, J. Dziedzic, P. D. Haynes, N. D. M. Hine, A. A. Mostofi, M. C. Payne and C.-K. Skylaris, J. Chem. Phys., 2020, 152, 174111.
47. Á. Ruiz-Serrano and C.-K. Skylaris, J. Chem. Phys., 2013, 139, 054107.
48. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77(18), 3865–3868.
49. P. C. Jennings, H. A. Aleksandrov, K. M. Neyman and R. L. Johnston, Phys. Chem. Chem. Phys., 2014, 16, 26539–26545.
50. A. G. Yohannes, K. Fink and I. Kondov, Nanoscale Adv., 2022, 4, 4554–4569.
51. P. Légaré, Surf. Sci., 2005, 580, 137–144.
52. Z. Gu and P. B. Balbuena, J. Phys. Chem. C, 2007, 111, 9877–9883.
53. J. M. Hawkins, J. F. Weaver and A. Asthagiri, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 79, 125434.
54. R. B. Getman, W. F. Schneider, A. D. Smeltz, W. N. Delgass and F. H. Ribeiro, Phys. Rev. Lett., 2009, 102, 076101.
55. F. Calle-Vallejo, Adv. Sci., 2023, 10, 2207644.
56. M. Tada, S. Murata, T. Asakoka, K. Hiroshima, K. Okumura, H. Tanida, T. Uruga, H. Nakanishi, S.-I. Matsumoto, Y. Inada, M. Nomura and Y. Iwasawa, Angew. Chem., Int. Ed., 2007, 46, 4310–4315.
57. R. K. Nomiyama, M. J. Piotrowski and J. L. F. Da Silva, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 84, 100101.
58. N. Seriani, W. Pompe and L. C. Ciacchi, J. Phys. Chem. B, 2006, 110, 14860–14869.
59. D. E. Bugaris, M. D. Smith and H.-C. zur Loye, Inorg. Chem., 2013, 52, 3836–3844.
60. A. F. Zurhelle, W. Stehling, R. Waser, R. A. De Souza and S. Menzel, Adv. Mater. Interfaces, 2022, 9, 2101257.
61. S. Takao, O. Sekizawa, G. Samjeské, T. Kaneko, K. Higashi, Y. Yoshida, X. Zhao, T. Sakata, T. Yamamoto, T. Gunji, T. Uruga and Y. Iwasawa, ACS Appl. Mater. Interfaces, 2018, 10, 27734–27744.
62. T. Fuchs, J. Drnec, F. Calle-Vallejo, N. Stubb, D. J. S. Sandbeck, M. Ruge, S. Cherevko, D. A. Harrington and O. M. Magnussen, Nat. Catal., 2020, 3, 754–761.
63. J. K. Nørskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard and H. Jónsson, J. Phys. Chem. B, 2004, 108, 17886–17892.
64. E. Sargeant, F. Illas, P. Rodríguez and F. Calle-Vallejo, J. Electroanal. Chem., 2021, 896, 115178.
65. D. Alfonso, D. Tafen and D. Kauffmann, Catalysts, 2018, 8, 424.
66. A. Farkas, D. Fantauzzi, J. E. Mueller, T. Zhu, C. Papp, H.-P. Steinrück and T. Jacob, J. Electron Spectrosc. Relat. Phenom., 2017, 221, 44–57.
67. J. Huang, J. Zhang and M. Eikerling, Phys. Chem. Chem. Phys., 2018, 20, 11776–11786.
68. M. Wakisaka, H. Suzuki, S. Mitsui, H. Uchida and M. Watanabe, Langmuir, 2009, 25, 1897–1900.
69. H. S. Casalongue, S. Kaya, V. Viswanathan, D. J. Miller, D. Friebel, H. A. Hansen, J. K. Nørskov, A. Nilsson and H. Ogasawara, Nat. Commun., 2013, 4, 2817.
70. W. Chen, J. Huang, J. Wei, D. Zhou, J. Cai, Z.-D. He and Y.-X. Chen, Electrochem. Commun., 2018, 96, 71–76.
71. H. You, D. J. Zurawski, Z. Nagy and R. M. Yonco, J. Chem. Phys., 1994, 100, 4699–4702.
72. H. You and Z. Nagy, Phys. B, 1994, 198, 187–194.
73. M. Teliska, W. E. O'Grady and D. E. Ramaker, J. Phys. Chem. B, 2005, 109, 8076–8084.
74. T. Zhu, E. J. M. Hensen, R. A. Van Santen, N. Tian, S.-G. Sun, P. Kaghazchi and T. Jacob, Phys. Chem. Chem. Phys., 2013, 15, 2268.
75. N. G. Hörmann, O. Andreussi and N. Marzari, J. Chem. Phys., 2019, 150, 041730.
76. R. Guidelli and W. Schmickler, Modern Aspects of Electrochemistry, Springer US, Boston, MA, 2005, vol. 38, pp. 303–371.
77. N. İnoğlu and J. R. Kitchin, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 82, 045414.
78. J. R. Kitchin, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 79, 205412.

---

Footnote

† Electronic supplementary information (ESI) available: Detailed description of the parameters used for both forcefields, Monte Carlo, molecular dynamics, and DFT calculations. Description of the experimental methods used for the synthesis and characterization of the nanoparticles. See DOI: https://doi.org/10.1039/d5cp00134j

---