分割統治(DC)電子状態計算プログラムのGAMESSへの実装

小林 正人, 赤間 知子, 中井 浩巳


Return

1 緒言

分子軌道(MO)法や密度汎関数理論(DFT)などの電子状態理論は,より大きな系を扱う方向と,より高精度に取り扱う方向の2つに沿って発展してきた.この2つの方向は互いが独自に発展してきたが,近年,大規模な系を高精度に取り扱う手法のニーズが高まりつつある.Hartree-Fock (HF)計算およびDFT計算のボトルネックは,2電子反発積分の計算とFock行列の対角化である.これらの形式的な計算コストは,それぞれ基底関数の数Nの4乗および3乗に比例して増加する(O(N4)やO(N3)と記す).2電子反発積分の計算に対しては,例えば,クーロン項の計算コストを減少させる高速多重極展開法(FMM) [1 - 3]などが提案され,実用的に用いられている.また,Fock行列の対角化に関しても様々なアルゴリズムが開発されており[4, 5],分割統治(DC)法[6, 7]もそのひとつである.しかしながら,これまでDC法は,DFTや半経験的MO法などのHF交換項を含まない計算に対してのみ適用されてきた.そこで,我々はDC法をHFおよびハイブリッドDFT計算に適用し,計算時間やエネルギー,状態密度などの再現性について精査を行った[8].また,空間分割に基づくO(N)法による計算は難しいとされているp共役系に対してもDC-HF計算を適用し,効率的な計算が可能なことを明らかにした[9].さらに,DC-HF法及びDC-ハイブリッド DFTに基づく解析的エネルギー微分法の開発も行っており,構造最適化計算も可能となっている[10].しかし,生体分子やナノ材料では,しばしばファンデルワールス力などの弱い相互作用が重要な役割を果たす.HF法や一般的なDFTではこのような相互作用を定量的に取り扱うことは困難なため,2次摂動(MP2)法やクラスター展開(CCSD, CCSD(T))法などのより高精度な計算手法が必要となる.しかし,高精度な計算手法ほど計算コストも増加し,スケーリングもMP2, CCSD, CCSD(T)法でそれぞれ形式的にはO(N5), O(N6), O(N7)と増加する.我々はDC法をMP2やCCSD法などのポストHF法に拡張し,O(N)を達成することに成功した[11 - 13].これにより,これまで従来の方法では計算が難しかった大規模系の電子相関計算が可能になった.
我々は,量子化学計算パッケージGAMESS [14]をベースとしてこれらの計算ルーチンを開発してきた.本論文では,理論的発展とあわせてDC計算プログラム(Figure 11)について報告する.まず2章で,DC法に基づく自己無撞着場(SCF)計算法(DC-HF法およびDC-DFT)と,DC法に基づく電子相関計算法(DC-MP2およびDC-CC法)の理論について述べる.3章ではこれらの手法のGAMESSへの実装について述べ,4章で数値検証の結果を報告する.


Figure 11. Capabilities of our DC program.

2 理論

2. 1 DC法に基づくSCF計算法(DC-HF法およびDC-DFT)

SCF計算に対するDC法の理論について概説する.DC法では,分割行列 pαを用いて,全系の密度行列 Dを部分系aの密度行列 Dαに分割する.

ここで,m およびn は原子軌道(AO)を表し,分割行列は次の形をとる.

分割行列の模式図をFigure 21に示す.部分系の密度行列は,電子数neの閉殻系の場合,Fermi準位eFとヘヴィサイド階段関数(h(x < 0) = 0, h(x >= 0) = 1 )を用いて次のように表される.

DC法では,全系の軌道エネルギーeiおよび軌道係数{Cmi}の代わりに部分系の軌道エネルギーeαiおよび軌道係数{Cαmi}を用いて,部分系の密度行列を近似する.

部分系のMOは,部分系aの領域にその周囲の領域をバッファとして加えた局在化領域の基底{fαm}を用いて次のように展開され,

eαiおよび{Cαmi}は,部分系のFock行列 Fαを対角化することにより決定される.

ここで,部分系のFock行列 Fαは全系のFock行列 Fから局在化領域に対応する部分を切りだした部分行列であり,Sαは局在化領域の基底からなる重なり行列を表す.実際に計算を行う際に扱いやすいように,h(x)をFermi関数fb(x) = [1+exp(-bx)]-1で置き換える.

b は温度の逆数である.Fermi準位はすべての部分系で共通の値を持ち,次の電子数保存の制約条件によって決定される.

全系の密度行列 DDCは,部分系の密度行列をすべて足し合わせることにより得られる.

この密度行列から通常のHF計算と同様の手続きにより全系のFock行列 Fを作成しなおし,密度行列が収束するまで繰り返す.DC-HF計算の全エネルギーは,DDCを用いて従来のHF計算と同様の手続きにより計算される.

DFT計算でも同様の手順で計算可能である.
DC法によるSCF計算の手続きを,Figure 22にまとめた.DC法ではFermi準位を介して電子数が決まっており,部分系毎に電子を割り当てる必要がない.従って,局在化領域の端の部分に水素原子をキャップするなどの人為的な処理を行うことなく,計算を実行することが可能である.


Figure 21. Schematic of partition matrix for DC calculation.


Figure 22. Schematic of the DC SCF procedure.

DC法に基づくHF計算では,式(2.6)を解く際に実行する部分系のFock行列の対角化にかかるコストは,局在化領域の大きさMを用いると形式的にはO(M3)と表され,全系の大きさNには依存しない.部分系の数はO(N)なので,すべての部分系で式(2.6)を解くためにかかる計算コストは理論上O(M3N)となり,十分に大きな系ではM3 << Nなので線形スケーリングが達成される.

2. 2 DC法に基づく電子相関計算法

従来のMP2やCCSD法などの電子相関計算法では,全系の軌道と軌道エネルギーを用いて相関エネルギーを計算する.DC法では,全系の軌道は得ることができないが,部分系にその周囲のバッファ領域を加えた局在化領域で構築された部分系の軌道が得られる.DC計算では軌道占有数にFermi関数を用いているため非整数になる場合があるが,勾配が十分に急なFermi関数を用いているため,部分系の軌道をFermi準位を境に占有軌道と仮想軌道に分類することは可能である.しかし,この軌道をそのまま用いて相関エネルギーを計算すると,バッファ領域に対応する相関エネルギーが二重に見積られることになる.そこで,部分系のみに対応する相関エネルギーを見積るために,著者らが提案したエネルギー密度解析(EDA)[15]を利用する.2.2.1節で相関エネルギーに対するEDAについて述べ,2.2.2節でDC電子相関計算法について述べる.さらに,2.2.3節では,より効率的な計算を可能とするためにHF計算と電子相関計算でバッファの大きさを変えるデュアル・バッファ法について述べる.

2. 2. 1 相関エネルギー密度解析

閉殻系に対する相関エネルギーは,占有軌道i, jと仮想軌道a, bを用いて以下のように表される.

ここで, は有効2電子励起係数を表しており,具体的にMP2法の場合には以下のように与えられる.

式(2.11)中の占有軌道と仮想軌道の割合を決めるパラメータwoccwvirは任意の値を取ることができるが,全エネルギーを保存するために以下の関係を満足する.

EDAでは,マリケンの電子密度解析[16]からの類推で,原子Aに属するAOについてのみ和をとることで,全エネルギーを原子の寄与に分割する.これに従って,原子Aに対する相関エネルギーDEAcorrは,

と定義される.

2. 2. 2 DC電子相関計算法

分割した相関エネルギーを用いることにより,DC法に基づく全系の相関エネルギーは近似的に次のように与えられる.

ここで,iαjαはFermi準位eFよりも小さい軌道エネルギーeαieαjを持つ部分系の占有軌道を表し,aαbαeFよりも大きな軌道エネルギーeαaeαbを持つ部分系の仮想軌道を表す. は部分系aの局在化領域で求め,具体的にMP2法の場合には部分系の軌道と軌道エネルギーを用いて以下のように計算される.

このようにして,DC法に基づいた電子相関エネルギーを得ることができる.特に,MP2の場合にはこの方法をDC-MP2法と呼ぶ.MP2法だけでなく,高次の摂動法や結合クラスター法の場合にも,同様の方法で電子相関エネルギーの見積りが可能である.たとえば,CCSD法の場合には部分系ごとにT1T2励起振幅(tαi,a, tαij,ab)を求め, を以下のように計算する.

たとえばDC-MP2計算では,式(2.17)により部分系の相関エネルギーを求めるための計算コストは,局在化領域の大きさMを用いると形式的にはO(M5)と表され,全系の大きさNには依存しない.部分系の数は全系の大きさNに比例するので,式(2.17)の計算にかかる全計算コストは理論上O(M5N)となり,十分に大きな系ではM5 << Nなので線形スケーリングが達成される.
特に電子相関計算では,計算時間だけでなく大容量のメモリを必要とする.たとえば,CCSD法では2電子励起振幅tij,abの総数はN4に比例して増加する.そのため,それらをメモリに保存するアルゴリズムでは,メモリサイズがボトルネックとなる.一方,DC法に基づく計算では,部分系の2電子励起振幅tαij,abの総数はM4に比例する.そのため,一定の大きさの部分系を用いればその総数は一定であり,必要とするメモリサイズも一定であるという特長を持つ.

2. 2. 3 デュアル・バッファ法

DC計算では,周囲との相互作用を取り込むために部分系を計算する際にバッファ領域をとる.しかし相関エネルギーのオーダーはHFエネルギーよりも3桁程度小さいので,DC計算で同じオーダーのエネルギー誤差を与えるバッファの大きさは,電子相関計算の方がHF計算よりも小さいであろうと予測される.そこで,我々はDC電子相関計算においてHF計算部分とポストHF計算部分で異なるバッファを用いるデュアル・バッファDC電子相関計算法を提案した[13].デュアル・バッファ法の計算手順は,以下の通りである.
  1. 大きいバッファを用いたDC-HF計算を密度とエネルギーが収束するまで解き,全系のFock行列とDC-HFによる全エネルギーを得る.(バッファ無限大の極限で,従来のHF計算を実行しても良い.)
  2. 小さいバッファを用いてDC-HF方程式を一度だけ解き,部分系の軌道と軌道エネルギーを得る.式(2.8)の制約条件から,電子相関計算で占有軌道と仮想軌道を分けるために用いるFermi準位を決定する.
  3. 部分系の軌道と軌道エネルギーを用いて,局在化領域の相関エネルギーを計算する.
  4. 相関エネルギー密度の考え方を用いて,部分系のみに対応する相関エネルギーを抽出する.全ての部分系の相関エネルギーを足し合わせ,全系の相関エネルギーを得る.
  5. HFエネルギーと相関エネルギーを足し合わせ,DC電子相関計算による全系のエネルギーを得る.
一般に相関部分の計算コストはHF部分に比べて非常に大きいため,デュアル・バッファ法を用いることによりさらなる計算コストの削減が期待される.

3 実装

我々は量子化学計算パッケージGAMESSをベースにDC法に基づく計算ルーチンの開発を行ってきた.具体的には,DC法に基づくSCF計算法としてDC-HF法およびDC-DFTとこれらのエネルギー勾配法を実装し,さらにDC法に基づく電子相関計算法としてDC-MP2およびDC-CC法を実装した.
本プログラムを用いてDC計算を行う場合のインプットファイルの例をFigure 31に示す.Figure 31 (a)はDC-HF計算のインプット例であり,太字の部分が通常のHF計算のインプットとは異なる部分である.DC法に基づくSCF計算を行うためには,新たに$DANDCネームリストを追加し,DC計算を実行するフラグ(DCFLG=.TRUE.)を指定する.$DANDCネームリストでは,部分系とバッファ領域を指定することができる.部分系の作り方はSUBTYPで指定することができ,1原子(ATOM),手動入力(MANUAL),自動決定(AUTOおよびAUTBND),カードからの読み込み(CARD)の中から選択可能となっている.SUBTYP=AUTOがデフォルトであり,この場合は原子間距離に従って入力座標を一意に回転させ,全系を格子状の空間に分割して部分系を決定している.さらに,SUBTYP=AUTBNDでは,以下のアルゴリズムに従って部分系の調整を行っている.


Figure 31. Sample input for (a) DC-HF and (b) DC-MP2 calculation. The words in boldface and those with underline are required to execute DC and DC-based correlation calculation, respectively.

  1. 結合距離から結合次数を判定
  2. 多重結合を切断しないように調整
  3. 水素原子とその最近接原子の結合を切断しないように調整
バッファ領域の作り方はBUFTYPで指定できる.デフォルトでは部分系に含まれる原子を中心とする球状領域に含まれる全ての原子を採用している(RADIUS).球の半径はBUFRADで指定する.また,原子ではなく部分系ごとに区切られるようにバッファ領域を決めるBUFTYP=RADSUBや,カードから読み込むBUFTYP=CARDも可能である.さらに,初期軌道の新たな構築法として,DC法に基づいた拡張ヒュッケル計算を行うGUESS=HUCSUBを$GUESSネームリストに加えた.
Figure 31 (b)には,DC法に基づく電子相関計算の場合の例として,DC-MP2計算のインプットを示した.MP2計算を実行することは,通常と同様に$CONTRLネームリストでMPLEVL=2を加えて指定する.太字の下線部で示した$DCCORRネームリストが,DC法に基づく電子相関計算に新たに必要な部分である.電子相関に対するバッファ領域の半径はRBUFCRで指定され,BUFRAD ($DANDCネームリスト)と異なる場合にはデュアル・バッファ法による計算が行われる.通常のHF計算の後に電子相関だけをDC法を用いて計算するには,DODCCR=.TRUE.を追加する.また,FZCORE=.TRUE.とすれば,電子相関計算で内殻軌道を凍結させることも可能である.
部分系とバッファの情報をカードで指定する(SUBTYP=CARD, BUFTYP=CARD)ためには,SCF計算では$SUBSCFグループを,電子相関計算では$SUBCORグループを追加する.Figure 32には$SUBSCFグループの指定例を示した.エクスクラメーションマーク(!)で始まる行はコメントであり,説明のために加えている.カードで指定する場合には,まず部分系の数を$DANDCネームリストのNSUBSで指定する.$SUBSCFグループでは,部分系1・部分系1のバッファ領域・部分系2・部分系2のバッファ領域・・・と順番に原子の番号を指定する.部分系の終わりには0を1つ,バッファ領域の終わりには0を2つ挿入する.この例では,部分系1は原子1から原子5であり,バッファ領域は原子6から原子13ということになる.また,番号が連続する場合,最初の原子の番号と最後の原子の番号に負号(-)を加えたものを入力することで指定することが可能である.Figure 32の部分系2の場合は,6と-9が入力されているので原子6から原子9が指定される.バッファ領域は原子1から5と10から17が指定されていることになる.この入力方法は組み合わせて用いることも可能である.また,電子相関計算で用いる部分系とバッファを指定する$SUBCORグループも同様のフォーマットである.さらに,SUBTYP=MANUALでは部分系だけを手動で入力することができる.この場合,カードから読み込みの場合と同様に$DANDCの中のNSUBSで部分系の数を入力し,配列LBSUBSで原子1から順に所属する部分系の番号を指定することができる.


Figure 32. Sample input for specifying subsystems and buffer regions at $SUBSCF group.

エネルギー以外のプロパティや解析手法としては,Mulliken電荷密度解析[16]をはじめ,EDA [15],1電子状態密度,双極子モーメントの計算などが可能である.また,DC法では全系の正準軌道を求めることはできないが,各軌道に相当する電子密度を求めることが可能である.

4 数値検証

GAMESSに実装したDCプログラムを用いて,DC計算の数値的検証を行った.すべての計算で6-31G基底関数系を用いた.

4. 1 DC-HF計算

4. 1. 1 バッファサイズ依存性

まず,バッファサイズを変えてb-ストランド型のグリシンペプチドH(NHCH2CO)30OH (以降 (gly)30と略す)のDC-HF計算を行い,バッファに対する依存性を検証した.部分系は1原子(SUBTYP=ATOM),バッファは原子を中心とする球状領域(BUFTYP=RADIUS)とし,半径BUFRADの値を変化させて計算を行った.Figure 41に,バッファサイズに対するDC-HF計算の全エネルギーの誤差(従来法基準)を示す.横軸はBUFRAD,縦軸はエネルギー誤差の絶対値を対数目盛で示したものである.バッファが大きくなるにつれて,DC計算によるエネルギーの誤差はほぼ指数関数的に減少することがわかった.この系の場合にはBUFRAD >= 7 Aで3 mHartree以下(BUFRAD >= 8 Aで1 mHartree以下)の誤差となり,従来のHF計算と同程度の精度で計算可能であることがわかった.


Figure 41. Buffer-size dependence of total energy error of b-strand glycine peptide (gly)30 in absolute value obtained by DC-HF/6-31G calculation.

4. 1. 2 部分系のサイズ依存性

次に,部分系に対する依存性を検証するため,b-ストランド型のグリシンペプチド(gly)30のDC-HF計算を,部分系を1原子,1残基から6残基,及び自動決定(SUBTYP= AUTBND, およそ2原子)と変化させて実行した.バッファは部分系の構成原子から半径6 Aの球状領域(BUFTYP=RADIUS, BUFRAD=6.0)とした.この場合の1繰り返しあたりに要する従来法の対角化ステップに相当する計算時間をFigure 42 (a)に,全エネルギーの誤差をFigure 42 (b)に示した.図の横軸は,部分系に含まれる平均原子数である.計算時間は1ノードのXeon 5060 (3.2 GHz) 1 CPUコアを用いて測定した.部分系の大きさに対して,計算時間は極小値をとり,この系では部分系が1残基(平均原子数が約7原子)の時に最も減少した.DC計算における対角化相当のステップの形式的な計算コストは,局在化領域の大きさMと部分系の数nsubを用いてO(M3nsub)と表される.バッファを固定して部分系のサイズを大きくすると,Mは増加しnsubは減少する.このトレードオフが原因となり,計算時間が極小となるような部分系の大きさが存在する.また,エネルギー誤差は部分系が大きくなるにつれて減少するが,バッファサイズが変化した場合よりも変化はゆるやかであることがわかった.


Figure 42. Subsystem-size dependence of (a) CPU time to solve SCF equation per iteration and (b) total energy error of DC-HF calculation of glycine peptide (gly)30 at 6-31G level.

4. 1. 3 軌道電子密度

DC計算では全系の正準軌道は得られないが,各軌道に相当する電子密度を求めることができる.閉殻系のHF軌道の場合,エネルギーの低い方からi番目の軌道に相当する密度行列は,以下のように表すことが可能である.

ここで

であり,明らかにDmv(eF) = Dmvである.この類推で,DC-HF計算では以下のように全系の軌道に相当する電子密度を求めることができる.

今回開発したプログラムでは,この密度行列を用いて,Gaussian cube file形式などで密度をマッピングすることが可能である.これを用いて,ポリエン鎖C30H32のDC-HF計算を行い,最高被占軌道(HOMO)に相当する電子密度を計算した結果をFigure 43に示す.従来のHF計算によるHOMOの電子密度も示した.部分系は1原子(SUBTYP=ATOM)とし,バッファは球状(BUFTYP=RADIUS)で半径(BUFRAD)は8, 10, 14 Aとした.DC-HF計算による軌道電子密度は,BUFRAD = 8 A程度でも従来のHF計算によるものとほぼ同じ描像を与えている.また,バッファが大きくなるにつれて軌道の節面などがHFのものに近くなっていくことも確認された.


Figure 43. Density map obtained by DC and conventional HF calculation corresponding to the HOMO orbital of polyene chain C30H32 (basis-set: 6-31G).

4. 1. 4 計算時間

次に,系の大きさに対する計算時間の依存性を検証した.系はポリエン鎖CnHn+2 (n = 30-150)とし,部分系はC6H6(端ではC6H7)のユニット,バッファ領域は最近接のユニットとした.用いた計算機は1ノードでCPUはXeon 5060 (3.2 GHz) 1コアである.Figure 44 (a)にDC-HFおよび従来のHF計算の1繰り返しあたりのFock行列生成に要する計算時間を,Figure 44 (b)に1繰り返しあたりの対角化相当の計算に要する時間を示す.横軸はポリエン鎖の炭素数nである.DC-HF計算では,Fock行列の生成・対角化の両方で従来のHF計算よりも計算時間が削減されており,高速化することに成功した.また,従来のHF計算の計算時間は,Fock行列生成と対角化に対してそれぞれO(n2.95)およびO(n3.37)であったが,DC-HF計算ではO(n2.25)およびO(n1.22)になり,スケーリングも改善されている.特に,Fock行列の対角化については,ほぼ線形の計算時間を達成している.Fock行列生成では,DC法により密度行列が疎となるため,交換項の計算時間が削減されている.従ってクーロン項計算を削減するFMMなどを用いることにより,さらなる高速化及びオーダーの減少が可能である.


Figure 44. CPU time for (a) forming Fock matrix and (b) solving SCF equation per iteration of polyene chains CnHn+2 at 6-31G level.

4. 2 DC-MP2,DC-CCSD計算

4. 2. 1 バッファサイズ依存性

DC-MP2計算のバッファサイズ依存性を検証するため,バッファの大きさを変えてポリエン鎖C60H62の計算を行った.部分系はC2H2(端ではC2H3)のユニットとし,バッファは最近接のnb個のユニットまでとるように$SUBSCFで指定した.また,占有軌道と仮想軌道の割合を決めるパラメータは,以前の検討結果[12]から(wocc, wvir) = (1, 0)とした.Figure 45に,バッファのユニット数nbに対するDC-HFエネルギーおよびDC-MP2電子相関エネルギー,DC-MP2全エネルギーの正準MP2エネルギーからの誤差を示す.バッファが大きくなるとエネルギー誤差は減少し,ゼロに近づく.また,DC相関エネルギーの誤差はDC-HFエネルギーの誤差に比べて1桁以上小さいことがわかった.このことからも,2.2.3節で提案したデュアル・バッファ法が有用であることが示唆された.


Figure 45. Buffer-size dependence of deviations of HF, MP2 correlation, and total MP2 energies in DC-MP2/6-31G calculation of polyene chain C60H62 from canonical MP2.

4. 2. 2 デュアル・バッファ法の相関バッファサイズ依存性

次に,デュアル・バッファ法のパフォーマンスを検証するため,相関計算に用いるバッファサイズのみを変化させてポリエン鎖C20H22の計算を実行した.先ほどと同様に部分系はC2H2(端ではC2H3)のユニットとし,HF計算のバッファは最近接の6ユニット(nbHF = 6)となるように$SUBSCFで指定した.相関計算で用いるバッファのユニット数(nbcorr)は1-4と変化させた.Table 41に,このときのDC-MP2およびDC-CCSD計算による相関エネルギーの従来法からの誤差を示す.相関計算に用いるバッファが大きくなると,DC計算による相関エネルギーは,MP2の場合もCCSDの場合も従来法によるものに収束することがわかった.また,DC-MP2計算とDC-CCSD計算の誤差は同程度の大きさであることがわかった.

Table 41. Buffer-size dependence of correlation energy errors [mHartree] of polyene chain C20H22 obtained by DC-MP2 and DC-CCSD calculations at 6-31G level.
nbcorrMP2CCSD
1+19.685+20.864
2+0.300-0.060
3+0.282+0.399
4+0.055+0.219
conv.-1785.601-1957.863

4. 2. 3 計算時間

最後に,DC法に基づく電子相関計算の計算時間を検証した.Figure 46に,ポリエン鎖CnHn+2 (n = 6-120)の炭素数nに対するDC-MP2およびDC-CCSD計算の計算時間を示す.用いた計算機は1ノードでCPUはPentium4/3.8 GHz,メモリ容量は2 GBである.部分系はC2H2(端ではC2H3)のユニットとし,HFおよび相関計算に対するバッファ領域のユニット数は(nbHF, nbcorr) = (6, 2)とした.DC-HF計算の実行に要した計算時間は含まれていない.双対数プロットであるので,傾きが各手法のスケーリングを表している.今回の系では,従来のMP2およびCCSD計算の計算コストはそれぞれO(n5.1)およびO(n5.3)(理論的にはそれぞれO(n5)およびO(n6))であるが,DC法に基づくこれらの計算ではO(n1.2)およびO(n1.3)となり,系の大きさに対してほぼ線形の計算時間を達成している.このように,DC法に基づく電子相関計算では従来法に比べて計算時間が大幅に削減されており,これまで計算が困難であった大規模な系の高精度計算が本プログラムを用いて実行可能となった.


Figure 46. System-size dependence of CPU time for DC and conventional MP2 and CCSD calculations of polyene chains CnHn+2 at the 6-31G level.

5 終わりに

我々は大規模系の電子状態計算を効率的に実行することが可能な分割統治(DC)法の開発を行ってきた.DC法はもともとYangにより半経験的MO計算やDFT計算に適用された手法であるが,我々はHF計算でのアセスメントを行い,また新しくポストHF電子相関計算手法も開発した.本稿ではこれらの手法の量子化学計算パッケージGAMESSへの実装について述べた.電子相関計算では,精度を保ったまま計算速度をさらに向上させることが可能であるデュアル・バッファ法を用いることも可能である.これらの計算を実行するのに必要なインプットファイルは,通常の計算を実行するインプットファイルにわずかな変更を施すだけで作成できることも示した.本プログラムを用いてアセスメントを行い,高精度かつ低コストに大規模系の計算を実行可能であることが検証された.
現在,このプログラムはGAMESS公開版[20]への実装を予定している.また,さらに高精度な計算手法で,現在の電子状態計算のgold standardと言われるCCSD(T)法への拡張も進めている[17 - 19].しかし,CCSD(T)法の3電子励起補正エネルギーは式(2.11)のような形ではなく,さまざまな分割が考えられるため,現在は予備的な計算にとどまっている.アセスメントが十分に行われれば,本手法で非常に高精度な大規模系計算が可能となるであろう.

本研究は,文部科学省次世代スーパーコンピュータプロジェクト,同省科学研究費・特定領域「実在系の分子理論」,同省Global COE「実践的化学知」および早稲田大学理工学研究所(RISE)の支援を受けて実施された.また,量子化学計算の一部は自然科学研究機構計算科学研究センター(RCCS)の計算機を利用して行った.赤間は日本学術振興会特別研究員奨励費の助成を受けて研究を行った.

参考文献

[ 1] C. A. White and M. Head-Gordon, J. Chem. Phys., 101, 6593 (1994).
[ 2] M. C. Strain, G. E. Scuseria, and M. J. Frisch, Science, 271, 51 (1996).
[ 3] C. H. Choi, K. Ruedenberg, and M. S. Gordon, J. Comput. Chem., 22, 1484 (2001).
[ 4] S. Goedecker, Rev. Mod. Phys., 71, 1085 (1999).
[ 5] S. Y. Wu and C.S. Jayanthi, Phys. Rep., 358, 1 (2002).
[ 6] W. Yang, Phys. Rev. Lett., 66, 1438 (1991).
[ 7] W. Yang and T.-S. Lee, J. Chem. Phys., 103, 5674 (1995).
[ 8] T. Akama, M. Kobayashi, and H. Nakai, J. Comput. Chem., 28, 2003 (2007).
[ 9] T. Akama, A. Fujii, M. Kobayashi, and H. Nakai, Mol. Phys., 19-22, 2799 (2007).
[10] H. Nakai, D. Sakura, A. Fujii, T. Akama, and M. Kobayashi, in preparation.
[11] M. Kobayashi, T. Akama, and H. Nakai, J. Chem. Phys., 125, 204106 (2006).
[12] M. Kobayashi, Y. Imamura, and H. Nakai, J. Chem. Phys., 127, 074103 (2007).
[13] M. Kobayashi and H. Nakai, J. Chem. Phys., 129, 044103 (2008).
[14] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis, and J. A. Montgomery, J. Comput. Chem., 14, 1347 (1993).
[15] H. Nakai, Chem. Phys. Lett., 363, 73 (2002).
[16] R. S. Mulliken, J. Chem. Phys., 23, 1833 (1955).
[17] 中井浩巳, 化学と工業, 61, 961 (2008).
[18] 小林正人,中井浩巳, 化学, 64(1), 38 (2009).
[19] 中井浩巳, 早稲田産学連携レビュー2009, 大学出版グループ編, 日経BP (2009), pp. 72-73.
[20] http://www.msg.chem.iastate.edu/gamess/


Return