Sep 4, 2020
37th CMD workshop@ online

Large-scale molecular dynamics simulationin conjunction with data-driven approach(10,368,000,000 iron atoms, 78h/500ps using 256GPUs)* A version for distribution
shibuta@material.t.u-tokyo.ac.jp
mse.t.u-tokyo.ac.jp

Yasushi Shibuta
Department of Materials Engineering
The University of Tokyo

Outline
1. Brief introduction of molecular dynamics (MD) simulation

2. Large-scale MD simulation of microstructure formation
・ Remaining problems in PFM and the role of MD
・ Large-scale MD simulation of nucleation and solidification
・ Evolution of multi-scale modeling

3. Inter-scale modeling: share of information
・ Conversion of atomistic configuration into phase-field profile
・ MPF simulation of grain growth starting from MD configuration

4. Cross-scale modeling: utilization of shared information
・ What is data-assimilation(DA)? - Bayes' theorem ・ Concept of data assimilation between MD and PFM
・ On-the-fly estimation of solid-liquid interfacial properties by DA

5. Summary with some ongoing works

Materials modelling laboratory– modelling for materials design –What is desiredfor materials modelling?

EtOH dissociation
on Fe cluster *1
CNT growth*2

Nucleation and solidification(10B atoms)*3

New materials designfrom computational simulations*1CPL 731 (2019) 136619*2CPL 382 (2003) 381*3MSMSE 27 (2019) 054002*4CMS, 152 (2018) 118.

Ab initio MD（electronic）
Grain growth*4

Classical MD（atomistic）
Convolutional Neural Network*5*5CMS, 184 (2020) 109880.

Microstructure formation
Heat and masstransformation
Phase-field method（Meso-scale）
Continuum dynamics（macro）

Bayesian inference

How to collaborate with
data-driven approach?

Molecular dynamics (MD) simulation
In the MD simulation, trajectories of all atoms in materials are determined bynumerically solving Newton's equations of motion for a system of interacting atoms.

Governing equation: Newton’s second law of motion

d 2r
m 2 Fdt

MD simulation is an deterministic method, in which position and velocity of all atoms
at every time step is determined explicitly, when initial position and velocity andas well as the force acting on all atoms is defined.
例） Snapshot of Ar216 cluster at 50 K and trajectory of sky-blue atom during 100 ps.
[y]
65

50ps
100ps

60

0ps

25ps

55
55

60

65 [x]

Force acting on atoms
1. Force from (time-independent) Schrödinger equation→ab initio MD

Hˆ   E

Under adiabatic approximation, the electronic state (wave functions) are calculated by
ab initio method. Then, the Hellmann–Feynman force is derived from the wave function.


H
F *

r

2. Force from the interatomic potential (force field)→classical MD
The interaction between atoms are approximatedby the potential function (the force-field).
0.01

7

F
  12   6 0
 (rij )  4      
i
F
m
rij
 rij   rij  
12
6

–0.01
   (rij )j    
f ij   242
  
 rij  rij
rij   rij 
  

3
The force is derived from the differential of the potential.

100dE/dr: Force (10 eV/m/atom)m

E: Potential energy (eV/atom)

Ex) Lennard-Jones potential

Lennard–Jones Potential for Ar

0
Repulsive

Attractive

4
Interatomic distance (Å)

Position and velocity are derived from the Newton's equationsof motion even if the force is given.
It is same process both for ab initio and classical MD.

5

d 2r
m 2 Fdt

Modelling of metal-catalyzed growth of carbon nanotube and graphene
Carbon cap formation on Ni cluster

Precursor (EtOH) dissociation
on transition metal (Fe, Co) cluster

20 nm

on Fe cluster

1.2nm
on Co cluster

Classical MD

Ab initio MD
S. Fukuhara, M. Misawa, F. Shimojo, Y. Shibuta,
Chem. Phys. Lett. 731 (2019) 136619.
Y. Shibuta, S. Maruyama,
Chem. Phys. Lett, 382 (2003) 381.

A proper method should be chosen
depending on the purpose of the study.

Outline
1. Brief introduction of molecular dynamics (MD) simulation

2. Large-scale MD simulation of microstructure formation
・ Remaining problems in PFM and the role of MD
・ Large-scale MD simulation of nucleation and solidification
・ Evolution of multi-scale modeling

3. Inter-scale modeling: share of information
・ Conversion of atomistic configuration into phase-field profile
・ MPF simulation of grain growth starting from MD configuration

4. Cross-scale modeling: utilization of shared information
・ What is data-assimilation(DA)? - Bayes' theorem ・ Concept of data assimilation between MD and PFM
・ On-the-fly estimation of solid-liquid interfacial properties by DA

5. Summary with some ongoing works

Phase-field, Phase-field model
Solid
=1

diffuse interface

Phase-field,   2 2  f  , c tSolid-Liquidinterface

Driving force, Gibbs-Thomson effectLiquid=0Concentration, cc M    , c tDiffusion in bulks, Stefan condition

Solid :  = 1
Temperature, T

Liquid :  = 0

T

t

Velocity, v

v

t

Remaining problems in PFM and the role of MD
Now, PFM can treat large-scale 3D simulation of solidification and grain growthRemaining problems in PFM
1. There is no explicit way to treat a nucleation event.
→Nuclei specified in advance as having a random (or a particular) distribution
or forcibly formed in line with an assumption
2. It is not easy to treat the anisotropy in grain boundary properties.
→Mostly the approximate model (such as Read-Shockley) employed.
Ex) Initial condition of PFM ofcompetitive dendrite growth.

Ex) Multi-PFM of grain growth

Figures from Prof. T. Takaki (KIT)Can we treat these problems by MD simulation?
・Spontaneous nucleation occurs when an enough spatiotemporal scale is given.
・The anisotropy in grain boundary properties is inherently taken into account.

Polycrystalline microstructures by Laguerre-Voronoi tessellation

Log-normal distribution

Voronoi division

Polycrystalline microstructure

Z. Fan, Y. Wu, X. Zhao, Y. Lu, Comp. Mater. Sci. 20 (2004) 301.

S. Falco, J. Jiang, F. De Cola, N. Petrinic, Comp. Mater. Sci. 136 (2017) 20.

How to include the orientational relationship between grains?
How to express the anisotropy in grain boundary properties?

Homogeneous nucleation and microstructure evolution
in million-atom MD simulation (quasi-2D)
1,037,880 atoms (53.4x53.4x4.3nm) 1400K(0.58Tm)

0.67

1600

0.58

1400

0.50

1200

0.42
0
10

1

1000

Classical nucleation theory
for homogeneous nucleation
Temperature [K]Normalized temperature [K/T m]

Nucleation rate

3
 16


SL

I  I 0 exp 2 3 ( s f T ) k BT 

Two competing effects in the exponent
result in the nose shape

10
–1
–1
Nucleation rate [(cell size) ns ]Y. Shibuta, K. Oguchi, T. Takaki, M. Ohno,
Scientific Reports, 5 (2015) 13534.

Thermally activated homogeneous nucleation occurs
in MD simulations without any inducing factor.

Billion-atom MD simulation by multi-GPU parallel computation
pure Fe (FS potential)
Y. Shibuta, S. Sakane, E. Miyoshi, S. Okita, T. Takaki, M. Ohno
Nature Communications, 8 (2017) 10.

TSUBAME2.5
512GPU(c.a.10ps/1h)(1,040,000,000 atoms)

Heterogeneity in homogeneous nucleation
Y. Shibuta, S. Sakane, E. Miyoshi, S. Okita, T. Takaki, M. Ohno
Nature Communications, 8 (2017) 10.

From experiments: the origin of peak from finemicrostructure, abundant twin boundary

Orientational relationshipbetween neighbor grains
Peak around 60 appears

Au-Cu-Ag

Au-Cu-Ag + 200ppm Ir

G. Kurtuldu, A. Sicco, M. Rappaz, Acta Mater 70 (2014) 240.

Micrometer-scale MD simulation of microstructure formation
Y. Shibuta, S. Sakane, E. Miyoshi T. Takaki, M. Ohno, Model. Simul. Mater. Sci. Eng., 27 (2019) 054002.

10,368,000,000 iron atoms in 1.03 x 1.03 x 0.11 m3
FS potential (pure iron)78h/500ps using 256GPUs(P100 on TSUBAME3.0)

Evolution of multi-scale modeling“Advent of cross-scale modeling: high-performance computing of solidification and grain growth”
Y. Shibuta, M. Ohno, T. Takaki, Advanced Theory and Simulations, 1 (2018) 201800065.

・Parameterization
Ab initio (Electronic)↓
MD/MC (Atomistic)↓
PFM (Continuum)
Averaged “static” information

・share of information
・validation of models

・utilization of shared information
Ex) Data-driven MD/PFM

Ex) direct comparison
between MD and PFM“On-the-fly” use of information

Outline
1. Brief introduction of molecular dynamics (MD) simulation

2. Large-scale MD simulation of microstructure formation
・ Remaining problems in PFM and the role of MD
・ Large-scale MD simulation of nucleation and solidification
・ Evolution of multi-scale modeling

3. Inter-scale modeling: share of information
・ Conversion of atomistic configuration into phase-field profile
・ MPF simulation of grain growth starting from MD configuration

4. Cross-scale modeling: utilization of shared information
・ What is data-assimilation(DA)? - Bayes' theorem ・ Concept of data assimilation between MD and PFM
・ On-the-fly estimation of solid-liquid interfacial properties by DA

5. Summary with some ongoing works

Connection between MD and PFM
Our large-scale MD (113,246,208 atoms) （50days by 128GPU）

S. Okita, E. Miyoshi, S. Sakane, T. Takaki, M. Ohno, Y. Shibuta, Acta Mater. 153 (2018) 108.

Multi-phase-field simulation（less than 1day by 1GPU）

e.g. E. Miyoshi and T. Takaki: Computational Materials Science, 120 (2016) 77.
E. Miyoshi and T. Takaki: Computational Materials Science, 112 (2016) 44.

Bridging MD and PFM for grain growth prediction
E. Miyoshi, T. Takaki, Y. Shibuta, M. Ohno, Comp. Mater. Sci. 152 (2018) 118.

MD

voxel

PFM

Atomistic configuration from MD is converted into phase-field profile.

Bridging MD and PFM for grain growth prediction
E. Miyoshi, T. Takaki, Y. Shibuta, M. Ohno, Comp. Mater. Sci. 152 (2018) 118.

MD
(atomistic)

MPF(continuum)

Conventional multiscale modeling:
Averaged values from smaller systemare extracted as parameters
and transfer to larger system.
New concept of multiscale modeling:
Direct mapping from bottom level

MPF simulation of grain growth starting from MD configuration
Y. Shibuta, S. Sakane, E. Miyoshi T. Takaki, M. Ohno, Model. Simul. Mater. Sci. Eng. 27 (2019) 054002.

Conversion of atomistic configuration into phase-field profile

Subsequent grain growth by MPF

Outline
1. Brief introduction of molecular dynamics (MD) simulation

2. Large-scale MD simulation of microstructure formation
・ Remaining problems in PFM and the role of MD
・ Large-scale MD simulation of nucleation and solidification
・ Evolution of multi-scale modeling

3. Inter-scale modeling: share of information
・ Conversion of atomistic configuration into phase-field profile
・ MPF simulation of grain growth starting from MD configuration

4. Cross-scale modeling: utilization of shared information
・ What is data-assimilation(DA)? - Bayes' theorem ・ Concept of data assimilation between MD and PFM
・ On-the-fly estimation of solid-liquid interfacial properties by DA

5. Summary with some ongoing works

Evolution of multi-scale modeling“Advent of cross-scale modeling: high-performance computing of solidification and grain growth”
Y. Shibuta, M. Ohno, T. Takaki, Advanced Theory and Simulations, 1 (2018) 201800065.

・Parameterization
Ab initio (Electronic)↓
MD/MC (Atomistic)↓
PFM (Continuum)
Averaged “static” information

・share of information
・validation of models

・utilization of shared information
Ex) Data-driven MD/PFM

Ex) direct comparison
between MD and PFM“On-the-fly” use of information

What is data-assimilation?
If the simulation model covers all observed phenomena completely,
the simulation can reproduce the observation exactly. However, it is impossible..
Therefore, the system noise is added to have a range in the answer
in order to compensate the incompleteness in the simulation modelor the numerical error due to the discretization of time differential.
→ called a system model
System model

Simulation modelxt  ft (xt 1 )p(xt | xt 1 )   (xt  ft (xt 1 ))xt
・Deterministic
・Delta function
as probability distribution

p(xt | xt 1 )x t  f t (x t 1 )  v t

xt

v t ：system noise・probabilistic・having a range
as probability distributionp ( A | B ) is the conditional probability: probability of A given B.
Here, p ( A, B )  p ( A | B ) p ( B )↑ joint probability for A and B

What is data-assimilation?
Moreover, we want to increase the accuracy of simulation results
by importing observation data (i.e., data assimilation).
p(ht (xt ))p(ht (xt ) | y t )PDF

System modelht (xt )ht (xt )yt ●y t ：observation data

Comparison with observation datayt ●Assimilation into observation data

In the system model, there is a probability function distribution (PDF) for the answer.
It is possible to increase the accuracy by revising the PDF.

This is basic concept of data assimilation. But how to assimilate?

Bayesian inference!
Practically, several filtering techniques are proposed:
Particle filter, Kalman filter, Ensemble Kalman filter, etc.

Bayes' theoremp ( A, B)  p ( A | B) p ( B)  p ( B | A) p ( A) derives the following relation.
Bayes’ theorem

likelihoodp ( B | A) p ( A)p( A | B) p( B)Thomas Bayes(1702-1761)

prior probability
Probability of A

posterior probability
Probability of A given B

Very simple but very important theorem→implying the derivation of cause (distribution) from the effect (observation)
For example, assuming A as cause and B as effect,p (effect | cause) p (cause)p (cause | effect ) p (effect )“From an effect obtained from the observation,
we can know the probability of the cause.”From 意味が分かるベイズ統計学：
一石賢，ベレ出版，２０１６

Bayes' theorem: example

Let us consider an airplane accident assuming following three reasons.
Event probability

Probability
to cause accident

Human error

0.6

0.02

Maintenance error

0.3

0.03

Error in control tower

0.1

0.01

When one accident happens, what probability is it for the cause as “Error in control tower”?
p (acc. | ect ) p (ect )p (ect | acc.) From Bayes' theoremp (acc.)ect: error in control tower
Acc.: accident
Maintenance error
Error in
Human error
control towerNow

0.01  0.1
 0.0450.02  0.6  0.03  0.3  0.01  0.1
p (acc.)  p (acc. | hue) p (hue) p (acc. | mae) p (mae)hue: human error
mae: maintenance error p (acc. | ect ) p (ect )Here,p (effect )  p(effect | cause) p(cause)cause

p (Ck | E ) 
0.03

0.02
0.6

0.3

0.01
0.1

p ( E | Ck ) p (Ck )p ( E | Ck ) p (Ck ) np( E )p ( E | Ci ) p (Ci )i 1Bayes’ expansion formulaSchematic image of data assimilationp (xt 1 | y1:t 1 )p (x t | y1:t 1 )p (xt | y1:t )p (xt 1 | y1:t )Predicted distribution: prediction of current state from observation data until t-1



p (x t | y1:t 1 )  p (x t | x t 1 ) p (x t 1 | y1:t 1 )dx t 1
obtained by system modelx t  f t (x t 1 )  v t v t ：system noise
Filtered distribution: prediction of current state from observation data until now (t)Likelihoodp (x t | y1:t ) 

p (y t | x t ) p (x t | y1:t 1 )p (y t | x t ) p (x t | y1:t 1 )dx t

Bayesian inference

Practically, several filtering techniques to define likelihood are proposed:
Particle filter, Kalman filter, Ensemble Kalman filter, etc.

Application example of data assimilation: weather forecast

How is the weather forecasted?

Referred from tenki.jp web(For Typhoon No. 18 at Sep, 2013.)

Forecast circle changes as the time goes.
→observation data is used for assimilation.

Interesting points
・Where the typhoon comes from now?
・Size of typhoon
Not so interesting (for most of people)
・Where and when the typhoon starts?
（i.e. the initial condition）
・model of forecast（NS eq.? What difference scheme？）

Data assimilation for estimation of grain boundary properties
Schematic image of EnKF

Parameter estimation for MPF simulation from 3D/4D experimental data.
Ensemble Kalman filter-based data assimilation for three-dimensional multi-phase-field model:
Estimation of anisotropic grain boundary properties
A. Yamanaka, Y. Maeda, K. Sasaki, Mater. Des. 165 (2019) 107577

Can we estimate interfacial parameters PFM dynamically from MD data?

Concept of data assimilation between MD and PFM
for estimation of solid-liquid interfacial properties
MD

50 nm

50 nm
LiquidLiquid

Solid

Solid

(T = 200K)
Isothermal solidification of Fe(T = 200K)
Isothermal solidification of Fe
M. Ohno et al.
Phys. Rev. E 101(2020) 052121.

Data assimilation (EnKF)Interfacial energy,  (J/m2)

Estimation of SL interfacial energy
0.6
0.5
0.4
0.3
0.2
0.1
0.0

0

PFM

50

100

150
Time, t (ps)

200

250

300

* Time-variable properties during
microstructure formation can be
estimated from MD profile.
* Prior knowledge is not required for
estimation of time-variable properties

New concept of collaboration
between MD and PFM~ cross-scale modelingBayesian inference of interfacial properties out of equilibrium
Data assimilation based on
Ensemble Kalman filter(EnKF)MD results employed as observation data

Simultaneous estimation of solid-liquid interfacial
energy, interfacial mobility and their anisotropy

M. Ohno, Y. Oka, S. Sakane, Y. Shibuta, T. Takaki,
Phys. Rev. E 101 (2020) 052121.

Procedure of on-the-fly parameter estimation
Y, Nagatsuma, in graduation thesis, The University of Tokyo, 2020 (unpublished)

A simple phase-field model (PFM)
𝜕𝛟
𝜕
𝜕𝒂 𝜕𝝓
𝜕
𝜕𝒂 𝜕𝝓
= 𝑴𝝓 [𝛻 𝒂2 𝛻𝝓 −𝒂+𝒂
+ 4𝑾𝝓(1 − 𝝓)(𝝓 − 0.5 + 𝜷)]
𝜕𝒕
𝜕𝒙
𝜕𝜽 𝜕𝒚
𝜕𝒚
𝜕𝜽 𝜕𝒙

𝑎 θ = 𝑎[1
ത + 𝜁𝑐𝑜𝑠{𝑘(𝜃 − 𝜃0 )}]Data assimilation based on
Ensemble Kalman filter (EnKF)

Anisotropy parameter ,inclination of solid crystal 0,change in the driving force, 
are simultaneously estimated.

1. Parameter estimation using synthetic
observation data (Twin experiment)

2. Parameter estimation using MD data
as observation data

Observation data prepared by PFM
simulation with known parameterswith some random noises- true value of parameter known in advance

MD data used as the observation data
- true values of parameters unknown in advanceOn-the-fly parameter estimation

Preparation of MD data for single Ni crystal growth
Y, Nagatsuma, in graduation thesis, The University of Tokyo, 2020 (unpublished)

Preparation of calculation system (pure Ni)

[010]
[110]
[100]

10 nm

insertLiquid(c.a. 2,700,000 atoms)

10 nm

3.5 nmSolid nucleus

Calculation condition and snapshots of MD results
1530 KInteratomic potential:
EAM* by Zope & Mishin*R.R. Zope & Y. Mishin,Phys. Rev. B 68, 024102 (2003).

Using LAMMPSVisualization: OVITO(using a-CNA)
T = 1530 and 1470 K

[010]
[110]
[100]
1480 K

200 ps400 ps600 ps800 psSimultaneous estimation of three parameters (，0，)Synthetic observation data from PFM(twin experiment)

MD data as observation dataestimated value:
0.249  0.01

True value: 0.2

0True value: 0.314

0estimated value:
0.028  0.004
estimated value:
0.371  0.004

True value: 0.35

βAll parameters converged to true valuesβ

All parameters converged to specific value

Observation data

Observation data

Estimated result

Estimated result

Succeeded to estimate three parameters simultaneously.
Y, Nagatsuma, in graduation thesis, The University of Tokyo, 2020 (unpublished)Outline
1. Brief introduction of molecular dynamics (MD) simulation

2. Large-scale MD simulation of microstructure formation
・ Remaining problems in PFM and the role of MD
・ Large-scale MD simulation of nucleation and solidification
・ Evolution of multi-scale modeling

3. Inter-scale modeling: share of information
・ Conversion of atomistic configuration into phase-field profile
・ MPF simulation of grain growth starting from MD configuration

4. Cross-scale modeling: utilization of shared information
・ What is data-assimilation(DA)? - Bayes' theorem ・ Concept of data assimilation between MD and PFM
・ On-the-fly estimation of solid-liquid interfacial properties by DA

5. Summary with some ongoing works

Calculation condition
Interatomic potential：EAM**T. T. Lau et al., PRL. 98 (2007) 215501

Fe：250 atoms C：1 atoms
300 K, NVT const.
CV: bond length

MD(10ns)

hypertime [ns]Collective variable-driven hyper dynamics (CVHD)
for acceleration of diffusion MD at room temperature

Original MD time [ns]MD＋CVHD(1.4×107ns)

S. Fukuhara, K.M. Bal, E.C. Neyts, Y. Shibuta, Comp. Mater. Sci., 177 (2020) 109581.

Keyword: metadynamics, hyperdynamics

Machine learning approach to automated analysis of
atomic configuration of MD simulation
Local atomic structures of solid-liquidbiphasic system are identified by 3D
convolutional neural network (CNN).

Accuracy of 3D-CNN is more than
90% for all structures.
T. Fukuya, Y. Shibuta, Comp. Mater. Sci., 184 (2020) 109880.

Summary
Based on large-scale MD simulations on microstructure formation,
Several new trials on multi-scale modeling are introduced.
Inter-scale modeling: share of information
・Direct mapping of atomistic configuration into phase-field profile
・PFM simulation staring from the microstructure by MD simulation

Cross-scale modeling: utilization of shared information
・On-the-fly estimation of SL interfacial properties based on data assimilation

Data-driven technique is opening new stage of multi-scale modeling

Acknowledgement
Tomohiro Takaki, Shinji Sakane, Eisuke Miyoshi (Kyoto Inst. Tech.)
Munekazu Ohno (Hokkaido Univ), Tetsuo Mohri (Tohoku Univ.)
Shin Okita, Kensho Ueno, Takuya Fujinaga, Yuhi Nagatsuma (Univ. Tokyo)
Financial support by
Grant-in-Aid for Scientific Research (B) (16H04490, 19H02415) from JSPS.
Support for computational environment by
JHCPN, HPCI, CDMSI (post K 7E)