1/23

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

Lecturer: Mitsuaki Kawamura (河村光晶)

2/23

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

3/23
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/23

Today’s Schedule

Kohn-Sham potential
Hartree term
(Pseudo) atomic potential
Exchange-correlation term
SCF loop (Solution of non-linear equation)(Quasi) newton method
Broyden's method
Visualization of grid data
Xcrysden structure format (XSF)5/23DFT loop
Unit lattice vector
atomic position
pseudopotential file𝑁𝜌 𝑟 =𝑉𝑢𝑐

Initial density isuniform𝑣𝐾𝑆 𝜌 𝒓
𝛁+𝒌 2−+ 𝑣𝐾𝑆 𝒓
2𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓∞

𝑁=1
න 𝑑 3 𝑘 ෍ 𝜃(𝜀F − 𝜀𝑛𝑘 )𝑉𝐵𝑍 𝐵𝑍
𝑛=1

∞

1
𝜌 𝒓 =
න 𝑑 3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍

2

End
Output 𝑣𝐾𝑆 𝒓
for post process

𝜃(𝜀F − 𝜀𝑛𝑘 )𝑛=1

No𝑣𝐾𝑆,𝑜𝑢𝑡 𝒓 − 𝑣𝐾𝑆,𝑖𝑛 𝒓

2

is small ?

Yes

6/23

Kohn-Sham potential

𝛁+𝒌 2−+ 𝑣𝐾𝑆 𝒓
2𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓𝑣𝐾𝑆 [𝜌] 𝒓 = 𝑣 𝒓 + 𝑣𝐻 𝒓 + 𝑣𝑋𝐶 𝜌 𝒓
𝑁𝐶 𝑁𝑠 𝑁𝑠,𝑎Atomic (pseudo) potential :𝑣 𝒓 = ෍ ෍ ෍ 𝑣𝑠 𝒓 − 𝝉𝑠,𝑎 − 𝑹𝐿
𝐿=1 𝑠=1 𝑎=1

Hartree potential :′𝜌
𝒓𝑣𝐻 𝒓 = න 𝑑3 𝑟 ′𝒓 − 𝒓′

𝛿𝐸𝑋𝐶 𝜌
Exchange-correlation potential : 𝑣𝑋𝐶 𝜌 𝒓 =
𝛿𝜌 𝒓

7/23

Hartree potential (hartree_pot@rho_v.F90)
𝑁𝐶

𝜌 𝒓′′
𝜌
𝒓3𝑟′𝑣H 𝒓 = න 𝑑3 𝑟 ′=෍න
𝑑′−𝑹𝒓 − 𝒓′𝒓−𝒓𝐿𝑢𝑐
𝐿=1𝑣෤H 𝑮 = න 𝑑3 𝑟 𝑒𝑢𝑐= න 𝑑3 𝑟 𝑒𝑒−𝑖𝑮⋅𝒓𝑉𝑢𝑐𝑢𝑐𝑢𝑐′𝑒 −𝑖𝑮⋅ 𝒓−𝒓 −𝑹𝐿𝑉𝑢𝑐𝑢𝑐= න 𝑑3 𝑟 𝑒𝑣𝐻 𝒓 = න 𝑑3 𝑟 𝑒′−𝑖𝑮⋅𝒓𝑒𝑉𝑢𝑐𝑣H 𝒓 = ෍ 𝑒 𝑖𝑮⋅𝒓 𝑣෤H 𝑮𝐺𝑒−𝑖𝑮⋅𝒓

𝑉𝑢𝑐𝑁𝐶′𝜌𝒓෍ න 𝑑3𝑟 ′′−𝑹𝒓−𝒓𝐿𝑢𝑐

𝐿=1𝑁𝐶′′ 𝑒 −𝑖𝑮⋅𝒓𝜌𝒓෍ න 𝑑3 𝑟 ′′−𝑹𝒓−𝒓
𝐿𝑢𝑐

𝐿=1−𝑖𝑮⋅𝒓𝑒4𝜋′3𝜌 𝒓 න𝑑 𝒓
= 𝜌෤ 𝑮𝒓𝑮2

8/23

Periodic atomic potential(generate_vps@rho_v.F90)
𝑁𝐶 𝑁𝑠 𝑁𝑠,𝑎𝑣 𝒓 = ෍ ෍ ෍ 𝑣𝑠 𝒓 − 𝝉𝑠,𝑎 − 𝑹𝐿
𝐿=1 𝑠=1 𝑎=1𝑣෤ 𝑮 = න 𝑑 3 𝑟 𝑒𝑢𝑐𝑁𝑠 𝑁𝑠,𝑎= ෍෍
𝑠=1 𝑎=1𝑒 −𝑖𝑮⋅𝒓
𝑉𝑢𝑐𝑒 −𝑖𝑮⋅𝝉𝑠,𝑎𝑉𝑢𝑐

𝑁𝑠= ෍ 𝐹𝑠 (𝑮) 𝑣෤𝑠 𝑮𝑣 𝒓 = න 𝑑3𝑟 𝑒𝑢𝑐𝑒 −𝑖𝑮⋅𝒓
𝑉𝑢𝑐𝑁𝐶 𝑁𝑠 𝑁𝑠,𝑎෍ ෍ ෍ 𝑣𝑠 𝒓 − 𝝉𝑠,𝑎 − 𝑹𝐿
𝐿=1 𝑠=1 𝑎=1𝑁𝐶෍ න 𝑑3 𝑟 𝑒 −𝑖𝑮⋅ 𝒓−𝝉𝑠,𝑎−𝑹𝐿 𝑣𝑠 𝒓 − 𝝉𝑠,𝑎 − 𝑹𝐿
𝐿=1 𝑢𝑐𝑁𝑠,𝑎𝑒 −𝑖𝑮⋅𝝉𝑠,𝑎𝐹𝑠 𝑮 ≡ ෍𝑉𝑢𝑐
𝑎=1𝑠=1∞𝑣෤𝑠 𝑮≡ න𝑑3 𝑟 𝑒 −𝑖𝑮⋅𝒓 𝑣𝑠𝒓= 4𝜋 න0𝑣 𝒓 = ෍ 𝑒 𝑖𝑮⋅𝒓 𝑣෤ 𝑮𝐺𝑑𝑟 𝑟 2sin 𝑮 𝑟𝑣𝑠 𝑟
𝑮𝑟9/23

Pseudopotential (generate_vps@rho_v.F90)∞

sin 𝑮 𝑟𝑣𝑠 𝑟 Pseudopotential (𝑟 ≤ 𝑟𝑐 )𝑣෤𝑠 𝑮 = 4𝜋 න 𝑑𝑟 𝑟
𝑣𝑠 𝑟𝑍𝑆𝑣𝑠 𝑟 = ቐ𝑮𝑟0−(𝑟 ≥ 𝑟𝑐 )
𝑟2𝑟𝑐

∞
sin𝑮𝑟
sin 𝑮 𝑟 𝑍𝑆22= 4𝜋 න 𝑑𝑟 𝑟
𝑣𝑠 𝑟 − න 𝑑𝑟 𝑟
𝑮𝑟𝑮𝑟 𝑟0𝑟𝑐

𝑍𝑆 cos 𝑮 𝑟𝑐
−𝑮2
Conpensation of 𝑣෤ 𝐺 = 0 and 𝑣෤H 𝐺 = 0𝑁𝑁𝑠,𝑎lim 𝑣෤𝑠 𝑮 + 𝑣෤H 𝑮
𝑉𝑢𝑐𝑮→0
𝑉𝑢𝑐𝑁𝑠𝑟𝑐
𝑍𝑆 cos 𝑮 𝑟𝑐
2
= 4𝜋 lim ෍ 𝐹𝑠 𝟎 න 𝑑𝑟 𝑟 𝑣𝑠 𝑟 −𝑮→0𝑮20𝑠=1

𝑁𝑠= 4𝜋 ෍ 𝐹𝑠 𝟎𝑠=1𝑟𝑐2𝑍𝑟
𝑆 𝑐න 𝑑𝑟 𝑟 2 𝑣𝑠 𝑟 +20Al (Zion=3)
𝜌෤ 0+𝑮2

10/23

Exchange-correlation potential in LDA(xc_pot@rho_v.F90)

𝐿𝐷𝐴𝑣𝑋𝐶𝜌 𝒓= 𝜀𝑋𝐶 𝜌 𝒓𝑑𝜀𝑋𝐶 𝜌+𝜌 𝒓
ቚ
= 𝜀𝑋𝐶 𝜌 𝒓𝑑𝜌
𝜌=𝜌 𝒓

𝜀𝑋𝐶 𝜌 = 𝜀𝑋 𝜌 + 𝜀𝐶 𝜌𝑟𝑠 𝑑𝜀𝑋𝐶 𝜌−3 𝑑𝑟𝑠

4𝜋 3 1𝑟𝑠 =3𝜌1/3

3 9𝜋1𝜀𝑋 𝜌 = −4𝜋 4𝑟𝑠𝜀𝐶 𝜌 = −0.0480 + 0.031 ln 𝑟𝑠 − 0.0116 𝑟𝑠 + 0.0020𝑟𝑠 ln 𝑟𝑠 (𝑟𝑠 ≤ 1)
0.1423=−(𝑟𝑠 ≥ 1)
1 + 1.0529 𝑟𝑠 + 0.3334𝑟𝑠
J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).

11/23Update 𝑣𝐾𝑆 (𝒓)𝛁+𝒌 2−+ 𝑣𝐾𝑆,𝑖𝑛 𝒓
2𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓∞

1𝑁=න 𝑑3 𝑘 ෍ 𝜃(𝜀F − 𝜀𝑛𝑘 )𝑉𝐵𝑍 𝐵𝑍𝑣𝐾𝑆,𝑖𝑛 𝑟

𝑛=1

∞

1
𝜌 𝒓 =
න 𝑑 3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍2 𝜃(𝜀 − 𝜀 )F𝑛𝑘𝑛=1𝑣𝐾𝑆 [𝜌] 𝒓Plane mixing

𝑣𝐾𝑆,𝑛𝑒𝑤 𝑟 = 𝑣𝐾𝑆,𝑖𝑛 𝑟 + 𝛽 𝑣𝐾𝑆,𝑜𝑢𝑡 𝑟 − 𝑣𝐾𝑆,𝑖𝑛 𝑟0<𝛽≤1𝑣𝐾𝑆,𝑜𝑢𝑡 𝑟12/23

Broyden's method𝐹 𝑣𝑖𝑛 (𝒓) = 𝑣𝑜𝑢𝑡 𝑣𝑖𝑛 (𝒓) − 𝑣𝑖𝑛 (𝒓)መ =𝟎𝑭 𝑣 + 𝐽Δ𝒗

𝑭(𝒗) = 𝒗𝑜𝑢𝑡 (𝒗) − 𝒗(𝑚)Δ𝒗

Solve 𝑭(𝒗) = 𝟎𝜕𝐹𝑖𝐽𝑖𝑗 ≡𝜕𝑣𝑗−1 (𝑚) (𝑚)መ=− 𝐽𝑭𝒗(𝑚+1) = 𝒗(𝑚) + Δ𝒗(𝑚)Newton's method
Approximate 𝐽መ−1 , and update it iteratively → Quasi Newton's method
𝑭(𝑚) − 𝑭 𝑚−1
Minimize−1𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥

𝑚−1𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥−1= 𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥−1= 𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥

𝑖𝑚𝑚−1

𝑚−1−−1𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥𝒗 𝑚 − 𝒗 𝑚−12

C. G. Broyden, Marh. Compur. 19, 577 (1965).

−1𝒗 𝑚+1 − 𝒗 𝑚 − 𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥+𝑚−1

𝑭 𝑚 − 𝑭 𝑚−1𝑭 𝑚 − 𝑭 𝑚−1 ⋅ 𝑭 𝑚 − 𝑭 𝑚−1

𝑭 𝑚 − 𝑭 𝑚−1

T13/23

Practical algorithm (scf_loop@scf.F90)

do iteration : 𝑚 = 1, 2, 3, ⋯Compute 𝑭(𝒗)if 𝐹 2 is small → exit
if 𝑚 > 1
Δ𝑭 = 𝑭 − 𝑭𝑜𝑙𝑑(𝑚)𝓙1

G. P. Srivastava, J. Phys. A: Math. Gen. 17, L317 (1984).
G. P. Srivastava J. Phys. A: Math. Gen. 17, 2737 (1984).
𝑚𝑚−1𝐽መ𝑎𝑝𝑝𝑟𝑜𝑥

𝑚−1𝑛(𝑛)= 𝛽Δ𝑭 + Δ𝒗 − ෍ 𝓙2 ⋅ Δ𝑭 𝓙1𝑛=2(𝑚)𝓙2 =(𝑛) 𝑛 𝑇= −𝛽 + ෍ 𝓙1 𝓙2𝑛=2Δ𝑭Δ𝑭 ⋅ Δ𝑭

Plane mixing(𝛽 = 0.1)

end if𝑚𝑛(𝑛)Δ𝒗 = 𝛽𝑭 − ෍ 𝓙2 ⋅ 𝑭 𝓙1
𝑛=2

𝑭𝑜𝑙𝑑 = 𝑭𝒗 = 𝒗 + Δ𝒗
end doW(001) surface
D. Singh, et al, Phys. Rev. B 34, 8391, (198).

14/23

Tutorial$ cd ~/pwdft/$ git checkout master$ git pull$ make clean; make$ cd sample/Al/$ ../../src/pwdft.x < scf.in
Iteration
1
Average LOBPCG steps :20
delta Vks [eV] :2.7895180896787383E-003
Iteration
2
Average LOBPCG steps :8delta Vks [eV] :1.9722767387705393E-003
Iteration
3
Average LOBPCG steps :9delta Vks [eV] :3.3737045076316851E-005
Iteration
4
Average LOBPCG steps :1delta Vks [eV] :2.3594012412709303E-005
Iteration
5
Average LOBPCG steps :2delta Vks [eV] :1.0001898370361127E-006
Converged ! iter = 5&CONTROL
calculation = 'scf'/&SYSTEM
nbnd = 5nat = 1ntyp = 1
ecutwfc =30.000000
ecutrho = 120.000000/&ELECTRONS
mixing_beta = 0.3
conv_thr = 1.000000e-5
electron_maxstep = 100/
CELL_PARAMETERS
0.000000 2.024700 2.024700
2.024700 0.000000 2.024700
2.024700 2.024700 0.000000
ATOMIC_SPECIESAl al.lda.lps
ATOMIC_POSITIONS
Al 0.000000 0.000000 0.000000
K_POINTS
8 8 8

15/23

Band structure plot

One-shot calculation with previously computed 𝑣𝐾𝑆 (𝒓)$ ../../src/pwdft.x < band.in$ gnuplot
gnuplot> set ylabel "Energy from E_F [eV]"gnuplot> unset xlabel
gnuplot> plot [][:40] "band.dat", ¥x**2*0.5*27.21138456-11.514956871425214 tit "k^2/2-E_F", ¥(x-1.42)**2*0.5*27.2-11.5 tit "(k-G)^2/2-E_F"

Band structure of Aluminum→ Free electron+ Brillouin-zone periodicity(folding)+ small potential effect16/23Visualization of grid data

Real-space grid data
• Potential 𝑣𝐾𝑆 𝒓 , 𝑣 𝒓 , 𝑣𝐻 𝒓 , 𝑣𝑋𝐶 (𝒓)• Charge density 𝜌(𝒓)• Kohn-Sham orbital 𝑢𝑛𝑘 (𝒓)• Etc.
Reciprocal-space grid data
• Kohn-Sham energy 𝜀𝑛𝒌• Fermi velocity 𝛁𝒌 𝜀𝑛𝒌• Berry connection 𝑢𝑛𝒌 𝑖𝛁𝒌 𝑢𝑛𝒌• Berry curvature 𝛁𝒌 × 𝑢𝑛𝒌 𝑖𝛁𝒌 𝑢𝑛𝒌• Etc.

• Isosurface• Section

17/23

Real space : XSF

XCrysDen : Crystalline structure viewer (like VESTA)
XCrysden Structure File (XSF) format

http://www.xcrysden.org/

http://www.xcrysden.org/doc/XSF.html

Unit cell + Atomic position (structure) + Volume (grid) data
Unit : Angstrom
The grid data should include data at all edge ofthe unit cell.

We need additional data at these points : copy from original points

18/23XSF format

Unit lattice vectors [Å]CRYSTAL
Number of atoms in unit cell
PRIMVEC
The second number is always "1"0.00000E+00
0.20247E+01
0.20247E+01
0.20247E+01
0.00000E+00
0.20247E+01
0.20247E+01
0.20247E+01
0.00000E+00
PRIMCOORD
Atomic position in [Å]1 1
Not fractional coordinate
Al
0.00000E+00
0.00000E+00
0.00000E+00
BEGIN_BLOCK_DATAGRID_3D
3D_PWSCF
Data grid number. FFT grid + 1
DATAGRID_3D_UNKNOWN(Including all edge of unit cell)
21 21 210.00000E+00
0.00000E+00
0.00000E+00
Grid shift from atomic position (usually 0)0.00000E+00
0.20247E+01
0.20247E+01
0.20247E+010.00000E+000.20247E+01
Cell of grid data (usually the same as u.c.
0.20247E+010.20247E+010.00000E+000.76311E+020.67851E+020.45909E+020.18866E+02-0.39980E+01-0.16986E+02-0.20120E+02-0.17490E+02-0.13890E+02-0.11749E+02-0.11130E+02-0.11749E+02-0.13890E+02-0.17490E+02-0.20120E+02-0.16986E+02-0.39980E+010.18866E+020.45909E+020.67851E+020.76311E+020.67851E+020.52697E+020.28399E+020.40345E+01-0.12896E+02-0.19773E+02-0.19031E+02-0.15398E+02-0.12470E+02:
END_DATAGRID_3D
END_BLOCK_DATAGRID_3D

Volume data19/23

Input/Output grid data (griddata.F90)

module griddata
contains
subroutine read_griddata(filename, gdata)
use gvec, only : g_rhcharacter(*),intent(in) :: filename
real(8),intent(out) :: gdata(g_rh%nft(1),g_rh%nft(2),g_rh%nft(3))
end subroutine read_griddata
subroutine write_griddata(filename, gdata)
use gvec, only : g_rhcharacter(*),intent(in) :: filename
real(8),intent(in) :: gdata(g_rh%nft(1),g_rh%nft(2),g_rh%nft(3))
end subroutine write_griddata
end module griddata

In pwdft.F90,
program pwdft
use griddata, only : read_griddata, write_griddata
if(calculation == "scf") thencall write_griddata("vks.xsf", Vks)elsecall read_griddata("vks.xsf", Vks)
end ifend program pwdft

20/23

Display 𝑣𝐾𝑆 (𝒓) with VESTA$ cd ~/pwdft/sample/Si/$ ../../src/pwdft.x < scf.in$ vesta vks.xsf• Semiconductor• Diamond structure (lattice is FCC)
• Each Si atom has 4 bonds.
"Objects" → "Properties" → "Isosurfaces..."21/23

Cross section"Objects" → "Properties" → "Sections...""Edit" → "Lattice Planes..."22/23

Report problem 2

Modify pwdft.F90 to output the grid data of charge density 𝜌(𝒓)as a file "rho.xsf".
Then plot the charge density of Si.
Explain the shape of 𝜌(𝒓) (where do the electrons locate ?).
Hint : 𝜌(𝒓) is stored in the following variable
module rho_v!
implicit none!
real(8),allocatable :: && Vks(:), & !< (g_rh%nr) Kohn-Sham potential [Htr]& Vps(:), & !< (g_rh%nr) Pseudopotential [Htr]& rho(:) !< (g_rh%nr) Charge density

Modified code should be submitted as a diff file generated by$ cd ~/pwdft/$ git add src$ git commit$ git show > 12142018.diff23/23

Today's summary Kohn-Sham potential
 Hartree potential is computed with FFT
 Charge compensation of Hartree and ionic potential
 Exchange correlation potential
 SCF loop = Solving non-linear equation Broyden's method
 Visualize grid data (Potential, charge density Kohn-Shamorbital) as a XSF data.