Development of Bond Energy Density Analysis with Informatics Technique

Abstract

The bond energy density analysis (bond-EDA), which is one of the two-body energy decomposition schemes to obtain bond energies, was extended by revisiting the formula and the constraint conditions, and using the informatics technique. Numerical assessments of the present schemes were performed for 44 chemical bonds in 33 small molecules with covalent and ionic bonds involving second and third row atoms. The results show that the present scheme qualitatively reproduces the bond dissociation energies obtained by the direct evaluation scheme in quantum chemical calculations.

1 はじめに

分子や分子集合体が示す構造，反応性，物性には，構成する原子間の結合の強さが密接に関係している．また，その様式は多種多様で，共有結合性やイオン結合性を有した強い力からファンデルワールス力など分子間に働く弱い力がある．結合の強さに対応する観測量は結合解離エネルギー(BDE)であるが，構成するすべての原子間の結合を切断し観測することは小分子を除いて現実的でない．また，環状化合物などのように部分系に分割するには複数の化学結合を同時に切断する必要がある分子も存在し，その場合，BDEの見積りは困難となる．

我々の研究室では，分子や分子集合体に対する電子状態計算の結果を利用して，全ての原子間の結合の強さを同時に見積る解析手法として，結合エネルギー密度解析(Bond-EDA) [1]を提案してきた．同様の試みは，Ichikawaら [2]やMayer [3]など他のグループ [4, 5]でも行われてきたが，多重結合の記述が困難であるなど，定性的にも満足のいく結果は得られなかった．2011年に我々が発表した改良では，結合エネルギー(BE)の微分に関する条件を課すことで，これまでの困難を克服できることに成功した [6]．一方で，弱い相互作用の部位を同時に扱うと，数値的不安定性が生じるという欠点も明らかとなった．最近，我々は機械学習の1つであるLASSO回帰を用いることでこの数値的不安定性の問題を回避できることを示した [7]．本稿では，このLASSO回帰を用いた解析手法とさらに非対称的な結合様式への拡張に関する予備的考察について解説する．

2 結合エネルギー(BE)の解析方法

分子内の化学結合A−Bが切断されることによって，分子が二つの異なるフラグメントX, Yに解離するとき，BDEは一般に以下の式で表される．   

$E_{\text{BDE}}^{AB}=E^{\text{tot}}-E^X-E^Y$(1)

ここでE^tot，E^X，E^Yはそれぞれ分子全体，フラグメントX，Yのエネルギーを示す．実験的に観測する場合，解離過程における構造緩和や振動モード変化なども含まれる．一方，理論的にはそのような効果を分けて解析することができる．しかし，上述の通り，環状化合物などのように複数の結合を同時に切断する必要がある場合，BDEを求めることができない．

我々の研究室では，分子を構成する原子間のBE (非観測量)を見積る解析手法，Bond-EDAを提案してきた．Bond-EDAでは，次式のように全エネルギーE^totは原子内一体成分E^AAおよび原子間二体成分E^ABに分割される．   

$E^{\text{tot}}=\sum _A{E}^{AA}+\sum _A\sum _{B\left(\ne A\right)}{E}^{AB}$(2)

E^ABの反対符号(−E^AB)はBEに対応する．さらに適切なBEの総和を取ることで以下のようにBDEを算出できる．   

$E_{\text{BDE}}^{AB}=\left(-E^{AB}\right)+\sum _{C\left(\ne A\right),C\in X}\left(-E^{CB}\right)+\sum _{D\left(\ne B\right),D\in Y}\left(-E^{AD}\right)$(3)

ここで，式(3)の二体成分は，第一項の該当する原子間だけでなく，フラグメントX, Y間のすべての相互作用を考慮している．この式は，BDEを推定できない原子間の結合エネルギーも評価することができる．

式(2)の一体成分および二体成分は，一般的に以下のように算出される．   

$\begin{array}{c}{E}^{AA}=T_{\text{S}}^{AA}+E_{Ne}^{AA}+E_{ee}^{AA}+E_{HFx}^{AA}+E_{DFTxc}^{AA}\\ +\frac{1}{2}\sum _{B\left(\ne A\right)}{\alpha }^{AB}\left(T_{\text{S}}^{AB}+E_{HFx}^{AB}+E_{DFTxc}^{AB}\right)\end{array}$(4)

  

$E^{AB}=E_{Ne}^{AB}+E_{ee}^{AB}+E_{NN}^{AB}+\left(1-{\alpha }^{AB}\right)\left(T_{\text{S}}^{AB}+E_{HFx}^{AB}+E_{DFTxc}^{AB}\right)$(5)

ここで，T_S, E_Ne, E_ee, E_NN, E_HFx, E_DFTxcはそれぞれ運動エネルギー，核―電子引力，電子―電子斥力，核―核斥力，HF交換エネルギー，DFT交換・相関エネルギー成分である．Mayerが提案した方法(M03 [3])はα^AB = 0に対応するが，一般的にBEを過小評価する傾向がある．オリジナルのBond-EDA (NK05 [1])では，(4)-(5)式の右辺括弧内の項，いわゆる量子力学的な項を一律に50%原子間から原子に分配することで，過小評価の低減に成功した．すなわち，α^AB = 0.5に対応する．NOIK11 [6]では，多重結合性を再現するため，平衡構造に対する拘束条件を用いることで，原子間ごとに異なるパラメータα^AB を決定することとした．   

$\frac{\partial E^{AB}}{\partial R^{AB}}=0$(6)

平衡構造に限らない場合，エネルギー勾配を用いる以下の条件によりα^AB を決定できる．   

$\sum _{C,D\left(>C\right)}\frac{\partial E^{CD}}{\partial R^{AB}}=\frac{\partial E^{\text{tot}}}{\partial R^{AB}}$(7)

この条件は，一体のエネルギーは結合長の微小な変化に対して不変と仮定する [5]．条件式(7)からなる連立方程式の数は分割表式(4)-(5)に現れるパラメータα^AB の数(N*(N−1)/2個，Nは原子数)と等しいため，一意にα^AB が決定できる．しかし，すべての原子間を考慮する場合，ほとんど相互作用がない原子間も存在するため，連立方程式の係数行列が非常に小さな値となる．その結果，数値的な不安定性が生じ，連立方程式の解であるパラメータα^AB が非物理的な挙動を示す．そこで，インフォマティクス手法の一つであるLASSO回帰を用いることで，問題回避する方法を提案した [7]．

本稿では，さらに一般的な分割表式について検討する．   

$\begin{array}{c}{E}^{AA}=T_{\text{S}}^{AA}+E_{Ne}^{AA}+E_{ee}^{AA}+E_{HFx}^{AA}+E_{DFTxc}^{AA}\\ +\sum _{B\left(\ne A\right)}{\alpha }^{AB,AA}\left(T_{\text{S}}^{AB}+E_{HFx}^{AB}+E_{DFTxc}^{AB}\right)\end{array}$(8)

  

$\begin{array}{c}{E}^{AB}=E_{Ne}^{AB}+E_{ee}^{AB}+E_{NN}^{AB}\\ +\left(1-{\alpha }^{AB,AA}-{\alpha }^{AB,BB}\right)\left(T_{\text{S}}^{AB}+E_{HFx}^{AB}+E_{DFTxc}^{AB}\right)\end{array}$(9)

ここで，パラメータα^AB,AAとα^AB,BBを用いることにより，原子間ABに分配された量子力学的な項を原子AとBに非対称に分配することができる．イオン結合性の大きい場合などに有効であると思われる．しかし，(4)-(5)式のパラメータα^ABの個数がN*(N−1)/2個であるのに対し，(8)-(9)式のパラメータα^AB,AAの個数は2倍のN*(N−1)個となる．条件式 (7) N*(N−1)/2個の連立方程式であるため，パラメータα^AB,AAは一意には決定できない．そこで，以下の条件式もあらわに用いることとした．   

$E^{AA}=E_0\left(A\right)$(10)

ここで，E₀(A)は自由原子Aのエネルギーを示す．(10)式はN個しか条件を与えないので，まだパラメータ数に比べて条件式が足りず，解の任意性が残っている．しかし，LASSO回帰による次元数削減を用いることで，パラメータα^AB,AAを決定することができると予想される．

3 計算条件

すべての計算はGAMESS [8]を用いて行い，計算レベルはB3LYP/6-311G (d,p)とした．BEの計算手法として分割表式(4)-(5)に条件式(7)を適用した方法(Model #1)，分割表式(8)-(9)に条件式(7)と(10)を適用した方法(Model #2)を，比較としてM03 [3]，NK05 [1]，NOIK11 [7]を用いた．Model #1, #2における連立方程式は独自に実装したLASSO回帰を用いて解いた．LASSO回帰におけるハイパーパラメータは，後に示す44結合に関する統計誤差が最小となる0.005を採用した．各解析手法におけるBDEは式(3)の定義にしたがって計算した．また量子化学計算によるBDEの参照値(QC)は構造緩和を考慮しない式(1)の定義に基づき算出した．

4 結果と考察

Table 1に第二周期までの原子を含む3原子分子の，共有結合(原子間を「−」で示す)およびthrough-space/bond相互作用(「…」で示す)に対する，各手法によるパラメータα^ABおよびα^AB,AAの値を示す．ここで，直接的な化学結合は存在しないが，クーロン相互作用など空間を介した長距離相互作用をthrough-space，いくつかの化学結合を介した間接的な相互作用をthrough-bondと呼んでいる．Model #2では各結合に対してそれぞれの原子に非対称に分配するため，α^AB,AAとα^AB,BBの和はNOIK11とModel #1のα^ABに対応する．NOIK11では，through-space/bond相互作用のαが1以上となり非物理的な挙動を示す．一方，Model #1, #2ではすべて0∼1の範囲になり，LASSO回帰により非物理的な挙動が回避されたことがわかる．Model #2では各結合に対してそれぞれの原子に非対称に分配され，α^AB,AAとα^AB,BBの和はModel #1のα^ABに近い値となっている．さらにα^AB,AAとα^AB,BBの大小と，それぞれの原子のもつ電気陰性度の大小に相関があることもわかる．たとえばHCNにおいて，C_A−H_B結合のパラメータ値は{α^AB,AA, α^AB,BB} = {0.50, 0.22}，C_A−N_B結合では{α^AB,AA, α^AB,BB} = {0.09, 0.41}となり，H < C < Nの順となり，それぞれの原子の電気陰性度の大小に対応している． このようにModel #2は，パラメータ数に比べて条件式が少ないにもかかわらず，すべてのパラメータが有意な挙動をしており，これはLASSOの次元削減が寄与していると考えられる．

Table 1. Parameter set {α^AB } in triatomic molecules by NOIK11, Models #1 and #2.

Molecule	Bond	NOIK11	Model #1	Model #2
		αAB	αAB	αAB,AA	αAB,BB
CH2	CA−HB	0.75	0.72	0.50	0.22
	HA…HB	1.54	1.00	0.50	0.50
NH2	NA−HB	0.58	0.59	0.50	0.08
	HA…HB	2.95	1.00	0.50	0.50
HCN	CA−HB	0.85	0.73	0.50	0.22
	CA−NB	0.55	0.51	0.09	0.41
	HA…NB	1.12	1.00	0.50	0.50
H2O	OA−HB	0.60	0.62	0.50	0.13
	HA…HB	−9.14	1.00	0.50	0.50
HCO	CA−HB	0.89	0.77	0.50	0.26
	CA−OB	0.67	0.62	0.00	0.64
	HA…OB	1.77	1.00	0.50	0.50
CO2	CA−OB	0.63	0.56	0.03	0.52
	OA…OB	1.71	1.00	0.50	0.50

Figure 1に，第三周期までの原子を含む，33個の典型的な小分子の44結合 (共有結合とイオン結合を含む) に対する，各手法によるBEから式(3)により算出したBDEの結果を示す．ここでMAD、MaxDおよびR²はそれぞれ，QCからの平均絶対偏差，最大偏差，y切片を0としたときの決定係数である．Model #2ではMADが10.6kcal/mol，R²が0.9051と最も良い相関を示したこれにより，本モデルの妥協性とLASSO回帰による過適合の解消および次元数削減の有効性を示している．

Figure 1.

 Correlation between the BDEs (kcal/mol) for 44 chemical bonds in 33 covalent and ionic molecules obtained by the quantum chemical calculations and those by the energy decomposition schemes.

5 まとめ

本研究では，bond-EDAにおける分割表式と最適化条件の拡張を行った．33分子44結合に対する数値検証から，本手法は様々な化学結合に対するBDEを定性的に正しく再現し，参照値からの誤差が小さくなることも確認された．

参考文献

• [1]   H. Nakai, Y. Kikuchi, J. Theor. Comput. Chem., 04, 317 (2005). doi:10.1142/S0219633605001465
• [2]   H. Ichikawa, A. Yoshida, Int. J. Quantum Chem., 71, 35 (1999).
• [3]   I. Mayer, Chem. Phys. Lett., 382, 265 (2003). doi:10.1016/j.cplett.2003.10.097
• [4]   H. Sato, S. Sakaki, J. Phys. Chem. B, 110, 12714 (2006). , doi:10.1021/jp061017s16800606
• [5]   S. F. Vyboishchikov, Int. J. Quantum Chem., 108, 708 (2008). doi:10.1002/qua.21545
• [6]   H. Nakai, H. Ohashi, Y. Imamura, Y. Kikuchi, J. Chem. Phys., 135, 124105 (2011). doi:10.1063/1.3636387
• [7]   H. Nakai, J. Seino, K. Nakamura, J. Phys. Chem., in press.
• [8]   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. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem., 14, 1347 (1993). doi:10.1002/jcc.540141112