# 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.

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:

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)

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.

```
[1]:
```

```
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
```

```
[1]:
```

```
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.

```
[2]:
```

```
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.

```
[3]:
```

```
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.

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]\).

```
[4]:
```

```
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.

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]\).

```
[5]:
```

```
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:

```
[6]:
```

```
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)
```