1/27

物質科学のための計算数理 II
Numerical Analysis for Material Science II
8th: Density Functional Theory (1)Nov. 30 (Fri)

Lecturer: Mitsuaki Kawamura (河村光晶)

2/27

Schedule (This semester W1, W2)

1. Sep. 28 (Fri) Guidance Y
2. Oct. 5 (Fri) Monte Carlo method O
3. Oct. 12 (Fri) Monte Carlo method O
4. Oct. 19 (Fri) Monte Carlo method O
5. Oct. 26 (Fri) Exact diagonalization Y
6. Nov. 2 (Fri) Exact diagonalization Y
7. Nov. 9 (Fri) Molecular dynamics O (1st report problem will be announced.)
8. Nov. 30 (Fri) Standard DFT code K
9. Dec. 7 (Fri) Density functional theory K
10. Dec. 14 (Fri) Density functional theory K
11. Dec. 21 (Fri) Density functional theory K
12. Dec. 25 (Tue) (2nd)Report problem K (遠隔講義室)
13. Jan. 11 (Fri) Density functional theory K※ Lecturers: Y … Yamaji, K … Kawamura, O… Ohgoe

3/27

Schedule in this section (DFT)1. Nov. 30 (Fri) Standard DFT code
• First-principles calculation and Density functional theory (Lecture)• One-body Schrödinger eq. for periodic system and Bloch theorem (L)• Numerical solution of Kohn-Sham (one-body Schrödinger) eq. (L)• Hands-on DFT code (Tutorial)• Version control system : Git (T)2. Dec. 7 (Fri) Kohn-Sham eq.
• Plane-wave basis and Pseudopotentials (L)• Iterative eigenvalue solution method (L & T)
3. Dec. 14 (Fri) Self-Consistent loop
• Hartree potential (Poisson eq.)
• Brillouin-zone integral (Tetrahedron method)• Visualization (T)4. Dec. 21 (Fri) Total Energy• Coulomb potential for periodic point charge (Ewald sum)
5. Dec. 25 (Tue) (2nd)Report problem K (遠隔講義室)
6. Jan. 11 (Fri) Density functional theory K

4/27

Today’s Schedule
What is the first-principles study
Density functional theory
Kohn-Sham method
Periodic system (Bulk crystal)
Bloch theorem
How to solve Kohn-Sham eq. with computer
DFT codeGit clone
Charge density
Band structure
Fermi surface
Density of states (DOS), Partial DOS

5/27

Electrons in materials
𝑁−෍𝑛=1

𝛁𝑛22𝑁+

1෍
2 ′𝑛,𝑛 =1

𝑁 𝑁atom

1𝑍𝐼+෍ ෍Ψ 𝒓1 , 𝒓2 , ⋯ , 𝒓𝑁 = 𝐸Ψ 𝒓𝑁
𝒓𝑛 − 𝒓𝑛 ′𝒓𝑛 − 𝑹𝐼
𝑛=1 𝐼=1Single H2O molecule : N = 10
Single C60 molecule : N = 360
Bulk crystal, liquid, glass, ・・・ :N > 1024

Both numerically and analyticallyunsolvable.

Information of electronic structure• Charge density• Magnetic moment• Superconductivity• Optical spectrum• etc.

Can we obtain a few quantities such as• Ground-state energy• Ground-state charge density
without solving this equation ?

6/27

Density functional theory
𝑁𝐸 = න 𝑑 3 𝑟1 ⋯ න 𝑑 3 𝑟𝑁 Ψ ∗ 𝒓𝑁−෍𝑛=1

𝐸 𝜌 = න 𝑑 3 𝑟𝜌 𝒓 𝑣 𝒓 + 𝐸𝑢𝑛𝑖𝑣 𝜌𝛁𝑛22𝑁+1෍2 ′𝑛,𝑛 =1

𝑁1+ ෍ 𝑣 𝒓𝑛𝒓𝑛 − 𝒓𝑛 ′Ψ 𝒓𝑁𝑛=1𝑁atom

𝑍𝐼𝑣 𝒓 ≡ ෍𝒓 − 𝑹𝐼
𝐼=1

Universal functional 𝐸𝑢𝑛𝑖𝑣 𝜌 exists.
𝐸𝜌
Functional depends on theoverall shape of 𝜌, e.g.,′1𝜌(𝒓)𝜌(𝒓)33′
𝐸𝐻 𝜌 = ඵ 𝑑 𝑟𝑑 𝑟
2𝒓 − 𝒓′

𝐸𝐺𝑆
𝜌𝐺𝑆

𝜌

Hohenberg-Kohn theorem
Phys. Rev. 136, B864 (1964).

7/27

Kohn-Sham method

W. Kohn and L. J. Sham,
Phys. Rev. 140, A1133 (1965).

′)1
𝜌(𝒓)𝜌(𝒓𝐸𝑢𝑛𝑖𝑣 𝜌 = 𝑇𝐾𝑆 𝜌 + ඵ 𝑑 3 𝑟𝑑3 𝑟 ′+ 𝐸𝑋𝐶 𝜌′
2𝒓−𝒓Kinetic energy of non-interacting system whose charge density is 𝜌(𝒓)
𝛁2−+ 𝑣𝐾𝑆 [𝜌] 𝒓
2𝜑𝑛 𝒓 = 𝜀𝑛 𝜑𝑛 𝒓

𝑁/2𝑇𝐾𝑆 = 2 ෍ න 𝑑 3 𝑟 𝜑𝑛∗ 𝒓𝑛=1

𝛁2−𝜑𝑛 𝒓2′𝜌𝒓𝑣𝐾𝑆 [𝜌] 𝒓 = 𝑣 𝒓 + න 𝑑 3 𝑟 ′+ 𝑣𝑋𝐶 𝜌 𝒓′𝒓−𝒓

𝑁/2

𝜌 𝒓 = 2 ෍ 𝜑𝑛 𝒓

2𝑛=1

𝛿𝐸𝑋𝐶 𝜌𝑣𝑋𝐶 𝜌 𝒓 ≡𝛿𝜌 𝒓

Self-consistent field (SCF)

The exact form of 𝐸𝑋𝐶 𝜌 is not known.

Approximate 𝐸𝑋𝐶 𝜌8/27

Local density approximation

𝐿𝐷𝐴𝐸𝑋𝐶 𝜌 ≈ 𝐸𝑋𝐶𝜌 ≡ න 𝑑3 𝑟𝜌 𝒓 𝜀𝑋𝐶 𝜌 𝒓

𝑣𝑋𝐶 𝜌 𝒓𝐿𝐷𝐴
≈ 𝑣𝑋𝐶
𝜌 𝒓

𝜀𝑋𝐶 𝜌 : Function (not functional)

𝐿𝐷𝐴
𝛿𝐸𝑋𝐶𝜌
≡
= 𝜀𝑋𝐶 𝜌 𝒓𝛿𝜌 𝒓𝑑𝜀𝑋𝐶 𝜌+𝜌 𝒓
ቚ𝑑𝜌
𝜌=𝜌 𝒓

𝜀𝑋𝐶 𝜌 : Use the XC energy of homogeneous-electron gas

Random phase approximation : 𝜀𝑋𝐶 𝜌 = −0.916+ 0.0622 ln 𝑟𝑠 − 0.0938
𝑟𝑠

M. Gell-Mann and K. A. Brueckner, Phys. Rev. 106, 364 (1957).

Monte-Carlo method

4𝜋 3 1
𝑟𝑠 =3
𝜌

0.916
𝜀𝑋𝐶 𝜌 = −+ 0.031 ln 𝑟𝑠 − 0.0480 − 0.0116 𝑟𝑠 + 0.0020𝑟𝑠 ln 𝑟𝑠 (𝑟𝑠 < 1)
𝑟𝑠[Monte-Carlo] D. M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980).
[Parameterize] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).

𝜀𝑋𝐶 𝜌, 𝛁𝜌 : Generalized gradient correction (GGA)9/27

Total energy and force
𝑁𝐸𝑡𝑜𝑡 = න 𝑑 3 𝑟1 ⋯ න 𝑑 3 𝑟𝑁 Ψ ∗ 𝒓𝑁

෍
𝑛=1

𝛻𝑛2
2𝑁+ 𝑉 𝒓𝑛 ; 𝑹𝑁atom+1෍
2 ′𝑛,𝑛 =11
Ψ 𝒓𝑁𝒓 𝑛 − 𝒓𝑛 ′𝑁atom

𝑍𝐼 𝑍𝐽1+ ෍2
𝑹𝐼 − 𝑹𝐽

Nuclear-nuclear interaction

𝐼,𝐽=1Finding the lowest-energy (stable) structure
Structure optimization : One of the most popular usage of DFT

𝜕𝐸𝑡𝑜𝑡𝑭𝐼 ≡= න 𝑑3 𝑟1 ⋯ න 𝑑3 𝑟𝑁 Ψ ∗ 𝒓𝑁
𝜕𝑹𝐼Hellman-Feynmann theorem𝑁𝜕𝑣 𝒓𝑛 ; 𝑹𝑁atom෍𝜕𝑹𝐼Ψ 𝒓𝑁

𝑛=1𝑁atom+ ෍
𝐽=1

𝑍𝐼 𝑍𝐽 (𝑹𝐼 − 𝑹𝐽 )𝑹𝐼 − 𝑹𝐽3

10/27Structure optimization

GGA-PBE functional : J. Perdew, K. Burke, and M. Ernzerhof Phys. Rev. Lett. 77, 3865 (1996)

Other usage : Formation (Cohesive) energy

A + B → C + Energy

𝐸𝐶 − (𝐸𝐴 + 𝐸𝐵 )Magnetic moment (Spin density functional theory) : 𝐸 𝑛↑ , 𝑛↓11/27

Kohn-Sham eq. for periodic system (1)𝛁2−+ 𝑣𝐾𝑆 𝒓
2𝜑 𝒓 = 𝜀𝜑 𝒓

Unit lattice vectors(Not unique)

𝑁×𝑁𝑢𝑐 /2

𝜌 𝒓 =2෍
𝑛=1𝒂1𝑣𝐾𝑆 𝒓 and 𝑣 𝒓 are periodic,i.e., 𝑣𝐾𝑆 𝒓 + 𝑛1 𝒂1 + 𝑛2 𝒂2 + 𝑛3 𝒂3 = 𝑣𝐾𝑆 𝒓

𝜑𝑛 𝒓

2

𝑁 electrons per unit cell
𝑁𝑢𝑐 cells → ∞Equation to solve in the whole region of bulk crystal

𝒂1𝒂2

Unit cell12/27Kohn-Sham eq. for periodic system (2)Bloch’s theorem

𝜑 𝒓 =𝜑(𝒓) can be written as1𝑒 𝑖𝒌⋅𝒓 𝑢𝑛𝒌 (𝒓)𝑁𝑢𝑐𝑢𝑛𝒌 𝒓 + 𝑛1 𝒂1 + 𝑛2 𝒂2 + 𝑛3 𝒂3 = 𝑢𝑛𝒌 𝒓
𝛁+𝒌 2−+ 𝑣𝐾𝑆 𝒓
2𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓

Equation to solve only in the unit cell

𝒂1𝒂2𝒂3

𝜀𝑛𝑘𝒃1 , 𝒃2 , 𝒃3 = 2𝜋𝐼መන 𝑑 3 𝑟𝜌 𝑟 = 𝑁𝑢𝑐∞1𝜌 𝒓 =2
න 𝑑 3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍
𝑛=12 𝜃(𝜀 − 𝜀 )F𝑛𝑘2𝜋 3𝑉𝐵𝑍 ≡Band structure

Unit reciprocal lattice vectors
𝜀𝑛𝑘 and 𝑢𝑛𝒌 𝒓 isperiodic with 𝒃𝛼 in the 𝒌 space.

𝒃1𝒌𝑉𝑢𝑐13/27

One-body energy level𝛁 + 𝑖𝒌 2−+ 𝑣𝐾𝑆 𝒓2e𝑢𝑛𝒌 𝒓 = 𝜀𝑛𝑘 𝑢𝑛𝒌 𝒓e𝜀𝑛𝑘 ≠True quasi-particle energy

Unit [eV]

LDA

GW

Exp.

Si

0.54

1.32

1.17

LiCl

6.07

9.34

9.40

AlP

1.52

2.59

2.50

AlAs

1.25

2.15

2.23

AlSb

0.99

1.64

1.68

GaP

1.82

2.55

2.39

GaAs

0.37

1.22

1.52

InP

0.57

1.44

1.42

X. Zhu and S. Louie, Phys. Rev. B 43, 14142 (1991).

𝜀𝑛𝑘e

Susceptibility

Qualitative analyzation• Angle resolved photo electronic spectroscopy• de Haas-van Alphen oscillation

eInversed magnetic field (1/H)

Band gap𝒌14/27

How to solve Kohn-Sham eq. : Basis𝛁 + 𝑖𝒌 2−+ 𝑣𝐾𝑆 𝒓
2𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓

𝑢𝑛𝒌 𝒓 → 𝑢𝑛𝒌 𝒓𝑖
Huge numerical cost !
Not efficient.

••••

Pseudopotential (plane wave, real space grid, FEM)
Local basis (atomic orbital, Gaussian)Augmented wave (complicated basis)
Mixed basis (PW + Local, augmented + Local) 𝑢 𝒓 = ෍ 𝑢෤
𝑛𝒌𝑛𝒌𝛼 𝜒𝛼 𝒓

𝛼𝑢𝑛𝒌 𝒓 = ෍ 𝑢෤ 𝑛𝒌 𝑮𝐺𝑒 𝑖𝑮⋅𝒓
𝑉𝑢𝑐

Plane waves below the
Cutoff frequency

𝑮+𝒌 2−+ 𝑣ො𝐾𝑆 𝑢෤ 𝑛𝒌 𝑮 = 𝜀 𝑢෤ 𝑛𝒌 𝑮2𝑣ො𝐾𝑆 𝑢෤ 𝑛𝒌 𝑮 = න 𝑑3 𝑟𝑢𝑐𝑒 −𝑖𝑮⋅𝒓
𝑉𝑢𝑐𝑣𝐾𝑆 (𝒓)𝑢𝑛𝒌 𝒓

15/27

Pseudopotential

Electronic state at core level (e.g. 1s state for Fe)
hardly affects the chemical properties(chemical reaction, conduction, superconductivity, etc.)
𝑢𝑛 𝑟
Pseudopotential is generated
in the calculation of single atom
PseudoTrue𝜀𝑛𝑢𝑛,𝑝𝑠 𝑟

valence

𝑟

2p2s

𝑣𝑝𝑠 𝑟
1s
only reproduce theenergy and long-range behaviorof valence state.
Similar to the potential of an ion.

𝑣 𝑟

Then,
used for the calculation ofmolecule, solid, etc.
Transferability• Norm conserving• Ultrasoft• Projector augmented wave

16/27

QuantumESPRESSO:
plane wave and pseudopotential based code

Before the explanation, we should start to install because it takes long time.

mac OSX in ECCS
$ mkdir -p ~/bin/$ echo 'export PATH=$HOME/bin:$PATH' >> ~/.profile$ git clone https://gitlab.com/QEF/q-e.git -b qe-6.3 ~/q-e$ cd ~/q-e/$ ./configure CC=gcc-8; make pw pp; cp bin/* ~/bin/

Prompt. We do not need type this.

17/27

Quantum ESPRESSO

Web page : https://www.quantum-espresso.org/
Developer : P Giannozzi (University of Udine, Italy), et al. (Many)
Language : Fortran (Mainly)
Feature:• Total energy, force, MD, chemical reaction path, charge/spin density• Band structure, density of states, Fermi surfaces• LDA, GGA, Van der Waals, exact exchange functional• MPI+OpenMP parallelism, GPGPU• Phonon and electron-phonon interaction• X-ray absorption spectra, Optical spectra• Many-body perturbation theory (GW, Bethe-Salpeter eq.)
• Wannier function, Berry's phase• Etc.
Manuals : https://www.quantum-espresso.org/resources/users-manual
Pseudopotential library• SSSP (https://www.materialscloud.org/discover/sssp/table/efficiency)• PSLibrary (https://www.quantum-espresso.org/pseudopotentials)

18/27Visualization tools

FermiSurfer : Fermi-surface viewer$ git clone https://scm.osdn.net/gitroot/fermisurfer/fermisurfer.git -b 1.10.1 ~/fermisurfer$ cd ~/fermisurfer/$ touch Makefile.in aclocal.m4 configure$ ./configure --prefix=$HOME; make; make install

VESTA : Crystal-structure viewer (Also charge density, etc.)$ wget https://jp-minerals.org/vesta/archives/3.4.5/VESTA.dmg$ hdiutil mount VESTA.dmg$ cp -rf /Volumes/VESTA/VESTA ~/$ echo 'alias vesta="open -a $HOME/VESTA/VESTA.app"' >> ~/.profile$ source ~/.profile

19/27

Tutorial : MgB2

$ wget https://osdn.net/projects/educational-pwdft/storage/113018.tgz$ tar xzvf 113018.tgz -C ~/$ cd ~/113018/$ vesta mgb2.xsf
 Metal Superconductivity : TC = 39 K
 Stacked Mg triangular lattice and B
honeycomb lattice20/27Structure optimization & input file$ mpirun -np 2 ~/bin/pw.x -npool 2 -in relax.in | tee relax.out

Variation-cell structure relaxation
Number of atoms in the unit cell&CONTROL
calculation = 'vc-relax'
pseudo_dir = './'/&SYSTEMibrav = 0nat = 3ntyp = 2
ecutwfc = 35.000000
ecutrho = 280.000000
occupations = 'tetrahedra_opt'/&ELECTRONS
/
&IONS
/
&CELL/

Number of kinds of atom

Cutoff energy of plane-waves forwave functions [Ry]𝑢𝑛𝒌 𝒓 = ෍ 𝑢෤ 𝑛𝒌 𝑮 𝑒 𝑖𝑮⋅𝒓𝐺

Cutoff depends on the pseudopotentialCutoff energy of plane-waves forthe charge density [Ry]Integration scheme for thenumerical integration in∞1𝜌 𝒓 =2
න 𝑑 3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍
𝑛=12 𝜃(𝜀 − 𝜀 )F𝑛𝑘21/27

Input-file format

CELL_PARAMETERS angstrom
3.200000 0.000000 0.000000-1.600000 2.771281 0.000000
0.000000 0.000000 3.700000
ATOMIC_SPECIES
B 10.811000 b_pbe_v1.4.uspp.F.UPF
Mg 24.305000 Mg.pbe-n-kjpaw_psl.0.3.0.UPF
ATOMIC_POSITIONS crystal
Mg 0.000000 0.000000 0.000000
B 0.333333 0.666667 0.500000
B 0.666667 0.333333 0.500000
K_POINTS automatic
12 12 8 0 0 0

Unit lattice vectors (Å unit)
Element
Relative atomic mass

Pseudopotential file nameAtomic position in thefractional coordinate
𝑹 = 𝑥1 𝒂1 + 𝑥2 𝒂2 + 𝑥3 𝒂3

𝒌-mesh for numerical integration in∞1𝜌 𝒓 =2
න 𝑑3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍2 𝜃(𝜀 − 𝜀 )
F𝑛𝑘

𝑛=1

Dense k-mesh improves accuracy and increases numerical cost.

22/27

Result$ grep volume relax.out
relax.out:
unit-cell volume=221.4261 (a.u.)^3
relax.out:
new unit-cell volume =204.22304 a.u.^3 (30.26273 Ang^3 )relax.out:
new unit-cell volume =193.53924 a.u.^3 (28.67956 Ang^3 )relax.out:
new unit-cell volume =194.50016 a.u.^3 (28.82195 Ang^3 )relax.out:
new unit-cell volume =194.38563 a.u.^3 (28.80498 Ang^3 )relax.out:
new unit-cell volume =194.25661 a.u.^3 (28.78586 Ang^3 )relax.out:
new unit-cell volume =194.25661 a.u.^3 (28.78586 Ang^3 )relax.out:
unit-cell volume=194.2566 (a.u.)^3$ grep -A 3 CELL_PARAMETERS relax.out
CELL_PARAMETERS (angstrom)
3.106677678
0.000000000 -0.000000000-1.553338839
2.690461490
0.000000000
0.000000000 -0.000000000
3.620638395
Exp. (XRD):
a = 3.085 Å, c = 3.523 Å-CELL_PARAMETERS (angstrom)
volume : 29.04 Å3
3.072422829
0.000000000
0.000000000-1.536211414
2.660795917
0.000000000
0.000000000
0.000000000
3.521167214

23/27

Band structure$ mpirun -np 2 ~/bin/pw.x -npool 2 -in nonscf.in$ mpirun -np 2 ~/bin/bands.x -npool 2 -in bands.in$ gnuplot band.gp$ mpirun -np 1 ~/bin/pp.x -npool 1 -in pp.in$ vesta tmp.pp_K001_B00*.xsf

𝑢𝑛𝑘 (𝒓)

AL

Γ

HKM

24/27

Density of states & Fermi surface$ mpirun -np 2 ~/bin/pw.x -npool 2 -in dense.in$ mpirun -np 2 ~/bin/projwfc.x -npool 2 -in pdos.in |tee pdos.out$ gnuplot pdos.gp
∞
1
𝐷 𝜀 =
න 𝑑3 𝑘 ෍ 𝛿(𝜀 − 𝜀𝑛𝑘 )𝑉𝐵𝑍 𝐵𝑍𝑛=1∞

𝐷𝑝 𝜀 =1
න 𝑑 3 𝑘 ෍ 𝛿(𝜀 − 𝜀𝑛𝑘 )𝑉𝐵𝑍 𝐵𝑍𝑛=12× න 𝑑 3 𝑟𝜓𝑝∗ 𝒓 𝜑𝑛𝑘 𝒓

Fermi surface :Isosurface at 𝜀𝑛𝑘 = 𝜀F
$ mpirun -np 1 ~/bin/fermi_proj.x -in proj_pz.in$ mv proj.frmsf pz.frmsf$ mpirun -np 1 ~/bin/fermi_proj.x -in proj_pxy.in$ mv proj.frmsf pxy.frmsf$ fermisurfer pz.frmsf$ fermisurfer pxy.frmsf

25/27

Why do we perform DFT calculation ?
Directly compute (predict)properties of materials.

DFT calculation

Force or total energy.

Compute effective single-body state,
Kohn-Sham orbitals.

Fit into model atomic force field.

Perform classical MD forlarge-size system or log-time scale.

Construct Hubbard model or
local spin model.
Perform Monte-Carlo simulation,
exact diagonalization for moreaccurate calculation of exotic
phenomena.

26/27

Today's summary Basics of DFT calculation
 Hohenberg-Kohn theorem Kohn-Sham method
 Exchange correlation functional
 Total energy and force
 Kohn-Sham eq. in periodic system (bulk crystal) Bloch theorem
 Band structure Numerical method to solve Kohn-Sham eq.
 Plane-wave (and other) basis
 Pseudopotential
 Usage of a DFT code : Quantum ESPRESSO

27/27

Reference books
Electronic Structure: Basic Theory and Practical Methods
Richard M. Martin
Cambridge University Press
Errata : https://es.polytechnique.fr/Electronic_Structure/errata

日本語訳
物質の電子状態 上下
寺倉清之、寺倉郁子、善甫康成
シュプリンガー

固体‐構造と物性
金森 順次郎, 川村 清, 米沢 富美子, 寺倉 清之
岩波
絶版

岩波オンデマンド