化学反応の成否が分子軌道計算で予測できるようになって三十数年が過ぎた．基底状態の分子の計算は容易であるが，遷移状態構造を求めるのは簡単ではない．そのために，パソコン本体より高価な専用ソフトを買ったという話を聞いたことがある．本稿では，フリーウェアを用いて遷移状態を求める方法について，O-allyl S-methyl dithiocarbonateの{3,3]-sigmatropy反応を例にとり，説明したい．具体的には，下図に於いて=CH2とS=C間で新たな結合が生成すると同時にC-O結合が切れる反応である．なお，分子軌道計算にはライセンスフリーのMOPAC2016を用いた．

MarvinSketchで二次元構造を描画

AvogadroでMarvinSketchファイルを開く

三次元構造ではない旨の警告が出るが，OKを押すと水素原子が付加された構造（共平面型）が表示される．原子に続き番号が表示されない場合は左側のメニューのラベルをONにする．

MOPAC入力データの生成

画面上部のExtentionsを開きMOPACを選択，MOPACインプットのFormatがデカルト座標（x,y,z) になっている場合は，ℤマトリックス（距離，角度，二面角）に変更する．

プレビューで確認し，生成を押し保存する．

ファイルの内容 （適当なテキストエディターで開く）

AUX LARGE CHARGE=0 SINGLET PM6

AOXASY

C 0.000000 1 0.000000 1 0.000000 1 0 0 0

C 1.340370 1 0.000000 1 0.000000 1 1 0 0

C 1.498610 1 123.364860 1 0.000000 1 2 1 0

O 1.443620 1 107.830791 1 180.037865 1 3 2 1

C 1.366873 1 123.152640 1 180.041160 1 4 3 2

S 1.676087 1 121.298443 1 359.918996 1 5 4 3 👈 変更個所

S 1.765786 1 113.006666 1 179.915522 1 5 4 3

C 1.801964 1 104.490613 1 179.961873 1 7 5 4

H 1.086113 1 120.814990 1 180.019453 1 1 2 3

H 1.085744 1 121.682069 1 0.015663 1 1 2 3

H 1.087530 1 119.087738 1 179.986526 1 2 1 3

H 1.094973 1 109.366200 1 61.265161 1 3 2 1

H 1.095369 1 109.336076 1 298.762845 1 3 2 1

H 1.093882 1 108.774770 1 180.005699 1 8 7 5

H 1.093456 1 111.216764 1 61.564329 1 8 7 5

H 1.093453 1 111.213330 1 298.446755 1 8 7 5

C=Sの硫黄原子の位置は３→４→５→６で作られているので，３→２→１→６によって作られるように変更する．そのためにはAvogadroの3次元構造から距離，角度，二面角を読みだす必要がある．

反応座標の設定

チオカルボニルの硫黄原子と二重結合の末端炭素原子の距離を少しづつ近づけながら各点に於いて構造最適化を行い，その生成熱を求める必要がある．まず，そのための初期構造を求める．Avogadro→Window→Toolbars→Toolsと移動し，赤←⇢を選択すると距離，角度，二面角が測定できるようになる．原子を二個選択すると原子間距離，３個選ぶと角度，４個で二面角が表示される．

AUX LARGE step=-0.1 point=40 SINGLET precise PM6 👈 書き換える

AOXASY

C 0.000000 1 0.000000 1 0.000000 1 0 0 0

C 1.339977 1 0.000000 1 0.000000 1 1 0 0

C 1.498243 1 123.352214 1 0.000000 1 2 1 0

O 1.443759 1 107.853634 1 180.175327 1 3 2 1

C 1.366667 1 123.227458 1 180.061596 1 4 3 2

S 5.486 -1 38.6 1 0.0 1 1 2 3 👈 書き替え

S 1.765481 1 113.034917 1 179.917283 1 5 4 3

C 1.801891 1 104.469350 1 179.978192 1 7 5 4

H 1.085698 1 120.839275 1 179.999882 1 1 2 3

H 1.086135 1 121.678703 1 359.988906 1 1 2 3

H 1.087707 1 119.098095 1 180.003869 1 2 1 3

H 1.095147 1 109.378674 1 61.444777 1 3 2 1

H 1.095137 1 109.393917 1 298.905647 1 3 2 1

H 1.093810 1 108.764962 1 179.996132 1 8 7 5

H 1.093464 1 111.217051 1 61.558308 1 8 7 5

H 1.093435 1 111.215010 1 298.435634 1 8 7 5

精密化フラグを−１と指定すると，5.486，5.386，5.286，5.186，5.086・・・・とstep幅の-0.1づつ距離を減らしながら，point数だけ構造最適化を行う．

二面角がゼロは不自然であるが，計算中に修正される．中途半端に収束するのを避けるため，preciseを指定して収束判定を厳しくした．

MOPAC計算（フリーソフト）

入力ファイルをドラッグ＆ドロップするとMOPAC計算が実行されて，４個のファイルが生成する．結果はoutあるいはxyzファイルで数値変化を追えば計算の成否が確認できる．outファイルの最後に各距離に対するエネルギーがまとめられた表が表示されるが見にくいので，xyzファイルを使用したアニメーションを勧めたい．MOPACについては別稿を参照．

アニメーション

Avogadroでxyzファイルを読み込めば，構造変化を動画で確認することができる．

ファイルを選択し，動的結合をオンにして再生▶を押すと６員環構造を経て生成物への変化が表示される．以下に５点の構造を画像カルーセルで示した．

STEP計算の結果はarcファイルに纏められているので，そのデータを表計算ソフトで読み込めば生成熱の変化を図示することができる．

抽出結果はGoogleのフリーアプリ（Google表計算）でグラフ化する．頂点がTS付近．黒文字の部分は共平面型から屈曲型に変化する過程のため使用しなかった．プロットしたのは青字（赤字を含む）の部分．

TS計算

頂点の座標を使用して，PM6, TS, XYZ, PRECISEを用いて精密化する．その際，すべての原子の精密可フラグは１にする．XYZを指定したのはデカルト座標の方がTS計算に適しているという経験上の理由によるものである．

TS構造の座標データはarcファイルから必要な部分を切り出せばよい．それを使って，FORCE計算，次いでIRC計算（PM6 irc=-1 LARGE=50）を行い，TS構造が出発物質，生成物に繋がっていることを確認する．IRC計算結果のxyzファイルはAvogadroのアニメーションで表示させることができる． irc=-1の指定で，TS→前駆体． irc=1の指定で，TS→転位体への構造変化に対応する結果を与える．

まとめ

Cope転位，Claisen転位反応などの遷移状態についても同様な手法で計算可能である．一連の計算や解析を統括した分子計算アプリ（Winmostar 等）が市販されているが，たいへん高価であり，リタイアした人間には手が出せない．しかし，計算の仕組みの概略が分かっていればフリーソフトで実行可能である．Winmostarのようなアプリが存在しなかった遥か昔，組み立て式分子構造模型と分度器を使ってZマトリックスを作成していた頃を思い出しながら，後期高齢者の備忘録として記した次第である．

フリーソフトのインストールについて

Macの場合， 分子計算関連ソフトウェアをダウンロードしてインストールしようとすると， 「開発元を検証できないため開けません」のメッセージが表示され， 先に進めないことが多い． そのような場合は， システム環境設定でセキュリティ設定を変える必要がある．

参考資料

64ビットMacOSのMOPAC計算

Marvin Sketchの3D化機能

Marvin Sketchによる分子描画 Ⅱ

Marvin Sketchによる分子描画 Ⅰ

(2022.07.17)