Origin of high oxygen reduction reaction activity of Pt12 and strategy to obtain better catalyst using sub-nanosized Pt-alloy clusters

In the present study, methods to enhance the oxygen reduction reaction (ORR) activity of sub-nanosized Pt clusters were investigated in a theoretical manner. Using ab initio molecular dynamics and Monte Carlo simulations based on density functional theory, we have succeeded in determining the origin of the superior ORR activity of Pt12 compared to that of Pt13. That is, it was clarified that the electronic structure of Pt12 fluctuates to a greater extent compared to that of Pt13, which leads to stronger resistance against catalyst poisoning by O/OH. Based on this conclusion, a set of sub-nanosized Pt-alloy clusters was also explored to find catalysts with better ORR activities and lower financial costs. It was suggested that Ga4Pt8, Ge4Pt8, and Sn4Pt8 would be good candidates for ORR catalysts.

sub-nanosized clusters 7 . Therefore, their electronic structures and the relationship between their structures and catalytic activities are unclear.
As experimental methods such as X-ray diffraction cannot be applied to investigate sub-nanosized Pt clusters, we have to investigate them using theoretical methods. In order to obtain a plausible theoretical geometry for each metallic cluster, however, it is necessary to use global optimization, which is a method that determines the most stable geometry in a multi-dimensional potential energy surface for a given system [10][11][12][13] . For an N-atom system, the dimension of the structural phase space is 3 N. The number of local minima is known to increase exponentially with respect to N 14,15 . Generally speaking, it is impossible to find all stable potential minima in 3N-dimensional space using a simple geometry optimization technique 16 . Furthermore, local optimizations from nonsensical initial geometries are often trapped in local minima with high energy; however, this is not important for discussions on catalysis. Moreover, the experimentally observable physicochemical properties of each cluster cannot be entirely understood through only a global minimum geometry owing to the existence of minima with similar potential energies 17,18 . For example, Rodríguez-Kessler et al. have mentioned that O/OH adsorption energies for Pt 12 /Pt 13 depend on the cluster geometries, i.e., fluctuation 19 . It should be noted that the adsorption energies have often been utilized as simple indicators of ORR activity [20][21][22] . Previous reports have clearly shown that geometrical flexibility and the accompanying electronic fluctuation should be considered when discussing ORR activities, at least in the case of sub-nanosized Pt clusters.
The aim of this study was to identify novel sub-nanosized clusters that have higher ORR catalytic efficiencies with lower financial costs from a theoretical point of view. Among the sub-nanosized Pt clusters, we focused on Pt 12 , as the atomicity can be controlled precisely using the method developed by Yamamoto, Imaoka et al. [6][7][8][9] . Since it is not completely understood why Pt 12 showed better ORR activity than Pt 13 did, we focused on Pt 13 as well. Thus, we investigated the difference between the electronic structures of Pt 12 and Pt 13 to reveal the cause of the difference in their ORR activities, taking the flexibility of their geometries into consideration. First, we performed ab initio molecular dynamics and simulated annealing 13 (AIMD-SA) calculations to find a global minimum and relatively low-energy geometries for Pt 12 and Pt 13 . The fluctuation behaviors of electronic structures within a set of stable geometries were compared. Second, to understand the difference between the electronic structures in terms of interactions with O/OH molecules that affect the rate-limiting steps of the ORR, we calculated adsorption energies for the fluctuating geometries using ab initio Monte Carlo (AIMC) simulations. Finally, based on a discussion for pure Pt clusters, we searched for several sub-nanosized alloy clusters containing Pt that can be novel catalysts with high activities and low financial costs.

Results and Discussion
An example of AIMD-SA, in which the initial geometry belongs to the highly symmetric I h point group, is shown in Fig. 1. It can be observed that the atomic configuration was well mixed to give an apparently different geometry from the initial one. The AIMD procedure succeeded in preventing the system from being trapped in local minima. It should be mentioned that the final geometry in Fig. 1 is ca. 3 eV more stable than the initial one in the I h point group, which has been considered a topological magic number. This result indicates well the importance of finding global minima in transition metal clusters whose electronic structures are governed by complex interactions between d-orbitals. Figure 2 shows final geometries of Pt 12 and Pt 13 obtained using AIMD-SA simulations started from nine different initial ones. The Cartesian coordinates of each geometry are written in the Supporting Information (SI). The energetically most probable geometries for Pt 12 and Pt 13 are (a) and (j), respectively. These geometries are consistent with previous global optimizations based on density functional theory (DFT) by Zhang 23 , Da Silva 24 , Wei 25 , and Zhai 26 et al., which focused on obtaining a static global minimum geometry in each cluster size. Although it is difficult to obtain a global minimum with a single simulated annealing (SA) run, we can increase the precision when obtaining a global minimum geometry by performing multiple independent AIMD-SA runs. In the present study, six AIMD-SA runs with different initial conditions gave the same geometry (j) for Pt 13 , which means that it should be the global minimum 16 . It was found that both Pt 12 and Pt 13 have isomers whose energies are within  13 . The vertical and horizontal axes are the energy relative to that for the initial geometry and simulation time (steps), respectively. The simulation temperature in each step is given above. The geometries depicted below represent snapshots in the AIMD-SA.
Scientific RepoRts | 7:45381 | DOI: 10.1038/srep45381 0.5 eV of those of the most stable ones. This indicates that the geometries of Pt 12 and Pt 13 fluctuate under ambient conditions. Our results clearly show the importance of the geometric fluctuations of the clusters. This agrees well with the previous report by Rodríguez-Kessler et al. 19 .
To investigate the differences in the electronic properties of Pt 12 and Pt 13 , first, we compared their atomic natural charges [27][28][29] and densities of states (DOSs) in the "static" global minimum geometries ((a) and (j)). Contrary to our expectations, as shown in Figures S1 and S2, we could find no apparent differences in their charge distributions. This is because we forgot that these clusters should have enough flexibility to allow their geometries to fluctuate under the experimental conditions. In order to discuss the differences between the electronic structures of Pt 12 and Pt 13 more precisely, it is mandatory to consider their atomic fluctuations. Next, we compared the thermal fluctuations of their electronic structures using wave functions at all atomic configurations obtained in our AIMD-SA simulations. Using this scheme, we can acquire thermally important geometries and accompanying electronic structures. Figure 3 depicts the relationship between atomic charge variance and effective coordination number (ECN) 30 at energies less than 0.8 eV. The vertical and horizontal axes in Fig. 3 indicate the atomic fluctuation and electronic polarization in a cluster, respectively. Upon comparing the data for Pt 12 and Pt 13 , we could observe no apparent difference in graph width in the vertical direction. This means that geometric softness, which can be viewed as atomic fluctuation, is similar for Pt 12 and Pt 13 . On the other hand, in the horizontal direction, we could observe a clearly wider distribution for Pt 12 . Thus, we could conclude that Pt 12 has a more flexible electronic structure than Pt 13 does even though their degrees of atomic fluctuation are similar.  Next, we shall move our focus onto O/OH adsorption interactions on the sub-nanosized Pt clusters. To date, interaction energies of O/OH and Pt have often been used to explain their catalytic activities [20][21][22] . The relationship between the interaction energies and corresponding ORR activities has been known to give a volcano plot. This is because O and OH are key intermediate species in rate-limiting steps of the ORR. The plot tells us that the adsorption energy should not be too strong or too weak in the ORR catalytic process. Basically, the idea behind volcano plots can be well applied to describe Pt-based ORR catalytic materials. Since the basic idea of the volcano plot is based on "static" interaction between adsorbent and adsorbed species, however, strictly speaking, it should be noted that discussions involving the plot should be applied only to well-fixed solid-state materials. As explained above, the sub-nanosized Pt clusters are geometrically flexible at room temperature. Thus, there should be an enormous number of possible rate-limiting complexes with different interaction energies for the sub-nanosized Pt catalysts. In order to explain the difference between the ORR activities of Pt 12 and Pt 13 , again, we have to consider the thermal fluctuation effect of the clusters on the adsorption interaction. In this study, we applied a set of AIMC simulations to deal with the energetic fluctuation in the adsorption interaction. The details of the AIMC simulations are described in the Methods section and SI. Using the AIMC scheme, 10 5 geometries for each complex were sampled and the effects of thermal fluctuation were taken into account.
In Fig. 4, it is shown that the O/OH adsorption energies to Pt 12 /Pt 13 have wide energy spreads extending over 1.5 eV. In the Figure, data for different adsorption sites (i.e., on-top, bridge, and three-hollow sites) were plotted with different colors. Recently, using the global minimum geometry of Pt 13 , Chaves et al. showed that the adsorption energies of OH in on-top and bridge sites are − 3.23 and − 3.60 eV, respectively 31 , which are similar to our data in Fig. 4. We could observe that there should be a large fluctuation for the interaction energy by thermal atomic motion even in the same kind of adsorption sites. Using a few optimized geometries, Rodríguez-Kessler et al. showed that the adsorption energies of O/OH and Pt 12 /Pt 13 depend on the cluster geometries 19 . Since the atomic charge of an O atom represents the amount of charge transfer between it and Pt-based ORR catalysts, it has often been considered that the atomic charge on an O atom should be a good indicator for the adsorption strength. However, in Fig. 4, where the flexibilities of the clusters are fully taken into account, we could not find such a strong correlation between the vertical and horizontal axes. Moreover, we could find no apparent difference between the degrees of atomic fluctuation (horizontal direction in Fig. 3) of Pt 12 and Pt 13 . These results mean that a factor other than atomic fluctuation, i.e., electronic flexibility, of the sub-nanosized cluster governs the difference between the catalytic activities of Pt 12 and Pt 13 . Table 1 shows the average and standard deviations of the O/OH adsorption energies. The standard deviations are 0.17, 0.27, 0.05, and 0.03 eV for Pt 12 -O, Pt 12 -OH, Pt 13 -O, and Pt 13 -OH, respectively. These results correspond well with the fact that Pt 12 has a more flexible electronic structure than Pt 13 even though they show similar magnitudes of atomic fluctuation (see Fig. 3). To put it briefly, Pt 12 undergoes greater electronic fluctuation in O/OH adsorption than Pt 13 does, which results in the prevention of catalyst poisoning. Our results explain well the experimental data reported by Imaoka et al., which indicated that Pt 12 shows higher ORR activity than Pt 13 does 6 .
One of the roles of current theoretical chemists is to predict novel functional materials. In this study, we also attempted to search for sub-nanosized alloy materials as candidates for ORR catalysts. As discussed above, it  has been found that electronic fluctuation of sub-nanosized clusters should play an important role in preventing catalyst poisoning. Although it has been clarified that better sub-nanosized catalysts for the ORR should have large electronic fluctuations to prevent poisoning effects, random atomic fluctuation is an inconvenience in the further design of catalysts based on sub-nanosized materials. Thus, a strategy for reducing uncontrollable atomic fluctuations must be used. The geometrical flexibility of the Pt clusters originates from complex hybridizations of d-orbitals. In this paper, we propose a novel approach for developing better ORR catalysts based on the sub-nanosized Pt clusters. This approach is the sub-nanosized alloying of the Pt clusters with p-block metals. Since the valence electronic structures of p-block metals are much simpler than those of the d-block ones, the geometries of d/p-alloyed clusters are expected to have moderately fixed atomic configurations and large electronic fluctuations. In this study, we explored the possibility of using M 4 Pt 8 (M = Al, Ga, Ge, and Sn) as ORR catalysts using the AIMD-SA and AIMC scheme mentioned above. Since Al, Ga, Ge, and Sn have lower electronegativities than that of Pt (1.61, 1.81, 2.06, 1.96, and 2.28, respectively), alloying them with Pt clusters is expected to enhance electronic polarization and fluctuation in the electronic structure of the alloy clusters. At the same time, as these metals are much cheaper than Pt is, these M 4 Pt 8 (M = Al, Ga, Ge, and Sn) alloy clusters are good candidates for ORR catalysts from an elemental strategy point of view. The number of heteroatoms here was chosen to suit the experimental synthesis technique developed by Yamamoto, Imaoka et al. [6][7][8][9] . The most stable geometries obtained for M 4 Pt 8 (M = Al, Ga, Ge, and Sn) are depicted in Fig. 5. For Al 4 Pt 8 , we could obtain only one geometry even though we performed several AIMD-SA runs with different initial conditions. This indicates that this geometry is most likely to be a global minimum 16 . For Ga 4 Pt 8 , Ge 4 Pt 8 , and Sn 4 Pt 8 , we found a few local minima within the energy range of 0.3 eV (data not shown). Figure 6 shows the relationships between geometric and electronic fluctuations for M 4 Pt 8 (M = Al, Ga, Ge, and Sn). It should be noted that the vertical axis of Fig. 6 is ten times larger than that of Fig. 3. The variances of atomic charges are 0.11, 0.03, 0.03, and 0.11 for Al 4 Pt 8 , Ga 4 Pt 8 , Ge 4 Pt 8 , and Sn 4 Pt 8 , respectively, which are much larger than those for Pt 12 (0.009) and Pt 13  (0.0075). This means that polarizations of the sub-nanosized alloy clusters are larger than those of pure Pt clusters. The averaged atomic charges of the p-block metal and Pt, respectively, are + 0.46 and − 0.23 for Al 4 Pt 8 , + 0.21 and − 0.11 for Ga 4 Pt 8 , + 0.22 and − 0.11 for Ge 4 Pt 8 , and + 0.47 and − 0.23 for Sn 4 Pt 8 . The negative values of the Pt atoms in the alloy clusters can be explained reasonably using Pauling's electronegativity. Alloying with p-block metals also enhanced charge polarization. That is, in the alloy system, electronic fluctuation is enhanced even though geometrical fluctuation is reduced.
Since we were able to control the atomic fluctuation, as expected, by introducing p-block metals, the next step was to examine the catalytic activities of the sub-nanosized alloy clusters. Figure S3 shows the adsorption energies between M 4 Pt 8 (M = Al, Ga, Ge, and Sn) and adsorbed species (O/OH). Compared to those in Fig. 4, the distributions in Figure S3 are wider. The energy-weighted average and standard deviation of each adsorption are listed in Table 1. From these data, first, we can find that all of the alloy clusters gave similar adsorption features to those of Pt 12 /Pt 13 . This means that the alloy clusters should be candidates for ORR catalysts. However, it seems that Al 4 Pt 8 binds with O/OH too strongly to catalyze the ORR. In the panel for Al 4 Pt 8 -OH in Figure S3, we can find adsorption sites whose adsorptions are too strong, with energies less than − 4 eV (even though the energetic fluctuation is as large as that of Pt 12 ). Poisoning by OH is expected to occur in Al 4 Pt 8 . For the other alloy clusters, no such poisoned site was found. Thus, it can be concluded that Ga 4 Pt 8 , Ge 4 Pt 8 , and Sn 4 Pt 8 should be good candidates for ORR catalysts. Ge 4 Pt 8 is especially promising because it has similar adsorption energy to Pt 12 with larger fluctuations.

Conclusions
The origin of the superior ORR activity of Pt 12 to that of Pt 13 can be explained by differences in electronic fluctuations. As Pt 12 and Pt 13 have numerous stable isomers, it was clarified that their physicochemical properties cannot be discussed without consideration of their atomic fluctuations. Taking their atomic fluctuations into account, we found that Pt 12 has a more flexible electronic structure than Pt 13 does, even though their degrees of atomic fluctuation are similar. Greater electronic fluctuation results in greater energetic fluctuation of O/OH adsorption for Pt 12 , i.e., stronger resistance against catalyst poisoning compared to Pt 13 . It was concluded that controlling the atomic/ electronic fluctuations is key in realizing higher ORR catalytic activities for sub-nanosized clusters containing Pt. Based on conclusions for pure Pt clusters, we attempted to predict whether sub-nanosized M 4 Pt 8 (M = Al, Ga, Ge, and Sn) alloy clusters, in which p-block metals were alloyed with Pt to enhance electronic flexibility and provide well-defined atomic configurations, would be good candidates for ORR catalysts. It was suggested that Ga 4 Pt 8 , Ge 4 Pt 8 , and Sn 4 Pt 8 would show better ORR catalytic activities than Pt 12 and Pt 13 do and with lower costs.

Methods
In the present study, all calculations were performed based on DFT using the TURBOMOLE 7.0 quantum chemical program package 32 . We employed the Perdew-Burke-Ernzerhof (PBE) functional 33 , which is a pure generalized gradient approximation (GGA) exchange-correlation functional, to describe the electronic structures of the target clusters. As Pt is a heavy element belonging to the fifth row in the periodic table, it is mandatory to take relativistic effects into account for sufficiently accurate quantum chemical discussions. Thus, we applied the def-SV(P) basis set with relativistic pseudopotentials 34,35 . Xiao et al. have reported that the relative energy of each isomer is not affected by spin-orbit coupling in Pt clusters 36 . Therefore, we did not include the spin-orbit coupling effect in our calculations. To reduce computational costs for two-electron integrals and to accelerate all calculations, the resolution of identity (RI) approximation was applied 37 . As our targets are sub-nanosized transition metal clusters, we had to consider their band-like (metallic) electronic structures with thermal excitation 38 . In the present study, the metallic character was described using the pseudo-Fermi smearing technique at 300 K, which is consistent with experimental conditions 39 .
To find the global and low-energy minima of Pt 12 , Pt 13 , and M 4 Pt 8 (M = Al, Ga, Ge, and Sn), simulated annealing (SA) using ab initio molecular dynamics simulations (AIMD) was performed, in which electronic potential energy and gradient were evaluated using the quantum chemical method mentioned above. Under the experimental conditions, the sub-nanosized clusters were supported on glassy carbon [6][7][8][9] . According to the report by Lim et al., the geometries of Pt 13 in the gas phase and on defective graphene supports are similar 40 . They also reported that ORR pathways are the same in both cases. Thus, in the present study, we used cluster models in gas-phase conditions. In the AIMD-SA simulations, the Leapfrog algorithm was applied. The details of the AIMD-SA simulations were as follows. First, a set of initial geometries was prepared. For Pt 12 and Pt 13 , nine geometries were chosen from the previous reports by Häkkinen, in which geometries of Au n (n = 4-14) were investigated using photoelectron spectroscopy and DFT calculations 41 . For the alloy clusters, M 4 Pt 8 (M = Al, Ga, Ge, and Sn), three sets of initial geometries were defined by randomly replacing four Pt atoms in stable structures of Pt 12 with M. Then, after generating the initial geometries, a set of the AIMD-SA simulations was performed with a time step of 10 fs for Pt 12 , Pt 13 , and M 4 Pt 8 (M = Ga, Ge, and Sn) and 5 fs for Al 4 Pt 8 . These time steps were chosen based on one-tenth of the vibration period of Pt 2 or PtM (M = Al, Ga, Ge, and Sn) dimers ( Table 2). The AIMD-SA simulations consisted of four stages to avoid being trapped in a local minimum. In the first stage of our simulations, a set of AIMD at 2000 K was performed for 5.2 ps. The object of this stage was to scramble cluster configurations. In the second and third stages, the temperature was lowered gradually to 1000 and 500 K, respectively. In these three stages, an NVT ensemble with Nosé-Hoover thermostats was applied 42 . Finally, in the fourth stage, a set of SA runs was performed to predict a global minimum structure for each cluster. The final geometries are depicted in Figs 2 and 5.
The ECNs in Figs 3 and 6, and S4 were calculated using Eqs. (1) and (2). The ECN is a fractional number defined as the sum of weights, representing the contribution to the coordination number of the target atom 30 . The contribution of atom j to atom i is calculated based on the distance between them (d ij in Eqs. (1) and (2)  The original definition of the ECN, Eqs (1) and (2), is not suitable for alloy systems, which contain multiple types of atoms. In the present study, we replaced d ij in Eqs (1) and (2)  The average ECNs and variances of atomic charges were calculated using Eqs. (4) and (5), respectively. A large variance of atomic charge means the cluster is strongly polarized in the geometry. To understand the differences in the adsorption interactions between ORR intermediates O/OH and each sub-nanosized cluster, we developed an AIMC simulation code by ourselves using Python with the NumPy module 44 . The detailed procedure of the AIMC simulation is as follows. First, we selected a set of sub-nanosized cluster geometries from pre-performed AIMD-SA simulations by random sampling with energy weight (Boltzmann factor). Second, an adsorbed molecule of O/OH was placed randomly near the surface of the cluster selected in the first step. In this step, the shortest distances between the adsorbed molecule and cluster were defined to be those in the 1.3-3.0 Å range. The geometries that did not satisfy the requirement were not adopted from the AIMC simulations. Then, using the accepted geometries, geometry optimizations were performed using the quasi-Newton-Raphson method with convergence criteria of 10 −6 a.u. For all sub-nanosized clusters considered in this study, we obtained 10 5 adsorption geometries with O/OH. Finally, the interaction energies between the sub-nanosized clusters and adsorbates were evaluated. In the geometry optimizations, all atoms were fully relaxed to allow geometry changes through adsorption interactions. In the final complex geometries, the interaction energies were calculated using Eq. (6).

ot tot
Further details of the AIMC procedure are described in the Supporting Information.