Phonon

1. Introduction

1.1 What is phonon?

Essentially, a phonon is the vibration mode of atoms in the crystal. Before we go to the concepts, let’s look at what is a phonon from this website.

Animation visualization of phonon. Cite from http://henriquemiranda.github.io/phononwebsite/phonon.html.

Let’s open the website -> select Si in the left-upper materials column -> click any point in right figure, i.e. phonon dispersion. Then, you can see every atom (yellow ball) in the middle window shows periodic motion, and the collection movement of all atoms seems like a wave that travels from one end to another.

Such a vibration mode in the crystal is called lattice vibration. Lattice vibration usually deals with quantum effects and it is called “Phonon”. In this tutorial explanation, we don’t distinguish between the term lattice vibration and phonon. Please refer appendix & reference books for more detail.

The direction of wave propagation is usually represented by a vector \(\vec q\) (called wave vector). The most important problem in the phonon analysis is to know what kind of vibration modes are allowed for a given wave vector \(\vec q\), and the frequency of atomic vibration in those vibration modes.

Phonon is a kind of quasiparticle. It has energy (\(E\)) and momentum (\(\mathbf{p}\)). The phonon in a solid is directly related to its thermal and acoustic properties. For example, the excitation of a phonon causes the increase of internal energy (\(U\)); the transportation of the phonon causes the conduction of sound and heat. As a consequence, phonon analysis has been a very important simulation tool to obtain the properties of solids.

1.2 How to calculate the phonon frequency?

Similar to the vibration analysis (section 4-1), the phonon is also calculated with the harmonic approximation. However, the atom interacts with those in the neighboring unit cell in the crystal. By considering the periodical boundary condition, the potential energy of a small distorted crystal structure can be expressed as:

\[V(\mathbf{r}) \approx V(\mathbf{r_0}) +　\frac{1}{2} \sum_{ij} \frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^a \partial r_j^b} \Delta \mathbf{r}_i^a \Delta \mathbf{r}_j^b\]

where, \(r_i^a\) indicates the coordination of the \(i\)-th atom at the \(a\)-th supercell. The force constant (explained in section 4-1) of a crystal should be a \(3NM \times 3NM\) matrix, where \(N\) is the number of atoms per unit cell and \(M\) is the number of the unit cell to be considered. The force constant can be a very large matrix.

We assume that the atomic displacement from its equilibrium position \(u\) in a vibration mode in the crystal should have the following form (plane-wave-like function)

\[u(\vec R, t) = A e^{ \vec{q} \vec{R} -i\omega_{m} t}\]

Then, the large force constant matrix can be converted into a \(3N \times 3N\) dynamical matrix by discrete Fourier transform in real space. The dynamical matrix depends on the wave vector \(\vec q\). The eigenvalue and eigenvector of the dynamical matrix are the vibration frequencies and vibration modes of the phonon. For a detailed explanation of the phonon calculation formula, please refer to the appendix.

In the next section, we show how to calculate the phonon properties using the matlantis-features package. Some important concepts of phonon are introduced at the same time.

2. Example: phonon calculation of Si

2.1 Material preparation

Here, we use crystal silicon as an example to run the phonon calculation. As the first step, we should relax the crystal structure. Usually, we use the primitive cell of the crystal to calculate phonon dispersion since it contains the least number of atoms.

import numpy as np
np.set_printoptions(precision=3)

from ase.build import bulk
from ase.optimize import BFGS
from ase.constraints import ExpCellFilter

import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode

estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v3.0.0")
calculator = ASECalculator(estimator)
atoms = bulk("Si")
atoms.calc = calculator

exp = ExpCellFilter(atoms)
BFGS(exp).run(fmax=0.01)

      Step     Time          Energy         fmax
BFGS:    0 07:07:45       -9.108037        0.4967
BFGS:    1 07:07:45       -9.113260        0.0002

True

2.2 Calculation of force constant matrix

Same with the vibration analysis, the force constant is necessary for calculating the vibration frequencies of phonons. Here, we obtain the force constant of the crystal using the ForceConstantFeature.

The primitive cell of crystal Si contains only 2 atoms. We should repeat the unit cell many times to obtain an accurate result. It is because the force constant element \(\frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^a \partial r_j^b}\) will decay with the distance between atoms \(i\) of unit cell \(a\) and the atoms \(j\) of unit cell \(b\). Only if the supercell is large enough, all non-negligible elements of the force constant matrix can be included. In this example, we expand the crystal Si structure by 10x10x10 time.

The ForceConstantFeature in matlantis-features is used to compute the force constant of crystal silicon.

from matlantis_features.features.phonon import ForceConstantFeature

fc = ForceConstantFeature(
    supercell = (10,10,10),
    delta = 0.1,
)
force_constant = fc(atoms)
print(force_constant.force_constant.shape)

(1000, 6, 6)

The calculation result, i.e. force constant \(C\), is a 1000x6x6 tensor.

print(force_constant.force_constant[1])

[[-0.196 -0.195 -0.115 -0.044 -0.026  0.031]
 [-0.195 -0.196 -0.115 -0.026 -0.044  0.031]
 [ 0.115  0.115  0.43   0.031  0.031 -0.136]
 [-3.125 -2.088  2.088 -0.196 -0.195  0.115]
 [-2.088 -3.125  2.088 -0.195 -0.196  0.115]
 [ 2.088  2.088 -3.125 -0.115 -0.115  0.43 ]]

The force_constant.force_constant[1] is a 6x6 tensor. It is the force constant between the first and second unit cell \(C^{1,2}\), i.e. $:nbsphinx-math:frac{partial^2 V(mathbf{r_0})}{partial r_i^1 partial r_j^2} $. Same as the vibration analysis, the element of \(C^{1,2}\) is approximated with the finite differences method.

\[C^{1,2}[i,j]=\frac{\partial^2 V(\mathbf{r_0})}{\partial r_i^1 \partial r_j^2} \approx \frac{F(\mathbf{r_0} + \Delta r_i^1)_j^2 - F(\mathbf{r_0})_j^2}{|\Delta{r_i^1}|}\]

Since we used 10x10x10 unit cells, we should have 1000x1000 such a 6x6 matrix. However, the crystal has translational invariance, which makes the \(C^{a,b} = C^{0, b-a}\). Thus, we can represent the force constant with a 1000x6x6 tensor.

2.3 Calculation of dynamical matrix

Unlike the molecular vibration, in which the vibration frequency is obtained from diagonalization of force constant \(C\), the vibration frequency of the crystal is the eigenvalue of the dynamical matrix \(D_{\vec q}\), which is the discrete Fourier transform of the force constant in the real space. To obtain the dynamical matrix, we introduced another parameter, wave vector \(\vec q\).

Then, we use the PhononFreqency provided by matlantis-features to obtain the dynamical matrix of crystal silicon for the wave vector \(\vec q = [0.1, 0.0, 0.0]\).

from matlantis_features.features.phonon.utils import PhononFrequency

phonon_frequency = PhononFrequency(
    force_constant.force_constant, force_constant.supercell, force_constant.unit_cell_atoms
)
dynamic_matrix = phonon_frequency._dynamic_matrix_q(np.array([0.1, 0.0, 0.0])).real
print("Dynamical matrix of q=[0.1, 0.0, 0.0]:")
print(dynamic_matrix)

Dynamical matrix of q=[0.1, 0.0, 0.0]:
[[ 0.456 -0.001 -0.001 -0.433 -0.013 -0.013]
 [-0.001  0.456  0.001 -0.013 -0.433  0.013]
 [-0.001  0.001  0.456 -0.013  0.013 -0.433]
 [-0.433 -0.013 -0.013  0.456 -0.001 -0.001]
 [-0.013 -0.433  0.013 -0.001  0.456  0.001]
 [-0.013  0.013 -0.433 -0.001  0.001  0.456]]

The dynamical matrix is a symmetric 6x6 matrix. It is calculated from the force constant \(C\) with the following formula.

\[\tilde{D}_{\vec{q}}[k\alpha, k'\beta] = \frac{1}{\sqrt{M_k M_{k'}}} \tilde{C}(\vec{q})[k\alpha, k'\beta]\]

\[\tilde{C}_{\vec{q}}[k\alpha, k'\beta] = \sum_{b} {\left( \frac{\partial ^{2} E}{\partial{u_{k}^{a\alpha}} \partial{u_{k'}^{b\beta}}} \right) e^{i\vec{q}(\vec{R_b} - \vec{R_a})}} = \sum_{b} {C^{0,b}[k\alpha, k'\beta] e^{i\vec{q}(\vec{R_b} - \vec{R_a})}}\]

Then, we can diagonalize this matrix and obtain the eigenvalue and eigenvector. Same as molecular vibration, the eigenvalue is the vibration frequency of the phonon, and the eigenvector is related to the displacement of atoms in the phonon mode. Please see the appendix for the deduction of above equations.

So, let’s get the phonon frequency of crystal silicon at \(\vec q = [0.1, 0.0, 0.0]\).

frequency, eigenvector = phonon_frequency._frequency_mode_q(np.array([0.1, 0.0, 0.0]), return_mode=True)
for i in range(6):
    print(f"Band {i}")
    print(f"    Frequency: {frequency[i]:.3f} meV")
    print(f"    Eigenvector:")
    print(eigenvector[i].real)

Band 0
    Frequency: 5.910 meV
    Eigenvector:
[[-0.075  0.028 -0.103]
 [-0.074  0.029 -0.103]]
Band 1
    Frequency: 5.912 meV
    Eigenvector:
[[ 0.052  0.087 -0.035]
 [ 0.048  0.083 -0.035]]
Band 2
    Frequency: 11.188 meV
    Eigenvector:
[[ 0.073 -0.073 -0.073]
 [ 0.077 -0.077 -0.077]]
Band 3
    Frequency: 60.955 meV
    Eigenvector:
[[-0.073  0.073  0.073]
 [ 0.077 -0.077 -0.077]]
Band 4
    Frequency: 61.389 meV
    Eigenvector:
[[-0.09  -0.014 -0.077]
 [ 0.088  0.011  0.077]]
Band 5
    Frequency: 61.390 meV
    Eigenvector:
[[-0.013 -0.09   0.077]
 [ 0.01   0.088 -0.077]]

wave number :math:`vec q`

The wave vector \(\vec q\) is introduced into the phonon calculation. The physical meaning of wave vector is the propagation direction of vibration mode. Theoretically, the \(\vec q\) is limited to the first Brillouin zone (please see the appendix for the detailed explaination).

Let’s use the function get_brillouin_zone_3d to see the first brillouin zone of crystal sillicon. Please compare with the figure in Wikipedia:

First Brillouin zone of FCC lattice, showing symmetry labels for high symmetry lines and points. Cite from https://en.wikipedia.org/wiki/Brillouin_zone

from matlantis_features.utils.visual_utils.brillouin_zone import get_brillouin_zone_3d

fig = get_brillouin_zone_3d(atoms.cell, show_reciprocal_basis=True, show_special_k_points=True)
fig.update_layout(width=600, height=500)