化学反応の反応速度をプロットしてみたよ~Arrheniusの式とEyringの式の比較~

# 化学反応の反応速度

化学反応の反応速度について、プロットして遊んでみます。当初の目論みとしてはアレニウスの式とアイリングの式を比較して、違いを指摘しようと思ってたんですが、妥当な値の範囲では傾向に大きな違いはないなと思いました。

### Arrheniusの式
化学反応の反応速度について、よく知られている式にアレニウスの式があります。

```math
k = A \exp(-\frac{E_a}{RT})
```
ここで、$k$は反応速度定数、$A$は頻度因子、$E_a$(kJ/mol)は活性化エネルギー、$T$(K)は温度、$R$は気体定数です。

### Eyringの絶対反応速度論
さらに、化学反応の反応速度について、統計熱力学の知見も加味して導出されたアイリングの式というものがあります。

```math
k = \frac{k_bT}{h} \exp(-\frac{\Delta G}{RT})
```
ここで、$k_b$はボルツマン定数、$h$はプランク定数、$\Delta G$(kJ/mol)はギブスの活性化エネルギーです。


## plot

アレニウスの式とアイリングの式をプロットしてみます。反応速度定数の温度依存性と活性化エネルギー依存性を確認するため、k vs Tプロット、k vs 活性化エネルギープロット、それらの二次元プロットを描画します。それぞれの範囲は温度が$T=[250, 350]$、活性化エネルギーが$[50000, 150000]$です。温度は0℃が約273 Kなので概ね室温周辺です。また、有機反応の$\Delta G$は40から200 kJ/molくらいなので両者まあまあ妥当かと思われます。

```Python
import matplotlib.pyplot as plt
import numpy as np

kb = 1.38e-23
h = 6.62e-34
R = 8.31

t = np.linspace(250, 350, 101)
dg = np.linspace(50*1000, 150*1000, 101)

T, dG = np.meshgrid(t, dg)

ka = np.exp(-dG / (R*T))
k = (kb*T / h)*np.exp(-dG / (R*T)) 

contourf_level = 15

### Arrhenius

## ka vs T plot
fig_ta, ax_ta = plt.subplots()

plt.plot(t, ka[0, :], "-o", label=str(dg[0]/1000) + " kJ/mol")
plt.plot(t, ka[50, :], "-o", label=str(dg[50]/1000) + " kJ/mol")
plt.plot(t, ka[100, :], "-o", label=str(dg[100]/1000) + " kJ/mol")

plt.legend(frameon=False,fontsize=14)

plt.xlabel("T (K)", fontsize=20, color='black')
plt.ylabel("k (/s)", fontsize=20, color='black')

plt.subplots_adjust(left=0.15, right=0.95, bottom=0.15, top=0.95)
plt.show()

## ka vs dG plot
fig_ga, ax_ga = plt.subplots()

plt.plot(dg/1000, ka[:,0], "-o",  label=str(t[0]) + " K")
plt.plot(dg/1000, ka[:,50], "-o",  label=str(t[50]) + " K")
plt.plot(dg/1000, ka[:,100], "-o",  label=str(t[100]) + " K")

plt.legend(frameon=False,fontsize=14)

plt.xlabel(r"$E_a$ (kJ/mol)", fontsize=20, color='black')
plt.ylabel("k (/s)", fontsize=20, color='black')

plt.subplots_adjust(left=0.15, right=0.95, bottom=0.15, top=0.95)
plt.show()


## 2D ka plot
# color map
figa, axa = plt.subplots()
contf = axa.contourf(T, dG/1000, np.log10(ka), contourf_level)

#contour
cont = axa.contour(T, dG/1000, np.log10(ka), levels=contourf_level, colors="k")
axa.clabel(cont, fmt='%1.2g', fontsize=12)

plt.xlabel("T (K)", fontsize=20, color='black')
plt.ylabel(r"$E_a$ (kJ/mol)", fontsize=20, color='black')
plt.colorbar(contf, label="k (/s)")

plt.subplots_adjust(left=0.15, right=0.95, bottom=0.15, top=0.95)

plt.show()



### Eyring
## k vs T plot
fig_t, ax_t = plt.subplots()

plt.plot(t, k[0, :], "-o", label=str(dg[0]/1000) + " kJ/mol")
plt.plot(t, k[50, :], "-o", label=str(dg[50]/1000) + " kJ/mol")
plt.plot(t, k[100, :], "-o", label=str(dg[100]/1000) + " kJ/mol")

plt.legend(frameon=False,fontsize=14)

plt.xlabel("T (K)", fontsize=20, color='black')
plt.ylabel("k (/s)", fontsize=20, color='black')

plt.subplots_adjust(left=0.2, right=0.95, bottom=0.15, top=0.95)
plt.show()

## k vs dG plot
fig_g, ax_g = plt.subplots()

plt.plot(dg/1000, k[:,0], "-o",  label=str(t[0]) + " K")
plt.plot(dg/1000, k[:,50], "-o",  label=str(t[50]) + " K")
plt.plot(dg/1000, k[:,100], "-o",  label=str(t[100]) + " K")

plt.legend(frameon=False,fontsize=14)

plt.xlabel(r"$\Delta G ^{\ddag}$ (kJ/mol)", fontsize=20, color='black')
plt.ylabel("k (/s)", fontsize=20, color='black')

plt.subplots_adjust(left=0.2, right=0.95, bottom=0.15, top=0.95)
plt.show()


## 2D k plot
# color map
fig, ax = plt.subplots()
contf = ax.contourf(T, dG/1000, np.log10(k), contourf_level)

#contour
cont = ax.contour(T, dG/1000, np.log10(k), levels=contourf_level, colors="k")
ax.clabel(cont, fmt='%1.2g', fontsize=12)

plt.xlabel("T (K)", fontsize=20, color='black')
plt.ylabel(r"$\Delta G ^{\ddag}$ (kJ/mol)", fontsize=20, color='black')
plt.colorbar(contf, label="k (/s)")

plt.subplots_adjust(left=0.2, right=0.95, bottom=0.15, top=0.95)

plt.show()
```

## 結果

### Arrheniusの式
・k vs Tプロット
![Figure_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/49e76cbd-82ec-46ab-8d7d-d70cfd30bab3.png)

・k vs 活性化エネルギープロット
![Figure_2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/c279a17e-9318-4799-aec9-6f32f1f07803.png)

・二次元プロット
![Figure_3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/1f64d852-f42a-4a39-acca-1c194e61e74e.png)


### Eyringの式
・k vs Tプロット
![Figure_4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/97e82c7f-d1a7-4b09-8d05-cabee3a7cce7.png)

・k vs 活性化エネルギープロット
![Figure_5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/fa0c9098-2aeb-47b6-be50-60b164751c4e.png)

・二次元プロット
![Figure_6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/453d3e20-b336-473f-8ba9-34dc12cea09c.png)

## 比較
アレニウスの式の頻度因子を1としているものあり絶対値は大きく異なりますが、傾向としては、活性化エネルギーが小さくて温度が高いと反応速度が速いということが見てとれます。この点、アレニウスの式とアイリングの式に大きな相違はありません。アイリングの式は、アレニウスの式の頻度因子に温度が含まれる形となっていますが、指数関数部分の挙動が主に反応速度に反映されています。

kの温度依存性を見てみると、活性化エネルギーが小さいときは温度に依存して速度定数が大きく変化します。よく"温度が10度上がると反応速度はおよそ2倍になる"と言われているものここから来ていそうです。しかし、活性化エネルギーが大きくなると反応速度定数の温度変化は小さくなります。

## 追記
活性化エネルギーをとても小さくしてプロットしてみました。活性化エネルギーを$[10, 50]$の範囲でプロットした二次元プロットを示します。

### Arrheniusの式
![Figure_7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/57b731c5-1841-4b39-8e37-eda2216eb80b.png)


### Eyringの式
![Figure_8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2804113/3e3a0a04-5b21-418b-b220-a01ccaacfecf.png)


ここまで活性化エネルギーが小さいと、アイリングの式の指数関数の前についた温度の寄与が大きくなって、温度依存性が顕著になっていることがわかります。とはいえ、実際の有機反応の活性化エネルギーが40から200 kJ/molなのに対して、描画範囲が10から50 J/mol(*kJではなくJです)であり、桁が3桁違います。現実には、なかなかそういう化学反応はないとは思われますので、室温付近の一般的な有機反応の反応速度を扱う際には、アレニウスの式でもアイリングの式でもどちらでも良さそうです。