# 6.1. Frozen phonon calculation

There are two main methods for phonon calculations: the density functional perturbation theory (DFPT) method and the frozen phonon method. The DFPT method is based on linear response theory; the frozen phonon method is based on calculating the force constant matrix (Hessian) by explicitly displacing atoms in real space. This tutorial is for the frozen phonon method in the NanoDCAL package. For DFPT, the user is referred to the RESCU package.

In this tutorial, we’ll show how to calculate the *Phonon band structure* in a silicon supercell (Fig. 6.1.1) using NanoDCAL code.

The calculation steps are given as follows:

Import the crystal;

Expand the primitive cell;

Self-consistent calculation;

Hessian matrix calculation;

Phonon band structure calculation;

## 6.1.1. Self-consistent calculation

In the previous tutorial (SOI), we already learned how:

Import the Si crystal using Device Studio software (step 1);

How to expand the crystal cell by

Redefine Crystal(step 2).

Note

Note that for crystals with small unit cells, it is often necessary to accurately account for long-range interactions in the dynamical matrix using periodically repeating supercells of the unit cell. Therefore it is important to set the number of repetitions along each lattice vector direction to at least 3 times the number of primitive cells.

In this exemple we’ll expanded the primitive cell 3 times. Thus, in the cell vector set **u = (3a, 0b, 0c)**, **v = (0a, 3b, 0c)**, and **w = (0a, 0b, 3c)**

The atomic coordenates for Si primitive cell can be download here:

Silicon primitive cell

Si.hzwfile (`Si-primitivecell.pbs`

).

Si.hzw

```
0
0.0000000000 2.7153500000 2.7153500000
2.7153500000 0.0000000000 2.7153500000
2.7153500000 2.7153500000 0.0000000000
2
Si 0.6788375000 0.6788375000 0.6788375000
Si 2.0365125000 2.0365125000 2.0365125000
```

From the step 2, we can generate the input file to peform the self-consistent calculation for the silicon supercell. The input file could be generated using Device Studio with NanoDCAL module simulator, as follows:

In the menu bar of Device Studio, click

**Simulator**➟**NanoDCAL**➟**SCF Calculation**to pop up the interface to set the parameters;On the

**Exchange correlation**choose**PBE_GGA96**;In the

**K-point Sampling**set**n1 = 5**,**n2 = 5**, and**n3 = 5**;

The generated input file for silicon-supercell could be downloaded here:

Silicon supercell

scf.inputfile (`Si-supercell.scf.input`

).

The keywords were also specified in the previous tutorial. It is can be check in keywords-scf.

Scf input file

```
%%What quantities should be calculated
calculation.name = scf
%Basic setting
calculation.occupationFunction.temperature = 300
calculation.realspacegrids.E_cutoff = 80 Hartree
calculation.xcFunctional.Type = GGA_PBE96
calculation.k_spacegrids.number = [ 5 5 5 ]'
system.centralCellVectors = [[0 8.14605 8.14605]' [8.14605 0 8.14605]' [8.14605 8.14605 0]']
system.spinType = NoSpin
%Iteration control
calculation.SCF.monitoredVariableName = {'rhoMatrix','hMatrix','totalEnergy','bandEnergy','gridCharge','orbitalCharge'}
calculation.SCF.convergenceCriteria = {1e-04,1e-04,[],[],[],[]}
calculation.SCF.maximumSteps = 200
calculation.SCF.mixMethod = Pulay
calculation.SCF.mixRate = 0.1
calculation.SCF.mixingMode = H
calculation.SCF.startingMode = H
%calculation.SCF.donatorObject = NanodcalObject.mat
%Basic set
system.neutralAtomDataDirectory = '../'
system.atomBlock = 54
AtomType OrbitalType X Y Z
Si PBE-DZP 0.67883750 0.67883750 0.67883750
Si PBE-DZP 3.39418750 3.39418750 0.67883750
Si PBE-DZP 6.10953750 6.10953750 0.67883750
Si PBE-DZP 3.39418750 0.67883750 3.39418750
Si PBE-DZP 6.10953750 3.39418750 3.39418750
Si PBE-DZP 8.82488750 6.10953750 3.39418750
Si PBE-DZP 6.10953750 0.67883750 6.10953750
Si PBE-DZP 8.82488750 3.39418750 6.10953750
Si PBE-DZP 11.54023750 6.10953750 6.10953750
Si PBE-DZP 0.67883750 3.39418750 3.39418750
Si PBE-DZP 3.39418750 6.10953750 3.39418750
Si PBE-DZP 6.10953750 8.82488750 3.39418750
Si PBE-DZP 3.39418750 3.39418750 6.10953750
Si PBE-DZP 6.10953750 6.10953750 6.10953750
Si PBE-DZP 8.82488750 8.82488750 6.10953750
Si PBE-DZP 6.10953750 3.39418750 8.82488750
Si PBE-DZP 8.82488750 6.10953750 8.82488750
Si PBE-DZP 11.54023750 8.82488750 8.82488750
Si PBE-DZP 0.67883750 6.10953750 6.10953750
Si PBE-DZP 3.39418750 8.82488750 6.10953750
Si PBE-DZP 6.10953750 11.54023750 6.10953750
Si PBE-DZP 3.39418750 6.10953750 8.82488750
Si PBE-DZP 6.10953750 8.82488750 8.82488750
Si PBE-DZP 8.82488750 11.54023750 8.82488750
Si PBE-DZP 6.10953750 6.10953750 11.54023750
Si PBE-DZP 8.82488750 8.82488750 11.54023750
Si PBE-DZP 11.54023750 11.54023750 11.54023750
Si PBE-DZP 2.03651250 2.03651250 2.03651250
Si PBE-DZP 4.75186250 4.75186250 2.03651250
Si PBE-DZP 7.46721250 7.46721250 2.03651250
Si PBE-DZP 4.75186250 2.03651250 4.75186250
Si PBE-DZP 7.46721250 4.75186250 4.75186250
Si PBE-DZP 10.18256250 7.46721250 4.75186250
Si PBE-DZP 7.46721250 2.03651250 7.46721250
Si PBE-DZP 10.18256250 4.75186250 7.46721250
Si PBE-DZP 12.89791250 7.46721250 7.46721250
Si PBE-DZP 2.03651250 4.75186250 4.75186250
Si PBE-DZP 4.75186250 7.46721250 4.75186250
Si PBE-DZP 7.46721250 10.18256250 4.75186250
Si PBE-DZP 4.75186250 4.75186250 7.46721250
Si PBE-DZP 7.46721250 7.46721250 7.46721250
Si PBE-DZP 10.18256250 10.18256250 7.46721250
Si PBE-DZP 7.46721250 4.75186250 10.18256250
Si PBE-DZP 10.18256250 7.46721250 10.18256250
Si PBE-DZP 12.89791250 10.18256250 10.18256250
Si PBE-DZP 2.03651250 7.46721250 7.46721250
Si PBE-DZP 4.75186250 10.18256250 7.46721250
Si PBE-DZP 7.46721250 12.89791250 7.46721250
Si PBE-DZP 4.75186250 7.46721250 10.18256250
Si PBE-DZP 7.46721250 10.18256250 10.18256250
Si PBE-DZP 10.18256250 12.89791250 10.18256250
Si PBE-DZP 7.46721250 7.46721250 12.89791250
Si PBE-DZP 10.18256250 10.18256250 12.89791250
Si PBE-DZP 12.89791250 12.89791250 12.89791250
end
```

*Run* the scf calculation. The output *NanodcalObject.mat* file is required for the Hessian matrix calculation.

## 6.1.2. Hessian matrix calculation

The Hessian matrix is calculated using a finite difference scheme, where each unique atom of the symmetry is displaced along all cartesian directions. It is a well-known approach called frozen phonons.

Frozen phonons are an alternative way to calculate phonons. The monochromatic perturbation is frozen with a finite amplitude in the system. The Fourier transform of force constants at the wave vector is calculated from finite differences of forces induced on all the atoms of the supercell by the monochromatic perturbation.

Once the Hamiltonian from the supercell are calculated (*NanodcalObject.mat*), we can move to Hessian matrix calculation.

The required input file for this procedure could be generated using Device Studio with NanoDCAL simulator module, as follows:

In the menu bar of Device Studio, click on

**Simulator**➟**NanoDCAL**➟**Analysis**to pop up the interface shown in Fig. 6.2.4.

Navigate to the

**Analysis**panel (left side) ➟ select**Hessian**➟ click on the*Right arrow*button to transfer**Hessian**matrix calculator to the**Calculation Selected**panel.Set the values according to the requirements of the atomic structure and click

**Generate files**.

In this tutorial, the default values are enough to describe the structure.
However, the user should note that more acurate calculations can be peformed by NanoDCAL code. Please, look at the documentation *Input Reference*.
Now, we are ready to start the Hessian matrix calculation.
Before doing that, we’ll briefly present the simulation parameters specified in the generated files.

Hessian matrix input file

```
system.object = NanodcalObject.mat
calculation.name = hessian
calculation.hessian.delta = 0.02
calculation.hessian.order = 2
calculation.hessian.primitiveCellVectors = [[ 0 8.14605 8.14605]' [8.14605 0 8.14605]' [8.14605 8.14605 0]']
calculation.control.xml = true
```

`calculation.hessian.delta`

This input parameter specifies that a small displacement of atomic position will be used to calculate the hessian matrix by numerical derivatives of forces. The default value is 0.03 au;`calculation.hessian.order`

The number of points used to calculate the derivative of force numerically. The possible values are 2, 3, or 5. The default value is 3 au.

*Run* the simulation with `hessian-matrix.input`

file.
After the calculation finish, the following output files are generated: *CalculatedResults.mat*, *log.txt*, and *Hessian.mat*.

The Hessian matrix is used to calculate the phonon frequencies at different points in the Brillouin zone. After this step we’ll calculate the phonon band structure, which provides information about how phonon frequencies vary with wave vector, and provides insights into the material’s vibrational properties.

## 6.1.3. Phonon band structure calculation

For the phonon band structure calculation the user should repeat the procedure discribed in here as following:

As shown in Fig. 6.2.4, navigate to the

**Analysis**panel (left side) ➟ select**PhononBandStructure**➟ click on the*Right arrow*button to transfer**PhononBandStructure**calculator to the**Calculation Selected**panel.The

**Symmetry k-points**and**Number of k-points**are show on the right side panel, click on**path**to see/edit the K-paths along the high-symmetry points in the Brillouin zone and press**Apply**. Returning to the**Analysis**window, click**Generate files**to create the*PhononBandStructure.input*file.

Phonon band structure input file:

```
system.object = NanodcalObject.mat
calculation.name = phononBandStructure
calculation.phononBandStructure.symmetryKPoints = {'U','X','G','X','W','K','G','L','U','W','L','K'}
%calculation.phononBandStructure.coordinatesOfTheSymmetryKPoints = [0.625 0.25 0.625;0.5 0 0.5;0 0 0;0.5 0 0.5;0.5 0.25 0.75;0.375 0.375 0.75;0 0 0;0.5 0.5 0.5;0.625 0.25 0.625;0.5 0.25 0.75;0.5 0.5 0.5;0.375 0.375 0.75]'
calculation.phononBandStructure.numberOfKPoints = 200
%calculation.phononBandStructure.plot = true
calculation.control.xml = true
```

The user should be able to *Run* the phonon band structure calculations with `phonon-bands.input`

file following the procedure previously described at run BS.

After the calculation is finished, the following output files are generated: *CalculatedResults.mat*, *log.txt*, *PhononBandStructure.mat*, *PhononBandStructure.xml*.
Those results could be plotted with MATLAB (*PhononBandStructure.mat*) or loaded in Device Studio software:

In menu bar of Device Studio, click on

**Simulator**➟**NanoDCAL**➟**Analysis Plot**➟ import*PhononBandStructure.xml*file.

Phonon band structure analysis

The phonon spectrum for the silicon supercell is shown in Fig. 6.1.3.

The key features and interpretations of the phonon band structure of a silicon supercell include:

This system has no imaginary frequency at \(\Gamma\)-point (G). In this case we can demonstrate that the Si structure is stable.

The horizontal axis of the band structure plot represents the crystal’s reciprocal lattice vectors (typically denoted as k-points), while the vertical axis represents the phonon frequencies. Different branches or curves correspond to different vibrational modes, and their slopes provide information about the speed or group velocity of the phonons.

Supercells may exhibit zone-folding effects in their phonon band structure. These effects arise due to the periodicity of the supercell, which affects the allowed phonon wave vectors and leads to folding of bands into the first Brillouin zone.

In summary, the phonon band structure of a silicon supercell provides insights into its vibrational properties, helping useres to understand thermal conductivity, lattice dynamics, and other related phenomena in the material.

Sergueev, D. Roubtsov, and Hong Guo Ab Initio Analysis of Electron-Phonon Coupling in Molecular Devices Phys. Rev. Lett 95 (2005) p. 146803.