# Dislocation modelling in Mg2SiO4 forsterite: an atomic-scale study based on the THB1 potential

S. Mahendran, P. Carrez, S. Groh, Philippe Cordier

► To cite this version:

S. Mahendran, P. Carrez, S. Groh, Philippe Cordier. Dislocation modelling in Mg2SiO4 forsterite: an atomic-scale study based on the THB1 potential. Modelling and Simulation in Materials Science and Engineering, 2017, 25 (5), pp.1-19. 10.1088/1361-651X/aa6efa. hal-01606755

HAL Id: hal-01606755

https://hal.science/hal-01606755

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Dislocation modelling in Mg₂SiO₄ forsterite: an atomic-scale study based on the THB1 potential

To cite this article: S Mahendran et al 2017 *Modelling Simul. Mater. Sci. Eng.* **25** 054002

View the [article online](#) for updates and enhancements.

## Related content

- Modelling [1 0 0] and [0 1 0] screw dislocations in MgSiO₃ perovskite based on the Peierls-Nabarro-Galerkin model
K Gouriet, Ph Carrez and P Cordier

- Peierls-Nabarro modelling of dislocation in MgO
Philippe Carrez, Denise Ferré and Patrick Cordier

- Dislocation motion, constriction and cross-slip in fcc metals
M S Duesbery
# Dislocation modelling in Mg₂SiO₄ forsterite: an atomic-scale study based on the THB1 potential

S Mahendran<fnref n="1" />, P Carrez<fnref n="1" />, S Groh<fnref n="2" />,<fnref n="3" /> and P Cordier<fnref n="1" />

E-mail: philippe.carrez@univ-lille1.fr

Received 16 January 2017, revised 29 March 2017

Accepted for publication 24 April 2017

Published 23 May 2017

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

Knowledge of the deformation mechanisms of (Mg,Fe)₂SiO₄ olivine is important for the understanding of flow and seismic anisotropy in the Earth's upper mantle. We report here a numerical modelling at the atomic scale of dislocation structures and slip system properties in Mg₂SiO₄ forsterite. Our study focuses on screw dislocations of [100] and [001] Burgers vectors. Computations are performed using the so-called THB1 empirical potential set for Mg₂SiO₄. Results of dislocation core structures highlight the primary importance of the (010) plane for [100] slip dislocations. For [001] dislocations, we confirm the occurrence of a stable narrow core that evolves into transient planar configurations to glide in (100) and (010). Such configurations suggest a locking-unlocking mechanism.

Keywords: atomistic simulation, ionic materials, core shell model, olivine, dislocations, crystal plasticity

(Some figures may appear in colour only in the online journal)

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 3.0 licence](https://creativecommons.org/licenses/by/3.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Univ. Lille, CNRS, INRA, ENSCL, UMR 8207, UMET—Unité Matériaux et Transformations, F-59000 Lille, France</footnote>

<footnote n="2">²Institute of Mechanics and Fluid Dynamics, TU Bergakademie Freiberg, D-09599 Freiberg, Germany</footnote>

<footnote n="3">³Department of Biomedical Engineering, University of Basel, CH-4123 Allschwill, Switzerland</footnote>
# 1. Introduction

It is nowadays well accepted that convection in the upper mantle is constrained by the plastic behaviour of olivine (Mg,Fe)₂SiO₄, one of its main constituents. Seismological observations, such as seismic anisotropy, can be related to the development of lattice-preferred orientations, and inform us about flow in the mantle (e.g. Panning and Romanowicz 2004). However, as pointed out by Mainprice and co-workers (e.g. Mainprice et al 2005) or more recently by Trampert and co-workers (de Wit and Trampert 2015), seismological observations remain as indirect evidence of the deformation mechanism and must be handled with care. Due to its primary importance in geology, the plastic deformation of olivine has been the subject of numerous experimental studies since the mid-seventies (Phakey et al 1972, Durham et al 1977, Evans and Goetze 1979, Darot and Gueguen 1981, Gaboriaud et al 1981, Gueguen and Darot 1982, Wang et al 1988, Barber et al 2010, Mussi et al 2014). Intracrystalline plastic deformation of olivine involves two types of dislocations characterised by either [100] or [001] Burgers vectors (here and throughout the study we refer to the *Pb* nm description for the crystal structure). The determination of slip planes is more uncertain, as several planes have been reported in deformation experiments involving either single crystals or poly-crystalline material. Under conditions of high stress and low temperature, glide of [001] dislocations is usually reported in (100), (010), or {110}. At high temperature, [100] dislocations gliding in (010), (001), {011}, and potentially {021} or {031}, dominate. A high lattice friction is suggested in olivine by dislocation microstructures since the dislocations tend to align along preferential crystallographic orientations. The dislocations with [001] Burgers vectors usually exhibit straight screw segments (Phakey et al 1972, Durham et al 1977, Evans and Goetze 1979, Darot and Gueguen 1981, Gaboriaud et al 1981, Wang et al 1988). At low temperature, [100] dislocations gliding in (001) planes are also characterised by long straight screw segments (Durham et al 1977). For temperature above 1200 °C, transmission electron microscope (TEM) investigations (Durham et al 1977, Gueguen and Darot 1982) show that, in (010), [100] dislocations exhibit long straight edge segments while the screw components appear now to be much shorter and curved.

More recently, numerical modelling of dislocations in pure-magnesium end-member forsterite Mg₂SiO₄ was carried out by Durinck et al (2007). In this work, dislocation core structures have been derived using the Peierls–Nabarro (PN) model (e.g. Peierls 1940 or Joós et al 1994). Although undertaken using *ab initio* calculations of generalised staking fault (GSF) energy, such core structures investigation suffers from the intrinsic limitations of the PN model (Schoeck 2005). In particular, using the PN model, Durinck et al (2007) could only study planar core configurations. Using an empirical potential parameterisation for Mg₂SiO₄ forsterite, Walker et al (2005) computed the atomic arrangement around a screw dislocation line of Burgers vector [001]. By comparison between the two approaches, it was later concluded that [001] dislocations exhibit a non-planar core (comparable to that found in bcc metals for instance), which is consistent with the occurrence of the long [001] screw segment observed in TEM images (Carrez et al 2008). Nevertheless, any conclusion regarding the dislocation slip systems in olivine based on these previous numerical modellings only is speculative.

In the following, we propose to model the core structures of dislocations with Burgers vector [100] and [001] in Mg₂SiO₄ olivine by performing atomic-scale calculations. As further described in the following section 2, we rely on the so-called THB1 potential (Price and Parker 1987) specifically derived to model forsterite. In order to account for polarisation effects in the material, the potential is parameterised with the inclusion of a core–shell (C/S) assumption. Since molecular static calculations are performed using the large-scale atomic/
<figure><img src="image_4.png" /><figcaption>Figure 1. Mg₂SiO₄ structure with the unit cell (black lines). The structure is described using the Pbnm space group. Magnesium atoms are shown in yellow, oxygen in red, and silicon in blue, located at the centre of the SiO₄ tetrahedra. The silicon and magnesium M2 lies on the mirror plane, and magnesium M1 on the inversion centre.</figcaption></figure>

molecular massively parallel simulator (LAMMPS) (Plimpton 1995), section 2 addresses how we implemented the THB1 potential for LAMMPS simulation. The atomic core structures and Peierls stresses for the two kinds of dislocations of interest are computed according to the methodology presented in section 3.

## 2. Implementation of the THB1 potential within the LAMMPS package

THB1 is an empirical potential developed to model Mg₂SiO₄ forsterite and its high-pressure polymorphs (Price and Parker 1987). It includes a long-range electrostatic term described through a classical Coulombic energy term, a short range repulsive term based on a Buckingham form, a harmonic three-body term to model the directionality of Si–O bonding, and a core–shell approximation to include electronic polarizability of ions. For computational efficiency, we implemented the THB1 potential into a new inter-atomic potential library freely available for the LAMMPS package. Our implementation currently relies on the following energy decomposition:

$$U_{\text{Total}} = U_{\text{Coulomb}} + U_{\text{Short}} + U_{\text{Spring}} + U_{\text{THB}} \quad (1)$$

The long-range Coulombic energy is calculated by application of a Wolf summation (Wolf et al 1999), which allows calculations on periodic or non-periodic systems. The electronic polarizability is modelled in the shell model approximation by considering oxygen ions as a massive core (i.e. the nucleus) connected by a harmonic spring to a massless shell (i.e. the electron cloud). Because of this approximation, the charge of the shell remains constant and independent of the deformation state. This additional degree of freedom enables the ion to polarise in response to an electric field due to surrounding ions. The coupling between a core and its corresponding shell is described by a harmonic spring term, $U_{\text{Spring}}$ defined as follows:

$$U_{\text{Spring}} = \frac{K_1 r_i^2}{2} + \frac{K_2 r_i^4}{24}, \quad (2)$$

where $K_1$ and $K_2$ are the shell spring constants, and $r_i$ is the separation distance between the core and its corresponding shell. The short-ranged Buckingham potential, $U_{\text{Short}}$, describes
<table><caption>**Table 1.** Parameters for the THB1 potential used to model Mg₂SiO₄ forsterite (Price and Parker 1987).</caption><thead><tr><th colspan="4">Charge (units of electronic charge)</th></tr><tr><th>Ions</th><th>Core</th><th>Shell</th><th>Core–shell spring constant K1 (eV Å<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Mg</td><td>2.0</td><td></td><td></td></tr><tr><td>Si</td><td>4.0</td><td></td><td></td></tr><tr><td>O</td><td>0.84819</td><td>-2.84819</td><td>74.92</td></tr><tr><th colspan="4">Buckingham potential</th></tr><tr><th></th><th>A(eV)</th><th>B (Å)</th><th>C (eV Å<sup>6</sup>)</th></tr><tr><td>Mg–O</td><td>1428.5</td><td>0.2945</td><td>0.0</td></tr><tr><td>Si–O</td><td>473.2</td><td>0.4157</td><td>0.0</td></tr><tr><td>O–O</td><td>22764.3</td><td>0.1490</td><td>60.08</td></tr><tr><th colspan="4">Harmonic three-body term</th></tr><tr><th></th><th>k<sup>B</sup>(eV rad<sup>-2</sup>)</th><th></th><th>θ<sub>0</sub> (degrees)</th></tr><tr><td>O–Si–O</td><td>2.09</td><td></td><td>109.47</td></tr></tbody></table>

<table><caption>**Table 2.** Elastic constant $C_{ij}$ tensor (in GPa) using either Wolf summation method (this study) or Ewald summation (Gale 1997). First-principles calculations (Durinck *et al* 2005) are given for comparison with the THB1 potential.</caption><thead><tr><th></th><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>33</sub></th><th>C<sub>44</sub></th><th>C<sub>55</sub></th><th>C<sub>66</sub></th><th>C<sub>12</sub></th><th>C<sub>13</sub></th><th>C<sub>23</sub></th></tr></thead><tbody><tr><td>Using Wolf summation</td><td>344.4</td><td>200.9</td><td>279.2</td><td>36.1</td><td>75.9</td><td>82.2</td><td>88.2</td><td>92.4</td><td>87.8</td></tr><tr><td>Using Ewald summation</td><td>344.9</td><td>201.7</td><td>279.4</td><td>37.0</td><td>76.1</td><td>81.7</td><td>88.5</td><td>92.3</td><td>88.5</td></tr><tr><td>First-principles calculations</td><td>300.7</td><td>195.7</td><td>224.1</td><td>65.8</td><td>78.3</td><td>77.3</td><td>61.1</td><td>62.2</td><td>64.7</td></tr></tbody></table>

the effect of the electron cloud on nearby neighbour ions. The Buckingham potential is parameterised by three constants $A_{ij}$, $B_{ij}$, and $C_{ij}$ between shells i and j. It's functional form is defined as follows:

$$U_{\text{Short}} = \sum_{ij} A_{ij} \exp\left(-\frac{r_{ij}}{B_{ij}}\right) - C_{ij} r_{ij}^{-6}. \qquad (3)$$

Finally, the directionality of Si–O bonding is included in THB1 through a harmonic three-body contribution given in equation (4):

$$U_{\text{THB}} = \sum_{ijk} k_{ijk}^{B} (\theta_{ijk} - \theta_{0})^{2}, \qquad (4)$$

where $k_{ijk}^B$ is a derivable spring constant, $\theta_{ijk}$ is the angle between the O–Si–O three-body term, and $\theta_0$ is the constant tetrahedral angle.

THB1 was previously derived by Price and Parker (1987), and table 1 summarises the parameterisation used in this work. Within this implementation, the lattice parameters of forsterite are $a = 4.7874$ Å, $b = 10.2717$ Å, and $c = 6.0227$ Å (as shown in figure 1). We also computed the nine elastic constants of forsterite, following the description of Beckstein *et al* (2001) for orthorhombic crystal structure. The calculated $C_{ij}$ are listed in table 2. It is worth pointing out that, compared to previously reported elastic constants, the Wolf summation method leads to identical results. Independent of the long-range summation techniques
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) GSF energy (J/m²) in (010). Calculations are performed by applying lattice translations in (010), where ideal relaxation conditions given by Durinck *et al* (2005) have been used (i.e. Mg and O can relax freely in every direction whereas Si atoms can move only perpendicular to the fault plane). γ-energies along [100] and [001] are shown in (b) and (c), respectively. One may notice that energies computed with the THB1 potential are higher than those of DFT calculations (Durinck *et al* 2005).</figcaption></figure>

(Ewald or Wolf), computed elastic constants are in perfect agreement with differences below 1%. Nevertheless, as commonly reported, empirical potential calculations do not fully render the first-principles elastic constants set, but correspond to the best compromise for forsterite (Price and Parker 1987). To further model dislocation structures at the atomic scale, one can also investigate the response of the potential to rather large and non-elastic displacements. This is commonly done by computing GSF. Therefore, prior to dislocation modelling, we perform GSF calculations and compare THB1 γ-energies with those computed from first-principles by Durinck *et al* (2005). For the GSF calculations, we follow the full periodic boundary conditions method described by Gouriet *et al* (2014) with Si atoms allowed to move perpendicularly to the stacking fault plane and, in order to increase the degrees of freedom of the system, with O and Mg atoms freely allowed to relax in any direction. Figure 2 shows the (010) γ-surface energy. Low-energy paths correspond to [100] and [001] (figures 2(b) and (c)), as expected from density function theory (DFT) calculations. The shapes of the GSF are
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Geometry of atomistic systems used to model the mobility of screw dislocations with cluster periodicity. The atoms in the shaded regions are frozen. The dotted lines in the middle of the simulation volume represent the trace of the glide plane.</figcaption></figure>

also consistent between the two methods, and one may notice that γ-surface energies computed with the THB1 potential are systematically slightly higher.

The analytical expressions of force fields implemented to the code and algorithm are given in the *appendix*.

**3. Dislocation core modelling methods**

The atomic systems used are designed to study dislocation core properties, and dislocation motion in different glide planes. We thus rely on a periodic system, defined as follows: the dislocation line is along *z*, and the normal to the glide plane is along *y* (see figure 3). Consequently, the motion of screw dislocation is along the *x* direction of the atomic system. Dealing with straight screw dislocation lines, the system is naturally periodic along the *z* direction (with a length equal to the magnitude of the Burgers vector, *b*, of the simulated dislocation). The periodicity along *x* is maintained by adding a tilt component of magnitude *b*/2 to this axis (Hirel *et al* 2014). For each investigated slip system, the Mg₂SiO₄ crystal structure is rotated to fulfil these two conditions. However, dealing with an orthorhombic symmetry material, the *y* direction (i.e. the normal to the glide plane) is not necessarily aligned with a crystallographic direction [*uvw*]. Therefore, the atoms located at the top and bottom regions of the supercell along the *y* axis are kept fixed to mimic an infinite perfect crystal. With such a setup, the systems investigated can be considered as periodic slabs. To further minimise spurious elastic interaction between periodic replicas, the system sizes are gradually increased to typically 480 Å between dislocation replicas, and 120 Å between the glide plane and artificially frozen regions of atoms. This corresponds to system sizes between 13 440 atoms and 53 760 atoms. Straight screw dislocations are finally introduced into the system thanks to the classical displacement field either relying on isotropic elasticity solutions or anisotropic elasticity as implemented in Atomsk (Hirel 2015). All as-built systems are then optimised using a conjugate gradient minimisation algorithm with a stopping tolerance force set at 10⁻¹² eV/Å. Following the minimisation procedure, the dislocation core structures are
<figure><img src="image_7.png" /><figcaption>**Figure 4.** [100] screw dislocation. (a) Atomic core configuration viewed along the dislocation line; yellow, blue, and red spheres represent Mg, Si, and O ions respectively; SiO$_4$ tetrahedra are shown in transparent blue; solid arrow describe the differential displacement along the [100] direction between two ions. (b) The corresponding disregistry function $S(x)$ in (010), and, shown by the dotted line, the density of [100] Burger vectors $dS/dx$ of the dissociated screw dislocation, calculated for the cationic sub-lattice.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Representation of local edge displacement around the dislocation core for the oxygen sub-lattice. Continuum displacement fields are potted on top of the relaxed atomic structure of the [100] dislocation. (a) Displacement field along the [010] direction. (b) Displacement field along the [001] direction.</figcaption></figure>

analysed using differential displacement maps (Vítek 1968) or through the computation of atomic disregistries across glide planes.

To determine the dislocations mobility, a simple shearing is applied to the system, which results in a resolved shear stress in the glide plane. For screw dislocations this can be achieved by applying shear strain increments of 1%, ensuring quasi-static loading. At each increment, the system is relaxed using a conjugate gradient method. Consequently, the aforementioned loading conditions allow for determining a critical strain, and thus a critical stress above which the system does not respond purely elastically. The dislocation starts moving when an inflection on the stress–strain behaviour is obtained. In the following, various cell sizes (along the x and y directions) have been considered to ensure system size effects. Using molecular statics, the computed critical shear stresses are the Peierls stresses for the considered slip system. By repeating such a process for various planes, the methodology can be applied to [100](010), [100](001), [001](100), [001](010), and [001](110) slip systems.

# 4. Results

## 4.1. Screw dislocations with [100] Burgers vector

Figure 4 shows the differential displacements map around a [100] screw dislocation core. Based on neighbours' cations analysis, the core of the [100] screw dislocation shows a clear tendency of spreading in (010) with a core centre located between SiO₄ tetrahedra at the M2 site (we used the classical distinction M1, M2 of cation sites in olivine group structure; see for instance Deer et al 1982). To investigate the spreading of the core, one can extract the disregistry functions within (010) from the atomic positions. We observe (figure 4) a dissociation of the core into two strongly correlated Burgers vector density peaks. The distance between these peaks, ~7 Å (table 3), is close to one lattice parameter in the [001] direction. Details of the atomic arrangements around the line indicate that the core of [100] dislocations involves some local edge displacements in the vicinity of (010) (figure 5). These local displacements do not produce any edge component of the dislocation Burgers vector but can be associated with a dilatation state of the core. Following Sun et al (2016), we use an atomistic-to-continuum cross over method to analyse the details of local edge displacements. The unrelaxed and relaxed configurations obtained during the atomistic calculations are used as the reference state and deformed state, respectively. The components of the displacement vector field, **u**, in the vicinity of the dislocation core are calculated for the three different sub-
<table><caption>**Table 3.** Dislocation disregistry analysis. Half widths $\zeta$ of the dislocation core are given according to a fit of the atomic disregistries with the following expressions: $S(x) = b/2 + b/\pi \arctan(x/\zeta)$ for a compact core, and $S(x) = b/2 + b_1/\pi \arctan(x-x_1/\zeta_1) - b_2/\pi \arctan(x-x_2/\zeta_2)$ for a dissociated core. In the case of a dissociated core, $b_i$, $x_i$, and $\zeta_i$ correspond to the partial Burgers vectors, the position of the partial cores, and their widths, respectively.</caption><thead><tr><th rowspan="2">Dislocation</th><th rowspan="2">Type of core configuration</th><th rowspan="2">Plane</th><th rowspan="2">Distance of separation of partials (Å)</th><th colspan="2">Half width</th></tr><tr><th>ζ₁(Å)</th><th>ζ₂(Å)</th></tr></thead><tbody><tr><td>[100]</td><td>Stable</td><td>(010)</td><td>7.48</td><td>2.23</td><td>2.23</td></tr><tr><td>[100]</td><td>Metastable</td><td>(010)</td><td>—</td><td>3.15</td><td>—</td></tr><tr><td>[001]</td><td>Stable</td><td>(010)</td><td>—</td><td>2.23</td><td>—</td></tr><tr><td>[001]</td><td>Stable</td><td>(100)</td><td>—</td><td>2.16</td><td>—</td></tr><tr><td>[001]</td><td>Stable</td><td>{110}</td><td>—</td><td>1.87</td><td>—</td></tr><tr><td>[001]</td><td>Transient</td><td>(010)</td><td>—</td><td>5.87</td><td>—</td></tr><tr><td>[001]</td><td>Transient</td><td>(100)</td><td>13.86</td><td>1.58</td><td>2.07</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Trace of strain tensor mapped around the [100] dislocation core region and calculated from the O sub-lattice. Positive values for the trace of the strain tensor quantify the dilatation of the dislocation core.</figcaption></figure>

However, since the compact anionic sub-lattice ensures continuity of the structure, only displacements corresponding to the oxygen sub-lattice are shown in what follows. The displacement vector field maps (figure 5) clearly indicate that the edge components are restricted to the vicinity of (010), and correspond to outward displacements of atoms. Thus, from the computed gradient of the displacement vector, **u**, one can compute the corresponding strain tensor giving access to the local dilatation taken as the trace of the strain tensor. As shown in figure 6, the core of the [100] screw dislocation is associated with positive values of the trace of the strain tensor, confirming the dilatation state of the dislocation core.

As mentioned in the methods section, a simple shearing was applied to the stable core to determine at which critical stress the dislocation starts moving in a given glide plane. We find two different critical stresses for glide in (010) and (001). For [100](010), the computed Peierls stress is 3.1 GPa. By contrast, when the [100] dislocation is forced to glide into (001) we find a higher value of the Peierls stress around 7 GPa. In the case of slip in (010), the [100] core structure evolves with stress. While it still spreads in (010), the core becomes narrower and does not show any evidence of dissociation (figure 7 and table 3). Computing the energy difference between the two [100] screw configurations, it turns out that the undissociated dislocation core corresponds to a metastable state with a higher energy of 0.4 eV $b^{-1}$ with
<figure><img src="image_10.png" /><figcaption>**Figure 7.** [100] screw dislocation. Metastable core configuration of the dislocation under a stress of 3.1 GPa. The disregistry function $S(x)$ and the [100] Burger vector density $dS/dx$ of the metastable core do not show any evidence of dissociation. Nevertheless, the core is spread in (010) with a core width of 6.3 Å.</figcaption></figure>

respect to the dissociated stable configuration (0.4 eV $b^{-1}$ is an average of energy differences calculated from simulation cells of different sizes).

## 4.2. Screw dislocations with [001] Burgers vector

Inserting a [001] dislocation into the structure leads to a stable core configuration located between four $SiO_4$ tetrahedra with a dislocation line centred on a Mg cation site (figure 8). As previously shown, the structure of the [001] stable core can be analysed with the differential displacement map. All significant displacement arrows around the dislocation line are restricted to a small area within one unit cell (figure 8). Thus, the core appeared to be compact with no particular evidence of spreading in any preferential plane. This is confirmed by computing the disregistry functions in several planes (100), (010), or {110} (figure 9). Independent of the plane, the core profile is narrow with a width smaller than 5 Å (table 3), preserving the shape of the $SiO_4$ tetrahedra.

We now investigate the critical stress needed to promote glide in these planes. We found that above a critical stress of 7.2 GPa, the dislocation moves over a large distance (with respect to the simulation cell size) in (010). Above 7.2 GPa of resolved shear stress, an analysis of the conjugate gradient minimiser steps shows that once the dislocation is set into motion, it exhibits a different core structure (figure 10), entirely planar, spread in (010). It is worth noting that such a planar configuration cannot be stabilised without applying stress in our simulations. Indeed, we verified that without stress a planar configuration always recombines into the compact stable core configuration (figure 8). However, once this transient state (figure 10) is reached, the dislocation core can travel over several lattice repeats before finally falling again into another stable core configuration. For a slip in (100), we found an onset for core motion corresponding to a critical stress of 6.9 GPa. Once again, the transient states of the core (figure 11) computed according to the conjugate gradient algorithm show that the core tends to adopt a planar configuration in (100) during motion. As for the other case, the transient planar configuration in (100) cannot be stabilised without stress. Finally, for the [001](110) slip, whereas we apply a simple shear component to maximise the resolved
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** [001] screw dislocation. (a) Atomic core configuration viewed along the dislocation line; yellow, blue, and red spheres represent Mg, Si, and O ions, respectively; SiO$_{4}$ tetrahedron are shown in transparent blue; solid arrows describe the differential displacement along the [100] direction between two ions. (b) The same core described within the oxygen sub-lattice only.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Disregistry function *S*(*x*) and the [001] Burger vector density d*S*/dx of the screw dislocation (dotted lines) for the cationic sub-lattices in (010), (100), and {110}, respectively. These disregistries are calculated along the BB', CC', and DD' profiles depicted in figure 7.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Transient core of the [001] dislocation as it glides in (010). (a) Atomic core configuration with differential displacement map superimposed; (b) disregistry function of the transient planar configuration.</figcaption></figure>

shear stress in (110), we find that above 7.8 GPa, the dislocation moves in (100) with a behaviour similar to that previously described. Indeed, due to the close orientation between (110) and (100), a simple shearing of the simulation cell cannot avoid inducing a significant resolved shear stress in (100).

**5. Discussion**

Our calculations of the atomic arrangements around the screw dislocation core of dislocations with the [001] Burgers vector show a compact core structure with no specific or preferential spreading in any given plane. As expected, differential displacement maps around the [001] core show typical features consistent with previous attempts to model screw dislocation cores
<figure><img src="image_17.png" /><figcaption>**Figure 11.** Transient core of the [001] dislocation as it glides in (100). (a) Atomic core configuration with differential displacement map superimposed; (b) disregistry function of the transient planar configuration.</figcaption></figure>

in forsterite directly from atomic-scale computation (Walker *et al* 2005, Carrez *et al* 2008). Nevertheless, the choice made for system set-ups are different. In the study by Walker *et al* (2005), computations were performed on a small system size with, consequently, a potentially strong influence on rigid periodic boundary conditions. In the present work, the investigated systems are much larger and the boundary conditions applied are supposed to be less restrictive (Woodward 2005, Walker *et al* 2010). Also, as discussed in the methods section, size effects were taken into account by gradually increasing the system size. One of the main achievements of the present study is to be able to compute Peierls stresses for glide in (100), (010), or (110). For [001](100) and [001](010), we find rather high comparable Peierls
stresses of, respectively, 6.9 GPa and 7.2 GPa. For [001](110), the determination of the Peierls stress is more ambiguous, but our results suggest that in this plane, the Peierls stress should be above 7.8 GPa. Whatever the glide plane, lattice friction opposed to [001] dislocation is therefore very high. The occurrence of high lattice friction can be viewed as a consequence of the compact core structure of [001] dislocation, and is therefore consistent with the experimental observations of long straight screw dislocation lines in deformed olivine samples (e.g. Kohlstedt et al 1976 or Gaboriaud et al 1981). One may also notice that whatever the glide plane, our results suggest that glide in (100), (010), or (110) may be activated at comparable stress levels. Indeed, glide of [001] dislocations in (100) and also possibly in {110} have been reported from TEM observations (see for instance Gaboriaud et al 1981) or using decoration techniques (Kohlstedt et al 1976). More recently, based on electron tomography techniques, the three possible slip planes, i.e. (100), (010), and {110}, have been reported simultaneously (Mussi et al 2014), but with a greater occurrence of {110}, which is not explained by our calculations of lattice friction.

At stresses higher than the Peierls stress, our calculations show that the [001] core adopts some planar transient state (figures 10 and 11). Since it is impossible to preserve the atomic configurations of these transient planar configurations when minimising the energy of the system if no external strain is applied, it is suggested that the planar core configurations for [001] dislocation are highly unstable. Since glide may involve high-energy configurations, we confirm a clear tendency of the [001] core to glide through a so-called locking-unlocking mechanism (Couret and Caillard 1989), as recently described at the atomic scale in hpc metals (Clouet et al 2015). This suggestion was initially proposed by Carrez et al (2008) based on unrelaxed core energy calculations and the PN model of core spreading. As for glide in (100) and (010), one can speculate that glide in (110) would involve a planar configuration. Although a comparison with PN model calculations may not be straightforward, it is however interesting to note that the transient core of the [001] screw dislocation in (010) does not show any evidence of dissociation whereas that in (100) indicates a tendency to dissociate (see figures 10(b) and 11(b)). Qualitatively, these two transient cores are similar to those computed using the PN model (Durinck et al 2007). This is another argument to suggest that glide in (110) would involve a planar dissociated transient core, as predicted by the PN model.

For [100] screw dislocations, our results show that the core dissociates in (010) into two collinear partial dislocations. This tendency to dissociation has never been reported; however, based on hcp lattice considerations, possible dissociations of dislocations in olivine were addressed by Poirier (1975). However, the computed width of dissociation is rather small with a distance between partial dislocations slightly larger than a lattice repeat along the [001] direction. Such a dissociation scheme can be viewed as a direct consequence of the occurrence of a stable stacking fault energy at ½[100] in the GSF calculations (figure 2(b)). Unfortunately, such a feature was not clearly observed on γ-energy computed from first-principles by Durinck et al (2005), and may be a function of the relaxation procedure used during the GSF calculation. Indeed, for DFT calculations, Durinck et al (2005) computed the GSF with non-periodic atomic systems, for which surfaces induce discontinuities and potentially oscillations in the γ-surface energies depending on cell sizes. As a consequence, the dissociation configurations of the [100] core may not be captured by the PN model used by Durinck et al (2007). More importantly, the THB1 potential calculations of the [100] dislocation core show local edge displacement of atoms within the vicinity of the core perpendicular to (010), which is in agreement with displacements of atoms outward from (010) observed by Durinck et al (2005) along the γ-energy path. Therefore, the present calculations confirm the occurrence of the dilatation state within the [100] screw core. In the context of the upper Earth's mantle (where olivine is one of the most important minerals until
a depth of ca. 410 km, where the pressure reaches 13 GPa), the dilatation state of the core may have strong implications on the activation of the [100] slip with respect to an increase in pressure. Indeed, as proposed by Durinck et al (2005), one effect of pressure is to work against this dilatation state, leading to an intrinsic hardening of the glide properties of [100](010).

## 6. Conclusion

Molecular static calculations relying on the THB1 potential have been undertaken to study dislocation core structures and slip in Mg₂SiO₄ forsterite. We analyse the dislocation core structures of screw dislocations with [100] and [001] Burgers vectors. For [100] dislocations, our results highlight the importance of the [100](010) slip system resulting from a clear spreading of the dislocation core in this plane. By contrast, the [001] screw dislocation is characterised by a compact core. Moreover, glide analysis of [001] reveals the occurrence of several planar core configurations depending on the slip planes. It should be emphasised that such configurations can be viewed as higher energy transient states of the core related to the dislocation motion typical of a locking-unlocking mechanism. The calculation of the energy of the transient planar configurations is out of the scope of this study, but could be addressed in the future. As a final remark, we believe that the methodology presented in this study can be used to further investigate the effect of pressure on the core structure in olivine, and further revisit the unsolved question of the change in main slip system activities from [100] to [001] in olivine, as reported over the last decades.

## Acknowledgements

This work was supported by funding from the European Research Council under the Seventh Framework Programme (FP7), ERC grant N°290424—RheoMan. Computational resources have been provided by the CRI-Université de Lille 1. The authors are grateful for fruitful discussions with Dr Xiao-Yu Sun and his help in the implementation of the atomistic-to-continuum cross over method. SG acknowledges the MIRACLE Project at the University of Basel funded by the Werner Siemens Foundation, Zug/Switzerland.

## Appendix. Implementation details

The THB1 potential library consists of an electrostatic term, a short-range Buckingham potential, a harmonic spring term, and a three-body term. The expressions required for the calculation of force fields are obtained from the negative gradient of the potential energy (1). In summary, the force can be calculated as

$$f = -\nabla U_{\text{Total}}. \qquad (5)$$

Analytical expressions for the first derivative of the electrostatic term have been described by Wolf et al (1999), and are given as follows:
$$f_{\text{Coulomb } \alpha} = \sum_{\substack{j \neq i \\ (r_{ij} < R_c)}} q_i q_j \left\{ \left( \frac{\text{erfc}(\alpha R_c)}{r_{ij}^2} + \frac{2\alpha}{\pi^{1/2}} \frac{\exp(-\alpha^2 r_{ij}^2)}{r_{ij}} \right) \times \frac{r_{ij\alpha}}{r_{ij}} - \left( \frac{\text{erfc}(\alpha R_c)}{R_c^2} + \frac{2\alpha}{\pi^{1/2}} \frac{\exp(-\alpha^2 R_c^2)}{R_c} \right) \times \frac{r_{ij\alpha}}{R_c} \bigg|_{r_{ij}=R_c} \right\} \quad (6)$$

with the force field corresponding to the two-body Buckingham potential expressed as

$$f_{\text{Short } \alpha} = \left( -6C_{ij}r_{ij}^{-6} + \frac{A_{ij}}{B_{ij}} \exp\left(-\frac{r_{ij}}{B_{ij}}\right) \times r_{ij} \right) \times \frac{r_{ij\alpha}}{r_{ij}^2}, \quad r_{ij} < R_B, \quad (7)$$

and the harmonic spring term as

$$f_{\text{Spring } \alpha} = -\left(K_1 + \frac{K_2 r_{ij}^2}{6}\right) r_{ij\alpha}, \quad r_{ij} < R_{\text{Spring}}. \quad (8)$$

It is not straightforward to write as a single term to express the force field arriving from the covalent three-body term. Hence, the force experienced by atom i due to interactions with atoms j and k are written as follows:

$$f_{\text{THB } i\alpha} = -f_{\text{THB } j\alpha} - f_{\text{THB } k\alpha}, \quad (9)$$

$$f_{\text{THB } j\alpha} = A_{11}r_{ij\alpha} + A_{12}r_{ik\alpha}, \quad (10)$$

$$f_{\text{THB } k\alpha} = A_{22}r_{ik\alpha} + A_{12}r_{ij\alpha}, \quad (11)$$

where

$$\Delta\theta = (\theta_{ijk} - \theta_0), \quad (12)$$

$$A = -2 \Delta\theta \sin \theta_{ijk}, \quad (13)$$

$$A_{11} = \frac{A \cos \theta_{ijk}}{r_{ij}^2}, \quad (14)$$

$$A_{12} = \frac{A}{(r_{ij} \times r_{ik})}, \quad (15)$$

$$A_{22} = \frac{A \cos \theta_{ijk}}{r_{ik}^2}. \quad (16)$$

## References

Barber D J, Wenk H R, Hirth G and Kohlstedt D L 2010 Dislocations in minerals *Dislocations in Solids* ed J P Hirth and L Kubin vol 19 (Oxford: Elsevier) pp 171–232

Beckstein O, Klepeis J, Hart G and Pankratov O 2001 First-principles elastic constants and electronic structure of α-Pt₂Si and PtSi *Phys. Rev. B* **63** 134112

Carrez P, Walker A M, Metsue A and Cordier P 2008 Evidence from numerical modelling for 3D spreading of [001] screw dislocations in Mg₂SiO₄ forsterite *Phil. Mag.* **88** 2477–85

Clouet E, Caillard D, Chaari N, Onimus F and Rodney D 2015 Dislocation locking versus easy glide in titanium and zirconium *Nat. Mater.* **14** 931–6

Couret A and Caillard D 1989 Glide in prismatic planes in beryllium: an *in situ* study *Strength of Metals and Alloys (ICSMA 8)* (Oxford: Pergamon) pp 211–6

Darot M and Gueguen Y 1981 High-temperature creep of forsterite single crystals *J. Geophys. Res.: Solid Earth* **86** 6219–34
Deer W A, Howie R A and Zussman J 1982 *Rock Forming Minerals*, vol. 1A *Orthosilicates* 2nd edn
(London: Longman)

de Wit R W L and Trampert J 2015 Robust constraints on average radial lower mantle anisotropy and consequences for composition and texture *Earth Planet. Sci. Lett.* **429** 101–9

Durham W B, Goetze C and Blake B 1977 Plastic flow of oriented single crystals of olivine: 2. Observations and interpretations of the dislocation structures *J. Geophys. Res.* **82** 5755–70

Durinck J, Carrez P and Cordier P 2007 Application of the Peierls–Nabarro model to dislocations in forsterite *Eur. J. Mineral.* **19** 631–9

Durinck J, Legris A and Cordier P 2005 Pressure sensitivity of olivine slip systems: First-principle calculations of generalised stacking faults *Phys. Chem. Miner.* **32** 646–54

Evans B and Goetze C 1979 The temperature variation of hardness of olivine and its implication for polycrystalline yield stress *J. Geophys. Res.* **84** 5505–24

Gaboriaud R J, Darot M, Gueguen Y and Woirgard J 1981 Dislocations in olivine indented at low temperatures *Phys. Chem. Miner.* **7** 100–4

Gale J 1997 GULP: a computer program for the symmetry-adapted simulation of solids *J. Chem. Soc. Faraday Trans.* **93** 629–37

Gouriet K, Carrez P and Cordier P 2014 Modelling [100] and [010] screw dislocations in MgSiO₃ perovskite based on the Peierls–Nabarro–Galerkin model *Modelling Simul. Mater. Sci. Eng.* **22** 25020

Gueguen P Y and Darot M 1982 Les dislocations dans la forstérite déformée à haute Température *Phil. Mag.* A **45** 419–42

Hirel P 2015 Atomsk: a tool for manipulating and converting atomic data files *Comput. Phys. Commun.* **197** 212–9

Hirel P, Kraych A, Carrez P and Cordier P 2014 Atomic core structure and mobility of [100](010) and [010](100) dislocations in MgSiO₃ perovskite *Acta Mater.* **79** 117–25

Joós B, Ren Q and Duesbery M S 1994 Peierls–Nabarro model of dislocations in silicon with generalized stacking-fault restoring forces *Phys. Rev. B* **50** 5890–8

Kohlstedt D L, Goetze C, Durham W B and Vander Sande J 1976 New technique for decorating dislocations in Olivine *Science* **191** 1045–6

Mainprice D, Tommasi A, Couvy H, Cordier P and Frost D J 2005 Pressure sensitivity of olivine slip systems and seismic anisotropy of Earth’s upper mantle *Nature* **433** 731–3

Mussi A, Cordier P, Demouchy S and Vanmansart C 2014 Characterization of the glide planes of the [001] screw dislocations in olivine using electron tomography *Phys. Chem. Miner.* **41** 537–45

Panning M and Romanowicz B 2004 Inferences on flow at the base of Earth’s mantle based on seismic anisotropy *Science* **303** 351–3

Peierls R 1940 The size of a dislocation *Proc. Phys. Soc.* **52** 34–7

Phakey P, Dollinger G and Christie J 1972 Transmission electron microscopy of experimentally deformed olivine crystals *Flow And Fracture Of Rocks* 16 (Washington DC: Amer. Geophysical Union) 117–38

Plimpton S 1995 Fast parallel algorithms for short-range molecular dynamics *J. Comput. Phys.* **117** 1–19

Poirier J P 1975 On the slip systems of olivine *J. Geophys. Res.* **80** 4059–61

Price G D and Parker S C 1987 The lattice dynamics of forsterite *Mineralogical Magazine* **51** 157–70

Schoeck G 2005 The Peierls model: progress and limitations *Mater. Sci. Eng. A* **400** 7–17

Sun X Y, Cordier P, Taupin V, Fressengeas C and Jahn S 2016 Continuous description of a grain boundary in forsterite from atomic scale simulations: the role of disclinations *Phil. Mag.* **96** 1757–72

Vítek V 1968 Intrinsic stacking faults in body-centred cubic crystals *Phil. Mag.* **18** 773–86

Walker A M, Gale J, Slater B and Wright K 2005 Atomic scale modelling of the cores of dislocations in complex materials part 2: applications *Phys. Chem. Chem. Phys.* **7** 3235–42

Walker A M, Carrez P and Cordier P 2010 Atomic-scale models of dislocation cores in minerals: progress and prospects *Mineralogical Magazine* **74** 381–413

Wang Y, Liebermann R C and Boland J N 1988 Olivine as an *in situ* piezometer in high pressure apparatus *Phys. Chem. Miner.* **15** 493–7

Wolf D, Keblinski P, Phillpot S R and Eggebrecht J 1999 Exact method for the simulation of Coulombic systems by spherically truncated, pairwise $r^{-1}$ summation *J. Chem. Phys.* **110** 8254

Woodward C 2005 First-principles simulations of dislocation cores *Mater. Sci. Eng. A* **400** 59–67