First-principles study of ground state properties of zirconium dihydride
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J u l First-principles study of ground state properties of ZrH Peng Zhang,
1, 2
Bao-Tian Wang,
3, 2
Chao-Hui He, and Ping Zhang
2, 4, ∗ Department of Nuclear Science and Technology,Xi’an Jiaotong University, Xi’an 710049, People’s Republic of China LCP, Institute of Applied Physics and Computational Mathematics,Beijing 100088, People’s Republic of China Institute of Theoretical Physics and Department of Physics,Shanxi University, Taiyuan 030006, People’s Republic of China Center for Applied Physics and Technology,Peking University, Beijing 100871, People’s Republic of China
Abstract
Structural, mechanical, electronic, and thermodynamic properties of fluorite and tetragonalphases of ZrH are systematically studied by employing the density functional theory within gen-eralized gradient approximation. The existence of the bistable structure for ZrH is mainly due tothe tetragonal distortions. And our calculated lattice constants for the stable face-centered tetrag-onal (fct) phase with c/a =0.885 are consistent well with experiments. Through calculating elasticconstants, the mechanically unstable characters of face-centered cubic (fcc) phase and fct structurewith c/a =1.111 are predicted. As for fct0.885 structure, our calculated elastic constants explicitlyindicate that it is mechanically stable. Elastic moduli, Poisson’s ratio, and Debye temperature arederived from elastic constants. After analyzing total and partial densities of states and valenceelectron charge distribution, we conclude that the Zr − H bonds in ZrH exhibit weak covalent fea-ture. But the ionic property is evident with about 1.5 electrons transferring from each Zr atom toH. Phonon spectrum results indicate that fct0.885 and fct1.111 structures are dynamically stable,while the fcc structure is unstable. PACS numbers: 61.50.Ah, 62.20.Dc, 71.15.Mb, 63.20.dk ∗ Author to whom correspondence should be addressed. E-mail: zhang [email protected] . INTRODUCTION Many transition metals react readily with hydrogen to form stable metal hydrides [1].The common metal hydrides are technologically attractive materials due to their abilityto store high densities of hydrogen safely. In the case of zirconium hydrides, another ap-pealing interest arises in nuclear industry, essentially used as a neutron moderator [2] andfuel rod cladding materials [3] in nuclear reactors. Additionally, in the fusion technologyapplications, zirconium hydrides are unavoidably precipitated in the excess of the hydrogensolubility, which may adversely affect several mechanical and thermal properties in somecrystal orientations [4] involving different stable and metastable hydride phases which maylead to a significant embrittlement. Furthermore, zirconium hydrides are also served asa part of the promising new type of actinide hydride fuels, such as U-(Th-Np-Am)-Zr-H[5]. Consequently, great attentions are needed to focus on the basic scientific research ofzirconium hydrides.Since 1952, there has occurred in literature large amounts of experimental investigationson zirconium hydrides. The electronic structures of ZrH x (1.63 ≤ x ≤ et al. [7] in 1990. And the transitionbetween the cubic and the tetragonal phases has been suggested by Ducastelle et al. [8] to beof Jahn-Teller type too. From theoretical point of view, the phase transition between cubicand tetragonal phases has been investigated by Switendick [9], where the non-self-consistentaugmented plane wave (APW) method was employed. Most importantly, the study of theelectronic structure and energetics for the tetragonal distortion in ZrH has been extensivelyperformed. Using the pseudopotential method with both local density approximation (LDA)and general gradient approximation (GGA), Ackland [10] obtained a double minimum in theenergy surface, one minimum for c/a < c/a > c/a > et al. [18] studiedthe accurate total energy and the electronic structure characteristic of ZrH by using LDAand GGA, respectively, within the full potential linear augmented plane wave (FP-LAPW)method. Both of these works obtained the correct ground state for c/a < . Moreover, althoughthe electronic properties as well as the chemical bonding in zirconium dihydride have beencalculated [17], the study of the bonding nature of Zr-H involving its mixed ionic/covalentcharacter is still lacking. As a consequence, these facts inhibit deep understanding of thezirconium dihydride. Motivated by these observations, in this paper, we present a first-principles study by calculating the structural, electronic, mechanical, and thermodynamicalproperties of zirconium dihydride. Our calculated results show that the bistable structureof ZrH is the fct phase at c/a =0.885 and c/a =1.111, respectively. The mechanical and dy-namical stability of fct and fcc phases of ZrH are carefully analyzed. In addition, throughBader analysis [19, 20] we find that about 1.5 electrons transfer from each Zr atom to Hatom for ZrH . II. COMPUTATIONAL METHOD
Our total energy calculations are self-consistently carried out using density functionaltheory (DFT) as implemented in Vienna ab initio simulation package (VASP) [21], which isbased on pseudopotentials and plane wave basis functions. All electron projected augmentedwave (PAW) method of Bl¨ochl [22] is applied in VASP with the frozen core approximation.The generalized gradient approximation (GGA) introduced by Perdew, Burke, and Ernz-erhof (PBE) [23] is employed to evaluate the electron exchange and correlation potential.Zirconium 4 s p d s and hydrogen 1 s electrons are treated as valence electrons. Theintegration over the Brillouin Zone (BZ) is performed with a grid of special k point-mesh de-termined according to the Monkhorst-Pack scheme [24]. After convergence test, 12 × × k point-mesh is chosen to make sure the total energy difference less than 1 meV per primitivecell. When exploring the electronic structure, a finer k point-mesh 24 × ×
24 is preferredto generate a high quality charge density. The expansion of the valence electron plane wavefunctions can be truncated at the cutoff energy of 650 eV, which has been tested to be fully3
IG. 1: Tetragonal unit cell in space group I /mmm (a) and cubic unit cell for ZrH in spacegroup F m ¯3 m (b) with larger cyan spheres for Zr atoms and the smaller white H. converged with respect to total energy. Relaxation procedures at ground state are carriedout according to the Methfessel-Paxton scheme [25] with a smearing width of 0.1 eV. Thegeometry relaxation is considered to be completed when the Hellmann-Feynman forces onall atoms were less than 0.001 eV/˚A. And the accurate total energy calculations are per-formed by means of the linear tetrahedron method with the Bl¨ochl’s corrections [26]. Theself-consistent convergence of the total energy calculation is set to 10 − eV. III. RESULTSA. Atomic structure and mechanical properties
At ambient condition, the stable zirconium dihydride crystallizes in a fct structure withspace group I mmm (No. 139). In its unit cell, there are two ZrH formula units with Zrand H atoms in 2 a :(0,0,0) and 4 d :(0, , ) sites, respectively [see Fig. 1(a)]. Each Zr atomis surrounded by eight H atoms forming a tetragonal and each H connects with four Zratoms to build a tetrahedron. The present optimized equilibrium lattice parameters ( a and c ) obtained by fitting the energy-volume data in the third-order Birch-Murnaghan equationof states (EOS) [27] are 5.030 ˚A and 4.414 ˚A (see Table I), in good agreement with theexperimental values. In addition, the fcc fluorite type structure with space group Fm ¯3 m (No. 225) is considered as the metastable phase for ZrH . The cubic unit cell is composedof four ZrH formula units with the Zr and H atoms in 4 a :(0,0,0) and 8 c :( , , ) sites,respectively [see Fig. 1(b)]. Here, each Zr atom is surrounded by eight H atoms forming acube and each H atom connects with four Zr atoms to build a tetrahedron. Our optimizedequilibrium lattice constant a for fcc ZrH is 4.823 ˚A. This value, although no experimental4 ABLE I: Calculated lattice constants, bulk modulus B , and pressure derivative of the bulkmodulus B ′ at ground state from EOS fitting. For comparison, available experimental values andother theoretical results are also listed.Compounds Property Present work Previous calculation Experimentfcc ZrH a (˚A) 4.823 4.817 a , 4.804 b - B (GPa) 133 136 a ,152 b - B ′ a (˚A) 5.030 5.021 a , 5.008 b , 5.000 c d , 4.975 e , 4.982 f c (˚A) 4.414 4.432 a , 4.419 b , 4.450 c d , 4.447 e , 4.449 f B (GPa) 127 - - B ′ a Reference [18], b Reference[17], c Reference [10], d Reference [11], e Reference[16], f Reference[12]. -68.0-67.9-67.8-67.7 E ne r g y ( e V / p r i m i t i v e c e ll ) c/a ZrH FIG. 2: Total energy as a function of the c/a ratio in the fct phase for ZrH . structural values are available for comparison, is in good agreement with previous theoreticalwork in the literature (see Table I).In present study, the stability of ZrH has been investigated by calculating the totalenergy as a function of c/a at the optimized constant volume of fcc structure, as illustratedin Fig. 2. The energy curve displays two local minima, one with c/a < c/a > c/a =0.885 has the lowest energy and corresponds tofct phase observed experimentally. The fct-fcc structure energy barrier obtained from thiscalculation is 0.1048 eV, in good agreement with the Ackland’s results [10] but much largerthan the FP-LAPW calculation results [17, 18].5 ABLE II: The strain combinations in the strain tensor [Eq. (1)] to calculate the elastic constantsof cubic and tetragonal ZrH .Phase Strain Parameters (unlisted e i =0) ∆ E/V in O( δ ) fcc ǫ e = e = e = δ ( C + 3 C ) δ ǫ e = e = δ ( C + C ) δ ǫ e = e = e = δ C δ fct ǫ e = e = e = δ ( C + C + 2 C + C ) δ ǫ e = e = δ ( C + 2 C + C ) δ ǫ e = e = δ ( C + C ) δ ǫ e = δ C δ ǫ e = e = δ C δ ǫ e = δ C δ Elastic constants can measure the resistance and mechanical features of crystal to externalstress or pressure, thus may evaluate the stability of crystals against elastic deformation.For small strain ε , Hooke’s law is valid. Therefore, we can enforce it onto the equilibriumlattice, determine the resulting change in the total energy, and from this information deducethe elastic constants. The crystal total energy E ( V, ǫ ) can be expanded as a Taylor series[28], E ( V, ǫ ) = E ( V ,
0) + V X i =1 σ i e i + V X i,j =1 C ij e i e j + O ( { e i } ) , (1)where E ( V ,
0) is the energy of the unstrained system with the equilibrium volume V , ǫ isthe strain tensor which has matrix elements ε ij ( i, j =1, 2, and 3) defined by ε ij = e e e e e e e e e , (2)and C ij are the elastic constants. Such a strain transforms the three lattice vectors definingthe unstrained Bravais lattice { a k , k =1, 2, and3 } to the strained vectors { a ′ k } [29] as definedby a ′ k = ( I + ε ) a k , (3)6 ABLE III: Calculated elastic constants in units of GPa for ZrH at zero pressure.Phase C C C C C C fcc 82.6 159.7 -19.5fct0.885 165.6 140.9 30.5 106.8 145.5 60.6fct1.111 125.7 145.5 30.9 115.0 190.6 42.0 where I is a unit 3 × a k or a ′ k is a 3 × C , C , and C ,which can be calculated through a proper choice of the set of strains { e i , i = 1 , ..., } (seeTable II). As for the fct phase, the six independent elastic constants, i.e., C , C , C , C , C , and C , can be obtained from six different strains listed in Table II. To avoid theinfluence of high order terms on the estimated elastic constants, we have used very smallstrains, i.e. within ± δ is varied insteps of 0.006 from δ = − E ( V, δ ) at these strain stepsare calculated and fitted through the strains with the corresponding parabolic equations of∆
E/V to yield the required second-order elastic constants.Our calculated elastic constants for fcc ZrH and fct ZrH at c/a =0.885 and c/a =1.111are collected in Table III. Obviously, fct ZrH at c/a =0.885 is mechanically stable due tothe fact that its elastic constants satisfy the following mechanical stability criteria [30] oftetragonal structure: C > , C > , C > , C > , ( C − C ) > , ( C + C − C ) > , [2( C + C ) + C + 4 C ] > . (4)Nevertheless, the fct phase at c/a =1.111 is mechanically unstable in light of the fact that C is smaller than C . As for the fcc structure of ZrH , it is also mechanically unstablebecause of the fact that its elastic constants cannot meet the following mechanical stabilitycriteria [30] of cubic structure: C > , C > , C > | C | , ( C + 2 C ) > . (5)After obtaining elastic constants, bulk modulus B and shear modulus G for fct phaseof ZrH at c/a =0.885 can be calculated from the Voigt-Reuss-Hill (VRH) approximations7 ABLE IV: Calculated bulk modulus, shear modulus, Young’s modulus, Poisson’s ratio, density,transverse sound velocity, longitudinal sound velocity, average sound velocity, and Debye temper-ature of the fct phase of ZrH at c/a =0.885. B (GPa) G (GPa) E (GPa) υ ρ (g/cm ) v t (km/s) v l (km/s) v m (km/s) θ D (K)130 29 80 0.397 5.5201 2.2823 5.5197 2.5825 364.9 [31–33]. Then the Young’s modulus E and Poisson’s ratio υ are calculated through E =9 BG/ (3 B + G ) and υ = (3 B − G ) / [2(3 B + G )]. The Debye temperature ( θ D ) is obtainedusing the relation [34] θ D = hk B (cid:20) n π (cid:18) N A ρM (cid:19)(cid:21) / v m , (6)where h and k B are Planck’s and Boltzmann’s constants, respectively, N A is Avogadro’snumber, ρ is the density, M is the molecular weight, n is the number of atoms in themolecule, and v m is the average wave velocity. The average wave velocity in the polycrys-talline materials is approximately given as v m = (cid:20)
13 ( 2 v t + 1 v l ) (cid:21) − / , (7)where v t = p G/ρ and v l = p (3 B + 4 G ) / ρ are the transverse and longitudinal elastic wavevelocity of the polycrystalline materials, respectively.All the calculated results for fct0.885 are collected in Table IV. Results of fcc and fct1.111can not be obtained due to their mechanically unstable nature. Note that we have alsocalculated the bulk modulus B by EOS fitting. The derived bulk modulus for fct0.885 turnsout to be very close to the one obtained from the above VRH approximation, which againindicates that our calculations are consistant and reliable. It is well known that the shearmodulus G represents the resistance to plastic deformation, while the bulk modulus B canrepresent the resistance to fracture. A high (low) B/G value is responsible for the ductility(brittleness) of polycrystalline materials. The critical value to separate ductile and brittlematerials is about 1.75. Using the calculated values of bulk modulus B and shear modulus G for fct0.885, the B/G value of 4.48 can be obtained. For element Zr, the
B/G value of2.63 can be derived from the elastic data in Ref. [35]. Therefore, fct phase of ZrH is ratherductile and its ductility is more predominant than that of element Zr. As for the Poisson’s8
10 -8 -6 -4 -2 0 2 4012 D O S ( s t a t e s / e v Z r a t o m ) Energy(eV) (a) fcc -10 -8 -6 -4 -2 0 2 4 (b) fct0.885 -10 -8 -6 -4 -2 0 2 4 total
Zr-s
Zr-p
Zr-d
H-s (c) fct1.111
FIG. 3: Total and orbital-resolved local densities of states for (a) fcc, (b) fct0.885, and (c) fct1.111structures of ZrH . The Fermi energy level is set at zero. ratio, the value is well within the range from 0.25 to 0.45 for typical metals. Unfortunately,no reliable experimental and theoretical results concerning on the mechanical properties inthe literature are available for comparison. We hope that our calculated elastic constantsand elastic moduli can be illustrative in the realistic application of the mechanical data forZrH . B. Electronic structure and charge distribution
Basically, all the macroscopical properties of materials, such as hardness, elasticity, andconductivity, originate from their electronic structure properties as well as the nature of thechemical bonding. Therefore, it is necessary to perform the electronic structure analysis ofZrH . The calculated total densities of states (DOS) and the orbital-resolved partial densitiesof states (PDOS) of fcc ZrH , fct ZrH at c/a =0.885 and c/a =1.111, are shown in Fig. 3.Wholly, the occupation properties of ZrH in fct and fcc phases are similar. The pseudogapbelow the Fermi level, located at around − s orbitaland Zr 4 d and Zr 4 p orbitals. States in this region have critical contribution to the Zr-H andH-H bonding state. However, from the pseudogap to the Fermi level, H 1 s states contributevery little to the total DOS, while the Zr 4 d states dominate in this energy range. The factthat the DOS occupation has obvious sharp peak at the Fermi level for fcc ZrH implies thatthe cubic phase should be unstable. In contrast, a strong reduction in the DOS at the Fermilevel is observed in the two fct structures. This evident difference between fcc and fct phase9 Zr H Zr
H Zr (a) (b) (c)
FIG. 4: Valence charge density of (a) fcc phase in (110) plane, (b) fct0.885 in (100) plane, and (c)fct1.111 in (100) plane. The contour lines are drawn from 0.0 to 1.0 at 0.05 e/˚A intervals. in the vicinity of the Fermi level is caused by the splitting of the band in this region. Asa consequence, the Fermi level in the two tetragonal phases is not located at a peak of theDOS but at a local minimum, which leads to a reduction in the electronic contribution to thetotal energy. This behavior has been referred to as a Jahn-Teller mechanism, as indicatedby previous experiments and theoretical works [6, 12, 17, 18].To analyze the ionic/covalent character of zirconium dihydride, in the following we willcarefully investigate the valence charge density distribution. The calculated valence chargedensity maps of the fct ZrH at c/a =0.885 and c/a =1.111 in (100) plane and fcc ZrH in(110) plane are plotted in Fig. 4. Obviously, the charge densities around Zr and H ions areall near spherical distribution with slightly deformed towards the direction to their nearestneighboring atoms. Covalent bridges, which represent the bonding nature of Zr − H bonds,is clearly shown in Fig. 4. In Fig. 5, we plot the line charge density distribution alongthe nearest Zr − H bonds. Clearly, three line charge density curves vary in almost the sameway along the Zr − H bonds. One can find a minimum value of charge density for each bondat around the bridge locus (indicated by the arrow in Fig. 5). These minimum values arelisted in Table V. Although these values are much smaller than 0.7 e/˚A for Si covalent bond,they are prominently higher than 0.05 e/˚A for Na − Cl bond in typical ionic crystal NaCl.Therefore, there are weak but clear weights in Zr − H bonds in both fcc and fct phases ofZrH .In addition, we have performed the Bader analysis [19, 20] for the three typical ZrH cells. The charge ( Q B ) enclosed within the Bader volume ( V B ) is a good approximation10 .0 0.5 1.0 1.5 2.001234567 oo HZr
Distance (A) C ha r ge D en s i t y ( e / A ) fcc fct0.885 fct1.111 FIG. 5: The line charge density distribution between Zr atom and the nearest neighbor H atomfor fcc, fct0.885, and fct1.111 structures of ZrH . to the total electronic charge of an atom. In this work, the default charge density gridsfor one primitive cell are 56 × ×
56, 56 × ×
56, 60 × ×
60 for fcc, fct0.885, and fct1.111cells, respectively. To check the precision, the charge density distributions are calculatedwith a series of n times finer grids ( n =2,3,4,5,6,7). The deviation of the effective chargebetween the six and the seven times finer grids is less than 0.2%. Thus we perform thecharge density calculations using the seven times finer grid (392 × × × × × × are almost equal to each other. This shows similar ionic character, througha flux of charge (about 1.5 electrons for each Zr atom) from cations towards anions, forzirconium dihydride in their stable phase and metastable phases. The slight difference incharge distribution is mainly due to the structure distortion. Thus the ionicity of Zr-Hbonds for ZrH is also evident either for its fcc structure or fct structures. C. Phonon dispersion curves
Phonon frequencies of the crystalline structure is one of the basic aspects when consideringthe phase stability, phase transformations, and thermodynamics of crystalline materials.11
ABLE V: Calculated charges Q B and volumes V B (˚A ) according to Bader partitioning as well asthe Zr − H distance (˚A) and relevant minimum values of charge densities (e/˚A ) along the Zr − Hbonds for ZrH .Compounds Q B (Zr) Q B (H) V B (Zr) V B (H) Zr-H distance Charge density min fcc 10.484 1.758 14.809 6.618 2.088 0.320fct0.885 10.512 1.744 14.953 6.544 2.095 0.316fct1.111 10.514 1.743 14.954 6.545 2.094 0.317 Employing the Hellmann-Feynman theorem [36, 37] and the direct method [38], we havecalculated the phonon curves along some high symmetry directions in the BZ, together withthe phonon DOS. For the phonon dispersion calculation, we use the 2 × × and the 4 × × k -point meshfor the BZ integration. In order to calculate the Hellmann-Feynman forces, we displace fourand eight atoms, respectively, for fcc and fct ZrH from their equilibrium positions and theamplitude of all the displacements is 0.03˚A. The calculated phonon dispersion curves alongthe Γ- X - K -Γ- L - X - W - L directions for fcc ZrH and along the Γ- N - P - X -Γ- Z directions forfct ZrH are displayed in Fig. 6. For both fcc and fct ZrH , there are only three atoms intheir primitive cell. Therefore, nine phonon modes exist in the dispersion relations. Due tothe fact that zirconium is much heavier than hydrogen, the vibration frequency of zirconiumatoms is apparently lower than that of hydrogen atoms. As a result, evident gap between theoptic modes and the acoustic branches exits and the phonon DOS of fcc (fct) ZrH can beviewed as two parts. One is the part lower than 7.4 (7.1) THz, where the main contributioncomes from the zirconium sublattice, while the other part range from 30.7 to 37.9 (29.6 to38.5) THz is dominated by the dynamics of the light hydrogen atoms. As shown by Fig.6(b) and 6(c), all frequencies are positive for the two fct structures. This assures that thetwo fct phases are both dynamically stable. In contrast, one can see from Fig. 6(a) thatthe fcc phase of ZrH is dynamically unstable. The transverse acoustic (TA) mode closeto Γ point becomes imaginary along the Γ- K and Γ- L directions. The minimum of the TAbranch locates at the Γ- K direction. Therefore, the fcc to fct phase transition probably canoccur along the [101] direction. This kind of dynamical instability of fcc phase compared tothe stable fct phase is well consistent with our above-discussed mechanical stability analysis12 (b) (c) L F r equen cy ( T H z ) X K L X W (a)
H Zr
ZXPN F r equen cy ( T H z ) H Zr F r equen cy ( T H z ) ZXPN
H Zr
FIG. 6: Calculated phonon dispersion curves along the high-symmetry directions (left panel) andthe corresponding PDOS (right panel) at ground state for (a) fcc ZrH , (b) fct ZrH at c/a =0.885,and (c) fct ZrH at c/a =1.111. of ZrH . IV. CONCLUSION
In summary, we have investigated the structural, mechanical, electronic, and thermody-namic properties of ZrH in its stable and metastable phases by means of first-principlesDFT-GGA method. Our optimized structural parameters are well consistent with previousexperiments and theoretical calculations. Tetragonal distortion leads to bistable structureswith c/a > c/a <
1, the latter of which corresponds to the experimentally observedlow-temperature phase. Our total energy calculations illustrate that the most stable phaseof ZrH is the fct structure with c/a =0.885. Elastic constants, various moduli, Poisson’sratio, and Debye temperature are calculated for fcc and two fct phases of ZrH . Mechani-cally unstable nature for fcc and fct1.111 structures are predicted. For fct0.885, mechanicalstability and dynamical stability are firstly predicted by our elastic constants results andphonon dispersion calculation. The occupation characters of electronic orbital also accordwell with experiments. Through Bader analysis, we have found that the Zr − H bonds ex-hibit weak covalent character while the ionic property is dominant with about 1.5 electronstransferring from each Zr atom to H. Our calculated phonon curves of fcc ZrH have shownthat the TA mode becomes imaginary close to Γ point along the [101] direction.13 . ACKNOWLEDGMENTS This work was partially supported by NSFC under Grants No. 90921003 and No.60776063. [1] Y. Fukai, Metal Hydrogen System. Basics Bulk Properties, Springer,Berlin,2005.[2] P.W. Bickel, T.G. Berlincourt, Phys. Rev. B 2 (1970) 4807.[3] L. Holliger, A. Legris, R. Besson, Phys. Rev. B 80 (2009) 094111.[4] C.E. Ellis, J. Nucl. Mater. 28 (1968) 129.[5] J. Huang, B. Tsuchiya, K. Konashi, M. Yamawaki, J. Nucl. Sci. Techol. 37 (2000) 887.[6] J.H. Weaver, D.J. Peterman, D.T. Peterson, A. Franciosi, Phys. Rev. B 23 (1981) 1692.[7] E. Zuzek, J.P. Abriata, A.S. Martin, F.D. Manchester, Bull. Alloy Phase Diagrams 11 (1990)385.[8] F. Ducastelle, R. Caudron, P. Costa J. Physique 31 (1970) 57.[9] A.C. Switendick, J. Less-Common Met. 101 (1984) 191.[10] G.J. Ackland, Phys. Rev. Lett. 80 (1998) 2233.[11] H.L. Yakel, Acta Crystallogr. 11 (1958) 46.[12] R.C. Bowman Jr., E.L. Venturini, B.D. Craft, A. Attalla, D.B. Sullenger, Phys. Rev. B 27(1983) 1474.[13] R.C. Bowman Jr., B.D. Craft, J. Phys. C 17 (1984) L477.[14] R.C. Bowman Jr., B.D. Craft, J.S. Cantrell, E.L. Venturini, Phys. Rev. B 31 (1985) 5604.[15] O.J. Zogal, B. Nowak, K. Niedzwiedz, Solid State Commun. 80 (1991) 601.[16] K. Niedzwiedz, B. Nowak, O.J. Zogal, J. Alloys Compd. 194 (1993) 47.[17] W. Wolf, P. Herzig, J. Phys.: Condens. Matter 12 (2000) 4535.[18] R. Quijano, R. Coss, D.J. Singh, Phys. Rev. B 80 (2009) 184103.[19] W.F.W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, New York,1990.[20] W. Tang, E. Sanville, G. Henkelman, J. Phys.: Condens. Matter 21 (2009) 084204.[21] G. Kresse, J. Furthm¨uller, Phys. Rev. B 54 (1996) 11169.[22] P.E. Bl¨ochl, Phys. Rev. B 50 (1994) 17953.
23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865.[24] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1972) 5188.[25] M. Methfessel, A.T. Paxton, Phys. Rev. B 40 (1989) 3616.[26] P.E. Bl¨ochl, O. Jepsen, O.K. Andersen, Phys. Rev. B 49 (1994) 16223.[27] F. Birch, Phys. Rev. 71 (1947) 809.[28] J.F. Nye, Physical Properties of Crystals, Their Representation by Tensors and Matrices,Oxford Press, Chap.VIII, 1957.[29] J.H Westbrook, R.L Fleischer (Ed.), Intermetallic Compounds: Principles and Practice, VolI:principles, Wiley, London, 1995, Chap. 9, p. 195.[30] J.F. Nye, Physical Properties of Crystals, Oxford University Press, 1985.[31] W. Voigt, Lehrburch der Kristallphysik, Teubner,Leipzig, 1928.[32] A. Reuss, Z. Angew, Math. Mech. 9 (1929) 49.[33] R. Hill, Phys. Soc. London 65 (1952) 350.[34] O.L. Anderson, J. Phys. Chem. Solids 24 (1963) 909.[35] W. Liu, B. Li, L. Wang, J. Zhang, and Y. Zhao, J.Appl. Phys. 104 (2008) 076102.[36] O.H. Nielsen, R.M. Martin, Phys. Rev. B 32 (1985) 3780[37] O.H. Nielsen, R.M. Martin, Phys. Rev. B 35 (1987) 9308[38] K. Parlinski, Z.Q. Li, Y. Kawazoe, Phys. Rev. Lett. 78 (1997) 4063.23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865.[24] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1972) 5188.[25] M. Methfessel, A.T. Paxton, Phys. Rev. B 40 (1989) 3616.[26] P.E. Bl¨ochl, O. Jepsen, O.K. Andersen, Phys. Rev. B 49 (1994) 16223.[27] F. Birch, Phys. Rev. 71 (1947) 809.[28] J.F. Nye, Physical Properties of Crystals, Their Representation by Tensors and Matrices,Oxford Press, Chap.VIII, 1957.[29] J.H Westbrook, R.L Fleischer (Ed.), Intermetallic Compounds: Principles and Practice, VolI:principles, Wiley, London, 1995, Chap. 9, p. 195.[30] J.F. Nye, Physical Properties of Crystals, Oxford University Press, 1985.[31] W. Voigt, Lehrburch der Kristallphysik, Teubner,Leipzig, 1928.[32] A. Reuss, Z. Angew, Math. Mech. 9 (1929) 49.[33] R. Hill, Phys. Soc. London 65 (1952) 350.[34] O.L. Anderson, J. Phys. Chem. Solids 24 (1963) 909.[35] W. Liu, B. Li, L. Wang, J. Zhang, and Y. Zhao, J.Appl. Phys. 104 (2008) 076102.[36] O.H. Nielsen, R.M. Martin, Phys. Rev. B 32 (1985) 3780[37] O.H. Nielsen, R.M. Martin, Phys. Rev. B 35 (1987) 9308[38] K. Parlinski, Z.Q. Li, Y. Kawazoe, Phys. Rev. Lett. 78 (1997) 4063.