LoginSignup
0
1
@chemweb000

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

Last updated at Posted at 2024-05-26

この記事は現代化学(東京化学同人)の2022年4月号から2024年3月号まで連載された安藤耕司先生の「Pythonによる化学シミュレーション入門」を読んでやってみた記事です。詳しい内容は上記の「Pythonによる化学シミュレーション入門」をご参照ください。今回は、2022年8月号の反応速度論(3)、2022年9月号の反応速度論(4)の内容です。

目標

いろいろな反応速度をプロットしてみよう

単分子n次反応

A \longrightarrow B

のような反応で、反応速度$v$が次のように表されるもの。

v = k A(t)^n = - \frac{dA(t)}{dt} = \frac{dB(t)}{dt}

ここで$A(t)$、$B(t)$は時刻$t$でのAとBの濃度、$k$は反応速度定数。

・解析解

$k A(t)^n = - \frac{dA(t)}{dt}$より、この微分方程式を解いて、

\frac{1}{A(t)^{n-1}} = \frac{1}{A(0)^{n-1}} + (n-1)kt

import numpy as np
import matplotlib.pyplot as plt

k, t_max, h = 0.5, 5.0, 0.01
t, a, b = 0.0, 1.0, 0.0
n_list = [0, 1, 2, 3, 4]

ttn = np.arange(t, t_max, h)

for i in range(len(n_list)):
  if n_list[i] == 1:
    a_exact = a*np.exp(-k*ttn)
  else:
    a_exact = (a**(1-n_list[i]) + (n_list[i]-1)*k*ttn)**(1/(1-n_list[i]))

  a_exact = [0 if j <= 0 else j for j in a_exact]
  plt.plot(ttn, a_exact, label=f'[A], n={n_list[i]}')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

n_reaction.png

二分子二次反応

A + B \longrightarrow C

のような反応で、反応速度$v$が次のように表されるもの。

v = k A(t)B(t) = - \frac{dA(t)}{dt} = - \frac{dB(t)}{dt} = \frac{dC(t)}{dt}

ここで$A(t)$、$B(t)$、$C(t)$は時刻$t$でのA、B, Cの濃度、$k$は反応速度定数。

・解析解

$k(t)A(t)B(t) = - \frac{dA(t)}{dt}$を解きたい。$A(t)= A(0) - x(t)$とすると、$B(t) = B(0) - x(t)$、$\frac{dA(t)}{dt} = -\frac{dx(t)}{dt}$となる。これを、最初の式に代入すると、微分方程式に含まれる$t$の関数が$x(t)$だけになるので解くことができます。

\displaylines{
A(t) = \frac{A(0)-B(0)}{A(0) - B(0) \exp\{-(A(0)-B(0))kt\} } A(0) \\

B(t) = \frac{A(0)-B(0)}{A(0) \exp\{-(A(0)-B(0))kt\} - B(0) }B(0) \\

C(t) = C(0) + A(0) - A(t)
}

・数値解(オイラー法)
オイラー法では$\frac{d y(t)}{dt}=x$を$y(t+h) \approx y(t) + x*h$と近似して計算します。$A(t)$、$B(t)$、$C(t)$は次のようになります。

\displaylines{
A(t+h) = A(t) - v*h \\

B(t+h) = B(t) - v*h \\

C(t+h) = C(t) + v*h
}

import numpy as np
import matplotlib.pyplot as plt

k, t_max, h = 1, 5.0, 0.01
t, a, b, c = 0.0, 1.0, 1.1, 0.0

tt, aa, bb, cc = [], [], [], []

while t <= t_max:
  tt.append(t)
  aa.append(a)
  bb.append(b)
  cc.append(c)

  t += h
  v = k*a*b
  a -= v*h
  b -= v*h
  c += v*h

plt.plot(tt, aa, "--", label='[A]')
plt.plot(tt, bb, "--", label='[B]')
plt.plot(tt, cc, "--", label='[C]')

ttn = np.array(tt)

a_exact = (aa[0] - bb[0])/(aa[0] - bb[0]*np.exp(-(aa[0] - bb[0])*k*ttn)) * aa[0]
b_exact = (aa[0] - bb[0])/(aa[0]*np.exp((aa[0] - bb[0])*k*ttn) - bb[0]) * bb[0]
c_exact = cc[0] + aa[0] - a_exact

plt.plot(tt, a_exact, label='[A] exact')
plt.plot(tt, b_exact, label='[B] exact')
plt.plot(tt, c_exact, label='[C] exact')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

bimol_reaction.png

反応次数の決定 単離法

二分子の反応で反応の次数、$v = k A(t)^mB(t)^n$において、$m$と$n$を知りたいとします。$B(0)$を$A(0)$に対して大過剰にすると反応速度はAの単分子の反応に近い振る舞いをします。このような方法を単離法もしくは分離法と呼びます。

ここでは、二分子二次反応と単分子一次反応を比較します。

import numpy as np
import matplotlib.pyplot as plt

k, t_max, h = 0.1, 5.0, 0.01
t, a, b, c = 0.0, 1.0, [1.1, 2.0, 5.0, 10.0], 0.0

ttn = np.arange(t, t_max, h)

for i in range(len(b)):
  a_exact = (a - b[i])/(a - b[i]*np.exp(-(a - b[i])*k*ttn)) * a
  b_exact = (a - b[i])/(a*np.exp((a - b[i])*k*ttn) - b[i]) * b[i]
  c_exact = c + a - a_exact

  plt.plot(ttn, a_exact, label=f'[A] bimol, B(0) = {b[i]}')
  #plt.plot(ttn, b_exact, label=f'[B] bimol, B(0) = {b[i]}')
  #plt.plot(ttn, c_exact, label=f'[C] bimol, B(0) = {b[i]}')

a_unimol = a*np.exp(-k*ttn)

plt.plot(ttn, a_exact, "--", label='[A] unimol')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

mono_bi.png

実線のbimolと名前がついているのが二分子二次反応の$A(t)$です。$B(0)$を変化させてプロットしています。それに対して、unimolというのが単分子一次反応の$A(t)$です。確かに$B(0)$が大きくなると二分子二次反応の$A(t)$が単分子一次反応の$A(t)$と似通った時間変化を示すことがわかります。

自己触媒反応

R \longrightarrow P

のような反応で、反応速度$v$が次のように表されるもの。

v = k R(t)P(t) = - \frac{dR(t)}{dt} = \frac{dR(t)}{dt}

ここで$R(t)$、$P(t)$は時刻$t$での$R$と$P$の濃度、$k$は反応速度定数。

・解析解

二分子二次反応と同様に、$R(t)= R(0) - x(t)$とすると、$P(t) = P(0) + x(t)$、$\frac{dR(t)}{dt} = -\frac{dx(t)}{dt}$となるので、これらを反応速度の式に代入し微分方程式を解けば良いです。

\displaylines{
R(t) = R(0) - R(0)P(0)\frac{\exp(R(0) + P(0))kt - 1}{R(0) + P(0) \exp(R(0) + P(0))kt} \\

P(t) = P(0) + R(0)P(0)\frac{\exp(R(0) + P(0))kt - 1}{R(0) + P(0) \exp(R(0) + P(0))kt}

}

・数値解(オイラー法)

\displaylines{
R(t+h) = R(t) - v*h \\

P(t+h) = P(t) + v*h
}

import numpy as np
import matplotlib.pyplot as plt

k, t_max, h = 1.0, 10.0, 0.01
t, r, p = 0.0, 1.0, 0.01

tt, rr, pp = [], [], []

while t <= t_max:
  tt.append(t)
  rr.append(r)
  pp.append(p)

  t += h
  v = k*r*p
  r -= v*h
  p += v*h

plt.plot(tt, rr, label='[R]')
plt.plot(tt, pp, label='[P]')

ttn = np.array(tt)

r_exact = rr[0] - rr[0]*pp[0]*(np.exp((rr[0] + pp[0])*k*ttn) - 1)/(rr[0] + pp[0]*np.exp((rr[0] + pp[0])*k*ttn))
p_exact = pp[0] + rr[0]*pp[0]*(np.exp((rr[0] + pp[0])*k*ttn) - 1)/(rr[0] + pp[0]*np.exp((rr[0] + pp[0])*k*ttn))

plt.plot(tt, r_exact, "--", label='[R] exact')
plt.plot(tt, p_exact, "--", label='[P] exact')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

self.png

次に$P(0)$を変化させてプロットしてみます。

import numpy as np
import matplotlib.pyplot as plt

k, t_max, h = 1.0, 10.0, 0.01
t, r, p = 0.0, 1.0, [0.0, 0.001, 0.01, 0.1]

ttn = np.arange(t, t_max, h)

for i in range(len(p)):
  r_exact = r - r*p[i]*(np.exp((r + p[i])*k*ttn) - 1)/(r + p[i]*np.exp((r + p[i])*k*ttn))
  p_exact = p[i] + r*p[i]*(np.exp((r + p[i])*k*ttn) - 1)/(r + p[i]*np.exp((r + p[i])*k*ttn))

  plt.plot(ttn, r_exact, label=f'[R] exact, P(0)={p[i]}')
  plt.plot(ttn, p_exact, label=f'[P] exact, P(0)={p[i]}')


plt.legend(fontsize=14, loc='upper left', bbox_to_anchor=(1, 1))
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

self2.png

逐次反応と律速段階

A \longrightarrow B \longrightarrow C

のような反応を逐次反応と呼びます。2段階の反応の両者とも一次反応で、それぞれの速度定数を$k_1$、$k_2$とすると濃度の時間変化を次のように表すことできます。

\displaylines{
\frac{dA(t)}{dt} = - k_1A(t) \\

\frac{dB(t)}{dt} = k_1A(t) - k_2B(t) \\

\frac{dC(t)}{dt} = k_2B(t)
}

ここで$A(t)$、$B(t)$、$C(t)$は時刻$t$でのAとBとCの濃度。

・解析解
$A(t)$については、単分子一次反応と同様です。$B(t)$は$\frac{dB(t)}{dt} = k_1A(t) - k_2B(t)$の一階非斉次線形微分方程式を解けば良いです。$C(t)$は残りの微分方程式を解くともとまります。

\displaylines{
A(t) = A(0)exp(-k_1t) \\

B(t) = \frac{k_1}{k_2 - k_1}A(0)\{\exp(-k_1 t) - \exp(-k_2 t) \} + B(0) \exp(-k_2 t) \\

C(t) = A(0)\{1 + \frac{k_1 \exp(-k_2 t) - k_2 \exp(-k_1 t)}{k_2 - k_1} \}+ B(0)\{ \exp(-k_2 t) - 1 \} + C(0)
}

・数値解(オイラー法)

\displaylines{
A(t+h) = A(t) - k_1A(t)h \\

B(t+h) = B(t) + (k_1A(t) - k_2B(t))*h \\

C(t+h) = C(t) + k_2B(t)h
}

import numpy as np
import matplotlib.pyplot as plt

k1, k2, t_max, h = 1.0, 1.1, 10.0, 0.01
t, a, b, c = 0.0, 1.0, 0.0, 0.0

tt, aa, bb, cc = [], [], [], []

while t <= t_max:
  tt.append(t)
  aa.append(a)
  bb.append(b)
  cc.append(c)

  t += h
  a -= k1*a*h
  b += k1*a*h - k2*b*h
  c += k2*b*h

plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')
plt.plot(tt, cc, label='[C]')

ttn = np.array(tt)

a_exact = aa[0]*np.exp(-k1*ttn)
b_exact = k1/(k2 - k1)*aa[0]*(np.exp(-k1*ttn) - np.exp(-k2*ttn)) + bb[0]*np.exp(-k2*ttn)
c_exact = aa[0]*(1 + (k1*np.exp(-k2*ttn) - k2*np.exp(-k1*ttn))/(k2 - k1) ) + bb[0]*(1 - np.exp(-k2*ttn)) + cc[0]

plt.plot(tt, a_exact, "--", label='[A] exact')
plt.plot(tt, b_exact, "--", label='[B] exact')
plt.plot(tt, c_exact, "--", label='[C] exact')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

chiku.png

ここで、1段階目の反応が2段階目の反応より十分早い、つまり$k1 >> k2$を考えます。この条件では$\exp(-k_1 t)$が$\exp(-k_2 t)$より十分小さく無視でき、$C(t)$は$k_2$のみに依存します。反応速度の遅い方が律速段階です。これを確かめるために、$C(t)$の解析解と$\exp(-k_1 t)$を無視する近似を行った式をプロットしてみます。

import numpy as np
import matplotlib.pyplot as plt

k1, k2, t_max, h = [1.1, 5.0, 10.0], 1.0, 10.0, 0.01
t, a, b, c = 0.0, 1.0, 0.0, 0.0

ttn = np.arange(t, t_max, h)

for i in range(len(k1)):
  a_exact = a*np.exp(-k1[i]*ttn)
  b_exact = k1[i]/(k2 - k1[i])*a*(np.exp(-k1[i]*ttn) - np.exp(-k2*ttn)) + b*np.exp(-k2*ttn)
  c_exact = a*(1 + (k1[i]*np.exp(-k2*ttn) - k2*np.exp(-k1[i]*ttn))/(k2 - k1[i]) ) + b*(1 - np.exp(-k2*ttn)) + c

  #plt.plot(ttn, a_exact, label=f'[A] exact, k1/k2 = {k1[i]/k2}')
  #plt.plot(ttn, b_exact, label=f'[B] exact, k1/k2 = {k1[i]/k2}')
  plt.plot(ttn, c_exact, label=f'[C] exact, k1/k2 = {k1[i]/k2}')

c_approx = a*(1 - np.exp(-k2*ttn)) + b*(1 - np.exp(-k2*ttn)) + c

plt.plot(ttn, c_approx,"--", label=f'[C] approx.')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

rissoku.png

$[C] approx.$というのが$\exp(-k_1 t)$を無視した式です。解析解は$k_1$が大きいほどこれに近づいていくことがわかります。

ここで、逐次反応において解析解と数値解の差を比較しておきます。

for i in range(len(h)):
  t, a, b, c = 0.0, 1.0, 0.0, 0.0
  tt, aa, bb, cc = [], [], [], []

  while t <= t_max:
    tt.append(t)
    aa.append(a)
    bb.append(b)
    cc.append(c)

    t += h[i]
    a -= k1*a*h[i]
    b += k1*a*h[i] - k2*b*h[i]
    c += k2*b*h[i]

  ttn = np.array(tt)

  a_exact = aa[0]*np.exp(-k1*ttn)
  b_exact = k1/(k2 - k1)*aa[0]*(np.exp(-k1*ttn) - np.exp(-k2*ttn)) + bb[0]*np.exp(-k2*ttn)
  c_exact = aa[0]*(1 + (k1*np.exp(-k2*ttn) - k2*np.exp(-k1*ttn))/(k2 - k1) ) + bb[0]*(1 - np.exp(-k2*ttn)) + cc[0]

  aaa = aa - a_exact
  bbb = bb - b_exact
  ccc = cc - c_exact

  plt.plot(tt, aaa, label=f'diff. [A], h = {h[i]}')
  plt.plot(tt, bbb, label=f'diff. [B], h = {h[i]}')
  plt.plot(tt, ccc, label=f'diff. [C], h = {h[i]}')


plt.legend(fontsize=14, loc='upper left', bbox_to_anchor=(1, 1))
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

chiku_diff.png

$h$を変化させて解析解と数値解の差分をプロットしてみました。$A(t)$と$B(t)$は$h=0.01$で$t$を大きくしていくと差分が$0$に収束していく傾向が見られます。一方で$C(t)$は差分が$0$に収束しなさそうですし、$h=0.001$にしても$A(t)$と$B(t)$に比べ差分が大きいです。理論的に正しいのは解析解です。数値解はあくまで近似値なので、常にどれくらいの誤差が生じているのか確認することが重要です。

対向反応と平衡定数

一般的に$A$が$B$になる反応では、$A \rightarrow B$となる正反応だけでなく、$B \rightarrow A$の逆反応が同時に進行します。

A \rightleftharpoons B

逆反応も考慮すると、反応速度は次のように表されます。

\frac{dA(t)}{dt} = - k_fA(t) + k_bB(t)

ここで$A(t)$、$B(t)$、は時刻$t$でのAとBの濃度、$k_f$と$k_b$は正反応と逆反応の速度定数です。

十分時間が経過した平衡状態の反応物と生成物の比を平衡定数$K$と呼び、速度定数と次のような関係があります。

K = \frac{B(\infty)}{A(\infty)} = \frac{k_f}{k_b}

・解析解

$B(0) = 0$とする時$B(t) = A(0) - A(t)$となり、これを代入して反応速度の微分方程式を解く。

\displaylines{
A(t) = \frac{k_b + k_f \exp\{-(k_f + k_b)t\} }{k_f + k_b}A(0) \\

B(t) = A(0) - A(t)
}

・数値解(オイラー法)

\displaylines{
A(t+h) = A(t) - k_fA(t)h + k_bB(t)h \\

B(t+h) = B(t) + k_fA(t)h - k_bB(t)h \\
}

import numpy as np
import matplotlib.pyplot as plt

kf, kb, t_max, h = 2.0, 1.0, 10.0, 0.01
t, a, b = 0.0, 1.0, 0.0

tt, aa, bb = [], [], []

while t <= t_max:
  tt.append(t)
  aa.append(a)
  bb.append(b)

  t += h
  a += -kf*a*h + kb*b*h
  b += kf*a*h - kb*b*h
  
plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')

ttn = np.array(tt)

a_exact = aa[0]*(kb + kf*np.exp(-(kf + kb)*ttn))/(kf + kb)
b_exact = aa[0] - a_exact

plt.plot(tt, a_exact, "--", label='[A] exact')
plt.plot(tt, b_exact, "--", label='[B] exact')

plt.legend(fontsize=14, loc='upper left', bbox_to_anchor=(1, 1))
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.text(7, 0.05, f"t = {format(tt[-1], '.3f')} \nkf/kb = {format(kf/kb, '.3f')} \nB/A = {format(b_exact[-1]/a_exact[-1], '.3f')}", fontsize=16)

plt.show()

eq.png

十分時間が経過して平衡状態になった時のB(t)とA(t)の比を計算すると、$k_f$と$k_b$の比と一致していることがわかります。この値が平衡定数です。

前駆反応

A \rightleftharpoons B \longrightarrow C

逐次反応の1段階目の反応に逆反応があるパターン。

\displaylines{
\frac{dA(t)}{dt} = - k_fA(t) + k_bB(t) \\

\frac{dB(t)}{dt} = k_fA(t) - k_bB(t) - k_rB(t) \\

\frac{dC(t)}{dt} = k_rB(t)
}

ここで$A(t)$、$B(t)$、$C(t)$は時刻$t$でのAとBとCの濃度、$k_f$と$k_b$と$k_r$は1段階目の正反応と1段階目の逆反応の2段階目の速度定数です。

もう、反応が複雑なので、数値解のみ示します。

・数値解(オイラー法)

\displaylines{
dA(t+h) = - k_fA(t)h + k_bB(t)h \\

dB(t+h) = k_fA(t)h - k_bB(t)h - k_rB(t)h \\

dC(t+h) = k_rB(t)h
}

import numpy as np
import matplotlib.pyplot as plt

kf, kb, kr, t_max, h = 1.0, 1.0, 1.0, 10.0, 0.01
t, a, b, c = 0.0, 1.0, 0.0, 0.0

tt, aa, bb, cc = [], [], [], []

while t <= t_max:
  tt.append(t)
  aa.append(a)
  bb.append(b)
  cc.append(c)

  t += h
  a += - kf*a*h + kb*b*h
  b += kf*a*h - kb*b*h - kr*b*h
  c += kr*b*h

plt.plot(tt, aa, label='[A]')
plt.plot(tt, bb, label='[B]')
plt.plot(tt, cc, label='[C]')

plt.legend(fontsize=14)
plt.ylabel("consentration", fontsize=16)
plt.xlabel("time", fontsize=16)

plt.show()

zenku.png

まとめ

いろいろな反応の反応速度をプロットしてみました。オイラー法を使えば複雑な反応も簡単にシミュレーションできることがわかります。特に、AとBの平衡の反応は$B(0)=0$の条件でしか解析解を示してませんし、前駆反応に関しては条件付きでも解析解がわかりません。(自分はこれらの反応の微分方程式が解くことが可能なのかもわかりませんでした。)しかし、そんな反応であってもそれっぽいシミュレーションができるということがオイラー法の効果的な点です。一方で、逐次反応で示した通り、あくまで近似値なので、理論値からの誤差をなくすことはできません。解析解が求められるときは解析解を用いましょう。

連載記事目次

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

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1