谈谈该从Gaussian输出文件中的什么地方读电子能量

sobereva

谈谈该从Gaussian输出文件中的什么地方读电子能量

Where should I read electronic energy from Gaussian output file?

文/Sobereva@北京科音

First release: 2019-Jun-1    Last update: 2023-Sep-16

1 前言


在网上答疑时，经常看到有Gaussian初学者问诸如“我要算xxx，是不是用HF能量”、“我应该用SCF Done能量对比么？”。这种问法被内行人看了很烦，完全是根本没入门、没常识的问法，回答他本来想问的问题前还得先给对方科普一下从输出文件里该读什么数据。鉴于每次都科普一遍太麻烦、甚至都懒得回复，索性这里专门写个小文章说一下对于常见类型的计算，该在什么地方读电子能量。

电子能量的计算属于最最基础的知识，若连这都还不懂的话非常建议参加北京科音初级量子化学培训班（http://www.keinsci.com/workshop/KEQC_content.html）系统、完整学学量子化学计算的各种基本知识。


2 什么叫电子能量、怎么得到

在介绍哪里读电子能量之前，笔者觉得有必要先给初学者简单科普一下什么叫电子能量。我发现很多人甚至连电子能量和自由能都分不清楚。

首先要知道量子化学程序能算出的描述体系的能量的量包括电子能量（electronic energy）以及一些热力学数据，即内能、焓、自由能。

一般量子化学书籍、文献、程序里所谓的“电子能量”包括四项：(1)电子的动能(2)电子与电子间的库仑互斥能(3)核与核之间的库仑互斥能(4)电子与核之间的库仑吸引能。所以一般所说的“电子能量”其实并不完全是它的字面意思，因为也把核之间的互斥能包含了进去。

电子能量的能量零点是假设核与电子都没有动能，体系中所有电子和原子核都被分离到无穷远的情况的能量。这个能量也可以被视为是体系的绝对能量。其数值本身并没有什么化学意义，只有通过求差来得到与物理/化学上研究的问题对应的量的时候才能体现出意义，如反应能、电离能。除了对微型体系使用极高精度的方法，否则电子能量的计算是不可能达到定量准确的，但由于求差的时候大部分误差都可以被抵消，因此目前常用的方法给出的有化学意义的数据的精度还是不错的。

任何量子化学计算任务都可以给出电子能量。量子化学任务中最简单的就是单点（Single point）计算任务，指的就是对于一个给定的几何结构去计算电子能量，这也叫单点能。所以“电子能量”和“单点能”这两个词有时候不加区分地使用。如果没有额外任务设定，量子化学程序默认就是算单点任务。

单点计算以外的任务也可以给出电子能量，只不过它们的目的不是专门算电子能量，电子能量只是作为一个副产物顺带产生的。比如做几何结构优化，最后一次输出的电子能量对应的就是优化后的几何结构下的电子能量。如果你用优化后的结构手动再做一次单点能计算任务，若此时用的计算级别以及各种影响能量的设定都与做几何优化时相同，那么此时得到的电子能量也和之前做优化时最后一次输出的相同（但也有不同的可能，因为Gaussian的几何优化每一步用的波函数初猜都是上一步收敛的波函数。直接算单点任务时，自动产生的初猜不一定总能收敛到与优化任务最后一步相同的波函数）。

做振动分析、NMR计算、极化率计算等等考察基态特征的任务，顺带给出的电子能量和计算单点任务给出的是完全相同的（前提是计算条件完全相同）。

如果想计算内能、焓、自由能这些热力学上的能量，必须要做振动分析才能得到（或者做热力学组合方法计算，如CBS-QB3）。它们与电子能量的一个关键差异就在于这些热力学量里也体现了核运动的贡献，而我们算电子能量的时候是完全忽略掉核的运动的。一般说的零点能（ZPE）对应的是0K下体系的内能/焓/自由能（此时这三者数值相同）与电子能量的差值，来自于原子核的振动运动。显然ZPE也必须通过做振动分析才能得到，电子能量里是不包含ZPE的。

关于上述热力学量的具体计算方法，看《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552）里面的例子。在文中介绍的Shermo程序手册的附录部分有详细的热力学数据计算的基础知识介绍，仔细看了就能彻底理清关系了。


3 单点任务该在哪里读电子能量？

老有初学者搞不清楚单点任务的输出文件中哪项才是自己需要的当前级别下的电子能量，这里专门说一下，这和你用的理论方法有直接关系。Gaussian输出文件里给出电子能量的时候单位都是Hartree，即原子单位（a.u.）制下的能量单位。不同版本的输出信息可能有异，下面的说法至少对于Gaussian 16 A.03是适用的。

3.1 Hartree-Fock

Hartree-Fock（HF）方法需要做SCF（自洽场）迭代，迭代过程中电子能量会一点一点逐渐逼近HF方法的最终结果，迭代收敛后会显示SCF Done，这后面的值显然就是这个结构下HF级别下给出的电子能量了。

另外，在输出文件末尾会有一大堆密密麻麻的输出，这叫archive段落，是当前计算过程中重要信息的汇总。其中HF=后面的数值和之前输出的SCF Done后面显示的完全相同，因此你直接从输出文件末尾来读电子能量也可以。

PS：HF已经彻底过时了，这年头还用HF绝对文章发不出去，绝对没有任何理由使用HF方法。这点我在《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）里专门强调了。


3.2 普通泛函的DFT计算

对于双杂化泛函以外的DFT泛函的计算，如B3LYP、M06-2X、PBE0等，读电子能量也是读SCF Done。

注意，在末尾的archive段落中输出电子能量时此时还是用HF=这个标签来输出，即此时HF=后面的值其实是DFT的结果，和SCF Done后面的值是一样的，而非指的是Hartree-Fock的结果。这是为什么我特别反感有人用诸如“是不是用HF能量”这种问法，谁知道他是真的做了巨垃圾的Hartree-Fock计算，还是指他做的是DFT计算但是指的实际上是HF=后面的那个DFT结果？


3.3 MP2

MP2是后HF方法，计算分为两个过程，先做常规的HF计算，然后基于HF波函数计算MP2能量。因此，SCF Done，以及archive段落的HF=后面的值，指的都是HF级别的电子能量。MP2级别的电子能量是输出文件里EUMP2 =后面的值，这个值和archive段落的MP2=后面的值是相同的。


3.4 双杂化泛函

双杂化泛函里面包含了基于DFT轨道按照MP2公式计算的一部分能量，因此双杂化泛函也分为两步，第一步和普通泛函DFT计算一样要进行SCF迭代，然后输出SCF Done，这个值和archive段落的HF=后面的值是相同的。然后程序会基于收敛的DFT轨道，再做类似MP2的那一部分。比如用的是B2PLYP双杂化泛函，则最终双杂化泛函的能量就是输出文件中E(B2PLYP) =后面的值，这个值在archive段落里也有，但是标签用的是MP2=。

因此，对于双杂化泛函计算，最省事的读取电子能量的做法就是直接从archive段落里读MP2=后头的值。有的初学者误将SCF Done后面的值当做双杂化泛函能量，此时相当于花了双杂化泛函的耗时，然而得到的却是一个垃圾数据。

另外，老有初学者非要用GaussView的Results - Summary界面来读取结果，我强烈不建议大家这样做，因为GaussView从Gaussian输出文件中解析出来的数据可能是有极度误导性的！！！比如用B2PLYP泛函计算，Summary界面里就给出一个能量数据E(RB2PLYP)，这个数据是前面提到的SCF Done后面的数据，根本不是B2PLYP级别的真正的电子能量！这个问题起码对于GaussView 5.0.9和6.0.16都有。显然，初学者盲目相信Summary这个很小儿科界面里显示的数据很容易被坑得头破血流，而只有从输出文件里自己读恰当的数据，才能万无一失。


3.5 CCSD(T)

CCSD(T)是如今最常用的高精度后HF方法。和MP2一样也是先做HF的SCF迭代，然后基于HF轨道再计算CCSD(T)的电子相关部分，最后给出CCSD(T)电子能量。CCSD(T)级别的电子能量是输出文件中CCSD(T)=后面的值，直接从archive段落里读取CCSD(T)=后面的值也是一样的。

CCSD(T)是高级别后HF方法，计算过程中也会利用中间数据顺带着把一些更低级别的后HF能量给算出来。因此在archive段落中你也会看到MP2能量、不同形式的MP4能量、CCSD能量。


3.6 CASSCF方法

CASSCF是个SCF迭代过程，用Gaussian做CASSCF计算在迭代收敛后，会输出各个被计算的电子态的能量和组态系数。可以直接搜EIGENVALUES AND  EIGENVECTORS OF CI MATRIX，下面的EIGENVALUE后面的值就是相应的电子态的电子能量了。archive段落的HF=后面的值也是CASSCF的电子能量，如果同时算了多个态的话，这个值是序号最高的那个态的值。


3.7 TDDFT计算

TDDFT的计算在《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）有非常详细的介绍，其中也专门说了root的设置，它默认为1。输出文件中，root设成了第几激发态，在输出激发态信息的段落中，第几激发态下面就会有一行Total Energy, E(TD-HF/TD-DFT) =，后面的值就是这个激发态的电子能量，相当于基态电子能量加上这个态的激发能。在archive段落里HF=后面的值是基态的电子能量，当前激发态的电子能量在archive段落里并没有输出。


3.8 半经验方法

半经验方法如AM1、PM7等，计算过程和HF一样要做SCF迭代，之后输出的SCF Done的值就是半经验方法的能量了，在archive段落里依然是以HF=作为标签来输出。

特别要注意的是，半经验方法在原作者拟合参数的时候大多都是朝着生成焓来拟合的，因此Gaussian做半经验方法计算给出的能量并不是一般意义的电子能量，而是当前体系在标况下的生成焓。由于生成焓是个相对能量，这也是为什么半经验方法计算得到的能量比起用前面那些方法的数量级都要小得多得多。


3.9 分子力学计算

Gaussian虽然是量子化学程序，但也支持基于AMBER、UFF等分子力场做分子力学计算。从输出文件中搜Energy=，后面的值就是分子力学能量，在archive段落以HF=标签来输出。

注意分子力学计算给出的能量也不是一般意义的电子能量，分子力学计算时根本都没有把原子核与电子单独进行描述，而是每个原子被整个当做一个质点考虑。分子力学的能量零点是每个几何变量都恰处于分子力场定义的能量最低位置。

由于半经验方法和分子力学方法给出的能量都不是体系的绝对能量，虽然你说成“电子能量”别人也知道你指的是什么，但是说成“单点能”或简单说成“能量”从语义上更合适。


3.10 热力学组合方法

G4、G4(MP2)、CBS-QB3等热力学组合方法都是几何优化、振动分析、高精度电子能量计算的过程的组合，目的是得到各种热力学量。如果你只想得到这些方法自动优化出的极小点结构处的电子能量，那就把程序输出的0 K下的内能减去ZPE即可。比如G4(MP2)任务末尾可看到
E(ZPE)=                     0.021064
G4MP2(0 K)=               -76.355852
即自动优化出的极小点结构下的电子能量为-76.355852-(0.021064) = -76.376916 Hartree。

如果你就是希望用这些热力学组合方法计算当前结构下的单点能，而不自动优化结构，也不做对当前没有意义的振动分析，就在热力学组合方法关键词后面写上=SP。比如关键词写G4MP2=SP，就是对当前结构用G4(MP2)方法定义的电子能量计算方式算单点能，此时在输出文件末尾看到的G4MP2 Energy=后面的数值就是了，也可以从archive段落部分读取G4MP2=后面的值，是一样的。

以上情况对于其它热力学组合方法也都是高度类似的，只不过输出的标签里面的名字不一样而已。


4 总结

希望看过本文后，读者能清楚知道当前用的理论方对应的电子能量应该从哪里读。也绝对不要再在提问时不给出任何前提就直接说什么“HF能量”、“SCF Done值”这种令内行人看着头疼、纠结的含糊不清的表达。正确、没有丝毫歧义的表达应当是比如“双杂化泛函xxx下算的电子能量”、“PM6方法下算的单点能”。另外，没有前提时，直接说“能量”的时候，内行人都一律会默认当做你说的是电子能量/单点能，所以如果你要问热力学量相关问题的时候必须明确说清楚是哪种热力学量。

如果你做的是几何优化任务，输出文件末尾archive段落输出的能量数据都是最后一步结构的。

星斗如盘

社长，想请教您一下从输出文件直接读取的电子能量和从Gview读取的能量应该是一致的吧？
————问题已解决，怪我没有仔细看博文，跳过了双杂化泛函的部分。

zjxitcc

星斗如盘 发表于 2020-6-2 11:47
社长，想请教您一下从输出文件直接读取的电子能量和从Gview读取的能量应该是一致的吧？

提问暴露了没仔细读博文。双杂化那一大段就被你忽略了

星斗如盘

zjxitcc 发表于 2020-6-2 13:03
提问暴露了没仔细读博文。双杂化那一大段就被你忽略了

是这样子。。。我只用B3LYP计算过，所以自动把这一部分忽略掉了，原来还真的不一样

a1063021280

受用 感谢社长

一条君

老师，双杂化泛函加色散校正后B2PLYPD3的单点能就是读取E(B2PLYPD3) =后面的值吧，谢谢~

sobereva

一条君 发表于 2023-8-6 11:44
老师，双杂化泛函加色散校正后B2PLYPD3的单点能就是读取E(B2PLYPD3) =后面的值吧，谢谢~

显然不是

一条君

sobereva 发表于 2023-8-6 17:36
显然不是

不是吗老师？那应该在mp2那里读？我仔细看了一下读E(B2PLYPD3) =处的和mp2数值一样呀

sobereva

一条君 发表于 2023-9-4 21:46
不是吗老师？那应该在mp2那里读？我仔细看了一下读E(B2PLYPD3) =处的和mp2数值一样呀

谈谈该从Gaussian输出文件中的什么地方读电子能量
http://sobereva.com/488（http://bbs.keinsci.com/thread-13450-1-1.html）
上面该说的都说了

如果有怀疑，说清楚你从哪看的，前后输出信息都贴出来

一条君

sobereva 发表于 2023-9-6 23:44
谈谈该从Gaussian输出文件中的什么地方读电子能量
http://sobereva.com/488（http://bbs.keinsci.com/th ...

老师，（1）请问双杂化泛函加色散校正后B2PLYPD3的单点能就是读取E(B2PLYPD3) =后面的值吧
我是在加粗部分读的The integrals were generated   1 times. E2(B2PLYPD3) =    -0.8155994207D-01 E(B2PLYPD3) =    -0.76407640512291D+02 或archive段落里MP2=-76.4076405
（2）普通泛函加dft色散校正后，读取能量的位置仍是scf done=后面吗
谢谢老师

sobereva

一条君 发表于 2023-9-7 10:24
老师，（1）请问双杂化泛函加色散校正后B2PLYPD3的单点能就是读取E(B2PLYPD3) =后面的值吧
我是在加粗部 ...

是
是

sobereva

给此文增加了热力学组合方法输出电子能量的说明

嘤嘤嘤

老师，我用gaussian软件采用ROCBS-QB3=sp计算单点能，那我的能量是读取最后的E(CBS-QB3)=  的能量吗？但是我通过看文章发现，有人是按照这个公式计算的单点能。请问这俩的结果是一样的吗？

sobereva

嘤嘤嘤 发表于 2024-1-17 16:48
老师，我用gaussian软件采用ROCBS-QB3=sp计算单点能，那我的能量是读取最后的E(CBS-QB3)=  的能量吗？但是 ...

是

你截图的公式里都带了ZPE，显然不是电子能量

嘤嘤嘤

sobereva 发表于 2024-1-18 03:21
是

你截图的公式里都带了ZPE，显然不是电子能量

请问老师，截图中的公式是什么意思啊？因为文章中说的很清楚，DFT条件下优化，CBS-QB3计算单点，结果给出的单点能计算公式是那个。难道这说明作者算错了？不过感觉截图中的公式有点合理，毕竟是热力学组合方法，最后单点能是几种方法的加和，但是后边那些加减的量不是很理解，比如the empirical higher-order correction −0.00579E(OIii）