谈谈不同量子化学程序计算结果的差异问题
谈谈不同量子化学程序计算结果的差异问题
On the differences in calculation results of different quantum chemistry programs
文/Sobereva@北京科音
First release: 2020-Nov-1 Last update: 2023-Jun-18
经常有人问，为什么他用的两个不同量子化学程序在相同级别计算下计算结果存在差异，也常有人问用不同的量子化学程序算的结果能否放在一起对比。在此文就专门说一下导致不同程序间，也包括相同程序的不同大版本间计算结果存在差异可能的因素。首先说最简单的情况，即真空下单点能的计算，之后再说其它任务、涉及额外因素的情况。由于Gaussian和ORCA是最流行的两个量子化学程序，所以对此二者说得相对多一些。
1 导致真空下计算的单点能存在差异的因素
1.1 基组的差异性
角动量>=D的Gauss型基函数分为球谐型和笛卡尔型两种，介绍见《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》（http://sobereva.com/51）和《球谐型与笛卡尔型Gaussian函数的转换关系》（http://sobereva.com/97）。对同一种基组，不同程序可能默认用不同形式的Gauss基函数。比如对于3-21G、6-31G系列等少数基组，Gaussian默认用笛卡尔型D基函数，而对>=F角动量的默认用球谐型；对于6-311G系列、def2系列、cc-pVnZ等系列，以及用gen的时候，Gaussian默认用球谐型。而ORCA程序只支持球谐型基函数，因此对6-31G系列基组的计算结果必然与Gaussian默认情况下的不同。而在NWChem中，不管是什么基组，都默认用笛卡尔型。
有些程序可能会对基组定义做轻微调整，Gaussian就是典型。在Chem. Phys. Lett., 260, 514 (1996)中Davidson指出，对于Dunning相关一致性系列这种（部分）广义收缩的基组，对于Gaussian等绝大多数没有为广义收缩优化积分代码的程序，可以利用类似高斯消元法的做法去除一些指数重复的GTF壳层来节约时间，而精度不会受到任何影响，下文称作“第一级简化”。如果想更进一步节约时间，还可以把收缩系数很小的GTF壳层给去掉，但会轻微影响计算结果，下文称作“第二级简化”。Gaussian默认情况下就会做这两级简化。例如碳的cc-pVTZ基组有两个广义收缩的S型基函数，Gaussian格式的原始定义为（第一列是指数，第二列是收缩系数）
S 8 1.00
8236.0000000 0.0005310
1235.0000000 0.0041080
280.8000000 0.0210870
79.2700000 0.0818530
25.5900000 0.2348170
8.9970000 0.4344010
3.3190000 0.3461290
0.3643000 -0.0089830
S 8 1.00
8236.0000000 -0.0001130
1235.0000000 -0.0008780
280.8000000 -0.0045400
79.2700000 -0.0181330
25.5900000 -0.0557600
8.9970000 -0.1268950
3.3190000 -0.1703520
0.3643000 0.5986840
Gaussian计算时如果用gfinput关键词输出基组定义，会发现变成了下面这样，可见这两个基函数的收缩度分别减少了1和2，在去除了部分指数相同的GTF的同时，收缩系数也发生了变化。
S 7 1.00 0.000000000000
0.8236000000D+04 0.5419783203D-03
0.1235000000D+04 0.4192873817D-02
0.2808000000D+03 0.2152216205D-01
0.7927000000D+02 0.8353432195D-01
0.2559000000D+02 0.2395828457D+00
0.8997000000D+01 0.4428528419D+00
0.3319000000D+01 0.3517995618D+00
S 6 1.00 0.000000000000
0.2808000000D+03 -0.5949224937D-04
0.7927000000D+02 -0.1148158310D-02
0.2559000000D+02 -0.1001913745D-01
0.8997000000D+01 -0.6121949230D-01
0.3319000000D+01 -0.1732698541D+00
0.3643000000D+00 0.1072915192D+01
当Gaussian计算时看到比如Ernie: 6 primitive shells out of 92 were deleted这种提示，就说明当前体系原本有92个GTF壳层，其中6个冗余的被自动去除了。Gaussian的这种基组简化行为不限于Dunning相关一致性基组，而是对任何基组，哪怕自定义基组，都会自动这么搞，这使得Gaussian的计算结果和其它程序在某些情况下有不可忽略的差异，有时候能量差异能达到0.0001 Hartree甚至更多。如果想让Gaussian完全不做基组简化，可以用int=nobasistransform关键词。如果只允许做完全不影响精度的第一级简化，可以用int=ExactBasisTransform关键词。做第二级简化用的收缩系数阈值默认是1E-4，可以通过int=BasisTransform=N关键词改为1E-N。注意，哪怕只做第一级这种精度无损的简化，也会令基函数的轨道展开系数发生变化，因此如果将Gaussian得到的波函数用于其它程序当做初猜，最好加上int=nobasistransform，这点在《将Gaussian等程序收敛的波函数作为ORCA的初猜波函数的方法》（http://sobereva.com/517）中特意提到了。顺带一提，在著名的BSE基组数据库里有个Optimized General Contractions选项，选了这个之后得到的基组就相当于是做了第一级简化之后的。
有的基组有新旧不同的版本，不同程序内置的可能不同。比如Chem. Phys. Lett., 825, 140575 (2023)中指出cc-pVnZ系列对于Li，在OpenMolcas/Molpro/PSI4和Gaussian/CFOUR/NWChem/Q-Chem里内置的定义存在差异，对结果有不可忽视的影响。
有些基组的定义有一定含糊性，要注意看手册里的说明，必要的时候让程序把基组定义输出出来进行检查、对比。值得一提的是，6-311G系列是Pople系列基组之一，但实际上Pople只是在J. Chem. Phys., 72, 650 (1980)中参与了这个系列基组的前两周期的定义，而Gaussian程序中的6-311G用于其它元素的定义是来自于其他不同研究者的多篇不同文章的，是为了让6-311G系列能完整涵盖前四周期所有元素而被Gaussian开发者凑在了一起（详见手册）；另外，6-31G、6-311G系列基组对于He的定义都没有公开发表，只在Gaussian程序里有。因此，其它程序里对某些元素用的Pople系列基组的定义可能和Gaussian程序有出入。BSE基组库上的6-311G甚至对碘有定义，这也是网站的维护者把其他人搞的碘的基组冠上了6-311G系列的名字。
有些程序中内置的某些非主流的基组的定义可能比较特殊。例如在《给基组以even-tempered方式增加弥散函数的工具adddiffuse》（http://sobereva.com/347）中笔者提到，Gaussian里的d-aug-cc-pVnZ系列基组的最外层弥散函数的指数来路不明，和BSE基组数据库上的基于even-tempered方式构建的不符。
有的程序在使用某些基组上有特殊的习俗。比如在GAMESS-US程序中，基组名写CCD、CCT，看起来似乎应当是分别对应cc-pVDZ、cc-pVTZ，但实际上用的是cc-pV(D+d)Z、cc-pV(T+d)Z。对于Al~Ar带不带+d的版本是明显不同的。所以使用时要注意看手册。
当基函数存在明显线性依赖的情况时，可能造成计算过程数值不稳定、结果出现异常，甚至任务失败。量子化学程序通常根据基函数的重叠矩阵的本征值来判断是否存在显著线性依赖，如果有的本征值过小，则说明基函数存在很大程度的线性依赖问题，此时会自动砍掉一些线性依赖基函数。通常对于体系较大且使用弥散函数较多的情况线性依赖问题会比较显著。不同程序去除去除线性依赖基函数的阈值不同（Gaussian里可以通过IOp(3/59)设置），这造成某些情况下，特别是用大量弥散函数时不同程序的结果可能会有一些出入，甚至可能带来毫Hartree数量级的差异。
有些非主流计算程序使用的不是主流的Gauss函数作为基函数，比如ADF用的是STO基函数计算，FHI-aims、Dmol3等用的是数值基组，这些程序的计算结果注定不可能与使用Gauss基函数的主流量子化学程序去严格对比。笔者在网上答疑时有时看到有些人非要试图去严格重复那些程序的数据，这根本没什么意义，本来就不可能精确重复。实际上只要自己用的基组档次和文章中用那些程序计算时候用的差不多，得到的结论就会是相同的，这便足够了。比如有人用ADF拿TZP基组算的，你用def2-TZVP也能得到差不多的研究结论。
1.2 理论方法的差异性
不同程序做后HF、双杂化泛函计算时用的冻核设定往往不同。比如Gaussian和ORCA默认是冻核的，CFOUR、PSI4、NWChem默认不冻核，GAMESS-US默认冻不冻核看具体模块。而且不同程序默认用的冻核规则也往往不同，比如Gaussian默认只保留价层轨道，但从K开始的碱/碱土金属原子还保留亚价层的s,p轨道，当使用6-31G/6-311G系列基组时还保留所有p族元素的亚价层d轨道。不同程序具体怎么冻核的大家应当注意看手册、对比程序输出的相关信息。
不同的程序用的理论方法在定义上可能存在差异，典型的是B3LYP定义的差异。在B3LYP定义之初没说清楚其VWN局域相关泛函用的具体形式，导致不同程序用了不同的形式，有人在Chem. Phys. Lett., 268, 345 (1997)中还做了专门讨论和对比测试。Gaussian中默认用的是VWN3，TURBOMOLE、GAMESS-US、CRYSTAL、ORCA等默认用的是VWN5。如果要用Gaussian的B3LYP，在ORCA里应当写B3LYP/G，在GAMESS-US里应当写B3LYPV1R（而非B3LYPV3，这点经过实测确认），在Dalton里应该写B3LYPg。
有的程序在实现某些方法的时候对原方法做了修改。比如在Gaussian里直接写PM7关键词，实际上用的是Throssel和Frisch修改的PM7，使得势能面更连续。如果要用Stewart原始定义的PM7，需要写PM7MOPAC（经过实测，两种关键词的能量结果并没区别，Gaussian在实现时具体改了什么，由于相应文章目前还没有发表，不得而知）。
有的程序由于一些莫名其妙的规则，实际计算的和你想象的可能明显不同。比如在Gaussian 16中，写PBE0关键词给你算的不是PBE0泛函，而是PBE0-DH双杂化泛函。这点笔者在http://bbs.keinsci.com/thread-13660-1-1.html中专门提醒了。
对于开壳层体系计算，不同程序默认的形式也可能不同。比如Gaussian、ORCA对于自旋多重度大于1的体系默认都是用U的形式，而Molpro默认用RO。Molpro里的开壳层耦合簇关键词也有一定迷惑性，此程序里的UCCSD不是像多数程序基于UHF参考态做的CCSD，而是基于ROHF做的（对于开壳层体系，Molpro里还有个RCCSD，这也是基于ROHF参考态，但在做CCSD的时候给振幅施加约束，使波函数线性部分是自旋平方算符的本征函数）。
限制性开壳层计算（指ROHF、ROKS计算）有不同的具体实现，不同程序往往采用不同算法，虽然总能量相同，但轨道能量截然不同，在J. Chem. Phys., 125, 204110 (2006)的表1的对比中就充分体现了这一点，这是值得注意的一个问题。另外，基于ROHF轨道的后HF方法也有不同的具体实现，因此不同程序给出的结果可能明显不同，比如基于ROHF的MP2有ROMP、RMP、OPT、IOPT、ZAPT、HCPT、OSPT等等一大堆。对比之前应当先看手册、看其中引用的具体文献，弄清楚不同关键词写法对应的是哪些文章提出的方法，不要太想当然。
DFT-D色散校正是如今经常使用的，之前笔者在《DFT-D色散校正的使用》（http://sobereva.com/210）、《谈谈“计算时是否需要加DFT-D3色散校正？”》（http://sobereva.com/413）、《乱谈DFT-D》（http://sobereva.com/83）中详细介绍过DFT-D3。要注意不同程序，哪怕相同的程序的不同版本，用的D3参数也有可能不同。参数比较混乱的典型是B2PLYP，这个泛函总共有4套D3参数：
(1)在Grimme的2010年DFT-D3的原文的表4中，B2PLYP的DFT-D3零阻尼参数是S6=0.5，Sr,6=1.332，S8=1.000
(2)在Grimme的PWPB95-D3的原文J. Chem. Theor. Comput., 7, 291 (2011)中，作者重新拟合了B2PLYP的DFT-D3零阻尼参数，由表3可见参数为S6=0.64，Sr,6=1.427，S8=1.022
(3)在Grimme的专门对比零阻尼和BJ阻尼形式的文章J. Comput. Chem., 32, 1456 (2011)的表2中，B2PLYP的DFT-D3(BJ)参数为S6=0.5，S8=1.0860，a1=0.3451，a2=4.7735
(4)在Grimme的GMTKN30测试集文章Phys. Chem. Chem. Phys.,13, 6670 (2011)的表S1中，B2PLYP的DFT-D3(BJ)参数为S6=0.64，S8=0.9147，a1=0.3065，a2=5.0570
可见哪怕同一种阻尼函数都有两套参数，发表年份还这么接近，不导致不同程序出现分歧才怪，所以互相对比时要看清楚输出信息中有无提示具体参数、看看手册里的参数引的哪篇文章。在Gaussian中，B2PLYP加DFT-D3(BJ)校正有两种写法，一个是B2PLYPD3，一个是B2PLYP em=GD3BJ，然而在G09中这两种写法的色散校正值结果不同，估计就是因为用的D3(BJ)参数不同所致，但手册里却没提。而在G16中这两种写法的结果则变得统一了，都等于G09中B2PLYPD3的写法，此时用的参数对应的是上面的(4)。
1.3 积分计算的差异性
笔者在《密度泛函计算中的格点积分方法》（http://sobereva.com/69）中详细介绍过几乎所有能做DFT的程序都用的交换-相关泛函的多中心格点积分算法。虽然算法是共通的，但是不同程序中的积分格点有所不同，导致DFT结果注定不可能严格对比。差异有这些方面：
(1)径向和角度的积分格点数。点数和元素所属的周期还有关系，很多程序中是对于周期越靠后的元素用的点数越多
(2)径向格点位置的产生算法。比如常见的有第二类Gauss-Chebyshev积分格点（ORCA用的）、Euler-Maclaurin积分格点（GAMESS-US默认的）等等
(3)角度格点位置的产生算法。这个在不同程序之间倒是差异不大，多数都是用Lebedev格点，但也有其它做法，比如GAMESS-US里也可以用Gauss-Legendre积分格点得到角度部分格点的极坐标位置
(4)格点的剪裁。为了降低耗时，可以让距离原子核较近的区域用较少的角度部分格点，可以去掉权重很小的格点。不同程序做法有所不同
不是说不同程序的DFT结果就一定不能比，而是要比的话必须都用很高质量的积分格点（比如Gaussian用int=superfine，ORCA用grid7），这样才能让积分格点的差异带来的影响最小化。相对于G16默认用的ultrafine精度的格点，ORCA默认的积分格点非常糙，因此默认设置下二者算的DFT能量没什么可比性。另外，也有人提出过标准化的DFT积分格点，典型的是SG-1，其作者希望这种积分格点就像标准化的基组一样，能让不同量子化学程序做DFT的结果有可比性。不过支持它的程序不算很多，Gaussian、Q-Chem等程序支持，Q-Chem还支持后来的SG-2和SG-3。
对于HF、普通泛函的计算，基函数间的双电子积分的计算是计算耗时的主要来源。一般程序都会为了节约时间使用积分屏蔽策略，若通过简单关系估计出某双电子积分数值非常小，或者结合密度矩阵判断某双电子积分对Fock矩阵的贡献非常小，就直接将之忽略而不去实际计算。在程序默认设置下，这种处理通常对于结果精度影响是微乎其微的，不过当使用大量弥散函数算大体系时影响可能不容忽略（阈值太宽松还容易导致SCF难收敛）。这种做法也可能明显影响不同程序计算的能量间的精确比较。不同程序用的忽略双电子积分的方式和阈值不同，在Gaussian里可以用int=acc2e=N来设置阈值为1E-N，G09是N=10，G16是N=12，N越大忽略得越少、计算耗时越高。在ORCA中的默认设置等价于%scf Thresh 1e-8 TCut 1e-10 end，其中Thresh控制忽略双电子积分的阈值，TCut控制忽略primitive双电子积分的前因子的阈值，详见手册的Integral Handling部分，二者设得越小忽略的积分越少。当使用弥散函数特别多时，在对比两个程序计算的能量的时候建议把积分忽略的阈值设得严一些。
值得一提的是，Terachem这个程序是专注于GPU加速做量子化学计算的程序，见《首个完全基于GPU的量化软件-TeraChem杂谈及真实性能测试》（http://sobereva.com/137）。由于目前消费级市场的GPU的单精度浮点性能远远超过双精度，但是量子化学计算对于精度要求较高，又不能像简单粗暴的基于经典力场的分子动力学模拟那样可以只用单精度浮点数，因此Terachem支持单-双精度的混合，在尽可能利用GPU的单精度性能的同时又让结果尽量接近于纯双精度计算的结果。具体来说，Terachem和一般量化程序一样，对于数值特别小的积分直接忽略，而对于数值不可忽略但又不太大的积分用单精度计算，而对于数值较大、对结果精度影响明显的积分用双精度计算。无疑，这种做法得到的结果肯定和其它量子化学程序算的有一定差异，不过考虑得当的话可以令差异在可接受范围内。
有的程序默认用密度拟合技术加速计算。比如ORCA做纯泛函计算的时候，默认就会用RIJ技术大幅节约耗时，这会给结果带来不可忽略的差异，影响有多大具体看辅助基组的大小，相关知识这里介绍过：《大体系弱相互作用计算的解决之道》（http://sobereva.com/214）。在ORCA中如果想关闭这种做法，可以写noRI。
1.4 收敛问题
不同程序在SCF迭代方面的算法不同、帮助收敛的策略不同、初猜不同，可能导致SCF收敛到不同波函数，导致最后能量没有可比性。如果你发现两个程序在严格可比的条件下（本文说的其它因素都考虑到了），一个算出来的能量比另一个更高，那能量较高的这个必定对应于不稳定波函数，即在轨道展开系数的空间下还可以找到能量更低的解。对于SCF做得不太好的程序，碰上电子结构比较复杂的情况，收敛到不稳定波函数的情况比较多见。为了解决这个问题，可以尝试程序支持的不同的初猜方法，可以尝试恰当设置片段组合波函数（见《谈谈片段组合波函数与自旋极化单重态》http://sobereva.com/82）或调整初始轨道的占据方式，也有的程序可以自动尝试优化能量更低波函数（比如Gaussian里的stable=opt），也可以试图把得到了能量更低解的另一个程序的收敛的波函数给当前程序当初猜用。比如利用Multiwfn，可以很容易地实现GAMESS-US、Gaussian、ORCA等程序得到的波函数的互用，见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）的“Multiwfn的文件格式转换功能”一节。
还要注意计算程序是否都真的收敛了。比如Gaussian 16 C.01和ORCA，做普通DFT泛函的单点任务，哪怕SCF没收敛，程序也不会报错。如果到最大迭代次数的时候能量还在显著震荡中，显然能量不能和其它程序比。
不同程序的SCF收敛限、CI计算的Davidson迭代收敛限、CCSD方程迭代的收敛限以及它们的判据都有所不同，如果有的程序收敛限特别宽松的话，结果也难以和其它程序准确对比。这里说一下Gaussian和ORCA的情况。Gaussian从09开始默认的SCF收敛限（tight）已经特别严了，虽然要求能量变化小于1E-6 Hartree，但由于对密度矩阵变化也有较严格的收敛要求，实际达到收敛的时候经常能量变化都小于1E-10 Hartree了。而ORCA算单点任务的收敛限默认是normalSCF，对总能量变化要求也是小于1E-6 Hartree，但并没有对波函数的收敛要求，所以默认设置下ORCA的SCF收敛要求远比Gaussian宽松。虽然看似能量收敛到1E-6 Hartree就已经挺精确了，但要注意仅达到这个条件的话，波函数可能仍未充分收敛，如果当前任务是后HF、双杂化泛函计算的话，可能因此导致最终的能量并不太准确，有可能和Gaussian算的结果有不可忽视的差异。所以我建议ORCA的用户做涉及到后SCF类型的计算都用tightSCF，比normalSCF的能量收敛限严了两个数量级。
2 其它方面的差异
计算涉及到更多因素的时候，导致不同程序间结果存在差异的来源就更多了，这里说一些常见的。
做TDDFT的时候，ORCA等程序默认开了TDA近似，这会对激发能、振子/转子强度等跃迁性质有明显影响，详见《乱谈激发态的计算方法》（http://sobereva.com/265）。此外，TDDFT一般都是通过Davidson迭代方式求解的，不同程序用的收敛阈值有所不同，而且不同程序默认求解出的态数往往不同，会导致结果有一定差异。
计算热力学数据的时候，ORCA会自动用Grimme提出的准RRHO方法考虑低频模式的贡献，如果你读的是这部分输出的热力学量，当体系存在低频模式时，会和一般程序基于RRHO模型得到的不同。什么叫RRHO、准RRHO，在这里有介绍：《使用Shermo结合量子化学程序方便地计算分子的各种热力学数据》（http://sobereva.com/552），其中提到的笔者的Shermo程序的介绍论文里有更充分说明。此外，对一个结构非常接近但不是严格满足点群对称性的结构，不同程序由于判断的点群有可能不同，对应的对转动对称数因此也可能不同，最终导致转动对热力学量的贡献可能不同。另外，不同程序计算热力学数据时用的同位素设置可能不同，比如Gaussian、NWChem、GAMESS-US都是用的丰度最大的同位素质量，但ORCA用的是丰度平均的元素的质量，这也会轻微影响结果。
不同程序计算的振动频率也会多多少少有所不同。算频率需要先计算Hessian矩阵，这牵扯求解CPHF方程，不同程序有不同的CPHF迭代求解的收敛限。对于DFT计算，CPHF计算还牵扯到积分格点，不同程序用的格点也有所差异。Gaussian算振动频率前会先把平动、转动模式对Hessian的影响投影掉，而GAMESS-US则没有这一步，这点也会导致得到的频率有所不同。频率的差异也影响振动对热力学量的贡献。另外，由于算振动频率需要对角化质量权重的Hessian矩阵，因此不同程序用的原子质量的不同也会给结果带来差异。
在相对论哈密顿下的计算时，不同程序默认的对核电荷的考虑的不同会给结果带来差异。比如做DKH2计算的时候，Gaussian、Dirac默认用高斯分布的有限核模型，即把原子核电荷用高斯函数描述，而GAMESS-US则用点电荷模型。
同一种隐式溶剂模型在不同程序、不同版本中有不同的实现，这导致不同程序，乃至同一个程序的不同大版本之间的结果是很难严格对比的。比如SMD模型，其极性部分在Gaussian中用的是IEFPCM第二版，而在ORCA（对于4.2版而言）中用的是CPCM，结果自然会不同。哪怕都是IEFPCM第二版，在G09用的是非对称形式，而G16是对称形式，结果也存在差异。哪怕都是CPCM，Gaussian和ORCA的实现细节也有所不同，比如孔洞构造方式、孔洞表面的离散化，等等。另外，Gaussian的CPCM模型对于中性体系用的参数还是错的，这点在J. Chem. Theory Comput., 11, 4220 (2015)中明确指出了。此外，ORCA较早版本用的是点电荷方式描述表面显著电荷，而4.2开始用的是数值更稳健的高斯函数形式描述，这使得同种溶剂模型在不同版本中的结果有所不同。另外，对于电子激发计算，不同程序对于溶剂的快、慢部分响应的考虑也有所不同，比如Gaussian和ORCA在隐式溶剂模型下做TDDFT算出来的某些态的激发能有时有不小差异，这点笔者在《使用ORCA在TDDFT下计算旋轨耦合矩阵元和绘制旋轨耦合校正的UV-Vis光谱》（http://sobereva.com/462）中提到过。
其它要素都严格相同的情况下，不同程序做几何优化任务得到的结构差异一般都在收敛限范围的差异内。但对于势能面复杂的大体系，以某些结构做初猜的时候，不同程序可能会优化到明显不同的极小点去，毕竟不同程序的优化算法、细节设置有别。
当输入、输出信息涉及到非原子单位的量的时候，由于不同程序或者同一个程序不同版本内置的物理常数、转换因子、元素质量可能有轻微差异，这也会导致结果可能有可查觉的不同。在Gaussian中专门有个关键词constants可以设置物理常数的版本，对于G16来说默认用的是constants=2010，而G09对应的是constants=2006。
顺带一提，G16有个关键词g09default，可以让结果和G09严格相同，这个关键词等价于int(fine,acc2e=10) constants=2006 scrf=G09Defaults。
程序bug也是导致不同程序间或者不同版本间结果无法对比的可能原因。如果经严谨的测试，且仔细阅读手册之后，可以确认这就是程序的问题的话，应当发邮件、去官方的邮件列表或论坛咨询程序开发者。
且不说不同程序之间对比，哪怕是同一个程序同一个版本、同一台机子，如果用了并行计算，两次计算的结果也会有非常细微的不同，原因在《数值误差对计算化学结果重现性的影响》（http://sobereva.com/88）中说了。这个问题对于程序间的能量可比性没什么直接影响，稍微了解一下就行。