第二回とめ研究所若手研究者懸賞論文

人工ニューラルネットワークに基づく
銀カルコゲナイドの原子間ポテンシャルの構築と熱伝導度計算への応用

熊本大学大学院 自然科学教育部
理学専攻 物理科学コース
学籍番号 202-D1207

竹下 雄輔

指導教員: 下條 冬樹 教授

0

1

目次
目次..................................................................................................................................... 2
1章

序論 .......................................................................................................................... 4

1.1 節

研究背景 ...................................................................................................................... 4

1.1.1 項

熱デバイスの重要性 ............................................................................................ 4

1.1.2 項銀カルコゲナイドの可能性 ................................................................................. 5

1.1.3 項

格子熱伝導度に注目する理由.............................................................................. 6

1.1.4 項

格子熱伝導度の計算方法とその問題点 ............................................................... 7

1.1.5 項

Artificial Neural Network (ANN) ポテンシャルの利用 ...................................... 8

1.1.6 項ANN ポテンシャルに関する課題 ........................................................................ 8

1.2 節

研究目的 .................................................................................................................... 10

1.3 節

本研究の流れ ............................................................................................................. 10

2章

計算手法 ................................................................................................................. 11

2.1 節

分子動力学法 ............................................................................................................. 11

2.1.1 項

力の計算 ............................................................................................................ 11

2.1.2 項時間発展 ............................................................................................................ 11

2.2 節

第一原理計算 ............................................................................................................. 122.2.1 項Born-Oppenheimer 近似 ................................................................................... 13

2.2.2 項

密度汎関数理論 ................................................................................................. 13

2.2.3 項Kohn-Sham 方程式 ............................................................................................ 14

2.2.4 項Self-Consistent Field 法 ..................................................................................... 15

2.2.5 項Hellmann-Feynman の定理 ............................................................................... 15

2.2.6 項交換相関エネルギー .......................................................................................... 16

2.2.7 項内殻電子の扱い ................................................................................................. 16

2.3 節

Artificial Neural Network (ANN) ポテンシャル ....................................................... 16

2.3.1 項ANN ポテンシャルとは .................................................................................... 16

2.3.2 項

学習.................................................................................................................... 18

2.3.3 項

ANN-MD........................................................................................................... 20

2.3.4 項Ewald 法との併用 .............................................................................................. 21

2.4 節

熱伝導度計算 ............................................................................................................. 22

2.4.1 項

平衡分子動力学法による計算方法 .................................................................... 222

2.4.2 項

非平衡分子動力学法による計算方法 ................................................................. 23

2.4.3 項熱流束の表式 ..................................................................................................... 243章

計算条件・結果・考察 ............................................................................................. 26

3.1 節

教師データ (DFT+U) の計算条件 ............................................................................ 26

3.2 節

ANN ポテンシャルに関する各種計算条件 ............................................................... 27

3.3 節

Virial 学習の比重調整過程 ........................................................................................ 28

3.4 節

Virial 学習の効果....................................................................................................... 33

3.5 節

Ewald 法併用の効果 .................................................................................................. 36

3.6 節

熱流の正誤による熱伝導度の変化 ............................................................................ 38

3.7 節

非平衡分子動力学法による計算結果 ........................................................................ 40

3.8 節

対称関数の変更による影響 ....................................................................................... 44

3.9 節

スーパーセルの形状が熱伝導度に与える影響 .......................................................... 47

4章

結論 ........................................................................................................................ 47

付録 A

計算手法、条件の補足 ...................................................................................... 50

A.1

対称関数 ........................................................................................................................ 50

A.2

対称関数のパラメータ .................................................................................................. 51

引用文献 ............................................................................................................................ 48

3

1章 序論本研究では熱デバイスの材料として注目されている銀カルコゲナイドのうち、Ag2Se に注目し、その格子熱伝導度をシミュレーションで調べた。手法としては、第一原理分子動力学 (First-Principles Molecular Dynamics, FPMD) 法のデータを学習した人工ニューラルネ
ットワーク (Artificial Neural Network, ANN) に基づく原子間ポテンシャル (ANN ポテンシ
ャル)を用いた MD (ANN-MD) シミュレーションの平衡状態のデータと Green-Kubo (GK)
公式を組み合わせて格子熱伝導度を計算する手法を用いた。特に本研究では、格子熱伝導度計算を十分な精度で行うために、ANN ポテンシャル構築手法や格子熱伝導度の計算方法の改良に重点をおいた。改良の末、十分な精度で格子熱伝導度を計算出来るようになり、
Ag2Se の格子熱伝導度は 300 K において 0.306±0.013 Wm−1 K −1 だと見積もられた。この値
は実験値(0.4~0.5 Wm−1 K −1 ) [1, 2]に近い。以下ではまず、本研究の背景 (1.1 節) と目的(1.2 節)、流れ (1.3 節) について紹介する。

1.1節

研究背景

持続可能な社会の実現の一助として未利用エネルギーの活用が重要であり、それを進めるためには実用性の高い熱デバイスの開発が必要である。このデバイスの実用性を高めてくれる材料候補として銀カルコゲナイドが注目されている。そしてこの材料を熱電デバイ
スとしてうまく制御したり、適した組成を探ったりするために熱伝導度についての理解が必要となる。我々は ANN-MD 法と GK 公式を組み合わせる手法によって格子熱伝導度を計算する手法に注目している。十分な精度で ANN-MD を行うためには、ANN ポテンシャ
ル構築のノウハウを蓄積することも重要なので、本研究ではこの点にも重点を置いている。
以下ではまず本研究の背景として熱デバイスの開発の重要性 (1.1.1 項)、銀カルコゲナイドの可能性 (1.1.2 項)、格子熱伝導度を計算する理由 (1.1.3 項)、格子熱電度の計算方法とその問題点 (1.1.4 項)、ANN-MD 法によるシミュレーションデータの収集 (1.1.5 項)、ANN
ポテンシャルに関する課題 (1.1.6 項)について説明する。

1.1.1項 熱デバイスの重要性
持続可能な社会の実現に向けて、未利用エネルギーの活用が重要である。未利用エネル
ギーとは工場排熱やエネルギー転換時の損失などのように有効活用出来ていないエネルギーのことである。例えば石油、天然ガス、太陽光などのような 1 次エネルギーの国内供給
を電気エネルギーなどに転換する際に約 1/3 のエネルギーが熱などとして失われるというデータが有る (図 1.1) [3, 4]。また、日本国内では年間 1 兆 kWh 程度のエネルギーが未利用熱として捨てられているという報告もある [5]。このように未利用エネルギーの量は莫4

大なので、それを一部でも有効活用できればエネルギー的に無駄の少ない社会の実現に役立つと考えられる。

図 1.1 一次エネルギーの国内供給に対する最終エネルギー消費の割合 [3, 4]
また、未利用熱からの発電は Internet of Things (IoT)やウェアラブルデバイスの普及を後押しする技術としても期待されている。例えば IoT の中には大量のセンサーの設置が必要なものもあり、それらの電源や配線をいかに確保するかが問題視されている。ここで未利用熱からの発電によって十分な電力を供給できれば、これらの問題は解決する。また、ウェアラブルデバイスについても、体温と気温との温度差で発電できるようになれば、充電の必要がなくなって活用の幅が広がると考えられる。
以上のような未利用のエネルギーの活用を可能にするデバイスが熱デバイスである。熱デバイスの種類としては熱電変換素子や熱整流素子、熱スイッチなどがある。これらがあれば、熱が作り出した温度差から電気エネルギーを抽出したり、熱源付近から熱が必要なところへ熱を輸送したり出来るようになる。未利用エネルギーの多くが熱であることからも、これらのデバイスの実用化によって可能になるエネルギー再利用の規模は大きいと考えられるので、熱デバイスの開発は重要である。

1.1.2項 銀カルコゲナイドの可能性
数ある材料の中でも銀カルコゲナイドは熱デバイス材料として期待されていて [6, 2, 7, 8,9]、その理由は以下のような銀カルコゲナイドの特徴にある。
⚫

熱伝導度が非常に低い。

⚫

一般的な結晶と異なり、温度が上昇しても熱伝導度が上昇しない [10]。

5⚫

非伝導相と超イオン伝導相 (銀イオンがカルコゲン副格子間を拡散する相 [11]) の間の相転移温度が室温付近にあり、この相転移温度と熱伝導度の温度依存性を組成によって変化させられる(図 1.2) [2]。

⚫

組成次第で毒性や延性も変化する

図 1.2 銀カルコゲナイド (Ag 2 S1−𝑥 Se𝑥 ) の熱伝導度の温度依存性。350~450 K の間に相転移温度がある。低温相は非伝導相で高温相は超イオン伝導相である [2]。
これらの特徴から熱デバイスに適した材料だと考えられている。例えば、熱電材料の性能
は無次元性能指数𝑧𝑇 [12]𝑆 2 𝜎𝑇
1.1.1
𝜅
によるので、熱伝導度𝜅が非常に低いという点は好都合である。ここで𝑆はゼーベック係数、𝑧𝑇 =𝜎は電気伝導度、𝑇は温度、𝜅は熱伝導度である。また、組成次第で室温付近にある相転移
温度や熱伝導度の温度依存性を変化させられることは、熱整流素子にとって好都合である。
以上の理由から銀カルコゲナイドには注目する価値が十分有る。

1.1.3項 格子熱伝導度に注目する理由
熱デバイス制御においては熱伝導度に対する深い理解が重要である。なぜなら1.1.2 項で説明した通りに熱伝導度は熱電変換素子や熱整流素子の性能に直接影響するからである。
熱伝導度𝜅は格子熱伝導度𝜅l と電子熱伝導度𝜅e に分けられる。
𝜅 = 𝜅l + 𝜅e

1.1.2

我々は熱伝導度の中でも特に格子熱伝導度𝜅l に注目している。その主な理由は次の２つで
6

ある。
⚫

実験で直接計測出来るものは𝜅であり、𝜅l と𝜅e を区別することは難しいため、理論的に計算することが重要であること。

⚫

電子熱伝導度は電気伝導度などと独立に制御するのが難しいが、格子熱伝導度は比較的独立に制御しやすいため、熱電材料の性能向上に役立つ知見を得られる可能性が有ること。

２点目の理由に関して補足する。多くの金属に対して当てはまる経験則の Wiedemann–Franz 則 [13, 14]
𝜅e = 𝐿𝜎𝑇

1.1.3

から電子熱伝導度と電気伝導度は比例関係であるため、電子熱伝導度の制御によって無次元性能指数𝑧𝑇 (1.1.1 式)を大きくするのは難しい。一方で、格子熱伝導度であれば比較的独立して制御出来るため無次元性能指数𝑧𝑇の向上に役立つと考えられる。
以上のような理由から我々は特に格子熱伝導度𝜅l に注目して研究を行った。以下では電子
熱伝導度𝜅e は登場しないので簡単のために、格子熱伝導度を𝜅として表記している。

1.1.4項 格子熱伝導度の計算方法とその問題点
数ある格子熱伝導度計の算方法の中で、銀カルコゲナイドのような複雑な振る舞いをす
る物質に適しているのは、平衡状態の MD (Equilibrium Molecular Dynamics, EMD) 計算と Green-Kubo (GK) 公式を用いる手法 [15, 16, 17]である。銀カルコゲナイドの超イオン伝導相は Ag イオンの振る舞いが明らかに複雑であり、更に低温の非イオン伝導相においても Ag 原子が異常な振る舞いを見せるということが知られている [7]。そのため、この系の熱伝 導度 を適 切に 見積 もれ る手 法は 限ら れて い る。例 えば 、格 子動 力学 (Lattice
Dynamics, LD) 法に基づいた計算手法 [18]は結晶に対して有効であるものの、超イオン伝
導相のように液体に近い振る舞いを含む系には使えない。一方 EMD 計算と GK 公式を用いる手法であれば、原子の振る舞いが結晶的であるかなどの制限なしに適用できるので、
銀カルコゲナイドに対しても使える。したがって本研究ではこの手法を主に用いた。
しかし、これらの手法ではシミュレーション時間を長くとる必要があり [19]、計算コス
トの高い第一原理分子動力学 (First-Principles Molecular Dynamics, FPMD) 法との相性は悪い。例えば、経験的ポテンシャルによって銀カルコゲナイドの熱伝導度を計算した先行研究によると、熱伝導度を十分収束させるのに 384 原子程度の系で 2.4 ns 程度のシミュレ
ーション時間が必要であり、この量の計算を私の使える計算資源で FPMD 法に基づいて行う場合にはおよそ１年程度計算し続ける必要がある。これを行うのは現実的ではない。

7

1.1.5項 Artificial Neural Network (ANN) ポテンシャルの利用
前項で触れた通りに EMD 法による長時間のシミュレーションデータを FPMD 法で収集するのは難しい。そこで人工ニューラルネットワーク (Artificial Neural Network, ANN) ポ
テンシャル [20, 21]を用いた MD (ANN-MD) 計算が有効だと考えられる。ANN ポテンシ
ャルについては2.3 節で詳細に記している。簡単に説明すると機械学習を応用した原子間ポ
テンシャルであり、計算コストは経験的ポテンシャル程度に低く、精度は教師データに依存する。つまり、FPMD 計算の結果を学習させれば、FPMD 計算と同程度の精度の結果を
経験的ポテンシャル程度の計算コストで与える原子間ポテンシャルになる。この ANN ポ
テンシャルを用いた場合、例えば、384 原子の銀カルコゲナイドに対して 2.4 ns 分の計算
を 12 時間程度でできるようになる。このように ANN ポテンシャルによって、EMD 計算と GK 公式を合わせた熱伝導度計算は十分可能なものになる。
また、ANN ポテンシャルによる恩恵としては計算コストの低減だけでなく、Ag2S1-xSexのような経験的ポテンシャルが報告されていない３元素以上から成る系に対して、熱伝導度などの計算が可能になることも挙げられる。Ag2Se などの２元系に対する経験的ポテン
シャルは報告されているので、それらで熱伝導度などの長時間のシミュレーションデータ
が必要な物理量を古典 MD 法で見積もることは可能である。しかし、３元素以上から成る
系には経験的ポテンシャルが報告されていない場合が多いので、古典 MD 法で見積もることが難しいという問題がある。この問題は ANN ポテンシャルで解決出来る。なぜなら、対象とする系の FPMD データさえ有れば ANN ポテンシャルを構築して、長時間のシミュ
レーションデータの収集が可能になるからだ。このように、ANN ポテンシャルには計算の高速化だけでなく、経験的ポテンシャルが報告されていない系に対しても熱伝導度などの長時間のシミュレーションデータが必要な物理量の計算を可能にするという利点もある。

1.1.6項 ANN ポテンシャルに関する課題
ANN ポテンシャルの精度は様々なハイパーパラメータに依存するので、高精度に熱伝導度を見積もるにはそれらを十分に調整しなければならない [22, 23]。ハイパーパラメータとは我々人が調節しなければならないパラメータであり、広義には以下のようなものが含まれる。学習対象とする物理量の選択、損失関数の係数のバランス、Ewald 法 (クーロン相互作用を計算する手法)との併用の有無、原子配置からの特徴抽出に用いる対称関数の種類
やパラメータ、隠れ層の層数やノード数などである。以下３段落では、学習対象とする物
理量、Ewald 法、対称関数を例にとってこれらのハイパーパラメータの調整の重要性について述べる。続けて、ANN ポテンシャルの初期重み依存性の考慮の必要性や ANN ポテン
シャルに対する厳密な熱流束の表式を用いる必要性についても述べる。
学習対象とする物理量に関しては、エネルギーだけを学習対象とする従来手法が良く知8られていた。しかし、Ag2Se のような振る舞いが複雑な系においては不十分であり、エネ
ルギーに加えて力も学習することで部分動径分布関数などの基本的な物理量を再現出来るようになることが知られている [24]。ところが、エネルギーと力を学習するだけでは圧力の精度が十分でないため、さらに圧力 (実際には Virial テンソル) の学習による影響を調べてみる必要がある (3.4 節)。また、この影響を調べるためには損失関数の係数のバランスを整えて圧力を含めた学習が十分な精度で行えるように調整する必要も有る (3.3 節)。
Ag2Se はイオン性の高い物質なのでクーロン相互作用を明示的に考慮するために Ewald 法が必要だと考えられる。しかし、Ewald 法は固定電荷でしかクーロン相互作用を計算出来ないため ANN ポテンシャルを教師データと異なる温度などの条件下で用いる場合に注意が必要であったり、そもそも電荷の任意性があったりするなどの問題点が有る。そのため、Ewald 法の併用の有無で熱伝導度などの結果が変わらないのであれば併用しない方が良い。
したがって、Ewald 法の併用の有無による熱伝導度などへの影響を調べる必要がある (3.5節)。
原子配置からの特徴抽出に用いる対称関数には Behler−Parrinello 型 [25]や Chebyshev
型 [26]などがあり、我々は以 前から Behler−Parrinello 型を用い ていた。しか し、
Chebyshev 型の方がパラメータ調整は容易であるため、今後他の物質に ANN ポテンシャ
ルを応用していく際に有利だとだと考えられる。したがって、この対称関数の種類の変更
による影響を調べる必要もある (3.8 節)。
このように十分な精度の ANN ポテンシャルを構築するために、様々なパラメータを調整しなければならない。パラメータ調整の詳細は3 章に示してある。本研究でパラメータ調整に関する知見がある程度は集まったが、まだ良くわからない部分もある。
以上のようなパラメータを固定したとしても、構築される ANN ポテンシャルはその初期重みに依存し、ANN-MD の結果も変わってくる。従って、上記のようなパラメータの調整や ANN-MD による物理量の見積もりを行う際には、初期重みの異なる複数の ANN ポ
テンシャルを構築してそれらによる結果を吟味する必要がある (ANN ポテンシャルの初期重み依存性の考慮)。またこれに伴って、対象とする物理量や系に対していくつの ANN ポ
テンシャルを用意すれば良いのかという問題も生じ、計算コストも増えてくる。本研究の場合 (Ag2Se の格子熱伝導度) は 5 つ用意すれば良いと経験的に分かったため、3 章の結果の殆どはそのようにして求めたものになっている。
また、熱伝導度計算に関しては、ANN ポテンシャルに対応した熱流束の表式を用いる事も重要である。従来一般的だった熱流束には圧力テンソルから自然に導かれる原子毎のVirial テンソルが使われており、これは ANN ポテンシャルの多体効果を考慮していないものであった。そのため、この熱流束から導かれる熱伝導度は実験値とかけ離れたものであり、エラーバーも大きかった。これに関して、ANN ポテンシャルの多体効果を考慮した熱流束とそうでない熱流束による熱伝導度の違いを確認する必要もある (3.6 節)。

9

1.2節

研究目的

以下ではまず背景を簡単にまとめた後に本研究の目的を述べる。
未利用熱の活用のために熱デバイスの開発が重要であり (1.1.1 項)、その材料として銀カルコゲナイドが注目されている (1.1.2 項)。銀カルコゲナイドの物性値の中でも熱伝導度、
特に格子熱伝導度についての深い理解は熱デバイスの制御に重要だと考えられるので、格子熱伝導度の温度や組成依存性を見積もることが重要である (1.1.3 項)。銀カルコゲナイドの非結晶的な振る舞いを考慮すると、格子熱伝導度を計算する手法としては EMD 法による長時間シミュレーションと GK 公式を用いる手法が有効だと考えられる (1.1.4 項)。長時間シミュレーションのデータを FPMD 法で計算するのは計算コスト的に現実的では無いため、ANN-MD 法によって FPMD 法と同程度の精度の結果を低い計算コストで得る手法を我々は用いることにした (1.1.5 項)。ただし、十分な精度の ANN ポテンシャルを構築する
ためには様々なハイパーパラメータの最適化が必要であり、これについての調査も必要である (1.1.6 項)。
以上から、我々の研究の主要な目的は、銀カルコゲナイドの格子熱伝導度の温度・組成依存性を ANN-MD 法を通して明らかにすることと、そのための手法を構築することである。これを達成するために、まずは１つの組成に対して高精度に格子熱伝導度を見積もれるようになる必要がある。そこで本研究では Ag2Se を研究対象として、以下のことを目的とした。
⚫

ANN-MD 法を用いて格子熱伝導度を高精度に見積もるための ANN ポテンシャル構
築手法や格子熱伝導度の計算手法の確立

1.3節

本研究の流れ

本研究の流れは次のようになっている。まず ANN ポテンシャルの精度を向上させるための調査として、エネルギーと力に加えて Virial テンソルも学習させることによる影響 (3.3
節と3.4 節) 、Ewald 法を併用して明示的にクーロン相互を計算することによる影響 (3.5 節)を調べた。また、熱伝導度計算の精度を上げるために、熱流束の表式による影響 (3.6 節)も調べた。さらに補助的な調査として、HNEMD 法によって計算する事による影響 (3.7節) や対称関数の変更による影響 (3.8 節)を調べた。

10

2章 計算手法
本研究で用いた計算手法を大別すると、MD 法 (2.1 節)、第一原理計算 (2.2 節)、ANN ポ
テンシャル (2.3 節)、熱伝導度計算 (2.4 節)の４つになる。以下ではこれらの詳細について
紹介する。

2.1節

分子動力学法

分子動力学 (MD) 法は古典 MD 法と FPMD 法に大別される。両者とも原子や分子に働く力に基づいて運動方程式を解き、それらの時間発展を追う手法である点は共通である。主に異なるのは原子や分子に働く力の求め方である。以下では MD 計算での力の計算 (2.1.1
項) と時間発展の方法 (2.1.2 項) について紹介する [27]。

2.1.1項 力の計算
古典 MD 法での原子間相互作用は、原子の位置を引数とするポテンシャル関数によって与えられる。このポテンシャル関数は物質ごとに物性データを再現するように経験的に決められるものであり、経験的ポテンシャルと呼ばれる。これには FPMD 法に比べて低い計算コストで原子間相互作用を計算出来るという利点がある。しかし、適切な関数形の決定は素人には難しく、電子状態に基づいた量子力学的性質を考慮していないので精度において FPMD 法には及ばないという欠点がある。
一方 FPMD 法での原子間相互作用は、原子の位置に基づいて電子状態を量子力学的に考慮した上で導かれるものである。そのため、古典 MD 法よりも現実の系に近い結果を得られるという利点がある。しかし、電子状態が収束するまで波動方程式を繰り返し解くという方法 (つじつまの合う場の方法, Self-Consistent Field 法, SCF 法) をとるので計算コストが非常に高いという欠点がある。そのため、熱伝導度のように収束まで時間のかかる物理量の計算には不向きである。ちなみに FPMD 法での原子間相互作用の計算は、第一原理(First-Principles, FP) 計算と呼ばれている。つまり FPMD 法とは、FP 計算に基づいた MD法だと言える。

2.1.2項 時間発展
粒子の位置と速度を時間発展させる手法でよく用いられるのは RESPA(reversible reference system propagator algorithm)法である。この手法は異なるΔ𝑡を用いても可逆性が保たれるように設計されている。ここではその時間発展の計算式を導く。
一般化座標𝒒 = (𝑞1 , 𝑞2 , ⋯ , 𝑞𝑛 )とそれに共役な一般化運動量𝒑 = (𝑝1 , 𝑝2 , ⋯ , 𝑝𝑛 )に依存する物

11

理量𝐴(𝒒, 𝒑)の時間微分は以下のように書ける。
𝑛𝑑𝐴
𝜕𝐴 𝑑𝑞𝑖 𝜕𝐴 𝑑𝑝𝑖
)
= ∑(
+
𝑑𝑡
𝜕𝑞𝑖 𝑑𝑡 𝜕𝑝𝑖 𝑑𝑡

2.1. 1𝑖=1

ここでリウビル演算子
𝑛𝑖ℒ = ∑ (
𝑖=1

𝑑𝑞𝑖 𝜕
𝑑𝑝𝑖 𝜕
)
+
𝑑𝑡 𝜕𝑞𝑖 𝑑𝑡 𝜕𝑝𝑖を導入すると、𝐴(𝒒, 𝒑)の時間微分は𝑑𝐴(𝑡)= 𝑖ℒ𝐴(𝑡)𝑑𝑡となり、形式的に解くと次のようになる。
𝐴(𝑡) = exp (𝑖ℒ𝑡𝐴(0))

2.1. 2

2.1.3

2.1.4

物理量𝐴(𝒒, 𝒑)が𝒒, 𝒑である時は𝑞𝑖 (𝑡) = exp(𝑖ℒ𝑡) 𝑞𝑖 (0)2.1.5

𝑝𝑖 (𝑡) = exp(𝑖ℒ𝑡) 𝑝𝑖 (0)2.1.6となり、時間の起点を変えると次のようになる。
𝑞𝑖 (𝑡 + Δ𝑡 ) = exp(𝑖ℒΔ𝑡) 𝑞𝑖 (𝑡)2.1.7

𝑝𝑖 (𝑡 + Δ𝑡) = exp(𝑖ℒΔ𝑡) 𝑝𝑖 (𝑡)2.1.8

𝑖番目の自由度に対するリウビル演算子は
𝑑𝑞𝑖 𝜕
𝑑𝑝𝑖 𝜕
𝜕
𝜕
𝑖ℒ 𝑖 =
+
= 𝑞𝑖̇
+ 𝐹𝑖
2.1.9
𝑑𝑡 𝜕𝑞𝑖 𝑑𝑡 𝜕𝑝𝑖𝜕𝑞𝑖𝜕𝑝𝑖と書けるので、exp(𝑖ℒ 𝑖 Δ𝑡)に対称化されたトロッタ分解を用いると、次式のようになる
Δ𝑡
𝜕
𝜕
Δ𝑡
𝜕
) exp (Δ𝑡𝑞𝑖̇
) exp ( 𝐹𝑖
) + 𝑂(Δ𝑡 3 )exp(𝑖ℒ 𝑖 Δ𝑡) = exp ( 𝐹𝑖
2.1.10
2 𝜕𝑝𝑖
𝜕𝑞𝑖
2 𝜕𝑝𝑖
Δ𝑡

これは異なる時間幅 2 とΔ𝑡に関して対称な形であり、時間発展に対する可逆性が保たれている。exp(𝑖ℒΔ𝑡)のうち𝑞𝑖 や𝑝𝑖 に作用するのはexp(𝑖ℒ 𝑖 Δ𝑡)部分だけなので、𝑞𝑖 (𝑡 + Δ𝑡 )と𝑝𝑖 (𝑡 + Δ𝑡)の形式解から次の式が導ける。
𝑞𝑖 (𝑡 + Δ𝑡 ) = 𝑞𝑖 (𝑡) + Δ𝑡𝑞̇ 𝑖̇ (𝑡) +Δ𝑡 2
𝐹 (𝑡) + 𝑂(Δ𝑡 3 )
2𝑚𝑖 𝑖

Δ𝑡 2
{𝐹𝑖 (𝑡) + 𝐹𝑖 (𝑡 + Δ𝑡)} + 𝑂(Δ𝑡 3 )2
これらの式に基づいて系を時間発展させていく。
𝑝𝑖 (𝑡 + Δ𝑡) = 𝑝𝑖 (𝑡) +2.2節

2.1.112.1.12

第一原理計算FPMD (First-Principles Molecular Dynamics, FPMD) 法とは、第一原理 (First-Principles,FP) 計算に基づいて計算された力を用いる分子動力学法のことであり、FP 計算による力は、

12量子力学的に電子状態を計算し、原子核や電子間の相互作用を考慮した上で与えられる。
この方法の利点としては、電子、原子核の質量、スピン以外のパラメータを基本的に必要としない非経験的手法である点と、量子力学に基づいているため精度の高い点が挙げられる。しかし、毎ステップ電子についての波動方程式を解く必要があるため古典 MD 法に比べて莫大な計算時間がかかってしまうという欠点がある。
FP 計算では「静止した原子核の分布に対して全電子が動いて基底状態をとる」という近似 (2.2.1 項) のもとで全電子から原子核への作用を計算する。その作用は全電子のエネル
ギーが分かれば求められる。全電子のエネルギーは基底状態において電子密度の汎関数として与えられることが知られており (2.2.2 項)、これを利用して基底状態のエネルギーを求める方法がコーン・シャム（Kohn-Sham, KS）方程式 [28]を SCF (Self-Consistent Field)法に従って解くという方法である (2.2.3 項、2.2.4 項)。この方法で求められた基底状態のエネルギーを用いれば Hellmann-Feynman の定理に従って当初の目的である全電子から原子核への作用を計算出来る (2.2.5 項)。以下ではこれらの詳細を紹介する。

2.2.1項 Born-Oppenheimer 近似
系は外部からの励起作用が無ければ安定な最低エネルギー状態 (基底状態) に留まると考えられる。そして、電子は原子核に比べて非常に軽いので、原子の動きに瞬時に反応して基底状態の分布を取ると考えられる。この時の電子の運動は静止した原子核によるクーロ
ン引力場と電子間のクーロン斥力場のもとでの運動として扱える。このような扱いをボルン・オッペンハイマー(Born-Oppenheimer, BO)近似という。電子から原子核への作用は
BO 近似のもとで得られる基底状態の全電子のエネルギーの原子核座標に対する勾配から計算できるので、以降ではこのエネルギーの計算に焦点を置く。

2.2.2項 密度汎関数理論
基底状態における全電子のエネルギー𝐸 𝐷𝐹𝑇 は電子密度𝑛(𝒓)の汎関数として次式で与えられる。
DFT [𝑛]𝐸 𝐷𝐹𝑇 [𝑛] = 𝑇0 [𝑛] + 𝐽[𝑛] + ∫ 𝑉 ext (𝒓) 𝑛(𝒓)𝑑𝒓 + 𝐸𝐱𝐜

2.2.1

これが密度汎関数理論 (Density Functioanl Theory, DFT) の主張である。ここで電子密度𝑛(𝒓)は次式
𝑛 (𝒓) = ∑ 𝜑𝑘∗ (𝒓)𝜑𝑘 (𝒓)

2.2.2

𝑘

で与えられ（𝜑𝑘 (𝒓)は𝑘番目の１電子軌道関数である。）、𝑇0 [𝑛]は分布𝑛(𝒓)に従った電子の自由 な 運 動 エ ネ ル ギ ー 、 𝐽[𝑛] は 電 子 間 の ク ー ロ ン 相 互 作 用 に よ る エ ネ ル ギ ー 、DFT [𝑛]は交換∫ 𝑉 ext (𝒓) 𝑛(𝒓)𝑑𝒓は原子核が作る外場からの作用𝑉 ext (𝒓)によるエネルギー、𝐸xc

13相関エネルギーの期待値である。特筆するべきは、電子の運動に付随する相互作用のように関数形が不明確な作用を全て交換相関エネルギーに盛り込むことで式の扱いを容易にしている点である。交換相関エネルギーの具体的な電子密度への依存性を決定する近似手法については後の Error! Reference source not found.で簡単に触れる。
2.2.1 式の𝐸 𝐷𝐹𝑇 [𝑛]は基底状態においてのみ物理的意味が有る。基底状態の𝐸 𝐷𝐹𝑇 [𝑛]を求める方法は、KS 方程式を SCF 法に従って解くという方法である。次項では KS 方程式について紹介する。

2.2.3項 Kohn-Sham 方程式
コーン・シャム（Kohn-Sham, KS）方程式 [28]は電子軌道関数の規格化直交条件のもとで全エネルギーの期待値𝐸 𝐷𝐹𝑇 [𝑛]に最小値を取らせるときに現れる条件式である。具体的にはラグランジュの未定乗数法に従って書き下した2.2.3 式が 1 電子軌道関数𝜑𝑘,𝜎 (𝒓)の変分に対して 0 になるようにすれば得られる。
′

𝜎𝜎∗𝛿 {𝐸 𝐷𝐹𝑇 [𝑛] − ∑ ∑ [𝜀𝑘𝑘′ ∫ 𝜑𝑘,𝜎 (𝒓) 𝜑𝑘 ′ ,𝜎 ′ (𝒓) 𝑑𝑟 − 𝛿𝑘𝑘 ′ 𝛿𝜎𝜎 ′ ]}2.2.3𝑘,𝜎 𝑘 ′,𝜎 ′′𝜎𝜎
ここで𝜑𝑘,𝜎 (𝒓)は𝑘番目の空間軌道上のスピン𝜎 の 1 電子軌道関数、𝜀𝑘𝑘′ はラグランジュの未定乗数で、𝛿𝑘𝑘 ′ はクロネッカーのデルタ記号である。これを解いた結果を対角化したもの
が2.2.4 式の KS 方程式である。
̂ 𝐷𝐹𝑇 𝜑𝑘,𝜎 (𝒓) = 𝜀𝑘𝜎 𝜑𝑘,𝜎 (𝒓)
𝐻

2.2.4

̂ 𝐷𝐹𝑇 は全電子系のハミルトニアンであり、2.2.5 式のように与えられ𝜀𝑘𝜎 は固有値である。𝐻
る。
1
𝜎 (𝒓)
̂ 𝐷𝐹𝑇 = − ∆ + 𝑉 𝑒𝑥𝑡 (𝒓) + 𝑉 𝐻 (𝒓) + 𝜇xc
2.2.5
𝐻
2
ここで𝑉 𝑒𝑥𝑡 (𝒓)は原子核による外場からのポテンシャル、𝑉 𝐻 (𝒓)は電子間のクーロン相互作
𝜎
用によるポテンシャルである。𝜇xc(𝒓)は交換相関ポテンシャルであり、唯一関数形の不明

なポテンシャルである。
交換相関ポテンシャルの形が決まれば、電子密度𝑛𝜎 (𝒓)を与える式∗ (𝒓)𝜑𝑛𝜎 (𝒓) = ∑ 𝜑𝑘,𝜎𝑘,𝜎 (𝒓)2.2.6

𝑘と連立して KS 方程式を解くことが出来る。ここで𝑛𝜎 (𝒓)は𝜎のスピンをもつ電子の密度である。
KS 方程式は、多電子問題を１電子問題に還元した式で有り、系が基底状態にある時のエネルギーを求めるには SCF 法に従って KS 方程式を繰り返し解く必要がある。次項ではSCF 法について紹介する。

14

2.2.4項 Self-Consistent Field 法
SCF (Self-Consistent Field) 法に従って KS 方程式を繰り返し解くことで基底状態のエネ
ルギーが得られる。以下でその流れを紹介する。まず原子配置から仮決めした電子密度
̂ 𝐷𝐹𝑇 を決定する。これを用いて KS 方程式を解き、𝑛𝜎 (𝒓)を用いて2.2.5 式のハミルトニアン𝐻
新たな 1 電子軌道関数𝜑𝑘,𝜎 (𝒓)を得る。さらにこれを用いてハミルトニアンを再構築し、再び KS 方程式を解く。これを KS 方程式の１電子エネルギー固有値が全て収束するまで繰り返すことによって、与えられた原子配置に対する基底状態の全電子のエネルギー固有値が得られる。以上が SCF 法の流れである。
𝑁

𝐸

𝐷𝐹𝑇

2.2.7

̂ 𝐷𝐹𝑇 |𝜑𝑘 ⟩= ∑⟨𝜑𝑘 |𝐻
𝑘=1

次項ではこうして得られた全電子のエネルギー固有値から原子核に働く力を計算する方法 (Hellmann-Feynman の定理) について紹介する。

2.2.5項 Hellmann-Feynman の定理
全電子から原子核𝛼への作用𝑭𝑒𝛼 は次式で計算できる。
𝑑𝐸({𝑹𝛼 })𝑭𝑒𝛼 = −
2.2.8
𝑑𝑹𝛼
ここで𝑹𝛼 は原子核 𝛼の座標、𝐸({𝑹𝛼 })は全原子核の配置が{𝑹𝛼 }である系での全電子のエネルギー固有値であり、2.2.7 式の𝐸 𝐷𝐹𝑇 に相当する。𝑭𝑒𝛼 は次のように変形できる。
𝑭𝑒𝛼 = −

𝑑𝐸({𝑹𝛼 })𝑑̂ 𝐷𝐹𝑇 |𝜑𝑘,𝜎 ⟩∑⟨𝜑𝑘,𝜎 |𝐻
=−
𝑑𝑹𝛼
𝑑𝑹𝛼

2.2.9

𝑘,𝜎

̂ 𝐷𝐹𝑇
𝜕𝜑
𝜕𝜑
̂ 𝐷𝐹𝑇 |𝜑𝑘,𝜎 ⟩ + ⟨𝜑𝑘,𝜎 | 𝜕𝐻
̂ 𝐷𝐹𝑇 | 𝑘,𝜎 ⟩}|𝜑𝑘,𝜎 ⟩ + ⟨𝜑𝑘,𝜎 |𝐻= − ∑ {⟨ 𝜕𝑹𝑘,𝜎 |𝐻
𝜕𝑹
𝜕𝑹
𝑘,𝜎

𝛼

𝛼

𝛼

∗
̂ 𝐷𝐹𝑇
𝜕𝜑
𝜕𝜑
̂ 𝐷𝐹𝑇 |𝜑𝑘,𝜎 ⟩ + ⟨𝜑𝑘,𝜎 | 𝜕𝐻
̂ 𝐷𝐹𝑇 |𝜑𝑘,𝜎 ⟩ }|𝜑𝑘,𝜎 ⟩ + ⟨ 𝑘,𝜎 |𝐻= − ∑ {⟨ 𝜕𝑹𝑘,𝜎 |𝐻
𝜕𝑹
𝜕𝑹
𝑘,𝜎

𝛼

𝛼

𝛼

∗
̂ 𝐷𝐹𝑇
𝜕𝜑𝑘,𝜎
𝜕𝜑𝑘,𝜎
𝜕𝐻
∗
= − ∑ {𝜀𝑘,𝜎 ⟨ 𝜕𝑹 |𝜑𝑘,𝜎 ⟩ + ⟨𝜑𝑘,𝜎 | 𝜕𝑹 |𝜑𝑘,𝜎 ⟩ + 𝜀𝑘,𝜎 ⟨ 𝜕𝑹 |𝜑𝑘,𝜎 ⟩ }
𝛼
𝛼
𝛼
𝑘,𝜎

̂ 𝐷𝐹𝑇
𝜕𝜑
𝜕𝜑
𝜕𝐻
= − ∑ {𝜀𝑘,𝜎 ⟨ 𝜕𝑹𝑘,𝜎 |𝜑𝑘,𝜎 ⟩ + ⟨𝜑𝑘,𝜎 | 𝜕𝑹 |𝜑𝑘,𝜎 ⟩ + 𝜀𝑘,𝜎 ⟨𝜑𝑘,𝜎 | 𝜕𝑹𝑘,𝜎 ⟩}𝛼𝑘,𝜎= − ∑ 𝜀𝑘,𝜎
𝑘,𝜎

𝛼

𝛼

̂ 𝐷𝐹𝑇
𝜕
𝜕𝐻
⟨𝜑𝑘,𝜎 |𝜑𝑘,𝜎 ⟩ − ∑ ⟨𝜑𝑘,𝜎 ||𝜑𝑘,𝜎 ⟩
𝜕𝑹𝛼
𝜕𝑹𝛼
𝑘,𝜎

= − ∑ ⟨𝜑𝑘,𝜎 |
𝑘,𝜎

̂ 𝐷𝐹𝑇
𝜕𝐻
|𝜑𝑘,𝜎 ⟩
𝜕𝑹𝛼2.2.1015̂ 𝐷𝐹𝑇 がエルミート演算子であることと、𝜑𝑘,𝜎 がその固有で状態であること、𝜀𝑘,𝜎 が途中で𝐻その固有値であり実数であること、⟨𝜑𝑘,𝜎 |𝜑𝑘,𝜎 ⟩ = 1であることなどを用いた。これらの仮定のもとで次式が成り立つことを Hellmann-Feynman の定理と言う。
𝑭𝑒𝛼 = − ∑ ⟨𝜑𝑘,𝜎 |
𝑘,𝜎

̂ 𝐷𝐹𝑇
𝜕𝐸({𝑹𝛼 })𝜕𝐻|𝜑𝑘,𝜎 ⟩ (= −
)
𝜕𝑹𝛼
𝜕𝑹𝛼

2.2.11

前項まで議論に従って全電子のエネルギー固有値を計算できていれば、この式で原子核に働く力を計算出来る。

2.3節

Artificial Neural Network (ANN) ポテンシャル

以下では ANN ポテンシャルの仕組み (2.3.1 項) や構築方法 (2.3.2 項)、ANN-MD 法での使われ方 (2.3.3 項)、Ewald 法併用時の仕組み (2.3.4 項)について紹介する。

2.3.1項 ANN ポテンシャルとはANN ポテンシャルとは、ANN という任意の関数形を近似できるネットワークを用いた
原子間ポテンシャル関数である。１原子周りの原子配置を基にした値を入力すると、その原子がもつポテンシャルエネルギーを出力する。FPMD のデータを学習させれば FPMD
法と同程度に高精度な結果を古典ポテンシャル程度の低い計算コストで与えてくれるとい
う大きな利点がある。つまり、役割としては FPMD 法と同程度の精度をもつ古典ポテンシ
ャルだと言える。以下では ANN ポテンシャルを構成するパーツや ANN ポテンシャルの入力値、出力値などについて紹介する。
そもそも ANN とは人間の脳神経回路を模した数理モデルのことであり、様々な関数形を近似的に再現できると言われている (万能近似定理)。ANN は汎化性能の高さから現在画像認識や商品の売上予想、チャットボットなど様々な場面で利用されており、その主な働きは次の２つに分けられる。１つは人間がしていた単純作業を代替するという働きである。
もう１つは人間が把握しきれないような膨大なデータから特徴を抽出して、新たな入力に対する適切な予想を与えるという働きである。我々が用いる ANN ポテンシャルはどちらかというと後者の働きをしている。
ANN と一言で言っても様々な種類が有るが、我々が用いる ANN は順伝搬型と呼ばれるものであり、以降の議論で ANN といった場合には基本的にこの形のものを指す。順伝播型は最も単純な形をした ANN であり、その名の通り入力された情報をワットワークの出力側に伝搬する仕組みになっている。具体的な形を図 2.1 に示してある。以下で我々が用いている ANN ポテンシャルの各パーツの説明などをする。

16

図 2.1 ANN の概念図
順伝播型の ANN は図 2.1 のように入力層と隠れ層、出力層からなっている。入力層は入力値が入る層である。隠れ層は ANN の複雑さに直結する層であり、数が多ければ多いほど ANN の汎用性能は上がる。出力層は ANN の最終的な出力値を与える層である。各層にはノードというスカラー量の入れ物が有り、ノード数の多さも ANN の汎化性能に影響する。本研究では、出力層のノードは１つで、それ以外の層のノードは複数に設定してある。
𝑘,𝑙
各ノードは隣接する層のノードと重み𝑤𝑖,𝑗で結ばれていて、数値のやり取りの際にこの重みが使われる。具体的には𝑘 層の𝑖番目のノードの数値を𝑙 (= 𝑘 + 1)層の𝑗番目のノードに送𝑘,𝑙る際に重み𝑤𝑖,𝑗をかけて送るという使われ方をする。また各隠れ層にはバイアスノード 𝑏𝑁𝑘 𝑘というものがあり、これは唯一前の層の影響を受けないノードとして利用されている。こ
こで𝑁𝑘 は𝑘層のノードの数であり、各層の最終ノードとしてバイアスノードが来る。このバイアスノードは簡単な関数𝑦 = 𝑎𝑥 + 𝑏の切片𝑏のように、次の層への入力値の底上げ、底𝑘,𝑙
下げをする働きがあると言える。以上から𝑙層𝑗番目のノードへの入力値は 𝑏𝑁𝑘 𝑘 + ∑𝑖 𝑤𝑖,𝑗
⋅ 𝑦𝑖𝑘
𝑘,𝑙 𝑘
となり、ノードが持つ値は𝑓𝑗𝑙 ( 𝑏𝑁𝑘 𝑘 + ∑𝑖 𝑤𝑖,𝑗𝑦𝑖 )となる。この𝑓𝑗𝑙 (⋅)は活性化関数と呼ばれるものであり、本研究では以下のような非線形関数を用いた [29]。
𝑓 𝑙 (𝑥) = a tanh (b𝑥)a = 1.7195 , b = 2/3

2.3.1
2.3.2

非線形関数を用いているのは二階微分しても傾きが消失しないからである。この点は後述
する力の学習のために重要である。
入力値としては１原子周りの原子位置を基にした対称関数 [25, 26]の値を用いる。対称関数とは系全体の並進や回転に対する対称性がある関数のことである。同様の対称性をエネ
ルギーも持っているので、対称関数と対応付ければ情報量の無駄が少なくて済む。一方で、原子座標をそのまま入力値として用いた場合には、系全体の並進や回転に対してエネルギ

17

ーが変わらないのに入力値が大きく変わってしまうため、情報量の無駄が大量に生じてしまう。この違いが原子座標ではなく対称関数を入力値として用いる大きな理由である。他にも対称関数を用いれば、教師データとサイズの異なる系に ANN ポテンシャルを利用出来るなどという利点も有る。また、対称関数には Behler−Parrinello 型 [25]と Chebyshev
型 [26]などの複数の種類が有る。各型の詳細については付録 Error! Reference source notfound.、Error! Reference source not found.で、またこれらによる結果の違いについては3.8
節で紹介している。
１つの ANN ポテンシャルからの出力値は構造𝜎中の𝑖原子がもつポテンシャルエネルギー
ANN
ANN
𝜀𝜎,𝑖
であり、構造𝜎中の全ポテンシャルエネルギー𝐸𝜎𝐴𝑁𝑁 は𝐸𝜎𝐴𝑁𝑁 = ∑𝑖∈𝜎 𝜀𝜎,𝑖で与えられる。
FPMD
FPMD で計算出来るものは𝜀𝜎,𝑖ではなく𝐸𝜎FPMD であるため、エネルギーの精度の検証にANN
は𝐸𝜎𝐴𝑁𝑁 が使われる。また、𝑖原子に働く力は𝜀𝜎,𝑖ではなく𝐸𝜎𝐴𝑁𝑁 を𝑖原子の座標で微分することによって得られる。
十分な精度の ANN ポテンシャルを構築するには様々なパラメータの最適化が必要である。
𝑘,𝑙
１番代表的なものは重み𝑤𝑖,𝑗の最適化 (学習) であり、これはアルゴリズムによって自動的に行われる。学習については次項で説明する。他には入力値を与える対称関数の種類やそれに含まれるパラメータ、隠れ層の層数やノード数、活性化関数の形、損失関数のパラメ
ータ、最適化の手法、学習のスケジュールなど我々人が試行錯誤して調整するべきパラメータ（ハイパーパラメータ）がある。十分な精度の ANN ポテンシャルを構築するにはこ
れらのハイパーパラメータも最適化する必要があるので、本研究ではそこにも重点を置いている。

2.3.2項 学習
ANN ポテンシャルを利用する前に、学習によって予測精度を高める必要がある。学習と𝑘,𝑙は ANN ポテンシャルによる出力と教師データとの誤差を小さくするために重み{𝑤𝑖,𝑗}をア

ルゴリズムに従って調整することである。ここで言う誤差は損失関数というもので測られ、損失関数の中には我々人が調整するパラメータ（ハイパーパラメータ）もある。このハイ
パーパラメータはエネルギーや力、圧力などのうちどの物理量を重点的に学習させるかを決めるものである。また、これを含めた全てのハイパーパラメータを固定したとしても、𝑘,𝑙
重み{𝑤𝑖,𝑗}の初期値が違えば異なる ANN ポテンシャルが構築される。したがって、初期重𝑘,𝑙み{𝑤𝑖,𝑗}ini の異なる複数の ANN ポテンシャルを構築して結果を議論するのが妥当である。

以下では、これら学習の詳細について紹介する。
学習とは ANN ポテンシャルによる出力と教師データとの誤差を小さくするために重み
𝑘,𝑙
{𝑤𝑖,𝑗
}を調整することである。この誤差は次式のような損失関数で表され、実質的には損失関数の値を小さくするように学習が進んでいく。
𝑝𝐸
𝑝𝐹
𝑝𝑊
(𝛥𝑊)2
ℒ = (𝛥𝐸)2 + (𝛥𝐹)2 +
2
2
2

2.3.3

18

𝑁

1
1
(𝐸 ANN − 𝐸𝜎FPMD )2
𝛥𝐸 = √ ∑ 2
𝑁𝑁atom𝜎 𝜎

2.3.4

𝜎=1

𝑁atom𝜎𝑁

1
1
FPMD )2∑ (𝑭ANN
𝛥𝐹 = √ ∑𝜎,𝑖 − 𝑭𝜎,𝑖𝑁
3𝑁atom𝜎
𝜎=1

2.3.5

𝑖=1

𝑁

6

2

ANN
FPMD
𝑊𝜎,𝑖
− 𝑊𝜎,𝑖1∑∑ (𝛥𝑊 = √)
6𝑁
𝑁atom𝜎2.3.6

𝜎=1 𝑖=1ANN
ここで𝐸𝜎ANN , 𝑭ANN
𝜎,𝑖 , 𝑊𝜎,𝑖 はそれぞれ構造値𝜎に対して ANN ポテンシャルにより計算されたポテンシャルエネルギー、原子𝑖 に働く力、 Virial テンソルの第 𝑖 成分を表しており、
FPMD
𝐸𝜎FPMD , 𝑭FPMD, 𝑊𝜎,𝑖はそれぞれ構造値𝜎に対して FPMD により計算されたポテンシャルエ
𝜎,𝑖ネルギー、原子𝑖に働く力、Virial テンソルの第𝑖成分を表している。また、𝑁は構造の数である。𝛥𝐸, 𝛥𝐹, 𝛥𝑊はそれぞれエネルギー、力、Virial テンソルの教師データからの誤差を表しており、RMSE (Root Mean Square Error) と呼ばれる。これが各物理量の学習精度の指標として使われる。Virial テンソルに関しては、圧力の RMSE である𝛥𝑃に変換して評価することが多い。この𝛥𝑃は圧力の中でもポテンシャルエネルギーの寄与に依る部分の RMSE で
ある。
𝛥𝑊
2.3.7
𝛺𝜎
ANN
FPMD
𝑊𝜎,𝑖
や𝑊𝜎,𝑖
は Virial の定理から計算されるものであり、アンサンブル平均ではない [30]。
𝛥𝑃 =また、𝑝𝐸 , 𝑝𝐹 , 𝑝𝑊 をそれぞれエネルギーの学習比重、力の学習比重、Virial の学習比重もし
くは圧力学習の比重と呼び、𝑝𝐸 , 𝑝𝐹 , 𝑝𝑊 それぞれが他より大きい場合の学習のことを、エネ
ルギー重視の学習、力重視の学習、Virial 重視の学習もしくは圧力重視の学習と呼ぶ。さらに3.3 節や3.4 節で頻出する「Virial 学習」という言葉は𝑝𝑊 ≠ 0の学習のことを指す。
𝑘,𝑙𝑘,𝑙重みの最適化は𝜕ℒ/𝜕𝑤𝑖,𝑗に基づいてℒが小さくなる方向に𝑤𝑖,𝑗を調整することによって行われる。この最適化の数学的手法には種類があり、我々は Levenberg-Marquardt 法 [31, 32]を用いている。経験的にこれが最も効率的な最適化を可能にすると分かっている。
学習対象が複数ある場合は損失関数の係数𝑝𝐸 , 𝑝𝐹 , 𝑝𝑊 の調整が必要である。この係数は、𝛥𝐸, 𝛥𝐹, 𝛥𝑊の次元とスケールを調整する働きがある。特にスケールを調整する働きによって、重点的に学習させる物理量を選べる。経験的には、適切な調整によって RMSE をより小さく出来るということが分かっている。この調整に関しては3.3 節で詳細に議論している。
他にも様々なハイパーパラメータによって学習の精度は左右される。ハイパーパラメー
タとは我々人間が直接調整するパラメータのことであり先述の𝑝𝐸 , 𝑝𝐹 , 𝑝𝑊 もその１つである。
他には学習対象とする物理量の選択 (3.4 節で議論している)、Ewald 法 (クーロン相互作用を計算する手法)との併用の有無 (3.5 節で議論している)、原子配置からの特徴抽出に用いる対称関数の種類やパラメータ (3.8 節で議論している)、隠れ層の層数やノード数 (私の学19

部の卒業論文 [24]で議論している)などである。
以上のようなハイパーパラメータを全て固定したとしても、得られる ANN ポテンシャル𝑘,𝑙𝑘,𝑙は初期重み{𝑤𝑖,𝑗}ini によって微妙に変わってくるので、初期重み{𝑤𝑖,𝑗}ini の異なる複数のANN ポテンシャルを用いて出した結果について議論するのが妥当である。このような初期𝑘,𝑙
重み依存性がある原因は、ANN ポテンシャルのパラメータの多さに起因して重み{𝑤𝑖,𝑗}の局所解が大量にあることである。局所解ではあるものの、十分な学習が行われていれば実用上問題ない精度の ANN ポテンシャルを構築出来る場合が殆どである。そう言っても、
初期重み依存性による ANN ポテンシャルの違いは RMSE にも現れるし、ANN ポテンシャ
ルを用いた MD (ANN-MD) 法の結果にも現れる。そのためハイパーパラメータの調整過𝑘,𝑙程においても、ANN-MD 法で物理量を計算する場合においても、初期重み{𝑤𝑖,𝑗}ini の異なる複数の ANN ポテンシャルについての結果を用いて議論すること（初期重み依存性の考慮）が重要である。本研究では一貫して初期重み依存性を考慮した結果を用いている。
ちなみに本研究では ANN ポテンシャルの構築に、Aenet [33]というソフトウェアパッケ
ージを用いている。しかし、そのままではエネルギーの学習しか出来なかったので、原子に働く力や Virial テンソルを学習できる形に拡張したものを用いている [23]。
以上のように学習とは ANN ポテンシャルによる出力と教師データとの誤差を小さくする
𝑘,𝑙
ために、ANN ポテンシャルの重み{𝑤𝑖,𝑗}を調整することであり、その精度を左右するハイパーパラメータが沢山ある。また、ANN ポテンシャルの初期重みによる違いも生じるので、初期重みの異なる複数の ANN ポテンシャルを用いて議論する必要が有る。

2.3.3項 ANN-MD
ANN-MD とは経験的ポテンシャルとして ANN ポテンシャルを用いた古典 MD のようなものであり、低い計算コストで FPMD 計算と同程度の精度の結果を得られるという長所がある。その一方で、ANN ポテンシャルの初期重み依存性を考慮する必要があるなどの短所も有る。以下では ANN-MD の簡単な仕組みや特徴についてについて紹介する。
ANN-MD の流れはほぼ古典 MD シミュレーションと同じであり、大きく異なるのはポテ
ンシャルエネルギーや各原子に働く力の計算に経験的ポテンシャルではなく ANN ポテン
シャルを用いている点である。ANN ポテンシャルによるポテンシャルエネルギーなどの計算の流れは2.3.1 項に記したとおりである。簡単にもう一度記すと次のようになる。ある時
刻での各原子周りの原子配置は対称関数に変換されて ANN ポテンシャルの入力値になる。
この入力値を基に ANN ポテンシャルで各原子がもつエネルギーを計算する。これらの総和として全系のポテンシャルエネルギーを得る。全系のポテンシャルエネルギーを各原子座標について微分することで各原子に働く力を計算できる。以上のようにして得られた力に基づいて、MD シミュレーションが行われる。
ANN-MD の利点は計算精度と計算速度にある。計算精度については、FPMD を教師デ20

ータとして十分な学習を行った場合には FPMD と同程度の精度の結果を出せるという点で優れている。計算速度については、経験的ポテンシャルを用いた MD (古典 MD) に比べたら遅いものの FPMD よりも非常に速いという点で優れている。例えば、Ag2Se の場合にANN-MD は FPMD より約１万倍速く計算できると分かっている。以上２つの利点から、大規模な系や長時間の計算を必要とする物理量の計算まで MD シミュレーションの適用範
囲を広げてくれる存在として期待されている。
ANN-MD から得られる物理量にも ANN ポテンシャルの初期重み依存性が反映される。
なぜなら１つの原子配置に対して計算される力は初期重みの異なる ANN ポテンシャル間で微妙に異なり、時間発展させる度にその微妙な違いが積算され、それが最終的な物理量にも影響するからである。例えば、3 章にて格子熱伝導度がプロットされたすべてのグラ
フでは、ANN ポテンシャルの初期重みによるばらつき表すためにエラーバーをつけてある。
他にも動径分布関数が微妙に異なったり、ANN-MD で破綻が起きるものとそうでないものができたりする。このように ANN-MD においても ANN ポテンシャルの初期重み依存性を考慮するために、初期重みの異なる複数の ANN ポテンシャルによる結果について議論するのが適当である。この分計算コストが高くなるので、これは短所であると言える。

2.3.4項 Ewald 法との併用
Ag2Se はイオン性が高いのでクーロン相互作用を考慮する必要があると考えられる。一方で ANN ポテンシャルは長くても原子周りの半径 9 Å 程度までの原子配置しか考慮しない
ので、ANN ポテンシャル単体では長距離相互作用であるクーロン相互作用を考慮できているかの懸念が残る。この問題の対処法としては ANN ポテンシャルと、クーロン相互作用を計算する手法である Ewald 法を併用するという方法が有る。以下では ANN ポテンシャ
ルが考慮出来る空間範囲の制限、ANN ポテンシャルと Ewald 法の併用法、併用する場合の学習法について紹介する。
ANN ポテンシャルにはカットオフ半径があり、それ以遠の原子配置は考慮されないため、
長距離相互作用をうまく再現できないと考えられる。カットオフ半径以内に入ってくる原子の数はおおよそ𝑟 3 に比例して増えるので、長くしすぎる対称関数などの計算コストが高くなる。そのため、カットオフ半径は5~9 Åの範囲で設定されることが多く、今回の例では 9 Å と設定している。このような制限のため、クーロン相互作用のような長距離相互作用を考慮するべき系に対しては ANN ポテンシャルに加えて補正が必要になると考えられる。
ANN ポテンシャルを用いつつ、クーロン相互作用を考慮する方法に ANN ポテンシャルと Ewald 法を併用する手法がある。ANN ポテンシャルで短距離的な相互作用を計算し、Ewald 法でクーロン相互作用を計算するという手法である。Ewald 法は各元素に与えた固定電荷にしたがってクーロン相互作用を計算してくれる。固定電荷を与えるため、局所的

21

な原子配置による各原子の電荷のゆらぎまでは考慮出来ないが、クーロン相互作用の近似的な考慮としては十分である。
Ewald 法を併用する際に ANN ポテンシャルに学習させるものは、もともとの教師データ
SRから Ewald 法による効果を引いた残りの{𝐸𝜎SR }, {𝑭SR𝜎,𝑖 }, {𝑾𝜎 }である。ここで SR は Short

Range の略である。これら SR の教師データは以下の手順に従って作られる。
1.

教師データから各元素の有効電荷を抽出し、Ewald 法で使う電荷として設定する。
我々は Mulliken 電荷 [34, 35, 36]を教師データから抽出した。

2.

その有効電荷に基づいた Ewald 法の計算によって、教師データの各構造{𝒓𝜎,𝑖 }に対す
る全ポテンシャルエネルギー{𝐸𝜎Ewald }と各原子に働く力{𝑭Ewald}、Virial テンソル𝜎,𝑖{𝑾Ewald
}を得る。
𝜎

3.

}, {𝑾FPMD}, {𝑾Ewald
もともとの教師データ{𝐸𝜎FPMD }, {𝑭FPMD}と{𝐸𝜎Ewald }, {𝑭Ewald
}との差
𝜎
𝜎
𝜎,𝑖
𝜎,𝑖
SR
分をそれぞれの物理量で取ることによって{𝐸𝜎SR }, {𝑭SR𝜎,𝑖 }, {𝑾𝜎 }を得る。

この点以外は通常の学習と同様である。

2.4節

熱伝導度計算以下では MD シミュレーションを用いて熱伝導度を計算する２つの方法について紹介する。１つは平衡分子動力学 (Equilibrium Molecular Dynamics, EMD) 法によるもので、もう １ つ は 非 平 衡 分 子 動 力 学 (Homogeneous Non-Equilibrium Molecular Dynamics,
HNEMD) 法 (今回用いた方法は特に摂動分子動力学 (perturbed MD, pMD) 法と言われる)によるものである。その後にそれらの計算に用いる熱流式について紹介する。

2.4.1項 平衡分子動力学法による計算方法
平衡分子動力学 (Equilibrium Molecular Dynamics, EMD) 法による熱伝導度計算は、長時
間 (Ag2Se の場合はおよそ 2 ns 以上) の平衡分子動力学計算のデータとグリーン・久保(Green-Kubo, GK) 公式を用いて行われる。GK 公式は以下のような形をしている。
𝜏𝛺∫⟨𝒋 (𝑡) ⋅ 𝒋𝑄 (0)⟩d𝑡
𝜅=2.4.1
3𝑘𝐵 𝑇 2 0 𝑄
ここで、𝛺 は体積、𝑘𝐵 はボルツマン定数、𝑇は温度、𝒋𝑄 (𝑡)時刻𝑡における熱流、〈⋅〉はアンサ
ンブル平均（シミュレーション時間に対する平均）、を表す。⟨𝒋𝑄 (𝑡) ⋅ 𝒋𝑄 (0)⟩は熱流の自己相関関数(Heat Current Auto-Correlation Function, HCACF)と呼ばれる。したがって、GK 公式というのは熱流の自己相関関数の時間積分によって熱伝導度𝜅を求める式だと言える。𝜏
は理想的には無限大としたいが、それは無理なので有限な値となっている。実際には 5 ps
程度で十分収束した値が得られる。
EMD 法においては熱流を
22

𝒋𝑄 = 𝒋Se + 𝒋Ag = 𝒋𝐸 + 𝒋𝐼のように分解することによって、式2.4.1 から熱伝導度を
𝜅Total = 𝜅Se−Se + 𝜅Se−Ag + 𝜅Ag−Se + 𝜅Ag−Ag = 𝜅𝐸−𝐸 + 𝜅𝐸−𝐼 + 𝜅𝐼−𝐸 + 𝜅𝐼−𝐼

2.4.2
2.4.3

のような元素のペアや力学的エネルギーの項 (E) と原子間相互作用の項 (I) のペアの寄与に分けられる。3 章での熱伝導度の成分はこの式に基づいて計算されている。
EMD 法による計算の長所は以下の通りである。
⚫

Ag2Se のような異常な振る舞いをする原子を含む結晶 [11]に対しても熱伝導度を計
算出来る。

⚫

結晶の異方性の有無に関わらず１回の計算で熱伝導度を計算できる。

⚫

平衡状態に対して MD を行えば良いので、ANN-MD でも比較的安定して MD を実行しやすい。

EMD 法による計算の短所は以下の通りである。
⚫

十分収束させるには 2 ns 以上の非常に長時間の MD シミュレーションが必要。

⚫

NVE アンサンブルでする必要が有るので、NVT アンサンブルなどに比べて温度がばらつき易い。

2.4.2項 非平衡分子動力学法による計算方法
非平衡分子動力学 (Homogeneous Non-Equilibrium Molecular Dynamics, HNEMD) 法[37]による計算では、系に摂動を加えながら MD を行う。摂動が𝑥方向にのみかかっている場合を考えると、非平衡状態の熱流〈𝑗𝑥 〉neは以下の式で計算できる。
𝜏𝐹〈𝑗𝑥 〉ne = 〈𝑗𝑥 〉eq + ext ∫ 〈𝑗𝑥 (𝑡)𝑗𝑥 (0)〉 𝑑𝑡
2.4.4
𝑘𝐵 𝑇 0
ここで、𝐹ext は系に加える摂動である。平衡状態の熱流の平均〈𝑗𝑥 〉eq は０であるため、GK
公式とこれを比べると熱伝導度𝜅を与える式が得られる。
𝛺𝜅=lim 〈𝑗 〉𝑇𝐹ext 𝑡→∞ 𝑥 ne
HNEMD 法においては熱流を
𝒋𝑄 = 𝒋Se + 𝒋Ag = 𝒋𝐸 + 𝒋𝐼のように分解することによって、式2.4.5 から熱伝導度を
𝜅Total = 𝜅Se + 𝜅Ag = 𝜅𝐸 + 𝜅𝐼

2.4.5

2.4.6
2.4.7

のような元素ごとの寄与(𝜅Se と𝜅Ag )や力学的エネルギーの項 (𝜅𝐸 ) と原子間相互作用の項(𝜅𝐼 ) の寄与に分けられる。3.7 節での熱伝導度の成分はこの式に基づいて計算されている。
HNEMD 法による計算の長所は以下の通りである。
⚫

摂動が大きいほど熱伝導度の収束が速いのでシミュレーション時間が短くて済む。
具体的には EMD 方と比べて 1/2 ~ 1/3 程度まで短く出来る。

23

⚫

NVT アンサンブル計算できるので温度が比較的安定しやすい。

HNEMD 法による計算の短所は以下の通りである。
⚫

摂動に対する熱伝導度の線形応答領域を見つけないといけない。

⚫

摂動が弱すぎると収束までに時間がかかる。しかし長くても EMD 方の 1/2 程度で済む場合がほとんど。

2.4.3項 熱流束の表式
元素が𝑁type種類有る系での熱流束は次式で与えられる [38]。
𝑁atom

𝑁atom

𝑁atom

𝑁type

𝑁𝜇

𝑖=1

𝑖=1

𝑖=1

𝜇=1

𝑗=1

1
1
𝒋𝑄 = ( ∑ 𝑡𝑖 𝒗𝑖 + ∑ 𝜀𝑖 𝒗𝑖 + ∑ 𝑾𝑖 𝒗𝑖 ) − ∑ ℎ𝜇 ∑ 𝒗𝑗
𝛺
𝛺

2.4.8

ここで𝛺は体積、𝑁atom は系に含まれる粒子の数、𝑡𝑖 は原子𝑖の運動エネルギー、𝑣𝑖 は原子𝑖の速度、𝜀𝑖 は原子𝑖のポテンシャルエネルギー、𝑾𝑖 は原子𝑖の Virial テンソル、ℎ𝜇 は元素𝜇の部分エントロピー、𝑁𝜇 は元素𝜇の粒子数である。ANN ポテンシャルを用いた場合には𝜀𝑖 が𝜀𝑖ANN になる。式2.4.2 の𝒋𝑄 = 𝒋𝐸 + 𝒋𝐼 という寄与の分解において、𝒋𝐸 は(𝑡𝑖 + 𝜀𝑖 )を含んだ項であり力学的エネルギーの項と呼ばれる。𝒋𝐼 は𝑾𝑖 を含んだ項であり原子間相互作用の項と呼ばれる。𝑾𝑖 はポテンシャルエネルギーが２体によるものか、多体によるものかで形が変わり、その形によって熱伝導度が大きく変わるということが知られている [39, 40, 41]。以下
では ANN ポテンシャルに対する𝑾𝑖 の表式を紹介する。
対称関数が Behler−Parrinello 型の場合、原子あたりのポテンシャルエネルギーが𝜀𝑖ANN = 𝜀𝑖ANN ({𝒓𝑖𝑗 } , {𝒓𝑗𝑘 }
)
2.4.9
𝑗≠𝑖

𝑗,𝑘≠𝑖

であることから次式が𝑾𝑖 の正しい表式であると導かれ、これによる熱流束は正しいものになる。
𝜕𝜀𝑗ANNrig∑∑
𝑾𝑖 =𝒓𝑖𝑗 ⨂ ()
𝜕𝒓𝑘𝑖𝑗≠𝑖 𝑘≠𝑖

2.4.10

ここで rig は rigorous の略である。一方、この対称関数に対して圧力テンソルから自然に導かれる𝑾𝑖 は次式であり [42, 25]、これによる熱流束は誤ったものになる。
𝜕𝜀𝑖ANN
∑
∑
𝑾wro
=
𝒓
⨂
(
)
𝑖𝑗
𝑖
𝜕𝒓𝑗𝑘

2.4.11𝑗≠𝑖 𝑘≠𝑖rig

ここで wro は wrong の略である。これら𝑾𝑖 と𝑾wroから得られる熱伝導度の違いについ𝑖ては3.6 節で紹介している。
また、対称関数が Chebyshev 型の場合のポテンシャルエネルギーは𝜀𝑖ANN = 𝜀𝑖ANN ({𝒓𝑖𝑗 } )2.4.12𝑗≠𝑖rig
なので𝑘についての和が消えて、上の𝑾𝑖 と𝑾wroに対応する表式がそれぞれ次のようにな
𝑖
る。

24

𝑁𝑎𝑡𝑜𝑚rig
𝑾𝑖 =∑ 𝒓𝑖𝑗 ⊗
𝑗≠𝑖
𝑁𝑎𝑡𝑜𝑚

𝑾wro
= ∑ 𝒓𝑖𝑗 ⊗𝑖𝑗≠𝑖

𝜕𝜀𝑗𝐴𝑁𝑁𝜕𝒓𝑗𝑖
𝜕𝜀𝑖𝐴𝑁𝑁𝜕𝒓𝑗𝑖

2.4.132.4.1425

3章 計算条件・結果・考察
目的「ANN-MD を通して格子熱伝導度を高精度に見積もるための ANN ポテンシャル構
築手法や格子熱伝導度の計算手法の確立」を達成するために、適切な ANN ポテンシャル
構築条件や熱伝導度の計算条件を探った。結果として 300 K において 0.306±0.013
Wm−1 K −1 という、実験値(0.4~0.5 Wm−1 K −1 ) [1, 2]に近くてモデルに対して妥当な熱伝導度を得られるようになった。また、次のような計算条件が適切だと分かった。
⚫

初期重みの異なる ANN ポテンシャルを５つ以上用意する。

⚫

エネルギーに加えて力と Virial テンソルを含めた学習をする (3.4 節)。

⚫

学習比重は次の通り (3.3 節)。
（ここでEpochとは学習回数のことである。）

𝑝𝐹
𝑝= 103, 𝑤 = 10−3
𝑝𝐸
𝑝𝐸
𝑝𝐹
𝑝𝑤
20~40 Epoch:= 1,= 10−2𝑝𝐸
𝑝𝐸

0~20 Epoch:⚫

Ewald 法の併用は不要 (3.5 節)。

⚫

ANN ポテンシャルの多体効果を考慮した正しい熱流束を用いる (3.6 節)。

2.4.1
2.4.2

他にも以下のようなことが分かった。
HNEMD 法と EMD 法とで熱伝導度の値はそれほど変わらないものの、系に異方性⚫

があるなら計算コストの観点で HNEMD 法よりも EMD 法が効率的である (3.7 節)。
対称関数が Behler-Parrinello 型でも Chebyshev 型でも結果はほとんど変わらない⚫

(3.8 節)。
以下では基本的な計算条件を紹介した後で、上記の設定が適切だと言える根拠を紹介す
る。

3.1節

教師データ (DFT+U) の計算条件

DFT+U の教師データの計算条件は以下のとおりである。
⚫

密度汎関数理論(Density Functional Theory, DFT) [43]⚫

一般化密度勾配近似(Generalized Gradient Approximation, GGA) [44]⚫

GGA+ ハバード補正(GGA+U) [45]⚫

ハバード補正：𝑈 − 𝐽= 6.0 eV [46, 47]⚫

PAW (projector augmented wave) 法 [48, 49]⚫

原子数：384 原子(Ag256 個、Se128 個)⚫

アンサンブル：粒子数、圧力、温度一定アンサンブル（NPT アンサンブル）⚫

圧力：0 GPa⚫

温度：300 K⚫

シミュレーション時間：7.2 ps (3000 step 分)
26

⚫

時間幅：2.4 fs

教師データとして採用したのは、3000 step 分全てである。速度スケーリング部分は教師
データに入っていない。

3.2節

ANN ポテンシャルに関する各種計算条件

本項では3.3 節以降での ANN ポテンシャルの構築や MD 計算、熱伝導度計算に関する計算条件を紹介する。なお、本章の各節では本条件に加えてそれ以前の項の結果で適切だと判断された条件を追加したものを用いている。例えば、3.4 節でエネルギーと力に加えてVirial テンソルを学習させることが適切だという結論に至っているので、3.5 節からの結果は Virial テンソル込の学習をした場合の結果について議論している。
教師データについて✓

3.1 節で紹介した 3000 step 分全てを用いている。

ANN ポテンシャルの設定✓

対称関数: Behler-Parrinello 型✓

隠れ層: 2 層, 各 10 ノード✓

初期重み依存性の考慮：それぞれの条件に対して、異なる初期重みから構築した５つの ANN ポテンシャルによる結果について議論した。
✓

学習比重：3.3 節の比重調整から得られた比重を使用した。具体的には以下の値。

ここでEpochとは学習回数のことである。

𝑝𝐹
𝑝= 103, 𝑊 = 10−3
𝑝𝐸
𝑝𝐸
𝑝𝐹
𝑝𝑊
20~40 Epoch:= 1,= 10−2𝑝𝐸
𝑝𝐸
0~20 Epoch:

✓

3.2.1
3.2.2

Ewald 法の併用：有り(3.5 節まで。3.5 節の検証で Ewald 法は不要と判断され、それ以降では用いられていない。)✓

Ewald 法の電荷：Mulliken 電荷 [34, 35, 36]を用いた。具体的には Ag が約 0.371𝑒、Se が約-0.742𝑒。
ANN ポテンシャルを用いた MD の設定✓

粒子数: Ag256, Se128✓

アンサンブル: NVE アンサンブル (粒子数、体積、エネルギー 一定)✓

体積: 平均 0 GPa になる体積✓

温度: 300 K✓

シミュレーション時間: 約 2.4 ns (∆𝑡 = 2.4 fs で約 1,000,000 step)27✓

速度 (温度) スケーリングに関するスケジュール:

NVT（粒子数、体積、温度一定）アンサンブルで 10000 step (20 step 置きに 500 回)の速度スケーリングした後、同じく NVT アンサンブルで 10000 step スケーリング無しの MD を行った。その後に NVE アンサンブルに移行して、GK 公式に用いる MD デー
タを集めた。
✓

スーパーセルの形: ユニットセルを 8×2×2 で並べた非等方的なもの。

EMD 法による熱伝導度計算の設定✓

GK 公式の積分の上端 (次式の𝜏) : 12 ps(実際には 12 ps までの様々な𝜏に対する𝜅が出力される。
次の条件で述べるが、実質 5ps 付近の値をグラフにしている)𝜏𝛺∫⟨𝒋 (𝑡) ⋅ 𝒋𝑄 (0)⟩d𝑡
𝜅=
3𝑘𝐵 𝑇 2 0 𝑄✓
以降の章でグラフに用いている値のとり方:

3.2.3

まず１つの ANN-MD に対して 𝜏 = 4.5~5.5 psの間の𝜅の平均を取る。次に、𝜅 の ANNMD 間の平均と標準偏差を計算し、この値をグラフに用いている。
𝜏は理想的には無限大で有るため、大きいほうが良いと考えられている。しかし、今回の計算では𝜏 = 12ps ではなく𝜏 = 4.5~5.5 ps程度の小さい値を用いている。これは実質的に
𝜅は𝜏 = 4.5~5.5 ps程度でほとんど収束していて、これ以上の𝜏に対する値は収束値にノイズが乗った値に過ぎないと判断したからである。また、𝜏 = 4.5~5.5 psで平均をとっているのは、𝜅の𝜏に対する微小な振動を考慮してのことである。

3.3節

Virial 学習の比重調整過程

本節では圧力の再現に不可欠な Virial 学習においての学習比重 (𝑝𝐸 , 𝑝𝐹 , 𝑝𝑊 ) 調整過程を紹介する。結論から言えば、

𝑝𝐹
𝑝= 103, 𝑤 = 10−3
𝑝𝐸
𝑝𝐸
𝑝𝐹
𝑝𝑤
20~40 Epoch:= 1,= 10−2𝑝𝐸
𝑝𝐸

0~20 Epoch:

3.3.1
3.3.2

という比重がおおよそ妥当であると分かった（ここでEpochとは学習回数のことである）。
学習対象ごとに調整は必要だと考えられる。
学習比重調整の基本方針としては、前半で力を重点的に学習して、後半で力の精度を下げない程度にエネルギーや圧力を重点的に学習するという方針をとった。なぜなら、１構造あたりの教師データの数は力が圧倒的に多く (表 3.1)、ANN-MD において原子の動きに直接影響するのは力だからである。
１構造あたりの教師データ数28

エネルギー

１

力

３×原子数

圧力

６

表 3.1 １構造あたりの教師データ数
仮に教師データの少ないものを前半に学習して教師データの多いものを後半に学習すれば、後半に学習したものに引っ張られて前半の学習の意味が無くなる可能性が高い。現にエネルギーを前半に重点的に学習して、力を後半に重点的に学習した場合はエネルギーのRMSE が後半で上昇することが確かめられている [24]。本研究では以上のような理由から、力を重点的に学習した後、力を妥協しない程度にエネルギー、圧力を重点的に学習するという方針のもとで学習比重を調整した。
力を重点的に学習する際の比重𝑝𝐹 /𝑝𝐸 と𝑝𝑤 /𝑝𝐸 の調整過程について紹介する。力重視の場
合は力の RMSE が最初の 10 Epoch 以下で大きくが下がり、その後ごく僅かずつ緩やかに下るという傾向がある [24] ので学習の回数は 20 Epoch で固定した。力を重視すればよいのだから、エネルギーと圧力の比重を力の学習の邪魔にならない程度に下げたらよい。具
体的には、𝑝𝐹 /𝑝𝐸 を十分大きくして、𝑝𝑤 /𝑝𝐸 を十分小さくすれば良い。𝑝𝐹 /𝑝𝐸 に関してはエネルギーと力の学習の比重調整を行った先行研究 [24]において𝑝𝐹 /𝑝𝐸 = 103 もあれば十分だということが分かっている。この数値は圧力を含めた学習でも変わらないと考えられる
ためそのまま採用する。次に𝑝𝑤 /𝑝𝐸 を調整するために、𝑝𝐹 /𝑝𝐸 = 103 固定下での𝑝𝑤 /𝑝𝐸 に対する力の RMSE の変化を見た（図 3.1）。

図 3.1 𝑝𝑊 /𝑝𝐸 = 10−6 ~102 での RMSE の平均値の変化。シンボルとエラーバーはそれぞれ初期重みの異なる５つの ANN ポテンシャルに対する RMSE の平均と標準偏差を示している。

29

初期値による依存性はあるものの𝑝𝑊 /𝑝𝐸 = 10−3 程度で十分収束していたのでこれを採用した。ちなみに、先程述べたように [24]を参考として学習回数を 20 Epoch に固定しているが、改めてこれは十分な学習回数であると分かる。なぜなら 𝑝𝑊 /𝑝𝐸 = 10−6 , 10の場合のRMSE の Epoch 変化（図 3.2 と図 3.3）の最後の 5 Epoch 間の力の RMSE の変化幅が、
ANN ポテンシャルの初期重みによって違いは有るものの図 3.1 で見ている値に対して小さいからである。結論としては、以下の比重で 0~20 Epoch 程度の学習をすれば十分だと分
かった。

𝑝𝐹
𝑝𝑊
= 103 ,= 10−3
𝑝𝐸
𝑝𝐸

3.3.3

図 3.2 𝑝𝑊 /𝑝𝐸 = 10−6 に対する RMSE の変化。初期重みの異なる５つの ANN ポテンシャル
に対する結果を示している (try01~try05)。挿入図では力の RMSE を拡大して表示してある。

30

図 3.3 𝑝𝑊 /𝑝𝐸 = 10に対する RMSE の変化。初期重みの異なる５つの ANN ポテンシャルに対する結果を示している (try01~try04)。挿入図では力の RMSE を拡大して表示してある。

力の学習後に圧力を重点的に学習する際の比重𝑝𝐹 /𝑝𝐸 、𝑝𝑊 /𝑝𝐸 の調整過程について紹介する。力の学習後の ANN ポテンシャルに対して、𝑝𝐹 /𝑝𝐸 = 103 と固定したまま𝑝𝑊 /𝑝𝐸 を
10−2 ~10の範囲で動かした結果、RMSE は図 3.4 通りに変化した。

31

図 3.4 𝑝𝑊 /𝑝𝐸 = 10~104 に対する RMSE の変化。挿入図では力の RMSE を拡大して表示してある。
圧力の RMSE は𝑝𝑊 /𝑝𝐸 の上昇に伴って下がっていき、𝑝𝑊 /𝑝𝐸 = 10で約2.5 × 10−3 GPaとなってほぼ底打ちになっている。一方で、力の RMSE は𝑝𝑊 /𝑝𝐸 = 10−2 では上昇しないが、それ以上の𝑝𝑊 /𝑝𝐸 に対しては上昇している。したがって圧力と力のどちらの精度を優先する
べきか決める必要がある。圧力の RMSE が小さいに越したことは無いが、0 GPa を再現できればよく、教師データも 0.5 GPa 程度のゆらぎは持っているのでΔ𝑃 = 1.0 × 10−2 GPa程度の精度もあれば十分である。この程度の精度は𝑝𝑊 /𝑝𝐸 = 101 でも達成できているので、この値を採用した。この比重であればエネルギーの RMSE も十分収束しているので圧力とエネルギーを同時に学習できる比重だと言える。しかし、エネルギーと力だけを学習する
際の比重に比べるとエネルギーとの比重が小さすぎて、ANN ポテンシャルの初期重みによってはエネルギーの RMSE が収束しない可能性があるので、以降の計算では𝑝𝐹 /𝑝𝐸 =10−2 ~10−3 とした。このようにエネルギーの比重を大きくしてもエネルギーと力だけを学習した場合 [24] と同様に、エネルギー以外の物理量の RMSE は大して上昇しなかった。
圧力を重点的に学習する際の比重調整においては ANN ポテンシャルの初期重み依存性を考慮していないものの、後の計算結果から今回の結論で問題ないことが分かっている。今回採用した比重を用いて５つの異なる初期重みに対する ANN ポテンシャルを構築したと
ころ、RMSE の平均と標準偏差は表 3.3 のようになった。標準偏差が十分小さく、ANN ポ
テンシャルの重みを変えても図 3.4 の傾向が変わらないと予想出来る。例えば、図 3.4 で𝑝𝑊 /𝑝𝐸 = 102 と𝑝𝑊 /𝑝𝐸 = 101 の最終的な力の RMSE の差が1.0 × 10−3であるのに対して、表32

3.2 での力の標準偏差は10−5程度の桁の数値であり十分小さい。

Train

Test

平均値

標準偏差

平均値

標準偏差

∆𝐸 (eV/atom)

6.86 × 10−4

1.06 × 10−5

6.75 × 10−4

3.69 × 10−5

∆𝐹 (eV/Å)1.93 × 10−2

3.60 × 10−5

1.93 × 10−2

5.63 × 10−5

∆𝑃 (GPa)7.32 × 10−3

9.79 × 10−5

7.41 × 10−3

1.97 × 10−4

表 3.3 今回採用した学習比重で 5 つの異なる初期重みに対する ANN ポテンシャルを構築したときの RMSE の平均と標準偏差
以上のような計算結果・考察から最終的に次の比重が最適であるとした。

𝑝𝐹
𝑝= 103 , 𝑊 = 10−3
𝑝𝐸
𝑝𝐸
𝑝𝐹
𝑝𝑊
20~40 Epoch:= 1,= 10−2𝑝𝐸
𝑝𝐸

0~20 Epoch:

3.3.4
3.3.5

これはあくまで Ag2Se に対して調べた結果なので、物質によって多少の調整が必要になる
可能性はある。以降の議論で出てくる ANN ポテンシャルの構築にはほぼこの比重を用い
ている。

3.4節

Virial 学習の効果

エネルギーと力だけを学習対象とした場合では圧力の再現精度が低いという問題があった。そこで圧力も含めた学習 (Virial 学習) を行った。学習比重としては前項で得られた比重を用いた。結果としては、Virial 学習をすることで圧力の RMSE が約 10-2 GPa となり、圧力の再現に十分な精度が得られた。また、圧力の総和は熱流束を通して熱伝導度と関係があるため、圧力の精度の向上には熱伝導度の ANN ポテンシャル間でのばらつきを小さくする効果が期待されている [23]。しかし今回調べている Ag2Se, 300 K の系では、そのよ
うな改善は見られなかった。以下では結果を詳細に紹介する。
エネルギーと力を学習させることで動径分布関数などの再現は出来ていたが、圧力を再現できていなかった。具体的には圧力の RMSE は約 1 GPa であり、実際に MD をした結果から得られる圧力も図 3.5 の通りに教師データのものから約 1 GPa ずれていた。これでは ANN-MD によって教師データと同じ熱力学的状況を再現できているとは言えない。そこで圧力 (実際には Virial テンソル) を含めた学習を行った。

33

図 3.5 圧力を学習していない ANN ポテンシャルによる圧力（実線）と教師データの圧力（点破線）の比較。黒が𝒙方向、赤が𝒚方向、緑が𝒛方向の圧力を表す。
前項の比重調整方法に従って Virial 学習を行った結果、他の物理量の精度を大きく下げる
こと無く圧力の精度を十分に改善出来た。具体的には圧力の RMSE が約 10-2 GPa と十分小さくなり (図 3.6)、MD 中に出力される圧力も教師データと同様の 0±0.5 GPa 程度になった。一方、Virial 学習によってエネルギーの RMSE が上昇したが (図 3.6)、部分動径分布関数、Mean Squared Displacement (MSD) などの基本的な物理量の再現は Virial 学習なしの場合と同様に出来ていた(図 3.7)。

図 3.6 Virial 学習による RMSE の変化。

34

図 3.7 Virial 学習の有無が異なる２つの手法 (ANN with Virial と ANN without Virial) に対する基本的物理量。左が部分動径分布関数。右が MSD。初期重みの異なる複数の ANN ポ
テンシャルのうち、力の RMSE が最も大きい ANN ポテンシャルによる結果をそれぞれの手法に対して示している。
Virial 学習にはもう一つ、ANN ポテンシャルの初期重みによる熱伝導度のばらつきを抑えるという効果も期待できる。なぜなら系全体のポテンシャルエネルギーによる圧力𝑷は
1
3.4.1
𝑷 = ∑ 𝑾𝑖𝛺𝑖というよう原子ごとの Virial テンソルに分解出来、この原子ごとの Virial テンソル𝑾𝑖 が2.4.8 式 (下に再掲) の熱流束𝒋𝑄 を通して熱伝導度に関与しているからである。
𝑁atom

𝑁atom

𝑁atom

𝑁type

𝑁𝜇

𝑖=1

𝑖=1

𝑖=1

𝜇=1

𝑗=1

1
1
𝒋𝑄 = ( ∑ 𝑡𝑖 𝒗𝑖 + ∑ 𝜀𝑖ANN 𝒗𝑖 + ∑ 𝑾𝑖 𝒗𝑖 ) − ∑ ℎ𝜇 ∑ 𝒗𝑗𝛺𝛺
このような繋がりから、Virial 学習によって熱伝導度が精度向上すると期待でき、実際にそのようなケースも確認されている [23]。図 3.8 は Virial 学習の有無が異なる両方で初期重みを変えた ANN ポテンシャル５つを用意し、その ANN ポテンシャルによる圧力のEP からのずれ|𝜅 ANN − 𝜅 EP |（縦軸）に対する分布
RMSE (横軸) と教師データの熱伝導度𝜅ave
ave
aveを示したものである [23]。Virial 学習をしていないものは、圧力の RMSE が大きいものが多く、熱伝導度の教師データからのずれも ANN ポテンシャルの初期重みによって大きく異なることが分かる。一方で Virial 学習をしたものは、全ての ANN ポテンシャルで圧力のRMSE が小さく、熱伝導度の教師データからのずれも小さいことが分かる。

35

図 3.8 圧力の RMSE と熱伝導度の誤差の相関 [23]しかし、今回の系では Virial 学習による熱伝導度のばらつき低減は確認できなかった。図3.9 には Virial 学習の有無が異なる２種類の ANN ポテンシャルによる熱伝導度が示してある。Total はエラーバーも含めてほとんど差がない。これは Virial 学習無しの時点で既に熱伝導度のばらつきが十分小さかったためだと考えられる。しかし、初期重みの異なる ANN
ポテンシャルなどでサンプル数を増やしていけば、Virial 学習をしていない系での熱伝導度のばらつきが大きくなるという可能性もある。

図 3.9 Virial 学習の有無による熱伝導度の変化。
以上の結果から、Virial 学習によって圧力の精度が十分に向上し、FPMD と同様の熱力学条件下で ANN-MD を行えるようになると分かった。したがって以下の議論では Virial 学
習をした ANN ポテンシャルによる結果を用いる。

3.5節

Ewald 法併用の効果

Ag2Se はイオン性が高く、クーロン相互作用を考慮する必要があると考えられる。一方
ANN ポテンシャルはカットオフ距離が設けられているため、ANN ポテンシャル単体では

36

長距離相互作用であるクーロン相互作用を考慮出来ないと予想される。そこでクーロン相互作用を計算する手法である Ewald 法と ANN ポテンシャルを併用した場合 （手法については2.3.4 項を参照） と、ANN ポテンシャルだけを用いる場合とで結果を比較した。結論から言えば動径分布関数などの基本的な物理量や熱伝導度の計算には Ewald 法は不要であった。以下ではこの結論に至る過程について説明する。
図 3.10 には計算して得られた格子熱伝導度を Ewald 併用の有無が異なるそれぞれの場合に対して示してある。見ての通り Total だけでなく各寄与においてもほとんど差は無い。
また、α-Ag2Se, 500 K の系に対しても同様のことが確かめられている [22, 23]。このことから少なくとも Ag2Se の熱伝導度計算には Ewald 法は不要だと言える。

図 3.10 Ewald 法の併用の有無による熱伝導度の変化
ANN ポテンシャルと Ewald 法を併用した場合 (ANN with Ewald と呼称する) と併用しなかった場合 (ANN without Ewald と呼称する) とで RMSE や部分動径分布関数の基本的な物理量も比較したところ、こちらも大きな違いは無かった。図には２つの場合の RMSE
(図 3.11) や動径分布関数、MSD (図 3.12) を示している。これらを見る限り Ewald 法併用の有無による違いはほとんど無い。このようになる理由として予想されるは、ANN ポテン
シャルのカットオフ半径内の原子配置にそれ以遠の原子配置に関する情報が間接的に含まれていることや、クーロン相互作用が他の作用に比べて十分小さいことなどである。

37

図 3.11 Ewald 法の併用の有無による RMSE の変化。

図 3.12 Ewald 法の有無が異なる２つの手法 (ANN with Ewald と ANN without Ewald) に対する基本的物理量。左が部分動径分布関数。右が MSD。初期重みの異なる複数の ANN ポ
テンシャルのうち力の RMSE が最も大きい ANN ポテンシャルによる結果をそれぞれの手法に対して示している。
以上のような結果から Ag2Se の動径分布関数などの基本的物理量の再現や熱伝導度計算には Ewald 法は不要だと考えられる。したがって、これ以降は基本的に ANN ポテンシャ
ル単体を使用した場合の結果を用いて議論する。

3.6節

熱流の正誤による熱伝導度の変化

熱伝導度計算に使う GK 公式には熱流束が含まれているため、正しい熱流束を用いる必38

要がある。本項では誤った熱流束 𝒋wrong 、正しい熱流束 𝒋rigorous それぞれによる熱伝導度を比較して、正しい熱流束𝒋rigorous を用いることの重要性を確認する。
GK 公式 (2.4.1 式) によると熱伝導度は熱流の自己相関係数に比例しているので、ANN ポ
テンシャルにあった正しい熱流束の表式𝒋rigorous を用いることが高精度な熱伝導度計算に必要である。ここで言う正しい熱流束というのは、ANN ポテンシャルの多体効果を考慮して
rig

導出された𝑾𝑖(詳細は2.4.3 項を参照) を用いた熱流束のことである。

一方私は研究の中盤まで誤った熱流束𝒋wrong を用いていた。誤った熱流束𝒋wrong というのは、圧力テンソルから自然に導ける𝑾wro(詳細は2.4.3 項を参照)を用いたものであり、こ𝑖の表式は 2 体ポテンシャルに対しては正しい結果を与えるが、多体効果を含んでいる ANN
ポテンシャルに対しては不適切である。
熱流束の正誤だけが異なる２つの場合の熱伝導度を図 3.13 に示している。誤った熱流束
𝒋wrong を用いた場合の熱伝導度は 0.449±0.046 Wm-1K-1 であったのに対して、正しい熱流束𝒋rigorous を用いた場合は 0.306±0.013 Wm-1K-1 になった。熱流の表式を正しくすることによって平均値が変化するだけなく、エラーバーが小さくなるという変化も見られた。これは ANN の初期重みへの依存性が減ったという点でも好ましい結果である。

図 3.13 熱流束の表式の正誤による熱伝導度の変化。シンボルとエラーバーはそれぞれ 5 つの ANN ポテンシャルによる平均と標準偏差を表している。
Ewald 法併用の場合についても同様の比較を行ったところ、やはり正しい熱流束を用いた場合の方が妥当な値であった (図 3.13)。具体的には誤った熱流束を用いた場合は 0.805±
0.269 Wm-1K-1 であるのに対して正しい熱流束を用いた場合は 0.306±0.013 Wm-1K-1 であった。このように Ewald 法併用時の方が顕著な変化を確認できた。

39

図 3.14 Ewald 法併用時の熱流束の表式の正誤による熱伝導度の変化。左が誤った熱流束、
右が正しい熱流束による結果。シンボルとエラーバーはそれぞれ 5 つの ANN ポテンシャ
ルによる平均と標準偏差を表している。
以上のような結果から、ANN ポテンシャルの多体効果を考慮した正しい熱流束を用いる
ことが適切な熱伝導度を得るために重要であると分かった。したがって以降は正しい熱流束𝒋rigorous を用いた場合の結果について議論する。

3.7節

非平衡分子動力学法による計算結果

2.4.2 項の「HNEMD 法による計算のメリット」の箇所で言及したように、非平衡分子動
力学 (HNEMD) 法の方が EMD 法よりも低い計算コストで熱伝導度を計算出来る。そこで
HNEMD 法でも EMD 法と同様な結果を得られるかを試した。HNEMD 法では系に加える摂動に対して熱伝導度が一定になる領域（線形応答領域）を探す必要がある。今回のβ相の線形応答領域はα相に比べると非常に狭かった。また、Ag2Se のβ相には異方性があり、摂動を加える方向が異なる３つの MD から得られる熱伝導度を平均する必要があった。このため、当初期待していた EMD 法から HNEMD 法へ変えることによる計算コストの低減は出来なかった。結果的に熱伝導度の値は EMD 法によるものとほぼ同じになった。以下
ではその詳細を紹介する。
HNEMD 法では系に加える摂動の大きさを調整して熱伝導度が一定になる領域（線形応答領域）を探る必要がある。摂動が小さければ EMD 法との違いが小さくなり、熱伝導度
収束までに時間がかかることになる。したがって、線形応答領域の中でできるだけ大きな摂動を用いたいところである。ただし、摂動が大きすぎれば ANN-MD が不安定になり易くなるので、実際には ANN-MD を安定して行える範囲内で出来るだけ大きな摂動を選ぶことになる。

40

摂動の強さを調整した結果𝐹ext /Ω = 0.001~0.005 bohr −1 程度が線形応答領域であると分かった (図 3.15)。これは Ag2Se のα相に対して経験的ポテンシャル [50]を用いて計算した
結果 (だいたい𝐹ext /Ω =0.003~0.06 bohr −1 が線形応答領域) [23]に比べると非常に狭い範囲であり、値自体も小さい。摂動が小さいため ANN-MD が不安定になる心配は無かった反面、熱伝導度を十分に収束させるのには 1.2 ns (500000 step) ほどのシミュレーション時間が必要であり、それ以上のシミュレーション時間でも若干のゆらぎが見られた (図 3.16)。

図 3.15 摂動の強さに対する熱伝導度の変化

41図 3.16 MDstep に対する熱伝導度の推移。それぞれ図中の摂動に対する結果を示している。
線形応答領域内だとみなせる𝑭𝐞𝐱𝐭 /𝛀 = 𝟎. 𝟎𝟎𝟑 𝐛𝐨𝐡𝐫 −𝟏においても、1000000 step まで若干の揺らぎが見られる。
β-Ag2Se の Se 副格子は直方晶で異方性があるため、方向ごとに熱伝導度の値が異なる可能性が有り、摂動を加える方向の異なる３つの MD から得られる熱伝導度を平均する必要があった。そのため必要なシミュレーション時間は 1.2 ns ×3=3.6 ns と長くなった。後述の通り、実際にこの系では熱伝導度が方向によって異なっていた。一方、Se 副格子が bcc構造であるα-Ag2Se のように対称性の高い系であれば、全方向の熱伝導度が同じだと考えられるので計算量が少なくて済む。
摂動の強さを 0.005 bohr −1として３方向の平均をとった結果、熱伝導度は 0.338±0.009
Wm−1 K −1 になり、EMD 法で得られた値 0.306±0.013 Wm−1 K −1 に近い値が得られた (図3.17)。これは EMD 法と HNEMD 法では結果がほとんど変わらないという先行研究 [51]の結果の通りである。細かいことを言えば各平均値がお互いの標準偏差の範囲外にあるが、そもそもこの物質は他の物質よりも熱伝導度が非常に小さいので、この 0.03 Wm−1 K −1 程
度の差は実質的にあまり意味が無いと考えられる。

42

図 3.17 EMD 法と HNEMD 法による格子熱伝導度の比較。右３つはぞれぞれ[100]、[010]、[001]の結晶方向の熱伝導度。左端は３方向の平均値を示す。
方向ごとの熱伝導度を見ると、[100]方向に大きいという結果が HNEMD 法でも EMD 法でも確認できた(図 3.17)。具体的には[100]方向は 0.4 Wm−1 K −1 に近く、他の方向は 0.3
Wm−1 K −1 以下だという傾向がどちらの計算手法からも得られている。このようになる原因としては、[100]方向は直方晶であるこの系の長辺方向であるため、他の方向に比べて原子の振動が邪魔されにくく熱の伝搬が起こり易いことが考えられる。
また、HNEMD 法においても熱流束の表式の正誤による熱伝導度の変化を見たところ、EMD 法の場合と同様の傾向が現れると分かった (図 3.18)。誤った表式を用いた場合には熱伝導度の平均値もエラーバー (標準偏差) も正しい表式の場合より大きくなっている。

図 3.18 HNEMD 法の場合の熱流束の表式の正誤による熱伝導度の変化
以上のように EMD 法と HNEMD 法の結果がほとんど変わらないことを確認できた。またβ-Ag2Se では結晶の異方性や線形応答領域を与える摂動の小ささや結晶の異方性などが43

原因で、EMD 法の計算コストの方が低いことも分かった。したがって以降の議論では前項までに引き続き EMD 法による結果を提示する。

3.8節

対称関数の変更による影響

ANN ポテンシャルの入力値を表現するための主要な対称関数としては Behler−Parrinello型 [25]と Chebyshev 型 [26]のものが有る (詳細は付録 Error! Reference source not found.
を参照)。以下では対称関数による熱伝導度などの変化について議論する。結論から述べると、どちらの対称関数でも物理量に差が出ないと分かった。そこで、以降の議論ではChebyshev 型を用いることにした。理由はとしてはパラメータ調整のし易さ (付録 Error!
Reference source not found.、Error! Reference source not found.を参照) が挙げられる。
計算の結果、力と圧力の RMSE に差が見られたものの (図 3.19)、基底の違いによる熱伝導度などの物理量の変化は見られなかった。動径分布関数 (図 3.20)、熱伝導度 (図 3.21)などを比較した図を以下に示している。

図 3.19 対称関数の種類による RMSE の変化。黒と赤はそれぞれ Behler−Parrinello 型とChebyshev 型に対する RMSE を表している。またそれぞれの色において、「＋」「×」はそれぞれ Train Data と Test Data に対する RMSE を示している。

44

図 3.20 対称関数の種類による部分動径分布関数の変化

図 3.21 対称関数の種類による熱伝導度の変化
RMSE は Behler−Parrinello 型の方が低かったものの、success rate (success rate = (2.4 ns
程度の MD シミュレーション中に破綻しない ANN ポテンシャルの数)/(ANN ポテンシャル数))は Chebyshev 型の方が高かった。具体的には Behler−Parrinello 型の success rate が16/28 であるのに対して、Chebyshev 型の方は 5/5 であった。統計量を揃えた場合に結果が変わる可能性は有るものの、最終の success rate が 16/28 以下になる学習条件でのデー
タ収集過程で success rate=5/5 を経る確率は非常に低いので、今回の結果を信じて良いと考えられる。RMSE の低い Behler−Parrinello 型の方の success rate が低い理由としては、過学習が考えられる。Train と Test の RMSE が変わらないので過学習が起きていないように見えるものの、教師データに入っていないデータ（FPMD で回れていない部分の位相空間に対応するデータ）に対する広義の過学習が起こっている事によって、未学習の構造に適切に対応出来ていない可能性があるということである。反対に Chebyshev 型の場合は用意されている教師データに対する学習精度が比較的低い（RMSE が高い）ことによって、広義の過学習を避けられている可能性がある。

45

以上のように Behler−Parrinello 型、Chebyshev 型どちらの基底関数を用いても物理量に差は出ず、むしろ Chebyshev 型の success rate が高くなる場合もあると分かった。

46

4章 結論本研究では Ag2Se に対して「ANN-MD を通して熱伝導度を高精度に見積もるためのANN ポテンシャル構築手法や熱伝導度計算手法の確立」を目的としていた。
手法としては、FPMD 法のデータを学習した ANN ポテンシャルによる EMD 法のデータと GK 公式を組み合わせて格子熱伝導度を見積もるという手法を用いた。
目的を達成するために、300 K の β 相を対象にして ANN ポテンシャルの構築条件や熱伝導度の計算方法を調整したところ格子熱伝導度は 0.306±0.013 Wm−1 K −1 と見積もられた。この値は実験値(0.4~0.5 Wm−1 K −1 ) [1, 2]に近い。また、次の設定が適切だと分かった。
⚫

初期重みの異なる ANN ポテンシャルを５つ以上用意する。

⚫

エネルギーに加えて力と Virial テンソルを含めた学習をする (3.4 節3.4 節)。

⚫

学習比重は次の通り (3.3 節)。 (ここでEpochとは学習回数のことである。)

𝑝𝐹
𝑝= 103 , 𝑊 = 10−3
𝑝𝐸
𝑝𝐸
𝑝𝐹
𝑝𝑊
20~40 Epoch:= 1,= 10−2𝑝𝐸
𝑝𝐸

0~20 Epoch:⚫

Ewald 法の併用は不要である (3.5 節)。

⚫

ANN ポテンシャルの多体効果を考慮した正しい熱流束を用いる (3.6 節)。

3.8.1
3.8.2

他にも以下のようなことが分かった。
⚫

HNEMD 法と EMD 法とで熱伝導度の値はそれほど変わらないものの、系に異方性があるなら計算コストの観点で HNEMD 法よりも EMD 法が効率的である(3.7 節)。

⚫

対称関数が Behler-Parrinello 型でも Chebyshev 型でも結果はほとんど変わらない(3.8 節)。

今後の展望としては本研究で分かった適切な計算手法を用いて銀カルコゲナイドの様々な組成に対して熱伝導度の温度依存性を見積もることが考えられる。

47

引用文献
[1] T. Takeuchi et al., private communication. .
[2] K. Hirata, et al., J. Electron. Mater. 49, 2895 (2020).
[3] 資源エネルギー庁 総務課戦略企画室. "平成３０年度（２０１８年度）における エネ
ル ギ ー 需 給 実 績 （ 確 報 ） ". 経 済 産 業 省 資 源 エ ネ ル ギ ー 庁 . 2020-04-30.
https://www.enecho.meti.go.jp/statistics/total_energy/pdf/stte_030.pdf, (参照 202202-17).
[4] 石田謙司. "捨てられた「廃熱」から電気エネルギーを作る：熱電変換材料". 神戸大学 . 2020-09-17. https://www.kobe-u.ac.jp/report/environmental/2020/6-3-3.html,(参照 2022-02-17).
[5] （一財）エネルギー総合工学研究所. "「IoE 社会のエネルギーシステム」 エネルギー
マネジメント研究会 産業・熱低炭素化検討グループ成果報告書（要約版）". 科学技術振興機構. 2020-05-07. https://www.jst.go.jp/sip/dl/p08/report2019_3.pdf, ( 参照
2022-02-17).
[6] T. Matsunaga, et al., MATERIALS TRANSACTIONS, 62, 16–19 (2021).
[7] P. Jood, et al., A Mater. Energy Sustain., 8, 13024–13037 (2020).
[8] Y. Ding, et al., Nat. Commun., 10, 1–7 (2019).
[9] F. Drymiotis, et al., Appl. Phys. Lett., 103, 143906 (2013).
[10] H. Wang, et al., J. Electron. Mater., 40, 624–628 (2011).
[11] F. Shimojo and M. Aniya, J. Phys. Soc. Jpn., 74, 1224–1230 (2005).
[12] A.F. Ioffe, et al., Phys. Today, 12, 42 (1959).
[13] “Theoretical Solid State Physics”, Courier Corporation, (1985).
[14] R. Franz and G. Wiedemann, Annalen Der Physik Und Chemie, 165, 497–531 (1853).
[15] S.G. Volz and G. Chen, Physical Review B, 61, 2651–2656 (2000).
[16] R. Kita and T. Dobashi, “ Nano/Micro Science and Technology in Biorheology:
Principles, Methods, and Applications”, Springer, (2015)..
[17] M.S. Green, J. Chem. Phys., 22, 398–413 (1954).
[18] W. Li et al., Comput. Phys. Commun. 185, 1747, (2014).
[19] M.E. Khaled, et al., Modell. Simul. Mater. Sci. Eng., 26, 055002 (2018).
[20] N. Artrith and A.M. Kolpak, Computational Materials Science, 110, 20–28 (2015).
[21] J. Behler and M. Parrinello, Phys. Rev. Lett., 98, 146401 (2007).

48

[22] K. Shimamura, Y. Takeshita, S. Fukushima, A. Koura, and F. Shimojo, J. Chem. Phys.,
153, 234301 (2020).
[23] K. Shimamura, Y. Takeshita, S. Fukushima, A. Koura, and F. Shimojo, Chemical Physics
Letters, 778, 138748 (2021)..
[24] 竹下雄輔、2019 年度

卒業論文、「Ag2Se の ANN ポテンシャルの 構築とその利用」(2019).
[25] J. Behler, J. Chem. Phys., 134, 074106 (2011).
[26] N. Artrith, et al., Phys. Rev. B Condens. Matter, 96, 014112 (2017).
[27] 岡崎進・吉井範行, 「コンピュータ・シミュレーションの基礎 [第 2 版]」, 化学同人,(2011).
[28] W. Kohn and L.J. Sham, Phys. Rev., 140, A1133–A1138 (1965).
[29] Y.A. LeCun, et al., Neural Networks: Tricks of the Trade: Second Edition, Springer
Berlin Heidelberg, Berlin, Heidelberg (2012), pp. 9–48.
[30] P.S. Branicio and D.J. Srolovitz, J. Comput. Phys., 228, 8467–8479 (2009).
[31] K. Levenberg, Quart. Appl. Math., 2, 164–168 (1944).
[32] D.W. Marquardt, Journal of the Society for Industrial and Applied Mathematics, 11,431–441 (1963).
[33] N. Artrith and A. Urban, Comput. Mater. Sci., 114, 135–150 (2016).
[34] R.S. Mulliken, The Journal of Chemical Physics, 23, 1833–1840 (1955).
[35] R.S. Mulliken, The Journal of Chemical Physics, 23, 1841–1846 (1955).
[36] F. Shimojo, et al., Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 77, 066103 (2008).
[37] D. Evans, Physica A: Statistical Mechanics and Its Applications, 118, 443–443 (1983).
[38] W.W. Wood, J. Stat. Phys., 57, 675–727 (1989).
[39] D. Surblys, et al., Phys Rev E, 99, 051301 (2019).
[40] P. Boone, et al., J. Chem. Theory Comput., 15, 5579–5587 (2019).
[41] Z. Fan, et al., Phys. Rev. B Condens. Matter, 92, 094301 (2015).
[42] M.R.G. Marques, et al., Phys. Chem. Chem. Phys., 21, 6506–6516 (2019).
[43] P. Hohenberg and W. Kohn, Phys. Rev., 136, B864–B871 (1964).
[44] J.P. Perdew, et al., Physical Review Letters, 78, 1396–1396 (1997).
[45] V.I. Anisimov, et al., J. Phys. Condens. Matter, 9, 767 (1997).
[46] D. Santamaría-Pérez, et al., Inorg. Chem., 51, 5289–5298 (2012).
[47] S. Fukushima, M. Misawa, A. Koura, and F. Shimojo, J. Phys. Soc. Jpn., 88, 115002

49

(2019).
[48] P.E. Blöchl, Phys. Rev. B Condens. Matter, 50, 17953–17979 (1994).
[49] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758, (1999).
[50] J.P. Rino, et al., J. Chem. Phys., 89, 7542–7555 (1988).
[51] S. Fukushima, K. Shimamura, A. Koura, and F. Shimojo, Phys. Status Solidi , 112,2000183 (2020).
[52] 上田顕, 分子シミュレーション(裳華房, 2003).

[53] M.C. Payne, et al., Rev. Mod. Phys., 64, 1045–1097 (1992).

50