Ab initio and atomistic study of generalized stacking fault energies in Mg and Mg–Y alloys

Magnesium–yttrium alloys show significantly improved room temperature ductility when compared with pure Mg. We study this interesting phenomenon theoretically at the atomic scale employing quantum-mechanical (so-called ab initio) and atomistic modeling methods. Specifically, we have calculated generalized stacking fault energies for five slip systems in both elemental magnesium (Mg) and Mg–Y alloys using (i) density functional theory and (ii) a set of embedded-atom-method (EAM) potentials. These calculations predict that the addition of yttrium results in a reduction in the unstable stacking fault energy of basal slip systems. Specifically in the case of an I2 stacking fault, the predicted reduction of the stacking fault energy due to Y atoms was verified by experimental measurements. We find a similar reduction for the stable stacking fault energy of the non-basal slip system. On the other hand, other energies along this particular γ-surface profile increase with the addition of Y. In parallel to our quantum-mechanical calculations, we have also developed a new EAM Mg–Y potential and thoroughly tested its performance. The comparison of quantum-mechanical and atomistic results indicates that the new potential is suitable for future large-scale atomistic simulations.

hand, other energies along this particular γ -surface profile increase with the addition of Y. In parallel to our quantum-mechanical calculations, we have also developed a new EAM Mg-Y potential and thoroughly tested its performance. The comparison of quantum-mechanical and atomistic results indicates that the new potential is suitable for future large-scale atomistic simulations.

Introduction
Magnesium alloys are promising candidate materials for lightweight industrial applications [1,2]. Unfortunately, a wider use of magnesium and many Mg-based alloys is impeded by their poor room-temperature ductility, which is caused by basal-type texturing and limited available deformation mechanisms. It is known [3] that addition of Y and/or rare-earth elements increases the ductility by texture weakening and higher activity of non-basal slip systems. Specifically, shear deformations with a component along the 0001 axis, including twinning and c + a slip systems, were shown (e.g. in [4]) to determine the mechanical behavior of hcp metals (figure 1). The limitation to mainly basal a slip and tensile twinning in pure Mg and most commercial wrought Mg alloys causes strain localization in a few shear bands and early fracture of the materials along these shear bands. It has been shown experimentally [5] that alloying Y leads to higher activation of non-basal deformation modes, i.e. non-basal dislocation slip and compression twinning, and consequently to much enhanced ductility at room temperature. Nucleation and mobility energies of dislocations, both perfect and dissociated dislocations, are related to the stable and unstable stacking fault energies which are available via atomistic calculations.
Different slip mechanisms in Mg have been intensively studied both experimentally [5][6][7][8] and theoretically [4,[9][10][11]. Theoretical studies generally focused on calculating so-called generalized stacking fault energy (GSFE) profiles obtained by sliding one half of a crystal over the second half across a given glide plane. As shown by Vitek [12], these energy profiles are closely connected with the motion of dislocations and the plasticity of materials in general. GSFEs of important planes along different slip directions were calculated employing both density functional theory (DFT) and the embedded atom method (EAM) for pure Mg [13][14][15][16][17]. It was shown, for example, that maxima occurring along GSFE profiles, so-called unstable stacking fault energies (USFEs), determine the behavior of slip modes and thus the ductility in fcc metals [18].
Our study is motivated by the fact that rather little is known about GSFEs in Mg-Y alloys, especially in the case of non-basal deformation modes (see e.g. table 1). Therefore, we study and compare the influence of Y additions on the behavior of both basal and non-basal slip systems. Our aim was to provide results that will contribute to the understanding of enhanced room temperature ductility of Mg-Y alloys [5] and to the design of ductile Mg-based materials in the future.
Due to the fact that it is experimentally very difficult to obtain a deeper insight into various mechanisms that are simultaneously active in Mg alloys at atomic scale, we employ theoretical methods as these allow for decomposing this complexity into individual processes and studying them separately. Unfortunately, even theoretical modeling is limited by available computer resources, especially in the case of state-of-the-art quantum-mechanical calculations. The predictive power and reliability of these calculations can nowadays be applied only to problems involving a few hundreds of atoms. As processes related to the plasticity of materials typically span over much larger scales, methods based on inter-atomic potentials are often used instead. For these simulations, well-tested and reliable potentials are necessary. Therefore, we have not only employed quantum-mechanical calculations but, in a complementary manner, we have also developed a new Mg-Y EAM potential. We first carefully validated our potential and only after this then used it to simulate processes involving so many atoms that quantum-mechanical calculations would be prohibitively long. In order to verify our theoretical predictions, an experimental TEM analysis of the I 2 stacking fault was performed and confirmed our theoretical findings.

γ -surface simulations
Employing theoretical methods, slip deformation modes can be studied by calculating generalized stacking fault energies, known as γ -surfaces [12]. Key properties of slip deformation modes such as the Peierls barrier as well as stacking fault energies can be deduced from the γ -surfaces. We determine the generalized stacking fault energies by incrementally shifting the upper half crystal along the slip direction and calculating the total energy of the system as a function of the applied shift vector. In order to facilitate the consideration of a single stacking fault, and so allow us to reduce the required system dimensions, we added the applied shift vector to the lattice vector along the glide plane normal (in our case the glide plane normal is denoted as the z-direction). For the minimization after each incremental shift, atomic positions are constrained along the lateral dimensions of the glide plane, but atoms can reduce the total energy of the system by relaxing in the direction perpendicular to the studied glide plane.
In order to eliminate any spurious interactions between the interface (i.e. the glide plane) and its periodic images, the supercell size must be sufficiently large. We therefore converged the GSFE with respect to the number of atomic planes parallel to the fault plane. For example, in order to determine an optimum number of atomic planes for the {0001} 1120 slip system, supercells with 2, 4, 6, 8 and 12 layers along the 1120 direction, the z-direction (parallel to the plane normal), were used to calculate GSFEs. The GSFE difference between 6-layer (2 × 2 × 6, 24 atoms) and 12-layer (2 × 2 × 12, 48 atoms) supercells is less than 1%. Therefore, supercells with more than six layers are considered large enough. In our case, a 48-atom supercell was employed to obtain a realistically low Y concentration. Converged supercell sizes of all the other slip systems considered were determined similarly (see table 2). In calculations of Mg-Y alloys, one Mg atom is replaced by a Y atom on the first plane (with the smallest value of the , where 0 x 1 and 0 y 1. More details about the calculation of the γ -surface can be found in [19].

Parameters of ab initio calculations
In our work, we calculate GSFE profiles using both DFT [20,21] and the EAM [22] (using the LAMMPS [23] package.). Our DFT calculations were carried out using the projector augmented wave method [39] and the electronic exchange-correlation effects were described by the generalized gradient approximation [40] as implemented in the Vienna Ab initio Simulation Package (VASP) [41,42]. A plane wave energy cut-off of 350 eV is used for both pure Mg and the Mg-Y alloys. The Brillouin zone was sampled using dense Monkhorst-Pack [43] k-point meshes that were chosen to ensure a convergence of the total energy to within 1 meV per atom. Atomic relaxations were performed until the energy and forces converged to 10 −7 eV and 10 −3 eV Å −1 , respectively. To test our computational parameters, we have calculated the lattice parameters of pure Mg. The results, a = 3.1886 Å and c/a = 1.6261 (table 3), are in good agreement with experimental data suggesting a = 3.21 Å and c/a = 1.624 [44].

Development of Mg-Y embedded atom method potential
Empirical potentials provide a means of exploring the physics of systems of atoms on lengthand time-scales currently inaccessible to more computationally expensive ab initio techniques. Such potentials seek to approximate the energy of a system of atoms, and the forces between those atoms, as a classical function of the relative positions of the atoms treated as point particles. The key attributes of such potentials are: (i) the level of physics they aim to capture, i.e. whether they are simple pair potentials or attempt to capture many-atom effects and whether they are functions only of atomic separations or also of bond angles; (ii) the nature of the functions used in the fitted potential, which is to say, the extent to which those functions have a (perhaps physically motivated) prescribed form versus being, say, piecewise spline fitted; and (iii) the data to which those functions are fitted, which may be energies and elastic properties from experiment and ab initio calculations or individual atomic forces from ab initio calculations for a set of atomic configurations. The key issue for empirical potentials, especially when they are to be used to explore systems beyond the reach of ab initio, is their transferability: it is one thing to correctly reproduce the properties to which a potential is fitted and quite another to provide accurate predictions of 'unseen' data. As no existing potential (e.g. [45]) has been developed specifically for GSFE calculations, we provide a new one.
In this study, an EAM interatomic potential was developed to describe the Mg-Y system, which has the following formalism [22,24]: and φ(r ), ρ(r ) and F(n i ) are the pair, density and embedding functions. In this work, the three functions φ(r ), ρ(r ) and F(n i ) are represented with quintic spine interpolations [25]. We used 15 equidistant spline knots for both the pair and the density functions (in the fitting range of 0.5-6.5 Å), and 6 spline knots for the embedding function. For the Mg-Y binary system, a total number of 83 parameters was fitted.
The EAM potential was developed based on the force-matching method [25,26] as implemented in the POTFIT package [27]. To this end, a first-principles database was first established to provide a coarse-grained potential energy surface (PES) of the Mg-Y system. The database was constructed to encompass a wide range of atomic configurations, including all crystalline phases in this alloy system, as well as their derivative structures such as defects, crystal equations of state, deformation paths, melting and cooling trajectories, etc. In addition to six crystal structures of Mg and Y (i.e. hcp, bcc, fcc, 9R, diamond and sc structures), three crystallographic types of Mg-Y intermetallic compounds [28] are considered in this work, including Mg 24 Y 5 (I43m), Mg 2 Y (P6 3 /mmc) and MgY (Pm3m).
Ab initio molecular dynamics (MD) simulation [29] was conducted to obtain liquid structures as well as their trajectories along the heating and cooling processes. Altogether, around 700 configurations (with each configuration typically containing 100 atoms) were selected and subjected to DFT calculations using the pseudopotential and plane-wave method implemented in the VASP [41,42]. The derived PES (the potential energy and stress tensors of each configuration, forces on each atom and elastic constants of two reference structures) was further modified to match the experimental values of the lattice parameter and cohesive energy of Mg and Y, and was then utilized to parameterize the EAM potential for the Mg-Y system. A similar practice can be found in [25]. During potential fitting, ad hoc EAM potentials were employed in classical MD to probe deeper potential basins on the PES, with new configurations added to the previously built potential database for a new round of EAM parameterization. Several iterations were performed until self-consistency between ab initio and EAM calculations was reached.
To demonstrate the overall performance of the as-developed EAM potential, we compare the present EAM model with previous models for Mg and Y in terms of the accuracy of predicting a set of material properties, as shown in tables 4 and 5. The performance of the EAM potential for Mg-Y alloys can be seen from figure 2, where the equations of state of the crystalline phases derived from the EAM model and the ab initio treatment are provided. The cohesive energies of the intermetallic compounds, with tabulated lattice parameters, are listed in table 6 for comparison. The general agreement between EAM and ab initio calculations is satisfactory. The as-developed EAM potential is available from https://sites.google.com/site/eampotentials/Home/MgY.

Results and discussion
Our quantum-mechanical and atomistic calculations of GSFE profiles include five slip systems for both pure Mg and Mg-Y. Specifically, we consider the three a -type slip systems  [30]. b Liu et al [31]. c Barrett and Massalski [34]. d Kittel [35]. e Simons and Wang [36]. f de Boer [37] g Average of basal plane and prism plane surface energies. h First-principles results as a part of this work. i Campbell [38].
({0001} 1120 , {1010} 1120 and {1011} 1120 ), and the second pyramidal slip system ({1122} 1123 ), which contributes to shear deformation out of the basal plane, as well as the {0001} 1010 slip system, which is related to the intrinsic stacking fault I 2 (I 2 SF). The I 2 stacking fault is formed when the [0001]-vector is altered by addition of 1 3 1010 . If we assign letters A, B and C to three possible stacking configurations of (0001) planes (in analogy to the stacking order of (111) planes in face-centered cubic lattices), the normal alternating hcpstacking (. . . ABABABAB . . . ) is locally changed to an (. . . ABABCACA . . . ) stacking. Thus the influence of Y on the energetics of the I 2 stacking fault can be studied by GSFE calculations. This part of our study represents a continuation of our previous experimental and theoretical study on the intrinsic stacking fault, I 1 (see [46] for details).     addition of Y. This finding is in agreement with previous DFT results [13] (0.288 J m −2 in Mg is reduced to 0.248 J m −2 in Mg-Y alloys). When studying the same GSFE employing atomistic EAM potentials, the trend is qualitatively the same but the reduction (from 0.233 J m −2 in Mg to 0.190 J m −2 in Mg 47 Y) is slightly smaller, 18.5%. The GSFE curves obtained by DFT and EAM calculations are similar with variations in the USFE ranging between 12% and 15% for Mg and the Mg-Y alloy. We can therefore conclude that both methods consistently predict the GSFE barriers to be reduced upon Y additions.  for elemental Mg is in good agreement with previous DFT studies for pure Mg that reported 0.218 J m −2 [16] and 0.225 J m −2 [17]. The corresponding EAM results are qualitatively very similar but the predicted reduction is lower. Therefore, despite the fact that the actual energy barriers obtained by DFT and EAM calculations differ, the changes induced by Y alloying are qualitatively very similar and thus independent of our selected computational method.   apparently lower than those obtained using DFT. However, the trend of a Y-induced reduction in the GSFE is the same. The DFT-calculated energies of local GSFE maxima and minima (I 2 SFE) are in quantitative agreement with previous results [13,16,17].

Pyramidal
In order to verify this theoretical prediction that Y additions lower the I 2 stacking fault, an experimental TEM study has been performed (for details related to TEM measurements, see our previous paper [5]). These TEM observations on slightly pre-deformed (1.5%) pure Mg and Mg-3Y (wt.%) were performed to measure the I 2 SFE. In pure Mg no I 2 stacking faults were observed, which indicates a relatively high stacking fault energy, i.e. a low probability to form these stacking faults or a dissociation width that is too small to be measured using conventional TEM. In Mg-3Y, frequent formation of the I 2 stacking fault was observed allowing the measurement of the I 2 SFE, figure 7. Here, the Burgers vectors (b) of the bounding partial dislocations of the stacking faults were determined according to the g · b invisibility criterion, Here, γ is the stacking fault energy, G the shear modulus, ν the Poisson's ratio, b the Burgers vector of the partials, β the angle between the partials and d the splitting width of the partials. Consequently, the I 2 stacking fault energy of Mg-Y amounts to 1.5 ± 0.5 mJ m −2 with an average dissociation width of about 20-30 nm.

Non-basal plane generalized stacking fault energies
As the {1122} 1123 slip system is an important non-basal slip system, its GSFE was also studied. The GSFE profiles of the slip system are calculated for pure Mg and Mg 47 Y (2.08 at.%) using both DFT and EAM calculations (see figure 8). The GSFE curves share similar shapes, with two maxima (one local and one global) and a local minimum along the shift vector. The local minimum corresponds to a stable stacking fault energy. The stable stacking fault energy predicted for Mg 47 Y by DFT (0.293 J m −2 ) is lower than that in pure Mg (0.318 J m −2 ). The EAM results suggest a smaller reduction in the stable SFE of only 0.008 J m −2 by the addition of Y. The energy difference corresponding to this reduction in the stable SFE is, however, very small and, in fact, quite close to the error bar of our EAM calculations.
The maximum of the GSFE (USFE) calculated by both DFT and EAM is higher in Mg 47 Y than in pure Mg. The DFT results are again very close to previous DFT results [14]. Both DFT and EAM data clearly indicate that the computed energies for Mg and Mg-Y crystals reverse their mutual order with increasing displacement (see figure 8). Specifically, the generalized stacking fault energies are first nearly equal for both Mg and Mg 47 Y for displacements lower than 0.2, then Y additions result in a clear lowering of stacking fault energies up to a displacement when GSFE curves of Mg and Mg 47 Y systems cross and the energies are beyond this point in the system containing the Y atoms. The displacement corresponding to the crossing point is ≈ 0.48 and 0.37 in the case of DFT and EAM results, respectively.
With regard to the local minima and maxima individual, the corresponding fractional displacements are not exactly equal when calculated by DFT and EAM methods.  [16] are the results of Sun EAM. The data taken from [17] are the results of KSDFT/LDA/TM-NLPS. a Muzyk et al [13]. b Yasi et al [16]. c Shin et al [17]. d Nogaret et al [14].
The DFT-calculated shallow minima (0.020 J m −2 for Mg and 0.015 J m −2 for Mg 47 Y) indicate a dissociation of the perfect 1 3 1123 dislocation family into partial dislocations according to the following reaction: where λ is a constant coefficient for the dissociation reaction. For λ = 0.4, the rows of atoms in the adjacent two gliding planes seem to be more regular than in any other case (see [19], figure 7, 0.6 c + a ). A dissociation of perfect c + a dislocations was also found by other authors, but with different λ values. In [14], based on DFT and EAM (Liu potential), values of λ = 0.33 and 0.25, respectively, were found for pure Mg. In [19], a λ value of 0.5 was reported using an EAM potential taking in-plane relaxations of the atoms into account. The different λ values are due to the different calculation methods and the specific conditions. Lower SFE and USFE values indicate a tendency toward deformation mechanisms including partial dislocations when compared with pure Mg.
Despite a few discrepancies between DFT and EAM results, it can be concluded that, from a qualitative point of view, the EAM gives a rather accurate description of the GSFE in the {1122} 1123 glide system. In particular, the trend that alloying Y in Mg reduces stable stacking fault energies is correctly captured. Further details related to the {1122} 1123 stable stacking faults will be described in the next subsection.
The calculated results obtained for the five studied slip systems are summarized in table 7 together with published data. The maxima along the GSFE profiles corresponding to the three basal slip systems indicate that the highest unstable stacking fault energy is expected for  The actual results are listed in table 7. Moreover, table 7 reveals that the ratio between the stable and unstable stacking fault energies decreases in the {0001} 1010 slip system upon alloying with Y, further supporting the conclusion that Y alloying lowers the SFE. For the {1122} 1123 slip system there is a more complex change of the positions and ratio of the stable and unstable stacking fault energies, and so the influence of Y atoms is more complex and a topic for future study.

Non-basal plane γ -surfaces
In order to deepen our understanding of 1123 non-basal slip processes and stable staking faults on the {1122} plane, the GSFEs in pure Mg and Mg-Y alloys were calculated employing our newly developed EAM potentials (see figure 9). There is no significant difference between the γ -surfaces for pure Mg and Mg-Y. In both figures 9(a) and (b), the lowest-energy slip path is a line section along y = 0, which is exactly the Burgers vector 1 3 1123 . Another path connects two maxima (marked as U1 and U2) and one minimum (marked S1) without consideration of the point where x = 0 and y = 0. This can be seen in figure 8. There are also another two minima (S2 and S3) in the γ surface, but their energies are both much higher than that of S1. So the S1 minimum is the most stable stacking fault on the {1122} plane. The calculation of the γ -surface for pure Mg was performed in [19] (see figure 4 in [19]), and quantitatively agrees with our results. Along [1010], there exists a curved minimum path across the S3 point. The fact that the energy at the S3 point is even higher than that calculated for U2 shows that this minimum path is not energetically favorable.

Conclusions
Using DFT and newly developed EAM potentials, we present GSFE profiles for five selected slip systems in both pure Mg and Mg-Y alloys. The results of the three a slip systems show that the USFEs decrease upon Y additions in Mg. From the slip system {0001} 1010 the I 2 SFE values for pure Mg and Mg-Y are obtained, and the I 2 SFE of Mg-Y alloy is again lower than that computed for pure Mg. Importantly, the generalized stacking fault energies associated with displacement along the 1123 direction in the non-basal plane {1122} are initially lower but with increasing displacement become higher than those of Mg. This is in contrast to the basal slip systems, for which the Mg-Y energies are consistently lower than for Mg. Lastly, after careful testing of our newly developed EAM Mg-Y potential, the generalized stacking fault energies for all glide systems within the {1122} plane are calculated using this potential and the GSFE profiles are visualized as a two-dimensional γ -surface. Our theoretical study has been complemented by an experimental TEM analysis that confirmed our theoretical prediction that the I 2 stacking fault energy is reduced due to Y additions.
To summarize, as a complement to previous experimental studies, we use theoretical methods to decompose the complex interplay of various mechanisms acting in Mg alloys in order to study some of them individually. Focusing solely on a few selected slip systems, we conclude that the impact of Y additions on generalized stacking fault energies is rather complex. On the one hand, the Y atoms reduce both stable and unstable stacking fault energies in the case of studied a slip systems. On the other hand, this influence may not be generalized to all slip systems as we predict also an increase of energies due to Y additions for most of the displacements along the 1123 direction in the non-basal plane {1122} and a shift in position of stable SF. Due to the fact that previous studies of fcc materials linked the stable and unstable stacking fault energies with the mobility and activity of dislocations, we speculate that our findings can have similar consequences in the case of Mg alloys. As a verification of this speculation will necessarily require further theoretical studies addressing the full complexity of these materials, we have also developed and carefully tested a new EAM potential. The comparison of first-principles and EAM results indicates that our EAM potential may be suitable for future larger-scale atomistic simulations.