A genetic algorithm for first principles global structure optimization of supported nano structures.

A genetic algorithm for first principles global structure optimization of supported nano structures Lasse B. Vilhelmsen and Bjørk Hammer Citation: The Journal of Chemical Physics 141, 044711 (2014); doi: 10.1063/1.4886337 View online: http://dx.doi.org/10.1063/1.4886337 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in A first-principles study on magnetocrystalline anisotropy at interfaces of Fe with non-magnetic metals J. Appl. Phys. 113, 233908 (2013); 10.1063/1.4811685 The role of van der Waals forces in water adsorption on metals J. Chem. Phys. 138, 024708 (2013); 10.1063/1.4773901 First-principles study of the noble metal-doped BN layer J. Appl. Phys. 109, 084308 (2011); 10.1063/1.3569725 First-principles study of metal–graphene interfaces J. Appl. Phys. 108, 123711 (2010); 10.1063/1.3524232 Characterization of methoxy adsorption on some transition metals: A first principles density functional theory study J. Chem. Phys. 122, 044707 (2005); 10.1063/1.1839552

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.102.42.98 On: Sun, 23 Nov 2014 17:34:18

THE JOURNAL OF CHEMICAL PHYSICS 141, 044711 (2014)

A genetic algorithm for first principles global structure optimization of supported nano structures Lasse B. Vilhelmsen and Bjørk Hammera) Interdisciplinary Nanoscience Center (iNANO) and Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark

(Received 25 April 2014; accepted 16 June 2014; published online 29 July 2014) We present a newly developed publicly available genetic algorithm (GA) for global structure optimisation within atomic scale modeling. The GA is focused on optimizations using first principles calculations, but it works equally well with empirical potentials. The implementation is described and benchmarked through a detailed statistical analysis employing averages across many independent runs of the GA. This analysis focuses on the practical use of GA’s with a description of optimal parameters to use. New results for the adsorption of M8 clusters (M = Ru, Rh, Pd, Ag, Pt, Au) on the stoichiometric rutile TiO2 (110) surface are presented showing the power of automated structure prediction and highlighting the diversity of metal cluster geometries at the atomic scale. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4886337] I. INTRODUCTION

Global structure optimization at the atomic scale is of great importance in many fields of computational physics and chemistry. One example of this is the enhanced reactivity and size dependent behavior of nano-sized particles arising from specific sites and non-bulk behavior.1, 2 Any detailed computational study of such phenomena must start from atomic configurations close in energy to the global minimum (GM) structure. Finding this GM structure within a reasonable amount of computational time and human effort is however a significant challenge.3 A set of N atoms require a 3N-dimensional vector for describing the atomic positions. The potential energy of these N atoms as function of their position is called the potential energy surface (PES). Finding the GM structure of a set of atoms is thus equivalent to finding the lowest energy point on the PES. The challenge associated with this is that the PES is a nonlinear function where the number of local minima scales exponentially with N, see Ref. 4, it is therefore only for a very limited number of systems that the PES can be mapped in a systematic manner.5 Several algorithms have been proposed for the unbiased search for the GM on the PES. These algorithms may be divided into thermodynamically motivated ones and evolutionary ones. The thermodynamically driven algorithms can be further divided into molecular dynamics ones and Monte Carlo type ones. The molecular dynamics algorithms follow a temporal evolution of the atoms at a given temperature T where the systematic lowering of the temperature can be used to condensate to the putative GM structure.6 The most popular Monte Carlo type algorithm is the Basin Hopping method which is based on a Metropolis-type selection criterion, random displacements, and relaxations.7 The evolutionary algorithm to be discussed in this paper is in contrast based on a

a) Electronic mail: [email protected]

0021-9606/2014/141(4)/044711/11/$30.00

population of candidates which are evolved through mutual pairing.8 Evolutionary genetic algorithms (GA’s) have successfully been employed at the atomic scale since the early 1990s9 with the largest conceptual developments made by Deaven and Ho who invented the cut-and-splice pairing operator and introduced local relaxations.10 The introduction of these concepts has allowed for the application of GA structure prediction both employing empirical potentials8, 11–14 and also in recent years employing first principles methods15, 16 as reviewed in Ref. 17. The present article discusses the implementation and application of a new GA developed with focus on the use of first principles structure relaxations for the optimization of both clusters in vacuum and supported structures. This implementation has been used successfully for the prediction of putative GM structures of supported metal clusters18–20 and alloys21 and for the reconstruction of extended defects on the rutile TiO2 (110) surface.22, 23 The implementation is expected to be applicable to many other nano-scale problems and it is now publicly available through the open source project “Atomic Simulation Environment” (ASE)24 which is a Python package for executing and analyzing atomic simulations. It should be noted that other publicly available GA codes exist for the optimization of molecular clusters25 and for 3D periodic structures.26, 27 Section II describes the recent successes of the implemented GA. Section III gives a detailed description and benchmark of the GA together with suggested best practices and input parameters. Finally in Sec. IV, the GA will be applied to a small trend study of the adsorption geometry of M8 clusters (M = Ru, Rh, Pd, Ag, Pt, Au) on a stoichiometric rutile TiO2 (110) surface. II. APPLICATIONS OF THE GA

The algorithm presented in this article has been developed as a general purpose GA for the optimization of both 141, 044711-1

© 2014 AIP Publishing LLC


044711-2

L. B. Vilhelmsen and B. Hammer

(a)

(b)

(a)

(b)

(c)

(d)

(e)

(f)

2003: -3.85 eV

1999: -3.59 eV (c)

J. Chem. Phys. 141, 044711 (2014)

(d)

2007: -4.13 eV

2012: -4.41 eV

FIG. 1. Four geometries of Au8 clusters supported on an F-center defect on MgO(100). The publication year and adsorption potential energy are given for each structure. Structures (a), (b), and (c) were reported as putative GM structures in Refs. 28,1, and 29 respectively. The structure in (d) was reported in Ref. 18 and found using the GA presented in this paper.

free standing and supported multi component atomic configurations using any total energy expression, i.e., energy calculator compatible with ASE. The largest structures investigated using a first principles representation of the PES are supported Pd32 clusters20 with most published applications considering between 8 and 20 atoms of one or two atomic types supported on either a surface or within a porous framework.18–23 The first application of the GA was for the determination of the structure of a Au8 cluster supported in an O vacancy on the MgO(100) surface. This investigation highlights the difficulties in manually determining GM structures because three previous publications published in the years from 1999 to 2007 presented three different putative GM structures for this supported cluster.1, 28, 29 Figure 1 illustrates each published structure together with their adsorption potential energies. Notice how the GA predicts a new putative GM structure which has a significantly lower potential energy than the previously proposed structures. Continuing with supported metal clusters, Figs. 2(a) and 2(b) show the result of applying the GA to the investigation of supported metal clusters in metal-organic-frameworks (MOFs) both as small clusters in alloy form (a) and as quite large clusters in single component form (b). Figure 2(c) shows how the GA can also predict the structure of larger clusters supported on a flat surface and Fig. 2(d) how the GA can operate on 1D periodic structures on a support. The studies described until now have all been of structures adsorbed at a support. The GA has however also been used to predict reconstructions of the support itself. Figure 2(e) shows the reconstruction of steps on the rutile TiO2 (110) surface and Fig. 2(f) the growth of 1D defects on the rutile TiO2 (110) surface. III. DESCRIPTION AND BENCHMARK OF THE GA IMPLEMENTATION

Section II established that GA’s can be used to gain insight into the structure and stability of supported nano-

FIG. 2. Examples of prior applications of the GA: (a) Au4 Pd4 supported in the metal-organic-framework MOF-74.20 (b) Pd28 in UiO-66.21 (c) Au24 supported onto the rutile TiO2 (110) surface.30 (d) A Au8 O3 rod along the [001] direction on the rutile TiO2 (110) support.19 (e) An energetically favored re¯ step on rutile TiO .22 (f) Strands grown along the construction of the [110] 2 trough on the rutile TiO2 (110) surface. Here, the GA has been used to predict the terminating structure of these strands (blue and green atoms).23

structures and of surface steps and imperfections. In this section, we want to give a detailed description of the implemented GA while illustrating the importance of different aspects of the algorithm through benchmark calculations. This section is divided into a number of subsections each focusing on different aspects of the algorithm. Figure 3 can be used as a roadmap for this section. For each step a reference is given to the subsection concerned with the given stage. The first step in the GA run is to generate a set of random starting candidates in which the atoms considered are distributed across the relevant part of the super cell. The random starting candidates are then relaxed and the main GA loop starts. The first part of the main loop is to generate the population from the relaxed candidates, each representing local minima in the PES. Next two candidates are selected from the active population and paired together with a possible subsequent mutation. After this offspring candidate has been created it is locally relaxed and eventually included into the set of relaxed candidates. This loop is repeated until the desired number of candidates have been tested or until the user identifies population stagnation. Due to parallelization considerations (see Sec. III F) the GA does not rely on generations where a set of candidates needs to be fully relaxed before a new generation can be generated. This is in contrast to, e.g., the Birmingham cluster genetic algorithm,8 but not unique to our implementation.25


044711-3


J. Chem. Phys. 141, 044711 (2014)

FIG. 4. Three different isomers of Ti6 O12 clusters as found with the GA using a DFTB energy expression. The left most structure is the global minimum and energies are given with respect to this structure.

FIG. 3. Overview of the steps (left) and data structure (right) in the genetic algorithm. For each step a reference is given to a section in the text with details about the given step.

A detailed description of each step in Figure 3 is given in the following together with a discussion of how tuning different optimization parameters can influence the GA performance. In the rest of this section, we will be concerned with benchmarking the different aspects of the GA in order to assess their relative importance. While the presented GA has been developed for use with first principles energy calculations, a systematic benchmark of the algorithm using, e.g., density functional theory (DFT) calculations would be too computationally demanding, and we therefore opt to employ Density Functional Tight Binding (DFTB) calculations.31 Even though DFTB calculations are vastly faster than DFT calculations they are still too costly to test the GA on supported structures. Consequently, a more tractable system of a two-component cluster in vacuum, specifically (TiO2 )6 aka Ti6 O12 , was chosen. The DFTB calculations are done with the DFTB+ program package32 and bond parameters are taken from Ref. 33. These parameters have been optimized for small organic molecules with transition metal centers, so full agreement with DFT calculations is not expected. This is, however, not an issue since the goal of this study is not to find the structure of these clusters, but rather to benchmark the performance of the GA with a potential having many well-defined local minima in the PES. Ti6 O12 clusters have a convenient size for this investigation since the global minimum on average can be found with the GA in about 500 tested candidates as will be documented

below. Figure 4 shows the global minimum and two other stable isomers of this cluster. To estimate how small a portion of the PES the global minimum occupies 10 000 random starting configurations were generated and locally relaxed. Only one of these 10 000 configurations relaxes into the GM structure. This indicates that the global minimum is not easily found and that the GA does a satisfying job of finding it. Figure 5 shows a histogram of the energy of these 10 000 randomly generated structures. Figure 5 also shows the energy of all structures encountered by the GA from 200 separate runs. These data will be discussed in the following. Notice how the distribution moved towards much lower energies for the GA data. Since the GA is based on a set of random starting candidates and on a stochastic evolution, each run of the GA will produce a unique set of candidates and each run will have a unique way to progress towards the global minimum. For this reason the GA can only be benchmarked by analyzing statistically across many separate runs. From a collection of separate GA runs the following information will be extracted: the average number of candidates needed to find the GM structure for the successful runs, success-rates, and the number of candidates needed before 80% of the GA runs have succeeded. Figure 6 shows the histogram obtained using default parameters for the optimization: a population size of 20 candidates, a set of 20 starting candidates, and no mutations. The histogram is based on 200 separate runs of the GA. Some of

FIG. 5. The energy distribution of all locally relaxed structures encountered from a random generation of structures (top) and from runs of the GA (bottom).


044711-4


J. Chem. Phys. 141, 044711 (2014)

processes although it is possible to do so. The candidates are generated by defining a simulation cell and distributing the atoms randomly within this cell. To ensure that no single atom is detached from the main cluster each atom is required to be closer than dbl,max to another atom. Similarly, no two atoms are allowed to be closer than dbl,min since overlapping atoms are unphysical and give rise to various computational difficulties. The distances dbl are dbl,max (t, o) = αmax · [rcov (t) + rcov (o)]

(1)

dbl,min (t, o) = αmin · [rcov (t) + rcov (o)],

(2)

and

FIG. 6. Histogram of the number of candidates needed before the global minimum is found for a GA run with 20 random starting configurations, a population size of 20, standard comparisons criteria, and no mutations. The red line marks 80th percentile of required candidates, i.e., the value at which 80% of the independent GA runs have found the GM structure. A successful run is when the GM is found with less than 3000 tested candidates.

the GA runs find the global minimum with very few candidate structures, whereas a few runs use more than 1500 candidates before the best structure is found. Only runs finding the GM within 3000 candidates will be considered successful and runs not fulfilling this are abandoned. Figure 6 reports a success rate of 98% and considering these runs alone, an average of λ = 468 ± 29 candidates are found to be needed in order to determine the global minimum structure. The red line marked in Fig. 6 marks the 80th percentile, γ 80 , of the distribution, i.e., the number of candidates needed before 80% of the GA runs have succeeded. This is an interesting number because it combines both the success-rate and the average number of structures needed for the successful runs to have completed. If one could estimate this number a priori one could opt to do a small number, β, of independent GA runs that are terminated when reaching γ 80 candidates. After the β independent runs the probability of having found the GM structure would be 1 − (1 − 0.8)β . With the present example, one would thus have a 99.9% chance of having found the GM structure after performing β = 5 independent runs of the GA, where each GA run had been stopped after γ 80 = 643 tested candidates. Histograms similar to the one in Fig. 6 will be provided in the following to gauge the impact of different parameter choices. When gauging the GA performance the reference will always be the one in Fig. 6.

where rcov (o) and rcov (t) are the covalent radii of atoms of types o and t, respectively. Typical values of αmin and αmax are 0.8 and 1.6, respectively. Before the main loop in Fig. 3 is initiated each random starting candidate is relaxed to its local minimum and the population always consists solely of relaxed candidates. It is important that the atoms of the starting candidates are distributed within a box of reasonable size. For the results in Fig. 6 the atoms of the starting candidates are distributed in a cube with an edge length of 20 Å. Figure 7 shows the consequence of instead generating the starting candidates from randomly positioned atoms in a cuboid of dimensions 4 Å × 4 Å × 20 Å (the GM structure has dimension of approximately 6 Å × 6 Å × 6 Å). Notice how the GA performance is significantly reduced by this. This is because the starting candidates are too elongated leading to an initial relaxed population consisting of structures far in shape from the GM structure. B. Maintaining a diverse population

While the GA is running, a population is maintained which consists of the N most stable and structurally different candidates. For the reference run in Fig. 6 the population size (a)

(b)

A. The random starting candidates

Before starting the main loop in Fig. 3 one needs to generate a set of random starting candidates. There are in principle no restrictions on the size and diversity of the set of starting candidates, but it is advantageous to choose a sufficiently large and diverse set of candidates such that many different types of candidates are represented from the start. A typical number of starting candidates is 20. All starting candidates are generated before initializing the main loop and new configurations are not introduced during the optimization

FIG. 7. The effect of generating the random starting population in a box of unreasonable proportions compared to the problem investigated. The result in (a) is reproduced from Fig. 6. In (b) the starting population is put in a 4 Å × 4 Å × 20 Å box.


044711-5


J. Chem. Phys. 141, 044711 (2014)

(a)

(b)

(c) FIG. 9. Two different structures of Au8 adsorbed on an MgO(100) F-site defect drawn together with the sorted distance list Di for each structure. The energy difference between the two structures is only 0.7 meV. Their structural properties are δrel = 0.05 and dmax = 1.03 Å. dmax is marked in the graph.

FIG. 8. The effect of changing the population size. The data in (b) are the same as in Fig. 6.

is fixed at 20 candidates. Figure 8 shows the effect of having either a smaller or larger population size. Notice how reducing the population size actually improves the speed at which the GM is found but at the cost of an increased failure rate. Doubling the population size increases the number of candidates needed to find the GM structure at the cost of slower convergence. The population size dependence in Fig. 8 shows how the GA performance is a trade off between speed and reliability. When the population is small, new improved configurations get to proliferate faster due to the smaller degree of competition and the GA can thereby progress more quickly. A small population does, on the other hand, not include as much diversity as a larger population, and there is therefore an increased risk that the population does not include the genetic information needed to generate the optimal structure. To help ensure that the population is genetically diverse, measures must be taken to avoid multiple occurrences of one structure in the population. A natural starting point for the comparison of two structures is their energy. Only employing an energy comparison can however be dangerous because two configurations can be very different and still be close in energy, or they can be similar in structure but have energies, which differ by more than some energy criterion, E. Figure 9 shows an example of two structures with an energy difference of only E = 0.7 meV but with very different structures. In the literature, several geometric comparisons are proposed; one can characterize clusters by their moment of inertia, radial distribution function, average and variance of atomic distances or in bulk by, for instance, the fingerprint function.34 Common for these comparisons is that they must be invariant under translations and permutations of identical atoms.

In the present GA implementation we ensure that two structures are considered different by employing an interatomic comparison criteria. A sorted list, Di , of all interatomic distances is calculated for all candidates i (see Fig. 9). Two configurations i and j are considered equal if k |Di (k) − Dj (k)| < δrel (3) 1 k Di (k) + Dj (k) 2 and maxk (|Di (k) − Dj (k)|) < dmax .

(4)

k indexes each entry in the list Di . Each new candidate i is compared to all relaxed candidates j in the population. The first criterion is the relative accumulated difference between the two candidates and the second criterion is the maximum difference between two distances in the two candidates. Typical values are δrel = 0.03 and dmax = 0.7 Å. Employing these criteria makes the two candidates in Fig. 9 distinct. If an optimization includes different atomic types, the structural criteria are calculated for each atomic type separately and added together by a weighted sum, depending on how many distances exist for each atom type. Figure 10 shows the effect of not including comparisons or only employing either energy or structure in the comparisons. As seen from the figure the GA performs very poorly if population diversity is not ensured (Fig. 10(b)), but it also shows that employing only one of the criteria already improves the GA performance greatly (Figs. 10(c) and 10(d)), with the biggest improvement being gained by the structural comparison. Comparing Figs. 10(a) and 10(d) it is evident that the energy comparison actually decreases both the speed and reliability of the GA showing that for this application of the GA, comparisons by energy should be avoided. The GA is implemented such that the user can easily employ or implement alternative comparisons. C. Selection of candidates

To impose a survival of the fittest scheme, the most stable candidates must be selected for pairing more frequently


044711-6


J. Chem. Phys. 141, 044711 (2014)

(a)

(a)

(b)

(b)

(c)

(c)

(d)

FIG. 11. (a) The same data as in Fig. 6. (b) The effect of not including Ui in the fitness function. (c) The effect of not including Fi in the fitness function.

FIG. 10. (a) Both comparisons are included. (b) No comparisons are included. (c) Only energy comparisons. (d) Only structural comparisons.

than their less stable counterparts. In Ref. 35, it is proposed to assign a fitness (Fi ) to each candidate in the population and to select candidates with a probability (Pi ) proportional to this fitness. To help ensure a diverse population we choose to multiply Fi with a uniqueness factor Ui that tends to zero the more times the ith candidate has participated in pairings or the more structures exist that are similar to the ith structure but outside of the population. Fi and Ui are defined as Fi =

1 Ei − Emin [1 − tanh(2ρi − 1)], with ρi = , (5) 2 Emax − Emin

Ui =

1 1 · , 1 + ni 1 + mi

(6)

with Pi = Fi · Ui . Emin (Emax ) is the minimum (maximum) energy of any structure in the population. ni is the number of times candidate i has participated in a pairing and mi is the number of candidates which are structurally similar to the ith candidate, but not included in the population because the ith candidate is already present. The functional form of Ui is constructed such that unique new structures do not have their fitness reduced and old candidates are not penalized too heavily.21 Figure 11 shows the GA performance if either Ui or Fi is left out of the fitness function. Notice how the success rate is reduced if Ui is left out and the performance is slightly reduced if Fi is left out. This shows that avoiding overuse of the same structures indeed improves the success-rate of the GA although at the cost of slightly worse performance. It is

surprising that leaving Fi out of the fitness function only has a minor detrimental impact on the overall GA performance.

D. Pairing of candidates

Arguably, the most important ingredient in a successful GA is the pairing operator. Two candidates need to be paired together in a manner which has a high likelihood of preserving the best part of each, while ensuring that a large and relevant portion of the PES is explored. In a general purpose GA, the data describing the candidates are represented as a string of information. When pairing two candidates together this string of characters is cut into smaller bits and parts from each candidate are selected. Deaven and Ho realized that an important property of a cluster of atoms is the local environment of each atom.10 When doing a pairing of two atomic clusters one should therefore not disregard this local behavior. They thus devised the cutand-splice operator depicted in Fig. 12. This pairing employs a randomly oriented dividing surface through the two candidates’ common center of mass which cuts the two clusters in halves. In its simplest form, the dividing surface is simply a cutting plane, but it may also be more elaborate, e.g., a sphere creating a core-shell splitting.36–38 When doing a pairing as depicted in Fig. 12 one often ends up with a candidate with too few or too many atoms, or with a structure where two atoms are prohibitively close. We choose the following procedure when ensuring the preserved stoichiometry of pairing candidate c1 and c2 consisting of N atoms: 1. Choose an orientation of the random cutting plane and position it in the common center of mass of c1 and c2 .


044711-7


J. Chem. Phys. 141, 044711 (2014)

(a)

(b)

FIG. 12. The cut-and-splice operator used on a supported Au8 cluster. The coloring is chosen to illustrate which configuration each atom is selected from. The black line indicates the cutting plane.

2. Choose all atoms from c1 which are to the left of the plane and all atoms from c2 which are to the right of the plane. 3. If the total number of atoms in the new cluster (Nn ) is larger than N remove the Nn − N atoms which are farthest away from the cutting plane. 4. If Nn < N randomly choose N − Nn atoms from the two candidates which were not previously included. 5. When a cluster with N atoms has been created check all interatomic distances in the new cluster. If any distance is smaller than dbl,min discard the new candidate and try with a new cutting plane. Otherwise the new candidate is accepted. If the cluster being optimized consists of M different types of atoms the above procedure is done separately for each atom type, and only step 5 is done with the coordinates of all atoms in the new candidate. When optimizing a cluster near a surface the surface is disregarded during the cut-and-splice operation. The cut-and-splice operator was developed to work on clusters in vacuum consisting of only one species. One could consequently think that the cut-and-splice operator would work poorly for multi-component and supported systems. In the literature, however, there are reports of GAs used successfully for the optimization of step reconstructions39 and the optimization of alloy clusters.35 Common for these previous uses is the fact that the cut-and-splice operator is not changed in any significant way. Figure 13 shows how one can easily implement a cutand-splice operator that performs sub-optimally. In Fig. 13(b), the atoms closest to the cutting plane are removed instead of the ones farthest away when too many atoms have been selected. This decreases the success rate of the GA to 90% from 98%. Figure 13(c) shows the effect of rotating one of the candidates being paired prior to the pairing. This clearly has a huge impact on the GA performance. It can be explained by the population evolving towards structures being aligned a certain way in the computational cell. If one of the candidates is rotated at each step this evolution is hindered. For clusters with a small degree of directional binding, like Lennard-Jones clusters, it has, however, been shown that rotating one of the candidates is advantageous.10

(c)

FIG. 13. (a) and (b) The performance when the atoms closest to the cutting plane are removed instead of the atoms farthest away from the cutting plane when too many atoms have been chosen. (c) The performance if one of the candidates used to generate a new structure is rotated randomly before being paired.

E. Mutation of candidates

Mutations are introduced to ensure diversity, and to avoid that the GA stagnates in a local minimum region of the PES with a too degenerate population. Whereas the cut-and-splice operator has not been changed significantly since its introduction, the number of mutations depending on application has expanded.8, 10, 18, 21, 26, 35, 36, 40 Mutations are applied after a new candidate is generated but before it undergoes local relaxation. Typically 30% of the generated candidates are mutated. The following mutations have been implemented in our code:

r Rattle: 40% of the atoms in the cluster are translated a random distance in a random direction. The translated distance is between 0 Å and 0.8 Å and is chosen differently for each atom.40 r Twist: For clusters adsorbed on a surface the cluster can be rotated a random angle with respect to the surface normal. The rotation point can be a point defect or the center of mass of the cluster depending on the type of support.8 r Permutation: The atomic number of one-third of the atoms in the cluster is permuted. This mutation is reported to strongly improve the convergence for multicomponent systems since it rearranges the atoms while preserving the local geometry.35 r Mirror: A randomly oriented plane is defined through the center of mass of a cluster. One-half of the cluster is selected and mirrored across the plane. To preserve the stoichiometry atoms are added or removed following the same principles as when two clusters are paired


044711-8


(a)

(b)

(c)

(d)

FIG. 14. (a)–(c) An increasing number of the generated structures are mutated. (d) Instead of mutating structures generated from a pairing a single structure is taken directly from the population.

together. A variation of this operator mirrors the selected half a second time through a plane perpendicular to the first. This is equivalent to rotating the selected half 180◦ instead of mirroring.21 Figures 14(a)–14(c) show the performance impact of mutating 0%, 30%, and 100% of the proposed candidates with either the rattle, permutation, or mirror mutation. For these runs of the GA a new candidate is first generated using the cut-and-splice pairing operator and the new candidate is then subsequently mutated. Figure 14(d) shows the performance if instead a candidate is taken directly from the population and mutated. Figure 14 shows that mutations only have a minor impact on the GA performance and reliability for this particular problem except if 100% of the generated structures are mutated. This can be understood by considering mutations as a means of introducing diversity to the population. If all candidates are mutated too much noise is introduced and the population therefore takes longer to converge to the GM structure. Since mutating 30% of the proposed candidates has a negligible effect on the GA performance this might be interpreted as the GA not needing additional measures to ensure diversity for this particular problem. The difference between applying mutations to candidates directly from the population (Fig. 14(d)) or to candidates generated from a pairing (Fig. 14(b)) is within the uncertainty of the data. It might be that this difference becomes larger for problems benefiting more from mutations.

J. Chem. Phys. 141, 044711 (2014)

F. Local relaxations: Parallelization and the computational cost

Since DFT calculations are computationally demanding to perform several strategies have been proposed in order to reduce the computational cost. For clusters in vacuum many studies employ empirical potentials for the structure optimization and only use DFT for the following data analysis.10, 40–43 In studies including a support, the modeling of the support usually requires the use of DFT, although some examples with empirical potentials exist for the interaction with defect free supports.44 Since our research has focused on supported structures and reconstructions in complex oxide systems we cannot avoid employing DFT optimizations for the entire GA process.18–23 We have therefore employed two different approaches which reduce the net computational cost. The first avenue to reducing the computational cost employs a two-step optimization technique where the proposed structures are first locally relaxed with the use of a localized basis set and subsequently with a full real-space grid or planewave basis.18, 21 A comparison of relaxing a structure directly or through such a two-step optimization scheme is shown in Fig. 15. Employing this technique gives a speed-up of around a factor three for one particular relaxation. Also by comparing the relaxed structure to all previously encountered structures before doing the costly relaxation one can further reduce the computational cost giving a net speedup, for this particular system, of a factor of 7 compared to doing all relaxations employing the full basis set. The other scheme for reducing the computational cost is to accept some errors in the structure optimization process by using smaller super cells and optionally worse basis sets.45, 46 In the studies described in Refs. 22 and 23 we investigated adsorbates and reconstructions on the rutile TiO2 (110) surface. Properly converged calculations on this support require modeling the support by at least a three tri-layer thick slab and using a setup for Ti describing 12 valence electrons. Such calculations are too expensive for an investigation of many structures. Instead we used a two tri-layer thick surface and a 4 valence electron setup for Ti together with a dzp LCAO ba-

FIG. 15. Comparison of locally relaxing a Au8 cluster on top an F-center on the MgO(100) surface. In the top row, the structure is first relaxed with a local basis set and subsequently with a full grid-basis. In the bottom row, the structure is relaxed directly with the full basis. The computational time for each step is indicated as the product of the wall time and the number of cores used.


044711-9


J. Chem. Phys. 141, 044711 (2014)

FIG. 16. Flow diagram showing how several structures are relaxed simultaneously across the computer cluster.

sis set. Running the GA with these lower precision parameters introduces some errors, but the correlation between the two levels of precision is sufficiently good (see Fig. 3 of Ref. 22), that a post relaxation of all structures, within some threshold of the most stable structure proposed by the GA, has a high chance of including the most stable structure. Because first principles structural relaxations do not achieve perfect scaling one further optimization that can be done is to relax many structures independently and in parallel. Figure 16 shows how the GA software monitors the number of running relaxations and resubmits new structures for local relaxation as needed. For a typical GA run, we perform structural optimization of up to 20 candidates concurrently on separate computing modes. This parallelization opportunity is the reason that generations are not employed in this GA implementation since waiting for a full generation to be relaxed would make the parallelization less effective. IV. GA ASSISTED TREND STUDIES

To conclude this article, we want to apply the GA to perform an automated trend study to show the power of the method. The beauty of automated structure optimization is that it minimizes the need for direct human influence in the structure prediction process. If sufficient computer resources are available it thus allows for detailed trend studies across different stoichiometries and sizes without much human effort. Figure 17 shows the result of applying the GA to the search for putative GM structures for six different M8 transition metal clusters adsorbed on a flat stoichiometric rutile TiO2 (110) surface. This study has been done using the same computational setup as described for the TiO2 studies in

FIG. 17. The putative GM structures of M8 clusters supported on a stoichiometric rutile TiO2 surface. The adhesion potential energies are with reference to metal clusters in vacuum.11, 47–51

Sec. III F and in further detail in Refs. 19 and 30 with the use of the PBE functional. Between three and five structures for each metal are selected and relaxed with the more precise setup. This is because of the possible discrepancies between the less precise setup used for the GA search and the more precise calculations presented in Fig. 17. The GA calculations are all done spin-paired due to computational constraints, but subsequent analysis has shown only small energy changes from applying spin polarization. The energetics in Fig. 17 are with respect to the corresponding cluster in its putative GM structure in vacuum.11, 47–51 The structure of adsorbed Pt8 is in agreement with Ref. 52 whereas our structure of adsorbed Ag8 differs from a previously manually proposed putative GM structure.53 The six putative GM structures shown in Fig. 17 have been established by allowing the GA to test N = 300 structures for each metal type. The choice of N was made based on five independent runs of the GA for adsorbed Pt8 . Figure 18(a) shows the energy evolution for these five independent GA runs. The figure shows that after 150 tested structures all five runs of the GA have encountered a structure with a very small energy difference to the putative GM structure. Allowing the GA to test N = 300 structures should therefore with a high probability be sufficient for finding the putative GM structure. This is confirmed from Fig. 18(b) for the other metals where after about 150 tested structures no new structures are encountered. The only exception to this is Au where


044711-10


(a)

(b)

FIG. 18. The energy of the most stable structure encountered relative to the putative GM structure as the GA is progressing. The energy is shown with respect to the energy of the most stable structure encountered at the end. (a) shows the evolution from five independent GA executions for Pt whereas (b) shows the data for the other five metals considered.

a new putative GM structure is found as late as at around 250 tested structures. There is however only a minute structural difference between this structure and the second most stable one. First doing a number of GA runs for one stoichiometry to establish N is the more computationally feasible alternative to the method described in the first part of Sec. III. When doing first principles calculations it is not possible to do enough runs that γ 80 can be determined within a reasonable precision allowing for a more systematic approach. From Fig. 17, it is clear that the six metals behave very differently at the atomic scale. Only Au8 and Ru8 maintain their vacuum structure when adsorbed on the surface whereas the other metals change their internal structure quite drastically in order to coordinate optimally to the rutile surface. We stress that the great structural diversity of the M8 clusters would be next to impossible to predict without the use of automated structure prediction methods. A further discussion of the electronic and geometric properties of the presented clusters is outside the scope of this article.

V. CONCLUSIONS

This article describes the implementation of a genetic algorithm for global structure optimization at the atomic scale. The algorithm has specifically been implemented for use with computationally demanding first principles local relaxation techniques and the implementation has proven useful for a number of applications. Section III described the algorithm in detail and benchmarked various parameters within the implementation using statistical averages from repeated independent runs. The benchmarks show a sensitivity to the choice of population size, comparison criteria, mutations, and pairing operator.

J. Chem. Phys. 141, 044711 (2014)

Sensible starting values are however given and discussed providing the user with a good starting point for employing the algorithm through the “Atomic Simulation Environment.”24 Shown in the final part of this article is a small trend study of the adsorption of transition metal clusters onto the rutile TiO2 (110) surface. For this study, we developed a technique for deciding when to stop the GA: First, we do multiple GA runs for one stoichiometry to estimate the complexity of the problem and the number of structures N to test. Subsequently, we do one single run for each other stoichiometry testing N structures for each. Such a trend study would have been a demanding task to undertake without the use of an automated structure prediction method, but it becomes a standard task with this GA at hand. In conclusion, we expect that global structure optimization algorithms will continue to aid researchers doing atomic scale modeling. We believe the use of such methods will increase in the coming years together with the ever increasing performance of computer clusters. ACKNOWLEDGMENTS

This work received support from the Lundbeck Foundation, The Danish Council for Independent Research, COST action CM1104, Nordforsk, and the Danish Center for Scientific Computing. 1 A.

Sanchez, S. Abbet, U. Heiz, W. D. Schneider, H. Häkkinen, R. N. Barnett, and U. Landman, J. Phys. Chem. A 103, 9573 (1999). 2 A. T. Bell, Science 299, 1688 (2003). 3 R. L. Johnston, Atomic and Molecular Clusters (Taylor & Francis, 2002). 4 F. Stillinger, Phys. Rev. E 59, 48 (1999). 5 D. J. Wales, M. A. Miller, and T. R. Walsh, Nature 394, 758 (1998). 6 R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 7 D. J. Wales and J. P. K. Doye, J. Phys. Chem. A 101, 5111 (1997). 8 R. L. Johnston, Dalton Trans. 22, 4193 (2003). 9 B. Hartke, J. Phys. Chem. 97, 9973 (1993). 10 M. D. Wolf and U. Landman, J. Phys. Chem. A 102, 6129 (1998). 11 B. Assadollahzadeh and P. Schwerdtfeger, J. Chem. Phys. 131, 064306 (2009). 12 F. Baletto and R. Ferrando, Rev. Mod. Phys. 77, 371 (2005). 13 K.-M. Ho, A. A. Shvartsburg, B. Pan, Z.-Y. Lu, C.-Z. Wang, J. G. Wacker, J. L. Fye, and M. F. Jarrold, Nature 392, 582 (1998). 14 A. Rapallo, G. Rossi, R. Ferrando, A. Fortunelli, B. C. Curley, L. D. Lloyd, G. M. Tarbuck, and R. L. Johnston, J. Chem. Phys. 122, 194308 (2005). 15 A. N. Alexandrova and A. I. Boldyrev, J. Chem. Theory Comput. 1, 566 (2005). 16 H. Xiang, S.-H. Wei, and X. Gong, J. Am. Chem. Soc. 132, 7355 (2010). 17 S. Heiles and R. L. Johnston, Int. J. Quantum Chem. 113, 2091 (2013). 18 L. B. Vilhelmsen and B. Hammer, Phys. Rev. Lett. 108, 126101 (2012). 19 L. B. Vilhelmsen and B. Hammer, J. Chem. Phys. 139, 204701 (2013). 20 L. B. Vilhelmsen and D. S. Sholl, J. Phys. Chem. Lett. 3, 3702 (2012). 21 L. B. Vilhelmsen, K. S. Walton, and D. S. Sholl, J. Am. Chem. Soc. 134, 12807 (2012). 22 U. Martinez, L. B. Vilhelmsen, H. H. Kristoffersen, J. Stausholm-Moller, and B. Hammer, Phys. Rev. B 84, 205434 (2011). 23 R. Bechstein, H. H. Kristoffersen, L. B. Vilhelmsen, F. Rieboldt, J. Stausholm-Moller, S. Wendt, B. Hammer, and F. Besenbacher, Phys. Rev. Lett. 108, 236103 (2012). 24 S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng. 4, 56 (2002). 25 J. M. Dieterich and B. Hartke, Mol. Phys. 108, 279 (2010). 26 A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006). 27 W. W. Tipton and R. G. Hennig, J. Phys.: Condens. Matter 25, 495401 (2013). 28 J. Wang and B. Hammer, Top. Catal. 44, 49 (2007). 29 H. Häkkinen, S. Abbet, A. Sanchez, U. Heiz, and U. Landman, Angew. Chem., Int. Ed. 42, 1297 (2003).


044711-11 30 L.


B. Vilhelmsen and B. Hammer, ACS Catal. 4, 1626 (2014). Koskinen and V. Mäkinen, Comput. Mater. Sci. 47, 237 (2009). 32 B. Aradi, B. Hourahine, and T. Frauenheim, J. Phys. Chem. A 111, 5678 (2007). 33 G. Zheng, H. A. Witek, P. Bobadova-Parvanova, S. Irle, D. G. Musaev, R. Prabhakar, K. Morokuma, M. Lundberg, M. Elstner, C. Köhler, and T. Frauenheim, J. Chem. Theory Comput. 3, 1349 (2007). 34 A. R. Oganov and M. Valle, J. Chem. Phys. 130, 104504 (2009). 35 L. D. Lloyd, R. L. Johnston, and S. Salhi, J. Comput. Chem. 26, 1069 (2005). 36 N. S. Froemming and G. Henkelman, J. Chem. Phys. 131, 234103 (2009). 37 Z. Chen, X. Jiang, J. Li, and S. Li, J. Chem. Phys. 138, 214303 (2013). 38 S. Lysgaard, D. Landis, T. Bligaard, and T. Vegge, Top. Catal. 57, 33 (2014). 39 R. M. Briggs and C. V. Ciobanu, Phys. Rev. B 75, 195415 (2007). 40 D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288 (1995). 41 A. Bruma, R. Ismail, L. Oliver Paz-Borbon, H. Arslan, G. Barcaro, A. Fortunelli, Z. Y. Li, and R. L. Johnston, Nanoscale 5, 646 (2013). 31 P.

J. Chem. Phys. 141, 044711 (2014) 42 D.

M. Daven, N. Tit, J. R. Morris, and K. M. Ho, Chem. Phys. Lett. 256, 195 (1996). 43 R. Ferrando, A. Fortunelli, and R. L. Johnston, Phys. Chem. Chem. Phys. 10, 640 (2008). 44 F. R. Negreiros, G. Barcaro, Z. Kuntová, G. Rossi, R. Ferrando, and A. Fortunelli, Surf. Sci. 605, 483 (2011). 45 G. Barcaro and A. Fortunelli, Chem. Phys. Lett. 457, 143 (2008). 46 G. Barcaro, E. Aprà, and A. Fortunelli, Chem. – Eur. J. 13, 6408 (2007). 47 V. Bonaci´ c-Koutecký, L. CeSpiva, P. Fantucci, and J. Koutecký, J. Chem. Phys 98, 7981 (1993). 48 G. Zanti and D. Peeters, Eur. J. Inorg. Chem. 2009, 3904. 49 A. Sebetci, Comput. Mater. Sci. 78, 9 (2013). 50 Y.-C. Bae, H. Osanai, V. Kumar, and Y. Kawazoe, Phys. Rev. B 70, 195413 (2004). 51 W. Zhang, H. Zhao, and L. Wang, J. Phys. Chem. B 108, 2140 (2004). 52 D.-e. Jiang, S. H. Overbury, and S. Dai, J. Phys. Chem. C 116, 21880 (2012). 53 V. E. Matulis, A. S. Mozheiko, and O. A. Ivashkevich, Russ. J. Gen. Chem. 80, 1068 (2010).


A Novel Hybrid Firefly Algorithm for Global Optimization.

A local stability supported parallel distributed constraint optimization algorithm.

Genetic Algorithm (GA)-Based Inclinometer Layout Optimization.

Optimization of genomic selection training populations with a genetic algorithm.

First-principles studies on graphene-supported transition metal clusters.

An adaptive hybrid algorithm based on particle swarm optimization and differential evolution for global optimization.

Protein structure optimization with a "Lamarckian" ant colony algorithm.

First-principles predicted low-energy structures of NaSc(BH4)4.

Observable-targeting global cluster structure optimization.

First principles molecular dynamics without self-consistent field optimization.

A Novel Hybrid Clonal Selection Algorithm with Combinatorial Recombination and Modified Hypermutation Operators for Global Optimization.

A Guiding Evolutionary Algorithm with Greedy Strategy for Global Optimization Problems.

A hybrid cuckoo search algorithm with Nelder Mead method for solving global optimization problems.

Genetic algorithm application in optimization of wireless sensor networks.

A Swarm Optimization Genetic Algorithm Based on Quantum-Behaved Particle Swarm Optimization.

A cuckoo search algorithm for multimodal optimization.

Multi-objective genetic algorithm for the optimization of a flat-plate solar thermal collector.

A newton cooperative genetic algorithm method for in silico optimization of metabolic pathway production.

A Combination of Genetic Algorithm and Particle Swarm Optimization for Vehicle Routing Problem with Time Windows.

Novel genetic algorithm search procedure for LEED surface structure determination.

An adaptive genetic algorithm for crystal structure prediction.

Research on laser marking speed optimization by using genetic algorithm.

A genetic algorithm based multi-objective shape optimization scheme for cementless femoral implant.

Accurate and scalable O(N) algorithm for first-principles molecular-dynamics computations on large parallel computers.