はじめに
身の周りに見られる現象は，突き詰めれば原子と原子の相互作用を通じて生じています。料理は素材の熱的変化と添加物質構成原子の侵入現象です。コップを叩くと鳴る音はコップの表面の振動が空気原子を揺らし，空気原子が波の変動パターンを作り，それが私たちの鼓膜にとどき，つぎに鼓膜を揺らします。
分子の動きをシミュレートできれば，原理的には多くのことを知ることができるでしょう。ここではPythonによって分子の運動をシミュレートする((分子動力学(MD)シミュレーション)ための簡易プログラムを作成します。
ここで扱うのは小さな計算セルサイズによるシミュレーション(NVE集団をとりあつかいます)ですが，MDシミュレーションで必要となる幾つかの大切なエッセンス(初期速度分布の設定，ビリアル定理による圧力計算方法，力のカットオフなど)が含まれいます。それらは大規模シミュレーションを行う際にも活かせると思います。また，ここに掲載しているのは2次元のシミュレーションコードですが，3次元への拡張は容易です。
内容
2次元原子系のMDシミュレーションを行い，粒子の時間発展を計算しアニメーションを作成します。
加えて，エネルギーや圧力といった各種物理量の時間発展の様子を調べます。
シミュレーションではレナード-ジョーンズポテンシャルによる2体相互作用を考慮します。
シミュレーション結果の後に各種テクニックをまとめています。
どんな人向け？
下記の方たちにとって少なからず有益だと思います。
● MDシミュレーションのエッセンスを手っ取り早くしりたい人
● Pythonを利用したMDコードのプロトタイプを作りたい人
計算コード
時間間隔: delta_t = 0.002
原子数: 25原子
格子は5x5の格子をとり，1辺の長さを1.291とした。
初期温度T_initial = 100 K, 質量m=1としています。
"""
レナード-ジョーンズ原子の2次元分子動力学シミュレーション
Sept. 2017
"""
%matplotlib nbagg
import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.animation as animation
fig = plt.figure(figsize = (8, 6))
ims = []
lat_x = 1.291 # 格子定数: x
lat_y = 1.291
Lx = 5 # 格子点: x
Ly = 5
N_atom = Lx*Ly
m=1 # 質量
delta_t = 0.002 # MDステップの間隔
T_initial = 50 # initial temperature [K]
T_initial = T_initial /119 # convert to natural unit
Max_step = 1000 #MD ステップ数
x = np.zeros([N_atom],float)
y = np.zeros([N_atom],float)
vx = np.zeros([N_atom],float)
vy = np.zeros([N_atom],float)
fx = np.zeros([N_atom,2],float)# Q: What is this "2"?
fy = np.zeros([N_atom,2],float)
# step毎の物理量の格納
T_lis = []
P_lis = []
KE_lis =[]
PE_lis =[]
Eot_lis =[]
Tav_lis = []
Pav_lis = []
KEav_lis =[]
PEav_lis =[]
Eotav_lis =[]
##
density = N_atom/(Lx*Ly*lat_x*lat_y) # 数密度: 2次元
print ("number density = ", density)
def twelveran():
s = 0
for i in range(1,13):
s += random.random()
return s/12.0 - 0.5
def Maxwell_velocity2D(TT): # マクスウェルの速度分布: ボックス-ミュラー法: self
for i in range(N_atom):
R1=random.random()
R2=random.random()
R3=random.random()
R4=random.random()
vx[i] = (np.sqrt(-2* (TT/m)*np.log(R1)))*np.cos(2*np.pi*R2)
vy[i] = (np.sqrt(-2* (TT/m)*np.log(R3)))*np.cos(2*np.pi*R4)
def initialposMaxwellvel():
i=-1
for ix in range(Lx):
for iy in range(Ly):
i += 1
x[i] = lat_x*ix
y[i] = lat_y*iy
Maxwell_velocity2D(T_initial)
def sign(a,b):
if b >= 0:
return abs(a)
else:
return -abs(a)
def Forces(t,w,PE,PEorW):
rcut = 4 # 力のカットオフ半径
PE = 0
r2cut = rcut**2
fx[:,t] = 0
fy[:,t] = 0
for i in range(N_atom-1):
for j in range(i+1, N_atom):
dx = x[i] - x[j]
dy = y[i] - y[j]
# 近接イメージ(鏡像)相互作用
if (abs(dx) > 0.5*lat_x*Lx):
dx = dx - sign(lat_x*Lx,dx)
if (abs(dy) > 0.5*lat_y*Ly):
dy = dy-sign(lat_y*Ly,dy)
r2 = dx**2+dy**2
if (r2 < r2cut):
if r2 == 0 :
r2 = 0.0001
invr2 = 1/r2
wij = 48*(invr2**3-0.5)*(invr2**3)
fijx = wij*invr2*dx
fijy = wij*invr2*dy
fx[i,t] += fijx
fy[i,t] += fijy
fx[j,t] -= fijx
fy[j,t] -= fijy
PE += 4*(invr2**3)*((invr2**3)-1)
w += wij
if (PEorW == 1) :
return PE
else :
return w
# メイン: 時間発展
def time_evolution():
avT = 0
avP = 0
Pavg = 0
avKE=0
avPE=0
t1 = 0
PE = 0.0
# 初期化
KE = 0
w = 0
initialposMaxwellvel() # 初期状態の構築
KE = (np.sum(vx**2)+np.sum(vy**2))/2
T = KE/N_atom
P = density*(2*KE+1.5*w)/(3*N_atom) # 圧力の計算(2次元系)
PE = Forces(t1, w, PE, 1)
time = 1
n_step=0
# 値の保存
T_lis.append(T*119) # K 単位
P_lis.append(P*2.382*1e-8) # Pa 単位
KE_lis.append(KE/N_atom)
PE_lis.append(PE/N_atom)
dumEtot = KE+PE
Eot_lis.append( dumEtot /N_atom)
while n_step < Max_step:
if n_step % 50 ==0:
print ("step=",n_step)
#print("T=",T*119,"K/N = ",KE/N_atom," Pot/N = ",PE/N_atom," Etot/N = ",(KE+PE)/N_atom, " P=",P)
#print("x[0] =",x[0] )
# ax.set_title("t =" + str("{0:.1f}".format(n_step*delta_t)))
if n_step % 5 ==0:
im=plt.scatter(x, y, s=120, c="red", alpha=0.9, linewidths="1.5",edgecolors="black")
ims.append([im])
n_step += 1
for i in range(N_atom):
PE = Forces(t1, w, PE, 1)
# 速度ベルレー法による更新
x[i] +=delta_t*(vx[i]+delta_t*fx[i,t1]/2)
y[i] +=delta_t*(vy[i]+delta_t*fy[i,t1]/2)
if x[i] <= 0:
x[i] += lat_x*Lx
if x[i] >= lat_x*Lx:
x[i] -= lat_x*Lx
if y[i] <= 0:
y[i] += lat_y*Ly
if y[i] >= lat_y*Ly:
y[i] -= lat_y*Ly
PE = 0.
t2=1
PE = Forces(t2, w, PE, 1)
KE = 0.
w = 0.
#for i in range(N_atom):
# 速度ベルレー法による更新
vx[:] += delta_t*(fx[:,t1]+fx[:,t2])/2
vy[:] += delta_t*(fy[:,t1]+fy[:,t2])/2
KE =(np.sum(vx**2)+np.sum(vy**2))/2
w = Forces(t2, w, PE, 2)
T = KE/N_atom
# P = density*(KE + w)
P = density*(2*KE+1.5*w)/(3*N_atom) # 圧力の計算
#平均値
avT += T
avP += P
avKE += KE
avPE += PE
time +=1
t = time
if (t == 0):
t = 1
# 時間平均
Pavg =avP/t
eKavg = avKE/t
ePavg = avPE/t
Tavg = avT/t
T_lis.append(T*119) # K 単位
P_lis.append(P*2.382*1e-8) # Pa 単位
KE_lis.append(KE/N_atom)
PE_lis.append(PE/N_atom)
dumEtot = KE+PE
Eot_lis.append( dumEtot /N_atom)
Tav_lis.append(Tavg*119 )# K 単位
Pav_lis.append(Pavg*2.382*1e-8)
KEav_lis.append(eKavg/N_atom)
PEav_lis.append(ePavg/N_atom)
##
time_evolution() #メイン
# プロット+動画の保存
plt.axis([0,lat_x*Lx,0,lat_y*Ly])
ani = animation.ArtistAnimation(fig, ims, interval =1)
plt.show()
ani.save("output.gif", writer="imagemagick")
結果
図. MDシミュレーション中の粒子の運動の様子。
MDシミュレーション中の運動エネルギー(KE)，相互作用ポテンシャルエネルギー(PE)，温度(T), 圧力(P)の変動を見てみましょう。
%matplotlib inline
import matplotlib.pyplot as plt
plt.close()
step=list(range(len(Eot_lis)))
plt.plot(step, KE_lis,'-', color='red', linewidth = 1, aa=True,label='KE/N')
plt.plot(step, PE_lis,'-', color='blue',linewidth = 1, aa=True,label='PE/N')
plt.plot(step, P_lis,'-', color='brown', linewidth = 1, aa=True, label='P')
#plt.plot(step, KEav_lis,'-', color='red', linewidth = 1, aa=True,label='KE/N')
#plt.plot(step, PEav_lis,'-', color='blue',linewidth = 1, aa=True,label='PE/N')
plt.plot(step, Eot_lis,'-', color='black', linewidth = 1, aa=True,label='Etot/N')
#plt.plot(step, Pav_lis,'-', color='brown', linewidth = 1, aa=True, label='P')
plt.legend(loc='best')
plt.xlabel('MD steps') # ｘ軸のラベル
plt.show()
plt.close()
plt.plot(step, T_lis,'-', color='purple', linewidth = 1, aa=True,label='T (K)')
#plt.plot(step, Tav_lis,'-', color='purple', linewidth = 1, aa=True)
plt.legend(loc='best')
plt.xlabel('MD steps') # ｘ軸のラベル
plt.show()
図. 粒子間相互作用を通じてKEとPEは常に変動しますがEtotは数値誤差内で一定です。
図. (即時)温度はMDシミュレーション中に変化していきます。
以上が，体積・エネルギー一定のMDシミュレーションの基本的な流れです。
シミュレーション手法
計算コードで使用した幾つかの手法について簡単に説明します。
速度ベルレ法による粒子の運動方程式
互いに力$F$を及ぼしあう粒子N個の座標($(x_i,y_i,z_i) \ i=1,2,...,N$)と速度($(vx_i,vy_i,vz_i)\ i=1,2,...,N$))に関する時間発展は，速度ベルレ法では以下のように計算します。
記述を簡単にするため，x座標だけを考えます。
$x_{n+1} = x_n +v_n \ \delta t +0.5 \ F(x_n,t_n)\ (\delta t)^2/M $
$v_{n+1} = v_n+ 0.5 \ \delta t\ (\ F(x_n, t_n) + F(x_{n+1}, t_{n+1}))/M$
実際にはx,y,zの3次元座標で計算します。粒子の数がN個なので，合計3N個の空間座標について計算することになります。また，速度も同様ですので，空間と速度をあわせた方程式の数は3N+3N = 6N個になります。
初期速度の決定
熱力学的に平衡状態にある系の速度分布は古典力学では"マクスウェルの速度分布"に従います。本シミュレーションではその分布にしたがった初期速度分布を用いています。
熱平衡状態において，粒子iの速度が$\bf v_i \sim v_i+dv_i$の範囲内にある確率を$f(v_i) \bf dv_i$は以下のように表されます。
(3次元の場合)
$$f(v_i)=(\frac{m}{2\pi k_BT})^{\frac{3}{2}}\exp\left(-\frac{mv_{ix}^2+mv_{iy}^2+mv_{iz}^2}{2kT}\right)$$
(2次元の場合)
$$f(v_i)=(\frac{m}{2\pi k_BT})\exp\left(-\frac{mv_{ix}^2+mv_{iy}^2}{2kT}\right)$$
$k_B$はボルツマン定数です。
これらはガウス分布関数の特殊形と見なすことができます。この分布関数を発生させるため本シミュレーションではボックス-ミュラー法を用います[3]。
[Pythonによる科学・技術計算] 与えられた確率密度関数を与える非一様乱数の生成，モンテカルロシミュレーション
この方法を用いると，マクスウェル速度分布に従う粒子の速度は，R1-R6を[0,1]に分布する一様乱数列として，
$$ v_x = (\frac{-2 k_B T}{m} ln(R_1))cos(2\pi R_2)$$
$$ v_y = (\frac{-2 k_B T}{m} ln(R_3))cos(2\pi R_4)$$
$$ v_z = (\frac{-2 k_B T}{m} ln(R_5))cos(2\pi R_6)$$
として生成できます。
力: レナード-ジョーンズポテンシャル
ここでは2体ポテンシャルとして中性原子系を良く記述することがしられている"レナード・ジョーンズポテンシャル"を用いて，$F = -\nabla U$により力を計算しています。
原子間距離を$r_{ij}$として，レナード・ジョーンズポテンシャルの関数形は以下のようになります。
$$U(r_{ij}) = 4β \left[ \left(\frac{\sigma}{r_{ij}}\right)^{12} - \left(\frac{\sigma}{r_{ij}}\right)^{6} \right] {\tag {13}}$$
$\beta$と$\sigma$は原子の種類に依存するパラメータです。
$r^{-12}$の項は核間距離が小さいときに生じる交換相互作用による強い反発力を，$r^{-6}$の項は双極子-双極子相互作用(ファンデルワールス相互作用)の主要項による引力を表していると考えることができます。
このとき，原子$j$が原子$i$に及ぼす力のベクトル$f_{ij}$は，
$$f_{ij}=\frac{48\beta}{r^2} \left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\frac{1}{2}\left(\frac{\sigma}{r_{ij}}\right)^{6} \right]\bf r_{ij} {\tag {14}} $$
となります。
圧力
圧力は，Vを体積として，
$$P= \frac{Nk_BT}{V}+\frac{1}{3V}\langle \ \sum_{i(\ne j)} \bf r_{ij} \cdot f_{ij} \ \rangle $$
として計算できます。
これは2体ポテンシャルに対するビリアル定理から導かれれます[2]。
周期境界条件
液体や固体はアボガドロ数($N_A \sim 6.02 \times 10^{23}/mol$ほどの原子から構成されいているのが普通です。その物質の性質を調べるには$N_A$個の原子を用いた大規模シミュレーションを行う必要がありますが，現実問題としてそれは不可能です。そのため，何らかの近似が必要です。
分子動力学シミュレーションで頻繁に使われる方法は，空間を一つのセルを周期的に並べ"表面を無くす"ことが行われます(下図)。シミュレーションセルが小さくなると表面効果が無視できなくなりますが，この方法によって表面を無くせば，実効的に"バルクの性質"をシミュレーションすることができます。これは境界条件として与えますので周期境界条件などともよばれます。
図. 周期的境界条件。真ん中のシミュレーションセル内の粒子の運動が周期的に繰り返されていると仮定する。
もし粒子が隣のセルに侵入したらどうしたらよいでしょうか。
その場合は原子を一周期分だけずらし，本物のシミュレーションセル内へと移動させます。
計算コード中の
if x[i] <= 0:
x[i] += Lx
if x[i] >= Lx:
x[i] -= Lx
if y[i] <= 0:
y[i] += Ly
if y[i] >= Ly:
y[i] -= Ly
がそれに対応します。
力の祐子範囲のカットオフ
上で周期境界条件を考えました。では力の計算はどうでしょうか。これは厄介です。
周期的に(無限に)セルが並ぶため，一つの粒子は無限個の粒子から力を受けます。力が減衰(収束)するにせよまともに計算することはできません。
そこでこの力の働く距離を有限の距離(カットオフ)で打ち切ることを考えます。
すなわち，粒子iとjの粒子間距離$r_{ij}$がある特定の距離$r_{cut}$より大きくなったときは力が働かないことにするのです。
$$f(r_{ij} = 0 \ for \ r_{ij} \ge r_{cut}$$
このようにすれば，無限個の粒子との相互作用を考える必要がなくなり，$r_{cut}$内の範囲にある粒子とだけ相互作用を考えれば良いのです。このようにするとエネルギー保存則が満たされなくなりますが，ある程度の大きさの$r_{cut}$ (だいたい2-3セル分)をとればその影響は小さいことが多いでです。
図. 着目している粒子(ピンク色)は$r_{cut}$内にある他の粒子からのみ力を受けることにします。
展望
このコードを起点として，以下のように様々な発展が考えられます。
- 3次元系への拡張
- 異核原子がある系
- 平衡化の作業 (NVTアンサンブル)
- 拡散係数，粘性係数，熱伝導率等の輸送係数の計算
- 自由エネルギー計算と状態方程式
など。
参考文献
[1]神山新一・佐藤明著，『分子動力学シミュレーション』，朝倉書店，1997.
[2]圧力の計算方法について:[Pythonによる科学・技術計算]二体ポテンシャル下で運動する系の圧力・ビリアルの計算，ビリアル定理，状態方程式，統計力学
[3]マクスウェル速度分布の生成について: [Pythonによる科学・技術計算] 与えられた確率密度関数を与える非一様乱数の生成，モンテカルロシミュレーション