基本的な反応機構の計算方法 [Gaussian]
私は学生時代が理論化学計算の研究室だったということもあり、
計算化学はとっても馴染み深いです。
学生さんから質問された事について、知識を整理するためにnoteに纏めます。
なるべく無料記事で書いていきたいと考えてますが、
難しい内容のものは一部有料にする場合があります。
ご了承ください。
反応機構の探索へのアプローチ
反応機構の解析を計算で行う場合、遷移状態（TS；Transition State）の探索をする必要があります。
遷移状態とはなんぞや、というレベルの方は別の記事に纏めておくので
そちらを読んでおいてください。
遷移状態の探索と言うと、必ず振動計算を行う必要があり、虚の振動が唯一でかつその虚振動の分子の動きが正しい必要があります。ものによってはすごく難しかったりします。
このTS構造の探索のアプローチ方法はよく用いられるもので3つあります。
以下に、難易度の順に示します。
①既に計算済の似た構造を用いて構造を作成する（TSモチーフ法）
②Minimun Energy Path（Minp）計算
③QST3計算
これらについてサンプルと一緒に解説していきたいと思います。
①TSモチーフ法[1]
参考になる論文（置換基違いとか）や研究室の先輩からデータを引き継いで計算を行う場合は、この方法を使うと遷移状態の計算の難易度がかなり下がります。
置換基効果等の影響を強く受ける場合を除いて、TS構造の反応部位の原子間距離はだいたいの値が決まっています。これを利用して、既にあるTS構造を鋳型として、置換基をターゲットの構造になるように変更して計算を行います。
計算を行う時には、反応部位の原子は固定しないと、全体のエネルギーが安定する方向にoptが進んでしまうので、気を付けてください。私は、なんだかんだで反応部位の原子は2～7個程度あるので、座標固定をよく用います。
計算を実際に行ってみて、TS構造に落とし込めない場合は、なるべくモチーフから変更する置換基の数を少なくして（なるべくモチーフに近い構造で）行うと上手くいくかもしれません。
サンプル1：ethyleneとbutadieneのDiels-Alder反応 [input]
サンプル2：butadieneに-OH基を付加したDiels-Alder反応 [input]
▼サンプル1
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample1.chk
# b3lyp/6-31g(d) opt=modredundant 6d optcyc=200 nosymm
Freeze Optimization
0 1
C 0.42739300 1.43895100 0.49122900
C 1.31628100 0.71368600 -0.28825400
C 1.32837700 -0.69329900 -0.28937400
C 0.45364400 -1.43558500 0.48964100
C -1.56507000 -0.70585400 -0.22840100
C -1.57334700 0.68510700 -0.22885300
H 0.11547800 1.06851000 1.46067500
H 0.36257700 2.51758000 0.37577300
H 1.86232300 1.22998200 -1.07560300
H 1.88345500 -1.19858400 -1.07757300
H 0.13462500 -1.07274400 1.45960300
H 0.40699600 -2.51482300 0.37127200
H -2.06494400 -1.24768000 0.56887000
H -1.46255000 -1.24676000 -1.16196500
H -1.47887900 1.22564700 -1.16353900
H -2.08273900 1.22083600 0.56656300
X 1 F
X 6 F
X 5 F
X 4 F
--Link1--
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample1.chk
# b3lyp/6-31g(d) 6d opt=(ts,noeigentest,calcfc,nofreeze) freq geom=check
Nonfreezed optimization
0 1▼サンプル2
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample2.chk
# b3lyp/6-31g(d) opt=modredundant 6d optcyc=200 nosymm
Freeze Optimization
0 1
C 0.42739300 1.43895100 0.49122900
C 1.31628100 0.71368600 -0.28825400
C 1.32837700 -0.69329900 -0.28937400
C 0.45364400 -1.43558500 0.48964100
C -1.56507000 -0.70585400 -0.22840100
C -1.57334700 0.68510700 -0.22885300
H 0.11547800 1.06851000 1.46067500
H 0.36257700 2.51758000 0.37577300
H 1.88345500 -1.19858400 -1.07757300
H 0.13462500 -1.07274400 1.45960300
H 0.40699600 -2.51482300 0.37127200
H -2.06494400 -1.24768000 0.56887000
H -1.46255000 -1.24676000 -1.16196500
H -1.47887900 1.22564700 -1.16353900
H -2.08273900 1.22083600 0.56656300
O 2.03369285 1.39201636 -1.32270465
H 2.81104516 1.81418729 -0.94975012
X 4 F
X 5 F
X 6 F
X 1 F
--Link1--
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample2.chk
# b3lyp/6-31g(d) 6d opt=(ts,noeigentest,calcfc,nofreeze) freq geom=check
Nonfreezed optimization
0 1②Minp計算
TSモチーフ法が使えない場合や、置換基効果等の作用によりTSモチーフ法で上手くTS探索できなかった場合に使用することが多いです。
変数に使用できるのは、距離、角度、二面角です。
変数の数は、距離の場合は2つ、角度は3つ、二面角は4つです。
これらの変数を変化させて、エネルギーを計算して、そのエネルギーの山（最も高くなっている点）からTS候補構造を探索します。このTS候補構造についいて、TS計算をした場合に虚の振動がない場合はTS構造とは言えないので（TS構造は虚の振動は唯一だから）、変数の変化量の幅をより小さくして探索する必要があります。
サンプル3：ethyleneとbutadieneのDiels-Alder反応のminp計算 [input]
▼サンプル3
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample3.chk
# b3lyp/6-31g(d) opt=modredundant 6d optcyc=200 nosymm
Minimum energy path calculation
0 1
C 0.42739300 1.43895100 0.49122900
C 1.31628100 0.71368600 -0.28825400
C 1.32837700 -0.69329900 -0.28937400
C 0.45364400 -1.43558500 0.48964100
C -1.56507000 -0.70585400 -0.22840100
C -1.57334700 0.68510700 -0.22885300
H 0.11547800 1.06851000 1.46067500
H 0.36257700 2.51758000 0.37577300
H 1.86232300 1.22998200 -1.07560300
H 1.88345500 -1.19858400 -1.07757300
H 0.13462500 -1.07274400 1.45960300
H 0.40699600 -2.51482300 0.37127200
H -2.06494400 -1.24768000 0.56887000
H -1.46255000 -1.24676000 -1.16196500
H -1.47887900 1.22564700 -1.16353900
H -2.08273900 1.22083600 0.56656300
B 4 5 S 20 0.1このサンプルは、4行目の原子を5行目の原子に0.1ずつ20ステップ近付けるというものになります。
③QST3計算
上記2つの方法でもTS構造が出せなかった場合に使用します。この方法では、(1) TS候補構造の座標だけでなく、(2) 反応物と (3) 生成物の座標も使用するため、そこそこに面倒です。
この計算のTS候補構造は、極めてTS構造に近い必要があります。じゃないと、余裕でエラーになるか自分が目的としているTS構造とは異なる挙動をする構造ができてしまいます。
学生時代に私も使ったことはありますが、なかなか上手くいかなかった記憶があります。確か上手く構造出すまで数日かかりました（①や②だとだいたい1～2日で出せます）。
そのため、難易度は上の2つよりも高いということにさせてもらいました。
QST2だと反応物と生成物の構造の2つからTS構造を求めることができる、とGaussianのマニュアル等読んでると書いてますが、まあ全然上手く出ないです。もう1構造与えたQST3の方がTSに近い構造から探索されるので、計算コストの面でも算出しやすいのかなと思います。
サンプル4：ethyleneとbutadieneのDiels-Alder反応のQST3計算 [input]
▼サンプル4
%mem=[メモリー数]
%nproc=[proc数]
%chk=sample4.chk
# b3lyp/6-31g(d) opt=qst3 6d scf=noincore optcyc=200 nosymm
reactant
0 1
C 1.36221 -1.5949 -0.61545
C 1.63921 -0.70868 0.35071
C 1.69667 0.75159 0.18774
C 0.9694 1.47342 -0.67529
C -2.69449 0.68165 0.56121
C -2.59731 -0.60382 0.22782
H 1.18039 -1.28547 -1.64166
H 1.32574 -2.66137 -0.41243
H 1.87072 -1.08177 1.34897
H 2.38725 1.27269 0.85176
H 0.22364 1.01612 -1.32048
H 1.08775 2.55054 -0.7518
H -3.65438 1.19025 0.61804
H -1.81431 1.27776 0.78948
H -1.63478 -1.10684 0.17317
H -3.4743 -1.2046 -0.00324
product
0 1
C 0.11709 -1.40009 -0.34168
C 1.2605 -0.64861 0.29417
C 1.23932 0.68811 0.29387
C 0.07315 1.40275 -0.343
C -1.27481 0.75782 0.06666
C -1.25007 -0.79737 0.06869
H 0.21556 -1.3504 -1.43784
H 0.14599 -2.46492 -0.08418
H 2.08559 -1.20292 0.73718
H 2.04653 1.26843 0.7365
H 0.1738 1.35527 -1.43905
H 0.06842 2.46817 -0.0864
H -2.06495 1.12864 -0.59701
H -1.52734 1.11248 1.07233
H -1.489 -1.15727 1.07583
H -2.02929 -1.19493 -0.5924
Candidate Structure of TS
0 1
C 0.46949 -1.42985 -0.4884
C 1.33318 -0.68241 0.29148
C 1.31117 0.7246 0.29084
C 0.425 1.44411 -0.49003
C -1.58822 0.66825 0.2247
C -1.56572 -0.71753 0.22648
H 0.14434 -1.06659 -1.45662
H 0.43346 -2.51087 -0.37731
H 1.88812 -1.18323 1.08354
H 1.85 1.24326 1.08254
H 0.11128 1.07008 -1.45792
H 0.35505 2.52355 -0.37973
H -2.09471 1.20194 -0.57445
H -1.49509 1.21405 1.15685
H -1.45463 -1.25751 1.16006
H -2.05516 -1.26953 -0.57086 QST3計算の結果（*.log）は、座標の保存ができないので
GaussViewで表示させた時に構造全体をコピーして
新規作成したWindowにペーストすると座標が抜けます。
この結果を①のTSモチーフ法と同様にTSモチーフを固定してTS計算をします。
④その他
参考までに上記の方法以外によく知られている手法も書いておきます。
■GRRM
前田先生、大野先生らによって開発されたプログラムで、自動で遷移状態や反応経路を導ける画期的なものです。私は使った事がないので解説できない・・・ごめんなさい(´・ω・｀)
興味がある方は検索してみてください。
■エネルギー等高線図法
反応に強く関係すると思われる2変数を選択して、片方の値を固定して、もう一方の値を変化させた計算からエネルギーの山を探索する手法。やり方としては、minp計算は1つの変数を変化させているところ、この手法では2つの変数について行えばOK。opt=modredundantで座標後の改行後に入力する変数を1つ増やして、固定する変数をfreezeさせて行います。難易度は高くないけど、計算量が多くなります。
■半経験的MO法の利用
私自身使ったことがないので、割愛させていただきます。詳しくは、堀先生の著書「Gaussianプログラムで学ぶ 情報化学・計算化学実験」に記載があります。
最後に
計算化学をやってみたいけど、わからない事だらけ！って方は↓こちらのテキストがおすすめです。
講義形式になっているので、ゼロから学ぶ場合に助けになります。
ある程度の知識がある場合や、↑のテキストの補助教材として購入するのであれば↓こちらがおすすめです。
カラーでとてもわかりやすく解説されています。
ページ下部のハートアイコンの「いいね」を押していただけると、励みになります。
引用
[1] 前山恵璃、山口徹、隅本倫徳、堀憲次、"最安定反応経路探索（MSRP）法開発とPinner Pyrimidine反応への適用"、JCAC、(2021).