動径分布関数の計算
動径分布関数
動径分布関数はある粒子からの距離に他の粒子が存在する確率で す。ある粒子からの距離との球殻の間にある粒子の数を としたとき、球殻内の粒子の密度は
となります。この量を系の平均密度をで割ると、ある粒子の動径 分布関数となります。
の粒子平均、時間平均を取ったものが系の動径分布関数となります。
動径分布関数はX線散乱実験から求めることができます。 下図は温度85Kでの液体Arの動径分布関数の 実験値 (Yarnell et al, Phys. Rev. A, 7, 2130-2144 (1973))です。
の小さいところで の値はゼロになります。 これは原子同士が重なり合わないことを意味しています。 3.7Å付近のピークはある粒子に隣接する粒子が、 その位置にいる確率が高いことを表しています。 最初のピーク位置の の2倍、3 倍、…、の位置のピークは隣の隣の粒子、 その隣の粒子が存在しやすい場所であることを意味しています。
動径分布関数を求めるプログラム
粒子の間の距離は力を求めるときに計算するので、その距離を利用して を求めるプログラムは以下の通りです。求まった と式(3)を使用すると、動径分布関数を求めることができます。
C言語enum { GR_SIZE = 100 }; int k, gr[GR_SIZE]; /* 動径分布関数計算用の配列 */ double dr = 0.05; /* 動径分布関数計算用の dr */ for(i = 0; i < GR_SIZE; i++){ gr[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(step > adjust_temp_step){ /* 動径分布関数を求めるための処理 */ k = sqrt(r2)/dr; if(k < GR_SIZE){ gr[k]++; } } if(r2 < cut_off2){ (略) } } } (略)
(略) gr_size = 100 gr = np.zeros(gr_size) # 動径分布関数計算用の配列 dr = 0.05 # 動径分布関数計算用の dr (略) # 力の計算 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 step > adjust_temp_step: # 動径分布関数を求めるための処理 k = int(math.sqrt(r2)/dr) if k < gr_size: gr[k] += 1 if r2 < cut_off2: (略) (略)
動径分布関数を求めるプログラムを作成し、 シミュレーションで得られる動径分布関数と 実験値 (Yarnell et al, Phys. Rev. A, 7, 2130-2144 (1973)) と比較してみてください。
このページではプログラミング言語として C 言語と Python を使用していますが、 Python で作ったプログラムの動作速度は C 言語に比べて遅いので、 MD シミュレーションのプログラムは C 言語のものを使うことを想定しています。
Python のプログラムは C 言語のプログラムを理解することを目的に作られており、 なるべく C 言語に近い書き方になっています。 そのため、 Python のプログラムでは、効率の悪い書き方になっている箇所もありますので注意してください。