粒子同士の相互作用(2)

分子間の引力、斥力

分子間の引力、斥力を計算に加えましょう。

ここからの数式展開は長いですが、最終目標はLJポテンシャルによる分子間の引力、斥力を計算し、それによって生じる加速度をpixel(長さ), frame(時間)単位で求めることです。 (少し下の方に青色で示されています。)

分子間の相互作用にはLennard-Jonesポテンシャルを用いることにします。

 \displaystyle u(r) = 4 \varepsilon \Bigl[ \Bigl( \frac{\sigma}{r} \Bigr)^{12} - \Bigl( \frac{\sigma}{r} \Bigr)^{6} \Bigr] 

粒子間に働く力は、ポテンシャルをr で微分して – を付ければ得られます。

 \displaystyle f(r) = -\frac{{\rm d}u(r)}{{\rm d}r} 

計算すると

 \displaystyle f(r) = 4 \varepsilon \Bigl[ 12 \Bigl( \frac{\sigma^{12}}{r^{13}} \Bigr) - 6 \Bigl( \frac{\sigma^6}{r^7} \Bigr) \Bigr] 

f がわかれば、運動方程式から物体の運動を記述できます。

 \displaystyle f = ma 

a は加速度なので

 \displaystyle a = \frac{{\rm d}v}{{\rm d}t} = \frac{f}{m} 

v の式に順次代入していきます。

 \displaystyle v = v_0 + \Bigl( \frac{{\rm d}v}{{\rm d}t} \Bigr) \cdot \Delta t \\ = v_0 + \Bigl( \frac{f}{m} \Bigr) \cdot \Delta t \\ = v_0 +  \frac{1}{m}\cdot 4 \varepsilon \Bigl[ 12 \Bigl( \frac{\sigma^{12}}{r^{13}} \Bigr) - 6 \Bigl( \frac{\sigma^6}{r^7} \Bigr) \Bigr]  \cdot \Delta t 

さらに質量は分子量M で、

 \displaystyle (m / {\rm kg})  = \frac{(M / {\rm g \cdot mol^{-1}})}{(N_A / {\rm mol^{-1}})} \cdot 10^{-3} 

ε をK単位で

 \displaystyle (\varepsilon / {\rm J}) = (\varepsilon_{\rm K} / {\rm K}) \cdot (k_{\rm B} / {\rm J \cdot K^{-1}}) 

表すことにすると

 \displaystyle v = v_0 +  \frac{N_A}{M \cdot 10^{-3}} \cdot 4 \cdot \varepsilon_{\rm K} \cdot k_{\rm B} \Bigl[ 12 \Bigl( \frac{\sigma^{12}}{r^{13}} \Bigr) - 6 \Bigl( \frac{\sigma^6}{r^7} \Bigr) \Bigr]  \cdot \Delta t 

となります。

単位合わせ

さて、上式を計算すれば分子間に働く力によって分子の運動(速度)がいかに変化するかがわかります。

具体的には、以前重力の効果を生じさせるときに

vy = vy + 0.04;

というように、ステップごとにv を変化させたのと同じように分子間の力によってv を変えていけばよい(加速度を生じさせてやればよい)わけです。

ここで単位が問題となります。通常、上式にはSI単位系に合わせた、時間はs (秒)、距離はm (メートル)の単位の数値を入れていきます。しかし、コンピュータ上で分子を動かすためには、時間はframe(1ステップあたりの時間)、距離はpixel(画面上のピクセル)単位での値が必要です。(分子の速度の単位は m s-1の代わりにpixel frame-1 となる)

まず、前ページと同じように、両者の対応を定義してやります。

 1 {\rm\ frame} = t_{\rm scale} {\rm\ ps} \\ 1 {\rm\ pixel} = s_{\rm scale} {\rm\ pm} \\ 

tscalesscaleは、それぞれ1 frame, 1 pixel が(今回のシミュレーションの世界の中で)何ps、何pmに相当するかを指定するパラメータで、プログラム冒頭で変数として宣言し、変更できるようにします。
frameとpixelはそれぞれ時間と距離の単位になります。上式に単位をつけた数値を代入していき、現れたs と m をframe とpixelに置き換えていってやります。

Lennard-Jonesのパラメータは窒素のものを入れてみます。

 \sigma = 370 {\rm\ pm} \\ \\ \varepsilon_{\rm K} = 95.1 {\rm\ K} \\ \\ M = 28.0 {\rm\ g \ mol^{-1}} 
 \displaystyle v = v_0 +  \frac{6.02 \times 10^{23} {\rm\ mol^{-1}}}{28.0 \times 10^{-3} {\rm\ kg \ mol^{-1}}} \cdot 4 \cdot (95.1 {\rm\ K}) \cdot (1.381 \times 10^{-23} {\rm\ kg\ m^2\ s^{-2}\ K^{-1}}) \Bigl[ 12 \Bigl\{ \frac{(370 \times 10^{-12} {\rm\ m})^{12}}{(d {\rm\ pixel})^{13}} \Bigr\} - 6 \Bigl\{ \frac{(370 \times 10^{-12} {\rm\ m})^6}{(d {\rm\ pixel})^7} \Bigr\} \Bigr]  \cdot (1 {\rm\ frame}) 

d はpixel 単位で表した粒子間の距離(の数値部分)です。
さて、現れたs と m を

 \displaystyle 1 {\rm\ s} = \frac{10^{12}}{t_{\rm scale}} {\rm\ frame}  \\ \\ \displaystyle 1 {\rm\ m} = \frac{10^{12}}{s_{\rm scale}} {\rm\ pixel}  \\ 

を使って、(あたかも代数記号のように)置き換えていきます。

 \displaystyle v = v_0 +  \frac{6.02 \times 10^{23} {\rm\ mol^{-1}}}{28.0 \times 10^{-3} {\rm\ kg \ mol^{-1}}} \cdot 4 \cdot (95.1 {\rm\ K}) \cdot \frac{t_{\rm scale}^2}{s_{\rm scale}^2}(1.381 \times 10^{-23} {\rm\ kg\ pixel^2\ frame^{-2}\ K^{-1}}) \Bigl[ 12 \Bigl\{ \frac{370^{12}}{d^{13} s_{\rm scale}^{12} {\rm\ pixel}} \Bigr\} - 6 \Bigl\{ \frac{370^6}{d^{7} s_{\rm scale}^6 {\rm\ pixel}} \Bigr\} \Bigr]  \cdot (1 {\rm\ frame}) 

単位を整理すると、下記のようにpixel / frame単位の式が得られます。

 \displaystyle v = v_0 \  + \  \frac{6.02 \times 10^{23} }{28.0 \times 10^{-3} } \cdot 4 \cdot (95.1) \cdot \frac{t_{\rm scale}^2}{s_{\rm scale}^2} \cdot (1.381 \times 10^{-23} ) \Bigl\{ \Bigl( \frac{12}{d^{13}} \cdot \frac{370^{12}}{s_{\rm scale}^{12}} \Bigr) - \Bigl( \frac{6}{d^{7}} \cdot \frac{370^6}{s_{\rm scale}^6 } \Bigr) \Bigr\}  \cdot {\rm\ pixel\ frame^{-1}} 

(d はpixel 単位で表した粒子中心間の距離)

Processingでは乗数はpow(a, b)を使って計算できます。
例えば1023はpow(10, 23)、37012はpow(370,12)、d7はpow(d,7)です。

プログラムへの組み込み

さて、いよいよ分子間の相互作用をプログラム中に組み込みましょう。
処理はParticleクラス内のcollide()で行うこととします。

  void collide() {
    float elast = 1;
    //他の粒子との衝突
    if(visible>0){ //自分自身が存在しているとき
      for (int i = id + 1; i < particle_number; i++) { //他の粒子との相互作用を計算
        //粒子間距離の算出
        float dx = particle[i].x - x;
        float dy = particle[i].y - y;
        float distance = sqrt(dx*dx + dy*dy); //粒子間距離

        float minDist = particle[i].r/2 + r/2; //粒子の半径の和

        if(particle[i].visible>0){ //相手の粒子が存在しているとき
          float angle = atan2(dy, dx);
          if (distance < 0.9 * minDist) { //半径の和の0.9倍より距離が小さく
            if(distance < old_distance[i]){ //かつ距離が縮まっている場合、衝突として処理する
              //粒子間のなす角を算出
              //速度を衝突方向の水平成分と垂直成分に分ける
              float v1h = vx*cos(angle)+vy*sin(angle);
              float v1v = vx*sin(angle)-vy*cos(angle);
              float v2h = particle[i].vx*cos(angle)+particle[i].vy*sin(angle);
              float v2v = particle[i].vx*sin(angle)-particle[i].vy*cos(angle);

              //垂直成分はそのまま
              //水平成分は運動量保存則と弾性係数から新しい値を計算する
              float nv1h = (v2h-v1h)*(1+elast)/(m/particle[i].m+1)+v1h;
              float nv2h = (v1h-v2h)*(1+elast)/(particle[i].m/m+1)+v2h;

              //速度の水平、垂直成分を画面上のx方向、y方向に書き戻す
              vx = nv1h*cos(angle)+v1v*sin(angle);
              vy = nv1h*sin(angle)-v1v*cos(angle);
              particle[i].vx = nv2h*cos(angle)+v2v*sin(angle);
              particle[i].vy = nv2h*sin(angle)-v2v*cos(angle);
            }
          } //衝突処理終わり

          else if(distance < 10 * minDist){ //粒子間距離が粒子直径の10倍より近い時
          //分子間の引力、斥力を処理
            float a;

            //f に両者間に働く力によって生じる加速度を入れる。ヒント:距離はdistanceという変数に入っている。(pixel単位)
            a = 0;

            vx = vx - a * cos(angle); //x成分と
            vy = vy - a * sin(angle); //y成分に分けて速度に足す
            particle[i].vx = particle[i].vx + a * cos(angle); //作用反作用で、相手の粒子にも(反対方向に)加速度を生じさせる
            particle[i].vy = particle[i].vy + a * sin(angle);
          } //分子間の相互作用の処理終わり
        }
        old_distance[i] = distance; //計算した粒子間距離は保存しておく
      }
    }
  }
  • 15行目~35行目まで
    粒子間距離がσ の0.9倍以下の時
    前回やったように、衝突(弾性衝突)として処理します。(後述)
  • 37行目~48行目まで
    ここが粒子間の相互作用を処理する部分となります。 
    • 上で導いた式を使って、42行目でa=0となっている右辺に、速度の増分(加速度; 1 pixel frame-1単位)を記入してください。
      (v = v0 + at。1 frameあたりなので、t = 1 frameとなりt はプログラム上は消える。)
    • 44~47行目では、42行目で定義されたaをx成分、y成分に分け、それまでの速度に加えています。相互作用は2つの粒子の対の両方に働くので、相手の粒子にも(加速度の方向を反対にして)同じ処理をしています。
    • この処理は全ての粒子対について行っています。つまり、複数の分子間で生じる相互作用は、単純に「足し合わせ」で表されると仮定しているわけです。

分子間のポテンシャルについて

以下、追加の話です。
ここで用いている分子間のポテンシャルについて少々補足します。

LJポテンシャルはこのような形状です。(青)

横軸はσ、縦軸は-ε で規格化してあります。
赤は斥力項(r の-12乗の項)、緑は引力項(r の-6乗の項)です。

 \displaystyle u(r) = 4 \varepsilon \Bigl[ \Bigl( \frac{\sigma}{r} \Bigr)^{12} - \Bigl( \frac{\sigma}{r} \Bigr)^{6} \Bigr] 

さて、遠い方から見ていくと、上の青線はr = 3 σ くらいまで行くとほぼ0になっていることがわかります。
上のプログラムでは38行目でσ の10倍より近い時だけ計算して、それ以遠については実質0なので、計算をしませんでした。
このような距離を「カットオフ半径」といい、これによって計算時間を短縮することができます。

分子を近づけていくと、r = 1.1 σ のあたり(正確には2の1/6乗≈1.122 σ)で青線は最小値となります。この距離より遠いと分子間には引力が働き(緑の効果が赤より大きくなる)、この距離より近いと分子間には斥力が働きます(赤の効果が緑より大きくなる)。
従って、分子の相対速度が小さい時(= 温度が低い時)、2つの分子はこの距離を中心に振動運動をするはずです。

r = σ は分子間のポテンシャルが0となる位置であり、分子同士が速い速度で衝突すると、分子はこの距離より近づきます。画面上で分子を直径σ で描いているなら、両者は「めり込む」ことがある、ということです。
計算(シミュレーション)を行う際、ここで問題が生じます。r < σ の領域では、少し近づいただけでやたらとポテンシャルが上昇します。
このプログラムは時間を(frame単位で)とびとびで計算しているため、速い相対速度で衝突すると、分子同士がめり込んだ後、ものすごい強い力で両者が反発してしまうのです。

これを避けるため、r = 0.9 σ で弾性衝突をする、という処理(15行目)を加えています。
先週やったような、r = σ で弾性衝突をする、という計算(剛体球モデル、ですね)は、ポテンシャルとして次のような仮定をしていることに相当します。

すなわち、rσ のとき U = ∞, r > σ のとき U = 0 です。
2つの粒子は、どんな速い速度で衝突してもrσ となることはありません。

上のプログラムではr = 0.9 σ で弾性衝突ですから

こんなポテンシャルということになります。
本当は(LJポテンシャルに従うとすると)点線のようなポテンシャルですが、r = 0.9 σ より近い時はLJ相互作用を計算せず、剛体球として処理しているということですね。
時間の刻み幅が分子の運動に対して十分に小さければ、r = 0.9 σ まで近づかず、r = 1 σ ~ 0.9 σ の間で、 LJ相互作用によって跳ね返るはずです。

追加の話のさらに追加になりますが、最初から分子がr = 0.9 σ より近い場合、(分子を生成させるとき、同じ場所でクリック連打などをしたときに起こる)
相互作用は0としていますから、より正確には

というようなポテンシャルを仮定していることになりますね。
(r < 0.9 σ では、分子間距離が近づいている場合剛体球として処理、遠ざかっている場合はF = 0。1)dU / dr) = 0

void collide() {
float elast = 1;
//他の粒子との衝突
if(visible>0){ //自分自身が存在しているとき
for (int i = id + 1; i < particle_number; i++) { //他の粒子との相互作用を計算
//粒子間距離の算出
float dx = others[i].x – x;
float dy = others[i].y – y;
float distance = sqrt(dx*dx + dy*dy); //粒子間距離 

float minDist = others[i].r/2 + r/2; //粒子の半径の和

if(others[i].visible>0){ //相手の粒子が存在しているとき
float angle = atan2(dy, dx);
if (distance < 0.9 * minDist) { //半径の和の0.9倍より距離が小さく
if(distance < old_distance[i]){ //かつ距離が縮まっている場合、衝突として処理する
//粒子間のなす角を算出
//速度を衝突方向の水平成分と垂直成分に分ける
float v1h = vx*cos(angle)+vy*sin(angle);
float v1v = vx*sin(angle)-vy*cos(angle);
float v2h = others[i].vx*cos(angle)+others[i].vy*sin(angle);
float v2v = others[i].vx*sin(angle)-others[i].vy*cos(angle);

//垂直成分はそのまま
//水平成分は運動量保存則と弾性係数から新しい値を計算する
float nv1h = (v2h-v1h)*(1+elast)/(m/others[i].m+1)+v1h;
float nv2h = (v1h-v2h)*(1+elast)/(others[i].m/m+1)+v2h;

//速度の水平、垂直成分を画面上のx方向、y方向に書き戻す
vx = nv1h*cos(angle)+v1v*sin(angle);
vy = nv1h*sin(angle)-v1v*cos(angle);
others[i].vx = nv2h*cos(angle)+v2v*sin(angle);
others[i].vy = nv2h*sin(angle)-v2v*cos(angle);
}
} //衝突処理終わり

else if(distance < 10 * minDist){ //粒子間距離が粒子直径の10倍より近い時
//分子間の引力、斥力を処理
float a;

//f に両者間に働く力によって生じる加速度を入れる。ヒント:距離はdistanceという変数に入っている。(pixel単位)
a = kr/pow(distance, 13)-ki/pow(distance, 7);

vx = vx – a * cos(angle); //x成分と
vy = vy – a * sin(angle); //y成分に分けて速度に足す
others[i].vx = others[i].vx + a * cos(angle); //作用反作用で、相手の粒子にも(反対方向に)加速度を生じさせる
others[i].vy = others[i].vy + a * sin(angle);
} //分子間の相互作用の処理終わり
}
old_distance[i] = distance; //計算した粒子間距離は保存しておく
}
}
}

脚注   [ + ]

1. dU / dr) = 0