「Pythonによる化学シミュレーション入門」をやってみたよ4~ランダムウォークと拡散現象~

この記事は現代化学(東京化学同人)の2022年4月号から2024年3月号まで連載された安藤耕司先生の「Pythonによる化学シミュレーション入門」を読んでやってみた記事です。詳しい内容は上記の「Pythonによる化学シミュレーション入門」をご参照ください。今回は、2022年11月号のランダム・ウォークと拡散現象(1)、2022年12月号のランダム・ウォークと拡散現象(2)の内容です。

ランダムウォーク

ランダムウォークとは、ある値が時間に対して確率的に変化していくようなモデルのことです。ランダムウォークは熱振動により微粒子が溶媒中で示すブラウン運動、拡散現象、株の値動きなど幅広い現象のモデルとして考察されています。

早速二次元のランダムウォークを図示します。

import matplotlib.pyplot as plt
import random

nstep = 100

x = 0.0
xx = []

y = 0.0
yy = []

for i in range(nstep):
  xx.append(x)
  yy.append(y)
  x += random.uniform(-1, 1)
  y += random.uniform(-1, 1)

plt.plot(xx,yy)

plt.legend(fontsize=14)
plt.ylabel("y", fontsize=16)
plt.xlabel("x", fontsize=16)

plt.show()

アニメーションにしてしてみると運動の軌跡がわかりやすいです。アニメーションもプロットしてみます。アニメーションは環境によって描画方法が異なるのですが、今回はGoogle colabで動くようにしています。

import numpy as np
import matplotlib.pyplot as plt
import random
from matplotlib import animation, rc

fig = plt.figure()
nsteps = int(input("Enter nsteps: "))

x, y = 0.0, 0.0

x_list, y_list, ims = [], [], []

for i in range(nsteps):
  x_list.append(x)
  y_list.append(y)

  x += random.uniform(-1, 1)
  y += random.uniform(-1, 1)

  ims.append((plt.plot(x_list,y_list, "b")[0],))

ani = animation.ArtistAnimation(fig, ims,interval=100)

rc('animation', html='jshtml')
ani

拡散の統計的性質

上記のように一点のランダムウォークのシミュレーションは簡単に実行することができます。しかし、実際の拡散現象においては、アボガドロ数のオーダーの粒子がそれぞれランダムウォークのような運動をしていると考えられます。したがって、拡散現象を議論するためには多数のランダムウォークの統計的な性質を知る必要があります。

まずはnrun個のランダムウォークの図示をしてみます。

import numpy as np
import matplotlib.pyplot as plt
import random

fig = plt.figure()
nsteps, nrun = 100, 20

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(121)
ax2 = fig.add_subplot(122)

for j in range(nrun):
  x, y = 0.0, 0.0
  x_list, y_list = [], []

  for i in range(nsteps):
    x_list.append(x)
    y_list.append(y)

    x += random.uniform(-1, 1)
    y += random.uniform(-1, 1)

  color = next(prop_cycle)['color']
  ax1.plot(x_list, color=color)
  ax2.plot(x_list, y_list, color=color)

ax1.set_xlabel("step")
ax1.set_ylabel("position")
ax2.set_xlabel("x")
ax2.set_ylabel("y")

plt.show()

ランダムウォークの統計的性質として次の二つがあります。

(1)変位の二乗の平均(分散)が時間に比例する

これは、変位x、分散$\sigma$、時間$t$を用いて次式で表されます。

\displaylines{
<x^2> = \sigma^2 = 2Dt
}

ここで、$<x^2>$が変位の二乗の平均を表し、$D$が拡散係数です。

(2)粒子の移動距離(変位)の分布が正規分布となる

これは、変位$x$、分散$\sigma$を用いて次式で表されます。

\displaylines{

f(x) = \frac{1}{\sqrt{2 \pi \sigma}} \exp{(-\frac{x^2}{2 \sigma})}

}

以下、これらを図示してみていきます。

変位の二乗の平均の時間依存性とその粒子数変化

縦軸$<x^2>$で横軸$t$に相当するグラフを書いてみましょう。粒子数に対応するnrunを変化させるとnrunが大きい方が直線に近づきます。この図の傾きが$2D$になるので、$y=ax$の直線近似から$D=a/2$を求めます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq
import random

fig = plt.figure()
nsteps, nrun = 100, [100, 1000, 10000]
xsq = []

for k in range(len(nrun)):

  xsq_k = np.zeros(nsteps)

  for j in range(nrun[k]):
    x = 0.0
    x_list = []

    for i in range(nsteps):
      x_list.append(x)
      x += random.uniform(-1, 1)

    xsq_k += np.array(x_list)**2

  xsq.append(xsq_k)
  plt.plot(xsq_k/nrun[k], "o", label=f"{nrun[k]} steps")

#fitting
xx = np.arange(nsteps)
yy = xsq[-1]/nrun[-1]

def func1(param,x,y):
  residual = y - (param[0]*x)
  return residual

param1 = [0]
prm = leastsq(func1, param1, args=(xx, yy))

fit_plot_range = np.max(xx) - np.min(xx)

x_func = np.arange(np.min(xx), np.max(xx), fit_plot_range/100)
y_func = prm[0][0]*x_func

plt.plot(x_func, y_func, "r--", label="<x^2> = {:10.4f} t ".format(prm[0][0]))

plt.legend(fontsize=14)

plt.xlabel("step", fontsize=16)
plt.ylabel("mean square displacement", fontsize=16)

plt.show()

変位の分布とその時間変化

変位の分布は、経過時間に対応するstep数を依存して次式のようになります。

\displaylines{

f(x) = \frac{1}{\sqrt{4 \pi D t}} \exp{(-\frac{x^2}{4 D t})}

}

さっきの直線近似から$2D \simeq 0.33$たっだので、$D=0.165$を上の式に代入して、実際に変位の分布が上記の式になっているか図示してみます。

import numpy as np
import matplotlib.pyplot as plt
import random

fig = plt.figure()
nsteps, nrun = 100, 10000
plot_steps = [9, 39, 99]
num_bins = []

trj, xx = [], []

for j in range(nrun):
  x = 0.0
  x_list = []

  for i in range(nsteps):
    x_list.append(x)
    x += random.uniform(-1, 1)

  trj.append(x_list)

trj_n = np.array(trj)

for k in range(len(plot_steps)):
  num_bins_k = int((np.max(trj_n[:,plot_steps[k]]) - np.min(trj_n[:,plot_steps[k]])) / 0.5)
  xx.append(np.linspace(np.min(trj_n[:,plot_steps[k]]), np.max(trj_n[:,plot_steps[k]]), num_bins_k))
  num_bins.append(num_bins_k)

def fx(xx, steps):
  fx = 1/(4 * np.pi * 0.165 * steps) * np.exp(-xx**2/(4 * 0.165 * steps))
  int_fx = sum(fx)
  fx = fx/int_fx*10000
  return fx

plt.hist(trj_n[:,plot_steps[0]], bins=num_bins[0], color="orange", rwidth=0.9, alpha=1.0) #10step目の位置の分布
plt.plot(xx[0], fx(xx[0], plot_steps[0]), "r--") #10step目の分布の理論式
plt.hist(trj_n[:,plot_steps[1]], bins=num_bins[1], color="green", rwidth=0.7, alpha=0.75) #40step目の位置の分布
plt.plot(xx[1], fx(xx[1], plot_steps[1]), "g--") #40step目の分布の理論式
plt.hist(trj_n[:,plot_steps[2]], bins=num_bins[2], color="blue" , rwidth=0.5, alpha=0.5) #最終stepの位置の分布
plt.plot(xx[2], fx(xx[2], plot_steps[2]), "b--") #最終stepの分布の理論式

plt.legend(["10th steps theoretical", "40th step theoretical", "final step theoretical", "10th step histogram", "40th step histogram", "final step histogram"])

plt.xlabel("position", fontsize=16)
plt.ylabel("counts", fontsize=16)

plt.show()

最後に、それぞれの粒子がランダムウォークを行い拡散していく様子を二次元的にプロットしてアニメーションにしてみます。

import numpy as np
import matplotlib.pyplot as plt
import random
from matplotlib import animation, rc

fig = plt.figure()
nsteps, nrun= 500, 10000
plot_steps = np.arange(0, nsteps-1, 10)

xend, yend, ims= [], [], []

for k in range(nrun):

  x, y = 0.0, 0.0
  x_list, y_list = [], []

  for i in range(nsteps):
    x_list.append(x)
    y_list.append(y)

    x += random.uniform(-1, 1)
    y += random.uniform(-1, 1)

  xend.append(x_list)
  yend.append(y_list)

xend_n = np.array(xend)
yend_n = np.array(yend)

plt.xlim(np.min(xend_n), np.max(xend_n))
plt.ylim(np.min(yend_n), np.min(yend_n))

plot_bins = np.arange(np.min(xend_n), np.max(xend_n), 0.5)

for l in range(len(plot_steps)):
  # Store the image (artist) returned by imshow, not the entire tuple
  _, xedges, yedges, image = plt.hist2d(xend_n[:,plot_steps[l]], yend_n[:,plot_steps[l]], bins=[plot_bins, plot_bins], cmap=plt.cm.jet)
  plt.xlim(np.min(xend_n), np.max(xend_n))
  plt.ylim(np.min(yend_n), np.max(yend_n))
  image.set_clim(0,20)
  ims.append([image])  # Wrap the image in a list for ArtistAnimation

plt.colorbar()
ani = animation.ArtistAnimation(fig, ims, interval=250)

rc('animation', html='jshtml')
ani

まとめ

一つ一つの粒子がランダムに動くランダムウォークを多数集めると拡散現象のモデルとなることがわかりました。拡散の様子はランダムウォークの統計的性質に従っていて、拡散係数がわかるとある時刻での拡散粒子の分布が予測できそうです。

連載記事目次

・「Pythonによる化学シミュレーション入門」をやってみたよ~単分子一次反応の反応速度論~

・「Pythonによる化学シミュレーション入門」をやってみたよ2~色々な反応の反応速度~

・「Pythonによる化学シミュレーション入門」をやってみたよ3~振動反応~