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

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

5.2 Simulation of Structural Transformations 109

As mentioned above, however, this kind of simulations suffers severely from

the timescale gap problem. Most structural transitions are of ﬁrst order (as a

consequence of the lack of group–subgroup relation between the symmetry

groups) and therefore proceed by nucleation and growth. In this process a free

energy barrier which separates the initial from the ﬁnal structure must be crossed.

If the barrier is large in comparison to the thermal energy k

T, the spontaneous

crossing of the barrier is a rare event and occurs on the time scale which is very

long compared to typical time scales of atomic vibrations. This occurs also in

experiment, in particular in case of the reconstructive phase transitions where

strong covalent bonds are broken across the transition, and some transitions would

not even occur on the experimental time scale unless the kinetics is accelerated by

increasing the temperature. In simulation, of course, the situation is much worse,

since the time scale over which we can observe the evolution of the system is many

orders of magnitude shorter than in experiment.

It is well known that simulation of such activated processes poses a problem if

the barrier height is substantially larger than the thermal energy k

T.Infact,in

simulation of structural transitions we encounter precisely this situation. While

in the real system nucleation is likely to start at structural defects (or at surface),

where the symmetry is already broken locally, in simulation one typically uses a

small system (at most few hundred atoms in the case of ab initio simulations)

arranged as supercell with periodic boundary conditions. Under these conditions

where neither surface nor defects are present it is very unlikely that a spontaneous

nucleation event will occur and structural transformations occur typically in a

collective manner, all atoms undergoing the same displacement at the same time.

The barrier associated with such a collective mechanism could be substantially

higher than the one corresponding to proper nucleation, which tends to make the

timescale gap problem even worse. The most simple solution here is to simply

increase the pressure beyond the thermodynamic transition pressure. This would

decrease the barrier and increase the probability that it will be spontaneously

crossed in the simulation time. When the pressure is high enough and the barrier

becomes close to zero, the system gets close to the point of mechanical instability

and the transition occurs on the fast time scale, typically in few picoseconds.

This time scale can be nowadays easily reached even in ab initio simulations of

systems consisting of several hundreds of atoms. The problem is, however, that the

pressure at which the system loses mechanical stability could be much larger (often

by an order of magnitude) than the equilibrium transition pressure. The transition

observed under the conditions of strong overpressurization is likely to be different

from the one observed at equilibrium, and might even result in the creation of a

different structure. This clearly has a negative impact on the predictive power of

the method, and even if the correct structure is created the mechanism might be

substantially different from the one that would operate at equilibrium conditions.

It is clear that any method aimed at improving this difﬁculty has to address the

timescale gap problem and introduce some mechanism accelerating the barrier

crossing.

110 5 Simulation of Structural Phase Transitions in Crystals: The Metadynamics Approach

We would like to stress that we do not aim here at simulating the structural

transitions in crystals in full complexity, including nucleation processes. Instead we

are looking for a practical technique allowing the simulation of the idealized model

represented by a relatively small supercell with periodic boundary conditions, yet

being able to observe the transitions that include barrier crossing without the need

for an excessive overpressurization.

5.3

The Metadynamics-Based Algorithm

The metadynamics algorithm invented by Laio and Parrinello [9] represents a

generic approach to the simulation of activated processes, such as chemical

reactions, ﬁrst-order phase transitions, or protein folding. It belongs to the family

of methods using biasing potential which is added to the original free energy

with the aim to lower the barrier and thus accelerate its crossing. This class

of techniques includes, e.g., local elevation [10], conformational ﬂooding [11],

or hyperdynamics [12]. The difference among these techniques consists in the

type of the space in which the biasing potential is applied as well as in the

manner in which it is constructed (the details are found in the original papers).

In metadynamics the biasing potential is constructed in the space of a suitably

chosen order parameter with low dimensionality. This important step reduces the

dimensionality of the search space in a substantial manner and makes the problem

more easily manageable. The idea of working in the space of order parameter is

based on the assumption of time separation. One assumes that it is possible to

split the microscopic degrees of freedom of the problem into two groups – the fast

ones which can be equilibrated on the simulation time scale, and the slow ones

which require excessive time. In the ideal case the order parameter – which can

consist of several components – should include all slow degrees of freedom. The

order parameter is then treated by metadynamics while the remaining fast degrees

of freedom are treated by a standard sampling method such as MD or MC.

While it is possible to invent various kinds of order parameters which could drive

the structural transitions between different crystalline phases, it is also possible to

use a simple and generic one which turns out to work well in many cases. The idea

is based on the separation of time scales between acoustic and optical phonons.

While optical phonons essentially represent internal vibrations of atoms inside the

unit cell, acoustic phonons are rather related to the ﬂuctuations of the unit cell of

the crystalline structure. Since the latter degrees of freedom are slower, it is possible

to construct an order parameter based on the unit cell components and use it in

metadynamics. This approach is in spirit close to that of the Parrinello–Rahman

(PR) method, where the variable unit cell is able to distinguish between different

crystalline structures. In the following we show how this idea can be implemented

in a practical scheme. We start with the original version proposed in Ref. [13] and

afterward present the improved version from Ref. [14].

5.3 The Metadynamics-Based Algorithm 111

In the PR method [6] one applies periodic boundary conditions deﬁned by the

simulation box edges



c and uses the 3 ×3matrixh = (



c) as collective

coordinate, or order parameter. The matrix h is treated as a dynamical variable

coupled to the atomistic degrees of freedom under the condition of constant average

pressure. Since the matrix h contains besides 6 degrees of freedom related to the

shape of the cell (3 cell lengths amd 3 angles) also 3 global rotations, it is useful

to exclude these 3 rotational degrees of freedom from the problem. In Ref. [15] it

was proposed to restrict the matrix h to the upper triangular form which effectively

freezes the box rotation.

In such a case we can deﬁne our six-dimensional order parameter as a vector

h = (h

, h

)

. Following the generic metadynamics algorithm [9]

we deﬁne a discrete dynamics by the equations

t+1

+ δh



(5.1)

Here, the driving force φ

=−

∂G

∂

is derived from a modiﬁed Gibbs potential G

which includes a history-dependent term

(

h) = G(

h) +



−

h−



2δh

(5.2)

The history-dependent term [10] in G

(

h) pushes the system out of the local

minimum. The ﬁrst derivative of the free energy with respect to the order parameter

can be expressed as

∂G

∂h

= V



−1

(σ + P)



(5.3)

and requires only an evaluation of the stress tensor

σ by means of an NVT ensemble

average over a relatively short MD (or Monte Carlo) simulation.

The above version of the algorithm has been successfully applied to several

systems [16–20]. In some cases, however, the use of the matrix h directly as

the order parameter may not be convenient since the shape of the free energy

well in the h coordinates may be quite anisotropic. Besides crystal anisotropy

resulting in different stiffness with respect to the different components of h,there

is also a strong coupling among the components of h. In particular, the energy

cost of a deformation where the system is compressed along all axes (volume

compression) is typically much higher than the cost of a deformation where

compression along one axis is compensated by expansion in the perpendicular

directions (volume-conserving deformation). The basin of attraction of a crystal

structure is therefore likely to be narrow in the direction of volume change and long

in the perpendicular directions. Proper exploration of such landscape would require

the Gaussian size to be small with respect to the narrow direction. Consequently,

ﬁlling the well would require a large number of metasteps, which is not practical.

If, on the other hand, the Gaussian size is not small enough, the metadynamics

algorithm does not guarantee that the system escapes from the basin of attraction

112 5 Simulation of Structural Phase Transitions in Crystals: The Metadynamics Approach

of the initial crystal structure via the lowest energy path. In the extreme case of

a large Gaussian (we note that the Gaussian size is equal to the metadynamics

step), the system may escape in a random direction. This is likely to be important

in systems with complex free energy surfaces where there are many competing

pathways leading out of the initial free energy basin. Some of these pathways

may lead to crystalline structures while others may result, e.g., in amorphization;

therefore, it is important to perform the exploration of the free energy surface in a

proper way so that the pathways crossing the low barriers are followed.

The shape of the bottom of the well is described by its curvature. Close to a given

equilibrium crystal structure characterized by a matrix

the Gibbs free energy

can be expanded to second order:

h) ≈ G(

h0) +

(

h −

h0)

h −

h0) (5.4)

The Hessian matrix

= ∂

h)/∂

∂



(5.5)

can be calculated from the ﬁnite differences of the stress tensor at different values of

h, making use of Eq. (5.3). At equilibrium the matrix A has positive real eigenvalues

{λ

} and can be diagonalized by an orthogonal matrix O.

For illustration we show here an example from ab initio calculation for silicon.

Table 5.1 shows the eigenvalues and eigenvectors of the Hessian matrix (5.5)

corresponding to a supercell consisting of 64 atoms in the cubic diamond structure

at a temperature T = 0 K and pressure p = 100 kbar. We note the large spread

of eigenvalues as well as the fact that the eigenvector corresponding to the largest

eigenvalue points along the direction (111 000). Since for an upper triangular matrix

the volume is simply expressed as V = det h =

, the direction of the volume

gradient in the

h space is ∇V = (

,0,0,0), which is approximately

parallel to the direction (111 000) if the diagonal elements of h are not too different.

Table 5.1 Eigenvalues (in units of (kbar

A)) and correspond-

ing eigenvectors of the Hessian matrix (5.5) for a sample of

silicon consisting of 64 atoms in the cubic diamond structure

at temperature T = 0 K and pressure p = 100 kbar.

Eigenvalues

6439.8 6487.9 11858.0 11914.0 12110.7 62638.7

Eigenvectors

0.7915 −0.2048 −0.0050 0.0048 −0.0015 0.5758

−0.5717 −0.5812 −0.0012 0.0058 0.0012 0.5791

−0.2161 0.7875 0.0069 −0.0074 −0.0032 0.5771

0.0012 0.0025 0.6297 0.7506 −0.2003 −0.0021

−0.0027 −0.0001 −0.3348 0.0296 −0.9418 −0.0018

0.0044 −0.0125 0.7009 −0.6601 −0.2699 0.0004

5.4 Practical Aspects 113

In our experience this is typical also for noncubic supercells; the largest eigenvalue

is mainly associated with volume change and is usually larger than the smallest

one by an order of magnitude and more.

In order to treat all degrees of freedom on equal footing it is convenient to

diagonalize the quadratic form (5.4). This can be achieved by expressing the

variables

h in terms of new collective variables s:

−



√

(5.6)

It is easily seen that in the new variables the well becomes spherical: G(s) ≈

h0) +



. The thermodynamic force ∂G/∂s

now reads

∂G

∂s



∂G

∂

√

(5.7)

and can be easily calculated from expression (5.3).

The metadynamics equations in the new variables read

t+1

= s

+ δs



(s) = G(s) +



−

|s−s



2δs

(5.8)

where φ

=−∂G

/∂s. We follow here the prescription of Ref. [14] which relates the

Gaussian height W to the Gaussian width δs,intheformW ∼ δs

. This prescription

is similar to the one proposed in Ref. [11].

We now comment about the calculation of the Hessian matrix. In principle, it

would be optimal to recalculate it and accordingly redeﬁne the collective coordinates

every time the system undergoes a transition to a new structure. However, the

eigenvector corresponding to the largest eigenvalue is in most cases approximately

parallel to the direction (111 000). Therefore a set of coordinates calculated for one

structure might still be usable for simulation of a series of transitions since the

separation between the direction of the volume gradient and the other degrees

of freedom spanning the orthogonal subspace remains approximately preserved.

When the bulk modulus of the system changes considerably which is often the case

when the system changes coordination, the Hessian eigenvalues may also change

substantially. In such a case it is preferred to recalculate the Hessian matrix and

continue metadynamics with new collective coordinates suitable for the new initial

structure.

5.4

Practical Aspects

Here we discuss some practical aspects related to the implementation and practical

use of the technique. The advantage of the algorithm is that it can be very easily

114 5 Simulation of Structural Phase Transitions in Crystals: The Metadynamics Approach

adapted for use with various MD codes. Since the exploration of the Gibbs free

energy landscape is performed by means of a series of NVT simulations (an

alternative continuous version of metadynamics is proposed in Ref. [21]), any MD

code is suitable if it allows for simulation of a crystal in a general (triclinic) supercell

and is able to calculate the stress tensor. The procedure can be performed by a

simple code acting as a driver which writes the input ﬁles for the MD simulation,

then calls the MD code executable via a system call and when the MD run is

completed, extracts the average value of the stress tensor, reads the MD output

ﬁles and writes new input ﬁles for the next run (Figure 5.1). The length of the MD

run should be sufﬁcient to allow both for a good equilibration of the system after

change of the supercell induced by metadynamics and for a good averaging of the

stress tensor. Even more important, however, is to allow sufﬁcient time for the

system to complete the transition, once this starts. Using too short simulation for

one metastep might result in creation of defective or even disordered structures,

since the structural transformation is disturbed by the rapidly changing supercell

before it is complete.

In the case of classical potentials the above requirements do not pose a practical

problem. For systems containing several hundreds of atoms an MD run of several

thousand MD steps, representing several ps, takes a wallclock time of the order

of 1 min. Since a typical metadynamics simulation takes of the order of hundred

metasteps, it can be readily performed in few hours on a PC. The situation is

different in the case of ab initio simulations. If the system is insulating and 

point description is sufﬁcient, for a typical system of about 100 atoms, running on

a parallel computer with 16 cores, one ionic step of ab initio Born–Oppenheimer

MD may take about 10 s. In such a case performing at each metastep an MD

run of several thousand steps might be prohibitively expensive. It is possible to

run metadynamics even performing at each metastep just few hundred MD steps

and in such case one metastep might take just few hours. A typical simulation of

hundred metasteps then could take about 10 days which is still feasible. In case of

systems requiring the use of k-points the situation is clearly more complicated and

here it is an advantage if the ab initio code is parallelized also over k-points and a

sufﬁcient number of computing cores are available.

Since during a structural transition the system might release a large amount of

heat, it is convenient to use some kind of temperature control in order to prevent

melting or creation of structural defects. To this end, simple velocity rescaling

might be sufﬁcient. If available in the MD code, more sophisticated thermostats

such as Berendsen thermostat [22] or Nose-Hoover thermostat [23] can be used.

An important issue is the choice of the Gaussian parameters. Since at the

beginning no information is available about the landscape of the Gibbs free energy

of the system, such as the height of the barrier or the distance between the minima

in the order parameter space, there is no obvious way to determine the appropriate

values. Using the relation W ∼ δs

mentioned in the previous section we need to

specify only one parameter. A convenient way to ﬁnd a starting value is based on

the assumption that the barrier for the transition in the supercell is typically of

the order of several 10 meV multiplied by the number of structural units N

5.5 Examples of Applications 115

start from an initial equilibrated structure h

, r

MD code (NVT simulation)

Metadynamics code

read h

, r

initial

equilibrate the system

average stress tensor σ

write h

, r

final

, σ

read h

, r

final

, σ

calculate h

t+1

get r

t+1

initial

by rescaling r

initial

to h

t+1

write h

t+1

, r

t+1

initial

Figure 5.1 Flow chart of a metadynamics

simulation showing the communication be-

tween the MD code and the metadynamics

driver.

the system. From this assumption one can take the Gaussian height per structural

unit, e.g., of the order of

∼ 10 meV. This represents a resolution in the free

energy space corresponding to a fraction of the total barrier corresponding to the

supercell. For illustration, in a system consisting of N

= 100 structural units,

the choice of a Gaussian width δs = 40 (kbar

)

and corresponding Gaussian

height W = δs

= 1600 kbar

≈ 1 eV corresponds to a resolution of

≈ 10

meV per structural unit and might be a good starting choice. Clearly this choice

is highly dependent on the nature of the interactions in the system as well as

on the simulation parameters, in particular pressure P. While at thermodynamic

transition pressure the barrier might be large, upon increasing the pressure it

becomes lower. The nature of the interactions determines the characteristic energy

scale of structural transformations and clearly a crystal with strong covalent bonds

and reconstructive transformations will have much higher barriers between the

different phases than, e.g., a molecular crystal bound by van der Waals forces. If

the Gaussian size is too large for a given system and conditions, metadynamics will

not perform a proper ﬁlling of the initial free energy well but instead quickly escape

from it in a random direction. Such situation is easily recognized by the fact that the

ﬁrst transition occurs within few metasteps. In such a case it is deﬁnitely necessary

to reduce the Gaussian size in order to obtain a proper exploration of the landscape.

5.5

Examples of Applications

In this section we review several recent applications of the metadynamics technique

to various kinds of crystals. The initial applications using a tight-binding scheme

116 5 Simulation of Structural Phase Transitions in Crystals: The Metadynamics Approach

and classical force ﬁelds were reviewed in Refs. [15, 24]. More recently metadynam-

ics studies using classical potentials were performed for MgSiO

enstatites in Refs.

[25, 26]. These materials are important crustal and upper mantle minerals and the

metadynamics studies succeeded in revealing microscopic details of complex trans-

formation mechanisms of the ortho to high-pressure clinoenstatite transition [25]

(Figure 5.2) and proto to high-pressure clinoenstatite transition [26]. Metadynamics

with a classical force ﬁeld has also been applied to the study of phase transitions

in CdSe [27, 28]. In this system the forward transition from four-coordinated

structures, wurtzite and zinc-blende, to six-coordinated rocksalt structure was stud-

ied. Besides that, the reverse transition induced by decompression of the rocksalt

structure was studied and it was found that the ﬁnal four-coordinated structure

represents a mixture of layers with wurtzite and zinc blende stacking (Figure 5.3).

This interesting outcome is caused by kinetic effects and is in agreement with a

metastep 1 700 800

14001150850

BOC

CC C

BOC

Figure 5.2 Selected relaxed structures from

metadynamics simulation in projections

along [010] (above) and [001] (below). The

SiO

tetrahedra are represented by polyhe-

dra and the Mg cations by spheres. a

∗

the projection of a perpendicular to b and

c. The position of the (100) slip planes and

the respective slip directions are indicated by

arrows. After Ref. [25].

5.5 Examples of Applications 117

(100) layers

[

(100) layers

[

(a)

(b)

(c)

[001]

[010]

[100]

Figure 5.3 Left and right: snapshots of

two different pairs of (001) layers in the

1728-atom CdSe sample during the transition

from six- to four-coordination. (a) During

metastep 353, alternating bonds have bro-

ken in adjacent (100) layers, such that all

the four-membered rings have been con-

verted into six-membered rings. In the pair

of (001) layers shown on the left, the bonds

broken in each (001) layer are the same,

while in those shown on the right alternating

bonds have been broken in each (001) layer.

(b) During metastep 368, each (001) plane

contains an array of hexagonal six-membered

rings. In the two (001) planes on the left

the six-membered rings are stacked directly

above one-another (‘‘eclipsed’’), while in the

planes shown on the right the stacking of

the rings alternates (‘‘staggered’’). (c) The

layers during metastep 368 shown along the

(010) direction. A puckering of each of the

(001) layers occurred. The eclipsed stacking

of the six-membered rings combined with

the puckering led to ABA WZ stacking, while

the staggered stacking of six-membered rings

gave the ABC ZB motif. After Ref. [27].

recent study [29] where both experiment and simulations using the transition path

sampling [30] technique were performed.

In the following we focus on the applications of the technique where meta-

dynamics was combined with density functional theory (DFT) calculations or

neural-network representation of DFT calculations. The ab initio calculations in

thecaseofMgSiO

and silica were performed with the Car–Parrinello method [31]

using the CPMD code [32] and in the case of CO

,C,andCawiththeprojector

augmented-wave (PAW) method [33] as implemented in the molecular-dynamics

program VASP (Vienna ab-initio simulation program) [34, 35].

118 5 Simulation of Structural Phase Transitions in Crystals: The Metadynamics Approach

(a) (b)

Figure 5.4 MgSiO

polytypes (after Ref.

[18]): a – perovskite (space group Pbnm),

d – postperovskite (Cmcm), b,c – newly

found structures 2 × 2(Pnma=Pbnm)and

3 ×1(P2

/m), respectively. Only silicate oc-

tahedra are shown; Mg atoms are omitted

for clarity. In the postperovskite structure,

the predicted plastic slip plane {110} is

shown by an arrow. Arrows also show the

likely slip planes in the other structures.

The ﬁrst application of ab initio metadynamics was the study of the perovskite to

postperovskite transition in MgSiO

[18], following the experimental discovery of

the latter phase [36, 37]. This phase transition is of high importance for geophysics

and seismology since it is related to the seismic anomaly in the D



layer of the

lowermost mantle. The perovskite structure was found to transform to the layered

structure called postperovskite under conditions that correspond to those of the

lowermost mantle (125–136 GPa, 2500–4000 K). The discovery of this transition,

previously unknown, offered a possibility to explain anomalous properties of the



layer. The metadynamics simulations of this transition were performed using

the original version of the algorithm [13] and revealed that the layered structure of

postperovskite and the original perovskite structure are actually end members of

an inﬁnite series of polytypes and the transition mechanism involves plane sliding

(Figure 5.4). Besides that it provided information about the dominant planes of

plastic slip in this system, which is important for the texture and anisotropic

properties of the D



layer.

Metadynamics has also been applied to a number of polymorphs of silica. Silica

is a material of high importance both for geophysics and for practical applications

and it has a large number of polymorphs. Due to the fact that most phase

transformations in this system are strongly reconstructive and require crossing

of high barriers, their investigation poses a problem for both experiment and

simulation. Kinetic effects play a strong role here, causing metastability of many

phases as well as dependence of the structure on the history of the changes of

pressure and temperature. We mention ﬁrst the simulation of transformation of

coesite which was performed using ab initio metadynamics [14, 38]. Experimentally

coesite was found to amorphize upon compression at room temperature [39]

and the same outcome was found in simulation [40]. Ab initio metadynamics was

applied in order to see whether coesite under pressure could possibly transform also

to another crystalline structure if amorphisation could be avoided. The simulation

started from a 48-atom supercell of coesite and resulted in transformation to