談談重復不出來計算化學文獻里的數據的可能原因
談談重復不出來計算化學文獻里的數據的可能原因
文/Sobereva@北京科音 2023-Aug-1
0 前言
經常有做計算化學的人在網上問為什么他沒法重復出來文獻里的計算數據。導致重復不出來的原因實在太多了,然而初學者提問此類問題時幾乎總是含糊其辭,提供的有效信息甚少,導致別人根本無從回復,或者需要別人羅列一大堆可能導致差異的原因。我感到每次回復此類問題太費事,于是在這里專門寫個文章說說這事。遇到這種問題的讀者仔細看完此文后,應該能料想出可能是什么原因導致自己和文獻的結果存在差異。就算還搞不明白,至少也應該能明白在網上問別人的時候應當提供什么信息,知道怎么把導致差異的可能原因的范圍盡可能縮小,免得別人打很多字羅列諸多當前不需要考慮的可能性,或者別人嫌你提供的有效信息完全不夠而根本不想回復。
這里首先要強調一點的是,如果你的目的就是為了精確重復文獻數據,那你就要盡最大可能使用和文獻完全相同的方式去做計算,而不要管人家的做法合不合理,哪怕你明知道人家的計算方式不當。下文說的內容都是為了讓你能盡可能準確重復出文獻里的數據,而不是講怎么合理地計算。如果你的實際目的是想獲得相同體系相同研究問題的合理的結果,那你只要確保你自己的計算是合理的就行了,沒必要非得把別人的數據重復出來,說不定人家的計算本來就有問題。如果最終驗證發現別人的計算有硬傷并因此導致結論錯誤,或者憑基本知識就能斷定別人的計算有問題,那么可以發一篇comment文章指正以免以訛傳訛,比如《我對一篇存在大量錯誤的J.Mol.Model.期刊上的18碳環研究文章的comment》(http://www.shanxitv.org/584)里介紹的文章。
1 結構不一樣
自己計算用的結構和文獻不一致,是導致數據無法重現的最常見原因之一。如果你用的結構就是文獻正文或補充材料里給的坐標,那就可以直接排除這個因素了。注意有些文獻的SI里鍵長、笛卡爾坐標單位是Bohr而不是埃,千萬別弄錯了。
如果你計算的基本結構特征都和文獻里不符,也即算的都不是同一個東西,能重現出文獻里的數據就怪了。很多體系有不同的形態,比如酚酞在酸性、中性、堿性環境下主要的存在形態明顯不一樣,氨基酸在真空和水環境下的質子化態明顯不一樣,姜黃素在水和有機溶劑中一種是烯醇式一種是酮式結構,等等。你必須確保你算的和文獻里算的確實是同一種狀態的結構,如果文獻里給了截圖,一定要仔細對照。
做計算要有基本的化學常識,并且要認真仔細,免得搭出來沒意義的結構去做計算,自然和文獻里對不上。比如文獻里給出Lewis結構式時通常是隱去碳上的氫的,如果計算的時候真不畫氫,或者少畫了個別的氫,自然計算毫無意義。再比如,文獻里算一個羧基配位的過渡金屬配合物,如果你在計算時羧基上還留著氫,明顯也是錯的。再比如,文獻里畫了一個帶一個小配體的過渡金屬配合物的Lewis式,若你算的是水溶液中的狀態并忽略了配位水分子,也自然是錯的。再比如,文獻里給出一個聯苯或者碗烯的二維結構式,看上去是在一個平面上,若你去計算其基態時真的擺成平面來優化,得到的也是大概率是錯誤的結構(有破壞平面的虛頻的結構)。還要注意一些結構的構建不要想當然,比如SiO2表面在水中時表面的O主要是以羥基及其去質子化形態存在的,比例和pH有關,如果都當成是烏禿禿的氧就錯了。結構構建的關鍵性信息在文獻里一般都會有具體說明,但也可能一些文獻中在描述上有疏漏。
有些體系在計算的時候結構構建方式不唯一,比如《要善用簇模型,不要盲目用ONIOM算蛋白質-小分子相互作用問題》(http://www.shanxitv.org/597)、《18碳環(cyclo[18]carbon)與石墨烯的相互作用:基于簇模型的研究一例》(http://www.shanxitv.org/615)和《使用量子化學程序基于簇模型計算金屬表面吸附問題》(http://www.shanxitv.org/540)里提到的簇模型,怎么構建簇明顯會影響結果,特別是當簇用得較小的時候。若要重復文獻里的數據,就要用盡可能相同的模型,哪怕其模型很昂貴。
很多分子有不止一種構象(勢能面極小點),柔性大分子甚至有一萬種以上的構象。幾何結構優化一般會優化到與初始結構最接近的勢能面極小點。而不同構象下計算的屬性可能有不小差異,因此你必須確保你用的構象和文獻里一致,也即結構優化使用的初猜結構要盡可能像文獻里給出的。如果文獻沒明確給出三維結構信息,按理說文獻應當用的是能量最低構象(確切來說是自由能最低構象)做的計算,但也有可能作者忽略了構象搜索并得到了不合理的結果。如果你對構象搜索一點概念都沒有,建議看下文了解些常識
《使用molclus程序做團簇構型搜索和分子構象搜索》(http://bbs.keinsci.com/thread-577-1-1.html)
《gentor:掃描方式做分子構象搜索的便捷工具》(http://bbs.keinsci.com/thread-2388-1-1.html)
《將Confab或Frog2與Molclus聯用對有機體系做構象搜索》(http://bbs.keinsci.com/thread-20063-1-1.html)
《使用Molclus結合xtb做的動力學模擬對瑞德西韋(Remdesivir)做構象搜索》(http://bbs.keinsci.com/thread-16255-1-1.html)
也不是所有情況都適合用能量最低構象來計算,比如某化學反應可能是從某個能量相對較高的構象發生的,可以對反應的IRC的反應物一端做幾何優化得到之。
團簇體系的構型搜索和柔性分子的構象搜索同等重要。哪怕水二聚體這么簡單的兩個小分子構成的團簇都都有不同的極小點構型,并且對應的結合能明顯不同。因此計算這種體系也必須確保和文獻用的構型相一致。如果文中沒明確交代,一般默認其用的是能量最低構型結構做的計算,但也可能作者忽略了這個問題,導致誤用了能量不是最低的構型并因此得到了錯誤的結果。如果你對構型搜索沒有概念,建議看此文:《genmer:生成團簇初始構型結合molclus做團簇結構搜索的超便捷工具》(http://bbs.keinsci.com/thread-2369-1-1.html)。
有的時候兩種極小點構型或構象可能結構相差很小。這種情況需要把可視化程序里的視角和文獻里的圖片擺得盡可能一致再去觀察,免得算的其實不是同一個結構。
上面說的是計算極小點結構的情況,自不必說,重復文獻里過渡態的計算結果時也必須嚴格確保和文獻里用的結構一致。
對于研究固體表面的情況,要注意文獻里的結構是怎么來的。表面是從體相直接切出來的未經優化的(unrelaxed),還是經過優化的(relaxed),還是使用的考慮表面重構的(reconstructed),這需要區分清楚。文獻中計算用的是具體哪個晶面也必須要搞清楚。固體表面吸附問題也有類似構象/構型搜索的勢能面存在多種極小點問題,一個小分子可能在表面上不同位點以不同方式吸附,實際算的結構必須和文獻相一致。文獻里算的吸附是物理吸附還是化學吸附也得搞清楚。對于吸附過程有勢壘的情況,初猜結構中你把被吸附物和表面擺得相距很近,一般優化完了的是形成新鍵的化學吸附,如果擺得較遠,一般得到的是弱相互作用維持的物理吸附,二者的結構、相互作用能相差甚大。研究表面吸附還有覆蓋率的問題,要注意文獻中對所用結構的相關描述,研究低覆蓋率和飽和覆蓋率時用的模型明顯不同。
用slab模型研究固體表面還要注意厚度問題。slab太薄的話,得到的功函數、吸附能等性質都不準確。對于二維層狀材料的研究,還要注意層數對計算的性質的影響。這些方面都應當和文獻嚴格一致。對表面類體系優化時還經常牽扯到對一部分原子凍結使之固定為體相的結構,之后振動分析也相應地牽扯凍結,凍結方式必須和文獻一致。
做某些體系的分子動力學模擬也同樣要注意初始結構的差異,否則現象可能和文獻不符。比如有文獻研究表明水中的蛋白質由于疏水效應會自發進入口徑足夠大的碳納米管,如果你一開始把蛋白質放得離碳納米管老遠,由于質量頗大的蛋白質本身整體運動就較慢,再加上碳納米管的隨意旋轉(沒凍結的話),可能文獻里說的那種蛋白質自發進入管內的現象跑挺長時間也出現不了。再比如想模擬結冰過程,哪怕一開始把溫度設得顯著低于冰點,但沒有在水中加入局部的冰的結構的話,跑很長時間也不會看到冰的結構的出現。至于在模擬的時間尺度內一定會出現的現象,初始結構就不重要了。比如水和乙醇的互溶,無論初始結構里把二者混合均勻,還是分成兩相,最后肯定都是互溶的。
做周期性體系的模擬還要注意盒子尺寸與文獻要一致。用第一性原理計算固體時,假設k點數是固定的,盒子尺寸的不同會影響精度。做動力學模擬時,盒子尺寸還間接影響原子運動行為以及模擬得到的屬性。
用Quantum ESPRESSO等平面波程序時,以及CP2K這種將平面波當輔助基函數的程序時,如果算的是某些維度為非周期性體系,還需要注意真空層的厚度,要盡可能和文獻一致。因為真空層可能對結果有不可忽視的影響,比如影響非周期性方向上與相鄰鏡像的作用是否可以忽略。而且比如CP2K里用MT的Psolver時還有非周期性方向上盒子尺寸大于相應方向有電子密度明顯分布的兩倍以上的要求,顯然真空區不能小了。
假設文獻用的是實驗測定的坐標直接進行性質的計算,還要注意實驗來源必須一致。比如文獻或晶體數據庫里對同一種晶體結構可能提供了不同溫度下測定的坐標,會有一定程度的不同(體現熱膨脹)。顯然壓力也直接影響測定的結構。哪怕是相同的測試環境,不同來源的晶體結構也會有差異,特別是氫的位置的確定,無疑這對氫鍵等問題的計算結果會有直接影響。做量子化學/第一性原理研究的人一般都會至少對氫的坐標重新進行優化來refine之,需要注意文獻里是怎么考慮的。蛋白質的高解析度和低解析度測出來的晶體結構里殘基側鏈的位置也可能相差甚大,甚至明顯影響到做酶催化、分子對接之類問題的結果。同樣的屬性,以不同類型實驗測定的結果注定也會有差異,比如氣相電子衍射和單晶衍射測的一個是氣相結構一個是晶體結構,對于柔性分子差異可能很大,這點參考《實驗測定分子結構的方法以及將實驗結構用于量子化學計算需要注意的問題》(http://www.shanxitv.org/569)。
上面列舉了諸多結構層面影響計算結果的可能性,當然不可能列舉得很完備,實際情況多種多樣。總之希望讀者能意識到嚴格保證和文獻里用的結構的一致對于重現數據有多么重要。
2 計算設置不一樣
為了精確重復文獻數據,要盡可能仔細閱讀文章及其補充材料,了解清楚文獻計算細節,使得自己的計算設置能最大程度與文獻相一致。本節里談到的很多可能導致差異的因素在《談談不同量子化學程序計算結果的差異問題》(http://www.shanxitv.org/573)里都有詳細說明,務必注意閱讀。
量子化學和第一性原理計算用的理論方法,以及理論方法中涉及的設定和參數(如范圍分離泛函的w參數、色散校正的經驗參數、CASSCF和多參考方法的活性空間的選擇、TDDFT計算是否用了TDA近似、DFT+U參數、CASPT2計算用的IPEA shift參數、DLPNO方法用的閾值、GW@DFT和DFT-SAPT具體基于什么泛函、sTDA計算用的經驗參數,等等)必須和文獻精確一致。有的理論方法在同一個程序里甚至有兩種定義,比如B3LYP在ORCA里寫成B3LYP還是B3LYP/G是不同的,要搞清楚文獻實際用的哪個(通常是前者)。還要注意理論方法名和關鍵詞的差異,比如里這里提及的:《Gaussian16里用PBE0關鍵詞計算的實際上是PBE0-DH雙雜化泛函》(http://bbs.keinsci.com/thread-13660-1-1.html)。有的文獻的寫法還可能存在含糊性,比如可能文獻就說對某個開殼層體系用了二階微擾理論,但實際上二階微擾理論用于開殼層的實現形式有一堆,比如UMP2、PUMP2、ROMP2、ZAPT2、ROMP、OPT2等等,如果文中給了引用的文獻,要根據文獻判斷具體是哪種。
對計算結果可能有明顯影響的各種條件必須一致,比如:
? 溶劑效應的考慮。如果用的是隱式溶劑模型,那么用的具體是什么方法?描述的是哪種溶劑?用CP2K的SCCS模型時α、β、γ參數是怎么設的?描述混合溶劑的時是怎么確定溶劑參數的?如果用的是顯式溶劑模型,那么溶劑分子排布是怎么考慮的?有沒有同時用隱式溶劑模型?
? 計算熱力學校正量的模型。比如《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》(http://www.shanxitv.org/552)里介紹的Shermo程序默認用的以及ORCA用的是Grimme的QRRHO模型,而Gaussian用的是較原始的不適合含低頻體系的RRHO模型
? 計算熱力學量時是否考慮了平動和轉動,顯然算周期性體系時不應當考慮(在Shermo程序中這由settings.ini里的imode決定)
? 計算轉動對熱力學量貢獻時對轉動對稱數的考慮,看Shermo手冊附錄部分了解相關常識
? 計算熱力學校正量、振動分析用的原子質量
? k點數目和分布
? 平面波相關計算的截斷能
? CP2K中,做雜化泛函時用的庫侖截斷半徑、DFT+U計算用的布居計算方法以及考慮的角動量、Poisson solver的選擇
? 對相對論效應的考慮
? 對外場/外勢的考慮
? DFT積分格點精度。參看《密度泛函計算中的格點積分方法》(http://www.shanxitv.org/69)
? 是否用了凍核、怎么凍的
? 雙電子積分的積分屏蔽設置
? 是否用了RI、DLPNO等加速計算方式
? 是否考慮了偶極校正
? QM/MM計算劃分的位置、MM與QM區之間的相互作用的考慮、link原子的設置
等等等等
基組必須和文獻嚴格一致,并注意上述http://www.shanxitv.org/573博文里說的許多細節問題。注意一個基組/贗勢名可能有不止一個指代,比如Stuttgart贗勢,有不同方式考慮相對論效應的版本、不同價電子數和氧化態的版本,Gaussian內置的和Stuttgart官網上的贗勢和贗勢基組定義往往又有所不同。再比如,GAMESS-US用戶寫CCD和CCT關鍵詞分別用的是cc-pV(D+d)Z和cc-pV(T+d)Z,而不是更常見的cc-pVDZ和cc-pVTZ。要弄清楚作者實際用的什么,仔細看文中的描述,并注意作者引的基組原文。
不光一般意義上的基組(叫主基組或者叫軌道基組)要和文獻計算時的一致,輔助基組也應當一致。輔助基組種類很多,比如RIJ用的/J輔助基組,RIJK用的/JK輔助基組,電子相關計算用的/C輔助基組,CP2K里加快雜化泛函計算用的ADMM輔助基組,CP2K里做LRIGPW計算用的輔助基組,等等。同一個主基組能夠或者適合搭配的某種目的的輔助基組可能不止一種,而且有的時候有含糊性和任意性,比如ORCA里用ma-def2-TZVP主基組時,可能有人用autoaux關鍵詞自動產生輔助基組,可能有的人用標準的def2/J輔助基組。
對于基于經典力場的模擬,力場參數顯然必須和文獻嚴格一致。如果參數是從文獻里拷來然后用在自己計算中的,注意單位轉換,以及1-4相互作用規則和LJ參數的混合規則等細節必須和文章相同。文獻給出的LJ參數有的是R0有的是σ需要區分。原子電荷也必須和文獻保持一致,如果原子電荷是作者自己算的,要搞清楚計算原子電荷的方法,并且如果是基于波函數算的,那么波函數是在什么幾何結構、溶劑環境、計算級別下得到的要弄清楚。這點很重要,比如對于RESP電荷,氣相下計算的和水環境下計算的,Hartree-Fock和B3LYP計算的,相差大了去了,參看這些討論:《RESP擬合靜電勢電荷的原理以及在Multiwfn中的計算》(http://www.shanxitv.org/441)、《RESP2原子電荷的思想以及在Multiwfn中的計算》(http://www.shanxitv.org/531)。
做動力學模擬的情況,必須確保模擬相關設置和文獻一致。比如步長、步數、坐標/速度/受力/能量保存頻率、參考溫度和壓力、控溫和控壓方法、熱浴和壓浴的時間常數、可壓縮系數、各項同性還是各向異性控壓、初速度怎么產生的、范德華和靜電作用的計算方式(如簡單截斷、twin-range、Ewald、PME、SPME等,以及其中涉及的截斷半徑,是否用了switching function等處理)、水模型是剛性還是柔性的、是否用了鍵/鍵角約束、用沒用能量/壓力色散校正,等等。
如果你用的計算程序和文獻里不同,或者和文獻里的相同但版本不同,那么由于默認設置不同,以及算法實現的不同,可能對結果帶來很大差異,一定要仔細閱讀《談談不同量子化學程序計算結果的差異問題》(http://www.shanxitv.org/573)。
在量子化學波函數分析領域也有很多細節設置會影響分析結果。這里以非常流行的波函數分析程序Multiwfn(http://www.shanxitv.org/multiwfn)為例:
? 按照《使用Multiwfn做電子密度、ELF、靜電勢、密度差等函數的盆分析》(http://www.shanxitv.org/179)、《使用Multiwfn的定量分子表面分析功能預測反應位點、分析分子間相互作用》(http://www.shanxitv.org/159)、《使用Multiwfn對靜電勢和范德華勢做拓撲分析精確得到極小點位置和數值》(http://www.shanxitv.org/645)等文章介紹的方法進行分析時,格點間距數值越小精度越高,但也越耗時
? 按照《使用Multiwfn繪制態密度(DOS)圖考察電子結構》(http://www.shanxitv.org/482)繪制PDOS和OPDOS圖時,在Multiwfn里選擇的計算軌道成份的方法會直接影響結果
? 按照《使用Multiwfn模擬掃描隧道顯微鏡(STM)圖像》(http://www.shanxitv.org/549)和《使用Multiwfn結合CP2K的波函數模擬周期性體系的隧道掃描顯微鏡(STM)圖像》(http://www.shanxitv.org/671)繪制STM圖時,偏壓的數值和常高模式掃描的平面的選擇都會明顯影響結果
? 計算《衡量芳香性的方法以及在Multiwfn中的計算》(http://www.shanxitv.org/176)和《使用AdNDP方法以及ELF/LOL、多中心鍵級研究多中心鍵》(http://www.shanxitv.org/138)里介紹的多中心鍵級時,基于原本的基函數和基于自然原子軌道(NAO)計算會得到不同的結果
? 計算Hirshfeld原子電荷時需要用到原子自由狀態下的密度。用自行提供的原子波函數文件來計算此密度,還是直接用Multiwfn內置的此密度,會有所不同。
本節提到的這些計算細節是理應充分體現在文章的computational details里或者補充材料里的。但由于一些文獻的作者為了省篇幅,或者用的很多設置和程序默認的一致,或者由于作者計算水平比較外行、意識不到很多細節對計算非常重要,就在文中沒特意說。如果有些要點不知道作者具體怎么考慮的,自己把各種可能性都試試就知道了,文章沒特意說的設置就姑且用程序默認的。
3 計算的電子結構狀態不一樣
凈電荷和自旋多重度的設置必須和文獻里的相同,否則相當于算的是其它電子態,自然算出來的各種性質都不同。比如計算基態的氧氣發生的反應、計算富勒烯的原子化能,如果氧氣和碳原子誤當成了單重態來算,顯然結果是錯的。此外,稍微了解點兩態反應性的概念就會知道,不少體系在不同自旋態下反應機理都明顯不一樣,如果自旋多重度或凈電荷設錯了,自然反應物/產物/中間體/過渡態等無法得到和文獻里一致的結果,而且相差懸殊。
如果文獻沒注明、從文獻中看不出相關信息,那么凈電荷一般為0。如果被研究的體系的自旋多重度不那么顯而易見,或者研究中涉及到了不止一個自旋態,文中一般都會說明自旋多重度。如果確實沒說,那就得查閱文獻或者數據庫,或者自己判斷、測試了。
很容易被忽略的是開殼層單重態或者叫自旋極化單重態。用KS-DFT等方法來計算的話應當做對稱破缺計算,參看《談談片段組合波函數與自旋極化單重態》(http://www.shanxitv.org/82)。比如整體為單重態的反鐵磁性耦合體系、雙自由基體系、共價鍵被拉長的均裂狀態就是最典型的例子。還有一些特殊的體系也是如此,比如18碳環的D18h點群的過渡態結構、乙烯二面角扭轉成90度的結構、六并苯。重復文獻里的這些體系的計算時別誤當成了閉殼層單重態來計算(文中如果也是做了對稱破缺計算,一般都會明確說。如果沒說,沒準文獻也誤當成了閉殼層計算)。
對于第一性原理計算磁性材料時必須謹慎,很多情況下不是光設個整體的自旋多重度就完事的,而必須恰當設置初始磁矩才能收斂到真正要考察的磁性狀態。為了重復文獻(或材料數據庫)中的數據,若它們給了收斂的波函數對應的各個原子的自旋磁矩/自旋布居信息的話,最好根據其來設置初猜波函數里的原子自旋磁矩,以盡最大可能收斂到和文獻相同的波函數狀態。
要注意波函數穩定性問題。即便凈電荷、自旋多重度的設置是正確的,也不代表SCF收斂后得到的波函數就是實際想研究的態、和文獻里算的狀態相一致。如果文獻里研究的是基態,你算的結果和文獻對不上,而且體系又不是簡單有機體系而是電子結構特殊的體系,那么應當做一下波函數穩定性檢測,比如Gaussian里可以用stable關鍵詞,發現不穩定的話可以用stable=opt嘗試搞出穩定波函數。不穩定的波函數可能意味著這是某個激發態,還可能意味著此波函數是沒有任何物理意義的虛假的狀態。
4 計算的東西或方式不一樣
重復不出文獻的數據時還要想清楚到底你算的東西、計算的方式是否和文獻一致、有可比性。下面列舉一些典型的需要注意的情況。
比如通過能量差E(AB)-E(A)-E(B)考察倆分子A和B間相互作用強度,復合物結構肯定是要優化的,而單分子是直接從復合物優化后的結構里直接取出來,還是也做優化,這在不同文中的做法就不一樣了。單分子不優化情況下得到的能量差一般稱為(復合物結構下的)相互作用能,單分子也優化情況下得到的能量差一般稱為結合能,二者之間相差各個單體的變形能之和。對于諸如《透徹認識氫鍵本質、簡單可靠地估計氫鍵強度:一篇2019年JCC上的重要研究文章介紹》(http://www.shanxitv.org/513)里面列舉的那些剛性小分子,相互作用能和結合能差異甚微,但如果考察的是兩個容易變形的分子間的作用強度,比如兩個氨基酸分子的結合,那么相互作用能和結合能就必須明確區分了。J. Chem. Phys., 158, 244106 (2023)還專門討論了什么情況變形能會較大。對于鍵能,也是類似地有不同計算方式。所以研究此類問題的時候一定要搞清楚文獻具體是怎么算的。
再比如做SAPT能量分解計算,會得到很多耦合項,比如交換與色散作用的耦合項,這是歸到色散項里,還是歸到交換互斥項里,不同文章的做法可能不一樣,需要注意文章怎么說的。再比如計算反應的自由能問題,是計算標準反應自由能,還是計算特定濃度下的反應自由能,這是不一樣,相差濃度校正項,要搞清楚文獻實際算的是什么。再比如電子親和能、電離能、激發能,都有垂直和絕熱之分,文獻里計算哪種都有,也必須確保和文獻算的東西對應。
文獻里有時候對熱力學量的標注比較含糊,要分清楚文獻算的是什么溫度下什么量。比如E,文獻里可能指電子能量,也可能指內能,需要結合語境判斷。
模擬ECD譜需要轉子強度,有長度形式和速度形式兩種,諸如Gaussian都會輸出,在Multiwfn里按照《使用Multiwfn繪制紅外、拉曼、UV-Vis、ECD、VCD和ROA光譜圖》(http://www.shanxitv.org/224)里說的繪制ECD時可以選采用哪種,一般推薦用速度形式。如果你用的和文獻不一致,會導致峰型存在一定差異。
基組的CBS外推的方法不唯一,有不同公式,有的情況外推方式里也涉及到參數,參看《談談能量的基組外推》(http://www.shanxitv.org/172)。不光是能量外推,還有比如GW計算的準粒子能量隨基組的外推、幾何結構隨基組的外推等。一些文獻對高級別理論方法結合大基組還用一些近似方式估計,比如一種常見的是CCSD(T)/大基組 = CCSD(T)/小基組 + (MP2/大基組-MP2/小基組),這在《各種后HF方法精度簡單橫測》(http://www.shanxitv.org/378)里專門說了。涉及到這些高精度數據估算時,需注意保證和文獻的做法一致。
還有一些量本來就沒有唯一定義,去重復文獻的數據更是得弄清楚文獻里用什么定義、什么方法算的。比如自由體積,我在《使用Multiwfn計算晶體結構中自由區域的體積、圖形化展現自由區域》(http://www.shanxitv.org/617)里介紹了一種合理的計算方式,但顯然這不是唯一方式。再比如分子半徑,我在《談談分子半徑的計算和分子形狀的描述》(http://www.shanxitv.org/190)就已經列舉了多種計算方式。還有諸如原子電荷,計算方法更是多了去了,如ADCH、RESP、Hirshfeld、NPA等,簡要介紹見《原子電荷計算方法的對比》(http://www.whxb.pku.edu.cn/CN/abstract/abstract27818.shtml)。
為了正確理解文獻具體算的是什么,必須把背景知識搞清楚,典型的就是gap這個概念,有HOMO-LUMO gap、band gap、optical gap、fundamental gap等,參見《正確地認識分子的能隙(gap)、HOMO和LUMO》(http://www.shanxitv.org/543)。如果沒搞清楚概念和定義方式,可能重復不出文獻的數據,或者實際上文獻的計算方式或用詞不當,你卻意識不到這一點。
某些問題的計算中可能涉及到一些實驗值,要搞清楚作者用的什么值。比如用量子化學程序計算含碳物質的生成焓,通過構造熱力學循環可知其中需要利用實驗的(或第一性原理計算的)石墨的升華焓,要搞清楚這個量是從哪來的,具體是多少。再比如, 量子化學計算某分子的標準氧化還原勢,其中利用到標準氫電極的電勢,其數值有不同來源,從4.05到4.44 V都有報道,Phys. Chem. Chem. Phys., 16, 15068 (2014)建議對理論計算的情況用4.28 V。
還有一些研究涉及到一些額外的參數需要注意,典型的就是頻率校正因子,見《談談諧振頻率校正因子》(http://www.shanxitv.org/221)。是否使用頻率校正因子,以及對某一特定級別使用什么目的的、哪個來源的校正因子,都會影響結果。
周期性體系的計算中,軌道能量的絕對數值是沒意義的,不同程序用的能級零點也不一樣,需要關心的只有相對數值,比如帶隙。所以不要問為什么DOS圖的形狀和文獻一致但橫坐標和文獻不一致。這也正是為什么文獻在繪制DOS圖時往往一般把價帶頂或者費米能級位置設為0點,因為這樣才便于公平對比。
還有些量在文獻里用的習俗和你用的可能不同。比如算結合能問題,有少數文章把結合過程的能量降低量作為正值記錄。再比如NBO給出的E2超共軛作用能,程序輸出時把正值作為超共軛作用造成的體系能量降低,但有的文獻報道E2值的時候取的是負值。
如以上列舉的情況可見,算某個東西的時候必須仔細閱讀文獻搞清楚作者是怎么算的、用的什么定義、算的到底是什么。
5 數據的可重現性本身就差
有些數據很容易重現,比如在RHF/6-31G**級別下優化水分子的基態結構,不管用什么主流量化程序、怎么優化,得到的都是相同結果,初始結構很隨意。有些數據本身就是難以精確重現,典型的就是凝聚相體系的分子動力學。由于分子動力學過程有混沌效應,而且實際計算中往往不可避免地引入隨機因素,會導致文獻中的動力學模擬的定量結果幾乎無法嚴格精確重現。關于這點我專門寫過文章,參看《數值誤差對計算化學結果重現性的影響》(http://www.shanxitv.org/88)。更具體來說也看重現文獻里的什么動力學模擬數據。比如小分子的液態性質相對來說是容易重復的。比如重復密度值,文獻里是801.92 g/L,你算出來801.88 g/L,這就可以認為是很好重復出來了。如果你算出來是810.21 g/L,而且模擬流程和設置是合理的,那就是有明顯不可忽視的差異了,已經超過了隨機性的影響了,需要根據前面所說的來考慮導致與文獻存在差異的可能原因。也有的量的重現相對比較困難,比如磷脂膜里插入了一個分子,去計算它的側向擴散系數,這就屬于相對比較難算準的,因為本身分子在磷脂膜這種擁擠的環境里擴散就困難,而且體系里就這么一個分子,沒法通過同類分子的平均來減小統計誤差。就算你自己跑的時間非常長以保證統計精度絕對足夠高,但文獻里給出的數據未必統計精度就夠高。像這種問題,對文獻里的數據的重現精度就別要求太高了。
還有一種情況是動力學模擬的現象在有必然性的同時本身也有一定隨機性,跑出什么現象有運氣成分。比如距離蛋白質的口袋一定距離處放一個小分子配體,在水環境下做動力學模擬,考察2 ns模擬時間內配體能否和蛋白質結合,初速度隨機生成。像這種模擬的結果就很有隨機性,如果跑10次模擬,可能有的模擬剛過幾百ps小分子就進入口袋了,也可能有的模擬開始后小分子恰好往蛋白質相反方向游離,跑到2 ns還沒進入口袋。對于這種問題的研究必須做大量的模擬來得到統計結果,較嚴謹的文獻普遍都會這么考察這種現象明顯有隨機性的問題,你若要重復文獻也必須以相同的方式去模擬。如果模擬的樣本數不夠多的話,也不要指望結果特別相符。比如前述問題,可能文獻里跑10次發現有5次進入口袋,你模擬10次發現有3次進入口袋,這未必是模擬設置和文獻有區別,而純粹是隨機性所致,要以正確的態度來看待這種情況。
順帶一提,量子化學、第一性原理做靜態的計算的可重復性比起分子動力學模擬要好得多得多,但也不排除極個別情況前者可能受到隨機性的明顯影響。比如對勢能面很復雜的大分子進行優化,特別是QM/MM計算時可能牽扯到幾千原子的情況的優化,有可能由于微小的隨機性因素(如使用了并行計算所帶來的)導致兩次優化分別收斂到了勢能面上可能相距挺遠的兩個極小點,且有肉眼可察覺的差異。
6 對相符程度要求過高
有人其實可能已經合理重復出了文獻中的數據,但由于他對重復的精度要求過高,導致誤以為沒重復出來。
頭腦里要有收斂限的概念。比如文獻里優化的結構二面角是123.64度,你優化出來的是123.69,這就可以認為是重現出來了,因為這點差異通常都已經小于計算程序的幾何優化的默認收斂限了,差異基本來自初始結構的任意性。再比如,文獻里DFT計算出的反應能是-8.12 kJ/mol,你算出來的是-8.18 kJ/mol,也完全可以接受了。如果本來你用的計算程序和文獻里都不同,在盡可能令設置和文獻相同的情況下,你算出來比如是-8.32 kJ/mol也可以算重復出了文獻,各種雞毛蒜皮不值得一提的因素足矣帶來這種程度的差異。但如果你算出來是比如-9.8 kJ/mol,那就明顯不是不值得一提的差異了。再比如TDDFT計算電子激發能、DFT計算HOMO能級,和文獻里相差0.02 eV,這都可以算一致,但如果相差到0.1 eV,這就不算一致了,畢竟TDDFT算激發能的精度都在0.3 eV左右(泛函選擇得當的前提下),當差異達到接近方法本身的誤差的數量級時明顯就算顯著差異了。
重復文獻數據的時候用的收斂限應當足夠嚴,起碼令收斂限顯著小于你期望的對文獻數據的重現精度。雖然文獻作者用的收斂限未必夠嚴,因此光是你把收斂限設嚴未必有意義,但至少先做自己這里能做的事。
7 自己的計算存在硬傷
重復不出文獻數據有可能是自己計算存在硬傷導致的。初學者缺乏常識、經驗和敏銳性,計算很容易存在硬傷,例如:
? 關鍵詞沒寫對或設置不當,導致數據沒意義,比如:Gaussian里想要用PBE0泛函但寫的關鍵詞是PBE0(此時做的是PBE0-DH泛函計算)、用贗勢基組時只定義了基組部分而沒定義贗勢、自定義基組時漏定義了某些原子或元素、Gaussian里用IOp自定義泛函時用opt freq關鍵詞(IOp沒法自動傳遞到freq這一步導致振動分析不是在與優化相同級別勢能面的極小點計算的)、Gaussian里用后HF方法算偶極矩時沒寫density結果得到的是HF級別的、計算UV-Vis光譜時算的態數太少導致光譜不夠完整、CP2K計算時用雜化泛函并且基于密度矩陣做積分屏蔽時沒讀取純泛函波函數當初猜(《解決CP2K中SCF不收斂的方法》http://www.shanxitv.org/665文中提到了)、CP2K算磁性體系時沒寫UKS結果當成了閉殼層算、Gaussian里按照特定方向加電場計算時不知道要寫恰當用nosymm導致電場加的方向和期望的不符(nosymm用處見《談談Gaussian中的對稱性與nosymm關鍵詞的使用》http://www.shanxitv.org/297)、用sobtop(http://www.shanxitv.org/soft/Sobtop)直接產生拓撲文件后沒有給里面添加原子電荷、用sobtop產生拓撲文件時誤把可旋轉二面角以諧振勢對待、用Multiwfn做空穴-電子分析時沒按照《使用Multiwfn做空穴-電子分析全面考察電子激發特征》(http://www.shanxitv.org/434)說的在Gaussian輸入文件里寫IOp(9/40=4),等等
? 缺乏基本知識瞎算,比如:計算小晶胞體系不知道考慮k點、算極化率時不知道要帶彌散函數、B3LYP算色散主導的弱相互作用卻不加色散校正、CP2K里用MOLOPT基組計算需要CUTOFF很高的Na的時候用的CUTOFF不夠、界面體系用了各向同性控壓而不知道應該用半各向同性控壓、算氣-液界面體系在垂直于界面方向開了控壓結果導致真空區消失、用Multiwfn基于CP2K產生的molden文件分析前沒在里面添加盒子信息(參見《詳談使用CP2K產生給Multiwfn用的molden格式的波函數文件》http://www.shanxitv.org/651)、做Mulliken布居分析用的基組帶了不該帶的彌散函數、算溶劑環境下溶質自由能的時候忘了考慮標準態濃度轉換問題(參看《談談隱式溶劑模型下溶解自由能和體系自由能的計算》http://www.shanxitv.org/327)、試圖用后HF方法計算分子軌道(卻渾然不知程序雖然給出了結果但實際上是HF部分產生的)、用背景電荷表現環境分子的勢場用的是對靜電勢重現性極爛的QEq電荷(相關討論看《基于背景電荷計算分子在晶體環境中的吸收光譜http://www.shanxitv.org/579),等等
? 計算流程/方式不對,比如:算高精度能量和性質之前結構沒經過幾何優化、直接把IRC兩端的結構當做反應物和產物的準確結構用于計算能量和屬性數據、振動分析和IRC計算之前結構沒有在嚴格相同級別下優化過、對動力學軌跡計算某個分子的RMSD衡量內部結構變化的時候不知道先要消除整體平動和轉動、算電子發射光譜之前沒優化激發態、模擬拉曼光譜時沒按照《使用Multiwfn計算特定方向的UV-Vis吸收光譜》(http://www.shanxitv.org/648)里說的將拉曼活性轉化成拉曼強度
? 數據讀取錯誤。比如Gaussian里做了雙雜化泛函計算,結果讀成了其中SCF Done后面的中間量(怎么正確讀看《談談該從Gaussian輸出文件中的什么地方讀電子能量》http://www.shanxitv.org/488)。再比如做幾何優化,Gaussian會對初始和最終結構都做布居分析,本來想讀最終結構下的,結果誤讀成第一次輸出的。再比如做開殼層體系的NBO分析時,NBO會對總密度矩陣、alpha密度矩陣和beta密度矩陣的分析結果依次輸出,如果讀錯地方就糟糕了
? 對計算數據質量缺乏把控。如果數據質量太糟糕,自然可能無法重現文獻數據。各個部分的計算都要保證數據質量。比如算穩定狀態必須確保波函數穩定、確保無虛頻。雖然做波函數穩定性測試和振動分析會增加耗時,但必要時必須做,否則可能得到錯誤的結果而無法重復文獻。不要輕易使用比如Gaussian里opt=loose這種關鍵詞放寬收斂限,否則很可能收斂限太松導致某些情況對文獻數據重復性太差,尤其是此時做振動分析得到的頻率的準確度很可能很糟糕。用取消SCF收斂情況檢查的IOp(5/13=1)這種危險關鍵詞的時候必須謹慎檢查最終SCF收斂情況,千萬別最后收斂精度很爛卻直接使用其數據。再比如對于通過分子動力學計算物質平衡狀態屬性時,必須保證采樣充分使統計誤差足夠小。再比如考察某種磷脂膜體系的特征時,需要模擬時間足夠長以使得體系充分達到平衡,本身這類體系達到平衡所需的時間就偏長,動輒需要100 ns以上
? 沒恰當考慮對稱性問題。對于有對稱性的體系,在Gaussian等支持對稱性的量子化學程序里,要盡可能利用對稱性,這樣計算又好又快,內行的研究者的文章里對于能利用對稱性的情況幾乎都會利用。不利用對稱性甚至可能還會得到和文獻不符的結論,比如文獻報道某個體系是平面的,但你一開始擺的結構是彎曲的,由于優化的收斂不可能無窮精確,最后你得到的結構可能輕微偏離平面,與文獻結論不一致。如果體系本身沒有那么高對稱性,你一開始當成高對稱性結構來計算,最后可能錯誤地得到了高對稱性結構,也和文獻不符(這種情況只要做一下振動分析看有沒有虛頻便知)。
? 其它低級錯誤:單位轉換因子沒用對或者用的不準確(比如百毒搜出來的大量轉換因子都是錯的)、該做構型/構象搜索的情況沒做或者做得不當、該考慮構型/構象權重平均的情況沒考慮、優化表面吸附問題時把被吸附分子放得太遠導致還沒優化到吸附狀態就滿足了收斂判據、處理數據用的Excel表格或者腳本存在問題、科學計數法記錄的數據沒讀對(如漏掉了指數部分),等等
8 文獻自身有問題
千萬不要只懷疑自己而不懷疑文獻,盡信文獻不如無文獻,很多文獻里(包括IF很高的期刊的文章里)的數據是有錯誤的,甚至是非常明顯的錯誤,上一節提到的各種犯錯的例子也都出現在很多文獻里。比如《我對一篇存在大量錯誤的J.Mol.Model.期刊上的18碳環研究文章的comment》(http://www.shanxitv.org/584)里我就指出一篇文章里從頭到尾的一大堆錯誤,諸如作者計算C18的時候大概率由于用的是非平面的初始結構,導致優化出的C18不是精確平面的,然后他們發現C18的極小點結構竟然有ECD信號,明顯這是虛假的結果。
還有些文獻里對數據的解釋、標注是錯的,是他們對數據理解有問題。比如直接把HOMO、LUMO能量的負值說成是氧化、還原勢是很多涉及電化學的文章里常見的事,顯然你如果以正確的方式來算,也即計算溶液下的電極反應的自由能變,結果肯定和文獻對不上。再比如有不少文獻用的費米能級是錯的,對有gap的體系直接就把價帶頂當做費米能級,這明顯不符合正確的確定有限溫度下費米能級的方式(正確定義參看wiki相應條目和Multiwfn手冊3.300.9節)。
還有些文獻里的計算方式從描述上就知道是錯的。比如有的文獻寫G=Eele+ZPE-TS,有本科化學知識的人都知道這明擺著不對,H(T)和U(0)之間的差異上哪去了?就算是作為近似不考慮這部分、假設這部分在求差的時候能極大抵消,至少也應該把等號寫成約等號。這種明顯不合理的數據就別去重復了,非要重復的話那就暫時性降智故意用錯誤的公式。
很多文獻作者也對數據很不負責任。比如之前我在IF很高的文獻里看到補充材料中的Gaussian輸入文件里居然全帶著IOp(5/13=1),這文獻的數據的可靠性沒法令人不放心(雖然作者可能檢查了最終收斂情況,但誰也不敢說作者確實這么做了)。現在有很多無良代算機構,沒一丁點搞理論計算的人的基本操守,為了能給出令做實驗的人滿意的與實驗相符的計算數據來交差,極有可能在算不出期望的數據的情況下捏造假數據,這樣的數據更別指望能被別人重復出來了。甚至有的代算方連計算都不計算,直接憑空胡編亂造數據,建議讀者看看《談談我對計算化學代算的看法》(http://www.shanxitv.org/505)里提及的情況。
有很多文獻對計算的關鍵要點、直接影響數據可重復性的因素交代得非常粗糙、簡略、語焉不詳,甚至根本不提,這無疑給讀者帶來了很多誤解、給別人重復計算造成了障礙。重復這些文章的數據只能是考慮各種可能情況去試了。另外,還有可能文章里用的參數、設置本來就不慎寫錯了,作者寫的程序本身存在bug等等。
對于已經恰當考慮了前面說的所有因素,竭盡全力去重復文獻,但還是重復不出來(且確實是有不可忽視的差異而不是在復現精度上過度較真),如果沒有絕對必要去重復文獻,那就別去重復了,極有可能文獻的數據就是有問題的,或者沒交代關鍵性計算細節。這種情況只要你能保證自己的計算是合理的就夠了。但如果有特殊原因就是非要重復文獻不可,顯然該做的不是在思想家公社QQ群或者計算化學公社論壇(bbs.keinsci.com)等公開場所提問,因為找誰都不如直接找作者,世界上還有誰比作者更清楚數據是怎么來的的人么?世界上還有誰比作者更有回復你的義務?發郵件給通訊作者,向作者索要相應數據的輸入輸出文件,或者向作者提供你復現文獻計算的輸入輸出文件問作者為什么沒重復出來即可。當然很有可能一直沒收到作者的回復,可能性有這些:
(1)作者看漏了你的郵件,或者郵件被誤歸入垃圾郵件,或者正忙著什么事沒來及回復你結果后來又忘了。可以過一個禮拜(且可以換個郵箱)再發一次試試。如果通訊作者有多個,也可以給另外的人發
(2)措辭不當,不夠禮貌和誠懇,而作者脾氣又大或者又比較忙,又覺得回復你也沒什么好處,就不回了。可以過一陣重新措辭再發一次。為了讓作者有動力花時間給你回復,最好表示你會在未來的文章中會引用他的文獻。
(3)文獻中的數據是作者找別人代算的,代算方又是非常不負責任的,相關文件也不提供、計算細節也不交代清楚,甚至代算之后就再也聯系不上了,因而作者也無能為力。
(4)做計算的是學生,學生畢業后找不到人了,導師又不清楚計算細節、沒原始文件。
(5)文獻中的計算本身就有嚴重錯誤或是假數據,可能收到你的郵件后作者才意識到,又或者可能之前就已經明知道用的是假數據,總之都不希望你知道他的數據存在問題,免得面子上掛不住,甚至怕文章最后被撤稿。