Counterpoise補償法によるBSSEの見積
分子間の相互作用は、計算化学が扱う主要な問題の一つです。単純に考えると、例えばAという分子とBという分子が相互作用してABという錯体を形成する場合、A,Bが単独で存在したときのエネルギーと錯体ABのエネルギーを比較して、その相互作用がどれぐらい強いかを議論するということになります。しかし、ここには思いがけない落とし穴が存在します。それが本稿で紹介するBSSE(Basis Set Superposition Error:基底関数重なり誤差)です。以下、GAMESSでのBSSEの計算方法とモデル化学による差について見てみたいと思います。
What's BSSE ?
先に述べた[A + B → A-B]における計算をもう少し詳しく見ると、Aを単独で計算するのと、A-Bという錯体で計算するのに違いがあることがわかります。即ち、錯体A-Bの中のAは、Bの基底関数を取り込むことができます。普通はこの「取り込み」はイコール相互作用ということになりますが、小規模の基底関数(例えばSTO-3Gや3-21Gなど)では、Bの基底関数をAの基底関数の分極関数やDiffuse関数のように使われてしまい、結果、相互作用のエネルギーだけでなく基底関数の改良による安定化エネルギーが取り込まれた形で数値が出ます。この「基底関数の改良による安定化エネルギー」こそがBSSEです。どんなに大規模な基底関数でもこの寄与が(大幅に縮小するとはいえ)完全に取り除けないことは、容易に想像できると思います。ですのでBSSEを求めることは、相互作用を扱う場合必須であると言えます。
(※これから紹介するBSSEの計算(counterpoise補償法)は、BSSEと完全にイコールではありません。よって、"BSSEを補正した"と言っても完全に除けたわけではなく、その比率が目的とする相互作用のエネルギーに比べて十分小さくするぐらいの意味合いで考えるべきでしょう。時には補正しすぎることもあるようです。)
An Estimate of BSSE
メチルイソシアニド(CH3-NC)と三フッ化ホウ素(BF3)の錯体を例に、実際にBSSEを計算してみましょう。
本来、2分子間の相互作用はGibbs自由エネルギーで取り扱うべきですが、BSSEは全電子エネルギーに影響を与えるということと、計算時間の都合で振動計算は省略します(手抜きじゃないですよ？)。原系も生成系も対称性が高い構造をしていますので、メチルイソシアニドと錯体はC3vで、三フッ化ホウ素はD3hで計算を行うことにします。
それぞれの入力構造($DATAグループ)は以下のとおりです。非常に単純な分子なので、モデリングソフトを使うよりも手計算で組んだほうが早いでしょう。
(1)CH3-NC
(2)BF3
(3)CH3-NC…BF3 complex
これらについて、理論モデルとしてRHFとRPBE0(hybrid-DFT)を、基底関数として3-21Gと6-31G(d)を使い(つまり4種類の異なるモデル化学で)構造最適化を行ってみました。最適化のオプションとして、METHOD=RFO, OPTTOL=1.0E-05に設定しています($STATPT内)。それぞれ見積もられた相互作用エネルギーΔEは以下の通りです。
| Eel (kcal/mol) | |
| CH3-NC | -82303.62308 |
| BF3 | -201722.8716 |
| complex | -284046.4195 |
| ΔE | -19.92481976 |
|---|
| Eel (kcal/mol) | |
| CH3-NC | -82764.96653 |
| BF3 | -202808.2373 |
| complex | -285579.2782 |
| ΔE | -6.074389726 |
|---|
| Eel (kcal/mol) | |
| CH3-NC | -82718.39100 |
| BF3 | -202364.8063 |
| complex | -285109.6244 |
| ΔE | -26.42704571 |
|---|
| Eel (kcal/mol) | |
| CH3-NC | -83177.16892 |
| BF3 | -203461.1812 |
| complex | -286652.6213 |
| ΔE | -14.27117145 |
|---|
モデル化学によって、相互作用エネルギー(ΔE)の値は大きく異なります。3-21Gを使ったときに絶対値が大きくなる傾向にありますが、これはBSSEが多く含まれているためです。問題のBSSEですが、上記計算だけでは求められません。BoysのCounterpoise補償法では、あと4つエネルギー計算をする必要があります。
Counterpoise補償法では、「ゴースト原子」と呼ばれる「原子核がない、空の軌道だけがある原子」を利用してBSSEを見積もります(下図)。まず、最適化された錯体から一方のみをそっくり取り出し、単独でエネルギー計算(RUNTYP=ENERGY)を行います(E1)。次に、錯体の中で相手の方をゴースト化し、エネルギー計算を行い(E2)、両者の差を取ります(E1-E2)。つまり構造は全く同じで、相互作用する相手分子の軌道がある時とない時のエネルギー差が「BSSE」であるとするわけです。
RPBE0/6-31G(d)を例に説明すると、まず、CH3-NC…BF3 complexの最適化構造を抜き出します。対称性は引き続き利用しますので、「COORDINATES OF SYMMETRY UNIQUE ATOMS」の方を使います。
これをCH3-NCとBF3とに分けます。これら2つのエネルギー計算の和がE1の計算に相当します。
次に、CH3-NCとBF3とをそれぞれゴースト化します。陽子数(2列目)を負の値にすればOKです。これら2つのエネルギー計算の和がE2の計算に相当します。
あとは、E1-E2でBSSE(counterpoise)が算出されます。前出の4つのモデル化学についてBSSEと補正ΔEを示します。(corr.ΔE = ΔE + BSSE : BSSE補正相互作用エネルギー)
| Eel (kcal/mol) | |
| ΔE | -19.92481976 |
| BSSE | 15.41135853 |
| corr.ΔE | -4.513461225 |
|---|
| Eel (kcal/mol) | |
| ΔE | -6.074389726 |
| BSSE | 4.158496483 |
| corr.ΔE | -1.915893243 |
|---|
| Eel (kcal/mol) | |
| ΔE | -26.42704571 |
| BSSE | 17.46429935 |
| corr.ΔE | -8.962746364 |
|---|
| Eel (kcal/mol) | |
| ΔE | -14.27117145 |
| BSSE | 5.530825435 |
| corr.ΔE | -8.740346015 |
|---|
3-21G基底を使った場合、BSSEが非常に大きくなっていることが一目瞭然です。6-31G(d)基底ではBSSEが大幅に減少しています。RPBE0ではBSSE補正がきちっと効いているのか、補正後のエネルギーが3-21Gと6-31G(d)でほとんど変わりません。比較のため、基底関数をさらに拡張した6-311+G(2d)と、6-31G(d)と同じスプリットバレンスDZ基底でDFTに最適化されているDZVP(6-41G(d)ぐらいに相当 : Can. J. Chem. 1992, 70, 560. EMSL Basis Set Exchangeから入手)でも同様に計算した結果を示します。
| Eel (kcal/mol) | |
| ΔE | -12.41652696 |
| BSSE | 0.932262563 |
| corr.ΔE | -11.4842644 |
|---|
| Eel (kcal/mol) | |
| ΔE | -14.86625195 |
| BSSE | 2.460337491 |
| corr.ΔE | -12.40591446 |
|---|
RPBE0/6-311+G(2d)の結果は流石に優れています。BSSEは1 kcal/mol以内。DZVPの結果も次いで良く、しかも6-311+G(2d)よりも3倍以上速く計算することができます。6-31G(d)とDZVPは所要時間はあまり変わらない(6-31G(d)が15.1 sec, DZVPが18.1 sec : ゴースト化錯体のエネルギー計算で比較)のですが、BSSEはDZVPが半分、最終的な補正エネルギーも6-311+G(2d)に近いということで、本系に対して非常に有効な基底関数のようです。
Morokuma Docomposition for RHF-BSSE
RHF法に限り、上記のような複数の計算を行わなくてもCounterpoise補償法によるBSSEを計算することができます。最適化した錯体の構造を使い、$CONTRLグループでRUNTYP=MOROKUMAとし、$MOROKMグループに「BSSE=.T. IATM(1)=n」と記述します。また、対称性は利用せずにC1に、Direct SCFはできないので$SCFグループでDIRSCF=.F.にしています(デフォルトが.F.ですので、何も書かなくても構いません)。(下記例はRHF/3-21G)
IATM(1)=nの「n」は、モノマーの区切りを示すもので、上記のケースでは6番目までがメチルイソシアニドですので「6」としています。もし3成分以上の錯体だった場合は、各モノマーの最後の原子をカンマ区切りの番号で指定します(最後の1つは必要ありません)。
出力は一番最後で、以下のような内容です(抜粋)。
一番最後のBSSE CORRECTIONが目的の値です。先の手計算とも合致しています。
コメント
ヨウ素原子の6-31G(d)相当の基底をEMSLで探していて、DZVPを見つけ、DZVPをgoogleで検索してこのページにたどり着きました。DZVPが6-311+G比較してあってとても参考になりました。
さて、EMSLから入手したGamess-US用DZVP基底には、DATAカードが3つ入っているのですが、どれを使えば良いのでしょうか？
Sasaki様の"inp"ファイルを拝見しましたが、external fileになっていて、判断できませんでした。
GamessによるDFT計算は初めてなので、非常に初歩的な事を質問しているかもしれませんが、ご教授いただければ幸いです。
よろしくお願いいたします。
投稿者: 藤原洋規 | 2008年02月03日 22:22
inpファイルではなくoutファイルの方を見ると分かります。「ATOMIC BASIS SET」の下を見ると、基底関数のデータが出力されていますので、是非ご確認を。
そこを見ると分かりますが、、私は縮約されている基底関数系(EMSLで言う"一番上"に出てくるデータ)を使っています。これが一番低いエネルギーを与えますし、縮約されている分計算も速いです。EMSLのサイトで「Get Basis Set」ボタンの下にある「more information...」から、その基底関数に関する情報を見ることができますが、そこでも基底関数の縮約が示されていますので、それを根拠に使っています。
実はお恥ずかしながら、なぜ3種類登録されているかは私は知りません。もともとDGaussで使われていた基底系で、そのときに非縮約型が用いられていたのでしょうか？誰かご存知の方にご教授頂ければと思います。
投稿者: s2k | 2008年02月03日 23:09
ご回答ありがとうございました。EMSLで配布されている3つの$DATAカードの内、1、2番目は、引用元論文(N. Godbout et al. Can. J. Chem. 70, 560 (1992))にありました。1番目は"Coefficients and exponetns for (63/5/1*) orbital basis sets"、2番目は"Exponetns for the (4,3;4,3) auxiliary sets for use with the (63/5) orbital set"というタイトルが付けられています。本文中に、2番目のカードは"elextron density"を近似するためのauxiliary basis setであると書かれています。
3番目のカードは"DGauss basis set library"からの引用したものであると思われます。Gaussian03のDGDZVP基底に関する文献(C. Sosa et al. J. Phys. Chem. 96 6630 (1992))によるとauxiliary basis setは"elextron density"を近似するためのものと"exchange-correlation potential and energy"を近似するためのものがあるそうです。
"In the analytical approach used in DGauss, the total electron
density is variationally fit to a set of auxiliary Gaussian-type basis
functions leading to exact Coulomb forces. The remaining
exchange-correlation energy term is a smooth function of the
density and can be accurately fit to another auxiliary set of
Gaussian-type functions using an adaptive set of grid points
as described below."
したがって3番目のカードは、"exchange-correlation potential and energy"を近似するためのauxiliary basis setだと思われます。
釈迦に説法かもしれませんが、auxiliary basis setは積分を近似計算するための基底関数だそうです。例えば"elextron density"をauxiliary basis setの線形結合で近似すると、電子同士の反発に対する静電力の寄与を表す積分の計算コストが、基底数の4乗に比例から基底数の自乗とauxiliary basis setの数の積に比例に軽減するそうです。(Wolfram Koch et al. A Chemist's Guide to Density Functional Theory, 2nd Ed., (2001) pp.102-103)
投稿者: 藤原洋規 | 2008年02月04日 13:45
>藤原さん
ref情報込みでの詳細なコメント、有難うございます。JPCについてはすぐに文献が見れるので、目を通してみます。
投稿者: s2k | 2008年02月12日 19:25