Ising-like models for stacking faults in a free electron metal

We propose an extension of the axial next nearest neighbour Ising (ANNNI) model to a general number of interactions between spins. We apply this to the calculation of stacking fault energies in magnesium—particularly challenging due to the long-ranged screening of the pseudopotential by the free electron gas. We employ both density functional theory (DFT) using highest possible precision, and generalized pseudopotential theory (GPT) in the form of an analytic, long ranged, oscillating pair potential. At the level of first neighbours, the Ising model is reasonably accurate, but higher order terms are required. In fact, our ‘ ANNNI model’ is slow to converge—an inevitable feature of the free electron-like electronic structure. In consequence, the convergence and internal consistency of the ANNNI model is problematic within the most precise implementation of DFT. The GPT shows the convergence and internal consistency of the DFT bandstructure approach with electron temperature, but does not lead to loss of precision. The GPT is as accurate as a full implementation of DFT but carries the additional benefit that damping of the oscillations in the ANNNI model parameters are achieved without entailing error in stacking fault energies. We trace this to the logarithmic singularity of the Lindhard function.

Ising-like models for stacking faults in a free electron metal Martina Ruffino, Guy C. G We propose an extension of the axial next nearest neighbour Ising (ANNNI) model to a general number of interactions between spins. We apply this to the calculation of stacking fault energies in magnesium-particularly challenging due to the longranged screening of the pseudopotential by the free electron gas. We employ both density functional theory (DFT) using highest possible precision, and generalized pseudopotential theory (GPT) in the form of an analytic, long ranged, oscillating pair potential. At the level of first neighbours, the Ising model is reasonably accurate, but higher order terms are required. In fact, our ' AN N NI model' is slow to converge-an inevitable feature of the free electron-like electronic structure. In consequence, the convergence and internal consistency of the AN N NI model is problematic within the most precise implementation of DFT. The GPT shows the convergence and internal consistency of the DFT bandstructure approach with electron temperature, but does not lead to loss of precision. The GPT is as accurate as a full implementation of DFT but carries the additional benefit that damping of the oscillations in the AN N NI model parameters are achieved without entailing error in stacking fault energies. We trace this to the logarithmic singularity of the Lindhard function.
2020 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited. 1

. Introduction
It is a 100 years since Wilhelm Lenz suggested a problem in magnetism to his student Ernst Ising, resulting in the now famous eponymous model [1] which has become a paradigm in statistical mechanics. Over 30 years ago, the Ising model and its extension to second neighbours, the axial next nearest neighbour Ising (ANNNI) model [2], found an application in the understanding of planar faults in semiconductors [3] and later in metals [4]. Moreover, in the study of semiconductors, the model was extended to third neighbours [5,6]. In addition to stacking fault energies (SFEs), the theory has been applied to antiphase boundaries [7,8]. In the simplest case of the comparison in stability of face centred cubic (fcc) versus hexagonal close-packed (hcp) lattices, fcc corresponds to the 'ferromagnetic' arrangement of spins, · · · ↑ · · · , and is stable in the Ising model for J > 0; and hcp is represented by 'antiferromagnetic' spin arrangement, · · · ↑↓ · · · , and is stable if J < 0. If the parameters of the ANNNI model differ in sign then the opportunity for competing interactions arises resulting in complexity in the phase diagram [2] and this has been suggested as an origin of the many high temperature modulated phases in binary alloys [9], first discovered in the same year as the Ising model by Johannson & Linde [10].
In this work, we study Ising-like models applied to the free electron metal magnesium. Planar faults in Mg are of great interest in problems of plasticity mediated by slip and twinning, effects of temperature and the role of alloying elements. It is curious that Mg alloys, having the lowest density among structural metals, are not as widely used as their low cost and high strength to weight ratio would suggest. One reason for this is difficulty in forging, as there is only a narrow temperature range in which this is possible; and their tendency during metallurgical processing to develop deleterious crystallographic textures. Present-day alloy design is targeted at formable, low cost alloys. From the viewpoint of fundamental theory, what singles out Mg in the present context is its electronic structure being free electron-like. This has the benefit that an analytic pair potential exists that emerges out of the density functional theory (DFT) [11,12]. On the other hand, this pair potential is long ranged as is to be expected from the logarithmic singularity in the Lindhard response function [13]. As a consequence, the natural question arises as to the validity of Ising and ANNNI models which are truncated after one or two neighbours. Our approach here is to propose an ' AN N NI model', which is an extension of the ANNNI model to interactions out to N neighbours. We develop this model with formulae for the four principal basal faults in hcp metals, up to N = 6. By making the most precise estimates of the parameters of the AN N NI model using electronic structure theory, we are able to assess the precision and the convergence of the Ising-like expansion.
The structure of this paper is as follows. We describe the basal planar faults peculiar to hcp metals in §2. Theory is presented in §3, with §3a on the Ising-like models; §3b on the total energy methods which we employ in this work. We present our results in §4, and Discussion and Conclusions will be found in § §5 and 6.

Planar faults in hcp metals
In hcp magnesium, there are four stable stacking faults in the basal plane, as described in Hirth & Lothe [14]. The first intrinsic stacking fault I 1 (or growth fault), is obtained through the removal of a basal plane, followed by the shear of the crystal above the fault of 1 3 a in the [1100] direction, where a is the hcp basal lattice constant (nearest neighbour distance). The second intrinsic stacking fault I 2 (or deformation fault) is formed by a shear of 1 3 a in the [1100] direction; it is this fault which separates partial dislocations of the a dislocation dissociated on the basal plane. The extrinsic stacking fault, E, is obtained by inserting an extra plane in the perfect sequence; the twin-like fault, T, possesses mirror symmetry about the fault plane [15] and is formed by replacing an A plane with a C plane [16] (see §3a). Experimental measurements of stacking fault energies are notoriously difficult and in Mg only upper bounds have been established. Hence, there has been a heavy reliance on first principles calculations based on DFT [17], with the main approach being the 'supercell' method [15,18]

Theory (a) Ising-like models
Many close-packed metals and alloys and many semiconductors based in the zincblende or wurtzite prototypes have crystal structures that correspond to the stacking of planes of closepacked spheres. Such hexagonal, or honeycomb, planes possess hollow sites labelled 'A' and 'B' which may receive a sphere from a subsequent plane, stacked atop the first, in such a way as to maintain close packing. If the receiving plane is designated 'C' then a close-packed repeating unit cell may be represented by a sequence of A, B and C designations which uniquely specify the crystal [19]. For example, · · · ABC · · · corresponds to the stacking of an fcc metal, while · · · AB · · · describes the ideal hcp crystal. In the case of the diamond lattice, the stacking is that in fcc, on the understanding that each close-packed 'sphere' is actually an atom plus a further basis atom. The scheme is applicable also to ordered alloys based on the fcc lattice and not necessarily restricted to the stacking of close-packed planes-for example, the L1 2 D0 22 and L1 0 superlattices [4,7,8].
If a crystal is constructed from a stacking sequence of atomic planes then its total energy may be written as the sum of the interaction energies between planes plus a contribution coming from the energy of the system when interactions are set to zero. A spin S = ±1 is assigned to the ith plane, according to whether the plane i+1 follows the 'correct' stacking sequence or not. Conventionally, we assert that the correct stacking sequence is A → B → C → A, whereas the incorrect one is A → C → B → A. Consequently, the energy of a crystal with an arbitrary stacking sequence is expressed as follows: with E 0 being the energy of the system when all interactions between planes are set to zero, M being the number of planes and N being the furthest neighbour interaction to be taken into account. The J n parameters are the interaction coefficients between neighbouring layers, such that e.g. J 1 represents the interaction energy between two layers which are nearest neighbours.
Here, we assume that the magnitude of the interaction energy between neighbours decreases as n increases. Truncated at the first term in the sum over n, equation (3.1) is the well-known Ising model [1]. In the ANNNI model [2], interactions up to J 2 are taken into account. In this work, we propose to observe the behaviour of the J parameters as further neighbours are added. We call these AN N NI models, with N = 1, 2, 3, 4, 5, 6 depending on how many neighbours are included. The sums in equation (3.1) can be evaluated for perfect structures straightforwardly; we have chosen seven structures, classified in terms of the type of unit cell as in the case of the polytypes of silicon carbide [5]. We label each structure with a number, representing the number of planes in the unit cell, followed by a letter, representing the Bravais lattice type, as in the Ramsdell classification [20]. The chosen polytypes and their stacking sequences are listed in table 1. It should be noted that the 2H structure is identical to hcp; while the 3R structure is identical to fcc, although it is no longer cubic if the material has non-ideal c/a ratio. In fact, the AN N NI model applied to an hcp metal need not employ the ideal axial ratio of c/a = √ 8/3, as long as the interplanar spacing of all the polytypes is adjusted to match that of the hcp, antiferromagnetic, ground state. In the present work, we use the axial ratio that is predicted by the electronic structure theory, see § §3bi and ii.
Evaluating the sums in equation (3.1) up to N = 6 and dividing by the number of planes M, we obtain the energies per plane of each polytype 2) 3)

I 1 fault
· · · ABABCBCB · · · · · · ↑↓↑↓↑↑↓↑↓↑↓ · · ·  Table 2. The energies of polytypes of pure Mg, calculated using the GPT; and LMTO-DFT with modified tetrahedron integration [25] and Fermi function smearing with kT = 0.03 Ry. The energy of the 2H structure is the binding energy relative to a free atom and is given in eV/atom, and the energies of the other structures are relative to that of 2H, i.e. given as differences in energy also in eV/atom. Note that the use of a large electron temperature leads to an absolute error in total energy of 0.12 eV/atom, while this systematic error clearly cancels in the energy differences between polytypes.
the structures containing the stacking faults, and thus eliminating one of the drawbacks of the supercell model. The SFE is then given as energy per unit area [3], such that γ = E/A with A being the area that defines the unit cell in one plane, A = √ 3/2 a 2 . SFEs calculated using this model are independent of cell size: each of the energies of the perfect crystals are calculated from the unit cell of each polytype, and the SFEs are then indirectly calculated as energies per plane, per atom. The cell size-dependency of the supercell method is thus avoided.

(b) Total energy
A key feature of the AN N NI model is that a large amount of information is extracted from a small number of total energy calculations on small unit cells. Specifically, calculations of the total energy of N + 1 polytypes are needed. Our data are displayed in table 2. It is vital that the greatest possible precision is used. In this work, we have used the DFT. We take two approaches to this. First, we use a standard procedure for obtaining solutions of the Kohn-Sham equations of DFT, here in a basis of linear muffin tin orbitals in a full potential implementation [22]; second, we exploit the fact that Mg is a free electron metal to apply the first principles pair potential formulation of the DFT total energy within the generalized pseudopotential theory (GPT) [11,12,23,24].

(i) LMTO-DFT
We use the full potential linear muffin tin orbital (LMTO) method of Methfessel and van Schilfgaarde [22] as implemented in the questaal code [26,27]. Kohn-Sham orbitals are expanded in a basis set of 'smooth Hankel functions' which are atom centred Hankel functions convoluted with Gaussians within a smoothing radius, R sm . The charge density is expanded in spherical and plane waves. For Mg, we use atomic spheres of radius 3 bohr, and a double energy basis set having Hankel energies of -0.1Ry in s, p and d channels and -0.9Ry in s and p channels. The smoothing radius is R sm = 2 bohr. The plane wave expansion cut-off is 8.2 Ry. Exchange and correlation are calculated using the GGA-PBE [28]. Brillouin zone integration is made using the modified tetrahedron method of Blöchl et al. [25]. To achieve the greatest precision, we are absolutely converged with respect to k-point sampling: typically dividing the Brillouin zone into 30 divisions along a reciprocal lattice vector of 2π/a. As an alternative, for reasons given below, for some calculations, we replace tetrahedron Brillouin zone integration with an effective electron temperature of kT = 0.03 Ry which acts to smear out the Fermi surface. The LMTO-DFT using tetrahedron integration predicts an equilibrium atomic volume, Ω 0 = 155.3 bohr 3 and an axial ratio of c/a = 1.627 and we use these in all LMTO calculations to avoid spurious stresses, especially in the supercells. The lattice constants are in good agreement with the measured lowtemperature values, Ω 0 = 153. 4   (ii) Generalized pseudopotential theory The GPT is also a DFT, but here the total energy is cast into a many-atom expansion comprising one, two and further body terms. The expansion is well known to be converged at the two-body term in the free electron metals [23,30]. This means that we can write the total energy as which is a function of the average atomic volume Ω and the bond lengths, R ij connecting atoms i and j. The pair potential v 2 is a function of Ω. We should note that since all our calculations in the present work are made at a fixed atomic volume, only the second term in the total energy (3.15) is involved in energy differences; as our structures are all closed packed there is no need for local volume corrections [24]. Detailed expressions for E vol and v 2 are given in [12]: v 2 (r) is proportional to where F is the energy-wavenumber characteristic [30,31]. Its oscillations in real space are evident from figure 1: the pair potential is long ranged and oscillating-a typical feature of a metal having a sharp discontinuity in the electron occupation number across the Fermi surface. The oscillations are asymptotically identical to those expected from the Lindhard screening function [13]. Exchange and correlation are treated in the screening of the perturbed free electron gas; and this allows the GPT to go beyond the local density approximation by using a local field correction which amounts to treating the electron-electron interaction at the level of the random phase approximation. In our case, we use an energy dependent non local pseudopotential for Mg in AHS form [32] within a small core, screened self consistently to first order using the dielectric function of Ichimaru & Utsumi [33] based on the correlation energy of Vosko et al. [34].
The GPT for Mg in this form has been thoroughly tested against experiment and standard DFT calculations [12,23,24]. For example, the GPT predicts that the atomic volume is Ω 0 = 160.9 bohr 3 and the axial ratio is c/a = 1. 621. In order to be comparing like with like, in our GPT calculations, we use the LMTO-DFT predicted lattice constants, see § §3a and 3bi.

(iii) Supercells
We build conventional supercells with the hcp c-axis normal to a rectangular area of basal plane. For Brillouin zone integrations in the LMTO calculations, we employ 60 divisions along the long axis of the Brillouin zone in the plane of the fault and 35 divisions along the short axis; there are six divisions along the reciprocal lattice vector normal to the basal plane and there are 12 (0002) planes separating the faults generated by the periodic boundary conditions. We use both modified tetrahedron integration [25], and smearing of the Fermi surface with a Fermi function at temperature equivalent to kT = 0.03 Ry: that is, an electron temperature of 4737 K. There is no Brillouin zone integration in the GPT; the total energy is simply given by summing equation (3.15) over all neighbours distant to R c = 28.3 bohr. The oscillating function v 2 is truncated smoothly to zero beyond R c by replacing v 2 (r) in the range R 0 < r < R c by a degree five Hermite interpolating polynomial designed to match continuously and twice differentiably to v(R 0 ) with R 0 = 26.7 bohr and zero at R c .

Results
As we mention in §3b, a rather modest number of total energy calculations (table 2) yields a very large body of data when entered into the formulae of the AN N NI model. This is because N + 1 total energies furnish us with the AN N NI model parameters J up to J N and these in turn lead to four stacking fault energies. The wealth of data extends beyond this since at any level of theory up to level five, we are free to choose which combinations of the seven total energies to use-that is to say, at level N we need to choose N + 1 total energies taken from the seven that we have available in table 2. This means that at level N we are free to make 7 N+1 combinations giving that many approximations to each of the four stacking faults at that level. In fact, a particular difficulty is how to display the data in a form that can be apprehended by the reader, so that analysis can be made and conclusions reached. Before taking that route, we display in table 3 the results from just the first neighbour, Ising ('ANNI') model and their comparisons to our supercell benchmarks and some from the literature. We will return to this table in the discussion, §5, but it is worthwhile to point out here that there is a considerable scatter in the published SFEs despite them all being calculated using standard DFT. This acts as an indication that such calculations are very demanding and require complete control over convergence with respect to supercell length, basis set and k-point mesh. To begin the presentation of the full data, we have assembled into table 4, in the appendix, a unique designation number, c, for each combination of polytypes that may be used to yield a unique estimate of each of the four-fault energies.   In addition to the GPT total energies, we have two total energies from the LMTO-DFT calculations. In figures 2-4, we show the AN N NI model J n parameters as functions of n plotted for all of our combinations from table 4. For the GPT total energies, the J-parameters for a particular n cluster quite closely together, while for the LMTO-DFT using the modified tetrahedron method there is considerable scatter. We find that this scatter can be largely removed at the cost of some precision by replacing the tetrahedon method in the Brillouin zone integration with a smearing of the Fermi cut-off with a Fermi function at temperature equivalent to kT = 0.03 Ry. We also observe that the J-parameters calculated using LMTO-DFT (figures 3-4) do not converge as rapidly towards zero as do the J-parameters calculated using GPT (figure 2). We will return to these observations in §5b. It is striking that the very evident contrast in the J-parameters using tetrahedron Brillouin zone integration or an electron temperature is concealed in the very small differences in polytype energies displayed in rows 2 and 3 in table 2. Although there seem to be very minor deviations in the energy differences in table 2 these clearly conspire to have dramatic consequences when equations (3.2)-(3.8) are inverted.
A third way to display the data is presented in figures 9-20 in the Appendix. Here, we show as ordinates each estimate of each fault energy as calculated from each of the 94 combinations from table 4 as abscissa to each plot. It will be possible for the reader to identify each calculation and in particular its excursion from the benchmark supercell SFE with one of the 94 combinations. In each of figures 9-20, we show as horizontal lines fault energies calculated using the supercell royalsocietypublishing.org/journal/rspa Proc. R     method described in [42]. There is evidently considerable scatter in these data and convergence is not convincingly demonstrated. This is to be expected in a free electron metal as will be discussed in §5. On the other hand much of the scatter, which is also evident in figures 2-4 can be traced to certain 'outliers' in figures 9-20. As we explain in the appendix, we can eliminate much of the noise by excluding results from those combinations in  longest period polytypes, 9R, 10H and 15R. We then plot the predicted SFEs for each fault and for each of the three total energy methods in figures 5-8. We defer discussion of these results to §5c.

Discussion
(a) Benefits of the AN N NI model, and the Ising model first neighbour approximation As far as we are aware, there are few or even no direct measurements of stacking fault energies in Mg. Normally only the intrinsic fault energy can be measured and we are aware of one experimental measure that places the energy of the I 2 fault at lower than 50 mJ m −2 [43]. This argues strongly for the need of theory to approach the question. Traditionally fault energies have been calculated using the supercell method in DFT and whereas this was expensive in the past, leading to the application of the Ising and ANNNI models to the problem, with current computers it is not particularly challenging to calculate fault energies (notwithstanding the considerable uncertainty in the literature-table 3). On the other hand, the approach we take here has a number of advantages. 1. Clear insight is obtained into the relation between structural energy differences between polytypes and stacking fault and twin boundary energies. 2. The method is in principle extensible to alloys, for which the supercell method becomes cumbersome. In addition, there is the promise of much greater insight since the theory of alloy stability and the effects of impurities on the electronic structure of metals may be transposed into an understanding of the role of alloying and impurities on fault energies and through that understanding to the plastic response of an engineering alloy. This is a matter for further work. 3. A possible disadvantage is that whereas using the AN N NI model does not require us to work with ideal axial ratio hcp metals, any relaxation of the interplanar spacing or atomic positions on creating the defect are neglected. In the present case, however, these are not limitations. The intrinsic stacking fault energy calculated with the supercell method and neglecting relaxation is 35 mJ m −2 (table 3), whereas using the same DFT program and lattice constants but allowing relaxation the result is 32.1 mJ m −2 [42]. 4. As an example of insight gained, it is clear from equations (3.11) to (3.14) that at the level of nearest neighbours we have the relation, and indeed this is reasonably well adhered to in the supercell calculations in table 3, at least in respect of I 1 and I 2 . We also have at the Ising (ANNI) level that γ I 2 = γ T . These relations expose the well known first approximation to the fault energies as proportional to energy differences between fcc and hcp crystals [3][4][5]18].

(b) Convergence of the J-parameters
We turn next to a discussion of the series (3.1) taken beyond second neighbours. It is rather clear from table 3 that the first neighbour Ising model leads to a faithful rendering of the four fault energies compared to supercell calculations, particularly if the GPT is employed in the total energy calculation. However, the ANNI relation γ I 2 = γ T is not supported by the supercell results. This is to be anticipated in view of the long range of atomic interactions expected in a free electron metal (figure 1), and it is important to examine the convergence of further neighbour interactions.
We are not aware from the literature of any study which makes such a thorough examination of the convergence of what we call the AN N NI model as we have done here. First, we have examined the convergence of the J-parameters; and second, we have considered every possible combination of up to seven polytypes in order to check the method for internal consistency. Our full sets of data are presented in the appendix. Whereas figures 2-4 and table 4 show that the J n for n > 1 are much smaller in magnitude than J 1 , the subsequent convergence is rather slow, as is to be expected in a free electron metal. One reason for this is that oscillations in the J n as functions of n are not sufficiently damped (table 4), especially in the case of the LMTO with tetrahedron integration; if J N has not approached zero at the point that we truncate the AN N NI model then neglect of the unknown slowly decaying J n , n > N will lead to the inconsistency across different combinations and scatter in the SFEs that we observe. A second reason is outside the control of the precision in total energy calculations and arises from the fact that the prefactors to the J-parameters do not get smaller with each order of approximation, both in the expression for the polytype energies (3.2)-(3.8) and for the fault energies (3.11)-(3.14). On the other hand, there is a very clearly increased scatter in figure 3 in the case of the LMTO-DFT with tetrahedron integration compared to the same parameters calculated using the GPT (figure 2). This scatter is removed in the LMTO-DFT if we replace the most precise Brillouin zone integration with the approximation of applying an electron temperature (figure 4). What the latter does, most significantly, is to change the abrupt discontinuity of occupancy of the one-electron eigenstates into a smooth transition. This is exactly equivalent to removing the logarithmic singularity in the energy-wavenumber characteristic in the theory of simple metals. The question then arises, whether and how this is evidently achieved in the GPT. One answer is that we implicitly smooth the Fermi surface discontinuity by applying a cut-off to the pair potential. This is indicated in figure 1 as two vertical lines: the pair potential is not taken to an infinite number of neighbours in the computation but is smoothly interpolated to zero using a polynomial which replaces v 2 between two distances which we choose to be 26.7 bohr and 28.3 bohr [24]. While applying an appropriate electron temperature in the LMTO-DFT has the desirable effect of removing the excessive scatter in the J-parameters there is a price to pay, namely  a loss of overall precision in the solution to the Kohn-Sham equations. This is reflected in table 3 which shows that the supercell fault energies using the finite electron temperature are not as well in accord with the most precise LMTO-DFT values which in turn are in very good agreement with the GPT. The conclusion must be that the GPT is not only as accurate as the standard DFT in predicting properties of Mg, it manages to deal with the screening of the free electron gas better than the standard DFT for which either at zero temperature there is scatter in the J n or at finite temperature there is a loss of precision. We cannot be absolutely certain from where this benefit of the GPT over the DFT bandstructure method arises; however, a comparison with the analytical pair potential of Pettifor and Ward is instructive [44]. These authors fitted a rational approximation to the Lindhard function so that their derived energy-wavenumber characteristic had no singularity at twice the Fermi wavevector. This effectively damps the pair potential, but unlike our v 2 which is only damped beyond R 0 , the Pettifor-Ward pair potential is damped at all r: the peaks are lowered and the troughs are filled in [30]. This circumstance, in our opinion, corresponds to the case of temperature smearing of the discontinuity in occupation number in our LMTO-DFT calculations: in that case all states in reciprocal space that are within kT of the Fermi surface are affected by the electron temperature, just as all interactions in real space are affected by the Pettifor-Ward sextic interpolation replacing the Lindhard screening function.

(c) Convergence of the fault energies
As long as we sensibly remove outliers (see appendix) then we clearly see the convergence of the J-parameters, as discussed in §5b, reflected in the mean and standard deviation of the stacking fault energies. Figures 5-8 show the SFEs calculated using GPT as functions of N in the AN N NI model. Convergence is rather clearly demonstrated and there is a useful, if not dramatic, benefit of increasing N, when compared to the supercell result. The diametrically opposite situation is found for the LMTO-DFT using tetrahedron integration. Our conclusions of §5b concerning the weak damping of the oscillations in J n , shown in figure 3, are confirmed in the large standard deviations and indeed negative predicted mean SFEs in the cases of the growth, I 2 , and twin-like faults. As expected in view of our discussion in §5b these issues are largely, but not fully, resolved with the use of a large electron temperature. The plots of SFEs calculated using LMTO-DFT and the electron temperature in figures 5-8 now show by and large a narrowing of the error bars as N is increased, and in the cases of the deformation and extrinsic faults significant convergence to the supercell result. But it must be recalled that the supercell result itself suffers from the error in Brillouin zone integration compared to the tetrahedron method.
It is important to point out that in the cases of the GPT and LMTO-DFT with tetrahedron integration the outcome of the N = 6 AN N NI model is extremely close to the supercell result for all stacking fault energies, see figures 9-16. This argues strongly for the correctness of the theory.

Conclusion
We have been motivated to calculate to the best available precision the energies of the four basal stacking faults in the hexagonal close packed, free electron metal, magnesium. While these are critical control parameters in the development of microstructure and mechanical properties, they are almost totally inaccessible to experiment so that recourse to theory is mandated.
A thorough investigation was made of the precision and convergence of the ' AN N NI model': an Ising model describing the total energy of crystals that may be described by a packing sequence  of close-packed planes, extended to a range of inter-plane interactions up to the sixth neighbour. By a suitable limiting process, the parameters obtained from the study of polytypes may be used to obtain the energies of isolated stacking faults [3,5]. We have chosen a particularly challenging case, that of the free electron metal magnesium, whose interatomic forces are known to be of a particularly long range due the prominence of the discontinuity of electron occupancy across the  nearly free electron Fermi surface. This results in a weak logarithmic singularity in the energywavenumber characteristic in reciprocal space: in turn leading to very weak damping of the screening of the pseudopotential (Friedel oscillation) in real space. These are well captured within the GPT by an analytic pair potential, v 2 , figure 1. Its range is particularly extended in the free electron metals, compared to the range of v 2 in transition metals [45].
A feature of the AN N NI model is that a small number of high precision total energy calculations of N + 1 polytypes provides the first N J-parameters of the model. In addition, at each level of the expansion, we have a choice of up to 7 N+1 combinations of polytype each in principle leading to equal values for the J-parameters. The deviation from this equality is a measure of the internal consistency of the model. This deviation was illustrated by points of the same colour in figures 9-20. These figures illustrated the cases of the three choices we made for finding the total energies of the polytypes, table 2, namely the GPT and the DFT as implemented in the full potential LMTO method [26]. We have validated our AN N NI model results against conventional 'supercell' calculations of fault energies. Our conclusions are as follows.
-Our benchmark calculations are the four fault energies calculated using supercells in the LMTO implementation of DFT [26] -As is evident from table 3, already at the level of first neighbours-the Ising model-fault energies are in reasonable agreement with supercell calculations using the same level of theory. In addition, relations between the fault energies within the ANNI model (5.1) are reasonably well adhered to in the supercell calculations (better using the GPT than the LMTO-DFT), although the extrinsic fault energy is significantly smaller in benchmark than predicted by equation (5.1). The Ising (ANNI) model prediction γ I 2 = γ T is not supported by the benchmark calculations. This means that corrections to the Ising model are desirable. -If such corrections are to be made then as seen in the remaining data of this paper, care must be taken. In our most precise calculations, the LMTO-DFT with tetrahedron integration, we find that the AN N NI model expansion up to N = 6 results in J-parameters that are oscillating and insufficiently close to zero to regard the expansion (3.1) to be converged ( figure 3). This is a principal conclusion of the this paper: the AN N NI model for a free electron metal is slow to converge. The symptom of this slow convergence is scatter in the values of the J-parameters obtained by using different choices of polytype combinations; scatter in the predicted fault energies, figures 13-16; and also in an error in the intrinsic fault energy in the Ising model, with large standard deviation as seen in column 4, table 3. We attribute this poor convergence to the sharp occupancy cut-off at the Fermi surface and confirm this supposition by making identical LMTO-DFT total energy calculations but with a finite electron temperature designed to smooth the zero temperature Heaviside step in the Fermi function. Convergence is thereby restoredfigure 4. However, a heavy price is to pay in the consequent loss of precision as seen in column 7, table 3, which shows that such a large smearing in the Brillouin zone integrals leads to errors in the supercell results for the fault energies, when compared to the tetrahedron method. -We find that the GPT offers good convergence of the J-parameters figure 2, with scatter in the fault energies consistent with the high electron temperature in LMTO-DFT, figures 9-12, but without the attendant loss of precision as is evident from the supercell results in table 3. This is a key conclusion of this paper. The damping of the logarithmic singularity in the energy-wavenumber characteristic is effected not in reciprocal space [44], but in real space [24] by a cut-off in the long-ranged part of the pair potential, figure 1; this does not compromise the precision of the total energy expression (3.15).
A matter for further work is the question of whether the inclusion of three-spin or even four or more spin terms in the AN N NI model may improve either precision or convergence. Our surmise is that it will not and that such an extension to the theory takes us beyond the simple context of the Ising description. Indeed, as we point out at the end of §5c, our longest ranged, N = 6, AN N NI model is in excellent agreement with our benchmark calculations, both the GPT and LMTO-DFT, in respect of all four stacking fault energies. Inclusion of three-spin terms would allow us to test the hypothesis that long-ranged interactions in the pairwise theory may be accounted for by first neighbour many-spin terms, and this remains a matter for future work.
Data accessibility. All input files used for the calculation of the total energies of the polytypes of magnesium as well as supercells can be found at a public GitHub repository http://doi.org/10.5281/zenodo. 4059479. The questaal suite, including the LMTO-DFT code, is available at https://www.questaal.org. The GPT code is available as part of the LAMMPS code at https://lammps.sandia.gov.