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のコンパイラーやライブラリーが勇気を与えてくれたと思うが、この様な煩わしい仕事をする研究者が、他にもいた事は驚きだった。彼らの研究は雑誌の表紙を飾ったようだ。彼らの研究では、我々と次の点が異なる。

  1. 彼らは「近い将来倍精度計算を高速にできるGPUが登場する」と予想し、単精度の問題、つまり誤差解析はしていない。他方私は回路規模から考えて「倍精度計算は単精度より相当遅い」と予想し、単精度で計算できる工夫をした。実際には倍精度計算速度は単精度の1/4〜1/10で、単精度計算は必須だろう。
  2. 彼らはlinear scaling (LS)法を使っていない。LS法を使うと計算すべき電子間反発積分の数が減るので、汎用CPUでは数倍速くなる。またshell 4つ組も不規則になり、streaming processorに不利になる。大型分子では普通LS法を使うので、この点は公正な比較でない。
  3. Threadへの色々なタスク割り当て方を比較検討している。

またVogtらはresolution of identity (RI) MP2法を、GPUで加速する方法を報告している。RI-MP2法の計算時間の5-8割を占める行列積の計算を、NVIDIAが公開したcuBLASライブラリーで単精度計算し、精度と計算時間を調べたものである。MP2エネルギー自体が小さいので、精度に殆ど問題が生じなかった。また2-4倍に加速できた。難しい電子間反発積分は汎用CPUで計算しており、私やMartinezらよりずっと安直だが、簡単に加速できるのは結構な事だ。

最後に、情報科学でも量子化学でも余り評価されない気がする事が、この研究を進める上での最大の障害だった。

前ページ 次ページ

参考文献

  • Two-Electron Integral Evaluation on the Graphics Processor Unit, Koji Yasuda, Journal of Computational Chemistry, 29, 334, Published Online: 5?Jul?2007
  • Accelerating Resolution-of-the-Identity Second-Order Moller-Plesset Quantum Chemistry Calculations with Graphical Processing Units, Leslie Vogt, Roberto Olivares-Amaya, Sean Kermes, Yihan Shao, Carlos Amador-Bedolla, and Aln Aspuru-Guzik, J. Phys. Chem. A, 112, 2049, 2008. Web Release Date: January 30, 2008
  • Quantum Chemistry on Graphical Processing Units. 1. Strategies for Two-Electron Integral Evaluation. Ivan S. Ufimtsev and Todd J. Martinez, Journal of Chemical Theory and Computation, February 2008
  • Accelerating Density-Functional Calculations with Graphics Processing Units, Koji Yasuda, Journal of Chemical Theory and Computation, July 2008