ToolBoX Seminar Lecture Notes

1 A Practical Introduction to Density Functional Theory

2 L. Rademaker1*

3 1 Department of Theoretical Physics, University of Geneva, 1211 Geneva, Switzerland

4 * louk.rademaker@gmail.com

5 March 19, 2020

6 Abstract
7 These notes contain a brief practical introduction to doing density functional
8 theory calculations for crystals using the open source Quantum Espresso soft-
9 ware. The level is aimed at graduate students who are studying condensed
10 matter or solid state physics, either theoretical or experimental.

11

12 Contents
13 1 Introduction: What is Density Functional Theory? 2
14 1.1 Born-Oppenheimer approximation 2
15 1.2 Hohenberg-Kohn theory 3
16 1.3 Approximating the functional 3
17 1.4 Kohn-Sham equation 4
18 1.5 Limitations 5

19 2 A Practical Guide 5
20 2.1 Wavefunction basis sets 6
21 2.2 Which functional? 6
22 2.3 Pseudopotentials 6
23 2.4 Quantum ESPRESSO 7
24 2.5 Materials Cloud 7

25 3 Example 1: Silicon 8
26 3.1 Self-consistent field 8
27 3.1.1 Writing the input file 8
28 3.1.2 Running the code 10
29 3.1.3 Reading the output 11
30 3.2 Homework: Find the lattice constant and bulk modulus of silicon 11
31 3.3 Bands 11
32 3.4 Try at home 13

33 4 Example 2: Graphene 13
34 4.1 Compute the band-structure 14

35 5 Example 3: WSe2 16
36 5.1 Relax 16
37 5.2 Spin-orbit coupling 18

38 6 Further reading 18

1
 ToolBoX Seminar Lecture Notes

39 References 19
40

41

42 1 Introduction: What is Density Functional Theory?

43 Any material on earth, whether in crystals, amorphous solids, molecules or yourself, con-
44 sists of nothing else than a bunch of atoms, ions and electrons bound together by electric
45 forces. All these possible forms of matter can be explained by virtue of one simple equation:
46 the many-particle Schrödinger equation,
 
N 2 2 N 2
∂ X h̄ ∂ X e Zi Zj 
ih̄ Φ(r; t) = − + Φ(r; t). (1)
∂t 2mi ∂r2i |ri − rj |
i i<j

47 Here Φ(r; t) is the many-body wavefunction for N particles, where each particle has its
48 own mass mi , charge Zi and position ri . The only interaction is the Coulomb interaction
49 e2 /r.
50 Despite its apparent simplicity, Eq. (1) is notoriously difficult to solve. This is where
51 density functional theory (DFT) comes in. Using a set of reasonable physical approxima-
52 tions we can simplify the many-particle Schrödinger equation to something that we can
53 actually solve numerically.

54 1.1 Born-Oppenheimer approximation

55 The first approximation arises from the physical problem we want to study: the ground
56 state of a collection of interacting ions and electrons. Because even the lightest ion is more
57 than a thousand times heavier than an electron, we will forget about the dynamics of the
58 ions all-together. This is known as the Born-Oppenheimer approximation. We then write
59 the time-independent Schrödinger equation for a collection of N electrons subject to the
60 electric potential created by the fixed ions,
 
N  2 2
 XN 2
X h̄ ∂ e
 − + V (ri ) +  Ψ(r) = E0 Ψ(r) (2)
2m ∂r2i |ri − rj |
i i<j

61 where ri are the positions of the electrons. The potential V (ri ) is created by the charged
62 ions,
X e2 Zj
V (ri ) = − (3)
|ri − Rj |
j

63 where R is the (static) positions of the ions and Zj their charge. Note that the above
64 Hamiltonian – the left hand side of Eq. (2) – contains three terms: the kinetic energy (T ),
65 the potential energy (V ) and the interaction energy (U ).
66 The electronic density is obtained by integrating out all electron degrees of freedom
67 exact one, Z
n(r) = d3 r2 · · · d3 rN |Ψ(r1 · · · rn )|2 . (4)

68 The total potential energy V is just given by the integral over the potential V (r) times
69 the density, Z
V = d3 rV (r)n(r). (5)

2
 ToolBoX Seminar Lecture Notes

70 1.2 Hohenberg-Kohn theory

71 Assume we found a solution of Eq. (2), with ground state energy E0 and a certain electronic
72 density n(r). The strength of the Coulomb interaction and the mass of an electron are
73 constants of nature, so the only input that can possibly influence the electronic density
74 n(r) and the energy E0 of our ground state is our choice of potential V (r). In other words,
75 the ground state energy is a functional of the input potential,

E0 [V (r)] = FE [V (r)] (6)

76 A functional is nothing else than a function whose input is another function; in this case
77 the functional F takes as input the electric potential generated by the ions and outputs
78 the ground state energy based on Eq. (2).
79 At first this results seems counterintuitive. After all, the ground state energy clearly
80 contains the kinetic energy T , the interaction energy U and the potential energy. Only
81 the latter term explicitly depends on the potential. We can thus write the ground state
82 energy in terms of a separate functional for the kinetic and interaction energy, and the
83 potential energy Z
E[n(r)] = FE0 [V (r)] + d3 rV (r)n(r) (7)

84 Hohenberg and Kohn [1] came to the elegant insight that the potential V (r) and elec-
85 tronic density n(r) are conjugate variables. Other conjugate variables you may know are
86 for example pressure and volume in thermodynamics or momentum and position in clas-
87 sical physics. The fact that the potential V (r) and the density n(r) are conjugate means
88 you can equally well describe any solution of Eq. (2) using the potential or the density.
89 Formally known as a Legendre transform (in the same way you go from the Hamiltonian
90 to the Lagrangian formulation of classical mechanics), we can change the functional of Eq. 7
91 to depend on the density n(r) rather than the potential V (r). This is the Hohenberg-Kohn
92 theorem: there exists a universal functional of electronic density, F[n(r)], such that for
93 the correct density n(r) it provides the ground state energy of Eq. (2),
Z
E[n(r)] = F[n(r)] + d3 rV (r)n(r). (8)

94 Knowing this functional, for any given potential V (r) we minimize the right hand side by
95 checking all possible electronic density distributions.
96 There are only two minor problems. We don’t know what this functional looks like.
97 And even if we did, we don’t know how to find the right electronic density.

98 1.3 Approximating the functional

99 The unknown functional F[n(r)] should describe the kinetic and interaction energy of a
100 system described by Eq. 2. Even though we cannot find its exact shape, we can look at
101 its shape in some limiting cases that we can solve.
102 We know that a free homogeneous electron gas with density n has a ground state
103 energy of

2/3
3h̄2 3π 2 5/3
E0 = n0 . (9)
10m
104 For a slowly varying electronic density, we can approximate the kinetic energy contribution
105 to the full functional F[n(r)] as the energy of Eq. (9) evaluated at each point separately,

2/3 Z
3h̄2 3π 2
T0 [n(r)] = d3 r(n(r))5/3 . (10)
10m

3
 ToolBoX Seminar Lecture Notes

106 Furthermore, we know from perturbation theory that the lowest order energy contri-
107 bution from Coulomb interactions is given by the Hartree term,

e2 n(r)n(r0 )
Z
UH [n(r)] = d3 rd3 r0 . (11)
2 |r − r0 |

108 It is natural to write out the full functional as containing the homogeneous electron gas
109 term and the Hartree term. The remaining terms, though still unknown, should be small.
110 This unknown part is conventionally called the exchange-correlation potential Exc [n(r)].
111 The full Hohenberg-Kohn functional, including the potential energy, is thus
Z
EHK [n(r)] = T0 [n(r)] + d3 rV (r)n(r) + UH [n(r)] + Exc [n(r)]. (12)

112 We will later discuss some general choices of exchange-correlation functionals in Sec. 2.2.

113 1.4 Kohn-Sham equation

114 We replaced an intractable problem (solving Eq. (2)) with the task of minimizing an
115 unknown functional F[n(r)] over infinitely many possible electronic densities n(r). In the
116 previous section we already gave some first suggestions for the functional. But once we
117 found it, how to find the right electronic density n(r)?
118 Because the correct density minimizes the functional, we can find the functional by
119 setting it’s derivative to zero,
δF[n(r)]
= 0. (13)
δn(r)
120 Using the functional Eq. (12), we write out1

n(r0 ) 3 0 δExc [n(r)]

Z
δT [n(r)]
+ V (r) + d r + = 0. (14)
δn(r) |r − r0 | δn(r)

121 The idea of Kohn and Sham [2] was to treat this as if it is a single-particle problem. The
122 first term represents the kinetic energy, and the remaining terms form the Kohn-Sham
123 potential
n(r0 ) 3 0 δExc [n(r)]
Z
VKS (r) = V (r) + d r + . (15)
|r − r0 | δn(r)
124 The Kohn-Sham equation is the single-particle Schrödinger equation with the potential
125 given by Eq. (15),
h̄2 ∂ 2
 
− + VKS (r) ψi (r) = i ψi (r). (16)
2m ∂r2
126 We solve these equations numerically, which is tractable because it’s just a linear differen-
127 tial equation. The electronic density is obtained by occupying the N solutions ψi (r) with
128 the lowest energy,
XN
n(r) = |ψi (r)|2 . (17)
i=1

129 Now the electronic density obtained this way can be used to calculate a new Kohn-Sham
130 potential following Eq. (15). We continue this iterative procedure until we reach conver-
131 gence.
1
There is a small subtlety included in this equation: the kinetic component of the functional T [n(r)]
should not be the one obtained for the free homogeneous noninteracting electron gas of Eq. (10), but the
one for a noninteracting gas subject to the Kohn-Sham potential.

4
 ToolBoX Seminar Lecture Notes

132 A final comment is in order: in the above derivation we completely ignored the spin
133 of electrons. Of course, real electrons have spin so that you need that degree of freedom
134 as well. This does not change anything fundamental about how to use the Kohn-Sham
135 equation.

136 1.5 Limitations

137 We have reached the end-goal: using the Born-Oppenheimer approximation, with an ap-
138 propriate choice of functional, we use the Kohn-Sham equations to find the ground state
139 energy and electronic density of a system of interacting electrons and ions. This combi-
140 nation of approximations and techniques is called density functional theory (DFT).
141 Despite is sometimes shaky assumptions, DFT turned out to be a resounding success. A
142 large majority of crystalline materials, many molecules and molecular structures have been
143 explained using DFT. Walter Kohn – the man who was involved in both the Hohenberg-
144 Kohn theory and the Kohn-Sham equations – received the Nobel Prize for DFT in 1998.
145 Mainly because of this success, I assume, you want to learn DFT in this ToolBoX.
146 However, let me briefly shed some clouds over DFT’s success.

147 • A major limitation of DFT is that it is impossible to tell you the size of the errors
148 that exist due to the assumptions. If you get a self-consistent solution, that is nice,
149 but only a comparison with the experimental system will tell you whether it is a
150 good solution.

151 • The electronic properties of many materials can be described using band theory,
152 meaning for every quasimomentum k we have a set of energies n (k). The solutions
153 of the Kohn-Sham equation Eq. (16) are commonly interpreted as these electronic
154 bands. However, it is important to bear in mind that in principle there is no con-
155 nection to the actual electronic energy levels in a material. It just turns out that, in
156 many materials, the Kohn-Sham energies happens to be a good approximation.

157 • As a corollary to the previous limitation: when using DFT for an insulator or semi-
158 conductor, you can also compute the Kohn-Sham energy gap between the highest
159 occupied and lowest unoccupied state. This is not the actual gap of the semicon-
160 ductor or insulator. In practice, it turns out that DFT typically underestimates the
161 real gap.

162 • Note that even if we had the exact functional, solving the corresponding Kohn-Sham
163 equations would not give you the exact solution for the ground state energy.

164 • Because the Kohn-Sham equations describe non-interacting electrons, many materi-
165 als with strong correlations cannot be described using DFT. In general, this is true
166 for materials with partially filled d or f -orbitals; or materials with localized electrons.
167 In particular, applying DFT to the class of high-temperature superconductors such
168 as cuprates, pnictides, and heavy fermions is relatively unsuccessful.

169 • Materials with ground state degeneracy – for example in the case of spontaneous
170 symmetry breaking – are known to be difficult to compute using DFT.

171 • There is some subtlety involved in choosing the right functional. Commonly, func-
172 tionals become widely accepted because of a good overlap with experiments.

173 A good book to learn more about the theoretical side of DFT is Ref. [3].

5
 ToolBoX Seminar Lecture Notes

174 2 A Practical Guide

175 DFT as introduced in the last section is a numerical technique – there are no analytical
176 expansions or solutions. And like many modern numerical techniques, it is better to use
177 an existing developed code than to write your own. There are numerous DFT software
178 packages, both open source and payed. In this section we will outline the most important
179 choices we need to make in order to start running some code.

180 2.1 Wavefunction basis sets

181 The core of any code consists of computing the solutions to the Kohn-Sham equation.
182 Because this is a linear differential equation, we need to choose a basis over which we
183 can expand the Kohn-Sham equation. In this way, we transformed the continuum differ-
184 ential equation into a matrix equation, for which there are many known techniques for
185 diagonalizing it.
186 The two most popular choices of basis are plane waves (PW) and Gaussian type orbitals
187 (GTO). If you are interested in the structure of molecules, the most logical basis set is
188 GTO where you take a certain set of polynomials multiplied by a Gaussian envelope, to
189 make sure the electronic density remains close to the ions.
190 For crystals, on the other hand, the most logical basis set is plane waves ψ(k) = eik·r
191 in a box with periodic boundary conditions. Because this course is aimed at condensed
192 matter physicists, we will use a plane-wave code.
193 A full list of existing DFT codes, with their preferred basis set, is maintained on
194 Wikipedia [4].

195 2.2 Which functional?

196 In Sec. 1.3 we introduced some basic ideas regarding the precise shape of the functional.
197 We separated the kinetic energy of a noninteracting gas and the Hartree-Fock interaction
198 energy from the exchange-correlation functional. Here we will discuss in more detail some
199 possible exchange-correlation functionals that have been proposed, without any attempt
200 at being encyclopedic.
201 It is known – see for example the textbook [5], chapter 5 – that the energy of the
202 homogeneous electron gas can be expanded in powers of the Fermi momentum kF ∼ n1/3 .
203 Kohn and Sham [2] suggested to use these analytical results for the exchange-correlation
204 functional, which is now known as the local density approximation (LDA).
205 A natural next step is to have a functional that not only depends on the density
206 n(r) but also on its derivative ∇n(r). Such functionals are known as generalized gradient
207 approximations (GGA). [6] A popular version of a GGA functional is the Perdew–Burke-
208 Ernzerhof functional (PBE) [7] - its publication is cited more than 100000 times! In this
209 work we will use the PBE functional, since it reproduces experimental band-structures
210 relatively accurate.

211 2.3 Pseudopotentials

212 Are we ready to start computing? Well, not yet. If we are interested in doing DFT
213 for a crystal, we would give our code the size of the unit cell, the type of atoms and
214 their positions. For example, silicon has an fcc crystal structure with two atoms per unit
215 cell. Silicon itself is element number 14, which means there are 28 electrons per unit cell.
216 However, both physically and numerically it is nonsensical to include all 28 electrons in
217 our calculation.

6
 ToolBoX Seminar Lecture Notes

218 Physically speaking, the electronic configuration of silicon is Ne 3s2 3p2. The core
219 electrons, given by the electronic configuration of neon, are completely irrelevant in the
220 physics of binding a silicon crystal. Numerically, the fact that a bare atomic core has
221 a diverging potential 1/r creates a lot problems. Both these problems can be solved by
222 putting a pseudopotential at the position of the silicon atoms. A pseudopotential is a
223 smeared-out potential that includes the charge of the physically irrelevant core electrons,
224 that is therefore more numerically stable than the diverging 1/r potential of a bare atomic
225 core. Using a pseudopotential for silicon, for example, means you now do DFT with only
226 the four outermost electrons, known as the valence electrons.
227 Like with the choice of functional, there are many different ways to compute a pseu-
228 dopotential. In fact, most pseudopotentials are tailored to work with certain functionals,
229 so in these notes we will use pseudopotentials that work well with the PBE functional.

230 2.4 Quantum ESPRESSO

231 Now we are ready to select a DFT implementation. In these notes we will use the open-
232 source plane-wave DFT code Quantum ESPRESSO (http://www.quantum-espresso.
233 org/). Its development was started in Trieste, Italy, but has by now many contributors
234 from all around the world. If you ever use Quantum ESPRESSO scientifically, make
235 sure you explicitly acknowledge the code and cite their original journal publications [8, 9].
236 Our first task is to install Quantum ESPRESSO on your own computer. The sim-
237 ulations we will do in these lecture notes are light enough that they can be done on a
238 standard laptop.
239 The full source package of Quantum ESPRESSO can be downloaded from
240 https://www.quantum-espresso.org/download
241 If you are comfortable doing so, you can install Quantum ESPRESSO using the source
242 package. Otherwise, on the above webpage you will also find links for stable binaries for
243 typical platforms such as Windows. For Mac OS X I recommend using MacPorts (https:
244 //www.macports.org/), which has a port called quantum-espresso. It automatically
245 takes care of dependencies, such as OpenMPI and Fortran libraries.

246 2.5 Materials Cloud

247 We will use the pseudopotential libraries collected by Materials Cloud, [10, 11] an online
248 tool and repository for doing DFT calculations developed by the EPFL and the ETH.
249 On their website, they have collected various pseudopotentials, benchmarked them, and
250 selected for each element the best choice of pseudopotential.
251 We will download their collection of pseudopotentials, which are all computed for the
252 PBE functional.
253 1. Go to the website https://www.materialscloud.org/, and navigate to the Dis-
254 cover page using the top menu.

255 2. Click on Standard solid-state pseudopotentials (SSSP).

256 3. You will see a periodic table (see Fig. 1). By clicking on the button that reads
257 Pseudo, you will download a tar.gz file containing all the pseudopotentials.

258 4. Once you downloaded the library, look into the folder. You will see files with names
259 like Si.pbe-n-rrkjus_psl.1.0.0.UPF and C.pbe-n-kjpaw_psl.1.0.0.UPF. These
260 are the pseudopotentials for silicon and carbon, respectively. Directly after that, you
261 can see that these pseudopotentials are computed for the PBE functional.

7
 ToolBoX Seminar Lecture Notes

Figure 1: The Standard solid-state pseudopotentials library from Materials Cloud [10, 11]
contains for all elements a choice of pseudopotential, which can be downloaded on https:
//www.materialscloud.org/.

262 3 Example 1: Silicon

263 Silicon is one of the most abundant materials on earth, and one of the most used in modern
264 technology. It is therefore logical that we start with silicon as the first material we are
265 going to study.
266 To start off, create a folder where we will do all of our calculations, with three sub-
267 folders: silicon, out and pseudo. Then move the pseudopotential file for silicon that we
268 downloaded from Materials Cloud into the pseudo directory.

269 3.1 Self-consistent field

270 3.1.1 Writing the input file
271 Quantum ESPRESSO works with input files. These are text-only files that contain all
272 the parameters that you want to give to your code. We will now build together an input
273 file for computing the ground state of a silicon crystal. Note that all possible parameters
274 are summarized in the file INPUT_PW.txt or INPUT_PW.html that came with the source
275 package. If you are unsure on what to write, please ***
276 So let’s get coding! Open your favorite plain text editor, and create a new file called
277 silicon.scf.in. The input file we will create is structured with cards, of the form
278 &NAME ... /.
279 The first card named CONTROL describes calculation parameters.
280
281 & CONTROL
282 calculation = ’ scf ’
283 prefix = ’ silicon ’
284 outdir = ’ ../ out / ’
285 pseudo_dir = ’ ../ pseudo / ’
286 tprnfor = . true .
287 verbosity = ’ high ’

8
 ToolBoX Seminar Lecture Notes

288
289
/

290 • calculation decides what calculation we will do. Choose ’scf’, which is short for
291 self-consistent field calculation. It is the default option of Quantum ESPRESSO.
292 In the scf-mode, the self-consistent loop described in Sec. 1 is done to find the ground
293 state energy and the electronic density.
294 • prefix is the name of your set of calculations. This allows you to later use the
295 output of this calculation in further calculations, for example to obtain the band
296 structure.
297 • outdir is the directory where all the detailed output files will be put, in files of the
298 form prefix.xml and prefix.save, with the exception of the command line output.
299 • pseudo_dir is the directory where you put the pseudopotentials.
300 • The last two flags are optional. We set tprnfor true, which means the code will
301 calculate and output the forces acting on the ions. This is a useful way to check
302 that your proposed crystal structure is stable. We also set verbosity to high, which
303 means the output file will contain a lot of extra information, such as the calculated
304 symmetry representations of the crystal. I can really advise to run both with high
305 and low verbosity and to compare the output files.

306 The next card describes the system we are studying. What kind of lattice, how many
307 atoms, and what are the cut-offs for our plane wave basis.
308
309 & SYSTEM
310 ibrav = 2
311 celldm ( 1 ) = 10 . 2
312 nat = 2
313 ntyp = 1
314 occupations = ’ fixed ’
315 ecutwfc = 30
316 ecutrho = 120
317 /
318 & ELECTRONS
319
320
/

321 • ibrav selects the type of Bravais unit cell. In the case of silicon, we choose 2, which
322 is the number for face-centered cubic (fcc). In this crystal structure, the lattice
323 vectors are

a1 = a/2(−1, 0, 1); a2 = a/2(0, 1, 1); a3 = a/2(−1, 1, 0), (18)

324 where a is the lattice constant given in celldm(1), see also Fig. 2, left side. The
325 parameter celldm(1) is given in atomic units, meaning Bohr radii. A full list of
326 possible Bravais lattices can be found in the file INPUT_PW.txt.
327 • nat is the total number of atoms per unit cell, ntyp is the total number of different
328 atoms per unit cell. In fcc silicon, there is only one type (silicon) but two atoms per
329 unit cell.
330 • occupations tells the code how to occupy the computed Kohn-Sham states. fixed
331 means we just occupy all the states below the Fermi level and empty all the states
332 above. Choose this option for insulators; however, for metals (such as graphene in
333 the next section) we need a different option.

9
 ToolBoX Seminar Lecture Notes

334 • Recall that Quantum ESPRESSO uses a plane-wave basis to describe the wave-
335 functions. We should tell the system how many plane waves should be included in
336 the calculation. This is characterized by a kinetic energy cut-off ecutwfc, which
2
|k|2
337 implies that we take all the plane waves with momenta such that h̄ 2m ≤ ecutwfc.
338 The parameter is given in Rydberg units (1 Ry = 13.606 eV). On the Materials
339 Cloud website, see Fig. 1, there is a suggested minimum wavefunction cut-off for
340 every pseudopotential, which for silicon is 30 Ry. Making it larger makes the code
341 slower, but should give a more accurate result.

342 • The code not only stores the plane waves of the electronic states, but also directly
343 the density distribution. Therefore a second cut-off is necessary, ecutrho. Because
344 density scales as the square of wavefunctions, and the kinetic energy scales as the
345 square of the momentum, we need to include a density cut-off at least four times
346 ecutwfc. For our case, exactly four times suffices.

347 • The card &ELECTRONS allows us to tell the system how to solve the Kohn-Sham
348 equation. In our case, we just use the default values. However, we need to include
349 the card!

350 The third part of the input contains explicit information about the atoms: their pseu-
351 dopotentials and positions.
352
353 ATOMIC_SPECIES
354 Si 28 . 086 Si . pbe - n - rrkjus_psl . 1 . 0 . 0 . UPF
355 ATOMIC_POSITIONS alat
356 Si 0 . 00 0 . 00 0 . 00
357
358
Si 0 . 25 0 . 25 0 . 25

359 • Below the line ATOMIC_SPECIES you list all the types of atoms that exist in your unit
360 cell. In our case, it is just silicon. For each type of atom, you give its name (Si), its
361 atomic mass in units of u (28.086), and the relevant pseudopotential file name.

362 • After ATOMIC_POSITIONS you list the positions of all the atoms. The flag alat
363 means the positions are given in cartesian coordinates in units of the lattice param-
364 eter a (celldm(1)). Alternatively, one can use angstrom (cartesian coordinates in
365 Angstrom) or crystal (multiples of the primitive lattice vectors). Here we have a
366 silicon atom at the origin, and a second silicon atom at a( 41 , 14 , 14 ).

367 In the final part of the input we tell the code at which momentum points we will do
368 the calculation.
369
370 K_POINTS automatic
371
372
6 6 6 1 1 1

373 The simplest option is to choose the flag automatic, which generates a Monkhorst-Pack
374 grid [12]. The first three numbers indicate the number of k-points in each of the three
375 directions ( 6 6 6). The last three provide a possible offset in each direction: 0 means
376 no offset and thus the inclusion of high-symmetry points like Γ; 1 means that you place
377 the momentum points exactly in between the points generated by 0. A finer momentum
378 mesh is generated if you choose 1, so that is what we will choose typically.

379 3.1.2 Running the code

380 Congratulations, you have now written your first DFT input file! To run it, simply use
381 the following command on the command line:

10
 ToolBoX Seminar Lecture Notes

382 pw.x -in silicon.scf.in > silicon.scf.out

383 The program pw.x is the main component of Quantum ESPRESSO. It takes the input
384 file silicon.scf.in and does the self-consistent calculation of the ground state energy
385 and density. On most modern computers, it should not take longer than a few seconds for
386 a system as simple as silicon.

387 3.1.3 Reading the output

388 The heavy part of the output is sent to the output directory ../out/. The human-readable
389 part is saved in the file silicon.scf.out. Let’s read through it together.
390 The beginning of the output file lists the properties of this calculation, many of them
391 were given by your input file. Some of them were implicit, yet properly picked up. For
392 example, on line 45 you can see that the program will use the PBE exchange-correlation
393 functional. On line 39 we read that we are going to use 4 Kohn-Sham states, with spin
394 degeneracy this corresponds to 8 electrons per unit cell. After that follows properties of
395 the crystal structure, including its symmetries, and the list of momentum points in our
396 grid.
397 Just before the actual calculation starts, the program estimates the amount of mem-
398 ory needed for this calculation, in the line Estimated max dynamical RAM per process.
399 This might be useful to know for more complicated crystal structures, but for silicon this
400 poses no problems.
401 Between the lines Self-consistent Calculation and End of self-consistent calculation the
402 self-consistent DFT loop is iterated until we have reached convergence. After this, the file
403 contains for each momentum point the resulting energies of the Kohn-Sham states. This
404 is directly followed by:
405
406 highest occupied level ( ev ) : 6 . 2506
407
408
409
! total energy = - 22 . 83862400 Ry

410 This is the main output of a DFT calculation: the ground state energy. For clarity, the
411 code also includes the energy of the highest occupied level. This is not exactly the same as
412 the Fermi level, because there might be occupied states with a higher energy at momentum
413 points that were not included in your momentum grid.

414 3.2 Homework: Find the lattice constant and bulk modulus of silicon
415 a. The ground state energy by itself is not a measurable quantity. However, one of the
416 ideas of DFT is that we can find the lattice constant of a crystal by calculating the ground
417 state energy as a function of the input lattice constant, E0 (a). The predicted actual lattice
418 constant is where this function is minimal. So can you predict with your DFT code the
419 lattice constant of silicon?
420 Hint: Write a code that automatically generates input files with different values of the
421 lattice constant a. Then extract the total energy for every value of a.
422 b. In the previous exercise you calculated the function E0 (a). The change of energy
423 under uniform compression is characterized by the bulk modulus. What is the value of the
424 bulk modulus you calculated?

425 3.3 Bands

426 Strictly speaking, the Kohn-Sham energies do not correspond to anything physical. How-
427 ever, it turns out that they are a pretty good approximation to the electron band energies

11
 ToolBoX Seminar Lecture Notes

Silicon (FCC) band structure

10

Dispersion [eV]
0

-5

L Γ X W Γ

Figure 2: Left: Crystal structure of face-centered cubic silicon. Right: Band structure
of silicon as computed using the steps of Sec. 3.3.

428 of weakly interacting materials. Therefore it is instructive to calculate the Kohn-Sham

429 energies for many points in the Brillouin zone, to get a sense of silicons band structure.
430 A calculation of the bands can only be done after you finished a self-consistent field
431 calculation with the same parameters . This is because the bands calculation takes the
432 electronic density (and thus the Kohn-Sham potential) obtained in an scf calculation,
433 and recomputes the Kohn-Sham energies for a new set of chosen momentum points. A
434 good starting point for the bands calculation input file is therefore to copy the scf input
435 file to a new file called silicon.bands.in. If you want to write your own bands file from
436 scratch, make sure it has the same prefix as the previous self-consistent field calculation.
437 At three points, we will change this new input file. First, we need to tell the code that
438 we want to calculate the band structure, so replace the calculation line with this:
439
440
441
calculation = ’ bands ’

442 For an insulator or semiconductor like silicon, the default setting is that only occupied
443 bands will be computed. If we are interested in the band gap, we also need to calculate
444 the unoccupied bands. To achieve this, add to the &SYSTEM card a line that says we want
445 to calculate 8 bands:
446
447
448
nbnd = 8

449 The final change we need is to the momentum point grid. A customary way to visualize a
450 bandstructure is to choose a path in the Brillouin zone and compute the bands along this
451 path. This can be done by the command tpiba_b, which means that the K_POINTS are in
452 units of 2π/a. The subscript _b indicates that we can define a path for a bandstructure
453 calculated. For a path L – Γ – X – W – Γ, we replace the old K_POINTS card by
454
455 K_POINTS tpiba_b
456 # tpiba_b = k - points in units of 2pi /a , in format for band calculation
457 # number of k - points ( use high - symmetry points only )
458 5
459 # kx , ky , kz , n . of points between this and next one
460 0 . 5 0 . 5 0 . 5 20
461 0 . 0 0 . 0 0 . 0 30
462 0 . 0 0 . 0 1 . 0 30
463 0 . 0 1 . 0 1 . 0 30
464
465
0.0 0.0 0.0 0

466 Notice that here we added some comment lines (the one starting with #). These will be
467 ignored by the code, and can be useful for our own understanding of the input files. Here,

12
 ToolBoX Seminar Lecture Notes

468 the comments explain us that we have a path with 5 momentum points, and that we have
469 20 momentum points in between L and Γ, and so forth.
470 We can run the bands calculation with this command,

471 pw.x -in silicon.bands.in > silicon.bands.out

472 It might take a few seconds longer than the scf calculation, because we have more mo-
473 mentum points.
474 The output file silicon.bands.out starts out with listing the parameters of the cal-
475 culation. After the line Band Structure Calculation it will calculate the Kohn-Sham
476 energies of each desired momentum point. At the end of the calculation, you will find lines
477 looking like this:
478
479 End of band structure calculation
480
481 k = 0 . 5000 0 . 5000 0 . 5000 ( 754 PWs ) bands ( ev ) :
482
483 - 3 . 3153 - 0 . 6624 5 . 1803 5 . 1803 7 . 9984 9 . 7300 9 . 7300 14 . 1551
484
485 k = 0 . 4750 0 . 4750 0 . 4750 ( 748 PWs ) bands ( ev ) :
486
487
488
- 3 . 3519 - 0 . 6101 5 . 1849 5 . 1849 8 . 0036 9 . 7355 9 . 7355 14 . 1635

489 It signals the end of the calculation, and then it will list for each momentum point all
490 the Kohn-Sham energies in eV. The momentum points themselves are given in units of 2π a
491 where a is the lattice constant (celldm(1)).
492 The Quantum ESPRESSO code itself comes with a set of post-processing tools, one
493 of them allows you to plot the band-structure thus calculated. You can also import the
494 output file into your favorite tool (Python, Mathematica, Matlab, GNUplot) and plot it
495 there. The resulting band-structure is shown in Fig. 2, right. Notice it is very similar to
496 the actual band-structure!
497 The band gap calculated using this code is about 0.8 eV, significantly smaller than the
498 actual band gap in silicon of about 1.1 eV. This is common among DFT calculations of
499 semiconductors. On the other hand, many qualitative features – including the fact that
500 silicon has an indirect band gap – are reproduced with our simple calculation!

501 3.4 Try at home

502 Congratulations, you have successfully predicted the properties of a material, purely from
503 first principles! To get comfortable with this technique, try some other three-dimensional
504 semiconductors, and calculate their lattice constant and band-structure.

505 1. First try other zincblende structures like C-diamond, β-SiC, and GaAs.

506 2. Next, calculate a different crystal structure: rock-salt NaCl. Can you find in
507 PW_input.txt how to program its simple cubic structure?

508 3. Finally, study hexagonal α-SiC. Which one has a lower ground state energy, α- or
509 β-SiC?

510 4 Example 2: Graphene

511 The study of silicon in the previous section allows you to calculate crystal and electronic
512 band structures of any three-dimensional semiconductor. In the remainder of these notes,

13
 ToolBoX Seminar Lecture Notes

513 we will switch gears and introduce a few new concepts: we discuss how to deal with metals,
514 with two-dimensional materials, and how we can optimize the structure within a unit cell,
515 and how to have more complicated unit cells.
516 Graphene is a perfect material to do this. It is, as you know, a two-dimensional
517 semi-metal with a hexagonal unit cell.

518 4.1 Compute the band-structure

519 As with the silicon, we need to write input files for a self-consistent calculation first. Make
520 a new folder graphene and start an input file named graphene.scf.in. The first card,
521 &CONTROL, is the same as in the silicon case but with only the prefix changed to graphene.
522 The &SYSTEM card will have some significant changes. Let’s write it out in its full
523 totality,
524
525 & SYSTEM
526 assume_isolated = ’ 2D ’
527 ibrav = 4
528 celldm ( 1 ) = 4 . 65
529 celldm ( 3 ) = 6
530 nat = 2
531 ntyp = 1
532 occupations = ’ smearing ’
533 smearing = ’ mv ’
534 degauss = 1 . 5000000000d - 02
535 ecutwfc = 45
536 ecutrho = 180
537 /
538 & ELECTRONS
539
540
/

541 • In the standard DFT implementation, we have periodic boundary conditions in all
542 three direction. The command assume_isolated = ’2D’ ensures that there is no
543 periodicity (neither in the charge density nor in the Coulomb interactions) in the
544 z-direction.

545 • Graphene has an hexagonal lattice, which has ibrav = 4. The lattice vectors are
546 given by √
a1 = a(1, 0, 0); a2 = a(− 12 , 23 , 0); a3 = a(0, 0, c/a), (19)
547 where as before a =celldm(1) in Bohr, and c/a =celldm(3) is the ratio between
548 the horizontal and vertical lattice size. The value for celldm(3) should be such
549 that, for our 2d set-up in graphene, the vertical unit cell size should be large than
550 the cut-off of the pseudopotentials – in this case at least 20 Bohr.

551 • In the case of a metal or semimetal, just computing the occupied Kohn-Sham ener-
552 gies is very numerically unstable. Tiny changes can lead to different shapes of the
553 Fermi surface. It is therefore necessary to smear the occupations of the Kohn-Sham
554 states, which is ensured by setting occupations = ’smearing’. We then also need
555 to set which type of smearing we will use; here we opted for Marzari-Vanderbilt (mv)
556 smearing [13], with a width set by degauss in Ry units.

557 • Notice we changed the wavefunction and density cut-offs.

558 The last part of the input file tells us where the carbon atoms are going to be, and
559 our choice of momentum points. The only subtleties are in the placement of the carbon

14
 ToolBoX Seminar Lecture Notes

Graphene band structure

Dispersion [eV]
-5

- 10

- 15
Γ M K Γ

Figure 3: Left: The honeycomb crystal structure of graphene. Right: Band structure of
graphene with the Dirac cone clearly visible.

560 atoms, see Fig. 3, left, and in the choice of K_POINTS: because we have a two-dimensional
561 system we only need one momentum point in the z-direction.
562
563 ATOMIC_SPECIES
564 C 12 . 0107 C . pbe - n - kjpaw_psl . 1 . 0 . 0 . UPF
565 ATOMIC_POSITIONS alat
566 C 0 . 000000 0 . 000000 0 . 000000
567 C 0 . 000000 0 . 5773503 0 . 000000
568 K_POINTS automatic
569
570
9 9 1 1 1 1

571 Run the self-consistent field calculation by

572 pw.x -in graphene.scf.in > graphene.scf.out

573 The next step is to make our bands calculation input file. Like before, we can just copy
574 the scf input file to a new file graphene.bands.in. Make sure you change the type of
575 calculation, and the list of K_POINTS. For the latter, I suggest a path Γ – M – K –
576 Γ. Because K and M are particularly easily expressed in terms of the reciprocal lattice
577 vectors, we write the K_POINTS in units crystal_b:
578
579 K_POINTS crystal_b
580 4
581 0 . 000000 0 . 000000 0 . 000000 40
582 0 . 500000 0 . 000000 0 . 000000 20
583 0 . 333333 0 . 333333 0 . 000000 40
584
585
0 . 000000 0 . 000000 0 . 000000 0

586 In the silicon case we explicitly asked the code to compute more than just the occupied
587 bands, using nbnd. Because we are computing a (semi)metal, using the smearing flag, the
588 code automatically calculates some unoccupied bands as well. We do not need to specify
589 the number of bands nbnd.
590 As before, the bands calculation can now be run by the command

591 pw.x -in graphene.bands.in > graphene.bands.out

592 The resulting band-structure, with the characteristic Dirac cone at the Fermi level, can
593 be seen in Fig. 3.

15
 ToolBoX Seminar Lecture Notes

WSe2 band structure

-3

Dispersion [eV]
-4

-5

-6

-7
Γ M K Γ

Figure 4: Left: Crystal structure of WSe2 . Right: Final band-structure with spin-orbit
coupling, with a sizeable spin-orbit splitting at K of ∆vSOC = 0.4 eV.

594 5 Example 3: WSe2

595 The third and final materials that allows us to learn some new features of DFT is the
596 two-dimensional material WSe2 . In-plane it has a honeycomb lattice structure, with on
597 one sublattice the W atoms, and on the other sublattice two Se atoms, displaced in the
598 positive/negative z-direction, as shown in Fig. 4, left. Using lattice relaxation calculations,
599 we will be able to find the exact displacement of the Se atoms. Furthermore, WSe2 is a
600 semiconductor with sizeable spin-orbit coupling, and we will show how to include that.

601 5.1 Relax

602 In a relaxation calculation, the DFT code not only computes the ground state energy but
603 also the derivative of the energy with respect to atomic displacements. This corresponds
604 to the forces acting on each atom. If your initial guess of atomic positions is not stable,
605 there will be nonzero forces. The code will suggest a new set of atomic positions based
606 on the direction of those forces. By repeating this until you have no more forces acting
607 on the atoms, you have relaxed the structure and minimized the ground state energy. We
608 will use this feature to calculate the position of the Se atoms in monolayer WSe2 .
609 As before, we start by making a new folder wse2 with in there an input file, which we
610 will call wse2.relax.in. We will calculate the position of the Se atoms, using a relax
611 calculation. The first card of the input file therefore contains the lines
612
613 calculation = ’ relax ’
614
615
prefix = ’ wse2 ’

616 In the &CONTROL card we can also include a force convergence threshold forc_conv_thr,
617 which determines how close to zero we want the final forces to be. In our simple calculation
618 we only need to use the default value, so we do not need include it in our input file.
619 The remainder of the input file looks like this:
620
621 & SYSTEM
622 assume_isolated = ’ 2D ’
623 ibrav = 0
624 nat = 3
625 ntyp = 2
626 occupations = ’ fixed ’
627 ecutwfc = 30
628 ecutrho = 120
629 /

16
 ToolBoX Seminar Lecture Notes

630 & ELECTRONS

631 /
632 & IONS
633 /
634 ATOMIC_SPECIES
635 W 183 . 840 W_pbe_v1 . 2 . uspp . F . UPF
636 Se 78 . 960 Se_pbe_v1 . uspp . F . UPF
637 ATOMIC_POSITIONS angstrom
638 W 0 . 000000 0 . 000000 0 . 000000 0 0 0
639 Se 0 . 000000 1 . 919689645 1 . 500000 0 0 1
640 Se 0 . 000000 1 . 919689645 - 1 . 500000 0 0 1
641 K_POINTS automatic
642 8 8 1 1 1 1
643 CELL_PARAMETERS angstrom
644 3 . 32500000 0 . 0000000000 0 . 0000000000
645 - 1 . 66250000 2 . 8795344676 0 . 0000000000
646
647
0 . 00000000 0 . 0000000000 32 . 0000000000

648 • In the &SYSTEM card, we reverted back to fixed occupations since WSe2 is a semi-
649 conductor. Notice how we changed the cut-offs and the number of atoms and types
650 of atoms.

651 • Because we are interested in the atomic positions, it is worthwhile to write out the
652 lattice vectors and the initial atomic positions explicitly in units of Angstrom. We
653 can do so by selecting ibrav = 0, meaning we have a free form of the unit cell.
654 We should then explicitly write out the three unit vectors in a new card called
655 CELL_PARAMETERS.

656 • We are interested in finding the z-position of the Se atoms. We know already that
657 their in-plane coordinates are given by the honeycomb lattice, which are given as
658 the second and third column of the lines after ATOMIC_POSITIONS angstrom. In
659 the third column we put the z-position. We set W at z = 0, and we guess an
660 initial distance of the Se atoms at z = ±1.5 A. The last three numbers indicate
661 which atomic position coordinates we will relax . For the Se atoms, 0 0 1 means
662 we keep the x, y coordinates fixed, and will minimize the ground state energy with
663 respect to the z-coordinate.

664 • We do need to include a new card &IONS, where we can specify the properties of
665 the atomic displacements the code will do. Having this card empty just means we
666 choose the default values, but for a relax calculation it has to be there!

667 Run the pw.x code as usual. In the output file wse2.relax.out you can see the two
668 self-consistent loops. Given a set of atomic positions, the ground state energy and forces
669 are calculated. If the forces are larger than the threshold, a new set of atomic positions is
670 proposed after the line saying ATOMIC_POSITIONS. After a few iterations, we have reached
671 convergence and we find the following lines containing the final coordinates:
672
673 Begin final coordinates
674
675 ATOMIC_POSITIONS ( angstrom )
676 W 0 . 000000000 0 . 000000000 0 . 000000000 0 0 0
677 Se 0 . 000000000 1 . 919689645 1 . 678711660 0 0 1
678 Se 0 . 000000000 1 . 919689645 - 1 . 678711660 0 0 1
679
680
End final coordinates

681 We predict the distance between the two Se atoms to be 3.36 Å.

17
 ToolBoX Seminar Lecture Notes

682 5.2 Spin-orbit coupling

683 In the above lattice relaxation calculation we did not take into account spin-orbit coupling.
684 In general, spin-orbit coupling becomes more important the heavier the element is, which
685 in our case applies to the tungsten (W). Because spin-orbit will likely not influence the
686 position of the Se atoms, we take the atomic positions from the previous calculations and
687 do the standard scf followed by bands to calculate the bandstructure of WSe2 . Only this
688 time, we will have spin-orbit coupling.
689 Make the wse2.scf.in and wse2.bands.in input files. You can combine the elements
690 from the wse2.relax.in and graphene.bands.in. Make sure you copy the atomic posi-
691 tions from the previous relax output file into our new input files. To turn on spin-orbit
692 coupling, we need to add the following two lines to the &SYSTEM card:
693
694 lspinorb = . true .
695
696
noncolin = . true .

697 The first term turns on spin-orbit coupling, and the second allows for noncollinear spins
698 (so not only up and down but also superpositions). Additionally, we need to have a
699 fully relativistic pseudopotential to study spin-orbit coupling . In the W pseudopotential
700 we have used so far, you can find the following line:
701
702
703
The Pseudo was generated with a Scalar - Relativistic Calculation

704 Because Scalar-Relativistic implies no spin-orbit coupling, we need to find a new pseu-
705 dopotential that is fully relativistic! Many different types of pseudopotentials can be down-
706 loaded from the Quantum ESPRESSO website. Go to https://www.quantum-espresso.
707 org/pseudopotentials/ps-library/, and download a full relativistic ultra-soft pseu-
708 dopotential (USPP) for tungsten (W) that works with the PBE functional. After that,
709 update the line in the input files where you give the pseudopotential:
710
711 ATOMIC_SPECIES
712
713
W 183 . 840 W . rel - pbe - spn - rrkjus_psl . 1 . 0 . 0 . UPF

714 Finally, the amount of valence electron per W is 14 and per Se is 6, meaning the code
715 computes 26 electrons. If you want also to see the conduction bands, I suggest putting
716 nbnd = 32 or higher.
717 As before, run the scf first, followed by a bands calculation. You can check in the
718 output files that the number of Kohn-Sham energies is now equal to the number of elec-
719 trons. Before we included spin-orbit coupling, the spin degeneracy meant we just needed
720 half the amount of Kohn-Sham energies.
721 The final band structure is shown in Fig. 4, right. As before, the band gap (here
722 about 1.3 eV) is smaller than experimentally detected (1.7 eV in monolayers). Notably,
723 the valence band at the K point is split due to the spin-orbit coupling, with a splitting of
724 ∆vSOC = 0.4 eV, comparable to what is measured in experiments. [14]

725 6 Further reading

726 You now have learned the basics of how to compute crystal structures and electronic
727 bands using density functional theory, implemented in the plane-wave tool Quantum
728 ESPRESSO. But there is much more to DFT than just this. There are some additional
729 tools that we haven’t discussed, such as the possibility to compute phonon dispersions,
730 Raman or optical spectra, and wannierization. We also haven’t looked into modern de-
731 velopments of functionals, such as the inclusion of Van der Waals interactions or strong
732 correlations (LDA+U, GW or DFT+DMFT).

18
 ToolBoX Seminar Lecture Notes

733 Luckily, there are many online courses available that are more in-depth than this short
734 ToolBox.

735 • The Materials Cloud webpage also contains a set of lectures, including notes, ex-
736 ercises and videos, on how to do DFT with Quantum ESPRESSO. You can find
737 them here: https://www.materialscloud.org/learn/.

738 • Many universities have their classes on density functional theory online, for example
739 MIT has https://ocw.mit.edu/courses/materials-science-and-engineering/
740 3-320-atomistic-computer-modeling-of-materials-sma-5107-spring-2005/labs/
741 sections

742 • The source manual of Quantum ESPRESSO named pw_user_guide.pdf, which

743 comes with downloading the source package, contains a lot of information on what
744 you can do with the code. The source package also contains examples on how to use
745 the code, see the folder PW/examples for input files and ideas for pw.x. You can also
746 look at examples of other parts of the code, such as PHonon/examples, that show
747 you how to compute phonon dispersions.

748 Acknowledgements
749 These notes would not exist without the enthusiasm of the ToolBoX organizers João
750 Ferreira and Michael Sonner. I would also like to thank Marco Gibertini for discussions.

751 Funding information L.R. acknowledges funding from the SNSF in the form of an
752 Ambizione grant.

753 References
754 [1] P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phys. Rev. 136, B864
755 (1964).

756 [2] W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Corre-
757 lation Effects, Phys. Rev. 140(4A), A 1133 (1965).

758 [3] G. Giuliani, G. Vignale and Cambridge University Press, Quantum Theory of the
759 Electron Liquid, Masters Series in Physics and Astronomy. Cambridge University
760 Press (2005).

761 [4] Wikipedia contributors, List of quantum chemistry and solid-state physics soft-
762 ware — Wikipedia, the free encyclopedia, https://en.wikipedia.org/w/index.
763 php?title=List_of_quantum_chemistry_and_solid-state_physics_software&
764 oldid=942835156 [Online; accessed 2-March-2020] (2020).

765 [5] G. D. Mahan, Many-Particle Physics, Kluwer Academic, New York, 3rd ed. edn.
766 (2000).

767 [6] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J.

768 Singh and C. Fiolhais, Atoms, molecules, solids, and surfaces: Applications of the
769 generalized gradient approximation for exchange and correlation, Phys. Rev. B 46,
770 6671 (1992), doi:10.1103/PhysRevB.46.6671.

19
 ToolBoX Seminar Lecture Notes

771 [7] J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made
772 simple, Phys. Rev. Lett. 77, 3865 (1996), doi:10.1103/PhysRevLett.77.3865.

773 [8] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli,

774 G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli et al., Quantum
775 espresso: a modular and open-source software project for quantum simulations of
776 materials, Journal of Physics: Condensed Matter 21(39), 395502 (19pp) (2009).

777 [9] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra,

778 R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo et al., Ad-
779 vanced capabilities for materials modelling with quantum espresso, Journal of Physics:
780 Condensed Matter 29(46), 465901 (2017).

781 [10] K. Lejaeghere, G. Bihlmayer, T. Bjorkman, P. Blaha, S. Blugel, V. Blum, D. Caliste,

782 I. E. Castelli, S. J. Clark, A. Dal Corso, S. de Gironcoli, T. Deutsch et al., Re-
783 producibility in density functional theory calculations of solids, Science 351(6280),
784 aad3000 (2016).

785 [11] G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet and N. Marzari, Precision and
786 efficiency in solid-state pseudopotential calculations, npj Computational Materials
787 4(1), 72 (2018).

788 [12] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys.
789 Rev. B 13(12), 5188 (1976).

790 [13] N. Marzari, D. Vanderbilt, A. De Vita and M. C. Payne, Thermal Contraction and
791 Disordering of the Al(110) Surface, Phys. Rev. Lett. 82(16), 3296 (1999).

792 [14] Y. Zhang, M. M. Ugeda, C. Jin, S.-F. Shi, A. J. Bradley, A. Martı́n-Recio, H. Ryu,
793 J. Kim, S. Tang, Y. Kim, B. Zhou, C. Hwang et al., Electronic Structure, Surface
794 Doping, and Optical Response in Epitaxial WSe 2Thin Films, Nano Lett. 16(4), 2485
795 (2016).

796 [15] G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao and D. Xiao, Three-band tight-binding model
797 for monolayers of group-VIB transition metal dichalcogenides, Phys. Rev. B 88(8),
798 085433 (2013).

20