フリーウエアで非経験的分子軌道計算を実行（Mac）
自動変換ファイルです．コピーにより移行した修正ファイルはこちらです．
経験的分子軌道計算法の最新バージョンPM6法はStewartによって，研究目的の場合，無料で利用できる．ところが，非経験的分子軌道計算となると高価なGaussian社のソフトを購入し，並列計算機を用いて計算処理するのが普通となっている．私の場合，大学ではPemtium ４ ３GHz PCを16台並列化した計算システムやItanium 2の並列処理計算システム（2ノード×４）を利用してきた．そのためのソフトの購入費用，ハードウエアを合わせると400万円程度は必要である．そのために予算獲得が必至であり，数年毎に悩まされてきた．
ところが，近年マルチコアのCPUが家庭用のパソコンに採用されるようになり，高価な並列計算機は必ずしも必要としないという報告を聞くようになった．また．ソフトウエアもマルチコア用のプログラム (GAMESS) がPC対応のフリー版で公開されていることがわかった．さっそく，Firefly Project TeamのDr. Alex A. Granovskyにメールを送り利用許可を得た．当初Windows環境で実行するつもりでいたが，Intel Mac上でも動くことが分かり，計算処理前後のプログラムを含めIntel Mac (2 core, 2GHz, RAM 4GB, 2009) で実行環境を整備してみた．
以下はそのメモである．
入力データ作成および結果解析プログラム
１）MacMolPlt
２）Avogadro
機能は両者ほとんど同じであるが，一長一短である．
前者は振動ベクトルの表示が可能，分子描画上のバグがある（実行例を参照）．
後者は簡単な構造最適化法が用意されているので適当に作成したラフな構造を整形することができる．
また，特定の結合を任意の距離に変更することも可能である．
ブタジエンとエチレンのDiels-Alder反応のような対称型の遷移状態はAvogadroでも試行錯誤しながら遷移状態に近い入力構造を作成することができるが，今回取り上げる[3,3]-sigmatropy 反応（２→３）のような非対称型の６員環遷移状態構造は簡単には作れない．Avogadroで適当に作成した６印環遷移状態をsaddle計算しても鞍点構造は見つからなかった．もっとも簡単な方法はZ-matrixを入力座標とするMOPAC計算で反応点間の距離（赤矢印で示した硫黄原子とγ位の炭素の距離を固定し，徐々に接近させる方法で生成熱のもっとも高い点の分子座標を入力座標とする方法である．
分子軌道計算プログラム
１）Firefly（実行形式）
PCのコア数によって選択可能
Run-Firefly-Job-1-Core
Run-Firefly-Job-2-Core
Run-Firefly-Job-4-Core
Run-Firefly-Job-8-Core
２）GamessQ（実行形式）
コア数を入力するように促されるが，認識しなかった （mac miniの場合）．
Fireflyで利用可能な基底関数
上で述べたように，ab initio計算のための入力座標を半経験的MO法で作り，その座標をもとに計算を実行した結果を以下に示す．
Fireflyの入力座標
PM6のTS構造をフロントエンドのMacMolPlt（フリーウエア）またはAvogadro（フリーウエア）に読み込んで，SubwindowのInput Builderを使って必要事項を慎重に入力する．
計算目的は，次図に示す連続ペリ環状反応の第一段階であるO-allyl S-methyl dithiocarbonateの[3,3]-sigmatropyの遷移状態を探すことである．
B3LYP/3-21G(d)とB3LYP/6-31G(d)の両 Basis Setで実行した．
左側の各項目について，必要事項を選択後 Write FileをクリックしてGAMESS計算に必要なinputファイルを生成させる．
修正する必要があるときは，テキストエディタで修正・追加の書き込みを行う．
MacMolPlt で生成させた入力データ
------------------------------------------------------------------------------------------------------------------------------
! File created by MacMolPlt 7.4.2
$CONTRL SCFTYP=RHF RUNTYP=SADPOINT DFTTYP=B3LYP MAXIT=30 MULT=1 $END
$SYSTEM TIMLIM=525600 MEMORY=1000000 $END
$BASIS GBASIS=N21 NGAUSS=3 NDFUNC=1 $END
$SCF DIRSCF=.TRUE. $END
$STATPT OPTTOL=0.0001 NSTEP=200 HESS=CALC $END
! $FORCE METHOD=ANALYTIC VIBANL=.TRUE. $END 使用しない
$DATA
Title
C1
C 6.0 0.00000 0.00000 0.00000
S 16.0 1.69510 0.00000 0.00000
C 6.0 2.07260 2.10660 0.00000
C 6.0 0.86940 2.61870 -0.46920
C 6.0 -0.29440 2.51670 0.33550
O 8.0 -0.80070 0.94180 0.27350
S 16.0 -0.74300 -1.54930 -0.39990
C 6.0 -2.51430 -1.23390 -0.33250
H 1.0 2.95420 2.12610 -0.64670
H 1.0 2.32040 2.15820 1.06490
H 1.0 0.76150 2.88280 -1.51540
H 1.0 -1.23270 2.93520 -0.04620
H 1.0 -0.20050 2.58100 1.42600
H 1.0 -3.02460 -2.16130 -0.03980
H 1.0 -2.89870 -0.91290 -1.30880
H 1.0 -2.77120 -0.45620 0.40200
$END
------------------------------------------------------------------------------------------------------------------------------
2行目はB3LYPの鞍点計算である．
3行目は3-21g(d)の指定である．
5行目のHESS=CALCは振動計算を行うことを指定．NSTEP=200は大きめに設定する，変更しなければ20になる．
ジョブの実行（遷移状態探索）
Mac mini (Core2 duo, 2GHz, RAM 4GB, 2009)であるので，Run-Firefly-Job-2-Coresを起動（インプットファイルをアプリケーションにかぶせるだけでよい）．
"SADDLE POINT LOCATED"が表示されジョブが終了する．
計算時間 ２０分
TIMING STATISTICS ON NODE 0:
CPU TIME: STEP = 0.00 , TOTAL = 1147.5 SECONDS ( 19.1 MIN)
WALL CLOCK TIME: STEP = 0.00 , TOTAL = 1208.2 SECONDS ( 20.1 MIN)
CPU UTILIZATION: STEP = 0.00%, TOTAL = 94.98%
AN INPUT FILE FOR -MOLPLT- HAS BEEN PUNCHED.
出力ファイルをMacMolPltで開き，各ステップの構造とEnergy Plotを表示させる．
この場合，初期構造から最終構造まで４０フレームの分子図を閲覧することができる．
遷移状態構造の候補をexportする．
注 途中結果を見て，エネルギーが一定になったらジョブを打ち切ってもよい．
SADDLE計算処理前後の結合距離の変化（黒→赤）
SADDLE探索失敗の場合
フロントエンドのビルダーを利用して適当に作った遷移状態構造？を使って計算させると，鞍点を見い出せず下記のメッセージを表示して終了した．
***** FAILURE TO LOCATE STATIONARY POINT, TOO MANY STEPS TAKEN *****
UPDATED HESSIAN, GEOMETRY, AND VECTORS WILL BE PUNCHED FOR RESTART
最後のフレームを入力座標としてPM3でSADDLE計算すると ***** SADDLE POINT LOCATED *****を表示して正常終了した．PM3の演算時間は2.3秒．半経験的分子軌道計算は非常に速いので，非経験的分子軌道計算を実行する前により合理的な初期座標を得るために試みた方がよい．
振動計算
入力データを作成するルーチンではHESSIANを選択する．
----------
! File created by MacMolPlt 7.4.2
$CONTRL SCFTYP=RHF RUNTYP=HESSIAN DFTTYP=B3LYP MAXIT=30 MULT=1 $END
$SYSTEM TIMLIM=525600 MEMORY=1000000 $END
$BASIS GBASIS=N21 NGAUSS=3 NDFUNC=1 $END
$SCF DIRSCF=.TRUE. $END
! $FORCE METHOD=ANALYTIC VIBANL=.TRUE. $END
$DATA
Title
C1
C 6.0 0.03766 -0.15939 0.32907
S 16.0 1.72757 -0.21898 0.51922
C 6.0 1.99340 2.18506 -0.15277
C 6.0 0.70743 2.54026 -0.51341
C 6.0 -0.28513 2.56567 0.47762
O 8.0 -0.75259 0.75874 0.78405
S 16.0 -0.65577 -1.47246 -0.67472
C 6.0 -2.44714 -1.09171 -0.56049
H 1.0 2.75950 2.02179 -0.90248
H 1.0 2.34450 2.36483 0.85619
H 1.0 0.40689 2.47824 -1.55404
H 1.0 -1.29985 2.84976 0.22659
H 1.0 0.00674 2.75021 1.50572
H 1.0 -2.96926 -1.89664 -0.03995
H 1.0 -2.85962 -0.96157 -1.56293
H 1.0 -2.52123 -0.16031 0.00574
$END
----------
MacMolPltによる解析
Avogadro (Version 1.0.3) で表示させた場合
振動計算で負のベクトルが一個だけ存在することが判る．
どのような振動かベクトルが図示されるので非常に分かりやすい．
この後，irc計算を行い，鞍点構造が原料および生成物へ繋がっていることを確認する．
計算は2CPUが100％稼働状態で実行された．
B3LYP/6-31G(d)レベルの計算は60.6minを要した（５２フレーム）．
一頃の大型計算機センターに匹敵する演算速度である．
参考資料
半経験的分子軌道法PM3を実行する場合の入力座標（演算時間 2.0秒, 1 coreにのみ対応）
! File created by MacMolPlt 7.4.2
$CONTRL SCFTYP=RHF RUNTYP=SADPOINT MAXIT=30 MULT=1 $END
$SYSTEM TIMLIM=525600 MEMORY=1000000 $END
$BASIS GBASIS=PM3 $END
$SCF DIRSCF=.TRUE. $END
$STATPT OPTTOL=0.0001 NSTEP=1200 HESS=CALC $END
! $FORCE METHOD=SEMINUM VIBSIZ=0.010000 VIBANL=.TRUE. $END
$DATA
pm
C1
C 6.0 0.00000 0.00000 0.00000
S 16.0 1.69510 0.00000 0.00000
C 6.0 2.07260 2.10660 0.00000
C 6.0 0.86940 2.61870 -0.46920
C 6.0 -0.29440 2.51670 0.33550
O 8.0 -0.80070 0.94180 0.27350
S 16.0 -0.74300 -1.54930 -0.39990
C 6.0 -2.51430 -1.23390 -0.33250
H 1.0 2.95420 2.12610 -0.64670
H 1.0 2.32040 2.15820 1.06490
H 1.0 0.76150 2.88280 -1.51540
H 1.0 -1.23270 2.93520 -0.04620
H 1.0 -0.20050 2.58100 1.42600
H 1.0 -3.02460 -2.16130 -0.03980
H 1.0 -2.89870 -0.91290 -1.30880
H 1.0 -2.77120 -0.45620 0.40200
$END
計算データの日付 2011年10月 (aoxsaddle, aox_hess etc)