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

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\}\langle 11\bar {2}0\rangle$ slip system, supercells with 2, 4, 6, 8 and 12 layers along the $\langle 11\bar {2}0\rangle$ 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 z-axis in the coordinate system). Therefore, in Mg₅₅Y the concentration of Y is 1.82 at.% and 2.08 at.% for Mg₄₇Y in all the other considered slip systems.

Table 2. Different stoichiometries with corresponding supercell sizes and slip systems in EAM calculations.

Slip systems	Supercell size	No. of atoms
$\{0001\}\langle 11\bar {2}0 \rangle$	2×2×6	48
$\{10\bar {1}0\}\langle 11\bar {2}0 \rangle$	2×2×7	56
$\{0001\}\langle 10\bar {1}0 \rangle$	4×4×6	192
$\{10\bar {1}1\}\langle 11\bar {2}0 \rangle$	4×4×3	192
$\{11\bar {2}2\}\langle 11\bar {2}3 \rangle$	12×6×2	1728

In the case of the γ-surface for the $\{11\bar {2}2\}$ plane, two perpendicular vectors ($\frac {1}{3}[11\bar {2}3]$ and $[10\bar {1}0]$) are selected on the plane. All displacements in the slip plane can be obtained by a combination of these two vectors. Suppose x and y are fractional displacements along the two directions; the GSF vector $\skew3\vec{f}$ may be written as $\skew3\vec{f} = x\frac {1}{3}[11\bar {2}3] + y~[10\bar {1}0],$ where 0 ⩽ x ⩽ 1 and 0 ⩽ y ⩽ 1. More details about the calculation of the γ-surface can be found in [19].

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⁻⁷ eV and 10⁻³ eV Å⁻¹, 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].

Table 3. DFT-computed lattice parameters for pure Mg and Mg–Y alloys with different stoichiometries with corresponding supercell sizes and slip systems.

Stoichiometry	a (Å)	c/a	Supercell size	Atom no.	Slip systems
Mg	3.1886	1.6261	2×2×2	16	All slip systems
Mg55Y	3.2004	1.6253	2×2×14	56	$\{10\bar {1}0\}\langle 11\bar {2}0 \rangle$
Mg47Y	3.2064	1.6215	2×2×12	48	$\{0001\}\langle 11\bar {2}0 \rangle$, $\{0001\}\langle 10\bar {1}0 \rangle$, $\{10\bar {1}1\}\langle 11\bar {2}0 \rangle$, $\{11\bar {2}2\}\langle 11\bar {2}3 \rangle$

Empirical potentials provide a means of exploring the physics of systems of atoms on length- and 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]: E_tot = Σ_i,jϕ(r_ij) + Σ_iF(n_i), where n = Σ_jρ(r_ij) 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₂₄Y₅ ($I\bar {4}3m$), Mg₂Y (P6₃/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.

Zoom In Zoom Out Reset image size

Table 6. Cohesive energies of Mg–Y intermetallic compounds predicted by EAM and ab initio calculations. The lattice parameters are taken from [28].

Compound	Lattice parameter	Cohesive energy ab initio (eV atom−1)	Cohesive energy EAM (eV atom−1)
MgY (Pm3m)	a=3.790 Å	3.06	3.05
Mg2Y (P63/mmc)	a=6.037 Å	2.55	2.64
	c=752 Å
Mg24Y5 ($I\bar {4}3m$)	a=3.790 Å	2.05	2.07

3.1.1. Basal-plane slip systems

The computed GSFEs of {0001} along $\langle 11\bar {2}0\rangle$ are plotted as a function of the shift vector for pure Mg and the Mg–Y alloy in figure 3. The USFEs based on DFT calculations for pure Mg and the Mg–Y alloy (Mg₄₇Y, 2.08 at.% Y) are 0.276 and 0.214 J m⁻², respectively. The USFEs are thus predicted to be reduced by 22.5% due to the addition of Y. This finding is in agreement with previous DFT results [13] (0.288 J m⁻² in Mg is reduced to 0.248 J m⁻² 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⁻² in Mg to 0.190 J m⁻² in Mg₄₇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.

Zoom In Zoom Out Reset image size

3.1.2. Prismatic $\{10\bar {1}0\}\langle 11\bar {2}0\rangle$ slip system

The GSFE curves of the prismatic 〈a〉 slip system $\{10\bar {1}0\}\langle 11\bar {2}0\rangle$ for pure Mg and the Mg–Y alloy (Mg₅₅Y, 1.82 at.% Y) are shown in figure 4. In the four GSFE curves, there is only one global maximum corresponding to the unstable stacking fault energy (USFE) of the studied slip system. Based on our DFT calculations, the USFE is 44% lower in the Mg–Y alloy (0.128 J m⁻²) than in pure Mg (0.231 J m⁻²). The latter value for elemental Mg is in good agreement with previous DFT studies for pure Mg that reported 0.218 J m⁻² [16] and 0.225 J m⁻² [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.

Zoom In Zoom Out Reset image size

3.1.3. Pyramidal $\{10\bar {1}1\}\langle 11\bar {2}0\rangle$ slip system

The GSFE results computed for the $\{10\bar {1}1\}\langle 11\bar {2}0\rangle$ pyramidal slip system are visualized in figure 5. The GSFE curves of the Mg–Y alloy calculated by DFT and EAM methods match very well. As in the case of the $\{0001\}\langle 11\bar {2}0\rangle$ and $\{10\bar {1}0\}\langle 11\bar {2}0\rangle$ slip systems, the GSFEs predicted for pure Mg differ slightly. However, irrespective of the computational method used, GSFE values are lower in Mg–Y than those determined for pure Mg. Hence the DFT-calculated USFEs in Mg–Y alloys are reduced by 9% (from 0.343 J m⁻² calculated for pure Mg). The corresponding EAM values indicate a reduction by 4% (from 0.323 J m⁻² in pure Mg). The DFT-based USFE of pure Mg is in very good agreement with previous DFT results [14, 17].

Zoom In Zoom Out Reset image size

3.1.4. The $\{0001\}\langle 10\bar {1}0\rangle$ slip system

The computed GSFE curves for the $\{0001\}\langle 10\bar {1}0\rangle$ slip system of pure Mg and Mg₄₇Y (2.08 at.%) are plotted as a function of the applied shift vector in figure 6. The four curves show similar characteristics, e.g. a local minimum corresponding to the displacement of ≈ $\frac {1}{3}$ and a global maximum at a shift of about $\frac {2}{3}$ of the perfect lattice vector. The local minimum corresponds to the local stacking configuration (...ABABCACACA ...), which is the I₂ SFE. The I₂ SFEs in pure Mg and Mg₄₇Y calculated by DFT are 0.037 and 0.030 J m⁻², respectively, and so they exhibit a reduction by 0.007 J m⁻² due to the Y atoms. The EAM values differ slightly but are qualitatively similar, with an I₂ SFE of 0.017 J m⁻² (in pure Mg) being reduced by 0.009 to 0.008 J m⁻² (in the Mg₄₇Y alloy).

Zoom In Zoom Out Reset image size

The global maximum of these GSFE curves corresponds to the stacking fault (...ABABABBCBCBC ...), which is characterized by the on-top position of the two atomic planes across the fault plane. The energy of this configuration in Mg and the Mg–Y alloy was found to be 0.475 and 0.445 J m⁻², respectively, using DFT. The EAM calculated energies of these maxima are 0.444 and 0.413 J m⁻², respectively, closely matching the corresponding Y-induced reduction in the USF energy predicted by DFT. Here we can conclude that irrespective of the material system (Mg or Mg–Y alloys), the GSFEs calculated by the EAM are 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₂ SFE) are in quantitative agreement with previous results [13, 16, 17].

In order to verify this theoretical prediction that Y additions lower the I₂ 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₂ SFE. In pure Mg no I₂ 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₂ stacking fault was observed allowing the measurement of the I₂ 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, figure 7. Both bounding partial dislocations are invisible under g = (0002), only one is visible under ${\bf g} = (1\bar {2}10)$ and both are visible under ${\bf g} = (\bar {1}010)$, ($01\bar {1}1$), ($\bar {1}101$) and ($11\bar {2}\bar {2}$), figure 7. Therefore, it is concluded that one partial dislocation has the Burgers vector ${\bf b} = 1/3 \ [10\bar {1}0$] and the other has the Burgers vector ${\bf b} = 1/3 \ [01\bar {1}0]$ building an I₂ stacking fault according to the dissociation reaction: 1/3 $[11\bar {1}0] \rightarrow 1/3 \ [10\bar {1}0] + 1/3$ [ $01\bar {1}0$]. The corresponding stacking fault energies are calculated based on the dissociation width of the partial dislocations according to equation

$$\begin{eqnarray*} &&\gamma = \frac{Gb^2}{8\pi d} \frac{2 - \nu}{1- \nu}\left( 1 - \frac{2\nu}{2-\nu} \cos{2\beta}\right). \end{eqnarray*}$$

Zoom In Zoom Out Reset image size

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₂ stacking fault energy of Mg–Y amounts to 1.5 ± 0.5 mJ m⁻² with an average dissociation width of about 20–30 nm.

As the $\{11\bar {2}2\}\langle 11\bar {2}3\rangle$ 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₄₇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₄₇Y by DFT (0.293 J m⁻²) is lower than that in pure Mg (0.318 J m⁻²). The EAM results suggest a smaller reduction in the stable SFE of only 0.008 J m⁻² 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.

Zoom In Zoom Out Reset image size

The maximum of the GSFE (USFE) calculated by both DFT and EAM is higher in Mg₄₇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₄₇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₄₇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. The DFT-calculated shallow minima (0.020 J m⁻² for Mg and 0.015 J m⁻² for Mg₄₇Y) indicate a dissociation of the perfect $\frac {1}{3}\langle 11\bar {2}3\rangle$ dislocation family into partial dislocations according to the following reaction:

$$\begin{equation} \frac{1}{3}\langle11\bar{2}3\rangle = \lambda~ \frac{1}{3}\langle11\bar{2}3\rangle + (1-\lambda) ~ \frac{1}{3}\langle11\bar{2}3\rangle ~(\lambda \simeq 0.4), \end{equation} \tag{ 1 }$$

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 $\{11\bar {2}2\}\langle 11\bar {2}3\rangle$ glide system. In particular, the trend that alloying Y in Mg reduces stable stacking fault energies is correctly captured. Further details related to the $\{11\bar {2}2\}\langle 11\bar {2}3\rangle$ 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 $\frac {1}{3}\langle 11\bar {2}0\rangle$ dislocations on the pyramidal plane and the lowest on the prismatic plane. The ratios between the stable and unstable stacking fault energies are calculated for the $\{0001\}\langle 10\bar {1}0\rangle$ and $\{11\bar {2}2\}\langle 11\bar {2}3\rangle$ slip systems. 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\}\langle 10\bar {1}0\rangle$ slip system upon alloying with Y, further supporting the conclusion that Y alloying lowers the SFE. For the $\{11\bar {2}2\}\langle 11\bar {2}3\rangle$ 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.

In order to deepen our understanding of $\langle 11\bar {2}3\rangle$ non-basal slip processes and stable staking faults on the $\{11\bar {2}2\}$ 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 $\frac {1}{3}\langle 11\bar {2}3\rangle$.

Zoom In Zoom Out Reset image size

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 $\{11\bar {2}2\}$ 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 $[10\bar {1}0]$, 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.