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.
Fig. 6.1.1 Silicon supercell \(3\times 3\times 3\) reconstruction.
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.hzw file (
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.input file (
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.
Fig. 6.1.2 The Hessian matrix module simulator.
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.deltaThis 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.orderThe 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.
Fig. 6.1.3 Phonon band structure of a silicon supercell of \(3\times 3\times 3\) reconstruction.
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.