力の計算

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  よりも くなるのでイメとの互作します

演習

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;