吸着の圧力温度依存性を色々プロットしてみたよ

三次元データって、どうプロットすればわかりやすいの？

科学やってたら、二種類のパラメータ$x$、$y$を変化させて関数$f(x,y)$の変化を見たいって時があるんですが、これってどうすればわかりやすくグラフにできるんだろうか、ってなりますよね。そんな疑問を解決するべく色々pythonでプロットしてみました。例として使うデータはなんでもいいんですが、今回は気体吸着の式を使ってみます。温度、圧力を変化させて吸着量をプロットしてみましょう。

温度圧力可変の吸着式

固体表面への単層吸着の物理吸着の場合を考えます。この等温条件での被覆率の圧力依存性は名前がついていて、ラングミュアの吸着等温式と呼ばれています。今回は変化させるパラメタが二つ欲しいのでラングミュアの吸着等温式より複雑な式を使います。被覆率($\theta$)の温度圧力依存性は次のように表されると言われています¹。

$ \theta = \frac{b(T)p}{1 + b(T)p}$
$ b(T) = \frac{s_0 \exp(\frac{q}{RT})}{v \sigma \sqrt{2 \pi m RT}} $

ここで、$s_0$, $q$, $v$, $\sigma$, $m$ はそれぞれ、吸着剤に到達した吸着分子の吸着確率、吸着熱[J/mol]、周波数[Hz]、サイトの密度、吸着物質の分子量を表します。$ R = 8.31$[J/K mol]は気体定数です。

プログラム

さて、早速プロットを作ります。このプログラムでは6種類のプロットが得られます。

import numpy as np
import matplotlib.pyplot as plt


T = np.linspace(200,400,100)
P = 1013.25*np.exp(np.linspace(0.1,7,100))
R = 8.13
q = 25000
b0 = 1e-8

theta_matrix = np.zeros((len(T), len(P)))

for i in range(len(T)):
    bt = b0*(np.exp(q/(R*T[i]))/np.sqrt(R*T[i]))
    
    for j in range(len(P)):
        theta_matrix[i, j] = bt*P[j]/(1 + bt*P[j])

##plot isobar
fig_isothem = plt.figure()
ax_isothem = fig_isothem.add_subplot(111)

for k in range(0, len(T), len(T)//5):
    plt.plot(P/1000000,theta_matrix[k, :],label= str(round(T[k], 1)) + " K")
    
ax_isothem.set_xlabel('Pressuer [MPa]', fontsize=20, color='black')
ax_isothem.set_ylabel('Coverage', fontsize=20, color='black')

plt.legend(frameon=False)

fig_isothem.subplots_adjust(left=0.2)

##plot isobar
fig_isobar = plt.figure()
ax_isobar = fig_isobar.add_subplot(111)

for l in range(0, len(P), len(P)//5):
    plt.plot(T,theta_matrix[:, l],label="{:.3}".format(P[l]/1000000) + " MPa") 
    
ax_isobar.set_xlabel('Temperature[K]', fontsize=20, color='black')
ax_isobar.set_ylabel('Coverage', fontsize=20, color='black')

plt.legend(frameon=False)

fig_isobar.subplots_adjust(left=0.2)


T_matrix, P_matrix = np.meshgrid(T, P)

contourf_level = 20
contour_level = 50

#color map
fig_contourf, ax_contourf = plt.subplots()
contf = ax_contourf.contourf(P/1000000, T, theta_matrix, contourf_level, cmap = "gray", vmin=-0.1, vmax=1.1)
plt.xlabel("Pressuer[MPa]", fontsize=20, color='black')
plt.ylabel("Temperature[K]", fontsize=20, color='black')
plt.colorbar(contf, label="Coverage")
plt.subplots_adjust(left=0.15, right=0.95, bottom=0.15, top=0.95)

plt.show()

#contour
fig_contour = plt.figure()
ax_contour = fig_contour.add_subplot(111)
ax_contour = plt.contour(P/1000000, T, theta_matrix, cmap="gray", vmin=-2, vmax=2)
ax_contour.clabel(fmt='%1.1f', fontsize=12) 
plt.xlabel("Pressuer[MPa]", fontsize=20, color='black')
plt.ylabel("Temperature[K]", fontsize=20, color='black')
plt.colorbar(contf, label="Coverage")
plt.subplots_adjust(left=0.15, right=0.95, bottom=0.15, top=0.95)

plt.show()

#surface
fig_surfase, ax_surfase = plt.subplots(subplot_kw={"projection": "3d"})
surf = ax_surfase.plot_surface(P_matrix/1000000, T_matrix, theta_matrix.T, cmap = "rainbow", vmin=-0.1, vmax=1.1, antialiased=False)
ax_surfase.set_zlim(np.min(theta_matrix)- 0.1, np.max(theta_matrix) + 0.1)
ax_surfase.zaxis.set_major_formatter('{x:.02f}')

cax = fig_surfase.add_axes((0.9, 0.25, 0.025, 0.4))
fig_surfase.colorbar(surf, cax=cax)

ax_surfase.set_xlabel("Pressuer[MPa]", fontsize=11, color='black')
ax_surfase.set_ylabel("Temperature[K]", fontsize=11, color='black')
ax_surfase.set_zlabel("Coverage", fontsize=11, color='black')

ax_surfase.view_init(elev=25, azim=150)
plt.show()


#wire
fig_wire = plt.figure()

ax_wire = fig_wire.add_subplot(111, projection="3d")
ax_wire.plot_wireframe(P_matrix/1000000, T_matrix, theta_matrix.T, color = "darkblue")

ax_wire.set_zlim(np.min(theta_matrix) - 0.1, np.max(theta_matrix) + 0.1)

ax_wire.set_xlabel("Pressuer[MPa]", fontsize=11, color='black')
ax_wire.set_ylabel("Temperature[K]", fontsize=11, color='black')
ax_wire.set_zlabel("Coverage", fontsize=11, color='black')

ax_wire.view_init(elev=25, azim=150)
plt.show()

グラフ

上記のプログラムを動かして得られるプロットを確認します。

最初二つは二つのパラメーターのうち片方を固定してプロットをしています。

1. 吸着等温線

いくつかの温度において被覆率の圧力依存性をプロットしています。圧力が上がると吸着量が増加するのがわかりやすいです。グラフの注釈とそれぞれの曲線の比較から温度が増加すると吸着量が減少することがわかります。

1. 吸着等圧線

一つ目のプロットと逆で、いくつかの圧力で被覆率の温度依存性をプロットしています。先ほどとは逆で、温度が上がると吸着量が減少するのがわかりやすいです。グラフの注釈とそれぞれの曲線の比較から圧力が増加すると吸着量が増加することがわかります。

次の二つは二次元平面へのプロットです。

1. 濃淡プロット

黒い箇所が被覆率が小さく、白い箇所が被覆率が大きいです、圧力が高くて温度が低い時に吸着量が大きいことが一目瞭然です。あと二次元プロットにすると、極値がないことがわかりますね。

1. 等高線

等高線では、線上に被覆率の値が記載されているのが特徴でしょうか。値が入っていればどちらが大きいかすぐにわかりやすいです。

残りの二つは三次元的なプロットになります。

1. 表面プロット

温度と圧力の二次元の底面に対して、被覆率を高さと色で表しています。カラフルなんで見栄えが良いです。傾斜が変化量に対応するので、変化量も把握しやすいかと思います。プログラムでプロットした際にはドラッグで動かすことができるのでとてもわかりやすのですが、画像としてデータを残す際には二次元になってしまうので、どの方向からを描画するかによってわかりやすさが大きく変わります。プログラム中では以下の箇所で描画向きを調整できます。

ax_surfase.view_init(elev=25, azim=150)

1. ワイヤーフレーム

表面プロットが線画になったものです。線が密なところが変化量の大きいところに対応します。線なので、表面プロットで見えなかった奥が少し見えることもポイントかもしれません。

まとめ

同じデータを複数の方法でプロットしてみました。プロットによって同じデータでも伝わるイメージが異なるなというのが私の所感です。みなさんはどれが一番良いなどありましたでしょうか。三次元データだからと言って三次元にプロットすることが最善だという固定概念にとらわれず、どういう印象で伝えたいかによって使い分けができると良いのではないでしょうか。

1. Weiss, Werner, and Wolfgang Ranke. "Surface chemistry and catalysis on well-defined epitaxial iron-oxide layers." Progress in Surface Science 70.1-3 (2002): 1-151. ↩