Phonon dispersion relation

Phonon Dispersion Relation

Yanming Wang, Saad Bhamla and Wei Cai

June, 2012

---

METHODS

Small displacement method for calculating the phonon dispersion relation

For a crystal lattice composed of a number of atoms bound by a specific potential, an equilibrium or minimum energy state can be reached by relaxing the structure. This is achieved using the conjugate gradient relaxation method in MD++. Once the local minimum is reached, a Taylor expansion is used around this state in terms of the atomic displacements which gives Equation 1

$V(\{{\vec {r}}_{i}\})=V(\{{\vec {r}}_{i}\})+{\frac {1}{2}}\sum _{ij}{\vec {u}}_{i}^{T}\cdot [{\frac {d^{2}V}{d{\vec {r}}_{i}{\vec {r}}_{j}}}]\cdot {\vec {u}}_{j}\qquad (1)$

where $V$ is the potential function, ${\vec {r}}_{i}$ is the coordinate of the $ith$ atom, ${\vec {R}}_{i}$ is the position of the $ith$ atom at emuilibrium and ${\vec {u}}$ is a small displacement from the equilibrium position: ${\vec {u}}_{i}={\vec {r}}_{i}-{\vec {R}}_{i}$.

In Equation 1, the linear terms in the Taylor expansion are at 0K with the minimum energy state and the higher order terms are neglected using the Harmonic approximation. The second derivative of the potential energy evaluated at the equilibrium position and is called the force constant matrix or the Hessian matrix and is obtained from MD++ using the calHessian function.

$K_{ij}=[{\frac {d^{2}V}{d{\vec {r}}_{i}{\vec {r}}_{j}}}]=\left[{\begin{array}{ccc}{\frac {\partial ^{2}V}{\partial x_{i}\partial x_{j}}}&{\frac {\partial ^{2}V}{\partial x_{i}\partial y_{j}}}&{\frac {\partial ^{2}V}{\partial x_{i}\partial z_{j}}}\\{\frac {\partial ^{2}V}{\partial y_{i}\partial x_{j}}}&{\frac {\partial ^{2}V}{\partial y_{i}\partial y_{j}}}&{\frac {\partial ^{2}V}{\partial y_{i}\partial z_{j}}}\\{\frac {\partial ^{2}V}{\partial z_{i}\partial x_{j}}}&{\frac {\partial ^{2}V}{\partial z_{i}\partial y_{j}}}&{\frac {\partial ^{2}V}{\partial z_{i}\partial z_{j}}}\\\end{array}}\right]\qquad (2)$

From Equation 2, we can calculate the force on the atoms, as shown in Equation 3 as follows, ${\vec {f}}_{i}=-{\frac {dV}{d{\vec {r}}_{i}}}{\bigg |}_{{\vec {R}}_{i}}=-\sum _{j}K_{ij}{\vec {u}}_{j}\qquad (3)$

This force acting on the atoms leads to vibrations in the lattice. Applying the Newtonian equation of motion, we obtain Equation 4,

${\vec {f}}_{i}=m{\ddot {u}}=-\sum _{j}K_{ij}{\vec {u}}_{j}\qquad (4)$

Obviously, equation 4 is a second order differential equation about ${\vec {u}}$,and ${\vec {u}}$ can be derived from a general solution of a plane wave shown in Equation 5,

${\vec {u}}({\vec {R}},t)={\vec {A}}\cdot exp[i({\vec {k}}\cdot {\vec {R}}-\omega t)]\qquad (5)$

where ${\vec {A}}$ determines the amplitude and direction of vibration, $\omega$ represents the vibration modes of the system, and ${\vec {k}}$ represents the wave vector. Substituting Equation 5 in Equation 4, we obtain the following eigenfunction expression:

${\vec {A}}\cdot \xi =Q({\vec {k)}}\cdot {\vec {A}}\qquad (6)$

where $\xi =m\omega ^{2}$ and $Q({\vec {k)}}$ is the dynamical matrix defined by:

$Q({\vec {k)}}=\sum _{j}K_{ij}\cdot exp[i({\vec {k}}\cdot {\vec {R}}-\omega t)]\qquad (7)$

Equation 6 is an eigenfunction where $\xi$ represents the eigenvalues, and ${\vec {A}}$ represents the eigenvectors of the dynamical matrix $Q({\vec {k}})$. We may notice each given ${\vec {k}}$ coressponds to a unique $Q({\vec {k)}}$ and by solving the eigenvalues of $Q({\vec {k)}}$, finally we can obtain the relation between ${\vec {k}}$ and $\omega$, which is the phonon dispersion relation.

TEST CASE

The tcl script Graphene_md.tcl can provide the atoms' configurations and the Hessian Matrix for graphene dispersion curve calculation. First, you need to compile the REBO potential in MD++ using the command like,

make rebo build=R SYS=mac

Please refer the details for compiling to Introduction to MD++.

Then the following command can be run to generate the data for the calculation.

bin/rebo_mac scripts/ME346/graphene_md.tcl 0 0

After the simulation, you can find two files named as perf.cn and hessian.out. The matlab script Graphene_phonon.m can extract data from the files and draw the phonon dispersion curve for graphene.

In the end, you will have the same figure shown below,

You may want to compare this result with DFT calculation. Here is the plot. The red curves are from this small displacement method using REBO potential and the black dots are the data from DFT calculation.

The three phonon dispersion branches, which originate from the $\Gamma$-point of the Brilliouin zone correspond to acoustic modes: an out-of plane mode (ZA), an in-plane transverse mode (TA), and in-plane longitudinal (LA), listed in order of increasing energy. The remaining three branches correspond to optical modes: one out-of plane mode (ZO), and two in-plane modes (TO) and (LO).