Lennard-Jonesポテンシャルを用いたMDシミュレーション

ターム指向記述とアトム指向記述

原子間に働く力を計算している部分

  for(i=0; i<N; ++i){
    for(j=i+1; j<N; j++){
      dxij = rx[i] - rx[j];
      dyij = ry[i] - ry[j];
      dzij = rz[i] - rz[j];
      rij  = sqrt(dxij * dxij + dyij * dyij + dzij * dzij);
      fij  = - 4.0 * epsilon *(-12.0 * pow(sigma/rij, 12.0) / rij
             + 6.0 * pow(sigma/rij, 6.0) / rij);
      fx[i] += fij * dxij / rij;
      fy[i] += fij * dyij / rij;
      fz[i] += fij * dzij / rij;
      fx[j] -= fij * dxij / rij;
      fy[j] -= fij * dyij / rij;
      fz[j] -= fij * dzij / rij;
      V += 4.0*epsilon*(pow(sigma/rij,12.0)-pow(sigma/rij,6.0));
   }
 }

は、次のように短く記述することもできる。

  for(i=0; i<N; ++i){
    for(j=0; j<N; j++){
      dxij = rx[i] - rx[j];
      dyij = ry[i] - ry[j];
      dzij = rz[i] - rz[j];
      rij  = sqrt(dxij * dxij + dyij * dyij + dzij * dzij);
      fij  = - 4.0 * epsilon *(-12.0 * pow(sigma/rij, 12.0) / rij
             + 6.0 * pow(sigma/rij, 6.0) / rij);
      fx[i] += fij * dxij / rij;
      fy[i] += fij * dyij / rij;
      fz[i] += fij * dzij / rij;
      if(i<j) V += 4.0*epsilon*(pow(sigma/rij,12.0)-pow(sigma/rij,6.0));
   }
 }

筆者は、前者の書き方を「ターム指向」、後者の書き方を「アトム指向」と呼んで区別している。Lennard-Jonesポテンシャルのように単純なポテンシャルでは大した違いはないが、配位数に依存する多体相互作用モデルや多元素混在系になると、実は前者の方が圧倒的に有利になる。プログラミング上の間違いが少なく、パラメータを格納しているメモリへのアクセス回数が減り、計算速度も向上するからである。

ターム指向記述とは、ポテンシャルを構成する各項(ターム)に注目し、各項に係わるすべての原子に対して、その項から生じる力を配分していく書き方である。ポテンシャルエネルギーの定義式のΣに従って大枠のforループを構成すればよい。

一方、アトム指向記述では、一つ一つの原子に注目し、その原子に加わる力を列挙して足し合わせていく。一番外側のforループは、系の全ての原子をスキャンするためのループであり、その内側のループは、各原子に加わる力の源となるポテンシャルの項をスキャンする。この方式では同じ項を何度も参照するので、パラメータの値をロードする回数が増え、効率が悪い。ポテンシャルエネルギーを足し合わせる条件にも気を配る必要があり、ミスの温床となる。

行数が少ないからといって、洗練された良いプログラムとは限らない。ポテンシャルエネルギーと力の計算ルーチンをインプリメントする際は、ターム指向記述を心がけよう。