多数のツール、ファイルの組み合わせで行ないますが、これは優れたユーザーインターフェースを持ったMacのファインダがあってこその芸当ですね。
バッチ処理できるような定型作業は、「研究」ではなく「作業」です。
Chem3Dを起動してみましょう。
上の方に、入力ボックスがあります。
ここに、“CH2CHOH”と入力して、リターンを押します。
すると、図のような構造が表示されるはずです。
これがビニルアルコールです。（以下VAと略）
「command+M」にて、エネルギーの最小化を行なって形を整えておきましょう。
官能基の向きを変えて、どの方向が一番エネルギーが低いか、もれなくチェックしてください。
エネルギーが「最小値」ではなく、「極小値」になっている事もあります。シビアに検討して下さい。
それぞれの原子に、きちんと番号をふって管理しやすくします。
Command+Aですべての原子を選択し、Buildメニューの「Show serial numbers」を選びましょう。
その後、入力BOXに“1”を入れます。そして１番に指定したい原子をクリックすると、それが一番に指定されます。
自動的に次の番号になっていますので、あとは番号順に原子をクリックしていけばいいだけです。
これをMOPAC形式で保存しましょう。ファイルメニューの「save as...」を選び、ファイル形式でMOAPCを選びましょう。それだけです。
拡張子など気にしなくてかまいません。日本語もＯＫ。
同様に、アセトアルデヒドも作ります。
作為的ですが、今回は後付けの解析なのでこれで良しとします。
ここをあなたの「仮説」とすれば、そして検証と考察を怠らなければ、それは立派な「研究」となります。
PM3 EF PRECISE XYZ GEO-OK ←一行目
VA ←二行目
←三行目
C 0.000000 0 0.000000 0 0.000000 0 0 0 0
H 1.101288 1 0.000000 0 0.000000 0 1 0 0
H 1.101608 1 119.340546 1 0.000000 0 1 2 0
C 1.343140 1 120.994659 1 -179.999619 1 1 2 3
H 1.102463 1 120.069687 1 0.000000 1 4 1 2
O 1.522858 1 123.168076 1 179.999619 1 4 1 2
H 0.975693 1 107.994354 1 1.379578 1 6 4 1
一行目には“キーワード”を入れます。これは、MOPACに与える命令です。ここでは説明しきれませんので、各自参考書を見てください。
PM3というのはパラメーターセットを、EFは計算方法、PRECISEは計算精度、XYZで座標系、GEO-OKで異常接近座標の無視、を指定しています。
二行目には、Chem3Dがファイル名を書き込んでおいてくれます。何を書いてもかまわないのですが、一般的にはファイル名です。
三行目も、何を書いてもかまわないのですが、一般にはメモ程度の内容を書き留めておきます。
何も書かなくとも、空けておいてください。
キーワードを指定したら、保存します。ここでも拡張子は必要ありません。
＊ダミー原子を使用した場合、出力ファイルには“XX”という原子名で記入されます。
EF Working.....
CYCLE: 1 TIME: 1.00 TIME LEFT: 99999.0 GRAD.: 130.744 HEAT:-17.96588
CYCLE: 2 TIME: 0.00 TIME LEFT: 99999.0 GRAD.: 117.469 HEAT:-25.10923
CYCLE: 3 TIME: 0.00 TIME LEFT: 99999.0 GRAD.: 73.993 HEAT:-27.98541
CYCLE: 19 TIME: 0.00 TIME LEFT: 99997.0 GRAD.: 0.031 HEAT:-31.56866
COMPUTATION TIME = 3.000 SECONDS
== MOPAC DONE ==
GRAD.というのは、エネルギー勾配です。これが十分小さい（もうエネルギーが落ちない=
安定である）ところまで来た、ということで計算が無事終わりました。
キーワードでPRECISEを指定しているので、0.05以下になるまで計算しています。
（しない場合は5以下。100倍厳しくなるのです、その分、時間もかかります。大きな分子では良く考えてから使いましょう。この場合は約３秒で計算が終わっています。買って良かったPowerMac）
さて、６と１１を開いて、エラーメッセージや不審な挙動が記録されていないかチェックしてください（判断は参考書などを参考にしてください）。
大丈夫でしたか？ならば、MOPACについてきた“MAKENEW”と同じフォルダに、「FOR012」を持っていきましょう。
そしてMAKENEWを起動します。そして出力したいファイルの名前を指定しリターンを押せば（何も指定しないと"NEW"という名前になります）MOPAC形式のファイルができます（Chem3Dでも読めます）。
ファイル書式が分かっている人は、エディタでコピー＆ペーストしてかまいません。
頭の中でZ-Matrixが構築できる人はいいのですが、私を含めた普通の人間はそんなこと出来ません。
＊ダミー原子（“XX”“X”“Cb”など）を使用した場合、Chem3Dに読み込ませる前に、“Du”と直しておく必要があります。
＊Chem3Dは、読み込んだ構造を、自動的に妥当な構造に直してしまうことがあります。
“Windows”メニューから“preferences”を選び、“Building”ウィンドウの中のチェックを“Center Model In Window”以外のチェックを外しておきます。
計算により得られた結果を検討してみます。
ビニルアルコール
-28.721972 Kcal/mol
アセトアルデヒド
-44.199987 Kcal/mol
アセトアルデヒドの方が、圧倒的にエネルギーが低い（＝安定）であることがわかります。
あとは、どのくらいのエネルギーで遷移状態を越えられるかが問題となります。
まず、エネルギー的な中間状態である「鞍点（saddle point）」を求めます。
二つのファイルから、原子の座標と電荷の部分を取り出します（ここはエディタで作業します）。
それぞれの座標部分を、順にペーストします。
その間に、すべて“0（ゼロ）”のカラムを入れます。
そして、キーワードも忘れずに入れます。
今回は、鞍点（saddle point）を求めるので、キーワード“SADDLE”を使います。
“EF”とは共存できないキーワードなので、EFは忘れずに外しておきます（外し忘れるとエラーで指摘される）。
こう言った計算の場合は、XYZ座標を用いるのが定石です。三つ異常の原子が直線上に並ぶと、Z-Matrixでは座標系が記述できなくなりますので、そのような可能性がある場合はわすれずに“XYZ”を使いましょう（普段から使っていても問題はありません）。
PM3 PRECISE XYZ SADDLE GEO-OK
VA-FA
C 0.0000000 0 0.000000 0 0.000000 0 0 0 0 -0.2426
H 1.0843430 1 0.000000 0 0.000000 0 1 0 0 0.0914
H 1.0858579 1 114.876261 1 0.000000 0 1 2 0 0.1020
C 1.3342365 1 121.552683 1 179.996476 1 1 2 3 0.0104
H 1.0957929 1 124.129405 1 0.000000 1 4 1 2 0.0883
O 1.3683073 1 118.487111 1 179.954189 1 4 1 2 -0.2377
H 0.9473633 1 106.740656 1 179.671012 1 6 4 1 0.1882
0 0 0 0 0 0 0 0 0 0 0
C 0.0000000 0 0.000000 0 0.000000 0 0 0 0 -0.1955
H 1.0978936 1 0.000000 0 0.000000 0 1 0 0 0.0622
H 1.0976857 1 107.716117 1 0.000000 0 1 2 0 0.0674
C 1.4985682 1 110.070664 1 124.235985 1 1 2 3 0.2784
H 1.1024178 1 116.655396 1 59.230609 1 4 1 2 0.0407
O 1.2097384 1 123.408302 1 -120.770163 1 4 1 2 -0.3154
H 1.0978891 1 107.538807 1 -115.853611 1 1 2 3 0.0622
上段がビニルアルコール、下段がアセトアルデヒド
MOPACで、同じように計算します。
HEAT OF FORMATION = 27.853705 KCAL
HEAT OF FORMATION = 28.440004 KCAL
上の方が小さいので、こちらの構造（ビニルアルコールから近づいた構造）を使います。
PM3 PRECISE XYZ TS FromVA C 0.0000000 0 0.000000 0 0.000000 0 0 0 0 -0.5592 H 1.0915401 1 0.000000 0 0.000000 0 1 0 0 0.1127 H 1.0918748 1 111.200282 1 0.000000 0 1 2 0 0.1293 C 1.4243168 1 120.815871 1 148.163793 1 1 2 3 0.3269 H 1.0892385 1 129.893323 1 -23.300824 1 4 1 2 0.0902 O 1.2670038 1 109.770395 1 150.312125 1 4 1 2 -0.3781 H 1.4438374 1 80.337494 1 -5.596066 1 6 4 1 0.2783
いつものようにMOPACで計算をかけましょう。
余談ですが、右の図と、上のSADDLEで求めた構造とを比べてください。
炭素−酸素が二重結合から単結合になっています。結合距離が変わっていることがわかります。
MOPACの画面に
RMS GRADIENT = XXXXX IS LESS THAN CUTOFF = 1.00000
GEOMETRY OPTIMISED USING EIGENVECTOR FOLLOWING (EF).
SCF FIELD WAS ACHIEVED
FINAL GEOMETRY OBTAINED
といった表示がされているかがポイントです。この計算の場合、PowerMac6100/80で３秒程度でした。
PM3 XYZ FORCE VibrationAnalysis C 0.0000000 0 0.000000 0 0.000000 0 0 0 0 -0.5929 H 1.0839061 1 0.000000 0 0.000000 0 1 0 0 0.1120 H 1.0944609 1 111.960983 1 0.000000 0 1 2 0 0.1417 C 1.4074935 1 121.910338 1 150.480663 1 1 2 3 0.2970 H 1.0899992 1 130.764118 1 -28.282621 1 4 1 2 0.1019 O 1.2856876 1 108.830602 1 148.022229 1 4 1 2 -0.3616 H 1.3561658 1 82.756661 1 -3.529235 1 6 4 1 0.3019結果が出ましたら、「FOR006」を開いて、チェックします。 以下のようなところで、負号がついた振動数（ FREQ.）を探してください。
VIBRATION 1 ATOM PAIR ENERGY CONTRIBUTION RADIAL
FREQ. -2303.12 O 6 -- H 7 -52.1% (-86.5%) 88.4%
T-DIPOLE 2.8137 C 1 -- H 7 -49.8% 47.8%
TRAVEL 0.1380 C 1 -- C 4 -0.6% 89.0%
RED. MASS 0.7685 C 4 -- H 5 -0.6% 0.2%
VIBRATION 2 ATOM PAIR ENERGY CONTRIBUTION RADIAL
FREQ. 363.43 O 6 -- H 7 45.6% (112.3%) 53.9%
T-DIPOLE 2.3894 C 4 -- O 6 40.6% 3.0%
TRAVEL 0.1165 C 1 -- C 4 37.8% 9.1%
RED. MASS 6.8287 C 1 -- O 6 31.8% 93.1%
C 4 -- H 5 19.3% 100.0%
C 1 -- H 3 11.2% 0.3%
.
.
.
最初の一個以外に、負の値を持っている振動がありましたでしょうか？
もし、たった一つの負の値（正確に言うと「虚の振動数」という状態を示している）を有しているのであれば、先ほど求めた構造が遷移状態であったことが証明できたことになります。
長い道のりでした（慣れれば楽ですが）
さて、MolCatの出番です。
と、その前に、前処理ツールである“VibAnal”を使います。
振動解析の結果で得られた「FOR006」ファイルを、“VibAnal”にドラッグ＆ドロップしてください（MolCatシリーズはアイコンがかわいいです）。
GAUSSIANか、MOPACか、形式を聞いて来ますので、「MOPAC」を選んでください。
しばらく処理した後、どの振動モードを可視化するか聞いてきます。とりあえず、全部選んでおけば間違いありません。
得られたファイル群のうち、表示したい振動モードの番号（振動解析の結果から調べてください）のファイルを“MolCat”にドラッグ＆ドロップします。
すると、振動モードが表示されます。右の座標やダイヤルで、分子を回転して視点をあわせることが出来ます。
右の図では、酸素原子から離れた水素原子が、隣の炭素原子に移って行こうとする様子がわかります。同じく、結合相手の炭素の側も、近づいています。図示化するとわかりやすいですね。
また、メニューから分子の表示サイズや原子の半径、ベクトルの長さ（相対比は保存されます）、ベクトルの反転などが行なえます。見やすいように指定してください。
ここで“PSI-88”が登場します。
まず、計算したいファイルを、キーワード“GRAPH”を使って計算し、電子雲の座標を書き出します。
そうして得られた出力ファイルの中の「FOR013」に目的の情報が書き込まれています。
これを、PSI-88と同じフォルダに移動し、右のような順で処理ます。ドラッグ＆ドロップには対応していないので、preplotをダブルクリックして起動します。
その後、psi1→psi2→psiconの順で起動していくと、最終的に「PLOT」というファイルが出来ます（ダイアログがいくつか出ますが、マニュアルを参照してください。HOMOとLUMOは何番です、と表示されるので、それを指定してください。あとはリターンのみでも行けます）。
これはヒューレット・パッカードのワークステーション形式の画像ファイルなので、「GraphicConverter」を使って開きます。ドラッグ＆ドロップしてください。
すると、図のように電子が表示されます。
残念ながら白黒画像で、ビットマップなので回転もできません。まぁ、表示が出来るだけ良しとしましょう。
MOPACまでの精度が必要ない場合は、「HMO-Plus」を使うともっと簡単に表示出来ます。
MOPAC93のそれに比べるとまだいろいろ問題が多いのですが、傾向を知るような目的には役立つことと思います。計算方法ですが、MOPAC7にて計算を行なうときに、キーワード“EPS=XXX”を使用するだけです。XXXXには数値が入ります。具体的には、目的とする溶媒の誘電率です。水ならば78.3です。
溶媒と電荷のやり取りがあるような場合には適用し難いなど、COSMO法自体にもいくつか弱点があります。参考書などでよく理解した上で使ってください。