GAMESSでかんたん、遷移状態探索。ただしMacに限る。

GAMESSで化学反応の遷移状態の探索を行う方法について。シクロペンタジエンとシクロペンタジエンからノルボルネンが生成する反応を例に


↑まずはモデリング。どんなモデリングソフトを使用してもいいのですが、ここでは"Avogadro"を使用。Avogadroを起動します。


↑ノルボルネンを適当に描き、ジオメトリ最適化（Avogadroの力場計算機能を使用）を行います。



↑ノルボルネンの、シクロペンタジエン部分とエチレン部分の距離を長くします（とりあえず最適化構造の距離の2倍）。ふたつ結合があるので、ひとつを切断し("Draw Setting"の"Adjust Hydrogens"のチェックを外すことをお忘れなく)


↑View → Properties → Bond Properties… から


↑シクロペンタジエン部分とエチレン部分をつなぐ結合の長さを示すところに…


↑数値(1.5 × 2 =3.0)をを入れると…


↑ノルボルネンはこうなります。


↑名前を付けて保存します。


↑拡張子が".gamin"なので、これを".inp"に変更します。


↑Finderがこんなことを聞いてくるので、「".inp"を使用」を選択。


↑こうなります。


↑"テキストエディット"に入力ファイルのアイコンをドラッグ＆ドロップします。


↑こんなふうに表示されます。1行目は消去しましょう。


↑GAMESSの構造最適化用のパラメータを入力します。予め作っておいたものをコピー＆ペーストすると楽です。ここではDFT, B3LYP/3-21G(d)レベルで計算。

"RUNTYPE=OPTIMIZE"に注目。いきなり遷移状態の探索はしません。とりあえず、Avogadroで作成したモデルを構造最適化し、シクロペンタジエン＋エチレンに最適化されるか、ノルボルネンに最適化されるかを調べます。
完全に最適化する必要はありませんので、探索の回数(NSTEP)は30程度でいいと思います。


↑ファイルを保存し、入力ファイルのアイコンを"GamessQ"にドラッグ＆ドロップします。

"Job Options"ウィンドウが出ますので、適当な数値を入力し、"OK"をクリック。構造最適化が始まります。


↑計算が終了(途中でも可)したら、ログファイルをAvogadro（MacMolPltでもいいです）で開きます。シクロペンタジエンとエチレンになっています。使用したモデルを初期構造とした場合、遷移状態を探索できず、シクロペンタジエンとエチレンの構造になってしまうでしょう（運良く遷移状態が見つかることもありますが…）。シクロペンタジエンとエチレンの距離が長過ぎるということです。

そこで、Avogadro（くどいようですがモデリングソフトならなんでもいいです）でGAMESS入力ファイル(.inp)を開き、シクロペンタジエン部分とエチレン部分をつなぐ結合を（消えてしまっている場合は描いて下さい）、1.5Aと3.0Aの間の2.25A（大体でいいです）にします。で、同じように入力ファイルを作成、構造最適化を行います。


↑得られた構造はノルボルネン。シクロペンタジエンとエチレンの距離が短過ぎるということです。

今度は2.25Aと3.0Aの間、その次は…と、構造最適化をしてシクロペンタジエン＋エチレンに最適化されるか、ノルボルネンに最適化されるかを調べます。シクロペンタジエン＋エチレンに最適化される初期構造、ノルボルネンに最適化される初期構造での、シクロペンタジエン部位とエチレン部位の距離の差が0.05Aくらいになるまでこの作業を続けます。

この作業が終わったら、シクロペンタジエン＋エチレンに最適化される初期構造の中で、シクロペンタジエン部位とエチレン部位の距離が最も短いモデルのGAMESS用入力ファイルをテキストエディットで開き、遷移状態探索のためのパラメータを入力します。


↑注目する箇所のひとつは"RUNTYPE=SADPOINT"。鞍点探索を行います。
もうひとつは”HESS=CALC"。振動解析を行ってから（必須です）鞍点探索を行います。予め振動解析を行っているのであれば"HESS=READ"となります（この場合は$HESS〜$ENDのデータが必要）。
NSTEPには大きな数値(1200)を入れています。

ファイルを保存し、入力ファイルのアイコンを"GamessQ"にドラッグ＆ドロップします。構造最適化の時とは異なり、"HESS=CALC"の場合、まず振動解析を行うので、鞍点の探索がスタートするのに時間がかかります。"HESS=READ"の場合は鞍点の探索がすぐにスタートします。

注意してほしいのは、GAMESSでは鞍点の探索がなかなか終わらないケースがしばしば見られることです。計算の進行具合を、ログファイルを開いて時々チェックする必要があります。計算が終了していなくても、分子全体がスライドするような最適化が続いているようであれば、鞍点探索ができている可能性大、です。
その場合、「スライドするような最適化」されている中で最もエネルギーの低い構造をファイルに書き出します。





↑"MacMolplt"だとFile → Export… から"GAMESS $DATA group (*.inp)"でエクスポートします。

遷移状態を思われる構造がピックアップできました。


↑次に、この構造が遷移状態であるかどうか、振動解析でチェックします。エクスポートしたファイルをテキストエディットで開き、振動解析用のパラメータを入力します。
"RUNTYPE=HESSIAN"に注目。

ファイルを保存し、入力ファイルのアイコンを"GamessQ"にドラッグ＆ドロップします。振動解析が始まります。


↑振動解析が終わり、吐き出されたログファイル(.log)を開き、振動数を表示します(↑MacMolPltの場合)。拡張子が".dat"のファイルを消さずに残しておくことをお忘れなく（私はよくやらかします）。


↑マイナスの数値（虚振動）がひとつだけあるのに注目。この構造が遷移状態である可能性があります。ない場合は遷移状態ではありません。MacMolPltの場合、振動数をクリックすると、モデルの振動の方向が表示されます。マイナスの数値をクリックすると、「それらしい」振動方向の矢印が表示されるはずです。

この構造が遷移状態であることを確認するために、もうひと作業必要です。
振動解析で使用したGAMESS入力ファイルを使ってIRC解析を行います。


↑テキストエディットで開き、パラメータを入力します。"RUNTYP=IRC"に注目。
"FORWRD=.TRUE."にも注目しておいて下さい。


↑次に、拡張子が".dat"のファイルをテキストエディットで開き、"$HESS"の行から"$END"の行までをコピーします（$HESSをキーワードに検索すると楽です）。$HESS行頭のスペースも含めてコピーしましょう。

コピーしたものを、GAMESS入力ファイルの最後のところにペーストし、ファイルを保存します。

このファイルを使って、もうひとつGAMESS入力ファイルを作成します。
入力パラメータ中の"FORWRD=.TRUE."を"FORWRD=.FALSE."に変更し、別名で保存します。

そして、作成した2つのファイルのアイコンを"GamessQ"にドラッグ＆ドロップします。計算開始。
".log"ファイルか、".trj"ファイル（こちらがおすすめ）を開き、構造がどう変化するか確認します。初期構造が遷移状態であれば、「シクロペンタジエンとエチレン」か「ノルボルネン」へと変化しているはずです。