Over the past decades, bimetallic nanoparticles have attracted continuing interests because of their unique catalytic, optical, and electronic properties. Their physical and chemical properties can be widely tuned by the size, composition, and the chemical element ordering [1-10]. For example, for a bimetallic A_mB_n nanoparticle, it could be uniform alloy, A-/B-surface segregated, or core-shell configuration, which directly affects the catalytic activity. Therefore, it is much pursued to rationally design and synthesize more superior bimetallic clusters. However, the atomic-scale structural and compositional information on the bimetallic clusters has been less understood and also remains a challenge to obtain in both the experimental and theoretical exploration, owing to ther extremely complex configuration space. In particular, among the current bimetallic systems, platinum-based clusters, such as Pt-Cu, Pt-Co, Pt-Fe and Pt-Ni, have been widely studied to reduce the amount of expensive Pt in diverge applications [11, 12], in which the Pt-Au bimetallic clusters has received great attention in many theoretical and experimental studies. Nevertheless, toward the controlled synthesis of Au-Pt nanoalloy, the reactivity and stability dependence on its composition and atomic arrange-ment need further exploration [5, 13-16]. Specifically, there are two critical questions to be answered for the bimetallic clusters: (ⅰ) For a specific composition of A_mB_n, what is the atom arrangement pattern at the thermodynamically stable state? (ⅱ) What is the composition (or A/B atomic ratio) for the thermodynamically most stable A_mB_n cluster given a specific size?

Unfortunately, even for small-sized bimetallic clusters, the global search for the stable atomic arrangement could not be performed efficiently at the first-principles density functional theory (DFT) level. For a N-atom bimetallic nanoparticle A_mB_N-m even at a given structural framework, the number of isomers is huge at the order of N!/(m!(N-m)!). Taking Pt₇Au₃₁ as an example, the isomer of Pt₇Au₃₁ is over than 10⁷ (if neglecting the symmetry). Also, those isomers could share structure but differ in composition, i.e., forming the so-called homotops, which further dramatically increase the complexity. Therefore, computational exploration of the potential energy surface (PES) and searching the global minima of the PES in bimetallic clusters would be much more difficult in comparison with the unary ones [15]. Currently, the empirical atomistic potentials, e.g., the embedded-atom method and Gupta potential, were usually employed but suffered from inadequate accuracy [17]. Overall, an efficient and accurate description of bimetallic clusters' stability is highly desirable.

Enlightened by the recent development of machine learning in building inference models that could learn the relationship among the composition, structures and properties of chemical materials based on materials databases [18-23], here we explored an advanced artificial neural network framework (SchNet) [20] coupled with the genetic-algorithm (GA) to systematically study the huge chemical space of bimetallic Pt-Au clusters. Significantly, a machine learning model was constructed based on bimetallic cluster dataset and iteratively trained under the guidance of GA. The trained neural network model enabled us to efficiently predict the PES of the whole space of Pt-Au bimetallic clusters and quantitatively identify the thermodynamically stable component and atoms' distribution.

Genetic algorithm: To explore the complex configuration space of Au-Pt bimetallic nanoparticles, we adopted a global search technique, i.e., the genetic algorithm, which is stochastic and usually allows us to search multidimensional parameter space more effectively than deterministic approaches such as hill climbing [24, 25]. In the present work, we carried out self-built GA according to the following steps:

(Ⅰ) Structure encoding. Each individual in GA, corresponding to a specific Pt_mAu_38-m configuration, is encoded by binary vector as P = c₁c₂ ··· c_N, in which c_i presents a specific type of atom; here 0 and 1 represents Pt and Au atom, respectively, as demonstrated in Scheme 1.

(Ⅱ) Constructing fitness function. The fitness function is defined as

where E_mixing is the formation energy of Pt_mAu_38-m cluster by reference to the pure Pt₃₈ and Au₃₈ cluster, F is the object fitness. The more stable Pt_mAu_38-m is, the more negative F is. Notably, F is able to describe the thermodynamic stability of Pt_mAu_38-m cluster with variable composition, and m_i is the number of Au for a specific-composition Pt_mAu_38-m and Cis the penalty coefficient; for example, to explore the invariable-composition clusters of Pt_miAu_38-mi (i.e., with a specific number of Pt of m_i), we set C to be extremely large (> 100 here), otherwise C = 0.

(Ⅲ) Roulette wheel selection. To effectively improve the population diversity and prevent the premature convergence of GA, the roulette wheel selection is employed to generate new generation. The probability of an individual can be computed by the fitness. Following Eq. 3, the individual (ⅰ) is selected for crossover if d_i < r, otherwise it will be eliminated. r is a randomly generated number belonging to [0, 1].

where F(Pi) is the fitness of individual P_i, and N_p is the number of individuals in the population.

(Ⅳ) Crossover and mutation. The single-point crossover is introduced to generate new children individuals. For example, as shown in Scheme 1, the binary code of 'parent1' and 'parent2' are chosen for crossover, then two new individuals are generated. In addition, an atom exchange operation is applied to keep the atomic composition ratio unchanged. Also, the mutation operator is used to increase the population diversity.

(Ⅴ) Workflow. We started the GA structural search with an initial population at a seed size of N_p= 200 randomly generated, i.e., 200 individuals, and then repeated step (Ⅰ-Ⅳ) until achieving the termination conditions. Finally, the individual with the highest fitness is the optimal structure. In our calculation, the probabilities of crossover and mutation were set 0.8 and 0.002, respectively.

DFT Calculation: In this work, all the first-principles-level calculations were performed using density functional theory (DFT) in the VASP code [26, 27]. The generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE) was employed to describe the exchange-correlation potential. The ion-electron interactions are described by the projector augmented wave (PAW) method. The 38-atom cuboctahedral clusters were calculated in a large 20×20×20 Å supercell, and a cutoff energy of 450 eV and a k-point grid of 2×2×2 were used.

Machine learning model (SchNet): To fit the PES of 38-atom Pt_mAu_38-m clusters to reduce the usage of the time-consuming DFT calculations, an advanced atomistic neural network (SchNet) was fully examined and adopted in this work. SchNet is end-to-end model, and has achieved excellent performance for predicting molecular properties on the QM9 database and for the predicting energies and forces on the MD17 database [20]. Here, we, for the first time, applied it into the bimetallic systems. The overview of SchNet architecture can be seen in Scheme 2, which is composed of atom embedding layer, interaction blocks and prediction blocks. The SchNet has been implemented in the open-source package (SchNetPack)[20, 21]. The primary components of the networks were as follows.

Atom embedding: In principle, a certain cluster with n toms is completely specified by its nuclear charges, Z = (Z₁, ···, Z_n) and the position of these nuclear charges, R = (r₁, ···, r_n). At the beginning of neural network, all atoms are represented by a tuple of features X^l = (x₁^l, ···, x_n^l), in which each x_i^l is a real-value vector of dimension F(= 128 in this work) describing the features of atom i. F is the feature map dimension, and x_i^l is the representation of atom i at the l^th NN-layer. The values of x_i⁰ is initiated by the embedding layer x_i⁰ = a_Zi, which is set randomly at the start of training and would be optimized during training.

Atomwise layer: Atomwise layers are fully connected layers applied to each of the atomic representations, x_i^l, defined as x_i^l+1 = W^lx_i^l+b^l. Different atoms share the same weights (W^l) and biases (b^l), which keep scalable with respect the size of clusters. Also, the weights and biases will be optimized during training.

Interaction block: To incorporate the influence of neighboring atoms, continuous-filter convolutions are employed, which are defined as follows:

where ⊙ is element-wise multiplication and nbh(ⅰ) are the neighbors of atom i. W^l(rj-r_i) is a filter-generating network between atoms i and j and will be mentioned briefly in the following.

Notably, the critical component of SchNet architecture over other convolutional neural network (CNN)-based approaches is the introduction of continuous filter convolutional layers (CFCONV) rather than the discrete filter convolutional layers commonly used in Euclidean data. As shown in Scheme 2, an interaction block is defined as an atomwise layer, followed by a continuous-filter convolutional layer, an atomwise layer, the application of the activation function to each feature element, and a final atomwise layer. The resulting interaction for atom i is then added to the representation of atom i. For the activation function, the original SchNet studies chose a shifted soft-plus function, f(x) = ln(1+ex)-ln(2), for all activations.

The schematic of the filter-generating network defined in Eq. (4) is shown in Scheme 2. The input of filter depends only on the distances between atoms, which gives the SchNet architecture invariance to both translations and rotations of the system. These distance are then used as input for a set of radial basis function (RBF), e_k, defined as

The number of radial basis functions, K (300), is a meta-parameter, and the center of each radial basis function, μ_k, is chosen such that the K radial basis functions are equally spaced between 0 and a distance cutoff, d_k (30.0 Å), another meta-parameter. In Eq. (5), γ is a scaling meta-parameter (10 Å ^-2). After the radial basis functions have been applied, a dense layer is used to align the output dimension equal to the feature dimension. The output would go through an activation layer (Softplus), then another dense layer and one more activation layer to obtain the final output of the filter-generating network. Finally, the blocks use a residual connection inspired by RseNet (Eq. (6)), enabling for a flexible residual mapping of atom i(v_i^l) in interaction layer l that incorporates interactions between atoms and feature maps, as seen in Scheme 2.

Overall, SchNet obtains a latent representation of atomistic system by first using an embedding layer to obtain feature X⁰. These features are then processed by L interaction blocks, which generate the latent representation X^L, which can be passed to the prediction block.

Prediction block: After interaction blocks, these representations are then processed by a prediction block to obtain the desired properties of clusters. Finally, a dense layer with an output dimension of one is applied, giving the final output value for each atom. After obtaining these final outputs, sum pooling or average is used to give the property of interest, which is similar to the method of Behler and Parrinello [28].

In this work, we used the recommended hyperparameters [20]. Six interaction blocks and the sum pooling in the prediction block were employed. The detail dimensions of vectors could be seen in Scheme 2. All models were trained with mini-batch stochastic gradient descent using the frequently-used optimizer (Adam) with mini-batches of 32 examples [20]. The initial learning rate of 0.001 was used for the training steps. The learning rate was reduced by exponential decay with a ratio of 0.96 every 100000 steps. The loss of the validation set was computed every 1000 batch gradient training steps for early stopping.

Scheme 3 illustrates the whole workflow to explore the configuration space of bimetallic clusters. Step 1: We carried out a certain amount of prior DFT calculations of different Pt_mAu_38-m configurations to construct the initial dataset (~4000 individual in this work), which were used to train the original SchNet until obtaining the desirable accuracy (DFT level). Step 2: We sampled clusters from existing dataset as seeds to perform structural evolution with GA, in which the energy of each new configuration generated by crossover or mutation was provided by the trained SchNet. In this work, the structural evolution was performed at a number of 1000 generations (200 seeds/generation) at least until that no more stable configuration appears within 50 consecutive generations. Step 3: The candidate stable clusters searched by GA were further validated by DFT calculation; the ones, whose energy errors between SchNet and DFT are unneglectable, would be feedback to the dataset for further training. Finally, repeating step 2 and 3 iteratively several times, which are all on-the-fly performed and thus efficient, the high-quality SchNet and globally best configuration would be determined (see discussions later). Notably, there may be overfitting for training the neural networks, and we performed early-stopping once the learning rate is less than the threshold value (10 ^-6).

In the proposed workflow, the configuration exploration and SchNet training is performed iteratively, which grantees the prediction ability and searching ability at less computational cost, thus greatly accelerating the exploration process. From the E_mixing comparison between SchNet and DFT (Fig. S1 in Supporting information), one can see that the energy prediction ability of SchNet can be promoted to the accuracy comparable with DFT level with four interactions. For example, at the beginning of iterations the SchNet-assisted GA search (SchNet-GA) indicates that Pt_mAu_38-m (m: 6, 7, 8) are possible candidates but over-estimates the E_mixing, which could result from the unbalanced data distribution; after four iterations, the SchNet-GA has already figured out the most stable cluster and accurately predicted the energy at a high level.

Specifically, it is worth mentioning that, to describe the PES of bimetallic clusters, SchNet was trained on scalable size of dataset and stopped at the criteria that the learning rate decays to the minimum value (10^-6). During the training process, data were randomly partitioned into training (80%), validation (10%), and test sets (10%). It is well-known that machine learning models may suffer from overfitting, which causes the models not generalized well to data outside the training set. The test set, independent of the validation set, can be used to enable the use of early stopping, in which the final trained parameters are taken to minimize the loss function of the validation set.

The overall energy predication accuracy of SchNet toward the E_mixing of the evolved Pt_mAu_38-m clusters was shown in Fig. 1a. It turns out that errors for training data, testing data, and validation data are nearly 0.05 eV for both MAE and RMSE (Fig. 1b). The training, validation and test data almost share the same error level, which indicates that the effect of over-fitting can be neglected and neural network have the excellent generalization ability. The learning curve during training SchNet was also illustrated in Fig. 1c, which indicates the SchNet model has converged. In other word, the trained SchNet neural network could enable us to efficiently explore the PES of bimetallic clusters combining with the GA global searching method.

With the genetic algorithm assisted by the trained SchNet, the complex configuration space of Pt_mAu_38-m at each composition (i.e., m = 1–38) was thoroughly searched. The located E_mixing for each fixed composition of Pt_mAu_38-m were listed in Fig. 2 (only the relatively negative values were shown); the most negative one of each composition, corresponding to the most stable configuration, constitute an upward parabolic curve as a function of the composition. Notably, the SchNet model exhibits a good agreement with DFT result for this parabolic curve.

For each composition of Pt_mAu_38-m cluster, the located thermodynamically most stable configuration is showed in Fig. 3, which clearly specifies the relative distribution of Pt and Au atoms. It was found that, when substituting an Au atom into Pt₃₈ clusters (m = 37), the Au atom would prefer to occupy the least-coordinated surface corner site; as the second Au atom enters, it also tends to be located at another surface corner site nearest the first Au atom. Interestingly, on the Pt-rich Pt_mAu_38-m cluster (m ≥ 7), a phenomena is that Au always forms the continuous patch-like shapes and keep at the surface layer (Fig. 3). In other word, there is a strong tendency for Au atoms to segregate to the surface of clusters as Au content increases.

Correspondingly, on the Au-rich Pt_mAu_38-m cluster, Pt atoms tend to occupy the core sites when mixing with less than six Pt atoms (m ≤ 6); when m = 6, it forms a perfect core-shell Pt₆Au₃₂ structure. When m ≥ 6, the extra Pt atoms are favorable to stay in the high-coordinated surface center site. For example, for the Pt₈Au₃₀ cluster, there are two Pt atoms are interestingly isolated at surface center sites besides the six Pt atoms at the core sites.

As Fig. 2 shows, from the parabolic curve between the E_mixing of the most stable configuration for each composition and the contained Pt number, one can see that as the number of Pt atom (m) increases from 1 to 6, the corresponding bimetallic clusters become more and more stable. At m = 6, i.e., the Pt₆Au₃₂, it form the thermodynamically most stable one with a perfect core-shell structure validated by DFT, in which all the core sites are Pt atoms and the surface are totally occupied by Au atoms. In other word, it indicates that the Pt prefer to keep in the core of clusters and Au prefer to distribute on the outer surface thermodynamically. When m > 6, more Pt involving in Pt_mAu_38-m would gradually weakens the relative stability. Noteworthily, when Pt is very rich, Pt_mAu_38-m (m = 34–37) exhibit even positive E_mixing values, which implies that they are less stable than their pure-phase Pt₃₈ and Au₃₈ clusters and thus low-content Au atoms is not favorable energetically mixed in the Pt metallic cluster.

Now several questions arise: Why do Pt atoms tend to preferentially occupy the core region and then the surface center site? Why do Au atoms prefer to be segregated at the surface layer and specifically the corner sites? Why is the core-shell Pt₆Au₃₂ thermodynamically most stable? To understand these interesting distribution behaviors, the structural features and the inherent mechanism were analyzed. Generally, within the 38-atom cuboctahedral cluster framework, each of three sites (core, center and corner) has its characteristic bonding environment, and the coordination number by counting the nearest metallic bonds are listed in Scheme S1 (Supporting information), corresponding to 12, 9 and 6 for the core, center and corner site, respectively. In other word, the core-region atom has the maximized coordination number, while the surface corner site is the least coordinated. Meanwhile, we estimated the bond strength of Pt-Pt, Au-Pt and Au-Au and the potential energy curve over different bond lengths (Fig. S3 in Supporting information), and it can be found that the bond strength always follows the order of Pt-Pt > Pt-Au > Au-Au within the bond length range spanning from 2.5 Å to 3.1 Å

Firstly, we used PtAu₃₇ as an example to illustrate the distribution behavior of Pt. The mixed atom (Pt) at different sites correspond to forming/destroying different number of Pt-Au bonds and Au-Au bonds. Considering the bond strength order of Pt-Pt > Pt-Au > Au-Au, it is obvious that Pt tend to form as many as bonds, thus preferentially occupying the core sites. Semi-quantitatively, by comparing Pt at the center versus corner site, or at the corner versus center site (Scheme S1), the corresponding E_mixing difference could be presented in Eqs. (7) and (8) by the first-order bond energy decomposition.

Thus, it can evidently give out that E_core-E_center and E_center-E_core is above 0, indicating that Pt is more stable at the core site, secondarily stable at the surface center site, and least stable at the surface corner site. Alternatively, by exploring Pt₃₇Au, it could be understood that Au would keep at least-coordinated corner site to make room for forming the most Pt-Pt bonds.

To deeply understand the dependence of atom distribution on the stability, the Pt₆Au₃₂ cluster was further employed as example. Its 28 typical configurations were selected and the corresponding structures can be seen in Fig. S4 (Supporting information). For each configuration, the corresponding bond-type number was summa-rized in Fig. 4a, which basically illustrates that the more stable it is, the less (more) the Au-Au (Pt-Pt) bonds are.

From the correlation between E_mixing and the number of different bond types depicted in Figs. 4b–d, one can interestingly see that the Au-Au bond number among those configuration has more close correlation with E_mixing in comparison with the other Pt-Pt and Pt-Au bonds, corresponding to R²= 0.82, 0.69 and 0.46, respectively. As it has been demonstrated that Au-Au is the weakest among the above three bond types, by referring to traditional cask theories, the number of Au-Au bond as the weakest part should largely decide a cluster's stability trend, which may qualitatively rationalize the higher correlation (comparing with that of Pt-Au and Pt-Pt) in Fig. 4.

More importantly, a new descriptor (the difference of Au-Au and Pt-Pt bond number) was found to exhibit excellent linear correlation with E_mixing. Specially, when considering a scaling coefficient α(= 0.50, roughly estimated from Fig. S3) that reflects the ratio of Au-Au and Pt-Pt bond strength, it gives a R²= 0.96. Therefore, it could be concluded that the most stable configuration with a fixed composition should minimize the number of weakest Au-Au bond and maximize the number of strongest Pt-Pt bond. In this sense, it is obvious that Au atoms tend to keep at the surface sites (9 and 6 bonds) while Pt stay in the core region (12 bonds).

Furthermore, it may be worth discussing that the origin determining the relative distribution of Pt and Au atom within these Pt_mAu_38-m could also be perceived from the SchNet itself. Beyond predicting E_mixing, the atomistic neural network architec-ture of SchNet enables us to readily get access to atom-wise latent variables, which provides a framework for atom-wise explanation out-of-the-black-box. The atom-wise saliency can be seen as the logical extension of the pixel-wise explanations used in the "imaging" strategy for describing the cluster structure [29]. An aggregation of these latent variable constitutes the final object property of interest, and thus these "latent variable" get assigned an inherent physical meaning. Here, with respect to the predicted E_mixing, these "latent variable" correspond to the energy contribu-tion relating to each atom involved, i.e., E_mixing =ƩPt_i + ƩAu_i. Fig. S5 (Supporting information) depicts the distribution of different atoms' energies obtained from the SchNet model. It could be seen that total Pt's contribution (ƩPt_i) basically determines and classifies the energy level of E_mixing in comparison with the total Au's contribution (ƩAu_i), as shown by the contour plot of E_mixing varying with ƩPt_i and ƩAu_i in Fig. S5a. Also, it could be seen from Fig. S5b that Pt's contribution is greater than Au's for most of the Pt_mAu_38-m with different compositions (the more negative the value, the greater the energy contribution). Significantly, these energy contribution of different elements (Pt, Au) obtained by SchNet is eventually in coincidence with the bond strength order (Pt-Pt > Au-Pt > Au-Au), demonstrating that the machine learning could automatically learn the potential nature of bonding mechanism. Thus, the neural network could give an excellent description on the potential energy surface of bimetallic clusters.

In summary, we have developed an approach by combining the machine learning model (SchNet neural network) and genetic algorithm to tackle the challenge problem for the high-throughput screening for the stable bimetallic clusters and determining the atomic distribution from a huge configuration space. The method enables us to quickly predict the energy of 38-atom Pt_mAu_38-m cluster at the DFT level and identify the most stable composition and the preferential relative distribution rule of Pt and Au atoms, respectively. The following significant results were obtained.

(ⅰ) By learning from ~15, 000 bimetallic cluster iteratively, the most stable bimetallic cluster was found and the correspond-ing atoms' distribution was identified. Thermodynamically, the most stable Pt_mAu_38-m cluster is the core-shell Pt₆Au₃₂, with all the core sites occupied by Pt.

(ⅱ) The relative distribution of Au and Pt atoms within each composition of Pt_mAu_38-m (m = 1 – 38) were identified, and it was found that Au atoms prefer to occupy the surface layer and specifically the least-coordinated corner sites, forming a continuous patch-like surface Au segregation; by comparison, Pt tend to accumulate at the high-coordinated core region.

(ⅲ) The stability trend of Pt_mAu_38-m was largely rationalized by the variation of bond number of Pt-Pt, Pt-Au and Au-Au. The number difference between the Pt-Pt bond and Au-Au bond can serve as a good descriptor to correlate the stability of the corresponding cluster; the more the strong Pt-Pt bond and the less the weak Au-Au bond, the more stable the Pt-Au cluster.

Overall, this is the first demonstration of SchNet applied to the PES exploration of bimetallic nanoparticles, which largely accelerate the discovery of bimetallic clusters, and the given global searching methods assisted by data-driven machine learning models could be generalized to other types of systems.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This project was supported by National Key R & D Program of China (No. 2018YFA0208602), NSFC (Nos. 21622305, 21873028, 21703067), National Ten Thousand Talent Program for Young Top-notch Talents in China, Shanghai ShuGuang project (No. 17SG30).

Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.cclet.2019.12.006.