一般化ランジュバン方程式を用いた局所拡散係数の導出
本記事では，一般化ランジュバン方程式を出発点にすることにより，局所拡散係数が
$$
\begin{align*}
D(z)&=\lim_{s\rightarrow 0}\frac{-\tilde{C}_{v}(s;z)\langle z^2\rangle\langle\dot{z}^2\rangle}{\tilde{C}_{v}(s;z)\left[s\langle z^2\rangle+\langle\dot{z}^2\rangle/s\right]-\langle z^2\rangle\langle\dot{z}^2\rangle}
\end{align*}\tag{1}
$$
もしくは
$$
\begin{align*}
D(z)&=\frac{\langle z^2\rangle^2}{\int_0^{\infty}{\rm d}tC_{z}(t)}
\end{align*}\tag{2}
$$
で計算できることを解説します。
$${\tilde{C}_{v}(s;z)}$$は位置$${z}$$における速度自己相関関数$${C_{v}(t)}$$のラプラス変換，$${C_{z}(t)}$$は位置$${z}$$における位置自己相関関数を表します。
つまり，(1)式は速度自己相関関数を用いた表式，(2)式は位置自己相関関数を用いた表式ということになります。
局所拡散係数とは？
溶質の拡散係数は周囲の溶媒環境に応じて異なる値となります。
例えば，均一の水溶媒や有機溶媒下における拡散係数を求めたい場合は，溶質が系内のどこに位置していても同一の拡散係数となるため，局所的な値を考える必要がありませんし，通常のMD（Molecular Dynamics，分子動力学）計算結果から速度相関関数や平均二乗変位を解析することで拡散係数を求めることができます。
一方，生体分子シミュレーションでは非均一の溶媒環境における拡散係数に興味があることがあります。創薬の観点では脂質二重膜への透過性が重要な性質の一つに挙げられますが，水溶媒と脂質二重膜からなる非均一溶媒での局所拡散係数が透過性に関連します。
原理的には局所拡散係数も通常のMD計算で得られたtrajectory（軌跡）から解析できますが，効率が著しく悪いため工夫が求められます。
そこで，局所拡散係数を求めたい位置$${z_0}$$を平衡点に持つ調和ポテンシャル$${U(z, z_0)}$$を溶質の重心座標$${z}$$に課すことを考えます。
$$
\begin{align*}
U(z, z_0)&=\frac{k}{2}(z-z_0)^2
\end{align*}\tag{3}
$$
$${z}$$の挙動に着目した場合，溶媒と相互作用する調和振動子と見ますことができます。
一般化ランジュバン方程式とは？
一般化ランジュバン方程式は，着目した自由度以外の影響を摩擦力，ランダム力として記述することにより，着目した自由度の挙動を調べるための方程式です。
今回の問題設定に従うと(4)式になります。
$$
\begin{align*}
m\ddot{z}(t)&=-k(z(t)-z_0)-\int_0^t{\rm d}\tau\dot{z}(\tau)\zeta(t-\tau)+R(t)
\end{align*}\tag{4}
$$
$${\zeta(t-\tau)}$$は摩擦カーネル，もしくは記憶関数と呼ばれる摩擦を表す関数，及び$${R(t)}$$はランダム力です。
摩擦カーネルがデルタ関数の場合，通常のランジュバン方程式に帰着します。
以後，表示の簡略化するため，$${z(t)-z_0\rightarrow z(t)}$$とすることにします。
局所拡散係数の表式(1)の導出
(4)式から速度自己相関関数を用いた表式である(1)式を導くことを考えます。
まず，(4)式の$${\ddot{z}(t)，z(t)}$$を$${\dot{z}(t)}$$を用いた表式に変形します。
$$
\begin{align*}
m\frac{{\rm d}}{{\rm d}t}\dot{z}(t)&=-k\int_0^t{\rm d}\tau\dot{z}(t)-\int_0^t{\rm d}\tau\dot{z}(\tau)\zeta(t-\tau)+R(t)\\
&=-\int_0^t{\rm d}\tau\left[k+\zeta(t-\tau)\right]\dot{z}(t)+R(t)
\end{align*}\tag{5}
$$
次に(5)式の両辺に$${\dot{z}(0)}$$をかけます。
$$
\begin{align*}
m\frac{{\rm d}}{{\rm d}t}\dot{z}(0)\dot{z}(t)&=-\int_0^t{\rm d}\tau\left[k+\zeta(t-\tau)\right]\dot{z}(0)\dot{z}(t)+\dot{z}(0)R(t)\\
\end{align*}\tag{6}
$$
(6)式の両辺の統計平均を取ることで速度自己相関関数$${C_v(t)}$$が従う積分微分方程式が得られますが，$${\dot{z}(0)R(t)}$$の統計平均が0（$${\lang\dot{z}(0)R(t)\rang = 0}$$）となるため(7)式になります。
$$
\begin{align*}
m\frac{{\rm d}}{{\rm d}t}C_v(t)&=-\int_0^t{\rm d}\tau\left[k+\zeta(t-\tau)\right]C_v(t)\\
\end{align*}\tag{7}
$$
積分微分方程式のままでは扱いづらいため，ラプラス変換します。
$$
\begin{align*}
m\left[s\tilde{C}_v(s)-C_v(0)\right]&=-\left[\frac{k}{s}+\tilde{\zeta}(s\right]\tilde{C}_v(s)
\end{align*}\tag{8}
$$
$${\tilde{\zeta}(s)}$$が左辺に来るように(8)式を変形します。
$$
\begin{align*}
\tilde{\zeta}(s)&=-\frac{m\left[s\tilde{C}_v(s)-C_v(0)\right]+\frac{k}{s}\tilde{C}_v(s)}{\tilde{C}_v(s)}\\
\therefore \frac{k_{\rm B}T}{\tilde{\zeta}(s)}&=-\frac{k_{\rm B}T\tilde{C}_v(s)}{m\left[s\tilde{C}_v(s)-C_v(0)\right]+\frac{k}{s}\tilde{C}_v(s)}
\end{align*}\tag{9}
$$
参考文献１によりますと，拡散係数$${D}$$と摩擦カーネル$${\zeta(t)}$$にはストークス・アインシュタインの式より(10)式が成立するとのことです。
$$
\begin{align*}
D&=\lim_{s\rightarrow 0}\frac{k_{\rm B}T}{\tilde{\zeta}(s)}&=\frac{k_{\rm B}T}{\int_0^{\infty}{\rm d}t\zeta(t)}
\end{align*}\tag{10}
$$
個人的に調べた限りでは，(10)式が一般的に成立することを示した参考書や論文にたどり着くことができませんでしたが，正しいと仮定して話を進めることにします。
ちなみに調和振動子の平衡点は$${z_0}$$なので$${D=D(z_0)}$$と書いた方が良いかもしれませんが，以後も$${D}$$の表記で進めます。
$$
\begin{align*}
D&=\lim_{s\rightarrow 0}\frac{k_{\rm B}T}{\tilde{\zeta}(s)}\\
&=\lim_{s\rightarrow 0}\frac{-k_{\rm B}T\tilde{C}_v(s)}{m\left[s\tilde{C}_v(s)-C_v(0)\right]+\frac{k}{s}\tilde{C}_v(s)}
\end{align*}\tag{11}
$$
ここで，ボルツマン統計における調和振動子の期待値がそれぞれ
$$
\begin{align*}
\lang z^2\rang&=\frac{k_{\rm B}T}{k}\\
\lang \dot{z}^2\rang&=\frac{k_{\rm B}T}{m}\\
\end{align*}\tag{12}
$$
となることを思い出すと，(11)式は(1)式に帰着します。
局所拡散係数の表式(2)の導出
ここでは，(1)式を変形して(2)式が得られることを導きます。
$${C_v(t)}$$は$${C_{z}(t)}$$を時間で2階微分したものに等しいことを利用します。
$$
\begin{align*}
C_v(t)&=\lang\dot{z}(0)\dot{z}(t)\rang& \\
&=\lang\dot{z}(t')\dot{z}(t+t')\rang& \\
&=\lang\dot{z}(t')\dot{z}(t)\rang& (t+t'\rightarrow t)\\
&=\frac{\partial^2}{\partial t\partial t'}\lang z(t') z(t)\rang\\
&=\frac{\partial^2}{\partial t\partial t'}\lang z(0) z(t-t')\rang\\
&=-\frac{\partial^2}{\partial t^2}\lang z(0)z(t-t')\rang\\
&=-\frac{\partial^2}{\partial t^2}C_{z}(t)\\
\end{align*}\tag{13}
$$
(13)式の左辺と最後の右辺をラプラス変換すると，
$$
\begin{align*}
\tilde{C}_{v}(s)&=-s^2\tilde{C}_{z}(s)+sC_{z}(0)+C_{z}^{(1)}(0)\\
&=-s^2\tilde{C}_{z}(s)+s\lang z^2\rang
\end{align*}\tag{14}
$$
になります（$${C_{z}^{(1)}(0)=\lang z\dot{z}\rang=0}$$）。
(14)式を(1)式に代入すると，
$$
\begin{align*}
D&=\lim_{s\rightarrow 0}\frac{(s^2\tilde{C}_{z}(s)-s\lang z^2\rang)\langle z^2\rangle\langle\dot{z}^2\rangle}{-(s^2\tilde{C}_{z}(s)-s\lang z^2\rang)\left[s\langle z^2\rangle+\langle\dot{z}^2\rangle/s\right]-\langle z^2\rangle\langle\dot{z}^2\rangle}\\
&=\lim_{s\rightarrow 0}\frac{(s\tilde{C}_{z}(s)-\lang z^2\rang)\langle z^2\rangle\langle\dot{z}^2\rangle}{-\tilde{C}_{z}(s)\left[s^2\langle z^2\rangle+\langle\dot{z}^2\rangle\right]+s\lang z^2\rang^2}\\
&=\frac{\lang z^2\rang^2}{\tilde{C}_{z}(0)}\\
&=\frac{\langle z^2\rangle^2}{\int_0^{\infty}{\rm d}tC_{z}(t)}
\end{align*}\tag{15}
$$
となり，(2)式に帰着します。
実用上の使い勝手の違い
理論上は(1)式と(2)式は同じ結果を与えますが，数値計算上では異なる性質を持ちます。
参考文献１では(1)式の優位性をアピールしているように見えますが，それ以外の論文では(2)式（もしくはそれ以外の手法）を採用している場合が殆どです。(1)式はラプラス変換をしなければならない点，及び$${s=0}$$が特異点となっている関係で外挿から求めなければいけない点が厄介な印象です。
(2)式も時間積分をどう計算するかは意見が分かれており，かつ計算手法に対して値が敏感に変わったりする難しさがありますが，それでも(1)式よりかは扱いやすい印象があります。
Thermostatの影響
通常のMD計算では熱浴と接している状況を取り入れることで有限温度の状態を実現します。今回考えている場合においては，調和ポテンシャル以外の粒子が調和ポテンシャルに対する熱浴の役割を果たします。更にnose-hooverやlangevinといった熱浴を設定したMD計算を実行することも可能ですが，その場合は厳密にはMD計算上の時間発展は物理的な意味を失いうることに注意する必要があります。相関関数がその影響を受けて拡散係数の誤差が大きくなってしまう可能性があります。
参考文献
この記事が気に入ったらサポートをしてみませんか？