文部科学省科学研究費新学術領域「コンピューティクスによる物質デザイン：複合相関と非平衡ダイナミクス」
平成23年度第１回研究会
東京大学工学部6号館セミナー室A/D 平成23年10月7日（金）－10月8日（土）電子状態計算の計算量
N 個の原子に対する密度汎関数計算( )規格直交条件の演算 ∝ O N 3

計算時間

O(N3)
O(N)

百原子を超えるとボトルネック

百原子以上の大規模計算をするためには計算量がNに線型に依存する方法が必要。
線形スケーリング法

O(N)法
NLinear-Scaling O(N ) DFT
Basic idea of the linear-scaling DFT
In the conventional DFT

φm = ∑ψ n Unmn

ψ n Eigenfunctions of H

φm Wannier functions“localized”“extended over the whole system”

If localization constraint is imposed for electron orbitals, eachorbital overlaps only with a finite number of neighboring orbitals,independent of the size of the system.
Linear-scaling method
Mauri and Galli, PRB50, 4316 (1994)
NN(band bandrrr* r* r
33
˜
E[{ri}, {φ n }]= 2 ∑ ∑ ∫ d rφ m (r )H φ n (r ) 2δ nm − ∫ d rφ n (r )φm (r )m =1 n=1)

O(N) 分割統治法（D-&-C法）
Ω=• 物理系 = 小さなドメインの集合体• 物理的性質 = ドメインの性質の線型結合(W. Yang, ‘91)∑α Ωααααρ (r ) = ∑α p (r )∑ n f (εn ) | ψ n (r ) |Domain Ωα• 演算量が N に線型に依存する計算が可能• 並列計算にも適する

2Buffer

計算法の有効性を十分に検証する必要がある。

Real-Space DFT on Hierarchical Grids• Robust convergence, controlled error bounds, &good energy conservation during MD simulation• On parallel computers,
we need an inter-processor communication of density only → scalability
DFT calculation within each domain: Real-space approach• High-order finite difference (Chelikowsky, Troullier & Saad, ‘94)
• Conjugate-gradient minimization with multigrid preconditioning(Bernholc et al., ‘96; Beck, ‘00)
• Double-grid method (Ono & Hirose, ‘99)
Multigrid accelerationof preconditioned
conjugate gradient

Adaptive fine grids forionic pseudopotentials

Boundary Condition
amorphous CdSe

バッファー長さの増加に対して周期境界条件の方が収束が早い
Domain ΩαBuffer

DC-DFT Computations
Repeat the self-consistent loop until convergence:
1 Calculate potential v(r) from electron density ρ(r), including Hartree potential bysolving Poisson eq. using multigrid method: ∇ 2 vH (r) = −4 πρ (r)2 For all domains αa Iteratively improve electronic wave functions {ψnα(r)} using multigridpreconditioned conjugate-gradient method (6-th order finite difference)⎡ ∇2⎤αα αg n (r) = ⎢− v(r) + ε n ⎥ψ n (r)2⎣⎦b Orthonormalize {ψnα(r)} with the Gram-Schmidt method (inner products)

ψαn ← ψαn − ∑ ψαm ψαm ψαnm<n3 Calculate new density: ρ(r ) =∑domain αα()
wavefunction np (r )∑ααf ε n ψ n (r )2

Scalability of DC-DFT Algorithm
alumina (Al2O3) systems
IBM BlueGene/LSGI Altix 3000
AMD Opteron cluster

The computation time ofDC-DFT algorithm scales
linearly with the numberof atoms.

The largest benchmark:19.2 million-atom DC-DFT calculations of Al2O3 on 131,072
BlueGene/L CPUs at LLNLLiquid Rb
MD simulations of liquid Rb at 350 K

432 Rb atoms

The grid lines show the domains in the DC-DFT calculations.
Atoms move smoothly without any effects of boundary between domains.

Liquid Rb: Energy Conservation
Time evolution of the total and potential energies

The total energy is conserved within 0.2 mHartree/psduring MD simulation.

Liquid Rb: Pair distribution function
DC-DFT432 atoms

The red line
Upper chart: DC-DFTLower chart: Conventional DFTg (r)The white circles
Experimental results( Matsuda and Tamura et al., Phys. Rev. Lett. ‘07)

circles (X-ray diffraction):
Tamura and Matsuda '07

Artificial effects due toperiodic boundary conditions

In contrast, the conventional DFT
result on a smaller system shows
discrepancy.

g (r)Conventional DFT
54 atoms

Extended size by DC-DFT :• Smooth profile over the entire
range of r• Excellent agreement withthe experiment resultsr (Å)

GraphenMD simulation at 300 K

512 C atoms

C atoms show an oscillating motion and the graphene sheet is stable.

Graphen: Vabrational density-of-states
DFT
8 atoms
DC-DFT
512 atoms

Density-functional perturbation theory lattice dynamic

Vabrational density-of-states(VDOS) of graphene
calculated fromvelocity autocorrelation
The 512-atom DC-DFT
represents low-frequency
modes that are absent in the8-atom DFT.
The frequency range and peak
positions for the large system areconsistent with VDOS obtained by ahigh-precision theory (densityfunctional perturbation theory
lattice dynamic).

CdSe Nanorods
CdSe nanorods

Tetrapods
P. Alivisatos(Berkeley)

216 Cd + 216 Se atoms

top view

side view

CdSe Nanorods: Energy Convergence
CdSe

Two parameters to control the data locality:
Domain size, Length of buffer layers

CdSe

DomainBuffer

The energy converges rapidly as afunction of these localization
parameters

CdSe Nanorods: Bandgap
Comparison of the IP, EA, and BG calculated by DC-DFT fordifferent system sizes with conventional DFT results

DFTDC-DFTIP

EA
BG

BG = IP − EA

Both results agree well for small systems.
DC-DFT allows us to study larger systems until the value
converges to the asymptotic value.

Thermite Reaction
http://www.ricen.hokkaido-c.ed.jp/221kagaku_t/2202ka-telmitto/tel.html

Thermite reaction: redox reaction
Fe2O3 + 2Al → Al2O3 + 2Fe
Reaction front propagates at cm/s
in conventional thermite materials.
(diffusion-based mechanisms)

A scene of chemical experiments
in a junior high school

J. Appl. Phys. 98, 064903 (‘05)

Thermal ignition of n‐Al/Fe2O3 mixture
A. Prakash, et al.,Chem. Mater. 16, 1466 (‘04)]• In nanocomposite thermites of Al nanoparticles
embedded in iron oxide, the combustion velocity isaccelerated to km/s from cm/s.
• Spallation (not diffusion) mechanism for the fastreactions?
Levitas, et al., J. Appl. Phys. 101, 083524 (‘07)

O(N) DC-DFT MD of Thermite Reaction
To study the initiation of such
Fast thermite reaction,

Al

Fe2O3/Al interface
1152 atoms with PBC( 144Fe2O3 + 432Al )Fe2O3Temperature 2000 KMD simulationΔt = 0.84 fs
12,000 step ( 10 ps )AlOFeAlO(N) DC-DFT MD of Thermite Reaction
AlReactants
Al & iron-oxide

Al2O3

Fe2O3FeAl

Al2O3Products
aluminum-oxide & Fe

Upper reaction
front position
Lower reaction
front position

OFeAl• The first layer of Al-O bonds
is formed rapidly within 1 ps.
• The propagation of thereaction fronts at a constantvelocity of 70 m/s follows therapid Al-O bond formation.

Time evolution of Mulliken Charges
0 ps 0.7 ps

4 ps

8 ps

Metals change their charges depending on surroundings,
while the charges of oxygen are nearly constant.

Bond-Overlap Population Analysis
The sum of bond overlaps with other Fe and Al atoms for one O atom
Reaction time = 0.25 ps

OiOiFeOiAl

As the O atom crosses theinterface, its character rapidlychanges from iron-oxide typeto aluminum-oxide type within
0.3 fs.

Oiα (t ) = ∑ j∈α Oij (t )Oi (t ) = OiFe (t ) + OiAl (t )AlFeOConcerted metal-oxygen
flip mechanism• Reaction time averaged over all O atoms = 1.7 ± 0.9 ps

Nonadiabatic Quantum MD• Electron transitions by the time-dependent density-functional theory &the fewest-switches surface-hopping method (TDDFT-FSSH)

Potential energy

Tully, J. Chem. Phys. 93, 1061 (’90); Craig et al., Phys. Rev. Lett. 95, 163001 (’05)

State I
Photoexcitationtransition
State JTime• Electronic transitions from the current state to another occur stochastically
based on the transition probability obtained by solving TDDFT equations
Ψ(r,t) = ∑ J C J(I ) (t)Φ J (r; R(t)) C J(I ) (0) = δ I,JJ-th adiabatic⎞ excited stated (I )∂(I ) ⎛C J (t) = −∑ C k (t)⎜iω K δ JK + Φ JΦ K ⎟ K-th excitation⎝⎠dt∂tkfrequency

デンドリマー内のエネルギー伝達過程
http://pltop.shocomarec.kumamoto‐u.ac.jp/index‐j.html

光電子とホールの再結合によるフォトルミネッセンス

光励起された電子とホールはコアへ移動
Zinc-porphyrin core

Hyper-branched
benzyl ether-type peripheries

モデル
Benzyl ether-type peripheries

periO
intercore

HNCZn

Zinc-porphyrin core

front view
91 原子(44C+6OZn+4N+36H)

振動子強度の分布
振動子強度
22m
f0n =
( En − E0 ) n d x 02(eh )d x = e∑ xkk

実験結果（吸収スペクトル）を定性的に再現する。

Akai et al.,phys. stat. sol. (c) 3, 3420 (2006)

光励起後の緩和過程のシミュレーション

LUMO

LUMO+1

LUMO+2

LUMO+3

LUMO+4HOMO

HOMO-1HOMO-2

HOMO-3

HOMO-4

HOMO-4 から LUMO+4 へ電子を励起

励起状態の時間発展： 準電子密度の空間分布

Excitedstatet (fs)The photoexcited electron is transferred to the coreby overlapping the orbitals (Dexter mechanism).

Ongoing Project• Nonadiabatic quantum molecular
dynamics simulation of chargetransfer & charge recombination
at a tetracene-C60 interface

Electronic density of states
LUMO

charge transfertetracene(donor)

LUMOEFHOMOcharge
recombination

C60(acceptor)
HOMO

まとめ• 分割統治法を用いた線形スケーリングDFT計算- 液体金属, グラフェン, CdSe ナノロッド- 酸化鉄(Fe2O3)/Al界面でのテルミット反応
現象の再現に成功、更なる検討が必要
課題
PAW法との結合より適切な境界条件の設定
非局所擬ポテンシャルの計算法の改良
適切なドメインサイズ、バッファーサイズの決定• 非断熱過程のシミュレーション-光機能性物質における高速エネルギー伝達過程
光捕集性デンドリマー内ではデクスター機構でエネルギー伝達
課題
分割統治法との結合

Electronic Excitation• Excited electron-hole pairs within the linear-response time-dependent
density functional theory (LR-TDDFT) [Casida, ’95]
δPijσ (t) = δ Φ(t) aˆ +iσ aˆ jσ Φ(t)δV (t) = δv klτ (t) aˆ +kτ aˆ lτ

χ ijσ ,klτ (t − t′) = δPijσ (t) /δv klτ (t′)

electronhole• Excitation energies from the poles of the response function χijτ,klσ(ω)Ωe(I ) = h2ω I2 e(I ) : 2NunoccupiedNoccupied × 2NunoccupiedNoccupied eigenequation
I-th excitation energy
Fermi distribution
Kohn-Sham energyΩ vοσ ,v ′o′σ ′ = δν ,ν ′δο ,ο ′δσ ,σ ′ (ε νσ − ε οσ ) 2 + 2 ( fοσ − fνσ )(ε νσ − ε οσ )K νοσ ,ν ′ο ′σ ′ ( fο ′σ ′ − fν ′σ ′ )(ε ν ′σ ′ − ε ο ′σ ′ )⎛ e2⎞
δ 2 E xc
K νοσ ,ν ′ο ′σ ′ = ∫∫ ψ οσ (r)ψ νσ (r)⎜+⎟ψ ο ′σ ′ ( r ′)ψ ν* ′σ ′ ( r ′)drdr ′⎝ | r − r ′ | δρσ (r)δρσ ′ ( r ′) ⎠*ν• I-th excited stateΦ I (r;R) =∑∑∑ groundo∈{occupied} v∈{unoccupied} σstateε −ε(I ) +aˆνσ aˆοσ Φ 0 (r;R)× νσ οσ eνοσhω Ielectron-hole pairEFΦ0οaˆν+ aˆο Φ 0

Non-self-consistent excited-state force
The I-th excited state

The excitation energy
Diagonal approximation for atomic forces

The non-self-consistent (NSC) force

Mixing scheme for electron density(Φ0 is the ground-state
Slater determinant)Non-self-consistent excited-state force
Excitation energy

Bond length

一電子波動関数の空間分布（基底状態）

LUMO

LUMO+1

LUMO+2

LUMO+3

LUMO+4HOMO

HOMO-1HOMO-2

HOMO-3

HOMO-4エネルギー固有値の時間変化

電子の状態密度電子は基底状態
300 K

アンテナとコアの両方に広がった状態

分子振動がエネルギー輸送に重要であることを示唆する。