分子軌道計算専用計算機用LSI (ERIC) の開発

中村 健太, 本田 宏明, 梅田 宏明, 小松 晴信, 村上 和彰


Return

1 はじめに

近年、化学や薬学などの理論や実験の多方面の研究現場において分子軌道法プログラムが頻繁に使用されており、分子の電子状態の計算のみならず、その成果の無機、有機、生化学や創薬などの広汎な分野への応用が進んでいる。しかしながら非経験的分子軌道法は、基底関数の数の4乗に比例する演算量と補助記憶量を必要とするために、取り扱う分子系が大きくなるに従い計算に必要な分子積分データの規模が急激に膨大なものとなり、現在望まれている生体系などの大規模分子系への適用は困難な状況にある。
実際、ここ数年のスーパーコンピュータや高性能ワークステーションのような大規模並列計算機の進歩により、ようやく(生体)高分子を意識した1,000基底を超える分子軌道計算が行われつつあるものの、これらの計算機は高価であり、かつ設備の維持管理も大変であるため、現在のところ研究者が研究室レベルで大規模計算をパーソナル・ユースに利用することは実質的に不可能である。またこれらの計算機上では、研究計算対象となる化学物質を他の研究者に対して秘密にしておくことが困難だという問題もある。そこで、我々は、「現実を反映した大規模分子系」の分子軌道計算を、「低コスト=パーソナル・ユース」で実現すべく、非経験的分子軌道法による分子軌道計算を高速に実行するEHPC (Embedded High Performance Computing) 専用並列計算機システムの開発に取り組んできた。
大規模分子系における分子軌道計算では2つの困難がある。1つは繰り返し行われる電子反発積分計算の計算回数の膨大さ、もう1つは個々の積分計算の処理の複雑さである。そこで我々が開発したEHPCシステムでは、木状に階層化されたシステム・アーキテクチャとリーフ・ノードに位置する多数のLSIにより膨大な電子反発積分計算を効率よく並列処理する。また、リーフ・ノードのLSIとして小規模な計算機システムへ多数の LSIを搭載できるように低消費電力化を図り、かつ、計算の特徴に最適化した高速な電子反発積分計算 (Electron Repulsion Integral Calculations) 専用 LSI ERICを開発することで、分子軌道計算プログラム全体の高速化を図る。
本稿では、電子反発積分計算の特徴解析に基づき我々が開発した電子反発積分計算専用LSI ERIC の構成について紹介するとともに、ERICを搭載したEHPCシステムの性能評価を行う。

2 電子反発積分計算アルゴリズムとその特徴

2. 1 分子軌道計算アルゴリズム

非経験的分子軌道法の解法には、Equation 1で表されるハートリー・フォック法を用いる[1]。

ここで Fi j(1), Fi j(2), C, e はそれぞれ、1電子フォック行列要素、2電子フォック行列要素、分子軌道の係数行列、固有値である軌道エネルギー行列を表しており、このうち Fi j(2) は以下の式により計算される。

なお (i j, k l), Pi jはそれぞれ、電子反発積分(ERI)、および、密度行列要素を表しており、{c}は基底関数組、N は基底関数の数を表す。全電子波動関数とエネルギーは最終的に {C} と必要な分子積分より計算される。F(2)P 行列と C から計算されるが、上記ハートリー・フォック方程式 (Equation 1)により再び C を求めるため、SCF (Self-Consistent Field) と呼ばれる繰り返し法により数値的に解かれ、この手続きは以下のようになる。
  1. 準備
  2. 分子座標や基底関数組のデータの読み込みと処理
  3. 1電子電子積分計算 (ONE-EIs)
  4. 電子反発積分より F(2) の生成 (ステップ 3 ~ 4 は繰り返し)
  5. フォック行列対角化、収束の有無のチェック、密度行列の更新
  6. 分子プロパティの計算
Table 1 は5つのペプチド分子G, GA,GAQ,GAQM,GAQMY (注:G: Glycine, A: Alanine, Q: Glutamine, M: Methionine, Y: Tyrosine) について、非経験的分子軌道計算におけるそれぞれの計算ステップの時間分布の割合を示している。SCF の計算においては F(2) 行列の生成の1計算ステップに 98% 以上の計算時間が費やされていることが分かる。
さらにこの傾向は対象とする分子が大きくなるに従い顕著となり、全計算時間の 99% 以上が F(2) 行列の生成に費やされる。そのため、この F(2) 行列生成の処理時間を大幅に削減することが非経験的分子軌道計算全体の高速化につながると言える。そこで、以下ではフォック行列生成アルゴリズムについて見ていき、その特徴に最適化した分子軌道計算専用LSIについて検討を行う。

2. 2 フォック行列生成アルゴリズム

前節でみたハートリー・フォック分子軌道法計算の律速段階ともいえるF(2) 行列生成のために Figure 1に示している積分駆動型アルゴリズム (Integral Driven Algorithm)[2] は有用な方法であり、現在種々の量子化学計算パッケージに組み込まれている。巨大分子系のフォック行列計算では必要な電子反発積分の量が莫大となるため、ハードディスクに保存できないという制約がある。そこで F(2) 行列の必要部分行列を計算するために、積分 (i j, k l)の添え字に関する対称性を考慮した上で、独立な積分計算を 1回のフォック行列計算において重複なく計算する手法がとられている。この電子反発積分計算の個数は N4 に比例しており、N の増加に従い計算量が爆発的に増大することがわかる。この単純な4重ループ構造が示している通り、i, j, k, l の一組の計算は他の組には依存していない。そのため、それぞれの電子反発積分計算には高い並列性があるといえる。

Table 1. Examples of ab initio MO Computation Timea (sec.)


Figure 1. Integral Driven Algorithm for the Fock Matrix Generation

2. 3 新小原積分法による電子反発積分の計算

前節で電子反発積分の N4 に比例する莫大な計算量について説明したが、1つの電子反発積分を計算すること自身が複雑で大変であることも広く知られている。電子反発積分のアルゴリズムは数種類あり、そのうち有名ないくつかのアルゴリズムについては現在広く利用されている。その中でも本研究では、縮約ガウス型関数にもとづく一般漸化表式の手法である新小原積分法[3, 4]を選んだ。この手法は積分の高速で高精度な計算に向いており、しかも電子反発積分のみならず、種々の分子積分法に適用できるものである。Figure 2はこの新小原積分法のアルゴリズムの電子反発積分におけるプログラム構造を示したものである。


Figure 2. Obara Algorithm for the Electron Repulsion Integral Calculation

新小原積分法の計算手法は初期積分計算部分と漸化計算部分に明確に分かれており、基底関数縮約のための縮約ループは初期積分計算内部にのみ現れることが特徴的である。

2. 3. 1 初期積分計算部分

初期積分計算は Figure 2に示す 4重のループ構造をとる。このループ1回あたりの演算数は非常に少なく、命令レベル並列性 (ILP: instruction-level parallelism) は低い。このため、複数の演算回路を用いた並列計算による高速化は期待できない。さらに、浮動小数点除算、開平逆数演算、指数関数演算、および、Equation 5で表される誤差関数演算といった複雑な倍精度の浮動小数点演算を含んでおり、かつ、それらの演算がクリティカル・パス上に存在するため、4重ループ全体としては非常に負荷の高い処理となっている。したがって初期積分計算部分を高速に処理するには、個々の複雑な演算を高速に処理することが重要となる。

2. 3. 2 漸化計算部分

漸化計算は初期積分計算が完了したのち、初期積分の計算結果をもとにして最終的に求める電子反発積分の値を計算するステップである。Figure 3 は漸化計算プログラムで実際に使われているコードの一部を表しているが、これから見て取れるように漸化計算は多数の積和演算命令のみから構成されている。また、非常に高い命令レベル並列性を持つことから、複数の積和演算回路を用いた並列計算による高速化が期待できる。


Figure 3. An Example of the Recurrence Calculation Program

3 分子軌道計算専用計算機システム・アーキテクチャ

EHPC (Embedded High Performance Computing) プロジェクトでは、非経験的分子軌道法による巨大分子の電子状態計算のためのEHPC専用並列計算機システムを構築するに際して、Figure 1で示した積分駆動型アルゴリズムによるフォック行列生成を効率よく処理するために、以下の方法論を採る。
  1. 木状に階層化されたシステム・アーキテクチャにより、多数の積分計算をリーフ・ノードに位置するプロセッサに割り当て、効率のよい並列処理を行う
  2. 積分計算の特殊性に特化した電子反発積分計算専用LSI ERICにより、個々の積分計算処理の高速化を図る
これを実現するために、我々は Figure 4 に示す EHPC 分子軌道計算専用並列計算機システム・アーキテクチャを開発した。
この専用並列計算機システムは、以下の4層からなる階層的なアーキテクチャを採る。
Bottom Layer:
積分計算の特徴に最適化された ERIC LSI で、木構造のリーフ・ノードに位置する。後述のSH-4プロセッサと共に、組み込みシステム業界標準CompactPCI規格に準拠した専用計算ボード(Figure 5)上に搭載される。
House Keeping Layer:
複数個の ERIC LSI を搭載した専用計算ボード上で、上位のSecond Layerからの初期データを受け取り、ERIC LSI上のメモリ空間に書き込む一方、専用LSIで得られた計算結果の吸い出しを行う。組み込みシステム業界で用いられている、株式会社ルネサス テクノロジ製のSH-4プロセッサを用いる。
Second Layer:
CompactPCI規格に準拠した汎用 CPUボードであり、PentiumIIIプロセッサを搭載したLinux PC である。Figure 6の様に、CompactPCI規格に準拠した筐体に最大 7枚の専用計算ボードと共に格納され、計算ボードと CompactPCIバスで接続される。
Top Layer:
Linuxを搭載したPCであり、下層に位置する複数の汎用 CPUボードと100Mbps Ethernetで接続される。CPUボードと併せてPCクラスタを構成し、専用並列計算機システム全体のユーザー・インタフェースとして使用する。


Figure 4. EHPC Platform Architecture


Figure 5. ERIC LSI Board


Figure 6. EHPC Parallel Computing System

木構造をもつシステム・アーキテクチャでは一般に、同一階層に位置し直接接続されていないリーフ・ノード間での通信効率が低いという問題がある。しかしながら、本並列計算機システムが対象としているフォック行列生成処理で必要なデータ通信は、主に以下の 4通りである。

  1. 初期の準備処理として全てのリーフ・ノードに対し、256KBytes以下の計算のための諸データと 200MBytes程度の密度行列をブロードキャスト
  2. 全計算の i, j, k, li, j, k, l それぞれが取り得る値は 1〜5,000)のうち、i, j のペア・データ(4*2Bytes以下)を House Keeping Layer に位置する SH-4 へ転送
  3. k, l のペア・データ(4*2Bytes以下)について、SH-4から各リーフ・ノードに対して転送
  4. 全ての計算が終わった段階で、各リーフ・ノードより 200MBytes程度のフォック行列を、SH-4を介して部分和を求めながら Top Layer に転送
ここで、Figure 1における F(2) 行列への累積計算は、データ通信量の削減を図るべくリーフ・ノード内で処理することとしている。これから判るように、大規模なデータ転送は計算処理の最初と最後のみでしか発生せず、実行中に頻繁に発生する通信は 4*2Bytes以下のショート・メッセージのみである。また、リーフ・ノード間相互での通信も発生しない。このため、高度な機能を持つ通信路は必要とせず、木構造を持ったシンプルなシステム・アーキテクチャが有効に機能する。

4 ERIC: 電子反発積分計算専用LSI

4. 1 設計方針

我々はまず最初に、電子反発積分計算の特徴が異なる2つの計算部分、初期積分計算と漸化計算の両方の特徴を考慮し、電子反発積分計算専用LSI ERIC の構成として以下のような基本方針を立てた。
  1. 初期積分計算に含まれる、倍精度浮動小数点除算、開平逆数演算、指数関数演算、誤差関数演算といった複雑な算術演算を高速に処理するため、専用算術演算回路を設ける。
  2. 漸化計算に含まれる非常に高い命令レベル並列性を活用するため、多数の積和演算器を設け、並列計算を行う。
しかしながら、初期積分計算は命令レベル並列性が低いことから、多数の積和演算器を設けても初期積分計算部分の処理の高速化は期待できない。一方、漸化計算は積和演算のみで構成されるため、漸化計算の処理中は複雑な算術演算専用回路を全く必要としない。そのため、この2つの機能を単一のプロセッサ上に実装するのは効果的ではない。そこで我々は、ERIC LSIのアーキテクチャとして、初期積分計算エンジンと漸化計算エンジンの 2つの計算エンジンにから構成される、チップ・マルチプロセッサ・アーキテクチャを採用することとした[5, 6]。

4. 1. 1 初期積分計算エンジン

初期積分計算は、整数演算、倍精度浮動小数点加減算、乗算、除算、開平逆数、指数関数演算、関数呼び出しや条件分岐といった、多種多様な命令を含む。そこで初期積分計算エンジンは、汎用のRISCプロセッサをベースとし、複雑な算術演算のための専用演算回路を付加することで、初期積分計算部分を高速に処理するアーキテクチャとした。

4. 1. 2 漸化計算エンジン

漸化計算は、多数の積和演算のみで構成されている。また、それぞれの積和演算の間には依存関係が少なく、同時並行的に処理可能なものも多い。そこで我々はこの特徴に最適なアーキテクチャを検討し、その結果として複数のマイクロ・エンジンを持つ、マイクロ・エンジン・ベースのチップ・マルチ・プロセッサ(mCMP)アーキテクチャ[7]を開発、漸化計算エンジンに採用することとした。ここで、マイクロ・エンジンとは、アプリケーションに特化し、機能を最小限に絞り込んだ小さいプロセッサのことである。

4. 2 ERIC LSIアーキテクチャ

前節で示した基本方針に加え、チップサイズを勘案しつつ必要なメモリサイズや搭載可能な演算器の個数などの検討を重ね、最終的なERIC LSIアーキテクチャは、1つの初期積分計算エンジン(IIC Engine)と4つのマイクロ・エンジンからなる漸化計算エンジン(RC Engine)の 2つのエンジンと、704Kbytesのオンチップ SRAMからなるFigure 7のようなアーキテクチャに決定した。
初期積分計算エンジンは、ERIC LSIの共同開発を行ったセイコーエプソン株式会社が持つ、独自の 32bit RISCマイコン S1C33 Familyをベースとし、倍精度浮動小数点数演算のための、積和演算回路、除算/開平逆数演算回路、指数関数演算回路、誤差関数演算回路を付加した。
漸化計算エンジンは、浮動小数点積和演算回路、演算に必要なデータを主記憶上から取得し、演算結果を主記憶に書き戻すためのロード/ストア・ユニット、アドレス計算用ALU(加減算のみ)、および、レジスタファイルのみで構成され、サポートする命令数も20種類程度のプロセッサを設計し、マイクロ・エンジンとして採用した。
また、設計したERIC LSIのチップ試作を行った。Figure 8は完成したERIC LSIのチップ写真である。

5 分子軌道計算専用並列計算機システムの性能

5. 1 ERIC LSIの主な仕様

今回、我々が設計したERIC LSIの諸元は Table 2に示すとおりである。また、ERIC LSI に搭載されている、倍精度浮動小数点演算ユニットにおける各演算処理時間を Table 3に示す。これらの演算回路のうち加算回路 (addition)、および、乗算回路(multiplication) ではパイプライン化を行っている。それ以外の演算回路については、積和演算の繰り返しにより値を求めるアルゴリズムを採用していることによる、実装上問題からパイプライン化を施していない。


Figure 7. Overview of the ERIC LSI Architecture

Table 2. Major Specifications of the ERIC LSI

Table 3. Execution Clock Cycles of each Floating-Point Arithmetic Operation

これらの浮動小数点演算命令は電子反発積分計算の高速な処理に重要である。例えば誤差関数計算を一般的なMIPSプロセッサ上でテイラー展開を用いたソフトウェアで実現すると 32命令、146クロック・サイクル掛かるが、ERIC LSIでは誤差関数 (error function) 専用演算回路を用いることで 1命令、16クロック・サイクルで完了することができる。この誤差関数計算は、フォック行列生成の際の電子反発積分計算を行う4重ループの最内側で高頻度に実行されるため、専用演算回路を設けることがフォック行列生成全体の高速化につながる。同様に、倍精度浮動小数点除算 (division) 回路、開平逆数 (square root) 演算回路、指数関数 (exponential function) 演算回路を専用のハードウェアとして実装したことにより、同様の処理をソフトウェアで実行する場合と比べ、大幅に実行時間を削減することができる。そのため ERIC LSI を用いることにより、効率のよい電子反発積分計算を行うことが期待できる。

5. 2 分子軌道計算専用並列計算機システムの見積もり性能

リーフ・ノードの計算エレメントとして SH-4を使用した分子軌道計算専用並列計算機システムのプロトタイプ・システムを使った並列性能評価では、63プロセッサを使用する 3ユニット・システムに於いて 62.6倍の高速化が報告されている[8]。このリニアな速度向上は、Figure 1のフォック行列生成アルゴリズムが備えている並列化の容易さという特色だけに起因するものではなく、EHPC専用並列計算機システム・アーキテクチャに依るところが大きい。


Figure 8. Chip Implementation of the ERIC LSI

一般に木構造を持つシステム・アーキテクチャでは直接通信ができないリーフ・ノード間での通信オーバーヘッドが大きく、並列化効率を高めるには計算負荷の分散とともに通信の局所化が必要となる。これを実現するには、木構造の中間ノード全てにジョブマネージャを置く多段の動的負荷分散機構が有効である。本システムで採用している、階層的なシステム・アーキテクチャはこれを狙ったものであり、さらにソフトウェアの最適化もあって高い並列化効率を実現している。これは、専用LSIを使用した場合でも同じことが実現できるため、最終的な分子軌道計算専用並列計算機システムでも 99%以上の高い並列化効率が期待できる。
前節の諸元より、ERIC LSI 1チップあたりの消費電力は最大 4.2W である。現在までに性能測定が行われている SH-4を搭載したEHPC並列計算機システムから見積もると、ERIC LSI をリーフ計算ノードとして置き換え、56個 のERIC LSI を搭載した4ユニット・システム構成での全消費電力は 1KW 以下になると予想される。また、Figure 9が SH-4を搭載したEHPC並列計算機システムの 4ユニット・システム構成の写真であり、この体積は1200×650×1010mm 程度であるが、ERIC LSIを搭載したシステムでも内部のボードが置き換わるのみの変更であるため、外形寸法は同じである。このことから本並列計算機システムは、電力ならびに体積の観点から小さな研究室でも設置し、維持管理できる規模の計算機システムとなっているという特徴もある。更には、Table 3 に示したように誤差関数回路や指数演算回路など電子反発積分頻出の特殊演算器を積んだ ERIC LSI の特徴から、大規模な分子軌道法計算に対し効率的で低コストな並列計算機システムを提供することができると言える。

6 まとめ

我々は、「現実を反映した大規模分子系」の分子軌道計算を、「低コスト=パーソナルユース」で実現すべく、EHPC専用並列計算機システムを開発した。非経験的分子軌道法において全計算時間の 98%以上を占めるのがフォック行列生成のステップであるのため、分子軌道計算全体の高速化を図るにはフォック行列生成アルゴリズムに着目し、並列処理を行うことにより高速化を図るのが効果的である。


Figure 9. EHPC Parallel Computing System (Four Units System)

本研究では、フォック行列生成に必要な膨大な数の電子反発積分計算を効率よく並列処理するため、木状に階層化されたアーキテクチャを採る EHPC専用並列計算機システムを構築すると同時に、個々の電子反発積分計算の高速化を目的とした電子反発積分計算専用LSI ERICを開発した。ERIC LSIは、積分計算アルゴリズムの解析結果に基づいた最適化を行い、倍精度浮動小数点除算、開平逆数演算、指数関数演算、誤差関数演算のための専用演算回路を持つと共に、多数の積和演算を高速に並列処理するための4つのマイクロエンジンを持つアーキテクチャを採っている。
ERIC LSIは、EHPC専用並列計算機システムに多数個組み込むことを前提とし、低消費電力な LSI としたため、外形寸法が 1200×650×1010mm のシステムに56個のERIC LSIを搭載することができ、その消費電力も 1KW以下になると予想される。これにより、大規模な分子軌道計算に対し効率的で、低コストな並列計算機システムを提供することが出来る。

本研究の一部は科学技術振興調整費 総合研究「科学技術計算専用ロジック組込み型プラットフォーム・アーキテクチャに関する研究」(代表 村上和彰 九州大学教授)によるものである。

参考文献

[ 1] A. Szabo and N. S. Ostlund, Modern Quantum Chemistry, Dover Publishers (1996).
[ 2] J. Almlof, K. Faegri, Jr., and K. Korsell, Principles for a Direct SCF Approach to LCAO-MO Ab Initio Calculations, J. Chem. Phys., 3, 385-399 (1982).
[ 3] H. Honda, T. Yamaki, and S. Obara, Molecular integrals evaluated over contracted Gaussian functions by using auxiliary contracted hyper-Gaussian functions, J. Chem. Phys., 117, 1457-1469 (2002).
[ 4] H. Honda and S. Obara, Molecular Integrals evaluated over Contracted Gaussian functions, Proceedings of the 11th International Congress of Quantum Chemistry (2003).
[ 5] K. Nakamura, H. Hatae, M. Harada, Y. Kuwayama, M. Uehara, H. Sato, S. Obara, H. Honda, U. Nagashima, Y. Inadomi and K. Murakami, Eric: A Special-Purpose Processor for ERI Calculations in Quantum Chemistry Applications, HPC-Asia 2002 (2002).
[ 6] 原田宗幸, 中村健太, 桑山庸史, 上原正光, 佐藤比佐夫, 小原繁, 本田宏明, 長嶋雲兵, 稲富雄一, 村上和彰, 二電子積分計算専用プロセッサ・アーキテクチャの開発, 情報処理学会論文誌:ハイパフォーマンスコンピューティングシステム, 44, No. SIG 1 (HPS6) pp.1-9 (2003).
[ 7] K. Nakamura, H. Honda, K. Inoue, H. Sato, M. Uehara, H. Komatsu, H. Umeda, Y. Inadomi, K. Araki, T. Sasaki, S. Obara, U. Nagashima and K. Murakami, A HighPerformance, LowPower Chip Multiprocessor for Large Scale Molecular Orbital Calculation, Proceedings of The Workshop on Unique Chips and Systems (UCAS-1), 87-94 (2005).
[ 8] H. Umeda, Y. Inadomi, H. Honda and U. Nagashima, Parallel Fock matrix construction on layered multi-processor system, Proceedings of The 229th ACS National Meeting (2005).


Return