Density Functional Perturbation Theory to Predict Piezoelectric Properties Density Functional Perturbation Theory to Predict Piezoelectric Properties

Among the various computational methods in materials science, only first-principles cal- culation based on the density functional theory has predictability for unknown material. Especially, density functional perturbation theory (DFPT) can effectively calculate the second derivative of the total energy with respect to the atomic displacement. By using DFPT method, we can predict piezoelectric constants, dielectric constants, elastic con - stants, and phonon dispersion relationship of any given crystal structure. Recently, we established the computational technique to decompose piezoelectric constants into each atomic contribution, which enable us to gain deeper insights to understand the piezo electricity of material. Therefore, in this chapter, we will introduce the computational framework to predict piezoelectric properties of polar material by means of DFPT and details of decomposition technique of piezoelectric constants. Then, we will show some case studies to predict and discover new piezoelectric material.


Introduction
In this chapter, we will introduce how recent computational techniques can successfully predict response properties, represented as piezoelectricity, by means of perturbation method. Piezoelectricity is the polarization change in response to external mechanical force. Inversely, if electrical field is applied to piezoelectric material, mechanical strain is induced (inverse piezoelectric effect). Therefore, piezoelectric materials are widely used as vibrational censors, surface acoustic wave devices, and actuators. Only the material having no inversion symmetry shows piezoelectricity. For example, Figure 1 shows schematic illustration of the piezoelectric effect. Positions of positively charged ion (cation) and negatively charged ion (anion) are represented as plus and minus symbols. Figure 1a shows the paraelectric phase, where ions are orderly located with inversion symmetry. On the other hand, ions are slightly displaced by δ with respect to those in paraelectric phase, as shown in Figure 1b. Such small displacement induces microscopic polarization P s along the ionic displaced direction.
Because ferroelectric phase is energetically more stable than paraelectric phase under low temperature, P s is frequently referred as the spontaneous polarization. Above Curie temperature, ferroelectric properties are disappeared since paraelectric phase becomes more stable than ferroelectric one. Figure 1c shows the schematic illustration of the principal of piezoelectricity, where external stress (red-colored arrows) increases the ionic displacement and resultant polarization. In this case, external stress increases the spontaneous polarization by ΔP s = P s ' −P s . Therefore, piezoelectric constant is defined as the derivative of the spontaneous polarization with respect to the external field. More detailed and comprehensive description of piezoelectricity is reviewed by Martin [1].
First-principles calculation based on density functional theory (DFT [2,3]) has been widely utilized as the computational method to predict the electronic properties of material under the ground state. Ideally, required information to conduct the first-principles calculation is only the crystal structure, including atomic species and position of periodic/nonperiodic structure unit. The most significant advantage of first-principles calculation is its predictability. Since King-Smith and Vanderbilt showed the theoretical methodology to calculate change in polarization per unit volume ΔP [4], dielectric and piezoelectric properties of wide range of materials in which electronic correlations are not too strong [5][6][7] have been accurately predicted. The derivative of total energy determines various properties. For example, determined forces, stresses, dipole moment (first-order derivatives), dynamical matrix, elastic constants, dielectric and piezoelectric constants (second-order derivative), nonlinear dielectric susceptibility, phonon-phonon interaction and Grüneisen parameters (third-order derivative), and so on. Thus, first-principles calculation has been made use of calculating the perturbed total energy of materials because of its accuracy. Although perturbations were made by hand up to the early 1980s, sophisticated methodology of density functional perturbation theory (DFPT) was proposed in 1987 by Baroni et al. [5]. They showed general formulation of total energy change with respect to atomic displacement and opened the way to efficiently compute the energy derivative with respect to the perturbation [5]. DFPT can compute response properties directly arising from the perturbations of strain, atomic displacement, and electric field by making use of linear response theory [8][9][10][11]. Numbers of ferroelectric materials are theoretically investigated on the origin of their ferroelectric properties (including piezoelectricity and dielectric properties) by using DFPT. Because of technological importance, such theoretical researches have been focused on Pb-based perovskite material (e.g., PbTiO 3 , PbZrO 3 , and their solid solution [12][13][14][15]) because they have excellent piezoelectric properties and are widely applied for actuators. However, due to the restriction of hazardous substance (RoHS) directive, researches on lead-free ferroelectric materials have gathered great attraction. By taking advantage of the predictability of DFPT, various lead-free ferroelectric oxide and nitride materials were theoretically investigated on their piezoelectric properties [16][17][18][19][20][21][22][23][24][25]. Moreover, DFPT calculations showed that piezoelectricity can be greatly enhanced by imposing isotropic stress for PbTiO 3 [26,27], uniaxial stress for SrHfO 3 [28], uniaxial and biaxial strain for AlN-GaN solid solution alloy [29], and two-dimensional epitaxial strain for doped ZnO [30]. As latterly explained, those enhancements of piezoelectric constant are thought to be closely related to the phase transition. In the next section, we will show the definition of piezoelectric constants within the framework of DFPT.

Formulation of piezoelectric constants
Formulation and calculation methodologies to obtain response properties of materials in the framework of DFPT have been developed in a step-by-step manner, because degrees of freedom by perturbations of atomic displacement, homogeneous electric fields, and strain are often strongly coupled. For example, piezoelectricity affects elastic and dielectric properties. Therefore, special care must be paid for the calculation of coupled properties. In 2005, Hamann et al. demonstrated that elastic and piezoelectric tensors can be efficiently calculated by treating homogeneous strain within the framework of DFPT [31]. At the same time, Wu et al. systematically formulated response properties with respect to displacement, strain, and electric fields [32]. In this section, we will briefly introduce how piezoelectric properties are formulated in the framework of DFPT. In each formulation, Einstein implied-sum notation is used. Cartesian directions {x, y, z} are represented as α and β. Subscription of j and k = 1, …, 6 is the standard Voigt notation (represents directions of xx, yy, zz, yz, zx, and xy). The subscripts m and n are the degrees of freedom in the cell. They range from 1 to 3i, where i is the number of irreducible atoms because each atom has three degree of freedom along x, y, and z directions.
Total energy of material under perturbation of atomic displacement u, electric field σ, and strain η, E(u,σ,η), is defined as follows: where E 0 is the total energy of material under the ground state, Ω 0 is volume of the unit cell (smallest repeat unit of crystal), Ω is deformed volume of the unit cell, and P is the electric polarization. Following response functional tensor can be obtained by second-order differential of Eq. (1): Clamped − ion term of electric susceptibility : Clamped − ion term of elastic tensor : Born effective charge tensor : Clamped − ion piezoelectric tensor : Clamped-ion term is a frozen quantity, which indicates that atomic coordinates are not allowed to relax as the homogeneous electric field or strain. Therefore, dynamical term should be added into the clamped-ion term in order to obtain proper response properties.
Simplest and physically well-understandable piezoelectric constant can be expressed as follows: In this expression, it is easily understood that piezoelectric e constant e αj is a measure of the change in polarization induced by the external strain. As the atomic positions are changed according to the strain, change of the polarization includes both electronic contribution (clamped-ion term) and dynamical contribution (internal-strain term). The internal-strain term of piezoelectric constant is represented as follows: Thus, proper piezoelectric constant can be obtained by Eqs. (7) and (9): Here, the first and second terms on the right-hand side in Eq. (10) are the clamped-ion term and internal-strain term, respectively. The former shows the electronic contribution ignoring the atomic relaxation effect, and the latter shows the ionic contribution including the response of the atomic displacement to the strain. The Born effective charge Z mα , force-constant matrix K mn , and internal-strain tensor Λ nj are the second derivatives of the energy with respect to the displacement and electric field, pairs of displacements, and displacement and strain, respectively. The internal-strain term of the piezoelectric stress constants can be further decomposed into the individual atomic contributions when the above second-derivative tensors are fully obtained.
On the other hand, the internal-strain term of the piezoelectric stress constant e αj is frequently described by the following equation, using the Born effective charge Z αβ and displacement u β of each atom in the calculation cell: shows the response of the first-order atomic displacement to the first-order strain. In this expression, the meaning of the piezoelectric stress constant, i.e., e j is a measure of the change in polarization induced by the external strain, is much more visible than in Eq. (9). In the DFPT formalism, ∂ u β / ∂ η j is implicitly calculated as a displacement-response internalstrain tensor Γ as follows [32]: Because the subscript n in Γ nj indicates the degrees of freedom, Γ nj can be decomposed into the individual atomic components, which also enables to calculate individual contribution of each atom for total piezoelectric constant.
Here, piezoelectric e constant defined as Eq. (9) is frequently referred as "piezoelectric strain constant." On the other hand, it is much more natural and easy to control the stress (electric field) than to control the strain in any case. In this case, the piezoelectric strain constant d αj is usually measured. It can be obtained from piezoelectric strain constant e αj using the following relation: where s jk is the elastic compliance, which is given by the inverse matrix of the elastic constants C jk .
Those formulations are implemented in specific first-principles simulation packages such as ABINIT [33] and Vienna ab initio simulation package (VASP) [34], and piezoelectric constants can be calculated on a daily basis. From the next section, we will show how DFPT calculation precisely gives piezoelectric properties of ferroelectric materials.

Introduction of target material
In this chapter, we selected LiNbO 3 as a target material to show the predictability of DFPT calculation. LiNbO 3 is one of ferroelectric materials and widely used as surface acoustic wave (SAW) and optical waveguide elements. Crystal structure of LiNbO 3 , which belongs to the space group of R3c, is frequently referred as "strained perovskite structure." Schematic illustrations of crystal structure of ABO 3 perovskite and LiNbO 3 are shown in Figure 2.
Crystal structures shown in the present chapter was visualized by using VESTA software [35]. Curie temperature of LiNbO 3 is quite high and ranges from 1140 [36] to 1210°C dependent on the quality of sample (variation of Li/Nb relation can shift Curie temperature [37]). Below Curie temperature, ferroelectric phase with R3c symmetry (crystal structure can be classified into 230 types of space group according to the symmetry group) shown in Figure 2a is stable. Paraelectric phase with R ¯ 3 c symmetry, shown in Figure 2b, becomes stable above Curie temperature. In the paraelectric phase, it can be seen that both Li and oxygen is positioned with the same height along c-axis, and the position of Nb is just the center between two oxygen layers along c-axis. On the other hand, both Li and Nb are shifted in ferroelectric R3c phase along downward direction of c-axis with respect to those in paraelectric R ¯ 3 c phase.
Due to the different bonding nature between Li-O and Nb-O, atomic positions of Li and Nb are off-centered within oxygen layers along c-axis. This structural characteristic is the ferroelectric nature of LiNbO 3 . One of the notable properties of LiNbO 3 is its high-curie temperature (~1400 K). However, piezoelectric properties of LiNbO 3 are not so much superior as compared with Pb-based perovskites. Crystal structure of piezoelectric ABO 3 perovskite is based on the cubic structure (of Pm3m symmetry), shown in Figure 3a.
Cubic lattice is symmetric and usually high-temperature phase, same as LiNbO 3 . The "strained perovskite structure" expression for LiNbO 3 means that LiO 6 and NbO 6 polyhedron are largely rotated with respect to the cubic perovskite structure. However, because of the simple atomic configuration of cubic structure, atoms can be displaced along various directions and change crystalline symmetry as shown in Figure 3a. Crystalline lattice is vibrated (referred as phonon) under finite temperature. Some lattice vibrations along specific directions are unstable. This specific phonon is called as soft mode with imaginary frequency. In such case, atoms are displaced along unstable phonon mode to lower the total energy. For example, cooperative atomic displacement along [001] direction shown in Figure 3b (referred as Γ 15 mode) changes symmetry from cubic to tetragonal (of P4mm symmetry), which leads polarization along [001] direction. Thus, polarization direction of perovskite is not restricted and allowed to be changed. This characteristic rotational polarization direction is favorable for piezoelectricity because grains in polycrystalline material are oriented along various directions. Thus, careful controlling of crystal structure is essential to obtain superior piezoelectric properties.
The most convenient way to control and drastically change the crystal structure is imposing high pressure. Many compounds have found to be possible to form LiNbO 3 -type structure under the high-pressure synthesis [38], and some of them were quenchable phase. For example, LiNbO 3 - Density Functional Perturbation Theory to Predict Piezoelectric Properties http://dx.doi.org/10.5772/intechopen.76827 type structured ZnSbO 3 was successfully synthesized [39] under high pressure, and improvement of the spontaneous polarization is suggested by enhancement of the covalency of Sn site from first-principles simulation [40]. Moreover, high-pressure synthesized research on LiNbO 3type structure is now extended to more complex compounds such as oxynitrides [41,42].
The crystal structure of ABO 3 compound is generally determined by the balance between the ionic radius of A and B element, which is frequently referred as tolerance factor. Due to the small size of the Li ion with respect to the tolerance factor of LiNbO 3 , this compound cannot form stably the popular perovskite structure under the ambient condition. On the other hand, we predicted the crystal structures of high-pressure phase of LiNbO 3 [43], which were not completely elucidated by experimental study [44]. Revealed structures are NaIO 3 -type structure (Pnma) as room temperature high-pressure phase and apatite-like structure (P6 3 /m) as hightemperature high-pressure phase. It should be noted that the NaIO 3 -type structure is closely related with the popular GdFeO 3 -type perovskite structure. The only difference between these structures is that A-site position and B-site position are inter-exchanged. Therefore, there seems to be a possible way to connect the perovskite structure and LiNbO 3 -type structure.
In our previous theoretical study on high-pressure phase, analysis was mainly concerned with phase transition mechanism only from the viewpoint of subgroup symmetry and energy barrier [43]. It will be instructive to deal with this phase transition phenomenon from the viewpoint of lattice instability as discussed in the field of the ferroelectric instability analysis.
In the following section, we will show investigation on the potential piezoelectric properties of LiNbO 3 with various hypothetical crystal structures by the method of DFPT, and possible phase transition mechanism will be discussed from the viewpoint of soft mode.

Computational methodology
First-principles calculation was performed by using VASP code [34]. Interactions between ion and electron were treated by projector augmented wave (PAW) method [45]. PBEsol functional [46] was used to approximate exchanges and correlate interactions of electrons, which can be used to reproduce the lattice constants of various materials [45]. Precise calculation on the lattice constant is essential to predict piezoelectric properties because they depend on volume of unit cell Ω as shown in Eq. (1). The kinetic energy cutoff for plane waves was set at 500 eV, and the k-point mesh was set at ~0.03/Å intervals to obtain the converged total energy at less than 0.1 meV/atom. Before calculating the piezoelectric constants, atomic positions and cell parameters were optimized until the forces on each atom and cell converged at below 5 × 10 −4 eV/Å.
Since VASP does not directly calculate Eq. (12), we added routine to calculate displacementresponse internal-strain tensor Γ nj and decompose piezoelectric constants into each atom. The sum of the decomposed piezoelectric constants was confirmed to accurately reproduce the total piezoelectric constants. Careful convergence tests with a higher energy cutoff and denser k-point mesh showed that the numerical accuracy of the calculated Γ nj was less than 0.01. It was confirmed that this error does not influence our discussion and conclusion. Moreover, it was confirmed that the values of Γ nj obtained by the DFPT method were consistent with those calculated by the direct method, in which the strain-displacement relation of each ion was explicitly calculated.
On the basis of cubic Pm3m phase, lattice instability analysis was performed by phonon calculation utilizing phonopy code [47]. Force constant matrix shown in Eq. (2) was constructed by DFPT calculation implemented in VASP code combined with supercell approach. Supercell was constructed by using unit cell so that orthogonal three axes of the supercell exceed 10 Å. Note that although supercell is not required in DFPT approach, the present VASP code implements perturbation at the zone center.

Calculated piezoelectric properties of LiNbO 3
Calculated piezoelectric properties of LiNbO 3 in ferroelectric phase are summarized in Table 1. Some experimentally measured values are also shown in Thus, Li vacancy is considered to have negligible influence on the piezoelectric properties.
Decomposed ionic contribution of piezoelectric strain constant e 33 is summarized in Table 2.
Although the Born effective charge of Nb is larger than its formal charge +5e, displacementresponse internal-strain constant of Nb is negative value. This indicates that piezoelectricity of LiNbO 3 is mainly dominated by displacement of Li. Born effective charge indicates a degree of polarization induced by atomic displacement and dominated by the change in the orbital hybridization. Although anomalously large Born effective charge is crucial for superior piezoelectric properties of perovskite ABO 3 materials [50], the present study of decomposition of piezoelectric constant shows that coupling degree between external strain and atomic displacement is also indispensable to understand the piezoelectric properties.

Piezoelectric properties of perovskite-LiNbO 3
Next, we will show how piezoelectric properties are affected by crystal structure, while chemical composition is kept as LiNbO 3 . Various hypothetical crystal structures common for perovskite-type structure were constructed, and their energetic stabilities were examined by calculating enthalpy H = U + PV (U is total energy obtained by first-principles calculation, P is external pressure, and V is equilibrium volume under pressure P) as a function of external pressure. Imposing high pressure is most convenient method to modify crystal structure and find unexpected stable phase. The following eight types of phases were considered: Cubic, Pm-3 m; tetragonal, P4mm; and rhombohedral, R-3 m.
where names of space groups are used to distinguish each structure. Crystal structure of each phase is shown in Figure 4a. Polyhedra shown in Figure 4a correspond to Nb-centered bonding structure of Nb-O bondings. Figure 4b shows the enthalpy difference of each phases as a function of external pressure. Here, external pressure is assumed to be isotropic. Standard of enthalpy was set to be the enthalpy of most stable R3c phase under ambient condition. At the positive (compressive) pressure region, P6 3 /m phase becomes stable above 21 GPa, which is close to the experimental phase transition pressure of 25 GPa [43]. Details of the phase transition behavior under high pressure are theoretically investigated in our previous work [44]. Unfortunately, P6 3 /m phase is highly symmetric and shows no piezoelectricity. At the negative (expansive) pressure region, enthalpy difference becomes smaller as there is an increase of negative pressure except for R ¯ 3 c and P6 3 /m phases.
Imposing negative pressure can be achieved by solid solution with parent phase of larger lattice constant. At −6 GPa, P4mm phase becomes stable, while R3m and Amm2 phases become stable at −9 GPa. However, bond breaking occurs in Nb-O bonding above −6 GPa for P4mm phase. The same bond breaking occurs in Amm2 and R3m phases at −11GPa and −14 GPa, respectively. Thus, those phase transitions occur just before bond breaking.
Within the eight phases shown in Figure 4a, only P4mm, R3m, R3c, and Amm2 phases show piezoelectricity. Piezoelectric stress constant, elastic constant, and dielectric constant of P4mm, R3m, and Amm2 phases are compared with those of R3c phase in Table 3. Various piezoelectric properties are observed by each phase. Especially for P4mm and Amm2 phases, high e 33 and relatively low C 33 values are predicted, which are advantageous for large piezoelectric strain constant d 33 . On the other hand, R3m phase was found to be unstable because following mechanical stability conditions of rhombohedral symmetry: C 11 + C 12 > 0, C 33 > 0, ( C 11 + C 12 ) * C 33 > 2 C 13 , C 11 − C 12 > 0, C 44 > 0, ( C 11 − C 12 ) * C 44 > 2 C 14 (14) are broken because of C 44 < 0. This giant piezoelectric constant is almost comparable to that of PZT material [51]. Giant piezoelectric constant is understood as a result of phase instability in morphotropic phase boundary [52]. The same as P4mm phase of LiNbO 3 , we revealed that ZnO also showed anomalously large piezoelectric constant just before phase transition [30]. Finally, we would like to show phase transition path between cubic perovskite structure and LiNbO 3 structure. Figure 6a shows the energy change of Pm3m phase as a function of Li displacement along <001>, <011>, and <111> directions. Pm3m phase is paraelectric phase. Because ferroelectricity and piezoelectricity of LiNbO 3 are dominated by off-centering and displacement of Li, respectively, phase transition from Pm3m phase is also expected to be occurred by Li displacement. Li displacement along <001>, <011>, and <111> directions induces tetragonal, orthorhombic, and rhombohedral phase transition from cubic phase. Figure 6a clearly shows that tetragonal phase transition from Pm3m phase to P4mm phase is the most energetically advantageous. Figure 6b shows the phonon dispersion curve of cubic Pm3m phase of LiNbO 3 . Horizontal axis corresponds to sampling path along high symmetric reciprocal point (q-point).
Within the whole Brillouin zone of reciprocal space, unstable phonon modes with imaginary phonon frequencies are observed. Here, imaginary phonon frequency is represented as negative value for convenience. Therefore, cubic Pm3m phase of LiNbO 3 is thermodynamically unstable and considered to show phase transition in accordance with specific phonon mode of imaginary frequency (referred as soft mode). Thus, modulated structures were constructed by imposing atomic displacement along normal modes at each symmetric q-points. Modulated structures were structurally relaxed, and their space group and energy change from cubic P4mm phase were investigated. Summary of such modulated structures are shown in Table 4. At Γ point, tetragonal phase transition along with Γ 15 soft mode of cubic phase shown in Figure 3b shows energy gain of −0.422 eV/formula unit (f.u.). On the other hand, it was found that modulation at R point gives more stable energy gain of −0.682 eV/f.u. In this case, R 25 soft mode induces phase transition from cubic Pm3m phase to R ¯ 3 c phase shown in Figure 2b.     shows schematic illustration of phase transition mechanism from cubic perovskite structure to LiNbO 3 structure. On the contrary to the result of Figure 6a, R 25 soft mode is represented as rotation of NbO 6 polyhedra. Then, Γ 15 soft mode of R ¯ 3 c phase leads R3c phase, which is ground state of LiNbO 3 . Although the present study shows that perovskite-structured LiNbO 3 is thermodynamically unstable while its piezoelectricity is excellent, it can be possible to control phase transition behavior by dopant substitution.

Summary and conclusion
In this chapter, we briefly introduced sophisticated method of density functional perturbation theory. DFPT can effectively calculate the second derivative of the total energy with respect to the atomic displacement within the framework of first-principles calculation. By using DFPT method, we can predict piezoelectric constants, dielectric constants, elastic constants, and phonon dispersion relationship of any given crystal structure. Moreover, we showed our established computational technique to decompose piezoelectric constants into each atomic contribution, which enable us to gain deeper insights to understand the piezoelectricity of material. By using LiNbO 3 as a model material, we showed the predictability of DFPT for piezoelectric properties. In addition, we showed that superior piezoelectric properties are hidden in perovskite-structured LiNbO 3 . Structural relationship and possible phase transition path between LiNbO 3 structure and perovskite structure were discussed and concluded that perovskite-structured LiNbO 3 is thermodynamically unstable. Further studies are expected to control relative phase stability between perovskite and LiNbO 3 structure by dopant substation and solid solution.

Author details
Kaoru Nakamura*, Sadao Higuchi and Toshiharu Ohnuma *Address all correspondence to: n-kaoru@criepi.denken.or.jp Central Research Institute of Electric Power Industry, Japan