Level 3 — 知道 MD 在模擬什麼What MD actually does
從 DFT 的世界跨到力場的世界,知道換到了什麼、失去了什麼。
你現在在哪
L0 到 L2 全部發生在生物與問題定義層,你一行程式都還沒跑。L3 是你第一次換世界:從 0 K 的電子結構,換到 300 K、有水、非諧的系綜。
這一關對你來說是最危險的一關,理由很違反直覺:你太熟悉計算了。 從 DFT 走過來的人第一個直覺是「MD 就是便宜版的 DFT」—— 用一個粗糙的位能面換來大系統與長時間,本質上還是在算能量、找穩定結構、比誰比較穩。這個直覺在每一個字上都是錯的,而且錯得很貴。
最乾淨的拆法是三個誤差軸:
- 數值軸 —— DFT 有(ENCUT、k-mesh),MD 也有(PME grid、real-space cutoff、Ewald tolerance)。兩邊都能靠算力壓到收斂。
- 系統誤差軸 —— DFT 是泛函,MD 是力場。兩邊都壓不動,加再多算力也只是更精確地逼近同一個近似的答案。
- 統計誤差軸 —— DFT 沒有,MD 有,而且你買的 GPU 時數幾乎全部落在這一軸上。
把三個軸混成一句「跑久一點就好」,就會出現「用算力解決力場問題」這種燒錢方式。
這一關結束之後,你能做兩件你現在做不到的事: 讀一篇 MD 論文的方法段,判斷它宣稱的東西撐不撐得住;以及在把計算送進佇列之前,自己列得出從一顆 PDB 到能跑的系統之間那十幾個決定,並指出哪一個最不可逆。
這關要學會什麼
- 你將能夠寫出力場位能函數的形式,並直接從式子讀出三個結構性限制:不能斷鍵、極化被平均進參數、參數是 transferable 的。這些不是精度問題,是能力邊界。
- 你將能夠把 DFT 世界的每一個動作對應到 MD(或明確指出「沒有對應」),特別是「ENCUT/basis set 收斂在 MD 裡沒有對應物」這一條。
- 你將能夠分辨哪些參數該做收斂測試(PME grid、cutoff、timestep),哪些做了也沒用(力場本身),以及檢驗力場的三條路(對實驗、對 QM、換獨立力場重跑)。
- 你將能夠說出起始結構為什麼是整條流程裡最不可逆的決定,並用 0/8 對 3/8 這組對照佐證。
- 你將能夠估算:在 300 K、87 ns 的軌跡裡,大約多高的 barrier 還跨得過去 —— 並用這把尺判斷某個構形轉變在你的尺度上有沒有機會發生。
- 你將能夠對每一個 observable 分別回答「它在我的取樣量下收不收斂」,而不是籠統地講「我的模擬收斂了」。
依序讀這些
順序的邏輯是:世界觀 → 起點 → 起點的替代品 → 建系統 → 跑 → 判斷夠不夠。 特別注意 cas9-structures 排在 system-setup 前面,這是刻意的 —— 起始結構是這條流程裡最貴、最不可逆的一個決定,把它排在「建系統的技術細節」後面,等於默許你先挑好水盒再想要模擬哪個狀態。
- 給 DFT 人的 MD — 換世界。這一頁不讀通,後面五頁全部會被讀成「DFT 的操作手冊」。
- 該用哪個結構起頭 — 你選哪個 PDB,決定的不是初始條件,是你在模擬哪一個反應步驟。
- 結構預測能不能用 — 承接上一頁:AlphaFold/Chai-1 到底能放在流程的哪一層。答案不是「不能用」,而是「不能坐在 production 的上游」。
- 系統建立 — 十幾個決定,每一個都能安靜地毀掉結論。
- Production 設定 — 每個參數都是一次換匯:物理正確性換牆鐘時間。
- 取樣與收斂 — 放最後,因為它是判斷前面五頁夠不夠的那把尺。先有東西可以量,尺才有意義。
三篇原始文獻(讀法比內容重要)
通關檢查表
通關測驗
你把 real-space cutoff 從 1.0 nm 加到 1.4 nm、PME grid 加密一倍,結果變了。這代表什麼?沒變又代表什麼?點開看參考答案
這兩個都是數值軸的參數,也是 MD 裡少數真的對應到你熟悉的 k-mesh 收斂測試的東西。
變了 → 你原本的設定沒有數值收斂,該收斂到平坦為止。 沒變 → 數值軸收斂了。
但兩種情況都完全沒有告訴你力場是對的。這是這一題真正的陷阱:把 cutoff 加大得到的是「同一個力場的更精確答案」,不是「更接近真實的答案」。要檢驗力場只有三條路 —— 對實驗量測比對、對 QM 比對、換一個獨立力場重跑看結論會不會翻。沒有第四條路,特別是「跑久一點」不是。
你的 4 條 87 ns 軌跡裡,HNH 從頭到尾沒有 docking。可以下什麼結論?點開看參考答案
只能下一個:在這個取樣量下,我們沒有觀察到 docking。
不能下的是「這個變體不會活化」。理由是量級的問題,把數字擺齊就很清楚:Palermo 等人累積約 15 μs 的 GaMD 才碰得到活化路徑,Nierzwicki 等人談 allosteric 訊號用的是 3 replicas × 約 1.2 μs(系統約 340,000 atoms)。你的單條 87 ns 對 Nierzwicki 的單條 1.2 μs 差約一個數量級,你的 4 × 87 ns 總量對 Palermo 的約 15 μs 差約一個半數量級 —— 而且 GaMD 帶 boost potential,等效尺度還要再拉開。差距落在一到兩個數量級之間,講的時候不要誇大也不要縮小。
這是所有錯誤裡最陰險的一種: 「它沒有發生」與「它不會發生」長得一模一樣。前者是取樣不足,後者是物理結論。判斷方式是去看你的 observable 在整條軌跡裡有沒有跨過任何一個 barrier;如果它從頭到尾待在起始 basin,你量到的是那個 basin 的寬度,不是自由能地景。
你檢查活性中心的 Mg²⁺,發現整條 87 ns 軌跡裡第一水合層一個水分子都沒換過,配位數是一條完美的直線。這是好消息嗎?點開看參考答案
不是,這是一個要換方法的訊號。
二價陽離子周圍的極化與部分電荷轉移貢獻很大,fixed-charge 力場把這些全部丟掉,結果是 Mg²⁺ 與水/磷酸的靜電吸引被系統性高估 —— 力場把它焊死了。你看到的不是「Mg 配位很穩定」這個生物學結論,你看到的是力場的一個已知失效模式。
這一題的一般化版本才是重點:古典力場的失效通常不會給你錯誤訊息,它會給你一條完美的軌跡和一個過度自信的結論。 需要成鍵、金屬配位化學、質子轉移的問題,都該回到 QM 或 QM/MM,而不是期待更長的軌跡把它救回來。
你用 Chai-1 的預測結構當 production 起點,跑完 4 條 replica 之後發現逐殘基 top-8 接觸殘基是 0/8 共享、ρ ≈ −0.03。同一套腳本換成 5F9R 晶體起始,變成 3/8 共享、Jaccard 0.48、Spearman +0.77。這說明了什麼?點開看參考答案
說明元兇是起始結構,不是 metric。metric 沒變、分析腳本沒變、力場沒變,變的只有起始座標。
背後的物理很簡單,而且可以估給你看:300 K 下 kBT ≈ 0.6 kcal/mol,用過渡態理論粗略回推,要在 87 ns 裡跨過去一次,barrier 大約不能超過 8 kcal/mol;要跨到能取樣兩邊的族群,門檻還要更低。構形重排的 barrier 通常遠在這之上。所以力場沒有任何機制把一個錯的起點拉回正確的 basin,它只會讓那個錯的起點在 300 K 下抖動。
有一個一般性教訓值得記到 L5:如果你只做了 Chai-1 那一組,你會得出「逐殘基接觸是雜訊,這個 metric 沒用」,然後把它丟掉。可靠性分析必須是二維的 —— 換種子(隨機誤差)與換起始結構(系統性誤差)要分開跑,缺一個就會把兩種完全不同的病灶誤診成同一種。
你的 replica 只是換了初速度的亂數種子。這算不算獨立取樣?點開看參考答案
算,但只算一半。
因為 MD 有 Lyapunov 不穩定性,兩條只差在初速度的軌跡在皮秒量級之後就完全分開(而你用的 Langevin 積分器本身還帶隨機噪音,換種子等於換一個噪音實現),所以換種子確實給了你統計上獨立的樣本 —— 這是最低要求,而且是有效的。
但它沒有涵蓋建系統本身的不確定性:loop 怎麼補、水怎麼放、離子放哪、質子化狀態怎麼決定。這些在所有 replica 裡是一致的,所以任何由它們造成的偏差,你的跨 replica CV 完全看不見。更好的做法是重新做溶劑化與離子放置,讓水與離子的初始構型也不同;要涵蓋建模層的不確定性,得從不同的 build 重來。
另外有一個一定要能明確否認的偷懶做法:把一條長軌跡切成四段不算四條 replica。 那些片段共享同一段歷史、同一個 basin,彼此高度相關,算出來的 CV 會系統性偏低 —— 等於自己騙自己。主動在方法段寫明這一點,比等 reviewer 問好。
用 DFT 的語言解釋一句話:「MD 的錢買的是取樣,不是精度。」點開看參考答案
加大 ENCUT、加密 k-mesh,你得到的是更精確逼近同一個泛函的答案;泛函選錯了,這些算力一分錢都救不回來。
MD 完全平行:多跑 87 ns、多跑一條 replica,你得到的是更精確逼近同一個力場的系綜平均;力場錯了,跑一年也是同樣的錯誤答案,只是誤差棒比較小。
差別在於 MD 多出來的那一軸。DFT 裡同一個結構永遠給同一個能量;MD 裡同一個起始結構、不同亂數種子會給不同軌跡,而且某些量的散度大到嚇人 —— 同一條軌跡吐出來的兩個數字,BSA 跨 4 條 replica 的 CV 是 2.1%,逐位點距離是 38 到 46%。一個可以拿來打分,一個只是雜訊。 分辨這兩者,就是 L4 與 L5 的全部工作。