はじめに

Step 3まででほぼ実用的なMDが書けた。チューニングはしないと書いたが、少しだけ高速化を考える。

今回のコードはgithub/kaityo256/mdstep/step4においてあるので、適宜参照されたい。

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

メッシュ探索

周期的境界条件の補正

perfなどのツールで見ると、adjust_periodicという関数が結構な割合を占めている。systemparam.cppに実体が書いてあるため、インライン展開されなかったようだ。これをsystemparam.hppに実体を書き、inline指定しよう。

systemparam.hpp

inline void adjust_periodic(double &dx, double &dy, double &dz){
  const double LH = L * 0.5;
  if (dx < -LH)dx += L;
  if (dx > LH) dx -= L;
  if (dy < -LH)dy += L;
  if (dy > LH) dy -= L;
  if (dz < -LH)dz += L;
  if (dz > LH) dz -= L;
}

systemparam.cppに書いてあった実体は消す。

# インライン指定前
$ time ./a.out > energy.dat
./a.out > energy.dat  5.26s user 0.15s system 99% cpu 5.433 total

# インライン指定後
$ time ./a.out > energy.dat
./a.out > energy.dat  4.34s user 0.14s system 99% cpu 4.494 total

5.26秒→4.34秒ということで、わりと早くなった。もちろん、周期境界条件については、もっとかしこく扱う方法があるが、とりあえずここではこれくらいで良いことにしよう。

相互作用相手によるソート

現在の力の計算は、ペアごとになっている。つまり、

まずk番目のペア情報を取って来る
ペア情報からi,jのインデックスを取って来る
i原子、j原子の運動量と座標の情報を取ってくる
力を計算する
i原子、j原子の運動量を書き戻す

という形になっている。データ構造を工夫することで、メモリアクセスの回数を大幅に減らすことができる。具体的には、i原子と相互作用しているj原子の情報を一箇所にまとめておいて、i原子を固定してその相互作用相手でループを回すことにする。以下のような二重ループになる。

まずi原子の情報を取ってくる
i原子と相互作用しているj原子の情報を取ってくる
力の計算をする
j原子の運動量を書き戻す。i原子にかかる力積はキャッシュしておく
まだ相互作用相手がいるなら2へ
i原子の力積を書き戻し、iをインクリメントして1へ

こうすると、i原子に関する情報がレジスタに乗り、相互作用相手の数が十分に多い場合は、i原子に関する情報の読み書きが無視できるため、メモリアクセスが半分になる。

上記の力計算を行うためには、ペアリストができた段階で、相互作用相手によるソートをしておけば良い。本当はペアリストを作りながら、相互作用相手によるソートの準備をしておくほうが効率が良いのだが、わかりやすさのために、既に完成したペアリストから相互作用相手のソートを行うことにする。

その前に、ペアリストpairsで、i原子のインデックスがj原子のインデックスより小さいようにしておかないといろいろ面倒なので、ペアリスト登録の際にチェックしておく。

      Pair p;
      if (i < j) {
        p.i = i;
        p.j = j;
      } else {
        p.i = j;
        p.j = i;
      }
      pairs.push_back(p);

MeshList::searchとMeshList::search_otherの二箇所あるので注意。

ペアリストを受け取ってソートするためのメソッド、Variables::make_neighbor_list(std::vector<Pair> &pairs)を実装しよう。さほど長くないのでメソッド全体を書いてしまう。std::vectorのメモリの確保/解放にかなり無頓着な書き方をしているが、よほどmallocが腐ってなければ速度に影響は無いと思う。無論、ここのmalloc/freeが見えるようであれば、何度もメモリの確保/解放をしないように書き換えなければならない。

void
Variables::make_neighbor_list(std::vector<Pair> &pairs) {
  const int pn = atoms.size();
  const int pp = pairs.size();
  neighbor_list.clear();
  neighbor_list.resize(pp);
  i_position.resize(pn);
  j_count.resize(pn);
  std::fill(j_count.begin(), j_count.end(), 0);
  for (auto &p : pairs) {
    j_count[p.i]++;
  }
  i_position[0] = 0;
  int sum = 0;
  for (int i = 0; i < pn - 1; i++) {
    sum += j_count[i];
    i_position[i + 1] = sum;
  }
  std::vector<int> pointer(pn);
  std::fill(pointer.begin(), pointer.end(), 0);
  for (auto &p : pairs) {
    int pos = i_position[p.i] + pointer[p.i];
    neighbor_list[pos] = p.j;
    pointer[p.i]++;
  }
}

後で気が向いたら図解するかもしれないが、やっていることは住所録作成と同じなので説明不要と思う。neighbor_listというのがソートしたj原子の情報で、ペアリストの長さの配列である。i_positionというのは、neighbor_listのどこに自分の相互作用相手がいるかのインデックスで、j_countが、相互作用相手が何個いるかの数である。例えばi原子の相互作用相手は、neighbor_listのi_position[i]番目から、j_count[i]個あることになる。

念のため、作成したソート情報が、元のペアリストを正しく反映しているか確認しよう。
もとのペアリストの情報は

    for(auto &p:pairs){
      printf("%03d %03d\n",p.i,p.j);
    }

と吐き、ソートした相互作用相手の情報を

    for(int i=0;i<pn;i++){
      for(int n=0;n<j_count[i];n++){
        int pos = i_position[i] + n;
        assert(pos < pp);
        int j = neighbor_list[pos];
        printf("%03d %03d\n",i,j);
      }

として吐く。これをsort書けてdiffを取り、違いがなければ正しく作成できている。

ペアリストを構築したら、このソートをしておかないといけない。MeshList::make_pairメソッドの最後に

  vars->make_neighbor_list(pairs);

と書いておこう。

この情報から力の計算は以下のように書ける。

void
MD::calculate_force_list(void) {
  Atom *atoms = vars->atoms.data();
  const int pn = vars->number_of_atoms();
  int *neighbor_list = vars->neighbor_list.data();
  int *i_position = vars->i_position.data();
  int *j_count = vars->j_count.data();
  for (int i = 0; i < pn; i++) {
    const double qix = atoms[i].qx;
    const double qiy = atoms[i].qy;
    const double qiz = atoms[i].qz;
    double pix = atoms[i].px;
    double piy = atoms[i].py;
    double piz = atoms[i].pz;
    const int ip = i_position[i];
    for (int k = 0; k < j_count[i]; k++) {
      const int j = neighbor_list[ip + k];
      double dx = atoms[j].qx - qix;
      double dy = atoms[j].qy - qiy;
      double dz = atoms[j].qz - qiz;
      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;
      pix += df * dx;
      piy += df * dy;
      piz += df * dz;
      atoms[j].px -= df * dx;
      atoms[j].py -= df * dy;
      atoms[j].pz -= df * dz;
    }
    atoms[i].px = pix;
    atoms[i].py = piy;
    atoms[i].pz = piz;
  }
}

i原子の情報が内側のループで変わらないことをわかりやすくするためにconstを使っているが、まっとうなコンパイラならこれが内側ループで変化しないことを見抜くと思われる。また、一時変数を使ってi原子の力積の足し込みをしておいて、最後に書き戻している。これは明示的に行わないとやってくれないコンパイラの方が多いと思う。

さて、この書き換えの効果だが、実際にはシングルプロセスではほとんど影響が無い。

$ time ./a.out > energy.dat
./a.out > energy.dat  4.42s user 0.13s system 99% cpu 4.569 total

むしろ余計な手間をかけた分、若干遅くなっている。これは、原子数が500と少なく全ての情報がキャッシュに乗ってしまい、かつカットオフが短いため、相互作用相手が少ないためである。手元のPCでは、この最適化の影響はほとんど無いが、別のIntel(R) Xeon(R) CPU E5-2680なマシンで、インテルコンパイラでビルドした結果では、

const double L = 100;
const double dt = 0.001;
const double CUTOFF = 3.5;
const double MARGIN = 0.5;
const int STEPS = 100;

といった条件、つまり50万原子をカットオフ3.5、時間刻み0.001で100ステップだけ時間発展させると、このソートしたほうが10%程度速かった。

この高速化が効くのは、原子の情報がキャッシュから溢れており、メモリアクセスをしている時である。一つの原子の情報は、6つのdouble、すなわち48バイトあるが、最近の石のL3なら20〜30MB程度はあるだろう。すると50万原子の情報程度ならキャッシュに乗ってしまう。従って、それより小さい原子数を扱う場合にはこの最適化の効果は薄い。

まとめ

一応少しだけ最適化してみた。前者はともかく、後者はシングルプロセスではほとんど恩恵が無いため、普通は気にしなくて良いと思う。後者が効いてくるのは、キャッシュに入りきらないくらいの原子数を扱う場合で、並列化などをしていると効いてくるが、それでも速度向上は10〜20%程度である。並列化して何億もするスパコンをぶん回すなら、10%〜20%の高速化は気にすべきであろうが、ローカルマシンを使ってる場合や、シングルコアの場合はあまり気にしなくて良いと思う。

本気で高速化するなら、Step 1でも述べたようにデータ構造をAoSにするかSoAにするかといったところから意識して、最初から全部組み直したほうが良いと思う。~~あとstd::vectorを全く使わないとかね。~~ただここで述べたかったことは、最近の計算機は演算能力に対してメモリアクセスがボトルネックになりやすいため、余計な計算をしてもメモリアクセスを減らすという最適化が有効であることが多い、ということ。

ここまでのコードはgithub/kaityo256/mdstep/step4においてある。