Diels-Alder反応をシミュレーション！量子化学計算ソフトで遊んでみよう-超初心者向け完全無料の量子化学計算入門２

第一回講座→計算ソフト導入編-超初心者向け完全無料の量子化学計算入門１
第二回講座→これ

はじめに

どもども
量子化学計算をゼロから勉強する連載記事第二弾です。

さっそく、前回の記事の内容ですが量子化学計算の流れをおさらいしましょう。

「モデリングソフト」で分子を組み立てる
↓
「プリ・ポストプロセッサ」で計算の細かい設定を行う。インプットファイルを生成。
↓
「計算ソフト」で実際に計算
↓
「可視化ソフト」で計算結果を表示

という4段階の流れでした。前回、計算の準備のために無料版のWinmostarを導入しました。Winmostar自体はモデリング・インプットファイル生成・計算結果の可視化ができます。また、計算ソフトとしてMOPACというソフトを同梱しています。

つまり、上記の量子化学計算のすべてのステップをWinmostar上で行えるということでした。便利ですね。

今回のゴール

今回は、実際にWinmostarを弄って遊んでみようと思います。数式使うような本格的な量子化学の内容はないです（笑）。
今回のゴールは以下のようなDiels-Alder反応のgifアニメーションを作ることです。
緑色の原子は炭素、黄色の原子は水素を表しています。アニメの上の”Length”というところは、C-C結合ができる部分の炭素間の距離を示しています。

ぬるぬる動いてめっちゃ気持ちいいｗ
反応のシミュレーションって結構雲の上の話だと思っていたのですが、思っていたより簡単にできるしめっちゃたのしいです
(* ´ ▽ ` *)

まずはお勉強

反応計算の流れ

１．反応前Aと反応後Bの、分子の安定な構造を計算(構造最適化)

２．反応の途中の構造TSを計算。(TS計算)

２．A→TS→Bという反応の流れをつなげてぬるぬるアニメーションにする。(IRC計算など)

1. 構造最適化(安定な構造を探す計算)

分子の安定な構造とは

私たちは紙の上で分子を表現しようとしたとき、2次元の構造式を書きます。構造式は分子の大まかな形が分かって便利なのですが、実際の原子は3次元の、立体的な構造を持っています。

じゃあ、実際の分子の3D構造ってどうなっているのでしょうか。たとえば原子と原子の距離や角度ってどのくらいか決まっているのでしょうか。

実は、現実の分子が取りやすい構造は決まっており、計算できます。というか今回はこれをまず計算します。
ポイントは分子の持つエネルギーに着目することです。分子の持つエネルギーはポテンシャルエネルギーとか内部エネルギーとか言います。分子は一般的に、今の構造よりもより低いポテンシャルエネルギーになるように構造を変化してさせていきます。最も低いエネルギーになったとき、構造を変化させてもエネルギーが高くなるしかないので、その構造に落ち着きます。こういう構造を安定な構造といい、これが分子がもっとも取りやすい構造ということです。
注意点として、安定になるようにいろいろ動いて構造が変化した結果、安定な構造に落ち着く可能性が極めて高いよねってだけで、別に分子が止まっているとか安定な構造以外の構造になりえないとかそういうわけではないです。分子は常に振動してます。

たとえばメタン(CH4)は正四面体の構造をとるといわれています。正四面体だと分子のエネルギーが最も低くなるような構造になりやすいです。

大学の化学に触れる機会がなかった方は、どうやったら安定な構造になるのか、つまりポテンシャルエネルギーを変化させる要因がなんなのかをイメージしにくいと思います。一番わかりやすい要因は”立体障害”です。

原子は原子核の周りに電子が広がっているような形をしている、という図を見たことがあるかと思います。、電子はマイナスの電荷をもっており、マイナスとマイナスは反発します。反発ってなんだよと思った方は、よくある磁石の反発をイメージしてください。このような反発のせいで、かさばった原子や分子同士は反発するためあまりに近づこうとはしません。逆に言うと、分子はあまりかさばらないような構造になろうとします。

（もちろん結合を保ったままです。究極にかさばらないのは原子同士が無限に離れているときです。当たり前ですがそんなのどうでもいいです（笑））

先ほど例に出したCH4の話を思い出してください。中学校の数学の幾何でやった方もいるかもしれませんが、正四面体の頂点は、中心と頂点の距離を決めたとき、4つの頂点が最も離れるような形になっています。つまり、一番かさばらない立体ということです。だから、H原子同士の距離が最も広がるのが正四面体構造なので、現実でもそのような形をとっていると考えられます。

“安定な構造”の分子のエネルギーを計算する方法 –MM法と第一原理計算

ということで、上の説明を読んで察しのいい方は
「これから立体障害をうまいこと計算のパラメータにして、エネルギーが低い安定な構造を探す計算を行うんだな！」
と思ったかもしれません。素晴らしい発想だと思います。でもごめんなさい。今回は違います。上の話は安定な構造を頭の中でざっくりとイメージする方法の話です。

確かに、立体障害をはじめとして、結合距離・結合角度・二面角（要するに面と面の角度）というイメージしやすいパラメータを使って分子全体のポテンシャルエネルギーを計算する方法はちゃんとあります。MM法と呼ばれる方法です。前回の記事でAvogadroを導入してみた方は、ボタンをぽちっと押すだけでめっちゃ簡単にMM法で遊べます。（下にgifアニメあとでいれる）

ＭＭ法は簡単ですし、計算量もめちゃくちゃ少ないのですぐに終わります。生物系のような計算をできるだけ少なくしたい分野の方だとなじみ深い計算方法かもですね。ネットで検索したこの記事とかけっこうおもしろかったです。でも、精度がいいとは言えません。またMM法のWikipediaにも欠点として書いていますが、前もって準備するパラメータが結構必要ですし、計算結果もそれに大きく依存します。論文でなんにも考えずに使ってたら怒られるらしいです。

ということで、もっと別の精度の高い方法が必要です。実は昔の偉い人たちが量子化学ででてくる式をこねこねしてもっと高精度で分子のエネルギーを求める方法を発明しています。第一原理計算という手法です。・・・一行で書きましたがどうやってそういう発明にたどり着いたのか私にはさっぱり理解できません。この分野勉強しようと思うと、「本当に誰が思いついたんだよこれ・・・」って思うような我々の理解を超越する内容が大量にでてくるので覚悟してください。歴史的にはたぶん相当面白い歴史があるんだと思います（笑）。

それはさておき、今回行う量子化学計算はまさにこれです。昔の偉い人が考えたことをコンピュータ使って遊んでみましょう。

分子のエネルギー求める量子化学の式ってなんなのよ

一応量子化学計算入門の記事のつもりなので、量子化学の式についてもちょっと触れときます。このサイトがわかりやすいです。

シュレディンガー方程式HΦ＝EΦという難しい等式が出てくるのですが、これがなんなのかは今はわからなくてもいいです。とにかく、この式の中にはエネルギーEが含まれているので、うまいこと解いてやることでEの数値が求まるというわけです。このようにシュレディンガー方程式を解くことを第一原理計算といいます。（第一原理計算は正確には、その手法の一つという意見もあるかと思いますが、人によって結構バラバラです・・・。シュレディンガー方程式を解くほかの表現方法がわからなかったのでとりあえず第一原理計算と書きます。もっと正しい言い回しあったらぜひ教えてください。）

”うまいこと”解くと書きましたが、ここに非常に多くのテクニックがあるわけです。計算時間が短い代わりに精度が低い結果にする計算方法や、計算時間は長いが精度がいい計算方法など、調べれば大量に出てきます。
このサイトに書かれている第一原理計算の分類は一応目を通しておいたほうがいいと思います。

2. TS計算

反応の途中の構造とは

反応とは、かさばらないように離れている(安定な構造をとっている)２つの分子が近づいて、くっつくことです。

くっつく途中は、どうしても分子同士は歪んでしまうため、反応前に比べて、反応途中の構造は不安定、つまりエネルギーが高い状態なんですね。

反応計算は、このくっつく途中の、エネルギーが高い状態の構造も計算しておく必要があります。

反応途中において、エネルギーが最も高い状態を日本語で遷移状態、英語でTS(Transition State)といいます。

くっついた後は、分子の歪みがなくなり、TSよりも安定になります。

不安定＝TS　じゃない！！

さて、ここでとても重要なことがあります。

「TS=最も不安定な状態」ではないということです。

「TS＝”反応の途中で”最も不安定な状態」です。

反応の途中というのがポイントです。

この制限がなかったら、原子核を重ねたとき最高に不安定になります。

そりゃ不安定です。２つの球が重なっているようなものです。核と核が重なったら核融合です。

つまり、この”反応の途中”という制限があるため、TS計算はTSぽい構造を手で入力する必要があります。これはセンスが問われます。

3. IRC計算

本当にTS構造なのかを確かめる計算

さて、TSは”反応の途中”であるということがポイントだと書きました。

TS計算が終わった後に、計算結果が本当に反応の途中の構造を表しているか確証を取りたいですよね？

これを行うのがIRC計算です。TS構造から安定な構造にたどっていけるかを確かめる計算です。

IRC計算に成功すると、安定構造を知ることができます。

手順

1.　反応前の構造を書いて、構造最適化

反応前を書きます

2.　反応後の構造を書いて、構造最適化

3．　TS計算

4.　IRC計算

5.　アニメーション作成

/