談談不同量子化學程序計算結果的差異問題
談談不同量子化學程序計算結果的差異問題
文/Sobereva@北京科音
First release: 2020-Nov-1 Last update: 2023-Jun-18
經常有人問,為什么他用的兩個不同量子化學程序在相同級別計算下計算結果存在差異,也常有人問用不同的量子化學程序算的結果能否放在一起對比。在此文就專門說一下導致不同程序間,也包括相同程序的不同大版本間計算結果存在差異可能的因素。首先說最簡單的情況,即真空下單點能的計算,之后再說其它任務、涉及額外因素的情況。由于Gaussian和ORCA是最流行的兩個量子化學程序,所以對此二者說得相對多一些。
1 導致真空下計算的單點能存在差異的因素
1.1 基組的差異性
角動量>=D的Gauss型基函數分為球諧型和笛卡爾型兩種,介紹見《談談5d、6d型d殼層基函數與它們在Gaussian中的標識》(http://www.shanxitv.org/51)和《球諧型與笛卡爾型Gaussian函數的轉換關系》(http://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/210)、《談談“計算時是否需要加DFT-D3色散校正?”》(http://www.shanxitv.org/413)、《亂談DFT-D》(http://www.shanxitv.org/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://www.shanxitv.org/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://www.shanxitv.org/137)。由于目前消費級市場的GPU的單精度浮點性能遠遠超過雙精度,但是量子化學計算對于精度要求較高,又不能像簡單粗暴的基于經典力場的分子動力學模擬那樣可以只用單精度浮點數,因此Terachem支持單-雙精度的混合,在盡可能利用GPU的單精度性能的同時又讓結果盡量接近于純雙精度計算的結果。具體來說,Terachem和一般量化程序一樣,對于數值特別小的積分直接忽略,而對于數值不可忽略但又不太大的積分用單精度計算,而對于數值較大、對結果精度影響明顯的積分用雙精度計算。無疑,這種做法得到的結果肯定和其它量子化學程序算的有一定差異,不過考慮得當的話可以令差異在可接受范圍內。
有的程序默認用密度擬合技術加速計算。比如ORCA做純泛函計算的時候,默認就會用RIJ技術大幅節約耗時,這會給結果帶來不可忽略的差異,影響有多大具體看輔助基組的大小,相關知識這里介紹過:《大體系弱相互作用計算的解決之道》(http://www.shanxitv.org/214)。在ORCA中如果想關閉這種做法,可以寫noRI。
1.4 收斂問題
不同程序在SCF迭代方面的算法不同、幫助收斂的策略不同、初猜不同,可能導致SCF收斂到不同波函數,導致最后能量沒有可比性。如果你發現兩個程序在嚴格可比的條件下(本文說的其它因素都考慮到了),一個算出來的能量比另一個更高,那能量較高的這個必定對應于不穩定波函數,即在軌道展開系數的空間下還可以找到能量更低的解。對于SCF做得不太好的程序,碰上電子結構比較復雜的情況,收斂到不穩定波函數的情況比較多見。為了解決這個問題,可以嘗試程序支持的不同的初猜方法,可以嘗試恰當設置片段組合波函數(見《談談片段組合波函數與自旋極化單重態》http://www.shanxitv.org/82)或調整初始軌道的占據方式,也有的程序可以自動嘗試優化能量更低波函數(比如Gaussian里的stable=opt),也可以試圖把得到了能量更低解的另一個程序的收斂的波函數給當前程序當初猜用。比如利用Multiwfn,可以很容易地實現GAMESS-US、Gaussian、ORCA等程序得到的波函數的互用,見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》(http://www.shanxitv.org/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://www.shanxitv.org/265)。此外,TDDFT一般都是通過Davidson迭代方式求解的,不同程序用的收斂閾值有所不同,而且不同程序默認求解出的態數往往不同,會導致結果有一定差異。
計算熱力學數據的時候,ORCA會自動用Grimme提出的準RRHO方法考慮低頻模式的貢獻,如果你讀的是這部分輸出的熱力學量,當體系存在低頻模式時,會和一般程序基于RRHO模型得到的不同。什么叫RRHO、準RRHO,在這里有介紹:《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》(http://www.shanxitv.org/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://www.shanxitv.org/462)中提到過。
其它要素都嚴格相同的情況下,不同程序做幾何優化任務得到的結構差異一般都在收斂限范圍的差異內。但對于勢能面復雜的大體系,以某些結構做初猜的時候,不同程序可能會優化到明顯不同的極小點去,畢竟不同程序的優化算法、細節設置有別。
當輸入、輸出信息涉及到非原子單位的量的時候,由于不同程序或者同一個程序不同版本內置的物理常數、轉換因子、元素質量可能有輕微差異,這也會導致結果可能有可查覺的不同。在Gaussian中專門有個關鍵詞constants可以設置物理常數的版本,對于G16來說默認用的是constants=2010,而G09對應的是constants=2006。
順帶一提,G16有個關鍵詞g09default,可以讓結果和G09嚴格相同,這個關鍵詞等價于int(fine,acc2e=10) constants=2006 scrf=G09Defaults。
程序bug也是導致不同程序間或者不同版本間結果無法對比的可能原因。如果經嚴謹的測試,且仔細閱讀手冊之后,可以確認這就是程序的問題的話,應當發郵件、去官方的郵件列表或論壇咨詢程序開發者。
且不說不同程序之間對比,哪怕是同一個程序同一個版本、同一臺機子,如果用了并行計算,兩次計算的結果也會有非常細微的不同,原因在《數值誤差對計算化學結果重現性的影響》(http://www.shanxitv.org/88)中說了。這個問題對于程序間的能量可比性沒什么直接影響,稍微了解一下就行。