遷移状態の構造最適化が止まったら[GaussianによるTS計算にまつわるTips]

遷移状態計算って、やってみると意外となんとかできるもんですが、それでも基底状態と比較して、やや面倒なことも多いと思います。

私はラジカル系の計算をすることが多いのですが、ラジカル的な解を求めようとしているのに、勝手にイオン的な（閉殻的な）遷移状態がでてきてしまったり、難しいことが多いです。

なかなか、ノウハウを共有している資料も多くないので、解決に向けたノウハウを少しずつ貯めて行きたいと思います。（Gaussian09を用いたDFT法です。）
（もしかしたら、こちらも役立つかも？）

[閉殻分子同士のラジカル性の反応を見たいのに、勝手に求核攻撃のような解に収束してしまう]
一重項分子の結合のホモリシスを伴いながら進行する系で、どう考えても遷移状態ではラジカル性がないとおかしいのに、求まったTSの軌道を見ると完全に閉殻、つまりイオニックな反応になってしまったことがありました。
三重項として計算をすすめて構造を見つけた後に、その構造を出発点に、もう一度一重項として計算を行う方法があるようです。とんでもないことをやっているようにも聞こえますが、αスピン、βスピンの距離が離れた状態だと、βスピンがひっくり返って三重項となった状態と、あまりエネルギーは変わらないようです。これを利用すると、閉殻的な、ニセの一重項にトラップされることなく、望みの電子配置に極めて近い構造、軌道、から計算を始められるようです。

（新たに求まった、開核性の一重項の電子配置にある遷移状態が、閉殻性の一重項よりも不安定だったら、、どうしたらいいんでしょうね。最適な汎関数を探す旅に出なければならないかもしれません）

Gaussian先生の構造最適化は、他の計算ソフトと比較して、とても優秀だそうです。ですが、その能力を過信してはいけません。（どんなアルゴリズムで最適化を進めているか、理解してませんよね？私もですが、、気になるあなたはBernyアルゴリズム、について検索すると幸せになれるかもしれません。）
他に、より安定な解があるのに、全然違う準安定な解を拾ってくる可能性については、常に気を配らなければなりません。
結合が安定にできる原子間距離にあるときの計算は、結構精度よく再現できているそうですが、原子間距離が離れてくると、誤差が大きくなってきます。
DFT法ではなく、HF法の例ですみませんが、名古屋工大和佐田先生の手引に良い図がありました。図４、図５を見ると分かるように、Ristrictedでの計算では、原子間距離が長くなるほどに大きな誤差が、Unristrictedでの計算では、無限遠に至る途中に、誤差が大きくなる部分があります。水素分子の解離ですらこれですから、より複雑な分子の計算は押して知るべしです。


[一度、エネルギー計算が終わったあとに、謎のエラーで計算が止まってしまう]
遷移状態の構造最適化を行う際、"それらしい"初期構造を作り、その構造での力の定数を計算し、それを基に遷移状態（n-1本の結合については正の振動数をもつが、たった1本の結合については負の振動数を持つ状態）を探索していきます。
状況について、私自身もあまり整理ができてませんが、他の計算（例えば、より単純な構造ではじめの計算を行う）で作った構造を流用して計算をしたりすると、初期構造についての計算のあと、以下のようなエラーを吐いてパタリと計算が止まる事例に悩まされていました。

Error termination request processed by link 9999.

Error termination via Lnk1e in /usr/local/g09/l9999.exe at Wed May 21 19:31:46 2014.

おなじ悩みを抱えた人もいるようで、検索していくとResearch Gateでお悩み相談している人がいました。(https://www.researchgate.net/post/Can_anyone_explain_a_consistent_link_9999_error_during_TS_with_Gaussian)

ここで、opt=(TS,calcfc) となっているところを、opt=(ts,noeigentest,calcfc) とするとよいのでは、というアドバイスがあり、私も試してみると、きちんと計算が回り始めました。この、構造最適化のオプションキーワード、"noeigentest"についてはHPCシステムズのページに、本家ガウシアンの解説ページの翻訳版があるのでそちらも参照してください。
かいつまむと、構造最適化を進める際に、反応座標のどんな位置にいるのか確認しながら効率的に計算を進めており、変なところにいるぞ！となったら計算を止めてしまうようです。が、分子量がおおきくなってくると、なかなかこのチェックの有用性が薄れてくるようで、公式のページでは計算のサイズが大きい場合に、Noeigentestのキーワードを利用することが推奨されています。HPCさんのサイトでは、Noeigentestで計算する場合はopt=(Maxcyc=5) 程度で、一旦計算を止めてあられも無い方角に、構造がぶっ飛んでいっていないことを確認しては？と、優しい提案をしてくれています。


[RMS Gradient Normが振動して収束の兆しが見えない！？]
再び、名古屋工業大和佐田先生の解説からの引用となりますが、単純な構造最適化でも起こることはありますが、遷移状態の探索では起こりやすいそうです。

• OPT=(TS,calcall) で、構造を変えるごとに力の定数を計算させる
• OPT=(TS,maxstep=5) で、一度に変形させる結合長、結合角を小さくする。

などが、提案されています。１つ目については、一般に、力の定数の計算は結構重たいので、計算のサイズと要相談のオプションです。
また２つ目のmaxstepですが、何も指定していないと(デフォルト値は)、maxstep=30 となっているそうです（元の解説では、IOP(1/8=5)のように指定できると説明されています）。もちろん、収束点から遠い場合にmaxstepを小さくすると、余分な計算が増えます。むやみに小さくするのは得策ではありません。

計算を回し続け、計算機資源、電気代を無駄につかうことは一旦やめて、計算の条件を少し変えてやったほうが良いかもしれません。
再び計算を始める際には、guess=read(guess=checkpointと同義), geom=checkpoint を指定し、先程まで使っていたchkファイルから情報を取り出しましょう（ファイルの名前をかえたらいけません）。


Enjoy、遷移状態探索！