Preparation

Preparation before phonon calculations

Based on the derivation of Phonon formalism below, the optimization of structure must be done before the beginning of phonon calculations (Phonon calculations: Pre-process). The structure optimization (i.e., optimizing the positions of constituent atoms into equilibrium) is carried out by relaxation calculation. It is generally recommended to use a strict force criterion for the relaxation convergence (e.g., corresponding to EDIFFG tag in VASP). This is because highly accurate values of forces acting on each atom are required for phonon evaluation.

Phonon formalism

In crystals, the equilibrium position of an atom is described by the labels of lattice (\(l\)) and basis (\(s\)):

\[R_{ls} = R_{l}+\tau_{s}\]

where \(R_{l}\) indicates a specific unit cell among infinite lattices and \(\tau_{s}\) is the position of the atom in the unit cell. In real situation where vibrational motions around their equilibrium positions {\(R_{ls}\)} with displacements {\({u_{ls}}\)} take place, the potential energy \(U\) of a crystal can be written by a Taylor expansion around the equilibrium positions up to the second order (harmonic approximation):

\[U = U_{0}+\frac{1}{2} \sum_{ls\alpha}\sum_{l's'\beta} \frac{\partial^{2}U}{\partial u_{ls\alpha}\partial u_{l's'\beta}} u_{ls\alpha}u_{l's'\beta}\]

where \(U_0\) is the zeroth term corresponding to the potential energy at equilibrium; and \(\alpha\) and \(\beta\) indicate the Cartesian directions. Note that the first derivative of the potential energy (corresponding to force) vanishes, since the force on each atom is zero near the equilibrium position. The second derivative, called the second order force constants (\(K_{ls\alpha,l's'\beta}\)), is numerically evaluated using the midpoint method within the FDM scheme:

\[K_{ls\alpha,l's'\beta} = \frac{\partial^{2}U}{\partial u_{ls\alpha}\partial u_{l's'\beta}} = -\frac{\partial F_{ls\alpha}}{\partial u_{l's'\beta}} = -\frac{F_{ls\alpha}(R_{l's'\beta}+u) - (R_{l's'\beta}-u)}{2u}\]

where \(F_{ls\alpha}(R_{l's'\beta}+u)\) is the force on the atom at \(R_{ls}\) along \(\alpha\) direction in response to the displacement of the atom at \(R_{l's'}\) by \(u\) along \(\beta\) direction.

The periodicity of crystals imposes the displacement to take the form of plane wave (\(u_{ls\alpha} = u^0_{s\alpha} e^{i(qR_{ls}-wt)}\)) as for a sound wave in solids. Finally, by making several transformations, the eigenvalue equation of the dynamical matrix with elements of \(D_{s\alpha,s'\beta}(q)\) is obtained:

\[w^2 v^0_{s\alpha} = \sum_{s'\beta} D_{s\alpha,s'\beta}(q) v^0_{s'\beta},\]

\[D_{s\alpha,s'\beta}(q) = \frac{1}{(M_s M_s')^{1/2}} \sum_{l'} K_{ls\alpha,l's'\beta} e^{iq(R_{l'}-R_{l})} e^{iq(\tau_{s'}-\tau_{s})}\]

where \(w\) is the square root of the eigenvalue corresponding to the normal frequency of lattice vibrations, and a vector with components \(v^0_{s\alpha}\) (\(=M^{1/2}_s u^0_{s\alpha}\)) is the eigenvector and is referred to as the mass-weighted normal mode.