給 DFT 人的 MDMD for the DFT-trained
力場不是便宜版的 DFT，它回答的是完全不同的問題。
01一句話只有一句
力場是一個帶參數的經典位能面，MD 給你的是系綜平均，不是單點能量。
02Why should I care?我為什麼要讀這頁
因為你即將把 H100 時數花在一個你以為你懂、但其實問錯問題的工具上。
從 DFT 走過來的人第一個直覺是：MD 就是便宜版的 DFT。用一個粗糙的位能面換來大系統和長時間，本質上還是在算能量、還是在找穩定結構、還是能比較誰比較穩。這個直覺在每一個字上都是錯的，而且錯的方式很貴。
這個誤解有一個很常見的具體後果：以平衡結合 ΔG 為排序主軸的變體篩選流程（不論那個 ΔG 來自 end-point 近似還是 alchemical 自由能），會系統性地漏掉專一性真正住的地方。文獻上的高保真變體 —— SpCas9-HF1、eSpCas9(1.1)、HypaCas9 —— 對 on-target 與 off-target 的結合親和力跟 WT 差不多。這不是算錯了，是問錯了：專一性住在 HNH（Cas9 兩個切割活性中心之一，會整塊翻轉去對準切位，見 Cas9 domain）構形檢查點的動力學裡，而不是住在結合自由能裡。而「該問系綜與速率，還是該問單點能量」正是這一頁的全部內容。
再講一個更直接的：你在 5F9R 上跑的是 535,265 個原子、4 條 replica、每條約 87 ns。如果你以為那些 GPU 時數買的是「精度」，你會把下一筆錢繼續花錯地方。MD 的錢買的是取樣，不是精度。 力場的系統誤差不會因為你多跑 87 ns 而變小一點點；能靠時間變小的只有統計誤差。
這裡要避開一個順手就會犯的誤讀：不是「DFT 的誤差可以靠算力消掉、MD 不行」。 DFT 一樣有兩個軸 —— 數值收斂（ENCUT、k-mesh）可以靠算力壓，泛函本身的系統誤差壓不動，加再大的 ENCUT 也只是更精確地逼近同一個泛函的答案。真正的差別是 MD 多了第三個軸：統計誤差，而且你買到的 GPU 時數幾乎全部落在這第三個軸上。三個軸各有各的解法：數值軸做收斂測試、統計軸加軌跡長度與 replica、系統誤差軸只能換力場或對外部參考驗證。把三個軸混成一個「跑久一點就好」，就會出現「用算力解決力場問題」這種燒錢方式。
03What這是什麼
一個力場就是一個寫死形式、參數 fit 出來的經典位能函數：
就這樣。沒有波函數、沒有電子密度、沒有自洽場。給定座標 R，這個式子直接吐出一個數字和一組力，成本是 O(N log N)（靠 PME 處理長程靜電）。這就是為什麼你能算五十萬個原子。
從這個式子直接讀出三個結構性的限制。它們不是精度問題，是能力邊界問題 —— 不管你怎麼調參數都跨不過去：
1. 沒有電子，所以不能斷鍵。 Σ_bonds 是諧振子，拉長只會單調上升，永遠不會解離。更關鍵的是：哪兩個原子之間有鍵、每個原子帶多少電荷、每個原子屬於哪個 atom type —— 這些在建系統那一刻就凍結了，整條軌跡都不會變。你不能有 charge transfer，不能改氧化態，不能有質子轉移，不能有共價中間態。
2. 極化不是被算出來的，是被「平均進去」的。 主流的蛋白／核酸力場（Amber、CHARMM）都是 fixed-charge、non-polarizable：每個原子的部分電荷是一個常數。真實分子在水裡會被鄰居極化，力場處理這件事的方式是在 fit 電荷的時候就故意讓它「已經被極化過」—— Amber 系列用 HF/6-31G* 的 RESP 電荷，正是因為這個組合會系統性地過度極化，剛好粗略模擬水環境的平均效果。
它跟 DFT+U 的精神是同一類的 —— 用一個 fit 出來的、狀態無關的參數，去替代一個你負擔不起的自由度。好處是便宜到荒謬，代價是它在「環境跟 fit 時不一樣」的地方會系統性地錯，而且錯的方向你事先不知道。
3. 參數是 transferable 的。 同一組 sp3 碳的參數用在整個蛋白的每一個 sp3 碳上。這是力場便宜的根源，也是它誤差的根源。
對照表：你在 DFT 裡的動作，在 MD 裡對應到什麼
| DFT 世界 | MD 世界的對應 | 真正的差別 |
|---|---|---|
| ENCUT / basis set 收斂 | 沒有對應 | 力場沒有「加大就變準」的旋鈕。它的誤差是 fit 的系統誤差，不能靠數值收斂消除 |
| k-mesh 收斂 | PME grid、real-space cutoff、Ewald tolerance | 這個是真的數值參數，可以也應該收斂測試 |
| SCF 收斂 | 無 | 沒有電子自由度要自洽；每步只是代入公式算力 |
| 幾何最佳化到 local minimum | minimization | 目的完全不同：MD 的 minimization 只是移除原子重疊，不是找極小點 |
| XC functional 的選擇 | 力場的選擇（ff14SB / ff19SB / CHARMM36m） | 兩者都是「同一個系統換一套近似就換一組答案」，都需要對外部參考驗證 |
| 同一個結構永遠給同一個能量 | 同一個起始結構、不同亂數種子 → 不同軌跡 | MD 有 Lyapunov 不穩定性，幾十 ps 內兩條軌跡就分道揚鑣。這不是 bug，是本質 |
| 答案是 0 K 的一個數字 | 答案是 300 K 系綜的一個分布 | 熵、溫度、溶劑全在裡面，代價是每個數字都帶統計誤差棒 |
04Why為什麼重要
因為 MD 真正在算的是一個機率分布，不是一個能量。
平衡態 MD 的目的是對 Boltzmann 分布取樣：
而 MD 用時間平均代替上面這個系綜平均（這就是 ergodic hypothesis，也是整個方法唯一的信仰跳躍）。所以每一個你從軌跡讀出來的數字 —— 接觸數、RMSF、氫鍵壽命、界面埋藏面積 —— 都是一個 ⟨A⟩ 的估計值，都應該附一個誤差棒。你的專案已經把這件事量化了：BSA 跨 4 條 replica 的 CV 是 2.1%，逐位點距離是 38 到 46%。同一條軌跡吐出來的兩個數字，一個可以拿來打分，一個只是雜訊。
能量的意義完全換了一個層次：
| 你要的量 | DFT 怎麼給 | MD 怎麼給 |
|---|---|---|
| 內能 | 直接就是 E(R) | ⟨U⟩，但這個數字本身幾乎沒有用（在五十萬原子的系統裡它是 10⁶ kcal/mol 量級，其中絕大部分是水，瞬時漲落遠大於你關心的差異） |
| 熵 | 事後補：算 phonon，諧近似積分 | 內建。你在取樣構形空間，構形多的狀態自動有高熵，不用額外算 |
| 自由能 | G = E + ZPE + 熱修正，一層層加 | 只能從分布回推：F(x) = −kBT ln P(x) + C，或用 alchemical／umbrella 之類的方法 |
| 結合能 | Eads = Eslab+ads − Eslab − Eads，乾淨 | 沒有這麼乾淨的東西。end-point 近似（MM-PBSA）方差極大，嚴謹做法是 alchemical FEP，貴得多 |
| 速率、時間相關 | 給不了 | 這是 MD 獨有的：residence time、correlation function、構形轉換路徑 |
最後一列是重點。你的專案的核心論點 —— 專一性由 kinetics 而非平衡 ΔG 決定 —— 在你手上的計算方法裡，只有 MD 這一族碰得到。DFT 給不了，序列模型給不了，結合親和力預測器也給不了（實驗端有單分子 FRET 這條路，但那不是你在做的事）。這就是為什麼你在跑 MD，也是為什麼把 MD 拿去算結合能，等於買了一台顯微鏡然後拿來當紙鎮。
DFT 回答「在這個固定幾何下，電子怎麼排、能量多少」。MD 回答「這個系統在 300 K 的水裡會擺出哪些構形、各佔多少比例、彼此之間切換多快」。這是兩個幾乎不重疊的問題集。把 MD 當便宜 DFT 用，會得到「便宜且錯」的答案。
05How怎麼做
每一步在做什麼： 算 U(R) 的梯度得到力，用 Verlet 類的積分器推進一個 timestep，重複幾千萬次。成本幾乎全部在 nonbonded 那兩項（LJ 與 Coulomb）。所以決定你速度的是原子數（主要是水的數量）和 cutoff，而不是「系統多複雜」。
該做收斂測試的參數（這些是真的數值參數）：
- real-space cutoff：1.0 nm 是常見值，你的專案就用這個。LJ 的截斷要不要加 dispersion correction 要說明。
- PME grid spacing / Ewald error tolerance：這相當於你的 k-mesh。預設通常夠，但如果你要報靜電相關的量，值得測一次。
- timestep：這不是「調快一點」，是動力學正確性的門檻，而且門檻跟積分器綁在一起。OpenMM 使用手冊的說法是：在
HBondsconstraint 之下，Verlet 類約 2 fs、Langevin 類約 4 fs；要再往上才需要靠 HMR 把氫加重。你的專案是LangevinMiddleIntegrator+HBonds+hydrogenMass=1.5 amu跑 4 fs。細節見 Production 設定。 - 軌跡長度與 replica 數：這是統計誤差軸，見 取樣與收斂。
不能靠收斂測試解決的（這是關鍵的認知差異）： 力場本身。你把 cutoff 加到 1.4 nm、PME grid 加密一倍，得到的是「同一個力場的更精確答案」，不是「更接近真實的答案」。要檢驗力場，只有三條路：對實驗量測比對、對 QM 比對、換一個獨立的力場重跑看結論會不會翻。沒有第四條路。
什麼時候必須回到 QM 或 QM/MM：
- 切 DNA 的化學步：磷酸二酯鍵的斷裂、RuvC 與 HNH 的 two-metal-ion 機制。這是共價鍵重組，力場結構上不能碰。
- 活性中心 Mg²⁺ 的配位化學：二價陽離子的極化與電荷轉移在 fixed-charge 力場裡被系統性描述錯，配位水的交換速率可以錯好幾個數量級。
- pKa 位移與質子轉移：protonation state 在 MD 裡是建系統時凍結的常數。要讓它會動，得用 constant-pH MD 或 QM/MM。
- 金屬中心的自旋態、電荷轉移、光激發：完全在力場的語言之外。
什麼時候留在 MD 就好： 構形、界面接觸、柔軟度、溶劑可及性、變體造成的構形分布位移 —— 也就是你的 workflow 目前 90% 的內容。
成本的量級感： 5F9R 完整三元複合體是 535,265 個原子，在 H100 上跑到約 87 ns 是一個 replica 的量級。同一筆錢在 DFT 裡，你連幾百個原子的單點能量都要斟酌。這個對比不是為了說 MD 比較好，是為了讓你知道兩者的錢花在完全不同的地方：DFT 的錢買精度，MD 的錢買取樣。
06When什麼時候用
- 你要問的是構形問題：這個接觸在 300 K 的水裡撐不撐得住？這個 loop 會不會塌？這個界面實際上有多少接觸是持久的（而不是晶體裡那一張快照恰好有的）？
- 你要比較 WT 與變體的構形分布，而且你接受這個比較必須建立在多條 replica 的聚合量上。
- 你要產生後續分析或 end-point 計算需要的構形集合（PCA、MSM、MM-PBSA 都是吃軌跡的下游）。
- 你關心的過程時間尺度落在 ns 到 μs，而且不涉及成鍵。
- 你需要一個能同時容納溫度、溶劑、離子、熵的描述，而且願意用力場的系統誤差換這些。
07When NOT什麼時候別用
-
任何涉及鍵斷裂或形成的問題，不要用。 為什麼會壞：
Σ_bonds是諧振子，鍵拉長只會爬能量、永遠不解離；topology 在建系統時就凍結了。怎麼看出它壞了：這個錯誤不會給你錯誤訊息，它會給你一條完美的軌跡和一個「這個鍵非常穩定」的假結論。判斷方式不是看輸出，是看問題本身 —— 只要你要回答的東西需要 bond order 改變，力場就結構上不能回答，這跟你跑多久、參數多好都無關。 -
活性中心 Mg²⁺ 的配位化學，不要用 fixed-charge 力場的結論。 為什麼會壞：二價陽離子周圍的極化與部分電荷轉移貢獻很大，固定點電荷把它們全部丟掉，結果是 Mg²⁺ 與水／磷酸的靜電吸引被系統性高估。怎麼看出它壞了：檢查第一水合層的水交換 —— 如果整條 87 ns 軌跡裡一個水分子都沒換過、配位數是一條完美的直線，那不是「Mg 配位很穩定」這個生物學結論，那是力場把它焊死了。
-
用單一 snapshot 或短軌跡的絕對能量做排序，不要。 為什麼會壞：在蛋白–DNA–RNA 這種高電荷系統裡，end-point 方法分解出來的靜電交互作用項是 10³ kcal/mol 量級，而且必須跟一個同量級、反號的去溶劑化項相消，而你想分辨的變體差異可能只有 1 到 2 kcal/mol。兩個大數相減的統計誤差比訊號大得多。怎麼看出它壞了：把同一個變體用不同 replica 各算一次，比較跨 replica 的 SD 與你宣稱的 WT–變體差異。如果 SD 比差異大，你報的是雜訊。 這正是拿 MM-PBSA 當排序軸時的結構性風險。
-
需要 μs 到 ms 的過程，不要用標準 MD 直接看。 HNH 完整翻轉並 docking、R-loop（sgRNA 把雙股 DNA 頂開、自己去配對 target 股後形成的三股結構）完整形成都在這個尺度。為什麼會壞：這是所有錯誤裡最陰險的一種 —— 你會看到「它沒有發生」，然後把它讀成「它不會發生」。前者是取樣不足，後者是物理結論，兩者長得一模一樣。怎麼看出它壞了：檢查你的 observable 在整條軌跡裡有沒有跨過任何一個 barrier。如果它從頭到尾待在起始 basin，你量到的是那個 basin 的寬度，不是自由能地景。要碰這個尺度得用 enhanced sampling（GaMD、metadynamics、umbrella），成本與詮釋難度都跳一級。
-
起始結構本身可疑時，不要用 MD 來「修」它。 為什麼會壞：力場只會讓一個錯的起點在 300 K 下抖動，它沒有任何機制把結構拉回正確的 basin —— 因為正確的 basin 隔著一個遠高於 kBT 的 barrier。給你一把尺：300 K 下 kBT ≈ 0.6 kcal/mol，用過渡態理論回推，要在 87 ns 裡跨過去一次，barrier 大約不能超過 8 kcal/mol；要跨到能取樣兩邊的族群，門檻還要更低。構形重排的 barrier 通常遠在這之上。你的專案有實測：用 Chai-1 預測結構起始時，逐殘基 top-8 接觸殘基跨 replica 是 0/8 共享、ρ ≈ −0.03（等同隨機）；換成 5F9R 晶體起始後變成 3/8 共享、pairwise Jaccard 0.48、Spearman +0.77。同一套 metric、同一套腳本。元兇是起始結構，不是 metric。 怎麼看出它壞了：跨 replica 一致性崩掉，而且是從逐殘基層級先崩。
-
不要跨力場、跨水模型比較絕對數值。 為什麼會壞：力場是一整套 fit 出來的，換掉任何一個零件，所有絕對值都會系統性平移。具體一點：TIP3P 的自擴散係數是實驗值的兩倍以上，剪切黏度只有實驗值的三分之一上下。也就是說 TIP3P 系統裡的「時間」是被系統性加速的。（這兩個比值是後續 MD 文獻的結論，不是 Jorgensen 1983 本身 —— 那篇是 Monte Carlo，方法上量不到動力學。確切數值隨 cutoff、PME 設定與有限尺寸修正而異，要寫進 limitation 就自己引一次原始量測，別引錯來源。）怎麼看出它壞了：任何把 TIP3P 軌跡的 ns 數直接對應到實驗速率的敘述都是錯的。能比的只有同一套設定下的相對排序。
「跑久一點應該就會準一點。」這個直覺來自 DFT：加大 ENCUT、加密 k-mesh，答案確實會更逼近該泛函的正確值，而且那條路你走過很多次。但在 MD 裡，時間買的是統計誤差，買不到力場的系統誤差（就像加大 ENCUT 也修不了泛函選錯）。如果你的問題出在力場，跑一年也是同樣的錯誤答案，只是誤差棒比較小。
08Project Lens跟我的 project 多相關
為什麼是五星： 因為這一頁不寫清楚，整張地圖上其他每一頁都會被讀錯。
你的專案的核心爭議 —— kinetics vs equilibrium —— 本質上根本不是 Cas9 生物學問題，是「MD 到底在回答什麼」的方法學問題。以 MM-PBSA 或 FEP 排序平衡結合 ΔG 的篩選流程，在技術上沒有錯；錯在把一個動力學的專一性問題交給一個平衡態的 observable。而高保真變體（SpCas9-HF1、eSpCas9(1.1)、HypaCas9）對 on/off-target 的結合親和力與 WT 差不多，這件事就是那個結構性盲區的直接證據。要把這件事講清楚，需要的不是更多 Cas9 文獻，是這一頁。
為什麼不是更低： 因為它零 GPU 時數卻能省下大量 GPU 時數，而且它管的是「這個 run 該不該投」，比 系統建立 管的「投出去要怎麼建」還要早一層。省錢的順序就是這個順序。
一個很常見、也很貴的浪費是這樣發生的：eSpCas9(1.1) 在 on-target 底物上是 null —— 而它本來就該 null，因為這個變體削弱的是對 non-target 股的靜電抓握，完全配對的底物上沒有東西可削，要看到專一性效應必須換帶 mismatch 的底物。這種 run 本身不是廢的（高保真變體在 on-target 上跟 WT 沒差，正好是一個乾淨的陰性對照，證明流程不會憑空生出差異），問題在於它是當成「找訊號」的 run 投出去的，而它在底物選擇上就結構性地給不出那個訊號。這不是模擬技術問題，是「這個 metric 在量什麼、在什麼條件下才會有訊號」沒有在按下 sbatch 之前寫下來。在「便宜」是硬約束的現實下，這種期待落差是負擔不起的。
具體怎麼用： 在 workflow 的每一層前面，強制問三個問題：
- 我要的量是系綜平均、時間相關量，還是單一構形的性質？（第三種通常不該用 MD 拿。）
- 這個量的時間尺度在不在標準 MD 碰得到的範圍？如果不在，我是要換 enhanced sampling，還是換一個代理量？
- 這個量跨 replica 穩不穩？你的 reliability 表已經先回答了：BSA 2.1%、總 protein–DNA 接觸 4.1%、target-DNA 接觸 9.4%、鹽橋 9.5% 可以打分；逐位點距離 38 到 46%、DNA RMSF 59.1%、groove 接觸計數 145.5% 不行。
09Reviewer Thinkingreviewer 會問什麼
你用 non-polarizable force field 模擬一個高電荷的 protein–DNA–RNA 界面，靜電項可信嗎？點開看參考答案
誠實的回答是：絕對值不可信，相對趨勢在同一套參數下可用。 蛋白–核酸界面是這類力場最吃力的場景之一，因為極化效應在高電荷密度處最大。可以防守的做法有三層：第一，不要報絕對結合能，只報接觸與幾何類的量；第二，報跨 replica 的統計而不是單點；第三，如果非要報能量，至少換一個獨立的力場重跑一組對照，證明結論的方向不變。打分層用聚合接觸量而非能量，正好繞開了這一問的正面攻擊。
你怎麼確定你看到的差異不是 force field artifact？點開看參考答案
唯一站得住的回答是外部參考：跟實驗量測比、跟 QM 比、或跟另一套獨立力場比。你沒有辦法用同一套力場的更多模擬來證明那套力場是對的。 實務上最便宜的一招是找一個已知答案的陽性對照（例如文獻上機制清楚的高保真變體），確認你的流程能把它排出來；排不出來就先別談新變體。這也是你的 ground truth panel 存在的理由。
你的結論需要多長的時間尺度？87 ns 夠嗎？點開看參考答案
關鍵是分開回答「你量的是什麼」。界面接觸、埋藏面積、鹽橋這類量在幾十 ns 就能取樣到合理的平均值，而你的 CV 表（BSA 2.1%、總接觸 4.1%）就是這件事的證據。但 HNH 的完整構形轉換不在這個尺度上：Palermo 等人累積約 15 μs 的 GaMD 才碰得到活化路徑，Nierzwicki 等人談 allosteric 訊號用的是 3 replicas × 約 1.2 μs。把數字擺齊：你的單條 87 ns 對 Nierzwicki 的單條 1.2 μs 差約一個數量級，你的 4 × 87 ns 總量對 Palermo 的約 15 μs 差約一個半數量級 —— 而且 GaMD 帶 boost potential，等效取樣尺度還要再拉開。差距落在一到兩個數量級之間，講的時候不要誇大也不要縮小。 結論是：任何從 87 ns 軌跡宣稱「HNH 沒有 docking，所以這個變體不會活化」的敘述都站不住。同一條軌跡對不同 observable 的收斂程度完全不同，方法段要一個一個講。
為什麼不做 QM/MM？點開看參考答案
因為你要回答的不是化學步。QM/MM 是切鍵那一步的工具，而你的專案問的是「HNH 有沒有走到能切的位置」—— 這是構形問題，在 QM/MM 之前。誠實的補充是：如果最後發現變體的差別出現在催化步而非 engagement 步，那 QM/MM 就變成必要的，而那會是另一個等級的計算成本。先講清楚你的 scope 到哪裡結束，比假裝什麼都涵蓋了要安全得多。
你的 replica 是怎麼獨立的？只換亂數種子算不算獨立取樣？點開看參考答案
只換初速度種子給的是「同一個建系統決定下的統計獨立」。因為 MD 有 Lyapunov 不穩定性，幾十 ps 後兩條軌跡就完全分開，所以它確實是有效的統計獨立樣本 —— 但它沒有涵蓋建系統本身的不確定性（loop 怎麼建、水怎麼放、離子放哪）。要涵蓋後者，得從不同的 build 重來。這一點值得主動在方法段寫明，而不是等 reviewer 問。
10Common Mistakes最常犯的錯
- 把 MM-PBSA 的輸出當成 ΔG_bind 報。 後果：報了一個方差比訊號大、而且系統性偏差未知的數字，reviewer 一問誤差棒就垮。正解：報跨 replica 的分布，明說這是 end-point 近似，並且不要用它做最終排序。
- 方法段只寫「用了 4 fs timestep」。 後果：4 fs 能不能撐得住，取決於積分器與 constraint，只寫步長等於沒寫；constraint 漏掉時高頻 X–H 振動被積分器解不動，能量漂移、結構慢慢變形，而且不一定會崩給你看。正解：底線是
HBondsconstraint，而且 OpenMM 手冊給的門檻是「HBonds 之下 Verlet 約 2 fs、Langevin 約 4 fs」；HMR 是在這之上多買的餘裕，不是省略 constraint 的替代品。積分器、constraint、hydrogen mass、步長四件事一起寫。 - 拿軌跡最後一幀當「模擬得到的結構」。 後果：那只是系綜裡隨機的一個樣本，跟平均構形沒有必然關係。正解：報群聚（clustering）後的代表結構，或直接報分布。
- 把 minimization 當成 DFT 的幾何最佳化，跑到力收斂為止。 後果：浪費時間，而且給人一種「已經找到穩定結構」的錯覺。正解：minimization 只是為了消除原子重疊，跑幾百到幾千步就夠，馬上進入受約束的分階段平衡。
- 跨 paper 比較絕對數值。 後果：力場、水模型、cutoff、長度全都不同，數字不可比。正解：只比同一套設定下的相對排序，要跨組比就自己重跑對照。
- 期待 MD 修好一個爛的起始結構。 後果：得到一條看起來很正常、結論卻是隨機的軌跡（你的 Chai-1 起始 0/8、ρ ≈ −0.03 就是實例）。正解：把起始結構當成獨立的品質關卡，見 系統建立 與 結構預測。
11Further Reading讀哪幾篇
12Summary三句話
力場是一個沒有電子、鍵不能斷、極化被平均進參數裡的經典位能面，它的誤差來源從基組與泛函換成了參數化與取樣。MD 的答案是 300 K 系綜的分布與時間相關量，不是單點能量，所以每個數字都要帶跨 replica 的誤差棒，而多跑時間只縮小統計誤差、不縮小力場的系統誤差。凡是需要成鍵、金屬配位化學或質子轉移的問題，都該回到 QM 或 QM/MM，而不是期待更長的軌跡把它救回來。