GAMESSでMEP計算。
久々の投稿、そして久々の計算化学の記事。
以前、GAMESSではMOPACのMEP計算のような計算ができないということを書きました。実は、今日CC-BBSで同じことが質問に出ていました。そこで再度挑戦したところ、なんとできてしまったわけです。
ちょっと裏技じみた(と思っているのは私だけ?)やり方ですが、こういうのを発見するのも、開発する立場ではなく使う立場の楽しみで…
基本はやはりSURFACE計算。
そしてRUNTP1=OPTIMIZEも変わりません。でも、直交座標系を使い、$STATPTのIFREEZで対象原子の座標(x,y,z)を固定することにより、擬似MEP計算へ導くことができます。
以下の入力は、メチレンジアミン(アミナール)がプロトン化された構造からアンモニアが脱離していく過程を、C-N(アンモニア)結合に沿った構造・エネルギー変化で追う計算です。
$CONTRL SCFTYP=RHF RUNTYP=SURFACE COORD=UNIQUE
MAXIT=200 NZVAR=0 ICHARG=1 $END
$SYSTEM TIMLIM=600 MWORDS=1 $END
$STATPT NSTEP=500 OPTTOL=0.0001 IFREEZ(1)=1,2,3,4,5,6 $END
$BASIS NGAUSS=3 GBASIS=STO $END
$SCF DIRSCF=.T. FDIFF=.F. DAMP=.T. $END
$GUESS GUESS=HUCKEL $END
$SURF RUNTP1=OPTIMIZE IVEC1(1)=1,2 IGRP1(1)=2,6,7,8 ORIG1=0.0
DISP1=0.2 NDISP1=5 $END
$DATA
test
C1
C 6.0 -0.110256 -0.472297 -0.384997
N 7.0 1.172596 -0.033181 0.386590
N 7.0 -1.140348 0.487881 -0.041832
H 1.0 0.182512 -0.471211 -1.441351
H 1.0 -0.307033 -1.495451 -0.044485
H 1.0 1.441267 0.934725 0.112543
H 1.0 0.984592 -0.013030 1.410551
H 1.0 1.983248 -0.663218 0.211685
H 1.0 -1.936309 0.034042 0.421060
H 1.0 -1.482556 0.983506 -0.873140
$END
太字&下線部が重要なキーワードです。これにより、平衡構造(C-N : 1.56Å)から0.2Åずつ1.66Åまで結合を伸ばしたときの構造とエネルギーの変化を追うことができます。この計算の出力は、Winmostarを使って各結合長での最適化の流れまで含めた形で可視化できます。これを使えるだけで、ずいぶんと反応解析はやりやすくなります。あとは、PC GAMESSで同じことをやる方法があるかどうか…(PC GAMESSでは$SURFにRUNTPnオプションがありません)