この項ではDiels-Alder反応を題材にしてIRC解析の方法を述べる。なぜDiels-Alder反応を選んだかというと、いろいろとIRC計算を試した中で最もやりやすかったからである。
正直IRC計算の敷居は中々高く、大概回りくどい解説が多いので1回生向きではない(これは仕方がない。というのもこのホームページに書いている目的が趣味でやる計算なのに対し、本来の目的は研究するためだからである)。それゆえ気楽に読めるものを目指して記すつもりである。
Diels-Alder反応
さて、IRC計算の前にDiels-Alder反応とは何かについて軽く述べる。この反応にはフロンティア軌道が深く関わる。下の図を見てほしい。
上側がcis-ブタジエンとそのHOMO、下側がエチレンとそのLUMOである。
HOMOとLUMOがバチッと合致する(?)ので反応が進行する。なおこの反応はcis-ブタジエンのLUMOとエチレンのHOMOを考えても同じ結果となる。生成物はシクロヘキセンとなる。
IRC計算の手順
さて本題のIRC計算である。手順は以下の通り。
(0)Hesse行列の計算 (RUNTYP=HESSIAN)
(1)遷移状態の計算 (RUNTYP=SADPOINT)
(2)正しい遷移状態を計算したかどうかの確認 (RUNTYP=HESSIAN)
(3)遷移状態から反応系または生成系に向かってIRC計算(計2回) (RUNTYP=IRC)
(1)遷移状態の作成
まず初めに遷移状態の構造を作る。ここが最大の難点である。遷移状態はWinmostarで作ることをお勧めする。「作る」というのもこの反応は簡単なものなので下の図のように適当に(ある程度の適切さで良いという意味)作って構わない。とにかく初めのうちは原子移動などを用いて作ってみる。
ここまで出来たら取り合えず構造をpdbファイルで保存し、Facioで開く。開いたらRUNTYP=SADPOINT、HESS=CALC、STO-3Gで計算する(今回ここでの計算にあまり精度は必要ない)。すると以下のような構造になってくれる(はず)。
成功すると横のパネルに
*** Saddle Point Located ***
と表示される(エネルギーの結果の上にあるはず)。この後IRC計算したい計算精度でもう一度同じことをやる。
なお失敗すると、
*** Failure to Locate Stationary Point ***
と表示される。
ここで注意したい点が何点かある。似たような構造になった人は下の注意点を無視して次に進むとよい。失敗した人は読むこと。
(※1)
*** Failure to Locate Stationary Point ***
と表示された場合、遷移状態が見つかりませんでした、ということである。この場合遷移状態を作り直さねばならない。この作業がIRC計算における最大の山場である。Diels-Alder反応は研究されまくって遷移状態がどのようなものかほぼ分かっているのであまり苦労はしないが、未知の反応解析となると胃が痛くなるに違いない。
ところで上の文章の意味は
停留点(stationary point)の探索に失敗しました
という意味で、停留点とは極大値、極小値、鞍点などのことである。鞍点って何?という人は検索するか、微積分学や解析学の教科書を開いてみること。鞍点は英語でsaddle pointと言う(自転車のサドルみたいな曲面の形)が、これがRUNTYP=SADPOINTとした理由である。とにかく遷移状態の探索に失敗した場合は微調整して計算を繰り返すことである。
(※2)
HESS=CALCとしたことについてである。これは「Hesse行列を読み込まず、遷移状態計算をする前にHesse行列の計算をしますよ」という意味である。読み込むときはHESS=READとする。簡単な計算ならばこの設定で良いが、複雑な反応系であればHesse行列は先に計算したほうが良い。この場合RUNTYP=HESSIANとして先に計算しておく。計算した後、FileよりHesse行列を読み込む。読み込むと、
*** $HESS group was successfully retrieved. ***
とでる。
(2)遷移状態の確認
さて遷移状態の計算が出来たらその構造が本当に遷移状態か確認するため、RUNTYP=HESSIANとして計算する。ここで
*** Firefly Output ***
Number of Imaginary Frequency : 1
**********************************************
There is Only ONE Imaginary Frequency at Normal Mode #1
**********************************************
と表示されれば先ほどの構造は遷移状態でした、ということになる。これで終わればめでたしめでたしだが、普通ここまですんなり行くことは少ない。もし
*** Firefly Output ***
**********************************************
This is NOT a stationary point on the molecular PES
The Vibrational Analysis is NOT Valid !!
**********************************************
Number of Imaginary Frequency : 2
などと表示されると、それは遷移状態ではないのでもう一度遷移状態を調整して計算しなおす必要がある。
※Number of Imaginary Frequencyとは虚振動数解の数のことである。固有振動数が虚数ということだが、意味が分からなければ解析学の教科書と量子化学の教科書のIRC計算の部分を読むこと。虚振動数が1つというのは遷移状態が一次の鞍点であるためである。
とにかくひたすら調整して、虚振動が1つになるまで努力すること。
(3) IRC計算
めでたく遷移状態が求まれば、あとは遷移状態からポテンシャルエネルギー面の山を下るだけである。下る方向は反応系と生成系の2つがあるので2回計算を行う。まず以下の赤線で囲った部分に注意する。
まずRUNTYP=IRC、HESS=READとする。この前にFileよりHesse行列を読み込んでおくこと。次に右下にあるMore Optionsを選択し、IRCの設定ウィンドウを開く。遷移状態から計算する向きをFORWRD=.TRUE.と.FALSE.の2つ選べるので、それぞれ計算する(計2回計算)。NPOINTは探索する点の数で多いほど計算量が増える。STRIDEは計算する点の間隔で、小さいほど計算量が増える。ここはパソコンのスペックと計算精度によって微調整する。
計算が終わると
*** IRC Calculation Performed ***
と表示されるはずである。もし、
*************************************
ADDRESS 0x00734840 HAS INITIATED PROGRAM ABORT BECAUSE OF FATAL ERROR(S)
*************************************
というエラーが出たら、何かを間違えているはずである。大概は、計算精度を途中で変えてしまっている(SADPOINT,HESS,IRCの計算はすべて同じ計算精度でなくてはならない)、電荷の間違い(-1を+1にしているなど)、スピン多重度の間違いなどである。
IRC計算の結果表示と分析
結果表示は
File > Load New GAMESS > Firefly IRC(*.irc) and WinGamess IRC(*.trj)
から*.ircのファイルを選択する。そんなファイルが見当たらない人は、Fireflyのフォルダの中にIRCDATAというファイルがあるはずなので、それを選択する。
ファイルを読み込んだら、
Tools > Viewers >IRC Data Viewer
を選択する。すると以下のような画面が出る。
Show Animationを選択するとパラパラ漫画のように反応の様子を見ることが出来る。
Show IRC Energy Plotを選択すると以下のような画面がでる。
RMS Gradientはエネルギー勾配の平均二乗偏差を意味し、RMSは二乗平均平方根(root mean square)を意味する。ここで下のSave As Textを選択すると、Fireflyのフォルダに
IRC_Energy.txt
IRC_Grad.txt
という2つのファイルが出るはずである。このテキストフォルダ内の数値をExcelやgnuplotなどの数値処理ソフトに入れて、エネルギー断面図を描くことができる。
こうして今回の結果を纏めたのが以下の図である。計算精度はHF/3-21G(d)で行った。
コメントをお書きください