はじめてのCUDAプログラミングで分子動力学計算

追記(2009/12/12): ここで示してるコードはかなりひどいので少し手直しして、CUDAで作成した分子動力学計算プログラムを書き直してみたという記事にまとめた。こちらも参照して欲しい。

最近、自分の周りではCUDAを使ったGPGPUが流行っているので、どんなものかこの週末にいじってみた。

一日目はNVIDIA GeForce 8800 GTSを搭載したWindows XPマシンにCUDAをインストールをしてリファレンスやサンプルプログラムをパラパラと読んでみた。二日目になんちゃって分子動力学(MD)計算プログラムを書いた。擬似原子を扱っており、全原子のvdW力とクーロン力のみを考慮して分子内結合については考慮しない、液滴中のMD計算だ。取り敢えず、半径20Åの液滴中に1229原子を入れて1000ステップ計算させてみた。C++で書いたプログラムでは4分30秒かかった計算がCUDAでは32秒で計算できた。約9倍ほどの速度が出ているようだ。

まあ、C++のコードも併せて一日で書いたので最適化や高速化については考慮していない。後でスレッドでの分割がやりやすいかもと思ってiとjをそのまま全部計算しているし。気が向いたらちゃんとしたプログラムを書いてみるかも。

CUDAアーキテクチャについてはまだまだ理解が足りない。効率のよいメモリの扱い方を覚えないとなぁ。あと、sharedメモリをうまく扱えてないっぽい。さすがに一日二日じゃキツイか。

初心者が一日で書いたコードなので参考になるかわからないが、以下にソースコードとサンプルの入力データを示しておく。入力データについては、上記で用いた1229原子ではブログ記事するには大きすぎるので154原子の小さなデータを用意した。また、記事の表示量が多くなりすぎるのでC++のソースは割愛させてもらった。Windows上でのコンパイルは下記の通り。

nvcc -lcutil32 -lkernel32 --host-compilation c++ -Xcompiler /EHsc,/W3,/nologo,/O2,/Zi,/MT simple_md_gpu.cu -o simple_md_gpu.exe

入力ファイルをmd.datとして、100ステップ計算させる場合は、以下のように実行する。

simple_md_gpu md.dat 100

最初に初期構造が表示され、次いでステップごとのエネルギー、最後に最終構造が出力される。

simple_md_gpu.cu

//////////////////////////////////////////////////////////////// // Simple Molecular Dynamics using GPU by nox, 2009.03.22. // // Perform molecular dynamics for pseudo-atoms without internal // interactions, such as bonds, angles, dihedrals. Only vdW and // Coulomb interactions. #include <iostream> #include <iomanip> #include <fstream> #include <sstream> #include <cmath> #include <cutil_inline.h> using namespace std; const int BLOCK = 256; const int WIDTH = 2048; __device__ float Distance(float* crd, int i, int j) { float dx = crd[i * 3] - crd[j * 3]; float dy = crd[i * 3 + 1] - crd[j * 3 + 1]; float dz = crd[i * 3 + 2] - crd[j * 3 + 2]; return sqrtf(dx * dx + dy * dy + dz * dz); } __device__ float Energy(float* ene, int num_atoms) { float sum = 0.0f; for (int i = 0; i < num_atoms; i++) sum += ene[i]; return sum / 2.0f; } __global__ void Center(float* crd,int num_atoms) { __shared__ float d[3]; for (int i = 0; i < 3; i++) d[i] = 0.0f; for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) d[j] += crd[i * 3 + j]; for (int i = 0; i < 3; i++) d[i] /= num_atoms; for (int i = 0; i < num_atoms; i++) for (int j = 0; j < 3; j++) crd[i * 3 + j] -= d[j]; } // Calculate energies and forces // Energy = vdW + Coulomb // vdW // U = eps * [(Rij / r[i])^12 - 2 * (Rij / r[i])^6] // F = -12 * eps / Rij * [(Rij / r[i])^13 - (Rij / r[i])^7] * r_xyz / r[i] // Coulomb // U = SUM_i>j qiqj / r[i] // F = SUM_j qiqj / r[i]^3 * r_xyz __global__ void Calc(float* crd, float* f, float* ene, float* chg, int num_atoms, float* tene) { int dtx = blockIdx.x * blockDim.x + threadIdx.x; const float eps = 0.2f; const float Rij = 2.5f; float r, r_xyz, df, r12, r6, q; for (int i = dtx; i < num_atoms; i += WIDTH) ene[i] = 0.0f; __syncthreads(); for (int i = dtx; i < num_atoms; i += WIDTH) { for (int j = 0; j < num_atoms; j++) { if (i == j) continue; r = Distance(crd, i, j); q = chg[i] * chg[j]; r12 = powf(Rij / r, 12); r6 = powf(Rij / r, 6); ene[i] += eps * (r12 - 2.0f * r6) + q / r; for (int k = 0; k < 3; k++) { r_xyz = crd[i * 3 + k] - crd[j * 3 + k]; df = (12 * eps / Rij) * (r12 * Rij / r - r6 * Rij / r) * r_xyz / r; df += q / powf(r, 3) * r_xyz; f[i * 3 + k] += df; } } } __syncthreads(); *tene = Energy(ene, num_atoms); } __global__ void Move(float* crd, float* f, int num_atoms, float cap_range) { int dtx = blockIdx.x * blockDim.x + threadIdx.x; float r, dr; for (int i = dtx; i < num_atoms; i += WIDTH) { for (int j = 0; j < 3; j++) { crd[i * 3 + j] += f[i * 3 + j] * 0.01f; f[i * 3 + j] = 0.0f; } r = 0.0f; for (int j = 0; j < 3; j++) r += crd[i * 3 + j] * crd[i * 3 + j]; dr = r - cap_range * cap_range; if (dr > 0.0f) for (int j = 0; j < 3; j++) f[i * 3 + j] -= crd[i * 3 + j] / cap_range * dr * 0.01f; } } class SimpleMD { private: int num_steps; int num_atoms; float cap_range; float *h_crd, *h_f, *h_ene, *h_chg; float *d_crd, *d_f, *d_ene, *d_chg; float tene; public: SimpleMD(int n, char*); ~SimpleMD(); void ReadCrd(char*); void PrintCrd(); void Run(); }; SimpleMD::SimpleMD(int n, char* fname) : num_steps(n) { ReadCrd(fname); h_f = new float[num_atoms * 3]; h_ene = new float[num_atoms]; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) h_f[i * 3 + j] = 0.0f; h_ene[i] = 0.0f; } cudaMalloc((void**)&d_crd, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_f, sizeof(float) * num_atoms * 3); cudaMalloc((void**)&d_ene, sizeof(float) * num_atoms); cudaMalloc((void**)&d_chg, sizeof(float) * num_atoms); } SimpleMD::~SimpleMD() { cudaFree(d_chg); cudaFree(d_ene); cudaFree(d_f); cudaFree(d_crd); delete[] h_ene; delete[] h_f; delete[] h_chg; delete[] h_crd; } void SimpleMD::ReadCrd(char* fname) { fstream fs(fname, ios_base::in); string line; stringstream ss; if (!fs.is_open()) { cerr << "File open error: " << fname << endl; exit(1); } getline(fs, line); cout << line << endl; getline(fs, line); ss.str(line); ss >> num_atoms; ss.clear(); getline(fs, line); ss.str(line); ss >> cap_range; ss.clear(); h_crd = new float[num_atoms * 3]; h_chg = new float[num_atoms]; for (int i = 0; i < num_atoms; i++) { getline(fs, line); ss.str(line); ss >> h_crd[i * 3] >> h_crd[i * 3 + 1] >> h_crd[i * 3 + 2] >> h_chg[i]; ss.clear(); ss.str(""); } fs.close(); } void SimpleMD::PrintCrd() { cout << endl; cout << num_atoms << endl; cout << cap_range << endl; for (int i = 0; i < num_atoms; i++) { for (int j = 0; j < 3; j++) cout << " " << fixed << setw(10) << setprecision(6) << h_crd[i * 3 + j]; cout << " " << fixed << setw(8) << setprecision(4) << h_chg[i]; cout << endl; } } void SimpleMD::Run() { float *d_tene; cudaMalloc((void**)&d_tene, sizeof(float)); cudaMemcpy(d_crd, h_crd, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_f, h_f, sizeof(float) * num_atoms * 3, cudaMemcpyHostToDevice); cudaMemcpy(d_ene, h_ene, sizeof(float) * num_atoms, cudaMemcpyHostToDevice); cudaMemcpy(d_chg, h_chg, sizeof(float) * num_atoms, cudaMemcpyHostToDevice); dim3 grid(WIDTH / BLOCK, 1, 1); dim3 threads(BLOCK, 1, 1); Center<<<1, 1>>>(d_crd, num_atoms); for (int i = 0; i < num_steps; i++) { Calc<<<grid, threads>>>(d_crd, d_f, d_ene, d_chg, num_atoms, d_tene); Move<<<grid, threads>>>(d_crd, d_f, num_atoms, cap_range); Center<<<1, 1>>>(d_crd, num_atoms); cudaMemcpy(&tene, d_tene, sizeof(float), cudaMemcpyDeviceToHost); cout << "Energy (" << setw(7) << i + 1 << "): "; cout << fixed << setw(15) << setprecision(5) << tene << endl; } cudaMemcpy(h_crd, d_crd, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_f, d_f, sizeof(float) * num_atoms * 3, cudaMemcpyDeviceToHost); cudaMemcpy(h_ene, d_ene, sizeof(float) * num_atoms, cudaMemcpyDeviceToHost); cudaFree(d_tene); } int main(int argc, char** argv) { if (argc != 3) { cerr << "Usage: " << argv[0] << " input_file number_of_steps" << endl; cerr << "Input: line 1 : title" << endl; cerr << " line 2 : number of atoms" << endl; cerr << " line 3 : radius of droplet" << endl; cerr << " line 4-: x-crd y-crd z-crd charge" << endl; exit(1); } stringstream ss; int nstep; cutilDeviceInit(1, argv); ss.str(argv[2]); ss >> nstep; SimpleMD md(nstep, argv[1]); md.PrintCrd(); md.Run(); md.PrintCrd(); return 0; }

サンプル入力ファイル: md.dat

simple_md_gpu sample data 154 10.0 -7.00000 -7.00000 -1.00000 -1.000 -7.00000 -4.00000 -4.00000 1.000 -7.00000 -4.00000 -1.00000 -1.000 -7.00000 -4.00000 2.00000 1.000 -7.00000 -4.00000 5.00000 -1.000 -7.00000 -1.00000 -7.00000 1.000 -7.00000 -1.00000 -4.00000 -1.000 -7.00000 -1.00000 -1.00000 1.000 -7.00000 -1.00000 2.00000 -1.000 -7.00000 -1.00000 5.00000 1.000 -7.00000 2.00000 -4.00000 -1.000 -7.00000 2.00000 -1.00000 1.000 -7.00000 2.00000 2.00000 -1.000 -7.00000 2.00000 5.00000 1.000 -7.00000 5.00000 -4.00000 -1.000 -7.00000 5.00000 -1.00000 1.000 -7.00000 5.00000 2.00000 -1.000 -7.00000 5.00000 5.00000 1.000 -4.00000 -7.00000 -4.00000 -1.000 -4.00000 -7.00000 -1.00000 1.000 -4.00000 -7.00000 2.00000 -1.000 -4.00000 -7.00000 5.00000 1.000 -4.00000 -4.00000 -7.00000 -1.000 -4.00000 -4.00000 -4.00000 1.000 -4.00000 -4.00000 -1.00000 -1.000 -4.00000 -4.00000 2.00000 1.000 -4.00000 -4.00000 5.00000 -1.000 -4.00000 -4.00000 8.00000 1.000 -4.00000 -1.00000 -7.00000 -1.000 -4.00000 -1.00000 -4.00000 1.000 -4.00000 -1.00000 -1.00000 -1.000 -4.00000 -1.00000 2.00000 1.000 -4.00000 -1.00000 5.00000 -1.000 -4.00000 -1.00000 8.00000 1.000 -4.00000 2.00000 -7.00000 -1.000 -4.00000 2.00000 -4.00000 1.000 -4.00000 2.00000 -1.00000 -1.000 -4.00000 2.00000 2.00000 1.000 -4.00000 2.00000 5.00000 -1.000 -4.00000 2.00000 8.00000 1.000 -4.00000 5.00000 -7.00000 -1.000 -4.00000 5.00000 -4.00000 1.000 -4.00000 5.00000 -1.00000 -1.000 -4.00000 5.00000 2.00000 1.000 -4.00000 5.00000 5.00000 -1.000 -4.00000 8.00000 -4.00000 1.000 -4.00000 8.00000 -1.00000 -1.000 -4.00000 8.00000 2.00000 1.000 -1.00000 -7.00000 -7.00000 -1.000 -1.00000 -7.00000 -4.00000 1.000 -1.00000 -7.00000 -1.00000 -1.000 -1.00000 -7.00000 2.00000 1.000 -1.00000 -7.00000 5.00000 -1.000 -1.00000 -4.00000 -7.00000 1.000 -1.00000 -4.00000 -4.00000 -1.000 -1.00000 -4.00000 -1.00000 1.000 -1.00000 -4.00000 2.00000 -1.000 -1.00000 -4.00000 5.00000 1.000 -1.00000 -4.00000 8.00000 -1.000 -1.00000 -1.00000 -7.00000 1.000 -1.00000 -1.00000 -4.00000 -1.000 -1.00000 -1.00000 -1.00000 1.000 -1.00000 -1.00000 2.00000 -1.000 -1.00000 -1.00000 5.00000 1.000 -1.00000 -1.00000 8.00000 -1.000 -1.00000 2.00000 -7.00000 1.000 -1.00000 2.00000 -4.00000 -1.000 -1.00000 2.00000 -1.00000 1.000 -1.00000 2.00000 2.00000 -1.000 -1.00000 2.00000 5.00000 1.000 -1.00000 2.00000 8.00000 -1.000 -1.00000 5.00000 -7.00000 1.000 -1.00000 5.00000 -4.00000 -1.000 -1.00000 5.00000 -1.00000 1.000 -1.00000 5.00000 2.00000 -1.000 -1.00000 5.00000 5.00000 1.000 -1.00000 5.00000 8.00000 -1.000 -1.00000 8.00000 -4.00000 1.000 -1.00000 8.00000 -1.00000 -1.000 -1.00000 8.00000 2.00000 1.000 -1.00000 8.00000 5.00000 -1.000 2.00000 -7.00000 -4.00000 1.000 2.00000 -7.00000 -1.00000 -1.000 2.00000 -7.00000 2.00000 1.000 2.00000 -7.00000 5.00000 -1.000 2.00000 -4.00000 -7.00000 1.000 2.00000 -4.00000 -4.00000 -1.000 2.00000 -4.00000 -1.00000 1.000 2.00000 -4.00000 2.00000 -1.000 2.00000 -4.00000 5.00000 1.000 2.00000 -4.00000 8.00000 -1.000 2.00000 -1.00000 -7.00000 1.000 2.00000 -1.00000 -4.00000 -1.000 2.00000 -1.00000 -1.00000 1.000 2.00000 -1.00000 2.00000 -1.000 2.00000 -1.00000 5.00000 1.000 2.00000 -1.00000 8.00000 -1.000 2.00000 2.00000 -7.00000 1.000 2.00000 2.00000 -4.00000 -1.000 2.00000 2.00000 -1.00000 1.000 2.00000 2.00000 2.00000 -1.000 2.00000 2.00000 5.00000 1.000 2.00000 2.00000 8.00000 -1.000 2.00000 5.00000 -7.00000 1.000 2.00000 5.00000 -4.00000 -1.000 2.00000 5.00000 -1.00000 1.000 2.00000 5.00000 2.00000 -1.000 2.00000 5.00000 5.00000 1.000 2.00000 5.00000 8.00000 -1.000 2.00000 8.00000 -4.00000 1.000 2.00000 8.00000 -1.00000 -1.000 2.00000 8.00000 2.00000 1.000 2.00000 8.00000 5.00000 -1.000 5.00000 -7.00000 -4.00000 1.000 5.00000 -7.00000 -1.00000 -1.000 5.00000 -7.00000 2.00000 1.000 5.00000 -7.00000 5.00000 -1.000 5.00000 -4.00000 -7.00000 1.000 5.00000 -4.00000 -4.00000 -1.000 5.00000 -4.00000 -1.00000 1.000 5.00000 -4.00000 2.00000 -1.000 5.00000 -4.00000 5.00000 1.000 5.00000 -1.00000 -7.00000 -1.000 5.00000 -1.00000 -4.00000 1.000 5.00000 -1.00000 -1.00000 -1.000 5.00000 -1.00000 2.00000 1.000 5.00000 -1.00000 5.00000 -1.000 5.00000 -1.00000 8.00000 1.000 5.00000 2.00000 -7.00000 -1.000 5.00000 2.00000 -4.00000 1.000 5.00000 2.00000 -1.00000 -1.000 5.00000 2.00000 2.00000 1.000 5.00000 2.00000 5.00000 -1.000 5.00000 2.00000 8.00000 1.000 5.00000 5.00000 -7.00000 -1.000 5.00000 5.00000 -4.00000 1.000 5.00000 5.00000 -1.00000 -1.000 5.00000 5.00000 2.00000 1.000 5.00000 5.00000 5.00000 -1.000 5.00000 8.00000 -1.00000 1.000 5.00000 8.00000 2.00000 -1.000 8.00000 -4.00000 -4.00000 1.000 8.00000 -4.00000 -1.00000 -1.000 8.00000 -4.00000 2.00000 1.000 8.00000 -1.00000 -4.00000 -1.000 8.00000 -1.00000 -1.00000 1.000 8.00000 -1.00000 2.00000 -1.000 8.00000 -1.00000 5.00000 1.000 8.00000 2.00000 -4.00000 -1.000 8.00000 2.00000 -1.00000 1.000 8.00000 2.00000 2.00000 -1.000 8.00000 2.00000 5.00000 1.000 8.00000 5.00000 -1.00000 -1.000 8.00000 5.00000 2.00000 1.000

追記:

特に記述しなかったが、出力エネルギーはポテンシャルエネルギーなので、念のため。また、ここでは時間積分について考慮していない。NVEアンサンブル(定エネルギー)であれば、時間積分は速度ベルレ法がもっとも実装しやすいように思う。dt/2で速度を求め、そこからdtの座標を計算し、その座標からさらにdt/2で速度を求める。