First principles global optimization of metal clusters and nanoalloys

ABSTRACT The global optimization of nanoparticles, such as pure or bimetallic metal clusters, has become a very important and sophisticated research field in modern nanoscience. The possibility of using more rigorous quantum chemical first principle methods during the global optimization has been facilitated by the development of more powerful computer hardware as well as more efficient algorithms. In this review, recent advances in first principle global optimization methods are described, with the main focus on genetic algorithms coupled with density functional theory for optimizing sub-nanometre metal clusters and nanoalloys. GRAPHICAL ABSTRACT


Introduction
In today's world, the global optimization (GO) of functions according to one or more criteria has become very popular in many different research fields and applications [1][2][3][4][5][6][7][8][9][10]. By the term 'global optimization' we mean finding the overall best solution for a given mathematically formulated problem, which usually corresponds to the global maximum or the global minimum (GM) of a function or a set of functions.
An important goal in modern nanoscience is the control of materials on the atomic scale, coupled with achieving a fundamental understanding of how physicochemical properties depend on the particle structure in order to tailor nanomaterials for real-world applications [11]. Thus, computeraided optimization can be used to predict a suitable atomic structure corresponding to an ideal (or near-ideal) value of the property of interest and to interpret experimental results, as well as guiding the development of future materials in the course of rational material design. Therefore, GO of nanosystems plays a very important role.
Metal clusters have received enormous interest due to their industrial applications in nanotechnology, electronical, biological, medical devices and catalysis . Furthermore, metal clusters are ideal test systems for probing physical theories. The limits of concepts for bonding theories can be evaluated particularly well with the help of high-quality experimental data. The range of properties of metallic clusters can be significantly enhanced by combining more than one metal, forming 'nanoalloys' [14]. These mixed metal clusters allow the modification of chemical and physical properties by changing, in addition to their size and shape, their composition and chemical ordering. The system complexity increases significantly for mixed clusters due to the presence of 'homotops' (inequivalent permutational isomers), which makes the GO a more challenging task. For nanoalloys, there are also several possible mixing patterns (chemical ordering). Common types of chemical ordering are: random or ordered mixed; core-shell and multi-layer (onion-like); and phase segregated (Janus) [14]. The preferred chemical ordering depends on several factors and is a balance between the cohesive energies, surface energies, relative atomic sizes and specific electronic and magnetic effects.
The possibility of performing GO on ever larger systems and applying more rigorous quantum chemical electronic structure methods has gone hand-inhand with the development of increasing computing power. For an accurate description of the smallest clusters (or sub-nanometre particles), a quantum chemical description is essential, due to the increasing importance of quantum size effects. The most rigorous methods are ab initio electronic structure methods, which only rely on the laws of nature without making use of empirically fitting parameters or additional assumptions [36]. The input to ab initio electronic structure methods are physical constants, atomic numbers, the positions of the nuclei and the number of electrons [37]. Popular ab initio methods are Hartree-Fock (HF) methods and post-HF methods such as Møller-Plesset-Perturbation theory (MPn), Configuration Interaction (CI) and Coupled Cluster (CC) methods. Another very widely used method in material science is Density Functional Theory (DFT). Though DFT is strictly not an ab initio electronic structure method, in this review we will describe both DFT and ab initio electronic structure methods by the term 'first principles', as in the review by Heiles and Johnston [10].
The main focus in this review is to introduce and present examples of GO methods, applied to sub-nanometre pure metal clusters and nanoalloys, using first principle methods.

Global optimization methods for metal clusters and nanoalloys
Within the Born Oppenheimer approximation, stable structures for a given cluster (e.g. different stable isomers) correspond to local minima on the multidimensional potential energy surface (PES). The PES of a cluster describes its energy as a function of its atomic coordinates. GO corresponds to finding the most stable atomic arrangement for a specific cluster (defined by the number and type of its constituent atoms as well as by the total charge) that is the lowest energy-point on the PES, the so-called 'global minimum'. The significant effort that has been expended in finding cluster GM [10,[38][39][40][41][42][43][44][45][46][47][48] is because the GM is usually the most likely structure to be formed in an experiment (though kinetic factors should not be forgotten). Even if the experimental structure is not the GM, typically it is a very low energy isomer [46]. In combined experimental and theoretical studies, the presence of one or more low lying isomers have been found to explain experimental results very well [49][50][51]. The approach in such combined studies is: (i) to measure a physico-chemical property of a particular cluster; (ii) to use a GO technique to search for the GM and other low energy structures and (iii) to calculate the same property for these isomers. Finally, the comparison between theory and experiment can lead to structural assignment discrimination of the cluster (see Figure 1) [51,52]. However, this approach is complicated by the fact that the number of possible stable structures (local minima on the PES) increases (at least) exponentially with the size of the cluster and, hence, so does the computational effort required. For mixed clusters, such as nanoalloys, the problem is compounded due to the combinatorial increase of possible minima (so-called 'homotops' [29]) resulting from swapping positions of different elements.
Deterministic location of the GM would be possible with an exact knowledge of the PES. However, this requires calculating the entire PES Figure 1. Combined experimental and theoretical approach for cluster structure discrimination. (a) Reproduced from Ref. [52] with permission from Wiley-VCH Verlag GmbH & Co: compares experimental electrostatic beam deflection measurements and calculated beam profiles for the globally optimized Sn 6 Bi 3 cluster and other low-lying isomers. (b) Reproduced from Ref. [51] with the permission of AIP Publishing: The experimental absorption spectra obtained by UV-Vis photodissociation spectroscopy is shown, along with the calculated absorption spectra for the lowest lying AgAu 3 + isomers.
to determine all its minima and saddle points, which is only possible for the simplest cases and is generally not feasible using relatively expensive ab initio (or even DFT) electronic structure methods, for all but the smallest systems. For these reasons, sophisticated GO methods are essential to perform an efficient and unbiased exploration of the PES. In general, a good GO strategy must combine both local and global aspects of PES searching. Whereas local searching addresses the sampling of single solutions, global searching deals with the efficient exploration of different regions of the multidimensional parameter space [53].
In constructing a GO for clusters, there are two issues to deal with. From a mathematical point of view this corresponds to the questions of how to choose good discrete argument values and how to calculate the target set. The first question pertains to the GO algorithm for generating cluster structures (defining atomic coordinates) in order to efficiently explore different regions on the PES. The second question concerns the level of theory used to calculate the energies (points on the PES) for the generated geometries.
Common classes of algorithm which are used for the GO of clusters are: Genetic algorithms (GA) [46]; basin hopping (BH) [45]; particle swarm optimization (PSO) [54]; artificial bee colony (ABC) [55]; simulated annealing [56] and threshold algorithms [57]. For a given cluster, with a fixed number of atoms and elemental composition, the computational time strongly depends on the selected level of theory. Therefore, a trade-off between accuracy and cost must be found. The choice of the appropriate level of theory for describing the chemical bonding between the atoms of the cluster depends on the system size as well as the available computational resources, which are dependent on the development of more efficient and ever faster computer hardware. Hence, in early GO studies of clusters, generally only empirical (or semi-empirical) potentials (EP) were used [58][59][60][61][62][63][64][65][66]. Even today, many studies employ this EP-GO approach, which is usually followed by geometry refinement with more rigorous first principle methods [67][68][69][70][71][72][73][74][75][76][77][78][79]. The main advantage of this approach is that it can also be applied to large clusters, of several hundreds or a few thousands of atoms. For larger nanoparticles (up to many tens or hundreds of atoms), which already show similar properties and chemical bonding as in the bulk phase (with properties typically scaling as the surface to volume ratio, i.e. as 1/R), the EP-GO approach yields reasonable results [38,58,[80][81][82][83][84][85][86][87]. Furthermore, EP-GO is not only beneficial when systems already have bulk-like behaviour, but also when larger clusters form unusual structural motifs so that it would be difficult to guess by chemical intuition alone. For instance, chiral icosahedral and decahedral nanoalloys in the size range between 100 and 200 atoms as well as pyritohedral structures and chiral decahedra for nanoalloys up to 586 atoms have been found with the EP-GO approach and subsequent first principle calculations [71,73]. Unfortunately in general, simple EPs do not reproduce subtle effects due to the intrinsic quantum chemical nature of chemical bonding, which become increasingly important with decreasing system size, i.e. in the region where 'every atom counts' [88]. With the term 'intrinsic quantum chemical nature of chemical bonding', we mean that in EPs no electronic interactions (as well as no explicit electronic correlation effects) are considered, in contrast to first principle methods [89,90]. For example, Lyons et al. examined with an EP the relative energetics of C 44 isomers [91]. They found an overall good agreement with DFT results, but they emphasized that electronic effects are not considered in the EP calculations and special stability due to resonance effects is not addressed in their study. Although in principle it is possible to build an EP for a proper description of some specific small clusters, such potentials may not be transferable to slightly larger clusters. There are several examples where first principle calculations and experimental findings disagree with EP computations regarding the cluster structure or relative energies of different isomers. For example, GO of trimetallic Ag k Au m Pt n (k + m + n = 13, 19, 33, 38) clusters result in different GMs for EP and DFT calculations [92]. GO of Si n (n = [16][17][18][19][20] clusters employing three different EPs, as well as DFT computations, provides different GMs for each theory level [64]. A combined EP and DFT study of 40-atom Pt-Au clusters shows only a partial agreement between EP and DFT results [80]. Moreover, a theoretical EP study of the structures of neutral Ag n clusters (n = 2-60) suggests a decahedral cluster growth and compared this result with experimental trapped ion electron (TIED) findings of silver cations of the same size, but the experimental data suggest an icosahedral growth [93]. Hence, in general, first principles electronic structure methods are the better choice for the accurate description of the chemical bonding (and hence for GO) of small subnanometre clusters [89]. When only dealing with the optimization of chemical ordering within the cluster using EP-based optimization, the present size limits are well above 1000 atoms for an unconstrained search and several thousand atoms for symmetry constrained searches [75,94]. Regarding merely the optimization of chemical ordering, a DFT optimization of a 309-atom AuCu cluster (an icosahedron) has been realized [74]. Another method that has been applied is the semi empirical tight-binding DFT (TBDFT) method. The TBDFT approach has empirical parameters, but can maintain close-to-DFT accuracy, for considerably reduced computational cost [95]. Hence, provided it is used within the region of parameterization, DFTB can be an efficient method for metal cluster GO, which was recently shown for TBDFT applications to silver and gold clusters [96]. A study has also been reported which combines TBDFT and DFT directly in the GO of gold clusters [97].
Usually, when performing first principles GO of clusters and nanoparticles, the GO algorithm is coupled to a quantum chemical (QC) software package. The choice of the particular electronic structure method and package depends on available computer resources and on the size and complexity of the cluster system (e.g. pure metal cluster vs. nanoalloy, presence or absence of ligands and/or supporting surface, importance of magnetism).
GO algorithms can be distinguished in a number of ways [10]. Firstly, they can be divided into biased and unbiased methods. In biased GO, prior knowledge (or guesses) about bonding in the system or preferred structural motifs are used to accelerate the exploration of the PES. In a strictly unbiased global search, however, no prior knowledge is used. Another classification concerns the energy calculation for a generated geometry. During the GO process, there are in general two possibilities for how to perform the local energy calculation: Either only the energy of the generated cluster geometry is calculated (single point calculation) or for each created cluster structure a local geometry optimization (energy minimization) takes place so that each cluster structure corresponds to a local energy minimum on the PES.
In the optimization of clusters, a cost-saving approach (which sometimes gives good results) involves first performing GO at a lower level of theory (EP or DFT with a small basis set) and looser convergence criteria. This is followed by refinement (reoptimization) of a set of low-energy and structurally distinct isomers using a more rigorous ab initio method, such as CCSD(T), or DFT with larger basis sets and/or more computationally expensive hybrid exchange-correlation (xc) functionals [10,50,51,[98][99][100][101][102][103]. Particular attention must be paid to the fact that the relative energy differences between the isomers, as well as their energetic order, can change (especially in DFT with the use of different xc functionals and basis sets) so that an originally energetically higher isomer can become the GM after refinement (see Figure 2) [51]. It is also possible that the high theory level GM is missed if the lower level PES is not a sufficiently good match to the 'true' PES in the low energy regions.
In the following sections, while several GO methods will be briefly mentioned, we will focus on unbiased GO of metal clusters and nanoalloys at the DFT level, principally using GA.

Genetic algorithms for optimizing cluster structures
The GA is a metaheuristic based on the principles of natural evolution and belongs to the class of evolutionary algorithms (EA), which are examples of the emerging area of artificial intelligence (AI). The GA can be seen as an intelligent search algorithm, which learns features of good solutions by the recognition of schemata or building blocks during the exploration of the multidimensional parameter (coordinate) space. The GA successively improves its solution during the search progress. We only give a brief description of the GA procedure here and refer the reader to the literature for more detailed explanations [46,104].
In general, GAs can be applied to any GO problem where the variables (genes) can be encoded as a string (chromosome), which represents a test solution of the problem. The overall task is the GO of the values of the variables ('alleles') in order to find the overall best solution [46]. To do this, the GA uses evolutionary processes such as mutation, mating (also known as crossover) and natural selection. The GA starts with an initial population of a number of individuals (with the population size generally being inversely related to the computational cost of the energy minimization step), which are usually generated randomly [10]. During the GA process, the population evolves for a specified number of generations or until energy or population convergence is reached, by applying the aforementioned evolutionary operators (crossover and mutation) to each generation. The crossover method mimics biological genetic crossover, combining genetic information from two parent strings (in a GA, more than two parents can be used), in order to produce one (or more) offspring. The most common crossover operators are one-point and twopoint crossover. In one-point crossover, the parent strings are cut at the same position and offspring are generated by adding complementary Figure 2. Lowest energy isomers for Ag 2 Au 2 + for different xc functionals for the same def2-TZVPP basis set. Different xc functionals lead to different relative energies and may also cause a different energy order of the isomers. Reproduced from Ref. [51], with the permission of AIP Publishing.
(parent) genes. In two-point crossover, the parent strings are cut at two different positions. To ensure a sufficient diversification of genetic information for each generation, a mutation operator is used which supplies new genetic material. This is to prevent stagnation of a population, corresponding to convergence on a non-optimal solution.
In each generation, the fitness (objective function in the GO problem) of every individual is evaluated and serves as a measure of the quality of the trial solution. The more fit individual members of the current population are chosen for subsequent crossover or mutation in order to breed a new generation. For this purpose, popular selection methods are roulette wheel and tournament selection. This selection process imitates the natural selection concept in biological evolution, which is also known as 'survival of the fittest'. Consequently, it is an elitist strategy, because the best (highest fitness) individual is guaranteed to survive into the subsequent generation, thereby enhancing the probability of passing on some of its 'good' genetic material. Thus, a gradual improvement of the trial solutions will be achieved by applying GAs to GO problems. There are additional strategies which have been introduced to improve the performance of the GA. One strategy which can help to accelerate the GM search is biasing the GA by introducing a seeded initial population [105]. Another strategy is to ensure a certain degree of diversity in each population in order to avoid stagnation (that is to avoid being trapped in a sub-optimal local minima). This can be accomplished with replacement or deletion (also known as predator) strategies. These strategies define which member(s) of the population will be replaced with a new one (in addition to the previously mentioned fitness-based replacement) or deleted. With these methods it is possible to remove identical (or at least too similar) members in one population in order to increase the diversity. Those strategies require one or several dissimilarity measures via the definition of descriptor operators to compare each member of the population. The results are usually a set of numbers or arrays of dissimilarity measures indicating how the individuals differ from each other. Thereafter a 'predator operator' removes certain members of the population, those which are too similar or identical according to the chosen descriptors [105][106][107][108][109][110][111][112]. In the case of clusters, commonly used descriptors for dissimilarity measures are: energies, moments of inertia, radial distribution functions, connection tables, variance of atomic distances, etc. [105,106,[109][110][111][112].
The first applications of GAs for finding the GM structure of clusters dates back to 1993 for atomic clusters (Si 4 ) by Hartke [41] and molecular clusters (benzene dimer, trimer and tetramer) by Xiao and Williams [113]. In these first publications, the cluster coordinates were encoded as binary strings and thus represented as a chromosome to which the evolutionary operators (crossover and mutation) were applied. In 1995, Zeiri introduced a GA which directly used real number cartesian coordinates of the cluster atoms and hence eliminated the representation problem (encoding and decoding the cluster coordinates to binary strings) [43]. The next milestone in GA development was due to Deaven and Ho who introduced a physical cut-and-splice method for the crossover operator, as depicted in Figure 3 [42]. This operation is highly effective for the GO of clusters since it considers the spatial distribution of the atoms and hence is sensitive to structural properties (this means it is able to pass 'good' structural features from the chosen parents to the new child cluster generation) as well as being highly effective. Thus, the cut and splice operator is widely used and implemented in many GAs for clusters and nanoparticles. The cut and splice operation (see Figure 3) applies random rotations of both parents, uses a horizontal cutting plane, cuts each parent and combines complementary parts of each of one of the parents to create descent(s). More recent studies have developed further crossover operators for clusters [114][115][116]. One new form is the generalized cut and splice (GenC&S) operator [114,115]. In contrast to the Deaven-Ho approach, the GenC&S operator does not make use of a plane or random rotations of the parents before mixing the genetic material. It uses the Euclidean distance as the criterion for selecting atoms of each cluster such that subsets of atoms that are close together in the parent clusters will form the new building blocks to create the offspring. In the GenC&S version there is no bias associated in the operator as in the Deaven-Ho approach (that is using a horizontal cutting plane and forcing the new sub-cluster to be parallel to this cutting plane) [115]. Chen et al. pointed out that the cut and splice operator of Deaven and Ho may fail in the case where the cutting plane lies near the major axis of one parent and near the minor axis of the other parent [116]. In this case, the resulting structure is unstable and after local minimization the final structure shows no apparent relationship to the parents. Hence, good (low energy) structural features of the parents would not be passed to the offspring. For this reason, Chen et al. proposed a new so-called 'spherecut-splice' crossover which employs a sphere instead of a cutting plane and generates the offspring by using atoms of one parent, which lie inside the sphere and atoms of the other parent, which lie outside the sphere. This method is able to preserve good schemata and has been shown revealed to be very suitable for the GO of larger clusters.
Another innovation of Deaven and Ho was the introduction of local energy optimizations for each generated cluster structure, so that each generated cluster is relaxed to the nearest local minimum. This approach was initially applied to fullerenes (C 20 to C 60 ) [42]. GAs which couple local minimizations to the global search are called 'Lamarckian' because, from a biological point of view, this process corresponds to Lamarckian rather than Darwinian evolution [46]. Doyle and Wales demonstrated that local minimizations effectively transform the PES into a stepped surface in which each step corresponds to a basin for the respective local minimum [62]. The Lamarckian approach, thereby, improves the efficiency of the GA tremendously.
The GA-DFT approach A widely used first principles approach for GO of clusters at the DFT level is the GA-DFT approach. The first GA-DFT approach was introduced in 1996 by Tomasulo and Ramakrishna who performed GO of AlP-clusters directly at the DFT level of theory, within the local density approximation [117]. Several years later, Alexandrova and Boldyrev firstly coupled a GA code to a quantum chemical software package and introduced the gradient embedded genetic algorithm (GEGA) [132]. Other recently developed first principles based GAs are, for example, OGOLEM [124], a versatile poolbased GA implementation for clusters and flexible molecules and a GA that uses machine-learning techniques to improve its performance [129]. Though the GA-DFT approach is the most widely used method for the GO with first principle methods, there is also a GA study which uses MP2 calculations for sodium potassium nanoalloy clusters [133]. The following briefly describes the development of the GA-DFT approach within the Birmingham Cluster Genetic Algorithm and its derivatives, developed by Johnston and collaborators.

The Birmingham Cluster Genetic Algorithm (BCGA)
The BCGA is a Lamarckian-GA, which is capable of performing GO of bare pure or mixed clusters at either the EP or DFT level of theory [46,98]. In this algorithm, clusters are represented as real-valued cartesian coordinates (as introduced by Zeiri) on which the genetic operators act directly. In general, such GAs are called 'phenotype algorithms' [111]. The Deaven-Ho cut-and-splice crossover method and several mutation methods (atom displacement, twisting, cluster replacement and (inequivalent) atom permutation) are employed. In order to perform GO using plane wave (PW) DFT, the BCGA was initially coupled to the software package Quantum Espresso. The first study with this BCGA-DFT approach was the investigation of the dopant-induced 2D-3D transition of 8-atom Au-Ag nanoalloys [98].
A major disadvantage of the traditional GA approach is that the local energy minimization steps are performed sequentially, which results in a bottleneck in the GO procedure and limits the cluster sizes which can be reliably optimized. This is especially important when using first principles methods in the local minimizations. For example, employing DFT calculations for GO, typically more than 99.9% of the computational time is spent in the local minimization steps [124,130]. Therefore, parallelization of the GA, taking advantage of larger computational resources, became absolutely necessary.
The Pool-BCGA code was developed [39] to allow several independent local optimizations to be performed at the same time during the GO. In the pool approach, multiple instances of the GA act on a global database (pool) in order to perform simultaneous local optimizations of cluster structures which are generated by crossover and mutation. These GA subprocesses share the GO workload (as depicted in Figure 4), which leads to parallelization of the algorithm. During GO, the pool database consists of the geometric and energetic information of a fixed number of isomers. Each individual GA instance uses the pool by applying crossover and mutation operators to the database members in order to generate new structure candidates. Each newly formed isomer competes with the current members of the pool and replaces the highest energy database structure if the new individual has a lower energy (higher fitness). At the beginning of the GO, an initial pool is created by random generation and local minimization of a pre-determined number of structures. After the database is filled, the pool-clusters are randomly selected according to their fitness using the roulette wheel method and are subjected to the crossover and mutation operators. The pool-GA was benchmarked and applied successfully to the GO of Au 10 and Au 20 clusters. An evolution plot of the GO of Au 10 is shown in Figure 5. In addition, the influence of the mutation rate on the evolutionary process was investigated for Au 10 Pd 10 and it was found that the higher the mutation rate the greater the number of higher-energy-isomers that were found and the more structures had to be generated to find the GM. This example shows that too high mutation rates reduce the efficiency of finding the GM.
The Fortran-based pool-BPGA was subsequently incorporated within a more efficient modern Python framework and was coupled to the Vienna Ab initio Simulation Package (VASP), which is a plane wave DFT code with PAW pseudo-potentials [134][135][136][137]. The resulting Birmingham Parallel Genetic Algorithm (BPGA) was first applied to the GO of iridium clusters with 10 to 20 atoms. The methodology of the BPGA is shown in Figure 6 [40]. The BPGA can also be used for the optimization of cluster geometries in the presence of a surface and applied to surface-grown or soft-deposited clusters in order to study catalytic active supported clusters [138,139].
The latest published GA-DFT code based on the BPGA is the Mexican Enhanced Genetic Algorithm (MEGA), written by Vargas, Beltran and coworkers [130]. The code is also written in Python and coupled to the VASP code. MEGA is more efficient and flexible and a number of new features have been introduced, such as progress documentation (a text file, which documents the applied evolutionary operator to each generated cluster structure) and several new mutation implementations, such as 'move' or 'twist'. The 'move' mutation performs a random displacement of 25% of the atoms, while 'twist' mutation rotates one half of the cluster atom by a random angle with respect to the other half. This algorithm can start from scratch with randomly generated structures for performing an unbiased GO, but it can also operate in biased mode, starting from predefined pool structures ('seeded pool'). One of MEGA's main advantages is the use of two criteria (energetic and structural comparison) for testing isomer diversity, deleting similar structures and adding new structures as required. Hence, the structural diversity of the pool is maintained during the whole GO, which decreases the probability of stagnation and  increases the GO efficiency. MEGA was firstly applied for the GO of neutral and negatively charged gas-phase Au n clusters (26 < n < 30) [130].
Current developments of the BPGA have the goal of making the GA more versatile, to allow GO of more complex systems. This includes: interfacing it with different QM packages (allowing diverse electronic structure methods); the explicit consideration and optimization of cluster spin; and the optimization of clusters in the presence of different chemical species (ligands or reactants). Further developments will include the simultaneous fitting of chemical and physical properties during the GO process.
Other first principles approaches for optimizing cluster structures In addition to the GA approach for the GO of clusters, several other effective methods exist, which cannot be discussed in detail here. A brief introduction is given below for a number of first principle-based GO techniques which have been applied for cluster optimization.
Particle Swarm Optimization (PSO) is another nature-inspired algorithm, which is based on social-psychological behaviour of birds flocking and belongs to the category of swarm intelligence methods [140]. PSO is a population-based algorithm, in which the behaviour of every individual (agent) is determined by swarm intelligence to probe promising regions of the PES [54,141]. The GO of Li clusters was successfully carried out with PSO in combination with DFT within the CALYPSO methodology [142].
One of the most popular GO methods is the Basin Hopping (BH) algorithm, which is a Monte Carlo (MC) method that combines random hopping moves with local minimization. Similar to a 'Lamarckian' GA, the PES is converted to a set of basins of attraction of all local minima [62,63]. One of the first examples of a direct DFT-based BH GO of metal cluster was carried out by Aprà et al. [143]. Usually the BH algorithm leads to random jumps between nearby minima, which corresponds to an exploration of vicinity basins on the PES around the starting geometry. A disadvantage of the traditional BH method is that there is a risk of being trapped in a very low local minimum ('funnel') if it is surrounded by high energy barriers. To overcome this limitation, several methods have been developed such as 'jumping' [144]. Other developments which uses descriptors to tackle this trapping problem are the Parallel Excitable Walkers (PEW) algorithm [145] and the population-based BH [112,146]. The use of first principles methods such as DFT for the local minimization step with the BH algorithm has become a well-established method in the GO of clusters [147][148][149][150][151][152][153].
Another search scheme for GO is the threshold algorithm, which uses MC methods to explore the PES. The threshold algorithm performs a stochastic investigation of the PES with the restriction to keep the energy below some threshold (lid) values [154,155]. The procedure can be considered as constrained random walk, which allows an ergodic sampling of the available PES. The threshold algorithm is also suitable to investigate energy barriers between local minima. New regions of the landscape become available by increasing the threshold and repeating the stochastic exploration of the PES. For example, DFT has been used in a combined threshold algorithm and GA study to investigate the structures and interconversions of low-lying isomers of Cu 4 Ag 4 [156].
The minima-hopping-algorithm (MHA) uses molecular dynamics (MD) steps to explore the configurational space of the cluster system. The MHA of Goedecker makes use of MD steps with subsequent local minimization and a historical list of the already visited minima to restrict the algorithm to search new regions of the PES [157]. MHA has been successfully used at the DFT level to investigate structures of Mg clusters [158].
Probably the first published GO method for cluster structures was simulated annealing (SA), but standard SA [56] is significantly outperformed by GA and BH methods [41][42][43]. The SA method consists of two steps: first the initial system is heated up to a high temperature, with the system evolving using MC or MD simulations; subsequently, the temperature is slowly decreased in a stepwise fashion. In order to perform first principles MD simulations, SA was used in combination with DFT-based Car-Parrinello simulations [159] and applied to predict the geometry of Se (Se 3 to Se 8 ) and sulphur clusters (S 2 to S 13 ) [160,161]. SA has also been coupled with HF calculations, with subsequent random quenching, to search for the GM of LiF clusters [162].

Example applications of first principles GO methods to metal clusters and nanoalloys
In the following section we describe several examples of GO studies of metal clusters and nanoalloys employing first principles methods. For a more complete and detailed summary, the reader is referred to the following references [10,47,53,163].

Au clusters
Clusters of noble metal elements are of considerable research interest, due to their intriguing optical properties and promising applications in catalysis and other technologies. Gold shows unique physical and chemical properties, which are significantly influenced by relativistic effects [11,29,88,90,92,96,103,107,109,111]. There have been several GO studies on gold clusters using EPs followed by DFT relaxation [164][165][166]. Assadollahzadeh and Schwerdtfeger performed a GO of neutral Au n (n = 2-20) clusters using a seeded GA-DFT approach. They predicted planar structures as the putative GM up to Au 10 and a tetrahedral structure (T d ) for Au 20 [167]. The experimentally confirmed (by farinfrared-multiple photon dissociation (FIR-MPD) spectroscopy) [168] highly symmetric tetrahedral structure (T d ) of Au 20 was also obtained with a BH-DFT approach [143] and with the pool-BPGA at the DFT level by Shayeghi et al., who also found the same planar geometry for the GM of Au 10 [39]. The BH-DFT(GGA-PBE) method was applied by Jiang and Walter for Au 40 , for which a twisted pyramid (C 1 symmetry), with a tetrahedral core, was obtained as the putative GM [116]. New GM structures, with core-shell structures, were proposed for Au n (n = 27-30) using MEGA-DFT(GGA-PBE) [130].
For anionic Au n -(n = 36, 37, 38) clusters a BH-DFT approach was used by Zeng et al. in a joint experimental (photoelectron spectroscopy) and theoretical study [169]. They found for the most stable geometries coreshell structures with a tetrahedral core. A combined experimental (trapped-ion electron diffraction: TIED) and theoretical study, using DFT MD (TPSS meta-GGA) with a quenching algorithm, of Au 12 found quasiisoenergetic 2D and 3D structures within the limits of the employed experimental and theoretical methods [170].
The structure of small Au 4 + clusters was elucidated by a combined experimental (photodissociation spectroscopy) and theoretical approach, using the BCGA-DFT method [103]. A rhombic D 2h structure was obtained in accordance with the previous study of Gilb et al. [171]. Moreover, the experimental depletion data indicates also an additional contribution of a Y-shaped isomer with C 2v symmetry, which was also experimentally observed by Dopfer et al. in a photodissociation experiment [172].

Ag clusters
Several joint experimental and theoretical studies have proven to be very powerful tools in the structure determination of small silver cluster cations [49,102,103]. For Ag 10 + , Ag 8 + , Ag 6 + and Ag 4 + , using the GA-DFT approach, it was possible to explain experimental UV-Vis photodissociation spectra in terms of one or more lowest lying energy isomers [49,102,103]. The GM structure of Ag 4 + is a rhombus with D 2h symmetry, while the putative GM geometry of Ag 6 + is a bicapped tetrahedron (C 2v ) and the lowest lying isomer of Ag 8 + was found to be a pentagonal bipyramid with a single face capped (C s ). The putative Ag 10 + GM structure is built up of two interpenetrating pentagonal bipyramids (D 2h ).
The GO of neutral Ag n (n = 5-12) clusters was carried out using the GA-DFT approach and compared to neutral gold and mixed AuAg clusters of the same size [173]. For neutral silver clusters, the 2D-3D transition was determined to lie between the hexamer and heptamer clusters, i.e. significantly smaller than for gold.

Pt clusters
Pure platinum clusters (trimers to hexamers) as well as Ti-and V-doped Pt clusters have been investigated by Jennings and Johnston with the GA-DFT approach [174]. It was found that varying the spin multiplicity has a strong effect on pure Pt clusters regarding energies and structures, whereas different effects are obtained for the doped ones, where partial spin quenching was observed. In general, higher spin states of PtTi clusters lead to higher energies. For singly V-doped Pt clusters the GM structures are all doublets whereas the doubly doped clusters have triplet ground states.

Ir clusters
Direct GO with spin polarized DFT was performed with the BPGA for Ir n (n = 10-20) clusters [40]. In agreement with previous CCSD(T) and CASSCF calculations [175], simple cubic structures (or derivatives) were found as the GM.

Sn clusters
There have been several first principle GO studies of main group metal clusters. For example, neutral Sn n (n = 6-20) clusters were examined with the DFT-GA approach in combination with electric beam deflection experiments [176]. The GM for n = 18-20 were found to be stacked prolate structures, with at least one trigonal prism motif. Singlet states were found to be most stable. Ahlrichs et al. examined the structures of tin cluster cations Sn 3 + to Sn 15 + , combining ion mobility measurements and unbiased GA-DFT GO [177]. All the small Sn n + clusters up to the heptamer exhibit a Jahn-Teller distortion. The GM structures for the larger clusters are mainly single or multiple capped prisms or antiprisms. Schooss et al. employed a joint GA-DFT approach with TIED and collision-induced dissociation to determine the structures of medium sized tin cluster anions Sn n − (n = 16-19) [178]. They found prolate cluster geometries, which are composed of at least one of three frequently occurring building blocks: pentagonal bipyramid; tricapped trigonal prism; bicapped tetragonal antiprism.

Pb clusters
Lead is a heavier homologue of tin. Its properties are strongly influenced by relativistic spin-orbit (SO) effects. Götz et al. elucidated the structure of neutral Pb n (n = 7-18), combining BCGA-DFT and molecular beam deflection studies [101]. The GM cluster geometries were mainly spherical and showed a strongly deviating growth pattern compared to tin clusters: whereas tin cluster adopts mainly prolate geometries, lead prefers a denser packing and forms more spherical structures. In a later study, Götz et al. investigated the structure of larger neutral Pb n (n = 19-25, 31, 36, 54) clusters with BPGA-DFT, finding that Pb n clusters with up to 36 atoms are not metallic [50]. They found oblate structures for the GM of Pb 19 and Pb 20 , while prolate structures were obtained for Pb 21 to Pb 25 .

Al clusters
A GM search for Al n clusters (n = 2-30) was performed with a BH-MC-DFT algorithm [126]. The 2D-3D transition was found to take place between n = 4 and n = 5, with larger Al clusters exhibiting closely packed structures. Al 14 to Al 19 display structures based on an underlying icosahedral Al 13 motif.

Nanoalloys
Mixing two or more metals in order to form nanoalloys leads to a tremendous increase in the tuning range of chemical and physical properties by adding two variables: composition and chemical ordering (the degree of mixing or segregation) [14]. However, as mentioned previously, this makes the GO procedure more complex.

Sn-Bi clusters
In a joint theoretical and experimental (electric field deflection method) approach, Sn m-n Bi n clusters (m = 5-13, n = 1-2) were studied using the BCGA-DFT approach [100]. GO calculations used plane wave DFT. Low energy structures were reminimized using orbital-based DFT and electric dipole moments and polarizabilities were calculated for comparison with the experiments. In nearly all cases, the lowest lying, metastable isomers were found to be homotops of the GM. The low energy structures are triangular faced polyhedral ('deltahedra'), which are also found for borane clusters and isoelectronic Zintl anions [179].
Au-Ag, Au-Cu and Ag-Cu clusters Nanoalloys of the coinage metals (Cu, Ag and Au) have been studied widely [14]. Heiles et al. employed BCGA-DFT to study the composition-dependence of the 2D-3D transition in neutral octameric gold-silver clusters (Au 8n Ag n , n = 0-8) [98]. They predicted the 2D-3D transition to lie between Au 6 Ag 2 and Au 5 Ag 3 . Subsequently, Heard et al. also used BCGA-DFT to find the GM for neutral and charged Au 8-n Ag n (with n = 0-8) clusters [180]. They found that the 2D-3D transition is highly sensitive to both charge and composition, with 3D structures becoming favoured at Au 8 + and Au 2 Ag 6 − for cationic and anionic clusters, respectively, which is consistent with previous studies that have shown that, for pure gold clusters, the 2D-3D transition occurs at smaller sizes for cations and larger sizes for anions, compared to neutral clusters [98,180]. The structure of mixed gold silver tetramers was elucidated by a combined experimental and theoretical approach. Using photodissociation spectroscopy and BCGA-DFT [51], Zhao et al. performed an unbiased GM search of all possible compositions of Au n Ag m in the size range 4 < n + m < 13 using GA-DFT(GGA-PBE) [173]. In agreement with previous studies, they found the same (octameric) GM structures and ascertained that planar binary Au-Ag clusters retain the geometry of the pure gold clusters, while 3D structures show a similar shape to the pure silver clusters.
The BCGA-DFT approach was used for the GM search of octameric AuCu and AgCu clusters [181]. All clusters except Au 8 and AuCu 7 were found to have compact 3D structures. Thus, the 2D-3D transition takes place after AuCu 7 .

Au-Ir clusters
Due to the promise of Ir-doped Au nanoparticles in catalysis (CO oxidation), the GM search of AuIr nanoparticles is of particular importance. Adding Ir to Au catalysts leads to an improvement of the catalytic activity and prevents sintering [201,202]. The GO of Au n Ir 4-n , Au n Ir 5-n and Au n Ir 6n clusters was performed with the BPGA, using spin-polarized DFT [138].
The GM structures for Au n Ir 4-n , Au n Ir 5-n and Au n Ir 6-n were found to be planar.

Au-Pd clusters
In a combined theoretical and experimental study (UV-vis depletion spectroscopy), the effect of doping Pd atoms into small cationic gold clusters (tetramer and pentamer) was investigated using the BCGA-DFT approach [203]. It was found that Pd doping leads to the formation of 3D structures. Johnston et al. used the BPGA-DFT approach to perform a systematic search for the GM of neutral Au n Pd m clusters (n + m = 4-10) [204], finding that Pd-rich clusters tend to adopt 3D geometries, whereas 2D geometries are only obtained for gold clusters doped with at most one Pd atom. In a subsequent investigation, medium-sized Au-Pd clusters with 11-18 atoms were globally optimized employing the BPGA with spin polarized DFT [205].

Au-Rh clusters
GO of small Au m Rh n (3 < m + n < 7) clusters was performed using BPGA-DFT [206]. Subsequently, MEGA-DFT calculation was used to investigate the GM structures of Au m Rh n (5 < m + n < 11) clusters [207]. Systems with a high gold concentration tend to form planar structures with Rh atoms in highly coordinated positions and electron transfer takes place from Rh to Au atoms. Clusters with a high concentration of Rh atoms form 3D structures, having an inner Rh core surrounded by Au atoms.

Au-Al clusters
Anionic Au n Al m clusters with (n + m = 7, 8) were studied in a joint approach consisting of experimental photoelectron spectroscopy and BH-DFT optimization [208]. The square bi-pyramidal structural motive of Al 6 − was preferably formed in Au n Al m − clusters, which was explained in terms of the dominance of the stronger Al-Al interactions compared to Au-Au.

Au-Bi clusters
In recent studies, Wang et al. investigated low-lying isomers of anionic AuBi n − (n = 4-8) clusters by combing BH-DFT GO and photoelectron spectroscopy [209]. The energies of the lowest energy DFT isomers were then recalculated at the ab initio CCSD(T) level of theory. Mostly 3D structures were found, with gold at higher coordination locations inside the cluster, and charge transfer from Bi to Au was observed.

Pd-Ir clusters
Noble metals, such as Pd, Pt, Ir and Ru and their alloys, are promising materials for chemical catalysis. GO studies of noble metal nanoalloys play an important role in the design of nanocatalysts for real-world applications. BCGA-DFT was used to explore the PES of Pd n Ir N-n (n = 8-10) clusters [210]. For Ir-rich clusters, cube-based structures are preferred, with more close-packed structures found for Pd-rich clusters. Due to the significantly stronger Ir-Ir bonds, there is strong tendency for Pd-Ir segregation.

Pd-Co clusters
The GO of small Pd n Co m (2 < n + m < 8) nanoalloy clusters was investigated using BCGA-DFT level [211]. Segregation of Pd atoms to peripheral positions in the cluster was observed, which results in a lower surface energy and maximization of strong Co-Co bonds. Pd-induced quenching of the magnetism was also noted.

Pt-Ru clusters
Subnanometre Pt n Ru m (2 < n + m < 9) clusters were studied using BPGA-DFT [212]. It was noted that Pt atoms prefer peripheral sites, while Ru atoms favour central sites, with maximization of the stronger Ru-Ru bonds. The putative GMs are depicted in Figure 7. The authors predicted Ru@Pt core-shell structures to be energetically favoured with increasing cluster size.

Conclusions and outlook
Recent progress in the GO of metal clusters and nanoalloys, applying first principles methods, has been reviewed here, with the main focus on GA-DFT. The use of first principles methods and sophisticated GO algorithms is absolutely necessary to determine the GM structure for very small clusters, because these ultra-small chemical particles tend to form unusual structural motifs, which one would not expect based on chemical intuition. The challenge is locating the GM structure scales with system size and system complexity (bi-or multimetallic, ligated, surface-supported etc.). Hence, the development of more versatile and efficient GO algorithms is essential. For example, the efficiency of GO methods can be increased by combining the advantages of different algorithms or combining GO algorithms with powerful machine-learning tools. The development of computer hardware will make a major contribution to the investigation of more complex systems.
In conclusion, the GO of clusters is a very important task for both basic research and industry, since it can be used to predict experimental structures and can lead to a deeper understanding of experimental data. GO of nanosystems can also be employed for accurate prediction of particle shape and properties and the design of new, improved catalysts, devices, etc. for real-world applications.