Physical factors controlling the observed high-strength precipitate morphology in Mg–rare earth alloys

Available online at www.sciencedirect.com

ScienceDirect Acta Materialia 65 (2014) 240–250 www.elsevier.com/locate/actamat

Physical factors controlling the observed high-strength precipitate morphology in Mg–rare earth alloys A. Issa ⇑, J.E. Saal, C. Wolverton Department of Materials Science and Engineering, Northwestern University, 2220, Campus Drive, Evanston, IL 60208-3108, USA Received 25 June 2013; received in revised form 29 October 2013; accepted 29 October 2013 Available online 21 November 2013

Abstract In an effort to understand the exceptional precipitation strength in Mg–RE (RE = rare earth) alloys, we use first-principles density functional theory calculations to study the energetic stability, elastic constants and coherency strain energy of Mg3 RE-D019 precipitate phases in a Mg matrix and make extensive comparisons with experimental and theoretical work, where available. We find the metastable b00 -D019 phases are energetically competitive with the stable Mg-rich phases for all RE elements. We also investigate the coherency strain energy of Mg–Mg3 RE binary systems using first-principles methods and harmonic elasticity theory. We find the two approaches to computing coherency strain energy in good agreement, indicating the validity of using harmonic elasticity equations, which we extend to hexagonal systems, to study the direction-dependent coherency strain energy of D019 precipitates in Mg–RE binary systems. From our coherency strain calculations, we find the D019 precipitates to strongly prefer prismatic, as opposed to basal, habit planes for all Mg–RE systems. This work thus provides an explanation for the observed prismatic plate-shaped morphology of many Mg–RE precipitates, which is ultimately responsible for their strengthening response. Ó 2013 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Magnesium alloys; Rare earths; Precipitate morphology; Density functional theory; Coherency strain

1. Introduction Magnesium alloys are attractive for structural applications because of their lightweight properties, making them suitable for aerospace and automotive applications. A key deficiency of Mg alloys, namely their low critical resolved shear stress (scrss ) in the basal plane [1], has inhibited them from gaining more widespread use. Increasing basal scrss may enable more widespread use of Mg alloys in structural applications. Precipitation hardening is one avenue for strengthening, where an element, or a combination of elements, form precipitates which impede dislocation movement in the basal plane, increasing basal scrss . Precipitation strengthening of Mg alloys has, therefore,

⇑ Corresponding author.

E-mail addresses: [email protected] (A. Issa), [email protected] (C. Wolverton).

increasingly been the subject of intense recent research (e.g. [2–18]). The extent of strengthening from precipitate formation is determined by several factors: the quantity of precipitates, their shape and size, their orientation with respect to the matrix, the degree of coherency with the matrix, and, particularly for coherent precipitates, their susceptibility to shear. Precipitates resistant to shearing improve alloy strength through the creation of Orowan loops as dislocations pass around the precipitate. Shearable precipitates can also improve alloy strength through the formation of hardening defects, such as anti-phase boundaries. Nie examined the effect of precipitate shape and orientation on strengthening in hexagonal close-packed (hcp) alloys [19] and showed that, for a given volume fraction of second-phase precipitates, the interparticle spacing in the basal plane is minimized by the formation of prismaticplate-shaped precipitates. One might therefore expect precipitates extended along prismatic planes to offer more

1359-6454/$36.00 Ó 2013 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.actamat.2013.10.066

A. Issa et al. / Acta Materialia 65 (2014) 240–250

strengthening, as opposed to precipitates extended along basal, pyramidal or other planes, regardless of whether the particle is deformable. Mg alloys with rare earth (RE) element additions [10– 18], in particular, have been the subject of numerous investigations focusing on their ability to provide precipitation hardening. Experimental investigations [15,18,20–25] have observed the presence of plate-like b00 precipitates that have prismatic habit planes, generally reported as f1010gMg [15,21,22]. Mordike et al. [25,24] have observed two distinct aging reactions, depending on the alloying elements. For Mg–heavy RE alloys (heavy RE = Tb, Dy, Ho, Er, Tm and Lu), the first phase to precipitate out of Mg solid solution is b0 , which has an orthorhombic structure. For Mg– (Ce, Nd, Gd) and WE (W = Y, E = Nd + a heavy RE element) alloys, a coherent hexagonal phase, thought to be in the D019 structure, called b00 , is the first phase to precipitate out of the solid solution, and, upon further aging, an orthorhombic phase, called b0 , forms with the {1 2 1 0} habit plane [26]. The composition of the b00 -D019 phase is Mg3 RE, while the composition of both b0 -orthorhombic phases is thought to be Mg7 RE [27]. While Mordike et al. [25,24] did not observe the formation of a b00 -D019 phase in Mg–HRE alloys, Rokhlin et al. found clusters of D019 in the Mg–(Y, Gd, Tb, Dy, Er) systems using transmission electron microscopy [28–31,18]. The hardness increase in these Mg–RE binary alloys is generally attributed to the formation of the metastable b00 and b0 precipitates during aging (see, for example, Ref. [4], Nie’s extensive review [32] and references therein). The emerging consensus, then, is that there are two primary phases in Mg–RE alloys, b00 and b0 , that increase the alloy strength. Both b00 and b0 appear in Mg–light RE alloys, while for Mg–heavy RE alloys only b0 is found. b00 is invariably observed to be fully coherent with the Mg matrix. On the other hand, b0 has been reported in some experimental work to be fully coherent [26], and in other work only semicoherent [33]. To achieve the goal of significantly strengthened Mg alloys, these hardening precipitates warrant an investigation into the underlying factors that give rise to their desirable mechanical properties. First-principles calculations lend themselves well to such an effort [34–42] for various reasons. First, the number of possible alloying elements and the composition ranges makes a comprehensive experimental study difficult. Second, while past common practices of alloying Mg with mischmetal (a combination of La, Ce, Nd and Pr) are known to produce a stronger alloy, the isolated element-by-element effect of rare earths on the properties of Mg alloys is not precisely known. First-principles methods are well suited to providing an understanding of each alloying element in isolation and in combination with other alloying elements. Third, first-principles calculations readily provide information that may be difficult to measure experimentally. For example, thermodynamic and elastic properties of metastable precipitate phases are readily calculated from first principles and are important for any study of precipitation in Mg–RE

241

systems. However, these quantities of metastable phases can be difficult to measure. Also, while our work here will only address Mg–RE binary systems, these calculations form a foundation from which to explore ternary or more complicated alloying schemes. Using first-principles calculations, we can explore equilibrium precipitate morphology which is, in general, controlled by a competition between two energetic contributions: interfacial energy between precipitate and matrix and the coherency strain energy of forming a coherent interface. Coherency strain calculations measure the energy penalty of elastically deforming the precipitate and matrix to form a coherent interface. While interfacial energy scales with the area of the interface, coherency strain energy scales with the volume of the precipitate, thus dominating as the size of the precipitate increases. In this work, we will focus solely on the coherency strain properties of these precipitates, with the understanding that a complete exposition of the morphology (e.g. determining the precipitate’s aspect ratio) will require knowing the interfacial energy contributions, which we will explore in a future paper. The main objective of this paper is to investigate, via density functional theory (DFT) calculations, whether we can explain the appearance and morphology, or at least habit plane, of D019 -Mg3 RE precipitates in Mg alloys and, if so, compute the preferred habit plane across all Mg–RE systems. We also calculate the formation energies and elastic constants of the D019 precipitates, which we use to assess their energetic and mechanical stability. This study will then provide a first step towards understanding the role of precipitate morphology in high-strength Mg alloys, and ultimately will accelerate the replacement of RE elements in these alloys, a long-term goal in the Mg alloy research effort. 2. Methodology 2.1. DFT calculations DFT calculations were performed using the Vienna Ab initio Software Package (VASP) [43–45]. We utilize the projector augmented wave method with the generalized gradient approximation [46,47] using the PBE exchange– correlation functional [48]. For all the D019 calculations, an energy cutoff of 350 eV and a gamma-centered k-point mesh of 10  10  13 were used. For the pure constituent elements in their ground state structure, the same energy cutoff was used, along with a k-point mesh size such that 10,000 k-points per reciprocal atom was approximately achieved. Tests of different k-point meshes indicate that our calculated total energies are converged to within 1 meV/atom. The f-electrons in the lanthanide RE elements exhibit unusual mixed valence/core characteristics, making a DFT description problematic [49–51]. The f-electron energies lie near the Fermi level, suggestive of valence

242

A. Issa et al. / Acta Materialia 65 (2014) 240–250

electrons, but these electrons are highly localized in space around the nucleus, thus making them somewhat core-like. These characteristics give rise to complicated magnetic ordering and make the RE metals difficult to model using either the local density approximation or the generalized gradient approximation. Whether the f-electrons are treated as valence or core states depends on the physical property of interest. Previous studies of ground state structures, formation energies and elastic constants of RE elements and compounds [52–60] indicate the f-core approach is the correct method to treat the RE elements when thermodynamic and elastic properties are of interest. In this work, pseudopotentials with the f-electrons frozen in the core were used throughout. The published Mg–RE phase diagrams [61–70] indicate the formation of a stable Mg3 RE in the BiF3 structure for the Mg–(La, Ce, Pr, Nd, Sm, Gd, Tb, Dy) binary systems and the formation of two-phase regions for Mg–(Eu, Ho, Er, Tm, Yb, Lu) systems. As mentioned earlier, the D019 -b00 phase also has a composition of Mg3 RE. The metastability of the D019 -b00 phase and the stability of the BiF3 structure indicate an energetic competition between these two structures at Mg3 RE stoichiometry and give us an opportunity to test the accuracy of our DFT approach. For all of the systems where BiF3 is observed, the energy of this phase must be lower than that of D019. Using DFT methods, we can calculate the formation energy of both phases, shedding light on the relative stability of the b00 phase for each Mg–RE system. The formation energy is calculated using the equation:   3 1 DH f ðMg3 REÞ ¼ EðMg3 REÞ  EðMgÞ þ EðREÞ ; ð1Þ 4 4 where E (X) is the total energy of structure X per atom, and DEf is the formation energy of a given structure. 2.2. Coherency strain calculations The coherency strain energy, DECS (^k; x), is defined as the strain energy required to maintain a coherent interface along direction ^k between bulk Mg and Mg3 RE precipitate. This quantity is dependent on both direction (^k) and composition (x). In this work, x is the mole fraction of Mg3 RE phase, making x ¼ 1 “pure” Mg3 RE. In cubic systems, the coherency strain energy is calculated from the energy change of the biaxial deformation of constituent phases A and B from their bulk lattice constants to a common in-plane lattice parameter (as ) and relaxing in the third direction [71]: epi ^ epi ^ ^ DEcubic CS ðk; xÞ ¼ min½ð1  xÞDE A ðk; as Þ þ xDE B ðk; as Þ; as

ð2Þ

epi where DEepi A and DEB are the energies required to deform material A and material B biaxially, respectively. In hexagonal systems, Eq. (2) is modified because there are two independent lattice constants to minimize over:

epi ^ epi ^ ^ DEhex CS ðk;xÞ ¼ min ½ð1  xÞDE Mg ðk; as1 ; as2 Þ þ xDE Mg3 RE ðk; as1 ;as2 Þ: as1 ;as2

ð3Þ

For the prismatic plane [1 0 1 0], the two in-plane substrate lattice constants are a and c, while for the basal plane [0 0 0 1], the in-plane lattice constants are identical: a. Consequently, for the high-symmetry case of the basal plane, Eq. (3) reduces to Eq. (2). Using DFT, we calculate the coherency strain energy, DEhex CS , as a function of composition (x) along basal [0 0 0 1] and prismatic [1 0 1 0] planes for Mg–Mg3 RE systems. For the basal plane coherency strain calculations, we use the Automated Theoretic Alloy Toolkit (ATAT) to set up the DFT calculations [71–73]. For the prismatic plane coherency strain calculations, a two-dimensional minimization over as and cs is necessary. We note here that these direct-DFT calculations capture the complete anharmonic elastic response of the constituent phases, in contrast to using harmonic elasticity along with elastic constants to calculate the coherency strain energies. We use both approaches to elucidate the nature, harmonic vs. anharmonic, of the elastic response of the precipitates. 2.3. Bulk moduli and elastic constant calculations The bulk moduli for the pure REs and D019 were calculated by performing a series of fixed-volume calculations, slightly varying the volume around the ground state value, V 0 . The energies and volumes are then fit to the fourparameter Birch–Murnaghan equation of state [74–76]: 8" #3  2=3 9V 0 B0 < V 0 EðV Þ ¼ E0 þ  1 B00 16 : V 9 "  #2 "  2=3 #= 2=3 V0 V0 þ 1 64 ; ð4Þ ; V V where E0 is the energy at the ground state, B0 is the bulk modulus at zero pressure, a0 is the equilibrium lattice constant and B00 is the derivative of the bulk modulus with respect to pressure. The four parameters in the fit are E0 ; B0 ; V 0 and B00 . We also compute the elastic constants, C ij , for RE metals in their ground states and the D019 structures by applying various deformations on the primitive cell, and extracting the elastic constants [77–80]. Since the D019 structures and most of the REs are of hexagonal symmetry, five independent elastic constants (C 11 ; C 12 ; C 13 ; C 33 and C 44 ) are needed to completely describe the elastic properties. Four deformation matrices are applied to the hexagonal unit cell to calculate the elastic constants. The lattice vectors for the primitive hcp unit cell are: pffiffi 01 1 a  23 a 0 2 pffiffi B1 C 3 ð5Þ @ a a 0 A: 2

0

2

0

c

A. Issa et al. / Acta Materialia 65 (2014) 240–250

The first 0 d B  ¼ @0 0

deformation to apply to the unit cell is: 1 0 0 C d 0 A; 0 0

ð6Þ

where d denotes a small deformation from the equilibrium cell vectors. The change in energy due to the deformation, neglecting higher than second-order terms, in d, is: DEðdÞ ¼ V ðC 11 þ C 12 Þd2 : The second deformation matrix is: 0 qffiffiffiffiffiffi 1 1þd 1 0 0 1d B C qffiffiffiffiffiffi B C ¼B 1d 0  1 0C @ A 1þd

ð7Þ

ð8Þ

0 0 0 and the associated energy change is: DEðdÞ ¼ V ðC 11  C 12 Þd2 :

ð9Þ

From the calculations of DE (d) along these two deformation paths (Eqs. (6)–(9)), we can compute C 11 and C 12 . To calculate C 33 , we use a third distortion of the form: 0 1 0 0 0 B C  ¼ @ 0 0 0 A; ð10Þ 0 0

d

which results in the equation: 1 DEðdÞ ¼ V  C 33 d2 : 2 The final 0 0 B  ¼ @0 0

ð11Þ

distortion matrix gives the C 44 elastic constant: 1 0 0 0 d C ð12Þ A; 1 2 d 4d

which gives the equation: 1 DEðdÞ ¼ V  C 44 d2 : ð13Þ 2 Since the bulk modulus is known along with four of the five elastic constants, the last elastic constant, c13 , can be determined through the relation: 2 B ¼ ðC 11 þ C 12 þ 2C 13 þ C 33 =2Þ: 9

243

various competing structures with the Mg3 RE composition: (1) the observed stable phases—BiF3 for the first 10 lanthanides (La–Dy), Mg24 RE5 –Mg2 RE two-phase regions for Mg–(Ho, Er, Tm and Lu) binary systems, and the Mg– Mg2 Yb two-phase region for the Mg–Yb system [61–70]; (ii) the metastable D019 ; and (iii) Mg3 RE hcp solid solution. Formation energies of all phases are shown in Fig. 1 for all the Mg–RE systems. Based on the observed phase diagrams described above we would expect BiF3 to be energetically more stable than D019 , which would in turn be more stable than Mg3 RE solid solution. Our DFT calculations show the correct phase stability for most Mg– RE systems. Specifically, we find BiF3 -type Mg3 RE compounds are lower in energy than D019 for La–Dy, in agreement with experiments. Both structures have negative formation energies and are fairly competitive. For Mg– (Er, Tm, Yb and Lu) systems, we also correctly predict the experimental two-phase region to be more stable than the D019 structure. For the Mg–Ho and Mg–Y systems, we predict D019 to be more stable than the two-phase region. This prediction agrees with previous DFT work [58], indicating either a difficulty in correctly describing those two systems from first principles, or perhaps that some entropic contribution, not considered here, changes the stability of these phases at finite temperatures. Out of the 17 systems calculated, our DFT calculations show the expected energetic order of stability at the Mg3 RE composition for 15 systems. Because the lattice mismatch between the Mg matrix and D019 precipitates is key to determining the magnitude of the coherency strain between the precipitates and the Mg matrix, we are interested in the volumes and lattice parameters of the metastable D019 phases. We are also interested in how well the DFT calculations can predict the volumes of the pure RE elements in comparison to experiment. It is

ð14Þ

For each deformation path, we used a mesh of five points including the equilibrium volume with a d spacing equal to approximately 1% of the lattice parameter. From our DFT calculations of the elastic constants, we examine the mechanical stability of the Mg3 RE D019 compounds at the DFT-relaxed lattice parameters and at the constrained lattice parameters of Mg. 3. Results and discussion 3.1. Energetic stability and volumes of D019 compounds We begin by assessing the phase stability of the Mg3 RE structures. We compute the formation energy, DEf , of

Fig. 1. Formation energies of D019 , experimentally stable phases [61–70] and solid solution. The experimentally stable phases are BiF3 for Mg–(La– Dy) systems, Mg24 RE5 –Mg2 RE two-phase regions for Mg–(Ho, Er, Tm, and Lu) and Mg–Mg2 Yb for the Mg–Yb system. The D019 phases are energetically competitive with the stable BiF3 phase. The DFT calculations correctly predict BiF3 to be the most stable stable, followed by the D019 , and finally the solid solution, as expected, validating our DFT approach.

244

A. Issa et al. / Acta Materialia 65 (2014) 240–250 Volume of RE and Mg3RE (D019)

Pure RE Bulk Modulus

50 70

(a)

DFT Expt.

60

Bulk modulus (GPa)

Volume (Å3/atom)

45

D019 Pure RE DFT Pure RE Expt.

40

35

30

25

50 40 30 20 10

20 Sc Y

0

La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

Sc Y

Fig. 2. Volume of D019 phases and pure REs from DFT and experiment [82]. The RE elements and compounds tend to decrease in size as one goes down the lanthanide series, with the exception of Eu and Yb, the larger size of which is due to their electronic structure.

3.2. Elastic properties The calculated bulk moduli for the pure REs and D019 structures are shown in Fig. 3. The bulk moduli of the pure REs are also compared to experiment [82], and show good agreement. An exception is pure Ce, a discrepancy that has been discussed in other work [53,83] and is thought to originate from the unusual magnetic properties of this element. The observed trend [82,81] of increasing stiffness with increasing number of f-electrons is also well reproduced in the calculations for pure REs. Eu and Yb are again the two notable exceptions to the trend, both theoretically and experimentally. The trend of increasing stiffness with increasing number of f-electrons is also present for D019 , but to a lesser degree than the pure REs. The weak trend of increasing stiffness in D019 structures is again interrupted with Mg3 Eu and Mg3 Yb being significantly softer than the rest of the D019 compounds. The elastic constants for the pure RE are given in Table 1 and compared to previously calculated values using ”f-core” potentials in VASP with the projector augmented wave (PAW) method within the GGA-PW91 parameterization [80,84] and experiment [85,86]. Two notable characteristics are immediately visible: the lanthanides get stiffer with increasing atomic number, and Eu and Yb are significantly softer than the neighboring RE metals, just as for the bulk moduli. The elastic constants for the

Mg3RE (D019) Bulk Modulus 55

50

Bulk modulus (GPa)

well known that REs go through the ”lanthanide contraction”, where the atomic volume gradually decreases with increasing number of f-electrons. The atomic volumes of the REs and the D019 compounds are shown in Fig. 2. Our DFT calculations correctly reproduce the lanthanide contraction for the pure REs, and we find a similar contraction across the series of Mg3 RE D019 compounds. Finally, we note that Eu- and Yb-containing systems do not follow the trends of the rest of the lanthanides. These deviations will be seen repeatedly throughout this work, and are well known and experimentally observed [81].

La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

(b)

D019 3/4 BMg+1/4 BRE

45

40

35

30

25

Sc Y

La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

Fig. 3. (a) Bulk modulus of pure REs using DFT calculations (red triangles) and experimental numbers (blue circles) [82]. The DFT values reproduce the experimental results well, including the very soft bulk moduli of Eu and Yb. (b) Bulk moduli of Mg3 RE compounds in the D019 structure. The compounds carry the same trend as the pure RE elements: increasing stiffness with increasing number of f-electrons, with the exception of Mg3 Eu and Mg3 Yb, which are significantly softer than the rest of the compounds. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Mg3 RE-D019 structures were calculated at two volumes: the ground-state DFT-relaxed volume (Veq ) and the DFT bulk Mg volume, in which the D019 structures are constrained to the DFT Mg lattice parameters (VMg ). We calculate the elastic constants of Mg3 RE-D019 at the lattice parameters of Mg to elucidate whether the elastic properties of the coherently constrained precipitate phase are significantly different from the relaxed geometry. The elastic constants for both volumes are shown in Table 2. In the constrained environment, the volume of the D019 precipitates is compressed to that of pure Mg, and the elastic constants are stiffer than in the relaxed state. We also use the Born stability criteria [87] to examine the mechanical stability of the D019 structure under constrained and relaxed conditions. For a hexagonal system, the stability criteria are: C 44 > 0; C 11  jC 12 j > 0; ðC 11 þ C 12 ÞC 33 

ð15Þ 2C 213

> 0:

A. Issa et al. / Acta Materialia 65 (2014) 240–250

245

Table 1 Elastic constant calculations (GPa) of the pure RE elements compared with previous DFT calculations and experiment, where available. Ce has the largest discrepancy between DFT and experiment, a known anomaly in DFT literature. Element

Structure

C 11

C 12

C 13

C 44

C 33

Ref.

Sc

hcp

106.2 104.8 99.3

38.3 37.7 39.7

30.3 29.3 29.4

33.9 31.6 27.7

106.5 105.2 107

Present DFT [84] Expt. [86]

Y

hcp

75.2 78.0 83.4

22.5 24.8 29.1

21.0 22.7 19.0

24.5 26.4 26.9

77.1 82.4 80.1

Present DFT [84] Expt. [86]

La

dhcp

47.9 51.4

17.2 17.27

10.1 10.4

15.3 13.92

56.7 54.6

Present DFT [80]

Ce

fcc

42.2 43.5 24.1

23.6 23.6 10.2

23.6 23.6 10.2

21.5 21.7 19.4

42.2 43.5 24.1

Present DFT [80] Expt. [85]

Pr

dhcp

58.3 60.8 50.8

22.9 25.4 24.6

15.5 17.9 14.8

17.7 17.4 14.9

64.2 67.3 55.7

Present DFT [80] Expt. [86]

Nd

dhcp

61.4 65.2 58.8

23.1 25.9 24.6

16.0 17.8 16.2

19.1 19.1 16.2

67.8 71.8 65.1

Present DFT [80] Expt. [86]

Pm

dhcp

66.0 70.3

29.6 24.6

17.6 18.6

20.5 21.0

81.4 77.2

Present DFT [80]

Sm

rhom

71.0 61.8

28.5 21.3

19.4 24.6

19.9 18.6

58.5 58.6

Present DFT [80]

Eu

bcc

17.5 16.5

10.9 10.6

10.9 10.6

17.1 16.3

17.5 16.5

Present DFT [80]

Gd

hcp

67.5 68.3 76.8

19.6 21.0 31.5

24.8 30.0 19.1

23.9 21.0 23.8

78.0 80.3 79.0

Present DFT [80] Expt. [86]

Tb

hcp

69.9 68.4 78.4

21.0 20.0 27.8

24.3 28.6 20.1

25.0 21.9 25.2

79.8 79.3 82.4

Present DFT [80] Expt. [86]

Dy

hcp

72.5 70.9 81.0

21.8 20.5 37.4

23.0 28.8 18.8

25.3 24.0 26.8

78.3 84.7 82.2

Present DFT [80] Expt. [86]

Ho

hcp

80.7 75.4 79.9

23.2 22.3 26.6

23.1 29.6 18.5

28.7 26.7 28.5

82.2 85.0 82.1

Present DFT [80] Expt. [86]

Er

hcp

82.5 81.5 86.8

26.0 24.3 28.5

22.4 28.3 24.2

28.2 28.9 26.5

85.5 88.0 82.1

Present DFT [80] Expt. [86]

Tm

hcp

88.1 88.4 92.5

27.4 25.6 33.5

23.4 28.0 25.0

31.0 30.2 28.2

90.5 94.2 81.5

Present DFT [80] Expt. [90]

Yb

fcc

21.6 23.2 18.6

14.0 12.9 10.4

14.0 12.9 10.4

15.3 17.4 17.7

21.6 23.2 18.6

Present DFT [80] Expt. [91]

Lu

hcp

91.6 95.8 91.0

26.6 31.4 31.9

21.0 31.4 28.8

29.7 32.7 29.1

90.5 98.1 84.0

Present DFT [80] Expt. [86]

Mg

hcp

62.6 67.5 63.5

26.0 24.8 25.9

20.9 24.1 21.7

13.3 24.0 18.4

64.9 72.4 66.4

Present DFT [80] Expt. [92]

We calculate and list the stability conditions in Table 3, where negative quantities denote a mechanical instability. We note that while all the D019 compounds are energeti-

cally metastable, the majority are mechanically stable. However, Mg3 Yb is mechanically unstable under relaxed conditions since it fails two stability criteria: its C 44 elastic

246

A. Issa et al. / Acta Materialia 65 (2014) 240–250

Table 2 DFT calculated elastic constants (GPa) of D019 structures constrained to bulk Mg lattice parameters (VMg ) and completely relaxed in DFT (Veq ). Element V Mg3 Sc Mg3 Y Mg3 La Mg3 Ce Mg3 Pr Mg3 Nd Mg3 Pm Mg3 Sm Mg3 Eu Mg3 Gd Mg3 Tb Mg3 Dy Mg3 Ho Mg3 Er Mg3 Tm Mg3 Yb Mg3 Lu

C 12

C 11 Mg

eq

78.2 94.1 94.8 104.0 96.7 96.0 96.0 95.7 71.5 94.5 93.3 104.3 90.9 103.7 90.2 71.4 89.1

Mg

V

V

77.8 69.1 57.0 60.6 62.3 63.9 66.0 66.8 44.0 68.5 69.1 68.1 68.6 64.8 63.1 30.2 65.9

36.7 44.6 51.3 43.7 47.8 46.9 45.8 45.1 50.0 44.5 45.0 35.7 46.3 33.8 44.5 49.2 43.2

C 13 V

eq

32.8 34.4 28.8 28.1 28.9 29.5 29.6 30.8 17.4 31.6 33.0 34.8 35.7 39.6 40.6 40.6 39.0

constant is negative and C 44 is smaller than the magnitude of C 12 . The mechanical instability of Mg3 Yb in the relaxed structure precludes its formation as a bulk phase, a prediction that is consistent with experiment [88]. Finally, we note that C 33 is greater than C 11 for all Mg3 RE compounds in both environments, indicating greater stiffness along the c or [0001] direction. 3.3. Coherency strain of Mg-D019 Fig. 4 shows the DFT-predicted coherency strain calculations for the Mg–Mg3 Ce system in both the basal and prismatic planes, an illustration of the calculations performed for all Mg–Mg3 RE systems. Fig. 5 summarizes

V

C 33

Mg

V

17.5 21.9 23.9 20.4 22.1 21.6 20.9 21.1 19.7 20.8 20.8 18.8 20.6 18.9 20.7 18.7 20.6

eq

18.8 18.8 18.8 17.1 17.9 17.9 17.6 18.1 14.7 17.9 18.3 18.3 18.4 18.7 18.8 14.2 19.5

C 44

Mg

V

V

eq

VMg

Veq

89.0 101.4 118.5 121.1 114.3 112.4 110.3 108.9 101.0 106.0 104.6 108.2 101.3 105.7 99.2 97.5 97.7

94.7 84.2 75.3 76.4 78.9 80.6 81.4 83.4 59.9 83.7 85.2 85.1 85.8 85.5 85.0 63.7 84.8

20.8 24.7 21.7 30.4 24.4 24.5 25.1 25.2 8.6 25.0 24.1 34.2 22.3 34.8 22.8 11.1 22.9

22.4 17.2 14.0 15.2 16.6 17.2 18.2 17.9 12.3 18.4 18.0 16.6 16.4 12.6 11.2 5.1 13.4

these calculations by showing the coherency strain energies at 10% D019 (2.5% RE) for all Mg–Mg3 RE systems. We selected this composition as it is representative of a Mg– RE alloy. The most significant observation is that the coherency strain energy in the basal plane is higher than that in the prismatic plane for all Mg–RE binary systems. This anisotropy in strain energy shows that prismatic habit planes are more energetically favorable for D019 phases in Mg–RE binary systems. Our prediction from DFT is fully consistent with experimental observation of Mg3 RE D019 precipitates adopting prismatic plane morphology [[32] and references therein]. We can now identify strain as the controlling factor in Mg3 RE D019 precipitate morphology, thus answering the long-standing question of why some

Table 3 Mechanical stability conditions for D019 structures when constrained to Mg lattice parameters (VMg ) and when relaxed in DFT (Veq ). For a compound to be stable, all three stability criteria must be greater than zero. Relaxed (equilibrium) Mg3 Yb is mechanically unstable, failing two of the three stability criteria.

Mg3 Sc Mg3 Y Mg3 La Mg3 Ce Mg3 Pr Mg3 Nd Mg3 Pm Mg3 Sm Mg3 Eu Mg3 Gd Mg3 Tb Mg3 Dy Mg3 Ho Mg3 Er Mg3 Tm Mg3 Yb Mg3 Lu

ðC 11 þ C 12 ÞC 33  2C 212 (GPa2 )

C 11  jC 12 j (GPa)

C 44 (GPa) VMg

Veq

VMg

Veq

VMg (103 )

Veq ð103 Þ

20.8 24.7 21.7 30.4 24.4 24.5 25.1 25.2 8.6 25.0 24.1 34.2 22.3 34.8 22.8 11.1 22.9

22.4 17.2 14.0 21.5 16.6 17.2 18.2 17.9 12.3 18.4 18.0 16.6 16.4 12.6 11.2 5.1 13.4

41.5 49.5 43.5 60.3 48.9 49.1 50.2 50.6 21.5 50.0 48.3 68.6 44.6 69.9 45.7 22.2 45.9

45.0 34.7 28.2 32.5 33.4 34.4 36.4 36.0 26.6 36.9 36.1 33.3 32.9 25.2 22.5 10.4 26.9

9.6 13.1 16.2 17.0 15.5 15.1 14.8 14.4 11.5 13.9 13.6 14.4 13.0 13.8 12.5 11.1 12.1

9.8 1.1 5.8 6.2 6.6 6.9 7.2 7.5 3.2 7.7 8.0 8.1 8.3 8.2 8.1 4.1 8.1

A. Issa et al. / Acta Materialia 65 (2014) 240–250

247

7.20

(a)

7.00 6.90 6.80 6.70 6.60

19

aD0 Lattice Parameter (Å)

7.10

6.50

aMg

6.40 6.30

Sc Y La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu 5.40

19

Fig. 4. Coherency strain energy as a function of composition for Mg– Mg3 Ce system. Basal strain energy is higher than prismatic for the entire composition range, a persistent fact for all the lanthanide systems. Lower prismatic strain energy indicates prismatic planes are energetically favored to be habit planes.

cD0 Lattice Parameter (Å)

5.35

Coherency Strain Energy (meV/atom)

Prismatic Versus Basal Coherency Strain: Mg/Mg3RE (D019) 10

Basal Prismatic

(b)

5.30 5.25 5.20 5.15

cMg

5.10 5.05 Sc Y La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

8

Fig. 6. (a) Lattice parameter a of the D019 structures, aD019 ¼ 2aMg . The corresponding lattice parameter of Mg is indicated as a horizontal line. (b) Lattice parameter c of the D019 structures, with the c lattice parameter of Mg indicated as a horizontal line.

6

4

2

0 Sc Y La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

Fig. 5. Coherency strain energy for basal and prismatic planes at 10% Mg3 RE (D019 ) (RE = name of the alloying element on the horizontal axis) using direct anharmonic DFT. Basal strain energy is higher than Mg for all RE systems, meaning prismatic planes will be energetically favored as habit planes.

precipitates in Mg–RE alloys adopt such morphologies. Identification of strain as the decisive factor in precipitate morphology is critical to both understanding the age hardening in Mg–RE alloys and to potential future alloy development using non-RE elements. We also note three results from Fig. 5. First, the coherency strain energies of Mg–Mg3 Eu and Mg–Mg3 Yb are higher than their neighboring lanthanide elements. Second, the basal coherency strain energy is monotonically decreasing for the lanthanide series, with the exception of Mg– Mg3 Eu and Mg–Mg3 Yb. Third, the prismatic coherency strain energy decreases in every lanthanide system up to Mg–Mg3 Tb, then stays fairly constant for the rest of the Mg–Mg3 RE systems, again with Mg–Mg3 Eu and Mg– Mg3 Yb deviating from the trend. All three observations

are explained by examining the lattice parameters of Mg– Mg3 RE systems, which are shown in Fig. 6. Since D019 is a superstructure of hcp, the lattice parameters of D019 can be described as multiples of the hcp structure lattice parameters: aD019 =bD019 =2ahcp and cD019 =chcp . From Fig. 6, the large coherency strain values of Mg–Mg3 Eu and Mg–Mg3 Yb can be accounted for by the large lattice mismatch between Mg and the D019 precipitates, over 10% mismatch in the a lattice parameter in the case of Mg–Mg3 Eu, for example. Similar reasoning qualitatively accounts for the monotonically decreasing basal coherency strain: aD019 gets smaller with increasing atomic number but is larger than 2aMg for all lanthanides, thus the gradual decrease in basal coherency strain energy with increasing atomic number. In Fig. 6b, we note that the mismatch between cD019 and cMg reaches a minimum of close to zero for the Mg–Mg3 Dy system, after which the mismatch starts to increase again as cD019 is now smaller than cMg . For Mg– Mg3 RE systems after Mg3 Dy, the lattice mismatch becomes smaller in a and larger in c, approximately compensating for each other, hence the approximately constant prismatic coherency strain energy for the Mg–Mg3 (Dy-Tm) systems. There remains the question of whether there is a habit plane with lower coherency strain energy than both prismatic and basal which was not explicitly calculated. Since,

248

A. Issa et al. / Acta Materialia 65 (2014) 240–250 (0001)

(a)

Energy (meV/atom)

8

6

4

2

_

(b)

2

4

6

8

(1010)

Energy (meV/atom)

Fig. 8. Coherency strain energy as a function of direction in Mg–Mg3 Gd at a composition of 50% Mg3 Gd (12.5% Gd) using harmonic DFT. Coherency strain is maximized for a basal habit plane and minimized for a prismatic one. Pyramidal planes are energetically between prismatic and basal. Harmonic DFT results in the basal and prismatic direction are in excellent agreement with the anharmonic DFT. The inset shows a threedimensional plot of the harmonic coherency strain energy as a function of direction.

(c)

Fig. 7. Basal and prismatic coherency strain curves using harmonic elasticity theory and DFT calculations. Solid lines indicate the coherency strain curves using the harmonic elasticity theory equations, and dashed lines indicate the curves obtained through DFT calculations. The close agreement between DFT and harmonic theory indicates that harmonic elasticity equations are sufficient to quantify the strain energy of a coherent interface in an arbitrary direction

technically, the number of possible habit planes between precipitate and matrix is infinite, it is not possible to calculate the coherency strain energy of every candidate habit plane. Laks et al. have shown [71], for cubic systems, how harmonic elasticity theory can be used to analytically express the coherency strain energy in terms of the elastic constants of the two phases. The extrema of this analytic expression then give the softest and hardest habit planes. However, it is not clear how far the harmonic theory can be applied, particularly in systems with very large size mismatch, such as some of the Mg–Mg3 RE systems desribed here. To test if harmonic elasticity theory can correctly describe Mg–Mg3 RE systems, we must compare the

coherency strain energy curves calculated using DFT, which accounts for harmonic and anharmonic strain energy, with analytically calculated harmonic elasticity theory curves (we derive these equations in the Appendix) which use DFT-calculated elastic constants. We perform such a comparison for the Mg–Mg3 La, Mg–Mg3 Eu and Mg–Mg3 Gd systems in Fig. 7. We chose these systems because of their large lattice mismatch with Mg, where we expect the largest anharmonic contributions to the coherency strain energy. For harmonic coherency strain calculations, we use the DFT-calculated elastic constants for Mg and Mg3 RE, each at their relaxed volumes. We show in Fig. 7 that the error in the prismatic and basal directions is 2 meV/atom or less. The favorable comparison gives us confidence that harmonic elasticity is the dominant contribution to the coherency strain energy. We can thus conclude that harmonic elasticity theory is sufficient to calculate the coherency strain energy and that the anharmonic contributions to the strain energy are small. The conclusion that harmonic elasticity is accurate means we can now use the elasticity equations to calculate the coherency strain energy for an arbitrary direction, which we do for basal, prismatic and pyramidal planes in the Mg–Mg3 Gd system at a concentration of 50% Mg3 Gd (12.5% Gd). Fig. 8 shows a polar plot of the coherency strain energy as a function of direction (^k) between the [1 0 1 0] prismatic plane and the [0 0 0 1] basal plane. Note that because elastic constants are isotropic with respect to rotations around the z direction [89], the quadrant shown suffices to investigate all possible habit planes. The transverse isotropy in hexagonal systems is seen in the inset in Fig. 8, where the prismatic plane coherency strain energy is isotropic with respect to rotations around the z axis, a direct result of

A. Issa et al. / Acta Materialia 65 (2014) 240–250

the harmonic elasticity equations. From Fig. 8 we see that the coherency strain energy is lowest along the horizontal axis, indicating that a prismatic habit plane would minimize the coherency strain energy of this system. 4. Conclusions In this paper, we have answered the long-standing questions about metastable precipitate morphology in high-strength Mg–RE alloys. We used first-principles calculations to assess the energetic stability, elastic constants and coherency strain of Mg3 RE precipitates in the D019 structure. We find: 1. The Mg3 RE-D019 structures all have negative formation energies that are competitive with the stable phase. The formation energies of the D019 structures become less negative with increasing RE atomic number, a possible reason why these b00 -D019 phases are not observed in Mg–heavy RE systems. 2. The lattice mismatch between Mg3 RE D019 precipitates and Mg explains the coherency strain behavior across the lanthanide series. The monotonic decrease in lattice mismatch in the basal plane leads to a consistent decrease in the basal coherency strain energy. The prismatic plane energy, on the other hand, is fairly constant for the last third of the lanthanides as a result of the competing effects of smaller mismatch in lattice parameter a and a bigger mismatch in lattice parameter c. 3. All the lanthanide D019 structures, and Mg3 Y, energetically prefer prismatic habit planes with hcp Mg over basal planes. This general tendency implies prismatic planes are indeed the precipitate habit planes, in agreement with experiment. The strain-dictated precipitate morphology explains the exceptional strength of Mg–RE alloys. Our core finding of prismatic habit plane morphology in the Mg–RE systems elucidates the experimentally observed morphology and points the way towards possible alloy design schemes. Attempts to replace RE alloying elements in Mg, while keeping the favorable prismatic plate precipitates, must take into account strain as a controlling factor in precipitate morphology. Acknowledgment The authors gratefully acknowledge the support of the Ford–Boeing–Northwestern (FBN) alliance, Award No. 81132882. Appendix A. Appendix Using continuum elasticity theory, we can derive equations that describe the coherency strain energy in an arbitrary plane. The general equation for crystal strain energy per unit volume, DE, is:

1 DE ¼ C ij i j : 2

249

ð16Þ

Following the standard Einstein notation of summing over repeated indices, Eq. (16) has nine terms on the right-hand side. To find the biaxial strain energy for a given crystal, DEepi , Eq. (16) is minimized with respect to strain along the epitaxial plane normal. Doing so for cubic systems produces [71]: ! 9 B epi ^ 2 DE ðkÞ ¼ V 0 B 1  0 ; ð17Þ 2 C 11 ð^kÞ where C 011 ð^kÞ is the out-of-plane elastic constant which is a function of direction ^k, (e.g. stiffness along the h1 1 1i direction for an epitaxial strain in the [1 1 1] plane). Because of reduced symmetry, the analogous epitaxial strain energy equation for a hexagonal crystal contains more terms:   0 0 0 0 V 0 0 02 0 0 C 12 2 þ C 13 3 epi DE ¼ C 22 2 þ C 12 2 þ C 032 02 03 2 C 011  0 0  0 0 0 0 C 12 2 þ C 13 3 0 0 0 0 02 þ C 13 3 þ C 32 2 3 þ C 33 3 : ð18Þ C 011 Similarly to Eq. (17), all the primed variables are a function of direction ^k. Again, C 011 is taken to be the out-of-plane elastic constant. The primed variables in Eqs. (17) and (18) are calculated using the tensor transformation laws [87]. Finally, to produce the coherency strain curves as a function of composition for a particular habit plane, as in Fig. 7, we use the equation: epi epi ^ DEhex CS ðk; xÞ ¼ ð1  xÞDE Mg þ xDE Mg3 RE ;

ð19Þ

Epi where DEEpi Mg and DE Mg3 RE are calculated using Eq. (18).

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

Partridge P. Metall Rev 1967;12:169–94. Zhu S, Gibson M, Easton M, Nie JF. Scripta Mater 2010;63:698–703. Li R, Nie JF, Huang G, Xin Y, Liu Q. Scripta Mater 2011;64:950–3. Gao X, He S, Zeng X, Peng L, Ding W, Nie JF. Mater Sci Eng A– Struct 2006;431:322–7. Honma T, Ohkubo T, Hono K, Kamado S. Mater Sci Eng A–Struct 2005;395:301–6. Oh-ishi K, Watanabe R, Mendis C, Hono K. Mater Sci Eng A–Struct 2009;526:177–84. Nie JF, Gao X, Zhu S. Scripta Mater 2005;53:1049–53. Honma T, Ohkubo T, Kamado S, Hono K. Acta Mater 2007;55:4137–50. Suzuki M, Kimura T, Koike J, Maruyama K. Scripta Mater 2003;48:997–1002. Zheng K, Dong J, Zeng X, Ding W. Mater Sci Eng A–Struct 2008;489:44–54. Shao X, Yang Z, Ma X. Acta Mater 2010;58:4760–71. Nie JF, Oh-ishi K, Gao X, Hono K. Acta Mater 2008;56:6061–76. Zhu YM, Morton AJ, Nie JF. Acta Mater 2010;58:2936–47. Yamada K, Hoshikawa H, Maki S, Ozaki T, Kuroki Y, Kamado S, et al. Scripta Mater 2009;61:636–9. Hisa M, Barry JC, Dunlop GL. Phil Mag A 2002;82:37–41. Nie JF. Acta Mater 2000;48:1691–703. Ping D, Hono K, Nie JF. Scripta Mater 2003;48:1017–22. Rokhlin L, Dobatkina T, Tarytina I, Timofeev V, Balakhchi E. J Alloys Compd 2004;367:17–9.

250

A. Issa et al. / Acta Materialia 65 (2014) 240–250

[19] Nie JF. Scripta Mater 2003;48:1009–15. [20] He S, Zeng X, Peng L, Gao X, Nie JF, Ding W. J Alloys Compd 2007;427:316–23. [21] Apps P, Karimzadeh H, King J, Lorimer G. Scripta Mater 2003;48:1023–8. [22] Antion C. Acta Mater 2003;51:5335–48. [23] Zhang M-X, Kelly P. Scripta Mater 2003;48:379–84. [24] Smola B, Stulikova I, Vonbuch F, Mordike B. Mater Sci Eng A– Struct 2002;324:113–7. [25] Vostry P, Smola B, Stuliokova I, von Buch F, Mordike BL. Phys Stat Sol A 1999;42:491–500. [26] Saito K, Hiraga K. Mater Trans 2011;52:1860–7. [27] Nishijima M, Hiraga K, Yamasaki M, Kawamura Y. Mater Trans 2006;47:2109–12. [28] Rokhlin L. Fiz Metal Metalloved 1985;59:1188–93. [29] Rokhlin L. Fiz Metal Metalloved 1983;55:733–8. [30] Rokhlin L. Fiz Metal Metalloved 1987;63:146–50. [31] Rokhlin L, Nikitina NI. Fiz Metal Metalloved 1986;62:781–6. [32] Nie J-f. Metall Mater Trans A 2012:1–96. [33] Smola B, Stulikova´ I. J Alloys Compd 2004;381:L1–2. [34] Shin D, Wolverton C. Acta Mater 2010;58:531–40. [35] Yasi JA, Nogaret T, Trinkle DR, Qi Y, Hector LG, Curtin WA. Model Simul Mater Sci Eng 2009;17:055012. [36] Yasi JA, Hector Jr LG, Trinkle DR. Acta Mater 2010;58:5704–13. [37] Zhang H, Shang S, Saal JE, Saengdeejing A, Wang Y, Chen L-Q, et al. Intermetallics 2009;17:878–85. [38] Wang Y, Chen L-Q, Liu Z-K, Mathaudhu S. Scripta Mater 2010;62:646–9. [39] Shin D, Wolverton C. Acta Mater 2012;60:5135–42. [40] Shin D, Wolverton C. Scripta Mater 2010;63:680–5. [41] Saal JE, Wolverton C. Scripta Mater 2012. [42] Saal JE, Wolverton C. Acta Mater 2012;60:5151–9. [43] Kresse G, Hafner J. Phys Rev B: Condens Matter 1993;47:558–61. [44] Kresse G, Furthmu¨ller J. Phys Rev B: Condens Matter 1996;54:11169–86. [45] Kresse G, Joubert D. Phys Rev B: Condens Matter 1999;59:1758–75. [46] Kresse G, Hafner J. J Phys: Condens Matter 1994;6:8245–57. [47] Vanderbilt D. Phys Rev B: Condens Matter 1990;41:7892–5. [48] Perdew JP, Burke K, Ernzerhof M. Phys Rev Lett 1996;77:3865–8. [49] Georges A. In: AIP conf. proc., Aip; 2004. p. 3–74. [50] Temmerman W, Petit L, Svane A, Szotek Z, Luders M, Strange P, et al. Handbook on the physics and chemistry of rare earths 2009;39:1–112. [51] Zhou F, Ozolins V. Phys Rev B: Condens Matter 2009;80:125127. [52] Duthie J, Pettifor D. Phys Rev Lett 1977;38:564–7. [53] Gao MC, Rollett AD, Widom M. Phys Rev B: Condens Matter 2007;75:174120:1–174120:16. ¨ nlu¨ N, Shiflet GJ, Mihalkovic M, Widom M. Metall [54] Gao MC, U Mater Trans A 2005;36:3269–79. [55] Mao Z, Seidman DN, Wolverton C. Acta Mater 2011;59:3659–66. [56] Gao M, Rollett a, Widom M. Calphad 2006;30:341–8. [57] Lu S, Hu Q-M, Yang R, Johansson B, Vitos L. Comp Mater Sci 2009;46:1187–91. [58] Tao X, Ouyang Y, Liu H, Feng Y, Du Y, He Y, et al. J Alloys Compd 2011;509:6899–907.

[59] Wu B, Zinkevich M, Aldinger F, Wen D, Chen L. J Solid State Chem 2007;180:3280–7. [60] Eriksson O, Brooks M, Johansson B. J Less-Common Met 1990;158:207–20. [61] Okamoto H. J Phase Equilib 1996;17:553. [62] Cacciamani G, Saccone A, De Negri S, Ferro R. J Phase Equilib 2002;23. [63] Saccone A, Delfino S, Maccio D, Ferro R. J Phase Equilib 1994;15:128. [64] Du Z, Yang X, Ling G, Liu H. Z Metallkd 2004;95:1115–9. [65] Guo C, Du Z. J Alloys Compd 2006;422:102–8. [66] Okamoto H. J Phase Equilib 1993;14:534–5. [67] Okamoto H. J Phase Equilib 1992;13:103–4. [68] Okamoto H. In: Massalski TB, editor. Binary alloy phase diagrams, 2nd ed.; 1990. p. 2549–50. [69] Nayeb-Hashemi A, Clark J. J Phase Equilib 1989;10:23–7. [70] Nayeb-Hashemi A, Clark J. In: Massalski T, editor. Binary alloy phase diagrams, 2 ed.; 1990. [71] Laks DB, Ferreira LG, Froyen S, Zunger A. Phys Rev B: Condens Matter 1992;46:12587–605. [72] van de Walle A, Ceder G. J Phase Equilib 2002;23:348–59. [73] van de Walle A, Asta M. Model Simul Mater Sci Eng 2002;10:521–38. [74] Shang S-L, Wang Y, Kim D, Liu Z-K. Comp Mater Sci 2010;47:1040–8. [75] Birch F. Phys Rev B: Condens Matter 1947;71:809–24. [76] Sholl DS, Steckel JA. Density functional theory: a practical introduction. Hoboken, NJ: John Wiley & Sons; 2009. [77] Guo X-Q, Podloucky R, Freeman A. J Mater Res 2011;6:324–9. [78] Mao Z, Chen W, Seidman D, Wolverton C. Acta Mater 2011;59:3012–23. [79] Wang A, Shang S, Du Y, Kong Y, Zhang L, Chen L-Q, et al. Comp Mater Sci 2010;48:705–9. [80] Ouyang Y, Tao X, Zeng F, Chen H, Du Y, Feng Y, et al. Physica B 2009;404:2299–304. [81] Rokhlin L. Magnesium alloys containing rare earth metals. New York; 2003. [82] Gschneidner KA. J Phase Equilib 1990;11:216–24. [83] Singh N, Singh SP. Phys Rev B: Condens Matter 1990;42:1652–8. [84] Shang S, Saengdeejing a, Mei Z, Kim D, Zhang H, Ganeshan S, et al. Comp Mater Sci 2010;48:813–26. [85] Stassis C, Gould T, McMasters OD, Gschneidner KA, Nicklow RM. Phys Rev B: Condens Matter 1979;19:5746–53. [86] Singh D, Varshni YP. Phys Rev B: Condens Matter 1981;24:4340–7. [87] Nye J. Physical properties of crystals. New York, NY: Oxford University Press; 1957. [88] Rokhlin L, Nikitina N, Tarytina I. Fiz Met Metalloved 1993;76:82–7. [89] Urusovskaya AA. In: Shuvalov LA, editor. Modern crystallography IV. Berlin: Springer-Verlag; 1988. p. 63. [90] Lim CM, Edwards C, Dixon S, Palmer SB. J Magn Magn Mater 2001;234:387–94. [91] Stassis C, Loong C-K, Theisen C, Nicklow RM. Phys Rev B: Condens Matter 1982;26:4106–10. [92] Slutsky LJ, Garland C. Phys Rev B: Condens Matter 1957;107:972–6.