分割統治量子化学計算におけるバッファ領域決定の自動化

Abstract

本稿では，著者らが開発してきた大規模量子化学計算手法である分割統治（DC）法において，バッファ領域を自動的に決定する手法について述べる．バッファ領域は，DC法の近似に伴い導入される誤差に直接関係し，その選択はエネルギー精度を決める要である．繰り返し計算であるDC Hartree-Fock法では，二層の階層構造を持つバッファ領域を用いる．外側バッファ領域の各原子からのエネルギー寄与を概算し，エネルギー閾値に基づいてバッファ領域をその方向に拡大すべきか否かを判断することで，バッファ領域を徐々に拡大していく．一方，繰り返し計算ではない2次Møller-Plesset摂動計算に対するDC法では，元のバッファ領域内の各原子に対してエネルギー寄与を概算し，与えられたエネルギー閾値以上の寄与を持つ原子のみをバッファ領域に残す．いずれの手法も，エネルギーに基づく1つの閾値だけをパラメータとして，ほぼ一定の精度でエネルギーを計算できることを実証した．

Translated Abstract

A scheme to automatically determine the buffer region in the divide-and-conquer (DC) large-scale quantum chemical method is introduced. The buffer region directly relates to the error introduced by the DC method. In the iterative DC Hartree-Fock procedure, the automatic scheme adopts two-layered buffer region and gradually enlarges the buffer region by evaluating the energy contribution from the outer buffer region and determining whether the buffer region should be enlarged or not based on the energy-based threshold. On the other hand, in the non-iterative DC second-order Møller-Plesset perturbation calculation, the energy contribution is approximately estimated for the atoms in the buffer region and only those atoms that contribute more than an energy-based threshold are left in the buffer region. We demonstrated that both methods achieve almost constant accuracy in the energy using only one energy-based threshold as a parameter.

1 緒言

まずは日本コンピュータ化学会設立20周年に際し，心より祝意を表したい．第一著者の小林が研究を始めてまもなく参加した学会が，2002年に東京都北区の北とぴあで開催された日本コンピュータ化学会の最初の春季年会であった．東京都北区は小林の地元であり，北とぴあは勝手知ったる会場であった．当時，早稲田大学の中井浩巳先生の研究室に配属されたばかりの学部4年生だった小林は，ポスター会場設置や受付等の手伝いに駆り出された．懇親会終了後には，長嶋雲兵副会長(当時は事務局長)が先頭に立って2次会会場の居酒屋へと繰り出し，店のワインを枯らしたことを鮮明に記憶している．それから20年，つまり小林の研究生活も20年を迎える，と思うと感慨深い．

これまでもたびたび，春・秋の年会で発表の機会をいただいた．ここには改めてリストしないが，これまでの発表タイトルを見てみると(最初の春季年会の翌年2003年以降，少なくとも共著者としての発表はほぼ毎年行っている)，我ながら自身の研究が右往左往しているさまがそのまま反映されていると感じる．しかし，発表のたびに貴重なコメントやヒントを得て次の研究につなげており，この学会に支えられてこれまでの研究を展開することができたことと強く実感する．2009年にはリニアスケーリング量子化学計算手法である分割統治(DC)法のGAMESSプログラムへの実装をまとめた日本コンピュータ化学会誌の論文 [1]に対して吉田賞(論文賞)を賜り，受賞講演の光栄に浴した．その後もDC法に関する開発は継続しているが，2011年春季年会の報告以来，第一著者としては本会年会では報告していないようである．

DC法は，考えている系をいくつもの部分系に分割して計算を行うことで系の大きさに対して線形の計算時間(リニアスケーリング)で量子化学計算を実行する手法で，いわゆるフラグメント分割型の高速計算法の一つである [2, 3]．DC法自身の提案は1991年デューク大学のYangら [4, 5]によるものであるが，同1991年には広島大学の今村らによりElongation法 [6, 7]が，1999年には産総研の北浦らによりフラグメント分子軌道(FMO)法 [8, 9]が提案されるなど，この分野は同時期から日本がリードをして研究を行っている．FMO法とDC法の大きな違いは，分子軌道(MO)の構築領域にある．FMO法では，全系を重なりなく分割したフラグメントの中だけでMOを構築して一体のエネルギーを計算し，フラグメントペアに対して同様の計算を行うことでこれを補正する．必要に応じて，三体 [10]，四体 [11]の寄与まで計算されることもある．誤解のないように言うと，各フラグメントの計算では，フラグメント外部からの静電相互作用も考慮している．一方DC法では，ペア以上の計算を行わない代わりに，重なりなく分割したフラグメント(これをDC法では中央領域と呼ぶ)に周囲の領域(バッファ領域と呼ぶ)を予め加えてMOを構築する(ここでも，バッファ領域の外側の静電相互作用は考慮されていることを注記しておこう)．DC法の計算精度は，このバッファ領域の形状，特に大きさに顕著に依存する．GAMESSには，中央領域の各原子から半径r_b以内の原子をバッファ領域とする1パラメータのオプションがあるが，適切なバッファ領域の大きさは系の性質によって異なる．

ところで，実践的な量子化学計算のほとんどが密度汎関数理論(DFT)ベースで行われるようになって久しいが，系統的で正しい理由により精度を改善できるpost Hartree-Fock (HF)波動関数理論の需要も決してなくなるものではない．特に大規模系のpost HF計算はDFT計算の正しさを保証するための試金石としての役目も果たす．我々もDC法に基づく大規模post HF電子相関計算法を開発してきた [12,13,14,15]．この方法では，バッファ領域まで広がったMOを用いて中央領域の相関エネルギーを見積もるが，やはり計算精度はバッファ領域の形状に依存する．さらに，考慮する相互作用の質が異なるため，HF計算とpost HF計算でも適切なバッファ領域の形状は異なってくる [16]．このため，精度とコストのバランスを取った適切なDC法による計算のためには，事前の評価が必要であった．

そこで我々は最近，求めるエネルギー誤差の上限から適切なバッファ領域を自動で決定する方法を開発してきた．HF法をはじめとする一体近似計算では，繰り返しを要する自己無撞着場(SCF)計算の過程でエネルギー閾値に基づきバッファ領域を拡大する方法を開発した [17]．一方，電子相関計算に対しては，予め設定した大きめの初期バッファ領域内の各原子に対して相関エネルギーの寄与を見積もり，閾値よりも小さい寄与と判定された原子を実際の部分系の計算から除外する方法を開発した [18]．これにより，エネルギーという最も重要なパラメータに基づく閾値でDC計算の精度を制御することが可能になる．また，距離をベースにしたバッファ領域の指定では難しい異方性の考慮も可能となる．本稿では，DC法の最近の展開として，このバッファ領域の自動決定法について概観するとともに，エネルギー以外のプロパティ値の所望誤差に基づいたバッファ領域の決定についても少し触れよう．

2 DC-SCF計算

2.1 理論

本稿では基本的に閉殻系を扱う．非制限開殻計算に対する拡張は文献 [19]等を参考されたい．

DC法では原子軌道(AO)を基底関数として用いることを前提としており，各基底関数は一つの原子に属している．全系をまずN_sub個の重なりのない中央領域に分割する．部分系αの中央領域に対応する基底関数の集合を以後S(α)で参照する．各中央領域にその周囲の領域をバッファ領域として加え，実際に部分系のMOを計算する局在化領域を構築する．DixonとMerzによる計算スキーム [20]では，内側と外側の2層のバッファ領域を利用する(Figure 1)．内側バッファ領域(この部分に対応する基底関数の集合をB_i(α)で参照する)は部分系のMOの構築に利用されるほか，密度行列にも寄与する．しかし，外側バッファ領域(この部分に対応する基底関数の集合をB_o(α)で参照する)は部分系のMOの構築にのみ利用され，密度行列には寄与しない．

Fig. 1.

 Structure of central, inner buffer, and outer buffer regions in the two-layer DC method.

DC-SCF法では，全系の一体密度行列を部分系の寄与の和として計算する．   

$D_{\mu \nu }\approx D_{\mu \nu }^{\text{DC}}=\sum _{\alpha }^{\text{subsystem}}{P}_{\mu \nu }^{\alpha }{D}_{\mu \nu }^{\alpha }$

(1)

ここで$D^{\alpha }$は以下で与えられる部分系αの密度行列である．   

$D_{\mu \nu }^{\alpha }=\sum _p{f}_{\beta }({\epsilon }_F-{\epsilon }_p^{\alpha })C_{\mu p}^{\alpha }{C}_{\nu p}^{\alpha *}$

(2)部分系の分子軌道$\left\{{\psi }_p^{\alpha }\right\}$は中央領域に内側バッファ領域と外側バッファ領域を加えた「外側局在化領域」の基底関数$\{{\varphi }_{\mu };\mu \in {\text{L}}_{\text{o}}(\alpha )\equiv \text{S}(\alpha )\cup {\text{B}}_{\text{i}}(\alpha )\cup {\text{B}}_{\text{o}}(\alpha \left)\right\}$を用いて展開する．   

${\psi }_p^{\alpha }\left(r\right)=\sum _{\mu \in {\text{L}}_{\text{o}}\left(\alpha \right)}{C}_{\mu p}^{\alpha }{\varphi }_{\mu }\left(r\right)$

(3)部分系のMO係数$\left\{C_p^{\alpha }\right\}$と軌道エネルギー$\left\{{\epsilon }_p^{\alpha }\right\}$は，以下の部分系のRoothaan方程式を解いて得られる．   

$F^{\alpha }\left[D^{\text{DC}}\right]C_p^{\alpha }={\epsilon }_p^{\alpha }{S}^{\alpha }{C}_p^{\alpha }$

(4)$F^{\alpha }\left[D^{\text{DC}}\right]$は密度行列$D^{\text{DC}}$に対して求められる有効ハミルトニアン行列(以下はHF法の場合)   

$F_{\mu \nu }\left[D^{\text{DC}}\right]=H_{\mu \nu }^{\text{core}}+\sum _{\lambda \sigma }{D}_{\lambda \sigma }^{\text{DC}}\left[2\left(\mu \nu \right|\sigma \lambda )-(\mu \lambda \left|\sigma \nu \right)\right]$

(5)の，$S^{\alpha }$は重なり積分   

$S_{\mu \nu }=\left({\varphi }_{\mu }\right|{\varphi }_{\nu })$

(6)の部分系αの外側局在化領域に対応する部分行列を表している．$\left(\mu \nu \right|\sigma \lambda )=\iint d{r}_{1}d{r}_2{\varphi }_{\mu }^{*}\left(r_1\right){\varphi }_{\nu }\left(r_1\right)r_{12}^{-1}{\varphi }_{\sigma }^{*}\left(r_2\right){\varphi }_{\lambda }\left(r_2\right)$は二電子積分である．(1)式の$P^{\alpha }$は分割行列であり，2層のバッファ領域を用いるDC法の場合は以下で与えられる．   

$P_{\mu \nu }^{\alpha }=\{\begin{array}{c}1\text{for }\mu \in S\left(\alpha \right)\wedge \nu \in S\left(\alpha \right)\\ 1/2\text{for (}\mu \in S\left(\alpha \right)\wedge \nu \in B_{\text{i}}\left(\alpha \right))\text{ and }viceversa\\ 0\text{otherwise}\end{array}$

(7)(2)式の$f_{\beta }\left(x\right)={[1+\exp (-\beta x\left)\right]}^{-1}$はFermi分布関数であり，ε_Fは以下の全電子数(n_e)保存の条件から決定されるFermi準位である．   

$n_{\text{e}}=2\text{Tr}\left(D^{\text{DC}}S\right)$

(8)DC-SCF計算では，(1)式の密度行列と(5)式の有効ハミルトニアン行列が自己無撞着となるように繰り返し計算を行う．

2.2 DC-SCFバッファ領域の自動決定

上述の2層バッファ領域の方法は，もともと主にSCF収束性の向上を目的に導入されたものである．半経験的MO計算では，DC法を用いない「通常計算」からのエネルギー誤差の減少が見られるという報告 [20]もあるため，我々がDC法をGAMESSに実装した際にもマニュアル非掲載のオプションとして2層バッファ領域が利用できるようになっていた．しかしHFやDFT計算では，外側バッファ領域を用いずに内側バッファ領域だけで計算した方がエネルギー誤差は小さくなることを経験的に知っていた．そこで，外側バッファ領域を内側バッファ領域に組み込んだ時のエネルギー変化を見積もろう．その際，密度行列がΔD変化するとし，SCFの影響は無視してHFエネルギーの変化を1次で近似すると   

$\Delta E=2\text{Tr}\left[\Delta DF\left[D^{\text{DC}}\right]\right]$

(9)と計算される．外側バッファ領域を内側バッファ領域に組み入れてもFermi準位が変化しないと仮定すると，ΔDは   

$\Delta D_{\mu \nu }\approx \Delta D_{\mu \nu }^{\text{DC}}=\sum _{\alpha }^{\text{subsystem}}\Delta D_{\mu \nu }^{\alpha }=\sum _{\alpha }^{\text{subsystem}}\Delta P_{\mu \nu }^{\alpha }{D}_{\mu \nu }^{\alpha }$

(10)で表される．ここで$\Delta P_{\mu \nu }^{\alpha }$は以下で与えられる．   

$\Delta P_{\mu \nu }^{\alpha }=\{\begin{array}{c}1/2\text{for (}\mu \in S\left(\alpha \right)\wedge \nu \in B_{\text{o}}\left(\alpha \right))\text{ and }viceversa\\ 0\text{otherwise}\end{array}$

(11)すると，(9)式のエネルギー変化は各部分系αの外側バッファ領域に存在する原子Aの寄与に分割できる．   

$\Delta E=\sum _{\alpha }\sum _{A\in {\text{B}}_{\text{o}}\left(\alpha \right)}\Delta E_A^{\alpha }$

(12)   

$\Delta E_A^{\alpha }=\sum _{\mu \in S\left(\alpha \right)}\sum _{\nu \in A}2D_{\mu \nu }^{\alpha }{F}_{\nu \mu }^{\alpha }\left[D^{\text{DC}}\right]$

(13)(13)式で用いられる密度行列の要素は，中央領域-外側バッファ領域という一部の非対角成分だけであることに注意されたい．密度行列の非対角成分は，導電体でなければ原子間距離に対して指数関数的に減少する．したがって，(13)式が閾値より大きければ，そこからさらにバッファ領域を拡大することでエネルギー誤差が小さくなる．これに対して(13)式が閾値よりも小さければエネルギーの寄与は十分に小さいため，それ以上拡大しなくてもよい，と判断することができる．

以上の考察のもと，DC-SCF計算においてバッファ領域を自動的に決定するスキームを開発した．具体的手順は以下のとおりである．

i). 各SCFサイクルで$F^{\alpha }\left[D^{\text{DC}}\right]$を構築後，(13)式にしたがって$\Delta E_A^{\alpha }$を評価する．

ii). 部分系αの外側バッファ領域に存在する全ての原子を内側バッファ領域に移動する．

iii). $\Delta E_A^{\alpha }$が閾値$e_{\text{thresh}}^{\text{SCF}}$以上の外側バッファ領域の原子を中心に，半径r_extの球内の原子を部分系αの新しい外側バッファ領域に含める．

iv). (4)式を用いて部分系のMOを計算し，(1)式と(2)式を用いて密度行列を構築する．上記の手続き後，i)に戻り，SCFサイクルを進める．SCFサイクルが進むと全ての$\Delta E_A^{\alpha }$が閾値以下になり，最終的に外側バッファ領域は消滅する．この手順によって，部分系当たりのエネルギー誤差を保持しながら，各部分系に対して適切なバッファ領域を選択することができると期待される．

2.3 結果

ここからは，上記のスキームでバッファ領域を自動決定したDC-SCF計算の結果について述べていこう．文献 [17]では主に半経験的分子軌道計算を行った結果を述べたので，本稿では小タンパク質のTrp-cage (PDB ID: 1L2Y)に対してDC-HF計算を行った結果を示そう．各アミノ酸残基を中央領域に取り，バッファ領域拡大を判定するエネルギー寄与の閾値を$e_{\text{thresh}}^{\text{SCF}}$ = 0.1$\text{μ}{E}_{\text{h}}{\text{atom}}^{-1}$，拡大の大きさをr_ext = 3.0 Åに設定した．Figure 2は，今回のスキームを適用したときにGly10部分系(図中でspace fillingで示した)のバッファ領域が拡大する様子を示した．Figure 2(a)が初期バッファ領域(内側バッファ領域の半径$r_{\text{b}}^{\text{in}}=3.5\text{Å}$の場合)を，Figure 2 (b)がDC-HF計算終了時のバッファ領域をball-stick表示している．拡大後の部分系には，Gly10から平均で11.38 Å離れているPro18の原子の約8割がバッファ領域に含まれている．しかし，Gly10からの平均距離が11.21 ÅであるIle4は約3割の原子しかバッファ領域に含まれていない．このように，本スキームでは部分系が適切な方向に異方性を持って拡大していることが確認できる．

Fig. 2.

 Expansion of the buffer region for Gly10 subsystem in Trp-cage system. (a) initial buffer region for $r_{\text{b}}^{\text{in}}=3.5\text{Å}$; (b) final DC-HF buffer region with $e_{\text{thresh}}^{\text{SCF}}=0.1\text{ μ}{E}_{\text{h}}{\text{atom}}^{-1}$; (c) final DC-MP2 buffer region with $e_{\text{thresh}}^{\text{corr}}=10\text{ μ}{E}_{\text{h}}{\text{atom}}^{-1}$.

本スキームが有効にはたらけば，DC法の計算結果は初期バッファ領域の取り方に依存せずに一定以上の精度を担保できるはずである．そこで，初期バッファ領域の大きさに対する結果の依存性を調査した．Table 1に，初期バッファ領域の大きさ(内側バッファ領域の半径$r_{\text{b}}^{\text{in}}$と外側バッファ領域の半径$r_{\text{b}}^{\text{out}}$)を変えてTrp-cageのDC-HF計算を行った時の全エネルギーを示した．通常のHF計算からの1原子当たりのエネルギー差を括弧内に示している．エネルギー誤差は初期バッファ領域の大きさに関わらず，高い精度を与えていることがわかる．閾値の0.1 μE_h atom⁻¹よりもエネルギー誤差が大きくなっているものもあるが，あくまでも見積もりは各部分系ごとに行っているため，ここに示したものと対応するわけではないことに注意されたい．また，各部分系の大きさの目安を表す主軸半径($l_{\text{local}}^{\text{HF,}\alpha }$)の平均値も表に示したが，これについてもどの場合もほぼ一定の値となっていることがわかる．このことからも，初期バッファ領域を深く吟味することなくDC法を適用することができるようになったことが確認できる．

Table 1.  Initial buffer-size dependence of the automated DC-HF energy for the Trp-cage system with 6-31G* basis set.

$r_{\text{b}}^{\text{in}}/\text{Å}$	$r_{\text{b}}^{\text{out}}/\text{Å}$	Energy /Eh	(Diff.)/μEh atom−1	$⟨l_{\text{local}}^{\text{HF,}\alpha }⟩$ /Å
4.0	5.0	−7439.551778	(+0.04)	9.79
4.5	5.5	−7439.551761	(+0.09)	9.65
5.0	6.0	−7439.551722	(+0.23)	9.64
Standard HF		−7439.551790

文献 [17]ではタンパク質のほかに，水分子クラスターやグラフェンといった系に対して適用した結果が示されている．また，本スキームを適用しても計算時間のリニアスケーリングが達成されることも確認しているので，参考されたい．

2.4 エネルギー勾配への適用

量子化学計算で求めたいプロパティ値はエネルギーに限らない．量子化学計算の用途として最も重要な構造最適化に不可欠なエネルギー勾配，すなわちエネルギーの核座標微分は，そのようなプロパティの最たる例である．ここまではDC法の近似に伴うエネルギー誤差の制御に焦点を当ててきたが，エネルギーを望みの精度で得られれば他のプロパティ値も同様に高い精度で得られる，と楽観的に考える理由はない．そこで我々は，DC-SCFエネルギーに加えてエネルギー勾配も望みの精度で得られるようにバッファ領域を調整できないか思案をしている．

まずはDC-HFエネルギー勾配の式を示しておこう．通常のHFエネルギーの核座標Qに関する微分は，   

$\frac{\partial E_{\text{HF}}}{\partial Q}=\sum _{\mu \nu }{D}_{\mu \nu }\left\{2\frac{\partial H_{\nu \mu }^{\text{core}}}{\partial Q}+\sum _{\lambda \sigma }{D}_{\lambda \sigma }\left[2\frac{\partial \left(\nu \mu |\sigma \lambda \right)}{\partial Q}-\frac{\partial \left(\nu \lambda |\sigma \mu \right)}{\partial Q}\right]\right\}+2\sum _{\mu \nu }\frac{\partial D_{\mu \nu }}{\partial Q}{F}_{\nu \mu }$

(14)で表される．ここで(14)式の第1項はHellmann-Feynman (H-F)項，第2項は原子核位置に基底関数中心を置くことで生じるPulay項である．Pulay項は密度行列の(あるいはMOの)核座標微分を必要とするように見えるが，Roothaan方程式の解になっていれば   

$\sum _{\mu \nu }\frac{\partial D_{\mu \nu }}{\partial Q}{F}_{\nu \mu }=-\sum _{\mu \nu }{W}_{\mu \nu }\frac{\partial S_{\nu \mu }}{\partial Q}$

(15)と変形できることが知られている [21, 22]．ここで，Wはエネルギー重み付き密度行列である．YangとLeeが提案したDC-HFエネルギー勾配式 [5]は，(15)式の右辺のPulay項におけるWをDC近似して求めるものである．つまり，   

$W_{\mu \nu }\approx W_{\mu \nu }^{\text{DC}}=\sum _{\alpha }^{\text{subsystem}}{P}_{\mu \nu }^{\alpha }{W}_{\mu \nu }^{\alpha }$

(16)   

$W_{\mu \nu }^{\alpha }=\sum _p{f}_{\beta }({\epsilon }_F-{\epsilon }_p^{\alpha }){\epsilon }_p^{\alpha }{C}_{\mu p}^{\alpha }{C}_{\nu p}^{\alpha *}$

(17)として(15)式の右辺によりPulay項を算出する．著者ら [22]は以前，DC法においては全系の密度行列が変分的に決定されていないために，(15)式の等式が成立しないことを指摘し，上記の方法で求めたエネルギー勾配が大きな誤差を有することを示した．代わりに，各部分系のRoothaan方程式(4)が成立し，変分的に決定されていると仮定すると，以下の改良DC-HFエネルギー勾配式を得ることができる [22]．   

$\frac{\partial E_{\text{HF}}^{\text{DC}}}{\partial Q}=\left(\text{H-F force}\right)-2\sum _{\alpha }\text{Tr}\left[\frac{\partial S^{\alpha }}{\partial Q}{X}^{\alpha }\right]$

(18)

ここで$X^{\alpha }=W^{\alpha }{S}^{\alpha \left[{\text{L}}_{\text{i}}\right(\alpha )\times \text{S}(\alpha \left)\right]}{D}^{\alpha \left[\text{S}\right(\alpha )\times {\text{L}}_{\text{i}}(\alpha \left)\right]}$である．行列の右上にある角括弧$\left[{\text{L}}_{\text{i}}\right(\alpha )\times \text{S}(\alpha \left)\right]$は${\text{L}}_{\text{i}}\left(\alpha \right)\equiv \text{S}\left(\alpha \right)\cup {\text{B}}_{\text{i}}\left(\alpha \right)$に対応する行成分と$\text{S}\left(\alpha \right)$に対応する列成分を取り出した部分行列を表している．

2.2節と同様に，部分系の増大に伴いエネルギー勾配がどれだけ変化するかをざっくりと見積もろう．ここでは，αヘリックス型のグリシン10量体Gly₁₀のDC-HF/6-31G計算を例として取り上げる．まずは，密度行列がΔDだけ変化したときの(14)式第1項(H-F項)の1電子積分部分の変化，すなわち   

$g_Q^{\text{OEI}}=2\sum _{\mu \nu }\Delta D_{\mu \nu }\frac{\partial H_{\nu \mu }^{\text{core}}}{\partial Q}$

(19)を(10)式のΔDに対して計算した(Table 2)．内側・外側バッファ領域は表の$r_{\text{b}}^{\text{in}}$と$r_{\text{b}}^{\text{out}}$で決定して固定し，SCF収束後に(19)式を見積もってエネルギー勾配のH-F項1電子積分部分の通常計算からの実際の誤差と比較した．表には原子数N_atomに対して3N_atom個ある勾配成分の最大絶対偏差(MaxAD)と平均絶対偏差(MAD)を典型的な値として示している．最もバッファ領域が小さい$r_{\text{b}}^{\text{in}}=3.5\text{Å}$のMaxADを除き，バッファ領域を拡大するにつれて実際のH-F項の誤差が減少する様子が確認できる．$r_{\text{b}}^{\text{in}}=3.5\text{Å}$の場合の方が$r_{\text{b}}^{\text{in}}=4.0\text{Å}$の場合よりもMaxADが小さくなっているが，これはlucky cancellationによるものと思われる．見積もられたH-F項の1電子積分部分の誤差($g_Q^{\text{OEI}}$)は実際の誤差よりも1~2桁程度小さくなっており，また，$g_Q^{\text{OEI}}$は実際のH-Fの誤差よりも急速に減少している．これはおそらくエネルギーよりもエネルギー勾配の方が密度行列の変化に鋭敏で，広範囲に影響が及ぶため，外側バッファ領域の厚みが十分でなかったことが大きく影響していると考えている．

Table 2.  Buffer-size dependence of the actual and estimated errors in the one-electron part of the Hellmann-Feynman term of the HF energy gradient of α-helix glycine oligomer Gly₁₀.

$r_{\text{b}}^{\text{in}}$	$r_{\text{b}}^{\text{out}}$	MaxAD/Eh bohr−1			MAD /Eh bohr−1
/Å	/Å	actual	estimated		actual	estimated
3.5	4.5	0.5021	0.2697		0.1032	0.0208
4.0	5.0	0.5795	0.0998		0.0808	0.0072
4.5	5.5	0.0664	0.0076		0.0107	0.0007
5.0	6.0	0.0201	0.0010		0.0037	0.0001
5.5	6.5	0.0160	0.0002		0.0031	0.0000

さらに，密度行列が(10)式のΔDだけ変化したときの(18)式のPulay項の変化，すなわち   

$\begin{array}{c}{g}_Q^{\text{Pulay}}=-2\sum _{\alpha }\sum _{\mu }{\left(\frac{\partial S^{\alpha }}{\partial Q}{W}^{\alpha }{S}^{\alpha \left[\text{L}{}_{\text{o}}\right(\alpha )\times \text{S}(\alpha \left)\right]}\Delta D^{\alpha \left[\text{S}\right(\alpha )\times {\text{L}}_{\text{o}}(\alpha \left)\right]}\right)}_{\mu \mu }\\ =-2\sum _{\alpha }\sum _{\mu }{\left(\frac{\partial S^{\alpha \left[{\text{B}}_{\text{o}}\right(\alpha )\times \text{L}{}_{\text{o}}(\alpha \left)\right]}}{\partial Q}{W}^{\alpha }{S}^{\alpha \left[\text{L}{}_{\text{o}}\right(\alpha )\times \text{S}(\alpha \left)\right]}{D}^{\alpha \left[\text{S}\right(\alpha )\times {\text{B}}_{\text{o}}(\alpha \left)\right]}\right)}_{\mu \mu }\end{array}$

(20)を収束したDC-HF計算の結果に対して見積もり，通常のHF勾配のPulay項と比較して誤差を求めた(Table 3)．Table 3にはエネルギー勾配の通常法からの誤差に加えて，エネルギーの誤差も示している．エネルギー，エネルギー勾配ともに，実際の誤差は初期バッファ領域を大きくするにつれて指数関数的に減少している．また，(20)式により推定したPulay項の誤差のオーダーは実際のPulay項の誤差とほぼ一致しており，(9)式によるエネルギー誤差の推定と同程度に正しく見積もることができている．(20)式のμは外側バッファ領域の原子軌道の和であるため，(12)式と同じように外側バッファ領域の寄与に分割することも可能である．エネルギーで0.0001 E_h(約0.25 kJ mol⁻¹)程度の誤差が許容できるとして，GAMESSの構造最適化計算におけるエネルギー勾配収束判定基準が0.0001 E_h bohr⁻¹であることを考えると，やはり構造を正しく議論するにはより大きなバッファ領域が必要となりそうである．

Table 3.  Buffer-size dependence of the actual and estimated errors in the HF energy and the Pulay term of the HF energy gradient ofα-helix glycine oligomer Gly₁₀. The standard HF energy is −2142.6879 E_h.

$r_{\text{b}}^{\text{in}}/\text{Å}$	$r_{\text{b}}^{\text{out}}/\text{Å}$	Energy error /Eh			MaxAD of Pulay term/Eh bohr−1			MAD of Pulay term/Eh bohr−1
		Actual	estimated		actual	Estimated		actual	estimated
3.5	4.5	+0.0553	+0.0540		0.0265	0.0317		0.0043	0.0027
4.0	5.0	−0.0397	−0.0511		0.0188	0.0209		0.0026	0.0014
4.5	5.5	+0.0027	+0.0022		0.0021	0.0007		0.0003	0.0002
5.0	6.0	+0.0002	−0.0007		0.0007	0.0007		0.0001	0.0001
5.5	6.5	+0.0001	−0.0003		0.0006	0.0003		0.0001	0.0000

議論の余地はまだまだあるものの，今回示したようなプロパティ値の精度保証付き計算を指向したDC計算の自動化も時期に完成すると思われる．勾配以外にも，動的分極率計算 [23]でこのような自動化手法が適用されると，励起周波数で分極率が発散することを利用して励起エネルギーと遷移強度を求められる [24]ため，励起状態計算の精度保証にも利用できると考えられる．

3 DC-MP2エネルギー計算

3.1 理論

ここからは，post-HF電子相関理論，特に2次のMøller-Plesset摂動(MP2)法に対するDC計算について述べていこう．DC post-HF計算では，部分系αに対応する相関エネルギー$E_{\text{corr}}^{\alpha }$を部分系αの分子軌道から計算し，全ての部分系の相関エネルギーを足し合わせることで全系の相関エネルギーE_corrを求める．すなわち，   

$E_{\text{corr}}^{\left(2\right)}\approx \sum _{\alpha }{E}_{\text{corr}}^{\alpha \left(2\right)}$

(21)ここで上付の(2)は本稿で取り上げるMP2相関エネルギーであることを表す．通常の相関エネルギーE_corrはNesbetの定理 [25]により，HF計算で求めた占有MO$\left\{{\psi }_i,{\psi }_j,\mathrm{...}\right\}$と仮想MO$\left\{{\psi }_a,{\psi }_b,\mathrm{...}\right\}$を用いて，以下のように表すことができる．   

$E_{\text{corr}}^{\left(2\right)}=\sum _{i,j}^{\text{occ}}\sum _{a,b}^{\text{vir}}\left(ai\right|bj)\left[2{\tilde{t}}_{ia,jb}^{\left(2\right)}-{\tilde{t}}_{ib,ja}^{\left(2\right)}\right]$

(22)ここで${\tilde{t}}_{ia,jb}^{\left(2\right)}$はアンプリチュードと呼ばれる係数で，MP2法の場合には   

${\tilde{t}}_{ia,jb}^{\left(2\right)}=\frac{\left(ia\right|jb)}{{\epsilon }_i+{\epsilon }_j-{\epsilon }_a-{\epsilon }_b}$

(23)である．DC法では，部分系の相関エネルギーを求めるため，上式のMOを部分系のMOに置き換えるが，各部分系はバッファ領域が互いに重なり合っているため，中央領域に相当する寄与だけを求める必要がある．これには，エネルギー密度解析(EDA) [26]の考え方を利用する．すなわち，(22)式における2電子積分の最後の積分変換の和を中央領域の原子軌道のみに限定し，下式で部分系αの相関エネルギーを計算する．   

$E_{\text{corr}}^{\alpha \left(2\right)}=\sum _{i^{\alpha },j^{\alpha }}^{\text{occ(}\alpha \text{)}}\sum _{a^{\alpha },b^{\alpha }}^{\text{vir(}\alpha \text{)}}\sum _{\mu \in \text{S}\left(\alpha \right)}{C}_{\mu i}^{\alpha }\left(a^{\alpha }\mu \right|b^{\alpha }{j}^{\alpha })\left[2{\tilde{t}}_{ia,jb}^{\alpha \left(2\right)}-{\tilde{t}}_{ib,ja}^{\alpha \left(2\right)}\right]$

(24)部分系αのアンプリチュードも部分系αのMOを用いて計算する．   

${\tilde{t}}_{ia,jb}^{\alpha \left(2\right)}=\frac{\left(i^{\alpha }{a}^{\alpha }\right|j^{\alpha }{b}^{\alpha })}{{\epsilon }_i^{\alpha }+{\epsilon }_j^{\alpha }-{\epsilon }_a^{\alpha }-{\epsilon }_b^{\alpha }}$

(25)DC法のMOはもともと占有MOと仮想MOに分けられてはいないが，たとえばDC-HF計算で決定されたFermi準位よりも軌道エネルギー${\epsilon }_p^{\alpha }$が低いMOを占有軌道，高いMOを仮想軌道とすればよい．また，小数占有数系に対するMP2エネルギーの表式を用いることも可能である [27]．もう一つ，DC-HF計算とDC-MP2計算で同じバッファ領域を使用しなくても良いことを付記しておこう．MP2法で取り込むことができる動的電子相関は局所的なものなので，DC-HF計算よりも小さなバッファ領域で十分な精度の計算を行うことができる [16]．

3.2 DC-MP2バッファ領域の自動決定

EDAに基づけば，部分系αのMP2相関エネルギー$E_{\text{corr}}^{\alpha \left(2\right)}$はさらにこの部分系の局在化領域内の原子の寄与に分割することができる．    (26)      

$E_{\text{corr},B}^{\alpha \left(2\right)}=\sum _{i^{\alpha },j^{\alpha }}^{\text{occ(}\alpha \text{)}}\sum _{a^{\alpha },b^{\alpha }}^{\text{vir(}\alpha \text{)}}\sum _{\mu \in S\left(\alpha \right)}\sum _{\nu \in B}{C}_{\mu i}^{\alpha }{C}_{\nu a}^{\alpha }\left(\nu \mu |b^{\alpha }{j}^{\alpha }\right)\left[2{\tilde{t}}_{ia,jb}^{\alpha \left(2\right)}-{\tilde{t}}_{ib,ja}^{\alpha \left(2\right)}\right]$

(27)MP2法で取り込まれるような動的電子相関の局所性を考慮すると，$E_{\text{corr},B}^{\alpha \left(2\right)}$は部分系αの中央領域から原子Bまでの距離が増大すると急速に減少することが期待される．そこで，$E_{\text{corr},B}^{\alpha \left(2\right)}$の絶対値がある閾値よりも小さいと判断された場合，部分系αのバッファ領域から原子Bを排除しても部分系相関エネルギーの変化は十分に小さくなると考えられる．大規模MP2計算法の一つである，寄与の大きさの判定を困難とする(27)式の軌道エネルギー分母を積分にして分子に移動するLaplace変換MP2法 [28, 29]を用いると，(27)式は   

$E_{\text{corr,}B}^{\alpha \left(2\right)}=-{\int }_0^{\infty }\sum _{\mu \in \text{S}\left(\alpha \right)}\sum _{\nu \in B}\sum _{\lambda \sigma }\sum _{\gamma \kappa \delta ϵ}{X}_{\mu \gamma }^{\alpha }\left(\tau \right)Y_{\nu \kappa }^{\alpha }\left(\tau \right)X_{\lambda \delta }^{\alpha }\left(\tau \right)Y_{\sigma ϵ}^{\alpha }\left(\tau \right)\left(\nu \mu \right|\sigma \lambda )\left[2\left(\gamma \kappa \right|\delta ϵ)-(\gamma ϵ\left|\delta \kappa \right)\right]d\tau$

(28)と変形できる．ここで，X^α(τ)とY^α(τ)は以下のエネルギー重み付き密度行列である(どちらも(17)式のエネルギー重み付き密度行列とは異なるので注意されたい)．   

$X_{\mu \nu }^{\alpha }\left(\tau \right)=\sum _i^{\text{occ}\left(\alpha \right)}{C}_{\mu i}^{\alpha }{C}_{\nu i}^{\alpha }{e}^{({\epsilon }_i^{\alpha }-{\epsilon }_{\text{F}})\tau }$

(29)   

$Y_{\mu \nu }^{\alpha }\left(\tau \right)=\sum _a^{\text{vir}\left(\alpha \right)}{C}_{\mu a}^{\alpha }{C}_{\nu a}^{\alpha }{e}^{-({\epsilon }_a^{\alpha }-{\epsilon }_{\text{F}})\tau }$

(30)正しい相関エネルギーを求めるには，(28)式を数値積分して計算する必要があるが，今回は事前評価が目的であるため，これをGauss-Laguerreの1点求積で近似する．   

$E_{\text{corr,}B}^{\alpha \left(2\right)}~-e\sum _{\mu \in \text{S}\left(\alpha \right)}\sum _{\nu \in B}\sum _{\lambda \sigma }\sum _{\gamma \kappa \delta ϵ}{X}_{\mu \gamma }^{\alpha }{Y}_{\nu \kappa }^{\alpha }{X}_{\lambda \delta }^{\alpha }{Y}_{\sigma ϵ}^{\alpha }\left(\nu \mu \right|\sigma \lambda )\left[2\left(\gamma \kappa \right|\delta ϵ)-(\gamma ϵ\left|\delta \kappa \right)\right]$

(31)   

$X_{\mu \nu }^{\alpha }=\sum _i^{\text{occ}\left(\alpha \right)}{C}_{\mu i}^{\alpha }{C}_{\nu i}^{\alpha }{e}^{({\epsilon }_i^{\alpha }-{\epsilon }_{\text{F}})}$

(32)   

$Y_{\mu \nu }^{\alpha }=\sum _a^{\text{vir}\left(\alpha \right)}{C}_{\mu a}^{\alpha }{C}_{\nu a}^{\alpha }{e}^{-({\epsilon }_a^{\alpha }-{\epsilon }_{\text{F}})}$

(33)近似式(31)が十分妥当であると仮定すると，その絶対値の上限はSchwarzの不等式等を用いて下式で見積もることができる [18]．   

$\left|E_{\text{corr,}B}^{\alpha \left(2\right)}\right|\lesssim e\left(\sum _{\lambda \sigma }\sum _{\delta ϵ}\left|X_{\lambda \delta }^{\alpha }\right|\left|Y_{\sigma ϵ}^{\alpha }\right|A_{\lambda \sigma }^{\alpha }\right)\sum _{\mu \in S\left(\alpha \right)}\sum _{\nu \in B}\sum _{\gamma \kappa }\left|X_{\mu \gamma }^{\alpha }\right|\left|Y_{\nu \kappa }^{\alpha }\right|A_{\mu \nu }^{\alpha }\left[2A_{\kappa \gamma }^{\alpha }\max \right(A^{\alpha }\left)\right]$

(34)ここで$A_{\mu \nu }^{\alpha }=\sqrt{\left|\left(\mu \nu \right|\mu \nu )\right|}$である．(34)式の右辺の丸括弧内は各部分系で定数であるので，次式の$e_B^{\alpha }$が原子Bからのエネルギー寄与の大きさを表すものと考えることができる．   

$e_B^{\alpha }=e\sum _{\mu \in S\left(\alpha \right)}\sum _{\nu \in B}\sum _{\gamma \kappa }\left|X_{\mu \gamma }^{\alpha }\right|\left|Y_{\nu \kappa }^{\alpha }\right|A_{\mu \nu }^{\alpha }\left[2A_{\kappa \gamma }^{\alpha }\max \right(A^{\alpha }\left)\right]$

(35)

上記の$e_B^{\alpha }$を用いて，次のようにDC-MP2計算におけるバッファ領域を自動決定するスキームを開発した．

i). 各部分系に対して，DC-MP2計算の初期バッファ領域を割り当てる．

ii). (35)式の$e_B^{\alpha }$を評価する．

iii). $e_B^{\alpha }$が閾値$e_{\text{thresh}}^{\text{corr}}$以下の原子を部分系αのバッファ領域から除外する．

iv). (4)式を用いて部分系のMOを再構築する．

v). (24)式により部分系の相関エネルギー$E_{\text{corr}}^{\alpha \left(2\right)}$を計算する．全ての部分系の初期バッファ領域の原子に対して$e_B^{\alpha }$を計算するのにかかる時間は，全部分系の数をN_sub，バッファ領域の大きさをmとするとO(N_subm³)である．DC-MP2計算にかかる時間はO(N_subm⁵)であるから，上記の手法によりバッファ領域の大きさをm’ (≤ m)にできれば，大幅な計算時間短縮が期待できる．また，本研究で示す結果では，i)におけるDC-MP2計算の初期バッファ領域は，2章で開発した自動制御型DC-HF計算を用いて決定した．これにより，HF計算からpost-HF計算までの一連の量子化学計算において，各部分系に対して適切なバッファ領域を決定できることになる．

3.3 結果

それでは，上記の方法でバッファ領域を自動決定したDC-MP2計算の結果を示していこう．特記しない限り，DC-HF計算に関する種々のパラメータは2章と同様，DC-MP2バッファ領域決定時のエネルギー閾値は$e_{\text{thresh}}^{\text{corr}}$= 0.1 μE_h atom⁻¹とした．まず，2.3節で示したFigure 2(c)に，DC-HF計算で自動決定されたTrp-cageのバッファ領域(Figure 2(b))をDC-MP2計算の初期バッファ領域として上記のスキームを適用した後の最終的なDC-MP2バッファ領域が示されている(ここでは$e_{\text{thresh}}^{\text{corr}}$= 10 μE_h atom⁻¹を使用した)．中央領域から離れた部分がバッファ領域から除外されていることが確認できる．ちなみに，この部分系に対してFigure 2(b)，Figure 2(c)のバッファ領域を用いて計算したMP2相関エネルギーはそれぞれ−597.328，−597.294 mE_hであり，その差は十分に小さいことも確認された．

次に本手法をポリアセチレン(C_2nH_2n+2)に適用し，全ての$e_B^{\alpha }$を計算するのに要した計算時間を分子サイズnに対してプロットした(Figure 3)．C₂H₂(もしくはC₂H₃)のユニットを中央領域とし，基底関数には6-31G*を，計算機はIntel Xeon E5-2667プロセッサを2基搭載したものを1台使用して計測した．初期バッファ領域は，自動化されたDC-HF法により決定した．両対数プロットを用いたスケーリング解析では，この計算時間は分子サイズn≧25に対してO(n^1.4)でスケールしている．なお，参考文献 [18]では立方体内に水分子をN_water個配置した3次元的な系で同様の測定を行っており，計算時間はO(N_water^1.5)であった．系の次元に関わらず，ほぼ線形スケーリングを達成している．

Fig. 3.

 System-size dependence of the CPU time of the evaluation of $e_B^{\alpha }$ for polyacetylene chain system containing 2n carbon atoms C_2nH_2n+2. The initial inner and outer buffer sizes for the DC-HF calculation were fixed at $r_{\text{b}}^{\text{in}}$=5.0 Å and $r_{\text{b}}^{\text{out}}$= 6.5 Å, respectively.

ここからは，本手法の有効性を異性化エネルギーの計算を例に示していこう．まず，Figure 4のようなポリアセチレンC₁₀₀H₁₀₂のcis-trans異性化エネルギーを従来のバッファ領域固定型のDC-MP2法(DC-MP2 methodの列に$r_{\text{b}}^{\text{MP2}}$でバッファ半径を示した)と今回開発した自動制御型DC-MP2法(同列に$e_{\text{thresh}}^{\text{corr}}$で閾値を示した)を用いて計算した結果をTable 4に示す．原子はすべて平面上に置き，C–C，C=C，C–H結合長は1.462，1.357，1.096 Åで固定した一点計算でcis体とtrans体のエネルギーを計算した．MP2計算の前に2.2節の自動制御型DC-HF計算を初期バッファ領域$r_{\text{b}}^{\text{in}}$= 5.0 Å，$r_{\text{b}}^{\text{out}}$= 6.5 Åで実行し，最終的なDC-HFバッファ領域をDC-MP2計算の初期バッファ領域とした．なお，本計算においては部分系の中で二重結合が切断されると大きな誤差を生むことが分かっているので，CH=CH(あるいはCH2=CH)の単位を中央領域とし，この単位ごとにバッファに含めるか否かを決定した．従来のDC-MP2計算では，バッファ領域を拡大することでtrans体，cis体ともにエネルギー誤差が小さくなり，結果として異性化エネルギーの誤差も小さくなる．これに対して自動制御型DC法では，trans体，cis体のそれぞれに対して適切な大きさのバッファ領域が自動的に選択され，どちらの誤差も同程度となって異性化エネルギーを精度よく見積もることができている．また，MP2バッファ領域の閾値を小さくするにつれて，異性化エネルギーの誤差も小さくなることが確認できる．

Fig. 4.

  cis-trans isomerization of polyacetylene C₁₀₀H₁₀₂.

Table 4.   cis-trans isomerization energy of polyacetylene C₁₀₀H₁₀₂ evaluated at (DC-)MP2/6-31G* level.

DC-MP2 method	Total energy /Eh		Isomerization energy
	trans	cis	/ kcal mol−1
$r_{\text{b}}^{\text{MP2}}$= 6 Å	−3858.088	−3857.721	229.9
$r_{\text{b}}^{\text{MP2}}$= 8 Å	−3858.094	−3857.721	234.0
$r_{\text{b}}^{\text{MP2}}$= 10 Å	−3858.098	−3857.722	235.8
$e_{\text{thresh}}^{\text{corr}}$= 10 μEh	−3858.094	−3857.721	234.1
$e_{\text{thresh}}^{\text{corr}}$= 0.1 μEh	−3858.098	−3857.722	235.9
Standard MP2	−3858.100	−3857.722	237.2

文献 [18]では小タンパク質やオリゴペプチド，水クラスターに適用した結果のほか，MP2相関エネルギーの局所性に関する基礎的な議論も展開しているので，参考されたい．

4 結言

本稿では，大規模量子化学計算手法のひとつである分割統治(DC)法において導入される誤差に深くかかわっているバッファ領域を，エネルギー誤差に基づく閾値パラメータを用いて自動的に決定する手法について紹介した．繰り返し計算であるDC-HF法においては，二層のバッファ領域を用いてバッファ領域を拡大すべきか否かを外側バッファ領域内で判定し，小さなバッファ領域から順次大きくする戦略をとった．一方，繰り返し計算ではないDC-MP2計算においては，動的電子相関の局所性からDC-HF法で必要とされるよりも小さなバッファ領域で十分な精度が得られる，という事実から，バッファ領域内の原子に対して簡単に計算できる指標を用いてエネルギー寄与を推定し，一定以上の寄与を与える原子のみをバッファ領域に残す戦略を採用した．いずれの方法においても，閾値の1パラメータだけでほぼ一定のエネルギー精度を得ることに成功した．

MP2アンプリチュードはより高精度な結合クラスター(CC)計算のアンプリチュードのよい近似になるので，DC-MP2バッファ領域の自動決定手法はDC-CC計算 [13,14,15]においてもそのまま用いることができると考えている．一方，DC-HF計算に対しては，エネルギー以外のプロパティ(今回の例ではエネルギー勾配)の精度を基準としたバッファ領域の構築法について予備的な結果を示した．こちらはまだまだ検討すべき課題が山積している．また，先に述べたように動的分極率から励起状態を求める方向など，応用も広がっている．これらの展開については，今後の日本コンピュータ化学会年会および本誌で発表していきたい．

謝辞

本研究は，JSPS科研費JP25810011及びJP17K05738，文部科学省元素戦略プロジェクト拠点形成型(京都大学) JPMXP0112101003，北海道大学機能強化プロジェクト「フォトエキサイトニクス研究拠点~光励起状態制御の予測と高度利用~」，JST-CREST JPMJCR1902，文部科学省 データ関連人材育成プログラム北海道大学事業(D-DRIVE-HU)の助成を受けて実施された．藤森は北海道大学リーディングプログラムの支援を受けた．文部科学省世界トップレベル研究拠点プログラム(WPI)により設置された北海道大学化学反応創成研究拠点(ICReDD)からも支援を受けた．計算の一部は岡崎の計算科学研究センターと九州大学情報基盤研究開発センターの計算機を利用して行った．この場を借りて謝意を表したい．

参考文献

• [1]   M. Kobayashi, T. Akama, H. Nakai, J. Comput. Chem. Jpn., 8, 1 (2009). doi:10.2477/jccj.H2027
• [2]   T. Akama, M. Kobayashi, H. Nakai, J. Comput. Chem., 28, 2003 (2007). doi:10.1002/jcc.20707 PMID:17455367
• [3]   M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications (Springer, Dordrecht), pp. 97-127 (2011).
• [4]   W. Yang, Phys. Rev. Lett., 66, 1438 (1991). doi:10.1103/PhysRevLett.66.1438 PMID:10043209
• [5]   W. Yang, T. S. Lee, J. Chem. Phys., 103, 5674 (1995). doi:10.1063/1.470549
• [6]   A. Imamura, Y. Aoki, K. Maekawa, J. Chem. Phys., 95, 5419 (1991). doi:10.1063/1.461658
• [7]   F. L. Gu, B. Kirtman, Y. Aoki, in Linear-Scaling Techniques in Computational Chemistry and Physics: Methods and Applications (Springer, Dordrecht), pp. 175-198 (2011).
• [8]   K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi, Chem. Phys. Lett., 313, 701 (1999). doi:10.1016/S0009-2614(99)00874-X
• [9]D. G. Fedorov and K. Kitaura, The fragment molecular orbital method: practical applications to large molecular systems (CRC Press, Baca Raton, 2009).
• [10]   D. G. Fedorov, K. Kitaura, Chem. Phys. Lett., 433, 182 (2006). doi:10.1016/j.cplett.2006.10.052
• [11]   T. Nakano, Y. Mochizuki, K. Yamashita, C. Watanabe, K. Fukuzawa, K. Segawa, Y. Okiyama, T. Tsukamoto, S. Tanaka, Chem. Phys. Lett., 523, 128 (2012). doi:10.1016/j.cplett.2011.12.004
• [12]   M. Kobayashi, Y. Imamura, H. Nakai, J. Chem. Phys., 127, 074103 (2007). doi:10.1063/1.2761878 PMID:17718602
• [13]   M. Kobayahsi, H. Nakai, J. Chem. Phys., 129, 044103 (2008). doi:10.1063/1.2956490 PMID:18681630
• [14]   M. Kobayashi, H. Nakai, J. Chem. Phys., 131, 114108 (2009). doi:10.1063/1.3211119 PMID:19778101
• [15]   T. Yoshikawa, M. Kobayashi, H. Nakai, Int. J. Quantum Chem., 113, 218 (2013). doi:10.1002/qua.24093
• [16]   M. Kobayashi, H. Nakai, Int. J. Quantum Chem., 109, 2227 (2009). doi:10.1002/qua.22111
• [17]   M. Kobayashi, T. Fujimori, T. Taketsugu, J. Comput. Chem., 39, 909 (2018). doi:10.1002/jcc.25174 PMID:29399822
• [18]   T. Fujimori, M. Kobayashi, T. Taketsugu, J. Comput. Chem., 42, 620 (2021). doi:10.1002/jcc.26486 PMID:33534916
• [19]   M. Kobayashi, H. Nakai, Phys. Chem. Chem. Phys., 14, 7629 (2012). doi:10.1039/c2cp40153c PMID:22513877
• [20]   S. L. Dixon, K. M. Merz, Jr., J. Chem. Phys., 107, 879 (1997). doi:10.1063/1.474386
• [21]   P. Pulay, Mol. Phys., 17, 197 (1969). doi:10.1080/00268976900100941
• [22]   M. Kobayashi, T. Kunisada, T. Akama, D. Sakura, H. Nakai, J. Chem. Phys., 134, 034105 (2011). doi:10.1063/1.3524337 PMID:21261328
• [23]   M. Kobayashi, T. Touma, H. Nakai, J. Chem. Phys., 136, 084108 (2012). doi:10.1063/1.3687341 PMID:22380033
• [24]   T. Yoshikawa, J. Yoshihara, H. Nakai, J. Chem. Phys., 152, 024102 (2020). doi:10.1063/1.5124909 PMID:31941302
• [25]   R. K. Nesbet, in Advances in Chemical Physics (John Wiley & Sons, West Sussex) vol. 14, pp. 1-34 (1969).
• [26]   H. Nakai, Chem. Phys. Lett., 363, 73 (2002). doi:10.1016/S0009-2614(02)01151-X
• [27]   M. Kobayashi, T. Taketsugu, Theor. Chem. Acc., 134, 107 (2015). doi:10.1007/s00214-015-1710-y
• [28]   M. Häser, Theor. Chim. Acta, 87, 147 (1993). doi:10.1007/BF01113535
• [29]   M. Kobayashi, H. Nakai, Chem. Phys. Lett., 420, 250 (2006). doi:10.1016/j.cplett.2005.12.088