Modeling defects and plasticity in MgSiO3 post-perovskite: Part 1—generalized stacking faults

In this work, we examine the transferability of a pairwise potential model (derived for MgSiO3 perovskite) to accurately compute the excess energies of the generalized stacking faults (GSF, also called γ-surfaces) in MgSiO3 post-perovskite. All calculations have been performed at 120 GPa, a pressure relevant to the D″ layer. Taking into account an important aspect of crystal chemistry for complex materials, we consider in detail all possible locations of slip planes in the post-perovskite structure. The γ-surface calculations emphasize the easiness of glide of slip systems with the smallest shear vector [100] and of the [001](010) slip system. Our results are in agreement with previous ab initio calculations. This validates the use the chosen potential model for further full atomistic modeling of dislocations in MgSiO3 post-perovskite.


Introduction
In 2004, the discovery of the so-called "last mantle phase transition" immediately became promising to shed light on the puzzling properties of the D″ layer and to provide new insights into our understanding of the dynamics of the lowermost part of the Earth's mantle. The transformation of the magnesium silicate perovskite (bridgmanite) into a denser, also affected. Atomic diffusion of Mg 2+ and Si 4+ in postperovskite is also extremely anisotropic, with almost eight orders of magnitude difference between the fast and slow directions (Ammann et al. 2010). Metadynamics trajectories identified preferential plane sliding involving the formation of stacking faults as the source of an easy pathway for the phase transition (Oganov et al. 2005; see also Zahn 2011 for the role of shear and stacking faults on the transition between the perovskite and post-perovskite structures). This raises the question of the response of the post-perovskite structure to shear. This question has profound implications on the plastic properties of the phase which in turn may determine seismic anisotropy (Nowacki et al. 2013) and affect the flow at the core-mantle boundary (Nakagawa and Tackley 2011).
The calculation of generalized stacking faults (GSF, also called hereafter γ-surface) represents a first powerful approach to crystal plasticity (Vitek 1968). By application of a rigid-body shear to a crystal structure followed by excess energy calculation, it is possible to identify easy shear paths. γ-Surfaces also represent a key ingredient in dislocation core modeling with the Peierls-Nabarro model. This approach has been applied to post-perovskite by Carrez et al. (2007) where γ-lines were calculated for several potential slip systems demonstrating a significant plastic anisotropy on post-perovskite revealed either by the γ-surfaces or by the dislocation modeling. However, both calculations were suffering limitations which prevented from fully clarifying the issue of plastic anisotropy in post-perovskite. These limitations are illustrated by the fact that distinct models are found for screw dislocations depending on the glide plane considered (see as an illustration the case of [001] dislocations in Carrez et al. 2007). This is a result of application of the Peierls-Nabarro model which searches for a solution corresponding to planar core spreading exclusively. Solving the issue of the actual (potentially 3D) core structure of screw dislocations (which may eventually constrain plastic anisotropy) requires further modeling techniques. In this series of papers, we propose to apply full atomistic modeling to establish the 3D structure of dislocations cores in postperovskite without any assumption. However, such an approach requires modeling of large atomic systems which are beyond computational capabilities as long as firstprinciple calculations are used. Alternatively, atomistic calculations can be performed on large systems containing defects if atomic interactions are calculated from empirical potentials. However, the use of these potentials fitted against equilibrium properties must be validated in case of plasticity studies which involve atomic configurations far from equilibrium (see for instance, Lu et al. 2000;Zimmerman et al. 2000;Godet et al. 2003;Carrez et al. 2008;Ryu et al. 2009).
Numerous theoretical studies have already showed interatomic potential modeling to be very effective for high-pressure and high-temperature simulations of the perovskite phase (Oganov et al. 2000;Ito and Toriumi 2010;Chen et al. 2012;Gouriet et al. 2014;Hirel et al. 2014). In the present work, the ability of this parameterization to reproduce the structure, elastic properties and γ-surfaces of the MgSiO 3 post-perovskite is examined (by comparing with previous calculations based on first principles). Having done that, we propose a detailed investigation of the influence of crystal chemistry on plastic shear properties of post-perovskite. This study is intended to serve as a basis for the full atomistic modeling of dislocations that will follow.

Computational methods
In this study, we focus on post-perovskite with the pure MgSiO 3 composition. Indeed, Metsue and Tsuchiya (2013) have shown that the introduction of iron has little effect on the γ-surfaces without changing the plastic anisotropy style of this phase. The classical semi-empirical approach is applied for the theoretical study of the MgSiO 3 postperovskite. Force field used in the following is based on an interatomic pairwise potential (1) which takes into account long-range and short-range interactions through Coulombic and Buckingham forms, respectively. The short-range interactions include repulsive and attractive van der Waals interactions: where r ij corresponds to the distance between ions with charges z i and z j ; b ij , ρ ij and c ij are constant parameters describing the short-range interactions. The parameterization used in this work (Table 1) was previously derived by Oganov et al. (2000) for MgSiO 3 perovskite (bridgmanite).
Molecular static simulations are carried out at P = 120 GPa using the program packages GULP (Gale and Rohl 2003) for optimization of ground-state structural parameters and elastic properties, and LAMMPS (Plimpton 1995) for calculation of γ-surfaces. Both codes rely on Ewald summation methods for Coulombic interactions. Optimization is performed using a conjugate gradient algorithm with a stopping tolerance for force on atoms of 10 −9 eV/Å (1.602×10 −1 -18 N). The athermal elastic constants C ij are calculated as the second derivatives of the energy density with respect to external strain using the procedure fully implemented in GULP code. The orthorhombic symmetry of post-perovskite results in the elastic stiffness tensor containing nine independent C ij coefficients. For the evaluation of the excess energies associated with a generalized stacking fault, we use fully periodic supercells which are built on three lattice vectors a 1 , a 2 and a 3 . The supercells are oriented in such a way that the (hkl) generalized stacking fault plane of interest is normal to the Cartesian z direction (Fig. 1) (110) and (011) calculations, monoclinic atomic configurations with angles 90.62° and 92.57° are designed. It should be noted that in complex materials there could be several non-equivalent geometric locations of a slip plane with a given (hkl) index which distort different atomic bonds. All possible configurations of slip planes in the post-perovskite structure are thus considered in detail in the Results section. The optimum size of a simulation supercell along the z direction is found to be about 100 Å, which is large enough to allow relaxation of atoms close to the slip plane and to prevent the effect of atomic distortions on the boundary condition. Typical orthorhombic simulation supercell contains about 2000 atoms, while a monoclinic supercell contains about 3000 atoms. The excess energies γ, corresponding to the energy cost resulting from a rigid-body shear, are calculated by displacing the upper part of a supercell over the lower part across a chosen (hkl) slip plane applying a shear displacement f = e 1 a 1 + e 2 a 2 ( Fig. 1). Using a 5 % increment of displacement e i , the resolution of the energy landscape corresponds to 400 calculation points. In order to keep a periodic boundary condition in an atomic system containing a stacking fault, the lattice vector a 3 is given an additional component along the same vector f (Fig. 1). After shearing the upper part, atomic relaxations are allowed along the direction perpendicular to the glide plane, i.e., along z, exclusively in order to avoid spurious recovery of the perfect crystal geometry during energy minimization at smallest shear step. The energy minimization is performed at constant volume, corresponding to the ground-state volume of a perfect crystal under the given external pressure field.
γ-Surface excess energies can be used to evaluate the ideal shear stress (ISS), i.e., the upper limit of the stress that a perfect crystal can sustain (Paxton et al. 1991). ISS is then computed as the absolute maximum of the γ-surface energy derivatives relative to the applied shear.

Structural parameters and elastic properties
An orthorhombic MgSiO 3 unit cell is optimized at P = 120 GPa and T = 0 K using GULP package (Gale and Rohl 2003). The calculated structural parameters and elastic constants are given in Table 2. The a and b parameters computed with the potential model are slightly overestimated in comparison with experimental data and DFT calculations, while the c parameter is a bit underestimated. In general, all structural parameters (including the unit cell volume) are in good agreement with previous experimental and theoretical data since differences do not exceed 3 %.
The elastic C ij tensor provided by the potential model (Table 2) is found to be in a reasonable agreement with the available literature data. Generally, all non-diagonal components are about 20 % stiffer, while all diagonal  (Table 2).

Generalized stacking faults
The C-lattice of the post-perovskite results in four potential shear vectors:

Location of the stacking faults
In contrast to metals, most minerals have complex crystal structures with interatomic bonds of significantly different nature. This fact requires a detailed consideration of bonding and makes the aspect of crystal chemistry to be dominant over the simple concept of close-packed atomic planes for predicting the most probable slip planes. Thus, for each (hkl) plane, we consider different possibilities of its location along z hkl 1 in order to analyze the effect of different types and number of bonds involved.
Atomic planes parallel to (100) are the most densely packed planes in the post-perovskite structure. All atoms in the structure are located at two levels along the shortest [100] direction: z 100 = 0 and z 100 = 0.5 (Fig. 2a). Consequently, for the (100), there is only one cutting level located at z 100 = 0.25 (equivalent to z 100 = 0.75 according to the symmetry of the structure).
The slip planes (010), (110) and (011) have three nonequivalent potential geometric locations for shear (Fig. 2b,d,e). In the (010) Table 2 Lattice parameters (in Å) and elastic constants (in GPa) for MgSiO 3 post-perovskite calculated in this study using the potential model by Oganov et al. (2000) with comparison with data from previous studies Differences with DFT studies are indicated in brackets. Superscripts refer to the following studies: 1 Oganov and Ono (2004) affected bonds per unit cell for each of the listed cutting levels is given in Table 3.

GSF excess energies
Previous ab initio studies of GSF in MgSiO 3 post-perovskite provided only γ-lines (Carrez et al. 2007;Metsue and Tsuchiya 2013). In this work, full γ-surfaces are calculated for all possible slip planes and cutting levels. The shape of the obtained γ-surfaces clearly reflects the symmetry of the structure in a given glide plane for all configurations considered in this work (Fig. 3a- (001) and (110) γ-surfaces ( Fig. 3a, b, h-j) such as all energy minimum valleys and maximum peaks are translated by ½ <110>. For (010) with z 010 = 0.7 and (110) with z 110 = 0.65 γ-surfaces, there are metastable stacking faults (Fig. 3c, h) whose presence is not caused by additional translation vectors of the lattice but by the favorable atomic arrangement in these configurations of faulted crystals. In case of (010), the stacking fault is located exactly in the middle of the γ-surface (Fig. 3c) and corresponds to a highsymmetry configuration where atoms are arranged in such a way that slip and unslip halfcrystals represent a mirror reflection of each other (due to the superposition of plane n 010 and ½<101> shear component). For the (110) plane (Fig. 3h), stacking faults are centered with respect to [001] and disposed at 1/3[110] with ½<110> periodicity. Rigid shear by 1/3[110] results in a stacking fault structure where Si-octahedra are connected by corners instead of edges along the slip plane. This geometry reproduces the faulted structure related to post-perovskite-perovskite transformation previously described by (Oganov et al. 2005 Fig. 3c, h. Among all examined γ-surfaces, the (010) plane with z 010 = 0.7 (Fig. 3c) is characterized by the lowest excess energies, while the highest energy barriers are observed for (100) and (001) at z 001 = 0.47 along the longest [010] direction (Table 3; Fig. 3a, k). Glide along the [100] and ½<110> directions are characterized by a single maximum peak, while glide for all    the other systems exhibit a camel hump shape with a local minimum at 50 % along the shear vector ( Fig. 4a- (Fig. 4e). As noted above, in (001) translation of atoms and, consequently, of GSF excess energy, maximum peaks and minimum valleys occur also by ½[110] (Fig. 3a, b) (Table 3).

Validity of the empirical potential
Considering the close lattice parameters obtained from interatomic potential with DFT and comparable results in major compressional and shear elastic constants (Table 2), we feel this set of parameters has the level of transferability from Pv to Ppv, although we do expect the differences between force field calculations and density functional calculations, since an insurmountable task is in practice to use simplified classical mechanics model to simulate quantum systems on bond formation energies and bond geometry. However, modeling plasticity and defects requires transferability of the potential to reproduce atomic configurations far from equilibrium where the potential is fitted. In order to validate the chosen parameterization for dislocation modeling at atomic scale, its ability to reproduce ab initio calculations of γ-lines is tested (Zimmerman et al. 2000;Godet et al. 2003;Ryu et al. 2009;Jelinek et al. 2012;Pizzagalli et al. 2013).
The γ-surface energies computed in this work for the (010), (001), (011) and (110) glide planes are in good agreement with previous ab initio simulations (Carrez et al. 2007;Metsue and Tsuchiya 2013) (Table 3; Fig. 4a-h). All DFT γ-lines for (010) and (001) are in agreement with those obtained with the semi-empirical potential at cutting levels leading to the lowest energy barriers, i.e., to z 010 = 0.7 and z 001 = 0.65, respectively (Fig. 4a, b, d-f). Discrepancies observed for [100](010) and [100](001) slip systems may be attributed to either the underestimation of C 55 and C 66 elastic coefficients predicted by the potential model (Table 2) or to an overestimation of excess energies in DFT calculations resulting from a size effect of the volume used in DFT approach. For (011) and (110), γ-lines from DFT correspond to different cutting levels. Thus, the [100](011) γ-line by Carrez et al. (2007) is related to z 011 = 0.42, while the γ-line by Metsue and Tsuchiya (2013) reproduces the shape of the high-energy barrier obtained with z 011 = 0.49 (Fig. 4c). Considering γ-lines for (110), the ½[110](110) asymmetric curve by Metsue and Tsuchiya (2013) is consistent with our lowest energy barrier obtained with z 110 = 0.65, while results by Carrez et al. (2007) agree with the high-energy symmetric curve obtained at z 110 = 0.42 (Fig. 4g). However, the easiest cutting level is not the same for ½[110] and [001] directions: along [001], the lowest energy barrier in (110) corresponds to z 110 = 0.42 (Fig. 4h). This phenomenon is considered below in detail.
For (100), the potential model leads to γ-lines with systematically higher energies and comparatively dissimilar shape (Table 3; Fig. 4i, j). The stacking fault along [001] (100) which appears in DFT curves is barely visible with the potential (Fig. 4j). Similarly, along [010](010), the shoulders at 25 % of shear found with the DFT are not reproduced by the potential (Fig. 4i). However, both DFT and potential model clearly show that (100) is one the most unfavorable planes for shear which involves very highenergy configurations.
Apart from the discrepancy observed for very unfavorable shear planes like (100), results on the γ-surfaces obtained from the potential model compare well with ab initio calculations. The empirical parameterization of pairwise potentials proposed by Oganov et al. (2000) is thus valid to model shear properties of post-perovskite and in particular crystal defects modeling involving atomic positions far from equilibrium.

Influence of crystal chemistry
Analyzing the computed γ-lines from different cutting levels for each (hkl) plane, it is found that the lowest and the highest energy barriers do not necessarily correspond to the same cutting levels for shear along different directions within the same plane. For instance, in the (010) (010), the highest energy barriers correspond to different cutting levels (Fig. 4a, d). In (010), shear is easier at z 010 = 0.55 compared to z 010 = 0.4 along [001], whereas along [100], the energy barrier is higher at z 010 = 0.55, although the difference with z 010 = 0.4 is very small (Fig. 4a, d; Table 3). A similar effect is observed for shear along ½ [110] and [001] in (110)  But in this case, the easiest cutting level is not the same for both directions (Fig. 4g, h; Table 3). This behavior illustrates that the number and nature of bonds is not enough to describe the easiness of shear. The ability for building up new bonds after shear is likely to play an important role.

Searching for the easiest slip systems
In order to further describe the relative resistance to shear of the potential slip systems, the ideal shear stresses (ISS) are calculated (Table 3). ISS represents the upper limit of stress that a perfect crystal can sustain before yielding (Paxton et al. 1991 (Table 3). In order to analyze an effect of the reduced C 55 and C 66 values predicted by the potential model on the calculated ISS for these slip systems, the ISS values can be normalized by the corresponding C ij components and compared with those from DFT. Both DFT and semi-empirical values provide the same ISS/C ij ratios (0.25, 0.3 and 0.19 for [001](010), [100](001) and [100] (010), respectively). This correlation between computed elastic constants and ISS indicates that underestimated C 55 and C 66 values inevitably result in smaller ISS values (~35%) for [100] slip systems predicted by the potential model (Table 3). However, even keeping in mind such an effect of reduced ISS for [100] systems, they clearly remain among the most probable slip systems together with [001] (010) what is consistent with the ab initio studies by Carrez et al. (2007) and Metsue and Tsuchiya (2013). The empirical potential model also shows that ½[110](110) (at z 110 = 0.65) and [001](110) (at z 110 = 0.42) are further possible slip systems. While considering all possible cutting levels in the post-perovskite structure, we show that the lowest and the highest energy barrier can correspond to different shear levels for different directions within the same glide plane. This effect was not taken into account in previous ab initio studies and therefore has affected the (011) and (110) γ-surface calculations. Possible relevance of the ½[110](110) system is in agreement with slip predicted by Oganov et al. (2005) based on first-principle metadynamics.

Implications
In MgSiO 3 post-perovskite, plastic anisotropy represents a potential cause for seismic anisotropy in the D″ layer through the development of crystal-preferred orientations (CPO) during plastic flow. Experimental CPO produced in the diamond anvil cell with MgSiO 3 post-perovskite and its close structural analog MgGeO 3 yield conflicting results. According to Merkel et al. (2006Merkel et al. ( , 2007, both MgSiO 3 and MgGeO 3 post-perovskite phases exhibit slip on (110) and/ or (100). High-pressure experiments on the same phases by Miyagi et al. (2010;2011) reveal strong (001) texture development during the perovskite → post-perovskite transformation and subsequent slip on (001) in the [100] direction in the transformed post-perovskite phase. However, the present study rather emphasizes the role of slip in (010). Dislocation modeling will be performed using the potential validated here to further constrain the easy slip systems in MgSiO 3 post-perovskite.

Conclusions
An effective empirical parameterization of pairwise potentials (Oganov et al. 2000) developed for MgSiO 3 perovskite (bridgmanite) was applied to the next highpressure phase-post-perovskite. We have examined its transferability to reproduce the structure, elastic properties and γ-surface excess energies of post-perovskite. The significant structural, elastic and energetic parameters computed with the potential model were shown to compare well with available theoretical and experimental data.
Calculations of the γ-surface energies were performed taking into account all possible non-equivalent geometric locations of slip planes breaking different bonds in the post-perovskite structure. The observed phenomenon of switching the cutting level with lowest energy barrier for different slip directions within one glide plane (see [001] (110) and ½[110](110) systems) highlights the necessity to consider all possible cutting levels for complex materials. In agreement with previous ab initio studies (Carrez et al. 2007;Metsue and Tsuchiya 2013), the slip systems involving the shortest [100] shear vector and the [001](010) system within the glide plane cutting only Mg-O bonds are characterized with very low γ-surface energies. Based on current γ-surface calculations, (110) could be considered as further possible slip plane.
Reasonable agreement with previous theoretical and experimental results proves the chosen potential parameterization (Oganov et al. 2000) to be valid for further fully atomistic modeling of dislocations and their mobility in MgSiO 3 post-perovskite.