Diamond & Related Materials 18 (2009) 632–636
Contents lists available at ScienceDirect
Diamond & Related Materials j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / d i a m o n d
Analytic equation of state and thermodynamic properties of diamond based on an analytic mean field approach Yang Wei ⁎, Jiu-xun Sun, Rong-gang Tian Department of Applied Physics, University of Electronic Science and Technology, Chengdu 610054, People's Republic of China
a r t i c l e
i n f o
Article history: Received 25 May 2008 Received in revised form 1 September 2008 Accepted 4 November 2008 Available online 27 November 2008 Keywords: Diamond crystal Mechanical properties Thermal properties Superhard materials
a b s t r a c t The analytic mean field potential (AMFP) method is applied to diamond. The analytic expressions for the Helmholtz free energy, internal energy and equation of state (EOS) are derived. The formalism for the case of the Morse potential is used in this work. The six potential parameters are determined by fitting the compression experimental data of diamond at ambient temperature. The thermophysical properties including the isothermals, thermal expansion, isochoric heat capacity, Helmholtz free energy and internal energy are calculated and analyzed. The theoretical results agree well with the experimental data available for diamond. These results presented in this paper verify that the AMFP method is a useful approach to consider the anharmonic effects at high temperature. Numerous reasonable predictions and the change trend of the properties for diamond at extreme conditions have been given. © 2008 Elsevier B.V. All rights reserved.
1. Introduction Diamond is a material with many striking and unique properties, e.g., the largest elastic moduli known for any material and correspondingly, the largest sound velocities [1,2]; a very large Debye temperature, making it a “quantum” crystal even at room temperature [1,3,4]. Its hardness and abrasive qualities, highly valued in technology and gem industry, are controlled by the large elastic moduli. Nowadays, the diamond anvil cell (DAC) technique, which is widely used as a powerful tool for generating high pressure in the laboratory [5], exploits the extreme hardness of diamond. Because of these remarkable properties, diamond has been extensively studied [6–17]. The possibility of producing under pressure high-density polymorphs of diamond, including metallic forms, has been discussed [6–8]. Structural changes have already been reported in diamond under nonhydrostatic pressures around 150 GPa and large deformation [9]. However, measurements [10,11] of the properties of diamond under hydrostatic pressure have been limited to below 40 GPa. Occelli et al. [12] report accurate measurements of the volume and of the optical phonon frequency of diamond under hydrostatic pressure up to 140 GPa. They show that diamond is more compressible than currently expected. Maezond et al.[13] have studied the equation of state and Ranman frequency of diamond using quantum Monte Carlo methods. Their EOS are in good agreement with the experimental data of Ref.[14], and with PBE-GGA results, but disagree with those of Ref. [12]. Zouboulis et al. have measured the elastic moduli of diamond by ⁎ Corresponding author. E-mail address:
[email protected] (Y. Wei). 0925-9635/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.diamond.2008.11.015
using Brillouin Scattering in the temperature range 300 to 1600 K [15]. Mounet et al. [16] provide a converged accurate determination of the structural, dynamical, and thermodynamic properties of diamond from first-principles calculations. Recently, Dewaele et al.[17] have measured the thermal expansion of diamond in the pressure– temperature range of 0–80 GPa and 300–900 K by X-ray diffraction in a resistively heated diamond-anvil cell. Although there exist many reports on the thermodynamic properties of diamond, most of them are studied in the harmonic or quasiharmonic approximation framework and are not in good agreement with each other. In order to consider the anharmonic effects at high temperature and determine the analytic results for diamond, in this paper, we apply an analytic approach to this material. Several years ago, Wang et al. [18–21] proposed the analytic mean field potential (AMFP) approach, and they applied it to many materials. Bhatt et al. [22,23] further applied the AMFP method to lead and alkali metals, and concluded that in comparison with other theoretical models the AMFP method is computationally simple, physically transparent and reliable in the study of thermodynamic properties at high pressures and high temperatures. Recently, Sun et al. have proven that the AMFP method is an analytic approximation of the free volume theory (FVT) [24]. The FVT is a mean field approximation to the thermal contribution of atoms to the Helmholtz free energy of crystalline phases. It is more valuable to directly use the strict FVT than the approximate AMFP, in the cases that the analytic equation of state can be derived based on the strict FVT. Nevertheless, in some cases the mean field integral and the equation of state (EOS) for the strict FVT are fairly complicated or cannot be analytically derived. Then it is convenient to develop simple analytic EOS through
Y. Wei et al. / Diamond & Related Materials 18 (2009) 632–636
the AMFP method, when the complete FVT fails. The AMFP method has been applied to solid C60 by using the Girifalco potential [25] by Sun et al. [26]; the numerical results agree well with the molecular dynamics (MD) simulations [27,28] and are superior to the CUSF of Zubov et al. [29,30]. This verifies that the AMFP method is a convenient approach to consider the anharmonic effects at high temperature. Thus, in this paper, we present the results of thermodynamic properties of diamond by using the AMFP method. The rest of this paper is organized as follows. In Section 2 we derive the analytic EOS based on the AMFP approach. In Section 3, the parameters of the Morse potential [31] are determined by fitting the compression experimental data of diamond, and the numerical results of thermodynamic properties are calculated and compared with the experimental data. In Section 4 the conclusion is presented.
Sun et al [24] have proved that the AMFP method is an analytic approximation of the FVT [32,33]. The fundamental spirit of the AMFP to use all cold energy to evaluate thermal properties is in contradiction with the embedded atom model (EAM) [24]. It is necessary to improve the AMFP in terms of the EAM by the replacement of all cold energy with cold energy from interaction between metallic atoms. According to this conclusion and the detail structure of diamond, we have modified the expression of the free energy in reference [33] as follows ð1Þ
where μ is the mass of per carbon atom, vf is the free volume [33], vf = 4π∫0rm exp½−g ðr; V Þ=kT r 2 dr
ð2Þ
rm is the largest displace, rm is approximate to the Wigner–Seitz radius, rm = (3a3/4πγ)1/3 ≈ a/2. To obtain an analytic approximation to the EOS, we propose the Einstein model. It must be pointed out that, except at very low temperatures, no significant improvement is gained by using the Debye model [34]. Fqm/NkT is the quantum modification [35], by using the Einstein model, we have Fqm Fq Fq = − lim = 3ln 1−e−ΘE =T −3lnðΘE =T Þ NkT NkT T→∞ NkT
ð3Þ
Ec (a) = E1(a) + E2(a), is the cohesive energy of an atom, a is the nearestneighbor distance. Based on the embedded atom model, Ec (a) can be divided into two parts: E1(a) and E2(a). For metal, the first part E1(a) represents the contribution of electron gas. As for diamond, we think that the first part E1(a) is the partial contribution of chemical bond. The second part E2(a) represents the contribution of the van der Waals interaction between atoms. By using the Morse potential [31], we have
E1 ðaÞuE1 ðy1 Þ = e1 fexp½2λ1 ð1−y1 Þ−2exp½λ1 ð1−y1 Þg E2 ðaÞuE2 ðy2 Þ = e2 fexp½2λ2 ð1−y2 Þ−2exp½λ2 ð1−y2 Þg
ð4Þ
y1, y2 is the reduced volume, (
y1 = a=r01 = ðV=V01 Þ1=3 ; V01 = Nðr01 Þ3 =γ : y2 = a=r02 = ðV=V02 Þ1=3 ; V02 = Nðr02 Þ3 =γ
The g (r, V) in Eq. (2) is the potential energy of an atom as it roams from the center atom to a distance r. In terms of the AMFP approach [18–21], g(r, V) can be expressed by E2(a) as follows r i 1 h r 1 + E2 ða + r Þ + 1− E2 ða−rÞ−2E2 ðaÞ ð6Þ g ðr; V Þ = 2 a a In order to consider the quantum effect, we develop g(r, V) into the quadratic function in the equilibrium position.
1 1 @ 2 E2 ðaÞ 2 @E2 ðaÞ r2 + ð7Þ g ðr; V Þ = μω2 r 2 = 2 2 2 a @a @a In terms of Eq. (7), the harmonic vibration frequency ω and Einstein temperature ΘE can be determined as
1=2 ħω ħ @ 2 E2 ðaÞ 2 @E2 ðaÞ = pffiffiffi + k k μ a @a @a2 " #1=2 ħ @ 2 E2 ðaÞ 2 @E2 ðaÞ = + pffiffiffi kr02 μ y2 @y2 @y22
ΘE =
2. Analytic equation of state
1 Fqm F 3 ½E1 ðaÞ + E2 ðaÞ−lnvf + = − ln 2πμkT=h2 + NkT 2 kT NkT
633
ð5Þ
where r01, r02 is the equilibrium distance, ε1, ε2 is the corresponding well depth, λ1, λ2 describe the decrease of potential as the distance increases. γ is structural constant, pffiffiffithe values for fcc structure have been given in Ref. [36], and γ = 2. The volume of the fcc solid is V = Na3/γ. Diamond is not the fcc solid. However, our calculations show pffiffiffi that the results are not sensitive to the values of γ, thus we take γ = 2 for diamond.
ð8Þ
And the Grüneisen parameter can be derived as V @ΘE y2 @ΘE =− ΘE @V 3ΘE @y2 1 @ 3 E2 =@y32 + ð2=y2 Þ@ 2 E2 =@y22 − 2=y22 @E2 =@y2 =− 6 @ 2 E2 =@y22 + ð2=y2 Þ@E2 =@y2
γG = −
ð9Þ
For simplicity, we introduce dimensionless reduced free volume v̄f, and reduced radial coordinate x as follows vf = 4πa3 v f = 4πγ ðV=NÞv f x = r=a;
rm 1 ≈ xm = a 2
ð10Þ ð11Þ
The reduced free volume v̄f and its derivatives with respect to temperature and reduced volume can be expressed as vf = ∫
xm 0 exp½−g ðx; y2 Þ=kT
x2 dx
ð12Þ
v fa = T
@ 1 xm ∫ exp½−g ðx; y2 Þ=kT g ðx; y2 Þx2 dx v = @T f kT 0
ð13Þ
v fb = −
@ 1 xm @ ∫ exp½−g ðx; y2 Þ=kT v = g ðx; y2 Þx2 dx @y2 f kT 0 @y2
ð14Þ
Here g(x, y2) ≡ g(r, V), combining Eqs. (4), (6) and (10), we have g ðx; y2 Þug ðr; V Þ =
1 ½ð1 + xÞE2 ðy2 + y2 xÞ + ð1−xÞE2 ðy2 −y2 xÞ−2E2 ðy2 Þ ð15Þ 2
The compressibility factor can be derived as Z=
Pf V Pqm V PV a @ F Pc V =− = + + NkT 3 @a NkT NkT NkT NkT P V
P
ð16Þ
V
Pc V f qm where NkT , NkT and NkT represent the contribution of cold energy, the contribution of the free volume and the quantum modification, respectively.
Pc V a @ y1 ′ y2 ′ =− E1 ðy1 Þ− E2 ðy2 Þ ½E1 ðy1 Þ + E2 ðy2 Þ = − NkT 3kT @a 3kT 3kT
ð17Þ
Pf V y2 v fb y2 @ =1+ v = 1− NkT 3v f @y2 f 3v f
ð18Þ
Pqm V γ G Uqm = NkT NkT
ð19Þ
where Pc is the cold pressure, Pf is the thermal pressure deduced from free volume.
634
Y. Wei et al. / Diamond & Related Materials 18 (2009) 632–636
The internal energy can be derived as U @ F 3 1 T @v f = −T = + ½E1 ðy1 Þ + E2 ðy2 Þ + NkT @T NkT 2 kT v f @T v fa Uqm 3 1 ½E1 ðy1 Þ + E2 ðy2 Þ + + = + 2 kT vf NkT
ð20Þ
where the first term represent the ideal gas, the second term, the third term and the fourth term represent the contribution of cold energy, the contribution of the free volume and the quantum modification, respectively. Uqm @ Fqm 3ΘE =T = −T = −3 @T NkT eΘE =T −1 NkT
ð21Þ
By using above equations, all other thermodynamic quantities can be analytically derived. The derivations are straightforward. However, the expressions for thermal expansion coefficient, compressibility coefficient and isochoric heat capacity are redundant, we would calculate these quantities by using numerical differentiation instead of the analytic expressions. The compressibility factor can be seen as function of variables T and V, Z = Z(T, V). In terms of the function, the formulas for thermal expansion coefficient, compressibility coefficient and isochoric heat capacity can be reduced to following form
−1 1 @V Z @Z @Z α= + = ð22Þ Z−V V @T P T @T V @V T
−1 1 @V V @Z Z−V = V @P T NkT @V T CV 1 @U U @ U +T = = @T NkT V Nk Nk @T V NkT
β=−
ð24Þ ð23Þ
In our calculations, it is found that following steps for the numerical differentiations in Eqs. (22–24) can reach stable numerical results, ΔT = 0.00001 × T and ΔV = 0.00001 × V. 3. Results and discussion In this section we apply the above formalism to diamond. Occelli et al. [12] have measured the isothermals of diamond under hydrostatic
Fig. 2. Potential energy E as function of the nearest-neighbor distance a calculated for diamond. E1(a) (dashed line) represents the partial contribution of chemical bond. E2(a) (dash–dotted line) represents the contribution of the van der Waals interaction between atoms. Ec(a) (solid line) represents the cohesive energy of an atom.
pressure up to 140 GPa at ambient temperature. Considering the structure of diamond is tetrahedron, we thus determined the parameters for the Morse potential in Eq. (4) by fitting the compression experimental data [12]. The experimental data and smoothed fitting curves are plotted in Fig. 1. The figure shows that the fitting precision is satisfactory and our results are in good agreement with the available experimental data. This suggests that the AMFP method is a useful approach to study the thermodynamic properties of diamond. Furthermore, we predict the variation of the pressure dependence of the molar volume at several higher temperatures (T = 1200 K, 2000 K and 3000 K), and also plot the results in Fig. 1. The results suggest that the molar volume of diamond decreases with pressure and increases with temperature. The determined values of the parameters are as follows:
Fig. 1. Isothermal curves at 296 K,1200 K, 2000 K, 3000 K calculated by using the parameters in Eq. (25). The open circles (o) are experimental data taken from Occelli et al. [12].
e1 = 188; 000K; λ1 = 1:66; V01 = 3:4206 cm3=mol; e2 = 159; 000K; λ2 = 2:622; V02 = 3:4133 cm3=mol;
ð25Þ
Substituting the determined parameter values (25) into Eq. (4), we obtain the specific form of E1(a), E2(a) and Ec(a) for diamond. We plot E1(a), E2(a) and Ec(a) versus the nearest-neighbor distance a of diamond in Fig. 2. The figure shows that the difference of tendency between E1(a) and E2(a) is fairly small. This means that the partial contribution of chemical bond and the contribution of the van der Waals interaction between atoms are identical. But the well depth of Ec (a) is fairly large. This is the reason that the strong, tetrahedrally coordinated, covalent bonds between nearest neighbors and the light mass of the constituent atoms lead to many striking and unique properties of diamond. In Figs. 3–6, the results of thermodynamic properties of diamond calculated by using the Morse potential are plotted. The variation of the bulk modulus versus pressure curves at 296 K, 1200 K, 2000 K, 3000 K calculated by using the parameters of (25) are plotted in Fig. 3. From this figure, we can see that the difference of the variation tendency among the four curves is not prominent at all. The bulk modulus of diamond increases with the pressure increasing when the temperature is
Y. Wei et al. / Diamond & Related Materials 18 (2009) 632–636
635
Fig. 3. Variations of bulk modulus BT versus pressure P at different temperatures for diamond. The curves represent results calculated by using the AMFP method. Fig. 5. The same as for Fig. 4, but for internal energy U.
changeless and the bulk modulus of diamond decreases with the temperature increasing when the pressure is changeless. The calculated free energy and internal energy as functions of the molar volume for diamond at four temperatures (296 K, 1200 K, 2000 K, 3000 K) are plotted in Figs. 4 and 5, respectively. The two figures show that both the free energy and internal energy of diamond are decreasing function of the molar volume and the two physical quantities have the same variation tendency as the molar volume increase at four different temperatures. On the other hand, the two figures also show that the free energy and internal energy elevate with the increasing temperature, and the values of both physical quantities are negative and quite little at ambient temperature. This suggests that the interatomic interaction of diamond is very strong. At the same time, the interaction decreases with the volume increasing at high pressure region and approaches to a flat at low pressure region under the condition of the fixed temperature. The temperature dependence of the molar volume, isochoric heat capacity, bulk modulus and thermal expansion coefficient for diamond at zero pressure from theoretical calculation in this paper are shown in Fig. 6. From the Fig. 6 (a), it can be found that the molar volume of diamond increases with temperature, at low temperature the variation is slow, but at high temperature the variation is fast. Fig. 6 (b) shows that isochoric heat capacity is an increasing function of temperature, the variation is fast at low temperature and slow at high temperature. Fig. 6 (c) presents the isothermal bulk moduli calculated by using the present AMFP method, ab initio calculations [16] and Raman-based Debye model calculations [17]. From this figure, it can be found that the bulk modulus of diamond is large at ambient temperature and
decreases slowly with the increasing temperature at zero pressure. The Raman-based Debye model [17] overestimates the bulk modulus at all temperatures, however, the results of ab initio calculations [16] are underestimated the bulk modulus up to 3000 K. At the same time, our calculation curve is lies between the curve of Raman-based Debye model and the curve of ab initio calculations. This means that our results based on the AMFP method are approach to the intrinsic behaviors of diamond. Fig. 6 (d) compares our predictions of linear thermal expansion coefficient versus temperature with the available experimental data [37] and the other authors' results [16,38,39] calculated by using different methods. From this figure, it is worth point out that our calculated results are in good agreement with the available experimental data [37] and the ab initio calculations [16]. This further verifies that the AMFP method is a very useful approach to study the thermodynamic properties of diamond. At the same time, our analytic calculation based on the AMFP method is obviously superior to the QHA-GGA thermal expansion calculated using the Grüneisen equation [16] and QHA-LDA [39]. The reason is that the AMFP method is considered the anharmonic effects at high temperature, but the QHA-GGA and QHALDA are in the quasi-harmonic approximation. In addition, the agreement between the results calculated by using different methods and the available experimental data is in fairly good as expected at low temperature, at high temperature our model tends to exhibit deviations. This phenomenon is often observed for most current EOS's, the deviations observed at high temperature confirm the fact that the thermal pressure is not a strictly linear function of temperature above the Debye temperature. 4. Conclusion
Fig. 4. Variation of free energy F versus molar volume V at different temperatures for diamond. Here the contribution of ideal gas to free energy is neglected, the curves represent results calculated by using the AMFP method.
In summary, the AMFP approach has been applied to diamond. The analytic expressions for the Helmholtz free energy, internal energy and EOS have been derived. By fitting the compression experimental data of diamond, we determined the parameters for the Morse potential and calculated the pressure dependence of the molar volume for diamond based on the AMFP method. The results obtained are in good agreement with the available experimental data. This suggests that the AMFP method is a very useful approach to study the thermodynamic properties of diamond. Additionally, we predicted the variation tendency of compression curves at high temperatures and present the potential function of diamond from theoretical calculation. Furthermore, the variation relationship of the free energy and internal energy versus the molar volume of diamond at various temperatures have been predicted and the temperature dependence of the molar volume, bulk modulus, thermal expansion coefficient and
636
Y. Wei et al. / Diamond & Related Materials 18 (2009) 632–636
Fig. 6. The molar volume (a), isochoric heat capacity (b), bulk modulus (c) and thermal expansion coefficient (d) of diamond as a function of temperature at zero pressure, respectively. (c). solid line, AMFP; dashed line, ab initio calculations [16]; dashed–dotted line, Raman-based Debye model calculations [17]. (d). solid line, AMFP; open circles, experimental data [37]; dashed line, ab initio calculations [16]; the signals “+”, a path integral Monte Carlo study using a Tersoff empirical potential [38]; dashed–dotted lined, QHALDA study by Pavone et al.[39]; dotted line, The QHA-GGA thermal expansion calculated using the Grüneisen equation [16].
isochoric heat capacity at zero pressure have been calculated by using the formalism developed in this paper. Acknowledgement This work was supported by the Support Programs for Academic Excellence of Sichuan Province of China under Grant No. 06ZQ026010, which of the Education Ministry of China under Grant No. NCET05-0799, and which of UESTC under Grant No. 23601008. References [1] M.H. Grimsditch, A.K. Ramdas, Phys. Rev., B 11 (1975) 3139. [2] R. Vogelgesang, A.K. Ramdas, S. Rodriguez, M. Grimsditch, T.R. Anthony, Phys. Rev., B 54 (1996) 3989. [3] D.L. Burk, S.A. Friedberg, Phys. Rev. 111 (1958) 1275. [4] J.E. Desnoyers, J.A. Morrison, Philos. Mag. 3 (1958) 42. [5] H.K. Mao, R.J. Hemley, Nature (London) 351 (1991) 721. [6] S. Fahy, S.G. Louie, Phys. Rev., B 36 (1987) 3373. [7] S.J. Clark, C.J. Ackland, J. Crain, Phys. Rev., B 52 (1995) 15035. [8] S. Serra, G. Benedek, M. Facchinetti, L. Miglio, Phys. Rev., B 57 (1998) 5661. [9] H.K. Mao, R.J. Hemley, Nature 351 (1991) 721. [10] M. Hanfland, K. Syassen, S. Fahy, S.G. Louie, M.L. Cohen, Phys. Rev., B 31 (1985) 6896. [11] P.h. Gillet, G. Fiquet, I. Daniel, B. Reynard, M. Hanfland, Phys. Rev., B 60 (1999) 14660. [12] F. Occelli, P. Loubeyre, R. LeToullec, Nat. Matters 2 (2003) 151. [13] R. Maezono, A. Ma, M.D. Towler, R.J. Needs, Phys. Rev. Lett. 98 (2007) 025701.
[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
H.J. McSkimin, P. Andreatch Jr., J. Appl. Phys. 43 (1972) 2944. E.S. Zouboulis, M. Grimsditch, A.K. Ramdas, S. Rodriguez, Phys. Rev., B 57 (1998) 2889. N. Mounet, N. Marzari, Phys. Rev., B 71 (2005) 205214. Agnès Dewaele, Frédéric Datchi, Paul Loubeyre, Mohamed Mezouar, Phys. Rev., B 77 (2008) 094106. Y. Wang, D. Chen, X. Zhang, Phys. Rev. Lett. 84 (2000) 3220. Y. Wang, Phys. Rev., B 62 (2000) 196. Y. Wang, Phys. Rev., B 63 (2001) 245108. Y. Wang, R. Ahuja, B. Johansson, Phys. Rev., B 65 (2001) 014104. N.K. Bhatt, A.R. Jani, P.R. Vyas, V.B. Gohel, Physica B 65 (2001) 014104. N.K. Bhatt, P.R. Vyas, A.R. Jani, V.B. Gohel, J. Phys. Chem. Solids 66 (2005) 797. J.X. Sun, L.C. Cai, Q. Wu, F.Q. Jing, Phys. Rev., B 71 (2005) 024107. L.A. Girifalco, J. Phys. Chem. 96 (1992) 858. J.X. Sun, Physica B 34 (2006) 381. M.C. Abramo, C. Caccamo, J. Phys. Chem. Solids 57 (1996) 1751. M.C. Abramo, C. Caccamo, D. Costa, et al., Phys. Rev., E 69 (2004) 031112. V.I. Zubov, N.P. Tretiakov, J.F. Sanchez, A.A. Caparica, Phys. Rev., B 53 (1996) 12080. V.I. Zubov, J.F. Sanchez-Ortiz, N.P. Tretiakov, I.V. Zubov, Phys. Rev., B 55 (1997) 6747. A.I. Karasevskii, W.B. Holzapfel, Phys. Rev., B 67 (2003) 224301. E. Wasserman, L. Stixrude, Phys. Rev., B 53 (1996) 8296. J.X. Sun, L.C. Cai, Q. Wu, F.Q. Jing, Phys. Rev., B 73 (2006) 155431. M. Taravillo, V.G. Baonza, J.M.F. Rubio, J. Núñez, M. Cáceres, J. Phys. Chem. Solids 63 (2002) 1705. J.X. Sun, Q. Wu, L.C. Cai, F.Q. Jing, Chin. Phys. 12 (2003) 6. N.X. Chen, Z.D. Chen, Y.C. Wei, Phys. Rev., B 55 (1997) R5. G.A. Slack, S.F. Bartram, J. Appl. Phys. 46 (1975) 89. C.P. Herrero, R. Ramírez, Phys. Rev., B 63 (2000) 024103. P. Pavone, K. Karch, O. Schütt, W. Windl, D. Strauch, P. Giannozzi, S. Baroni, Phys. Rev., B 48 (1993) 3156.