PowerPoint プレゼンテーション

社会人再教育プログラム（コース１-ナノマテリアル・ナノデバイスデザイン学）
励起状態MDシミュレーション１
励起状態のシミュレーション手法
宮本良之
NEC基礎･環境研究所
© NEC Corporation 2004
1
材料設計の方針
1. ある特定元素を特定位置に配置→安定性、電子構造
2. 反応経路の特定→活性化障壁の最小パスの探査
3. 効率のよい反応探査→光励起反応
© NEC Corporation 2004
2
普段、理論屋がなにげなくおこなっていること。。。
波動関数|ψ>が求められている場合。
1
ħ∂
電流密度 J= m <ψ| i
|ψ> + i <ψ|[Vnl, r]|ψ>
∂r
ħ
何故だったか？期待値の時間依存方程式
d
J= d t <ψ| r |ψ>=
i
ħ<ψ| [H, r]|ψ>
何故だったか？時間依存Schrödinger方程式
d
iħ dt |ψ>=H| ψ >
しかし、通常は永年方程式H|ψ>=ε|ψ>を解いている。
密度汎関数理論で時間依存計算をして良いのか？
© NEC Corporation 2004
3
Motivation (Nobody is doing!)
Ab-initio approach on non-equilibrium electron dynamics
excited states
ehν
transport
ee
ground states
1. Optical spectroscopy
2. Electron-phonon coupling (adiabatic/nonadiabatic)
3. Photo-induced dynamics (Ultrafast dynamics)
© NEC Corporation 2004
4
Ultrafast dynamics in femtosecond time constant.
An example: photo-isomerization of retinal
|e>
|g>
F. Gai et al. Science
279, 1886 (1998)
Why fast? Absence/lowering of the activation barrier
Diffusion on |g> takes ps,
while on |e> takes fs.
© NEC Corporation 2004
5
Can we do CP-like MD on |e> with the same manner as MD on |g>?
Sometimes, yes.
Mauri,Car [PRL 74, 3166 (1995)]:
Diamond graphitization by optical excitation.
Why yes?
HOMO-LUMO long living excitation and no level alternation
among differently occupied states throughout the simulation.
But these are NOT always the case.
© NEC Corporation 2004
6
密度汎関数理論→全エネルギー最小の原理
∫∫
ρ(r)ρ(r’)
Etot=∑ifi<ψi|-Δ|ψi>+1/2
drdr’+ Exc{ρ}ρ(r)dr
| r-r’|
ρ(r)Z(R)
dr +∑ifi<ψi|Vnl|ψi>+Eion
|r-R|
δEtot
=0 → HKSψi=εiψi 永年方程式
δψi
∫
∫
これを時間発展方程式に拡張するには量子力学での作用最小原
理の考え方に戻る。 L=<ψ|iħ d ψ> - <ψ|H|ψ>
dt
ただしDFTでは<ψ|H|ψ>の換わりにEtotが来る。
d
LDFT＝ =∑ifi <ψ|iħ ψ> - Etot
dt
δLDFT
d ψ =H ψ
=0
→
iħ
δψi
dt i KS i
© NEC Corporation 2004
7
Hohenberg-Kohn :Physical Review 136, B864 (1964)
基底状態において電子状態と電荷密度は1対1対応
Runge-Gross, Physical Review Letters 52, 997 (1984)
時間発展する電子状態は時間発展する電荷密度と1対1対応
正しい交換相関作用を示す汎関数さえ
見つかれば時間発展問題も精度良く計
算できる。
© NEC Corporation 2004
8
First principles molecular dynamics under electronic excitation.
Potential
|e>
|g>
Reaction coordinate
How to do this?
(Pseudopotentials, Plane-waves, LDA)
Time-dependent DFT for dynamics
under electronic excitation.
More details:
Sugino & Miyamoto PRB 59, 2579 (1999) ;
ibid, B66, 89901(E) (2002).
TM
First-Principles Simulation tool for Electron-Ion Dynamics
© NEC Corporation 2004
9
MD under electronic excitation [Sugino, Miyamoto,PRB 59, 2579, (1999)]
t = 0: Promote the electronic occupations to mimic the
excited states. Then perform the static SCF calculation.
t > 0: Solve ψn(t+Δt)=exp{- i ΔtH(t)} ψn(t).
Hellmann-Feynman
theorem works
|e>
Is Matrix of HKohn-Sham diagonal?
Do MD.
Yes
|g>
1.
No need ofmonitoring
level assignment
for a hole and
an excited
electron
2. Automatic
of the nonradiative
decay
(lifetime,
except
at thewithout
beginning.
decay path)
experiences.
No
lifetime,
Observation of the nonradiative decay!
decay path
© NEC Corporation 2004
10
時間依存Schrödinger方程式
d
i dt|ψ>=H|ψ> をどう解く？
両辺を時間積分！ ψ (t0)がわかっている場合
∫
t+dt
ψ(t0+dt)=ψ(t0)+ (-idt)H(t1)ψ(t1)dt1
t0
t1
t+dt
= ψ(t0)+ (-idt) H(t1){ψ(t0)+ (-idt) H(t2)ψ(t2)dt2} dt1
t0
t0
∫
∫
= …….
© NEC Corporation 2004
11
Integral of the time-dependent Schrödinger equation, time
evolution of wave functions is obtained as follows.
Here, ignorance of T means ignoring the following commutation
relation.
t
Potential V(t) depends on the charge density.
So
must be smaller than typical period of plasma oscillation.
© NEC Corporation 2004
12
How to treat exp(-idtH)ψ?
1
1
1
2
3
Taylor expansion: 1 -idtH + (-idtH) + (-idtH) + 4! (-idtH)4 + ….
3!
2!
Numerical instability beyond 100 fs.
Use of eigen vectors and eigen values:
†
U
exp(-idtε1)
exp(-idtε2)
0
..
.
0
exp(-idtεn)
U
Diagonalization at every timestep is needed
© NEC Corporation 2004
13
Split operator method (Trotter formula)
When H = Δ + V
exp (-idtH)≈exp(-idtΔ)exp(idtV)exp( -idt
Δ)
2
2
Accurate in the second order.
1. Needed more higher order of accuracy.
2. Pseudopotentials: H = Δ + ∑τVnl(τ) + V
© NEC Corporation 2004
14
Suzuki-Trotter formula (split operator for e-iΔtH)
M. Suzuki, and T. Yamauchi, J. Math. Phys. 34, 4892 (1993).
Suppose H consists of non-commutable operators as
.
The second order decomposition of exponential of H is
.
The fourth order decomposition is
.
and
Here,
© NEC Corporation 2004
15
.
< Application for the plane-wave method >
Adapted Hamiltonian(Kohn-Sham - Pseudopotentials)
Wavefunction expression is
or
FFT
.
Exponential of each operator of H acts as
,
,
???????? .
© NEC Corporation 2004
16
Here the nonlocal pseudopotentials are adopted as
(separable K-B type) .
Then we know
and
.
So the exponential is
,
which can directly be operated to
© NEC Corporation 2004
.
17
∑Gρ(G)exp(iGr)=∑G,G’ψ*(G)ψ（G‘)exp{-i(G-G’)r}: convolution
Gcut
Gcut
or
Gcut: for Wfs
Full grid is necessary for timeevoluving WFs to avoid
numerical noise.
FFT
exp(-iΔ/2dt)exp(-iV/2dt).....ψ(r)
© NEC Corporation 2004
Gcut: for charge
18
Car-Parrinello-like smoothing technique in G-space
Ekin
For kinetic energy operator
- /2 = max { G2/2, Gcut2/2 }.
For charge density and local potential
(G) → (G)f(G,Gcut)
and
V H xc(G) → V H xc(G) f(G,Gcut).
Gcut
f
G2
Gcut2
f(G,Gcut)= 1/[1 + exp(
© NEC Corporation 2004
G
19
G2-Gcut2
Gcut2/10
)].
Railway curve interpolation scheme for SCF predictorcorrector routine.
Repeat this part until the
self-consistency between
n,k and VHxc is achieved.
t-Δt
t
t+Δt
VH xc(s)=(
s -t-Δt
Δt
.
s -t-Δt
)2 [3VHxc(t)+ΔtVH xc(t)+
Δt
.
s
.
( 2VHxc (t)+ΔtVHxc (t) )]
.
s-t
( 2VHxc (t+Δt) - ΔtVHxc (t+Δt) )]
Δt
+( s - t )2 [ 3VHxc (t+Δt)-ΔtVHxc (t+Δt) Δt
Eq. 14 in Sugino Miyamoto (PRB) was mistyped!
Iri and Fujimoto, Common sense for numerical calculations(Kyoritsu, Tokyo, ‘85), p116 (Japanese)
H. Akima, J. Assoc. Comput. Machinery, Vol. 17, 589 (1970).
© NEC Corporation 2004
20
CPU 1:
Ψ1(t+Δt)=exp{-i/ħΔtHKS(t)} Ψ1(t)
HKS (t+Δt)
|Ψ1(t+Δt)|2
HKS (t+Δt) CPU 2:
Ψ2(t+Δt)=exp{-i/ħΔtHKS(t)} Ψ2(t)
CPU 0:
ρ(t+Δt)=Σ|Ψi(t+Δt)|2
|Ψ2(t+Δt)|2
HKS (t+Δt)=HKS{ρ(t+Δt)}
Simple communication!
|Ψn(t+Δt)|2
HKS (t+Δt)
Parallel computing
© NEC Corporation 2004
CPU n:
Ψn(t+Δt)=exp{-i/ħΔtHKS(t)} Ψn(t)
21
Advantage of the split-operator method
1. Ψn(t+Δt)=exp{-i/ħHKS(t)Δt}Ψn(t) can be parallelized.
(No need of orthonormalization.)
2. No need to identify occupied and unoccupied
bands throughout the simulation, except the
beginning.
Disadvantage of the split-operator method
1. Full grids of FFT should be used for extending
plane waves to express Ψn(t). (Large core memor
is needed.)
2. Many FFT operations are necessary.
Suitable computer system:
Many processors, each of which has large memory size and high
FFT performance. (Vector-parallel computer system.)
© NEC Corporation 2004
22
Applications of the TDDFT-MD scheme
1. Re-activation of Si donor in GaAs by electronic excitation –
Miyamoto, Sugino, and Mochizuki, Appl. Phys. Lett. 75, 2915
(1999).
2. H-desorption from Si(111) by electronic excitation - Miyamoto,
and Sugino, Phys. Rev. B 62, 2039 (2000).
3. Si-H bond weakening around defects in SiO2 by electronic
excitation – Yokozawa and Miyamoto, J. Appl. Phys. 88, 4542
(2000).
4. Br-desorption from Si(111) by photoemission – Miyamoto,
Solid State Comm. 117, 727 (2001). (TDLDA and TDGGA,
Cutoff 8 Ry, t=0.5 a.u.=1.21x10-2 fs)
5. Self-healing of vacancies in nanotubes, Miyamoto, Berber,
Yoon, Rubio, and Tománek, Chem. Phys. Lett. 392, 209 (2004).
6. Induction of localized vibration on SW defects in nanotubes,
Miyamoto, Rubio, Berber, Yoon, and Tománek, Phys. Rev. B69,
121413(R ) (2004).
© NEC Corporation 2004
23
御清聴ありがとうございました
次回12月1日（水）：励起状態MDシミュレーション２
「ナノスケール物質の高速現象（事例）」
© NEC Corporation 2004
24

Download Report