力の計算

LJポテンシャル

シミュレションでは 2 互作しますアルゴン ファンデルワルス互作くのでLennard-Jonse(LJ) ポテンシャル使して互作 しますLJ ポテンシャル
          {( σ )12  ( σ )6}
φ(r) = 4ɛ    --   -   --
             r        r
(1)

されアルゴンパラメσ =  3.40  (Å)ɛ∕kB = 120  (K) 使します シミュレションではこのした使します
         {              }
          ( 1)12   ( 1)6
φ(r) = 4    r    -   r
(2)

無次元化したLJポテンシャルの図示すると以下のようになります。

LJポテンシャル

LJポテンシャルを計算するプログラムは

#include <stdio.h>
#include <math.h>
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);
  }
}

となります。

演習

上記のプログラムで得られた数値をgnuplotを使用してグラフにせよ。 また、プログラムを変更しLJポテンシャルの微分を求め、そのグラフを作成せよ。

力の計算

i  ベクトルr
  i  x ,y ,z
 i  i  i  としてj  i  ぼす F ij  そのXij,Yij,Zij  きますこの交座x, y,z  ベク トルex,ey,ez  とすれば
ri = xiex + yiey + ziez
(3)

であり
F ij = - ∇i φ(rij)
(4)

となります∇ ナブラばれる
∇  = e  -∂--+ e  ∂--+  e -∂-
  i    x∂xi    y ∂yi    z∂zi
(5)

されますi  j  r
 ij
                 ∘ ----------------------------------
rij = |ri - rj| =   (xi - xj)2 + (yi - yj)2 + (zi - zj)2
(6)

けるので
∂rij = xi---xj
∂xi      rij
(7)

となりますこのいるとF ij  x, y,z
        dφ-(rij)xi---xj
Xij = -   drij    rij
(8)

        dφ(r  )y -  y
Yij = - ----ij--i----j
          drij    rij
(9)

Z   = - dφ-(rij)zi --zj
  ij       drij    rij
(10)

となります

力の計算のプログラム

上記の式を使った力の計算プログラムは以下のようになります。

#include <stdio.h>
#include <math.h>
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]);
  }
}

このプログラムを実行して図を作ると以下のようになります。 LJポテンシャルの図と比較すると、 LJポテンシャルの微分になっていることも確認できます。

LJポテンシャルの力

2つの粒子が(1.0,1.0,1.0)と(2.0,2.0,2.0)にある場合、 それぞれの粒子か働く力は、以下のプログラムで計算できます。

#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]);
  }
}

このプログラムをコンパイルし、実行すると

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

相互作用の数

シミュレションではN  互作します互作I
I = N (N  - 1)
(11)

となりますN 2   にならないのはF
  ii  くためです
F ij = - F ji
(12)

なのでする互作となりI
I =  1N (N -  1)
     2
(13)

となり水色部分のみすればよいことになります

計算する相互作用の数

演習

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

周期境界条件を考慮した力の計算

またがあります

力の計算の周期境界条件の処理

上図のように j  r  シミュレションセルL  よりもきい イメジセルにあるj  イメとのr′ r  よりも くなるのでイメとの互作します

以下では簡単のため 軸方向のみで、 粒子、粒子の位置をとし、 粒子 のイメージ粒子 の位置を として、

(a)
(b)

の場合を考えます。 ただし、粒子と粒子の距離をとし、 粒子と粒子の距離を とします。

(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ポテンシャルはrが大きくなると値が急速に小さくなるので、 rの大きいところでの相互作用の計算を省略することができます。 例えばLJポテンシャルでは、 以下のようにrが大きくなるとLJの微分値(=力)が急激に小さくなっていきます。

 r  LJの微分値
1.5 1.158029
2.0 0.181641
2.5 0.038999
3.0 0.010944
3.5 0.003726
4.0 0.001464

通常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
force0.c force1.c force2.c force3.c

力の計算のプログラム

以上のことを考慮に入れると以下のようなプログラムになります。

pot_erg = 0;
for(i = 0; i < pn*3; i++){
  fc[i] = 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;
    rd  = dd[0]*dd[0] + dd[1]*dd[1] + dd[2]*dd[2];
    if(rd < cut_off2){
	r6 = rd*rd*rd;
	r12 = r6*r6;
	phi = 4.0*(1.0/(r12) - 1.0/(r6));
	dphi = 4.0*(12.0/(r12*rd) - 6.0/(r6*rd));
	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];
    }
  }
}

argon_fc.c

演習

上記のプログラムを参考に、実際のアルゴンのプログラムで、 初期位置での原子間相互作用を計算し、 各原子に働く力を求めよ。また、力の和がゼロになることを確認せよ。

補足
sidehはsideの半分、cut_off2はsidehの2乗とする。
sideh = side*0.5;
cut_off2 = sideh*sideh;