Density Functional Theory Calculations and Thermodynamic Analysis

Subscriber access provided by University of Winnipeg Library

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Density Functional Theory Calculations and Thermodynamic Analysis of the Forsterite MgSiO(010) Surface 2

4

Ming Geng, and Hannes Jonsson J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.8b09047 • Publication Date (Web): 06 Dec 2018 Downloaded from http://pubs.acs.org on December 6, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Density functional theory calculations and thermodynamic analysis of the forsterite

M g2SiO4(010) surface Ming Geng∗ † ‡ and Hannes Jónsson¶ § , ,

†Key

,

Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China

‡Institutions

of Earth Science, Chinese Academy of Sciences, Beijing 100029, China

¶Faculty §Dept.

of Physical Sciences, University of Iceland, 107 Rekjavík, Iceland

of Energy Conversion and Storage, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark

E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 37

Abstract The stability of possible termination structures for the (010) surface of forsterite,

M g2 SiO4 ,

is studied using a density functional theory (DFT) based thermodynamic

approach. The DFT calculations are used to estimate the surface Gibbs free energy of various surface structures and compare their stability as a function of the chemical environment. Among 9 possible terminations, the SiO-II, M2, O-II terminations are found to be most stable as conditions range from Mg-poor to Mg-rich. This relative stability order remains the same at elevated temperature.

The surface phase

diagram obtained provides ground for further theoretical studies of chemical processes on forsterite surfaces in terrestrial planets.

Introduction Silicate minerals are the building blocks of terrestrial planets. Within a wide variety of silicate minerals, olivine is the predominant mineral in both Earth's upper mantle and interstellar media in space. Consequently, olivine plays a fundamental role in dening the properties and inuencing the physiochemical processes of the interior of terrestrial planets. Knowledge of the physical and chemical propertied of olivine is of great geophysical and astrophysical interest because of its role in many important processes. Olivines are silicate solid solutions with the general formula (M gx F e2−x )SiO4 . Forsterite is at the Mg rich end of the range corresponding to x = 2 and fayalite corresponds to the iron rich end with x = 0. Both on Earth and in space, the most common olivines are richer in magnesium than in iron. Olivine crystals have orthorhombic structure with space group

Pbmn and are characterized by [SiO4 ]4− tetrahedra linked by the divalent metal cations as indicated in Fig.1. The exposed surface structure or termination of forsterite M g2 SiO4 surfaces determines its ability to adsorb volatile chemicals, an important aspect in many natural processes. For example, the chemisorption of water molecules onto forsterite dust grain has been con2


Page 3 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sidered the possible source of water on terrestrial planets. 14 Despite its low porosity and permeability, the formations of the ultramac rock peridotite provides attractive sites for permanent carbon dioxide sequestration because of the high olivine content and thus signicant potential for carbon dioxide mineralization. 5,6 During the subduction of a tectonic plate, carbonated peridotite (Olivine + CO2 ) and serpentinization (Olivine + H2 O) can transport carbon and water into the deep interior of the Earth. This process is an important pathway of the volatile element cycle on Earth and aects many important features of our planet. Therefore, a fundamental understanding of the surface structure of M g2 SiO4 under various reaction conditions is a prerequisite for atomic scale understanding of various processes on Earth. Several computational studies using empirical force elds 7,8 have been conducted of the main low-index surface of forsterite. Recently,

ab initio

calculations 912 with various basis

sets have also been performed. Both the force eld and ab initio calculations have shown that the (010) surface of forsterite is the most stable one among the various low-index surfaces. Several studies have also been conducted of the absorption of volatile atoms and molecules, such as H atom 11 and water molecule. 13 All these studies have been based on the non-polar surfaces. Polar surface terminations may, however, also be important. An important issue is surface termination, i.e. which atoms are exposed at the (010) surface. Previous studies have not been based on the same surface termination and apparently it has been selected rather arbitrarily. In recent experimental work based on high-resolution X-ray reectivity (HRXR) 14 measurements, various surface terminations and dierent levels of hydrations of the forsterite(010) surface were observed depending on the sample preparation such as the method used for polishing the surface. The termination is, therefore, an important consideration when determining the relative stability of forsterite surfaces and in studies of molecular adsorption. In this article, we present simulations of M g2 SiO4 (010) surfaces with all possible terminations and analyze their structure and surface stability as a function of oxygen partial 3



pressure at high and low temperature using density functional theory (DFT) calculations and thermodynamics. Section II describes the computational methodology and the surface structures. In section III the results obtained from the calculations are presented along with some discussion. Conclusions are summarized in Section IV.

METHODOLOGY Surface structure There are two dierent metal ion sites, three oxygen sites and one silicon site in forsterite,

M g2 SiO4 . In the [010] surface orientation, the atomic layers are stacked 2M g(M 1)−O(O1)− SiO(O3) − 2O(O2) − M g(M 2) − M g(M 2) − 2O(O2) − SiO(O3) − O(O1) − 2M g(M 1)−. Crystal planes with [010] orientation can have 9 dierent surface terminations. In the present study, slabs are constructed with these nine structurally dierent (010) surface terminations. The slabs are labeled after the cut positions of the atomic stack as follows: M1 (metal site 1),

SiO (O3 and Si site), O (O1 site), O2 (O2 site) and M2 (metal site 2). Because M1 (metal site 1) is on a symmetry line, it can have only one surface structure. The other four cut sites (SiO, O, O2 and M2) can form two dierent surface structures in opposite directions. We use a "-II" mark to distinguish between dierent slabs in the same position. Each slab has the same termination on both sides. Of these nine slabs, only the M2 terminated one is stoichiometric.

Cleavage Energy When a forsterite crystal is cleaved, two complementary surfaces are created. There are ve complementary pairs of surfaces for the forsterite [010] orientation. The M2-terminated surface is complementary to itself, and the created slab is stoichiometric. In the present study, asymmetric slabs terminated with dierent surfaces on two sides are used. Since the combination of two slabs with complementary terminations is a has twice the thickness of 4


Page 4 of 37

Page 5 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the crystal unit cell, the cleavage energy per area of a unit cell is calculated as follows:

1 i j bulk + Eslab − 2EM Ecleave (i, j) = − [Eslab g2 SiO4 ] 2

(1)

j bulk i where EM g2 SiO4 is the total energy for a unit cell of M g2 SiO4 , and Eslab and Eslab are the

values of the total energy for complementary surface slabs. The cleavage energy density per unit area is equal to

cleave (i, j) =

Ecleave (i, j) A

(2)

where A is the surface area of a unit cell.

Computational method The DFT calculations make use of the Perdew-Burke-Ernzerhof (PBE) approximation to the exchange-correlation functional 15 with a 500 eV kinetic energy cuto in the plane wave basis set. The projector-augmented wave (PAW) method 16 is used to describe the eect of the inner electrons. Monkhorst-Pack k-point mesh schemes were used in the optimization of crystal phases (8×8×8) and surface terminations (8×8×1). To minimize the eect of the periodic boundary conditions, a 30 Å-thick vacuum space is added on top of the surface slab models. The ground-state atomic geometries of the bulk and surface are obtained by minimizing the forces on each atom to below 0.01 eV/Å. The Vienna Ab-initio Simulation Package (VASP) 1720 is used to perform the DFT calculations.

Thermodynamic stability The results of the DFT calculations are used to estimate the thermodynamic stability of the various forsterite surfaces in equilibrium with O2 and the various oxides. By evaluating the Gibbs free energy, the exchange of atoms between the crystal, its surface and the gas phase is taken into account in the analysis of surface stability. The surface Gibbs free energy is 5



Page 6 of 37

a measure of the excess energy of a semi-innite crystal in contact with matter reservoirs with respect to the bulk crystal. The most stable surface is the one with the smallest surface Gibbs free energy. The surface Gibbs free energy is a function of the chemical potentials of the various atomic species. For example, the surface Gibbs free energy of a forsterite slab is:

1 Ωi = (Gislab − NM g µM g − NSi µSi − NO µO ) 2

(3)

where NM g ,NSi , NO denote the numbers of Mg, Si and O atoms in the slab, Gislab is the Gibbs free energy of the symmetric slab with identical surface terminations on both top and bottom sides, and µM g , µSi and µO are the chemical potentials for the Mg, Si and O atomic species in the forsterite crystal. The surface Gibbs free energy per unit area is

ω=

Ω A

(4)

Because the surface of each slab is in equilibrium with forsterite crystal, the chemical potential of M g2 SiO4 is equal to the crystal Gibbs free energy, which results in the following expression:

bulk µM g2 SiO4 = gM g2 SiO4

(5)

where µM g2 SiO4 is the chemical potential of forsterite, which is equal to the sum of the chemical potentials of all atom types in the M g2 SiO4 crystal:

µM g2 SiO4 = 2µM g + µSi + 4µO .

(6)

The equation for the surface Gibbs free energy as a function of variations of the Mg and O chemical potentials can be simplied as

6


Page 7 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 i slab gM g2 SiO4 ) − ΓiSi,M g µM g − ΓiSi,O µO Ωi = (Gislab − NSi 2

(7)

where the parameters ΓiSi,M g and ΓiSi,O are the excesses in the i terminated surface of Mg i i i (NM g )and O (NO ) atoms with respect to the number of Si atoms(NSi ) in the slab:

bulk NM 1 i 1 i g i i − N ) = (NM ΓiSi,M g = (NM g Si g − 2NSi ) bulk 2 2 NSi

(8)

bulk 1 1 i NO i ΓiSi,O = (NOi − NSi ) = (NOi − 4NSi ) bulk 2 2 NSi

(9)

In order for the Mg and Si atoms and O2 molecule not to have a thermodynamical driving force to leave the M g2 SiO4 crystal and precipitating on the surface, their chemical potential must be less than the Gibbs free energy of the corresponding bulk: bulk µM g ≤ gM g

(10)

bulk µSi ≤ gSi

(11)

Similarly, the precipitation of M gO and SiO2 will not be thermodynamically favored if bulk µM g + µO ≤ gM gO

(12)

bulk µSi + 2µO ≤ gSiO 2

(13)

The vibrational contribution to the surface Gibbs free energy is often neglected, 21 but it can make a signicant contribution when surfaces with similar stabilities are compared. We explicitly introduce the vibrational contribution to the Gibbs free energy. The Gibbs free energy g can be approximated as the sum of the energy and F vib (T ):

g(T ) = E + F vib (T ) − T S conf + pV ≈ E + F vib (T ) 7


(14)


Page 8 of 37

Here, the terms F vib (T ) ,T , p, V and S conf correspond to the vibrational free energy, temperature, pressure, volume, and congurational entropy, respectively. In addition, the crystal volume per molecule is approximately 73.174 Å3 , so the largest pV term is only several meV for 100 atm pressure, which is negligible compared with typical error bars in the DFT calculations. T S conf is also approximated to be zero because of its negligible energetic contributions, as has been done in previous studies. Thus, g in Equation 14 is simplied to have contributions from only E and F vib (T ). The F vib (T ) term can be expressed using the phonon density of states (DOS) as follows:

F vib (T ) =

X 1X hωqj + kB T ln[1 − exp(−hωqj /kB T )] 2 qj qj

(15)

where q is the wave vector, j is the band index, and ωqj is the phonon frequency of the phonon mode labeled by a set {q,j}. The phonon calculations are conducted using density functional perturbation theory (DFPT) as implemented in the VASP/Phonopy software. 22 For a slab containing a single unit cell, some of the calculated phonons have a negative eigenvalue with a small magnitude (represented as negative frequency in the phonon dispersion shown in gure Fig. 4). By increasing the cell size, up to 4x4x1, this problem is reduced. Only phonon modes with positive frequency are included in the evaluation of F vib(T ) and the vale obtained appears to be well converged as a function of the slab size. The variation in

F vib(T ) is no more than 0.02 eV which can easily be ignored in comparison to DFT error bars. Some previously reported studies 23,24 have met with similar problems in the phonon calculation and ascribed them to artifacts in the numerics. Since the phonon calculations are computationally intensive, we take the (2 × 2 × 1) slab to give a good compromise for the estimation of the Gibbs free energy of the slabs studied here. By introducing the deviation in the Mg and Si chemical potentials bulk ∆µM g = µM g − gM g

8


(16)

Page 9 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and bulk ∆µSi = µSi − gSi

(17)

2∆µM g + ∆µSi + 4∆µO = ∆gf (M g2 SiO4 )

(18)

eq. (5) can be rewritten as

where ∆gf (M g2 SiO4 ) is the formation Gibbs free energy of M g2 SiO4 from Mg, Si and O2 in their standard states. bulk bulk bulk ∆gf (M g2 SiO4 ) = gM g2 SiO4 − 2gM g − gSi − 2EO2

(19)

The above boundary conditions can be transformed into:

∆µM g ≤ 0

(20)

2∆µM g + 4∆µO ≥ ∆gf (M g2 SiO4 )

(21)

∆µM g + ∆µO ≤ ∆gf (M gO)

(22)

1 ∆µM g + ∆µO ≥ (∆gf (M g2 SiO4 ) − ∆gf (SiO2 )) 2

(23)

where ∆gf (M gO) and ∆gf (SiO2 ) are the Gibbs free energies of formation for M gO and

SiO2 . The oxygen atoms in M g2 SiO4 are in equilibrium with oxygen gas in the atmosphere above the crystal surface:

1 . µO = µgas 2 O2

(24)

The chemical potential of oxygen gas depends on the gas temperature and partial pressure. This dependence can be used to express the Gibbs free energy of the M g2 SiO4 surface in terms of the temperature and oxygen gas partial pressure. The oxygen gas is approximated as an ideal gas, so the dependence of the chemical potential on pressure can be written as 9



gas µgas O2 (T, P ) = µO2 (T ) + kT ln fO2 (T, P )

=

µgas O2 (T )

P + kT ln( 0 ) P

Page 10 of 37

(25)

where k is the Boltzmann constant, µgas O2 (T ) is the standard chemical potential of oxygen gas at temperature T, evaluated from the NIST-JANAF data. 25 The variation of the oxygen atom chemical potential can be written as

1 ∆µO (T, P ) = µO (T, P ) − EO2 2 1 P gas = [∆GO2 (T, P 0 ) + kT ln( 0 )] + δµ0O 2 P

(26)

0 0 where ∆Ggas O2 (T, P ) is the change in oxygen gas Gibbs free energy at standard pressure P

and temperature T with respect to its Gibbs energy at 298.15 K. The calculated values for 0 0 ∆Ggas O2 (T, P ) are presented in Table 3. δµO is a correction to match the DFT calculation of

O2 (which is known to have a relatively large error) to experimental data. This correction was estimated from the calculations of metal oxides and metals similarly to what has been done previously. 21 In our calculation, δµ0O is 0.1997 eV and 0.4208 eV for SiO2 and M gO, respectively. We use an average value of 0.3102 eV for the calculation hereafter. Finally, we write Equation 3 as

Ωi = φi − ΓiSi,M g ∆µM g − ΓiSi,O ∆µO

(27)

1 i 1 bulk i bulk gM g2 SiO4 ) − ΓiSi,M g gM φi = (Gislab − NSi g − ΓSi,O EO2 2 2 1 i i bulk vib ≈ (Eslab + Fslab − NSi gM g2 SiO4 ) 2 1 i bulk − ΓiSi,M g gM g − ΓSi,O EO2 2

(28)

and φi as:

The surface Gibbs free energy for all surfaces described in Section II A was calculated. 10


Page 11 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The region where the M g2 SiO4 surface is stable with respect to the precipitation of Mg, Si or their oxides is determined using the inequalities.

RESULTS AND DISCUSSION Atomic and electronic structure The optimized lattice constant of the crystal obtained from the DFT calculation is close to what has been reported from previous calculation 10,26 and experimental measurements, 27 as summarized in Table 1. After the relaxation of the slab surface, the topmost atoms are displaced from the ideal crystal positions towards the interior of the slab. These displacements are listed in Table 5, and are evaluated as

drel

origin relaxed ∆Zi,j − ∆Zi,j = . ∆Z0

(29)

origin relaxed is the is the distance between atoms i and j in the unrelaxed slab, ∆Zi,j Here, ∆Zi,j

distance between atoms i and j in the relaxed slab, and ∆Z0 is the length of the unit cell of the forsterite crystal. For drel , a positive sign implies a contraction between the surface layers, whereas a negative sign implies an expansion. The slab models are thicker than twice the crystal unit cell (more than 40 layers), and the displacement of atoms from the crystal position is conned only to layers close to the surface. Figure 5 shows that the displacement becomes negligible beyond 10 layers. A unit cell of bulk crystal forsterite has 19 layers in the [010] direction. In M1, O2 and O2-II terminations, the topmost two layers have the same types of atoms. These atoms remain in the same layer after relaxation. In almost all terminations, the displacement between O and Mg atoms is larger than the displacement between O and Si or O layers. For example, in the M1 and M2 terminations, the Mg atoms in the topmost layers are displaced towards the bulk as indicated by the large positive drel . The O atom in the topmost layer of the O and O-II termination is less constrained. An

11



interlayer expansion appears in the rst two layers of these two terminations. Si and O are also split into two layers at the topmost SiO and SiO-II terminations. The Bader charges, calculated using the grid algorithm, 28 for the atoms in the outermost few layers of all nine terminations are listed in Table 6. For reference, the charges of Mg, Si, and O atoms in M g2 SiO4 crystal are also listed in Table 6. The M2 , O-II and SiOII terminations have the smallest charge dierence among these nine terminations. The variations are all not larger than 0.05 e in the rst 10 layers, as shown in the Table. 6. The other six terminations have similar charge variations of approximately 1.50 e in the rst 10 layers.

Cleavage energy The energy required to create complementary surface pairs for forsterite (010) were calculated using equations (1) and (2), and the results are collected in Table 2. The largest cleavage energy was obtained for a pair of O2- and SiO-terminated (010) surfaces. The lowest energy is the creation of an M2-M2 terminated (010) surface pair. This slab is also the only stoichiometric slab. Therefore, one should expect the formation of these surfaces, when an

M g2 SiO4 crystal is cleaved perpendicular to the [010] direction. Thus, the M2 termination is the one previous forsterite surface-related studies have focused on. 11,29,30 Evolution of these surfaces from one termination to another depends on the possibility of ion exchange between crystal and surface, surface and environment, and mobility of atoms on the surface.

Surface stability of various terminations The cleavage energy is the energy required to split a crystal into two parts with complementary terminations. Therefore, the cleavage energy itself does not give directly the stability of the dierent surfaces formed. The Gibbs surface free energy is a measure of the excess energy of a semi-innite crystal in contact with matter reservoirs, and can be used to analyze the stability of the various surface terminations. The DFT calculations have been used to 12


Page 12 of 37

Page 13 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


determine all parameters needed to calculated the surface free energy of a variety of surfaces and formation energy of the M g2 SiO4 , M gO and SiO2 crystals. The excess of oxygen and Mg atoms with respect to Si atoms in the simulated slabs and values are calculated and listed in Table 4. A spontaneous surface formation line of every surface with a specic termination, presented on phase diagrams (shown on the left part of Fig. 6 with a label denoting the termination), can be determined by solving the equation

Ωi (∆µM g , ∆µO ) = 0

(30)

The direction indicated by the arrow on each spontaneous formation line, shows where the Gibbs free energy becomes positive, suggesting that these terminations are stable and may be exposed on M g2 SiO4 particles under the given conditions. The most stable surface for any particular value of the Mg and O chemical potentials is the surface with the smallest, positive surface Gibbs free energy. The boundaries among the regions of stability for dierent terminations are marked by red dashed lines. The boundary between stability regions for surfaces with terminations i and j is determined by the solution of the equation

Ωi (∆µM g , ∆µO ) = Ωj (∆µM g , ∆µO )

(31)

In the present study, only three surface terminations satisfy the two condition at the same time. Thus, we only show two boundaries of these three surfaces. On the left part of Fig.6, each color block represents the region of stability for a dierent surface termination. Pure M g2 SiO4 exists when the conditions of Equations (20) - (23) are all satised. All of these conditions are shown in Fig. 6 as solid lines. The formation energy of M g2 Si4 , magnesium and silicon oxides determine the position of respective precipitation lines. Precipitation of silicon occurs below the Si precipitation line, magnesium metal precipitates on the right from the magnesium precipitation line. M gO crystal will grow on the right and above the M gO precipitation line, and SiO2 will grow on the left from the SiO2 precipitation

13



line. A pure M g2 SiO4 can only be obtained in the narrow striped region between the M gO precipitation line on the left and the SiO2 precipitation line on the right. At the bottom of the diagram, the strip is limited by the Si precipitation line. The calculated stability diagram is presented in Fig. 6. It shows regions of oxygen and magnesium chemical potentials where the surface free energy [Equation(27)] of the various surface terminations is minimal. Every color block is a stability zone of the corresponding surface slab. At each point on the diagram, the M g2 SiO4 surface is in equilibrium with oxygen gas. The equilibrium is characterized by the oxygen chemical potential. Since a single value of the chemical potential can correspond to a wide range in temperature and pressure, the dependency of the oxygen chemical potential on temperature is shown for a number of values for the gas pressure on the right side of Fig.6. These functions were calculated from experimental data, taken from the thermodynamic tables 25 following the approach described earlier by Eq. (26). The design used for constructing the diagrams in Fig. 6 make it possible to determine the conditions for the oxygen environment that correlates with the points on the phase diagrams on the left side of the gures. To illustrate the calculated results for the surface Gibbs free energy of the various

M g2 SiO4 surfaces, graphs are shown for several specic conditions. The surface Gibbs free energy dened in Eq. 27 for an oxygen gas pressure equal to 1 bar is shown in Fig. 7 and Fig. 8 representing ambient conditions and near-melting conditions, 31 respectively. For smaller ∆µM g (corresponding to Mg-poor conditions), the most stable surface is SiO-II terminated. When the environment is Mg-rich, the M2,+ and O-II terminated surfaces become more stable. As the temperature is increased, the ordering of the stability of these surface terminations is unchanged. It is, however, conceivable that other surfaces may become more stable under extreme conditions and this would need to be addressed in further calculations.

14


Page 14 of 37

Page 15 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CONCLUSIONS We have used DFT calculations combined with thermodynamics to analyze nine possible surface terminations of the forsterite (010) surface. The relative stability of the various surface terminations has been estimated as a function of magnesium and oxygen chemical potentials. While most previous studies have focused on the M2 terminated surface, we nd that the SiO-II and O-II surfaces can also be stable under some conditions. In nature, olivine is a solid solution of forsterite (M g2 SiO4 ) and faylaite (F e2 SiO4 ), so Fe can be present as an impurity in forsterite and this can be important for the surface chemistry of olivine. The M2 terminated surface only exposes the metal ions in the M2 site. Previous studies 32 have shown that Fe atoms are preferably placed in the M1 site of the forsterite crystal, so it is not likely that Fe impurity atoms will be exposed on the forsterite surface. Surfaces that expose the M1 site, such as the O-II termination, are more likely to expose Fe atoms on the surface. Since the presence of Fe on the surface will greatly increase the surface reactivity, the emergence of O-II termination as the most stable one can lead to increased reactivity. This point represents an extension of previous studies 33 of the surface chemistry of forsterite. A stability diagram for comparing the surface Gibbs free energy for the various possible surface terminations of the forsterite(010) surface has been generated. In the Mg-poor condition, the SiO-II termination has the smallest surface Gibbs free energy. When the environment is Mg-rich, the most stable termination will be M2 or O-II. Projections of the phase diagram for ambient and high-temperature conditions have also been generated. The stability ordering of the various terminations does not change as the temperature is increased from room temperature to 2170 K. These calculations are based on the harmonic approximation for vibration, which may, however, not be accurate at such high temperature. The calculated results presented here provide a foundation for further theoretical studies of forsterite surfaces, in particular non-stoichiometric surfaces, possible surface reconstructions, and chemical processes that can take place on forsterite (010) surfaces.

15



Page 16 of 37

Acknowledgement This work was supported by National Science Foundation of China (Grants #41503060 and #41590620), Strategic Priority Research Program (B) of Chinese Academy of Sciences (#XDB18000000 and #XDB10020301). MG thanks Dr. Javed Husssain for help with the VASP calculations and useful discussions at an early stage of this research. The computations were performed on resources provided by the Computer Simulation Lab, IGGCAS and the Icelandic High Performance Computing Centre at the University of Iceland.

References (1) Stimp, M.; Walker, A. M.; Drake, M. J.; de Leeuw, N. H.; Deymier, P. An ångströmsized window on the origin of water in the inner solar system: Atomistic simulation of adsorption of water on olivine.

Journal of Crystal Growth

2006, 294, 8395.

(2) Leeuw, N. H. d.; Catlow, C. R. A.; King, H. E.; Putnis, A.; Muralidharan, K.; Deymier, P.; Stimp, M.; Drake, M. J. Where on Earth has our water come from? Chemical Communications

2010, 46, 89238925.

(3) Muralidharan, K.; Deymier, P.; Stimp, M.; de Leeuw, N. H.; Drake, M. J. Origin of water in the inner Solar System: A kinetic Monte Carlo study of water adsorption on forsterite.

Icarus

2008, 198, 400407.

(4) King, H. E.; Stimp, M.; Deymier, P.; Drake, M. J.; Catlow, C. R. A.; Putnis, A.; de Leeuw, N. H. Computer simulations of water interactions with low-coordinated forsterite surface sites: Implications for the origin of water in the inner solar system. Earth and Planetary Science Letters

2010, 300, 1118.

(5) Trias, R.; Ménez, B.; le Campion, P.; Zivanovic, Y.; Lecourt, L.; Lecoeuvre, A.; SchmittKopplin, P.; Uhl, J.; Gislason, S. R.; Alfredsson, H. A. et al. High reactivity of deep

16


Page 17 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


biota under anthropogenic CO2 injection into basalt.

Nature Communications

2017,

, 1063.

8

(6) Aradóttir, E. S. P.; Sonnenthal, E. L.; Björnsson, G.; Jónsson, H. Multidimensional reactive transport modeling of CO2 mineral sequestration in basalts at the Hellisheidi geothermal eld, Iceland.

International Journal of Greenhouse Gas Control

2012,

,

9

2440. (7) Watson, G. W.; Oliver, P. M.; Parker, S. C. Computer simulation of the structure and stability of forsterite surfaces.

Physics and Chemistry of Minerals

1997, 25, 7078.

(8) de Leeuw, N. H.; Parker, S. C.; Catlow, C. R. A.; Price, G. D. Proton-containing defects at forsterite 010 tilt grain boundaries and stepped surfaces.

American Mineralogist

2000, 85, 11431154. (9) de Leeuw, N. H. Density Functional Theory Calculations of Hydrogen-Containing Defects in Forsterite, Periclase, and α-Quartz. The Journal of Physical Chemistry B

2001,

, 97479754.

105

(10) Bruno, M.; Massaro, F. R.; Prencipe, M.; Demichelis, R.; De La Pierre, M.; Nestola, F. Ab Initio Calculations of the Main Crystal Surfaces of Forsterite (M g2 SiO4 ): A Preliminary Study to Understand the Nature of Geochemical Processes at the Olivine Interface.

The Journal of Physical Chemistry C

2014, 118, 24982506.

(11) Garcia-Gil, S.; Teillet-Billy, D.; Rougeau, N.; Sidis, V. H Atom Adsorption on a Silicate Surface: The (010) Surface of Forsterite.

The Journal of Physical Chemistry C

2013,

, 1261212621.

117

(12) Demichelis, R.; Bruno, M.; Massaro, F. R.; Prencipe, M.; De La Pierre, M.; Nestola, F. First-principle modelling of forsterite surface properties: Accuracy of methods and basis sets.

Journal of Computational Chemistry

17

2015, 36, 14391445.



Page 18 of 37

(13) Asaduzzaman, A. M.; Laref, S.; Deymier, P. A.; Runge, K.; Cheng, H.-P.; Muralidharan, K.; Drake, M. J. A rst-principles characterization of water adsorption on forsterite grains.

Philosophical Transactions of the Royal Society A: Mathematical, Physical and

Engineering Sciences

2013, 371 .

(14) Yan, H.; Park, C.; Ahn, G.; Hong, S.; Keane, D. T.; Kenney-Benson, C.; Chow, P.; Xiao, Y.; Shen, G. Termination and hydration of forsteritic olivine (0 1 0) surface. Geochimica et Cosmochimica Acta

2014, 145, 268280.

(15) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple.

Physical Review Letters

1996, 77, 38653868.

(16) Blöchl, P. E. Projector augmented-wave method.

Physical Review B

1994, 50, 17953

17979. (17) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. B

Physical Review

1993, 47, 558561.

(18) Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of the liquid-metalamorphous-semiconductor transition in germanium.

Physical Review B

1994,

,

49

1425114269. (19) Kresse, G.; Furthmüller, J. Eciency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set.

Computational Materials Science

1996, 6, 1550. (20) Kresse, G.; Furthmüller, J. Ecient iterative schemes for ab initio total-energy calculations using a plane-wave basis set.

Physical Review B

1996, 54, 1116911186.

(21) Reuter, K.; Scheer, M. Composition, structure, and stability of RuO2 (110) as a function of oxygen pressure.

Physical Review B

18

2001, 65, 035406.


Page 19 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22) Togo, A.; Tanaka, I. First principles phonon calculations in materials science. Materialia

Scripta

2015, 108, 15.

(23) Zhang, Q.; Schwingenschlögl, U. Emergence of Dirac and quantum spin Hall states in uorinated monolayer As and AsSb.

Phys. Rev. B

2016, 93, 045312.

(24) Sun, X.; Liu, X.; Yin, J.; Yu, J.; Li, Y.; Hang, Y.; Zhou, X.; Yu, M.; Li, J.; Tai, G. et al. Two-Dimensional Boron Crystals: Structural Stability, Tunable Properties, Fabrications and Applications. (25) Chase, M. W.

, 1603300.

Advanced Functional Materials 27

NIST-JANAF Thermochemical Tables, 4th Edition

; American Institute

of Physics, 1998. (26) José, M. S.; Emilio, A.; Julian, D. G.; Alberto, G.; Javier, J.; Pablo, O.; Daniel, S.-P. The SIESTA method for ab initio order- N materials simulation. Condensed Matter

Journal of Physics:

2002, 14, 2745.

(27) Kirfel, A.; Lippmann, T.; Blaha, P.; Schwarz, K.; Cox, D. F.; Rosso, K. M.; Gibbs, G. V. Electron density distribution and bond critical point properties for forsterite, M g2 SiO4 , determined with synchrotron single crystal X-ray diraction data. istry of Minerals

Physics and Chem-

2005, 32, 301313.

(28) Henkelman, G.; Arnaldsson, A.; Jónsson, H. A fast and robust algorithm for Bader decomposition of charge density.

Computational Materials Science

2006, 36, 354360.

(29) Kerisit, S.; Weare, J. H.; Felmy, A. R. Structure and dynamics of forsterite-scCO2 /H2 O interfaces as a function of water content.

Geochimica et Cosmochimica Acta

2012, 84,

137151. (30) Navarro-Ruiz, J.; Ugliengo, P.; Rimola, A.; Sodupe, M. B3LYP Periodic Study of the Physicochemical Properties of the Nonpolar (010) Mg-Pure and Fe-Containing Olivine Surfaces.

The Journal of Physical Chemistry A

19

2014, 118, 58665875.



Page 20 of 37

(31) Davis, B. T. C.; England, J. L. The melting of forsterite up to 50 kilobars. Geophysical Research

Journal of

1964, 69, 11131116.

(32) Chatterjee, S.; Sengupta, S.; Saha-Dasgupta, T.; Chatterjee, K.; Mandal, N. Site preference of Fe atoms in FeMgSiO4 and FeMg(SiO3 )2 studied by density functional calculations.

Physical Review B

2009, 79, 115103.

(33) Navarro-Ruiz, J.; Ugliengo, P.; Sodupe, M.; Rimola, A. Does Fe2+ in olivine-based interstellar grains play any role in the formation of H2? Atomistic insights from DFT periodic simulations.

Chemical Communications

2016, 52, 68736876.

Table 1: Experimental and Calculated Values for Cell Parameters of Bulk Forsterite (in Å) Lattice Constant a b c

This work 4.7654 10.2396 5.9984

Experiment 27 4.756 10.207 5.980

Gaussian Basis 26 4.804 10.280 6.032

Bruno et al. 10 4.7892 10.2539 6.0092

Table 2: Calculated Cleavage Energies Created pairs of surfaces M1 + O M2 + M2 M2-II + O2-II O-II + SiO-II O2 + SiO

Cleavage energy (eV/unit cell) 11.765 3.619 9.841 5.586 11.904

20


Cleavage energy (J/m2 ) 6.594 2.029 5.516 3.130 6.672

Page 21 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Variation in Gibbs free energy for oxygen gas at standard pressure with respect to its value at 0 K. Data are taken from NIST-JANAFtable 25 T(K) 100 200 298.15 300 400 500 600 700

0 ∆Ggas O2 (T, p )(eV ) -0.146 -0.338 -0.539 -0.545 -0.762 -0.987 -1.219 -1.457

0 ∆Ggas O (T, p )(eV ) -0.073 -0.169 -0.270 -0.273 -0.381 -0495 -0.610 -0.729

T(K) 800 900 1000 1100 1200 1300 1400 1500

gas 0 0 ∆Ggas O2 (T, p )(eV ) ∆GO (T, p )(eV ) -1.699 -0.850 -1.945 -0.973 -2.196 -1.098 -2.450 -1.225 -2.708 -1.354 -2.968 -1.493 -3.232 -1.616 -3.498 -1.749

Table 4: Excesses of O and Mg atoms in slabs with respect to Si atoms and the free energy of formation for various surfaces Surface i M1-term. M2-term M2-II-term O-term O-II-term O2-term O2-II-term SiO-term SiO-II-term

NSi 8 8 8 10 8 8 10 10 10

ΓiSi,M g 1 0 1 -1 1 1 -1 -1 -1

ΓiSi,O 0 0 0 0 1 2 0 -2 -1

φi (eV /unitcell) φi (J/m2 ) 2.645 1.483 1.810 1.014 2.163 1.212 9.154 5.131 -3.057 -1.713 -0.839 -0.469 9.878 5.537 12.815 7.186 8.642 4.844

Table 5: Atom displacement of the topmost layers of nine terminated surfaces (drel (%)) M1 term. M2 term. M2-II term. 4.45 (M1-O1) 4.64 (M2-O2) -1.44 (M2-M2) 0.94 (M2-O2) -1.53 (O1-SiO) 0.22 (O2-SiO) -1.35 (SiO split) -1.51 (SiO split) -0.36 (O2 split) 0.03 (O2-SiO) -1.73 (SiO-O2) 0.16 (SiO-O1) 1.20 (O2-M2) 0.74 (O1-M1) -0.54 (SiO split) 0.36 (M1-O1) -1.16 (SiO-O1) -2.09 (M2-M2) -0.16 (O1-SiO) 1.32 (O1-M1) 0.88 (M2-O2) -0.05 (O2-SiO) -0.99 (SiO split) -1.67 (M1 split) 1.51 (M1-O1) -0.41 (SiO-O1) 0.13 (SiO-O2) 0.62 (O2-M2) -0.58 (O1-SiO) 0.35 (O1-M1)

O term. O-II term. O2 term. -8.73 (O1-SiO) -1.75 (O1-M1) 0.70 (O2-M2) 4.15 (SiO-O2) 3.84 (M1-O1) 1.15 (M2-M2) -5.17 (O2 split) -1.42 (O1-SiO) -0.12 (M2-O2) 1.89 (O2-M2) -0.23 (SiO split) -0.13 (O2 split) -0.48 (M2-M2) -0.66 (SiO-O2) -0.25 (O2-SiO) -0.61 (M2-O2) 1.42 (O2-M2) -0.72 (SiO split) -1.60 (O2 split) -2.27 (M2-M2) 0.04 (SiO-O1) -0.65 (O2-SiO) 0.75 (M2-O2) 0.27 (O1-M1) -0.71 (SiO-O1) -0.51 (O2-SiO) 0.34 (M1-O1) 0.65 (O1-M1) -0.40 (SiO split) -0.53 (O1-SiO)

21


O2-II term. SiO term. SiO-II term. -0.12 (O2 split) -3.09 (SiO split) -3.38 (SiO split) -1.99 (O2-SiO) 0.47 (Si-O1) 3.01 (Si-O2) -1.57 (SiO split) -0.32 (O1-M1) -0.78 (O2-M2) -0.11 (SiO-O1) 2.76 (M1-O1) 2.33 (M2-M2) 1.55 (O1-M1) -1.75 (O1-SiO) -1.28 (M2-O) -0.13 (M1 split) -0.22 (SiO split) -0.14 (O2-Si) -1.20 (M1-O1) -0.68 (SiO-O2) -0.18 (SiO split) 0.29 (O1-SiO) 0.70 (O2-M2) 0.33 (O-O1) -0.43 (SiO split) -0.51 (M2-M2) -0.20 (O1-M1) -0.21 (SiO-O2) 0.30 (M2-O2) 0.07 (M1-O1)


Table 6: Bader charges of the bulk crystal and slab with dierent terminations Unit Cell atom charge Mg(M1) 1.67 Mg(M2) 1.69 O(O1) -1.63 O(O2) -1.61 O(SiO) -1.62 Si(SiO) 3.1

M1 term. atom charge Mg(M1) 1.59 O(O1) -1.60 O(SiO) -1.62 Si(SiO) 1.62 O(O2) -1.61 Mg(M2) 1.68 Mg(M2) 1.68 O(O2) -1.61 Si(SiO) 3.08 O(SiO) -1.60

M2 term. atom charge Mg(M2) 1.66 O(O2) -1.60 Si(SiO) 3.09 O(SiO) -1.58 O(O1) -1.62 Mg(M1) 1.66 O(O1) -1.61 O(SiO) -1.61 Si(SiO) 3.08 O(O2) -1.60

M2-II term. atom charge Mg(M2) 0.44 Mg(M2) 1.45 O(O2) -1.63 O(O2) -1.70 O(SiO) -1.66 Si(SiO) 3.09 O(O1) -1.64 Mg(M1) 1.66 O(O1) -1.62 O(SiO) -1.60

O term. atom charge O(O1) -0.72 O(SiO) -0.89 Si(SiO) 3.08 O(O2) -1.56 O(O2) -1.58 Mg(M2) 1.68 O(O2) -1.59 O(O2) -1.62 Si(SiO) 3.10 O(SiO) -1.62

O-II term. atom charge O(O1) -1.60 Mg(M1) 1.64 O(O1) -1.59 O(SiO) -1.62 Si(SiO) 3.08 O(O2) -1.61 Mg(M2) 1.69 Mg(M2) 1.68 O(O2) -1.60 O(SiO) -1.59

O2 term. atom charge O(O2) -0.86 O(O2) -0.88 Mg(M2) 1.67 Mg(M2) 1.68 O(O2) -1.58 Si(SiO) 3.08 O(SiO) -1.60 O(O1) -1.61 Mg(M1) 1.66 O(O1) -1.61

Graphical TOC Entry 2P

M

re c

gO

ip

ita

0

-4

O2-II term. Slab

∆µO(eV)

itat

-2

n

ip

O2 te

ita

tio

rm . Sla

n

b

Oxygen Releasing Mg Precipitation

Si P rec ip

Pr ec

tio

ion

SiO

Si O -II m. Sla te b rm

ter

.S

O

la

-II

b

-6

te

rm .S

la

O-II term. Stable Zone

b

SiO-II term. Stable Zone M2 term. Stable Zone Mg2SiO4 Crystal Stable Zone

-10

-8

-6

-4

∆µMg(eV)

22

-2

0


2

M1 term. Slab

O

M2-II term. Slab

Si

2

O term. Slab

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 37

O2-II term. atom charge O(O2) -1.06 O(O2) -1.07 O(SiO) -1.49 Si(SiO) 3.07 O(O1) -1.35 Mg(M1) 1.67 O(O1) -1.55 O(SiO) -1.59 Si(SiO) 3.08 O(O2) -1.54

SiO term. atom charge Si(SiO) 1.58 O(SiO) -1.59 O(O1) -1.59 Mg(M1) 1.66 O(O1) -1.63 O(SiO) -1.62 Si(SiO) 3.08 O(O2) -1.61 Mg(M2) 1.68 -1.6 O(O2)

SiO-II term. atom charge O(SiO) 1.50 Si(SiO) 2.95 O(O2) -1.56 Mg(M2) 1.68 O(O2) -1.61 Si(SiO) 3.08 O(SiO) -1.61 O(O1) -1.62 Mg(M1) 1.66 O(O1) -1.62

Page 23 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Structure of the forsterite unit cell and M1 termination

23



M2 termination

Page 24 of 37

M2-II termination

Top View Side View Mg O2 termination

O

O2-II termination

Si

Figure 2: M2 and O2 termination surface structure

24


Page 25 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


O termination

O-II termination

Top View Side View SiO termination

Mg

SiO-II termination

O Si

Figure 3: O and SiO termination surface structure

25



1×1×1 Slab

35

Fvib(0K) = 2.343 eV

30

25

25

20

20

Frequency (THz)

30

15 10

0

0

35

3×3×1 Slab

X(0.5,0,0)

Wave Vector

Γ (0,0,0)

-5 Γ (0,0,0)

Y(0,0.5,0)

35

Fvib(0K) = 2.362 eV

30

30

25

25

20

20

15 10

0

0

Wave Vector

Γ (0,0,0)

-5 Γ (0,0,0)

Y(0,0.5,0)

Wave Vector

Γ (0,0,0)

Y(0,0.5,0)

Γ (0,0,0)

Y(0,0.5,0)

Fvib(0K) = 2.362 eV

10 5

X(0.5,0,0)

4×4×1 Slab

X(0.5,0,0)

15

5

-5 Γ (0,0,0)

Fvib(0K) = 2.361 eV

10 5

-5 Γ (0,0,0)

2×2×1 Slab

15

5

Frequency (THz)

Frequency (THz)

35

Frequency (THz)

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 37

X(0.5,0,0)

Wave Vector

Figure 4: Phonon dispersions of dierent sizes M2 termination surface slabs

26


0

0

8

8

16

16

24

24

32

32

40

40

48

48

56

56

M1 M2 M2-II

64 72

0.0%

O O-II

4.3% -8.50%

0.00% -7.90%

64

SiO SiO-II

O2 O2-II -8.50%

0.00%

Displacement drel (%)

72

0.00%

Figure 5: Displacement oscillation pattern of the relaxed atom layers

M

gO

re cip

ita

0

-4

ita

tio

ter

m.

n

Sla b

Oxygen Releasing

itat

-2 O2-II term. Slab

∆µO(eV)

rec ip

n

O2

ip

ion

SiO

ter

Si O

m. Sl

ab

-II

te rm

.S

la

O

b


-6

P=1010 bar

P=0.2

Mg Precipitation

Si P

tio

2 Pr ec

bar

P=1 b

P=1

0 -10

-II

te

rm

.S

P=

r

bar

10 -20

P=10 -5

la b

r

0

Mg2SiO4 Bulk Stable Zone

-10

-8

-6

-4

∆µMg(eV)

-2

0

2

0

250

500

750

1000

Τ(K)

1250

-2

-4 ba

P= 10 -

M2 term. Stable Zone

bar

ba

P= 10 -3

SiO-II term. Stable Zone

0

P=10 5 ba

ar

r

∆µO(eV)

2P

M1 term. Slab

O

M2-II term. Slab

Si

2

O term. Slab

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Interatom Layers

Page 27 of 37

-6

40

ba r

1500

1750

2000

Figure 6: Phase diagram: the regions of stability of the M g2 SiO4 (010) surfaces with dierent terminations as functions of the chemical potential variations for Mg and O atoms

27



T = 300K, P = 1 bar

10

M M 1t 2er II m te . rm .

4

O Si

2t er

rm

te

.

m

.

O

-II

te

rm

.

MgO precipitation

6

O

m

er It

SiO2 precipitation

Ωi (eV/Unit Cell)

8

Si precipitation

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 37

O

.

I 2O . m r te

. m

-I

O

Si

er It

2

M2 term.

0 -10

-8

-6

-4

∆µMg(eV)

-2

0

Figure 7: Surface Gibbs free energy as a function of ∆µM g at T=300 K and pO2 = 1bar

28


Page 29 of 37

10

M

Ωi (eV/Unit Cell)

8

1t

er

m . te rm .

O

-II

te

rm

.

O

2t

. r te

O

Si

M2 term.

.

-II O i S

0 -10

er m

-8

. rm . e t m II er 2- O t O

. m

4

2

MgO precipitation

2II

SiO2 precipitation

M

6

T = 2170K, P = 1 bar

Si precipitation

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-6

rm te

∆µMg(eV)

-4

-2

0

Figure 8: Surface Gibbs free energy as a function of ∆µM g at T=2170 K and pO2 = 1bar

29


The Journal of Physical Chemistry 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 34 35 36 37 38 39 40 41 42 43 44 45 46


Page 30 of 37

Page 31 of 37 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 34 35 36 37 38 39 40 41 42 43 44 45 46

M2 termination


M2-II termination

Top View Side View Mg O2 termination

O Si


O2-II termination

O termination 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 34 35 36 37 38 39 40 41 42 43 44 45 46


O-II termination

Top View Side View SiO termination

Mg O Si


SiO-II termination

Page 32 of 37

1×1×1 Slab

The Journal of Physical35Chemistry

Fvib(0K) = 2.343 eV

2×2×1 Slab

30

25

25

20

20

Frequency (THz)

30

15 10

10 5

0

0

35

3×3×1 Slab

X(0.5,0,0)

Wave Vector

Γ (0,0,0)

-5 Γ (0,0,0)

Y(0,0.5,0)

35

Fvib(0K) = 2.362 eV

30

25

25

20

20

Frequency (THz)

30

15 10

0

0

Wave Vector

Γ (0,0,0)

Wave Vector

Γ (0,0,0)

Y(0,0.5,0)

Γ (0,0,0)

Y(0,0.5,0)

Fvib(0K) = 2.362 eV

10 5

X(0.5,0,0)

4×4×1 Slab

X(0.5,0,0)

15

5

-5 Γ (0,0,0)

Fvib(0K) = 2.361 eV

15

5

-5 Γ (0,0,0)

Frequency (THz)

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 34 35 36 37 38 39 40 41 42 43 44 45 46

Frequency (THz)

35 Page 33 of 37

ACS Paragon Plus Environment -5 Y(0,0.5,0) Γ (0,0,0)

X(0.5,0,0)

Wave Vector


Interatom Layers

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 34 35 36 37 38 39 40 41

Page 34 of 37

0

0

8

8

16

16

24

24

32

32

40

40

48

48

56

56

64 72

M1 M2 M2-II 0.0%

4.3% -8.50%

O O-II 0.00% -7.90%


SiO SiO-II

O2 O2-II 0.00%

Displacement drel (%)

-8.50%

0.00%

64 72

Page 35 of 37

M

ip

ita

tio

0 rec ip

O2

n

ip

ita

tio

ter

m.

n

Sla

b

Oxygen Releasing

itat

-2

-4

Pr ec

ion

SiO

ter

m.

Si

O

Sla

b

-II

te

rm

.S

P=0.2

te

rm

.S

r

ar

0 -10

-II

b

bar

P=

bar

10 -20

P=10 -5

r

-4

10 -3

la

b

0

P=

ba

r

10 -

M2 term. Stable Zone

-6

-2

P=

SiO-II term. Stable Zone

-8

bar

ba

-6

40

ba

r

Mg2SiO4 Bulk Stable Zone

-10

0

P=10 5 ba

P=1 b

P=1

O

la


-6

P=1010 bar

Mg Precipitation

Si P

2

-4

∆µMg(eV)

-2

0

2 ACS Paragon Plus Environment

0

250

500

750

1000

Τ(K)

1250

1500

1750

2000

∆µO(eV)

re c

gO

M1 term. Slab

2P

M2-II term. Slab

O

O term. Slab

∆µO(eV)

Si

2

O2-II term. Slab

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 34 35 36 37 38 39 40 41 42 43


M

2-

8

O

2t

4

er

.

m

S

iO

II

r te

1t

te

rm

er

.

.

.

O

te

rm

.

t I I -

O

2 O . m r te

.

rm

e

t I I -

.

rm

e

m

-II

Page 36 of 37

m

MgO precipitation

6

M

SiO2 precipitation

Ωi (eV/Unit Cell)

10

Si precipitation

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47

of Physical Chemistry TThe = Journal 300K, P = 1 bar

O

Si

2

M2 term.

0 -10

-8

-6

-4


∆µMg(eV)

-2

0

Page 37 of 37

M

2-

II

8

6

1t

er

te

m

.

rm

.

O

-II

te

rm

.

O

2t

er

4

2

m

.

O

Si

. m er m. t II ter 2 O O

.

rm

te

M2 term.

.

rm

e

t I I -

O

Si

0 -10

MgO precipitation

M

SiO2 precipitation

Ωi (eV/Unit Cell)

10

T = 2170K, P = 1 bar

Si precipitation

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 34 35 36 37 38 39 40 41 42 43 44 45 46 47


-8

-6

-4


∆µMg(eV)

-2

0

Density Functional Theory Calculations and Thermodynamic Analysis

Recommend Documents