この記事は現代化学(東京化学同人)の2022年4月号から2024年3月号まで連載された安藤耕司先生の「Pythonによる化学シミュレーション入門」を読んでやってみた記事です。詳しい内容は上記の「Pythonによる化学シミュレーション入門」をご参照ください。今回は、2023年1月号のブラウン運動とランジュバン方程式の内容です。
ランダムウォークからブラウン運動へ
最初に今回の内容を説明します。以前の記事でランダムウォークについて解説しました。ここでは、ランダムウォークで拡散現象をおおよそ再現できることを見てみました。ここで、ランダムウォークにおいては粒子の位置は乱数によってのみ決定されています。しかし、実際のブラウン運動では粒子は周囲の媒質から力を受けて運動しています。そこで、ランダムウォークでは乱数として考えていなかった粒子の運動を、ブラウン運動の物理的なモデルに置き換えて、実際の物理現象に拡張してみましょうというのが今回の内容です。
ランジュバン方程式
ランジュバン方程式とはブラウン運動をしている粒子に成り立つ次のような運動方程式です。
\displaylines{
m\frac{d^2x}{dt^2} = -\gamma\frac{dx}{dt} + R(t)
}
左辺は$F = ma = m\frac{d^2x}{dt^2}$です。つまり右辺がブラウン運動をしている粒子にかかる力であると言っています。右辺の第一項はブラウン運動をしている粒子が媒質から受ける摩擦力です。この摩擦力の大きさは粒子の運動速度($\frac{dx}{dt}$)に比例していて、その比例定数が摩擦係数$\gamma$です。さらに、力の向きは運動方向と逆向きなので値にマイナスがついています。
次に、右辺の第二項はブラウン運動をしている粒子が媒質の粒子にランダムに衝突されるときに受ける力です。これについて考えます。ブラウン運動をしている粒子の周りにはmol数のオーダーの媒質の粒子が存在し、それらが熱振動によって絶え間なくぶつかってきます。従って、$R(t)$は確率変数としましょう。式中に確率変数が入ってくるせいで、後々積分したりするのに難儀しますが、今は$R(t)$統計的に満たすべき条件をまとめます。この力の統計平均$ < R(t) > $は0として良いでしょう。
\displaylines{
< R(t) > = 0
}
また、各時間に受ける力に相関はないでしょう。さらに、温度$T$が大きくなると$R(t)$が大きくなると考えられるし、媒質の摩擦係数$\gamma$が大きいと$R(t)$が大きくなるでしょう。これを式で表すと以下の式になります。
\displaylines{
< R(t_1)R(t_2) > = 2k_BT\gamma\delta(t_1-t_2)
}
ここで、$< R(t_1)R(t_2) >$は時間相関関数と呼ばれ、$k_B$はボルツマン定数、$\delta(t_1-t_2)$はデルタ関数です。
ランジュバン方程式のシミュレーション
シミュレーションの定法としては、以前の記事で紹介したように、微分方程式を数列に書き換えてオイラー法で計算していきたいところです。しかしながら、確率変数を含む微分方程式では、確率変数$R(t)$が決められないので、そのままではオイラー法を実行することができません。このような確率微分方程式を数値解析する方法としてオイラー・丸山法というのがあります。
確率微分方程式であるランジュバン方程式をオイラー・丸山法によって導いた式を次に示します。後でこれに数値を入れてシミュレーションをします。
\displaylines{
x(t+h) = x(t) + v(t)h \\
v(t+h) = v(t) + \frac{1}{m}(-\gamma v(t) h + \sqrt{2k_BT\gamma h}N(0,1)) \\
}
ここで$N(\mu,\sigma^2)$は平均$\mu$,分散$\sigma^2$の確率分布$f(x)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp(-\frac{(x-\mu)^2}{2\sigma^2})$です。
オイラー・丸山法の概説
証明は先人の解説に頼り、導出を概要だけ示します。
オイラー法では$ v(t) = \frac{dx}{dt} \simeq \frac{ x(t + h) - x(t)}{h} $と近似するのでした。$x(t+h) = x(t) + v(t)h$はこれを式変形したものです。次にランジュバン方程式を、次にように$t$から$t+h$の範囲で積分したものを考えます。
\int^{t+h}_{t} m\frac{dv(t)}{dt} dt= \int^{t+h}_{t}(-\gamma v(t) + R(t)) dt
ここから$\int^{t+h}_{t}(-\gamma v(t))dt \simeq -\gamma v(t)h$を用いて次式を得ます。
v(t+h) - v(t) =\frac{1}{m}( -\gamma v(t) h + \int^{t+h}_{t} R(t) dt)
ここで、$ < \hat{R}(t_1)\hat{R}(t_2) > = \delta(t_1-t_2) $を満たす$\hat{R}(t)=\frac{1}{\sqrt{2k_BT\gamma}}R(t)$を考えます。$\hat{R}(t)$は白色雑音と呼ばれ、その積分$\int^{t+h}_{t} \hat{R}(t) dt$を$N(0,h)=\sqrt{h}N(0,1)$とできることが知られています。以上代入して、次式を得ます。
\displaylines{
v(t+h) = v(t) + \frac{1}{m}(-\gamma v(t) h + \sqrt{2k_BT\gamma h}N(0,1))
}
シミュレーション
nrun個の粒子の一次元のランジュバン方程式のシミュレーションを実装します。
import numpy as np
import matplotlib.pyplot as plt
def langevin(m, gamma, kBT, tmax, h, x0, v0):
t, x, v = 0, x0, v0
t_list, x_list, v_list = [], [], []
while t <= tmax:
t_list.append(t)
x_list.append(x)
v_list.append(v)
t += h
x += v*h
v += (-gamma*v*h + np.sqrt(2*gamma*kBT*h)*np.random.normal(0,1))/m
return t_list, x_list, v_list
nrun = 20
m, gamma, kBT = 1, 1, 1
tmax, h = 100, 0.1
x0, v0 = 0, np.sqrt(kBT/m)
fig = plt.figure()
plt.rcParams["axes.prop_cycle"] = plt.cycler("color", plt.get_cmap("cool")(np.linspace(0,1,nrun)))
prop_cycle = plt.rcParams["axes.prop_cycle"]()
ax1 = fig.add_subplot()
for i in range(nrun):
t_list, x_list, v_list = langevin(m, gamma, kBT, tmax, h, x0, v0)
color = next(prop_cycle)['color']
ax1.plot(t_list, x_list, color=color)
ax1.set_xlabel("time")
ax1.set_ylabel("position")
plt.show()
各パラメーターは、粒子数が$20$、$m, \gamma, k_B T$が$1, 1, 1$、
初期位置$x_0$, 初期速度$v_0$が$0, \sqrt{k_B/m}$としました。初期速度$v0$は、エネルギー等分配則$\frac{1}{2}mv^2=\frac{1}{2}k_B T$を用いてます。実行結果がこちらです。
分布関数$f(x)$と平均二乗変位$ < x^2 > $はランダムウォークと同じく、
\displaylines{
f(x)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp(-\frac{x^2}{2\sigma^2}) \\
<x^2> = \sigma^2 = 2Dt
}
です。ブラウン運動では拡散係数$D$に以下のアインシュタインの関係式が成り立ちます。
\displaylines{
D = \frac{k_B T}{\gamma}
}
さらにここで摩擦係数$\gamma$は溶媒の粘性率$\eta$、溶質の微粒子の半径を$a$とすると次の関係が知れらています。
\displaylines{
\gamma = 6 \pi \eta a
}
これらの関係を使って、物性値を変更したときにブラウン運動がどう変化するかプロットしてみます。
import numpy as np
import matplotlib.pyplot as plt
def langevin2d(density, eta, a, kBT, tmax, h):
gamma = 6*np.pi*eta*a
m = density*(4*np.pi*a**3)/3
t, x, y, vx, vy = 0, 0, 0, np.sqrt(kBT/m), np.sqrt(kBT/m)
t_list, x_list, y_list , vx_list, vy_list = [], [], [], [], []
while t <= tmax:
t_list.append(t)
x_list.append(x)
y_list.append(y)
vx_list.append(vx)
vy_list.append(vy)
t += h
x += vx*h
y += vy*h
vx += (-gamma*vx*h + np.sqrt(2*gamma*kBT*h)*np.random.normal(0,1))/m
vy += (-gamma*vy*h + np.sqrt(2*gamma*kBT*h)*np.random.normal(0,1))/m
return t_list, x_list, y_list , vx_list, vy_list
tmax, h = 10, 0.1
fig = plt.figure(figsize=(8,8))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)
nrun= 20
for j in range(nrun):
t_list1, x_list1, y_list1, vx_list1, vy_list1 = langevin2d(1, 1, 1, 5, tmax, h)
t_list2, x_list2, y_list2, vx_list2, vy_list2 = langevin2d(1, 1, 5, 5, tmax, h)
t_list3, x_list3, y_list3, vx_list3, vy_list3 = langevin2d(1, 1, 1, 1, tmax, h)
t_list4, x_list4, y_list4, vx_list4, vy_list4 = langevin2d(1, 1, 5, 1, tmax, h)
ax1.plot(x_list1, y_list1, color="red")
ax2.plot(x_list2, y_list2, color="orange")
ax3.plot(x_list3, y_list3, color="green")
ax4.plot(x_list4, y_list4, color="blue")
ax1.set_xlabel("x",fontsize=12)
ax1.set_ylabel("y",fontsize=12)
ax1.set_xlim([-5,5])
ax1.set_ylim([-5,5])
ax2.set_xlabel("x",fontsize=12)
ax2.set_ylabel("y",fontsize=12)
ax2.set_xlim([-5,5])
ax2.set_ylim([-5,5])
ax3.set_xlabel("x",fontsize=12)
ax3.set_ylabel("y",fontsize=12)
ax3.set_xlim([-5,5])
ax3.set_ylim([-5,5])
ax4.set_xlabel("x",fontsize=12)
ax4.set_ylabel("y",fontsize=12)
ax4.set_xlim([-5,5])
ax4.set_ylim([-5,5])
plt.show()
シミュレーションに用いるパラメタを物性値に近づけてみました。このプログラムで指定する値は粒子の密度($density$)、溶媒の粘性率($eta$)、粒子の半径($a$)、ボツルマン定数と温度の積($kBT$)と、時間($tmax, h$)です。ここで粒子の半径と温度を変更した時、ブラウン運動がどう変化するか確認するために、[$density, eta, a, kBT$]を[$1, 1, 1, 5$] (赤)、[$1, 1, 5, 5$] (黄)、[$1, 1, 1, 1$] (緑)、[$1, 1, 5, 1$] (青)と指定してプロットしてみました。
結果は、赤＞緑＞黄色＞青の順番に拡散がしやすいという傾向が見られました。これはブラウン運動する粒子が小さいほど、またその時の温度が高いほど早く拡散が進行するということに対応します。