Implementation of Wang-Landau Algorithm for Probing Thermodynamic Stable Configuration of Multi-Element Materials and Application to Multinary Alloy Nanoparticles

Abstract

With an increase in an element, the configurational entropy effect stabilizes the multi-element materials. However, the expansion of configuration space and heterogeneous surfaces such as nanoparticles preclude the analytical evaluation of configurational entropy. Then, we implemented the Wang-Landau algorithm, which is one of the Monte-Carlo algorithms, for evaluating the configurational entropy and probing the thermodynamic stable configuration of multi-element materials. The regression equation obtained by density functional theory calculation and multiple regression analysis is used in the energy estimation in the sampling. This method was applied to binary alloys in the bulk and ternary alloy nanoparticles and the obtained features of the stable configuration were discussed.

1 INTRODUCTION

High-entropy alloy nanoparticles (NPs) have been attracted due to their excellent catalytic activity and durability [1,2,3]. With an increase in the constituent elements, the configurational entropy effect enhances the alloy NP stability. Solid-solution (SS) binary alloy NPs consisting of two elements, that are immiscible in the bulk state, have been synthesized [4, 5]. Although PdRu alloy NPs, which are one of the examples, showed excellent catalytic activity [6, 7], they are unstable under the practical use environment. Adding third elements to binary alloy NP is one of the methods for stabilization [8]. The high entropy effect is effective for not only alloy NP but also other multi-element materials. However, with increases in element and size, the configurational space becomes more expanded as shown in Figure 1. Evaluating the configurational entropy and probing the stable configuration becomes more difficult.

Figure 1.

 Number of configurations of alloy NP models consisting of 201, 405, 711, and 1289 atoms.

A simple calculation of entropy is based on the Bragg-Williams approximation [9]. The estimation of the entropy is used in the stability assessment of alloy NPs through the genetic algorithm [10]. However, the configuration in this approximation includes all types such as the SS, ordered, and segregated types. The real configurational entropy is lower than the ideal configurational entropy [11]. To estimate the configurational entropy of a specific configuration, the cluster variation method was introduced [12,13,14,15]. The combinations of bonds, tetrahedrons, and octahedrons divide into specific configurations. This method is specified when the structural features of the consistent elements such as coordination number (CN) are homogeneous. The CN of atoms in the surface layer is different from that of the bulk state. The inhomogeneous feature should be introduced. When the configurational entropy of SS and skin layer configurations were compared, the components of the surface and interior were considered [16, 17]. However, the analytical estimation with the cluster variation and inhomogeneous surface is extraordinary for the multi-element materials.

In the numerical estimation of configurational entropy, the multicanonical Monte Carlo algorithm [18,19,20,21,22] is considered. The transition probability from state A with the energy level E_A to state B with energy level E_B is expressed as   

$p\left(E_A⟶E_B\right)= min\left(1.0,\frac{g\left(E_A\right)}{g\left(E_B\right)}\right)$(1)

.

The configurational density of states (DOS) g(E), meaning the number of possible configurations for an energy level E, determines the transition probability, which is different from the Metropolis sampling determined by the temperature and energy difference (E_A − E_B). However, g(E) must be known beforehand. Then, Wang and Landau suggested a new algorithm to solve this problem [23, 24]. The configurational DOS is asymptotically estimated. This method was applied to the spin system, and the critical temperature was revealed. Besides, two-dimensional DOS g(E, M) with an energy level E and a magnetization M serving as the variables provides the expected value of the magnetization [25, 26]. The stable configuration of the spin system was revealed. The obtained configurational entropy results in the partition function. Various thermodynamic variables can be derived using the partition function. Therefore, the Wang-Landau algorithm is useful to evaluate a large number of configurations for multi-element materials.

Herein, we implemented the Wang-Landau algorithm to evaluate the configurational entropy and stable configuration of the multi-element materials. Wang-Landau algorithm was combined with the density functional theory (DFT) calculation and machine learning. The approach was successful in assessing stable configuration for ABC ternary alloy NPs (A = Pd, B = Ru, and C = Cu, Rh, Ir, or Au) [27]. However, the previous study has focused on specific systems based on PdRu alloy NP. It must be clarified that our method is effective in the systems without one or both of Pd and Ru. Our findings indicate that the surface energy of the third element C, relative to those of Pd and Ru, has a significant influence on the stable configuration and the temperature dependence of alloy NPs. On the other hand, designing functional and durable materials requires a comprehensive understanding of their thermodynamic stability. To this end, the implemented method was applied to the corresponding binary and ternary alloys in the bulk and NP models, providing more detailed insights into the thermodynamic stability.

2 IMPLEMENTATION

Energy estimation is a crucial step in the Wang-Landau sampling process. For a spin system, the Ising model was considered, and the energy was estimated by the formula [23,24,25,26]. In the context of alloy models, the energy is calculated by the interaction between the constituent elements instead of the up- and down-spin. We calculated the excess energies of various configurations using DFT. The calculated excess energies were subjected to multiple regression analysis (MRA). The regression equation based on the DFT calculations was used to estimate the energy during sampling.

2.1 DFT calculation

Spin-polarized calculations were conducted using the Vienna ab initio simulation package code based on DFT [28, 29]. Perdew–Burke–Ernzerhof exchange-correlation functionals [30] with the generalized gradient approximation were used. A projector-augmented wave method was used to model the interactions between valence and core electrons [31, 32]. The cutoff energy of the plane wave basis set was set to 400 eV. Monkhorst−Pack k-point grids were 9 × 9 × 9 in the bulk model, whereas the k-point grid was 1 × 1 × 1 in the NP model. The convergence criteria for the self-consistent field and geometry optimization were set to 10⁻⁵ and 10⁻⁴ eV, respectively. The volume, shape, and atom in the bulk models were relaxed, while the atom in the NP model was relaxed. In the NP models, the vacuum space of the unit cell was assumed to be greater than 12Å to avoid interaction with the neighbor NPs.

We prepared several types of alloy models, including those with a SS configuration. In the SS alloys, the constituent elements are randomly mixed. The Warren-Cowley parameter was used to evaluate the homogeneity of the SS in the binary alloy models [33]. This parameter is defined as   

$\alpha =1-\frac{P_{AB}}{c_B}$(2)

, where P_AB and c_B represent the probability of finding the B atom surrounding the A atom and the overall composition of the B atom, respectively. When Warren-Cowley parameter α is close to zero, the two elements in binary alloy models are mixed. In the ternary alloy models, the homogeneity of the SS state was evaluated using pairwise short-range order [34] as follows:   

${\alpha }_{AB}=\frac{P_{AB}-c_B}{{\delta }_{AB}-c_B}$(3)

, where δ_AB = 1 (A = B) and δ_AB = 0 (A ≠ B). The SS configuration forms when the value of the short-range order is close to zero.

To assess the stabilities of the binary alloy models, we calculated the excess energy expressed as   

${\epsilon }_{excess}=\frac{1}{N}\left({\epsilon }_{AB}-\frac{N_A}{N}{\epsilon }_A-\frac{N-N_A}{N}{\epsilon }_B\right)$(4)

, where ε_A and ε_B represent the total energies of A_N and B_N models, respectively. For the ABC ternary alloy NP (N = N_A + N_B + N_C), we calculated the excess energy expressed as   

${\epsilon }_{excess}=\frac{1}{N}\left({\epsilon }_{ABC}-\frac{N_A}{N}{\epsilon }_A-\frac{N_B}{N}{\epsilon }_B-\frac{N_C}{N}{\epsilon }_C\right)$(5)

, where ε_A, ε_B, and ε_C represent the total energies of A_N, B_N, and C_N NP models, respectively.

2.2 Machine Learning

MRA in supervised learning (SL) was performed. For the bulk model, leave-one-out validation was used. We calculated the predictive squared correlation coefficients and mean absolute error (MAE) by the following equations:   

$R^2=1-\frac{{\sum }_i^n{\left(y_{i }-Y_i\right)}^2}{{\sum }_i^n{\left(y_{i }-\stackrel{-}{y}\right)}^2}$(6)

and   

$MAE=\frac{1}{n}{\sum }_i^n\left(y_{i }-Y_i\right)$(7)

, where Y and y denote the excess energy predicted by SL and calculated by DFT, respectively. The bar symbol represents the average value. In the NP model, the hold-out method was used to validate the analysis, and this was repeated five times. In each validation, 1/4 and 3/4 of the prepared configurations were randomly selected as the test and training sets.

2.3 Wang-Landau sampling

Figure 2 shows the algorithm suggested by Wang and Landau. The state with energy E_A moves to that with energy E_B by the operation such as displacing and exchanging atoms. The transition probability in the multicanonical Monte-Carlo algorithm was evaluated by the configurational DOS g(E). When a trial move is accepted, the configurational DOS corresponding to the new configuration is updated by multiplying the existing value by a modification factor f. If a trial move is rejected, the configurational DOS corresponding to the old configuration is updated. In this study, the only exchanges of atoms were considered to probe the configurations. Figure 3 (a) shows the example of the increment of configurational DOS, where an arrow indicates a transition from state A to state B. The block in the Figure represents the configurational DOS, which refers to the number of possible configurations for an energy level. In the ln (g(E)) representation, multiplying by f( = e¹) is equivalent to the increment of one block. When the configurational DOS of state A is larger than that of state B, a trial move is accepted. When the configurational DOS of state A is smaller than that of state B’, a trial move is determined by the transition probability. With an increase in the difference of configurational DOS, the transition probability becomes lower. For example, the transition probabilities for transition with a difference of ln (g(E)) of 1 and 2, indicated by blue and magenta arrows in Figure 3 (a), are 0.368 and 0.135, respectively. Together with the update of configurational DOS, a histogram is incremented. The configurational DOS was asymptotically obtained by repeating the calculation until convergence was reached. The flatness criterion of the histogram was set to 0.8. When the histogram was flat, the histogram entries were reset to 0 and the modification factor was reduced from f to f^1/2. The resulting increment in configurational DOS in Figure 3 (b) is half of that in Figure 3 (a). When the configurational DOS of state A is larger than that of state B, a trial move is unchanged. On the other hand, the transition probability from state A to state B’ with larger configurational DOS varies according to the increment of configurational DOS. For instance, when the difference of ln (g(E)) is 1/2 for the transition (Figure 3 (b), magenta arrow), the transition probability is 0.607. The modification factor used in the simulation was determined with reference to previous studies [23,24,25,26]. The initial modification factor f₀ was assumed to be Napier’s constant. When the modification factor was smaller than exp (10⁻⁸), the sampling ended. The obtained configurational DOS g(E) results in the partition function and Helmholtz free energy expressed as   

$Z= {\sum }_{i}g\left(E_i\right)e^{-\beta E_i}$(8)

,   

$F=-k_{B}Tln\left(Z\right)$(9)

, where β and k_B represent the thermodynamic beta ( = 1/k_BT) and the Boltzmann constant. The two-dimensional DOS g(E, σ) was calculated to determine the stable configuration, where the second variable was the structural parameter σ instead of magnetization M in a spin system. The structural parameter σ in the bulk model corresponds to the bond fraction for the alloys. In the case of the alloy NP, the composition of the constituent elements in the vertex (CN = 6), ridge (CN = 7), and facet (CN = 8&9) are considered as well as the bond fraction in the whole NP and surface layer. Satisfying the flatness criterion of the two-dimensional histogram is difficult [24, 25]. The less stringent was used; The number of entries that are larger than 2000 remains unchanged for N × 10⁶ trials. The obtained configurational DOS g(E, σ) results in the expected value of structural parameter expressed as   

$⟨\sigma ⟩=\frac{1}{Z}{\displaystyle \sum }_{i,j}{\sigma }_{j}g\left(E_i,{\sigma }_j\right)e^{-\beta E_i}$(10)

.

Figure 2.

 Algorithm of Wang-Landau sampling

Figure 3.

 Increment in configurational DOS for (a)f = e¹ and (b)f = e^1/2. The red and green arrows indicate a transition to state B with lower or equal confirmational DOS, while the blue and magenta arrows represent a transition to state B’ with larger configurational DOS. The colored block represents an increment in configurational DOS resulting from the acceptance or rejection of a transition. Transition probability for the operation is expressed by p.

Several expected values were obtained for the structural parameters. In this study, one of the possible configurations with structural parameters close to the expected values is shown.

3 RESULTS

3.1 Bulk model

We took Pd−M (M = Rh, Ag, Ir, or Au) and Ag−Rh alloys as examples of binary alloys, where all constituent elements show the fcc structures. The excess energies of Pd−M and Ag−Rh alloys were calculated to obtain the regression equation of energy estimation in Wang-Landau sampling. The bulk models consist of 32 ( = N) atoms as shown in Figure 4 (a), where their composition x/N of A_xB_N−x alloy is 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, and 0.875. We prepared five configurations (α = α_min, α_min/2, 0, α_max/2, and α_max) for each composition. Figure 4 (b) shows the compositional dependence of excess energies in Pd−M and Ag−Rh alloys. The excess energies of Pd−Rh, Pd−Ir, and Ag−Rh alloys are positive, and increase in the order: Pd−Rh < Pd−Ir < Ag−Rh. On the other hand, the excess energies of Pd−Ag and Pd−Au alloys are negative. Ag and Rh atoms are immiscible [35], while the SS Pd_0.5Rh_0.5 and Pd_0.5Ir_0.5 can be formed above 1200 and 1500 K [36, 37]. Pd−Au alloys show the ordered state [38]. The trend of the excess energy is consistent with that of the transition temperature in the phase diagram.

Figure 4.

 (a) Bulk models, (b) excess energies, (c) predictive squared correlation coefficients as a function of the number of samples, and (d) comparison between excess energies determined by MRA and by DFT calculation for A_xB_32−x alloys.

The excess energies of alloy models were subjected to MRA, where the A−A and B−B bond fractions were used as the descriptor. A total of 35 configurations were considered in this study. Figure 4 (c) shows the predictive squared correlation coefficients as a function of the number of samples. The predictive squared correlation coefficient of the Pd−Ag alloy was lower than those of the other alloys. The low predictive squared correlation coefficient is caused by the small difference between the maximum and minimum values of the response variable [39]. On the other hand, the MAE of the Pd−Ag alloy was much the same as those of the other alloys. With an increase in the samples, the predictive squared correlation coefficients improved. When the number of samples exceeded 30, the predictive squared correlation coefficients were consistently above 0.70 and converged for each combination. Therefore, a total of 35 configurations were found to be sufficient for expressing the excess energies in the bulk model through MRA. Figure 4 (d) shows comparisons between excess energy predicted by SL and calculated by DFT in Pd−M and Ag−Rh alloys. The excess energies predicted by SL matched those calculated by DFT. The MAE of Pd−Rh, Pd−Ag, Pd−Ir, Pd−Au, and Ag−Rh are 2.96, 6.42, 8.76, 5.03, and 7.18 meV/atom. The regression equation is enough to estimate the excess energy in the sampling. Table 1 shows the regression coefficients of the excess energies of Pd−M and Ag−Rh alloys. The order of the intercept β₀ is consistent with that of the excess energy in A_0.5B_0.5. The sign of the intercept β₀ is opposite to those of the regression coefficients β₁ and β₂. With a decrease in the A−B bond fraction, the absolute values of excess energies become smaller. As shown in Figure 4 (b), the compositional dependence of the excess energy at A-rich composition is different from that at B-rich composition, which results in the difference between the regression coefficients β₁ and β₂. Thus, the obtained regression equation is reflected by the characteristics of the calculated excess energies.

Table 1.  Intercept (β₀) and regression coefficients of A−A and B−B bond fractions (β₁ and β₂) for the excess energy of the A−B bulk model

	β0	β1	β2
Pd−Rh	0.1451	−0.1510	−0.1549
Pd−Ag	−0.080	0.1001	0.0574
Pd−Ir	0.2228	−0.2299	−0.2430
Pd−Au	−0.1359	0.1563	0.1159
Ag−Rh	0.4729	−0.5103	−0.4696

In Wang-Landau sampling, the bulk model was expanded. We considered the bulk model consisting of 256 atoms. When the same calculation conditions were used with the exception of k-point grids (1 × 1 × 1), the error between the excess energy calculated by DFT and estimated by the regression equation is close to MAE for A_xB_32−x bulk models. Thus, the regression equations of the excess energies were used in the energy estimation for the A_xB_256−x bulk model in the Wang-Landau algorithm. The sampling procedure is specifically shown in Pd_xRh_256−x alloy as an example. In the Pd−Rh alloy, the SS configuration forms at high temperatures, while Pd and Rh are segregated at low temperatures. According to the phase diagram, the transition temperature between SS and segregated regions exhibits a maximum (1188 K) around 0.5 of Pd composition and varies in an inverse parabolic shape with respect to Pd composition. Figure 5 (a) shows the configurational DOS of Pd_xRh_256−x, which results in the partition function and excess free energy (equations (8) and (9)). Excess free energy depends on the composition and temperature as shown in Figure 5 (b). The excess free energies of Pd_xRh_256−x alloys are unchanged at low temperatures. When the slope of excess energy for temperature changes, the stable configuration of the alloy model changes. The excess free energy of Pd_0.5Rh_0.5 is negative above 1000 K, indicating that Pd−Rh alloy forms the SS configuration. The temperature showing the negative excess free energy becomes lower with an increase/decrease in Pd composition. Although the temperature of the miscibility gap is a little high, the stability of the Pd_xRh_256−x alloys is close to the phase diagram. Next, we reveal the stable configuration of binary alloys. The expected value of the A−B bond fraction is an easy comparison for understanding the stable binary alloy. Figure 6 shows the composition and temperature dependence of the expected A−B bond fraction in Pd−M and Ag−Rh alloys together with the stable configurations of A₆₄B₁₉₂, A₁₂₈B₁₂₈, and A₁₉₂B₆₄. The stable configuration changes by combination, composition, and temperature. The variation depends on the sign of excess energy. Pd−Rh, Pd−Ir, and Ag−Rh alloys, in which the excess energies are positive, show fewer A−B bonds. With an increase in temperature, the A−B bond fractions of Pd−Rh and Pd−Ir alloys become higher. The energy gain from entropy effect overcomes the energy loss of the A−B bond. The random structure of the binary alloy model shows 0.5 of the A−B bond fractions. However, the Pd−Rh bond fractions of Pd_0.5Rh_0.5 at 1000 K are about 0.33, which suggested that the stable configuration of Pd−Rh alloys in the high-temperature region is not completely random structure. The stable configuration of Ag−Rh alloy did not show a large variation. The energy gain from the entropy effect did not overcome the energy loss of Ag−Rh bonds. On the other hand, the alloys with negative excess energy (Pd−Ag and Pd−Au) show the ordered configuration at 0 K. The A−B bond fractions in the ordered configuration are maximum. With an increase in temperature, the A−B bond fractions become lower. The energy gain from the entropy effect overcomes the energy gain from the A−B bond. The Pd−Au bond fraction of Pd_0.5Au_0.5 alloy at 1000 K is about 0.55, which is close to the bond fraction of the random structure. The Wang-Landau algorithm with DFT calculation and MRA revealed the detailed temperature variation of the stable configuration of the binary alloy in the bulk state.

Figure 5.

 (a) Configurational DOS as a function of excess energy and (b) temperature and Pd compositional dependence of excess free energy.

Figure 6.

 Temperature and compositional dependence of A−B bond fraction and stable configurations of A_xB_256−x alloys.

3.2 Nanoparticle model

The enthalpy of the alloy in the bulk state is determined by the bond fraction. On the other hand, the surface affects the stability of the alloy NP. The elements with low surface energy preferred the low-CN atomic positions [40, 41]. In our previous study, the comparison between PdRuIr and PdRuAu alloy NPs showed that the stable configuration and the temperature variation depend on the third element [27]. In PdRuAu alloy NPs, Au atoms are located at low-CN atomic positions. With an increase in CN (vertex → ridge → facet), Pd composition increases. With an increase in temperature, the Pd composition in the vertex increase while those in the facet decrease. The characteristic indicates that Pd and Au atoms are mixed at the vertex and ridge. On the other hand, Pd occupied the vertex and ridge of PdRuIr alloy NP, while the facet is mainly occupied by Ir atoms. With an increase in temperature, the Ru composition and Pd−Ru bond fractions in the surface layer increase. The contribution of the configurational entropy effect overcomes the energy loss due to the Ru surface atoms. The surface energy becomes higher in the order of Au, Pd, Ir, and Ru. While the value of the surface energy for Pd is close to that for Au, that for Ru is close to that for Ir. With an increase in the temperature, elements with similar surface energies become mixed at positions with similar CN. The surface energy of the constituent elements is important to determine the stable configuration and temperature variation. On the other hand, the bond affinity of the constituent elements is different according to the phase diagram of Pd−Ru and Pt−Ru alloys. In the Pd−Ru alloy, Pd and Ru are immiscible with the exception of Pd- and Ru-rich composition [42]. Conversely, the Pt−Ru alloy forms the SS at any composition although the structure varies depending on the composition [43]. Therefore, while the Pd−Ru bonds in the PdRuIr alloy NPs cause energy loss, the Pt−Ru bonds in PtRuIr alloy NPs lead to energy gain. This study compared the temperature dependence of the stable configuration between PdRuIr and PtRuIr alloy NPs. The order of the surface energy in PtRuIr alloy NPs is the same as that of PdRuIr alloy NPs. The influence of bond fractions may result in the difference between the stable configuration and the temperature variation of PdRuIr and PtRuIr alloy NPs.

We made the regression equation of the excess energy for fcc PtRuIr alloy NP before Wang-Landau sampling. We considered the NP model comprising 201 ( = N) atoms to calculate the excess energies of alloy NPs. Note that high-symmetry hexagonal close-packed (hcp) Ru₂₀₁ NP was not prepared. The cohesive energies were assumed to scale to approximately N^−1/3. We estimated the total energy of hcp Ru₂₀₁ NP using the cohesive energies from the previous study [44]. The compositional and configurational variations were considered. We calculated the excess energies of SS, subcluster segregated, and vertex segregated configurations. The SS states of ternary alloy NP models were evaluated by pairwise short-range order, while the Warren-Cowley parameter was used in the binary alloy NP models. In vertex segregated configurations, the vertex is assumed to be occupied by Pt atoms. When the number of Pt atoms is less than 24, Ir atoms occupied the vertex. In PtRuIr alloy NPs, the third element C ( = Pt, Ru, or Ir) was added into A−B ( = Ru−Ir, Pt−Ir, or Pt−Ru) binary alloy NP, where the A composition is equal to the B composition. Besides, the additional 40 configurations at the random composition were prepared to enhance the accuracy of SL, where the atoms were exchanged without the specific condition.

Figure 7 (a) shows the compositional dependence of the excess energy in the binary alloy NPs. The excess energies of SS Pt−Ru and Ru−Ir alloy NPs were negative and lower than that of the subcluster segregated configuration. In the bulk system, the SS configuration of Ru−Ir alloy is observed below the melting point, where the structure of the alloy depends on the composition [45]. On the other hand, the SS configuration of Pt−Ir alloy is observed at high temperatures [46]. In Pt−Ir alloy NP, the SS configurations show positive excess energy, while the excess energy of the subcluster segregated configuration is almost zero. The relation between the SS and subcluster segregated configurations is consistent with the formation of the SS in the bulk. The excess energies of vertex segregated configurations were lower than those of SS configurations. The occupation of low-CN atomic positions affects the stability of alloy NPs. Figure 7 (b) shows the excess energies of PtRuIr alloy NPs as a function of composition. When the C composition was low, the signs of excess energies of PtRuIr alloy NPs were the same as those of A−B binary alloy NPs. With an increase in C composition, the signs of the excess energies become opposite. Various trends of the excess energies for the composition are shown in PtRuIr alloy NPs.

Figure 7.

 Excess energies of (a) A−B binary and (b) A−B−C ternary alloy NPs as a function of composition (A, B, C = Pt, Ru, or Ir).

In the bulk model, the bond fractions were used as the descriptor. In the alloy NPs, the elements with low surface energy are located at low-CN atomic positions. The segregation of the composition in the surface layer from the overall composition was added. In this descriptor, the surface was divided into the vertex, ridge, and facet for the Ru and Ir composition. In the validation process, 64 and 192 configurations were selected as the test and training sets. Figure 8 shows one of the comparisons between the excess energies predicted by SL and calculated by DFT. The MAE in the test set is lower than 5.0 meV/atom for five cases. The excess energies predicted by SL matched those calculated by DFT. Table 2 shows the regression coefficients of the excess energy of PtRuIr alloy NPs. While the regression coefficients of the Pt−Ru and Ru−Ir bond fractions were negative, that of the Pt−Ir bond fraction was positive. The regression coefficients of the bond fractions have a trend similar to the excess energy of binary alloy NPs (Figure 7 (a)). The Ru and Ir surface segregations showed positive regression coefficients. The regression coefficients of Ru surface segregations are larger than those of the Ir surface segregations. These characteristics of regression coefficients for surface segregation are consistent with the order in the surface energy of the constituent elements [27, 47]. The obtained regression coefficients were used in the energy estimation of the Wang-Landau sampling.

Figure 8.

 Comparison of excess energy values determined by MRA and DFT calculation.

Table 2.  Intercept (β₀) and regression coefficients of bond fractions (β_1i) and surface segregations (β_2j) for the excess energy of the PtRuIr alloy NP model

β0		−0.0003
β1	Pt−Ru	−0.0042
	Pt−Ir	0.0945
	Ru−Ru	−0.0120
	Ru−Ir	−0.0819
	Ir−Ir	−0.0082
β2	Ru (vertex)	1.0602
	Ru (ridge)	0.9381
	Ru (facet)	0.7319
	Ir (vertex)	0.8732
	Ir (ridge)	0.6860
	Ir (facet)	0.5396

Figure 9 shows the overall composition and temperature dependence of bond fractions and compositions in the surface layer of PtRuIr alloy NPs, where the probability distribution of the composition in the surface layer of Pt₆₇Ru₆₇Ir₆₇ alloy NPs at 0, 1000, and 2000 K was Figure 10. The surface energy becomes higher in Pt, Ir, and Ru. The Pt atoms occupied the vertex and ridge, while the facet was mainly occupied by Ir atoms. The Ru atoms were located in the core. In the surface layer, the Pt−Ir bond fraction was higher than the Pt−Ru and Ru−Ir bond fractions. With an increase in the temperature, the Ru composition in the facet increase. The Ru surface atoms resulted in the formation of the Pt−Ru bonds in the surface layer. On the contrary, the Pt−Ir bonds were deformed with an increase in temperature. The contribution of the configurational entropy effect overcomes the energy loss due to the Ru surface atom.

Figure 9.

 Temperature and compositional dependence of structural parameters in PtRuIr alloy NPs.

Figure 10.

 Probability distribution of the composition in the surface layer of Pt₆₇Ru₆₇Ir₆₇ alloy NPs at 300, 1000, and 2000 K. Ternary plot represents the probability of Pt, Ru, and Ru.

We compared the stable configuration and the temperature variation in PtRuIr alloy NPs with those in PdRuIr alloy NPs. The details of the regression equation for PdRuIr alloy NPs are shown in Ref 27. Figure 11 shows the overall composition and temperature dependence of bond fractions and compositions in the surface layer of PdRuIr alloy NPs. The order of the surface energy in PdRuIr alloy NPs is the same as that in PtRuIr alloy NPs. Thus, the trends of composition and bond fractions in the surface layer of PdRuIr NPs are similar to that of PtRuIr NPs. However, the temperature dependence of PdRuIr alloy NP is small compared to PtRuIr alloy NPs. In the surface layer, the Pd−Ru bond fractions are lower than the Pt−Ru bond fractions, while the Pd−Ir bond fractions are higher than the Pt−Ir bond fractions. The Ru composition in the facet of PdRuIr alloy NPs were lower than those of PtRuIr alloy NPs. The key factor is the characteristics of Pd−Ru and Pt−Ru bonds as mentioned above. In the PdRu alloy, the SS configuration is observed in unique composition and high temperature. The excess energies of PdRu alloys were positive [27]. On the contrary, the excess energies of PtRu alloy NPs were negative as shown in Figure 7. The energy loss is caused by the Ru surface atoms. The Ru surface atoms form the Pd−Ru and Pt−Ru bonds when Pt atoms occupied the low-CN position. The energy loss due to the Ru surface atoms is suppressed by the formation of Pt−Ru bonds and is enhanced by the formation of Pd−Ru bonds. The combinations of the constituent elements give a significant influence on the stable configuration and the temperature dependence of alloy NPs. In the bulk state, the ordered configuration was observed, when the excess energy is negative. On the other hand, the ordered configuration was not observed in the considered ternary alloy NPs even if the excess energies of the two combinations were negative. It indicates that the contribution of the surface may be larger than that of the bond.

Figure 11.

 Temperature and compositional dependence of structural parameters in PdRuIr alloy NPs.

The algorithm used in this study was combined with the DFT calculation and MRA. Specifically, the regression equation was used in the energy estimation of the sampling. If the regression equations are obtainable, the method can be applied to various materials. In this study, we treated the binary alloys in the bulk and ternary alloy NP model. With an increase in the constituent elements, the descriptors of the bond fraction increased. In the NP model, the descriptors for the composition of the surface (or CN) were added to the regression equation. For the descriptors of the multinary alloy NP, the kind of bond fractions and the composition of CN may increase. In the complicated structure, additional descriptors may be required. For example, the A−B bond fractions in the distorted system are divided into the longer- and shorter-bond types. The successful preparation of the descriptor is the key factor for the application to other systems.

4 CONCLUSION

We implemented the Wang-Landau algorithm to probe the stable configurations of the alloy NPs. In this algorithm, the configurational DOS in the multicanonical Monte Carlo algorithm is asymptotically obtained, which results in the partition function and free energy. We combined the energy estimation in the sampling with the DFT calculation and MRA. If the regression equation can be obtained, our idea is applicable to an increase in the element and the system variation. We applied this method to the bulk model of binary alloy and the NP model of ternary alloy. For the subject, the procedure of the DFT calculation, MRA, and Wang-Landau algorithm are shown. The stable configuration and temperature variation of alloy in the bulk and alloy NP can be calculated. The formation of SS alloy in bulk state occurs when the entropy effect contribution overcomes the energy loss or gain due to A−B bonds, and the transition temperature varies with the bond affinity of the constituent elements. In the case of alloy NP, on the other hand, the contribution of the surface plays a significant role in determining the stable configuration. At low temperatures, the low surface energy element prefers to occupy the low CN atomic position. As the temperature increases, the entropy effect contribution overcomes the energy loss due to the transfer of high surface energy elements to the surface layer, leading to a change in the stable configuration. The energy loss depends on the bond affinity of constituent elements. The energy gain due to A−B bonds facilitated the change in stable configuration due to the entropy effect of the alloy NP. Therefore, our approach could capture the trend of thermodynamics stability for various constituent element combinations, and it revealed that the stable configuration and temperature variation are affected by the surface energy and bond affinity of the constituent elements.

Acknowledgments

This work was supported by ACCEL, JST (JPMJAC1501), and JSPS KAKENHI Grant Number 20H05623. The computation was performed using MASAMUNE-IMR at Center for Computational Materials Science, Institute for Materials Research, Tohoku University.

REFERENCES

• [1]   K. Kusada, D. Wu, H. Kitagawa, Chem. Eur. J., 26, 5105 (2020). , doi:10.1002/chem.201903928 PMID:31863514
• [2]   D. Wu, K. Kusada, T. Yamamoto, T. Toriyama, S. Matsumura, I. Gueye, et al., Chem. Sci. (Camb.), 11, 12731 (2020). , doi:10.1039/D0SC02351E PMID:34094468
• [3]   D. Wu, K. Kusada, T. Yamamoto, T. Toriyama, S. Matsumura, S. Kawaguchi, et al., J. Am. Chem. Soc., 142, 13833 (2020). , doi:10.1021/jacs.0c04807 PMID:32786816
• [4]   H. Kobayashi, K. Kusada, H. Kitagawa, Acc. Chem. Res., 48, 1551 (2015). , doi:10.1021/ar500413e PMID:25993560
• [5]   K. Kusada, M. Yamauchi, H. Kobayashi, H. Kitagawa, Y. Kubota, J. Am. Chem. Soc., 132, 15896 (2010). , doi:10.1021/ja107362z PMID:20979361
• [6]   K. Kusada, H. Kobayashi, R. Ikeda, Y. Kubota, M. Takata, S. Toh, et al., J. Am. Chem. Soc., 136, 1864 (2014). , doi:10.1021/ja409464g PMID:24455969
• [7]   K. Sato, H. Tomonaga, T. Yamamoto, S. Matsumura, N. D. B. Zulkifli, T. Ishimoto, et al., Sci. Rep., 6, 28265 (2016). , doi:10.1038/srep28265 PMID:27340099
• [8]   K. Kusada, D. Wu, Y. Nanba, M. Koyama, T. Yamamoto, X. Q. Tran, et al., Adv. Mater., 33, 2005206 (2021). , doi:10.1002/adma.202005206 PMID:33751709
• [9]   W. L. Bragg, Proc. R. Soc. Lond. A Math. Phys. Sci., 145, 699 (1934). doi:10.1098/rspa.1934.0132
• [10]   J. Dean, M. J. Cowan, J. Estes, M. Ramadan, G. Mpourmpakis, ACS Nano, 14, 8171 (2020). ,doi:10.1021/acsnano.0c01586 PMID:32515581
• [11]   Y. Yao, Z. Liu, P. Xie, Z. Huang, T. Li, D. Morris, et al., Sci. Adv., 6, eaaz0510 (2020). , doi:10.1126/sciadv.aaz0510 PMID:32201728
• [12]   R. Kikuchi, Phys. Rev., 81, 988 (1951). doi:10.1103/PhysRev.81.988
• [13]   R. Kiuchi, J. Chem. Phys., 60, 1071 (1974). doi:10.1063/1.1681115
• [14]   J. M. Sanchez, D. de Fontaine, Phys. Rev. B Condens. Matter, 17, 2926 (1978). doi:10.1103/PhysRevB.17.2926
• [15]   T. Mohri, M. Ohno, Y. Chen, J. Phase Equilibria Diffus., 27, 47 (2006). doi:10.1361/105497106X92790
• [16]   D. S. R. Rocabado, Y. Nanba, M. Koyama, Comput. Mater. Sci., 184, 109874 (2020). doi:10.1016/j.commatsci.2020.109874
• [17]   Y. Nanba, M. Koyama, Comput. Mater. Sci., 203, 111132 (2022). doi:10.1016/j.commatsci.2021.111132
• [18]   B. A. Berg, T. Neuhaus, Phys. Lett. B, 267, 249 (1991). doi:10.1016/0370-2693(91)91256-U
• [19]   B. A. Berg, T. Neuhaus, Phys. Rev. Lett., 68, 9 (1992). , doi:10.1103/PhysRevLett.68.9 PMID:10045099
• [20]   U. H. E. Hansmann, Y. Okamoto, J. Comput. Chem., 14, 1333 (1993). doi:10.1002/jcc.540141110
• [21]   U. H. E. Hansmann, Y. Okamoto, F. Eisenmenger, Chem. Phys. Lett., 259, 321 (1996). doi:10.1016/0009-2614(96)00761-0
• [22]   N. Nakajima, H. Nakamura, A. Kidera, J. Phys. Chem. B, 101, 817 (1997). doi:10.1021/jp962142e
• [23]   F. Wang, D. P. Landau, Phys. Rev. Lett., 86, 2050 (2001). , doi:10.1103/PhysRevLett.86.2050 PMID:11289852
• [24]   D. P. Landau, F. Wang, Comput. Phys. Commun., 147, 674 (2002). doi:10.1016/S0010-4655(02)00374-0
• [25]   S. H. Tsai, F. Wang, D. P. Landau, Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 75, 061108 (2007). , doi:10.1103/PhysRevE.75.061108 PMID:17677221
• [26]   D. P. Landau, F. Wang, S. H. Tsai, Comput. Phys. Commun., 179, 8 (2008). doi:10.1016/j.cpc.2008.01.004
• [27]   Y. Nanba, M. Koyama, Bull. Chem. Soc. Jpn., 94, 2484 (2021). doi:10.1246/bcsj.20210199
• [28]   G. Kresse, J. Furthmüller, Phys. Rev. B Condens. Matter, 54, 11169 (1996). , doi:10.1103/PhysRevB.54.11169 PMID:9984901
• [29]   G. Kresse, J. Furthmüller, Comput. Mater. Sci., 6, 15 (1996). doi:10.1016/0927-0256(96)00008-0
• [30]   J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett., 77, 3865 (1996). , doi:10.1103/PhysRevLett.77.3865 PMID:10062328
• [31]   P. E. Blöchl, Phys. Rev. B Condens. Matter, 50, 17953 (1994). , doi:10.1103/PhysRevB.50.17953 PMID:9976227
• [32]   G. Kresse, D. Joubert, Phys. Rev. B Condens. Matter, 59, 1758 (1999). doi:10.1103/PhysRevB.59.1758
• [33]   J. M. Cowley, Phys. Rev., 77, 669 (1950). doi:10.1103/PhysRev.77.669
• [34]   D. de Fontaine, J. Appl. Cryst., 4, 15 (1971). doi:10.1107/S0021889871006174
• [35]   I. Karakaya, W. T. Thompson, Bull. Alloy Phase Diagr., 7, 362 (1986). doi:10.1007/BF02873023
• [36]   S. N. Tripathi, S. R. Bharadwqj, J. Phase Equilibria, 15, 208 (1994). doi:10.1007/BF02646369
• [37]   S. N. Tripathi, S. R. Bharadwaj, M. S. Chandrasekharaiah, J. Phase Equilibria, 12, 603 (1991). doi:10.1007/BF02645078
• [38]   H. Okamoto, T. B. Massalski, Bull. Alloy Phase Diagr., 6, 229 (1985). doi:10.1007/BF02880404
• [39]   Y. Nanba, M. Koyama, ACS Omega, 6, 3218 (2021). , doi:10.1021/acsomega.0c05649 PMID:33553938
• [40]   K. E. A. Batista, J. L. F. Da Silva, M. J. Piotrowski, J. Phys. Chem. C, 122, 7444 (2018). doi:10.1021/acs.jpcc.7b05714
• [41]   P. C. D. Mendes, S. G. Justo, J. Mucelini, M. D. Soares, K. E. A. Batista, M. G. Quiles, et al., J. Phys. Chem. C, 124, 1158 (2020). doi:10.1021/acs.jpcc.9b09561
• [42]   S. N. Tripathi, S. R. Bharadwaj, S. R. Dharwadkar, J. Phase Equilibria, 14, 638 (1993). doi:10.1007/BF02669157
• [43]   H. Okamoto, J. Phase Equilibria Diffus., 29, 471 (2008). doi:10.1007/s11669-008-9366-0
• [44]   Y. Nanba, T. Ishimoto, M. Koyama, J. Phys. Chem. C, 121, 27445 (2017). doi:10.1021/acs.jpcc.7b08672
• [45]   H. Okamoto, J. Phase Equilibria, 13, 565 (1992). doi:10.1007/BF02665768
• [46]   S. R. Bharadwaj, S. N. Tripathi, M. S. Chandrasekharaiah, J. Phase Equilibria, 16, 460 (1995). doi:10.1007/BF02645356
• [47]   R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson, et al., Sci. Data, 3, 160080 (2016). , doi:10.1038/sdata.2016.80 PMID:27622853