Записки Российского минералогического общества, 2021, T. 150, № 5, стр. 92-102

INTRODUCTION

Grain boundaries and interfaces are known to strongly influence the mechanical and transport behavior of polycrystalline materials (Sutton, Balluffi, 1995). Many efforts have been made in the last decades (Panchal et al., 2013) to establish the relation between local interfacial structure and chemistry and various phenomena such as microstructure evolution (Mott, 1948; Beck, Sperry, 1950; Burke,Turnbull, 1952), segregation (Lejček et al., 2017), creep (Langdon, 1970; Raj, Ashby, 1971; Watanabe, 1982; Chen, Schuh, 2007), fatigue (Tanaka et al., 1986; Sangid et al., 2011; Sangid, 2013; Musinski, McDowell, 2016), fracture (Watanabe, Tsurekawa, 1999, 2004) and corrosion ( Lin et al., 1995; Lehockey et al., 1997; Shimada et al., 2002). The intrinsic complexity of both homo- and heterophase interfaces requires at least five parameters to uniquely describe their misorientation angle and boundary plane orientation (Sutton, Balluffi, 1995). Thus, various approaches including Coincidence Site Lattice (CSL) and Displacement Shift Complete (DSC) lattice (Sutton, Balluffi, 1995), dichromatic patterns and complexes (based on Shubnikov groups) (Shubnikov, Koptsik, 1972; Pond, Bollmann, 1979; Pond, Vlachavas, 1983), and Bollman’s 0-lattice theory (Bollmann, 1967; Smith, Pond, 1976) were introduced to provide this description. Recent experimental studies show that the chemical and structural diversity of grain boundaries (GBs) can be even greater. According to Dillon et al. (Dillon et al., 2007; Dillon, Harmer, 2007), and Cantwell et al. (Cantwell et al., 2014, 2016), GBs can demonstrate phase-like behavior similarly to bulk crystals and undergo first-order transitions called complexion transitions at the same macroscopic crystallographic parameters. This makes the study of grain boundary structure and its relation to the properties of materials even more difficult.

Experimental observation of grain boundaries and complexion transitions at the atomic scale is extremely complicated since it requires a thorough sample growth and preparation followed by accurate in situ high-resolution microscopy analysis in ultrahigh vacuum (Chen et al., 2008; Cantwell et al., 2020). On the other hand, theoretical simulations of GBs are relatively cheap since they are often carried out within classical molecular dynamics with empirical interatomic potentials, where CSL theory is used for creating the initial orientation of crystallites (Béré et al., 2002; Uberuaga, Bai, 2011; Uberuaga et al., 2013; Frolov et al., 2013, 2015; Fujii et al., 2019). However, this approach cannot provide a general description of all GB phases and complexion transitions for given GB crystallographic parameters, being limited only to manually prepared structural patterns. Furthermore, such calculations can suffer from inaccuracy of both CSL approximation and errors of empirical interatomic potentials. Despite the non-trivial nature of the problem, an effective method of grain boundary structure prediction based on evolutionary algorithm USPEX (Oganov, Glass, 2006; Oganov et al., 2011; Lyakhov et al., 2013), was recently proposed and successfully applied to the study of GB phase behavior in elemental metals (Frolov, Setyawan et al., 2018; Frolov, Zhu et al., 2018; Zhu et al., 2018; Meiners et al., 2020).

EVOLUTIONARY ALGORITHM

Here, we extend this method to grain boundaries in compounds, and boost its efficiency and reliability by using accurate machine learning interatomic potentials. USPEX was previously used to study bulk crystals at zero and finite temperatures (Kruglov et al., 2021), free surfaces (Zhu et al., 2013), two-dimensional materials (Zhou et al., 2016), epitaxial thin films (Mazitov, Oganov, 2021), nanoclusters (Lepeshkin et al., 2018), and grain boundaries in single-component systems (Zhu et al., 2018). The interface system in our method is represented by a thick slab immersed in vacuum and consisting of three parts: two bulk regions and an interfacial region located between them (Fig. 1).

Our algorithm optimizes the structure of the interface keeping bulk regions unchanged. The first generation of individuals is produced with a topological random structure generator (Bushlanov et al., 2019), while all subsequent generations are created by variation operators using the best representatives of the previous generation as parents. There are four variation operators implemented in our method: (1) heredity, (2) softmutation, (3) transmutation, and (4) addition/deletion. The first one creates offspring from two randomly sliced parents by combining their fragments. The remaining operators act on the structure of a single parent by shifting the atoms along with the softest vibrational modes (2), changing the chemical identity of some atoms (3), and adding (or deleting) the atoms according to their coordination number (4). A certain percentage of random structures are included in each generation to ensure diversity of the population. All structures are relaxed and ranked by fitness function based on their energy. This process is repeated iteratively until the set of the fittest individuals remains unchanged for a sufficiently large number of generations.

MACHINE LEARNING INTERATOMIC POTENTIALS

The evolutionary search typically requires hundreds or even thousands of structure relaxa-tions, which should be performed with high accuracy to properly determine the ground state of the system. In contrast to crystal structures, the atomic structures of computationally studied grain boundaries are usually larger by an order of magnitude, which makes their prediction with the usage of ab initio relaxation practically unfeasible. Previous studies of grain boundaries in simple metals (Frolov et al., 2015, 2016; Frolov, Zhu et al., 2018; Meiners et al., 2020; Zhu et al., 2018) show that classical interatomic potentials in principle can be utilized for their simulation and structure prediction. However, the accuracy of these potentials is not always sufficient to reliably determine the ground state, both for bulk crystals and for grain boundaries. In this work, we replace classical potentials with an accurate machine learning (ML) interatomic potential using the MTP package (Shapeev, 2016), and utilize the two-stage relaxation scheme with two separately trained potentials. The first potential was actively trained on relaxation trajectories of random crystal structures generated by USPEX. It is designed to operate in a wide region of the phase space and perform the initial crude relaxation. The second potential was trained on molecular dynamics trajectories of bulk supercells with the addition of point defects and random deformations of the unit cell. Since there were no structures with free surfaces and grain boundaries in such a training set, we extended it with GB structures randomly generated with USPEX and preliminarily relaxed with the first potential. The resulting training set represents a relatively narrow region of the phase space, while the corresponding interatomic potentials turn out to be more accurate and can be used for a final structure relaxation. Data on energies, interatomic forces and stresses of the structures were obtained with density functional theory (DFT) with the generalized gradient approximation (Perdew et al., 1996) and projector-augmented wave method (Bloechl, 1994; Kresse, Joubert, 1999), as implemented in VASP code (Kresse, Furthmuller, 1996). Accuracy of the second interatomic potential in prediction of equilibrium lattice constant ${{a}_{0}}$ and elastic constants ${{C}_{{11}}}$, ${{C}_{{12}}}$, and ${{C}_{{44}}}$ is given in Table 1. For more information on ML interatomic potentials, the reader is referred to work (Podryabinkin, Shapeev, 2017).

RESULTS

We applied our method to study the atomic structure of ${{\Sigma }}5\left( {310} \right)\left[ {001} \right]$ grain boundaries in mineral periclase (MgO). This boundary is a symmetric tilt boundary resulting from the simultaneous rotation of two (310) surface slabs by an angle of 36.9° around the [001] axis. Our choice of this particular GB orientation was conditioned by a presence of a sufficient number of studies in the literature suitable for comparison of the results. Being probably the most stu-died metal oxide, MgO is commonly utilized as a model system for more complex oxides. Grain boundaries are known to significantly affect various properties of polycrystalline metal oxides in a wide range of practically important applications, such as MOSFETs (Kukli et al., 2002; Yanev et al., 2008), fuel cells (Maier, 2000; Suzuki et al., 2001), gas sensors (Kosacki et al., 2005; Dey, 2018), varistors (Clarke, 1999), SQUIDs and high-Tc superconductors (Hilgenkamp, Mannhart, 2002). According to a large number of experimental and theoretical studies (Duffy, 1986; Yan et al., 1998b; Yan et al., 1998a; Parker, Harris, 1999; Harris et al., 2001; Harding, 2003; Wynblatt et al., 2003), point defects in MgO tend to segregate at grain boundaries and diffuse along them, which makes the properties of polycrystalline MgO depend considerably on GB structure. The evolutionary search was carried out for up to 100 generations with 30 individuals in each generation. The initial population of 30 GB structures with 4–16 MgO units (up to 32 atoms in the interface region) was produced by a topological random generator. The structures in subsequent generations were produced by heredity (40%) and softmutation (20%) operators, while the remaining 30% were produced randomly to diversify the population. The thickness of the grain boundary region was automatically adjusted according to the number of MgO units in the structure and the average atomic density of bulk MgO, while both bulk blocks were 50 Å thick. We also considered reconstructions of the grain boundary unit cell with a cell area up to 4 times larger than the original interface unit cell area. Structure relaxations and energy calculations were performed in two stages with the two ML potentials described above, using the LAMMPS package (Plimpton, 1995). Each relaxation stage consisted of a short finite-temperature molecular dynamics run followed by conjugate gradients minimization, where only the atoms of the GB region were allowed to move. Finally, for each GB structure, we calculated the value of the interface energy according to the formula

where $E$ is the energy of whole atomic block with grain boundary, two bulk regions and two free (310) surfaces, ${{E}_{s}}$ is the (310) surface energy, ${{N}_{{{\text{MgO}}}}}$, ${{\mu }_{{{\text{MgO}}}}}$ are the number of MgO units in the structure and chemical potential of MgO, and is the area of the grain boundary.

Results of the evolutionary search are illustrated in Fig. 2. We represented each structure as a point on the phase diagram in ($n$,$\gamma $) space, where $n$ is the atomic density on the grain boundary plane. To derive $n$, we first calculate the number of atoms in the system $N$ and the number of atoms in one (310) plane of the MgO bulk region $N_{{{\text{plane}}}}^{{{\text{bulk}}}}$. Finally, we calculate $n$ as the ratio ${{\left( {N\,{\text{modulo}}\,N_{{{\text{plane}}}}^{{{\text{bulk}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {N\,{\text{modulo}}\,N_{{{\text{plane}}}}^{{{\text{bulk}}}}} \right)} {N_{{{\text{plane}}}}^{{{\text{bulk}}}}}}} \right. \kern-0em} {N_{{{\text{plane}}}}^{{{\text{bulk}}}}}}$. Physical meaning of this quantity is the fraction of atoms from the ideal (310) plane located in the grain boundary plane. This approach, recently proposed by Zhu et al. (Zhu et al., 2018), separates structures with different atomic densities, considering them as grain boundary complexions. Our algorithm successfully found a widely known CSL-type structure (ID 587) with $n = 0$, which was extensively studied in previous works (Fujii et al., 2019; Harris et al., 2001; Uberuaga, Bai, 2011; Yan et al., 1998b) to be the most stable. Predicted value of its interface energy is ${{\gamma }} = 87{\text{\;}}{{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}.{\text{\;}}$ Moreover, in addition to grain boundaries with standard atomic density ($n = 0$), various GB complexions with different atomic densities (ID 1570, ID 434, ID 252) were found. All these structures are essentially modifications of the ground state, which can be observed during grain boundary segregation and diffusion or complexion transitions. Interface energies of these structures, ${{\gamma }} = 102{\text{\;}}{{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}$ (ID 1570), ${{\gamma }} = 105{\text{\;}}{{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}$ (ID 434), and ${{\gamma }} = 102{\text{\;}}{{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}$ (ID 252), are also fairly close to the ground state.

The accuracy of MTP interatomic potential in prediction of the interface energy was subsequently tested on a subset of the best GB structures from the evolutionary search. We selected 25 structures for each value of $n$ with various predicted $\gamma $ including both stable and unstable structures, and calculated their interface energies with VASP, at DFT level of theory. Results of this accuracy test are presented in Fig. 3. MTP demonstrates an outstanding performance in prediction of the interface energy with a root-mean-squared error (RMSE) of 8.4 ${{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}$or 6.7% with respect to mean value of $\gamma $ in this set of structures. Moreover, the predicted interface energy of the ground state (ID 587) differs from DFT value only by 4.8 ${{{\text{meV}}} \mathord{\left/ {\vphantom {{{\text{meV}}} {{{{\AA}}^{2}}}}} \right. \kern-0em} {{{{\AA}}^{2}}}}$ while the atomic configuration of the DFT ground state is almost indistinguishable from that predicted by MTP.

Interface energies in ${{\Sigma }}5\left( {310} \right)\left[ {001} \right]$ GBs are rather ordinary for ionic crystals. In Table 2, we compare our results with typical values of grain boundary formation energies in different substances: copper, tungsten, diamond, tausonite, cubic zirconia and rutile.

It is worth noting that all GB structures found in the calculation have a considerable excess atomic volume. Large excess volume is usually considered proportional to the degree of segregation (Aaron, Bolling, 1972; Frolov, Mishin, 2012a, 2012b). In case of MgO, this tendency can be implicitly confirmed by experimental studies of grain boundary segregation (Yan et al., 1998b; Wang et al., 2011), where the cavities inside the GB plane act as sinks for point defects. Figure 4 shows the distribution of excess volume ${{\upsilon }_{{ex}}}$. For each structure, we calculated the GB excess volume as ${{\upsilon }_{{ex}}} = {1 \mathord{\left/ {\vphantom {1 A}} \right. \kern-0em} A}\left( {\mathcal{V} - \mathcal{N}{{v}_{{{\text{MgO}}}}}} \right)$, where $\mathcal{V}$ is the volume of arbitrary chosen region of the structure containing GB plane, $\mathcal{N}$ is the number of atoms in this volume, and ${{v}_{{{\text{MgO}}}}}$ is the volume per atom in bulk MgO. Defined in this way, ${{\upsilon }_{{ex}}}$ represents the deviation of atomic volume in GB region from the bulk. Phase diagram constructed in (${{\upsilon }_{{ex}}},\gamma $) space allows to study the structural transitions resulting from application of stress ${{\sigma }_{{33}}}$ in the direction normal to the GB plane. Construction of the convex hull in such a diagram gives information about stable phases at each value of ${{\sigma }_{{33}}}$, the stress at which a transition occurs is equal to the negative slope of the corresponding section of the convex hull: ${{\sigma }_{{33}}} = {{ - \partial \gamma } \mathord{\left/ {\vphantom {{ - \partial \gamma } {\partial {{\upsilon }_{{ex}}}}}} \right. \kern-0em} {\partial {{\upsilon }_{{ex}}}}}$.

The values of ${{\sigma }_{{33}}} > 0$ correspond to compression and result in stabilization of dense ID 1421 structure. Negative sign of ${{\sigma }_{{33}}}$ corresponds to tensile stress. The application of increasingly negative stresses induces complexion transitions, where structures with increasing excess volume are stabilized (ID 402, ID 1270) and finally, the system fractures (ID 605). In other words, analysis of stable GB phases in (${{\upsilon }_{{ex}}},\gamma $) space may reveal the fracture mechanics at the atomic scale and provide the thermodynamic description of related structural transitions. In particular, one can estimate the value of tensile strength by calculating the slope of the convex hull section preceding the formation of a crack. The resulting value of ${{\Sigma }}5$(310)[001] grain boundary tensile strength is 1.14 GPa. Experimental observations (albeit for a different grain boundary) of fracture strength in MgO bicrystals (Ku, Johnston, 1964) range from 0.05 to 0.3 GPa, depending on temperature ($T = 300{\kern 1pt} - {\kern 1pt} 1400$ K) and grain size ($D = 1{\kern 1pt} - {\kern 1pt} 6$ mm). The agreement is good, especially given that tensile strength should significantly decrease with temperature, and should differ for different GBs in the same material.

To sum up, here we demonstrate the application of our newly extended method to the prediction of the atomic structure of grain boundaries and interfaces. Utilizing the natural principles of Darwinian evolution, it is capable of automatically exploring grain boundary configuration space, using only the knowledge of the structure, composition, and orientation of adjacent crystallites. The application of machine learning interatomic potentials makes the study of interfacial structure possible for an arbitrary system of interest almost at ab initio level of accuracy. We tested our algorithm on grain boundary complexions in MgO for ${{\Sigma }}5\left( {310} \right)\left[ {001} \right]$ orientation. To perform structure relaxation during the evolutionary search, two machine learning interatomic potentials were trained for two successive relaxation stages. The first one was trained on relaxation trajectories of randomly generated MgO crystal structures, and the second one – on molecular dynamics trajectories of MgO supercells with the addition of point defects and random deformations of the cell. Our results confirm that the well known CSL-type structure of ${{\Sigma }}5\left( {310} \right)\left[ {001} \right]$ GB is the ground state for this orientation, at least in zero-temperature case. In addition to the ground state, several GB complexions with various atomic densities on the boundary plane were found. Moreover, the analysis of grain boundary excess volume revealed the atomistic mechanism of fracture under applied stress. The proposed methodology allows one to thoroughly investigate structures of any crystalline interfaces, which may significantly deepen our knowledge and understanding of this type of systems.