Computational Materials Science 48 (2010) 859–865
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Thermo-elastic and lattice dynamical properties of Rh3Hf compound G. Surucu a,*, K. Colakoglu a, E. Deligoz b, N. Korozlu a, H. Ozisik b a b
Gazi University, Department Of Physics, Teknikokullar, 06500 Ankara, Turkey Aksaray University, Department Of Physics, 68100 Aksaray, Turkey
a r t i c l e
i n f o
Article history: Received 18 November 2009 Received in revised form 31 March 2010 Accepted 8 April 2010
Keywords: Ab initio calculations Electronic structure Phonons
a b s t r a c t The lattice dynamical calculations have been performed on the L12-type (space number 221) of intermetallic compound Rh3Hf using the ab initio density-functional theory within the local density approximation (LDA) and the generalized gradient approximation (GGA). Beside the basic physical parameters such as lattice constant, bulk modulus, elastic constants, shear modulus, Young’s modulus, and Poison’s ratio; the phonon dispersion curves and corresponding one-phonon density of states (DOS) are also calculated for the same compound. The temperature and pressure variations of the volume, bulk modulus, thermal expansion coefficient, heat capacity, and Debye temperature in a wide pressure (0–180 GPa) and temperature (0–2000 K) ranges are presented in this study. In particular, our structural parameters (the lattice constant and bulk modulus) are consistent with the available experimental and other theoretical data. Ó 2010 Elsevier B.V. All rights reserved.
1. Introduction Ni-based super alloys are important materials for various technological applications, especially for gas turbine engines and aircraft engines [1] due to their high-temperature physical properties. The recent studies [2–14] on these materials have shed light on the thermo-mechanical properties of platinum metal-base alloys, such as Ir-and Rh-based alloys, because it is thought that these are more promising materials than Ni-base superalloys due to their high-melting temperatures, good high-temperature strengths and good oxidation resistances [13,14]. So, the results of the detailed calculations on the mechanical and thermo-elastic properties of Rh3Hf are not only important for academic interest but also serve as a valuable guide for high-temperature structural applications. Crystallographically, the rhodium based L12 intermetallic compound Rh3Hf crystallizes in Cu3Au type structure with the Pm3m space group symmetry (space number 221) [15]. The constituent atom Hf and the Rh atom are localed at (0, 0, 0) and (0, 1/2, 1/2) positions, respectively. There are only a few theoretical [4,5] and experimental [3,14] studies dealing with the structural, elastic, and electronic properties of Rh3Hf in the literature. Chen et al. [4] investigated the elastic and mechanical properties of this compound based on the ab initio density-functional theory within the GGA. Rajagopalan and Sundareswari [5] reported the structural and electronic properties using the self-consistent tight binding linear muffin tin orbital
* Corresponding author. Tel.: +90 312 2021233; fax: +90 312 2122279. E-mail address:
[email protected] (G. Surucu). 0927-0256/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2010.04.008
(TB-LMTO) method for Rh3Hf. Experimentally, Terada et al. [3] measured the thermal conductivity and thermal expansion as a function of the temperature for the same compounds in the temperature range of 300–1100 K. Yamabe et al. [14] investigated the microstructural evolution and high-temperature strengths of Rh-based alloys. Up to now –to our best knowledge–, no systematic study on the vibrational, and thermodynamical properties of Rh3Hf has been reported. Therefore, we have aimed to provide some additional information about the physical properties of Rh3Hf to the existing data using the ab initio total energy calculations. Especially, we focus our attention on the important bulk properties such as mechanical, thermodynamical, and the lattice dynamical behavior of this material. We have also discussed the temperature effect on the structural parameters; bulk modulus, thermal expansion coefficient, heat capacity, and Debye temperature in a wide pressure (from 0 to 180 GPa) and temperature (from 0 to 2000 K) ranges.
2. Method of calculation All calculations have been carried out by using the Vienna ab initio simulation package (VASP) [16–19] based on the densityfunctional theory (DFT). The electron–ion interaction was considered in the form of the projector-augmented-wave (PAW) method with the plane wave up to energy of 400 eV [18,20]. This cut-off energy was found to be adequate to study the structural and lattice dynamical properties. We have not observed statistically significant changes in the key parameters when the energy cut-off is increased from 300 to 500 eV.
860
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865
We have used both the local density approximation (LDA) parameterized by Perdew and Zunger [21] and the generalized gradient approximation (GGA) parameterized by Perdew–Burke–Ernzerhof (PBE) [22] for the exchange and correlation terms in the electron–electron interaction for k-space summation which was 12 12 12 Monkhorst and Pack grid of k-points [23]. The relaxation on the atomic coordinates has been carefully checked in each step of calculations for structural and elastic properties. The thermodynamical properties of Rh3Hf are computed using the quasi-harmonic Debye model [25] in which the non-equilibrium Gibbs function G(V; P, T) takes the form:
G ðV; P; TÞ ¼ EðVÞ þ PV þ Avib ½hðVÞ; T:
ð1Þ
In Eq. (1), E(V) is the total energy for per unit cell of Rh3Hf, PV corresponds to the constant hydrostatic pressure condition, h (V) is the Debye temperature and AVib is the vibrational Helmholtz free energy and can be written as [26–30]:
AVib ðh; TÞ ¼ nkT
9h h h þ 3 ln 1 eT D 8T T
ð2Þ
where n is the number of atoms per formula unit, D Th is the Debye integral. For an isotropic solid, h is expressed as [28]
hD ¼
rffiffiffiffiffi i1=3 h h Bs f ðrÞ 6pV 1=2 n k M
ð3Þ
where M is the molecular mass per unit cell and BS is the adiabatic bulk modulus, which is given approximately by the static compressibility [30]: 2
Bs BðVÞ ¼ V
d EðVÞ dV
2
ð4Þ
The values of f ðrÞ are taken from Refs. [28,29], and for the Poisson’s ratio our calculated value of 0.28 is used. n and M are taken to be 4 and 281.3 a.u, respectively. Therefore, the non-equilibrium Gibbs function G(V; P, T) as a function of (V; P, T) can be minimized with respect to volume V as:
@G ðV; P; TÞ ¼0 @V P;T
ð5Þ
By solving Eq. (5), one can obtain the thermal equation-of-state (eos) V(P, T). The heat capacity at constant volume CV and the thermal expansion coefficient (a) are given [31] as follows:
h 3h=T h=T C v ¼ 3nk 4D T e 1
ð6Þ
h S ¼ nk 4D 3 ln 1 eh=T T
ð7Þ
a¼
cC v BT V
ð8Þ
Here c represents the Grüneisen parameter and its general expression is given as:
c¼
d ln hðVÞ d ln V
based on the same Murnaghan’s eos and the results are listed in Table 1 along with the experimental [3,6,14] and other theoretical [4,5] values. As expected, the lattice constant obtained from LDA (3.86 Å) is underestimated from the experimental values of 3.91–3.90 Å, and overestimated from the experimental ones, and is in agreement with the other theoretical LDA result of 3.88 Å in Ref. [5]. In the case of GGA, the lattice constant is found to be 3.94 Å which is in agreement with the other theoretical GGA value of 3.95 Å [4]. The bulk modulus is a fundamental physical property of solids and can also be used as a measure of the average bond strengths of atoms for the given crystals [25]. While our bulk modulus obtained from LDA about 20% lower than the other theoretical value given in Ref. [5], the present result obtained from GGA is about 10% higher than the theoretical value given in Ref. [4]. The calculated bulk modulus of Rh3Hf compares with that of Rh3Ti (263 GPa) [12], Rh3V (289 GPa) [12], Rh3Ta (326.3 GPa) [5], Rh3Nb (315.7 GPa) [5], Rh3Zr (277 GPa) [5]. These results show that the Rh3Hf is a highly compressible material. The elastic constants of solids provide a link between the mechanical and dynamical behavior of crystals, and give important information concerning the nature of the forces operating in solids. In particular, they provide information about the stability and stiffness of materials, and their ab initio calculation requires precise methods. Since the forces and the elastic constants are functions of the first-order and second-order derivatives of the potentials, their calculation will provide a further check on the accuracy of the calculation of forces in solids. Here, the elastic constants are computed using both the ‘‘volume-conserving” technique [33] and the ‘‘stress–strain” relations [34], and the results are listed in Table 2 along with the other theoretical values of Ref. [4]. Our results for elastic constants, obtained from both methods are consistent with each other, but especially, our C11 and C12 values significantly higher than those of other theoretical results of Ref. [4]. These deviations may stem from the used lower number of K-POINTS (8 8 8) and energy cut-off (290 eV) values in their work because the correct values of the elastic constants and phonon frequencies are closely related to the number of K-POINTS. The traditional mechanical stability conditions in cubic crystals on the elastic constants are known as C11 – C12 > 0, C11 > 0, C44 > 0, C11 + 2C12 > 0, and C12 < B < C11. Our results for elastic constants in Table 2 satisfy these traditional stability conditions. The Zener anisotropy factor (A) is an indicator of the degree of anisotropy in the solid structures. For a completely isotropic material, the factor A takes the value of 1. When the value of A is smaller or greater than unity, it is a measure of the degree of elastic anisotropy. In the hardness investigations of the polycrystalline materials, Poisson’s ratio ðtÞ, shear modulus (G), and Young’s modulus (E), which are the most interesting elastic properties for applica-
ð9Þ
Table 1 The calculated equilibrium lattice constant (a0), bulk modulus (B), and pressure derivative of bulk modulus (B0 ) for Rh3Hf.
3. Results and discussion 3.1. Structural and elastic properties Firstly, the equilibrium lattice parameters have been computed by minimizing the crystal total energy calculated for different values of lattice constant by means of Murnaghan’s eos [32]. The bulk modulus, and its pressure derivative have also been calculated
a b c d
Material (Rh3Hf)
a0 (Å)
B (GPa)
B0
Present-LDA Present-GGA Theory-LDAa Theory-GGAb Experimentalc Experimentald
3.860 3.940 3.880 3.950 3.911 3.910
259 215 274 204 –
4.56 4.79 –
Ref. Ref. Ref. Ref.
[5]. [4]. [7,8]. [9].
–
861
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865 Table 2 The calculated elastic constants (in GPa unit), Poisson’s ratio (m), Young’s modulus (E), isotropic Shear modulus (G), and Zener anisotropy factor (A) for Rh3Hf.
a
Material (Rh3Hf)
C11
C12
C44
G
Gv
Gr
t
A
E
Volume-conserving-LDA (present) Stress–strain-LDA (present) Average-LDA (present) Volume-conserving- GGA (present) Stress–strain-GGA (present) Average-GGA (present) Stress–strain GGA (theory a)
380 379 380 321 319 320 296
191 199 195 157 163 160 158
168 169 169 142 144 143 140
133 131 132 114 112 113 –
138 137 138 118 117 118 112
128 125 127 110 107 109 99
0.28 0.28 0.28 0.27 0.28 0.28 0.28
1.78 1.88 1.83 1.73 1.84 1.79 –
340 336 338 290 286 288 270
Ref. [4].
tions, are often measured. We have calculated these quantities, in terms of the computed data, using the following relations [35]:
A¼
2C44 C11 C12
ð10Þ
t¼
" # 1 B 23 G 2 B þ 13 G
ð11Þ
9GB G þ 3B
ð12Þ
and
E¼
where G = (GV + GR )/2 is the isotropic shear modulus, GV is Voigt’s shear modulus corresponding to the upper bound of G values, and GR is Reuss’s shear modulus corresponding to the lower bound of G values, and they can be written as GV = (C11 C12 + 3C44)/5, and 5/GR = 4/(C11 C12) + 3/C44. Our results are given in Table 2 along with the theoretical ones [4]. The present average value of A shows that Rh3Hf is an elastically anisotropic material, and it is slightly smaller than that for Rh3Ti and Rh3V given in Ref. [12]. The bulk modulus is a measure of resistance to volume change by applied pressure, whereas the shear modulus is a measure of resistance to reversible deformations upon shear stress [36]. Therefore, isotropic shear modulus is a better predictor of hardness than the bulk modulus. The calculated average isotropic shear modulus is about 123 GPa for this structure. We conclude that Rh3Hf is a sort of low compressible compound. According to a com-
mon criterion [37,38] in the literature, a material is brittle (ductility) if the B/G ratio is less (high) than 1.75. The present average value of the B/G is higher than 1.75, thus this material will behave in a ductility manner. The typical value of Poisson’s ratio is about t = 0.1 for covalent materials and 0.25 for ionic materials [38]. In the present case the value of t is 0.28 for Rh3Hf, and it is consistent with the theoretical value of Chen et al. [4] found for this compound. Thus, the ionic contributions to the atomic bonding are dominant for Rh3Hf. The Young’s modulus is defined as the ratio of the tensile stress to the corresponding tensile strain, and is an important quantity for technological and engineering applications. It provides a measure of the stiffness of a solid, and the material is stiffer for the larger value of Young’s modulus. The computed value of Young’s modulus obtained from the both functional LDA and GGA is given in Table 2, and the result obtained from GGA is in agreement with the theoretical value of Ref. [4]. 3.2. Thermodynamic properties The thermal properties are determined in the temperature range from 0 to 2000 K for Rh3Hf compound using the results obtained from LDA, where the quasi-harmonic model remains fully valid. The pressure effect on the some thermodynamical quantities is studied in the range of 0–180 GPa. The variation of volume with the temperature and/or pressure are shown in Figs. 1a and 1b, respectively. The temperature dependence of volume almost re-
Fig. 1a. The variation of volume with temperature for various pressures.
862
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865
Fig. 1b. The variation of volume with the pressure of Rh3Hf for various temperatures.
main constant at higher pressures, but its slope slightly increases at the lower pressures (see Fig. 1a). It is seen from the Fig. 1b that when the pressure increases from 0 GPa to 180 GPa, the volume almost exponentially decreases. This behavior can be attributed to the strengthen of atomic interaction due to becoming closer of the atoms in the interlayer. The variation of bulk modulus with temperature for the bulk modulus at lower pressures remains almost constant, but relatively at higher pressures the bulk modulus slightly decreases with temperature (see Fig. 2a). It is seen from Fig. 2b that the bulk modulus rapidly increases-almost linearly-with pressure for all temperatures. In the quasi-harmonic Debye model, the Debye temperature hðTÞ and the Grüneisen parameter cðTÞ are two key quantities, and they are very sensitive to the vibrational modes. Their values
at various temperatures (100, 300, 1000, 1500, and 2000 K) and pressures (0, 45, 90, 135, and 180 GPa) are given in Tables 3 and 4, respectively. According to the Tables 3 and 4, when the temperature increases, the Debye temperature decreases and the Grüneisen parameter increases, but when the pressure increases the opposite situation takes place for this compound. The heat capacity at constant volume Cv at different temperatures T and pressures P is shown in Fig. 3 for Rh3Hf. It is realized from the figure that when T < 800 K, the Cv increases very rapidly with the temperature; when T > 800 K, the Cv increases slowly with the temperature, and it almost approaches a constant value called as Dulong–Petit limit for this compound. The variation of the thermal expansion coefficientðaÞ with temperature with various pressures is shown in Fig. 4 for Rh3Hf. It is clearly seen that a exhibits similar trend for all isobars, and it rap-
Fig. 2a. The variation of bulk modulus with temperature for various pressures.
863
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865
Fig. 2b. The variation of bulk modulus with pressure for various temperatures.
Table 3 The calculated Debye temperature H (K) over a wide temperature and pressure range for Rh3Hf compound. T (K)
100 300 1000 1500 2000
Rh3Hf P (GPa)
H H H H H
Table 4 The calculated the Grüneisen parameter over a wide temperature and pressure range for Rh3Hf compound. T (K)
0
45
90
135
180
543.85 539.98 519.41 502.28 483.88
699.68 697.77 686.35 677.52 668.34
800.46 799.25 791.80 785.98 779.80
880.90 879.85 873.71 868.91 863.81
947.69 946.66 941.21 936.92 932.44
100 300 1000 1500 2000
Rh3Hf P (GPa)
c c c c c
0
45
90
135
180
2.129 2.142 2.213 2.275 2.344
1.700 1.704 1.730 1.751 1.773
1.496 1.498 1.512 1.523 1.534
1.363 1.365 1.374 1.382 1.390
1.269 1.270 1.277 1.283 1.289
Fig. 3. The variation of Cv with temperature at various pressures for Rh3Hf.
864
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865
idly increases with the temperature T at lower temperatures (about T < 500 K) and the above this value (about T > 500 K) its increasing rate gradually decreases for lower pressures. The thermal expansion coefficient possesses the highest values for lowest (P = 0 GPa) pressure in all temperature range considered. 3.3. Phonon dispersion curves Although the superiority of GGA on LDA is generally accepted, but not for all physical properties and all types of materials, we have chosen the local density approximation (LDA) instead of general gradient approximation (GGA) because Grobowski et al’s [39] have obtained more accurate results using LDA than GGA for phonon dispersion curves of elemental Rh. Recently, Barrera et al. [40] have also observed the superiority of LDA on GGA for anharmonic properties of alkali metal hydrides compounds at high-temperatures. These results have encouraged us to use LDA instead of GGA for lattice dynamical properties of Rh3Hf alloy. The present LDA phonon frequencies of Rh3Hf compound in L12 phase are calculated using the PHON program [41] utilizing the
Helmann–Fynman forces obtained from the VASP. The PHON code calculates force constant matrices and phonon frequencies using the ‘‘Small Displacement Method” that is described in Refs. [42,43]. The present phonon dispersion curves and one-phonon density of state in high symmetry directions have been calculated using a 2 2 2 cubic supercell of 32 atoms. The obtained results along the high symmetry direction are illustrated in Fig. 5 for Rh3Hf. Unfortunately, for the lattice dynamics of this compound, there is neither experimental nor other theoretical data for comparing with the present data. Fig. 5 shows the phonon dispersion curves and corresponding density of states for Rh3Hf. These curves are very similar to those obtained for other platinum-based alloys [1,12] in the same structure. The calculated phonon dispersion curves do not contain soft mode at any vectors, which confirms the stability of Rh3Hf in L12 structure. The unit cell of Rh3Hf has four atoms, which give rise to a total of 12 phonon branches, which contains three acoustic modes and nine optical modes. It is well known that the mass difference between anions and cations strongly affects the maximum and minimum values of the acoustic and optic branches. The
Fig. 4. The variation of the thermal expansion coefficient with temperature at various pressures.
Fig. 5. The phonon dispersions and corresponding density of states for Rh3Hf.
G. Surucu et al. / Computational Materials Science 48 (2010) 859–865
maximum value of optical branches (117.5 cm1) for Rh3Hf are higher than corresponding values for Rh3V (106.86 cm1) and Rh3Ti (102.5 cm1) in Ref. [12]. Thus, the maximum values of the phonon frequencies for optical branches decrease on going from Hf to V and Ti atom, and a clear gap between the acoustic and optic branches is not observed for Rh3Hf. 4. Summary and conclusion We have performed the first principles total energy calculation for Rh3Hf using the plane-wave pseudopotential approach within the local density approximation and generalized gradient approximation. The calculated lattice parameter and bulk modulus are, reasonably, consistent with the literature values. The present elastic constants satisfy the traditional mechanical stability conditions and their values from both ‘‘volume-conserving” and ‘‘stress– strain” methods are consistent with each other. The analysis of the elastic properties shows that Rh3Hf is a soft and ductile material. Beside the other contributions, the original aspect of the present paper is concerned with the phonon dispersion curves and thermo-elastic results, which have not been considered so far. Consequently, we conclude that our theoretical predictions on the considered properties of Rh3Hf compound would be helpful for the further experimental and theoretical investigations in the future. Acknowledgments This work is supported by Gazi University Research-Project Unit under Project No: 05/2009-55. References [1] E.I. Isaev, A.I. Lechtenstein, Y.K. Vekilov, E.A. Smirnova, I.A. Abrikosov, S.I. Simak, R. Ahuja, B. Johansson, Solid State Commun. 129 (2004) 809–814. [2] M. Sundareswari, M. Rajagopalan, Eur. Phys. J. B 49 (2006) 67. [3] Y. Terada, K. Ohkubo, S. Miura, J.M. Sanchez, T. Mohri, J. Alloys Comp. 354 (2003) 202. [4] K. Chen, L.R. Zhao, J.S. Tse, J.R. Rodgers, Phys. Lett. A 331 (2004) 400. [5] M. Rajagopalan, M. Sundareswari, J. Alloys Comp. 379 (2004) 8.
865
[6] S. Miura, K. Honma, Y. Terada, J.M. Sanchez, T. Mohri, Intermetallics 8 (2000) 785. [7] A.E. Dwight, P.A. Beck, Trans. AIME 215 (1959) 976. [8] D.L. Ritter, B.C. Giessen, N.T. Grant, Trans. Met. Soc. AIME 230 (1964) 1250. [9] P. Villars, L.D. Calvert, Pearson’s Handbook of Crystallographic Data for Intermetallic Phases, ASM International, Metals Park, OH, 1985. [10] S.M.L. Sastry, R. Mahapatra, A.W. Davis, Mater. Sci. Eng. A329-331 (2002) 486. [11] K. Chen, L.R. Zhao, J. Appl. Phys. 93 (2003) 2414. [12] G. Surucu, K. Colakoglu, E. Deligoz, H. Ozisik, Intermetallics 18 (2010) 286. [13] Y. Yamabe, Y. Koizumi, H. Murakami, Y. Ro, T. Maruko, H. Harada, Scr. Metall. Mater. 35 (1996) 211. [14] Y. Yamabe, Y. Koizumi, H. Murakami, Y. Ro, T. Maruko, H. Harada, Scr. Metall. Mater. 36 (1997) 393. [15] W.B. Pearson, A Handbook of Lattice Spacings and Structures of Metals and Alloys, Vol. 2, Pergamon Press, Oxford, 1967, p. 1. [16] G. Kresse, J. Hafner, Phys. Rev. B 47 (1994) 558. [17] G. Kresse, J. FurthmÄuller, Comp. Mater. Sci. 6 (1996) 15. [18] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [19] G. Kresse, J. FurthmÄuller, Phys. Rev. B 54 (1996) 11169. [20] P.E. Blochl, Phys. Rev. B 50 (1994) 17953. [21] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048. [22] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. [23] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188–5192. [25] A.A. Maradudin, E.W. Montroll, G.H. Weiss, I.P. Ipatova, Theory of Lattice Dynamics in the Harmonic Approximation, Academic Press, 1971. [26] M.A. Blanco, PhD Thesis, Universidad de Oviedo, 1997.
. [27] M.A. Blanco, A. Martín Pendás, E. Francisco, J.M. Recio, R. Franco, J. Molec. Struct. Theochem. 368 (1996) 245. [28] E. Francisco, J.M. Recio, M.A. Blanco, A. Martín Pendás, J. Phys. Chem. 102 (1998) 1595. [29] E. Francisco, G. Sanjurjo, M.A. Blanco, Phys. Rev. B 63 (2001) 094107. [30] M. Flórez, J.M. Recio, E. Francisco, M.A. Blanco, A. Martín Pendás, Phys. Rev. B 66 (2002) 144112. [31] R. Hill, Proc. Phys. Soc. London A 65 (1952) 349. [32] F.D. Murnaghan, Proc. Natl., Acad. Sci. USA 30 (1994) 5390. [33] M.J. Mehl, Phys. Rev. B 47 (1993) 2493. [34] O.H. Nielsen, R.M. Martin, Phys. Rev. Lett. 50 (1983) 697. [35] B. Mayer, H. Anton, E. Bott, M. Methfessel, J. Sticht, P.C. Schmidt, Intermetallics 11 (2003) 23. [36] I.R. Shein, A.L. Ivanovskii, J. Phys.: Condens. Matter 20 (2008) 415218. [37] S.F. Pugh, Phil. Mag. 45 (1954) 833. [38] V.V. Bannikov, I.R. Shein, A.L. Ivanovskii, Phys. Stat. Sol. (RRL) 1 (2007) 89. [39] B. Grabowski, T. Hickel, J. Neugebauer, PRB 76 (2007) 024309. [40] G.D. Barrera, D. Colognesi, P.C.H. Mitchell, A.J. Ramirez-Cuesta, Chem. Phys. 317 (2005) 119. [41] . 1998. [42] D. Alfè, G.D. Price, M.J. Gillan, Phys. Rev. B 64 (2001) 045123. [43] G. Kresse, J. Furthmüller, J. Hafner, Europhys. Lett. 32 (1995) 729.