- Email: [email protected]
Analytic equation of state for FCC C60 solid based on analytic mean-field potential approach
ARTICLE IN PRESS
Physica B 381 (2006) 34–39 www.elsevier.com/locate/physb
Analytic equation of state for FCC C60 solid based on analytic mean-field potential approach Sun Jiu-xun Department of Applied Physics, University of Electronic, Science and Technology, Chengdu 610054, PR China Received 21 October 2005; received in revised form 30 November 2005; accepted 1 December 2005
Abstract The analytic mean-field approach (AMFP) was applied to the FCC C60 solid. For the intermolecular forces the Girifalco potential has been utilized. The analytic expressions for the Helmholtz free energy, internal energy and equation of state have been derived. The numerical results of thermodynamic quantities are compared with the molecular dynamic (MD) simulations and the unsymmetrized selfconsistent field approach (CUSF) in the literature. It is shown that our AMFP results are in good agreement with the MD data both at low and high temperatures. The results of CUSF are in accordance with the AMFP at low temperature, but at high temperature the difference becomes prominent. Especially the AMFP predicted that the FCC C60 solid is stable upto 2202 K, the spinodal temperature, in good agreement with 2320 K from the MD simulation. However, the CUST just gives 1916 K, a temperature evidently lower than the MD data. The AMFP qualifies as a useful approach that can reasonably consider the anharmonic effects at high temperature. r 2006 Elsevier B.V. All rights reserved. PACS: 64.10.+h; 61.48.+c; 50.70.Ce; 62.50.+p Keywords: Equation-of-state; Fullerene; Thermodynamic properties; Anharmonicity
1. Introduction Since their discovery [1] and especially after their preparation in solid state [2,3], fullerenes have been the focus of attention of many scientific groups [3–31]. In the research of thermodynamic properties of fullerites, the Girifalco potential [4] has been extensively applied [5–31]. The C60 molecule consists of 60 carbon atoms symmetrically arranged on a sphere. Girifalco [4] pointed out that the molecules rotate rather freely and therefore can be approximated by spheres for many purposes. This is true even in the solid, at least for temperatures above 300 K [4]. Under these conditions, the hollow molecular cages can be assimilated to spheres whose surface consists of a uniform distribution of carbon sites. The overall interaction between two fullerene particles is then obtained by an integral of the interaction between pairs of sites on different cages. Using the standard (12-6) Lennard–Jones Tel.: +86 28 83202586.
E-mail address: [email protected]. 0921-4526/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2005.12.268
(LJ) function for the atom–atom interaction between two molecules and averaging, Girifalco obtained following intermolecular potential for the C60 fullerene [4]: 1 1 2 ðsÞ ðr=dÞ ¼ 2 þ sðs 1Þ9 sðs þ 1Þ9 s10 1 1 2 1 þ , ð1Þ sðs 1Þ3 sðs þ 1Þ3 s4 1 ¼ N 20 A=12d 6 ;
2 ¼ N 20 B=90d 12 ;
A ¼ 32:0 1060 erg cm6 ; d ¼ 0:71 nm,
d ¼ 0:71 nm,
(2)
B ¼ 55:77 10105 erg cm6 , ð3Þ
where s ¼ r=d, and N 0 ¼ 60. After the Girifalco potential was proposed, it has been extensively used in many works involving physical–chemical properties of C60 system [5–31]. The properties of fluid phase have been extensively researched by using computer simulations [5–31] and model equations of state [14–20]. For example, the molecular dynamics (MD) simulations on the properties
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
of fluid phase have been reported by Ohno et al. [5], Cheng et al. [6] and Rey et al. [7]. The Monte Carlo (MC) simulations have been reported by Hagen et al. [8,9], Caccamo et al. [10,11] and Hasegwa and Ohno [12,13]. At the same time, the density functional theory (DFT) [14–16], modified hypernetted-chain theory [17], hierarchical reference theory [18], Percus–Yevick approximation [19] of the Ornstein–Zernike equation and a variational series mean spherical approximation [20] have been used to research the thermodynamic properties of C60 fluid by the aid of double Yukawa fitting of Girifalco potential. However, for the solid state, the vast majority of experimental and theoretical work on pure and doped solid fullerenes (fullerites and fullerides) deals with their magnetic properties and superconductivity [21–25]. Up to now, just several works have been reported for thermodynamic properties [26–31]. Cheng et al. [26] predicted the phase diagram by using an integral equation approach combined with MD simulations based on the Girifalco potential. Abramo group published MD simulation results [27,28]. As for model equations of state, we just found two works. Girifalco firstly researched thermodynamic properties of FCC C60 under harmonic Debye–Gru¨neisen (DG) theory [29]. Zubov et al. [30,31] pointed out that the anharmonicity of the lattice vibrations of the C60 fullerite is strong at T4700 K. As this temperature is reached, the main anharmonic terms are no longer small corrections and must be included. Zubov et al. [30,31] considered anharmonic modification up to fourth order by using their unsymmetrized self-consistent field (CUSF) approach. However, Shukla et al. [32,33] have shown that the anharmonic expansion is divergent at high temperature. Although Zubov et al. [30,31] improved the convergence of the CUSF approach by using a self-consistent procedure and including anharmonic modification up to fourth order, this may not be enough at higher temperature. For example, the MD simulations [26] by using the Girifalco potential show that FCC C60 solid is stable up to about 2300 K. However, the CUSF predicted the stable temperature just is 1916 K. Thus, it is necessary to use other approaches that can include anharmonic effect soundly to study the system. The free volume theory (FVT) is a mean-field approximation to the thermal contribution of atoms to the Helmholtz free energy of crystalline phases [34–40]. The model assumes that each atom is confined to the Weigner–Seitz cell formed by its nearest neighbors. Although interatomic correlations are neglected, the cell model soundly includes anharmonic terms which are important for high temperatures, but which are ignored in DG model and quasiharmonic lattice dynamics. Many authors have demonstrated that the cell model matches successfully the thermodynamic properties of the FCC LJ crystal [34–36], sodium chloride [37] and exponential-6 model solid [38] calculated from Monte Carlo simulations. Wasserman and Stixrude further applied the cell model to a metallic solid iron [39]. The calculated properties are in
35
good agreement with available static and shock-wave experimental measurements. Sun et al. [40] extended the FVT to a generalized LJ solid including the quantum modification. The numerical results are in good agreement with experimental data of solid xenon, both at low and high temperature. Recently, Wang et al. [41–44] proposed an analytic mean-field potential (AMFP) approach, and applied it to many materials. Bhatt et al. [45,46] further applied the AMFP to lead and alkali metals, and concluded that in comparison with other theoretical models the AMFP is computationally simple, physically transparent and reliable to study the thermodynamic properties in the high pressures and high temperatures environment. However, the approach has been proven equivalent to the FVT, and just is an analytic approximation of FVT [47]. 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. However, in some cases the mean-field integral and the equation of state for the strict FVT are fairly complicated or cannot be analytically derived, the AMFP approach becomes useful. The Girifalco model for the FCC C60 solid just is the case, so we would develop analytic equation of state for the Girifalco model solid by using the AMFP approach. As was pointed out by Girifalco [29], in C60 molecule, the atoms are held together by strong chemical bonds while the interactions among these molecules in the FCC solid arise from van der Waals-type forces are much weaker. To a good approximation, therefore, the thermal properties of FCC C60 can be treated as the sum of those resulting from intermolecular vibrations and those arising from intramolecular vibrations [29]. Here we just consider the intermoleculear contributions, as having been done by Girifalco [29]. In Section 2 we derive analytic equation of state based on the AMFP approach. In Section 3 the numerical results are compared with results of Zubov et al. and MD data of Abramo et al. In Section 4 the conclusion is presented.
2. Analytic equation of state In terms of the FVT, the free energy can be expressed as [34–40] F 3 uð0Þ ¼ lnð2pmkT=h2 Þ þ ln vf , (4) NkT 2 2kT where m is mass of C60 molecule, uð0Þ is the potential energy of a molecule as the lattice is static, vf is the free volume. X Ri X a X uð0Þ ¼ zi z i di zi ðdi yÞ, (5) ¼ ¼ d d ia0 ia0 ia0 Z
rm
vf ¼ 4p 0
exp½gðr; V Þ=kTr2 dr,
(6)
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
36
where Ri ¼ di a is the distance of molecules in ith shell from the central molecule ði ¼ 0Þ. zi and di are structural constants, the values for FCC structure have been given in Ref. [48]. gðr; V Þ is the potential energy of a molecule as it roams from the center to a distance r. In terms of the AMFP approach [41–47], gðr; V Þ can be expressed by the static energy E c ðaÞ of a molecule: gðr; V Þ ¼ 12½ð1 þ r=aÞE c ða þ rÞ þ ð1 r=aÞE c ða rÞ 2E c ðaÞ,
(8)
where rm is the largest displace. For molecule with hardcore of diameter d, rm can be expressed as d a d 1 rm ¼ rw , (9) 2 2 a 1
where rw ð3a3 =4pgÞ3 a=2, is the Wigner–Seitz radius. a the nearest-neighbor distance, and pffiffiffi g is structure constant, for FCC structure it equals to 2. The volume of the FCC solid V ¼ Na3 =g. For simplicity, we introduce dimensionless reduced free volume v¯ f , reduced volume y, and reduced radial coordinate x: 3
vf ¼ 4pa v¯ f ¼ 4pgV v¯ f ,
x ¼ r=a;
(10) vd ¼ Nd 3 =g.
(11)
rm 1 1 1 . xm ¼ 2 y a
(12)
The reduced free volume v¯ f and its derivatives with respect to temperature and reduced volume can be expressed as Z xm v¯ f ¼ exp½gðx; yÞ=kTx2 dx, (13) 0
q 1 v¯ f ¼ v¯ fa ¼ T qT kT
Z
xm
exp½gðx; yÞ=kTgðx; yÞx2 dx,
0
(14) v¯ fb ¼
q 1 v¯ f ¼ qy kT
Z
xm
exp½gðx; yÞ=kT 0
q gðx; yÞx2 dx, qy (15)
1 xm 2 v¯ fc ¼ exp½gðxm ; yÞ=kT 0. 2 y
(16)
1X gðx; yÞ gðr; V Þ zi ½ð1 þ xÞðdi y þ di yxÞ 4 ia0 ð17Þ
ð18Þ
1 9 1 þ þ 9 10 2 sðs 1Þ s ðs þ 1Þ9 1Þ 9 20 1 3 þ 11 þ 1 þ 10 3 2 s sðs þ 1Þ s ðs 1Þ sðs 1Þ4 1 3 8 þ þ . s2 ðs þ 1Þ3 sðs þ 1Þ4 s5
0 ðsÞ ¼ 2
s2 ðs
ð19Þ
The compressibility factor and internal energy can now be derived: PV y q F Pc V Pf V ¼ ¼ þ , (20) Z¼ NkT 3 qy NkT NkT NkT Pc V y q ¼ uð0Þ, NkT 6kT qy
(21)
X q uð0Þ ¼ zi di 0 ðdi yÞ, qy ia0
(22)
Pf V y q y y¯vfb ¼1þ v¯ f ¼ 1 ð¯vfb v¯ fc Þ 1 , NkT 3¯vf qy 3¯vf 3¯vf (23) U q F 3 uð0Þ T q¯vf 3 uð0Þ v¯ fa ¼ T ¼ þ þ þ ¼ þ , NkT qT NkT 2 2kT v¯ f qT 2 2kT v¯ f (24) where Pc is the cold pressure, Pf is the thermal pressure deduced from free volume. By using above equations, all other thermodynamic quantities can be analytically derived. The derivations are straightforward. However, the expressions for thermal expansivity, 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 y and T, Z ¼ Zðy; TÞ. In terms of the function, formulas for thermal expansivity, compressibility coefficient and isochoric heat capacity can be reduced to following form: 1 qV 3 qy a¼ ¼ V qT P y qT P " # yqZ 1 Z qZ þ Z ¼ , ð25Þ T qT y 3 qy T qV 3 qy ¼ qP T y qP T 1 3 Vdy y qZ ¼ , Z 3 qy T NkT CV 1 qU U q U þT ¼ ¼ , qT NkT y Nk Nk qT V NkT b¼
Here gðx; yÞ gðr; V Þ, combining Eqs. (5), (8) and (10), we have
þ ð1 xÞðdi y di yxÞ 2ðdi yÞ,
þ ð1 xÞ0 ðdi y di yxÞ 20 ðdi yÞ,
ð7Þ
1 1X zi ðdi yÞ, E c ðaÞ ¼ uð0Þ ¼ 2 2 ia0
y ¼ a=d ¼ ðV =V d Þ1=3 ;
q 1X gðx; yÞ zi di ½ð1 þ xÞ0 ðdi y þ di yxÞ qy 4 ia0
1 V
ð26Þ
(27)
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
37
for many problems far from the critical point. Since the thermodynamic property of solids is a problem far from the critical point, it is not surprising that the AMFP can give good results [36–46]. However, the CUST is an approach to consider the anharmonic effects by using the method of series expansion [30,31]. Shukla et al. [32,33] have systematically studied the convergence of the anharmonic expansion, it has been shown that the numerical results would oscillate for odd and even order expansions, the convergence is acceptable at low temperature, slows down at slightly high temperature, and diverges at high temperature. In order to improve the convergence of the anharmonic expansion, Zubov et al. [30,31] developed their CUST approach by introducing consistent procedure and retaining the anharmonic terms up to fourth orders. Although the CUST indeed improves the convergence of anharmonic expansion, it may not be enough at high temperature as suggested by the numerical results in present work. As for the correlation effect, it has been shown contributing little to the thermodynamic properties for realistic potential models [36,37], and can be neglected in most cases except for the critical phenomena. In fact, both AMFP and CUST have neglected the correlation effect, so the correlation effect is not responsible to the discrepancies of the two approaches at high temperatures. In Fig. 1, the results for the internal energy and the pressure as functions of the temperature at various lattice spacings are compared with the MD results [28]. The isothermal curves at 300 K and 1900 K are plotted in Fig. 2. The two figures show that our results are in fairly good agreement with MD simulations [28]. Fig. 3 gives the variation of isochoric heat capacity C V , thermal expansivity a and bulk modulus BT versus pressure. The figure shows that C V is an increasing function of pressure, and the variation is slow. The a is a decreasing function of pressure, at low pressure the variation is fast, but at high pressure, the variation slows down and shows some saturation effect. The BT is a linear increasing function of pressure, but the slope
3. Results and discussion The thermodynamic properties calculated by using formalism developed above are listed and compared with the CUSF approach of Zubov et al. [30] in Table 1. The spinodal point T s is the temperature satisfying the condition, BT ðT s Þ ¼ 0. The system is unstable for temperatures above T s . The MD simulations performed by Cheng group [26] and Abramo group give T s ¼ 2320 K [27]. The CUSF approach gives T s ¼ 1916 K [30], a value lower than the MD result about 20%. The reason maybe that the CUSF approach cannot consider anharmonic effects soundly. Although it can consider anharmonic terms up to fourth order, it is not enough at high temperature. So Zubov et al. [30] just computed thermodynamic properties up to 1900 K. However, our calculation gives T s ¼ 2202 K, is in good agreement with MD simulations. So our calculations have been done up to 2200 K. The results in this paper verify the viewpoint that the AMFP is a useful approach to consider anharmonic effects at high temperatures [32–48]. Table 1 shows that our results are in agreement with the CUST of Zubov et al. at low temperature, at high temperature the difference becomes prominent. Especially near the sinodal temperature of CUSF 1900 K, the thermal expansivity coefficient becomes divergent, and the compressibility coefficient tends zero. However, both coefficients take finite values at this temperature, as the spinodal point of AMFP is 2202 K. Although some difference exists for the linear thermal expansivity at ambient temperature, our calculation gives 1:148 105 K1 at 300 K, is in slightly better agreement with the MD data [28] 1:137 105 K1 than the CUSF result 1:084 105 K1 [30]. Now we explain in more detail the different features of CUSF and AMFP responsible for the discrepancies at high temperatures. The AMFP is a mean field approximation to the thermal contribution of atoms to the Helmholtz free energy of crystalline phases [36–46]. It is well known that the mean field approximation can give fairly good results
Table 1 Properties of the FCC phase of C60 T
261.4
400
600
800
1000
1200
1400
1600
1800
a
1.0068 1.0067
1.0084 1.0082
1.0108 1.0106
1.0136 1.0132
1.0166 1.0161
1.0201 1.0194
1.0241 1.0233
1.0290 1.0282
1.0354 1.0354
a
1.130 1.070
1.198 1.126
1.313 1.223
1.456 1.347
1.640 1.512
1.884 1.745
2.227 2.103
2.750 2.763
3.667 4.710
1900 1.0394 1.0426 4.469 12.713
2000
2100
1.0445
1.0517
5.853
9.073
BT
136.04 140.51
124.19 128.89
107.87 112.88
92.40 97.54
77.71 82.73
63.77 68.22
50.53 53.71
37.95 38.61
25.95 21.11
20.11 7.40
14.33
8.41
CV
23.50
23.33
22.41
22.90
21.87
21.30
20.69
20.02
19.26
18.83
18.32
17.69
2200 1.0692 71.86 0.857 16.49
The nearest-neighbor distance in nm, linear thermal expansion coefficient in 105 K1 , the heat capacity in kJ mol1 K1 , the bulk modulus in kbar. The values in the first lines have been calculated in this work, in the second lines are from Zubov et al. [30].
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
38 100
27 CV (kJ/mol K)
P (kbar)
80 60 40 20 0 200
400
600
800
1000 1200 T (K)
1400
1600
1800
26 25 24
2000
0
100
200 P (kbar)
300
400
0
100
200 P (kbar)
300
400
0
100
200 P (kbar)
300
400
-50
0.9 -100 -150 -200 200
400
600
800
1000 1200 T (K)
1400
1600
1800
2000
α (10-5K-1)
U (kJ/mol)
0
Fig. 1. Pressure (top) and internal energy (bottom) as functions of temperature for several lattice spacing a. The symbols refer to MD simulation data for a ¼ 1:36 nm (circles), 1.38 nm (triangles), and 1.42 nm (squares) [28]. The lines are corresponding numerical results calculated in this work.
P (kbar)
300
6000 4000 2000
T = 300 K 200
0 100 0 1.45
1.5
1.55
1.6
1.65
1.7
1.75
ρ (nm-3) 400 300 P (kbar)
0.3 0
BT (kbar)
400
0.6
T = 1900 K 200 100 0 1.45
1.5
1.55
1.6
1.65
1.7
1.75
ρ (nm-3)
Fig. 2. Variation of pressure versus density at two temperatures, 300 K (top) and 1900 K (bottom). The circles refer to MD simulation data [28]. The lines are corresponding numerical results calculated in this work.
Fig. 3. Variation of thermodynamic quantities versus pressure calculated in this work. The heat capacity is in kJ mol1 K1 , linear thermal expansion coefficient in 105 K1 , the bulk modulus is in kbar.
monic effects available at present, can be applied to high temperature, our comparison shows that some discrepancy exists at high temperature between the CUST with MD data and our FVT results. However, the numerical results from the AMFP are in good agreement with the MD simulations available. This further verifies the viewpoint that the AMFP is a useful approach to consider the anharmonic effects at high temperature. Acknowledgment
is fairly large and the increase is fast. This means the Girifalco model solid for the FCC C60 is difficult to be compressed. 4. Conclusion In summary, the AMFP has been applied to the Girifalco model solid of C60. All thermodynamic quantities can be analytically evaluated based on the formalism developed in this paper. Although Zubov et al. claimed their CUST approach, another theory to consider anhar-
This work was supported by the Natural Science Foundation of China under Grant no. 19904002, by the program of middle-aged and young academic cadres of UESTC. References [1] H.W. Kroto, J.R. Heath, S.C. O’Brien, et al., Nature 318 (1985) 162. [2] W. Krat¨schmer, L.D. Lamb, K. Fostinoupolos, R.D. Huffman, Nature 347 (1990) 354. [3] O. Zhou, D.E. Cox, J. Phys. Chem. Solids 53 (1992) 1373.
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
L.A. Girifalco, J. Phys. Chem. 96 (1992) 858. K. Ohno, Y. Maroyama, Y. Kawazoe, Phys. Rev. B 52 (1996) 4078. A. Cheng, M.L. Klein, C. Caccamo, Phys. Rev. Lett. 71 (1993) 1200. C. Rey, L.J. Gallego, J.A. Alonso, Phys. Rev. B 49 (1994) 8491. M.H.J. Hagen, E.J. Meijer, G.C.A.M. Mooij, et al., Nature 365 (1993) 425. M.H.J. Hagen, D. Frenkel, J. Chem. Phys. 101 (1994) 4093. C. Caccamo, D. Costa, A. Fucile, J. Chem. Phys. 106 (1997) 255. D. Costa, G. Pellicane, M.C. Abramo, C. Caccamo, J. Chem. Phys. 118 (2003) 304. M. Hasegawa, K. Ohno, J. Chem. Phys. 111 (1999) 5955. M. Hasegawa, K. Ohno, J. Chem. Phys. 113 (2000) 4315. M. Hasegawa, K. Ohno, J. Phys.: Condens. Matter 9 (1997) 3361. L. Mederos, G. Nawascue´s, Phys. Rev. B 50 (1994) 1301. M. Hasegawa, K. Ohno, Phys. Rev. E 54 (1996) 3928. C. Caccamo, Phys. Rev. E 51 (1995) 3387. M. Tau, A. Parola, D. Pini, L. Reatto, Phys. Rev. E 52 (1995) 2644. H. Guerin, J. Phys.: Condens. Matter 10 (1998) L527. M. Khedr, M.S. Al-Busaidy, S.M. Osman, J. Phys.: Condens. Matter 17 (2005) 4411. J.F. Fischer, P.A. Heiney, J. Phys. Chem. Solids 54 (1993) 1725. A.V. Eletsky, B.M. Smirnov, Usp. Fiz. Nauk 163 (1993) 33; A.V. Eletsky, B.M. Smirnov, Phys. Usp. 36 (1993) 202. H.W. Kroto, J.E. Fischer, D.E. Cox (Eds.), The Fullerenes Pergamon, Oxford, 1993. W. Kra¨tschmer, Nuovo Cimento A 107 (1994) 1077. V. Buntar, Phys. Lett. A 184 (1993) 131. A. Cheng, M.L. Klein, C. Caccamo, Phys. Rev. Lett. 71 (1993) 1200. M.C. Abramo, C. Caccamo, J. Phys. Chem. Solids 57 (1996) 1751.
39
[28] M.C. Abramo, C. Caccamo, D. Costa, et al., Phys. Rev. E 69 (2004) 031112. [29] L.A. Girifalco, Phys. Rev. B 52 (1995) 9910. [30] V.I. Zubov, N.P. Tretiakov, J.F. Sanchez, A.A. Caparica, Phys. Rev. B 53 (1996) 12 080. [31] V.I. Zubov, J.F. Sanchez-Ortiz, N.P. Tretiakov, I.V. Zubov, Phys. Rev. B 55 (1997) 6747. [32] R.C. Shukla, E.R. Cowley, Phys. Rev. B 31 (1985) 374. [33] D.J. Lacks, R.C. Shukla, Phys. Rev. B 54 (1996) 3266. [34] Z.W. Salsburg, W.W. Wood, J. Chem. Phys. 37 (1962) 798. [35] F.H. Ree, A.C. Holt, Phys. Rev. B 8 (1973) 826. [36] K. Westera, E.R. Cowley, Phys. Rev. B 11 (1975) 4008. [37] E.R. Cowley, J. Gross, Z. Gong, G.K. Horton, Phys. Rev. B 42 (1990) 3135. [38] A.C. Holt, M. Ross, Phys. Rev. B 1 (1970) 2700. [39] E. Wasserman, L. Stixrude, Phys. Rev. B 53 (1996) 8296. [40] J.X. Sun, H.C. Yang, Q. Wu, L.C. Cai, J. Phys. Chem. Solids 63 (2002) 113. [41] Y. Wang, D. Chen, X. Zhang, Phys. Rev. Lett. 84 (2000) 3220. [42] Y. Wang, Phys. Rev. B 62 (2000) 196. [43] Y. Wang, Phys. Rev. B 63 (2001) 245108. [44] Y. Wang, R. Ahuja, B. Johansson, Phys. Rev. B 65 (2001) 014104. [45] N.K. Bhatt, A.R. Jani, P.R. Vyas, V.B. Gohel, Physica B 65 (2001) 014104. [46] N.K. Bhatt, P.R. Vyas, A.R. Jani, V.B. Gohel, J. Phys. Chem. Solids 66 (2005) 797. [47] J.X. Sun, L.C. Cai, Q. Wu, F.Q. Jing, Phys. Rev. B 71 (2005) 024107. [48] N.X. Chen, Z.D. Chen, Y.C. Wei, Phys. Rev. B 55 (1997) R5.
Physica B 381 (2006) 34–39 www.elsevier.com/locate/physb
Analytic equation of state for FCC C60 solid based on analytic mean-field potential approach Sun Jiu-xun Department of Applied Physics, University of Electronic, Science and Technology, Chengdu 610054, PR China Received 21 October 2005; received in revised form 30 November 2005; accepted 1 December 2005
Abstract The analytic mean-field approach (AMFP) was applied to the FCC C60 solid. For the intermolecular forces the Girifalco potential has been utilized. The analytic expressions for the Helmholtz free energy, internal energy and equation of state have been derived. The numerical results of thermodynamic quantities are compared with the molecular dynamic (MD) simulations and the unsymmetrized selfconsistent field approach (CUSF) in the literature. It is shown that our AMFP results are in good agreement with the MD data both at low and high temperatures. The results of CUSF are in accordance with the AMFP at low temperature, but at high temperature the difference becomes prominent. Especially the AMFP predicted that the FCC C60 solid is stable upto 2202 K, the spinodal temperature, in good agreement with 2320 K from the MD simulation. However, the CUST just gives 1916 K, a temperature evidently lower than the MD data. The AMFP qualifies as a useful approach that can reasonably consider the anharmonic effects at high temperature. r 2006 Elsevier B.V. All rights reserved. PACS: 64.10.+h; 61.48.+c; 50.70.Ce; 62.50.+p Keywords: Equation-of-state; Fullerene; Thermodynamic properties; Anharmonicity
1. Introduction Since their discovery [1] and especially after their preparation in solid state [2,3], fullerenes have been the focus of attention of many scientific groups [3–31]. In the research of thermodynamic properties of fullerites, the Girifalco potential [4] has been extensively applied [5–31]. The C60 molecule consists of 60 carbon atoms symmetrically arranged on a sphere. Girifalco [4] pointed out that the molecules rotate rather freely and therefore can be approximated by spheres for many purposes. This is true even in the solid, at least for temperatures above 300 K [4]. Under these conditions, the hollow molecular cages can be assimilated to spheres whose surface consists of a uniform distribution of carbon sites. The overall interaction between two fullerene particles is then obtained by an integral of the interaction between pairs of sites on different cages. Using the standard (12-6) Lennard–Jones Tel.: +86 28 83202586.
E-mail address: [email protected]. 0921-4526/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2005.12.268
(LJ) function for the atom–atom interaction between two molecules and averaging, Girifalco obtained following intermolecular potential for the C60 fullerene [4]: 1 1 2 ðsÞ ðr=dÞ ¼ 2 þ sðs 1Þ9 sðs þ 1Þ9 s10 1 1 2 1 þ , ð1Þ sðs 1Þ3 sðs þ 1Þ3 s4 1 ¼ N 20 A=12d 6 ;
2 ¼ N 20 B=90d 12 ;
A ¼ 32:0 1060 erg cm6 ; d ¼ 0:71 nm,
d ¼ 0:71 nm,
(2)
B ¼ 55:77 10105 erg cm6 , ð3Þ
where s ¼ r=d, and N 0 ¼ 60. After the Girifalco potential was proposed, it has been extensively used in many works involving physical–chemical properties of C60 system [5–31]. The properties of fluid phase have been extensively researched by using computer simulations [5–31] and model equations of state [14–20]. For example, the molecular dynamics (MD) simulations on the properties
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
of fluid phase have been reported by Ohno et al. [5], Cheng et al. [6] and Rey et al. [7]. The Monte Carlo (MC) simulations have been reported by Hagen et al. [8,9], Caccamo et al. [10,11] and Hasegwa and Ohno [12,13]. At the same time, the density functional theory (DFT) [14–16], modified hypernetted-chain theory [17], hierarchical reference theory [18], Percus–Yevick approximation [19] of the Ornstein–Zernike equation and a variational series mean spherical approximation [20] have been used to research the thermodynamic properties of C60 fluid by the aid of double Yukawa fitting of Girifalco potential. However, for the solid state, the vast majority of experimental and theoretical work on pure and doped solid fullerenes (fullerites and fullerides) deals with their magnetic properties and superconductivity [21–25]. Up to now, just several works have been reported for thermodynamic properties [26–31]. Cheng et al. [26] predicted the phase diagram by using an integral equation approach combined with MD simulations based on the Girifalco potential. Abramo group published MD simulation results [27,28]. As for model equations of state, we just found two works. Girifalco firstly researched thermodynamic properties of FCC C60 under harmonic Debye–Gru¨neisen (DG) theory [29]. Zubov et al. [30,31] pointed out that the anharmonicity of the lattice vibrations of the C60 fullerite is strong at T4700 K. As this temperature is reached, the main anharmonic terms are no longer small corrections and must be included. Zubov et al. [30,31] considered anharmonic modification up to fourth order by using their unsymmetrized self-consistent field (CUSF) approach. However, Shukla et al. [32,33] have shown that the anharmonic expansion is divergent at high temperature. Although Zubov et al. [30,31] improved the convergence of the CUSF approach by using a self-consistent procedure and including anharmonic modification up to fourth order, this may not be enough at higher temperature. For example, the MD simulations [26] by using the Girifalco potential show that FCC C60 solid is stable up to about 2300 K. However, the CUSF predicted the stable temperature just is 1916 K. Thus, it is necessary to use other approaches that can include anharmonic effect soundly to study the system. The free volume theory (FVT) is a mean-field approximation to the thermal contribution of atoms to the Helmholtz free energy of crystalline phases [34–40]. The model assumes that each atom is confined to the Weigner–Seitz cell formed by its nearest neighbors. Although interatomic correlations are neglected, the cell model soundly includes anharmonic terms which are important for high temperatures, but which are ignored in DG model and quasiharmonic lattice dynamics. Many authors have demonstrated that the cell model matches successfully the thermodynamic properties of the FCC LJ crystal [34–36], sodium chloride [37] and exponential-6 model solid [38] calculated from Monte Carlo simulations. Wasserman and Stixrude further applied the cell model to a metallic solid iron [39]. The calculated properties are in
35
good agreement with available static and shock-wave experimental measurements. Sun et al. [40] extended the FVT to a generalized LJ solid including the quantum modification. The numerical results are in good agreement with experimental data of solid xenon, both at low and high temperature. Recently, Wang et al. [41–44] proposed an analytic mean-field potential (AMFP) approach, and applied it to many materials. Bhatt et al. [45,46] further applied the AMFP to lead and alkali metals, and concluded that in comparison with other theoretical models the AMFP is computationally simple, physically transparent and reliable to study the thermodynamic properties in the high pressures and high temperatures environment. However, the approach has been proven equivalent to the FVT, and just is an analytic approximation of FVT [47]. 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. However, in some cases the mean-field integral and the equation of state for the strict FVT are fairly complicated or cannot be analytically derived, the AMFP approach becomes useful. The Girifalco model for the FCC C60 solid just is the case, so we would develop analytic equation of state for the Girifalco model solid by using the AMFP approach. As was pointed out by Girifalco [29], in C60 molecule, the atoms are held together by strong chemical bonds while the interactions among these molecules in the FCC solid arise from van der Waals-type forces are much weaker. To a good approximation, therefore, the thermal properties of FCC C60 can be treated as the sum of those resulting from intermolecular vibrations and those arising from intramolecular vibrations [29]. Here we just consider the intermoleculear contributions, as having been done by Girifalco [29]. In Section 2 we derive analytic equation of state based on the AMFP approach. In Section 3 the numerical results are compared with results of Zubov et al. and MD data of Abramo et al. In Section 4 the conclusion is presented.
2. Analytic equation of state In terms of the FVT, the free energy can be expressed as [34–40] F 3 uð0Þ ¼ lnð2pmkT=h2 Þ þ ln vf , (4) NkT 2 2kT where m is mass of C60 molecule, uð0Þ is the potential energy of a molecule as the lattice is static, vf is the free volume. X Ri X a X uð0Þ ¼ zi z i di zi ðdi yÞ, (5) ¼ ¼ d d ia0 ia0 ia0 Z
rm
vf ¼ 4p 0
exp½gðr; V Þ=kTr2 dr,
(6)
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
36
where Ri ¼ di a is the distance of molecules in ith shell from the central molecule ði ¼ 0Þ. zi and di are structural constants, the values for FCC structure have been given in Ref. [48]. gðr; V Þ is the potential energy of a molecule as it roams from the center to a distance r. In terms of the AMFP approach [41–47], gðr; V Þ can be expressed by the static energy E c ðaÞ of a molecule: gðr; V Þ ¼ 12½ð1 þ r=aÞE c ða þ rÞ þ ð1 r=aÞE c ða rÞ 2E c ðaÞ,
(8)
where rm is the largest displace. For molecule with hardcore of diameter d, rm can be expressed as d a d 1 rm ¼ rw , (9) 2 2 a 1
where rw ð3a3 =4pgÞ3 a=2, is the Wigner–Seitz radius. a the nearest-neighbor distance, and pffiffiffi g is structure constant, for FCC structure it equals to 2. The volume of the FCC solid V ¼ Na3 =g. For simplicity, we introduce dimensionless reduced free volume v¯ f , reduced volume y, and reduced radial coordinate x: 3
vf ¼ 4pa v¯ f ¼ 4pgV v¯ f ,
x ¼ r=a;
(10) vd ¼ Nd 3 =g.
(11)
rm 1 1 1 . xm ¼ 2 y a
(12)
The reduced free volume v¯ f and its derivatives with respect to temperature and reduced volume can be expressed as Z xm v¯ f ¼ exp½gðx; yÞ=kTx2 dx, (13) 0
q 1 v¯ f ¼ v¯ fa ¼ T qT kT
Z
xm
exp½gðx; yÞ=kTgðx; yÞx2 dx,
0
(14) v¯ fb ¼
q 1 v¯ f ¼ qy kT
Z
xm
exp½gðx; yÞ=kT 0
q gðx; yÞx2 dx, qy (15)
1 xm 2 v¯ fc ¼ exp½gðxm ; yÞ=kT 0. 2 y
(16)
1X gðx; yÞ gðr; V Þ zi ½ð1 þ xÞðdi y þ di yxÞ 4 ia0 ð17Þ
ð18Þ
1 9 1 þ þ 9 10 2 sðs 1Þ s ðs þ 1Þ9 1Þ 9 20 1 3 þ 11 þ 1 þ 10 3 2 s sðs þ 1Þ s ðs 1Þ sðs 1Þ4 1 3 8 þ þ . s2 ðs þ 1Þ3 sðs þ 1Þ4 s5
0 ðsÞ ¼ 2
s2 ðs
ð19Þ
The compressibility factor and internal energy can now be derived: PV y q F Pc V Pf V ¼ ¼ þ , (20) Z¼ NkT 3 qy NkT NkT NkT Pc V y q ¼ uð0Þ, NkT 6kT qy
(21)
X q uð0Þ ¼ zi di 0 ðdi yÞ, qy ia0
(22)
Pf V y q y y¯vfb ¼1þ v¯ f ¼ 1 ð¯vfb v¯ fc Þ 1 , NkT 3¯vf qy 3¯vf 3¯vf (23) U q F 3 uð0Þ T q¯vf 3 uð0Þ v¯ fa ¼ T ¼ þ þ þ ¼ þ , NkT qT NkT 2 2kT v¯ f qT 2 2kT v¯ f (24) where Pc is the cold pressure, Pf is the thermal pressure deduced from free volume. By using above equations, all other thermodynamic quantities can be analytically derived. The derivations are straightforward. However, the expressions for thermal expansivity, 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 y and T, Z ¼ Zðy; TÞ. In terms of the function, formulas for thermal expansivity, compressibility coefficient and isochoric heat capacity can be reduced to following form: 1 qV 3 qy a¼ ¼ V qT P y qT P " # yqZ 1 Z qZ þ Z ¼ , ð25Þ T qT y 3 qy T qV 3 qy ¼ qP T y qP T 1 3 Vdy y qZ ¼ , Z 3 qy T NkT CV 1 qU U q U þT ¼ ¼ , qT NkT y Nk Nk qT V NkT b¼
Here gðx; yÞ gðr; V Þ, combining Eqs. (5), (8) and (10), we have
þ ð1 xÞðdi y di yxÞ 2ðdi yÞ,
þ ð1 xÞ0 ðdi y di yxÞ 20 ðdi yÞ,
ð7Þ
1 1X zi ðdi yÞ, E c ðaÞ ¼ uð0Þ ¼ 2 2 ia0
y ¼ a=d ¼ ðV =V d Þ1=3 ;
q 1X gðx; yÞ zi di ½ð1 þ xÞ0 ðdi y þ di yxÞ qy 4 ia0
1 V
ð26Þ
(27)
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
37
for many problems far from the critical point. Since the thermodynamic property of solids is a problem far from the critical point, it is not surprising that the AMFP can give good results [36–46]. However, the CUST is an approach to consider the anharmonic effects by using the method of series expansion [30,31]. Shukla et al. [32,33] have systematically studied the convergence of the anharmonic expansion, it has been shown that the numerical results would oscillate for odd and even order expansions, the convergence is acceptable at low temperature, slows down at slightly high temperature, and diverges at high temperature. In order to improve the convergence of the anharmonic expansion, Zubov et al. [30,31] developed their CUST approach by introducing consistent procedure and retaining the anharmonic terms up to fourth orders. Although the CUST indeed improves the convergence of anharmonic expansion, it may not be enough at high temperature as suggested by the numerical results in present work. As for the correlation effect, it has been shown contributing little to the thermodynamic properties for realistic potential models [36,37], and can be neglected in most cases except for the critical phenomena. In fact, both AMFP and CUST have neglected the correlation effect, so the correlation effect is not responsible to the discrepancies of the two approaches at high temperatures. In Fig. 1, the results for the internal energy and the pressure as functions of the temperature at various lattice spacings are compared with the MD results [28]. The isothermal curves at 300 K and 1900 K are plotted in Fig. 2. The two figures show that our results are in fairly good agreement with MD simulations [28]. Fig. 3 gives the variation of isochoric heat capacity C V , thermal expansivity a and bulk modulus BT versus pressure. The figure shows that C V is an increasing function of pressure, and the variation is slow. The a is a decreasing function of pressure, at low pressure the variation is fast, but at high pressure, the variation slows down and shows some saturation effect. The BT is a linear increasing function of pressure, but the slope
3. Results and discussion The thermodynamic properties calculated by using formalism developed above are listed and compared with the CUSF approach of Zubov et al. [30] in Table 1. The spinodal point T s is the temperature satisfying the condition, BT ðT s Þ ¼ 0. The system is unstable for temperatures above T s . The MD simulations performed by Cheng group [26] and Abramo group give T s ¼ 2320 K [27]. The CUSF approach gives T s ¼ 1916 K [30], a value lower than the MD result about 20%. The reason maybe that the CUSF approach cannot consider anharmonic effects soundly. Although it can consider anharmonic terms up to fourth order, it is not enough at high temperature. So Zubov et al. [30] just computed thermodynamic properties up to 1900 K. However, our calculation gives T s ¼ 2202 K, is in good agreement with MD simulations. So our calculations have been done up to 2200 K. The results in this paper verify the viewpoint that the AMFP is a useful approach to consider anharmonic effects at high temperatures [32–48]. Table 1 shows that our results are in agreement with the CUST of Zubov et al. at low temperature, at high temperature the difference becomes prominent. Especially near the sinodal temperature of CUSF 1900 K, the thermal expansivity coefficient becomes divergent, and the compressibility coefficient tends zero. However, both coefficients take finite values at this temperature, as the spinodal point of AMFP is 2202 K. Although some difference exists for the linear thermal expansivity at ambient temperature, our calculation gives 1:148 105 K1 at 300 K, is in slightly better agreement with the MD data [28] 1:137 105 K1 than the CUSF result 1:084 105 K1 [30]. Now we explain in more detail the different features of CUSF and AMFP responsible for the discrepancies at high temperatures. The AMFP is a mean field approximation to the thermal contribution of atoms to the Helmholtz free energy of crystalline phases [36–46]. It is well known that the mean field approximation can give fairly good results
Table 1 Properties of the FCC phase of C60 T
261.4
400
600
800
1000
1200
1400
1600
1800
a
1.0068 1.0067
1.0084 1.0082
1.0108 1.0106
1.0136 1.0132
1.0166 1.0161
1.0201 1.0194
1.0241 1.0233
1.0290 1.0282
1.0354 1.0354
a
1.130 1.070
1.198 1.126
1.313 1.223
1.456 1.347
1.640 1.512
1.884 1.745
2.227 2.103
2.750 2.763
3.667 4.710
1900 1.0394 1.0426 4.469 12.713
2000
2100
1.0445
1.0517
5.853
9.073
BT
136.04 140.51
124.19 128.89
107.87 112.88
92.40 97.54
77.71 82.73
63.77 68.22
50.53 53.71
37.95 38.61
25.95 21.11
20.11 7.40
14.33
8.41
CV
23.50
23.33
22.41
22.90
21.87
21.30
20.69
20.02
19.26
18.83
18.32
17.69
2200 1.0692 71.86 0.857 16.49
The nearest-neighbor distance in nm, linear thermal expansion coefficient in 105 K1 , the heat capacity in kJ mol1 K1 , the bulk modulus in kbar. The values in the first lines have been calculated in this work, in the second lines are from Zubov et al. [30].
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39
38 100
27 CV (kJ/mol K)
P (kbar)
80 60 40 20 0 200
400
600
800
1000 1200 T (K)
1400
1600
1800
26 25 24
2000
0
100
200 P (kbar)
300
400
0
100
200 P (kbar)
300
400
0
100
200 P (kbar)
300
400
-50
0.9 -100 -150 -200 200
400
600
800
1000 1200 T (K)
1400
1600
1800
2000
α (10-5K-1)
U (kJ/mol)
0
Fig. 1. Pressure (top) and internal energy (bottom) as functions of temperature for several lattice spacing a. The symbols refer to MD simulation data for a ¼ 1:36 nm (circles), 1.38 nm (triangles), and 1.42 nm (squares) [28]. The lines are corresponding numerical results calculated in this work.
P (kbar)
300
6000 4000 2000
T = 300 K 200
0 100 0 1.45
1.5
1.55
1.6
1.65
1.7
1.75
ρ (nm-3) 400 300 P (kbar)
0.3 0
BT (kbar)
400
0.6
T = 1900 K 200 100 0 1.45
1.5
1.55
1.6
1.65
1.7
1.75
ρ (nm-3)
Fig. 2. Variation of pressure versus density at two temperatures, 300 K (top) and 1900 K (bottom). The circles refer to MD simulation data [28]. The lines are corresponding numerical results calculated in this work.
Fig. 3. Variation of thermodynamic quantities versus pressure calculated in this work. The heat capacity is in kJ mol1 K1 , linear thermal expansion coefficient in 105 K1 , the bulk modulus is in kbar.
monic effects available at present, can be applied to high temperature, our comparison shows that some discrepancy exists at high temperature between the CUST with MD data and our FVT results. However, the numerical results from the AMFP are in good agreement with the MD simulations available. This further verifies the viewpoint that the AMFP is a useful approach to consider the anharmonic effects at high temperature. Acknowledgment
is fairly large and the increase is fast. This means the Girifalco model solid for the FCC C60 is difficult to be compressed. 4. Conclusion In summary, the AMFP has been applied to the Girifalco model solid of C60. All thermodynamic quantities can be analytically evaluated based on the formalism developed in this paper. Although Zubov et al. claimed their CUST approach, another theory to consider anhar-
This work was supported by the Natural Science Foundation of China under Grant no. 19904002, by the program of middle-aged and young academic cadres of UESTC. References [1] H.W. Kroto, J.R. Heath, S.C. O’Brien, et al., Nature 318 (1985) 162. [2] W. Krat¨schmer, L.D. Lamb, K. Fostinoupolos, R.D. Huffman, Nature 347 (1990) 354. [3] O. Zhou, D.E. Cox, J. Phys. Chem. Solids 53 (1992) 1373.
ARTICLE IN PRESS S. Jiu-xun / Physica B 381 (2006) 34–39 [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]
L.A. Girifalco, J. Phys. Chem. 96 (1992) 858. K. Ohno, Y. Maroyama, Y. Kawazoe, Phys. Rev. B 52 (1996) 4078. A. Cheng, M.L. Klein, C. Caccamo, Phys. Rev. Lett. 71 (1993) 1200. C. Rey, L.J. Gallego, J.A. Alonso, Phys. Rev. B 49 (1994) 8491. M.H.J. Hagen, E.J. Meijer, G.C.A.M. Mooij, et al., Nature 365 (1993) 425. M.H.J. Hagen, D. Frenkel, J. Chem. Phys. 101 (1994) 4093. C. Caccamo, D. Costa, A. Fucile, J. Chem. Phys. 106 (1997) 255. D. Costa, G. Pellicane, M.C. Abramo, C. Caccamo, J. Chem. Phys. 118 (2003) 304. M. Hasegawa, K. Ohno, J. Chem. Phys. 111 (1999) 5955. M. Hasegawa, K. Ohno, J. Chem. Phys. 113 (2000) 4315. M. Hasegawa, K. Ohno, J. Phys.: Condens. Matter 9 (1997) 3361. L. Mederos, G. Nawascue´s, Phys. Rev. B 50 (1994) 1301. M. Hasegawa, K. Ohno, Phys. Rev. E 54 (1996) 3928. C. Caccamo, Phys. Rev. E 51 (1995) 3387. M. Tau, A. Parola, D. Pini, L. Reatto, Phys. Rev. E 52 (1995) 2644. H. Guerin, J. Phys.: Condens. Matter 10 (1998) L527. M. Khedr, M.S. Al-Busaidy, S.M. Osman, J. Phys.: Condens. Matter 17 (2005) 4411. J.F. Fischer, P.A. Heiney, J. Phys. Chem. Solids 54 (1993) 1725. A.V. Eletsky, B.M. Smirnov, Usp. Fiz. Nauk 163 (1993) 33; A.V. Eletsky, B.M. Smirnov, Phys. Usp. 36 (1993) 202. H.W. Kroto, J.E. Fischer, D.E. Cox (Eds.), The Fullerenes Pergamon, Oxford, 1993. W. Kra¨tschmer, Nuovo Cimento A 107 (1994) 1077. V. Buntar, Phys. Lett. A 184 (1993) 131. A. Cheng, M.L. Klein, C. Caccamo, Phys. Rev. Lett. 71 (1993) 1200. M.C. Abramo, C. Caccamo, J. Phys. Chem. Solids 57 (1996) 1751.
39
[28] M.C. Abramo, C. Caccamo, D. Costa, et al., Phys. Rev. E 69 (2004) 031112. [29] L.A. Girifalco, Phys. Rev. B 52 (1995) 9910. [30] V.I. Zubov, N.P. Tretiakov, J.F. Sanchez, A.A. Caparica, Phys. Rev. B 53 (1996) 12 080. [31] V.I. Zubov, J.F. Sanchez-Ortiz, N.P. Tretiakov, I.V. Zubov, Phys. Rev. B 55 (1997) 6747. [32] R.C. Shukla, E.R. Cowley, Phys. Rev. B 31 (1985) 374. [33] D.J. Lacks, R.C. Shukla, Phys. Rev. B 54 (1996) 3266. [34] Z.W. Salsburg, W.W. Wood, J. Chem. Phys. 37 (1962) 798. [35] F.H. Ree, A.C. Holt, Phys. Rev. B 8 (1973) 826. [36] K. Westera, E.R. Cowley, Phys. Rev. B 11 (1975) 4008. [37] E.R. Cowley, J. Gross, Z. Gong, G.K. Horton, Phys. Rev. B 42 (1990) 3135. [38] A.C. Holt, M. Ross, Phys. Rev. B 1 (1970) 2700. [39] E. Wasserman, L. Stixrude, Phys. Rev. B 53 (1996) 8296. [40] J.X. Sun, H.C. Yang, Q. Wu, L.C. Cai, J. Phys. Chem. Solids 63 (2002) 113. [41] Y. Wang, D. Chen, X. Zhang, Phys. Rev. Lett. 84 (2000) 3220. [42] Y. Wang, Phys. Rev. B 62 (2000) 196. [43] Y. Wang, Phys. Rev. B 63 (2001) 245108. [44] Y. Wang, R. Ahuja, B. Johansson, Phys. Rev. B 65 (2001) 014104. [45] N.K. Bhatt, A.R. Jani, P.R. Vyas, V.B. Gohel, Physica B 65 (2001) 014104. [46] N.K. Bhatt, P.R. Vyas, A.R. Jani, V.B. Gohel, J. Phys. Chem. Solids 66 (2005) 797. [47] J.X. Sun, L.C. Cai, Q. Wu, F.Q. Jing, Phys. Rev. B 71 (2005) 024107. [48] N.X. Chen, Z.D. Chen, Y.C. Wei, Phys. Rev. B 55 (1997) R5.