動的モンテカルロ法のすゝめ

概要

粒子の運動を追跡するモンテカルロ法の一種である動的モンテカルロ法について紹介する．基礎的な観点から，動的モンテカルロ特有のアルゴリズムに起因して，移動イベントの失敗が発生しないために効率的な発展が期待できる点や，シミュレーションにおける経過時間を見積もれる点を挙げる．一方で，実践的な課題として，固体中での不純物拡散など限られた問題にしか現実的な適用が難しいという弱点もある．これらを踏まえた処方箋や発展的な課題についても紹介する．

1　 はじめに

ミクロな計算からマクロな現象の説明へ繋げることは，分子シミュレーションの研究者にとって共通の目的であろう．ただし，原子・分子のような計算から，一足飛びに流体・連続体へと繋げるといった問題は少ない．現在取り組まれていることの多くは，ミクロな系の中から重要な情報は残しつつも，少しでも大きい時空間スケールの系へ拡張していく粗視化の手法が色々と考えられている．何が重要な情報となるかは扱う現象によって千差万別であり，汎用的な手法とはなりにくく，問題に応じて様々な粗視化があり得る．

本稿で取り扱う動的モンテカルロ法（KMC）は，時間方向の粗視化を行う手法といえる．空間スケールについてはむしろ分子動力学（MD）と同程度の数nmから数百nm程度で良く，計算で扱うオブジェクトも原子・分子で良いという問題に向いている．空間スケールはミクロのままで，時間スケールだけを拡大し，長時間の現象を追いかけたいというのが，少なくとも私がKMCを利用するモチベーションである［1］．このような要請は特殊なことではなく，現在のナノ材料の形成では，空間サイズがnmスケールの物質を作るにあたって，その生成にかかる時間は数十分から数時間に及ぶことが少なくない．私の専門であるプラズマ－物質相互作用においては，例えば半導体加工などはプラズマからのイオン照射で素材をnmスケールに加工する．これに時間がかかる要因は，プラズマからイオン粒子が飛来する頻度がせいぜい${10}^{20}-{10}^{22}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}$程度と小さいことにある．100 nm²程度の表面にとってはμsに1度しかイオン粒子が到達しない．また，もっと一般的には，nmスケールの材料だからこそ，ゆっくりと加工しないとエネルギーを散逸しきれず，すぐに壊れてしまうという側面もあろう．このような現象のシミュレーションでは，原子・分子の特徴は残したまま，何とかして生成にかかる時間を追いかけたいのである．

ところで，動的モンテカルロ（kinetic Monte-Carlo）という名前はあまりに一般的すぎる．実際に1990年代に名前が定着するまでは，dynamic Monte-Carloやtime-dependent Monte-Carloという呼ばれ方もされていたらしい［2］．よって，何を最初の研究とするかは難しいが，VoterによればBeelerによる固体材料の放射線損傷時のアニーリングへの適用であるらしい［3］．そのためか，核融合炉材料の研究現場では，KMCの利用は比較的多い印象がある．

さて，KMCを用いれば，もちろん条件によりけりだが，MDでは到達できなかったようなミリ秒や秒・分といった時間に到達できる可能性がある．また，KMC特有のアルゴリズムも，シミュレーション手法開発の観点からは非常に興味深い．モンテカルロ法にも関わらず，シミュレーションの経過時間が見積れるのである．それらKMCの魅力を共有することを目指し，基本的なアルゴリズムから，実践的な課題と処方箋，実用上の問題についても紹介する．

2　 基本アルゴリズム

KMCの基本的なアルゴリズムについて紹介する．ここで分子シミュレーションにおいて代表的なモンテカルロ法であるメトロポリス法［4］と比較することで，KMCのアルゴリズムの特徴を分かりやすくする．

まずはメトロポリス法について簡単に復習する．図1(a)に概念図を示すが，粒子集団の中で1ステップに1つの粒子が移動するようなシミュレーションを考える．メトロポリス法では，基本的に次のような手順で1ステップの更新を行う．

図1  (a)メトロポリス法の概念図．(b)格子系における動的モンテカルロ法の概念図．灰色の球は移動可能な粒子を，太い矢印は移動イベント候補を示す．

1．乱数を生成して移動対象となる粒子$i$を選ぶ．

2．乱数を生成して移動先となる候補地を選ぶ．

3．現在地から候補地への移動が成功する確率$p_i$を求める．一般的なMDのポテンシャルモデル関数などから粒子配置に応じたエネルギーを求め，詳細釣り合いの原理を満たすように確率を決める．よく用いられる方法としては，移動前の状態のエネルギーを$U_i$，移動後の状態のエネルギーを$U_i^{\prime }$としたとき，$p_i=\min [1,\exp \{-\beta (U_i^{\prime }-U_i\})]$とする．

4．区間$[0,1)$の乱数$\xi$を生成し，$\xi \le p_i$であればイベント成功と見做して粒子を移動させる．そうでなければ，イベント失敗として粒子は移動させない．

5．次のステップへ進むため，手順1へ戻る．

メトロポリス法の強みはその汎用性の高さにある．粒子配置からエネルギー$U_i$や$U_i^{\prime }$が評価できれば，この手順を繰り返すことで粒子集団の運動を追いかけることができる．ただし，問題となるのは，手順4でイベントが失敗し得ることである．エネルギー差が大きかったり，温度$k_{\mathrm{B}}T=1/\beta$が低かったりすると成功確率$p_i$が小さくなってしまい，粒子集団の運動が殆ど進まなくなることが現実的にあり得る．

イベント発生の効率という点に注目すると，KMCのアルゴリズムには原理的にイベントの失敗が起こらないという強みがある．すなわち，1ステップに1つの粒子が必ず移動する．KMCでも粒子集団の運動を考えることができるが，図1(b)のように，各々の粒子は決まった格子点上だけを動くものとしよう．簡単のため，1度に動ける距離は隣接点までとする．このような条件の下で，KMCのアルゴリズムは次の手順となる．

1．最初に全ての粒子に対して移動先とその発生確率$p_k$をリストアップしておく．

2．全イベントの発生確率の和を求める．

3．乱数を生成して，イベントを選択する．選ばれたイベントは成功と見做し，粒子を移動させる．

4．乱数を生成してこのステップでの時間刻み$\tau$を求める．

5．移動した粒子，および，周辺の粒子について，移動イベントの発生確率のリストを更新する．

6．次のステップに進むため，手順2へ戻る．

それぞれの手順の詳細は以後に述べるが，まずは手順3においてイベントを必ず成功させているところに注目してほしい．発生確率$p_k$の大小に関わらず，必ず何か1つの粒子が移動することになる．これによって効率的な発展が期待できる．

それでは，各手順における詳細について説明する．

手順1は最初に一度だけ行う処理であるが，全ての粒子に対して，全ての移動候補をリストアップする必要がある．図1(b)を例にとると，サイトf，c，nに3つの粒子が存在する．ここで1番の粒子の移動先の候補は隣接するサイトb，e，g，jである．同様に，2番の粒子の移動先はサイトb，d，gであり，3番の粒子の移動先はサイトk，m，oである．より一般的に，粒子$i$の今いるサイトを$s_i$とし，そこからの移動先（すなわち$s_i$の隣接サイト）を$n\in N(s_i)$とする．$N(s_i)$は全ての隣接サイトであり，そのうちのある一つを$n$と書くという意味である．このとき，粒子が移動イベントの発生確率を$p(i,n)$と書くことにする．ただし，後の議論を見やすくするため，粒子と隣接サイトのインデックスを纏めて$k=(i,n)$とし，$p_k=p(i,n)$と略記するので注意されたし．

発生確率も手順1で全て求めておく必要がある．最も単純な方法は，メトロポリス法と同様に，移動前後のエネルギー$U(s_i)$, $U(n)$の差から

  

$$\begin{array}{c}\hfill p(i,n)=p_0\min [1,\exp \{-\beta (U(n)-U(s_i))\}]\end{array}$$(1)

とすることである．$p_0$はプレファクターである．

手順2では，全てのイベントとその発生確率が列挙された状態から，発生確率の総和として全確率$P$を求める．

  

$$\begin{array}{c}\hfill P=\sum _i^{}\sum _n^{N\left(s_i\right)}p(i,n)=\sum _k^{}{p}_k\end{array}$$(2)

ここで第二式の1つ目の和の記号は全ての粒子についての和，2つ目の和の記号は粒子$i$のいるサイト$s_i$から見た隣接サイト$n$についての和である．

手順3では一様乱数$\eta$を範囲$[0,P)$から発生させる．範囲の上限を式(2)で求めた全確率とすることが重要である．そして，次の関係を満たすイベント$k^{\prime }=(i^{\prime },n^{\prime })$を選択する．

  

$$\begin{array}{c}\hfill \sum _{k=1}^{k^{\prime }-1}{p}_k\le \eta <\sum _{k=1}^{k^{\prime }}{p}_k\end{array}$$(3)

すなわち，イベントの発生確率を順に積算していき，その和が乱数$\eta$を初めて超えたところが$k^{\prime }$番目のイベントであるという意味である．そして，選ばれた$k^{\prime }=(i^{\prime },n^{\prime })$に相当する粒子$i^{\prime }$をサイト$n^{\prime }$へ実際に移動させる．

手順4ではこのステップのイベントの時間刻み$\tau$を見積もる．実はここまで発生確率と呼んできた$p_k$は，KMCにおいては時間の逆数の次元をもっており，正確には発生頻度（英語ではratioではなくrate）と呼ぶ方が相応しい．詳細は後の節で記述するが，範囲$[0,1)$から生成した一様乱数$\xi$を用いて，このステップの時間刻みは次のように与えられる［5］．

  

$$\begin{array}{c}\hfill \tau =-\frac{\ln (1-\xi )}{P}\end{array}$$(4)

時間刻みは毎ステップ異なる値となる．ともかく，これを積算していけばシミュレーション全体での経過時間が得られる．

手順5では，次のステップに進むためにイベントリストの更新を行う．手順1で列挙したイベントのうち，移動させた粒子$i^{\prime }$については現在位置が$s_{i^{\prime }}\to n^{\prime }$と変更されたために，サイト$s_{i^{\prime }}$を起点としたイベントは全て無効となりリストから外す．新たにサイト$n^{\prime }$を起点とした移動先を列挙し，それらについての移動確率$p(i^{\prime },m)$，ただし，$m\in N(n^{\prime })$をリストに加えることになる．また，これ以外に，移動しなかった周辺粒子に対してもリストを更新する必要がある．例えば，粒子1がサイトfからgに移動したとすると，粒子2にとっても移動先の候補だったサイトgが埋まってしまうことになる．よって移動確率$p(2,g)=0$と更新してリストから外す必要がある．モデルを詳細にした場合には，周辺粒子の配置によってエネルギー$U_i$が変わるので，それらに合わせて発生確率も$p_k\to {p^{\prime }}_k$へと再計算する必要がある．これで1ステップの処理はすべて終了し，手順2へ戻る．一部の発生確率が手順5で行使されたので，式(2)に従って全確率$P$も更新する必要がある．

手順5の処理は一見すると入り組んでいて面倒だが，要は粒子配置の変化を発生確率に反映させたいのである．よって，手順1に戻って全ての洗い出しを行っても良いが，実際上は手順1の処理がKMCでは最も計算量の大きい処理である．そこで手順5では，イベントの更新対象を必要最低限に留めた結果，移動粒子とその周辺だけに限ることができ，無駄な計算が増えるのを抑えている．

ともかく，手順1および手順5によって，常に全てのイベント候補が列挙された状態にあり，全確率$P$が分かっている．このことが，手順3において選んだ確率を必ず発生させられるというKMCの利点に繋がっている．

3　 経過時間の見積もり

さて，ここまででKMCの1つ目の特徴であるイベントの失敗がない点について説明した．ここではKMCの2つ目の特徴である経過時間を見積もれる点について説明する．

KMCは広い意味でのマルコフ過程に従っており，経過時間をも積もる際にはポアソン過程に倣った議論が行われる［6］．先に導入した全確率$P$は，あらためて単位時間当たりのイベント発生回数である．ここで定常だと仮定すると時刻$t$から$t+dt$の間に少なくとも1つのイベントが起こる確率は

  

$$\begin{array}{c}\hfill Pdt+O(dt)\end{array}$$(5)

で与えられる．このとき，$O(dt)$は時刻$t$から$t+dt$の間に2つ以上のイベントが起こる確率であり，これは$dt\to 0$において第一項よりも早く0に収束するとする．$P$はポアソン過程の用語で言えば推移確率に相当する．

あるイベントの発生直後，基準時刻$t=0$に系がある状態$\mathrm{\Omega }$に遷移したとする（ここまでの例では状態$\mathrm{\Omega }$は粒子配置）．そこから時間$t$だけ経過した時に，まだ次のイベントが発生しておらず，系が状態$\mathrm{\Omega }$のままである確率を$R(t)$とする．微小時間$dt$の間のイベントの発生確率は$Pdt+O(dt)$であったので，時間$t+dt$になっても次のイベントが発生していない確率$R(t+dt)$は

  

$$\begin{array}{c}\hfill R(t+dt)=(1-Pdt-O(dt))R(t)\end{array}$$(6)

である．これより，次の微分方程式が得られる．

  

$$\begin{array}{c}\hfill \frac{dR(t)}{dt}=-PR(t)\end{array}$$(7)

この解は，

  

$$\begin{array}{c}\hfill R(t)=\exp (-Pt)\end{array}$$(8)

である．ただし，初期値として$R(0)=1$とした．数式の形から$t\to \infty$で$R(t)\to 0$であることもわかり，直感に反しない．

ここで，1ステップあたりの時間刻み$\tau$を考えるために，発生時間の密度分布関数$f(\tau )$を導入する．すなわち，経過時間が$\tau$から$\tau +d\tau$の間にちょうど次のイベントが起こる確率を$f(\tau )d\tau$とする．系が状態$\mathrm{\Omega }$のまま留まっている確率$R(t)$とは次の関係で結ばれる．

  

$$\begin{array}{c}\hfill R(t)=1-{\int }_0^{t}f(\tau )d\tau =-{\int }_{\infty }^{t}f(\tau )d\tau \end{array}$$(9)

ここに$R(t)$の解である式(8)を代入して両辺を微分すれば，発生時間に関する密度分布関数$f(\tau )$が次のように得られる．

  

$$\begin{array}{c}\hfill f(\tau )=P\exp (-P\tau )\end{array}$$(10)

これは，ポアソン過程における待ち行列分布の密度関数に相当する．すなわち$f(\tau )$は，あるイベントが起こってから次のイベントが起こるまでに掛かる時間$\tau$の分布である．

ここで念のため，イベントの全確率$P$と発生確率密度分布関数$f(\tau )$の関係を整理しておく．状態$\mathrm{\Omega }$にあった系が時刻$t+dt$の間に別の状態へ遷移する確率は，式(6)のように書くと$(P+O(dt))R(t)dt$である．一方で，発生時間の密度分布関数を使って書けば，$f(t)dt$である．$dt\to 0$の極限で両者が一致することは，式(8)および式(10)からわかる．

次のイベントが起こるまでに必要な時間$\tau$の密度分布が分かったので，期待値$\mathcal{⟨}\tau \mathcal{⟩}$を計算することができる．

  

$$\begin{array}{ccc}\mathcal{⟨}\tau \mathcal{⟩}& =\hfill & {\int }_0^{\infty }\tau f(\tau )d\tau \hfill \\ & =\hfill & {\int }_0^{\infty }\tau P\exp (-P\tau )d\tau \hfill \\ & =\hfill & \frac{1}{P}\hfill \end{array}$$(11)

これは最初の乱数を使った定義の期待値と一致する．

続いて，式(4)のように乱数を用いて時間刻み$\tau$を決める方法を考える．一般論として，分布関数$f(x)$に従う乱数列$x$は，一様乱数$\xi \subset [0,1)$から生成することが可能である．変数$u$が$[0,1)$の範囲で一様分布$g(u)=1$に従うものとする．ここで変数が$u$から$u+du$の間の値を取る確率$g(u)du=du$と，変数$x$が$x$から$x+dx$の間の値を取る確率$f(x)dx$が等しくなるような変換$x=\varphi (u)$を考える．すると，逆関数${\varphi }^{-1}$を用いて，

  

$$\begin{array}{c}\hfill du=\frac{d{\varphi }^{-1}(x)}{dx}dx=f(x)dx\end{array}$$(12)

である．これより，分布関数$f(x)$の不定積分を$F(x)={\int }_0^{x}f(x^{\prime })d{x}^{\prime }$とすれば，

  

$$\begin{array}{c}\hfill \varphi (u)=F^{-1}(u)\end{array}$$(13)

である．すなわち，変数$u$の代わりに一様乱数$\xi \subset [0,1)$の列を生成し，目的の分布関数$f(x)$の不定積分の逆関数を通して得た値$x=F^{-1}(\xi )$の列は，分布$f(x)$を成す乱数列となる．

これを応用すると，時間刻み$\tau$は式(10)の分布関数$f(\tau )$に従う$\tau \subset [0,\infty )$の範囲の乱数列と見做すことができる．すると，

  

$$\begin{array}{c}\hfill F(\tau )={\int }_0^{\tau }f({\tau }^{\prime })d{\tau }^{\prime }={[-\exp (-P{\tau }^{\prime })]}_0^{\tau }=1-\exp (-P\tau )\end{array}$$(14)

である．一様乱数$\xi \subset [0,1)$に対して$\xi =F(\tau )$であることから，逆関数を考えると

  

$$\begin{array}{c}\hfill \tau =F^{-1}(\xi )=-\frac{\ln (1-\xi )}{P}\end{array}$$(15)

となって式(4)に一致する．念のため，対数部分の発散は起こらないことは，一様乱数の範囲から明らかである．

このように，ポアソン過程に従うと仮定することで，KMCの経過時間を見積もることができる．ただし，ここまで1ステップの時間刻みに対してだけの議論を行った．複数ステップ経過した時の経過時間は，本当に各ステップで求めた$\tau$の積算で良いのかという点については，詳細は割愛するが，アーラン分布を考えてやればよい．全確率$P$が殆ど変化しないと仮定できる場合には，$s$ステップ後の経過時間の期待値$\mathcal{⟨}{t}_s\mathcal{⟩}$は，

  

$$\begin{array}{c}\hfill \mathcal{⟨}{t}_s\mathcal{⟩}=s\mathcal{⟨}\tau \mathcal{⟩}=\frac{s}{P}\end{array}$$(16)

となる．これより，複数ステップでの経過時間は，1ステップの時間刻みの積算で良さそうである．

4　 実践的な課題

ここまでで，KMCのアルゴリズムの基本的な部分について述べた．ここからは実際のシミュレーションを行うための実践的な課題を挙げる．

4.1　 適用可能な系

先にも述べたが，KMCのアルゴリズムの最大の特徴であると同時に最大の問題は，粒子の移動先を全て列挙できなけらばならないという点である．つまり，移動対象と移動先が有限でなければならない．この制限から，粒子の移動先を無限に考えられる気体や液体のような系には適さない．

歴史的にも適用例の多いのは，固体結晶中の不純物拡散であったり空孔拡散が多い．これらは，不純物原子の位置できる場所を格子間サイト（安定サイト）に限定したり，空孔の位置できる場所を原子サイトに限定することで，移動イベントの数を有限に抑え込むことができる．現実的な固体結晶では，粒子当たりの移動先はせいぜい10以下であり，粒子数$N$に対して$O(N)$の列挙で済む．

同様の適用先として，固体結晶表面のアドアトム拡散にも応用される．

4.2　 エネルギー差か障壁エネルギーか

次の問題は，粒子の移動イベントに対して，移動確率$p_k$をどのように定めるかということである．先の例ではメトロポリス法と同様に移動前後のエネルギー差から確率を設定したが，そのエネルギーをどのように見積もるかが，KMCのモデル化において重要な差異を生む．

規則的な格子上での運動を条件にした最もシンプルなモデルは，粒子$i$ごとにの隣接サイトに存在する粒子数（結合数）$b_i$を数えて，エネルギーを$b_i$の関数として$U(b_i)=c_1{b}_i+c_2\sqrt{b_i}$のように定めてしまうことである．ここで係数$c_1$，$c_2$は実験値や分子動力学（MD）のポテンシャル，密度汎関数理論（DFT）計算などと合うように決めてやる．結合数以外の変数を加えて，もう少し複雑な関数でモデル化することもある．

また，移動確率$p_k$はエネルギー差以外から決めてもよい．詳細釣り合いの原理を満たすように設定すればよいので，より詳細な物理過程を反映するには障壁エネルギー$\mathrm{\Delta }{E}_k$によって，$p_k=p_0\exp (-\beta \mathrm{\Delta }{E}_k)$とすることができる．特に昨今では，DFTが普及したことから，Nudged Elastic Band（NEB）法［7］などによって，固体中の不純物粒子移動のエネルギー障壁を精度良く見積もることが容易になった．格子が規則的であることから原子配置とエネルギー障壁の間の関係をモデル化する．関数としてモデル化でき無い場合には，テーブル化・データベース化することになる．

一方で，ここのエネルギー算出に掛かる計算時間はKMCの計算の大部分を占めることになる．そのため，毎ステップMDのポテンシャルから算出するという方法を選択するには注意が必要である．

4.3　 プレファクターの決定

さて，経過時間を見積もることができる根拠は，イベントの発生確率$p_k$は，実際には単位時間当たりのイベント発生回数であるためであった．そのための時間の逆数の次元は，プレファクター$p_0$に押し付けられている．そのため，プレファクターにどのような値を設定するかが，KMCの再現時間を決めるカギとなる．大別すると，実験から決める方法と，遷移状態理論を用いて第一原理的に決める方法の2つが挙げられる．

実験から決める方法としては，粒子の拡散係数から算出することができる．確率$p_k$で発生する移動イベントにおいて，移動距離を$l_k$，移動障壁エネルギーを$\mathrm{\Delta }{E}_k$とすると，単位時間当たりの移動距離の二乗平均は

  

$$\begin{array}{c}\hfill \frac{\mathcal{⟨}{l}^2\mathcal{⟩}}{\mathcal{⟨}\tau \mathcal{⟩}}=\sum _k^{}(l_k{)}^2{p}_0\exp (-\beta \mathrm{\Delta }{E}_k)\end{array}$$(17)

となる．一方で，拡散係数$D(\beta )$との間にアインシュタインの関係式

  

$$\begin{array}{c}\hfill \mathcal{⟨}{l}^2\mathcal{⟩}=2dD(\beta )\mathcal{⟨}\tau \mathcal{⟩}\end{array}$$(18)

を課す．ただし，$d$は空間の次元である．これらの比較から，$l_k$と$\mathrm{\Delta }{E}_k$の具体的な値が分かれば$p_0$が求まる．実践的には，基本的にプレファクター$p_0$は系の状態や場所に依らず，一度定めたら常に共通の定数である．そこで，系の中で最もシンプルな状況で設定してやればよい．より単純化できる状況として，1粒子だけ存在するものとし，規則正しい格子上の移動とすれば$l_k$および$\mathrm{\Delta }{E}_k$も共通の値$l$，$\mathrm{\Delta }E$とできるだろう．各サイトの隣接数を$n$とすれば，$p_0=2dD(\beta )/n{l}^2\exp (-\beta \mathrm{\Delta }E)$と簡単に見積もることができる．具体的な例として，タングステンの(110)表面のアドアトム拡散を取り上げると，$l$は格子定数となり，隣接数は$n=4$，次元は$d=2$となる．アドアトム拡散の実験値は，例えば$D(\beta )=D_0\exp (-\beta \mathrm{\Delta }{E}^{\prime })$という形式で報告されており，$D_0=2.6\times {10}^{-7}[{\mathrm{m}}^2\mathrm{/}\mathrm{s}]$，$\mathrm{\Delta }{E}^{\prime }=0.92[\mathrm{e}\mathrm{V}]$である［8］．

注意点として，原子レベルの実験値として複数の報告がある場合，それらは桁で違ってくることが多い．また，表面拡散からプレファクターを決めた場合と，内部の空孔拡散から決めた場合とでは，異なる値になってしまうこともしばしばある．この辺りは見たい現象に合わせて塩梅を決める必要も出てくる．

また，実験から報告される拡散係数は，アレニウスプロットされた温度の関数として報告されるため，プレファクターも温度の関数になる．それを嫌って，KMCのモデル中で採用する$\mathrm{\Delta }E$と，実験値として報告される$\mathrm{\Delta }{E}^{\prime }$を一致させてやればよいかと言えば，必ずしもそう単純ではない．実験で測定される量は，多体効果や自由エネルギーの効果を含んだエフェクティブな量であるため，ミクロからモデル化したKMCの設定値と一致しないのは当然という考え方もできる．

より第一原理的な観点でプレファクターを決める方法として，遷移状態理論（TST）を用いる方法がある．基本的な考え方として，KMCで粒子が位置するサイトは，ミクロにはポテンシャル極小点に対応し，粒子は極小点付近のポテンシャル場の上で振動している．近似のレベルに応じていくつかの手法が提案されているが［9–14］，この振動の振動数がプレファクターに対応する．

4.4　 計算上のボトルネックと二分木探索

今度は，コード開発の観点から，実践的な処方箋を紹介する．

先に述べたKMCの手順1から6のアルゴリズムのうち，シミュレーションにおける計算時間のボトルネックになるのは，以外にも確率の総和を取る手順2とイベントを選択する手順3である．他の部分が非常に高速に処理できるために，式(3)に従ってイベントの発生確率の足し算と比較処理を行う部分に相対的に時間が掛かる．

手順1は初回に一度だけ行えばよいので，手順2から手順5の計算量を簡単に見積もってみる．各々の粒子の移動先はせいぜい数か所であるので，粒子数$N$に対してイベント数は比例する．手順2では式(2)の総和を取るために単純な逐次処理を行うと，足し算の回数はイベント数に比例した$O(N)$となる．さらに手順3の式(3)の処理でも，逐次的に$p_k$を積算しつつ乱数$\eta$との比較を行う．この足し算と比較の回数も$O(N)$となる．手順4は乱数から時間を見積もるだけであり$O(1)$である．手順5では，エネルギー$U_i$や移動障壁エネルギー$\mathrm{\Delta }{E}_k$を算出しなおす必要がある．この計算は，モデル関数やテーブル参照となり，計算量は実装次第ではあるものの，算出しなおす対象となる粒子は，移動した粒子とその周辺の粒子だけであってせいぜい数個である．よってエネルギー算出の計算量は$O(1)$と言える．

このように，粒子数が増えるほど，式(2)や式(3)といった単純な計算が大きな負荷となる．これを解決する大変有効な方法として，足し算や比較の回数を$O({\log }_{2}N)$で抑える方法がBlue［15］らによって提案されている．その方法をさらに効率よくするため，二分木探索を応用した方法をここで紹介する［16］．

図2のような二分木を考える．各ノードはKMCの1つのイベントを指す．さらに，各ノードは左右に計2つの子ノードを持つ．このノードの親子関係は，あくまでデータ構造の話であり，イベント自体に優劣はない．また，ステップが進むにつれてイベントの追加や削除が起こると親子関係は入れ替わることがある．左右片方だけのノードを持つこともあるし，末端となる子ノードはそれ以上子ノードを持たない．さて，各ノードは情報としてイベントを指すインデックス$k$と，そのイベントの発生確率$p_k$を保持する．加えて，自身を含めて子ノード以下にぶら下がる全ての発生確率の和$q_k$を持つ．あるノード$k$の左右の子ノードをそれぞれ$k_L,k_R$としたとき，$q_k$は次の関係を満たす．

図2  二分木を用いたKMCのイベント選択の概念図．この場合ノード1がルートであり，直接の子はノード2およびノード8である．各々のノードは発生確率$p_k$と自分以下の子ノード全ての発生確率の和$q_k$を保持する．矢印はイベント選択時の探索経路であり，発生確率の更新時はその逆の経路を辿る．

  

$$\begin{array}{c}\hfill q_k=p_k+q_{k_L}+q_{k_R}\end{array}$$(19)

ただし，左右のノードを持たないときは，$q_{k_L}$や$q_{k_R}$を0とみなす．

このような二分木構造を構築すると，最上位の親（ルート）$k=r$の持つ総和$q_r$は，KMCの全確率$P$に等しい．つまり二分木の構築によって自然と，式(2)の全確率の計算が完了する．手順5において，あるノード$k^{\prime }$のイベント確率$p_{k^{\prime }}$が更新されたときは，すぐに$q_{k^{\prime }}$も更新する．このとき子ノード以下の再計算は不要である．そして$q_{k^{\prime }}$の更新を親ノード$k^{″}$に伝える．更新の伝えられた親ノード$k^{″}$は自身の$q_{k^{″}}$を更新して，さらに親に伝える．これをルート$r$に到着するまで繰り返せば，ルートの持つ全確率$q_r$が更新される．ルートに到着するまでの処理回数は$O({\log }_{2}N)$となっている．イベントのノードを追加したり削除したりした時も，同様に$q_k$を更新しながら親に向かって遡るだけで良く，処理回数は$O({\log }_{2}N)$ですむ．

続いて，この二分木を使って手順3のイベント選択を行う．乱数$\eta$が与えられたとき，まず最初にルートのノード$k=r$を起点とする．ここで，(i)そのノードの発生確率$p_k$と比較し，$\eta <p_k$ならばイベント$k$を選択して完了する．(ii)これを満たさない場合，$\eta <p_k+q_{k_L}$ならば，$\eta \leftarrow \eta -p_k$と更新してから左の子ノード$k_L$へ移り，(i)へ戻る．(iii)これも満たさない場合は，$\eta \leftarrow \eta -p_k-q_{k_L}$と更新してから右の子ノード$k_R$へ移り，(i)へ戻る．これを繰り返すことで，やがて条件を満たすイベントが見つかる．この時の比較回数はせいぜい$O({\log }_{2}N)$となっている．

こうして，単純な逐次処理では$O(N)$となった計算回数を$O({\log }_{2}N)$に削減することができる．ただし，これは理想値である．実用的には，二分木構造がどのノードから見ても左右に連なる子ノードの個数がバランスが取れている必要がある．もしバランスが取れておらず，例えばルートノード以下が全て右の子ノードしか持っていない一次元鎖状になったとすると，全ての処理は最悪の$O(N)$となる．左右のバランスをとるために，ノードの追加や削除に応じて，定期的に親子関係の入れ替えを行うなどの対策が講じられる．我々のコードでは二分木構造のアルゴリズムのうちAVL木［17,18］を用いた実装を行っている．

この二分木探索の効果は非常に大きく，エネルギー計算のモデルが簡単であれば現在の計算機において1コアで${10}^{11}$ steps/dayの計算が可能である．

4.5　 MDと比べて本当に速いか

さて，KMCの1ステップ当たりの計算速度は$O({\log }_{2}N)$と非常に速くなった．比較対象としてMDを取り上げるならば，固体のような短距離相互作用だけの計算であっても，1ステップ当たりの計算量は$O(N)$である．これだけ見るとKMCは速いように見える．しかしながら，状況によってはMDを下回ることがあることも注意しておかなければならない．

今行った比較は実はフェアではない．KMCの場合，1ステップ当たりの経過時間は粒子数や温度に応じて小さくなる．式(11)にあるように1ステップ当たりの時間刻み$\mathcal{⟨}\tau \mathcal{⟩}$は全確率$P$に反比例するが，ここで全確率$P$は粒子数$N$に比例する．つまり，粒子数が増加すると，目標時間$t_{\mathrm{tgt}}$に達するのに必要なステップ数が粒子数に比例して増加することになる（$t_{\mathrm{tgt}}/\mathcal{⟨}\tau \mathcal{⟩}\sim N$）．すると，$t_{\mathrm{tgt}}$に達するまでに必要な計算量は$O(N{\log }_{2}N)$となってしまう．

また，温度が上がると，確率$p_k$は指数関数的に1に近づくが，逆に時間刻み$\mathcal{⟨}\tau \mathcal{⟩}$は指数関数的に減少することになる．よって高温になると時間が進みにくくなる．

一方でMDの場合には，タイムステップは粒子数や温度に依存しない．KMCの計算量$O(N{\log }_{2}N)$とMDの計算量$O(N)$を比べれば，粒子数を増やしていくとKMCとMDで計算時間が逆転することになる．また，MDの場合には領域分割法により粒子数の増加をカバーできるが，KMCの場合はアルゴリズムが逐次的ゆえに並列化が難しいという欠点もある．もちろん，KMCならば，トラップサイト内での粒子の長時間の振動を無視して，1ステップで別のサイトへ移動できるので，これだけでMDと比べて時間的にかなり有利ではあるが，粒子数が多い場合や高温の場合にはMDと比べて本当に早いか検討してから利用した方が良い．全てKMCにしようとせず，MDや他の手法とのハイブリッド化［19,20］を検討することも良い選択肢である．

4.6　 On-the-fly KMC

さて，KMCのアルゴリズムでは系の中の全イベントを事前に列挙出来る必要があり，そのために規則的な格子状のモデルが好まれる．固体結晶で構造が大きく変わらない現象を扱っているときはこれでも良い．しかし，固体であっても結晶粒界やアモルファス構造の場合には，規則的な格子を組むことができない．ましてや，不規則な格子において，サイト間のエネルギー差や障壁エネルギーをシンプルにモデル化するのは困難である．そこで，不規則な格子に対してもKMCを適用するための幾つかの方法が提案されている．

我々が開発した手法もその一つで，不規則な格子中の不純物拡散現象を対象として，手順1での移動サイトの探索を自動的に行い，合わせて障壁エネルギーも算出するものである［21］．移動サイトを探すには，系の中から半径1～数nm程度の局所的な領域を切り出してきて，不純物原子を配置し，MDによる構造緩和を行ってやればよい．続いて，見つかったサイトの中から距離の近い2つを選び，再び局所的に切り出した領域でのMDでNEB法を行えば，サイト間の移動障壁エネルギーも算出できる．この計算を系全体に対して網羅的に行えば，原理的には全てのイベント見つかる．ここで，各々のMD計算は場所ごとに独立に行うことができるので，マルチプロセス環境での並列化（MPMD［22］）が可能である．

さらに進んで，手順5においても粒子の移動に伴うイベントの再計算をMDで行えばよい．このとき，移動させる粒子を選択された粒子に限らず，それに伴って周辺の粒子の構造緩和もMD行ってやれば，母材の構造が変わっていく様子も扱えるようになる．そのような手法をOn-the-fly KMCと呼ぶ．代表的なコードとしては，kART［23］やSEAKMC［24］が存在する．

我々も上述の方法を拡張してFlyAMというコードを開発しているが，それゆえに課題も見えてきた．On-the-fly KMCは汎用的なフレームワークではあるものの，MDのポテンシャルモデルに強く依存する．そして，固体系のMDのポテンシャルモデルは十分な精度で汎用的なものがあるとは言い難い状況である．ただし，最近流行りの機械学習ポテンシャルがより使いやすくなれば（ユーザーが自前のコードと結合しやすくなれば），On-the-fly KMCが真の意味で活躍する日が来るであろう．

5　 まとめ

本稿ではKMCの基本的なアルゴリズムを紹介するところから始まり，実践的な高速化の方法や，実用上の課題を紹介した．分子シミュレーション分野のメトロポリス法から始まる多くのモンテカルロ法と比べると，KMCではアルゴリズムの違いに起因する異なる課題に出会えて興味深い．開発者や利用者も相対的に多くないために，まだまだ新しいKMCの発展や応用があり得るだろう．ぜひ一度，現在の研究対象にKMCを利用できないか検討してみて欲しい．より詳しい解説としてVoterによるレビュー［2］が大変参考になる．また，上記で紹介したAVL木を用いたKMCのサンプルコードをGithub［16］に公開しているので，興味のある方はぜひご覧いただきたい．

 謝辞

筆者がKMCの研究を立ち上げるにあたり，核融合研（当時）のCOE研究員であった小田泰丈博士に基本的なアルゴリズムから応用例まで，多くのことをご教授いただきました．また，MDとのハイブリッド化やプレファクターの設定などでは東北大の村島隆浩博士のご助言に助けられました．この場を借りてお礼申し上げます．本研究はJSPS科研費JP23K17679，24K00617，JP24H02251，JP24H02246の助成を受けたものです．

参考文献

• [1]   A. M.  Ito, JSAP Review, 2023, 230412 (2023).
• [2]   A. F.  Voter, “Introduction to the Kinetic Monte Carlo Method”, Radiation Effects in Solid, Mathematics, Physics and Chemistry Vol. 235, edited by K. E. Sickafus, E. A. Kotomin, and B P. Uberuaga, Springer, pp. 1–24 (2004).
• [3]   J. R.  Beeler, Jr., Phys. Rev., 150, 470 (1966).
• [4]   N.  Metropolis,  A. W.  Rosenbluth,  M. N.  Rosenbluth,  A. H.  Teller and  E.  Teller, J. Chem. Phys., 21, 1087–1092 (1953).
• [5]   A. B.  Bortz,  M. H.  Kalos and  J. L.  Lebowitz, J. Compt. Phys., 17, 10–18 (1975).
• [6]   K. A.  Fichthorn and  W. H.  Weinberg, J. Chem. Phys., 95, 1090 (1991).
• [7]   G.  Mills and  H.  Jónsson, Phys. Rev. Lett., 72, 1124 (1994).
• [8]   W.  Xu and  J. B.  Adams, Surf. Sci., 319, 58–67 (1994).
• [9]   H.  Mehrer, Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion-Controlled Processes, Springer, Berlin, Chap. 4, pp. 55–67 (2007).
• [10]  H. メーラー著，藤川辰一郎訳，固体中の拡散：基礎と方法，異種物質中の拡散，拡散律速過程，Springer，東京 (2007).
• [11]   H.  Eyring, J. Chem. Phys., 3, 107 (1935).
• [12]   H.  Eyring, Trans. Faraday Soc., 34, 41 (1938).
• [13]   E.  Wigner, Z. Phys. Chem. Abt. B, 19, 203 (1932).
• [14]   E.  Wigner, Trans. Faraday Soc., 34, 29 (1938).
• [15]   J. L.  Blue,  I.  Beichel and  F.  Sullivan, Phys. Rev. E, 51, R867 (1995).
• [16]  https://github.com/atsushi-m-ito/
• [17]   G.  Adelson-Velsky and  E.  Landis, “An algorithm for the organization of information”, Proceedings of the USSR Academy of Sciences, 146, 263–266 (1962); English translation by Myron J. Ricci in Soviet Mathematics - Doklady, 3, 1259–1263 (1962).
• [18]   R.  Sedgewick, “Balanced Trees”, Algorithms, Addison-Wesley, p. 199 (1983).
• [19]   A. M.  Ito,  A.  Takayama and  H.  Nakamura, Materials Research Express, 10, 125002 (2023).
• [20]   A. M.  Ito,  A.  Takayama and  H.  Nakamura, Plasma and Fusion Research, 13, 3403061 (2018).
• [21]   A. M.  Ito,  S.  Kato,  A.  Takayama and  H.  Nakamura, Nucl. Mater. & Energy, 12, 353–360 (2017).
• [22]   A.  Takayama,  K.  Shimizu,  Y.  Tomita and  T.  Takizuka, J. Plasma Fusion Res. SERIES, 9, 604–609 (2010).
• [23]   N.  Mousseau and  G. T.  Barkema, Phys. Rev. E, 57, 2419 (1998).
• [24]   H.  Xu,  Y. N.  Osetsky and  R. E.  Stoller, Phys. Rev. B, 84, 132103 (2011).

著者紹介

伊藤 篤史（博士（理学））

〔経歴〕2009年名古屋大学大学院理学研究科博士課程修了，2010年核融合科学研究所に入所（現所属）．〔専門〕プラズマ－物質相互作用．〔趣味〕ガンダム，木工，F1鑑賞．