明 細 書
並列計算方法及び装置
技術分野
[0001] 本発明は、分子シミュレーションなどに適した並列計算方法及び装置に関し、特に 、 RHF (制限ハートリ^ ~ ·フォック; Restricted Hartree- Fock)法による分子軌道計算 に適合した並列計算方法及び装置に関する。
背景技術
[0002] 量子化学理論の発展や計算機技術の進歩によって、計算により、分子の構造、物 性や、分子内の化学結合や分子軌道、電子状態などを精度よくシミュレーションでき るようになってきた。そのような手法は一般に分子軌道法と呼ばれる。分子軌道法の 中でも、経験的パラメータに原則として依存しない非経験的分子軌道法は、近似計 算であるにせよシュレーディンガー方程式を解くために莫大な計算量を必要とするた め、これまでは、小さな分子にしか適用することができなかった。し力しながら、近年 のコンピュータの性能の著しい向上により、生体関連物質を含む比較的大きな分子 に対しても非経験的分子軌道法計算を実行できるようになり、例えば、生理活性を有 する物質の解析や探索に使用されるようになってきて!ヽる。
[0003] 非経験的分子軌道法にはいくつかの手法がある力 分子の全エネルギーを得るた めに最も使用される方法として、ハートリー 'フォック(HF)法がある。 HF法は、シユレ 一ディンガー方程式に対して 1電子近似、線形近似を行なったロータン方程式 (式 (1) )を解く方法として定式化されて ヽる。
[0004] [数 1]
FC = SCe CD
[0005] 式 (1)にお!/、て、解析対象である分子における原子軌道 (AO ; Atomic Orbital)の数 を N、分子軌道(MO ; Molecular Orbital)の数を Mとすると、 F, Sは、いずれも N X N 行列、 Cは M X N行列、 εは Μ次元ベクトルである。 Fはフォック行列と呼ばれ、式 (2) で与えられる。
[0007] ここで Dは密度行列と呼ばれ、式 (3)のように MO係数 Cによって定義されている c [0008] [数 3]
OCC
Dtu― ^ ^ 。 リ uj (3)
[0009] ここで、記号∑での" occ"は、電子が占有する分子軌道にっ 、ての和であることを表 している。 S, H, (rs I tu)は、各々、重なり積分、 Hコア積分、 2電子積分と呼ばれる 物理量であり、原子軌道 φを用いて、式 (4)一 (6)のように表される。
[0010] [数 4]
00
(4)
'OO
(rs\tu) = I φ;(1)φ3(1) ^ φ;(2) φα(2)ά3χ (6)
ノ一oo r12
[0011] ここで、記号∑における" core"は、原子核について和であることを表す。ところで式 (1) は、線型方程式の型に書かれてはいる力 フォック行列 Fが原子軌道 φに依存して 決まるため、実際には非線形方程式であって、完全に解くことはできない。そこで、こ の方程式を解くのに用いられるのが、自己無撞着場の方法(SCF(Self-consistent field)法)である。
[0012] SCF法では、よく知られているように、
[1] MO係数 Cの初期推測値を求める;
[2] 密度行列 Dを MO係数 Cから求める;
[3] 得られた密度行列 Dを使、、フォック行列 Fを求める;
[4] フォック行列 Fを対角化し、固有値 εと固有ベクトルを求める;
[5] 得られた固有ベクトルから、新しい ΜΟ係数 Cと密度行列 Dを求める; [6] 密度行列 Dが変化しなくなるまで、 [3]から [5]までをくり返す;
t 、う手順で計算が行われる。
[0013] SCF法の計算において、その費やす時間の大部分は [3]のフォック行列 Fの計算 である。その理由は、 2電子積分 (rs I tu)の計算を N4回繰り返す力 である。一度計 算した 2電子積分の結果をディスクなどのストレージ装置に保存しておく方法も考えら れるが、大規模計算、例えば Nが数万程度の場合であると、必要なディスク容量が膨 大になってしまうため、多くの場合は 2電子積分を毎回計算するダイレクト方式が取ら れる。したがって、フォック行列 Fの計算を高速ィ匕することが、 SCF法計算全体の高 速化に直結する。
[0014] この分子軌道法計算を高速に実行する方法として、例えば、特開 2000— 293494 号公報、特開 2001— 312485号公報及び特開平 9— 50428号公報に開示された方 法がある。これらの方法は、 1台のべ外ル計算機などのホスト機で行列計算などを行 ない、その一方で、フォック行列 Fの計算または 2電子積分の計算を並列計算機や計 算機クラスタに行なわせる方法である。
[0015] し力しこれらの方法では、対角化などの行列計算をいわゆるホスト機で行なうため、 ホスト機のメモリ容量を超える大きさの行列を扱うことはできない、という問題点がある 。これに対してはホスト機を複数台並列に用意するなどの対処法は存在するが、この 対処法は、高価なため、容易に行なうことはできない。近年は、汎用計算機の低価格 化や性能向上もあり、計算機クラスタで高速で廉価なシステムを作れるようになった。 汎用計算機クラスタでは、計算機 1台あたりのメモリ容量は高性能計算機に比べ小容 量である力 接続する計算機の台数を増やせばシステム全体の容量は大容量になる 。また、接続する計算機の台数を増やすことは、高速化にもつながる。ただし、汎用 計算機クラスタでは、計算機 1台ごとのメモリ容量は小さいため、大型分子の行列計 算を 1台の計算機で行なうことはできず、分子軌道数が数万軌道になると、 1つの行 列を 1台の計算機のメモリ上に確保することもできない、という問題点がある。
特許文献 1:特開 2000—293494
特許文献 2:特開 2001—312485
特許文献 3:特開平 9-050428
発明の開示
発明が解決しょうとする課題
[0016] 上述したように、計算機クラスタのような並列計算機システムを用いたとしても、大規 模かつ高速な分子軌道計算を行うためには解決すべき多くの課題が残されている。 大規模かつ高速な分子軌道計算を実現するために必要なのは、第 1に、全ての行列 要素を常に分散してメモリ領域に保存し、行列計算の際にも 1台の計算機上に集め ない計算方法の確立である。これにより、計算規模はシステム全体の総メモリ容量の みに依存し、計算機 1台のメモリ容量には依存しなくなり、また高価なホスト機を用意 する必要もなくなる。すなわち、大規模計算を行なうためには、接続される廉価な計 算機の個数を多くすればよいことになる。第 2に、上記の方法において、計算量を削 減する手法の確立である。 1台の計算機上で計算する場合には計算量を 1Z16にす る方法が知られている力 上記のように全ての行列要素を分散してメモリ領域に保存 した場合の計算量を削減する方法は知られて 、な 、。並列台数を増やせば高速ィ匕 できるが、計算量を削減できればさらなる高速ィ匕につながるので、このような手法の 確立は重要である。
[0017] そこで本発明の目的は、密度行列を分割しても、アルゴリズムの工夫によって計算 を高速に行えるようにでき、さらに、大規模計算を可能にし、従来は計算できないで いた生体高分子などの分子軌道計算を行うことができる並列計算方法を提供するこ とにある。
[0018] 本発明の別の目的は、密度行列を分割しても、アルゴリズムの工夫によって計算を 高速に行えるようにでき、さらに、大規模計算を可能にし、従来は計算できないでい た生体高分子などの分子軌道計算を行うことができる並列計算装置を提供すること にめる。
課題を解決するための手段
[0019] 本発明の並列計算方法は、複数の計算機力 なる計算機クラスタを使用する段階 と、密度行列を複数の部分密度行列に分割してそれらの複数の部分密度行列を複 数の計算機上に分散して格納する段階と、複数の部分密度行列を複数の計算機間 で順番に転送させながら、各計算機において部分密度行列に対する演算処理を実 行する段階と、を有する。
[0020] この並列計算方法は、上述した RHF法の計算に適したものである。本発明によれ ば、密度行列を複数の計算機上に分散して保存した場合でも RHF法の計算が可能 となり、計算時間を短縮することができる。本発明では、密度行列の複製を用いて積 分計算を半数に減らし、計算時間を短縮するようにしてもよい。また、密度行列とその 複製を計 4個用意し、 2電子積分における (rs I tu) ^ (tu I rs)の対称性を使うことで 、さらに計算時間を短縮するようにしてもよい。
[0021] 本発明の並列計算装置は、分子軌道法におけるハートリー 'フォック法の計算を実 行するための並列計算装置であって、複数の計算機を備える計算機クラスタを有し、 各計算機は、密度行列を分割した部分密度行列を格納する行列格納部と、計算機ク ラスタ中の他の計算機に対して部分密度行列を送受信する転送制御部と、行列格納 部に格納された部分密度行列に関する演算を実行する演算処理部と、を有し、複数 の部分密度行列を複数の計算機間で順番に転送させながら、各計算機にぉ 、て前 記部分密度行列に対する演算処理を実行する。
[0022] 本発明は、密度行列を分割したノヽートリー'フォック法の計算方法に関するものであ つて、分割された部分密度行列を順次転送することによりフォック行列の生成を可能 にし、作業領域を増やし異なる転送順序を組み合わせることにより計算量を減らし、 転送方法を一回跳ばすことにより転送量を削減したものである。本発明の方法は、計 算規模がシステム全体のメモリ容量にのみ依存する方法であるので、多数のコンビュ ータを並列接続することで大規模計算を可能にするととともに、計算時間を短縮する ことを可能にする。
図面の簡単な説明
[0023] [図 1]分散メモリ型並列計算機すなわち PCクラスタとして構成された本発明の実施の 一形態の並列計算装置の構成を示すブロック図である。
[図 2]各計算機の論理的な構成を示すブロック図である。
[図 3]本発明の並列計算方法に基づぐ密度行列を分割した場合におけるフォック行 列生成アルゴリズムを示すフローチャートである。
[図 4]ノード数 4の場合における領域 (RS) (行歹 (RS))に対する転送順序を表した図 である。
[図 5]ノード数 4の場合における領域 (TU) (行列 D(TU))に対する転送順序を表した 図である。
[図 6]ノード数 4の場合における領域 (RU) (行列 K(RU))に対する転送順序を表した 図である。
[図 7]ノード数 4の場合における領域 (TS) (行列 D(TS))に対する転送順序を表した図 である。
[図 8]図 6に示した転送順序を分力りやすく書き直した図である。
[図 9]ノード数 4のとき、 2電子積分 (RS I TU)を (m | n)と書き直し、行列であるとみ なして、 2電子積分の計算実行の有無を説明する図である。
[図 10]ノード数 9のとき、 2電子積分 (RS I TU)を (m | n)と書き直し、行列であると みなして、 2電子積分の計算実行の有無を説明する図である。
発明を実施するための最良の形態
[0024] 本実施形態では、図 1に示すような計算機クラスタすなわち分散メモリ型並列計算 機を想定している。図 1に示す本発明の実施の一形態の計算機システムは、同様の 性能を持った複数台の計算機を、通信機器によって接続したものである。従来の分 子軌道計算アルゴリズムでは、大容量のメモリを必要とするため、ホスト計算機として いわゆるスーパーコンピュータなどの高性能計算機を用意することが多い。そして、 高性能計算機は一般に高価であるため、費用面の問題となる場合がある。ところが本 発明の方法を用いれば、そのようなホスト計算機を必要としないため、費用を軽減で きる。
[0025] 図 1に示した例では、計算機システムは、複数台 (ここでは n台)の計算機 1
1一 1と、 n これら計算機 1一 1が接続するハブ 2とを備えた計算機クラスタとして構成されている
1 n
。計算機 1一 1としては、典型的には、パーソナルコンピュータ (PC)が使用されるの
1 n
で、この計算機クラスタは、 PCクラスタということになる。ここでは 1台のハブ 2を用いて 複数の計算機を接続しているが、各計算機を接続する形態はこれに限定されるもの ではなぐ例えば、リング型あるいはバス型のネットワークに接続するようなものであつ てもよい。
[0026] 図 2は、上述した計算クラスタを構成する各計算機における論理的な機能を示した
ブロック図である。本実施形態では、フォック行列の算出に必要な密度行列などの各 行列を部分行列に分割して各ノード (計算機クラスタを構成する各計算機)に分割し 、各ノードすなわち計算機ではその格納して 、る部分行列に対してあるいはその部 分行列に基づいて演算を行うとともに、ノード間でそのような部分行列を転送し、この ような演算と転送とを繰り返すことによって、最終的な結果 (例えばフォック行列)を得 るようにしている。このように各計算機が機能するために、各計算機 (ノード)は、図 2 に示すように、部分行列を格納する行列格納部 11と、行列格納部 11に格納されて V、る部分行列を別のノードに転送し、別のノードから部分行列を受け取るための転送 制御部 12と、例えば 2電子積分などの演算を行って行列格納部 11内の部分行列に 関する演算処理を実行する演算処理部 13と、演算と転送との繰り返しの制御を行う 制御部 14と、が設けられている。
[0027] 次に、本実施形態におけるこのような並列計算装置を用いた、密度行列 Dを分割し た場合のフォック行列生成アルゴリズムの手順について、図 3を用いて説明する。本 実施形態の方法では、密度行列を複数の部分密度行列に分割し、それをノード間で 転送しつつ 2電子積分を繰り返し、最後に Hコア行列との和を足し合わせることによつ て、フォック行列を生成する。図 3では、複数のノード間での部分密度行列の転送を 説明するために、 PC1と PC2の 2つのノードでの処理が並列して示されている。
[0028] まず、カウンタ nを用意して、ステップ 101において、カウンタの初期値を n=0とする 。次に、ステップ 102において、 nに 1をカ卩算し、ステップ 103において、 2電子積分を 部分的に計算し、計算した計算した 2電子積分値を用い、上述の式 (2)の計算を行う 。そして、ステップ 104において、次に必要な部分密度行列を得るために、密度行列 の転送を行う。ここでノード PC1からノード PC2に部分密度行列が転送される。ノード PC1には別のノード力 部分密度行列が送られて来るので、ノード PC1は、ノード PC 2に転送した部分密度行列の代わりに、その送られてきた部分行列を行列格納部内 に格納する。またノード PC2は、格納していた部分密度行列をさらに別のノードに転 送する。
[0029] そして、ステップ 105において、カウンタ nが、ノード数 (すなわち計算機数)に満た ないかかどうかを判定する。 nくノード数の場合には、ステップ 102— 104の処理を繰
り返すためにステップ 102に戻り、カウンタ nを 1だけインクリメントして再び 2電子積分 を行う。このとき行う計算では、転送されてきた部分密度行列に合わせ、前回とは異 なる部分を計算させる。一方、ステップ 105においてカウンタ nがノード数と等しくなる 力または超えたときは、全ノードで全ての部分密度行列に対して計算が行われたこと になるので、ステップ 106において、 Hコア行列等を足し合わせ、計算を終了する。
[0030] 本発明の並列計算方法は上述したものに限定されるものではない。転送される部 分行列や転送の態様を選択することによって、以下に説明するように、各種の実施形 態で本発明を実施することができる。以下の説明では、説明の簡単のために、ノード 数が 4であり、また、行列を部分行列に分割する際の分割数が 4であるとする。
[0031] 《例 1:巡回密度行列法》
ます、巡回密度行列法と呼ぶ実施形態について説明する。再びフォック行列 Fを式 (7)のように表す。
[0032] [数 5]
N
^rs hTS + tu (rs\tu)― - (rt\su) (7)
ム
t Σ,u=l
[0033] 式 (7)において、行列を表しているのは、 F , H , (r, s, t, u= l, · ··, N)である
。行列の分割数 4より、式 (8)に示すように、 r, s, t, uを各々 2つの領域に分割する。
[0034] 園
Rl, Sl,Tl, Ul = 1, . .. , ΛΓ/2, R2, 52,T2, C/2 = ΛΓ/2 + 1, .. . , ΛΓ, (8)
[0035] このように領域を分割すると、行列は 4つに分割される。分割された部分行列を式
(9)のように呼ぶことにする。
[0036] [数 7]
F(ll) = ¾i5i H(H) = msi D{aa) = Ότιυι \
F{21) = Fmsi H(21) = H職 D{ba) = DTWl
(9)
F{12) = Fms2 H{12) = Hms2 D{ab) = DT1U2
F{22) = ¾252 H(22) = HmS2 D{bb) = Drwi /
[0037] また、これらをまとめて表す場合には、式 (10)のように表すこととする。
[0038] [数 8]
F(RS) = { (11), (21), (12), (22)}
H{RS) = {H(11),H(21), H(12), H(22)} (10)
D{TU) = {D{aa),D{ba),D{ab),D{bb)}
[0039] この部分行列を、 4つのノード (計算機)に分配する。ノード名を Pll, P21, P12, P22として、部分行列が各ノードに式 (11)に示すように分配されるとする。このように行 列が分配された状態を初期分配状態 (S— 0)と呼ぶことにする。
[0040] [数 9]
P11 F{ll),H(ll),D(aa)
P21 F{2l),H{2l),D{ba)
(11)
12 F{12),H{12),D{ab)
P22 F{22),H{22),D{bb) ノ
[0041] さて、式 (7)における t, uについての和は、全領域について計算しなければならない 力 各ノードには密度行列 Dの一部が保持されているだけなので、このままでは計算 できない。そこで、密度行列 Dを転送し、その都度、式 (7)の一部を計算していくことに する。この方法がうまく行くための条件は、各ノードが常に行列 Dの各々異なる部分 行列を持つことである。そのような方法のひとつとして、巡回的に部分行列を転送す る方法がある。これは式 (7)を式 (12)— (14)のように計算することである。
[0042] [数 10]
P11 F{11) = H(ll) + D{aa)G{llaa) + D{bb)G{llbb) + (a6)G(llo6) + D{ba)G{llba) P21 F(2l) = H(21) + D{ba)G(2lba) + D{aa)G{2laa) + D{bb)G{2lbb) +Z>(a6)G(21a6) P12 (12) = H(12) + D{ab)G(l2ab) + D{ba)G(12ab) + D{aa)G{Uab) + D(bb)G{12bb) P22 F(22) = H(22) + D(bb)G{22bb) + D(ab)G{22ab) + D{ba)G(22ba) + D(aa)G{22aa)
■■■(12)
G{RSTU) = (RS\TU) - ^{RU\TS) (13)
R = {Rl, 2}, 5 = {51, 52}, T = {Tl,T2}, U = {U1,U2}, (14)
[0043] ここで、 1, 2は R, Sの部分領域を表し、 a, bは T, Uの部分領域を表すので、それら の領域全てについて和を取ることとし、和の記号は省略した。式 (12)は、それぞれノー ド PI 1— P22上での計算を表しており、初期分配状態 (S - 0)では第 2項まで計算で きる。第 3項以降は、密度行列 Dを転送した後でなければ計算できない。そこで、密 度行列 Dにのみ注目すると、その転送順序は図 5に示すようになつていることが分か る。ここで、図の(S-0)—(S— 3)は、転送によって変化した分配状態の名前である。
部分行列 D(aa)に注目すれば、 Pl l, P21, P12, P22の順に巡回することが分かる 。逆に、ノード P21に注目すると、 D(ba), D(aa), D(bb), D(ab)の順に部分密度行列 が巡って来ることが分かる。すなわち、全ての部分密度行列がノード P21に順次転送 されてきたことになる。全ノードについて同様なので、結局、式 (7)を計算するには、こ の順序に部分密度行列を転送すればよいことになる。この転送を実現するには、ノー ドの名前を式 (15)のように読みかえると便利である。
[0044] [数 11]
P(l) = ll, P(2) = P21, P(3) =尸 12, P(4) = P22 (15)
[0045] このようにすると、以下のようなアルゴリズムによって、上記の転送による計算は実現 できる。
[0046] [1] Fの一部分を計算する;
[2] ノード P(i)は、自らが格納して 、る部分密度行列を P(i+ 1)へ送信する; [3] ノード P(i)は、 P(i— 1)が格納していた部分密度行列を受け取る;
[4] [1]一 [3]をノードの数だけ繰り返す。
[0047] ただし、 4ノードで考えているので、ノード番号 iには周期的境界条件 i=i+4を課す 。図 5のように、ノードの数だけの回数 (ここでは 4回)転送を行うと、部分密度行列が 一周し、計算の初期分配状態 (S - 0)に戻り、計算が終了する。
[0048] 《例 2.二重巡回密度行列法》
次に、二重巡回密度行列法と呼ぶ形態について説明する。ここでは、密度行列 D の複製を用いることにより、積分計算を半数に減らし、計算時間を短縮するようにして いる。ここでの記号は、上述した巡回密度行列法と同様のものを使用する。この方法 では、式 (7)を式 (16)のように分解する。
[0049] [数 12]
[0050] ここで J, Kは、各々、クーロン積分と交換積分の和を示す行列である。これを、巡回 密度行列法の場合と同様に、領域とノードごとの計算に書き直すと、式 (17),(18)を得 る。
[0051] [数 13]
P11 F{11)
P21 (21)
P12 F{12)
P11 J(ll) = Ι?(αβ) (11|αα), J(ll) = Z)(66) (ll|6b), J(ll) = £»(β6)(11|α&), J(ll) = £)(6a)(ll|ba) P21 J(21) = D(6a)(21|6a), J(21) = Z) (oo) (21 |ao), J(21) = i3(66)(21|6 &) , J(21) = i3(a6)(21|o6) P12 J(12) = jD(ab)(12|a6), J(12) = jD(6a)(12|6a), J(12) = jD(aa)(12|ae), J(12) = JD(&6)(12|6 &) Ρ22 J(22) = D(66)(22|&6), J(22) = D(a6)(22|a6), J(22) = D(6a)(22|6e), J(22) = D(oa)(22|a )
P11 K(la) = D(al)(ll|aa), K(lb) =
K(la) = £)(61)(ll|6a) 21 Jf(2o) = D(61)(21|bo), Jf(2o) = D(ol)(21|oo), ίΓ(26) = D(61)(21|b6)
I Jf(26) = D(al)(21|ab) P12 K{lb) = D(o2)(12|o6), JT(lo) = D (& 2)(12|6a), K{la) = D(a2)(L2|o ), = D(62) (12|66) Ρ22 K{2b) = D(62)(22|66), K(2b) = D(o2)(22|o6), K{2a) = D(62)(22|6a), Jf(2o) = D(o2) (22|oe)
■■■ (18)
[0052] 二こでは、各ノードで計算される 2電子積分 (rs | tu)が、 J, Kの分配状態ごとに等し くなるように並べられている。言い換えると、 2電子積分 (rs | tu)を一回計算するだけ でよいように並べてある。そのために、 Kの部分行列も転送しなければならない。行列 J(RS), D(TU), K(RU), D(TS)の転送は、それぞれ、図 4一図 7に示すようになる。 K(RU) (R, S, T, Uは領域)の転送は、一見複雑に見える力 図 8のように書き直す と分かりやすい。すなわち、 K(RU)については、 P11と P12、 P21と P22の間でのみ 転送がおき、一方 D(TS)は、図 7より、 P11と P21、 P12と P22の間でのみ転送がおこ る。これは、ノードを式 (19)のように読みかえることで分力りやすくなる。
[0053] [数 14]
Ρίΐ, ΐ) = Ρ11, Ρ{2, 1) = P21, P(l, 2) = P12, Ρ(2, 2) = Ρ22 (19)
[0054] このようにしたとき、
[1] K(RU)の転送は、 P(i, j), (i, j = { l, 2})の iが等しいノードとの間でのみ行わ れ、また、
[2] D(TS)の転送は、 P(i, j), (i, j = { l, 2})の jが等しいノードとの間でのみ行わ れる。
[0055] K(RU)は、転送される場合と転送されな!ヽ場合とが存在するが、これは分配状態(
S— n)の状態番号 n ( = 0, 1, 2, · ··)と巡回周期 t (= 1, 2, 3, · ··)で判定できる。まず 、巡回周期 tは、ノード数を mとして、式 (20)で与えられる。
[0057] これを用いることにより、 K(RU)の転送が起きる条件は具体的には式 (21)のようになる
[0058] [数 16] i = {n mod ί) + 1 (21)
[0059] これが成り立つノード P(i, j)でのみ、次回に転送を起こさなければならない。以下に おいて、式 (21)を転送条件と呼ぶ。これら転送を実現する計算アルゴリズムは、以下 のようになる。
[0060] [1] J, Kの一部分を計算する;
[2] P(i, j)は、転送条件が成立するとき、自らが格納する K(RU)を P(i, j + 1)へ送 信する;
[3] P(i, j)は、転送条件が成立するとき、 P(i, j 1)が格納していた K(RU)を受け 取る;
[4] P(i, j)は、自らが格納する D(TS)を P(i+ 1, j)に送信する;
[5) P(i, j)は、 P(i— 1, j)が格納していた D(TS)を受け取る;
[6] P(k)は、自らが格納する D(TU)を P(k+ 1)へ送信する;
[7] P(k)は、 P(k— 1)が格納していた D(TU)を受け取る;
[8] [1]一 [7]をノードの数だけ繰り返す;
[9] フォック行列 Fを計算する。
[0061] ただし、 k=(j-Di +iである。この方法でも、全ての転送が終了したとき、分配状 max
態は初期状態 (S - 0)に戻る。
[0062] 《例 3.四重巡回密度行列法 その 1》
次に、四重巡回密度行列法と呼ぶ形態について説明する。ここでは、密度行列 Dと その複製とを合計 4個用意し、さらに、 2電子積分における (rs I tu) ^ (tu I rs)の対
称性を用いることで、二重巡回密度行列法よりもさらに積分計算を減らし、計算時間 の短縮を図っている。ここでの記号は、上述した巡回密度行列法と同様のものを使用 する。この方法では、式 (7)を式 (22)— (24)のように分解する。
[0063] [数 17]
m = {s - 1)N + r, n = {u - 1)N + 1, (23)
d(m,n) = I for m n, 1/2 f or m = n (24)
[0064] これは、 r, sによって一意に決まる数 mと、 t, uによって一意に決定される数 nを用 いて、クーロン積分 J及び交換積分 Kを 2つの行列に分けて計算することを表している o 2電子積分 (rs I tu)は、 m, nを用いて (m | n)と表される。 2電子積分は (rs | tu) = (tu I rs)という対称性を持つので、条件 m≤nより、 (m | n)を行列と考えると、 J1, K1は (m I n)のその行列の下三角部分、 J2, K2は (m | n)の上三角部分を計算し ていることになる。条件 m≤nをひとまず無視して、式 (25)ズ26)に示すように、式 (22)を 領域とノードごとの計算に書き直す。
[0065] [数 18]
ΛΓΙ(ΙΙ) + ϋΓ2(11) 、
Kl(21) + K2(21) Kl(12) + K2(1 ) (25)
Kl(22) + ϋΓ2(22) ノ ll Jl(ll)
P21 Jl(21)
P12 Jl(12)
Pll J2(aa) = 0(ll)(ll|afl)/2, J2(&&) = = ΰ(11)(11|α&), J2(&a) = = Ζ)(11)(11|&α) P21 J2(6 ) -D(21)(21|6a)/2, J2(aa) - D(21)(21|oo), - £)(21)(21|&&), J2(a&) - : )(21)(21|a&) P12 J2(a&) =D(12)(12|a6)/2, J2(6a) = D(12)(12|6a), = D(12)(12|a ), J2(66) = (12)(12|W>) P22 J2(66) = = (22)(22|a )
Pll iCl(la)
P21 (2t¾5
P12 Kl(lb)
P22 Kl(2b)
Pll iT2(al)
P21 K2(bl)
12 K2{a2)
P22 K2(b2)
[0066] この計算を実現する各行列の転送順序としては、二重巡回密度行列法と同じもの を使用する。行歹 1(RS), D(RS)については図 4、行列 D(TU), J2(TU)については 図 5、行列 K1(RU), D(RU)については図 6、行列 D(TS), K2(TS)については図 7に したがった転送を行う。すなわち領域 (RS), (TU), (RU), (TS)によって転送方法が 異なる。次に、条件を領域ごとに適用できるように、数 m, nを式 (27)のように書き直す
[0067] [数 19]
m = (J-l)imBx+i {i = {1,2}), n = {l-l)krmx+k (k,l = {a, b}) (27)
[0068] ここで、 i, jは領域 R, Sの番号 {1, 2}を表し、 k, 1は領域 T, Uの番号 {a, b}を表す。
ただし、 a=l, b = 2と換算する。この m, nを使用して、行歹 を Jl, J2に、行列 Κを Κ 1, Κ2に分けて計算する条件は、以下のようにいくつか考えられる。
[0069] [数 20]
条件 (C_1)
m> n
条件 (C_2)
m < n
条件 (C_3)
m = n, m +n = odd for m <n,, m + n = even for m > n 条件 (C_4)
m = n, m+n = odd for m> n,, m + n = even for m <
[0070] 例えば、条件 (c 1), (c—3)を適用した場合は、 Jl, J2, Kl, K2の計算は次のよう になる。
[0071] 条件 (c 1)を適用した場合、 Jl, J2, Kl, Κ2の計算は式 (28)で表される。
[0072] [数 21]
[0073] 条件 (c 3)を適用した場合、 Jl, J2, Kl, K2の計算は式 (29)で表される。
[0074] [数 22]
ll Jl(ll)
P21 Jl(21)
P12 Jl(12)
Pll J2(a ) = D(ll)(ll|a )/2} J2( & 6) = D(ll)(ll\bb) J2(&a) = Z)(ll) (ll|6a)
P21 J2{ba) = Χ)(21) (21|&ο)/2} J2(a&) = D{2l) (2l\ab)
P12 J2(ab) = Γ>(12) (12|α&)/2} J {aa) = D{12) {12\a )J J2{bb) = Z)(12) (12|66)
• · · (29)
[0075] 条件 (c 3) , (c-4)の方が、各ノードの計算量が比較的平均化される。結局、計算 アルゴリズムは以下のようになる。
[0076] [1] 条件(c 1)一(c 4)の何れかを満たすとき、 Jl, J2, Kl, K2の一部分を計 算する;
[2] P(i, j)は、転送条件が成立するとき、自らが格納する Kl(RU), D(RU)を P(i, j + 1)へ送信する;
[3] P(i, j)は、転送条件が成立するとき、 P(i, j 1)が格納していた K1(RU), D(R U)を受け取る;
[4] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 1, j)に送信する; [5] P(i, j)は、 P(i— 1, j)が格納していた K2(TS), D(TS)を受け取る;
[6] P(k)は、自らが格納する J2(TU), D(TU)を P(k+ 1)へ送信する;
[7] P(k)は、 P(k— 1)が格納してい^ J2(TU), D(TU)を受け取る;
[8] [1]一 [7]をノードの数だけ繰り返す;
[9] J, Kを計算する;
[10] Fを計算する。
[0077] ただし、ノード番号は、式 (30)のように読みかえて 、る。
[0079] 《例 4.四重巡回密度行列法 その 2》
次に、四重巡回密度行列法と呼ぶ形態の別の例について説明する。
[0080] 上述した四重巡回密度行列法では、ノード数が奇数個の場合に限り、条件 (c 3) ,
(c 4)下の計算を別の方法で行なうことができる。ノード数が 4, 9で条件 (c— 3)を満 たす場合について、領域ごとの 2電子積分 (rs I tu) = (m | n)を、図 9、図 10のよう に図示する。白四角の部分は、計算する 2電子積分 (m I n)を表し、黒四角の部分 は、条件 (c— 3)によって計算しない (m I n)を表す。また、各ノード、 P (i, j) (i, j = { l, 2}または { 1, 2, 3})については、横一列を計算する。計算の順序は、対角成分 (m I m)から計算が始まり、転送により(m I m 1)に移動する。図 9、図 10の四角形内 部の数字が個の計算順序を表している。ノード数 9の場合は、 2, 4, 6, 8に対応する 四角形部分は全ノードで黒四角形であり、 2電子積分を計算しないことを示している。 ノード数が奇数個の場合は、転送を繰り返すと、黒四角形部分と白四角形部分とを 交互に迪ることになる。ノード数 4の場合は、どの番号の四角形部分においても必ず 一つは白四角形があり、ノード数が偶数個の場合にも一般ィ匕できる。したがって、ノ ード数が奇数個の場合には、 1, 3, 5, 7, 9の四角形部分に対応して計算しさえすれ ばよぐ転送を一回分跳ばして行なうことができる。ただし、最後の一回だけは一回分 のみ転送する。この場合の計算アルゴリズムは次のようになる。
[0081] [1] Jl, J2, Kl, K2の一部分を計算する;
[2] P(i, j)は、転送条件が成立するとき、自らが格納する Kl(RU), D(RU)を P(i, j + 1)へ送信する;
[3] P(i, j)は、転送条件が成立するとき、 P(i, j 1)が格納していた K1(RU), D(R U)を受け取る;
[4] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 2, j)に送信する; [5] P(i, j)は、 P(i— 2, j)が格納していた K2(TS), D(TS)を受け取る;
[6] P(k)は、自らが格納する J2(TU), D(TU)を P(k + 2)へ送信する;
[7] P(k)は、 P(k— 2)が格納してい^ J2(TU), D(TU)を受け取る;
[8] [ 1]一 [7]を (ノード数 Z2)に相当する回数繰り返す;
[9] Jl, J2, Kl, K2の一部分を計算する;
[10] P(i, j)は、転送条件が成立するとき、自らが格納する Kl(RU), D(RU)を P(i , j + 1)へ送信する;
[11] P(i, j)は、転送条件が成立するとき、 P(i, j— 1)が格納していた K1(RU), D( RU)を受け取る;
[12] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 1, j)に送信する; [13] P(i, j)は、 P(i— 1, j)が格納していた K2(TS), D(TS)を受け取る;
[14] P(k)は、自らが格納する J2(TU), D(TU)を P(k+ 1)へ送信する;
[15] P(k)は、 P(k-l)が格納してい 2(TU), D(TU)を受け取る;
[16] J, Kを計算する;
[17] Fを計算する。
[0082] 条件 (c 4)を満たす場合は、転送順序を逆にするか、最後の一回分のみの転送を 最初に行なえばよい。すなわち、上記のアルゴリズムを [9] , · ··, [15] , [1] , · ··, [8] , [16] , [17]と行なえばよい。
[0083] 《例 5.—般化》
上述した各種の巡回密度行列法は、ノード数 N、分割数 M ( >N)の場合でも計算 することができるが、各行列の分割方法は揃えておく必要がある。すなわち、領域尺と Tの個数が一致し、領域 Sと Uの個数が一致していなければならない。「四重巡回密 度行列法その 2」では、一般ィ匕すると次のようなアルゴリズムになる。
[0084] [1] 条件 (c 3)を満たすとき、 Jl, J2, Kl, K2の一部分を計算する;
[2] P(i, j)は、転送条件が成立するとき、自らが格納する Kl(RU), D(RU)を P(i, j + 1)へ送信する;
[3] P(i, j)は、転送条件が成立するとき、 P(i, j 1)が格納していた K1(RU), D(R U)を受け取る;
[4] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 2, j)に送信する;
[5] P(i, j)は、 P(i— 2, j)が格納していた K2(TS), D(TS)を受け取る;
[6] P(k)は、自らが格納する J2(TU), D(TU)を P(k + 2)へ送信する;
[7] P(k)は、 P(k— 2)が格納してい^ J2(TU), D(TU)を受け取る;
[8] Jl, J2, Kl, K2の一部分を計算する;
[9] P (i, j)は、転送条件が成立するとき、自らが格納する K1(RU), D(RU)を P(i , j + 1)へ送信する;
[10] P(i, j)は、転送条件が成立するとき、 P(i, j-1)が格納していた K1(RU), D( RU)を受け取る;
[11] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 2, j)に送信する;
[12] P(i, j)は、 P(i-2, j)が格納していた K2(TS), D(TS)を受け取る; [13] P(k)は、自らが格納する J2(TU), D(TU)を P(k + 2)へ送信する; [14] P(k)は、 P(k— 2)が格納してい 2(TU), D(TU)を受け取る;
[15] [8]— [14]を { (ノード数 Z2)—l }に相当する回数繰り返す;
[16] Jl, J2, Kl, K2の一部分を計算する;
[17] P(i, j)は、転送条件が成立するとき、自らが格納する Kl(RU), D(RU)を P(i , j + 1)へ送信する;
[18] P(i, j)は、転送条件が成立するとき、 P(i, j 1)が格納する K1(RU), D(RU) を受け取る;
[19] P(i, j)は、自らが格納する K2(TS), D(TS)を P(i+ 1, j)に送信する;
[20] P(i, j)は、 P(i— 1, j)が格納する K2(TS), D(TS)を受け取る;
[21] P(k)は、自らが格納する J2(TU), D(TU)を P(k+ 1)へ送信する; [22] P(k)は、 P(k— 1)が格納してい 2(TU), D(TU)を受け取る;
[23] J, Kを計算する;
[24] Fを計算する。
[0085] 条件 (c 4)を満たす場合は、転送を逆に巡回させる力、最初の転送を一回分にす ればよい。
[0086] 以上、本発明の好ましい実施形態を説明したが、本発明は、計算機クラスタでの実 現を前提とするものである。したがって、計算機クラスタを構成する各計算機は、上述
した各ノードでの処理を実行するものでなくてはならない。各計算機は、ノードとして の処理を実現するための計算機プログラムを読み込み、そのプログラムを実行するこ とによって、上述した各処理を実行するようになる。そのようなプログラムは、磁気テー プゃ CD— ROMなどの記録媒体によって、あるいはネットワークを介して、計算機に 読み込まれる。
[0087] 具体的にはそのプログラムは、複数のノードから構成される計算機クラスタにおける 各ノードの計算機を、密度行列を分割した部分密度行列を格納する行列格納部、計 算機クラスタの他のノードに対して前記部分密度行列を送受信する転送制御部、行 列格納部に格納された部分密度行列に関する演算を実行する演算処理部、として 機能させ、それによつて、複数の部分密度行列が複数のノード間で順番に転送され ながら、各ノードにおいて部分密度行列に対する演算処理が実行されるようにする。
[0088] さらには本発明の範疇には、上述したプログラム力 なるプログラムプロダクト、この プログラムを格納した機械可読記録媒体、このプログラムを伝送する伝送媒体も含ま れる。