第３１回コンピュテーショナル・マテリアルズ・デザイン(CMD)ワークショップ
先端事例講義2017年9⽉11-15⽇、⼤阪⼤学基礎⼯学研究科

電⼦ダイナミクス第⼀原理計算と光科学への応⽤
First-Principles calculation for electron dynamics
and its applications to optical sciences
筑波⼤学計算科学研究センター
⽮花⼀浩
Center for Computational Sciences, University of Tsukuba
Kazuhiro Yabanar2

c 2 @t2 A =4⇡J/c

Bulk material(Electromagnetics)

i~@t

=Ĥ

Crystalline unit cell
(Quantum mechanics)

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACS

SALMON Project(Univ. Tsukuba + Institute for Molecular Sciences)

Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience
Software for first-principles real-time electron dynamics calculation-

Numerical experiment platform simulating forefront optical experiments
Software usable by experimentalists and company researchers

・ポスト「京」重点課題７「次世代の産業を⽀える新機能デバイス・⾼性能材料の創成」
サブ課題B「光・電⼦融合デバイス」
・JST-CREST 研究課題「光・電⼦融合第⼀原理ソフトウェアの開発と応⽤」

http://salmon-tddft.jp/
Open source, Apache 2.0

SALMON

Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience
http://salmon-tddft.jp/

Density functional theory (DFT)First-principles calculaton for electronic ground states"%' ! 2 "" "" n(r')" )' ""
2
2∇ + ∑Vion (r − Ra ) + e ∫ dr' " " + µ xc (n(r ))* ϕ i (r ) = εiϕ i (r )&−r − r'
'( 2m
'+
a
!
! 2
n(r ) = å fi (r )iTime-dependent density functional theory (TDDFT)
Ab-initio description for electron dynamics!
ìï " 2 ! 2!
! !
! n( r ' , t )!
! üï !
¶
2i" y i ( r , t ) = í Ñ + å Vion (r - Ra ) + e ò dr ' ! ! + µ xc (n(r , t )) + Vext (r , t )ýy i (r , t )¶tr - r'
ïî 2m
ïþ
a
!
! 2
n( r , t ) = å y i ( r , t )iReal-time, real-space computational method
Time-dependent Kohn-Sham equation

!
!
!
!
! 2
¶
y i (r , t ) = hKS [n (r , t )]y i (r , t )n (r , t ) = å y i (r , t )
¶t
i
2
2
!
!
! e
!
!
"
hKS [n (r , t )] = Ñ2 + Vion (r ) + ò dr ' ! ! + µ xc [n (r , t )] + Vext (r , t )2mr - r'

i"

Norm-conserving pseudopotential
High-order finite-difference scheme using 3D Cartesian grid

Algorithm for time evolution: Explicit method using 4th order Taylor expansionN
1 æ hKS (t )Dt öé hKS (t )Dt ù()
y i (t + Dt ) = expêyt»ç
÷y i (t ),
å
i
ú
i! ø
ë i! û
k = 0 k!èN =4

SALMON

Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience
http://salmon-tddft.jp/

Calculations feasible by SALMON: Interaction between pulsed light and matters

101

Basic optical quantities:
polarizability, dielectric function,
nonlinear response,…
Various applications:
measurements, communications,energy, lighting, processing,chemistry, medical applications,…

Nonlinear optics:
saturable absorption

Excitation energy (eV/atom)

10-0
10-1

Saturable
regionDamageregion

10-2
10-3
10-4
10-5
10-6
109

Linearregion
1010

1011

1012

1013

Intensity (W/cm2)

Attosecond science:measuring electron
dynamics in time
domainE. Goulielmakis et.al,
Nature 466, 739 (2010)

Femtosecond nonthermal processing

B.N. Chichkov et.al, Appl. Phys. A63, 109 (1996)

1014

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACSLinear response calculation in real time:
Response after impulsive force
applied to electrons

Dipole moment!
!
p (t ) = ò dr zr (r , t ) µ a (t )!
Vext (r , t ) µ d (t )ztime

Fourier transf.
Oscillator strength µ w Im a (w )p-p* oscillations electron
emission

Ethylene (C2H4) molecule

Excitation energy

Photoabsorption of molecules

calc.
exp.

K. Yabana, Y. Kawashita, T. Nakatsukasa, J.-I. Iwata, Charged Particle and Photon Interactions with Matter:
Recent Advances, Applications, and Interfaces Chapter 4, Taylor & Francis, 2010.

N2

C6H6

H2HO2O

C60

9

Optical responses of nano-particles (by IMS group)
Surface plasmon of Au nano-particles

4 nm

40 6

Oscillator strength (1 / eV)

Au1414Au560

35 4 146
30 2

1414
54

25 02.0 2.5 3.0
2.0

922

3.0

20
560

15

308

10
5
0
1.5

2.0

2.5

3.0

3.5

4.0

Energy (eV)15,554 electrons

0.08

strength (1 / eV)

Au54

0.07 0.02
largest simulation
in 0.2real nanostructures
0.06

0.010.12.1
K. Iida, M. Noda, K. Ishimura, K. Nobusada,
JPCA,2.5118, 11317 (2014)
2.3

0.05

2.1

0.04

54
146

2.3

2.5

10

At present, not available in SALMON!
Vext (r , t ) µ d (t )z
Hit a molecule with chirality,
and measure the angular momentum
gives us Circular Dichroism (CD) spectrum
K. Yabana, G.F. Bertsch, Phys. Rev. A60, 1271(1999)

Hit a molecule in a magnetic field,
and measure the dipole moment
gives us Magnetic Circular Dichroism (MCD)
K.-M. Lee, K. Yabana, G.F. Bertsch,
J. Chem. Phys. 134, 144106 (2011)

MCD of C60

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACS

Electron dynamics in solids: theory and computational method
Crystalline SiliconTime-dependent band calculation

First-principles band calculation
For electronic ground state(Density functional theory = DFT)(Time-dependent density functional theory = TDDFT)⎡ 1 ⎛ ! ! e ! ⎞2
!
! ⎤
!
∂
i! unk ( r, t ) = ⎢⎜ p + k − A (t ) ⎟ +V ( r, t )⎥ unk ( r, t )⎠∂tc⎢⎣ 2m ⎝
⎥⎦

" " 2
⎧ 1
!⎫
!
"
−i!∇ + k +V ( r )⎬ unk! (r ) = ε nk! unk" (r )
⎨
⎩ 2m
⎭

(

)

t
!
!
A (t ) = −c ∫ dt 'E (t ')Electron density change by pulsed light!ω = 1.55eV
I = 1.0 ×1014 W/cm 2

Dielectric function calculation in real time
(Si, LDA)
Apply an impulsive field to all electrons
E (t ) = kd (t ),A(t ) µ q (t )Induced current = conductivity in time

Fourier transformation → dielectric function
14pis (w )s (w ) = ò dt eiwt J (t ) e (w ) = 1+kw

J (t ) = ks (t )t » 2fs

It is possible to calculate c(2), c(3)14Using improved energy functional, it is possible to get rid of ‘band gap problem’
S.A. Sato et.al, J. Chem. Phys. 143, 224116 (2015)

meta-GGA (Tran and Blaha)Hybrid functional, HSE03

F. Tran, P. Blaha, PRL102, 226401 (2009)

HSE03: J. Heyd, G.E. Scuseria, M. Ernzerhof,
J. Chem. Phys. 118, 8207 (2003)

Imaginary part of dielectric function of Si
TB-mBJ

LDA
Exp
15

Real-time calculation for nonlinear polarizability: χ(3)From calculations using electric field
Of two different intensities, one may
Extract nonlinear polarization.

For a given electric field, we can calculate
polarizability in a unit cell of solid.

⌘2
@
1 ⇣
e
ub,k (r, t) =i~r + ~k + A(t) + vion + vH + vXC ub,k (r, t)
@t
2m
c
Z
⇣
⌘2
e X1
e
J(t) =
Re u?b,k (r, t)
i~r + ~k + A(t) ub,k (r, t) dr + JNL
m
⌦ ⌦
c

i~

P(t) =p(1) (t)E + p(2) (t)E 2 + p(3) (t)E 3 + · · ·p(1) (t)b,kZ tJ(t0 ) dt0

0.15
pz(1)(t) (au)P(t) =

a-quartz
SiO2

I 1, I 3
I 2, I 3
E z(t)

0

-0.15
0

5

1.0
0.5
0.0
-0.5
-1.0
5

10
Time (fs)15

20

Polarization (C/m2)

Polarizability

15

20

5

0

-6
0

5

10
Time t (fs)𝑝(#) (t) =𝜒 (#) 𝐸(t)3

0.010
0.005
0.000
-0.005-0.010020

6Electric field015

p(3) (t)pz(3)(t) (au)Electric field (GV/m)

1.55eV (800nm)

10Time t (fs)10Time (fs)15

20

𝜒 (#) Value, Exp. 2.5x10-22 V2/m2, Cal. 4.3x10-22 V2/m2(LDA)

Saturable absorption of graphite(Collaborator: S. Kurata, N. Kawaguchi (IHI))

k-space- At the beginning of the pulse, most valence electrons are excited.
Then remaining part of the pulse are not absorbed.
- Saturable absorption of graphene is applied to Q-switch of mode-locked laser.

Energy transfer from a pulsed light to electrons in graphite

101

Excitation energy (eV/atom)

10-0
10-1

Damage region• I < 109 W/cm2
Energy transfer proportional topulse duration (linear regime)

Saturation region

10-2
10

• I > 1010 W/cm2
Energy transfer becomes
independent of the pulse
duration (saturated regime)

-3

10-4
10-5
10-6
109

Linear region

1010

1011

1012

1013

1014

Intensity (W/cm2)
18

TDDFT calculation
Electron current induced by pulsed light of 109 – 1011 W/cm2

Applied electric field (blue) vs induced current (green)

Atomic unit (10-4 au)

I0=109 W/cm2

I0=1011 W/cm2

I0=1010 W/cm2

Time (fsec)

Time (fsec)

I < 109 W/cm2 Ohmic current (proportional between current and field)
I = 1010 W/cm2 Phase change during irradiation
I = 1011 W/cm2 Metal to insulator immediately after irradiation

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACS

Interactions between pulsed light and matters𝐼)*+, : 3.51×1045 W/cm8

Maximum power

Solar constant

0.1366 W/cm 2
LinearNonlinear optics

1013 - 1015 W/cm2
Rescattering, Laser damage
Nonthemal laser processing

10 22 W/cm 2
Laser acceleration
Vaccum breakdown𝐼)*+, : 2.32×108: W/cm8- Perturbative -e (w )c (2 ) (w1 , w2 ), c (3) (w1 , w2 , w3 )- Extreme nonlinear response High harmonic generation
Nonthermal laser processing

Ultrashort laser pulse and various ultrafast phenomena

Picosecond

10

-12

s

Optical phonon(Si, 64 fs)

Femtosecond

Attosecond
-18

-15

10 s

10 s
Period of electron
In hydrogen atom,
150as
Ti-sapphire
Laser period

Attosecond measurement of laser pulse shape

-23

10 s

Time

Motion of protons
And neutrons in nuclei

Shortest laser pulse
available ~ 50 as

Real-time observation ofvalence electron motion
E. Goulielmakis et.al,
Nature 466, 739 (2010).

22

Electron rescattering and high harmonic generation
Electrons emitted from atoms and moleculesby light electron field are accelerated by the light electric field
and recollide with the atoms and molecules,
to produce extremely nonlinear and high frequency field.

I (w ) µ ò dte d A (t )! !
dD (t ) d
= ò dr zr (r , t )
dt
dt

iwt

2

9

10

0.8
0.6

8

10
0.4

7

10

0.2
0.0

6

10
-0.2

5

10

-0.4
-0.6

4

10
-0.8

0

5

10

15

20
x10

3

25

30

35
3

10

2

10

0

10

20

30

40

50

Ethylene (C2H4 )Irradiating further strong laser pulse
produces highly-charged ion and Coulomb expansion…
N2 molecule
I=3.35x1015W/cm2, 27fs

Irradiation of high intensity pulsed light on Silicon crystal
First-principles band calculation（Density functional theory＝DFT）

Time-dependent band calculation（Time-dependent density functional theory =TDDFT)⎡ 1 ⎛ ! ! e ! ⎞2
!
! ⎤
!
∂
i! unk ( r, t ) = ⎢⎜ p + k − A (t ) ⎟ +V ( r, t )⎥ unk ( r, t )⎠∂tc⎢⎣ 2m ⎝
⎥⎦

" " 2
⎧ 1
!⎫
!
"
−i!∇ + k +V ( r )⎬ unk! (r ) = ε nk! unk" (r )
⎨
⎩ 2m
⎭

(

)

t
!
!
A (t ) = −c ∫ dt 'E (t ')Electron density change by optical electric field!ω = 1.55eV
I = 1.0 ×1014 W/cm 225

Nonthermal laser processing using femtosecond laser pulse

Shorter pulse

Longer pulse

B.N. Chickov et.al, Appl. Phys. A, Mat. Sci. Proc. 63, 109 (1996)‘Crator’ created on a surface of glass by single-shot laser pulse
Depth is about 100nm (10-7m)

B. Chimier et al, Phys. Rev. B84, 094104 (2011)26

Initial process of laser processing(Only electron dynamics in first few tens of femtosecond)Time

femtosecond

10-15 spicosecond

10-12 sOur approach
Nonthermal processing
Multiphoton or tunnelingionization(single electron motion)

Thermal processing
Avalanche mechanism(electron-electron collision)
Light intensity» 1014 W/cm 2
1 - 2 J/cm

2

27

Two time scales in laser processing simulation
using femtosecond pulsed light＜100 femtosecond
Energy transfer frompulsed light to electrons
Electron density change after irradiating
10fs pulsed Light on silicon!ω = 1.55eV
I = 1.0 ×1014 W/cm 2＞Picosecond
Energy transfer fromelectrons to atomic motions
and ablation process
Molecular dynamics simulation after irradiating 100fs
pulsed light on metal surface, taken from a web page
Of Gauss Center for Supercomputing (Germany).
Simulation done by Dr. Johannes Roth (Univ. Stuttgart).

Application to initial process of nonthermal laser processingof transparent materials
Crator created on a glass surface by single-shot laser pulse- At which intensity of the laser pulse, glass becomes opaque?
- How deep an intense pulsed light penetrated inside the material?

Multiscale, large-scale simulation coupling
29
Maxwell equation (light propagation) and real-time TDDFT (electron dynamics)

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACSA pulsed light irradiating normally on silicon surface

Pulsed light
λ = 800nm,!ω = 1.55eV(below direct gap)

A/c

Vacuum
e=1

Si
e=16

Dielectric constant
!# 16
(Si)
ε ( z) = "#$ 1 (vacuum)
Wave equatione (z ) ¶ 2
c2¶2
A( z , t ) - 2 A( z , t ) = 0
2¶t¶z

For strong pulsed light,
use real-time TDDFT instead of dielectric function.

[mm]

31

Light propagation description: Coupling with Maxwell equations‘Multi-scale’ coupling
Macroscopic gridfor Maxwell equationµm scale

Microscopic gridfor RT-TDDFT electron dynamics

ｚ

1 ¶2
¶2
4p
(
)
(
)
A
t
A
t
=
J Z (t )ZZ
c 2 ¶t 2¶Z 2
cless than nm scale

!
!
!
!
∂
⎡
i! uZ,nk! ( r, t ) = h ⎣ AZ (t ), nZ ( r, t )⎤⎦ uZ,nk! ( r, t )
∂t
!
1 ! !
J Z (t ) = ∫ jZ ( r, t )V V

At each macroscopic grid point,
we solve real-time electron dynamics in parallel
K. Yabana et.al, Phys. Rev. B85, 045134 (2012).

Propagation of weak pulse(Linear response regime, ordinary linear optics applies)
I=1010W/cm2
A/c
Laser frequency：1.55eV：
lower than direct bandgap
2.4eV(LDA）

VacuumSi

200 grid points in Si region.
Representative 3 points are shown.
[µm]

Z=0 [µm]

Z=0.8 [µm]Z=1.6 [µm]33At which intensity a glass thin film starts to absorb light?

SiO2thin film
10µm

Maxwell + TDDFT multiscale simulation : 10 mm SiO2
Laser electric field, red (strong), blue (weak)!ω = 1.55eV

λ =800nm
I=7 ×1013 W/cm 2

Energy transfer from laser pulse to electrons [eV/atom]
We may calculateenergy transferfrom laser pulse
to electrons.

80,000 cores, 20 hours
at K-Computer, Kobe

35

Energy deposited to electrons in SiO2
increases rapidly at a certain laser intensity.

0.1 eV/atom

I=5 ×1013 W/cm 2>> Damage threshold
0

Multiphoton/Tunneling excitation

0.1 eV/atom

I=7 ×1013 W/cm 2

Conduction bandbandgap
10 eV

10mm

1.55eV*7
0.1 eV/atom

I=1×1014 W/cm 2

Valence band36

Intensity at which the glass becomes opaque:
A comparison with experiments
Energy transfer from laser pulse to SiO2 at mid point (5µm).
Calc.

Laser intensity
×1014W/cm2

Exp.

Laser intensity
×1014W/cm2

A. Sommer et.al, Nature 534, 86 (2016)
Exp: Attosecond streaking method
By Max Planck Inst. For Quantum Optics

37

Evalutate the depth of crator by laser pulse

a-quartz (SiO2)!w = 1.55eV (l = 800nm), TFWHM = 7fs
E gap = 9.0 eV (mBJ)

Estimate ablation depth
from a comparison
with cohesive energy[eV/atom]

Ablation depth
100 – 150 nm
Melthing depth
600 - 800 nm

Maxwell + TDDFT calculation
Energy transfer from pulsed light to electrons in SiO2

Ecohesion = 6 eV

Laser ablation depth: Comparison with experiment

Crater formation: measurement, B. Chimier et al, Phys. Rev. B84, 094104 (2011)

Maxwell + TDDFT calculation
provides reasonable estimation
for the ablation threshold
and depth.

S.A. Sato et.al, Phys Rev. B92, 205413 (2015)39

CONTENTS

・About SALMON
・Real-time calculation for molecular polarizability
・Electron dynamics calculations in solids
・A frontier in optical sciences: Intense and ultrashort light・First-principles multiscale simulation for light propagation
・A computational challenge: full nodes calculation at Oakforest-PACS1992 計算物理学研究センター発⾜2004 計算科学研究センターに改組拡充

コ・デザイン：
計算科学と計算機科学の
密接な協⼒
41

SALMON Project(Univ. Tsukuba + Institute for Molecular Sciences)

Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience
Software for first-principles real-time electron dynamics calculation-

Numerical experiment platform simulating forefront optical experiments
Software usable by experimentalists and company researchers

・ポスト「京」重点課題７「次世代の産業を⽀える新機能デバイス・⾼性能材料の創成」
サブ課題B「光・電⼦融合デバイス」
・JST-CREST 研究課題「光・電⼦融合第⼀原理ソフトウェアの開発と応⽤」

http://salmon-tddft.jp/
Open source, Apache 2.0

SALMON is optimized for various processors (Prof. Boku, Mr. Hirokawa, CCS, U. Tsukuba)

Max: 3Tflops

Max: 128Gflops

43

Oakforest-PACS: All-node calculation
U. Tokyo + U. Tsukuba
1st in Japan
7th in the world (2017.6 Top 500 list)
Machine time 2017.3.25-28 (72 hours)

CPUIntel Xeon Phi 7250 with 68 cores, 1.4 GHz base clock# of nodes

8208 (use up to 8192)

Memory16 GB MCDRAM and 96 GB DDR4-2400

Interconnect

Intel Omni-Path Architecture with 100Gbps link

Network topology

Full bisection bandwidth of Fat-Tree

Cooling

Water cooling for CPU and Air cooling for others

OS

CentOS7 (for compute node)Compiler and MPIIntel compiler 17.0.1 and Intel MPI 2017 update 1

Example: pulsed light (1.55eV, 5fs) irradiate on silicon sphere (530nm)

Silicon sphere is expressed as 32,752 points.
At each point, ab-initio electron dynamics
calculations are carried out.
530nm

45

3D Maxwell + 3D real-time TDDFT simulation(a) Macroscopic system(b) Microscopic system (Silicon)

Picture by M. Uemoto46

Laser pulse (1012W/cm2, 1.55eV, 5fs) irradiated on Si nano-sphere2.5 hours by Oakforest-PACS, 8188 nodes

Energy density ofElectromagnetic
fields 𝐸2 + 𝐵2

Electronic excitationenergy by multi-phtonabsorption

47

We expect to build numerical experimental platform for nano-optics

Weak scaling
4 PFLOPS
Hamiltonian performance
4096

450

2048

300
250
200
150
100
50

Graphite

1024
512

Higher is
BetterPerformance [TFLOPS]

350

Lower isFasterDynamics time / Iteration [msec]

400

256
128
64

Silicon0Graphite

Silicon

32
128

512
2048
# of compute node

8192

128

512
2048
# of compute node8192

Summary
Theory and Method:- Real-time, real-space ab-initio simulation for electron dynamics- Interaction of pulsed light with atoms, molecules, nano-materials, and solids- Both linear (polarizability, dielectric function ) and nonlinear regime- Couple electron dynamics calculation with Maxwell equation for light propagation-> enables us to make ‘numerical experiments’, faithfully simulating experiments
Applications:- femtosecond and attosecond solid state sciences- nonthermal laser processing
Code development:- SALMON: Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience- Close collaboration with computer science researchers- Code applicable to various materials and phenomena, andusable for non-experts including experimentalists and company researchers
49

Acknowledgement
Financial supports

Collaborators
Univ. Tsukuba
Mitsuharu Uemoto
Yuta Hirokawa
Taisuke Boku
QST(KPSI)
Tomohito Otobe
Univ. Tokyo
Yasushi Shinohara
Max-Planck Institute for
Structure and Dynamics
of MatterShunsuke Sato

IHI Corporation
Shinichiro Kurata
Norihito Kawaguchi
Univ. Washington
George F. Bertsch
Max Planck Institute forQuantum Optics
Annkatrin Sommer
Martin Schultze
Ferenc Krausz

Supercomputers50