Direct correlation of oxygen adsorption on platinum-electrolyte interfaces with the activity in the oxygen reduction reaction

See allHide authors and affiliations

Science Advances  09 Jun 2021:
Vol. 7, no. 24, eabb1435
DOI: 10.1126/sciadv.abb1435


The oxygen reduction reaction (ORR) on platinum catalysts is essential in fuel cells. Quantitative predictions of the relative ORR activity in experiments, in the range of 1 to 50 times, have remained challenging because of incomplete mechanistic understanding and lack of computational tools to account for the associated small differences in activation energies (<2.3 kilocalories per mole). Using highly accurate molecular dynamics (MD) simulation with the Interface force field (0.1 kilocalories per mole), we elucidated the mechanism of adsorption of molecular oxygen on regular and irregular platinum surfaces and nanostructures, followed by local density functional theory (DFT) calculations. The relative ORR activity is determined by oxygen access to platinum surfaces, which greatly depends on specific water adlayers, while electron transfer occurs at a similar slow rate. The MD methods facilitate quantitative predictions of relative ORR activities of any platinum nanostructures, are applicable to other catalysts, and enable effective MD/DFT approaches.


The oxygen reduction reaction (ORR) has many applications including hydrogen fuel cells for energy conversion (1, 2), biological respiration (3), and materials dissolution (4). Platinum catalysts and platinum alloy catalysts (57) have shown promising catalytic activity for ORR and are widely used to increase slow reaction rates and address high overpotential (8, 9). At the same time, the high cost of platinum requires improvements in catalyst composition and performance. A major problem for designing ORR catalysts are reports of conflicting reaction mechanisms (1017), difficulties to rationalize reported studies for the development of catalysts with desired specifications, and lack of predictions of reaction rates. Therefore, trial-and-error search in experiments continues to be the main mode of discovery. A major roadblock has been coarse approximations in prior theoretical studies, especially the neglect of solution conditions (1428). Oxygen has a very low O2 solubility in water under standard conditions (0.00025 M) (29), which greatly limits contact with the surface; the acid concentration is high (0.1 M); and solvent layers on the surface of Pt catalyst particles play a critical role for ORR activity. This study illuminates the mechanisms and predictions of initial oxygen adsorption, which is the first step of ORR, using molecular dynamics (MD) simulations of the electrode-electrolyte interfaces in high accuracy. We identified a direct correlation of oxygen access to Pt surfaces with ORR activity, which augments our understanding of the ORR process and enables predictions of relative activities between different Pt catalysts. We explain the use of computational MD methods to support the design of Pt catalysts for ORR and similar catalysts and electrode materials from the atomic scale.

The ORR is a multistep, four-electron reduction with several intermediates, O2 (aq) + 4 H3O+ (aq) + 4e → 6 H2O (shown in acid environment) (fig. S1, A and B) (30). The initial O2 adsorption to the catalyst surface in solution has not been studied before since suitable experimental and computational techniques have not been available. The subsequent elementary reaction steps and possible intermediates have been extensively explored using experimental and computational approaches. In situ/operando techniques such as synchrotron x-ray spectroscopy recently made large advances and support the existence of bonded oxygen intermediates (31). The precise chemistry and lifetime of the intermediates on the surface, however, could not yet be identified. Many quantum mechanical (QM) and theory studies using density functional theory (DFT) have been carried out to gain insights into possible reaction intermediates and energy barriers (1417).

Nevertheless, the role of initial oxygen adsorption has remained unclear. QM calculations commonly neglect explicit solvents, electrolytes, and specific concentrations that play a vital role for ORR to occur (14, 16, 1823). These simplifications are a major challenge since ORR does not occur in vacuum and the use of implicit water models or few water molecules (clusters) in DFT simulations cannot approximate realistic electrolytes and dynamics. Many DFT simulations also assume an unrealistic temperature of 0 K (water would be frozen to ice in the fuel cell), and the time scale of the reaction at the nanometer scale of milliseconds is not accessible (2428). Large uncertainties also arise from the choice of exchange-correlation functionals, including up to 40% error in Pt surface energies and in adsorption energies (table S3) (32). Together, these limitations explain that DFT calculations disregard the initial adsorption-desorption dynamics of O2 and lead to limited understanding of the entire ORR process.

As a result, existing studies diverge with regard to the ORR mechanism and several conceivable pathways, and rate-determining steps have been reported (33). For example, the first electron transfer from Pt to O2 and the dissociation of O─O bonds in peroxide species have been suggested as rate-limiting (slowest) steps, while experimental evidence is still not sufficient to back up either hypothesis (1013). The formation of *OH groups and *OH/*OOH intermediates bound to Pt surfaces appears likely, and binding energies from DFT calculations have been proposed (1417). However, the occurrence of bound species depends on the external voltage, confirmation by experiments is not available under typical operating conditions of 0.9 V versus reference hydrogen electrode (RHE), and the lifetime may be short.

Verifiable data include the range of ORR activities of Pt catalysts in experiments, which is approximately 1 to 30 times relative to Pt/C catalysts (k1k2=30) (5, 34, 35). Considering the Arrhenius equation to describe the reaction kinetics (k1k2=eERT), the difference in the effective activation energy ΔE for ORR is then less than 0.09 eV (= 2.0 kcal/mol = 3.4 RT at 298.15 K). The differences are thus smaller than typical uncertainties in DFT calculations, which are on the order of 0.1 eV (= 2.3 kcal/mol = 4 RT at 298.15 K) (36, 37) and suffer from the additional limitations described above. The reliability of DFT methods is therefore too coarse for quantitative predictions of relative ORR rates even if limitations to include explicit solvents, temperatures, and dynamics could be overcome.

Here, we use a to-date unused tool for ORR analysis, MD simulations with the Interface force field (IFF) that reaches the necessary accuracy of 0.005 eV (= 0.1 kcal/mol = 0.2 RT) for metal-water-gas interfaces and can describe the dynamics of entire nanoscale Pt catalyst particles over tens of nanoseconds. The MD simulations uncover the role of solvent and solute adlayers, free energies of O2 binding, and statistics of O2 adsorption on the metal surfaces. Simulations using IFF recently generated insights into nanocrystal growth, solvation effects, and activity of other catalysts in combination with experimental data and local DFT calculations due to accurate representations of chemical bonding, interpretable parameters, and validated intermolecular energies (23, 3841). We used new force field parameters for Pt/electrolyte/O2 interfaces that exceed the accuracy of prior models by a factor of 10, including several times higher reliability than embedded atom models and density functionals in metal surface energies and interfacial energies (tables S1 to S4 and fig. S1C). We overcome the following major problems: (i) The force field for O2 is compatible with common energy expressions for water and metals, enabling the simulation of mixtures, which was not possible with incompatible prior models; (ii) we provide validation of structural and thermodynamic properties, including critical interactions of O2 gas and Pt surfaces with water, reducing errors or previous model in excess of 100 to <5% (table S2 and section S1); and (iii) the new polarizable Pt parameters capture image charges and responses to electric fields on the fly, which was not possible before (fig. S1C and tables S1, S3, and S4) (36, 42, 43).

The dynamics and adsorption of molecular O2 to the metal surfaces could thus not be studied previously, especially under realistic electrolyte conditions, and deviations higher than 100% in key interfacial properties are reduced to ~5% with the new IFF models. Knowledge of initial oxygen adsorption is essential to examine subsequent ORRs and intermediates (4449). The models and MD simulations here add this critical capability, facilitate the simulation of entire nanoscale catalysts, include more realistic conditions in fuel cells than reported before (16, 17, 5053), and reach high accuracy that is necessary to predict relative activities. Here, oxygen adsorption is only part of the overall ORR process. We explain that O2 adsorption scales with the relative ORR rate, while the slowest (rate-determining) step likely involves initial chemisorption to Pt, O─O dissociation, formation of water, or another subsequent step in the mechanism. The rate of this yet uncertain slowest step appears to be approximately constant and independent of the specific Pt surface structure, such as (hkl) facet, defective, or irregular surface, and renders the relative ORR rate proportional to oxygen adsorption to the Pt surface, which we can predict by MD simulation.


We find that the oxygen affinity to differently structured Pt surfaces, e.g., in the form of free energy curves or contact time, directly correlates with the ORR activity measured in experiments (Fig. 1A). Solvent layers on the Pt surface structures play a critical role as they limit access of molecular oxygen to Pt surface features to a different extent, commensurate with the expected small differences in effective activation energies (<2.3 kcal/mol) that correspond to rate differences up to ~50 times in experiments. The effective contact time with the surface correlates with relative ORR activity, while the slowest reaction steps, i.e., electron transfer or subsequent steps during the formation of water, seem to occur at the same rate independent of the Pt surface geometry. The affinity of oxygen molecules to Pt surfaces increases in the order Pt(100) < Pt(111) < Pt(110) and even further for regular Pt nanowires (R-NWs), Pt nanoplates, and jagged Pt nanowires (J-NWs), which matches the trend of higher ORR activity observed in experiments (Fig. 1B) (5, 34, 35, 5456). This order was observed in the absence of site-blocking electrolytes, for example, in perchloric acid solution. According to experiments, the rate slows down in the presence of hydrogen sulfate or phosphoric acid that inhibit adsorption sites (9, 12, 54, 5658), consistent with the proposed correlation of oxygen adsorption and Pt surface accessibility (table S5).

Fig. 1 Correlation of oxygen adsorption on Pt(hkl) surfaces and Pt nanostructures in acid solution with ORR activity.

(A) Oxygen molecules encounter barriers to Pt surface contact through superficial molecular layers of water. The number of red arrows indicates the affinity, and the depth of blue color above the surfaces represents the affinity of water adlayers to Pt surfaces. Corrugation and irregularity of the Pt surfaces in nanowires and nanoplates cause disruptions in the water layers, locally stronger O2 adsorption, O2 coordination numbers with Pt (CN), and increased ORR activity. (B) Relative order oxygen adsorption in simulations and ORR activity in experiments [data from (5, 9, 12, 34, 35, 5458)]. Further support comes from competition by electrolytes (Fig. 3). RRDE = Rotating ring-disk electrode, HMRDE = hanging meniscus rotating disk electrode.

Oxygen adsorption onto platinum in water is thus surface sensitive and electrolyte sensitive, and small differences in the involved energy barriers are decisive for ORR. First, we describe the interactions of O2 with well-defined low-index Pt(hkl) surfaces in water at the atomic level, including binding free energies and an examination of site-blocking role of electrolytes such as HSO4 ions and H3PO4 using MD simulation. Second, we discuss the role of time scale and applied potentials. Third, we examine possible initial reaction pathways with DFT following adsorption using realistic initial configurations from MD simulations. Fourth, we examine oxygen adsorption and the correlation with ORR activity for Pt nanowire and Pt nanoplate surfaces, which feature more realistic irregular shapes found in ORR catalyst particles, including grain boundaries and other defects. The paper ends with conclusions.

Oxygen affinity to Pt(hkl) surfaces and binding free energies

The free energy profiles of oxygen adsorption from steered molecular dynamics (SMD) simulations indicate adsorption free energies of oxygen molecules in direct contact with Pt(110), Pt(111), and Pt(100) surfaces of −1.54, −1.03, and +0.83 kcal/mol, respectively (Fig. 2, A to C). Negative adsorption free energies of O2 on Pt(110) and Pt(111) indicate thermodynamically favorable and direct contact. On the contrary, oxygen molecules avoid contact with the Pt(100) surface atomic layer, although a preference for oxygen adsorption is seen within the second adlayer of water. This free energy minimum of O2 on Pt(100) amounts to −0.73 kcal/mol, and adsorbed O2 molecules can still directly approach Pt atoms in the top atomic layer to initiate electron transfer, although likely less effectively. Accordingly, the differences in binding free energies of molecular oxygen are small, on the order of 0.81 kcal/mol, and similar to RT (1 RT = 0.59 kcal/mol = 0.026 eV at 298.15 K). These differences imply subtle differences in surface contact time by a factor of 2 to 4 using equilibrium constants (K=eGRT) and semiquantitatively match the changes in ORR activity by a factor of 2 to 3 observed in experiment (Fig. 1B) (5). The magnitude is comparable to thermal motion and statistically significant over time, similar to other equilibrium and kinetic processes. For example, rotation states in alkanes and proteins differ by 0.6 kcal/mol between anti and gauche conformations with well-known implications on folding and reactivity (59, 60). If more energy-intensive changes of bond-breaking and bond formation processes would determine the relative ORR rate, as suggested in earlier theoretical studies (1417), then we would expect much larger differences in ORR activity. For example, a difference of 2.3 or 23 kcal/mol (0.1 or 1 eV) in barriers would result in 50 or 1017 times different rate, respectively, while ORR experiments only show rate differences by two to three times on low-index (hkl) surfaces (according to k1k2=eERT).

Fig. 2 Oxygen adsorption free energy and competition with water on Pt(hkl) surfaces.

(A to C) Free energy profiles of oxygen interaction with Pt(110), Pt(111), and Pt(100) surfaces in solution. The distance for direct contact with the surface is indicated by black triangles and varies slightly for each surface. (D) Water density profiles on Pt(110), Pt(111), and Pt(100) surfaces. (E to G) Representative snapshots of water-oxygen coadsorption on Pt(110), Pt(111), and Pt(100) surfaces in top view and side view. The adsorbed conformation on the Pt(100) surface (G) is thermodynamically unfavorable (seen infrequently) and shown for representation purposes only. Water molecules are shown as line representation in side view and CPK representation in top view for adlayer. Oxygen molecules are shown in VDW representation as defined in the VMD program.

We also examined adsorption energies of O2 on the Pt surfaces in vacuum that are often used for simplicity in QM calculations (fig. S1D). Adsorption energies between −7 and −8 kcal/mol are many times stronger than in solution, whereby vacuum allows uninhibited contact of O2 with the Pt surfaces. Therefore, rate-modifying energy barriers by water adlayers and electrolytes are excluded, and simulations in vacuum add little or no insight related to ORR activity.

The free energy barriers of oxygen access to the Pt surface in aqueous solution are associated with networks of hydrogen bonds that must be broken (Fig. 2, A to C). Hydrogen bonding between superficial water molecules thereby depends on the (hkl) surface as can be seen by the different location of the peaks in the water density profiles and the respective intensities (Fig. 2D). The location and height of these peaks correlate with the barriers for O2 adsorption and regulate O2 penetration (Fig. 2, A to C). The intrinsic roughness of the Pt(110) surface promotes corrugation of the adlayer of water molecules (Fig. 2E) and enables many hydrogen bonds between the first water adlayer and bulk water molecules. The higher adlayer mobility on the corrugated Pt(110) surface results in continuity of the water density profile past the adlayer and eases exchange with oxygen molecules (Fig. 2D). The free energy minimum is found at 2.1 Å (mogul in Fig. 2A), and O2 molecules break the H-bonded network before replacing adsorbed water molecules (Fig. 2, A and E). Conversely, on the even Pt(111) and Pt(100) surfaces, the flat-on structured water adlayer is somewhat disjointed from bulk water and contains a stronger network of hydrogen bonds (Fig. 2, F and G). When the relatively bulky O2 molecules enter the low-density region around 4 Å (Fig. 2D), higher free energy barriers arise from simultaneous breaking of hydrogen bonds and replacement of water molecules. Oxygen adsorption on Pt(100) involves a particularly high energy barrier associated with penetration into a near-perfectly ordered network of hydrogen bonds in the water adlayer (Fig. 2G) (38). The Pt(100) water adlayer is strongly stabilized by a square pattern of epitaxial sites with a side length of 2.78 Å, which attracts the oxygen atoms in water molecules with O···O distances very close to those found in bulk water (~2.8 Å) and involves four shared hydrogen bonds per water molecule (6163). This epitaxial pattern and preferential alignment of O─H bonds on bridge sites explain the high barrier for O2 adsorption (4 kcal/mol) and unfavorable close contact overall (Fig. 2G). Predictions of the structure and thermodynamics of metal-electrolyte-gas interfaces in a high accuracy of 0.1 kcal/mol require the quality of IFF models, as previously shown for metal-water and metal-organic interfaces (23, 39, 64, 65).

Influence of competing electrolytes

The correlation of O2 adsorption to Pt with ORR rates is further supported by experimental data for competing adsorption of other electrolytes (Fig. 3, A and B). While we focus on aqueous solutions and perchloric acid solutions that are believed to involve minor coadsorption, experiments have shown that ORR activity decreases on all low-index Pt surfaces in sulfuric acid solution and in phosphoric acid solution at pH ≤ 1 (table S5). Site-blocking effects are known to be particularly strong on Pt(111) surfaces (9, 12, 54, 5658). MD simulations of the affinity of HSO4 ions and H3PO4 molecules on the low-index Pt(hkl) surfaces confirm these trends at the atomic scale (Fig. 3, C to E). SMD simulations show strong direct adsorption of HSO4 ions with free energies of binding of −8.0, −3.6, and −5.3 kcal/mol, respectively, on Pt(111), Pt(100), and Pt(110) surfaces in water. H3PO4 behaves similar with free energies of binding of −7.3, −2.9, and −5.4 kcal/mol, respectively, and has largely the same effect on oxygen adsorption as HSO4 ions according to experiments and simulations (Fig. 3, C to E). The two electrolytes thus adsorb multiple times stronger to Pt surfaces than oxygen molecules (3 to 12 RT more), are typically present in 400 times higher concentration (100 mM versus 0.25 mM), and inhibit Pt(111) surfaces in particular.

Fig. 3 Influence of electrolytes on ORR activity and underlying mechanism.

(A) Effect of bisulfate ions (HSO4) and phosphoric acid (H3PO4) in acidic solution (pH 1) on overall oxygen adsorption and ORR rate. Strong binding to Pt(111) surfaces strongly diminishes oxygen adsorption and leads to the lowest ORR rate. (B) Relative order of ORR activity according to experiment and simulation. Perchloric acid is known to have only minor competing effects on adsorption and ORR activity and was not explicitly considered. HSO4 and H3PO4 modify the relative availability of (100) and (111) surfaces for ORR compared to perchloric acid solution, resulting in the activity order (111) < (100) < (110) in agreement with multiple experimental studies (9, 12, 54, 5658). The (110) surface maintains the highest relative activity. (C to E) Free energy profiles and representative snapshots of HSO4 ion and H3PO4 adsorption on Pt(111), Pt(100), and Pt(110) surfaces in solution (blue subpanels, bisulfate; red subpanels, phosphoric acid). Preferred distances, measured for the S and P atoms, and associated adsorption free energies are highlighted by black triangles. Tripod-like adsorption geometries on the Pt(111) surface explain the particularly large negative adsorption energies and inhibit oxygen adsorption more than on other (hkl) surfaces. The electrolytes and Pt atoms are shown in VDW representation. Water molecules are shown in line representation.

Adsorption of HSO4 ions and H3PO4 molecules increases in the order Pt(100) < Pt(110) < Pt(111), and the energy barriers to reach the surfaces are smaller compared to O2 molecules (Fig. 3, C to E). The high internal polarity and presence of hydrophilic ─OH groups permit easier integration into the superficial water layers on Pt. The ORR activity in sulfuric acid and phosphoric acid at pH values equal to or lower than 1 is then modified to (111) < (100) < (110) from (100) < (111) < (110) in water or perchloric acid (Fig. 3B and table S5). The results support the proposed correlation of O2 adsorption with ORR activity, explaining how inhibition of adsorption sites on the Pt surface by other species decreases the ORR rate via reduced contact time [especially on Pt(111)]. Summarizing so far, the simulation data explain the relative ORR activity on Pt(hkl) surfaces in four experimental studies in water/HClO4 (Fig. 2) (5, 5456) and in six experimental studies in HSO4 and H3PO4 solution due to electrolyte competition (Fig. 3) (9, 12, 54, 5658). We will discuss additional evidence for irregular Pt nanostructures (Fig. 1) (34, 35, 66) and the role of oxygen partial pressure (67) further below, after in-depth discussion of site-specific adsorption, time scale, applied potential, and initial chemical reactions following adsorption.

Local site-specific adsorption

The analysis of O2 interaction with regular Pt(hkl) surfaces shows that O2 adsorption occurs on specific local sites of the (hkl) surfaces (Fig. 4, A to C). Preferred adsorption sites include epitaxial and threefold sites on (110) surfaces (Fig. 4, A and D) and an even distribution of binding sites without a strong preference on (111) surfaces (Fig. 4, B and E). (100) surfaces favor desorption (Fig. 2C), and intermittent contacts occur on epitaxial and bridge sites (Fig. 4, C and F). Adsorption in the second water layer on Pt(100) surfaces allows proximity to top sites (Fig. 2G) and involves higher surface mobility than in direct contact. In addition to local site preferences on every (hkl) surface, the linearly shaped oxygen molecules display preferred orientation angles with respect to the Pt(hkl) surfaces (Fig. 4, D to F). These orientations could play a role in the onset of electron transfer in ORR as discussed below. O2 molecules exhibit a tilt angle near 35° on both epitaxial and threefold sites on the most-active Pt(110) surfaces (Fig. 4D) and relatively flat orientations <15° on less-active Pt(111) surfaces (Fig. 4E). More tilted orientations of about 40° are seen on Pt(100) surfaces when in contact (Fig. 4F).

Fig. 4 Quantitative analysis of preferred local oxygen adsorption sites and molecular orientation on the Pt surfaces.

(A to C) Scatter plot of the position of oxygen molecules on Pt(110), Pt(111), and Pt(100) surfaces. The points represent adsorption sites, and the color indicates the adsorption angle θ, defined between the surface plane and orientation of the O═O bond (schematic and color bar on the right). (D to F) Local site preferences for adsorbed oxygen molecules on Pt(110), Pt(111), and Pt(100) surfaces, given as a percentage of time adsorbed for each site. On the right-hand side vertical axis of each chart, the average tilt angle of adsorbed oxygen molecules is shown for each common site on the Pt(110), Pt(111), and Pt(100) surfaces. Schematics of the orientation are indicated above each column (upward arrows indicate desorption).

Role of time scale and applied potential for oxygen adsorption and ORR activity

In addition to electrolyte interfaces, it is critical to consider the time scale of the reaction to explain the ORR mechanism and relative activities. We use experimental reference data for reliable benchmarks and MD and DFT results to add details at shorter time scales. The time for Pt oxide layer formation in air is minutes to hours (6870), and the typical current density in ORR measurements in solution (6 mA/cm2) at 0.9 V versus RHE corresponds to milliseconds to observe a reaction event on the nanometer scale (= number of electrons transferred per unit time and surface area) (Fig. 5, A and B) (71). The oxidation reaction in air is slow although O2 is abundantly available and adsorption to Pt is strong (fig. S1D). The ORR reaction in a rotating disk electrode (RDE) is orders of magnitude faster. The applied negative potential and high solution acidity (0.1 M) accelerate the reaction, although O2 availability to the Pt surface is much more limited in solution due to the very low solubility (0.00025 M) and weaker binding to Pt (Fig. 2, A to C). Nevertheless, the time interval for a single electron transfer in ORR is 2.7 ms/nm2 of Pt surface area (Fig. 5B and see calculation details in section S2) (71), which remains orders of magnitude slower compared to adsorption-desorption equilibria on the nanosecond scale in MD simulation (Fig. 2, A to C). Nanosecond time scales for adsorption are also typical for other small molecules (23, 3841). While we use high O2 supersaturations in the MD simulations (0.05 to 0.36 M) to obtain good statistics and to avoid working with excessively large model systems, adsorption events at 0.00025 M O2 concentration would occur at a rate of microseconds per square nanometers of Pt surface area. Therefore, adsorption-desorption equilibria are approximately three orders of magnitude faster than the electron transfer rate in ORR. When we consider further that surface contact is necessary to initiate an ORR reaction, limited O2 contact with the surface and a low probability for electron transfer are the limiting factors for ORR activity. DFT calculations also help to estimate time scales. Large energy barriers for changes in chemical bonding (>> 0.1 eV) suggest electron transfer to be distinctly slower than adsorption-desorption equilibria (1417), consistent with the experimental data.

Fig. 5 Oxygen adsorption under an electrical potential and analysis of the transition from adsorption to chemisorption on Pt(hkl) surfaces.

(A) Setup of a three-electrode system for ORR activity measurements by cyclic voltammetry. (B) Schematic of the electrode configuration in the simulation with electrode colors corresponding to (A). The reference electrode was not included in MD simulations under applied potentials. Voltage, strength of the electric field (E), charge density (σ), and average time for electron transfer (te) are indicated. (C) Free energy profiles of oxygen molecule adsorption on Pt(110) at different electrical potentials from SMD simulation using the polarizable Pt models, as well as nonpolarizable Pt models at zero potential. (D and E) Initial physisorbed (left) and corresponding chemisorbed structures (right) on Pt(110) surfaces according to DFT calculations (optPBE functional), starting with the threefold and epitaxial adsorption sites from MD simulations. The energy difference for chemisorption for the threefold site was set to zero as a reference, and the dissociative reaction on the epitaxial site (E) leads to products that are −1.27 eV lower in energy. (F to I) Initial physisorbed and corresponding chemisorbed structures on Pt(111) from DFT calculations, starting with four equilibrium configurations from MD simulations. (J and K) Initial physisorbed and corresponding chemisorbed structures on Pt(100) surfaces for two possible initial states. Initial adsorption sites and corresponding angles are listed on the bottom left of each panel. Insets at the top of each panel show top views before and after initial reaction (water molecules are hidden for clarity). (L) Summary of the transition from oxygen adsorption to the Pt surface and formation of Pt─O bonds. Oxygen molecules are shown in blue color. Chemisorbed structures with a maximum of four Pt─O bonds on the Pt(110) surface tend to be most stable.

To examine how an external electric potential affects the dynamics of O2 adsorption, we carried out SMD simulations of O2 adsorption with a polarizable Pt model (tables S1, S3, and S4, and fig. S1C). We applied various external potentials, mimicking the three-electrode system used in RDE measurements (Fig. 5, A to C). The onset potential is typically reported between the working electrode and an RHE (72) so that the exact potential between the working electrode and the counter electrode is not known. The potential difference does not exceed a few volts, however (Fig. 5B). For a potential difference of 1 V, we note that the induced surface charge density σ on the two electrodes is only ±10−7 electrons (e)/nm2 (σ = q/A = Eε = Uε/d according to the law of Gauss and the definition of voltage; see section S3). A charge density this small has a negligible effect on oxygen adsorption (Fig. 5C) and requires thousands of encounters with oxygen molecules to induce a single reduction reaction, qualitatively consistent with the experimental rate (1e per 2.7 ms and nm2; Fig. 5B). A visible enhancement of the adsorption free energies of molecular oxygen from −1.54 to −4.50 kcal/mol on the Pt(110) surface requires the potential difference to increase to the megavolt range, then enabling a surface charge density of ~−1 e/nm2, or ~−0.1e per surface atom (Fig. 5C and figs. S1E and S2). External potentials of this high magnitude are not used in fuel cells and would interfere with the standard reduction potential of O2. The results clearly show that typical potentials on the order of 1 V have no impact on oxygen adsorption. Details and full validation of the polarizable IFF for platinum are available in the Supplementary Materials (section S1; tables S1, S3, and S4; and fig. S1C). The parameterization followed the protocol for polarizable Au (39).

In conclusion, the gap in time scales between fast oxygen adsorption/desorption cycles and slow electron transfer rate explains the role of adsorption-desorption equilibria preceding Pt oxidation. The adsorption/desorption dynamics of O2 controls surface access and modulates the otherwise slow and constant activity. Adsorption and desorption of O2 are unaffected by applied voltages in the typical operating range.

ORR mechanism and transition from adsorption to first electron transfer

The applied negative electrode potential thus affects the rate of the slowest reaction steps in ORR only after adsorption occurs. Experiments and simulation data suggest that (i) the contact time of O2 with the Pt surface scales with the relative ORR rate, (ii) initial chemisorption (or somewhere else along the reaction coordinate) is the slowest step, (iii) and the slowest step has about the same barrier regardless of the Pt(hkl) surface environment and electrolyte. Once an electron transfer takes place upon a random adsorption event every few milliseconds, subsequent reactions are likely faster because of the negative potential (below the O2/H2O standard potential) and 103 times higher hydronium concentration (0.1 M) in comparison to the oxygen saturation concentration (0.00025 M). The negative electrode potential adds a driving force to attract hydronium ions to negatively charged Pt-oxygen species (such as Pt─O2) and restore partially oxidized Pt atoms to neutral Pt metal. Difficulties to fully identify intermediates in experiments due to short lifetime support this view (73).

To gain further understanding, we examined the transition from adsorption to chemisorption using electronic structure calculations with the optPBE functional and with simpler DFT functionals (PBE) (Fig. 5, D to L). We used realistic start conformations from high-accuracy MD simulations, different from prior approaches, by including the local solvent environment. We selected several preferred adsorbed conformations of O2 on each Pt(hkl) surface (Fig. 4, D to F) with typical orientation angles θ and several water molecules near the adsorption site. Specifically, we chose the epitaxial and threefold sites on the Pt(110) surface; the top, face-centered cubic (fcc), hexagonal close-packed (hcp), and bridge sites on the Pt(111) surface; and the epitaxial and bridge sites of intermittent contact with the Pt(100) surface (Figs. 4, D to F, and 5, D to K). Oxygen molecules transition from nontop sites, preferred for adsorption, to the nearest top Pt atoms to initiate chemical bonding on all surfaces (Fig. 5, D to K). The only top-site start structure on Pt(111) surface also results in a top-site binding conformation (Fig. 5F). The changes upon chemisorption involve shortening of Pt─O distances from ~3 Å in the adsorbed state to Pt─O bond lengths near 2 Å in the chemically bound state (Fig. 5, D to K). The DFT calculations hereby shortcut a process that takes hours in experiment (6870). Depending on the initial structure and adsorbed O2 precursor orientations, we obtained a variety of chemisorbed structures. On the Pt(110) surface, starting on the threefold site, an associated bonded Pt─O─O─Pt structure was formed (Fig. 5D). This associated structure was by +1.27 eV less stable than a dissociated structure with two sets of Pt─O─Pt bonds that developed from O2 bound to a common epitaxial site (Fig. 5E). Epitaxial adsorption on the Pt(110) surface thus induced the oxygen molecule to approach Pt atoms across the long bridge (wide groove), elongating the O─O bond and promoting O─O dissociation (insets in Fig. 5E). The dissociated structure with four Pt─O bonds of ~2.0 Å in length was strongly favored for chemisorption, and similar dissociative adsorption was reported earlier on Ag(110) surfaces compared to Ag(111) and Ag(100) surfaces (74). On the Pt(111) surface, originating from a top site, a tilted single-bonded associated structure (Pt─O─O) was obtained (Fig. 5F), which is approximately +0.56 ± 0.03 eV higher in energy than three other parallel, vicinal associated bonded conformations originating from fcc, hcp, and bridge sites (Pt─O─O─Pt) (Fig. 5, G to I). Accordingly, two covalent Pt─O bonds are favorable over one covalent Pt─O bond on the Pt(111) surface. On the Pt(100) surface, an epitaxial initial conformation with a high 42° tilt angle transformed into a tilted, twofold geminal bonded associated structure (Pt > O─O) (Fig. 5J). The reaction process leading to this bound structure was much higher in energy, by +0.86 eV, compared to the transition from adsorption on a bridge site to a parallel, vicinal associated bonded structure (Pt─O─O─Pt) (Fig. 5K).

In summary, the oxygen molecules first need to overcome barriers from water adlayers to reach the surface atomic layer and tend to adsorb on nontop sites (steps 1 and 2 in Fig. 5L). Adsorption-desorption dynamics continues for many times without a reaction (step 3 in Fig. 5L). Proportional to the statistical amount of surface contact, chemisorption eventually occurs at the nearest top-layer Pt atoms (step 4 in Fig. 5L). The formation of covalent bonds reorients adsorbed O2 molecules toward top sites, shortens the Pt─O distance, and favors dissociated O2 molecules with up to four Pt─O bonds over two Pt─O bonds and one Pt─O bond if the local (hkl) geometry allows (Fig. 5, E and L). The stability of bonded species increases approximately in the order (Pt > O─O) < (Pt─O─O) < (Pt─O─O─Pt) < (2 Pt─O─Pt). An associative or dissociative mechanism of reduction may then follow as previously proposed (fig. S1, A and B) (33).

The results from DFT calculations (optPBE) also show that chemisorption energies of the covalently bound structures differ from −0.5 to −1.3 eV, equivalent to −11 to −30 kcal/mol. These variations are clearly too large to explain changes in specific ORR activity of only ~3 times on low-index Pt surfaces (Fig. 1, A and B) (5), which only equal changes in barrier heights of 0.025 eV (= 1.0 RT = 0.6 kcal/mol) and include contributions by adsorption. The DFT data therefore support that the rate is regulated by oxygen adsorption, i.e., contact time of oxygen with the Pt surface, which varies with surface structure and electrolyte composition in the experimentally observed range of ORR activities. We repeated the DFT calculations with the simpler PBE density functional for comparison, which changed the computed energy differences by up to 0.4 eV (table S6), while adsorption geometries and qualitative trends remained the same. Convergence testing as a function of energy cutoff and k-space grid points showed that an energy cutoff of at least 400 eV and k-point meshes above 4 × 6 × 1 yield steady energies (fig. S3, A and B).

The proposed mechanism agrees with available in situ spectroscopic data. Only small amounts of Pt─O and Pt─OH intermediates have been detected at 0.9 V versus RHE and remain difficult to identify (31, 73, 7577). In situ surface-enhanced Raman spectroscopy revealed the formation of intermediates on Pt(100), Pt(111), and Pt(110) surfaces during ORR at voltages higher than 1.1 V versus RHE and suggested low coverage of oxygen species at 0.9 V versus RHE (76). In situ x-ray spectroscopy and calculations on Pt nanoparticles showed that some adsorption was observed at 0.74 to 0.80 V versus RHE and O or OH groups bind to Pt top sites at low coverage (75), consistent with the geometries obtained after chemisorption in DFT simulations (Fig. 5, D to K). In situ x-ray absorption spectroscopy of ORR on Pt/Rh(111) catalysts and calculations indicated very low incidence of chemisorbed oxygen at 1 V versus RHE (77). While modification of Pt with Rh is expected to reduce binding energies and coverage, very low coverage is more consistent with the proposed mechanism than high coverage (not observed). Bisulfate electrolytes were also found to impede oxygen adsorption (75), consistent with the site-blocking effects seen in MD simulation (Fig. 3, C to E). In summary, in situ spectroscopy suggests low amounts of adsorbed oxygen species under ORR operating conditions (0.9 V versus RHE) and mainly records whether intermediates are present or not. It is not possible to quantify the coverage per square nanometer of surface area, the full chemical identity, or monitor adsorption-desorption cycles preceding the formation of chemisorbed oxygen species. The experimental and computational data suggest that bonded intermediates likely form after many adsorption-desorption cycles in a probabilistic manner regardless of the Pt surface structure. Sparse O2─Pt surface contact and weak attraction (Fig. 2, A to C) lead to the correlation of oxygen adsorption with ORR activity under typical operating conditions.

Correlation of O2 adsorption with ORR activity on irregular Pt nanostructures

Pt nanocatalysts used in practice are structurally complex, consist of multiple Pt(hkl) facets and defects, and often reconstruct during reactions (34, 35, 66, 78). The analysis of O2 adsorption on Pt nanostructures, which are covered by multiple (hkl) surfaces, therefore provides further insights under more realistic reaction conditions. We illustrate the correlation of O2 adsorption with ORR activity for several realistic, less-regular Pt nanostructures that have been well characterized with regard to structure and ORR activity in recent experiments (Fig. 6, A and B). (34, 35, 66) In comparison to flat (100), (111), and (110) surfaces (Fig. 6, C to E), these structures include R-NWs without grain boundaries (Fig. 6F), R-NWs with grain boundaries (Fig. 6G), J-NWs (Fig. 6H), and Pt nanoplates (Fig. 6I). Some of these systems have shown ultrahigh ORR activity in experiments (Fig. 1). The molecular models closely reflect the dimensions, surface corrugation, and grain boundary densities as observed in experiments (Fig. 6, F to I, and details in Methods) (34, 35, 66). We analyzed the interaction of O2 with the irregular Pt nanostructures in water by MD simulation, focusing on the O2─Pt coordination numbers, local Pt─Pt coordination numbers, the percentage of time adsorbed for every oxygen molecule (Fig. 6A), as well as the residence time during adsorption events (Fig. 6B and see fig. S3C with a logarithmic scale). These indicators are more broadly applicable because of the roughness and irregularity of the Pt surfaces and superficial water layers in comparison to free energy profiles that become strongly location dependent and difficult to sample across irregular surfaces. The coordination number of adsorbed O2 on Pt was defined as the number of Pt atoms within 3.5 Å in distance of either O atom in O2 (3.5 Å is the approximate cutoff of the first coordination shell as seen in Fig. 2, A to C). The Pt─Pt coordination number equals the number of Pt atoms within 3.5 Å in distance from a given Pt atom.

Fig. 6 Oxygen adsorption on Pt nanowires and Pt nanoplates.

(A) Coordination numbers and percent of time O2 molecules adsorb on Pt(hkl) surfaces and Pt nanostructures from unrestrained MD simulations over 50 ns (45 ns in equilibrium used for analysis). Higher coordination numbers and stickiness to the surface correlates with higher ORR activity observed in experiments. (B) Coordination numbers and residence time for every O2 adsorption event on Pt(hkl) surfaces and Pt nanostructures. Higher coordination numbers and residence time on the surface correlates with higher ORR activity observed in experiments. (C to E) Equilibrium models of the Pt(100), Pt(111), and Pt(110) surfaces, showing adsorbed conformations of O2, the coordination number of O2 with the Pt surface, and Pt─Pt coordination numbers (3.5-Å cutoff). (F to I) Models of R-NWs without and with grain boundaries, J-NW, and a Pt nanoplate with grain boundaries, as well as adsorbed conformations of O2 molecules. The models represent nanostructures in recent experiments (34, 35) and highlight several effective coordination patterns with high O2─Pt coordination numbers that lead to strong adsorption and extraordinary ORR activity.

Movies S1 and S2 visualize the adsorption-desorption dynamics in solution with 88 oxygen molecules over a period of 50 ns of simulation time. Although we use a high O2 supersaturation in solution (~0.36 M), the surface coverage of O2 on Pt in equilibrium is low on all surfaces (Fig. 6, C to I) related to weak adsorption free energies (Fig. 2, A to C). The low solubility of O2 in water of 0.00025 M (29) would further decrease the O2 coverage on a real catalyst surface by ~1000 times compared to the simulation. At the same time, some sites on the corrugated nanowires and nanoplates have remarkably higher O2 retention and adsorption than on any of the flat surfaces (Fig. 6, C to I). O2─Pt coordination numbers exceed seven for J-NWs and nanoplates, while typical coordination numbers are three to five on (110), (111), and (100) surfaces. Maximum O2─Pt coordination numbers of nine occur at buried Pt sites on nanoplates, or potentially on J-NWs depending on the model build (Fig. 6, A, B, H, and I). Some of these sites function as “oxygen traps,” and the percentage of contact time of several O2 molecules is 100%, which exceeds the typical range of 10 to 70% for regular Pt(hkl) surfaces (Fig. 6, A and B). The most attractive sites often also have higher Pt─Pt coordination numbers. On the coordination number–contact time percentage spectrum, regular nanowires are found in between flat surfaces and J-NWs and nanoplates (Fig. 6A), consistent with the ORR activity observed in experiments (Fig. 1B and table S7). R-NWs without grain boundaries show lower coordination numbers and lower contact time than R-NWs with grain boundaries, also consistent with experiments (66). Jagged nanowires and nanoplates show the highest coordination numbers, percent of time adsorbed, and residence time of O2, consistent with the highest ORR activities seen in experiments (34, 35). Here, experiments favor J-NWs over nanoplates, whereas simulations favor nanoplates over J-NWs. The difference could be related to some competition by ClO4 ions, which were not included in the simulation and may inhibit some sites on extended nanoplates, such as HSO4 ions, rather than on highly curved nanowires (Fig. 3, C to E).

Accordingly, we confirm the relationship between adsorption of molecular oxygen and ORR activity for irregular Pt nanostructures with some minor uncertainties. Higher O2─Pt coordination numbers, higher percentage of time adsorbed, and longer residence time correlate with increased ORR activity (Fig. 6, A and B, and table S7). Recent operando spectroscopic evidence supports that defective steps can act as templates to disturb interfacial water networks, similar to corrugated water adlayers (Fig. 2, A, D, and E), and enhance ORR activity (31). An interesting follow-on question is how to connect small potential shifts in cyclic voltammograms to the Pt surface–electrolyte structure, oxygen adsorption, residence times, and ORR pathways under realistic solution conditions.

Further evidence and application to other operating conditions and electrocatalysts

The small differences in oxygen adsorption energies on the order of 0 to 2 kcal/mol correlate with the relative range of known ORR activity (1 to 30) via the Boltzmann factor. Accordingly, pressure and temperature are also expected to play a role. The proposed mechanism is consistent with a linear increase in ORR rate for higher oxygen partial pressure in fuel cells (67), which increases the aqueous saturation concentration of O2 and, thereby, the contact time with the Pt surface. Higher temperature, such as 350 K in fuel cells, tends to slightly increase the overall ORR rate (67, 79). The exact temperature dependence involves kinetic and thermodynamic contributions. The kinetic contribution accelerates the rate of the slowest step of ORR at higher temperature at constant potential versus RHE, which has been observed in experiments (67, 79, 80) and is not the focus of our contribution. The thermodynamic contribution decreases the ORR rate at higher temperature associated with lower O2 solubility and modified adsorption equilibria, which can be quantified by MD simulation. On balance, the overall ORR rate may increase up to ~350 K and decrease at a higher temperature again (67, 79). Quantifying the two effects is still an open question. In the simulation, the free energy of O2 adsorption on Pt(110) slightly increased at 350 K relative to 298 K (−1.85 kcal/mol versus −1.54 kcal/mol), assuming same solubility. Simulation results are in the expected range of magnitude, and the relative activity of different Pt(hkl) surfaces and irregular surfaces at different temperatures can be further investigated.

For mixed-metal catalysts, such as Pt─Ni and other alloys, the ORR activity is likely influenced by two contributions, namely, by O2 adsorption and by changes in the electronic structure of the catalyst. The contribution of O2 adsorption to the rate is quantitatively predictable with MD simulation. The contribution of changes in the electronic structure to the rate can be estimated by a comparative analysis of the reaction path with and without alloying elements using DFT calculations at the local scale to estimate the overall change in activity. For example, Pt3Ni surfaces have shown an increase in ORR activity by a factor of 3 to 10 relative to pure Pt surfaces under otherwise identical conditions (5). Thereby, changes in the electronic structure of the surface accelerate the electron transfer. Follow-on simulations including solution conditions can yield interesting total quantitative correlations, for example, using combined MD and DFT methods such as quantum mechanics/molecular mechanics (QM/MM) to capture sensitive O2 adsorption and insight into the role of substituting elements and defects. Start configurations from realistic MD simulations could also much improve the quality of QM/DFT studies that often assume structures at 0 K without explicit water molecules and electrolytes, which greatly differs from the conditions in fuel cells and in rotating disk experiments. The MD methods will also be beneficial to examine the activity of metal-free catalysts for ORR and other reactions [nitrogen reduction reaction (NRR), carbon dioxide reduction reaction (CO2RR), oxygen evolution reaction (OER), hydrogen evolution reaction (HER)].


We have elucidated mechanisms of O2 adsorption on Pt surfaces and found a correlation of oxygen adsorption with ORR activity, which allows predictions of the relative activity of Pt catalysts under realistic operating conditions. O2 adsorption under solution conditions is the first step of ORR, preceding the formation of chemisorbed oxygen species, and has not been studied before. Accurate models for Pt and molecular oxygen in the IFF were derived and used to explain oxygen-platinum interactions in electrolyte solutions, quantifying inhibition by superficial water layers and coelectrolytes, site-specific oxygen binding, and molecular orientations. The methods can be used to estimate ORR activities for realistically sized Pt nanocatalysts up to the large nanometer scale using MD simulation at a computational expense more than 103 times lower than with DFT methods (see Methods).

A key role to modify oxygen contact with the surface is played by water networks on the Pt surfaces. Oxygen has more contact time when approaching certain Pt surface environments, such as corrugated, slightly curved, or locally concave surfaces that feature corrugated water adlayers with weaker hydrogen bonding. Stronger oxygen adsorption with longer residence time increases the probability to initiate chemisorption. Initial oxygen adsorption and time in contact with the Pt surface scale with the ORR rate, while initial chemisorption (or somewhere else along the reaction coordinate) is the rate-determining (slowest) step. According to the available data for pure Pt surfaces at 0.9 V versus RHE, this slowest step exhibits a constant rate regardless of the surface environment, and O2 adsorption regulates accessibility to the Pt surface. The differences in adsorption free energy are small (<2 kcal/mol) and correlate with the experimentally observed relative ORR activity (1 to 50 times), explaining quantitatively the range of ORR activity via equilibrium constants or the Arrhenius equation, which was not possible with prior models. Very low O2 saturation concentration in the electrolyte (0.00025 M = 8 parts per million) prevents high surface coverage under common operating conditions. On the basis of the available data, after the slow initiation of electron transfer, the remaining reaction sequence leading to water is likely faster, supported by the high concentration of hydronium ions of 0.1 M and the instability of transient bound oxygen species. As a limitation, there is still no experimental or computational technique that can monitor the activation barrier in solution and the reaction intermediates in real time in atomic resolution. The relationship between oxygen adsorption and ORR activity for flat Pt surfaces, irregular Pt nanostructures, and competition by electrolytes invites follow-on studies, for example, to identify Pt(hkl) surfaces and irregular surface geometries with highest ORR activity and uncover yet unknown correlations with ORR activity. We also recommend the inclusion of coadsorbing electrolytes in MD simulations.

Evidence for the correlation of oxygen adsorption with ORR activity on pure Pt catalysts includes multiple experimental and computational data on flat surfaces, studies on coadsorption of electrolytes, irregular structures, and the correlation with oxygen partial pressure. When extending rate predictions for other metals and alloys, we need to consider two major factors: adsorption of oxygen and the electronic structure. We expect the relative ORR rates to be regulated by the availability of oxygen on different surfaces, which is modulated by the solvent and competing electrolytes, and by the slowest step in ORR that involves initial electron transfer or a later step in the reduction reaction. The development of simulation tools for semiquantitative predictions of ORR rates at catalyst-solution interfaces is a major improvement over the lack of methods previously. Several aspects are still unexplored and require further experiments and simulations at multiple levels to clarify and further develop mechanistic insights and correlations with ORR activity. Understanding of catalyst-solution interfaces at this level of detail can benefit a multitude of reactions in electrocatalysis (OER, HER, NRR, and CO2RR) to find the role of catalyst composition, surface structure, and electrolytes for activity and stability. The methods can connect MD and QM approaches in new ways to study changes in local composition in alloy, single atom, and metal-free catalysts for quantitative predictions of the performance.


Computational models and methods summary

Models for all systems were prepared in all-atomic detail using the Materials Studio graphical user interface (81). MD simulations were carried out using the Large-scale Atomic/Molecular Massively Parallel Simulator program (82) at a temperature of 298.15 K and a pressure of 101.3 kPa consistent with common experimental conditions for tens of nanoseconds. We used a spherical cutoff at 12 Å for the summation of Lennard-Jones interactions and long-range electrostatics with high accuracy (10−5) using the Particle-Particle Particle-Mesh (PPPM) solver. Simulations were carried out using periodic boundary conditions in three dimensions except systems with applied potential. The periodic boundary condition in z dimension (direction of electric field) was then disabled to represent the electric potential.

New IFF models of oxygen with unprecedented performance were derived and used. Errors in critical interfacial properties are up to 10 times less than before (tables S1 and S2). The models reproduce the liquid density, vaporization energy, and free energy of hydration within 3% deviation from experimental measurements. We further developed a polarizable IFF model for platinum to represent electric polarization under applied fields in the simulation, including similar large improvements in accuracy over earlier models, embedded atom method (EAM), and DFT methods (tables S1, S3, and S4). The Pt-e model achieves an excellent match with lattice parameters, surface energies, mechanical properties (tables S3 and S4), as well as with the image potential (fig. S1C). The equivalent model for go4ld with available experimental data on monolayer hydration energies showed ±5% agreement without any additional fit parameters (39).

SMD simulations were used to determine adsorption free energy profiles of molecular oxygen with sufficiently slow, constant pulling velocity, ranging from 1 to 0.05 Å/ns using up to three independent trajectories. Adsorption conformations were analyzed from equilibrated MD trajectories up to 50 ns.

DFT calculations were carried out with the Vienna Ab initio Simulation Package (VASP) program (8386) using the optPBE and PBE functionals, an energy cutoff of 400 eV, 4 × 4 × 1 k-point grid for Pt(100) systems, and 4 × 6 × 1 k-point grid for Pt(110) and Pt(111) systems. The models included an adsorbed oxygen molecule, the first solvation shell of water, and a six-layer Pt slab. The initial coordinates were used from frames of the MD trajectory in equilibrium. Details are given in the following sections.

Calculation of adsorption free energy profiles

The free energy profiles of oxygen adsorption were obtained by SMD simulations (using the “fix smd” command in the LAMMPS program). Pt(100), Pt(110), and Pt(111) slabs constructed for SMD simulations were prepared from multiples of the metal unit cells, had a thickness of ~15 Å, and were infinite in the lateral dimensions using three-dimensional (3D) periodic boundary conditions. Pre-equilibrated water layers of 50 Å in height containing ~1000 explicit water molecules, represented by the flexible SPC model, were used above the Pt slabs. Single-target molecules (oxygen molecule, bisulfate ion, and phosphoric acid molecule) were initially positioned at least 12 Å away from the platinum surface atomic layer. A force was applied to the target molecule so that the target molecule was able to move toward to the surface at a constant velocity. When the SMD was applied to charged electrolyte molecules (bisulfate ions), no hydronium ions were involved to represent the acid conditions; instead, sodium ions were used to neutralize the simulation system. This change was made because hydronium ions interact with the steered molecule during the pulling process and result in poor sampling. The free energy was sampled using 1 Å/ns of velocity in bulk. As the target molecule approached the surface (~6 Å or less from the surface atomic layer), a slower pulling velocity of 0.05 to 0.5 Å/ns was selected to achieve sufficient sampling within the metal-water interface region. The most suitable velocity was selected as a tradeoff according between the different metal-water-adsorbate interfacial interactions and computation costs [0.05, 0.5, and 0.1 Å/ns for oxygen adsorption on Pt(100), Pt(110), and Pt(111) surfaces, respectively; and 0.1 Å/ns for electrolyte molecules on Pt(100), Pt(110), and Pt(111) surfaces]. The selected velocity was slow enough for every system to reach equilibrium at every time step to simulate a reversible process. The results do not change for smaller pulling velocity. For water adlayers with a discontinuity to bulk water such as Pt(100) and Pt(111), a slow pulling velocity as given above is necessary to obtain an accurate free energy profile. To avoid movement of the platinum slab, the bottom atomic layer of the slab was constrained in the pulling direction of the target molecules while mobile in the other directions. The simulation systems were first pre-equilibrated for 1 ns in the isobaric-isothermal (NPT) ensemble and then in the isochoric-isothermal (NVT) ensemble to carry out the SMD runs. The trajectory frames and the corresponding free energies were recorded every 1 ps in bulk and every 0.5 to 0.05 ps within 6 Å from the surface atomic layer. The Visual Molecular Dynamics (VMD) program was used to determine the distance between the geometric center of the target molecule and surface atomic layer at each frame. The reported free energy profiles were generated by converting the free energies from all frames into a function of the distances between the target molecule and surface atomic layer and then smoothed using the percentile filter method. The geometric centers of oxygen molecules, sulfur atoms, and phosphorus atoms were used as reference points to calculate the distance from the top atomic layer.

Quantitative analysis of adsorption probability and molecular orientation at each site

The site preference of oxygen adsorption on Pt(hkl) surfaces was analyzed using equilibrium MD trajectories. The models contained a Pt slab, water molecules, and a single O2 molecule, as used in the SMD simulation. Unrestrained MD simulations were carried out in the NPT ensemble at a temperature of 298.15 K and a pressure of 1 atm for 50 ns. Snapshots (coordinates) were collected every 10 ps. Frames with O2 adsorbed on Pt in the equilibrium trajectories were distinguished from frames depicting nonadsorbed O2 using a cutoff distance between the surface atomic layer and the oxygen molecule of 3.0 Å. The lateral coordinates (x and y coordinates in Fig. 4, A to C) of adsorbed oxygen molecules were then wrapped into a unit cell to visually superimpose all adsorbed positions [Pt(100): a = b = 3.9236 Å; Pt(110): a = 3.9236 Å, b = 2.7744 Å; Pt(111): a = 2.7744 Å, b = 4.8054 Å]. The geometric centers of the oxygen molecules were used as reference points for the position of the molecules. Each adsorbed position was assigned to a surface site according to the distance between the adsorbed position and the center position of the nearest labeled surface sites. After assigning all adsorbed orientations to a certain site, the adsorption probability of each site was then determined by the ratio of number of adsorbed frames pertaining to this site to the total number of adsorbed frames. The adsorption angle of each site was calculated by averaging all adsorption angles of adsorbed oxygen molecules assigned to the site. The trajectories were visualized by the VMD program and analyzed using customized new MATLAB code and Python code with the MDAnalysis package (87, 88).

Analysis of ORR mechanism using DFT

All-atom models for DFT calculations were chosen from equilibrated MD trajectories using representative oxygen adsorption sites and orientations on each Pt surface. Because of the high cost of DFT calculation, only six layers of Pt atoms and the first hydration shell around adsorbed oxygen molecules were included in the simulation boxes. A five-angstrom cutoff from adsorbed oxygen molecules was used to determine the first hydration shell and results in 13 to 16 water molecules included in simulation boxes. Box dimension for Pt(110), Pt(111), and Pt(100) systems are 15.6944 Å by 11.0960 Å by 30 Å, 14.4162 Å by 11.0976 Å by 30 Å, and 11.7708 Å by 11.7708 Å by 30 Å, respectively. Periodic boundary conditions were used with ~10-Å vacuum above the hydration shell in z direction. Geometry optimization was first conducted until the formation of Pt─O bonds was observed. No proton transfer occurred during this process. To compare the energy differences between different adsorption conformations on the same Pt(hkl) surface, hydration shells were removed to maintain model systems with the same number of atoms. Then, geometry optimization was performed again for all systems until convergence, and systems with the highest energy were used as a reference state for other adsorbed conformations on the same Pt(hkl) surface. All DFT calculations were carried out using the VASP program with the projector augmented wave pseudo-potentials (89, 90). Both generalized gradient approximation (GGA) optPBE and PBE functionals were applied, and the results suggested negligible effect on the calculated energies according to the reference states (table S6). The convergence tests were performed on the energy cutoff and k-point grid using the final structure in Fig. 4E with optPBE functional, and an energy cutoff of 400 eV with 4 × 6 × 1 k-point grid showed sufficient accuracy (fig. S3, A and B). The 4 × 4 × 1 k-point grid was applied for Pt(100) systems due to the equal cell dimensions of 11.7708 Å.

Analysis of O2 adsorption on irregular Pt nanostructures

Models of R-NW without grain boundaries, R-NW with grain boundaries, J-NW, and Pt nanoplates were prepared using dimensions, grain boundary density, and surface corrugation as described and tested for ORR in previous experiments (34, 35, 66). Models of R-NW without grain boundaries were built starting with an elongated bulk Pt supercell of 157.5 Å in length with (111) surfaces at the ends, which was radially cut to yield a cylinder of 1.8 nm in diameter. For the model of Pt nanowires with grain boundaries, the R-NW was cut into ~20-Å-long segments, which were rotated around the long axis at random angles to create grain boundaries. The surface atoms were annealed and cooled down to relax the surface structure. Last, the Pt nanowire with grain boundaries was subjected to 2 ns of MD simulation at 298.15 K in vacuum for equilibration before use. Building models of J-NW started with Pt─Ni alloy as a precursor, as reported in synthesis (34). A 30,600-atom Pt supercell with a dimension of 52.714 Å by 45.651 Å by 163.102 Å was used, and 26,400 Pt atoms were randomly replaced with Ni atoms to match the overall ratio of Ni to Pt to experiment (34). The Pt─Ni structure was equilibrated for 5 ns of MD simulation in vacuum at 443.15 K using the NPT ensemble with pressure applied in z direction, followed by the removal of all Ni atoms. The resulting structure formed a J-NW with a diameter of ~3 nm upon 5 ns of further MD simulation. Subsequently, further core atoms were randomly removed, and the model equilibrated in the NVT ensemble until a J-NW of 2.2 nm in diameter was obtained, consistent with transmission electron microscopy (TEM) images (34). The model of the nanoplates was assembled from 36 individual cuboctahedral Pt nanoparticles with 3.5 nm in diameter, following the mechanism reported in experiments (66). The constituting nanoparticles were placed at close distance next to each other at about the same height in a simulation box with 188 Å by 200 Å cross-sectional area in vacuum, bringing the nanoparticles in a position to fuse into a nanoplate. Then, MD simulation in the NPT ensemble (with pressure only in the two lateral dimensions) was carried out at 1 atm and 298.15 K to transform the dispersed nearby nanoparticles into a nanoplate with an equilibrium dimension of 146.2 Å by 157.4 Å. The model of the nanoplate was cut down in size into a cross-sectional area of approximately 94 Å by 104 Å, followed by 2 ns of further equilibration in the 2D NPT ensemble in vacuum at 298.15 K. The final dimension of the equilibrated nanoplate was 93.46 Å by 103.21 Å. The grain boundary density (total length of gain boundaries divided per surface area) was 0.60 nm−1 in the nanoplate model and 0.58 nm−1 in TEM images (66).

For the analysis of O2 adsorption, the Pt nanostructures, 13,476 water molecules, and 88 oxygen molecules were included in the model system. The solution was supersaturated with oxygen molecules (~0.36 M) due to limitations in system size for good statistics. The oxygen molecules had no pairwise interactions in solution or on the surface that would affect results besides having more frequent surface contact for analysis. Larger slabs were rebuilt for the Pt low-index surfaces [(110), (111), and (100)] to enable unbiased comparisons with the irregular Pt nanostructures (see models in Fig. 5, C to I). All models were initially equilibrated for 200 ps in the NPT ensemble and then in NVT ensemble for 55 ns. All systems reached steady-state adsorption-desorption equilibria after 10 ns, and snapshots from the equilibrium trajectories were collected every 10 ps to obtain 4500 trajectory frames for analysis.

Python code was developed using the MDAnalysis package (87, 88) to determine the percentage of time adsorbed and the residence time. The number of Pt neighbor atoms of each oxygen molecule was calculated on every frame with a radial distance cutoff at 3.5 Å (O2─Pt coordination number). A 4500 × 88 matrix was generated, with rows suggesting the number of neighboring Pt atoms for every O2 at a certain time step and columns, suggesting the number of neighboring Pt atoms for certain O2 during the 45-ns simulation. An O2─Pt coordination number of zero indicated oxygen desorption, and nonzero coordination numbers indicated adsorption. From this matrix, the percentage of time adsorbed could be determined for every O2 molecule as the ratio of the number of frames indicating adsorption to the total number of frames (4500). The residence time for every adsorption event was analyzed by counting the number of frames with continuous nonzero neighbor Pt atoms (= continuous nonzero O2─Pt coordination number). Adsorption events of <100 ps (less than 10 frames) were eliminated to avoid statistical noise caused by the fixed 3.5-Å distance cutoff. The Pt─Pt coordination numbers were computed using a cutoff at 3.5 Å.

Computational advantage

DFT simulations take 106 to 108 times longer than MD simulations for running a system with the same size for the same amount of time (91). MD simulations can therefore access a larger length and time scale including irregular surfaces and accurate electrolyte interactions and remain a factor of ~103 times faster. If DFT simulations are added locally, then the total computational cost remains effectively the same as standard DFT.

For example, simplified DFT calculations may use model systems of 200 atoms compared to 10,000 atoms in MD. The simulation time is then 10 ps with DFT while it is 100 ns in MD, overall leading to a speed up of 103 times or more. Major differences are the computed properties: in case of MD, the dynamic adsorption-desorption equilibria and ion competition, and in case of DFT, qualitative insights into electron transfer and follow-on binding conformations. DFT is not necessary to analyze adsorption of molecular oxygen.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. Liu, S. Hoff, and K. Kanhaiya for discussions. We also thank J. Fan for coding support. Funding: This work was supported by the NSF (DMREF 1623947, CBET 1530790, and OAC 1931587). The allocation of computational resources at the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357, and at the Summit supercomputer supported by the NSF (award numbers ACI-1532235 and ACI-1532236) is acknowledged. Author contributions: S.W. performed the simulations, collected the data, analyzed the data, and wrote the manuscript. E.Z. and Y.H. analyzed the data and wrote the manuscript. H.H. designed the study, analyzed the data, and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article