Temperature dependence of the stacking-fault Gibbs energy for Al, Cu, and Ni

The temperature-dependent intrinsic stacking fault Gibbs energy is computed based on highly converged density-functional-theory (DFT) calculations for the three prototype face-centered cubic metals Al, Cu, and Ni. All relevant temperature-dependent contributions are considered including electronic, vibrational, magnetic, and explicit anharmonic Gibbs energy contributions as well as coupling terms employing state-of-the-art statistical sampling techniques. Particular emphasis is put on a careful comparison of different theoretical concepts to derive the stacking fault energy such as the axial-next-nearest-neighbor-Ising (ANNNI) model or the vacuum-slab approach. Our theoretical results are compared with an extensive set of previous theoretical and experimental data. Large uncertainties in the experimental data highlight the necessity of complementary parameter-free calculations. Speciﬁcally, the temperature dependence is experimentally unknown and poorly described by thermodynamic databases. Whereas CALPHAD derived data shows an increase of the stacking fault energy with temperature for two of the systems (Cu and Ni), our results predict a decrease for all studied systems. For Ni, the temperature induced change is in fact so strong that in the temperature interval relevant for super-alloy applications the stacking fault energy falls below one third of the low temperature value. Such large changes clearly call for a revision of the stacking fault energy when modeling or designing alloys based on such elements.


I. INTRODUCTION
Metals and alloys-while crystalline-contain defects which largely determine their mechanical behavior, most importantly, plasticity and fracture. Examples include point defects such as vacancies, line defects such as dislocations, and planar defects such as grain boundaries. Of particular importance are planar defects known as stacking faults, more concretely the intrinsic stacking fault, which arises when a dislocation dissociates into so-called Shockley partials [1]. The extent of such partial dissociation is influential in determining the kinematics of dislocation flow in the presence of a stress field during plastic deformation [2].
It is the energy of the intrinsic stacking fault that is of importance in determining the extent of dissociation. When the stacking fault energy (SFE) is low, dislocation dissociation is favored and vice versa. Critical phenomena such as the kinetics of dislocation motion, work hardening, recrys-tallization, and crystalline texture-all of which influence the properties of these materials-are strongly dependent on the SFE size. Early experimental measurements of the SFE were prone to large scatter [3,4]. The advent of transmission electron microscopy enabled direct measurements of the extent of dislocation dissociation which somewhat improved the situation [4,5]. However, depending on the SFE magnitude, large uncertainties remain. Moreover, the intrinsic temperature dependence of the SFE has remained experimentally inaccessible for most elements.
Computational modeling is an attractive approach to overcome these difficulties. In particular, first principles based methods such as density functional theory (DFT) are now contributing a great deal to understanding and controlling materials properties [6][7][8]. In the past, an accurate determination of high temperature properties using DFT has been a great challenge due to the difficulty of capturing all the relevant entropic contributions. The majority of SFE calculations was thus restricted to T = 0 K or utilized approximate temperature dependencies (see the compilation in Appendix C). Meanwhile, methods have been put forward that enable an accurate calculation of finite temperature properties up to the melting point for perfect bulk and point defects [9][10][11][12][13]. An example is the DFT study of the temperature dependence of the Gibbs energy of formation of a vacancy, which mediates diffusional flow at high temperatures. These studies revealed that commonly assumed extrapolations from high temperature experimental data to low temperatures may fail qualitatively [9,12].
In this paper, we extend these methods to enable an accurate study of the SFE. We apply them to predict the temperature dependence of the SFE, with emphasis placed on the face-centered cubic (fcc) metals Al, Cu, and Ni. These elements have been selected since their alloys form the basis of many structural components which are used at elevated temperatures. Due to the critical impact of small free energy changes (in the range of meV/atom) on the temperature dependence of the SFE, special emphasis is placed on a detailed analysis of all relevant technical parameters.

A. General approach
Within DFT, the SFE can be computed by performing explicit stacking fault calculations [14][15][16], in which the stacking fault is modeled in an appropriate supercell. Our focus here is specifically on the temperature dependence of the intrinsic SFE at ambient pressure, given by the excess Gibbs energy where G sf (P , T ) and G fcc (P , T ) represent the Gibbs energies at pressure P and temperature T of the supercell including one stacking fault (sf) and the reference supercell of a perfect fcc structure, respectively. Further, A sf (P , T ) denotes the pressure and temperature dependent interface area over which the stacking fault extends in the (111) plane.
In the framework of first-principles calculations, it is convenient to calculate first the Helmholtz free energy, F α (V , T ) as a function of volume V , for each of the supercells (α=sf, fcc). The Gibbs energy is then obtained by a Legendre transformation: Note that the volume V on the right hand side of this equation is not a free variable but needs to be adjusted to correspond to the given pressure P (i.e., the negative derivative of the free energy surface with respect to volume). To account for the different temperature-dependent excitations, the free energy is decomposed as where E α denotes the conventional T = 0 K total energy, F el α the electronic free energy for the static lattice, F qh α the quasiharmonic free energy, F ah α the explicitly anharmonic free energy due to phonon-phonon interactions, F el-vib α the coupling contribution between electrons and atomic vibrations, and where (for Ni) with the magnetic contribution F mag α on the static lattice, its coupling to electrons F mag-el α and to atomic vibrations F mag-vib α . In Sec. II B we will discuss the methodological details related to the computation of the different free energy contributions. Sections II C to II E focus on a few particularly important technical aspects of the calculations. The treatment of the magnetic contribution for Ni will be detailed in Sec. II F. Note that the individual excitation mechanisms can usually not be treated independently. Rather, they can sensitively depend on each other as expressed by the various coupling terms in Eq. (3). For example, electronic and magnetic excitations can be affected by explicit atomic vibrations [17][18][19][20] and we will thus investigate these effects carefully.
Calculating the SFE according to Eq. (1) has the advantage of capturing local effects in the vicinity of the stacking fault such as, e.g., local atomic relaxations, local magnetic moments, or impurity interactions with or within the defect structure [15,21]. A typically quoted disadvantage is the increased computational costs due to larger system sizes and due to a reduced number of symmetries. Therefore, an often employed approximate alternative is the axial-next-nearestneighbor-Ising (ANNNI) model [22,23], in which the intrinsic SFE is expressed by a series expansion of different stacking sequences. In first and second order, the SFE can be formulated as where F hcp , F dhcp , and F fcc denote the Helmholtz free energies of the perfect hexagonal-close-packed (hcp; ABAB stacking along the [111] direction), double-hcp (dhcp; ABAC stacking), and fcc (ABC stacking) bulk structures. All these free energies are calculated at the fcc equilibrium volume V fcc := V fcc (P , T ) at a given pressure P and temperature T . For the hcp and dhcp structures the c/a ratio is set to the ideal value of √ 8/3 ≈ 1.633. The interface area of the mimicked stacking fault is derived from the fcc phase by A fcc (T ) = √ 3/4 [a(T )] 2 with a(T ) = [4V fcc (T )] 1/3 . For more details we refer to Appendix A.
Note that in Eqs. (5) and (6) only the free energies of perfect crystal structures are required. Thus, when considering only contributions on the perfect lattice in Eq. (3) [e.g., E α (V )], the ANNNI approach is computationally significantly less expensive than the explicit calculation in Eq. (1). This makes the ANNNI model attractive, e.g., for studying chemical trends [24,25]. Note however that the argument does not apply to the computation of free energy contributions that require per se large supercells, as the quasiharmonic or anharmonic one. In fact, in such a case the second order ANNNI model, Eq. (6), requires even more computational effort than Eq. (1) (hcp, dhcp, fcc vs sf, fcc). Nevertheless, due to its wide-spread usage, we will investigate the performance of the ANNNI model in detail.

B. Methodological details
The DFT total energies entering the free energies in Eq. (3) were calculated with VASP [26,27] employing the projectoraugmented wave (PAW) method [28] within the local density approximation (LDA) and within the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) parametrization [29]. The provided PAW potentials [30] were employed, treating the 3s 2 3p 1 , 3d 10 4s 1 , and 3d 8 4s 2 orbitals as valence electrons for Al, Cu, and Ni, respectively. T = 0 K total energies for the explicit SFE calculations were obtained employing a two-step procedure. First, the atomic positions in the stacking-fault supercell were optimized by performing an atomic relaxation parallel to the [111] direction. In a second step, the tetrahedron method in combination with the correction scheme proposed by Blöchl [28] was used to obtain accurate total energies. The T = 0 K energies for the perfect bulk structures were calculated likewise with the tetrahedron method. All energy-volume curves were parametrized by the Vinet equation of state [31]. Table I in Appendix B summarizes the relevant parameters for the total energy calculations.
Electronic free energies, F el α (V , T ), were computed based on finite-temperature DFT [32] for a dense mesh of temperature and volume points. The temperature-volume dependence was parametrized as discussed in Ref. [33]. We used in particular a fourth order fit for the temperature expansion of the density of states and a third order polynomial for the volume dependence of the free energy. Table II in Appendix B summarizes the relevant parameters for the electronic free energy calculations on the static lattice. Changes in the electronic free energy due to explicit vibrations [19] were fully included at the stage of the anharmonic calculations as mentioned below. This coupling contribution is important for Ni as will be shown in Sec. III A.
The phonon calculations required for F qh α (V , T ) were performed via the finite-displacement method as implemented in SPH/IN/X [34]. A comparison with linear response calculations will be given in Sec. II D. Supercells for the bulk structures included 256 and 288 atoms for fcc and (d)hcp. Only for the Al fcc calculations with LDA we increased the supercell size to 500 atoms to remove artificial imaginary frequencies close to the point. The supercell for the explicit SFE calculations included 180 atoms (3 × 5 × 2 supercell in terms of the primitive stacking fault cell, cf. Sec. II C) for both the stacking fault supercell and the perfect fcc reference supercell. The impact of choosing different supercell sizes on finite-temperature SFE computations will be discussed in Sec. II E. As discussed in Ref. [35], a dense grid for the augmented charges is crucial for the here considered metallic systems and was therefore carefully adjusted. All phonon calculations were carried out for a dense set of volumes and the volume dependence of F qh α (V , T ) was parametrized employing a third-order fit. Table III in Appendix B summarizes the relevant parameters for the quasiharmonic calculations.
The anharmonic free energies F ah α (V , T ) were computed using the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) approach [10,36]. The required interatomic potentials were parameterized by the embedded atom method (EAM) utilizing the MEAMfit code [37]. The DFT molecular dynamics (MD) simulations used as input for the potential fits were performed at the respective (experimental) melting temperature of the elements and at several volumes for at least 1000 MD steps. The potentials obtained for the fcc phase were also used for the stacking fault supercell calculations. Separate potentials were obtained for the hcp phase and also used for the dhcp phase. Overall the potentials showed a very good performance, in the sense that the standard deviation of the energy difference between DFT and potential was only 1-3 meV/atom across all volumes and temperatures. Due to this good performance, the thermodynamic integration from the potentials to DFT (with low converged parameters) could be efficiently performed, using five λ values, i.e., 0, 0.15, 0.5, 0.85, and 1, and a cubic fit. The thermodynamic integration between the quasiharmonic and the optimized potentials was instead performed on a dense grid of coupling parameters λ, i.e., at least 20 λ values, to capture the nonlinear dependence of the integrand. All MD simulations were run until a standard error of well below 1 meV/atom was reached. The Langevin thermostat was used with a friction parameter of 0.01 fs −1 and a time step of 5 fs. The impact of the atomic vibrations on the electronic free energy was fully taken into account at the upsampling step of the TU-TILD method (i.e., high converged calculations) by adjusting the electronic temperature (i.e., Fermi-Dirac broadening) to the temperature of the atomic vibrations. Following the analytical formula proposed in Ref. [38], the explicitly computed anharmonic free energies for a set of volume and temperature points were used to fit a smooth anharmonic free energy surface F ah α (V , T ) in terms of a renormalized frequencyω ah (V , T ) = a 0 + a 1 V + a 2 T , where a 0 , a 1 , and a 2 are fitting coefficients. Table IV in Appendix B summarizes further relevant parameters for the anharmonic calculations. The calculation strategy for the magnetic contribution F mag α (V , T ) (for Ni) as well as an analysis of the impact of magnetism on other finitetemperature contributions such as vibrations will be discussed in Sec. II F and Sec. III B, respectively.

C. Supercell geometry and relaxation effects at zero temperature
Different geometrical setups exist to compute the SFE within a periodic boundary approach [39]. Two frequently chosen setups are (i) a tilted supercell geometry [ Fig. 1(a)] and (ii) a slab+vacuum geometry [ Fig. 1(b)]. Within the first setup a tilt is introduced to one of the lattice vectors in order to create the stacking fault while preserving periodic boundary conditions. An argument against this setup is that the tilt leads to supercell geometries which are difficult to converge to a good precision [39]. Therefore in setup (ii), a technically simpler supercell geometry is preserved at the expense of creating a vacuum to comply with the boundary conditions [40][41][42]. The unavoidable surfaces accompanying the slab+vacuum approach lead, however, to spurious finite size effects on the wave function and to a modification of the SFE. We show in the following that setup (i) is superior to setup (ii) in efficiently obtaining well converged SFE results.
We first concentrate on the results for the tilted supercell according to setup (i) shown by the black lines in Fig. 2. The SFE is shown as a function of the separation between neighboring stacking faults [cf. Fig. 1(a)]. The critical question is, how far the stacking faults need to be separated to avoid interaction between them. We observe from Fig. 2 that for all three considered elements a separation of about 7 Å (corresponding to a 1 × 1 × 1 stacking fault supercell with six atoms) is sufficient to reach an SFE to within 10% of the converged value. A twice as large separation (1 × 1 × 2 supercell with 12 atoms) gives SFEs to within 1% of the converged value and can be thus considered as a well converged geometry. An interesting observation is that relaxation (difference between dashed black and solid black lines in Fig. 2), which can be responsible for a reduction in the SFE of up to 3%, is not well captured by the smaller supercell geometry with only six atoms. The reason is that the smaller supercell introduces constraints on the atoms, i.e., the too close periodic image atoms prevent relaxation. We turn now to the results for the slab+vacuum approach. We restrict the discussion on Al. A converged vacuum region of ≈14 Å is used. The corresponding unrelaxed SFE is shown in Fig. 2(a) by the red-dashed line as a function of the separation between the stacking fault and the surface [see Fig. 1(b)]. This energy should be compared with the unrelaxed SFE for the tilted supercell approach given by the black-dashed line. We observe that for a separation of about 14 Å, for which the tilted setup was very well converged, the slab+ vacuum approach shows an error of more than 20%. This discrepancy indicates that there is a strong interaction between the surface and the stacking fault. Only by further increasing the separation between the stacking fault and the surface, the SFE converges to the same value as for the other approach.
The plot as a function of separation between the respective defects shadows the actual computational time involved in the SFE calculations. For the same separation, the actual number of atoms involved in the slab+vacuum setup is in fact two times the number of atoms for the tilted supercell setup as shown in Fig. 2(b). Since it is the number of atoms that mainly determines computational time, we clearly see from Fig. 2(b) that the tilted supercell setup is by far more efficient than the slab+vacuum one.
There are in fact other problems with the slab+vacuum approach related to atomic relaxation and vibrations. Our test calculations have revealed that atoms close to the surface strongly relax either into the vacuum or into the bulk depending on the actual lattice constant. Fixing the surface atoms can be used to circumvent this problem but it introduces an artificial constraint on the system that might affect the final SFE. Furthermore, when calculating finite-temperature vibrations, a surface acts as a hard wall reflecting phonons back and causing an interference with other phonons and possibly standing waves. These artifacts might again have an effect on the SFE, in particular on its temperature dependence. All of our explicit SFE results shown in the following are therefore based on the tilted-supercell approach.

D. Finite displacements versus linear response
The dynamical matrix, from which phonons, the quasiharmonic free energy, and thus quasiharmonic SFE are derived, can be computed by finite displacements on the Born-Oppenheimer energy surface or alternatively by linear response calculations within perturbation theory. For both approaches very high convergence parameters are required to obtain a precision allowing us to resolve the temperature changes of the SFE. The energy cutoff, k-mesh sampling, and augmentation grid used for our calculations are detailed in Appendix B in Table III. Here, we discuss the influence of the convergence criterion for the electronic loop showing that it needs to be extremely narrow and that the finite displacement method converges better than the linear response approach.
For the finite displacement method we used a displacement of 0.02 Bohr radius which provides a decent compromise between numerical accuracy and harmonic displacement. For the linear response calculations we used the implementation in VASP which employs a virtual displacement of ≈10 −4 Bohr radius for the perturbation of the wave function. It should be noted that the linear response theory in VASP applies the perturbation in real space rather than in reciprocal space in contrast to the frozen phonon method (or sometimes also called linear response method) [43].
Convergence results for the harmonic free energy and SFE at the melting temperature are shown for both methods in Fig. 3 as a function of the energy convergence criterion for the electronic loop for the example of fcc Ni. The latter has been chosen as it displays the largest absolute changes in SFE with temperature. We observe that the finite displacement method (black lines) converges quicker with the electronic-loop break condition than linear response theory (red lines). For the finite displacement method, a break condition of 10 −6 eV provides already free energies to within 1 meV/atom accuracy, whereas for the linear response calculation a value of 10 −7 eV should be rather chosen to obtain a similar precision. The faster convergence can be understood by comparing the magnitudes of the displacements in both methods. Within the linear response method a very small displacement is required for the perturbation theory to function properly. The small displacement leads however to a significant sensitivity of the results with respect to numerical noise. For the finite displacement method a larger displacement can be chosen and this increases numerical stability leading thus to the observed faster convergence. The fact that both approaches converge to the same free energies and SFE ensures that the displacement of the finite displacement method is still in the harmonic regime.
As a general rule, we suggest to always employ an energy break condition of 10 −7 eV for quasiharmonic calculations regardless of the method, system, and supercell size. We have observed that this value guarantees a well-converged dynamical matrix and resulting quasiharmonic properties throughout our calculations with only a minor increase in computational time as compared to a break condition of 10 −6 eV. It is important to stress that the suggested convergence criterion should not be rescaled with the number of atoms. The reason is that the energy derivatives determining the dynamical matrix (both within the finite displacement and the linear response method) depend on the absolute energy convergence. 1 1 Note that in the implementation of the linear response approach within VASP, specifically in ILINEAR_RESPONSE.F, an exit condition  Figure 4 shows the impact of the supercell size on the quasiharmonic SFE. Except for the first order ANNNI SFE of Ni, we observe a rather quick convergence. Note the different energy and temperature scales for the three considered elements in Fig. 4. Deviations between the different shown supercell sizes are at most 7 mJ/m 2 regardless of the temperature, which translates into about 1 meV/atom when considering the ANNNI SFE. This result shows that quasiharmonic SFEs can be calculated accurately already in rather small supercells with a few tens of atoms. A more detailed analysis of our data reveals that the underlying reason for this are rather short-ranged force constants of the involved phases.

E. Impact of supercell size at finite temperatures
A clear exception is the first order ANNNI SFE for Ni. At the melting temperature of 1728 K, the SFE is reduced for the linear response cycle is included when the initial norm of the residual vector, RMS, does not change within three iteration cycles by more than 10%. This coincidentally caused for some virtual displacements the linear response cycles to stop before the chosen minimization criteria has been reached (e.g., at 10 −3 eV instead of 10 −6 eV). This could be circumvented by commenting out the corresponding part in the routine. The harmonic curve (gray) corresponds to a 'usual' finite displacement calculation of the dynamical matrix, i.e., with single atoms displaced away from their perfect hcp positions. The oscillations indicate long ranged interactions, which are however compensated by the explicitly anharmonic contribution (red curve), i.e., by explicit vibrations of all atoms. To obtain the anharmonic value for the largest supercell size (490 atoms) we assumed the same change in the free energy between the optimized EAM potential and DFT as for the 288 atoms supercell, since this change was already well converged (below 1 meV/atom) with supercell size. by −55 mJ/m 2 (−9 meV/atom) when going from the 32/36 atoms to the 108/96 atoms calculation, and then it increases again by 36 mJ/m 2 (6 meV/atom) when going to the 256/288 atoms calculation; cf. the green and black arrows in Fig. 4. Additional (harmonic) calculations for an even larger supercell combination of 500/490 atoms were necessary to clarify that the SFE value obtained for the 256/288 atoms calculation is eventually converged (to about 2 mJ/m 2 ).
The main reason for the slow convergence of the first order ANNNI SFE for Ni with supercell size are long-ranged interactions in the dynamical matrix of the hcp phase. This is exemplified in Fig. 5 by the gray line. Interestingly, the corresponding explicitly anharmonic free energy shows a similarly slow convergence behavior (red line) but with a compensating dependence. This means that if we do not separate the harmonic and anharmonic contributions, we obtain a comparatively smooth full vibrational free energy as a function of the supercell size (black curve). In fact, the full vibrational free energy is converged down to 1 meV/atom already for the smallest supercell size with 36 atoms. We thus observe that long-ranged interactions, which are present in the (quasi)harmonic dynamical matrix derived from a perfectly ordered lattice, are removed by explicit vibrations. In this sense, including the explicitly anharmonic contribution renders the calculations easier.
However, there is another issue which requires in some cases the application of larger supercell sizes for the DFT MD anharmonic calculations. We have observed that for MD runs performed at high temperatures and large volumes, there is a significant chance that during the MD a whole plane of atoms shifts and then vibrates around new positions. This collective movement of a whole atomic plane is visualized in Fig. 6 where one (0001) B plane of the hcp structure transforms into a C plane of dhcp (cf. atomic trajectories shown in black and gray). Since such a collective shift requires a strictly correlated, simultaneous movement of all atoms in the plane, the probability of the shift is strongly suppressed for larger supercell sizes. Our calculations reveal that a 4 × 4 × 3 supercell with 96 atoms for hcp Ni and a minimum size of 2 × 3 × 2 with 72 atoms for the explicit stacking fault calculations is sufficient to render the calculations for the required number of MD steps stable [see Table IV (column "cell" under the "potential → DFT" category) for the safe supercell size for each phase].

F. Magnetic treatment
Although modeling of the magnetic free energy from ab initio has been significantly advanced over the last years (see, e.g., Refs. [44][45][46]), it still remains a challenging task. This applies specifically to complex geometries and materials which are prone to longitudinal spin fluctuations (LSFs), i.e., temperature induced variations in the magnitude of the local magnetic moments. The key issues are: (i) to construct and solve a magnetic model, (ii) to include the effect of magnetism on other degrees of freedom (e.g., atomic vibrations), and (iii) to incorporate the reverse effect from the other degrees of freedom. The situation becomes even more complicated when (i)-(iii) have to be solved self-consistently, i.e., when the degrees of freedom cannot be adiabatically decoupled. Some progress has been made by combining spin dynamics and MD [20], but such approaches are not directly applicable to Ni because of LSFs.
In more general terms, approaches utilizing directionally constrained local magnetic moments cannot be directly applied to materials prone to LSFs. In the case of Ni the underlying reason is that a standard supercell calculation with disordered magnetic moments will converge into the nonmagnetic solution [47,48]. Attempts to tackle this challenge utilize, e.g., an effective Heisenberg model fitted to the experimental Curie temperature [49] or an extended Heisenberg model with LSF terms fitted to magnitude-constrained spin calculations [47]. These approaches are, as of now, not applicable to complex geometric structures (such as the stacking fault) and cannot be coupled with MD simulations.
To overcome this difficulty and still capture the relevant physics, we employ here a semiempirical approach informed with spin-polarized DFT MD calculations, which has been recently shown to resolve subtle energy differences and phase stabilities [50,51]. The approach exhibits a qualitatively correct temperature dependence below and above the Curie temperature, T C , and also in the high temperature limit. Specifically, we employ an empirical function for the magnetic specific heat developed within the CALPHAD approach [52], [1120] with 105, and p = 2, and the magnetic entropy where k B is the Boltzmann constant and m a magnetic moment. The magnetic free energy is readily obtained via integration of Eq. (7). The two parameters determining the magnetic free energy in the above approach are the Curie temperature, T C , and the magnetic moment, m. For the Curie temperature we take the experimental value [53] and discuss in Sec. III B the possible impact of the stacking fault on T C . The influence of atomic vibrations on T C is neglected which can be justified by a recent study [54] which showed that T C derived from magnetic exchange parameters from MD simulations is hardly affected in the case of Ni. The magnetic moment m entering Eq. (8) is obtained from spin-polarized DFT MD simulations in the ferromagnetic state. In this way the impact of atomic motion on the site-averaged local magnetic moments, i.e., m MD (V , T ), is included as well as the impact of the thermal expansion, cf. the coupling term F mag-vib α in Eq. (4). The usage of ferromagnetic calculations in the paramagnetic regime is justified by recent studies for Ni [12,47,48,54]. The local magnetic moments obtained for a paramagnetic state stabilized by longitudinal spin fluctuations are close to the ferromagnetic value [47]. Further, the impact of magnetism turns out to be negligible on other finite temperature contributions such as, e.g., lattice vibrations [48] or vacancy formation energies [12]. Specifically, in Ref. [48] the impact of magnetic fluctuations on the force constants for fcc Ni was extensively studied including longitudinal and transverse spin fluctuations. It was found that both magnetic degrees of freedom hardly affect the force constants and hence have a negligible impact on vibrations in fcc Ni. Ferromagnetic and non-spin-polarized MD simulations for fcc Ni [12] further corroborated this finding, i.e., that the magnetic fluctuations apparently do not alter the sampled atomic phase space for Ni. The results of the present study discussed in Sec. III B will likewise support these findings. Note, however, that these statements are element specific and that the situation is different, e.g., for Fe [17,55] and Fe-based alloys [56]. Figure 7 shows our main results for the temperature dependence of the SFEs for Al, Cu, and Ni in the first, second, and third row, respectively. The different lines and colors emphasize the various finite-temperature contributions. The thick black solid lines are the final SFEs for Al and Cu. For Ni the red dashed-dotted lines are the final SFEs, highlighting additionally the (small) contribution from magnetic excitations. The first two columns show a comparison between LDA and GGA-PBE. Columns two to four present a comparison between the first and second order ANNNI model as well as the explicit calculations. We note the overall small energy scale in the range of a few meV/atom for which temperature-induced changes in the SFE need to be resolved (see red markers/numbers). This reveals the genuine importance of a high-accuracy approach as described in the previous sections.

A. Temperature dependence of the SFE: Mechanisms and approximations
A general feature that we observe irrespective of the functional, the SFE model, and the element is a strong decrease of the SFE with increasing temperature. Consequently, at the respective melting temperature the initial T = 0 K SFEs are significantly reduced. For Al, the reduction can be up to 60%, for Cu it is stronger with up to 70%, and for Ni even up to 85%. Thus, it becomes clear that it is critical to consider carefully the temperature dependence of the SFE, for example, when using it as an input to higher level modeling.
The main contribution causing this reduction is, for all cases, due to the quasiharmonic phonons (i.e., blue and green shaded regions taken together). Since the SFE is computed from a difference where the (perfect) fcc crystal structure is always the reference structure irrespective of the model [cf.
Eqs. (1), (5), and (6)], we can conclude that the phonons of the fcc structure are stiffer than the ones of the other involved structures. Explicit analysis of our data reveals that this is indeed the case, and that the averaged phonon frequencies in fcc are by 1-2 meV higher in energy. This finding can be intuited as follows. All considered elements have a positive SFE at T = 0 K, i.e., the (perfect) fcc crystal structure is the ground state at zero K and energetically more stable than the geometrically similar, also closed-packed hcp and stacking fault containing supercells. This suggests a stronger electronic binding (alas force constants) for the fcc ground state and thus stiffer phonons. Hence, the competing hcp and stacking fault supercells have softer phonons and thus a larger vibrational Gibbs energy contribution resulting in a decreasing SFE, i.e., competing with the initial T = 0 K fcc stability.
The finite temperature excitations beyond the quasiharmonic phonons constitute only a small correction to the SFEs, in particular for Al and Cu. Ni is somewhat special as it shows a stronger electronic contribution and likewise a compensation effect from the electron-vibrational coupling. These points will be analyzed in detail below. Before, we turn to the impact of the exchange-correlation functional and the SFE model.
Comparing the first column with the second, we see a similar temperature dependence of the SFE for LDA and GGA for each element. This statement does not only apply to the final SFE but also separately to each contribution of the various excitation mechanisms. The main difference between LDA and GGA lies in the T = 0 K SFE which is always higher in LDA, by about 10 to 20 mJ/m 2 . The observation that the difference between LDA and GGA is mainly due to the T = 0 K energy surface whereas the temperature dependence is similar is consistent with previous findings for other properties, such as thermal expansion and heat capacities [35] and vacancy formation energies [9]. 16   Comparing columns two to four we observe, overall, that the first and second order ANNNI models produce similar temperature dependent SFEs as the explicit calculations. The contribution of the different excitation mechanisms is, mostly, similar as well. A general feature is visible for the T = 0 K SFEs which are highest for the first order ANNNI model and which consistently converge to the lower explicit SFE values when going from the first via the second order approximation to the explicit model. For Al, the convergence can be seen also for the temperature dependence: The reduction of the SFE at the melting point is −79 mJ/m 2 for the first order ANNNI model, −72 mJ/m 2 for the second order ANNNI model, and −70 mJ/m 2 for the explicit SFE. However, some features are not resolved by the ANNNI approximations. For example, in Cu the anharmonic and coupling contributions show a qualitatively different behavior for the explicit SFE.
We come back to the discussion of the finite temperature contributions beyond the quasiharmonic phonons, focusing first on anharmonicity (i.e., difference between dashed orange lines and thin solid black lines in Fig. 7). The underlying anharmonic free energy surfaces of the various contributing structures are shown in Fig. 8. A general observation is that the overall magnitude of the anharmonic free energy is relatively small with only a few meV/atom even at the highest temperatures and largest volumes. Further, comparing the different structures for each element with each other we find a very similar anharmonic free energy dependence, resulting in the already mentioned, small contribution to the SFEs. Hcp Ni, in particular within GGA, falls slightly out of this picture showing a positive temperature dependence of the anharmonic free energy in contrast to the other structures. The resulting anharmonic contribution to the SFE in the first order ANNNI model is thus strongest with 18 mJ/m 2 at the melting point. This behavior is likely related to the compensation of the long-ranged interactions in the dynamical matrix as discussed in Sec. II E.
Electronic excitations (orange shading in Fig. 7) decrease the SFE of all three metals, only little in the case of Al and Cu, whereas a significant reduction is found for Ni. The underlying physics is explained in Fig. 9 for the example of the first order ANNNI model. The electronic contribution is determined by (i) the considered temperature, which provides the energy window of the excited electron-hole pairs (gray broadening function in Fig. 9) as well as by (ii) the electronic density of states (DOS) of the considered structure (green and red lines). For Al and Cu, we observe in general a small electronic DOS in the relevant energy region at the Fermi level (0.4 states/eV atom for Al and 0.3 states/eV atom for Cu) leading to small electronic free energies for each structure (a few meV/atom at the melting point). Changes between Ni behaves differently because, for the minority spin channel (Fig. 9), the relevant energy window falls together with a sharp peak in the electronic DOS. The stability of the fcc structure at T = 0 K is achieved by the possibility to shift some of the states at the Fermi level to an energetically lower peak, about 1 eV below the Fermi level as indicated by the red arrow. This is not the case for the hcp structure, and thus the electronic DOS of the hcp at the Fermi energy is significantly larger than the one of fcc (difference of about 0.4 states/eV atom). Therefore the electronic contribution compensates the T = 0 K stability of the fcc phase and favors the hcp structure with increasing temperature.
The electronic DOS itself depends, however, also on the atomic temperature, i.e., it can change due to atomic fluctuations [F el-vib α term in Eq. (3)]. Note that we have already taken the impact of the electronic temperature on the electronic DOS into account as discussed in Ref. [19]. The coupling between vibrational and electronic degree of freedom has so far only implicitly been taken into account via the volume dependence of the electronic DOS. Now, by switching "on" the atomic motion at finite temperatures employing molecular dynamics simulations, the electronic DOS becomes also explicitly temperature dependent, not only via the change of atomic volume but also by the atomic temperature.
Interestingly, atomic vibrations do not only have a strong impact on the electronic contribution to the SFE of Ni (change between thin and thick solid black lines in Fig. 7), this coupling contribution in fact compensates to a large extent the just discussed electronic contribution for the static lattice. The underlying reason for this is a strong broadening of the electronic DOS as illustrated in Fig. 10. Specifically, the broadening results in electronic DOS' that are extremely similar for fcc and hcp, an observation consistent with recent results for fcc and bcc Fe [18] and other d-transition elements in the fcc, bcc, and hcp structures [19]. The free energy change due to the coupling turns out to be much larger in the case of hcp than fcc (9 meV/atom vs 1 meV/atom at the melting point). This strong impact for hcp can be traced back to a stronger suppression of the DOS at the Fermi level with temperature as compared to fcc [cf. Fig. 10(b)]. Noticeably, even an artificial fcc-hcp transition (indicated by a negative value for the first ANNNI model) would be obtained if this higher-order contribution was not taken into account. We can therefore conclude that a proper treatment of anharmonic contributions and inherent coupling terms is critical for the consideration of phase stabilities [18,19].
We have also evaluated the reverse impact of this mutual interaction to address the question, if a finite electronic temperature can cause a different sampling of the phase space. We find, however, that the electronic temperature has no significant impact on the phase trajectory of the MD runs themselves which is consistent with our inherent assumption of the UP-and TU-TILD methods of a robust phase space sampling employing computationally efficient parameters.

B. Special considerations for Ni: Estimating the impact of LSFs and coupling terms
For all scenarios considered in Fig. 7 (LDA, GGA, ANNNI, explicit SFE), the impact of magnetism in Ni iswithin the utilized magnetic model-very small (compare black solid and red dash-dotted lines). The reason are compensating contributions for the defect/hcp and ideal fcc structures, i.e., similar local magnetic moments and critical temperatures. This can be verified in Fig. 11(a)  the magnetic free energy contributions for both structures are almost identical, with a maximum difference of less than 1 meV at the melting temperature [blue line in Fig. 11(b)]. This indicates that the influence of magnetism on the SFE in Ni is negligible and thus that the SFE temperature dependence is well captured by the present treatment.
As discussed in Sec. II, longitudinal spin fluctuations are only implicitly included in the applied empirical approach (via the self-consistently determined volume-dependent local moments including the impact of atomic vibrations). However, in general they are known to be of importance for an accurate computation of the high-temperature paramagnetic state of fcc bulk Ni [47]. While a full consideration of LSFs [47] on the SFE is at present not feasible, we provide here several supporting arguments that these effects are unlikely to modify the obtained SFE temperature dependence.
For the magnetic contribution, the key parameters are the local magnetic moments, m, and the Curie temperature, T C . One might assume that in the vicinity of the defect the local magnetic structure shows a different tendency for magnetic ordering than in the bulk due to different magnetic interactions within or close to the defect. Within the first order ANNNI model, such a scenario would correspond to a different virtual Curie temperature of hcp Ni. As pointed out in Sec. II F, a reliable calculation of T C for Ni requires the inclusion of LSFs as revealed for fcc Ni by Ruban et al. [47]. Similar computations for hcp Ni are beyond the present scope and we therefore estimate the virtual Curie temperature of hcp Ni, T hcp C , by resorting to a mean field approach [51] Here, T hcp,MF C and T fcc,MF C denote the conventional mean field expressions [44,46,57,58] for the Curie temperature of hcp and fcc which, apart from physical constants, depend only on the sum of all magnetic exchange interactions, i J 0i [44]. We have computed the latter utilizing the magnetic force theorem and the exact muffin-tin orbitals (EMTO) approach [59,60] employing the Lyngby version [61] of the code. Using the experimental Curie temperature for fcc Ni, T fcc C = 633 K [53], we find a slightly reduced T hcp C ≈ 609 K indicating overall slightly reduced magnetic interactions in hcp compared to fcc Ni. Inserting these values into Eq. (7), the impact of the slightly different T C values turns out to be no more than 1 meV up to the melting temperature. Note that this estimation provides an upper limit of the effect of an altered Curie temperature because the two-dimensional stacking fault defect embedded in the fcc host will likely adapt to the global magnetic ordering of the fcc host.
The second quantity entering the magnetic model are the local magnetic moments, m. Our FM calculations reveal very similar magnetic moments for all considered structures (fcc, hcp, dhcp as well as explicit stacking fault). This is one of the main reasons why the magnetic contributions within the current treatment cancel each other providing negligible impact on the SFE. In principle the inclusion of the LSF degree of freedom on the local magnetic moments might change this picture, i.e., the local magnetic moments in the paramagnetic state could differ more strongly between fcc and hcp compared to the ferromagnetic state. Indeed, as already mentioned, LSFs are crucial for the stabilization of local magnetic moments in paramagnetic Ni, i.e., neglecting LSFs, a standard disordered-local-moment (DLM) calculation (e.g., by setting up a large cell with random magnetic moments) will converge into a nonmagnetic solution. To estimate the local magnetic moments in the paramagnetic state of Ni we therefore employed a recently developed LSF theory [62] and computed the effective local moments of Ni at high temperatures (1600 K) in the hcp and fcc structure. These calculations reveal that the difference between the fcc and hcp derived local moments is almost unaffected and similar to the FM calculations. Based on this we can estimate the change of the magnetic energy contribution due to the thermally LSF enhanced local magnetic moments to be <0.3 meV. This further indicates that the inclusion of LSFs is not decisive for a correct description of the temperature dependence of the SFE of Ni.
We now give further supporting arguments that the vibrational contribution in fcc Ni is hardly affected by magnetism. As already discussed above, vibrations contribute the main part to the temperature dependence of the SFE. However, as outlined in Sec. II F, for fcc Ni the application of an extended version (taking into account LSFs) [48] of the spin-space averaging technique [55] has revealed that neither transverse nor longitudinal spin disorder impact the force constants and phonon spectra [48]. Even nonmagnetic, i.e., non-spinpolarized calculations turn out to provide similar results for phonon spectra of fcc Ni compared to the ferromagnetic ones. This suggests that for Ni the impact of magnetism including LSF on vibrational contributions to the SFE temperature dependence is also negligible.
In order to corroborate this assumption we have performed nonmagnetic (NM) calculations for the SFE. Such calculations are sometimes employed to model the paramagnetic state but they often provide unsatisfactory quantitative results. Qualitatively this approach can be nevertheless used to estimate the impact of magnetism on different quantities, e.g., vibrations [48]. To obtain a similar estimate for Ni, we have performed the complete set of calculations for all Gibbs energy contributions entering Eq. (3) for the NM scenario. Consistently with our phonon calculations of paramagnetic fcc Ni [48], we find that the overall temperature dependence within the first order ANNNI model is similar for the NM and FM scenarios as shown in Fig. 12. The individual contributions differ over the whole temperature range by not more than 20 mJ/m 2 (≈2 meV). Such a close temperature dependence for the NM and FM vibrational and electronic contributions indicates that magnetism, in particular also the coupling to vibrations, cancels each other for the computation of the SFE in Ni.
Finally we have investigated the combined impact of vibrations and LSFs on the electronic free energy contributions, cf. coupling terms in Eq. (4). For that purpose we selected snapshots from our FM MD runs for fcc (32 atom cell) and hcp (36 atom cell) at the largest considered temperature (i.e., 1728 K) and volume. These snapshots were employed in subsequent EMTO calculations [63] in combination with the recently developed LSF theory [62] as also discussed in Ref. [54]. The electronic free energy difference obtained in the FM and LSF-stabilized paramagnetic state turned out to be less than 0.5 meV/atom. The underlying reason can be observed in Fig. 13 which shows that the thermally averaged DOS' including the effect of LSFs are almost indistinguishable between fcc (red line) and hcp (green line).

C. Comparison with experiment and thermodynamic databases
In order to assess our theoretical prediction, we have conducted an extensive literature survey of previous experimental as well as theoretical data (see Appendix C). A representative set of data is summarized in Fig. 14 and, focusing on experiment, we observe a very strong scatter for all three metals. Among the various techniques available to measure the SFE, the more recent weak beam and high-resolution transmission-electron-microscopy techniques are considered as most reliable [4,5]. The corresponding data points in Fig. 14 have been given a black outline. Interestingly, these more accurate experimental SFE predictions give typically lower values than the earlier methods, in much better agreement with DFT calculations. Nevertheless, depending on the SFE magnitude, large uncertainties remain as evidenced by the error bars for the Al SFE.
The scatter in between theoretical SFE values obtained with different empirical potentials-mostly based on the EAM approach-is enormous. For example, for Ni an SFE as low as 13 mJ/m 2 is predicted at T = 0 K by one potential (OJ Zimmermann 2000 [66]), whereas another potential results in an SFE of 304 mJ/m 2 [67] (above the scale shown in Fig. 14; cf. Table VII), i.e., more than an order of magnitude larger. One has to be aware that such differences in the SFE will definitely change the dominant deformation mechanism when employing the potentials in large scale MD simulations.
The scatter in between previous DFT based SFE predictions is significantly smaller. Our zero K derived DFT values agree well with the previously reported ones. Small deviations are caused by the employed flavor of the exchange-correlation functional (LDA vs GGA) or other technical details. For example, employing the vacuum slab approach, previously reported values, e.g., for Al, are higher than the ones obtained in setups not introducing a vacuum [41,42]. This can be attributed to the inherent deficiencies of the vacuum slab approach as discussed in Sec. II C. In fact, for our smallest considered vacuum-slab supercell (Fig. 2), the obtained value of 151 mJ/m 2 for Al is in good agreement with previous vacuum-slab based calculations [41,42]. But note that this higher value is not converged with respect to the chosen supercell as exemplified in Fig. 2.
Only very few high temperature experimental data for Cu exist. It needs to be stressed that the corresponding technique ("2× coherent twin energy" in Table VI) is an indirect method relying on model assumptions. The scatter in these high temperature data is similar as for the low temperature values and it is impossible to extract any meaningful temperature dependence of the SFE. For none of the considered elements a true temperature dependence has been measured so far. This situation clearly reveals the necessity of highly accurate DFT calculations at finite temperatures as performed in the present study.
Our final data for the explicit SFEs including all finite temperature excitation mechanisms are shown in Fig. 14 by the thick solid lines. We observe a significant decrease of all SFEs with temperature, with the strongest decrease for Ni. At the melting temperature the SFE of Ni has a very small value of 10 mJ/m 2 and it becomes clear that the room temperature value of 110 mJ/m 2 (i.e., an order of magnitude higher) is not at all representative of high temperatures that are relevant for example for applications of Ni based super alloys. A proper consideration of the temperature dependence of the SFE is thus crucial. Although not available for the considered elements, a temperature dependent SFE has been experimentally derived for pure Ag [68] and the measured decrease of the SFE with temperature is consistent with our theoretical findings here for Al, Cu, and Ni.
Our results provide further insights: In Fig. 14 we have included additionally SFEs extracted from the SGTE unary thermodynamic database [64]. This database is the fundamental basis of the CALPHAD approach. These SFEs are based on the lattice stabilities for fcc and hcp and correspond thus to the first order ANNNI model. For Al we observe a very similar temperature dependence of the SFE as for our DFT calculations, although the absolute values from CALPHAD are by more than a factor of two higher. For Cu and Ni, the situation is very different and we observe a qualitatively different temperature dependence. Whereas our DFT calculations clearly reveal a decrease of the SFE with temperature, CALPHAD indicates an increase with temperature.
In order to exclude that the difference between DFT and CALPHAD is due to the comparison of the explicit SFE with the first order ANNNI model employed in the CALPHAD approach, Fig. 15 provides an additional comparison. Here, the CALPHAD curves (red lines) are the same as in Fig. 14. The solid black lines show now also for DFT the first order ANNNI results. We can draw a similar conclusion, i.e., that there is a qualitative difference in the temperature dependence of the SFE for Cu and Ni. Additional calculations also show that calculating the hcp phase at constant pressure has a negligible effect. Relaxing the c/a ratio can also not account for the observed difference between DFT and CALPHAD as it could only decrease the DFT SFE further downward.
It should be noted that, since there is no experimental data for the metastable hcp structures of these metals available, the corresponding CALPHAD based Gibbs energies are derived from extrapolations in different alloy systems. Our results reveal that, due to the subtle involved energy differences, such extrapolations may introduce qualitatively wrong temperature dependencies of the lattice stabilities and thus approximated SFEs. Preliminary results show that our DFT-derived lattice stabilities with the modified temperature dependence can be used to consistently parametrize Ni binary phase diagrams [70].

IV. SUMMARY AND CONCLUSION
We have conducted a highly-precise finite-temperature ab initio study of the temperature dependence of the stacking fault energy (SFE) of three prototype fcc metals Al, Cu, and Ni. All relevant temperature induced excitations have been considered, i.e., harmonic, quasiharmonic, explicitly anharmonic, electronic, as well as magnetic contributions for Ni. In all cases the SFE decreases with temperature. The largest Gibbs energy contribution to the SFE is captured by the quasiharmonic approximation. For Ni, electronic excitations are significant but are largely compensated by the coupling to explicit vibrations. Explicitly anharmonic as well as magnetic contributions (Ni) play a minor role. Due to the similarities of the magnetic properties of the various structures of Ni (fcc, hcp, dhcp as well as explicit stacking fault), the longitudinal spin fluctuations, being important for fcc Ni, are found to be largely compensated. Non-spin-polarized as well as ferromagnetic SFE calculations for Ni provide very similar temperature dependencies further supporting this finding.
The performance of the commonly employed ANNNI model in first and second order has been evaluated in comparison to explicit stacking fault calculations at elevated temperatures. Consistently with previous studies we find that the largest deviations occur at T = 0 K, whereas the absolute change in temperature is well captured by the ANNNI model. The impact of choosing a different approximation for the exchange-correlation functional (LDA or GGA-PBE) is similar, i.e., considerable at T = 0 K, but moderate for the temperature dependence.
We have evaluated the technical (supercell size convergence etc.) as well as conceptional treatments for practical SFE calculations. Our results show that the technical parameters such as the electronic convergence criterion have to be chosen extremely carefully to ensure converged SFE data. For the same level of accuracy, the vacuum approach requires by construction a larger number of atoms compared to the tilted periodic boundary supercell method and is therefore computationally less efficient. Self-consistent and perturbation approaches within the finite-displacement method are practically equivalent, but both techniques are sensitive to the employed technical convergence parameters. Long-ranged interactions, visible specifically in the quasiharmonic force constants of hcp Ni, are compensated by the explicitly anharmonic contribution that renders the interactions local, i.e., already a small supercell size (36 atoms) is well converged. Yet, we found that larger supercells may be necessary to prevent correlated displacements of close-packed atomic planes in the explicit anharmonic molecular dynamics simulations.
The present SFE predictions have been compared to a large set of previous experimental and theoretical data. The scatter within experimental as well as theoretical empirical potential based approaches is significant. No temperature dependence of the SFE has been measured so far for the investigated metals. Comparing the SFEs derived in this study with empirical CALPHAD derived data reveals qualitative differences for the temperature dependence in case of Cu and Ni, whereas the results for Al are in qualitative agreement.
We expect important consequences of our results for the field of Ni based superalloys. In designing such alloys, a key strategy is the reduction of the SFE by alloying elements (e.g., Co or Cr) in order to increase the resistance of dislocation movement [71]. Our calculations predict that already the SFE of pure Ni exhibits a strong intrinsic reduction from a room temperature value of 110 mJ/m 2 to only 10 mJ/m 2 for high temperatures.
Having now the ab initio based tools to accurately predict SFEs will allow us to accurately model the impact of alloying elements on the mechanical behavior at high temperatures. While in the present study we focused on three selected elements, the approach outlined here and the underlying concepts and convergence studies are general and can be thus applied not only to unary metals but also to alloys. We show here the relation between the pressure and temperature dependent explicit SFE, Eq. (1), and the approximate SFEs within the ANNNI model, Eqs. (5) and (6). The ANNNI approximations can be derived within the axial interaction model (AIM) [22,23]. Within the AIM the energy of closepacked structures is mapped onto an Ising-type Hamiltonian with interaction parameters (J 1 , J 2 ,...) describing the interaction between (111) layers. The distance between the layers is assumed equal in the different structures (fcc, hcp, dhcp, stacking fault) and in particular to be determined by the equilibrium geometry of the fcc phase. Thus, the SFE within the AIM model is expressed as (see, e.g., Ref. [69]) with the fcc equilibrium volume V fcc = V fcc (T ). The interface area is defined as before for Eqs. (5) and (6)  (A2) Depending up to which order the interaction parameters are included, one can obtain the ANNNI approximations Eqs. (5) and (6) (see, e.g., Ref. [69]): Second order, J 1 , In order to relate the constant volume expression Eq. (A1) to the constant pressure SFE [Eq. (1)], we start with  (3). The supercell size is given in terms of the conventional cubic unit cell for fcc (four atoms) and in terms of the primitive cell for hcp (two atoms), dhcp (four atoms), and stacking fault ["sf"; six atoms; cf. Fig. 1(a)]. "sf-ref" denotes the fcc reference supercell used to determine the explicit SFE, with the same geometry as the "sf" supercell but without the tilt to produce the stacking fault. The plane wave cutoff is given in eV and the lattice constants in Å. The notation for the lattice constants (a start ...a end , δa) means that the mesh is going from a start to a end in steps of δa. The four columns to the right of the table display results of parametrizing the E α (V ) curves using the Vinet equation of state [31]. a eq is the equilibrium lattice constant in Å, V eq the equilibrium volume in Å 3 , B eq the bulk modulus in GPa, and B eq the derivative of the bulk modulus with respect to pressure. For hcp and dhcp an ideal c/a ratio was assumed.   (3), on the static lattice (the coupling to vibrations was included during the anharmonic calculations, cf. Table IV). The supercell size is given in terms of the conventional cubic unit cell for fcc (four atoms) and in terms of the primitive cell for hcp (two atoms), dhcp (four atoms), and stacking fault ["sf"; six atoms; cf. Fig. 1(a)]. "sf-ref" denotes the fcc reference supercell used to determine the explicit SFE, with the same geometry as the "sf" supercell but without the tilt to produce the stacking fault. The plane wave cutoff is given in eV, the lattice constants in Å, and the temperatures in eV/k B (k B = Boltzmann constant). The notation for the lattice constants (a start ...a end , δa) means that the mesh is going from a start to a end in steps of δa. The same notation is used for the temperatures. The latter correspond specifically to the electronic temperature, i.e., to the Fermi-Dirac broadening employed in the finite temperature DFT calculations. where V sf and V fcc are implicitly defined by and where δV = V sf − V fcc . Now, we Taylor expand F sf (V fcc + δV , T ) up to first order in δV , with δP defined by Hence we see that the constant volume SFE γ (V fcc , T ) and thus the ANNNI approximations Eqs. (A3) and (A4) [or equivalently Eqs. (5) and (6)] differ from the constant pressure SFE γ (P , T ) only by smaller correction terms. Note in particular that, as long as V fcc (T ) is properly determined by Eq. (A7), this statement extends also to (arbitrary) finite pressures P . This is the reason why in Eqs. (5) and (6) γ 1 and γ 2 have been expressed with a dependence on P . Tables V to VII display a comprehensive list of previous experimental and theoretical (DFT and EAM) SFE data for Al, Cu, and Ni from the literature. The theoret-ical values refer to T = 0 K. For the experimental values the relevant temperature(s) is stated if available from the reference. The symbol τ 3 in the strain-rate dependence method refers to the stress at the onset of dynamical recovery in a tensile test. For further details on the experimental methods and a discussion on the credibility of experimental data we refer to Refs. [3][4][5]. The abbreviations used in the theory parts of the tables have the following meaning: APW=augmented plane wave, ASA=atomic sphere approximation, ASW=augmented spherical wave, EMTO=exact muffin tin orbitals, FP=full potential, LAPW=linearized augmented plane wave, LCAO=linear combination of atomic orbitals, LKKR=layer Korringer-Kohn-Rostoker, LMTO=linear muffin tin orbitals, MEAM=modified EAM, NC=norm conserving, PP=pseudopotential, PW91=Perdew-Wang 91, TB=tight binding, US=ultrasoft. IV. Parameters employed in the calculation of the anharmonic free energies, i.e., F ah α (V , T ) in Eq. (3), within the TU-TILD method. The supercell sizes are given in terms of the conventional cubic unit cell for fcc (four atoms) and in terms of the primitive cell for hcp (two atoms), dhcp (four atoms), and stacking fault ["sf"; six atoms; cf. Fig. 1(a)]. "sf-ref=fcc" means that the fcc supercell calculation was used as the reference for the explicit SFE. The two columns labeled "qh→potential" correspond to the thermodynamic integration from the quasiharmonic to the optimized EAM potentials. The following four columns correspond to the thermodynamic integration from the EAM potentials to DFT with low converged parameters, with the low plane wave cutoff and low k points denoted by E low and "kp low," respectively. The high converged cutoff and k points for the upsampling are denoted by E high and "kp high," respectively. The plane wave cutoffs are given in eV, the lattice constants in Å, and the temperatures in K.   Year Method