PIMDを用いた機械学習分子シミュレーション
精度保証された機械学習シミュレーション手法：自己学習ハイブリッドモンテカルロ法の紹介
機械学習ポテンシャルの精度に依らずに第一原理計算精度を保証する計算手法
原子力機構システム計算科学センターシミュレーション技術開発室
理研
1

革新知能統合研究センター (AIP)

永井佑紀

2

アウトライン
イントロダクション
原子・分子シミュレーションと機械学習
機械学習分子シミュレーションの問題点と解決策：自己学習ハイブリッドモンテカルロ法
PIMDによる原子分子シミュレーション：SiO2とYNi2B2Cを例にしてまとめ

3

イントロダクション

4

今日の講演
機械学習を何に使うか？
猫画像判別、自動運転、囲碁、自動お絵描きetc
何かを判断したり、推測したりを機械が行う -> 代わりにやってもらうA

B:人間には定式化できない何かシミュレーションに機械学習を使うには？
これまで計算が大変だった部分を機械が行う

A

B

定式化済みだけど計算コスト大A

B
計算コストの軽い機械

シミュレーションの高速化に機械学習を用いる

5

今日の講演
シミュレーションに機械学習を使うには？
これまで計算が大変だった部分を機械が行う

AA

B

定式化済みだけど計算コスト大B
計算コストの軽い機械

シミュレーションの高速化に機械学習を用いる
物質の情報x

何らかのアウトプットyf(x):場の量子論、密度汎関数理論、流体力学etc

物質の情報x

何らかのアウトプットyg(x): ニューラルネットワーク

機械学習を使ったシミュレーションの妥当性はどう検証するのか？
-> 自己学習ハイブリッドモンテカル
ロ法の提案

6

今日の講演
PIMD：CCSEが開発している原子分子シミュレーション
パッケージ(メイン開発者：志賀基之氏)
オープンソース:無料

原子同士に働く力を計算する様々な方法に対応している
古典ポテンシャル、第一原理計算、機械学習ポテンシャル
量子効果も考慮できる（経路積分MD）ので、水素、重水素、三重水素の違いも扱える

自己学習ハイブリッドモンテカルロ法が実装ずみ
1. 機械学習ポテンシャルをシミュレーションしながら学習できる
液体、固体、合金、様々な系で「最適な」ポテンシャルを構築できる

2. 機械学習ポテンシャルの精度に依らずに、第一原理MDの精度を達成できる

7

今日の講演
PIMD：CCSEが開発している原子分子シミュレーション
パッケージ(メイン開発者：志賀基之氏)•••

アンサンブル

構造最適化、基準振動解析、フォノン計算、弾性定数•

反応経路探索：ストリング法、最急降下法、最緩上昇法

エネルギー一定（NVE）、温度一定（NVT）、圧力一定（NPT）、応力一定（NtT）

境界条件

古典統計：レプリカ交換法••

量子統計：経路積分分子動力学法•

古典力学：分子動力学法、ハイブリッド・モンテカルロ法•

半古典動力学：セントロイドおよびリング・ポリマー分子動力学法

自由境界、周期境界（立方セル、平行六面体セル）

ポテンシャル•

第一原理的手法：SMAS•

古典力場：ADP(new!)、AMBER、CHARMM、CLAYFF(new!)、EAM、OPLS、•

非断熱動力学：平均場動力学法、サーフェス・ホッピング法

Tersoff(new!
• 分極力場：誘起双極子とTholeダンピング補正 (new!

•自由エネルギー計算：メタ・ダイナミクス法、拘束分子動力学•C)

M

TURBOMOLE、VAS

T

機械学習：自己学習・ハイブリッドモンテカルロ法(new!)H

GAMESS、GAUSSIAN、MOLPRO、NTCHEM、ORCA、QUANTUM ESPRESSO、P

法、平均力ダイナミクス)•

他のコードとのインターフェース（第一原理的手法）：ABINIT-MP、CP2K(new!)、•

他のコードとのインターフェース（半経験的手法）：DFTB、MOPA•

マルチスケール法：ONIOM、QM/M•

機械学習ポテンシャル：AENE•

ユーザー定義のサブルーチン

8

物理学と機械学習は相性が良い
物理学では何をするか観測する

機械学習では何をするか法則を見つける

ニューラルネットワークを作る

ニュートンの法則？

y

空気抵抗を無視してみる

猫?

大体合う？

yes or no
データに合う箱を作る

もっと精緻に考える？
形状も考える？
ボールの落下地点を理論を構築し、記録する

解いて予測データを集める

箱を作り、入れて予測9

機械学習をどう使うか
機械学習

学習によって作られたニューラルネットワーク

猫?
yes or no

人間は「猫を識別する理論」を知らないが「猫を識別する関数」は作れる知っている理論に適用すると何が起こる？

物理学で、「インプットxから何かを計算する理論f(x)」を知っているとする「インプットxから何かを計算する関数g(x)」をニューラルネットで作る
f(x)=g(x)となるg(x)を作ると何が嬉しいか？

10

機械学習をどう使うか物理学で、「インプットxから何かを計算する理論f(x)」を知っているとする「インプットxから何かを計算する関数g(x)」をニューラルネットで作る
f(x)=g(x)となるg(x)を作ると何が嬉しいか？
物理学
違いは何か
物質の情報x

何らかのアウトプットyf(x): 非常に計算が大変な場合がある

f(x):場の量子論、密度汎関数理論、流体力学etc

物質の情報x

何らかのアウトプットyg(x): ニューラルネットワーク

g(x): 学習が終われば計算は一瞬計算の高速化に使うことができる！

11

原子・分子シミュレー
ションと機械学習

12

量子力学と物質
我々は量子力学を知っている⃗⃗⃗⃗̂Hψ( r ⃗ , ⋯, r ⃗ , R , ⋯, R ) = Eψ( r ⃗ , ⋯, r ⃗ , R , ⋯, R )1N

1M

1N1M

注：時間依存がない場合

電子位置rと原子核の位置Rに関する多体のシュレーディンガー方程式を知っている

これが解ければ物質に関する殆どのことがわかる（はず）

密度汎関数理論の登場(1998年ノーベル化学賞Walter Kohn)によってE⃗
R 1⃗ ,⋯, R Mは計算できるようになった（原子核位置Rを固定し電子の量子性を考慮）

R2
R1

F1R3

あらゆるものは原子からできている

F2
F3

固体、液体、気体、電池、生体、創薬、ウイルスetc.
計算ができればなんでもシミュレーションできる?
-> 計算が重たい!

分子動力学法 (MD)
原子や分子間に働く力がわかればシミュレーションできる

r2
r1

エネルギーから力は計算できる

F2∂E(r1, ⋯, rN )∂E(r1, ⋯, rN )∂E(r1, ⋯, rN )
F1 =
F2 =
F3 =
∂r1
∂r2
∂r3

F1r3

Lennard-Jonesポテンシャル等

F3

精度の良いエネルギー評価ができれば精度の良い力がわかる

第一原理計算 -> 精度の良いエネルギーの評価が可能
密度汎関数理論(DFT)

各時刻ごとにDFT計算を行えば、精度の良い力が計算できる

典型的には、1fsごとにDFT計算を行う
f -> p -> n -> μ -> m -> 1s

ニューラルネットに置き換えよう14

Behler-Parrinello型ニューラルネットワーク
第一原理計算 -> 高精度な力の評価 とても重たい -> ニューラルネットワークでフィット“Neural network models of potential energy surfaces
T. B. Blank et al., J. Chem. Phys. 103, 4129 (1995)

2007 年 革新的な論文が登場“Generalized Neural-Network Representation of High-Dimensional Potential-Energy
Surfaces”, J. Behler and M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007)

原子位置DFT計算

エネルギー

原子位置ニューラルネットワーク

エネルギー”

1000倍以上高速に

15

やってみたい！そんな時には原子・分子シミュレーション用のニューラルネットワークを作るソフトウェア

MDを行うソフトウェア

オープンソースソフトウェア

LAMMPS

The Atomic Energy Network (aenet)

A neural network potential package (n2p2)

PIMD：我々（原子力機構）が開発

市販ソフトウェア

ニューラルネットワーク分子動力学システム

Advance/NeuralMD

汎用原子レベルシミュレータMatlantis

新たなプログラミングなしに機械学習分子シミュレーションができる環境が整っている16

機械学習シミュレーションの問題点と解決策：
自己学習ハイブリッドモンテカルロ法

17

機械学習の(本質的な？)問題点
機械学習は何かを模倣する

g(x) ∼ f(x)与えられたインプットから、オリジナルによく似たアウトプットを出力する
機械学習をシミュレーションに使う場合ちゃんと計算できる部分を機械学習に置き換えた事による誤差は？
物質の情報x

何らかのアウトプットyf(x):場の量子論、密度汎関数理論、流体力学etc

物質の情報x

何らかのアウトプットyg(x): ニューラルネットワーク

機械学習を使ったシミュレーションの妥当性はどう検証するのか？
完全に一致することはないが、どこまで精度を高めれば良いのか？

機械学習MDの問題点について訓練データはどうやって集めるのか？
1. 圧縮したり引き伸ばしたり歪めたりして作った色々な構造のデータを集める
様々な構造で第一原理計算を行う ->訓練データが集まる
問題点 1 : どんな構造を入れるべき？（構造候補は無限にある）
問題点 2 : どのくらいの数の訓練データを集めれば、やりたいシミュレーションを再現するに足る十分な数だと結論づけられるのか？

1000？10000？100000?

2. 第一原理MDを行なって得られる構造データを集める
毎ステップごとに新しい構造データが得られる
問題点 1 : どのくらいの長さのMDを行えばいいのか？どのくらいの数の訓練データを集めればやりたいシミュレーションを再現するに足る十分な数だと結論づけられるのか？
1000steps? 10000steps? 100000steps?

問題点 2: そもそも第一原理MDの計算コストが大きいために機械学習MDを使いたいわけで第一原理MDでデータを集めていたら時間がかかってしまう機械学習MDの問題点について構築したニューラルネットワークの精度はどうやって評価するのか
機械学習ポテンシャルはDFTポテンシャルの近似

悪いニューラルネットワークは悪い力を生む

大きな系を取り扱いたい
大きな系は教師データがないが、使っているNNの精度をどうやって知ればいいのか？
(第一原理MDができないような大きな系で使いたいので、教師データがない)

機械学習MDは原子種の数が増える（4原子種以上）と不安定になりがち
複雑なインプットから作ったニューラルネットワークの精度は悪くなりがち
ニューラルネットワークの精度の良し悪しにに振り回されたくない！

->自己学習ハイブリッドモンテカルロ法を使おう20

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

機械学習MD

機械学習MD

エネルギー計算

エネルギー計算

エネルギー計算−βH({r},({p})−βH({r},({p})−βH({r},({p})

eee従来のHMC: MD部分は第一原理MD
自己学習HMC: MD部分を機械学習MDに置き換えておく

MLMDとMCを合わせた手法 自己学習HMCの計算精度は第一原理MDと等しい
ニューラルネットワークを使っているのに、計算精度が元の手法と同じになっている！
ニューラルネットワークの置き換えによる高速化と、第一原理計算の精度を両立

YN et al., Phys. Rev. B 102, 041124 (2020)

SLHMC法の計算について
機械学習MD: 近似されたポテンシャルを利用 ->厳密ではないSLHMC: 近似されたポテンシャルを利用-> それでも厳密！
機械学習MD

機械学習MDDFT

DFTMetropolis-Hastings法の一種なので機械学習MD

DFT

on-the- yに改善
訓練データ

訓練データ

訓練データ

新しい機械学習ポテンシャルを構築

訓練データはターゲット
アンサンブルから取得

良い教師

fl

機械学習ポテンシャルがシミュレーション中に逐次的に改善される

YN et al., Phys. Rev. B 102, 041124 (2020)

自己学習ハイブリッドモンテカルロ法 (SLHMC)
1. 厳密: 第一原理MDと同じ精度を持つが高速に計算ができる
2. On-the- y: 機械学習ポテンシャルはシミュレーションの最中に自動的に改善
される

DFTポテンシャルに対するメトロ
ポリスーヘイスティングス法
ハイブリッドモンテカルロ法の一種

統計的に厳密な結果
以後の計算には我々が開発したPIMDを使用

PIMD

https://ccse.jaea.go.jp/software/PIMD/index.en.html

fl

ニューラルネットにはaenetを使用

23

PIMDによる原子分子シミュレーション

様々な分子シミュレーション手法
Accuracy

Computational cost

Classical MDlow

very lowFirst-principles MDvery highvery highMachine-learning MD

highlow

Self-learning HybridMonte Carlovery highmiddleSLHMCを使うと嬉しい時：
1. DFTレベルの精度の計算をしたいがDFT-MDの計算コストが大きすぎて断念してしまっている時
2. 機械学習MDのための高精度なニューラルネットワークを構築したい時

どちらの場合でも、ニューラルネットワークは自動的に改善される

YN et al., Phys. Rev. B 102, 041124 (2020)

逐次的訓練と厳密な結果
α-quartz (SiO2)
2x2x2 スーパーセル 300K NVT
DFT：VASP
ニューラルネットワーク：隠れ層2層、各15ユニット
PIMD with aenetを使用-785008DFT [K/atom]∆V/NAtoms [meV/atom]

10-78550
64

ニューラルネットワークが2自動的に改善されていく-78600-786500-78700−2−40200

400

600

800

10001200

1400

1600

Num. of DFT calc.

1meV/atomの精度

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

YN et al., Phys. Rev. B 102, 041124 (2020)

逐次的訓練と厳密な結果
α-quartz (SiO2)
2x2x2 スーパーセル 300K NVT
DFT：VASP
ニューラルネットワーク：隠れ層2層、各15ユニット
PIMD with aenetを使用

SLHMCはDFTポテンシャルに対する35DFT−MD: Si−Si
SLHMC: Si−Si
DFT−MD: Si−O
SLHMC: Si−O
DFT−MD: O−O
SLHMC: O−O

30RDF

25

Metropolis-Hastings法
1. 訓練データ: 対象とするアンサンブルそのもの

20

効率的なサンプリングと訓練が可能

1510

2. ニューラルネットワークの精度に依らず、常に厳密

5001234r [Å]動径分布関数56

結果の精度が保証されている！

YN et al., Phys. Rev. B 102, 041124 (2020)

逐次的訓練と厳密な結果
α-quartz (SiO2)
2x2x2 スーパーセル 300K NVT
DFT：VASP
ニューラルネットワーク：隠れ層2層、各15ユニット
PIMD with aenetを使用

計算コストがDFT-MDと比べて劇的に減少

87MSD [×10−3 Å2]6543
SLHMC: ∆tDFT = 5fsSLHMC: ∆tDFT = 12.5fs
SLHMC: ∆tDFT = 25fs
SLHMC: ∆tDFT = 50fs
SLHMC: ∆tDFT = 250fs
DFT−MD: ∆t = 1fs
DFT−MD: ∆t = 0.25fs

2100100

200

300

Num. of DFT calc.

平均二乗変位400

機械学習MD

DFT計算の回数が減っている
500

YN et al., Phys. Rev. B 102, 041124 (2020)

高精度な機械学習ポテンシャルによる長時間MD
SLHMCで機械学習ポテンシャルを訓練

YNi2B2C:
フォノン由来の非従来型s波超伝導体(Tc ~ 15K)
DFT: VASPPIMD with aenetを使用
1000K HMC -> 高精度な機械学習ポテンシャルを構築
1fs 間隔で100,000 ステップで MD x 10
通常のDFT-MDではこの長さの計算は大変
2x2x2

YN et al., Phys. Rev. B 102, 041124 (2020)

高精度な機械学習ポテンシャルによる長時間MD
SLHMCで機械学習ポテンシャルを訓練
フォノン状態密度
7.06.0Phonon DOS [Arb. unit]

1. SLHMCで訓練

YNi
B
CTotal5.0
4.0

2. ML-MDを実行

60K

3.0

1000K HMC -> 高精度な機械学

2.0
1.0
0.00

20406080

100120140160

Frequency [meV]7.0Phonon DOS [Arb. unit]

1fs 間隔で100,000 ス
テップで MD x 10YNi
B
C
Total

6.0
5.0
4.0

300K

3.02.0強い温度依存性を発見

1.02x2x20.0

習ポテンシャルを構築

020406080

100

Frequency [meV]

120140160

実際の計算時間
SLHMC: SiO2

240CPUコア並列：JAEAスパコンを使用

SLHMC+MD: YNi2B2C

SLHMCは訓練時間込みでも速い！

K. Kobayashi, YN, M. Itakura, and M. Shiga, J. Chem. Phys. 155, 034106 (2021)

液体シリカ
NVTアンサンブル（粒子数、体積、温度一定） 以外でも使いたい
NPTアンサンブル（粒子数、圧力、温度一定）の手法を開発
72原子、4000Kで、50ps後の構造

MLMD
8600構造で学習->結晶構造が壊れる

SLHMC
悪い構造はreject
MD部分を可変にしたしているため安定約200fsに一回だけ第一原理計算を行う
第一原理計算精度を保った1ns オーダーのシミュレーションが可能に->安定なシミュレーションが可能

K. Kobayashi, YN, M. Itakura, and M. Shiga, J. Chem. Phys. 155, 034106 (2021)

液体シリカ
NVTアンサンブル（粒子数、体積、温度一定） 以外でも使いたいNPTアンサンブル（粒子数、圧力、温度一定）の手法を開発

80 psの第一原理MDの結果との比較 (4000K)

ポテンシャルが完全に一致していなくても、物理量は同じ精度

K. Kobayashi, YN, M. Itakura, and M. Shiga, J. Chem. Phys. 155, 034106 (2021)

液体シリカ
NPTアンサンブル（粒子数、圧力、温度一定）の手法を開発
SLHMC

平均二乗変位：三つの時間領域がある
ballistic, plateau and di usive regimes
実験での融点は1983K
DFT-MDでは3500Kでもdi usive regimeがよく見えない

ff

SLHMCだとMDを長くすることで実効的に長時間平均が計算可能 -> diffusiveregimeが観測できる

ff

DFT-MDK. Kobayashi, YN, M. Itakura, and M. Shiga, J. Chem. Phys. 155, 034106 (2021)

液体シリカ
NPTアンサンブル（粒子数、圧力、温度一定）の手法を開発

平均二乗変位：三つの時間領域がある
ballistic, plateau and di usive regimes
実験での融点は1983K
DFT-MDでは3500Kでもdi usive regimeがよく見えない
SLHMCだとMDを長くすることで実行的に長時間平均が計算可能 -> diffusiveregimeが観測できる-> 実験値をよく再現する結果が得られる
ff

ff

x線による構造因子

35

まとめ

36

まとめ
精度保証された機械学習シミュレーション手法：自己学習ハイブリッドモンテカルロ法の紹介
機械学習ポテンシャルの精度に依らずに第一原理計算精度を保証する計算手法

PIMD：CCSEが開発している原子分子シミュレーションパッケージ(メイン開発者：志賀基之氏)

オープンソースソフトウェア
機械学習をシミュレーションに使う不安模倣なので結果が正しいか不安 -> 自己学習ハイブリッドモンテカルロ(SLHMC)法SLHMCを使うと嬉しい時：

YN et al., Phys. Rev. B 102, 041124(R) (2020)

1. DFTレベルの精度の計算をしたいがDFT-MDの計算コストが大きすぎて断念してしまっている時
2. 機械学習MDのための高精度なニューラルネットワークを構築したい時