計算化学 ソフトウェア

反応の進んでいく過程をシミュレーションしたい人にとって、反応が進む瞬間の構造(遷移構造)を求めるのは一つの重要なテクニックです。 そのやり方の一つである、QST法による遷移構造探索をGaussianとGaussViewで行う方法を、連載で紹介していきます。

構造最適化計算を、分子のポテンシャルエネルギーを極小化する構造を探索する計算ととらえると、遷移構造を構造最適化計算で求めるのが 困難であるのがわかります。反応前の構造あるいは反応後の構造に少しでも近いと、そちらの方がエネルギーが小さいので、それらの構造に落ち着いてしまうためです。
これを、ポテンシャルエネルギー曲面において考えてみると、遷移構造は鞍点(あんてん)という概念に対応します。 鞍点の図がWikipediaのページ右にありますので、 初耳の人は見てみてください。鞍点とは、文字通り鞍のようになっている曲面において、ある一方向についてのみ極大で、それ以外の方向については 全て極小になっている箇所と言えます(この場合、構造最適化は全方向への極小化計算となりますね)。 そこで通常の構造最適化計算とは違った構造探索計算が必要になるわけですが、 今回紹介するQST法では、反応の前後の安定構造(計2つ)をもとに、 その間の遷移構造を探索します。簡単な題材として、アセトアルデヒドとビニルアルコールの間の遷移構造を求めていきます。

QST法は2つの平衡体(reactantとproductと呼びます)の間の遷移構造を求めます。これには次のステップを踏みます。

1. reactantの安定構造を求める
2. productの安定構造を求める
3. reactantとproductの両方で各原子のTagを揃える
4. Transition State Optimization計算する
5. 遷移構造が出力されているかチェックする

まずは、reactant(今回はアセトアルデヒドとします)の安定構造を求めましょう。 GaussViewでアセトアルデヒドをモデリングします(詳細は割愛します。GaussViewのHelpをご参照ください)。

• 図1．GaussViewでアセトアルデヒドのモデリング結果

ここでステップ3.に関係してくることですが、気を付けることがあります。それは原子のTag(番号)です。周期表の原子番号のことではなく、 分子内で便宜上付けられた番号で、Gaussianが内部的に気にする情報です。まずは見てみましょう。[Edit]メニューから[Atom List…]を選びます。 すると、Atom List Editorダイアログが現れます。

• 図2．Atom List Editorダイアログ

このダイアログでは各原子のシンボルや座標、角度等を編集できます。ダイアログ上部にRow,Highlight,Display,Tag,Symbol,NA,…などと 書いてありますが、HighlightとTagの列が今回気を付ける部分です。Highlight列の任意の●項目をクリックすると、対応する原始が明るい水色で ハイライトされ、その原始の中に番号が表示されます。この番号がTag列に記載されている数値です。

• 図3．1原子だけハイライトされたアセトアルデヒド

つまり、分子内の各原子に対応するTagがあるのです。またこのTagの値はAtom List Editor内で変更(他の原子と入れ替え)できます。ここでは、全原子にどのTagが振られているのかを(PrintScreenキーで スクリーンショットを撮ってペイントで名前を付けて保存する等して)保存しておきましょう。その後、Atom List Editorを閉じて、アセトアルデヒドの構造最適化計算を 実行します。これには[Calculate]メニューから[Gaussian Calculation Setup…]を選んで、Job Typeコンボボックスで[Optimization]を指定したら、 [Submit…]ボタンをクリックして計算実行です。

続いて、product(ビニルアルコール)の構造を作りましょう。product分子を作る際、前回述べた、原子のTagをreactantと対応するように作っていきましょう。

• 図1. GaussViewでビニルアルコールのモデリング結果

上図のように作っていたら、前回撮ったスクリーンショットと比較して、各原子とTagが、前回のreactantのものとズレていることが Atom List Editorを見てわかりました。そこでAtom List Editorで各原子のTag値を変更して(入れ替えて)Enterキーで適用します。

• 図2. Atom List EditorダイアログでTagを変更

productとreactantにおいて対応する原子のTagが全てそろったら、構造最適化です。これには[Calculate]メニューから [Gaussian Calculation Setup…]を選んで、Job Typeコンボボックスで[Optimization]を指定したら、 [Submit…]ボタンをクリックして計算実行です。

productとreactantの安定構造が求まったでしょうか？それぞれの構造最適化結果(.logファイル)をGaussViewで開いて表示しておきましょう。 これらを基に今回遷移状態計算に用いるQST2法のインプットを作っていきます。

まず、[File]メニューから[New]⇒[Create Molecule Group]と選択し、新規インプットを作成するウィンドウを表示します。 つづいて、product(ビニルアルコール)の構造最適化結果ウィンドウを最前面にし、[Edit]メニューから[Copy]を選びます。 その後、[Edit]メニューから[Paste]⇒[Add to Molecule Group]を選択します。これで新規インプットにproductの構造が貼り付けられました。

同様に、reactant(アセトアルデヒド)の構造最適化結果をコピーし、新規インプットに貼り付けます。この状態で新規インプットのウィンドウ上部にある 表示モード切替ボタンをクリックし、multiview modeに切り替えます。

• 図1. multiview modeへ切り替え

multiview modeでは、編集中のMolecule Groupがタイル状に並べて表示されます。今回の例ではproductとreactantが左右に表示されます。 ここで、[Edit]メニューから[Connection…]を選択します。すると、Connection Editorが表示されます。

• 図2. Connection Editorとmultiview mode

Connection Editorが表示されている間、multiview modeのウィンドウでは、Molecule Groupごとに、各原子に割り振られたTagが表示されます。 ここで、双方の構造で各原子の対応が合うよう、Tagを調整します。そのためには、productの原子を左クリックで選択して(ハイライトされます)、 それに対応するreactantの原子を右クリックで選択します(ハイライトされます)。なお、productとreactantの作業エリアのフォーカスを切り替える際に 切り替え先のウィンドウを１回クリックする必要がありますので注意してください。

• 図3. 調整後のTag

各原子のTag対応が合ったら、[OK]ボタンをクリックしてConnection Editorを閉じ、QST2法で構造最適化計算を行います。 それには、[Calculate]メニューから[Gaussian Calculation Setup]を選択します。[Job Type]にOptimizationを選択し、 [Optimize to a]欄で TS (QST2) を選択します。また、[General]タブにて[Write Connectivity]チェックボックスをオフにします。 そして[Submit…]ボタンをクリックして計算実行です。

前回、QST2法で構造最適化計算を実行しました。今回は、その計算結果で遷移構造が出力されているかどうかをチェックします。

チェックすべきポイントは次の2つです。

(1)収束判定「Converged?」が全部YESになっていること。
(2)虚の固有値が唯一つ存在していること。

それぞれのチェック方法を以下に説明します。

QST2計算の過程を直接確認する必要があるので、QST2法の構造最適化計算の出力ファイル(.logもしくは.out)をテキストエディタ(Windowsのメモ帳等)で開きます。

(1)収束判定「Converged?」が全部YESになっているか？

出力ファイル中を「Converged?」で検索します。たいていの場合、ファイル中で複数見つかります。そのうちで、ファイル中で一番後ろに出現する Converged? を探して、 その下4行を確認します。次のような行が見つかります。

 Item                       Value     Threshold  Converged?
 Maximum Force            0.000114     0.000450     YES
 RMS     Force            0.000044     0.000300     YES
 Maximum Displacement     0.000668     0.001800     YES
 RMS     Displacement     0.000341     0.001200     YES

これらの4行全ての Converged? の項目が YES になっていることを確認します。 もし一つでも YES でない行があれば、それは遷移構造ではありませんので、 入力ファイルをこれまでの手順に沿って正しく作成したかどうかを見直してください。

(2)虚の固有値が唯一つ存在しているか？

(1)で見つけた行の70行程前にある、「 Eigenvalues —」で始まる複数の行を探します。次のような行が見つかります。

Eigenvalues ---   -0.21139   0.01293   0.03911   0.05263   0.08490
Eigenvalues ---    0.11446   0.16127   0.16377   0.17601   0.24041
Eigenvalues ---    0.36010   0.36275   0.36462   0.52903   0.60070

これらの値が、一つだけ負の値をとっていて、かつ負の値が複数あらわれていない ことを確認します。負の値があらわれなかったり、あるいは負の値が複数あらわれた 場合は、遷移構造ではありません。

上記(1)(2)が全て満足できている場合に限り、最後のチェックに進みます。

(1)収束判定「Converged?」が全部YESになっているか？
(2)虚の固有値が唯一つ存在しているか？

前回のように上記(1)(2)が全て満足した結果となったでしょうか？満足できた場合は、最後に求まった構造で振動数計算を実行します。

[Calculate]メニューから[Gaussian Calculation Setup]を選択します。

[Job Type]にFrequencyを選択します。計算理論や基底関数系をそろえた後、[Submit…]ボタンをクリックして計算実行です。

振動数計算が正常に完了したら、振動数計算にて虚の固有値が唯一つとなったなら遷移状態と認めます。

具体的には、出力ファイル(.log)をテキストエディタ(Windowsのメモ帳など)で開きます。 ファイルの末尾から、「 Eigenvalues —」で始まる複数の行を探します。次のような行が見つかります。

Eigenvalues ---   -0.39161   0.02340   0.04121   0.05702   0.07807
Eigenvalues ---    0.09233   0.15318   0.19388   0.24331   0.28561
Eigenvalues ---    0.32554   0.69021   0.89937   1.05213   1.31076

これらの値が、一つだけ負の値をとっていて、かつ負の値が複数あらわれていないことを確認します。 負の値があらわれなかったり、あるいは負の値が複数あらわれた場合は、遷移構造ではありません。

虚の固有値が唯一つとなっていれば、その構造が遷移状態です。(次回配信は2012年7月17日予定)

Gaussian入門者の皆様に、Gaussianを使いこなして高度な計算化学者に飛躍していただけるよう、Gaussianに関する様々な情報を発信しております。メールニュースの配信をご希望の方は、以下のフォームよりお願いいたします（お問い合わせ内容の項に「Gaussian入門メールニュース配信希望」とご記入ください）。

メールにてお申し込みされる方は、以下の必要事項をご記入の上、送信してください。

◆電子メール件名 : Gaussian入門メールニュース配信希望
◆記入項目：
1. 氏名 :
2. E-mail：
3. 法人名/部署名：
4. 住所：
5. 電話番号：