processingの備忘録 -LJポテンシャル-

こんにちは。「processingの備忘録」シリーズ、第3回目です。
これまでに、「クラス」、「衝突」の話をしましたが、今回は「相互作用」を入れてみます。

ここでの「相互作用」は、Lennard-Jonesポテンシャル（LJポテンシャル）から説明できる引力と斥力です。
Lennard-Jonesポテンシャルは分子間の引力・斥力を表したもので、 $$ u = 4 \epsilon \left( \left( \frac{σ}{r} \right)^{12} - \left( \frac{σ}{r} \right)^{6} \right) $$ という形で表されるものが一般的で、$\epsilon$は最も安定な状態でのポテンシャル、σは分子の直径に対応する値をとります。

さて、エネルギーを微分して-1倍すると力となり、運動方程式に当てはめると、加速度として表せます。
加速度は速度の変化量なので、プログラム上では速度の値に加速度を足していけばよいことになります。
これより加速度は、 $$ a = \frac{4 \epsilon}{m} \left( 12 \frac{σ^{12}}{r^{13}} - 6 \frac{σ^{6}}{r^{7}} \right) $$ と表せるので、プログラムに組み込んでみます。
（$\epsilon$、σの部分はそれっぽい値に調節しました。）

int  n = 50;
Circle[] circles = new Circle[n];
 
void setup(){
    size(500, 500);
 
    for(int i = 0; i < n; i ++){
        circles[i] = new Circle(i);
    }
}
 
void draw(){
    background(0);
 
    for(int i = 0; i < n; i ++){
        circles[i].display();
        circles[i].move();
        circles[i].wall();
        circles[i].collide();
    }
}
 
class Circle{
    int id;
    float x, y, m, d, v, vx, vy, epsilon;
    float[] odistance = new float[n];
 
    Circle(int _id){
        id = _id;
        d = 50;
        x = id%10*(width/10)+d/2;
        y = id/10*(height/10)+d/2;
        m = 1;
        v = 3;
        vx = random(-v, v);
        vy = random(-v, v);
        epsilon = 1;
        for(int i = 0; i < n; i ++){
            odistance[i] = width + height;
        }
    }
 
    void display(){
        ellipse(x, y, d, d);
    }
 
    void move(){
        x += vx;
        y += vy;
    }
 
    void collide() {
        for (int i = id+1; i < n; i ++){
            float dx = circles[i].x - x;
            float dy = circles[i].y - y;
            float distance = sqrt(dx*dx + dy*dy);
            float angle = atan2(dy, dx);
 
            if (distance < d && distance < odistance[i]){
                float v1h = vx*cos(angle)+vy*sin(angle);
                float v1v = vx*sin(angle)-vy*cos(angle);
                float v2h = circles[i].vx*cos(angle)+circles[i].vy*sin(angle);
                float v2v = circles[i].vx*sin(angle)-circles[i].vy*cos(angle);
 
                float v1ha = (m*v1h+2*circles[i].m*v2h-circles[i].m*v1h)/(m+circles[i].m);
                float v2ha = (2*m*v1h+circles[i].m*v2h-m*v2h)/(m+circles[i].m);
 
                vx = v1ha*cos(angle)+v1v*sin(angle);
                vy = v1ha*sin(angle)-v1v*cos(angle);
                circles[i].vx = v2ha*cos(angle)+v2v*sin(angle);
                circles[i].vy = v2ha*sin(angle)-v2v*cos(angle);
            }
            if(distance > d){ //丸同士が接していないとき
                float a = 4*epsilon/m*(12*pow(d, 12)/pow(distance, 13)-6*pow(d, 6)/pow(distance, 7));
 
                vx -= a*cos(angle);
                vy -= a*sin(angle);
                circles[i].vx += a*cos(angle);
                circles[i].vy += a*sin(angle);
            }
            odistance[i] = distance;
        }
    }
 
    void wall(){
        if(x < d/2){
            x = d/2;
            vx *= -1;
        }
        if(x > width-d/2){
            x = width-d/2;
            vx *= -1;
        }
        if(y < d/2){
            y = d/2;
            vy *= -1;
        }
        if(y > height-d/2){
            y = height-d/2;
            vy *= -1;
        }
    }
}

これを実行してみると、50個の丸が互いに衝突し合うような動きが見えると思います。
速度を遅くしてみる(v = 0.5など)と、丸同士は引き合い、集まっていきます。
これを実際の分子と考えると、これは相転移にあたると言えます。
丸の速度は温度に対応していると考えれば、v = 5のときは気体として存在しますが、温度を低くすることで丸（分子）同士の引き合う力の方が強くなり、集まります（液体、固体）。
液体と固体は集まりの度合いによって区別できます。

このままでは、速度を変えるにはいちいちプログラムを書き換えなければなりませんが、実行した状態で速度を変えられるようにしてみます。
ここは工夫の余地がありますが、とりあえず以下のようなコードを書き足してみます。

void mousePressed(){
    for(int i = 0; i < n; i ++){
        if(mouseX < width/2){
            circles[i].vx *= 0.9;
            circles[i].vy *= 0.9;
        }
        if(mouseX > width/2){
            circles[i].vx *= 1.1;
            circles[i].vy *= 1.1;
        }
    }
}

これを追加することで、画面の左半分をクリックするたびに速度が0.9倍になり、右半分をクリックするたびに速度が1.1倍になるようになります。
実行してすぐは気体状態ですが、左半分を15～20回ほどクリックするとなんとなく集まってきて液体状態となり、さらに10回ほどクリックするとその集まりはある程度並び方を保って固体状態となります。

ここまでで分子の衝突、相互作用を追加でき、分子の速度を変えることで相転移する様子が見えるようになりました。
目標はひとまず達成したので、以上で「processingの備忘録」シリーズは一区切りとしたいと思います。
この後、コンピュータ世界のスケールと現実世界のスケールを対応させたり、重力を追加したりといろいろと加えると面白いこともあるのですが、とりあえず基本の骨格のようなものはつくることができたので、このあたりにしたいと思います。
最後まで読んでいただき、ありがとうございました。

おまけ
さらに、大きな丸を一つ描くようなプログラムを加え、その他諸々の機能を付けたり消したりした結果、ブラウン運動っぽい感じになりました。
クリックすれば丸を増やせます。青色の線がブラウン運動の軌跡で、「c」キーを押せば軌跡をリセットできます。