kivantium活動日記

プログラムを使っていろいろやります

NWChemをコンパイルして動かす

なにも分からない。

NWChemは計算化学ソフトで、GSoCのOpenChemistryに参加していたので動かしてみたかった。

とりあえず動いた方法

環境はUbuntu 18.04。

ドキュメントのCompiling NWChemにあるNWChem 6.6 on Ubuntu 14.04 (Trusty Tahr)の項目を見ればいいと思っていたが微妙に動かなかった。そもそも現在の最新リリースバージョンは6.8.1なのでドキュメントが更新されてなさそう。以下のように改変したらとりあえず動くようになった。

必要なパッケージのインストール

python-dev gfortran libopenblas-dev libopenmpi-dev openmpi-bin tcsh make 

最新リリースをクローンする

git clone -b hotfix/release-6-8 https://github.com/nwchemgit/nwchem.git nwchem-6.8.1
cd nwchem-6.8.1/src

環境変数の設定

ドキュメントでは

export USE_MPI=y  
export NWCHEM_TARGET=LINUX64  
export USE_PYTHONCONFIG=y  
export PYTHONVERSION=2.7  
export PYTHONHOME=/usr 
export BLASOPT="-lopenblas -lpthread -lrt"  
export BLAS_SIZE=4  
export USE_64TO32=y

とするように指示されているが、BLASOPTを指定するとmake時にLAPACK_LIBを指定しろというエラーが発生する。しかし、その正しい値についてはドキュメントに言及がなく分からなかった。そこで

unset BLASOPT

を実行したところとりあえずエラーが消えたので良しとした。(絶対良くないんだよなぁ)

コンパイル

make nwchem_config NWCHEM_MODULES="all python"
make 64_to_32
make

定石どおり make install しようとしたら動かなかったのでコンパイル後にできたbin/以下をPATHに加えてインストール終了とした。

動作チェック

Getting Started にあるコードが動くことを確認した。(aptでインストールしたバージョンではこれが動かなかったのがそもそもコンパイルしないといけなくなった原因なのでNWChemのメンテナンスが全般的に間に合ってないっぽい)

ついでに閉殻Hartree-Fock法によるHeH+のエネルギー計算 - kivantium活動日記と同じ設定で走らせて結果がどうなるか見てみた。

title "HeH+ STO-3G SCF energy calculation"
charge 1
geometry  
  H 0 0 0
  He 0 0 1.4632
end
basis
  H  library sto-3g
  He library sto-3g
end
task scf energy

結果は

       Final RHF  results 
       ------------------ 

         Total SCF energy =     -2.825194209471
      One-electron energy =     -4.549711559836
      Two-electron energy =      1.001202357200
 Nuclear repulsion energy =      0.723314993165

        Time for solution =      0.0s

となり、前回の結果のTOTAL ENERGY = -2.860662163500と近い値が出ることが確認できた。

閉殻Hartree-Fock法によるHeH+のエネルギー計算

量子化学計算を勉強するために新しい量子化学―電子構造の理論入門〈上〉を読んでいます。
サンプルコードがFortranだったので勉強がてらC言語に移植しました。

新しい量子化学―電子構造の理論入門〈上〉

新しい量子化学―電子構造の理論入門〈上〉

理論

本で200ページくらいかけて説明されている内容をブログに書くのは大変なので省略します。

ちなみに、以下のツイートの「100ページかけて説明してる別の本」はこの本を指します。
難しい内容を易しく簡潔に説明しようとするのは大変で、難しい内容は難しいまま長い時間をかけて説明する必要があるのだなぁという学びがありました。


実装

本のFortran実装を素直にC言語で書き直しました。
本のコードを打ち込むのが面倒だったので、検索して見つけたページにあったコードと見比べながら書きました。

最終的に得られたエネルギーはFortran実装の結果と小数第5位まで一致しました。
この差がどこで生じているのか謎ですが、実装は正しいんじゃないでしょうか……

#include <math.h>
#include <stdio.h>

// プロトタイプ宣言
void HFCALC(int, int, double, double, double, double, double);
void INTGRL(int, int, double, double, double, double, double);
void COLECT(int, int, double, double, double, double, double);
void SCF(int, int, double, double, double, double, double);
double calcS(double, double, double);
double calcT(double, double, double);
double F0(double);
double calcV(double, double, double, double, double);
double TWOE(double, double, double, double, double, double, double);
void FORMG(void);
void DIAG(double [][2], double [][2], double [][2]);
void MULT(double [][2], double [][2], double [][2], int, int);
void MATOUT(double [][2], int, int, int, int, char*);

double S12, T11, T12, T22, V11A, V12A, V22A, V11B, V12B, V22B,
       V1111, V2111, V2121, V2211, V2221, V2222;

// Slater型関数に対する最小2乗近似の短縮係数と軌道指数
double COEF[3][3] = {{1.0, 0.678914, 0.444635}, 
                     {0.0, 0.430129, 0.535328},
                     {0.0, 0.0, 0.154329}};
double EXPON[3][3] = {{0.27095, 0.151623, 0.109818},
                      {0.0, 0.851819, 0.405771},
                      {0.0, 0.0, 2.22766}};
double D1[3], A1[3], D2[3], A2[3];
double PI = 3.1415926535898;

double S[2][2], X[2][2], XT[2][2], H[2][2], F[2][2], G[2][2], C[2][2],
       FPRIME[2][2], CPRIME[2][2], P[2][2], OLDP[2][2], TT[2][2][2][2], E[2][2];

int main(void) {
  int IOP = 2;
  // 基底関数STO-NG
  int N = 3;
  // 基底状態の平衡結合長(a.u.)
  double R = 1.4632;
  // Heに対するSTO-3G軌道指数値
  // He原子に対する最良の軌道指数は1.6875
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA1 = 2.0925;
  // Hに対するSTO-3G軌道指数値
  // H原子に対する最良の軌道指数は1.0
  // HeH+では電子雲が縮んでいることを反映して1.24倍
  double ZETA2 = 1.24;
  double ZA = 2.0;
  double ZB = 1.0;
  HFCALC(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/*
 * 1s基底関数に対応するSTO-NG最小基底を用いて2電子系2原子分子に対する
 * Hartree-Fock計算を行う関数
 *
 * IOP=0 何も表示しない
 * IOP=1 収束の結果のみ表示
 * IOP=2 全繰り返しの結果を表示
 * N STO-NG を使う (N=1, 2, 3)
 * R 結合長 (a.u.)
 * ZETA1 原子Aの基底関数のSlater軌道指数値
 * ZETA2 原子Bの基底関数のSlater軌道指数値
 * ZA 原子Aの原子番号
 * ZB 原子Bの原子番号
 */
void HFCALC(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  if (IOP != 0) {
    printf("   STO-%dG for atomic numbers %5.2lf and %5.2lf\n\n\n", N, ZA, ZB);
  }
  INTGRL(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  COLECT(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
  SCF(IOP, N, R, ZETA1, ZETA2, ZA, ZB);
}

/* SCF計算に必要な積分を計算する関数*/
void INTGRL(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  double R2 = R * R;
  // 軌道指数のスケーリングと短縮係数の規格化
  for (int i = 0; i < N; ++i) {
    A1[i] = EXPON[i][N-1] * pow(ZETA1, 2.0);
    D1[i] = COEF[i][N-1] * pow(2.0 * A1[i] / PI, 0.75);
    A2[i] = EXPON[i][N-1] * pow(ZETA2, 2.0);
    D2[i] = COEF[i][N-1] * pow(2.0 * A2[i] / PI, 0.75);
  }
  S12 = 0.0;
  T11 = 0.0;
  T12 = 0.0;
  T22 = 0.0;
  V11A = 0.0;
  V12A = 0.0;
  V22A = 0.0;
  V11B = 0.0;
  V12B = 0.0;
  V22B = 0.0;
  V1111 = 0.0;
  V2111 = 0.0;
  V2121 = 0.0;
  V2211 = 0.0;
  V2221 = 0.0;
  V2222 = 0.0;

  // 1電子積分の計算
  // V12A = 原子核Aへの各引力の非対角項
  // RAP2 = 中心Aと中心Pの距離の2乗
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      double RAP = A2[j] * R / (A1[i] + A2[j]);
      double RAP2 = pow(RAP, 2.0);
      double RBP2 = pow(R - RAP, 2.0);
      S12 += calcS(A1[i], A2[j], R2) * D1[i] * D2[j];
      T11 += calcT(A1[i], A1[j], 0.0) * D1[i] * D1[j];
      T12 += calcT(A1[i], A2[j], R2) * D1[i] * D2[j];
      T22 += calcT(A2[i], A2[j], 0.0) * D2[i] * D2[j];
      V11A += calcV(A1[i], A1[j], 0.0, 0.0, ZA) * D1[i] * D1[j];
      V12A += calcV(A1[i], A2[j], R2, RAP2, ZA) * D1[i] * D2[j];
      V22A += calcV(A2[i], A2[j], 0.0, R2, ZA) * D2[i] * D2[j];
      V11B += calcV(A1[i], A1[j], 0.0, R2, ZB) * D1[i] * D1[j];
      V12B += calcV(A1[i], A2[j], R2, RBP2, ZB) * D1[i] * D2[j];
      V22B += calcV(A2[i], A2[j], 0.0, 0.0, ZB) * D2[i] * D2[j];
    }
  }
  // 2電子積分の計算
  for (int i = 0; i < N; ++i) {
    for (int j = 0; j < N; ++j) {
      for (int k = 0; k < N; ++k) {
        for (int l = 0; l < N; ++l) {
          double RAP = A2[i] * R / (A2[i] + A1[j]);
          //double RBP = R - RAP;
          double RAQ = A2[k] * R / (A2[k] + A1[l]);
          double RBQ = R - RAQ;
          double RPQ = RAP - RAQ;
          double RAP2 = RAP * RAP;
          //double RBP2 = RBP * RBP;
          //double RAQ2 = RAQ * RAQ;
          double RBQ2 = RBQ * RBQ;
          double RPQ2 = RPQ * RPQ;
          V1111 += TWOE(A1[i], A1[j], A1[k], A1[l], 0.0, 0.0, 0.0) 
                   * D1[i] * D1[j] * D1[k] * D1[l];
          V2111 += TWOE(A2[i], A1[j], A1[k], A1[l], R2, 0.0, RAP2)
                   * D2[i] * D1[j] * D1[k] * D1[l];
          V2121 += TWOE(A2[i], A1[j], A2[k], A1[l], R2, R2, RPQ2)
                   * D2[i] * D1[j] * D2[k] * D1[l];
          V2211 += TWOE(A2[i], A2[j], A1[k], A1[l], 0.0, 0.0, R2)
                   * D2[i] * D2[j] * D1[k] * D1[l];
          V2221 += TWOE(A2[i], A2[j], A2[k], A1[l], 0.0, R2, RBQ2)
                   * D2[i] * D2[j] * D2[k] * D1[l];
          V2222 += TWOE(A2[i], A2[j], A2[k], A2[l], 0.0, 0.0, 0.0) * D2[i] *
                   D2[j] * D2[k] * D2[l];
        }
      }
    }
  }
  if (IOP != 0) {
    printf("   R          ZETA1      ZETA2      S12       T11\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", R, ZETA1, ZETA2, S12, T11);
    printf("   T12        T22        V11A       V12A        V22A\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", T12, T22, V11A, V12A, V22A);
    printf("   V11B       V12B       V22B       V1111        V2111\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf%11.6lf\n\n\n", V11B, V12B, V22B, V1111,
           V2111);
    printf("   V2121      V2211       V2221     V2222\n\n");
    printf("%11.6lf%11.6lf%11.6lf%11.6lf\n", V2121, V2211, V2221, V2222);
  }
}

// 重なり積分
double calcS(double A, double B, double RAB2) {
  return pow(PI / (A + B), 1.5) * exp(-A * B * RAB2 / (A + B));
}

// 運動エネルギー積分
double calcT(double A, double B, double RAB2) {
  return A * B / (A + B) * (3.0 - 2.0 * A * B * RAB2 / (A + B)) *
         pow(PI / (A + B), 1.5) * exp(-A*B*RAB2/(A+B));
}

// F0関数
double F0(double ARG) {
  double ret;
  if (ARG < 1.0e-6) {
    ret = 1.0 - ARG / 3.0;
  } else {
    ret = sqrt(PI / ARG) * erf(sqrt(ARG)) / 2.0;
  }
  return ret;
}

// 核引力積分
double calcV(double A, double B, double RAB2, double RCP2, double ZC) {
  double v = 2.0 * PI / (A + B) * F0((A + B) * RCP2) * exp(-A * B * RAB2 / (A + B));
  return -v * ZC;
}

// 2電子積分
double TWOE(double A, double B, double C, double D, double RAB2, double RCD2,
            double RPQ2) {
  return 2.0 * pow(PI, 2.5) / ((A + B) * (C + D) * sqrt(A + B + C + D))
         * F0((A + B) * (C + D) * RPQ2 / (A + B + C + D))
         * exp(-A * B * RAB2 / (A + B) - C * D * RCD2 / (C + D));
}

// 積分の結果から行列を求める
void COLECT(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
            double ZB) {
  // 核ハミルトニアン cf. 式(3.153)
  H[0][0] = T11 + V11A + V11B;
  H[0][1] = T12 + V12A + V12B;
  H[1][0] = H[0][1];
  H[1][1] = T22 + V22A + V22B;

  // 重なり行列 cf. 式(3.136)
  S[0][0] = 1.0;
  S[0][1] = S12;
  S[1][0] = S[0][1];
  S[1][1] = 1.0;

  // 正準直交化 cf. 式(3.169), p.191
  X[0][0] = 1.0 / sqrt(2.0 * (1.0 + S12));
  X[1][0] = X[0][0];
  X[0][1] = 1.0 / sqrt(2.0 * (1.0 - S12));
  X[1][1] = -X[0][1];

  // 変換行列の転置
  XT[0][0] = X[0][0];
  XT[0][1] = X[1][0];
  XT[1][0] = X[0][1];
  XT[1][1] = X[1][1];

  // 2電子積分の行列
  TT[0][0][0][0] = V1111;
  TT[1][0][0][0] = V2111;
  TT[0][1][0][0] = V2111;
  TT[0][0][1][0] = V2111;
  TT[0][0][0][1] = V2111;
  TT[1][0][1][0] = V2121;
  TT[0][1][1][0] = V2121;
  TT[1][0][0][1] = V2121;
  TT[0][1][0][1] = V2121;
  TT[1][1][0][0] = V2211;
  TT[0][0][1][1] = V2211;
  TT[1][1][1][0] = V2221;
  TT[1][1][0][1] = V2221;
  TT[1][0][1][1] = V2221;
  TT[0][1][1][1] = V2221;
  TT[1][1][1][1] = V2222;

  if (IOP != 0) {
    MATOUT(S, 2, 2, 2, 2, "S   ");
    MATOUT(X, 2, 2, 2, 2, "X   ");
    MATOUT(H, 2, 2, 2, 2, "H   ");
    printf("\n");
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        for (int k = 0; k < 2; ++k) {
          for (int l = 0; l < 2; ++l) {
            printf("   (%2d%2d%2d%2d )%10.6lf\n", i+1, j+1, k+1, l+1, TT[i][j][k][l]);
          }
        }
      }
    }
  }
}

// SCFを実行する
void SCF(int IOP, int N, double R, double ZETA1, double ZETA2, double ZA,
         double ZB) {
  // 密度行列の収束規準
  double CRIT = 1.0e-4;
  // 繰り返しの最大回数
  int MAXIT = 25;
  // 核ハミルトニアンを最初の推定に用いる(つまりPを零行列で初期化)
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      P[i][j] = 0.0;
    }
  }
  if (IOP == 2) {
    MATOUT(P, 2, 2, 2, 2, "P    ");
  }
  int ITER = 0;
  double EN;
  for (;;) {
    ITER++;
    if (IOP == 2) {
      printf("\n    START OF ITERATION NUMBER = %2d\n", ITER);
    }
    // Fock行列の2電子部分を計算
    FORMG();

    if (IOP == 2) {
      MATOUT(G, 2, 2, 2, 2, "G    ");
    }

    // 核ハミルトニアンを足してFock行列を得る
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        F[i][j] = H[i][j] + G[i][j];
      }
    }

    // 電子エネルギーを計算する
    EN = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        EN += 0.5 * P[i][j] * (H[i][j] + F[i][j]);
      }
    }

    if (IOP == 2) {
      MATOUT(F, 2, 2, 2, 2, "F   ");
      printf("\n\n\n    ELECTRONIC ENERGY = %20.12lf\n", EN);
    }

    // Gを用いてFock行列を変換
    MULT(F, X, G, 2, 2);
    MULT(XT, G, FPRIME, 2, 2);
    // Fock行列の対角化
    DIAG(FPRIME, CPRIME, E);
    // 固有ベクトルから行列Cを得る
    MULT(X, CPRIME, C, 2, 2);
    // 新しい密度行列を計算する
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        OLDP[i][j] = P[i][j];
        P[i][j] = 0.0;
        for (int k = 0; k < 1; ++k) { // ここはk<1で本当にいいのか不明
          P[i][j] += 2.0 * C[i][k] * C[j][k];
        }
      }
    }
    if (IOP == 2) {
      MATOUT(FPRIME, 2, 2, 2, 2, "F'   ");
      MATOUT(CPRIME, 2, 2, 2, 2, "C'   ");
      MATOUT(E, 2, 2, 2, 2, "E   ");
      MATOUT(C, 2, 2, 2, 2, "C   ");
      MATOUT(P, 2, 2, 2, 2, "P   ");
    }
    // deltaを計算する
    double DELTA = 0.0;
    for (int i = 0; i < 2; ++i) {
      for (int j = 0; j < 2; ++j) {
        DELTA += pow(P[i][j] - OLDP[i][j], 2.0);
      }
    }
    DELTA = sqrt(DELTA / 4.0);
    if (IOP != 0) {
      printf("\n    DELTA(CONVERGENCE OF DENSITY MATRIX) = %10.6lf\n", DELTA);
    }

    // 収束判定
    if (DELTA < CRIT) break;
    if (ITER < MAXIT) continue;
    // 収束しなかった場合
    printf("    NO CONVERGENCE IN SCF");
    return;
  }

  double ENT = EN + ZA * ZB / R;
  if (IOP != 0) {
    printf("\n\n    CALCULATION CONVERGED"
           "\n\n    ELECTORONIC ENERGY = %20.12lf"
           "\n\n    TOTAL ENERGY =       %20.12f", EN, ENT);
  }
  if (IOP == 1) {
    MATOUT(G, 2, 2, 2, 2, "G   ");
    MATOUT(F, 2, 2, 2, 2, "F   ");
    MATOUT(E, 2, 2, 2, 2, "E   ");
    MATOUT(C, 2, 2, 2, 2, "C   ");
    MATOUT(P, 2, 2, 2, 2, "P   ");
  }
  MULT(P, S, OLDP, 2, 2);
  if (IOP != 0) {
    MATOUT(OLDP, 2, 2, 2, 2, "PS  ");
  }
}

void FORMG(void) {
  for (int i = 0; i < 2; ++i) {
    for (int j = 0; j < 2; ++j) {
      G[i][j] = 0.0;
      for (int k = 0; k < 2; ++k) {
        for (int l = 0; l < 2; ++l) {
          G[i][j] += P[k][l] * (TT[i][j][k][l] - 0.5 * TT[i][l][k][j]);
        }
      }
    }
  }
}

// Fを対角化し、Cに固有ベクトル、Eに固有値を入れる
void DIAG(double F[][2], double C[][2], double E[][2]) {
  double theta;
  if (fabs(F[0][0] - F[1][1]) > 1.0e-20) {
    theta = 0.5 * atan(2.0 * F[0][1] / (F[0][0] - F[1][1]));
  } else {
    theta = PI / 4.0;
  }
  C[0][0] = cos(theta);
  C[1][0] = sin(theta);
  C[0][1] = sin(theta);
  C[1][1] = -cos(theta);
  E[0][0] = F[0][0] * pow(cos(theta), 2.0) + F[1][1] * pow(sin(theta), 2.0)
            + F[0][1] * sin(2.0 * theta);
  E[1][1] = F[1][1] * pow(cos(theta), 2.0) + F[0][0] * pow(sin(theta), 2.0)
            - F[0][1] * sin(2.0 * theta);
  E[1][0] = 0.0;
  E[0][1] = 0.0;
  // 固有値と固有ベクトルを整列する
  if (E[1][1] > E[0][0]) return;
  double temp = E[1][1];
  E[1][1] = E[0][0];
  E[0][0] = temp;
  temp = C[0][1];
  C[0][1] = C[0][0];
  C[0][0] = temp;
  temp = C[1][1];
  C[1][1] = C[1][0];
  C[1][0] = temp;
}

// AとBの積をCに入れる
void MULT(double A[][2], double B[][2], double C[][2], int IM, int M) {
  for (int i = 0; i < M; ++i) {
    for (int j = 0; j < M; ++j) {
      C[i][j] = 0.0;
      for (int k = 0; k < M; ++k) {
        C[i][j] += A[i][k] * B[k][j];
      }
    }
  }
}

void MATOUT(double A[][2], int IM, int IN, int M, int N, char *LABEL) {
  int IHIGH = 0;
  for (;;) {
    int LOW = IHIGH + 1;
    IHIGH = IHIGH + 5;
    IHIGH = (IHIGH < N) ? IHIGH : N;
    printf("\n\n\n    THE %s ARRAY\n", LABEL);
    printf("               ");
    for (int i = LOW; i <= IHIGH; ++i) {
      printf("          %3d      ", i);
    }
    printf("\n");
    for (int i = 1; i <= M; ++i) {
      printf("%10d     ", i);
      for (int j = LOW; j <= IHIGH; ++j) {
        printf(" %18.10lf", A[i - 1][j - 1]);
      }
      printf("\n");
    }
    if (N <= IHIGH) break;
  }
}

最終的な出力は

    CALCULATION CONVERGED

    ELECTORONIC ENERGY =      -4.227529304014

    TOTAL ENERGY =            -2.860662163500


    THE PS   ARRAY
                           1                  2      
         1            1.5296370381       1.1199277981
         2            0.6424383868       0.4703629619

でした。参照したFortranコードの出力は

    CALCULATION CONVERGED

    ELECTRONIC ENERGY =  -0.422752913203D+01

    TOTAL ENERGY =       -0.286066199152D+01



    THE PS   ARRAY
                           1                  2
         1        0.1529637121D+01   0.1119927796D+01
         2        0.6424383099D+00   0.4703628793D+00

でした。

FortranからCに移植するときのtips

詰まったことを挙げておきます。

配列の添字

Fortranでは添字が1から始まるので配列にアクセスするときに添字をうまく調整する必要があります。

多次元配列の並び順

このページに詳しく書いてありますが、C言語の多次元配列はメモリ上でa[0][0], a[0][1], a[1][0], a[1][1], a[2][0], a[2][1]のような順番で並びますが、Fortranではa(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2)のような順番に並びます。多次元配列の初期化を書き直すときに気をつける必要があります。

整合配列

Fortranではサブルーチンに配列を渡した配列のサイズをサブルーチン側で定義することができます。(整合配列という)

      dimension a(6)
      a(1) = 1
      a(2) = 2
      a(3) = 3
      a(4) = 4
      a(5) = 5
      a(6) = 6
      call sub(a,3,2)
      end

      subroutine sub(a,n,m)
      dimension a(n,m) <---- a(n,*) でも同じ結果を得る
      do i=1,n
         write(6,*) (a(i,j),j=1,m)
      end do
      return
      end

コードはFortran プログラミングの基礎知識から引用しました。

C99で導入された可変長配列の機能を使うと同じことができるそうです。cf. 3.3 整合配列、可変長配列
今回は配列のサイズが決め打ちできたのでvoid MULT(double A[][2], double B[][2], double C[][2], int IM, int M) { のように関数を定義しました。

C言語のabs

特にFortranとは関係がないのですが、C言語のabsはint型に対して使う関数で、double型に対してはfabsという別の関数が用意されています。(これはC言語で関数オーバーロードができないことが原因です。C11で導入されたGenericという機能を使うとそれっぽいことが実現できるらしいです。cf. Generic - C言語入門

Docker上で開発を行う

ライブラリの開発をするときに、現在のバージョンと開発中のバージョンの性能比較をすることがあります。同じ名前のライブラリが1つの環境に複数存在すると面倒なことになるのでDockerを使って環境を分けることにしました。そのときに使ったコマンドについて記録します。

Dockerの仕組みやインストール方法については触れません。

Ubuntu 16.04のコンテナのダウンロード

Ubuntu 16.04を使うには

docker pull ubuntu:16.04

とします。

Dockerイメージの実行

Dockerイメージを入手したら

docker run -it ubuntu:16.04 /bin/bash

を実行するとコンテナ上でシェルが実行されます。オプションの意味は以下の通りです。

  • -i: 標準入力を開いたままにする
  • -t: 擬似ttyに接続する。

コンテナの実行中にCtrl+PCtrl+Qを続けて押すとデタッチされます。

アタッチするには、

docker ps

で現在実行中のコンテナのCONTAINER IDを確認してから

docker attach <CONTAINER ID>

を実行します。

別の実行方法

上の方法でコンテナを実行するとCtrl+Dをうっかり押したときにコンテナが終了してしまいます。以下のようにすると回避できます。

まずバックグラウンド実行のオプションをつけた状態でコンテナを起動します。

docker run -it -d --name mycontainer ubuntu:16.04

シェルを開くときには

docker exec -it mycontainer /bin/bash

のようにコンテナ上でbashを実行します。

docker stop mycontainer

とすればコンテナが停止します。

参考 - Dockerコンテナ内で操作 attachとexecの違い - Qiita

サーバーの変更

Ubuntuイメージの初期設定では海外サーバーになっているのでapt-getする前に国内サーバーに変更するとよいです。

sed -i.bak -e "s%http://[^ ]\+%http://ftp.jaist.ac.jp/pub/Linux/ubuntu/%g" /etc/apt/sources.list

参考: apt-getの利用リポジトリを日本サーバーに変更する - Qiita

Dockerイメージの作成

ライブラリの開発に必要なパッケージをインストールしたDockerイメージを用意しておくと便利です。こういうときにはDockerfileというファイルを作ります。

OpenBabel 2.4.1をmakeしてインストールしてあるDockerイメージを作るときに使ったDockerfileの例を示します。

FROM ubuntu:16.04
RUN apt-get update && apt-get install -y \
  build-essential \
  cmake \
  libcairo2-dev \
  libeigen2-dev \
  wget \
  zlib1g-dev
RUN wget https://github.com/openbabel/openbabel/archive/openbabel-2-4-1.tar.gz && \
  tar xf openbabel-2-4-1.tar.gz && \
  cd openbabel-openbabel-2-4-1/ && \
  mkdir build && cd build/ && \
  cmake .. && make -j ${nproc} && make install

これをDockerfileという名前で保存し、

docker build -t openbabel2.4.1 .

を実行すると openbabel2.4.1 という名前のDockerイメージが作成されます。

RDKitのDockerfileの場所: rdkit_containers/Dockerfile at master · rdkit/rdkit_containers · GitHub

ファイルのコピー

ファイルをコピーするにはホスト側で

docker cp hoge.txt <CONTAINER ID>:/huga

のようにすればいいです。

停止したコンテナの扱い

停止したものも含めたコンテナ一覧を見るには

docker ps -a

とします。

docker start <CONTAINER ID>

でコンテナを起動することができます。

docker rm <CONTAINER ID>

でコンテナを削除することができます。

docker ps -a -qでIDの一覧を出せるので、コンテナを一括で削除するには

docker rm `docker ps -a -q`

とすればよいです。(参考: Dockerイメージとコンテナの削除方法 - Qiita

イメージの削除

イメージ一覧は

docker images

で確認できます。

削除するには

docker rmi <イメージ名>

とすればよいです。

他に有用なコマンドを学んだら追記します。

Vue.jsとFlaskの間でフォーム情報をやり取りする

Pythonで作ったプログラムにWebインターフェースをつけることができると便利です。
せっかくなのでフロントエンドも勉強しようと思って見た目簡単そうなVue.jsを試しました。

ソースコード

バックエンド (Python)

PythonのシンプルなウェブアプリケーションフレームワークであるFlaskを使いました。
/にアクセスするとindex.htmlが表示されます。
/apiにPOSTでデータを送りつけると、送られたフォームデータをsplit()したものをjsonで返します。

from flask import Flask, request, render_template, jsonify
app = Flask(__name__)

@app.route('/')
def index():
    return render_template('index.html')

@app.route("/api", methods=['POST'])
def hello():
    result = {
        "results": request.json["text"].split()
    }
    return jsonify(result)

if __name__ == '__main__':
    app.run(debug=True)

フロントエンド (Vue.js)

フロントエンドにはVue.jsを使いました。
axiosを使って/apiにフォームの情報を投げています。
検索した範囲ではVue.jsではフォーム情報をこうやって取るのが正しい雰囲気があったのですが、本当にこれが正解なのかはよく分かりません。
結果が返ってきたら結果表示欄にリストで表示します。
FlaskのレンダーとdelimiterがかぶっているのでVue.js側のdelimiterを変更する必要がありました。

<html>
<head>
<script src="https://cdn.jsdelivr.net/npm/vue/dist/vue.js"></script>
<script src="https://unpkg.com/axios/dist/axios.min.js"></script>
</head>
<body>
<div id="app">
<form v-on:submit.prevent="submit">
<textarea v-model="formData.text"></textarea>
<input type="submit" value="submit">
</form>
<div v-if="seen">
<h2>Results</h2>
<ul>
<li v-for="result in results"> [[ result ]] </li>
</ul>
</div>
</div>
</body>
<script>
const vm = new Vue({
  el: '#app',
  delimiters: ['[[', ']]'], // Flaskのdelimiterとの重複を回避
  data: {
    seen: false,
    results: [],
    formData: {
      text: ''
    }
  },
  methods: {
    submit: function() {
      axios.post('/api', this.formData)
        .then(response => {
          this.results = response.data.results;
          this.seen = true;
        })
      .catch(error => {
        console.log(error);
      })
    }
  }
});
</script>
</html>

動作例

submitする前

f:id:kivantium:20180723173258p:plain

submitした後

f:id:kivantium:20180723173339p:plain


もっと良い書き方があったら教えてください

C/C++でのメモリリーク検出方法 〜AddressSanitizer, Valgrind, mtrace〜

C/C++でプログラムを書いているときに遭遇する厄介なバグの一つがメモリリークです。
今回はメモリリークを検出するのに使えるツールの使い方について書きます。

AddressSanitizer

コンパイルオプションをつけるだけで使えて出力も見やすいのでおすすめです。

AddressSanitizerはGCC 4.8以降かLLVM 3.1以降で使うことができます。
コンパイル時にオプションをつけるだけでメモリリークを検出してくれます。(若干実行時間が長くなります)

以下のメモリリークのあるプログラム leak.cpp を例に使い方を説明します。

int main() {
  int *a = new int[10];
}

newで作った動的配列をdeleteしていないのでメモリリークになります。

g++ -fsanitize=address -fno-omit-frame-pointer -g leak.cpp
./a.out

のようにオプションをつけてコンパイルして実行すると以下のような出力を得ます。(clangの場合もオプションは同じです)

=================================================================
==6027==ERROR: LeakSanitizer: detected memory leaks

Direct leak of 40 byte(s) in 1 object(s) allocated from:
    #0 0x7f606cd7a6b2 in operator new[](unsigned long) (/usr/lib/x86_64-linux-gnu/libasan.so.2+0x996b2)
    #1 0x4006f7 in main /home/kivantium/leak.cpp:2
    #2 0x7f606c93782f in __libc_start_main (/lib/x86_64-linux-gnu/libc.so.6+0x2082f)

SUMMARY: AddressSanitizer: 40 byte(s) leaked in 1 allocation(s).

leak.cppの2行目のメモリリークの存在が分かりました。

詳細はAddressSanitizer · google/sanitizers Wiki · GitHubが詳しいです。

Valgrind

仮想機械の上でプログラムを実行することでメモリリークの検出などを行うツールです。実行ファイルだけあれば使える点が便利ですが、実行はけっこう遅くなります。

インストールは基本的には

sudo apt install valgrind

でできるはずですが、僕の環境ではlinux-tools-4.13.0-45-genericのインストールも必要でした。

int main() {
  int *a = new int[10];
}

をValgrindでチェックするには、

g++ -g leak.cpp
valgrind --leak-check=full ./a.out

を実行します。結果は、

==6911== Memcheck, a memory error detector
==6911== Copyright (C) 2002-2015, and GNU GPL'd, by Julian Seward et al.
==6911== Using Valgrind-3.11.0 and LibVEX; rerun with -h for copyright info
==6911== Command: ./a.out
==6911== 
==6911== 
==6911== HEAP SUMMARY:
==6911==     in use at exit: 72,744 bytes in 2 blocks
==6911==   total heap usage: 2 allocs, 0 frees, 72,744 bytes allocated
==6911== 
==6911== 40 bytes in 1 blocks are definitely lost in loss record 1 of 2
==6911==    at 0x4C2E80F: operator new[](unsigned long) (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==6911==    by 0x400607: main (leak.cpp:2)
==6911== 
==6911== LEAK SUMMARY:
==6911==    definitely lost: 40 bytes in 1 blocks
==6911==    indirectly lost: 0 bytes in 0 blocks
==6911==      possibly lost: 0 bytes in 0 blocks
==6911==    still reachable: 72,704 bytes in 1 blocks
==6911==         suppressed: 0 bytes in 0 blocks
==6911== Reachable blocks (those to which a pointer was found) are not shown.
==6911== To see them, rerun with: --leak-check=full --show-leak-kinds=all
==6911== 
==6911== For counts of detected and suppressed errors, rerun with: -v
==6911== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 0 from 0)

のようになり、2行目にメモリリークがあることが分かります。

他にも機能があるようなので機会があったら追記します。

mtrace

メモリリークの位置を絞って探すときに使えます。glibcに入っているので使える環境が多いですが出力が見にくいです。

mtraceを使うにはmcheck.hをインクルードして、調べる範囲にmtrace()muntrace()を書く必要があります。

例としては

#include <stdlib.h>
#include <mcheck.h>
 
int main() {
  mtrace();
  int *a = (int*) malloc(sizeof(int)*10);
  muntrace();
}

のようになります。

実行するには

gcc -g leak.c
export MALLOC_TRACE=./mtrace.log
./a.out
mtrace a.out mtrace.log

のようにログ出力先を環境変数に指定したあとにプログラムを実行し、バイナリとログファイルを引数に与えてmtraceを実行します。

出力は

Memory not freed:
-----------------
           Address     Size     Caller
0x0000000001456450     0x28  at /home/kivantium/leak.cpp:6

のようになります。

C++

#include <mcheck.h>
 
int main() {
  mtrace();
  int *a = new int[10];
  muntrace();
}

のように書いた場合、出力は

Memory not freed:
-----------------
           Address     Size     Caller
0x00000000008b5060     0x28  at 0x7fe1589bbe78

のようになってしまい、どの行が問題だったのか分かりません。addr2lineを使って

addr2line -e a.out 0x7fe1589bbe78

とすればどの行か分かる場合もあるようですが、僕の環境では分かりませんでした。
C++ではかなり使いにくいようです……

制約付き最適化問題をD-Waveと同じ方法で解く

まだ古典コンピュータで消耗してるの?

イジングモデルとQUBO

D-Waveが提供しているのは量子アニーリングを使った量子コンピュータです。
量子アニーリング方式では次の形の最適化問題を高速に解くことができます。
 {\displaystyle \text{min.}\ \sum_{i \lt j} J_{ij} s_{i} s_{j} - \sum_{i=1}^{N} h_i s_i }
ここで J_{ij}, h_{i} \in \mathbb{R} s_i \in \{1, -1\}です。

1と-1よりは0と1の方が扱いやすいので、 x_i = (s_i+1)/2と変数変換して得られる以下のQUBO (Quadratic unconstrained binary optimization) という形に変形してから使うことが多いようです。
 {\text{min. } \displaystyle \sum_{i=1}^{N} c_i x_i + \sum_{i=1}^{N} \sum_{j=1}^{i} q_{ij} x_i x_j }

制約付き最適化問題とQUBO

QUBOでは変数間の制約がありませんでしたが、実際には  x_1+x_2=1 のような制約が必要になることがあります。
その場合、最適化したい関数 f(x)と十分大きい \lambda > 0を使って
 {\displaystyle \text{min.} f(x) + \lambda (x_1+x_2-1)^2 }
を解けば、(実行可能解が存在するなら)最小値を取るときには自動的に制約が満たされることになるので、QUBOの形で制約付き最適化問題が解けることになります。
なお、QUBOには x_i^2 の項がありませんが、 x_i は0または1なので  x_i^2=x_i が成り立ち、2乗の項が出てきてもQUBOの形に変形することができます。

qbsolv

qbsolv はD-Waveが開発しているQUBOのソルバです。
設定すればD-Waveの実機上でも動かせるようですが、デフォルトではローカル環境でのシミュレーションを使って解く設定になっているようです。
ここではコマンドライン環境を使います。

コンパイル

git clone https://github.com/dwavesystems/qbsolv.git
cd qbsolv
mkdir build; cd build
cmake ..
make

qbsolvは.quboファイルを読み込んでその解を出力します。
.quboファイルの仕様は以下の通りです。(READMEより引用)

qbsolv "qubo" input file format

   A .qubo file contains data which describes an unconstrained
quadratic binary optimization problem.  It is an ASCII file comprised
of four types of lines:

1) Comments - defined by a "c" in column 1.  They may appear
    anywhere in the file, and are otherwise ignored.

2) One program line, which starts with p in the first column.
    The program line must be the first non-comment line in the file.
    The program line has six required fields (separated by space(s)),
    as in this example:

  p   qubo  topology   maxNodes   nNodes   nCouplers

    where:
  p         the problem line sentinel
  qubo      identifies the file type
  topology  a string which identifies the topology of the problem
            and the specific problem type.  For an unconstrained problem,
            target will be "0" or "unconstrained."   Possible, for future
            implementations, valid strings might include "chimera128"
            or "chimera512" (among others).
  maxNodes   number of nodes in the topology.
  nNodes     number of nodes in the problem (nNodes <= maxNodes).
            Each node has a unique number and must take a value in the
            the range {0 - (maxNodes-1)}.  A duplicate node number is an
            error.  The node numbers need not be in order, and they need
            not be contiguous.
  nCouplers  number of couplers in the problem.  Each coupler is a
            unique connection between two different nodes.  The maximum
            number of couplers is (nNodes)^2.  A duplicate coupler is
            an error.

3) nNodes clauses.  Each clause is made up of three numbers.  The
            numbers are separated by one or more blanks.  The first two
            numbers must be integers and are the number for this node
            (repeated).  The node number must be in {0 , (maxNodes-1)}.
            The third value is the weight associated with the node, may be
            an integer or float, and can take on any positive or negative
            value, or zero.

4) nCouplers clauses.  Each clause is made up of three numbers.  The
            numbers are separated by one or more blanks.  The first two
            numbers must be different integers and are the node numbers
            for this coupler.  The two values (i and j) must have (i < j).
            Each number must be one of the nNodes valid node numbers (and
            thus in {0, (maxNodes-1)}).  The third value is the strength
            associated with the coupler, may be an integer or float, and can
            take on any positive or negative value, but not zero.  Every node
            must connect with at least one other node (thus must have at least
            one coupler connected to it).

Here is a simple QUBO file example for an unconstrained QUBO with 4
nodes and 6 couplers.  This example is provided to illustrate the
elements of a QUBO benchmark file, not to represent a real problem.

        | <--- column 1
        c
        c  This is a sample .qubo file
        c  with 4 nodes and 6 couplers
        c
        p  qubo  0  4  4  6
        c ------------------
        0  0   3.4
        1  1   4.5
        2  2   2.1
        3  3   -2.4
        c ------------------
        0  1   2.2
        0  2   3.4
        1  2   4.5
        0  3   -2
        1  3   4.5678
        2  3   -3.22

最短路問題をQUBOで解く

制約付き最適化問題の例として、最短路問題をQUBOに変換して解きます。
良い子のみんなはダイクストラなどを使いましょう。

f:id:kivantium:20180617163419p:plain
今回考えるグラフ

最短路問題を最適化問題として解く時は、辺  e_i を通るときは  x_i=1、通らない時は  x_i=0 として、重みの和  \sum w_i x_i が最小になるようにする問題と見ます。( w_i は辺  e_i の重み)
これだけだと何も辺を選ばないのが最適になってしまうので、辺の選び方が経路になるように各頂点に制約を加えます。

  •  v_0 からは1辺が出ていくので  x_0 + x_2 = 1
  •  v_1 に入って来る辺の数と出て行く辺の数は等しいので  x_0 - x_1 = 0
  •  v_2 には1辺が入ってくるので  x_1 + x_2 = 1

 (w_0, w_1, w_2) = (1, 1, 3) とすると、この最短路問題を解くためのQUBOは、
 {\displaystyle \text{min. } x_0 + x_1 + 3x_2 + \lambda \{(x_0+x_2-1)^2 + (x_0-x_1)^2+(x_1+x_2-1)^2\} }
となります。
 \lambda=100として、 x_i^2=x_iに気をつけながら変形すると
 {\displaystyle \text{min. } x_0 + x_1 +  (-197)x_3 + (-200) x_0x_1+200x_0x_2+200x_1x_2 }
となります。(間違ってたらごめんなさい)

これを.quboファイルにすると

p  qubo  0  3  3  3
c ------------------
0  0   1
1  1   1
2  2   -197
c ------------------
0  1   -200
0  2   200
1  2   200

となります。

./qbsolv -i test.qubo

のように実行すると

3 bits,  find Min, SubMatrix= 47, -a o, timeout=2592000.0 sec
110
-198.00000 Energy of solution
0 Number of Partitioned calls, 1 output sample 
 0.02363 seconds of classic cpu time

となり、 (x_0, x_1, x_2) = (1, 1, 0) という正しい答えが得られました。

 w_0=3にすると答えが

3 bits,  find Min, SubMatrix= 47, -a o, timeout=2592000.0 sec
001
-197.00000 Energy of solution
0 Number of Partitioned calls, 1 output sample 
 0.02023 seconds of classic cpu time

に変わるので正しそうな挙動をしています。

Open Babelのビルドとインストール

Open Babel は化学で使われる分子フォーマット間の変換によく使われるソフトです。

開発されたばかりの機能を使うためには自分でビルドしてインストールする必要があるのですが、Pythonバインディングコンパイルするところにハマったのでメモしておきます。 基本的にはInstall Open Babelを見ればいいです。 環境はUbuntu 16.04を使いました。

Anacondaでリリース済みバージョンを使う

Anacondaをインストールしていれば

conda -c openbabel openbabel

とすればよいです。

自分でコンパイルする

Anacondaを使っていない場合

必要なライブラリをインストールします。

sudo apt install cmake python-dev libxml2-dev zlib1g-dev libeigen2-dev libcairo2-dev swig2.0

ビルドとインストールは以下のようにします。

git clone https://github.com/openbabel/openbabel.git
cd openbabel
mkdir build
cd build
cmake .. -DPYTHON_BINDINGS=ON 
make
sudo make install
echo 'export PYTHONPATH=/usr/local/lib:$PYTHONPATH' >> ~/.bashrc
source ~/.bashrc

pythonを起動して import openbabel を実行してエラーが起きなければ成功です。

(2019/11/22 追記) Open Babel 3.0のリリース時点ではSWIGがないとPythonバインディングコンパイルできないようです。 Python bindings don't build in the 3.0 release tarball due to incorrect path · Issue #2081 · openbabel/openbabel · GitHub

sudo apt install build-essentials cmake python3-dev libeigen3-dev swig wget
wget https://github.com/openbabel/openbabel/archive/openbabel-3-0-0.tar.gz
tar xf openbabel-3-0-0.tar.gz
cd openbabel-openbabel-3-0-0/
mkdir build; cd $_
cmake ..  -DPYTHON_BINDINGS=ON -DRUN_SWIG=ON
make
sudo make install
echo 'export PYTHONPATH=/usr/local/lib:$PYTHONPATH' >> ~/.bashrc
source ~/.bashrc

コンパイルできました。

Anacondaを使っている場合

Pybel import crashes Anaconda Python3.6 Ubuntu · Issue #1805 · openbabel/openbabel · GitHub に書いてあります。

conda install cmake lxml swig

cmake .. -DPYTHON_BINDINGS=ON -DRUN_SWIG=ON -DCMAKE_INSTALL_PREFIX=~/anaconda3 -DPYTHON_INCLUDE_DIR=~/anaconda3/include/python3.6m/ -DCMAKE_LIBRARY_PATH=~/anaconda3/lib -DSWIG_DIR=~/anaconda3/share/swig/3.0.12/ -DSWIG_EXECUTABLE=~/anaconda3/bin/swig -DPYTHON_LIBRARY=~/anaconda3/lib/libpython3.6m.so

make && make install

C++で使う

インストールすると obabelなどの通常のコマンドラインツールが使えるようになります。

より複雑なことをするにはC++Pythonなどでコードを書くことになります。 以下に標準入力からSMILESを読んで、3次元座標を計算した後、標準出力にSDFを吐くコード例を示します。 エラー処理は省略しています。

#include <iostream>

#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/builder.h>

int main(int argc,char **argv) {
    OpenBabel::OBConversion conv(&std::cin, &std::cout);
    OpenBabel::OBMol mol;
    OpenBabel::OBBuilder builder;

    conv.SetInAndOutFormats("SMI", "SDF");
    conv.Read(&mol);
    builder.Build(mol);
    conv.Write(&mol);
}

コンパイルは以下のように行います。ディレクトリは環境によって異なるかもしれません。

g++ -I /usr/local/include/openbabel-2.0/  -L /usr/local/lib/openbabel/2.4.90/ hoge.cpp  -lopenbabel

Makefileの例を示します。ディレクトリ名は適宜変更してください。

CC = g++
CXXFLAGS = -I /usr/local/include/openbabel-2.0/
LDFLAGS = -L /usr/local/lib/openbabel/2.4.90/
LDLIBS = -lopenbabel

all: hoge fuga

SMILESをstringで指定する例

#include <iostream>
#include <string>
#include <sstream>

#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/builder.h>

int main(int argc, char** argv){
  std::string smiles = "FCl(=O)=O";
  std::stringstream ss(smiles);
  OpenBabel::OBConversion conv(&ss, &std::cout);
  conv.SetInAndOutFormats("smi", "sdf");
  OpenBabel::OBMol mol;
  conv.Read(&mol);

  OpenBabel::OBBuilder builder;
  builder.Build(mol);

  conv.Write(&mol);
}

広告コーナー