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

1. Import the crystal;

2. Expand the primitive cell;

3. Self-consistent calculation;

4. Hessian matrix calculation;

5. 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)

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:

1. In the menu bar of Device Studio, click SimulatorNanoDCALSCF 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 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.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.

1. 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 SimulatorNanoDCALAnalysis to pop up the interface shown in Fig. 6.2.4.

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

3. 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:

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

2. 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 SimulatorNanoDCALAnalysis 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.

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