複雑で大規模な不均一系における分子輸送を記述するための新規動的モンテカルロ法

概要

大規模不均一系における溶質分子の長時間にわたる大域的な拡散係数の解析は通常のMD計算では事実上不可能であるが，これを可能にする新たな粗視化の考え方と，そこで用いる新規動的MC法について，現実系への応用例とともに述べる．この動的MC法が生成する分子の軌跡は，位置に依存する拡散係数と自由エネルギーを有するn粒子拡散方程式を統計力学的に正しく満たし，MD計算と同様に，個々の分子軌跡から動的性質の詳細な解析が可能となる．

1　 はじめに

溶液などの均一媒体中における物質の拡散については，長い伝統の下に実験的，理論的に様々な取り組みがなされ，方法論的にもまた物理化学的描像としても十分に理解が進み，サイエンスとしてほぼ完成していると言えるだろう［1–4］．しかしながら，空間的に複雑な形状を持った不均一系における低分子の大域的な拡散は，その複雑さから十分な理解には至っておらず，解析手法についても理論，シミュレーションいずれも十分には確立されていない．

計算科学における不均一系の取り扱いの難しさは，ひとえに系の大きさと軌跡を追跡すべき時間の長さにある．つまり，ここで興味のある動的性質とは，物質が局所的な不均一性を数多く乗り越えて統計的に大きな空間を運動した後の大域的な拡散であり，この解析には非常に長い距離の運動を非常に長い時間をかけて直接追跡する必要がある．したがって，拡散係数が小さな物質に対しては通常のMD計算による解析は困難である．

たとえば，高分子電解質膜中でミクロに相分離している高分子相と水相の不均一性の空間スケールは3～5 nmである（図5，6参照）．このような不均一構造で自由エネルギー障壁を5個含む構造を作成するためには，1辺が20 nmの基本セル中に10⁶個オーダーの原子数が必要である．一方で，高分子電解質膜中の水素の拡散係数は実験的に10⁻⁶–10⁻⁷ cm²/sであると報告されており，これらの障壁を乗り越えてセルの端から端までの運動を追跡するためには，1–10 μsのシミュレーションが要求される．そして，このような個々の運動から拡散的性質を得るためには10⁴あるいはそれ以上の数の軌跡を統計的に平均して解析しなければならない．このような計算は，事実上不可能である．

幸いにも，我々が興味を持っている系は，高分子電解質膜中の物質輸送に代表されるように，輸送される小分子は空間を運動していくが，その環境を作っている高分子構造は大きくは変化しないといった系である．このことを利用して，小分子の運動は場の中の拡散運動に粗視化し，拡散方程式で近似する．

一方で，得られた拡散方程式を統計的に満たす分子軌跡を生成する新規動的MC法を開発し［5,6］，これを用いることにより，小分子に対してミリ秒もしくはそれ以上の長時間のMC計算を実現する．この計算から得られた軌跡の解析から，大域的な拡散係数が評価される．ここで得られる分子軌跡は，拡散的挙動という意味では物理的に正しいものであり，得られた分子軌跡に対してMD計算と同様な詳細な解析を施すことができる．

2　 粗視化の考え方

現実系における拡散分子に注目すると，分子は高分子など環境を構成する原子の局所的な運動により，常時ランダム力を受けていると考えることができる．これが駆動力となって，分子は環境が作る平均化されたポテンシャル場の中を拡散運動するとみなして，運動を拡散方程式で近似する．これは物理的に正当な近似であり，環境分子が拡散分子に及ぼす影響のすべてを，拡散方程式に含まれる拡散係数とポテンシャルに粗視化したことに相当する．

その際，たとえば高分子電解質膜中では，分子が高分子相にいるか水相にいるかによってランダム力は大きく異なる．また，同じ高分子相にいるとしても，Cの近傍にいるかFの近傍にいるか，またSやOの近傍にいるかによってもランダム力の性質は異なったものとなる．つまり，不均一系において駆動力，つまるところ拡散係数は大きく位置に依存する．これを$D_0(\mathit{r})$と表すこととする．ここで，$D_0(\mathit{r})$は以下に示すような自由エネルギー曲面の中を運動した結果としての拡散係数Dではなく，平坦な場の中を環境からの揺動力のみにより拡散した際の拡散係数である．そのため，Dとは区別してD₀と書く．

一方で，拡散分子は環境に応じて空間的に特定の分布を持つ．たとえば，高分子相と水相では存在確率が異なる．また，同じ高分子相でも，環境原子が大きくは動かないため，環境原子との引力，斥力相互作用により，これらを反映して異なる空間分布を持つ．つまり，この節での詳細な解析を待たなくても，不均一系においては拡散分子の持つ自由エネルギーは，明らかに位置に依存する．この位置に依存する自由エネルギーを$F(\mathit{r})$と書く．

これらにより，環境分子も含む極めて多数の原子に対するニュートンの運動方程式は，比較的少数の拡散分子の拡散方程式へと粗視化されたことになる．簡単な例として，拡散分子が1個の場合，1粒子拡散方程式は

  

$$\begin{array}{c}\hfill \frac{\partial P(\mathit{r},t)}{\partial t}=\{\frac{\partial }{\partial \mathit{r}}{D}_0\left(\mathit{r}\right)\frac{\partial }{\partial \mathit{r}}+\frac{\partial }{\partial \mathit{r}}{D}_0\left(\mathit{r}\right)\beta \frac{\partial F\left(\mathit{r}\right)}{\partial \mathit{r}}\}P(\mathit{r},t)\end{array}$$(1)

で与えられる．この拡散運動を分子軌跡として統計力学的に正しく追跡することができれば，その軌跡をMD計算同様に解析することができる．

2.1　 環境分子の粗視化

ここではより一般的に，系は高分子など極めて大きな数N個の原子${\mathit{R}}^N=({\mathit{R}}_1,{\mathit{R}}_2,\dots ,{\mathit{R}}_N)$で構成される環境と，${\mathit{r}}^n=({\mathit{r}}_1,{\mathit{r}}_2,\dots ,{\mathit{r}}_n)$で表される比較的少数のn個の拡散分子からなり［6］，また，分配関数は全相互作用$V({\mathit{r}}^n,{\mathit{R}}^N)$を用いて

  

$$\begin{array}{c}\hfill Z={\iint }_{}^{}\exp \left\{-\beta V({\mathit{r}}^n,{\mathit{R}}^N)\right\}d{\mathit{r}}^{n}d{\mathit{R}}^N\end{array}$$(2)

で表されるとする．ここで，全相互作用を拡散分子間，拡散分子－環境，環境分子間に分割して

  

$$\begin{array}{ccc}V({\mathit{r}}^n,{\mathit{R}}^N)\hfill & =\hfill & \sum _{i>j}^n\sum _j^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+\sum _i^n\sum _J^N{V}^{iJ}({\mathit{r}}_i,{\mathit{R}}_J)\hfill \\ & & +\sum _{I>J}^N\sum _J^N{V}^{IJ}({\mathit{R}}_I,{\mathit{R}}_J)\hfill \\ & =& \sum _{i>j}^n\sum _j^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+\sum _i^n{V}_{env}^i({\mathit{r}}_i;{\mathit{R}}^N)+V_{env}\left({\mathit{R}}^N\right)\end{array}$$(3)

とし，先に環境変数${\mathit{R}}^N$に関して積分を実行して

  

$$\begin{array}{c}\hfill Z={\int }_{}^{}z\exp \left\{-\beta \sum _{i>j}^n\sum _j^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)\right\}d{\mathit{r}}^n\end{array}$$(4)

と書くと，zは

  

$$\begin{array}{c}\hfill z={\int }_{}^{}\exp \left\{-\beta \sum _i^n{V}_{\mathit{\text{env}}}^i({\mathit{r}}_i;{\mathit{R}}^N)\right\}\exp \left\{-\beta V_{env}\left({\mathit{R}}^N\right)\right\}d{\mathit{R}}^N\end{array}$$(5)

であり，この中では拡散分子同士は相互作用していない．n = 1のとき

  

$$\begin{array}{c}\hfill z={\mathcal{⟨}\exp \left\{-\beta V_{env}^1({\mathit{r}}_1;{\mathit{R}}^N)\right\}\mathcal{⟩}}_{env}=\exp \left\{-\beta F_1^1\left({\mathit{r}}_1\right)\right\}\end{array}$$(6)

また，n = 2のとき

  

$$\begin{array}{ccc}z\hfill & =\hfill & {\mathcal{⟨}\exp \left\{-\beta (V_{env}^1({\mathit{r}}_1;{\mathit{R}}^N)+V_{env}^2({\mathit{r}}_2;{\mathit{R}}^N))\right\}\mathcal{⟩}}_{env}\hfill \\ & =\hfill & \exp \left\{-\beta F_2^{12}({\mathit{r}}_1,{\mathit{r}}_2)\right\}\hfill \end{array}$$(7)

である．$F_1^1\left({\mathit{r}}_1\right)$は拡散分子1の真空中から系中${\mathit{r}}_1$への移行のヘルムホルツ自由エネルギーであり，$F_2^{12}({\mathit{r}}_1,{\mathit{r}}_2)$は互いに相互作用しない拡散分子1と2を同時に真空中から系中${\mathit{r}}_1$，${\mathit{r}}_2$へと移行したときの自由エネルギー変化である．これらは，二体の平均力ポテンシャル${\varphi }^{12}({\mathit{r}}_1,{\mathit{r}}_2)$と

  

$$\begin{array}{c}\hfill {\varphi }^{12}({\mathit{r}}_1,{\mathit{r}}_2)=F_2^{12}({\mathit{r}}_1,{\mathit{r}}_2)-\{F_1^1\left({\mathit{r}}_1\right)+F_1^2\left({\mathit{r}}_2\right)\}+V^{12}({\mathit{r}}_1,{\mathit{r}}_2)\end{array}$$(8)

の関係にある．n = nのときも同様に

  

$$\begin{array}{c}\hfill z={\mathcal{⟨}\prod _i^n\exp \left\{-\beta V_{env}^i({\mathit{r}}_i;{\mathit{R}}^N)\right\}\mathcal{⟩}}_{env}=\exp \left\{-\beta F_n\left({\mathit{r}}^n\right)\right\}\end{array}$$(9)

と表される．

ここで，環境は高分子や無機固体，カーボン，MOFなどsolidなものであり，拡散分子との弱い相互作用は環境に大きな構造変化をもたらさないと考えると

  

$$\begin{array}{c}\hfill F_n\left({\mathit{r}}^n\right)\approx \sum _i^n{F}_1^i\left({\mathit{r}}_i\right)\end{array}$$(10)

と近似することができる．これを用いると，我々が興味を持っている本来の分配関数(2)もしくは(4)は

  

$$\begin{array}{c}\hfill Z\approx \int \exp \left\{-\beta (\sum _{i>j}^n\sum _{}^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+\sum _i^n{F}_1^i\left({\mathit{r}}_i\right))\right\}d{\mathit{r}}^n\end{array}$$(11)

で表される．式(11)は，環境変数を先にintegrate outして粗視化すると，全系の相互作用$V_n\left({\mathit{r}}^n\right)$は

  

$$\begin{array}{c}\hfill V_n\left({\mathit{r}}^n\right)=\sum _{i>j}^n\sum _{}^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+\sum _i^n{F}_1^i\left({\mathit{r}}_i\right)\end{array}$$(12)

と表され，位置に依存する一体の自由エネルギーと比較的少数の拡散分子間の直接の相互作用の和で得られることを示している．この計算負荷は本来の全原子MDと比較すると格段に小さく，手元のPCで充分間に合うものであり，極めて効率の良い粗視化が実現されている．

これまでに用いた近似は式(10)と拡散分子の運動は拡散方程式に従うということの二つだけであり，多くの粗視化モデルで用いられているような大胆な近似は含まれていない．また，式(10)が成り立たない場合でも，$F_n\left({\mathit{r}}^n\right)$を二体の自由エネルギー$F_2^{ij}({\mathit{r}}_i,{\mathit{r}}_j)$の和で表すなど，より良い近似を模索することもできる．

これまでの粗視化に対する議論をまとめると，n粒子系に対して解くべきn粒子拡散方程式は

  

$$\begin{array}{c}\frac{\partial P_n({\mathit{r}}^n,t)}{\partial t}=\hfill \\ \sum _{i=1}^n\{\frac{\partial }{\partial {\mathit{r}}_i}{D}_0^i\left({\mathit{r}}_i\right)\frac{\partial }{\partial {\mathit{r}}_i}+\frac{\partial }{\partial {\mathit{r}}_i}{D}_0^i\left({\mathit{r}}_i\right)\beta \frac{\partial (\sum _{j\ne i}^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+F_1^i\left({\mathit{r}}_i\right)}{\partial {\mathit{r}}_i}\}\hfill \\ P_n({\mathit{r}}^n,t)\hfill \end{array}$$(13)

である．

2.2　 パラメータの計算

拡散係数$D_0^i\left({\mathit{r}}_i\right)$と自由エネルギー$F_1^i\left({\mathit{r}}_i\right)$が得られれば，式(13)を基本方程式として大域的拡散係数の解析が可能となる．しかしながら，現実系の$F_1(\mathit{r})$と$D_0(\mathit{r})$を実験的に求める手法はなく，また理論もない．したがって，これらはMD計算に基づいて“実測”するしかない．

この節ではこれらのパラメータがMD計算からどのように求められるかについて簡単に議論しておく．

 2.2.1　$F_1(\mathit{r})$の計算

1個の分子を真空中から系の中の位置rへ移動させる移行の自由エネルギー$F_1(\mathit{r})$は，環境系の全空間にわたる三次元マップを得るという大量の計算を前提にすると，Widomの式に基づく粒子挿入法で，より高い精度が必要な場合はoverlapping distribution methodを用いて計算することになる［7］．位置依存性は，分子の挿入先や削除する分子の位置をrに指定することにより得られる．しかしながら，位置依存性を得るにはやはり相当に大きなサンプリングを実施しなければならず，一般的に計算コストは高い．

このため，我々のグループでは，Widomの式に基づいて粒子挿入法と同等なqualityでありながら圧倒的に高いサンプリング効率を持ち，計算コストの格段に低い新たな方法論を開発中である［8］．この方法に従うと，本稿で議論している拡散分子の解析に必要な分解能を超えた，はるかに精密な三次元自由エネルギー図の作成も可能である．

 2.2.2　${\mathit{D}}_0(\mathit{r})$の計算

位置に依存した拡散係数$D(\mathit{r})$は，従来溶解拡散法に基づいた浸透係数の計算の一環として発達してきており，代表的な方法にBerendsenの方法［9］，Woolf-Rouxの方法［10］がある．前者はSHAKEを用いて拡散分子をrにデルタ関数的に拘束した上で，環境から分子に作用する力の相関関数を求めて拡散係数を計算する．後者は遷移状態における反応速度の計算のために開発されたStraub-Berneの方法［11］の拡散係数版であり，拡散分子を調和振動子的にrに拘束し，分子の変位の相関関数を求めて拡散係数を得るものである．ここで必要としている揺動力のみからの拡散係数$D_0(\mathit{r})$は，拡散方程式の移流項である$F_1(\mathit{r})$を通常のポテンシャルから差し引いた上でMD計算を実行すれば得られる．

しかしながら，これらの方法では拡散分子が大きく，かつ溶液中で小さな溶媒と充分頻繁に衝突しているような場合には高い精度で拡散係数が見積もられるが，本稿で興味の対象としているような系，つまり拡散分子が小さくかつ環境分子があまり位置を変えず，拡散分子から衝突していかなければ力が働かないような場合には，値を大きく過小評価する．

このため我々は，環境の中で拡散分子を位置rに設定したメッシュ状のwall内に拘束しながら，wall内では環境と相互作用しながら自由に運動させ，得られた軌跡から拡散係数を評価するflat-bottom法と呼んでいる方法を開発した［12］．軌跡からは非直線的なMSDが得られるが，あらかじめ求めておいた形状因子を用いて補正する．この方法によると拡散分子の位置分解能は低くなるものの，拡散係数は十分精度よく求まる．この場合，MSDの評価には一定の空間が必要なので自動的に格子モデルとなる．

3　 $F_1\mathbf{(}\mathit{r}\mathbf{)}$と$D_0\mathbf{(}\mathit{r}\mathbf{)}$からの拡散運動の解析

ひとたび$F_1(\mathit{r})$と$D_0(\mathit{r})$が得られると，小分子のダイナミックスはMC法に限らず，拡散方程式に基づいたいくつかの道筋により記述することができる．ここでは従来用いられてきた代表的な解析手法をいくつか紹介し，MC法と比較しながらその長所短所についても議論しておく．

3.1　 従来の解析について

最も簡単な記述は溶解拡散モデル［13］である．しかしながら，反応座標が1次元の場合には有効に解析できるが，多次元の場合には経路探索に大きな困難を伴い，事実上適用できない．また，Kramersモデル［2–4］などの近似解析解も反応座標が1次元の場合には有効であるが，多次元の解析解は無い．

さらに，この種の解析で現在最もしばしば用いられている方法は，主にタンパク質の構造変化の研究に適用されているMSMモデル［14,15］であろう．この方法は拡散方程式を直接扱う手法とは多少異なり，$F_1(\mathit{r})$と$D_0(\mathit{r})$ではなく，軌跡の解析から状態間の遷移確率行列を抽出し，それをベースに解析する．行列の対角化からは$F_1(\mathit{r})$と同等な平衡分布が得られ，さらには平均最速到達時間などの解析が行われる．解析は厳密であるが，行列の扱いに起因して大規模系への適用が困難である．また遷移確率行列の解析手法が限られており，得られる情報は少ない．

3.2　 分子軌跡の生成による詳細解析

$F_1(\mathit{r})$と$D_0(\mathit{r})$から個々の分子の長時間の大域的な軌跡を得ることができれば，これらは拡散に関するすべての物理情報を含み，軌跡に対する多様で厳密な解析は，我々の化学的，物理的直観と整合した豊富な議論を可能とする．そこで，$F_1(\mathit{r})$と$D_0(\mathit{r})$から分子軌跡を直接得る方法について検討する．

分子や系の運動を粗視化して記述する際に最も重要なことは，本研究の目的からするとそれが拡散方程式を満たすことである．個々の粒子の運動を統計的に取り扱ったときに，それが拡散方程式に従うものであれば拡散運動を正しく記述できていると考えるのが自然であり，得られる平衡分布も保証される．

 3.2.1　Langevin方程式

Langevin方程式の数値解，つまりBrownian dynamicsは，上記の条件をすべて満足する有効な伝統的方法である［3,4］．しかしながら，Brownian dynamicsは連続座標系には有効であるが，粗視化した系つまり格子モデルには適用できず，格子座標系を何らかの方法で連続座標系に変換する必要がある．また，微分方程式である運動方程式の数値解であるため，解の安定性にも注意が必要である．

 3.2.2　マスター方程式

位置に依存する拡散係数とポテンシャルを持つ拡散方程式の離散的な近似解から出発して，マスター方程式がSzaboらにより導かれている［16］．このマスター方程式を例えばGillespie algorithm［17］により解けば，系の拡散的な軌跡を得ることができる．この方法は，化学反応の解析によく用いられている．しかしながら，Gillespie algorithmはBrownian dynamicsとは逆に，格子座標系に対しては有効であるが，連続座標系に対しては計算コストが極めて高く，実用的でない．

 3.2.3　MC法

MC法は格子座標系に対しても，また連続座標系に対しても安定に精度よく解を得ることのできる優れた方法である［18,19］．周知のように，MC法は歴史的には配置積分の実行，つまりボルツマン分布を生成する平衡量の計算手法としてとらえられ，その動的な意味については少数の研究を除いて積極的には語られてこなかったというきらいがある．

その中で，Metropolis法［20］は最もよく用いられているMCアルゴリズムのひとつであり，特に分子の拡散係数が位置に依存しない場合は，その動的な意味についてもある程度積極的に議論されてきている［21,22］．研究対象の多くは互いに相互作用する凝集系であるが，自由エネルギーは位置に依存しない系である．このような場合について，厳密な証明は省きながらも，凝集系における粒子の拡散が研究されてきている．さらに，1粒子系の場合は，位置に依存するポテンシャル場の中でのMC計算は，位置に依存するポテンシャルを持つ拡散方程式を満足することが証明されている［23］．しかしながら，拡散係数も位置に依存するn粒子系に対する明示的な議論は，少なくとも我々が調べた範囲では見つかっていない．

4　 新規動的MC法

本稿では，以上の議論を踏まえながら，我々が最近開発した新規動的MC法［5,6］について，導出も含めて多少詳細に解説する．この方法論に従って得られた分子の軌跡は．$F_1(\mathit{r})$と$D_0(\mathit{r})$を持つn粒子拡散方程式(13)を統計的に満たす．

4.1　 位置に依存する拡散速度を記述する遷移確率

Metropolis法を含む従来のMC法について先に述べておこう．MC法におけるi番目の分子の${\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_{\mathit{i}}$の遷移確率は，一般的な形で

  

$$\begin{array}{c}\hfill p_{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_{\mathit{i}}}=\frac{1}{M_i}{p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}\end{array}$$(14)

と表される．ここで，M_iは分子iがr_iから遷移できる状態の数，${p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}$はMC法におけるacceptの確率であり，Metropolis法の場合は

  

$$\begin{array}{c}\hfill {p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=min(1,\frac{{\pi }_i({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{{\pi }_i\left({\mathit{r}}_i\right)})\end{array}$$(15)

が用いられている．${\pi }_i\left({\mathit{r}}_i\right)$は分子が${\mathit{r}}_i$に見出される確率で，式(11)，(12)に従うと

  

$$\begin{array}{cc}{\pi }_i({\mathit{r}}_i;{\mathit{r}}^{n-1})\hfill & \propto \exp \left[-\beta \{\sum _{j\ne i}^n{V}^{ij}({\mathit{r}}_i,{\mathit{r}}_j)+F_1^i\left({\mathit{r}}_i\right)\}\right]\hfill \\ & =\exp \left\{-\beta V^i({\mathit{r}}_i;{\mathit{r}}^{n-1})\right\}\hfill \end{array}$$(16)

と書ける．これらに対し式(14)を用いてMC計算を行うと，拡散分子のカノニカル分布が形成される．

これに対し，新規動的MC法では，単位時間当たりの分子の${\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_{\mathit{i}}$への遷移確率を

  

$$\begin{array}{c}\hfill p_{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_{\mathit{i}}}=\frac{1}{M_i}{d}_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i){p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}\end{array}$$(17)

と選ぶ．$d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)$は分子の拡散速度を制御するための因子で，この方法で中心的な役割を果たす．

  

$$\begin{array}{c}\hfill d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{r}_i)=f({\mathcal{D}}_i\left({\mathit{r}}_i\right),{\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i))\end{array}$$(18)

であり，${\mathcal{D}}_i\left({\mathit{r}}_i\right)$はすべての空間にわたる$D_0^i\left({\mathit{r}}_i\right)$の最大値$D_0^{\max }=\max \left(D_0^i\left({\mathit{r}}_i\right)\right)$を用いて換算した，無次元化した拡散係数である．

  

$$\begin{array}{c}\hfill {\mathcal{D}}_i\left({\mathit{r}}_i\right)=\frac{D_0^i\left({\mathit{r}}_i\right)}{D_0^{\max }}\mathrm{\text{where}}0\le {\mathcal{D}}_i\left({\mathit{r}}_i\right)\le 1\end{array}$$(19)

4.2　 Fokker-Planck方程式

Fokker-Planck方程式を得るために，式(17)にマスター方程式を適用する．導出の詳細は文献［5,6］を参照されたいが，次の5つの条件が満たされていれば，式(13)の拡散方程式を統計的に満たすことが証明される．

  

$$\begin{array}{c}\hfill \text{条件1}\pi \left({\mathit{r}}_i\right){p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_{\mathbf{i}}){p\prime }_{{{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i\to {\mathit{r}}_i}^{}}\end{array}$$(20)

  

$$\begin{array}{c}\hfill \text{条件2}{d}_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)=d_i({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i,{\mathit{r}}_i)\end{array}$$(21)

  

$$\begin{array}{c}\hfill \text{条件3}{p\prime }_{{{\mathit{r}}_i-\mathrm{\Delta }{\mathit{r}}_i\to {\mathit{r}}_i}^{}}={p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}+O\left(\mathrm{\Delta }{\mathit{r}}_i^2\right)\end{array}$$(22)

  

$$\begin{array}{c}\hfill \text{条件4}\underset{}{\lim }{p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=1\end{array}$$(23)

  

$$\begin{array}{c}\hfill \text{条件5}{d}_i({\mathit{r}}_i,{\mathit{r}}_i)={\mathcal{D}}_i\left({\mathit{r}}_i\right)\end{array}$$(24)

条件1は微視的詳細つり合い，条件2はdの可換性，条件3はp′のtranslational invariance，条件4，5はp′とd_iの連続性である．

MC計算1ステップあたりの実時間は

  

$$\begin{array}{c}\hfill {\mathrm{\Delta }t}_{MC}=\frac{{\mathrm{\Delta }\mathit{r}}_i^2}{M_i{D}_{0MD}^{\max }}\end{array}$$(25)

であり，これを用いて実時間でのMSDが計算され，拡散係数を評価することができる．ここで，$D_{0MD}^{\max }$はMD計算から実測した$D_0(\mathit{r})$に基づいたものである．

4.3　 5つの条件を満足する遷移確率

条件1から5を満足する${p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}$と$d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)$は数多くあるが，その中のいくつかは直ちに見出すことができる．たとえば${p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}$に対しては，従来用いられてきていたものとして

  

$$\begin{array}{c}\hfill {p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=\min (1,\frac{\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{\pi \left({\mathit{r}}_i\right)}):\mathrm{\text{Metropolis}}\end{array}$$(26)

  

$$\begin{array}{c}\hfill {p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=\sqrt{\frac{\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{\pi \left({\mathit{r}}_i\right)}}:\mathrm{\text{Szabo}}\end{array}$$(27)

  

$$\begin{array}{c}\hfill {p\prime }_{{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}^{}}=\frac{2\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{\pi \left({\mathit{r}}_i\right)+\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}:\mathrm{\text{Heat bath}}\end{array}$$(28)

が挙げられる．ここでは，これらをそのまま用いる．拡散を制御する$d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)$に対しては，直ちに

  

$$\begin{array}{c}\hfill d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)=\frac{{\mathcal{D}}_i\left({\mathit{r}}_i\right)+{\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathbf{\alpha }}_i)}{2}\end{array}$$(29)

  

$$\begin{array}{c}\hfill d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)=\sqrt{{\mathcal{D}}_i\left({\mathit{r}}_i\right){\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathbf{\alpha }}_i)}\end{array}$$(30)

  

$$\begin{array}{c}\hfill d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)=\frac{2{\mathcal{D}}_i\left({\mathit{r}}_{\mathbf{i}}\right){\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathbf{\alpha }}_i)}{{\mathcal{D}}_i\left({\mathit{r}}_i\right)+{\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathbf{\alpha }}_i)}\end{array}$$(31)

が思い浮かぶ．式(29)，(30)，(31)はそれぞれ代数平均，幾何平均，調和平均の形をしている．これらはすべて，課せられた連続性や対称性等に起因するものである．式(26)～(28)と式(29)～(31)の組み合わせで作られる9つのMC遷移確率を表1に示す．

表1 式(26)～(28)と式(29)～(31)の組み合わせから得られる9つのMC遷移確率．

これらは$O\left(\mathrm{\Delta }{\mathit{r}}^2\right)$までの精度で統計力学的にはすべて等価であり，正しく拡散方程式(13)を満たす分子軌跡を生成する．したがって，組み合わせのどれを選んでも構わない．ただし，式(27)と(28)を含む式は確率が1を超える場合があり，これに対しては特別な処理が必要である．この煩雑さを避けるために，問題の生じないMetropolis形式の式(26)を用いたものを選択するのが賢明である．$d_i({\mathit{r}}_i,{\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)$としてたとえば式(29)を選ぶと

  

$$\begin{array}{c}\hfill p_{{\mathit{r}}_i\to {\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i}=\frac{1}{M_i}\frac{{\mathcal{D}}_i\left({\mathit{r}}_i\right)+{\mathcal{D}}_i({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{2}\min (1,\frac{\pi ({\mathit{r}}_i+\mathrm{\Delta }{\mathit{r}}_i)}{\pi \left({\mathit{r}}_i\right)})\end{array}$$(32)

となる．

表1の遷移確率をKramers-Moyal展開に代入すると，拡散方程式(13)が直接得られるので，確かめて見られるとよい．

5　 新規方法論の検証と計算例

式(32)に対して，これらが拡散方程式(13)を統計的に満たす分子軌跡を生成することを，一次元モデル系ならびに二次元モデル系を用いて検証する．その上で，ここで提案している手法を高分子電解質膜に適用し，燃料である水素の大域的な拡散の解析を計算例として示す．

5.1　 検証1：周期場の中の一次元格子上の拡散

一次元格子系における一粒子ダイナミックスに対して，式(32)による動的MC計算と拡散方程式(13)の数値解を比較することにより，方法論の検証を行った［5］．モデルは，周期的な存在確率と周期的な拡散係数

  

$$\begin{array}{c}\hfill \pi (i)\propto \frac{1}{2a+1}[1+a\{\cos \left(\frac{2m\pi i}{i_0}\right)+1\}]\end{array}$$(33)

  

$$\begin{array}{c}\hfill \mathcal{D}\mathcal{(}i)=\frac{D(i)}{D_{\max }}=\frac{1}{2b+1}[1+b\{\cos \left(\frac{2n\pi i}{i_0}\right)+1\}]\end{array}$$(34)

の中を粒子が拡散していくものであり，初状態は正規分布した粒子とした．MC計算の粒子軌跡から求めたMSDに対して正規分布を重みにして平均し，拡散方程式の数値解である分布の分散の時間発展と比較した．

二つの型の遷移確率を用いたMC計算の結果と，拡散方程式の数値解の結果を図1に比較する．図から明らかなように，MC計算と拡散方程式の結果は，極めて良い一致を示している．

図1  動的MC計算と拡散方程式の数値解の比較．MC計算における拡散係数は，いずれも代数平均型を用いた．

5.2　 検証2：二次元連続座標上で互いに相互作用するn粒子ダイナミックス

低密度極限もしくは低濃度極限を仮定するときは一粒子拡散を考えれば十分であるが，拡散分子濃度が高い場合には拡散する分子間の相互作用も考慮する必要がある．ここでは，拡散粒子としては二酸化炭素に相当するLJ粒子を仮定し，二次元系に対して検討を行った［6］．

モデル系は，高分子を模したある不均一な構造を持つ固体粒子中のn個の拡散粒子であり，これらの粒子は互いに相互作用する．あらかじめMD計算を実施し，密度分布から求めたF(r)と位置に依存する拡散係数D₀(r)を用いて，n個の拡散粒子に対して動的MC計算を行い，その軌跡から50 nmを越える大域的な拡散係数を求めた．用いた密度分布を図2，拡散係数を図3に示す．いずれも強い不均一性を示しており，モデル系としては十分である．

図2  MC計算で仮定した参照系の密度分布．位置に依存する自由エネルギーの形では，exp{-βF(r)}に相当する．

図3  MC計算で仮定した位置に依存する拡散係数D₀(r)．一辺3.56 nmのメッシュ内では一定値を取る．

得られた大域的拡散係数の結果を，長時間計算のMD計算から求めた拡散係数の結果とともに図4に示す．n = 1は低密度極限，n = 1500は環境分子も含めると液体に近い密度に相当する．両者は低密度から高密度に至るまで良好な一致を示し，方法論の有効性を実証している．

図4  大域的拡散係数のMC計算とMD計算結果の比較．希薄密度から液体に相当する密度まで，極めて良好な一致を示している．

5.3　 計算例：高分子電解質膜中の水素拡散

はじめにの節でも述べたが，高分子電解質膜は高分子相と水相がミクロに相分離した複雑で不均一な構造を取っている．図5，6は，含水率10（ひとつのスルホン酸基に対して何個の水が存在しているかという指標）の電解質膜における高分子相と水相のスナップショットを示したものである．水相の構造は含水率に大きく依存するが，共通的にスピノーダル様の複雑な構造を有している．水素分子は，このような構造の中で相を越えて大域的に拡散していく．

図5  含水率10の高分子電解質膜．茶色が高分子相，水色が水相．二相はミクロに相分離している．

図6  含水率10の高分子電解質膜の水相だけを描いたもの．複雑なスピノーダル様の構造をしている．

ここでは，この系に動的MC法を適用して行った研究［24］について，本稿の主テーマである方法論に関連したトピックスを中心に紹介する．拡散機構やその物理化学的意味等，サイエンスサイドの議論については共同研究者の永井哲郎氏が本誌に美しい図とともにまとめて紹介しているので，そちらを参照されたい［25］．

図7に，MC計算で得られたMSDを示す．まず，横軸と縦軸の数値に注目していただきたい．横軸の単位はμs，また縦軸の単位はnm²である．つまり，このMSDは4 μs，90 nmにも及ぶ水素拡散を記述している．基本セルの長さを越えたMSDに特段の意味はないが，ここでは本方法論を用いると何ができるかを理解していただくために，あえて示した．これは，通常のMD計算からは到底不可能な時間領域，空間領域の解析がなされていることを表している．さらに，この計算は富岳のようなスーパーコンピュータではなく，普通のメニコアのPCで計算されたものであることにも注意していただきたい．

図7  MC計算で得られた拡散分子の軌跡から計算したMSD．横軸の単位はμs，縦軸はnm²．

図7は，本方法論の威力を如実に物語るものであると勝手に信じているが，研究としては長距離の軌跡が大きな統計量で生成され，サイエンス研究の準備ができた，ということだけである．サイエンス研究としては，永井氏の総説を参照していただきたい．

 謝辞

ここで紹介した方法論の研究は，定年前後から福岡大の永井哲郎氏，新潟大の吉森明氏との共同研究の一環として行ってきたものである．電解質膜など実在系への展開は，ほぼ永井氏に依存している．吉森氏には，我々が泥臭く直観的に導出してきた式を系統的な解析から美しく導出する手法を教示していただいた．これにより，当初はあまり考えていなかった大きな展開に結びついた．ここで，あらためて両者に感謝の意を表したい．

参考文献

• [1]   R.  Zwanzig, Nonequilibrium Statistical Mechanics, Oxford University Press (2001).
• [2]   A.  Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer and Reactions in Condensed Molecular Systems, Oxford University Press (2013).
• [3]   H.  Risken, The Fokker-Planck Equation: Methods of Solution and Applications, 2nd ed., Springer-Verlag (1996).
• [4]   N. G.  van Kampen, Stochastic Processes in Physics and Chemistry, 3rd ed., Elsevier (2007).
• [5]   T.  Nagai,  A.  Yoshimori and  S.  Okazaki, J. Chem. Phys., 156, 154506 (2022).
• [6]   S.  Okazaki, J. Chem. Phys., 160, 174111 (2024).
• [7]   T.  Nagai,  K.  Fujimoto and  S.  Okazaki, J. Chem. Phys., 156, 044507 (2022).
• [8]   藤本  和士， 湯  之也， 永井  哲郎， 岡崎  進，第38回分子シミュレーション討論会，318S，姫路 (2024).
• [9]   S.-J.  Marrink and  H. J. C.  Berendsen, J. Phys. Chem., 98, 4155 (1994).
• [10]   T. B.  Woolf and  B.  Roux, J. Am. Chem. Soc., 116, 5916 (1994).
• [11]   J. E.  Straub,  M.  Borkovec and  B. J.  Berne, J. Phys. Chem., 91, 4995 (1987).
• [12]   T.  Nagai,  S.  Tsurumaki,  R.  Urano,  K.  Fujimoto,  W.  Shinoda and  S.  Okazaki, J. Chem. Theory Comput., 16, 7239 (2020).
• [13]   E.  Awoonor-Williams and  C. N.  Rowley, Biochim. Biophys. Acta, 1858, 1672 (2016).
• [14]  An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation, edited by G. R. Bowman, V. S. Pande and F. Noé, Springer (2014).
• [15]   B. E.  Husic and  V. S.  Pande, J. Am. Chem. Soc., 140, 2386 (2018).
• [16]   D. J.  Bicout and  A.  Szabo, J. Chem. Phys., 109, 2325 (1998).
• [17]   D. T.  Gillespie, Annu. Rev. Phys. Chem., 58, 35 (2007).
• [18]  The Monte Carlo Method in Condensed Matter Physics, edited by K. Binder, Springer-Verlag (1992).
• [19]   M.  Newman and  G. T.  Barkema, Monte Carlo Methods in Statistical Physics, Clarendon Press (1999).
• [20]   N.  Metropolis,  A. W.  Rosenbluth,  M. N.  Rosenbluth,  A. H.  Teller and  E.  Teller, J. Chem. Phys., 21, 1087 (1953).
• [21]   D. L.  Ermak, J. Chem. Phys., 62, 4189 (1975).
• [22]   P.  Meakin,  H.  Metiu,  R. G.  Petschek and  D. J.  Scalapino, J. Chem. Phys., 79, 1948 (1983).
• [23]   K.  Kikuchi,  M.  Yoshida,  T.  Maekawa and  H.  Watanabe, Chem. Phys. Lett., 185, 335 (1991).
• [24]   T.  Nagai and  S.  Okazaki, J. Chem. Phys., 157, 054502 (2022).
• [25]   永井  哲郎， 吉川  信明， 陣内  亮典， 木村  将之， 岡崎  進，アンサンブル, 26, 281 (2024).

著者紹介

岡崎 進（工学博士）

［経歴］1982年京都大学大学院工学研究科博士課程修了，名古屋大学名誉教授，分子科学研究所名誉教授，2024年から現職．［専門］分子シミュレーション．［趣味］畑，山歩き，古典文学．