An atomistic model of MgSiO3 perovskite and post-perovskite phases

Computational Materials Science 126 (2017) 351–359

Contents lists available at ScienceDirect

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

An atomistic model of MgSiO3 perovskite and post-perovskite phases C. Pinilla a,b,⇑, M. Acuña-Rojas a, N. Seriani c, S. Scandolo c,d a

Departamento de Física y Geociencias, Universidad del Norte, Km 5 via Puerto Colombia, Barranquilla, Colombia Earth Sciences Department, University College London, Gower Street, WC1E 6BT London, United Kingdom c The Abdus Salam International Centre for Theoretical Physics, Strada Costiera 11, I-34151 Trieste, Italy d INFM/Democritos National Simulation Center, Via Beirut 2-4, I-34014 Trieste, Italy b

a r t i c l e

i n f o

Article history: Received 1 August 2016 Received in revised form 22 September 2016 Accepted 23 September 2016

Keywords: MgSiO3 phases Force fields Molecular dynamics Materials at extreme conditions

a b s t r a c t A classical force field for MgSiO3 polymorphs is presented and tested for the perovskite, post-perovskite and enstatite phases. In this force field each ion has a fixed partial charge and a fixed polarizability and the potential energy is minimised with respect to the ions’ dipole moments at each step of molecular dynamics by iterating them to self-consistency. The potential parameters are obtained by fitting to forces, stresses and energies calculated by density functional theory. The phase transition from perovskite to post-perovskite is predicted to take place at a pressure of 60 GPa at zero temperature. The potential reproduces well structural and thermodynamic properties as well as defects formation in a wide range of pressures, and therefore will be useful for the investigation of MgSiO3 based materials in a geophysical context. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction In recent years computational modelling has become an important tool to characterise properties of geo-materials and study their behaviour at geological conditions. The fact that materials can be analysed at the atomic level allows for a better understanding of many of the features observed in samples from nature. Computational methods make it possible to study materials at high pressure and temperature where experimental techniques fall short. In this line, properties such as equation of state, deformation at high pressures, defects and trace element incorporation for several materials have been performed in recent years [1]. However, the impact of atomistic simulations on Earth science is constrained by two factors: the accuracy of the model and the length and time scales of the simulation. In recent years most of the modelling has been based on very accurate quantum methods based on density functional theory (DFT). These methods make it possible to calculate properties of a wide range of materials with different compositions ab initio, i.e. without a specific assumption for the atomic interactions. Unfortunately, DFT is currently limited to systems containing a few hundreds of electrons and to dynamics with timescales of a few tens of picoseconds. In most of the cases this is enough to study relative stabilities [2], equations of state ⇑ Corresponding author at: Earth Sciences Department, University College London, Gower Street, WC1E 6BT London, United Kingdom. E-mail address: [email protected] (C. Pinilla). http://dx.doi.org/10.1016/j.commatsci.2016.09.032 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

[3,4], crystal defects [5–7] and even aggregation and dynamics of melts [8,9]. For large scale problems other methods such as empirical atomistic force fields are often used. They are based on a description of the system in terms of atoms (i.e. without explicit electrons) interacting with each other through a potential described by an analytic expression. The advantage of such methods is that they can allow for simulations containing up to millions of atoms and for simulation times on the order of nanoseconds. However, the accuracy of the method strongly depends on the functional form and on the choice of the parameters used to describe the interatomic interactions. Parameters are usually optimised to reproduce known experimental properties [10] or high level DFT calculations [11,12]. Simple additive potentials are often used; however, they fail to describe non-central forces [13]. In the case of oxides, the high polarisability of O2
cannot be described by a pair potential. Of particular interest for geophysical research is the study of MgSiO3, which is believed to form 60% of the lower mantle volume. For this oxide, modelling of processes such as heat transport, melt dynamics and atomic fractionation are relevant for the understanding of the Earth’s interior [14]. The complexity of this task is underlined by the recent discovery of a previously unknown phase of MgSiO3, which forms only under the extreme conditions found deep in the lower mantle [15,16]. In the case of MgSiO3, Matsui [17] has developed a pair potential that gives good results for the structural properties of

352

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

the perovskite phase. For the investigation of polymorphs of these oxides it is usually necessary to employ more advanced functional forms, like the shell model, where the polarisation of the species is taken into account by a Drude-like model [13]. To overcome the problem of the transferability to different polymorphs and extreme conditions, Jahn and Madden [18] have recently presented a series of potentials for several solid phases of geological materials containing Ca, Mg, Si, Al and O. Electrostatics are described within their model using multipole expansions of the ions’ charge densities and short-range repulsive forces between ions also take account of aspherical distortions of ions. These potentials describe with good accuracy a wide range of materials and phases but there is considerable computational expense associated with describing dipolar and quadrupolar distortions of ions. A simpler approach with only point charges and dipoles has been followed by Tangney and Scandolo [12] for several oxides such as SiO2, MgO, TiO2, and Al2O3 [12,19,20] as well as for water [21] and oxide interfaces [22]. These potentials allow ions to dipole polarise in response to their environments. The dipole moment of each ion interacts electrostatically with all other ions. For a given set of ion positions and dipole moments, the potential energy of interaction between ions is a sum of pairwise terms. The manybody character of the potential is in the dependence of the set of dipole moments on the positions of all ions [23,24]. These potentials can be reasonably transferable between pressures, temperatures, and phases. They are also relatively inexpensive, computationally, allowing for simulations containing many thousands of atoms. Following this trend, in this paper we present a potential for MgSiO3 capable of reproducing most of the structural properties of the perovskite and post-perovskite phases. This potential is completely based on the fitting of ab initio properties of MgSiO3. We have used it to calculate a range of physical properties of both phases, including equilibrium structures, equations of state, phonons, elastic constants, and heat capacities. In general our potential shows good transferability, making it possible to study MgSiO3 at high temperature and pressure conditions. 2. Computational methods

A classical polarizable potential with fixed point charges and atomic polarizabilities, as used in Ref. [12], was employed in this work. The total potential energy function U used in this work can be divided into two contributions: short range and long range electrostatic. The short range contribution U sr is described by a pairwise potential of the Morse stretch form:

U sr ¼

"

h

cij 1
ðrij =r0ij Þ

Dij e

i

#
2e

ðcij =2Þ½1
ðr ij =r 0ij Þ

;

ð1Þ

i>j

where r ij ¼ jri
rj j is the distance between atoms i and j and Dij , cij and r 0ij are the parameters that specify the strength of the interactions and are part of the set of variables found through fitting. This interaction is truncated at an interionic separation of 18 bohr. The total electrostatic interaction includes charge-charge, charge-dipole and dipole-dipole interaction terms:

U es ¼

# rij rTij 1 Tij ¼ I
3 2 : rij r 3ij

" # X p2 1 XX qi qj rij  p i
2qj 3 i þ pi  Tij  pj þ 2 i j–i r ij 2 a rij i i

ð2Þ

where qi and ai are parameters corresponding to the charge and polarisation of atom i respectively. Tij is the dipole propagator, a second rank tensor

ð3Þ

With rTij the transposed of the vector rij . Here we use the approach by Madden et al. [25,26] to describe an additional contribution to dipoles from the deformation of ions by their near neighbours. In this model the ‘‘short range” polarisation is given by

psr i ¼ ai

X q ri i

i–j

r3ij

f ij ðr ij Þ;

ð4Þ

where

f ij ¼ cij e
bij rij

k 4 X ðbij r ij Þ k¼0

k!

ð5Þ

and bij and cij are parameters to be fit. For each atomic configuration the electrostatic potential U es is minimised with respect to the set of dipole moments. Due to the inter-dependence of the electric field and the dipole moments, the complete set of dipole moments is calculated selfconsistently following the equation m
1 pm gi–j ; frj gj–i Þ þ psr i ¼ ai Eðri ; fpi i ;

ð6Þ

where pm i is the dipole moment of ion i at iteration m of the selfconsistent cycle and Eðri Þ is the total electric field at position ri . 2.2. Fitting the force field parameters. We have performed density functional calculations with an LDA exchange-correlation functional [27] and norm-conserving pseudopotentials as implemented in the Quantum-ESPRESSO code [28]. In this work we have used the set of pseudopotentials developed by de Gironcoli et al. [29] and known to give a good description of many structural properties of solid silicates [30,29]. For the fitting of the potential we used the variant used by Tangney and Scandolo [12,31] of the force-matching method [11]. Force field parameters (fgij g) are fitted by minimizing the function Cðfgij gÞ described as:

C ¼ xf DF þ xs DS þ xE DE;

2.1. Force field model

X

"

ð7Þ

where DF; DS and DE are normalised root-mean-squared differences between classical and quantum forces, stresses and energies respectively [12,32] and xf ¼ 1:0; xs ¼ 0:01 and xE ¼ 0:1 are the weights assigned to these differences during the fitting. As described in Ref. [31], our parameterisation procedure is an iterative one. At each iteration, parameters are fit to DFT forces, energies and stress tensors calculated on ‘‘snapshot” configurations from classical molecular dynamics simulations using the potential parameters from the previous iteration. The iterations stop when a given parameter set fits DFT data on snapshots generated with it as closely as it does the DFT data to which it was originally fit. In this case, selfconsistency was achieved with the fourth update of the parameter set. The fit was converged to within the precision of the fit ðDF; DS; DEÞ that is reported below. Parameterisation was performed on 360-atom supercells constructed by a 3  3  2 repetition of the primitive 20-atom orthorhombic unit cell. The final parameters were fit to DFT calculations on 22 snapshots taken from three different long MD simulations each of which was equilibrated, and then run for at least 10 ps. Each MD simulation was started in the perovskite crystal structure. Of these 22 snapshots, 6 were taken from a constant pressure simulation at P ¼ 120 GPa and T ¼ 3700 K (Simulation 1); 6 were taken from a constant volume simulation at a density of q ¼ 4:95 g cm
3 and T ¼ 8000 K, in which the average pressure was P ¼ 140 GPa (Simulation 2); and 10 were taken from a

353

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

constant volume simulation at a density of q ¼ 4:88 g cm
3 and T ¼ 6000 K, in which the average pressure was P ¼ 116 GPa (Simulation 3). Simulations 2 and 3 were highly disordered and diffusive and were included to ensure that a wide range of configurations were sampled. DFT calculations were performed on these snapshots with an energy cutoff for the planewave basis set used to represent the wavefunctions of 130 Ry. Such a large energy cutoff is necessary to adequately converge the stress tensor. The Brillouin zone was only sampled at the C
point. The final fit to the 22 snapshots was ðDF; DS; DEÞ ¼ ð0:11; 0:004; 0:04Þ and the individual fits to the configurations from simulations 1, 2, and 3 were ð0:09; 0:007; 0:13Þ; ð0:14; 0:003; 0:3Þ, and ð0:10; 0:003; 0:2Þ, respectively. To test the parameter set on the post-perovskite phase, we used it to simulate a 120-atom cell of post-perovskite at a constant pressure of P ¼ 120 GPa and T ¼ 4500 K. We then performed DFT calculations on 50 equally-spaced snapshots from a 50 ps MD simulation and we found that the fit of the parameters to the DFT data was ðDF; DS; DEÞ ¼ ð0:10; 0:007; 0:12Þ. This is our first indication of the transferability of our force field between the perovskite and the post-perovskite phases. The parameter set is presented in Table 1.

Table 2 Lattice parameters calculated with the new force field (FF). Data from first principles calculations and experiments are given for comparison. Perovskite and Enstatite data have been estimated at 0 GPa. Data for post-perovskite are taken at 120 GPa. FF

LDA [35]

Exp. [36]

Enstatite (P = 0 GPa) 17.92 18.06 8.77 8.62 5.24 5.10 823.0 793.9

a0 (Å) b0 (Å) c0 (Å) V0 (Å3)

FF a0 (Å) b0 (Å) c0 (Å) V0 (Å3)

4.78 4.90 6.90 161.2

18.26 8.82 5.19 836.1

LDA [37]

Exp. [38,39]

Perovskite (P = 0 GPa) 4.79 4.92 6.89 162.5

FF

4.77 4.93 6.90 162.5

LDA [33]

Exp. [16]

Post-Perovskite (P = 120 GPa) 2.49 2.46 8.05 8.05 6.14 6.11 123.0 121.1

a0 (Å) b0 (Å) c0 (Å) V0 (Å3)

2.46 8.04 6.09 120.9

3. Results and discussion 4

In order to check the reliability and transferability of the new force field, we have computed the structural properties of three crystalline polymorphs of MgSiO3, namely the perovskite, postperovskite, and enstatite structures. Structural optimizations for these three crystals were performed via the steepest-descent method. The resulting lattice parameters are presented in Table 2, together with experimental data [16] and results of ab initio calculations [33,34]. The calculated lattice parameters are in excellent agreement both with experimental data and ab initio calculations. This is even more significative given that no structures from enstatite were used in the fitting of the potential and, even for this structure, the discrepancies are below 2%. Since we are particularly interested in the behaviour of MgSiO3 at high pressure, we have also tested the behaviour of these phases under compression. In Fig. 1 the energy as a function of volume for perovskite and post-perovskite phases is reported. The atomic positions have been fully relaxed in these calculations. As expected, the perovskite phase has the lowest total energy. However, our model seems to describe in both cases a less compressive solid than the one from reference LDA calculations. The relationship between energy and volume can be described by the Birch-Murnaghan equation of state [40], which we fit to our calculations to give the state parameters at zero pressure reported in Table 3. The resulting bulk moduli, reported in Table 3, are in fair agreement with experiments and somewhat higher than our DFT estimates. Although our estimates are closer to experiments it is important to remember that experiments are usually carried out at conditions of pressure and

Energy (eV/formula unit)

3.1. MgSiO3 polymorphs: structure and enthalpy

Perovskite Post-Perovskite

2

0 30

35

40

45

50

V (Å) Fig. 1. Energy as a function of volume for the three MgSiO3 polymorphs. Quantities are per formula unit. Squares: perovskite; Circles: post-perovskite. Dashed lines correspond to ab initio results.

temperature where thermal expansion can be important and so the comparison with these results should be made with care. As will be shown later, using the quasi-harmonic approximation we have estimated the thermal expansion coefficient of perovskite and post-perovskite at 300 K with this force field to be 2:5  10
5 and 2:2  10
5 K
1. In the case of perovskite this value compares well with the ab initio [37] and experimental [41,43] results of 2:1  10
5 and 1:6
2:2  10
5 K
1 respectively. With the total energies at different volumes it is possible to calculate the static

Table 1 Fitted parameters for the polarizable potential (in a.u.). Parameter

Mg-Mg

Mg-Si

Mg-O

Si-Si

Si-O

O-O

Dms

1:6093  10
7 16.9120 13.3412 – – – 1.3862

2:4395  10
4 4.9446 30.6176 – – – –

3:1869  10
3 10.7814 4.6453 2.0659
1.9366 – –

4:7468  10þ8
40.4770
23.0339 – – – 2.7996

2:3503  10
3 10.4644 4.8625 2.0866
1.3921 – –


2:4500  10
1 13.2057 2.8236 1.4720 0.8044 4.1351
1.3953

cms rms Bij Cij

ai qi

354

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

Table 3 Calculated equation of state parameters at zero pressure and temperature according to the Birch-Murnaghan equation. FF

LDA

Exp. [41–43] 162.5 246–272 3.9–4.0

Perovskite V (Å3) B0 (GPa) B00

161.2 266.9 5.50

162.5 240 3.9

V (Å3) B0 (GPa) B00

158 247 4.98

Post-Perovskite 165 220.5 4.07

enthalpy difference between perovskite and post-perovskite as a function of pressure (Fig. 2). A simple calculation of enthalpy differences between perovskite and post-perovskite indicates that the latter phase becomes more stable above pressures of 60 GPa (Fig. 2). Transition pressures calculated by LDA and GGA by Oganov and Ono [15] are 83.7 GPa and 98.7 GPa which compare fairly well with the estimate from this potential. Furthermore, this difference corresponds to a discrepancy between the present force field calculations and density functional theory of around 0.1 eV per formula unit, which is comparable to typical expected errors of DFT.

3.2. Dielectric properties We have calculated the high-frequency dielectric tensor 1 , which is a diagonal tensor, and the Born effective charge tensor Z for perovskite and post-perovskite. These properties have been calculated directly by looking at the polarisation responses to a small electric field and small atomic displacements [44], respectively. The dielectric constant and the effective charges will be used in the next section to calculate non-analytic corrections to the dynamical matrix in the long wavelength limit. This is important for correctly reproducing the frequency difference between longitudinal optical (LO) and transverse optical (TO) phonons at the Brillouin zone center. As discussed in Refs. [32,20], the screened effective charge tenpffiffiffiffiffiffi sor (Z ab = 1 cd ) is a good indicator of the accuracy of long range electrostatic interactions between ions. Individually, the tensors Z and 1 are much less meaningful tests of the quality of atomistic force fields, whose role is to reproduce forces on atoms, rather than electronic properties.

0.2

ΔΗ [eV]

0 -0.2 -0.4 -0.6 -0.8

50

100

150

Pressure [GPa] Fig. 2. Difference in enthalpy per formula unit at 0 K between perovskite and postperovskite as a function of pressure. Dots represent results from the new force field and squares are data from LDA by Oganov et al. [15]. The dashed and dotted line is there to show the transition point.

The screened effective charge tensors are compared to DFT in Tables 4 and 5. Overall, we have found a good agreement between the calculated values with the new force field and those found from DFT calculations and so we suggest that the new force field should describe reasonably well long range interactions. 3.3. Vibrational properties In the Earth’s mantle the materials are exposed to both high pressure and high temperature and, in fact, the post-perovskite phase is only experimentally observed when both conditions are present. For a potential for MgSiO3 to be useful, it is therefore necessary that it describes vibrational modes correctly, as these play a fundamental role in thermodynamics, in thermal transport and in the propagation of seismic waves. To have a further insight into the performance of the new force field we have calculated the phonon frequencies for the perovskite and post-perovskite phases of MgSiO3. These are calculated, using the PHONOPY package [47], by diagonalising the dynamical matrix. The diagonal matrix is constructed from the forces on all the atoms when particular atoms are displaced slightly from equilibrium. Thermodynamic properties for perovskite and postperovskite are evaluated by calculating the free energy, entropy and phonon specific heat under constant volume in the quasiharmonic approximation. For both perovskite and post-perovskite we have calculated the dynamical matrix at a large number of phonon wavevectors (q-points) by using a 8  8  8 supercell of the primitive cells of perovskite (20-atom) and post-perovskite (10-atom). The atomic displacements were of 0.001 Å. The calculated phonon dispersion curves of the perovskite and post-perovskite phases are shown in Figs. 3 and 4, respectively. At the C point the long-wavelength correction to the LO phonon frequency is added a posteriori. We have plotted only frequencies corresponding to q-points commensurate with the supercell. We have excluded q-points at small but finite q because these were not directly calculated with the supercell we used, and their interpolation is far from trivial [49]. For comparison we have included frequencies calculated for some branches and symmetry points from DFT-LDA by Karki et al. [45] and Tsuchiya et al. [46] for perovskite and post-perovskite respectively. Overall the potential shows a fair description of the vibrational properties. It can be seen that frequencies for perovskite at 0 GPa are in good agreement with experiments and DFT-LDA calculations. However, an increment in pressure gives a set of frequencies that are somewhat lower than the ones from DFT-LDA calculations. In the case of post-perovskite the agreement is similar, with FF frequencies generally being lower that the ones calculated by DFT methods. For both structures, the phonon frequencies increase with increasing pressure, and this has consequences for the thermodynamic properties. It is important to remark here that, despite the fitting process being carried out using high pressure configurations, the force field’s description of high frequency optical vibrations is better at low pressures than at high pressures. As discussed in detail in Ref. [20], for a number of other materials, force fields of the form presented in this paper have underestimated the screened effective charges and the frequencies and dispersions of the highest energy optical phonons [31,50,20]. The parameterisation process seeks the parameter set that gives the best overall fit to the DFT potential energy surface. Presumably, the charges (q) and polarisabilities (a) which provide the best overall description of the phonon band structure, are too small and/or too large, respectively, to correctly reproduce high energy optical vibrations, particularly in the long wavelength limit. Tables 4 and 5 and Figs. 3 and 4 indicate that our potential is also afflicted by this overscreening, at high pres-

355

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359 Table 4 Diagonal components of the screened effective charge tensor calculated with the new force field and from DFT [45] (in italic) for perovskite. Mg Z xx =

pffiffiffiffiffiffi

pffiffiffiffiffiffi

? ¼ Z yy = ?

Z zz =

1.14 1.13 1.14

pffiffiffiffiffi

k

Si Perovskite (P = 0 GPa) 2.16 2.14 2.16

1.18 Z xx =

pffiffiffiffiffiffi

pffiffiffiffiffiffi

? ¼ Z yy = ?

Z zz =

1.02 1.10 1.04

pffiffiffiffiffi

k

2.16 Perovskite (P = 100 GPa) 1.98 2.11 1.97

1.15

2.13

O1

O2


1.15
1.21
0.96


0.99
0.94
1.36


0.94


1.51


1.04
1.16
0.90


0.92
0.97
1.21


0.97


1.42

Table 5 Diagonal components of the screened effective charge tensor calculated with the new force field and from DFT [46] (in italic) for post-perovskite.

Z xx =

pffiffiffiffiffiffi

?

Z yy = Z zz =

?

pffiffiffiffiffi

k

Z xx =

pffiffiffiffiffiffi

Z yy = Z zz =

pffiffiffiffiffiffi

?

pffiffiffiffiffiffi

?

pffiffiffiffiffi

k

Mg

Si

O1

O2

1.09 1.06 1.21 1.31 1.13

Post-perovskite (P = 0 GPa) 2.26 2.11 2.02 1.63 2.16


1.21
1.20
1.12
1.04
0.98


0.93
0.77
0.99
0.86
1.32

1.18

2.00


0.90


1.39

0.96 1.02 1.01 1.22 1.22

Post-perovskite (P = 120 GPa) 2.03 2.09 1.92 1.71 1.94


1.08
1.15
1.02
1.02
0.90


0.83
0.82
0.89
0.88
1.16

1.14

2.02


0.91


1.34

sure, where it was fit. Presumably, a better fit to these phonons could be achieved by reducing a or increasing the q’s, however this would be at the expense of a worse description of other phonon modes. From Table 4 we see that, for perovskite, the screened effective charges are in better agreement with DFT at low pressure. This may be fortuitous, but it does help to explain the improved agreement with DFT on the highest optical phonons frequencies that we see in Fig. 3(a). For post-perovskite (Table 5) it is much less clear that screened charges are in closer agreement with DFT at low pressures and the improvement of the phonon frequencies (Fig. 4) is also less dramatic. 3.4. Elastic properties The elastic constants play a fundamental role in explaining the seismic behaviour of the lower mantle, and are therefore an important benchmark for a classical potential for MgSiO3. The elastic constants for the perovskite and post-perovskite phases have been calculated at zero temperature by introducing one component of the strain at a time and calculating the corresponding stress after relaxation of the atomic positions. The absolute value of the strain was 0.005 in all cases. The results for perovskite, reported in Fig. 5, compare well with the values calculated by DFT-LDA [51]. At zero pressure, the best agreement is found for the elastic components C11, C33, C44 and C66, for which the discrepancy is at most 10% with respect to both experiment and DFT calculations [51], and in some cases substantially lower (Table 6). For the shear components C12 and C13 the discrepancies are also less than 40 GPa, but the small values of these components lead to higher relative errors. At high pressure, the components C11, C33, C44 and C66 tend to be slightly underestimated with respect to the calculated values of Ref. [51]. The shear components, on the contrary,

are overestimated by 15–20% at a pressure of 100 GPa. For the post-perovskite phase, the elastic constants at zero temperature have been calculated for a pressure of 100 and 120 GPa (Table 7). Also in this case, two different behaviours are observed. The elastic constants C11, C22, C33, C44 and C55 differ from those calculated in Ref. [34] by less than 5%. On the contrary, C12, C13 and C66 differ by 30–40%. 3.5. Thermodynamic properties We have calculated the Helmholtz free energy within the quasiharmonic approximation (QHA) to study the thermal expansion of perovskite and post-perovskite. In this case the free energy takes the following form:

FðV; TÞ ¼ U 0 ðVÞ þ

X   1X hxj ðq; VÞ þ kB T ln 1
eð
hxj ðq;VÞ=kB TÞ ; 2 q;j q;j ð8Þ

where the first, second and third term correspond to the internal, zero point and thermal energy contributions respectively and the summations are carried out on a number of atoms and q-points. We have used 34 q-points (29 q-points for post-perovskite) in the Brillouin zone to guarantee full convergence of the free energy. Subsequently, we have estimated the free energy as a function of volume for each temperature to produce equation of state isotherms as shown in Fig. 6. For comparison we provide experimental and DFT results for the volume of perovskite and post-perovskite at various pressures and temperatures. In the case of perovskite the diverse experimental data have been taken from works using samples under pressures of 24–30 GPa and temperatures ranging from 300 K to 3000 K [53,41,54]. For post-perovskite a unique point at a pressure above 125 GPa and 300 K [16] is shown as well the

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

1000

800

800

Frequency (cm )

1000

-1

-1

Frequency (cm )

356

600

400

200

400

200

Γ

X

S

Y

Γ

Z

T

R

0

U

1500

1500

1000

1000

-1

Frequency (cm )

500

0

Γ

X

S

Y

Γ

Z

T

R

0

U

Fig. 3. Phonon band structure of the perovskite phase at pressures of 0 GPa (top) and 100 GPa (bottom). The black filled circles correspond to frequencies calculated with the present potential, the red empty circles represent DFT calculations [45] and the blue full squares at C point are experimental data from Ref. [48]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

results of DFT calculations by Tsuchiya et al. Our results are found to be slightly higher than experiments but overall the potential is able to describe, within a small error, the variation of the volume with pressure and temperature. The calculated thermal expansion coefficients at different temperatures and pressures can be found in Table 3. It is important to note that the validity of the quasiharmonic approximation is questionable at high temperatures, where thermal expansion is dominated by the anharmonicity of the potential energy surface. Typically, thermal expansivities are significantly overestimated by the QHA at high temperatures. The DFT results of Tsuchiya et al. also make use of the QHA. In Fig. 7 the heat capacities at constant volume of the two phases as a function of temperature are reported. The difference in heat capacities between the two phases at the same temperature is very small, and at temperatures above 1500 K the heat capacities are very near to the Dulong-Petit value. Still, the hardening of the phonon frequencies with pressure leads to a distinct change, with the free energies shifting to higher values with pressure, and the heat capacities becoming smaller. It should be noted however that the effect on the heat capacities is noticeable only below 1500 K, i.e. at temperatures not thought to be present in the lower mantle. The ambient values for perovskite of CV is 79.0 Jmol
1 K
1 which compares well with the 81.85 Jmol
1 K
1 found from DFT calculations [45]. Similarly the calculated value of 62.6 Jmol
1 K
1 compares favourably with the DFT and measured values of 58.72 Jmol
1 K
1 and 57.2 Jmol
1 K
1 respectively.

Γ

Z

T

Y

Γ

S

R

Z

Γ

Z

T

Y

Γ

S

R

Z

500

Fig. 4. Phonon band structure of the post-perovskite phase at pressures of 0 GPa (top) and 120 GPa (bottom). The black filled circles correspond to frequencies calculated with the present potential, the red empty circles represent DFT calculations [46]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

1000

Elastic constants (GPa)

-1

Frequency (cm )

0

600

C11 C12 C13 C33 C44 C66

800 600 400 200 0

0

20

40

60

80

100

Pressure (GPa) Fig. 5. Elastic constants of the perovskite phase at zero temperature. Solid line: C11; dashed line: C12; dotted line: C13; dot-dashed line: C33; double-dotted-dashed line: C44; dotted-double-dashed line: C66.

3.6. Defect energetics In order to explore further the performance of this force field we have estimated the energetics of vacancies as a function of pressure. In this case we have defined the defect energy as:

Edef ðVÞ ¼ EðN
1; VÞ
EðN; VÞ
Eint ðVÞ

ð9Þ

357

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

Table 6 Elastic constants (in GPa) for the perovskite phase calculated with the present potential (FF) at zero pressure. For comparison, the DFT calculations from Ref. [51] and the experimental results from Ref. [52] are reported as well.

FF DFT Exp.

C11

C12

C13

C33

C44

C66

470 487 482

124 128 144

165 144 147

389 456 485

170 203 204

112 145 147

Table 7 Elastic constants for the post-perovskite phase calculated with the present potential (in GPa). For comparison, the DFT calculations from Ref. [34] are reported in brackets. Pressure (GPa)

C11

C12

C13

C22

100 120

1171 (1175) 1269 (1270)

520 (366) 603 (425)

440 (289) 503 (329)

840 (862) 923 (937)

Pressure (GPa)

C33

C44

C55

C66

100 120

1151 (1157) 1263 (1264)

277 (264) 304 (291)

245 (242) 269 (264)

590 (368) 660 (412)

180

6

3000 K 1000 K 300 K Static

5

2000 K

V (Å)

Cv [meV/ K]

160

4 3 2

140

0 GPa 100 GPa

1

120

0 0

100

50

150

0

500

1000

1500

2000

1500

2000

Temperature [K]

P (GPa) 6

180

3000 K 2000 K 1000 K 300 K

5

Static

V (Å)

Cv [meV/K ]

160

140

4 3 2 1

120

0 0

100

50

150

0

500

P (GPa) Fig. 6. Volume as a function of pressure for the static lattice and at a range of temperatures within the quasi-harmonic approximation, for perovskite (top) and post-perovskite (bottom). Experimental data for perovskite at 300 K, 1000 K and 2000 K are denoted by circles [53], squares [41] and crosses [54] respectively. In the case of post-perovskite, the large square denotes experimental data at 300 K [16] and other symbols correspond to DFT computer simulations results at 300 K (black dots), 1000 K (crosses) and 2000 K (triangles) [46].

where EðN
1; VÞ is the energy of a system containing one vacancy and EðN; VÞ the energy of the perfect crystal containing N lattice sites at a given volume V. Eint ðVÞ is the subtracted correction energy of the defect-defect interaction defined as [5,55]

1000

Temperature [K] Fig. 7. Heat capacity of the perovskite phase as a function of temperature for perovskite (top) at 0 GPa (circles) and 100 GPa (stars) and for post-perovskite (bottom) at 0 GPa (circles) and 120 GPa (stars).

Eint ¼


cq2



:

ð10Þ

The values of c used in this work range from 1.46 at zero pressure to 1.57 at 200 GPa.  corresponds to the static dielectric constant as has been discussed earlier. We have used a constant volume approach to simulate the vacancies. In this case the full supercell containing 360 (320) atoms

358

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359

Table 8 Calculated defect energies (in eV) for vacancies in perovskite and post-perovskite as a function of pressure (in GPa). 0

25

11.52 53.90 102.70

11.35 56.55 102.20

60

100

60

100

120

11.12 57.92 101.89

10.68 58.48 100.48

10.33 59.01 99.46

120

150

10.28 58.10 12.25

9.90 58.35 11.86

Perovskite VMg VSi VO

Post-Perovskite 10.95 57.25 12.94

VMg VSi VO

10.56 57.86 12.49

perovskite, post-perovskite and enstatite phases, from zero pressure up to 120 GPa. Given the relatively simple potential form and the limited number of parameters employed, we are confident that this potential will be useful for the investigations of several properties of these phases of relevance for the understanding of the behaviour of the inner mantle of the Earth, such as thermal conductivity and viscosity.

Schottky Enthalpy (eV)

50

40

Acknowledgments

30

20 0

100

50

150

P (GPa) Fig. 8. Schottky formation enthalpy for perovskite (squares) and post-perovskite (plus) from calculations using this force field. Circles and diamonds are the respective ab initio results taken from a previous work by Karki and Khanduja [6].

for pervoskite (post-perovskite) was optimised with respect to its lattice parameters as a function of pressure. Subsequently, the defect calculation has been carried out keeping the volume constant. The estimated defect energies with this force field can be found in Table 8. The total energies of the perfect and defective supercell are used to calculate Schottky defect enthalpies, i.e. the formation enthalpy needed to remove a whole MgSiO3 unit from the lattice and add it to the surface. In this way, the Schottky defect energy is only the sum of the extraction energies of one Mg, one Si and three O vacancies minus the energy per formula unit of the perfect lattice. Fig. 8 shows the estimated Schottky defect enthalpies for perovskite and post-perovskite as estimated from this force field. We have also included the results from Karki and Khanduja [6] for comparison. Defect enthalpies estimated with the classical force field are in good agreement with ab initio estimates and increase as a function of pressure. In addition, the estimated defect enthalpies and the effect of pressure are very similar for the two crystalline phases.

4. Conclusions A polarizable potential for the simulation of MgSiO3 polymorphs has been developed and tested. The potential has 28 free parameters which have been obtained through a fit to forces, stresses and energies calculated from density functional theory. Although only highly disordered configurations have been employed in the fit, the potential properly describes a variety of structural, energetic and vibrational properties of MgSiO3 in the

CP and NS acknowledge financial support by the European Commission under the EU-FP7-NMP grant 229205 ‘‘ADGLASS”. CP acknowledges funding from the French Agence National de la Recherche (ANR, Project No. 11-JS5-001 ‘‘CrIMin”) and from an internal agenda I+D+I of Uninorte (2015-005). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

[29] [30] [31] [32] [33]

T.S. Duffy, Phil. Trans. R. Soc. 366 (2008) 4273. L. Li, J.P. Brodholt, D. Alfè, Phys. Earth Planet. Int. 177 (2009) 103. B. Li, J. Zhang, Phys. Earth Planet. Int. 151 (2005) 143. R.M. Wentzcovitch, B.B. Karki, M. Cococcioni, S. de Gironcoli, Phys. Rev. Lett. 92 (2004) 018501. J.P. Brodholt, Am. Mineral. 82 (1997) 1049. B.B. Karki, G. Khanduja, Earth Planet. Sci. Lett. 260 (2007) 201. Z. Du, N.L. Allan, J.D. Blundy, J.A. Purton, R.A. Brooker, Geochim. Cosmochim. Acta 72 (2008) 554. L. Stixrude, B. Karki, Science 310 (2005) 297. J.T.K. Wan, T.S. Duffy, S. Scandolo, R. Car, J. Geophys. Res. 112 (2007) B03208. R.L.C. Vink, G.T. Barkema, W.F. van der Weg, N. Mousseau, J. Non-Cryst. Solids 282 (2001) 248. F. Ercolessi, J.B. Adams, Europhys. Lett. 26 (1994) 583. P. Tangney, S. Scandolo, J. Chem. Phys. 117 (2002) 8898. S. Yip, Handbook of Materials Modeling, Springer, 2005. D. Anderson, New Theory of the Earth, Cambridge University Press, 2002. A.R. Oganov, S. Ono, Nature 430 (2004) 445. M. Murakami, K. Hirose, K. Kawamura, N. Sata, Y. Ohishi, Science 304 (2004) 855. M. Matsui, Phys. Chem. Mineral. 16 (1988) 234. S. Jahn, P.A. Madden, Phys. Earth Planet. Int. 162 (2007) 129. P. Tangney, S. Scandolo, J. Chem. Phys. 131 (2009) 124510. J. Sarsam, M.W. Finnis, P. Tangney, J. Chem. Phys. 139 (2013) 204704. C. Pinilla, A. Irani, N. Seriani, S. Scandolo, J. Chem. Phys. 136 (2012) 114511. N. Seriani, C. Pinilla, S. Scandolo, J. Phys. Chem. C 116 (2012) 11062. D. Herzbach, K. Binder, M.H. Müser, J. Chem. Phys. 123 (2005) 124711. Y. liang, C.R. Miranda, S. Scandolo, Phys. Rev. Lett. 99 (2007) 215504. A.J. Rowley, P. Jemmer, M. Wilson, P.A. Madden, J. Chem. Phys. 108 (1998) 10209. P.A. Madden, M. Wilson, Chem. Soc. Rev. 25 (1996) 339. J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G.L. Chiarotti, M. Cococcioni, I. Dabo, J. Phys.: Condens. Matter 21 (2009) 395502. . B.B. Karki, R.M. Wentzcovitch, S. de Gironcoli, S. Baroni, Science 286 (1999) 1705. B.B. Karki, R.M. Wentzcovitch, S. de Gironcoli, S. Baroni, Geophys. Res. Lett. 28 (2001) 2699. P. Tangney, S. Scandolo, J. Chem. Phys. 119 (2003) 9673. X.J. Han, L. Bergvist, P.H. Dederichs, H. Müller-Krumbhaar, J.K. Christie, S. Scandolo, P. Tangney, Phys. Rev. B 81 (2010) 134108. T. Tsuchiya, J. Tsuchiya, K. Omemoto, R.M. Wentzcovitch, Earth Planet. Sci. Lett. 224 (2004) 241.

C. Pinilla et al. / Computational Materials Science 126 (2017) 351–359 [34] [35] [36] [37] [38] [39] [40] [41] [42]

T. Iitaka, K. Hirose, K. Kawamura, M. Murakami, Nature 430 (2004) 442. S. Jahn, Am. Mineral. 93 (2008) 528. S. Sasaki, Y. Takeuchi, K. Fujino, S. Akimoto, Z. Kristallogr. 156 (1982) 279. A.R. Oganov, J. Brodholt, G. Price, Earth Planet. Sci. Lett. 184 (2001) 555. A. Yeganeh-Haeri, Phys. Earth Planet. Int. 87 (1994) 111. N. Ross, R. Hazen, Phys. Chem. Miner. 16 (1989) 415. F. Birch, Phys. Rev. 71 (1947) 809. E. Knittle, R. Jeanloz, G. Smith, Nature 319 (1986) 214. H. Mao, R. Hemely, Y. Fei, J. Shu, L. Chen, A. Jephcoat, J. Geophys. Res. 96 (1991) 8069. [43] Y. Wang, D. Weidner, R. Liebermann, Y. Zhao, Phys. Earth Planet. Inter. 83 (1994) 40. [44] Y. Liang, C. Miranda, S. Scandolo, J. Chem. Phys. 125 (2006) 194524. [45] B.B. Karki, R.M. Wentzcovitch, S. de Gironcoli, S. Baroni, Phys. Rev. B 62 (2000) 14750.

359

[46] T. Tsuchiya, J. Tsuchiya, R.M. Wentzcovitch, J. Geophys. Res. 110 (2005) B02204. [47] A. Togo, F. Oba, I. Tanaka, Phys. Rev. B 78 (2008) 134106. [48] R. Lu, A.M. Hofmeister, Y. Wang, J. Geophys. Res. 99 (1994) 11795. [49] A. Togo, F. Oba, I. Tanaka, Phys. Rev. B 78 (2008) 134106. [50] J.R. Kermode, S. Cereda, P. Tangney, A.D. Vita, J. Chem. Phys. 133 (2010) 094102. [51] B.B. Karki, L. Stixrude, S.J. Clark, M. Warren, G.J. Ackland, J. Crain, Am. Mineral. 82 (1997) 635. [52] A. Yeganeh-Haeri, Phys. Earth Planet. Int. 87 (1994) 111. [53] N. Funamori, T. Yagi, W. Utsumi, T. Kondo, T. Uchida, J. Geophys. Res. 101 (1996) 8257. [54] H. Mao, P. Bell, J. Geophys. Res. 84 (1979) 4533. [55] M. Leslie, M. Gillan, J. Phys. C. 18 (1985) 973.