DFT study on elastic and piezoelectric properties of tetragonal BaTiO3

Computational Materials Science 49 (2010) S372–S377

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

DFT study on elastic and piezoelectric properties of tetragonal BaTiO3 Xiangying Meng a,b, Xiaohong Wen b, Gaowu Qin a,* a b

Key Laboratory for Anisotropy and Texture of Materials (MOE), Northeastern University, Shenyang 110004, PR China College of Sciences, Northeastern University, Shenyang 110004, PR China

a r t i c l e

i n f o

Article history: Received 21 October 2009 Received in revised form 14 April 2010 Accepted 15 April 2010 Available online 15 May 2010 Keywords: First-principles calculations Elastic constants Piezoelectric effects BaTiO3 crystals

a b s t r a c t A theoretical method for studying the mechanism of the piezoelectric effects in crystals from first-principles calculations is described and applied for tetragonal BaTiO3 in the present work. The elastic constants of tetragonal BaTiO3 are calculated by using the homogeneous deformation method in the scheme of density-functional theory and the generalized gradient approximation. The bulk modulus and the Debye temperature are also calculated from the elastic constants and found to be in good agreement with corresponding experimental values. Ó 2010 Elsevier B.V. All rights reserved.

1. Introduction Piezoelectric materials can convert mechanical energy into electrical energy and vice versa. There has been considerable interest in this electromechanical property arising from the coupling of spontaneous polarization with lattice strain, and many devices such as ultrasonic transducers and piezoelectric actuators are designed making use of this property of piezoelectric materials. Now, although experimental techniques are available for the measurement of piezoelectric effect, such as low-frequency equivalent circuit method and Brillouin scattering, material designers may still be eager to predict the piezoelectric effect of a potential, or a designed piezoelectric crystals without any experimental input due to the difficulties in preparing suitable specimens for these materials. Piezoelectricity is a macroscopic property of many materials that display a macroscopic electric polarization upon application of a stress. The calculation of polarization and piezoelectric response within the framework of first-principles methods has proven to be a rather subtle problem. In a landmark paper, Martin [1] showed that the piezoelectric tensor is well defined as a bulk quantity in a crystalline insulator. Following the formalism developed by Martin, Press and Ellis [2] investigated piezoelectricity in copper halides using a localized basis set and a molecular-cluster-embedding method. The quality of their calculation is difficult to assess because of its lacking in the value of dielectric constants. McKitterick [3] dealt with GaAs within the plane-wave pseudopo* Corresponding author. Tel.: +86 24 83683772. E-mail address: [email protected] (G. Qin). 0927-0256/$ - see front matter Ó 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2010.04.026

tential scheme, and predicted a value for the piezoelectric constant that is about one-half of the experimental one. Then Resta et.al [4] showed that by using the linear-response theory (LRT) it is possible to estimate GaAs piezoelectric tensors quite accurately. This excellent agreement is strong dependent on the internal strain parameter n. Nielsen and Martin [5] have shown that the disagreement between theoretical value of the internal strain parameter and experiments can be up to 27%. So this apparently excellent agreement may vary with each individual because the theoretical value of the internal strain parameter is not very certain. In 1993, the geometrical Berry-phase method [6,7] for computing the electric polarization was proposed by King-Smith and Vanderbilt. Later Resta gave a review on this theory [8]. Within the Berry-phase approach the piezoelectric coefficients are accessible via finite differences and then the full piezoelectric stress tensors of PbTiO3 were determined based on this model using the all-electron linear augmented plane wave method [9]. The results are found to be in relatively good agreement with experimental data measured on single crystal material. However, the authors also mentioned that there is a significant spread between various measurements and some subtleties remain regarding the computation of the piezoelectric tensor components by finite differences [10,11]. Vanderbilt [12] further clarified this divergence in that the variations of the Barry-phase with deformation determined the ‘‘proper” piezoelectric response, instead of ‘‘improper” piezoelectric response. Nevertheless, the exact theoretical prediction of piezoelectric constants is still a hard work due to the experimental systematic errors and theoretical itself drawbacks. Here we present an efficient scheme for studying the mechanism of the piezoelectric effects in crystals from first-principles

S373

X. Meng et al. / Computational Materials Science 49 (2010) S372–S377

calculations. The key point of this method for the piezoelectric crystals lies in the fact that the main contribution to the piezoelectric effects is the large relative displacement of cationic and anionic sublattices caused by the macroscopic external strain. This displacement changes the distribution of valence charges in the crystal through relaxations of the positions of the atoms and rehybridization of the electronic orbitals, which then induces a change of polarization in crystals. BaTiO3, with its distinct properties of spontaneous polarization, is an excellent ferroelectric piezoelectric material. There have been considerable works involving both experimental and theoretical methods on BaTiO3 compound. In view of symmetry relations, elastic constants and piezoelectric constants (in Voigt notations) for tetragonal BaTiO3 are:

2

c11 6c 6 12 6 6 c13 C¼6 6 0 6 6 4 0

c12 c11

c13 c13

0 0

0 0

c13

c33

0

0

0

0

c44

0

0

0

0

c44

0

0

0

0

0

3 0 0 7 7 7 0 7 7 0 7 7 7 0 5

2. Calculations Piezoelectricity is the ability of some materials to generate an electric potential in response to applied mechanical stress. Therefore the tensor of the elastic constants at constant electric field must be first known to describe the response of the crystal to external forces. Then, piezoelectricity is associated with the lattice deformations and basically depends on the polarization in crystals. Our method for calculating elastic constants draws upon finite strain continuum elasticity theory [15,16]. Isothermal elastic constants are defined by expanding the Helmholtz free energy per unit mass as a Taylor series in elastic strain at constant temperature. We perform a first-principles calculations at zero temperature, so the internal energy E, instead of the Helmholtz free energy F, is used. To calculate elastic constants using the method of homogeneous deformation, a specified Lagrangian strains gij is applied to the equilibrium lattice as it removes irrelevant crystal rotations, then the energy E (gij) for the strained crystal is calculated based on density-functional theory (DFT). The elastic constants are extracted as proportional to the second order coefficient in a polynomial fit of the total energy as a function of Lagrangian strains gij [15–17]:

c66

q0 Eðgij Þ ¼ q0 Eð0Þ þ and

2

0 6 d¼4 0 d31

3

0

0

0

d15

0

0

d15

0

7 0 5:

d31

d33

0

0

0

0

All six independent elastic constants, as well as three independent piezoelectric constants of tetragonal BaTiO3 single crystals have been determined experimentally [13,14]. This motivated us to calculate all six independent elastic constants and related piezoelectric properties of tetragonal BaTiO3 from first principle electronic structure calculations. By comparing theoretical and experimental data we can judge the applicability of the presented method in searching for other piezoelectric crystals. In the present work we have no intention of giving the theoretical prediction of piezoelectric constants for crystals although the dipole moment of the cell can indeed be calculated with the method presented. There are two main reasons. First, experimental piezoelectric values are often lacking or affected by systematic errors, especially for a potential or designed piezoelectric crystal. In the case of tetragonal BaTiO3 crystal which has been widely investigated, different experimental values of d15 [13,14], 28.2 and 58 (in the unit of 1011 C/N), have been reported. Consequently the theoretical predictions of piezoelectric constants can make only a weak power due to the uncertain experimentations. Second, from the viewpoint of application, it is not necessary for material designer to know the exact piezoelectric coefficients in the beginning of searching a novel material. In most cases a prediction of good piezoelectric phenomenon is enough. In this paper we discuss the influence of the lattice distortions on polarization properties for tetragonal BaTiO3 in detail. As fundamental parameters, elastic constants, bulk modulus and Debye temperature of tetragonal BaTiO3 are also calculated from firstprinciples calculations. Furthermore, the mechanism of the piezoelectric effect has been basically understood through the calculations and the differences in piezoelectric behavior along the different crystallographic directions have also been explained. The numerical accuracy and predictive power of presented scheme allow us to make quantitative assertions about the improvement in the piezoelectric properties of novel materials.

  1X C ijkl gij gkl þ O g3ij ; 2 ijkl

ð1Þ

where q0 is the initial mass density of the material and the initial state is assumed to be stress-free. The expansion coefficients Cijkl is elastic constants defined as [15–17]:

C ijkl ¼ q0

@2E j : @ gij @ gkl g¼0

ð2Þ

We make use of the Voigt notation (replaces 11, 22, 33, 23, 13 and 12 by 1, 2, 3, 4, 5, and 6, respectively) for the tensor indices to write Cijkl as Cab. To account for the symmetry of deformation, we introduced the factor ni, which takes the value 1 if the Voigt index is 1, 2, or 3 and the value 2 if the Voigt number is 4, 5, or 6. Thus, the Taylor expansion of the total energy under Eulerian strain, d, can be done in powers of the strain tensor with respect to the initial energy of the unstrained crystal in the following way [18]:

EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0

X i

si ni di þ

1X cij ni di nj dj 2 ij

! þ Oðd3 Þ:

ð3Þ

In the above equation, V0 is the volume of the unstrained system, si is an element in the stress tensor. The Taylor expansion of the total energy should be done in terms of Lagrangian coordinates g, while we use the Eulerian coordinates d here. Since we shall consider only small lattice distortions in order to remain within the elastic limit of the crystal, the approximation g = d is valid in the Taylor expansion. However, the homogeneous deformation lowered the symmetry of the crystal in the calculations of total energy for the strained crystals. To obtain the strained unit cell, the symmetric distortion matrices D was applied to the lattice vectors R according to the rule R0 = RD. The R0 is a matrix containing the components of the distorted lattice vectors. The symmetric distortion matrices D was determined from the Lagrangian strains gij by inverting the equation [19]:

gij ¼

1X ðDki Dkj  dij Þ: 2 k

ð4Þ

For tetragonal BaTiO3 there are six independent elastic constants, c11, c33, c44, c66, c12 and c13. We first apply the distortion matrices

S374

X. Meng et al. / Computational Materials Science 49 (2010) S372–S377

0 B D11 ¼ @

1

1þd 0 0

0 C 1 0A

0

0 1

and

0

D33

1 1 0 0 B C ¼ @0 1 0 A; 0 0 1þd

to the crystal, and correspond to straining the lattice along x, and zaxis, respectively. In these distortions, the symmetry of the strained lattice still remains tetragonal. However, the volume is slightly changed by the small distortion. The energy associated with these distortions can be obtained by putting the values of the matrices D11 and D33 in Eq. (3) as [18]:

EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0



s1 d þ

c11 2  d 2

ð5Þ

s3 d þ

c33 2  d ; 2

ð6Þ

and

EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0



respectively. From above relations we obtain the elastic constants c11 and c33. Similarity, using the following volume-conserving monoclinic shear distortions,

2 6 6 D44 ¼ 6 4 and

2

1 ð1d2 Þ1=3

0

0

0

1 ð1d2 Þ1=3

d ð1d2 Þ1=3

0

d ð1d2 Þ1=3

1 ð1d2 Þ1=3

1

d ð1d2 Þ1=3

0

1 ð1d2 Þ1=3

0

0

1 ð1d2 Þ1=3

2 1=3

D66

6 ð1d Þ 6 d ¼ 6 ð1d2 Þ1=3 4 0

3 7 7 7 5 3 7 7 7; 5

P¼

we obtain c44 and c66 from

EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0 ð2s4 d þ 2c44 d2 Þ

ð7Þ

and 2

EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0 ð2s6 d þ 2c66 d Þ;

ð8Þ

respectively. Also, applying the volume-conserving orthorhombic distortions

2 D12

and

6 ¼6 4

2

6 6 D13 ¼ 6 4

1þd ð1d2 Þ1=3

0

0

3

0

1d ð1d2 Þ1=3

0

0

0

1 ð1d2 Þ1=3

7 7 5

1þd ð1d2 Þ1=3

0

0

0

1 ð1d2 Þ1=3

0

0

0

1d ð1d2 Þ1=3

Next, the piezoelectric effect is determined by calculating the dependence of macroscopic polarization on the local strains. After a fully optimization of the ions system, we calculate the changes of the dipole moments along the different crystallographic directions in the subsequent strained lattices. The analysis of polarization upon strains allows us to explain the piezoelectric behavior of tetragonal BaTiO3 single crystals. The geometry optimization and the total energy calculations have been carried out using the ab initio total-energy and molecular dynamics program VASP (Vienna ab initio simulation package) based on the density-functional theory (DFT) [20,21]. The current calculations make use of the VASP implementation of ultrasoft pseudopotentials [22] and an expansion of the electronic wave functions in plane waves with a kinetic-energy cutoff of 600 eV. All calculated results were derived employing the generalized gradient approximation (GGA) for exchange and correlation due to Perdew and Wang [23]. Brillouin-zone integrations were performed using Monkhorst–Pack k-point meshes [24] and the Methfessel–Paxton technique [25] with a 0.1 eV searing of the electron levels. For each structural optimization, tests were carried out using different k-point meshes to ensure absolute convergence of the total energy with respect to the structural degrees of freedom to a precision of better than 2 meV/atom. The total energy is minimized with respect to the volume, the unit cell-external degree of freedom (i.e. the unit shape) and unit cell-internal degree of freedom (i.e. the Wyckoff positions of all atoms). The data sets of DE  d were obtained by varying d from the equilibrium lattice in steps of 0.002. The change of dipole moments and polarization property are analyzed using the software utility ‘‘LEV00” [26]. We still define the polarization P of a unit cell in a classical way:

3 7 7 7; 5

to the lattice, c12 and c13 can be calculated from

  1 EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0 ðs1  s2 Þd þ ðc11 þ c22  2c12 Þd2 2

ð9Þ

1

"

v

e

X

Z l Rl þ

Z

# rqðrÞ dr ;

ð11Þ

l

where e is the electron charge, v is the sample volume, the l summation is over the ionic sites, eZl are bare ionic charges. q(r) is the electronic charge density and can be calculated by

qðrÞ ¼

Occ X

jwnk ðrÞj2 :

ð12Þ

n;k2BZ

Although Vanderbilt [7] and Resta [8] have defined polarization as a quantum phenomenon, the dipole moment of the cell can indeed be calculated with the method presented as the choice of the integration region is not important for calculating the cell dipole moment difference. Here it is possible to integrate the charge density over one or several spheres of different radii. The point is that the notion of atomic charge is not well defined. To analyze charge transfer during atom adsorption, we calculated the charge of the atom in spheres of several different radii and compare them. In principle, calculating the dipole moment of the density in the cell is only meaningful in the cases of molecules contained inside big cells. We choose a sphere which has the whole unit cell inside it with the center somewhere nearby the cell center. The smallest and largest radii both included the whole BaTiO3 molecule. Since the charge is not properly normalized, a ‘‘non-conserving” algorithm is used for the charge integration. After that, a reasonable estimate of dipole moments can be obtained.

and

1 EðV; dÞ ¼ EðV 0 ; 0Þ þ V 0 ½ðs1  s3 Þd þ ðc11 þ c33  2c13 Þd2 ; 2

ð10Þ

respectively. The already calculated elastic constant c11, which is equal to c22, and c33 are used in the calculations of c12 and c13. Therefore, all six independent elastic constants of tetragonal BaTiO3 are determined.

3. Results and discussion 3.1. Elastic constants, bulk modulus and Debye temperature The tetragonal BaTiO3 is well known as a simple ABO3-type perovskite crystal with P4mm symmetry [27]. Its primitive cell

X. Meng et al. / Computational Materials Science 49 (2010) S372–S377 Table 1 The equilibrium structural parameters for BaTiO3 obtained from GGA calculations. The lattice parameters a, and c are in Å, and the equilibrium unit cell volume (V0) is in Å3.

GGA Exp.a a

From Ref. [27].

a

c

c/a

V0

3.984 4.000

4.066 4.024

1.021 1.006

64.537 64.384

S375

has one BaTiO3 molecule and the lattice parameters are a = b = 4.000 Å and c = 4.024 Å. In the tetragonal phase Ba ions occupy the eight apexes of the cube, Ti ions are located at the position 0.201 Å above the body center of the octahedron along the C4-axis and O ions are located in the face center of the cube to form an octahedron (TiO6). Before performing calculations to obtain the elastic constants we have first optimized the structural parameters of BaTiO3 using GGA. We first adopted the experimental a, and c and optimized the

Fig. 1. Changes in the energy (DE) under various strains: (a) D11, (b) D33, (c) D44, (d) D66, (e) D12, and (f) D13. The black square represent the calculated values and the line is the polynomial fit.

S376

X. Meng et al. / Computational Materials Science 49 (2010) S372–S377

equilibrium volume. We then used the theoretical equilibrium volume and lattice parameters to calculate elastic constants. In Table 1 we list the GGA and experimental values of the equilibrium lattice parameters and volume. The theoretically estimated equilibrium volume is overestimated by 0.3% in our GGA calculations. We have also found that the equilibrium lattice parameters obtained from GGA are in good agreement with the experiments. In Fig. 1 we have plotted the calculated changes in total energy (DE) versus strain d for tetragonal BaTiO3 for the six types of distortions discussed above. The elastic constants cij are then obtained by means of polynomial fits truncated after the third-order term. We used very small distortion to minimize the errors coming from higher order terms. The elastic constants obtained from our GGA calculations are given in Table 2. We can found that our theoretical results provide better agreement with experimental data than results from previous phenomenological calculations. The bulk modulus B, which represents the resistance and influences the speed of sound and other mechanical waves in the material, is an important material parameter. From elastic theory, a formula for calculating the bulk modulus and its components along different crystallographic axes is provided with the following result [18]:

K ; ð1 þ a þ bÞ2 K Ba ¼ ; ð1 þ a þ bÞ Ba Bb ¼ ;

B¼

ð13Þ ð14Þ ð15Þ

a

and

Bc ¼

Ba ; b

ð16Þ

where K = c11 + 2c12a + c22a2 + 2c13b + c33b2 + 2c23ab. For tetragonal crystals a = 1 and b is defined as

b¼

c11 þ c12  2c13 : c33  c13

ð17Þ

The Debye temperature of a material is a suitable parameter to describe phenomena of solid-state physics which are associated with lattice vibrations. It basically depends on the elastic constants. H. Siethoff presents a Debye-temperature-elastic-constants relationship for tetragonal symmetry crystals in the following way [18]:

hD ¼

ahC t Gt 1=6

kMs

als belonging to the tetragonal crystal system, Gt is associated with the elastic constants and defined as

Gt ¼

 1=3 c44 c66 ðc11  c12 Þ : 2

ð19Þ

Therefore the present calculated Debye temperature hD = 503.5 K. All the calculated material parameters are listed in Table 3 and found to be in good agreement with the experimental values. 3.2. Piezoelectric effects To investigate the effect of lattice distortions on the changes of polarization in tetragonal BaTiO3 crystals, the dipole moments have been calculated from first-principles for each strained lattice. Since we have three independent piezoelectric constants, in Voigt notations, d31, d33, and d15. The meaning of the notation, for example d15, which is equal to d24 in view of the symmetry, can be explained as the change of the polarization along y-axis (DPy) with respect to the lattice distortion D44. Similar explain can be made for other two piezoelectric constants, d31 and d33 as well. So we apply the distortion, D11, D33, and D44 to the BaTiO3 lattice to calculate the corresponding changes of the dipole moments along y-axis and z-axis, respectively. The changes of polarization along different directions versus strain d are plotted in Fig. 2. Our result shows that change of polarization is linear over a wide rage of strains. We can see from Fig. 2 that the changes of polarization DP are almost in direct proportion to the lattice distortions in small strain d. As the piezoelectric effect is determined by the dependence of macroscopic polarization on the local strains, the rates DP/Dd can present the values of piezoelectric effect in different directions in tetragonal BaTiO3. Our calculation results show that for the same Table 3 Some material parameters of tetragonal BaTiO3 crystals: the piezoelectric constants (dij in 1011 C/N), the bulk modulus and its components (B in 1011 N/m2), and the Debye temperature (hD in K).

GGA Exp. a b c

d31

d33

d15

B

Ba

Bb

Bc

hD

2.325 3.34a

7.072 9.00a

25.14 28.2a

1.354 1.41b

6.394 –

6.394 –

2.350 –

503.5 480.16c

From Ref. [13]. From Ref. [31]. From Ref. [28].

ð18Þ

;

where h is Plank constants and k is Boltzmann constants. s is the number of atoms and M is the weighted arithmetical average of the masses of the species in the crystallographic unit cell. Ct is a constant, a fit to the data sets of all analyzed tetragonal materials yielded Ct = 5.05  1011 [28]. a is the lattice parameter. For materiTable 2 The GGA results of elastic constants (cij in 1011 N/m2) for tetragonal BaTiO3 crystals.

a b c

cij

GGA

Phenomenological methoda

Exp.b

Exp.c

c11 c33 c44 c66 c12 c13

2.548 1.585 0.680 1.186 1.014 1.041

1.67 1.26 0.52 0.93 1.25 1.21

2.22 1.51 0.61 1.34 1.08 1.11

2.11 1.60 0.56 1.27 1.07 1.14

From Ref. [30]. From Ref. [13]. From Ref. [32].

Fig. 2. The changes of the dipole moment in different directions with respect to the lattice distortions.

X. Meng et al. / Computational Materials Science 49 (2010) S372–S377

deformation d the change of polarization along y-axis caused by strain D44 is largest and along z-axis caused by strain D33 is next, while the change of polarization along z-axis caused by strain D11 is smallest. The result hints that the piezoelectric constants d15 is the largest in magnitude and d31 is the smallest. To test our numerical accuracy and predictive power of presented scheme, we calculated the piezoelectric constants of tetragonal BaTiO3 crystals. In describing the response of a piezoelectric crystal to external mechanical fields one can use materials constants depending on deferent sets of independent variables. Here we use the elastic constants at constants electric field strength cijkl and the piezoelectric stress tensor eijk. The materials constants in such a set can be related as [29]:

eijk ¼ clmjk dilm :

ð20Þ

The tensor dilm is the piezoelectric constants. By solving Eq. (20), we can get the explicitly analytical solutions for the piezoelectric constants dijk in terms of the piezoelectric stress tensor eijk and the elastic constants cijkl. In Voigt notations, they are [13]:

c33 e31  c13 e33 ; ðc11 þ c12 Þc33  2c213 ðc11 þ c12 Þe33  c13 e31 ¼ ; ðc11 þ c12 Þc33  2c213 e15 ¼ : 2c55

d31 ¼ d33 d15

S377

were simulated from first-principles calculations. The main idea of this method lie in the fact that the external forces induced lattice distortions change the distribution of the valence charge in the crystal through relaxations of the positions of the atoms and rehybridization of the electronic orbitals, which then induces a change of polarization in crystals. Our calculations successfully explained the differences in piezoelectric behavior along different crystallographic axes in BaTiO3 and agree well with the experimental results. We hope that our results will be beneficial for the design and search for other good piezoelectric materials in the future. Acknowledgments The authors would like to thank Prof. Chuangtian Chen for helpful discussions. This work was supported by the Chinese Research Fund for the Dectoral Program of Higher Education (No. 200801451021) and the Fundamental Research Fund for the Central Universities (No. N090405016). References

ð21Þ ð22Þ ð23Þ

Previously measured piezoelectric stress tensors are e31 = 0.7, e33 = 6.7, and e15 = 34.2 (in the unit of C m2) [13]. The constants are given at room temperature at the wavelength 633 nm. Both the theoretical and experimental values of the piezoelectric constants are also shown in Table 3. We find that the ordering of DP/Dd is completely in agreement with that of the magnitude of the three independent piezoelectric constants. The same sequence in magnitude of the differential variation DP/Dd and dij can clearly help us judge the piezoelectric behavior along the different directions in tetragonal BaTiO3 crystals before the experimental results. By calculating the changes of the polarization DP under the strain, we can also make quantitative assertions about the improvement in the piezoelectric properties of new materials by comparing it to the familiar piezoelectric crystals. 4. Conclusions The strain energies for six different distortions of tetragonal BaTiO3 are calculated using GGA in the theoretically optimized crystal structure in order to calculate elastic constants. As fundamental parameters, bulk modulus and Debye temperature of tetragonal BaTiO3 are also calculated here from the elastic constants and are compared with experimental values. Furthermore, We discuss the influence of the lattice distortions on polarization properties for tetragonal BaTiO3 in detail. An efficient computational scheme has been designed and the piezoelectric effects in BaTiO3 crystals

[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32]

R.M. Martin, Phys. Rev. B 5 (1972) 1607–1613. M.R. Press, D.E. Ellis, Phys. Rev. B 38 (1988) 3102–3111. J.B. McKitterick, Phys. Rev. B 28 (1983) 7384–7386. S. de Gironcoli, S. Baroni, R. Resta, Phys. Rev. Lett. 62 (1989) 2853–2856. O.H. Nielsen, R.M. Martin, Phys. Rev. Lett. 50 (1983) 697–700. R.D. King-Smith, D. Vanderbilt, Phys. Rev. B 47 (1993) 1651–1654. D. Vanderbilt, R.D. King-Smith, Phys. Rev. B 48 (1993) 4442–4455. R. Resta, Rev. Mod. Phys. 66 (1994) 899–915. G. Sághi-Szabó, R.E. Cohen, H. Krakauer, Phys. Rev. Lett. 80 (1998) 4321–4324. A. Dal Corso, M. Posternak, R. Resta, A. Baldereschi, Phys. Rev. B 50 (1994) 10715–10721. F. Bernardini, V. Fiorentini, D. Vandefbilt, Phys. Rev. B 56 (1997) 10024–10027. D. Vanderbilt, J. Phys. Chem. Solid. 61 (2000) 147–151. } nter, M.H. Garrett, D. M. Zgonik, P. Bernasconi, M. Duelli, R. Schlesser, P. Gu Rytz, Y. Zhu, X. Wu, Phys. Rev. B 50 (1994) 5941–5949. Z. Li, S. -K. Chan, M.H. Grimsditch, E.S. Zouboulis, J. Appl. Phys. 70 (1991) 7327– 7332. R.N. Thurston, Physical Acoustics Principles and Methods, Academic, New York, 1964. D.C. Wallace, Solid State Physics, Academic, New York, 1970. K. Brugger, Phys. Rev. 133 (1964) A1611–A1612. P. Ravindran, Lars Fast, P.A. Korzhavyi, B. Johansson, J. Wills, O. Eriksson, J. Appl. Phys. 84 (1998) 4891–4904. Jijun Zhao, J.M. Winey, Y.M. Gupta, Phys. Rev. B 75 (2007) 094105– 094111. G. Kresse, Phys. Rev. B 54 (1996) 11169–11186. G. Kresse, Comput. Mater. Sci. 6 (1996) 15–50. D. Vanderbilt, Phys. Rev. B 41 (1990) 7892–7895. John P. Perdew, Yue Wang, Phys. Rev. B 45 (1990) 13244–13249. H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188–5192. M. Methfessel, A.T. Paxton, Phys. Rev. B 40 (1989) 3616–3621. L.N. Kantorovich. . S. Aoyagi, Y. Kuroiwa, A. Sawada, I. Yamashita, T. Atake, J. Phys. Soc. Jpn. 71 (2002) 1218–1221. Hans Siethoff, Karl Ahlborn, J. Appl. Phys. 79 (1996) 2968–2974. J.F. Nye, Physical Properties of Crystals, Clarendon, Oxford, 1957. A. Khalal, D. Khatib, B. Jannot, Physica B 271 (1999) 343–347. A. Schaefer, H. Schmitt, A. Dorr, Ferroelectrics 69 (1986) 253–266. Z. Li, M. Grimsditch, C.M. Foster, S.-K. Chan, J. Phys. Chem. Solid. 57 (1996) 1433–1438.