レッドムーン法：一つの“モンテカルロ”法，何なら元祖進化形！？

概要

複雑な反応現象を伴う分子凝集系である複合反応（complex reaction, CR）系を扱うためのレッドムーン（Red Moon, RM）法を紹介する．まず，スタニスワフ・ウラムらによって考案された，「統計的アプローチ」としての「モンテカルロ法」を概観しながら，CR系を扱うレッド－ムーン法との類似点に触れる．次に，CR系を確率過程としてモデル化するために「分子状態（molecular state）」を定義し，そのマルコフ過程の数学的取り扱いとして基本マスター方程式の基礎を概説して，その解が「定常状態」に漸近することを確認する．さらに，RMシミュレーション（レッドムーン法による分子シミュレーション）の概要を記した後，その本質的特徴が，媒質ダイナミクスを含めた内在MDシミュレーションの実行にあり，その結果，反応過程の非格子サンプリングが実現できることを述べる．また，RMシミュレーションで得られる定常状態が「準安定状態」として現れる実例として，リチウムイオン電池の負極被膜形成を挙げ，電解質分子や添加剤濃度によって決まる微細構造とその力学的安定性の解析が電池特性向上に重要であることを指摘する．最後に，レッドムーン法という「統計的アプローチ」は，ウラムらが始めた，元祖「統計的アプローチ」の進化形（拡張版）として，メソ・マクロスケールの物質状態を探求する有効な計算分子技術の一つになりうると結論づける．

1　 はじめに

2013年，著者らは，材料科学［1–4］や工業的精密合成［3］の深化に貢献しようと夢見て，反応のような希少事象を有効的に駆動する仕組みを備えた混合MC/MD（Monte Carlo/Molecular Dynamics）反応法であるレッドムーン（Red Moon, RM）法を発表しました［5–8］．この反応法は，従来の反応法［9–14］に加えて，より汎用的で，現実的な問題解決のための分子シミュレーションの新しい選択肢を提供しました．本稿では，レッドムーン法を紹介すると共に，その統計的アプローチの理論的取り扱いと，その特徴に着目して，「モンテカルロ法の新展開」について考えてみたいと思います．

1.1　 元祖モンテカルロ法とは

モンテカルロ（Monte Carlo method, MC）法とは，「シミュレーションや数値計算を，乱数を用いて行う手法の総称」で，媒質中での核粒子の運動を研究するために，1946年に，スタニスワフ・ウラムが，その基本的アイディアを考案し¹［15–17］，ジョン・フォン・ノイマン（あるいはニコラス・メトロポリス）により，その名が命名されました²（図1）［16,18］．

図1  (a) スタニスワフ・マルチン・ウラム（Stanisław Marcin Ulam, 1909–1984）(b) ジョン・フォン・ノイマン（John von Neumann，1903–1957）(c) ニコラス・コンスタンティン・メトロポリス（Nicholas Constantine Metropolis, 1915–1999）(https://ahf.nuclearmuseum.org/voices/oral-histories/nicholas-metropolis-interview/)

まさに“The Monte Carlo Method”と題された論文［15］で，ウラムらは，MC法の開発に繋がった問題群（中程度の数の粒子からなる対象系）を考察して，ボルツマン方程式やフォッカー・プランク方程式を例に挙げた後，より一般的に，点 (x, y, z) における粒子密度u(x, y, z; t) に関する，非常に簡単化された放物型偏微分方程式（反応拡散方程式）を示しています［15］．ただ，現実の物理的状況では，独立変数の多さから，“これらの方程式を扱う古典的な方法は，非常に手間がかかり「閉じた形式」での解が得られないという意味で不完全である”と述べています．そのために考案した「統計的アプローチ」がMC法という訳です．ちょうど初の電子式コンピュータENIACが開発され，弾道研究のための射撃表計算への応用が検討されていた頃のことです．

具体的には，“媒質中の核粒子が示す連鎖反応（複雑な過程（complicated process））を「マルコフ連鎖」の例と捉えれば，有限回の「（計算）実験」から，可能性のある系統のクラスまたはサンプルを得ようと試みることができるでしょう”と述べています（“one might try to perform … and obtain …”ととても控えめな表現！）［15］．実際，“各過程の確率（遷移確率）が与えられていると仮定すれば，想定された確率分布に対応する確率で，多数の偶然のゲームを行うことができるので，こうした複雑な過程の遷移確率として解釈される正の成分をもつ冪乗行列の漸近的（asymptotic）性質を経験的に研究することができる”と述べています［15］．

フォン・ノイマンは，ウラムから新しい「統計的アプローチ」のアイディアを聞き，直ちにその重要性を感得し，1947年3月，その詳細な説明を，暫定的な「計算シート」（“a tentative ‘computing sheet’”）を含めた手紙にして，当時，ロスアラモス国立研究所の理論部門のリーダーであったロバート・リヒトマイヤーに送り，この「統計的アプローチ」のデジタル処理との適合性と核分裂装置における中性子拡散と増殖の問題解決への寄与を説きました．驚いたことに，彼は“この「計算シート」に記載された制御命令は，ENIACの論理的能力を超えるものではない可能性が高いと思われる”などとENIACのスペックを想定した計算量の見積もりにも言及しています．これが電子計算機のためのMC法的手続きの最初のアウトラインでした［17］．

複雑な過程をマルコフ過程と捉えると，その確率分布の発展は，マスター方程式で記述することができます．しかし，彼らは，問題を“決定論的過程と確率論的過程の混合からなる流れの研究として数学的に捉え”ます［15］．そして（偏微分方程式やマスター方程式の）“確率行列の多重積分や乗算を避け，代わりに単一の事象の連鎖をサンプリングする”という“本質的な特徴”をもつ「統計的アプローチ」（すなわちMC法）を提唱したのです［15］．そのサンプリングの手続きは，“実際の過程に必要な回数だけ繰り返されるか，または定常分布が存在すると考えられる問題では，「（計算）実験」分布が大幅な変化を示さなくなるまで繰り返される”としました［15］．

1.2　 レッドムーン法とその似たところ

“媒質中の核粒子が示す連鎖反応”の問題では，媒質は，中性子減速材や反射材などであり，その性質は，連鎖反応を構成する吸収，散乱，分裂といった希少過程の反応断面積の値や中性子の飛行距離の指数分布に，間接的に反映しますが，基本的に「真空」です．分子科学では，媒質が分子という「粒子」集合体（分子凝集系）となり，そこには，個々の分子の振動運動（または秤動運動）が生み出す，分子（群）の動的過程が溢れています．（気体状態では，並進，回転，振動運動が区別されますが，凝集状態では，並進も回転も“振動（ときに秤動）”です）．

このとき，化学反応や拡散は，頻繁には起こりませんが，ときどき起こる希少事象です．ところが，分子の立体化学や会合構造，均一あるいは不均一な準安定構造などは，（分子振動の時間スケールで見たときの）長い時間スケールで生まれ，度重なる複雑な分子過程の結果として形成されるため，このような現象を従来の分子シミュレーション手法を直接適用して，計算で実現することは困難です．

レッドムーン法で扱う希少過程は，媒質中の核粒子が示す吸収，散乱，分裂に対応して，溶質・溶媒分（群）が生み出す，非常に多くの分子振動群（孤立・集団モード）の中に埋もれながら，希れに起こる化学反応や拡散です．分子凝集系で起こる反応では，「媒質」の捉え方と粒子の構造が，さらに複雑になっていますが，ウラムらの対象系と似ています．

素反応（elementary reaction, ER）系においてさえ，化学反応や拡散（以下，単に反応）自体は，多くの場合，とても希れにしか起こらない出来事です（文献6,19などを参照してください）．しかし，実際には，複合反応（complex reaction, CR）系のように，複数の種類の反応が，同時に競合して，生成物の組成やその空間分布に反映したり，どの反応が速く（あるいは早く）起こるかによって，立体化学に関連する幾何学的問題として現れたりします．

化学平衡にある，相空間のもっともらしい状態を実現する計算効率のよい方法として考えられるのは，遷移状態領域によって分割されて存在する安定相状態の，比較的安定な領域間で，系全体を試行的に遷移させ，確率論的アプローチを適用して，希れな遷移の発生後の状態が実現されるべきか否かを判断するという方法です．20世紀最後の10年間に，レッドムーン法が想定するような似た考え方が，気相化学平衡を扱うために導入され［11］，その後より複雑な化学系に応用されました［12–14］．しかし，溶液中の化学反応系では，さらに，このような希少過程に関連した現象を，反応分子の「原子論的」な複雑さと柔軟性をそのまま残して計算するという問題に取り組む必要がありました．レッドムーン法の開発に関連した，2000年以前，またそれ以降の研究動向については，文献6, 19を参考にしてください．

2　 複合反応系の「統計的アプローチ」の基礎

2.1　 確率過程としての複合反応系

従来，計算科学で，分子の反応を考えるとき，単独のERや多段階の連続したER群に注目して，孤立系あるいは，無限希釈系としてモデル化することが主流でした．しかし，CR系を扱う際には，それをたくさんの反応分子群が含まれた広がりをもった対象系と捉え，そのあちこちで，複数の希少過程を別々に示す，多数の反応分子群からなる分子集合体としてモデル化することが便利です．実際には，反応分子群そのものに加えて，触媒反応では，さらに触媒分子や金属原子が，凝集系では溶媒分子群がCR系には含まれることになります［6,19］．

 分子状態

分子状態（molecular state）という用語は，例えば，一つのER

aA + bB → cC + dD　（始状態i → 終状態f）　(R)

のみを考える場合には，反応「系」（すなわちaA + bB）あるいは，生成「系」（cC + dD）に対して使用される用語の「系」と，また，始「状態」iまたは終「状態」fに対して使用される用語の「状態」［e.g.,20–22］［23］と同義です［6,19］．ただ，分子状態は，一般に，多種類の成分分子が複数の反応を起こすCR系の場合には，全系のあちこちに存在する，或る反応が可能な或る同種の成分分子（群）を，それぞれ全数含んだ，独立した「系」または「状態」を区別して指すために用います．

分子状態とは，或る反応(R)に関与する分子種AとB，分子種CとDを，単に濃度[A]などとスカラー的に捉えずに，a個の分子Aとb個の分子Bの組が，全系に空間的に広がって分散して存在している様子を思い浮かべて捉えようとした，或る反応に着目したときの，その反応に関与する同種分子（群）達に対する「状態」の名称です．また反応(R)の（aA + bB）におけるAという同種分子群とBという同種分子群のように，同種分子群が複数の場合や，或る分子状態に含まれる成分分子（群）が同時に別の分子状態に含まれる場合もあります．分子状態では，従来の「系」や「状態」に，空間的に点在する個々の分子群の様態を埋め込んだように思い描いています．

 風景類似性の仮定（assumption of landscape similarity）

「分子凝集状態」における反応拡散現象は，微視的な分子運動の時空間スケールのもとで，反応や拡散などの各希少過程そのものが，互いに無相関と見做すことができれば，マルコフ過程と見做せます．つまり，或る希少事象が発生しても，その周辺では激しい分子運動によって相空間情報の混合（mixing）が急速に起こり，それまでの履歴はすぐに失われてしまい，平衡状態が生まれます［24］．これは反応系の特徴であり，遷移状態（transition state (TS)）理論［25,26］が非常にうまく成り立つ事実の起源です．従って，同じ種類の希少事象，例えば同じ反応が起こったときの，反応分子群を取り巻く「風景（landscape）」の類似性は，TS理論と同様にレッドムーン法でも仮定されていると言えます．

 分子状態のマルコフ過程

こうして，CR系に存在するM個の分子状態を区別して，各分子状態をm (= 1, 2, ..., M)で表します．それぞれの成分分子（群）の数x_mをベクトル的に表現してCR系全体の全分子状態をベクトル変数x = (x_m) = (x₁, x₂, ..., x_M)で記述します．より具体的な例は，文献6,19あるいは，4.1小節を参照してください．

xは，反応などが起こる度に，時間的に変化しますが，レッドムーン法では，そのダイナミクスを確率的に起こる変化と見做し，その，すべての分子状態の確率ダイナミクスを，マルコフ過程x(t)（ベクトル）として捉えます［27,28］．化学反応などにより分子状態mがnへ確率的に遷移するとして，ひとまず，遷移は離散的な時刻$t=s\cdot \mathrm{\Delta }t$（ここで，整数s = 0, 1, …, ∞）で順番に起こると仮定します．ここで，$\mathrm{\Delta }t$は連続するステップの間の基本的な時間間隔とします．こうして，離散的な時刻sで順番に起こる遷移x_sをマルコフ過程

  

$$\begin{array}{c}\hfill {\mathbf{x}}_{s=0}\to {\mathbf{x}}_1\to {\mathbf{x}}_2\to {\mathbf{x}}_3\to \cdots \to {\mathbf{x}}_s\to \cdots \end{array}$$(1)

と見做します．

2.2　 分子状態の確率分布の基本マスター方程式

マルコフ過程の標準的取り扱いに従って，チャップマン－コルモゴロフ方程式を基礎に，CR系の基本マスター方程式を求めることができます［29,30］．

  

$$\begin{array}{c}\hfill \frac{\partial P(n,t)}{\partial t}=\sum _{m=1}^M[w_{nm}\cdot P(m,t)-w_{mn}\cdot P(n,t)]\end{array}$$(2a)

  

$$\begin{array}{c}\hfill \iff \frac{\partial \mathbf{P}(t)}{\partial t}=\mathbf{W}\cdot \mathbf{P}(t).\end{array}$$(2b)

この方程式は，確率分布（ベクトル）$\mathbf{P}$，すなわち時刻tにおけるx = (x_n)の実現確率の集まり

  

$$\mathbf{P}={}^t(P(n))={}^t(P(1)P(2)\cdots P(M))$$(3a)

  

$$={}^t\left(P_n\right)={}^t\left(P_1{P}_2\cdots P_M\right),$$(3b)

の時間発展を支配します．ここで$P_n$は分子状態nの部分確率分布を表します．また，式(2b)において，$\mathbf{W}$は遷移行列（時間発展演算子），

  

$$\mathbf{W}=\left(W_{nm}\right)\equiv (w_{nm}-{\delta }_{nm}\sum _l^{}{w}_{lm}),$$(4)

を表し，対象CR系の反応スキーム（その例については4.1小節を参照）に応じて決まる遷移率行列$\mathbf{w}$を用いて定義されます．ここで，$w_{nm}$は分子状態mからnへの確率的遷移，すなわち反応m→nの単位時間あたりの遷移率（transition rate）を表し，$\mathbf{W}$は，次の共通の性質をもちます．

  

$$W_{nm}\ge 0(n\ne m),$$(5a)

  

$$\sum _n^{}{W}_{nm}=0(\text{各}m\mathrm{に対して})．$$(5b)

ここでは，或る整数Nに対して${\mathbf{W}}^N$のすべての要素が正になると考えておきます．このとき$\mathbf{W}$は正則（またはエルゴード的）であると言います（これは任意の分子状態間で遷移可能という意味です）［29,30］．

今後，反応m→nを，単に反応nmと表記します［6,19］．nmにおけるnとmの順序は，化学反応式の通常の順序とは異なります．この「反応」を，振動などのありふれた分子運動よりも遅い時間スケールで起こる，一般的な分子過程と見做しても，ここでの理論的な一般性は失われません．また，特に断らずに，頻繁に起こらない，時間のかかる過程をすべて「反応」という用語を使うことがあります．

2.3　 基本マスター方程式の解とその特徴

基本マスター方程式(2)については，物理学や化学におけるマルコフ過程の数学的取り扱いが，よく知られています［29,30］．

 Wの固有値と固有ベクトル

CR系のモデル化において，遷移行列$\mathbf{W}$は一般に非対称実行列となります³［29,30］．このとき，一般に，$\mathbf{W}$の左固有ベクトル｛${\mathit{\chi }}_j$｝と右固有ベクトル｛${\mathit{\phi }}_j$｝とは異なります．今，右固有値方程式

  

$$\mathbf{W}{\mathit{\phi }}_j={\lambda }_j{\mathit{\phi }}_j(j=1,2,\dots ,M)$$(6)

を考えて，初期確率分布$\mathbf{P}(0)$を，

  

$$\mathbf{P}(0)=\sum _{j=1}^M{c}_j{\mathit{\phi }}_j$$(7)

とスペクトル分解して表すと，基本マスター方程式の解，すなわち，時刻tにおける確率分布$\mathbf{P}(t)$を，

  

$$\mathbf{P}(t)=\sum _{j=1}^M{c}_j{e}^{{\lambda }_{j}t}{\mathit{\phi }}_j$$(8)

と表すことができます．

他方，性質(5b)は，左固有値方程式

  

$${\mathit{\chi }}_j\mathbf{W}={\mathit{\chi }}_j{\lambda }_j(j=1,2,\dots ,M)$$(9)

が，ゼロ固有値（${\lambda }_1$とします）とその左固有ベクトル${\mathit{\chi }}_1$ = (1, 1, 1, …, 1)をもつことを意味しています．従って，同じゼロ固有値${\lambda }_1$をもつ右固有ベクトル${\mathit{\phi }}_1$も存在します［29,30］．

 確率分布の定常分布への緩和

一般に，性質(5a)と(5b)をもつ$\mathbf{W}$の固有値${\lambda }_j$は，非正値条件

  

$${\lambda }_j\le 0(j=1,2,\dots ,M)$$(10)

を満たすことが示せます［29,30］．従って，式(8)において，ゼロ固有値以外の右固有ベクトル成分は，十分長い時間が経つと減衰するので，確率保存条件を考慮すると，右固有ベクトル${\mathit{\phi }}_1$に対して，

  

$$c_1=1$$(11)

が成り立ち，

  

$$\mathbf{P}(t)\to {\mathit{\phi }}_1\text{（}t\text{→長い時間）}$$(12)

となることが分かります［29,30］．

こうして，“「（計算）実験」分布が大幅な変化を示さなくなるまで繰り返せ”ば［15］，CR系で扱うような，ハミルトン運動方程式と確率的遷移で記述される「孤立」系の問題［15］では，定常分布（もちろん平衡分布を含みます）が存在するであろうことが分かります［29］．

ウラムらは，彼らの「統計的アプローチ」を用いると，エルゴード性と定常分布を純粋に適用して，“遷移確率として解釈される正の係数（coefficients）をもつ冪乗行列の漸近的性質を経験的に研究することができる”と述べています［15］．

3　 レッドムーン法：決定論的過程と確率論的過程が混合した一つの「統計的アプローチ」

RMシミュレーション（レッドムーン法による分子シミュレーション）では，MD法とMC法の一連のステップを組み合わせたサイクル（RMサイクル）を繰り返して，CR系の十分長い時間に亘る分子シミュレーション（「反応」という希少過程に伴う分子凝集系の微視的構造の効率的サンプリング）を実現して，その漸近的性質を調べます［6,19］．

3.1　 内在MDシミュレーションと非格子サンプリング

CR系の確率ダイナミクスを支配する基本マスター方程式(2)を，実際に多自由度系について解くことは，ウラムらが述べたように“不完全”で，現実的ではありません．そのために分子シミュレーションを使って，有効的にそれを解くことを考えます．このような場合，マスター方程式に基礎をおいた運動学的シミュレーション（KMC法［31,32］と呼ばれています）が，古くから用いられています．

 内在MDシミュレーション

レッドムーン法でも，式(2)を解いて$\mathbf{P}(t)$を求めるのではなく，背景となる（それゆえ基礎となる）内在（intrinsic）分子動力学（iMD）シミュレーションを行いながら，全体として$\mathbf{P}(t)$を構成する確率的軌道のサンプルを求めます［6,19］．iMDシミュレーションは，分子群の振動運動を追跡し続け，一旦起これば瞬時に起こるものの，頻繁には起こらない希少過程（「反応」など）に対して，適時適切にMC処理できるように備えます．

レッドムーン法の本質的な特徴は，RMシミュレーション全体のバックグラウンドで，自由度の削減や座標サンプリング空間の簡略化，そしてTSの探索などを介さずに，媒質ダイナミクスも含んだiMDシミュレーションを実行することだと言えます．RMシミュレーションでは，iMDシミュレーションが，CR系の全原子に対する決定論的過程を受け持ちます．これは，ウラムらの“媒質中の核粒子が示す連鎖反応”の問題では想定していない（というか，その必要がなかった）取り扱いでした．

 反応候補の非格子サンプリング

こうして，iMDシミュレーションを通して，複数の反応群｛nm｝の中から選ばれた或る反応nmを起こしうる（つまり，構造的に反応を引き起こすことが可能な）分子状態mの反応候補種（つまり候補分子（群））を追跡し，その中から，適切な候補種を選択して，反応nmに伴う3次元的変化を許容できる分子状態nの分子（群）へと構造変化を試行します（確率論的過程）．

その際，反応確率は，当然，反応候補種数に依存します．iMDシミュレーションに沿って，特定の型の反応nmの種類とその候補種の数｛$N_{nm}^{\mathrm{can}}$｝が，随時更新され，候補リストに保存されます．$N_{nm}^{\mathrm{can}}$は，iMDシミュレーション中に，原理的に，CR系全体で反応nmを示すことが可能な注目種の，その時刻の濃度に応じて合理的に得ることができます［6,19］．

このiMDシミュレーションに基づいた，完全な実空間非格子処理を通してなされる反応の試行遷移は，メトロポリス因子［32,33］を含んだ遷移率

  

$$w_{nm}\propto {\Omega }_{nm}\cdot \min \{1,\exp \left[-\beta \mathrm{\Delta }{U}_{nm}\right]\}(>0)$$(13)

  

$$=\{\begin{array}{cc}\hfill {\Omega }_{nm}\hfill & \hfill \mathrm{\text{(}}\mathrm{\Delta }{U}_{nm}\le 0),\hfill \\ \hfill {\Omega }_{nm}\cdot \exp \left[-\beta \Delta U_{nm}\right]\hfill & \hfill \mathrm{\text{(}}\mathrm{\Delta }{U}_{nm}>0),\hfill \end{array}$$(14)

を使って，RMシミュレーションの遷移過程（反応過程）として採択するか棄却するかを判断します［33］．ここで，$\mathrm{\Delta }{U}_{nm}$は，分子状態mにある，或る分子（群）を始状態として，分子状態nにある，試行的に構造変化させた，対応する分子（群）（試行的な終状態）とのポテンシャルエネルギー差

  

$$\mathrm{\Delta }{U}_{nm}=\mathrm{\Delta }{U}_{nm}^{\mathrm{\text{MM}}}+\mathrm{\Delta }{U}_{nm}^{\mathrm{\text{reac}}}$$(15)

を表し，$\mathrm{\Delta }{U}_{nm}^{\mathrm{\text{MM}}}$は，標準的MM力場の枠組みで得られた分子状態mの構造配置とnのそれにおける局所ポテンシャルエネルギーの素朴な差です．他方，$\mathrm{\Delta }{U}_{nm}^{\mathrm{\text{reac}}}$はMM力場の枠組みでは通常考慮されない，上の始状態と終状態のそれぞれの配置に近い，分子内力場で表された2つの平衡構造配置（エネルギー「ゼロ」点）の間の断熱ポテンシャル補正で，近似的な反応ポテンシャルエネルギーです［6,19］．また，式(13)と(14)で，反応nmの相対的選択重み${\Omega }_{nm}$は，

  

$${\Omega }_{nm}\sim N_{nm}^{\mathrm{can}}•\exp \left(-{\beta •\mathrm{\Delta }E}_{nm}^{\mathrm{a}}\right),$$(16)

のように，反応候補種数$N_{nm}^{\mathrm{can}}$とアレニウス型の指数関数因子（$\mathrm{\Delta }{E}_{nm}^{\mathrm{a}}$は，反応nmの活性化自由エネルギー）との積によって見積もられる反応確率です．

このように，レッドムーンアプローチは，「決定論的過程」と「確率論的過程」の混合からなる，一つの「統計的アプローチ」で，現実の“分子”集合体に対して，成分分子の数x_mのサンプリングのみならず，実空間で微視的に蠢く3次元構造をもサンプリングできる，実空間非格子KMCアプローチと言えます．

3.2　 媒質ダイナミクスの遷移行列への影響

レッドムーン法では，各反応過程の個別性を保ちながら，iMD軌跡上に現れる個々の反応過程について，反応種類ごとの事前データ（$\mathrm{\Delta }{E}_{nm}^{\mathrm{a}}$と$\mathrm{\Delta }{U}_{nm}^{\mathrm{\text{reac}}}$）を用いて，確率的に採択棄却を判定します．この近似の，計算コストに対する「経済」効果により，不均質な化学状態の不均一な微細構造が，分子シミュレーションで調製可能になるのです．その際，CR系全体（媒質分子（群）も含む）のiMDシミュレーションが反応候補種を探索するため，分子（群）の対称性や反応条件，各成分種のどの分子（群）が，いつどこにいるのかに応じて，$\mathrm{\Delta }{U}_{nm}$の値に影響し，$W_{nm}$の値は試行毎に変わります．従って，$\mathbf{W}$を暗に定数行列とした2.2，2.3小節の議論は，媒質ダイナミクスの影響も含めたランダム行列などを用いて拡張する必要があります．

また，iMDシミュレーションの実空間の「軌跡」は，従来のMDシミュレーションの「軌跡」のように，本来の意味で，時間的にも空間的にも滑らかに繋がっている訳ではありません．しかし，「RM軌跡の集合」，すなわち多数のRMシミュレーションによって得られた軌跡全体は，全部集まって初めてマスター方程式の解である確率分布$\mathbf{P}$と等価になると言えます．この関係は，ファインマン経路積分［34］の各経路の集合と，時間依存シュレーディンガー方程式［34–36］の解である波動関数で表される時間相関関数との関係に似ています．

4　 応用：リチウムイオン電池の負極被膜形成

リチウムイオン電池において，非プロトン性有機溶媒，例えば炭酸エチレン（EC）や炭酸プロピレン（PC）からなる液体電解質では，充放電中における電解質分解による劣化を防ぐための被膜層を形成することが非常に重要です［37–41］．これは初期充電中に負極電解質界面で，複雑な化学反応によって形成され固体電解液相間層（solid-electrolyte interphase (SEI) layer）（あるいはSEI膜）と呼ばれています．SEI層は，Liイオンなどのキャリアイオン伝導性をもつ一方で，電子伝導性をもたないことが理想的で，二次電池の寿命や安全性に深く関わっているため，その形成機構や微細構造の理解は，高性能化を実現するために不可欠です．そこで，実験研究によるSEI層形成過程の直接観察を補完するために，電解液の溶媒選択性とSEI層形成に対する添加剤効果についてレッドムーン法を用いて調べました［37–41］．

4.1　 リチウムイオン電池のモデル化

この研究[37,38]では，グラファイト様炭素系負極と，ヘキサフルオロリン酸リチウムLiPF₆を溶解させたEC系およびPC系電解液（塩濃度～1.1 mol/L）を想定しました．有機溶媒の中でもPCは，可燃性をもつものの，広い液温範囲（218～513 K）をもつ魅力的な非水電解質溶媒です．これらのモデル系（図2）として，負極は，グラフェンシートの積層体，つまり10層の炭素原子層でモデル化し，層間距離は，LiG（lithium-intercalated graphite）の層間距離に対応させて3.7 Åと取りました．

図2  （上）リチウムイオン電池の初期状態における負極界面付近のモデル（炭素系電極Cと電解液（炭酸エチレンC₃H₄O₃（EC）（あるいは，炭酸プロピレンC₄H₆O₃（PC）），とリチウム塩LiPF₆）の混合物）．（下）初期充電過程における反応スキーム（開環：21，分解：32と65，重合：42）（数字は各分子状態（1～6）を表す）．

また，SEI膜形成のRMシミュレーションに必要十分な反応スキーム（開環，分解，重合からなる主要反応リスト）（図2）を考慮しました［6,19］．実際，反応スキームは，6個の分子状態m = 1～6に対して，EC（/PC）分子の還元を含む反応群｛nm｝=｛21, 32, 42｝と塩起源のPF₆⁻イオンの分解反応｛65｝の計4つから成り立ち，それらを表す遷移率行列$\mathbf{w}$と遷移行列$\mathbf{W}$は図3のようになります．

図3  (a) 遷移率行列$\mathbf{w}$と (b) 遷移行列$\mathbf{W}$（式(4)）．各遷移率w_nmの数値は，文献37のTable 1を参照せよ．

4.2　 SEI膜という準安定状態の形成

SEI膜形成の二次電池の性能上の議論については，文献6,19で引用した，著者らの論文や実験化学者との共著論文に記されています．ここでは，分子シミュレーション的視点に立って，基本マスター方程式で裏付けられた「統計的アプローチ」RMシミュレーションが生み出す定常状態としてのSEI膜について述べます．図4は，温度298 KのECあるいはPCを主成分とした2種類の電解質溶液に対するRMシミュレーションの結果で，5つのRMサイクルにおける反応生成物の状態変化の典型的なスナップショット（側面図）です．RMサイクルが進むに連れて電解質の還元生成物がSEI層として成長する様子が分かります（文献37の補足資料動画参照）．さらに，図5では，表面分子数密度のRMサイクル変化から，2000 RMサイクルには，各成分量が一定値に達し，定常状態に至ることが分かります．EC系電解質では，二量体Li₂BDC（二炭酸ブチレン二リチウム）の数が大幅に増加することが判明しました．その際に得たSEI膜成分の負極表面鉛直方向の質量密度分布（図6）からは，無機塩は負極近くに分布し，有機塩は主に外層に分布する傾向が分かり，実験観察を支持する準安定状態（あるいは非平衡定常状態）が形成されていることが示せました．

図4  (a) EC系，および (b) PC系電解質におけるSEI膜形成過程．RMシミュレーション中の反応生成物の凝集状態の変化を示す典型的な5つのスナップショット（バルク電解質は非表示）［37］．Copyright 2014, American Chemical Society.

図5  (a) EC系，および (b) PC系電解質における，SEI膜形成過程における溶媒分子と反応生成物の表面分子数密度のRMサイクルに対する変化［37］．Copyright 2014, American Chemical Society.

図6  (a) EC系，および (b) PC系電解質における，RMシミュレーションの定常状態（2000 RMサイクル）における反応生成物の質量密度分布［37］．Copyright 2014, American Chemical Society.

4.3　 僅かな構造差が生む構造安定性の大きな差

図7は，負極表面近傍のSEI層の界面構造を正面から見たスナップショットです．EC系とPC系電解質のそれぞれから形成された2つのSEI層を比較すると，PC系SEI層の界面構造に点在する孔径は，EC系SEI層の孔径よりも明らかに大きいことが分かります．これは，EC分子に対してPC分子には余分なメチル基が存在し，この僅かな分子構造の差が，PCの分解生成物間の凝集を阻害し，SEI層の構造安定性に負の影響を与えたと考えられます．

図7  (a) EC系，および (b) PC系のSEI膜界面構造の正面スナップショット（深さ0≦z≦20Åをバルク電解質側から見て可視化）．これらの10個の構造は，異なる10個の初期配置（初期分布）から得られた［37］．Copyright 2014, American Chemical Society.

また，SEI層の「安定性」は，そのトポロジカルな特性（空隙率や膜厚など）に大きく関係し，SEI層の機械的安定性（強度と言っても良いかもしれません）が，充電容量（可逆容量）とサイクル寿命の向上に繋がることが示唆されました［6,19,41］．実際，SEI層の「安定性」には，必ずしも，エネルギー的安定性は直接関係しません．

例えば，添加剤炭酸フルオロエチレン（fluoroethylene carbonate, FEC）を添加したNaPF₆/PC系電解液においては，実験的には，1.0 vol%程度の低FEC濃度で，サイクル特性が最適化されることが知られています．一方，レッドムーン法で得たSEI層は，1.0 vol%で，エネルギー的には最も不安定であることが分かりました．SEI層は，FEC濃度に依存して，サンプル毎に独立に決まる化合物成分比をもつ分子集合体ですが，この場合，その構造は，エネルギー的な安定性と引き換えに，むしろ機械的な安定性を優先して実現した構造であると言えます［41］．従って，この場合，エネルギー的安定性のFEC濃度依存性の比較はあまり有意義ではありません．これは，シリカガラスが結晶シリカよりもエネルギー的に不安定であるにも関わらず，十分に「安定」に存在している状況に似ていると言えるかもしれません．

こうして，レッドムーン法で調製したSEI層で観察されたエネルギー的安定性の一見「逆」と感じるFEC濃度依存性は，SEI層が準安定状態として形成されるためであることから理解できます．結果として，二次電池の被膜「安定性」には，機械的安定性が不可欠であると結論づけることができます［6,19］．

5　 考察

5.1　 「モンテカルロ法」とは…

平衡状態（あるいは定常状態）の性質を求めたいときに，よく用いられる“MC法”ですが，初期濃度で混ぜた始状態の化学物質が，反応を繰り返して終状態（定常状態）に到達する過程は，一種の緩和過程であり非定常状態にあります．

MC法を，広く「シミュレーションや数値計算を，乱数を用いて行う手法の総称」と解すれば，レッドムーン法は，MC法の一種である，あるいはMC法を使っている，と言えます．しかし，MC法を,もう少し絞って，マルコフ連鎖モンテカルロ（Markov-Chain Monte Carlo, MCMC）法などと解すると，途端に「えっ？モンテカルロ法とは言えないかも…」と分からなくなってしまいます．

レッドムーン法を始めた初期の頃，その研究成果を学会等で発表すると，「詳細釣り合い条件（detailed balance condition, DBC）［29,30］

  

$$w_{nm}\cdot P(m,t)=w_{mn}\cdot P(n,t)$$(17)

は成り立っているのですか？」という質問を頻繁に受けました．ちょうどDBCを満たさない新しいMCMC法（諏訪藤堂法［42］）が発表された直後で，MC法に大きな進展があったことも影響したのかもしれません．私自身が不勉強であったことも手伝って「（この人は一体何を尋ねているのだろう？）」と内心とても困惑したものでした．その理由は，緩和過程のような非定常状態においては，DBCはもちろん，釣り合い条件（balance condition, BC）［29,30］

  

$$\sum _{m=1}^M{w}_{nm}\cdot P(m,t)=\sum _{m=1}^M{w}_{mn}\cdot P(n,t)$$(18)

も満たされる訳はないからです．

また，RMシミュレーションでは，QuenchMDやRelaxMDのステップを入れることがあるのですが，この点についても指摘されることがあります．もし「RMシミュレーション中に，DBCが成り立つような状態に一旦遷移させる」のだとしたら（実際は違います！），DBCが満たされなくなってしまう“根拠”として，これらのステップを入れる点を指摘されるのです．もし，もともとの確率的遷移過程がDBCを満たすようになされたのなら，その状態から，Quenchしたり，Relaxさせたりしたら，元も子もないという指摘なのかもしれません⁴．

文献等で「メトロポリス法［33］などの多くのMC法では…」と言うとき，このMC法は，DBCを仮定しています．また，「シミュレーテッド・アニーリング法［43］では，“温度”パラメータT依存性があるため，Tが変化する際に詳細釣り合い条件が破れるので，MC法とは異なり…」などと書かれていたりします．ここで言っているMC法は，絞った意味での統計力学的MC法です⁵［44］．

現在，物質科学の多くの文献で，一般に「モンテカルロ法」というと，ほぼMCMC法を指しているように思います．つまり，「カノニカル分布を定常状態にもつマルコフ過程を用いて平衡状態の統計的性質を調べる統計力学的MC法」［44］です．その際，対象系の動的な性質は遷移確率の具体的な取り方に依存し，一般性を失ってしまいます．しかし，そこで用いた遷移確率が何らかの物理的な過程に対応していると考えられるとその緩和過程も意味をもってくるのです［44］．

5.2　 遷移率が変化したり分子状態が結合したり

図3(b)の遷移行列$\mathbf{W}$を思い出してください．$\mathbf{W}$に限らず，一般に行列が適当な相似置換によって，ブロック対角行列

  

$$\mathbf{W}=\left(\begin{array}{cc}\hfill \mathbf{A}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{B}\hfill \end{array}\right),$$(19)

に変換できる場合，完全可約あるいは分解可能と呼ばれます［29,30］．ここで，$\mathbf{A}$と$\mathbf{B}$は，$\mathbf{W}$より低次の2つの正方行列です．この場合，$\mathbf{W}$の状態は，その間で遷移が不可能な2つの副状態に分解されます．ここでは$\mathbf{A}$も$\mathbf{B}$も，それぞれ，性質(5a)と(5b)を満たすとして，それぞれゼロ固有値の右固有ベクトル${\phi }_1^{\mathrm{A}}$と${\phi }_1^{\mathrm{B}}$が存在します．従って，$\mathbf{W}$は，ゼロ固有値の線形独立な固有ベクトルを少なくとも2つもちます［29,30］．

  

$$\left(\begin{array}{cc}\hfill \mathbf{A}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{B}\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\mathit{\phi }}_1^{\mathrm{A}}\hfill \\ \hfill \mathbf{0}\hfill \end{array}\right)=\mathbf{0},$$(20a)

  

$$\left(\begin{array}{cc}\hfill \mathbf{A}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill \mathbf{B}\hfill \end{array}\right)\left(\begin{array}{c}\hfill \mathbf{0}\hfill \\ \hfill {\mathit{\phi }}_1^{\mathrm{B}}\hfill \end{array}\right)=\mathrm{𝟎}.$$(20b)

この場合，$\mathbf{W}$の状態は，$\mathbf{A}$と$\mathbf{B}$を遷移行列にもつ2つの基本マスター方程式に支配され，全系は，互いに遷移し合わない2つの非相互作用系（部分系）から成り立っていることを意味します［29,30］．この状況は，全CR系の温度が十分低かったり，各反応に関与する分子濃度が十分低く，互いに全く影響を及ぼし合わなかったりする場合などに対応するでしょうか．

ところが，一般に，RMシミュレーションでは，iMDシミュレーションをしながら反応候補種を探索するため，3.2小節でも触れたように，媒質ダイナミクスが，$W_{nm}$の値に陰に影響を与えたり，形の上では完全可約(19)であるにも関わらず，$\mathbf{W}$の各副状態の成分種同士，すなわち，非相互作用系であった分子状態同士が，実空間配置を介して相互に結合したりします（隠れた結合）．こうした事実が，分子科学的CR系における独特な準安定“定常状態”が実現することに繋がります．その結果として，4節で紹介したような特徴的な多孔性のSEI膜が出現するのです．

6　 おわりに　さらなる温故知新の展開へ

レッドムーン法による研究は，MC法の新展開というか，「ルネッサンス」のように見えなくもないなと思います．というのも，レッドムーン法では，ウラムらの“媒質中の核粒子が示す連鎖反応”の問題では想定する必要がなかった媒質ダイナミクスが，分子科学的CR系を念頭にした場合，複雑多様かつ不可避必須になり，媒質分子（群）と反応分子（群）との微視的相互作用を介して，生成物の組成や時空間分布に大きな影響をもたらすからです．その際，媒質分子（群）そのものの微視的特徴自体も，CR系の準安定な定常状態の構造に反映します［6,19］．

レッドムーン法は，ウラムらが始めた，元祖「統計的アプローチ」の進化形と言えるかもしれません．なぜなら，それは，この80年の間に飛躍的に進展した情報技術と，その結果実現した，電子状態，振動状態，粗視状態のモデリングを時間空間的に繋ぐマルチスケールモデリング［45］とが相まって，初めて可能になった計算分子技術だからです．今後は，さらに可能な温故知新を展開して，メソ・マクロスケールの物質状態を探っていきたいと思います．

 謝辞

本稿で紹介した内容は，名古屋大学情報学部・大学院情報学研究科と未来社会創造機構において，学部生・大学院生を初めとした多くの共同研究者と一緒に進めた研究教育活動を通して考えてきた事柄と考察から生まれたものです．皆さんとの議論に心から感謝いたします．また，私達の理論的研究を常に温かく見守り，ときに厳しく，ときに寛大に指導してくださった諸先輩，多くの先生方にも，この場をお借りして厚くお礼申し上げます．最後に，本稿をまとめる機会を与えてくださった「アンサンブル」担当編集委員の森俊文先生，そして編集委員会の皆様に深く感謝申し上げます．あらためて色々と考えを巡らし，それを深める貴重な機会となりました．読者の皆様には，ご意見ご感想など，是非，お聞かせ頂けたらと思います．宜しくお願いします．

注

1  文献17，特に，ウラムが1946年の事として1983年に語った未発表の言葉：“The first thoughts and attempts I made to practice [the Monte Carlo method] were suggested by a question which occurred to me in 1946 as I was convalescing from an illness and playing solitaires. … I wondered whether a more practical method than “abstract thinking” might not be to lay it out say one hundred times and simply observe and count the number of successful plays.”などや，フォン・ノイマンが，1947年にロバート・リヒトマイヤーに宛てた手紙の中の言葉：“This is the letter I promised you in the course of our telephone conversation on Friday, March 7th. (改行) I have been thinking a good deal about the possibility of using statistical methods to solve neutron diffusion and multiplication problems, in accordance with the principle suggested by Stan Ulam. …”，（図1中の写し）．

2  文献18，論文中，次の言葉があります：“It was at that time that I suggested an obvious name for the statistical method—a suggestion not unrelated to the fact that Stan had an uncle who would borrow money from relatives because he “just had to go to Monte Carlo.” The name seems to have endured.”（ IはMetropolis ）．

3  この特徴は，量子力学では，ハミルトニアンがエルミートであり，基準振動解析では，ヘシアンが対称であることと対照的です．

4  CR系全体の開始状態（非平衡状態）から試行過程が繰り返される緩和過程の途中で，メトロポリススキームを適用しているのは，確率的軌道を効率的に定常状態に到達させようとするためです．ここでの定常状態は，局所的な平衡状態と言ってもよいかもしれません．ただし，本来，統計力学的MC法では，平衡状態における微視的遷移が満たすべき条件として，十分条件であるDBCを満たす遷移率の取り方の一つとして，メトロポリススキームが適用されます．

5  例えば，文献44．本書では，統計力学的MC法を超えて，“時間的に非定常な系の巨視的動的性質”の解明をも視野に入れて解説されていて，とても勉強になります．

参考文献

• [1]   M.  Nagaoka,  M.  Takayanagi,  N.  Takenaka,  Y.  Suzuki,  K.  Matsumoto,  N.  Koga,  S. K.  Sankaran,  P.  Uppula and  Y.  Kitamura, In Molecular Technology, Vol. 3: Materials Innovation, H. Yamamoto, T. Kato, Eds. Wiley-VCH, New York, 201 (2019).
• [2]  Research Project (JST-CREST): “Establishment of Computational Molecular Technology towards Macroscopic Chemical Phenomena -Molecular Control of Complex Chemical Reactions, Stereospecificity and Aggregate Structures-” https://pubs.aip.org/aip/acp/article-abstract/1702/1/090029/998219/Computational-molecular-technology-towards?redirectedFrom=fulltext Accessed 25 Feb. 2025.
• [3]  Research Area of JST-CREST “Establishment of Molecular Technology towards the Creation of New Functions” https://www.jst.go.jp/kisoken/crest/research_area/completed/bunyah24-4.html Accessed 25 Feb. 2025.
• [4]  Elements Strategy Initiative for Catalysts and Batteries (ESICB) Interplay Between Experimental and Theoretical Studies. https://www.mext.go.jp/elements-strategy/en/digest/p9/ Accessed 25 Feb. 2025.
• [5]   M.  Nagaoka,  Y.  Suzuki,  T.  Okamoto and  N.  Takenaka, Chem. Phys. Lett., 583, 80 (2013).
• [6]   M.  Nagaoka, WIREs Comput. Mol. Sci., 14, e1714 (2024).
• [7]   Y.  Suzuki and  M.  Nagaoka, J. Chem. Phys., 146, 204102 (2017).
• [8]   T.  Fujie,  N.  Takenaka,  Y.  Suzuki and  M.  Nagaoka, J. Chem. Phys., 149, 044113 (2018).
• [9]   P. G.  Bolhuis and  D. W. H.  Swenson, Adv. Theory Simul., 4, 2000237 (2021).
• [10]   H.  Fujisaki,  M.  Shiga and  A.  Kidera, J. Chem. Phys., 132, 134101 (2010).
• [11]   M. S.  Shaw, J. Chem. Phys., 94, 7550 (1991).
• [12]   W. R.  Smith and  B. J.  Triska, Chem. Phys., 100, 3019 (1994).
• [13]   J. K.  Johnson,  A. Z.  Panagiotopoulos and  K. E.  Gubbins, Mol. Phys., 81, 717 (1994).
• [14]   N.  Hansen,  S.  Jakobtorweihen and  F. J.  Keil, J. Chem. Phys., 122, 164705 (2005).
• [15]   N.  Metropolis and  S.  Ulam, J. Amer. Stat. Assoc., 44, 335 (1949).
• [16]   J.  Von Neumann and  S.  Ulam, Bull. Amer. Math. Soc., 53, 1120 (1947).
• [17]   R.  Eckhardt, LA-Science, 15, 131 (1987).
• [18]   N.  Metropolis, LA-Science, 15, 125 (1987).
• [19]   長岡 正隆，Molecular Science, 18, A0136 (2024).
• [20]  P. Atkins, J. De Paula and J. Keeler, Atkins’ Physical Chemistry, 12th Ed. Oxford University Press, London (2022).
• [21]  G. M. Barrow, Physical Chemistry, 6th Ed. (International Student Edition), McGraw-Hill, New York (1996).
• [22]  W. J. Moore, Basic Physical Chemistry, Prentice-Hall, Englewood Cliffs (1983).
• [23]   J. I.  Steinfeld,  J. S.  Francisco and  W. L.  Hase, Chemical Kinetics and Dynamics, 2nd Ed., Prentice Hall, New York (1998).
• [24]   P.  Gaspard and  M.  Nagaoka, J. Chem. Phys., 111, 5668 (1999).
• [25]   S.  Glasstone,  K. J.  Laidler and  H.  Eyring, The Theory of Rate Processes: the Kinetics of Chemical Reactions, Viscosity, Diffusion and Electrochemical Phenomena, McGraw-Hill, New York (1941).
• [26]   P.  Hänggi,  P.  Talkner and  M.  Borkovec, Rev. Mod. Phys., 62, 251 (1990).
• [27]   W.  Feller, An Introduction to Probability Theory and Its Applications Vol. I, Wiley, New York (1957).
• [28]   W.  Feller, An Introduction to Probability Theory and Its Applications Vol. II, Wiley, New York (1966).
• [29]   N. G.  Van Kampen, Stochastic Processes in Physics and Chemistry, 3rd Ed., North-Holland (2007).
• [30]   L. E.  Reichl, A Modern Course in Statistical Physics, 4th Ed., Wiley-VCH (2016).
• [31]   K. A.  Fichthorn and  W. H.  Weinberg, J. Chem. Phys., 95, 1090 (1991).
• [32]   D. W.  Heerman, Computer Simulation Methods in Theoretical Physics, 2nd Ed. Springer, Berlin (1990).
• [33]   N.  Metropolis,  A. W.  Rosenbluth,  M. N.  Rosenbluth,  A. H.  Teller and  E.  Teller, J. Chem. Phys., 21, 1087 (1953).
• [34]   R. P.  Feynman and  A. R.  Hibbs, Quantum Mechanics and Path Integrals, Emended ed., Dover, New York (2010).
• [35]   M.  Andersen,  C.  Panosetti and  K.  Reuter, Front Chem., 7, 1 (2019).
• [36]   M. F.  Weber and  E.  Frey, Rep Prog Phys., 80, 046601 (2017).
• [37]   N.  Takenaka,  H.  Sakai and  M.  Nagaoka, J. Phys. Chem. C, 118, 10874 (2014).
• [38]   N.  Takenaka,  H.  Sakai,  Y.  Suzuki,  P.  Uppula and  M.  Nagaoka, J. Phys. Chem. C, 119, 18046 (2015).
• [39]   A.  Bouibes,  N.  Takenaka,  S.  Saha and  M.  Nagaoka, J. Phys. Chem. Lett., 10, 5949 (2019).
• [40]   A.  Bouibes,  N.  Takenaka,  K.  Kubota,  S.  Komaba and  M.  Nagaoka, RSC Adv., 12, 971 (2022).
• [41]   A.  Bouibes,  N.  Sakaki and  M.  Nagaoka, ACS Omega, 8, 16570 (2023).
• [42]   H.  Suwa and  S.  Todo, Phys. Rev. Lett., 105, 120603 (2010).
• [43]   S.  Kirkpatrick,  C. D.  Gelatt Jr. and  M. P.  Vecchi, Science, 220, 671 (1983).
• [44]  宮下精二，熱・統計力学，培風館 (1993).
• [45]  2013年ノーベル化学賞「for the development of multiscale models for complex chemical systems」https://www.nobelprize.org/prizes/chemistry/2013/summary/

著者紹介

長岡 正隆（博士（工学））

〔経歴〕1988年京都大学大学院工学研究科博士課程修了，2002名古屋大学教授．2024年から現所属．〔専門〕分子シミュレーション，化学反応理論．〔趣味〕アート鑑賞．