First principles phonon calculations in materials science

Phonon plays essential roles in dynamical behaviors and thermal properties, which are central topics in fundamental issues of materials science. The importance of first principles phonon calculations cannot be overly emphasized. Phonopy is an open source code for such calculations launched by the present authors, which has been world-widely used. Here we demonstrate phonon properties with fundamental equations and show examples how the phonon calculations are applied in materials science.


Introduction
Application of first principles calculations in condensed matter physics and materials science has greatly expanded when phonon calculations became routine in the last decade. Thanks to the progress of high performance computers and development of accurate and efficient density functional theory (DFT) codes, a large set of first principles calculations are now practical with the accuracy comparable to experiments using ordinary PC clusters. In addition to electronic structure information, a DFT calculation for solids provides energy and stress of the system as well as the force on each atom. Equilibrium crystal structures can be obtained by minimizing residual forces and optimizing stress tensors. When an atom in a crystal is displaced from its equilibrium position, the forces on all atoms in the crystal raise. Analysis of the forces associated with a systematic set of displacements provides a series of phonon frequencies. First principles phonon calculations with a finite displacement method can be made in this way. The present authors have launched a robust and easy-to-use open-source code for first principles phonon calculations, phonopy [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]. The number of users is rapidly growing world-wide, since the information of phonon is very useful for accounting variety of properties and behavior of crystalline materials, such as thermal properties, mechanical properties, phase transition, and superconductivity. In this article, we show examples of applications of the first principles phonon calculations.
In Sections 2-4, we take FCC-Al as examples of applications of first principles phonon calculations. For the electronic structure calculations, we employed the planewave basis projector augmented wave method [17] in the framework of DFT within the generalized gradient approximation in the Perdew-Burke-Ernzerhof form [18] as implemented in the VASP code [19,20,21]. A plane-wave energy cutoff of 300 eV and an energy convergence criteria of 10 −8 eV were used. A 30 × 30 × 30 k-point sampling mesh was used for the unit cell and the equivalent density mesh was used for the supercells together with a 0.2 eV smearing width of the Methfessel-Paxton scheme [22]. For the phonon calculations, supercell and finite displacement approaches were used with 3×3×3 supercell of the conventional unit cell (108 atoms) and the atomic displacement distance of 0.01Å.

Harmonic approximation
In crystals, it is presumed that atoms move around their equilibrium positions r(lκ) with displacements u(lκ), where l and κ are the labels of unit cells and atoms in each unit cell, respectively. Crystal potential energy Φ is presumed to be an analytic function of the displacements of the atoms, and Φ is expanded as where α, β, · · · are the Cartesian indices. The coefficients of the series expansion, Φ 0 , Φ α (lκ), Φ αβ (lκ, l κ ), and, Φ αβγ (lκ, l κ , l κ ), are the zero-th, first, second, and third order force constants, respectively. With small displacements at constant volume, the problem of atomic vibrations is solved with the second-order terms as the harmonic approximation, and the higher order terms are treated by the perturbation theory. With a force F α (lκ) = − ∂Φ ∂uα(lκ) , an element of secondorder force constants Φ αβ (lκ, l κ ) is obtained by Crystal symmetry is utilized to improve the numerial accuracy of the force constants and to reduce the computational cost. The more details on the calculation of force constants are found in Ref. [4,5].
As found in text books [23,24,25,26], dynamical property of atoms in the harmonic approximation is obtained by solving eigenvalue problem of dynamical matrix D(q), where m κ is the mass of the atom κ, q is the wave vector, and j is the band index. ω qj and e qj give the phonon frequency and polarization vector of the phonon mode labeled by a set {q, j}, respectively. Since D(q) is an Hermitian matrix, its eigenvalues, ω 2 qj , are real. Usually D(q) is arranged to be a 3n a × 3n a matrix [26], where 3 comes from the freedom of the Cartesian indices for crystal and n a is the nubmer of atoms in a unit cell. Then e qj becomes a complex column vector with 3n a elements, and usually e qj is normalized to be 1, i.e., ακ |e ακ qj | 2 = 1. e qj contains information of collective motion of atoms. This may be understood as a set of atomic displacement vectors, where A is the complex constant undetermined by Eq. (3), and e κ qj T = (e xκ qj , e yκ qj , e zκ qj ). As a typical example, the phonon band structure and phonon density of states (DOS) of Al are shown in Fig. 1. The phonon DOS is defined as where N is the number of unit cells in crystal. Divided by N , g(ω) is normalized so that the integral over frequency becomes 3n a . The phonon band structure can be directly comparable with experimental data by neutron or X-ray inelastic scattering. They often show reasonable agreements [16,27,28]. Frequency data by Raman and infrared (IR) spectroscopy can also be well reproduced [8,29]. Irreducible representations of phonon modes, which can be used to assign Raman or IR active modes, are calculated from polarization vectors [8,30]. Atom specific phonon DOS projected along a unit direction vectorn is defined as This g κ (ω,n) can be directly compared with that measured by means of nuclear-resonant inelastic scattering using synchrotron radiation. In Ref. [13], phonon calculations of L1 0 -type FePt projected along the c-axis and basal plane are well comparable to experimental 57 Fe nuclearresonant inelastic scattering spectra measured at 10 K in the parallel and perpendicular geometries, respectively.  Once phonon frequencies over Brillouin zone are known, from the canonical distribution in statistical mechanics for phonons under the harmonic approximation, the energy E of phonon system is given as where T , k B , andh are the temperature, the Boltzmann constant, and the reduced Planck constant, respectively. Using the thermodynamic relations, a number of thermal properties, such as constant volume heat capacity C V , Helmholtz free energy F , and entropy S, can be computed as functions of temperature [26]: and The calculated F , C V , and S for Al are shown in Fig. 2.

Mean square atomic displacements
With the phase factor convention of the dynamical matrix used in Eq. (4), an atomic displacement operator is written as, whereâ † andâ are the creation and annihilation operators, respectively. Distributions of atoms around their equilibrium positions are then obtained as the expectation values of Eq. (12). The mean square atomic displacement projected alongn is obtained as Eq. (13) is lattice-point (l) independent since the phase factor disappears. Anisotropic atomic displacement parameters (ADPs) to estimate the atom positions during thermal motion can also be computed and compared with experimental neutron diffraction data. Thermal ellipsoids may be discussed using mean square displacement matrix B(κ) defined by The shape and orientation of an ellipsoid is obtained solving eigenvalue problem of this matrix. The method has been applied to show the ORTEP (Oak Ridge Thermal Ellipsoid Plot)-style drawing of ADPs [14]. Ref. [7] shows an example for a ternary carbide Ti 3 SiC 2 having a layered structure known as MAX phases, in which we can see good agreement between calculated and experimental ADPs.

Quasi-harmonic approximation
By changing volume, phonon properties vary since the crystal potential is an anharmonic function of volume. In this article, the term "quasi-harmonic approximation (QHA)" means this volume dependence of phonon properties, but the harmonic approximation is simply applied at each volume. Figure 3a shows calculated phonon frequencies of Al at X and L points with respect to ten different unit-cell volumes. Typically phonon frequency decreases by increasing volume, and the slope of each phonon mode is nearly constant in the wide volume range. The normalized slope is called mode-Grüneisen parameter that is defined as Once dynamical matrix is known, γ qj is easily calculated from the volume derivative of the dynamical matrix [25], The quantity can be related to macroscopic Grüneisen parameter γ using mode contributions to the heat capacity C qj found in Eq. (9), γ = qj γ qj C qj /C V [24,26].
Silicon is known as a famous exception to have large negative mode-Grüneisen parameters. Mode-Grüneisen parameter is a measure of anharmonicity of phonon modes and is related to third-order force constants directly [25]. Therefore crystals that possess large anharmonic terms beyond third-order terms in Eq. (1) can show non-linear change of phonon frequency with respect to volume. This is often observed in crystals that exhibit second-or higherorder structural phase transitions [2].
The phonon frequency influences the phonon energy through Eq. (8). The thermal properties are thereby affected. Using thermodynamics definition, thermodynamic variables at constant volume is transformed to those at constant pressure that is often more easily measurable in experiments. Gibbs free energy G(T, p) at given temperature T and pressure p is obtained from Helmholtz free energy F (T ; V ) through the transformation, where the right hand side of this equation means finding a minimum value in the square bracket by changing volume V . We may approximate F (T ; V ) by the sum of electronic internal energy U el (V ) and phonon Helmholtz free energy F ph (T ; V ), i.e., F (T ; V ) U el (V ) + F ph (T ; V ). U el (V ) is obtained as total energy of electronic structure from the first principles calculation, and the first principles phonon calculation at T and V gives F ph (T ; V ). The calculated U el (V ) + F ph (T ; V ) are depicted by the filled circle symbols in Fig. 3b, where the ten volume points chosen are the same as those in Fig. 3a. The nine curves are the fits to equation of states (EOS) at temperatures from 0 to 800 K with 100 K step. Here the Vinet EOS [34] was used to fit the points to the curves though any other functions can be used for the fitting. The minimum values at the temperatures are depicted by the cross symbols in Fig. 3b and are the Gibbs free energies at the temperatures and the respective equilibrium volumes are simultaneously given. Volumetric thermal expansion coefficient, ∂T , is obtained from the calculated equilibrium volumes V (T ) at dense temperature points. β(T ) for Al is shown in Fig. 3c, where we can see that the calculation shows reasonable agreements with the experiments. In thermodynamics, heat capacity at constant pressure C P is given by In Eq. (18), the second term of the second equation is understood as the contribution to heat capacity from thermal expansion. C P for Al is shown in Fig. 2. The agreement of the calculation with the experiment is excellent. At high temperatures, the difference between C P and C V is not negligible in Al. Therefore it is essential to consider thermal expansion for heat capacity. QHA is known as a reasonable approximation in a wide temperature range below melting point except for temperatures very close to melting point where higher-order terms in Eq. (1) become evident [35].

Stability condition and imaginary mode
At equilibrium, ∂Φ ∂rα(lκ) = 0, a crystal is dynamically (mechanically) stable if its potential energy always increases against any combinations of atomic displacements. In the harmonic approximation, this is equivalent to the condition that all phonons have real and positive frequencies [25]. However under virtual thermodynamic conditions, imaginary frequency or negative eigenvalue can appear in the solution of Eq. (3). This indicates dynamical instability of the system, which means that the corrective atomic displacements of Eq. (5) reduce the potential energy in the vicinity of the equilibrium atomic positions. Imaginary mode provides useful information to study displacive phase transition. A typical example is shown in Figs. 4a to 4c [12]. Imaginary modes can be found only for β-Ti, that has BCC structure, at both P and N points. This indicates that β-Ti is unstable at low temperature. Such imaginary modes cannot be seen for either ω-Ti whose crystal structure is shown in Fig. 4d or α-Ti (HCP). Experimentally β-Ti is known to occur above 1155K. At such high temperatures, large atomic displacements can stabilize the BCC structure. In such a case, the perturbation approach is invalid. Phonons with large atomic displacements may be treated by self-consistent phonon method [25,36] or by a combination of molecular dynamics and lattice dynamics calculation [37,38,39], which is not discussed in this article.
A given structure having imaginary phonon modes can be led to alternative structures through continuous atomic displacements and lattice deformations. The present authors systematically investigated the evolution of crystal structures from the simple cubic (SC) structure [6]. The inset of Fig. 5 shows the phonon band structure of SC-Cu (P m3m). Imaginary modes can be found at M (1/2, 1/2, 0) and X(1/2, 0, 0) points. Then the SC structure is deformed along these directions. In order to accommodate the deformation in the calculation model with the periodic boundary condition, the unit cells are expanded by 2×2×1 for the M point and 2×1×1 for the X point. The M point deformation breaks the crystal symmetry of SC (P m3m) to P 4/nmm. The doubly degenerated instability at the X point leads to P mma and Cmcm as highest possible symmetries. The deformed crystal structures are relaxed by first principles calculations imposing the corresponding space-group operations. After these procedures, body-centered tetragonal (BCT), simple hexagonal (SH), and FCC are respectively formed. The whole procedure finishes when all crystal structures at the end-points are found to be dynamically stable. Finally a treelike structure of crystal structure relationships was drawn as shown in Fig. 5, where thick lines indicate phase transition pathways (PTPs). The space-group type written near a line is a common subgroup of initial and final structures. The presence of the line indicates that the energy decreases monotonically with the phase transition. In other words, the transition can take place spontaneously without any energy barrier. The line is terminated when the final structure is dynamically stable. In the line diagram, ω is located at the junction of two pathways, i.e., ω → BCC → HCP and ω → FCC. The instability of ω at the Γ point leads to BCC, which is still dynamically unstable and eventually leads to HCP. Another instability at the M point leads to FCC. The other instability at the K point, which is doubly degenerate, leads to FCC. On the path from ω to BCC, the crystal symmetry of ω having the space-group type of P 6/mmm is once lowered to P3m1 and then becomes Im3m (BCC) after the geometry optimization. Both ω-and BCC-Cu are dynamically unstable, which can be formed only under crystal symmetry constraints. FCC-Cu is, of course, dynamically stable. Several PTPs between BCC and FCC have been proposed in literature. However, they are mostly based only upon investigation of continuous lattice deformation. For example in the classical Bain path, formation of BCT in between BCC and FCC is considered. Formation of SC is taken into account in the "trigonal Bain path." Normal modes of phonon, which should be most adequate to describe the collective atomic displacements, have not been considered. The presence of ω as the lowest energy barrier in the BCC-FCC pathway had not been reported before Ref. [6]. The situation is the same for the BCC-HCP transition, known as the Burgers path. The Burgers path was thought to be quite complicated from the viewpoint of the lattice continuity. On the basis of the present study, it can be easily pointed out that the BCC-HCP transition pathway is along the space-group type of Cmcm.
Evolution diagram was constructed in the same way for NaRTiO 4 (R: rare-earth metal) with Ruddlesden-Popper type structure [9]. Inversion symmetry breaking by oxygen octahedral rotations was unambiguously demonstrated. The mechanism is expected to lead to many more families of acentric oxides.

Interactions among phonons and lattice thermal conductivity
Using the harmonic phonon coordinates, anharmonic terms in Eq. (1) are transformed to a picture of phononphonon interactions [4,24]. Lattice thermal conductivity can be accurately calculated by solving linearized Boltzmann transport equation with the phonon-phonon interaction strength obtained using first principles calculation [5,40,41]. Although the computational cost for such calculations is many orders of magnitudes higher than the ordinary DFT calculations of primitive cells, such calculations have already been applied for many simple compounds and computed lattice thermal conductivities show good agreements with experimental data [5,41]. Calculations with special focus on searching thermoelectric materials have also been made [10,16,41].