動径分布関数の計算

動径分布関数

分布g(r)  はあるからのr  するですあるからの r  r + dr  にあるn(r)  としたとき
-n(r)--
4πr2dr
(1)

となりますこのρ  るとあるi  分布g (r)
 i  となりま
gi(r) = --n(r)--
        4πr2dr ρ
(2)

g(r)  ったものが分布となります
                 ⟨n(r)⟩-
g(r) = ⟨gi(r)⟩ = 4πrdrρ
(3)

距離rとr+dr内にある粒子

動径分布関数はX線散乱実験から求めることが出来ます。 下図は温度85Kでの液体Arの動径分布関数の 実験値 (Yarnell et al, Phys. Rev. A, 7, 2130 (1973))です。

Arの動径分布関数の実験値

r の小さいところでg(r) の値はゼロになります。 これは原子同士が重なり合わないことを意味しています。 3.7Å付近のピークはある粒子に隣接する粒子が、 その位置にいる確率が高いことを表しています。 最初のピーク位置のr の2倍、3 倍、...、の位置のピークは隣の隣の粒子、 その隣の粒子が存在しやすい場所であることを意味しています。

動径分布関数を求めるプログラム

粒子の間の距離は力を求めるときに計算するので、その距離を利用してn(r) を求めるプログラムは以下の通りです。求まったn(r)  からと式(3)を使用すると、動径分布関数を求めることが出来ます。

  enum {
    GR_NUM = 100
  };
  int gr[GR_NUM];
  double r,dr = 0.1;

  中略

    /* 力の計算 */
    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];
        r = sqrt(rd);
        k = r/dr
        if(k < GR_NUM){
          gr[k]++;
        }
        if(rd < cut_off2){
        中略
        }
      }
    }

動径分布関数を求めるプログラムを作成し、 シミュレーションで得られる動径分布関数と実験値 と比較してみてください。