從 PDB 到軌跡From PDB to trajectory
建系統與跑模擬這一段的每個決定與它的代價。
01一句話只有一句
建系統的每個決定都很便宜，但幾乎每一個都不可逆。
02Why should I care?我為什麼要讀這頁
因為這一段是唯一你還能用零成本改變結果的地方。軌跡一旦跑出來，起始結構、力場、盒型、底物設計全部凍結在裡面，之後所有分析都只是在讀這些決定的後果。
一個具體到殘忍的例子：同樣的分析腳本、同樣的建系統流程，只換起始結構。Chai-1 預測結構起始 → 逐殘基 top-8 接觸殘基跨 replica 0/8 共享、ρ ≈ −0.03；5F9R 晶體起始 → 3/8 共享（450、695、765）、pairwise Jaccard 0.48、Spearman +0.77。從隨機變成有訊號，代價是零 GPU-hour，只是在 wget 的時候打了不同的四個字元。
反過來的例子一樣殘忍。菱形十二面體盒相對於同樣 image distance 的立方盒，體積只有 71%，GROMACS 手冊逐字寫著這會省下約 29% 的 CPU 時間。注意這個百分比的方向：29% 是「用菱形十二面體比用立方盒省下的」；反過來從菱形十二面體換成立方盒，溶劑量是 1/0.71 ≈ 1.4 倍。你的系統是 535,265 atoms，其中絕大多數是水；換成立方盒會落到七十幾萬原子的量級（這是幾何換算，不是實測）。整條流程最貴的一步，被一個下拉選單省掉了將近三成。
這件事對你不陌生。菱形十二面體就是 FCC 晶格的 Wigner–Seitz cell——你在 VASP 裡選 primitive cell 而不是 conventional cell，做的是完全一樣的事：同樣的物理、同樣的 minimum image 距離、少算一堆重複的體積。差別只在於材料人天生會這樣做，而生物 MD 的教學經常從立方盒開始講，然後就沒有人回頭改了。
03What這是什麼
從 PDB 到能跑的軌跡，中間是一串有順序的決定。順序重要，因為後面的決定會被前面的鎖死。
D1 — 起始座標。 這是全流程槓桿最大的一個決定。選項：實驗晶體／cryo-EM 結構、預測結構（AlphaFold、Chai-1）、或同源建模。代價：零算力，但決定了逐殘基層級的結論能不能存在。
D2 — 系統邊界。 apo 蛋白？蛋白 + sgRNA？完整三元（蛋白 + sgRNA + target/non-target DNA）？你手上的參考點是完整 R-loop 三元 535,265 atoms；apo 蛋白少了 sgRNA、雙股 DNA 以及它們帶進來的水與抗衡離子，會小上一截（確切數字去查你自己那個 run 的 log，不要憑印象引用）。代價：原子數大致線性反映在每步的計算量上。
D3 — 缺失片段與端點。 晶體結構一定有缺的 loop，核酸一定有截斷的端點。你補或不補、怎麼補，會直接製造出高柔軟度的假訊號。
D4 — 質子化與金屬。 His 的三種質子化態、活性中心的 Mg²⁺ 數量與配位。古典力場不會幫你決定，也不會在錯的時候報錯。
D5 — 力場組合。 蛋白一套（ff14SB 或更新版）、DNA 一套（OL15）、RNA 一套（OL3／χOL3）、水一套（TIP3P 或 OPC）、離子參數一套。這是這條流程的 exchange-correlation functional：五套獨立發展的參數要湊在一起用，而它們的相容性沒有變分原理保證。
D6 — 盒型與 padding。 菱形十二面體對近球形溶質是預設答案。padding 的慣例下限大約在 1.0 到 1.2 nm 之間，關鍵判準是 minimum image 距離要大於兩倍 cutoff，而且要留給溶質在模擬中可能的伸展。
D7 — 中和與鹽濃度。 Cas9 三元複合體帶大量核酸負電荷，中和需要的離子數不少。慣例是中和之後再加到約 0.15 M 的生理鹽濃度。
D8 — Minimize → 加熱 → NPT → 逐步鬆綁約束。 這一串是有階梯的：先固定重原子只讓水鬆弛，再逐步放鬆蛋白側鏈、主鏈、核酸。跳階梯會在 production 的前幾 ns 看到不該有的大幅重排。
D9 — Integrator、timestep、約束。 有 H 原子鍵長約束時 2 fs 是安全值；配合 hydrogen mass repartitioning（HMR）可以推到 4 fs，等於同樣的 wall time 多拿一倍取樣。這是整條流程裡投報率最高的單一設定，但它不是免費的預設值——原始文獻只在胜肽與蛋白上驗過，沒有核酸，所以要自己做一次對照（見 When NOT 第 5 條）。
D10 — 恆溫恆壓與長程靜電。 Langevin 或 Nosé-Hoover、Monte Carlo barostat、PME。PME 對你來說就是 Ewald，跟 VASP 處理帶電系統的邏輯同源。
D11 — Replica 策略。 幾條、怎麼產生獨立性、seed 怎麼記錄。這一項決定了你最後能不能算誤差棒，也就是決定了整篇 paper 能不能防守。
D12 — 輸出頻率與儲存。 家目錄配額通常只有百 GB 量級，這不是細節，是會讓 job 在第 60 ns 死掉的硬約束。
04Why為什麼重要
因為這一段決定的不是「精度」，是你的結論存在於哪個尺度。
D1 的例子已經給了：起始結構把逐殘基尺度的訊號從有變成無。原因不難懂——Chai-1 對 HNH 的 pLDDT 只有 62.9，target DNA 64 到 65，最低的一段是 L1 linker 接進 HNH 的區間（res 765 到 809，pLDDT 46；UniProt 標 HNH 從 770 起、RuvC-II 到 765 止，所以這一段橫跨 linker 與 HNH 前緣）；而 REC lobe 有 78 到 85、RuvC-I 有 80。你從一個在檢查點區域信心最低的模型出發，那個區域的逐殘基接觸當然在每條 replica 裡跑去不同的地方。低信心不代表那裡是隨機噪音，而是代表模型對那裡的擺法沒把握，你很可能被放進一個結構上自洽、動力學上錯誤的盆地。
D3 的例子更直接。你補建的 5′ 單股末端（res 1 到 11）RMSF 有 6 到 8 Å，而晶體本體只有 3.9 Å。這些原子沒有任何實驗座標約束，它們的高柔軟度反映的是「這裡本來就沒資料」。如果它們留在 domain 平均值裡，就是用建模的自由度污染物理結論。
D11 最不直觀但最貴。你的 CV 分層結果——BSA 2.1%、總 protein–DNA 接觸 4.1%、target-DNA 接觸 9.4%、鹽橋 9.5% 是綠燈；protein RMSD 12.2%、non-target 接觸 17.7% 是黃燈；interface iRMSD 23.9%、DNA RMSF 59.1%、groove 距離 48%、逐位點距離 38 到 46%、groove 接觸計數 145.5% 是紅燈——全部是 replica 之間的變異，不是時間上的變異。這代表你買 replica 買到的是誤差棒，買長度買到的主要是同一個盆地裡更細的分布。兩者不能互換。
建系統階段的錯誤不會在建系統階段被發現。它們會偽裝成分析階段的「這個 metric 不可靠」。所以判斷一個 metric 好不好用之前，先確定你不是在評價自己的建模決定。
05How怎麼做
最小可行配置（MVP）
這是「我今天要開始，最少要決定什麼」的版本。目標是拿到一組能算出誤差棒的 WT 軌跡。
- 起點：5F9R 完整 R-loop 三元晶體結構（Jiang 等人 2016）。這個結構捕捉的是 R-loop 已形成、非 target 股被推到 RuvC 附近、準備切割的狀態，正好是你想談檢查點的那個階段，而且 deposition 沒有標註去活突變（是催化活性型序列）。同時先知道它缺什麼：解析度 3.40 Å；non-target 股 30 個核苷酸只有 res 12 到 30 有電子密度（19/30），res 1 到 11 完全沒有——而那正是要走進 RuvC 的那一段。你補出來的那 11 個殘基會主導所有 RuvC 相關的結論，這是這個起點最大的建模風險，必須寫在方法段裡。
- 系統邊界：完整三元。不要為了省錢先跑 apo——apo 的 HNH RMSF（1.86 Å）跟三元的（1.61 Å）雖然都指向同一個結論，但 apo 沒有介面，你的綠燈 metric（BSA、protein–DNA 接觸、鹽橋）一個都算不出來。
- 缺失片段：補，但全部記錄下來，並在分析腳本裡預設排除。做法是產生一個 residue mask 檔案，跟軌跡放在一起。
- 力場：ff14SB（蛋白）＋ OL15（DNA）＋ OL3／χOL3（RNA）＋ TIP3P（水）＋ 對應的離子參數。這是目前生物分子 MD 對「蛋白 + 核酸」系統最常見的組合。若你之後要談水合結構的細節，OPC 是更好的水模型，但換水模型等於整批重跑，要在這一步決定。
- 盒型：菱形十二面體，padding 取 1.0 到 1.2 nm 這個量級。
- 離子：中和後加到約 0.15 M。
- 平衡階梯：能量最小化 → 位置約束下加熱到 300 K → NPT 下逐步鬆綁約束 → 無約束 NPT 平衡。每一階段都看能量、體積、密度是否平了再往下走。
- Production：2 fs 加 H 鍵長約束，或 4 fs 加 HMR。300 K、1 bar、PME。
- Replica：4 條。這是能報 CV 的實務下限，而你自己的 n = 2 → n = 4 經驗（逐殘基一致性 6/8 掉到 3/8）就是不能只跑 2 條的理由。但要知道下限的代價：從 4 個點估出來的標準差，本身的相對誤差就有 1/√(2(n−1)) ≈ 40%。所以你的 CV 表要當成粗略的三檔紅綠燈，不要當成小數點後一位的排序——12.2% 與 17.7% 之間在 n = 4 下其實分不太開。
- 長度：以 ~87 ns（你已有的參考點）為第一個里程碑，然後看目標量的後半段有沒有單調漂移再決定要不要延長。
- 輸出：只存 solute 的高頻軌跡 + 低頻的全系統快照。先算好 GB 數再送 job。
成本量級：五十萬原子等級的 explicit-solvent 系統在單張現代資料中心 GPU 上，公開 benchmark 大致落在每天數十 ns，所以一條 ~87 ns 是單卡數天。這只是量級——ns/day 對 timestep、precision 模式與輸出頻率極度敏感，光是 2 fs 換 4 fs 就差兩倍，所以看到別人的 benchmark 沒附這些設定就當它不存在。請務必自己跑一個 1 ns 的計時 job 校準，之後所有排程都用實測值推；引用別人的數字排自己的程，是被佇列卡死的標準原因。
若算力加倍，先加在哪
按投報率排序，不要跳過前面直接做後面：
- 加 replica，4 → 8。 理由不是偏好，是你的 CV 資料：主要誤差來源在 replica 之間。但要講清楚它買到什麼、買不到什麼——CV 不會因為多跑 replica 而變小，它是母體的散布，不是平均值的誤差；縮小的是平均值的標準誤，而且只按 1/√n 縮。以 non-target 接觸（CV 17.7%）為例：n = 4 時平均值的相對誤差約 8.9%，n = 8 降到約 6.3%。所以黃燈的量不會「變綠」，但你能可靠偵測的最小效應量會縮小約 √2 倍——這才是它有機會被拉進打分層的實際理由，也是唯一能擴大打分公式可用範圍的投資。
- 加條件，特別是 mismatched 底物。 eSpCas9(1.1) 在 on-target 底物上是 null，而且它本來就該 null。你要看到專一性效應，必須有錯配底物這個維度。加一個底物條件比把單一條件跑到兩倍長有用得多。
- 才是加長度。 只有當你的目標量在後半段還在單調漂移時，加長度才是對的答案。否則你買到的是同一個盆地裡更平滑的分布。
- 最後才是加系統大小或加解析度（極化力場之類）。 這會讓你同時失去統計力與可比性，是投報率最差的花法。
砍長度，不要砍 replica。一組 4 × 40 ns 能算誤差棒，一組 2 × 80 ns 不能。前者能寫成 methodology paper，後者只能寫成觀察報告。
怎麼產生「獨立」的 replica
從同一個平衡構形分支、只換初速度亂數種子，是可接受的下限，但你要知道它們共享同一個溶劑構形與同一段 equilibration 歷史，軌跡早期是相關的。更穩健的做法是從溶劑化那一步就分開，四個盒子各自獨立加水加離子、各自平衡。成本只多了幾小時的建系統時間，換來的是 reviewer 問「你的 replica 真的獨立嗎」時有答案。
絕對不要做的事：把一條長軌跡切成四段當成四條 replica。段與段之間高度相關，CV 會被系統性低估，你會以為自己的量比實際上可靠。
06When什麼時候用
- 每一次要產生新的 production 軌跡時。 這一頁是 checklist，不是讀完就算。
- 當你要換任何一個上游設定時（力場、水模型、盒型、起始結構）——換了就等於開一個新的比較群組，舊軌跡不能混進來比。
- 當你要把 WT 的設定複製到變體上時。 逐項核對，任何一項不同，Δ(變體 − WT) 就同時包含了設定差異。
- 當算力預算改變時。 上面那個排序就是決策表。
- 當 reviewer 要求補對照組時。 你會需要知道哪些設定可以只重跑一部分、哪些必須整批重來。
07When NOT什麼時候別用
- 不要用低信心的預測結構當 production 起點。 為什麼會壞： 預測結構在動態區域的信心最低（你的 Chai-1 模型：HNH 62.9、target DNA 64 到 65、L1 linker 進 HNH 那一段只有 46），而那恰好是你想量的區域。最可能的解釋是模型把起點放進一個結構上自洽、但動力學上不對的盆地，每條 replica 從那裡各自漂走；還有一個競爭解釋是那個區域本來就柔軟，任何起點都會散。這兩者你目前的資料分不開——但不必分開也能做決定，因為換成晶體起點就把訊號救回來了。怎麼看出它壞了： 逐殘基層級的 top-N 在 replica 之間幾乎沒有交集（0/8、ρ ≈ −0.03），而且逐殘基 RMSF 與 pLDDT 呈反相關（ρ = −0.36）——後者告訴你模型自己已經在示警。替代做法： 預測結構拿來做 QC 與挑目標區域，production 用晶體座標；真的沒有晶體時，把所有結論限制在聚合量層級。
- 不要把建模補出來的片段留在分析裡。 為什麼會壞： 這些原子沒有實驗約束，它們的高柔軟度是建模自由度不是生物學。你的 5′ 端 res 1 到 11 RMSF 有 6 到 8 Å，本體只有 3.9 Å——少數殘基就能把整個 domain 的平均值拉走。怎麼看出它壞了： RMSF 曲線的高峰全部落在鏈端與你補的 loop 上；把它們排除後結論就翻掉。
- 不要為了省錢把高電荷的核酸複合體改用 implicit solvent。 為什麼會壞： implicit solvent 沒有顯式離子，對 Cas9 三元這種高電荷密度系統的靜電屏蔽處理不足，而 R-loop 的穩定性本質上就是屏蔽問題。怎麼看出它壞了： 核酸骨架彼此靠得比晶體近、R-loop 有塌陷傾向、鹽橋數量異常偏高。這幾個徵狀會同時出現。注意： implicit solvent 對小的、電中性的體系（例如做教學用的 Trp-cage）仍然完全合理，這條限制是針對帶大量電荷的核酸複合體。
- 不要在同一批比較裡混用不同的盒型、水模型或力場版本。 為什麼會壞： 這些設定會系統性地平移你的 metric，差異會被誤讀成突變效應。怎麼看出它壞了： 你的 Δ(變體 − WT) 與「這個變體是哪一批跑的」相關性，比它與突變位置的相關性還高。預防方式： 每條軌跡的 metadata 存成一個 JSON 跟軌跡放一起，分析時強制比對。
- 不要把 timestep 推到 4 fs 卻沒做 HMR，也不要推上去之後不檢查。 為什麼會壞： 含 H 的高頻振動週期大約在 10 fs 量級，timestep 太大時積分器會系統性地把能量灌進這些自由度。怎麼看出它壞了： 在 NVE 測試裡總能量單調漂移；或在 NVT 裡溫度需要 thermostat 持續移除能量。做法： 換 timestep 之後跑一段短 NVE 檢查能量守恆，這是幾分鐘的事，但幾乎沒人做。
- 不要在「同時只允許少數並行 job、家目錄配額只有百 GB 量級」的環境下，一次送出 8 條 replica 的長軌跡。 為什麼會壞： 不是物理錯誤，是工程錯誤——job 會排隊、磁碟會滿、checkpoint 會寫不進去，然後你在第 60 ns 損失整批。怎麼看出來： 這個要在送出之前算，不是看出來的。先算
輸出頻率 × 原子數 × 軌跡長度 × replica 數，再看你的配額夠不夠。
08Project Lens跟我的 project 多相關
為什麼不是更低： 因為這一段是你在整條流程裡唯一能用零成本改變結論的地方，而你的專案剛好有一個乾淨的因果證據可以證明這件事。
- 一個決定，把訊號從無變有。 只換起始結構（Chai-1 → 5F9R），逐殘基接觸的跨 replica 一致性從 0/8、ρ ≈ −0.03 變成 3/8、Jaccard 0.48、Spearman +0.77。零 GPU-hour 的改動，救回了整個逐殘基尺度。這是全站最強的「上游決定支配下游結論」實例。
- 一個下拉選單，省掉整條流程最大開銷的三成。 菱形十二面體盒相對同 image distance 的立方盒體積只有 71%（GROMACS 文件：約省 29% CPU 時間）。在「便宜」是硬約束的 project 裡，這種等級的免費午餐必須吃。
- 這一段決定了誤差棒存不存在。 你所有可用的打分量（BSA 2.1%、總接觸 4.1%、target-DNA 接觸 9.4%、鹽橋 9.5%）之所以能被判定為可用，前提是有 4 條獨立 replica。這個前提是在 D11 決定的，不是在分析時決定的。
風險要誠實講： 這一段目前有一個沒解掉的方法學缺口——力場組合的相容性沒有針對 Cas9 這類蛋白–RNA–DNA 三元系統做過系統性驗證。ff14SB、OL15、OL3 各自都有紮實的基準，但它們湊在一個數十萬原子的三股介面上表現如何，文獻上沒有給你一個乾淨的答案。這必須寫進 limitations。
09Reviewer Thinkingreviewer 會問什麼
你為什麼選這個起始構形？換一個構形會不會給不同答案？點開看參考答案
會，而且你有數據可以證明會。誠實的答法是兩層：第一，說明 5F9R 捕捉的是 R-loop 已形成、準備切割的狀態，這正是你想討論檢查點的階段；第二，主動報告你用預測結構跑過同一套分析，逐殘基層級變成隨機（0/8、ρ ≈ −0.03）。把這個對照當成結果而不是失敗，它直接支撐了「起始結構是逐殘基指標的支配因素」這個 methodology 論點。
你的蛋白、DNA、RNA 三套力場相容嗎？有沒有交叉驗證？點開看參考答案
這一問沒有完美答案，不要硬掰。可以講的是：三套都是同一個 AMBER 家族、共用同一套 Lennard-Jones 混合規則與同一個水模型，這是目前的社群慣例而不是我的創新。可以提供的證據是內部一致性檢查——晶體可解析區的 RMSD 是否維持在合理範圍、鹼基配對是否維持、螺旋參數有沒有系統性漂移。明說「這是慣例，不是驗證過的最佳組合」，比宣稱它經過驗證安全得多。
你的 replica 是怎麼產生的？它們真的獨立嗎？點開看參考答案
要準備確切的答案：是從同一個平衡構形換初速，還是從溶劑化就分開。如果是前者，主動說明軌跡早期是相關的，並展示你丟掉了多少 equilibration。可以補的定量證據是：把每條 replica 的前半與後半分開算目標量，看 replica 之間的差異是否大於 replica 內部前後半的差異——如果是，那 replica 之間的變異是真的。
補建的區段占多少？你怎麼確定結論不是來自它們？點開看參考答案
給出確切的殘基清單與比例，並展示排除前後的結果。你的 5′ 端 res 1 到 11（RMSF 6 到 8 Å，本體 3.9 Å）是最好的示範案例：把它們留著會怎樣、排除後結論是否不變。這一問答得好，圖的可信度會整個提高一級。
你的 timestep 和約束設定有沒有影響到你關心的自由度？點開看參考答案
如果用了 HMR 到 4 fs，reviewer 有理由問這會不會影響含氫的動力學。可以給的證據是一段短 NVE 的能量守恆檢查，以及同一段用 2 fs 重跑的對照——後者只要幾 ns 就足以說明目標量沒有系統性偏移。這個對照很便宜，但它把一個潛在的攻擊點變成一張補充圖。
10Common Mistakes最常犯的錯
- 用預測結構跑 production 然後報逐殘基結論。 後果：結論在統計上是隨機的（0/8、ρ ≈ −0.03）。正解：晶體座標當起點，預測結構只做 QC 與區域挑選。
- 忘了記錄哪些殘基是補建的。 後果：半年後你自己也分不出哪些高 RMSF 是真的。正解：建系統時就輸出 residue mask，跟軌跡一起版本控制。
- 用立方盒。 後果：多算三成的水，整條流程最大的開銷憑空變貴。正解：近球形溶質用菱形十二面體，這是 FCC 的 Wigner–Seitz cell，跟你在 VASP 選 primitive cell 是同一個道理。
- 平衡階梯跳級。 後果：production 前幾 ns 出現大幅重排，你以為那是生物學。正解：約束逐步鬆綁，每一階段確認能量與密度平了再往下。
- WT 與變體用了不同批次的設定。 後果：Δ 值裡混著設定差異，而且無法事後拆開。正解：把設定寫成一個 config 檔，變體只改 mutation 欄位，其餘由腳本繼承。
- 把一條長軌跡切段當 replica。 後果：CV 被低估，你會高估自己的可靠性。正解：真正獨立的初速與（最好是）獨立的溶劑化。
- 輸出頻率沒算就送 job。 後果：job 在 60 ns 因為磁碟滿而死，checkpoint 也寫不進去。正解：送出前算一次
輸出頻率 × 原子數 × 長度 × replica 數，對照你的家目錄配額。 - 先跑 apo 省錢。 後果：拿到一組漂亮的 RMSF，但你的綠燈打分量（BSA、protein–DNA 接觸、鹽橋）一個都算不出來。正解：如果目標是介面打分，第一天就跑完整三元。
11Further Reading讀哪幾篇
12Summary三句話
從 PDB 到軌跡這一段的所有決定都幾乎不花算力，卻把後面每一個數字的尺度與誤差棒都鎖死了。最小可行配置是 5F9R 完整三元、AMBER 家族力場、菱形十二面體盒、HMR 4 fs、4 條真正獨立的 replica，而算力加倍時的正確花法依序是加 replica、加錯配底物條件、最後才加長度。判斷一個 metric 好不好用之前，先確定你不是在評價自己的建模決定。