Accounts
o f c h e m i c a l r e s e a r c h
March 2011
Volume 44
number 3
pubs.acs.org/accounts
Published on the Web 03/01/2011 www.pubs.acs.org/accounts Vol. 44, No. 3 ’ 2011 ’ 227–237 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 227
10.1021/ar1001318 & 2011 American Chemical Society
How Evolutionary Crystal Structure Prediction
Works;and Why
ARTEM R. OGANOV,*, †, ‡ ANDRIY O. LYAKHOV,† AND
MARIO VALLE§
†Department of Geosciences and Department of Physics and Astronomy, Stony
Brook University, Stony Brook, New York 11794-2100, United States, ‡Geology
Department, Moscow State University, 119992 Moscow, Russia, and §Data
Analysis and Visualization Group, Swiss National Supercomputing Centre
(CSCS), via Cantonale, Galleria 2, 6928 Manno, Switzerland
RECEIVED ON OCTOBER 2, 2010
CONS P EC TU S
O nce the crystal structure of a chemical substance is known, many prop-erties can be predicted reliably and routinely. Therefore if researchers
could predict the crystal structure of a material before it is synthesized, they
could significantly accelerate the discovery of newmaterials. In addition, the
ability to predict crystal structures at arbitrary conditions of pressure and
temperature is invaluable for the study of matter at extreme conditions,
where experiments are difficult.
Crystal structure prediction (CSP), the problem of finding the most
stable arrangement of atoms given only the chemical composition, has long
remained amajor unsolved scientific problem. Two problems are entangled
here: search, the efficient exploration of the multidimensional energy
landscape, and ranking, the correct calculation of relative energies. For
organic crystals, which contain a few molecules in the unit cell, search can
be quite simple as long as a researcher does not need to include many
possible isomers or conformations of the molecules; therefore ranking becomes the main challenge. For inorganic crystals,
quantum mechanical methods often provide correct relative energies, making search the most critical problem. Recent
developments provide useful practical methods for solving the search problem to a considerable extent. One can use simulated
annealing, metadynamics, random sampling, basin hopping, minima hopping, and data mining. Genetic algorithms have been
applied to crystals since 1995, but with limited success, which necessitated the development of a very different evolutionary
algorithm. This Account reviews CSP using one of the major techniques, the hybrid evolutionary algorithm USPEX (Universal
Structure Predictor: Evolutionary Xtallography).
Using recent developments in the theory of energy landscapes, we unravel the reasons evolutionary techniques work for
CSP and point out their limitations. We demonstrate that the energy landscapes of chemical systems have an overall shape and
explore their intrinsic dimensionalities. Because of the inverse relationships between order and energy and between the
dimensionality and diversity of an ensemble of crystal structures, the chances that a random search will find the ground state
decrease exponentially with increasing system size. A well-designed evolutionary algorithm allows for much greater
computational efficiency.
We illustrate the power of evolutionary CSP through applications that examine matter at high pressure, where new,
unexpected phenomena take place. Evolutionary CSP has allowed researchers to make unexpected discoveries such as a
transparent phase of sodium, a partially ionic form of boron, complex superconducting forms of calcium, a novel superhard
allotrope of carbon, polymeric modifications of nitrogen, and a new class of compounds, perhydrides. These methods have also
led to the discovery of novel hydride superconductors including the “impossible” LiHn (n = 2, 6, 8) compounds, and CaLi2. We
discuss extensions of the method to molecular crystals, systems of variable composition, and the targeted optimization of
specific physical properties.
228 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 227–237 ’ 2011 ’ Vol. 44, No. 3
Evolutionary Crystal Structure Prediction Oganov et al.
1. Combinatorial Complexity of the Problem
Following a simple combinatorial argument,1 the number of
possible distinct structures can be evaluated as
C ¼ V=δ
3
N
� �Y
i
N
ni
� �
(1)
where N is the total number of atoms in the unit cell of
volume V, δ is a relevant discretization parameter (for
instance, 1 Å) and ni is the number of atoms of ith type in
the unit cell. Already for small systems (N ≈ 10-20), C is
astronomically large (roughly 10N if one uses δ = 1 Å and
typical atomic volume of 10 Å3).
It is useful to consider the dimensionality of the energy
landscape:
d ¼ 3Nþ3 (2)
where 3N - 3 degrees of freedom are the atomic posi-
tions, and the remaining six dimensions are lattice para-
meters. For example, a systemwith 20 atoms/cell poses a
63-dimensional problem! We can rewrite eq 1 as C ∼
exp(Rd), where R is some system-specific constant. With
such high-dimensional problems, simple exhaustive
search strategies are clearly unfeasible.
The global optimization problem can be greatly simpli-
fied if combined with relaxation (local optimization). During
relaxation, certain correlations between atomic positions set
in: interatomic distances adjust to reasonable values and
unfavorable interactions are avoided to some extent. The
intrinsic dimensionality of this reduced energy landscape
consisting only of local minima (Figure 1) is now
d� ¼ 3Nþ3-K (3)
where κ is the (noninteger) number of correlated dimen-
sions. d* depends both on system size and on chemistry.
We found2 d* = 10.9 (d=39) for Au8Pd4, d*=11.6 (d=99)
for Mg16O16, and d* = 32.5 (d = 39) for Mg4N4H4. The
number of local minima is then
C�∼ exp(βd�) (4)
with β < R, d* < d, and C* , C, implying that efficient
search must include local optimization. Even simple ran-
dom sampling, when combined with relaxation,3 can
deliver correct solutions for systems with N< 8-10.With
USPEX, the limit is much higher, but the exponential
increase of C* with system size means that CSP is an
NP-hard problem and for sufficiently large sizes CSP will
always be intractable. In most cases, we are interested in
systems with N < 20-200, and systems with N < 100 are
tractable, while the range 100 < N < 200 may become
accessible in the foreseeable future.
2. How the Method Works
Evolutionary algorithms work best when the energy (or,
more generally, fitness) landscape has an overall shape, as
in Figure 1. Analysis2 suggests suchoverall shape in the energy
landscapes of chemical systems and implies that evolu-
tionary algorithms are highly appropriate for CSP. Such overall
structure is also expected for landscapes of many physical
properties. In evolutionary simulations, a population of struc-
tures evolves, gradually “zooming in” on the most promising
regionsof the landscape and leading to further reductionof d*.
The evolutionary algorithm USPEX (Universal Structure
Predictor: Evolutionary Xtallography1,4,5), unlike a previous
FIGURE 1. Energy landscape:2 (a) 1D scheme showing the full landscape (solid line) and reduced landscape (dashed line joining local minima); (b) 2D
projection of the reduced landscape of Au8Pd4, showing clustering of low-energy structures in one region.
Vol. 44, No. 3 ’ 2011 ’ 227–237 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 229
Evolutionary Crystal Structure Prediction Oganov et al.
genetic algorithm for crystals6 but similarly to an algorithm
for clusters,7 includes local optimization and treats structural
variables as physical numbers, instead of nonintuitive bin-
ary “0/1” strings (the latter is the defining difference between
“genetic” and more general “evolutionary” algorithms; the
former use binary strings). Other important considerations
are as follows:
(1) The algorithm incorporates “learning from history”
(i.e., offspring structures bear resemblance to the
more successful of the previously sampled structures),
which is done through selection of the low-energy
structures to become parents of the new generation,
survival of the fittest structures, and variation opera-
tors (i.e., recipes for producing child structures from
parents). Acting upon low-energy structures, variation
operators lead, with high probability, to yet other low-
energy structures. Four variationoperators are used in
our method:1,4,5
(i) heredity (creating child structures from planar slabs
cut from two parent structures8)
(ii) lattice mutation (large random deformation ap-
plied to the unit cell shape)
(iii) permutation (swaps of chemical identity in pairs of
chemically different atoms)
(iv) special coordinate mutations (displacements of
the atoms, but not in a fully random way, see
below). For molecular crystals, where the structure
is assembled from entiremolecules (of a particular
isomer), rigid or flexible, the above variation op-
erators act on molecular centers, and additional
variation operators must act on orientation and
conformation of the molecules.
(2) The population should remain diverse, allowing very
different solutions to be produced throughout the
simulation. Diversity can be measured by the collec-
tive quasientropy, Scoll:
Scoll ¼ - (1-Dij)ln(1-Dij)
� �
(5)
where Dij are abstract cosine distances between all
pairs of structures (these distances measure struc-
tural dissimilarity and can only take values be-
tween 0 and 12). Figure 2 shows that in a good
simulation quasientropy retains large values and
can exceed quasientropy of the first random gen-
eration, that is, evolutionary search not only is
more efficient in finding low-energy structures
but also can have more structural diversity than
random search, thus depriving the latter of any
potential advantages.
Initialization of the first generation can be random for
small systems (N < 20). For large systems, most of the
structures produced by random samplingwill be very similar
(Figure 3), disordered and with high energies.2 It will be hard
to produce good structures from such a population. There is
an inverse relationship between the intrinsic dimensionality
and the mean μ of the distance distribution,
μ
- � (d�)-m (6a)
and variance of this distribution,
σ � (d�)-n (6b)
where positivem and n depend on the distance measure
used (cosine vs Euclidean distances).
FIGURE 2. Evolutionary simulation of the binary Lennard-Jones A5B16
(the potential model used here is known to yield low-energy quasi-
crystalline structures9). The insets show the lowest energy as a function
of generation number, and the lowest-energy structure.
FIGURE 3. Distribution of distances between randomly sampled local
minima in a binary Lennard-Jones system AB2.
230 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 227–237 ’ 2011 ’ Vol. 44, No. 3
Evolutionary Crystal Structure Prediction Oganov et al.
To obtain a diverse population, one should reduce the
number of degrees of freedom in the first generation by (i)
assembling initial structures from ready-made building
blocks (molecules, coordination polyhedra, and low-energy
seed structures) or (ii) generating the initial population using
symmetry and/or pseudosymmetry. Since variation opera-
tors break symmetry, structures with different symmetries
will have a chance to emerge. Consider splitting the unit cell
into S subcells. When all ni/S are integers, splitting is done
into identical subcells, introducing additional translational
symmetry.When ni/S is a noninteger, random vacancies are
created to maintain the correct number of atoms (Figure 4),
introducing pseudosymmetry and leading to nontrivial or-
dered structures that are difficult to create otherwise. Such
structures are well-known in nonstoichiometric compounds
and even in the elements, for example, complex high-pressure
phases Cs-III and Rb-III can be represented as supercells of
thebody-centered cubic structurewithadditional atoms (which
can be thought of as additional partially occupied sublattices).
After a generation is completed, locally optimized struc-
tures are compared using their fingerprints,10 and all non-
identical structures are ranked in order of their free energies.
The probability P of selecting a structure to be a parent is
determined by its fitness rank i, e.g. in a linear scheme:
P(i) ¼ P1- (i-1)P1c ,
Xc
i¼1
P(i) ¼ 1 (7)
where c is a selection cutoff. This scheme is superior to
Boltzmann-type selection, because it is insensitive to
peaks and gaps in energy distributions and does not
require an additional parameter (“temperature”) needed
for defining Boltzmann probabilities; a quadratic analo-
gue of (7) often works even better.
Niching (i.e., removal of identical structures using
fingerprints5,11) allows a large number of lowest-energy
structures to be carried over into the next generation,
increasing the learning power, retaining diversity, and en-
abling a more thorough exploration of low-energy meta-
stable structures.
The current algorithm is efficient for systems with <300
degrees of freedom and can be enhanced.5 Directed moves
that have a higher likelihood of leading to lower-energy
structures or new promising areas of the energy landscape
will be essential. For instance, moving the atoms along the
eigenvectors of lowest-frequency phonon modes typically
leads to low-energy structures. For this, one has to construct
and solve the dynamical matrix, which is computationally
extremely demanding at the ab initio level. We have effi-
ciently solved this problem5 by constructing an approximate
dynamical matrix using the bond hardnesses computed
from bond lengths and atomic covalent radii and
electronegativities.
Given the usually good correlation between the energy
and degree of order,2 one could preferentially use the more
ordered pieces of parent structures in the heredity operator.
Figure 5 shows local order in a hypothetical defective
structure of SiO2; clearly, defective regions correspond to
low-order atoms. Giving low-order atoms larger displace-
ments while preserving positions of high-order atoms leads
to a very effective coordinate mutation operator;5 note that
“blind” indiscriminate displacement of all atoms is much
more likely to destroy than to create good structural motifs.4
Figure 6 shows different types of optimization. These
examples show that with very minor adaptations, this
method is powerful for solving a wide range of problems.
The limits of applicability of this approach are not fully
known. Using a partly successful reimplementation of
the method,1,4 Trimarchi and Zunger14 failed to predict
FIGURE 5. Illustration of the concept of local order for defective SiO2.
Low-order atoms are blue; high-order atoms are red. Low-order regions
correspond to the planar defect.
FIGURE 4. Pseudo-subcells for composition A3B6 (atoms A, large black
circles; atoms B, small filled circles; vacancies, empty circles). The true
cell (thick lines) is split into four pseudo-subcells (thin lines).
Vol. 44, No. 3 ’ 2011 ’ 227–237 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 231
Evolutionary Crystal Structure Prediction Oganov et al.
fcc-ordered structures of the Au8Pd4 alloy that would be
competitive with structures predicted by cluster expansion;
however, this failure was probably due to inadequate Bril-
louin zone sampling. With this in mind, we found12 a P21/m
ordering, more stable than any structure found to be stable
by cluster expansion. Thus, atomic ordering can be effi-
ciently predicted. What are the real limitations then? First,
like for any NP-hard problem, dimensionality is a great
restriction. The current version of the method was found to
work very well for dimensionalities below 100-200. Re-
cently, we could determine, at DFT level of theory, a very
complex high-pressure structure of methane with 105
atoms/cell (Q. Zhu, unpublished), but to enable this, we
had to operate with molecular, rather than atomic, entities,
which considerably lowers the dimensionality of the
problem.
Topology of the landscape is an important factor -
single-funnel landscapes (as in Figure 1) are much more
amenable for evolutionary algorithms than multifunnel or,
even worse, featureless landscapes. Energy landscapes
usually have an overall shape with a small number of
funnels,2 but landscapes of some properties may be more
erratic.
Below we consider some this method's recent applica-
tions to high-pressure chemistry; most of these findings
could not be expected by chemical intuition and required
a powerful CSP method, nonempirical, unbiased and cap-
able of arriving at completely unexpected solutions.
3. Applications
3.1. The High-Pressure Chemistry of “Inert” Atomic
Cores and of the Electron Gas. Highly compressible atoms
of alkali and alkali earth metals enter a chemically interest-
ing regime at strong compression when their cores begin to
overlap:15 valence electrons get increasingly “trapped” in
the interstitial space, and valence bandwidth decreases on
compression. Vacant in the free atoms, p- and d-orbitals
become dominant at strong compression, eventually mak-
ing K, Rb, Cs, Ca, Sr, and Ba d-metals, while Li becomes a
predominantly p-element, even adopting the diamond
structure above 483 GPa.16 The most interesting picture
occurs for Na: in this element at megabar pressures, valence
s-, p-, and d-orbitals are populated nearly equally.17
Sodium has a deep minimum on the melting curve at
∼118 GPa and 300 K, belowwhich extremely complex (and
not resolved) crystal structures were observed.18 We
predicted17 that above 250 GPa an insulating hP4 structure
becomes stable, and this prediction was experimentally
FIGURE 6. Predictions of (a) stable crystal structure of MgSiO3 with 80
atoms in the unit cell,12 (b) stable compounds and their structures in a
binary Lennard-Jones system,12 and (c) the hardest structure of TiO2.
13
232 ’ ACCOUNTS OF CHEMICAL RESEARCH ’ 227–237 ’ 2011 ’ Vol. 44, No. 3
Evolutionary Crystal Structure Prediction Oganov et al.
verified,17 although at lower pressures (>200 GPa).19 The
band gap, estimated using the state-of-the-art GW
calculations,17 turned out to be remarkably wide, from ∼2
eV at 200 GPa to over 5 eV at 500 GPa. This implies that the
hP4 phase should be optically transparent already at 200
GPa (this was experimentally confirmed17) and even color-
less above ∼320 GPa. The band gap is due to strong
interstitial localization of electron pairs. The hP4 structure
(Figure 7) can be described in several equivalent ways as
(i) NiAs-type structure where both sites are occupied by
Na atoms.
(ii) Ni-sublattice of the Ni2In structure, In sublattice being
occupied by the interstitial electron pairs. The Ni2In
structure is known to be remarkably dense, consid-
ered to be the highest-pressure phase for AB2 com-
pounds20 until a post-Ni2In transitionwas discovered.
21
(iii) Double hexagonal close packed structure. The stack-
ing of close-packed layers ofNa atoms is CACBCACB...
(underlined layers contain interstitial electron pairs,
Figure 7) and is squeezed by a factor of >2 along the
c-axis, while the interstitial electron pairs form a
nearly ideal hcp ABAB-stacking (c/a ≈ 1.3-1.6).
The hP4 structure of Naminimizes core-valence overlap
and maximizes packing efficiency of the interstitial electron
pairs. It can be called an “electride”,22 that is, a “compound”
made of ionic cores and localized interstitial electron pairs;
electride formshavealsobeenpredicted inLi underpressure,16
and Li was experimentally shown to become a semiconductor
under pressure.23 There are surprisingly faint hints of electride
behavior in K in anarrowpressure range,24while Rb andCs do
not showelectridebehavior at all. This is due to thepresenceof
shallow d-orbitals in heavy alkali and alkali earth elements;
when d-states get populated under pressure, the atom be-
comes more compact
本文档为【ACR-USPEX-2011】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑,
图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
该文档来自用户分享,如有侵权行为请发邮件ishare@vip.sina.com联系网站客服,我们会及时删除。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。