GPUによる量子シミュレーションの加速現在量子系のシミュレーションに最もよく使われているのは、密度汎関数法である。価電子ドブロイ波長よりはるかに大きい分子では、離れた部分の量子相関は無視できるので、近似的な独立性を使った効率良い解法（linear-scaling法）がほぼ確立されている。他方、近い部分の量子相関は、Gauss型基底関数を使い、波動関数を展開する方法が、不均一系でも効率の良い方法として確立されている。密度汎関数法を解く際必要になる、電子間反発積分や交換相関項の計算法もほぼ確立されている。つまり密度汎関数法では、理論やアルゴリズムが成熟し、これらの改良の余地は少ないため、計算機自体を工夫して高速化する強い動機がある。逆に、明らかに改良の余地ある理論を、効率良く実行するアルゴリズムや実装を考えても、理論の改良によりその努力が無に帰する可能性が大きい。 そこで私は業界標準であるGaussian03プログラムを改造し、Gauss基底を使ったlinear-scaling 密度汎関数法を加速する事にした。Gaussianの密度汎関数法の計算速度や、分子の安定構造探索手法などは、改造に良く使われるGAMESSより優れている。元々遅いプログラムやアルゴリズムをGPUで加速する事には意味が無いので、我々は最良と思われる、高速多重極子法とdirect J engineを使う事にした。その結果Gaussian03の密度汎関数計算を、最速の汎用CPUに比べて10倍に加速できた。 電子間反発積分密度汎関数計算には多数の複雑なステップがあるが、最も良く使われる500〜1000基底関数で計算時間を食う部分は、(1) 近距離Coulombポテンシャルの計算（15%）、(2) 交換相関ポテンシャルの計算（80%）、(2) Fock行列対角化（5%）なので、この部分のみ加速すれば良い。この内streaming processorが苦手なのは(1)である。ここではGauss関数間の電子間反発積分を計算する。漸化公式を使い短縮Gauss関数を同時処理する方法が、理論計算量が最小で、汎用CPUでは最も高速である。この方法では1個の電子間反発積分を計算するのに、数千語の作業領域が必要な事が、第一の問題である。Streaming processorでは僅かな作業領域しかないが、direct J法で静電ポテンシャルを計算する時に限り、primitive Hermite Gauss 関数を使っても計算速度は下がらず、かつ作業領域は数十語に減る。そこで最も作業領域の少ないGauss-Rys求積法を使った。 2つめの問題は単精度である。大型分子で計算される電子間反発積分を調べると、ごく少数で大きな値の積分が、誤差の大部分を作る事が分かった。現在の積分計算アルゴリズムでは、Schwarz不等式を用いて積分上限値を求め、十分小さい値は計算しない。そこで単精度計算での誤差が無視できる項のみGPUで計算し、残りはCPUで計算した。9割の積分を単精度で計算しても、精度に問題が生じない事が分かった。積分計算の各ステップで必要精度を調べ、数表の補間公式を与えた。精度上危険な部分が分かったので、この部分のみ倍精度計算する事もできる。 3つめの問題はCPU-グラフィックボード間の通信の遅さである。1積分あたりの計算量は数十FLOPSなので、計算した積分をグラフィックボードから回収する事はできない。そこで静電ポテンシャル行列に縮約して回収した。これは古典2体力計算を加速するGRAPEシリーズと同じアイデアである。 最後に実装上の問題として、十分な細粒度並列性を確保するため、original programでshell4つ組を作る部分を、相当修正する必要があった。またNVIDIAのコンパイラーの吐くkernel codeは最適化が不十分だったので、手探りで最適化した。 交換相関項これは密度汎関数計算で、8割の時間を占める。交換相関項は電子密度の複雑な関数なので、3次元グリッドで求積する必要が有る。グリッド上での基底関数値の計算、グリッド上での電子密度の計算、交換相関ポテンシャルの行列要素の計算が必要である。無視できる程小さい項は計算しないので、不規則なデータが多く、SIMD processorで効率良く処理する方法が問題だった。十分な細粒度並列性を確保するため、original programを相当書き換えた。また各手順は単純だが、短縮長や基底関数の数など、可変長ループが複数あり、作業領域もやはり足りない。どの項を保存し、どの項を再計算すると最も高速か、何通りかプログラムを作り比較した。kernel codeは手探りで最適化した。ハードウエアの制約に由来する、一般性が乏しく、論文に書けないし故に評価されない苦労が、沢山あった。単精度の問題は、計算すべき項を、解析計算する項と、数値計算する項に分け、後者のみ単精度で数値計算した。その結果精度の問題は生じない事が分かった。 関連研究NVIDIAがCUDAをリリースする1年前には、私はGRAPE-DRで密度汎関数計算を行うプログラムを既に完成し、シミュレータでテストしていた（GRAPE-DRボードはまだ開発中だった）。またGeForce 8800 GTXの構造は、GRAPE-DRとかなり似ていた。そのためCUDAリリース後、直ちにGPUで密度汎関数計算を行えた。論文出版の半年後に、MartinezらもGPUで電子間反発積分を計算する試みを発表した。NVIDIAのコンパイラーやライブラリーが勇気を与えてくれたと思うが、この様な煩わしい仕事をする研究者が、他にもいた事は驚きだった。彼らの研究は雑誌の表紙を飾ったようだ。彼らの研究では、我々と次の点が異なる。
またVogtらはresolution of identity (RI) MP2法を、GPUで加速する方法を報告している。RI-MP2法の計算時間の5-8割を占める行列積の計算を、NVIDIAが公開したcuBLASライブラリーで単精度計算し、精度と計算時間を調べたものである。MP2エネルギー自体が小さいので、精度に殆ど問題が生じなかった。また2-4倍に加速できた。難しい電子間反発積分は汎用CPUで計算しており、私やMartinezらよりずっと安直だが、簡単に加速できるのは結構な事だ。 最後に、情報科学でも量子化学でも余り評価されない気がする事が、この研究を進める上での最大の障害だった。 参考文献
|