Machine Learning Accelerated Genetic Algorithms for Computational Materials Search

A machine learning (ML) model is trained on-the-ﬂy as a computationally inexpensive energy predictor before analyzing how to augment convergence in Genetic Algorithm (GA)-based approaches by using the ML model as a surrogate. This leads to a machine learning accelerated genetic algorithm (MLaGA) combining robust qualities of the GA with rapid learning of the ML. The MLaGA is used to search for stable, compositionally variant nanoparticle alloys to illustrate its capability for accelerated materials discovery, e.g. , nanoalloy catalysts. The MLaGA, in this case, yields a 50-fold reduction in the number of required energy calculations compared to a traditional “brute force” GA. This makes searching through the space of all compositions of a binary alloy particle in a given structure feasible, using density functional theory calculations.


I. INTRODUCTION
The current rate of discovery of clean energy materials remains a key bottleneck in the transition to renewable energy, and computational tools enabling accelerated prediction of the chemical ordering and structure of such materials, e.g., nanoparticle alloys and catalysts, are in high demand.
Genetic algorithms (GAs) are metaheuristic optimization algorithms inspired by Darwinian evolution. Performing crossover, mutation and selection operations, the algorithm progresses a population of evolving candidate solutions. Selecting well designed operators and optimal parameters, GAs have exhibited a high degree of robustness in terms of finding ideal solutions to difficult optimization problems 1,2 . The robustness results from the evolutionary process being able to advance solutions that would have been very difficult to predict a priori. Though, GAs often require a large number of function evaluations, resulting from typical offspring not being very "fit" solutions. Modern machine learning (ML) methods have the capacity to fit complex functions in high-dimensional feature spaces while controlling overfitting 3,4 . However, the high-dimensional feature space means that finding an optimum in an ML model is not a simple task. The robustness of the GA is analyzed while accelerating its convergence through integration with an on-the-fly established Gaussian process (GP) regression model of the feature space.
For materials applications, GAs have typically employed (semi-) empirical potentials (EP) 5-11 to describe the potential energy surface (PES). [12][13][14][15] The utilization of more accurate methods to describe the PES, such as a) Electronic mail: teve@dtu.dk b) Electronic mail: bligaard@stanford.edu density functional theory (DFT) has been limited, due to computational cost. To account for the increased computational cost of searching the PES directly with DFT, studies have often been limited in size, 16 though these methods have successfully been used in a number of investigations. [17][18][19][20][21][22][23][24][25] This study focuses on utilizing the GA to gain an understanding of chemical ordering within larger particles. Searching across a range of compositions is particularly important in the field of materials discovery, where composition can have a profound effect on the desired property e.g. catalytic activity. 26,27 Further, the optimal composition may vary with the size of the nanoparticle. Therefore, the accurate description of chemical ordering is important; where, for certain motifs, the ordering is very complex. 28 Focus is placed on expediting a fast unbiased homotop search by reducing the number of energy evaluations needed to explore the PES and locate the putative global minimum (GM) for a given template structure.

II. RESULTS AND DISCUSSION
The chemical ordering of atoms is optimized for a 147atom Mackay icosahedral template structure. 29 Searches elucidate the full convex hull of possible Pt x Au 147−x for x ∈ [1, 146] compositions. A small number of PtAu compositions will preferentially distort to form rosette-icosahedral instead of the Mackay icosahedral structures. 30 The GA locates these structures in a number of cases, though as structure optimization is not the focus of these benchmarks, when the rosette distortion occurs, the calculations are prevented from entering the population preserving the template structure.
The excess energy is used to assign fitness within the GA, as in Equation 1.
For particles containing a total of N atoms, n A and n B are the number of atom type A and B, respectively. E AB is the total energy of the mixed particle, while E A and E B are the energies of the pure particles. The number of homotops for each particle rises combinatorially toward the 1:1 composition. The number of possible homotops is given by Equation 2.
There are a total of 1.78 × 10 44 homotops for all 146 compositions. The total number homotops for each composition is shown in Fig. 1 as well as an example of a randomly ordered icosahedral structure under consideration in this study. We first run a traditional GA (described in details in the Methods section) to baseline our benchmark and then describe the ML extensions and their results. When using the traditional GA, it is possible to locate the hull of local minima with ∼16,000 energy evaluations. This is significantly lower than the total number of homotops that are present, and thus the number of energy calculations required if a brute-force method was used (1.78 × 10 44 ). However, this is still typically above the number of energy calculations one would wish to perform if a more expensive energy calculator were being employed. To overcome inefficiencies in this method, the underlying search algorithm is optimized and coupled with machine learning selection. A GP regression model is used to predict excess energies of nanoparticles before employing electronic structure calculations. The squared exponential kernel was utilized for the mapping function, as in Equation 3.
The kernel is applied to determine relationships between the fingerprint vectors (x) of two candidates, where x is the Euclidean L 2 -norm and w denotes the kernel width. The training dataset is comprised of unique numerical fingerprint vectors, with features representing distinct chemical ordering within a particle, based on a simple measure, the averaged number of nearest neighbors, as in Equation 4.
Where, e.g. #A-A accounts for the number of homoatomic bonds between atom type A. The summed mass (M ) is appended to account for compositional changes. The ML model is trained on relaxed nanoparticles, though predictions are based on features generated for the unrelaxed structure. The set of descriptors generated in the fingerprint vector are invariant to small changes to the coordinate system, such as a small expansion or contraction of the lattice resulting from the geometry relaxation. A similar ∆-learning method, has been discussed by von Lilienfeld et al. 31 Within the ML accelerated GA (MLaGA) implementation, exists two tiers of energy evaluation, one by the ML functions giving a predicted fitness and the other by the energy calculator providing the actual fitness. A nested GA has been implemented to search the surrogate model representation, generated by the ML. This acts as a high throughput screening function based solely on predicted fitness, running in the "master" GA. The nested surrogate GA takes the current population and is able to progress through additional search iterations, where evaluation and selection are based only on the current model of the PES. The final population from the nested GA returns unrelaxed candidates to the master GA. This is well suited to making large steps on the PES without performing expensive energy evaluations. A difficulty when searching with the MLaGA is that convergence criteria typically used in these studies is no longer suitable. The MLaGA methodology is specifically implemented to limit the number of energy evaluations that are performed. Therefore, every candidate in the generation typically progresses the population. This progression within the population continues until the ML routine is unable to find new candidates that are predicted to be better, essentially stalling the search. For this reason, convergence is considered to have been achieved by the point at which the ML routine prevents new candidates from being evaluated. The general MLaGA methodology is shown in Figure 2.
The GA can be run with a pool or generational population. When running the MLaGA with a generational population, a ML model is trained and utilized to search for a full generation of e.g. 150 candidates. When combining the MLaGA with the generational nested GA, a greater number of candidates are generated in total, compared with the traditional GA. However, the majority of candidates generated in the nested GA routine are discarded prior to the expensive energy evaluation step. Therefore, the MLaGA with a nested search, is able to locate the full convex hull of minima in an average of 1200 candidates. It is possible to reduce the total number of energy calculations through the employ of different acceptance criteria. Tournament acceptance was particularly efficient at reducing the number of required energy calculations, reducing to fewer than 600 for the search.
Tournament acceptance is able to improve search efficiency by restricting the number of candidates passed from the nested, to the master GA. To exploit this further, the MLaGA can also be run with a pool based population where the surrogate model is trained for each new data point resulting from an electronic structure calcu-lation. In this case, the search must progress in serial. Despite the potential for further reduction in the number of calculations required, this may end up being time consuming. This is because, performing the electronic structure calculations cannot be parallelized, as would be possible with the generational population. When this methodology is utilized the number of energy calculations required to search the convex hull is approximately 310.
When training a new model for every energy calculation, it is also possible to exploit uncertainty within the variance distribution on the predicted mean, as in Equation 5. 32 Where a new candidate has the fingerprint vector x * , k is the covariance vector between a new data point and the training data set, K is the covariance, or Gram, matrix for the training data and λ is the regularization hyperparameter. In order to progress the search as efficiently as possible, the cumulative distribution function (cdf), as in Equation 6, is used as the fitness of a candidate.
When the fitness function also accounts for the variance, it is possible to utilize the inherent uncertainty within a prediction to either exploit the current known information in the model, or to explore unknown regions of the search space. 33 The cdf is calculated up to the current known fittest candidate in the composition. The pool based MLaGA is able to locate the convex hull of stable minima in approximately 280 energy calculations. A comparison of the different methods is in Fig. 3. There are clear advantages to performing the search with the augmented ML method.
To ensure that advantages of the methodologies discussed above were not an artifact of utilizing the less accurate EMT calculator, the MLaGA was tested searching directly on the DFT PES. As a significant reduction in the number of energy calculations is likely to be achieved and parallelization of calculations is favorable, the search is performed with the generational population setup. Utilizing the MLaGA methodology, while allowing the nested search to run for a greater number of generations, it is possible to locate the convex hull of minima with approximately 700 DFT calculations. When optimizing geometries with the DFT calculator, there was a 0 eV barrier to structural rearrangement for a small range of the Au deficient compositions.
The convex hull located for the DFT search is in Fig.  4. The shaded region shows the difference in stability between the distorted structures and the most stable icosahedral structures located. The complete core-shell Au 92 Pt 55 structure is located as the most stable for both the EMT and DFT searches. There is good general agreement between the structures obtained elsewhere on the hull, aside from the region of distortion. Further there is broad agreement in the efficiency of the search routines based on the benchmarking and actual searches. Fig 4 also shows the convergence profile as a function of each subsequent DFT calculation. The abrupt bend after around 150 calculations corresponds to a particularly favorable chemical ordering that is distributed to all compositions in the following calculations. This is of course an effect of similar chemical ordering across the whole Au-Pt composition range.
Coupling ML with the GA provides significant advantages in accelerating searches. Performing a search on the surrogate model provides a cheap energy descriptor without requiring expensive electronic structure calculations to assess stability of these nanoparticles. The exact method should be optimized based on the advantages of parallelizing the execution of energy calculations and reducing the total CPU hour cost of the search. A hierarchy of methods have been utilized to reduce the total number of energy calculations required to fully search the convex hull of local minima from 16,000 to around 300.

A. Computational details
The effective-medium theory (EMT) potential 12 is used as the energy calculator for initial benchmarking. The fast inertial relaxation engine 34 optimization routine is utilized to relax the structures, with forces on all individual atoms minimized to at least 0.1 eVÅ -1 . DFT calculations are performed using GPAW with a real space implementation of the projector-augmented wave method. 35 GPAW is run in the linear combination of atomic orbitals mode 36 with a double zeta basis set and RPBE exchange correlation functional. 37 Calculations are run spin-polarised with a Fermi smearing of 0.05 eV in a non-periodic 32 × 32 × 32Å unit cell.

B. Traditional Genetic Algorithm
The GA implemented within the Atomic Simulations Environment (ASE) software package 38 has been utilized. A niching fitness function is employed to efficiently search across the full compositional convex hull. 39 When initializing the traditional GA, the population size is set to 150 candidates. The method for selecting parents is handled by roulette wheel selection. Selection probabilities are directly related to the ascribed fitness, which accounts for the stabilities of the nanoparticle. Offspring are created by either mating two parents, or by mutating a single candidate. The mating and mutation routines are mutually exclusive and thus are not allowed to stack i.e. performing crossover and mutation before evaluation. Cut and splice crossover functions, described by Deaven and Ho, 5 are used to generate new candidates with a call probability of 0.6. Random permutation mutations are utilized with a call probability of 0.2, e.g. swapping the positions of two random atoms of different elemental species. A random swap mutation is also employed with a call probability of 0.2, where one atom type is swapped for another. Convergence criteria is assigned through a lack of progression in the population e.g. the fitness of the population does not change for a number of generations. 6 The GA is run with relatively loose convergence criteria, when there has no observed change in the population for two generations, the search is concluded.

IV. DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.