電子相関計算用積分変換モジュールの開発

望月 祐志, 長嶋 雲兵


Return

1 はじめに

分子軌道計算において単一電子配置のHartree-Fock(HF)近似を超えて電子相関を導入するには、単一配置の枠組みのままで密度汎関数理論(DFT: density functional theory)によって相関エネルギー項を導入するか、配置間相互作用(CI: configuration interaction)や摂動論(PT: perturbation theory)、あるいは結合クラスター展開(CC: coupled cluster)のように励起電子配置を含めた多配置の線形結合から相関エネルギーを求めるか、二つに大別される[1 - 4]。
分子軌道を展開する基底関数の総数をNとすれば計算量の形式的なオーダーは、よく知られているように、DFTではHFと同じく4乗、2次のPT(MP2)では5乗、CIやCCでは6乗であり、計算コスト的にはDFTが最も有利である。しかし、電子構造が複雑な系や弱い結合状態の反応中間体などではDFTの精度が低下する場合がある。PTは、弱い相互作用は記述出来るが、3d遷移元素などの局在電子の相互作用が大きい系では十分ではないこともあり、汎用性と信頼性の観点からはCIの方が本質的には好ましい。CI計算では、励起配置の生成が通常は1,2電子励起に限られるため、相関電子数が増えると記述が次第に悪化する"大きさ矛盾性"の欠陥があるが、結合電子対近似(CEPA: coupled electron pair approximation)的な考え方に基づき、元のCIコードの修正を最低限に留めつつ、多電子励起の寄与を取り込んで是正する近似法が提案されているので、中型以上の分子に対しても計算コスト的な問題を上手く解決していけるのであれば、CIの魅力は再評価されていくと期待される。
本稿では、続く2章でCI法の概要をまとめた後、3章で望月が進めているLCIプログラムの開発についてふれ、4章ではLCIの中で積分変換に関するモジュール群について個別に機能とアルゴリズムを記し、Linux上での評価について報告していく。

2 CI計算の概略

準備として、以下にCI法の概略をまとめる。文献としては教科書の[1 - 4]、CIに特化したレビューとして[5]を参考にされたい。
CIでは、分子の基底状態を表わすHF配置を参照して、占有軌道から仮想軌道へと電子を励起させた励起配置群を式(1)のように線形結合

して多電子波動関数を記述する。係数のセット{TI}は、ハミルトニアン行列(2)

を対角化する固有ベクトル、すなわちCIベクトルとして得られ、固有値がエネルギーとなる。ハミルトニアン行列は、非零率が高々数%程度で非常に疎ではあるが、長さのオーダーは容易に105となり、107に達することもある。
式(1)中の個々の配置は、系のスピン状態を満たすように励起Slater行列式群を束ねて作られる配置関数(CSF: configuration state function)によって与えられる。例えば、1電子励起によって占有軌道と仮想軌道に1つづつ開殻が生じる配置で1重項状態を満たすには、開殻部分でabbaの2つのスピン並びを持つSlater行列式をまとめて1つのCSFとする。後述するハミルトニアンの行列要素の評価は、CSFを使う方が複雑になる。そのため、軌道と電子数を決めて可能なあらゆる励起とスピン結合を含めるFull-CI計算では、簡単さを最優先してSlater行列式があらわに使われることが多い。
HF配置からの1,2電子励起CI(CISD: CI with singles and doubles)では、最低固有値のエネルギーとHFエネルギーとの差が電子相関のエネルギーに対応する。他方、1電子励起(CIS)に限って多数の根を解き、分子の励起エネルギースペクトルを見積ることもGaussian[6]ではよくなされる。
目的とする状態の記述精度を上げたり、オゾン等のジラジカルの記述や化学結合の解裂を扱う際には、参照配置を複数とする必要(MRCISD: multi-reference CISD)があるが、この場合、入力する分子軌道にはHF法ではなくCASSCF(complete active space self-consistent-field)法で求めたものを使う。CASSCFは、活性空間と称する小数の軌道群の中でFull-CI展開を行いつつ、占有軌道全体を最適化する。なお、多参照のPTは近年、開発と応用が進んできているが、多参照CCについては、未だ開発途上にあるというのが現状である。
CIでは展開長が104以上となること、またエネルギーの低い方の根だけを求めればいいことから、ハウスホルダー法等の直接対角化ではなく、試行的なCIベクトルを逐次的に改良する反復解法が用いられる。反復法では、試行ベクトルとハミルトニアン行列の積である、いわゆるsベクトル

が繰り返し計算される。改良は、残差を分子に、対角要素とエネルギーを分母に持つ補正ベクトル

を規格化を行いつつCIベクトルに繰り返し加えて行う。反復により、式(4)の分母が閾値の範囲で零となれば、エネルギーとCIベクトルが収束したと見なす。多根を解く際は、Davidson法などが用いられる。
式(2)のハミルトニアン行列の要素は、基底関数の添字から分子軌道の添字に変換

された1電子積分hijと2電子積分gij,klに対し、CSF間の軌道占有のパターンとスピン結合の組み合わせによって決まる結合係数を掛けた積の和

で書ける。結合係数のセット{AIJij, BIJij,kl}は、ハミルトニアン行列のエネルギー表現とも呼ばれる。
式(7)の和の添字は、ハミルトニアンが2電子の相互作用までしか含まないことから、配置の占有パターンの差が3電子以上では寄与が無くなるため、実際にはフルに走る必要はない。また、積分の添字に関する交換等値性

を用いれば、式(7)の非零の和もコンパクトな正準形の添字の関係

を使ってユニークな寄与だけを取れる。ここで、2電子積分に関する節約は1/8である。
正準形の和としても、分子軌道の添字を4つ持つ2電子積分gij,klの式(7)での処理はCI計算のコストを支配する。CISDの場合、CSF添字のIJの各々を指定するのに最大で4つ(占有軌道で2、仮想軌道で2)づつ励起に係わる軌道の添字が現われるので、合わせて8つとなるが、HIJで2電子積分が寄与を持つためには互いに2つは共通でなくてはならない。そのため2つ減り、結果としてコストのオーダーは基底数Nに対して6乗となる。
gij,klのリストを得る変換処理にも注意が要る。変換の式(6)には基底の添字が4つ、軌道の添字が4つ現われており、そのまま単純に8重のループで回すことはコスト的に現実的ではない。そこで、5乗オーダーのループ処理

で1/4づつの変換を続けて行う。基本的なアルゴリズムは幾つか知られてはいるが、プログラムとしての実装にあたっては改良する余地が残っている(本論文も、まさに2電子積分の変換に関するものである)。
ここで、式(3)のsベクトルの計算に話を戻す。式(7)から非零のハミルトニアン行列要素を分子積分と結合係数から一度あらわにつくってファイルに書き出し、sベクトル計算の際に繰り返し読み出して処理するCI計算の仕方をConventional-CIと呼ぶ。ここで、非零率を別とすれば、ハミルトニアン行列は本質的にNの6乗オーダーのデータ量であることに注意して欲しい。ファイルへの収納には要素の値とIJの添字で合わせて1つあたり16バイトが必要なので、非零率1%の場合、展開長が5×105の場合で20GBもの容量になる。そのためConventional-CIでは、相関エネルギーに重要な寄与を与える励起CSF群をPT計算によって予め選択して限定的なCIを行う。実際に扱える対角化の規模としては、現在でも106程度が限界となっている。フルの対角化によって得るべき"真のエネルギー"を求めるのであれば、選択の閾値を変えながら計算を繰り返して外挿によって推定する他ない。また、選択によって相関電子対の扱いが均等でなくなり得るため、CPF(coupled pair functional)[7]やACPF(averaged CPF)[8]のようにCEPA流に"大きさ矛盾性"を是正したCIがやり難いという弱みがConventional-CIにはある。
Conventional-CIに対し、ハミルトニアン行列要素をあらわにはつくらないで、個々の分子積分の寄与を直接sベクトルの中に加算

していく仕方が、Roosによって提案されたDirect-CIである。Direct-CIでファイルに保存すべき主なデータは、Nの4乗オーダーの2電子積分、結合係数BIJij,kl、ならびに式(4)の補正ベクトルの計算で分母に必要なハミルトニアン行列の対角要素であり、要求ファイルの容量はConventional-CIに比してオーダー違いで少ない。実のところ、Direct-CIでは107規模の大きなCI展開をそのまま扱うことが可能であり、今日では汎用CI計算の標準解法となっている。また、CPF/ACPF計算も容易である。
Direct-CIのプログラムは、2電子積分の処理に係る部分がConventional-CIに比してずっと複雑になる。その複雑さは「2電子積分を占有軌道の添字と仮想軌道の添字の組合せと並び方に応じて分類し、各々の積分タイプに応じて随伴するIJをグループ化してまとめ、sベクトルへの寄与を効率よく処理するループ群を個別に場合分けしてつくらなくてはいけない」ことにある。また、積分に乗じるべき結合係数の扱いにも工夫が要る。実装としては様々なやり方があり、MOLPROやCOLUMBUSなどのプログラムとして利用出来るようになっている[9 - 14]。しかし残念ながら、著者らが知る限り、国産の汎用Direct-CIのプログラムは未だ開発されてはいない。

3 LCIプログラム

LCIは、望月が提案した代表化エネルギー表現法[15]に基づき、自身で開発しているDirect-CIプログラムシステムの名称である。語頭のLには、{large, liberal,...}などが込められている。
LCIの基本的な考え方は、
であるが、「自在に手を入れられる国産のDirect-CIプログラムを確保する」ことも目的としては重要である。LCIの開発最終段階では、MOLPRO等の他のプログラムが得手としない中大型の分子の様々なCI問題を6-31G**に代表されるSVP(split-valence plus polarization)程度の基底関数で扱えることを主に想定しており、具体的にはバイオテクノロジー関係で小ペプチドやポルフィリン等の生体機能分子、あるいはナノテクノロジー関係で遷移金属の複核錯体や表面吸着系を考えている。
基底関数による積分の生成とHF/CASSCFによる軌道の最適化をGAMESS等の外部の分子軌道コードに拠るLCIシステムでは、式(5)や(6)の積分変換の処理は入力インターフェース部分に相当する。入力データとして書けば{hpq , gpq,rs , cpi}を受け取ることになる。積分を変換した後は、ハミルトニアン行列のエネルギー表現の生成、Direct-CI法によるCIベクトルの求解、その後で得られる1電子密度行列

とそれを対角化して決まる自然軌道の生成といった一連の処理はLCIシステムの内部で完結する。HF軌道の占有数は{0, 1, 2}の整数だが、自然軌道では占有軌道空間からの励起によって仮想軌道空間に電子が若干量移るために0-2の間の小数値となる。自然軌道は系の電子状態の特徴を捕らえるのに便利であるし、自然軌道を入力とした別のCI計算を行うことも可能である。従って、LCIの出力インターフェースとしては、自然軌道を他の分子軌道プログラム、あるいは可視化プログラムに渡すことが考えられる。
LCIの開発は現時点では第一段階にあり、モジュールとして開発が済んでいるのは、次章で述べる積分変換に関するステップである。生体分子や錯体では、多くの場合に対称性はC1、すなわち対称性無であることから、コーディングに際して対称性は課していない。基底関数積分と軌道係数を提供してもらうインターフェース相手としては、オープンソースの代表的なコードであるGAMESS[12]となっている。もちろん他のプログラムも可能であり、同じくオープンソースではDALTON[14]が有望な相手であるし、国内では分割分子軌道(FMO: fragmented MO)によりペプチド類の高速HF計算を可能としているABINIT-MP[16]に対して電子相関計算の機能を提供する形での連携もあり得る。

4 積分変換に係わるモジュール

GAMESSから供給される基底関数による積分と分子軌道を使って行う積分変換は、実際には3つのモジュールが関与している。以下の3節で個々の機能とアルゴリズムについて述べ、最後の1節で評価を行う。

4. 1 1電子積分用IFSORT

名称のIFSORTは、Integral-Fock-Sortから取ったものであり、モジュールの機能としては、Fock型の1電子積分の生成、ならびに2電子積分の変換前にソーティングを行うようになっている。これらについて順に記す。
LCIでは文献[15]にあるように、仮想的な閉殻電子配置を"原点"として、そこからの差分によってあらゆるCSFセットを記述する粒子-空孔の定式化によってエネルギー表現を計算し、これによって閾値性に因らず2電子積分の処理量を根本的に減らす工夫をしている。このため、1電子積分については、hijに2電子積分の寄与を予め含めた閉殻型のFock演算子の分子軌道表現

を使う必要がある。ここで、和の添字klは全ての閉殻軌道について走る。例えばCSF間で占有軌道mと仮想軌道aが違う場合のハミルトニアン行列の要素は次式

となるので、既に埋め込まれている2電子積分の和を扱う必要はない。
式(19)でfijをつくる際も、分子軌道添字に変換された2電子積分をあらわに扱うことは得策ではない。そこで、普通のHF計算と同様に、基底関数の添字で"原点"に関するFock行列

を先ず計算し、軌道添字に変換

を行う。
式(21)中で、密度行列(式(22)で定義されている)に掛け合わされる括弧内の2電子積分部分は、いわゆるSuper Matrixに相当する。しかし、Fock行列計算のためにだけSuper Matrixを改めてつくるのはコストが大きい。GAMESSでは閾値以上の個々の積分を正準番地のラベルを付け、その寄与を直接扱っており、LCIでも同じアルゴリズムを使うことにした。つまり、gpq,rsについて次の6通りの寄与を計算

する。実装では、正準添字の同値性(式(10)〜(12)で添字をpqrsに読み替える)を重み因子で適宜考慮しつつ、読み込んだ積分バッファ長で回る単一ループで処理を行う。また、下述する積分のソーティングも同時に行うので、GAMESSの積分ファイルの読み込みは一度だけで済ませている。
"原点"の仮想閉殻配置の電子エネルギーは次式

で書ける。閉殻HF配置を参照するCISDの場合、このエネルギーはHFの電子エネルギーそのものになる。また、入力された軌道が正準軌道の場合、Fock行列の対角要素fiiは軌道エネルギーeiに対応する。開殻軌道が1つの2重項HFを参照する場合、電子数が1つ少ない閉殻を"原点"とし、HF配置の参照CSFは「開殻軌道に関するaスピンの粒子生成演算子を作用させる」ものとして記述され、"原点"に対する追加のエネルギーを持つことになる。軌道xyの2つの開殻を持つ3重項HF配置を参照する場合は、エネルギーの低い方の軌道xに閉殻を作って"原点"とすると、HFのエネルギーは次式

となる。これは、「"原点"の配置に対し、軌道xbスピンの空孔を生成し、軌道yaスピンの粒子を生成する」ことで得られる。同様に、あらゆる励起CSF群の行列要素は粒子と空孔の生成演算子の積によって評価され、対応する分子積分{fij, gij,kl}と結合係数{AIJij, BIJij,kl}の積和で与えられる。こうしたエネルギー表現の詳細は、第二量子化演算子の代数について十分な準備が必要になるので、ここではこれ以上は踏み込まない。
式(13)〜(16)の2電子積分変換は、次節で述べる山本-長嶋のアルゴリズム[17]に基づいたITRANS(Integral-Transformationの意)モジュールで行う。このアルゴリズムでは、前半2つの変換pq®ijのために予め「ある1つのrs対について閾値以上のgpq,rsを持つ全てのpq対を集めておく」ソート処理が必要である。IFSORTでは、前述のFock型1電子積分fijの計算と同時にソートしている。ここで、pqの並びについてはソートする必要がないのでロジックは直截に実装出来るのだが、注意が要るのは正準添字でpqrsの場合だけなく、pqrsの場合も考慮、すなわちpqrsを置き換えた処理も含めなくてはいけないことである。これにより、正準性による節約が1/8から1/4に落ちるため、ソートされた2電子積分のファイルは元のGAMESSのファイルに比して2倍程の大きさになる。
メモリ上でソートする、すなわちインコアのソーティングの場合、第一次元をスタックしていくpq、第二次元をソーティング変数のrsとするW(N(N+1)/2, N(N+1)/2)を積分値のための配列、ならびにpqの添字値を収納する配列IW(2, N(N+1)/2, N(N+1)/2)の2つを確保出来ればよく、rsについて順に順編成ファイルに書き出せばソーティング処理は完了する。インコアで処理できない場合は、配列Wの第一次元(pqの添字値に関する配列IWでは第二次元)を適当なIOバッファ長Lで取り、GAMESSの2電子積分ファイルを読み込んで処理しながら、スタックがLになる毎に作業用の直接編成ファイルにレコードとして書き出す。この時、吉嶺の連鎖[18]に従って"1つ前のレコードの番地"を同時に書いておく。元の積分ファイルの読み込みが終わったら、スタックに残っているデータも書き出してしまう(フラッシュ操作)。
順編成ファイルへのソート済み積分の書き出しでは、直接編成ファイルを「あるrsについての最後のレコードから逆向きに次々に読んでいく」ことになる。バッファ長Lは、IFSORTの既定値では820となっている。Lを大きくすればIO頻度は減るが、その分だけメモリ要求は大きくなる。セミコア的な扱いとして、「あるrsまではインコアで、残りは直接ファイルを使ってソートする」ことも考えられるが、元の基底関数添字の2電子積分のファイルを複数回読まなくてはならないこと、またメモリ利用の最適化を含めたロジックが煩雑になることから、現時点ではシンプルさを優先して実装していない。基底関数の総数と利用出来るメモリ量に応じて、インコア式とファイル式をプログラムが自動的に選択する。CPU時間はどちらの処理も大きくはないが、直接編成ファイルでのソーティングの場合、マシン環境に因るが経過時間はかかる。

4. 2 2電子積分用ITRANS

LCIで採用した山本-長島の2電子積分変換アルゴリズム[17]には以下の特徴がある。
要は、積分のsparsenessを活かして無駄な演算を避けられる点が一番の魅力で、中大型分子、特に分子形状が広がった生体分子の扱いに有利である。実のところ、山本らが今から十五年も前に、基底総数N=232の鉄ポリフィリン錯体酵素モデルのCASSCF計算を当時の最新鋭ベクトル機日立S-810上で実現出来た[19]のも、この変換アルゴリズムをJASON2プログラムに実装したからこそであった。JASON2では、当時の計算機環境の中で最適動作させるために、変換処理の核心部分に随伴してIO処理やメモリ割当に相当の工夫がこらされていたが、今日では普及価格のパーソナルコンピュータでさえも数百MBのメモリは容易に利用出来るので、LCIの変換モジュールITRANSでは、こうした付随処理を省き、元のアルゴリズムをシンプルにコーディングすることを先ず考えた。その際、最深部には、世界標準の基本線形代数ライブラリであるBLAS[20]の中で、データ移動量が1次であるサブルーチンDAXPYとDDOTを導入することにした。BLASはソースで使うことはもちろんだが、計算機毎に最適化されたバイナリライブラリとしてリンクすることも可能なので、パーソナルコンピュータからベクトル型スーパーコンピュータまでカバー出来る。また、JASON2は直接編成の作業ファイルを常に使う仕様であったが、ITRANSではgij,rsがメモリ上に展開出来れば、IFSORTと同様にインコアでの変換処理が自動的に行われるように書かれている。
これから、個々の1/4変換について順に見ていく。処理のコストが最も大きいのは、式(13)の最初の1/4変換のステップである。変換すべき分子軌道の総数をM(≦N)とすると、演算数のオーダーはζN4M/4となる。ここで、ζは基底積分の生存率で、因子1/4は添字正準性に因る。ITRANSでは、この1/4変換の最深部分は以下のようになっている。
DO IWPQ=1,ICOUNT
  IP=IBUFRS(1,IWPQ)
  IQ=IBUFRS(2,IWPQ)
  V=BUFRS(IWPQ)
  IF(IP.EQ.IQ) THEN
    CALL DAXPY(M,V,TACOEF(1,IP),1,GIQRS(1,IP),1) 
  ELSE
    CALL DAXPY(M,V,TACOEF(1,IQ),1,GIQRS(1,IP),1)
    CALL DAXPY(M,V,TACOEF(1,IP),1,GIQRS(1,IQ),1)
  END IF
END DO
このpqループの外側には、IFSORTでソートされた積分バッファを読み込んで駆動される rsループがあるので配列GIQRSはrsの次元を持つ必要はない。言い換えると、上のループ処理の前にGIQRS領域のクリア操作(零で初期化する)が要る。TACOEFは変換に係る軌道係数cpiを転置し、第一次元を軌道、第二次元を基底とした配列c'ipである。これにより、配列アクセスのストライドは1となり、DAXPY内でのループアンロールによる加速が有効となる。IPとIQが違う場合、正準性から二つの処理が並んでいる。
続く式(14)の2段目の変換は次の通り。
IJ=0
DO I=1,M
  DO J=1,I
    IJ=IJ+1
   S=DDOT(N,ACOEF(1,J),1,GIQRS(I,1),M)
    IF(DABS(S).GT.THIJRS) THEN
      IZ=IZ+1
      GIJRS(IZ,IJ)=S
      LGIJRS(1,IZ,IJ)=IR
      LGIJRS(2,IZ,IJ)=IS
      IF(IZ.EQ.LBUFRS) THEN
        IREC=IREC+1
        WRITE(Dir.) IREC,IRECPR,GIJRS(IJ),LGIJRS(IJ)
        IZ=0
      END IF
    END IF
  END DO
END DO
IRとISは外側で走っているrsループの変数、ACOEFは軌道係数(cpi)を基底-軌道の順に含む配列である。基底数Nに関するDDOT処理のGIQRSのストライドは、軌道数Mとなっている。
前半、2/4変換の段階でgij,rsとなった時点で閾値判断があり、生き残ったものはスタック変数IZをインクメントして直接編成ファイルへのバッファに入れる。IZがバッファ長LBUFRSに達すると、前のレコード番号IRECPRと共にそのIJについてファイルに書き出す。これが、吉嶺の連鎖アルゴリズム[18]の要である。変換と書き込みの処理を最後のrsまで行い、バッファをファイルにフラッシュすれば、後半のrs®kl変換に必要なijに関するソート処理も済んでいる。
変換の後半は、直接編成ファイルからij 毎にgij,rsを吉嶺の連鎖により読み出して駆動されるループが外側になる。式(15)の3/4段目の変換は、DAXPYを使って次のように書ける。
DO IWRS=1,ICOUNT
  X=GIJRS(IWRS)
  IR=LGIJRS(1,IWRS)
  IS=LGIJRS(2,IWRS)
  IF(IR.EQ.IS) THEN
    CALL DAXPY(M,X,ACOEF(IR,1),N,GIJKS(IR,1),N)
  ELSE
    CALL DAXPY(M,X,ACOEF(IS,1),N,GIJKS(IR,1),N)
    CALL DAXPY(M,X,ACOEF(IR,1),N,GIJKS(IS,1),N)
  END IF
END DO
演算コストは、gij,rsの閾値生存率をlとしてlN2M3/2となる。式(13)の最初の1/4変換の時と同様、上のループの前に配列GIJKS領域のクリアが必要である。
式(16)、すなわち最終4/4段目の変換は、DDOTによって配列アクセスのストライドを1として処理出来る。
KL=0
DO K=1,M
  DO L=1,K
    KL=KL+1
    IF(IJ.GE.KL) THEN
      SS=DDOT(N,GIJKS(1,K),1,ACOEF(1,L),1)
      IF(ABS(SS).GT.THIJKL) THEN
        IZZ=IZZ+1
        GIJKL(IZZ)=SS
        LGIJKL(1,IZZ)=I
        LGIJKL(2,IZZ)=J
        LGIJKL(3,IZZ)=K
        LGIJKL(4,IZZ)=L
        IF(IZZ.EQ.LR2BUF) THEN
          WRITE(Seq.) LR2BUF,GIJKL,LGIJKL
          IZZ=0
        END IF
      END IF
    END IF
  END DO
END DO
閾値以上のgij,klは、添字を付けてバッファに入れてパッキングを行いつつ、順編成ファイルに書き出す。ijのループを抜け、残っているバッファ領域をファイルにフラッシュすれば、ITRANSモジュールによる2電子積分の変換処理は終了する。
付記しておけば、gij,klで添字の範囲をik を仮想軌道に関するabに、jl を占有軌道に関するmnとすれば、MP2-PTによる電子相関のエネルギー

の計算に必要な積分は揃う(MP2のコストのオーダーがNの5乗であるのは、積分変換に拠る)。言い換えると、ITRANSに僅かの修正を施すだけでMP2計算のプログラムとして使える。

4. 3 密度局在化用PLOCAL

HF計算で通常得られる正準分子軌道は、Fock演算子の固有関数(固有値が軌道エネルギー)となっているが、振幅は分子全体に広がっている。この状況は、小型分子だけでなく中大型の分子でも共通しており、生体分子も例外ではない。
閉殻のHF配置を表わすSlater行列式はユニタリー不変性[1 - 5]を持っているので、占有軌道と仮想軌道のそれぞれの空間内で局在化、すなわち分子の特定領域に振幅が集中している局在化軌道に変換することがエネルギー不変のままで可能である。距離の離れた局在化軌道間の重なりは小さくなり、こうした組に関する2電子積分gij,klは十分小さな値となるはずである(中間的なgij,rsの段階でも小さくなると期待される)。従って、局在化軌道を使うことで、電子相関の記述に有効な積分数、さらには励起配置のタイプを電子対の遠近に応じて低減出来る。これが、局在化相関計算の考え方である[21]。
幸い、CISD、それにCPF/ACPF[7, 8]もユニタリー不変性を持っているので、局在化軌道を使った計算は、積分駆動型のsベクトル計算を行うLCIでは、処理すべき積分の数を減らして計算コストを下げる一つの方策となり得る。デメリットとしては、Fock演算子が対角でなくなるため、結果としてハミルトニアン行列の対角主導性が低下し、sベクトルの反復計算の回数が若干増えてしまうことである。
PLOCAL(Population-based localizationの意)は、IFSORTとITRANSに対する補助的なモジュールとして、PipekとMezey[22]の処方箋、すなわちMulliken電子密度を極大化して軌道局在化を行う機能を持っている。この方法での評価量は、軌道による電子密度の和

である。Aは原子核について走る添字で、式(34)中で片方の和はAに限定されていることに注意して欲しい。実装にあっては、2×2の軌道対毎のユニタリー変換を繰り返していく。基底関数の重なり積分 spqは、1電子積分hpqと同様にGAMESSからファイル経由で得ている。
PLOCALでは、局在化の変換効率を上げるために2×2変換の角度については、補角を含めた幾つかの可能性を一度に評価して、利得が最大となるものを選び出すようにしている。変換の範囲は、占有軌道群と仮想軌道群の各々で、エネルギーの低い内殻を正準軌道のまま保存する指定も出来る。変換された軌道群は、原子密度<iPAi>の大きさ順にソートされ、GAMESSからの正準軌道の係数と同一の形式で別のファイルに書き出される。

4. 4 Linux上での評価

LCIプログラムの開発は、Intel Pentium-III (800MHz) / 384MB-DRAMのパーソナルコンピュータ(Dell製Dimension XPS800R)にLinux RedHat-8[23]を乗せた環境で進めている。デバッグはLinux標準のg77で行い、動作確認後はIntel製のFortranコンパイラであるifc(Version7/Free版)[24]でバイナリを作っている。
Table 1に、テストした系と基底/変換軌道、それにCPU時間(秒)をまとめた。分子は、環状の水の多量体とアミノ酸で、Gaussian[6]を使いHFレベルで構造最適化を行っている。参考比較に用いたGAMESSは、配布時の標準指定に従った最適化のオプション
g77 -O2 -malign-double -fautomatic
でコンパイルしている。GAMESSの積分変換は、基底の数が増えるとファイルによる変換になり、パスが複数回になる(つまり、適当にブロック化して行われる)ので、パスの数をTable 1中、GAMESSのタイミングに付随する括弧内に示している。LCI(ITRANS)は、GAMESSと同じg77オプションとifcの両方でコンパイルしたものを並べている。DAXPYとDDOTも、サイト[20]からソースとしてダウンロードしたものを同時にコンパイルした。また、作業用の直接編成ファイルのバッファサイズLBUFRSは820、変換された積分を書き出すバッファのサイズLR2BUFは1148に取った。積分の閾値は、基底積分gpq,rsに対しては10-10、2/4変換のgij,rsと4/4変換のgij,klの閾値(THIJRSとTHIJKL)はGAMESSに合わせて10-9としている。
g77とifcでタイミングを比較すると、やはりIntel-CPUに最適化された後者での実行は優位に高速である。水の5量体の6-31G**の場合、ifcではg77よりも89秒も速く、DAXPY/DDOTとifcによるコンパイルの相性が良いことが伺える。グリシン2量体の6-31G*の変換は、GAMESSではテストマシン上ではメモリが不足して実行出来なかった。ITRANSはifcコンパイルでは699秒で行えた(前処理のIFSORTでのCPU時間は52秒)が、g77では300秒以上も遅く、水の5量体の場合より差が開いている。ただ、gij,rsを収納する直接編成ファイルと全変換されたgij,klの書き出しの順編成ファイルが1つのドライブに集中したこと、そもそもマシンがやや旧型でディスクが低速であることから、ifcのコンパイルでもCPU時間の倍近い1375秒も経過時間がかかっている。

Table 1. Timings of integral transformation for water polymers and amino acids ona Pentium-III (800MHz) personal computer (refer texts).
Moleculea(H2O)2(H2O)3(H2O)4(H2O)5GlyAlaGly2Gly2
Basis setb6-31G**6-31G**6-31G**6-31G**6-31G6-31G6-31G6-31G*
No. Basis5075100125556897151
No. Orbital466992115506291133
Timingc
LCI / g774.532.0107.2f315.3f8.122.295.01015.3f
LCI / ifc3.423.578.1f225.6f6.016.769.7698.9f
GAMESSd3.828.390.7f(4)240.6f(6)6.5 18.581.6f(3)unable
a The structures of all molecules included here were optimized at HF level with the Gaussian program [6]. The shape of water polymers is of the cyclic type.
b The s-contaminant of cartesian d-polarization functions was projected out from the orbital space.
c CPU timing in seconds. The symbol "f" means a direct-access file-based transformation. The number of passes (or blocked-processing) is indicated in the parentheses for GAMESS.
d GAMESS was compiled with g77 (see texts for options).

GAMESSとITRANSでCPU時間を比べるのは、実装アルゴリズムが異なるのであくまでも参考ではあるが、g77でコンパイルした場合にはGAMESSが速いが、ifcを使えばGAMESSと同等か若干速いというところで、"及第レベル"に達していると言えるだろう。直接編成ファイルを使う場合、経過時間を向上させるには、IOのバッファサイズを最適化する必要があるが、マシン毎にパラメータは違ってくる。Pentium-4でメモリが1GB、7200/rpmといった高速回転ディスクを複数台(RAID0化するなどして)持つような新しいパーソナルコンピュータなら、CPU時間、経過時間共に、現版でも大幅に短縮されるはずであるし、ワークステーションや大型計算機でも同様の改善が見込まれる。もちろん、直接編成ファイルを使う変換の経過時間をきちんと効率化するには、JASON2[17]と同様にIO周りの最適化を行うべきで、今後の課題であろう。
次に、PLOCALによる密度局在化軌道を使った積分数の低減を水とアミノ酸類について調べてみた。占有軌道、仮想軌道共に局在変換を行っている。低減は、軌道範囲から決まる正準添字の積分数に対する比である。Table 2に結果を示す。相関エネルギーを"化学的精度"の下限となる0.1kcal/mol、つまり10-4の桁で見積もる場合、式(4)の分子の収束判定は一桁下の10-5が安全である。式(17)に示すように、積和係数とCIベクトルを乗じてsベクトルの要素に加算されるgij,klの閾値としては、さらに二桁下の10-7は十分なマージンがある値である。この閾値では、6-31G**による水の5量体では局在化軌道では63%まで積分数が落ちる。水多量体のSTO-3G計算では、あまり効果はない。しかし、STO-3Gによるグリシン類の計算では、鎖が伸びていくと局在化軌道を使うと積分数の減少は速く、5量体では27%になっている。トリプトファンは、共役環を持つためか大きさの割に低減の効果は大きくはない。密度局在化の場合、直交性の要求から仮想軌道は占有軌道ほどコンパクトな局在性は得られないことが知られている。特に、6-31++G**のようにSVPの組に広がった関数を補っていくと、その傾向は強くなる。仮想軌道については、別の局在化法を使うことも将来的には検討すべきかもしれない。

Table 2. Reduction with population-based localized orbitals [22] for water polymers andamino acids relative to the formal number of transformed integrals
Moleculea(H2O)4((H2O)4(H2O)5(H2O)5Gly2Gly3Gly4Gly5Trp
Basis setSTO-3G6-31G**STO-3G6-31G**STO-3GSTO-3GSTO-3GSTO-3GSTO-3G
No. Basis281003512553769912287
No. Orbital24923011544638210172
Reductionb
Can. / 10-799.999.8100.0100.0100.099.999.698.999.9
Loc. / 10-999.999.798.697.299.893.076.559.297.8
Loc. / 10-792.285.375.063.590.363.341.027.073.6
a All the geometries of additional molecules to those in Table 1 were similarly optimized by the Gaussian program [6].
b Reduction in percentage relative to the formal number of final transformed integrals with canonical indices. There is no reduction for canonical orbitals with threshold of 10-9.

5 まとめ

本稿では、開発中の汎用Direct-CIプログラムシステムLCIの中で、基底関数から分子軌道の添字へと2電子積分を変換するITRANSモジュール、変換の前処理として基底関数の積分をソートしつつ、Fock型の有効1電子演算子を生成するIFSORTモジュール、及び軌道局在化を行うPLOCALモジュールについて報告した。IO処理の最適化などに改良の余地は残されているが、積分変換の処理速度は、GAMESS[12]と比較して"及第レベル"にあると考えている。
現在、MRCISDエンジン開発前の"試金石"として、閉殻HF配置参照のCISD/CPF[7]を行う限定版を開発しているところである。今後もLCIの開発を鋭意進め、経過報告を行っていきたい。

本研究の一部は、文部科学省ITプログラム「戦略的基盤ソフトウェアの開発(FSIS)」において実施されたものであり、関係各位に深謝する[25]。LCIプログラムは、FSISのプロジェクト成果物として公開される予定である。 JASON2プログラムの2電子積分変換部分の参照を承諾いただいた中京大学の山本茂義助教授、GAMESSの積分周りについてメールで情報をいただいたGAMESS管理者のM. Schmidt博士に感謝する。また、CI計算全般について議論や助言を継続していただいている北海道大学の田中皓教授にも謝意を表しておきたい。

参考文献

[ 1] Methods of Electronic Structure Theory, H. F. Schaefer III ed., Plenum, New York (1977).
[ 2] A. Szabo, N. S. Ostlund, Modern Quantum Chemistry, MacMillan, New York (1982).
[ 3] Lecture Notes in Quantum Chemistry I, B. O. Roos ed., Springer-Verlag, Berlin (1992).
[ 4] Lecture Notes in Quantum Chemistry II, B. O. Roos ed., Springer-Verlag, Berlin (1994).
[ 5] T. Helgaker, P. Joergensen, J. Olsen, Molecular Electronic-Structure Theory, Wiley, Chichester (2000).
[ 6] C. D. Sherrill, H. F. Schaefer III, Advan. Quant. Chem., 34, 143 (1999).
[ 7] R. Ahlrichs, P. Scharf, C. Ehrhardt, J. Chem. Phys., 82, 890 (1985).
[ 8] R. J. Gdanitz, R. Ahlrichs, Chem. Phys. Lett., 143, 413 (1988).
[ 9] H. -J. Werner et al., MOLPRO.
http://www.molpro.net/
[10] K. Andersson et al., MOLCAS.
http://www.teokem.lu.se/molcas/
[11] H. Lischka et al., COLUMBUS.
http://www.itc.univie.ac.at/~hans/Columbus/columbus.html
[12] M. Schmidt et al., GAMESS.
http://www.msg.ameslab.gov/GAMESS/GAMESS.html
[13] M. F. Guest et al., GAMESS-UK.
http://www.dl.ac.uk/CCP/CCP1/gamess.html
[14] T. Helgaker et al., DALTON.
http://www.kjemi.uio.no/software/dalton/dalton.html
[15] Y. Mochizuki, N. Nishi, Y. Hirahara, T. Takada, Theor. Chim. Acta, 93, 211 (1996).
[16] T. Nakano et al., ABINIT-MP.
http://moldb.nihs.go.jp/abinitmp/
[17] S. Yamamoto, U. Nagashima, T. Aoyama, H. Kashiwagi, J. Comp. Chem., 9, 627 (1988).
[18] M. Yoshimine, J. Comp. Phys., 11, 449 (1973).
[19] S. Yamamoto, J. Teraoka, H. Kashiwagi, J. Chem. Phys., 88, 303 (1988).
[20] BLAS.
http://www.netlib.org/blas/
[21] Correlation and Localization, P. R. Surjan ed., Springer-Verlag, Berlin (1999).
[22] J. Pipek, P. G. Mezey, J. Chem. Phys., 90, 4916 (1989).
[23] RedHat Linux.
http://www.redhat.com
[24] IFC.
htttp://www.intel.com/software/products/compilers/flin/
[25] FSIS.
http://www.fsis.iis.u-tokyo.ac.jp/


Return