圧力の計算

ビリアル定理

シミュレションするときにはビリアル使しますビリアル
                  ⟨ N       ⟩
P =  N-kBT--+ -1-  ∑  ri ⋅ F i
       V      3V   i=1
(1)

となりますただしいるF i  イメからのってい るのでこの使することはできませんって以下のようにしま

F
  i  2 ポテンシャルからられる
⟨∑N        ⟩     ⟨ ∑N     ∂Φ ⟩      ⟨ N∑  N∑      ∂φ(rij)⟩
     ri ⋅ F i = -     ri ⋅---- =  -         ri ⋅-------
  i                 i     ∂ri          i j⁄=i     ∂ri
(2)

となりますこのi  j  ∂∕∂rj =  - ∂ ∕∂ri  すると(2)
     ⟨                       ⟩        ⟨                            ⟩
   1- ∑N ∑N           ∂-φ(rij)       1- ∑N ∑N            ∂rij-∂φ(rij)
-  2        (ri - rj) ⋅ ∂r      = - 2        (ri - rj) ⋅ ∂r  ∂r
       i j⁄=i               i            i j⁄=i             i     ij
(3)

となりますさらに
           ∂rij
(ri - rj) ⋅----=  rij
           ∂ri
(4)

いると
⟨ N        ⟩       ⟨  N  N           ⟩
  ∑  r ⋅ F   =  - 1- ∑  ∑  r  dφ(rij)
   i  i   i       2   i j⁄=i ij ∂rij
(5)

となりますって
                  ⟨                 ⟩
     N--kBT-   -1-  ∑N ∑N    dφ(rij)
P  =    V   -  3V         rij ∂r
                     i j>i       ij
(6)

よりすることがますいるときの(1) ではなく(6) 使するがあります

圧力を求めるプログラム

式(6)の< >内の値(ビリアル)は、 力を求めるループの中で計算することが出来ます。

  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;
        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));
          virial -= rd*dphi;
         中略
        }
      }
    }

上記のプログラムで求めたビリアルと式(6)を使用すると、 シミュレーション中の圧力の計算を行うことが出来るので、 実際に計算してみてください。