シミュレーションの産業利用
新しい原子論的シミュレーション法/アルゴリズムでチャレンジする自動車関連部品の高性能化

対象物理

1-102nm電子構造，原子配置，化学反応

0.1-10μm

材料の微細構造， 高分子系/固体系/流体系毎に，調べたい物理量
相構造，相変化 に応じて，熱統計力学(自由エネルギー)も援用

尾形 修司
名工大院 創成シミュレーション工学専攻

電子状態，イオン/原子の動力学，分子構造（密度汎関数法，分子軌道法，分子動力学法…）
à基礎手法は，ある程度固まっているが…
規模/時間 vs. 精度の相反する需要à混沌とした状況

名工大 オプトバイオテクノロジー研究センター> 10μm

16:15-16:50, 2015年1月22日(木); 名工大 講堂; バイオスーパーコンピューティング名古屋2015

計算法，現況

対象scale流れ，変形，振動，空力

連続体(流体，構造体)方程式の様々な数値解法（差分法，波数空間法，有限要素法, 粒子法…)

à特定分野の企業では，プレゼンツールでなくCAEとして，重要技術と認知されている

本日の内容

本日の内容

１．電子系のための，オーダーN型の実空間グリッドDFTの開発

１．電子系のための，オーダーN型の実空間グリッドDFTの開発

è 固体電解質皮膜-電解液界面でのLiイオン透過 è２次電池

è 固体電解質皮膜-電解液界面でのLiイオン透過 è２次電池

２．稠密固体にも適用できる，ハイブリッド量子古典法の開発
è無機有機界面の熱伝導，物理化学的強度 è放熱材料，塗装

３．高速な，剛体分子動力学アルゴリズムの開発
è氷と水，不凍タンパク質 è○○液 (with hybrid sim.)

２．稠密固体にも適用できる，ハイブリッド量子古典法の開発
è無機有機界面の熱伝導，物理化学的強度 è放熱材料，塗料

３．高速な，剛体分子動力学アルゴリズムの開発
è氷と水，不凍タンパク質 è○○液 (with hybrid sim.)1

実空間グリッドを用いる Real-space-Grid (RG) DFT

電子群の取扱い法
e+e-

Kohn-Sham Hamiltonian for electrons in DFT	

単純化	

Self-consistent solution for eigen-orbitals: φi (r ) …real variables+e-◉ 複数の電子が，互いにクーロン反発しながら，イオンに引き寄せられてイオン周囲に存在する状態 è 多体問題を直接解くのは大変!
! " e 2 ρ (r!) !
! % !
!
!
!
!2−Δφi (r ) + $$ ∫ ! ! dr! + Vion (r )''φi (r ) + VXC (r )φi (r ) = εiφi (r )2m
# r −r'&à Use finite
difference(e.g., 9pt/dir)Classical	  electronic	  Coulomb	  Classical	  ionic	  Coulomb	  Quantum	  electronic	  Ex.-­‐Corr.	  pot.	  

 Occuppied 2ρ (r ) = ∑ 2 φi (r )i ∫ φ (r )φ (r )dr = δijij

φi (rmax ) = 0

Real-space Grid for orbitals & potential field for KS-DFT-Hamiltonian◉ 定まったポテンシャル場に存在する1電子状態なら容易

Grid: 106 ~ 107 points in spherical shape	èただし，ポテンシャル場に，他の電子が作るポテンシャルも自己無撞着な平均場(Self-Consistent Field)として取り入れる

ion	  

◉	◉	
Orbital-data on a single grid-point (r0): 	
103 ~ 104 components as	 {φ1 (r0 ), φ 2 (r0 ), φ3 (r0 ),...}
è In total, 109 ~ 1011 variables to be determined	

Kohn-Shamの密度汎関数理論 (DFT)RGDFT èOrder-N3 �

Divide-&-Conquer type, Real-space Grid DFT (DC-RGDFT) code
Order-N Real-space Grid DFT	Our formulation to realize:	

Comp. Phys. Comm. 183 (2012) 1664 	

Real-space à Finite difference methodfor derivatives	

(1) Density-continuity at domain-boundary(2) Little artificial cutting-effect on orbitals(from ionic, kinetic, exchange, & correlation)(1) Domain decomposition of a target system: can be up to ~ 27	
(2) Spatial decomposition of grid points in a domain: efficient up to 5x5x5 	
(3) Level parallelization of orbitals: efficient up to 5	ρ	Total	  system	

domain-­‐I	Parallel coding of DC-RGDFT
Parallelization in four points	" ρ for r ∈ Iρ tot = # I$ ρ II for r ∈ II

ρI	bond	ion	
real	
real	
real	
buffer	(4) OpenMP parallelization of various do-loops: efficient up to 4 threads
domain-I	r	

domain-II	àscalable up to 50,000 CPU-cores	For domain-I:	

E( ρ tot ) = E( ρ I ) + [E( ρ tot ) − E( ρ I )]vDT (r) = [ ρ I (r) − ρ tot (r)] / α Density-­‐
template	  = [ ρ I (r) − ρ II (r)] / α poten8al	−[δ (TS + E XC ) / δρ ]ρ =ρ

€localext∑
iIKohn-Sham eq. of domain-I:	

n−localext+v+ vemb ]ψi = εiψiKohn-Sham- [− 2 + vxc ( ρ I )+ vDT + vH,I + v
type SCF eq.
2for domain-I:	 ρ I = 2 ψi …up to global (=all domain) Fermi level	
Δlevel-parallel
computation	

vemb (r) = [δ (TS + E XC ) / δρ ]ρ =ρ
Embed-­‐
tot
poten8al	local[− Δ2 + vxc ( ρ I ) + vH ,I + vext

domain-III	domain-IV	

non−local+vext+ vembed-DT ]ψi = εiψi2

現在のMDシミュレーションで扱える対象系の規模

Benchmark tests of DC-RGDFT code
K-computer	 (Grid size=0.55 aB)	Strong-scaling test à very good!

Weak-scaling test à order-N !
NCPU-CORE

0
2500

1600
wall-clock time (sec)wall-clock time (sec)

1400
1200
1000
800

Si

1728

cluster

600
400
200
00

SCF-loop = 10grid size = 0.55 a.u.
100

200

300

400

500

600

N

2000

4000

6000

8000

Coulomb multipoles (far) = 8
Coulomb multipoles (far) = 1◉ 2大 分子動力学(MD)シミュレーション法＊密度汎関数法(DFT)で電子状態計算 è イオンに働く力で第一原理MD＊古典原子間ポテンシャルによる古典MD◉ スパコン活用の為，2大シミュレーション法のコードを高度並列化

2000

è 粒子間の(実効的)相互作用距離が有限であることを利用する
空間分割法による並列化 è Order-N化

1500
1000

Various Si-clusters
SCF-loop = 10grid size = 0.55 a.u.

500
00

10000

20000

30000

◉ 実用的な，対象系規模
è 数千個のCPUコアを用いて，電子系でも，古典原子系でも，1010個の自由度系 まで

N

CPU-CORESi

Scalable up to 50,000 CPU-cores= 27 domains x 125 spatial decom. x 5 orbitals x 4 threads	

First-principles MD: Modeling and settings�

Target system: Li-ion battery�
Li+-ion transport through SEI-electrolyte boundary	

Molecules

Simulation systems … about 2400 atomsSEI è LEFT

à 72 EDC2- + 144 Li+Electrolyteè RIGHT� à 120 EC vs. 112 EC + 8 LiPF6	
◉ At T= 825 & 1100K	
◉ Various settings for surf. mols for stability; No PBC	→Discharge	Charging← Li+Negativeelectrode(graphite)

Electrolyte

PBE-GGA;
Grid size=0.70aB with 0.23aB	

Positive
electrode(metaloxide)

SEI (Solid Electrolyte Interface)◉ Blocks electrons, but transports Li+◉ 101 nm only in depth◉�A major resistance part for Li+-ions�3

First-principles MD: diffusive behavior of Li+-ions�

Change in Li+-ion detaching energy at the boundary�◉ At time=0, insert 36 Li+-ions at the left-end	From separate DFT analyses	

◉ Every 0.3ps, remove Li+-ions at the right-end	◉ Li+-flow by “pressure” of electric field (Li+) and high Li-density at the left	

T=825K without salt

◉ Case of without salt:
Li+-detaching as Li+-EC
Detaching energy = 1.7 eV

T=825K WITH 1.0 mol/l salt (LiPF6)◉ Case of WITH salt (PF6-):
Li+-detaching as Li+-PF6Detaching energy = 0.9 eV

Followed by dissociation of LiPF6
and solvation (溶媒和) by EC’s:
LiPF6 + 4EC à LiEC4 + PF6-Li-passing: 3.5 Li-ions/5ps�

Li-passing: 6.5 Li-ions/5ps�

本日の内容

１．電子系のための，オーダーN型の実空間グリッドDFTの開発

スパコンでの，DFTと古典MDの計算時間の違い
MD = 古典分子動力学
DFT = 密度汎関数法
オーダーN-DFTè 固体電解質皮膜-電解液界面でのLiイオン透過 è２次電池

２．稠密固体にも適用できる，ハイブリッド量子古典法の開発

スパコンの1024 CPU-coreを使った計算è無機有機界面の熱伝導，物理化学的強度 è放熱材料，塗料

３．高速な，剛体分子動力学アルゴリズムの開発
氷と水，不凍タンパク質 è○○液 (with hybrid sim.)◉ DFTと古典MDとでは，同じ原子数なら7桁も計算時間が異なる◉ 第一原理MD è ３千原子まで
古典MD
è 1010原子まで

4ハイブリッド量子古典シミュレーション法

放熱材料稠密な固体内部での，電子と古典領域の力学的接続 (稠密固体に適したハイブリッド法)◉ 発熱部品（パワー素子）からの熱を筐体へ逃がす放熱材料◉ SiC (〜2020年度）利用に向け，更なる放熱性の要求◉ 密着性を保ちつつ熱伝導性向上したい è エポキシ樹脂に無機フィラを混入◉ 量子(QM)領域 à オーダーN割の実空間グリッド 密度汎関数法 (2千原子まで)◉ 古典(CL)原子領域 à 様々な経験的分子動力学ポテンシャル◉ 独自開発したbuﬀered cluster法による固体内部での接続
対象系 (半導体Si, C，金属Al，セラミックスalumina) およびQM領域形状に依らず適用可

例: グラファイト中のLiイオン群の熱拡散パワー素子の放熱経路と放熱材料
発熱部品
放熱材料

熱伝導筐体

例: 原子スケールの摩擦

フィラ充填率，表面結合剤の効果
表面結合剤と熱伝導率（実験）

冷却放熱材の熱伝導性が重要
放熱材料
C�Li(上)�

Si�

発熱部品

O�

熱伝
導H�

Li(下)�C�

フィラ（Al2O3など）
フィラ-樹脂界面
樹脂筐体（冷却ﾁｭｰﾌﾞ）

明色球
=QM原子

高フィラ充填率での，「アルミナ-樹脂-アルミナ」領域の熱伝導非平衡MDを用いたアルミナ-樹脂界面の熱伝導率
アルミナと樹脂(bisA)との界面シミュレーション
high-T heat flux

low-T

heat flux

放熱性能upの理論限界は？	

樹脂にフィラを入れて熱伝導性を向上

bisA
エポキシ分子high-Tfix

high-T

low-T

heat flux

アルミナ間に挟まれた樹脂部分の実効的な熱伝導率：

Λ=

fix

y0.15z alumina

界面での“温度”分布370K

w/o SCbisA(エポキシ)

vacuumyz

D

440

T (K)

bisA

SC molecules

0.09!T10.06

400

10
380

!T2T360

結合剤分子の導入により，界面熱抵抗が低減し，熱伝導率が40%向上

0.12T

420

Aluminaw/o SC
w/ SCvacuum

SC (16mols.)x

alumina

430K

aluminabisA

SC(表面結合剤)! (W/m•K)

aluminax[Heat flux](T1 − T2 ) / D

0102030203040 50
D (Å)

60

70

FIG.3

40

z-position (Å)
FIG.2

5060

樹脂部が薄くなるとΛが低下する理由：
高分子の配向具合が変化しΔT 増大

5固体材料の皮膜系，有機・無機材料との界面系

本日の内容

金属の塗装・接着部の物理化学的耐久性
金属: Al, Mg, Fe, Zn…
樹脂: エポキシ樹脂, PP…

１．電子系のための，オーダーN型の実空間グリッドDFTの開発
è 固体電解質皮膜-電解液界面でのLiイオン透過 è２次電池

アルミニウムエポキシ樹脂酸化膜樹脂
OHOHOO
自然酸化膜

２．稠密固体にも適用できる，ハイブリッド量子古典法の開発
è無機有機界面の熱伝導，物理化学的強度 è放熱材料，塗料

H2Oと反応?

３．高速な，剛体分子動力学アルゴリズムの開発

Al基金属

è氷と水，不凍タンパク質 è○○液 (with hybrid sim.)

QM領域
Work in progress…

New rigid-body MD algorithm: idea�

Simulation results of Quasi-Liquid-Layer (QLL) of ice�

◉ Rigid-body molecules: TIP4P potential with a long time-step = 5.0 fs◉ Fast time-reversible algorithm for rigid-body dynamics:
Fastest algorithm with explicit time-reversibility

H 2O

Ideas to ENFORCE time-reversibility�Angular position (quaternion): q
Angular velocity in lab-frame: ω#2(n)% q (n+1) = q (n) + Δtω (n) + Δt dω%2 dt$(−Δt)2 dω (n+1)% (n)q = q (n+1) + (−Δt)ω (n+1) +%2dt&◉ 1.31 million Water molecules◉ Top surface: (0001) surface	èHeight-difference (=bump height) ~ 12 ÅAngular velocity in body-fixed-frame: ΩMoment of inertia: IdΩi= (I j − I k )Ω j Ωk + torqueidtω = A[q]Ω

Redàbumps(凸)	

Tm-1K	

Tm: bulk melting temp.	
Tm-13K	

Tm-23K	Ii⇒ Get complex non-linear Eq. for (n +1)-th ω (n+1)⇒ Simplify the Eq. without losing time-reversibility èiteration free!	

6

Time evolution of Quasi-Liquid-Layer (QLL)-bumps of ice�

氷の擬似液体層(QLL)： 大規模シミュレーションでわかったこと

QLL-bumps at Tm-1K	

◉ 融点直下の温度では，凸部の擬似液体層とその下のシート状の疑似液体層が共存する

Red à Bumps	

QLL-bumps grow up & down

◉融点より6%程度低い温度では，凸部が移動しその直下で再結晶化する◉ 再結晶化の際，液体中の分子を取込む(環境から取込)可能性もある不凍タンパク質の作用メカニズムから新しい”産業品”へ

まとめ◉ 特殊なタンパク質を作り出し，極寒環境での凍死に耐えている変温生物◉ 水の固化点低下のメカニズムは，氷の曲率が鍵？

~25nm◉ 新しい○○液に展開可能かも？

Antifreeze proteins
H Nada and Y Furukawa

è 固体電解質皮膜-電解液界面でのLiイオン透過 è２次電池

２．稠密固体にも適用できる，ハイブリッド量子古典法の開発

ハイブリッド原子・剛体分子シミュレーション（100万個の水分子）

è無機有機界面の熱伝導，物理化学的強度 è放熱材料，塗料

タンパク質とその近傍：原子領域

３．高速な，剛体分子動力学アルゴリズムの開発
è氷と水，不凍タンパク質 è○○液 (with hybrid sim.)

遠方：剛体分子領域

Work in progress…

１．電子系のための，オーダーN型の実空間グリッドDFTの開発「低温科学」より

ater molecules forming a structure resembling ice at 4 ns for all conformations. The dashed lines show the initial positions

connecting a pair of water molecules
It is clear that AFP was partiallygrown during the simulation. As wehe surrounding ice had an importantding of AFP to the interface.
displacement Drz2 of the center of massn the z direction from the beginning ofshows the shift in rotation angle Df
FP as a function of t. These two panels
y of binding of AFP to the interface forn other words, initial conformation A
al position during the growth of ice,
B became stably bound to the interface
nformation A. Initial conformations C
initial positions, indicating unstable
A was deemed to be the most stable atce.

grown around AFP was estimated forgure 8). A large decrease in R occurred
nformations A and B, whereas R for
constant over the entire simulationicate that a large decrease in the Rding AFP that was stably bound to thedecrease in R likely corresponded to ice

surrounding AFP for conformations A, we suspected that the decrease in R7