2023年10⽉31⽇

2023年度第1回計算科学フォーラム

富岳を⽤いた⾼度計測むけデータ駆動科学：2次元物質解析への応⽤
核融合研究所，高エネ研 (客員), 東大物性研(客員) 星健夫（出⾝：東京⼤学⼯学部物理⼯学科）

・背景1：21世紀における計算科学のパラダイムチェンジ
・背景 2：⾼度計測におけるデータ解析によるブレークスルー
・計測データ解析フレームワーク「2DMAT」
・超並列モンテカルロ法によるベイズ推定
・全反射⾼速陽電⼦回折(TRHEPD)への応⽤
・まとめ・⼤局的展望：デジタルツインなど

謝辞：和田 健・望月出海（高エネ研低速陽電子）,
吉見一慶・本山裕一（東大物性研），福島孝治（東大総合文化），花⽥貴（東北⼤⾦研），
⾼草⽊達・朝倉清⾼（北⼤触媒研），⼭本有作・⼯藤周平（電通⼤），岩本晴道，⽊下直希（⿃取⼤） ，(*)謝辞
中野陽⽃，寺地雄真（核融合研）
・富岳は，HPCI課題(hp210267, hp220248, hp230304)などで実施した．
・ 2DMAT(https://github.com/issp-center-dev/2DMAT/)は，「2020・2021年度 東⼤物性研 ソフトウェア開発・⾼度化プロジェクト」などの⽀援を受け開発・公開され，東⼤物性研スーパーコンピュータにプリインストールされている．
・2DMATの⼀部データは，以下ディレクトリにアップロードされている；
2DMAT Gallery （https://datarepo.mdcl.issp.u-tokyo.ac.jp/repo/17）

背景 1 ：21世紀における計算科学のパラダイムチェンジ
21世紀

20世紀
コンピュータの誕生→数値計算アルゴリズムの整備

パラダイムチェンジ（古いプログラムは「死ぬ」？）

計算機の質的変化：
超並列コンピュータ(「京」「富岳」など）
並列（「手分け」）アルゴリズムによるブレークスルー

参考：20世紀のトップ10アルゴリズム(F. Sullivan and J. Dongarra, Comp. Sci. Eng. 2, 22 (2000))
・ The Metropolis Algorithm (1946)
・ Simplex Method (1947)
・ Krylov Subspace Method (1950)
・ The Decompositional Approach
to Matrix Computations (1951)
・The Fortran Optimizing Compiler (1957)

・ QR Algorithm (1959)
・ Quicksort (1962)
・ Fast Fourier Transform (1965)
・ Integer Relation Detection (1977)
・ Fast Multipole Method (1987)

解説書籍「20世紀のトップ10アルゴリズム」，2022年→

背景 2：高度計測におけるデータ解析によるブレークスルー
データ解析発展で「⾒えないモノを⾒る」（計測限界突破）（従来計測⼿法では検出不可能だったモノを検出する）
例：クライオ電⼦顕微鏡（2017年ノーベル賞）
本研究の対象：先端計測・ ⼤型実験施設や⾼度実験技術が必要
・ データが「断⽚的」「⼀⾒⽭盾」し，「死蔵」されたデータや未解決問題が多数．
→ データ解析による計測⾰新への期待

研究開発サイクル
合成計測解析
❌
ボトルネックデータ駆動科学の立場から開発された，
新しい計測データ解析フレームワーク「2DMAT」
https://github.com/issp-center-dev/2DMAT/; Motoyama, et al., Comp. Phys. Commun. 280, 108465 (2022)
特徴1： 高信頼な計測データ解析を，高速に実行
・「見えないものを，見る」(測定限界突破）→ 不確かさ・ノイズに隠された「真の情報」を抽出．
・人的試行錯誤を伴わない→ 大域解析（解のinitial guessは，不要)・並列（ 「手分け」できる）アルゴリズム(PCからスパコンまで）
特徴2 ：汎用性：さまざまな計測手法に対応
・当座対象：2次元物質むけ計測

TRHEPDRHEED(1) 全反射高速陽電子回折(TRHEPD，トレプト)（本発表）(2) 反射高速電子回折(RHEED)(3) 表面エックス線回折(SXRD)(4) 低速電子回折(LEED)(5) (予定) オペランド偏光全反射蛍光(PTRF)
X線吸収微細構造(XAFS)(高草木・ 朝倉(北大触媒研))・他実験にも対応可．

2DMATSXRD
LEED
オペランド
PTRF-XAFSデータ駆動科学の立場から開発された，新しい測定データ解析フレームワーク「2DMAT」
https://github.com/issp-center-dev/2DMAT/; Motoyama, et al., Comp. Phys. Commun. 280, 108465 (2022)◎ ⼤域型逆問題解析
・ 内容：測定量𝐷 ≡ (𝐷! , . . 𝐷" )(観測データ)から⽬的状態量(欲しい量) 𝑋 ≡ (𝑋! , . . 𝑋# )を得る．
・ 解析原理：1) 物理学により，測定量𝐷を⽬的量𝑋の関数(順問題ルーチン)で与える
𝐷$%& ≡ 𝐷$%& 𝑋(1)→「順問題ルーチン」を取り替えることで様々な実験に対応
2) 計算値𝐷$%& 𝑋 と実測値𝐷の残差を⽬的関数𝐹 𝑋; 𝐷 とする．
例えば，𝐹(𝑋; 𝐷) ≡ ∑' 𝐷' − 𝐷$%&,' 𝑋→ 最適化・不確かさ解析を実⾏)

図 ⽬的関数𝐹 𝑋; 𝐷 の概念図

(2)・5種解析法が実装(1) 局所探索型最適化(Nelder-Mead法）(2) グリッド型探索(3) 並列化ベイズ最適化(PHYSBOライブラリ(*)を利⽤)(*)https://www.pasums.issp.u-tokyo.ac.jp/PHYSBO/(4) レプリカ交換(RE)モンテカルロ(MC)法(5) ポピュレーションアニーリング(PA)モンテカルロ(MC)法→「富岳」などの超並列計算機に最適

𝐹 𝑋; 𝐷

不確
かさ

不確
かさ
局所解

真解

ベイズ推定𝑋

対象実験:全反射高速陽電子回折(Total-reflection high-energy positron diffraction; TRHEPD, トレプト)
総説 ：Fukaya, et al., J. Phys. D 52, 013002 (2019)
解析例：[1] Hoshi, et al, Comp. Phys. Comm. 271, 108186 (2022)

・2次元的構造⾮破壊計測の⾰新⼿法（図1）．⽇本発．KEKで実⽤的測定可(全国共同利⽤)．
反射⾼速電⼦回折(Reflection high-energy positron diffraction; RHEED)の陽電⼦版
陽電⼦を⽤いることで，表⾯領域（厚さ< 1nm)のみをプローブできる．
・回折データ（各スポットの強度）から表⾯原⼦配列を決定する
2次元物質(表面) 構造(原子座標) 𝑋の決定例

実験装置模式図(https://www2.kek.jp/imss/spf/)(𝐷 ⟹ 𝑋)H
Ce2DMAT(解析ソフト)𝐷
触媒系：⽔素吸着したセリア(CeO2)表⾯,( Mochizuki et al. SLOPOS, 2023 )O

カーボンニュートラル向け(CO2再資源化)触媒開発
例：CO2+⽔+太陽光で灯油⽣成(厚み)≤1nm

逆解析

出⼝（の１つ）：(側⾯図)
灯油

Zoller et al., Joule 6,1606 (2022)

対象実験:全反射高速陽電子回折(Total-reflection high-energy positron diffraction; TRHEPD, トレプト)
総説 ：Fukaya, et al., J. Phys. D 52, 013002 (2019)
解析例：[1] Hoshi, et al, Comp. Phys. Comm. 271, 108186 (2022)

・2次元的構造⾮破壊計測の⾰新⼿法（図1）．⽇本発．KEKで実⽤的測定可(全国共同利⽤)．
⽇本語解説記事：
反射⾼速電⼦回折(Reflection high-energy
positron diffraction; RHEED)の陽電⼦版[a]⼀宮彪彦，「我が国における反射⾼速電⼦回折の発展と全反陽電⼦を⽤いることで，表⾯領域（厚さ<1nm)のみをプローブできる．
射陽電⼦回折への展開」, ⽇本物理学会誌，2015年
・回折データ（各スポットの強度）から表⾯原⼦配列を決定する[b] 深⾕有喜，「全反射⾼速陽電⼦回折（TRHEPD），陽電⼦学出⼝（の１つ）：
会誌，2019年
実験装置模式図

図(https://www2.kek.jp/imss/spf/)

2次元物質(表面) 構造(原子座標) 𝑋の決定例

TRHEPD装置模式図[b](𝐷 ⟹ 𝑋)

HCe2DMAT(解析ソフト)𝐷
触媒系：⽔素吸着したセリア(CeO2)表⾯,( Mochizuki et al. SLOPOS, 2023 )O(厚み)≤1nm

逆解析

カーボンニュートラル向け(CO2再資源化)触媒開発
例：CO2+⽔+太陽光で灯油⽣成

(側⾯図)
灯油

Zoller et al., Joule 6,1606 (2022)

例：[1] Hoshi, et al, Comp. Phys. Comm. 271, 108186 (2022)

TRHEPDにおけるデータ解析

解析原理：測定量（回折強度）における，実験データ𝐷と計算データ𝐷!"# 𝑋 の数値的⼀致→ 「R因⼦（残差）」 R 𝑋; 𝐷 ≡ 𝐷!"# 𝑋 − 𝐷 の2乗が，⽬的関数(𝐹(𝑋; 𝐷) ≡ R 𝑋; 𝐷 $ )図 TRHEPDにおけるデータ解析の概要図

00スポット強度 I

ロッキングカーブ(回折強度)

D実験
計算
𝑅 𝑋; 𝐷 ≈ 0 となる
Xを探索

実験

Dcal(X) 計算（試行座標X ）（多数）

視射⾓θ（度）順問題：量⼦散乱問題(偏微分⽅程式)求解(*)

z: 半無限系, xy⽅向のみ周期的
⼀宮, JJAP22, 176 (1983)(RHEED)(*)利⽤ルーチン：SIM-TRHEPD-RHEED (https://github.com/sim-trhepd-rheed/; Hanada, et al., Comp. Phys. Commun. 277, 108371 (2022)基礎：計測データ解析におけるベイズ推定・ 計測データ𝐷 ≡ (𝐷% , . . 𝐷& )が与えられた時に，⽬的状態量(欲しい量) 𝑋 ≡ (𝑋% , . . 𝑋' )の確率分布（ベイズ事後確率分布）を，ベイズの定理（下記）に基づいて得る

𝑃 𝑋|𝐷 ＝! "|$ !($)!(")

・例：病気Aには検査Bがあるが，誤判定する場合もある．

計測データ𝐷

今，検査Bを3回⾏ったところ，3回とも陽性と診断された．
実際に， 病気 Aにかかっている確率は，何パーセントか？
（計算に必要な情報：病気Aの患者の割合は，⽇本国⺠の10%．
検査Bは，病気Aに罹患している⼈のうち，30％を陰性と誤判定する．
検査Bは，病気Aに罹患していない⼈のうち，20％を陽性と誤判定する）

使用解析手法 : ポピュレーションアニーリングモンテカルロ法(PAMC法)
K. Hukushima and Y. Iba., AIP Conf. Proc. 690, 160 (2003).

PAMC法の概念図→ レプリカと呼ばれる多数サンプルデータを，乱数により⽣成・更新．
→ 100万以上の並列⾃由度．「富岳」に理想的 (右図)・パラメータ𝛽（𝛽 (% が観測データの不確かさ強度）に依存した事後確率分布𝑃(𝑋|𝐷; 𝛽)をヒストグラムとして表現(下図）．

パラメータ 𝛽!" 軸・並列化されたモンテカルロ法の1種で，⼤域解析を実現．

．
注：ベイズ事後確率分布𝑃(𝑋|𝐷)
をベイズの定理から計算．
𝑃 𝑋 𝐷; 𝛽𝑃 𝐷 𝑋; 𝛽 𝑃(𝑋)𝑃 𝑋 𝐷; 𝛽 =𝑃(𝐷; 𝛽)𝑃 𝐷 𝑋; 𝛽 ∝ exp −𝛽𝐹(𝑋, 𝐷)（ボルツマン型尤度関数）
𝐹(𝑋, 𝐷) ≡ 𝑅(𝑋, 𝐷)$ :（測定量についての実験と計算の2乗残差）
𝑃 𝑋 ：定義域で⼀様分布

サンプル⽣成軸ベイズ事後確率分布の模式図
𝑃 𝑋 𝐷; 𝛽𝛽値の⼤きさ⼤⼩

⽬的状態量 𝑋

⽬的状態量 𝑋

「富岳」での高速計算
・ 「富岳」で，PAMC法をコード調整中・ [暫定データ]並列ベンチマーク例(計算ノード数)=3,072- 82,944(全体の約1/2)で，強スケーリング型並列効率α=0.91 （注：理想値は α= 1）
詳細：高速化の論点(a) プログラム言語は，MPI並列化されたPython．MPI通信は，汎用Pythonオブジェクト(例： comm.send, comm.reduce ）ではなく，高速なbuffer-likeオブジェクト（例： comm.Send, comm.Reduce ）を利用．
(b) 順問題ルーチン(𝐷./0 𝑋 の計算)はFortranで書かれており， iso_c_bindingを通じて，外部ライブラリとして，Pythonコードにリンク．
(c) 計算途中でリサンプリング操作（レプリカ再配分）が実施される．
→ Walker’s tableアルゴリズムを利用（各MPIプロセスで，理想的ロードバランス）(d) 「富岳」特有の論点： Lightweight Layered I/O Accelerator （LLIO） の利用（データI/Oにおけるデータ通信を最適化）

TRHEPD計測におけるベイズ事後確率分布計算のデモ：Ge(001)-c4x2表面構造解析
Motoyama, et al., Comp. Phys. Commun. 280, 108465 (2022)

・最表⾯2原⼦のz座標 𝑧% , 𝑧$ に対するベイズ事後確率分布𝑃 𝑧% , 𝑧$ 𝐷
→ ⼤域解析(initial guess不要)．不確かさパラメータ値𝛽により，真解・局所解を検出・区別できる．
(a) 原子構造

9.0z

局所解

10!"

8.0

7.0

7.0

10!#

10!#

6.0

6.0

10!$

5.0

5.0

10!$

4.0

4.0

真解

3.0

2.0

5×10!"(c)不確かさ小9.010!"8.0𝑧) (Å)

側⾯図

5×10!"(b)不確かさ大2.0

3.0

4.0

5.0

6.0𝑧! (Å)

7.0

10!%

8.0

2.0 2.0

9.0!&

10

真解

3.0

<3.0

4.0

5.0

6.0𝑧! (Å)

7.0

10!%8.0

9.010!& <注：実測定データの代わりに，⼈⼯データを利⽤結果：Ge(001)-c4x2表面構造解析 2変数探索
ベイズ事後確率分布𝑃 𝑧! , 𝑧) 𝐷 の不確かさパラメータ依存性（不確かさ：⼤→⼩）

𝑍) (単位：0.1 nm)

等⾼線図の値：R因⼦(残差)R 𝑧% , 𝑧$ ×100

𝑍! (単位：0.1 nm)

6Ge(001)-c4x2構造への応⽤：
表⾯1~6層のz座標の⼤域解析

543・真解が既知である構造を使った，⼿法テスト
・⼀波条件(z座標のみの決定）・側⾯図（右図）で，0 < 𝑧 < 9 Å での⼤域解析→ 事後確率分布𝑃 𝑋 𝐷; 𝛽 をヒストグラムとして⽣成．( 𝑋 ≡(𝑧4, 𝑧5 … 𝑧6) )2
1𝑧'𝑧(𝑧) , 𝑧*𝑧+𝑧,

0-1 単位: ÅGe(001) -c4x2（側⾯図）(これまでの）まとめ：
・ 「2DMAT」: データ駆動科学の⽴場から開発された，データ解析フレームワーク「2DMAT」(https://github.com/issp-center-dev/2DMAT/)
・ 超並列モンテカルロ法(PAMC法)によるベイズ推定（ベイズ事後確率分布計算）→「富岳」などの超並列計算機に理想的．
・全反射⾼速陽電⼦回折実験(TRHEPD,トレプト)で実践．
→回折データ𝐷から，表⾯付近（厚さ：サブナノメートル）の原⼦座標 𝑋 を推定．
→⼤域解析(解へのinitial guess不要)で，真解が得られた．
・展望：(i) 「見えないものを，見る」(測定限界突破）
例：表⾯⽔素の検出（→次世代触媒・電池）(望⽉，招待講演，SLOPOS-16, France, 2023年7⽉16⽇)

⽔素吸着したセリア(CeO2)表⾯,HCe(ii) 種々実験（ TRHEPD, RHEED, SXRD, LEED,
オペランドPTRF-XAFS)，DFT計算との連携→マルチモーダル（複数⼿法）計測における．統合的データ解析（関連：KEK物構研での試料ケース共通化）

O⼤局的展望 : デジタルツインなど
・計測データ解析は汎⽤→物質科学以外への展望（例：核融合科学）
・⼤局的展望→⾼エネ研・核融合研などでの， 次世代実験装置開発にむけたデジタルツイン開発（要素パラメータが向上したときに，装置性能はどの程度向上するか?）
外挿型予測の例：
スパコンの計算性能予測（ベイズ推定）

現⾏実験装置

学習データ
物理ベースモデル

デジタルツイン（外挿的性能予測など）
設計

次世代実験装置[予測：メジアン，信頼性区間の上下端をプロット]
学習データ

テストデータ上端

メジアン
下端

Kohashi et al., JSIAM Letters 14, 13 (2022)