2019.09.06

Applications of Materials Informatics
For High-productive Research:
Machine-learning potential for amorphousand spectrum-data analysis

Research Center for Computational Design of Advanced Functional
Materials (CD-FMat), AIST
Senior Researcher

Yasunobu Ando1

Self-Introduction
Name：Yasunobu ANDO（36）CV✓ 2012. 3.
Ph.D in physics @ Dept. Physics, UTokyo.
✓ 2012. 4 ~ 2013. 4
Posdoc @ AIST✓ 2013. 5 ~ 2016. 3
Assistant professor @ Dept. Materials Eng., UTokyo✓ 2016.4 ~Researcher @ CD-FMat, AIST

Major topics

Materials x Informatics(Physics・Chemistry・Informatics)!2

Ab initio MD
in Solid-Liquid Interfaces

ute of Advanced Industrial Science and Technology - Google Maps

2017/11/16 13)52

Place where I come fromAIST Tsukuba Central, National Institute of Advanced Industrial Scienceand Technology

Employees and Budget● Researchers (foreign nationals)……2,284(116)［ Permanent ］［ Fixed term ］

Employees and Budget[ 1,925 ][ 359 ]

● Administrative employees (foreign nationals)………686(1)Total number of employees : 2,970(117)● Executives (full time) …………………………13● Visiting researchers ………………………185● Postdoctoral researchers …………………190● Technical staﬀ ……………………………1,487

Num
indu
●C
●U
●O

(As of July 1, 2016)● Researchers (foreign nationals)……2,284(116)［ Permanent ］［ Fixed term ］[ 1,925 ][ 359 ]

● Administrative employees (foreign nationals)………686(1)Total number of employees : 2,970(117)● Executives (full time) …………………………13
Kartendaten © 2017 Google,ZENRIN
20 km● Visiting researchers ………………………185● Postdoctoral researchers …………………190
3

Number of researchers accepted throughindustry/academia/government
Composition of partnerships
researchers by d
● Companies …………………………………1,856● Universities National………………………………1,924
D
Metrology Institute of Japan● Other organizations ………………………936(foreign nationals : 456)
(Total number of researchers accepted in FY 2015)14 ％Geological Survey of Japan

10 ％

18

Major research topics

My interests:• Electric-double layer• Reaction on surface and interfaces.

Tool: Density-functional-theory based (MD) simulation.
4

73
145

We need “NEW” tools

TABLE I. Explanatory variables: The
mechanical computations from the str

Spectra-data sets including 600 samples
Now we “can” obtain hundreds
thousands spectra from simulation
FPMD of electrochemical Interfaces

Efficient Analysis is necessary.
5

ace open innovation. Rapid
ional modeling and data
nced algorithms for modelingbe developed to supplementnd a data exchange systemhers to index, search, andmplemented to allow greater
ation.

data sharing systems and more integrated engineering
teams will allow design, systems engineering, andmanufacturing activities to overlap and interact(see Figure 2).

Milestone project with Informatics

This new integrated design continuum — incorporating
greater use of computing and information technologies
coupled with advances in characterization andexperiment — will significantly accelerate the time andinuum are necessarily linear
number of materials deployed by replacing lengthy
occur before systems design),
and costly empirical studies with mathematical models
benefit from increased dataand computational simulations. Now is the ideal time
mmunication. Currently, noto enact this initiative; the computing capacity
low
necessary to achieve theseamsadvances exists and related
Future Materials
Number ofels.
technologies
suchasContinuum
New Materials
to Market
may
nanotechnology and bioact
technology have matured toeenenable us to make great
act
progress in reducing time
Materials Continuumonto market at a very low cost.
TodayofMultiple international entities
uct
have
recognizedtheseaissues and a number of
tain
Time
foreigncountries
havebealready
embarkedon
has
Figure 2: Initiative acceleration of the materials continuum
programs to address them.6
nedThe
National
Researchno1. DevelopingaMaterials
Innovation
Infrastructure
Councilofthe
National
Academiesof
Sciences,in itsrch for it. Data transparency
reporton
Integrated
Computational
Materials
stryandmultidisciplinary
2. Acheiving
National Goals With Advanced Materials
Engineering, describes the potential outcome:Materials Genome Initiative（MGI from 2011）

Main purpose

the Next-Generation
Materials
Workforce
Integrating
materials computational tools
andmaterial does not 3.Equipping
end with
information with sophisticated computational andhat is coming more to theanalytical tools already in use in engineering
consumers is the recyclability
fields… [promises] to shorten the materials
aterials. Materials engineersdevelopment cycle from its current 10-20 years to
er-changing parameters and6
2 or 3 years.7
heir initial intended purpose;

Using Database and machine-learning

Our research in MI
1. Clustering and correlation analysis of transmission spectrum of molecular junction system
2. A new descriptor of perovskites for performance estimation of fuel cells
3. Materials Search for a well-worked 2D substrate for Germanene and Stanene
4. Unsupervised clustering of PDOS in Surface system
5. Development of machine-learning potentials for Amorphous research
6. Yield prediction in experiments from Simulation (Catalyst informatics)
7. High-throughput peak fitting on many XPS spectra
8. Model selection of equivalent circuits on impedance spectroscopy
9. Model selection of preferred orientation distribution of tourmaline-grains
10.Parameter optimization of equation of states for an inner earth environment!7

So, what is informatics?
What is machine-learning?

8

Short introduction
Functions of ML✓ Prediction✓ Characterization✓ Classification✓ Pattern Recognition

Linear Regression (PCA)Y =X

ai Xi + bi

Support Vector Machine（SVM）

What they do is just“Putting points andDrawing lines”
9Representing Complex Situation
NOT work lines ?
JUST draw the “curve”✓ Basis expansion (polynomial, spline etc.)✓ Kernel Regression, Gaussian Process✓ Neural Network

Fitting process by Neural network

Applicable to high-dimensional data10

Exploration and Exploitation
Bayesian Optimization✓ Based on the observation, predicting values and its “credibility” at unobserved points✓ Automatically searching observable space considering predicted value and credibility.
11

4

1.56

PbICl

Pnma

0.18

2.72

pounds with high Z
51.56
PbClBr
Pnma
0.09
3.44
which are neighboring
7
1.44
PbI2
0.29
2.42
R 3̄m
Bi2 Te3 is a famous
P63 mc
0.29
2.45
8
1.43
PbI2
s known that some121
0.39
K2 CdPb
Ama2
0.45
0.18
LTC. This may be
144
0.29 Cs2 ½PdCl4 "I2 I4=mmm
0.31
0.88
ng dataset composed
nts and three kinds ofwords, the training
conductivity, electrical conductivity, and the Seebecks unavoidable at thecoefficient, should also be optimized. Although they can
C calculations are still
be tuned by elemental doping, the band gap, Eg , should be aly unbiased training
week ending
simple measure of the electronic structure and allows us to
P
H
Y
S
I
C
A
L
R
E
V
I
E
W
L
E
T
T
E
R
S
of data to
cover
the
13 NOVEMBER 2015
PRL 115, 205901
discriminate(2015)
in a simple way between materials that can be
Jpn. J. Appl. Phys. 55, 045502 (2016)
S. Kiyohara et al.
nd crystal structures.
good thermoelectrics or not. All of the 221 compounds
set will be verified to(a) showing a positive Z score are(b)listed in the SM [33]elemental descriptors is found to improve the robustness of
50
ls. Because of the use
together with Eg (DFT-PBE) given in the MPD library.
P
b
S
e
L
i
I
may not be able tothe efficient search.
Among them, only 19 compounds satisfy 0.1 <rials in the library.
Eg < 1.0 eV. First-principles
LTCs
areevaluated
for
them.
Jpn. J. Appl. Phys.Kriging
55, 045502(2016)
S. Kiyohara et al.
with
2 descriptors
As described in the Supplemental
Material (SM) [33],part of them.
Crystal structures and LTC for two of
them
are
shown
in
+ 34 elemental descriptors
he Z score does not Kriging
better correlations with LTC can be found for parameters that
descriptors
Fig. 4 withand 2Table
I. Both K2 CdPb and
2 ½PdCl4 "I2 are(Nav =Cs19)
st-principles ranking. (N=11)
predicted
to exhibit LTC of less than 0.5 W=mK (at 300 K)
av
10
are obtained from the phonon density of states. However,the candidates of thetogether with a band gap of smaller thanRandom
1 eV. Thel screening is one of
we do not use such phonon parameters as predictors in the
Kriging with
2 descriptors
discoveryof suchcompounds may open a(Ngate
toward= 50)
avring” low-LTC 5
com+
34
elemental
descriptors
designing new thermoelectric materials with an exceptionpresent study because there is no data library available foruated for the top eight= 19)
ally(Nhighfigure of merit.
av
Kriging
2 descriptors
ng. All of them aresuch phonon parameters for a wide range of compounds.
In this Letter,
we first report the theoretical
LTCwithof 101
Random(N=
65)
LTC calculation wasav
compounds by(Nfirst-principles
anharmonic lattice-dynamics
Hereafter, we show results only with the predictor set= 55)
av
e presence of imagicalculations. Using
these data, the virtual screening of acomposed of 34 elemental descriptors on top of V and ρ.
ell used in the present
library containing 54 779 compounds is performed by
1
TC of the rest of theBayesian optimization using the kriging method based
Screening for low LTC compounds over compounds in a
All of the top five 0 on20the Gaussianwith 80 100
40 60 process
80 regressions.
100 0 221
20 materials
40 60
large library is carried out using a GPR prediction model.
W=mK (at 300 K),
very lowof LTC
are found from thisNumber
screening.
A final
Numberobservationsofobservations
e rocksalt PbSe, i.e.,filtering of those low-LTC compounds is made using theSuch a screening based on a prediction model is called as the powerfulness of“virtual screening” in biomedical communities [34]. For
FIG.discov2 (color online).
Lowest LTC(b)values
at each iteration in theor efficiently(a)I
structuresof highly
the virtual screening, we adopt all 54 779 compounds in the
krigingsearch for finding (a) PbSe Cland (b) LiI. Those random
PbRb4 Br6 , and PbI2
K
searches
When performing a
MPD library [30,35], which is composed mostly of crystal
Pd
PbbICl and PbClBr
have are also shown for Cdcomparison.
Cs
krigingfor finding LiI, PbSe and PbTe are intentionally
PbI2 (R 3̄m)
and PbIsearch
2structure data available in ICSD [28]. This means that mostcking sequences.
Allomitted and the rest of the 99 compounds are used.
ofthem have been synthesized experimentallyat least once.
sixfold or eightfold
Fig. 4. (Color online) (a) Trajectory of the calculated GB energy to theFig. 2. (Color online) Stable atomistic arrangements of Σ5[001]=(210)convergence.
The “conventionalwith (a)thethe conventional
methodand (b)GPR
the krigingprediction
method. Dashed line modelare of stoichiometricOn
basisofGB.Selected
themadeby V, ρmethod”
andmeans the computation of all
K CdPbCs [PdCl ]I
conﬁgurations. To understand the calculation eﬃciency easily, the calculated
represents
the (Color
positionof the(a)Polyhedronshowsduring
the six-memberedFig.
3.
online)data
pointsthekrigingsearch.
valent.
resultby the conventional
method (black)structureunit. correspondsdescriptors
Theelemental
numberto the trial number for the
kriging.
The 101
red number
34
for
the
LTC
data,
a rankingforwas plotted by order of the GB
varianceFIG.
are4 described
Gaussiankernelfuncnsidered its
for thermoenergy, whereas it was plotted by order of the trial number in the kriging
indicates that the initial calculations. (b) Magniﬁed image around (x = 5.0,
(color online). using
Crystal thestructuresof K2 CdPb
andmethod (red). (b) Magniﬁed image around the area pointed by blue arrow in
y = 1.2, z = 1.8) in(a). The number “69” (indicatedby pink arrow)denotes
d to electronic
struclow-LTC
compounds
is madeaccordingto the Z score ofCs2 ½PdCl4our"I2 predicted
to showthelow LTCof <0.5 W=mK
tion.
Therefore,
GPR
has
two
free
parameters,i.e.,the (a).
the converged structure, i.e., the stable structure.
ribution of thermal(at 300 K) and narrow band gap of < 1 eV.
3. Results and discussion

Minimum of LTC (W/m·K)

Material Search with Bayesian Opt.

2

2

4

2

the 54 779 compounds.
variances of the Gaussian kernel and prior distribution. Here,A GB
of copperSince
wasdepends
chosenonas thethe ﬁrst
test random
casecalculations.
its eﬃciency
been reported
The kriging method was applied
Figure3Σ5[001]=(210)
shows
the
distributionof Z scores
forpreviously.
the 54 779
andselection
the stablestructurewasinvestigated.
Using
the
conventhey are
given as 20 and 0.1, respectively.
205901-4of data points, we repeated the kriging operation by to search for the most stable conﬁguration, just as in the casetionalapproach,allsame
conﬁgurations
were Vcalculated
thecompounds
along
with
ρ.space
The magnitudeofcoppertheGBZdescribed above. The result
1. (Color online) (a) Search space. Each grid point corresponds to the74 times
for theGB. Theresultsofallandtrials and
are
shown
of the Σ5[001]=(210)
Figure 2(a) shows the Fig.
resultofthekrigingsearchinmoststablepoint
wasdeterminedfromthesearchpossible rigid body translations. (b) GB energies in the search spaceas histograms in Figs. 5(a) and 5(b). The 42 kriging opera- was that we succeeded in determining the stable structure afterscore
is completed
plottedin panels
corresponding
to constituent
shown
Fig.
1(b).
The obtained
shown
determined by comprehensive calculations.
comparison to the random
search of the lowest LTC
tionsinwere
by fewerstable
than structure
50 time iscalculation,
105 calculations. The calculated structure is shown in Fig. 6
in whereas
Fig. 2(a),more
and thecalculatedGB
energy
was
0.96
J=m
.
than 200 time calculations
wererequiredin elements(purple circles). The
stablestructure determined by the kriging
elements.
(Transition
metal
and
are
shown
As thecan4bekriging
seen fromthe ﬁgure,GB
isof composed
ofotherancompounds within the 101 compounds.
The
average
numJpn.
J. Appl.
Phys.
(2016)
operations.
The thenumberrequired operation
approachagrees wellwith that55of theprevious 2–6.
study [Fig. 6PRL
115,
205901(2015).
array
with
a
six-memberedstructure
unit;
this
structureagrees(3) Based on the joint probability distribution, the possible intheSM[33]).
The
Z
score
is
relativetorocksalt
PbSe,
becomes
large
when
the
structuresatalimitedarea
inthe(orangecircles)].
Since
17,466
conﬁgurations
are present forbers of compounds required forthe optimization
using
the
the previously
reported
structures.by theHowever,GB energies
are estimated at all points
in thesearch with"whole
data space
was selected
random 17,983
structure Science
this !3½110"=ð
111Þ
GB
bcc-Fegrainboundary,
the kriging
Adv.,2, twoe1600746-1-7.
space.
complete
calculations
werenecessary
toin reachthisspacestableshowingthelowestLTC
among
themethod101
compounds.
selection.
Suchnon-uniform
selection
the
data
leads
again
achieves
ordersofmagnitudehighereﬃkriging and random searches,
N
,
are
11
and
55,
respec12
av calculated by the following equation, are structure by the conventional method.
(4) Z-scores,
to
the
small
varianceofthedescriptors.
Astheresult,theciencythantheconventional
method.
This
clearly
indicates
54for779
compounds,
221 compounds, which
estimated at all points in the search space. Z ! score ¼ Among
To acceleratethe
the search
the stablestructure, the kriging

Thermoelectric materialInterface matching
26)

2

26)

18,20)

Energy storage: Battery materials

from manufacturing equipment and processes.46 Indeed,
Energy storage materials such as in Li ion batteries rep“…the United States is the Saudi Arabia of waste heat.”47
Appl.
Phys.
Lett.
91,
132102(2007). Appl. Phys. Rev. 4, 011105 (2017).
resent another opportunity for the HTE approach. The expoRecent studies have shown that for the United States alone,
nential growth of computer processing power, combined
annual potential for electrical energy recovery from waste48
with the laws of physics expressed through quantum
heat could be in the multi-terawatt range. Although thermechanics, has made it possible to design new materials
moelectric devices have significant potential to recoverfrom first principle physics using supercomputers. In thewaste heat from industrial processes, commercially available
mid-2000s, the development of high-throughput computadevices are only about 5% efficient. Therefore, discovery oftional methods and software infrastructure was pioneered
higher efficiency thermoelectric materials using HTE is critiand applied to the discovery of novel energy storage materical to enabling the practical recovery of waste heat.
Materials that exhibit a large Seebeck coefficient, highals.60,61 Importantly, HTE synthesis and measurement

High-Throughput Observation

High-Throughput Experimental (HTE) methodologies

FIG. 2. (a) Electrical conductivity, (b) Seebeck coefficient, and (c) power factor of the composition-spread (Ca1#x#ySrxLay)3Co4O9 film (0 < x < 1/3 and
0 < y < 1/3). Reproduced with permission from Appl. Phys. Lett. 91, 3 (2007). Copyright 2007 AIP Publishing LLC.56

Materials “Library” : Dispersing the several compositions on a single sheet

Observing Big Materials Space at once, Finding optimum one
13

Y. N.-Ohno, M. Haze, Y. Yoshida, K. Fukushuma, Y. Hasegawa, and M. Okada, J. Phys. Soc. Jpn. 85, 093702 (2016).

Downsampling method for quasi-particle Ohs.
J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by on 10/26/17

J. Phys. Soc. Jpn.
Downloaded from journals.jps.jp by on 10/26/17

Letters
Quasi-particle interference (QPI)
on Ag(111) surface

Y. Nakanishi-Ohno et al.
Letters

These patterns are attributed
to the aliasing eﬀect
that makes(a)(b)
Conventionalsome diﬀerent wavenumber components indistinguishable
1.4
owing to the periodicity of sampling. Figures 2(c) and 2(f )1.3
show the case of data points in a small central region
1.2
of Fig. 1(a). As shown in the ﬁgures, the QPI pattern is
1.1
1
considerably deteriorated. The inner ring disappears because
0.9
long-wavelength components can hardly be detected in the
0.8
small region. In addition, the intensity of the outer ring is
0.7
very weak because (d)the region is too distant fromdefects to
(e)
LASSO
0.6
the QPI pattern. On the basis of the above results, we
360x180 (64,800) pixelsform
conclude that CS performs well with random sampling in a(b)
broad region.
15
To quantitatively evaluate the eﬀects of LASSO, we
examine line sections of the FTs. Figures 2(g) and 2(h) show
the radially averaged line sections of Figs. 2(a) and 2(d),Fourier10
respectively. Each of the ﬁtting curves is obtained by the(g)Trans. least-squares method, where,
in the same way as Sessi
30)
et al., the following model function is employed:
Random!
Every
5" thirdX
1fðqÞ sampling¼
ak exp %ðq % "kpixels
Þ2 þ c;
ð7Þ
2
2!
k
k¼1;2
(a)

−3

0
−1
qx [nm ]

3

6

9

9
6

0

−9
−9

2.5

9

0

1

−3

6

0.5

Intensity [arb. unit]

−6

−3
−6

−6

−3

−3

0
−1
q [nm ]

3

6

9

3

6

9

−1

q [nm ]

3

6

9

y

Intensity [arb. unit]

−9
−9

2.5

9

3

1.5

0

1

−3

−6

−3

x

0
−1
q [nm ]

3

6

9

−6

2.5

2

2

1.5

1

0.5

0
−1
qx [nm ]

3

6

9

0

2.5

6

2

3
1.5
0
1
−3
0.5

−6
−9
−9

0

−6

x

2.5

−3

(f)

0.5

−9
−9

0.5

−6

2

−6

0

1
−3

0

−3

0
−1
q [nm ]

3

6

9

0

x

(h)

Using
centralregion1.5

1128 x128 (16,384) pixels
where f"1 ; "2 ; !1 ; !2 ; a1 ; a2 ; cg is a set of model parameters.

0q [nm−1]
x

−6

0
−1
qx [nm ]Intensity [arb. unit]

y

0

−9
−9

Intensity [arb. unit]

3

−3

−1

−1
qy [nm ]9q [nm−1]

1.5

−6

6

qy [nm ]

3

1.5
0

0.5

−6

2

3

y

−6

1

2

−1

−9
−9

0

−3

0.5

−6

1.5

y

1

−3

6

2

3

q [nm−1]0yIntensity [arb. unit]

q [nm−1]

1.5

Intensity [arb. unit]

3

Q space

6

2

2.5

9

q [nm ]

6

(c)

2.5

9

Intensity [arb. unit]

2.5

9Intensity [arb. unit]

R space−9−9

Y. Nakanishi-Ohno et a

Intensity [arb. unit]

J. Phys. Soc. Jpn. 85, 093702 (2016)

Intensity [arb. unit]

J. Phys. Soc. Jpn. 85, 093702 (2016)

0.5

0Fig. 1. (a) dI=dV map of Ag(111) surface. (b) FT of (a) obtained byconventional method.

0
0

2

3

4

5

6

7

8

9

0
0

1

2

3

4

5

6

7

8

9

q [nm ]
Assume without loss of generality that "1 < "2 . Note that theFig. of
2. (a)–(festimatedfrom has
diﬀerentapartsof Fig. 1(a). Subsets of data areshapecomposed whenof 7200 pixels. Results of the conventional method and LASSO
peakthe) FTs
outer
ringdelta-function-likeare shown in the top and middle rows, respectively. (a) and (d) are obtained by using randomly chosen pixels. (b) and (e) are obtained by using every thir
thepixel
phase
noise
that
arisesfromtherandomdistributionof in the r space. (g) and (h) are radially averaged lin
in both thehorizontal
and verticaldirections.
(c) and(f ) areobtained by using
only the central region
sections corresponding to (a) and (d), respectively. Each of the ﬁtting curves is composed of two Gaussian functions and a background constan
scatterers
is removed. The results show that LASSO provides(g) ! ¼ 1:3 nm , ! ¼ 3:2 nm , " ¼ 0:91 nm , " ¼ 0:18 nm , a =c ¼ 0:92, a =c ¼ 0:96. (h) ! ¼ 1:2 nm , ! ¼ 3:2 nm , " ¼ 0:84 nm , " ¼
0:15 nm higher, a =c ¼ 14,signal-to-noisea =c ¼ 18.
a muchratio (ak =c) and determines the
14 with less uncertainty (!k ). Overall, the synergy
peak locations−1

q [nm ]

−1

Q space is really sparse
Sampling in R space can be reduced.

1

pixels is 360 ! 180, namely, M ¼ 64800. There is a QPI

1

−1

−1

1

−1

2

2

1

−1

2

−1

1

2

1

−1

2

−1

1

−1

2

Basic flow chart of Application of MI
The most important

Always Start with an Issue✓ What do you want to know from data?

✓ What benefits are obtained by applying machine-learning?
✓ Can you breakdown the issue enough to solve?
✓ NEVER just USING the machine-learning.

Basic Flow
Clustering

Nα

Regression
Optimization

Precond.

Obs.

Raw Data

Low-dim.
Descriptors

ML!15

Postcond.

Analysis

Visualization

Machine-learning potentials

!16

–1.0
Low
–1

0

High
(nm)

1

Low
–0.5

0.5

High
(nm)

8

12

16

Temperature (K)Molecular dynamics simulations| Surface topographies and superconducting critical temperatures. (a,b) STM topographic images of as-deposited thin film at substrate
ure of 600 !C (a) and 400 !C (b). (c,d) STM images after post-deposition annealing (PDA) for films deposited at 400 !C (c) and 300 !C (d).
t b,c are taken with using the same film. a–c are obtained at 77 K and d is obtained at 4.2 K (all the STM images were observed at a sample-bias
f þ 300 mV and a tunnelling current is about 10 pA). Scale bar, 80 nm (a–d). (e) Growth temperature dependence of root mean square of surfaces (RRMS) values: as-grown samples (blue symbols) and after PDA (red symbols). The value of RRMS is evaluated from topographic images
at a sample-bias voltage of þ 300 mV and a tunnelling current of 10 pA (scan area of 400 nm). (f) Temperature dependence of the field-cooled
etic susceptibility for the LiTi2O4 films in a magnetic field of 50 Oe, which was applied parallel to the (111) plane. Clear diamagnetism is observed.

障壁計算
Normal modes and transition
states NEB

Structure Optimizationacbvinyl alcohol to acetaldehyde (NEB method)m~0.6 n52.14 kcal/mol

Low

High
–0.05

0.00
(nm)

Low

High

Low

–0.05 0.00
(nm)

High
–0.05

0.00
(nm)

DMol3 TS search (LST/QST) 51.473 cal/mol, -10.851 cal/mol (PBE)!19| Typical topographic images on a terrace. (a) Filled-state STM image of LiTi2O4(111) surface (11.6 nm $ 11.6 nm, sample-bias voltage Vs ofV, tunnelling current Iset of 30 pA). (b) Empty-state STM image (4 nm $ 4 nm, Vs ¼ þ 30 mV, Iset ¼ 30 pA). (c) Zoomed-up image
1.7 nm, Vs ¼ þ 30 mV, Iset ¼ 30 pA) of b. The image shows three-fold symmetry representing the spinel crystal structure. Scale bars, 2 nm (a),b) and 0.3 nm (c).

Sampling for Free energy

Determining theLTO
surfacesingle
structurevortex core by applying an external magnetic field of 1.5 T
samples range from 3.5 to 4.0 in
point-contact

using 8,perpendicular
to the surface at 4.2 K (Fig. 4a–e). At Vs ¼ " 8 mV
copy33,34 and Andreev reflection36. A recent(Nat.report
Commun.
15975 (2017))
LiTi2O4 thin films claims 2D/kBTc ¼ 4.07 from point- and þ 8 mV, we observed uniform conductance over the scanned
spectroscopy22, and thus the obtained value of 3.0 is much region (Fig. 4a,e), whereas conductance values were depressedthan that of all the previous reports. Furthermore, the around the centre of images at Vs ¼ " 4 mV and þ 4 mV
value is even smaller than that for the weak coupling limit (Fig. 4b,d). This depressed conductance is a consequence ofFrom stable structure to dynamical behaviors
ave BCS superconductivity of 3.52. We discuss later the suppressed coherence peaks. In contrast, the conductance map at
Vs ¼ 0 mV clearly represents enhanced conductance in the centre
originsMustof this unexpectedly
2D/kBTc. in nano sciences
tools forsmall
research
region (Fig. 4c). This enhanced zero-bias conductance around thecentre region is because of pair breaking. These energy evolutions
nce length on the surface. To further study the super- of conductance map indicate signatures of a vortex core (Fig. 4a–e),
ivity on the surface, we investigated the value of x from and the evolution of tunnelling spectra along line A–B in Fig. 4atronic structures around a magnetic vortex core. We first clearly shows a typical spatial evolution of spectral shape across ad the Vs dependent conductance (dI/dV) map around a vortex core (Fig. 4f).

••

We need to modelinteraction between atoms

OMMUNICATIONS | 8:15975 | DOI: 10.1038/ncomms15975 | www.nature.com/naturecommunications!17

Phys. Chem. Chem. Phys., 2018, 20, 11586-11591
3

a =0. 118, 0.424) ' for silicon and (9s5p)/[3s2p] with p
unctions (with a =0.059) ' for negative-ion states ofed in Fig. 1, in which we stretch all the Si-0 bonds keeping the Td symmetry of the cluster. We have also obtained the results for the C3„, mode, in which only one
xygen. We assume no electron orbitals around theSi-0 distance is changed, and the D2d mode, in which
oint charges.
0-Si-0 angles are varied with constant Si-0 distance.
The potential energy surface is then obtained by our
VOLUME 61, NUMBER 7
REVIEW
PHYSICALLETTERS
15 AUGUST 1988
It has turned out that these three potential energy surhanging the Si-0 distances or 0-Si-0 angles with three
ifl'erent modes. One potential energy surface is depictfaces can be fitted well by a sum of pairwise interatomic't=0.
we stretch
the Si-0
bonds keepwhichemployed
ed in Fig. 1,
a118, 0.424) for silicon and (9s5p)/[3s2p] with p
the allfunction
Wein havepotentials.
form,
the Td symmetry of the cluster. We have also obfunctions (witha =0.059) ' for negative-ion states ofing—c;cj/r',(1)
U/(r) =U;/c'""~ b(r)+f o(b;+b J)exp[( a+a J r)/(—)]
b;+bitained the results for the C3„, mode, in which only one
oxygen. We assume no electron orbitals around theParametrizingbyfitting
on Hartree-FockofSiO
44- + 4 point charge model
Si-0
which
in disperdistance
isBorn-Mayer-type
the repulsion,
point consists
charges. of Coulomb interactionchanged, and calculations
D2d mode,and
whichwith some potentials
corrections discussed
below,
0-Si-0
Si-0 distance.
is thenbetween
are varied
constant
potential Hereobtainedatoms
energyr issurface
ive Theinteraction.
by ourthe distanceand a; (b;)
is theangles
effective
radius with
(softnessof the ith
parameter)
Phys.
Rev.
Lett.
61,
86
(1988).
Si-0 distances
It has Coulombthat these withthree thewith .threeturned out
or 0-Si-0
potentialchanging
anglesmolenergy surwith the standard
tomforce fo
1 kcal A
We also include
interactions
point charges.
209(1989).
difl'erent
well bybecause
potential
surfaceis depictfaces can bein fitted
a sum ofA cautionmodes.
must One
beNature
madein339,
theenergy
evaluationof theCoulomb interaction
the cluster,
thepairwise
effectiveinteratomic
charge int•charge,
Dividing
the electrostaticthehave employedform,
he bulk, Q;, is different from that in a Si04
terms of theWecluster, Q;: In potentials.
fractionaltransferredfrom a Siinteraction to longAn, function—2/)ne andtom U/(r)to an=U;/c'""~atom b(r)+f
we a+ahave Qp=bond,have short-range
andper Si —Qs;=4/)ne, while we
Qo= —(I+An)e
range,
parts.
o(b;+b(1)J)exp[(
b;+bi )] —c;cj/r',
J r)/(—=4/)neWe.
theinteraction
Coulomb
inexpressthe cluster as a sum of long-range • andshort-range
parts asmodel, applying bulk.
Qs;
Fitting
on cluster
which consists of Coulomb interaction with some corrections discussed below, Born-Mayer-type repulsion, and disper• Reproducing SiO2 polymorphs
Uij~'"" '=Q;Q, [1 —sive interaction.
Here rg~)(r)]/r+Q;Q)gii(r)/r,is the distance between atoms and a; (b;) is the effective radius (softness parameter) of the ith
atom with the standard force fo 1 kcal A
mol. We also include Coulomb interactions
with the point charges. (2)——=[1+11((r)/8+3((r)/4+((r)in /6]exp(
gso(r)
A
caution=(I+fr)exp(
must be made 2/r),in the goo(r)evaluationof the Coulomb interaction
the cluster,2(r).
because the effective charge inthe bulk,
different in fromthat in a Si04 (theQ;: In terms of the fractional charge, An, transferred from a Si
Q;, is g&(r),Thecorrection,
the —long-range
first —
part cluster,I
atom involving
to an
atom per Si2/)neand Qs;=4/)ne, while we have Qo= —(I+An)e andbond, we have Qp=
erm)
hn 0 7)-for. the equilibrium bond length, is shown to be
Q reflects the distribution of the excess
We. express
the Coulombinteraction in the cluster as a sum of long-range
and short-range parts
Qs; =4/)neof oxygen,for which
we assume a hydrogenlike
harge
Theas featurea function of the Si-0 distance (Fig. 1).
with a'=Q;Q,
rbitalUij~'""radius [1
The radius is of the orderthat the atomic charge varies with the bond length clearI/(—here.
g~)(r)]/r+Q;Q)gii(r)/r,f the ionic radius of 0, so that we employ I/(=1. 4 A
ly indicates a many-body character of interatomic forces
(2)
'
—
Since
will
be
—
ollowing
theshort-range
Pauling.
mechanism
remaining
the
The
in
covalent
system.
physicalgso(r) =(I+fr)exp( 2/r), goo(r) =[1+11((r)/8+3((r) /4+((r) /6]exp( 2(r).
art is expected to be insensitive to the environment, we
discussed in more detail elsewhere. To concentrate on
Thethecorrection,in theg&(r),Oncesebulk Q; there.
the long-range
parameters partand first
the pair-potential approach, however, we have used cona;, b;,(the
c;
I
term)
involvingfrom
the calculation,
distribution we
excess-for. theof theis shown
to bebond results
0 7)
equilibrium
length, arehn Q;
rethe cluster
switchinsufficient
optimized
small-cluster
Since
Q reflects
stant
(Q;).
Q
for whichSi-0 distance
we assumecharge1). The featurehydrogenlike
function ofthe
the long-range(Fig. interaction,
ack intoofQoxygen,simula(i.e. , UPi""in thea bulk
andCoulomb=Q;Qi/r)
toa determine
clearwiththea radiusorbitalhere. The haveradiustheis sameof the funcorderwith Mulliken
the bond chargevarieslength itself
atomicthat the theI/( potentials
chargeon. Thus
final pairvalueof theabsolute
because'ofionicas radiusof 0,that wecharacterof interatomic
ina many-bodyonaltheform
Ida. I/(=1. 4 A
Gilbertandemploy
ly indicates
we have notforcessuggestedfunction,choice of basison theby sodepends'
Since thetheremaining
followingshort-range
in the covalentsystem.
physical mechanism
From the Pauling.
cluster calculation,
triedbecharge obtained
Instead we willprocedure.
cludedAn in thefitting The
by—is expected analysis,
to be insensitive
concentrate
part Mulliken
more detail
discussedhewhich isto the
1.7e (i.e.we
Mulliken
Qp-environment,
charge.onaround the To
of An elsewhere.
severalfixedin values,
conuse the bulk Q; there. Once the parameters a;, b;, and c;
we fitted
have used
the pair-potential
approach,paramestudiedhowever,
the trial values
here, theAmong
are optimized from the cluster calculation, we switch Q
insufficient
small-cluster
stantwith(Q;).0.Since
the bestresults
crystalareparameters.
6 reproduce
ters
Q; An
-1.
0
back-589.
into 4 Q (i.e. , UPi"" =Q;Qi/r) in the bulk simulaand
Coulomb
long-range
to determine
for a;,
Since
the fittingtheprocedurec; is nonlinear,b;, and interaction,tion. Thus the final pair potentials have Sithe same funcitself
value of theareMulliken
becauseobtained.charge
We haveset of parametersmorethan theoneabsolute
Discovery
of
new
SiO
2 high-pressure phase'Otionalform
not inbasis function,
choice I)of whichthe (Table
structuredependstheonone~ -589.
optimizes wethehave
chosen1.2 O
5- as suggested by Gilbert and Ida.
p
!18 fitting
From the cluster calculation, the charge obtained by
we tried
Insteadsimulation
procedure.
cludedAn in the ofa-quartz'
in a static
andcompressibilitye'Potential modeling by fitting

TTAM (Tsuneyuki-Tsukada-Aoki-Matsuda) Potential

0

0

0

0

I
C

~

[1] J. Chem. Phys. 103 (10), 8 (1995) .

Problem of NNP modeling
In 1995, NN fitting for interatomic potential has been reported [1]. But…

Construction of proper and simple descriptorse.g.) Descriptor based on Fourier expansion on fcc(111) surface structure[2][2] J. Behler, S. Lorenz, and K. Reuter, J. Chem. Phys. 127, 014705(2007)

HOWEVER… Such an expansion is needed to construct one by one
and is not is not applicable to a complex surfaces and.

3-dimensional coordinates are not appropriate.

• Input layer depends on the total amount of atoms• Impossible to apply for extended systems• Permutation symmetry for the same particles are broken• The other symmetries are not considered

Descriptor, symmetry are the keywords.
!19

Behler-Parrinello Ansatz
1. Introduction of Symmetry Functions describing local environments
G1i =

all
X

e

⌘(Rij

R s )2fc (Rij ),

(2)

j6=i

G2i = 21 ⇣

all
X

⇣

(1 + cos ✓ijk ) e

2
2
2
⌘(Rij
+Rik
+Rjk
)

fc (Rij )fc (Rik )fc (Rjk ),(3)

j,k6=i• Indexing the local environments based on lengths and angles• They have a invariance for translational and rotational operations

2. Dividing the total energy
energiesP H Y S I Cinto
A L atomic
REVIEW
LETTERSPRL 98, 146401 (2007)

nput and the output layer are one or more ‘‘hidden layers,’’
each witha certain numberof nodes.
All nodesfor
in theeachpermutation
symmetry• Keepingayer are connectedto the nodes in the adjacent layers bysame particle
eal-valuedEasy
weightwhich initiallyare chosen
to parameters,
extend a systemsize (or•andomly. For a given set of coordinates the output of thenumber of particles) by adding subnets
NN is then given by the expression
!
"
#$
3
2
X
X
E i ! fa2 w201 "w2j1 fa1 w10j "w1!j G!
: (1)
i
j!1

!!1

week ending
6 APRIL 2007

FIG.
!20 2. Structure of the neural network as applied in this Letterto a system containing three atoms. The Cartesian coordinates ofFlow of constructing NNP
Describing the local environments
2

1

Symmetry function vector Gk={Gk1,Gk2,…}

4
3

Atomic configurations

DFT

G1

NN1

E1

G2

NN1

E2

G3

NN1

E3

G4

NN2

E4

E

Total E• Total INPUT dimension: input of subnetworks x number of atoms• Same parameters on the subnetworks are used for the same atomic species• It is easy to expand the number of atoms by just adding the subnetwork.
!21

Total flow chart of NNP
予測

Nα
Precond.

Sampling

Raw data
DFT results（X）atomics configurations（Y）total energy (or forces)

Analysis

Descriptor

Modeling

Prediction• Symmetry Function• Structural fingerprint• SOAP• Coulomb Matrix• MBTR

Neural network
Polynomial exp.
Gaussian process
Kernel ridge
Mixture models

High-Throuput MD• Structure optimization• Materials screening• Large-scale simulation

Type of ML potentials:（１）Descriptors（２）modeling
Keeping in mind their a merit, a purpose, and a issue!22

normalized to have unit
carefully
controlledduringtraining.
!.
from the full training database (Nt ≪ Nfull ). The number ofPHYS. REV.
" B(MB)# (2017)95,094203(MB)The
kernelfunction
Both for 2b and 3b contributions, we use a squared
k qa ,qt=
prepresentative points differs for each descriptor and must be
product
exponential kernel [35]nn′ lcarefully controlled during training.
⎡⎤" (MB) (MB) #" (d )#
2(d)
toBoth for" 2b and3b contributions,we
use a squared and we findkit qadvantageous
!
=
,q
q
−q
#
t
1
a
ξ,t VOLKER
ξ,i
L. DERINGER AND GÁBOR CSÁNYI(d ) (d ) (d )⎣⎦for a sharper distinction betwe
K
q
=
exp
−
,q
,
(3)
exponential kernel
[35]
t
i
2
2 ξ
θξ
gives the final kernelmo⎡⎤" (d )#fun(d ) 2 30
DFT it advantage
and
we
find
GAP: Gaussian
Approximation
Potential
!
"
q
−q
" (d
1 the components(MB)sch
2b only (MB)(MB)
ξ,t ofξ,i) index(d ) # runningwhereKξ(d )is
an
over
the
⎣
⎦
K
q
,q
for a sharper
distinctiont an
qi ,qt (d ) = exp −, (3)a2b+3b only2
20 hasthe case2Regressionof pairs, theθdescriptor
descriptor
q . In Process
we
SOAP only
by Gaussian✓ Fittingvector
ξ
gives
the
final
kernelξ
2b+3b+SOAPThis dot productkernel is aonesinglescalarcomponent(namely,thedistancer
between12
Decomposing
2-,
3-,
many-body
terms
ap
✓
"
10
power spectrum descriptor,asthe✓two
atoms
involved):
tra(MB) (MB)
where
ξ
is
an
index
runningoverthecomponentsofthe2- and 3-body: Gaussian kernelKqa ha(up to normalization) to thein(d(2b)
)
.
In
the
case
of
pairs,
thedescriptor
has
descriptor
vectorqthe0 (4)✓ Many-body: Simple|r2 product−r1 | ≡rkernel
q =dot12 ;neighbor densities
Fu∆Energy (eV)GAP and SOAP

This dot product kern
hig
/
0
0
power spectrum descrip
/ mo
dis
/ to
ρa
d
R̂
13
23
(up
to
normalization)
0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0/
to
(2b)
permutationoftheneighbor
3 [52]:= |ratoms≡r12;(4) C–Cneighborqdistance (Å) densities
2 −r12| and2-body:⎞⎛Theforexpression
for Thisthetota
FIG. 2. Potential-energy scans
an isolated carbon dimer.
+ rthe0
13 natural coordinates,
for triplets, we do not directlyr12user
12as reference (blue), allows us to assess the use
plot, with DFT data
therefore
given
byare needed for ab
(3b)
2
⎝
⎠
of
different
structuraldescriptors:
all
three
combined
(r
−r
)
q
.
(5)
=
12
13
r13 ,3-body:and r23 , but a different form to enforce symmetryoverddriR̂
a high-quality fit (see text).
!
!
r23 2 and 3 [52]:permutation of the neighbor atomsare
E = (δ (2b) )2αt(2b) Kyie
To illustrate the role of the combined descriptors, we use⎞⎛tiNote that with this choiceof descriptors, the firstterm(andinincreasingly complex)Thetdifferent
GAPexpression
models to computefor do
Many-body:SOAP
r12 + r13!!
(m
the potential-energy
curve for an isolated carbon
dimer;
these (3bEq. (1) is equivalent
to
a
pair
potential,
andthesecondisa(3b)2
therefore
given
by
(3b)
2⎠
⎝(r
(Smoothoverlapof12atomic
position)(δand surface) training αt wa
models have been
fitted to the full+
bulk
−r
)
q
.
(5)
=
13
fin
generic three-body potential, but in the GAP framework
bothbelow that additionally incorporates DFT!
set describeddata !
tjwar23(2b)The2 results
between 0.8 and 3.7 Å in small increments.
do not impose constraints on the specificfunctional points
form.
E = (δ )one single scalar component (namely, the distance r12 between
fortwotriplets,
weinvolved):
do not directly use the natural coordinates–10 r12 ,theatoms
Descriptor
r , and r , but a different form to enforce symmetry over
241714-8

Fujikake et al.

J. Chem. Phys. 148, 241714 (!!
are summarized in Fig. 2: GAP models using 2b descriptors
c.f.：Li
Diﬀusioningraphites(MB)2
For the many-body term, we use the recently introduced
only, or a combination of 2b+3b, reproduce the minimum
+ (δ

)

ite

(M
α
t ab
t

i
Note that with this choice of descriptors, the andfirstterm intherepulsion
at
small
C-C
distances
reasonably
well,
but
Fujikake
et al.,
J. Chem.
148,241714 (2018).
smoothoverlapof Phys.
atomicpositions(SOAP) [53] descriptor,
!
!
a
t
wa
longer-range
is not yet correctly described.
Eq. (1) is equivalent to a pair potential, and the thesecond
is behaviora(3b) 2 Anthe+ (δ descriptor)
interestingfor
result is seen when using a many-bodywhich has proven successful in generating!23 GAP models
liq
where
allregion
fittingcoefficientsαgeneric three-body potential, but in the GAP framework
only: the fit both
is very goodfor thewheredata points

FIG. 6. Molecular dynamics simulations of Li diffusion in a graphite-like framework at 1000 K. Results of GAP-MD have been structurally benchm
against DFT-MD data. (a) Exemplary GAP-MD trajectory visualized by plotting coordinates of Li atoms at equally spaced time intervals as purple dots, w

eter Q6 [36], is additionally reported in the Supplemental
Phys.
Lett.
120, 143001(2018).
Material andshowsRev.
excellentagreement
between DPMD
and AIMD trajectories. In the case of the molecules, weFI
perform DPMD at the same temperature of the originalec
data, using a Langevin thermostat with a damping timelin
τ ¼ 0.1 ps. The corresponding distributions of interatomic
distances are very close to the original data (Fig. 4).
PHYSICAL REVIEW LETTERS
120,and143001(2018)
Scalability
computationalcost.—All the physicalsoquantities in DPMD are sums of local contributions. Thus,th
TABLE III. The equilibrium energy and density, Ē and ρ̄, of[modeling]
after training on a relatively small system, DPMD can bewater and ices, with DPMD and AIMD. The numbers in squareacbrackets are the AIMD results. The numbers in parenthesesdirectlyareapplied to much larger systems. The computaco
Framework of Behler-Parrinello
NNP
statistical uncertainties in the last one or two digits. The training
tional cost of DPMD scales linearly with the number ofsm
AIMD trajectories for the ices are shorter and more correlated
atoms. Moreover, DPMD can be easily parallelized due tothan in the water case.
ev
its local decomposition and the near-neighbor dependenceth[Descriptors]
SystemĒ (eV=H2 O)
ρ̄ (g=m3 )of its atomic energies. In Fig. 5, we compare the cost ofin
Liquid (invariance−467.678ð2Þ [−467.679ð6Þ]
1.013(5) [1.013(20)]・Inner coordinatesforrotation)
DPMD fixed-cell simulations (NVT) of liquid water withmwater
that of equivalent simulations with AIMD and the empiriIce Ih (b) −467.750ð1Þ [−467.747ð4Þ] 0.967(1) [0.966(6)]Ex・Based on Inverse
Ice Ih distance(c) −468.0478ð3Þ [−468.0557ð16Þ] 0.950(1) [0.949(2)]cal FF TIP3P (transferable intermolecular potential with 3
in
TERS 120, 143001 (2018) Ice Ih (d) −468.0942ð2Þ [−468.1026ð9Þ] 0.986(1) [0.985(2)]points) [41] in units of CPU core seconds/step/molecule.
While in principle the environmental dependence of Ei
atomj
is analytical, in our implementation, discontinuities are
RDFs of ice are reported in the Supplemental Material.
A higher-order correlation function, the probability distripresent in the forces, due to adoption of a sharp cutoff
bution function of the O─O bond orientation order paramradius, limitation of angular information to a fixed number
Rijeter Q6 [36], is additionally reported in the Supplementalez
of atoms, and abrupt changes in the atomic lists due toMaterial and shows excellent agreement between DPMD

DeePMDand AIMD trajectories. In the case of the molecules, weof the organic molzij DPMD at the same temperature of the original3.5 FIG. 4. Interatomic distance distributions
ey
perform
0.8
ecules. The solid lines denote the DPMD results. The dashed
DPMD O−O
DPMD
data, using a Langevin thermostat with a damping timelines denote theDPMD
AIMD
results.
0.7
O−H
DFT
3.0
ex
τ
¼
0.1
ps.
The
corresponding
distributionsofinteratomic
DPMD
H−H
yij
0.6
DFT O−O
2.5
distances are very close to the original data (Fig. 4).
DFT O−H
0.5
DFT H−H
Scalability
and computational cost.—All the physical2.0 sorting. These discontinuities
are similar in magnitude toxijquantities in DPMD are sums
or of local contributions. Thus,1.5 those present in the AIMD forces0.4due to finite numericalafter training on a relatively small system, DPMD can be
accuracy in the enforcement of 0.3
the Born-Oppenheimer
directly applied to much larger systems. The computa-1.0 condition. In both cases, the discontinuities
0.2
are muchtional cost of DPMD scales linearly with the number of0.5 smaller than thermal fluctuations0.1and perfect canonicalatoms. Moreover, DPMD can be easily parallelized due toevolution is achieved by coupling
the systems to a
0.0
FIG. 1. Schematic plot of theits local
neuralnetwork andinput
for the dependence0.0 0 thermostat.
decomposition
the near-neighbor
1
2
3
4
5
6
0.5
1
1.5 Coulomb
2
2.5
3
We further note that long-rangeof its atomic
In Fig.
5, wethe cost ofr [Å] not treated explicitly in the current
ψ [rad]environment of atom i, taking water
as an energies.
example.
Atomj iscompareainteractions are
impleDPMD fixed-cell simulations (NVT) of liquid water withmentation, although implicitly present in the training data.
generic neighbor of atom i, ðex; thatey; eofz Þequivalentis the local
frame of atom
simulations with AIMD and the empiriFIG. 3.Explicit
Correlationof liquidwaterfrom DPMD
treatmentfunctionsof Coulombic
effects may
be necessary!24i, ex is along the O─H bond, ezcal
is FF
perpendicular
to theplane of potential
TIP3P (transferable
intermolecularwith3
in some applications
and(Right)
deservesThefurtherstudy.ADF within aand PI-AIMD.
(Left) RDFs.
O-O-O

It works well in molecular systems.

P(ψ)

RDF g(r)atom i

FI
D
pe
CP
co

Diﬀusion pathways in amorphous
Problem：Too expensive to calculate with DFT（N3 times optimization x simulation time T) N~50, T > 1 h => 13 year
NEB calculation is also necessary to estimate the activation energy…

Etot = Eamorph +ECu +Cu insertion term

Eopt

Relaxation term

Construction with Behler-Parrinello ansatz✓ Considering amorphous matrix as “a field”✓ Learning only energy of Cu atom with its surroundings✓ Able to simplify the NN even for ternary systems✓ The calculation cost is cheaper than that of full NN✓ Impossible to execute the molecular dynamics!25

Found diffusion pathways

upercells with a vacancy have been generated by removing
Li atom and then relaxed the structure. Since the supercell
ontains 12 Li atoms that are not equivalent to one another,
we can obtain 12 different Li vacancy formation energies.
igure 4 shows the comparison of the vacancy formation

tion of Li vacancy. In reality, a vacancy-interstitial pair would
be formed rather than an isolated vacancy. Since the purposeof this section is the verification of the applicability of the NN
potential in various calculations, we do not discuss the caseswith complex defects here.

Small system to Large system
C. Diffusion paths and networkNN

IG. 3. The final structure of the amorphous Li3 PO4 model (Li12 P4 O16 ) after
elt-quenching simulation. The structure was fully optimized with DFT.

DFT for small cell

Here we examine the Li vacancy diffusion network(including the Li migration paths, barrier energies, and the
topology connection of paths) using the NEB method. Thismethod is frequently used to find migration paths and theirenergies. In this method, a number of intermediate images,
which are constrained by imaginary “spring” force, are interposed between the starting and ending equilibrium positionsof the target atom. By optimizing the atomic coordinates ineach image according to the forces given by both imaginary spring and real potential energy surface, we can find theminimum energy path of atom migration between the two equilibrium sites. The optimization can be based on the potential
energy surface provided by the DFT calculation or classicalinteratomic potential (e.g., NN potential).
In the present study, we assumed that a Li atom candirectly
hop to the neighborwhen the distance
Simulation
forvacancy
largesitecell

Small, accurate
training data-sets

Scale up with ML potential

Reproduce the RDF and Diffusion constant estimated by DFT!26

W. Li, and Y. Ando, Phys, Chem. Chem. Phys. accepted.

Vibration property depending on structure
Amorphous model

Phonon density of state

The phonon density of state (DOS) calculated from ML force fields.

1728 atoms[Science, 1997, 275, 1925–1927; Europhys. Lett., 2002, 60, 269–275; Phys. Rev.
Lett., 1985, 54, 441–443]

It is possible to anneal slowly for large-scale system by ML potential（0.01K/fs）
Improving ring statistics of Si-O networks/Phonon DOS also agrees well.
!27

High-throughput
Spectrum analysis!28Collaborated with 永村直佳（NIMS）

Automatic estimation of peak position
Experimental data-sets of spectrum

By hand

D

S


x

Peak estimation
by EM algorithm

X-ray spectrum imaging
29

Collaborated with 永村直佳（NIMS）, 松村太郎次郎（PD）, 永田賢二(AIST), 赤穂昭太郎(AIST)

Auto-estimation of peak position
Issue：”Automatic fitting to finding peaks in many spectra”
14

E
x
y
S

301
302
303
304


305

D
0E6
B /.4
e
r

26 6
m

E6 G 6(e

VVaf

X-ray Spectromicroscopy

c

b

m

x

Need to extract peak positions of whole spectra
Fitting by gaussian mixture model?
30

Character mappingDifficulty of parameter fitting
Hard to use non-linear fitting scheme.
non-linear LS fitting= searching better initial guess

Searching based on their experience

Handling each data manually

Impossible to analyze big-data

y = A sin(Bx + C) + D

解決策
Eﬀective way to find（EM algorithm）
Stochastic sampling（monte-carlo method）!31

Even though there are no noise,
It is not work with bad initial guess.

Maximum likelihood approach
ML estimation For Gaussian distributionp({x1 , · · · , xN }) =

N
Y

n=1

N (xn |µ,

2

)=

✓

p

1
2⇡ 2

◆N

exp

1

N
X

2 2 n=1

(xn

µ)2

!

Likelihood: probability of obtaining observed data
Maximizing it

ML estimation for Gaussian mixture modelp(x) =

K∑
k=1• Mixing ration(sum. 1)• Summation are in log.
• Diﬃcult to solve analyticallyπk N(x | μk, σk)ln p({x1, ⋯, xN}) =

N

∑
n=1

ln

K

πk N(xn | μk, σk)
∑
(
)
k=1

32

Collaborated with 永村直佳（NIMS）, 松村太郎次郎（PD）, 永田賢二(AIST), 赤穂昭太郎(AIST)

Expectation-Minimization algorithm
Gaussian mixture model(1)

Estimation
1. mean μk , variance σ2k, mixing ration πk
2. Latent variable rnk for data n(2)

E (Expectation)-step：Estimate (2) from (1).
expectation of rnk（: responsibility γ(znk)）

πk N(xn, θk)
E[znk] == γ(znk)∑j πk N(xn | θk)Depend on (1) and xn

図３­１ 混合ガウス分布の例

Nk =M (Maximization)-step：Estimate (1) from (2).

N

∑
n=1

γ(znk)

Maximizing expectation of log. Likelihood from complete data（pairs of xn and rnk）

μknew =N1γ(znk)xn,
Nk ∑n=1

σk2 new =N

1γ(znk)(xn − μknew)2,
Nk ∑
n=1
33

πknew =

Nk
N

Likelihood monotonically increase!

Spectrum modeling of EM algorithm
Conventional EM algorithm: the cost depends total data numberTime [s]

0.8

481.0

RMSE

8707.4

933.8⇦

improved！
.

19748.5
9PN (dszr
66086.0c

B U U LK SL OUK c

A KPTH 8 HRNU P OS
cost.
• Robust for selection of initial guess, cheap computational{ m d• Multi-trial of initialization makes improve the accuracy.
• Noisy and peak overlapping case is diﬃcult.

現在はこの手法を拡張して「ピーク本数」「pseudo Voigt関数による自動フィット」まで可能
34

Key consideration for using Informatics“Finding” does not mean “Understanding”✓ Many researchers have Internet of “Mechanism”.
✓ Machine-learning does not take into account the “physical law” generally

Analysis and Interpretation is must to do.
Trade of “Accuracy” and “White box”
Non-linear formalism

Linear formalismY =X

ai Xi + bi
NOT so high-accurate
BUT WHITE BOX

High-Accurate one
BUT BLACK BOX
35

Several problem on Materials Informatics“Small data” rather than “Big data”✓ Conventional Experiment: less than 100, maximally 1,000.
✓ Systemic error depending on labs. (human, equipment, etc.)✓ Spectra and Simulation data is close to “Big data”.

Sequencer

NO universal representation✓ object: polymer・semiconductor・metal etc.
✓ Scale: 1 nm 〜 1 m（9 order!!!）

Genomics✓ Diversity of representation for systems and scales

Key aspects ofMaterials sciences“controllability (reproducibility)”and ”theory (prior knowledge)”

Prior cases are important36

Always starts with an issue
Not necessary to be an expert of ML.
✓ Algorithm development is too diﬃcult for materials researchers.
✓ Basic knowledge helps us to follow the cutting-edge algorithms.
✓

Important things is communication with experts.

Making a issue in materials science IS YOUR WORK
✓ ML experts cannot make an issue from the aspect of Materials science.
✓ Our duty is making an issue by ourselves.
✓ Applying basic algorithm by ourselves initially if possible.

Understanding basics of ML makes us to communicate with experts.
Communication makes us to solve our issue!
37

Take-home messages
The most important

Always Start with an Issue✓ What do you want to know from data?

✓ What benefits are obtained by applying machine-learning?
✓ Can you breakdown the issue enough to solve?
✓ NEVER just USING the machine-learning.

Basic Flow
Clustering

Nα

Regression
Optimization

Precond.

Obs.

Raw Data

Low-dim.
Descriptors

ML!38

Postcond.

Analysis

Visualization