Bond-type embedded CGCNN (BE-CGCNN)
Previous studies using the ML approach to predict surface adsorption energies mostly focused on single adsorbates27, 28, 29, 30, 31, 32, 33, 34. However, they performed poorly in predicting adsorption energies for cases involving multiple adsorbates and were thus not successful in reflecting surface coverage effects. In this regard, an ML model that can cover various adsorbates and coverages must be developed. For high surface coverage cases, adsorbates on NP surfaces are close enough that several types of interactions can occur, including intermolecular and intramolecular interactions. To reflect such complexity, we propose the bond-type embedded CGCNN (BE-CGCNN) model in Fig. 1, in which edge vectors are uniquely designed.
BE-CGCNN overall follows the schemes of the original CGCNN approach, in particular, those for graph construction, node vectorization, and selection of convolution functions. For each NP and slab structure, atoms and bonds are encoded into node vectors and edge vectors to construct a graph of the corresponding structure. As shown in Fig. 1b, for edge vectors (representing bonds), we classify them into four types: covalent bond within an adsorbate (e.g., O-H), metallic bond within an NP (e.g., Pt-Pt), chemisorption between an NP and an adsorbate (e.g., Pt-O), and, lastly, nonbonded interaction between different adsorbates (e.g., H…O), in which the edge vectors are encoded in a one-hot manner with four categorized vectors. The last term (nonbonded interaction) is only valid when the atomic distance is larger than 1.25 Å. In the previous CGCNN and its modified versions, long-distance interactions were not treated fairly. However, in many cases, there could be strong nonbonded interactions, such as hydrogen bonds, between OH species, and these were treated as important in our graph constructions. Note that the edge vectors are designed to be distance-insensitive.
The processes of node vector construction and optimization are basically identical to those of the original CGCNN. To select appropriate features for the node vector, we calculated the mean absolute error of the adsorption energy with increasing number of features, as shown in Fig. 1c. Candidates for the features include elemental properties available in the periodic table of elements as follows: group number, period number, atomic number, radius, electronegativity, ionization energy, electron affinity, volume, atomic weight, melting temperature, boiling temperature, density, Zeff, polarizability, resistivity, capacity, number of valence electrons, and number of d-electrons. The value ranges and categories are provided in Supplementary Table 1. For cost efficiency, the best feature combination set was fixed with increasing number of features. The best result was obtained for the feature set comprising group number, atomic radius, and electron affinity. Because the adsorption of adsorbates is closely related to the electronic interaction between adsorbates and NPs, the model with selected features related to the number of valence electrons (group number and atomic radius) and the energy of valence electrons (electron affinity) likely shows the best performance.
Adsorption energy dataset for various NP surface coverages.
In this work, we aim to build surface Pourbaix diagrams of Pt NPs, which are well-known materials in various catalytic applications. The adsorption of O and OH species on Pt NPs was calculated for training set construction. Adsorption of O and OH on Pt13 and Pt55 NPs was included. Two NP structures, cuboctahedron (Coh) and icosahedron (Ih), which are known as the stable morphologies of NPs with 13 and 55 atoms, were considered. NPs with a truncated octahedron (Toh) structure with 38 atoms were also included. Slabs with (100), (110), (111), and (211) exposed surfaces were also included. The structures of these Pt NPs are available in Supplementary Fig. 1. For each NP and slab structure, adsorption configurations of O and OH up to 1 monolayer (ML) were modeled. O and OH were adsorbed on either bridge or top sites, which were the most stable adsorption sites on the Pt(111) slab. The total adsorption energy (\(\varDelta \text{E}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]\)) and the adsorption energy per adsorbate (\(\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]\)) were computed by the following equations:
$$\varDelta \text{E}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]=\text{E}\left[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}\right]-\text{E}\left[\text{N}\text{P}\right]-n\text{E}\left[\text{A}\right]$$
$$\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]=\frac{\varDelta \text{E}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]}{n}$$
where E[NP-(Aads)n], E[NP], and E[A] are the total energies of the structures of the NP including n adsorbates, the NP only, and adsorbate A only, respectively. A refers to an adsorbate, such as O and OH.
Unfortunately, neither the total adsorption energy (\(\varDelta \text{E}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]\)) nor the adsorption energy per adsorbate (\(\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]\)) is suitable for accurate ML training for the following reasons. For the former (the total adsorption energy), the data range was estimated to be too large (> 120 eV; Supplementary Figs. 2c and 2d) because of the cases involving several tens of adsorbates, and thus, the absolute errors of the ML models are also very large. On the other hand, for the latter case (the adsorption energy per adsorbate), the data range is much smaller (< 3 eV; Supplementary Figs. 2a and 2b), and the ML errors seemingly look small. However, in the end, surface Pourbaix diagrams are fed with total adsorption energy inputs, and thus, the predicted values for many adsorbate cases (relatively large surface coverage) would be very erroneous and misleading. To overcome this limitation, we introduced a different metric, namely, the adsorption energy difference (\(\varDelta \varDelta \text{E}\left[\text{n}\text{p}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}\right]\)), which served as a much more suitable form for accurate ML training and prediction and could be computed as follows:
$$\varDelta \varDelta \text{E}\left[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}\right]=\varDelta \text{E}\left[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}\right]-n\overline{\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]}$$
$$\overline{\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]}= \frac{\sum _{}^{\text{M}}\varDelta {\text{E}}_{ads}\left[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}\right]}{\text{M}}$$
where ∆∆E[NP-(Aads)n] is the adsorption energy difference, \(\overline{\varDelta {\text{E}}_{ads}[\text{N}\text{P}-{\left({\text{A}}_{\text{a}\text{d}\text{s}}\right)}_{n}]}\) is an averaged value of adsorption energies per adsorbate at each n, and M is the number of each n adsorbate case in our dataset. The value range of this metric (∆∆E[NP-(Aads)n]) is much smaller (~ 25 eV; Supplementary Figs. 2e and 2f) than that of the total adsorption energy (> 120 eV); thus, the ML errors were also estimated to be much smaller.
BE-CGCNN training results with bond-type embedding.
In Fig. 2, we show the BE-CGCNN model training results obtained using the dataset comprising 736 adsorption energy difference data points (both O and OH species on Pt NPs), in which 80% of the dataset was used for training and the remaining 20% was used for the test. For both adsorbates, the BE-CGCNN with bond-type embedding greatly outperforms the original CGCNN, exhibiting a mean absolute error (MAE) of 0.33 eV for the O adsorbate and an MAE of 0.07 eV for OH. These values are much smaller than those of 0.86 eV and 0.49 eV for the O and OH cases from the original CGCNN model (‘without bond-type embedding’ results in Fig. 2), proving the effectiveness of bond-type embedding in our ML models. Note that these MAE values are remarkably small given that the data range of the adsorption energy difference (∆∆E[NP-(Aads)n]) is as large as ~ 25 eV. The error values of the original CGCNN are not sufficiently small to warrant reliable construction of Pourbaix diagrams, and thus, implementation of BE-CGCNN with much enhanced accuracy is required.
Different bond-type embeddings lead to different error levels, as shown in Table 1. For both O and OH species, the original CGCNN without bond-type differentiation exhibits the worst accuracy. Upon the introduction of three bond types (covalent bond, metallic bond, and chemisorption) into the bond vector encoding, the error is significantly reduced to 0.44 eV and 0.30 eV for each O and OH dataset. Finally, by adding a nonbonded interaction term to the bond types, the ML model reaches even lower-level errors of 0.33 eV and 0.07 eV, respectively. These successful error reductions indicate that the nonbonded interaction is a critical term and should be treated as important, particularly for high surface coverage cases where van der Waals interactions are present.
Bond embedding | Test MAE (O DB) [eV] | Test MAE (OH DB) [eV] |
Without bond type embedding | 0.86 | 0.49 |
Bond type (C, M, CH) | 0.44 | 0.30 |
Bond type (C, M, CH, NB) | 0.33 | 0.07 |
C, covalent bond; M, Metallic bond; CH, Chemisorption; NB, Nonbonded interaction |
Table 1. MAE values with the variations of bond type embedding.
Constructing Pt surface Pourbaix diagrams using BE-CGCNN.
Now, using the well-trained BE-CGCNN model, we are ready to build surface Pourbaix diagrams without DFT inputs. Our model systems are Pt NPs whose surfaces are under the competition between O and OH adsorption. The first validation process is to construct reliable Pourbaix diagrams of icosahedron-shaped (termed Ih hereafter) Pt55 NPs using the BE-CGCNN model trained on the data of smaller systems, including slabs and smaller NPs (Pt13 and Pt38). In Fig. 3, the surface Pourbaix diagram of a Pt55 NP (Ih) constructed solely by the BE-CGCNN model is compared to the ground truth diagram obtained by DFT computations. The DFT and ML results are observed to be very similar, exhibiting differences in the y-intercepts of boundary lines of even less than 0.1 eV on average. In both diagrams, as the pH and U (applied potential) increase, transitions from bare Pt NPs to OH-covered Pt NPs and finally to O-covered Pt NPs are observed. In addition, Pt dissolution regions are found at low pH (approximately less than 3) and high U (approximately greater than 0.7 V). This test reveals the feasibility of reliable construction of Pourbaix diagrams of larger-size Pt NPs solely based on BE-CGCNN.
Next, we expand the study of ML-based Pourbaix diagrams to Pt NPs of various sizes (Pt55 and Pt147 NPs) or shapes (Coh and Ih), as shown in Fig. 3. Here, the BE-CGCNN model was trained on the dataset of slabs and relatively small NPs (Pt13, Pt38, and Pt55). In Fig. 3a-3d, the surface Pourbaix diagrams of the Pt55 (Ih) NP and Pt55 (Coh) NP are shown. The diameter of Pt55 NPs is approximately 0.5 nm. Qualitatively, as the pH and given potential U increase, the phase transition from bare Pt to the OH- and O-adsorbed surface occurs. The DFT-based and ML-based predictions appear quantitatively very similar for both shapes. For example, in the Ih shape cases, the phase boundaries between the bare Pt NP and Pt-(OH)0.29ML appear at almost identical y-intercepts (0.52 V for DFT versus 0.51 V for ML). Other boundary lines within comparative diagrams of Ih cases appear at very similar positions, with the differences in terms of the y-intercept positions being much less than 0.1 V on average. The diagrams for Coh NP cases are more complicated due to the appearance of many more phases. Similar to Ih cases, the ML-based predictions overall also worked great for Coh cases, except for the quantitatively erroneous description of the boundary line between Pt-(O)0.67ML and Pt-(O)1ML (y-intercept difference of approximately 0.24 V).
We further explore larger-size NP systems, namely, Pt147 (Ih) and Pt147 (Coh), as shown in Figs. 3e and 3f. The diameter of the Pt147 NPs is approximately 0.8 nm. For both Ih- and Coh-shaped NPs, Pt-(OH) phases are substantially destabilized over other phases of bare Pt and Pt-(O) 35 in comparison to diagrams of smaller Pt NPs (e.g., Pt55) or Pt slab systems2 available in Supplementary Fig. 4. This observation is consistent with the tendency of relatively strong OH adsorption on Pt surfaces for smaller Pt NPs compared to bulkier NPs36.
A prominent difference between Pt147 (Ih) and Pt147 (Coh) is the relative amount of the fully O-covered phase (Pt-(O)1ML). The Pt-(O)1ML phase stands out for the Ih case, whereas this phase shrinks for the Coh structures due to the appearance of partial O coverage phases, such as Pt-(O)0.53ML. Reduction of the Pt-(O)1ML phase was also similarly found for smaller Pt55 NPs (Coh), where partial O coverage phases (Pt-(O)0.42ML, Pt-(O)0.56ML, Pt-(O)0.67ML) preferentially appear over the full O coverage phase. This difference comes from the shape effect of the NPs. Unlike the Ih structures comprising only (111) surface planes, Coh NP structures have mixed surfaces of (111) and (100) planes. The adsorption of O species on the Pt(100) surface is well known to be stronger than that on the (111) surface due to the presence of more dangling bonds37. As a result, for the Coh NP structures, partial O phases (Pt-(O)<1ML) are likely to be stable when oxygen is dominantly adsorbed on (100) surfaces compared to the full O phase (Pt-(O)1ML) where oxygen is adsorbed on both (100) and (111). This phenomenon would be unlikely for the Ih NP structures without (100) surfaces.
Pourbaix diagrams for real-scale Pt NPs.
The NPs explored thus far are NPs involving up to 147 atoms. This size corresponds to only 0.8 nm in diameter, which is unfortunately far from the real scale. Typically, the diameters of experimentally synthesized NPs are over 3–4 nm, which involves thousands of atoms. Using BE-CGCNN, we demonstrate the construction of Pourbaix diagrams of real-scale Pt NPs, including Pt561 (Coh), Pt3871 (Coh), and Pt6535 (Coh) NPs (approximately 2.8, 3.9, and 4.8 nm in diameter), in Figs. 4a-4c. Since these diagrams are for several nanometer-size Pt NPs, we can now compare the results with the available experimental reports. For example, there are studies reporting the measured onset potentials for surface oxide (or fully O-covered phase) generation on Pt NPs, including 0.96 V on 1.2 nm-size NPs from Merte et al. 38 and 0.9–1.15 V on 4.0 nm-size NPs from Mom et al. 39. These values are marked on the Pourbaix diagrams in Figs. 4a and 4b and are found near the boundary lines of Pt-(O)1ML in each diagram. Such quantitative agreement well supports that the constructed Pourbaix diagrams are highly reliable and that the size dependences are greatly reflected.
As the NP size increases, the OH-covered phases are observed to be reduced compared to the O-covered phase. This is consistent with the experimental report that the OH surface coverage decreases with increasing NP size35. The Pourbaix diagrams asymptotically converge above the size of 3,871 atoms. The converged case can be compared to the Pourbaix diagram of the Pt(111) slab system, available in Supplementary Fig. 5: the compositions of the O-covered phases are quite different, whereas those of the bare Pt phases and OH-covered phases are very similar. This difference arises because the Coh-shaped NPs have (100) terraces, unlike the Pt(111) slabs.
The increasing O- to OH-covered phase ratio with increasing NP size can be understood in terms of the relative adsorption strength of O and OH species. Figure 4d compares the adsorption energies per adsorbate of the fully O-covered (Pt-(O)1ML) and OH-covered (Pt-(OH)1ML) phases, and their difference becomes larger with increasing NP size, which is consistent with the shrunken OH phases for larger NPs. The NP size dependence of the adsorption energies can be adequately explained by the relative adsorption strength on vertex sites, edge sites, and terrace sites. Taking Pt55 (Coh) NPs as an example, the adsorption of OH on a vertex site is 0.68 eV stronger than that on a terrace site, whereas the adsorption of O on a vertex site is only 0.27 eV stronger than that on a terrace site. The smaller the Pt NPs become, the larger the ratio of vertex and edge sites to terrace atom sites becomes, and thus, the adsorption energy difference between O and OH species (Eads, O 1ML – Eads, OH 1ML in Fig. 4d) is gradually reduced. These trends are clearly presented in Fig. 4d.
In addition to the O- to OH-covered phase ratio, the Pt dissolution phase is also an interesting spot to focus on. As the NP size increases, the Pt dissolution phases shrink. The Pt dissolution area is closely related to the stability of NPs operating in electrochemical catalysis. For example, Pt NPs are one of the most widely used catalysts for the oxygen reduction reaction (ORR) at the cathode of proton exchange membrane fuel cells (PEMFCs). However, the operating conditions of the PEMFC cathode are approximately an applied potential of 0.8 V and a pH of 140. This point is inside the Pt dissolution domain for relatively small Pt NPs (Pt55 and Pt147), as shown in Fig. 3, whereas it resides outside of the Pt dissolution region for larger Pt NPs (Pt561, Pt3871, and Pt6525), as shown in Figs. 4a-4c. Very interestingly, an experimental report41 confirmed that Pt NPs with diameters smaller than 2.0 nm were easily dissolved, and thus, the specific electrochemical surface area was greatly decreased over potential cycling compared to larger NPs.
We importantly note that the experimental operating conditions (0.8 V, pH = 1) are near the fully O-covered phase, which is a surface Pt oxide. Although this position is outside of the Pt dissolution region, the surface Pt oxide layer induces nonequilibrium transient dissolution of Pt42, 43, 44. In these situations, Pt dissolution may occur even for large-size Pt NPs, and the ORR performance of the catalyst would be degraded under long-term working conditions. This indicates that to fully understand the electrochemical stability of NP-based catalysts, taking into account kinetic factors, not only the surface Pourbaix diagram in which only equilibrium states are considered, is also important.
The construction of Pourbaix diagrams of several nanometer-size NPs was enabled by the extremely fast prediction speed of BE-CGCNN compared to DFT. Thus, discussing the computing time taken for the task of predicting the adsorption energy of each structure by the BE-CGCNN and DFT methods should be worthwhile. For Pt147 NPs, the total computing time (both training and prediction) for BE-CGCNN was estimated to be approximately 150 seconds based on a personal computer implemented with an NVIDIA GPU of GTX 2070. In contrast, the computing time for the same task was approximately 90 hours (2,160 times longer than in the BE-CGCNN case) based on a high-performance computing node implemented with a 2.3 GHz 20-core CPU, as shown in Fig. 5. Because DFT theory follows the computational scaling of O(N3), where N is the number of electrons45, 46, the computing time differences between BE-CGCNN and DFT will be extremely large for NPs involving thousands of atoms. Following the extrapolation lines based on O(N3) for DFT, for the example of the NP of 6,535 Pt atoms, the computing time of DFT will be 2,200 days, which is approximately 1.9 × 108 times longer than that in the BE-CGCNN case.