3.10. Vibrations in Solids - Phonons

download:pdf

3.10.1. Overview

MedeA Phonon is a computational tool for deriving vibrational properties of solids from first principles. This section covers the underlying physical theorems and equations, provides a comprehensive overview on methods and numerical techniques.

  • Optionally impose translational and rotational invariance on forces to improve accuracy.
  • Modeling longitudinal optical branches close to LO/TO splitting with an approximate scheme based on Born effective charges and dielectric constants for polar crystals.
  • Analyze irreducible representations of phonon modes at the \(\Gamma\)-point as well as Infrared and Raman activities.
  • Debye-Waller factors and incoherent inelastic neutron scattering based on off-diagonal phonon density of states
  • Thermodynamic functions decomposed into contributions of individual atoms and polarization directions: internal energy E, free energy F, entropy S and heat capacity CV as a function of temperature.
  • Intensity (form factors) of coherent neutron scattering in different Brillouin zones and different geometries.

Doubly differential incoherent inelastic neutron scattering cross sections of one Phonon processes for mono-crystals and of up to five Phonon processes for poly-crystals.

3.10.2. Background

A detailed knowledge of lattice vibrations is critical for the understanding and quantitative prediction of a wide variety of physical properties of solids. The interpretation of optical, infrared, Raman and synchrotron spectroscopic methods as well as neutron scattering experiments rely on an accurate theoretical approach to lattice dynamics. The fundamental thermodynamic functions of internal and free energy, entropy, heat capacity as well as non-linear properties such as thermal expansion and heat conduction are to a considerable extent determined by the vibrations of the constituent atoms in the lattice. Further phenomena such as electrical resistivity of metals, superconductivity and the temperature dependence of optical spectra are governed by electron-Phonon coupling. Fortunately, the quantum theory of lattice dynamics is well developed and has proven to be one of the most successful theories of solid state physics and therewith provides an impressive reassurance of the validity and applicability of the concepts of quantum theory in general.

From a historical point of view, the fundamental theory of lattice dynamics was developed in the 1930’s and was reviewed in the book of Born and Huang [1], which is still regarded as a standard textbook in the field. In those early days, the main interest focused on general properties of the dynamical matrix, e.g. its symmetry and analytical properties, but less emphasis was put on its underlying physical origin based on the interactions of the constituent ions and electrons. The first approaches to theoretically determine phonon dispersion relations were based on lattice-dynamical models with adjustable parameters to be fitted to experimentally measured frequencies, e.g. from inelastic neutron scattering experiments, and were mainly used as a means to interpolate between experimental data points. These model potentials either did not account for electronic polarization at all, like in the force-constant model, Born-von Kármán model, Born-Mayer potentials for alkali halides or Lennard-Jones potentials for noble gas crystals [1], or polarization effects were approximated by electric dipoles generated by relative displacements of the shells of valence electrons relative to the ion cores (shell models [2]). In addition, to tackle systems with highly anisotropic electron distributions as found in covalently bonded crystals, alternative models have been developed, relying on angular forces (Keating models [3]) or bond charges (bond charge model [4]). The predictive capabilities of these methods turned out to be quite limited. Although experimental phonon dispersions are accurately reproduced close to the data used for fitting - especially along high symmetry directions - predictions along lower-symmetry directions frequently failed due to the lack of a rigorous description of the interactions among constituent atoms. The first attempts to investigate the connections between the electronic structure and the dynamical properties were made as late as in the early seventies and were based on dielectric matrices or the electron polarizability [5].

With the advent of density functional theory [6] and the progress with numerical methods for solving quantum physical equations together with the emergence of more and more powerful computers made it feasible to accurately describe the interatomic interactions in crystals and molecules based on quantum mechanics. It is thus possible, nowadays, to calculate a large variety of properties of simple materials without relying on input data from experiments. These so-called first-principles or ab-initio methods also revolutionized the practical approaches to lattice dynamics and accurate theoretical phonon dispersion curves can be calculated completely independent of any experimental knowledge.

Three different techniques for ab initio evaluation of vibrational properties have been developed, namely

  • direct methods based on total energy changes or forces calculated for atoms displaced from their equilibrium position,
  • analytical calculation of force constants based on a perturbative expansion around the equilibrium geometry,
  • Fourier transform of the atomic velocity autocorrelation function obtained from a molecular dynamics trajectory. [7]

The third option usually suffers from difficulties to reach equilibrium within reasonable simulation times, although several technical tricks have been developed to address these problems (e.g. Nosé-Hoover thermostat [8] and multiple signal classification (MUSIC) algorithm [9].

The second option makes use of the fact that the lattice distortion associated with a phonon is a static perturbation acting on the electrons, causing a linear response of the electron density which determines the energy variation up to third-order (\(2 \nu + 1\) theorem [10]). The linear variation of the electron density is determined in terms of the linear response of the Kohn-Sham orbitals applying first-order density functional perturbation theory [11]. This method is assessed as quite accurate, efficient and elegant, however, requiring highly specialized ab initio computer codes and considerable implementation efforts. For a review of the formalism and implementation of density functional perturbation theory and its application see Baroni, de Gironcoli, Corso & Giannozzi [12].

Direct methods (option (1) above) require the evaluation of total energy and forces for the equilibrium geometry as well as of several distorted geometries from which the force constant matrix can be assembled. Several different techniques are subsumed under the title direct methods: in the early eighties, the first attempts for ab initio lattice dynamics simulation were based on frozen phonon total energy calculations to derive accurate phonon frequencies for specific high symmetry points in the Brillouin zone (e.g. Yin & Cohen [13]). Phonon dispersion curves along specific high symmetry directions in reciprocal space were determined by the method of interplanar force constants [14], where planes perpendicular to these directions are displaced within an elongated supercell. Since these one-dimensional force constants between high-symmetry planes are linear combinations of the three-dimensional interatomic force constants, the phonon dispersion can be derived, in principle, for general directions [15].

The most general direct approach to lattice dynamics is based on the ab initio evaluation of forces on all atoms produced by a set of finite displacements of a few atoms within an otherwise perfect crystal. The perfect crystal environment has to be sufficiently large to ensure that interactions of the perturbation with all its translational symmetry equivalent copies are small, which usually requires construction of suitable supercells and atomic displacements, assembling force constant matrices from the calculated forces and calculating phonon dispersion relations via Fourier transforms [16]. Phonon dispersion curves determined from direct approaches are very accurate provided the studied crystal area (supercell) is larger than the interaction range of all constituent atoms. The main advantages of direct methods consist in the fact that no specialized ab initio codes are required, as long as forces can be derived, and that anharmonic effects can be treated straightforwardly. The necessity to consider supercells in the case of simple crystal structure with small unit cell extensions is one of the main drawbacks with regards to the linear response methods.

The MedeA Phonon module is based on the general direct approach to lattice dynamics and is designed to work independently of a specific underlying code for deriving forces and total energies. However, together with VASP, a fully automatic and highly parallel procedure is provided within MedeA.

3.10.3. Theory of Lattice Dynamics

In this section the basic theory of lattice dynamics is summarized and the essential analytical and numerical aspects of the direct approach to the calculation of phonon dispersions are discussed in detail. Furthermore, relevant equations for calculating form factors, thermodynamic functions, thermal displacements, and neutron scattering cross-sections as implemented in the Phonon module are provided. For further details of the theory as well as applications in diverse research fields, we refer to several recent review articles [17].

3.10.3.1. Crystal Symmetry

The crystal structure is subject to space group symmetry

(1)\[G = \lbrace h_{1} \vert \tau_{1} \rbrace, \lbrace h_{2} \vert \tau_{2} \rbrace, ..., \lbrace h_{g} \vert \tau_{g} \rbrace\]

\(h_{i}\) are the point group symmetry elements and \(\tau_{i}\) are partial translations. The point group of the crystal \(\hat{G}\) contains all point symmetry elements \(h_{i}\). All atomic positions can be generated from non-equivalent sites applying these symmetry operations. For each non-equivalent atom \(m\) with site point group \(P_{\mu}{p_{i}}\) the space group can be decomposed into left cosets with respect to \(P_{m}\)

(2)\[G = \lbrace h_{1}^{\mu} \vert \tau_{1}^{\mu} \rbrace P_{\mu} + \lbrace h_{2}^{\mu} \vert \tau_{2}^{\mu} \rbrace P_{\mu} + ... + \lbrace h_{g/l}^{\mu} \vert \tau_{g/l}^{\mu} \rbrace P_{\mu}\]

Applying the generating symmetry elements \(\lbrace h_{1}^{\mu} \vert \tau_{1}^{\mu} \rbrace\) all equivalent atom positions are found.

3.10.4. Basic Relations and Approximations

The ground state energy E of a crystal as a function of atom positions R(n,m) (where n denotes the unit cell and m the atom index) can be written as a Shubham expansion in terms of atomic displacements around the equilibrium positions

(3)\[E \left( ...R(n,\mu)...R(m,\nu) \right) = E_{0} + \dfrac{1}{2} \sum_{n,\mu,m,\nu} \Phi(n,\mu;m,\nu) U(n,\mu) U(m,\nu)+ O(U^{3})\]

\(E_{0}\) is the equilibrium energy. Terms linear in U are absent due to the condition that at equilibrium all forces on the atoms vanish. The elements of the atomic force constant matrix are defined by

(4)\[\Phi_{i,j}(n,\mu;m,\nu) = \left[ \dfrac{\partial ^{2}E}{\partial R_{i} (n,\mu) \partial R_{j}(m,\nu)} \right] _{0}\]

with the gradients taken at the minimum energy configuration for which all first-order derivatives vanish. The frequently used harmonic approximation consists in retaining in (3) only terms up to quadratic order in the displacements.

Within the harmonic approximation the classical equations of motion for each atom are

(5)\[M_{\mu} \ddot{U} (n,\mu) = \sum_{n,\mu,m,\nu} \Phi (n,\mu;m,\nu)U(m,\nu)\]

where \(M_{\mu}\) is the atom mass of atom \(\mu\). Solutions are required to exhibit Bloch-wave form because of translational invariance.

(6)\[U(n,\mu) = \dfrac{1}{M_{\mu}} e(k) e^{ikR(n,\mu)-i \omega t}\]

The k vectors are chosen to fulfill Born-von Karman periodic boundary conditions.

The equation of motion can thus be cast into the simple form

(7)\[D(k)e(k,j) = \omega^{2}(k,j)e(k,j)\]

where for each mode j the phonon frequencies \(\omega^{2}(k,j)\) are the eigenvalues and the polarization vectors \(e(k,j)\) are the eigenvectors of the dynamical matrix \(D(k)\), which has been introduced as the discrete Fourier transform

(8)\[D(k;\mu,\nu) = \dfrac{1}{\sqrt{M_{\mu}M_{\nu}}} \sum_{m} \Phi(0,\mu;m,\nu) exp \left( 2 \pi i k \left( R(0,\mu) - R(m,\nu) \right) \right)\]

The summation \(m\) runs over all atoms of the crystal. The forces \(F(n,\mu)\) on all atoms generated by the displacement \(U(m,\nu)\) of atom \(\nu\) can be written as

(9)\[F_{i}(n,\mu) = \sum_{m,\mu,j} \Phi_{i,j} (n,\mu;m,\nu) U_{j} (m,\nu)\]

which relates the generated forces with the force constant matrices and atomic displacements. This is the central relation of the direct method.

The complete quantum mechanical description for a system of ions and electrons requires additional approximations. Since the nuclear mass is much larger than that of an electron, it is reasonable to consider the nuclei in their equilibrium position while dealing with the electronic motion. In mathematical terms, this so-called adiabatic (or Born-Oppenheimer) approximation yields a separation of the general Schrödinger equation into one for the motion of nuclei and one for the electrons. The nuclear motion is determined by the potential field generated by the average motion of the electrons and the corresponding nuclear Schrödinger equation yields formally the same results as the classical theory. The above equations remain valid, however, with the exact effective potential. The purely electronic Schrödinger equation can be solved within density functional theory.

3.10.4.1. Symmetry Properties of Force Constants

The atomic force constants exhibit certain symmetry properties related to the crystal symmetry. Because of the translational invariance of the crystal, force constants depend only on the difference \(R(n,m)-R(m,n)\) and satisfy the acoustic sum rule

(10)\[\sum_{n,\mu,n,\nu} \Phi (n,\mu;m,\nu) = 0\]

which expresses the invariance of the potential energy on uniform translations of the whole crystal.

Therefore, a force constant \(\Phi(m,\mu,n,\nu)\) represents a bond between atoms \(m,\mu\) and \(n,\nu\). The bond point group \(B= \lbrace b_{1}b_{2}...b_{b} \rbrace\) leaves this specific bond invariant and its elements \(b_{i}\) span a subgroup of the point group of the crystal. The symmetry of force constants, therefore, is set by transformations by \(b_{i}\) which reduces the number of independent parameters. Furthermore, each force constant can be decoupled into a part \(A\) determined exclusively by symmetry and a part \(P\) which is dependent on potential strength. The potential parameters contained in \(P\) are no longer dependent on any other element. This decoupling allows the use of fitting procedures on a minimal set of independent parameters. The crystal space group can be decomposed with respect to the bond point group \(B\)

(11)\[G = \lbrace h_{1}^{(b)} \vert \tau_{1}^{(b)} \rbrace B + \lbrace h_{2}^{(b)} \vert \tau_{2}^{(b)} \rbrace B + ... + \lbrace h_{g/l}^{(b)} \vert \tau_{g/l}^{(b)} \rbrace B\]

Therefore, knowing the symmetry of one force constant, all remaining equivalent ones can be created by applying the generating symmetry elements \(\lbrace h_{i}^{(b)} \vert \tau_{i}^{(b)} \rbrace\).

3.10.4.2. Supercell Force Constants for the Direct Method

The ab initio calculations for the direct approach are performed for a supercell constructed as a multiplication of the original primitive cell with periodic boundary conditions. There are symmetry implications for the force constants resulting from the supercell geometry. The symmetry elements of the supercell point group \(\hat{S} = \lbrace s_{1}, s_{2}, ... s_{s} \rbrace\) leave the supercell as parallelepiped (not as crystal) invariant. The point group symmetry of the supercell crystallite \(\hat{H}\) consists of the common point symmetry elements of the supercell point group \(\hat{S}\) and the point group of the crystal \(\hat{G}\), i.e. \(\hat{H} = \hat{S} \cap \hat{G}\). The supercell crystallite space group H (including translations and partial translations) can be equal to the crystal space group \(G\), or it can be a subgroup of \(G\). In the second case, some equivalent atoms become non-equivalent, and site and bond point groups are reduced to corresponding subgroups. As a consequence, the symmetry of phonon modes derived from these supercell calculations might be reduced and degeneracy of modes may be decreased.

The direct method for calculating vibrational properties starts with an optimization of the atom positions (and lattice parameters) of the unperturbed supercell introduced in the previous paragraph. Any computational method providing atomic forces is suitable for the direct approach, but usually first-principles methods are employed. After the initial total energy minimization process, all forces should vanish, and the resulting structure is the reference and starting point for all further steps. In order to obtain complete information on all force constants, it is necessary to displace each non-equivalent atom of the supercell along three non-equivalent directions, and to calculate for each of these perturbed supercells the forces on all other atoms generated by this displacement. The non-equivalent atoms and directions are defined with respect to the supercell crystallite space group H. Further reduction of the number of displacements is possible: if the site point group of the atom is cubic, a single displacement along a single fourfold axis is adequate. If the site point group is tetragonal, one displacement along the fourfold symmetry axis and one perpendicular to it are sufficient.

The use of supercells with periodic boundary conditions rather than interaction shells has some consequences on the definition of force constants. The displacement \(U(m,\nu)\) of an atom \((m,\nu)\) reappears for all atoms \((m+L,\nu)\) in all images L of the supercell due to periodicity. \(L=(L_{a},L_{b},L_{c})\) are the lattice constant indices of the supercell. Thus, the displacement of a single atom \((m,\nu)\) in the original supercell generates on atom \((n,\mu)\) a net force

(12)\[F_{i}(n,\mu) = - \sum_{L} \Phi_{i,j}(n,\mu;m+L,\nu)U_{j}(m,\nu) = -\Phi_{i,j}^{(\Sigma)}(n,\mu;m,\nu)U_{j}(m,\nu)\]

where the cumulative force constant is defined as

(13)\[\Phi_{i,j}^{\Sigma} (n,\mu;m,\nu)) = \sum_{L}\Phi_{i,j} (n,\mu;m+L,\nu)\]

and the summation L runs overall supercell images. In order to make sure that all neighbors of a given interaction shell are taken into account for constructing the dynamical matrix, the atom \((n,\mu)\) under consideration is positioned at the center of an extended supercell, which has the same size as the original one but includes all atoms on its surface planes, edges, and corners. The extended supercell, in general contains more atoms than the conventional one. If the displaced atoms \((M+L,\nu)\) for several \(L\) are located at the same distance from the considered atom \((n,\mu)\), i.e. if they are located at the surface of the extended supercell, the corresponding force constant has to be scaled by the factor \(1/n_{m}\), where \(n_{m}\) is the number of equivalent atoms at the surface. Thus, the supercell force constants are defined as

(14)\[\Phi_{i,j}^{(SC)} (n,\mu;m,\nu) = \dfrac{1}{n_{m}} \sum \Phi_{i,j}^{\Sigma} (n,\mu;m+L,\nu)\]

and the net force can be written as

(15)\[F_{i}(n,\mu) = -\sum_{m,\nu \in SC} \Phi_{i,j}^{(SC)} (n,\mu;m,\nu) U_{j}(m,\nu)\]

where the summation is limited to atoms of the extended supercell. The supercell force constants may have higher symmetry than the conventional ones, provided by the symmetry relations among the equivalent surface atoms.

The above equation has to be solved in the direct method, yielding the supercell force constants from the known forces and displacement variables. As stated in section Symmetry Properties of Force Constants for force constants in general, also the supercell force constants can be decoupled into a matrix A determined by symmetry and the matrix P of independent potential parameters. The net force is

(16)\[\begin{split}F(n,m) & = - \sum_{m,\nu \in SC} A(n,\mu;m,\nu) P^{(SC)} (n,\mu;m,\nu) U(m,\nu) \\ & = C(n,\mu;m,\nu)P^{(SC)}(n,\mu;m,\nu)\end{split}\]

where the matrix \(C\) is defined as

\[C(n,\mu;m,\nu) = - \sum_{(m,\nu) \in SC} A(n,\mu;m,\nu) U(m,\nu)\]

The main equation of the direct method is cast into the simple global form

(17)\[F = C \cdot P^{(SC)} \Rightarrow P^{(SC)} = C^{-1} \cdot F\]

Usually, the number of forces is larger than the number of independent parameters, and the system is over-determined. Equation (17) is solved by an algorithm for the singular value decomposition of the sparse matrix \(C\) [18]. All supercell force constants can be expressed in terms of the potential parameters \(P^{(SC)}\).

3.10.4.3. Translational-Rotational Invariance

The forces should satisfy automatically the translational-rotational invariance conditions guaranteeing that all acoustic modes start at zero frequency at the G-point. Due to numerical errors of the forces as obtained in actual ab initio calculations, this condition is not exactly fulfilled. It is possible, however, to enforce the invariance conditions in the derivation of the force constants.

For the case of supercell force constants with an interaction range confined inside the supercell, the translational-rotational invariance conditions are

(18)\[\sum_{(m,\nu) \in (SC)} \Phi_{i,j}^{(SC)} (n,\mu;m,\nu) = 0\]
(19)\[\sum_{(m,\nu) \in (SC)} \left[ \Phi_{i,j}^{(SC)} (n,\mu;m,\nu) R_{k}(m,\nu) - \Phi_{i,j}^{(SC)} (n,\mu;m,\nu) R_{j}(m,\nu) \right] = 0\]

where the summation is limited to the atoms of the extended supercell, \(R(m,n)\) are position vectors and \(i,j,k\) are Cartesian indices. In the simple global form, the invariance conditions can be written as \(O = M \cdot P^{(SC)}\), where \(M\) is a \((18n \times p)\) matrix. Here, \(n\) is the number of atoms in the primitive unit cell, \(p\) is the number of potential parameters and 18 is the sum of 9 translational and 9 rotational equations. This equation may be added to the main equation of the direct method and the resulting equations can be solved by singular value decomposition. The \(b\) factor specifies the strength of enforcement of the invariance conditions.

(20)\[\begin{split}\begin{pmatrix} F \\ 0 \\ \end{pmatrix} = \begin{pmatrix} C \\ \beta M \\ \end{pmatrix} \cdot P^{(SC)}\end{split}\]

3.10.4.4. Supercell Dynamical Matrix

The direct method yields supercell force constants from the forces and by discrete Fourier transform the supercell dynamical matrix is obtained:

(21)\[D^{(SC)} (k;\mu,\nu) = \dfrac{1}{\sqrt{M_{\mu}M_{\nu}}} \sum_{m \in SC} \Phi^{(SC)} (0,\mu;m,\nu) e^{-2 \pi i k \cdot \left( R(0,\mu) - R(m,\nu) \right)}\]

Here, the atom \((0,m)\) is located at the center of the extended supercell and summation is limited to all neighbors of this cell. The summation over images of the supercell is alcreated included in the supercell force constant.

The supercell dynamical matrix becomes identical to the conventional dynamical matrix for all \(k\) vectors if the supercell is large enough such that the interaction range is confined inside. In this case, all phonon frequencies are calculated exactly. If the interaction exceeds the extended supercell size, both dynamical matrices are identical only for special reciprocal vectors \(k_{s}\) fulfilling the condition

(22)\[exp \lbrace -2 \pi i k_{s} \cdot a_{SC} \rbrace = 1 \ \ exp \lbrace -2 \pi i k_{s} \cdot b_{SC} \rbrace = 1 \ \ exp \lbrace -2 \pi i k_{s} \cdot c_{SC} \rbrace = 1\]

where \(a_{SC}\), \(b_{SC}\) and \(c_{SC}\) are the basis vectors of the extended supercell. For these special vectors \(k_{i}\) the phonon frequencies are calculated exactly. Furthermore, the accuracy of phonon frequencies mainly depends on the mutual symmetries of the supercell and the wave vector. The precision of phonon frequencies for a discrete wave vector \(k_{i}\) can be assessed by considering the point groups \(\hat{G}\), \(\hat{S}\) and \(\hat{G}_{ks}\) of the crystal space group \(G\), the supercell and the wave vector \(k_{i}\), respectively. The accuracy is

(23)\[\textrm{Exact if} \quad \hat{S}=\hat{G} \textrm{ and } \hat{G}_{k_{z}} \subseteq \hat{S}\]
(24)\[\textrm{High precision} \quad \hat{S} \subset \hat{G} \textrm{ and } \hat{G}_{k_{z}} \subseteq \hat{S}\]
(25)\[\textrm{Low precision} \quad \hat{S}=\hat{G} \textrm{ and } \hat{G}_{k_{z}} \not\subset \hat{S}\]

For wave vectors fulfilling the requirements for exact solutions, the phonon frequencies and degeneracy are exactly calculated from the supercell. For wave vectors satisfying the requirement for high precision, the supercell approach usually yields highly precise phonon frequencies but the exact mode degeneracy is not guaranteed. As a consequence, modes along equivalent directions with respect to \(\hat{G}\) but non-equivalent directions with respect to \(\hat{S}\) might have different frequencies and degeneracy. Modes for \(k_{s}\) vectors fulfilling the criterion for low precision suffer from the same deficiencies but in addition may also have poorer frequency values, because the list of neighbors is far from complete. All intermediate phonon branches are interpolated between these special vectors \(k_{s}\). Supercells providing the maximal number of exact wave vectors are: \(n \times n \times n\) for cubic, \(n \times n \times m\) for tetragonal, \(n \times n \times n\) for hexagonal (rhombohedral), and \(n \times m \times l\) for orthorhombic, monoclinic and triclinic symmetry, respectively. The shape of the supercell should be as close as possible to a cube to obtain accurate frequencies for all directions; elongated supercells may be applied to achieve high accuracy in the direction of elongation.

3.10.4.5. Phonon Dispersion and Polarization Vectors

The frequencies \(w^{2}(k,j)\) of phonon modes \(j\) are calculated by diagonalization of the supercell dynamical matrix for each wave vector \(k\) along a specified path through the Brillouin zone, thus creating phonon dispersion curves.

(26)\[D(k) \cdot e(k,j) = \omega ^{2} (k,j) e (k,j)\]

The irreducible representations of all phonon modes at the \(G\) point can be calculated, providing in addition Raman and infrared activities of the modes. The complex polarization vectors satisfy the orthonormality relations

(27)\[\sum e_{i}^{\ast} (k,j;\mu) \cdot e_{i} (k,j;\nu) = \delta_{i,j} \delta_{\mu,\nu}\]
(28)\[\sum e_{i}^{\ast} (k,j;\mu) \cdot e_{i} (k,j;\mu) = \delta_{i,j}\]

The polarization vectors \(e(k,j;)\) defined for the wave vector \(k\) centered at the origin of reciprocal space differ from the conventional polarization vector as defined for the wave vector \(k_{tau}\) pointing from the center of a given Brillouin zone labeled by the reciprocal vector \(\tau\). Because of \(k=\tau+k_{tau}\) the relation between these differently defined polarization vectors is

(29)\[e(k,j;\mu) = e(k_{\tau},j;\mu) exp (-2 \pi \tau \cdot r_{\mu})\]

Using the polarization vectors, the displacements caused by a particular phonon and its intensity can be calculated. Assuming amplitude \(Q_{k}\) and phase \(0 \leq \Phi_{k} \leq 1\) of the displacement wave, the displacements \(U(n,m)\) of atoms \((n,\mu)\) for a given wave vector \(k\) and phonon branch \(j\) are given by the equation

(30)\[\begin{split}U(n,\mu) = \dfrac{Q_{k}}{2 \sqrt{M_{\mu}}} \left[ \begin{matrix} \textrm{Re} e (k,j;\mu) cos \left[ 2 \pi \left( k \cdot R(n,\mu) - \Phi_{k} \right) \right] - \\ \textrm{Im} e (k,j;\mu) sin \left[ 2 \pi \left( k \cdot R(n,\mu) - \Phi_{k} \right) \right] \\ \end{matrix} \right]\end{split}\]

The intensity of phonon modes is obtained from the form factors. The form factor projected on the wave vector is defined as

(31)\[F^{(p)}(k,j) = \dfrac{1}{k^{2}} \left( \sum_{\mu k} \dfrac{k \cdot e (k,j;\mu)}{\sqrt{M_{\mu}}} \right)^{2}\]

However, the intensity of a phonon mode is represented by the simple form factors

(32)\[F^{(s)}(k,j) = \dfrac{1}{k^{2}} \left( \sum_{\mu} \dfrac{e (k,j;\mu)}{\sqrt{M_{\mu}}} \right)^{2}\]

which may be applied to remove unessential phonon branches originating from backfolding, or to estimate relative intensities of all modes in varying Brillouin zones.

3.10.4.6. Phonon Density of States

The phonon density of states \(g(\omega)\) provides the frequency distribution of normal modes and is obtained as a histogram plot from

(33)\[\begin{split}g(\omega) = \dfrac{1}{nd \Delta \omega} \sum_{k,j} \left( \omega - \omega(k,j) \right) \textrm{ where } \delta_{\Delta \omega)}(x) = \biggr\lbrace \begin{matrix} & 1 & \frac{\Delta \omega}{2} < x \leq \frac{\Delta \omega}{2} \\ & 0 & otherwise \end{matrix}\end{split}\]

The summation covers all \(n\) wave vectors homogeneously distributed in the first Brillouin zone and all phonon branches. The frequency interval of the histogram is denoted \(\Delta \omega\) and \(d\) is the dimension of the dynamical matrix (equal to the number of phonon branches and the number of degrees of freedom in the unit cell). The phonon density of states is normalized

(34)\[\int g(\omega) d\omega = 1\]

The partial phonon density of states provides the contribution of one specific atom \(m\) vibrating along a particular Cartesian direction

(35)\[g_{i,\mu}(\omega) = \dfrac{1}{nd\Delta \omega} \sum \vert e_{i}(k,j;\mu) \vert ^{2} \delta_{\Delta \omega} \left( \omega - \omega(k,j) \right)\]

satisfying the normalization condition (\(d\) being the number of degrees of freedom, again).

(36)\[\int g_{i,\mu} (\omega) d \omega = \dfrac{1}{d}\]

The off-diagonal partial phonon density of state is defined as

(37)\[g_{il,\mu} (\omega) = \dfrac{1}{n d \Delta \omega} \sum_{k,j} e_{i} (k,j;\mu) e_{l}^{\ast} (k,j;\mu) \delta_{\Delta \omega} \left( \omega - \omega(k,j) \right)\]

where \(e_{i}(k,j;\mu)\) and \(e_{l}(k,j;\mu)\) are the ith and jth Cartesian components of the polarization vectors for particle \(m\) and phonon frequency \(\omega(k,j)\). The diagonal elements are the partial phonon density of states.

3.10.5. Thermodynamic Functions

The contribution of lattice vibrations to thermodynamic functions such as internal energy, free energy, entropy, heat capacity, and thermal displacements are derived from the integrated phonon density of states. For each thermodynamic function, the contributions of particular atoms and polarization directions can be separated making use of partial phonon densities of states. This section provides a brief summary of the equations involved. Since thermodynamic functions are sensitive to phonons of low frequencies, it is advised to check for correct low and high-temperature limits. These limits are also summarized below.

3.10.5.1. Internal Energy

The internal energy of a crystal of N unit cells is given by \(E_{tot} = NE\) where E is the internal energy of the unit cell in the harmonic approximation

(38)\[E = \sum E_{i,\mu}\]
(39)\[E_{i,\mu} = \dfrac{1}{2} d \int_{0}^{\infty} g_{i,\mu} (\omega) \hbar \omega coth \dfrac{\hbar \omega}{2 k_{B} T} d \omega\]

Herein, \(E_{i,\mu}\) is the contribution of atom \(\mu\) and direction \(i\) to the internal energy, \(g_{l,u}(w)\) is the corresponding partial phonon density of states, \(d\) is the number of degrees of freedom in the unit cell, \(\hbar\) is the Planck constant, \(k_{B}\) is the Boltzmann constant and \(T\) is the temperature. The low and high-temperature limits are

(40)\[\lim_{T \rightarrow 0} E_{i,\mu} = \dfrac{1}{2} d \int_{0}^{\infty} g_{i,\mu} \hbar \omega d \omega\]
(41)\[\lim_{T \rightarrow \infty} E_{i,\mu} = k_{B} T\]

3.10.5.2. Free Energy

The Helmholtz free energy of a crystal of \(N\) unit cells at a given volume is given by \(A_{tot}=NA\), where \(A\) is the Helmholtz free energy of the unit cell in harmonic approximation

(42)\[\begin{split}A = \sum A_{i,\mu} \\ A_{i,\mu} = dk_{B} T \int_{0}^{\infty} g_{i,\mu}(\omega) ln \left( 2 sinh \dfrac{\hbar \omega}{2 k_{B}T} \right) d \omega\end{split}\]

\(A_{i,m}\) is the contribution of atom \(m\) and direction \(i\) to the free energy. The low-temperature limit of the free energy is equal to the mean energy of a given degree of freedom.

(43)\[\lim_{T \rightarrow 0} A_{i,\mu} = \dfrac{1}{2} d \int g_{i,\mu} (\omega) \hbar \omega d \omega\]

3.10.5.3. Entropy

The entropy of the crystal is \(S_{tot} = N S\) and the entropy of the unit cell \(S\) is described by

(44)\[S = \sum_{i,\mu} S_{i,\mu}\]
(45)\[S_{i,\mu} = dk_{B} \int_{0}^{\infty} g_{i,\mu} \left( \dfrac{\hbar \omega}{2 k_{B}T} \left( coth \left( \dfrac{\hbar \omega}{2 k_{B} T} \right) - 1 \right) - ln ( 1 - exp \left( \dfrac{\hbar \omega}{k_{B}T} \right) \right) d \omega\]

where \(S_{i,\mu}\) is the contribution of atom \(m\) and direction \(i\) to the entropy. The low-temperature limit of the entropy is

(46)\[\lim_{T \rightarrow 0} S_{i,\mu} = 0\]

3.10.5.4. Heat Capacity

Within the harmonic approximation the heat capacity at constant volume and at constant pressure are equal. The heat capacity of the crystal is \(C_{V,tot} = C_{P,tot} = N C\), and the heat capacity of the unit cell \(C\) is represented in terms of contributions of atoms and directions \(C_{i,m}\) by

(47)\[C = \sum_{i,\mu} C_{i,\mu}\]
(48)\[C_{i,\mu} = dk_{B} \int g_{i,\mu} (\omega) \left( \dfrac{\hbar \omega}{2 k_{B} T} \right) ^{2} \dfrac{exp \left( \dfrac{\hbar \omega}{k_{B} T} \right) }{ \left( exp \left( \dfrac{\hbar \omega}{k_{B}T} \right) - 1 \right) ^{2}} d \omega\]

The low and high-temperature limits are

(49)\[\lim_{T \rightarrow 0} C_{i,\mu} = 0\]
(50)\[\lim_{T \rightarrow \infty} C_{i,\mu} = k_{B}\]

3.10.5.5. Thermal Displacements, Debye-Waller Factor

Thermal vibrations have a considerable effect on neutron scattering data. The form factor describing diffraction scattering contains the Debye-Waller factor \(exp[W_{\mu}(k)]\), where the argument \(W_{\mu}(k)\)

(51)\[W_{\mu}(k) = \dfrac{1}{2} (2 \pi k) \cdot B(\mu)(2 \pi k)\]

is described in term of the 3x3 matrix \(B(\mu)\) representing the static correlation function of displacements \(U(\mu)\) of atom \(m\) from the equilibrium position

(52)\[B_{ij}(\mu) = \left( U_{i}(\mu)U_{j}(\mu) \right)\]

The symmetric matrix \(B(\mu)\) represents the square mean displacements of an atom \(m\), and is determined by the off-diagonal partial phonon density of states

(53)\[B_{il}(\mu) = \dfrac{\hbar d}{2 M_{\mu}} \int_{0}^{\infty} g_{il,\mu} (\omega) \dfrac{1}{\omega} coth \dfrac{\hbar \omega}{2 k_{B}T} d\omega\]

Debye-Waller factors are the main ingredient for analyzing cross sections for phonon creation in neutron scattering experiments.

3.10.6. Neutron Scattering

There are coherent and incoherent contributions to neutron scattering depending on the properties of the nuclei involved. To obtain sensible results, for each constituent nucleus some neutron scattering input data are required: the coherent scattering length \(a_{coh}(\mu)\), the incoherent scattering cross section \(\sigma_{incoh}(\mu)\), and the total scattering cross section for polycrystals \(\sigma_{tot}(\mu)\). The main equations for obtaining neutron scattering cross sections are summarized below.

3.10.6.1. Coherent Neutron Scattering

The doubly differential coherent scattering cross section for the creation of one phonon is given by the expression

(54)\[\dfrac{d^{2} \sigma_{coh}^{(+)}}{d \Omega dE} = \dfrac{K}{K_{0}} \sum_{j} \int \dfrac{\hbar F(k,j)}{2 \omega (k,j)} \left( n \left( \omega(k,j)+1 \right) \times (\kappa - k) \times (\epsilon - \omega(k,j)) \right) d \omega\]

where the coherent form factor, also called dynamical structure factor, is defined as

(55)\[F(k,j) = \biggr\vert \sum \alpha_{coh}(\mu) exp (-W_{\mu (k)}) \dfrac{2 \pi k \cdot e(k,j;\mu)}{\sqrt{M_{\mu}}} \biggr\vert ^{2}\]

and

(56)\[n(\omega(k,j)+1) = \dfrac{exp \frac{\hbar \omega}{k_{B}T}}{exp \frac{\hbar \omega}{k_{B}T} -1} \quad \epsilon = E_{0} - E \quad \kappa = K_{0} - K\]

Here, \(exp[W_{\mu}(k)]\) denotes the Debye-Waller factor (see section Entropy and \(a_{coh}(m)\) the coherent scattering length. Neutron scattering spectra are sensitive to the direction of incidental and scattered neutron beams and orientation of the monocrystal. This is taken into account in the equation above by the difference \(\epsilon\) of energy of incident neutrons \(E_{0}\) and scattered neutrons \(E\), and by the neutron momentum change \(\kappa = K_{0} - K\), where \(K_{0}\) indicates the direction of the incident neutron beam and \(K\) the direction of scattered neutrons. Both types of neutron scattering experiments, direct geometry (fixed energy of incident neutrons and spectrum given as a function of the energy of scattered neutrons) and inverse geometry (energy of scattered neutrons fixed and spectrum is shown as a function of energy of incident neutrons) may be simulated according to these equations.

3.10.6.2. Incoherent Neutron Scattering on Monocrystals

The one-Phonon double differential incoherent neutron scattering cross section is calculated as a sum of contributions from all atoms of the unit cell

(57)\[\dfrac{\partial ^{2} \sigma_{incoh}}{\partial \Omega \partial E} = \sum \dfrac{\partial ^{2} \sigma_{incoh}}{\partial \Omega \partial E} (\mu)\]

The contribution of one atom \(\mu\) then is

(58)\[\dfrac{\partial ^{2} \sigma_{incoh}}{\partial \Omega \partial E} = \sigma_{incoh} (\omega) \dfrac{K}{K_{0}} \sum_{k,j} \dfrac{\hbar}{2 M_{\mu} \omega (k,j)} (\kappa \cdot e(k,j;\mu))^{2} \cdot n(\omega(k,j)+1) e^{-2W_{\mu}(k)} \delta (\epsilon - \omega(k,j))\]

\(\sigma_{incoh}(\mu)\) is a materials constant and denotes the incoherent scattering cross section for each atom \({\mu}\), and all other quantities in this equation are explained in the previous section. In terms of the off-diagonal partial phonon density of states this equation reads

(59)\[\dfrac{\partial ^{2} \sigma_{incoh}}{\partial \Omega \partial E} = \sigma_{incoh} (\mu) \dfrac{K}{K_{0}} \sum_{k,j} d \dfrac{\hbar}{2M\mu} \int n(\omega(k,j)+1) \cdot e^{-2W_{\mu}}(\kappa) \sum_{i,j=1}^{3} \kappa_{i} \kappa_{l} g_{il,\mu} (\omega) \delta (\epsilon - \omega) d \omega\]

3.10.6.3. Incoherent Neutron Scattering on Polycrystalline Materials

For analyzing neutron scattering on polycrystalline materials contributions from both, incoherent and coherent scattering are considered. Since the coherent contribution disregards the wave vector conversation law due to orientational averaging, it is treated as effective incoherent scattering. For polycrystals, the double differential scattering cross section can be expressed for one-Phonon as well as for multi-Phonon processes.

The n-Phonon orientationally averaged double differential neutron scattering cross section is written as a sum of contributions from all atoms in the unit cell

(60)\[\dfrac{\partial ^{2} \sigma_{incoh}^{(av)} (n)}{\partial \Omega \partial E} = \sum \dfrac{\partial ^{2} \sigma_{incoh}^{(av)}}{\partial \Omega \partial E}\]

The contribution of atom \(\mu\) is

(61)\[\begin{split}\dfrac{\partial ^{2} \sigma_{incoh}^{(av)} (n,\mu)}{\partial \Omega \partial E} = \sigma_{tot}(\mu) \dfrac{K}{K_{0}} exp \left( -W_{\mu}^{(av)} (\omega) \right) \dfrac{1}{n!} \left( \dfrac{\hbar^{2} \omega^{2}}{2M_{\mu}} \right)^{n} \cdot \\ \quad \cdot \int_{-\infty}^{\infty} f_{\mu}(\omega_{1}) d \omega_{1} \int_{-\infty}^{\infty} f_{\mu}(\omega_{2}) d \omega_{2}... \int_{-\infty}^{\infty} f_{\mu}(\omega_{n}) d \omega_{n} \delta (\epsilon - \omega_{1} - \omega_{2} - ... - \omega_{n})\end{split}\]

where

(62)\[W_{\mu}^{(av)} (\kappa) = \dfrac{1}{6} \kappa^{2} T r \left( B(\mu) \right)\]

and the orientationally averaged Debye-Waller factor is

(63)\[f_{\mu}(\omega) = \dfrac{d exp \frac{\hbar \omega}{k_{B}T}}{3 \omega \left( exp \frac{\hbar \omega}{k_{B}T} - 1 \right)} \left( g_{x,\mu}(\omega) + g_{y,\mu}(\omega) + g_{z,\mu}(\omega) \right)\]

3.10.7. LO/TO Splitting

Polar crystals exhibit splitting of infrared active optical modes into longitudinal and transversal modes at the \(\Gamma\)-point (so-called LO/TO split). The frequency of LO modes is raised above those of the TO modes as a consequence of the long-range part of the Coulomb interaction, i.e. the macroscopic electric field arising from the displacements of entire ionic sublattices. This field breaks the Born-von K{‘a}rm{‘a}` n periodic boundary conditions and leads to a non-analytical term of the dynamical matrix at the \(\Gamma\)-point. As a consequence, without any further approximation the LO/TO mode splitting at the \(\Gamma\)-point cannot be treated with the direct method and only the TO mode only is obtained.

Applying elongated supercells the phonon frequency at \(G\) will still be degenerate, however, it is possible to recover the \(k \rightarrow 0\) limit of the LO branch by extrapolating the LO curve along the elongated direction from \(k \pi 0 \textrm{ to } k=0\).

As an alternative to this procedure, e.g. in case the elongation of the cell would raise the system size beyond computational feasibility for the ab initio calculations, one may take into account the non-analytical contribution of the dynamical matrix in an approximate form [19]. Introducing phenomenological charges in the form of the Born effective charge tensor \(Z^{\ast}(\mu)\) and the dielectric constant \(\epsilon_{\infty}\), this approximate additional term is written as

(64)\[\begin{split}D_{\alpha,\beta}^{M}(k;\mu,\nu) = D_{\alpha,\beta}(k;\mu,\nu) + 4 \pi \dfrac{e^{2}}{V \epsilon_{\infty} \sqrt{M_{\alpha} M_{\beta}}} \dfrac{(k \cdot Z^{\ast}(\mu))_{\alpha} (k \cdot Z^{\ast} (\mu))_{\beta}}{|k|^{2}} \times \\ \times e^{-2 \pi i g \cdot \left( r(\mu) - r(\mu) \right)} \times d (q) e^{\left( -\pi^{2} \left( \left( \frac{k_{x}}{\rho_{x}} \right)^{2} + \left( \frac{k_{y}}{\rho_{y}} \right)^{2} + \left( \frac{k_{z}}{\rho_{z}} \right)^{2} \right) \right) }\end{split}\]

where the damping factor \(d(q)\) is

(65)\[\begin{split}d(q) = \Bigg\lbrace \begin{matrix} & \dfrac{1}{2} + cos \left[ \pi \left( \dfrac{\sqrt{q_{1}^{2}+ q_{2}^{2} + q_{3}^{2}}}{\sqrt{b_{1}^{2}+ b_{2}^{2} + b_{3}^{2}}} \right)^{n} \right] \quad & \textrm{ if n } \geq \textrm{ 1 } \\ & \dfrac{1}{2} + cos \left[ \pi \left( 1- \dfrac{\sqrt{q_{1}^{2}+ q_{2}^{2} + q_{3}^{2}}}{\sqrt{b_{1}^{2}+ b_{2}^{2} + b_{3}^{2}}} \right)^{n} \right] \quad & \textrm{if 0 < n < 1} \\ \end{matrix}\end{split}\]

In the non-analytical term the wave vector \(k=(k_{x},k_{y},k_{z})\) is counted from the closest Brillouin zone center given by the vector \(g\), \(q=(q_{1},q_{2},q_{3})\) is the lattice vector in reciprocal lattice coordinates, and \(B=(b_{1},b_{2},b_{3})\) is the lattice vector from the Brillouin zone center to the Brillouin zone surface in the direction specified by \(q\). The volume of the primitive cell is denoted \(V\); \(M_{i}\) and \(r_{i}\) are the atomic masses and positions within the primitive cell. The index \(n\) is the power of the interpolation function. It allows modeling of the longitudinal phonon dispersion curve between the Brillouin zone center and the Brillouin zone surface. For \(n \geq 1\) the longitudinal dispersion curve is closer to the value of at the Brillouin zone center for most wave vectors, except in the close vicinity of the Brillouin zone surface. The opposite behavior is obtained for \(0<n<1\), where the dispersion curve reaches the longitudinal phonon mode only quite close to the Brillouin zone center.

The vector \(\rho_{i} = \rho \kappa_{i}\) is determined by the scalar quantity \(\rho\) and the vector \(\kappa_{i}\), the wave vector distance from the Brillouin zone center to the Brillouin zone surface. The macroscopic electric field range factor \(\rho\) is a free parameter that can further suppress the influence of the second term in Eq. (57). The fixed default value is 10.0, which eliminates the Gaussian damping factor from suppressing dispersion curves. The electric charge conservation law must be satisfied, and that can be achieved either by setting appropriate parameters of the Born effective charge tensor, or the program may take care of this by manipulating the charge tensor.

[1](1, 2) Max Born and K Huang, Dynamical Theory of Crystal Lattices, Clarendon Press, Oxford, 1966.
[2]B Dick and A Overhauser, “Theory of the Dielectric Constants of Alkali Halide Crystals,” Physical Review 112, no. 1 (October 1, 1958): 90; W Cochran, “Theory of the Lattice Vibrations of Germanium,” Physical Review Letters 2, no. 12 (June 15, 1959): 495; U Schröder, “A New Model for Lattice Dynamics (“Breathing Shell Model”),” Solid State Communications 4, no. 7 (1966): 347-349.
[3]P Keating, “Effect of Invariance Requirements on the Elastic Strain Energy of Crystals with Application to the Diamond Structure,” Physical Review 145, no. 2 (May 13, 1966): 637.
[4]Richard Martin, “Dielectric Screening Model for Lattice Vibrations of Diamond-Structure Crystals,” Physical Review 186, no. 3 (October 15, 1969): 871. Werner Weber, “Adiabatic Bond Charge Model for the Phonons in Diamond, Si, Ge, and A-Sn,” Physical Review B 15, no. 10 (May 15, 1977): 4789.
[5]Robert Pick, Morrel Cohen, and Richard Martin, “Microscopic Theory of Force Constants in the Adiabatic Approximation,” Physical Review B 1, no. 2 (January 15, 1970): 910.
[6]Pierre Hohenberg and Walter Kohn, “Inhomogeneous Electron Gas,” Physical Review B 136, no. 3 (1964): 864-871.
[7]A Rahman, “Correlations in the Motion of Atoms in Liquid Argon,” Physical Review 136, no. 2 (October 19, 1964): A405.
[8]Glenn J Martyna, Michael L Klein, and Mark Tuckerman, “Nos\({\acute{e}}\) -Hoover Chains: the Canonical Ensemble via Continuous Dynamics,” Journal of Physical Chemistry 97, no. 4 (1992): 2635.
[9]S Lawrence Marple, “Digital Spectral Analysis with Applications,” Englewood Cliffs, NJ, Prentice-Hall (1987). Jorge Kohanoff, Wanda Andreoni, and Michele Parrinello, “Zero-Point-Motion Effects on the Structure of C 6 0,” Physical Review B 46, no. 7 (August 15, 1992): 4371. Jorge Kohanoff, “Phonon Spectra From Short Non-Thermally Equilibrated Molecular Dynamics Simulations,” Computational Materials Science 2, no. 2 (1994): 221-232.
[10]X Gonze and J-P Vigneron, “Density-Functional Approach to Nonlinear-Response Coefficients of Solids,” Physical Review B 39, no. 18 (June 15, 1989): 13120.
[11]Stefano Baroni, Paolo Giannozzi, and Andrea Testa, “Green’s-Function Approach to Linear Response in Solids,” Physical Review Letters 58, no. 18 (May 4, 1987): 1861. Paolo Giannozzi et al., “Ab Initio Calculation of Phonon Dispersions in Semiconductors,” Physical Review B 43, no. 9 (March 15, 1991): 7231.
[12]Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi, “Phonons and Related Crystal Properties From Density-Functional Perturbation Theory,” Reviews of Modern Physics 73, no. 2 (July 6, 2001): 515.
[13]M Yin and Marvin L Cohen, “Theory of Lattice-Dynamical Properties of Solids: Application to Si and Ge,” Physical Review B 26, no. 6 (September 15, 1982): 3259.
[14]K Kunc and Richard Martin, “Ab Initio Force Constants of GaAs: a New Approach to Calculation of Phonons and Dielectric Properties,” Physical Review Letters 48, no. 6 (February 8, 1982): 406.
[15]S Wei and M Chou, “Ab Initio Calculation of Force Constants and Full Phonon Dispersions,” Physical Review Letters (1992). Siqing Wei and M Chou, “Phonon Dispersions of Silicon and Germanium From First-Principles Calculations,” Physical Review B 50, no. 4 (July 15, 1994): 2221.
[16]Georg Kresse, J Furthm ller, and J rgen Hafner, “Ab Initio Force Constant Approach to Phonon Dispersion Relations of Diamond and Graphite,” Europhysics Letters (EPL), 1995. W Frank, Christian Elsässer, and M Fähnle, “Ab Initio Force-Constant Method for Phonon Dispersions in Alkali Metals,” Physical Review Letters 74, no. 10 (March 6, 1995): 1791. K Parlinski, ZQ Li, and Yoshiyuki Kawazoe, “First-Principles Determination of the Soft Mode in Cubic ZrO 2,” Physical Review Letters 78, no. 21 (1997): 4063-4066.
[17]Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi, “Phonons and Related Crystal Properties From Density-Functional Perturbation Theory,” Reviews of Modern Physics 73, no. 2 (July 6, 2001): 515. G J Ackland, M C Warren, and S J Clark, “Practical Methods in Ab Initio Lattice Dynamics,” Journal of Physics: Condensed Matter, 1997. Pasquale Pavone, “Old and New Aspects in Lattice-Dynamical Theory,” Journal of Physics: Condensed Matter, 2001. K Parlinski, ZQ Li, and Yoshiyuki Kawazoe, “First-Principles Determination of the Soft Mode in Cubic ZrO 2,” Physical Review Letters 78, no. 21 (1997): 4063-4066.
[18]WH Press, SA Teukolsky, WT Vetterling, and BP Flannery, Numerical Recipes in FORTRAN Example Book: the Art of Scientific Computing, Cambridge University Press, 1992.
[19]Peter E Blöchl, “Projector Augmented-Wave Method,” Physical Review B 50, no. 24 (December 1994): 17953-17979.
download:pdf