多成分密度汎関数理論のための電子-核相関汎関数の開発

Abstract

最も軽い水素の原子核の運動においては，量子効果が顕著に現れる. 原子核運動の量子効果を直接的に取り込んだ理論として多成分分子理論を開発してきた．また近年では，計算時間のかからない密度汎関数理論が分子理論の中心理論となりつつある．我々は，多成分密度汎関数理論のための簡便な電子-核相関汎関数を開発した. この汎関数は，相関波動関数から簡単かつ正当な物理条件のみにもとづいて導出した汎用性の高い汎関数であり，パラメータを1つしか含まない．にもかかわらず，水素原子を含む小分子の電子-プロトン相関エネルギーを化学的精度(1∼2 kcal/mol程度の誤差)で再現する．

1 研究目的

水素原子は多くの分子系に存在し，水素結合や水素移動を介して分子の構造や物性に決定的に重要な役割を果たす．また，近年の実験技術，理論計算手法の進展により，プロトントンネリングのような水素の原子核の持つ量子効果が顕著な寄与を示す現象の重要性が確認され始めている．しかし，量子化学計算手法では，電子に比べて圧倒的に質量の重い原子核の運動はあらわに考慮せず，固定された原子核の作る力場の中の電子の運動状態のみを取り扱うBorn-Oppenheimer (BO)近似に基づくことがほとんどである．密度汎関数理論(DFT) [1,2]は，少ない計算時間で小分子から大規模系まで高精度な構造や物性を与える手法として量子化学計算手法の中心理論になりつつあるばかりでなく，原子核物理学や物性科学•生体物質科学といったきわめて幅広い分野の応用計算で成功を収めている．しかし，DFTも一般的にBO近似に基づいており，プロトントンネリングのような原子核の量子効果が重要な役割を果たす現象は取り扱えなかった．

電子のみならずプロトンのような軽い原子核自身の量子効果を直接考慮可能な第一原理計算手法として多成分分子軌道(Multicomponent molecular orbital: MC_MO)法 [3]がある．この理論は，様々な分子や反応に対するH/D同位体効果の解析に利用されてきた．しかし，これまでのMC_MO法は定量性が十分ではなく，さらなる適用性の向上には理論の効率のよい高精度化が必要であった．多成分分子理論の高精度化には，電子間の相関効果のみならず，電子と原子核との間の相関効果(電子-核相関)をも取り込む必要がある．電子-核相関はこれまで，摂動論 [4]や配置間相互作用法 [5]を用いて評価されてきたが，計算時間がかかりすぎること，収束解が得られないことといった問題を抱えていた．この問題を解決するため，多体効果をDFTで評価する多成分DFT (MC_DFT)に着目した．

MC_DFTは，CapitaniとParrらにより1982年に証明された多成分系に対するHohenberg-Kohnの定理に基づく理論 [6]であり，これまで多くの研究がなされてきた [7,8,9]．我々も，電子相関を考慮したMC_DFTにより，多成分Hartree-Fock法では定性的にすら正しく再現できないポルフィセンおよびポルフィリン分子におけるH/D同位体効果を正しく評価できることを示した [10]．しかし，その結果はまだ定量的と言えるレベルではなかった．定量性を確保するためには電子-核相関を取り込む必要があったが，実用的な電子-核相関汎関数は存在していなかった．

本研究では，最も物理的裏付けのある電子相関汎関数の1つの型でありLee-Yang-Parr (LYP)汎関数やOne-parameter Progressive (OP)汎関数の基となって広く利用されているColle-Salvetti (CS)型相関汎関数 [11]を電子-核相関へと拡張し，パラメータを一つしか含まない汎用的な電子-核相関汎関数を開発することを目的としている [12]．

2 理論

電子-核相関のためのCS型相関波動関数は，電子-核間に相関のない波動関数，すなわち電子どうしの相関効果のみを考慮したKohn-Sham (KS)波動関数(Ψ_KS)にCS型の電子-核相関カスプ(相関領域の核近傍の尖頭)を表現した因子を乗じた形をもち，以下のように定義される．   

${\Psi }_{\text{CS}}^{\text{en}}={\Psi }_{\text{KS}}\prod _{i,n}\left[1-{\varphi }^{\text{en}}\left(r_i, r_n\right)\right]$(1)

  

${\varphi }^{\text{en}}\left(r_i, r_n\right)=\text{exp}\left(-{\beta }_{\text{en}}^2{r}^2\right)\left\{1-{\text{Φ}}^{\text{en}}\left(R\right)\left[1+{\eta }_{\text{cusp}}^{\text{en}}r\right]\right\}$(2)

ここで，ϕ^enは電子-核相関のためのCS型相関カスプ因子であり，r_iとr_nはそれぞれ電子と核の座標ベクトルである．rとRはそれぞれ，$r=\left|r_i-r_n\right|$, $R=\alpha r_i+ \left(1-\alpha \right)r_n\approx r_n$と定義される．式(2)中の${\eta }_{\text{cusp}}^{\text{en}}$は，電子-核相関カスプ条件から決定される定数であり，${\eta }_{\text{cusp}}^{\text{en}}=-Z_n$(Z_nは核電荷)である．また，式(2)のΦ^en(R)については，電子-核に対する相関因子(式(2))を用いて，電子相関の場合と同様の手順で求めることができる．本研究では，Φ^en(R)の二次方程式を近似せずに解いた．

このCS型電子-核相関波動関数は，CS型電子-電子相関波動関数が満たすべき2つの物理条件を電子-核相関の場合に拡張した，以下の2つの物理条件を満たす．

(i) 電子と核が十分離れた場合には，電子-核に相関のない波動関数，つまり通常のKS波動関数に漸近する．

(ii) 電子と核が近接する場合は，電子-核相関カスプ条件を満足する．

(i)の条件は，CS型相関因子に含まれる指数因子により達成され，(ii)の条件は，CS相関因子中に含まれる${\eta }_{\text{cusp}}^{\text{en}}$により達成される．

求められた式(2)の電子-核CS相関因子の適用性を評価するため，超高精度計算である多成分量子モンテカルロ(MC_QMC)法 [13]により得られる電子-核相関Jastrow因子との比較検証を行った．MC_QMC法では，本手法と同様に非相関波動関数に相関因子(Jastrow因子)を掛けた形で相関波動関数が定義されるため，単純にMC_QMC法のJastrow因子とCS型相関カスプ因子を比較することで，正当性を評価することが可能である．Figure 1に，水素原子に対してMC_QMC法により最適化した電子-核Jastrow因子と，それにフィッティングしたCS型電子-核相関カスプ因子を比較した．Figure 1 (a)より，電子-核距離(r_ep)が約1.0 Åよりも大きくなると，式(2)のCS型相関カスプ因子ではJastrow因子を再現できなくなることがわかる．そこで，相関カスプ因子の有効領域を広げるため，式(2)を3次のTaylor展開まで以下のように拡張した．   

$\begin{array}{c}{\varphi }^{en}\left(r_i, r_n\right)=\text{exp}\left(-{\beta }_{en}^2{r}^2\right)\\ \times \left\{1-{\text{Φ}}^{en}\left(R\right)\left[1+{\eta }_{cusp}^{en}r+\frac{1}{2}{\eta }_{cusp}^{en}{}^2{r}^2+\frac{1}{6}{\eta }_{cusp}^{en}{}^3{r}^3\right]\right\}\end{array}$(3)

Figure 1.

 (a) The electron-nucleus correlation factor given by MC_QMC calculations and the original CS correlation factor for H atom, and (b) that given by MC_QMC calculations and new CS correlation factor (eq. (3)) for H atom [12].

その結果，Figure 1 (b)に示したように，拡張した電子-核CS相関因子(式(3))は式(2)のCS相関因子に比べ，より遠い電子-核間距離においても，MC_QMC法により得られた電子-核Jastrow因子を再現することに成功した．

この電子-核相関波動関数に基づく応用計算を行うには，相関カスプ因子に含まれるβ_enを決定する必要がある．β_enは相関領域の広さを決定する変数であり，適切な物理モデルに基づいて構築される．電子どうしの相関に対してはいくつかのモデルが提案されている．ColleとSalvettiは，相関領域の体積がWignerの排除体積 [14,15]に比例すると仮定し，${\beta }_{\text{ee}}=q{\rho }_{\text{e}}^{\frac{1}{3}}$ (qはパラメータ，ρ_eは電子密度)と求めた [11]．一方，常田らのOP相関汎関数では，Beckeの相関距離の定義を利用し，電子交換孔の体積が電子相関孔の体積に比例すると仮定し， ${\beta }_{\text{ee}}=q\frac{{\rho }_{\alpha }^{\frac{1}{3}}{\rho }_{\beta }^{\frac{1}{3}}{K}_{\alpha }{K}_{\beta }}{{\rho }_{\alpha }^{\frac{1}{3}}{K}_{\alpha }+{\rho }_{\beta }^{\frac{1}{3}}{K}_{\beta }}$と求めている [16]．ここで，K_σは電子交換エネルギー汎関数($E_{\text{x}}=-\frac{1}{2}\sum _{\sigma }{\rho }_{\sigma }^{\frac{4}{3}}{K}_{\sigma }{d}^{3}r$)により定義される．結果的に，Colle-Salvetti相関汎関数は常田らのOP相関汎関数の局所密度近似に相当する．

電子どうしの相関と異なり，電子-核相関と電子交換との間に直接的な関係性はないため，上記のβ_eeの関係式をそのまま適用することはできない．我々は，MC_QMC法によるテスト計算により，水素アニオンH⁻の電子-核相関領域の広さが，電子-核相関のみを考慮した場合と，電子-電子相関と電子-核相関の両者を考慮した場合とで著しく異なること，すなわち両者に関係性があることを見出した [12]．この結果を元に，電子-核相関領域の体積が電子-電子相関領域の体積に比例すると仮定した．この仮定に基づくと，電子-核相関汎関数におけるβ_enは   

${\beta }_{\text{en}}=q_{\text{en}}\frac{{\rho }_{\alpha }^{\frac{1}{3}}{\rho }_{\beta }^{\frac{1}{3}}{K}_{\alpha }{K}_{\beta }}{{\rho }_{\alpha }^{\frac{1}{3}}{K}_{\alpha }+{\rho }_{\beta }^{\frac{1}{3}}{K}_{\beta }}$(4)

と与えられる．β_eeと同様にこのβ_enも電子交換汎関数に   

$\begin{array}{c}{E}_{\text{C}}^{\text{en}}=-4\pi \int d_{}^3{r}_n{\rho }_e{\rho }_n\{-\frac{3}{4{\beta }_{\text{en}}^3}+\frac{32{\beta }_{\text{en}}^3+\left(8\sqrt{2\pi }-32\sqrt{\pi }\right)Z_n{\beta }_{\text{en}}^2+24Z_n^2{\beta }_{\text{en}}+\left(\sqrt{2\pi }-8\sqrt{\pi }{Z}_n^3\right)}{64{\beta }_{\text{en}}^5}{\Phi }_{}^{\text{en}}\left(r_n\right)\\ +\frac{384{\beta }_{\text{en}}^6-192\sqrt{2\pi }{Z}_n{\beta }_{\text{en}}^5+384Z_n^2{\beta }_{\text{en}}^4-96\sqrt{2\pi }{Z}_n^3{\beta }_{\text{en}}^3+112Z_n^4{\beta }_{\text{en}}^2-115\sqrt{2\pi }{Z}_n^5{\beta }_{\text{en}}^{}+8Z_n^6}{1536{\beta }_{\text{en}}^8}{\Phi }_{}^{\text{en}}{\left(r_n\right)}_{}^2\}\end{array}$(5)

より定義されるK_σを含んでいるが，それは電子-電子相関を介して電子-核相関領域と電子交換領域にも比例関係が与えられるからである．これらのCS型電子-核相関波動関数，相関因子，β_enを用いてOP相関汎関数 [17]に倣って電子-核相関エネルギー汎関数を導出すると，上の定式(式(5))が与えられる．ここで，ρ_nは核密度である．

電子-交換汎関数において唯一の未確定のパラメータは，相関領域の大きさを決めるq_enであり，これは用いる電子交換汎関数に依存する．このパラメータを決めるため，水素の原子核のみを量子力学的に取り扱い，電子-プロトン相関のみを考慮して全エネルギーのフィッティングを行った．BOP (Becke電子交換+OP電子相関)汎関数を用いたMC_DFTエネルギー(E^MC_BOP)と式(5)の電子-プロトン相関エネルギー($E_{\text{c}}^{\text{en}}$)の和が，BOP汎関数を用いた従来のDFTエネルギーと零点振動エネルギー補正の和(E^BOP + E^ZPE)と合うようにq_enをフィッティングした．計算系には，中性およびアニオンとカチオンの両方の環境を考慮するため，H₂, LiH, HeH⁺分子を選んだ．その結果，電子-プロトン相関汎関数に対する最適な値としてq_en = 18.0を得た [12]．

3 結果

開発した式(5)の電子-核相関エネルギー汎関数を用い，水素原子を含む18種類の小分子の電子-プロトン相関エネルギーを計算した．構造はDFT (BOP)/6-31G法により最適化した．電子の基底関数には6-31G基底を，量子的に取り扱った水素の原子核の基底関数にはs, p, d型のガウス型関数を1つずつ用いた．Table 1に零点振動エネルギー補正したDFT (BOP)のエネルギー(E^BOP + E^ZPE) [Hartree]と，BOP汎関数を用いたMC_DFT (MC_BOP)と式(5)により評価した電子-プロトン相関エネルギーの和 (E^MC_BOP + $E_{\text{c}}^{\text{en}}$) [Hartree]，およびそれらの差(E^Dev) [mHartree]をまとめた．その結果 18種類の分子に対する電子-プロトン相関エネルギー汎関数の平均絶対誤差は2.8 mHartreeであることが分かった．これは，開発した電子-核相関汎関数によって，きわめて精度良く分子の電子-プロトン相関エネルギーを評価できていることを示している [12]．また，それぞれの分子の誤差に着目すると，水素原子のアニオン性もしくはカチオン性が強いほど，誤差が大きい傾向があることがわかった．

Table 1.  The sum of the Kohn-Sham BOP energies and the corresponding zero-point vibrational energies (E^BOP + E^ZPE) [Hartree], sum of MC_BOP and electron-proton correlation energies (E^MC_BOP + $E_{\text{c}}^{\text{en}}$) [Hartree], and deviations (E^Dev)[mHartree] for 18 small hydrogen-containing molecules [12].

	H2							HeH+
EBOP + EZPE	–1.16148							–2.96621
EMC_BOP + $E_{\text{c}}^{\text{en}}$	–1.16130							–2.97605
EDev	0.1							–9.8
	LiH	BeH2	BH3	CH4	NH3	H2O	FH	NeH+
EBOP + EZPE	–8.07197	–15.88387	–26.55590	–40.43552	–56.47276	–76.34735	–100.37958	–128.95422
EMC_BOP + $E_{\text{c}}^{\text{en}}$	–8.06673	–15.87819	–26.55287	–40.44337	–56.47904	–76.35237	–100.38158	–128.95672
EDev	5.2	2.8	1.0	–2.0	–2.1	–2.5	–2.0	–2.5
	NaH	MgH2	AlH3	SiH4	PH3	SH2	HCl	ArH+
EBOP + EZPE	–162.82428	–201.19665	–244.12415	–291.76988	–343.03396	–399.29519	–460.72476	–527.58919
EMC_BOP + $E_{\text{c}}^{\text{en}}$	–162.81851	–201.18783	–244.11501	–291.76319	–343.03200	–399.29530	–460.72569	–527.59295
EDev	5.8	4.4	3.0	1.7	0.7	–0.1	–0.9	–3.8

本研究では，電子-プロトン相関エネルギーについてのみ実例を挙げて報告した [12]が，本研究で開発した電子-核相関汎関数は汎用的な汎関数であり，各原子に合うようなパラメータq_enを考えることで，同位体まで含めた種々の電子-原子核相関エネルギーを評価することができると考えられる．原子核の重さに依存するような形で適切にβ_enを構築できれば，原子種ごとにパラメータをフィッティングすることなく，すべての原子核を網羅した汎用的な電子-核相関汎関数の開発も可能になるであろう．

Acknowledgment

本研究を行うにあたり，多成分量子モンテカルロ法による計算をご協力いただきました北幸海准教授(横浜市大)に御礼申し上げます．

参考文献

• [1]   P. Hohenberg, W. Kohn, Phys. Rev., 136, B864 (1964).
• [2]   W. Kohn, L. J. Sham, Phys. Rev., 140, A1133 (1965).
• [3]   M. Tachikawa, K. Mori, H. Nakai, K. Iguchi, Chem. Phys. Lett., 290, 437 (1998).
• [4]   M. Tachikawa, H. Sainowo, K. Iguchi, K. Suzuki, J. Chem. Phys., 101, 5925 (1994).
• [5]   M. Tachikawa, Chem. Phys. Lett., 360, 494 (2002).
• [6]   J. F. Capitani, R. F. Nalewajski, R. G. Parr, J. Chem. Phys., 76, 568 (1982).
• [7]   Y. Imamura, H. Kiryu, H. Nakai, J. Comput. Chem., 29, 735 (2008).
• [8]   T. Kreibich, R. van Leeuwen, E. K. U. Gross, Phys. Rev. A, 78, 022501 (2008).
• [9]   A. Sirjoosingh, M. V. Pak, S. Hammes-Schiffer, J. Chem. Phys., 136, 174114 (2012).
• [10]   T. Udagawa, M. Tachikawa, J. Chem. Phys., 125, 244105 (2006).
• [11]   R. Colle, O. Salvetti, Theor. Chim. Acta, 37, 329 (1975).
• [12]   T. Udagawa, T. Tsuneda, M. Tachikawa, Phys. Rev. A, 89, 052519 (2014).
• [13]   Y. Kita, R. Maezono, M. Tachikawa, M. Towler, R. J. Needs, J. Chem. Phys., 131, 134310 (2009).
• [14]   E. Wigner, F. Seitz, Phys. Rev., 43, 804 (1933).
• [15]   E. Wigner, Phys. Rev., 46, 1002 (1934).
• [16]   T. Tsuneda, K. Hirao, Chem. Phys. Lett., 268, 510 (1997).
• [17]   T. Tsuneda, T. Suzumura, K. Hirao, J. Chem. Phys., 110, 10664 (1999).