ガラスのプログラム

理系ネタ

ガラスは本当に固体か？ −コンピュータシミュレーションと情報理論による予測−
http://www.kyoto-u.ac.jp/ja/research/research_results/2014/150122_1.html

元論文を読むと、１３００個のビリヤード球の衝突の計算で結果を出されてます。腕に憶えのあるプログラマは挑戦してみては？
その道筋を書いときます。

1:玉の衝突の計算

半径１の玉を箱の中で動かします。壁があるとしてもいいし、壁に来たら反対から出てくるドラクエ式でもいい。
まずは時間DTだけ進めて速度＊DTだけ玉を動かし、DTのうちに壁にぶつかるかどうかを判定し、ぶつかるなら、衝突する時間まで時間を進めて玉を壁にひっつけ、そこで速度を反転させます。

#define DIM 2
double pos[DIM], v[DIM];
double min[DIM],max[DIM];
double pos_new[DIM];
double dt;
...

  while(1){

  double t_hit= dt;
  int hit_flag= -1;

  for(m=0;m<DIM;m++)
  { 
     pos_new[m]=pos[m]+v[m]*dt;
     if(pos_new[m]<min[m])
     {
          double tt = -(pos[m]-min[m])/v[m];
          if(hit_flag==-1 || tt<t_hit)
          {
              t_hit = tt;
              hit_flag = m;
          }
     }
     if(pos_new[m]>max[m])
     {
          double tt = (max[m]-pos[m])/v[m];
          if(hit_flag==-1 || tt<t_hit)
          {
              t_hit = tt;
              hit_flag = m;
          }
     }
  }//m

  for(m=0;m<DIM;m++)
    pos[m] += t_hit * v[m];
  
  if (hit_flag!=-1)
     v[hit_flag] = -v[hit_flag]; 


  }//while

2：玉を増やす

二個以上に玉を増やします。玉同士の衝突も判定します。玉の中心同士の距離が２以下になるかどうかで判定。
玉同士、玉と壁などそれぞれ判定し、ぶつかるのがあるなら、一番早い時間にぶつかるイベントを見つけます。
その時間まで時間を進めて、衝突による速度の反転処理を行います。
玉iが場所Riにいて速度Viの時（R,Vはベクトル）、玉iとjが衝突する時間tは(Ri-Rj +t(Vi-Vj))^2=2^2 という二次方程式を解いてtの小さいほうの解からわかります。解がないときは衝突なし。
これで完成

#define N 100
#define Sq(x) ((x)*(x))
  // position, velocity, radius
  double pos[N][DIM], v[N][DIM], r[N];
  ...
  double tmin = dt;
  int imin=-1, jmin=-1;
  for(i=0;i<N;i++)
  {
    for(j=0;j<i;j++)
    {
        double a = 0.0;
        double b = 0.0;
        double c = -Sq(r[i]+r[j]);
        for(m=0;m<DIM;m++)
        {
           double dp = pos[i][m] - pos[j][m];
           double dv = v[i][m] - v[j][m];
           a += dv*dv;
           b += 2*dv*dp;
           c += dp*dp;
        }
        double det = b*b - 4*a*c;
        if(det<0.0) continue;
        if(a < 1.0e-10)continue;
        //課題:sqrtの計算回数を極力減らす工夫をせよ
        double t_c = (-b -sqrt(det))/(2*a);
        if(t_c > 0.0 && t_c < tmin)
        {
          tmin = t_c;
          imin=i;
          jmin=j;
        }
    }//j
  }//i

3実験

二次元でも三次元でもいいけど、箱に玉を隙間無くぎっしり詰めます。そして玉の半径を0.99倍にして、各玉にランダムな速度を与えます。
ガチガチガチガチとぶつかりまくります。しかし回りを玉に囲まれてるんで、同じ場所で振動する感じです。
速度を変更せず、更に半径を少し縮めます。動ける範囲は少し広がるけど、やっぱりガチガチです。
どんどん小さくしていくと、ある所で急にブワーっとグチャグチャになるはずです。これが固体・液体の相転移。
これを実際の高分子の玉でやったのがこの実験
http://physicsworld.com/cws/article/news/2012/oct/05/tiny-spheres-simulate-crystal-melting

実験２

液体の状態から、逆に玉の半径を徐々に増やしていきます。減らすのより若干別の処理がいるかも。どこかで、かたまって固体になってガチガチになるはず。

実験３

液体の状態で、玉の半径を若干ばらつかせて、大きいのや小さいのがあるようにします。ここから玉の半径を、互いの比率は保ったままで徐々に大きくしていきます。どこかでガチガチになるけど、きれいに揃えない。これがガラス。
動画を発見。固まるところまでいってませんが。