今回は、遷移状態の候補構造を探索するMinimum Energy Path計算Fireflyで実行してみましょう。反応座標の計算は、反応部位の結合距離や角度などを一方向に変化させながら、繰返し計算をおこない遷移状態を探索する方法です。
mepg00

はじめに知っておきたいこと

基本的には以前紹介したMOPACを使った方法をマスターしておけば簡単だと思います。分子座標は、MOPACの内部座標形式(ZMTMPC)を使うと便利です。
MoCalc2012を使用できる環境であれば、入力ファイルの作成も比較的簡単です。ただし、分子座標の結合関係をきちんと定義しておかないと、反応座標の計算を行っても異なる目的物に収束したり、発散したりするので注意が必要です。

入力ファイルの作成方法は?

臭化メチルとアンモニアの2分子求核置換反応を例に計算してみましょう。理論/基底関数はHF/3-21Gを使用します。入力ファイルの内容は以下の通りです。赤枠で囲った部分が新しいキーワードです。
megg01
$CONTRLグループ
RUNTYP=RSURFACE: 計算方法をRelaxed Surface Scanに指定します。各点で構造最適化を行いながらポテンシャルエネルギー面(PES)の探索を実行します。変化させられるのは距離のみ、角度は利用できません。

*) GAMESS(US)の場合、RUNTYP=SURFACEを使用します。しかし、FireflyのRSURFACEとは異なるので注意が必要です。GAMESS(US)のSURFACEは、構造最適化は行われない仕様になっています。例えば、原子間の距離を変更させて相互作用エネルギーの変化を確認する場合など、各構造で最適化が必要でない場合に適しています。

$SURFグループ
VECT1(n)=1: カッコ内にIFREEZで指定した値を指定します。
ORIG1=n: 座標の初期値を指定します。ゼロの場合は$DATAで指定した距離に対応します。
DISP1=n: 座標の変化量(ステップサイズ)を指定します。
NDISP1=n: 何回計算するか(座標に必要なステップ数)を指定します。

$STAPTグループ
IFREEZ(1)=n: 固定する構造定数を指定します。詳しい使い方は、こちらの記事を参照してください。
METHOD=GDIIS: GDIISアルゴリズムを使用するフラグです。GDIIS法は、構造最適化の収束を加速および安定化させるために用います。この手法は、以前の最適化サイクルで辿った座標の線形結合で新たな座標を計算することで、探索表現を豊かにし計算を加速する意図があります。比較的大きい系や、条件が厳しい場合、平坦なポテンシャルエネルギー面を持つ分子などを計算する際に推奨されます。

つまり、ここでは反応座標を10番目の窒素原子で指定し(-1)、2.501551Åから-0.05Åづつ一番目の炭素原子との距離を変化させながら最適化計算を20ステップ(計-1.00Å)行うことを指定したことになります。

後は、これまで紹介してきたキーワードなので過去の記事を参照してください。MAXIT(SCF計算の繰り返し回数)はデフォルトの30回ではこころもとないので少し大きくしておきましょう。

MoCalc2012を使用する場合、詳しい使い方は「部分構造最適化をしてみよう」の記事を参考にしてみてください。
megg05
まず、初期構造を読み込んだらZ-Matrixエディッタ-で編集します。窒素原子を反応座標に選び、最適化フラグを「S」(Scan=反応座標(-1))にします。ScanのStartとStopで座標の初期値と終点を、#Stepsでステップ数を指定します。
megg06
メイン画面でJob TypeをReaction Coordinate/Grid Calculationに設定すればOKです。後は、計算を実行するだけです。反応座標値毎に繰返し計算を行うので相応のCPU時間を要するので気長に待ちましょう。

計算結果から遷移状態(候補)を利用するには?

計算で得られたエネルギープロファイルの極大点が、遷移状態近傍の構造であると推測されるので、構造を抜き出し遷移状態の初期構造として利用できます。

MoCalc2012の場合は、計算終了後にCoordinateボタンをクリックしアニメーションをプレビューすることができます。後は、エネルギープロファイルの極大点の構造をXYZ形式などで保存すればOKです。
megg04
MEPGG00
MaSKの場合は、ツールバーのEボタンをクリックすると表示される、エネルギーダイアグラムの極大点の棒グラフをクリックします。エネルギーの極大点の構造が表示されたら、Fire >Export ToからXYZファイルに保存して、遷移状態計算用の入力ファイルに利用します。
megg03
MacユーザーはMoCalc2012が使用できないのでMaSKやWebMOを利用してみてください。ちなみに、候補構造をHF/3-21Gで遷移状態計算した結果、上手く遷移状態が見つかりました。計算はFirefly(PC GAMESS)で行いました。
megg07

おわりに

反応座標の手法を使った遷移状態探索は、一回の計算で上手く求まれば良いのですが大半はトライアンドエラーの繰返しです。特別な理由がない限り、GAMESSよりもMOPACで遷移状態の候補化合物を探索する方が入力ファイルの作成の簡便さや速度の面からおすすめです。もちろん、探索の段階から半経験的手法ではなくab initioやDFTの手法を用いたい場合は仕方ありません。
初めのうちは、分子座標に関する反応座標の結合関係が上手く設定できずに、エラーや求めたい構造に収束しない場合が多々あります。その場合、MEP計算が成功したファイルを再利用し原子のみを置き換えて計算してみると良いでしょう。

この記事で使用した計算ファイルを以下に置いておきます。
http://pc-chem-basics.blog.jp/NH3+CH3Br(MEP).inp
http://pc-chem-basics.blog.jp/NH3+CH3Br(MEP).out
http://pc-chem-basics.blog.jp/NH3+CH3Br(TS).out

※) 当ブログの連載が書籍になりました!!興味がある方は、お手に取っていただければ幸いです!!
(Amazon購入ページ)(書籍の紹介記事
ICam2
【関連記事】
遷移状態を「MEP法」で探索しよう
部分構造最適化をしてみよう