WO2005029352A1

WO2005029352A1 - 並列計算方法及び装置

Info

Publication number: WO2005029352A1
Application number: PCT/JP2004/013734
Authority: WO
Inventors: Toshikazu Takada; Jun-Ichi Yamamoto; Kazuto Nakata
Original assignee: Nec Corporation; Nec Soft, Ltd.
Priority date: 2003-09-22
Filing date: 2004-09-21
Publication date: 2005-03-31
Also published as: US20060271301A1; US7885796B2; JP4612546B2; JPWO2005029352A1

Description

明細書

並列計算方法及び装置

技術分野

[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) で与えられる。

[0006] [数 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)

^F— 00

'OO

(rs\tu) = I φ;(1)φ₃(1) ^ φ;(2) φ_α(2)ά³χ (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)の計算を N⁴回繰り返す力である。一度計算した 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 h_TS + 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) = F_msi H(21) = H職 D{ba) = D_TWl

(9)

F{12) = F_ms2 H{12) = H_ms2 D{ab) = D_T1U2

F{22) = ¾2₅₂ H(22) = H_mS2 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)

Ρ22 (22)

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)で与えられる。

[0056] [数 15]

[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)

P22 Jl(22)

= =

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)

P22 Jl(22)

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)のように読みかえて、る。

[0078] [数 23]

[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の一部分を計算する；

[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の一部分を計算する；

[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] さらには本発明の範疇には、上述したプログラム力なるプログラムプロダクト、このプログラムを格納した機械可読記録媒体、このプログラムを伝送する伝送媒体も含まれる。

Claims

請求の範囲

[1] 分子軌道法におけるハートリー 'フォック法の計算を実行するための並列計算方法であって、

複数の計算機からなる計算機クラスタを使用する段階と、

密度行列を複数の部分密度行列に分割して該複数の部分密度行列を前記複数の計算機上に分散して格納する段階と、

前記複数の部分密度行列を前記複数の計算機間で順番に転送させながら、前記各計算機において前記部分密度行列に対する演算処理を実行する段階と、を有する並列計算方法。

[2] 前記密度行列を複製を用い、前記密度行列及び前記複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させることにより、積分計算を減らす、請求項 1 に記載の方法。

[3] 前記密度行列と前記密度行列の複製とを計 4個使用し、前記密度行列及び前記複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させ、さらに、 2電子積分における (rs I tu) <-Xtu I rs)の対称性を使うことにより、積分計算を減らす、請求項 1に記載の方法。

[4] 前記各計算機において部分的に 2電子積分を行い、該 2電子積分の結果に基づいて、格納している部分密度行列を更新する段階を有する、請求項 1に記載の方法

[5] 前記密度行列を複製を用い、前記密度行列及び前記複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させることにより、積分計算を減らす、請求項 4 に記載の方法。

[6] 前記密度行列と前記密度行列の複製とを計 4個使用し、前記密度行列及び前記複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させ、さらに、 2電子積分における (rs I tu) <-Xtu I rs)の対称性を使うことにより、積分計算を減らす、請求項 4に記載の方法。

[7] 分子軌道法におけるハートリー 'フォック法の計算を実行するための並列計算装置であって複数の計算機を備える計算機クラスタを有し、

前記各計算機は、密度行列を分割した部分密度行列を格納する行列格納部と、前記計算機クラスタ中の他の計算機に対して前記部分密度行列を送受信する転送制御部と、前記行列格納部に格納された部分密度行列に関する演算を実行する演算処理部と、を有し、

前記複数の部分密度行列を前記複数の計算機間で順番に転送させながら、前記各計算機において前記部分密度行列に対する演算処理を実行する、並列計算装置

[8] 前記密度行列及び前記密度行列の複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させることにより、積分計算を減らす、請求項 7に記載の装置。

[9] 前記密度行列と前記密度行列の複製との計 4個をそれぞれ部分密度行列に分割して前記各計算機間で転送させ、さらに、 2電子積分における (rs I tu) <-Xtu I rs)の対称性を使うことにより、積分計算を減らす、請求項 7に記載の装置。

[10] 前記演算処理部は部分的に 2電子積分を行い、該 2電子積分の結果に基づいて前記行列格納部内の部分密度行列が更新される、請求項 7に記載の装置。

[11] 前記密度行列及び前記密度行列の複製をそれぞれ部分密度行列に分割して前記各計算機間で転送させることにより、積分計算を減らす、請求項 10に記載の装置。

[12] 前記密度行列と前記密度行列の複製との計 4個をそれぞれ部分密度行列に分割して前記各計算機間で転送させ、さらに、 2電子積分における (rs I tu) <-Xtu I rs)の対称性を使うことにより、積分計算を減らす、請求項 10に記載の装置。

[13] 複数のノードから構成される計算機クラスタにおける各ノードの計算機を、

密度行列を分割した部分密度行列を格納する行列格納部、前記計算機クラスタの他のノードに対して前記部分密度行列を送受信する転送制御部、前記行列格納部に格納された部分密度行列に関する演算を実行する演算処理部、として機能させ、それによつて、前記複数の部分密度行列が前記複数のノード間で順番に転送されながら、前記各ノードにおいて前記部分密度行列に対する演算処理が実行されるようにする、プログラム。

[14] 計算機が読み取り可能な記録媒体であって、請求項 13に記載のプログラムを格納した記録媒体。