分散メモリ型並列計算機に適した新しい大規模フォック行列生成アルゴリズム
- 積分カットオフとの関連 -

高島 一, 山田 想, 小原 繁, 北村 一泰, 稲畑 深二郎, 宮川 宣明, 田辺 和俊, 長嶋 雲兵

1 はじめに

Figure 1. MOEngine, a special purpose machine for accelerating molecular orbital calculations.

現在開発中の専用計算機MOEngineでは、安価なコストで高速計算を実現するため搭載メモリ量（や設備の規模）を小さくしている。そのため各並列プロセッサに全ての行列要素を搭載できないので、MOEngineで大規模フォック行列を計算するためには、ホストと各プロセッサの間で行列要素の通信を行なう並列分散アルゴリズムが必要になる。その際、大規模分子軌道計算で重要な問題となる二電子積分のカットオフやそれに伴う行列要素転送のカットオフに関する考察が十分でなければ実用性に乏しくなる。これまでにフォック行列の並列分散アルゴリズムはいくつか[3, 4]提唱されてはいるが、いずれもアーキテクチャを特定していたり通信面やカットオフの面の考察が十分でなかったりするため、蛋白質の様な大規模計算に適用するのは難しい。
そこで我々は、各専用計算機MOEngineのような各プロセッサのメモリが小容量（数MByte）しかない分散メモリ型並列アーキテクチャのもとで大規模フォック行列を計算するための新規アルゴリズムを考案した。本論文ではまず、大規模分子系で重要となる積分カットオフについて、以下の順序で説明する。最初に、分子軌道計算及び二電子積分とフォック行列の計算方法について概要を示す。次に、基底関数について、「縮約(contracted)」や「原始(primitive)」など以後の議論の理解に必要となるいくつかの言葉の定義を説明しておく。そして、二電子積分のカットオフの方法、並列フォック行列アルゴリズムにおける行列要素転送のカットオフとの関係、及び、大規模計算における積分カットオフの重要性、について詳細に述べる。なお、続く第二報[17]では、我々がRT並列アルゴリズムと呼んでいる新規フォック行列計算アルゴリズムの詳細について述べる。

2 非経験的分子軌道法について

Figure 2. Typical flow charts for an ab initio molecular orbital calculation.

Table 1. Computing times in eeach SCF calculation step for some peptide molecules(sec).

	G	GA	GAQ	GAQM	GAQMY
Number of dimensions	55	110	207	316	427
Initial guess	0.1	0.7	4.4	15.0	43.7
One-Electron Integral	0.1	0.5	2.3	7.8	14.4
Two-Electron Integral and Fock Matrix	21.9 (96.9%)	268.5 (98.8%)	2241.3 (99.1%)	10994.7 (99.4%)	32624.7 (99.4%)
Diagonalization	0.1	1.0	6.9	28.4	103.3
Properties	0.1	0.5	2.8	9.7	24.2
Others	0.3	0.6	3.0	9.5	25.1
Total	22.6	271.8	2260.7	11065.1	32835.4

G: Gly, GA:Gly-Ala, GAQ: Gly-Ala-Gln, GAQM: Gly-Ala-Gln-Met, GAQMY: Gly-Ala-Gln-Met-Tyr

3 基底関数

4 積分カットオフ

4. 1 概論

4. 2 並列分散フォック行列アルゴリズムとの関係

for(I=1;I≦N;I++) { for(J=1;J≦N;J++) { for(K=1;K≦N;K++) { for(L=1;L≦N;L++) { compute (IJ,KL), (IK,JL); FIJ += PKL*(IJ,KL); FIJ -= 0.5*PKL*(IK,JL); }}}}

Figure 3. Index-driven Fock matrix algorithm

これは、式(2.10)に示したフォック行列二電子項をそのまま素直に展開した方法であり、添字駆動型アルゴリズムと呼ばれる。このアルゴリズムを並列分散化する場合について考える。必要に応じてホスト計算機との間で行列要素データの送受信を行なう事になるので、データの転送回数を減らし同じ行列要素を使いまわして利用するためには、行列要素の転送をできるだけ外側のループで行なうようにすればよい。一例をFigure 4に示す。
Sendはホスト計算機からプロセッサエレメントへのデータ送信を、receiveはプロセッサエレメントからホスト計算機へのデータ受信をそれぞれ表す。プロセッサエレメントで必要な行列要素はF_IJとP_KLだけである。よって、Figure 4に示すように、四重ループの順序をI,J,K,LをからI,L,J,Kに入れ替えることにより、行列転送をJ,Kの二重ループの外側に配置し、O(N²)回の積分計算に対しO(N)個の行列転送量に抑える事ができる。
次に、Figure 4のアルゴリズムを用いて計算時間の短縮化のために二電子積分のカットオフを行なうことを考える。計算すべき二電子積分は(IJ,KL)と(IK,JL)であり、行列情報を送受信する時点ではすでにIとLの値は定まっている事から

for(I=1;I≦N;I++) { for(L=1;L≦N;L++) { send PKL(1≦K≦N) for(J=1;J≦N;J++) { for(K=1;K≦N;K++) { compute (IJ,KL), (IK,JL); FIJ += PKL*(IJ,KL); FIJ -= 0.5*PKL*(IK,JL); }} receive FIJ(1≦J≦I) }}

Figure 4. Index-driven Fock matrix algorithm for distributed memories.

・(IJ,KL)からは、Iに対するカットオフで生き残るJの組と、Lに対するカットオフで生き残るKの組が、
・(IK,JL)からは、Iに対するカットオフで生き残るKの組と、Lに対するカットオフで生き残るJの組が
それぞれ分かる。この場合、上記二つの条件から、F_IJの生き残る添字Jの組や、P_KLの生き残る添字Kの組はいずれも、Iに対するカットオフで生き残るものとLに対するカットオフで生き残るものの論理和で与えられる。従って、JとKのループでは、1〜Nの全てをではなく、カットオフ操作で生き残る添字だけを実行すればよい。
カットオフ操作でどの添字が生き残るかについては、あらかじめ、全てのインデックスに対してカットオフでどの添字が生き残るかというリストを示した「カットオフテーブル」を作成しておくと便利である。この計算量はO(N²)であり、またSCF計算の間で変化する事はないので、式(4.6)を用いて最初に一度ホスト計算機で作成しておけばよい。カットオフテーブルの具体例をFigure 5に示す。この例では、基底数は6であり、基底番号1〜6のそれぞれに対してカットオフ操作で生き残る添字のリストが与えられている。
カットオフ操作を行なう場合にFigure 4のアルゴリズムに沿ったフォック行列作成方法をFigure 6に示す。また、その手順を、Figure 5のカットオフテーブルを用いながら具体的に説明する。
0)ホスト計算機でカットオフテーブルを作成しておく。
1)プロセッサエレメントにI,Lを割りあてる。
例えば、I=3,L=2とする。

Figure 5. Cutoff survival index

Figure 6. Flow charts of parallel Fock matrix algorithm for distributed memories.

2)割り当てられたIとLに対して生き残るJとKの組をカットオフテーブルから読み出す。
カットオフテーブルから、
・I=3に対して生き残る添字Jは{2,3,5,6}、
・L=2に対して生き残る添字Kは{2,3,6}、
であることが分かる。よって、これらの論理和である{2,3,5,6}が、Figure 4のアルゴリズムで実際に取りうるJとKの値になる。
3)読み出したリストに従ってカットオフされた密度行列P_KLの情報を作成する。
P_*2(1≦*≦6)から、P₂₂, P₃₂, P₅₂, P₆₂を取り出す。
4)作成したP_KLの情報をプロセッサエレメントに転送する。
5)プロセッサで二電子積分を求め、フォック行列要素F_IJを計算する。
作成されるフォック行列要素は、F₃₂, F₃₃, F₃₅, F₃₆である。
6)計算されたフォック行列要素F_IJの情報をプロセッサエレメントから受け取る。
7)全てのI,Lペアが終了していなければ、1)〜6)を繰り返す。
このように、積分カットオフの判断基準に従ったカットオフテーブルを参照することによって、カットオフを考慮した必要最小限の密度行列およびフォック行列をホスト計算機側で準備してプロセッサエレメントと通信する事が可能になる。これにより、ホスト計算機とプロセッサエレメント間の転送量やプロセッサエレメントのメモリ量を大幅に減らす事ができる。
カットオフで生き残るKの添字はKLカットオフとIKカットオフの論理和で与えられるので、カットオフ生き残り率をa（0以上1以下の実数）とすると、転送すべき密度行列D_KLの要素数は2aNで与えられる。よって、本アルゴリズムを用いた場合の、二電子積分１個当りの行列要素の平均転送量（ここでは仮に[Word/ERI]とする）は、

となる。これが、仮に行列転送のカットオフができないとすると、二電子積分一個当たりの行列転送量は 4/(a²N) [Word/ERI]となり、カットオフ可能な場合に比べて 1/a倍の転送量（速度）とメモリ量が要求されることになる。4.3節で明らかにするように、大規模系になるほどaの値はどんどん小さくなるので、大規模分子を計算する場合にはカットオフの効果を考慮する事が重要であることが分かる。
このように、分散小メモリ型アーキテクチャにおける大規模並列フォック行列計算アルゴリズムを考えるときには、従来のように積分の対称性や積分のカットオフを有効に利用して計算量の低減化を図る以外に、
・行列データの送受信は二重ループの外側で行ない、O(N²)回の積分計算にO(N)個の行列要素転送を行なう、
・転送する行列要素もカットオフを行なう、
ことが大事な点である。つまり、一度送信した行列要素は何回も計算に用いると同時に、計算に用いられない行列を転送するような無駄を省くことが、通信量の軽減とメモリ量の縮小の両方の点において重要である。次節では、蛋白質のような大規模分子において、積分カットオフ及びそれに伴う行列転送カットオフの割合が実際どの程度の大きさになるかついて検討する。
なお、Figure 4に示すアルゴリズムの欠点は、式(2.13)に示す積分の対称性が利用されていないため、重複した積分計算回数が多くなることである。この点を改善したアルゴリズムとして、我々はRT並列アルゴリズムを開発した。その詳細は第二報[17]で示す。

4. 3 大規模計算におけるカットオフの実例

Figure 7. Energy Changes. (a) Relationship between density change (DP) and energy change (DE; kcal/mol) in SCF iteration process. Fill mark denotes that the energy change is positive.(b) Changes in total energies with various integral cutoff thresholds.

従って、10^-15程度を閾値とした二電子積分のカットオフを行なうことにより、全エネルギー値は十分計算精度を保つ事ができるのがわかる。これは、分子軌道計算が有効桁数16桁の倍精度浮動小数点演算で行われており、フォック行列の値がおおよそ0.01から100程度の範囲であることを考え合わせると妥当な結果であるといえる。
次に、積分カットオフ閾値を10^-15とした場合に、実際の大規模分子でどの程度カットオフが行われるのかについて調べた。縮約基底関数で表現される二電子積分(IJ,KL)はN⁴個あるが、カットオフを考慮した後に計算すべき二電子積分の数は、

になる。そこで、小分子から蛋白質程度の大きさまでさまざまなサイズの分子で、カットオフ生き残り率aおよび生き残った基底数Mを見積もった。結果をFigure 8とTable 2に示す。
縮約基底関数で表現される二電子積分(IJ,KL)は、式(4.1)で示す通り原始基底関数で表現される二電子積分g(i,j,k,l)の和である。一方、積分がカットオフできるかどうかは原始基底関数（もしくは原始シェル）に与えられる軌道指数及び軌道中心で決定される。そこで、ある縮約基底二電子積分(IJ,KL)を構成するすべての原始基底二電子積分g(i,j,k,l)がカットオフされるなら、(IJ,KL)はカットオフされると判断できる。aはこの縮約基底二電子積分に対する生き残り率を表す。

Figure 8. Effects of Cutoff. (a) Cutoff survival Ratio a. (b) Number of survival basis M.

計算に用いた分子は、H₂O、ベンゼン(C₆H₆)、アミノ酸であるArginine、GAQ(Gly-Ala-Gln)、GAQMY(Gly-Ala-Gln-Met-Tyr)、ヌクレオチドd(CGC)、Gramicidin A二量体（30残基）、BPTI（Bovine Pancreatic Trypsin Inhibitor、57残基）、Human Immunoglobulin IGG1（99残基）、Lysozyme（129残基）、の10種類である。基底関数による違いを調べるため、STO-3G、6-31G、6-31G**、分極関数やdiffuse functionを取り込んだ拡張基底(Extended Basis)として6-311++G(3d,2p)、の4種類を用いた。積分カットオフの閾値Kは10^-15を用いた。
Figure 8(a)は、基底数Nの分子におけるカットオフ生き残り率a（全ての添字での平均値）を両対数プロットした図である。小分子で基底数Nも小さい時（100以下）は、aはおおよそ1であり縮約基底関数で表現される積分のカットオフは殆どなされない (注：一つの縮約基底の中にはさまざまな軌道指数を持つ原子基底が含まれているため、原子基底で表現される二電子積分g(i,j,k,l)で見ると多くの積分がカットオフされている。) 。しかしながら、分子サイズが大きくなり基底数が増えるにつれて急速に積分カットオフがなされるようになる。これは、より中心間距離が離れたガウス関数のペアが多くなるためである。10,000次元程度になると、STO-3G、6-31G、6-31G**クラスの基底関数では、カットオフ生き残り率aは0.1以下になっている。一方、拡張基底の場合、diffuse関数や分極関数のような空間的に広がった基底関数が多数用いられているため、10,000次元でもカットオフ生き残り率aが0.3程度と、他の基底関数に比べて少し大きな値である。このように、カットオフ生き残り率aは基底関数の種類によって異なり、より精度の高く空間的に広がった軌道をもつ基底関数ほどカットオフされにくく、最小基底であるSTO-3Gのような空間的広がりを持たない基底関数系ではカットオフされる割合が大きい（カットオフ生き残り率aが小さい）事が分かる。

Table 2. Number of basis N and the number of cutoff survival basis M.

molecule	STO-3G		6-31G		6-31G**		Extended
	N	M	N	M	N	M	N	M
H2O	7	7	13	13	25	25	55	55
C6H6	36	36	66	66	120	117	270	262
Arginine	74	67	136	119	250	196	560	457
GAQ	113	105	207	180	375	263	845	671
GAQMY	235	188	427	314	769	455	1733	1176
d(CGC)	668	288	1203	479	2106	635	4797	2069
Gramicidin	1587	448	2910	710	5295	869	11910	3525
BPTI	2572	462	4694	739	8504	891	19138	4121
IGG	4334	507	7940	803	14438	962	32476	4806
Lysozyme	6004	552	10967	881	19850	1044	44705	5578

Figure 8(b)は、基底数Nに対してカットオフで生存する平均基底数Mをプロットしたものである。図中の対角線は、カットオフが全くなされない場合を表している。Table 2に数値を示す。基底数Nが小さい場合にはカットオフが殆どなくどの基底でも対角線上に近いが、基底数が大きくなると、精度の低い（広がった基底の少ない）基底関数系から順番に対角線から外れ出す。そして、5,000基底程度以上になると、カットオフ後に生存する基底数Mは拡張基底の場合を除いて殆ど増加しなくなり、ほぼ一定値になっていることが分かる。このように、蛋白分子のように基底数Nが大きい場合では、式(4.13)から判るように、計算しなければならない二電子積分の数はO(N⁴)からO(N²)になる。しかも、精度の低い基底関数ほどO(N²)になりやすい傾向が見える。先に述べた通り、精度の低い基底は電子雲の歪みを表現するための空間的に広がった関数が少ないので、他の中心を持つガウス関数との重なりが小さくなりカットオフされやすいのがその理由である。
なお、ここで示したカットオフ生き残り率aや生き残った基底数Mは分子全体での平均値であり、原子の種類や位置、縮約シェルの種類（例えば、内殻に相当する1s軌道かvalenceのp軌道か）、等によって、同じ分子内でもカットオフされる割合がかなり異なることは注意すべきである。6-31G基底におけるLysozyme分子の場合、カットオフで生き残る基底数Mの平均値は881であるが、最も少ない場合にはM=91であり、逆に最も多い場合にはM=3409であった。このように、原子の種類や位置、軌道の種類などによって、カットオフで生き残る基底数Mはばらつきがある。
では、ある軌道に対してカットオフ操作を行なった時に生き残る軌道の空間的配置はどのようなものであろうか？一例として、6-31G基底でのLysozyme分子において、カットオフ操作で生き残る軌道の空間的広がりを示したのがFigure 9である。
分子のほぼ中央に位置する56番目Leucineのd炭素原子（大きな球で表している）の2p軌道に対してカットオフを考慮した場合に、生き残る軌道が空間的にどのように分布しているかを調べた。カットオフ閾値は10^-15である。生き残る軌道を一つでももつ原子を小さな球で示しているが、これらの多くは半径10Å以内に集まっている。しかし、原子の種類によっては、例えば10Åより離れていても生き残る軌道を持つものもあるし、逆に近くにあっても生き残る軌道を全く持たないものもある。さらに、小さな球で示されている原子であっても、その原子に含まれる全ての軌道が生き残るわけではなく「少なくとも一つの軌道が」生き残っているだけであり、多くの軌道はカットオフされている。

Figure 9. Cutoff scheme in Lysozyme molecule.

このように、ある添字Iの軌道に対してカットオフにより生き残る軌道JはIの近傍の原子に集中しており、遠くの原子には殆ど存在しない。従って、蛋白質のような大規模系では、カットオフされる割合は非常に大きくなり二電子積分の数は実質上O(N²)になることが理解できる。
大規模フォック行列を計算するための並列分散アルゴリズムでは、ホスト−プロセッサ間の転送量（速度）及びプロセッサエレメントのメモリ量の観点から、積分のカットオフ及びこれに伴う行列転送のカットオフを考慮する事が大事である。実際、大規模分子でのカットオフの割合を調べる事で、カットオフを考慮する場合としない場合においては10倍以上転送量（速度）やメモリ量が異なることが分かった。
そこで我々は、上記の点を取り入れて、しかも積分の対称性もある程度考慮した新規並列フォック行列アルゴリズムを考案した。このアルゴリズムはI,Kループが最外二重ループになり、ループ内を一つのジョブ単位として各プロセッサに振り分けている。但し、実際に用いるときには、縮約基底の単位ではなく縮約シェルR,S,T,Uのインデックスを用いてループを回し、縮約シェルの添字R,Tが最外ループに位置する。そこで、このアルゴリズムを我々はRT並列アルゴリズムと呼んでいる。[17]

5 まとめ

本研究に当たりいつも有益な議論を頂いている島根大学助教授網崎孝志博士、九州大学助教授村上和彰博士に感謝する。本研究の一部は、産業技術情報基盤整備研究開発「ヘテロジニアスネットワークコンピューティング環境における分子シミュレーションシステムの構築」、科学技術振興調整費「第一原理計算手法に基づくシミュレーション技術に関する研究」、文部省科学研究費試験研究B「非経験的分子軌道計算の超高速化に基づくハードウェア及びソフトウェアの開発」、科学技術振興事業団計算科学技術活用型特定研究開発推進事業（短期集中型）「大規模分子軌道計算を高速かつ低コストで実現するための研究開発」に基づくものである。

参考文献