取樣與收斂Sampling and convergence
你的數字是物理，還是這條軌跡剛好長這樣？
01一句話只有一句
你要證明的不是軌跡跑得夠久，而是換個亂數種子答案不會變。
02Why should I care?我為什麼要讀這頁
因為你的 methodology paper 的主張，整篇就是從這一頁長出來的。
專案已經定案的那句 headline —— 「打分要用聚合量，不要用逐殘基／逐位點」 —— 不是從理論推出來的，是從 4 條 replica 的 CV 表推出來的。BSA 2.1%、總 protein–DNA 接觸 4.1%、target-DNA 接觸 9.4%、鹽橋 9.5% 落在綠燈；protein RMSD 12.2%、non-target 接觸 17.7% 是黃燈；介面 iRMSD 23.9%、DNA RMSF 59.1%、groove 距離 48%、逐位點距離 38 到 46%、groove 接觸計數 145.5% 全在紅燈。這張可靠性分層表（完整版見指標可靠性）成立的前提，是你能回答「87 ns 夠不夠、4 條夠不夠」。
用你熟的世界講，這件事的嚴重性就很清楚了：
只跑一條軌跡就報結果，等於只測一個 k-point 就宣布能量收斂。 你在 VASP 裡絕對不會這樣做 —— 你會掃 k-mesh，看能量差落到 meV/atom 等級才停手，而且會把那張收斂表放進補充材料。
這座橋只能搭到這裡，剩下兩件事要自己接住。第一，MD 沒有 ENCUT 的對應物：力場的解析度是別人配死的，你調不動，它的誤差不會因為你跑更久而變小。第二，k-mesh 收斂是單調的，replica 不是：k 點加密，能量朝一個極限靠；多跑一條 replica，逐殘基 RMSF 不會變穩，你只是把它本來就有的散度量得更清楚（這一點在指標可靠性展開）。真正對得起來的只有一件事 —— replica 數就是你的獨立取樣點數目。而 k-mesh 收斂測試是 DFT 的文化，沒人敢不做；MD 的 replica 收斂測試不是文化，多數 Cas9 論文只跑一到兩條，而且不報散布。 這個落差就是你的機會，也是你唯一不能自己犯的錯。
03What這是什麼
「收斂」這兩個字在 MD 裡被拿來指三件完全不同的事，混在一起講是大部分錯誤的起點：
- 系統還在鬆弛嗎（equilibration） —— 你的觀測量是不是還在往某個方向單調漂。這個可以量，而且很便宜。
- 我的平均值誤差多大（statistical error） —— 就算已經在平衡態，有限取樣仍然給你一個有誤差的平均值。這個也可以量，但要先處理相關性。
- 我有沒有取樣到別的構形盆地（ergodicity） —— 分子有沒有可能在別的地方待著，只是你 87 ns 沒去。這個原則上量不到，只能反證：跑更多、跑更久、從不同地方出發，看有沒有跑出新東西。
Ergodicity 的正式敘述是「時間平均 = 系綜平均，在 t → ∞ 的極限下」。 那個極限就是問題所在。你有的是 87 ns，而 Cas9 的 HNH（負責切 target 股的那個催化 domain）要翻轉並 docking 到催化構形，文獻上是微秒以上的事。所以你不是在做遍歷取樣，你是在做一個盆地內的局部取樣，而且你必須把這句話明白寫進論文，不然就是在誤導。
有效樣本數才是你真正的 n。 相鄰的軌跡幀高度相關，所以 4350 幀不等於 4350 個獨立樣本。標準寫法是統計低效率（statistical inefficiency）：
τ 是該觀測量的自相關時間（以幀為單位）。注意 τ 是「每個觀測量各自的」，不是軌跡的性質：能量的 τ 可能是幾皮秒，一個大 loop 的開合可能是幾十奈秒。所以同一條軌跡對某些量收斂得很好，對另一些量根本沒開始 —— 這正是你那張 tier 表在講的事。
Block averaging 是不用算 τ 就能拿到誠實誤差棒的辦法。 把軌跡切成長度 b 的塊，每塊算一個平均，再算這些塊平均的標準誤。b 小的時候塊之間仍然相關，標準誤被低估；b 加大到超過自相關時間，標準誤會爬升到一個 plateau。那個 plateau 值才是真的誤差。如果曲線一直爬、看不到 plateau，那就是你的軌跡短於自相關時間 —— 這本身就是「還沒收斂」的證據。
跨 replica 的 CV 是這一切之上的最終仲裁。 定義是標準差除以平均，用百分比表示。它的優點是不需要任何關於相關性的假設 —— 每條 replica 是一個獨立樣本，就這麼簡單。缺點是樣本數很小（你有 4 個），所以 CV 本身的估計也很粗。
單條軌跡的內部誤差棒會系統性低估真實不確定度，因為它看不見自己沒去過的地方。跨獨立 replica 的散布才是誠實的誤差棒。 你所有的收斂宣稱，最終都要落回這一句。
04Why為什麼重要
因為 MD 的所有輸出都是統計量，而統計量沒有誤差棒就不是數據。
更具體的理由是：取樣不足的失敗模式長得跟真實訊號一模一樣。 一條軌跡裡某個殘基特別柔軟、某個接觸特別穩定、某個距離特別短 —— 這些都可以是物理，也都可以是那一條軌跡剛好落進的一個角落。兩者在單條軌跡裡沒有任何辦法分辨。唯一的分辨方式是重複。
你的專案有一個直接的實證，而且是最貴的那種：
在 2 條 replica 的時候，逐殘基 top-8 接觸殘基共享 6/8，看起來漂亮到可以直接寫進打分軸。跑到 4 條，掉到 3/8（只剩 450、695、765）。
這不是統計波動 —— 而且機制要講對，因為 reviewer 一定會指出來：「所有 replica 都共享」是一個單調的統計量，交集只會愈交愈小，多加一條 replica 在數學上永遠不可能讓它變大。所以 n=2 的 6/8 天生就比 n=4 的 3/8 好看，這一半是定義造成的，不是發現。真正的問題在另一半：只有兩條的時候，你沒有第三條可以否證它，於是分不出「真的一致」與「碰巧一致」，而偏誤的方向是單向的 —— 小 n 幾乎只會讓你太樂觀，不會讓你太悲觀。
同一批數據還告訴你這件事的邊界在哪，這一點要一起講才誠實：逐殘基的排序並不是純亂數。 四條 replica 兩兩之間的 top-8 Jaccard 平均 0.48、接觸佔據率的 pairwise Spearman 平均 +0.77 —— 也就是「哪些殘基比較常接觸」是有一致傾向的，垮掉的是**「四條一致同意的前 8 名」這種嚴格共識**。對打分軸來說這個區別無關緊要（你要的正是嚴格共識），但寫進論文時把兩者混為一談，是會被抓的。
6/8 → 3/8。 這是同一個系統、同一套設定、只是多跑了兩條 replica。如果當時停在 n=2，你會把一份撐不住四條一致檢驗的清單寫進打分軸，然後帶著它跑完下游所有的變體排序。你日後最可能自己騙自己的地方就在這裡，因為 n=2 的結果總是先到，而且因為交集是單調的，它在數學上不可能比較差。
05How怎麼做
五個檢查，全部是既有軌跡的後處理，零額外 GPU 時數。
1. 先切掉平衡段。 目視畫你要報的量對時間的曲線，找出它停止單調漂的位置。要更客觀就用 pymbar.timeseries.detect_equilibration，它會回傳三個東西：t₀（production 起點）、g（統計低效率）、Neff。它的判準是選一個 t₀ 使得剩下那段的有效樣本數最大 —— 直觀而且不需要對分布做假設。丟太少會被鬆弛污染，丟太多會浪費昂貴的取樣，這個方法直接幫你找那個折衷點。
2. 前後半漂移。這是你專案採用的主檢查，也是最容易解釋給 reviewer 聽的一個。 把 production 切成前半後半，各算一次你要報的量，看差幾個百分比。你的實測結果是：聚合量的前後半漂移落在 2 到 17%，判定為大致收斂；rep1、rep3、rep4 走平（plateau），rep2 尾段還在微升。 注意最後這一句 —— rep2 沒有走平這件事必須報出來，不能因為「四條裡有三條漂亮」就把它省略。
3. Block averaging。 對每一條 replica 畫「標準誤 vs 塊長」，找 plateau。看到 plateau 才有資格報單條軌跡的誤差棒；看不到 plateau，就只能報跨 replica 散布。
4. 自相關時間與有效樣本數。 對你要報的每一個量各算一次 τ，換成 Neff。這裡有個殘酷的算術值得先做一次：如果某個量的自相關時間是 10 ns，那 87 ns 只給你不到 10 個獨立樣本；如果是 30 ns，你手上只有兩三個。 這就是為什麼跨 replica 的散布常常大於單條軌跡自己算出來的誤差棒 —— 後者假設了它其實沒有的獨立性。
5. 跨 replica CV，最終仲裁。 你專案用的分層是：
| CV（4 replica） | 判定 | 你的實測落在這一層的量 |
|---|---|---|
| < 10% | reproducible（可進打分軸） | BSA 2.1%、總 protein–DNA 接觸 4.1%、target-DNA 接觸 9.4%、鹽橋 9.5% |
| 10 到 20% | marginal（只能定性） | protein RMSD 12.2%、non-target 接觸 17.7% |
| > 20% | noise（不可用） | 介面 iRMSD 23.9%、逐位點距離 38 到 46%、groove 距離 48%、DNA RMSF 59.1%、groove 接觸計數 145.5% |
這三條界線是你自己定的操作性門檻，不是文獻標準。 這件事一定要在方法段說清楚，並且說明為什麼是 10% 而不是 15% —— 最誠實的理由是它跟下游決策掛鉤：只有綠燈的量才進最後的打分軸。
多條短 vs 一條長，怎麼選：
- 要探索構形空間 → 多條短贏。 四組不同的初速與噪音會走到四個不同的地方，一條長軌跡只會在它落進的那個盆地裡待著。而且多條短可以平行：在共享叢集上「同時可跑的 job 數與單次 walltime 都有上限」的現實下，四條 87 ns 分兩波跑完；同樣的總長度換成一條，你一樣得靠 checkpoint 串好幾段，卻只換到一個盆地。
- 要單一盆地內的平均值精度 → 一條長也可以，前提是它遠長於該量的自相關時間。
- 要把多條短縫成完整的動力學圖像 → 這就是 Markov state model 的正規做法，把短軌跡當成轉移矩陣的採樣，估出比任何單條軌跡都長的時間尺度。這是你目前還沒走、但值得知道存在的那條路。
實務上怎麼讓 replica 真的獨立： production 腳本把同一個 SEED 同時餵給 integrator、MonteCarloBarostat、以及 setVelocitiesToTemperature，然後從同一個 minimize 過的座標各自重新加熱平衡。要講精確：光換 integrator 種子就已經是不同的 Langevin 噪音實現，軌跡會分開 —— 三個一起換不是為了物理，是為了方法段能一句話講完、不留爭議的縫。真正沒解決的是更上游那一層：四條共用同一個起始座標與同一次溶劑化。更嚴格的做法是連溶劑化與離子放置都重做，見指標可靠性。
06When什麼時候用
- 每一條 production 軌跡跑完，第一件事就做前後半漂移檢查。 這是最便宜、最容易解釋的收斂證據。
- 任何要進論文的數字，都先過跨 replica CV。 不論它看起來多漂亮。
- 要決定「該再跑一條，還是該跑更長」的時候。 若跨 replica 散布 ≫ 單條內部誤差棒，你缺的是 replica；若單條的 block averaging 看不到 plateau，你缺的是長度。
- 要在方法段寫「我們的模擬收斂了」之前。 這句話沒有配套的量化證據就不要寫。
- 審別人的 Cas9 MD 論文的時候。 先翻方法段找 n 和軌跡長度；找不到就知道那篇的結論該打幾折。
- 要向共同作者或審稿人說明為什麼打分只用聚合量的時候。 這一頁就是那個論證本身。
07When NOT什麼時候別用
-
不要用「RMSD 走平了」當收斂證據。 為什麼會壞：RMSD plateau 是必要條件不是充分條件。RMSD 把整個結構壓成一個數字，對正交的自由度完全不敏感 —— 系統可以在同一個 RMSD 值上換過好幾個構形，曲線卻是一條平線。怎麼看出它壞了：把 RMSD 曲線跟你真正要報的那個量的曲線疊在一起看。你的數據就是現成的反例：protein RMSD 的跨 replica CV 是 12.2%（看起來還好），但同一批軌跡的 DNA RMSF CV 是 59.1%。RMSD 平不代表別的東西也平。
-
不要把一條長軌跡切成 4 段當成 4 條 replica。 為什麼會壞：它們共享同一個起始構形、同一次平衡的記憶、同一條相空間路徑。段與段之間的差異只反映該量的自相關，完全不反映「從別的地方出發會不會得到別的答案」，因此系統性低估散布。怎麼看出它壞了：段間 CV 明顯小於獨立 replica 的 CV。你的四條是各自從
state_min用不同種子重新加熱平衡的，這一點要在方法段寫清楚，因為它決定了你的 CV 在量什麼。 -
n = 2 的時候，不要下任何逐殘基或逐位點的結論。 為什麼會壞：兩個樣本的重疊率是一個變異極大的統計量，兩條碰巧同意的機率不低，而你沒有第三條可以否證它。更根本的是，共識型統計量（交集、共享數）隨 n 增加只會下降，所以 n=2 量到的值是一個上界，不是估計值。怎麼看出它壞了：加第三條、第四條，看數字往哪邊走。你的實測是 6/8 掉到 3/8。唯一安全的用法是把 n=2 的結果當 pilot，不進表、不進 abstract、不進任何對外簡報。
-
收斂宣稱不要跨觀測量套用。 為什麼會壞：自相關時間是每個量各自的。同一批軌跡的 BSA 跨 replica 只差 2.1%，groove 接觸計數卻是 145.5%。寫「我們的模擬已收斂」而不指明是哪個量，等於什麼都沒說。（後者要小心讀：groove 接觸是一個平均值接近零的稀疏計數，CV 的分母一小就會發散，所以那個 145.5% 至少混進了分母造成的假象，不能直接讀成「取樣完全沒開始」—— 診斷方式見指標可靠性。） 怎麼看出它壞了：把所有你要報的量的 CV 排成一張表 —— 若跨了兩個數量級，你就知道單一句「已收斂」是謊話。
-
不要在收斂還沒釐清之前，把高變異歸因給生物學。 為什麼會壞：噪音跟訊號長得一樣。你的 DNA RMSF CV 59.1% 看起來像「DNA 比蛋白質柔軟得多」這個很誘人的生物學結論，但拆開來看，主因是建模補出來的 5′ 單股末端在甩：res 1 到 11 的 RMSF 有 6 到 8 Å，晶體本體只有 3.9 Å。這不是 DNA 的物理，這是「這一段本來就沒有實驗座標」—— 而且這件事在 PDB 條目上就查得到：5F9R 的 non-target 股（chain D）只有 residues 12 到 30 有座標，1 到 11 完全無電子密度，是 PDBFixer 幫你補出來的。怎麼看出它壞了：把無實驗座標的區段標色排除，重算一次；數字掉下來就知道剛才差點寫錯結論。
-
不要指望純 MD 取樣到大構形轉變。 為什麼會壞：HNH 的構形檢查點翻轉在微秒以上，你的 4 × 87 ns 差了好幾個數量級。怎麼看出它壞了：數那個轉變在你的軌跡裡發生幾次。0 次不是證據，1 次也不是證據 —— 要從轉變估任何速率或自由能，你需要的是很多次，而那要靠 enhanced sampling（metadynamics、GaMD、umbrella sampling）或 Markov state model，不是靠再跑久一點。
-
不要把「收斂」誤當成「答案對」。 為什麼會壞：收斂只保證你把這個力場、這個起始結構、這個質子化狀態下的分布取樣清楚了，它對系統誤差完全沉默。你的專案有最直接的證據：Chai-1 預測結構起始給出逐殘基 0/8 共享、ρ≈−0.03（隨機）；5F9R 晶體起始給出 3/8 共享、Spearman +0.77。兩邊各自都可以宣稱「收斂」，答案卻不同 —— 元兇是起始結構，不是取樣。 收斂是誤差的下界，不是上界。
08Project Lens跟我的 project 多相關
為什麼是五星，不是四星： 因為這一頁就是這個專案的產品。
目前多數變體篩選流程都由數道 gate 與一連串預測／打分工具串接而成 —— 從序列設計模型、蛋白–核酸親和力預測器，一路到 MM-PBSA 與 FEP —— 但很少有哪一層在問「這個分數換條軌跡還會不會是同一個」。可靠性分層表要補的正是這個缺口。整個 methodology paper 的骨架是：先證明哪些量在獨立 replica 之間可重現，再讓篩選流程只用那些量打分。 沒有這一頁，你手上就只是「又一組 Cas9 MD 軌跡」；有了這一頁，你手上是一個別人可以照著做的判準。
再加三個理由：
- 它已經救過這個專案一次。 6/8 → 3/8 那件事若沒被抓到，逐殘基接觸就會變成打分軸，下游每一個變體排序都會建立在隨機數上。這種「已經產生實際價值」的方法學，是可防守的五星。
- 成本為零。 在「便宜」是硬約束的 workflow 裡，這是唯一一個不用多花任何 GPU 時數、卻直接決定其他所有結果可信度的環節。CV 表是既有軌跡的後處理。
- 它是 reviewer 的第一個攻擊點，也是你唯一準備好的防區。 多數 Cas9 MD 論文答不出「你跑了幾條、散布多大」，你答得出，而且答案本身有結論。
為什麼還不能宣稱做完了： 目前的證據是一個系統、一個力場、一個起始結構、4 條 replica。要把「聚合量可重現、逐位點量不可重現」升格成一個通用的方法學主張，至少還需要在第二個系統上重現同一組 tier —— 最自然的候選是換一組帶 mismatch 的底物條件，因為完全配對的底物在結構上就給不出高保真變體的專一性訊號（on-target 上的 eSpCas9(1.1) 在設計上本來就該是 null）。另外，n=4 算出來的 CV 本身也很粗，tier 邊界上的量（protein RMSD 12.2%、non-target 接觸 17.7%）不該被當成確定的分類。
具體怎麼用： 把 CV 表做成 workflow 的准入閘門：任何要進打分軸的量，必須先通過跨 replica CV < 10%。這條規則寫進工作流驗證，並在每次換系統時重跑一次 —— 因為 tier 是系統相依的，不是普世的。
09Reviewer Thinkingreviewer 會問什麼
4 條 87 ns 你憑什麼說收斂？點開看參考答案
分三層答，而且要主動講出不漂亮的部分。第一，前後半漂移：聚合量落在 2 到 17%，rep1、rep3、rep4 走平，rep2 尾段仍在微升 —— 這句要自己說。第二，跨 replica CV：綠燈的四個量在 2.1 到 9.5%。第三，也是最重要的，收斂宣稱只涵蓋聚合量；逐位點與 DNA RMSF 明確標為未收斂，並且沒有任何結論建立在它們上面。把收斂宣稱的邊界畫清楚，比宣稱得更大聲有用得多。
你的 10% 和 20% 門檻哪裡來的？點開看參考答案
誠實回答：這是本研究自訂的操作性分層，不是文獻標準。可防守的理由是它跟下游決策直接掛鉤 —— 只有綠燈的量會被用在打分軸上，所以門檻的意義是「這個量能不能承擔排序責任」。要更硬一點，就做敏感度測試：把門檻改成 15%，結論（聚合量可用、逐位點不可用）會不會翻？答案是不會，因為兩群之間有明顯的間隙 —— 綠燈最高 9.5%、紅燈最低 23.9%，中間只有兩個黃燈。那個間隙本身就是最好的辯護。
只有 4 個樣本算出來的 CV，它自己的誤差多大？點開看參考答案
這是最尖銳的一問，而且沒有辦法四兩撥千斤。n=4 的標準差估計本來就很不精確。誠實做法有兩個：一是報 CV 的同時給 bootstrap 信賴區間；二是承認 tier 邊界上的量不該被當成確定分類，明說 protein RMSD 12.2% 與 non-target 接觸 17.7% 是「未定」而非「marginal 這一類」。注意這一問不會摧毀你的主結論，因為主結論靠的是綠燈與紅燈之間一個數量級的差距，而不是靠某個量落在 9% 還是 11%。
你怎麼排除「四條 replica 只是同一個錯誤重複了四次」？點開看參考答案
排除不了，而且要直接承認。四條共享起始座標、力場、質子化狀態，所以跨 replica CV 量的是「同一個起始盆地內的取樣散布」，不是「起始結構不確定性」。後者需要換起始結構才能測 —— 而你剛好做過：Chai-1 起始給 0/8、ρ≈−0.03，5F9R 晶體起始給 3/8、Spearman +0.77。這組對照本身就是答案：最大的系統誤差源是起始結構，而不是取樣。 把這兩件事分開報，是這個專案最強的方法學論點之一。
為什麼不直接用 enhanced sampling 把構形空間掃乾淨？點開看參考答案
因為目標和預算都不允許。這個 workflow 的設計約束是「便宜、可大量篩選」，而 metadynamics 或 GaMD 每個變體的成本會跳一個數量級，還要選 collective variable —— 選錯 CV 的失敗模式比取樣不足更難察覺。目前的定位很清楚：平衡態的介面參與量用便宜的 unbiased MD 加多 replica 就夠，需要跨越 HNH 檢查點的問題留給後續的 enhanced sampling 或 Markov state model。 把這條界線寫進 limitations，比假裝已經涵蓋要安全得多。
10Common Mistakes最常犯的錯
- 只跑一條軌跡，用它自己的時間序列標準差當誤差棒。 後果：誤差被系統性低估，因為那條軌跡看不到自己沒去過的地方。正解：至少 3 到 4 條獨立 replica，報跨 replica 散布；單條的內部誤差棒只能當輔助。
- 在 n=2 的時候下逐殘基結論。 後果：6/8 這種漂亮數字會直接進簡報和下游的打分流程，然後在 n=4 崩掉。正解：n=2 只當 pilot；任何逐位點宣稱都等到 n≥4。
- 拿 RMSD plateau 當唯一收斂證據。 後果：對正交自由度完全盲目，收斂宣稱沒有保障。正解：對每一個要報的量各做一次前後半漂移，不要靠 RMSD 代言。
- 把一條長軌跡切段當獨立 replica。 後果：散布被低估，誤差棒假性收窄。正解：獨立 replica 必須換初始速度種子並各自重新平衡，方法段要寫明。
- 報「模擬已收斂」而不指定觀測量。 後果：同一條軌跡對 BSA 收斂、對 groove 接觸完全沒收斂，一句話同時是真也是假。正解：收斂宣稱一律綁定到具體的量，並附 CV 表。
- 把建模補出來的區段留在統計裡，然後解讀成生物學。 後果：5′ 端在甩會被讀成「DNA 很軟」。正解：標色排除無實驗座標的區段，重算後再解讀。
- 在軌跡裡看到一次構形轉變就開始講機制。 後果：n=1 的事件沒有任何統計意義，reviewer 一眼看穿。正解：報發生次數；次數個位數就只能寫成觀察，不能寫成結論。
- 存太密的幀，以為幀數就是樣本數。 後果：硬碟爆掉但有效樣本數沒增加。正解：先估自相關時間再決定輸出節奏；Neff = N/g 才是你的 n。
11Further Reading讀哪幾篇
12Summary三句話
收斂在 MD 裡是三件不同的事：系統還在不在鬆弛、平均值的誤差多大、以及有沒有整片構形空間沒去過，而最後一件原則上只能反證不能證明。你的 4 條 replica × 約 87 ns 給出的答案很明確 —— 聚合量前後半漂移 2 到 17%、跨 replica CV 在 10% 以下，可以承擔排序責任；逐位點量與 DNA RMSF 的 CV 在 24 到 146%，只能當定性敘述。這一頁最貴的一課是逐殘基 top-8 共享從 n=2 的 6/8 掉到 n=4 的 3/8，小 n 的偏誤是單向的，它只會讓你太樂觀。