1/26

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

Lecturer: Mitsuaki Kawamura (河村光晶)

2/26

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/26

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/26

Today’s Schedule
Kohn-Sham eq. in periodic system
Plane-wave representation
Iterative diagonalization
Matrix-vector product
Kinetic energy term
Potential term
Visualize

5/26

last week

Kohn-Sham method

𝐸 𝜌 = න 𝑑3 𝑟𝜌 𝒓 𝑣 𝒓 + 𝐸𝑢𝑛𝑖𝑣 𝜌𝑁atom

𝑍𝐼𝑣 𝒓 ≡ ෍𝒓 − 𝑹𝐼
𝐼=1′
1𝜌(𝒓)𝜌(𝒓)33′𝐸𝑢𝑛𝑖𝑣 𝜌 = 𝑇𝐾𝑆 𝜌 + ඵ 𝑑 𝑟𝑑 𝑟+ 𝐸𝑋𝐶 𝜌′
2𝒓−𝒓

Kinetic energy of non-interacting system whose charge density is 𝜌(𝒓)
𝛁2−+ 𝑣𝐾𝑆 [𝜌] 𝒓2

𝑁/2𝜑𝑛 𝒓 = 𝜀𝑛 𝜑𝑛 𝒓

𝜌 𝒓 = 2 ෍ 𝜑𝑛 𝒓

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

Self-consistent field (SCF)

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

6/26

last week

Kohn-Sham eq. for periodic system (1)𝛁2−+ 𝑣𝐾𝑆 𝒓
2𝑣𝐾𝑆 𝒓 is periodic,i.e., 𝑣𝐾𝑆 𝒓 + 𝑛1 𝒂1 + 𝑛2 𝒂2 + 𝑛3 𝒂3 = 𝑣𝐾𝑆 𝒓
𝜑 𝒓 = 𝜀𝜑 𝒓

Unit lattice vectors(Not unique)𝑁 electrons per unit cell
𝑁𝑐 cells → ∞𝑁×𝑁𝑐

𝜌 𝒓 = ෍ 𝜑𝑛 𝒓
𝑛=1𝒂12Equation to solve in the whole region of bulk crystal

𝒂1𝒂2Unit cell

7/26

Kohn-Shameq.
forperiodicsystem(2)last week
𝜑(𝒓) can be written as

𝛁+𝒌−21 𝑖𝒌⋅𝒓
𝜑 𝒓 =𝑒 𝑢𝑛𝒌 (𝒓)𝑁𝐶න 𝑑 3 𝑟𝜌 𝑟 = 𝑁2+ 𝑉𝐾𝑆 𝒓

𝑢𝑛𝒌 𝒓 = 𝜀𝑢𝑛𝒌 𝒓

Equation to solve only in the unit cell

𝒂1𝒂2𝒂3Bloch’s theorem

𝜀𝑛𝑘𝑈.𝐶.
∞𝜌 𝒓 =21
න 𝑑 3 𝑘 ෍ 𝑢𝑛𝒌 𝒓
𝑉𝐵𝑍 𝐵𝑍2 𝜃(𝜀 − 𝜀 )F𝑛𝑘𝑛=1Band structure

𝒃1 , 𝒃2 , 𝒃3 = 2𝜋𝐼መUnit reciprocal lattice vectors
𝜀𝑛𝑘 and 𝑢𝑛𝒌 𝒓 isperiodic with 𝒃𝛼 in the k space.

𝒃1𝒌

8/26

last week

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

෍𝑮′𝑢𝑛𝒌 𝒓 = 𝜀𝑛𝑘 𝑢𝑛𝒌 𝒓𝑒 𝑖𝑮⋅𝒓

𝑉𝑢𝑐

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

9/26

Unit in this lecture

Input and output :Length : Angstrom
Energy : Electron volt

Useful to read by VESTA

Inside program code:
Length : Atomic unit (Bohr radius). 1 [Bohr] = 0.529177249 [Å]Energy : Atomic unit (Hartree). 1 [Eh] = 27.21138456 [eV][Bohr-1]

෍𝑮′𝑮+𝒌 2
𝛿𝑮𝑮′ + 𝑣෤𝐾𝑆 (𝑮 − 𝑮′ ) 𝑢෤ 𝑛𝒌 𝑮 = 𝜀𝑛𝑘 𝑢෤ 𝑛𝒌 𝑮
2
[Eh]

[Eh]

10/26

Fourier Transformation−𝑖𝑮⋅𝒓𝑒𝑣෤𝐾𝑆 𝑮 = න 𝑑 3 𝑟𝑣𝐾𝑆 (𝒓)𝑉𝑢𝑐𝑢𝑐

Discretize (considering one dimension)
𝑁𝑃𝑊𝑛2𝜋1 −𝑖2𝜋𝑛𝑚×𝐿
𝐿𝑁𝑃𝑊 𝑣 𝐿
𝑁𝑃𝑊 → ∞ : Exact𝑣෤
𝑚 = ෍𝑒𝐿𝑁𝑃𝑊𝑁𝑃𝑊𝑛=1𝑁𝑃𝑊𝑁𝑃𝑊
𝑚=−, ⋯ , −1,0,1, ⋯ ,22
𝑁𝑃𝑊 × 𝑁𝑃𝑊𝑉𝑉෨

2𝑂 𝑁𝑃𝑊numerical cost

Fast Fourier Transformation (FFT) : 𝑂 𝑁𝑃𝑊 ln 𝑁𝑃𝑊 numerical cost
We do not detail in this lecture. We just use FFT numerical library.

11/26

Cutoff frequency
𝑮22 <𝐸𝑐𝑢𝑡𝑢𝑛𝒌 𝒓 =෍𝑢෤ 𝑛𝒌 𝑮𝑮2𝜋𝑢෤
𝑁𝑃𝑊 − 𝑚𝐿𝑁𝑃𝑊𝑒 𝑖𝑮⋅𝒓𝑉𝑢𝑐

1 −𝑖2𝜋𝑛 𝑁𝑁𝑃𝑊−𝑚𝑛2𝜋෨
𝑃𝑊= ෍𝑒𝑉 𝐿=𝑉 − 𝑚𝑁𝑃𝑊𝑁𝑃𝑊𝐿
𝑛=1

12/26

Very simple plane-wave DFT program$ cd ~$ git clone git://git.osdn.net/gitroot/educational-pwdft/pwdft.git$ cd pwdft&CONTROLcalculation = 'direct'$ cp make.inc.mac make.inc/$ make&SYSTEM$ cd sample/Al/
nbnd = 10$ ../../src/pwdft.x < direct.in
nat = 1Aluminum• Face centered cubic (fcc)• Metal• Almost free electron• keep"ecutrho" ≧ 4*ecutwfc

Only perform first SCF step𝑁
Initial density is𝜌 𝑟 =𝑉𝑢𝑐uniform

ntyp = 1
ecutwfc = 60.000000
ecutrho = 240.000000/&ELECTRONS/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
1
0.0 0.0 0.0

13/26

direct@diag_direct.F90

http://mitsuaki1987.github.io/pwdft/index.html
lwork = -1
allocate(work(1))
call zheev('V', 'U', npw, ham, npw, eval_full, work, lwork, rwork, info)
lwork = nint(dble(work(1)))
deallocate(work)
allocate(work(lwork))
call zheev('V', 'U', npw, ham, npw, eval_full, work, lwork, rwork, info)
deallocate(work)

Workspace (memory size) query : Do not need to care.
Diagonalize Hamiltonian matrix "ham"3Computational cost : 𝑂 𝑁𝑃𝑊
3Memory size : 𝑂(𝑁𝑃𝑊)

keep ecutrho ≧ 4*ecutwfc
Changing "ecutrho" and "ecutwfc"
and see "Kohn-Sham Time" in Standard output
Large numerical cost for large 𝑁𝑃𝑊
We only need lower energy (occupied) band.

14/26

Locally Optimal Block Preconditioned
Conjugate Gradient (LOBPCG) method

One of the iterative eigen solution methods
A. V. Knyazev, SIAM J. Sci. Compute. 23, 517 (2001).
山田進 他, 日本計算工学会論文集, 20060027 (2006).

For computing 𝑁𝑏 eigenvector, at each step
Approx. eigenvector 𝒙1 , 𝒙2 , ⋯ 𝒙𝑁𝑏
෡ 𝒙𝑖 − 𝜀𝑖 𝒙𝑖
Residual vector 𝒘𝑖 = 𝐻

Construct subspace Hamiltonian (3𝑁𝑏 dim.)𝒙𝑖†
෡ 𝒙𝑖 , 𝒘𝑖 , {𝒑𝑖 })𝒘† 𝐻(𝑖𝒑†𝑖

CG vector 𝒑1 , 𝒑2 , ⋯ 𝒑𝑁𝑏

Diagonalize with direct method

Take lowest 𝑁𝑏 vectors as approximate eigenvectors at next step
Loop until all𝑟𝑖

become smaller than the threshold

Computational cost : 𝑂 𝛼𝑁𝑏3 + 𝛽𝑁𝑏 𝑁𝑃𝑊 ln 𝑁𝑃𝑊 + 𝛾𝑁𝑏2 𝑁𝑃𝑊Memory : 𝑂 𝑁𝑏 𝑁𝑃𝑊

15/26

Algorithm of LOBPCG method
lobpcg_main@lobpcg.F90

𝒑𝑖 = 𝟎𝑷𝑖 = 𝟎Initial guess 𝒙𝑖 (Random or atomic)
෡ 𝑖𝑿𝑖 = 𝐻𝒙𝜀𝑖 = 𝑿†𝑖 𝒙𝑖
do iteration
All 𝒘𝑖 are small enough ? → Exit
𝒘𝑖 = 𝑿𝑖 − 𝜀𝑖 𝒙𝑖
෠ 𝑖𝒘𝑖 = 𝑃𝒘
𝒘𝑖 = 𝒘𝑖 / 𝒘𝑖෡ 𝑖
𝑾𝑖 = 𝐻𝒘
෡𝑠𝑢𝑏 =𝐻𝑂෠𝑠𝑢𝑏 =𝑾𝑖 , 𝑿𝑖 , 𝑷𝑖†

𝒘𝑖 , 𝒙𝑖 , 𝒑𝑖

𝒘𝑖 , 𝒙𝑖 , 𝒑𝑖†

𝒘𝑖 , 𝒙𝑖 , 𝒑𝑖

෡𝑠𝑢𝑏 𝑥𝑖,𝑠𝑢𝑏 = 𝜀𝑖,𝑠𝑢𝑏 𝑥𝑖,𝑠𝑢𝑏Solve 𝐻

Take lowest 𝑁𝑏 : 𝑥𝑖,𝑠𝑢𝑏 𝑙𝑜𝑤𝒙𝑖=𝒘𝑖 , 𝒙𝑖 , 𝒑𝑖𝑥𝑖,𝑠𝑢𝑏 𝑙𝑜𝑤𝒑𝑖=𝒘𝑖 , 𝟎 , 𝒑𝑖𝑥𝑖,𝑠𝑢𝑏 𝑙𝑜𝑤𝑿𝑖=𝑾𝑖 , 𝑿𝑖 , 𝑷𝑖𝑥𝑖,𝑠𝑢𝑏 𝑙𝑜𝑤

𝑷𝑖=𝑾𝑖 , 𝟎 , 𝑷𝑖𝑥𝑖,𝑠𝑢𝑏 𝑙𝑜𝑤

𝒙𝑖 = 𝒙𝑖 / 𝒙𝑖
end do iteration

𝑿𝑖 = 𝑿𝑖 / 𝑿𝑖

𝒑𝑖 = 𝒑𝑖 / 𝒑𝑖

𝑷𝑖 = 𝑷𝑖 / 𝑷𝑖

16/26

Hamiltonian vector producth_psi@hamiltonian.F90

෡ 𝑖
𝑾𝑖 = 𝐻𝒘
𝑮+𝒌 2
෡ 𝑢෤ 𝑮 = −H+ 𝑣ො𝐾𝑆 𝑢෤ 𝑮2𝑣ො𝐾𝑆 𝑢෤ 𝑮 = න 𝑑3 𝑟𝑒 −𝑖𝑮⋅𝒓𝑢𝑐

𝑢𝑛𝒌 𝒓 = ෍ 𝑢෤ 𝑛𝒌 𝑮
𝑮𝑒 𝑖𝑮⋅𝒓

𝑉𝑢𝑐𝑣𝐾𝑆 𝒓 𝑢 𝒓

𝑉𝑢𝑐 −𝑖𝑮⋅𝒓≈෍𝑒𝑣𝐾𝑆 𝒓 𝑢 𝒓
𝑁𝒓𝒓

𝑉𝑢𝑐Done in subroutine hpsi in hamiltonian.F90

17/26

Running iterative method$ cd ~/pwdft/sample/Al/$ ../../src/pwdft.x < iterative.in

See "Kohn-Sham Time" inStandard output
keep ecutrho ≧ 4*ecutwfc&CONTROL
calculation = 'iterative'/&SYSTEM
nbnd = 10
nat = 1ntyp = 1
ecutwfc = 60.000000
ecutrho = 240.000000/&ELECTRONS/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
10.0 0.0 0.0

18/26

Report problem 1(1) Compare the computational time of direct method andLOBPCG method by changing ecutrho and ecutwfc.

Time [sec]keep ecutrho ≧ 4*ecutwfc

ecutwfc [Ry](2) Plot the cutoff-energy dependence of Kohn-Sham energy 𝜀𝑛𝑘 for
𝑘 = (0,0,0) and 𝑛 = 1, 2 (they are outputted to a file band.dat).

19/26

What is the command "git"$ git clone ssh://user-name@133.11.72.58:/home/Student/Public/pwdft

•••
•••

Backup code
Port program in any system
Merge modifications
Compare diffLog each modification with commentEtc.

20/26

Considered case 1

System C
Laptop PC

PC cluster

Labo PC
System B

Home PC・ Where is the latest source code ?
・ We need to update source code easily.

21/26

Considered case 2
Program

System A

System B

Speed up with symmetry
Program

System A

OpenMP

Program

System A

OpenMP+MPI

Program

System A

System B

・ We need to trace when was it broken.
・ We need to be available to return the stable (no error) code.

22/26

Version control
Program (commit 1)

Repository
server$ git pull
update modification

Program (commit 2)
Program (commit 3)⋮Program (commit 10) Latest
Program (commit 11)$ git commit$ git push

$ git pull
Update latest version
local repository

Logging :Who
WhenWhatmodified

local repository
Modify23/26

When modified simultaneously
Program (commit 1)
Repository
server$ git push

Program (commit 2)
Program(commit 3)⋮
Program (commit 10)
Program (commit 11) Latest

commit 10+$ git pull

Program (commit 12)

commit 11+
local repo
commit$ git pushlocal repo

24/26

Tutorial$ cd ~/pwdft/$ git branch -a$ git config --global user.name "user-name"$ git config --global user.email "user-name@bkks"$ git checkout -b user-name$ echo "Hello, I am user-name." > user-name.txt$ git add user-name.txt$ git commit (vi is used)$ git push --set-upstream origin user-name$ git branch -a25/26

Example of branches
MasterV0.1

V0.2

Hotfix (Bug fix)
ReleaseDevelopFeature (new features)"A successful branch in git"

V1.0

26/26

Today's summary Plane-wave representation of Kohn-Sham eq.
 Direct method and LOBPCG method Hamiltonian-vector product with FFT Git