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言語入門

2019年2月

2月もあっという間に終わってしまいました。

今月のトピック

今月読んだ本

今月は論文を書くのに忙しくてあまり本を読めませんでした(大嘘)

天冥の標X Part3

全ての力を尽くして天冥の標シリーズをオススメする - 基本読書を読んだのをきっかけに追い始めたシリーズの完結篇。
2013年1月の記事なので読み始めてから6年間経ったことになる。

2013年というと高校2年生で、当時の自分は医学部に行こうと思っていたはずのだが、人類の軌跡を宇宙規模の争いと絡めながら語っていくスケールの大きなお話を読んでしまった結果、自分がやりたいのは人類の科学を進めて歴史を動かすことであってどうせ死ぬ人間をほんの少し延命させることではないななどとうっかり思ってしまい(このシリーズの主人公は医者なのに!)、医学部志望をやめてよく分からない学科を受験してしまったので、私の人生は天冥の標に狂わされたと言っていいだろう。
そんな思い入れの深いシリーズが予定通り10巻で終わってしまうことになり完結篇を読むのが非常に怖かったのだが、楽しみにしていてよかったと思える結末に辿りついてくれていた。

1巻は導入としてだけではなく植民星での革命劇として純粋に面白かったし、2巻の感染症との戦いは素晴らしく、3巻の世界観や酸素要らずの価値観には人生を破壊され、4巻はおいといて、5巻のスケールの大きさには度肝を抜かれ、6巻での戦いに圧倒され、7巻での絶望に心を打たれ、8巻での伏線回収に感心し、9,10巻での物語の畳み方には舌を巻いた素晴らしいシリーズであった。最終巻のあとがきに「読者の何割かに多少のエネルギーと二、三個の発見を提供できたかもしれない」とあったが、少なくとも私には大変なエネルギーを与えてくれる作品だった。

かぐや様は告らせたい 1〜13巻

アニメ放映に合わせてBookwalkerで80%オフキャンペーンをやっていたので全巻購入した。「天才たちの恋愛頭脳戦」をやっていた最初の頃は正直しょうもないことやってるなと思っていまいちだったのだが、5巻収録の「花火の音は聞こえない」以降、かぐや様が天才設定を捨てたベタ惚れポンコツ乙女になってからは毎回素晴らしかった。恋愛ものはあまり好きじゃないと思っていたが、単に面白い作品に出会ってないだけだったのではと心変わりした。

五等分の花嫁 1〜8巻

かぐや様を読んでラブコメもいいなと思っていたところに80%オフキャンペーンをやっていたのを見たので購入。性格悪いヒロインが多くてあまり好きになれなかったのだが、7巻あたりで皆がデレ出してからは面白く、8巻での二乃の一転攻勢は素晴らしくてここまで読んでよかったと思える出来だった。(ここまで書いてきて思ったが単に自分はデレデレヒロインが好きなだけなのでは……)

山田エルフ大先生の恋する純真ごはん

kivantium花嫁修業の元ネタである山田エルフ先生が料理を作る漫画。連載されると聞いてからはマガジンWALKERに契約して毎月電撃大王を読んでいたので全て読んだ話ではあったがお布施。

みらいのふうふですけど?

みらいのふうふですけど?(2) (百合姫コミックス)

みらいのふうふですけど?(2) (百合姫コミックス)

中学生百合。二人のコミカルなやりとりが素晴らしい作品だったが2巻で完結してしまった。作者が途上のIPコンテンツと言っているので打ち切りっぽい。悲しい。

お兄ちゃんはおしまい!

お兄ちゃんはおしまい!   (2) (IDコミックス)

お兄ちゃんはおしまい! (2) (IDコミックス)

妹に女の子にされてしまったニートのお兄ちゃんの話。1巻を買ったので2巻も。

GUNSLINGER GIRL

セールをやっていたので1巻だけ購入。全体に漂う暗い雰囲気や諦念が好みじゃなかったので2巻以降はまたの機会に。

量子化学

量子化学 上巻

量子化学 上巻

ラボの人に勧められたので読んだ。仮定と近似がはっきり書かれていており、説明も省略が少なかったので分かりやすかった。

2019年1月

2018年12月分も含みます。

今月のトピック

行列の微分を求めるサイト

NIPS 2018で発表された成果らしい: Computing Higher Order Derivatives of Matrix and Tensor Expressions

Vimプラグイン ale

flaskなどのlinterをVim上で走らせることができる上に、保存したときにインデントを自動で修正することもできる。便利。

その他

https://twitter.com/tmaehara/status/937385784440582144



読んだ本

量子計算理論 量子コンピュータの原理

量子計算理論 量子コンピュータの原理

量子計算の理論に関する本。古典確率計算を量子計算のようなベクトル表示で表すことで対比する導入は面白かった。
物理的な実現方法やShorのアルゴリズムのような具体的な量子アルゴリズムは扱っておらず、「量子コンピュータを持っていると主張しているQ-Wave社が本当に量子計算を行っていると検証するにはどうすればよいか」のような理論について扱っていた。論文のダイジェストという感じですんなり読める感じではなかった。
以前読んだ量子コンピュータ入門(第2版)のほうが一冊目としては良いだろう。

ここまでの9巻で広げた風呂敷を回収していくのがすごく良い。

まほろまてぃっくの作者の新作。ママ(小学生)がかわいい。

三ツ星カラーズ6 (電撃コミックスNEXT)

三ツ星カラーズ6 (電撃コミックスNEXT)

ハミングガール (百合姫コミックス)

ハミングガール (百合姫コミックス)

2018年11月

今月のトピック

PRMLの無料公開

PRML以外で無料公開されている有名本・PDFを挙げておく

いろいろな法則とか

こういうのは名前を覚えておかないとドヤ顔で披露できないので記憶する必要がある。

Content-Disposition レスポンスヘッダー

ブラウザでPDFを開いたときに、ブラウザで開かれる場合とダウンロードされる場合があるが、この挙動はHTTPヘッダーで指定できる。
HTML5のdownload属性で指定される場合やJavaScriptでダウンロードされる場合もある)

その他



2018年10月

今月知った有益情報

Multi-service image search

複数サービスを横断して画像検索ができる

UCL Course on RL

http://www0.cs.ucl.ac.uk/staff/d.silver/web/Teaching.html
強化学習のオンライン講座。誰かが分かりやすいと言っていた。

The Architecture of Open Source Applications

著名なOSSプロジェクトの内部について簡単に説明していく本らしい。

古地図

いろんな古地図を見ることができる

Exploit-Exercises: Protostar (v2)

CTFのpwn問題の勉強におすすめの問題とのこと

電子サイコロの回路出典

4015を1つ使うだけでサイコロを作れる回路の出典は「初歩のラジオ1980年10月号に松本悟さんが発表された記事」とのこと

読んだ本

量子コンピュータ入門(第2版)

量子コンピュータ入門(第2版)

ゲート型量子コンピュータの解説本。今まで読んだ量子コンピュータ本の中で一番わかりやすかった。

青葉ちゃんと紅葉ちゃんのコンペ対決。人々が頑張るのはよい。

まちカドまぞく (4) (まんがタイムKRコミックス)

まちカドまぞく (4) (まんがタイムKRコミックス)

4コマオブザイヤーで3年連続1位になっている作品。不思議な魅力がある。4巻ではギャグ描写だと思っていたものが展開上重要なアイテムになる展開に感心した。

3巻までは能力を持たないナナが能力者たちと対決する話として非常によくできており、4巻では少し方向性を変えつつナナの過去や能力者との友情が描かれていた。面白かったが、展開や新キャラからそこはかとない車輪の国っぽさを感じてしまった。

零號琴

零號琴

飛 浩隆の16年ぶりの長編。巨大楽器の竣工を記念した仮面舞踏会という面白くなさそうなあらすじだったが、高い描写力や緻密な世界構築に圧倒された。さすがにこれはふざけすぎではとか展開に納得感がなかったなぁなどと思わないでもない。

RISC-V原典  オープンアーキテクチャのススメ

RISC-V原典 オープンアーキテクチャのススメ

RISC-Vの命令セットについての本。ちょっとだけお気持ちが書いてあるが基本的にデータシートと同じ内容が書いてあるので買うほどではないかもしれない。

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 <イメージ名>

とすればよいです。

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

特定商取引法に定められた事項は請求により遅滞なく提供する