サイエンティフィック・システム研究会HPCフォーラム2017「計算科学の新たな可能性」
汐留シティセンター 2017年8⽉30⽇

量⼦ダイナミクスの第⼀原理計算と光科学フロンティアへの展開
⽮花⼀浩
筑波⼤学計算科学研究センター

r2c 2 @t2 A =4⇡J/c

Bulk material(Electromagnetics)

i~@t=Ĥ

Crystalline unit cell
(Quantum mechanics)

⾃⼰紹介
1982.4 – 1987.3
京都⼤学⼤学院理学研究科物理学第⼆専攻
原⼦核物理学の理論（博⼠学位取得）
1987.4 – 1988.2
京都⼤学基礎物理学研究所（学振特別研究員）
1988.2-1999.7
新潟⼤学理学部助⼿、助教授
1999.8-2004.3
筑波⼤学物理学系助教授、教授
2004.4-現在
筑波⼤学計算科学研究センター教授
・量⼦物性研究部⾨（部⾨⻑）・原⼦核物理研究部⾨

原⼦核物理学

量⼦（波動）ダイナミクス・シミュレーション

光科学

・量⼦（波動）ダイナミクス計算とは・原⼦核物理学からナノサイエンス、そして光科学へ
・光科学フロンティア：強い光・短い光・電磁気学と量⼦⼒学の多階層連結シミュレーション
・Oakforest-PACS全ノード計算の試み

⾼校物理の最後に学ぶこと

物質と光は、ともに粒⼦でもあり波動でもある。

しかし、現実の問題どのような状況で粒⼦になるか波動になるかは、結構複雑。

陽⼦

中性⼦

原⼦核の中の陽⼦と中性⼦＝波動

物質：電⼦＝波動
原⼦＝粒⼦

光電磁場：
（多くの場合）動DFT（波動）・MD（粒⼦）FDTD法（波動）物質の波動計算の難しさ
粒⼦：古典⼒学（ニュートン⽅程式）
波動：量⼦⼒学（シュレディンガー⽅程式）Né N #2!
! ! ù ! !
!
! !
!
¶()
()
(D+Vr+v(rr)
Yr,r,",r,t=i#Yr,r,",rêåå
iiij ú
1 2N1 2N ,t)
2m
¶
ti< j
ë i =1
û

ニュートン⼒学であれば、!
!
d Ri (t )Mi= å Fij
2dtj2

⾃由度の数は、粒⼦数をNとして3N

量⼦⼒学では、座標を100点に分割すると、⾃由度の数は、粒⼦数をNとして! !
!
Y(r1 , r2 ,", rN , t )1003 N

直接計算が可能なのは、⾼々数粒⼦系（２，３，４，…）まで多数の粒⼦を含む物質をどうやって扱うか？

密度汎関数理論
Ψ 𝑟# , 𝑟% , 𝑟& , ⋯ , 𝑟(波動関数!
密度 n(r ) を扱えば良い。

ではなくて、W. Kohn
1998年ノーベル化学賞

実際の計算は、「平均場近似」
それぞれの粒⼦が、平均的な場の中を独⽴に運動する。
Ψ 𝑟# , 𝑟% , 𝑟& , ⋯ , 𝑟( ≈ 𝐴𝜓# 𝑟# 𝜓% 𝑟% ⋯ 𝜓( 𝑟(1003 N 元々の⾃由度

𝜓, 𝑟

⾃由度の数は、N×1003

厳密な計算に⽐べて、圧倒的に計算量は少ないが、⼤胆な近似をすることになる。
「物質科学の第⼀原理計算は、経験的パラメータは含まないが、⼤胆な近似を伴う。」

波には定在波と進⾏波がある。

これまでの、物質科学の計算：
電⼦波動（軌道）関数の定在波計算

𝐻𝜙, 𝑟 = 𝜀, 𝜙, 𝑟（密度汎関数理論）

4

この講演で考える量⼦（波動）ダイナミクス計算：
𝐻𝜓, 𝑟, 𝑡 = 𝑖ℏ 𝜓, 𝑟, 𝑡45
電⼦波動（軌道）関数の進⾏波計算（時間依存密度汎関数理論）

光が照射した物質の中で、電⼦の運動を記述できる。

原子核衝突のシミュレーション（時間依存密度汎関数理論）
私が⼤学院に⼊学 (1982) する数年前

H. Flocard, S.E. Koonin, M.S. Weiss,
Phys. Rev. 17(1978)1682.

Time evolution of proton andneutron orbitals.
!
!
!
¶y i (r , t ) = h[n (r , t )]y i (r , t )¶t!
! 2n( r , t ) = å y i ( r , t )i#i! !
!
!
!
!
y (r1 , r2 ,", rN , t ) = A{y 1 (r1 , t )y 2 (r2 , t )"y N (rN , t )}16O-16O collision
Spatial grid 30x28x16, time step 4x102

当時のスパコン
CRAY-1 (80MFLOPS)

原⼦核物理の最近の話題
原⼦番号１１３番： 新元素合成（ニホニウム）
2015.12.30 IUPAC, IUPAPが、新元素Z=113, 115, 117, 118を認定。
Z=113元素の命名権が理研グループに。2016.11正式名称決定。

超重原⼦核 Z=120 ⽣成の時間依存密度汎関数理論によるシミュレーション
K. Sekizawa, K. Yabana, Phys. Rev. C93, 054616 (2016)

238U + 64Ni

(proton number: 92 + 28 = 120）

全部で238+64=302個の陽⼦と中性⼦（フェルミ粒⼦）に対する量⼦シミュレーション!
!
!
¶y i (r , t ) = h[n(r , t )]y i (r , t )¶t!
! 2n( r , t ) = å y i ( r , t )i"i空間：フェムトメートル(10-15m)
時間：ゼプトセカンド(10-21s)

10・量⼦（波動）ダイナミクス計算とは・原⼦核物理学からナノサイエンス、そして光科学へ
・光科学フロンティア：強い光・短い光・電磁気学と量⼦⼒学の多階層連結シミュレーション
・Oakforest-PACS全ノード計算の試み

量⼦ダイナミクス計算が、原⼦核物理から光科学・物質科学に展開するきっかけ

1980年代：ナノサイエンスの誕⽣

質量分布：どのサイズのナノ粒⼦ができやすいか？

⾦属ナノ粒⼦（ナトリウム）の魔法数

⾦属ナノ粒⼦： 2,8,20,40,58,92,138、...
原⼦核： 2,8,20,28,50,82,126、...

炭素のナノ粒⼦：
フラーレンC60の発⾒

ナノサイエンスは、様々な分野の学際科学として発展
異分野交流が盛んな時代
物性物理学（固体）、量⼦化学（分⼦）、原⼦核物理学（原⼦核）・・・

ナノ粒⼦の光応答の時間依存密度汎関数理論計算(1996-)
K. Yabana, G.F. Bertsch, Phys. Rev. B54, 4484 (1996).

Li147+

光吸収断⾯積[eV-1]

×Exp:Li139+
TDDFT

全ての電⼦に撃⼒を加えた後の電⼦運動
正20⾯体配位を仮定

13計算法：実時間・実空間解法

YXZ

３次元デカルト座標を離散化：⾼次差分法
擬ポテンシャル
J. Chelikowsky et.al, Phys. Rev. B50(1994)11355

RSDFT (2013 Gordon Bell Prize)と同じ計算⽅法

時間発展のアルゴリズム: 陽解法(テイラー展開）
K. Yabana, G.F. Bertsch, Phys. Rev. B54, 4484 (1996).
N1 æ hKS (t )Dt öé hKS (t )Dt ù()y i (t + Dt ) = expêyt
»ç÷y i (t ),åú ii!
k!
i!
ëû
è
ø
k =0

N =4

物質科学では、これまで平⾯波を⽤いた計算が主流。
この頃から並列化に優れた実空間解法が⼤規模計算で⽤いられるようになる。
14

分⼦にパルス光を照射する

エチレン分⼦( C2H4 )電⼦密度の時間変化
電⼦再散乱
⾼次⾼調波発⽣

結晶にパルス光を照射する

シリコン結晶中の電⼦の密度分布

・量⼦（波動）ダイナミクス計算とは・原⼦核物理学からナノサイエンス、そして光科学へ
・光科学フロンティア：強い光・短い光・電磁気学と量⼦⼒学の多階層連結シミュレーション
・Oakforest-PACS全ノード計算の試み

光と物質の相互作⽤を特徴づけるのは・・・ 屈折率と誘電率n= ε

分散関係（振動数依存性）

⾮線形性（強度依存性）n (ω ), ε (ω )n = n0 + n2 E 218

もう少し、詳しく書くと・・・

光と物質の相互作⽤には、２つの物理学が関わっている。

電磁気学

量⼦⼒学

物質中の光の伝播を記述するのは巨視的マクスウェル⽅程式

光電場と電⼦の相互作⽤に摂動論を⽤いると、誘電率の式が得られる。

D = εE19しかし、それだけでは済まないこともある。
例えば、レーザーを⽤いて、「透明」なガラスを「加⼯」することができる。

Im ε = 0Im ε ≠ 0

つまり、強いレーザーを照射すると、ガラスは透明ではなくなる。
これは、通常の電磁気学では記述できない。

20

⾼強度パルス光と物質の相互作⽤

太陽定数0.1366 W/cm 2最強のレーザー
10 22 W/cm 2 光強度

1013 - 1015 W/cm2

⾮線形光学
χ(2), χ(3)

極限⾮線形光学⾼次⾼調波発⽣

原⼦核反応
真空崩壊

レーザー加⼯⾼次の⾮線形光応答で、可視光からX線を作る。

超短パルス光（フェムト秒）を使うと微細加⼯が可能になる。

超短パルス光と、さまざまな現象の時間スケール

ピコ秒

10-12s

光学フォノン(Si, 64 fs)

アト秒

フェムト秒

-18-1510 s

10 s
⽔素原⼦の電⼦軌
道周期, 150as

Ti-sapphire
レーザー周期-2310 s

時間

原⼦核中の陽⼦・中性⼦運動

最も短いレーザーのパルス⻑ 64 as

アト秒科学：光の波形が測定可能に
http://www.attoworld.de (Max Planck Inst. for Quantum Optics)

電⼦が動く様⼦を撮影する試み
E. Goulielmakis et.al,
Nature 466, 739 (2010).

22

未来の「ペタヘルツ光デバイス」実現に向けて
Attosecond nanophotonics
G. Vampa, H. Fattahi, J. Vuckovic, F. Krausz, Nature Photonics 11, 210 (2017)

・量⼦（波動）ダイナミクス計算とは・原⼦核物理学からナノサイエンス、そして光科学へ
・光科学フロンティア：強い光・短い光・電磁気学と量⼦⼒学の多階層連結シミュレーション
・Oakforest-PACS全ノード計算の試み

結晶表面に入射する強いパルス光

パルス光
λ = 800nm,!ω = 1.55eV(below direct gap)

A/c

Vacuum
e=1

Si
e=16

誘電率を仮定!# 16(Si)ε ( z) = "#$ 1 (vacuum)
波動方程式e (z ) ¶ 2c2¶2A( z , t ) - 2 A( z , t ) = 02¶t¶z[µm]

弱い光は、これで良い。
強い光の場合、誘電率を使わない記述が必要となる。
25

電磁気学と量⼦⼒学を繋ぐ多階層シミュレーション

マクスウェル⽅程式（電磁気学）µm scale

時間依存密度汎関数理論（量⼦⼒学）ｚ1 ¶2¶24p()()AtAt=
J Z (t )ZZc 2 ¶t 2
¶Z 2c

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

マクスウェル⽅程式の各格⼦点で、電⼦ダイナミクス計算を並列実⾏。
K. Yabana et.al, Phys. Rev. B85, 045134 (2012).

弱いパルス光がSi結晶に⼊射する場合(線形領域、誘電率による記述が有効)

I=1010W/cm2A/cVacuumSi

振動数：1.55eV
直接バンドギャップ：2.4eV(LDA）[µm]

Z=0 [µm]

Z=0.8 [µm]200点で電⼦ダイナミクスを計算
代表的な３点を⽰す。

Z=1.6 [µm]27

ガラスの薄膜は、どのくらいの光強度から不透明になるのか？

SiO2薄膜10µm

Laser pulse propagation: 10 µm SiO2Laser electric field, red (strong), blue (weak)!ω = 1.55eV

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

Energy transfer from laser pulse to electrons [eV/atom]

京コンピュータ
10,000ノード, 20時間29あるレーザー強度から、光から電⼦へのエネルギー移⾏が急激に増加

0.1 eV/atom

I=5 ×1013 W/cm 2010µm

0.1 eV/atom

多光⼦/トンネル励起
I=7 ×1013 W/cm 2伝導帯バンドギャップ

10eV

７光⼦吸収

1.55eV*7

0.1 eV/atom

I=1×1014 W/cm 2

価電⼦帯
30

レーザーアブレーションの深さ
パルスレーザーをガラス表⾯に１回照射すると、およそ100nmの⽳があく。

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

電磁気学＋量⼦⼒学を多階層連結した第⼀原理シミュレーションは、レーザーの加⼯閾値と、⽳の深さを、おおよそ再現する。

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

31

SALMON プロジェクト （筑波⼤学＋分⼦科学研究所）
Scalable Ab-initio Light-Matter simulator for Optics and Nanoscience- 先端の光科学実験を精密に模倣する数値実験室を実現する。
- 実験家や企業研究者でも使って役に⽴つソフトウェアを整備する。
・ポスト「京」重点課題７「次世代の産業を⽀える新機能デバイス・⾼性能材料の創成」
サブ課題B「光・電⼦融合デバイス」
・JST-CREST 研究課題「光・電⼦融合第⼀原理ソフトウェアの開発と応⽤」
Open-source software, License: Apache 2.0

http://salmon-tddft.jp/
計算機科学者（朴泰祐教授のグループ）との⽇常的な協働（コ・デザイン）
筑波⼤D2 廣川祐太⽒（４年⽬）

・量⼦（波動）ダイナミクス計算とは・原⼦核物理学からナノサイエンス、そして光科学へ
・光科学フロンティア：強い光・短い光・電磁気学と量⼦⼒学の多階層連結シミュレーション
・Oakforest-PACS全ノード計算の試み

Oakforest-PACS 全ノード計算

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

３次元マクスウェル＋３次元電⼦ダイナミクスシミュレーション(a) Macroscopic system(b) Microscopic system (Silicon)

Picture by M. Uemoto35例：シリコン球 (直径530nm)に、1.55eV、5fsのパルス光が⼊射

シリコン球を、32,752 点で表現
各点で、電⼦の量⼦ダイナミクス
第⼀原理計算を実施
530nm

36

レーザーパルスをナノシリコン球に照射(1012W/cm2, 1.55eV, 5fs)
2.5 hours by Oakforest-PACS, 8188 nodes

電磁場のエネルギー密度
𝐸2 + 𝐵2

励起電⼦分布(２光⼦吸収)37

Weak scaling

4 PFLOPS
Hamiltonian performance
4096

450

2048

300
250
2001501024

Higher is
BetterPerformance [TFLOPS]

350

Lower isFaster

Dynamics time / Iteration [msec]

400

512
256
128

100
6450

3201285122048# of compute node

8192

1285122048# of compute node8192

Strong scaling
Lower is Faster

Dynamics time / Iteration [msec]

1000

100105121024

2048# of compute node40968192

考察: ノード間に⽣じる性能ギャップ

Hamiltonian

Current

Normalized elapse time / Iteration

Lower isFasterMisc. computation Communication
1.4
1.2
10.8

0.127
0.095

0.126
0.109

• 差のほとんどがハミルトニアンの計算時間
1. プロセス間通信は⼀切含まれない
2. 全プロセスの問題サイズは⼀定

0.6
0.4

0.764

0.925

0.20Best

Worst※Bestで正規化• Non-algorithmic load-imbalancing
1. プロセッサを2ノードペアで冷却
2. Turbo Boostによる動的クロック調整

まとめ

量⼦（波動）ダイナミクス計算
原⼦核物理学の分野で⽣まれ、ナノサイエンスを経て光科学に展開
光科学フロンティアと計算光科学
強い光・短い光を⽤いた科学・技術が発展
計算科学による階層連結が役に⽴つ
SALMONプロジェクト http://salmon-tddft.jp
量⼦論の第⼀原理計算による光科学の新しいシミュレータを開発
計算科学者と計算機科学者のインハウス連携

41

謝辞
共同研究者
筑波⼤学
植本光治
廣川祐太
朴泰祐
関⻄原研
⼄部智仁
東京⼤学
篠原康
マックスプランク構造
動⼒学研究所佐藤駿丞

研究費
分⼦科学研究所
信定克幸
野⽥真史
ワシントン⼤学
George F. Bertsch
マックスプランク
量⼦光学研究所
Annkatrin Sommer
Martin Schultze
Ferenc Krausz

スーパーコンピュータ

42