DC法をブラックボックス的に使う落とし穴?

[本記事は私が以前書いていたblogからサルベージしたものです]

ブラックボックス的に使えれば何よりですが、そうもいかないようで。

一例として、抗うつ薬フルオキセチン(Fluoxetine, SSRI)の2つの配座間のエネルギー差を調べる計算を示しますと…

GAMESSに実装されたDC法では、プログラムにSubsystemとbuffer領域を構築してもらうのが最も手っ取り早く計算を行う方法です。例えばこんな感じで入力を組みます。

$CONTRL SCFTYP=RHF RUNTYP=ENERGY COORD=UNIQUE　$END
$SYSTEM TIMLIM=36000 MWORDS=10 $END
$BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 $END
$GUESS GUESS=HUCSUB $END
$DANDC DCFLG=.T. SUBTYP=AUTBND
BUFTYP=RADSUB BUFRAD=8.0 $END
$DATA
fluoxetine, conf-1 : DC-RHF/6-31G(d) energy
C1
C 6.0 5.707825 2.291899 -3.172497
H 1.0 6.812142 2.291895 -3.172497
…
$END

どのキーワードがどういう意味かはマニュアルや中井先生の論文(邦文, JCCJ)などをご覧頂きたいと思いますが、ともかく、SUBTYP=AUTBNDを指定すると、分子を適当な格子(この目の細かさも別途指定できる)で区切った上で、結合を考慮して(結合長から判断し、二重結合をぶった切ったりしない)Subsystemを設定してくれるわけです。
で、このありがたいキーワードを使ってエネルギー計算を行った結果が以下の表です(6-31G(d)基底)。

	RHF, E(Conf-1)-E(Conf-2)	MP2, E(Conf-1)-E(Conf-2)
通常の計算	1.122	1.110*
分割統治法	0.954	-0.606

(数値はkcal/mol, *Firefly 7.1.Eの計算値) ※DC-MP2のHF計算は非分割で実施(DCFLG=.F.)

RHFでは両者のずれは大きくありませんが、MP2では符号すら違う結果になります。これはなぜか？？
その答えはおそらくSubsystemの切り方の違いにあります。

2つの配座では、同じ目の細かさの空間格子で分子を区切っても、Subsystemの構成が異なる部分があるのです。察しの良い方はすぐに気づくと思いますが、メチルアミノエチル側鎖の部分に違いが出ます。

Conf-2の分割法をConf-1に適用し(Conf-2のPunchファイルから $SUBSCF をコピーし、SUBTYP=CARD, BUFTYP=CARDとする)、再度計算を行うと、このようになります。

	RHF, E(Conf-1)-E(Conf-2)	MP2, E(Conf-1)-E(Conf-2)
通常の計算	1.122	1.110*
分割統治法	0.954	-0.606
分割統治法 (修正)	1.118	1.759

RHFはかなりよい一致になりますし、MP2も符号が合います(絶対値はまだ少し差がありますが)。
ブラックボックス的にプログラムにおまかせすると、思わぬ落とし穴に嵌ることもあり得ますので、計算結果は慎重に見つめなければなりませんね。

おまけ

フルオキセチンのDC-CCSDも実施してみたのだが、CCSD計算に入ったところで計算が落ちる…何のエラーメッセージも出ずに…指定メモリ量も足りてるし…なぜだ。

09.02.08 追記
DC-CC計算はConventional SCFで実行すればOK。Direct SCFをいつもの癖で指定していたせいだった。気づくまで時間がかかったなぁ…(マニュアルにちゃんと書いてあるのに)