Optimizing the Shape and Chemical Ordering of Nanoalloys with Specialized Walkers

New algorithms for the optimization of alloy nanoparticles (nanoalloys) are presented. The new algorithms are based on the concept of multiple basin‐hopping walkers running in parallel, each with its own specialized task—the flying walker, exploring the energy landscape at high temperatures to sample different geometric structures, and the landing and hiking walkers mainly refining the optimization of chemical ordering at low temperatures. These algorithms are referred to as flying‐landing (FL) and flying‐landing‐hiking (FLH). The algorithms are tested against several benchmarks (AuCu and AuRh clusters of 400 atoms and PtNi clusters of 38 and 55 atoms). In all cases, both FL and FLH are shown to perform very well compared to previous results in the literature. In addition, the algorithms are applied to the optimization of larger AgCu nanoparticles, with sizes up to 4000 atoms, in order to establish the behavior of the mixing energy and to compare full global optimization of shape and chemical ordering with optimization of chemical ordering alone at a fixed shape. In general, the results show that the simultaneous optimization of shape and chemical ordering is necessary in many cases, and that the FLH approach is especially efficient for that purpose.


Introduction
Metallic and multi-metallic nanoparticles (the latter being referred to as nanoalloys) have been widely studied in recent years for their surprisingly different properties from their bulk counterpart. [1] These differences at the nanoscale led to a wealth DOI: 10.1002/adts.202300268 of potential applications in important areas such as catalysis, [2,3] biomedicine, [4] and optics. [5,6] The experimental determination of many properties is very often aided by numerical simulations, which play a key role in the whole characterization of these systems. In fact, numerical simulations have many advantages. It stands out for the importance of its ability to examine a nanoparticle at the atomic level and on very short time scales-features that are not always possible in experiments. In fact, the experimental resolution now available is very often not sufficient for understanding the atomic-level mechanisms that determine the resulting structures of the clusters. In this regard, numerical simulations are very useful since they may allow a thorough understanding of these mechanisms. [7,8] A central aspect of the joint experimental-theoretical investigation of nanoparticles is the determination of their configurations of lowest potential energy. This is a fundamental step in the study of nanoparticle properties, [9] since the structure with the lowest potential energy corresponds to the thermodynamically stable configuration at low temperatures. This problem is known as global optimization, since it consists in determining the lowest minimum of a function of many variables-the potential energy surface (PES) which depends on all atomic coordinates.
The search of the global minimum of a function of many variables is a long-standing open problem. Due to its intrinsically complicated nature, there is no choice of attempting any analytical calculation; in practice, this problem can be tackled by numerical techniques only. The two key questions that any global optimization technique always needs to answer to, are the following: i) How broad is the exploration of the high-dimensional space in which the function is defined? ii) How fast is the convergence to the global minimum?
It is in fact crucial, for any algorithm that is in charge of finding the global minimum of any function, to be able to look for as many different relevant portions of that space as possible that is to perform a wide exploration of the energy landscape. [10] In addition, the algorithm must be able to go down deep in the relevant regions of the landscape, thus performing an efficient exploitation. [10] Several approaches have been recently proposed to deal with these issues. [11][12][13][14][15][16] It should be noted that in general, for question (ii), any global minimum found by any optimization is really worth that name only if it is possible to compare it to all the local minima of the function. Therefore, since it is virtually impossible to find all local minima of a function with more than a few variables in a reasonable amount of time, the global minimum found after a global optimization search is only a putative one.
The global optimization of a function of many variables has received notable attention due to the large variety of problems in which it is directly involved, [17][18][19][20] two notable examples beyond metal nanoparticles being the traveling salesman problem [21] and protein folding prediction. [22] Even focusing on the study of the global optimization of metal nanoparticles and nanoalloys only, there has been intense research activity and much effort, because of the importance of these systems both in fundamental science and in practical applications. [23][24][25][26][27] Different methods have been recently developed to tackle the problem of global optimization: basin hopping [28] and evolutionary-based algorithms, [29][30][31][32][33][34][35][36][37] lattice Monte Carlo, [38] mirror-rotation optimization algorithms, [39] curvilinear coordinate, [40] biminima optimizations, [41,42] and particle swarm optimization. [43] Global optimization is a very challenging task already for one-component (elemental) nanoparticles due to the enormous number of local equilibrium configurations (i.e., of local minima) of the PES, which is expected to exponentially increase with the number of atoms N. [44][45][46] These local minima correspond to different geometric structures which in principle can present quite different properties. The global optimization problem is even more complex for binary systems. In fact, in binary systems, not only the geometric shape has to be optimized, but also the chemical ordering. Chemical ordering is the pattern in which the two elements are arranged in the cluster or nanoparticle. Its optimization is challenging because of the enormous number of permutational isomers (often referred to as homotops [47,48] ).
In this article, we introduce two new algorithms specifically designed for the global optimization of bi-metallic nanoparticles. The new algorithms improve upon the basin hopping algorithm previously developed in refs. [49,50] which were favorably tested for the optimization of nanoalloys of sizes up to a few hundred atoms. [51] The basic principle behind these new algorithms is to perform the optimization process by means of different basin hopping walkers exploring the PES in parallel. Each walker is designed with its own specific features, in such a way that the combination of the walkers allows at the same time a broader exploration of the PES and a faster convergence to low-energy minima. Here we propose two algorithms based on this principle, the flying-landing (FL) algorithm and the flying-landing-hiking (FLH) algorithm. We show that these algorithms perform very well when compared to our previous multi-temperature basinhopping (BH) scheme [51] and to other algorithms proposed in the recent literature. [31] In addition, we clearly show that physically meaningful nanoalloy optimization requires in most cases the simultaneous optimization of shape and chemical ordering, while the optimization of chemical ordering with constraints on shape and/or symmetry can lead to misleading results in many cases.

Simulation Code
The simulation code [52] builds up on a standard Basin Hopping code and is able to perform the simulations by the new algorithms described below, together with the algorithms of refs. [49,50] that make use of multiple BH walkers interacting in a suitable order-parameter space. The code is written in C++ language. Besides binary systems, the code can optimize also multicomponent systems.
The set of elementary moves implemented in the code comprises shape-changing moves and atomic swap moves. [1,50] Among shape-changing moves, we recall: the shake move, in which all atoms are moved by randomly selected displacements within spheres centered in their present position; the Brownian move, in which the atoms are displaced by a short moleculardynamics run at high temperature, the bonds move, in which low-coordinated surface atoms are displaced to sites of higher coordination. More details on the settings of Brownian moves are given when treating the specific examples. In the bonds move, a list of the atoms with the lowest coordination n low (list 1) and a list of atoms of coordination n low + 1 (list 2) are made. An atom from list 1 is randomly chosen and it is displaced within nearestneighbor distance in the vicinity of another atom randomly chosen between those of list 1 and list 2. Exchange moves comprise both random exchanges, in which the atoms of the pair are randomly selected, and tailored exchanges, in which the atoms of the pair are chosen according to their local environment to better achieve either phase separation or intermixing. [53] This means that, for each exchange move, a list of A atoms (list A) and a list of B atoms (list B) are made. Then we select randomly one atom from list A and one atom from list B and the coordinates of the two atoms are exchanged. In random exchanges, the probability of choosing an atom from a list is uniform, while in tailored exchanges is weighted on the basis of the local environment of the atoms. [53] The parameters employed in the moves are specified when treating the specific case studies.

Search Algorithms
In the following, we describe the FL and FLH algorithms.
In the FL algorithm, there are two basin-hopping walkers, the flying walker (fw) and the landing walker (lw). The fw evolves by a majority of shape-changing moves [50] accepted at high temperatures and a smaller percentage of exchange moves, which are accepted at a lower temperature. Temperature values are given when treating the specific case studies. As shown in ref. [51], this is the most efficient scheme for a standard BH algorithm in the optimization of nanoalloys. On the other hand, the lw evolves by a majority of exchange moves and a smaller percentage of shapechanging moves, with all moves accepted at low temperatures. Examples of acceptance temperatures and percentages for different moves can be found in the Supporting Information. We note that the most efficient settings for BH searches in nanoalloys, as determined in ref. [51], imply acceptance rates in the range 20-40% for moves aiming at the wide exploration of the landscape, while the best refinement of chemical ordering (i.e., for www.advancedsciencenews.com www.advtheorysimul.com exploitation purposes [10] ) is obtained with much lower acceptance rates, below 5%.
The basic idea of the FL algorithm is that the lw periodically takes configurations from the fw (i.e., it is restarted) and refines them to go deep down in energy, while the fw continues its broad range exploration of the landscape at high temperatures. To do so, a recipe to periodically restart the lw is necessary. Here we propose the following scheme.
Let us consider a simulation of N BH steps, numbered by n = 0, 1, … , N − 1. With C fw (t n ) and C lw (t n ) we indicate the configurations of the fw and lw, respectively, whereas E fw (t n ) and E lw (t n ) are their locally minimized energies. Let t n 0 be the time of the last restart of the lw. At t n 0 both the fw and the lw have the same energy E(t n 0 ) and then they evolve independently. The two walkers are updated sequentially one by one, by some kind of move, according to percentages and acceptance temperatures as explained above. Now, let us assume that at t̃n > t n 0 the fw has improved its own best energy. When this happens, the lw may be restarted. One possibility is to restart the lw each time the fw improves its own best energy. More generally, the fw may be restarted by a probabilistic criterion. By definition, we have that and, since we expect that the landing walker is decreasing faster in energy than the flying walker, we have also so that the quantity satisfies 0 ≤ ≤ 1 in most cases. We now restart the lw with probability p r where ≥ 0, by extracting a random number r fl uniformly distributed between 0 and 1. If r fl ≤ p r the lw is restarted, and takes the configuration of the fw at t̃n, that is, C lw (t̃n +1 ) = C fw (t̃n). If the condition of Equation (2) is not fulfilled, the lw is restarted with probability 1. For = 0 the lw is restarted every time the fw finds a new lowest energy, since p r = 1 always. Other reasonable choices of are = 1∕2 or = 1. The number of restarts decreases at increasing . Figure 1 reports a block diagram of the FL algorithm. Representative behavior of the time evolution of the fw and lw walkers is shown in Figure 2, where = 1 is used, for an example of optimization of an Au 200 Cu 200 nanoalloy (the optimization of this system will be thoroughly treated in a following chapter). The fw evolves at high temperatures, finding local minima whose energies correspond to the orange line. Every time the fw improves its best local minimum (see the red squares), the algorithm decides whether to restart the lw according to the criterion of Equations (3) and (4). When the lw restarts, its energy suddenly increases (see the light blue line in Figure 2) and then quickly drops down thanks to the exchange and shape-changing moves accepted at low temperatures. Each time the lw finds a better local minimum, the structure is stored together with its energy (these events correspond to the purple circles).
A further refinement of the algorithm is obtained by adding a third walker, denoted as the hiking walker (hw), whose configurations and energies are C hw (t n ) and E hw (t n ). In the following, this algorithm is referred to as FLH algorithm. The hw is used for continuing the refinement of the configurations at low temperatures even when the lw jumps up in energy to restart its search. In order to do that, the hw uses the same type of moves and acceptance temperatures as the lw. Again, examples of acceptance temperatures and percentages for different moves can be found in the Supporting Information. The FLH algorithm works as follows. Let us assume that at time t̂n the lw reaches energy which is lower than the lowest energy reached by the hw walker so far, that is, E lw (t̂n) < E best hw . If this condition holds, the hw is restarted by taking the configuration of the lw and its energy, so that E hw (t̂n +1 ) = E lw (t̂n). At the same time, the lw is restarted too, that is, it takes configuration and energy from the fw, so that E lw (t̂n +1 ) = E fw (t̂n). In this way the lw goes up and down in energy between the fw and the hw, with the aim of feeding the hw with new good configurations to be further refined. Also in this case the three walkers are updated sequentially one by one, according to percentages and acceptance temperatures for different moves as explained above.
The block diagram of the FLH algorithm is given in Figure 3 with a representative behavior of the time evolution of the walkers is shown in Figure 4. The fw evolves in the same way as in the FL algorithm, whereas the lw is restarted in two cases: i) according to the criterion of Equations (3) and (4), as in the FL algorithm and ii) every time it finds a new putative global minimum (i.e., when its energy is lower than the best energy found by the hw so far). In case (ii), the hw takes configuration and energy from the lw with probability 1.
The FL algorithm can be used also with multiple flying walkers running in parallel, each with its individual landing walker. In this scheme, referred to as FL-PEW algorithm (with PEW standing for Parallel Excited Walkers), the flying walker may interact with each other in the space of a suitable order parameter, as in the PEW algorithm, [49,50] while the landing walkers do not interact with each other. Hiking walkers can be added too (FLH-PEW algorithm).

Results and Discussion
In this section we present some case studies to test the efficiency of the different algorithms. In all cases, we consider binary metallic systems whose interactions are described by the second-moment tight binding potential, also known as Gupta potential. [54][55][56][57] The form and parameters of the potential are given in the Supporting Information.
The AuCu, AuRh, and PtNi systems are chosen because they have been the subject of previous optimization studies so a comparison can be made. The AgCu system is chosen to deal with the optimization of large nanoalloys, containing up to a few thousand atoms, also comparing the outcomes of unseeded and seeded optimizations. AgCu is a widely studied system from the experimental point of view, [58][59][60][61][62][63][64] and, as such, it has been the subject  In this simulation, one fw is used together with its associated lw, which is restarted according to the probability of Equations (3) and (4) with = 1. When the lw is restarted its energy suddenly increases to match the energy of the fw. Then the energy sharply decreases due to the exchange and surface rearrangement moves accepted at low temperatures. of many optimization studies, [65][66][67][68][69][70][71] but systematic unseeded searches for sizes of few hundred atoms and above are still missing.

Benchmarks: Optimization of Au 300 Cu 100 , Au 200 Cu 200 , and Au 300 Rh 100
The Au 300 Cu 100 , Au 200 Cu 200 , and Au 300 Rh 100 systems were studied in Rossi et al. [51] by the BH algorithm in which different shape-changing moves and exchange moves were accepted at different temperatures. Here we compare the BH optimization algorithm as used in ref. [51] with the FL and FLH algorithms.
In agreement with ref. [51], we find that for all systems the global minimum is a Mackay icosahedron (Ih) with an incomplete outer shell as shown in Figure 5, where the best structures belonging to the other relevant motifs are shown too. These motifs are the Marks decahedron (Dh), [72] the face-centered-cubic truncated octahedron (fcc), and the icosahedron with incomplete anti-Mackay shell (Ih-aM). [67] Here below we analyze the performance of the BH, FL, and FLH algorithms separately for each system.  In this simulation, one fw is used together with its associated lw and hw, which are restarted according to the rules explained in the text.

Au 300 Cu 100
For Au 300 Cu 100 we performed searches of N = 300 000 steps for each algorithm. Note that this implies that in the FL algorithm both the fw and the lw make 150 000 steps, whereas in the FLH algorithm, each walker makes 100 000 steps. In all simulations, the same types of moves were used, namely the Brownian and the bonds moves, which mainly affect the cluster shape, and the random exchange move, which mainly affects chemical ordering. [50] The Brownian move consists of a molecular dynamics run of 200 steps at a temperature of 2000 K.
For each algorithm, two settings of the moves were used, and, for each setting, 100 independent searches were run to accumulate 200 simulations per algorithm. The settings of the simulations are reported in Tables S5 and S6, Supporting Information, for AuCu and AuRh, respectively.
All algorithms found lowest-energy structures with exactly the same geometric shape but slightly differing in chemical ordering, with the best and second best chemical ordering found by the FLH and FL algorithms, respectively, as can be seen in Table 1.
Here we display the difference in energy between the global minimum found for any of the three system and the best structure found for the other relevant structural motifs, for all three algorithms. For example, ΔE best Dh in first row is the difference between the best decahedron found by the standard BH algorithm and the global minimum, which is a Mackay icosahedron found by the FLH algorithm. The main difference between the performances of the three algorithms is shown in Figure 6, and concerns the frequency at which the searches were able to converge close to the lowest-energy structure. In 300 000 steps, FLH was able to get within 0.25 eV from the lowest-energy structure 27 times,   Table S5, Supporting Information) whose local minimization is much faster than that of shape-changing moves. At the same CPU time, the number of searches that ended within 0.25 eV was 27, 17, and 5 for FHL, FL, and BH. On the other hand, if we compare the efficiency to explore other structural motifs beyond the dominant one, FL obtains better results than FLH and BH (see Table 1).

Au 200 Cu 200
For Au 200 Cu 200 we adopted the same scheme as for Au 300 Cu 100 . The lowest-energy Mackay icosahedral shape was found by all algorithms, but FLH was able to perform a better optimization of chemical ordering than FL and BH, thus obtaining the lowestenergy structure (see Table 1). The histograms of Figure 7 show that, in 300 000 steps per simulation, FLH was able to optimize within 0.25 eV from the best structure in 11 simulations over 200, compared to the nine of FL and to the two of BH. These numbers become 11, 4, and 0 when the comparison is made at the same CPU time.
As regards other structural motifs FL performed slightly better than FLH, whereas BH was less efficient than both.

Au 300 Rh 100
For Au 300 Rh 100 we performed 100 searches of 300 000 steps for each algorithm. Also in these simulations, the same types of moves were used, namely the Brownian and the bonds moves, which mainly affect the cluster shape, and the tailored exchange move, which mainly affects chemical ordering and favors the convergence toward phase-separated configurations. [53] The Brownian consisted of a molecular dynamics run of 200 steps at a temperature of 3000 K.
The settings of the searches are reported in, Table S6, Supporting Information. For these systems, whose most favorable chemical ordering is core@shell Rh@Au, [51] all algorithms were able to find the lowest-energy structure. However, FLH was able to find low-energy structures more frequently than FL and BH, as shown in Figure 8. In fact, when comparing the algorithms for the same number of steps, FLH searches found isomers within 0.25 eV from the global minimum in 28 simulations over 100, against 17 and 9 of FH and BH. The superior performance of FLH was even more evident in the comparison for the same CPU time (28 against 14 and 7 searches within 0.25 eV). On the other hand, FL was more efficient in exploring other motifs than the Mackay Ih (see Table 1).

Benchmark: Optimization of Pt 19 Ni 19 and Pt 27 Ni 28
For PtNi, the parameters of the Gupta potential were taken from Yang et al., [31] who optimized all compositions for sizes 38 and Adv. Theory Simul. 2023, 6, 2300268 www.advancedsciencenews.com www.advtheorysimul.com  55 with a new improved algorithm. Here, we concentrate on Pt 19 Ni 19 and Pt 27 Ni 28 , for which complete data on the CPU times are available. [31] For both Pt 19 Ni 19 and Pt 27 Ni 28 we run 100 independent searches by the FLH algorithm. The settings of the moves for these searches are given in Table S7, Supporting Information. We note that at variance of the previous cases, shake moves instead of Brownian moves were used for shape changes. In the shake moves every atom was displaced randomly around its present position with a sphere with radius of 1.3 Å.
Each run was of 100 000 steps, and took on average 234 s for Pt 19 Ni 19 and 312 s for Pt 27 Ni 28 , running on a single core of a workstation equipped with a quad-core Intel(R) Core(TM) i7-6700K CPU with 16 GB of RAM. The FLH algorithm was able to locate the same global minima as in ref. [31], with a success rate of 85/100 for Pt 19 Ni 19 , which is a very high rate compared to the 12/40 of ref. [31]. For Pt 27 Ni 28 the success rate was 100/100. The average CPU times of successful searches were very short: 31 s for Pt 19 Ni 19 and 39 s for Pt 27 Ni 28 , corresponding to 13 300 and 12 500 steps on average, respectively. These times are considerably shorter than those in the previous literature. The complete data about the time distribution of successful searches are given in Figure 9. Finally, we note that the relatively easier optimization of Pt 27 Ni 28 compared to Pt 19 Ni 19 is due to the single-funnel character of the PES of the former and the double-funnel character of the PES of the latter, [73] where the global minimum is in competition with a regular truncated octahedral structure whose energy is higher by 0.10 eV.

Optimization of Larger Nanoparticles: The Case of AgCu
In this section, we tackle the problem of optimizing larger binary nanoparticles. We choose the AgCu system because it has been the subject of several optimization studies from very small sizes to a few hundred atoms, in which both shape and www.advancedsciencenews.com www.advtheorysimul.com  chemical ordering were fully optimized starting from random configurations. [65,[67][68][69][74][75][76][77][78][79] Optimization of chemical ordering alone at a fixed shape has been performed for AgCu nanoparticles of a few thousand atoms. [70,71] Here we extend the size of full global optimization well above 10 3 atoms, with the aim of treating the following issues: i) In the size range up to a few hundred atoms, the best structures of AgCu nanoparticles can be significantly different from those of the pure clusters Ag and Cu of the same size. [65,67,68,79] Is this fact still true for larger nanoparticles of a few thousand atoms? If yes, this would mean that the bulklike behavior is attained at a much larger size in nanoalloys than in the corresponding pure systems. ii) The mixing energy in AgCu nanoalloys of a few hundred atoms is negative for almost all compositions (see the cases of 100 and 200 atoms in ref. [79]). This is at variance with bulk AgCu, where mixing energy is always positive. [80] What is the sign of the mixing energy for nanoparticles of several thousand atoms? iii) Is the optimization of chemical ordering alone at a fixed shape able to produce low energy structures resembling those of full global optimization?
In order to treat these issues, we optimize AgCu nanoparticles of sizes N = 586, 1289, 2406, and 4033 atoms, corresponding to the magic sizes of regular truncated octahedra. [81] For N = 586, 1289, and 2406 we fully optimized the nanoparticles starting from random configurations and consider compositions with N Ag = 0.1 N, 0.2 N, … ., 0.9 N, whereas for N = 4033 we only perform seeded optimizations for two compositions only. We calculate the mixing energy E mix as a function of the composition defined as (5) www.advancedsciencenews.com www.advtheorysimul.com where N = N Ag + N Cu and E(N Ag , N Cu ) is the energy of the global minimum. E(N Ag , 0) and E(0, N Cu ) are the energies of the reference structures of pure clusters, in our cases fcc regular truncated octahedra. [1,81] Representative lowest-energy structures for the different sizes are shown in Figures 10-12.
For N = 586, the substitution of 10% Ag in Cu or of 10% Cu in Ag is already sufficient to cause a transition from the fcc truncated octahedron to the icosahedron, whose magic number is at N = 561, so that the 25 additional atoms are placed in an island with anti-Mackay stacking. [67] Truncated octahedral structures however reappear for intermediate compositions, [76] but they are different from the regular truncated octahedron of the pure Ag and Cu clusters. An example of these truncated octahedral structures is given in Figure 10 for composition Ag 234 Cu 352 . This structure is core@shell Cu@Ag. It presents an Ag surface and only one Ag internal atom. In addition, there is an Ag adatom on top of its surface. Therefore this structure is magic, with a monolayer Ag shell, for N = 585 and composition Ag 232 Cu 353 . It is asymmetric, with two 3 × 3, two 4 × 3, one 5 × 4, and one 4 × 4 square facets. For comparison, the regular truncated octahedron for N = 586 has six symmetric 4 × 4 facets.
For N = 1289 and 2406 (see Figure 11), icosahedral structures are still dominant for a wide spectrum of compositions, especially in the intermediate range. Truncated octahedral or decahedral structures were found for Ag-rich and Cu-rich compositions.
Finally, for N = 4033 calculations are very cumbersome. Therefore we consider only two compositions, Ag 2823 Cu 1210 and Ag 3000 Cu 1033 , and in both cases, we run seeded optimizations. As seeds, we used regular truncated octahedra or incomplete icosahedra with random chemical ordering. The optimization runs allowed for both shape-changing and exchange moves. Also for this size (results in Figure 12), icosahedral nanoparticles were lower in energy than truncated octahedral ones.
In relation to point (i), our result indeed shows that icosahedral structures are still dominant for intermediate compositions in a size range where for elemental nanoparticles truncated octahe-dra are much lower in energy. [82] This effect is due to the atomic stress release effect allowed by size mismatch between Cu and Ag atoms, [83] since Cu atoms are about 13% smaller. The transition from truncated octahedra to icosahedra experimentally observed in AuPd clusters has been attributed to stress release allowed by size mismatch. [84] These results demonstrate that the bulk-like behavior in nanoalloys can be reached at much larger sizes than in elemental nanoparticles.
Regarding point (ii), we calculated E mix as a function of composition for N = 586, 1289, and 2406. The results are shown in Figure 13. In all cases, E mix shows a rather regular behavior as a function of N Ag ∕N. For N = 586 and 1289, E mix is always negative at least in the range of 0.1 ≤ N Ag ∕N ≤ 0.9. For N = 2406, there is an Ag-rich range, roughly above N Ag ∕N > 0.8 in which E mix is positive. For N = 4033, the mixing energies of the two lowest-energy icosahedral nanoparticles, whose compositions are in the range 0.70 ≤ N Ag ∕N ≤ 0.75, are slightly negative. From the results of Figure 13, we note also that E mix has a rather sharp minimum whose position shifts to smaller and smaller N Ag ∕N as size increases.
The change of sign of E mix with composition should indeed be expected for any size N of the reference truncated octahedron. In fact, in the extreme Cu-rich limit, E mix will always be negative because, in a pure Cu truncated octahedron, the substitution of one Cu vertex atom by an Ag atom produces a negative mixing energy for all sizes of the truncated octahedron. On the contrary, in the extreme Ag-rich limit, the substitution of one Ag atom with a Cu atom, whose best position is subsurface [85] also according to DFT calculations, [59] produces a positive mixing energy. Therefore, since the two extremes have opposite signs, there must be at least one change of sign in between. With increasing size, the Cu-rich composition range in which E mix is negative shrinks. We expect that the position of the minimum of E mix shifts approximately as N surf ∕N, where N surf is the number of surface atoms because we find the lowest values of E mix in correspondence to the structures with a Cu core surrounded by a monolayer-thick Ag shell. N surf can be estimated as the number of surface atoms in an elemental truncated octahedron rescaled by (r Cu ∕r Ag ) 2 to account for size mismatch. This estimate gives N surf ∕N = 0.364, 0.294, and 0.245 for N = 586, 1289, and 2406, which correspond well to the minima of E mix in Figure 13.
To deal with point (iii) we considered the optimization of N = 586 by three different approaches. The first approach is full global optimization of shape and chemical ordering from random initial configurations, whose results have been already discussed before. The second approach is a seeded optimization which starts from a truncated octahedron with random chemical ordering and used only exchange moves. In this second approach, the cluster shape can change only because of local minimization. In the third approach, the seeded optimization uses mainly exchange moves but with a relatively small percentage (10-20%) of shapechanging moves. The results of the three approaches are compared in Figure 14. In the seeded searches with exchange moves only, E mix is negative below N ≃ 400, in qualitative agreement with the results of the symmetry-constrained optimizations of refs. [70,71], but with somewhat lower values, indicating that the absence of symmetry constraints may allow a somewhat better optimization. However, this is not the main issue, because this type of optimization itself is often not satisfactory at all. In fact,   the inclusion of a relatively small percentage of shape-changing moves in the seeded searches can lead to major improvements, especially in the intermediate composition range, producing fcc structures with much more regular shape and better chemical ordering, whose energy is often very close to that obtained by full global optimization. For example, structures (a) and (b) of Figure 14 at composition Ag 234 Cu 352 were obtained by exchanges only and by exchanges and Brownian moves together, respectively. Structure (a) presents a quite rough surface exposing some Cu patches. Structure (b) has a smooth surface containing Ag atoms only. As a result, structure (b) is lower in energy than (a) by almost 13 eV. For Cu-rich and Ag-rich compositions, the results of chemical ordering optimization are closer to those of full optimization. [86,87] Therefore the answer to point (iii) is no in many cases. The optimization of chemical ordering alone at fixed geometrical structures is often unable to lead to low-energy structures in AgCu because the search often remains trapped in unphysical highenergy shapes. Even if the optimization of chemical ordering alone starts with the correct motif, the final results may turn out to be quite unsatisfactory. The equilibrium statistical weight of structures such as that of Figure 14a is completely negligible, and its shape is quite far from any physically sensible structure which may be produced by the kinetics of a true growth process. We believe that these findings can be generally extended to systems in which atomic species present relevant size mismatch, for example AgNi, AgCo, AuCo, AuCu, AuRh, AuPt, and AuNi.

Conclusions
In this work, we have proposed new algorithms (FL and FLH algorithms) for the full global optimization of nanoalloys. In these algorithms, there are different BH walkers with specific rolesthe flying walker, the landing walker, and (in FLH only) the hiking walker. These walkers act together to ensure a wide exploration of different shapes together with a deep refinement of chemical ordering.
The effectiveness of FL and FLH has been tested in the optimization of a few systems, for sizes from a few ten to a few thousand atoms. The tests have revealed that both FL and FLH are more efficient than other approaches available in the literature. FLH shows some advantages over FL for what concerns the exploitation task, since it is able to better optimize chemical ordering within the global minimum geometric structure. On the other hand, FL shows a somewhat better performance in exploration of different funnels than that of the global minimum. These differences are due to the different balance between exploration and exploitation moves in the two algorithms. We believe that an effective strategy when dealing with a new system is to run searches by both algorithms.
Further improvement in efficiency may be expected by using multiple flying walkers running in parallel, each with its individual landing (and eventually hiking) walker. In systems with tendency to intermixing, such as AuCu and PtNi, a further improvement in efficiency may be expected by implementing more efficient schemes for homotop search [71] instead of random exchanges.
Our results have also shown that the optimization of shape and chemical ordering together is indeed necessary in many cases, since the optimization of chemical ordering alone at a fixed shape (apart from local relaxations) may lead to unphysical results. This is especially true for intermediate compositions in systems with size mismatch between the atomic species.
Finally, the study of AgCu nanoalloys for sizes up to 4000 atoms has shown that the bulk-like regime, in which the optimal structures are fragments of the fcc lattice common to both metals, is not yet attained for sizes of a few thousand atoms, since icosahedral structures are found as the lowest in energy for several compositions. This behavior is quite different from that of elemental Ag and Cu nanoparticles, in which icosahedra are already largely unfavorable for sizes above a few hundred atoms.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.