Complexions and stoichiometry of the 60.8°//[100](011) symmetrical tilt grain boundary in Mg₂SiO₄ forsterite: A combined empirical potential and first-principles study

Grain boundaries (GBs) play a key role in numerous physical processes associated with mechanical, chemical diffusion, or electrical conductivity behavior of polycrystalline materials (Balluffi and Sutton 1996). Understanding their properties is, therefore, particularly important for a mineral like olivine, which is the principal constituent of the rocks of the Earth’s upper mantle. GBs are involved in several first-order processes such as creep (Hirth and Kohlstedt 1995; Maruyama and Hiraga 2017a, 2017b), grain boundary migration (Bollinger et al. 2018; Furstoss et al. 2020) and diffusion (Demouchy 2010; Fei et al. 2016). At the mesoscopic scale, GBs are often viewed as interfaces with particular properties, like their mobility or their diffusion coefficient. Unfortunately, their structure at the atomic level is challenging to characterize. At the moment, atomically resolved observations of GB in olivine are scarce (Marquardt and Faul 2018), and the few high-resolution transmission electron microscope (HRTEM) micrographs available in literature (Marquardt and Faul 2018; Fei et al. 2016; Heinemann et al. 2005) do not allow determination of the chemical composition of GB, i.e., whether or not the GB is stoichiometric.

As a complement to experimental efforts, simulations at the atomic scale can provide insights into the structure of GBs. Yet even that is challenging due to the complex crystallography of olivine and the many possible GB configurations to explore. Nonetheless, people have used density functional theory (DFT) (Ghosh and Karki 2014) or classical molecular dynamics (MD) (Adjaoud et al. 2012; Mantisi et al. 2017) to determine the physical properties (energy, excess volume) of some select GBs in olivine. Ultimately, the comparison of numerical models with high-resolution observations will be necessary to assess the relevant GB structures and infer their physical properties.

Description of a GB can be done at multiple scales involving different degrees of complexity and precision. At the mesoscopic scale, a GB is generally described by the disorientation between the two adjacent grains, and its energy is almost always described through this variable using extended Read-Shockley type models (initially designed only for low-angle GBs) (Read and Shockley 1950; Gui-Jin and Vitek 1986). However, a more complete “macroscopic” description should also specify the rotation axis and the crystallographic planes in contact at the GB. At smaller scales, this description can be enriched by indicating the translation vector between the neighboring crystals and finally by describing the atomic arrangement of the GB (e.g., bonds, stoichiometry, charge). Using these lower scale descriptions introduces additional degrees of freedom, which also impact the physical properties of GB (Han et al. 2016) and question the common assumption that a given disorientation is related to a single GB structure. The multiplicity of GB structures even for a single GB disorientation and contact plane has been considered by several studies on metals (Rittner and Seidman 1996; Oh and Vitek 1986), but rarely for geologically relevant materials (Hirel et al. 2019). Recently the concept of GB complexion was introduced (Cantwell et al. 2014) to account for the multiplicity of GB atomic structures. In spite of growing evidence that GBs in a polycrystal may exist in a wide variety of complexions, numerical studies often only account for the complexion of lowest energy, discarding a great number of other possible complexions.

In the present work, we use a combination of molecular statics and ab initio calculations to study the low-energy complexions of a particular GB in the magnesium member of olivine, forsterite, or Mg₂SiO₄. Forsterite is an orthorhombic crystal in which magnesium ions are in two distinct sites while oxygen ions (three distinct sites) are arranged in tetrahedra with a silicon ion at the center. We use the Pbnm space group (where the cell parameters are ordered, such as b > c > a) to describe the forsterite crystal. The studied GB corresponds to a symmetrical tilt grain boundary (STGB) with a disorientation of 60.8° about the [100] direction and (011) planes in contact at GB. This orientation corresponds to a high symmetry of the quasi hexagonal close-packed (hcp) oxygen sub-lattice. In such a GB, the hcp stacking of oxygen ions is almost identical to the bulk crystal, which might explain its over-representation in natural olivine aggregates. Indeed, Marquardt et al. (2015) have shown from EBSD measurements on annealed forsterite polycrystals that among GB with 60.8° disorientation, already highly represented in the misorientation distribution function, the rotation about the [100] axis with a (011) GB plane is the more ubiquitous. Although first principles (Ghosh and Karki 2014) and MD (Adjaoud et al. 2012) works have examined this special STGB, both have focused on a single complexion, which raises the question of its relevance and occurrence in natural systems. Most recent studies of GB in metallic systems have revealed that many metastable configurations are possible and should be accounted for, arguing that materials are rarely in their lowest energy state (Han et al. 2016, 2017). In this work, we follow this hypothesis to explore multiple possible complexions of the GB in forsterite.

Here, we first review the low-energy complexions of this STGB by determining their energies, excess volumes and by describing their atomic-scale features. Classical MD simulations allow us to probe the energy landscape and identify low-energy configurations. Then, we perform ab initio calculations to characterize more accurately the properties of these configurations.

Deviation from stoichiometry at GB has been noticed by several experimental and numerical studies (Farkas 2000; Baker et al. 1990). In fact, GBs are known to be sources and wells of vacancies and, therefore, it seems interesting to explore the stoichiometry as another degree of freedom for complexion by deviations from the bulk composition at GB. In ionic compounds such as forsterite, an easy way to deviate from stoichiometry, while still keeping the system electrically neutral, is to incorporate neutral vacancy pairs. In forsterite, the energetically most favorable Schottky defect is the MgO one (Brodholt 1997). Thus, we investigate, in a second step, the effect of stoichiometry deviation at GB by inserting MgO vacancy pairs.

Atomistic calculations are performed with LAMMPS (Plimpton 1995). Interactions between ions are modeled with a rigid-ion potential accounting for the long-range Coulomb interactions, and short-range interactions are described by a Morse function and a repulsive r⁻¹² term. Potential parameters were optimized by Pedone et al. (2006), where ions have partial charges of 1.2e, 2.4e, and −1.2e for Mg, Si, and O ions, respectively. This allows considering neutral vacancy pairs such as MgO. Moreover, this potential was shown to offer a very good description of several properties of forsterite, including bulk, surface, and defect properties (Hirel et al. 2021). Coulomb interactions are computed using the particle-particle particle-mesh (pppm) method (Eastwood et al. 1980). Ion positions are optimized by means of the conjugate-gradient algorithm until the maximum force is smaller than 10⁻¹⁰ eV/Å. The pressure is maintained at 0 GPa by rescaling the simulation cell.

To consolidate our results and to investigate possible changes in ionic charge state, we also performed first-principles calculations on some systems containing grain boundaries. We use the Vienna Ab initio Simulation Program (VASP) (Hafner 2008) with the projector augmented wave (PAW) method (Kresse and Joubert 1999) in the generalized gradient approximation (GGA) (Perdew et al. 1996). To treat explicitly only valence electrons, the interactions between valence and core electrons with nuclei are described by the pseudopotential developed by Perdew and Wang (1992). A cut-off of 500 eV is applied to the plane wave basis and the first Brillouin zone is sampled with Monkhorst-pack grid (Monkhorst and Pack 1976) using a 2 × 1 × 1 mesh.

First-principles calculations give access to the electron density. We perform a Bader analysis, which decomposes the charge density into atomic basins based on a change of sign in the electron density gradient (Bader 1990). The calculations are performed using the software “Bader Charge Analysis” developed by the Henkelman group (Henkelman et al. 2006; Tang et al. 2009). We use it to quantify the changes in ion charges in GBs with respect to a perfect crystal environment.

We focus on the 60.8°//[100](011) symmetric tilt grain boundary (STGB) in forsterite, i.e., where two crystals share a common [100] axis and meet along (011) planes with a disorientation of 60.8°. In forsterite, two different types of (011) planes, with different stoichiometries, are possible (Fig. 1). One type of plane is terminated by a silicon ion and the oxygen ions forming the edge of tetrahedra, which we will refer to as edge-planes, and the second type is terminated by Mg ions and an oxygen ion belonging to the tip of a tetrahedron, referred to as tip-planes. We begin by constructing cells of forsterite terminated by the same type of plane (Fig. 1). Such cells preserve the stoichiometry but must be truncated by ½[001] (which removes one formula unit) to obtain the same type of surface.

The GB is constructed by stacking two crystals meeting across (011) planes. To avoid any spurious effect due to charged surfaces, we use 3D periodic boundary conditions (PBC) and ensure that the system remains charge-neutral. Moreover, we also make sure that SiO₄ tetrahedra remain unaffected because breaking Si-O bonds is energetically unfavorable in this system. Since the GB studied is a STGB, the two crystals are mirror images of each other. Two types of reference GBs are constructed, one by stacking the “edge-plane” system with its mirror image (“edge-to-edge”), and one by stacking the “tip-plane” system with its mirror image (“tip-to-tip”).

Due to PBC, the same GB is formed at the center of the cell and at the edges (Fig. 2). This allows computing the reference GB energy density (γ_ref in J/m²) and excess volume (Ω_ref):

where E_tot and V_tot are the total energy and volume of the relaxed system containing the GBs, E_bulk and V_bulk are the total energy and volume of the equivalent defect-free crystal, and A is the area of the GB.

We computed the reference GB energy for different system sizes and found that a height of 180 Å is necessary for GB energy and excess volume to converge, which corresponds to a system containing 980 atoms. When performing DFT calculations, we use smaller systems of 420 atoms due to computational costs. Although GB energies do not converge for this system size, our convergence tests show that the atomic structure is the same as in larger systems.

The method described above does not guarantee an optimal atomic GB structure. An important degree of freedom is the relative translation of the two grains, quantified by a vector $τ→$ = (τ_x,τ_y,τ_z). To explore this degree of freedom, we impose the translation vector (τ_x,τ_y) in the GB plane and relax ions in the direction normal to the GB to determine the optimal translation τ_z. Computing the GB energy for each translation vector, one obtains a map of the energy density γ along the GB plane. Similar to the classical generalized stacking fault approaches (Vitek 1968; Mishin and Farkas 1998), this method allows the identification of low-energy configurations of the GB.

The use of 3D PBC requires special attention to the simulation cell construction and the relaxation scheme involved during the relaxation process (Fig. 2). Indeed, if one crystal was displaced with respect to the other by keeping the same cell, then the GBs at the center and at the edges of the cell would not be equivalent anymore, thus preventing a direct calculation of the interface energy. For this reason, when the top crystal is translated, we tilt the box so that the interface at the edges of the cell remains identical, with the energy E_ref. For the given relative translation $τ→$⁠, the energy and excess volume of the central GB are then:

For a given imposed translation vector (τ_x,τ_y), we relax atom positions along the direction z normal to the GB. To allow for rotational relaxation of tetrahedra, O ions are relaxed in all directions. Within 25 Å in the vicinity of the central GB, all ions are left free to relax (see Fig. 2).

After the identification of the lowest energy configurations (Fig. 3), the selected systems are fully relaxed (no ionic positions fixed), which provides the final GB complexions and allows computing their energies and excess volumes.

From the atomistic simulation outputs, we compute the density profile along a direction normal to the GB plane. In practice, the mass of each atom is expanded through Gaussian distributions (with a standard deviation of 2 Å) and the density profile ρ(z) is computed from the contribution of all masses contained in an elementary volume (grid size of 0.2 Å). These calculations are performed with Atomsk (Hirel 2015).

The GB excess volumes can be expressed as a function of the density profile along the direction normal to the GB plane. Indeed, the density can be related to an infinitesimal change in volume dV(z):

where ρ_UC is the density of a unit cell. The GB excess volume is then computed by integrating dV around the GB in a length F, and dividing by the GB area:

We verified that this methodology yields the same values as the ones given by Equations 1 and 2 in the case of stoichiometric GB structures. For nonstoichiometric GBs, Equations 1 and 2 are no longer valid and Ω is computed using Equation 4.

As explained in the previous section, the first step toward the calculation of the GB energy landscape is the construction of a reference configuration containing two equivalent GBs. After full relaxation, the edge-to-edge GB (labeled E1) has an energy γ_E1 = 0.9 J/m² and an excess volume Ω_E1 = 0.39 Å. Its atomic configuration is shown in Figure 4a. Atomic displacements after relaxation are small (about 1.15 Å), and tetrahedra do not rotate with respect to the [100] direction. Computation of the mass density profile shows that the GB is less dense than the surrounding crystals, which is consistent with the positive excess volume.

The tip-to-tip GB (labeled T1) has an energy γ_T1 = 1.22 J/m² and an excess volume Ω_T1 = 0.76 Å. Displacements are also moderate (about 1.5 Å) but are associated with significant tetrahedra rotation, as shown in Figure 4d. This is probably due to strong repulsion of O ions at the tip of tetrahedra.

Following the methodology presented in “Conservative sampling” section, we computed the GB energy landscape for the two types of GB, which are reported in Figure 3. These landscapes are irregular and present abrupt energy variations, as it was already noticed by Adjaoud et al. (2012). Nonetheless, both surfaces exhibit a central symmetry, which is consistent with the symmetries of the (011) plane in forsterite. In addition to the reference GBs E1 and T1, each energy landscape reveals two more low-energy basins, labeled E2, E3, T2, and T3. By symmetry, equivalent configurations are found, respectively E2′ = E2, E3′ = E3, T2′ = T2, and T3′ = T3 (see Fig. 3).

Each low-energy configuration found from the energy landscape sampling is fully relaxed (i.e., without constraint). We obtain six distinct GB complexions with energies ranging between 0.66 and 1.25 J/m². Reference GB complexions E1 and T1 were already presented above.

The complexion E2 is shown in Figure 4b. Contrary to the reference E1, it shows significant relaxation with tetrahedra sharing an O ion (orange arrow in Fig. 4b), which causes an O ion to become free of any bond with a Si ion (blue arrow in Fig. 4b). Although the presence of this free O ion seems unfavorable, the interface energy is γ_E2 = 0.66 J/m², significantly lower than the reference E1. This complexion is also more compact, as evidenced by the smaller excess volume (i.e., 0.13 Å) and the density profile.

The complexion E3 is shown in Figure 4c. Like E2, it contains joined tetrahedra and free O ions, however, its energy is the highest with a value γ_E3 = 1.25 J/m². This complexion is the least favorable for the stoichiometric GBs. It is also associated with a large expansion (Ω_E3 = 0.60 Å).

Concerning the tip-to-tip configurations, the complexion T2 (Fig. 4e) also displays joined tetrahedra and unbound O ions. Its energy γ_T2 = 1.08 J/m² is lower than the reference T1. Finally, the configuration T3 also contains joined tetrahedra and free O ions. Its energy is the lowest for the tip-to-tip GB type, γ_T3 = 0.95 J/m². It is also the most compact of all complexions, with an excess volume Ω_T3 = 0.13 Å (see Table 1).

In summary, the most favorable complexions are E2 for the edge-to-edge type and T3 for the tip-to-tip type. Both are characterized by tetrahedra connected by an O ion and free O ions remaining in the GB. We also note a positive correlation between GB energy and excess volume: the lower the GB energy, the lower its volume. This is consistent with previous studies on similar GBs (Adjaoud et al. 2012; Ghosh and Karki 2014). From the density profiles, we estimate the structural width of GBs to be ~10 Å.

To complement the results obtained with the interatomic potential, we performed first-principles calculations. Starting from our previous relaxed configurations, we performed a full ionic relaxation by means of DFT calculations.

Ionic displacements are negligible, so the GB atomic structures remain the same for all complexions. The GB energies and excess volumes obtained by DFT are reported alongside those presented previously in Table 1. We obtain very good agreement between the two methods, thus confirming the suitability of the potential for modeling GBs. The ordering of GB energies is identical in both methods, except for the complexions E1 and E2. This small discrepancy may be attributed to the small size of simulation volumes involved in the DFT calculations.

Using the electron density from our DFT calculations, we performed Bader analysis to evaluate the charges of ions. Figure 5b shows the distribution of the charge of O ions in a perfect bulk environment (blue curve), with a median value of about −1.61e, which we consider as the reference charge here. Charge analysis of oxygen ions in GBs shows that they have little effect on this charge distribution. As an example, we show the charges of O ions in the complexion E2 in Figure 5b (orange curve). Their distribution is slightly wider than in the bulk but remains within 1% of the reference charge. Two noticeable exceptions are indicated with arrows. The first one (orange S arrow) corresponds to the O ion bonding two tetrahedra, which has a reduced electron charge. Visualization of its isosurface of charge density (Fig. 5a) shows that it is elongated in the direction of the two Si ions, which indicates mixed ionic and covalent bonding. The second exception is the unbound O ion (blue F arrow in Fig. 5), which has a greater charge of about −1.65e. Its isosurface is almost spherical, indicating that this ion is free of covalent bonding. Overall, we find that using an interatomic potential with fixed charges is a reasonable approximation for modeling GBs in forsterite.

We investigate the effect of stoichiometry by introducing MgO vacancy pairs in the six GBs presented above. We construct atomic systems corresponding to all possible combinations of Mg and O vacancies, which represent about 800 different systems. After relaxation, we obtain several GB complexions that can be compared with the parent GB complexion. Contrary to the stoichiometric case, the interface energy cannot be computed unambiguously, which will be discussed below. For that reason, in this section, we focus on the segregation energies of the vacancies, ΔE, defined as the difference between the energy of a system where a vacancy pair is inside a grain $ErefN−2$⁠, and one where it lies inside the GB $EtotN−2$⁠:

This quantity indicates if it is more energetically favorable for the vacancy to be in the bulk phase (i.e., ΔE > 0) or within the GB (i.e., ΔE < 0). Since in forsterite, two Mg sites and three O sites are possible in the perfect crystal, $ErefN−2$ can take multiple values. We choose to take the one that minimizes ΔE, so that it vanishes in the bulk.

Figure 6 presents the segregation energy ΔE as a function of the location of the MgO vacancy pair for the complexion E1. We can see that it takes negative values everywhere within the GB, which corresponds to the preferential segregation of vacancies at the GB. We conclude that when such vacancies are present in the system, it is energetically favorable for them to be incorporated in the GB rather than inside a grain. Similar results are obtained for all six complexions and are available in the Online Materials¹.

For each initial GB complexion, among all the defective configurations investigated, we select the one with the lowest energy. This leads to six final, non-stoichiometric GB complexions. These complexions are associated with segregation energies in the order of −2 eV, as reported in Table 1.

Their atomic configurations are presented in Figure 7. In the reference complexions E1 and T1, the removal of a MgO pair causes one tetrahedron to miss an O ion and to share one with a neighboring tetrahedron.

The stoichiometric complexions E2, E3, T2, and T3 already contained connected tetrahedra and free O ions, as presented above. One would expect the most favorable configuration to be the one where the free O ion would be removed. We find that this is true only in the complexion T3n, as shown in Figure 4f. On the contrary, for complexions E2n, E3n, and T2n, we find that the most favorable systems are the ones where the free O ions are untouched, and three tetrahedra are connected, as shown in Figures 4b, 4c, and 4e.

Although vacancies were introduced in the GBs, we find that they have a negligible impact on the formation volume, as evidenced by the density profiles shown in Figure 7. As a result, we expect that a simple measurement of the density would not allow discrimination between the different GB complexions. Here again, we estimate the GB structural widths to be ~10 Å.

As done for the stoichiometric complexions presented previously, we performed DFT calculations on the six nonstoichiometric GBs. Starting from the configurations relaxed with the interatomic potential, we perform a full ionic relaxation with DFT. Again, we find that this relaxation does not change the atomic structure of GBs, indicating that the interatomic potential is efficient for determining GB configurations. In particular, DFT confirms the stability of tetrahedra sharing O ions, as well as the excess volumes (see Online Materials¹).

We apply Bader analysis to the non-stoichiometric GBs. The results are presented for the non-stoichiometric complexion E2 in Figure 8. As seen previously, the shape of isosurfaces allows us to recognize mixed ionic and covalent Si-O bonds in tetrahedra, as well as free O ions. Where tetrahedra are connected, shared O ions (arrows labeled S) have a decreased charge, while other O ions belonging to those tetrahedra (arrow labeled T) have a greater charge. As before, free O ions (arrows labeled F) have a higher charge.

To the best of our knowledge, we have identified for the first time six stoichiometric complexions of the 60.8°//[001](011) GB in forsterite, with formation energies ranging from 0.66 to 1.25 J/m². For the same GB, Adjaoud et al. (2012) reported energy of 1.30 J/m² and an excess volume of 0.35 Å using an interatomic potential. First-principles calculations (Ghosh and Karki 2014) reported energy of 1.15 J/m² and an excess volume of 0.37 Å. We find that our own values are in good agreement with both previous studies. Both stoichiometric and non-stoichiometric GBs that we modeled exhibit structural widths of ~10 Å, a value commonly observed by transmission electron microscopy for GBs in olivine as reviewed by Marquardt and Faul (2018).

Experimental determination of GB energies is scarcer. Duyster and Stöckhert (2001) measured GB energies from dihedral angles in an equilibrated coarse-grained natural peridotite. Their values span between 1.12 and 1.47 J/m², which is in good agreement with our results. The values reported by Cooper and Kohlstedt (1982) are significantly lower, 0.9 ± 0.35 J/m², probably due to the presence of melt at the interfaces in the olivine-basalt system considered.

The linking of tetrahedral units and the free O ions at forsterite GBs, as observed in this work, has never been noticed by prior studies. The finding of Mantisi et al. (2017) that GBs are almost structureless above 3 Å cannot be confirmed here. Although differing in their structures, the complexions presented here exhibit a short-range order that can be seen in Figures 4 and 7. Some structural characteristics highlighted by Ghosh and Karki (2014) are seen in our simulations, such as the conservation of Si ions in a tetrahedral environment and modification of Si-O bond lengths. A continuous treatment of this particular GB by Sun et al. (2016) shows that the O and Si sub-lattices accommodate the main part of the crystal disorientation, which is in good agreement with the structure of the complexions presented here.

The combination of first principles and MD calculations has permitted us to improve our confidence in the interatomic potential of Pedone et al. (2006) to the study of GBs in forsterite. All stable complexions found with the potential have been confirmed by ab initio calculations with good agreement concerning the GB atomic structures, energies, and excess volumes. Charge analysis has shown that the use of fixed partial charges (the so-called rigid ion potential) is justified. The small deviations from the reference charge also justify the use of the interatomic potential.

The low-energy basins in the energy landscapes are separated by high-energy barriers. It implies that GBs cannot change complexion conservatively (e.g., E1 → E2) as it would require too much energy. Adsorption of vacancies at the GB may provide a more favorable path toward a change in GB complexion. We have shown that it is more favorable for MgO vacancy pairs to be within the GB than in the bulk phase. This finding is in good agreement with the expected vacancy sink behavior of GBs (Uberuaga et al. 2015).

As explained above, the intrinsic GB energy cannot be defined unambiguously when it is non-stoichiometric. One must account for the chemical potential of the removed ions:

The chemical potential depends on environmental conditions. One possible choice is to define it as the lattice energy of MgO periclase, which is μ_MgO = –16.57 eV using the same interatomic potential by Pedone. Nevertheless, other approximations can be performed, for example, μ_MgO = μ_Mg2SiO4 – μMgSiO₃ (where μ_Mg2SiO4 and μMgSiO₃ are the lattice energies of forsterite and enstatite, respectively), which gives μ_MgO = –17.00 eV; or 2μ_MgO = μ_Mg2SiO4 – μ_SiO2 (where μ_SiO2 is the lattice energy of quartz), which gives μ_MgO = –16.44 eV. Using these values, the GB energy of the E3n complexions (which has the lowest segregation energy) should range between 1.60 and 1.76 J/m₂, which is higher than the one of the parent GB (i.e., 1.25 J/m²). More generally, the energies of the non-stoichiometric complexions are higher than that of the parent GB.

A comparable numerical study (Chua et al. 2010) has also noticed that for different GBs in SrTiO₃ the most stable GBs (i.e., with lowest energies) were the stoichiometric ones, which is in contradiction with the non-stoichiometric GBs observed in SrTiO₃ (McGibbon et al. 1994; Kim et al. 2001; Yang et al. 2013).

The GB can be viewed as a phase at equilibrium with other phases (e.g., neighboring grains) (Cantwell et al. 2014), and the minimization of the whole system energy does not necessarily lead to the minimization of the GB energy itself. Our calculations show that if vacancies are present in the system, they will segregate at GBs even if it increases the GB energy. The GB structure and its energy, therefore, depend on external factors, which emphasizes the lower importance of the GB energy on the stable GB structures.

The atomic structure of GBs is a key input parameter, as all effective physical properties derive from it. For the one disorientation that we studied here (i.e., 60.8°), we expect each GB complexion to have different properties in terms of mobility, diffusion, segregation, and so forth. We want to stress the importance of considering all possible complexions of GBs when performing atomic-scale simulations instead of focusing on a single particular complexion.

Our work also demonstrates the attractiveness of GBs for vacancies (Tschopp et al. 2012; Uberuaga et al. 2015), highlighting the exchange capacities between a crystal and GBs. At finite temperatures, these exchanges should be of first importance for the structure of GBs.

A remarkable observation from our study is the presence of an anomalous non-silicate oxygen site in some GBs, which may have several implications. First, these free O ions may diffuse easier than the others (which are bound in a tetrahedral environment), which could contribute to the high diffusivity of oxygen in GBs compared to bulk (Yurimoto et al. 1992). It may therefore increase the O vacancy concentration in forsterite GBs.

Second, the occurrence of such an unbound O ion is likely to have implications on water storage. The same feature is found in the structure of wadsleyite where it has been shown to strongly favor protonation (Smyth 1994; Jacobsen et al. 2005). Our calculations suggest the possibility of a strong segregation of hydrogen at GBs with a maximum concentration of 1.74.10⁻² Å⁻² on the basis that each unbound O ion is protonated. Assuming spherical grains, the contribution of such a mechanism to hydrogen storage should be of the order of 390, 40, and 4 H/10⁶ Si for grain sizes of 1, 10, and 100 μm, respectively.

Finally, some stable GBs can have relatively high free volumes, and we can expect these sparsely dense boundaries to be preferential zones for incompatible elements. Further simulations are required to properly assess these hypotheses.

This project has received funding from the French government through the Programme Investissement d’Avenir (I-SITE ULNE/ANR-16-IDEX-0004 ULNE) managed by the Agence Nationale de la Recherche, under the project name NuMoGO, and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program under grant agreement no. 787198—TimeMan. Computational resources have been provided by the DSI at Université de Lille.