1/24

物質科学のための計算数理 II
Numerical Analysis for Material Science II
11th: Density Functional Theory (4)Dec. 21 (Fri)

Lecturer: Mitsuaki Kawamura (河村光晶)

2/24

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) Density functional theory 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) Density functional theory K (遠隔講義室)
13. Jan. 11 (Fri) (2nd)Report problem K※ Lecturers: Y … Yamaji, K … Kawamura, O… Ohgoe

3/24
1.

2.

3.

4.

5.

6.

Schedule in this section (DFT)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)Dec. 7 (Fri) Kohn-Sham eq.
• Plane-wave basis and Pseudopotentials (L)• Iterative eigenvalue solution method (L & T)
Dec. 14 (Fri) Self-Consistent loop
• Hartree potential (Poisson eq.), Atomic potential, XC potential• Update (Broyden's method)• Visualization of grid data (T)Dec. 21 (Fri) Total Energy• Total energy• Brillouin-zone integral (Tetrahedron method)• Coulomb potential for periodic point charge (Ewald sum)Dec. 25 (Tue) Advanced subjects for productive calculation (遠隔講義室)• Generalized gradient correction• Non-local pseudopotentials (Norm-conserving, ultrasoft, PAW)
• ProcedureJan. 11 (Fri) (2nd) Report problem, Question time

4/24

Today’s Schedule

Total energy
Hartree term
Exchange correlation term
Brillouin-zone integration
BasicsSmearing method
Tetrahedron method and its improvement
Fermi surface plot

5/24

Total energy functional for DFT𝑁𝑐𝑁1𝑍𝝉 𝑍𝝉′ 1 − 𝛿𝑅𝑅′ 𝛿𝝉𝝉′𝐸𝑡𝑜𝑡 = 𝐸 𝜌 + ෍ ෍2 ′ ′ 𝝉 + 𝑹 − 𝝉′ − 𝑹′𝑅,𝑅 𝝉,𝝉

𝐸 𝜌 = න 𝑑 3 𝑟𝜌 𝒓 𝑣 𝒓 + 𝐸𝑢𝑛𝑖𝑣 𝜌′1𝜌(𝒓)𝜌(𝒓)33′𝐸𝑢𝑛𝑖𝑣 𝜌 = 𝑇𝐾𝑆 𝜌 + ඵ 𝑑 𝑟𝑑 𝑟+ 𝐸𝑋𝐶 𝜌′
2𝒓−𝒓

𝑑3 𝑘
∗
𝑇𝐾𝑆 𝜌 = 𝑁𝑐 න෍ න 𝑑3 𝑟𝜑𝑛𝒌𝑟
𝐵𝑍 𝑉𝐵𝑍𝑛𝛁2−𝜑𝑛𝒌 𝑟
2𝑑3 𝑘𝒌+𝑮 2= 𝑁𝑐 න෍෍𝑢෤ 𝑛𝒌 𝑮
2
𝐵𝑍 𝑉𝐵𝑍𝑛2

1 𝑖𝒌⋅𝒓𝜑𝑛𝒌 𝒓 =
𝑒 𝑢𝑛𝒌 (𝒓)𝑁𝐶𝑒 𝑖𝑮⋅𝒓
𝑢𝑛𝒌 𝒓 = ෍ 𝑢෤ 𝑛𝒌 𝑮𝑉𝑢𝑐𝐺𝑮Local density approximation

𝐸𝑋𝐶 𝜌 = න 𝑑3 𝑟𝜌 𝑟 𝜀𝑋𝐶 𝜌 𝑟= 𝑁𝑐 න 𝑑3 𝑟𝑟𝜌 𝑟 𝜀𝑋𝐶 𝜌 𝑟𝑢𝑐6/24

Tutorial for total energy$ cd ~/pwdft/$ git checkout master$ git pull$ make clean; make$ cd sample/Al/$ ../../src/pwdft.x < scf.in

Energy per u.c. [eV]
Kinetic energy :22.571740309226847
Hartree energy :0.10954359264556501
rho*V :
14.630190122503379
Ewald eta [Bohr^-1] :1.4221682587100493
Ewald grid (G) :25
25
25
Ewald grid (R) :222
Eward energy :-73.366857049902464
XC energy :-21.802851254893593
Total energy :-57.858234280420255$ bash vol.sh7/24

Hartree and atomic potential
energy.F90

1
ඵ 𝑑 3 𝑟𝑑 3 𝑟 ′2𝜌(𝒓)𝜌(𝒓′ )𝒓 − 𝒓′

1= ෍ 𝜌෤ 𝑮 𝜌෤ 𝑮′ ඵ 𝑑3 𝑟𝑑 3 𝑟 ′2 ′′ ⋅𝒓′𝑖𝑮⋅𝒓𝑖𝑮𝑒 𝑒𝒓 − 𝒓′

𝑮,𝑮𝑖𝑮⋅ 𝒓−𝒓′

1𝑒′ ⋅𝒓′′3′𝑖𝑮+𝑮
3= ෍ 𝜌෤ 𝑮 𝜌෤ 𝑮 න 𝑑 𝑟 𝑒
න𝑑 𝑟
2 ′𝒓 − 𝒓′𝑮,𝑮

𝑁𝑐 𝑉𝑢𝑐4𝜋
𝑁𝑐 𝑉𝑢𝑐=෍ 𝜌෤ 𝑮 𝜌෤ −𝑮=෍ 𝜌෤ 𝑮𝟐2𝑮2𝑮𝑮

Atomic potential termන 𝑑 3 𝑟𝜌 𝒓𝑣 𝒓

2= ෍ 𝜌෤ 𝑮 𝑣෤ 𝑮′𝑁𝑐 𝑉𝑢𝑐 4𝜋 𝑁
lim
𝑮→0 2𝑮 𝟐 𝑉𝑢𝑐

4𝜋𝑮𝟐2′ ⋅𝒓′
3𝑖𝑮⋅𝒓𝑖𝑮න 𝑑 𝑟𝑒 𝑒= 𝑁𝑐 𝑉𝑢𝑐 ෍ 𝜌෤ 𝑮 𝑣෤ −𝑮

𝑮,𝑮′Compensating
divergent term

𝑮,𝑮′

4𝜋 𝑁
− lim 𝑁𝑐 𝑉𝑢𝑐
𝑮→0𝑮 𝟐 𝑉𝑢𝑐

2

8/24

Ion-ion Coulomb repulsion

𝑁𝑐

𝑁𝑁𝑐

𝑁1𝑍𝝉 𝑍𝝉′ 1 − 𝛿𝑅𝑅′ 𝛿𝝉𝝉′𝑁𝑐
𝑍𝝉 𝑍𝝉′ 1 − 𝛿𝑅𝟎 𝛿𝝉𝝉′෍෍
=෍෍′′′+𝑹
2 ′ ′ 𝝉+𝑹−𝝉 −𝑹2𝝉−𝝉′𝑅 𝝉,𝝉

𝑅,𝑅 𝝉,𝝉

𝑁𝑁𝑐𝝉,𝝉

𝑅−𝑖𝑮⋅𝒓𝑒න 𝑑3 𝑟𝒓

𝑁𝑐
1 − 𝛿𝑅𝟎 𝛿𝝉𝝉′=෍ 𝑍𝝉 𝑍𝝉′ ෍2 ′
𝝉 − 𝝉′ + 𝑹
𝑁𝑐

𝑁𝑐3 ′1𝑑 𝑟 −𝑖𝑮⋅𝒓′
14𝜋𝑖𝑮⋅𝒓𝑖𝑮⋅𝒓෍= ෍𝑒න𝑒෍ ′= ෍𝑒𝒓+𝑹𝒓 +𝑹𝑉𝑢𝑐 𝑮 2𝑢𝑐 𝑉𝑢𝑐
𝑅𝐺𝑅𝑁𝐺

𝑁𝑐

𝑁𝑐
4𝜋
𝛿𝑅𝟎 𝛿𝝉𝝉′′𝑖𝑮⋅
𝝉−𝝉=෍ 𝑍𝝉 𝑍𝝉′ ෍ 𝑒−෍
22 ′𝑉𝑢𝑐 𝑮𝝉 − 𝝉′ + 𝑹
𝝉,𝝉

𝐺

𝑁𝑐 𝑉𝑢𝑐 4𝜋 𝑁lim
𝑮→0 2
𝑮 𝟐 𝑉𝑢𝑐

𝑅2Slow decayin sum1
෍𝑹𝑅𝑒 𝑖𝑮⋅𝒓
෍𝑮2𝐺

9/24

Ewald's method (1)12 ∞2𝑥2−𝒓+𝑹=න 𝑑𝑥𝑒
𝒓+𝑹𝜋 0𝑁𝑐

3 −𝑮2𝜋 2 𝑒 4𝑥 2𝑁𝑐

෍𝑒− 𝒓+𝑹 2 𝑥 2𝑹3 ′𝑑𝑟 −𝑖𝑮⋅𝒓′′ +𝑹 2 𝑥 2𝑖𝑮⋅𝒓−𝒓= ෍𝑒න𝑒෍𝑒= ෍ 𝑒 𝑖𝑮⋅𝒓
3𝑉𝑉𝑥𝑢𝑐
𝑢𝑐
𝑢𝑐
𝑮𝑁𝑁𝑐

𝝉,𝝉

𝑅𝑹𝑮

𝑁𝑐
1 − 𝛿𝑅𝟎 𝛿𝝉𝝉′
෍ 𝑍𝝉 𝑍𝝉′ ෍2 ′
𝝉 − 𝝉′ + 𝑹𝑁𝑁𝑐

𝝉,𝝉

𝑅

∞
∞
2 2𝑁𝑐2′=
෍ 𝑍𝝉 𝑍𝝉′ ෍න 𝑑𝑥𝑒 − 𝝉−𝝉 +𝑹 𝑥 − 𝛿𝑅𝟎 𝛿𝝉𝝉′ න 𝑑𝑥2 ′𝜋 𝜂0𝑁2𝜂𝑮𝑁𝑐
2𝜋𝑑𝑥′−+ ෍ 𝑍𝝉 𝑍𝝉′ ෍ 𝑒 𝑖𝑮⋅ 𝝉−𝝉න 3 𝑒 4𝑥 22 ′𝑉𝑢𝑐 0 𝑥𝝉,𝝉𝑮cos 𝑮 ⋅ 𝝉 − 𝝉′𝑮2− 22𝑒 4𝜂

𝑮2

10/24

Ewald's method (2)
𝑁𝑠𝑁𝑐

𝝉,𝝉

𝑅𝑁𝑐
=
෍ 𝑍𝝉 𝑍𝝉′ ෍ 1 − 𝛿𝑅𝟎 𝛿𝝉𝝉′2 ′𝑁

𝜂𝑁𝑐
22− ෍ 𝑍𝝉න 𝑑𝑥2𝜋 0𝝉𝑁𝑮≠𝟎𝑁

𝜂𝑍𝝉2−𝑁𝑐 ෍𝜋𝑁𝑐+ ෍ 𝑍𝝉 𝑍𝝉′ ෍ cos 𝑮 ⋅ 𝝉 − 𝝉′2 ′𝝉,𝝉2 ∞′ +𝑹 2 𝑥 2−
𝝉−𝝉
න 𝑑𝑥𝑒
𝜋 𝜂erfc 𝝉 − 𝝉′ + 𝑹 𝜂𝝉 − 𝝉′ + 𝑹

𝝉

𝑮2− 24𝜋𝑒 4𝜂Complementaryerror function

𝜋𝑁 2𝑁𝑐𝑁𝑐 𝑉𝑢𝑐 4𝜋 𝑁−+ lim
22𝑮→0 2𝑉𝑢𝑐 𝑮
2 𝑉𝑢𝑐 𝜂𝑮 𝟐 𝑉𝑢𝑐2

11/24

Optimize lattice constant($ bash vol.sh)$ grep "Total energy" scf*.out | ¥sed -e 's/scf//g' -e 's/.out://g' > ene.dat$ gnuplot
gnuplot> set xlabel "Lattice constant [Ang]"
gnuplot> set ylabel "Total energy [eV]"gnuplot> plot "ene.dat" u 1:5 w lp

12/24

Brillouin-zone integration𝜌 𝑟 = ෍ ෍ 𝜑𝑛𝒌 𝑟𝒌2𝜃 𝜀 − 𝜀𝐹𝑛𝑘𝑛

Born-von Karman's boundary condition : 𝜑𝑛𝒌 𝑟 + 𝐿 = 𝜑𝑛𝒌 𝑟 ∶ 𝐿 → ∞2𝜋𝐿

The size of each state in the reciprocal space :32𝜋 3=
𝑁𝑐 𝑉𝑢𝑐𝑑3 𝑘𝑑3 𝑘෍⋯ = න⋯ = 𝑁𝑐 න⋯
3𝐵𝑍 2𝜋 /(𝑁𝑐 𝑉𝑢𝑐 )𝐵𝑍 𝑉𝐵𝑍𝒌𝑉𝑢𝑐 =𝒂1 × 𝒂2 ⋅ 𝒂3𝑑3 𝑘
𝜌 𝑟 = 𝑁𝑐 න
෍ 𝜑𝑛𝒌 𝑟𝑉𝐵𝑍 𝐵𝑍
𝑛2𝜋 3≡ 𝑉𝐵𝑍 =𝑉𝑢𝑐2𝜃 𝜀 − 𝜀𝐹𝑛𝑘𝒃1 × 𝒃2 ⋅ 𝒃3𝑑3 𝑘=න෍ 𝑢𝑛𝒌 𝑟𝑉𝐵𝑍 𝐵𝑍𝑛2𝜃 𝜀 − 𝜀𝐹𝑛𝑘

13/24𝒌-integration with 𝜃 𝜀𝑘 or 𝛿 𝜀𝑘Charge density

Density of states𝑑3𝑘
𝜌 𝑟 =න
෍ 𝑢𝑛𝒌 𝑟𝐵𝑍 𝑉𝐵𝑍𝑑3𝑘
𝐷 𝜀 =න෍ 𝛿 𝜀 − 𝜀𝑛𝑘
𝐵𝑍 𝑉𝐵𝑍2𝜃 𝜀 − 𝜀𝐹𝑛𝑘𝑛𝑛Using 𝜃(𝜀𝐹 − 𝜀𝑛𝑘 ) as isToo many k points are required.
Broadening

𝑁𝒌M. Methfessel, A. T. Paxton, Phys. Rev. B 40, 3616 (1989).

Replace 𝜃(𝜀𝐹 − 𝜀𝑛𝑘 ) with a smeared function such as the error function.
Tetrahedron
O. Jepsen and O. K. Andersen, Solid State Commun. 9, 1763 (1971).

14/24

Why we need to reduce 𝑘-points?

E.g. 1 DFT calculation of crystals with plane wave and pseudopotencial.
𝑣𝐾𝑆 𝒓 𝑢𝑛𝒌 𝒓 (for all orbitals) O 𝑁𝒌 × 𝑁𝑏 × 𝑁𝑃𝑊 log 𝑁𝑃𝑊
Subspace diagonalization O 𝑁𝒌 × 𝑁𝑏3
Update 𝑢෤ 𝑛𝒌 𝑮 O 𝑁𝒌 × 𝑁𝑏2 × 𝑁𝑃𝑊

E.g. 2 Hartree-Fock, hybrid DFT, TC, GW𝑣ො𝐹𝑜𝑐𝑘 𝑢𝑛𝒌 𝒓 (for all orbitals) O 𝑁𝒌2 × 𝑁𝑏2 × 𝑁𝑃𝑊 log 𝑁𝑃𝑊※ 𝑁𝑏 ～ The number of electrons
The reduction of k-points without loss of accuracy isthe best way for efficient calculation.

15/24

Broadening method
We can calculate only on discrete k-points.

𝜀𝒌Replace 𝛿(𝑥) with smeared function
such as (Hermite) Gaussian.
∞−1 𝑛𝑥 −𝑥 2 /𝜎2𝛿 𝑥 =෍𝐻2𝑛𝑒𝑛
𝜎𝑛! 4 𝜋𝑛=0

M. Methfessel and A. T. Paxton,Phys. Rev. B 40, 3616 (1989).

This method becomes worse whensmeared function not cover allenergy range.

𝜎𝒌True band structure16/24

Interpolation method
libtetrabz_occ_main@libtetrabz/libtetrabz_occ_mod.F90
𝜀1𝑉𝑇
𝜀2

O. Jepsen and O. K. Andersen, Solid
State Commun. 9, 1763 (1971)

𝜀4𝑉𝐵𝑍
𝑉𝑇 =6𝑁𝑘𝜀3

When we use linear interpolation, the next is satisfied.
If 𝜀1 ≤ 𝜀 ≤ 𝜀2 ≤ 𝜀3 ≤ 𝜀4 ,න 𝜃 𝜀 − 𝜀𝑘𝑇𝜀 − 𝜀1 𝜀 − 𝜀1 𝜀 − 𝜀1= 𝑉𝜀𝑘≤𝜀 ≈ 𝑉𝑇𝜀2 − 𝜀1 𝜀3 − 𝜀1 𝜀4 − 𝜀1𝜕න 𝛿 𝜀 − 𝜀𝑘 = න 𝜃 𝜀 − 𝜀𝑘𝜕𝜀 𝑇𝑇
If 𝜀1 ≤ 𝜀2 ≤ 𝜀 ≤ 𝜀3 ≤ 𝜀4 , ⋯If 𝜀1 ≤ 𝜀2 ≤ 𝜀3 ≤ 𝜀 ≤ 𝜀4 , ⋯If 𝜀1 ≤ 𝜀2 ≤ 𝜀3 ≤ 𝜀4 ≤ 𝜀, ⋯These expressions are notapplicable for higher order (2nd ,3rd ・・・) interpolation !

17/24

Systematic error in tetrahedron method(Traditional) tetrahedron method linearly interpolate thefunction of k (such as wave function, KS energy, etc. ).

TT

If the sign of curvature ofthese function is always
negative, we overestimete.

Not Interpolating
But Fitting

18/24

Optimized tetrahedron method

M. Kawamura et al., Phys. Rev. B 89, 094515 (2014).
Conventional tetrahedron method

Optimized tetrahedron method

TT𝒌𝒌
1. Linear interpolation of 𝜀𝑘 , etc.
2. Tetrahedron method withinterpolated linear function

1. Higher order interpolation of 𝜀𝑘 ,etc.
2. Fit this function into a linear
function in the tetrahedron
3. Tetrahedron method with fittedlinear function

19/24

Library of tetrahedron method

LibTetraBZ

http://libtetrabz.osdn.jp/
MPI + OpenMP (Thread safe)
Fortran/C/C++
3𝑑 𝑘෍න𝛿 𝜀 − 𝜀𝑛𝒌 𝑋𝑛𝒌
𝐵𝑍 𝑉𝐵𝑍𝑛3𝑑 𝑘෍න𝜃 𝜀 − 𝜀𝑛𝒌 𝑋𝑛𝒌
𝐵𝑍 𝑉𝐵𝑍𝑛𝑑3𝑘
෍න𝛿 𝜀𝐹 − 𝜀𝑛𝒌 𝛿 𝜀𝐹 − 𝜀𝑛′ ′𝒌 𝑋𝑛𝑛′ 𝒌
𝑉′ 𝐵𝑍 𝐵𝑍𝑛,𝑛𝑑3𝑘
෍න
𝜃 𝜀𝑛′ ′ 𝒌 − 𝜀𝐹 𝜃 𝜀𝐹 − 𝜀𝑛𝒌 𝑋𝑛𝑛′𝒌
𝑉′ 𝐵𝑍 𝐵𝑍𝑛,𝑛𝑑3 𝑘 𝜃 𝜀𝑛′ ′𝒌 − 𝜀𝐹 𝜃 𝜀𝐹 − 𝜀𝑛𝒌෍න
𝑋𝑛𝑛′𝒌′𝑉𝜀𝑛′𝒌 − 𝜀𝑛𝒌′ 𝐵𝑍 𝐵𝑍𝑛,𝑛𝑑3 𝑘෍න𝜃 𝜀𝐹 − 𝜀𝑛′ ′𝒌 𝜃 𝜀𝐹 − 𝜀𝑛𝒌 𝛿 𝜀𝑛′ ′𝒌 − 𝜀𝑛𝒌 + 𝜔 𝑋𝑛𝑛′ 𝒌 (𝜔)𝑉′ 𝐵𝑍 𝐵𝑍𝑛,𝑛𝑑3 𝑘 𝜃 𝜀𝐹 − 𝜀𝑛′ ′ 𝒌 𝜃 𝜀𝐹 − 𝜀𝑛𝒌෍න𝑋𝑛𝑛′ 𝒌 (𝜔)′𝑉𝜀𝑛′𝒌 − 𝜀𝑛𝒌 + 𝑖𝜔
′ 𝐵𝑍 𝐵𝑍𝑛,𝑛20/24

LibTetraBZ$ cd ~/$ git clone git://git.osdn.net/gitroot/libtetrabz/libtetrabz.git$ cd libtetrabz$ touch Makefile.in aclocal.m4 configure$ ./configure; make$ cd example$ ./dos.x
Which tetrahedron method ?(1 = Linear, 2 = Optimized): 2⏎
k-point mesh ?: 10⏎$ gnuplot
gnuplot> set xlabel "Energy"
gnuplot> set ylabel "DOS"gnuplot> plot "dos.dat" w l

𝜀𝑘 = − cos 𝑘𝑥 − cos 𝑘𝑦 − cos 𝑘𝑧

21/24

Fermi energy with bisection method

libtetrabz_fermieng@libtetrabz/src/libtetrabz_occ_mod.F90

𝑑3 𝑘
𝑁 = ෍න
𝜃 𝜀𝐹 − 𝜀𝑛𝒌
𝐵𝑍 𝑉𝐵𝑍𝑛𝜀𝑢𝑝 = max{𝜀𝑛𝑘 }𝜀𝑙𝑜𝑤 = min{𝜀𝑛𝑘 }do iteration
𝜀𝐹 = 𝜀𝑢𝑝 + 𝜀𝑙𝑜𝑤 /2𝑑3 𝑘෩ = ෍න𝑁
𝜃 𝜀𝐹 − 𝜀𝑛𝒌
𝐵𝑍 𝑉𝐵𝑍𝑛෩ < 𝑁 is small → Finishif 𝑁෩<𝑁else if 𝑁𝜀𝑙𝑜𝑤 = 𝜀𝐹
෩>𝑁else if 𝑁
𝜀𝑢𝑝 = 𝜀𝐹end ifend do iteration

22/24

Iso surface with tetrahedron method

FermiSurfer

Input file format
http://fermisurfer.osdn.jp/en/_build/html/fermisf_input_en.html

23/24

Report problem 3(1) Modify dos.F90 or dos_c.c to compute DOS with Gaussian broadening
method.
𝑁𝑘3𝑑 𝑘1
1 − 𝜀−𝜀𝒌 2𝜎න𝛿 𝜀 − 𝜀𝒌 ≈෍𝑒𝑁𝑘𝜎 𝜋
𝐵𝑍 𝑉𝐵𝑍𝑘Plot DOS at each 𝑁𝑘 and 𝜎, then compare with tetrahedron methods(Linear and optimized)
Note:Replacecall libtetrabz_fermieng(ltetra,bvec,nb,nge,eig,ngw,wght,ef,nelec)

toef = 0.0d0(2) Output 𝜀𝒌 and 𝛁𝒌 𝜀𝒌 in the FermiSurfer format (vf.frmsf).
Modified code should be submitted as a diff file generated by$ cd ~/libtetrabz/$ git add example$ git commit$ git show > 12212018.diff24/24

Today's summary

 Total energy Hartree, atomic, XC, kinetic, Ewald energy
 Coulomb repulsion of periodically aligned ions is computed
efficiently with Ewald's method.
 Brillouin zone integration
 Tetrahedron method
 Optimized tetrahedron method
 Library Fermi surface