Acta Materialia 53 (2005) 3729–3736 www.actamat-journals.com
Prediction of structural stabilities of transition-metal disilicide alloys by the density functional theory G. Shao
*
BCAST - Brunel Centre for Advanced Solidification Technology, Brunel University, Uxbridge, Middlesex UB8 3PH, UK Received 14 March 2005; received in revised form 13 April 2005; accepted 14 April 2005 Available online 3 June 2005
Abstract Transition metal silicide based refractory composite alloys have great potential for applications as the next generation of turbine airfoil materials, due to their significantly higher operating temperatures than advanced Ni-based superalloys. In order to achieve the most of refractory silicide based materials, multi-component alloying has been considered to be the most useful route for developing alloys of a good balance of creep resistance, fracture toughness, oxidation resistance, and room-temperature ductility. In this work, structural stabilities of disilicides of the group IVB to VIB transition metal elements are studied using the ab initio density functional theory, shedding light on alloying/bonding behaviour at the electronic level. The results are also useful for developing full thermodynamic databases using the CALPHAD approach, where lattice stabilities of counter-phases are critical for accurate thermodynamic descriptions of multi-component systems. 2005 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Ab initio electronic theory; CALPHAD; Density functional; Thermodynamics; Phase transformation
1. Introduction Intermetallic alloys that have high modulus and good creep resistance at high temperatures have recently attracted a great deal of attention in the international community of aerospace materials researchers. Of particular interest are composite materials based on niobium disilicide (C40 structure) or molybdenum disilicide (C11b structure), owing to their significantly higher operating temperatures than the current advanced Ni-based superalloys. These composite silicide materials show great promise as candidates for the next generation of turbine airfoil materials [1–4]. The potential applications of these alloys at high temperature require a balance of high creep strength, good oxidation
*
Tel.: +44 01895 266 407; fax: +44 01895 269 758. E-mail address:
[email protected].
resistance and acceptable low-temperature fracture toughness. Multi-alloying of NbSi2 and MoSi2 with both Al and group IVB to VIB transition metal elements such as Cr, V, Ta, Hf, Ti has been employed to improve the properties of the composites [4–8]. For the purpose of advanced design of disilicide based composites, one needs a good understanding of phase equilibria and thermodynamic properties of competing structures in the corresponding multi-component systems. It has been demonstrated that a most effective tool for designing multi-component alloys is the CALPHAD approach, which relies on assessed thermodynamic databases [9]. A major difficulty for assessing a multi-component system lies in the uncertainty in determining lattice stabilities (enthalpy and entropy) of counter phases, which cannot be measured directly due to their metastable or even unstable nature. For example, to assess the thermodynamic properties of the C40 structure of NbSi2 in the
1359-6454/$30.00 2005 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2005.04.025
3730
G. Shao / Acta Materialia 53 (2005) 3729–3736
Nb–Mo–Si system, one needs to have the lattice stability data for the same structures in both the equilibrium stoichiometry NbSi2 and the metastable stoichiometry MoSi2, with the latter being denoted as the counterphase (or hypothetical structure) of the former. A thermodynamic database in the framework of the CALPHAD methodology would be reliable only if the lattice stabilities of the counter-phase structures are accurately defined, this being particularly true when the database is to be extrapolated into metastable regimes that are important for materials processing. There are four crystal structures for refractory disilicides of interest as aerospace materials, the C40 (hP9), C11b (tI6), C54 (oF24) and C49 (oC12) [10]. Table 1 shows the experimentally observed crystal structures for the group IV to VI transition metal disilicides, revealing a strong correlation between the periodic table and the experimentally observed equilibrium structures of the phases. This suggests a dominant electronic origin for the evolution of stable structures in these disilicides. It is therefore expected that the fundamental reasons for structural transitions in refractory silicides will be understood at the electronic level, this being the task of the current work. The author aims to provide some theoretical insight into alloying refractory disilicides with transition metals, with the objectives being to: (a) calculate heats of formation of different MSi2 structures that would be difficult to determine otherwise, (b) probe the physical ground for structural competitions and bonding nature in MSi2, and (c) compare the results from different DFT methods.
2. Methodology Calculations in this work have been performed within the framework of the density functional theory (DFT) [11], which states that all ground state properties are functionals of the charge density q. Specifically, the total energy Et is Et ¼ T ½q þ U ½q þ Exc ½q;
ð1Þ
where T [q] is kinetic energy, U[q] the potential field, and Exc[q] the exchange-correlation energy. Self-consistent solution of Eq. (1) is equivalent to treating the manybody problem of interacting electrons and nuclei by
Table 1 Experimentally observed equilibrium MSi2 structures TiSi2 C54 (oF24) ZrSi2 C49 (oC12) HfSi2 C49 (oC12) a
VSi2 C40 (hP9) NbSi2 C40 (hP9) TaSi2 C40 (hP9)
High temperature equilibrium structure is C40.
CrSi2 C40 (hP9) MoSi2 C11b(tI6)a WSi2 C11b(tI6)a
solving a series of one-electron Kohn–Sham (KS) equations [12]. The final term in Eq. (1), the exchangecorrelation energy, requires some approximation for this method to be computationally tractable. Unless defined otherwise, the generalised gradient approximation (GGA) functionals have been used throughout this work, as they provide a better overall description of the electronic subsystem than the local density approximation (LDA) functionals [13,14]. Two well-tested computing codes have been used for the calculations in this work. The CASTEP code employs the plane-wave basis set to treat valence electrons and pseudopotentials to approximate the potential field of ion cores (including nuclei and tightly bond core electrons) [15], which is an efficient DFT tool especially suitable for calculations involving necessary geometric optimisation especially for metastable structures. The other code, DMol3 [16], which uses the numeric basis set to treat all electrons (with each basis function that corresponds to an atomic orbital being given numerically on an atomic-centred spherical-polar mesh), has been used to compare with enthalpies of formation from CASTEP.
3. Results and discussion 3.1. Heat of formation The heat of formation of a compound is calculated by subtracting the total energies of pure elements in their stable crystal structures from the total energy of the compound. For an ApBq compound, for example, the heat of formation H is defined as A Bq
DH ¼ Etotp
B pEA tot qE tot ;
ð2Þ
where the superscripts of the total energies Etot define the compound and constituent elements. Geometric optimisation was carried out for minimising the total energies of compounds and reference elemental structures against unit cell sizes and delocalised atomic coordinates, using the CASTEP code. Ultrasoft pseudo-potentials and PBE GGA functionals [14] were used for all CASTEP calculations in this work. All data were obtained for the relaxed structures via geometrically optimised unit cell sizes and delocalised atomic coordinates. The energy cutoff was set to be above 330 eV for both elements and silicides. Parameters for convergence control of geometry optimisation were 5 · 106 ˚ for maximum force, eV/atom for total energy, 0.01 eV/A ˚ for displacement. 0.02 GPa for stress, and 5 · 104 A The Monkhorst–Pack scheme [25] was used to produce a uniform grid of k-points along the three axes in reci˚ 1 procal space, and an actual separation of 0.04 A was chosen to guarantee the accuracy of the Brillouin
G. Shao / Acta Materialia 53 (2005) 3729–3736 -10 -20
C49
3d_series
-30
C54
-40
C40 C11b
-50 -60 -10
Hf (kJ / mole_atom)
zone sampling [15]. Table 2 shows the calculated unit cell sizes together with experimental data, showing overall errors smaller than 1%. The CASTEP predicted heats of formation of competing structures C11b, C40, C49 and C54 of MSi2 (with M as group IVB to VIB transition metals) are shown in Fig. 1. This figure predicts that the stable ground state structures for the group IVB to VIB disilicides are C49, C40 and C11b, being consistent with experimental observations (Table 1), except for TiSi2 (C49 instead of C54) and CrSi2 (C11b instead of C40). The predicted ground-state stable C11b structure for CrSi2 is consistent with previous workersÕ findings [17,18]. The calculated heats of formation for all four crystal structures are summarised in Table 3, together with available experimental measurements [19–22], and DFT calculated data from other workers [17]. Pank-
3731
-20
2 Ti
3
V
4 Cr
C49
4d_series
-30
C54 C40 C11b
-40 -50 -60 -70 2
Zr
3
Nb
4 Mo
-10
Table 2 CASTEP predicted MSi2 lattice parameters compared with available experimental data from [10] in parentheses M
Structure
Lattice parameter (0.1 nm) a
Ti
V
Cr
Zr
Nb
Mo
Hf
Ta
W
a
C40 C11b C54 C49 C40 C11b C54 C49 C40 C11b C54 C49 C40 C11b C54 C49 C40 C11b C54 C49 C40a C11b C54 C49 C40 C11b C54 C49 C40 C11b C54 C49 C40a C11b C54 C49
4.6764 3.1578 8.1790 3.5069 4.5239 3.1093 7.9434 3.3649 4.3544 3.0437 7.5926 3.0698 4.9878 3.2497 8.6384 3.6493 4.7793 3.2457 8.4301 3.5059 4.5676 3.1863 7.9950 3.3935 4.9230 3.2483 8.5665 3.6404 4.6867 3.2061 8.2815 3.4527 4.5582 3.2013 7.8389 3.1803
b
(8.252)
4.6050 (4.783) 13.4397
(4.571) 4.4981 13.0153 (4.431) 4.3346 11.7473
(3.721) (4.819)
5.1099 14.9186 (14.68)
4.7806 14.2725 (4.642) (3.202) 4.5656 11.9023
(3.672) (4.7821)
5.0212 14.5612 (14.570)
4.6687 13.8824 (4.614) (3.211)
High temperature structure.
4.5247 12.2309
5d_series
C49
-20 -30
C54 C40
-40
c 6.5317 8.3500 8.4696 3.5311 6.2900 7.7275 8.3163 3.4995 6.2762 7.4348 8.3504 3.8760 6.8668 9.1798 8.9058 3.6494 6.5320 8.2159 8.5901 2.5390 6.5388 7.7885 8.6566 3.3847 6.8146 8.9837 8.8482 3.6257 6.4336 8.0108 8.5060 3.5129 6.6037 7.7695 8.9595 4.1512
C11b -50
(8.540) (6.372)
-60
2 Hf
3 Ta
4W
Fig. 1. Heats of formation for MSi2 of different structures (M = Ti, V, Cr, Zr, Nb, Mo, Hf, Ta, W) at 0 K, predicted by geometric optimisation using the CASTEP code.
(4.364)
(3.683) (6.592)
(6.529) (7.851)
(3.641) (6.5695)
(6.414) (7.868)
hurst et al. [17] calculated the heats of formation of competing C11b and C40 structures of the CrSi2 and MoSi2, using the full-potential linear augmented plane waves (FP-LAPW) method, which are in excellent agreement with the CASTEP prediction of the present work. The energetic difference between the theoretical stable structures and the experimentally observed ones for TiSi2 and CrSi2 are very small (2.3 and 1.2 kJ per mole of atoms, respectively). To test if such small disparity were due to the pseudopotential approximation of core electrons, the all-electron code DMol3 was also used to calculate heats of formation in this work, on the basis of the CASTEP predicted unit cell sizes and subsequent structural relaxation for the delocalised atoms. The allelectron calculation predicts the same ground state stable structures as the CASTEP for both TiSi2 and CrSi2 (Table 3). As was pointed out by Pankhurst et al. [17], such small energy differences between the stable and metastable structures at ground-state could be outperformed by entropy differences at temperatures above zero Kelvin, leading to the observation of the metastable structure in metallurgical experimental temperature ranges.
3732
G. Shao / Acta Materialia 53 (2005) 3729–3736
Table 3 Calculated heats of formation of MSi2 compared to the literature data, with ground state stable structures emboldened M
Ti
Structure
C40 C11b C54
C49
V
C40 C11b
C54 C49 Cr
C40
C11b
C54 C49 Zr
C40 C11b C54
C49
Nb
C40
C11b
C54 C49 Mo
C40
C11b
Table 3 (continued) M
Structure
Hf (kJ per mole of atoms) Calc.
Method
Experimental
53.42 45.18 34.56 28.53 54.68 46.47 62.12 57.0 48.43 63.47
CASTEP DMol GGA CASTEP DMol GGA CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA
48.69 36.83 41.18 30.85 47.49 46.26 35.72 31.56 22.58
CASTEP DMol GGA CASTEP DMol GGA DMol LDA CASTEP DMol GGA CASTEP DMol GGA
37; 42; 50 [19]
38.19 26.26 44.67 34.65 [17] 39.37 29.17 47.16 36.03 [17] 34.17 22.70 15.83 6.30
CASTEP DMol GGA DMol LDA FP-LAPW CASTEP DMol GGA DMol LDA FP-LAPW CASTEP DMol GGA CASTEP DMol GGA
27 [19]
49.61 44.44 35.83 19.63 53.46 48.42 61.22 64.87 58.25 70.25
CASTEP DMol GGA CASTEP DMol GGA CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA
53.29 45.06 58.39 42.67 36.03 49.73 51.16 43.46 40.91 34.20
CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA CASTEP DMol GGA CASTEP DMol GGA
46.03 41.59 57.52 46.20 [17] 48.50 45.21 60.67
CASTEP DMol GGA DMol LDA FP-LAPW CASTEP DMol GGA DMol LDA
C54 C49
60; 59; 57 [19]
Hf
C40 C11b C54
C49
Ta
C40
C11b
C54 C49 W
C40
C11b
C54 C49
60.3 ± 1.7 [21]
53.7 ± 1.6 [21]
47.9 ± 2.1 [21]
Hf (kJ per mole of atoms) Calc.
Method
48.90 [17] 41.06 37.52 15.43 0.93
FP-LAPW CASTEP DMol GGA CASTEP DMol GGA
45.82 47.15 28.94 34.10 47.66 50.49 59.39 53.34 60.33 73.05
CASTEP DMol GGA CASTEP DMol GGA CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA
57.74 43.01 55.54 49.04 32.49 46.25 54.66 41.97 37.60 23.03
CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA CASTEP DMol GGA CASTEP DMol GGA
42.29 38.12 54.41 48.88 41.78 60.01 38.0 33.16 23.84 16.99
CASTEP DMol GGA DMol LDA CASTEP DMol GGA DMol LDA CASTEP DMol GGA CASTEP DMol GGA
Experimental
37; 39; 40 [19] 46 ± 2.5 [22]
27.3 ± 2.0 [22]
The predicted heats of formation using CASTEP and DMol3 are compared with experimental data of corresponding structures in Fig. 2. Both LDA and GGA (PW91 potentials [13]) were used for the DMol3 calculation. The trends of predicted heats of formation are largely in agreement. Except for the 5d disilicide HfSi2 and TaSi2, the CASTEP values are between the lower and upper limits defined by the DMol3 data. This is very interesting as it is well recognised that LDA tends to overbind atoms, so that the bond lengths and the cell volume are usually underestimated by a few per cent and the bulk modulus is correspondingly overestimated. GGA corrects this error but may underbind instead, leading to slightly long bond lengths and smaller bulk modulus. It should be pointed out that large discrepancies exist in experimental data of different workers (see Fig. 2 for VSi2 and TaSi2 data). This is consistent with the very slow transformation kinetics of
G. Shao / Acta Materialia 53 (2005) 3729–3736 0
60 ZR
NB
MO
HF
TA
CASTEP DMol GGA DMol LDA
H F ( kJ mole_at
-1
)
-20
TI
W
-30 -40 -50 -60 -70 -80
Fig. 2. Comparison of calculated heats of formation of MSi2 with corresponding experimental data. GGA was used for the CASTEP calculation and both GGA and LDA were used for the DMol3 calculation. The experimental data are from the review of de Boar et al. [19] (triangles) and the work of KleppaÕs group [20–22] (squares with error bars).
these refractory materials (due to low diffusion coefficient) [4,8], so that it was difficult for one to achieve the equilibrium condition during measurements. Therefore, it is reasonable to observe that more recent measurements (data with error bars) tended to lead to more negative heats of formation, largely due to improved techniques and better efforts on equilibration (cf. experimental data for TaSi2). Overall, there is good agreement between the theoretical and experimental data of heats of formation, except for the case of WSi2, for which the measured heat of formation for WSi2 (C11b) is outside the band of predicted data. Further efforts in experimental measurement could be helpful, considering that the melting point of WSi2 is higher than that of MoSi2 and the CASTEP predicted heat of formation similar to that of MoSi2 seems to be reasonable, as there is positive correlation between the melting points and heats of formation of silicides [23]. From the theoretical point of view, it will be useful to study if the 5d-series disilicides such as the C11b WSi2 are of abnormal lattice dynamics with respect to the 3d- and 4d-series. In Fig. 3, heats of formation with respect to those of the C40 structure are compared. The agreement between the CASPEP (GGA) and DMol3 (GGA) predictions is extremely good. Such good agreements also hold using other structures as references. This suggests that higher accuracy is achievable for energy differences between different structures of the same stoichiometry than for total energy or heat of formation, possibly due to the cancellation of errors when structures of similar states are used as references. This is useful for efficient use of ab initio input for CALPHAD modelling, as energy dif-
50
)
CR
-1
V
H F vs. C40 ( kJ mole_at
TI
-10
3733
V
CR
ZR
NB
MO
HF
TA
W
C11b C49
40
C54 C11b DMol
30
C49 DMol C54 DMol
20 10 0 -10 -20
Fig. 3. Heats of formation of MSi2 with respect to that of the C40 structure, calculated using the CASTEP and DMol3. For the DMol3 calculations, structural relaxation was carried out by geometric optimisation of the delocalised atoms using the unit cell sizes predicted by CASTEP. Notice both methods predict almost the same trends for energetic differences between different structures.
ferences between different structures are very important data for the latter. 3.2. Density of states Electron density of states (DOS) have been calculated in this work to elucidate the nature of bonding in transition metal disilicides of the C11b and C40 structures, using the CASTEP code with k-point separation of ˚ 1. Fig. 4 shows the DOS for the C40 and C11b 0.04 A CrSi2 structures. Both structures exhibit a major peak due to p–d hybridisation below zero energy. The semiconducting C40 structure is characterised by a stronger d–d hybridisation shoulder near the zero energy than the C11b structure, being in agreement with Pankhurst et al. [24]. In addition, a pseudo-gap exists above the d–d peak for both structures. DOS of various C40 and C11b MSi2 phases are shown in Figs. 5 and 6, respectively, being arranged in the same way as the positions of metals in the periodic table. It is seen firstly that the lower the number of delectrons, the higher the energy of the major p–d peak and the larger the value of DOS at zero energy (cf. the 3d-series). Secondly, the major p–d peak shifts to lower energy from the 3d- to 5d-series, this being attributed to corresponding enhanced p–d hybridisation. While such effects exist for DOS of both structures, the effects for the latter structure (C11b) are less pronounced. For example, the p–d peak energies of both structures are summarised in Fig. 7. This is attributable to a rigid band nature of these transition metal silicides, so that reduced overall electron concentration shifts the Fermi energy below the
3734
G. Shao / Acta Materialia 53 (2005) 3729–3736
15
p-d
5
(a)
p-d
(b)
C40_CrSi2
C11b_CrSi2
4 10 3
d-d
d-d 2 5 1
d p
0 -15
-10
-5
0
5
10
0 -20
p -15
-10
-5
d 0
5
10
15
20
Fig. 4. DOS for the C40 (a) and C11b (b) structures of CrSi2. Notice the major peak below the zero energy being attributed to the p–d hybridisation. The shoulders above this major peak are due to the d–d hybridisation. The semiconducting C40 structure is characterised with the d–d shoulder being more pronounced than that of the C11b structure.
15
15
15
p-d
p-d
C40 Ti
V
10
Cr 10
10
d-d d-d 5
5
5
p
d
0
-15
-10
-5
0
5
0 -15
Zr
-10
-5
0
5
0 -15
10
10
5
5
5
-15
-10
-5
0
5
0 -15
Hf
-10
-5
0
5
0 -15
Ta 10
10
5
5
5
-15
0
-10
-5
0
5
-15
0
5
10
-10
-5
0
5
10
-10
-5
0
5
10
W
10
0
-5
Mo
Nb
10
0
-10
0
-10
-5
0
5
-15
E (eV) Fig. 5. DOS of various MSi2 with the C40 structure. Notice how the p–d peak shifts for disilicides of the Cr-group. The emboldened metals correspond to experimentally observed stable C40 MSi2.
pseudo-gap [24]. These characteristics in DOS are consistent with the trends of heats of formation, with enhanced binding corresponding to lower energy of the
p–d and d–d peaks. Also, it makes these silicides more metallic to shift the Fermi energy below the pseudogap via lowering the overall electron concentration.
G. Shao / Acta Materialia 53 (2005) 3729–3736
5
5
3735
p-d
5
C11b
4
4
Ti
4
V
3
3
3
2
2
2
1
1
1
0
-20
4
-15
-10
-5
0
5
10
1
0
-20
4
Zr
3
-15
-10
-5
0
5
10
1
p-d
Nb
0
-20
4
3
3
2
2
2
1
1
1
p-d
0
-20
4
-15
-10
-5
0
5
10
1
0
-20
4
Hf
-15
-10
-5
0
5
10
1
0
-20
4
Ta
3
3
3
2
2
2
1
1
1
0 -20 -15 -10 -5
0
5
0
10 15 -20 -15 -10 -5
0
5
d-d
Cr
p -15
d
-10
-5
0
5
10
1
20
-10
-5
0
5
10
1
20
0
5
10 15 20
Mo
-15
W
0
10 15 -20 -15 -10 -5
E (eV) Fig. 6. DOS of various MSi2 with the C11b structure. Notice the slight p–d peak shifts for disilicides of the Cr-group. The emboldened metals correspond to experimentally observed stable C11b MSi2.
p-d peak positi on (eV)
-1
-2
3d
MSi2 C40 -3
4d 5d 2
3
4
-1
3d 4d 5d
-2
MSi2 C11b -3 2
3
4
The theoretically predicted lattice parameters are in excellent agreement with available experimental measurements. There is also very good overall agreement between the calculated and experimental heats of formation. The results are very useful for developing full thermodynamic databases using the CALPHAD approach, where lattice stabilities of counter-phases are critical for accurate thermodynamic descriptions of multi-component systems.
Acknowledgements
d-electrons Fig. 7. The energy of the major p–d peak for the C40 and C11b structures, plotted against the number of d-electrons per metal atom.
The author thanks Dr. D.A. Pankhurst, Dr. D. Nguyen-Manh and Prof. D.G. Pettifor for useful discussions.
4. Conclusions Structural stabilities of disilicides of the group IVB to VIB transition metal elements are studied using the ab initio density functional theory, resulting in a full database for the MSi2 lattice stabilities of the stable and metastable structures.
References [1] Jackson MR, Bewlay BP, Rowe RG, Skelly DW, Lipsitt HA. JOM 1996;48:39. [2] Subramanian PR, Mendiratta MG, Dimiduk DM. JOM 1996;48:33.
3736
G. Shao / Acta Materialia 53 (2005) 3729–3736
[3] Zhao JC, Peluso LA, Jackson MR, Tan L. J Alloy Compd 2003;360:183. [4] Zhao JC, Jackson MR, Peluso LA. Acta Mater 2003;51:6395. [5] Liu YQ, Shao G, Tsakiropoulos P. Intermetallics 2000;8:593. [6] Liu YQ, Shao G, Tsakiropoulos P. Intermetallics 2001;9:125. [7] Shao G. Intermetallics 2004;12(6):655. [8] Shao G. Intermetallics 2005;13:69. [9] Saunders N, Miodownik AP. CALPHAD – A comprehensive guide. Oxford: Pergamon Press; 1998. [10] Villars P, Calvert LD, editors. PearsonÕs handbook of crystallographic data for intermetallic phases. 2nd ed. Materials Park (OH): ASM International; 1991. [11] Hohenberg P, Kohn W. Phys Rev 1964;B136:864. [12] Kohn W, Sham LJ. Phys Rev 1965;A140:1133. [13] Perdew JP, Wang Y. Phys Rev 1992;B45:13244. [14] Perdew JP, Burke K, Ernzerhof M. Phys Rev Lett 1996;77:3865.
[15] Segall MD, Lindan PJD, Probert MJ, Pickard CJ, Hasnip PJ, Clark SJ, et al. J Phys: Condens Matter 2002;14:2717. [16] Delley B. J Chem Phys 2000;113:7756. [17] Pankhurst A, Nguyen-Manh D, Pettifor DG. Phys Rev 2004;B69:075113. [18] Carlsson AE, Meschter PJ. J Mater Res 1991;6:1512. [19] Boer FR de, Boom R, Mattens VCM, Miedema AR, Niessen AK. Cohesion in metals, tranistion metal alloys. Amsterdam: NorthHolland; 1988. [20] Meschel SV, Kleppa OJ. J Alloy Compd 1998;267:128. [21] Meschel SV, Kleppa OJ. J Alloy Compd 1998;274:193. [22] Meschel SV, Kleppa OJ. J Alloy Compd 1998;280:231. [23] Chart TG. High Temp–High Press 1973;5:241. [24] Pankhurst DA, Yuan Z, Nguyen-Manh D, Abel ML, Shao G, Watts JF, et al. Phys Rev 2005;B71:075114. [25] Monkhorst HJ, Pack JD. Phys Rev 1977;B16:1748.