分割統治法に基づく有限温度型単参照静的相関手法

Abstract

We have developed the divide-and-conquer (DC) based self-consistent-field (SCF) method with finite-temperature (FT) scheme including the electronic entropy contribution, as denoted by FT-DC-SCF. This could approximately involve the static correlation effect appearing in bond-breaking systems, diradicals, and conjugated polymers, within small additional computational costs. This letter reports the energy errors for the conventional and FT-type DC-SCF calculations of polyethylene and polyacetylene. Furthermore, the singlet-triplet separation in polyacene is examined by adopting the conventional and FT-type DC-SCF methods. The numerical data confirm the effectiveness and accuracy of the FT-DC-SCF method.

1 はじめに

エネルギーを電子密度の汎関数から算出する密度汎関数理論(DFT)は，計算コストと精度のバランスが良く，非常に多くの化学者に利用されている．DFTにおいて，エネルギー汎関数の厳密形は知られていないため，近似的なエネルギー汎関数が用いられる．短距離の交換・相関エネルギーについては，近似的なエネルギー汎関数(局所密度近似，一般化勾配近似等)により十分な精度で見積ることができる．また，混成汎関数，領域分割型汎関数 [1,2,3,4]や分散力補正 [5,6,7]を行うことにより，長距離効果についてもある程度の補正が可能である．しかし，DFTは単参照型の計算理論であるので，長距離効果の一つである静的相関効果を取り込むことは非常に困難である．静的相関の効果を取り込むためには多参照法を導入する，あるいは，乱雑位相近似型汎関数 [8, 9]を用いる必要があるが計算コストが高く適用できる分子サイズに制限がある．そこで，静的電子相関を単参照の枠組みで簡便に記述する，非整数占有数を用いる有限温度(FT)法 [10, 11]やThermally assisted occupation法 [12]などが開発され，拡張が期待されている．

筆者らの研究グループでは，大規模計算理論の一つである分割統治(DC)法 [13,14,15,16]の開発を行ってきた．DC法では，分割した部分系が非整数の電荷やスピンをもった場合も取扱い可能とするため，電子温度の概念を導入することにより電子数を非整数で取り扱っている．そのため，有限温度法との高い適応性が期待される．筆者らは最近，電子のエントロピー項まで考慮したFT法をDC法に適用したFT-DC 自己無撞着法(SCF)法 [17]を開発した．本稿ではFT-DC-SCF法をπ共役分子の一重項-三重項エネルギー差に適用し，大規模分子に対する精度検証結果を報告する．

2 理論的背景

FT-HF法では，一電子積分行列H^core，Fock行列Fを用いてエネルギーは以下のようになる．   

$\begin{array}{c}E=\frac{1}{2}\text{tr}\left\{D\left(H^{\text{core}}+F\right)\right\}\\ +2\theta \sum _p\left\{f_p\ln f_p+\left(1-f_p\right)\ln \left(1-f_p\right)\right\}.\end{array}$(1)

ここで{p}は一般分子軌道を表す．通常のHFと異なり，密度行列Dは占有数にFermi関数を用いる．   

$\begin{array}{c}{D}_{\mu \nu }=2\sum _p{f}_p{C}_{\mu p}{C}_{\nu p}^{*},\\ f_p^{}=\left\{1+\exp \left[-\left({\epsilon }_{\text{F}}-{\epsilon }_p^{}\right)/\theta \right]\right\}{.}^{-1}\end{array}$(2)

ここで，Cは軌道係数，εは軌道エネルギー，ε_FはFermi準位，θは擬似電子温度パラメータであり，{μ, ν}は原子軌道を表す．Fermi準位は全電子数を保存するように決定する．   

$N_{\text{e}}=\text{tr}\left(DS\right)=\sum _{p}2f_p^{}.$(3)

式(1)の第二項の電子エントロピー項は，占有数と擬似温度パラメータから算出される．また，我々は，擬似温度パラメータについてChaiらによって提案された系ごとに適した値を決定する手法 [18]を実装した．

DC法では全系を重なりのない中央領域αに分割する．中央領域αの周りのバッファを含めた局在化領域の軌道を構築することで，環境の効果を取り込んでいる．全体の密度行列Dは各局在化領域の分割行列P^α, 分子軌道C^α, 軌道エネルギーε^αと共通のFermi準位ε_Fを用いて構築する．   

$D_{\mu \nu }^{\text{DC}}=2\sum _{\alpha }{P}_{\mu \nu }^{\alpha }\sum _p{f}_p^{\alpha }{C}_{\mu p}^{\alpha }{C}_{\nu p}^{\alpha *}$(4)

,   

$f_p^{\alpha }={\left\{1+\exp \left[-\left({\epsilon }_{\text{F}}-{\epsilon }_p^{\alpha }\right)/\theta \right]\right\}}^{-1}$(5)

.

DC法の各局在化領域の電子数が非整数となるため，FT法と同様にFermi関数を用いている．Fermi準位は全電子数保存の条件から以下の式で決定する．   

$\begin{array}{c}{N}_{\text{e}}=\text{tr}\left(DS\right)=\sum _{\alpha }\sum _p\sum _{\mu \nu }2P_{\mu \nu }^{\alpha }{f}_p^{\alpha }{C}_{\mu p}^{\alpha }{C}_{\nu p}^{\alpha *}{S}_{\mu \nu }^{}\\ =\sum _{\alpha }\sum _{p}2f_p^{\alpha }{w}_p^{\alpha },\end{array}$(6)

  

$w_p^{\alpha }=\sum _{\mu \nu }{P}_{\mu \nu }^{\alpha }{C}_{\mu p}^{\alpha }{C}_{\nu p}^{\alpha *}{S}_{\mu \nu }^{}.$(7)

w^αは重み付き行列であり，部分系の占有数f^αに対する中央領域の寄与を表している．FT法へ拡張する場合，電子エントロピー項は重み付き行列を用いて次のように近似する．   

$\begin{array}{c}{E}_S\approx \sum _{\alpha }{E}_S^{\alpha },E_S^{\alpha }\\ \equiv \sum _{p}2\theta w_p^{\alpha }\left[f_p^{\alpha }\ln f_p^{\alpha }+\left(1-f_p^{\alpha }\right)\ln \left(1-f_p^{\alpha }\right)\right]\end{array}$(8)

DFTの場合は，Hartree-Fock-Roothaan方程式からKohn-Sham-Roothaan方程式に置き換えることにより，同様の手順で計算ができる．

3 結果と考察

FT-DC法を用いてポリエチレンC₄₀H₈₂とポリアセチレンC₄₀H₄₂に対する数値検証を行った．計算条件はBLYP/6-31G** [19,20,21,22]とした．DC計算において，部分系は炭素2個からなるユニット，バッファサイズは左右n_bユニットとし，初期の電子温度はθ = 5 mHとした．Figure 1に通常のDC法とFT-DC法に対するエネルギー誤差のバッファサイズ依存性を示す．ここでエネルギー誤差はDC-BLYP法については整数占有数を用いたBLYP法との差，FT-DC-BLYP法についてはFT-BLYP法との差として定義した．C₄₀H₈₂は，単結合で構成されているためHOMO-LUMOギャップが大きい．HOMO-LUMOギャップが大きい場合は，FT法を用いた場合でも占有数がほぼ整数となり，DC法とFT-DC法は同様の結果を示す．一方，C₄₀H₄₂は二重結合と一重結合が交互に結合しておりHOMO-LUMOギャップが小さい．HOMO-LUMOギャップが小さい場合は，非整数占有数を持つ軌道が多くなり，通常の整数占有数の解とFT法の解は大きく異なる．FT-DC法の誤差は，C₄₀H₈₂の場合と同様に，バッファサイズの増大に対して速やかに減少していく．通常のDC法は，バッファサイズの増大に対してエネルギー誤差は減少する傾向を示すものの，数mH程度の誤差が残る．これは，通常のDC法においても式(4)に示す通り占有数を非整数として取り扱っているため，整数占有数の解ではなく非整数占有数解に収束しやすくなることに起因する．

Figure 1.

 Buffer-size n_b dependences of DC-BLYP and FT-DC-BLYP energy errors with respect to conventional BLYP energy and FT-BLYP energy, respectively.

同様の計算条件で，ポリアセンの一重項-三重項エネルギー差に対する検証を行った．部分系はC₄H₂またはC₆H₄，バッファサイズは左右4ユニットである．Figure 2に炭素数を増大させた際のポリアセンC_4n+2H_2n+4の一重項-三重項エネルギー差を示す．ポリアセンはπ共役の増大に従いHOMO-LUMOギャップが小さくなるため一重項-三重項エネルギー差は単調減少する傾向にある．通常のDFT計算ではn > 10において非物理的な挙動を示す．このようなHOMO-LUMOギャップが小さい場合，静的相関効果が大きく，通常のDFT計算では表現することが困難である．一方FT-DFT法では，静的相関効果が取り込まれることによって非物理的な挙動が改善されている．

Figure 2

 Singlet-triplet gap of polyacene C_4n+2H_2n+4 (n = 1-20).

4 まとめ

本稿ではFT法をDC-SCF法に拡張したFT-DC-SCF法について解説した．電子のエントロピー項の導入と電子温度の自動決定法を用いることにより，静的相関効果を簡易的に取り込むだけでなく，通常のDC-SCF法のエネルギー誤差を改善できる．また，FT-DC-SCF法は，静的相関効果が重要なポリアセンの一重項-三重項エネルギー差のような大規模π共役系を取り扱うことができる．

Acknowledgment

本研究で行った量子化学計算の一部は，自然科学研究機構(NINS)・計算科学研究センター(RCCS)の計算機を利用した．

参考文献

• [1]   P. J. Stephens, F. J. Devlin, C. F. Chabalowski, M. J. Frisch, J. Phys. Chem., 98, 11623 (1994). doi:10.1021/j100096a001
• [2]   Y. Zhao, D. G. Truhlar, Theor. Chem. Acc., 120, 215 (2008). doi:10.1007/s00214-007-0310-x
• [3]   H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys., 115, 3540 (2001). doi:10.1063/1.1383587
• [4]   T. Yanai, D. P. Tew, N. C. Handy, Chem. Phys. Lett., 393, 51 (2004). doi:10.1016/j.cplett.2004.06.011
• [5]   S. J. Grime, J. Comput. Chem., 25, 1463 (2004). , doi:10.1002/jcc.2007815224390
• [6]   J. D. Chai, M. Head-Gordon, Phys. Chem. Chem. Phys., 10, 6615 (2008). , doi:10.1039/b810189b18989472
• [7]   T. Sato, H. Nakai, J. Chem. Phys., 131, 224104 (2009). , doi:10.1063/1.326980220001021
• [8]   S. Kimmel, L. Krunik, Rev. Mod. Phys., 80, 3 (2008). doi:10.1103/RevModPhys.80.3
• [9]   D. Peng, S. N. Steinmann, H. van Aggelen, W. Yang, J. Chem. Phys., 139, 104112 (2013). , doi:10.1063/1.482055624050333
• [10]   N. D. Mermin, Phys. Rev., 137, 5A, A1441 (1965). doi:10.1103/PhysRev.137.A1441
• [11]   C. A. Bauer, A. Hansen, S. Grimme, Chemistry, 23, 6150 (2017). , doi:10.1002/chem.20160468227906486
• [12]   J. D. Chai, J. Chem. Phys., 136, 154104 (2012). , doi:10.1063/1.370389422519312
• [13]   W. Yang, Phys. Rev. Lett., 66, 1438 (1991). , doi:10.1103/PhysRevLett.66.143810043209
• [14]   W. Yang, T. S. Lee, J. Chem. Phys., 103, 5674 (1995). doi:10.1063/1.470549
• [15]   T. Akama, M. Kobayashi, H. Nakai, J. Comput. Chem., 28, 2003 (2007). , doi:10.1002/jcc.2070717455367
• [16]   M. Kobayashi, Y. Imamura, H. Nakai, J. Chem. Phys., 127, 074103 (2007). , doi:10.1063/1.276187817718602
• [17]   T. Doi, T. Yoshikawa, H. Nakai, in preparation.
• [18]   C. Y. Lin, K. Hui, J. H. Chung, J. D. Chai, RSC Advances, 7, 50496 (2017). doi:10.1039/C7RA10241K
• [19]   A. D. Becke, Physical Review A, 38, 3098 (1988). doi:10.1103/PhysRevA.38.3098
• [20]   C. Lee, W. Yang, R. G. Parr, Phys. Rev. B, 37, 785 (1988). doi:10.1103/PhysRevB.37.785
• [21]   W. J. Hehre, R. Ditchfield, J. A. Pople, J. Chem. Phys., 56, 2257 (1972). doi:10.1063/1.1677527
• [22]   J. D. Dill, J. A. Pople, J. Chem. Phys., 62, 2921 (1975). doi:10.1063/1.430801