Phonons: Theory

Phonons are the collective excitation of nuclei in an extended periodic system.

Taylor expansion in ionic displacements

To compute them we start by Taylor expanding the total energy ($E$) around the set of equilibrium positions of the nuclei ($\{\mathbf {R} ^{0}\}$)

$E(\{\mathbf {R} \})=E(\{\mathbf {R} ^{0}\})+\sum _{I\alpha }{\frac {\partial E(\{\mathbf {R^{0}} \})}{\partial R_{I\alpha }}}(R_{I\alpha }-R_{I\alpha }^{0})+\sum _{I\alpha J\beta }{\frac {\partial E(\{\mathbf {R} ^{0}\})}{\partial R_{I\alpha }\partial R_{J\beta }}}(R_{I\alpha }-R_{I\alpha }^{0})(R_{J\beta }-R_{J\beta }^{0})+{\mathcal {O}}(\mathbf {R} ^{3}),$

where $\{\mathbf {R} \}$ the positions of the nuclei. The first derivative of the total energy with respect to the nuclei corresponds to the forces

$F_{I\alpha }(\{\mathbf {R} ^{0}\})=-\left.{\frac {\partial E(\{\mathbf {R} \})}{\partial R_{I\alpha }}}\right|_{\mathbf {R} =\mathbf {R^{0}} }$,

and the second derivative to the second-order force-constants

$\Phi _{I\alpha J\beta }(\{\mathbf {R} ^{0}\})=\left.{\frac {\partial E(\{\mathbf {R} \})}{\partial R_{I\alpha }\partial R_{J\beta }}}\right|_{\mathbf {R} =\mathbf {R^{0}} }=-\left.{\frac {\partial F_{I\alpha }(\{\mathbf {R} \})}{\partial R_{J\beta }}}\right|_{\mathbf {R} =\mathbf {R^{0}} }.$

Changing variables in the Taylor expansion of the total energy with $u_{I\alpha }=R_{I\alpha }-R_{I\alpha }^{0}$ that corresponds to the displacement of the atoms with respect to their equilibrium position $R_{I\alpha }^{0}$ leads to

$E(\{\mathbf {R} \})=E(\{\mathbf {R} ^{0}\})+\sum _{I\alpha }-F_{I\alpha }(\{\mathbf {R} ^{0}\})u_{I\alpha }+\sum _{I\alpha J\beta }\Phi _{I\alpha J\beta }(\{\mathbf {R} ^{0}\})u_{I\alpha }u_{J\beta }+{\mathcal {O}}(\mathbf {R} ^{3})$

Dynamical matrix and phonon modes

If we take the system to be in equilibrium, the forces on the atoms are zero and then the Hamiltonian of the system is

$H={\frac {1}{2}}\sum _{I\alpha }M_{I}{\dot {u}}_{I\alpha }^{2}+{\frac {1}{2}}\sum _{I\alpha J\beta }\Phi _{I\alpha J\beta }u_{I\alpha }u_{J\beta },$

with $M_{I}$ the mass of the $I$-th nucleus. The equation of motion is then given by

$M_{I}{\ddot {u}}_{I\alpha }=-\Phi _{I\alpha J\beta }u_{J\beta }.$

We then look for solutions of the form of plane waves traveling parallel to the wave vector $\mathbf {q}$, i.e.

$\mathbf {u} _{I\alpha }^{\mu }(\mathbf {q} ,t)={\frac {1}{2}}{\frac {1}{\sqrt {M_{I}}}}\left\{A^{\mu }(\mathbf {q} )\xi _{I\alpha }^{\mu }(\mathbf {q} )e^{i[\mathbf {q} \cdot \mathbf {\mathbf {R} } _{I}-\omega _{\mu }(\mathbf {q} )t]}+A^{\mu *}(\mathbf {q} )\xi _{I\alpha }^{\mu *}(\mathbf {q} )e^{-i[\mathbf {q} \cdot \mathbf {\mathbf {R} } _{I}-\omega _{\mu }(\mathbf {q} )t]}\right\}$

where $\xi _{I\alpha }^{\mu }(\mathbf {q} )$ are the phonon mode eigenvectors and $A^{\mu }(\mathbf {q} )$ the amplitudes. Replacing it in the equation of motion we obtain the following eigenvalue problem

$\sum _{J\beta }D_{I\alpha J\beta }(\mathbf {q} )\xi _{J\beta }^{\mu }(\mathbf {q} )=\omega ^{\mu }(\mathbf {q} )^{2}\xi _{I\alpha }^{\mu }(\mathbf {q} )$

with

$D_{I\alpha J\beta }(\mathbf {q} )={\frac {1}{\sqrt {M_{I}M_{J}}}}\Phi _{I\alpha J\beta }e^{i\mathbf {q} \cdot (\mathbf {R} _{J}-\mathbf {R} _{I})}$

the dynamical matrix in the harmonic approximation. Now by solving the eigenvalue problem above we can obtain the phonon modes $\xi _{I\alpha }^{\mu }(\mathbf {q} )$ and frequencies $\omega ^{\mu }(\mathbf {q} )$ at any arbitrary q point.

We can write the positions of the atoms in the supercell $\mathbf {R} _{I}$ in terms of integer multiples of the lattice vectors of the unit cell $\mathbf {l}$ such that $\mathbf {R} _{I}\rightarrow \mathbf {R} _{li}=\mathbf {l} +\mathbf {r} _{i}$ with $\mathbf {r} _{i}$ being the position of the ion in the unit cell. The force constants then become $\Phi _{I\alpha ,J\beta }\rightarrow \Phi _{li\alpha ,l'j\beta }$. The dynamical matrix is then given by

$D_{i\alpha j\beta }(\mathbf {q} )={\frac {1}{\sqrt {M_{i}M_{j}}}}\sum _{l}\Phi _{li\alpha ,0j\beta }e^{-i\mathbf {q} \cdot (\mathbf {l} +\mathbf {r} _{i}-\mathbf {r} _{j})},$

with $\mathbf {l}$ chosen using the minimal image convention.

This allows us to compute the phonons in the unit cell using the following equation

$\sum _{j\beta }D_{i\alpha j\beta }(\mathbf {q} )\xi _{j\beta }^{\mu }(\mathbf {q} )=\omega ^{\mu }(\mathbf {q} )^{2}\xi _{i\alpha }^{\mu }(\mathbf {q} )$

with the dynamical matrix having dimension $3n$ with $n$ the number of atoms in the unit cell.

Long-range interatomic force constants (LO-TO splitting)

For semiconductors or insulators, the electronic screening of the ions is incomplete which leads to long-range (LR) interatomic force constants. To compute them explicitly would require infinitely large supercell calculations. For practical calculations, a finite size truncation is needed which leads to Gibbs oscillations in the phonon dispersion. Fortunately, this long-range behavior can be modeled by looking at the analytic form of the ion-ion contribution to the total energy.

For that, follow the approach outlined in Ref. and start by splitting the second-order force constants into short-range and long-range parts,

$\Phi _{I\alpha J\beta }=\Phi _{I\alpha J\beta }^{\text{SR}}+\Phi _{I\alpha J\beta }^{\text{LR}}$

with the long-range part being obtained from the analytic derivative of the long-range part of the ion-ion contribution to the total energy ${\textstyle E_{\text{ion-ion}}}$. This contribution is typically evaluated using an Ewald sum technique in which we separate the ion-ion contribution to the total energy into two part, one is evaluated in real space and captures the short-range part and the other one in reciprocal space which captures the long-range part ${\textstyle E_{\text{ion-ion}}=E_{\text{ion-ion}}^{\text{SR}}+E_{\text{ion-ion}}^{\text{LR}}}$. The separation is governed by a Ewald parameter ${\textstyle \alpha }$ which represent a truncation length.

This leads to the following analytical expression for the long-range interatomic force constants,

$\Phi _{I\alpha J\beta }^{\text{LR}}={\frac {4\pi e^{2}}{\Omega _{0}}}\sum _{\mathbf {G} }{\frac {(\mathbf {G} \cdot \mathbf {Z} _{I\alpha }^{*})(\mathbf {G} \cdot \mathbf {Z} _{J\beta }^{*})}{\mathbf {G} \cdot \epsilon ^{\infty }\cdot \mathbf {G} }}e^{i\mathbf {G} \cdot (\mathbf {R} _{J}-\mathbf {R} _{I})}e^{\frac {-\mathbf {G} \cdot \epsilon ^{\infty }\cdot \mathbf {G} }{4\alpha ^{2}}}$

with ${\textstyle \epsilon ^{\infty }}$ the clamped ion dielectric function, ${\textstyle \mathbf {Z} _{I\alpha }^{*}}$ the Born effective charges, ${\textstyle \alpha }$ the Ewald parameter which is chosen such that the contributions from ${\textstyle e^{\frac {-\mathbf {G} \cdot \epsilon ^{\infty }\cdot \mathbf {G} }{4\alpha ^{2}}}}$ are negligible within a certain ${\textstyle \mathbf {G} }$ vector cutoff sphere PHON_G_CUTOFF.

This also allows us to separate the dynamical matrix into short and long-range parts

$D_{I\alpha J\beta }(\mathbf {q} )=D_{I\alpha J\beta }^{\text{SR}}(\mathbf {q} )+D_{I\alpha J\beta }^{\text{LR}}(\mathbf {q} ),$

with the long-range part of the dynamical matrix

$D_{i\alpha j\beta }^{\text{LR}}(\mathbf {q} )={\frac {4\pi e^{2}}{\Omega _{0}}}\sum _{\mathbf {G} }\sum _{l}{\frac {{\big [}(\mathbf {G} +\mathbf {q} )\cdot \mathbf {Z} _{i\alpha }^{*}{\big ]}{\big [}(\mathbf {G} +\mathbf {q} )\cdot \mathbf {Z} _{j\beta }^{*}{\big ]}}{(\mathbf {G} +\mathbf {q} )\cdot \epsilon ^{\infty }\cdot (\mathbf {G} +\mathbf {q} )}}e^{i(\mathbf {q} +\mathbf {G} )\cdot (\mathbf {l} +\mathbf {r} _{i}-\mathbf {r} _{j})}e^{\frac {-(\mathbf {G} +\mathbf {q} )\cdot \epsilon ^{\infty }\cdot (\mathbf {G} +\mathbf {q} )}{4\alpha ^{2}}}$

The equations above give us the practical method for computing the phonon dynamical matrices including the long-range force constants using a moderately sized supercell calculation with the steps:

• Compute ${\textstyle \Phi _{I\alpha J\beta }}$ using a finite size supercell
• Compute ${\textstyle \Phi _{I\alpha J\beta }^{\text{SR}}=\Phi _{I\alpha J\beta }-\Phi _{I\alpha J\beta }^{\text{LR}}}$
• Compute ${\textstyle D_{i\alpha i\beta }^{\text{SR}}(\mathbf {q} )}$ using ${\textstyle \Phi _{I\alpha J\beta }^{\text{SR}}}$
• Compute Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle D^\text{LR}_{i\alpha i\beta}(\mathbf{q})} in the unit cell and add to $D_{i\alpha i\beta }^{\text{SR}}(\mathbf {q} )$

The treatment is done automatically inside VASP using the LPHON_POLAR=.TRUE. tag and specifying the dielectric tensor with PHON_DIELECTRIC and the Born effective charges with PHON_BORN_CHARGES.

Finite differences

The second-order force constants are computed using finite differences of the forces when each ion is displaced in each independent direction. This is done by creating systems with finite ionic displacement of atom $a$ in direction $i$ with magnitude $\lambda$, computing the orbitals $\psi _{\lambda }^{u_{I\alpha }}$ and the forces for these systems. The second-order force constants are then computed using

$\Phi _{I\alpha J\beta }={\frac {\partial ^{2}E}{\partial u_{I\alpha }\partial u_{J\beta }}}=-{\frac {\partial F_{I\alpha }}{\partial u_{J\beta }}}\approx -{\frac {\left(\mathbf {F} [\{\psi _{\lambda /2}^{u_{J\beta }}\}]-\mathbf {F} [\{\psi _{-\lambda /2}^{u_{J\beta }}\}]\right)_{I\alpha }}{\lambda }},\quad {I=1,..,N_{\text{atoms}}}\quad {J=1,..,N_{\text{atoms}}}\quad {\alpha =x,y,z}\quad {\beta =x,y,z}$

where $u_{I\alpha }$ corresponds to the displacement of atom $I$ in the cartesian direction $\alpha$ and $\mathbf {F} [\psi ]$ retrieves the set of forces acting on all the ions given the $\psi _{n\mathbf {k} }$ KS orbitals.

Similarly, the internal strain tensor is

$\Xi _{I\alpha l}={\frac {\partial ^{2}E}{\partial u_{I\alpha }\partial \eta _{l}}}={\frac {\partial \sigma _{l}}{\partial u_{I\alpha }}}\approx {\frac {\left(\mathbf {\sigma } [\{\psi _{\lambda /2}^{u_{I\alpha }}\}]-\mathbf {\sigma } [\{\psi _{-\lambda /2}^{u_{I\alpha }}\}]\right)_{l}}{\lambda }},\qquad {l=xx,yy,zz,xy,yz,zx}\quad {\alpha =x,y,z}\quad {J=1,..,N_{\text{atoms}}}$

where $\mathbf {\sigma } [\psi _{n\mathbf {k} }]$ computes the strain tensor given the $\psi _{n\mathbf {k} }$ KS orbitals.

Density functional perturbation theory

Density-functional-perturbation theory provides a way to compute the second-order linear response to ionic displacement, strain, and electric fields. The equations are derived as follows.

In density-functional theory, we solve the Kohn-Sham (KS) equations

$H(\mathbf {k} )|\psi _{n\mathbf {k} }\rangle =e_{n\mathbf {k} }S(\mathbf {k} )|\psi _{n\mathbf {k} }\rangle ,$

where $H(\mathbf {k} )$ is the DFT Hamiltonian, $S(\mathbf {k} )$ is the overlap operator and, $|\psi _{n\mathbf {k} }\rangle$ and $e_{n\mathbf {k} }$ are the KS eigenstates.

Taking the derivative with respect to the ionic displacements $u_{I\alpha }$, we obtain the Sternheimer equations

$\left[H(\mathbf {k} )-e_{n\mathbf {k} }S(\mathbf {k} )\right]|\partial _{u_{I\alpha }}\psi _{n\mathbf {k} }\rangle =-\partial _{u_{I\alpha }}\left[H(\mathbf {k} )-e_{n\mathbf {k} }S(\mathbf {k} )\right]|\psi _{n\mathbf {k} }\rangle$

Once the derivative of the KS orbitals is computed, we can write

$|\psi _{\lambda }^{u_{I\alpha }}\rangle =|\psi \rangle +\lambda |\partial _{u_{I\alpha }}\psi \rangle .$

where $\lambda$ is a small numeric value to use in the finite differences formulas below.

The second-order response to ionic displacement, i.e., the force constants or Hessian matrix can be computed using the same equation used in the case of the finite differences approach

$\Phi _{I\alpha J\beta }={\frac {\partial ^{2}E}{\partial u_{I\alpha }\partial u_{J\beta }}}=-{\frac {\partial F_{I\alpha }}{\partial u_{J\beta }}}\approx -{\frac {\left(\mathbf {F} [\{\psi _{\lambda /2}^{u_{J\beta }}\}]-\mathbf {F} [\{\psi _{-\lambda /2}^{u_{J\beta }}\}]\right)_{I\alpha }}{\lambda }},\quad {I=1,..,N_{\text{atoms}}}\quad {J=1,..,N_{\text{atoms}}}\quad {\alpha =x,y,z}\quad {\beta =x,y,z}$

where again $\mathbf {F} [\{\psi \}]$ yields the forces for a given set of $\psi _{n\mathbf {k} }$ KS orbitals.

Similarly, the internal strain tensor is computed using

$\Xi _{I\alpha l}={\frac {\partial ^{2}E}{\partial u_{I\alpha }\partial \eta _{l}}}={\frac {\partial \sigma _{l}}{\partial u_{I\alpha }}}\approx {\frac {\left(\mathbf {\sigma } [\{\psi _{\lambda /2}^{u_{I\alpha }}\}]-\mathbf {\sigma } [\{\psi _{-\lambda /2}^{u_{I\alpha }}\}]\right)_{l}}{\lambda }},\qquad {l=xx,yy,zz,xy,yz,zx}\quad {\alpha =x,y,z}\quad {J=1,..,N_{\text{atoms}}}$

where $\mathbf {\sigma } [\psi _{n\mathbf {k} }]$ computes the strain tensor given the $\psi _{n\mathbf {k} }$ KS orbitals. The Born effective charges are then computed using Eq. (42) of Ref. .

$Z_{I\alpha \gamma }^{*}=2{\frac {\Omega _{0}}{(2\pi )^{3}}}\int _{\text{BZ}}\sum _{n}\langle \partial _{u_{I\alpha }}\psi _{n\mathbf {k} }|({\vec {\beta }}_{n\mathbf {k} })_{\gamma }\rangle d\mathbf {k}$

where $I$ is the atom index, $\alpha$ the direction of the displacement of the atom, $\gamma$ the polarization direction, and $|{\vec {\beta }}_{n\mathbf {k} }\rangle$ is the polarization vector defined in Eq. (30) in Ref. . The results should be equivalent to the ones obtained using LCALCEPS and LEPSILON.