Pure and Hybrid Evolutionary Computing in Global Optimization of Chemical Structures: from Atoms and Molecules to Clusters and Crystals

2015·Arxiv

Abstract

Abstract

The growth of evolutionary computing (EC) methods in the exploration of complex potential energy landscapes of atomic and molecular clusters, as well as crystals over the last decade or so is reviewed. The trend of growth indicates that pure as well as hybrid evolutionary computing techniques in conjunction of DFT has been emerging as a powerful tool, although work on molecular clusters has been rather limited so far. Some attempts to solve the atomic/molecular Schrodinger Equation (SE) directly by genetic algorithms (GA) are available in literature. At the Born-Oppenheimer level of approximation GA-density methods appear to be a viable tool which could be more extensively explored in the coming years, specially in the context of designing molecules and materials with targeted properties.

Kanchan Sarkar Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, Minnesota, USA e-mail: pcksarkar@gmail.com&ksarkar@umn.edu

S. P. Bhattacharyya Department of Chemistry, IIT Bombay, Powai: 400076, Mumbai, INDIA. e-mail: pcspb@iacs.

1 Introduction

An N-atom system can be arranged in many different ways to form chemically meaningful structures. These structures represent minima on the 3N dimensional potential energy surface (PES) on which the nucleii move. One of these minima is the global one and search for the global minimum is a computationally challenging job. Optimization strategies play a very important role in theoretically elucidating the minimum energy structures supported on a potential energy surface. The ideas of PES and minimum energy structures are rooted in the Born Oppenheimer (BO) theory of molecules. BO theory starts by separating the slow nuclear motion from the much faster electronic motion which allows the electronic charge distribution in a molecule to adjust itself practically instantaneously as the nucleii move. The separation results in the creation of the so called adiabatic potential energy surface En(x1,x2x3N) where Ens are the eigenvalues of a molecular electronic Hamiltonian H(r,x1,x2x3N) in which the nuclear displacement coordinates (x1,x2x3N) appear as parameters. The PES can have a number of stationary points. At all such points

If in addition, the eigenvalues of the matrix of second derivatives of En with respect to 3N nuclear displacement coordinates (xi) are > 0, the stationary point is identified with a minimum energy structure. Out of 3N eigenvalues 3 representing transnational and 3 representing rotational motion of the center of mass are discarded as motion along these coordinates do not define any new chemical structures. It is immediately seen that the problem of locating minimum energy structures on an N-atom PES is a problem of minimizing a 3N dimensional function En(X1,X2X3N). The function itself may be provided on the fly by solutions of the appropriate molecular Schrodinger equation or may be generated by superposing calibrated pair or higher body potentials.

In addition to the minima on the PES, chemists are also interested in locating first order saddle points and constructing minimum energy path (the reaction path) that connects the saddle point to two neighboring minima. At the first order saddle

and all but one of the eigenvalues (excluding the six) of the 2nd derivative matrix are positive, the special one being negative.

The present review focuses attention on the ways of exploring PES of molecules, clusters, crystals, etc. by invoking techniques of evolutionary computing for locating minima. In addition, we also review attempts to solve the atomic or molecular electronic Schrodinger equation (MESE) by similar techniques. In the case of molecular Schrodinger equation, the objective may be either to construct the PES by solving the MESE at a series of configuration and locating minima on the surface or to look for the deepest minimum in one go.

2 The Optimization Problem

Any problem can in principle, be reduced to an equivalent optimization problem. The domain of the function to be optimized is called the search space (S). Goodness of any solution in S is measured with a mathematical function, called the objective function or fitness function (O : S R). A fitness landscape (mapping from a configuration space into the real numbers) may be considered as a triple (S,O,D), where D is a metric defined on S.

In general, optimization problems involve setting a vector X of free parameters of a system in order to optimize (maximize or minimize) some objective function O(X) subject to the satisfaction of inequality constraints gk(X), equality constraints hk(X), as well as upper and lower bounds on the range of allowable parameter values. Any constrained nonlinear optimization problem that deals with the search for a minimum of a nonlinear function O(X) of m variables can be formulated as follows:

Global optimization methods differ from local optimization methods in that they attempt to find not just any local optimum, but the smallest (largest) local optimum in the search space S. Global optimization is a difficult problem since no general criterion exists for determining whether the global optimum has been reached.

3 Types of Optimization Algorithms

The goal of any optimization method is to design mathematical and computational infrastructure in order to locate the extremum (minimum or maximum) of a function (or functional). Automated optimization algorithm designing (i.e. the algorithm can learn adaptively to tune all its parameters by its own) can be very costly and not always straightforward. Expertise is needed to tackle with the parameter sensitivity, efficiency, robustness and solvability of the algorithms. Even mathematically equivalent formulations often differ substantially in efficiency and solvability. So it requires careful thoughts along with mathematical details while designing these optimization methods. Based on the use of random numbers, optimizers can be divided into two major categories, deterministic and stochastic. Deterministic algorithms are in general unidirectional (there exists at most one way to proceed, otherwise, the algorithm gets terminated) and do not use random numbers in any step of execution. State Space Search, Branch and Bound, Gradient based methods (such as Newton’s method, Gauss-Newton method, Steepest-Descent method, Levenberg-Marquardt algorithms), Conjugate – Direction methods (such as Conjugate-Gradient method, Fletcher-Reeves method, Powell’s method, Partan method), Quasi-Newton methods (such as Davidon-Fletcher-Powell method, Broyden-Fletcher-Goldfarb-Shanno method, Hoshino method) are some examples of deterministic optimization methods [1]. The gradient-based methods are the most used classical techniques for solving optimization problems. However, these methods can only be applied to the objective functions that are continuous and differentiable. Some of them, such as the Newton method, even require the knowledge of the second derivatives of objective functions as well. These techniques are also limited to the optimization problems having a sole optimum, as they are gradient-based. However, most of the objective functions in real-world problems are not differentiable and they have a large number of local optima, which makes these methods inapplicable. The unavailability of analytical gradients can be tackled by using numerical gradients. But this is only feasible for relatively low dimensional problems. So gradient-free optimization procedures are often much sought after techniques in fairly large scale optimization.

Stochastic algorithms, on the other hand, incorporate the concept of probability, employs at least one instruction or at least one operation that makes use of random numbers and do not use the gradient or Hessian matrix. The function to be optimized need not even be continuous or differentiable. If the relation between a candidate solution and its “fitness” are not so obvious or too complicated, or the dimensionality of the search space is very high, it becomes harder to solve a problem deterministically. Apart from these reasons there can be many other pitfalls and booby traps [2, 3, 4, 5, 7, 8, 6] (such as conflicting objectives, heavily constrained fitness function, oversimplification, overfitting, non-separability, scalability, rugged and deceptive fitness landscape, neutrality, epistasis etc.) that make the optimization problems difficult to handle and can lead an optimizer to a suboptimal region in the search space. Simulated annealing, Monte-Carlo sampling, Stochastic tunneling, Parallel tempering, Stochastic Hill Climbing, PSO, GA, Evolution Strategies, Memetic Algorithms, Differential Evolution are some examples of stochastic algorithms [9, 10, 11, 12, 13, 2, 3, 4, 5, 14, 6, 7, 8, 15]. Our focus in this review has been mainly on pure evolutionary computing techniques and their hybrids in which power of several methods are combined to produce a superior algorithm for global search. Averaged over the set of all problems, any pair of optimizers perform equivalently, which essentially implies that designing of a general purpose optimizer is bootless. Independent of the fitness function, one cannot (without prior domain knowledge) successfully choose between two algorithms based on their success in a different set of problems.

4 Evolutionary Computation

Biological life appeared on the earth probably as a culmination of many chance events involving chemical and physical interactions of molecules. Ever since the appearance of unicellular life, progressively complex life-forms have gradually evolved by the process of genetic evolution. The process of genetic evolution has elements of adaptive learning built into it. Evolutionary Computing methodologies are mathematical models of natural evolution implemented on a computer and are the most versatile complex problem-solvers. These techniques operate on a set of individuals or chromosomes (population) simultaneously. Each individual represents a potential solution to the problem being solved. The cardinality of population depends on the complexity of the problem. A chromosome is a sequence of genes (essentially the system parameters). Goodness of an individual is defined with a mathematical function, called fitness function (f : S R). Individuals having lower fitness value are slowly washed out by the dominant competitors. There is thus a natural process of screening in the course of genetic evolution which finds expression in the Darwinian dictat of the survival of the fittest [16]. Chromosomal crossover and mutation produce new features in the chromosome. The evolution is a continuous process by which a species continuously strives to attain a genetic structure of the chromosomes that maximizes their probability of survival in a given environment. The evolutionary algorithms continue their iterations to improve the fitness of the individuals until an fitness maximal solution is found. The way individuals in the population are distributed has a major influence on the search. In the panmictic model, there is no structure in the population, all individuals are potential partners. Each individual can interact with every other one in the population during the evolutionary process. However, it is possible to define some structure providing the algorithm with higher exploration capabilities of the search space with respect to panmictic populations. There are two main canonical kinds of structured populations, namely distributed or coarse-grained and cellular or fine-grained [17, 18]. The choice must be made carefully.

Basic Ingredients of Evolutionary Computation

In a diploid organism like the human being a pair of parental chromosomes form the chromosome-pair of the child, each chromosome containing millions of genes. One gene from the father and the corresponding gene from the mother constitute a gene-pair for the child, the gene-pair (genotype) determines the specific attributes or

Fig. 1 Flowchart for simple evolutionary computation

characteristics called the phenotypes. Most gene values (alleles) are inherited by the child from the parents unaltered; but on rare occasions, one or more of them may undergo change, by a process called mutation. If such changes are beneficial for the individual for better adaptation to the environment, they have higher probabilities of survival, and have correspondingly higher probabilities for producing their offsprings. Over generations those beneficial changes tend to stay on while the non-beneficial ones tend to disappear. In Evolutionary Computation, new individuals are generated in the search space by applying certain genetic operators to the current population. The dominant individuals will mate more often creating descendants with similar or better fitness in accordance with the natural selection process and statistically moving toward more optimal places in the search space. A simple picture of process steps in sequence for general Evolutionary Algorithms is shown in Figure 1. The three basic components of EC are A) Selection, B) Crossover, C) Mutation. Each component can be realized in a number of different ways. Let us start by looking at different ways of implementing a selection process on an evolving population of potential solutions

4.1 Selection

Selection operator selects a proportion of the existing population with a specified probability to create the basis for the next generation. The operator is designed to ensure that promising solutions of the population have a greater probability of being selected for mating. The selection process is controlled by selection pressure, which is defined as the degree to which the better individuals are favored. So, selection process is essentially biased towards the more fit individuals and drives the Evolutionary Computation (EC) to improve the population fitness over the successive generations. Higher selection pressure can result in higher convergence rate, but with perhaps a higher chance of premature convergence. On the contrary, EC with low selection pressure will take unnecessarily longer time to find the optimal solution. In addition to the selection pressure, selection schemes should also try to preserve population diversity, to avoid premature convergence. The selection mechanism consists of two steps, the selection probability calculation and the sampling algorithm.

4.1.1 Generational replacement:

Entire set of parents are replaced by their descendants or the n worst parents are replaced with n best offspring.

4.1.2 Truncation Selection:

A proportion (p) of the population are selected based on their fitness value and reproduced 1p times so that the population size is maintained. Less fit individuals are not given the opportunity to evolve.

4.1.3 Roulette Wheel Sampling Algorithm:

Individuals are selected according to their fitness values [9, 19]. The higher the fitness, the higher the probability of being selected. For each individual Ci in a population P the selection probability (ps(Ci)) is given by

The population is then mapped onto a roulette wheel, where each chromosome Ci is given a slot that proportionally corresponds to ps(Ci). To select an individual, a random number is generated and the individual whose slot spans the random number is selected. The wheel is then spun N times, N being the cardinality of the population) to choose individuals until next generation mating pool is fully populated.

4.1.4 Stochastic Universal Sampling:

Population is represented by a pie chart, in which each zone is allocated for an individual of the population. The areas of the zones are directly proportional to the fitness of the representative individuals, exactly as in roulette-wheel selection. But unlike roulette-wheel selection, here a single spin of the wheel is required to construct the mating pool. N (cardinality of the population) number of equally spaced pointers are put around the pie chart. The position of first pointer is generated randomly in the interval (0,1/N). Individuals for mating pool are then selected by generating the N pointers, starting with randomly generated first pointer, spaced by 1/N, and selecting the individuals whose fitness spans the positions of the pointers. Stochastic Universal Sampling provides low spreading over the desired distribution of individuals and is bias-free. This approach is better than roulette wheel, because it keeps the diversity and prevents best individual from dominating the population.

4.1.5 Rank-based Selection:

Individuals are ranked according to the fitness values in ascending order (least fit individual has rank=1, the fittest individual rank=N). Rank based selection uses a function to map the indices of individuals in the sorted list to their selection probabilities. The performance of the selection scheme depends greatly on this mapping function.

For linear rank-based selection, the bias is controlled by adjusting the selection pressure (SP), such that 2SP 0, the expected sampling rate of the best individual being SP, while that of the worst individual is (2SP). The selection pressure for all other individuals in the population can be obtained by linear interpolation of the selection pressure according to the rank. The fitness value for ith individual is calculated as:

where F(i) and ri are the fitness and rank of the ith individual, respectively. Nonlinear ranking permits higher selection pressures (N SP 0) than what is used in the linear ranking method. The corresponding fitness values can be com-

puted by using

where X is the root of the polynomial equation

Rank-based selection schemes help prevent premature convergence, but are computationally expensive because of the need to sort populations, and slow convergence.

4.1.6 Boltzmann Selection:

A temperature like selection parameter controls this selection procedure. Initially the temperature has been kept high and then gradually lowered (selection pressure increases gradually) throughout the EA run. The probability of accepting an individual k:

where Fk and Fa are the fitness of kth individual and average fitness of the population respectively. KB is Boltzmann constant and T is the temperature.

4.1.7 Tournament Selection:

A tournament is held among n randomly picked competitors from the population. The individual with the best fitness value of those random n tournament competitors (winner of the tournament) is then copied into the mating pool. A tournament size of n = 2 is commonly used in practice. Selection pressure can be controlled by controlling the tournament size n. Tournament selection is superior to the proportionate selection as the latter is found to be significantly slower in terms of convergence time. Linear ranking and stochastic binary tournament selection have identical performance in expectation, but binary tournament selection is preferred because of its more efficient time complexity [20].

4.1.8 Sexual Selection:

Crossover takes place only between chromosomes of opposite sex. The sex of the individuals in the current generation are determined either randomly or based on some specific criterion of the individuals. Generally all the females will get to reproduce only once regardless of their fitness level to facilitate exploration of the search space. The male selection is fitness biased. In feminine selection, the female chooses one of the males to mate based on a problem dependent attraction function [21, 22].

4.1.9 Clonal Selection:

A proportion of the fitter individuals are selected for cloning. Each of these individuals receive a number of copies proportional to its position in the ranking. The clones then undergo the maturation process. A given individual and its maturated clones forms a subpopulation and the best of each subpopulation is allowed to pass to the next generation [23].

4.2 Crossover

It is designed for sharing information between individuals by swapping or intermingling the genetic materials of two randomly chosen parent chromosomes, with the possibility that good chromosomes may generate promising descendants. A user de-

Fig. 2 Simple Crossover involves exchanging genetic information, beyond a randomly chosen cross site, by swapping the gene values within the parent’s chromosome. Traditionally chromosomes are crossed at a single point. However, some problems could benefit from using multiple crossover points.

fined crossover probability (pc) determines the crossover frequency, i.e., how often crossover will take place. A general crossover scheme is shown in Figure 2. In the discussion of different crossover techniques that follows, we have represented the jth parent and offspring, of cardinality d, in a population as:

4.2.1 Simple crossover:

A single crossover point i is randomly chosen (say, i), where the two mating chromosomes are cut and the sections after the cuts are swapped to create two new individuals. N point crossover is a generalization of this simple crossover applied to N different segments.

4.2.2 Cut and splice crossover:

This operator is an extension of simple crossover and designed specifically for cluster geometry optimization problems (see section Metal Clusters & Nanoalloys).

4.2.3 Uniform crossover:

It is assumed to be more disruptive in nature than simple crossover. So a lower crossover rate (say, 0.50) is used. Descendants are created by swapping the genes of two parent chromosomes with a predefined probability [24].

4.2.4 Cell crossover:

Cells (a set of genes) of two random parent individuals are swapped with predefined probability to produce two descendants [25].

4.2.5 Heuristic crossover:

A single offspring is generated by linear extrapolation of two parent individual using equation 14 that may lie outside of the range of the two parent vectors [26]. It is also possible for this operator to generate an unfeasible offspring vector. In such a case another offspring created with different choice of r, r being a uniform random variable belonging to the interval [0,1].

4.2.6 Simplex crossover:

Simplex is constructed with randomly selected N parents from the population. Centroid (C) of the simplex is calculated by

The simplex is then expanded by a small degree

and within this expanded simplex, one or more new individuals are sampled [27].

randomly chosen crossover point(s)

Fig. 3 Arithmatic Crossover involves intermingling genetic information, beyond a randomly cho- sen cross site and creates descendants that are the weighted arithmetic mean of two parents. For simplicity all the gene values in the parent strings are shown by same color.

4.2.7 Linear crossover:

Three offspring are generated (one in the exploitation zone and two in the exploration zone), with the two most promising solutions being chosen as descendants by means of fitness evaluation [26].

4.2.8 Arithmetic crossover:

Like simple crossover, chromosomes are crossed at a single or multiple crossover points. All positions arithmetic crossover are also useful. It was introduced to prevent invalid crossover and to increase searching space and convergence speed. Two complimentary linear combinations of the parents are generated on the basis of the arithmetic mean [28].

Arithmetic crossover operates in a way that every gene in descendants is a convex combination of the cooresponding genes from the two parents. is constant in uniform arithmetical crossover or it may vary with generations in non-uniform arithmetical crossover. Figure 3 schematically shows arithmatic crossover.

4.2.9 BLX-α crossover:

To expand the range of arithmetic crossover Eshelman and Schaffer suggested blend crossover [29], which generates a single offspring by blending two floating point parent vectors as follows:

is a function generating a uniform random number in the range (a,b). The user-defined parameter is usually set to a value of 0.5. Setting 0 and 0.25 will give flat crossover [30] and extended intermediate crossover [31], respectively.

4.2.10 Similarity crossover:

In arithmetic crossover, there has always been a problem in selecting a perfect weighting factor. The weight factor () in arithmetic crossover (equation 18) is replaced by similarity measurementbetween parents [32]

In addition to these, there are other variants of crossover schemes such as parentcentric BLX-crossover [33], Fuzzy recombination [34], SBX [35], PCX [36], XLM [37], Laplace crossover [38], differential evolution [39], partition [40, 41], linear BGA [42], UNDX [43], fuzzy connectives based [44], direction-based [45], multiple crossover [46], Ring crossover [47], hybrid crossovers [48] are available in the literature.

4.3 Mutation

The mutation operator perturbs one or more components (genes) of a selected chromosome, regulated by a predefined mutation probability, pm, so as to increase entropy or to decrease the mutual information in the population. Mutation simply restores lost information or import unexplored genetic material into the population in order to distribute solutions widely across the search space and thus avoid premature convergence.

Fig. 4 Mutation causes small alterations at randomly chosen gene(s) (with a predefined mutation probability) in an individual within available range of possible values.

The performance of any evolutionary search algorithm heavily depends on the choice of its main control parameters: population size, mutation rate, and recombination rate and among these mutation rate is the most sensitive control parameter. Mutation-only-EA is possible, but, crossover-only-EA would not work with finite population size. There should be always a proper balance between exploration (discovering promising areas in the search space) and exploitation (optimizing within a promising area already found) abilities of a search algorithm. Crossover exploits existing genetic material, while mutation explores for newer material. Mutation does its best to avoid premature convergence and explore newer areas of search space. High mutation rate increases the probability of searching over large areas of search space; however, it prevents the population from converging to an optimum solution. On the other hand, too small a mutation rate may result in premature convergence. The best value of mutation rate is strongly problem specific and depends on the nature and implementation of the algorithm. There is no universal best mutation rate for most of the real world problems. The crossover and mutation rates (probability) are critical to the success of any EA, and are usually determined by means of trial-and-error. Automating the parameter settings i.e., making the search algorithm self-adaptive has been one of the most fashionable area of current interest among the EA community. A general mutation scheme is shown pictorially in Figure 4. We shall use specific notations for an individual before and after mutation, right through the discussions of different mutation schemes. The chromosome before and after mutation are expressed as follows:

4.3.1 Uniform mutation:

Genes in an individual are selected with a predefined probability (pm) and replaced with a uniform random value according to [28].

is a function generating a uniform random number in the range (a,b), LBi and UBi represents lower bound and upper bound to the ith gene. Uniform mutation tends to forestall premature convergence. In the neighborhood of the optimal region, however, uniform mutation rate may disrupt the fine-tuning of the search to the optimal point by inducing random oscillations.

4.3.2 Nonuniform mutation:

Due to its dynamical nature, non-uniform mutation [28] is one of the most commonly used mutation scheme in real coded EC. Gene(s), selected with a predefined mutation probability (pm), are mutated as

where r is a random number, and t is the generation number. The user defined parameter b takes care of the non-uniform nature of the mutation step sizes.

4.3.3 Gaussian mutation:

A noise with a mean of 0 and a standard deviation of 0.1 times the maximum value of the gene is added to the gene-value. A new individual is obtained by adding a

random value to each element of the selected parent.

The mutation function that describes the transformation of an element x into y [49]

is the scale parameter. Luteevy-type mutation [50, 51], Real Number Creep [52] are some variants of Gaussian mutation. The scheme of Real Number Creep:

Based on the PBX-recombination operator and Gaussian distribution function, Dorronsoro et al. designed GPBX[53] mutation operator.

4.3.4 Dynamic Random Mutation:

The scheme enables the mutation size to be adjusted dynamically [54].

where is a random perturbation vector within the interval of (0,1]. The mutation step size is dynamically tuned by

The parameter b > 0 controls the decay rate of s and tmax denotes the maximum number of generations over which the search takes place. The parameter b governs the shape of the allowable mutation region. The mutation range dynamically decreases as the number of generation increases. The idea is adopted from the annealing procedure in metallurgy.

4.3.5 Directed Random Mutation:

The value of a randomly selected gene is modified according to

where r is a random number in the interval 0 < r < 1. The + or sign is chosen with a probability of 0.5. is the mutation intensity, dynamically adjusted on the basis of the degree of acceptability of mutation over a number of past generations [55, 56].

where rg is a Gaussian random number in the range (0,1) and Naccept is the number of accepted generation over past N generations. Nlower and Nupper are the user defined lower and upper bound to Naccept.

4.3.6 Wavelet mutation:

A continuous-time function is called a “mother wavelet” if it satisfies the following properties:

One example of such mother wavelet is Morlet wavelet:

In order to control the magnitude and the position of is defined as

where a is the dilation parameter and b is the translation parameter. Ling et al. [57], based on the wavelet theory, proposed wavelet mutation where each gene in a chromosome have a chance to mutate according to a preset mutation probabilty pm . The mutated gene is represented by

For Morlet wavelet is randomly generated.

4.3.7 BGA mutation:

The mutated gene yi is computed according to [31]

The + or sign is chosen with a probability of 0.5, rangei is normally set to 0(UBi LBi) and is computed as follows.

is randomly generated and takes the value 1 with a probability p. The scheme ensures that at least one variable will be mutated in an individual. Discrete and continuous modal mutations [58] are generalization of BGA mutation. They only differ in how is computed.

4.4 Supplement to Selection: Elitism

Sometimes crossover and mutation operations destroy potentially valuable genetic information which has evolved during the search. There is no guarantee whether EA will re-discover these lost information or not. If EA rediscovers those previously discarded information, it does so at the cost of some unnecessary time and computational power. To prevent this, the concept of elitism is introduced. Elitism ensure that a small proportion (typically 1 or 2) of the fittest individuals are passed onto the new generation so that the valuable information remains intact within the population and thus reduces the genetic drift. Higher proportion (causes rise in selection pressure) may lead to premature convergence. So the degree of elitism should be adjusted carefully. Elitism can speed up the performance of the GA significantly.

Apart from these operators, initial population also has major influence on the convergence of the algorithms. It is difficult to hit the global optimum in large systems starting from a fully randomized initial population as it leads to nearly identical glassy structures that have similar (high) energies and low degree of order [59]. Often random but symmetric structures may be helpful as initial guess for systems with large dimensions [60]. Presence of constraints often make an optimization problem difficult to solve. Several techniques such as adding penalty term to the objective function (the penalty increases with constraint violation), adjustment of the crossover or mutation operation, inducing instabilities to the infeasible solutions of the populations without modifying the objective function, have been proposed [61].

4.5 Evolutionary Computing Methodologies

Let us consider several variants of EC methodologies that have been in use.

4.5.1 Random Mutation Hill Climbing:

It is the simplest algorithm in the class of Evolutionary Computing methodologies. The algorithm works with one randomly generated string containing all the system parameters, (say, Cartesian position variables for geometry optimization) required for optimization and it involves mutation, as the only evolutionary process in order to climbing up the fitness landscape. Thus the underlying principle is to randomly choose a solution in the neighborhood of the current solution by means of mutating the string and retain this new solution if it improves the fitness function. A modified Random Mutation Hill Climbing was proposed and implemented by Sarkar et al. [56, 62, 63, 64]. The modification lies in additional built-in features for enforcing adaptive control on all the parameters of the search heuristic.

4.5.2 Genetic Algorithms:

GAs are a class of stochastic heuristic search methods based on simplifications of evolutionary processes observed in Nature. They maintain a finite population size and are typically good at both the exploration and exploitation of the search space. They are invoked particularly in large, complex, and poorly understood problem domain. The search begins with randomly initializing a population of individuals followed by selection operation to create a mating pool. The mating process, implemented by crossing over the genetic material from two parents to form the genetic material for two new solutions. Random mutation(s) is(are) applied to promote diversity. If the new solutions are better than those in the population, the individuals in the population are replaced by the new solutions. The key idea is to maintain a population of chromosomes, that evolves over time through a process of competition and controlled variation. GA have had a great measure of success in search and optimization problems due to their ability to exploit the information accumulated about an initially unknown search space in order to bias subsequent searches into useful subspaces. The population scalability of performance in evolutionary computation has been an issue pursued since 1990s. The population sizing equation of Goldberg [65], M¨uhlenbein’s study on minimal population size that can ensure convergence to the minimum [66], the work relating to adaptive control of population size by Arabas et al [67] are the early examples of interest in this area. However, GA is trivially parallelized and its parallel efficiency could often outweigh the higher computational demand in every generation [68].

4.5.3 Evolutionary Programming:

Both Evolutionary Programming (EP) and Evolution Strategy (ES) are known as phenotypic algorithms, whereas the GA are genotypic algorithms. Phenotypic algorithms operate directly on the parameters of the system itself, whereas genotypic algorithms operate on chromosomes or individuals representing the system. In recent work however the distinction has become blurred as ideas from the ES and EP have been applied also to strings or chromosomes that constitutes a population in GA. The primary difference between EP and the other approaches (GA, GP, and ES) is that no exchange of material between individuals in the population is made. Evolution is caused solely by appropriate mutation operators.

4.5.4 Evolution Strategy:

The earliest form of ES is based on populations that consist of two competing individuals and the evolution process consists of the two operations mutation and selection. Only the fittest of the two individuals is allowed to produce descendants in the following generation. Later, general versions of the algorithm were developed through increasing the number of members in a population (-parents and -descendants), among which the most common strategies are () and (). In () the parents are simply replaced with the descendants and in (), the parents compete with their children.

4.5.5 Differential Evolution:

DE is a population based technique [69, 70, 71]. The members of the population are chosen randomly to start with (X1,X2Xn). In a Differential Evolution (DE), mutation plays the dominant role. Mutation is used to produce new individual by adding the vector-difference between two randomly selected members (say, Xk and Xl) of the population to another (say, Xm) producing a trial vector Vi, where,

f being the mutation intensity. A total of n such individuals are generated (V1, V2, , Vn). A crossover operator is then used to produce new candidate Ui from the trial vector Vi by setting the jth component of Vi as follows:

The candidate is evaluated, and the better of the current individual and the candidate solution is selected.

4.5.6 Genetic Programming:

GP is an extension of the GA in which the structures in the population are not fixed-length strings that encode candidate solutions to a problem, but programs expressed as syntax trees. So essentially Genetic Programming is a method to evolve computer programs. The genetic operators allow the syntax of resulting expressions to be preserved.

4.5.7 Memetic Algorithms:

MAs are extensions of EA that repeatedly apply a local-search algorithm to exploit the search space in between generations to reduce the likelihood of the premature convergence.

4.5.8 Particle Swarm Optimization:

PSO was introduced by Kennedy and Eberhart [72] as a gradient free evolutionary computation technique that utilizes the ‘collective intelligence’ of a cooperative swarm of particles. Each particle representing a candidate solution is randomly assigned variables to be optimized and then allowed to “fly” with random velocities (vmin vmax) in the search space. The task is to reach a point in the Ndimensional search space where the fitness has the maximum value and which is hopefully the optimal solution as well. Each particle keeps track of its coordinates in the search space, associated with the previous best fitness value (pbest). The overall best fitness value that the swarm has achieved so far is stored as gbest. At t + 1 generation, velocity (v) and position (x) of ith particle are updated based on its current velocity (vti), the distance from its previous best position (pti), and the distance from the global best position (gti) according to

The updating rule applies to all the components of velocities and positions. C1,C2 are the acceleration coefficients usually set to a value of 2.05, and are uniformly distributed random numbers within the interval of (0,1). The inertia weight decreases linearly from 9 to 4 over the whole run to control the momentum of the particles. The coefficient C1 is called particle’s self-confidence due to its a contribution towards the self-exploration of a particle, whereas C2 is called swarm confidence due to its contribution towards motion of the particles in global direction, which takes into the collective wisdom [73]. The velocities are assigned cutoffs Vmin and Vmax to prevent the swarm from disintegrating. The coordinates evolve continuously, and have no cutoffs.

5 Monte Carlo Algorithms

The main principles of Monte Carlo algorithm, named after a famous casino city in Monaco, are ergodicity and repeated random sampling to find the optimal solution. It was originally practiced under more generic names such as statistical sampling. In the 1940s, physicists working on nuclear weapons projects at the Los Alamos National Laboratory coined the present name. The random search, direct Monte Carlo sampling, the random hill climbing methods, Simulated Annealing, Quantum Annealing and Parallel Tempering are some examples that belongs to the class of Monte Carlo algorithms.

Simulated Annealing:

The term “annealing” is commonly used in metallurgy to define a process in which the metal is heated to a high temperature, followed by slow and controlled cooling. If molten metal is quickly cooled, it does not reach the perfect crystalline state of minimum free energy. The essence of the process lies in slow cooling which allows, time for redistribution of the atoms as they slowly lose mobility. Metropolis et al. [74] proposed an algorithm for finding the equilibrium configuration of a collection of interacting atoms at a given temperature. The transition to a new state ( j) from a state (i) is accepted with the probability

Kirkpatrick et al. [13] introduced the concept of Simulated annealing (SA) in combinatorial optimization following Metrpolis’s idea. In SA, energy is regarded as the cost function and a set of parameters {xi} is used to define the configuration at some effective temperature. The target is to find the global minimum of the cost function using temperature as the control parameter. The SA first allows the system to ‘melt’ at a high temperature and the temperature is lowered in small steps allowing it to spend sufficient time at each temperature until the system freezes. At each step of temperature, the sampling must be carried out for a sufficiently long time for the system to reach a steady state. The sequence of temperatures and the number of samplings allowed to reach ‘equilibrium’ at each temperature constitute an “annealing schedule”. At each reconfiguring step perturbations are applied to generate a new point in the search space. The new point is accepted based on the Metropolis algorithm mentioned above (equation 44). Thus the algorithm differs from the HillClimbing method in the criterion used to replace the existing solution with a new one. The advantage of SA is the built-in local optimization scheme, since in the limit of the vanishing control parameter (T), the algorithm is a standard downhill simplex method. A performance comparison with EAs shows different strengths and weaknesses. While the short-range exploitation of SA based algorithms is generally considered to be better than that of pure EAs, EAs generally provide a better long-range exploration. Employing hybrid EAs with enabled local optimization steps in general makes the EAs at least on par with SA-based methods.

An inverse approach (Reverse Monte Carlo) to the classical Monte Carlo algorithms successfully reproduces many salient features of glass and metallic glass systems obtained from molecular dynamics simulation [76, 75]. Along with the success in different optimization problems SA has also been effectively used in quantum mechanical calculations, especially in the optimization of parameters of wave functions, constructing reaction paths for transformations on model potential energy surfaces and for rearrangement in Lennard Jones clusters [77, 78]. We have so far summarized the important aspects of evolutionary computing along with some relevant natural computing methods. We hope the readers will be able to exploit the information to form a well founded basis to choose appropriate schemes which maintain a good balance between exploration and exploitation for specific problems. In the next section we will provide an overview of applications of evolutionary computing and its hybrid methods in the general context of computing electronic structure of atoms, molecules and clusters.

6 Evolutionary Computing in Electronic Structure Calculation

The applications of EC to problems of calculation of electronic structure of atoms, molecules, clusters and crystals have been interesting and numerous. We will review the applications under these categories:

• Solution of Schrodinger equation (SE) by EC. • Locating minima on the complex potential energy surfaces of molecules, clusters, crystal structure prediction.

6.1 GA and Electronic Structure Calculation

The first attempt to marshal the power of GA to tackle the problem of electronic structure calculation can be traced to a paper by Zeiri et al. [79] on the application of GA to the calculation of bound state within the framework of a local density approximation almost 30 years after GA came into being. GA was invoked for directly solving the time-independent SE [80] by Chaudhury and Bhattacharyya in 1998. They tested the workability of their algorithm in finding the ground state of a particle moving in screened coulomb potential and that of an oscillator with quartic anharmonicity. They went on to examine the applicability of their method to the problem of solving the inhomogeneous differential equations of the RSPT with special reference to the ground states of two-electron atoms. Nakanishi and Sugawara (2000) proposed a ANN based way of solving the SE in which the wavefunction was represented on a perceptron type NN, and the weights and the biases were optimized by a micro-GA [81].

Saha and Bhattacharyya (2001) proposed a stable and generalizable basis set free strategy of solving the SE by GA [82]. Test calculations on the ground and excited states of H-atom, coupled harmonic oscillators, and symmetric double well, demonstrated the robustness and workability of the algorithm proposed. Saha et al. (2003) parallelized [83] the fitness evaluation steps and the parallel GA method was successfully invoked to solve an exactly solvable coupled oscillator problem and the H-atom problem. The kinetic energy was computed by numerical FFT while the potential energy was obtained by quadrature. Saha and Bhattacharyya (2004) further explored GA in solving the radial SE for the ground state of two electron atoms and ions [84]. The approach is basis set free and the results were better than HF results. The same authors extended the theory to calculate the ground state of Hydrogen molecule ion (H). In one realization, the internuclear distance (R) was allowed to evolve along with amplitudes of the electronic wave function at different points on a 3-d uniform grid. The equilibrium internuclear separation and the energy was obtained in a single run. Sahin et al. (2006) explored the GA based approach for computing the energy levels of hydrogenic impurity in a quantum dot and then went on to further extend the method for the self-consistent electronic structure calculation of many electron quantum dots [85].

GA has been deftly used to compute the molecular electronic structure of doped as well as undoped polythiophene and polyselenophene oligomers [56, 86] within the framework of a modified SSH method. The density based approach in which the GA operators act only on the molecular geometry strings (arrays of the geometrical parameters) that uniquely define the corresponding fixed nuclei Hamiltonian H(R). H(R) in turn, defines a unitary transformation U(R) which transforms a trial density P(Ro) into a new density PRo)U(R) at the new genetically modified geometry. Thus along with the geometry strings, the density P(R) is forced to coevolve, till [H(R),P(R)] = 0. The technique appears to be sufficiently flexible and could find use in ab-initio calculations of molecular electronic structure.

6.2 GA in the Exploration of PES of Clusters

Traditional optimization methods perform well on strongly convex response surfaces, whereas EAs are particularly appropriate for most of the real-world NP-hard problems in which the objective is multimodal and expensive to evaluate, the gradient of the objective function is not analytically definable, the problem domain is large in dimension and/or multiple linear and non-linear constraints have to be applied simultaneously. Traditional hard computing methods are occasionally fused with the evolutionary methodologies to develop computationally intelligent hybrid systems with moderate computational cost.

Atomic or molecular clusters are aggregates of similar or dissimilar particles and generally have physical and chemical properties different from the bulk. Geometry optimization of atomic or molecular clusters belongs to the class of nondeterministic polynomial hard (NP-hard) problems. Exhaustive search throughout the PES in any reasonable amount of time is almost an impossible task. By the early 1990s EAs were applied to these classes of problems. Deaven et al. [87] applied GA to find fullerene structures up to C60 starting from random atomic coordinates. Empirical Brenner potential was used to model the PES. Dugan et al. [88] had proposed Monte Carlo type local optimization between GA steps for moderately sized carbon clusters. Chaudhury et al. [89] have used EA methods on empirical potentials to obtain optimized structures of halide ions in water clusters, which they then coupled with quantum-chemical (AM1) calculations for simulation of vibrational spectra. By suitably defining the objective functions, EC methods can also be used to locate additional features on the PES rather than just low-energy local and global minima. Chaudhury et al. [90, 91] have implemented methods for finding first-order saddle points and reaction paths, on PES of Lennard-Jones (LJ) clusters up to n=30.

Clusters may have properties which are different from those of discrete molecules or bulk matter: for example, some metals (e.g. palladium) which are non-magnetic in the solid state, are magnetic, with relatively high local magnetic moments in both neutral and anionic PdN clusters [92]. The Coulomb explosion of large raregas clusters release huge energy (approaches the energy of nuclear processes). At very low temperatures (< 2K for 4He), He clusters may form a superfluid droplet. Metal nanoparticles play a crucial role as heterogeneous catalysts in petrochemical, pharmaceutical and clean energy sectors. Novel and interesting nanostructures, having useful chemical and physical properties could be obtained by controlling the nanoparticle size, composition, surface site preferences, and degree of segregation between the metal constituents [93]. Therefore, computing ground state configu-ration of atomic and molecular clusters is a challenging task from the theoretical perspective. The number of local minima of the clusters rises exponentially with the growth in cluster size. It increases further with composition heterogeneity because structures more complex than their homogeneous counterparts are possible. It is due to the existence of isomers with the same geometry and composition but showing different distributions of the constituent particles that the optimization task becomes notoriously difficult [95, 100, 101, 96, 97, 98, 99, 94, 102]. For a cluster made of nA and nB atoms, with nA +nB = N, one finds that a single geometrical isomer can have N!nA!nB! different homotopes. Various attempts are there in the literature to find the global minimum energy structure of these systems, as for example the dynamic lattice searching method [103], basin-hopping (BH) approach [104, 105, 106, 107], adaptive immune optimization algorithm (AIOA) [108, 109, 110, 111] and evolutionary algorithms [87, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128]. Hybrid approaches combining evolutionary algorithm and a local search procedure, have been increasingly used to handle these problems. We review recent work in the area in what follows.

6.2.1 LJ and Morse Clusters

Ab initio quantum mechanical calculation for finding the global minimum geometry of a cluster of atoms or molecules is often extremely expensive [129, 130]. Computational cost increases with the accuracy of the calculation and size of the system. Generally empirical or semiempirical pair potentials are employed for large systems. LJ and Morse functions are the two most widely used and relatively simple pairwise additive potential models describing interactions between atoms which have been regularly used as benchmark to assess the performance of global search methods for cluster geometry optimization. The LJ and Morse potentials are respectively given by

where ri j stands for the Euclidean distance between atoms i and j, r0 is the equilibrium bond length, is the well depth of the pair-potential, is the separation at which the pair-potential between the atoms goes through zero and is the interaction range scaling parameter. For the prediction of the ground state structure of crystals, the USPEX method [131, 135, 136, 137, 138, 134, 139, 140, 132, 133] has turned out to be extremely powerful and already guided materials scientists in finding interesting and unexpected crystal structures. Zhu et al. [141] enhanced the efficiency of USPEX method by designing additional variation operators and constraints for partially or completely fixed molecules. Goedecker et al. [142, 143, 144] offered improved minima hopping with a softening method and a stronger feedback mechanism to predict structures of homoatomic and binary clusters with LJ interaction as well as structures of silicon and gold clusters described by force fields. Cheng et al. [97] employed funnel hopping algorithm coupled with GA to locate the putative global minima of the LJ clusters and the Morse clusters up to N = 160. Wales et al. showed that the use of approximate symmetry provides a more productive way to explore the configuration space and substantially improves the efficacy of global optimization for the atomic clusters [145, 146, 147, 148]. Froltsov et al. [149] developed the “cut and splice” GA augmented with twinning mutation moves for the structural optimization of both a series of single-funnel LJ clusters up to 70 atoms and the double-funnel LJ38 cluster. Binary Lennard Jones (BLJ) cluster optimization is an even more challenging problem from the point of view of combinatorial complexity. They are also interesting because catalytic properties of such clusters depend on the composition and structure. Comparatively little work has been done on the mixed clusters [150, 151, 152] as the search space is dramatically enlarged with the inclusion of more atom types. Mixed LJ clusters of widely varying compositions offer highly interesting additional perspectives on how variations in preferred structures emerge. Hybrid approaches like combining a global optimizer with a local search procedure have proved to be useful for mixed clusters [153, 154, 155, 156, 157, 158, 159]. 17 new putative global minima for BLJ clusters in the size range of 90-100 particles have been predicted by coupling hidden-force algorithm and non-Markovian parallel Monte Carlo search [155]. Dor et al., proposed a multi-swarm based algorithm called PSO-2S [156] and tested it successfully on a number of benchmark functions. It uses charged particles in a partitioned search space and applies a repulsion heuristic on particles to increase the diversity for solving both unimodal and multimodal problems. Using GA-based global structure optimization framework and new set of fitted parameters for LJ potential from high-end ab initio calculations Hartke et al. optimized strongly mixed binary to quinary rare gas clusters [157, 158]. An average offspring method [142] designed for cluster structure prediction was shown to perform better in systems with compact optimal structures. Two individuals are randomly selected and for each atom of the first cluster, the closest lying atom of the other cluster is identified. The corresponding atom of the child is placed randomly on the connecting line between the two parent atoms. The randomness of this operator is necessary to prevent the algorithm from producing a lot of identical offsprings. Leitao et al. [160] claimed that they applied for the first time an Island Model to the optimization of Morse clusters ranging from 41 to 80 atoms, combined with a hybrid steady-state evolutionary algorithm and a local optimization method. The performance was slightly more robust than that of the sequential approach. Dieterich et al. [158, 161, 162] designed EA program suite OGOLEM for structure optimization of mixed clusters. The OGOLEM framework provides both an MPI fronted for MPP parallelization and a threading fronted for SMP parallelization omitting unnecessary MPI overhead. They demonstrated the possibility to design molecules with targetted properties in the area of photochemistry using the OGOLEM framework [161]. There are prescriptions for the inclusion of Taboo-search features [163, 164] into the evolutionary algorithms which might help in reducing the amount of time spent in local optimizations rediscovering already known minima. Daskin et al. [165] presented Group Leaders optimization algorithm in which the influence of the leaders in social groups is the inspiration for the evolutionary technique. The method is applied to locate the geometric structures of LJ clusters as well as to the quantum circuit design problems. Deep et al. [166] made an attempt to solve LJ problem by incorporating a multi-orbit dynamic neighborhood topology in PSO. In multi-orbit topology, the swarm has heterogeneous connectivity with some subsets of the swarm strongly connected while with the others are relatively isolated. This heterogeneity of connections balances the exploration-exploitation trade-off in the swarm. The dynamic neighborhoods topology helps avoid entrapment in local optima.

6.2.2 Ionic clusters

Oleksy et al. [167] calculated the equilibrium geometries and dissociation energies of electronic ground-state of Heclusters (N = 3-35) via an extended GA method and employing a semiempirical valence-bond model of intracluster interactions. They found that for the cluster sizes of N = 3, 4, 7-9, and 14-35, the positive charge of Hedelocalizes over a trimer ionic core, which is more or less linear and centrosymmetric. For N = 2, 10-13, a dimer ionic core develops, while for N = 5 the positive charge is delocalized over the whole cluster; for N = 6, a tetramer ionic core is formed. Atoms outside the ionic core are neutral and are distributed among several solvation rings. The electronic and geometrical structures of sodium cluster anions [168] [Na, n=20-57] and cationic Naclusters [169] were determined by apply- ing GA – density functional theory (DFT). Structures of Indium oxide nanoclusters and zirconia nanoclusters have been predicted by a method combining a robust evolutionary algorithm with classical interatomic potential and quantum chemical models [170, 171]. Kim et al. [172] used GA to design Nanoporous TiO2 for LowTemperature Processable Photoanodes of Dye-Sensitized Solar Cells. Darwinian and Lamarckian schemes within evolutionary algorithms have been compared in the context of structure prediction of the titania (TiO2) polymorphs. Lamarckism in natural evolution regards the effects of “inheritance of acquired characters” as the motive force of evolution, while the Darwinism claims that evolution is nothing but the cumulative processes of natural selection with random mutation and denies the possibility of inheritance of acquired characters. Although the mainstream of today’s evolutionary theory follows Darwinism, it has been found that the Lamarckian scheme is more successful and efficient at generating the target structures [173]. A combination of Buckingham and LJ potential functions describe the interactions between different atomic species in the system

The parameters Ai j, Bi j, Ci j, and are dependent on the species involved in the interaction. To perform a global geometry optimization of clusters resulting from the microsolvation of alkali metal ions (i.e., Na, K, and Cs) with benzene molecules Marques et al. [174] used modified LJ function for nonelectrostatic contributions. Zhu et al. [175] explored all the possible stoichiometries for Mg-O system at pressures up to 850 GPa. They found that two extraordinary compounds MgO2 and Mg3O2 became thermodynamically stable at 116 GPa and 500 GPa, respectively. They predicted the existence of thermodynamically stable Xe-O compounds [176] at high pressures (XeO, XeO2 and XeO3 became stable at pressures above 83, 102 and 114 GPa, respectively). Woodley et al. [177] predicted a greater stability of tetrahedral and trigonal coordinations compared to the tetragonal one for zinc oxide clusters [(ZnO)n, n=1-32] using an evolutionary algorithm with polarizable shell interatomic potentials. Wang et al. [178] successfully applied PSO for the prediction of elemental (Li, C, Mg, Si), binary (SiO2, SiC, ZnO, TiH2, TiB2, MoB2), and ternary (MgSiO3 and CaCO3) compounds in various chemical-bonding environments (metallic, ionic, and covalent bonding). Li et al. [179] obtained the geometrical structures of (Al2O3)n (n = 1-7) clusters via GA coupled with the DFT. Avaltroni et al. [180] identified peculiar low energy isomers of small clusters (C2Al4 and CB2) using both standard random search and GA procedures. Pal et al. [181, 182] determined the lowest energy structure of ZnS quantum dots of different sizes by using a search based on GA coupled with the density-functional tight-binding method (DFTB) and found a new ring-like configurations of ZnS quantum dots which had higher HOMO-LUMO gaps compared to other ZnS quantum dot structures. Pereira et al. developed an evolutionary algorithm for the global minimum energy search for water clusters up to (H2O)20, benzene clusters up to (C6H6)30 and (C6H6)with n=2-20 [183]. Do et al. [184] identified the energy minima of water, methanol, water + methanol, protonated water, and protonated water + methanol clusters with DFT combined with basin hopping. Addicoat et al. [185] optimized a GA in order to identify the minimum energy structure of arbitrarily hydrogenated and hydroxylated fullerenes, and implemented the method for exhaustive calculations on all possible isomers. They suggested that crossover operator did not have signifi-cant effect on the search efficiency of GA and hence the most efficient version of their GA does not employ crossover, thereby reducing it to an Evolutionary Algorithm (EA). Cai et al. [186] used GA to optimize structures of silicon clusters. Using the GA, Zhang et al. [187] predicted the most stable structures and a number of low-energy metastable structures for Si[001] symmetric tilted grain boundaries with various tilt angles which are found to be in very good agreement with the results of first-principles calculations. Yao et al. [188] investigated the high-pressure structures of solid nitrogen through GA combined with first-principles electronic structure calculations. Hooper et al. [189] studied generic defect association complexes in metal oxide materials at experimentally relevant dopant concentrations by a specialized GA-inspired search procedure for doped metal oxides. The search algorithms had been tested on lanthanide-doped ceria (L = Sm, Gd, Lu) with various dopant concentrations.

6.2.3 Metal Clusters & Nanoalloys

Chen et al. [190] proposed a parallel differential evolution for cluster optimization (PDECO) with triangle mutation and migration operators and applied it to Pt clusters with great efficiency. Rogan et al. [191] implemented a GA based methody on small Pd Clusters. The lowest-lying isomers of the copper cluster, Cu9, have been obtained by combining a GA driven approach with DFT [192]. Assadollahzadeh et al. [193] employed a seeded GA technique using DFT together with a relativistic pseudo-potential to search for global and energetically low-lying minimum energy structures of neutral gold clusters Aun (n = 2-20). GA coupled with a tight-binding potential is employed by Ping et al.[194] to optimize neutral lead clusters Pbn (n = 2-20). Pereira et al. [195] made a study on the effectiveness of different crossover operators in the global optimization of atomic clusters. They came up with a ‘Cut and Splice’ crossover operator which proved to be extremely useful. The first parent string is “Cut” randomly along a random horizontal cutting plane into two complementary parts and the second parent string is “Cut” in such a way that the number of the atoms beneath and above the cutting plane are equal to the number of atoms of the two complementary parts of the first parent. “Splice” joins the head of one string with the tail of the other one. Pereira et al. [195] modified this operator by the way it determines the sub-clusters to be exchanged. Here a random atom from the parent string is selected and placed on the offspring. A random number (M ) is generated (N is the number of atoms in the cluster). M atoms from the first par-

Fig. 5 To make the operation consistent and produce descendants with same number of atoms: the center of mass of the parent structures have to be translated to the origin of the coordinate system and have to rotate the plane of the parent clusters until the generated offspring contains the correct number of atoms.

ent string closer to the first selected atom are chosen and placed on the offspring. The remaining (N-M-1) atoms are taken from another parent string closer to the first selected atom. Those atoms which are already occupied are skipped. Figure 5 schematically depicts ‘Cut and splice crossover’. The energetic ground states of gold clusters, for some magic numbers, with up to 318 atoms were obtained by minima hopping method [196]. A GA–DFT was invoked by Sai et al. [197] who explored the size evolution of structural and electronic properties of neutral gallium clusters of 20-40 atoms. The atomic structures and electronic states of Gan clusters significantly differ from the solid but resemble solid and liquid to a certain extent. Meng et al. [198] found that NbCl5 and ZnMg intercalation in graphite results in pand n-type doping, respectively. To model clusters they used Gupta semiempirical potential and obtained globally optimal structures of those metal clusters employing GA with “cut-and-splice” crossover and triangle mutation. The optimized clusters were then placed between two layers of graphite, each of which contained 12 12 primitive cell and 288 atoms. The whole system was optimized with Universal potential and the electronic structures at the optimized doped cluster geometries were finally calculated using DFT. The atomic arrangement has an immense contribution to the catalytic activities of a nanoalloy. So, the development of a reliable method for predicting atomic arrangements has become increasingly important. Hong et al. [199] combined GA with first-principles electronic structure calculations to find the most stable configurations of elementary Aum and Agn clusters, as well as binary clusters AumAgn (5 12). The search led to the critical Au:Ag ratios for which the 2D-3D transition takes place. The study showed that Ag atoms prefer peripheral positions with lower coordination number while Au atoms tend to occupy central sites with higher coordination numbers. Johnston et al. [200] investigated small Sn-Bi clusters employing a DFT/GA method. Oh et al. [201] developed a combination of GA and classical molecular dynamics (MD) simulations to evaluate the optimal arrangement of multilayered core-shell structure for the Pt-Cu nano alloy with 1,654 atoms, regardless of the composition ratio. Using Tabu search in descriptor space and DFT Orel et al. [202] found the global minima of the neutral binary SnmPbn atomic clusters, 7 12, for all the possible stoichiometric ratios. Fournier et al. [203] found the minimum-energy structures of AgnRbn (n = 2-10) clusters by a combination of DFT and Tabu search. Chu et al. [204] described a new fragment-based evolutionary algorithm (EA) for de novo optimization, specifi-cally developed to handle organometallic and transition metal compounds. Au atoms preferentially segregates on the surface of Pd-Au nanoalloys. DFT calculations on fcc-type cubo-octahedral Pd-Au nanoparticles have indicated that the surface Pd atoms occupy (111) rather than (100) faces, thereby maximizing the number of relatively strong surface Pd-Au bonds [205]. Larger cohesive energy, lower electronegativity and smaller atomic radius (minimize bulk elastic strain) of Pd favors a Pd core and Au surface shell by lowering surface energy of the cluster. Besides electron transfer from Pd to more electro-negative Au, greater order of the Pd-Au than the Pd-Pd and Au-Au bonds favors Pd-Au mixing [96]. Ferrando et al. had considered a methodology based on extensive global-optimization within empirical potential models and subsequent DFT based local relaxation of low-energy structures pertaining to different structural motifs (or basins on the energy landscape) of gas-phase alloy nanoclusters [117, 207, 208, 209, 210, 211, 206] modelled by the many-body Gupta potential. Pittaway et al. [96] performed global optimization of Pd-Au bimetallic clusters using a GA, coupled with the Gupta many-body empirical potential (EP) to model inter-atomic interactions. Johnston and coworkers [100, 101, 94, 102, 117, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219] did a series of development on theoretical simulation of alloys and bimetallic clusters (e.g., Cu-Ag, Ni-Ag, Cu-Au, Ag-Pd, Ag-Au, Au-Pd) with GA coupled with the Gupta many-body potential. AgcorePdshell nano catalysts could be important in the development of micro polymer electrolyte membrane fuel cells for portable devices, and could also be applied to the promotion of other catalytic reactions [220]. Ag3Pd10 cluster was identified as a wedge-shaped nano-shell with Cs point group symmetry by a combination of GA global optimization and DFT calculation. Wu et al. [109] performed Global optimization of AgmPdn(m + n = 15) and Ag3mPd3812) clusters and came to the conclusion that silver atoms had a strong tendency towards segregating at the surface. Combined experimental/theoretical studies have indicated the preference for PdcoreAushell 5nm) and an “onion-like” Pd Au Pd configuration for larger Pd Au nanoalloys 12nm). Johnston et al. [101] have studied mixed coinage metal clusters, using the Gupta potential.

6.3 Crystal Structure Prediction

Evolutionary algorithms [60, 132, 133, 134, 135, 137, 138, 141, 142, 167, 170, 171, 173, 178, 221, 222, 223, 224, 225, 226] have provided reliable and powerful means of exploring the PES of crystals leading to the identification of the most stable, interesting and sometimes unexpected crystal structures. There are several other methods such as basin hopping [229, 227, 228, 142, 143], metadynamics [230, 231], PSO [178, 232], and simulated annealing [233] for performing similar search. For large dimensional problems, due to the complex nature of PES and the presence of high-energy disordered structures (random initial population) these algorithms often encounter serious challenges. Lyakhov et al. introduced an additional variation operator - the coordinate mutation or soft mutation [60] for the purpose of crystal structure prediction. Instead of using complete randomness in the mutation operator concerted mutation which directs the system to choices that have a higher probability to improve the fitness of the solution string is introduced. The idea of soft-mode mutation or softmutation [60, 134] is to move the atoms along the softest modes i.e., the lowest frequency eigenmodes that correspond to directions of lowest curvatures of the energy surface. The eigenvector corresponding to lowest non-zero eigenvalue determine the direction of softmutation. To calculate the softest modes one has to construct the dynamical matrix. It is often enough to have an approximate direction and sufficiently large mutation amplitude to arrive at a new low-energy structure. There is no need therefore of doing computationally expensive ab initio dynamical matrix calculation. Cheaper method for example dynamical matrix computed from bond hardness coefficients can be used. The soft mutation favor atoms with higher local order, Lyakhov et al. [134] constructed the first generation using pseudosubcells with fractional atomic occupancies to improve the order and diversity of the structures. They used a fingerprint function that improves the selection process through removing clones. Superconducting high pressure phase of germane [135], structural characterization of compressed silane [234], rhombohedral structure of superhard BC2N [235], superconducting structures of BC5 [236], superhard monoclinic polymorph of carbon [237], high-pressure orthorhombic polymorph of MgB2 stable above 190 GPa [238], high-pressure phases of CaLi2 [239], metallic structures of oxygen at pressures in the range of 100-250 GPa [137], two new hexagonal ultra-hard phases of WN2 [240], two unique high-pressure metallic phases of Stannane (SnH4) at superconducting temperatures of 15-22 K for the Ama2 phase at 120 GPa and 52-62 K for the P63/mmc phase at 200 GPa [241] and SiH4(H2)2 at high superconducting temperatures of 98-107 K at 250 GPa [242, 241] are some of the systems explored through ab initio evolutionary methodology. By combining PSO with first principles calculations, Xiang et al. [243] discovered a new metastable Si phase (Si20-T has a quasidirect gap of 1.55 eV) which could be a promising solar energy absorber. Oganov et al. [138] found a new boron phase (comprised of icosahedral B12 clusters and B2 pairs in a NaCl-type arrangement), stable between 19 and 89 GPa, and exhibiting evidence for charge transfer. In the context of reliably predicting ground state geometry, the Universal Structure Predictor: Evolutionary Xtallography (USPEX) [131, 132, 133, 141, 244] method has turned out to be extremely successful. Oganov et al. using USPEX methodology first demonstrated the possibility of computing hardness (at least for insulators and semiconductors) just from the crystal structure opens up the possibility of global optimization of hardness, aimed at the computational discovery of new superhard materials. They had introduced local measures of the quality of structure to locate defective regions in the crystal and used fingerprint niching, soft-mode mutation, symmetry- and pseudosymmetry-enabled generation of structures to greatly speed up the search for the global minimum [245, 246]. In a recent article [247] they have summarized the different applications of the USPEX method as a tool for crystal structure prediction and showed that this method has an enormous applications in both computational materials design and studies of matter at extreme conditions. In a recent application Oganov et al. showed in conjunction with the first-principles calculations, evolutionary algorithm will lead to the discovery of novel dielectrics [248]. In another application they successfully predicted energetically lower novel 2D boron structure than the -sheet structure [249], stable hafnium carbides [250]. Woodley et al. predicted ground state geometry for a series of transparent conducting indium oxide [170] and zirconia nanoclusters [171] by combining evolutionary algorithm with classical interatomic potential and quantum chemical models. Wu et al. showed that structure exploration by classical potentials with the accuracy of density functional theory is an efficient scheme for complex crystal structure prediction [251]. In another work they favorably compared Lamarckian schemes within evolutionary algorithms with Darwinian schemes in the context of minimum energy structure prediction of titania phases [173]. O’Keeffe [252] describes several difficulties encounted in the crystal structure prediction. Study on GA and minima hopping method based crystal structure prediction reveals that for relatively smaller sizes both methods show comparable efficiency while for larger systems GA becomes advantageous over minima hopping [229]. Oganov et al. explained how and why evolutionary crystal structure prediction works the way they do [253]. Using adaptive GA Zhao et al. [254] studied the structures and stabilities of the alkaline earth metal peroxide XO2. They also predicted complex crystal structures of the orthorhombic, rhombohedral, and hexagonal polymorphs close to the Zr2Co11 intermetallic compound [255]. Li et al. [256] explored the high-pressure crystal structures of Mg by the PSO algorithm. Wang et al. have developed a software package titled Crystal Structure Analysis by Particle Swarm Optimization (CALYPSO), enforcing symmetry constraints on structure generation, bond characterization matrix for elimination of similar structures, introducing partial random structures per generation for enhancing structural diversity, and penalty function, etc. [178, 232] for predicting crystal structure from random initial starting geometries. They have applied CALYPSO to crystal structure prediction, earth and planetary materials [257, 258, 259, 260, 261, 262, 263, 264]. Luo et al. [265] predicted new stable structures of 2D boron-carbon compounds for a wide range of boron concentrations by the PSO algorithm implemented in the CALYPSO code. XtalOpt [266] and EVO [223] are the two other packages that have proved to be useful in this context. Zhu et al. proposed a method based on metadynamics and evolutionary algorithms [230] and found stable and metastable states for Al2SiO5, SiO2, MgSiO3, and carbon clusters [231] from reasonable initial structures providing insight into the mechanisms of phase transitions. They designed an evolutionary algorithm based method [267] to automatically explore low energy surface reconstructions with variable surface atoms and reconstruction cells and illustrated it by the identification of N3 trimeric reconstruction of GaN(10¯11) surfaces. Liu [268] constructed a Multi-algorithm-collaborative Universal Structureprediction Environment (Muse) to efficiently find the stable and metastable structures of materials under given conditions. In Muse the evolutionary algorithm was coupled with the simulated annealing and the basin hopping algorithms. With the inclusion of two new variation operators, slip and twist the performance of Muse was greatly enhanced. Fadda et al. [269] described an evolutionary algorithm based on symmetry-preserving and symmetry-breaking mutations for the exploration of the space of the conjugacy classes of crystal structures. Johnston et al. have implemented and parameterized both Darwinian and Lamarkian GA search to identify structure of clusters from experimental STEM images [270]. Meredig et al. [271] had given considerable theoretical and computational efforts for the development of first-principles-assisted structure solution (FPASS) that automatically solve crystal structures. Nguyen et al. [272] performed genetic algorithm to find two meta-stable Si-IX phase with good experimental agreement in the lattice parameters.

7 Future Directions

The evolutionary computing techniques and their hybrids have registered impressive success in the elucidation of minimum energy structures of atomic clusters and crystalline solids. These methods have not benn explored that extensively for weakly bonded molecular clusters and crystals. We anticipate rapid growth of hybrid EC based exploration of such structures. Since these surfaces are dominated by shallow minima separated by small barriers, special mutation and crossover operators would have to be designed for exhaustive exploration of the PES and locating global as well as local close-lying minima. The information on such structures and their energies could then lead to the correct prediction of thermally averaged properties of such species.

The use of EC for directly solving the molecular Schrodinger equation has been very limited. The density based method described in reference [56] has the potential to be a viable, even a superior alternative for computing the equilibrium molecular structure and the corresponding one electron density resolved in a given basis. We anticipate extensive exploration along that line in the future years. Intelligent interfacing of the standard electronic structure codes and EC codes could prove instrumental for further growth of EC techniques in exploration of molecular structures. One area of great promise where EC could decisively prove superior is in designing molecules and materials with targetted properties. Such problems can be cast in the mould of multimodal optimization problem [63] which can be efficiently handled by the EC techniques using for example, the prey-predator model for multi-objective constrained optimization. Alongwith artificial neural network, EA can meet tate-of-the-art methods for powerful quantitative structure-property relationships modeling [273, 274]. We anticipate rapid growth in such studies in the next decade.

8 Conclusions

The last decade has seen impressive growth of EC techniques for exploring PES of clusters and crystals. The review highlights the important role played by soft-computing techniques, especially the EC methods in the elucidation of large scale structures on complex potential energy landscapes. When combined with DFT, such methods can often produce surprising results. The evolutionary computing techniques or for that matter the whole gamut of Soft-computing methods have not yet been explored exhaustively in the context of solving Schr¨odinger equation for even few electron systems although their potentials as viable tools have been demonstrated. One anticipates that the cross-talk between EC and computational chemistry would be productive and new hybrid algorithms for exploring electronic structures of atoms, molecules and clusters will evolve out of it.

Acknowledgement:

S.P.B. thanks the DAE, Government of India, for the award of Raja Ramanna Felloship, the authorities of IIT, Bombay for hosting the fellowship, and Professor B. L. Tembe, Department of Chemistry, IIT, Bombay for many helpful discussions. We thank Dr. Pinaki Chowdhury (Calcutta University) for helpful discussions.

References

1. Antoniou, A., Lu, W. S. Practical Optimization: Algorithms and Engineering Applications, Springer: 2010.

2. Deb, K. Evolutionary Computation 1999, 7, 205–230.

3. Weise, T. Global Optimization Algorithms – Theory and Application –, 2008.

4. Weise, T., Zapf, M., Chiong, R., Nebro, A. Why Is Optimization Difficult?. In NatureInspired Algorithms for Optimisation, Vol. 193, Chiong, R., Ed., Springer Berlin Heidelberg: 2009.

5. Weise, T., Chiong, R., T´ang, K. Journal of Computer Science and Technology (JCST) 2012, 27, 907–936.

6. Jin, Y., Branke, J. Evolutionary Computation, IEEE Transactions on 2005, 9, 303–317.

7. Coello, C. A. C., Lamont, G. B., van Veldhuizen, D. A. Evolutionary Algorithms for Solving Multi-Objective Problems, Genetic and Evolutionary Computation Springer: 2007.

8. Deb, K. Multi-Objective Optimization Using Evolutionary Algorithms, Wiley paperback series Wiley: 2009.

9. Holland, J. H. Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence, Bradford Books, MIT Press: 1992.

10. Mitchell, M. An Introduction to Genetic Algorithms, A Bradford book, MIT Press: 1998.

11. Brownlee, J. Clever Algorithms: Nature-Inspired Programming Recipes, Lulu Com: 2011.

12. Tang, Z., Bagchi, K. K. Computer and Information Science 2010, 60–71.

13. Kirkpatrick, S., Gelatt, C. D., Vecchi, M. P. Science 1983, 220, 671–680.

14. Talbi, E. G. Metaheuristics: From Design to Implementation, Wiley Series on Parallel and Distributed Computing Wiley: 2009.

15. Hartke, B. Wiley Interdisciplinary Reviews: Computational Molecular Science 2011, 1, 879– 887.

16. Darwin, C. On the Origin of the Species by Means of Natural Selection: Or, The Preservation of Favoured Races in the Struggle for Life, John Murray: 1859.

17. Alba, E., Dorronsoro, B. Cellular Genetic Algorithms, Operations Research/Computer Science Interfaces Series Springer: 2009.

18. Tomassini, M. Spatially Structured Evolutionary Algorithms: Artificial Evolution in Space and Time, Natural Computing Series Springer: 2010.

19. Goldberg, D. E. Genetic Algorithms, Pearson Education: 2013.

20. Goldberg, D. E., Deb, K. A comparative analysis of selection schemes used in genetic algo- rithms. In Foundations of Genetic Algorithms, Morgan Kaufmann: 1991.

21. Goh, K. S., Lim, A., Rodrigues, B. Artificial Intelligence Review 2003, 19, 123–152.

22. Agrawal, A. F. Nature 2001, 411, 692–695.

23. Campelo, F., Guimaraes, F. G., Igarashi, H., Ramirez, J. A. Magnetics, IEEE Transactions on 2005, 41, 1736–1739.

24. Syswerda, G. Uniform Crossover in Genetic Algorithms. In Proceedings of the 3rd International Conference on Genetic Algorithms, Morgan Kaufmann Publishers Inc.: San Francisco, CA, USA, 1989.

25. Bao, Z., Watanabe, T. A novel Genetic Algorithm with cell crossover for circuit design opti- mization. In IEEE International Symposium on Circuits and Systems, 2009.

26. Wright, A. H. Genetic Algorithms for Real Parameter Optimization. In Foundations of Genetic Algorithms, Morgan Kaufmann: 1991.

27. Tsutsui, S., Yamamura, M., Higuchi, T. Multi-parent recombination with simplex crossover in real coded genetic algorithms. In Proceedings of the genetic and evolutionary computation conference, Vol. 1, 1999.

28. Michalewicz, Z. Genetic Algorithms + Data Structures = Evolution Programs, Artificial intelligence Springer: 1996.

29. Eshelman, L. J., Schaffer, J. D. Real-Coded Genetic Algorithms and Interval-Schemata.. In FOGA, Whitley, L. D., Ed., Morgan Kaufmann: 1992.

30. Radcliffe, N. J. Complex Systems 1991, 5, 183–205.

31. M¨uhlenbein, H., Schlierkamp-Voosen, D. Evol. Comput. 1993, 1, 25–49.

32. Guvenc, U. Scientific Research and Essays 2010, 5, 2451–2456.

33. Lozano, M., Herrera, F., Krasnogor, N., Molina, D. Evolutionary Computation 2004, 12, 273–302.

34. Voigt, H.-M., M\ddotuhlenbein, H., Cvetkovic, D. Fuzzy recombination for the Breeder Genetic Algorithm. In Proc. Sixth Int. Conf. on Genetic Algorithms, Morgan Kaufmann Publishers: 1995.

35. Deb, K., Agrawal, R. B. Complex Systems 1995, 9, 115–148.

36. Deb, K., Anand, A., Joshi, D. Evolutionary computation 2002, 10, 371–395.

37. Takahashi, O., Kobayashi, S. An adaptive neighboring search using crossover-like mutation for multi modal function optimization. In IEEE International Conference on Systems, Man, and Cybernetics, Vol. 1, 2001.

38. Deep, K., Thakur, M. Applied Mathematics and Computation 2007, 188, 895–911.

39. Storn, R., Price, K. Journal of global optimization 1997, 11, 341–359.

40. Whitley, D., Hains, D., Howe, A. Tunneling between optima: partition crossover for the traveling salesman problem. In Proceedings of the 11th Annual conference on Genetic and evolutionary computation, GECCO ’09 ACM: New York, NY, USA, 2009.

41. Whitley, D., Hains, D., Howe, A. A Hybrid Genetic Algorithm for the Traveling Salesman Problem Using Generalized Partition Crossover. In Parallel Problem Solving from Nature, PPSN XI, Vol. 6238, Schaefer, R., Cotta, C., Kolodziej, J., Rudolph, G., Eds., Springer Berlin Heidelberg: 2010.

42. Schlierkamp-voosen, D., M\ddotuhlenbein, H. Strategy Adaptation by Competing Subpopulations. In Parallel Problem Solving from Nature (PPSN III), Springer-Verlag: 1994.

43. Ono, I., Kita, H., Kobayashi, S. A Real-coded Genetic Algorithm using the Unimodal Normal Distribution Crossover. In Advances in Evolutionary Computing, Ghosh, A., Tsutsui, S., Eds., Natural Computing Series Springer Berlin Heidelberg: 2003.

44. Herrera, F., Lozano, M., Verdegay, J. L. Artificial Intelligence Review 1998, 12, 265–319.

45. Arumugam, M. S., Rao, M. V. C., Palaniappan, R. Applied Soft Computing 2005, 6, 38–52.

46. Chang, W.-D. Computers & Mathematics with Applications 2006, 51, 1437–1444.

47. Kaya, Y., Uyar, M., Tekin, R. CoRR 2011, abs/1105.0355,.

48. Gwiazda, T. D. Crossover for single-objective numerical optimization problems, volume 1 of Genetic algorithms reference Tomasz Gwiazda: 2006.

49. Baeck, T., Fogel, D. B., Michalewicz, Z. Handbook of Evolutionary Computation, Taylor & Francis: 1997.

50. Iwamatsu, M. Computer Physics Communications 2002, 147, 729–732.

51. Tin´os, R., Yang, S. Soft Computing 2011, 15, 1523–1549.

52. Davis, L. Handbook of Genetic Algorithms, International Thomson Computer Press: 1996.

53. Dorronsoro, B., Bouvry, P. IEEE Transactions on Evolutionary Computation 2011, 15, 67–98.

54. Chuang, Y.-C., Chen, C.-T. A study on real-coded genetic algorithm for process optimization using ranking selection, direction-based crossover and dynamic mutation. In IEEE Congress on Evolutionary Computation (CEC), 2011.

55. Sharma, R., Bhattacharyya, S. P. Direct search for wave operator by a Genetic Algorithm (GA): Route to few eigenvalues of a Hamiltonian. In IEEE Congress on Evolutionary Computation (CEC), 2007.

56. Sarkar, K., Sharma, R., Bhattacharyya, S. P. Journal of Chemical Theory and Computation 2010, 6, 718–726.

57. Ling, S. H., Leung, F. H. F. Soft Computing 2007, 11, 7–31.

58. Voigt, H.-M., Anheyer, T. Modal mutations in evolutionary algorithms. In IEEE World Congress on Computational Intelligence., Proceedings of the First IEEE Conference on Evolutionary Computation, 1994.

59. Oganov, A. R., Valle, M. The Journal of Chemical Physics 2009, 130, 104504.

60. Lyakhov, A. O., Oganov, A. R., Valle, M. Computer Physics Communications 2010, 181, 1623–1632.

61. Vlachos, D. S., Parousis-Orthodoxou, K. J. Population Induced Instabilities in Genetic Algo- rithms for Constrained Optimization. In International Conference on Mathematical Modelling in Physical Sciences, Vol. 410, 2013.

62. Sarkar, K., Bhattacharyya, S. P. The Journal of Chemical Physics 2013, 139, 074106.

63. Sarkar, K., Sharma, R., Bhattacharyya, S. P. International Journal of Quantum Chemistry 2012, 112, 1547–1558.

64. Sarkar, K., Bhattacharyya, S. P. AIP Conference Proceedings 2013, 1512,.

65. Goldberg, D. E., Deb, K., Clark, J. H. Complex systems 1992, 6, 333–362.

66. M¨uhlenbein, H., Schlierkamp-Voosen, D. Evolutionary Computation 1993, 1, 335–360.

67. Arabas, J., Michalewicz, Z., Mulawka, J. GAVaPS-a genetic algorithm with varying pop- ulation size. In Evolutionary Computation, 1994. IEEE World Congress on Computational Intelligence., Proceedings of the First IEEE Conference on, 1994.

68. Nandy, S., Sharma, R., Bhattacharyya, S. P. Applied Soft Computing 2011, 11, 3946–3961.

69. Dorronsoro, B., Bouvry, P. IEEE Transactions on Evolutionary Computation 2011, 15, 67–98.

70. Mininno, E., Neri, F., Cupertino, F., Naso, D. Evolutionary Computation, IEEE Transactions on 2011, 15, 32–54.

71. Das, S., Suganthan, P. N. Evolutionary Computation, IEEE Transactions on 2011, 15, 4–31.

72. Kennedy, J., Eberhart, R. Particle swarm optimization. In Neural Networks, 1995. Proceedings., IEEE International Conference on, Vol. 4, 1995.

73. Das, S., Abraham, A., Konar, A. Particle Swarm Optimization and Differential Evolution Algorithms: Technical Analysis, Applications and Hybridization Perspectives. In Advances of Computational Intelligence in Industrial Systems, Vol. 116, Liu, Y., Sun, A., Loh, H., Lu, W., Lim, E.-P., Eds., Springer Berlin Heidelberg: 2008.

74. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E. The Journal of Chemical Physics 1953, 21, 1087–1092.

75. Muller, C. R., Kathriarachchi, V., Schuch, M., Maass, P., Petkov, V. G. Physical Chemistry Chemical Physics 2010, 12, 10444–10451.

76. Fang, X. W., Huang, L., Wang, C. Z., Ho, K. M., Ding, Z. J. Journal of Applied Physics 2014, 115, 061905.

77. Biring, S. K., Chaudhury, P. Chemical Physics 2010, 377, 46–53.

78. Biring, S. K., Chaudhury, P. Chemical Physics 2012, 400, 198–206.

79. Zeiri, Y., Fattal, E., Kosloff, R. The Journal of Chemical Physics 1995, 102, 1859–1862.

80. Chaudhury, P., Bhattacharyya, S. P. Chemical Physics Letters 1998, 296, 51–60.

81. Nakanishi, H., Sugawara, M. Chemical Physics Letters 2000, 327, 429–438.

82. Saha, R., Chaudhury, P., Bhattacharyya, S. P. Physics Letters A 2001, 291, 397–406.

83. Saha, R., Bhattacharyya, S. P., Taylor, C. D., Zhao, Y., Cundari, T. R. International Journal of Quantum Chemistry 2003, 94, 243–250.

84. Saha, R., Bhattacharyya, S. P. Journal of Theoretical and Computational Chemistry 2004, 3, 325–337.

85. S¸AH˙IN, M., Atav, ¨U., Tomak, M. Turkish Journal of Physics 2006, 30, 253–275.

86. NANDY, S., CHAUDHURY, P., SHARMA, R., BHATTACHARYYA, S. P. Journal of Theoretical and Computational Chemistry 2008, 07, 977–987.

87. Deaven, D. M., Ho, K. M. Physical Review Letter 1995, 75, 288–291.

88. Dugan, N., Erkoc, S. Computational Materials Science 2009, 45, 127–132.

89. Chaudhury, P., Saha, R., Bhattacharyya, S. P. Chemical Physics 2001, 270, 277–285.

90. Chaudhury, P., Bhattacharyya, S. P. Chemical Physics 1999, 241, 313–325.

91. Chaudhury, P., Bhattacharyya, S. P., Quapp, W. Chemical Physics 2000, 253, 295–303.

92. Moseler, M., H¨akkinen, H., Barnett, R. N., Landman, U. Physical Review Letter 2001, 86, 2545–2548.

93. Nø rskov, J. K., Bligaard, T., Rossmeisl, J., Christensen, C. H. Nature chemistry 2009, 1, 37–46.

94. Ferrando, R., Jellinek, J., Johnston, R. L. Chemical Reviews 2008, 108, 845–910.

95. Ferrando, R. Computational Methods for Predicting the Structures of Nanoalloys. In Nanoalloys, Alloyeau, D., Mottet, C., Ricolleau, C., Eds., Engineering Materials Springer London: 2012.

96. Pittaway, F., Paz-Borbon, L. O., Johnston, R. L., Arslan, H., Ferrando, R., Mottet, C., Barcaro, G., Fortunelli, A. The Journal of Physical Chemistry C 2009, 113, 9141–9152.

97. Cheng, L., Feng, Y., Yang, J., Yang, J. The Journal of Chemical Physics 2009, 130, 214112.

98. Huang, W., Lai, X., Xu, R. Chemical Physics Letters 2011, 507, 199–202.

99. Wilcoxon, J. P., Abrams, B. L. Chem. Soc. Rev. 2006, 35, 1162–1194.

100. Paz-Borb´on, L. O., Johnston, R. L., Barcaro, G., Fortunelli, A. The Journal of Physical Chemistry C 2007, 111, 2936–2941.

101. Paz-Borbon, L. O., Johnston, R. L., Barcaro, G., Fortunelli, A. The Journal of Chemical Physics 2008, 128, 134517.

102. Chen, F., Johnston, R. L. ACS Nano 2008, 2, 165–175.

103. Shao, X., Cheng, L., Cai, W. Journal of Computational Chemistry 2004, 25, 1693–1698.

104. Grosso, A., Locatelli, M., Schoen, F. Mathematical Programming 2007, 110, 373–404.

105. Rossi, G., Ferrando, R. Journal of Physics: Condensed Matter 2009, 21, 84208.

106. Oakley, M. T., Johnston, R. L. Journal of Chemical Theory and Computation 2013, 9, 650– 657.

107. Cerbelaud, M., Ferrando, R., Barcaro, G., Fortunelli, A. Physical Chemistry Chemical Physics 2011, 13, 10232–10240.

108. Wu, X., Cai, W., Shao, X. Journal of Computational Chemistry 2009, 30, 1992–2000.

109. Wu, X., Wu, Y., Kai, X., Wu, G., Chen, Y. Chemical Physics 2011, 390, 36–41.

110. Wu, X., Huang, C., Sun, Y., Wu, G. Chemical Physics 2013, 415, 69–75.

111. Wu, X., Wu, G., Chen, Y., Qiao, Y. The Journal of Physical Chemistry A 2011, 115, 13316– 13323.

112. Pereira, F. B., Marques, J. M. C., Leitao, T., Tavares, J. Analysis of Locality in Hybrid Evolutionary Cluster Optimization. In IEEE Congress on Evolutionary Computation (CEC), 2006.

113. Pereira, F. B., Marques, J., Leitao, T., Tavares, J. Designing Efficient Evolutionary Algo- rithms for Cluster Optimization: A Study on Locality. In Advances in Metaheuristics for Hard Optimization, Siarry, P., Michalewicz, Z., Eds., Natural Computing Series Springer Berlin Heidelberg: 2008.

114. Pereira, F. B., Marques, J. M. C. A self-adaptive evolutionary algorithm for cluster geometry optimization. In Hybrid Intelligent Systems, 2008. HIS’08. Eighth International Conference on, 2008.

115. Pereira, F., Marques, J. Evolutionary Intelligence 2009, 2, 121–140.

116. Marques, J. M. C., Pereira, F. B., Leitao, T. The Journal of Physical Chemistry A 2008, 112, 6079–6089.

117. Ferrando, R., Fortunelli, A., Johnston, R. L. Phys. Chem. Chem. Phys. 2008, 10, 640–649.

118. Schwerdtfeger, P., Gaston, N., Krawczyk, R. P., Tonner, R., Moyano, G. E. Physical Review B 2006, 73, 64112.

119. Marques, J. M. C., Pais, A. A. C. C., Abreu, P. E. Journal of Computational Chemistry 2010, 31, 1495–1503.

120. Cassioli, A., Locatelli, M., Schoen, F. Optimization Methods and Software 2009, 24, 819– 835.

121. Alberto Fernandez-Lima, F., Vilela Neto, O. P., Silva Pimentel, A., Pacheco, M. A. C., Ponciano, C. R., Nascimento, M. A. C., da Silveira, E. F. The Journal of Physical Chemistry A 2009, 113, 15031–15040.

122. Fernandez-Lima, F. A., VilelaNeto, O. P., Pimentel, A. S., Ponciano, C. R., Pacheco, M. A. C., Nascimento, M. A. C., da Silveira, E. F. The Journal of Physical Chemistry A 2009, 113, 1813–1821.

123. Pullan, W. Unbiased geometry optimisation of Morse atomic clusters. In IEEE Congress on Evolutionary Computation (CEC), 2010.

124. Oakley, M. T., Wales, D. J., Johnston, R. L. The Journal of Physical Chemistry B 2011, 115, 11525–11529.

125. Buck, U., Pradzynski, C. C., Zeuch, T., Dieterich, J. M., Hartke, B. Physical Chemistry Chemical Physics 2014, 16, 6859–6871.

126. Larsson, H. R., Duin, A. C., Hartke, B. Journal of computational chemistry 2013, 34, 2178– 2189.

127. HORNTON, A. R. T., WEINHART, T., OGARKO, V., LUDING, S. Computer methods in materials science 2013, 13, 197 – 212.

128. Li, Y., Hartke, B. ChemPhysChem 2013, 14, 2678–2686.

129. Alexandrova, A. N. The Journal of Physical Chemistry A 2010, 114, 12591–12599.

130. Doll, K., Schon, J. C., Jansen, M. The Journal of Chemical Physics 2010, 133, 24107.

131. Glass, C. W., Oganov, A. R., Hansen, N. Computer Physics Communications 2006, 175, 713–720.

132. Oganov, A. R., Glass, C. W. Journal of Physics: Condensed Matter2008, 20, 64210.

133. Woodley, S. M., Catlow, R. Nature materials 2008, 7, 937–946.

134. Lyakhov, A. O., Oganov, A. R., Stokes, H. T., Zhu, Q. Computer Physics Communications 2013, 184, 1172–1182.

135. Gao, G., Oganov, A. R., Bergara, A., Martinez-Canales, M., Cui, T., Iitaka, T., Ma, Y., Zou, G. Physical Review Letters 2008, 101, 107002.

136. Oganov, A. R., Ono, S., Ma, Y., Glass, C. W., Garcia, A. Earth and Planetary Science Letters 2008, 273, 38–47.

137. Ma, Y., Oganov, A. R., Glass, C. W. Physical Review B 2007, 76, 64101.

138. Oganov, A. R., Chen, J., Gatti, C., Ma, Y., Ma, Y., Glass, C. W., Liu, Z., Yu, T., Kurakevych, O. O., Solozhenko, V. L. Nature 2009, 457, 863–867.

139. Valle, M., Oganov, A. R. Acta Crystallographica Section A 2010, 66, 507–517.

140. Ma, Y., Eremets, M., Oganov, A. R., Xie, Y., Trojan, I., Medvedev, S., Lyakhov, A. O., Valle, M., Prakapenka, V. Nature 2009, 458, 182–185.

141. Zhu, Q., Oganov, A. R., Glass, C. W., Stokes, H. T. Acta Crystallographica Section B 2012, 68, 215–226.

142. Sch¨onborn, S. E., Goedecker, S., Roy, S., Oganov, A. R. The Journal of chemical physics 2009, 130, 144108.

143. Amsler, M., Goedecker, S. The Journal of Chemical Physics 2010, 133, 224104.

144. Sicher, M., Mohr, S., Goedecker, S. The Journal of Chemical Physics 2011, 134, 44106.

145. Oakley, M. T., Johnston, R. L., Wales, D. J. Physical Chemistry Chemical Physics 2013, 15, 3965–3976.

146. Wales, D. J., Carr, J. M. Journal of Chemical Theory and Computation 2012, 8, 5020–5034.

147. Wales, D. J., Head-Gordon, T. The Journal of Physical Chemistry B 2012, 116, 8394–8411.

148. Chakrabarti, D., Wales, D. J. Phys. Chem. Chem. Phys. 2009, 11, 1970–1976.

149. Froltsov, V. A., Reuter, K. Chemical Physics Letters 2009, 473, 363–366.

150. de Souza, V. K., Wales, D. J. The Journal of Chemical Physics 2009, 130, 194508.

151. Parodi, D., Ferrando, R. Physics Letters A 2007, 367, 215–219.

152. Marques, J. M. C., Pereira, F. B. Journal of Computational Chemistry 2013, 34, 505–517.

153. Marques, J. M. C., Pereira, F. B. Chemical Physics Letters 2010, 485, 211–216.

154. Pereira, F. B., Marques, J. M. C. Towards an effective evolutionary approach for binary Lennard-Jones clusters. In IEEE Congress on Evolutionary Computation (CEC), 2010.

155. Kolossv´ary, I., Bowers, K. J. Physical Review E 2010, 82, 56711.

156. El Dor, A., Clerc, M., Siarry, P. Computational Optimization and Applications 2012, 53, 271–295.

157. Dieterich, J. M., Hartke, B. Journal of Computational Chemistry 2011, 32, 1377–1385.

158. Dieterich, J. M., Hartke, B. Molecular Physics 2010, 108, 279–291.

159. Dzhurakhalov, A. A., Atanasov, I., Hou, M. Physical Review B 2008, 77, 115415.

160. Leitao, A., Pereira, F. B., Machado, P. Enhancing cluster geometry optimization with Island Models. In IEEE Congress on Evolutionary Computation (CEC), 2012.

161. Carstensen, N. O., Dieterich, J. M., Hartke, B. Physical Chemistry Chemical Physics 2011, 13, 2903–2910.

162. Dieterich, J. M., Hartke, B. arXiv preprint arXiv:1207.4318 2012, .

163. Stepanenko, S., Engels, B. The Journal of Physical Chemistry A 2009, 113, 11699–11705.

164. Stepanenko, S., Engels, B. Journal of Computational Chemistry 2008, 29, 768–780.

165. Daskin, A., Kais, S. Molecular Physics 2011, 109, 761–772.

166. Deep, K., Madhuri, Liquid-Drop-Like Multi-Orbit Topology Versus Ring Topology in PSO for Lennard-Jones Problem. In Proceedings of Seventh International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA 2012) SE - 20, Vol. 202, Bansal, J. C., Singh, P., Deep, K., Pant, M., Nagar, A., Eds., Springer India: 2013.

167. Oleksy, K., Karlick`y, F., Kalus, R. The Journal of chemical physics 2010, 133, 164314.

168. Huber, B., Moseler, M., Kostko, O., V.Issendorff, B. Physical Review B 2009, 80, 235425.

169. Huang, X., Sai, L., Jiang, X., Zhao, J. The European Physical Journal D 2013, 67, 1–7.

170. Walsh, A., Woodley, S. M. Physical Chemistry Chemical Physics 2010, 12, 8446–8453.

171. Woodley, S. M., Hamad, S., Catlow, C. R. A. Physical Chemistry Chemical Physics 2010, 12, 8454–8465.

172. Kim, S., Sohn, K.-S., Pyo, M. ACS Combinatorial Science 2011, 13, 101–106.

173. Woodley, S. M., Catlow, C. R. A. Computational Materials Science 2009, 45, 84–95.

174. Marques, J. M. C., Llanio-Trujillo, J. L., Albert´ı, M., Aguilar, A., Pirani, F. The Journal of Physical Chemistry A 2012, 116, 4947–4956.

175. Zhu, Q., Oganov, A. R., Lyakhov, A. O. Physical Chemistry Chemical Physics 2013, 15, 7696–7700.

176. Zhu, Q., Jung, D. Y., Oganov, A. R., Glass, C. W., Gatti, C., Lyakhov, A. O. Nature chemistry 2013, 5, 61–65.

177. Al-Sunaidi, A. A., Sokol, A. A., Catlow, C. R. A., Woodley, S. M. The Journal of Physical Chemistry C 2008, 112, 18860–18875.

178. Wang, Y., Lv, J., Zhu, L., Ma, Y. Physical Review B 2010, 82, 94116.

179. Li, R., Cheng, L. Computational and Theoretical Chemistry 2012, 996, 125–131.

180. Avaltroni, F., Corminboeuf, C. Journal of Computational Chemistry 2012, 33, 502–508.

181. Pal, S., Sharma, R., Goswami, B., Sarkar, P. Chemical Physics Letters 2009, 467, 365–368.

182. Pal, S., Sharma, R., Goswami, B., Sarkar, P., Bhattacharyya, S. P. The Journal of chemical physics 2009, 130, 214703.

183. Llanio-Trujillo, J. L., Marques, J. M. C., Pereira, F. B. The Journal of Physical Chemistry A 2011, 115, 2130–2138.

184. Do, H., Besley, N. A. The Journal of chemical physics 2012, 137, 134106.

185. Addicoat, M. A., Page, A. J., Brain, Z. E., Flack, L., Morokuma, K., Irle, S. Journal of Chemical Theory and Computation 2012, 8, 1841–1851.

186. Cai, Z., Cao, Q., Du, S., Ji, X. The Genetic algorithms in optimization of silicon clusters. In Natural Computation (ICNC), 2012 Eighth International Conference on, 2012.

187. Zhang, J., Wang, C.-Z., Ho, K.-M. Physical Review B 2009, 80, 174102.

188. Yao, Y., Tse, J. S., Tanaka, K.