Asaf's Course Project

Computational Physics Course Project - Calculation of Phonon Related Properties in PbTe, Using DFT with Quantum Espresso

Image
PbTe lattice - image made with XCRYSDEN

In this project I explore some possibilities of calculating phonon related properties in PbTe.

I am interested in PbTe because it is a well known material for thermoelectric applications, and I've encountered it in another project done in Prof. Yaron Amouyal's group , in the Materials Science&Engineering faculty. Thermoelectrics are materials that, based on the Seebeck effect (or inversely the Peltier effect), convert a temperature gradient to an electrical potential. Using this effect, scientists and engineers hope to harness otherwise unused waste heat for helpful purposes (e.g. using the waste heat produced by a vehicle's motor to power the air conditioning).

The efficiency of thermoelectric materials is defined by a dimensionless number, zT, known as the electrical figure of merit: $$zT=\frac{\sigma S^{2}T}{\kappa}$$ where $\sigma$ is the electric conductivity, $S$ the Seebeck coefficient, $T$ is the temperature and $\kappa$ is the thermal conductivity. A good thermoelectric should ideally have high electrical conductivity combined with relatively low thermal conductivity. Thermal conductivity in semiconductors is dominated by the lattice (phononic) thermal conductivity (electronic transport being only secondary), and hence the interest in researching phonon related properties.

Quantum Espresso is a widely used and freely distributed DFT package (under public GNU license). The main QE package, for basic DFT calculations such as structure relaxation and total energy is called PWscf (Plane-Wave self consistent Field). The QE distribution also includes some subpackages for more specific applications or post processing, as described in the user guide. The sub-package I will use is called Phonon, which calculated vibrational properties using Density functional perturbation theory (DFPT).

QE is installed on Tamnun, the Technion's parallel computer. Small tasks can be run directly from the command line, but it recommended to run all tasks through the PBS queuing system, is described in the Manual section.

My project consists of 3 bash script to be run on Tamnun, and some auxillary files.

The 3 scripts are:

  1. mpi_PbTe-vcrelax.sh : Stucural relaxation for finding the correct lattice parameter.
  2. run_ph_PbTe.sh : The main part of the project - calculates total engery (self consistent DFT calcultions) and phonon frequencies (the dynamical matrix) on a 4*4*4 grid in the recipocal lattice.
  3. run_and_plot_dispersion_PbTe.sh - Manipulates the raw data, Fourier transforming the k-space dynamical matrix, to real-space IFC's (Interatomic Force Constant), interpolates back to dynamical matrices over a select path in the reciprocal lattice. Plots the resulting phonon dispersion and density of states (DOS) with gnuplot.


All the files are available, together with output for reference here: Asaf_Project.tar.gz. Please download the archive and unpack it's contents using the following command (in Linux)

tar -zxvf Asaf_Project.tar.gz

In order to run the files, you must either have a user on Tamnun, or you need to install Quantum Espresso on your PC (my scripts won't work on PC's though, you'll have to edit them a little).


This project is based on tutorials provided by the QE team (specifically the chapter of DFPT by Nicola Bonini), and a few articles:

  1. Tian, Z., Garg, J., Esfarjani, K., Shiga, T., Shiomi, J., & Chen, G. (2012). Phonon conduction in PbSe, PbTe, and PbTe_{1?x}Se_{x} from first-principles calculations. Physical Review B, 85(18), 184303. doi:10.1103/PhysRevB.85.184303
  2. An, J., Subedi, A., & Singh, D. J. (2008). Ab Initio Phonon Dispersions for PbTe. . Retrieved from http://arxiv.org/abs/0806.2727v2
  3. Baroni, S., De Gironcoli, S., Dal Corso, A., & Giannozzi, P. (2001). Phonons and related crystal properties from density-functional perturbation theory . Reviews of Modern Physics.
  4. Baroni, S., Density Functional Theory - a different elementary introduction' , Internazionale, S., & Avanzati, S. (2007).
  5. Zhang, Y., Ke, X., Chen, C., Yang, J., & Kent, P. (2009). Thermodynamic properties of PbTe, PbSe, and PbS: First-principles study. Physical Review B, 80(2), 024304. doi:10.1103/PhysRevB.80.024304
  6. Cochran, W., Cowley, R. a., Dolling, G., & Elcombe, M. M. (1966). The Crystal Dynamics of Lead Telluride. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 293(1435), 433451. doi:10.1098/rspa.1966.0182
  7. A. Kokalj, XCrySDen--a new program for displaying crystalline structures and electron densities, J. Mol. Graphics Modelling, 1999, 17, 176--179.

Acknowledgements

I'd like to thank Dr. Arik Landu for his assistance with using Quantum Espresso, and Dr. Joan Adler for letting me use her group's computational resources (queue Simphony_q on Tamnun).
Go to Top

Theory

I present her some (very) brief background on DFT and phonons. The part on DFT is just for demystifying what DFT is, while the part on phonons is more relevant to the results of this project. Reference 3 and 4 (which I heavily rely upon) have lot's more information for those who are interested.

Density Functional Theory


Density Functional Theory (DFT) is an approximate method for solving Schrodinger's equation for many electron problems.
A generic hamiltonian of a system of many interacting electrons in some external potential is (using units in which $\hbar=m=e=1$ for simplicity, and under the Born-Openheimer approximation (i.e. Ions are stationary). The index MB indicates many-body system): $$\hat{H}_{MB}=-\frac{1}{2}\sum_{i}\nabla_{r_{i}}^{2}+\frac{1}{2}\sum_{i\neq j}\frac{1}{\left|\vec{r}_{i}-\vec{r}_{j}\right|}+\sum_{i}V\left(\vec{r}_{i}\right)$$ An exact solution for with system would require to solve the Schrodinger equation for the ground state wave function and energy, $ \Psi\left(\vec{r}_{1},\ldots,\vec{r}_{N}\right)$ and $E^{0}$, respectively: $$\hat{H}_{MB}\Psi\left(\vec{r}_{1},\ldots,\vec{r}_{N}\right)=E^{0}\Psi\left(\vec{r}_{1},\ldots,\vec{r}_{N}\right)$$ But, the whole point of DFT is to find the ground state energy without solving the full many body equation. We want to somehow rewrite the problem in terms of the electron charge density (that's where the D in DFT comes from)

So naturally, we begin by defining the electron charge density: $$n\left(\vec{r}\right)=N\int\left|\Psi\left(\vec{r}_{1},\ldots,\vec{r}_{N}\right)\right|d\vec{r}_{2}\ldots d\vec{r}_{N}$$ Let's cell the kinetic energy operator $-\frac{1}{2}\sum_{i}\nabla_{r_{i}}^{2}\to\hat{K}$ , and the electron-electron interaction term $\frac{1}{2}\sum_{i\neq j}\frac{1}{\left|\vec{r}_{i}-\vec{r}_{j}\right|}\to\hat{W}$. Using this notation (and the definition of the electron charge density) we can state: $$E=\left\langle \Psi\left|\hat{H}\right|\Psi\right\rangle =\left\langle \Psi\left|\hat{K}+\hat{W}\right|\Psi\right\rangle +\int n\left(\vec{r}\right)V\left(\vec{r}\right)d\vec{r}$$
At this point we consult the basic theorem behind DFT: the Hohenberg and Kohn (HK - 1964) Theorem, which makes (together with some auxillary theorems) the following statements:
  1. No two different (external) potentials ($V\left(\vec{r}\right)$ ) acting on a given system of electrons can give rise to a same ground state electronic charge density, $n\left(\vec{r}\right)$: $$V\left(\vec{r}\right)\Longleftrightarrow n\left(\vec{r}\right)$$
  2. A universal functional (one that does not depend on the specific potential $V\left(\vec{r}\right)$) of the electron charge density $F\left[n\left(\vec{r}\right)\right]$ exists such that the energy functional $$E\left[n\right]=F\left[n\right]+\int n\left(\vec{r}\right)V\left(\vec{r}\right)d\vec{r}$$ is minimized by the electron charge density of the ground state. Using the variational principle, this can be alternatively stated in a differential form: $$\frac{\delta F}{\delta n\left(r\right)}+V\left(r\right)=\mu$$ where $\mu$ is a Lagrange multiplier.
  3. The value of this minimum for $E\left[n\right]$ coincides with the with the ground state energy.


Using the HK theorem, we can identify the term $\left\langle \Psi\left|\hat{K}+\hat{W}\right|\Psi\right\rangle $ in equation (5) as that 'universal functional' ,$F\left[n\right]$ , and we can rest assure that it exists. We don't, however, know its exact form!

This is where the machinery of finding the ground state charge density kicks in: the Kohn-Sham equations (1965). Kohn and Sham's strategy is to 'localize' our ignorance about the functional form of $F\left[n\right]$ by identifying in it all the sensible contributions, that are both physically motivated and easily calculable (in terms of the electron density). The leads to following equation for $F\left[n\right]$ : $$F\left[n\right]=T_{0}\left[n\right]+E_{H}\left[n\right]+E_{XC}\left[n\right]$$ $T_{0}$ denotes the kinetic energy of the electrons. $E_{H}=\frac{1}{2}\int\frac{n\left(\vec{r}\right)n\left(\vec{r}\right)}{\left|\vec{r}-\vec{r}'\right|}d\vec{r}d\vec{r}'$ is called the Hartree Energy and is the classical electrostatic coulombic contribution we would expect from a cloud of charge density. The last term, $E_{XC}$ , is called the exchange and correlation energy, and is really just all the rest: everything about the electron-electron interaction we don't know how to calculate in terms of $n\left(\vec{r}\right)$ . The actual implementation of DFT relies on making good approximations for this term. For now let's assume that such an approximation exists. Using the definition for $F\left[n\right]$ we made in equation (9), applying the variational principle to the energy functional $E\left[n\right]$ gives: $$\frac{\delta T_{0}}{\delta n\left(r\right)}+V\left(r\right)+V_{H}\left(r\right)+V_{XC}\left(r\right)=\mu$$ where $V_{H}$ and $V_{XC}$ are the functional derivatives (with respect to $n\left(r\right)$ ) of $E_{H}$ and $E_{XC}$ , respectively.
Now for the 'trick' of Kohn&Sham: if we define an effective potential (denoted $V_{KS}$ ), $$V_{KS}\left(\vec{r}\right)=V\left(r\right)+V_{H}\left(r\right)+V_{XC}\left(r\right)$$ equation (10) takes the form: $$\frac{\delta T_{0}}{\delta n\left(r\right)}+V_{KS}\left(r\right)=\mu$$ This is formally equivalent to equation (8), but for a system of non-interacting electrons (because of there is no electron-electron interaction, all that is left of $F\left[n\right]$ is the kinetic energy term). Since it is a non-interacting system, we can define an auxiliary (fictitious) Schrodinger equation for a single electron, that's relatively easy to solve: $$\left(-\frac{1}{2}\nabla^{2}+V_{KS}\left(\vec{r}\right)\right)\phi_{n}\left(\vec{r}\right)=\epsilon_{n}\phi_{n}\left(\vec{r}\right)$$ Using these fictitious single particle wave functions, $\phi_{n} \left(\vec{r}\right)$ (called the Kohn-Sham orbitals) we can recalculate the electron density of the multi-particle system (summation is over the total number of electrons in the system) $$n\left(\vec{r}\right)=\sum_{n}\left\langle \phi_{n}\left(\vec{r}\right)|\phi_{n}\left(\vec{r}\right)\right\rangle$$ Knowing the electron charge density we can now calculate the total energy of the system: $$E=\underbrace{-\frac{1}{2}\sum_{n}\int\phi_{n}^{\star}\frac{\partial^{2}\phi_{n}\left(\vec{r}\right)}{\partial\vec{r}^{2}}d\vec{r}}_{T_{0}\left[n\right]}+E_{H}\left[n\right]+E_{XC}\left[n\right]+\int V\left(\vec{r}\right)n\left(\vec{r}\right)$$ It is important to notice that while the result of the single particle Schrodinger equation (13) is the electron density, the effective potential $V_{KS}\left(\vec{r}\right)$ depends itself on $n\left(\vec{r}\right)$! This means that the Kohn Sham equations must be solved in an iterative or self-consistent way: repeating the calculation loop until the resulting charge density is unchanged (up to some convergence threshold).

All in all, DFT enormously simplifies the task of finding the ground state properties of a system of interacting electrons : instead of trying to find a wave function $\psi\left(\vec{r}_{1}\ldots\vec{r}_{N}\right)$ that depends on the 3N (N electrons times 3 spatial dimensions) variables, we seek a solution based on the electron charge density $n\left(\vec{r}\right)$s , which depends only on 3 spatial variables. In recent years, new and improved exchange and correlation functionals have made DFT computations accurate for many applications. The relative speed and accuracy of DFT has made it the de facto standard for modeling and simulating materials on the atomic scale.

Phonons



Phonons are the quantum mechanical descriptions of vibrational motion in which a lattice of atoms or molecules oscillate uniformly and at a single frequency. They are the quantum analog on normal modes. The lattice-dynamical properties of a system are defined by the eigenvalues and eigenfunction of the lattice hamiltonian (taking \hbar=1 for simplicity): $$\left(-\sum_{l}\frac{1}{M_{I}}\nabla_{\vec{R}_{I}}^{2}+E\left(\vec{R}\right)\right)\Phi\left(\vec{R}\right)=\varepsilon\Phi\left(\vec{R}\right)$$ $\vec{R}_{I}$ is the coordinate of the Ith nucleus and $M_{I}$ is its mass. $\vec{R}$ denotes the set of all nuclear coordinates($\vec{R}=\left\{ R_{I}\right\}$ ). $E\left(R\right)$ , called the Born-Oppenheimer energy, is the energy of a system of interacting electrons moving in the field of the fixed nuclei (i.e. the result of a typical self consistent DFT calculation). The (electronic) Hamiltonian by which $E\left(\vec{R}\right)$ is computed depends parametrically on $\vec{R}$ , meaning that the variables of the wave function of the electronic Hamiltonian are the electron positions only. In analogy to a harmonic oscillator, the equilibrium positions of the nuclei are given by the condition that the forces on the nuclei vanish: $$\vec{F}_{I}=-\frac{\partial E\left(\vec{R}\right)}{\partial\vec{R}_{I}}=0$$ And the vibrational frequencies of motion around these equilibrium positions are related to the second derivatives of the potential, known as the matrix of Interatomic Force Constants (IFC's) $$\frac{\partial^{2}E\left(\vec{R}\right)}{\partial\vec{R}_{I}\partial\vec{R}_{J}}$$ This is a matrix (of size 3, for spatial degrees of freedom, times the number of atoms in the cell), whose eigenvalues give the vibrational frequencies $\omega$ (scaled by the nuclear masses), by the equation: $$\det\left|\frac{1}{\sqrt{M_{I}M_{J}}}\frac{\partial^{2}E\left(\vec{R}\right)}{\partial\vec{R}_{I}\partial\vec{R}_{J}}-\omega^{2}\right|=0$$ Just to show that this makes sense, the simplest (classical) example is the 1-D harmonic oscillator. The potential is $E\left(R\right)=\frac{1}{2}kR^{2}$ , and according to equation (NUMBER): $$\frac{1}{m}\frac{\partial^{2}}{\partial R^{2}}\left(\frac{1}{2}kR^{2}\right)-\omega^{2}=0\Rightarrow\omega=\sqrt{\frac{k}{m}}$$ A well known result.

In the more general case, the matrix of IFC's is written $$C_{st}^{\alpha\beta}\left(\vec{R}\right)=\frac{\partial^{2}E}{\partial u_{s}^{\alpha}\partial u_{t}^{\beta}}$$ where $u_{s}^{\alpha}$ is the $\alpha$ -th Cartesian component of the deviation of the s-th atom from its equilibrium position.
The Fourier transform of the real-space matrix of IFC's, denoted $\tilde{C}_{st}^{\alpha\beta}\left(q\right)$ is called the Dynamical matrix: $$\tilde{C}_{st}^{\alpha\beta}\left(\vec{q}\right)=\sum_{\vec{R}}e^{-i\vec{q}\cdot\vec{R}}C_{st}^{\alpha\beta}\left(\vec{R}\right)$$ The physical meaning of the dynamical matrix is better understood as the second derivative of the energy with respect to deviations resulting from a lattice distortion in the form a wave with wave vector $\vec{q}$ . That is, when the deviation vector of the s-th atom is $\vec{u}_{s}=\vec{u}_{s}\left(\vec{q}\right)\cdot e^{i\vec{q}\cdot\vec{R}}$ : The phonon frequencies $\omega\left(\vec{q}\right)$ are found by diagonalizing the dynamical matrix, or solving the secular equation $$det\left|\frac{1}{\sqrt{M_{s}M_{t}}}\tilde{C}_{st}^{\alpha\beta}\left(\vec{q}\right)-\omega^{2}\left(\vec{q}\right)\right|=0$$

In this project, the computational scheme used to calculate the phonon dispersion involves:

  1. Calculation (and diagonalization) of the dynamical matrix on a grid of points reciprocal lattice.
  2. Calculation of the real-space IFC matrix by FFT $(N_{q}$ is the total number of q points): $$C\left(\vec{R}_{l}\right)=\frac{1}{N_{q}}\sum C\left(\vec{q}_{n}\right)e^{i\vec{q}\cdot\vec{R}_{l}}$$ .
  3. Calculation of the dynamical matrices in required q-points (denoted $\vec{q}'$) , along a selected path in the reciprocal lattice, by Fourier interpolating of the real space IFC's : $$C\left(\vec{q}'\right)=\sum_{l}C\left(\vec{R}_{l}e^{-\vec{q}'\cdot\vec{R}_{l}}\right)$$

Go to Top

Setup

Step 1 - Structure Relaxation

Unless otherwise specified, DFT calculations try to emulate materials at a temperature of 0 K. Naturally, structural parameters found in the scientific literature were recorded at some specific temperature, most commonly room temperature (300 K). On top of this, DFT also inevitably includes artificial (yet necessary) corrections to the model introduced by using pseudopotentials and the infamous exchange-correlation functional.

All this means that as a starting point for phonon calculations, using structural parameters (namely the lattice constant) straight from some reference might lead to inaccuracies. As a partial remedy, it is good practice to start all calculations with a structural relaxation.

What does structural relaxation mean? Well we first need to remember that DFT calculations are done using the Born-Oppenheimer approximation, which means that the nucleuses of atoms (we'll call them just 'atoms' for here on) don't move - only the electrons do. Structural Relaxation is in iterative calculation, which allows the atoms move (a little) in order to find a structural configuration that is more stable - i.e. has lower energy. The procedure can be simply described in the following figure:

Client Logo

In Quantum Espresso - a structural relaxation is done be specifying the 'calculation' flag. It can be either 'relax' or 'vc-relax' (vc meaning variable cell), the difference between the two is that in 'vc-relax' the cell parameter itself is allowed to change, whereas in 'relax' the cell is fixed but the atoms can move within the cell. Since it is known that PbTe maintains its Rock-Salt structure at low temperatures, we are interested in how the lattice parameter changes, and therefore we'll be using 'vc-relax'. Let's have a look at the input script for the relaxation run (notice the comments in green, after the ! mark):

 &control
    calculation='vc-relax', ! This means we're doing a variable cell relaxation
    restart_mode='from_scratch', ! This is where all the data (LARGE files!) is stored - devote some folder for this purpose
    prefix='PbTe', ! Name of the folder to be opened in the outdir. It's useful to use the same prefix for sequential calculations, as the data produced by this calculation be be used by others
    outdir = '/u/gery/outs//', ! This is where all the data (LARGE files!) is stored - devote some folder for this purpose
    pseudo_dir = '/u/gery/pseudo//', ! Folder where the pseudopotential files are located.
 /
 &system    
    nat = 2 ! Number of atoms in the unit cell - PbTe has two atoms in it's Rock-Salt structure
    ntyp = 2 ! Number of types of atoms
    ibrav = 2 ! type of crystal lattice - ibrav=2 means that we're using a FCC (Face Centered Cubic) cell.
    celldm(1) =12.7, ! Specifies the cell dimension - only one value is needed for a cubic lattice. I'm using an initial value of 12.7 bohr (1 bohr = 0.529177249 Angstrom)
    ecutwfc =60 ! This is the kinetic energy cutoff for wave functions (in Ry units, 1 Ry = 13.6 eV) - 
    ! the higher this is the more accurate the total energy value will be, but the run will take more time
 /
 &electrons
 /
 &ions
 /
&cell
 /
ATOMIC_SPECIES
Pb  207.2  Pb.pbe-d-hgh.UPF ! Specify the mass of atom type 1 (Pb), and the potential we're using - the potential must be place in the pseudo directory
Te  127.6  Te.pbe-hgh.UPF
ATOMIC_POSITIONS crystal
Pb 0.00 0.00 0.00
Te 0.5 0.5 0.5
K_POINTS AUTOMATIC
10 10 10 0 0 0 ! 

So, our initial guess is 12.7 bohr (6.72 Angstrom). This is an intentionally high guess (experimental values are around 6.4 Angstrom), so the algorithm will have something to do. Run the relaxation from the command line (don't forget to edit the PBS input to match your authorization and mail and also change WORK_DIR, PSEUDO_DIR and OUT_DIR to the ones you chose):

qsub mpi_PbTe-vcrelax.sh

We can fish out the new lattice parameter from the out log, with an awk command:

[gery@tamnun relaxation]$ awk '/CELL_/{x=NR+3}(NR<=x){print}' PbTe.vcrelax.out
CELL_PARAMETERS (alat= 12.70000000)
  -0.497353022  -0.000000000   0.497353022
   0.000000000   0.497353022   0.497353022
  -0.497353022   0.497353022   0.000000000
CELL_PARAMETERS (alat= 12.70000000)
  -0.493382555   0.000000000   0.493382555
  -0.000000000   0.493382555   0.493382555
  -0.493382555   0.493382555   0.000000000
CELL_PARAMETERS (alat= 12.70000000)
  -0.487426855  -0.000000000   0.487426855
  -0.000000000   0.487426855   0.487426855
  -0.487426855   0.487426855   0.000000000
CELL_PARAMETERS (alat= 12.70000000)
  -0.487117744  -0.000000000   0.487117744
  -0.000000000   0.487117744   0.487117744
  -0.487117744   0.487117744   0.000000000
CELL_PARAMETERS (alat= 12.70000000)
  -0.487117744  -0.000000000   0.487117744
  -0.000000000   0.487117744   0.487117744
  -0.487117744   0.487117744   0.000000000

We see here the new lattice vectors, in terms on the old one (not so straight forward but that's what it is). an original lattice vector belongs to the family (0.5 0.5 0), so we can compute our new lattice vector by $\frac{0.487117744}{0.5}\cdot12.7\approx12.372\ \left[bohr\right]\approx6.547\left[Angstrom\right]$. This is actually still a bit higher than most experimental values, which would require further investigation if we would want to publish this result in some journal. However, for our purposes it's enough that it minimizes the energy of the system, and from now one we'll use this value.

Step 2 - Phonon Calculation


This is the main step of the calculation. Let's have a look at the shell script (which both generates the input data for QE and executes the mpirun. it's called run_ph_PbTe.sh):

#!/bin/sh
#PBS -q  simphony_q
#PBS -l select=8:ncpus=12:mpiprocs=12 -l place=scatter
#PBS -M sasafhad@t2.technion.ac.il

PBS_O_WORKDIR=/u/gery/Phonon/week1/day3/PbTe/Final/ # you MUST update this to the right place!!!
cd $PBS_O_WORKDIR


PSEUDO_DIR="/u/gery/pseudo/"
OUT_DIR="/u/gery/outs/"
PW_COMMAND="/usr/local/espresso/bin/pw.x"
PH_COMMAND="/usr/local/espresso/bin/ph.x"

NAME="PbTe"
alat=12.372 # This is the lattice parameter that the relaxation gave (we started out with an initial guess of 12.7)
cutoff=60 # This is a relatively high cutoff - the simulation will be long but without such a cutoff the phonon calculation will be erroneous
 
# self-consistent calculation
cat > $NAME.scf.in << EOF
 &control
    calculation='scf' ! we're running a self consistent calculation (not vcrelax as we previously did)
    restart_mode='from_scratch',
    prefix='$NAME',
    outdir = '$OUT_DIR/'
    pseudo_dir = '$PSEUDO_DIR/',
 /
 &system    
    nat = 2
    ntyp = 2
    ibrav = 2
    celldm(1) =$alat, 
    ecutwfc =$cutoff
    lspinorb= .true. ! include the spin-orbit term in the calculation - in PbTe this term is not negligible
    noncolin=.true.
 /
 &electrons
    conv_thr =  1.0d-10
 /
ATOMIC_SPECIES
Pb  207.2  Pb.pbe-d-hgh.UPF
Te  127.6  Te.pbe-hgh.UPF
ATOMIC_POSITIONS crystal
Pb 0.00 0.00 0.00
Te 0.5 0.5 0.5
K_POINTS AUTOMATIC
10 10 10 0 0 0
EOF

mpirun -np 96 $PW_COMMAND < $NAME.scf.in > $NAME.scf.out

# phonons calculation on a grid - this will take a *LONG* time
cat > $NAME.phGrid.in << EOF
phonons on a grid
 &inputph
   prefix='$NAME',
   outdir = '$OUT_DIR/'
   ldisp=.true., ! the flag means we are doing the phonon calculation on a grid of points in the reciprocal lattice.
   fildyn='dyn',
   nq1=4, ! specifies the density of points in one direction in the reciprocal lattice
   nq2=4,
   nq3=4,
   tr2_ph=1.0d-14,
 /
EOF
mpirun -np 96 $PH_COMMAND < $NAME.phGrid.in > $NAME.phGrid.out

So we have two stages - a self consistent calculation (type scf) and a phonon calculation (which uses the PHonon executable, ph.x).

The scf calculation looks almost the same as the relaxation, the only difference being that now the lattice parameter is fixed (at our relaxed value), and that we're including a spin-orbit term in the hamiltonian. As suggested by reference 1 , in PbTe the spin-orbit term is not negligible when it comes to phonon calculations. Also, note that we're using relatively high cost parameters, namely a cutoff of 60Ry and a dense 10*10*10 k-point grid. This will take time, but is necessary because the phonon calculation will produce meaningless results (negative frequencies) otherwise.

The phonon calculation computes the dynamical matrix on each point on a 4*4*4 grid in q-space (a total of eight points). This is the raw data for all the subsequent fourier interpolation. Each point will be stored in xml format in files call dyn#, # ranging from 1 to 8. Let's have a look at a piece of the output log concerning the calculation on a single point in the grid, to see what is actually going on:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
     Computing dynamical matrix for 
                    q = (  -0.2500000   0.2500000  -0.2500000 )

      6 Sym.Ops. (no q -> -q+G )

     G cutoff =  930.5320  (    310 G-vectors)     FFT grid: ( 45, 45, 45)
     number of k points=   440
.
.
.
.
<--Lets zip to the interesting part-->
.
.
.
.


     End of self-consistent calculation

     Convergence has been achieved 

     Number of q in the star =    8
     List of q in the star:
          1  -0.250000000   0.250000000  -0.250000000
          2   0.250000000  -0.250000000  -0.250000000
          3   0.250000000  -0.250000000   0.250000000
          4   0.250000000   0.250000000   0.250000000
          5  -0.250000000  -0.250000000  -0.250000000
          6  -0.250000000  -0.250000000   0.250000000
          7  -0.250000000   0.250000000   0.250000000
          8   0.250000000   0.250000000  -0.250000000

     Diagonalizing the dynamical matrix

     q = (   -0.250000000   0.250000000  -0.250000000 ) 

 **************************************************************************
     omega( 1) =       0.983445 [THz] =      32.804184 [cm-1]
     omega( 2) =       0.983445 [THz] =      32.804184 [cm-1]
     omega( 3) =       1.587818 [THz] =      52.963920 [cm-1]
     omega( 4) =       2.477828 [THz] =      82.651447 [cm-1]
     omega( 5) =       2.477828 [THz] =      82.651447 [cm-1]
     omega( 6) =       3.391280 [THz] =     113.120929 [cm-1]
 **************************************************************************

     Mode symmetry, C_3v (3m)   point group:

     omega(  1 -  2) =         32.8  [cm-1]   --> E    L_3           
     omega(  3 -  3) =         53.0  [cm-1]   --> A_1  L_1           
     omega(  4 -  5) =         82.7  [cm-1]   --> E    L_3           
     omega(  6 -  6) =        113.1  [cm-1]   --> A_1  L_1    

I've added line numbers for reference. Some things to notice:


line 1-2 : This is the calculation for point (-0.25 0.25 -0.25)

lines 23-32 : QE recognizes the symmetry of this bravais lattice, and lists 8 equivalent points in the first Brillouin zone

lines 39-44 : the vibrational frequencies (eigenvalues on the dynamical matrix) - this is what we want!

lines 47-52 : QE recognizes the symmetry of the normal modes (we won't be going into this subject, but QE can do it.)


Run the calculation with the command:
qsub run_ph.PbTe.sh

This will take some time (about 8 hours, using 8 nodes). If you don't want to wait, you may use the reference data available in the folder reference/phonon_generation to proceed to the next step.


Step 3 - Fourier Interpolation and Plotting


As explained in the Theory section, once the dynamical matrices are calculated over a grid in q-space, we want to convert them to real space Interatomic Force Constants, and then convert back to dynamical matrices over a path of selected q-points..

This part uses the shell script called run_and_plot_dispersion_PbTe.sh, which does all of that and plots the results as well.


Run the script with:
qsub run_and_plot_dispersion_PbTe.sh

Let's have a look at the script (its name is run_and_plot_dispersion_PbTe.sh):
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#!/bin/sh
#PBS -q  simphony_q
#PBS -l select=1:ncpus=12:mpiprocs=12 -l place=scatter
#PBS -M sasafhad@t2.technion.ac.il

PBS_O_WORKDIR=/u/gery/Phonon/week1/day3/PbTe/Final/ # you MUST update this to the right place!!!
cd $PBS_O_WORKDIR


PSEUDO_DIR="/u/gery/pseudo/"


MATDYN_COMMAND="/usr/local/espresso/bin/matdyn.x"
PLOTBAND_COMMAND="/usr/local/espresso/bin/plotband.x"
Q2R_COMMAND="/usr/local/espresso/bin/q2r.x"

NAME="PbTe"

# run q2r to calculate the IFC's (in realspace) from the dynamical matrices (in q space)

cat > q2r.in <<EOF
 &input
  fildyn='dyn.xml',
  zasr='crystal',
  flfrc='$NAME.444.fc'
 /


EOF

mpirun -np 12 $Q2R_COMMAND < q2r.in > q2r.out


#  recalculate omega(q) from C(R) using matdyin
 
cat > matdyn_disp.in <<EOF
 &input
    asr='crystal',  	amass(1)=207.2,
	amass(2)=127.6, q_in_cryst_coord=.TRUE.,
    flfrc='$NAME.444.fc.xml', flfrq='$NAME.freq'
 /
         398
        0.0000000000    0.0000000000    0.0000000000     
        0.0000000000    0.0054347826    0.0054347826     
    <-----lots of points....--->
        0.0063291139    0.5000000000    0.5000000000     
        0.0000000000    0.5000000000    0.5000000000     
EOF


mpirun -np 12 $MATDYN_COMMAND < matdyn_disp.in > matdyn_disp.out

rm freq*

#  make freeuency files for plotting the phonon dispersions
cat > plotband.in <<EOF
$NAME.freq
0 200
freq.plot
freq.ps
0.0
0.0 
0.1
EOF

mpirun -np 12 $PLOTBAND_COMMAND < plotband.in > /dev/null

# unify freq.plot.* files to one file (to be used by gnuplot)

cat freq.plot* > freq.plot

#make gnuplot for the dispersion

cat > plot_dispersion.gnu <<EOF
set terminal post eps enhanced solid color "Helvetica" 20
set output "gnuplot_dispersion.eps"
set key off

set xrange [0:4.3065]
set yrange [0:150]
set arrow from 1,0. to 1,150 nohead  lw 3
set arrow from 2.585,0. to 2.585,150 nohead  lw 3
set arrow from 1.5,0. to 1.5,150 nohead  lw 3
set arrow from 3.4142,0. to 3.4142,150 nohead  lw 3
set ylabel "{/Symbol w} (cm^{-1})"
unset xtics
set label "{/Symbol G}" at -0.05,-4
set label "X" at 0.95,-4
set label "W" at 1.45,-4
set label "{/Symbol G}" at 2.585,-4
set label "L" at 3.37,-4
set label "X" at 4.3065,-4

plot "freq.plot" u 1:2 w p pt 3 

EOF

# run matdyn for DOS calculation

cat > matdyn_ex.in <<EOF
 &input
  asr='crystal',
	amass(1)=207.2,
	amass(2)=127.6, 
	flfrc='$NAME.444.fc.xml',
	flfrq='$NAME.dos.freq'
	dos=.true.,
	fldos='$NAME.dos'
	deltaE=1.d0,
	nk1=20, nk2=20, nk3=20,
 / 

EOF

mpirun -np 12 $MATDYN_COMMAND < matdyn_ex.in > matdyn_ex.out


# make gnuplot for DOS

cat > plot_dos.gnu <<EOF
set terminal post eps enhanced solid color "Helvetica" 20
set output "gnuplot_dos.eps"
set key off

set xrange [0:150]
set ylabel "VDOS"
set xlabel "{/Symbol w} (cm^{-1})"

plot "$NAME.dos" u 1:2 w l lw 2


EOF


#plot everything and convert it to png

gnuplot plot_dispersion.gnu
gnuplot plot_dos.gnu
sh eps2png gnuplot_dispersion.eps gnuplot_dispersion.png
sh eps2png gnuplot_dos.eps gnuplot_dos.png


What the script does:

Lines 19-31 : Calculate the IFC's in real space using a Quantum Espresso module named q2r.x

Lines 34-58 : Calculate the dynamical matrices over select q-points (lots of them), using a QE module name matdyn.x. I created the path of Q points with a visualization program named XCRYSDEN, as explain in the Manual section

Lines 55-66 : Create a data file of the phonon frequencies for plotting

Lines 72-96 : Create a gnuplot file that has plotting instructions for plotting the phonon dispersion relation $\omega\left(\vec{q}\right)$

Lines 98-115 : Use the real-space IFC's and calculates phonons over a 4*4*4 grid for a plot of the vibrational density of states (DOS)

Lines 118-132 : Create a gnuplot file with instructions for plotting the phonon density of states

Lines 135-140 : Plot both graphs with gnuplot, and convert the resulting .eps files to .png with a script called eps.png

The script should run pretty fast (a few seconds). Let's review the results.
Go to Top

Results

The resulting phonon dispersion is:

Client Logo
Calulated phonon dispersion - image made with gnuplot.


The Q-point path the reciprocal lattice over which the dispersion is mapped

Client Logo
Q-point path used - image from XCRYSDEN

While the dispersion may appear like just a blurb of lines, it actually has some physically correct features. For one, around the $\Gamma$ point the acoustic modes converge to 0. Another feature is the appearance of optical modes (modes with non-zero frequency at $\Gamma$), as we would expect from a lattice with more that one atom type. The total number of bands is also correct - while sometimes the frequencies are degenerate (i.e. appear more than once) it's possible to see that there are a total of 6 bands, which is good because we're expecting 3N=6 bands.

The results are very similar to those calculated by reference 2. Compare the two:

this is image #1 this is image #2


There are a few differences, particularly between point W and $\Gamma$, which is reasonable since the lattice parameter used for reference 2.


Also, compare to experimental result from reference 6, mainly for a qualitative comparison, as the Q-point path isn't that same and the results were obtained at room temperature:
this is image #1 this is image #2

As for the density of states, my result (on the left) is brought together with another computational result from reference 5:

this is image #1 this is image #2
Go to Top

Manual

Sending tasks to the PBS queue:

A QE task is run by sending a shell script to the queueing system. Here is an example script (you can either copy the text or download it here):

#!/bin/sh

#------Script Parameters-----

#PBS -q  simphony_q 

# 		This is the queue you're asking to use (I'm asking for 'simphony_q'). 
# 		You must have authorization to the queue. 
#		Everyone has authorization at least to the general 'all_l_p' queue.

#PBS -l select=2:ncpus=12:mpiprocs=12  

# 		Here you're specifying the computer resources. 
# 		In the example use M chunks (nodes), use N (=< 12) CPUs for N mpi tasks on each of M nodes.

#PBS -M sasafhad@t2.technion.ac.il 

# 		Specify you're emal address (must be a technion address) - you'll
#		receive notification on the execution, completion and errors of each run to your mail.


#------Script Body------

PBS_O_WORKDIR=/u/gery/PSU/PbTe/ 
# 		Specifies the working directory, were the 
#		input files are and where output files will be sent to.
cd $PBS_O_WORKDIR

  
#------Actual Running Command (run MPI executable with M*N processes)------
mpirun -np 24   /usr/local/espresso/bin/pw.x < scf.in > scf_new.out  
#		The "np" must be equal the number of chunks multiplied by the number of "ncpus"

The task is then sent to the queue with the command:

qsub scriptname.sh

You can follow the execution status with the command:

qstat -u your_user_name

You can also check who's in line in your queue with the following command:

qstat | grep queue_name



XCRYSDEN - Visualization for Quantum Espresso

XCYRSDEN is a crystalline visualization program, that runs on linux. It's useful for Quantum Espresso because it can read the input and output files of the main suite, pw.x (PWscf).

To install XCrysden, download the program and unzip the binary package. I used XCRYSDEN for visualizing the unit cell (the picture in the Overview), and for creating a Q-point path. Using the program is self explanatory.

If you plan to use it a lot, it's useful to set an alias command (a command line shortcut). Do this by adding the following line to your .bashrc file (located in your user's root directory):

alias xcrysden='/./PATH_TO_FILE/xcrysden'

Some tips on debugging input scripts:

If a run fails for some reason, the reason for failure will usually be written at the end of the log file. To demonstrate this, in the following example I intentionally made a mistake in the address of the pseudopotential files. From reading the .out file we can see QE complain:

[gery@tamnun PbTe]$ cat scf_new.out

    ...

     Current dimensions of program PWSCF are:
     Max number of different atomic species (ntypx) = 10
     Max number of k-points (npk) =  40000
     Max angular momentum in pseudopotentials (lmaxx) =  3
     Waiting for input...
     Reading input from standard input

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                   %%%%%%%%%%%%%%%%
     Error in routine readpp (48):
     file /u/gery/WRONG_FOLDER_EXAMPLE/Pb.rel-pz-dn-kjpaw_psl.0                                                   .2.2.UPF not found
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                   %%%%%%%%%%%%%%%%

     stopping ...
APPLICATION TERMINATED WITH THE EXIT STRING: Hangup (signal 1)

The Quantum Espresso homepage has useful documentation. Particularly useful are the input data descriptions, where each flag is explained. The best resource on specific questions is the QE Forum, - on the page there is a search bar (make sure you check the 'forum' option). Most questions have already been asked and answered

In General, the documentation assumes a very good command of solid state physics and the subtle jargon of DFT, which can be annoying for new users.

Go to Top