WO2007071095A1 - Method for crystal structure determination - Google Patents

Abstract

The disclosure pertains to a method for determining a crystal structure comprising the steps of initialisation, local optimisation, selection, production of the next generation with the help of variation operators wherein as variation operators genetic crossover, mutation and permutation are used, wherein the simulation is terminated after it converges to a global minimum.

Description

SPECIFICATION

TITLE Method for Crystal Structure Determination

TECHNICAL FIELD

The disclosure pertains to a method for determining a crystal structure comprising the steps of initialisation, local optimisation, selection, production of a next generation, wherein the simulation is terminated after it converges to a global minimum.

BACKGROUND OF THE INVENTION

The stable crystal structure is arguably the greatest bearer of information on a given material at given conditions. By crystal structure we mean shape of the unit cell (lattice) and the position of each atom within the unit cell. While in many cases it is possible to solve crystal structure from experimental data, the possibility to predict it purely on the basis of theory is crucially important for several reasons:

1) When experimental data are of insufficient quality (poor powders, small amounts of the sample, especially at high pressures and temperatures) theory provides the last resort;

2) Theory is the only way of investigating matter under extreme conditions that cannot be studied with today's experimental techniques, e.g. at ultrahigh pressures;

3) The ability to predict crystal structures will open up new ways for materials design, providing the complete information about a series of compounds - their crystal structures and properties - entirely by computation and thus guiding experimental efforts to synthesise the most promising phases.

Crystal structure prediction is a major unsolved scientific problem, which attracts significant interest. The computational requirements in this respect are exceedingly demanding. To get some feeling of the number of possible structures, let us consider a simplified case of a fixed cubic cell with volume V, within which one has to position N identical atoms. For further simplification let us assume that atoms can only take discrete positions on the nodes of a grid with a certain resolution. This discretisation makes the number of combinations of atomic coordinates finite. If the resolution is chosen to be a significant fraction of the characteristic bond length (e.g., resolution = 1 A), the number of combinations would be a reasonable estimate of the number of local minima of the free energy. If there is more than one type of atoms, the number of different structures significantly increases. Assuming a typical atomic volume ~10 A³, and taking into account Stirling's formula, the number of possible structures for an element A (compound AB) is 10¹¹ (10¹⁴) for a system with 10 atoms in the unit cell, 10²⁵ (10³⁰) for a system with 20 atoms in the cell, and 10³⁹ (10⁴⁷) for the case of 30 atoms in the unit cell. One can see that these numbers are enormous and practically impossible to deal with even for small systems with the total number of atoms N ~ 10. Even worse, these numbers increase exponentially with N. It is clear then, that point-by- point exploration of the free energy surface (going through all possible structures) is not a viable solution, except for the simplest systems with ~l-5 atoms in the unit cell.

However, one does not have to explore the entire free energy surface in order to locate the global minimum - it would suffice to explore just the most promising regions of that surface, and there are a number of methods that can do this job. Clearly, to solve the crystal structure prediction problem one either has to start already in a good region of configuration space (so that no effort is wasted on sampling bad regions), or use a "self- improving" method that locates, step by step, the best structures - and in doing so it should ideally reduce the search space exponentially as a function of time in order to counteract the "exponential wall". The first group of methods includes meta-dynamics, simulated annealing, "basin hopping" and "minima hopping" methods. The second group essentially includes only evolutionary algorithms. The strength of evolutionary simulations is that such simulations do not require any system specific knowledge (except chemical composition) and are self-improving, i.e. in subsequent generations increasingly good structures are found and used to generate new structures. This allows a 'zooming in' on promising regions. Furthermore, due to the flexible nature of the variation operators, it is very easy to incorporate features of other methods into an evolutionary algorithm.

So several methods have been applied and developed in order to solve this problem, e.g. simulated annealing and meta-dynamics. These methods proceed via neighbourhood search, thus depending on a starting structure. Even for good starting structures - which are usually very hard to come by - these methods are computationally extremely expensive, while yielding modest success. Evolutionary algorithms have previously been used, but in a very different implementation which yielded low success rates even for simple structures and could only work for fixed experimental cell parameters. The problems of that algorithm are the lack of local optimisation and the use of bit-strings to represent structural variables and exchange genetic information; in addition a very simple heuristic fitness function is normally used, which limits applicability of that algorithm.

Further methods for computational crystal structure determination but with different methodological approaches are for example protected in EP 0 943 131 Bl, EP 0 992 011 B 1 or EP 0 592 421 B 1 and the corresponding patent families worldwide.

SUMMARY OF THE INVENTION

The method discussed here is completely different from those discussed above. Its design allows it to be independent from prior assumptions/knowledge (e.g. experimental results, starting structures) while having a much higher success rate and being much faster. It is capable of finding the most stable structure and discovering a large set of meta-stable structures on the way. Furthermore, the generated structures and their free energies yield information on the chemical regime of the system at given conditions. The most important achievement done here is the purely ab initio crystal structure prediction, allowing one to design new materials on the computer. According to many experts, crystal structure prediction is an insoluble task - this disclosure shows that this view must be changed.

The most non-trivial aspects are the hard constraints (which greatly decrease the search space and help local optimisation), cell-vectors cross-over, and the use of planar slices for genetic crossover with addition/removal of atoms to maintain correct chemical composition (without this trick it should be very difficult to study any system containing more than one type of atoms, and to study large systems), and random slice shifts introducing benign mutations. Results of this disclosure include: user-friendly and stable USPEX code and a number of new crystalline materials predicted.

Brief technical description of the USPEX (Universal Structure Predictor: Evolutionary Xtallography) method:

This is an evolutionary method for crystal structure prediction. The simulation is targeted to provide the stable crystal structure and a set of robust meta-stable structures at given pressure-temperature conditions.

1. The method uses ab initio free energy as the fitness function.

2. Structural parameters (cell vectors, atomic coordinates) are represented by real numbers rather than bit-strings. 3. Local optimisation is performed for each generated structure.

4. Calculations require only the number of atoms in the unit cell (for each atomic species), cell parameters are found as a result of the simulation.

As a result, one can predict crystal structures at any given pressure-temperature conditions without having any experimental information. This allows one to design new materials computationally, as well as predict hitherto unknown crystalline phases. When cell parameters are known from experiment, they can be fixed, thus simplifying the simulation.

Technically, there are several stages in the process:

Initialisation: the first population is preferably created randomly. However, all structures must satisfy a set of hard constraints: for example, the interatomic distances should be greater than 0.5 A, cell angles between 60 and 120 degrees, distances between cell coordinate planes should be greater than 2 A.

Local optimisation: performed on every individual at each generation. During local optimisation, the k-point grid changes in accordance with cell changes. This enables strict comparability between all obtained free energies while keeping the computational costs low. Selection: at each generation, only a selected percentage of the best structures (e.g., 50%) are used to produce the next generation while the worst structures are discarded.

Production of the next generation with the help of variation operators. By variation operators we mean genetic crossover (heredity), mutation and permutation. The probability that a certain structure is chosen for such an operation is derived from its rank of their fitness. A certain fraction of the new generation is created by cross-over, where we represent atomic positions by real numbers, not by the often used, but unphysical, bit-strings, and new structures are produced by matching planar slices (chosen in random directions and with random positions and shifts in the unit cell) of the parent structures. When slices put together do not contain the correct number of atoms, some atoms are added or deleted to maintain the correct chemical composition. Large number of random factors incorporated in cross-over can be viewed as "coordinate mutation". Importantly, we always use the fractional coordinates - unlike Cartesian coordinates, these have good invariance properties under cell shape transformations. Similar cross-over is done for the lattice vectors matrix elements (for this matrix we always use the upper-triangular form in order to avoid rotations of the cell). A certain fraction of the new generation is created by permutation (switching identities of two or more atoms in a structure) and mutation (random Gaussian change of the cell vectors). The best structure in a given generation always survives, mates and competes in the next generation.

The simulation is terminated after it converges to a global minimum. For systems with up to ~ 12 atoms in the cell the global minimum is often found within the first few generations, for larger systems this usually takes ~ 20 generations.

One of the objective problems underlying the present invention is therefore to provide an improved method, particularly computer-based, for the determination of crystal structures. The method for determining a crystal structure comprises the steps of initialisation, wherein the population is either created by seeding the first generation with complete structures and if necessary filling up the ranks with random individuals, wherein the structures can be chosen from any external source or can be taken from any other runs of the method or is created randomly, application of hard constraints, all structures must satisfy a set of hard constraints preferentially selected from the group of: interatomic distances, which preferably should be greater than 0.5 A; cell angles, which preferably should be between 60 and 120 degrees; distances between cell coordinate planes, which preferably should be greater than 2 A; local optimisation, wherein this step is performed on each individual at each generation and wherein first the atoms are relaxed, keeping the lattice fixed and then both lattice and atoms are relaxed simultaneously, and wherein during local optimisation, the k-point grid changes in accordance with cell changes; selection, wherein at each generation, only a selected percentage of the best structures, preferably in the range of 50%, are used to produce the next generation while the worst structures are discarded, wherein as an evaluation function preferably the ab initio free energy of the locally optimized individual is used as fitness value; production of the next generation with the help of variation operators wherein as variation operators genetic crossover, mutation and permutation is used and wherein the probability that a certain structure will participate in a given operation is derived from its rank of their fitness, wherein for this cross-over operation atomic positions are represented by real numbers, and new structures are produced by matching planar slices, wherein further preferably fractional coordinates are used, and wherein similar cross-over is done for the lattice vectors matrix elements, wherein preferably for this matrix the upper-triangular form is used, wherein a certain fraction of the new generation is created by permutation, preferably by switching identities of two or more atoms in a structure, and mutation, preferably by random Gaussian change of the cell vectors, and wherein the best structure in a given generation always survives, mates and competes in the next generation; wherein the simulation is terminated after it converges to a global minimum.

Above reference is made to structure, meaning the crystal structure. Reference is also made to atoms within that structure. Of course this concept applies for any kind or type of atoms, and independent of whether there is e.g. a specific type of bonding interaction between the atoms. As a matter of fact the proposed method also can be used for molecular systems or partial molecules (fragments), even for large molecular systems including large numbers of atoms which are bonded by chemical bonds or other types of interactions.

If in the following reference is made to atoms specifically, this can also include molecules and/or fragments.

For molecular systems, instead of - or besides - atoms, correspondingly whole molecules or fragments can be treated. In this case it is possible to regard every molecule/fragment as an 'atom' in the above procedures with a few exceptions and/or amendments and/or adaptations: Every molecule/fragment in addition to its position also has an orientation in space (for example expressed as three angles, e.g. three Euler angles) which preferably is be carried over from parents to offspring;

Two extra operators can additionally and beneficially be incorporated in case of treatment of molecules/fragments: Rotation of molecules/fragments (rotational mutation) and inversion of molecules/fragments (inversion mutation - only useful where the inversion gives a different molecule from the original). These possible additional mutations can be treated like any mutation, by applying them to some selected structures. The interatomic distance constraint, which in case of molecules/fragments is actually an intermolecular/interfragmental distance, needs to be bigger than an actual interatomic distance constraint and in addition to that it might be useful take into account the molecule's shape, such as to avoid overlap of different molecules/fragments. So the intermolecular distances to be taken into account for the constraints depend on molecular size and shape, however preferably such that no overlaps between molecules are possible (interatomic distances still need to fulfil minimal distances for computational stability).

For big, flexible molecules/fragments, some internal degrees of freedom (of every molecule/fragment), like e.g. torsion angles, can be to be carried over from parents to offspring. For very large and very flexible molecules/fragments, an extra operator can be useful in this case: molecular/fragmental distortion, which includes torsional angles, changes in bond angles or angles between fragments and the like.

The matching planar slices can be chosen in random directions and with random positions. Particularly the matching planar slices can be chosen with shifts in the unit cell, of the parent structures. The actual positioning of the respective reference system of the cells of the two parents can therefore be chosen independently for each parent when crossing the parents.

The hard constraints do not have to be regarded as hard constraints which simply lead to an exclusion if a certain criterion is not met and to continuation if the criterion is met. Individuals violating certain or all hard constraints can e.g. not be discarded but a fitness penalty for the violating individual can be applied.

Further embodiments of the present invention are outlined in the dependent claims.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

In the following the specifics of the proposed method are presented. The description shall however only serve to illustrate and document how to carry out the invention, it shall not be construed to limit the scope of the invention as defined in the appended claims. The proposed method allows the prediction of crystal structures starting from chemical composition. By chemical composition we mean number of atoms of every type involved (per unit cell). It is thus completely independent from experimental data. Calculations can be performed at any P-T condition.

It is based on an evolutionary algorithm or method, mainly featuring ab initio calculations, local optimization and spatial cross-over.

Each structure characterized by a complete set of values is called an individual. A set of individuals is called a population and each population a generation. The method generates an initial population. Then it uses variation operators and local optimization on selected individuals of this population to produce new individuals. A new population is chosen from the old population and/or the created individuals.

Input - Output:

The minimal input (of course there is further input required, being the parameters of the method and the P-T condition at which the calculations shall be performed) relies purely on theory, being: • Number of atoms of each type

• Initial guess of the lattice volume

• Hard constraints The initial guess of the lattice volume evolves during the run (see below) and the lattice of every individual is scaled so that its volume matches this value. This is not essential but effectively reduces the search space and helps avoiding pathological lattices.

The method is therefore independent from experimental data and if the initial volume is omitted, no prior assumptions/knowledge is required. If lattice parameters are available from experimental data, they can be incorporated. This reduces the search space significantly. In this case manipulations of the lattice are usually prohibited. However, relaxation of the lattice can be performed if necessary (e.g. if the lattice parameters have a considerable error). Every locally optimized individual is stored with the corresponding free energy. Representation:

Since variables in crystal structures are continuous numbers, the method represents every variable by a floating-point value to acquire sufficient resolution and facilitate the use of variation operators. Other representations (e.g bit-strings) are possible, but not practically advantageous. Two types of variables are discriminated: Lattice parameters and atom coordinates.

There are 6 lattice parameters, three angles (α, β, γ) - coded as the fraction of 2*π - and the length of the three lattice vectors. Alternatively lattice parameters are encoded as a three by three matrix in a lower-triangular form. Each atom has three coordinates, coded as the fraction of the length of the corresponding lattice vector.

Initialization:

There are two ways to initialize the method (initialize means to generate the first generation). The first possibility is to generate the individuals contained in the first generation randomly, only accepting the ones fulfilling the hard constraints. The second possibility is to seed the first generation with complete structures and if necessary fill up the ranks with random individuals. Again hard constraints apply. These structures can be chosen from any external source, they can be taken from many other runs of the method with purposefully small populations (with very small populations the method will quickly converge. By doing many runs one often has a nice selection of promising structures, defining a good region to start a big calculation by seeding) or the method can perform runs with smaller unit cells, using those to produce seeds by blowing up the best resulting structures to the desired unit cell size.

The first possibility is the most robust. It assumes nothing and can be universally applied. The second possibility works well if good candidate structures are a priori known, the system has many competitive regions or the scaling from smaller to desired unit cell size works well.

Local optimisation:

By locally optimizing individuals, a superior comparability between individuals and great reduction of iterations can be achieved. The method discussed locally optimizes every individual within every generation and replaces each by its respective local minima, before comparing fitness values and applying variation operators.

For further speed up and precision of the method, local optimization is performed in three steps. First the atoms are relaxed, keeping the lattice fixed. Then both lattice and atoms are relaxed simultaneously (Provided the lattice parameters are treated as variables).

It is not necessary to optimize every individual, but preferably this should be done. Evaluation function:

The method preferably uses the ab initio free energy of the locally optimized individual as fitness value. This is the most accurate and universal measure of quality. For systems where good approximate functions for quality exist, the ab initio free energy can be substituted by them. The advantage of ab initio calculations is accuracy and no required prior assumptions on the system. Surrogate evaluation functions are computationally cheaper but typically work only for a narrow range of system types and otherwise usually fail completely.

Producing a new population:

Candidates for the new population are produced by using variation operators on selected individuals of the old population. The method features three variation operators - crossover, mutation and permutation. The influence of each (relative distribution of usage) can be controlled via input parameters. The produced candidate structures need to fulfil the hard-constraints to be accepted and afterwards the unit cell is scaled so that its volume matches the currently assumed volume (this scaling is used only to produce initial structures. During local optimization cell volume will change.). The assumed volume can change during a run due to parameter control. Once a set of candidates was produced, the new generation can be acquired by different means. The old population is deleted, a specific number of best individuals are carried over from the old generation or the old generation and the candidates are mixed and a new generation is chosen. The choice depends on the current goal. For a more exploratory search (e.g. for randomly initialized runs), the old generation is preferably deleted or the best 1-2 individuals are carried over. For better exploitation (e.g. a run with many good seeds) the mixing is preferably used.

Selection: The method uses a stochastic selection scheme, where individuals with better fitness value ranking have a higher probability of being chosen. There are two probability assignments from which can be chosen, linear and quadratic. In both cases the initial number is the fitness rank counted from the worst up to the best. In the linear case this number is left as it is, in the quadratic case it is squared. These numbers are scaled so that the sum is equal to one. Now each number is the probability of the corresponding individual being chosen for a specific operation.

Furthermore a specified amount of worst individuals are discarded from selection. However, any selection mechanism where better individuals procreate more often has a good chance of working. Cross-over:

Two individuals are selected and used to produce one off spring. This is achieved by taking a fraction of each individual and combining those to a new solution. However, the fraction of each parent should carry as much information of the parent as possible. The information within crystal structures is the relative position of atoms. Thus to conserve information, the fraction of a parent is selected spatially, by taking slabs of each parent. The two slabs are fitted together and the result thereafter made feasible by adjusting the number of atoms of each type to the requirements.

In more detail this works as follows. One lattice dimension is picked, preferably randomly. Before atoms are selected, atoms may be shifted. For each parent and each dimension to be shifted one random number between zero and one is generated and added to the respective coordinate of the atoms. Since the unit cell is periodically repeated, the atoms ending up with a coordinate value above one (outside the unit cell) are adapted so as to lie within the unit cell again. This is achieved by subtracting one where required (in a periodic system, values shifted by a period or a multiple of a period are identical). Original and shifted systems are physically identical. Both for the chosen dimension and the other dimensions the user can specify in how many percent of the cases this should be done. This operation increases diversity and enhances the power of the method, if set correctly. Typically a setting close to 100% for the chosen dimension and close to 10% for the remaining dimensions works well. If shifting was performed, all following operations apply to the shifted individual.

Now a value X between zero and one is determined, again preferably randomly. From the first parent, every atom is taken which has a coordinate value in the chosen dimension between zero and X. From the second parent every atom is taken with a respective coordinate value between X and one. Now the atoms of both parents are taken together.

The total number of atoms of each type is counted and compared with the required number. If there are too many atoms of a type, atoms are removed randomly. If there are too few, the following procedure is repeated until the number is correct. A sector (zero - X or X - one) is picked randomly. The probability of a sector being chosen is either equal to the size of the sector (for sector zero - X the probability is X. For sector X - one it is 1-X) or it is inverse proportional to the atom density in the sector. Then one atom (having coordinates falling in this sector) is chosen randomly from the parent not having provided the atoms within this sector. Alternatively, the probabilities of placing/deleting atoms can be subject to a distribution, making placing/deleting atoms more/less probable, the closer to the boundary of the slab.

The cross-over of lattice parameters is achieved by taking the weighted average of the lattice parameters of both parents, where the weight is chosen randomly.

The power of this cross-over operation lies in the spatial selection of atoms and the immanent diversification. In this way, whole substructures of an individual can be transferred to any location in the off spring. Many details of this procedure can be slightly altered without loosing the functionality of the operation. Random choices can to a large extent be replaced by e.g. fitness proportional choices, fixed procedures or fixed values. Mutation:

One individual is selected. The lattice values are changed by applying a strain matrix, and the atomic positions are mutated by adding a zero-mean Gaussian random variable.

Mutation is important for the diversity within the population. Due to local optimization and diversification within the cross-over operation, mutation of the atomic positions is not important and can be left out. Mutation of the lattice should be present for optimal performance, to prevent a possibly premature convergence towards a certain lattice.

Mutation can be achieved in different ways, e.g. by multiplying each lattice vector by a matrix [ I + εij ] (where strains εy are random zero-mean Gaussian variable between -1 and 1 and re-scaling to its original volume. Permutation:

One individual is selected. Two atoms of different types are exchanged, a variable number of times. Optionally the user can define which atom types are allowed to be interchanged. Preferably one here specifies chemically similar atom types to be allowed.

Obviously, permutation is possible only for systems with different types of atoms. Parameter control:

The method features a parameter control of the ab initio calculations. Depending on the lattice dimensions of an individual, the number of k_poi_nts are calculated and adapted. This greatly enhances accuracy and speed of the calculations, rendering the method much faster overall. The number of k_poi_nts for a given reciprocal space dimension is calculated as follows: kpoints =l/(interPlanarDistance * k_pointsResolution) The resulting k_poi_nts is rounded to next higher integer.

Furthermore adaptation of the originally guessed unit cell volume is possible and enhances speed for systems where this guess is not sufficiently accurate. For each new generation, the new assumed unit cell volume is a weighted (weight is user-defined) average between the old volume and the average volume of the few best individuals of the old generation. Thus depending on the weight, the method adapts more or less fast to the currently most suited volume.

Hard constraints:

Every solution produced by the variation operators is tested against three hard constraints: • Minimal inter-atomic distances, specified for every combination of different atoms.

• Minimal value of the angles α, β and γ

• Minimal lattice vector length.

Solutions violating at least one of these constraints are disqualified.

Hard constraints fulfil two purposes. First the minimal inter-atomic distances must be sufficient to ensure stability of the ah initio calculations (or for the surrogate function if utilized). Second, they reduce the search space and allow for inclusion of system specific knowledge (e.g. if one knows a certain inter-atomic distance to be large in reality, one can set that distance to a large value). Different sets of constraints are possible, however ensuring stability is important and the set of constraints mentioned above work very well in reducing the search space.

Possible is also a routine to make individuals violating the inter-atomic distance constraint feasible by shifting atoms. This is important for large systems, since with increasing number of atoms, the probability of fulfilling this constraint becomes increasingly difficult, even practically impossible. Molecules:

For molecular systems, instead of - or besides - atoms, whole molecules or partial molecules (fragments) can be treated. So every molecule/fragment is regarded as an 'atom' in the above procedures with a few exceptions:

- Every molecule/fragment has in addition to a position also an orientation (three angles) which needs to be carried over from parents to offspring

- Two extra operators need to be incorporated: Rotation of molecules/fragments (rotational mutation) and inversion of molecules/fragments (inversion mutation - only where the inversion gives a different molecule/fragment from the original). These new mutations can be treated like any mutation, by applying them to some selected structures.

- The interatomic distance constraint needs to be bigger (since a molecule/fragment is much bigger than an atom) and take into account its shape, to avoid overlap of different molecules/fragments.

- For big, flexible molecules/fragments, some internal degrees of freedom (of every molecule/fragment, e.g. torsion angles) need to be carried over from parents to offspring. For very large and very flexible molecules/fragments, an extra operator is necessary: molecular/fragmental distortion. Internal degrees of freedom of some molecules/fragments in a structure are changed.

Parameters:

For a system of 20 atoms in the unit cell with unknown lattice and a reasonable guess at the unit cell volume, a reasonable setting for the parameters would look as follows.

Population size: 40-60 Percentage cross-over: 85

Percentage lattice mutation: 10 Percentage permutation: 5 Percentage shifting in chosen direction: 100

Percentage shifting in other directions: 10

Percentage of worst individuals to be discarded from selection: 40

Standard deviation of number of two-atom-exchanges for permutation: 2-3 Standard deviation of lattice strain matrix: 0.7 Resolution for kpoints determination: 0.12 AT¹

Weight for unit cell volume adaptation from generation to generation: 0.5 Number of best individuals to be averaged for unit cell adaptation: 4 Number of best individuals to be carried over from one generation to the next: 1-2 Hard constraints are system specific. Parallelization:

Evaluating the fitness of different solutions within one generation is independent and can thus be processed in parallel. Parallelization is however limited by the fact that in order to generate a new population, all fitness values of the old population need to be known.

Summarized Results:

The method is not only capable of finding the most stable structure, it also discovers a large set of meta-stable structures on the way. Furthermore, the generated structures yield information of the chemical regime at given conditions. This method can thus be used for materials design, both in finding promising structures to synthesize and in giving information on what conditions would be best suited for production.

The method has been tested on various systems: C at 0, 100, 300, 500, 1000, 2000 GPa C at exp.cell of diamond Xe at 100-1000 GPa Si at 10, 14, 20 GPa Fe at 350 GPa TiO₂-anatase at exp.cell SiO₂ at 0 GPa Al₂O₃ at 300 GPa MgSiO₃ at 80,120 GPa MgSiO₃at exp.cell of post-perovskite Si₂N₂O at 0 GPa SrSiN₂ at 0 GPa (NH₂)₂CO at exp.cell

For all these systems (with up to 20 atoms/unit cell) calculation were performed with minimal input or providing the lattice where this is specified and the method has never failed in predicting the stable structure. Additionally, as discussed above, it yielded numerous meta-stable structures.

Furthermore it is extremely efficient. For systems of 20 atoms/unit cell with minimal input, it took around one week running on four processors of the type HP- itaniuni2(1.35GHz).

Already numerous previously unknown new structures have been found by this method:

O at exp.cell

O at 25, 50, 130 GPa

N at 100, 250 GPa F at 100 GPa

MgCO₃ at 150 GPa

CaCO₃ at 50, 80,150 GPa

H at 200, 600 GPa

CO₂ at 50 GPa Up to now, one structure has been verified by experiment. The following modifications are possible:

First, the method also allows as input molecules instead of atoms. For molecular systems this reduces the degrees of freedom significantly, thus improving the power for molecular crystals enormously. This input will be a definition of the molecules and the method will operate on whole molecules (e.g. by using as internal variables a point in the lattice and orientation angles for each molecule). Furthermore, the definition of variable torsion angles is possible and if so as internal variable. Since van-der-Waals interaction is poorly captured by ab-initio calculations, the method will optionally calculate the system under slight pressure - thus simulating the van-der-Waals forces. As mentioned above a different approach towards hard constraints is possible. The quality of fractions used for cross-over can be adapted. This allows to choose not only fractions of good structures, but even good slabs of good structures. Of course, incorporation of this functionality needs to respect diversity - thus not always preferring better slabs. Promising candidates for evaluation functions are Brown's bond valence model for zero pressure and/or Pauling's rules.

Furthermore analysis of likeness of different individuals can be incorporated. This allows optional extension of the fitness including exclusiveness of a given individual. Rare and simultaneously good individuals can be treated with higher priority than slightly better common individuals. Control of the parameters defining the population size, how many individuals should be produced by mutation/permutation and standard deviation of number of swaps per permutation can also be included. The control is e.g. based on free energy distribution within a population (lower differences - smaller population size/bigger mutation and permutation rates. And vice versa.) and on success of permutation (many successful permutations - more swaps per permutation. And vice versa).

Summary and further results:

Our method has been implemented in the USPEX (Universal Structure Predictor: Evolutionary Xtallography) code. The code has a minimal input: the number of atoms of each sort, pressure-temperature conditions and method parameter values. Optionally, calculations can be performed under fixed lattice parameters (if these are known from experiment). Our procedure is the following:

1. The first generation is produced by a random-number generator and/or seeded with structures (only those structures which satisfy the pre-specified hard constraints are allowed) and then locally optimised using powerful conjugate- gradients or steepest-descent optimisers available in many first-principles and atomistic simulation codes. Currently, USPEX can use VASP and SIESTA (Sanchez-Portal D., Ordejon P., Canadell E. (2004). Computing the properties of materials from first principles with SIESTA. Struct. Bond. 113, 103-170) for first-principles simulations, and GULP (Gale J.D. (2005). GULP: Capabilities and prospects. Z. Krist. 220, 552-554) for atomistic simulations. During optimisation with ab initio methods, the k-point grid changes in accordance with cell changes. This enables strict comparability between all obtained free energies while keeping the computational costs low. 2. Among the locally optimised structures, a certain percentage of the worst ones are rejected, and the remaining structures participate in creating the next generation through cross-over, permutation and mutation. Selection probabilities are derived from the rank of their fitness (i.e. their calculated free energies).

For the cross-over operation we represent atomic positions by real numbers, not by the often used, but unphysical, bit-strings, and new structures are produced- by matching slices (chosen in random directions and with random positions and shifts in the unit cell) of the parent structures. Importantly, we always use the fractional coordinates - unlike Cartesian coordinates, these have good invariance properties under cell shape transformations. Similar cross-over is done for the lattice vectors matrix elements (for this matrix one may always use the upper- triangular form in order to avoid rotations of the cell).

A certain fraction of the new generation is created by permutation (i.e. switching identities of two or more atoms in a structure) and mutation (random change of the cell vectors). The best structure in a given generation always survives, mates and competes in the next generation. 3. For systems with up to ~ 12 atoms in the cell the global minimum is often found within the first few generations, for larger systems this usually takes ~ 20 generations. Among the important results of the simulation are the stable crystal structure and a set of robust meta-stable ones at given pressure-temperature conditions.

Usually, simulations are performed at 0 K (so, the free energy reduces to the enthalpy) - while this makes simulations much more affordable, it is in principle possible to deal with non-zero temperatures and the related entropic effects in the framework of USPEX, since it is possible to calculate the entropies from first principles. We note that the method is efficient even in problematic situations, such as finding correct atomic ordering (thanks to atomic permutations) and polytypism (thanks to spatial cross-over building the structures of children from slices of parent structures).

Detailed tests and applications:

In this section we discuss several applications of USPEX — both tests (i.e. where the stable structure is known) and real applications for finding hitherto unknown structures.

The systems discussed here will include simple "prototypic" cases, some physically, chemically or geophysically interesting compounds, or cases where other methods did not succeed. All the calculations discussed below were done ah initio, in the framework of density functional theory (DFT) within the generalised gradient approximation (GGA).

Elemental carbon was put as one of the cases where crystal structure prediction is "beyond the mortals' ken". The prediction of graphite as the stable structure at 1 atm could be problematic, because the GGA does not fully describe van der Waals forces between graphite layers. Nevertheless, the simulations did indeed produce graphite as the ground state at 1 atm. Simulations at high pressures (100, 300, 500, 1000, 2000 GPa) did not show graphite (it did not appear even among the metastable structures) and yielded other stable structures - diamond below 1000 GPa and "bc8" structure above 1000 GPa. The bc8 structure was known for a long time as a meta-stable structure of Si, and was long suspected to become the stable phase of carbon at pressures above 1000 GPa - our simulations strongly support this idea. Simulations at 1000 GPa and 2000 GPa show very diverse meta-stable structures - among those are diamond and other tetrahedral structures, as well as metallic structures. At 100 GPa we see lonsdaleite ("hexagonal diamond") among meta-stable structures, this structure is experimentally known as a metastable state of carbon at high pressures. Another good meta-stable structure found at 100 GPa and containing 5- and 7-fold rings. All these simulations were performed with 8 atoms in the unit cell, and no information on the cell dimensions or shape was used.

Several simulations were performed for the elements. For Xe at 100-1000 Mbar we predict the hep-structure, in agreement with experimental data. Xenon metallises at 130 GPa - our results confirm that this metallization occurs without a change of structure type. For Fe at pressures of the Earth's inner core (~ 350 GPa) we also predict the hep- structure, in agreement with most experimental and theoretical studies. At high temperatures, however, a transition to the bcc- or fcc-structure may be possible - to investigate this possibility, finite-temperature USPEX simulations are needed. For MgSiO₃ we correctly predict the CaIrO₃-type post-perovskite structure at 120 GPa: according to ab initio simulations and experiments this structure is more stable than perovskite above 100 GPa. USPEX simulations at 80 GPa with 20 atoms/cell and no experimental information correctly produced the perovskite structure stable at those conditions. Interestingly, at 120 GPa simulations most configurations found on the way to the post-perovskite structure contained chains of edge-sharing SiOβ-octahedra (like in the post-perovskite structure), which indicates that such structures become favourable at very high pressures. At 80 GPa, most structures contained corner-sharing octahedra, although edge sharing was also occasionally observed. This demonstrates one of the additional powers of our methodology: it allows one to "scan" the chemistry of the system at given conditions by looking at the stable structure and (importantly) at a large set of meta-stable structures.

For Al₂O₃, experiment and simulations demonstrated that the CaIrO₃-type structure is stable above ~ 130 GPa. Indeed, this pressure corresponds to the conditions where shock-wave experiments evidenced a sudden drop of the electrical resistivity and a small increase of density. Further (and much greater) density anomaly was observed at 250 GPa in literature and was tentatively attributed to a solid-solid phase transition. At least for some materials, e.g., MgSiO₃ and NaMgF₃, the CaIrO₃-type structure seems to be the ultimate high-pressure structure before decomposition, but decomposition can be ruled out for chemically simple Al₂O₃. In search of a post-CalrO₃ type phase, we performed variable-cell evolutionary simulations at 300 GPa, but the result was that the CaIrO₃-type structure is stable at these conditions. This leads us to suggest that the phase transition seen in shock-wave experiments at 250 GPa (and very high temperatures) was melting, rather than any solid-solid transition.

For SiO₂ at 1 atm (simulations with 9 atoms/cell), our variable-cell simulations correctly find quartz as the most stable structure. As meta-stable structures, we find various quartz-like structures, at higher free energies layered tetrahedral structures, then structures containing pairs of edge-sharing SiO₄-tetrahedra, then structure containing isolated trigonal chains of tetrahedra and another structure containing simultaneously SiO₄-tetrahedra and SiOβ-octahedra. All reasonable meta-stable structures (within 1-1.5 eV/cell of the global minimum) are based on corner-sharing SiO₄-tetrahedra, in accordance with Pauling's rules. Structures contradicting Pauling's rules (e.g. those containing edge-sharing tetrahedra or those with octahedral coordination at ambient pressure) come out to be energetically very unfavourable.

Previous attempts to use evolutionary algorithms within the fixed-cell implementation found that it is extremely difficult to reproduce the structure of anatase (12 atoms/cell): the success rate of such calculations was only ~ 10%. With USPEX we were able to easily find this structure.

A good illustration of the power of USPEX is the relatively complex and hitherto unknown structure of the post-aragonite phase of CaCO₃- experimental studies found this phase to be stable above 40 GPa, but could not solve its structure. Post-aragonite is an important mineral phase, likely to be a major host of carbon in the Earth's lower mantle. Finding this structure in not easy, as it contains three sorts of atoms with very different chemical bonds and molecular CO₃ ^2" ions that could rotate with little energy penalty. The structure found by USPEX closely reproduces experimental measurements; in agreement with experiment calculations show that this structure is more stable than aragonite above 42 GPa. The same structure is stable for SrCO₃ at high pressure. With USPEX, one can find not just the most stable structure at given pressure- temperature conditions, but also a set of robust meta-stable phases. The possibility to scan a large number of meta-stable structures and the stable one allows one to extract rich chemical information.

It can be expected that the most challenging materials for this methodology are molecular crystals. There are two reasons: first, the availability of large "empty" space (i.e. the chance of putting an atom by mistake in this "empty" space is very high), second, the possibility of creating different molecules from the same set of atoms - e.g., in the case of urea, the correct (NH₂)₂CO molecules or a set of N₂, H₂, O₂, CO₂, NH₃ molecules etc. Often these different possibilities would be energetically competitive - in such cases one could expect that the overall shape of the free energy landscape would contain several nearly equal attractor basins and the simulation could end up in any of those. However, the case of CaCO₃, close to the case of molecular crystals and discussed above, shows that the method still works. For simple molecular solids the method is also very efficient and reliable (we have done simulations on hydrogen, oxygen, nitrogen, chlorine, fluorine). An interesting test case is ice, H₂O, which at atmospheric pressure adopts the Ih structure (where the oxygen sublattice is isostructural with lonsdaleite polymorph of carbon). Together with graphite, prediction of this structure is, according to literature, "beyond the mortals' ken". One can expect some difficulties for this material - low density and relatively weak bonding between the molecules make the problem quite challenging. More importantly, the structure of ice is frustrated, since it cannot simultaneously fulfil the "ice rules" and keep the full symmetry of the oxygen sublattice. As a result, the structure is disordered and characterised by partial occupancies of hydrogen positions. Disordered structures cannot be treated by USPEX - instead, the method produces the lowest-free energy ordered structure. Fixing the unit cell parameters at experimental values for a 12-atom cell, we correctly obtain the ice Ih structure with correctly fulfilled ice rules. Interestingly, with variable cell parameters for a 12-atom cell we obtain both ice Ih and a slightly less stable cubic ice Ic structures. For an even more challenging case of urea, (NH₂) ₂CO, with 16 atoms in the unit cell, we correctly and easily find the structure when the unit cell parameters were fixed at experimental values. This success is far from trivial: during the simulation we observed a very large number of different molecules and with different orientations in the cell. These structures had similar energies, but the lowest-energy structure was (correctly) corresponding to the tetragonal urea structure. Without any experimental information, the simulation is still successful in identifying the correct molecules of urea, but it takes much longer to find their correct orientations in the unit cell. Indeed, the energies of rotation of the entire molecules are quite small compared to the intramolecular bond energies, thus making the process of finding correct molecular orientations difficult. To overcome this difficulty, one would need to fix the molecular geometry and use only molecular rotations and translations as search parameters.

Discussion:

The reasons why this method is so successful are e.g. rooted in the "good" properties of free energy landscapes mentioned above. Random generation of the first population ensures unbiased sampling of the landscape, and thanks to local optimisation these structures are projected into the chemically interesting part of the landscape (property 1). Local optimisation is very effective for reducing the noise of the landscape. Relatively large basin area of the global minimum further increases the chance of sampling the global minimum. For systems with up to ~ 12 atoms/cell the global minimum is often reached already in the first few generations. Combining the fragments of the selected optimised structures is likely to produce other reasonable structures. By selecting only the best structures (usually, the best ~ 50% of the structures are selected to produce the next generation) introduces an exponential "zooming in" on the most promising area of the configuration space until the global minimum is found (property 4). Since no symmetry constraints are imposed during simulations, symmetry is one of the results of our method. Moreover, this ensures that the resulting structures are mechanically stable and do not contain any unstable F-phonons.

Our approach enables crystal structure prediction without necessitating experimental input. Essentially, the only input is the chemical formula (then, we typically perform simulations for different numbers of formula units in the cell). However, in some cases, one would like to be able also to predict possible stoichiometries. This is especially important for metallic alloys, where it is very difficult to rationalise (and even more to predict a priori) stable stoichiometries. We propose to search for the structure and composition simultaneously. It is possible to incorporate variable composition in our method in order to find stable structures and likely stoichiometries hi a given system (e.g., for MgO-CaO-CO₂ system, it is possible to predict such stoichiometries as MgCO₃, CaCO₃, CaMg(COa)₂, MgO, CaO₅ CO₂).

Finally, the quality of the global minimum found by USPEX depends on the quality of the ab initio description of the system. Present-day DFT simulations (e.g., within the GGA) are adequate for most situations, but it is known that these simulations do not fully describe van der Waals bonding and the electronic structure of Mott insulators. In both areas there are significant achievements, which can be used for calculating ab initio free energies and evaluation of structures hi evolutionary simulations.

Claims

1. A method for determining a crystal structure comprising the steps of initialisation, wherein the population is either created by seeding the first generation with complete structures and if necessary filling up the ranks with random individuals, wherein the structures can be chosen from any external source or can be taken from any other runs of the method or is created randomly, application of hard constraints, wherein all structures must satisfy a set of hard constraints preferentially selected from the group of: interatomic and/or intermolecular distances and/or molecular shape considerations, which in case of interatomic distances preferably should be greater than 0.5 A; cell angles, which preferably should be between 60 and 120 degrees; distances between cell coordinate planes, which preferably should be greater than 2 A; local optimisation, wherein this step is performed for every individual at each generation and wherein preferably first the atoms and/or molecules are relaxed, keeping the lattice fixed and then both lattice and atoms and/or molecules are relaxed simultaneously; selection; production of the next generation with the help of variation operators wherein as variation operators genetic crossover and/or mutation and/or permutation in case of atomic systems, and in case of molecular systems optionally additionally rotational mutation and/or inversion and/or molecular distortion, is used and wherein the probability that a certain structure will participate in a given operation is derived from the rank of its fitness, wherein for this cross-over operation atomic and/or molecular positions are represented by real numbers, and new structures are produced by matching planar slices, wherein further preferably fractional coordinates are used, and wherein similar cross-over is done for the lattice vectors matrix elements, wherein preferably for this matrix the upper-triangular form is used, wherein a certain fraction of the new generation is created by permutation, preferably by switching identities of two or more atoms and/or molecules in a structure, and mutation, preferably by random Gaussian change of the cell vectors, and wherein structures in a given generation can survive and compete/mate in the next generation, preferably only the best one-two.

2. A method according to claim 1, wherein the matching planar slices are chosen in random directions and with random positions and shifts in the unit cell, of the parent structures, wherein preferably, when slices put together do not contain the correct number of atoms and/or molecules, some atoms and/or molecules are added or deleted to maintain the correct chemical composition.

3. A method according to any of the preceding claims, wherein the selection method uses a stochastic selection scheme.

4. A method according to any of the preceding claims, wherein for cross-over two or more individuals are selected and used to produce new individuals, wherein preferably this is achieved by taking spatial and/or functional fractions and combining those.

5. A method according to any of the preceding claims, wherein for mutation the lattice values are changed by applying a strain matrix and/or the atomic positions are mutated by adding random variables, wherein preferably lattice mutation is achieved by multiplying each lattice vector by a matrix [ I + εy ], wherein preferably strains εy are random zero-mean Gaussian variable between -1 and 1.

6. A method according to any of the preceding claims, wherein for permutation at least two atoms and/or molecules of different types are exchanged a variable number of times.

7. A method according to any of the preceding claims, wherein prior to free energy calculations the k-point grids adapt according to the lattice shape.

8. A method according to any of the preceding claims, wherein a certain fraction of every generation is completely discarded and not use for further operations, preferably due to bad fitness values.

9. A method according to any of the preceding claims, wherein as an evaluation function the ab initio free energy of the locally optimized individual is used as fitness value.

10. A method according to any of the preceding claims, wherein the lattice volume of individuals gets rescaled to a constant or a self-evolving value, preferably before local optimisation.

11. A method according to any of the preceding claims, wherein certain or all hard constraints are left out and/or individuals violating certain or all hard constraints are not discarded, wherein in the latter case a fitness penalty for the violating individual can be applied.

12. A method according to any of the preceding claims, wherein local optimisation is applied to only a selection of individuals.

13. A method according to any of the preceding claims, wherein the fitness value is used as basis of selection probability, not the fitness rank.

14. A method according to any of the preceding claims, wherein variables are represented by bits.

15. A method according to any of the preceding claims, wherein the fitness is derived from surrogate functions, evaluating the enthalpy.

16. A method according to any of the preceding claims, wherein parameters of the method are controlled during the run, wherein the parameters are preferably selected from the group population size, permutation/mutation/inversion/rotation/distortion rate, fraction of permuted/mutated/inverted/rotated/distorted/crossed-over individuals per generation, precision of local optimisation.

17. A method according to any of the preceding claims, wherein similarity between individuals is included in the fitness function.

18. A method according to any of the preceding claims, wherein full or partial molecules are used instead of individual atoms.

19. A method according to claim 18, wherein the information on full or partial molecules, carried over from parents to offspring, comprises besides coordinates such as position coordinates, also orientation variables and/or values of internal variables such as bond lengths, torsion angles, deformation parameters.