ORCAによる反応解析の基本

ここでは、ORCAを使った反応解析のやり方について、簡単にまとめておきたいと思います。
題材としては、座標が組みやすくシンプルな反応ということで、マロンジアルデヒド(エノラート型)の分子内プロトン移動を取り上げます(下図)

Fig.1

1.安定構造の最適化

まず安定構造の最適化から始めます。上図の番号に沿ってGaussianタイプのz-matrixを組むと、例えば以下のような入力ファイルになります。

ORCAのキーワード行は「!」で始まりますが、Gaussianと同様に整理されたキーワードで完結に計算内容を指定できます。ここでは「Opt」が構造最適化、「Freq」が(解析的な)振動解析、「B3LYP」がDFT計算の汎関数指定、「6-31G(d)」が基底関数の指定です。B3LYP/6-31G(d)というモデル化学で、構造最適化と(最適化構造の)振動解析を行うという意味になります。
構造データは「*」で挟まれた部分になりますが、「gzmt」はGaussian形式のz-matrixを示し、「0」は系の電荷、「1」は系の多重度を示します。閉殻系で中性の分子ですので、このような表記になります。
ちなみに「#」はコメント行を表します。

計算を実行すると、7ステップで構造最適化が終了します。最適化された座標はxyz形式で出力されます。

振動解析の結果は出力ファイルの後半に記録されています。

虚振動(負の振動数)が見られないことから、安定構造であると推察されます。
ファイル末尾には各種熱力学データが出力されており、デフォルトでは298.15K (25 ℃)での計算値が示されます。全電子エネルギー(Electronic energy)やゼロ点エネルギー(Zero point energy)、Gibbsの自由エネルギー(ORCAではGibbs free enthalpyと表記)といったよく使われる値はここから見つけることができます。

 

2.遷移状態の探索①

得られた安定構造をベースに、遷移状態を探索していきたいと思います。OHのプロトンが反対側のカルボニル酸素に移動していく反応ですので、この辺の結合距離を変化させて、エネルギーが高いところを探索するのがわかりやすい方法です。
これを実現する入力ファイルは、例えば以下のような形になります。

キーワード行は構造最適化と全く同じですが、下に「%geom」「end」で挟まれたブロックが挿入されています。「%geom」は構造最適化計算の細かい設定をするブロックで。更に「scan」「end」というサブブロックで、特定の構造パラメータ(結合長や結合角、二面角)を変化させながら構造最適化を行うことを指定します。「B 4 8」は4番と8番の原子間距離を変化させるという意味で、「1.65, 0.95, 8」は1.65 Åから0.95 Åまで8 stepで変化させることを意味します(つまり1.65, 1.55, 1.45, 1.35, 1.25, 1.15, 1.05, 0.95という風に0.1 Åずつ縮めます)。結合角なら「A 4 8 2」、二面角なら「D 4 8 2 1」という感じで指定します。
ちなみに、ORCAでは原子を指定する番号が「0」から始まりますので注意が必要です。xyz座標で一番上にある炭素が0番、一番下の水素が8番です。

計算は構造最適化を8回繰り返すため、少し時間がかかります。終わると、各ステップの最適化構造が全てxyz形式で出力されています。logファイルの末尾を見ると、どこが一番エネルギーが高いかわかります。

これを見るに、最もエネルギーが高いのは1.25 Åの構造です。つまり、この構造が一番遷移状態に近い可能性があります。上から5番目ですので、「ファイル名.005.xyz」がその構造のxyzファイルです。

 

3.遷移状態の探索②

得られた候補構造を基に、遷移状態の最適化を行います。典型的な入力ファイルは以下のような形です。

最適化のキーワードは「OptTS」で、%geomブロックで特殊な指定をします。「Calc_Hess」は構造最適化前に振動計算を行うかどうか、「Hybrid_Hess」はいくつかの原子を指定して正確な振動解析を行い、他は推算で補うこと(これをHybrid Hessianと呼びます)を指定します。Hybrid Hessianでは、遷移状態に深く関与する原子を指定することが重要で、ここでは2つの酸素原子と挟まれた水素原子をすれば良いと思います。

振動解析を先にしなくてもOptTS計算はできるのですが、遷移状態は不安定な状態ですので、正確な道標が無いと安定構造に落ちていってしまうことも多いです。Hybrid Hessianで計算コストに配慮しつつ精度の高い道標を用意することで、効率よく遷移状態最適化ができるようになります。今回のように小さい分子であれば、Hybrid_Hessの行を削除して、全原子の振動計算をしても構わないでしょう。

計算が終わると、通常の最適化と同様に構造がxyz形式で出力されます。

logファイルで振動解析結果を確認すると、一つだけ負の振動数があります。

これで、確かに遷移状態であることがわかります。しかしながら、求めているプロトン移動の遷移状Avogadroで振動可視化態かどうかは、振動数を見たところでわかりません。可視化にはいくつか方法がありますが、無償ツールだけで行うならば、「orca_pltvib」というORCAの同梱ツールで虚振動のアニメーションを20フレーム分xyz座標に変換し、これをAvogadroで可視化することで確認することができます。有償ツールであれば、ChemCraftがとても簡便です。
右に、Avogadroで振動解析を表示した際の様子を示しました。ここに至る手順を完結に記載すると、①hessファイルからorca_pltvibで虚振動のアニメーションデータを切り出す(コマンドラインで「orca_pltvib <hessファイル> 6」) ②出力された006.xyzファイルをAvogadroで開き、Extensions→Animationで確認 の2ステップで行います。

見事に、プロトン移動の遷移状態であることが確認できると思います。

 

4.IRCの確認

単純な遷移状態であれば、手順3までで十分解析が出来ていると言えなくもないのですが、複雑な反応になっていくと、得られた遷移状態から安定構造に向けて峠をそーっと下っていくと、思ってもみない安定構造に行き着く(=求めた遷移状態は、最初に期待したものとは異なる可能性がある)ということがあります。
きちんと想定通りの麓(安定構造)と繋がった遷移状態かどうか、峠をゆっくり下って確かめる計算はIRC(Intrinsic Reaction Coordinate)計算と呼ばれますが、ORCAではこれを指定するキーワードが存在しません。構造最適化計算を使って擬似的に再現することで、ほとんど同じことができます。

手順3の最後に作成したアニメーション用ファイルは以下のような内容です。

数字が6列並んでいますが、左側3列がxyz座標に相当します。また、一番上と一番下のデータは遷移状態そのものに相当します。
この中で、上から2番めと下から2番めのxyzデータを取り出します。これらは、遷移状態からちょっと進んだ構造とちょっと戻った構造に相当します。そして、これらを特定の条件で構造最適化することで、IRCを再現することができます。

%geomブロックの設定ですが、inhess, update, stepはこれで固定します。構造最適化アルゴリズムとしてQuasi-Newton法を使い、構造変化を追跡するためのhessianとしてunit hessianを使い、hessianのアップデートをしない設定にすることで、最急降下法的にエネルギーの山を下っていくことができます。
この計算の最適化過程や最適化構造を見ていくことで、求めた遷移状態が確かに想定していたものなのかが確認されます。

 

本ページで紹介した計算の入出力ファイル:反応解析の基礎(ORCA)

Be the first to comment

Leave a Reply