Structure prediction of titania phases: Implementation of Darwinian versus Lamarckian concepts in an Evolutionary Algorithm

Computational Materials Science 45 (2009) 84–95

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Structure prediction of titania phases: Implementation of Darwinian versus Lamarckian concepts in an Evolutionary Algorithm S.M. Woodley *, C.R.A. Catlow Davy Faraday Research Laboratory, Department of Chemistry, 3rd Floor, Kathleen Lonsdale Building, Gower Street, London WC1E 6BT, United Kingdom

a r t i c l e

i n f o

Article history: Received 1 November 2007 Received in revised form 6 February 2008 Accepted 14 February 2008 Available online 10 September 2008 PACS: 2.70.c 02.70.Uu 61.50.Ah 61.66.f Keywords: Global optimisation Genetic algorithms Semi-classical simulations Interatomic potentials Titania polymorphs

a b s t r a c t Darwinian and Lamarckian schemes within evolutionary algorithms have been implemented and optimised. We compare the performance of these two approaches applied to the problem of structure prediction of the titania polymorphs. A number of well-known phases have thus been reproduced as well as several plausible novel microporous and dense structures. Two different potential parameter sets, within the Born model description of a solid, were employed. Following the Lamarckian concept in a genetic or more generally in an evolutionary algorithm, all new candidate structures are immediately relaxed (analogous to the ageing process in nature); consequently, competition within any current population to procreate only occurs between mature candidate structures, which correspond to local minima on the energy landscape. In the Darwinian scheme, no local optimisation is performed, which should result in significant saving in CPU time per candidate structure considered. We show, however, that the Lamarckian scheme (which ultimately searches for the global minimum on a simplified landscape) is more successful and efficient at generating the target structures. Analysis of why the Lamarckian scheme produced a perfect success rate uncovered a weakness in the Darwinian approach when diversity of the population is allowed to reduce, and further methodological developments are suggested. Ó 2008 Elsevier B.V. All rights reserved.

1. Introduction Evolutionary algorithms have been developed and applied in many different fields, as demonstrated in the Genetic Algorithm Symposium in the 2007 EMRS autumn meeting, with applications ranging from optimising the process of bending steel plates for hulls of ships [1] and reduction of iron ore coal mixture in a packed bed reactor [2], to finding the most energetically favoured arrangement of water molecules [3] and the atomic structure of steps on stable crystal surfaces [4]. In the application to predicting organic crystal structures [5], drug companies are very interested in knowing what polymorphs an organic molecular solid can adopt as the different polymorphs may have different properties and, for example, a drug with the wrong polymorph that is administered as a pain killer could potentially be a source of a dangerous overdose. In another application, that of small sized nanoparticles [6–8], current experimental techniques cannot yield the information needed for a detailed understanding of (a) the structure, (b) their growth mechanisms, or (c) their specific properties. Computer simulations provide a complementary

* Corresponding author. Tel.: +44 (0) 20 76790315. E-mail address: [email protected] (S.M. Woodley). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.02.033

tool, in which key, low energy, configurations are found, and the relation between size of particle and their properties established. In particular, different structure prediction techniques have been employed to shed light on why certain sized clusters are more readily formed during their creation from either laser ablation or nucleation in water. In this work we focus on the problem of crystal structure solution and prediction. In the case of bulk inorganic [9–11], organic [12,13], and inorganic–organic hybrid [14,15] materials that can only be synthesised in a powder form, only the unit cell dimensions, and possibly symmetry operations that the structure must obey, can be readily extracted and, importantly, an approximate structure must be deduced before the structure can be solved using powder refinement techniques [16–19]. We are interested in developing an efficient and reliable approach for predicting the low energy polymorphs of inorganic, crystalline materials. We assume that the dimensions of the unit cell and the ionic contents are known, from diffraction data, and by chemical analysis and knowledge of the preparation procedure for example. We have therefore developed a computational approach to predict the possible structures within a fixed, predefined cell unit and constituent ions [20]. In this article, we report two evolutionary approaches and compare their success

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

in predicting several target structures. In line with previous authors [20–22], we choose, as our test system, to model the possible structures of titania. However, we do not restrict our study to the common phases, rutile [23], anatase [23] and brookite [24], but also include the metastable microporous structures, Hollandite [25] and Ramsdellite [26]. In Fig. 1, we show the structure of each of the target phases. During our search (when the cell parameters were also relaxed), other high-pressure phases of TiO2 were found, including columbite, baddeleyite, fluorite and cotunnite, although some of these phases have not been synthesised in the laboratory. We have developed our implementation of an evolutionary algorithm as a module within the General Utilities Lattice Program (GULP) [27]. Our method is based upon a multi-stage approach, whereby a genetic or evolutionary algorithm [28] is employed to find feasible approximate structural solutions that are then subsequently relaxed, using the main modules already available in GULP, so that the lattice energy, as defined using the Born model of a solid, is minimised. The code has been developed so that a range or variants of global optimisers are available; examples include, a genetic algorithm (GA) that uses genotype operators, a more general evolutionary algorithm (EA) that uses phenotype operators, and other Monte Carlo approaches, selectable via the input file. Using a Darwinian approach within the GA, the multistage scheme proved successful in the prediction of a previously unknown structure for lithium ruthenate, Li3RuO4 [29]. Simulated annealing, Monte Carlo approaches (SA-MC), have also been developed for predicting inorganic crystalline materials [30–32]. Likewise, in the field of predicting the most stable structures for small-sized inorganic nanoparticles, both GA and SAMC, or SA-Molecular Dynamics, have been employed, but more importantly it is reported that hybrid methods, whereby each new candidate structure is immediately relaxed, have proved to be more reliable in locating the global minimum structures [33]; this approach is adopted by the so-called hybrid-GA and

85

Monte Carlo basin hopping (MCBH) methods. The latter technique has been popularised by Wales, who has predicted the global minimum structures for Lennard–Jones clusters containing 1– 110 particles [34]. In both cases, one can say (although perhaps less convincingly for MCBH) that Lamarckian concepts have been used. For a GA, Darwinian evolution is simulated whereby the genetic makeup of a candidate structure, the unknown ionic coordinates, does not change during its lifetime so, ignoring mutations created during procreation, offspring will resemble their parents, grandparents and so forth. For a hybrid-GA, Lamarckian evolution is simulated; the genetic makeup can change during the life of a candidate structure, when directly minimising the candidate’s energy, and the modified genetic makeup is passed on to the next generation so offspring may not resemble their grandparents. The latter, therefore, combines global and local optimisation techniques during the global search. Below we present more details of our implementation to predicting the global minimum structures [35–37], albeit for bulk materials. The success of hybrid approaches is not confined to clusters; Turner et al. [38] showed that a hybrid-GA gave significant improvements in efficiency and reliability when used for organic crystal structure solution. In another recent application of crystal structure prediction to the problem of carbon polymorphism, Abraham and Probert [39] adapted a variable number approach of Chuang et al. [40] within their hybrid-GA so that, during the search, the number of carbon atoms within the candidate’s supercell need not be constrained to a predefined value. Further hybrid methods include using MC, rather than a local optimiser during the mature phase of the GA [41] and the use of MCBH as steps within a set of MD runs, with temperature decreasing [36], both applied successfully to predicting low energy structures of carbon and titania clusters, respectively. Here we report and compare the performance, in predicting inorganic crystal structures, of an EA incorporating either Darwinian or Lamarckian concepts.

Fig. 1. The titania target phases, (a) rutile, (b) anatase (IP2), (c) distorted anatase (IP1), (d) brookite, (e) Ramsdellite (IP2), (f) distorted Ramsdellite (IP1), (g) Hollandite (IP2) and distorted Hollandite (IP1).

86

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

sible values for structures containing short interatomic distances, we added an extra repulsive term so that:

2. Methodology 2.1. Comparing the quality of candidate structures: interatomic potentials

Vðrij Þ

Global optimisation techniques require a quick and robust way to compare candidate solutions. As we are interested in finding the global minimum structure, and possibly other low-lying, meta-stable structures, the lattice energy is an ideal ‘‘measuring stick” when comparing how good each particular candidate structure is within a set, or ‘‘population”, of candidate structures. Moreover, as discussed below, the set of lattice energy values within the current population of solutions governs, or weights, the ‘‘selection probability” so that better candidate structures are more likely to procreate. We note, however, as used by Chakraborti when applying multiobjective optimisations [42], if we are trying to predict a stable energy phase with certain desired properties, then the Pareto front should be optimised rather than just the energy or a fixed weighted sum of the components. As many candidate structures need to be assessed, we choose to use a two-body, rigid ion model to describe the solid [43]. The lattice energy per unit cell is simply:

E¼

X

kij Eij ;

ð1Þ

ij

where ij labels all the pair-wise interactions and kij is 1, ½ or 0 when both, either or neither ions i and j are within the defined unit cell, respectively. Typically, Eij is reported as the sum of two terms,

Eij ¼ k

qi qj r ij

Vðr ij Þ;

þ

ð2Þ

a long-range interaction, or the Coulomb contribution between the point charges, qi and qj, that are a distance rij apart, and the shortrange interactions, V(rij), between ions i and j. The standard Ewald method [44,45] is employed to compute the sum of Coulomb contributions, whereas a cut-off of 15 Å was used for contributions to V(rij). We have chosen to use two sets of interatomic potential (IP) parameters [46,47], which were empirically fitted to reproduce the observed bulk structure and a number of physical properties, including the elastic constants, of the rutile phase for TiO2. One set, used previously in a study by Freeman et al. [46] to predict rutile, anatase and brookite phases of TiO2, is based on formal charges and another, by Matsui and Akaogi [47], which was subsequently employed to model the surfaces of rutile and anatase phases [48], used partial charges. We will refer to the former set as IP1 and the latter set as IP2. Buckingham short-range terms,

Aij expðrij =qij Þ 

C ij ; r6ij

ð3Þ

are included in both potential sets, where A, q, C and q, have values given in Table 1. In order to make these expressions for the lattice energy more robust, so that our procedure should also provide sen-

Table 1 Potential parameters, IP1 [46] and IP2 [47], for modelling titania; see Eq. (4)

IP1

IP2

i

j

qi |e|

Ti O Ti

Ti O O

4.000 2.000

Ti O Ti

Ti O O

2.196 1.098

Aij (eV)

qij (Å)

Cij (Å6)

22,764.3 754.2

0.1490 0.3879

27.88

31,120.1 11,782.7 16,957.5

0.1540 0.2340 0.1940

5.25 30.22 12.59

¼

Aij expðr ij =qij Þ þ

Bij C ij  6: r12 r ij ij

ð4Þ

A value of 1.0 Å12, for all interactions, was chosen for Bij so that when a candidate structure is relaxed, using standard optimisation procedures to minimise the lattice energy, ions that are too densely packed would separate to produce structures that have plausible interatomic distances. Moreover, the inclusion of this additional term should cause minimal changes for plausible structures and, beneficially, removes the possibility of simulating unphysical structural collapse. To obtain an insight into the quality of the chosen interatomic potentials, we computed two lattice energies, Ev and Ep, for the most commonly known dense phases of titania, namely rutile, anatase and brookite, as well as two framework phases, namely Ramsdellite and Hollandite – see Table 2 and Fig. 1. Ev is the local minimum on the lattice energy landscape obtained by relaxing the ionic coordinates, whilst constraining the lattice parameters, and where all structural parameters are initially set to those observed [23–26,49], whereas Ep is the calculated enthalpy at ambient pressure, i.e. lattice parameters are also relaxed. For the dense phases, an acceptable change in the unit cell parameters (less than a change of 0.02 in the values shown in the last column of Table 2) and a negligible difference in our calculated values for Ev and Ep (less than 0.005%) was found when V(rij) is defined by Eq. (3) rather than (4). Comparing values for enthalpies, the order of stability of the three phases is correctly predicted when IP1 is employed. However, the average change, of about 5%, in the lattice parameters is very poor (although fortunately we will be conducting constant volume simulations when testing our algorithms). The lattice energy values obtained when employing IP2 are somewhat artificial, as the Coulomb contribution does not use the formal oxidation states of the ions, and erroneously locates brookite as the most stable of the three phases. However, the reproduction of the crystal structures for all three dense phases is far superior to that predicted when employing IP1, with an average change in the lattice parameters that is less than 1%. After relaxation, but with the dimensions of the unit cell fixed, we only obtained a good reproduction of the microporous framework structures of Ramsdellite and Hollandite when using IP2. The relaxation when using IP1 is pronounced; compare Fig. 1e with f, and g with h. In the Hollandite structure we obtain parallel 1 D rods, or tubes, where the cations are at the corners of the square cross section of these rods. Each TiO6 octahedron has one elongated Ti–O bond of 2.89 Å (not shown in Fig. 1h), which connects neighbouring rods. When the cell parameters are also relaxed, the cell volume for both phases expands by over 40%, where we obtain the microporous frameworks shown in Fig. 2. The change in the

Table 2 Calculated lattice energies for TiO2 and the difference in calculated and observed [23,24] values for the unit cell volume at zero pressure Phase

Ev eV/TiO2

Ep eV/TiO2

DVolume (%)

IP1

Rutile Anatase Brookite Ramsdellite Hollandite

112.133 112.135 112.333 111.330 111.766

112.50 112.761 112.422

4.25 6.48 5.27 >40 >40

IP2

Rutile Anatase Brookite Ramsdellite Hollandite

39.759 39.495 39.617 39.208 39.161

39.800 39.496 39.622 39.395 39.190

2.72 0.17 1.26 16.89 1.52

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

87

Fig. 2. The local minimum energy (IP1), microporous, titania structure obtained by relaxing the (a) Hollandite framework and (b) Ramsdellite framework, under constant pressure.

cell volume is much better when using IP2: 1.5% for Hollandite and 17% for Ramsdellite – a reasonable agreement given that the potential parameters were fitted to model dense phases of titania and that interstitial atoms/ions are normally observed in the framework structures. Even though the use of IP1 is poor for modelling the framework structures, we still can use the structures shown in Fig. 1 as our targets. We note that Density Functional Theory (DFT) approaches to calculating the energy for each new candidate structure during a global search have also been used, for example, to predict the surface structure of oxides [50]. However, as the computational task of DFT methods is much more expensive, their application is more frequently used in a later stage, when refining only the plausible or hypothetical structures that are predicted by the global search algorithm. In some cases, the assignment of the structure which corresponds to the global minimum may change. Thus, the additional stage simultaneously provides a check on the reliability of the IP parameters initially employed to model the chosen system, which is useful when observed data are not available, for example, when predicting the structure of clusters using IP parameters, particularly if the IP parameters were originally fitted to reproduce the observed structure of the bulk phase [51]. 2.2. Global optimisation: Darwinian and Lamarckian concepts in Evolutionary Algorithms Genetic or evolutionary algorithms evolve a population of solutions (candidate structures), where hopefully one candidate will best fit the chosen fixed environment (i.e. definition of lattice energy, unit cell dimensions and contents), by simulating competition to procreate and survive. There are many good textbooks on genetic or evolutionary algorithms, for example see reference [28], and general details of our implementation are also published elsewhere [9,20,37], so here we present only key aspects and a few important details. In a previous study, concerned with the structure prediction of microporous silica frameworks, we found our algorithm worked more efficiently when phenotype, rather than genotype, operators were used during the procreation of new candidate structures. Thus, we have selected to use only phenotype operators, or the more general EA. The algorithms used in this study, shown in Fig. 3, have evolved themselves over the last decade; in particular, following Johnston [38] we now use natural selection with elimination of duplicate candidates. The main differences between the two algorithms are (i) smaller population size and (ii) new candidate structures are immediately relaxed when the Lamarckian scheme, as opposed to Darwinian scheme, is used. Throughout this article we will refer to each method as D-EA and L-EA, respectively. The elimination of duplicates helps to keep the diversity of the population, which is particularly important for the L-EA where we use a smaller population size. To simulate competition to procreate, a two-member

tournament is used rather than the roulette wheel. As pointed out in our workshop [52], as well as being more expensive to use when compared to an N-member tournament, the roulette wheel should be used carefully as it can easily kill diversity. Inherited from our original GA, a random structure is created so that each ion is constrained to a different site on a 64  64  64 grid. 3. Results and discussion Ten cycles of the L-EA were first carried out in a search for the lowest energy structures with the cell dimensions and constituents of rutile (with the number of formula units, Z = 2), anatase (Z = 4), Ramsdellite (Z = 4), brookite (Z = 8) and Hollandite (Z = 8). Fifty runs were made for each system, whereby the initial population contained a different set of 20 random candidate structures. For each run, the best 10 candidate structures in the twentieth and final population were stored for further refinement in a second stage. During the ageing process of each new candidate in stage one, whereby its structure is relaxed using the method of conjugate gradients, tolerances for the accuracy of the coordinates, energy and forces of 105, 107, 105, respectively, were set. Although there are more sophisticated and reliable approaches, the criterion used to distinguish possible duplicate structures was simply that the difference in lattice energy should be less than 0.01 eV/unit cell. Five different structural refinement procedures, which the lattice energy is minimised, were tried as our second stage; in the first four, constant lattice parameters were maintained throughout the optimisations, whilst in the final method, they were allowed to relax. For convenience, we refer to these approaches as BFGS, RFO, SWITCH, RFO2 and CONP. In the first approach, the BFGS scheme was used to update the Hessian (which is used to determine the direction of the linear line searches). In the RFO approach, a rational function optimiser [53] is used to obtain the local minimum of the lattice energy with respect to atomic coordinates. When the initial candidate structure is not within a quadratic region (typically candidates produced from the D-EA after a few cycles), we sometimes found that the RFO approach became stuck on a non-stationary point of the energy landscape. With the RFO2 scheme, by restarting the RFO procedure from the last point on the energy landscape (and tightening the tolerances), we found that the optimiser was able to find a local minimum. Alternatively, we also tried using a SWITCH, whereby the BFGS scheme was replaced with that of RFO once the gradient norm fell below 0.1. With all, apart from CONP (where we choose to use the RFO optimiser), the cell parameters were constrained to those of the target structure. Normally, CONP would be used as a stage 3. Rutile (Z = 2): For the IP1 landscape, the average time (as measured when using a modern PC with a single processor) for each run was 256 s and, on average, the subroutine for evaluating the energy was called 99,853 times, and 22,472 times when the deriv-

88

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

first cycle

Generate 2M random candidate structures (create population)

Expand population: create L random candidate structures

Procreation: Select 2M candidates such that the better candidates are more likely to be chosen

Elitism: all candidates selected from population

For each pair of successful candidates create two new candidate structures using crossover and mutation routines

next cycle

Evaluate cost function for candidates in population to determine likelihood of each candidate to procreate (become a parent)

Natural Selection: From 2M+N candidate structures select the top 2M non-duplicate candidates as the new population On last cycle - store the best structures in current population for stage 2

first cycle

Generate 2M random candidate structures (create population)

Expand population: create L random candidate structures

Lamarckian or Ageing Process: Relax all candidate structures: growing up

Procreation: Select 2M candidates such that the better candidates are more likely to be chosen and create two new candidate structures using crossover & mutation routines + Relax Children Structures

Elitism: all candidates selected from population

next cycle

Evaluate cost function for candidates in population to determine likelihood of each candidate to procreate (become a parent)

Natural Selection: From 2M+N candidate structures select the top 2M non-duplicate candidates as the new population On last cycle - store the best structures in current population for stage 2

Fig. 3. Schematic diagrams of our evolutionary algorithms where (a) a Darwinian approach and (b) a Lamarckian approach are used in the search for the target structures of titania.

atives were also required. Similar numbers, 301 s, 115,410 and 25,702 times, were found for the IP2 landscape. For IP1, the 10 best candidates in each population, for all fifty runs, were found to be within the basin, on the energy landscape, that contains the rutile structure. In fact, the rutile structure was always obtained after just one L-EA cycle. For IP2, slightly different results were obtained depending on whether the BFGS or RFO approach was used for stage 2. Considering all the best 10 candidates from each run, we found 403 or 405 candidates (including all the best candidates found on each run) adopted the target structure of rutile, while 97 or 95 candidates adopted a metastable, layered structure shown in Fig. 4c. For a few candidate structures, BFGS and RFO obtained different structures because these candidates were not near a true local minimum on the IP2 landscape.

In an attempt to make a fair comparison, we set up D-EA runs so that each took, on average, approximately the same time to run as one L-EA cycle. In more detail, for each run we evolved a population of 100 candidates over 140 GA cycles, which, on average took 30 s (evaluating the cost function 7519 times). For IP1, 479 and 497 of the 500 candidates relaxed to the rutile structure using the BFGS and RFO approach, respectively. The remaining structures also relaxed to the target phase after further refinement and after we had removed an imaginary frequency mode and/or a better tolerance was used. For IP2, 489 (495 when imaginary frequency mode removed or better tolerance used) and 499 candidates adopted the rutile phase and the remaining, including a candidate from one run that was labelled the best candidate after stage 1, adopted the layered structure shown in Fig. 4. If success was measured by how

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

89

Fig. 5. The fluorite structure as (a) generated when targeting rutile and (b) viewed within the full cell.

Fig. 4. Two different layered IP2 structures where (a) and (b) are stable under ambient pressure and (c) and (d) are stable when the unit cell dimensions are constrained to that of rutile.

Table 3 The distribution of candidate structures generated using an evolutionary algorithm, followed by either BFGS or RFO (see text), where rutile (Z = 2) is the target phase Polymorph

L-EA (BFGS)

L-EA (RFO)

D-EA (BFGS)

D-EA (RFO)

Energy (eV/uc)

Rutile Rutile Fig. 4c and d Fig. 4a and b

500 402 (403) 96 1

500 405 95 0

479 (500) 489 (495) 5 0

500 497 (499) 1 0

224.5650 79.5176 75.0864 75.0195

often the best candidate generated on each run was the target, then the D-EA scored a 100% (IP 1) and 98% (IP 2) success rate, compared to the 100% (IP 1 and 2) success rate found for the LEA. The distributions of the different structures found in the final populations are tabulated in Table 3, which draws a picture of the IP1 and IP2 landscapes being dominated by the basin containing the GM. Unsolved structures, typically having many more degrees of freedom, would prove more challenging than the rutile phase, which does not need the application of either a D-EA or LEA as it can trivially be generated by optimising several random structures. We note that the addition of the repulsive short-range r12 term between all ions modifies the landscapes so as to increase the success rate of finding the target by simply relaxing random structures. Previously, Freeman et al. [46], generated the rutile, anatase and brookite phases on 43, 4 and 3 out of 50 attempts, respectively, (86%, 8% and 6%) when using a Metropolis MC approach on the IP1 landscape, although here the repulsive short-range r12 term was only used in a pre-stage to ensure that ions within each initial random structure are adequately spaced. Hence, although global optimisation techniques are not necessary for rutile, we may find anatase and brookite more challenging for our EAs. More interestingly, the results from the GA, using IP1 and CONP, generated more than just the rutile phase (although, of course, if CONP were used as a stage 3, then all candidates would have the rutile structure). In fact, for IP1, 248 candidates adopt the rutile structure (226.376 eV/uc), 74 candidates – anatase (225.521 eV/triclinic unit cell with a = c = 3.614 Å, b = 6.116 Å, a = 72.82° and c = 107.18°) and five – fluorite (223.828 eV/uc; full  cell has a = 4.915 Å with space group Fm3m, as shown in Fig. 5). Clearly, the candidate structures obtained from the D-EA runs are not sufficiently close to the minimum for the rutile phase on the enthalpy landscape in that (as defined by IP1) the more stable

phase of anatase and the less stable phase of the fluorite structure are generated. Considering only the best candidate from each run, 26 out of the 50 adopted the target structure, or rutile, and seven – anatase. 17 (or 169/500) candidates adopted a structure not yet discussed, which is composed of five coordinated titanium ions, with an energy of 226.376 eV/uc (space group Imma and lattice parameters of the full cell, shown in Fig. 6, of 3.625, 5.858, 9.413 Å) and, thus, the global minimum (GM) on the enthalpy landscape (at ambient pressure). As shown in Fig. 6d and c, the structure has hexagonal channels and anatase-style ladders. For IP2, when using CONP as the second stage after either the LEA or D-EA, although the target rutile structure has been found, a similar zoo of structures was generated as seen for IP1 when DEA was followed by CONP. Specifically, for the L-EA, 400, 1, 26 and 12 candidates adopted the rutile, anatase, fluorite and the layered structures shown in Fig. 4a, respectively. Two further structures were also found; with lattice energy of 78.376 eV, 60 adopted the structure shown in Fig. 7 (full cell has parameters a = 3.376, b = 4.376, c = 3.376 Å, space group Pmn21) and, with a higher lattice energy of 77.376 eV, one adopted the structure shown in Fig. 8 (full cell has parameters a = 2.376, b = 6.376, c = 5.376 Å, space group Cmcm). For the D-EA, 362, 2, 119 and 3 candidates adopted respectively the structure of rutile, anatase, that shown in Fig. 7 and fluorite. Seven other stable structures, shown in Fig. 9, with lattice energies of 77.376, 77.488 (Ramsdellite), 77.404, 77.289, 76.946 (titania (II) phase, see below), 76.846, 76.791 eV/uc, were also generated. Anatase (Z = 4): The distribution of structures found in the final populations for the anatase system is shown in Table 4. With twice the number of atoms within the chosen unit cell, compared with

Fig. 6. The IP1 GM structure found when (a) initially targeting rutile and (b–d) viewed within the full cell along three different directions.

90

S.M. Woodley, C.R.A. Catlow / Computational Materials Science 45 (2009) 84–95 Table 4 The distribution of candidate structures after stage two, where an evolutionary algorithm is used in stage one and where anatase (Z = 4) is the target phase

Fig. 7. A higher energy (IP2) structure (a) found when targeting rutile and (b–d) as viewed within the full cell from three different directions.

Fig. 8. A higher energy (IP2) structure (d) found when targeting rutile and (a–c) as viewed within the full cell from three different directions.

Polymorph

L-EA (BFGS)

Anatase buckled Cu2O (Fig. 10b) Split phase (Fig. 10c) Fig. 10d Fig. 10e Fig. 10f Other

227 + 5 (232) 232 (232) 96 (93)

Anatase Fig. 11a Layered Fig. 11b Fig. 11c Figs. 11d and e Fig. 11f Other

262 (262) 101 (101)

268 (268)

L-EA (RFO)

D-EA (BFGS)

D-EA (RFO)

148 (148) 448.7552

268 (268) 359 (373) 310 (310) 448.4814 4 (0) 22 (17) 12 (11) 5 (1) 2 (5)

50 (50) 52 (52) 1 (27) 27 (3) 8 (5)

Energy (eV/uc)

262 (262) 143 (145) 101 (101) 46 (50) 4 (0) 50 (50) 29 (26) 52 (52) 57 (55) 0 (27) 0 (50) 27 (3) 50 (57) 8 (5) 171 (117)

447.4160 13 (13) 26 (26) 3 (3)

447.0540 447.0427 446.7681 < 444.08410

159 (165) 157.9790 87 (28) 155.6858 155.0711 28 (30) 154.8568 57 (62) 154.7481 0 (28) 154.6813 27 (30) 154.6784 142 (95) < 154.6166

stable (448.755 eV/uc) than the unbuckled structure, and is therefore considered as the target structure. Note that, once the cell parameters are no longer constrained, the (unbuckled) anatase phase is lower in energy and can immediately be obtained during most structural relaxation procedures. With BFGS, all 50 runs had at least one candidate that adopted the buckled anatase structure in the final population. With a slightly higher energy, 448.713 eV/uc, 5 out of the final 500 candidates were found to adopt an unstable form of the anatase structure that was buckled only in one direction, as shown in Fig. 10a. As found in previous studies by Freeman et al. [46] and Woodley et al. [20], the Cu2O structure, shown in Fig. 10b, was also readily generated (268 out of the 500 candidates). Using RFO, the 5 partially buckled structures relaxed to the target buckled anatase structure and again, even if we had only considered the candidates within the first population of each run, a 100% success rate was found. Using CONP as a third stage, all buckled anatase structures straighten out and the

Fig. 9. Metastable (IP2) structures found when targeting rutile, with the space group (a) Cmc21, (b) C2/m, (c) P4/nmm, (d) P4/mmm, (e) C2221, (f) Amm2, (g) Pmc21.

rutile, one can expect more possible metastable phases and, perhaps, either a smaller success rate or a longer search time during stage 1. For IP1, with the unit cell parameters fixed to those observed, a buckled form (along the c-axis in both the a and b directions) of the anatase phase is erroneously predicted to be more

Fig. 10. Metastable (IP1) structures found when targeting anatase: (a) anatase buckled in the a-axis along c-axis, (b) Cu2O, (c) alternating phases, (d–e) two defective Cu2O and (f) Cu2O with square-planar Ti.

91

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

 with a = 4.172 Å. On Cu2O structure adopted the space group Pn3m, average, the L-EA took 32 min, i.e. about eight times longer than the runs for rutile, and evaluated 322,830 candidates (72,563 times when derivatives were also calculated). We estimated that an L-EA run for 1 population of 20 candidates is equivalent to a D-EA run that generates 450 populations of 100 candidates. On average, a D-EA run for the anatase system took just under 3½ min. For D-EA (CONV) and D-EA (RFO), from the final populations, 18 and 28 unique structures were found, which increased to 35 and 46 when we searched the landscape defined by IP2. Upon a closer examination, we found that many of these structures were obtained because the optimiser in stage 2 failed to converge properly. In an attempt to reduce the number of unique structures (and to remove all saddle point structures), which would need to be examined, BFGS and RFO procedures were replaced with that of SWITCH and RFO2. The results are promising: 10 and 6 unique candidates on the landscape of IP1, and 23 and 20 unique candidates for IP2 were obtained. The new distributions of candidates structures found are shown within the parentheses in Table 4. For D-EA (SWITCH, IP1), although 93 candidates adopted the buckled anatase structure, only 8 of the 50 best candidates relaxed to the target – again emphasising that the better candidates produced from a D-EA must be examined rather than just the best candidate. Other low energy structures that were generated are shown in Fig. 10; for example, the structure in Fig 10c is composed of two alternating phases that are possible in half the unit cell of the target; Cu2O and another where Ti is 5-fold coordinated as in Fig. 11f. Another example is that in Fig. 10f, where the cations adopt the Cu2O structure and the anions are arranged so that the corner sharing TiO4 structural units are square planar (rather than tetrahedra). The lowest energy structures found for IP2 are shown in Fig. 11. The first of these structures resembles that of anatase but has two oxygen ions that are two- rather than threefold coordinated (highlighted by use of a larger sphere), and the second structure is that of TiO2 (II), as previously found [20,46] when targeting the brookite phase. The distributions of structures generated, shown in the lower half of Table 4, clearly show that the anatase phase is readily found. The success rate of the D-EA (how often the target was found at least once at the end of a run) was 76% (SWITCH) and 96% (RFO2) for IP1 and 88% (SWITCH) and 98%

(RFO2) for IP2, and thus not as good as that found for the L-EA (always 100%). Although probably not a fair comparison, based on CPU time used, the EA approaches appear to be more efficient than the simulated annealing approach with an 8% success rate [46]. Ramsdellite (Z = 4): On average, the L-EA took 33 (IP1) and 29 (IP2) min, i.e. about the same time as the runs for anatase. Thus, we used the same number of cycles, as used when targeting anatase, for the D-EA. The distributions of structures generated when targeting Ramsdellite are shown in Table 5. Both L-EA and D-EA readily found the global minimum. Unfortunately, even with the cell parameters constrained to those of Ramsdellite, the rutile phase was found to be the structure with the lowest energy for both IP1 and IP2. In fact for IP1, only four out of the 500 candidates in the final populations of the L-EA runs were that of the target, which reduced to zero and one for D-EA (SWITCH) and D-EA (RFO2), respectively. On each run, with a sufficient number of cycles, the search should find up to 10 different structures. Although four structures had a lower energy than that of the target, it is the dominant basin containing the rutile phase that causes the failure of the search to find the Ramsdellite phase readily with the fixed number of cycles employed. We note that the difference in energies, rutile is over 2 eV lower in energy than all other structures found when targeting Ramsdellite. The energy landscape defined by IP2 is less problematic, as the L-EA is readily able to generate the target with a success rate of 96%, and also found another porous structure shown in Fig. 12d. Comparing L-EA with D-EA, with the limited number of cycles, the former also does better at finding more key (low energy) structures. Brookite (Z = 8): The average time for the completion of 10 cycles for the L-EA was 73 min for IP1 and 66 min for IP2, and as all runs obtained the GM after just one cycle, we chose to use 1000 cycles for our D-EA runs with Z = 8 (which took, each run on average, just over 6½ min). The distribution of candidate structures found (see Table 6) shows that the D-EA also readily finds the GM, however, as with Ramsdellite, the GM is not the target structure of brookite, but, in fact, the TiO2 (II) phase shown in Figs. 13a and b. Viewed from one direction (as shown in Fig. 13b and the upper view in Fig. 1d), brookite and the TiO2 (II) phase appear to have the same structure. However, as seen in Fig. 13a, the GM structure, like the Cu2O structure found within the anatase unit cell, has a periodicity that is half that of the longest length of the unit cell for brookite (as seen in Fig. 9e). As the two shorter cell parameters are similar but not equal, there are also further strained versions of each phase, so that when relaxed under constant pressure, these are essentially the same. The success rates of finding the target structure for IP1 are 100% (L-EA), 54% (D-EA(SWITCH)) and 48% (D-EA(RFO2)), but the D-EA figures improve to 82% (SWITCH) and 72% (RFO2) if the strained version of the target is counted.

Table 5 The distribution of candidate structures generated using an evolutionary algorithm, followed by either SWITCH or RFO2 (see text), where Ramsdellite (Z = 4) is the target phase

Fig. 11. Metastable (IP2) structures found when targeting anatase. Note that, for clarity, in (a) we have used a larger ball to indicate oxygen anions that are undercoordinated, and we have redefined the unit cell for the structure (b), which has the space group Cmc21.

Polymorph

L-EA (SW)

L-EA (RFO2)

D-EA (SW)

D-EA (RFO2)

Energy (eV/uc)

Rutile (Fig 12a) Layered (Fig. 12b) Low energy structures Target (Fig. 1f) Porous phase (Fig. 12d) Other structure

469 26 0 4 1 0

469 26 0 4 1 0

487 6 2 (2) 0 4 1

498 0 1 0 0

448.3201 446.1659 < 446.0015 445.3084 444.8659 430.5990

Rutile Target (Fig. 1e) Porous phase (Fig. 12d) Other structures

357 134 9 0

357 134 9 0

495 1 0 4 (2)

491 0 0 9 (3)

157.9287 156.8220 156.5734 < 155.8143

Numbers in parentheses relate to the number of different structures.

92

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

Fig. 12. Structures found when targeting Ramsdellite: (a) rutile, (b) layered phase (relaxed using IP1), (c) Ramsdellite and (d) porous framework, where half the pores are twice as wide as that in Ramsdellite.

Table 6 The distribution of candidate structures generated using an evolutionary and genetic algorithm, followed by either SWITCH or RFO2 (see text), where brookite (Z = 8) is the target phase

Table 7 The distribution of candidate structures generated using an evolutionary algorithm, followed by either SWITCH or RFO2 (see text), where Hollandite (Z = 8) is the target phase

Polymorph

L-EA (SW)

L-EA (RFO2)

D-EA (SW)

D-EA (RFO2)

Energy (eV/uc)

Polymorph

L-EA (SW)

L-EA (RFO2)

D-EA (SW)

D-EA (RFO2)

Energy (eV/uc)

TiO2 (II) (Fig. 13a and b) Target (Fig. 1d) Target (a M b) TiO2 (II) (a M b) (Fig. 13c) Other structures

213 132 58 60 31 6 (2)

213 132 58 60 31 6 (2)

245 37 40 64 29 85 (62)

235 41 44 62 26 92 (67)

899.6037 898.6637 898.3873 897.7961 896.9099 < 896.6038

Layered (Fig 14a) Target (Fig. 1h) Double layer (Fig. 14c) Layered (Fig. 14b) Porous (Fig. 14e) Other structure

115 77 153 71 55 29

115 77 153 71 55 29

5 6 20 3 10 456

3 8 15 3 8 463

894.1369 894.1069 894.0942 893.3668 893.2951 < 892.8156

TiO2 (II) (Fig. 13a and b) Target (Fig. 1d) Target (a M b) TiO2 (II) (a M b) Other structures

186 132 22 50 110 (6)

186 132 22 50 110 (60)

72 26 21 28 353 (*)

64 24 24 29 359 (*)

317.3131 316.9384 316.1744 315.3733 < 315.3605

Target (Fig. 1g) Porous (Fig. 14e) Porous (Fig. 14d) Porous (Fig. 12d) Other structures

11 62 30 51 346

11 62 30 51 346

1 19 2 9 469

0 21 2 4 473

313.2744 312.3811 311.8256 311.6100 < 311.3921

Numbers in parentheses relate to the number of different structures.

Fig. 13. Structures found when targeting brookite: (a) brookite, (b) TiO2 (II) phase, and (c–d) two higher energy phases.

The success rates of finding the target structure for IP2 are 100% (LEA), 44% (D-EA(SWITCH)) and 38% (D-EA(RFO2)), and, as with IP1, the D-EA figures improve to 62% (SWITCH) and 64% (RFO2) if the strained version of the target is counted. As with the results when targeting anatase, these success rates are much better than those obtained previously by Freeman et al. [46]. Hollandite (Z = 8): The average time for the completion of 10 cycles of the L-EA was 242 min for IP1, and 212 min for IP2. Unlike the previous systems, where a perfect success rate at generating the GM after just one cycle was found, the microporous framework of Hollandite proved much more difficult. For IP1, a 100% success rate for generating the GM structure was obtained only after an average of 2.9 cycles and from the final distribution of candidate structures (Table 7), the target, which is the second lowest energy structure, was generated with a 98% success rate. The success rate

dropped dramatically, to 12% (SWITCH) and 16% (RFO2), when the D-LA was employed for the same number of cycles used when targeting brookite, i.e. 1000 cycles. However, more D-EA cycles should have been used for a fairer comparison, as the L-EA Hollandite runs took much longer than the L-EA brookite runs due to the greater CPU cost of the maturing process (relaxing newly created dense candidate structures is less demanding than that of porous structures), whereas average time for D-EA brookite and D-EA Hollandite runs were similar. For IP2, the target structure of Hollandite was also found to be the GM. Rather than favouring 1D or layered structures, as shown in Fig. 14a–c, as well as the target, Ramsdellite (within a 2  1  1 supercell, see Fig. 14d) and other feasible, metastable structures (see Figs. 12d and 14e) on the IP2 energy landscape was also (and more readily) predicted. There were also an abundance of higher energy structures, symptomatic of both the larger energy landscape (Z = 8) and the low density, or due to employing too few cycles. The success rate of finding the target after 10 L-EA cycles or 1000 D-EA cycles was 20%, 2% (SWITCH) and 0% (RFO2), respectively. Another approach to increasing the success rate of this microporous material would be to employ exclusion zones [37,54,55]. Number of line searches (m): When considering how to reduce the computational cost of the L-EA, one key parameter is the maximum number of line searches allowed when maturing, or relaxing, each new candidate structure. Obviously for smaller values of m, less computational effort is required to evaluate each candidate, where in the limit m = 0 the L-EA becomes identical to the D-EA. Already we have seen that the L-EA is more successful than the D-EA, so m = 0 is not favourable. However, just a few lines searches, when the gradient norm can typically be large and initial searches tend to separate ions within regions that are erroneously too densely packed, maybe beneficial whereas any more line searches

93

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

Fig. 14. Structures found, other than the target, when targeting Hollandite using (a–c) IP1 and (d) IP2.

Table 8 The distribution of candidate structures generated using L-EA, where anatase (Z = 4) is the target phase Polymorph/m

1000

50

25

5

1

Energy (eV/uc)

Target

148 (98%) 310 13 26

459 (100%) 41

275 (94%) 225 0 0

81 (46%) 85 (54%) 381 23 14 1 0

116 (86%) 151 (94%) 354 319 17 9 4 18 20 63

448.7552

173 (90%) 186 (90%) 80 77 41 30 50 65 40 29 27 22

132 (90%) 187 (98%) 49 86 43 28 63 63 69 27 48 18 95

Cu2O (Fig. 10b) Fig. 10d Fig. 10e Fig. 10f Other Anatase Fig. 11a Fig. 11b Fig. 11c Fig. 11d and e Fig. 11f Other

3 165 (98%) 28 30 62 28 30 95

0 500 (100%)

498 (100%) 2

448.4814 447.0540 447.0427 446.7681 < 444.08410 157.9790 155.6858 154.8568 154.7481 154.6813 154.6784 < 154.6166

could be wasteful. To remove small structural refinements, a more suitable tolerance for the accuracy in atomic coordinates (i.e. less accuracy than that used above) could be chosen. Above, the criterion used for determining whether two candidate structures are

within the same energy basin, by comparing lattice energies, is less reliable once candidate structures are no longer local minimum structures. Alternatively, the tolerance for the accuracy found for the lattice energy could be changed; however, we chose to examine the effect of varying m. In the L-EA results discussed above we effectively used an unlimited number of line searches (m = 1000). In Table 8 we present the distribution of structures when targeting anatase for when m = 1, 5, 25, 50 and 1000. The average time for each L-EA run was 21 s, 61 s, 239 s, 452 s and 1932 s for IP1, and 22 s, 63 s, 245 s, 467 s, 2043 s for IP2. Typically for m < 6, not all the energies of the candidate structures in the 10th and final population are near enough to LM so that distribution is dependent upon the chosen method used in the second stage. For both landscapes, one expects the proportion of the final population of candidate structures, which adopt the target structure, P, to increase with m. However, some populations will have more than one copy of the target even though we try to kill off duplicates during stage 1, which is evident in the larger diversity for m = 1000 compared to 50 and 25. Therefore, particularly as more duplicates will survive for smaller m, P may be a misleading measure of the success. Hence, we also report the success rate, as a percentage of final populations that contain the target structure. Ignoring data for m < 6, as these are dependent upon stage

Table 9 The distribution of candidate structures generated using D-EA (SWITCH), when anatase (Z = 4) was the target phase Polymorph/n

50,000

450

200

100

10

Energy (eV/uc)

Target Cu2O (Fig. 10b) Fig. 10d Fig. 10e Fig. 10f Other

127 (78%) 348 10 14 1 0

93 (76%) 373 17 11 1 5

113 (88%) 347 14 21

119 (98%) 323 12 4

5

100 (86%) 344 17 22 5 12

448.7552 448.4814 447.0540 447.0427 446.7681 < 444.08410

Anatase Fig. 11a Fig. 11b Fig. 11c Fig. 11d and e Fig. 11f Other

185 (90%) 49 23 57 48 66 72

145 (88%) 50 26 55 50 57 117

142 (94%) 69 28 56 55 36

165 (98%) 57 32 60 38 47

177 (100%) 60 26 54 49 34

157.9790 155.6858 154.8568 154.7481 154.6813 154.6784 < 154.6166

94

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

Table 10 The distribution of candidate structures generated using D-EA (RFO2), when anatase (Z = 4) was the target phase Polymorph

50,000

450

200

100

10

Energy (eV/uc)

Anatase buckled Cu2O (Fig. 10b) Fig. 10d Fig. 10e Other

85 (63%) 385 10 20 0

148 (98%) 310 13 26 3

140 (94%) 318 11 26 5

135 (90%) 314 9 34 8

167 (100%) 297 11 20

448.7552 448.4814 447.0540 447.0427 < 444.08410

Anatase Fig. 11a Fig. 11b Fig. 11c Fig. 11d and e Fig. 11f Other

195 (100%) 36 26 70 58 35 80

165 (98%) 28 30 62 28 30 95

199 (100%) 76 20 53 35 19

200 (100%) 71 22 61 36 22

189 (100%) 97 18 49 30 25

157.9790 155.6858 154.8568 154.7481 154.6813 154.6784 < 154.6166

Fig. 15. 1D landscape, E(x), containing two minima and marked evenly along x. The dotted lines are at the height of the local minimum and the lowest mark within the basin containing the global minimum. The thirteen energy marks are the random structures within the first population. Most (six) of the candidates within the LM basin, if a local optimisation is not performed, are deemed better than those in the GM basin and are therefore more likely to procreate. Moreover, if natural selection allows only six candidates to survive then all the marks within the GM basin will be lost in the second population. If elimination of duplicates was extended so that two candidates must be 3l apart then four of these six LM candidates would survive – thus preventing the population’s diversity from collapsing (to solutions within just one basin) and allowing the GM candidates to survive long enough to procreate a candidate further down in the GM basin.

2, the L-EA is very successful at generating the target structure. Although for IP1 the success rates increase to 100% when m increases to 50, as we find a 98% success rate for m = 1000, we can conclude that our statistical error is of the order of 2%. Number of D-EA cycles (n): In Tables 9 and 10, we present the distribution of structures when targeting anatase for n = 10, 100, 200, 450 and 50,000. The average time for each D-EA run was approximately 0.47n s, with the longest sets taking approximately 6½ h per run. For the much longer runs, with the exception of one set, the success rate, as expected, improves, but we also find that the diversity of structures within the final populations drops and that results still depend on the choice of optimiser in the second stage; SWITCH proves to work much better than straight RFO. Surprisingly, the success rate initially falls as more D-EA cycles are used. From the results discussed above, it has been found that the anatase system, with Z = 4, is trivial enough so that merely quenching (relaxing) several randomly generated structures may generate the target structure. However, without the ageing phase, which is used within the L-EA, the landscape favours the discovery of the Cu2O structure. Fig. 15 shows an landscape, where an LM is initially favoured over a GM when using a D-EA, which indeed probably explains the cause of only an 8% success rate of the simulated annealing approach [46]. Notice that points high up the walls of the basin containing the GM are less likely to survive and procreate than points low down within the basin containing the LM. However, this loss can be avoided if an L-EA is used instead of the D-EA, where linear minimisations take points high up on the basin walls down to the GM. With a much-increased number of D-EA cycles, as the population evolves, the D-EA will find the GM, or target, provided the diversity of the population does not collapse. One way to ensure this was discussed by Lee [56].

4. Conclusion Comparing the efficiency of the two approaches, the Darwinian-Evolutionary Algorithm (D-EA) relies on the simple fact that a single point on the energy landscape is far faster to evaluate than performing a structural relaxation. The efficiency of the Lamarckian-EA (L-EA), where every new candidate structure is immediately relaxed, results from searching the relaxed landscape (where every point has the value of its local minimum), which is much less complex and has many fewer possible states to sample than the D-EA. Moreover, candidates within the same local region (e.g. basin) of the energy landscape are much easier to find (and therefore eliminate) in the L-EA approach; the efficiency of any EA will deteriorate dramatically if the population has a low diversity. In the work reported here, the L-EA proved more successful and efficient at finding the chosen crystalline target structures of titania, in agreement with experience, on other systems, of the other delegates at the symposium on genetic algorithms. We have used two sets of interatomic potential parameters, and thus defined two different energy landscapes to search: the first, which is based on formal charges and the second based on partial charges. The local minima on the latter energy landscape correspond to structures that are in better agreement with the observed polymorphs of titania. For example, using the first set of potential parameters, one of the Ti-O bonds in each of the TiO6 octahedra in the framework structures is erroneously stretched; in Hollandite the Ti–O bond stretches to 2.89 Å to create parallel 1D tubes. However, the formal charge model is better at reproducing the relative energetics of the meta-stable structures and has fewer non-targeted, high-energy structures, although, on the negative side, we found that the formal charge model also erroneously predicted

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

the GM structure to have titanium cations that are five-, rather than sixfold coordinated by the oxygen anions. Interestingly, although the unit cells had an identical number of constituent ions, the less dense phase of Hollandite was much harder to generate than that of brookite. Compared to typical unsolved structures, where the number of unknown variables is larger, most of these phases proved to be too trivial as quenching random structures may generate them. Moreover, earlier development work on D-EAs that were applied to searching the landscape of these trivial systems should not be discarded because the target structure momentarily was lost. For future testing of the performance of the two approaches, larger supercells of these systems will be used, thus we are unlikely to obtain the currently observed high success rate in the earlier cycles of the D-EA (before the diversity of the population falls). We will further investigate how the success rate of finding the target varies with the parameters of the maturing phase of the L-EA, but with the constraint that the total CPU time is fixed. Acknowledgements We would like to thank Alexey Sokol for valuable comments, EPSRC for financial support (Portfolio Grant EP/D504872/1), Acceryls for the provision of visualisation software, colleagues, and particularly the symposium organisers Wojciech Paszkowicz and Kenneth Harris, all of whom made the Genetic Algorithm Symposium a success. References [1] [2] [3] [4] [5]

[6] [7] [8] [9] [10] [11] [12] [13] [14]

N.R. Mandal, Sci. Technol. Weld. Join. 4 (1999) 290. A. Kumar, G.G. Roy, Metallurg. Mater. Trans. B 36 (2005) 901. B. Bandow, B. Hartke, J. Phys. Chem. A 110 (2006) 5809. R.M. Briggs, C.V. Ciobanu, Phys. Rev. B 75 (2007) 195415. G.M. Day, W.D.S. Motherwell, H.L. Ammon, S.X.M. Boerrigter, R.G. Della Valle, E. Venuti, A. Dzyabchenko, J.D. Dunitz, B. Schweizer, B.P. van Eijck, P. Erk, J.C. Facelli, V.E. Bazterra, M.B. Ferraro, D.W.M. Hofmann, F.J.J. Leusen, C. Liang, C.C. Pantelides, P.G. Karamertzanis, S.L. Price, T.C. Lewis, H. Nowell, A. Torrisi, H.A. Scheraga, Y.A. Arnautova, M.U. Schmidt, P. Verwer, Acta Cryst. B 61 (2005) 511. R.L. Johnston, Dalton Trans. (2003) 4193. E. Flikkema, S.T. Bromley, J. Phys. Chem. B 108 (2004) 9638. M.A. Millar, D.J. Wales, J. Phys. Chem. B 109 (2005) 23109. S.M. Woodley, in: R.L. Johnston (Ed.), Structure and Bonding, vol. 110, Springer-Verlag, Heidelberg, 2004, p. 95. J.C. Schön, M. Jansen, Z. Kristallogr. 216 (2001) 307. A. Corma, M. Moliner, J.M. Serra, P. Serna, M.J. Diaz-Cabanas, L.A. Baumes, Chem. Mater. 18 (2006) 3287. Z.G. Pan, E.Y. Cheung, K.D.M. Harris, E.C. Constable, C.E. Housecroft, Powder Diff. 20 (2005) 345. A.J. Hanson, E.Y. Cheung, K.D.M. Harris, J. Phys. Chem. B 111 (2007) 6349. A. Dolbecq, C. Mellot-Draznieks, P. Mialane, R. Marrot, G. Ferey, F. Secheresse, Eur. J. Inorg. Chem. 15 (2005) 3009.

95

[15] G. Ferey, C. Mellot-Draznieks, C. Serre, F. Millange, Acc. Chem. Res. 38 (2005) 217. [16] K.D.M. Harris, M. Tremayne, P. Lightfoot, P.G. Bruce, J. Am. Chem. Soc. 116 (1994) 3543. [17] K.D.M. Harris, S. Habershon, E.Y. Cheung, R.L. Johnston, Z. Kristallogr. 219 (2004) 838. [18] S. Habershon, Z. Zhou, G.W. Turner, B.M. Kariuki, E.Y. Cheung, A.J. Hanson, E. Tedesco, R.L. Johnston, K.D.M. Harris, EAGER, Cardiff University and University of Birmingham. [19] W. Paszkowicz, Mater. Sci. Forum 228 (1996) 19. [20] S.M. Woodley, P.D. Battle, J.D. Gale, C.R.A. Catlow, Phys. Chem. Chem. Phys. 1 (1999) 2535. [21] C.M. Freeman, C.R.A. Catlow, J. Chem. Soc., Chem. Commun. 2 (1992) 89. [22] L. Reinaudi, E.P.M. Leiva, R.E. Carbonia, J. Chem. Soc. Dalton 23 (2000) 4258. [23] C.J. Howard, T.M. Sabine, F. Dickson, Acta Cryst. B 47 (1991) 462. [24] E.P. Meagher, G.A. Lager, Can. Mineral. 27 (1979) 77. [25] L.D. Noailles, C.S. Johnson, J.T. Vaughey, M.M. Thackeray, J. Power Sources 81– 82 (1999) 259. [26] Y. Takahashi, N. Kijima, J. Akimoto, Chem. Mater. 18 (2006) 748. [27] J.D. Gale, J. Chem. Soc., Faraday Trans. 93 (1997) 629. [28] D.A. Coley, An Introduction to Genetic Algorithms for Scientists and Engineers, World Scientific, Singapore, 1999. [29] T.S. Bush, C.R.A. Catlow, P.D. Battle, J. Mater. Chem. 5 (1995) 1269. [30] J.C. Schön, M. Jansen, Comput. Mater. Sci. 4 (1995) 43. [31] L. Reinaudi, E.P.M. Leiva, R.E. Carbonio, J. Chem. Soc. (2000) 4258. [32] M.W. Deem, J.M. Newsam, J. Am. Chem. Soc. 114 (1992) 7189. [33] B. Hartke, in: R.L. Johnston (Ed.), Structure and Bonding, vol. 110, SpringerVerlag, Heidelberg, 2004, p. 35. [34] D.J. Wales, J.P.K. Doyle, J. Phys. Chem. A 101 (1997) 5111. [35] S.M. Woodley, A.A. Sokol, C.R.A. Catlow, Zeit. Anorg. Allge. Chem. 630 (2004) 2343. [36] S. Hamad, C.R.A. Catlow, S.M. Woodley, S. Lago, J.A. Mejías, J. Phys. Chem. B 109 (2005) 15741. [37] S.M. Woodley, Phys. Chem. Chem. Phys. 9 (2007) 1070. [38] G.W. Turner, E. Tedesco, K.D.M. Harris, R.L. Johnston, B.M. Kariuki, Chem. Phys. Lett. 321 (2000) 183. [39] N.L. Abraham, M.I.J. Probert, Phys. Rev. B 73 (2006) 224104. [40] F.C. Chuanga, C.V. Ciobanub, V.B. Shenoyb, C.Z. Wanga, K.M. Ho, Surf. Sci. 573 (2004) L375. [41] N. Dugan, S. Erkoc, Symposium G, EMRS, autumn 2007. [42] S. Datta, F. Pettersson, S. Ganguly, H. Saxen, N. Chakraborti, Isij Inter. 47 (2007) 1195. [43] For example, Molecular mechanics across chemistry, A.K. Rappé, C.J. Casewit, University Science Books, Sausalito, California, 1997. [44] P.P. Ewald, Ann. Phys. 64 (1921) 253. [45] M.P. Tosi, Solid State Phys. 16 (1964) 1. [46] C.M. Freeman, J.M. Newman, S.M. Levine, C.R.A. Catlow, J. Mater. Chem. 3 (1993) 531. [47] M. Matsui, M. Akaogi, Mol. Simul. 6 (1991) 239. [48] P.M. Oliver, G.W. Watson, E.T. Kelsey, S.C. Parker, J. Mater. Chem. 7 (1997) 563. [49] We have chosen more recent reported observations than that used in the potential fitting procedure. [50] M. Sierka, T.K. Todorova, J. Sauer, S. Kaya, D. Stacchiola, J. Weissenrieder, S. Shaikhutdinov, H.J. Freund, J. Chem. Phys. 126 (2007) 234710. [51] E. Flikkema, S.T. Bromley, Chem. Phys. Lett. 378 (2003) 622. [52] By P. Collet in the opening presentations in the Genetic Algorithms for Beginners workshop held before the start of the 2007 EMRS Fall meeting. [53] A. Banerjee, N. Adams, J. Simons, R. Shepard, J. Phys. Chem. 89 (1985) 52. [54] S.M. Woodley, P.D. Battle, J.D. Gale, C.R.A. Catlow, Phys. Chem. Chem. Phys. 6 (2004) 1815. [55] S.M. Woodley, Phys. Chem. Chem. Phys. 6 (2004) 1823. [56] J. Lee, I.H. Lee, J. Lee, Phys. Rev. Lett. 91 (2003) 080201.