Global minimization of gold clusters by combining neural network potentials and the basin-hopping method†

Neural network potentials trained by first-principles density functional theory total energies were applied to search for global minima of gold nanoclusters within the basin-hopping method. Using Au58 as an example, we found a new putative global minimum which has a core–shell structure of Au10@Au48 and C4 symmetry. This new structure of Au58 is 0.24 eV per formula more stable than the best previous model that has C1 symmetry. This work demonstrates that neural network potentials combined with the basin-hopping method could be very useful in global minimization for medium-sized metal clusters which might be computationally prohibitive for first principles density functional theory.

Gold nanoclusters and nanoparticles have attracted great research interest due to their novel physical and chemical properties in contrast with their bulk counterpart. Pioneered by Haruta, catalysis by nanogold has been one of the most active research areas in heterogeneous catalysis. [1][2][3] Understanding the structures of gold nanoclusters and nanoparticles can help elucidate their catalytic actives. Because of their small sizes, gold nanoclusters and nanoparticles are usually stabilized by an oxide support or a ligand monolayer. 2,[4][5][6] Synthesis and structure determination of ligand-protected gold nanoclusters have advanced greatly in the past decade, 6-10 while structural understanding of oxide-supported gold nanoparticles has been hampered by their complexity (size dispersion, metal-support interaction, oxygen vacancies, etc.). 2,11,12 Au clusters in the gas phase, on the other hand, offer a simpler and interesting playground for both theorists and experimentalists to explore the evolution of their intrinsic structures and properties with size. 13 For example, it was found that very small Au n clusters with n < 13 for anion, 14 <8 for cation, 15 and <12 for neutral 16 adopt planar configurations; clusters from Au 16 to Au 18 have cage configuration, while Au 19 and Au 20 have a pyramid shape. [17][18][19][20][21] From Au 21 to Au 35 the structure evolves from pyramidal to tubular, and then to core-shell structure. 13 The global minimum of Au 40 has been predicted by first principles basin hopping to have a twisted pyramid structure with a tetrahedral Au 4 core. 22 One of the largest Au clusters explored by DFT is Au 58 and the putative global minimum was found to have a double-shell structure of Au 10 @Au 44 of C 1 symmetry. 23 It becomes increasingly difficult to do global minimization for larger size Au clusters (say, n > 50) using first principles methods such as density functional theory (DFT) due to the exponential increase of the local minima number with size 24,25 and the nonlinear scaling of the computational cost (roughly ∼n 3 for DFT). In this context, alternative methods such as empirical potentials were often employed, [26][27][28][29][30] such as the Rosato-Guillope-Legrand potentials for Au clusters up to 318 atoms, 26 Sutton-Chen embedded atom potentials for Pt 55 and Au 55 , 27 and effective-medium-theory potentials for Au 147 and Au 309 . 30 Most empirical potentials, however, are not accurate enough, and sometimes could lead to unreliable predictions. 13,25,[27][28][29] One promising way to overcome this difficulty is to use the artificial neural network (NN) potentials, trained by large data sets of first principles total energies. As a major approach of machine learning, an artificial NN is a group of interconnected nodes mimicking how neurons in the brain work. In the context of NN potentials for a chemical system (for example, a nanocluster), the goal is to construct a parameterized analytical expression for the potential-energy surface (PES) of the system by using an artificial NN. The details about how to construct such a NN PES based on DFT datasets can be found from the recent literature. [31][32][33] The NN potentials combine the advantages of the speed of empirical potentials and the accuracy of the first principles methods. A number of recent papers have already shown the success of NN potentials in geometry optimization and molecular dynamics of large systems, such as free standing metal particles, oxide supported metal particles, and solution systems. 32,34,35 However, NN-potentialbased global minimization has not been demonstrated.
In this work, we apply the DFT-trained NN potentials to the global minimum search of Au nanoclusters for the first time. The Au 58 was chosen here because it was reported to be highly stable and robust as observed from experiment. 36,37 The basinhopping method has been demonstrated to be a highly efficient algorithm for global minimization of clusters. 38,39 Therefore, the basin-hopping method combined with the NN potentials (the NN-BH method for short) was employed in this work for global minimization of Au 58 . Stable structures found by NN-BH were then validated by DFT. For DFT calculations, we employed the Vienna Ab Initio Simulation (VASP) code, 40 with projector augmented wave 41 for the description of core-valence electronic interaction, Perdew-Burke-Ernzerhof 42 (PBE) functional for electron exchange-correlation, and plane wave for the basis set. Unless mentioned otherwise, calculations were done in normal precision with a cubic box of size 20 × 20 × 20 Å for the Au 58 . Only Γ-point was used to sample the Brillouin zone. Conjugate-gradient algorithm was employed for geometry optimization with force convergence criterion of 0.03 eV Å −1 .
The flow chart to construct accurate NN potentials is shown in Fig. 1. The construction was started from the generation of 400 random initial structures of Au 58 , which were created by randomizing the 58 Au atoms within a sphere of certain radius and with a constraint of minimal Au-Au bond length of 2.4 Å. These random structures were then optimized by DFT calcu-lations, and the structures and energies during all the relaxations (totally 23 358 configurations spanning 10 eV) were then used as initial reference data to train the NN potentials. A feed-forward NN was used which is a nested sum of activation functions and contains many weight parameters. [31][32][33] Construction of NN potentials is to fit the weight parameters to known reference potential surfaces. The DFT data set was split randomly into 80% for the training and 20% for the test. The NN architecture of 56-30-30-1 (one input layer of 56 nodes, two hidden layers of 30 nodes each, and one output layer of 1 node which is just the energy) was adopted, which is reasonable to avoid inadequate or over fitting. 32,35 The constructed NN potentials were then used to calculate energy and force of given structures. To further refine the NN potentials, we used the NN-BH method to generate some new structures. Their energies were then checked by DFT: if the NN energies were in poor agreement with the DFT energies, these new structures would then be incorporated into the reference data set to refine the NN potentials. The final adopted NN potentials have a root mean square error (RMSE) of 0.60 meV per atom for the training set and 0.66 meV per atom for the test set, indicating that the quality of the NN potentials is acceptable. Fig. 2 shows the validation of the NN potentials on new structures from BH. Although some deviation can be as large as 0.6 eV, the overall agreement is good between NN potentials and DFT.
Armed with the NN potentials, we applied the BH method to search the global minimum of Au 58 . 20 NN-BH jobs (that is, Monte-Carlo walkers) were started from 20 random initial structures, each running for 6500 Monte-Carlo steps. The energy landscape from the 20 NN-BH jobs totaling 120 000 Monte-Carlo steps is plotted together in Fig. 3. Fig. 4 shows an example NN-BH run. It can be seen that combining the BH runs and many random initial structures is an effective way to explore the complex configuration space of the Au 58 cluster. Fig. 3 labels the seven low energy isomers of (a)-(g), whose energy order has been checked by DFT (Table 1).  Isomer (a) was confirmed by DFT to be the most stable structure of Au 58 found from this work: we further optimized it with a tighter force tolerance (0.001 eV Å −1 ) and finite-difference analysis of the normal modes found no imaginary frequencies (see the ESI † for the coordinates and the calculated IR spectrum). The structure of isomer (a) is shown together with the other six low-energy isomers in Fig. 5. Isomer (a) has C 4 symmetry and is of a prolate spheroid shape with two squares at the top and bottom. It comprises a 48-atom outer shell and a 10-atom inner core. The other isomers also have a core-shell structure but they are less symmetric and less stable by at least 0.30 eV than isomer (a).
The HOMO and LUMO orbitals of the most stable Au 58 structure (Fig. 5a) are shown in Fig. 6. One can see that the HOMO orbital primarily distributes at the core and the waist of the shell, whereas the LUMO orbital populates the top and bottom of the outer shell. This indicates the spatial specific reactivity [43][44][45] of the Au 58 cluster, which could be interesting to explore in further studies. Isomer (a) has the largest HOMO-LUMO gap of 0.82 eV, while isomers (b)-(g) have gaps between 0.4-0.7 eV.
It is interesting to compare our best Au 58 structure with previous models. In an earlier theoretical work, Dong and Gong (DG) performed global minimum search using genetic algorithm based on empirical potential, and then examined the Fig. 3 The energy landscape of 20 BH runs based on NN potential. Adjacent BH runs are indicated by different colour (blue or black). (a)-(g) are the local minima with their structure shown in Fig. 5. Fig. 4 Energy landscape of an example single BH run based on NN potential with 6500 steps. The two red dots denote the initial random structure and the most stable structure found by the BH.  found structures by DFT. 23 Their best model also has the construction of a 10-atom core and 48-atom shell. We compared the energy between their structure and ours by DFT. We tested various functionals (PBE, 42 PBEsol, 46 and TPSS 47 ) and also employed a large box size of 30 × 30 × 30 Å in order to accurately evaluate the energy difference. We found that our structure is more stable than theirs for all functionals; for the TPSS functional which is known to be good for gold, our structure is 0.24 eV more stable. Geometrically, the two structures are quite similar, having a similar Au 10 core, but the Au 58 from this work is of C 4 symmetry, whereas the DG structure is of C 1 symmetry, which may explain the higher stability of our structure. The main difference lies in the shell structure as shown in Fig. 7. The outer shell of isomer (a) from this work has two squares at the two ends (top and bottom) of the prolate spheroid (Fig. 7a) and eight five-coordinated Au atoms symmetrically distributed on the shell (Fig. 7a and 7c). In contrast, the outer shell of the DG structure has a diamond at one end and a square at the other end (Fig. 7b); there are 10 five-coordinated Au atoms on the shell (Fig. 7b and 7d). We think that the higher symmetry of our structure leads to its higher stability. The electronic density of states (DOS) of the two structures near the Fermi level also shows small differences (Fig. 8): the higher symmetry of our structure leads to sharper peaks in the DOS. In summary, neural network (NN) potentials have been used to search for the global minimum of Au 58 in combination with the basin-hopping (BH) algorithm. NN potentials were trained and refined by DFT data sets. Then 20 NN-BH runs with 20 random initial structures were performed with a total of 120 000 Monte-Carlo steps. The putative global minimum of Au 58 was found to have C 4 symmetry and a core-shell configuration of a 10-atom core and 48-atom shell. Its stability was confirmed by DFT and found to be 0.24 eV more stable than the best previous model of a similar construction but of lower symmetry (C 1 ). This work shows that the neural network potentials can be used for global minimization of nanoclusters beyond the capability of the conventional first principles method.