Untitled
Untitled
Untitled
Email: [email protected]
Web: www.jennystanford.com
All rights reserved. This book, or parts thereof, may not be reproduced
in any form or by any means, electronic or mechanical, including
photocopying, recording or any information storage and retrieval
system now known or to be invented, without written permission from
the publisher.
Preface ix
1 Introduction 1
1.1 How to read and use the book? 1
1.2 What do we need to run a program? 3
1.3 What we get, and what we do not get? 3
1.4 Organization of the book 4
2 Software Installation 5
2.1 Preparing the operating systems 5
2.1.1 Ubuntu Linux 6
2.1.2 Windows 8
2.1.3 macOS Catalina 14
2.2 Installation of Quantum ESPRESSO and its supporting
software 14
2.3 VirtualBox approach 18
2.4 Processing input and output files 20
2.4.1 Basic execution of Quantum ESPRESSO
commands 20
2.4.2 Choice of plotting software 21
2.4.3 Obtaining example files for hands‐on tutorials 23
Bibliography 343
Index 357
Preface
Introduction
1 It is noted that we have many principles in the first‐principles. Thus, we do not say
You will be surprised if you know what we can get from Quantum
ESPRESSO. The following concepts are what we get from Quantum
ESPRESSO: (1) lattice structure of solid and X‐ray diffraction spectra,
(2) electronic band structure and density of states, (3) phonon
dispersion, (4) optical absorption spectra, (5) Raman spectra,
(6) electronic transport properties such as carrier mobility, and
(7) thermal properties such as thermal expansion and thermal
conductivity. These concepts are subjects of solid‐state physics that
you do not need to know now. For the readers who are not familiar
with physics, we will give minimum information about solid‐state
physics by using some equations to understand the special words.
Although solid‐state physics contents are minimal, you can use this
book as the first step of understanding solid‐state physics. It is highly
desired for you to read some solid‐state physics textbooks, too.
It is important to note that running Quantum ESPRESSO does
NOT mean that (1) we can explain the phenomena of solid‐state
4 | Introduction
Software Installation
This chapter will guide the readers to install all necessary software
for running Quantum ESPRESSO. We expect that the readers already
have a computer with an operating system (OS) such as Linux,
Windows, and macOS. Firstly, we will show you the minimum
requirements for each OS to avoid start‐up problems in installing
Quantum ESPRESSO. Next, we will show the main method to install
Quantum ESPRESSO and its supporting software. Finally, we will
show you how to run the simplest running instance of Quantum
ESPRESSO.
In Ubuntu, please check that you can open a Terminal window from
the Activities panel. Terminal is a place in the OS where we
can interact with the “shell”. Simply saying, the shell is a program
that takes commands from the keyboard and returns them to the
operating system to perform.
In the old days, the shell was the only user interface available
on a Unix‐like system. Nowadays, we have graphical user interfaces
(GUIs) in addition to the shell, which is a kind of command‐line
interfaces (CLIs). To interact with the shell in the terminal, we can
type a line of “command” or (lines of commands) corresponding to
some computer operations.
The basic structure of Ubuntu Terminal is shown in Fig. 2.1.
We can see there is a prompt showing the login username, machine
hostname, and a current working directory (“∼” sign, indicating the
so‐called “home” directory of the user). The prompt ends with a
dollar ($) sign. The cursor next to the $ sign is the starting point
Quantum ESPRESSO Course for Solid‐State Physics |7
where we can type the actual Linux commands. There are so many
commands for various purposes, but we will only use a few of them to
install software and manage our Quantum ESPRESSO jobs and files.
To update all software in Ubuntu, on the terminal, type the
following commands (next to the $ sign) and press enter after each
line:
You will be asked for “your password” (not root password) when
executing the above commands because the administrator (often
called “root”) privileges are required when installing or updating the
software in Ubuntu. The root privileges are represented by the sudo
command.
The apt command corresponds to the “Advanced Packaging
Tool”, a command‐line tool that helps in handling packages in Ubuntu.
Its main task is to retrieve the information and packages from the
authenticated sources (mostly from the internet) for installation,
upgrade, and removal of packages along with their dependencies.
In the above example, the update and upgrade texts are a set of
commands to ensure your Ubuntu software is all updated.
We also need some development tools and libraries, such as
Git (version control system), GNU Wget (file retrieval software),
GCC/G++(C/C++ compiler), GFortran (Fortran compiler), LAPACK
(linear algebra package), FFTW (Fourier transform library), and
Open MPI (parallel computing implementation). All of them are
8 | Software Installation
If you are prompted by a “Yes/No” question, just type “y” (or “Yes”)
and press enter to agree with the installation.
2.1.2 Windows
Figure 2.2 Windows Update. You can click Check for updates and
Download and install now to guarantee your Windows is in the latest
version.
Figure 2.3 Windows specifications in About your PC. Make sure that you
are running at least Windows 10 version 1903. Again, for a native support
of Linux graphical applications in WSL, we recommend at least Windows 11
build 22000.
from BIOS, which can be accessed before booting your OS. The key
you should press on the keyboard to access BIOS will depend on
the manufacturer of the computer. You can check from the computer
screen when you turn it on. Assuming the key is F2 (or DEL), below
are the steps to enable the virtualization from the BIOS settings:
• Press the F10 key and select Yes and press the Enter key to save
changes and Reboot into Windows.
wsl --install
Figure 2.6 First Ubuntu instances through WSL and creation of a user
account.
$ explorer.exe .
You can work with files normally from here (see Fig. 2.7). Use
drag and drop, copy and paste them, or even open them directly in
Windows applications.
Quantum ESPRESSO Course for Solid‐State Physics | 13
Besides typing the above command, you can also open the $HOME
directory from Windows Explorer by directing the explorer to the
following address:
\\wsl$\Ubuntu\home\quantum
is the same as that in the “real” Linux OS, which will be explained in
Sec. 2.2.
The command is quite long, and it may look broken into several lines
in this book or your terminal, but it is still a single line of command.
You should ensure typing the command properly to the terminal as a
unified entity, without pushing “Enter” key.
By having Homebrew in macOS, the prerequisites for the
Quantum ESPRESSO installation can be obtained by another single
line of command:
Again, no matter how the above (long) command looks broken into
separate lines, it is a single line of command.
With the various OSs that the readers have, we are ready to install
Quantum ESPRESSO using the same approach on the terminal.
Quantum ESPRESSO Course for Solid‐State Physics | 15
We will install Quantum ESPRESSO and an additional Wannier90
package, either in “real” Ubuntu Linux, WSL Ubuntu, or macOS.
For beginner users of Ubuntu Linux and WSL Ubuntu, the
following commands are sufficient to install Quantum ESPRESSO and
Wannier90:
$ mkdir opt
$ cd opt
In the opt folder, we can download the latest source files of Quantum
ESPRESSO using the Git utility as follows:
After the source files are downloaded, we can enter the q-e folder:
$ cd ~/opt/q-e/
$ ./ configure
16 | Software Installation
$ make all
$ make w90
The purpose of the first line above is to make executable files for
all standard Quantum ESPRESSO capabilities, such as electronic
structure, phonon dispersion, and optical properties calculations.
The second line is to make the additional Wannier90 package
available in the binary folder. In addition to this compilation process,
our OS must be “taught” how to find the executable files and include
them in its variable path. For this purpose, we can execute two more
commands to update the path:
$ pw.x
We will learn how to use this command briefly in Sec. 2.4 and
more extensively in Chapter 3. To exit the Quantum ESPRESSO
environment generated from the above command line, push CTRL+c
keystroke.
We can install some additional software to support our work‐
flow when performing first‐principles calculations with Quantum
ESPRESSO. Executing the following commands line by line on the
terminal will give us all necessary software:
Quantum ESPRESSO Course for Solid‐State Physics | 17
Some notes:
$ bash QEinstall.sh
18 | Software Installation
1. Quantum Mobile:
https://quantum-mobile.readthedocs.io
2. MateriApps LIVE!:
https://cmsi.github.io/MateriAppsLive/
3. QE‐2021: http://qe2021.ijs.si
Each virtual machine has its username and password, which may
change over time. Therefore, it is better to check the documentation
on one of the websites above, depending on which virtual machine
you like. However, all of them use the same capability of VirtualBox to
host a virtual machine. Here we just explain how to install Virtual Box
and add the virtual machine. The steps to install and use the virtual
machine is as follows:
Figure 2.8 Adding a virtual machine in VirtualBox. Click the Add button,
either from the Machine menu or from the toolbar.
Command
Text file Text file
Command Purpose
pw.x SCF and NSCF calculations
.................................................................................................................................................................
$ jupyter -lab
Note that some traditional Python users may prefer to directly access
the Jupyter notebook from the command line by typing:
$ jupyter -notebook
Chapter 3. To download the files from the terminal, we can use the
one‐line git clone utility:
1 $ cd ~/QE -SSP/gr/scf/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
1 If you do not find ‘mpirun’ command, you need to install libopenmpi‐dev (see
Chapter 2).
Quantum ESPRESSO Course for Solid‐State Physics | 27
‐ Line 1: Go to the scf directory that includes the input files.
‐ Line 2: Run pw.x (pw = plane‐wave, x = executable file) by using in
parallel with 4 processors (-np 4) by the command mpirun (running
a program with parallel processors) with the input file is scf.in and
the output file is scf.out. The symbol & makes the command run in
the background. Here, we select 4 processes for parallel calculations,
but the readers can run with serial calculations (without mpirun, that
is pw.x < scf.in > scf.out &) or any number of processes (e.g.,
-np 8 for Intel Core i7 or Core i9), depending on your computer. For
detailed commands for running in parallel, the readers can find on
Sec. 6.2.3.
How to check: To check whether or not the output file exists, the
readers can use the ls command by typing:
$ ls
scf.in scf.out
If the readers see the scf.out, the readers can check the contents of
scf.out by tail command as follows:
$ tail scf.out
The tail command prints the last few lines of the scf.out file. If the
calculation normally finishes, a message JOB DONE is written at the
end of this file as
=------------------------------------------------------=
JOB DONE.
=------------------------------------------------------=
Input file: Now, let us explain the details of the input file. To see
the input file, the readers can use vi editor by typing:
$ vi scf.in
28 | Hands‐On Tutorials of Quantum ESPRESSO
To quit the vi editor without saving any changes, the readers must
press : (colon). Then the cursor should reappear in the lower‐left
corner of the screen beside a colon prompt. After that, the readers
need to enter q! to quit the file without saving. Since vi editor
is usually available in all Linux distributions, we would like to use
vi. However, the readers can use any text editor, such as Emacs or
Notepad, to open scf.in directly. When the readers open scf.in file,
the readers will see the input variables as
QE‑SSP/gr/scf/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 2.4623
10 c = 10.0
11 nat = 2
12 ntyp = 1
13 occupations = 'smearing '
14 smearing = 'mv '
15 degauss = 0.02
16 ecutwfc = 60
17 /
18 & ELECTRONS
19 mixing_beta = 0.7
20 conv_thr = 1.0D-6
21 /
22 ATOMIC_SPECIES
23 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
24 ATOMIC_POSITIONS (crystal)
25 C 0.333333333 0.666666666 0.500000000
26 C 0.666666666 0.333333333 0.500000000
27 K_POINTS (automatic)
28 12 12 1 0 0 0
Continued
30 | Hands‐On Tutorials of Quantum ESPRESSO
$ xcrysden &
Atomic positions
c = 10.0 Å
a a = 2.4623 Å
4 2
Trigonal P v3 = a(0, 0, ac )
....................................................................................................................................................................
v1 = ( 2a , 2b , 2c ), v2 = (− 2a , 2b , 2c ),
11 Orthorhombic
v3 = (− 2a , − 2b , 2c )
....................................................................................................................................................................
Continued
34 | Hands‐On Tutorials of Quantum ESPRESSO
Monoclinic
v1 = ( 2a , 0, − 2c ), v2 = (b cos γ, b sin γ, 0),
13 base‐centered
v3 = ( 2a , 0, 2c )
(unique axis c)
....................................................................................................................................................................
Monoclinic
base‐centered v1 = ( 2a , 2b , 0), v2 = (− 2a , 2b , 0, 0),
(unique axis v 3 = (c cos β, 0, c sin β)
−13
b)
....................................................................................................................................................................
ATOMIC_POSITIONS (alat)
C 0.0000000 0.5773503 2.0306218
C 0.5000000 0.2886751 2.0306218
QE‑SSP/gr/scf/scf.out
If the readers want to see only the total energy from scf.out, the
readers can use the grep command to find the lines containing the
symbol ! from scf.out as follows:
$ grep ! scf.out
Try It Yourself
1 $ cd ~/QE -SSP/gr/ecut/
2 $ ./run.sh &
‐ Line 1: Change directory to the ecut that includes the input files.
‐ Line 2: Run a bash script file run.sh, which generates and runs
many jobs with changing the cut‐off energies.
How to check: To check whether or not the output file exists, the
readers can use the ls command by typing:
$ ls
The ls command displays a listing of the many input and output files
including run.sh in the ecut directory as
QE‑SSP/gr/ecut/run.sh
1 #!/bin/bash
2 # Convergence test of cut -off energy.
3 # Set a variable ecut from 20 to 80 Ry.
4 for ecut in 20 22 24 26 30 35 40 45 50 55 60 65 \
5 70 75 80; do
6 # Make input file for the SCF calculation.
7 # ecutwfc is assigned by variable ecut.
8 cat > ecut.$ecut.in << EOF
9 & CONTROL
10 calculation = 'scf '
11 pseudo_dir = '../ pseudo/'
12 outdir = '../tmp/'
13 prefix = 'gr '
14 /
15 &SYSTEM
16 ibrav = 4
17 a = 2.4623
18 c = 10.0
19 nat = 2
20 ntyp = 1
21 occupations = 'smearing '
22 smearing = 'mv '
23 degauss = 0.02
24 ecutwfc = ${ecut}
25 /
26 & ELECTRONS
27 mixing_beta = 0.7
28 conv_thr = 1.0D-6
29 /
30 ATOMIC_SPECIES
31 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
32 ATOMIC_POSITIONS (crystal)
33 C 0.333333333 0.666666666 0.500000000
34 C 0.666666666 0.333333333 0.500000000
35 K_POINTS (automatic)
36 12 12 1 0 0 0
37 EOF
38 # Run SCF calculation.
39 mpirun -np 4 pw.x <ecut.$ecut.in > ecut.$ecut.out
40 # Write cut -off and total energies in calc -ecut.dat.
41 awk '/!/ {printf "%d %s\n",'$ecut ',$5}' ecut.$ecut.
out >> calc -ecut.dat
42 # End of for loop
43 done
38 | Hands‐On Tutorials of Quantum ESPRESSO
$ vi calc -ecut.dat
QE‑SSP/gr/ecut/calc‑ecut.dat
1 20 -23.58222934
2 22 -23.74580589
3 24 -23.82861718
4 ...
5 70 -23.91018738
6 75 -23.91023390
7 80 -23.91024577
QE‑SSP/gr/ecut/plot‑ecut.ipynb
Try It Yourself
System Shifting
Cubic, tetragonal, orthorhombic, and monoclinic systems 1
....................................................................................................................................................................
(a) b2 b2
Γ b1 Γ b1
shift = 0 shift = 1
(b) b2 b2 b2
Γ b1 Γ b1 Γ b1
shift = 1 &
shift = 0 shift = 1
symmetrization
Figure 3.4 Schematic diagram of the k‐point sampling in a cubic cell (a) or a
hexagonal cell (b). We also illustrate how to shift the k‐points by specifying
“shift=1” and “symmetrization”. Red points denote the inequivalent k‐points
in the first Brillouin zone.
Next, we will learn how to perform the convergence test for this
k‐points grid.
How to run: To run this tutorial, the readers should type the
following command lines:
1 $ cd ~/QE -SSP/gr/k-point/
2 $ ./run.sh &
$ ls
QE‑SSP/gr/k‑point/run.sh
1 #!/bin/bash
2 # Convergence test of k-points grid.
3 # Set a variable k-point from 4 to 14.
4 for k in 4 5 6 7 8 9 10 11 12 13 14; do
5
6 # Make input file for the SCF calculation.
7 # k-points grid is assigned by variable n.
8 cat > kpoint.$k.in << EOF
9 & CONTROL
10 calculation = 'scf '
11 pseudo_dir = '../ pseudo/'
12 outdir = '../tmp/'
13 prefix = 'gr '
14 /
15 &SYSTEM
16 ibrav = 4
17 a = 2.4623
18 c = 10.0
19 nat = 2
20 ntyp = 1
21 occupations = 'smearing '
22 smearing = 'mv '
23 degauss = 0.02
24 ecutwfc = 40
25 /
26 & ELECTRONS
27 mixing_beta = 0.7
28 conv_thr = 1.0D-6
29 /
30 ATOMIC_SPECIES
31 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
32 ATOMIC_POSITIONS (crystal)
33 C 0.333333333 0.666666666 0.500000000
34 C 0.666666666 0.333333333 0.500000000
44 | Hands‐On Tutorials of Quantum ESPRESSO
35 K_POINTS (automatic)
36 ${k} ${k} 1 0 0 0
37 EOF
38
39 # Run pw.x for SCF calculation.
40 mpirun -np 4 pw.x <kpoint.$k.in >kpoint.$k.out
41 # Write the number of k-points (= k*k*1) and
42 # the total energy in calc -kpoint.dat
43 awk '/!/ {printf "%d %s\n",'$k*$k ',$5}' kpoint.$k.out
>> calc -kpoint.dat
44 # End of for loop.
45 done
$ vi calc -kpoint.dat
QE‑SSP/gr/ecut/calc‑kpoint.dat
1 16 -23.91186965
2 25 -23.90766485
3 36 -23.90599780
4 ...
5 144 -23.90913362
6 169 -23.90963897
7 196 -23.90953385
QE‑SSP/gr/k‑point/plot‑kpoint.ipynb
Try It Yourself
minf(x), x ∈ Rn , (3.1)
SCF calculation
NO
BFGS algorithm
Check convergence criteria
YES
Optimal
atomic positions
xk+1 = xk − H−1
k ∇f(xk ),
sk = xk+1 − xk ,
yk = ∇f(xk+1 ) − ∇f(xk ), (3.2)
Hk sk sTk Hk yk yTk
Hk+1 = Hk − , with k = 0, 1, 2, …
sTk Hk sk yTk sk
+
1 $ cd ~/QE -SSP/gr/relax/
2 $ mpirun -np 4 pw.x <relax.in > relax.out &
$ vi relax.in
QE‑SSP/gr/relax/relax.in
1 & CONTROL
2 calculation = 'relax '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 etot_conv_thr = 1.0D-5
7 forc_conv_thr = 1.0D-4
8 /
9 &SYSTEM
10 ibrav = 4
11 a = 2.4623
12 c = 10.0
13 nat = 2
14 ntyp = 1
15 occupations = 'smearing '
16 smearing = 'mv '
17 degauss = 0.02
18 ecutwfc = 60
19 /
Quantum ESPRESSO Course for Solid‐State Physics | 49
20 & ELECTRONS
21 mixing_beta = 0.7
22 conv_thr = 1.0D-9
23 /
24 &IONS
25 ion_dynamics = 'bfgs '
26 /
27 ATOMIC_SPECIES
28 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
29 ATOMIC_POSITIONS (crystal)
30 C 0.333333333 0.666666666 0.500000000 0 0 0
31 C 0.666666666 0.333333333 0.400000000
32 K_POINTS (automatic)
33 12 12 1 0 0 0
$ vi relax.out
To search forward for final coordinates, press Esc and then enter
/final, vi editor will search the word 'final' in relax.out.
QE‑SSP/gr/relax/relax.out
Forces acting on atoms (cartesian axes , Ry/au):
ATOMIC_POSITIONS (crystal)
C 0.3333333330 0.6666666660 0.5000000000 0 0 0
C 0.6666666660 0.3333333330 0.4999977356
Quantum ESPRESSO Course for Solid‐State Physics | 51
End final coordinates
< 1.0 × 10−5 Ry and the force < 1.0 × 10−4 Ry/Bohr. The
second carbon atoms is moved from the initial position
(0.6666666660 0.3333333330 > 0.4000000000) to the final
position (0.6666666660
>
0.3333333330 0.4999977356).
Note: The Total force in the output file is the square root
of the sum of all the squared force components rather than the
sum of the magnitudes of the individual forces on the atoms. If the
Total SCF correction is large or comparable to the Total force,
the output file will show a message "SCF correction compared to
forces is large: reduce conv_thr to get better values".
It usually means the readers need to try with a relatively smaller
conv_thr.
Visualizing optimized ionic positions: The readers can visualize
the optimized ionic positions by using XCrySDen as
$ xcrysden &
Then the readers can go through the following steps: File → Open
PWscf → Open PWscf Output File and select relax.out. In the box
[PWSCF Input: “relax.out”], we click on the OK button, then in the
box [Question], we select Display All Coordinates as Animation
and click on the Continue button. The initial or optimized ionic
positions can be showed by choosing Current slide in the box
[Animation Control Center], as shown in Fig. 3.7.
Try It Yourself
1 1
Quantum ESPRESSO Course for Solid‐State Physics | 53
$ vi vc -relax.in
QE‑SSP/gr/vc‑relax/vc‑relax.in
1 & CONTROL
2 calculation = 'vc -relax '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 etot_conv_thr = 1.0D-5
7 forc_conv_thr = 1.0D-4
8 /
9 &SYSTEM
10 ibrav = 4
11 a = 2.5
12 c = 15.0
13 nat = 2
14 ntyp = 1
15 occupations = 'smearing '
16 smearing = 'mv '
17 degauss = 0.02
18 ecutwfc = 60
19 /
20 & ELECTRONS
21 mixing_beta = 0.7
22 conv_thr = 1.0D-9
23 /
24 &IONS
25 ion_dynamics = 'bfgs '
26 /
27 &CELL
28 cell_dynamics = 'bfgs '
29 press_conv_thr= 0.05
30 cell_dofree = '2Dxy '
54 | Hands‐On Tutorials of Quantum ESPRESSO
31 /
32 ATOMIC_SPECIES
33 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
34 ATOMIC_POSITIONS (crystal)
35 C 0.333333333 0.666666666 0.500000000
36 C 0.666666666 0.333333333 0.500000000
37 K_POINTS (automatic)
38 12 12 1 0 0 0
$ vi vc -relax.out
To search forward for final coordinates, press Esc and then enter
/final, vi editor will display as
QE‑SSP/gr/vc‑relax/vc‑relax.out
Computing stress (Cartesian axis) and pressure
ATOMIC_POSITIONS (crystal)
C 0.3333333330 0.6666666660 0.5000000000
C 0.6666666660 0.3333333330 0.5000000000
End final coordinates
Try It Yourself
C.rel-pbe-n-rrkjus_psl.1.0.0.UPF
http://pseudopotentials.quantum-espresso.org/legacy_
tables/hartwigesen-goedecker-hutter-pp.
A recommended way of selecting a good PP is based on testing
and comparing the preliminary results to experimental data. Note
that when comparing Quantum ESPRESSO calculation results, you
should not directly compare the energy values by using different PPs.
Since each PP has a different option or condition, using a different
type of PP results in different energy values, which does not have
any meaning. It is recommended to compare the properties such as
lattice constant, energy band gap, etc., with experiments for observed
materials. Since we showed how to obtain the lattice constant of
graphene in Sec. 3.1.5, in this tutorial, we will compare the lattice
constant for several PPs.
How to run: To run this tutorial, the readers should type the
following command lines:
1 $ cd ~/QE -SSP/gr/pps/
2 $ ./run.sh &
$ ls
run.sh
vc -relax.C.pbe -hgh.UPF.in
vc -relax.C.pbe -hgh.UPF.out
vc -relax.C.pbe -n-kjpaw_psl .1.0.0. UPF.in
vc -relax.C.pbe -n-kjpaw_psl .1.0.0. UPF.out
vc -relax.C.pbe -n-rrkjus_psl .0.1. UPF.in
vc -relax.C.pbe -n-rrkjus_psl .0.1. UPF.out
vc -relax.C.pz -hgh.UPF.in
vc -relax.C.pz -hgh.UPF.out
vc -relax.C.pz -n-kjpaw_psl .0.1. UPF.in
vc -relax.C.pz -n-kjpaw_psl .0.1. UPF.out
vc -relax.C.pz -n-rrkjus_psl .0.1. UPF.in
vc -relax.C.pz -n-rrkjus_psl .0.1. UPF.out
Quantum ESPRESSO Course for Solid‐State Physics | 59
Input file: All input files are generated automatically by a bash
script file run.sh as
QE‑SSP/gr/pps/run.sh
1 #!/bin/bash
2 # Set a variable pp for 6 PPs.
3 for pp in C.pbe -hgh.UPF \
4 C.pz -hgh.UPF \
5 C.pbe -n-rrkjus_psl .0.1. UPF \
6 C.pz -n-rrkjus_psl .0.1. UPF \
7 C.pbe -n-kjpaw_psl .1.0.0. UPF \
8 C.pz -n-kjpaw_psl .0.1. UPF; do
9 # Make input file for the vc -relax calculation.
10 cat > vc -relax.$pp.in << EOF
11 & CONTROL
12 calculation = 'vc -relax '
13 pseudo_dir = '../ pseudo/'
14 outdir = '../tmp/'
15 prefix = 'gr '
16 etot_conv_thr = 1.0D-5
17 forc_conv_thr = 1.0D-4
18 /
19 &SYSTEM
20 ibrav = 4
21 a = 2.4639055825
22 c = 15.0
23 nat = 2
24 ntyp = 1
25 occupations = 'smearing '
26 smearing = 'mv '
27 degauss = 0.02
28 ecutwfc = 80
29 /
30 & ELECTRONS
31 mixing_beta = 0.7
32 conv_thr = 1.0D-9
33 /
34 &IONS
35 ion_dynamics = 'bfgs '
36 /
37 &CELL
38 cell_dynamics = 'bfgs '
39 press_conv_thr= 0.05
40 cell_dofree = '2Dxy '
41 /
42 ATOMIC_SPECIES
43 C 12.0107 $pp
44 ATOMIC_POSITIONS (crystal)
60 | Hands‐On Tutorials of Quantum ESPRESSO
where the index i runs over all states and fi are the occupation
numbers, which can be either 1 or 0 depending on that the i‐th state
is occupied or not, respectively, as shown in Fig. 3.9 (a). In Quantum
ESPRESSO, the sum of Eq. (3.3) is replaced by an integration over the
Brillouin zone (BZ) of the system as
(a) (b)
f fFD Determined by σ
1 1
0 0
EF E EF E
Figure 3.9 (a) The step function f and (b) the smooth Fermi‐Dirac function
fFD are plotted as a function of the energy E. The width of the Fermi‐Dirac
function is determined by σ.
1 Eik − EF
fFD = , with x = , (3.5)
exp(x) + 1 σ
1
fG = (3.6)
2
erfc(x),
in Fig. 3.10, the width of the Gaussian smearing is smaller than that
of the Fermi‐Dirac smearing. This means that in order to get similar
k‐points convergence as for the Fermi‐Dirac smearing, the Gaussian
smearing can use a smaller value of σ.
The Methfessel‑Paxton smearing: Methfessel and Paxton
proposed a more sophisticated approximation to the step function
based on the Gaussian smearing [Methfessel and Paxton (1989)].
The first‐order Methfessel‐Paxton approximation is expressed as
1 1
fMP = erfc(x) − √ x exp −x2 , (3.7)
2 2 π
1 2 1
fMV exp(−x − 2x − 1/2) + erfc x + √
2
. (3.8)
√
2 2
=
π
Marzari-Vanderbilt
Fermi-Dirac
f (x) Gaussian
Methfessel-Paxton
x
Figure 3.10 The several smearing functions f(x) are plotted as a function of
x = (Eik − EF )/σ.
1 $ cd ~/QE -SSP/gr/smearing/
2 $ ./run.sh &
$ ls
QE‑SSP/gr/smearing/run.sh
1 #!/bin/bash
2 # Set a variable sf for the smearing functions.
3 for sf in fd gauss mp mv; do
4
5 # Set a variable se for the smearing energies.
6 for se in 0.005 0.010 0.015 0.020 0.025 0.030 \
7 0.035 0.040; do
8
9 # Make input file for the scf calculation.
10 cat > $sf.$se.in << EOF
11 & CONTROL
12 calculation = 'scf '
13 pseudo_dir = '../ pseudo/'
14 outdir = '../tmp/'
15 prefix = 'gr '
16 /
17 &SYSTEM
18 ibrav = 4
19 a = 2.4639055825
20 c = 15.0
21 nat = 2
22 ntyp = 1
23 occupations = 'smearing '
24 smearing = '$sf '
25 degauss = $se
26 ecutwfc = 40
27 /
28 & ELECTRONS
66 | Hands‐On Tutorials of Quantum ESPRESSO
29 mixing_beta = 0.7
30 conv_thr = 1.0D-6
31 /
32 ATOMIC_SPECIES
33 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
34 ATOMIC_POSITIONS (crystal)
35 C 0.333333333 0.666666666 0.500000000
36 C 0.666666666 0.333333333 0.500000000
37 K_POINTS (automatic)
38 12 12 1 0 0 0
39 EOF
40
41 # Run pw.x for SCF calculation.
42 mpirun -np 4 pw.x <$sf.$se.in > $sf.$se.out
43
44 # Write the name of the smearing functions , the
smearing energies , and total energies
45 awk -v var="$sf" '/!/ {printf "%-6s %1.3f %s\n",var ,'
$se ',$5}' $sf.$se.out >> calc -smearing.dat
46 # End of for sf loop.
47 done
48 # End of for se loop.
49 done
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error in routine electrons (1):
charge is wrong: smearing is needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
$ vi calc -smearing.dat
QE‑SSP/gr/smearing/calc‑smearing.dat
1 fd 0.005 -23.90915383
2 fd 0.010 -23.90934651
3 fd 0.015 -23.90954805
4 ...
5 mv 0.030 -23.90923457
6 mv 0.035 -23.90927836
7 mv 0.040 -23.90931451
QE‑SSP/gr/smearing/plot‑smearing.ipynb
20 se[i]. append(float(tmp1))
21 ener[i]. append(float(tmp2))
22
23 # Create figure object
24 plt.figure ()
25 # Plot the data
26 plt.plot(se[0], ener [0],'o-', label='Fermi -Dirac ')
27 plt.plot(se[1], ener [1],'s-', label='Gaussian ')
28 plt.plot(se[2], ener [2],'v-', label='Methfessel -
Paxton ')
29 plt.plot(se[3], ener [3],'^-', label='Marzari -
Vanderbilt ')
30 # Add the legend
31 plt.legend(loc='lower left ')
32 # Add the x and y-axis labels
33 plt.xlabel('Smearing energy (Ry)')
34 plt.ylabel('Total energy (Ry)')
35 # Set the axis limits
36 plt.xlim (0.003 , 0.042)
37 plt.ylim ( -23.914 , -23.908)
38 # Save the figure.
39 plt.savefig('plot -smearing.pdf')
40 # Show the figure
41 plt.show ()
Try It Yourself
1 $ cd ~/QE -SSP/gr/rho/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pp.x < pp.in > pp.out &
$ vi pp.in
QE‑SSP/gr/rho/pp.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'gr '
4 plot_num = 0
5 /
6 &PLOT
7 iflag = 3
8 output_format= 6
9 fileout = 'gr_rho.cube '
10 nx = 64
11 ny = 64
12 nz = 12
13 /
$ VESTA gr_rho.cube
1 $ cd ~/QE -SSP/gr/bands/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 bands.x < bands.in > bands.out &
Input file: The detail of the scf.in file is given in Sec. 3.1.1. For
input file nscf.in, the readers can use vi editor by typing:
$ vi nscf.in
QE‑SSP/gr/bands/nscf.in
1 & CONTROL
2 calculation = 'bands '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 2.4639055825
10 c = 15.0
11 nat = 2
12 ntyp = 1
13 nbnd = 16
14 occupations = 'smearing '
15 smearing = 'mv '
16 degauss = 0.020
17 ecutwfc = 40
18 /
19 & ELECTRONS
20 mixing_beta = 0.7
21 conv_thr = 1.0D-6
22 /
23 ATOMIC_SPECIES
24 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
25 ATOMIC_POSITIONS (crystal)
26 C 0.333333333 0.666666666 0.500000000
27 C 0.666666666 0.333333333 0.500000000
28 K_POINTS (crystal_b)
29 4
30 gG 40
31 K 20
32 M 30
33 gG 0
ibrav = 3 b3 ibrav = 4
b3
P H A L
Γ b2
H
Γ K
b1 N b2 M
b1
SYSTEM, which means that 16 bands are calculated for the case of
graphene.
We also need to specify the k‐points in the card
K_POINTS (crystal_b) (line 28), in which crystal_b allows
us to use the labels of high‐symmetry points in the Brillouin‐
zone. In Fig. 3.13, we show the labels for several structures with
ibrav = 1, 2, 3, and 4. Note that the Γ‐point is denoted by gG
(line 30) in the input file. The labels of all structures can be found
in Ref. [Setyawan and Curtarolo (2010)]. It is noted that the k‐point
path must be continuous in the Brillouin‐zone.
To see the input file bands.in, the readers can type:
$ vi bands.in
QE‑SSP/gr/bands/bands.in
1 &BANDS
2 outdir = '../tmp/'
3 prefix = 'gr '
4 filband = 'gr.bands '
5 /
76 | Hands‐On Tutorials of Quantum ESPRESSO
Note: outdir and prefix must be the same as for scf.in for
pw.x calculation. By setting lp = .true. in the namelist &BANDS,
the readers can obtain the square of the absolute value of the
matrix elements of the momentum operator between valence and
conduction bands in the output file p_avg.dat.
Output file: The band structure of graphene is given by gr.bands,
which is output of bands.x. The readers can use vi editor to see the
gr.bands file as follows:
$ vi gr.bands
QE‑SSP/gr/bands/gr.bands
QE‑SSP/gr/bands/plot‑bands.ipynb
20
10
Energy (eV)
0
10
20
K M
Figure 3.14 Electronic band structure of graphene. The blue lines plot the
calculated electron dispersion. Dirac point is the transition between the
valence and conduction bands at the K point and zero energy. The ARPES
intensity of graphene appears on the electron dispersion of the valence band,
which is reproduced with permission [Ohta et al. (2007)].
42 plt.show ()
Try It Yourself
1 $ cd ~/QE -SSP/gr/edos/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 dos.x < dos.in > dos.out &
$ vi nscf.in
QE‑SSP/gr/edos/nscf.in
1 & CONTROL
2 calculation = 'nscf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 2.4639055825
10 c = 15.0
11 nat = 2
12 ntyp = 1
13 nbnd = 16
14 occupations = 'tetrahedra '
15 ecutwfc = 40
16 /
17 & ELECTRONS
18 mixing_beta = 0.7
19 conv_thr = 1.0D-6
20 /
21 ATOMIC_SPECIES
22 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
23 ATOMIC_POSITIONS (crystal)
24 C 0.333333333 0.666666666 0.500000000
25 C 0.666666666 0.333333333 0.500000000
26 K_POINTS (automatic)
27 48 48 1 0 0 0
$ vi dos.in
Quantum ESPRESSO Course for Solid‐State Physics | 81
QE‑SSP/gr/edos/dos.in
1 &DOS
2 outdir = '../tmp/'
3 prefix = 'gr '
4 fildos = 'gr.dos '
5 /
Note: outdir and prefix must be the same as for scf.in for
pw.x calculation.
Output file: The output file for the DOS of graphene is given by
gr.dos, which is output of step (3). The readers can use vi editor to
see the gr.dos file as follows:
$ vi gr.dos
QE‑SSP/gr/edos/gr.dos
3'552ITGFQURNQVFQUKR[PD
Try It Yourself
Valence configuration:
nl pn l occ Rcut Rcut US E pseu
2S 1 0 2.00 1.000 1.300 -1.010676
2P 2 1 2.00 1.000 1.450 -0.388489
1 $ cd ~/QE -SSP/gr/pdos/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 projwfc.x < projwfc.in > projwfc.out
&
84 | Hands‐On Tutorials of Quantum ESPRESSO
$ vi projwfc.in
QE‑SSP/gr/pdos/projwfc.in
1 &projwfc
2 prefix = 'gr '
3 outdir = '../tmp/'
4 /
QE‑SSP/gr/pdos/gr.pdos_atm#1(C)_wfc#1(s)
QE‑SSP/gr/edos/plot‑dos.ipynb
DOS
PDOS of 2p orbital
PDOS of 2s orbital
In Fig. 3.16, we show the DOS and the PDOS of the 2s and 2p
orbitals as functions of energy E. For the valence bands (E < 0),
the<DOS is equal to the total PDOS since all valence bands mainly
come from the 2s and 2p orbitals. For the conduction bands, > with
E > 2.5 eV, the DOS is not equal to the sum of PDOS because the
high energy bands include not only the 2s and 2p orbitals but also
more delocalized wavefunctions far from the atoms, which are called
interlayer‐bands [Fretigny et al. (1989)].
Quantum ESPRESSO Course for Solid‐State Physics | 87
Try It Yourself
D̃Iα,Jβ (q)ũJβ (q) = ω2q ũIα (q) with (Iα = 1, …, 3N), (3.9)
Jβ
where N is the number of atoms in the unit cell and α (β) corresponds
to x, y, z components of oscillation. ũIα (ũJβ ) is a phonon eigenvector
(or the amplitude of the phonon) of the I‐th (J‐th) atom. Summation
on Jβ is taken for neighbor atoms from a given I‐th atom up to several
nearest neighbors. D̃Iα,Jβ (q) is called the dynamical matrix, which is
3We use k‐points to refer to electron wavevectors and q‐points to refer to phonon
wavevectors
88 | Hands‐On Tutorials of Quantum ESPRESSO
defined by
1
D̃Iα,Jβ (q) = F̃Iα,Jβ (q), (3.10)
MI MJ
∂2 Etot
F̃Iα,Jβ (q) = . (3.11)
∂ũIα (q)∂ũJβ (q)
∗
Nq
1
FIα,Jβ (R) = F̃Iα,Jβ (qijk )eiqijk ·R , (3.12)
Nq
i=1
4 The DFPT is a method that can directly calculate the second‐order derivatives of the
Eqs. (3.12) and (3.14) allow the interpolation of the dynamical matrix
at arbitrary q′ , by a few IFCs. It is noted that the Fourier interpolation
works well if the IFCs decay rapidly in the real space. Therefore, the
Fourier interpolation might fail in some special cases when IFCs do
not decay. For example, when the Kohn anomaly occurs in a metal, the
dynamical matrix is not a smooth function of q, and the IFCs decay
slowly in the real space.
The calculation of the phonon dispersion of graphene has the
following four steps:
1. Perform the SCF with pw.x calculation.
2. Calculate the dynamical matrix in Eq. (3.10) on a uniform mesh of
q‐points by using ph.x.
3. Perform the inverse Fourier transform of the dynamical matrix in
Eq. (3.12) to obtain a set of the IFCs in the real space by using
q2r.x.
4. Perform the Fourier transform of the real space IFCs in Eq. (3.14)
to obtain the dynamical matrix at any q′ by using matdyn.x. This
allows us to calculate the phonon frequencies for a set of points
along lines between high‐symmetry points in the Brillouin zone,
which is similar to the electronic band structure in Sec. 3.2.2.
How to run: The readers can run the following command lines:
1 $ cd ~/QE -SSP/gr/phonon/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 ph.x < ph.in > ph.out &
4 $ mpirun -np 4 q2r.x < q2r.in > q2r.out &
5 $ mpirun -np 4 matdyn.x < matdyn.in > matdyn.out&
$ vi scf.in
QE‑SSP/gr/phonon/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 2.4639055825
10 c = 15.0
11 nat = 2
12 ntyp = 1
13 occupations = 'smearing '
14 smearing = 'mv '
15 degauss = 0.02
16 ecutwfc = 80
17 assume_isolated = '2D'
18 /
19 & ELECTRONS
20 mixing_beta = 0.7
21 conv_thr = 1.0D-6
22 /
23 ATOMIC_SPECIES
24 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
25 ATOMIC_POSITIONS (crystal)
26 C 0.333333333 0.666666666 0.500000000
27 C 0.666666666 0.333333333 0.500000000
28 K_POINTS (automatic)
29 12 12 1 0 0 0
$ vi ph.in
QE‑SSP/gr/phonon/ph.in
1 phonon calc.
2 & INPUTPH
3 outdir = '../tmp/'
4 prefix = 'gr '
5 tr2_ph = 1d-14
6 ldisp = .true.
7 nq1 = 6
8 nq2 = 6
9 nq3 = 1
10 fildyn = 'gr.dyn '
11 /
$ vi q2r.in
QE‑SSP/gr/phonon/q2r.in
1 &INPUT
2 fildyn = 'gr.dyn '
3 zasr = 'crystal '
92 | Hands‐On Tutorials of Quantum ESPRESSO
$ vi matdyn_disp.in
5 The ASR is a requirement for IFCs by translational and rotational invariance of the
1 &INPUT
2 asr = 'crystal '
3 flfrc = 'gr.fc '
4 flfrq = 'gr.freq '
5 flvec = 'gr.modes '
6 loto_2d = .true.
7 q_in_band_form = .true.
8 /
9 4
10 gG 40
11 K 20
12 M 30
13 gG 0
$ vi gr.freq
QE‑SSP/gr/phonon/gr.freq
1 &plot nbnd= 6, nks= 91 /
2 0.000000 0.000000 0.000000
3 -0.0000 -0.0000 -0.0000 876.7549 1488.7595 1488.7595
4 0.016667 0.000000 0.000000
5 7.4052 31.1526 49.6068 876.5584 1489.4337 1490.5041
6 ...
QE‑SSP/gr/phonon/plot‑phonon.ipynb
TO
ZO
LA
TA
ZA
7 The Kohn anomaly is an anomaly in the phonon dispersion of a metal. For a specific
Try It Yourself
1 $ cd ~/QE -SSP/gr/phdos/
2 $ cp ../ phonon/gr.fc ./
3 $ mpirun -np 4 matdyn.x < matdyn.in > matdyn.out &
$ vi matdyn.in
QE‑SSP/gr/phdos/matdyn.in
1 &INPUT
2 asr = 'crystal '
3 flfrc = 'gr.fc '
4 flfrq = 'gr.freq '
5 flvec = 'gr.modes '
6 loto_2d = .true.
7 fldos = 'gr.dos '
8 dos = .true.
9 nk1 = 48
10 nk2 = 48
11 nk3 = 1
12 /
$ vi gr.dos
QE‑SSP/gr/phdos/gr.dos
QE‑SSP/gr/edos/plot‑dos.ipynb
Figure 3.19 Phonon density of states (solid line) and phonon PDOS of a C
atom (dashed line) of graphene.
Try It Yourself
h̄
1/2
gνq (k, i, j) = (3.15)
2Mωνq
⟨ψi,k |∂νq VSCF |ψj,k+q ⟩,
where M is the total mass of the atoms in the unit cell, ωνq is the
phonon frequency, ∂νq VSCF is the derivative of the self‐consistent
potential VSCF by ionic displacement by phonon mode ν at q, and
ψi,k and ψj,k+q are the wavefunctions at initial state i and final state
j, respectively.
The phonon linewidth γ νq , which is the imaginary part of the
phonon self‐energy, is defined by gνq in Eq. (3.15) as
(3.16)
where ΩBZ is the volume of the Brillouin zone (BZ), ϵF is the Fermi
energy, and ϵi,k and ϵj,k+q are the energies of the initial state i and
final state j, respectively.
The electron‐phonon coupling constant λνq is defined by γ νq in
Eq. (3.16) as
, (3.17)
γ νq
λνq =
πh̄D(ϵF )ω2νq
1 $ cd ~/QE -SSP/gr/elph/
2 $ mpirun -np 4 pw.x < scf -dense.in > scf -dense.out &
3 $ mpirun -np 4 pw.x < scf.in > scf.out &
4 $ mpirun -np 4 ph.x < ph.in > ph.out &
5 $ mpirun -np 4 q2r.x < q2r.in > q2r.out &
6 $ mpirun -np 4 matdyn.x < matdyn.in > matdyn.out &
$ vi scf -dense.in
QE‑SSP/gr/elph/scf‑dense.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'gr '
6 /
Quantum ESPRESSO Course for Solid‐State Physics | 103
7 &SYSTEM
8 ibrav = 4
9 a = 2.4639055825
10 c = 15.0
11 nat = 2
12 ntyp = 1
13 occupations = 'smearing '
14 smearing = 'mv '
15 degauss = 0.020
16 ecutwfc = 80
17 assume_isolated = '2D'
18 la2F = .true.
19 /
20 & ELECTRONS
21 mixing_beta = 0.7
22 conv_thr = 1.0D-6
23 /
24 ATOMIC_SPECIES
25 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
26 ATOMIC_POSITIONS (crystal)
27 C 0.333333333 0.666666666 0.500000000
28 C 0.666666666 0.333333333 0.500000000
29 K_POINTS (automatic)
30 36 36 1 0 0 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
task # 3
from elphsum : error # 2
q is not a vector in the dense grid
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The input file scf.in is same as Sec. 3.3.1. For the input file
ph.in, the readers can open by typing:
104 | Hands‐On Tutorials of Quantum ESPRESSO
$ vi ph.in
QE‑SSP/gr/elph/ph.in
1 phonon calc.
2 & INPUTPH
3 outdir = '../tmp/'
4 prefix = 'gr '
5 tr2_ph = 1d-14
6 ldisp = .true.
7 nq1 = 6
8 nq2 = 6
9 nq3 = 1
10 fildyn = 'gr.dyn '
11 fildvscf = 'grdv '
12 electron_phonon = 'interpolated '
13 el_ph_nsigma = 10
14 el_ph_sigma = 0.02
15 /
$ vi elph_dir/elph.inp_lambda .1
QE‑SSP/gr/elph/elph_dir/elph.inp_lambda.1
$ vi elph.gamma .1
QE‑SSP/gr/elph/elph.gamma.1
QE‑SSP/gr/elph/plotband.in
1 elph.gamma .5
2 -30 300
3 elph.gamma .5. gnu
4 elph.gamma .5.ps
5 0.0
6 0.0 0.0
QE‑SSP/gr/elph/gamma.ipynb
1.00
0.50
TO
0.25
LO
0.00
Γ K M Γ
Figure 3.20 Phonon linewidth of graphene along the high symmetry‐lines of
the Brillouin zone.
Try It Yourself
1
α2 F(ω) =
(3.18)
dq
2 ν
ωνq λνq δ(ω − ωνq ).
BZ ΩBZ
α2 F(ω)
λ=2 λνq . (3.19)
dω =
ω νq
2
α F(ω) log(ω) , (3.21)
dω 2
ωln = exp
λ ω
1 $ cd ~/QE -SSP/gr/alpha/
2 $ cp -rf ../ elph/elph_dir/ ../ elph/gr.fc ./
3 $ mpirun -np 4 matdyn.x < matdyn.in > matdyn.out &
$ vi a2F.dos1
QE‑SSP/gr/alpha/a2F.dos1
QE‑SSP/gr/alpha/alpha.ipynb
For the parameters λ and ωln , they are given by the output file
lambda for each Gaussian broadening. The readers can open lambda
as follows:
$ vi lambda
QE‑SSP/gr/alpha/lambda
Try It Yourself
1 $ cd ~/QE -SSP/mos2/optic/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 epsilon.x < epsilon.in > epsilon.out
&
$ vi scf.in
QE‑SSP/mos2/optic/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'mos2 '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 3.1378055413
10 c = 20.0
11 nat = 3
12 ntyp = 2
13 ecutwfc = 80.0
14 /
15 & ELECTRONS
16 mixing_beta = 0.7
17 conv_thr = 1.0d-10
18 /
19 ATOMIC_SPECIES
20 Mo 95.94 Mo.pz -hgh.UPF
21 S 32.065 S.pz -hgh.UPF
22 ATOMIC_POSITIONS (crystal)
Quantum ESPRESSO Course for Solid‐State Physics | 115
$ vi nscf.in
116 | Hands‐On Tutorials of Quantum ESPRESSO
QE‑SSP/mos2/optic/nscf.in
1 & CONTROL
2 calculation = 'nscf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'mos2 '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 3.1378055413
10 c = 20.0
11 nat = 3
12 ntyp = 2
13 nbnd = 40
14 ecutwfc = 80.0
15 nosym = .true.
16 noinv = .true.
17 /
18 & ELECTRONS
19 mixing_beta = 0.7
20 conv_thr = 1.0d-10
21 /
22 ATOMIC_SPECIES
23 Mo 95.94 Mo.pz -hgh.UPF
24 S 32.065 S.pz -hgh.UPF
25 ATOMIC_POSITIONS (crystal)
26 Mo -0.0000000000 -0.0000000000 0.5000000000
27 S 0.3333333333 0.6666666667 0.4218571774
28 S 0.3333333333 0.6666666667 0.5781427066
29 K_POINTS (automatic)
30 24 24 1 0 0 0
9 Irreducible Brillouin zone is the first Brillouin zone, which is reduced by all
$ vi epsilon.in
QE‑SSP/mos2/optic/epsilon.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'mos2 '
4 calculation = 'eps '
5 /
6 & ENERGY_GRID
7 smeartype = 'gauss '
8 intersmear = 0.15
9 intrasmear = 0.0
10 wmin = 0.0
11 wmax = 10.0
12 nw = 500
13 shift = 0.0
14 /
$ vi epsr_mos2.dat
QE‑SSP/gr/optic/epsr_mos2.dat
QE‑SSP/mos2/optic/eps.ipynb
Nevertheless, peaks A and B show excitation energy. This is because the correction to
the transition energy from the exciton binding energy is largely offset by many‐body
corrections to the band gap in the DFT. Therefore, the predicted transition energy
within a single‐particle calculation is closer to the experiment than might be
expected [Li et al. (2014)].
Quantum ESPRESSO Course for Solid‐State Physics | 121
(a) (b)
C C
C
AB
AB
AB
Figure 3.23 In‐plane optical properties of monolayer MoS2 . (a) Real part
Re(εxx ) and imaginary part Im(εxx ) of the dielectric tensor are plotted by the
solid and dashed lines, respectively. (b) Absorption coefficient in x‐direction,
αxx=yy , of the monolayer MoS2 . It is noted that εxx = εyy and αxx = αyy . The
peaks labeled A and B correspond to transition energies at the K point of the
Brillouin zone, and C peak corresponds to higher‐lying interband transitions,
including the transitions between the K and Γ points.
Try It Yourself
1 $ cd ~/QE -SSP/mos2/optic/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 epsilon.x < epsilon -jdos.in > epsilon
-jdos.out &
$ vi epsilon -jdos.in
11 The optical absorption intensity is proportional to the product of JDOS and the
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'mos2 '
4 calculation = 'jdos '
5 /
6 & ENERGY_GRID
7 smeartype = 'gauss '
8 intersmear = 0.15
9 intrasmear = 0.0
10 wmin = 0.0
11 wmax = 10.0
12 nw = 500
13 shift = 0.0
14 /
QE‑SSP/mos2/optic/jdos.ipynb
Try It Yourself
1 $ cd ~/QE -SSP/mos2/raman/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 ph.x < ph.in > ph.out &
4 $ mpirun -np 4 dynmat.x < dynmat.in > dynmat.out &
$ vi ph.in
QE‑SSP/mos2/raman/ph.in
1 phonon calc.
2 & INPUTPH
3 outdir = '../tmp/'
4 prefix = 'mos2 '
5 fildyn = 'mos2.dmat '
6 tr2_ph = 1d-14
7 lraman = .true.
8 epsil = .true.
9 trans = .true.
10 asr = .true.
11 /
12 0.0 0.0 0.0
$ vi dynmat.in
QE‑SSP/mos2/raman/dynmat.in
1 &INPUT
2 fildyn = 'mos2.dmat '
3 asr = 'crystal '
4 /
$ vi dynmat.out
QE‑SSP/mos2/raman/dynmat.out
QE‑SSP/mos2/raman/raman.ipynb
Try It Yourself
1 $ cd ~/QE -SSP/mos2/soc/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 pw.x < nscf.in > nscf.out &
4 $ mpirun -np 4 bands.x < bands.in > bands.out &
$ vi scf.in
QE‑SSP/mos2/soc/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'mos2 '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 3.1825188839
10 c = 20.0
11 nat = 3
12 ntyp = 2
13 ecutwfc = 60.0
14 noncolin = .true.
15 lspinorb = .true.
16 /
17 & ELECTRONS
18 mixing_beta = 0.7
19 conv_thr = 1.0d-6
20 /
21 ATOMIC_SPECIES
22 Mo 95.94 Mo.rel -pbe -spn -rrkjus_psl .1.0.0. UPF
23 S 32.065 S.rel -pbe -nl -rrkjus_psl .1.0.0. UPF
24 ATOMIC_POSITIONS (crystal)
25 Mo 0.0000000000 0.0000000000 0.5000000000
26 S 0.3333333333 0.6666666667 0.4217548051
27 S 0.3333333333 0.6666666667 0.5782450789
28 K_POINTS (automatic)
29 6 6 1 0 0 0
$ vi nscf.in
QE‑SSP/mos2/soc/nscf.in
1 & CONTROL
2 calculation = 'bands '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'mos2 '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 3.1825188839
10 c = 20.0
11 nat = 3
12 ntyp = 2
13 nbnd = 60
14 ecutwfc = 60.0
15 noncolin = .true.
16 lspinorb = .true.
17 /
18 & ELECTRONS
19 mixing_beta = 0.7
20 conv_thr = 1.0d-6
21 /
22 ATOMIC_SPECIES
23 Mo 95.94 Mo.rel -pbe -spn -rrkjus_psl .1.0.0. UPF
24 S 32.065 S.rel -pbe -nl -rrkjus_psl .1.0.0. UPF
25 ATOMIC_POSITIONS (crystal)
26 Mo 0.0000000000 0.0000000000 0.5000000000
27 S 0.3333333333 0.6666666667 0.4217548051
28 S 0.3333333333 0.6666666667 0.5782450789
29 K_POINTS (crystal_b)
30 4
31 gG 40
32 K 20
33 M 30
34 gG 0
Δ = 0.15 eV
Figure 3.26 Electronic band structure of monolayer MoS2 (a) with SOC and
(b) without SOC.
$ vi bands.in
QE‑SSP/mos2/soc/bands.in
1 &BANDS
2 outdir = '../tmp/'
3 prefix = 'mos2 '
4 filband = 'mos2.bands '
5 /
Try It Yourself
Interlayer
distance
s6 CIJ6
EvdW = − fd (RIJ ), (3.24)
2
I̸=J
R6IJ
136 | Hands‐On Tutorials of Quantum ESPRESSO
1
fd (RIJ ) = , (3.25)
1+ e−d(RIJ /Rr −1)
1
Enl n(r)Φ(r, r′ )n(r′ )drdr′ ,
(3.26)
vdW [n(r)] =
2
$ ls
QE‑SSP/bi‑gr/vdw/run.sh
1 #!/bin/bash
2 # Set a variable vdw for 1 GGA and 6 vdW functionals
.
138 | Hands‐On Tutorials of Quantum ESPRESSO
3 for vdw in pbe vdw -df vdw -df2 vdw -df3 -opt1 \
4 vdw -df3 -opt2 vdw -df -C6 rvv10; do
5 # Make input file for the vc -relax calculation.
6 cat > vc -relax.$vdw.in << EOF
7 & CONTROL
8 calculation = 'vc -relax '
9 pseudo_dir = '../ pseudo/'
10 outdir = '../tmp/'
11 prefix = 'bi -gr '
12 etot_conv_thr = 1.0D-5
13 forc_conv_thr = 1.0D-4
14 /
15 &SYSTEM
16 ibrav = 4
17 a = 2.4639055825
18 c = 20.0
19 nat = 4
20 ntyp = 1
21 occupations = 'smearing '
22 smearing = 'mv '
23 degauss = 0.02
24 ecutwfc = 60
25 input_dft = '$vdw '
26 /
27 & ELECTRONS
28 mixing_beta = 0.7
29 conv_thr = 1.0D-9
30 /
31 &IONS
32 ion_dynamics = 'bfgs '
33 /
34 &CELL
35 cell_dynamics = 'bfgs '
36 press_conv_thr= 0.05
37 cell_dofree = '2Dxy '
38 /
39 ATOMIC_SPECIES
40 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
41 ATOMIC_POSITIONS (crystal)
42 C 0.000000000 0.000000000 0.412500000
43 C 0.333333333 0.666666666 0.412500000
44 C 0.000000000 0.000000000 0.587500000
45 C 0.666666666 0.333333333 0.587500000
46 K_POINTS (automatic)
47 8 8 1 0 0 0
48 EOF
49 # Run pw.x for vc -relax calculation.
50 mpirun -np 4 pw.x <vc -relax.$vdw.in > vc -relax.$vdw.
out
Quantum ESPRESSO Course for Solid‐State Physics | 139
51 # End of for loop.
52 done
Try It Yourself
Purpose: The bilayer graphene in Sec. 3.5.3 has a zero band gap.
However, the band gap can be tuned by applying an external electric
field perpendicular to the graphene layer, in which a band gap is
opened at the Dirac point. In this tutorial, we show how to obtain the
electrostatic potential and the band structure of the bilayer graphene
under the external electric field.
Background: In Quantum ESPRESSO, an external electric field E
can be applied for a 2D system by using a saw‐tooth potential Vext , as
shown in Fig. 3.29. Assuming that the 2D material is set to be in the
Quantum ESPRESSO Course for Solid‐State Physics | 141
2D material
E* E
Vext(z)
Vacuum Vacuum
c Z
emaxpos ℓ
emaxpos + eopreg
Figure 3.29 Schematics of saw‐tooth potential. The dark area represents the
2D material and the light gray spaces are the vacuums in the unit cell. The
solid lines with labels E and E∗ denote the linear external electric potential
and its counterpart, respectively. c is the size of the periodic cell along z
direction, and ℓ is the length of external electric potential. emaxpos and
eopreg are input parameters for Quantum ESPRESSO.
E∗ = −E , (3.28)
c−ℓ
ℓ
where ℓ and c are the length of the external electric field E and the
size of the unit cell in the z direction, respectively.
We note that the unit cell of the 2D material contains a dipole
moment in the vacuum space. In particular, when the bottom and top
layers of the 2D material are not equivalent to each other, there is a
142 | Hands‐On Tutorials of Quantum ESPRESSO
z 1
Vdip (z) = 4πm , (3.29)
c 2
−
$ vi scf.in
Quantum ESPRESSO Course for Solid‐State Physics | 143
QE‑SSP/bi‑gr/elec‑field/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'bi -gr '
6 tefield = .true.
7 dipfield = .true.
8 /
9 &SYSTEM
10 ibrav = 4
11 a = 2.4857910097
12 c = 20.0
13 nat = 4
14 ntyp = 1
15 occupations = 'smearing '
16 smearing = 'mv '
17 degauss = 0.02
18 ecutwfc = 80
19 input_dft = 'vdw -df2 '
20 edir = 3
21 emaxpos = 0
22 eopreg = 0.05
23 eamp = 0.0097234527
24 /
25 & ELECTRONS
26 mixing_beta = 0.7
27 conv_thr = 1.0D-9
28 /
29 ATOMIC_SPECIES
30 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
31 ATOMIC_POSITIONS (crystal)
32 C 0.0000000000 0.0000000000 0.4115083725
33 C 0.3333333333 0.6666666667 0.4114773543
34 C 0.0000000000 0.0000000000 0.5884916275
35 C 0.6666666667 0.3333333333 0.5885226457
36 K_POINTS (automatic)
37 8 8 1 0 0 0
in which the positions of the C atoms are located at the center of the
unit cell. Thus, as shown in Fig. 3.29, the position of the maximum
of the saw‐tooth potential will be set to 0 by emaxpos = 0 (line
21). We set the width of the window, where the saw‐tooth potential
decreases, about 1 Å by eopreg = 0.05 (line 22). We note that
emaxpos and eopreg have the unit of c = 20 Å. The magnitude of the
electric field is set to 0.5 V/Å by eamp = 0.0097234527 (line 23), in
which eamp has the unit of a.u. (1 V/Å= 0.0194469054 a.u.).
Note: For the external electric field, the readers should use a
large value of ecutwfc (line 18). For example, if the readers use
ecutwfc = 60, the readers might get an error in scf.out as follows:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error in routine electrons (1):
charge is wrong
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
For the input file pp.in, the readers can open as follows:
$ vi pp.in
QE‑SSP/bi‑gr/elec‑field/pp.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'bi -gr '
4 filplot = 'bi -gr.dat '
5 plot_num = 11
6 /
$ vi average.in
Quantum ESPRESSO Course for Solid‐State Physics | 145
QE‑SSP/bi‑gr/elec‑field/average.in
1 1
2 bi -gr.dat
3 1.0
4 300
5 3
6 3.0
$ vi nscf.in
QE‑SSP/bi‑gr/elec‑field/nscf.in
1 & CONTROL
2 calculation = 'bands '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'bi -gr '
6 tefield = .true.
7 dipfield = .true.
8 /
9 &SYSTEM
10 ibrav = 4
11 a = 2.4857910097
12 c = 20.0
13 nat = 4
14 ntyp = 1
146 | Hands‐On Tutorials of Quantum ESPRESSO
15 nbnd = 30
16 occupations = 'smearing '
17 smearing = 'mv '
18 degauss = 0.02
19 ecutwfc = 80
20 input_dft = 'vdw -df2 '
21 edir = 3
22 emaxpos = 0
23 eopreg = 0.05
24 eamp = 0.0097234527
25 /
26 & ELECTRONS
27 mixing_beta = 0.7
28 conv_thr = 1.0D-9
29 /
30 ATOMIC_SPECIES
31 C 12.0107 C.pbe -n-rrkjus_psl .0.1. UPF
32 ATOMIC_POSITIONS (crystal)
33 C 0.0000000000 0.0000000000 0.4115083725
34 C 0.3333333333 0.6666666667 0.4114773543
35 C 0.0000000000 0.0000000000 0.5884916275
36 C 0.6666666667 0.3333333333 0.5885226457
37 K_POINTS (crystal_b)
38 4
39 gG 40
40 K 20
41 M 30
42 gG 0
For the input file nscf.in, the readers can open as follows:
$ vi bands.in
QE‑SSP/bi‑gr/elec‑field/bands.in
1 &BANDS
2 outdir = '../tmp/'
3 prefix = 'bi -gr '
4 filband = 'bi -gr.bands '
5 /
$ vi avg.dat
QE‑SSP/bi‑gr/elec‑field/avg.dat
QE‑SSP/bi‑gr/elec‑field/plot‑potential.ipynb
0.3 eV
Try It Yourself
1 $ cd ~/QE -SSP/gr/w90/
2 $ mpirun -np 4 pw.x < scf.in > scf.out &
3 $ mpirun -np 4 open_grid.x < open_grid.in >
open_grid.out &
4 $ wannier90.x -pp gr &
5 $ mpirun -np 4 pw2wannier90.x < pw2wan.in > pw2wan.
out &
6 $ wannier90.x gr &
Input file: The input file scf.in is similar to scf.in in Sec. 3.1.1
for graphene, but we add nbnd = 16 in the namelist SYSTEM. For
the input file open_grid.in, the readers can open with vi editor as
follows:
$ vi open_grid.in
QE‑SSP/gr/w90/open_grid.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'gr '
$ vi gr.win
Quantum ESPRESSO Course for Solid‐State Physics | 153
QE‑SSP/gr/w90/gr.win
1 # CONTROL
2 num_bands = 16
3 num_wann = 5
4
5 dis_win_min = -22
6 dis_win_max = 11
7 dis_froz_min = -22
8 dis_froz_max = 0.5
9
10 num_iter = 20
11
12 # TIGHT -BINDING
13 write_hr = .true.
14
15 # PLOTTING
16 wannier_plot = .true.
17 bands_plot = .true.
18
19 wannier_plot_supercell = 3
20
21 begin kpoint_path
22 G 0.0000 0.0000 0.0000 K 0.3333 0.3333 0.0000
23 K 0.3333 0.3333 0.0000 M 0.5000 0.0000 0.0000
24 M 0.5000 0.0000 0.0000 G 0.0000 0.0000 0.0000
25 end kpoint_path
26
27 # SYSTEM
28 begin unit_cell_cart
29 ang
30 2.463906 0.000000 0.000000
31 -1.231953 2.133805 0.000000
32 0.000000 0.000000 15.000000
33 end unit_cell_cart
34
35 begin atoms_frac
36 C 0.333333333 0.666666666 0.500000000
37 C 0.666666666 0.333333333 0.500000000
38 end atoms_frac
39
40 # PROJECTIONS
41 guiding_centres = .true.
42
43 begin projections
44 C: pz
45 f= 0.5, 0.5, 0.5: s
46 f= 0.5, 0.0, 0.5: s
47 f= 0.0, 0.5, 0.5: s
154 | Hands‐On Tutorials of Quantum ESPRESSO
20
Energy (eV) 10
dis_win
0
froz_win
10
20
K M
Figure 3.32 An optimized selection for the frozen froz_win and the
disentanglement dis_win for the band structure of graphene. Solid lines
are the energy bands that are calculated by the DFT. Solid lines with circles
are “fitted energy bands” to the WFs, while only solid lines denote “not‐
fitted energy bands”. froz_win should be taken as large as possible with the
condition that we should not include “not‐fitted bands” (only solid lines).
dis_win should be taken as small as possible with the condition that we
should include all “fitted bands” (lines with circles).
48 end projections
49
50 # KPOINTS
51 mp_grid = 12 12 1
52
53 begin kpoints
54 0.00000000 0.00000000 0.00000000 0.0069444
55 0.00000000 0.08333333 0.00000000 0.0069444
56 ...
57 -0.083333333 -0.08333333 0.00000000 0.0069444
58 end kpoints
Continued
156 | Hands‐On Tutorials of Quantum ESPRESSO
Exiting .......
dis_windows: Energy window contains fewer states
than number of target WFs
$ vi pw2wan.in
QE‑SSP/gr/w90/pw2wan.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'gr_open '
4 seedname = 'gr '
5 write_unk = .true.
6 /
Figure 3.33 The Wannier functions of graphene includes pz (left) and s (right)
orbitals.
5). This option also produces a set of wavefunctions unk (r) on a real‐
space grid in the output files UNK00001.1, UNK00002.1, …
Note: If the k‐points in gr.win and open_grid.out are not the
same or prefix in pw2wan.in is not correct, the readers will get an
error in pw2wan.out as follows:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Error in routine pw2wannier90 (144):
Wrong number of k-points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Output file: The main output files for this tutorial are gr_*.xsf,
gr_band.dat, and gr_hr.dat. We will visualize the WFs from
gr_*.xsf, plot the band structure from gr_band.dat, and extract the
tight‐binding parameters from gr_hr.dat.
Visualizing the Wannier functions: The gr_*.xsf file can open
directly by using VESTA. For XCrySDen, we can open as File →
Open Structure → Open XSF (XCrySDen Structure File) and select
gr_*.xsf file. Then, we select Tools → Data Grid, click on the OK
button, enter a value in the box Isovalue (2 is a good choice for
graphene), click on the Submit button. In Fig. 3.33 (a) and (b), we
show the pz and s orbitals for the output files gr_00001.xsf and
gr_00003.xsf, respectively.
Plotting energy dispersion: By running step (5), the band
structure based on the Wannier interpolation (see Sec. 5.14.3) is
written in the output file gr_band.dat. In a similar way to plot
band structure of graphene in Sec. 3.2.2, the readers can run the
file plot-bands.ipynb to plot the band structure from output file
gr_band.dat. The band structure is plotted in Fig. 3.32, in which
160 | Hands‐On Tutorials of Quantum ESPRESSO
$ vi gr_hr.dat
...
0 0 0 5 1 -0.000004 0.000014
0 0 0 1 2 -2.910519 0.000006
0 0 0 2 2 -1.391550 -0.000000
...
0 0 0 5 5 -12.557218 0.000000
0 1 0 1 1 0.224064 -0.000001
0 1 0 2 1 0.019864 0.000001
...
Ty
Tx Tz n n’
Translation vector WF Real & imaginary parts
of unit cell indices of hopping parameters
the open dots and solid lines denote the band structure obtained
from wannier90.x and pw.x, respectively. For five lowest energy‐
bands, which correspond to three s‐ and two pz ‐orbitals, the Wannier
interpolation can reproduce the band structure from SCF calculation.
The advantage of the Wannier interpolation is that the band structure
can be accurately calculated for a dense k‐mesh based on the SCF
calculation for a coarse k‐mesh. This is helpful for the case that the
SCF calculation is time‐consuming for a dense k‐mesh, such as the
hybrid functional calculation, which will be explained in Sec. 3.6.2.
Extracting tight‑binding parameters: The tight‐binding param‐
eters of graphene are given by the output file gr_hr.dat, as shown
in Fig. 3.34. The readers can use the vi editor to open this file.
Each line in gr_hr.dat includes 7 values as follows: the first three
values give the translation vector T of the unit cell in x, y, and z
directions, respectively, the next two values give the WF indices n
and n′ , and the least two values are the real and imaginary parts of the
hopping parameters in units of eV tnn′ = ⟨wn,0 |H|wn′ ,T ⟩ (see Eq. 5.80).
For graphene, we consider the hopping parameters up to the 1st
nearest neighbor unit cell, i.e., t12 and t11 for T = (0, 0, 0) and (0, 1, 0),
respectively. Then, we obtain t12 = −2.91 eV and t11 = 0.224 eV, as
shown in Fig. 3.34.
Quantum ESPRESSO Course for Solid‐State Physics | 161
Try It Yourself
1 $ cd ~/QE -SSP/mos2/w90/
2 $ mpirun -np 8 pw.x < scf.in > scf.out &
3 $ mpirun -np 8 open_grid.x < open_grid.in >
open_grid.out &
4 $ wannier90.x -pp mos2 &
5 $ mpirun -np 8 pw2wannier90.x < pw2wan.in > pw2wan.
out &
6 $ wannier90.x mos2 &
$ vi scf.in
QE‑SSP/mos2/w90/scf.in
1 & CONTROL
2 calculation = 'scf '
3 pseudo_dir = '../ pseudo/'
4 outdir = '../tmp/'
5 prefix = 'mos2 '
6 /
7 &SYSTEM
8 ibrav = 4
9 a = 3.1825188839
10 c = 20.0
11 nat = 3
12 ntyp = 2
13 nbnd = 30
14 ecutwfc = 60.0
15 input_dft = 'hse '
16 exx_fraction = 0.25
17 nqx1 = 2
18 nqx2 = 2
19 nqx3 = 1
20 /
21 & ELECTRONS
22 mixing_beta = 0.7
23 conv_thr = 1.0d-6
24 /
25 ATOMIC_SPECIES
26 Mo 95.94 Mo.pbe -spn -rrkjus_psl .1.0.0. UPF
27 S 32.065 S.pbe -nl -rrkjus_psl .1.0.0. UPF
28 ATOMIC_POSITIONS (crystal)
29 Mo 0.0000000000 0.0000000000 0.5000000000
30 S 0.3333333333 0.6666666667 0.4217548051
31 S 0.3333333333 0.6666666667 0.5782450789
32 K_POINTS (automatic)
33 6 6 1 0 0 0
Quantum ESPRESSO Course for Solid‐State Physics | 163
Explanation of scf.in: The HSE functional is calculated by
setting input_dft = 'hse' (line 15) in the namelist SYSTEM.
exx_fraction (line 16) is the mixing coefficient, which is 0.25 for
the HSE functional, as shown in Table 4.3. nqx1, nqx2, and nqx3
(lines 17, 18, and 19), respectively, are dimensions of the q‐grid for
the Hartree‐Fock exchange interaction. The readers should check
the convergence of the energy band by changing the q‐grid. The
calculation of the q‐grid convergence will be very time‐consuming
compared to the k‐grid convergence in Sec. 3.1.3. Therefore, we
recommend increasing the number of processes in parallel, which
might require a workstation. Here, we select a coarse 2 × 2 × 1 grid
for q‐points, which is appropriate to run with a PC with 8 processes,
such as Intel Core i7 or i9.
For the input file open_grid.in, the readers can open as follows:
$ vi open_grid.in
QE‑SSP/mos2/w90/open_grid.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'mos2 '
4 /
For the input file mos2.win, the readers can open as follows:
$ vi mos2.win
QE‑SSP/mos2/w90/mos2.win
1 # CONTROL
2 num_bands = 30
3 num_wann = 11
4
5 num_iter = 20
6
7 dis_win_min = -8
8 dis_win_max = 5
9 dis_froz_min = -8
10 dis_froz_max = 3.5
164 | Hands‐On Tutorials of Quantum ESPRESSO
11
12 # PLOTTING
13 bands_plot = .true.
14
15 begin kpoint_path
16 G 0.0000 0.0000 0.0000 K 0.3333 0.3333 0.0000
17 K 0.3333 0.3333 0.0000 M 0.5000 0.0000 0.0000
18 M 0.5000 0.0000 0.0000 G 0.0000 0.0000 0.0000
19 end kpoint_path
20
21 #SYSTEM
22 begin unit_cell_cart
23 ang
24 3.182518 0.000000 0.000000
25 -1.591259 2.756142 0.000000
26 0.000000 0.000000 20.000000
27 end unit_cell_cart
28
29 begin atoms_frac
30 Mo 0.0000000000 0.0000000000 0.5000000000
31 S 0.3333333333 0.6666666667 0.4217548051
32 S 0.3333333333 0.6666666667 0.5782450789
33 end atoms_frac
34
35 # PROJECTIONS
36 guiding_centres = .true.
37
38 begin projections
39 Mo: dxy; dyz; dxz; dx2 -y2; dz2
40 S: px; py; pz
41 end projections
42
43 # KPOINTS
44 mp_grid = 6 6 1
45
46 begin kpoints
47 0.00000000 0.00000000 0.00000000 0.0277777
48 0.00000000 0.16666666 0.00000000 0.0277777
49 ...
50 -0.16666666 -0.16666666 0.00000000 0.0277777
51 end kpoints
Figure 3.35 Energy band structure of monolayer MoS2 . Solid and dashed
lines denote the band structure with the GGA and HSE functionals,
respectively.
For the input file pw2wan.in, the readers can open as follows:
$ vi pw2wan.in
QE‑SSP/mos2/w90/pw2wan.in
1 & INPUTPP
2 outdir = '../tmp/'
3 prefix = 'mos2_open '
4 seedname = 'mos2 '
5 /
Try It Yourself
1. Calculate the energy gap of the bulk Si with the LDA, GGA,
and HSE functionals, and compare to the experimental
value (1.12 eV).
2. Plot the Wannier function of the bulk Si with the HSE
functional.
Chapter 4
Density‐Functional Theory
(a)
Quantum
Input Output
ESPRESSO
Structures Properties
of materials of materials
HΨ = Etot Ψ, (4.1)
where H is the Hamiltonian of the system, Etot is the total energy of the
electrons and nuclei, and Ψ(r1 , …, rNe , R1 , …, RNn ) is the wavefunction
of many particles (electrons and nuclei). By using the atomic units
listed in Table 4.1, a general Hamiltonian in Eq. (4.1) is given by
kinetic energies T and potential energies V of the electrons and
nuclei as
H = Tn + Vn + Te + Ve + Ven
ZI ZJ
Nn Nn
∇2RI 1
2MI 2 |RI − RJ |
=− +
I=1 I̸=J
(4.2)
nuclei
Ne Ne Ne Nn
∇2r 1 1 −ZI
.
i
2 2 |ri − rj | |ri − RI |
− + +
i=1 i̸=j i=1 I=1
electrons mixed
is expressed by
h̄2 ∇2ri e2
Ne Ne
1
. (4.3)
2me 2 4πϵ0 |ri − rj |
− +
i=1 i̸=j
e–
r1 |r1 – r2|
e–
r2
2e+
x y
1 1 1 2 2
− ∇2r1 − ∇2r2 + ψ(r1 , r2 )
2 2 |r1 − r2 | |r1 | |r2 | (4.5)
− −
= Eψ(r1 , r2 ).
Ne Ne
Nn
∇2r
ψ = Eψ.
−ZI
(4.6)
i
|ri − RI |
(Te + Ven )ψ =
2
− +
i=1 i=1 I=1
N
∇2r n
h(r) = −
−ZI
, (4.7)
2 |r − RI |
+
I=1
From Eq.(4.8) we can say that, for the case of the non‐interacting
electrons, the Hamiltonian of the system can be written as Ne
independent single‐particle Hamiltonian, and the wavefunction
is written by ψ(ri ). Eq. (4.8) is used to obtain the energy of
electrons for the following two cases: (1) distinguishable4 and (2)
indistinguishable5 electrons. For simplicity, we consider the two
electrons in a helium atom as follows:
(1) Two distinguishable electrons: Let us consider electron 1
in one state ϕa (r1 ) and electron 2 in another state ϕb (r2 ). Since two
electrons are independent of each other, the probability of finding
(electron 1 at r1 ) and (electron 2 at r2 ) can be given by the product
as |ψ(r1 , r2 )|2 = |ϕa (r1 )|2 |ϕb (r2 )|2 . In this case, the wavefunction can
4 If the electrons are distinguishable, the physical properties of the system are
[h(r1 ) + h(r2 )] ϕa (r1 )ϕb (r2 ) = Eϕa (r1 )ϕb (r2 ). (4.9)
where ϵα is the energy of one electron, the left‐hand side of Eq. (4.9)
is rewritten as
E = ϵa + ϵb . (4.12)
Thus, for the case of the non‐interacting system with the distinguish‐
able electrons, the total energy of the electrons is the sum of the
energy of each electron.
(2) Two indistinguishable electrons: In quantum mechanics,
two electrons are indistinguishable (i.e., ϕa (r1 )ϕb (r2 ) is equivalent to
ϕa (r2 )ϕb (r1 )). Thus, one possible shape of the wavefunction ψ(r1 , r2 )
is a linear combination of two states as [Kittel (1976)]
1
ψ(r1 , r2 ) = √ [ϕa (r1 )ϕb (r2 ) − ϕb (r1 )ϕa (r2 )], (4.13)
2
where P1 (r) = |ψ(r, r2 )|2 dr2 and P2 (r) = |ψ(r1 , r)|2 dr1 are prob‐
By substituting Eq. (4.13) into Eq. (4.16) and using the normalization
and orthogonality conditions, respectively, as
ϕ∗a (r)ϕa (r)dr = ϕ∗b (r)ϕb (r)dr = 1 (4.17)
and
ϕ∗b (r)ϕa (r)dr = ϕ∗a (r)ϕb (r)dr = 0, (4.18)
Z3/2 Z2
ϕ(r) = √ exp(−Z|r|) and ϵ = − , (4.20)
π 2
2Z3
n(r) = exp(−2Z|r|) and E = −Z2 , (4.21)
π
16
n(r) = exp(−4|r|) (4.22)
π
He+ + e− −
I I
He −
→1
→2
He++ + 2e− , (4.23)
1
h(r1 ) + h(r2 ) + ϕ (r1 )ϕb (r2 ) = Eϕa (r1 )ϕb (r2 ). (4.24)
|r1 − r2 | a
with
ϵHb = E − ϵa , (4.28)
n(r′ )
dr′ . (4.33)
|r − r′ |
VH (r) =
occupied states. We repeat the procedure until n(r) does not change
180 | Density‐Functional Theory
No End
Figure 4.4 Flow‐chart of the self‐consistent field method for the solution of
the single‐particle Schrödinger equation with the Hartree potential VH .
1 ∂ ∂
2
r 2
r ∂r ∂r
∇ = 2
(4.34)
1 ∂ ∂ 1 ∂2
.
r2 sin θ ∂θ ∂θ r2 sin θ ∂φ2
+ sin θ + 2
1 ∂ 2 ∂VH
r (4.35)
r2 ∂r ∂r
= −4πn(r).
r
1 r′
C
r′′ n(r′′ )dr′ dr′′ −
. (4.36)
2
r′2 r
VH (r) = −4π
0 0
It is noted that C/r in Eq. (4.36) is a solution of Eq. (4.35) for any value
of C. Then by substituting the electron density of the helium atom in
Eq. (4.21) into Eq. (4.36) and using
x
1 x′
2
x′′ exp(x′′ )dx′ dx′′ = (1 −
(4.37)
2
x′2 x
) exp(x),
0 0
1 C
VH (r) = −2 Z + . (4.38)
r r
exp(−2Zr) −
1
r n(r′ )dr′ = Zr . Applying this boundary condition for Eq. (4.38), we
obtain C = −Z, and Eq. (4.38) can be rewritten for the helium atom
182 | Density‐Functional Theory
Figure 4.5 The Hartree potential VH (r) (solid line), the Coulomb attraction
between electrons and nuclei Ven (dotted line), and the total potential Ven +
VH (r) (dashed line) of helium atom are plotted in atomic unit as functions of
r.
as
1 Z
VH (r) = −2 Z + exp(−2Zr) + . (4.39)
r r
1
h(r1 ) + h(r2 ) +
|r1 − r2 |
[ϕa (r1 )ϕb (r2 ) − ϕb (r1 )ϕa (r2 )]
(4.40)
= E[ϕa (r1 )ϕb (r2 ) − ϕb (r1 )ϕa (r2 )].
By multiplying ϕ∗b (r2 ) in Eq. (4.40) and integrating over r2 with the
normalization and orthogonality conditions (Eqs. (4.17) and (4.18)),
we obtain:
where ϵHF
a is defined by
γ(r2 , r1 )
a ϕa (r1 ), (4.43)
ϕ (r2 )dr2 = ϵHF
|r1 − r2 | a
[h(r1 ) + VH (r1 )] ϕa (r1 ) −
γ(r2 , r1 ) = ϕ∗a (r2 )ϕa (r1 ) + ϕ∗b (r2 )ϕb (r1 ). (4.44)
γ(r1 , r2 )
b ϕb (r2 ), (4.45)
ϕ (r1 )dr1 = ϵHF
|r1 − r2 | b
[h(r2 ) + VH (r2 )] ϕb (r2 ) −
in which ϵHF
b is symmetrically given by
b = E − ϵa .
ϵHF (4.46)
γ(r2 , r1 )
Vx (r1 , r2 ) = − . (4.47)
|r1 − r2 |
Vx exists only for two electrons at r1 and r2 , which have the same
quantum state ϕa or ϕb including spin, i.e., the exchange potential
is the effective Coulomb interaction between the electrons with
parallel spins in the system. The origin of Vx is that two‐electrons
with the same spin can not exist at the same position (r1 = r2 ), in
which VH should be small by the restriction of the Pauli principle. By
considering both VH and Vx in Eq. (4.39), two electrons are interacted
not only by their electronic charge but also by their spins. Therefore,
the quantum behavior of electrons is included by Vx term in the
Hartree‐Fock equation. It is important to note that Vx is a non‑local
potential on r1 since it depends on an integration over the additional
variable r2 . Therefore, the practical solution of the Hartree‐Fock
equation becomes more complicated and time‐consuming.
Quantum ESPRESSO Course for Solid‐State Physics | 185
Let us discuss the energy that is obtained by the Hartree‐Fock
equation. From Eq. (4.43), we can write the energy by integrating the
wavefunctions as
ϵHF
a = ϕ∗a (r1 )[h(r1 ) + VH (r1 )]ϕa (r1 )dr1
γ(r2 , r1 ) (4.48)
ϕ (r2 )dr1 r2
|r1 − r2 | a
+ ϕ∗a (r1 )
=ϵa + J − K,
and
K=
dr1 dr2 . (4.50)
ϕ∗a (r1 )ϕb (r1 )ϕ∗b (r2 )ϕa (r2 )
|r1 − r2 |
E = ϵa + ϵb + J − K. (4.51)
186 | Density‐Functional Theory
On the other hand, from Eq. (4.40), the total electron energy can be
obtained, too, as
1
ψ∗ (r1 , r2 ) ψ(r1 , r2 )dr1 r2
(4.52)
|r1 − r2 |
+
|ψ(r2 , r1 )|2
dr1 r2 .
|r1 − r2 |
= ϵa + ϵb +
|ψ(r1 , r2 )|2
J−K=
dr1 dr2 ≥ 0. (4.53)
|r1 − r2 |
Total energy of helium atom: In the case of the helium atom, two
electrons at the 1s state ϕ must have different spin directions, up ↑
and down ↓ because of the Pauli principle. Let us consider the spin
states of two electrons, respectively, as
1 0
and ϕ↓ (r) = ϕ(r) , (4.54)
0 1
ϕ↑ (r) = ϕ(r)
1
ϕ∗↑ (r)ϕ↑ (r) = ϕ2 (r) 1 0 = ϕ2 (r), (4.55)
0
and
0
ϕ∗↑ (r)ϕ↓ (r) = ϕ2 (r) 1 0 = 0. (4.56)
1
Quantum ESPRESSO Course for Solid‐State Physics | 187
By substituting Eqs. (4.55) and (4.56) into Eqs. (4.49) and (4.50),
respectively, we obtain:
2 2
J=
drdr′ ,
ϕ∗↑ (r)ϕ↑ (r)ϕ∗↓ (r′ )ϕ↓ (r′ ) |ϕ(r)| |ϕ(r′ )|
|r − r′ |2 |r − r′ |2
drdr′ =
(4.57)
and
K=
drdr′ = 0. (4.58)
ϕ∗↑ (r)ϕ↓ (r)ϕ∗↓ (r′ )ϕ↑ (r′ )
|r − r′ |2
1 |ϕ(r)|2 n(r′ ) 1
J=
|ϕ(r)|2 VH (r)dr. (4.59)
2 |r − r′ |2 2
drdr′ =
1
E = 2ϵ +
|ϕ(r)|2 VH (r)dr. (4.60)
2
In order to apply Eq. (4.60) for the helium atom, we define an effective
atomic number Z, i.e., ϕ(r, Z), which gives minimized E when subject
to the Hartree potential VH (r). Eq. (4.60) can be obtained by the
following steps:
• Step 1: Using a beginning value Z = 2, from Eq. (4.39), the
Hartree potential VH (r) is given by
1 2
VH (r) = −2 2 + exp(−4r) + . (4.61)
r r
188 | Density‐Functional Theory
E = 2ϵ + J − K
Eexp = − 2.9 Ha
E = 2ϵ
6 The variational principle states that the energy of any approximate wavefunction is
higher than or equal to the energy of ground‐state wavefunction. The lower the
energy is the better the approximation to the ground state.
Quantum ESPRESSO Course for Solid‐State Physics | 189
could be done by recalculating the Hartree potential with new state
ϕ(r, Zeff = 1.826) in Step 1, then we recalculate Steps 2 and 3 until
the value of Zeff does not change significantly. Finally, we can obtain
the Eeff = −2.85 Ha, which is not too far from the exact value from the
experiment about −2.9 Ha. Slater derived a formula to determine Zeff
as Zeff = Z − S, where S is the shielding constant, which is known as
Slater’s rule [Slater (1930a)]. For the 1s state in the helium atom,
S = 0.3 and Zeff = 1.7, which is close to the result of the SCF method
(Zeff = 1.826).
In Sec. 4.6, we show that due to the Pauli exclusion principle, the
electrons with parallel spins have an exchange term of the Coulomb
interaction Vx . The two electrons with the parallel spin move to avoid
each other because the electrons can not overlap. Avoiding picture
should exist for the two electrons with anti‐parallel spins because
these electrons also move with keeping apart to lower the Coulomb
repulsion. This behavior of the electrons with anti‐parallel spins is
missing in the Hartree approximation as well as in the Hartree‐Fock
approximation. A complete picture of the motion of an electron in a
system is shown in Fig. 4.7. In general, the correlation interaction,
Vc , is defined by the remaining term that is missing in the Hartree‐
Fock method.7 Thus, an effective potential in the single‐particle
Schrödinger equation is defined by
∇2r
(4.64)
2
− + Veff [n(r)] ϕ(r) = ϵϕ(r).
Vx
VH Vc
Vc
When we know the value of Vxc [n(r)], Eq. (4.64) can be solved by
using the SCF method as discussed in Sec. 4.5. However, we do not
know the exact shape of Vxc [n(r)]. Many people proposed the many
functionals for Vxc [n(r)] over the past few decades. The efforts to
develop the accurate approximations for Vxc [n(r)] led to a theory,
the so‐called density‐functional theory (DFT). We will discuss these
approximations in the next section.
Let us consider how to obtain Vxc [n(r)] for a free‐electron gas, where
the contribution of ions is treated as a uniform‐positive‐charge
background, which is called the Jellium model [Brack (1993)], as
shown in Fig. 4.8. Since the potential is uniform, the electronic states
are expressed by plane waves as (see Sec. 5.3)
1
ϕk,σ (r) = √ exp(ikr)χ σ , (4.65)
V
Quantum ESPRESSO Course for Solid‐State Physics | 191
Positive ion charges Uniform background
+ + + + + + +
+
+ + + + + + + +
+ + + + + + + +
+ + + + + + +
+
Since the unit cell can not be defined by the uniform background,
the value of k is taken not from the Brillouin zone but all the k
space. Since ϕk,σ (r) is a one‐electron wavefunction, ϕk,σ (r) can be
a solution of the Hartree‐Fock equation in Eq. (4.43). Moreover, in
the free‐electron gas, the electron density of the ground state is
also a uniform‐negative‐charge background. Therefore, the electron
density cancels the positive charge background, i.e., Ven + VH = 0
and only the exchange term Vx survives in Eq. (4.43). Thus, the single‐
particle Schrödinger equation of a free‐electron gas can be written as
Ix = −
ϕ∗k′ ,σ′ (r′ )ϕk′ ,σ′ (r)
|r − r′ |
ϕk,σ (r′ )dr′
k
(4.67)
σ ′ ′
∗
ϕk′ ,σ′ (r′ )ϕk,σ (r′ ) ′
dr .
|r − r′ |
=− ϕk′ ,σ′ (r)
′ σ′ k
192 | Density‐Functional Theory
V3/2 |r − r′ |
′
dr′
k ′
f(k − k) ϕk,σ (r),
′
= −
k′
1 1 4π
f(q) =
, (4.69)
exp(−iqx)
V V q2
dx =
|x|
1 4π
Ix = −
V ′ |k′ − k|2
ϕk,σ (r)
k
1 4π
kF
(4.70)
(2π)3 0 |k′ − k|2
= − dk ϕk,σ (r)
′
kF k
=− F
kF
ϕk,σ (r),
π
1 − x2 1 + x k
F(x) = 1 + , with x = . (4.71)
2x 1 − x kF
ln
k = kF
kF k
− F
∇r
(4.72)
kF
ϕk,σ (r) = ϵk ϕk,σ (r).
2
−
π
k2 kF k
− F . (4.73)
kF
ϵk,σ
2
=
π
By taking the sum of all occupied states, we obtain the total energy
as
k2 1 kF
k
E= F
2 2 σ kF
−
σ |k| < kF
π
|k|<kF
2V k kF V k
kF 2 kF
F
(2π)3 >0 2 π (2π)3 0 kF
= dk − dk
V k kF V k
kF 4 kF
4π dk − 4 4π k F
2
(4.74)
4π3 2 8π k
= dk
0 0 F
V 5 V 4 1 2
k − k x F(x)dx
10π2 F 2π3 F 0
=
V 5 V
=1/2
kF − 3 k4F ,
10π 2 4π
=
194 | Density‐Functional Theory
with the factor of 1/2 in the second term of the right‐hand side in the
first line of Eq. (4.74) is needed to compensate for double counting
of the exchange interaction.
For the free‐electron gas, the number of electrons Ne is equal to
the number of occupied states in the systems:
kF
V kF
k3F
Ne = k2 dk = V
= 2V . (4.75)
dk
(2π)3 3π2
= 2
σ k 0 π 0
Ne k3
n= = F2 . (4.76)
V 3π
By substituting Eqs. (4.75) and (4.76) into Eq. (4.74), the total energy
per electron is expressed by
E
= C1 n2/3 + C2 n1/3 , (4.77)
Ne
3 3 3
1/3
C1 = (3π2 )2/3 = 2.871 and C2 = − = −0.738. (4.78)
10 4 π
The first term (C1 n2/3 ) in Eq. (4.77) represents the kinetic energy, and
the second term (C2 n1/3 ) represents the effective electron‐electron
interaction due to exchange potential. We can see that the exchange
potential reduces the total energy of the system from the kinetic
energy.
For a system with the free electron gas, n is a constant. However,
for a system with the non‐uniform‐electron charge, n becomes
a function of r. Slater [Slater (1951)] introduced the exchange
potential as a function of n(r) as
3 3
1/3
Vx [n(r)] = 2C2 n1/3 (r) = − n1/3 (r), (4.79)
2 π
where α usually taken in the range 2/3 < α < 1. For example, α =
0.772 and 0.978 [Schwarz (1972)] for He and H atoms, respectively.
Eq. (4.80) is called the Xα potential. > >
E Z 1 n(r′ )
= C1 n2/3 + C2 n1/3 −
dr′ , (4.81)
Ne |r| 2 |r − r′ |
+
where Z/|r| is the external potential that each electron feels the
presence of the ions and the last term is the Coulomb repulsion
between an electron at position r and all other electrons at r′ , which
is expressed by the density n(r′ ). The factor of 1/2 of the last term
196 | Density‐Functional Theory
Ne = n(r)dr.
(4.82)
By substituting Eq. (4.82) into Eq. (4.81), the total energy of the
system is expressed as
(4.83)
n(r) 1 n(r)n(r′ ) ′
−Z
2 |r − r′ |
dr + dr dr,
|r|
where the first and second terms of the right‐hand side are the
kinetic and exchange energies of the free‐electron gas, respectively,
the third term is the Coulomb attraction between the ions and the
electron density, and the last term describes Coulomb repulsion
between the electrons. n(r) for the ground‐state and the total
energy of the system can be obtained by minimizing the energy
functional E[n(r)]. Although the Thomas‐Fermi‐Dirac approximation
is not accurate good‐enough compared with the present electronic
structure calculation, the approach shows a prototype expression of
the DFT, in which n(r) is a crucial physical quantity to calculate the
ground‐state properties of a many‐electron system.
where F includes all the terms in the Hamiltonian except for the
external potential. That means that F contains the kinetic energy and
electron‐electron interaction terms. Therefore, F has the same shape
for all Ne ‐electrons systems with any external potentials.
The total energies of the ground states of the Hamiltonians is
given by
(r)] n(r)dr
⟨ψ |Ven −
′
Ven
′
|ψ′ ⟩ = [Ven (r) − Ven
′
(4.89)
and
(r)] n(r)dr.
⟨ψ|Ven − Ven
′
|ψ⟩ = [Ven (r) − Ven
′
(4.90)
potential Ven (r) and hence the Hamiltonian H. Because the ground‐
state wavefunction ψ is obtained by solving the Schrödinger equation
for H, ψ must be a unique functional of n(r). Therefore, we obtain the
following equation:
1 n(r)n(r′ )
drdr′ + Exc [n(r)], (4.94)
|r − r′ |
F[n(r)] = Ke [n(r)] +
2
i
200 | Density‐Functional Theory
∇2
Ke [n(r)] = ϕi − r ϕi , (4.96)
i
2
1 n(r)n(r′ )
E[n(r)] = Ke [n(r)] +
2 |r − r′ |
drdr′
(4.97)
+ Exc [n(r)] + Ven (r)n(r)dr,
∇2r
+ Veff [n(r)] ϕi (r) = ϵi ϕi (r), (4.99)
2
−
n(r′ )
δExc [n(r)]
, (4.100)
|r − r′ |
Veff [n(r)] = Ven (r) + dr′ +
δn(r)
where the first term Ven (r) is the external potential due to the ions,
the second term is the Hartree potential VH (see Eq. (4.33)) and the
Quantum ESPRESSO Course for Solid‐State Physics | 201
last term is the variational functional derivative of the exchange‐
correlation interaction Exc [n(r)]. The last term in Eq. (4.100) is
defined as the exchange‑correlation potential:
δExc [n(r)]
Vxc [n(r)] = . (4.101)
δn(r)
Multiplying the Kohn‐Sham equation in Eq. (4.99) by ϕ∗i (r) from the
left and summing over all occupied states, we obtain the Kohn‑Sham
energy as
(4.102)
ϵi = Ke [n(r)] + Veff (r)n(r)dr.
i
202 | Density‐Functional Theory
By substituting Veff (r) from Eq. (4.100) into Eq. (4.102), then
subtracting the total energy E[n(r)] in Eq. (4.97), we find:
1 n(r)n(r′ )
E[n(r)] =
|r − r′ |
ϵi −
2
drdr′ − ΔExc [n(r)]
i
(4.103)
1
ϵi − VH (r)n(r)dr − ΔExc [n(r)],
2
=
i
Note that the ratio of the exchange potential in Eq. (4.109) to the
Slater exchange potential in Eq. (4.79) is 2/3.
Since Eq. (4.107) is based on the free‐electron gas (see Sec. 4.8),
the correlation interaction is needed to capture accurately the many‐
body system. Therefore, the energy of correlation interaction ϵc (r) is
added to Eq. (4.107) as
LDA
Exc [n(r)] = [ϵx (r) + ϵc (r)]n(r)dr. (4.110)
4πr3s 1 3
1/3
, or rs = , (4.112)
3 n(r) 4πn(r)
=
ci c0 c1 c2
+ 3/2 + 2 + … (rs ≫ 1), (4.113)
∞
rs rs
ϵc (rs ) =
i=0 rs rs
i/2+1
=
n↑ (r) − n↓ (r)
, (4.114)
n(r)
ζ=
b
0)
X(x0 ) X(x)
3.72744 7.06042
ln (x−x +
VWN
c = 12.9352
=
18.0578
]}
2(b+2x0 )
arctan 2x+b , where
Q
Q x0 = −0.10498 −0.32500
x = rs , Q = 4c − b2 , and
√ √
X(x) = x + bx + c 2
....................................................................................................................................................................
A = 0.0311 0.0311
B = −0.048
A ln(rs ) + B + Crs ln(rs ) + Drs ,
C = 0.0020
−0.048
0.0007
where rs < 1 and
PZ D = −0.0116
γ/(1 + β1 rs + β2 rs ) with
√ −0.0048
rs ≥ 1
γ = −0.1423 −0.0843
β1 = 1.0529 1.3981
0.3334 0.2611
>
β 2 =
....................................................................................................................................................................
A = 0.031091 0.015545
α1 = 0.21370 0.20548
[
−2A(1 + α1 rs ) ln 1 + β1 = 7.5957 14.1189
PW ] β2 = 3.5876 6.1977
1 β3 = 1.6382 3.3662
β4 = 0.49294 0.6251)
p+1
2A β1 rs +β2 rs +β3 rs +β4 rs
1/2 3/2
p = 1.0 1.0
where n↑ (r) and n↓ (r) are the electron densities of up‐spin and
down‐spin states. ζ = 0 and 1 are the spin‐unpolarized and spin‐
polarized systems, respectively.
In Fig. 4.10, we show −ϵc by VWN, PZ, and PW parameters
as a function of rs for the uniform electron system with both
the spin‐unpolarized and spin‐polarized systems by using the
parameters in Table 4.2. All functions fit very well with the
numerical results (the circles and diamonds symbols for the spin‐
unpolarized and spin‐polarized cases, respectively) of Ceperley and
Alder [Ceperley and Alder (1980)]. For rs < 5, which is relevant for
condensed metals, the correlation energy becomes more important.
On the other hand, the correlation energy > becomes a small
206 | Density‐Functional Theory
WVN PZ PW
spin-unpolarized spin-unpolarized spin-unpolarized
spin-polarized spin-polarized spin-polarized
Fx (s) = 1 + κ − , (4.117)
κ
1 + μs2 /κ
where κ = 0.804 and μ = 0.21951 are constants. Note that the range
of s for the real systems is 0 ≤ s ≤ 3. Eqs. (4.115) and (4.117) show
that
1 + At2
H[n(r), ζ] = γϕ3 ln 1 + t2 . (4.120)
β
γ 1 + At2 + A2 t4
9 The repulsive Coulomb energy (exchange plus correlation energy) has a lower
bound of the form C2 n4/3 (r)dr, where n(r) is the single‐particle electron density.
208 | Density‐Functional Theory
LDA-PW
Expt.
GGA-PW91 GGA-PBE
μ (N+1)
Δxc
ϵg Eg
μ (N)
k k
Figure 4.12 Illustration of the contribution of Δxc , the discontinuity in
Vxc [n(r)]. The energies ϵ of the Kohn‐Sham equation are shown in the form
of a band structure for the (N)‐ and (N + 1)‐ electron systems. The two
differ in a uniform increase of the eigenvalues by Δxc . The quasiparticle
band‐gap Eg is the difference between the two eigenvalues as Eg = ϵN+1 −
(N+1)
ϵg = ϵN+1 − ϵN .
(N) (N)
GGA always underestimate the band gap. First, let us explain the
reason why both the LDA and GGA underestimate the band gaps
for semiconductors and insulators even though they are the basis of
good approximations for the ground‐state properties.
The origin of the band‐gap problem is the discontinuity
of chemical potential in the exchange‐correlation potential
Vxc [Sham and Schlüter (1985)]. This is because the density
is given by the summation up to the chemical potential μ as
n(r) = i Θ(μ − ϵi )|ϕi |2 , where Θ is the Heaviside step function.10
<
10 The Heaviside step function Θ(x) is a discontinuous function, whose value is zero
for negative arguments x < 0 and one for positive arguments x > 0.
>
210 | Density‐Functional Theory
Eg = ϵg + Δxc , (4.122)
where ϵg is the gap obtained from the Kohn‐Sham equation for the
ground state. In Fig. 4.12, we show the contribution of Δxc to Eg =
ϵN+1 − ϵN , where ϵJ denotes J‐th energy for (N)‐electron system.
(N+1) (N) (N)
As shown in Secs. 4.11.1 and 4.11.2, both the LDA and GGA do not
consider the discontinuity in Vxc [n(r)] since they are the smooth and
local functions of n(r). Therefore, the LDA and GGA would give zero
Δxc , which is the reason why they always underestimate the value of
Eg even they give a good approximation to ϵg .
Several hybrid functionals have been proposed to obtain non‐
zero value of Δxc . The idea of the hybrid functional is based on a
non‑local functional to the exchange‐correlation functional Exc [n(r)].
As discussed in Sec. 3.6, the exchange potential of the Hartree‐
Fock equation is a non‑local potential (see Eq. (4.47)). The exchange
energy can be obtained from the Hartree‐Fock exchange potential as
1
ϕ∗i (r1 )ϕ∗j (r2 )ϕj (r1 )ϕi (r2 )
ExHF dr1 dr2 , (4.123)
2 i,j |r1 − r2 |
=−
1 1 − erf(ηr) erf(ηr)
, (4.124)
r r r
= +
SR LR
1
EH [n(r)] = VH (r)n(r)dr, (4.126)
2
where a factor of 1/2 takes into account the double counting and
VH (r) is the Hartree potential, which is given in Eq. (4.33).
For a periodic solid, the total electron density and the potentials
can be expressed in the terms of the plane waves, eiGr , with G is the
reciprocal lattice vectors (see Sec. 5.3) as
G G
n(G) = n(r)e
−iGr
dr, and VH (G) = VH (r)e−iGr dr. (4.128)
For the neutral case, the total negative charge of the electrons is
canceled by the total positive charge of the ions. Therefore, the
average potential is zero. Since G = 0 in Eq. (4.127) corresponds
to the average potential over all the space, we can omit G = 0 in
the summation due to charge neutrality.12 Then, by substituting
12 There might be a contribution from the G = 0 term if the systems have an intrinsic
1
ei(G+G )r VH (G)n(G′ )dr
EH [n(r)] =
2 ′
′
GG
Ω
δG+G′ ,0 VH (G)n(G′ )
2
=
Ω (4.129)
VH (G)n(−G)
2
=
G̸=0
Ω
VH (G)n(G),
2
=
G̸=0
where Ω is the volume of the unit cell. Here we assume that n(+G) =
n(−G). On the other hand, by substituting Eq. (4.127) into the Poisson
equation (see Sec. 4.5), we have
n(G)
∇2 VH (r) = −4πn(r) ⇐⇒ VH (G) = 4π . (4.130)
|G|2
n2 (G)
EH = 2πΩ . (4.131)
|G|2
G̸=0
where ϵxc (r) = ϵx (r) + ϵc (r) can be expressed in the terms of the
plane waves as
(4.134)
Exc [n(r)] = Ω ϵxc (G)n(G).
G̸=0
Since the external potential term Ven (r) in Eq. (4.137) is costly
computation with the plane waves, Ven (r) is replaced by a pseudopo‑
tential Vps (r), which is related to replacing the effects of the core
electrons with an effective potential.
A local form of the pseudopotential of a single ion, vps (r),
is proposed by Ashcroft‐Heine‐Abarenkov [Ashcroft (1966),
Heine and Abarenkov (1964)] as
<
v0 r < rcut
vps (r) = , (4.138)
−r Z
r > rcut
where Z is the charge of the ionic core and rcut is the effective radius,
>
v
Core electrons Shell electrons
rcut
r
v0
vps
Z
−
r
1
vIps (|r − T − RI |)e−iGr dr
V
Vps (G) =
I T
Ncell
e vIps (|r|)e−iGr dr
V
−iGRI
=
I
(4.140)
1
e vIps (|r|)e−iGr dr
−iGRI
=
I
Ω
SI (G)FI (G),
=
I
Quantum ESPRESSO Course for Solid‐State Physics | 217
in which SI (G) and FI (G) are defined by
1
SI (G) = e−iGRI , and FI (G) = vIps (|r|)e−iGr dr,
(4.141)
Ω
factor13 and the atomic form factor14 of the I‐th atom, respectively.
For a given pseudopotential such as Eq. (4.138), the external
energy can be expressed as
(4.142)
Eext [n(r)] = Vps (r)n(r)dr = Ω Vps (G)n(G).
G
Finally, we will discuss the last term EEwald in Eq. (4.125). In Quantum
ESPRESSO, EEwald denotes the ion‐ion Coulomb energy, which is given
by
ZI ZJ
N
1n
. (4.144)
2 |RI − RJ |
EEwald =
I̸=J
13 The structure factor contains all the information about the lattice structure. The
structure factor gives the extinction rule for X‐ray diffraction (see Sec. 5.2).
14 The atomic form factor contains the information of each atom, which can be
calculated or obtained by fitting experimental data. The form factor is important for
obtaining the intensity of X‐ray scattering.
218 | Density‐Functional Theory
ZJ 1 − erf(η|R − RJ |)
Nn n N
ZI
|R − RJ | |R − RJ |
=
J J
SR
(4.145)
erf(η|R − RJ |)
Nn
ZJ ,
|R − RJ |
+
J
LR
ZJ ZJ . (4.146)
|R − RJ | V J |G|2
=
J G̸=0
Then the both SR and LR parts can converge quickly with increasing
|G| and |R| for a given value of η. Therefore, the ion‐ion Coulomb
energy can be obtained by a few terms in the summations over |G|
and J.
When we insert Eqs. (4.145) and (4.146) into Eq. (4.144), we
must subtract separately the term for R = RJ in Eq. (4.146), to avoid
the divergence. If we adopt the following equation:
erf(η|R − RJ |)
= 2√ , (4.147)
η
|R − RJ |
lim
R→RJ π
Quantum ESPRESSO Course for Solid‐State Physics | 219
Eq. (4.144) can be rewritten as
1 − erf(η|RI − RJ |)
N
1n
ZI ZJ
2 |RI − RJ |
EEwald =
I̸=J
ZI ZJ (4.148)
V I,J |G|2
+
G̸=0
Nn
Z2J √ .
η
−
J
π
The ionic forces are used to study the dynamics of ions such as
optimizing ionic positions (see Sec. 3.1.4). Using the total energy in
Sec. 4.12, we can calculate the forces that act to the ion I, which is
given by taking the derivative of the total energy with respect to
individual ionic position RI as
There are two contributions to the ionic forces, one from the
ion‐ion interaction energy EEwald in Eq. (4.148), and one from the
ion‐electron interaction energy Eext in Eq. (4.142). For the ion‐ion
interaction energy, the force of ion I is given by Eq. (4.148) as follows:
∂EEwald ZI ZJ
n N
Fion (RI − RJ ), (4.150)
I
∂RI |RI − RJ |3
=− =
J̸=I
∂Eext ∂Vps
Fion‐e ϕ . (4.151)
I
∂RI ∂RI i
=− =− ϕi
i
∂Vps (r)
Fion‐e n(r)dr.
(4.152)
I
∂RI
=−
15 The Hellmann‐Feynman theorem states that the forces are given by the
1 d ℓ(ℓ + 1)
r2 R(r) = ϵR(r), (4.154)
d
2r2 dr 2r2
− + + V(r)
dr
1 d2
u(r) = ϵu(r). (4.156)
2 dr2
− + V(r)
d2 u
(4.157)
dr2
= −k(r)u(r),
c0 = 1 − 12 h kn (r)
5 2 2
c1 = 1 + 12 h kn−1 (r)
1 2 2 , (4.160)
c = 1 +
12 h kn+1 (r)
1 2 2
2
u2 (r)dr = 1.
(4.161)
Then, the electron density n(r) of a single electron for the 1s state
(Y(θ, φ) = 1/ 4π) is given by
√
R2 (r) u2 (r)
n(r) = . (4.162)
4π 4πr2
=
Note that, for the helium atom, the total electron density is 2n(r) since
we have two electrons.
How to run: To run this tutorial, the readers will do the following
command lines:
equations of second order in which the first‐order term does not appear.
Quantum ESPRESSO Course for Solid‐State Physics | 223
Input file: The input files include a subroutine file
(schrodinger.py) and a main file (check-schrodinger.ipynb).
The subroutine contains the Numerov algorithm in Eq. (4.159) and
the radial Schrödinger equation in Eq. (4.157), and the main file
calculates the ground‐states energy and wavefunction for the helium
atom.
SSP‑QE/dft‑he/schrodinger.py
38 else:
39 assert num_nodes > 0, 'expect #nodes >0
since u[0]<0 while u[infty]>0'
40 eps_max = eps
41 raise Exception('Not converged after %d
iterations.' %( maxiter))
SSP‑QE/dft‑he/check‑schrodinger.ipynb
Figure 4.15 Electron density of helium. Solid and dashed lines denote the
electron densities, which are obtained by the Numerov algorithm and
analytical methods, respectively.
that the electron density from the numerical method (the Numerov
algorithm) reproduces the exact electron density of the helium atom
given by Eq. (4.21).
u2 (r)
∇2 UH (r) = − . (4.164)
r
17 The Verlet algorithm is a simple method for integrating second order differential
SSP‑QE/dft‑he/poisson.py
SSP‑QE/dft‑he/check‑poisson.ipynb
Figure 4.16 The Hartree potential of helium. Solid and dashed lines denote
the Hartree potential, which are obtained by the Verlet algorithm and
analytical methods, respectively.
4
ϵx (r) = C2 n1/3 (r) and Vx [n(r)] = C2 n1/3 (r), (4.166)
3
>
Quantum ESPRESSO Course for Solid‐State Physics | 229
where the parameters A, B, C, D, γ, β1 , β2 are given in Table 4.2, too.
Then, the correlation potential is given by
∂[ϵc (r)n(r)]
. (4.168)
∂n(r)
Vc [n(r)] =
rs d
Vc (rs ) = 1− ϵc (rs ). (4.169)
3 drs
(4.170)
>
1 d2
en H xc u(r) = ϵu(r), (4.171)
2 dr2
− + V (r) + V (r) + V [n(r)]
where Vxc [n(r)] = Vx [n(r)] + Vc (rs ). Then, the total energy of the
helium atom in Eq. (4.103) can be rewritten by
E = 2ϵ −
VH (r)u2 (r)dr − ΔExc [n(r)], (4.172)
where ϵxc [n(r)] = ϵx [n(r)] + ϵc (rs ). Note that we have to use the
electron density n(r) and VH (r) for the helium atom.
How to run: To run this tutorial, the readers will do the following
command lines:
SSP‑QE/dft‑he/dft‑lda‑he.ipynb
3
rj = sij ai . (5.1)
>
i
It should be careful to set sij in the unit cell when the unit cell is not
cubic or cuboid. A lattice vector which is defined by
where R represents a position of the unit cell in the crystal. Thus, the
coordinate of any atom in a crystal is given by the sum of R + rj .
Reciprocal lattice vector: Any wave in the crystal2 that
propagates in the direction of k is expressed by f(k, r) ≡ exp(ik · r)
which we call a plane wave. The product of k · r in exp(ik · r) is called
a phase of the wave that has a periodicity of 2π. If we take r = ai and
k = bi so that ai · bi = 2π, we get the same amplitude of f for r and
r + ai , that is f(bi , r) = f(bi , r + ai ). Further, if we take r = ai , we get
f(k, ai ) = f(k + bi , ai ). It means that the wave has a periodicity not
only in the real space but also in the k space, which we call reciprocal
lattice.3 The vector bi is a reciprocal lattice vector. The unit box in
the k space which consists of the three bi is called the Brillouin
zone. Similar to the lattice vector R, we can define reciprocal lattice
vector, G,
converted to each other. In physics, we call the extended 6 dimensional space of {r, k}
as a phase space. Any motion of a particle is expressed by a curve in the phase space.
Quantum ESPRESSO Course for Solid‐State Physics | 237
where G represents the position of the Brillouin zone in the k space.
The integers of (ℓ, n, m) are called the Miller indices which are used
in X‐ray analysis (Sec. 5.2).
The bi is obtained in the output file of Quantum ESPRESSO since
bi can be calculated by the ai in the input file as follows:
1 if i = j
ai · bj = 2πδij , . (5.4)
0 if i ̸= j
δij ≡
In fact, we have the following formula for bi which satisfies Eq. (5.4):
(a) λ (b)
kin kout
θ θ Δk
θ θ d
Figure 5.1 (a) Bragg’s condition in the real space. (b) Bragg’s condition is the
k space.
have one peak at a 2θ, which we can label by a Miller index (ℓ, n, m)
(see Sec. 5.1), which we explain below.
Bragg’s condition: For a periodic lattice, we can define a crystal
plane as is shown by thin lines in Fig. 5.1 (a) on which many atoms
exist. Since the lattice is periodic, we have many equivalent crystal
planes parallel to the crystal plane. In Fig. 5.1 (a), we show the side‐
view of the crystal planes as parallel lines. When the incident X‐ray
with the wavelength λ enters the crystal with an angle θ (solid dots
in Fig. 5.1 (a)) measured from the plane, the diffracted X‐ray that
propagates in the direction of θ, too, interfere to each other from
two adjacent planes as shown in Fig. 5.1 (a). The condition of the
constructive interference of the diffracted light is given by
nλ
2d sin θ = nλ, , (5.6)
2d
sin θ =
6 If you know diffraction grating, the concept of the diffracting grating comes from
nλ 2π
Δk = 2k sin θ = 2k · · n, (5.7)
2d d
=
2π 2π
d= . (5.8)
(ℓb1 + mb2 + nb3 )2
=
|G|
CG (k)ei(k+G)r , (5.9)
Ψk (r) =
G
8 The Schrödinger equation is given by the differential equation, and thus we solve
h̄2
(5.10)
2m
H(r) = − Δ + V(r),
where − 2m
h̄
Δ and V(r) are kinetic and potential energy operators,
2
1
VG e , where VG =
iGr
dr′ V(r′ )e−iGr . (5.11)
Nu
′
V(r) =
G
h̄2
(k + G)2 CG ei(k+G)r + VG′ CG′′ ei(k+G +G )r
2m
′ ′′
G G′ ,G′′ (5.12)
=E CG ei(k+G)r .
f(x) = fp ei2πpx/L
∞
p=0
which is called the Fourier series expansion. Here 2π/L is the reciprocal lattice
vector for L.
10 The Laplacian Δ are invariant for transformation of x′ = x + a and
V(r + a) = V(r).
242 | Solid‐State Physics for Quantum ESPRESSO
In the second term, we put G′ + G′′ = G and take the term of ei(k+G)r ,
we get the following equation,11
h̄2
(k + G)2 CG + VG−G′′ CG′′ = ECG . (5.13)
2m ′′ G
h̄2
HGG′ = 2m (k + G) , if (G = G ) .
2
(5.14)
′
V ′,
G−G if (G ̸= G′ )
h̄2 2
Ecut‐off = G , (5.15)
2m max
where Gmax is the largest reciprocal lattice vector. We use the cut‐
off energy because the values of N depend on the size of the unit
11 Since ei(k+G)r are independent functions for different values of G, the coefficients
of the both‐hand sides should be identical. Further, we change the summation on G′
and G′′ to the sum on G and G′′ . Note that G′ = G − G′′ .
12 LAPACK is the name of a famous software package for linear algebra calculation
Figure 5.2 Cut‐off energy. Dots in the k space represent the reciprocal lattice
vectors G. When we define a sphere whose diameter is |Gmax |, we will use
G’s in the sphere. Then cut‐off energy is given by Eq. (5.15). (a) When the
unit cell is large, G becomes small, and thus we have more G’s in the sphere.
(b) When the unit cell is relatively small, the number of G’s becomes small
even for the same cutting energy. Reprinted with permission from Asakura
Publishing Co. Ltd.
cell. When the size of the unit cell is relatively large because of many
atoms in the unit cell, the corresponding absolute value of G becomes
small as shown in Fig. 5.2 (a). In order to obtain the same accuracy for
calculating the atomic potentials in the smaller unit cell (Fig. 5.2 (b)),
a larger Gmax is needed. On the other hand, when we use the same cut‐
off energy for two different sizes of unit cells, we expect the accuracy
of the calculation to be similar to each other.
The selection of cut‐off energy also depends on the selection of
pseudopotential for each atom, which we explain below. As shown
in Eq. (5.12), the off‐diagonal matrix element of the Hamiltonian
matrix is the Fourier transform of the crystal potential V, which
is given by the sum of atomic potentials. If we adopt the real
atomic potential for an atom, we expect a large Coulomb potential,
−Ze2 /(4πϵ0 r) (Z is the atomic number of the atom). Fourier
transform of 1/r is 1/|G|2 which is a slowly decreasing function of
|G| compared with exp(−C|G|) (C is a constant) in three dimensions.
Thus, it is not a good idea that we do not adopt the real atomic
potential in the plane wave expansion.
Pseudopotential: For discussing most of the solid‐state proper‐
ties, inner‐core electrons in the 1s, 2s, 2p atomic orbitals for a heavy
element are not so important. Only the valence electrons near the
highest occupied orbitals contribute to the properties. The valence
244 | Solid‐State Physics for Quantum ESPRESSO
easily excite the unoccupied states within the same energy band. This
situation is called “metal”. We can say that metal is defined by a finite
density of states at the Fermi energy, D(EF ) ̸= 0. When the valence
band is fully occupied, and the next higher energy band, that is
conduction band is unoccupied, we usually have an energy gap. An
energy gap, Eg , is defined by an energy difference between the highest
energy of valence band and the lowest energy of conduction bands.13
Since there are no states in the energy gap region, the occupied
electron in the valence band needs at least Eg for exciting the electron.
This situation corresponds to “semiconductor” or “insulator”. The <
difference between “semiconductor” or “insulator” is the value of Eg .
If Eg <3eV, we can say that the material is semiconductor, while Eg > 5
eV, we call “insulator”.14
When Eg = 0 or Eg is a small value compared with the thermal
energy kB T (kB is the Boltzmann constant), part of electrons in
>
13 If you are a chemist, you can imagine that an energy gap is a kind of HOMO‐LUMO
A hole is an unoccupied state in the almost occupied energy band. Like a going‐up
beer bubble,
> the > unoccupied states carry a positive charge in the solid in the
direction opposite to that of the electron.
Quantum ESPRESSO Course for Solid‐State Physics | 247
(a) D(E) | M |2 (b)
E
A
I
0
EF + eV −I tip −V
EF
sample tip sample
−V 0
Figure 5.3 Scanning tunneling spectroscopy. (a) Tunneling current I (Eq.
(5.16)) flows from the metallic materials (left) to a scanning metallic tip
(right) in the energy region between EF and EF + eV, in which electrons
occupied in the metallic side while not occupied in the tip side. (b) In the
experiment, we fix the position of the tip, we change the voltage. Then the
value dI/dV is proportional to the density of state D(EF + eV) (Eq. (5.17)).
surface in the k space on which the energy is the Fermi surface) which
costs more CPU times. Total energy depends on the number of the
sampling points, which should be checked by increasing the sampling
points.
800
400
0 -2
0.0 1.0 x 10
Γ M K Γ -1
State/1C-atom/cm
Figure 5.4 (a) Phonon dispersion of graphene. The unit of frequency is cm−1
(1 eV = 8065 cm−1 ). (b) Phonon density of states.
transverse optic (TO) modes. The assignment of LA, TA, LO, and TO is
possible by looking at the direction of oscillation of each atom at the
Γ point (or near the Γ point for acoustic modes by the Material Cloud
(see Sec. 3.3.1.)).
When the material is two‐dimensional, such as graphene as
shown in Fig. 5.4, TA modes are further classified by in‐plane TA
(iTA) and out‐of‐plane TA (oTA), and TO modes are classified by in‐
plane TO (iTO) and out‐of‐plane TO (oTO) phonon modes. In LA and
LO modes, since the propagation direction is in‐plane, we do not
usually say iLA or iLO. On the other hand, we sometimes call oTA and
oTO modes as ZA and ZO modes which shows that the oscillational
direction is in the direction of z axis for xy plane of two‐dimensional
materials.
Phonon dispersion of a given material is calculated by solving
the so‐called dynamical matrix. In the dynamical matrix, we define
force constants between two atoms in the crystal. In order to
define the force constant, any atom should have a restoring force
for all directions of the displacement of the atom. This situation
corresponds to the ground state in which the total energy has a
minimum at the optimized position as a function of displacement of
each atom. The restoring force can be calculated by the gradient of the
total energy by displacing an atom in the unit cell with keeping other
250 | Solid‐State Physics for Quantum ESPRESSO
V(r)
rmin
r
in which the second term of the right‐hand side of Eq. (5.18) is called
deformation potential. When we adopt the adiabatic approximation
δR
t t t t t
v(r−R)
v(r﹣R)
Figure 5.6 Deformation potential. Solid lines are atomic potentials. When an
atom shifts position by δR, the potential energy for an electron (showing
solid circles) of the nearest neighbor atom increases (or decreases) by
moving away (or approaching) the atom. Further the absolute value of the
transfer integral |t| decreases or increases by moving away (or approaching)
the atom.
(a) (b)
R k′
R′′
R1
R2 k
q
R′
Figure 5.7 (a) Three‐centers integral M(R′ , R, R′′ ) in Eq. (5.20) is a function of
the positions of three atoms, R′ , R, R′′ . The integration is given as a function
of the relative coordinate R1 = R − R′ , R2 = R′′ − R′ . (b) The momentum of
an electron h̄k is scattered to h̄k′ = h̄k′ − h̄q by emitting a phonon with the
crystal momentum h̄q.
that the atomic wavefunction follows the position of the atom, the
atomic wavefunction of the distorted atom is given by
1
Vk′ ,k = exp(−ik′ R1 + ikR2 + i(−q − k′ + k)R′ )M(R1 , R2 )
N
R1 ,R ,R2
= δ(−q − k′ + k) exp(−ik′ R1 + ikR2 )M(R1 , R2 )
′
R1 ,R2
(5.21)
Figure 5.8 When light comes to a material, the light is transparent (left),
reflected (center), or scattered (right). Even for the case of transparent, the
direction of light changes after the transmission. Reprinted with permission
from Asakura Publishing Co. Ltd.
When light hits a material (see Fig. 5.8), the light is transparent (left),
reflected (center), or scattered (right). Even in the case of transparent
properties, the direction of light changes after the transmission. It
means that all phenomena are related to the interaction of a photon
and materials.
Light or a photon interacts with an electron (or electrons) in
the materials. The optical properties of the materials consist of
(1) absorption, (2) emission, and (3) scattering of light. In the
case of semiconductors, optical absorption (or emission) occurs
for the photon whose energy is more than the energy gap of the
semiconductor. Suppose the momentum and energy of a photon are
matched to the difference of momentum and energy between the
initial and final electronic states. In that case, the optical absorption
occurs by annihilating the photon and by exciting the electron from
the initial to the final states. This type of excitation for the electron
is called single‐particle excitation. On the other hand, for metal
or heavily doped semiconductors, an electric field of a photon as
Quantum ESPRESSO Course for Solid‐State Physics | 255
an electro‐magnetic wave can excite an electric current of metal.
The optical absorption occurs by the Joule heat by the current. In
this case, since the current is the collective motion of electrons,
this type of excitation is called collective excitation. Usually, the
single‐particle excitation and collective excitation occur exclusively
or simultaneously for a given energy of the photon.
Scattering of light in the semiconductor is defined as a sequential
process of optical absorption and emission, in which the emitted light
propagates in any direction of solid angle. The scattering amplitude
depends on materials as a function of photon frequency, which
determines the color of the materials. On the other hand, in metal, the
reflection of light occurs by screening the electric field of the photon
by the induced electric current, which is the origin of metallic luster.
For the region of the wavelength of the visible light, from 400 to
800 nm, there are two kinds of scattering of light: elastic and inelastic
scattering of light, which we call the Rayleigh and Raman scattering,
respectively. In the Rayleigh (Raman) scattering, the scattered light
does not (does) change the energy from that of the incident light.
The scattering amplitude becomes strong if the energy difference of
the initial and the final electronic states is the same as the photon
energy for either the incident or scattered photon. This effect is
called resonant Rayleigh (Raman) scattering. The resonant Rayleigh
scattering can be one of the reasons the materials have an intrinsic
color.17
Optical absorption can be explained by the perturbation theory
for electronic structure in which the perturbation Hamiltonian is
given by the vector potential of the electromagnetic field as follows:
1
{−ih̄∇ − eA(t)} + V(r)
2
2m
H=
1 2
−h̄ Δ + ieh̄A(t) · ∇ + ieh̄∇ · A(t) + e2 A(t)2 + V(r).
2m
=
(5.22)
17 Other reasons are (1) absorption spectra by induced current, (2) structural color,
and so on.
256 | Solid‐State Physics for Quantum ESPRESSO
ieh̄
A(t) · ∇. (5.23)
m
H′ =
dP
= 2 | ⟨f|H′ |i⟩ |2 δ(Ef − Ei − h̄ω), (5.24)
π
dt h̄
where δ is the Dirac delta function. The delta function tells us that
the optical transition occurs when the energy difference between
the initial state Ei and the final state Ef is matched to the energy
of the phonon h̄ω, which we call resonant optical transition. In the
mathematics, the delta function δ(x) describes a singular function at
x = 0. However, in physics, because of the uncertainty principles for
the energy, ΔEΔt ∼ h̄, the strict condition of x = 0 can be relaxed for
the fast optical transition with a small Δt. For example, if Δt = 1 ps
(or 1 THz), ΔE becomes about 1 meV. ΔE = 1 meV means that the
optical absorption spectra have a spectral width of 1 meV. For the
case of Raman scattering, a typical value of scattering is Δt = 10 fs,
ΔE becomes 0.1 eV.20
⟨f|H′ |i⟩ in Eq. (5.24), which is called transition dipole moment,
is calculated by putting the wavefunction of valence and conduction
18 The gauge is an additional condition for a vector potential combined with a scalar
electro‐static potential that does not change the value of electric and magnetic fields.
If we adopt the Coulomb gauge, the scalar electro‐static potential satisfies the
Poisson equation, while if we adopt the Lorentz gauge, the scalar electro‐static
potential satisfies the wave equation.
19 The Fermi golden rule is a general formula of the rate of transition from the initial
to the final states for a given time‐dependent perturbation that has a time
dependence of exp −iωt. See detail in any textbook on quantum mechanics.
20 ΔE = 0.1 eV means that the resonance Raman scattering occurs within 0.1 eV
energy width as a function of incident laser energy, which we call resonance Raman
profile. It is noted that the spectral width of Raman spectra comes from the lifetime
of a phonon.
Quantum ESPRESSO Course for Solid‐State Physics | 257
bands for the initial and final states, respectively. Here we can express
the vector potential of the light by polarization vector P that is the
unit vector in the direction of A (or E),
i I0
A= exp {i(kopt · r ± ωt)}P, (5.25)
ω cε0
eh̄ I
exp {i(ωf − ωi ± ω)t}Dfi · P. (5.27)
mω cε0
′
⟨f|H |i⟩ =
G G′
=i Cf∗
G′ (k)CG (k)(k + G)δ G,G
i
′
G G′
=i Cf∗
G (k)CG (k)(k + G).
i
(5.28)
|ε(ω)| + Re(ε(ω))
n(ω) = , (5.29)
2
258 | Solid‐State Physics for Quantum ESPRESSO
and
|ε(ω)| − Re(ε(ω))
. (5.30)
2
κ(ω) =
and
2ωκ(ω)
, (5.32)
c
α(ω) =
V
JDOS(E) = dkδ {Ecσ (k) − Evσ (k) − E},
(5.34)
σ
(2π)3
L
R=ρ , (5.35)
W2
(Ω)
W2
G=σ , (5.36)
L
(S)
nc e2 τE
J = nc ev̄ = , (5.37)
m
n c e2 τ
. (5.38)
m
σ=
1
τ n (k) = , (5.41)
Γn (k)
2π dq
|Mmn (k + q, k)|2 ×
h̄ mμj
Γn (k) =
ΩBZ
(5.42)
1+j
− jfm,k+q + nq.μ δ (ΔEk.q − jh̄ωqμ ),
0
2
(a) (b)
ℏω′ ℏω
q′ q
ℏω ℏω′′
q q′′
ℏω′′ ℏω′
q′′ q′
1
nqj ≡ , (5.43)
e h̄ωj (q)/kB T −1
1
γ j (q, T) ≡
τ j (q, T)
2
Vqj,q′ j′ ,q′′ j′′ δ(−q + q′ + q′′ + G) (5.44)
π
h̄ Nq ′ ′′ ′
(3)
= 2
q ,q ,j ,j′′
×[(1 + nq′ j′ + nq′′ j′′ )δ(ωqj − ωq′ j′ − ωq′′ j′′ )
+2(nq′ j′ − nq′′ j′′ )δ(ωqj + ωq′ j′ − ωq′′ j′′ )] ,
264 | Solid‐State Physics for Quantum ESPRESSO
When there are hot and cold points in material, heat flows from the
hot to the cold points. The heat flux J in units of W/m2 at a point is
proportional to the gradient of temperature ∇T, (K/m), that is,
which is called the Fourier law, and the coefficient κ > 0 is called
thermal conductivity (W/m/K). In Quantum ESPRESSO, we can
calculate the thermal conductivity by Thermal2 code (https://
anharmonic.github.io/thermal2/), which comes bundled with
the D3Q code. Both Thermal2 and D3Q codes can usually install
directly inside Quantum ESPRESSO by doing "make d3q". The
thermal conductivity κ becomes large when the amplitude of the
oscillation is large compared with a few % of the bond length.
Diamond and carbon nanotubes are known to be materials with a
large thermal conductivity since (1) the mass of carbon atoms is
relatively small for a given spring constant and (2) the anharmonicity
of the potential is large compared with other materials. It is noted
here that the thermal conductivity of metal is generally large, too,
because phonons are interacted by electron‐phonons interactions.
Thus, the phonon scattering process is relevant to the thermal
conductivity.
266 | Solid‐State Physics for Quantum ESPRESSO
dQ
(5.47)
dt
= −divJ,
Q = ρCT. (5.48)
dQ dT
(5.49)
dt dt
= ρC = −divJ = κΔT.
h̄2 2
ω τ j (q, T)nqj (nqj + 1)|vqj |2 . (5.50)
2VkB T2 qj qj
κ=
h̄2 V
Cvqj = ω2 nqj (nqj + 1), C= Cvqj . (5.51)
VkB T2 qj 2N qj
qj
scattering.
268 | Solid‐State Physics for Quantum ESPRESSO
The energy difference between the incident light and the scattered
light is called the Raman shift, whose unit is cm−1 (1 eV = 8065
cm−1 ). In the experiment, they observe the intensity of the scattered
light as a function of the Raman shift, which is called Raman spectra.
The peak positions of Raman spectra correspond to the phonon
modes, which are Raman active modes. The Raman active modes
are defined by the optical phonon modes in which the photo‐excited
electron emits a phonon by the electron‐phonon interaction. The
Raman active modes have the symmetry of quadratic form of x,
y, z such as x2 + y2 , x2 − y2 , xz, etc. The quadratic function for
each phonon mode is obtained as functions for the irreducible
representation in the character table of point group for the unit cell
if we know the information of the irreducible representation of the
phonon mode in the point group.23
In Quantum ESPRESSO, we can calculate the so‐called non‐
resonant Raman spectra to obtain the phonon mode’s frequency
and symmetry and the relative intensity in the Raman spectra,
which we call density‐functional perturbation theory (DFPT).
The non‐resonant Raman intensities are calculated within the
DFPT [Lazzeri and Mauri (2003)] as
nν + 1
I(ν) ∝ |⃗es · R(ν) · ⃗ei | , (5.53)
2
ων
23 Number of the Raman active modes and their irreducible representation are
a 0 0 0b0
∂3 Eel uν
Rαβ (ν) = √ ζδ , (5.55)
∂Eα ∂Eβ ∂rζδ Mδ
ζδ
where Mδ the mass of the δth atom, rζδ the position of the δth atom
along the direction ζ (ζ = x, y, z), and uζδ the atomic displacement of
the δth atom along the ζ direction for the eigenmode ν.
For the calculated Raman tensor and polarization vectors, we
can calculate the relative Raman intensities by |⃗es · R(ν) · ⃗ei | in
2
Eq. (5.53). For example, when the polarization of the incident light
lies in the x direction, ⃗ei =t (1, 0, 0), we get a non‐zero intensity for
the x or y polarizations of the scattered light ⃗es = ex ≡t (1, 0, 0), or
⃗es = ey ≡t (0, 1, 0) for R(E2g : x2 − y2 ) or R(E2g : xy), respectively, as
follows:
a 0 0 1
t
ex R(E2g : x2 − y2 )ex = (1, 0, 0) 0 −a 0 0 = a, (5.56)
0 0 0 0
or
0b0 1
t
ey R(E2g : xy)ex = (0, 1, 0) b 0 0 0 = b. (5.57)
000 0
is called zone‐center phonon mode. The reason why only the Γ point
phonon is observed in the first‐order Raman spectra is understood
by the three Raman sub‐processes: (1) optical absorption of an
electron with the wavevector k in which the electron can excite
almost “vertically” in the k space, (2) the photo‐excited electron can
emit a phonon with a phonon wavevector q, and (3) the scattered
electron is recombined with a hole at the original k. In order for the
scattered electron to recombine with a hole, the phonon wavevector
should be q = 0.
If we can not find the corresponding phonon frequency in
the observed Raman spectra, the Raman spectra might be two‐
phonon Raman spectra whose frequency is the sum of two‐phonons.
In the two‐phonon Raman scattering, the restriction of q = 0 in
the case of one‐phonon scattering is relaxed, and a pair of q ̸= 0
and −q can be possible for the wavevectors of the emitted two
phonons, which enables the photo‐excited electron to recombine
with the hole. Thus we expect that the two phonon spectra are
generally broad and weak. When two of the three intermediate states
are resonant to the electronic states, the Raman intensity of the
two‐phonon Raman scattering becomes comparable or even larger
than the intensity of one‐phonon Raman scattering. This situation
is called “double‐resonance Raman scattering” [Saito et al. (2002b),
Saito et al. (2003), Saito et al. (2002a)]. The double resonance Ra‐
man occurs, too, when one of the two scattering processes with
q and −q is an elastic scattering of a photo‐excited electron by
impurity. One famous example of defect‐oriented, double reso‐
nance Raman spectra is the D‐band of graphite at 1350 cm−1
[Pimenta et al. (2007)] in which iTO phonon mode at the K point is
relevant to the Raman spectra. The two‐phonon Raman spectra of the
iTO phonon mode at the K point appear at 2700 cm−1 which we call
G′ band (or 2D band). It is noted that the G′ band is an intrinsic Raman
spectra of graphite without any defects.
When the intermediate state of the photo‐excited electrons
is a real electronic state, the Raman intensity is enhanced sig‐
nificantly compared with non‐resonant Raman spectra, which we
call resonant Raman spectra. In order to calculate the resonant
Raman intensity by first‐principles calculation, we need to calculate
electron‐phonon and electron‐photon matrix elements. See detail in
the references [Tatsumi and Saito (2018), Tatsumi et al. (2018)].
Quantum ESPRESSO Course for Solid‐State Physics | 271
5.14 Warnier functions
1 −ik·R
wnR (r) = e (5.59)
Nk
ψnk (r),
k
and
nR
272 | Solid‐State Physics for Quantum ESPRESSO
k2
0 R1 R2
Unit cell
Figure 5.10 (a) The Bloch functions ψnk (r) of a single band n in 1D system for
three different values of the wave vector k. (b) The Wannier functions wnR (r)
for the same band n for the three different values of the lattice vector R.
where unk (r) is a periodic function (unk (r) = unk (r − R)). By substi‐
tuting Eq. (5.62) into Eq. (5.59), we obtain:
1 ik(r−R)
wnR (r) = e unk (r)
Nk
k
1 (5.63)
ψnk (r − R)
Nk
=
k
= wn,0 (r − R).
Therefore, wnR (r) has the same periodicity as the atomic structure of
the crystal, as shown in Fig. 5.10 (b).
It is important to note that the WFs are not unique because the
Bloch functions are not unique. For each ψnk (r), we can obtain the
same electron density by multiplying a different phase α in front of
the Bloch function (i.e., eiα ψnk (r)). However, by choosing a different
phase of the Bloch functions, we get new Wannier functions that
may exhibit very different shapes in real space. When we consider
Ne energy bands, Ne Bloch functions at a given k can be rotated by
Quantum ESPRESSO Course for Solid‐State Physics | 273
the unitary operator as
Ne
Ukmn ψmk (r), (5.64)
ψnk (r) =
m=1
Ne Ne
⟨(r − r̄n )2 ⟩n = ⟨r ⟩n − r̄2n , (5.66)
2
Ω=
n=1 n=1
where r̄n and ⟨r2 ⟩n are the center and the second moment of the
average for WF at R = 0 that are defined by
For the 1D system, the spread of the WFs is shown in Fig. 5.10 (b).
Since wn,0 is a function of Ukmn (see Eq. (5.65)), Ω is a functional
of Ukmn . Thus, it allows us to minimize the spread Ω by variation
of these unitary matrices. To perform this minimization, steepest‐
descent or conjugate‐gradient algorithms can be applied if dΩ/dUkmn
is known [Marzari et al. (2012b)].
274 | Solid‐State Physics for Quantum ESPRESSO
1 ∂
r̄n = ⟨wn,0 |r|wn,0 ⟩ = unk i unk . (5.68)
Nk ∂k
k
1 ∂ 2
⟨r2 ⟩n = ⟨wn,0 |r2 |wn,0 ⟩ = − unk u . (5.69)
Nk ∂k nk
k
i
r̄n = wb b[⟨unk |un,k+b ⟩ − 1], (5.70)
Nk
k,b
and
1
⟨r2 ⟩n = wb [2 − 2Re⟨unk |un,k+b ⟩], (5.71)
Nk
k,b
Mk,b
mn = ⟨umk |un,k+b ⟩. (5.72)
nn − 1 = iIm ln Mnn ,
Mk,b k,b
(5.73)
and
nn = 1 − Mnn
+ Im ln Mk,b
k,b 2 2
2 − 2ReMk,b nn . (5.74)
Quantum ESPRESSO Course for Solid‐State Physics | 275
By substituting Eqs. (5.73) and (5.74) into Eqs. (5.70) and (5.71), we
obtain:
1
r̄n = − wb bIm ln Mk,b
nn , (5.75)
Nk
k,b
and
1
wb 1 − Mk,b M
2 k,b 2
. (5.76)
⟨r2 ⟩n =
Nk nn nn
+ Im ln
k,b
From Eqs. (5.66), (5.75), and (5.76), we can see that the spread
Ω with respect to the unitary matrix Ukmn can be expressed as a
mn . Therefore, Mmn is needed
function of the overlap matrix element Mk,b k,b
1 k
|w̃n,0 ⟩ = Amn |ψmk ⟩, (5.77)
Nk
mk
where tnn′ (R, R′ ) are the hopping parameters that depend on only the
distance vector connecting the MLWFs as
Productivity Tools
Figure 6.2 AFLOW search page with some highlighted steps for obtaining the
GeTe information. (1) Push Ge, (2) push Te and then we get Ge, Te in (3) and
we select “2” for the number of species.
Figure 6.3 Sample search results for GeTe. We click the first entry to get the
information of bulk GeTe in its cubic phase with the Fm3̄m (#225) space
group.
specified “Ge, Te” keyword. The different entries mean that various
structures and parameters for the GeTe calculations are reported to
be explored.
For example, in the “space group” column, we can find Fm3̄m
(#225) and R3m (#160) groups that correspond to the cubic and
rhombohedral lattice systems, respectively, as shown in the first and
the second rows of Fig. 6.3. The same space groups may appear
several times, indicating the other details of the parameters used in
the structures. Here we take the first row of the search results as
our example and click the corresponding GeTe entry as highlighted
in Fig. 6.3. This GeTe entry gives information about the bulk GeTe in
its cubic phase.
By opening the page after clicking the first entry of GeTe, we
will see several structural details and physical properties of the
bulk cubic GeTe. The first information available on that page is the
structure visualization. The screenshot is given in Fig. 6.4, from
which we can obtain the crystallographic information file (CIF) of the
material. It is recommended to set the relaxed structure as calculated
in AFLOW by clicking the corresponding button indicated by number
1 in Fig. 6.4. Then, we can save it in the CIF format by clicking the
button number 2 in Fig. 6.4. We may save it as any filename we like.
In this example, we use the default one given by AFLOW:
Quantum ESPRESSO Course for Solid‐State Physics | 281
Ge1Te1_ICSD_638014.cif
Note that AFLOW, in this case, already gives us the relaxed geometry
of GeTe. If we are not satisfied with the relaxed structure, we may
repeat the structural optimization process in Quantum ESPRESSO by
following the related hands‐on in Sec. 3.1.5. However, for simplicity,
in this chapter, we will not do such a process again. Instead, we rely
on the relaxed structure given by AFLOW.
On the same AFLOW page of Fig. 6.4 below the structure visual‐
ization, we can explore physical properties such as thermodynamic,
magnetic, and electronic properties calculated by AFLOW. When
we further scroll down the page, there are also AFLOW calculation
282 | Productivity Tools
Figure 6.5 The default view of Quantum ESPRESSO input generator and
structure visualizer in MaterialsCloud (materialscloud.org).
Figure 6.6 Options selected to generate the SCF input file for the bulk cubic
GeTe in MaterialsCloud.
$ unzip PWscf.zip
After unzipping the file, we will see that there is a folder named
PWscf containing a PWscf.in file and a pseudo directory with the
suggested pseudopotential files from MaterialsCloud. To match with
the material’s name (GeTe), it is better if we rename the folder and
the SCF input file by the following set of command lines:
$ mv PWscf gete
$ cd gete
$ mv PWscf.in gete.scf.in
Quantum ESPRESSO Course for Solid‐State Physics | 285
The first line means that we change the name of the PWscf folder
to gete. The second line is to enter the gete folder (which was
previously named PWscf). In this folder, by using the third command
line, we rename one more file, i.e., PWscf.in to gete.scf.in. Now
we can see the contents of gete.scf.in by either opening it with a
text editor or by typing the cat command in the terminal:
$ cat gete.scf.in
QE‑SSP/gete/gete.scf.in
1 & CONTROL
2 calculation = 'scf '
3 outdir = './out/'
4 prefix = 'gete ' # originally 'aiida '
5 pseudo_dir = './ pseudo/'
6 verbosity = 'high '
7 /
8 &SYSTEM
9 ecutrho = 3.2000000000d+02
10 ecutwfc = 4.0000000000d+01
11 ibrav = 0
12 nat = 2
13 ntyp = 2
14 /
15 & ELECTRONS
16 conv_thr = 4.0000000000d-10
17 mixing_beta = 4.0000000000d-01
18 /
19 ATOMIC_SPECIES
20 Ge 72.64 ge_pbesol_v1 .4. uspp.F.UPF
21 Te 127.6 te_pbesol_v1.uspp.F.UPF
22 ATOMIC_POSITIONS crystal
23 Ge 0.0000000000 0.0000000000 0.0000000000
24 Te 0.5000000000 0.5000000000 0.5000000000
25 CELL_PARAMETERS angstrom
26 4.2603675143 0.0000000000 0.0000000000
27 2.1301837572 3.6895864968 0.0000000000
28 2.1301837572 1.2298621656 3.4785755089
29 K_POINTS automatic
30 13 13 13 1 1 1 # originally 10 10 10 0 0 0
Next, let us prepare the input files for DOS and band structure
calculations from the information already available in the previous
subsection combined with the new information we will gather here.
In the following example, we use gete.scf.in as the starting
point for the DOS and band structure calculations, which are both
categorized as the NSCF calculations. We can copy the SCF input file
to make new input files for the NSCF calculations of DOS and band
structure:
$ cp gete.scf.in gete.nscfdos.in
$ cp gete.scf.in gete.nscfbands.in
QE‑SSP/gete/gete.nscfdos.in
1 & CONTROL
2 calculation = 'nscf '
Quantum ESPRESSO Course for Solid‐State Physics | 287
3 outdir = './out/'
4 prefix = 'gete '
5 pseudo_dir = './ pseudo/'
6 verbosity = 'high '
7 /
8 &SYSTEM
9 ecutrho = 3.2000000000d+02
10 ecutwfc = 4.0000000000d+01
11 ibrav = 0
12 nat = 2
13 ntyp = 2
14 nbnd = 30 # same as in *. nscfbands.in
15 occupations = 'tetrahedra_opt ' # only for DOS
16 /
17 & ELECTRONS
18 conv_thr = 4.0000000000d-10
19 mixing_beta = 4.0000000000d-01
20 /
21 ATOMIC_SPECIES
22 Ge 72.64 ge_pbesol_v1 .4. uspp.F.UPF
23 Te 127.6 te_pbesol_v1.uspp.F.UPF
24 ATOMIC_POSITIONS crystal
25 Ge 0.0000000000 0.0000000000 0.0000000000
26 Te 0.5000000000 0.5000000000 0.5000000000
27 K_POINTS automatic
28 26 26 26 1 1 1 # twice the SCF 's case
29 CELL_PARAMETERS angstrom
30 4.2603675143 0.0000000000 0.0000000000
31 2.1301837572 3.6895864968 0.0000000000
32 2.1301837572 1.2298621656 3.4785755089
QE‑SSP/gete/gete.dos.in
1 &DOS
2 outdir = './out/'
3 prefix = 'gete '
4 fildos = 'gete.dos '
5 emin = -5
288 | Productivity Tools
Figure 6.8 Obtaining the k‐path and high‐symmetry points of bulk cubic GeTe
in MaterialsCloud.
6 emax = 25
7 /
Figure 6.9 The high‐symmetry points of bulk cubic GeTe in the scaled units
or crystal_b coordinates given by the SeeK‐path tool in MaterialsCloud.
QE‑SSP/gete/gete.nscfbands.in
1 & CONTROL
2 calculation = 'bands '
3 outdir = './out/'
4 prefix = 'gete '
5 pseudo_dir = './ pseudo/'
6 verbosity = 'high '
7 /
8 &SYSTEM
9 ecutrho = 3.2000000000d+02
10 ecutwfc = 4.0000000000d+01
11 ibrav = 0
12 nat = 2
13 ntyp = 2
14 nbnd = 30 # same as in *. nscfdos.in
15 /
16 & ELECTRONS
17 conv_thr = 4.0000000000d-10
18 mixing_beta = 4.0000000000d-01
19 /
290 | Productivity Tools
20 ATOMIC_SPECIES
21 Ge 72.64 ge_pbesol_v1 .4. uspp.F.UPF
22 Te 127.6 te_pbesol_v1.uspp.F.UPF
23 ATOMIC_POSITIONS crystal
24 Ge 0.0000000000 0.0000000000 0.0000000000
25 Te 0.5000000000 0.5000000000 0.5000000000
26 CELL_PARAMETERS angstrom
27 4.2603675143 0.0000000000 0.0000000000
28 2.1301837572 3.6895864968 0.0000000000
29 2.1301837572 1.2298621656 3.4785755089
30 K_POINTS crystal_b # was automatic in SCF file
31 6
32 0.5000000000 0.2500000000 0.7500000000 50 !W
33 0.5000000000 0.5000000000 0.5000000000 30 !L
34 0.6250000000 0.2500000000 0.6250000000 20 !U
35 0.5000000000 0.0000000000 0.5000000000 50 !X
36 0.0000000000 0.0000000000 0.0000000000 50 !G
37 0.3750000000 0.3750000000 0.7500000000 1 !K
QE‑SSP/gete/gete.bands.in
1 &BANDS
2 outdir = './out/'
3 prefix = 'gete '
4 filband = 'gete.bands '
5 /
Having created the above file, we have completed all the steps in
the preparation of Quantum ESPRESSO input files for the DOS and
band structure calculations.
In addition to the SCF, NSCF, and DOS or band processing input files,
we may want to generate the Wannier90 input file directly from the
provided CIF file. For this purpose, there is a useful tool in the forms
Quantum ESPRESSO Course for Solid‐State Physics | 291
of Python script named as cif2qewan, which can be downloaded by
the following command:
$ cd cif2qewan
$ cp ~/QE -SSP/ Ge1Te1_ICSD_638014 .cif .
$ python cif2qewan.py Ge1Te1_ICSD_638014 .cif
cif2qewan.toml
With the above set of command lines, we can obtain the Quan‐
tum ESPRESSO and Wannier90 input files automatically (scf.in,
nscf.in, pw2wan.in, and pwscf.win). In order to understand
and run these files, please check Sec. 3.6.1. Make sure that the
pseudo_dir in \cif2qewan.toml file has been edited properly to
be the location of “PS Library” pseudopotentials as explained in
Sec. 3.1.6, and the name of the pseudopotential of each atom in
pp_psl_rrkj.csv file is correct.
The most basic way that we have explained in Chapters 2 and 3 for
running the Quantum ESPRESSO calculations is using the following
Linux commands typed in the terminal:
292 | Productivity Tools
Table 6.1 File‐ and directory‐related commands in Linux, where file1 or dir1
are the name of file or directory, respectively.
Command Purpose
pwd show present directory
....................................................................................................................................................................
$ pwd
The pwd command displays the absolute (full) path of our present
working directory. Below is an example of pwd execution along with
its output:
$ pwd
/home/quantum/QE -SSP
$ cd
$ pwd
/home/quantum
$ cd QE -SSP
or
$ cd QE -SSP/
For better clarity of the relative path, we may also add a period
character (.) before the directory name we want to enter since
the period character emphasizes the path of the current working
296 | Productivity Tools
$ cd ./QE -SSP/
$ cd /home/quantum/QE -SSP
$ cd ~/QE -SSP
$ cd ../
$ cd ../../
Quantum ESPRESSO Course for Solid‐State Physics | 297
Note that if the directory that we want to change has spaces in
its name, we should either surround the path with quotes (single or
double quotes):
$ ls
gete.bands.in gete.nscfdos.in
gete.dos.in gete.scf.in
gete.nscfbands.in pseudo
$ ls gete
gete.bands.in gete.nscfdos.in
gete.dos.in gete.scf.in
gete.nscfbands.in pseudo
298 | Productivity Tools
We see that the default output of the ls command shows only the
names of the files and directories. We can use the -l option to print
files in a “long” listing format:
$ ls -l
$ ls -l
total 24K
-rw -r--r-- 1 user group 74 Jul 10 06:39 file1
-rw -r--r-- 1 user group 100 Jul 10 06:39 file2
-rw ------- 1 user group 1.2K Jul 10 06:39 file3
-rw ------- 1 user group 892 Jul 10 06:39 file4
-rw ------- 1 user group 829 Jul 10 06:39 file5
drwxr -xr -x 2 user group 4.0K Jul 10 06:39 dir1
drwxr -xr -x 1 user group 4.0K Jul 10 06:39 dir2
The first line of the output is the total size of files (and directory)
in the directory. Then, the next lines list the full information of the
files and directories. The leftmost character in each line will be either
a dash (-), which represents a file or a letter (d), which represents
a directory. The set of following letters (rwxrwxrwx) indicate the
permissions given to the particular file or directory.
The permissions are broken into three roles, i.e. rwx, rwx, and
rwx for a user, a group, and others, respectively. The name of the user
and the group are listed in user and group columns in the output
of ls -al. Then, each role has specified permission for the file or
directory in the following order: read (r), write (w), and execute (x).
A dash character (-) is shown if there is no permission to read, write,
or execute. Therefore, the following line,
Similarly,
Note that the ls command does not list the hidden files by
default. A hidden file is any file that begins with a period character.
To display all files, including the hidden files, use the -a option:
ls -al
$ ls -al
$ ls -al /home/quantum/QE -SSP/gete
The first command above will list all files and folders in the current
working directory, including the hidden files therein in the long
listing format. Meanwhile, the second command will output the list
for the specified directory. An example of the output is given below:
$ ls -al
total 36K
drwxr -xr -x 3 user group 4.0K Jul 10 06:39 .
drwxr -xr -x 5 user group 4.0K Jul 10 06:39 ..
-rw -r--r-- 1 user group 74 Jul 10 06:39 file1
-rw -r--r-- 1 user group 100 Jul 10 06:39 file2
-rw ------- 1 user group 1.2K Jul 10 06:39 file3
-rw ------- 1 user group 892 Jul 10 06:39 file4
300 | Productivity Tools
If the file already exists, the touch command will change the file’s
last access and modification times to the current time. If the file does
not exist yet, the command creates an empty new-file.txt with a
zero‐byte size.
$ ln -s sourcefile symboliclink
$ cd ~
$ ln -s ~/QE -SSP/gete gete
The above commands will give a symlink named as gete in the home
directory so we can access gete folder directly from the home
directory.
$ cd ~/ gete
$ ln -s ~/QE -SSP/gete
$ mkdir /tmp/newdirectory
mkdir can take one or more directory names as its arguments. If the
argument is a directory name, without the full path, the new directory
is created in the present working directory.
To create directories along with the parent directory, we should
use the -p option. For example,
$ mkdir -p QE -SSP/gete/out
$ rm file.txt
$ rm -i file.txt
...
rm: remove regular empty file 'file.txt'?
Only when we answer the above question with “y” (without quotes),
the rm command will perform.
Use the ‐d option to remove one or more empty directories:
$ rm -d dirname
$ rm -rf dirname
The -f option does not prompt the user, and rm -rf command
deletes all files and arguments under the directory without asking
remove or not. It should be carefully noted that the rm -rf ./*
command is considered as the most dangerous command in Linux
since all files and folders under the current working directory will
Quantum ESPRESSO Course for Solid‐State Physics | 303
completely disappear. Imagine if we execute this command in our
home directory, all our works will be wiped out, and if we execute
the command in the root (/), our OS can be broken. Therefore, the
readers should not use this command or use it after to back up the
contents of the directory before performing the command.
$ cp file.txt filebackup.txt
$ mv file.txt /tmp
304 | Productivity Tools
$ mv file.txt filerenamed.txt
$ cat gete.scf.in
$ more gete.scf.out
head (or tail): Display the beginning (or end) of the file
The head and tail commands can be used to display the beginning
and end of the file, respectively. By default, the number of lines shown
by the two commands is limited to 10. If we want to change the
number of lines printed on the screen, we can use the -n option
followed by the desired number. We can try performing the head and
tail commands on the gete.scf.out file. For the head command,
the example execution along with the resulting output is as follows:
$ head -n 13 gete.scf.out
For the tail command (this is useful to see JOB DONE message):
$ tail -n 8 gete.scf.out
PWSCF : 10m33 .89s CPU 16m 8.73s WALL
=-----------------------------------------------=
JOB DONE.
306 | Productivity Tools
=-----------------------------------------------=
The above command will print the last 8 lines of gete.scf.out into
the log.txt file, which can be opened later by a text editor.
The tail command also has a special option: -f, which is useful
for displaying the contents of the file as it grows, starting with the
last 10 lines. We can benefit from such an option, in particular, when
monitoring a job running on the background. For example, when we
execute
$ tail -f gete.scf.out
Command Purpose
whoami display who we are logged in as
....................................................................................................................................................................
$ whoami
quantum
$ w
Quantum ESPRESSO Course for Solid‐State Physics | 309
$ du -h
392K ./out
1.0M ./ pseudo
6.7M .
From the output above, we can see the usage of that directory
(denoted by the period symbol) is around 6.7 MB. The output
obtained by the readers may be different with this example
depending on the readers’ activities in the directory.
$ whereis app
$ whereis pw.x
pw:
/apps/tools/qe -6.4.1/ bin/pw.x
/home/ahma079/qe -6.7/ bin/pw.x
The output indicates that there are two locations of pw.x. The first
one is /apps/tools/qe-6.4.1/bin/pw.x and the second one is
/home/ahma079/qe-6.7/bin/pw.x. This output can vary depending
on the system used by the readers.
The question now is, “Which app is used by default when we call
pw.x?” To get the answer, we can execute the which command with
the following format:
$ which app
$ which pw.x
/apps/tools/qe -6.4.1/ bin/pw.x
^Z
[1]+ Stopped pw.x < gete.scf.in > gete.scf.out
Now the command prompt returns to the original state with the
cursor waiting for our action. What happens to the pw.x calculation
above is that it stops until we enter either the fg command or bg
command. First, let us try the fg command. Its execution with the
example output is as follows.
$ fg
pw.x < gete.scf.in > gete.scf.out
$ bg
[1]+ pw.x < gete.scf.in > gete.scf.out &
$ top
Tasks: 16 total , ...
...
PID USER ... COMMAND
850 quantum ... pw.x
...
Quantum ESPRESSO Course for Solid‐State Physics | 313
There are actually further details (…) of the output but we just show
the number of tasks, the process number (PID), the user who has
the process, and the command related to the process. Note that the
output of PID and USER may vary. However, if we run the same pw.x
command for the SCF calculation of GeTe, we can see it somewhere
below the COMMAND information of the top output. If we think that the
calculation was run wrongly or if we just want to stop it, we should
write down its PID, press q to exit the top command mode, and run
the kill command as will be explained next.
$ kill pid
$ kill 850
$ kill -9 850
The second way is by using the killall with the process (or
command) name as its argument.
$ killall pw.x
$ lscpu
Architecture: x86_64
CPU op -mode(s): 32-bit , 64-bit
Byte Order: Little Endian
CPU(s): 8
On -line CPU(s) list: 0-7
Thread(s) per core: 2
Core(s) per socket: 4
Socket(s): 1
...
$ free
total used free ...
Mem: 11337068 143932 10964292 ...
Swap: 3145728 0 3145728 ...
nohup: No hang up
When exiting the terminal or shell of a Linux System, all running
processes are usually terminated or hang up. If we want to keep the
processes running when exiting the terminal, we can use the nohup
command. It should be noted, however, that if the computer (not only
a physical computer but also a virtual machine or WSL) is shut down,
the nohup command and the processes to which it is assigned will
also stop. Therefore, if we want to run very long Quantum ESPRESSO
jobs, in addition to using the nohup command, we should make sure
that the computer can survive running within the calculation time
frame.
The nohup command can be combined with the background
operator & as follows
$ exit
The majority of the codes that come with the Quantum ESPRESSO
package can run in parallel similarly to the above way.
In a more sophisticated HPC system (let us say a supercom‐
puter), the result of lscpu should not be taken into consideration
for deciding the value of the -np option because when we login
to that supercomputer, we are only brought to the “control node”
that will distribute our job to the other nodes dedicated for parallel
processing. In this case, we should consult the system administrator
or the manual of the supercomputer for the best setting of the CPU
number we can use for our calculations. On the other hand, if we are
just running Quantum ESPRESSO calculations on our own laptop, it
is safe to use a smaller CPU number than what is indicated by lscpu.
We should also be aware that plane‐wave DFT calculations (like
that in Quantum ESPRESSO) do not scale linearly. Our calculation
will get faster to a certain value of N, after which the calculation
will become slow if we add more parallel processes. The situation
can vary depending on the PC hardware (chipset, memory size, etc.),
software, and type of calculation we are doing, but usually we will see
a reduction in the speedup after around 8–32 processes.
If we do the generic way of mpirun above (such as mpirun -np
8 pw.x < gete.scf.in > gete.scf.out), we accept the default
“strategy” for parallelizing the calculation. Using this default strategy
for the pw.x command, in most cases, is already a good choice.
However, if we are not satisfied with the speedup given by the default
Quantum ESPRESSO Course for Solid‐State Physics | 317
strategy, we can use other parallelization schemes as discussed
in detail in the Quantum ESPRESSO documentation, especially for
pw.x (“PWscf”) and ph.x (“PHonon”) applications, in the sections of
performances and parallelism, respectively.
Command Purpose
grep pattern files search for pattern in files
....................................................................................................................................................................
6.2.5 Searching
For searching files and patterns, Linux has some powerful commands
such as grep, find, and locate. For our purpose of analyzing and
understanding Quantum ESPRESSO output files, it is sufficient for
us to focus on the grep command. The usage summary of the grep
command is given in Table. 6.3. Let us practice this command for
obtaining some important information about our calculation.
The most basic syntax for grep is
In this case, since the SCF calculation performed iterations for finding
the minimum total energy, we see the words “total energy” several
times according to the number of iterations. The one with the
exclamation mark (!) is the final result of total energy. Therefore,
if we just want to know this final result, we can change the search
pattern to the exclamation mark. The grep execution with its output
in this case is
$ grep ! gete.scf.out
! total energy = -239.63636259 Ry
We can practice this syntax for finding the Fermi energy in the gete
working directory (~/QE-SSP/gete/).
$ grep -r Fermi .
Remember that the period (.) sign denotes the current directory,
so that the above syntax will look for all occurrences of “Fermi” in
all files of the current directory. The example of the output is given
below:
320 | Productivity Tools
The output tells us that the grep command could find the word
“Fermi” in the gete.nscfdos.out and gete.dos files in the current
directory.
The grep command can also be combined with the other Linux
commands through the “piping” method with the pipe (|) character.
The syntax is
where command is any command or app. For example, we can list the
current directory and find any file or directory which contains “out”
pattern. The example execution and some parts of its output are given
below.
There are some useful shortcuts available in the Linux terminal. The
summary of these keyboard shortcuts is given in Table. 6.4. You can
try them one by one following the table.
Shortcut Purpose
CTRL+c halt the current command
....................................................................................................................................................................
about how to use the command. The way we write it in the terminal
is as follows.
$ command --help
For example,
322 | Productivity Tools
$ cp --help
$ man command_name
For example, to open the man page of the command to change the
directory (cd), you would type
$ man cd
To navigate the man pages, use the Arrow, Page Up, and Page
Down keys. We can also press the Enter key to move one line at a time,
the Space bar to move to the next screen, and the letter b key to go
one screen back. To exit the man page, press the letter q key.
Shell scripts are useful for automating processes and running many
jobs simultaneously because shell scripts allow us to string together
sets of commands. With the shell scripts, we could write a script
that runs calculations, parses the important results to another file,
and even generates a plot. In this section, we will learn basic shell
scripting to be able to automate Quantum ESPRESSO calculations and
to submit batch jobs as a background job with nohup command.
There are often several different shells installed on a Linux
system, and bash is typically the default one in most Linux
distributions. The bash shell is running in the terminal, and it
interprets and executes the commands we type in.
Quantum ESPRESSO Course for Solid‐State Physics | 323
6.3.1 Environment
$ echo $PATH
processes rather than local to the script. Note when we are naming
the variable we want to assign a value to we do not use a $. Also we
just want to add a directory to the existing PATH so we add :$PATH to
keep all the existing directories.
We can practice updating the PATH variable to include not
only the ~/bin directory but also the Quantum ESPRESSO binary
directory. Suppose we install Quantum ESPRESSO latest version
(q-e) under ~/opt directory, we can add the export statements for
the PATH variable as follows:
~/.bashrc
export PATH="~/bin:$PATH"
export PATH="~/opt/q-e/bin:$PATH"
When you login the next time, you do not need this.
6.3.2 Scripting
$ nano helloworld.sh
QE‑SSP/scripts/helloworld.sh
#!/ bin/bash
# Simple script that outputs " Hello World !"
Quantum ESPRESSO Course for Solid‐State Physics | 325
The first line above tells the system what command is used to
interpret the contents of the script. This line should always begin
with #! and then the path to the executable. For bash the executable
will almost always be /bin/bash. In the script, comments are
represented by #. Linux commands can be typed the same as we
would enter them to the terminal. After writing the script, we will
need to make the script executable with chmod u+x hello.sh.
$ ./ helloworld.sh
Hello World!
$ mv helloworld.sh ~/bin/
$ helloworld.sh
Hello World!
QE‑SSP/scripts/helloworld‑rev.sh
#!/ bin/bash
# Example using variables .
var1="Hello"
var2="World!"
326 | Productivity Tools
We can also read a value from standard input using the read
command as follows:
QE‑SSP/scripts/enterinput.sh
#!/ bin/bash
# Example showing how the read command is used.
Command substitution
Command substitution allows the result of a command to replace the
command itself, acting much like a variable in practice. For example:
QE‑SSP/scripts/substitute.sh
#!/ bin/bash
# Example of command substitution .
Conditional statements
The conditional if statements can be used in bash, and many types
of conditional tests are possible. We can test if the value stored in a
variable equals something as follows:
QE‑SSP/scripts/if1.sh
#!/ bin/bash
# Example using if statements
QE‑SSP/scripts/if2.sh
#!/ bin/bash
# Check if the directory " tmpdir " exists , and if not
, create it , then check
# if the file " tmpdir / testfile " exists , and if not ,
create it.
dirname="tmpdir"
filename="testfile"
if [ ! -d "$dirname" ]
# The "!" is logical negation
# -d tests that a file exists and is a directory .
328 | Productivity Tools
then
mkdir $dirname
fi
if [ ! -f "$dirname/$filename" ]
# -d tests that a file exists and is a regular file
(i.e., not a directory ).
then
touch "$dirname/$filename"
fi
QE‑SSP/scripts/arithmetictest.sh
#!/ bin/bash
# Example showing the use of the (( )) construction
var1 =4
var2 =5
var3 =8
Looping
In bash shell, we can use the for loops to iterate a variable over a
range of values. For example:
Quantum ESPRESSO Course for Solid‐State Physics | 329
QE‑SSP/scripts/loop1.sh
#!/ bin/bash
for i in 1 2 3
do
echo "Iteration number $i"
done
QE‑SSP/scripts/loop2.sh
#!/ bin/bash
In the above example, basename is a tool to strip the suffix from a file.
Here we use it to construct a new filename with the .dat extension.
Bash also supports numeric loop ranges using the syntax
{start..finish..increment). For example:
#!/ bin/bash
# Example of for loop with numeric increments
for i in {0..6..2}
do
echo "The value of i is $i"
done
330 | Productivity Tools
In addition to the for loops, the bash shell also supports while
loops, where a set of commands are repeated continuously when a
given condition is true. For example:
QE‑SSP/scripts/while1.sh
#!/ bin/bash
# Example of while loop in bash
counter =0
QE‑SSP/scripts/while2.sh
#!/ bin/bash
# Example of while loop reading lines from a file
linenumber =0
Arguments
We can pass arguments to bash scripts from the command line. For
example:
QE‑SSP/scripts/argument1.sh
#!/ bin/bash
# Example using command line arguments .
Quantum ESPRESSO Course for Solid‐State Physics | 331
QE‑SSP/scripts/argument2.sh
#!/ bin/bash
# Example using a for loop to iterate over command
line arguments .
for arg in $@
do
echo $arg
done
Now let us make a simple script for submitting batch jobs of Quantum
ESPRESSO calculations. The script runEnDOS.sh includes Quantum
ESPRESSO calculations of SCF, NSCF for obtaining DOS, and NSCF for
obtaining the energy dispersion.
QE‑SSP/gete/runEnDOS.sh
#!/ bin/bash
# Simple QE job script
# This script can be re -used for other materials
# Later we just need to change
# the following 3 variables :
# first is the Quantum ESPRESSO binary location
export QEBIN =/ home/quantum/opt/q-e/bin/
# second is the number of CPUs
export NCPU =4
# last is the prefix of the calculations
export PREFIX=gete
# QE binaries that we want to use
export PW=$QEBIN/pw.x
export DOS=$QEBIN/dos.x
332 | Productivity Tools
export BND=$QEBIN/bands.x
# define variables for calculations of
# SCF , DOS , and energy dispersion
export SCF=$PREFIX.scf
export NSCFD=$PREFIX.nscfdos
export POSTD=$PREFIX.dos
export NSCFB=$PREFIX.nscfbands
export POSTB=$PREFIX.bands
echo "start SCF"
mpirun -np $NCPU $PW < $SCF.in > $SCF.out
echo "SCF done"
echo "start NSCF DOS"
mpirun -np $NCPU $PW < $NSCFD.in > $NSCFD.out
echo "NSCF DOS done"
echo "start Postprocessing DOS"
mpirun -np $NCPU $DOS < $POSTD.in > $SCF.out
echo "Postprocessing DOS done"
echo "start NSCF BANDS"
mpirun -np $NCPU $PW < $NSCFB.in > $NSCFB.out
echo "NSCF BANDS done"
echo "start Postprocessing BANDS"
mpirun -np $NCPU $BND < $POSTB.in > $POSB.out
echo "Postprocessing BANDS done"
echo "All jobs are finished"
QE‑SSP/scripts/plot.ipynb
respectively. We can load this data using the loadtxt function from
the NumPy package that we imported as np.
QE‑SSP/scripts/plot.ipynb
QE‑SSP/scripts/plot.ipynb
Figure 6.12 A default plot of gr.dos, which is not suitable for a scientific paper.
We can see that the plot shown in Fig. 6.12 has a few meanings
already, but we notice that the default Matplotlib settings do not give
the figure with sufficient publication quality. For example, there is no
axis label, the font size is too small, and we may also not like the tick
marks outside the main frame. As we change some of the parameters
described in the next subsections, we can get a better looking plot
that is good enough for publication as shown in Fig. 6.13.
It is safe to say that the appearance of Fig. 6.13 is better than
Fig. 6.12 in many ways. In particular, the plot now has axis labels
for variables along with their units, readable fonts with a reasonable
size, standard colors, clear tick marks for both the major and minor
336 | Productivity Tools
There are at least three general parameters that we need to set at the
beginning of the Python plotting script: (1) font, (2) font size, and (3)
axis line width. These three are essentially global parameters defined
at the beginning of the Python script that we do not change later, i.e.,
we do not have to explicitly set font/sizes for each label down the
line). We can add the following code before we generate any figures:
QE‑SSP/scripts/plot.ipynb
QE‑SSP/scripts/plot.ipynb
1 import matplotlib.font_manager as fm
2
3 # Collect all the available fonts
4 font_names = [f.name for f in fm.fontManager.ttflist
]
5 print(font_names)
We should not forget to add a label to each axis. The tick widths and
lengths should be edited to match our axis parameters. If we have
minor ticks, we can also specify by a parameter which='mirror'.
QE‑SSP/scripts/plot.ipynb
After all basic settings are ready, we can plot the dataset and give
annotations or legends to the plots.
QE‑SSP/scripts/plot.ipynb
At the very last, do not forget to save the plot. The PDF format
is recommended because it is a vector image that will not lose the
resolution.
When the readers want to plot many figures with the same style,
making your Matplotlib style is helpful to save time. You do not
need to copy/paste settings in Matplotlib whenever you create a new
figure. By importing a style file, you can ensure consistency while still
maintaining the ability to override settings as you wish within the
individual scripts.
The readers can find a style file for this book at
QE-SSP/matplotlib/sci.mplstyle. The sci.mplstyle file
can open and edit by any text editor for example, nano:
$ nano sci.mplstyle
QE‑SSP/matplotlib/sci.mplstyle
1 # Figure properties
2 figure.figsize : 6.5, 5
3
4 # Font properties
5 font.family : Arial
6 font.size : 22
7
8 #### LaTeX
9 mathtext.default : regular
10
11 # Axes properties
12 axes.titlesize : medium # fontsize of the axes
title
13 axes.titlepad : 10 # pad between axes and
title in points
14 axes.titleweight : normal # font weight for axes
title
15 axes.linewidth : 2 # edge linewidth
340 | Productivity Tools
The readers can modify the file freely. The details for each setting
can be found at https://matplotlib.org/stable/tutorials/
Quantum ESPRESSO Course for Solid‐State Physics | 341
Figure 6.14 Energy dispersion (left) and density of state (DOS) (right) of
GeTe.
QE‑SSP/gete/plot‑gete.ipynb
8 k = np.unique(data[:, 0])
9 bands = np.reshape(data[:, 1], (-1, len(k)))
10 # Load data from gete.dos
11 ener , dos , idos = np.loadtxt('gete.dos', unpack=True
)
12
13 # Fermi energy from ./out /*. xml :
14 HOMO = 2.635710373417968e-1
15 LUMO = 2.758546355831248e-1
16 unitE = 27.2114
17 E_F = (HOMO+LUMO)* unitE /2
18
19 # Set high -symmetry points from *. nscf.in file
20 rW = k[0]; rL = k[50]; rU = k[80]; rX = k[100]; rG =
k[150]; rK = k[200]
21
22 # Create figure object
23 fig = plt.figure(figsize =(6, 3))
24 # Add x and y-axes
25 axE = fig.add_axes ([0.00 , 0.0, 0.60, 1])
26 axD = fig.add_axes ([0.70 , 0.0, 0.30, 1])
27
28 # Plot band structure
29 for band in range(len(bands)):
30 axE.plot(k, bands[band , :]-E_F , c='b')
31 # Plot dotted line at Fermi energy
32 axE.axhline (0, c='gray ', ls=':')
33 # Plot dotted lines at high -symmetry points
34 axE.axvline(rL , c='gray ')
35 axE.axvline(rU , c='gray ')
36 axE.axvline(rX , c='gray ')
37 axE.axvline(rG , c='gray ')
38 # Add labels for high -symmetry points
39 axE.set_xticks ([rW , rL , rU , rX , rG , rK])
40 axE. set_xticklabels (['W', 'L', 'U', 'X', '$\Gamma$ ',
'K'])
41 # Hide x-axis minor ticks
42 axE.tick_params(axis='x', which='both ', length =0)
43 # Set the axis limits
44 axE.set_xlabel('')
45 axE.set_ylabel('Energy (eV)')
46 # Set the axis limits
47 axE.set_xlim(rW , rK)
48 axE.set_ylim (-2.2, 2.2)
49
50 # Plot the DOS
51 axD.plot(dos , ener -E_F , c='r')
52 # Set the axis limits
53 axD.set_xlim (0, 2)
Quantum ESPRESSO Course for Solid‐State Physics | 343
54 axD.set_ylim (-2.2, 2.2)
55 # Add the x label
56 axD.set_xlabel('DOS (state/eV/u.c.)')
57
58 # Save figure to the pdf file
59 plt.savefig('GeTeEnDOS.pdf')
60 # Show figure
61 plt.show ()
Bibliography