Oganov A.R. (Ed.) Modern Methods of Crystal Structure Prediction

Подождите немного. Документ загружается.

7.1 Theory 149

2) Initialization of the ﬁrst generation, that is, a set of points in the search space

that satisfy the problem constrains.

3) Determination of the quality for each member of the population using the

so-called ﬁtness function.

4) Selection of the ‘‘best’’ members from the current generation as parents from

which the algorithm creates new points (offspring) in the search space by

applying specially designed variation operators to them.

5) Evaluation of the quality for each new member of the population.

6) Selection of the ‘‘best’’ offspring to build the new generation of the population.

7) Repeat steps 4–6 till some halting criteria is achieved.

Onecanseethatsomebiologicaltermsareusedtodescribethisapproach:

•

Population. A set of points in the search space (which we will further call solutions)

that are analyzed as possible candidates for the optimal solution.

•

Parents. A set of solutions that are used to create new candidates for the optimum

solution.

•

Offspring. A set of solutions created from parents using variation operators.

•

Selection. A process that divides solutions into ones that have to ‘‘die’’ and ones

that have to ‘‘survive’’ to build the next generation.

Two main classes of variation operators also have biological names:

•

Heredity operators use a few parent solutions to build one offspring solution.

•

Mutation operators use a single-parent solution to produce a single child.

Why are evolutionary algorithms so effective? No one actually knows, and as we

wrote earlier no one could prove that the best solution found by the algorithm is

in general close to the optimum one. However there are two things that inspired

people to develop this class of algorithms and made us believe that this approach

is indeed a good one. First of all in a correctly designed evolutionary algorithm

the quality of the best solution(s) in each new generation (quantiﬁed by ﬁtness

function) is at least not worse than in the previous one. We can simply keep the best

solution or any given number of best solutions, if none of the offspring solutions

are better than them. And, secondly, you could see the success of the algorithm

executed by Mother Nature if you look into the mirror or hug your pet.

There is one thing to keep in mind during the design and execution of ev-

ery evolutionary algorithm. You have to deal with the trade-off: diversity of the

population versus convergence to the optimum solution. Higher diversity of the

population means that you can explore your search space better. However it slows

down the process of ﬁnding the minimum even if you have a few structures in its

basin of attraction. On the other hand, using elitist approaches and reducing the

diversity helps your population collapse into a local minimum faster while the risk

of omitting the global optimum is higher. The desired balance between diversity

and convergence can usually be achieved only if you can tune it ‘‘on the ﬂy’’ when

you reveal more information about the system and your search space.

150 7 Crystal Structure Prediction Using Evolutionary Approach

7.1.1

Search Space, Population, and Fitness Function

Crystal structure prediction requires us to ﬁnd the global minimum on the

free-energy landscape (that we will call the search space) for system of a given

stoichiometry. Each point on this landscape (solution) represents a crystal structure

with certain atomic positions and lattice vectors. The set of locally optimized

solutionswewillcallapopulation.

One of the features of evolutionary algorithms, which is very helpful for the

crystal structure prediction problem, is their ability to ﬁnd metastable states – good

local minima on the energy landscape that are clearly separated from the global

minimum.

Fitness function describes the quality of each solution and allows us to compare

them. Naturally, the free energy would be the relevant ﬁtness function for a crystal

structure prediction algorithm. Lower free energy will correspond to better solution

and the most stable structure under given conditions will have the lowest ﬁtness

function value.

7.1.2

Representation

As we have already mentioned, the choice of the right representation is crucial

for the effectiveness of the algorithm. One of the reasons why ﬁrst attempts

to design an evolutionary algorithm for crystal structure prediction ([4, 5]; see

also [6]) had small success was counterproductive representation choice. In these

approaches, a discrete grid of atom positions was used. The structural variables of

the crystal were represented by a binary string, and standard evolution operators

for binary strings were applied. These operators were not physically meaningful

and, therefore, the algorithm was basically performing a random walk in the

search space. Obviously it could reliably ﬁnd the global minimum only in the

simplest systems, and even in this case the number of optimizations it has to do is

comparable with the random search methods. Later algorithms use the real-number

representation for atom positions and lattice parameters [7, 8]. This representation

requires more sophisticated variation operators that are better ‘‘optimized’’ for

their task and allow researchers to build powerful methods for crystal structure

prediction.

In the rest of this chapter, we will describe our evolutionary algorithm Universal

Structure Predictor: Evolutionary Xtallography (USPEX) and results that were

achieved with it [8, 9]. USPEX represents the coordinates of atoms in the unit cell

and lattice vectors by real numbers. Therefore, the search space is continuous and

not discrete like in algorithms with binary string representation. The difﬁculty of

the problem is increased, but this is more than compensated by the possibility to

construct physically intuitive and powerful variation operators. Using real number

representation, we also have less risk to omit some shallow local minima.

7.1 Theory 151

Free energy

Order parameter(s)

Figure 7.2 Reducing the noise in the complex search space after local optimization.

7.1.3

Local Optimization and Constrains

Each structure produced by USPEX and similar modern algorithms is locally

optimized. Local optimization is a crucial part of the technique because otherwise

it would be impossible to ﬁnd the lowest free-energy structure for any reasonable

complex system. It also makes the search space much more ‘‘smooth’’ (see

Figure 7.2) and reduces the intrinsic dimensionality of the energy landscape,

see

Table 7.1.

Since complexity of the problem has exponential dependence on the dimension-

ality, local optimization greatly simpliﬁes the task of ﬁnding the global minimum

on the energy landscape. USPEX uses external tools for local optimization. So far

it can use VASP [12] and SIESTA [13] for ﬁrst-principles optimization and GULP

[14] for semiclassical simulations. It is also relatively easy to add support of other

energy-computing engines, if needed.

We would like to mention that not every structure on the energy landscape is a

chemically and physically feasible solution. If we try, for example, to optimize the

structure where two atoms are located in the same spot then the calculation may

1) Intrinsic dimensionality is the minimum

number of variables that are needed to rep-

resent the search space. Local optimization,

by creating correlations between atomic po-

sitions (such as favoring the formation of

some bonds and disfavoring others), de-

creases the intrinsic dimensionality and

dramatically simpliﬁes the crystal structure

prediction problem.

152 7 Crystal Structure Prediction Using Evolutionary Approach

Table 7.1 Intrinsic vs. extrinsic dimensionality

System Extrinsic dimensionality Intrinsic dimensionality

39 10.9

99 11.6

39 32.5

The numbers were produced using the molecular visualization toolkit STM4 [10] using the

Grassberger–Procaccia algorithm with a suitably adapted Camastra sampling correction [11].

‘‘explode’’ the unit cell and we would not get any meaningful results. Therefore, we

have to apply constrains to discard unfeasible solutions. The structure is considered

unfeasible if:

1) The distance between any two atoms is smaller than threshold determined by

user (e.g., there are no known bonds shorter than 0.5

A). One can set different

thresholds for different pairs of atom types separately; for example, the sum of

correspondent atom radii.

2) One of the lattice vectors is too small. User can determine the threshold value;

for example, it can be set to the diameter of the largest atom present in the

system.

3) The angle between two lattice vectors is too small or the angle between the

lattice vector and the diagonal of the parallelogram formed by other two lattice

vectors is too small, see Figure 7.3. One can always choose the lattice vectors

in such a way that the angle between any of them is in the (60

◦

, 120

◦

)range.

To do this, one can replace the longer vector (let’s say a) in a pair violating the

constrain by a



= a − ceil



|a · b|/b



· sign(a · b) ·b see [20] for more details.

Our experience shows that this procedure speeds up local optimization and

improves the efﬁciency of the algorithm as a whole

7.1.4

Initialization of the First Generation

Common way to create the ﬁrst generation in evolutionary approach is uniform

random sampling, aimed at achieving high diversity of the population. In fact, if

you do not have any aprioryinformation about how the possible optimal solution

could look like, random sampling is the most reasonable way to do the unbiased

2) The 60–120 degree conditions are a sim-

pliﬁed version of the cell constraints – more

complete conditions are that a unit cell vec-

tor (and also cell diagonals) should not have

projections onto other cell vectors that are

greater by absolute value than half of the

latter vector [20].

7.1 Theory 153

(a)

(b)

(c)

Figure 7.3 Cell angles constrains [20]. (a) All constrains

satisﬁed, (b) angle between lattice vectors is too small, and

allelogram formed by other two lattice vectors is too small.

initialization. You just pick up random points in the search space and check them

for feasibility. Feasible solutions are added to the ﬁrst generation. This process is

repeated until the desired number of trial solutions is reached. In our case it means

that we randomly choose the lattice vectors and then randomly drop the atoms

into the unit cell. If some information about the optimum solution is known, for

example, unit cell volume or lattice parameters, this can be used as constraints

for optimization. In this case one would, for example, ﬁx the lattice parameters

and vary only atomic positions within the unit cell. User can also ‘‘seed’’ the ﬁrst

generation with the structures that seem reasonable (e.g., those known for similar

compounds, or for the same compound at different conditions, or coming from

previous structure prediction runs) and ﬁll the rest of it with random ones.

High diversity of the ﬁrst generation is the key to the success of the algorithm. If

we do not have a structure in the basin of attraction (so-called funnel, see Figure 7.4)

of the global minimum, then the probability to ﬁnd that optimal solution can be

low.

Random initialization, however, poses a problem relevant for large systems.

When the number of atoms in the unit cell rises, randomly generated structures

become more and more similar [15] from chemical point of view to each other

154 7 Crystal Structure Prediction Using Evolutionary Approach

Funnel of minimum 2

Funnel of minimum 1

Figure 7.4 Basin of attraction (funnel) of a good local mini-

mum is a set of surrounding local minima. Different funnels

typically contain very distinct types of structures and are

separated by high-energy barriers.

and to a disordered system. This can be visualized as trying to build a crystal

by replicating small volumes of liquid. Thus purely random sampling cannot be

used for effective crystal structure prediction in this case. However, there is a good

method that combines relatively high diversity of the ﬁrst population with a high

degree of randomness. We call it unit cell splitting: large cell is split into smaller

subcells that are ﬁlled with atoms randomly and then replicated to ﬁll the full

cell, see Figure 7.5a. Such structures have a higher translational symmetry and are

more ordered, usually leading to a more diverse population. To help to break this

Figure 7.5 Schematic illustration of the (a) subcell and (b)

pseudo-subcell splitting for compositions (a) A

and

(b) A

(atoms A – large black circles, B – medium

dark-gray circles, C – small gray circles, empty circles – va-

cancies). Thick lines show the true unit cell, split into four

(pseudo-)subcells.

7.1 Theory 155

undesired ‘‘additional’’ symmetry by variation operators (described below) as well

as further increase the diversity of the ﬁrst population, different structures have to

be split into different number of subcells.

For cells where the number of atoms (e.g., a prime number) is not good

for splitting into a small number of identical subcells, the algorithm creates

random vacancies to keep the correct number of atoms in the whole unit cell,

see Figure 7.5b. In this case, no additional symmetry is induced and nontrivial

solutions can be found. We also believe that this pseudo-subcell method is able to

improve conventional random sampling methods [16, 17] when dealing with large

systems.

7.1.5

Variation Operators

In general, the choice of variation operators follows naturally from the representa-

tion. Mutation operators usually randomly distort the numbers from the set that

represents the solution, while heredity operators combine different parts of these

sets from different parent solutions into one child solution. For real-number crystal

structure representation, we use two different types of mutation operators – lattice

mutation and atom permutation.

Lattice mutation applies strain matrix with zero-mean Gaussian random strains

to the lattice vectors:



= (I + ε)L





1 + ε

/2 ε

/21+ε

/2 ε

/21+ε





Lattice mutation is shown in Figure 7.6. The position of atoms (their fractional coor-

dinates within the lattice) remains unchanged. This operator allows the algorithm

parent structure

mutant structure,

optimized

Figure 7.6 Lattice mutation applies strain to the lattice vectors. From Oganov et al.[9].

156 7 Crystal Structure Prediction Using Evolutionary Approach

parent structure permutated structure

Figure 7.7 Permutation operator swaps identities of a few

atom pairs. From Oganov et al.[9].

to investigate the neighborhood of good individuals. Also sometimes structures

of similar quality differ essentially in lattice and that is why the premature con-

vergence of the lattice to some ‘‘optimal’’ value reduces the effectiveness of the

algorithm. Lattice mutation operator increases the diversity of the lattices in the

population.

Atom permutation operator swaps chemical identities of atoms in randomly

selected pairs, see Figure 7.7, while lattice remains unchanged. Such swaps

provide the algorithm with steps in the search space that are very far and difﬁcult in

Euclidean distance. This operator is especially useful for systems, where chemically

similar atoms are present.

Heredity operators are vital part of any evolutionary algorithm. If only mutation

operators are present, then the algorithm is identical to a sophisticated version of a

random search (such as D. Wales’s basin hopping method [18]). Heredity operators

are responsible for utilizing and reﬁning the information about the system that we

gather during the execution of the algorithm. Since properties of the crystal are

determined by spatial arrangement of atoms in the unit cell, the most physically

meaningful way to build a heredity operator is to conserve the information from

parents by using spatially coherent pieces (spatial heredity).

To create a child from two parents, the algorithm ﬁrst randomly chooses the

lattice vector and a point on that vector. Then the unit cells of the parent structures

are cut by the plane parallel to other vectors that goes through this point. Planar

slices are matched, see Figure 7.8, and the number of atoms of each kind is

adjusted. For big cells, it is possible to use more than two structures as parents

and combine slices from all of them into a single child structure. The lattice of the

child structure is a weighted averaged of parent lattices.

Altogether three variation operators described above explore the search space

while preserving and reﬁning the good spatial features through generations. It can

be visualized by comparing the best structures from different generations, see, for

example, Figures 7.9 and 7.10.

7.1 Theory 157

slice from parent #1 slice from parent #2

unoptimized offspring optimized offspring

Figure 7.8 Heredity operator combines spatial slices from

different parent structures to form an offspring structure.

From Oganov et al.[9].

7.1.6

Survival of the Fittest and Selection of Parents

When the offspring structures are produced, we have to create a new generation.

Usually it is a good idea to select the best few structures from the previous

generation into the new one. This will make sure that we would not lose the best

structures found during the algorithm execution and is also helpful if we want

to search for low-energy metastable states, in addition to the global minimum.

New generation is ﬁlled with the offspring structures, which are obtained from

parent structures using variation operators. But which structures are to be chosen

as parents? We rank the solutions by their free energy, worst 40% of structures

are discarded (this threshold can be changed by user), and the probability for

structure among best 60% to be chosen as a parent is proportional to its ﬁtness

rank. Stochastic selection is usually slower than deterministic one, when we choose

only the best offspring to pass selection. However, it maintains a relatively high

diversity of the population, which saves the simulation from being ‘‘trapped’’ in

the basin of attraction of a good local minimum that is not the global one.

158 7 Crystal Structure Prediction Using Evolutionary Approach

generation 1, 0.22 eV/atom generation 2, 0.18 eV/atom

generation 3, 0.11 eV/atom generation 11, 0 eV/atom

Figure 7.9 Best structures in generation

from the evolutionary search for boron struc-

ture at 1 atm. The energy per atom shows

the deviation from the ground state. One

can see that parts of the icosahedra that

build the most stable structure started to

appear in the very ﬁrst generations and were

slowly reﬁned during the run. From Oganov

et al.[19].

7.1.7

Halting Criteria

After the new generation is produced and structures are relaxed, algorithm applies

selection and variation operators to create the next one. This process is repeated

until some halting criteria are achieved. Structure with the best enthalpy is

considered as a candidate for the ground state of the system under given conditions.

One can repeat the simulation a few times to increase the conﬁdence of the

results.

Obvious halting criterion would be reaching a certain number of subsequent

generations where the best structure is not changing. However, sometimes we

may need more sophisticated criteria. For example, we can continue calculations

till the certain level of diversity is maintained and stop only when the population is

ﬁlled with structures from a around the best structure. Alternatively, we may know

something about the desired structure, for example, symmetry group or some other