1. 概要
LJポテンシャルをpythonで計算、可視化してみます。
2. LJポテンシャルの定義
詳しくはWikipediaを参照ください。
式は以下の通りです。
U(r)=4\varepsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right)
ここで、$\varepsilon$, $\sigma$ は原子の種類によって決まる定数です。異なる原子間でのポテンシャルを計算する場合にはそれぞれ相乗平均、相加平均をとります。
\varepsilon_{ij}=\sqrt{\varepsilon_i\varepsilon_j},\quad\sigma_{ij}=\frac{\sigma_i+\sigma_j}{2}
プログラムを書くときには
U(r)=\frac{A}{r^{12}}-\frac{B}{r^6}
というふうに定数A, Bを予め計算しておいてコードをクリーンにすることもあります。
空気中にもある程度存在する希ガスArで試してみます。Arのパラメーターは
\varepsilon=1.67\times10^{-21}[J],\quad\sigma=3.40\times10^{-10} [m]
でいきましょう。これでなくても良いのですが。
3. 実際に書いてみる
コード
eps = 1.67e-21
sigma = 3.40e-10
def U(r):
return 4 * eps * ((sigma / r) ** 12 - (sigma / r) ** 6)
print(U(1))
出力
-1.0319293498880002e-77
1m離れていると-77乗オーダーのエネルギーになるんですね。
大気圧での平均自由行程は68nm程度らしい(Wikipedia)ので、それだとどうでしょう。
コード
print(U(68e-9))
出力
-1.043749999999984e-34
-34乗くらいです。まあ単位がジュールなのでこんなもんでしょう。
4. 力
ポテンシャルから力を導きます。
\frac{d}{dr}U(r)=-F(r)
ですから、素直にrで微分して
\begin{align}
F(r)&=-4\varepsilon\left(-12\left(\frac{\sigma^{12}}{r^{13}}\right)+6\left(\frac{\sigma^6}{r^7}\right)\right)\\
&=24\varepsilon\left(2\left(\frac{\sigma^{12}}{r^{13}}\right)-\left(\frac{\sigma^6}{r^7}\right)\right)
\end{align}
を得ます。ちょっと汚くなりますね。
コード
def F(r):
return 24 * eps * (2 * (sigma ** 12 / r ** 13) - (sigma ** 6 / r ** 7))
print(F(68e-9))
出力
-9.209558823529126e-27
単位はN(ニュートン)です。
ここいらで定数を$A, B, C, D$とおいてしまって式を綺麗にしましょう。
\begin{align}
U(r)&=\frac{A}{r^{12}}-\frac{B}{r^6}\qquad\left(A=4\varepsilon\sigma^{12},B=4\varepsilon\sigma^{6}\right)\\
F(r)&=\frac{C}{r^{13}}-\frac{D}{r^7}\qquad\left(C=48\varepsilon\sigma^{12},D=24\varepsilon\sigma^{6}\right)\\
\end{align}\\
eps = 1.67e-21
sigma = 3.40e-10
A = 4 * eps * sigma ** 12
B = 4 * eps * sigma ** 6
C = 48 * eps * sigma ** 12
D = 24 * eps * sigma ** 6
def U(r):
return A / r ** 12 - B / r ** 6
def F(r):
return C / r ** 13 - D / r ** 7
可読性が上がりました！やったね！
5. 可視化
LJポテンシャルがどういう形をしているのか見てみましょう（ググったら出てきますが・・・）
数値計算ライブラリnumpy、グラフ描画ライブラリmatplotlibを使います。入っていない方はpip installしておいてください。
コード
import numpy as np
import matplotlib.pyplot as plt
r_arr = np.linspace(3e-10, 1e-9, 100)
plt.plot(r_arr, U(r_arr))
plt.show()
出力
馴染み深いLJポテンシャルが姿を現しました。
力の方も合わせて、ちょっと体裁を整えて可視化しましょう。
コード
r_min = 3e-10
r_max = 1e-9
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
r_arr = np.linspace(r_min, r_max, 100)
ax[0].plot(r_arr, U(r_arr), color='k')
ax[0].set_xlabel('r [m]')
ax[0].set_ylabel('Energy [J]')
ax[0].hlines(0, 0, r_max, color='r')
ax[0].set_xlim(r_min, r_max)
ax[0].set_ylim(-0.25e-20, 1e-20)
ax[0].set_title('Potential Energy')
ax[1].plot(r_arr, F(r_arr), color='k')
ax[1].set_xlabel('r [m]')
ax[1].set_ylabel('Force [N]')
ax[1].hlines(0, 0, r_max, color='r')
ax[1].set_xlim(r_min, r_max)
ax[1].set_ylim(-0.25e-10, 1e-10)
ax[1].set_title('Force')
plt.show()
出力
0.38nmくらいが安定な距離になるわけですね。
それより距離が短くなると爆発的にポテンシャルが高く、不安定になります。
MDで時間の刻み幅が大きすぎると、安定点を通り越して分子同士が近づきすぎてしまい、莫大な力が働いて分子が飛んでいってしまう事故がよくあります。
-1. 次回予告
次回は数値積分について学びます。
リンクはこちら
https://qiita.com/PlusF/items/22ed49534b6fc6be174c
前回の記事はこちらhttps://qiita.com/PlusF/items/48c98bae2437c28efc0e