力の計算
LJポテンシャル
分子動力学シミュレーションでは2体の粒子間に働く相互作用を仮定します。ア ルゴンの原子間にはファンデルワールス相互作用が働くので、 Lennard-Jonse(LJ)ポテンシャルを使用して相互作用を計算します。LJポテンシャ ルの式は
で表され、アルゴンの場合のパラメータは(A)、 (K)を使用します。 実際のシミュレーションではこの式を無次元化した形で使用します。
無次元化したLJポテンシャルの図示すると以下のようになります。
LJポテンシャルを計算するC言語のプログラムは以下のようになります。
C言語#include <stdio.h> #include <math.h> int main() { double r; double phi; for(r = 0.9; r < 4; r += 0.01){ phi = 4.0*(pow(r,-12)-pow(r,-6)); printf("%f %f\n",r,phi); } return 0; }
上記と同じ出力を行う Python のプログラムは以下のようになります。
Pythonimport numpy as np for r in np.arange(0.9,4.0,.01): phi = 4.0*(r**(-12)-r**(-6)) print(r,phi)
演習
上記のプログラムで得られた数値をgnuplotを使用してグラフにせよ。 また、プログラムを変更しLJポテンシャルの微分を求め、そのグラフを作成せよ。
Python では以下のプログラムでグラフを直接表示できます。
Pythonimport numpy as np import matplotlib.pyplot as plt r = np.arange(0.9,4.0,.01) phi = 4.0*(r**(-12)-r**(-6)) fig, ax = plt.subplots() ax.plot(r,phi) ax.set(xlabel='r',ylabel='phi') plt.show()
力の計算
粒子の位置ベクトルの座標成分をとして、 粒子が粒子に及ぼす力を、その座標成分を と置きます。この直交座標系の軸方向の単 位ベクトルをとすれば、
であり、力は
となります。上式のはナブラと呼ばれる記号で
と定義されます。粒子と粒子の距離は
と書けるので、
となります。この式を用いるとの成分は
となります。
力の計算のプログラム
上記の式を使った力計算のC言語のプログラムは以下のようになります。
C言語#include <stdio.h> #include <math.h> int main() { int i,j; double cd[6], df[3], dd[3]; double rd,dphi; for(i = 0; i < 6; i++){ cd[i] = 0; } i = 0; j = 3; for(cd[j] = 1.0; cd[j] < 4.0; cd[j] += 0.01){ dd[0] = cd[i] - cd[j]; dd[1] = cd[i+1] - cd[j+1]; dd[2] = cd[i+2] - cd[j+2]; rd = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2]; dphi = 4.0*(12.0*pow(rd,-7) - 6.0*pow(rd,-4)); df[0] = dphi*dd[0]; df[1] = dphi*dd[1]; df[2] = dphi*dd[2]; printf("%f % f\n",cd[j],df[0]); } return 0; }
上記と同じ出力を行う Python のプログラムは以下のようになります。
Pythonimport numpy as np cd = np.zeros(6) df = np.zeros(3) dd = np.zeros(3) i = 0 j = 3 for v in np.arange(1.0,4.0,0.01): cd[j] = v dd[0] = cd[i] - cd[j] dd[1] = cd[i+1] - cd[j+1] dd[2] = cd[i+2] - cd[j+2] rd = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2] dphi = 4.0*(12.0*rd**(-7) - 6.0*rd**(-4)) df[0] = dphi*dd[0] df[1] = dphi*dd[1] df[2] = dphi*dd[2] print(cd[j],df[0])
このプログラムを実行して図を作ると以下のようになります。 LJポテンシャルの図と比較すると、 LJポテンシャルの微分になっていることも確認できます。
Python では以下のプログラムでグラフを直接表示できます。
Pythonimport numpy as np import matplotlib.pyplot as plt r = np.arange(1.0,4.0,0.01) rd = r*r dphi = 4.0*(12.0*rd**(-7) - 6.0*rd**(-4))*(-r) fig, ax = plt.subplots() ax.plot(r,dphi) plt.show()
2つの粒子が(1.0,1.0,1.0)と(2.0,2.0,2.0)にある場合、 それぞれの粒子に働く力は、以下のプログラムで計算できます。
C言語#include <stdio.h> #include <math.h> enum { N = 2 }; main() { int i,j; int n = N; double cd[N*3], fc[N*3]; double df[3], dd[3]; double rd,dphi; cd[0] = 1.0; cd[1] = 1.0; cd[2] = 1.0; cd[3] = 2.0; cd[4] = 2.0; cd[5] = 2.0; for(i = 0; i < n*3; i++){ fc[i] = 0; } i = 0; j = 3; dd[0] = cd[i] - cd[j]; dd[1] = cd[i+1] - cd[j+1]; dd[2] = cd[i+2] - cd[j+2]; rd = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2]; dphi = 4.0*(12.0*pow(rd,-7) - 6.0*pow(rd,-4)); df[0] = dphi*dd[0]; df[1] = dphi*dd[1]; df[2] = dphi*dd[2]; fc[i] += df[0]; fc[i+1] += df[1]; fc[i+2] += df[2]; fc[j] -= df[0]; fc[j+1] -= df[1]; fc[j+2] -= df[2]; for(i = 0; i < n; i++){ printf("%d % f % f % f\n",i,fc[i*3],fc[i*3+1],fc[i*3+2]); } }
import numpy as np n = 2 cd = np.zeros(n*3) fc = np.zeros(n*3) df = np.zeros(3) dd = np.zeros(3) cd[0] = 1.0; cd[1] = 1.0; cd[2] = 1.0 cd[3] = 2.0; cd[4] = 2.0; cd[5] = 2.0 i = 0 j = 3 dd[0] = cd[i] - cd[j] dd[1] = cd[i+1] - cd[j+1] dd[2] = cd[i+2] - cd[j+2] rd = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2] dphi = 4.0*(12.0*rd**(-7) - 6.0*rd**(-4)) df[0] = dphi*dd[0] df[1] = dphi*dd[1] df[2] = dphi*dd[2] fc[i] += df[0] fc[i+1] += df[1] fc[i+2] += df[2] fc[j] -= df[0] fc[j+1] -= df[1] fc[j+2] -= df[2] for i in range(n): print(f'{i} {fc[i*3]: .6f} {fc[i*3+1]: .6f} {fc[i*3+2]: .6f}')
プログラムの実行結果は以下のようになります。
0 0.274348 0.274348 0.274348 1 -0.274348 -0.274348 -0.274348
演習
上記のプログラムを参考に、 2つの粒子が(1.0,1.0,1.0)と(2.0,3.0,1.0)にある場合、 それぞれの粒子に働く力が以下のようになることを確認せよ。
0 0.037786 0.075571 0.000000 1 -0.037786 -0.075571 0.000000
相互作用の数
実際のシミュレーションでは個の粒子の相互作用を計算します。相互作用 の数は、
となります。にならないのは、同じ粒子同士の力の計 算を除くためです。また、
なので、実際の計算する相互作用の数は半分となり、は
となり、下図の水色の部分のみ計算すればよいことになります。
演習
4つの粒子が(1.0,1.0,1.0), (2.0,3.0,1.0), (4.0,1.0,1.0), (5.0,5.0,1.0) にある場合、それぞれの粒子に働く力が以下のようになることを確認せよ。
0 0.048821 0.075663 0.000000 1 -0.023594 -0.085565 0.000000 2 -0.022330 0.012822 0.000000 3 -0.002897 -0.002920 0.000000
周期境界条件を考慮した力の計算
また、力の計算を行う場合は周期境界条件の処理も行う必要があります。
上図のように粒子と粒子の距離が、シミュレーションセルの長 さ の半分よりも大きい場合は、隣のイメージセルにある粒子のイメー ジ粒子との距離が実際の距離よりも小さくなるのでイメージ粒子と の相互作用を計算します。
以下では簡単のため軸方向のみで、 粒子、粒子の位置を、とし、 粒子のイメージ粒子の位置を として、
(a)
(b)
の場合を考えます。 ただし、粒子と粒子の距離をとし、 粒子と粒子の距離を とします。このとき が成り立ちます。
周期境界処理では のとき原子の力は、 原子とではなく、イメージセル内の原子と計算します。
の場合、 、 なので はとなります。 このとき、粒子と粒子の距離はイメージ粒子との距離の方が近く、 距離計算で粒子の位置はでなく、イメージ粒子の位置を使用します。 ここで、であることを用いると、 は以下のように書けます。
つまり、この条件が成り立つときは とする必要があり、このとき となるため の置き換えを行うことになります。
一方、の場合、 、 なので はとなります。 このとき、粒子と粒子の距離はイメージ粒子との距離の方が近く、 距離計算で粒子の位置はでなく、 イメージ粒子の位置を使用します。 であることを用いると、 は以下のように書けます。
つまり、この条件が成り立つときは とする必要があり、このとき となるため の置き換えを行うことになります。
軸方向と軸方向でも同様の処理を行うことで、 周期境界条件下で粒子と粒子の間に働く力を求めることができます。
演習
4つの粒子が、一辺の長さが6.0のセル内の(1.0,1.0,1.0), (2.0,3.0,1.0), (4.0,1.0,1.0), (5.0,5.0,1.0)にある場合、それぞれの粒子に働く力を周期境界 条件下で計算した場合、以下のようなることを確認せよ。
0 0.037056 0.063898 0.000000 1 -0.023594 -0.085565 0.000000 2 0.015169 -0.063898 0.000000 3 -0.028631 0.085565 0.000000
カットオフ
LJポテンシャルはが大きくなると値が急速に小さくなるので、 の大きいところでの相互作用の計算を省略することができます。 例えばLJポテンシャルでは、 以下のようにが大きくなるとLJの微分値(=力)が急激に小さくなっていきます。
通常LJポテンシャルを使用する場合は3.0(無次元量)程度でカットオフします。
演習
4つの粒子が、一辺の長さが6.0のセル内の(1.0,1.0,1.0), (2.0,3.0,1.0), (4.0,1.0,1.0), (5.0,5.0,1.0)にある場合、 それぞれの粒子に働く力を周期境界条件下で計算した場合、 以下のようなることを確認せよ。 ただし、粒子間の距離3.0で力をカットオフすること。
0 0.026113 0.063898 0.000000 1 -0.026113 -0.087244 0.000000 2 0.026113 -0.063898 0.000000 3 -0.026113 0.087244 0.000000
力計算のプログラム
以上のことを考慮に入れると力計算のプログラムは以下のようになります。
C言語pot_erg = 0; for(i = 0; i < npa*3; i++){ fc[i] = 0; } for(i = 0; i < npa*3; i += 3){ for(j = i + 3; j < npa*3; j += 3){ dd[0] = cd[i] - cd[j]; dd[1] = cd[i+1] - cd[j+1]; dd[2] = cd[i+2] - cd[j+2]; if(dd[0] < -sideh) dd[0] += side; if(dd[0] > sideh) dd[0] -= side; if(dd[1] < -sideh) dd[1] += side; if(dd[1] > sideh) dd[1] -= side; if(dd[2] < -sideh) dd[2] += side; if(dd[2] > sideh) dd[2] -= side; r2 = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2]; if(r2 < cut_off2){ r6 = r2*r2*r2; r12 = r6*r6; phi = 4.0*(1.0/(r12) - 1.0/(r6)); dphi = 4.0*(12.0/(r12*r2) - 6.0/(r6*r2)); pot_erg += phi; df[0] = dphi*dd[0]; df[1] = dphi*dd[1]; df[2] = dphi*dd[2]; fc[i] += df[0]; fc[i+1] += df[1]; fc[i+2] += df[2]; fc[j] -= df[0]; fc[j+1] -= df[1]; fc[j+2] -= df[2]; } } }
pot_erg = 0 for i in range(npa*3): fc[i] = 0 for i in range(0, npa*3, 3): for j in range(i+3, npa*3, 3): dd[0] = cd[i] - cd[j] dd[1] = cd[i+1] - cd[j+1] dd[2] = cd[i+2] - cd[j+2] if dd[0] < -sideh: dd[0] += side if dd[0] > sideh: dd[0] -= side if dd[1] < -sideh: dd[1] += side if dd[1] > sideh: dd[1] -= side if dd[2] < -sideh: dd[2] += side if dd[2] > sideh: dd[2] -= side r2 = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2] if r2 < cut_off2: r6 = r2*r2*r2 r12 = r6*r6 phi = 4.0*(1.0/(r12) - 1.0/(r6)) dphi = 4.0*(12.0/(r12*r2) - 6.0/(r6*r2)) pot_erg += phi df[0] = dphi*dd[0] df[1] = dphi*dd[1] df[2] = dphi*dd[2] fc[i] += df[0] fc[i+1] += df[1] fc[i+2] += df[2] fc[j] -= df[0] fc[j+1] -= df[1] fc[j+2] -= df[2]
演習
上記のプログラムを参考に、実際のアルゴンのプログラムで、 初期位置での原子間相互作用を計算し、 各原子に働く力を求めよ。また、力の和がゼロになることを確認せよ。
補足
sidehはsideの半分、cut_off2はsidehの2乗とする。
sideh = side*0.5;
cut_off2 = sideh*sideh;
このページではプログラミング言語として C 言語と Python を使用していますが、 Python で作ったプログラムの動作速度は C 言語に比べて遅いので、 MD シミュレーションのプログラムは C 言語のものを使うことを想定しています。
Python のプログラムは C 言語のプログラムを理解することを目的に作られており、 なるべく C 言語に近い書き方になっています。 そのため、 Python のプログラムでは、効率の悪い書き方になっている箇所もありますので注意してください。