量子化学計算により出力されるデータの中で、必ずと言っていいほど利用されるのがエネルギー値です。その値と出力パターンは計算の方法(プログラム・計算目的・モデル化学)の違いで変わってきます。
ここでは、GAMESS(WinGAMESS/PC GAMESS)の出力ファイルのどこにエネルギーが記載され、その値はどのように処理すべきかについて、ごく簡単に述べたいと思います。
Single Point / Geometry Optimization
エネルギー計算(Single Point, RUNTYP=ENERGY)も構造最適化計算(Geometry Optimization, RUNTYP=OPTIMIZE)も、エネルギー値の出力に関しては共通しており、いずれも全電子エネルギーEel(全エネルギーから熱力学的な因子を除いた部分)が算出されます。構造最適化では、入力構造と最適化構造の2つについて詳細な値が出力され、その点のみエネルギー計算と異なります。以下、具体的に出力ファイルの一部を紹介します。
Eelの出力例 (RHF)
------------------------------
properties for the RHF density
------------------------------
-----------------
ENERGY COMPONENTS
-----------------
WAVEFUNCTION NORMALIZATION = 1.0000000000
ONE ELECTRON ENERGY = -123.2240816654
TWO ELECTRON ENERGY = 37.9247670254
NUCLEAR REPULSION ENERGY = 9.2885681244
------------------
TOTAL ENERGY = -76.0107465155
ELECTRON-ELECTRON POTENTIAL ENERGY = 37.9247670254
NUCLEUS-ELECTRON POTENTIAL ENERGY = -199.0837581392
NUCLEUS-NUCLEUS POTENTIAL ENERGY = 9.2885681244
------------------
TOTAL POTENTIAL ENERGY = -151.8704229894
TOTAL KINETIC ENERGY = 75.8596764739
VIRIAL RATIO (V/T) = 2.0019914406
これがRHF計算でのエネルギーの出力です。UHFやROHFでも変わりません。下線部の「TOTAL ENERGY」が実際に比較等に用いるEelに当たります。DFT計算の場合も、タイトルの所が「properties for the R-B3LYP density」のように汎関数の名称に置き換わるだけで、出力フォーマットは全く変わりません。
MP2の場合もフォーマットは変わりませんが、同じようなエネルギー出力が2回出ます。即ち、上記のような「RHF density」の結果が印字された後、それにMP2補正が加わった結果が印字されます。当然、「properties for the MP2 density」の方のTOTAL ENERGYが見るべき値です。
Eelの出力例 (RMP2)
-------------------------------------------------
SCF PROPERTIES...FOR THE UNPERTURBED WAVEFUNCTION
-------------------------------------------------
------------------------------
properties for the RHF density
------------------------------
-----------------
ENERGY COMPONENTS
-----------------
WAVEFUNCTION NORMALIZATION = 1.0000000000
ONE ELECTRON ENERGY = -122.8608464837
TWO ELECTRON ENERGY = 37.7640363230
NUCLEAR REPULSION ENERGY = 9.0870142835
------------------
TOTAL ENERGY = -76.0097958772
ELECTRON-ELECTRON POTENTIAL ENERGY = 37.7640363230
NUCLEUS-ELECTRON POTENTIAL ENERGY = -198.6385330316
NUCLEUS-NUCLEUS POTENTIAL ENERGY = 9.0870142835
------------------
TOTAL POTENTIAL ENERGY = -151.7874824250
TOTAL KINETIC ENERGY = 75.7776865478
VIRIAL RATIO (V/T) = 2.0030630300
.....................(中略)..................................
-----------------------------------------------
MP2 PROPERTIES...FOR THE 1ST ORDER WAVEFUNCTION
-----------------------------------------------
MP2 NATURAL ORBITAL OCCUPATION NUMBERS ARE
2.0000 1.9876 1.9739 1.9706 1.9697 0.0238 0.0218 0.0175 0.0102 0.0052
0.0048 0.0044 0.0042 0.0041 0.0008 0.0007 0.0004 0.0002 0.0001
THERE ARE 9.9018 ELECTRONS IN OCCUPIED RHF ORBITALS.
THERE ARE 0.0982 ELECTRONS IN VIRTUAL RHF ORBITALS.
MP2 NATURAL ORBITALS HAVE BEEN PUNCHED.
------------------------------
properties for the MP2 density
------------------------------
-----------------
ENERGY COMPONENTS
-----------------
WAVEFUNCTION NORMALIZATION = 1.0000000000
ONE ELECTRON ENERGY = -122.6679970253
TWO ELECTRON ENERGY = 37.3841349458
NUCLEAR REPULSION ENERGY = 9.0870142835
------------------
TOTAL ENERGY = -76.1968477960
ELECTRON-ELECTRON POTENTIAL ENERGY = 37.3841349458
NUCLEUS-ELECTRON POTENTIAL ENERGY = -198.5993485420
NUCLEUS-NUCLEUS POTENTIAL ENERGY = 9.0870142835
------------------
TOTAL POTENTIAL ENERGY = -152.1281993127
TOTAL KINETIC ENERGY = 75.9313515167
VIRIAL RATIO (V/T) = 2.0034965304
MP4(PC GAMESSのみ)の出力はMP2とは異なります。以下のような出力です。この例ではMP4(SDQ)計算なので、一番下の「E(MP4-SDQ)」の値を参照します。
Eelの出力例 (RMP4(SDQ))
RESULTS OF MOLLER-PLESSET 4TH ORDER CORRECTION ARE
E(RHF) = -76.0481301848
E(D,2) = -0.2411794348
E(MP2) = -76.2893096197
E(D,3) = -0.0055041932
E(D,2+3) = -0.2466836281
E(MP3) = -76.2948138129
E(S,4) = -0.0011282054
E(D,4) = -0.0039993863
E(D,2+3+4) = -0.2506830144
E(R+Q,4) = 0.0026702882
E(SDQ,4) = -0.0024573035
E(SDQ,2+3+4) = -0.2491409316
E(MP4-SDQ) = -76.2972711164
Coupled-Cluster(WinGAMESSのみ)の出力は以下のようになります。この例はCCSD(T)のもので、最後の行が見るべき値です。
Eelの出力例 (CCSD(T))
SUMMARY OF RESULTS
REFERENCE ENERGY: -76.0481301849
MBPT(2) ENERGY: -76.2893095943 CORR.E= -0.2411794093
CCSD ENERGY: -76.2971790713 CORR.E= -0.2490488864
CCSD[T] ENERGY: -76.3029917488 CORR.E= -0.2548615639
CCSD(T) ENERGY: -76.3028580851 CORR.E= -0.2547279002
THE FOLLOWING METHOD AND ENERGY WILL BE CONSIDERED THE HIGHEST LEVEL RESULT:
COUPLED-CLUSTER ENERGY E( CCSD(T)) = -76.3028580851
Hessian
振動計算(Frequency/Force, RUNTYP=HESSIAN)の出力では、全電子エネルギー以外に熱力学エネルギーが出力されます。全電子エネルギーについては前掲のエネルギー計算と同じですので、熱力学エネルギーの出力について紹介します。
熱力学エネルギーの出力例 (RHF/6-31G(d))
-------------------------------
THERMOCHEMISTRY AT T= 298.15 K※1
-------------------------------
USING IDEAL GAS, RIGID ROTOR, HARMONIC NORMAL MODE APPROXIMATIONS.
P= 1.01325E+05 PASCAL.※2
ALL FREQUENCIES ARE SCALED BY 0.91350※3
THE MOMENTS OF INERTIA ARE (IN AMU*BOHR**2)
2.10181 4.09291 6.19472
THE ROTATIONAL SYMMETRY NUMBER IS 2.0
THE ROTATIONAL CONSTANTS ARE (IN GHZ)
857.87288 440.53944 291.06843
THE HARMONIC ZERO POINT ENERGY IS (SCALED BY 0.914)
0.020989 HARTREE/MOLECULE 4606.656594 CM**-1/MOLECULE
13.171105 KCAL/MOL 55.107905 KJ/MOL※4
Q LN Q
ELEC. 1.00000E+00 0.000000
TRANS. 3.00431E+06 14.915558
ROT. 4.13176E+01 3.721289
VIB. 1.00032E+00 0.000319
TOT. 1.24170E+08 18.637166
E H G CV CP S
KJ/MOL KJ/MOL KJ/MOL J/MOL-K J/MOL-K J/MOL-K
ELEC. 0.000 0.000 0.000 0.000 0.000 0.000
TRANS. 3.718 6.197 -36.975 12.472 20.786 144.800
ROT. 3.718 3.718 -9.225 12.472 12.472 43.412
VIB. 55.114 55.114 55.107 0.172 0.172 0.024
TOTAL 62.551 65.030 8.907 25.115 33.429 188.236※5
E H G CV CP S
KCAL/MOL KCAL/MOL KCAL/MOL CAL/MOL-K CAL/MOL-K CAL/MOL-K
ELEC. 0.000 0.000 0.000 0.000 0.000 0.000
TRANS. 0.889 1.481 -8.837 2.981 4.968 34.608
ROT. 0.889 0.889 -2.205 2.981 2.981 10.376
VIB. 13.173 13.173 13.171 0.041 0.041 0.006
TOTAL 14.950 15.543 2.129 6.003 7.990 44.989※6
......END OF NORMAL COORDINATE ANALYSIS......
振動計算の出力はどのモデル化学でも同じです。まず最初に設定温度が示され(※1，デフォルトは25℃=298.15K)、次に気圧(※2，1 atmに相当)が示されます。温度は$FORCEグループのTEMP(1)=xxx,xxx,...,xxxとカンマ区切りで最大10個の値を指定できますが、気圧は1 atmに固定されています。
スケール因子が示されていますが(※3)、これは$FORCEグループのSCLFACで指定できます。この例はRHF/6-31G(d)なので、そのZPE用スケール因子を用いています。ちなみに、ここで指定されたスケール因子は振動数(IRスペクトルデータ)には適用されません。
中段やや上に零点エネルギー(ZPE)が出力されます(※4)。4種類の単位で出力され、値はいずれもスケールされたものです。下段の2つのテーブルはどちらも同じ内容で、上がジュール(※5)，下がカロリー(※6)で換算された各種熱力学補正値です(Eは内部エネルギー補正(Ecorr),Hはエンタルピー補正(Hcorr),Gはギブズ(自由)エネルギー補正(Gcorr),CVは定積モル比熱,CPは定圧モル比熱,Sはエントロピー)。実際によく使うのは、ZPEと熱力学補正値です。
それぞれの補正値は、以下のように使います。まず、ZPEは0 Kにおける内部エネルギー補正で、並進および回転の寄与が0 Kでは0なので、0 KではZPE=Ecorr=Hcorr=Gcorrです。即ち、0 Kにおける全エネルギーは、
E0 = Eel + ZPE
となります。一方、温度がT Kの時内部エネルギー(E)，エンタルピー(H)，ギブズエネルギー(G)は以下のようになります。
ET = Eel + Ecorr
HT = Eel + Hcorr
GT = Eel + Gcorr
よく、CCSD(T)/6-311+G(2df,2p)//RHF/6-31G(d)といった記述がありますが、これはRHF/6-31G(d)で構造最適化と振動計算(つまり熱力学補正値の計算)を行い、最適化構造についてCCSD(T)/6-311+G(2df,2p)で一点計算(Eelの計算)を行い、両者の計算値を統合して結果を求めたという意味です。熱力学補正と全電子エネルギーが別のモデル化学による計算、というのは一見違和感を覚えますが、振動計算はスケーリングによってモデル化学の差がある程度均されていますし、多くの場合全電子エネルギー計算の方がモデル化学による差が大きいため、全電子エネルギーのみを高度なモデル化学で求め直すことが一般的に認められているようです。
以上、一通り出力されたエネルギー値の見方について簡単に書きました。他にも様々な出力パターンがありますが、それらについてはここで網羅するのではなく、そういった計算のトピックスの中で説明できればと思います。