IRC 計算がうまくいかない時の対処方法

反応経路探索において、もっともフラストレーションがたまるのが IRC 計算がうまくいかない時です。

opt=ts がうまくいき、虚振動も一つ!もう少しで反応エネルギーなどの情報が手に入る!

、、、といったところで、いつも足を引っ張ってくるのが IRC 計算です!

IRC 計算がうまくいかない時には試行錯誤が必要ですが、それらの努力が必ずしも報われるとは限りません、、、。でも、最低限試した方が良いと思うことをまとめてみました。

この記事に書いてある方法以外にも良い解決方法をご存知の方がいらっしゃいましたら、ぜひコメント欄に書いてください。そうすればみんなで情報が共有できると思います。

このようなエラー対処方法を、研究者 or 研究室特有の技術として秘密にしている方々もたくさんいらっしゃると思いますが、そういうのは科学の進歩にとって害にしかならないと思います。どんどん公開しましょう。

英語版When the IRC calculation doesn’t work…

stepsize を変えてみよう!

中間体の構造最適化では maxstep=30、IRC 計算では maxstep=20、遷移状態構造最適化では maxstep=10 が良いとされています。
IRC 計算がうまくいかない時は、step size を少し小さくしてみるとうまく行くという場合があります。

L115 を使う

うまくいかないときは、古いバージョンの gaussian を使ってみることも一つの手です。IRC は default では L123 ですが、use=L115 というオプションを用いることにより gaussian03 のプロトコルをしようすることが可能です。ちなみに、maxstep の上限は gaussian03 で 400、gaussian09 で 1,000 となっています。

recorrect オプションを付ける

IRC が 2 点だけ進んで終了してしまうエラーのときに、irc=(recorrect=never) というキーワードを使うとうまく行く場合があります。IRC の各ステップは、閾値と比較してクリアしていれば次の correction に行くという流れになっています。recorrect = never は閾値テストを行わないオプションです (supress thethreshold test)。ちなみに、irc=tightverytight を付けると閾値の値を下げられます。default は recorrect=YES です。

rcfc を使わない

遷移状態構造最適化で使用した chk ファイルの情報を読み込む rcfc を使うと IRC 計算がうまくいかない時があります。
このような場合は、calcfc を使いましょう!

ちなみに、calcfc で うまくいかない時に calcall を使うという良いという都市伝説がありますが、管理人はうまくいったことがありません。また calcall は非常に時間がかかります。

calcall を試してみろと言ってくる人に限って calcall の使用経験が無かったりします。。。

追記)本記事の公開後、calcall でうまくいったことあるよ!という意見が寄せられました。うまく行くこともあるようですが、非常に時間がかかります。

gaussian の違うバージョンを使う

Gaussian09 の IRC 計算には一部バグがあることが知られています。
Gaussian03 を使ってみましょう!

また、Gaussian 09 の場合でも g09 B.01 や g09 D.01 など revision ごとに差があります。違う revision を使ってみましょう。

Gaussian16 で IRC が改良されることを期待していましたが、どうやら大きな変化は無いようです。

GRRM を使う

この際はっきり言いましょう、Gaussian の IRC はクソです。
Gaussian でいろいろと試行錯誤するよりも GRRM を使ったほうが早いです。GRRM の IRC でうまくいかなかったことはほとんどありません。

他の IRC のアルゴリズムを使う

Gaussian には様々なIRC のアルゴリズムが搭載されています。
うまくいかない場合はアルゴリズムを変えてみるのも一つの方法です。

参考: Gaussian 社の IRC ページ

g16 では、EulerPC がデフォルトとして設定されていますが、LQADVV なども選択可能です。
ちなみに、GRRM では PageMcIver の方法で IRC をトレースしています (参考文献 1,2)。

IOp(2/12=3) の指定も役に立つかもしれません。

どうしてもうまくいかない時

上記の方法を試してみてうまくいかない時は、妥協しましょう。
QM 以外の方法で反応経路を探索してみるのも一つの方法かもしれません。MD とか NEB とかでも良いんじゃないでしょうか?

また、微妙に構造の違う遷移状態構造が真の遷移状態かもしれません。TS の探索からやり直してみるというのも、一考の価値ありです。

次のような言い訳を使っている論文もあります (参考文献 3)。

In nearly all of these cases, attempted IRC calculations only ran a few steps before stopping at spurious minima, behavior which we often experience when working with flat regions of potential energy surfaces.

“Pentalenene formation mechanisms redux”
Michael W. Lodewyk, Dan Willenbring and Dean J. Tantillo Org. Biomol. Chem. 2014, 12, 887–894.

それでも諦められない人へ

網羅的計算により PES を作るのも、人を納得させる一つの方法だと思います。
虚振動のベクトル方向へ原子を動かすスクリプトを書いて opt & freq すれば、それっぽい PES が作れると思います。

管理人のデフォルト設定

ちなみに、管理人は以下の keyword 指定をいつも使っています。

#p B3LYP/6-31+G** int=(grid=ultrafine) scf=(maxcyc=300,direct,tight) irc=(lqa,forward,rcfc,noeig,maxcyc=30,maxpoints=30,stepsize=15) geom=check guess=read IOp(2/15=3) IOp(2/12=3)

ポイントは、lqa を使用しているところと、IOp(2/15=3) IOp(2/12=3) を指定しているところです。

もっとこうした方が良いよ!という意見がありましたら、是非コメント欄にお願いします。

参考文献

  1. “On evaluating the reaction path Hamiltonian” M. Page and J. W. McIver Jr., J. Chem. Phys. 1988, 88, 922-35. DOI: 10.1063/1.454172
  2. “Following steepest descent reaction paths – the use of higher energy derivatives with ab initio electronic-structure methods” M. Page, C. Doubleday Jr., and J. W. McIver Jr. J. Chem. Phys. 1990, 93, 5634-42. DOI: 10.1063/1.459634
  3. “Pentalenene formation mechanisms redux”
    Michael W. Lodewyk, Dan Willenbring and Dean J. Tantillo Org. Biomol. Chem. 2014, 12, 887–894.

関連する記事

是非フォローしてください

最新の情報をお伝えします

2 comments

  1. IRCは反応座標を刻むために各ステップごとに何らかの最適化が必要になります。(この方法だと美しいPESが得られる。)単に遷移状態の確認をするだけなら IOP(1/8=xx)を小さめにしてから遷移状態を虚の振動モードの方向に歪ませて OPTすることで両端につながるのを確認するだけでいいんじゃないかなと最近は思ってます。

コメントを残す(投稿者名のみ必須)