対称性と深層学習にみる
ポテンシャルフィッティングの近代史
産業技術総合研究所 研究員
安藤 康伸

1

–1.0
Low
–1

0

High
(nm)

1

Low
–0.5

0.5

High
(nm)

8

12

16

Temperature (K)ure 1 | Surface topographies and superconducting critical temperatures. (a,b) STM topographic images of as-deposited thin film at substratemperature 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).
te that 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-biastage of þ 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 surface
ghness (RRMS) values: as-grown samples (blue symbols) and after PDA (red symbols). The value of RRMS is evaluated from topographic images
served 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-cooledmagnetic susceptibility for the LiTi2O4 films in a magnetic field of 50 Oe, which was applied parallel to the (111) plane. Clear diamagnetism is observed.

原子・分子シミュレーションの応用例

安定構造探索振動モード・反応座標

a

b

c

障壁計算 NEB

vinyl 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

ure 2 | 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 of
900 mV, 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
7 nm $ 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),nm (b) and 0.3 nm (c).

統計サンプリング(自由エネルギー・エントロピー)

LTO表面構造の同定

Ti2O4 samples range from 3.5 to 4.0 in point-contact single vortex core by applying an external magnetic field of 1.5 T(Nat.report
Commun.
15975 (2017))
using 8,perpendicular
to the surface at 4.2 K (Fig. 4a–e). At Vs ¼ " 8 mV
ectroscopy33,34 and Andreev reflection36. A recent
itaxial LiTi2O4 thin films claims 2D/kBTc ¼ 4.07 from point- and þ 8 mV, we observed uniform conductance over the scanned
ntact spectroscopy22, and thus the obtained value of 3.0 is much region (Fig. 4a,e), whereas conductance values were depressedmaller than that of all the previous reports. Furthermore, the around the centre of images at Vs ¼ " 4 mV and þ 4 mV
esent value is even smaller than that for the weak coupling limit (Fig. 4b,d). This depressed conductance is a consequence ofs-wave BCS superconductivity of 3.52. We discuss later the suppressed coherence peaks. In contrast, the conductance map at
ssible origins of thisunexpectedly small 2D/kBTc.
Vs ¼ 0 mV clearly represents enhanced conductance in the centre
静的な安定状態から動的な挙動まで
region (Fig. 4c). This enhanced zero-bias conductance around thecentre region is because of pair breaking. These energy evolutions
構造と特性・観測データを結びつけるキー
oherence length on the surface. To further study the super- of conductance map indicate signatures of a vortex core (Fig. 4a–e),nductivity on theナノスケール物理では必須のツール
surface, we investigated the value of x from and the evolution of tunnelling spectra along line A–B in Fig. 4ae electronic structures around a magnetic vortex core. We first clearly shows a typical spatial evolution of spectral shape across aalysed the Vs dependent conductance (dI/dV) map around a vortex core (Fig. 4f).

•••

TURE COMMUNICATIONS | 8:15975 | DOI: 10.1038/ncomms15975 | www.nature.com/naturecommunications!2

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

by

2)
isaccord
discovered
the following:
full
scalar distance is obvious; similarly, f3 must possess (2.
withthe
These lattertwo
reasonable
criteria
experiment.
here/(r) isforce
betweenmolecules when se
r.
rk!o)
/o,
rk)=Ef3(r;v3(r;,~/o,rj,
translationalardiamondbe our point of view in the folstructurewas
indeedThat will symmetry.
themosttials, v2 and andstableperiodicv3. rotational
iniserved
tocontribute
to
our
of
severalofourrejection=7.
3
049 556277, B =0.602 224 558 4,
selectionofbetweenIn theInchosenbehindof obviously
as thatwe
rangementlow
the context,same(1.1), toparticles atsets.
spiritleast
lowing.
pressure,toat——1,interactions
tial
choicesparameterand forisrelated(r) amongby the equation •and o Eq.
is chosen
e ispresent
where
to giveq depth
o.
=1
task
than
that
more
Si
atoms
is
considerably
demandingalternatives.
a
a
first
units
and
then
lattice-sum(2.cric
and
after
this
introduceandsimpleHowever,7) /•Q
length
80,
p=4,
q=0,
energy
,
the
most
set
thus
far
Overall,satisfactory
parameterO
make f2(2'~ ) vanish. That f2 must be a function only of
dissmalltodescribe
confronted
when
terion
was
only
it
was
still
attempting
satisfied,that
thenecessarywrite
melting
is the following:
discovered
scalar distance is obvious; similarly, f3 must possess full
=21.0, liquid
20 .
r) = —in- I
y =1.structureThe
tetrahedral geometry.
fromthe ideal symmetry.
inferred for the n(placements
tentative
pointA, and thetranslationalandrotational
Downloadedfrom
http://rspa.royalsocietypublishing.org/
on May22, 2018=Efq(r~/cr),
)
v2(rj
=7.
=0.
3
049
B
602
224
558
556277,Jr
4,and phonondescribe selection
elastic properties
lattersuffices to context,
molecular-dynamics
teractionthesimulation
also
in
be
by
is
which
all
results
inThis
thechoicethereported
upon
interactionsofbetweenIn thepresent(2.
re-2)
restricted
In such atask
ofisthecrystalline material.
spectrawith
These latter criteria
reasonable
experiment.
remainder
ofaccord
this paper
based.
are
than
that
more
Si atoms
considerably
demanding
a
a =1
(2.7)
80,
p=4,
q=0,
r.
rk!o) potentials of the/o, ~/o,v3(r;,
rj, rk)=Ef3(r;463Theidentityofthetwoformulaeforutilizetoitbeappropriate
gime
may
ini- B is easily estab
servedtocontribute
to
our
of
severalofourrejection
1
shows
the
corresponding
binding
Figureenergyperdissmalltodescribe
confronted whenonlyattempting
'
is
to
=1.
it
=21.
however,
20 . units)
span
Keating—1,mandatory
fortial A,choices
sets. versus number 二体バネモデルなど
0&/X 0,reducedatom(in yparameter
density p
ando is chosen
where e type.
is chosen
toTheto
idealq depth
tetrahedral
from Now,
thegiveplacements
geometry.
piation(2.03)
byparts.
The
proof
requires
that
in
order
achievetoavailable
configuration
theentire
space
the
most
set
thus Alfar
Overall,satisfactory
parameterin reduced units) for several simple lattices.
f2(2'~ ) to(also
Thatelastic
be a function
vanish.
make suffices
only ofand phonondescribe
latterproperties
f2 must
is
which
all
results
inThis
thechoicethereportedupon
the short-range
orderin theofPotential
description
satisfactory
ascalaris thepressure
discovered
following:zero
diamond structure (DIA) is
thebased.
at of
full
is obvious;material.
thoughpossess
distance
similarly,
reIn f3
suchmust
a restricted
Lennard-Jonesof the crystalline
spectraremainder
this paper are.
r3body(e*j”{r)—1)
diffusive
and
of
the
atom-exchanging
On
the
DeterminationofMolecular
Fields.—II.
From =
the0,liquidphase,
cubicLtother cases=0.simple
thepreferred,translationalrotational symmetry.
(sc),of theutilize
toSoc.
it may andbe appropriate
potentials
gime
=7.
1049
shows
the
corresponding
binding3
B
602
224
Figure558
556277,energy
per
4,
Proc.
R.
Lond.
A,
106
463(1924).
in
occur
the
motionsthat
continuouslyr00liquidphase.
' Now,
cubic
arecentered
selectionit ofbetween(fcc)
In the type.
(bcc),
context,
present
is interactions
however,mandatory to
span
Keating
oface-centeredf Statef cubic
a Gas.
reducedand units)
versus onumber
atom 0&/X(inEquation
density p
✓
✓
◆◆
in
for
our
was
use
The
reducedpotential
study
pair
⇣ more⌘12demanding=1for80,severalis importasktothanthatclose.
We
that thissimple
closenessreasonably
Si atoms
a considerably
a believe(2.
q=0,
achieveavailable
configuration
theentire isspace in 6orderAl-7)
in reduced
units)
lattices.
(also p=4,
hichiseasilyseen
tobeaconditionthat
B shall remafrom when
the following
five-parameterfamily:selected
By
J.
E.
J
o
n
e
s
,
D
.S
c
.,
1851
Exhibition
SeniorResearch
Student,
discolsmall
tant
in
that
to
a
to
describe
thermodynamic
melting
confronted
assuring
only
attempting
a satisfactory
ULJ (r) =description
4✏ . of the short-range order in theisdiamondstructuretheat zeropressure(DIA)
thoughA, =21.
=1.
200,structurey Trinitybe. College,
temTheideal
tetrahedral r geometry.
at cubica reasonable bodypossible Cambridge.
placements
lapsed liquidthe
r atom-exchanging
diffusive
and thetheliquidphase, fromother oncases~ of r~)exp[(rsimplepreferred,
restriction
themolecularmodelswhich are possible
—
(sc),
—
A
(Br
r
a
a)
~
and
.
describeelastic
interaclatter
sufficestophononthesumofthree-bodyproperties
construction,
perature.
Bymotions
that occur continuously in the liquid phase. (2.3)
is
which
all
results
inThis
thechoicethereported
f2(r)=
upon
centered cubic (bcc), and face-centered cubic (fcc) are
reAlIn suchmaterial.
a restrictedof the crystalline
for
diamond structure.
tions
spectranearly of vanishes
in ourfor use
The reduced
study wasremainder
based.
are thepair potential
paper(CommunicatedbythisProf.
S.believe
Chapman,
F.R.S.—Received
April 22, 1924.)
close.
We
that thiscloseness is imporreasonably
latto utilize potentials
sum substantiallydestabilizes
the other colthat1 shows
be following
thoughappropriatefive-parametergime it frommaythefamily: of theselected
the thermodynamic
corresponding
binding
気体粘性の温度依存性を「斥力」「引力」相互作用で説明する事が目的
tantFigure
in assuring
thatto a permelting energy
•
'
auare
where
and
a
This
form
positive.
generic
3,
B,p,alternativecompensate
tices,
partiallynumber
Keating type.
Now, however, it is mandatory to span
byp
0&/X
reducedstructures
atom these
density
temstructure
beunits)
liquid (inlapsedC
O possible
N T Eversus
N TatS a. reasonable
I" =a without discontinuities
気体の状態方程式の２次のビリアル展開を解析的に議論
off
in anytomaticallycuts
at• available
~
—
in
order
achieve
bonds
to
—
more
of
than
does
configuration
thestabilizing
the entire
space
incorporating
pair
r ~a
r~)exp[(r a)Alunits) fortheseveral
lattices.
(also in reduced. A (Br
Page
interacsum ofsimple
three-bodyconstruction,
perature.
By lattice.
ra satisfactory
which is aof distinct
derivatives,
advantageorderin in(2.
f2(r)=
any
the short-range
thedescription
まだ分子の内部構造は十分分かっていない時代のこと
diamonditisachieve3)
theimportanttoEvidently,
•
is
zero
diamondstructuretheatpressure(DIA)
thoughlie
conditionrequiresthat
neither
n
nor
m
shall
Alvanishes
forthediamondstructure.
tions
11.
Introduction..........................................................................................................................
463 have anearly
本質をつき且つ解析的に扱える点がすごい
application.
molecular-dynamics
diffusiveof . the atom-exchanging
in orderbetweenthese competing
effects
apreferred,balanceliquid phase, and simulation
proper that
cubic
other
bodythecases
simple(sc),sum substantially destabilizes the other latthough Virial Coefficient...........................................................................................
2. The Second
465
The same
cutoff continuously
extended
to j theadvantage can thebeliquidoccur
motionsthat
5264
phase.
FRANK
modelThe
H. STILLINGER
tocenteredSicubic
successfully.
AND THOMAS
difficultyofexpressingBasan
explicit
functionwhere
and a are positive.in This
form auface-centered
cubic
and
are
generic(fcc)
3,(bcc),
B,p,
alternativethese
structurescompensate
tices,
partially
byJwas
3. The Determinationof the Force
Constants.................................................................
468
three-bodyinteractions
them
weourassign
Specifically,
f3.
in
for
use
The
reducedpotentialstudy
theminimuminthediamond
of
curve
At
pair
Fig.
1,
I"
is importomatically cuts off at =a without discontinuities in any
Webelievethat thispairclosenessreasonably close.
bonds thanmoreof thedoes
stabilizing
incorporating
arrangeF. Stillinger 5
and
A. Weber,Phys. Rev.
B 31,
5262
(1985).
4. T.The
Molecular
Fieldof
Argon...........................................................................................
469 talline
form from the following five-parametertheStllinger-Weberfamily:
0(ly
of
this
paper.
The
later
valueshavebeen treatedindepecolrselectedwhich is a Potential
tant
in
that
to
a
derivatives,
distinct advantage
in anythermodynamic
meltingassuring
menta.
The—
lattice.
itisachievethe pdiamond
importantto=Evidently,
460o.
=0.
1.999 993m. .
(2.8)
, @/X
lead
to
explictemstructure
at a reasonableapplication.
bethesepossiblemolecular-dynamics
liquid= h simulationh ( r~;,
HJ' k ) +J'k )f thispaper,
leading
toresults
more ininorderline w ith those enter
already
betweeneffects
alapsedbalance3 ( r;, rj', rk )competingproperthe New~ r;k—berjk extended—(Br ( r;~',r ~ a to ther~)exp[(rinteracthe
sumofthree-bodyconstruction,
perature.
Thef2(r)=same . A
cutoffcan a)
advantage
1.
Introduction.
integrated usi
By
model
to
Si
successfully.
(2.4)3)
leasurements(seeforAppendix).
mental time(2.
6;z~Al0, +h (r„;,vanishes
the diamondtions
nearlythree-body interactions
we assign them
f3.r&~,Specifically,theminimuminthe diamondofstructure.
curvehasAtFig.
1, that the temperature1s
The
investigationofaprecedingpaper*
shownF
sum substantiallythe cit.
other latKeesom,
Core, J. E.destabilizes
Jones,
loc.
though * thatthe form
subtended
and rkat• Siなど共有結合性結晶でのstandard.
is theandbetweenwhereangle
r=o(m/sauwhere OJ;ka areformvariation
positive.rj This
generic
3,
—
B,
=
p,
460o.
=0.
1.
993m.
.
@/X
999
(2.
8)
alternativethese
structures Theory
tices,
partially ofcompensate
by, Dynamicalofviscosity,asdetermined
experimentally,canbesatisfactorilyp
two-parameter
ThefunctiontoatJeans,‘
Gases,’
3rd
edn.,
p.
132,
eqns.
vertexetc.
hbelongsi,
Under these=aHJ' withouth ( atdiscontinuities
in any• ３体項を加えて角度成分を表現
tomaticallyr;~', I"r;kJ'k )rk ) =offrjkk ) + h ( r~;,
3 ( r;, rj', cuts
incorporating more of the stabilizing pair bonds than does
constant to a
less
(A, , y&0). Provided that both r,j and r;k are explained
family
in
many
gases
on
the
assumption
that
therepulsive
andattractiver derivatives, which is a distinct advantage in any• カットオフ半径の導入
diamondlattice.
itisachievetheimportanttotypically of s
Evidently,(2.4)
introduced
cutoffthan
the previously +h(r„;,r&~, 6;z~a, it has the followed by the usu
simulation
application.
molecular-dynamics
parts of athemolecular
field
are
each
accordingtoaninverse
power
of
the
between these competing
effects in orderproper balance
mon factor.
fcc
form:same cutoff advantage can be extended to the• シンプルな記述
bcc
ing TheN8- — 2r/0
atwhere OJ;k is the angle between rj and rk subtended distance.
to modelSi successfully.
In somecases, in argon, for example, it was further
shown that the During som
cussed below
three-body
2. Reduced—a) '+y(rk—assignpair potential
f3. h Specifically,
two-parameter
function
to a wevertexbelongs
in the diamondFIG.curve
ofbeyondh (rj,a) '] themi, etc.
At the minimumFig.
1,1.80. vs distance. This function
rk,interactions;k) =kexp[y(»;,
8,The
—1.4
rl(t), . . . , r~vanishes
identicallyr =molecularexperimental
facts
can
be
explained
by
more
than
one
model,
from
form
the
family (A, , y&0). Provided that both r,j and r;k are lessuration rlq, .
(2.
5)
1.999 993m. .
(2.8)
X ( cosign'y +cutoff
did this mapp
, @/X = —7 )h ( a, it has )the followp =0.460o.
which
we inferredthatthan the previously= hintroduced—1.6 viscosity results alone are insufficient to determine pre~
数値計算に向いた設計がされている
HJ'
(
)
(
)
+
rj',
r;~',
J'k
rk
r;,
r;k
ics; it was sir~;, rjk
k3
bcc toofcorresponding form: h vanishes identically. The "ideal" tetrahedral
2
fororder fcc
theseto the observed
cisely the nature ofN8- —molecularfields. In The
objectthe present
paperlattice
is to The mapping
otherwise
method to-1.
andatomization
6
spacingenergy of crystalline Si at 0 K, it
Q
(2.
4)
+h
(r„;,
r&~,
6;z~
such ;k)angleh 0,
steepest-descen
wouldchoosea molecular
model
of
the
same
type
will
also
explainavail
benecessary
to!3=kexp[y(»;, —a) '+y(rk —a) '] ascertain whether(rj,isrk,
8, that—1.4
subtended atandis thebetweenwheref

原子間ポテンシャル
f

—

r)a

—

'],

0, r)af—'],0, r)a

—

—

),

—

'],

0.0

f

),

—1.5

&.

f

),

&.

0.6

0.8

I

I

1.2

1.4

I

1.8

1.8

8.0

COMPUTER SIMULATION OF LOCAL ORDER IN. . . SILICON

31

tials, v2 and v3. That will be our point of view in the following. In the same spirit as that behind Eq. (1.1), we
first introduce energy and length units c and o., and then
write

5263

diamond structure was indeed the most stable periodic arrangement of particles at low pressure, at least among
simple alternatives.
However, after this lattice-sum criterion was satisfied, it was still necessary that the melting
point and the liquid structure inferred for the tentative interaction by the molecular-dynamics
simulation also be inreasonable accord with experiment. These latter criteria
served to contribute to our rejection of several of our initial choices for parameter sets.
Overall, the most satisfactory parameter set thus far
discovered is the following:

ポテンシャルのパラメータ決定
v2(rj ) =Efq(r~/cr),v3(r;, rj, rk)=Ef3(r; /o, r.~/o, rk!o)Stllinger-Weber Potential

f

(2.2)

—1, andRev.oB is31,chosentowhere e is chosen
to giveand qT.depth
F. Stillinger
A. Weber, Phys.
5262 (1985).
make f2(2'~ ) vanish. That f2 must be a function only of
満たすべき条件（物性・特性）
scalardistance is obvious; similarly, f3 must possess fulltranslational and rotational symmetry.
3 =7.049 556277, B =0.602 224 558
低温低圧でdiamond構造を取るように調整
In the present context, selection of interactions between
Si atoms is a considerably more demanding task than thata =1 80,(2.7)
融点や高温での液体構造を再現するように
confronted when attempting only to describe small dis20 .
A, =21.
y
調整from the ideal tetrahedral geometry. The
placements
latter suffices to describe elastic properties and phonon
This is the choice upon which all results reported in the
spectra of the crystalline material. In such a restricted reremainder of this paper are based.
gime it may be appropriate to utilize potentials of the
Figure 1 shows the corresponding binding energy per
Keating type. ' Now, however, it is mandatory to span
atom 0&/X (in reduced units) versus number density p
the entire available configuration space in order to achieve
several(also in reduced units) forLJ
simple lattices. Alpotential
ポテンシャルの物理モデルを理論的に作成
a satisfactory description of the short-range order in thethough at zero pressure the diamond structure (DIA) is
diffusiveliquid phase, and of . the atom-exchanging
SWsimplepotential
cubic (sc), bodypreferred, the other cases
motions that occur continuously in the liquid
phase.
face-centeredmodel
cubic (fcc) are
分子動力学の実行
centered cubic (bcc), and
Gay-Berne
The reduced pair potential for use in our study wasreasonably close. We believe that this closeness is imporパラメータ空間の探索
from the following five-parameter family:
selectedetc. melting to a coltant in assuring that thermodynamic

4,

•
•

p=4, q=0,0, =1.

ポテンシャル構成の流れ——

be possible at a reasonable temthe sum of three-body interacモデル選択と多目的最適化を
for the diamond structure. Altions nearly vanishes
though that sum substantially destabilizes the other lat満たすべき条件１を満たすパラメータを収集
where 3, B, p, and a are positive. This generic form au専門家が独自に行なっている
structures partially compensate bytices, these alternative
tomatically cuts off at I" =a without discontinuities in any
incorporating more of the stabilizing pair bonds than doesr derivatives, which is a distinct advantage in any
ようなもので、パラメータ発
the diamond lattice.
Evidently, it is important to achieve
simulation application.
molecular-dynamics
a proper balance between these competing effects in order
別の条件も満たすベストなパラメータセットを選択
The same cutoff advantage can be extended to the
見も非常に高難度（逆問題）
to model Si successfully.
three-body interactions f3. Specifically, we assign themAt the minimum in the diamond curve of Fig. 1,
the form!4=0.460o.
@/X = —1.999 993m. .
(2.8)
A (Br
f2(r)= .

~0, r)a=lapsed liquid

structure—'], r ~ ar~)exp[(r —a) 物性・構造を計算
perature. By construction,
(2.3)

p

,

逆問題からフィッティングへ
第一原理計算A first principle is a basic, foundational, self-evident proposition or assumptionthat cannot be deduced from any other proposition or assumption. (@Wikipedia)• 基礎方程式と基礎物理量に基づき実験事実に合わせ込んだ経験的パラメータを用いない手法• 実験結果を仮定せず扱いたい系を支配する物理モデルのうち最も基本的なもののみで物理量を計算する手法• 密度汎関数法（DFT）がもっとも主流（many-body Schrödinger eq.を有効一体問題にマップして解く）• 固体振動のバネモデルはバネ定数を実験値から決めなければいけないので第一原理的ではない• DFTからバネ定数を求めて振動モードを計算する方法は第一原理計算• ハートリー・フォック法, QCDなども第一原理計算

ポテンシャル設計の流れ
観測量

観測量に合うよう設計
試行錯誤が必要

High-throughput

経験的ポテンシャルモデル

パラメータが決まれば計算可能
計算コストが低いので種々の検討が可能

個別の実験事実を考慮しなくても計算可能

第一原理計算第一原理計算を再現するよう設計
Vmodel ~ VFP というフィッティングをすれば良い

計算コストが高いので検討範囲が必然的に狭まる!5

PHYSICAL REVIEW LETTERSVOLUME 61, NUMBER 7

15 AUGUST 1988

a =0. 118, 0.424) ' for silicon and (9s5p)/[3s2p] with pfunctions (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
oxygen. We assume no electron orbitals around theSi-0 distance is changed, and the D2d mode, in which
point 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
PHYSICAL REVIEW
15 AUGUST sur1988
It hasLETTERS
turned out that these three potentialchanging the Si-0 distances or 0-Si-0 angles with threeenergydifl'erent modes. One potential energy surface is depictfaces can be fitted well by a sum of pairwise interatomic
we stretchtheallfunction
the Si-0
bonds keepwhich employed
ed in Fig. 1,We
in4-have
a =0. 118, 0.424) ' for silicon and (9s5p)/[3s2p] with p t potentials.
form,
Hartree-Fock法でSiO
4 + 4 point chargeモデルを元に、ポテンシャルを作成'=0.
the Td symmetry of the cluster. We have also obfunctions =U;a b(r)+f
for negative-ion states059) o(b;+bofing—c;cj/r',/c'""~(1)
U/(r) (with)]
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 theSi-0 distancewhich
and the repulsion,
point charges.
which•changed,
Si-O間の電荷移動量をモデル化
D2d mode,andin disperconsists of Coulomb interaction with some corrections discussed
below,is Born-Mayer-type
Phys.
Rev.
Lett.
61,
86between(1988).
0-Si-0
Si-0 distance.
is thenare varied
constant
potential
obtainedanglesenergysiveTheinteraction.
r issurface
by ourHerethedistance
and a; (b;)
atoms
is theeffective
radius with
(softnessof the ith
parameter)• interactions
静電相互作用をlong-range,
Si-0
distances
It has Coulomb
that
these threewith threeturned out
or 0-Si-0
potentialchanging
energy sur- short-rangeに分離
with the339,(1989).
atom
theNature
standard
force209
1 kcal angles
A
mol. We also include
with thepoint charges.
fodifl'erent
modes.
well by because
potential
surfaceof isthedepictcan beinfitted
a sum ofthepairwise
A cautionmust Onebe madein theenergy
evaluation
Coulomb faces
interaction
the cluster,effectiveinteratomic
charge in• クラスターモデルから固体系へ移植
haveform,from a Si
employed
the bulk, Q;, is different from that in a Si04
terms of theWefractional
cluster, Q;: t In potentials.
transferredAn, function
charge, the—2/)ne and— bond,atomU/(r)to an=U;/c'""~atom b(r)+f
we a+ahave Qp=
have Qo= —(I+An)e andper Si o(b;+b
Qs;=4/)ne, • while
SiO2weの様々なpolymorphの安定構造を再現(1))] —c;cj/r',
J)exp[(b;+bi
J r)/(—=4/)neWe.
theinteraction
Coulomb
inexpresstheclustersumof
as along-range and short-range parts as
Qs;
which consists of Coulomb interaction with some corrections discussed below, •Born-Mayer-type
repulsion, and disper6配位, ４配位の差に軌道を仮定しなくてもよかった'=Q;Q,Here[1 —sive Uij~'""interaction.
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)—goo(r) =[1+11((r)/8+3((r)— because the effective charge in/4+((r)in /6]exp(
Agso(r)caution=(I+fr)exp(
must be made 2/r),in the evaluationof the Coulomb interaction
the cluster,2(r).
the bulk, Q;, is different from that in a Si04
cluster, Q;: In terms of the fractional charge, An, transferred from a Si
The correction,in the—long-range part (the first—
g&(r),peratom to an
atom
Si
2/)ne
and Qs;=4/)ne, while we have Qo= —(I+An)e andbond, we have Qp=
Iterm)=4/)ne
involving Q reflects the distribution of the excess
the equilibrium bond length, is shown to behn7)-forWe. express the Coulomb interaction in the cluster as a 0sum
of. long-range
and short-range parts as
Qs;charge of oxygen, for which we assume a hydrogenlike
a function of the Si-0 distance (Fig. 1). The feature—g~)(r)]/r+Q;Q)gii(r)/r,
Uij~'""
with '=Q;Q,
orbitala radius [1I/(
here. The radius is of the orderthat the atomic charge varies with the bond length clearof the ionic radius of 0, so that we employ I/(=1. 4 A(2)ly indicates a many-body character of interatomic forces'
Since 2/r),
followingthe —short-range
remaining=(I+fr)exp(
The
in the /4+((r)
covalent /6]exp(
system. —gso(r)Pauling.
2(r).physical mechanism will begoo(r) =[1+11((r)/8+3((r)
part is expected to be insensitive to the environment, we
discussed in more detail elsewhere. To concentrate on
The correction, g&(r), in the long-range part (the firstuse the bulk Q; there.
Once the parameters a;, b;, and c;
approach, however, we have used conI the pair-potentialterm) involving Q reflects the distribution of the excessforbond length, is shown to be0. the equilibriumhn7)
are optimized fromthe cluster calculation, we switch Q
stant Q; (Q;). Since small-cluster results are insufficientfor which we assume a hydrogenlike
Si-0
charge of oxygen, UPi""distancefunction
theof(Fig. 1). The featureaback into Q (i.e. ,=Q;Qi/r) in the bulk simulato determine the long-range Coulomb interaction, andorbital with a radius I/( here. The radius is of the orderwith the bond length clearvaries
theatomicthatchargetion. Thus the final pair potentials have the same funcbecause the absolute value of the Mulliken charge itselfof the ionic radius of 0, so that we employ ' I/(=1. 4 A
ly indicates a many-body character of interatomic forces
tional form as suggested
by Gilbert and Ida.
depends on the choice of basis function, we have not in'Since the remaining short-range
following Pauling.
in the covalent system. The physical mechanism will be
From the cluster calculation, the charge obtained byInstead we triedcluded An in the fitting procedure.
wepart is expected to be insensitive to the environment,
discussed in more detail elsewhere. To concentrate on
the Mulliken
analysis, which is Qp- —1.7e (i.e.,several fixed values of An around the Mulliken charge.
use the bulk Q; there. Once the parameters a;, b;, and c;
the pair-potential approach, however, we have used conparameAmong the trial values studied here, the fitted
are optimized from the cluster calculation, we switch Q
stant Q; (Q;). Since small-cluster results are insufficientparameters.
best
thecrystal
with
An
0.
6
reproducetersback -589.
simulainto 4Q (i.e. , UPi"" =Q;Qi/r) in the bulk -1.
to determine the long-range Coulomb interaction, and0and c; is nonlinear,
Since the fitting procedure for a;, b;,tion. Thus the final pair potentials have the same funcbecause the absolute value of the Mulliken charge itself
Si 'obtained.
We have
are
set
of
more
than
one
tionalform as suggested by Gilbert and Ida.
of basis function, we have not independs on the choiceparametersO
structure
~From
the
which
optimizes(Tablechosentheone
I)
1.2 Oby5- cluster calculation, the charge-589.thep obtained
Instead we triedcluded An in the fitting procedure.
in a static simulationand compressibility of a-quartz'—e'theC Mulliken
1.7e (i.e.,analysis, which is ~Qpseveral fixed values of' An around the Mulliken charge.
Z0
the program wMtN.
by
2
Among the trial values studied here, the fitted parame-1.4—589.6have parameters.
We
performed
Now
we turn !6to the MD study.
ters with An 0.6 reproduce the best crystal-1.0
of silica
test forfor four polymorphs
thedynamicalg) -589.4
procedure
Since
the fittingstabilitya;, b;, and c; is nonlinear,

フィッティングによるモデル作成

Tsuneyuki Potential

0

0

0

0

I

高圧下でのSiO 新構造を発見！

surface dimer buckles slightly away from the molecule. The
Previous calculations by Wu etconstant orientation of H2 along the unstable mode is consisJing and Whitten have found mu
tent with a cool rotational energy distribution and large
with the prepairing mechanism. Ouchanges in H–H bond length suggest that the product will be
due to different approximations us
vibrationally excited.
tion and in the electronic structure
The B-LYP transition state is 23 kcal/mol higher in ention by Nachtigall et al.10 is mostergy than the desorption products. Since most H2 molecules
Constraints that they imposed on thdesorb without significant vibrational, rotational, or translaatoms(10),prevent
anticorrelated
J. Chem. Phys. 103
8 the(1995).
tional excitation, the excess energy at the transition state
above. We have used their basis se（おそらく）初めてのNNによるフィッティング論文
must be deposited in surface modes. Somecan models
be of potential
et al.: Neural network
surfaces
Blankinsighttional~Becke3LYP! with4133the GUA
obtained from classicalrequired
estimatesof hidden
the energy
vibrational
only one morenode than thein
networkusedcalculate the energy at geometries
in the two-dimensional case, though the prediction error ismodes at the transition
state
and
the
assumption
that
theClassical modelの入力を低次元化する目的
with no constraints. The resulting
not as low. Nevertheless, it is remarkable that such a small
to the net topology
can achieve noreasonableaccuracyafter
surface and hydrogenchange
molecule
exchangeenergy74 kcal/mol ~without ZPE! at theon a substantially more complicated problem.
The predictionsof potentialfrom a cubic predict
spline fit to that
• CO/Ni(111)のCO吸着に関する研究
reaching the transition state.
B-LYP
calculations
etries, respectively, compared todata on the 73737 grid is an order of magnitude more aca hydrogen molecule stretched
as
in
our
transition
state
will
curate than the 3-7-1sl network. Increasing the density ofNachtigall et al. The difference re• COの３次元配置で与えられるポテンシャルを３つの特徴量でモデル化
data in the z direction allows the inclusion of morebe 15 kcal/mol abovetraining
the nodes
potential
minimum.
Thehowever.
first exhiddenin the network
without overtraining,straints. Increasing the basis to 6-3
training data on a 737313 grid and fifteen hidden• 7 x 7 x 11 samples, 3-15-1という構成
cited vibrational stateUsingis18
kcal/molabove
theminimum,
nodes, the neural net is comparable to the spline in accuracy,
the Becke3LYP prediction to 72.3
a MADdesorb
below 0.02 kcal/moleanda maximum
error
oftheronly
in
this
state
if
it
is
suggesting that H2 canwith
ometry. We have also done B-LYP
only 0.14 kcal/mole. This absolute error is as small as that
• Extend Kalman filterを使って最適化
achievedin the two-dimensionalcase, butfraction
as a fraction ofmally excited. This may
rationalize
the smallofthemolbasis at the same geometry, and fi
energy range, the error in fitting the three-dimensional model• 低次元NNで、path-integralを計算し、量子効果を見積もっている
ecules that desorb withvibrational
excitation.
the vast
is smaller
by half ~0.15%!.
Figure 2 shows For
the prediction
mol, only 0.5 kcal/mol below the v
for this network. The largest errors occur at energies
majority of molecules errors
whichin thevibrational
groundclose
to the desorb
minimum, where
the curvature
with respecttoson of the B-LYP and Becke3LY

ニューラルネットによるフィッティング

the z coordinate is greatest.
DFTの結果をモデル化することも試みられている
These training sets include an order of magnitude more

4136

Blank et al.: Neural network models of plarge basis suggests that differenc
calculatexamples than in the two-degree-of-freedom case. However,negligibtionals account for about 5 kcal/mothe x, z, u network can be trained in only 80 epochs, which
programis 50–70% of the number of epochs needed to train the x, u
dessimpledicted E a . This raises the questi
network. This demonstrates one of the strengths of the neuralnetwork
network: Each of the network weights depends on many data
more accurate. Nachtigall et al.
hatypically• Si(100)上のH2発生のクラスター
points, and training can be effected either by presenting moredata comdata, or by multiple presentations of the same data. The trainparisons of several DFT metho
training
ing data contain redundant information, and the neural net is
モデル計算
37
network
an efficient method for reducing a large number of datachemistry.
Unfortunately,
there
Our
points to a smaller set of weights. In this sense, the neural net
structure• 750 sample, 12-7-1という構成
has more in common with fits to assumed functional forms
analog of the electronic structuretechnolo~which use a small number of parameters to determine athousanwhich could provide an experimen
global potential! than with a spline ~which forces a local fit,
• 誤差2.1 kcal/mol.とやや大きい
ometryusing as many parameters as there are data points!. This
Since the initial submission
calculat
redundancy also can help to make the network less sensitive38
the info
to noise in the input training data than a spline.
learned of DFT slab calculations
network
The time required to train one of the larger networks is
S. Pai and D. Doren, J. Chem. Phys. 103, 1232 (1995).
glectedless than 90 min, or an average of one minute per epoch, in
E desZ
FIG. 1. Configurations along
the unstable normal mode at the transition
a of 55– 60 kcal/mol ~without
train or
MATLAB on a Sun Sparcstation 2/Weitek mP. In applicastate. Selected bond lengthstions
forwhere
the the
LSD-optimized
transition@shownMost ofpotential will be evaluatedmany statetimes, thetion of whether cluster methods
states, a
training timeis negligible
compared Å;
to thetime –H
spent !51.67
evaluÅ; r~Si
–H
!51.70r~SiÅ;
in 1~b!# include: r~Si1 –Si2!52.37
des
ただし、そもそものシミュレーション手法がそこまで進んでいない。。
11
12
be theE aerror.distribution
The forslab
calculations
Bm
ating the function. Under such circumstances, much longer
FIG. 2. Prediction
a three-dimensionalmodel of CO/ use
Å. beCorresponding
values for the
r~Si2 –H2!52.18 Å; r~H1 –H
priate f
Ni~111!. The potential minimum is at 231 kcal/mole. The network has a
2!51.00
training
times could
tolerated.
trained from a 737313
training set. Results
are showngeneralized
gradient
approximatio
where reB-LYP-optimized geometry are:Å; structure,!7 r~Si1 –Si2!52.47 Å; r~Si1 –H1!51.69 3-15-1sl
for a test set of 2000 points. The MAD over the test set is 0.022 kcal/mole.
43
IV. PATH INTEGRAL TRANSITION STATE THEORY~a!
Scatterplotoftheprediction
error
againstthetrueenergy.
~b!
Histogramused the GAUSSIAN94 programstates.tor~Si –H !51.68 Å; r~Si –H !52.18 Å; r~H –H !50.99 Å.

NNポテンシャル発展の鍵
表面化学反応のシミュレーション

方法論に求められること
1.電子状態変化が取り扱えること
2.エネルギー超曲面上でのダイナミクス
3.ダイナミクスの統計性• 電子状態変化を伴うため、DFT計算が必須• DFTのみでは計算が大変で統計が集まらない

複合的な方法論を構築しないとどうにもならない

Pd表面でのH2分解反応!8

PHYSICAL REVIEW B

VOLUME 57, NUMBER 4

15 JANUARY 1998-II

Ab initio quantum and molecular dynamics of the dissociative adsorption of hydrogen on Pd„100…
Axel Gross and Matthias Scheffler

Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, D-14195 Berlin-Dahlem, Germany
AB INITIO QUANTUM AND MOLECULAR DY
57~Received 2 May 1997!
The dissociative adsorption of hydrogen on Pd~100! has been studied by ab initio quantum dynamics and abinitio molecular-dynamics calculations. Treating all hydrogen degrees of freedom as dynamical coordinates
implies a high dimensionality and requires statistical averages over thousands of trajectories. An efficient andaccurate treatment of such extensive statistics is achieved in a three-step approach: In a first step we evaluate
要請（１）the ab initio potential-energy surface ~PES! for a number of appropriate points in configuration space. Then ~asstep 2! we determine an analytical representation that serves as an interpolation between the actually calculated
Pd上のH
2を様々な配置でそのエネルギーを
points.
In an independent third step dynamical calculations are performed on the analytical representation of
the PES. Thus the dissociation dynamics is investigated without any crucial assumption except for the BornDFTで計算（250点）
Oppenheimer
approximation which is anyhow employed when density-functional-theory calculations are performed. The ab initio molecular dynamics is compared to detailed quantum-dynamical calculations on exactly
the same ab initio PES. The occurence of quantum oscillations in the sticking probability as a function of
kinetic energy is addressed. They turn out to be very sensitive to the symmetry of the initial conditions. At low
表面構造の「対称性」を持つ基底関数で上のkinetic energies sticking is dominated by the steering effect, which is illustrated using classical trajectories. The
結果をフィッティング
steering
effect depends on the kinetic energy, but not on the mass of the molecules, as long as no energy
transfer to the substrate atoms is considered. The comparison between quantum and classical calculations ofe.g.
Symmetry-adaptedfunction
the sticking
probability shows the importanceof zero-point effects in the hydrogen dynamics. The dependenceofsticking probability
the angle
of incidence is analyzed; it is found to be in good agreement with
[J.theChem.
Phys., on104,
8(1996)]experimental data. The results show that the determination of the potential-energy surface combined withhigh-dimensional dynamical calculations, in which all relevant degrees of freedom are taken into account,
要請（２）leads to a detailed understanding of the dissociation dynamics of hydrogen at a transition metal surface.
@S0163-1829~98!05204-7#
フィットしたモデルポテンシャル上で、H2のGross-Sheﬄerの方法

量子ダイナミクスを計算

FIG. 1. Contour plot of the PES along a two-dimensional cut
through the six-dimensional coordinate space of H 2 /Pd ~100!. The
inset shows the orientation of the molecular axis and the lateral H 2center-of-mass coordinates, i.e., the coordinates X, Y , u , and f . The
coordinates in the figure are the H 2 center-of-mass distance fromney of the particles,
which
on the distanceorder of
100–1000
the surface
Z and requires
the H-H interatomicr. Thedashed line is
optimum reaction
are in eV
per Htheself-consistentthe
calculations
for path.
eachEnergies
trajectory.
Since
de2 molecule.
The contour spacing is 0.1 eV.

要請（３） I. INTRODUCTION
The 繰り返して統計性を担保・実験と比較
dissociative adsorption of molecules on surfaces istermination of total energies and forces is still a heavy com-

one of the fundamental reaction steps occuring in catalysis,
putational task,degrees
the numberof trajectories obtainable in suchof freedom. They can be represented, e.g., by the
corrosion, and the hydrogen gas storage in metals. In prin- !9 a ‘‘during the center-of-massjourney’’ abcoordinates
initio molecular-dynamics
X,Y ,Z, the interatomicsimudistancer, andpolar andwellazimuthal
angle
of1,2
theThesemolecular
axis ulation is limited
to thenumbers
below
100.
numciple, and this is the approach emphasized and taken in this

1
rot

V 5

(

m50

1

(

n

and

G52 p /a isurface reci
between Pd
note that W
equivalent a
H 2 /Cu~100!
coordinates
nates.
It turns odifferent coeters necess
potential pobtained41 u
generalizedexchange alinear augm
54 and refev (s) are decalculations
The intro
term V vib @E
minimum ethrough th
H 2 /Pd~100!
however, itmetry. It is important to note that the terms x1, . . . ,x8
represent only the minimum symmetry employed in
the data set, any other symmetry can then be learned
by the neural network.
Without such an approach, the number of training
points necessary for an accurate fit becomes much larger. A similar eﬀect can also be achieved by including
the gradients of the potential in the optimization procedure. A 8–24–18–1sl neural network, i.e., with two hidden layers, each with 24 or 18 nodes, respectively, and685 weights, a sigmoidal transfer function in the hidden
layers and a linear function in the output layer was used
Pd上のHfor2を様々な配置でそのエネルギーを
the fit.
In an activated system such as H2/K(2 · 2)/Pd(1 0 0)
DFTで計算the configurations
close to the valley of each elbow plot
are most crucial for the reaction dynamics, as revealed
by molecular dynamics simulations. During the fit, forthese configurations the weighting factors ui in the costfunction Eq. (2) were chosen to be 10 times larger than
表面構造の「対称性」を持つ関数値を
for other regions where the potential energy is significantly elevated. Thus, the region accessible for molecu「INPUT」としてNNでフィッティング
lar dynamics simulations is represented with a higher
accuracy.
The training and test root mean squared errors forthe neural network fit as a function of the number of
要請（２）so-called epochs are shown in Fig. 2. An epoch corresponds to the process of presenting the network one
フィットしたモデルポテンシャル上で、H2のwhole cycle of the data. The training and test errors
dropvery quickly and start to saturate at around 100
量子ダイナミクスを計算
epochs. Both error functions continue to go down withlonger training time, but after 100 epochs the test error

below the desired accuracy of 0.1 eV.
Fig. 3 illustrates the accuracy of the neural network
Chem.
Phys. Lett.
fit. Onetwo-dimensionalcut 395,
through210
the (2004).
six-dimensional potential energy surface, a so-called elbow plot,
obtained by the DFT calculations (Fig. 3a) is comparedto the corresponding neural network elbow plot (Fig.
3b). This cut corresponds to the dissociation above thehollow site between two K atoms which is hindered by

より柔軟なフィッティング

要請（３）繰り返して統計性を担保・実験と比較

H2 distance from the surface (A)

5.0
K

Pd

4.0

Pd

0.1

3.5
0.1

3.0
2.5
2.0

0.2

0.3
0.3
0.2

1.0

0.0

0.0
0.5

0.7

0.5

1.5

0.5

(a)

K

4.5

-0.2

-0.6

0.2

0.0

-0.5
-0.7

-0.7

-0.4

-0.3
0.7

0.5

1.0

1.5

2.0

H–H distance (A)

2.5

0.5

1.0

1.5

2.0

(b)

H–H distance (A)2.5

Fig. 3. Two-dimensional cuts through the six-dimensional PES of H2/K(2 · 2)/Pd(1 0 0). The molecular configuration within the surface unit
cell is depicted in the insets: (a) two-dimensional spline interpolation of
ab initio results; (b) two-dimensional cut through a six-dimensional
neural network interpolation. The dots in (b) mark the calculated gridof ab initio energies for this specific cut.

表面が持つ対称性の情報を含めた状態での柔軟なフィッティングが可能に!10

NNフィッティングの諸問題
Symmetry-Adapted Functionの設計
fcc(111)の表面形状に合わせてフーリエ変換を元に構成する方法を提案[1][1] J. Behler, S. Lorenz, and K. Reuter, J. Chem. Phys. 127, 014705 ︎(2007︎)

ただし問題に合わせて複雑な関数を設計することは決して容易ではなくもっと簡便な表現が必要に。

３次元配置を入力にしたNNの問題• 入力層が原子数に依存するため、異なる系には適用できない（汎用性）• 同種粒子の入れ替えに対する対称性が破れてしまう• 同様にポテンシャルが満たすべき対称性が保証されない

入力表現・汎用性・対称性がキーワード!11

Behler-Parrinelloのアイデア
1. 入力として新しいSymmetry Function（対称性関数）を導入
G1i =

all
X

e ⌘(Rij

R s )2fc (Rij ),

(2)

j6=i

G2i = 21 ⇣

all
X

⇣

2

2

2

(1 + cos ✓ijk ) e ⌘(Rij +Rik +Rjk ) fc (Rij )fc (Rik )fc (Rjk ),(3)

j,k6=i• 各原子周辺環境を２体・3体の指標（Symmetry function）で数値化• これらの指標は系の回転・並進で普遍

week ending2.
全エネルギーを原子ごとの寄与に分割、原子種ごとに同じNNを構成する
P
H
Y
S
I
C
A
L
R
E
V
I
E
W
L
E
T
T
E
R
S
6 APRIL 2007
PRL 98, 146401 (2007)

input and the output layer are one or more ‘‘hidden layers,’’同種原子の入れ替え対称性を担保
each with •a certain number of nodes. All nodes in eachlayer are connectedto the nodes in the adjacent layers by• 原子数が増えてもシステムサイズで入
real-valued weight parameters, which initially are chosen
力層が変らず学習済みNNのサブネッ
randomly. For a given set of coordinates the output of the
NN is then givenby the expression
トを増やすだけで良いため
!
"
#$
3
2
X
X
!
2
2
2
1
1
1
E i ! fa w01 "wj1 fa w0j "w!j Gi : (1)
j!1

!!1

FIG. 2. Structure of the neural network as applied in this Letter!12
to a system containing three atoms. The Cartesian coordinates of汎用NNの構成の流れ
2

1

Symmetry function vector Gk={Gk1,Gk2,…}で原子の周辺環境を記述

4
3

G1

NN1

E1

原子配置

G2

NN1

E2

G3

NN1

E3

G4

NN2

E4

DFT計算

E

全エネルギーE• 正味のINPUTは、サブネットの入力 x 原子数• 学習すべきNNのパラメータは、原子種ごとに共通。
• また入力ベクトルGの次元は固定できるので原子数を増やすことは容易!13sum to range over a number of “representative points” drawn productnormalized to have unit lecarefully
controlledduringtraining.
from the full training database (Nt ≪ Nfull ). The numberof" (MB)# !
PHYSICALREVIEWB95,094203(2017).
The
kernelfunction
fo
(MB)
(a
Both
for
2b
and
3b
contributions,weuseasquared
k
q
,q
=
p
t
a
nn
representative points differs for each descriptor and must be
product
exponential kernel [35]′nn lcarefully controlled during training.
⎡⎤" (MB) (MB) # !
"#2(d ) we (duse)
Both for" 2b and3b contributions,a squared and we findkitqadvantageous
to ra
,q
=
!
q
−q
#
t
a
1
ξ,t
ξ,i
) (d )⎦ , (3)for a sharper distinction between
K (d ) kernelq(d= exp ⎣ −exponential
i ,qt [35]n22VOLKERL.
DERINGER
AND
G
ÁBOR
CS
ÁNYI
θ
ξ
gives the final kernel
⎡
⎤
ξ
"
#
2
(d )
(d )
GAP: Gaussian
Approximation
Potential
and we find it" advantageo!
#
mostatq−q"#
1
ξ,t ofξ,i(MB) (MB) (MB)(d ) index(d )whereKξ(d )isqan
runningoverthecomponents
the⎣⎦
KDFT distinction
qa ,qt function
30
for a sharperb, (3)i ,qt (d ) = exp −2
scheme2b
only
descriptor
vector q . In the case2of pairs, theθdescriptor
has✓ ガウス過程回帰で相互作用を記述
ξ
gives
the
final
ξ
2b+3b
only kernel and a cu
This
dot
productkernel wereis ausone single
scalarcomponent(namely,thedistancer
between20
12
SOAP only✓ Cの相互作用を2-, 3-, many-bodyと展開"
2b+3b+SOAP
Exchpower spectrum
descriptor,
asi(MB) (MB)the twoξ atoms
involved):
where
is
an
index
runningoverthecomponents
of
the
K
q
,q
approxi
kernelを利用✓ 2- and 3-body: Gaussiana inte10(up to normalization) to the)training
descriptor vector qq(d(2b)
. In
the
case
of
pairs,
thedescriptor
has
=
|r
−r
|
≡r
;
(4)
2
1 kernel
12
has long✓ Many-body: Simple dot productneighbor
densitiesThis
dot product kerne
one single scalar component (namely, the distance r120 between
the de fa/
0Further0 descripto
fortwotriplets,
weinvolved):
do not directly use the natural coordinates r12 ,power spectrum
Descriptor
the
atoms
/ higherr13 , and r23 , but a different form to enforce symmetry
over
ρaex
d R̂ // more
–10
(up
to
normalization)
to
t(
(2b)
dispersi
permutation= |ratoms≡r312[52]:;(4)
2 −r12| and2-body: of the qneighborneighbor
densities
0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
to the G⎛⎞
expression
for the
C–CThedistance(Å)+ rthe0 total/13 natural coordinates r12 ,for triplets, we do not directlyr12use/
therefore
given
bydimer. This
(3b)
2⎠
⎝
FIG. 2. Potential-energyscans for an isolated
carbon
(r
−r
)
q
=
.
(5)
12
13 enforce symmetry overr13 3-body:, and r23 , but a different form
to
d Our
R̂ //
plot, with DFT data as reference (blue), allows us to assess
the
use
!
!
r23 2 and 3 [52]:of different structural descriptors:Eall =three combined(2b) 2 are needed for (2b)(2
initiopermutation of the neighbor atoms
(δ
)
αt ab
K
a high-quality fit (see text).
driven⎛⎞thentiNote that with this choice of descriptors, the first term in
The expression foraretheMany-body:
r
+
r
yield
SOAP
12
13
! (3b)w
To illustrate the role of the combined descriptors,!
we useEq. (1) is equivalent(3b)to a pairpotential,
andthesecondisa(3b)us
therefore
given
by αdone2⎠
different (and increasingly
complex)GAP⎝(r+ (δmodels)to2 compute
K
−r
)
q
=
.
(5)
t
(Smoothoverlapofatomicposition)
12
13(maximthe potential-energy
generic three-body potential, but in the GAP frameworkboth curve for an isolated carbon dimer; thesewave ext !
j2 !
models have been fitted to the full bulk and surface
trainingr(2b)23
do not impose constraints on the specific functional form.
finite-b
E = (δ )α∆Energy (eV)

そのほかの表現とモデル

set described below that additionally incorporates DFT
data!
!
t
was(MBapp(MB)
2
between 0.8 and 3.7 Å in small+increments.
The
results
For the many-body term, we use the recently points
introduced
(δ
)
α
t t
i
iteration
are summarized in Fig. 2: GAP models using 2b descriptors
smooth overlap of atomic positions (SOAP) [53] only,
descriptor,
Initi
a
t
or a combination of 2b+3b, reproduce the minimum
(3b)
2
ab
14 GAP and
the repulsion
which has proven successful in generating
modelsfor at small C-C distances reasonably well, but wasinitio
fitt

対称性を保った局所表現による予測が主流に。
Note that with thischoice of descriptors, the first term in!!
Eq. (1) is equivalent to a pair potential, and the second is a
+ (δ )

現状のまとめ
機械学習ポテンシャル設計の流れ
観測量

経験的ポテンシャルモデル
High-throughput

個別の実験事実を考慮しなくても計算可能

機械学習ポテンシャル

だが数は少ない

NN, GAP, LMM• 全エネルギー/Forceをフィット• ３次元座標ではなく、symmetry

第一原理計算
計算コストが高いので検討範囲が必然的に狭まる

invarianceな環境特徴ベクトルを入力• モデルの関数系を仮定しないで良い
注意点• 機械学習ポテンシャルは、もはやブラックボックス• ポテンシャルの性能を保証できる領域が曖昧（学習データに依存・補間）• 結局、学習データの取得コストが高いので学習コストは下げたい!15

ここから私たちがやるべき仕事
Throughputを上げて初めてわかる課題設定• 「精度よくフィットができた」は結局意味がない• モデル構築自体の価値は下がるので、アプリケーションが重要• 小数個の計算では不十分で、高速にたくさん計算できてわかること• 分布に対する統計性の問題や、網羅的な探索が必要な問題
Deepにしすぎないモデル構築• 学習データ取得コストを劇的には下げることができない• 少ないデータ数でもある程度性能が出る機械学習技術が必要• NN自体をシンプルにする• 入力ベクトルを何かしらの方法で減らす、etc.
!16

W. Li, Y. Ando, ando S. Watanabe, J. Phys. Soc. Jpn. 86, 104004 (2017).

アモルファス内の原子拡散経路の網羅探索
大問題：高精度に計算できるDFTではコストがかかりすぎる（N3 回の構造最適化 x 計算時間 T) N~50, T~1 hour以上 = 5.7 year
さらに拡散障壁の計算も必要…

Etot = Eamorph +

ECu +

Eopt

Cu挿入による

構造緩和によるエネルギー変化

エネルギー変化Neural Networkで予測✓ アモルファス母体を「場」として取り扱う✓ 学習対象は拡散粒子（Cu）のみ✓ ３元系でありながら小規模なネットワークで記述可能✓ 事実上１粒子の計算なので、エネルギー評価も容易✓ アモルファスは動かせないのでMDはできない!17

拡散ネットワークの全体図

T. Yayama, W. Li, Y. Ando in prep.

非一様拡散経路のモデル化
原子の準安定点が作る拡散ネットワークを可視化
解析の概要

53
35

41

42

40

39
47

36
38

43

37

45

✓ ホッピング確率を使ってマル

57

46

44

51

29

52

コフ過程の定常状態を計算

50

32

48

✓ 重なりの弱い複数の定常状態

49

1
33

14
6

12

13

11

9

56

55

54

24

10
8

21
18
17

16

23

15

るノードの特徴を同定✓ アモルファス内部のイオン拡

222

散を阻害する準安定点の原子

0

26

25

31

30
19

を発見
✓ 定常状態間の接続を妨げてい

28

7

3

34

構造的な特徴を発見

20
5

4

27

第一原理計算結果から複雑ネットワークモデルまで抽象化することが可能に!18

we examined several NN potentials constructed with different
choices of angular SFs. In all the cases, the NN potentials weretrained with the Broyden–Fletcher–Goldfarb–Shanno (BFGS)
第1表
algorithm,42 using the training set as reference data, while
the testing set was used to validate them. Table I shows thehidden layers with 15 hidden nodes in each layer. The hyperbolic tangent activation functions were employed at the nodesof the two hidden layers, and a linear activation function was
employed at the output node. This NN architecture contains
3123 parameters in total.

コンパクトなNNの設計

TABLE I. The root mean squared errors (RMSEs) of the NN energies obtained for the training set and the testing set of the Li3 PO4 system using different input
symmetry function sets. The number of parameters in the NN potential is also given for each architecture.

Species

Types of angular symmetry function set

1

Li
P
O

Li–LiLiP–LiLiO–LiLi

Li–PP
P–PP
O–PP

Li–OO
P–OO
O–OO

2

Li
P
O

Li–LiLiP–LiLiO–LiLi

Li–PP
P–PP
O–PP

Li–OO
P–OO
O–OO

3

Li
P
O

Li–OO
P–OO
O–LiLi

4

Li
P
O

5

Li
P
O

Li–LiP
P–LiP
O–LiP

Li–LiO
P–LiO
O–LiO

O–PP

P–OO
O–PP

Training

Testing156-15-15-1
156-15-15-1
156-15-15-1

7833

4.9

5.7

90-15-15-1
90-15-15-1
90-15-15-1

4863

5.0

5.5

46-15-15-1
46-15-15-1
68-15-15-1

3213

5.0

5.6

24-15-15-1
46-15-15-1
46-15-15-1

2553

5.7

6.0

24-15-15-1
24-15-15-1
24-15-15-1

1893

8.3

8.6

Network
Li–PO
P–PO
O–PO

RMSE (meV/atom)

No. ofparameters

パラメータ数を半分以下にしても、性能はほとんど変わらない
入力ベクトルの特徴選択は学習効率を上げるポイント!19

It should be noted that the large formation energies(⇠5 eV) obtained here mean very small probability of formation 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.

Li3 PO4 supercell, and µLi is the chemical potential of Li
obtained from DFT calculation of a perfect Li bcc crystal. The
supercells with a vacancy have been generated by removing
a Li atom and then relaxed the structure. Since the supercell
contains 12 Li atoms that are not equivalent to one another,
we can obtain 12 different Li vacancy formation energies.
Figure 4 shows the comparison of the vacancy formation

NNの移植性の確認

C. Diffusion paths and network

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 neighbor vacancy site when the distanceNN

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

DFT for small cell

Simulation for large cell

DFTと同様の同型分布・拡散係数の実験値を再現
!20

まとめ
• 原子間相互作用は物質シミュレーションの最重要テーマ• 数式として扱いやすい・計算機に乗りやすい・共有結合性の表現など
時代の変遷に応じて相互作用の形式に変化• 実験値に合わせ込む最適化から第一原理計算結果のフィッティングへ• 対称性適応関数で表面系のポテンシャルをフィットする試みからNNのインプットである「対称性関数」のアイデアへ繋がる• ３元系のアモルファス問題に対して、２成分を「外場」として扱うことで低次元のNNを作成、アモルファス内部の原子拡散経路を全て探索することに成功• 物質の特徴から、重要な対称性関数を選択することでNNを簡略化す
ることに成功!21