Quasiharmonic approach to a second-order phase transition

PHYSICAL REVIEW B 70, 104109 (2004)

Quasiharmonic approach to a second-order phase transition J. Łażewski,* P. T. Jochym, P. Piekarz, and K. Parlinski Institute of Nuclear Physics, Polish Academy of Sciences, Radzikowskiego 152, 31-342 Kraków, Poland (Received 17 December 2003; revised manuscript received 23 April 2004; published 17 September 2004) We present an ab initio based method of constructing the phase diagram of the second-order phase transitions with a significant volume dependence of the soft mode. This approach is based on the assumption that in the investigated crystal the phonon frequencies are mainly dependent on crystal volume and geometry. A critical volume for the transition at T = 0 K is derived from the crystal symmetry breaking point. The quasiharmonic approximation is used to find the free energy and the equation of state at finite temperature. The phase transition line is then the locus of points in the p–T diagram corresponding to the critical volume. As an example we determine the p–T phase diagram for CaCl2, which undergoes a rutile type—orthorhombic secondorder phase transition. DOI: 10.1103/PhysRevB.70.104109

PACS number(s): 64.70.Kb, 61.50.Ks, 78.30.Hv, 63.20.Dj

I. INTRODUCTION

In a first-order phase transition the free-energy curves FA and FB of two phases A and B cross at a phase transition T0 and volume V0. The continuation of free energies on either side of the transition point represent metastable phases, which correspond to overheated or undercooled states. An estimate of the transition point can then be obtained by computing the free energies including the extension to the metastable regions, using for example the quasiharmonic approximation. A number of papers have dealt with the ab initio determination of thermodynamic properties, phases stability limits, phase transitions, and phase diagrams for crystals with first-order transitions, among them GeO2,1 AgBr,2 NiTi,3 and CuInSe2.4 Combined with the ab initio technique, the quasiharmonic approximation was successfully used to study thermal properties of simple metals,5,6 semiconductors,7 and more complex, geologically important materials, such as MgO,8 MgSiO3,9 Al3Li,10 and Mg2SiO4.11 The quasiharmonic approach was also applied to study the temperature-driven structural transition of tin12 and boron nitride13 as well as the melting curve of the face-centeredcubic phase of aluminium.14 In a continuous (second-order) phase transition the situation is rather different.15–17 The phases do not exist in metastable forms, which means that no overheating, or undercooling is possible. Hence, the transition point is determined either by detecting a symmetry breaking, or by singularities observed in second-order derivatives of the free energy, caused by critical fluctuations. Using ab initio methods it is easy to find the symmetry breaking point at T = 0 K, but at finite temperature it is hardly possible without additional assumptions. Generally, there exist a number of kinds of second order phase transitions (displacive, order-disorder, magnetic). Some of them occur owing to volume changes under pressure or/and with increased temperature. However, thermal expansion should not be considered as a transition mechanism but rather as an accompanying effect of the phase transition, an external agent evoking the phase transition. In this paper we propose a way to determine the phase diagram of a crystal exhibiting a second-order phase transition connected 1098-0121/2004/70(10)/104109(7)/$22.50

with significant volume changes. The approach is based on the observation that in the harmonic crystal at constant volume regime the phonon frequencies do not alter. Therefore, the symmetry breaking points could occur in the same critical volume Vc regardless of temperature. Using this idea, together with first-principles densityfunctional-theory (DFT) calculations, we describe the rutile type—orthorhombic phase transition in calcium dichloride, CaCl2. At T = 0 K this phase transition can be induced by applying negative pressure. We compute the phonon spectra in both phases, tetragonal and orthorhombic, and study the soft mode behavior as a function of pressure in the vicinity of the critical point. Using the quasiharmonic approximation, we calculate the total free energy up to 1000 K and determine the finite-temperature equation of state. Finally, the assumption of constant phonon frequencies at constant crystal volume regime allows us to draw the p–T phase diagram. II. CONSTANT-CRITICAL-VOLUME QUASIHARMONIC APPROACH

The thermodynamic properties of the crystals are to a large extent determined by phonons. Such properties can be satisfactorily described in general within the quasiharmonic approximation, in which a change of crystal volume due to finite temperature is mapped to the change of crystal volume at T = 0 K (typically as a function of pressure). The thermodynamic functions are calculated using the formulae for the harmonic crystal, but anharmonic effects are nevertheless included, at least partly, through the volume (pressure) dependence of phonon frequencies. When the crystal exhibits a second-order phase transition, one should describe the critical volume at which the symmetry of the crystal changes. With an ab initio method at T = 0 K one may easily determine it. For example, two hitherto equal lattice parameters become different upon a change in volume, or a 90° angle starts to change. Another opportunity would be to observe a volume Vc at which a given mode’s frequency vanishes, ␻soft = 0. Since the practical application of ab initio methods is limited to T = 0 K, these criteria cannot be used for finite temperatures without additional assumption.

70 104109-1

©2004 The American Physical Society

PHYSICAL REVIEW B 70, 104109 (2004)

ŁAŻEWSKI et al.

In the strictly harmonic potential the phonon frequencies do not vary when the crystal volume, symmetry, lattice parameters and atomic positions remain constant, even if the temperature is changing. Indeed, the temperature influences the occupation of harmonic phonon states and the phonon amplitudes. Thus, at any temperature, but at constant volume and geometry of the crystal, harmonic phonon frequencies remain the same. Such properties are certainly obeyed further away from the phase transition point at T = 0 K and at not too high temperatures. At higher temperatures anharmonic effects not included in quasiharmonic approximation will enhance frequency renormalization. We would like to suggest that one may trace the symmetry breaking point of the phase transition from T = 0 K to finite temperatures along the constant volume regime as well. Taken together, these criteria mean that the phase transition instability corresponds to some definite critical volume Vc independent of temperature. Therefore, the task of computing the phase diagram starts by determining the critical volume Vc by means of T = 0 K first-principles calculations. Then the equation of state p共V , T兲 is obtained using the quasiharmonic approximation, and formally inverted to give V = V共p , T兲. The phase transition line in the p–T diagram is then determined from condition Vc = V共pc , Tc兲 = constant. The critical volume Vc can be derived from the soft-mode behavior. However, in a pseudoproper ferroelastic phase transition the Brillouin zone-center optic soft mode does not completely soften due to the linear coupling with elastic shear of the same symmetry and causes acoustic instability at the phase transition. Therefore, the acoustic soft mode or elastic shear vanishing at the symmetry breaking point can be used as an indicator of the phase transition. Often the critical point can be also easily determined from changes of structural parameters when—due to symmetry breaking— lattice constants cease to be identical or some atoms no longer occupy high symmetry points.

III. RUTILE-TYPE—ORTHORHOMBIC PHASE TRANSITION IN CaCl2 A. Present knowledge

Calcium dichloride 共CaCl2兲 is a ferroelastic material which, at ambient pressure, exists in two phases: a hightemperature (above Tc = 491 K) tetragonal (rutile-type, space group P42 / mnm) and a low-temperature orthorhombic phase (space group Pnnm).18 The orthorhombic phase is a result of a deformation which involves the rotation of nearly-rigid Cloctahedra around the tetragonal axis. The phase transition is characterized by a one-dimensional order parameter of symmetry B1g, and an optical mode of the same symmetry shows a strong softening in the tetragonal phase.18 In addition an elastic instability follows from the bilinear coupling of the acoustic shear mode B1g to the order parameter. This is an example of a pseudoproper ferroelastic phase transition. The transition in CaCl2 was studied with ab initio methods by Válgoma et al.19,20 For the rutile-type phase, their calculations confirmed the instability of the B1g optic mode at the ⌫ point. All other zone-center modes turned out to

have real frequencies and agreed fairly well with the experimental data, and the lattice parameters of the tetragonal phase agreed with the experimental values extrapolated to 0 K. The energy map as a function of the amplitude of the B1g mode and the orthorhombic strain were also obtained and interpreted in the framework of a Landau free energy expansion. Many other compounds with the general formula AB2 exist in the rutile-type phase at 0 K and transform into the orthorhombic structure at high pressure. To this class of compounds belong materials of geophysical importance, such as stishovite—high pressure phase of silica 共SiO2兲, and germanium dioxide, GeO2. The rutile-type—CaCl2-type phase transition was studied experimentally21,22 and theoretically.1,23,24 The local density approximation (LDA) was used to study SiO2.23 Analysis of the zone-center phonons showed the softening of the B1g mode in the rutiletype phase near the critical pressure 47 GPa. This soft mode transforms to an Ag mode in the orthorhombic phase, according to the group theory prediction. The coupling between the shear elastic constants and the B1g optic mode was studied in the rutile-type structure of SnO2,24 and it was shown that beyond the ⌫ point the B1g soft branch interacts with the transverse acoustic (TA) mode. As a result, part of the TA mode along ⌫-M direction becomes imaginary. This destabilizes the tetragonal crystal and leads to the orthorhombic phase. Łodziana et al.1 studied the pressure-induced phase transitions in GeO2 and confirmed the soft mode behavior in both the rutile-type and CaCl2-type phases. They found that the shear modulus of both phases goes to zero at about 19 GPa. B. Calculation details

The total energy calculations within DFT were performed using the VASP program.25 This program solves the KohnSham equations by an iterative diagonalization scheme. Ions were represented by ultrasoft pseudopotentials26 with 4s2 and 3s23p5 valence-electron configurations for Ca and Cl, respectively. The exchange-correlation energy functional was derived using the generalized gradient approximation (GGA).27 The electronic wave functions were represented in a planewave basis with 274 eV energy cutoff. Both the tetragonal and orthorhombic phases were optimized in 2 ⫻ 2 ⫻ 3 supercells containing 72 atoms (see later). Brillouin zone integrations were performed by sampling the wave functions at k points generated by the Monkhorst-Pack scheme28 with 2 ⫻ 2 ⫻ 2 grids. After the optimization of the electronic degrees of freedom, Hellmann-Feynman (HF) forces on all atoms and stresses were calculated. The minimization procedure ended when the residual forces were less than 0.01 meV/ Å. Phonon frequencies were calculated using the direct method,29,30 in which atomic force constants up to an appropriately large interaction range are derived from the HF forces. These forces are induced by small displacements of nonequivalent atoms from their equilibrium positions in supercells. The force-constant matrices are evaluated from the HF forces using the singular value decomposition method to

104109-2

PHYSICAL REVIEW B 70, 104109 (2004)

QUASIHARMONIC APPROACH TO A SECOND-ORDER…

solve an overdetermined linear system, and phonon frequencies are computed at each k point by the diagonalization of the corresponding dynamical matrix. Since the interaction range was confined to the interior of our large supercells we expect to obtain rather precise phonon frequencies for all wave vectors. One important feature of those supercells is that their shapes are nearly cubic (for the orthorhombic phase at p = 0 kb the size was: 12.95 Å ⫻ 12.74 Å ⫻ 12.49 Å). This ensures the same interaction range in all directions in the computation of the dynamical matrix. Finite systems, however, cannot exhibit the longitudinal-transverse (LO/TO) splitting of optical modes at k = 0, which results from the long-range polarization field induced by the atomic displacements. To extrapolate the LO phonon frequency at the ⌫ point we use elongated supercells (1 ⫻ 1 ⫻ 8 and 1 ⫻ 1 ⫻ 10) of the rutile-type structure.31 A more elegant method for the calculation of the LO/TO splitting involves the determination of the effective charges.32 We computed these in the tetragonal phase using the Berry-phase method33 as implemented in the SIESTA34 program. We obtained the following values of Ca Ca = 1.798, Z Ca the effective charge tensor: Z xx xy = 0.133, Z zz Cl Cl Cl = 1.936, Z xx = −0.901, Z xy = −0.135, and Z zz = −0.968. From the calculated LO/TO splitting and the effective charges, following the procedure described in Ref. 35, the electronic part of the dielectric constant ␧⬁ has been estimated to be 1.69. In further calculations we assumed the same values of effective charges and dielectric constant for both phases at all pressures, but in any case we found that the LO/TO splitting has negligible influence on the values of thermodynamic functions. The equation of state of the crystal can be obtained from the expression p=−

冋

⳵F共V,T兲 ⳵V

册

共1兲

, T

where the free energy F共V , T兲 contains two main contributions: the ground state energy of the crystal E共V兲 and the free energy of the phonons Fph共V , T兲: F共T,V兲 = E关V共T兲兴 + Fph共V,T兲.

共2兲

Any purely electronic contribution is neglected. The groundstate energy E varies according to temperature-induced changes in crystal volume and structure. In the (quasi)harmonic approximation the phonon free energy is written in the form ⬁

Fph共V,T兲 = kBT

冕 0

冋 冉 冊册

g共␻兲ln 2 sinh

␻ 2kBT

d␻ ,

共3兲

where ␻ denotes phonon frequency, g共␻兲 is the phonon density of states (implicitly depending on the volume), ប is the Planck constant, and kB is the Boltzmann constant. g共␻兲 is obtained by the summation over wave vectors k randomly selected from the first Brillouin zone. Equation (3) includes the zero-point contribution of phonons.

TABLE I. Lattice parameters of the orthorhombic CaCl2. The positions of Cl atoms are given in fractional coordinates.

a共Å兲 b共Å兲 c共Å兲 xCl y Cl

Present

Experiment18

Calculated20

6.474 6.368 4.162 0.316 0.291

6.446 6.167 4.137 0.325 0.275

6.429 6.054 4.088 0.347 0.255

C. Crystal structure

We begin the structural studies with the orthorhombic phase, which is stable at ambient conditions. This phase was previously considered by Válgoma et al.,20 whose calculated crystal parameters agree very well with the experimental values extrapolated to 0 K. In Table I, the lattice constants of our optimized crystal are compared with the experimental data and the results obtained in Ref. 20. The present lattice constants are larger than the values of Ref. 20, as a result of the different exchange-correlation functional used: the LDA, used in Ref. 20, gives usually too small lattice constants, while the GGA lattice constants are larger. The experimental points lie between the LDA and GGA values.36 In the same Table I we present also the Cl coordinates 关xCl , y Cl兴. The positions of cations are fixed by the crystal symmetry. In the next step, we applied an external pressure to the supercell. For each pressure the lattice constants and internal degrees of freedom were reoptimized. In Fig. 1, we plot the pressure dependence of the lattice constants, the fractional coordinates, volume, energy and enthalpy. Beyond the positive-pressure we consider also the negative-pressure regime, where the structural transformation is expected. The results show that at the T = 0 K critical point occurs at a small negative pressure. At this critical pressure (about −5 kb) the crystal relaxes spontaneously to the tetragonal symmetry. Because of the smaller unit cell volume in the orthorhombic phase, atoms are more densely packed in the orthorhombic than in the rutile-type phase. Consequently, the CaCl6 octahedra are more distorted in the orthorhombic phase, and a larger volume allows the Cl ions to reach the diagonal positions in the rutile-type structure. The amplitude of the optical mode connected with the rotation of octahedra is actually the order parameter characterizing this phase transition.20 The crystal structure of the tetragonal phase for p ⬍ pc was also obtained by the minimization procedure with the imposed P42 / mnm space group symmetry. In Fig. 1, these points are marked by filled squares. The full agreement between these two approaches is evident. An exact determination of pc is not possible due to the limited accuracy of the calculations, but it can be confidently stated that the transition takes place between −5 and −6 kb, a range which corresponds to volumes of 174.8 and 175.5 Å3 (volume changes in the vicinity of the critical point are clearly smooth to the transition). Also the ground state energy, E, and the enthalpy, H = E + pV, are continuous functions of pressure. For large negative pressures (below ps = −32 kb) the lattice constants attain very unrealistic values and the crystal structure becomes unstable.

104109-3

PHYSICAL REVIEW B 70, 104109 (2004)

ŁAŻEWSKI et al.

TABLE II. A comparison of calculated and experimental phonon frequencies (in THz) at the ⌫ point in the orthorhombic phase of CaCl2.

FIG. 1. Dependence of the lattice constants, fractional parameters, specific volume, ground state energy, and enthalpy on volume and pressure for the orthorhombic (open circles) and rutile-type (filled squares) phases of CaCl2. D. Phonon modes

Using the optimized orthorhombic crystal structure at p = 0 kb we have obtained the HF forces from six independent atomic displacements: three for Ca and three for the Cl atoms in each crystallographic direction. At the zone center, there are 15 optic modes described by the following irreducible representations:

Symmetry

Activity

Present

Experiment18

Ag Au B3u B2u B1g B2g B3g Ag B2u B3u B1u Au B3u B2u B1g

R

1.34 2.35 2.65 2.92 3.38 4.45 4.46 5.92 6.04 6.07 6.33 6.48 6.97 7.00 7.25

1.56

IR IR R R R R IR IR IR IR IR R

3.42 4.68 6.26

7.52

volume increases. In agreement with the experimental result and the previous theoretical studies, the Ag mode strongly softens in the vicinity of the phase transition. The computed frequencies of the Ag mode are presented in Table III. Across the transition, a part of the TA mode along the ⌫-S direction becomes unstable (plotted with negative values) as a result of the bilinear coupling of the optical soft mode and the elastic strain. Since the acoustic soft mode approaches zero before the B1g optic mode does, the softening of the optical mode is not complete. The frequency of the soft optic mode at the critical point was estimated using Landau theory,20 and the obtained value 0.51 THz is in good agreement with the experimental value 0.42 THz.18 Our calculations give an upper bound for the soft mode frequency: ␻ ⬍ 0.86 THz. Below

⌫ortho = 2Ag + 2Au + B1g + B2g + B3g + B1u + 3B2u + 3B3u . 共4兲 Six of them are Raman (R) active, seven infrared (IR) active, and two are silent. In Table II we compare the calculated frequencies at the ⌫ point with the experimental data taken from Raman scattering experiments.18 For most frequencies the agreement between theory and experiment is very good (differences are less than 5%). The largest deviation is observed for the lowest Ag mode and could be explained by its soft nature. The underestimation of all the phonon frequencies results from the systematic overestimation of the lattice constants in the GGA. In Fig. 2 (solid curves) the phonon dispersion is plotted along high symmetry directions. All phonon modes have real frequencies, confirming the stability of the orthorhombic phase. In order to study the soft mode behavior, we calculated phonon dispersion curves at several pressures. In Fig. 2, the dispersion curves at p = −7 kb are plotted with dashed lines. All phonon modes decrease their frequencies because the

FIG. 2. Calculated phonon dispersion curves for the orthorhombic phase of CaCl2 at p = 0 kb (solid lines) and p = −7 kb (dashed lines).

104109-4

PHYSICAL REVIEW B 70, 104109 (2004)

QUASIHARMONIC APPROACH TO A SECOND-ORDER… TABLE III. Frequencies of the soft optic mode as a function of pressure.

Pressure (kb)

Ag (THz) orthorhombic

B1g (THz) tetragonal

0.0 −5.0 −6.0 −12.0 −14.0 −15.0 −16.0 −18.0

1.34 0.86 0.91 ¯ ¯ ¯ ¯ ¯

¯ ¯ ¯ 1.58 1.87 1.97 2.05 2.09

−7 kb the crystal relaxes spontaneously to the tetragonal symmetry. In the rutile-type phase of CaCl2, the zone-center phonon frequencies were computed in Ref. 20, and it was demonstrated that the B1g mode is unstable (with imaginary frequency). One can expect that by increasing the crystal volume (or, equivalently, decreasing the pressure beyond the critical value), this instability vanishes. We have repeated calculations of phonon dispersion curves for the rutile-type structure in the interval from p = −12 to − 30 kb. (In this tetragonal symmetry, the HF forces were computed from four independent atomic displacements: two for Ca and two for Cl in x and z directions.) In Fig. 3 we plot the dispersion curves at p = −30 kb: as expected, all phonon modes, including the B1g soft mode, have real frequencies. E. The phase diagram

Using the quasiharmonic approximation we computed for all stable structures which do not show imaginary phonon frequencies, the total free energy as a function of volume for temperatures from 0 up to 1000 K. In Fig. 4 we fit a third-

FIG. 4. Dependence of free energy F on crystal volume for a number of selected temperatures. The bold line connects the F minima for the different temperatures. The crystal symmetry changes between the two dashed vertical lines, which represent then the lower and upper estimates of the critical volume.

order polynomial curve to the calculated points for the tetragonal and orthorhombic phases, and we indicate with a bold curve the locus of the free-energy minima at each temperature. These minima correspond, according to Eq. (1), to zero pressure. The bold curve is then a representation of the temperature dependence of the crystal volume at ambient pressure V共T , p = 0兲. The V共T , p兲 function for nonzero pressures can be obtained from the appropriate tangents to the F共V , T兲 curves, as in Eq. (1). The two dashed lines in Fig. 4 indicate the critical-volume range determined by the structural and soft-mode analysis considered in the preceding sections at T = 0 K. On the basis of the discussion presented in Sec. II we expect the critical volume Vc to remain constant along the second-order phase transition boundary, or at least to be confined between the upper and lower bounds shown in Fig. 4. We could thus evaluate the critical boundary Tc共p兲 between the orthorhombic and tetragonal phases from the equation Vc = V共Tc , p兲. We plot the phase diagram in Fig. 5. There are actually two lines, upper and lower bounds of the phase transition boundary, corresponding to the uncertainty we can draw the criticalvolume estimate. The experimental critical temperature at ambient pressure falls within the predicted boundaries. Notice the linear dependence of the critical temperature as a function of positive pressure. At higher temperatures the anharmonic effects could contribute more to the thermodynamic functions, and the linear Tc共p兲 dependence could change. In Fig. 5 we draw the phase boundary for negative pressure as well, but of course this interval is usually inaccessible to the experiment. IV. CONCLUSIONS

FIG. 3. Calculated phonon dispersion curves in the rutile-type phase of CaCl2 at p = −30 kb.

Phonon frequencies change with temperature due to volume changes and/or pure temperature renormalization inde-

104109-5

PHYSICAL REVIEW B 70, 104109 (2004)

ŁAŻEWSKI et al.

FIG. 5. The phase diagram of CaCl2. A filled square indicates the experimental critical point at T = 491 K.

proach to study a phase transition line. This method can be used, provided that we stay within the applicability limits of the quasiharmonic approximation. At T = 0 K the critical point is described by the singularities related to a breaking of crystal symmetry and/or by the vanishing of a soft-mode frequency. We have assumed that the harmonic phonon frequencies are only a function of the crystal volume (at fixed geometry), and depend on temperature only through the accompanying volume changes. This means that there is a single critical volume for the phase transition which is either constant, or at most very weakly dependent on temperature. On this basis we have evaluated from first-principles the upper and lower boundaries for the orthorhombic to tetragonal phase transition on the p–T phase diagram of calcium dichloride CaCl2. In the process we have obtained the phonon dispersion relations for both tetragonal and orthorhombic phases as a function of volume (pressure). We find good agreement of the lattice parameters, Raman frequencies, and phase diagram with the available experimental data.

pendent of volume. Together with temperature increase the contribution of the later process rises, however, for relatively small temperatures it is negligible for a broad class of materials. In the second-order phase transition with a soft mode, the critical temperature is directly connected with the phonon frequency renormalization in the vicinity of the critical point. For phase transitions in which volume agent plays a dominant role, we have proposed a constant-critical-volume ap-

The authors thank Alberto Garcia for fruitful discussions and critical reading of the manuscript. This work was partially supported by the EU Project No. NMP4-CT-2003001516 and the State Committee of Scientific Research (KBN), Grant No. 5 P03B 069 20.

*Author to whom correspondence should be addressed. Electronic

17 K.

address: [email protected] 1 Z. Łodziana, K. Parlinski, and J. Hafner, Phys. Rev. B 63, 134106 (2001). 2 P. T. Jochym and K. Parlinski, Phys. Rev. B 65, 024106 (2002). 3 K. Parlinski and M. Parlinska-Wojtan, Phys. Rev. B 66, 064307 (2002). 4 J. Łażewski, H. Neumann, K. Parlinski, G. Lippold, and B. J. Stanbery, Phys. Rev. B 68, 144108 (2003). 5 A. A. Quong and A. Y. Liu, Phys. Rev. B 56, 7767 (1997); J. Xie, S. Gironcoli, S. Baroni, and M. Scheffler, ibid. 59, 965 (1999). 6 S. Narasimhan and S. de Gironcoli, Phys. Rev. B 65, 064302 (2002). 7 S. Biernacki and M. Scheffler, Phys. Rev. Lett. 63, 290 (1989). 8 B. B. Karki, M. C. Warren, L. Stixtrude, G. J. Ackland, and J. Crain, Phys. Rev. B 55, 3465 (1997). 9 B. B. Karki, R. M. Wentzcovich, S. de Gironcoli, and S. Baroni, Phys. Rev. B 61, 8793 (2000). 10 Z. Li and J. S. Tse, Phys. Rev. B 61, 14531 (2000). 11 P. Piekarz, P. T. Jochym, K. Parlinski, and J. Łażewski, J. Chem. Phys. 117, 3340 (2002). 12 P. Pavone, S. Baroni, and S. Gironcoli, Phys. Rev. B 57, 10421 (1998). 13 G. Kern, G. Kresse, and J. Hafner, Phys. Rev. B 59, 8551 (1999). 14 L. Vocadlo and D. Alfe, Phys. Rev. B 65, 214105 (2002). 15 K. Parlinski and Y. Kawazoe, Eur. Phys. J. B 16, 49 (2000). 16 J. Łażewski, P. Piekarz, and K. Parlinski (unpublished).

ACKNOWLEDGMENTS

Parlinski, Y. Kawazoe, and Y. Waseda, J. Chem. Phys. 114, 2395 (2001). 18 H.-G. Unruh, D. Mühlenberg, and Ch. Hahn, Z. Phys. B: Condens. Matter 86, 133 (1992). 19 J. A. Válgoma, J. M. Perez-Mato, and A. Garcia, Ferroelectrics 237, 73 (2000). 20 J. A. Válgoma, J. M. Perez-Mato, A. Garcia, K. Schwarz, and P. Blaha, Phys. Rev. B 65, 134104 (2002). 21 K. J. Kingma, R. E. Cohen, R. J. Hemley, and H.-K. Mao, Nature (London) 374, 243 (1995); J. Haines and J. M. Leger, Phys. Rev. B 55, 11 144 (1997). 22 H. Hellwig, A. F. Goncharov, E. Gregoryanz, H. Mao, and R. J. Hemley, Phys. Rev. B 67, 174110 (2003). 23 B. B. Karki, R. M. Wentzcovitch, S. de Gironcoli, and S. Baroni, Phys. Rev. B 62, 14 750 (2000). 24 K. Parlinski and Y. Kawazoe, Eur. Phys. J. B 13, 679 (2000). 25 G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993); 49, 14 251 (1994); G. Kresse and J. Furthmüller, Software VASP, IMP-UW, Vienna, Austria, 1999; 54, 11 169 (1996); Comput. Mater. Sci. 6, 15 (1996) 26 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990); G. Kresse and J. Hafner, J. Phys.: Condens. Matter 6, 8245 (1994). 27 J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992). 28 H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). 29 K. Parlinski, Z. Q. Li, and Y. Kawazoe, Phys. Rev. Lett. 78, 4063

104109-6

PHYSICAL REVIEW B 70, 104109 (2004)

QUASIHARMONIC APPROACH TO A SECOND-ORDER… (1997). Parlinski, software PHONON, Cracow, Poland, 2003. 31 K. Parlinski, J. Łażewski, and Y. Kawazoe, J. Phys. Chem. Solids 61, 87 (2000). 32 W. Zhong, R. D. King-Smith, and D. Vanderbilt, Phys. Rev. Lett. 72, 3618 (1994). 33 R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 49, 5828 (1994). 30 K.

34 E.

Artacho, J. Gale, A. Garcia, J. Junquera, P. Ordejon, D. Sanchez-Portal, and J. Soler, Software SIESTA, Spain, 2001. 35 J. Łażewski, K. Parlinski, W. Szuszkiewicz, and B. Hennion, Phys. Rev. B 67, 094305 (2003). 36 Within the quasiharmonic approximation, the difference between the LDA and GGA volumes seems to increase with temperature. See Ref. 6.

104109-7