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

1、請(qǐng)下載支持！1 前言曾經(jīng)有很多人問，諸如為什么明明輸入文件一樣，但是動(dò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)。本文將對(duì)這個(gè)問題做一些分析探討。首先先討論數(shù)值算法、軟件、硬件因素是如何導(dǎo)致誤差（錯(cuò)誤）和隨機(jī)性的。然后再看這些問題對(duì)計(jì)算化學(xué)會(huì)產(chǎn)生何等的影響。2 數(shù)值誤差（錯(cuò)誤）、隨機(jī)性的根源2.1 內(nèi)存因?yàn)椴簧偃藢?duì)內(nèi)存有誤解，所以下面說得多一些。內(nèi)存在數(shù)據(jù)讀、寫的過程中不可避免地會(huì)

2、發(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+RDRAM的日子，不裝個(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ì)含

3、量不同）。(4)制程。一般來說制程越精細(xì)可靠性越低，比如22nm制程的比起0.35微米制程的更容易出問題。但不要忽視工藝質(zhì)量也在逐步改進(jìn)，新的精細(xì)制程的產(chǎn)品出錯(cuò)幾率不一定比以往的粗糙制程的要大。(5)宇宙輻射。這看起來比較微妙，但是當(dāng)嚴(yán)重的時(shí)候，如太陽風(fēng)暴劇烈時(shí)，對(duì)電子設(shè)備的影響就明顯了。這也是為什么在一些宇航領(lǐng)域(接受的輻射更嚴(yán)重)的芯片仍然用著設(shè)計(jì)過時(shí)、制程古舊的產(chǎn)品，沒準(zhǔn)兒出現(xiàn)一個(gè)小錯(cuò)大家就沒命了。機(jī)箱內(nèi)電磁輻射也對(duì)內(nèi)存穩(wěn)定性有影響。(6)內(nèi)存參數(shù)。所設(shè)的工作頻率越高、延遲越低出問題的幾率越高，當(dāng)然也要看內(nèi)存類型。(7)內(nèi)存容量、讀寫次數(shù)：內(nèi)存容量越大、讀寫越頻繁，相同時(shí)間內(nèi)出錯(cuò)的幾率越

4、大。(8)老化。在越惡劣的環(huán)境使用越長(zhǎng)時(shí)間，出錯(cuò)的可能性越大。要注意以上的因素不是分立、絕對(duì)的，而是互相關(guān)聯(lián)的、相對(duì)的，不能單獨(dú)拿出來對(duì)比?？梢姡c內(nèi)存數(shù)據(jù)出錯(cuò)相關(guān)的因素很多，也很難衡量?jī)?nèi)存平均多長(zhǎng)時(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)致的可能性微乎其微。但是，對(duì)于數(shù)據(jù)中心，內(nèi)存出錯(cuò)幾率乘上相應(yīng)內(nèi)存數(shù)目，再考慮數(shù)據(jù)出錯(cuò)的損

5、失，這就是大問題了。內(nèi)存的ErrorCheckingandCorrecting(EC唾得十分必要，這是一種數(shù)據(jù)錯(cuò)誤檢測(cè)和校正的技術(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í)相差甚微?，F(xiàn)今服務(wù)器內(nèi)存幾乎都帶ECC，它比普通內(nèi)存略貴，多了額外內(nèi)存顆粒用于儲(chǔ)存校驗(yàn)位，而現(xiàn)今服務(wù)器主板也都支持ECC（有些主板則只支持ECC內(nèi)存）。個(gè)人計(jì)算沒有絕對(duì)的必要使用ECC內(nèi)存。導(dǎo)致內(nèi)存出錯(cuò)的原因

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

7、數(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ù)生成器，比如MonteCarlo過程、動(dòng)力學(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)

8、行結(jié)果的隨機(jī)化，比如利用當(dāng)前系統(tǒng)的時(shí)間，對(duì)這種情況用戶不要指望結(jié)果會(huì)重現(xiàn)。2.5 庫文件、操作系統(tǒng)、程序版本不同的數(shù)學(xué)庫中相同的功能在內(nèi)部會(huì)使用不同的算法或者不同的編寫方式，即便是同一個(gè)庫，在哪怕相近版本中，也可能由于對(duì)代碼進(jìn)行的優(yōu)化、除蟲，導(dǎo)致不同結(jié)果。有些庫函數(shù)的運(yùn)行結(jié)果甚至取決于運(yùn)行期的情況，比如快速傅里葉變換庫FFTW會(huì)在啟動(dòng)時(shí)檢測(cè)哪種算法在當(dāng)前系統(tǒng)下效率最高，然后會(huì)在之后一直使用，這就可能在檢測(cè)的時(shí)候由于正在運(yùn)行的其它任務(wù)對(duì)系統(tǒng)CPU資源占用不同，使檢測(cè)結(jié)論不同，影響后續(xù)算法。有些程序是通過靜態(tài)鏈接產(chǎn)生的，庫文件會(huì)被合并在程序可執(zhí)行文件當(dāng)中，或者有些程序雖然使用動(dòng)態(tài)鏈接，但是自帶了動(dòng)

9、態(tài)庫文件，這樣不同系統(tǒng)下結(jié)果的重現(xiàn)性還好說。而有些程序需要依賴于當(dāng)前系統(tǒng)中的庫文件，不同系統(tǒng)自帶版本不同，差異就容易出現(xiàn)了。若平臺(tái)不同，如一個(gè)程序的Linux版本與Windows版本，差異會(huì)更明顯。不要指望不同版本號(hào)的程序能獲得相同結(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)編譯器類型及具體版本號(hào)、編譯參數(shù)對(duì)結(jié)果影響明顯。不同編譯器會(huì)對(duì)代碼進(jìn)行不同形式的優(yōu)化，有些優(yōu)化是以犧牲精度為代價(jià)

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

11、2.6662977E-005CVF6.5debug/release-2.6424300E-0052.7 并行執(zhí)行時(shí)的運(yùn)算順序雖然諸如a+b+c=a+(b+/數(shù)學(xué)上是成立的，但是在計(jì)算機(jī)執(zhí)行中未必成立。比如在CVF6.5里執(zhí)行這樣的語句a=sin(4D0)b=sin(99D0)c=sin(2D0)write(*,*)a+b+c-(a+(b+c)彳導(dǎo)至U的結(jié)果不是0.000000000000000E+00Q而是1.4625157E-016,如果結(jié)果被乘上一個(gè)比較大的數(shù)，那么這個(gè)差異就不能被忽略了，可見計(jì)算機(jī)運(yùn)算未必滿足結(jié)合律。在并行運(yùn)算中，經(jīng)常要把一部分工作拆成幾個(gè)線程以分配給多個(gè)處理器并行執(zhí)行，

12、然后將算出來的結(jié)果累加到一起。通常哪個(gè)線程先算完，哪個(gè)結(jié)果就先被累加上。然而實(shí)際計(jì)算中各個(gè)線程完成的先后順序是不確定的，每次程序執(zhí)行時(shí)都可能不同，因此數(shù)據(jù)累加的順序不一樣，這就造成結(jié)果出現(xiàn)一定不確定性。另外，有的程序中使用動(dòng)態(tài)均衡負(fù)載以獲得更高的效率，這導(dǎo)致每個(gè)線程執(zhí)行的內(nèi)容也是有隨機(jī)性的，使結(jié)果的重現(xiàn)更加困難。舉個(gè)實(shí)例，在Gaussian09中用4核計(jì)算一個(gè)體系單點(diǎn)能，執(zhí)行三次，最后一次迭代的能量分別為-234.3147-234.3143-234.3144結(jié)果在末位有所差異。如果只用nproc=1來串行計(jì)算，則結(jié)果每次都能重現(xiàn)。盡管并行時(shí)每次的結(jié)果差異看似微不足道，以a.u.為單位一般精確到

13、小數(shù)后5、6位足矣，然而這個(gè)差異在其它一些任務(wù)中會(huì)被逐漸放大，甚至可能達(dá)到影響定性結(jié)論的地步。有些人以為通過增加積分精度、收斂精度之類能減輕結(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)變量以

14、二進(jìn)制形式精確保存，有效數(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的過程，即重新算軌跡中每一幀的勢(shì)能的值會(huì)與之前動(dòng)力學(xué)過程中輸出的一致。有些程序在導(dǎo)出/導(dǎo)入的過程中還可能會(huì)對(duì)結(jié)構(gòu)、導(dǎo)數(shù)等等進(jìn)行坐標(biāo)變換，對(duì)單位進(jìn)行轉(zhuǎn)換等操作，也會(huì)引入數(shù)值誤

15、差。2.9 操作過程一些操作過程的細(xì)節(jié)差異往往不被注意，對(duì)最終結(jié)果的影響可能卻不小。例如對(duì)晶體測(cè)得的結(jié)構(gòu)加氫，很多程序加氫結(jié)果表面看似都一樣，位置沒有區(qū)別，但是查看文件中記錄的坐標(biāo)，卻總有細(xì)微差別。再比如模擬一個(gè)NMR測(cè)得的體系，pdb文件通常會(huì)有很多幀，往往每一幀差異很小，似乎取哪幀是隨意的，但實(shí)際上用于模擬會(huì)得到不同結(jié)果。還比如Gaussian中設(shè)定不同內(nèi)存大小，看似和計(jì)算結(jié)果無關(guān)，但實(shí)際上Gaussian會(huì)自動(dòng)根據(jù)分配的內(nèi)存容量選擇不同的積分計(jì)算、儲(chǔ)存方式，也會(huì)引入差異?？偨Y(jié)一下，上面的討論中，數(shù)值差異可以歸納為兩部分，一種是可重現(xiàn)性差異，也可以叫靜態(tài)差異，包括硬件配置、編譯過程、庫文件

16、、文件操作、操作過程、程序版本，只要運(yùn)行條件完全一致就可以徹底避免；另一種是不可重現(xiàn)性差異，或者叫動(dòng)態(tài)差異，它取決于運(yùn)行期的實(shí)際條件，有的可以通過一定對(duì)策避免，比如把并行運(yùn)算改為串行，但有的很難完全消除，尤其是硬件偶發(fā)錯(cuò)誤。3 數(shù)值誤差、隨機(jī)性對(duì)計(jì)算化學(xué)結(jié)果的影響從前面簡(jiǎn)單例子看到，誤差或者隨機(jī)性對(duì)結(jié)果的影響甚微，在默認(rèn)的輸出精度下甚至看不出任何差異。然而對(duì)于某些任務(wù)情況卻不是如此，數(shù)值的微小擾動(dòng)(下面用擾動(dòng)這個(gè)詞代表前述各種因素造成的微小數(shù)值差異)可能造成結(jié)果有明顯定量乃至定性的差異，這是因?yàn)樵谶@些任務(wù)中擾動(dòng)的效應(yīng)會(huì)被放大。放大的方式主要可以歸結(jié)為兩類，一類是數(shù)值運(yùn)算方式，另一類是體系和任務(wù)

17、的混沌效應(yīng)。3.1 數(shù)值運(yùn)算方式對(duì)擾動(dòng)的放大一些任務(wù)中，由于算法的原因，初始的擾動(dòng)會(huì)在結(jié)果中被放大，典型的例子是數(shù)值差分。在幾何優(yōu)化、頻率計(jì)算、超極化率計(jì)算等任務(wù)中，若理論方法在當(dāng)前程序中不支持能量解析導(dǎo)數(shù)計(jì)算，就必須要通過數(shù)值差分來近似計(jì)算導(dǎo)數(shù)。x處一階導(dǎo)數(shù)計(jì)算方法為：d1(x)=(E(x+-E儀必)/(2A)其中是差分步長(zhǎng)。假設(shè)最糟的情況，計(jì)算E(x+邪E(x-時(shí)都引入了大小為k的擾動(dòng)，且符號(hào)相反而無法被抵消，則此時(shí)的一階導(dǎo)數(shù)寫為d1'(x)=(E(x+AE(x-A)+2k)/(2那么由于擾動(dòng)，造成一階導(dǎo)數(shù)誤差為d1'(x)-d1(x)=k/。由于必須足夠小以保證導(dǎo)數(shù)計(jì)算準(zhǔn)

18、確，比如為0.0001,那么初始擾動(dòng)k就會(huì)被放大10000倍，對(duì)結(jié)果造成了10000k的影響！再來看差分計(jì)算二階導(dǎo)數(shù)d2(x)=(d1(x+-dlXX-A)/(2A)由于計(jì)算di的時(shí)候已經(jīng)引入了k/加勺誤差，則實(shí)際計(jì)算出的二階導(dǎo)數(shù)為d2'(x)=(d1(x+-di(X-A)+2k/A)/(2A)二階導(dǎo)數(shù)的計(jì)算誤差為d2'(x)-d2(x)=k/若=0.0001則初始擾動(dòng)k就會(huì)被放大0倍！可見，差分計(jì)算n階導(dǎo)數(shù)時(shí)由于初始擾動(dòng)值k最多可能導(dǎo)致結(jié)果產(chǎn)生k/八時(shí)勺誤差。可想而知，計(jì)算三階導(dǎo)數(shù)和四階導(dǎo)數(shù)時(shí)，比如有限場(chǎng)方法純數(shù)值計(jì)算一階和二階超極化率時(shí)，誤差已經(jīng)被放大到不可忽略的地步了。所

19、以，(取得越小，由于誤差會(huì)越大，結(jié)果未必越精確，必須選擇最合適的值。通過選擇合理的算法、良好的編程習(xí)慣，可以使一些任務(wù)中誤差的放大盡可能小，比如避免大數(shù)除小數(shù)之類的情況。這就是編程者的事情了。3.2 體系、任務(wù)的混沌性對(duì)擾動(dòng)的放大混沌效應(yīng)對(duì)擾動(dòng)的放大是特別需要關(guān)注的?；煦缧?yīng)粗俗地說就是捷克一個(gè)人從體內(nèi)釋放一股不圣潔的氣體可能造成韓國(guó)一場(chǎng)臺(tái)風(fēng)，用古話來講則是“失之毫厘，謬以千里”?；煦缧?yīng)對(duì)數(shù)值計(jì)算巨大影響的發(fā)現(xiàn)，可追溯到1959年洛倫茲用計(jì)算機(jī)進(jìn)行氣象模擬(其實(shí)更早之前有人在計(jì)算天體運(yùn)行軌道時(shí)也多多少少意識(shí)到了這點(diǎn))。他原先模擬時(shí)一直使用十幾位精度的浮點(diǎn)數(shù)保存數(shù)據(jù)，后來有一次為了更細(xì)致地分析

20、某段時(shí)間的數(shù)據(jù)，他用以前的數(shù)據(jù)為初值重算了一遍，為了打印的數(shù)據(jù)美觀，這次他將小數(shù)精度只保留三位。正是這微小誤差的引入，導(dǎo)致模擬出的兩個(gè)月之后的氣象結(jié)果和原先大相徑庭！在計(jì)算化學(xué)中，有兩類任務(wù)混沌效應(yīng)很明顯，一個(gè)是分子動(dòng)力學(xué)，另一個(gè)是幾何優(yōu)化。3.3 擾動(dòng)對(duì)分子動(dòng)力學(xué)的影響分子動(dòng)力學(xué)的混沌效應(yīng)，尤其是對(duì)蛋白質(zhì)折疊過程的影響被研究得比較多，如PROTEINS,29,41。7動(dòng)力學(xué)過程中初始擾動(dòng)會(huì)隨著模擬時(shí)間指數(shù)型放大，這并不是由于數(shù)值計(jì)算問題所導(dǎo)致的，而是體系內(nèi)在性質(zhì)的體現(xiàn)。若在t=0時(shí)對(duì)體系坐標(biāo)引入微擾，則模擬到t時(shí)軌跡與未引入微擾時(shí)的差異可寫為diff（t尸A*exp（入t）其中人是依賴于體系

21、的Lyapunov指數(shù)。擾動(dòng)的放大速度相當(dāng)快，研究發(fā)現(xiàn)對(duì)蛋白質(zhì)結(jié)構(gòu)引入RMSD=1*10八-9埃的擾動(dòng)，經(jīng)過2Ps的模擬就能被放大到RMSD=1埃。假設(shè)在模擬中觀察到某時(shí)刻有一個(gè)氨基酸側(cè)鏈發(fā)生了扭轉(zhuǎn)，如果在此10Ps前引入10次RMSD=0.001的不同方式的擾動(dòng)，并分別跑10次動(dòng)力學(xué)，則只有一次能夠重現(xiàn)側(cè)鏈扭轉(zhuǎn)?？梢?，動(dòng)力學(xué)過程對(duì)初始條件十分敏感，各種因素帶來的微小數(shù)值擾動(dòng)都足矣使動(dòng)力學(xué)軌跡發(fā)生很大分歧。更何況實(shí)際模擬中擾動(dòng)并不僅僅是在初始條件中出現(xiàn)，在模擬過程中還可能會(huì)不斷納入。注意對(duì)于蛋白質(zhì)體系，擾動(dòng)不是隨模擬的進(jìn)行而無限放大的，因?yàn)榈鞍椎脑邮艿搅?chǎng)的束縛，且模擬時(shí)的溫度有限，因此勢(shì)能

22、面可及的范圍是有限的，故初始擾動(dòng)造成的差異在迅速放大后會(huì)遇到上限。對(duì)于柔性體系，由于勢(shì)能面可及范圍會(huì)加更大，可以預(yù)見差異放大的上限會(huì)比剛性體系更高。一些辦法可以提高軌跡重現(xiàn)性。提高浮點(diǎn)數(shù)精度是重要策略之一，這可以減少擾動(dòng)量的數(shù)量級(jí)。有人抱怨gromacs軌跡重現(xiàn)性不如Amber，這很大程度上是由于gromacs默認(rèn)是單精度編譯，而Amber默認(rèn)用雙精度。不過，提高精度并不解決根本問題，只不過是讓軌跡的歧化慢一些罷了，對(duì)于長(zhǎng)時(shí)間的動(dòng)力學(xué)模擬結(jié)果的偏差仍然是無法預(yù)測(cè)的，除非能讓浮點(diǎn)精度無窮高，但這顯然是不可能的。若想讓軌跡完全重現(xiàn)，通用的辦法是將第一節(jié)涉及的所有可能帶來擾動(dòng)的因素都去除。而有些程序

23、也專門考慮了軌跡重現(xiàn)問題，在gromacs里的mdrun有-reprod選項(xiàng)以幫助軌跡重現(xiàn)。另外還有人弄了個(gè)固定精度的分子動(dòng)力學(xué)方法，不僅可以完全重現(xiàn)軌跡，還可以將時(shí)間反轉(zhuǎn)從末態(tài)得到初始結(jié)構(gòu)，完全體現(xiàn)出分子動(dòng)力學(xué)的決定性特點(diǎn)。動(dòng)力學(xué)軌跡的重現(xiàn)究竟有無意義？如果只是出于計(jì)算目的，比如一段軌跡文件弄丟了，或者希望以更高的頻率輸出體系能量，那么重新跑出一次完全一樣的軌跡是有意義的。然而，從理論角度上來講，軌跡的嚴(yán)格的重現(xiàn)性是毫無意義的。動(dòng)力學(xué)過程就是在體系相空間（坐標(biāo)+動(dòng)量空間）中游歷的過程，如果跑兩次得到了截然不同的軌跡，只不過是說明在相空間不同位置區(qū)域游歷罷了，更沒有誰對(duì)誰錯(cuò)之分。對(duì)于一個(gè)平衡態(tài)

24、的模擬，體系的性質(zhì)（比如內(nèi)能、擴(kuò)散系數(shù)）就是動(dòng)力學(xué)過程時(shí)間的平均，兩次跑出來軌跡相不相同完全無所謂，只要時(shí)間足夠長(zhǎng)，它們都足以充分遍歷相空間，結(jié)果會(huì)是等價(jià)的。對(duì)于有限時(shí)間的模擬，為了加強(qiáng)某些區(qū)域的采樣，實(shí)際上還有人專門利用動(dòng)力學(xué)軌跡的混沌效應(yīng)，給初始速度分配不同的值，分別做動(dòng)力學(xué)模擬并將結(jié)果合并。而對(duì)于非平衡過程的模擬，比如說蛋白質(zhì)折疊過程，軌跡（折疊路徑）的重現(xiàn)性同樣毫無意義。因?yàn)榈鞍椎恼郫B過程本來就是混沌的，初始的微小擾動(dòng)會(huì)使折疊路徑相差迥異。在現(xiàn)實(shí)中每次蛋白的初始條件絕不可能完全一致，因此同樣的折疊路徑在現(xiàn)實(shí)中無法重現(xiàn)，或者說現(xiàn)實(shí)中沒有可能發(fā)生兩次完全相同的折疊過程，那么苛求在多次模擬中

25、重現(xiàn)相同的軌跡有何意義？實(shí)際上不僅不應(yīng)該重現(xiàn)，折疊路徑的多樣性才是真正感興趣的，也很適合通過概率統(tǒng)計(jì)來研究，這就需要做很多次動(dòng)力學(xué)模擬，以得到不同的折疊路徑并進(jìn)行分析。正好，由于各種數(shù)值因素的介入，哪怕相同的初始條件都會(huì)產(chǎn)生不同的折疊路徑，都省得修改輸入文件了，直接反復(fù)跑多次就行了。注意盡管每次的軌跡十分不同，但最終卻一定會(huì)折疊到同樣的最穩(wěn)定結(jié)構(gòu)，這并不違背混沌效應(yīng)，而是由于每一次折疊過程都會(huì)感受到相同的“吸引子”，即趨向成為最穩(wěn)結(jié)構(gòu)的一種特殊勢(shì)場(chǎng)。折疊過程可以比喻為站在一個(gè)坑坑洼洼的大坑（對(duì)應(yīng)自由能面的形狀）旁邊往坑里扔石子，每次石子掉落的軌跡都很不同，但最終都掉到坑底。這里順便談?wù)勔粋€(gè)相關(guān)

26、問題，也就是單個(gè)動(dòng)力學(xué)模擬的軌跡是否足以說明問題。上面說了，一條折疊軌跡不能反映全部的折疊過程，同樣地，研究其它過程，原則上也需要用多條軌跡來統(tǒng)計(jì)地說明問題。比如研究在結(jié)合配體后激酶的活性位點(diǎn)打開的過程，在軌跡中發(fā)現(xiàn)經(jīng)過530ps模擬后活性位點(diǎn)打開了，那么就真的能說明在現(xiàn)實(shí)中配體結(jié)合后經(jīng)過約530Ps后活性位點(diǎn)就能打開？非也，模擬10次，會(huì)得到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)致動(dòng)力學(xué)過

27、程中各種“事件”發(fā)生的時(shí)刻的波動(dòng)范圍極大，或者說每種事件在不同的時(shí)間范圍都有一定出現(xiàn)概率，只分析單一軌跡就造成了一葉障目，結(jié)論的不可靠性會(huì)較大。然而目前絕大多數(shù)模擬都是基于單一軌跡的，很少有人對(duì)同一問題做多軌跡統(tǒng)計(jì)分析，主要是多軌跡的計(jì)算量太大或者研究的內(nèi)容不需要那么深入詳細(xì)，其次是很多人沒意識(shí)到這個(gè)問題，還有少數(shù)人是因?yàn)橥稒C(jī)取巧。比如論證一個(gè)抑制劑的作用，可能模擬反復(fù)跑了7次都無法論證抑制劑管用，當(dāng)跑到第8次時(shí)慶幸地發(fā)現(xiàn)能說明抑制劑管用，于是他只用這一次的軌跡說事，另外跑的7次雖然沒成功卻都避而不談，可想而知這抑制劑的實(shí)際功效如何了.3.4擾動(dòng)對(duì)幾何優(yōu)化的影響相較分子動(dòng)力學(xué)的混沌效應(yīng)，幾何優(yōu)

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

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