Dalton Transactions

Evolving better nanoparticles: Genetic algorithms for optimising cluster geometries

Roy L. Johnston
School of Chemistry, University of Birmingham, Edgbaston, Birmingham, UK B15 2TT

A review is presented of the design and application of genetic algorithms for the geometry optimisation of clusters and nanoparticles, where the interactions between atoms, ions or molecules are described by a variety of potential energy functions. A general introduction to genetic algorithms is followed by a detailed description of the genetic algorithm program that we have developed to identify the lowest energy isomers for a variety of atomic and molecular clusters. Examples are presented of its application to model Morse clusters, ionic MgO clusters and bimetallic nanoalloy clusters. Finally, a number of recent innovations and possible future developments are discussed.

Roy L. Johnston

Roy Johnston was born in Aberdeen in 1961 and received his BA (1983) and DPhil (1986) from the University of Oxford, where he was a member of St. Catherine's College. His DPhil, on theoretical aspects of cluster molecules, was carried out under the supervision of Professor D. M. P. Mingos. From 1987–1989, he was a SERC/NATO postdoctoral fellow in the USA, spending a year each at Cornell University (working with Professor Roald Hoffmann on theoretical solid state chemistry) and the University of Arizona (working with Professor Dennis Lichtenberger on gas phase and surface photoelectron spectroscopy). He returned to the UK in 1989 as a Royal Society University Research Fellow at the University of Sussex, where he collaborated with Professor John Murrell on the development and application of many-body potentials for modelling clusters, surfaces and solids. He was appointed to a lectureship in inorganic chemistry at the University of Birmingham in 1995 and was promoted to Senior Lecturer in 2002 and to a Reader in Computational Chemistry in 2003. His research interests lie in theoretical and computational chemistry, with particular recent interests being: the application of genetic algorithms and other evolutionary techniques to optimization problems in cluster chemistry, crystallography and protein folding; and the study of the geometries and segregation properties of bimetallic "nanoalloy" clusters. He is a Fellow of the Royal Society of Chemistry and is currently Secretary of the RSC Theoretical Chemistry Group.

Clusters are aggregates of between a few and many millions of atoms or molecules. They may consist of identical atoms, or molecules, or two or more different species and can be studied in a number of media, such as molecular beams, the vapour phase, in colloidal suspensions and isolated in inert matrices or on surfaces.^1,2 Interest in clusters arises, in part, because they constitute a new type of material which may have properties which are distinct from those of discrete molecules or bulk matter: for example, some metals (e.g. palladium) which are non-magnetic in the solid state can give rise to non-zero magnetic moments in discrete clusters.^3,4 Another reason for the interest in clusters is the size-dependent evolution of cluster properties.^1,5,6 Examples of some different types of clusters (fullerenes, metal clusters, molecular clusters and ionic clusters) are shown in Fig. 1.

	Fig. 1 Examples of some cluster types.

Theory plays an important role in cluster science because many cluster properties are difficult to measure directly and spectroscopic and mass spectrometric data often need to be interpreted in terms of theoretical models.^1,7 Since, for large clusters (of hundreds or thousands of atoms) ab initio electronic structure calculations are still, at present, computationally expensive, there has been much interest in developing empirical atomistic potentials for the simulation of such species.⁸

Whether one is using empirical atomistic potentials, ab initio Molecular Orbital (MO) theory or Density Functional Theory (DFT) to describe the bonding in clusters, one of the prime objectives is to find, for a given cluster size, the arrangement of atoms (or ions or molecules) corresponding to the lowest potential energy—i.e. the global minimum (GM) on the potential energy hypersurface.^9,10 As it is possible that a number of different states (isomers or local minima) will be thermally populated at a given finite temperature, one may pose the question Why is it important to locate the global minimum? Clusters corresponding to global minima (or low-lying local minima) are the most likely candidates for the most probable structure formed in a cluster experiment,¹¹ though, depending on experimental conditions, the structures observed experimentally may be kinetic (metastable), rather than thermodynamic products.¹² It is also important, however, to locate the global minimum from the point of view of testing the accuracy and how physically reasonable a particular potential energy function (or other theoretical model) is—for example, it is no use trying to model carbon clusters with a potential that predicts close packed structures!

As the number of minima rises exponentially with increasing cluster size, finding the global minimum can be very difficult. It has been found that traditional Monte Carlo (MC) and Molecular Dynamics (MD) Simulated Annealing (SA) approaches often encounter difficulties finding global minima for particular types of interatomic interactions (such as the short ranged Morse potential discussed below).¹³ For this reason, genetic algorithms have found increasing use in the area of finding global minima for clusters—i.e. cluster geometry optimisation.

The genetic algorithm (GA)^14–16 is a search technique based on the principles of natural evolution. It uses operators that are analogues of the evolutionary processes of mating (or crossover at the gene level), mutation and natural selection to explore multi-dimensional parameter spaces. The GA belongs to the class of evolutionary algorithms, which also includes evolution strategies, differential evolution and genetic programming.¹⁷

1.2.1 Terminology and definitions of GA operators.

In principal, a GA can be applied to any problem where the variables to be optimised (genes) can be encoded to form a string (chromosome). Each string represents a trial solution of the problem. By analogy with biology, the values of the individual variables are known as alleles. The relationship between GA chromosomes, genes and alleles is shown schematically in Fig. 2. In a typical GA, a population of individuals evolves for a certain number of generations—which is either fixed in advance or may depend on some convergence criterion. Further details on generic GAs and some specific implementations can be found in the standard texts.^14–16

	Fig. 2 Schematic representation of an individual in a generic optimisation GA. The genes (yellow circles) represent the variables in terms of which the optimization problem is defined and the alleles (symbols) are the values of these variables. The individual is represented by a chromosome or string of variables.

The initial population corresponds to the starting set of individuals which are to be evolved by the GA. These individuals are usually generated randomly, though it may sometimes be beneficial to use any available prior knowledge or chemical intuition in the generation of the initial population (taking care not to bias the search too much).

Fitness is an important concept for the operation of the GA. The fitness of a string is a measure of the quality of the trial solution represented by the string with respect to the function being optimised. Thus, high fitness corresponds to a high value (in a maximisation problem) or a low value (in a minimisation problem) of the function. If the upper and lower limits of the function being optimised are known, then absolute fitness may be used—where fitness values may be compared from generation to generation. Otherwise (as in most GA applications), dynamic fitness scaling can be adopted, where, in each generation the fitnesses of all the individuals are scaled relative to the best and worst members of the current population. Fitness is important in determining the likelihood of an individual taking part in crossover and also in deciding which individuals will survive into the next generation.

Selection refers to the way in which individual members of the population are chosen for subsequent crossover. There are a number of selection methods, the two most common being roulette wheel and tournament selection. In roulette wheel selection, a string is chosen at random and selected for crossover if its fitness value (f_i) is greater than a randomly generated number between 0 and 1 (i.e. if f_i > R[0, 1]), otherwise another string is chosen and tested. The analogy with a roulette wheel is apparent if one envisages a roulette wheel (see Fig. 3) with a slot for each member of the population, where each slot has a width (and therefore a probability of the ball dropping into it) which is proportional to the fitness of the corresponding member of the population. The tournament selection method picks a number of strings at random from the population to form a tournament pool. The two strings of highest fitness are then selected as parents from this tournament pool.

	Fig. 3 Schematic representation of roulette wheel selection. The probability of selecting an individual (i) for crossover is proportional to its fitness (fi)—i.e. proportional to the width of the corresponding slot on the roulette wheel.

Crossover is the way in which genetic information from two (or sometimes more than two) parent strings is combined to generate offspring. Fig. 4 shows two examples of commonly used crossover operators for a generic GA. In one-point crossover (Fig. 4a), two parent strings are cut at the same point and offspring are formed by combining complementary genes from the parents (i.e. the first part of parent 1 with the second part of parent 2 and vice versa). In two-point crossover (Fig. 4b), the two parents are cut at two points and offspring are formed by inserting the central sequence from parent 1 into parent 2 and vice versa. Other types of crossover are possible, such as uniform crossover, where offspring are generated by taking a certain number of genes from each parent, with no restriction on where these genes occur in the string.

	Fig. 4 Schematic representation of the crossover operation in a generic GA. (a) One-point crossover—two parents are cut at the same point and offspring are formed by combining complementary genes from the parents. (b) Two-point crossover—two parents are cut at the same two points.

Mutation—While the crossover operation leads to a mixing of genetic material in the offspring, no new genetic material is introduced, which can lead to lack of population diversity and eventually stagnation—where the population converges on the same, non-optimal solution. The GA mutation operator helps to increase population diversity by introducing new genetic material. As shown in Fig. 5, this can be accomplished by making a random change to one or more randomly chosen genes in an individual. In static mutation, the mutated gene is assigned a completely random value, while in dynamic mutation its value is changed by a small, random amount about its original value.

	Fig. 5 Schematic representation of the mutation operation in a generic GA. In this example a single variable (gene) has its value changed (represented here by its colour being changed.

Natural selection—In biological evolution the concept of the survival of the fittest (or best adapted to the environment) is a strong evolutionary driving force. In the case of a GA, although the selection is clearly not natural, individuals (be they parents, offspring or mutants) are likewise selected to survive into the next generation on the basis of their fitness (their quality with regards to the quantity being optimised). There are many modifications of the selection step, however: for example all mutants may be accepted, none of the parents may be accepted, or only the best parents may survive (the so-called elitist strategy—adopted so that the best member of the population cannot get worse from one generation to the next).

Schemata—The GA operators essentially exchange information between individual strings to evolve new and better solutions to the problem being optimised. A critical feature of the GA approach is that it operates effectively in a parallel manner, such that many different regions of parameter space are investigated simultaneously. Furthermore, information concerning different regions of parameter space is passed actively between the individual strings by the crossover procedure, thereby disseminating genetic information throughout the population. The GA is an intelligent search mechanism that is able to learn which regions of the search space represent good solutions, via the recognition of schemata. Fig. 6, shows three strings of seven variables (genes). The values of the variables (alleles—represented by blue symbols) are such that the three strings have identical values for genes 2 and 4. These three strings, therefore, contain a common schema (or building block, shown in red at the bottom of the figure) corresponding (in this example) to a star at position 2 and a snow flake at position 4. In a real application of a GA, a good schema (corresponding to a set of optimal or near-optimal variables) may be implicitly recognised and propagated within the population, if it leads to individuals of relatively high fitness.

	Fig. 6 Illustration of the concept of schemata (building blocks) in GA. The three strings depicted share the same values (alleles—blue symbols) of two of their variables (genes 2 and 4)—this constitutes a common schema (shown in red) for the three individuals.

The preceding section was an introduction to genetic operators for a generic GA. For a real application, however, both the crossover and mutation operators will generally have some problem-specific characteristics, to facilitate global optimisation in the particular problem being studied. Similarly, the extent to which schemata can be identified depends on how separable or correlated the variables are, which, in turn depends on the specific application. These considerations should become evident in the discussion (below) of the GA that we have developed for cluster geometry optimisation.

1.2.2 Applications of GAs in chemistry.

In this Perspective, there will be a brief discussion of critical stages in the development of GAs for studying the problem of cluster geometry optimisation, followed by a detailed treatment of our own cluster geometry optimisation program and its application to a number of cluster types. Finally, new and possible future directions for the application of GAs and related techniques in cluster science will be introduced.

The use of GAs for optimising cluster geometries was pioneered in the early 1990s by Hartke (for small silicon clusters)³³ and Xiao and Williams (for molecular clusters).³⁴ In both cases the cluster geometries were binary encoded, with the genetic operators acting in a bitwise fashion on the binary strings. Hartke has subsequently published the results of GA geometry optimisations for a number of different types of cluster, including silicon clusters,³⁵ water clusters³⁶ and mercury clusters.³⁷ An important stage in the evolution of GAs for cluster optimisation occurred when Zeiri³⁸ introduced a GA that operated on the real-valued cartesian coordinates of the clusters. This approach allowed for a representation of the cluster in terms of continuous variables and removed the requirement for encoding and decoding binary genes. The next signifficant step in the development of GAs for cluster optimisation was due to Deaven and Ho³⁹ who performed a gradient driven local minimisation of the cluster energy after each new cluster was generated. As Doye and Wales have pointed out,⁴⁰ the introduction of local minimisation effectively transforms the cluster potential energy hypersurface into a stepped surface, where each step corresponds to a basin of attraction of a local minimum on the potential energy surface, as shown in Fig. 7. This simplification of the surface greatly facilitates the search for the global minimum by reducing the space that the GA has to search. This principle also underpins the Basin Hopping MC method developed by Doye and Wales⁴⁰ and the MC + energy minimisation approach of Li and Scheraga.⁴¹ These related methods have proved very efficient for the structural optimisation of clusters, crystals and biomolecules.⁴² In the GA context, such local minimisation corresponds to Lamarckian, rather than Darwinian evolution, as individuals pass on a proportion of the characteristics that they have acquired to their offspring. In the case of clusters, these acquired characteristics are the geometries after local minimisation, rather than the characteristics they themselves inherited. Such hybrid (Lamarckian) GAs, which couple local minimisation with GA searching, have been found to improve GA efficiency for a number of different applications of GAs in global optimisation.⁴³

	Fig. 7 Simplification of the potential energy surface in the Lamarckian GA. Initially generated clusters which are in the same basin of attraction (e.g. A and A) are minimized to the same structure (A0) while cluster B minimizes to B0.

Another significant development in cluster optimisation GAs, also due to Deaven and Ho, was the introduction of the 3-dimensional cut and splice crossover operator.³⁹ This operator, which has been employed in most subsequent cluster GA work, gives a more physical meaning to the crossover process. In this crossover mechanism, good schemata correspond to regions of the parent clusters which have low energy local structure. Deaven and Ho applied their GA to the optimisation of carbon clusters bound by a tight binding potential³⁹ and Lennard-Jones clusters with 2–100 atoms.⁴⁴ Their work on Lennard-Jones clusters yielded many more low energy minima than had previously been found by alternative geometry optimisation methods.

Since the mid 1990s there have been a significant number of genetic algorithm programs which have been developed to study a wide range of cluster types. These GAs have introduced a number of new pseudo-genetic operators and ways of handling populations, carrying out crossover and mutation etc. A systematic discussion of this work (which is beyond the scope of the present Perspective) and further references can be found in a number of recent reviews^45–48 and a list of applications of GAs in the cluster field is included in Clark's bibliography of chemical applications of evolutionary algorithms.³²

A flow chart representing the operation of our cluster geometry optimisation GA program⁴⁹ is shown in Fig. 8. (Note: the term mating, used in this figure, can be regarded as synonymous with crossover.) Specific features of the GA are described below.

	Fig. 8 Flow chart for the Birmingham cluster genetic algorithm program.

For a given cluster size (nuclearity N),a number of clusters, (population size N_pop, typically ranging from 10 to 30) are generated at random to form the initial population (the zeroth generation). We have followed the approach of Zeiri³⁸ in using the real-valued cartesian coordinates of the cluster atoms as the genes. The x, y and z coordinates are chosen randomly in the range [0,N^1/3]. This ensures that the cluster volume scales correctly with cluster size (i.e. linearly with N). All of the clusters in the initial population are then relaxed into the nearest local minima, by minimising the cluster potential energy as a function of the cluster coordinates, using the quasi-Newton L-BFGS routine.⁵⁰ This routine utilises analytical first derivatives of the potential.

The GA operators of crossover, mutation and selection (on the basis of fitness) are performed to evolve one generation into the next.

As the cluster GA is being used to minimise the cluster potential energy (V_clus, a negative quantity), the lowest energy (most negative V_clus) clusters have the highest fitness and the highest energy (least negative V_clus) clusters have the lowest fitness. The cluster GA uses dynamic fitness scaling: the fitness of the lowest energy cluster in the current population is equal to (or close to) one and the highest energy cluster has a fitness equal to (or close to) zero. Dynamic scaling is achieved by using a normalised value of the energy, , in the fitness calculations:

i = (Vi – Vmin)/(Vmax – Vmin)

(1)

The most common fitness functions that we have used are:

Exponential:

fi = exp(–i)

(2)

Linear:

fi = 1 – 0.7i

(3)

Hyperbolic tangent:

fi = ½[1 – tanh(2i – 1)]

(4)

The choice of fitness function controls how rapidly fitness falls off with increasing cluster energy—i.e. the relative weighting given to good vs. mediocre and bad structures.

The selection of parents is accomplished using either roulette wheel or tournament selection. In both of these selection schemes, low energy clusters (with high fitness values) are more likely to be selected for crossover and therefore to pass their structural characteristics on to the next generation. Once a pair of parents have been selected, they are subjected to the crossover operation.

Crossover is carried out using a variant of the cut and splice crossover operator of Deaven and Ho,³⁹ as shown schematically in Fig. 9. In our implementation, random rotations (about two perpendicular axes) are performed on both parent clusters and then both clusters are cut horizontally about one or two positions (corresponding to one-point and two-point crossover, respectively) parallel to the xy plane, and complementary fragments are spliced together. For the single cut method, the cutting plane can be chosen at random, it can be defined to pass through the middle of the cluster, or weighted according to the relative fitnesses of the two parents (with more atoms taken from the fittest parent). For the double cut method, the cutting planes are chosen at random. (Note that in our GA only one offspring is produced from each crossover event, though, in principle two offspring can be produced.)

	Fig. 9 Schematic representation of the Deaven–Ho cut and splice crossover operation, as implemented in our GA program. (Single cut crossover is shown.)

Crossover continues until a predetermined number of offspring (N_off) have been generated. The number of offspring is generally set to approximately 80% of the population size. Each offspring cluster is subsequently relaxed into the nearest local minimum, as described above. The local minimisation step obviously changes the structure of the offspring cluster, and this structural rearrangement will be greatest in the region of the join between the two fragments donated by its parents. As the clusters get larger, however, the perturbation due to the local minimisation should become relatively smaller and confined to the join region. In this way, the principle of schemata should apply, as parents with high fitness are more likely to have fit offspring by passing on fragments with low energy arrangements of atoms.

In an attempt to avoid stagnation and to maintain population diversity, a mutation operator is introduced, whereby each individual has a probability (P_mut) of undergoing mutation. The mutation perturbs some or all of the atoms (or molecules) within the cluster.

We have adopted a number of mutation schemes, with the choice as to which scheme(s) to use depending on the type of cluster being studied:

Atom displacement—The cluster is mutated by replacing the atomic coordinates of a certain number of the atoms with randomly generated values.

Twisting—The cluster is mutated by rotating the upper half of the cluster about the z axis by a random angle relative to the bottom half.

Cluster replacement—An entire cluster is removed and replaced by a new, randomly generated cluster. The cluster is generated in an identical way to that used for the generation of the initial population.

Atom permutation—The cluster is mutated by swapping the atom types of one or more pairs of atoms without perturbing the structure of the cluster. This type of mutation is used for hetero-elemental clusters, such as ionic clusters and bimetallic clusters.

After mutation, each mutant cluster is subsequently relaxed into the nearest local minimum, using the L-BFGS minimisation routine, as described above.

The program contains an option of removing clusters from the population that have a difference in energy less than a value E (typically 10^–6 eV). If two or more clusters have energies less than E apart, then the lowest energy cluster is retained and the other clusters are discarded. Use of this operator ensures that population diversity is maintained.

The new population (the next generation) is selected from the N_pop lowest energy (highest fitness) clusters selected from the set containing the old population the new offspring clusters and the mutated clusters. The inclusion of clusters from the previous generation makes the GA elitist, ensuring that the best member of the population cannot get worse from one generation to the next.

Once the new generation has been formed, the potential energies of the best (V_min,) and worst (V_max) members of the population are recorded and the fitness values calculated for the entire population. The whole process of crossover, mutation and selection is then repeated for a specified number (N_gen) of generations or until the population is deemed to have converged. The population is considered to be converged if the range of cluster energies in the population has not changed for a prescribed number of generations.

A considerable amount of effort has been expended in optimising the GA operations and parameters described above.⁵¹ It should be noted that, because of the stochastic nature of the GA, the GA program is run several times for each cluster nuclearity and for each set of operations/parameters. It is generally found that there is not a great dependence of the success rate (of finding the global minimum structure) on the type of fitness function used, though even small improvements are useful, especially for larger cluster sizes. The best types of crossover and mutation operators to use can be system dependent (e.g. using atom permutation mutation in heteronuclear clusters) and the absolute and relative rates of crossover and mutation can be critical in determining the frequency with which the GM is found and the average number of generations required. Similarly, for larger clusters, a larger population size and maximum number of generations is generally required. While a detailed study of the effect of the GA operators and the values of the GA parameters is beyond the scope of this Perspective, the interested reader can find such information in our previous publications (referenced below).

In the next section, details of the application of our GA to a number of cluster types will be presented. Values of the GA parameters, and any other options adopted are listed where appropriate.

Our GA (or modifications of it) has been applied to the study of a wide variety of clusters, ranging from model Morse clusters⁵² to fullerenes,^53,54 ionic clusters,⁵⁵ water clusters,⁵⁶ metal clusters^57,58 and bimetallic nanoalloy clusters^59–63 Further details are given in these papers and in a recent review.⁴⁵ Many of the structures of the global minima found by our GA, for a number of cluster types, can be found on the Birmingham Cluster Web site.⁶⁴

The first type of potential that was studied with our cluster GA was the Morse potential,⁵² because, although Morse clusters had not previously been studied with GAs, Doye and Wales had described the structural consequences of varying the range of the Morse potential.¹³ Furthermore, using the Basin Hopping Monte Carlo approach,⁴⁰ they had found global minima for Morse clusters with different range parameters and noted that the short range Morse potential (which has many local minima and a very noisy potential energy surface) presents a particular challenge for global optimisation techniques.¹³ We, therefore, decided to apply our GA to find global minima for Morse clusters and to compare our results with Doye and Wales' tabulated coordinates and energies, which are available on the Cambridge Cluster Database website.⁶⁵

The Morse potential is a pair-wise additive potential which depends only on the separations r_ij between pairs of atoms:

VMij = De[e(l – rij/re)(e(l – rij/re) – 2)]

(5)

As in the work of Doye and Wales,¹³ we have adopted a simplified, scaled version of the Morse potential with D_e and r_e both set to one:

VMij = De[e(l – rij)(e(l – rij) – 2)]

(6)

The total potential energy of a cluster of N atoms is then obtained by summing over all atom pairs:

(7)

We have compared short range ( = 14) and medium range ( = 6) Morse potentials. Firstly, the efficiency of the GA was compared with that of a simple random search algorithm (RSA) for Morse clusters with N = 20, 30, 40 and 50 atoms. The RSA generates a number of structures using an identical method to that used to generate the initial population of the GA and then minimises the energy of these clusters using the L-BFGS local minimisation routine. Table 1 lists the success or failure of the RSA in finding the global minimum energy structures of Morse clusters bound with a Morse potential with = 6. It can be seen that the RSA is only successful in locating the GM for clusters of 20 atoms, even when 30,000 random geometries are minimised, as for the 50-atom clusters.
Table 1 Comparison of the success of the random search algorithm and the genetic algorithm in locating the global minima of Morse (= 6) clusters with N= 20, 30, 40 and 50 atoms. (N_search is the number of minimisations carried out by the RSA and <N_min> is the average number of minimisations required by the GA to find the global minimum)

N	Nsearch	GM found?	Nmin
20	5000	Yes	31
30	10000	No	301
40	20000	No	323
50	30000	No	472

The GA was used to search for the GM of Morse clusters of the same nuclearities, with the GA being run from 40 different initial populations for each nuclearity. The GM was found during each run of the GA for all 4 cluster nuclearities, which is a considerable improvement over the RSA where the GM was found only for the smallest cluster. The average number of energy minimisations N_min required to locate the GM for these clusters, measured over the 40 runs for each nuclearity, are listed in Table 1. The GA requires an average of 472 energy minimisations to locate the GM for the 50-atom Morse cluster, which is significantly less than the 30,000 energy minimisations performed by the unsuccessful RSA algorithm, so the success of the GA, as compared with the random search algorithm, is obvious. Table 1 also reveals that it is particularly challenging to find the GM for N = 30, which takes nearly as many minimisations (on average) as for N = 40.

The cluster geometry optimisation GA program was subsequently used to study Morse clusters (for = 6 and 14) with N = 19–50 atoms. Preliminary calculations were performed on a number of trial clusters and the following optimal values were obtained for the GA program parameters: N_pop = 10; N_off = 8; P_mut = 0.1; and N_gen = 10–300 (increasing with cluster size). The exponential fitness function was adopted, with roulette wheel selection, random single cut crossover and atom replacement mutation.

The GA located all of the previously published GM¹³ for Morse clusters with 19 to 50 atoms, both for medium ( = 6) and short range ( = 14) Morse potentials. The GA found a lower energy structure for the N = 30 ( = 14) Morse cluster than was initially reported by Doye and Wales,¹³ though this structure has indeed been found by their Basin Hopping algorithm.⁶⁵

The structures of some of the GM for = 6 and = 14 are shown in Fig. 10. These structures have been discussed in detail by Doye and Wales,¹³ the most obvious difference between the two potentials being that the longer range ( = 6) potential tends to favour poly-tetrahedral, icosahedral geometries, while the short range ( = 14) potential favours decahedral and fcc-like packing (such as the truncated octahedral cluster which is the global minimum for N = 38).

	Fig. 10 GM for Morse clusters (MN) with 19, 38 and 50 atoms for = 6 and = 14.

One way of monitoring the progress of the cluster GA is to construct an Evolutionary Progress Plot (EPP),⁶⁶ of the lowest energy (V_min,), highest energy (V_max) and the average energy (V_ave), as a function of generation number. (As well as EPPs, a number of other methods for analyzing the progress of GAs have been developed, such as the use of family trees, where the ancestry of a particular individual can be mapped out.⁶⁷) Such EPPs are shown in Fig. 11 for N = 38, with = 6 and = 14. The EPPs show that there is a rapid improvement in the population (a sharp drop in V_min, V_max and V_ave) in the early generations, relative to the randomly generated initial population (generation 0). This early improvement is mainly due to the crossover process. Subsequent, less dramatic improvement occurs in a stepwise fashion and may be due to crossover or mutation.

	Fig. 11 Evolutionary Progress Plots for 38-atom Morse clusters with = 6 and = 14.

Fig. 11 shows that, in both cases, the population converges on a single structure—as evidenced by V_min, V_max and V_ave becoming equal (diversity checking was not used in this case). The converged structure was found to be the global minimum in each case. The fcc-like truncated octahedral geometry of the 38-atom ( = 14) Morse cluster is difficult to find with most global optimisation techniques,¹³ but is here found before the 100th generation—even for a small population size of 10. Comparison with EPPs for larger clusters confirms that the GA must, in general, be run for a greater number of generations for larger cluster sizes. Similarly, comparison of the EPPs for = 14 with those for = 6 confirms that the short range Morse potential is more difficult to search,¹³ typically taking 2–3 times as many generations to find the GM.

The bonding in ionic clusters (clusters of NaCl, MgO etc.) is dominated by electrostatic interactions and the simplicity of modelling ionic clusters makes them ideal systems in which to study size-dependent properties. Ionic clusters have a number of practical applications, such as: silver halides in the photographic process; ZnS clusters as gas sensors; CdS and CdSe clusters as photodetectors; and polar semiconductors made from GaAs and InP clusters. The study of NaCl and NaBr clusters is also important in determining the mechanisms of ozone depletion and pollution in marine environments. This section concerns the simplest ionic oxide species, magnesium oxide, which crystallises in the rock salt structure. The suitability of the GA for investigating the structure of ionic clusters was investigated using a rigid ion potential to describe the electrostatic bonding. The study of finite MgO clusters may lead to an improved understanding of (for example) epitaxial growth of MgO surfaces. MgO clusters are also of interest in their own right: for example they have been invoked as possible nucleation sites for the formation of particulate matter in oxygen-rich regions of space.⁶⁸

The rigid ion potential is a simple model of the bonding in ionic solids and their clusters, in which the ions are of fixed size and carry fixed charges. The interaction between a pair of ions is given by the sum of the long range electrostatic Coulomb energy and a repulsive Born–Mayer potential (which reflflects the short range repulsive energy due to the overlap of the electron density of the ions):

Vij = (qiqje2/40rij) + Bijexp(–rij/ij)

(8)

_ij

_ij

Although the potential parameters (listed in Table 2) were derived by Lewis and Catlow for formal charges of Mg²⁺ and O^2–,⁶⁹ we have used the cluster GA to optimise the geometries of stoichiometric (MgO)_N clusters (N = 10–35), with formal ionic charges of Mg⁺O^–, as well as Mg²⁺O^2–.⁵⁵ The following GA parameters were used: N_pop = 20; roulette wheel selection; tanh fitness function; weighted single cut crossover, N_off = 16; twist mutation, P_mut = 0.05; N_gen = 200.
Table 2 Rigid ion potential parameters for MgO

Parameter	Value
Bij (Mg–O)	821.6 eV
Bij (O–O)	22764 eV
ij (Mg–O)	0.3242
ij (O–O)	0.1490

Using the GA and the rigid ion potential, the lowest energy Mg⁺O^– clusters are found to be compact cubic structures based on the NaCl structure (the structure of bulk MgO), whereas the Mg²⁺O^2– clusters are more spherical, cage like structures with square, hexagonal and octagonal faces. The cation–cation and anion–anion repulsions are much larger in the Mg^q+O^q– clusters with q = 2, as compared with those with q = 1. This causes the structures to open out from highly coordinated lattice structures (q = 1) to cage like structures (q = 2). Ziemann and Castleman have studied (MgO)⁺_N and (MgO)_NMg⁺ clusters experimentally, using laser-ionisation time of flight mass spectrometry.⁷⁰ The mass spectra showed magic number clusters (corresponding to relatively intense peaks) for compact cuboidal (rock salt) structures, as predicted by the rigid ion potential with ionic charges of +1 and –1. Similar structures are known to be preferred for alkali halide clusters.⁷¹

The rigid ion potential does not contain any terms to describe ionic polarisation. As the O^2– ions are highly polarisable and the Mg²⁺ ions are strongly polarising, the rigid ion potential, therefore, provides a poor description of Mg²⁺O^2– clusters. The effect of polarisation in the clusters will be to reduce the effective charges of both ions, thereby reducing the repulsion felt between like ions. This will lead to the more highly coordinated (bulk-like) cuboidal structures being preferred over the cage like structures, due to the higher degree of bonding present. The rigid ion potential calculations used singly charged ions in order to mimic, partially the effect of polarisability. It should be noted that ab initio calculations^72–74 and calculations using model potentials including explicit ion polarisabilities⁷⁵ confirm that MgO clusters should adopt rock salt-like structures.

As noted above, MgO clusters with formal ionic charges of ±1 and ±2 have different structures. To investigate the dependence of structure on charge further, we have studied the effect, on the lowest energy structures of small (Mg^q+O^q–)_N clusters, of varying the magnitude of the charge (q) in increments of 0.1 from 1.0 to 2.0. The results for N = 8 and 9 are summarised below.

Fig. 12a shows that the (MgO)₈ cluster has three possible structures, a 4 × 2 × 2 lattice for q = 1.0–1.3, a cage structure for q = 1.4–1.6 and a structure consisting of two stacked octagonal rings for q = 1.7–2.0. The two more open structures can both be generated by lengthening different Mg–O contacts in the 4 × 2 × 2 cuboid. The stacked ring structure is the ground state structure predicted by the Hartree–Fock (HF) calculations of de la Puente and Malliavian, while their correlated HF (CHF) calculations (which include intra-ion correlation) predict a structure of two hexagonal MgO rings with an Mg₂O₂ square capping one of the square faces.⁷³ In contrast to the N = 8 cluster, a single geometry was found for (MgO)₉ in the range q = 1.0–2.0: corresponding to three stacked hexagonal Mg₃0₃ rings. This geometry, which is shown in Fig. 12b, is the ground state structure found by de la Puente and Malliavian using both the HF and CHF methods.⁷³

	Fig. 12 (a) The structures of (Mgq+Oq–)8 clusters with q in the range 1.0–2.0. (b) The structure of (Mgq+Oq–)9 clusters with q in the range 1.0–2.0. (Mg cations are shown in red and O anions in yellow.)

Examples of the cage structures predicted for q = 2 include the pair of 60 ion enantiomers, shown in Fig. 13, which are predicted as the GM for (MgO)₃₀, for q = 2.

	Fig. 13 Enantiomers of (Mg2+O2–)30.

Although these pseudo-fullerene hollow structures, comprised of 4-, 6- and 8 membered rings are inconsistent with experimental and ab initio theoretical studies of (MgO)_N clusters, it is possible that such structures will be stabilised for polar III–V semiconductor nanoparticles, where the formal ionic charges are likely to be higher, even taking polarisation into account. DFT calculations by Sedlmayr and co-workers for (A1N)_N, (GaAs)_N and (InSb)_N clusters (N = 12, 24, 60) confirm that such structures should have high energetic stabilities.⁷⁶

EPPs showing the convergence of the GA on the lowest energy (MgO)₃₀ clusters, with formal charges of ±1 and ±2 are shown in Fig. 14. The GA required 74 generations to find the lowest energy structure of (Mg²⁺O^2–)₃₀ but only 7 generations to find the GM for (Mg⁺O^–)₃₀. In both cases the minimum, maximumum and average energies converge, indicating that the GA has converged on a single solution. The GA converges in 27 generations for q = 1 and in 100 generations for q = 2. These results are typical of the other cluster nuclearities studied. The rigid ion potential with formal charges of ±2 is shorter ranged and the potential energy surface is therefore likely to have more local minima for the GA to search, leading to greater difficulty in finding the global minimum.

	Fig. 14 EPPs for (Mg+O–)30 and (Mg2+O2–)30 clusters.

There is continuing interest in metal clusters because of potential applications in fields such as catalysis and nano-electronics (e.g. in single electron tunnelling devices).¹ It is known that alkali metal clusters, with sizes of up to thousands of atoms, conform to the jellium model, in that certain nuclearities are relatively stable (the so-called magic numbers) due to their having filled electronic shells.⁷⁷ By contrast, clusters of transition metals and some main group metals (e.g. Al, Ca and Sr) exhibit magic numbers which correspond to clusters consisting of concentric polyhedral shells (geometric shells) of atoms, where the relative stability of a given cluster is determined by the competition between packing and surface energy effects.⁷⁸

The range of properties of metallic systems can be greatly extended by taking mixtures of elements to generate intermetallic compounds and alloys. In many cases, there is an enhancement in specific properties upon alloying, due to synergistic effects and the rich diversity of compositions, structures and properties of metallic alloys has led to widespread applications in electronics, engineering and catalysis. The desire to fabricate materials with well defined, controllable properties and structures, on the nanometre scale, coupled with the flexibility afforded by intermetallic materials, has generated interest in bimetallic alloy clusters—or nanoalloys.^79–81 One of the major reasons for interest in nanoalloy particles is the fact that their chemical and physical properties may be tuned by varying the composition and atomic ordering, as well as the size of the clusters. Their surface structures, compositions and segregation properties^82,83 are of interest as they are important in determining chemical reactivity, and especially catalytic activity.^84,85 Nanoalloy clusters are also of interest as they may display structures and properties which are distinct from those of the pure elemental clusters. There are also examples of pairs of elements (such as iron and silver) which are immiscible in the bulk phase but which readily mix in finite clusters.⁸⁶

There have been a number of MD and MC studies of alloy clusters using empirical many-body potentials.^12,79,87 To date, however, there have been few studies of nanoalloys using electronic structure methods. Some notable applications of electronic structure calculations to small nanoalloy clusters include the work of Fortunelli and Velasco (Extended Hückel: Pt–Fe clusters⁸⁸), Calleja et al. (DFT: Ni–Al clusters⁸⁹) and Bromley et al. (DFT: Ru–Cu clusters⁹⁰). Where comparisons are possible between DFT calculations and empirical potentials, it is often found that the empirical potentials cannot reproduce specific electronic effects such as Jahn–Teller distortions, which are observed for small clusters. For larger clusters, DFT studies on homonuclear clusters (for example Pd clusters⁹¹) show that the lowest energy structures found using empirical potentials, while not always corresponding to DFT global minima, may lie close in conformation space to the DFT global minima, or may correspond to low-lying metastable isomers. In future, it should prove possible to combine empirical potentials with electronic structure calculations in order to explore conformation space more fully—i.e. to select regions for study at the higher level of theory (see Section 5.4). Such an approach has indeed been applied by Bromley et al. to propose a structure for a Ru₁₂Cu₄ cluster formed by thermal decomposition of an organometallic precursor cluster in mesoporous silicon.⁹⁰ Alternatively, the results of DFT calculations on selected benchmark clusters can be used to test the parameters of an empirical potential (as in a recent Gupta potential-DFT study of Co–Cu nanoalloy clusters by Wang et al.⁹²) or to facilitate reparameterisation.

A brief review is presented here of the application of our cluster GA to study Cu–Au,⁵⁹ Ni–Al⁶¹ and Pd–Pt clusters.⁶² In these studies, metal–metal bonding is modelled by the many-body Gupta potential,⁹³ which can be expressed as the sum of repulsive pair (V^r) and attractive many-body (V^m) terms:

(9)

(10)

(11)

In the above equations, r_ij is the distance between atoms i and j in the cluster and the parameters A, r₀, , p and q are fitted to experimental values of the cohesive energy, lattice parameters and independent elastic constants for the bulk metals and alloys at 0 K.⁹³ For A–B alloy clusters, the parameters take different values for each of the different types (A–A, B–B and A–B) of interaction. (Note: in the above equations, a and b are the atom types of atoms i and j, respectively.)

Jellinek has introduced the term homotops to describe A_aB_b nanoalloy isomers, for fixed number of atoms (N = a + b) and composition (a/b ratio), which have the same geometrical arrangement of atoms, but differ in the way in which the A and B-type atoms are arranged⁷⁹ (i.e. they are related by atom label permutation). As the number of homotops rises combinatorially with cluster size, global optimisation (in terms of both geometrical isomers and homotops) is a difficult task.

In the calculations reported below for Cu–Au and Pd–Pt clusters, the GA operators included atom permutation mutation for the alloy clusters—which was found to greatly improve the reproducibility of the results and the likelihood of finding the GM. The following GA input parameters were used: N_pop = 30; N_gen = 400; N_off = 24; P_mut = 0.1.

4.3.1 Cu–Au nanoalloys.

We have made a study of stoichiometric clusters with the compositions of the common bulk Cu–Au alloy phases: (CuAu₃)_N; (CuAu)_N; and (Cu₃Au)_N and have compared them with the pure copper and gold clusters.⁵⁹ The Gupta parameters used to study these clusters are listed in Table 3.⁹³
Table 3 Gupta potential parameters for Cu–Au nanoalloys

---

Parameter	Cu–Cu	Cu–Au	Au–Au
A/eV	0.0855	0.1539	0.2061
/eV	1.2240	1.5605	1.7900
p	10.960	11.050	10.229
q	2.2780	3.0475	4.0360
r0/	2.556	2.556	2.884

---

Using our GA, pure copper clusters were found to adopt regular, symmetric structures based on icosahedral packing, while gold clusters showed a greater tendency towards amorphous structures, as found previously by Garzón et al.⁹⁴ In many cases (e.g. for 14, 16 and 55 atoms), the replacement of a single Au atom by Cu was found to change the GM structure to that of the pure Cu cluster, which has also been observed by López et al.⁸⁷

As an example of the results obtained for Cu–Au nanoalloys, the global minima found by the GA for 40-atom clusters of varying composition—Cu₄₀; Cu₃₀Au₁₀; Cu₂₀Au₂₀; Cu₁₀Au₃₀; and Au₄₀—are shown in Fig. 15. From the figure, it is apparent that the clusters Cu₄₀, Cu₃₀Au₁₀ and Cu₂₀Au₂₀ have the same decahedral geometry. In the alloy clusters the Au atoms generally lie on the surface, while the Cu atoms are encapsulated. Although the Au₄₀ cluster has a low symmetry, amorphous structure, the gold-rich Cu₁₀Au₃₀ cluster has a structure which is more symmetrical than the Au₄₀ GM, but it is not decahedral. The Cu₁₀Au₃₀ GM has a flattened (oblate) topology, in which all of the Au atoms, except one, lie on the surface of the cluster and 7 of the Cu atoms occupy interior sites. The observed atomic segregation in the alloy clusters is favoured because of the lower surface energy of Au and the shorter Cu–Cu bonds which can allieviate the internal strain of the cluster. Similar tendencies have been found in a study of larger Cu–Au clusters.⁹⁵

	Fig. 15 GM structures of selected 40-atom Cu (red), Au (yellow) and Cu–Au clusters.

4.3.2 Ni–Al nanoalloys.

The cluster optimisation GA has been used to find the GM for Ni–Al nanoalloys with the approximate composition of the bulk alloy Ni₃A1.⁶¹ The Gupta potential parameters for the Ni–Al system are listed in Table 4.⁹³
Table 4 Gupta potential parameters for Ni–Al nanoalloys

---

Parameter	Ni–Ni	Ni–Al	Al–Al
A/eV	0.0376	0.0563	0.1221
/eV	1.070	1.2349	1.316
p	16.999	14.997	8.612
q	1.189	1.2823	2.516
r0/	2.4911	2.5222	2.8637

---

An important aspect of a reliable search method is that it enables one to have confidence when comparing stuctures and considering trends in structures and homotop stability as a function of size and composition. Fig. 16 shows, for example, that the GM structure type changes from truncated octahedral (fcc packing) to pseudo-icosahedral on adding a single Ni atom, going from 38-atom Ni₂₈A1₁₀ to 39-atom Ni₂₉A1₁₀. This icosahedral growth continues until the complete 2-shell icosahedral 55-atom cluster Ni₄₁A1₁₄ (also shown in Fig. 16) is reached.

	Fig. 16 GM for Ni28A110, Ni29A110 and Ni41A114. (Ni atoms are shown in red and Al in green.)

For Ni–Al clusters in general, we have found that the Al atoms lie preferentially on the surface of the cluster. In the Ni-rich clusters the Al atoms tend to be separated (i.e. there are few Al–Al bonds), as can be seen in Fig. 16. These observations, which are consistent with the empirical potential calculations of Rey et al.⁹⁶ and Rexer et al.,⁹⁷ can be rationalised by the lower surface energy of Al, the smaller size of Ni and the greater strength of Ni–Al bonding, as compared with Ni–Ni and Al–Al interactions.⁶¹

4.3.3 Pd–Pt nanoalloys.

A detailed study has been made of Pd, Pt and Pd–Pt bimetallic clusters, (PdPt)_M, with up to 56 atoms, modelled by the many-body Gupta potential.⁶² The Gupta parameters used to study these clusters are listed in Table 5.⁹³
Table 5 Gupta potential parameters for Pd–Pt nanoalloys. (See text for discussion of Pd–Pt parameter sets I, II and III)

---

Parameter	Pd–Pd	Pt–Pt	Pd–Pt(I)	Pd–Pt(II)	Pd–Pt(III)
A/eV	0.1746	0.2975	0.23	0.35	0.23
/eV	1.718	2.695	2.2	2.2	3.0
p	10.867	10.612	10.74	10.74	10.74
q	3.742	4.004	3.87	3.87	3.87
r0/	2.7485	2.7747	2.76	2.76	2.76

---

Pure clusters of Pd and Pt are predicted to adopt a variety of structures, depending on the cluster size. Many of the structures are regular (ordered), though there is a tendency (which is greater for Pt than for Pd) towards forming disordered structures. Another difference between the two metals is the increased tendency for Pd to form clusters based on icosahedral packing.

Heteronuclear Pd–Pt interaction parameters have not previously been derived for the Gupta potential. As the stoichiometric 1 : 1 Pt–Pd alloy is a solid solution (Pd_0.5Pt_0.5), with a small exothermic enthalpy of formation (–4 kJ mol^–1),⁹⁸ rather than an ordered intermetallic, the Pd–Pt parameters were generated by taking averages of the Pd–Pd and Pt–Pt parameters. These parameters are listed in Table 5 as Pd–Pt(I) and were used for most of the calculations on Pd–Pt clusters.

Using parameter set I, it was found that the lowest energy structures for stoichiometric (PtPd)_M nanoalloy clusters generally have different geometric structures than the corresponding pure Pd or Pt clusters: with a reduced tendency to display icosahedral packing and a larger number of capped decahedral structures. Compared with Pd, there is also an increase in the number of disordered structures for the Pd–Pt clusters. Shell-like atomic segregation is favoured for these Pt–Pd clusters, with the surface becoming richer in Pd and the core becoming richer in Pt. This segregation is consistent with experimental studies on Pd–Pt particles^99,100 and can be explained in terms of the lower surface energy of Pd and the greater cohesive energy of Pt.⁶²

An example of a nuclearity for which pure Pd and Pt clusters have different geometries is found for 14 atoms. Pd₁₄ has a capped icosahedral structure, while Pt₁₄ is a bicapped hexagonal antiprism. Fig. 17 shows that doping a single Pd atom into Pt₁₄, changes the structure to a capped icosahedron (where the Pd atom adopts the low-coordinate capping site), whereas doping a Pt atom into Pd₁₄ leaves the structure unchanged (though the Pt atom adopts the highly-coordinated interstitial site). These findings are consistent with our results for Cu–Au clusters.⁵⁹ For intermediate compositions the GM is usually the capped icosahedron, with the exception of Pd₁₀Pt₄–Pd₆Pt₈, where alternative nonicosahedral structures are found—as shown in Fig. 17 for Pd₇Pt₇.

	Fig. 17 GM structures for selected 14-atom Pd (blue), Pt (green) and Pd–Pt clusters, using Pd–Pt parameter set I.

While the averaged Pd–Pt parameters (set I) give a good qualitative description of Pd–Pt clusters, we have also investigated the effect of varying the Pd–Pt parameters on the structures of a limited set of Pd–Pt clusters, with the Pd : Pt ratio fixed at 1 : 1. The modified parameter sets were obtained from the averaged set I: by increasing the 2-body energy scaling parameter (A), such that A(Pd–Pt) > A(Pt–Pt) > A(Pd–Pd) (denoted Pd–Pt(II) in Table 5); and by increasing the many-body energy scaling parameter (), such that (Pd–Pt) > (Pt–Pt) > (Pd–Pd) (denoted Pd–Pt(III) in Table 5). In both cases, all other parameters were left unchanged.

Fig. 18 compares the lowest energy isomers obtained for (PdPt)₁₄ clusters for parameter sets I, II and III. Cluster I has a capped decahedral structure, based on a central 15-atom stack of two pentagonal prisms, with 2 interstitial atoms, 2 caps on the 5-fold axis and 9 of the 10 square faces capped (giving a star-shaped cross section). In agreement with our general findings for Pd–Pt clusters, using parameter set I, all 7 Pd atoms lie on the surface of the cluster, showing a preference for the low-coordinate capping sites. For parameter set II, Pd–Pt interactions are destabilised relative to Pt–Pt and Pd–Pd, so that the Pd–Pt clusters segregate into Pd and Pt sub-units held together by a small number of Pd–Pt bonds. In cluster II, both the Pd₁₄ and Pt₁₄ sub-units have the bicapped hexagonal antiprismatic structure previously observed for Pt₁₄, even though the preferred structure for Pd₁₄ is the capped icosahedron (see Fig. 17). Presumably this alternative structure is adopted because it enables stronger inter-fragment bonding. For parameter set III, ordered mixing of the Pd and Pt atoms is favoured, so as to maximise the number of strong Pd–Pt bonds. The clusters adopt body centred cubic (bcc) packing (the -brass structure), with each Pt atom surrounded by a cube of Pd atoms and vice versa, as shown for cluster III. Of the three parameter sets tested here, set I clearly yields structures and segregation patterns which are in closest agreement with experiment.^99,100

	Fig. 18 Comparison of the GM for (PdPt)14 clusters, for Pd–Pt parameter sets I, II and III.

4.3.4 Atomic mixing in nanoalloy clusters.

Maximisation of the number of the strongest interatomic interactions—this may favour mixing or segregation, depending on the system.

Minimisation of the cluster surface energy—this favours segregation, with the cluster surface becoming richer in the element which has the lower surface energy.

Minimisation of bulk strain—this favours the location of the smaller atom at the centre of icosahedral clusters, for example.

Three main types of mixing patterns of nanoalloy clusters can be identified, which are shown in Fig. 19, with the observed atomic arrangement for a particular A–B system depending critically on the balance of the above factors. Firstly, there are structures (Fig. 19a) which can be described ideally as shell-segregated, with a shell of one type of atoms (A) surrounding a core of another (B)—though there may be some mixing between the shells. Secondly, there are segregated structures (Fig. 19b) consisting of A and B sub-clusters, which may share a mixed interface (top) or which may only have a small number of A–B bonds (bottom). Finally, there are mixed A–B clusters (Fig. 19c), which may either be random (i.e. a solid solution—top) or ordered (bottom). Other intermediate mixing arrangements are possible (for example layered structures). Cluster generation under kinetic growth conditions, for example, can lead to metastable mixing patterns, such as onion-like alternating –A–B–A– shells.¹⁰³

	Fig. 19 Schematic representation of idealised atomic mixing/segregation found for nanoalloy clusters. (a) Shell segregation. (b) Sub-cluster segregation (top = with some mixing, bottom = complete segregation). (c) Mixed (top = random mixing, bottom = ordered mixing). (The two elemental components of the nanoalloy are denoted by blue and yellow colouring.)

As far as I am aware, to date there have been no electronic structure calculations on nanoalloy clusters for the systems we have studied in the size ranges reported here. It is to be hoped that our results (more details and references to lists of cluster coordinates can be found in our previous publications) will provide a useful starting point for DFT (or other electronic structure) calculations of Ni–Al, Cu–Au and Pd–Pt nanoalloys.

There have been a number of important recent developments and advances in the GA methodology, which have already been applied to the problem of cluster geometry optimisation—for example the introduction of predators, parallel GAs, hybrid GAs and GAs combining ab initio and empirical calculations.^45,36 Some of these approaches are briefly discussed below.

GAs have been very successful in determining global minima, but in a number of physical applications, structures corresponding to higher local minima may be more important than the GM. For example, carbon cluster ions formed in laser-ablation experiments¹⁰⁴ are observed in several different geometries, distinguished by their relative mobilities. In some experiments, kinetically favoured, higher energy isomers may be formed, and the distribution of and interconversion between isomers is also of great interest.¹³ Finally, it is worth noting that the biologically active forms of proteins may not always correspond to global minima.¹⁰⁵

In recent work, we have extended the analogy with natural evolution by considering the use of predators to remove unwanted individuals or traits from the population.⁵³

Sometimes unwanted members of a population can be removed by imposing a constraint on the fitness function, however, in seeking minima other than the GM, a suitable modification of the fitness function is not always possible. In principle, predation can be carried out using any property of the cluster, for example a shape selective predator could be used to remove individuals from the population (with a certain probability) if they show (or fail to show) certain topological features, such as sphericity, ring size or adjacent/non-adjacent pentagons (in the case of fullerene clusters).

The simplest predator is the energy predator, in which clusters with energies at or below a certain value are removed (predated) from the population. The energy predator can thus be used to search for low energy minima other than the GM, or to enhance the efficiency of the GA by removing specific low energy non-global minima that the GA may be incorrectly converging towards. When searching for low energy isomers, the GA is first run without the predator to find the GM cluster and its energy. Then the predator is invoked to remove the GM, so that the GA converges on the first metastable isomer. This is then repeated, to remove both the GM and the first metastable isomer, and the cycle is continued until the required number of low-lying isomers have been found. So far, we have applied this energy predator to find low-lying metastable isomers of aluminium clusters⁵³ and carbon fullerenes^45,54 and to explore the stability of the evolutionary trajectory when removing certain intermediate structures, which may either be key ancestors or evolutionary dead ends on the path to finding the global minimum.⁶³

The genetic algorithm has an intrinsically parallel nature: the production of each new cluster during a generation is independent of the production of any other new cluster. In this way the GA can study different regions of the potential energy surface in parallel. This parallel nature makes it easy to create a parallel version of the GA code that will utilise many processors simultaneously, thereby spreading the computational effort and reducing the time taken for the run of the algorithm.⁴³ The availability of bigger, faster and cheaper multi-processor computers, Beowulf clusters etc. means that parallel programming is an important and growing area in modern computational science.

There are a number of parallel GA paradigms which may be applied—for example master-slave, distributed and sub-population models.^45,51 The master-slave and distributed parallel algorithms work analogously to the normal, serial GA but they spread the work (carrying out the genetic operations) over a number of processors. They differ in whether each processor holds a copy of the entire population (distributed) or a part of the population (master-slave). The sub-population algorithm, however, differs from the other two parallel algorithms in that it does not have a single population. Each processor runs what is effectively the serial GA on a population of reduced size.

Because of the small size of the sub-populations, sub-population stagnation can be a problem, so it is essential that a mechansim exists to exchange individuals between the sub-populations, thereby helping to maintain genetic diversity within each sub-population. The nature or topology of this exchange can be treated in a number of ways. In the island model, for example, exchange is allowed between any pair of sub-populations, while in the ring or stepping-stone model, exchange may only take place between neighbouring populations.

Preliminary work with a parallel version of our GA showed that the sub-population and distributed parallel algorithms afford essentially linear speed-up with the number of processors.⁴⁵ Subsequent testing has shown that, as well as being faster for a fixed total population size and number of generations, the sub-population algorithm tends to find the global minimum more quickly and more reliably than the other algorithms, provided that the sub-populations are not too small.⁵¹

A sequential multi-population approach (sometimes known as a micro-GA) has also been applied to cluster optimisation.^56,60 In this method, a population is evolved over a certain number of generations (epochs), or until the population converges on a single structure. At the end of each epoch, an anihilation (mass extinction) operator removes all or most of the individuals, though high fitness structures from previous epochs may later be seeded back into the population via a memory operator.

In future, it is probable that combining GAs with other search techniques will enable one to combine the best features of stochastic, deterministic and heuristic search methods, involving global as well as local searching. The development of problem-specific, rather than generic search algorithms may also be desireable. In the cluster field, Zacharias et al. have developed a combined SA/GA approach for optimising the structures of silicon clusters (SA = simulated annealing).¹⁰⁶ They report that, for this application, the SA/GA method outperformed individual SA or GA by an order of magnitude (in terms of the CPU time required for convergence).

Hartke has introduced the concept of using an empirical potential to guide an ab initio calculation towards the GM on the ab initio, rather than the empirical hypersurface, without having to perform a large number of costly ab inito calculations.^35,107 His method involves on-the-fly reparameterisation of the empirical potential to a limited number of ab initio calculations, within a GA framework, and has proved successful in the geometry optimisation of small silicon¹⁰⁷ and water clusters.¹⁰⁸ With the advent of faster computers, this technique promises to open up the possibility of global optimisation of moderately sized clusters at high levels of computational sophistication.

Future developments may include the application of self-optimising GAs¹⁵ and Genetic Programming,¹⁰⁹ where the functional form, as well as the parameter values, defining a cluster potential energy function can be determined in an evolutionary fashion, in order to fit experimental or ab initio calculated data. In addition to GAs, recent years have witnessed the development of a large number of computational algorithms based on natural evolution, as well as methods based on Artificial Intelligence (such as Artificial Neural Networks, Expert Systems and Fuzzy or Soft Computing).^18,20 Many of these methods, together with other computational techniques based on nature,¹¹⁰ have already been applied to a number of chemical problems and it is likely that, before long, they will also find application in geometry optimisation or some other aspect of cluster science.

The genetic algorithm is a powerful global optimisation technique, which (together with other, related evolutionary methods) is currently being applied to many problems in the physical and biological sciences. In this Perspective, as well as discussing generic features of GAs and their implementation, I have shown how GAs have been developed to search for low energy isomers for a variety of cluster types. There are a number of search methods which have been applied to global optimisation problems, including simulated annealing and basin hopping.^9,10 It has not been my intention to present a critical comparison of these methods, as it is known that some methods work better than others for specific applications and much depends on the degree to which any given search method has been optimised for a particular problem. Instead, I have aimed to introduce readers to what may be, for them, a new class of optimisation methods which may be of application in their own research.

The importance of GAs (and other reliable search methods) lies not just in their ability to find global minima and low-lying metastable isomers (note that the use of a GA to find transition states and to map out reaction pathways has also been reported¹¹¹) but also because they enable us to test the suitability of a particular potential energy function for simulating a given class of clusters. Of course, the best that a search method can do is to find the global minimum for the potential energy function being tested. If this function is unphysical, then the search method will find a global minimum which is not realistic. As has been alluded to above, however, reliable global optimisation methods (such as the GA) can also facilitate the reparameterisation of potential energy functions by fitting to electronic structure calculations on benchmark structures. The GA can also aid in the use of empirical potentials to guide more sophisticated electronic structure methods towards relevant regions of conformation space.

The current state of the art concerning global cluster geometry optimisation using genetic algorithms is that for most potential energy functions the GM can generally be found reliably, in a reasonable length of time (minutes or hours, rather than days), for homonuclear atomic clusters with of the order of 100–200 atoms. For nanoalloy clusters, because of the occurrence of permutational (homotops) as well as geometric isomers, the size range of confidence is reduced somewhat to around 100. The situation is similar with molecular clusters, where in addition to the molecular positions, their orientation and (for flexible molecules) internal degrees of freedom must be simultaneously optimised.¹¹² It is very likely, however, that the future will see important developments, both in hardware and in genetic algorithm design, which will enable the study of ever-larger clusters, at increasing levels of sophistication.

Acknowledgements

I would like particularly to thank Dr Christopher Roberts, Thomas Mortimer-Jones, Mark Bailey, Sarah Darby and Claire Massen who carried out the work reported here. I also wish to thank Dr Lesley Lloyd, Dr Nicholas Wilson and Freddy Fernandes-Guimaraes who have also been involved in our work on GAs for cluster optimisation and to acknowledge collaborations in this field with Professor Jadson Belchior (Belo Horizonte, Brazil), Dr Alvaro Posada Amarillas (Hermosillo, Mexico), Dr Fred Manby (Bristol) and Dr Said Salhi (Birmingham).

I wish to acknowledge Professor Kenneth Harris and his group for many useful discussions on GAs and for a long and fruitful collaboration in their application to solving structures from powder diffraction data. I thank Professors Peter Edwards, Bernd Hartke, Julius Jellinek and Richard Palmer for helpful discussions on GAs and/or clusters and Dr Hugh Cartwright for introducing me to GAs, through his very readable Oxford Chemistry Primer Applications of Artificial Intelligence in Chemistry. I am indebted to Professors Mike Mingos and John Murrell for introducing me to the topics of cluster chemistry and many-body potential energy functions, respectively.

Finally, I am grateful to the EPSRC, The Leverhulme Trust, The Royal Society and The University of Birmingham for funding this research.

References

1		R. L. Johnston, Atomic and Molecular Clusters, Taylor and Francis, London, 2002.
2		Clusters of Atoms and Molecules, ed. H. Haberland, Springer-Verlag, Berlin, 1994.
3		B. V. Reddy, S. N. Khanna and B. I. Dunlap, Phys. Rev. Lett., 1993, 70, 3323 [Links].
4		M. Moseler, H. Häkkinen, R. N. Barnett and U. Landman, Phys. Rev. Lett., 2001, 86, 2545 [Links].
5		J. Jortner, Z. Phys. D, 1992, 24, 247 [Links].
6		R. L. Johnston, Philos. Trans. R. Soc. London, Ser. A, 1998, 356, 211.
7		Theory of Atomic and Molecular Clusters, ed. J. Jellinek, Springer, Berlin, 1999.
8		S. Erkoç, Phys. Rep., 1997, 278, 80.
9		J. P. K. Doye, in Global Optimization Selected Case Studies, ed. J. D. Pintér, Kluwer, Dordrecht, in press.
10		D. J. Wales, Energy Landscapes, Cambridge University Press, Cambridge, to be published, 2004.
11		B. Hartke, J. Comput. Chem., 1999, 20, 1752 [Links].
12		F. Baletto, C. Mottet and R. Ferrando, Phys. Rev. B, 2001, 63, 155 408.
13		J. P. K. Doye and D. J. Wales, J. Chem. Soc., Faraday Trans., 1997, 93, 4233 [Links].
14		J. Holland, Adaptation in Natural and Artificial Systems, University of Michigan Press, Ann Arbor, MI, 1975.
15		D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading, MA, 1989.
16		M. Mitchell, An Introduction to Genetic Algorithms, MIT Press, Cambridge, MA, 1998.
17		H. M. Cartwright, in Applications of Evolutionary Computing in Chemistry, ed. R. L. Johnston, Struct. Bonding (Berlin) , in press.
18		H. M. Cartwright, Applications of Artificial Intelligence in Chemistry, Oxford University Press, Oxford, 1993.
19		Evolutionary Algorithms in Molecular Design, ed. D. E. Clark, Wiley-VCH, Weinheim, 2000.
20		Soft Computing Approaches in Chemistry: Studies in Fuzziness and Soft Computing, ed. H. M. Cartwright and L. M. Sztandera, Springer-Verlag, Heidelberg, 2003, vol. 120.
21		K. D. M. Harris and R. L. Johnston, Chem. Soc. Rev., in preparation.
22		Applications of Evolutionary Computing in Chemistry, ed. R. L. Johnston, Struct. Bonding (Berlin) , in press.
23		J. T. Pedersen, in Evolutionary Algorithms in Molecular Design, ed. D. E. Clark, Wiley-VCH, Weinheim, 2000, ch. 11, p. 223.
24		R. Unger, in Applications of Evolutionary Computing in Chemistry, ed. R. L. Johnston, Struct. Bonding (Berlin) , in press.
25		V. J. Gillet, in Soft Computing Approaches in Chemistry: Studies in Fuzziness and Soft Computing, ed. H. M. Cartwright and L. M. Sztandera, Springer-Verlag, Heidelberg, 2003, vol. 120, p. 1.
26		K. D. M. Harris, R. L. Johnston and B. M. Kariuki, in Evolutionary Algorithms in Molecular Design, ed. D. E. Clark, Wiley-VCH, Weinheim, 2000, ch. 9, p. 159.
27		K. D. M. Harris, R. L. Johnston, E. Y. Cheung, G. W. Turner, S. Habershon, D. Albesa-Jové, E. Tedesco and B. M. Kariuki, CrystEngComm, 2002, 4, 356 [Links].
28		T. S. Bush, C. R. A. Catlow and P. D. Battle, J. Mater. Chem., 1995, 5, 1269 [Links].
29		S. M. Woodley, P. D. Battle, J. D. Gale and C. R. A. Catlow, Phys. Chem. Chem. Phys., 1999, 1, 2535 [Links].
30		S. M. Woodley, in Applications of Evolutionary Computing in Chemistry, ed. R. L. Johnston, Struct. Bonding (Berlin) , in press, and references therein.
31		B. C. Sanctuary, in Evolutionary Algorithms in Molecular Design, ed. D. E. Clark, Wiley-VCH, Weinheim, 2000, ch. 10, p. 195.
32		D. E. Clark, Evolutionary Algorithms in Computer-aided Molecular Design, URL http://panizzi.shef.ac.uk/cisrg/links/ea_bib.html.
33		B. Hartke, J. Phys. Chem., 1993, 97, 9973 [Links].
34		Y. Xiao and D. E. Williams, Chem. Phys. Lett., 1993, 215, 17 [Links].
35		B. Hartke, Chem. Phys. Lett., 1996, 258, 144 [Links].
36		B. Hartke, Phys. Chem. Chem. Phys., 2003, 5, 275 [Links].
37		B. Hartke, H.-J. Flad and M. Dolg, Phys. Chem. Chem. Phys., 2001, 3, 5121 [Links].
38		Y. Zeiri, Phys. Rev. E, 1995, 51, 2769 [Links].
39		D. M. Deaven and K. M. Ho, Phys. Rev. Lett., 1995, 75, 288 [Links].
40		J. P. K. Doye and D. J. Wales, J. Phys. Chem. A., 1997, 101, 5111 [Links].
41		Z. Li and H. A. Scheraga, J. Mol. Struct. (THEOCHEM), 1988, 179, 333 [Links].
42		D. J. Wales and H. A. Scheraga, Science, 1999, 285, 1368 [Links].
43		A. Tuson and D. E. Clark, in Evolutionary Algorithms in Molecular Design, ed. D. E. Clark, Wiley-VCH, Weinheim, 2000, p. 241.
44		D. M. Deaven, N. Tit, J. R. Morris and K. M. Ho, Chem. Phys. Lett., 1996, 256, 195 [Links].
45		R. L. Johnston and C. Roberts, in Soft Computing Approaches in Chemistry: Studies in Fuzziness and Soft Computing, ed. H. Cartwright and L. Sztandera, Springer-Verlag, Heidelberg, 2003, vol. 120, p. 161.
46		B. Hartke, in Proceedings of the Genetic and Evolutionary Computation Conference, GECCO-2001, ed. L. Spector, E. Goodman, A. Wu, W. B. Langdon, H.-M. Voigt, M. Gen, S. Sen, M. Dorigo, S. Pezeshk, M. Garzon and E. Burke, Morgan Kaufmann, San Francisco, 2001, p. 1284.
47		B. Hartke, in Global Optimization Selected Case Studies, ed. J. D. Pintér, Kluwer, Dordrecht, in press.
48		B. Hartke, in Applications of Evolutionary Computing in Chemistry, ed. R. L. Johnston, Struct. Bonding (Berlin) , in press.
49		R. L. Johnston and C. Roberts, Cluster Geometry Optimization Genetic Algorithm Program, University of Birmingham, 1999.
50		R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, SIAM J. Sci. Comput., 1995, 16, 1190.
51		C. Roberts, Genetic Algorithms for Cluster Optimization, PhD Thesis, University of Birmingham, 2001.
52		C. Roberts, R. L. Johnston and N. T. Wilson, Theor. Chem. Acc., 2000, 104, 123 [Links].
53		F. R. Manby, R. L. Johnston and C. Roberts, MATCH, 1998, 38, 111.
54		R. L. Johnston, T. V. Mortimer-Jones, C. Roberts, S. Darby and F. R. Manby, in Applications of Evolutionary Computation, ed. S. Cagnoni, J. Gottlieb, E. Hart, M. Middendorf and G. R. Raidl, Lect. Notes Comput. Sci., 2002, 2279, 92.
55		C. Roberts and R. L. Johnston, Phys. Chem. Chem. Phys., 2001, 3, 5024 [Links].
56		F. F. Guimarães, J. C. Belchior, R. L. Johnston and C. Roberts, J. Chem. Phys., 2002, 116, 8327 [Links].
57		L. D. Lloyd, R. L. Johnston, C. Roberts and T. V. Mortimer-Jones, Chem. Phys. Chem., 2002, 3, 408 [Links].
58		A. Posada Amarillas, M. F. Ortíz, C. Roberts, T. V. Mortimer-Jones and R. L. Johnston, in preparation.
59		S. Darby, T. V. Mortimer-Jones, R. L. Johnston and C. Roberts, J. Chem. Phys., 2002, 116, 1536 [Links].
60		R. A. Lordeiro, F. F. Guimarães, J. C. Belchior and R. L. Johnston, Int. J. Quantum Chem., 3003, 25, 41.
61		M. S. Bailey, N. T. Wilson, C. Roberts and R. L. Johnston, Eur. Phys. J. D, in press.
62		C. Massen, T. V. Mortimer-Jones and R. L. Johnston, J. Chem. Soc., Dalton Trans., 2002, 4375 [Links].
63		L. D. Lloyd, R. L. Johnston and S. Salhi, in preparation.
64		R. L. Johnston and N. T. Wilson, Birmingham Cluster Web , URL: http://www.tc.bham.ac.uk/bcweb.
65		D. J. Wales, J. P. K. Doye, A. Dullweber, M. P. Hodges, F. Y. Naumkin and F. Calvo, The Cambridge Cluster Database , URL: http://www-wales.ch.cam.ac.uk/CCD.html.
66		K. D. M. Harris, R. L. Johnston and B. M. Kariuki, Acta Crystallogr., Sect. A, 1998, 54, 632 [Links].
67		S. Habershon, K. D. M. Harris, R. L. Johnston, G. W. Turner and J. M. Johnston, Chem. Phys. Lett., 2002, 353, 185 [Links].
68		T. M. Köhler, H.-P. Gail and E. Sedlmayr, Astron. Astrophys., 1997, 320, 553.
69		G. V. Lewis and C. R. A. Catlow, J. Phys. C: Solid State Phys., 1985, 18, 1149 [Links].
70		P. J. Ziemann and A. W. Castleman Jr., J. Chem. Phys., 1991, 94, 718 [Links].
71		R. L. Whetten, Acc. Chem. Res., 1993, 26, 49 [Links].
72		E. de la Puente, A. Aguado, A. Ayuela and J. M. López, Phys. Rev. B, 1997, 56, 7607 [Links].
73		M.-J. Malliavin and C. Coudray, J. Chem. Phys., 1997, 106, 2323 [Links].
74		J. M. Recio, R. Pander, A. Ayuela and A. B. Kunz, J. Chem. Phys., 1993, 98, 4783 [Links].
75		M. Wilson, J. Phys. Chem. B, 1997, 101, 4917 [Links].
76		Ch. Chang, A. B. C. Patzer, E. Sedlmayr, T. Steinke and D. Sülze, Chem. Phys. Lett., 2001, 350, 399 [Links].
77		W. D. Knight, K. Clemenger, W. A. de Heer, W. A. Saunders, M. Y. Chou and M. L. Cohen, Phys. Rev. Lett., 1984, 52, 2141 [Links].
78		T. P. Martin, Phys. Rep., 1996, 273, 199 [Links].
79		J. Jellinek and E. B. Krissinel, in Theory of Atomic and Molecular Clusters, ed. J. Jellinek, Springer, Berlin, 1999, p. 277, and references therein.
80		S. Giorgio, H. Graoui, C. Chapan and C. R. Henry, in Metal Clusters in Chemistry, ed. P. Braunstein, L. A. Oro and P. R. Raithby, Wiley-VCH, Weinheim, 1999, vol. 2, p. 1194.
81		B. Pauwels, G. Van Tendeloo, E. Zhurkin, M. Hou, G. Verschoren, L. Theil Kuhn, W. Bouwen and P. Lievens, Phys. Rev. B, 2001, 63, 165 406.
82		A. V. Ruban, H. L. Skriver and J. K. Norskov, Phys. Rev. B, 1999, 59, 15 990.
83		G. Bozzolo, J. Ferrante, R. D. Noebe, B. Good, F. S. Honecy and P. Abel, Comput. Mater. Sci., 1999, 15, 169 [Links].
84		A. M. Molenbroek, S. Haukka and B. S. Clausen, J. Phys. Chem. B, 1998, 102, 10680 [Links].
85		G. Schmid, in Metal Clusters in Chemistry, ed. P. Braunstein, L. A. Oro and P. R. Raithby, Wiley-VCH, Weinheim, 1999, vol. 3, p. 1325.
86		M. P. Andrews and S. C. O'Brien, J. Phys. Chem., 1992, 96, 8233 [Links].
87		M. J. López, P. A. Marcos and J. A. Alonso, J. Chem. Phys., 1996, 104, 1056 [Links].
88		A. Fortunelli and A. M. Velasco, J. Mol. Struct. (THEOCHEM), 1999, 487, 251 [Links].
89		M. Calleja, C. Rey, M. M. G. Alemany, L. J. Gallego, P. Ordejón, D. Sánchez-Portal, E. Artacho and J. M. Soler, Phys. Rev. B, 1999, 60, 2020.
90		S. Bromley, G. Sankar, C. R. A. Catlow, T. Maschmeyer, B. F. G. Johnson and J. M. Thomas, Chem. Phys. Lett., 2001, 340, 524 [Links].
91		P. Nava, M. Sierka and R. Ahlrichs, Phys. Chem. Chem. Phys., 2003, 5, 3372 [Links].
92		J. Wang, G. Wang, X. Chen, W. Lu and J. Zhao, Phys. Rev. B, 2002, 66, 14419.
93		F. Cleri and V. Rosato, Phys. Rev. B, 1993, 48, 22 [Links].
94		I. L. Garon, K. Michaelian, M. R. Beltraan, A. Posada-Amarillas, P. Ordejón, E. Artacho, D. Sánchez-Portal and J. M. Soler, Phys. Rev. Lett., 1998, 81, 1600 [Links].
95		N. T. Wilson and R. L. Johnston, J. Mater. Chem., 2002, 12, 2913 [Links].
96		C. Rey, J. García-Rodeja and L. J. Gallego, Phys. Rev. B, 1996, 54, 2942 [Links].
97		E. F. Rexer, J. Jellinek, E. B. Krissinel, E. K. Parks and S. J. Riley, J. Chem. Phys., 2002, 117, 82 [Links].
98		F. R. de Boer, R. Boom, W. C. M. Mattens, A. R. Miedama and A. K. Niessen AK, Cohesion in Metals: Transition Metal Alloys, Elsevier, Amsterdam, 1988.
99		A. J. Renouprez, J. L. Rousset, A. M. Cadrot, Y. Soldo and L. Stievano, J. Alloys Compd., 2001, 328, 50 [Links].
100		J. L. Rousset, L. Stievano, F. J. Cadete Santos Aires, C. Geantet, A. J. Renouprez and M. Pellarin, J. Catal., 2001, 202, 163 [Links].
101		J. M. Montejano-Carrizalez, M. P. niguez and J. A. Alonso, Phys. Rev. B, 1994, 49, 16649 [Links].
102		A. M. Schoeb, T. J. Raeker, L. Q. Yang, X. Wu, T. S. King and A. E. DePristo, Surf. Sci., 1992, 278, L125 [Links].
103		F. Baletto, C. Mottet and R. Ferrando, Phys. Rev. Lett., 2003, 90, 135 504 [Links].
104		G. von Helden, M.-T. Hsu, N. Gotts and M. T. Bowers, J. Phys. Chem., 1993, 97, 8182 [Links].
105		A. R. Leach, Molecular Modelling: Principles and Applications, Adison-Wesley Longman, Harlow, 1996.
106		C. R. Zaharias, M. R. Lemes and A. Dal Pino, J. Mol. Struct. (THEOCHEM), 1998, 430, 29 [Links].
107		B. Hartke, Theor. Chem. Acc., 1998, 99, 241 [Links].
108		B. Hartke, M. Schütz and H.-J. Werner, Chem. Phys., 1998, 239, 561 [Links].
109		J. R. Koza, Genetic Programming: On the Programming of Computers by Means of Natural Selection, MIT Press, Cambridge, MA, 1992.
110		Computing With Biological Metaphors, ed. R. Paton, Chapman and Hall, London, 1994.
111		P. Chaudhury, S. P. Bhattacharyya and W. Quapp, Chem. Phys., 2002, 253, 295 [Links].
112		B. Hartke, Z. Phys. Chem., 2000, 214, 1251 [Links].