導(dǎo)致每次計(jì)算結(jié)果不同的原因

1、1 前言       曾經(jīng)有很多人問，諸如為什么明明輸入文件一樣，但是動力學(xué)軌跡每次跑出來的都不一樣、為什么幾何優(yōu)化結(jié)果相差甚遠(yuǎn)等等，甚至懷疑計(jì)算科學(xué)是否遵循決定論。 實(shí)際上結(jié)果的差異是由于運(yùn)算開始或運(yùn)算過程中各種形式的數(shù)值誤差、隨機(jī)性引起的。這個(gè)是很重要卻不被多數(shù)人重視的問題，它和理論本身、算法、軟件環(huán)境、硬 件環(huán)境等都有密切關(guān)聯(lián)。本文將對這個(gè)問題做一些分析探討。首先先討論數(shù)值算法、軟件、硬件因素是如何導(dǎo)致誤差（錯(cuò)誤）和隨機(jī)性的。然后再看這些問題對計(jì)算化學(xué)會產(chǎn)生何等的影響。      2 數(shù)值誤差（錯(cuò)誤）、隨機(jī)性

2、的根源      2.1 內(nèi)存因?yàn)椴簧偃藢?nèi)存有誤解，所以下面說得多一些。內(nèi)存在數(shù)據(jù)讀、寫的過程中不可避免地會發(fā)生一些錯(cuò)誤，輕則數(shù)據(jù)出現(xiàn)異常（比如一個(gè)變量為1卻成了0），重則死機(jī)、重啟。和內(nèi)存錯(cuò)誤（包括各種類型錯(cuò)誤）相關(guān)的主要有以下幾個(gè)方面：      (1)供電質(zhì)量差。如電壓偏離標(biāo)準(zhǔn)值過多、電流不穩(wěn)。這不僅取決于電源質(zhì)量，還取決于主板好壞（內(nèi)存供電電路）。      (2)過熱。隨著溫度升高內(nèi)存出錯(cuò)的幾率總是增加的，但只有高溫下增加的幅度才最明顯。遙想筆者用i860+RDR

3、AM的日子，不裝個(gè)小風(fēng)扇吹著內(nèi)存的話高負(fù)載下20分鐘左右就死機(jī)。燙手的話應(yīng)該裝上散熱片乃至小風(fēng)扇。      (3)內(nèi)存質(zhì)量。一方面是內(nèi)存設(shè)計(jì)，即布線、PCB板數(shù)、貼片元件等。另一方面取決于內(nèi)存顆粒制造商的品質(zhì)，也有RP的關(guān)系，比如顆粒來自晶圓的邊角還是中心（雜質(zhì)含量不同）。      (4)制程。一般來說制程越精細(xì)可靠性越低，比如22nm制程的比起0.35微米制程的更容易出問題。但不要忽視工藝質(zhì)量也在逐步改進(jìn)，新的精細(xì)制程的產(chǎn)品出錯(cuò)幾率不一定比以往的粗糙制程的要大。     

4、0;(5)宇宙輻射。這看起來比較微妙，但是當(dāng)嚴(yán)重的時(shí)候，如太陽風(fēng)暴劇烈時(shí)，對電子設(shè)備的影響就明顯了。這也是為什么在一些宇航領(lǐng)域（接受的輻射更嚴(yán)重） 的芯片仍然用著設(shè)計(jì)過時(shí)、制程古舊的產(chǎn)品，沒準(zhǔn)兒出現(xiàn)一個(gè)小錯(cuò)大家就沒命了。機(jī)箱內(nèi)電磁輻射也對內(nèi)存穩(wěn)定性有影響。     (6)內(nèi)存參數(shù)。所設(shè)的工作頻率越高、延遲越低出問題的幾率越高，當(dāng)然也要看內(nèi)存類型。     (7)內(nèi)存容量、讀寫次數(shù)：內(nèi)存容量越大、讀寫越頻繁，相同時(shí)間內(nèi)出錯(cuò)的幾率越大。     (8)老化。在越惡劣的環(huán)境使用越長

5、時(shí)間，出錯(cuò)的可能性越大。要注意以上的因素不是分立、絕對的，而是互相關(guān)聯(lián)的、相對的，不能單獨(dú)拿出來對比?？梢?，與內(nèi)存數(shù)據(jù)出錯(cuò)相關(guān)的因素很多，也很難衡量內(nèi)存平均多長時(shí)間出一次錯(cuò)。一個(gè)令我比較認(rèn)可的、由大規(guī)模 統(tǒng)計(jì)得到的說法是一般的使用環(huán)境下（外界環(huán)境正常，計(jì)算機(jī)運(yùn)行穩(wěn)定），平均每年每條內(nèi)存有3%幾率出錯(cuò)（包括各種形式的錯(cuò)誤）。實(shí)際上這個(gè)幾率是很小的， 換句話說，等到機(jī)子被淘汰了也未必能趕上一次內(nèi)存出錯(cuò)。所以不要把內(nèi)存問題過分夸大，比如持續(xù)幾天的計(jì)算，計(jì)算結(jié)果無法重現(xiàn)由內(nèi)存錯(cuò)誤導(dǎo)致的可能性微乎其 微。但是，對于數(shù)據(jù)中心，內(nèi)存出錯(cuò)幾率乘上相應(yīng)內(nèi)存數(shù)目，再考慮數(shù)據(jù)出錯(cuò)的損失，這就是大問題了。內(nèi)存的Err

6、or Checking and Correcting (ECC)顯得十分必要，這是一種數(shù)據(jù)錯(cuò)誤檢測和校正的技術(shù)，每次在寫入數(shù)據(jù)時(shí)同時(shí)寫入額外的數(shù)據(jù)位用作校驗(yàn)，在讀出數(shù)據(jù)時(shí)也同時(shí)把校驗(yàn)位讀出，由此根據(jù) 一定算法分析讀出的數(shù)據(jù)是否和寫入的數(shù)據(jù)相同，若不同則進(jìn)行修復(fù)。這使得內(nèi)存數(shù)據(jù)出錯(cuò)幾率大大降低，但也并非能解決所有問題。有些人誤認(rèn)為ECC內(nèi)存性能 比普通同類型內(nèi)存要低很多，其實(shí)相差甚微。現(xiàn)今服務(wù)器內(nèi)存幾乎都帶ECC，它比普通內(nèi)存略貴，多了額外內(nèi)存顆粒用于儲存校驗(yàn)位，而現(xiàn)今服務(wù)器主板也都支持 ECC（有些主板則只支持ECC內(nèi)存）。個(gè)人計(jì)算沒有絕對的必要使用ECC內(nèi)存。導(dǎo)致內(nèi)存出錯(cuò)的原因?qū)τ贑PU/G

7、PU運(yùn)算出錯(cuò)也基本是一樣的，就不累述了。      2.2 處理器類型       對于不同架構(gòu)的CPU、擁有不同指令集的CPU、乃至GPU，不能期望它們能完全重現(xiàn)同樣的結(jié)果。尤其是浮點(diǎn)運(yùn)算的標(biāo)準(zhǔn)可能不同，比如末位的舍入、極小值 的處理。有些指令在不同處理器中執(zhí)行方式不同，比如先乘一個(gè)數(shù)再加上一個(gè)數(shù)，這兩條指令在有的處理器中會被融合為結(jié)果更精確的“乘加”一條指令。也有一些很特殊的情況，數(shù)值錯(cuò)誤是由CPU設(shè)計(jì)bug引起的。老玩家們應(yīng)該還記得，較早一批的Pentium 60/66Mhz的浮點(diǎn)除法指令有錯(cuò)誤。 

8、     2.3 網(wǎng)絡(luò)數(shù)據(jù)傳輸     對于集群，尤其是分布式計(jì)算，需要考慮網(wǎng)絡(luò)連接引入的錯(cuò)誤。網(wǎng)絡(luò)傳輸都會有一定誤碼率，即在多少位數(shù)據(jù)中出現(xiàn)一位差錯(cuò)，其大小取決于網(wǎng)絡(luò)規(guī)范、傳輸 介質(zhì)、傳輸距離、外界干擾、網(wǎng)絡(luò)適配器質(zhì)量等因素。主流的以太網(wǎng)、infiniband等網(wǎng)絡(luò)標(biāo)準(zhǔn)本身都帶校驗(yàn)幀以保證數(shù)據(jù)完整、準(zhǔn)確性，已很大程度上減 輕了網(wǎng)絡(luò)傳輸引發(fā)的數(shù)據(jù)錯(cuò)誤。     2.4 隨機(jī)數(shù)生成器      在一些程序中要用到隨機(jī)數(shù)生成器，比如Monte Car

9、lo過程、動力學(xué)的初始速度設(shè)定。通常隨機(jī)數(shù)是以不同的隨機(jī)數(shù)算法產(chǎn)生的，常用的比如乘同余法。這些算法都需要一個(gè)稱為“種子”的數(shù)，不同的種子將 產(chǎn)生出不同的隨機(jī)數(shù)序列。也就是說，如果程序中利用了隨機(jī)數(shù)生成器，想重現(xiàn)結(jié)果就必須用相同的種子。但種子往往不允許由用戶設(shè)定，也不是固定不變的，而是 根據(jù)一些規(guī)則產(chǎn)生來保證每次運(yùn)行結(jié)果的隨機(jī)化，比如利用當(dāng)前系統(tǒng)的時(shí)間，對這種情況用戶不要指望結(jié)果會重現(xiàn)。      2.5 庫文件、操作系統(tǒng)、程序版本      不同的數(shù)學(xué)庫中相同的功能在內(nèi)部會使用不同的算法或者不同的編寫方式，即便是同

10、一個(gè)庫，在哪怕相近版本中，也可能由于對代碼進(jìn)行的優(yōu)化、除蟲，導(dǎo)致不同 結(jié)果。有些庫函數(shù)的運(yùn)行結(jié)果甚至取決于運(yùn)行期的情況，比如快速傅里葉變換庫FFTW會在啟動時(shí)檢測哪種算法在當(dāng)前系統(tǒng)下效率最高，然后會在之后一直使用， 這就可能在檢測的時(shí)候由于正在運(yùn)行的其它任務(wù)對系統(tǒng)CPU資源占用不同，使檢測結(jié)論不同，影響后續(xù)算法。有些程序是通過靜態(tài)鏈接產(chǎn)生的，庫文件會被合并在程序可執(zhí)行文件當(dāng)中，或者有些程序雖然使用動態(tài)鏈接，但是 自帶了動態(tài)庫文件，這樣不同系統(tǒng)下結(jié)果的重現(xiàn)性還好說。而有些程序需要依賴于當(dāng)前系統(tǒng)中的庫文件，不同系統(tǒng)自帶版本不同，差異就容易出現(xiàn)了。若平臺不同， 如一個(gè)程序的Linux版本與Windo

11、ws版本，差異會更明顯。不要指望不同版本號的程序能獲得相同結(jié)果，不僅所用的庫往往不同，默認(rèn)參數(shù)、算法也經(jīng)常不同，比如 Gaussian03的幾何優(yōu)化默認(rèn)用的berny方法，而09版后來默認(rèn)用GEDIIS，同一個(gè)輸入文件有時(shí)前者能優(yōu)化成功而后者卻不收斂。即便算法和 參數(shù)都不變，代碼細(xì)節(jié)也可能改變。      2.6 編譯器與編譯選項(xiàng)       編譯器類型及具體版本號、編譯參數(shù)對結(jié)果影響明顯。不同編譯器會對代碼進(jìn)行不同形式的優(yōu)化，有些優(yōu)化是以犧牲精度為代價(jià)的。如果調(diào)用了程序語言標(biāo)準(zhǔn)中定義 的標(biāo)準(zhǔn)數(shù)學(xué)函數(shù)，比如開方、求

12、三角函數(shù)、求對數(shù)等，不同編譯器會用不同的內(nèi)建數(shù)學(xué)庫，顯然也會帶來差異。不同編譯器對于語言標(biāo)準(zhǔn)中沒有定義的情況的處理也 不同，比如未初始化的變量，有些編譯器一律將之初始化為0，而有些編譯器則不這么做，相應(yīng)內(nèi)存地址位置目前是什么數(shù)值就是什么數(shù)值，由于這個(gè)數(shù)值無法預(yù) 測，每次運(yùn)行時(shí)可能造成結(jié)果不同。來看一個(gè)簡單的比較，用Multiwfn 2.02對一個(gè)水分子全空間積分電子密度拉普拉斯值的結(jié)果。Intel visual fortran與CVF結(jié)果有所不同，在不同優(yōu)化參數(shù)下結(jié)果也有所不同。IVF12  0d/O1       

13、60;-2.603320357104659E-005IVF12  O3            -2.603320356662977E-005CVF6.5 debug/release -2.603320356424300E-005      2.7 并行執(zhí)行時(shí)的運(yùn)算順序雖然諸如a+b+c=a+(b+c)在數(shù)學(xué)上是成立的，但是在計(jì)算機(jī)執(zhí)行中未必成立。比如在CVF6.5里執(zhí)行這樣的語句a=sin(4D0)b=sin(99D0)c=sin(2D0)write(*,

14、*) a+b+c-(a+(b+c)得 到的結(jié)果不是0.000000000000000E+000，而是1.110223024625157E-016，如果結(jié)果被乘上一個(gè)比較大的數(shù)，那么這 個(gè)差異就不能被忽略了，可見計(jì)算機(jī)運(yùn)算未必滿足結(jié)合律。在并行運(yùn)算中，經(jīng)常要把一部分工作拆成幾個(gè)線程以分配給多個(gè)處理器并行執(zhí)行，然后將算出來的結(jié)果累 加到一起。通常哪個(gè)線程先算完，哪個(gè)結(jié)果就先被累加上。然而實(shí)際計(jì)算中各個(gè)線程完成的先后順序是不確定的，每次程序執(zhí)行時(shí)都可能不同，因此數(shù)據(jù)累加的順序 不一樣，這就造成結(jié)果出現(xiàn)一定不確定性。另外，有的程序中使用動態(tài)均衡負(fù)載以獲得更高的效率，這導(dǎo)致每個(gè)線程執(zhí)行的內(nèi)容也是有隨機(jī)性

15、的，使結(jié)果的重現(xiàn)更加 困難。      舉個(gè)實(shí)例，在Gaussian09中用4核計(jì)算一個(gè)體系單點(diǎn)能，執(zhí)行三次，最后一次迭代的能量分別為-234.579551993147-234.579551993143-234.579551993144      結(jié)果在末位有所差異。如果只用nproc=1來串行計(jì)算，則結(jié)果每次都能重現(xiàn)。盡管并行時(shí)每次的結(jié)果差異看似微不足道，以a.u.為單位一般精確到小數(shù) 后5、6位足矣，然而這個(gè)差異在其它一些任務(wù)中會被逐漸放大，甚至可能達(dá)到影響定性結(jié)論的地步。有些人以為通過增加積分精度、收斂精度之

16、類能減輕結(jié)果的隨 機(jī)性，實(shí)際上從原理上講根本于事無補(bǔ)，根本不是一碼事（僅在極個(gè)別情況下能有改善）。      2.8 文件格式      不同文件格式的精度，以及格式轉(zhuǎn)換過程中的信息丟失是需要注意的。浮點(diǎn)數(shù)據(jù)在計(jì)算機(jī)中通過二進(jìn)制形式記錄，而在磁盤上記錄的文件往往是文本格式，這便于 人類閱讀和第三方程序處理。在內(nèi)存中占8字節(jié)的雙精度變量若想完整以文本方式記錄下來起碼要用22字節(jié)，這不僅太占空間，閱讀也不方便，因此往往故意減少 有效數(shù)字位數(shù)。比如Gaussian的chk文件是二進(jìn)制文件，浮點(diǎn)變量以二進(jìn)制形式精確保存，有

17、效數(shù)字位數(shù)約16。而轉(zhuǎn)換成fch后，浮點(diǎn)數(shù)有效數(shù)字就 被截為9位，若用unfchk再將之轉(zhuǎn)化回chk，則此時(shí)的chk的數(shù)據(jù)精度比起之前的chk已經(jīng)損失了。即便是二進(jìn)制文件也并非都是完整精度形式記錄數(shù) 據(jù)，比如Gromacs的xtc文件也是二進(jìn)制形式，而變量記錄精度可控（xtc_precision參數(shù)），默認(rèn)的記錄精度低于單精度浮點(diǎn)變量（這也是 為什么Gromacs軌跡文件比較小的原因），因此不要指望rerun的過程，即重新算軌跡中每一幀的勢能的值會與之前動力學(xué)過程中輸出的一致。有些程序 在導(dǎo)出/導(dǎo)入的過程中還可能會對結(jié)構(gòu)、導(dǎo)數(shù)等等進(jìn)行坐標(biāo)變換，對單位進(jìn)行轉(zhuǎn)換等操作，也會引入數(shù)值誤差。 

18、;      2.9 操作過程       一些操作過程的細(xì)節(jié)差異往往不被注意，對最終結(jié)果的影響可能卻不小。例如對晶體測得的結(jié)構(gòu)加氫，很多程序加氫結(jié)果表面看似都一樣，位置沒有區(qū)別，但是查看 文件中記錄的坐標(biāo)，卻總有細(xì)微差別。再比如模擬一個(gè)NMR測得的體系，pdb文件通常會有很多幀，往往每一幀差異很小，似乎取哪幀是隨意的，但實(shí)際上用于 模擬會得到不同結(jié)果。還比如Gaussian中設(shè)定不同內(nèi)存大小，看似和計(jì)算結(jié)果無關(guān)，但實(shí)際上Gaussian會自動根據(jù)分配的內(nèi)存容量選擇不同的積分 計(jì)算、儲存方式，也會引入差異。  &#

19、160;   總結(jié)一下，上面的討論中，數(shù)值差異可以歸納為兩部分，一種是可重現(xiàn)性差異，也可以叫靜態(tài)差異，包括硬件配置、編譯過程、庫文件、文件操作、操作過程、程 序版本，只要運(yùn)行條件完全一致就可以徹底避免；另一種是不可重現(xiàn)性差異，或者叫動態(tài)差異，它取決于運(yùn)行期的實(shí)際條件，有的可以通過一定對策避免，比如把并 行運(yùn)算改為串行，但有的很難完全消除，尤其是硬件偶發(fā)錯(cuò)誤。      3 數(shù)值誤差、隨機(jī)性對計(jì)算化學(xué)結(jié)果的影響      從前面簡單例子看到，誤差或者隨機(jī)性對結(jié)果的影響甚微，在默認(rèn)的輸出精度下甚至看

20、不出任何差異。然而對于某些任務(wù)情況卻不是如此，數(shù)值的微小擾動（下面 用擾動這個(gè)詞代表前述各種因素造成的微小數(shù)值差異）可能造成結(jié)果有明顯定量乃至定性的差異，這是因?yàn)樵谶@些任務(wù)中擾動的效應(yīng)會被放大。放大的方式主要可以 歸結(jié)為兩類，一類是數(shù)值運(yùn)算方式，另一類是體系和任務(wù)的混沌效應(yīng)。     3.1 數(shù)值運(yùn)算方式對擾動的放大     一些任務(wù)中，由于算法的原因，初始的擾動會在結(jié)果中被放大，典型的例子是數(shù)值差分。在幾何優(yōu)化、頻率計(jì)算、超極化率計(jì)算等任務(wù)中，若理論方法在當(dāng)前程序中不支持能量解析導(dǎo)數(shù)計(jì)算，就必須要通過數(shù)值差分來近

21、似計(jì)算導(dǎo)數(shù)。x處一階導(dǎo)數(shù)計(jì)算方法為：    d1(x)=( E(x+)-E(x-) )/(2)    其中是差分步長。假設(shè)最糟的情況，計(jì)算E(x+)和E(x-)時(shí)都引入了大小為k的擾動，且符號相反而無法被抵消，則此時(shí)的一階導(dǎo)數(shù)寫為    d1'(x)=( E(x+)-E(x-)+2k )/(2)那么由于擾動，造成一階導(dǎo)數(shù)誤差為d1'(x)-d1(x)=k/。由于必須足夠小以保證導(dǎo)數(shù)計(jì)算準(zhǔn)確，比如為0.0001，那么初始擾動k就會被放大10000倍，對結(jié)果造成了10000k的影響！    &

22、#160;再來看差分計(jì)算二階導(dǎo)數(shù)     d2(x)=( d1(x+)-d1(x-) )/(2)     由于計(jì)算d1的時(shí)候已經(jīng)引入了k/的誤差，則實(shí)際計(jì)算出的二階導(dǎo)數(shù)為    d2'(x)=( d1(x+)-d1(x-)+2k/ )/(2)    二階導(dǎo)數(shù)的計(jì)算誤差為d2'(x)-d2(x)=k/2，若=0.0001，則初始擾動k就會被放大100000000倍！可見，差分計(jì)算n階導(dǎo)數(shù) 時(shí)由于初始擾動值k最多可能導(dǎo)致結(jié)果產(chǎn)生k/n的誤差。可想而知，計(jì)算三階導(dǎo)數(shù)

23、和四階導(dǎo)數(shù)時(shí)，比如有限場方法純數(shù)值計(jì)算一階和二階超極化率時(shí)，誤差已 經(jīng)被放大到不可忽略的地步了。所以，取得越小，由于誤差會越大，結(jié)果未必越精確，必須選擇最合適的值。通過選擇合理的算法、良好的編程習(xí)慣，可以使一些任務(wù)中誤差的放大盡可能小，比如避免大數(shù)除小數(shù)之類的情況。這就是編程者的事情了。      3.2 體系、任務(wù)的混沌性對擾動的放大混沌效應(yīng)對擾動的放大是特別需要關(guān)注的?；煦缧?yīng)粗俗地說就是捷克一個(gè)人從體內(nèi)釋放一股不圣潔的氣體可能造 成韓國一場臺風(fēng)，用古話來講則是“失之毫厘，謬以千里”?；煦缧?yīng)對數(shù)值計(jì)算巨大影響的發(fā)現(xiàn)，可追溯到1959年洛倫茲用計(jì)算機(jī)

24、進(jìn)行氣象模擬（其實(shí)更早之 前有人在計(jì)算天體運(yùn)行軌道時(shí)也多多少少意識到了這點(diǎn)）。他原先模擬時(shí)一直使用十幾位精度的浮點(diǎn)數(shù)保存數(shù)據(jù)，后來有一次為了更細(xì)致地分析某段時(shí)間的數(shù)據(jù)，他 用以前的數(shù)據(jù)為初值重算了一遍，為了打印的數(shù)據(jù)美觀，這次他將小數(shù)精度只保留三位。正是這微小誤差的引入，導(dǎo)致模擬出的兩個(gè)月之后的氣象結(jié)果和原先大相徑 庭！在計(jì)算化學(xué)中，有兩類任務(wù)混沌效應(yīng)很明顯，一個(gè)是分子動力學(xué)，另一個(gè)是幾何優(yōu)化。      3.3 擾動對分子動力學(xué)的影響分子動力學(xué)的混沌效應(yīng)，尤其是對蛋白質(zhì)折疊過程的影響被研究得比較多，如PROTEINS,29,417。 動力學(xué)過程中初始擾

25、動會隨著模擬時(shí)間指數(shù)型放大，這并不是由于數(shù)值計(jì)算問題所導(dǎo)致的，而是體系內(nèi)在性質(zhì)的體現(xiàn)。若在t=0時(shí)對體系坐標(biāo)引入微擾，則模擬 到t時(shí)軌跡與未引入微擾時(shí)的差異可寫為diff(t)=*exp(t)，其中是依賴于體系的Lyapunov指數(shù)。擾動的放大速度相當(dāng)快，研究發(fā)現(xiàn) 對蛋白質(zhì)結(jié)構(gòu)引入RMSD=1*10-9埃的擾動，經(jīng)過2ps的模擬就能被放大到RMSD=1埃。假設(shè)在模擬中觀察到某時(shí)刻有一個(gè)氨基酸側(cè)鏈發(fā)生了扭 轉(zhuǎn)，如果在此10ps前引入10次RMSD=0.001的不同方式的擾動，并分別跑10次動力學(xué)，則只有一次能夠重現(xiàn)側(cè)鏈扭轉(zhuǎn)。可見，動力學(xué)過程對初始條 件十分敏感，各種因素帶來的微小數(shù)值擾動都足矣使

26、動力學(xué)軌跡發(fā)生很大分歧。更何況實(shí)際模擬中擾動并不僅僅是在初始條件中出現(xiàn)，在模擬過程中還可能會不斷納 入。      注意對于蛋白質(zhì)體系，擾動不是隨模擬的進(jìn)行而無限放大的，因?yàn)榈鞍椎脑邮艿搅龅氖`，且模擬時(shí)的溫度有限，因此勢能面可及的范圍是有限的，故初始擾 動造成的差異在迅速放大后會遇到上限。對于柔性體系，由于勢能面可及范圍會加更大，可以預(yù)見差異放大的上限會比剛性體系更高。一些辦法可以提高軌跡重現(xiàn)性。提高浮點(diǎn)數(shù)精度是重要策略之一，這可以減少擾動量的數(shù)量級。有人抱怨 gromacs軌跡重現(xiàn)性不如Amber，這很大程度上是由于gromacs默認(rèn)是單精度編譯

27、，而Amber默認(rèn)用雙精度。不過，提高精度并不解決根本問 題，只不過是讓軌跡的歧化慢一些罷了，對于長時(shí)間的動力學(xué)模擬結(jié)果的偏差仍然是無法預(yù)測的，除非能讓浮點(diǎn)精度無窮高，但這顯然是不可能的。若想讓軌跡完全 重現(xiàn)，通用的辦法是將第一節(jié)涉及的所有可能帶來擾動的因素都去除。而有些程序也專門考慮了軌跡重現(xiàn)問題，在gromacs里的mdrun有-reprod 選項(xiàng)以幫助軌跡重現(xiàn)。另外還有人弄了個(gè)固定精度的分子動力學(xué)方法，不僅可以完全重現(xiàn)軌跡，還可以將時(shí)間反轉(zhuǎn)從末態(tài)得到初始結(jié)構(gòu)，完全體現(xiàn)出分子動力學(xué)的決 定性特點(diǎn)。      動力學(xué)軌跡的重現(xiàn)究竟有無意義？如果只是出于

28、計(jì)算目的，比如一段軌跡文件弄丟了，或者希望以更高的頻率輸出體系能量，那么重新跑出一次完全一樣的軌跡是 有意義的。然而，從理論角度上來講，軌跡的嚴(yán)格的重現(xiàn)性是毫無意義的。動力學(xué)過程就是在體系相空間（坐標(biāo)+動量空間）中游歷的過程，如果跑兩次得到了截然 不同的軌跡，只不過是說明在相空間不同位置區(qū)域游歷罷了，更沒有誰對誰錯(cuò)之分。對于一個(gè)平衡態(tài)的模擬，體系的性質(zhì)（比如內(nèi)能、擴(kuò)散系數(shù)）就是動力學(xué)過程時(shí)間的平均，兩次跑出來軌跡相不相 同完全無所謂，只要時(shí)間足夠長，它們都足以充分遍歷相空間，結(jié)果會是等價(jià)的。對于有限時(shí)間的模擬，為了加強(qiáng)某些區(qū)域的采樣，實(shí)際上還有人專門利用動力學(xué)軌 跡的混沌效應(yīng)，給初始速度分配不

29、同的值，分別做動力學(xué)模擬并將結(jié)果合并。而對于非平衡過程的模擬，比如說蛋白質(zhì)折疊過程，軌跡（折疊路徑）的重現(xiàn)性同樣毫無意義。因?yàn)榈鞍椎恼郫B過 程本來就是混沌的，初始的微小擾動會使折疊路徑相差迥異。在現(xiàn)實(shí)中每次蛋白的初始條件絕不可能完全一致，因此同樣的折疊路徑在現(xiàn)實(shí)中無法重現(xiàn)，或者說現(xiàn)實(shí) 中沒有可能發(fā)生兩次完全相同的折疊過程，那么苛求在多次模擬中重現(xiàn)相同的軌跡有何意義？實(shí)際上不僅不應(yīng)該重現(xiàn)，折疊路徑的多樣性才是真正感興趣的，也很適 合通過概率統(tǒng)計(jì)來研究，這就需要做很多次動力學(xué)模擬，以得到不同的折疊路徑并進(jìn)行分析。正好，由于各種數(shù)值因素的介入，哪怕相同的初始條件都會產(chǎn)生不同的 折疊路徑，都省得修改輸

30、入文件了，直接反復(fù)跑多次就行了。注意盡管每次的軌跡十分不同，但最終卻一定會折疊到同樣的最穩(wěn)定結(jié)構(gòu)，這并不違背混沌效應(yīng)，而是 由于每一次折疊過程都會感受到相同的“吸引子”，即趨向成為最穩(wěn)結(jié)構(gòu)的一種特殊勢場。折疊過程可以比喻為站在一個(gè)坑坑洼洼的大坑（對應(yīng)自由能面的形狀）旁 邊往坑里扔石子，每次石子掉落的軌跡都很不同，但最終都掉到坑底。       這里順便談?wù)勔粋€(gè)相關(guān)問題，也就是單個(gè)動力學(xué)模擬的軌跡是否足以說明問題。上面說了，一條折疊軌跡不能反映全部的折疊過程，同樣地，研究其它過程，原則上 也需要用多條軌跡來統(tǒng)計(jì)地說明問題。比如研究在結(jié)合配體后激酶的活性位點(diǎn)

31、打開的過程，在軌跡中發(fā)現(xiàn)經(jīng)過530ps模擬后活性位點(diǎn)打開了，那么就真的能說明 在現(xiàn)實(shí)中配體結(jié)合后經(jīng)過約530ps后活性位點(diǎn)就能打開？非也，模擬10次，會得到10個(gè)不同的結(jié)果，恰好都出現(xiàn)在530ps是不可能的。經(jīng)過10次模擬 后，假設(shè)根據(jù)統(tǒng)計(jì)結(jié)果可說明位點(diǎn)打開的時(shí)間主要分布在300ps1.0ns區(qū)間內(nèi)，且最容易在500ps附近打開，則這樣的統(tǒng)計(jì)性結(jié)論比起從單個(gè)軌跡中 下的結(jié)論有意義、可信得多?；煦缧?yīng)導(dǎo)致動力學(xué)過程中各種“事件”發(fā)生的時(shí)刻的波動范圍極大，或者說每種事件在不同的時(shí)間范圍都有一定出現(xiàn)概率，只分析單 一軌跡就造成了一葉障目，結(jié)論的不可靠性會較大。然而目前絕大多數(shù)模擬都是基于單一軌跡的，

32、很少有人對同一問題做多軌跡統(tǒng)計(jì)分析，主要是多軌跡的計(jì)算量太 大或者研究的內(nèi)容不需要那么深入詳細(xì)，其次是很多人沒意識到這個(gè)問題，還有少數(shù)人是因?yàn)橥稒C(jī)取巧。比如論證一個(gè)抑制劑的作用，可能模擬反復(fù)跑了7次都無法 論證抑制劑管用，當(dāng)跑到第8次時(shí)慶幸地發(fā)現(xiàn)能說明抑制劑管用，于是他只用這一次的軌跡說事，另外跑的7次雖然沒成功卻都避而不談，可想而知這抑制劑的實(shí)際 功效如何了.      3.4 擾動對幾何優(yōu)化的影響相較分子動力學(xué)的混沌效應(yīng)，幾何優(yōu)化的混沌效應(yīng)被關(guān)注的相對較少，有興趣者可參看J Comput Aided Mol Des,22,39，討論數(shù)值誤差對分子力場

33、優(yōu)化結(jié)果的影響。之所以幾何優(yōu)化也有明顯混沌效應(yīng)，不妨將幾何優(yōu)化看作是0K下的分子動力學(xué)，尤其對于最陡下降法優(yōu)化，這個(gè) 比喻最為恰當(dāng)。最陡下降法優(yōu)化過程類似把一個(gè)玻璃珠輕輕放到陡峭的山上某斜坡處，讓它自動滾落。雖然0K下沒有動能（準(zhǔn)確來說尚有振動能），通常不會導(dǎo)致 引入擾動后幾何優(yōu)化過程的軌跡像動力學(xué)軌跡那樣隨步數(shù)增長分歧得那么快，但是，由于沒有了動能以越過勢壘，體系將沒機(jī)會感受到全局的吸引子（全局最穩(wěn)結(jié) 構(gòu)），而只能感受到離初始位置最近的局部吸引子（勢能面局部極小點(diǎn)），這導(dǎo)致了幾何優(yōu)化的最終結(jié)果的重現(xiàn)性可能比動力學(xué)更差。對于小體系還好，勢能面結(jié)構(gòu) 比較簡單、光滑，然而極小點(diǎn)、鞍點(diǎn)的數(shù)目是隨著原

34、子數(shù)目增加呈指數(shù)型增長的，對于大體系，勢能面的復(fù)雜程度是難以想象的，微小的擾動可能就會使結(jié)構(gòu)收斂到 很不同的極小點(diǎn)。之所以平時(shí)在量子化學(xué)的幾何優(yōu)化中優(yōu)化的結(jié)果重現(xiàn)性都很好，沒有顯示出太大不可重現(xiàn)性，主要是因?yàn)榱炕幚?的體系一般都比較小，通常少于50個(gè)原子，所需優(yōu)化步數(shù)一般也較少，通常少于100步。不僅誤差來不及充分放大，而且由于相對簡單的勢能面的形狀以及它起 到的束縛作用，誤差也不容易放大。然而因數(shù)值誤差帶來的分歧偶爾仍可能顯著，記得有人在一臺機(jī)子上某個(gè)體系優(yōu)化幾十步成功了，換了臺機(jī)子卻不收斂了。重現(xiàn)性好壞與初猜坐標(biāo)關(guān)系很大，在山腳下釋放一個(gè)珠子，和在山尖上釋放一個(gè)珠子，明顯珠子最終位置在后者的 情況中受初始位置影響更大。所以優(yōu)化時(shí)如果能給出一個(gè)好的初猜，比如在極小點(diǎn)附近的二次區(qū)域內(nèi)（勢能面可以近似用二次函數(shù)描述的區(qū)域），那么優(yōu)化幾乎一定 收斂到那個(gè)極小點(diǎn)，數(shù)值誤差因素不會導(dǎo)致收斂到別處去。然而如果在很接近過渡態(tài)的位