數值誤差對計算化學結果重現性的影響
1 前言
曾經有很多人問，諸如為什么明明輸入文件一樣，但是動力學軌跡每次跑出來的都不一樣、為什么幾何優化結果相差甚遠等等，甚至懷疑計算科學是否遵循決定論。實際上結果的差異是由于運算開始或運算過程中各種形式的數值誤差、隨機性引起的。這個是很重要卻不被多數人重視的問題，它和理論本身、算法、軟件環境、硬件環境等都有密切關聯。本文將對這個問題做一些分析探討。
首先先討論數值算法、軟件、硬件因素是如何導致誤差（錯誤）和隨機性的。然后再看這些問題對計算化學會產生何等的影響。
2 數值誤差（錯誤）、隨機性的根源
2.1 內存
因為不少人對內存有誤解，所以下面說得多一些。內存在數據讀、寫的過程中不可避免地會發生一些錯誤，輕則數據出現異常（比如一個變量為1卻成了0），重則死機、重啟。和內存錯誤（包括各種類型錯誤）相關的主要有以下幾個方面：
(1)供電質量差。如電壓偏離標準值過多、電流不穩。這不僅取決于電源質量，還取決于主板好壞（內存供電電路）。
(2)過熱。隨著溫度升高內存出錯的幾率總是增加的，但只有高溫下增加的幅度才最明顯。遙想筆者用i860+RDRAM的日子，不裝個小風扇吹著內存的話高負載下20分鐘左右就死機。燙手的話應該裝上散熱片乃至小風扇。
(3)內存質量。一方面是內存設計，即布線、PCB板數、貼片元件等。另一方面取決于內存顆粒制造商的品質，也有RP的關系，比如顆粒來自晶圓的邊角還是中心（雜質含量不同）。
(4)制程。一般來說制程越精細可靠性越低，比如22nm制程的比起0.35微米制程的更容易出問題。但不要忽視工藝質量也在逐步改進，新的精細制程的產品出錯幾率不一定比以往的粗糙制程的要大。
(5)宇宙輻射。這看起來比較微妙，但是當嚴重的時候，如太陽風暴劇烈時，對電子設備的影響就明顯了。這也是為什么在一些宇航領域（接受的輻射更嚴重）的芯片仍然用著設計過時、制程古舊的產品，沒準兒出現一個小錯大家就沒命了。機箱內電磁輻射也對內存穩定性有影響。
(6)內存參數。所設的工作頻率越高、延遲越低出問題的幾率越高，當然也要看內存類型。
(7)內存容量、讀寫次數：內存容量越大、讀寫越頻繁，相同時間內出錯的幾率越大。
(8)老化。在越惡劣的環境使用越長時間，出錯的可能性越大。
要注意以上的因素不是分立、絕對的，而是互相關聯的、相對的，不能單獨拿出來對比。
可見，與內存數據出錯相關的因素很多，也很難衡量內存平均多長時間出一次錯。一個令我比較認可的、由大規模統計得到的說法是一般的使用環境下（外界環境正常，計算機運行穩定），平均每年每條內存有3%幾率出錯（包括各種形式的錯誤）。實際上這個幾率是很小的，換句話說，等到機子被淘汰了也未必能趕上一次內存出錯。所以不要把內存問題過分夸大，比如持續幾天的計算，計算結果無法重現由內存錯誤導致的可能性微乎其微。但是，對于數據中心，內存出錯幾率乘上相應內存數目，再考慮數據出錯的損失，這就是大問題了。內存的Error Checking and Correcting (ECC)顯得十分必要，這是一種數據錯誤檢測和校正的技術，每次在寫入數據時同時寫入額外的數據位用作校驗，在讀出數據時也同時把校驗位讀出，由此根據一定算法分析讀出的數據是否和寫入的數據相同，若不同則進行修復。這使得內存數據出錯幾率大大降低，但也并非能解決所有問題。有些人誤認為ECC內存性能比普通同類型內存要低很多，其實相差甚微。現今服務器內存幾乎都帶ECC，它比普通內存略貴，多了額外內存顆粒用于儲存校驗位，而現今服務器主板也都支持ECC（有些主板則只支持ECC內存）。個人計算沒有絕對的必要使用ECC內存。
導致內存出錯的原因對于CPU/GPU運算出錯也基本是一樣的，就不累述了。
2.2 處理器類型
對于不同架構的CPU、擁有不同指令集的CPU、乃至GPU，不能期望它們能完全重現同樣的結果。尤其是浮點運算的標準可能不同，比如末位的舍入、極小值的處理。有些指令在不同處理器中執行方式不同，比如先乘一個數再加上一個數，這兩條指令在有的處理器中會被融合為結果更精確的“乘加”一條指令。
也有一些很特殊的情況，數值錯誤是由CPU設計bug引起的。老玩家們應該還記得，較早一批的Pentium 60/66Mhz的浮點除法指令有錯誤。
2.3 網絡數據傳輸
對于集群，尤其是分布式計算，需要考慮網絡連接引入的錯誤。網絡傳輸都會有一定誤碼率，即在多少位數據中出現一位差錯，其大小取決于網絡規范、傳輸介質、傳輸距離、外界干擾、網絡適配器質量等因素。主流的以太網、infiniband等網絡標準本身都帶校驗幀以保證數據完整、準確性，已很大程度上減輕了網絡傳輸引發的數據錯誤。
2.4 隨機數生成器
在一些程序中要用到隨機數生成器，比如Monte Carlo過程、動力學的初始速度設定。通常隨機數是以不同的隨機數算法產生的，常用的比如乘同余法。這些算法都需要一個稱為“種子”的數，不同的種子將產生出不同的隨機數序列。也就是說，如果程序中利用了隨機數生成器，想重現結果就必須用相同的種子。但種子往往不允許由用戶設定，也不是固定不變的，而是根據一些規則產生來保證每次運行結果的隨機化，比如利用當前系統的時間，對這種情況用戶不要指望結果會重現。
2.5 庫文件、操作系統、程序版本
不同的數學庫中相同的功能在內部會使用不同的算法或者不同的編寫方式，即便是同一個庫，在哪怕相近版本中，也可能由于對代碼進行的優化、除蟲，導致不同結果。有些庫函數的運行結果甚至取決于運行期的情況，比如快速傅里葉變換庫FFTW會在啟動時檢測哪種算法在當前系統下效率最高，然后會在之后一直使用，這就可能在檢測的時候由于正在運行的其它任務對系統CPU資源占用不同，使檢測結論不同，影響后續算法。
有些程序是通過靜態鏈接產生的，庫文件會被合并在程序可執行文件當中，或者有些程序雖然使用動態鏈接，但是自帶了動態庫文件，這樣不同系統下結果的重現性還好說。而有些程序需要依賴于當前系統中的庫文件，不同系統自帶版本不同，差異就容易出現了。若平臺不同，如一個程序的Linux版本與Windows版本，差異會更明顯。
不要指望不同版本號的程序能獲得相同結果，不僅所用的庫往往不同，默認參數、算法也經常不同，比如Gaussian03的幾何優化默認用的berny方法，而09版后來默認用GEDIIS，同一個輸入文件有時前者能優化成功而后者卻不收斂。即便算法和參數都不變，代碼細節也可能改變。
2.6 編譯器與編譯選項
編譯器類型及具體版本號、編譯參數對結果影響明顯。不同編譯器會對代碼進行不同形式的優化，有些優化是以犧牲精度為代價的。如果調用了程序語言標準中定義的標準數學函數，比如開方、求三角函數、求對數等，不同編譯器會用不同的內建數學庫，顯然也會帶來差異。不同編譯器對于語言標準中沒有定義的情況的處理也不同，比如未初始化的變量，有些編譯器一律將之初始化為0，而有些編譯器則不這么做，相應內存地址位置目前是什么數值就是什么數值，由于這個數值無法預測，每次運行時可能造成結果不同。
來看一個簡單的比較，用Multiwfn 2.02對一個水分子全空間積分電子密度拉普拉斯值的結果。Intel visual fortran與CVF結果有所不同，在不同優化參數下結果也有所不同。
IVF12 0d/O1 -2.603320357104659E-005
IVF12 O3 -2.603320356662977E-005
CVF6.5 debug/release -2.603320356424300E-005
2.7 并行執行時的運算順序
雖然諸如a+b+c=a+(b+c)在數學上是成立的，但是在計算機執行中未必成立。比如在CVF6.5里執行這樣的語句
a=sin(4D0)
b=sin(99D0)
c=sin(2D0)
write(*,*) a+b+c-(a+(b+c))
得到的結果不是0.000000000000000E+000，而是1.110223024625157E-016，如果結果被乘上一個比較大的數，那么這個差異就不能被忽略了，可見計算機運算未必滿足結合律。在并行運算中，經常要把一部分工作拆成幾個線程以分配給多個處理器并行執行，然后將算出來的結果累加到一起。通常哪個線程先算完，哪個結果就先被累加上。然而實際計算中各個線程完成的先后順序是不確定的，每次程序執行時都可能不同，因此數據累加的順序不一樣，這就造成結果出現一定不確定性。另外，有的程序中使用動態均衡負載以獲得更高的效率，這導致每個線程執行的內容也是有隨機性的，使結果的重現更加困難。
舉個實例，在Gaussian09中用4核計算一個體系單點能，執行三次，最后一次迭代的能量分別為
-234.579551993147
-234.579551993143
-234.579551993144
結果在末位有所差異。如果只用nproc=1來串行計算，則結果每次都能重現。盡管并行時每次的結果差異看似微不足道，以a.u.為單位一般精確到小數后5、6位足矣，然而這個差異在其它一些任務中會被逐漸放大，甚至可能達到影響定性結論的地步。有些人以為通過增加積分精度、收斂精度之類能減輕結果的隨機性，實際上從原理上講根本于事無補，根本不是一碼事（僅在極個別情況下能有改善）。
2.8 文件格式
不同文件格式的精度，以及格式轉換過程中的信息丟失是需要注意的。浮點數據在計算機中通過二進制形式記錄，而在磁盤上記錄的文件往往是文本格式，這便于人類閱讀和第三方程序處理。在內存中占8字節的雙精度變量若想完整以文本方式記錄下來起碼要用22字節，這不僅太占空間，閱讀也不方便，因此往往故意減少有效數字位數。比如Gaussian的chk文件是二進制文件，浮點變量以二進制形式精確保存，有效數字位數約16。而轉換成fch后，浮點數有效數字就被截為9位，若用unfchk再將之轉化回chk，則此時的chk的數據精度比起之前的chk已經損失了。即便是二進制文件也并非都是完整精度形式記錄數據，比如Gromacs的xtc文件也是二進制形式，而變量記錄精度可控（xtc_precision參數），默認的記錄精度低于單精度浮點變量（這也是為什么Gromacs軌跡文件比較小的原因），因此不要指望rerun的過程，即重新算軌跡中每一幀的勢能的值會與之前動力學過程中輸出的一致。有些程序在導出/導入的過程中還可能會對結構、導數等等進行坐標變換，對單位進行轉換等操作，也會引入數值誤差。
2.9 操作過程
一些操作過程的細節差異往往不被注意，對最終結果的影響可能卻不小。例如對晶體測得的結構加氫，很多程序加氫結果表面看似都一樣，位置沒有區別，但是查看文件中記錄的坐標，卻總有細微差別。再比如模擬一個NMR測得的體系，pdb文件通常會有很多幀，往往每一幀差異很小，似乎取哪幀是隨意的，但實際上用于模擬會得到不同結果。還比如Gaussian中設定不同內存大小，看似和計算結果無關，但實際上Gaussian會自動根據分配的內存容量選擇不同的積分計算、儲存方式，也會引入差異。
總結一下，上面的討論中，數值差異可以歸納為兩部分，一種是可重現性差異，也可以叫靜態差異，包括硬件配置、編譯過程、庫文件、文件操作、操作過程、程序版本，只要運行條件完全一致就可以徹底避免；另一種是不可重現性差異，或者叫動態差異，它取決于運行期的實際條件，有的可以通過一定對策避免，比如把并行運算改為串行，但有的很難完全消除，尤其是硬件偶發錯誤。
3 數值誤差、隨機性對計算化學結果的影響
從前面簡單例子看到，誤差或者隨機性對結果的影響甚微，在默認的輸出精度下甚至看不出任何差異。然而對于某些任務情況卻不是如此，數值的微小擾動（下面用擾動這個詞代表前述各種因素造成的微小數值差異）可能造成結果有明顯定量乃至定性的差異，這是因為在這些任務中擾動的效應會被放大。放大的方式主要可以歸結為兩類，一類是數值運算方式，另一類是體系和任務的混沌效應。
3.1 數值運算方式對擾動的放大
一些任務中，由于算法的原因，初始的擾動會在結果中被放大，典型的例子是數值差分。在幾何優化、頻率計算、超極化率計算等任務中，若理論方法在當前程序中不支持能量解析導數計算，就必須要通過數值差分來近似計算導數。x處一階導數計算方法為：
d1(x)=( E(x+Δ)-E(x-Δ) )/(2Δ)
其中Δ是差分步長。假設最糟的情況，計算E(x+Δ)和E(x-Δ)時都引入了大小為k的擾動，且符號相反而無法被抵消，則此時的一階導數寫為
d1'(x)=( E(x+Δ)-E(x-Δ)+2k )/(2Δ)
那么由于擾動，造成一階導數誤差為d1'(x)-d1(x)=k/Δ。由于Δ必須足夠小以保證導數計算準確，比如為0.0001，那么初始擾動k就會被放大10000倍，對結果造成了10000k的影響！
再來看差分計算二階導數
d2(x)=( d1(x+Δ)-d1(x-Δ) )/(2Δ)
由于計算d1的時候已經引入了k/Δ的誤差，則實際計算出的二階導數為
d2'(x)=( d1(x+Δ)-d1(x-Δ)+2k/Δ )/(2Δ)
二階導數的計算誤差為d2'(x)-d2(x)=k/Δ^2，若Δ=0.0001，則初始擾動k就會被放大100000000倍！可見，差分計算n階導數時由于初始擾動值k最多可能導致結果產生k/Δ^n的誤差。可想而知，計算三階導數和四階導數時，比如有限場方法純數值計算一階和二階超極化率時，誤差已經被放大到不可忽略的地步了。所以，Δ取得越小，由于誤差會越大，結果未必越精確，必須選擇最合適的值。
通過選擇合理的算法、良好的編程習慣，可以使一些任務中誤差的放大盡可能小，比如避免大數除小數之類的情況。這就是編程者的事情了。
3.2 體系、任務的混沌性對擾動的放大
混沌效應對擾動的放大是特別需要關注的。混沌效應粗俗地說就是捷克一個人從體內釋放一股不圣潔的氣體可能造成韓國一場臺風，用古話來講則是“失之毫厘，謬以千里”。混沌效應對數值計算巨大影響的發現，可追溯到1959年洛倫茲用計算機進行氣象模擬（其實更早之前有人在計算天體運行軌道時也多多少少意識到了這點）。他原先模擬時一直使用十幾位精度的浮點數保存數據，后來有一次為了更細致地分析某段時間的數據，他用以前的數據為初值重算了一遍，為了打印的數據美觀，這次他將小數精度只保留三位。正是這微小誤差的引入，導致模擬出的兩個月之后的氣象結果和原先大相徑庭！在計算化學中，有兩類任務混沌效應很明顯，一個是分子動力學，另一個是幾何優化。
3.3 擾動對分子動力學的影響
分子動力學的混沌效應，尤其是對蛋白質折疊過程的影響被研究得比較多，如PROTEINS,29,417。動力學過程中初始擾動會隨著模擬時間指數型放大，這并不是由于數值計算問題所導致的，而是體系內在性質的體現。若在t=0時對體系坐標引入Δ微擾，則模擬到t時軌跡與未引入微擾時的差異可寫為diff(t)=Δ*exp(λt)，其中λ是依賴于體系的Lyapunov指數。擾動的放大速度相當快，研究發現對蛋白質結構引入RMSD=1*10^-9埃的擾動，經過2ps的模擬就能被放大到RMSD=1埃。假設在模擬中觀察到某時刻有一個氨基酸側鏈發生了扭轉，如果在此10ps前引入10次RMSD=0.001的不同方式的擾動，并分別跑10次動力學，則只有一次能夠重現側鏈扭轉。可見，動力學過程對初始條件十分敏感，各種因素帶來的微小數值擾動都足矣使動力學軌跡發生很大分歧。更何況實際模擬中擾動并不僅僅是在初始條件中出現，在模擬過程中還可能會不斷納入。
注意對于蛋白質體系，擾動不是隨模擬的進行而無限放大的，因為蛋白的原子受到力場的束縛，且模擬時的溫度有限，因此勢能面可及的范圍是有限的，故初始擾動造成的差異在迅速放大后會遇到上限。對于柔性體系，由于勢能面可及范圍會加更大，可以預見差異放大的上限會比剛性體系更高。
一些辦法可以提高軌跡重現性。提高浮點數精度是重要策略之一，這可以減少擾動量的數量級。有人抱怨gromacs軌跡重現性不如Amber，這很大程度上是由于gromacs默認是單精度編譯，而Amber默認用雙精度。不過，提高精度并不解決根本問題，只不過是讓軌跡的歧化慢一些罷了，對于長時間的動力學模擬結果的偏差仍然是無法預測的，除非能讓浮點精度無窮高，但這顯然是不可能的。若想讓軌跡完全重現，通用的辦法是將第一節涉及的所有可能帶來擾動的因素都去除。而有些程序也專門考慮了軌跡重現問題，在gromacs里的mdrun有-reprod選項以幫助軌跡重現。另外還有人弄了個固定精度的分子動力學方法，不僅可以完全重現軌跡，還可以將時間反轉從末態得到初始結構，完全體現出分子動力學的決定性特點。
動力學軌跡的重現究竟有無意義？如果只是出于計算目的，比如一段軌跡文件弄丟了，或者希望以更高的頻率輸出體系能量，那么重新跑出一次完全一樣的軌跡是有意義的。然而，從理論角度上來講，軌跡的嚴格的重現性是毫無意義的。動力學過程就是在體系相空間（坐標+動量空間）中游歷的過程，如果跑兩次得到了截然不同的軌跡，只不過是說明在相空間不同位置區域游歷罷了，更沒有誰對誰錯之分。
對于一個平衡態的模擬，體系的性質（比如內能、擴散系數）就是動力學過程時間的平均，兩次跑出來軌跡相不相同完全無所謂，只要時間足夠長，它們都足以充分遍歷相空間，結果會是等價的。對于有限時間的模擬，為了加強某些區域的采樣，實際上還有人專門利用動力學軌跡的混沌效應，給初始速度分配不同的值，分別做動力學模擬并將結果合并。
而對于非平衡過程的模擬，比如說蛋白質折疊過程，軌跡（折疊路徑）的重現性同樣毫無意義。因為蛋白的折疊過程本來就是混沌的，初始的微小擾動會使折疊路徑相差迥異。在現實中每次蛋白的初始條件絕不可能完全一致，因此同樣的折疊路徑在現實中無法重現，或者說現實中沒有可能發生兩次完全相同的折疊過程，那么苛求在多次模擬中重現相同的軌跡有何意義？實際上不僅不應該重現，折疊路徑的多樣性才是真正感興趣的，也很適合通過概率統計來研究，這就需要做很多次動力學模擬，以得到不同的折疊路徑并進行分析。正好，由于各種數值因素的介入，哪怕相同的初始條件都會產生不同的折疊路徑，都省得修改輸入文件了，直接反復跑多次就行了。注意盡管每次的軌跡十分不同，但最終卻一定會折疊到同樣的最穩定結構，這并不違背混沌效應，而是由于每一次折疊過程都會感受到相同的“吸引子”，即趨向成為最穩結構的一種特殊勢場。折疊過程可以比喻為站在一個坑坑洼洼的大坑（對應自由能面的形狀）旁邊往坑里扔石子，每次石子掉落的軌跡都很不同，但最終都掉到坑底。
這里順便談談一個相關問題，也就是單個動力學模擬的軌跡是否足以說明問題。上面說了，一條折疊軌跡不能反映全部的折疊過程，同樣地，研究其它過程，原則上也需要用多條軌跡來統計地說明問題。比如研究在結合配體后激酶的活性位點打開的過程，在軌跡中發現經過530ps模擬后活性位點打開了，那么就真的能說明在現實中配體結合后經過約530ps后活性位點就能打開？非也，模擬10次，會得到10個不同的結果，恰好都出現在530ps是不可能的。經過10次模擬后，假設根據統計結果可說明位點打開的時間主要分布在300ps~1.0ns區間內，且最容易在500ps附近打開，則這樣的統計性結論比起從單個軌跡中下的結論有意義、可信得多。混沌效應導致動力學過程中各種“事件”發生的時刻的波動范圍極大，或者說每種事件在不同的時間范圍都有一定出現概率，只分析單一軌跡就造成了一葉障目，結論的不可靠性會較大。然而目前絕大多數模擬都是基于單一軌跡的，很少有人對同一問題做多軌跡統計分析，主要是多軌跡的計算量太大或者研究的內容不需要那么深入詳細，其次是很多人沒意識到這個問題，還有少數人是因為投機取巧。比如論證一個抑制劑的作用，可能模擬反復跑了7次都無法論證抑制劑管用，當跑到第8次時慶幸地發現能說明抑制劑管用，于是他只用這一次的軌跡說事，另外跑的7次雖然沒成功卻都避而不談，可想而知這抑制劑的實際功效如何了...
3.4 擾動對幾何優化的影響
相較分子動力學的混沌效應，幾何優化的混沌效應被關注的相對較少，有興趣者可參看J Comput Aided Mol Des,22,39，討論數值誤差對分子力場優化結果的影響。
之所以幾何優化也有明顯混沌效應，不妨將幾何優化看作是0K下的分子動力學，尤其對于最陡下降法優化，這個比喻最為恰當。最陡下降法優化過程類似把一個玻璃珠輕輕放到陡峭的山上某斜坡處，讓它自動滾落。雖然0K下沒有動能（準確來說尚有振動能），通常不會導致引入擾動后幾何優化過程的軌跡像動力學軌跡那樣隨步數增長分歧得那么快，但是，由于沒有了動能以越過勢壘，體系將沒機會感受到全局的吸引子（全局最穩結構），而只能感受到離初始位置最近的局部吸引子（勢能面局部極小點），這導致了幾何優化的最終結果的重現性可能比動力學更差。對于小體系還好，勢能面結構比較簡單、光滑，然而極小點、鞍點的數目是隨著原子數目增加呈指數型增長的，對于大體系，勢能面的復雜程度是難以想象的，微小的擾動可能就會使結構收斂到很不同的極小點。
之所以平時在量子化學的幾何優化中優化的結果重現性都很好，沒有顯示出太大不可重現性，主要是因為量化處理的體系一般都比較小，通常少于50個原子，所需優化步數一般也較少，通常少于100步。不僅誤差來不及充分放大，而且由于相對簡單的勢能面的形狀以及它起到的束縛作用，誤差也不容易放大。然而因數值誤差帶來的分歧偶爾仍可能顯著，記得有人在一臺機子上某個體系優化幾十步成功了，換了臺機子卻不收斂了。
重現性好壞與初猜坐標關系很大，在山腳下釋放一個珠子，和在山尖上釋放一個珠子，明顯珠子最終位置在后者的情況中受初始位置影響更大。所以優化時如果能給出一個好的初猜，比如在極小點附近的二次區域內（勢能面可以近似用二次函數描述的區域），那么優化幾乎一定收斂到那個極小點，數值誤差因素不會導致收斂到別處去。然而如果在很接近過渡態的位置作為初猜，一開始的坐標或者計算受力稍微有點擾動，很可能結構就會收斂到另一側極小點去了。重現性好壞也與體系特點緊密相關，同樣是720個原子，若說多肽的勢能面是坑坑洼洼的山丘，那么C720這個碳球的勢能面就是一個光滑的碗（至少是指勢能面中平衡結構附近的很大區域），顯然C720優化的結果幾乎不受初猜擾動的影響。
前面說到動力學軌跡的重現性沒有物理意義，那么優化結果的重現性是否有物理意義？或者，一次優化的結果是否足夠說明問題？對于結構簡單的體系，且初猜較好時，結果應當能重現，否則說明計算條件、軟件可能存在問題，對這些重現性好的情況做一次優化足矣。而對于柔性大體系，在確保軟硬件沒問題的情況下，追求重現性是沒意義的，而且有時還要特意利用混沌性。例如優化蛋白質，初始結構肯定不是完全準確的，X光衍射實驗本身測出來的結構都有誤差，加氫也會引入誤差，所以既然初始結構不完全精確，那么就有必要人為地引入不同的擾動（對原子坐標隨機做微小移動）以反映出初始結構誤差分布，并由此獲得許多不同的優化后的結構，然后從這些結構中挑選最佳的（比如能量最低的一個）用于后續分析。大部分情況，比如做動力學之前的優化，沒有必要搞得這么精細，畢竟熱運動開始后結構又亂了。然而對于研究比較小的能量差異問題時，如計算受體配體相互作用能，分析多次優化的結果是有意義的，因為數值擾動對優化后能量的影響可能會達到與相互作用能相同的級別。
盡管分子動力學、幾何優化都是迭代過程，SCF計算也是迭代過程，但SCF任務的重現性卻比它們好得多，甚至不同量子化學程序算出的Hartree-Fock能量都十分接近。這主要是因為SCF過程的參數空間（能量相對于軌道系數的空間）結構比較簡單，不像復雜體系勢能面那樣有很多極小點和鞍點，因而對數值擾動敏感性低，即便刻意引入一些擾動，其影響一般不僅不會隨迭代的進行而被放大，相反地還會逐漸降低。很多情況只要能收斂，不管初猜如何，得到的都是參數空間全局最低能量。
4 總結
計算化學在原理上是嚴格遵循決定論的，而之所以在實際研究中時常遇到結果無法重現的問題，一方面是諸多因素對計算產生了擾動，另一方面是執行的任務和體系本身對擾動敏感，呈現混沌效應。筆者常看到有人因為結果無法重現，比如兩次軌跡跑出來RMSD很不一致，于是就苦悶、茫然，并懷疑自己的模擬是錯誤的。只有當了解了任務重現性差的數值原因以及其體現的內在物理意義，才能認清這個問題的本質，甚至利用這個效應獲得更可靠的結果。