分子動力学法ステップ・バイ・ステップ その2

はじめに

Step1では、O(N^2)のコードを作った。今回はO(N)に落とす準備として、ペアリストを導入する。

今回のコードはgithub/kaityo256/mdstep/step2においてある。

• その1　 $O(N^2)$のルーチンまで
• その2　 ペアリストの構築とBookkeeping法 ← イマココ
• その3　 メッシュ探索
• その4 　少しだけ高速化
• その5 　温度制御法を三種類実装してみる
• その6 　圧力測定ルーチンの実装

ペアリスト

短距離相互作用をする系では、相互作用が届かないところまで力の計算をする意味はない。そこで、カットオフ距離を導入し、その距離より遠いところでは相互作用を計算しないことで計算の加速をはかる。

そのために、まずどの原子が相互作用をしているか、相互作用ペアリストを作る必要がある。最終的にはこのペアリストもO(N)で作ることにするが、まずはペアリストをO(N^2)で作ることにしよう。

C++にはペアを扱うためのテンプレートも用意されているが、とりあえず自分で作ってしまうことにしよう。二体相互作用をする原子は、よく「i粒子/原子」「j粒子/原子」という呼び方をするので、インデックスiとjを持つ構造体を作る。

struct Pair{
  int i,j;
};

ペアリストを管理するクラスを作っても良いが、今回はMDクラスに突っ込んでしまうことにしよう。ペアリストを管理するためのvectorも定義する。

  std::vector<Pair> pairs;

まず、ペアを作る関数、make_pairを実装する。といってもO(N^2)なら作るのは簡単。

void
MD::make_pair(void){
  pairs.clear();
  const int pn = vars->number_of_atoms();
  Atom *atoms = &(vars->atoms[0]);
  for (int i = 0; i < pn - 1; i++) {
    for (int j = i + 1; j < pn; j++) {
      double dx = atoms[j].qx - atoms[i].qx;
      double dy = atoms[j].qy - atoms[i].qy;
      double dz = atoms[j].qz - atoms[i].qz;
      adjust_periodic(dx, dy, dz);
      double r2 = (dx * dx + dy * dy + dz * dz);
      if (r2 > CL2)continue;
      Pair p;
      p.i = i;
      p.j = j;
      pairs.push_back(p);
    }
  }
}

そして、作ったペアリストから力を計算する、calculate_force_pairも実装する。

void
MD::calculate_force_pair(void) {
  const int pp = pairs.size();
  Atom *atoms = &(vars->atoms[0]);
  for(int k=0;k<pp;k++){
    const int i = pairs[k].i;
    const int j = pairs[k].j;
    double dx = atoms[j].qx - atoms[i].qx;
    double dy = atoms[j].qy - atoms[i].qy;
    double dz = atoms[j].qz - atoms[i].qz;
    adjust_periodic(dx, dy, dz);
    double r2 = (dx * dx + dy * dy + dz * dz);
    if (r2 > CL2)continue;
    double r6 = r2 * r2 * r2;
    double df = (24.0 * r6 - 48.0) / (r6 * r6 * r2) * dt;
    atoms[i].px += df * dx;
    atoms[i].py += df * dy;
    atoms[i].pz += df * dz;
    atoms[j].px -= df * dx;
    atoms[j].py -= df * dy;
    atoms[j].pz -= df * dz;
  }
}

二重ループが、ペアに対する一重ループになる。

まずは、ペアリストを使うバージョンのテストをしよう。具体的には、MD::calculateをこう書き換える。

void
MD::calculate(void) {
  update_position();
  make_pair();
  calculate_force_pair();
  update_position();
  periodic();
  vars->time += dt;
}

実行して、step1と結果が変わらないことを確認する(重要)。

$ time ./a.out > energy.dat
./a.out > energy.dat  23.11s user 0.12s system 99% cpu 23.246 total

$ diff energy.dat ../step1/energy.dat

実行時間は手元のPCで23秒。もともとのコードが22秒なので、やや遅くなった。

Bookkeeping法

さて、ペアリストを導入しただけでは速度があがらない。ペアリストを使う「うまみ」は、Bookkeeping法にある。Bookkeeping法とは、カットオフにマージンを設け、カットオフよりも長い距離でペアを登録して、ある程度の時間、そのペアリストを使い回す方法である。今回はカットオフが2.0なので、マージンを0.5くらいにしよう。後で使うので、マージンと、マージンのカットオフの和の平方をsystemparam.hppに定義しておく。

const double MARGIN = 0.5;
const double ML2 = (CUTOFF+MARGIN)*(CUTOFF+MARGIN);

ペアリスト構築のところで、CL2と書いてあったところをML2とすればマージン付きでペアリストを構築したことになる。問題はどれくらい安全にペアリストを使いまわせるかだが、これを「系の最速の原子が、マージンを半分通過するまで」とする。半分にするのは、最速の原子と、二番目に速い原子が正面衝突をする場合が最悪の場合、つまりペアリストの寿命が最も短い場合だからである。

さて、ペアリストの「寿命」を長さで表現する。これをmargin_lengthとして、MDクラスのメンバに入れておこう。ペアリストを作った直後は、margin_length = MARGINである。しかし、時間が経過するにつれて、マージンは減っていく。減る量は「系最速の原子の速度×2×タイムステップ」である。従って、以下のようなチェックルーチンを作る。

void
MD::check_pairlist(void){
  double vmax2 = 0.0;
  for (auto &a : vars->atoms) {
    double v2 = a.px*a.px + a.py*a.py + a.pz*a.pz;
    if (vmax2 < v2) vmax2 = v2;
  }
  double vmax = sqrt(vmax2);
  margin_length -= vmax*2.0*dt;
  if(margin_length < 0.0){
    margin_length = MARGIN;
    make_pair();
  }
}

つまり、系最速の原子を探し、その速度vmaxを求め、margin_lengthをvmax2.0dtだけ減らし、もしmargin_lengthを使い切ったらペアリストを作るルーチンである。作りなおす際にペアリスト寿命(マージン)を回復させている。

計算ルーチンでは、make_pairの代わりにcheck_pairlistを呼べば良い。

void
MD::calculate(void) {
  update_position();
  check_pairlist();
  calculate_force_pair();
  update_position();
  periodic();
  vars->time += dt;
}

また実行して、step1と結果が変わらないことを確認する(重要)。

$ time ./a.out > energy.dat
./a.out > energy.dat  5.62s user 0.11s system 99% cpu 5.736 total

$ diff energy.dat ../step1/energy.dat

結果は全く変わらず、速度が4倍になった。

観測ルーチンのO(N)化

せっかくペアリストを導入したので、観測ルーチンもペアリストを使うように修正しよう。Observer::potential_energyの引数に、 std::vector<Pair> &pairsを追加して、力の計算と同様にペアリストでなめるようにする。

double
Observer::potential_energy(Variables *vars, std::vector<Pair> &pairs) {
  double v = 0.0;
  const int pp = pairs.size();
  const int pn = vars->number_of_atoms();
  Atom *atoms = vars->atoms.data();
  for(int k=0;k<pp;k++){
    const int i = pairs[k].i;
    const int j = pairs[k].j;
    double dx = atoms[j].qx - atoms[i].qx;
    double dy = atoms[j].qy - atoms[i].qy;
    double dz = atoms[j].qz - atoms[i].qz;
    adjust_periodic(dx, dy, dz);
    double r2 = (dx * dx + dy * dy + dz * dz);
    if (r2 > CL2)continue;
    double r6 = r2 * r2 * r2;
    double r12 = r6 * r6;
    v += 4.0 * (1.0 / r12 - 1.0 / r6) + C0;
   }
  v /= static_cast<double>(pn);
  return v;
}

なお、シミュレーションの最初で、ペアリストを作っておかなければ正しくポテンシャルエネルギーが評価されないので、MD::runで、初期条件を作った直後にペアリストを作っておこう。

void
MD::run(void) {
  makeconf();
  make_pair(); // ←追加

これで、ペアリスト構築以外は全てO(N)となった。最終的な実行速度は

$ time ./a.out > energy.dat
./a.out > energy.dat  5.45s user 0.11s system 99% cpu 5.570 total

となった。

まとめ

今回はペアリストとBookkeeping法を導入した。ここまでのコードはgithub/kaityo256/mdstep/step2においてある。

ペアリストを作る際、カットオフ距離よりも長い距離のペアを登録してペアリストをしばらく使い回すことで、ペアリストを構築する手間が減る。ペアリストの構築の手間はO(N^2)だが、その寿命チェックと、ペアリストを用いた力の計算はO(N)なので、その効果は500粒子程度でもかなり大きい。

なお、マージンを長く取れば取るほどペアリストの寿命が伸びるが、その分ペアリスト構築の手間と、力の計算で相互作用していないペアを扱う時間が増えるので、あまり長くすると逆効果となる。最適なマージンは系の密度、温度、時間刻み等に依存するが、経験的にカットオフ距離の10%程度、もしくはペアリストの使い回しが10回程度になるあたりがちょうど良いことが多い気がする。

原子数が少ない場合はこれでよいが、さらに大きな系を計算したい場合はメッシュ探索を用いるなどしてペアリスト作成の手間をO(N)に落とす必要がある。