水素結合は、分子間または分子内の水素を介して形成される弱い相互作用です。水素結合の強さは1~7kcal/mol程度と共有結合やイオン結合よりはるかに弱いです。しかし、水素結合は物性や化学的性質に大きな影響を与えていることは有名です(沸点、溶解度、立体異性、酸解離度etc.)。
hyd00
DNAが二重らせん構造を形成する際にも水素結合が重要な役割を果たしています。例えば、DNAの中でグアニンとシトシンそしてアデニンとチミン同士が互いに水素結合を介して分子認識し合う原因も分子軌道法の計算できちんと説明できます。今回は、水の二量体について分子間に働く相互作用を計算してみましょう。基本的な方法がマスターできれば、分子間相互作用に由来する様々な物性を説明できるようになります。
 

分子間相互作用の計算方法は？

分子間の相互作用エネルギーは、複合体モデルのエネルギー値から、各分子単体(孤立系)のエネルギー値を引いたものです。つまり、複合体モデルでは相互作用に因る安定化エネルギーを含みますが、孤立系のエネルギーはこれを含まないので、単純にこのエネルギー差が相互作用による安定化エネルギーに相当すると言う訳です。
bs3
ただし、前回説明したように弱い相互作用では特にBSSE(基底関数重なり誤差)を排除する必要があります。ここでは、水の2量体を複合体モデルとして分子間相互作用(主に水素結合)に因り分子がどの程度安定化エネルギーを得ているかを計算してみましょう。BSSEの計算は、Counterpoise法(CP法)で求めてみましょう。

BSSEを求めよう！！

前回の説明で、CP法によるBSSEの計算の流れはつかめたかと思います。実際に手を動かして計算してみることが上達の近道です。

1) まずは、水の2量体を分子モデリングソフトで作成します。ここでは、Avogadroを使って説明します。まずは水分子をある程度近づけて2分子描画させます。分子力場法にMMFF94(UFFよりも、極性の高い官能基を持つ分子で良い再現性を示します)を使うと酸素原子ともう1分子の水素原子が直線状に並んだ分子座標を作成できると思います。水2分子程度であれば、ある程度適当でも次のab initio計算で上手く最適化できるのでそれらしい構造をつくっておけばOKです。
hyd1
2) ここでは、Firefly(PC GAMESS)を使って構造最適化してみましょう。$STATPTグループのNSTEP(探索回数)は200程度まで増やし、HSSEND=.t.を追加して振動解析を行うオプションを追加しておきましょう。ここでは、理論/基底関数にHF/6-31G＊を指定します(もちろん、余裕がある方は基底関数の数を増やしてもかまいません)。
hyd2
3) あとは、Firefly(PC GAMESS)で計算して構造最適化されたことを確認してください(※振動解析の結果、虚数振動がないか確認してください)。計算結果をMaSKで表示させると分子間で水素結合が表示されていることがわかります。
hyd3
典型的な水の水素結合の距離は1.97A程度なので概ね近い値が得られています。分子間距離はツールバーの赤枠で囲ったボタンをクリックして原子を選択すれば、画面下に計測した距離が表示されます。
※結合距離が表示されない場合は、View >Status Barにチェックマークを入れてください

4) Outputファイルをテキストエディッタで開き”LOCATED”でキーワード検索してください。最適構造が以下の箇所に出力されるのでこれを抜き出します(赤枠で囲った部分)。
hyd4
5) Avogadroを開き、ツールバーのBuild >Cartesian Editorを開き、先程コピーしたCartesian座標をサブウインドウに貼り付けます。今度は、最適化した2量体のうち一分子ずつ(孤立系)のエネルギーを求める必要があります。まずは、青枠の①の分子座標だけ残し、②の部分は消した状態でApplyボタンを押しましょう。そうすると、2量体の片側の水分子のみ描画されるはずです。
hyd5
6) 次に、Extensions >GAMESS >Input Generatorで振動解析用のInputファイルを作成します。$CONTRL グループでRUNTYP=HESSIANを指定します。Avogadroで作成すると、自動で$FORCEグループが追加されますが、この行は消去してもかまいません。後は、もう片側の水分子(青枠の②の分子座標)についても同様にInputファイルを作成して計算してみましょう。
hyd6
一応、ここで新たに出てきたキーワードを説明しておきます。これらのキーワードは何も記述しなくても、デフォルト値が採用されます。オプションを変更したい場合のみ追記してください。

$FORCEグループはRUNTYP=HESSIAN指定時のオプションになります。
METHOD=ANALYTIC：通常の振動解析計算を実行します。SCFTYP=RHF, ROHFなどではデフォルトでこの解析計算が実行されます。別のオプションでMETHOD=SEMINUMにすると、解析的に計算された一次導関数の数値微分を行います。UHF、MCSCFを指定する場合これがデフォルトになります。
VIBANL=.TRUE.：振動解析を有効にするフラグです。RUNTYP=HESSIANの場合、デフォルトがTRUEです。
※) 他のオプションについてもっと詳しく知りたい場合は、GAMESS付属のマニュアルで$FORCEグループを参照してみてください。

7) これで2量体のエネルギー値及び、そこから抜き出した各分子単体でのエネルギー値の3つのファイルが揃いました。次は、片方の水分子をGhost軌道に置き換えて計算します。
hyd7
先程の2量体の構造最適化後のCartesian座標を準備します。Ghost軌道の指定は、原子番号に「-」(マイナス)を追加すればOKです(青枠で囲った部分)。1分子ずつなので2つのinputファイルを作成して、振動計算を実行します。

8) これで合計5つのエネルギー計算の結果が出揃いました。後は、結果を抽出してBSSEを算出するだけです。ここでは、結果の抽出にMaSKを使用して表計算ソフトでBSSEを計算してみましょう。
hyd9
OutputファイルをMaSKで開き、ツールバーのresults >Extend Summaryをクリックします。大きな赤枠で囲った部分が、熱力学的諸量になるのでこれをコピーして表計算ソフトにペーストします。

9) 表計算ソフトで作成した計算ファイルは、記事の最後でダウンロードできるようにしているので、是非活用してください。BSSEは、孤立系の分子のエネルギー値の合計から、Ghost軌道を使った疑似的な2量体のエネルギー値を引くことで算出されます。
tet01
つまり、水の2量化をHF/6-31G＊で計算した場合、CP法では0.93kcal/molの基底関数重なり誤差が生じるということが算出できました(※今回はエネルギー補正を行いません)。

記事の最後でダウンロードできる計算ファイルは、赤枠で囲った部分にMaSKでコピーした値をペーストすれば後は、自動で計算してくれます。また、MaSKの仕様上Ghost軌道を用いた計算結果の抽出は、各補正値が表示されないので手動で入力する必要があります。
tet02

BSSEを排除して相互作用エネルギーを求めよう！！

あとは、BSSEを補正値として加えるだけです。水の2量体のエネルギー値から孤立系の各分子のエネルギー値を引いたものが相互作用エネルギーに相当します。BSSEが0.93kcal/molだったので、さらにこの補正値を加えれば水の2量体の相互作用エネルギーが4.71kcal/molであることが無事算出できました。
tet03

おわりに

今回は、水の2量体における相互作用エネルギーの計算方法を紹介しました。CP法によるBSSEの計算方法は前回からの続きなので、わからなくなったら前回の記事を見返してみてください。理論/基底関数にHF/6-31G*を用いているので定量性はそこまで期待できません。時間がある方は、6-311++Ｇ**など基底関数の数を増やして計算してみてください。
一方、HF/3-21Gのようなレベルの低い基底関数で計算するとBSSEの値が大きくなり、B3LYP(DFT法)を用いた場合は、水素結合のエネルギーを過大評価してしまうなど、これまで述べてきたことが実際に『体験』を通して実感できると思います。

CP法の他にGAMESSにはBSSEを求める方法として「KITAURA-MOROKUMA ANALYSIS」というオプションがあります。これは、分子間相互作用をエネルギー分割法によりBSSEを算出する方法です。CP法に比べ、Ghost軌道や孤立系の計算が必要なく簡単なのですが、色々と制約があります。
次回はKITAURA-MOROKUMA ANALYSISの実行方法について紹介します。 

この記事で使用した計算ファイルを以下に置いておきます。
http://pc-chem-basics.blog.jp/BSSE%20Calc..xlsx
http://pc-chem-basics.blog.jp/BSSE.zip

【関連記事】
BSSE(基底関数重なり誤差)について知ろう！！
熱力学的諸量をみてみよう  
KITAURA-MOROKUMA法でBSSEを計算してみよう