The Birmingham parallel genetic algorithm and its application to the direct DFT global optimisation of Ir N ( N = 10 – 20) clusters †

A new open-source parallel genetic algorithm, the Birmingham parallel genetic algorithm, is introduced for the direct density functional theory global optimisation of metallic nanoparticles. The program utilises a pool genetic algorithm methodology for the e ﬃ cient use of massively parallel computational resources. The scaling capability of the Birmingham parallel genetic algorithm is demonstrated through its application to the global optimisation of iridium clusters with 10 to 20 atoms, a catalytically important system with interesting size-speci ﬁ c e ﬀ ects. This is the ﬁ rst study of its type on Iridium clusters of this size and the parallel algorithm is shown to be capable of scaling beyond previous size restrictions and accurately characterising the structures of these larger system sizes. By globally optimising the system directly at the density functional level of theory, the code captures the cubic structures commonly found in sub-nano-metre sized Ir clusters.


Introduction
Nanosized materials are currently being investigated for potential use in a variety of applications.This is because the nanosizing effects seen in such materials result in properties different from those of the bulk material.These properties can also be tuned, normally through altering the size and shape of the cluster.
Metallic nanoparticles are such materials, with potential optical, magnetic and catalytic applications. 1Small Ir nanoparticles, in particular, are currently used as catalysts for a range of organic reactions, including olefin hydrogenation, oligomerisation, and ring-opening of cycloalkanes. 2Ir has been shown both experimentally 3 and theoretically 4 to exhibit significant nanosize-induced hydrogen adsorption capacity.Larger Ir nanoparticles have been shown to be active in C-C bond hydrogenolysis. 5Selective molecular recognition has also been seen in supported Ir cluster-based catalysts. 6 key step in rationalising properties, such as the catalytic activity of nanoparticles, is their structural characterisation.To achieve this it is necessary to sample comprehensively the potential energy landscape (PES) of the nanoparticle.A wide variety of methods is available for the exploration of the PES.These methods include statistical mechanical methods, such as the CBEV/FCEM method, 7,8 basin-hopping methods, 9 such as GMIN, 10 and genetic algorithms, such as the Birmingham cluster genetic algorithm (BCGA). 11The choice of method largely depends on the size and complexity of the system.
It is necessary to decide the level of theory required to replicate accurately a particular PES of a system.For example, the electronic structure of larger nanoparticles is thought to resemble closely that of the bulk material.This means the use of empirical potentials, such as the Gupta potential, 12 is suitable for the accurate representation of the PES.Statistical mechanical methods, such as CBEV/FCEM, may be best suited to these larger systems. 84][15] A variety of methods has been developed to achieve this, many of which have been outlined by Heiles and Johnston. 16he cubic structures adopted by sub-nanometre Ir clusters have been previously shown up to the CCSD(T) level of theory, [17][18][19][20][21][22][23][24] differing from the fcc structure of the bulk material.It is therefore vital that any global optimisation of Ir N (N = 10-20) structures is performed directly at the density functional theory (DFT) level of theory at least.
A quantum description of the PES greatly increases the cost of exploring it comprehensively, limiting the size of the cluster it is possible to investigate. 16It is therefore necessary that efficient parallel methodologies, with the ability to utilise greater computational resources, are developed.][27][28][29] This work presents the global optimisation of Ir N (N = 10-20) clusters directly at the DFT level of theory.This is achieved using the Birmingham parallel genetic algorithm (BPGA), a new open-source genetic algorithm available via Bitbucket. 30The BPGA utilises a pool genetic algorithm methodology combined with the evaluation of potential cluster geometries in parallel. 25This combination ensures highly efficient scaling when compared with generation based genetic algorithms and allows the structural characterisation of larger and more complex systems.The pool methodology has been recently applied to metallic clusters 31 and was benchmarked and applied successfully to the global optimisation of the much studied Au 10 and Au 20 clusters. 31Its predictions for Ag 10 + have also been shown to be accurate when compared with spectra from molecular beam experiments. 32The BPGA incorporates this highly efficient algorithm within a flexible Python framework.Due to the 5d 7 6s 2 ground state and other low lying states originating from its 5d 8 6s 1 configuration, 33 the spin of the Ir N clusters must be considered in the calculations. 33To account for this, spin-polarised DFT global optimisations are performed.The use of spin-polarised local minimisations effectively doubles the computational cost and can only now be performed due to the parallelism of the BPGA.
The BPGA is also capable of globally optimising bimetallic nanoalloys, whose PES is complicated by the presence of homotops. 34It is hoped that this work demonstrates that scaling capability of the BPGA and its ability to utilise massively parallel architectures, which enable the program to predict accurately the geometries of metallic nanoparticles.

Birmingham parallel genetic algorithm
The BPGA is a parallel genetic algorithm for the structural characterisation of nanoparticles.The program is written in object-oriented Python.This allows greater flexibility and the ability to utilise the large existing libraries of Python code, such as the atomic simulation environment. 35Python is well suited to job submission, required by the DFT interface, on a large shared HPC resource, such as ARCHER. 36The program is open-source and available via BitBucket. 30he BPGA utilises a pool methodology, shown in Fig. 1. 25,31 This differs from a generation-based code, where structures belong to and are evaluated generation by generation.When executed in parallel, multiple instances of the BPGA are started sequentially within a run.Each instance is a separate BPGA run with its own set of processes.The BPGA initially generates a fixed-size pool of n random geometries and places them in a central database file which is available to the other instances of the program.In the present study the pool is set to n = 15 random geometries.These initial geometries are fixed so that no two atoms are overlapping.
In the local minimisation the energy of a structure is minimised with respect to its coordinates.This transforms the PES into a simpler stepped surface, greatly reducing the search space.If an instance becomes free and all structures in the pool are being or have been minimised the instance will continue to evaluate further random structures.If one of these new structures is lower in energy than the highest energy cluster in the pool, the new lower energy structure will replace it.
Once the initial pool of structures has been minimised, offspring and mutants are produced through crossover and mutation.The choice of producing either an offspring or mutant is based on the mutation rate, which is set to anywhere between 0-100% of the fixed pool size.In the present work the mutation rate is set to 10%.
Mutation is defined as the selection of a cluster at random from the pool and the displacement of two of its atoms by up to 1 Å.Other mutations schemes are available in the code, including generating a new random geometry or, for bimetallic systems, swapping unlike atoms.
Selection for crossover is carried out using the tournament method.Once selected, clusters undergo crossover according to the Deaven and Ho cut and splice method. 37The cutting plane is weighted based on the fitness of each of the clusters selected.A higher fitness represents a lower energy.
A local minimisation of the offspring is performed and its energy is checked against those of the other structures in the pool.If the offspring's energy is lower than that of the highest energy structure in the pool, the offspring structure replaces it.Convergence is achieved when the energies of the structures in the pool differ by no less than 10 −3 eV.For larger clusters, convergence may not be achievable.In this case the lowest energy structure in the pool after a run of around 500 separate minimisations is taken as the putative global minimum.
The BPGA also has the ability to perform a DFT-level global optimisation of a cluster supported on a surface.This method, within a generation based code, has been demonstrated previously. 38

Energetics
The binding energies per atom were calculated using where N is the total number of atoms, E Ir N is total energy of an N-atom Ir cluster and E Ir is the energy of a spin-polarised Ir atom.
The stability of the clusters, relative to their N + 1 and N − 1 neighbours, is given by their second-order differences Δ 2 E, calculated using

Results
The BPGA calculations were performed on the UK's national supercomputer ARCHER. 36Each was run in parallel with eight instances of the code operating on the pool.The theoretical scaling of this parallel pool methodology has been shown previously. 31Around 500 structures were evaluated for each cluster size.Due to the high computational cost of the calculations multiple runs are not possible for each system size.
The parallelism within the code, and the scaling capabilities of VASP on ARCHER, allows for spin-polarised calculations to be carried out during the global optimisation.The binding energies E b , point groups and spin multiplicities of the putative global minimum are given in Table 1.The coordinates of the global minima and the additional minima discussed are supplied in the ESI.2][23] Some structures have been previously characterised and give a good indication of the ability of the BPGA to find the putative global minimum at a given level of theory.Other structures are reported here for the first time.
The BPGA successfully finds the C 2v dimer-capped ("house") structure, shown in Fig. 2, as the putative global minimum for Ir 10 and scales successfully well beyond this previous 10-atom limit.The putative global minimum structures of Ir N (N = 11-20) clusters are shown in Fig. 3.
The overall global minimum structure for Ir 11 is a trianglecapped cube.This structure, together with a second highly competitive low lying minimum, an edge-bridged structure based on the Ir 10 "house", are shown in Fig. 4. The two structures differ by 0.05 eV.The global minimum structure is a high spin structure with a spin multiplicity of 4, compared with the competitive minimum's multiplicity of 2. The spin polarised DFT global optimisation has allowed this lower energy putative global minimum to be reported for the first time.
The additional Ir in Ir 12 now makes it possible to complete a third cubic face and the 3 × 2 × 2 D 2h cuboid structure becomes the global minimum.For Ir 13 the extra Ir bridges an edge on one of the cubes.
It was thought that the global minimum structure for Ir 13 may be a C 4v structure, with an Ir atom capping an end of the cuboid, and that the cubic bounding cell used to generate the initial random geometries may have biased the search against any elongated structures.Local minimisations were carried out on C 2v centre edge-bridged, C s top edge-bridged and C 4v topcapped cuboid structures, shown in Fig. 5.The C 2v centre edge-bridged structure locally minimises into a face-capped C s structure.The structures were found to lie 0.26, 0.33 and 0.97 Fig. 2 The dimer-capped ("house") structure of Ir 10 .
eV higher in energy than the putative global minimum structure from the BPGA search, respectively.This structural preference is also confirmed by previous work on Ir 13 . 19,20he structure of Ir 14 is a dimer-capped 3 × 2 × 2 cuboid, analogous to the Ir 10 structure with the dimer lying perpendicular to the long edges of the Ir 12 cuboid structure.This structure has been previously shown. 21Upon addition of an extra Ir, the structure of Ir 15 becomes a trimer-capped cuboid, with the Ir 3 trimer lying parallel to the long sides of the rectangular face.
The preference of Ir 16 is shown to be a slightly deformed L-shaped cubic structure, so that two elongated 3.1 Å bonds can form between two cubes of the structure.Cubic bounding may also have affected Ir 16 , as it follows from the structure of Ir 8 , a cube, and Ir 12 that the global minimum could be a 4 × 1 × 1 cuboid.This structure was assessed alongside two other low-lying minima, T-capped and square-capped cuboid structures.The structures and relative energies of these minima are shown in Fig. 6, with the 4 × 1 × 1 cuboid found to lie 0.88 eV higher in energy than the BPGA global minimum.
The structure of Ir 17 shows the additional Ir in between two cubes of the Ir 16 cuboid, the start of a complete 3 × 3 × 2 cuboid.The additional Ir of Ir 18 sits between two cubes of the L-shaped Ir 17 cluster and forms the complete 3 × 3 × 2 cuboid.
The extra Ir in Ir 19 caps a cubic face on the 3 × 3 × 2 cuboid.This putative global minimum was tested against an edgebridged structure, which was believed to be the more likely global minimum structure following the trend seen in smaller   clusters.The C 1 edge-bridged structure, shown in Fig. 7, was found to have a spin multiplicity of 6 and to be just 0.008 eV lower in energy than the BPGA global minimum.
The energetic differences between the low lying minima of Ir 11 and Ir 19 are far smaller than that seen for the various Ir 13 minima and smaller than the error in the current DFT calculations.To determine accurately the global minimum, higher level calculations, coupled with experiment, will be required in the future.
The structure of Ir 20 is analogous to that of Ir 13 with an Ir For several N-atom clusters E b is higher (less stable) than for the N − 1 cluster: in particular Ir 13 , Ir 15 and Ir 19 .Ir 20 , despite Ir forming a dimer capping a face, has a higher E b than Ir 18 .The increased stability displayed by the even numbered clusters is more clearly shown by their second-order difference plot shown in Fig. 9.Each odd numbered cluster clearly shows decreased stability, indicated by a positive Δ 2 E, relative to its even numbered neighbours.
The spin multiplicities of the clusters are given in Table 1.All global minima are found to have multiplicities of between 2 and 7.This is a clear indication of the importance of magnetic effects in these clusters.In particular, high spin structures are shown to be the global minimum structures of Ir 11 and Ir 19 .It is likely that these and other low energy minima would have been excluded from the searches if non-spinpolarized calculations had been carried out.

Conclusions
The BPGA has successfully coupled the computation resources of ARCHER 36 with the scaling capability of a pool genetic algorithm methodology.This has allowed the direct DFT global optimisation of Ir N (N = 10-20) clusters, with some    global minimum structures being reported here for the first time.
The code has captured the cubic nature of the sub-nanometre Ir clusters with the putative global minima evaluated and compared with previous results.For Ir 11 and Ir 19 their spin has been shown to determine their global minimum structures, which would have otherwise been missed in a lowspin search.The use of spin polarized calculations has been made possible because of the BPGA's ability to utilise greater computational resources.The structural characterisation of a system is a vital first step in exploring its catalytic properties.The structures of the Ir N clusters will form the basis of future studies of the catalytic properties of the system, including modelling their interaction with small molecules.
The cubic structures found here are in agreement with higher level CCSD(T) and CASSCF calculations on Ir 8 , reported by Dixon et al. 18 In future work, we will investigate the effect of changing the functional in DFT calculations on Ir clusters; in particular we will explore the use of meta-GGA functionals. 47,48he development of the BPGA will continue.This will include the implementation of new features, such as the global optimisation of a system in the presence of a ligand or directly on a variety of surfaces.Further improvements to the efficiency of the parallel scheme will also be made.The code will be applied to a variety of new supported and ligated mono-and bimetallic cluster systems.The number of interfaces to common quantum chemistry programs within the BPGA will be expanded and applied to the global optimisation of systems at theory levels beyond DFT.

Fig. 4
Fig. 4 Putative global minimum (left) and competitive low-lying minima (right).Relative energies and spin multiplicities (in brackets) are shown below.

Fig. 5
Fig. 5 Global minimum (top-left), top edge-bridged (top-right), centre edge-bridged (bottom-left) and top-capped (bottom-right) cubic structures assessed for Ir 13 .The centre edge-bridged structure is shown after local minimisation to the face-capped cuboid structure.
2 dimer lying perpendicular to the long edges of the 3 × 3 × 2 cuboid structure.The binding energies of the clusters give an indication of the relative stability of the clusters, with more negative values indicating greater stability.E b values are shown in Fig. 8. Overall, E b tends to decrease as N increases, with Ir 18 having the lowest E b .

Table 1
Binding energies E b , point groups and multiplicities (2S + 1) for the putative global minimum of Ir N (N = 10-20) clusters