精度保証された機械学習分子動力学法：
自己学習ハイブリッドモンテカルロ法
原子力機構システム計算科学センター
理研AIP

YN, M. Okumura, K. Kobayashi, and M. Shiga,“Self-learning Hybrid Monte Carlo: A First-principles Approach”,
arxiv:1909.02255
1

永井佑紀

2

このセミナーについて
国内ワークショップ2018年6月 場所：大阪大学
Y. Nagai : Self-learning Monte Carlo Method

国際ワークショップ2019年10月31日~ 場所：京都大学“物理学者,機械学習を使う ー機械学習・深層学習の物理学への応用ー ”
朝倉書店(2019/10)
第5章”自己学習モンテカルロ法”

ディープラーニングと物理学2020

オンライン

場所：インターネット

最新の進展の紹介

3

本日のお話
原子や分子の運動をシミュレーションしたい！
日常に現れるあらゆる物質（固体、液体、気体、生体）をシミュレーション

分子動力学

原子や分子に働く力とは？

原子は電子と原子核からなる->量子力学と電磁気学を使えば良い
電子や原子の数はアボガドロ数：無理！
3N次元シュレーディンガー方程式(N~1023)

第一原理計算（密度汎関数理論）
3N次元波動関数と3次元電子密度の一対一対応：固体の計算が可能に“第一原理計算”：原子位置を与えればエネルギーが計算できる
実用上は近似が入る…化学の人は”第一原理計算”と呼ばないetc

力はエネルギーの微分！
第一原理分子動力学計算：高精度なシミュレーション

めちゃめちゃ計算が重いのでなんとかしたい！

機械学習使ってみよう！

本研

究

学習結果に依存しない精度が欲しい！

4

アウトライン
機械学習分子動力学法
自己学習ハイブリッドモンテカルロ法
適用例
まとめ
シミュレーション中にニューラルネットワークを改善し、かつ、第一原理分子動力学計算精度で計算を行う
YN, M. Okumura, K. Kobayashi, and M. Shiga,“Self-learning Hybrid Monte Carlo: A First-principles Approach”,
arxiv:1909.02255

5

機械学習分子動力学法

6

クイズ1：この動物は何？

猫ですBasile Morin / CC BY-SA (https://creativecommons.org/licenses/by-sa/4.0)

色々な猫を見たことがあるので、猫とわかる7

クイズ2：この動物は何？

サイガです
ウシ科サイガ属に分類される偶蹄類（モンゴル）：絶滅危惧種
https://ja-jp.facebook.com/9GAGCute/photos/saiga-antelope-a-priority-species-for-conservationthe-break-up-of-the-former-uss/816439751883761/

見たことがなければわからない機械学習とは
猫画像その他関数f(x)y猫?

猫です

yes or no

学習プロセスx

大量の入力データxを用いて、f(x)=yを満たす関数f(x)を決める

例：囲碁棋譜データ→勝利の方程式
囲碁プロ棋士に勝利

機械学習分子動力学計算

原子配置

ある原子配置におけるエネルギーを使って高速にシミュレーション

機械学習とは教師あり学習
大量の入力データxを用いて、f(x)=yを満たす関数f(x)を決める一番シンプルな例

y = ax+ b直線で近似（線形回帰）

インプットが複数→xがベクトル

y = Wx+ b
もっと表現力を高めたいW

yy = W2 f(W1x+ b1)+bf:非線形関数

y = W3f(W2 f(W1x+ b1)+b2)+b：x

深層学習
…

猫です
ニューラルネットワーク

10

機械学習の物理学への応用
例：相転移の検出
フェルミオン負符号問題のある量子系での量子モンテカルロ法による相転移検出

イジング模型の相転移検出
A.Tanaka and Y. Tomita, J. Phys. Soc. Jpn.
86, 063001 (2017)

スピン分布の解析

P. Broecker et al. Scientiﬁc Reports 7 8823 (2017)
volume7

ハニカムハバード模型

同時刻グリーン関数の解析

多くの論文が画像認識技術をベースにしている11

機械学習の物理学への応用
画像認識的な方法：相転移の検出等
別のアプローチはないか？
自然界は複雑すぎる
単純なモデルを作って解析する
例：ボールを投げたらどこに落ちる？
風は吹いていない「とする」 -> まあそんなに間違っていない結果
Done is better than perfect

物理学者が長年やってきた方法： 現象を記述する簡単なモデルを作る

機械が簡単なモデルを作ることも可能では？

12

分子動力学法 (MD)
原子や分子に働く力を仮定して時間発展を追うシミュレーション

r2
r1

Lennard-Jones ポテンシャルetc.

F2F1
r3

F3

力はエネルギーを偏微分すれば得られる∂E(r1, ⋯, rN )∂E(r1, ⋯, rN )∂E(r1, ⋯, rN )
F1 =
F2 =
F3 =
∂r1
∂r2
∂r3
正確なエネルギーがあれば正確な力を計算できる
原子分布{R}の関数としてのエネルギーE({R})

第一原理計算：多体系における高精度なエネルギー計算
密度汎関数理論:Hohenberg‒Kohn theorems

3N次元波動関数と3次元電子密度が一対一対応

多電子系の基底状態の性質が3つの空間座標だけに依存する電子密度によって一意に決定される固体の電子状態の計算が可能に！

時間ステップごとに第一原理計算して力を求める
非常に計算が重い！

機械学習を使えないか？

13

機械学習をどう使う？
インプット
原子の配置{ri}

アウトプット
第一原理計算etc

r2

r1

ポテンシャルエネルギーE({ri})
原子間に働く力∂E(r1, ⋯, rN )
F1 =
∂r1

r3

…

インプットを入れればアウトプットが出てくるニューラルネットを作ればいいのね

ニューラルネットワーク

そう簡単ではない14

機械学習をどう使う？
第一原理分子動力学で機械学習が使いたい時

学習で作ったニューラルネット
従来では重すぎて実行できないようなは原子数の異なるより大きな系大きな系の計算を行いたい時で使いたい！

例：H2O 32個で学習 -> H2O 1000個で実行

インプット

アウトプット…

どんなインプット、どんな形状のニューラルネットなら可能か？

Behler-Parrinello型ニューラルネットワーク
第一原理計算：高精度な力の評価->高精度シミュレーション
時間ステップ1fsのMD -> 1nsのシミュレーションに106回の第一原理計算が必要

とても重たい -> ニューラルネットでフィットしたい“Neural network models of potential energy surfaces”
T. B. Blank et al., J. Chem. Phys. 103, 4129 (1995)

解決すべき問題
1. エネルギーE(r1,r2,…,rN) は3N次元空間の関数で複雑
2. Portability: 作ったニューラルネットを他で使えるか？
3. 原子が持つ対称性（回転対称性、並進対称性）を持つか？

2007 (13年前)“Generalized Neural-Network Representation of High-Dimensional Potential-EnergySurfaces”, J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007)

原子配置 -> 第一原理計算 -> エネルギー
原子配置 -> ニューラルネット -> エネルギー

1000倍以上の高速化

16

Behler-Parrinello型ニューラルネットワーク

基本アイディア
1. インプットは原子ごとに局所的な位置情報だけから作る：原子ごとの”指紋”
2. 同じ原子種は同じニューラルネットワークを使う

原子数増減に対応可能

Behler-Parrinello型ニューラルネットワーク(r1,r2,r3)の代わりに (G1(r1,r2,r3),…,G3(r1,r2,r3)) をインプットに使う

S1: 原子周りでのニューラルネットワーク
回転対称性や並進対称性が自動的に考慮されている

Behler-Parrinello型ニューラルネットワーク
教師あり学習：インプットとアウトプット

Schematic figure from Y. Ando, Deep-learning and physics 2018 @Osaka Univ.

同じ原子種には同じニューラルネットワークを使用
Gの数が固定-> 原子数が増えても同じニューラルネットで対応が可能

オープンソースソフトウェアaenet
LDA+U 16,047 LiMO2

N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96, 014112 (2017).

インプットベクトルの長さ（周辺環境を記述するGの数）とニューラルネットワークの精度

隠れ層2層(各15ユニット)のニューラルネットワーク

エネルギーの精度: 3mev/atom

The Atomic Energy Network (ænet)http://ann.atomistic.net
近年：オープンソースソフトウェアが充実
2007年から10年近く、弟子入りしないと使えないような状態だった…

20

自己学習ハイブリッド
モンテカルロ法

機械学習MDの改善点
訓練データをどう集めるか？
2通り

1. つぶしたりのばしたりして構造を系を沢山作る
色々な原子位置に対してエネルギーを第一原理計算->教師データを集める
問題点：どういう系を集めてくれば良い？（配置は無限にある）
問題点：何個集めれば十分に精度の良いニューラルネットができる？
1000個？10000個？100000個?

2. 第一原理分子動力学計算を実行して、構造を集める
1ステップごとに訓練データが増える
問題点：何個集めれば十分に精度の良いニューラルネットができる？
1000ステップ？10000ステップ？100000ステップ?

問題点：これを高速化したいのに…。時間がかかりすぎる機械学習MDの改善点作ったニューラルネットワークの精度は？
悪いニューラルネットワークは悪い力を返す

より大きな系に適用したい！
その系でも精度がいい保証は？

多成分系(4元素以上)に弱い
入力データが複雑になり特徴が捉えづらくなる->悪いポテンシャルにニューラルネットワークの質に踊らされたくない！
->自己学習モンテカルロ法を使おう23

自己学習モンテカルロ法
モンテカルロ法を機械学習を使って高速化する手法“物理学者,機械学習を使う ー機械学習・深層学習の物理学への応用ー ”
朝倉書店(2019/10)
第5章”自己学習モンテカルロ法”
Learn（オリジナル模型を使用）

Earn（SLMC法による本シミュレーション）xi
x1 = (1,-1,1,1,-1,1,1)

xj

H(xi)

H(xj)

Heff

x2 = (1,-1,1,-1,-1,1,1)

採択する？

…

機械学習xi+1xi+2xi
xi = (1,-1,1,-1,-1,1,1)

Heff(xi)

xj-1
xj

提案マルコフ連鎖

Heff(xj)

機械学習で有効模型を作り、計算量を削減
固体物理、特にスピン系や電子系の模型に適用されてきた結果は統計的に厳密に元のモンテカルロ法に等しい
機械学習MDにも使えるのでは？

24

SLMCと機械学習MDの比較
自己学習モンテカルロ法(SLMC)
モンテカルロ配位C
スピン配位etc

モンテカルロ重み W(C)

機械学習分子動力学法
原子位置 R=(r1,r2,r3,…)
第一原理計算によるエネルギー E(R)

有効ハミルトニアン

ニューラルネットででフィッティング

フィッティング

有効モンテカルロ重み Weﬀ(C)有効模型は配位提案にのみ使う
厳密！

有効エネルギー Eeﬀ(R)
有効模型のMDで時間発展：近似MCを使おう！

厳密!

25

ハイブリッドモンテカルロ法
分子動力学法もモンテカルロ法も統計平均をとるという意味では同じ分配関数と物理量の期待値を計算する手法

ハイブリッドモンテカルロ法

格子QCDで使われている

S. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, and R. L. Sugar, Hybrid-molecular-dynamics algorithms forthe numerical simulation of quantum chromodynamics, Phys. Rev. D 35, 2531 (1987).

固体物理でも考案されている
B. Mehlig, D. W. Heermann, and B. M. Forrest, Hybrid Monte Carlo method for condensed-matter systems,
Phys. Rev. B 45, 679 (1992).

初期原子配置{Ri}MDで十分だったからかあまり使われてこなかった

提案原子配置{Riʼ} {Riʼ}を受け入れる？

機械学習MDで時間発展
MDで時間発展

メトロポリス法に従って受け入れ確率を計算

MD+MC = ハイブリッドMCYN et al., arxiv:1909.02255

自己学習ハイブリッドモンテカルロ法(SLHMC)
機械学習ポテンシャルを利用して、第一原理MDと統計的に厳密に等しい結果を得ることができる

機械学習=近似計算
ハイブリッドモンテカルロ法の一種

ではない

原子の運動量{p}と位置{r}をモンテカルロ法で変更

1. Trial conﬁgurationを機械学習MDで生成

第一原理計算第一原理計算第一原理計算

(DFT)

(DFT)

(DFT)

機械学習MD

2. 提案した{pʼ},{rʼ}を以下の確率で受け入れる

機械学習MDDFTのエネルギーをボルツマン重みとした
メトロポリスヘイスティングス法

機械学習ポテンシャルは配位提案に使う->第一原理HMCと結果が厳密に等しい自己学習ハイブリッドモンテカルロ法(SLHMC)
1. Trial conﬁgurationを機械学習MDで生成

2. 提案した{pʼ},{rʼ}を以下の確率で受け入れる

第一原理計算第一原理計算第一原理計算

(DFT)

(DFT)

(DFT)

機械学習MD

機械学習MD

MDの時間ステップΔt が微小ならばPacc(p, r) → {p′, r′}) ∼ min (1, exp [−β(VDFT({r′}) − VML({r′}) − (VDFT({r}) − VML({r})))])
ポテンシャルの差

ポテンシャルの差

機械学習ポテンシャルが完全に第一原理計算のポテンシャルを模倣すると、アクセプト率は1となる

YN et al., arxiv:1909.02255

自己学習ハイブリッドモンテカルロ法(SLHMC)
機械学習ポテンシャルを利用して、第一原理MDと統計的に厳密に等しい結果を得ることができる

機械学習=近似計算

ではない

1. 学習に使うデータを自動的効率的に収集第一原理MDと同じアンサンブルを出力 -> 効率的なデータ集め
機械学習MDを行う目的の一つは第一原理MDの代わりに計算することターゲットの第一原理MDと同じアンサンブル -> 欲しいデータそのもの
大量の学習データをあらかじめ集める必要がない
ポテンシャルはon-the-ﬂyで逐次的に改善されていくYN et al., arxiv:1909.02255

自己学習ハイブリッドモンテカルロ法(SLHMC)
機械学習ポテンシャルを利用して、第一原理MDと統計的に厳密に等しい結果を得ることができる

機械学習=近似計算

ではない

2. 作ったポテンシャルは近似でも結果は厳密+トレーニングしながら物理量計算
第一原理計算第一原理計算第一原理計算第一原理計算第一原理計算

(DFT)

(DFT)

(DFT)

(DFT)

(DFT)

機械学習MD

機械学習MD

機械学習MD

機械学習MDon-the-ﬂy
学習データセット

学習データセット

学習データセット

学習データセット

新しい機械学習ポテンシャルを生成

改善

30

分子シミュレーション手法の比較
計算精度

計算コスト

古典分子動力学低い

非常に低い

第一原理分子動力学

非常に高い

非常に高い

機械学習分子動力学

高い（学習結果に依存） 低い

自己学習ハイブリッドモン 非常に高い（学習結果に 中間
テカルロ法依存しない）

要求される精度と計算コスト（サイズ、原子数）で手法を選ぶ

SLHMCは二つの利用法がある
1. シミュレーション：第一原理分子動力学法ができない領域で精度保証で計算したい
2. ニューラルネット作成：機械学習分子動力学用のニューラルネットを高精度に構築したい

どちらにせよ、ニューラルネットは自動的に構築される->手間も大幅減31

実際にシミュレーションしたい方へ
第一原理分子動力学パッケージ

PIMD

https://ccse.jaea.go.jp/software/PIMD/index.en.html
経路積分第一原理分子動力学（水素の量子効果考慮等）
も可能
VASP,QuantumEspresso等様々な第一原理計算パッケージに対応

機械学習MDが可能

The Atomic Energy Network (ænet)自己学習ハイブリッドモンテカルロ法が実装済み！
実習開催の実績あり（東京大学2019年）

関連記事

PIMDとAENETを使って機械学習分子シミュレーションをやってみよう

https://qiita.com/cometscome_phys/items/660a84119b5d4acce903

適用例

YN et al., arxiv:1909.02255

α-quartz (SiO2)
α-quartz (SiO2)

2x2x2スーパーセル 300K NVT
第一原理計算：VASP
機械学習ポテンシャル：artiﬁcial neural networks (ANN)The Atomic Energy Network (ænet)ニューラルネットワークは標準的な構造(15ユニットの隠れ層2層)ハイブリッドモンテカルロ法：
PIMD: An Open Source Software for Parallel Molecular Simulation
https://ccse.jaea.go.jp/ja/download/pimd/index.en.html
10

1. 第一原理MDで300構造集めて初期ANNを作成∆V/NAtoms [meV/atom]

86

2. 0.25fsステップで機械学習MD4

3. 5fsごとにで第一原理計算してアクセプト率計算
し、学習データ(構造とエネルギー)も収集

20

4. 同時に物理量の計算も行う−2

5. 500fsごとにでANNをアップデート−40200

400

600

800

1000

1200

Num. of DFT calc.

1400

1600

ANNポテンシャルの精度が改善されていく->1meV/atom以下に

YN et al., arxiv:1909.02255

α-quartz (SiO2)
α-quartz (SiO2)

2x2x2スーパーセル 300K NVT
第一原理計算：VASP
機械学習ポテンシャル：artiﬁcial neural networks (ANN)The Atomic Energy Network (ænet)ニューラルネットワークは標準的な構造(15ユニットの隠れ層2層)ハイブリッドモンテカルロ法：
PIMD: An Open Source Software for Parallel Molecular Simulation
https://ccse.jaea.go.jp/ja/download/pimd/index.en.html
10-78500

DFT [K/atom]∆V/NAtoms [meV/atom]

8-78550
6-78600

42-78650
0-78700

−2−40

200

400

600

800

1000

1200

Num. of DFT calc.

1400

1600

ANN [K/atom]-78750-78750 -78700 -78650 -78600 -78550 -78500

YN et al., arxiv:1909.02255

α-quartz (SiO2)
α-quartz (SiO2)

2x2x2スーパーセル 300K NVT
第一原理計算：VASP
機械学習ポテンシャル：artiﬁcial neural networks (ANN)The Atomic Energy Network (ænet)ニューラルネットワークは標準的な構造(15ユニットの隠れ層2層)ハイブリッドモンテカルロ法：
PIMD: An Open Source Software for Parallel Molecular Simulation
https://ccse.jaea.go.jp/ja/download/pimd/index.en.html
35

DFT−MD: Si−Si
SLHMC: Si−Si
DFT−MD: Si−O
SLHMC: Si−O
DFT−MD: O−O
SLHMC: O−O

30

RDF

25

動径分布関数（原子の位置の分布）
第一原理MDの結果を完全に再現している
DFTのエネルギーをボルツマン重みとした
メトロポリスヘイスティングス法

201510

機械学習ポテンシャルは配位提案に使う

5->第一原理HMCと結果が厳密に等しい

00123r [Å]456

注：トレーニングしながら計算が可能

YN et al., arxiv:1909.02255

α-quartz (SiO2)
α-quartz (SiO2)

2x2x2スーパーセル 300K NVT
第一原理計算：VASP
機械学習ポテンシャル：artiﬁcial neural networks (ANN)The Atomic Energy Network (ænet)ニューラルネットワークは標準的な構造(15ユニットの隠れ層2層)ハイブリッドモンテカルロ法：
PIMD: An Open Source Software for Parallel Molecular Simulation
https://ccse.jaea.go.jp/ja/download/pimd/index.en.html
8平均自乗変位: 位置相関
(固体では一定値に収束する)

7MSD [×10−3 Å2]

6

機械学習MDの長さを変更

543
SLHMC: ∆tDFT = 5fs
SLHMC: ∆tDFT = 12.5fs
SLHMC: ∆tDFT = 25fs
SLHMC: ∆tDFT = 50fs
SLHMC: ∆tDFT = 250fs
DFT−MD: ∆t = 1fs
DFT−MD: ∆t = 0.25fs

2100100

200300

Num. of DFT calc.

400

機械学習MD

DFT計算の回数を減らすことができる
500

高速化

YN et al., arxiv:1909.02255

YNi2B2C超伝導体(Tc ~ 15K)

フォノン起源の非従来型超伝導体(強い異方性をもつ)

中性子散乱実験：フォノンDOSに強い温度依存性->非調和振動のフォノンが重要な系
第一原理計算でフォノンの温度依存性を出すには？

->第一原理分子動力学法
フォノンを計算するには長時間シミュレーションが必要->第一原理MDでは厳しい機械学習MDを使えばフォノンDOSを計算できる？

->4元素系->安定な機械学習MDは難しい

2x2x2

SLHMCで「良いポテンシャル」を作成し、機械学習MDをやってみる

YN et al., arxiv:1909.02255

YNi2B2Cフォノン起源の非従来型超伝導体(強い異方性をもつ)超伝導体(Tc ~ 15K)
7.0
6.0

Phonon DOS [Arb. unit]

1000KでSLHMC
well-trained ANNを作成

YNi
BC
Total

5.0
4.060K

3.0

1step 1fsでANN-MD

2.0

10万ステップ->100ps

1.0
0.0

0

20406080100120140160

Frequency [meV]
7.0Y
Ni
BC
Total

Phonon DOS [Arb. unit]

6.0
5.0

フォノンDOSが強い温度依存性を持つ

4.0

300K

3.02.0通常のDFT(Phononのバンド分散

1.02x2x20.0

10回繰り返し、計算を計算してDOS)では計算できない
020406080

100

Frequency [meV]

120140160

効果39まとめYN et al., arxiv:1909.02255

まとめ
自己学習ハイブリッドモンテカルロ法
1. 機械学習ポテンシャルを利用して、第一原理MDと統計的に厳密に等しい結果を得ることができる
第一原理計算第一原理計算第一原理計算

(DFT)

(DFT)

(DFT)

機械学習MD

機械学習MD

ポテンシャルはon-the-ﬂyで逐次的に改善されていく
トレーニングしながら計算が可能

DFTのエネルギーをボルツマン重みとした
メトロポリスヘイスティングス法

DFT計算の回数を減らすことができる

機械学習ポテンシャルは配位提案に使う

高速化->第一原理HMCと結果が厳密に等しい

2.学習で得られた機械学習ポテンシャルは、ターゲットMDからサンプルされたもの->well-trained なポテンシャル
得られたポテンシャルで長時間安定にML-MDができる