圧力の計算
ビリアル定理
分子動力学シミュレーションで圧力を計算するときには、ビリアル定理を使用し ます。ビリアル定理の式は
(1)
となります。ここで、は粒子 に働く力であり、 これは他の粒子 との相互作用による力の和であるので
(2)
です。周期境界条件を用いる場合ははイメージ粒子との 力となることがあり、位置が実際ものとは違う値となるため、 式(1)をそのまま使って圧力を計算することはできません。 従って以下のように式を変形します。
(3)
ここでを用いると
(4)
となります。従って、式(1)は
(5)
とするができ、とは 周期境界条件下でも常に対応することになります。 そのため、周期境界条件を用いるときの圧力計算は式(1)ではなく、 式(5)を使用します。
圧力を求めるプログラム
式(5)の内の値(ビリアル)は、 力を求めるループの中で計算することができます。
C言語double virial; (略) virial = 0; /* 力の計算 */ for(i = 0; i < pn*3; i += 3){ for(j = i + 3; j < pn*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)); virial += r2*dphi; (略) } } }
virial = 0.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 virial += r2*dphi (略)
上記のプログラムで求めたビリアルと式(5)を使用すると、 シミュレーション中の圧力の計算を行うことができるので、 実際に計算してみてください。
このページではプログラミング言語として C 言語と Python を使用していますが、 Python で作ったプログラムの動作速度は C 言語に比べて遅いので、 MD シミュレーションのプログラムは C 言語のものを使うことを想定しています。
Python のプログラムは C 言語のプログラムを理解することを目的に作られており、 なるべく C 言語に近い書き方になっています。 そのため、 Python のプログラムでは、効率の悪い書き方になっている箇所もありますので注意してください。