系統建立System preparation
從一個 PDB 到能跑的系統，中間十幾個決定每一個都能毀掉結果。
01一句話只有一句
PDB 不等於系統：中間十幾個決定，每一個都能安靜地毀掉結論。
02Why should I care?我為什麼要讀這頁
因為建系統的錯誤幾乎都不會炸給你看。
這一步相當於你在 VASP 裡寫 POSCAR、選 POTCAR、設 INCAR。差別在於：VASP 的錯誤通常會用力砸在你臉上 —— 收斂不了、能量荒謬、離子步發散。而 MD 的建系統錯誤會安靜地跑完 87 ns，吐給你一條漂亮的軌跡、一張好看的圖、和一個錯的結論。這是這一頁存在的全部理由。
這一步有兩種典型的失敗，這個專案兩種都遇過：一種真的付了 GPU 時數，另一種在 build 階段就攔下來了。
真的付了錢的那種是起始結構。 用 Chai-1 預測結構起始時，逐殘基 top-8 接觸殘基跨 replica 是 0/8 共享、ρ ≈ −0.03（等同隨機）。換成 5F9R 晶體起始後，變成 3/8 共享 [450, 695, 765]、pairwise Jaccard 0.48、Spearman +0.77。同一套 metric、同一套分析腳本、同樣的 replica 數。元兇是起始結構，不是 metric。 那幾條 Chai 起始的 replica 是真的跑完才知道的。
攔下來的那種是盒子。 5F9R 完整 R-loop 三元複合體第一次 build 出來是 762,027 個原子，改了盒型之後降到 535,265。這一次因為在 build report 就看到原子數不對勁，還沒投長跑就改掉了 —— 省下的 30% 不是靠犧牲任何科學嚴謹性換來的，是靠改一行程式。在算力配額與「便宜」都是硬約束的現實下，這 30% 直接換成你能跑幾條 replica，而 replica 數又直接決定你的 CV 表可不可信。下面 ## How 第 9 節會完整拆這個決定。
兩者的差別只有一件事：有沒有人在投長跑之前看 build 的輸出數字。
03What這是什麼
建系統就是把 U(R) 的全部內容與邊界條件一次決定完。在 DFT 裡這件事分散在 POSCAR（原子與晶格）、POTCAR（贗勢）、INCAR（數值設定）三個檔案；在 MD 裡它們全部擠在一個 build script 裡，而且大部分決定之後不能改。
這是你必須做的決定清單。它同時就是你之後方法段要逐條交代的東西：
- 起始座標 —— 哪個 PDB？晶體、cryo-EM，還是預測結構？
- 保留哪些 heterogen —— 結晶助劑、金屬離子、結晶水，各留不留；還要分辨出哪些被登記成 heterogen 但其實是共價接在鏈上的非標準殘基（刪了等於改序列）。
- 缺失殘基與 loop —— 補哪些、用什麼補、怎麼標記。
- protonation state —— pH、每個 His 的三態、任何位在活性中心或界面的可解離基團。
- 力場 —— 蛋白、DNA、RNA、離子，四套要互相相容。
- 水模型 —— TIP3P、OPC、TIP4P-Ew…，而且它跟第 5 項是綁在一起的。
- 盒型、padding、離子強度 —— 這一項直接決定成本。
- 約束、hydrogen mass、timestep —— 決定你一天能跑幾 ns。
- 最小化與分階段平衡 —— 決定前幾百 ps 的「構形變化」是真的還是人造的。
其中第 1、3、4 項是單向門：跑下去之後沒有辦法用更多模擬把它們救回來。第 7、8 項是成本旋鈕。第 5、6 項介於兩者之間 —— 可以換，但換了就得整套重跑。
04Why為什麼重要
因為這一步的錯誤有一個共同特徵：它們製造的假訊號，長得跟真訊號一模一樣。
補出來的 loop 會甩，甩出來的高 RMSF 跟真正的柔軟區在同一張圖上、同一個顏色。protonation 選錯，氫鍵網路就整個換掉，但軌跡照跑、能量照收斂。padding 開太小，蛋白透過 PME 感受到自己的鏡像，靜電被人工週期性汙染 —— 而你在任何常規輸出裡都看不到這件事。
這跟你熟悉的 DFT 失敗模式很不一樣。DFT 裡選錯 POTCAR、k-mesh 太疏、真空層太薄，通常會在能量或收斂行為上留下明顯痕跡。MD 的建系統錯誤更像是選錯了 slab termination：計算本身完全正常，只是你從頭到尾在算另一個系統。
所以這一步的正確心態是：把 build 當成一個要被 review 的產出，不是一個要趕快跑過去的步驟。
05How怎麼做
以下每一節對應上面的決定清單。用你的專案的實際做法當骨架。
1. 起始座標：晶體優先，預測結構要當作假設。
5F9R 是 catalytically-active WT SpCas9（沒有 H840A 這類失活突變）＋ sgRNA ＋ 完整 R-loop dsDNA，X-ray 解析度 3.40 Å，論文標題就叫「primed for DNA cleavage」—— 它捕捉的正是切之前那個構形，DNA 螺旋被折彎 30°。這是為什麼它比一般的三元結構更適合你的問題。「完整」是相對於 non-target 股被截斷的三元結構（例如 4UN3）而言：兩條股都在。但它自己的 non-target 股仍有 11 個 nt 沒有電子密度，見下面第 3 節 —— **「完整」不等於「全部解出來」，這兩件事在方法段要分開講。**預測結構不是不能用，但它必須被當成一個待驗證的假設，而不是起點（實測見上面的 0/8 vs 3/8）。詳見 Cas9 結構 與 結構預測。
2. 清掉什麼、留下什麼。
fixer.removeHeterogens(keepWater=False)。這一行會把 PDB 裡所有非標準殘基清掉，所以你必須先知道自己清掉了什麼。5F9R 的非標準項目只有兩種，而它們的性質完全不同：
- SO₄²⁻（三個：A201、B1401、C101） —— 結晶條件帶進來的硫酸根，不是生物學配體，拿掉沒有爭議。
- GTP（chain A residue 1） —— 這個不是結晶助劑。 打開 PDB 檔會看到
SEQRES 1 A 118 GTP G C G C ...，還有一行LINK O3' GTP A 1 → P G A 2：它是 sgRNA 的 5′ 端第一個核苷酸，因為 sgRNA 是 T7 in vitro transcription 做出來的，5′ 端保留了三磷酸，PDB 因此把它登記成非標準殘基。所以removeHeterogens在這裡做的事是：把 sgRNA 從 118 nt 截成 117 nt（順帶拿掉那個 5′ 三磷酸）。
removeHeterogens 是靠殘基名稱查標準表，它分不出「泡在旁邊的結晶助劑」和「共價接在鏈上的非標準殘基」。5F9R 的 GTP 就是後者。這件事對你的結論影響很小（少一個 5′ 核苷酸、遠離 protein–DNA 界面），但方法段必須寫成「sgRNA 5′ 端的 GTP 一併移除，模擬的是 117 nt 的 sgRNA」，而不是「移除結晶助劑 SO₄ 與 GTP」 —— 後者的說法會讓任何看過 5F9R 的 reviewer 知道你沒有逐項檢查過自己刪了什麼。通用做法：build 前把 HET 與 LINK 記錄各印一遍，有 LINK 接到聚合物鏈上的，就不是助劑。
Mg²⁺ 是最需要想清楚的一個：如果你要看催化，非加不可，而且要接受 fixed-charge 力場對二價離子的已知缺陷；如果你在看 engagement 與界面，可以不加。你的專案就是後者 —— 這一組 run 問的是 engagement 與界面接觸，不是催化，而且 5F9R 結構裡本來就沒有任何 Mg（非標準項目只有上面那三個 SO₄ 與 sgRNA 的 5′ GTP）。這個決定要寫進方法段，不要讓 reviewer 自己發現。
3. 缺失殘基：補，然後標記。
PDBFixer 的 findMissingResidues() → findMissingAtoms() → addMissingAtoms()。5F9R 要補 19 個 disordered residues，而這 19 個不是一團模糊的數字，PDB 的 REMARK 465 直接列給你：
- 蛋白 chain B：M1–D2（N 端 2 個）、L1365–D1368（C 端 4 個）
- sgRNA chain A：C117–G118（3′ 端 2 個）
- non-target DNA chain D：residues 1–11（11 個）
2 + 4 + 2 + 11 = 19。做這個加法本身就是一道品質關卡：如果 PDBFixer 報的數字跟 REMARK 465 對不起來，代表你的輸入檔被動過（截鏈、改編號、少了 header），先查清楚再往下。
5F9R 的 non-target 股是 30 nt，但只有 residues 12 到 30 有座標 —— 19/30，觀測率 63.3%。也就是說，R-loop（sgRNA 把雙股 DNA 頂開、自己去配對 target 股後留下的三股結構）裡那條被置換出來的單股，有三分之一在晶體裡根本沒有電子密度，是你自己補出來的。這不是「補個鬆散的末端」，這是你正在分析的那個東西本身有三分之一是模型。任何關於 non-target 股的結論都必須帶著這個前提講。
順帶一個看起來很像因果、但目前只是一致的觀察：你的 non-target 接觸 CV 是 17.7%（marginal），target-DNA 接觸只有 9.4%（reproducible）。補出來的 11-nt 是這個落差的合理嫌疑犯，但你沒有做過「把 res 1–11 排除後重算 CV」這個對照，所以現在只能寫成假說。要把它變成證據，就跑那個對照 —— 成本是零，因為軌跡已經在硬碟上。
真正關鍵的動作不是「補」，是標記。你的 build script 的做法值得抄成標準流程：先只輸出 solute（不含水），把每個原子的 B-factor 欄位寫成 1.0 = modelled / 0.0 = crystal，存成一個 prepared_solute PDB，先請熟悉這個結構的人（實驗端或共同作者）確認，確認過了才進入溶劑化。這樣做有三個回報：對方可以直接在 PyMOL 裡按 B-factor 上色看你補了什麼；後續每一張分析圖都能一鍵排除補出來的區段；reviewer 問「哪些是你自己畫的」時你有一個檔案可以直接給。
這件事的實際重量：補出來的 5′ 單股末端（res 1 到 11）RMSF 有 6 到 8 Å，晶體本體只有 3.9 Å。留在平均值裡就是純汙染。
4. protonation state：pH 只在加氫那一瞬間存在。
mod.addHydrogens(ff, pH=7.0)。這一行看起來人畜無害，但它做的事是：對每個可解離基團查一個標準 pKa 表，決定加不加氫，然後永遠凍結。His 有三態（HID / HIE / HIP），三者的重原子幾何幾乎一樣，力場不會抱怨，但氫鍵網路與金屬配位完全不同。
- 一般表面殘基：預設猜法通常沒問題。
- 埋在口袋裡、參與金屬配位、或位在界面上的 His：預設猜法經常錯。 Cas9 的 H840 是 HNH 活性中心的催化 His —— 這種位置不能交給預設值。
- 工具：PROPKA、H++、pdb2pqr 可以算結構相關的 pKa 位移。真的關鍵的位點，最好用兩個工具交叉檢查，不一致就在方法段誠實寫出來。
橋回你熟悉的世界：這相當於你在做吸附計算前選錯了吸附物的質子化態，或選錯 slab termination。之後所有的能量、所有的排序，全部建在那個沒被檢驗的選擇上。差別是 MD 更容易忘記自己做過這個選擇，因為它藏在一行 addHydrogens 裡面。
5. dangling 5′-phosphate：一個會直接擋住 build 的技術細節。
Amber 的核酸 residue template 對懸空的 5′-phosphate 會 template match 失敗。做法是把每條核酸鏈第一個殘基的 P、OP1、OP2 刪掉，變成 5′-OH。這是純技術修補，但要寫進方法段，因為它改變了那個末端的電荷。
注意這一步跟上面第 2 節是連動的：sgRNA 的原始 5′ 端是 GTP，被 removeHeterogens 拿掉之後，接棒的 G2 才是這裡被剝掉磷酸的那個殘基。也就是說 sgRNA 的 5′ 端經歷了兩次改動（GTP 移除、再變成 5′-OH），方法段要一次講完，不要拆成兩處。
6. 力場：整套用，不要混搭零件。
你寫的是 ForceField("amber14-all.xml", "amber14/tip3p.xml")，而 amber14-all.xml 其實是四個檔案的捷徑：protein.ff14SB + DNA.OL15 + RNA.OL3 + lipid17。也就是說那一行同時選了三套力場。知道自己選了什麼，是能討論它的前提：
- 蛋白 ff14SB（Maier 2015）：配 TIP3P 是被 validate 過的組合。Amber 現在的建議組合是 ff19SB 配 OPC，其中 ff19SB 的 backbone 參數是對溶液中的 QM 能量面逐胺基酸訓練的。
- DNA OL15（Zgarbová 2015）：OL15 這個名字指的是 χOL4 + ε/ζOL1 + βOL1 三項扭轉角修正疊在 ff99bsc0 之上的組合，2015 那篇貢獻的是 βOL1。它被驗證的對象是 B-DNA 與 Z-DNA 這類相對規則的雙股結構，沒有涵蓋被酵素撐開、解旋、彎折的非典型骨架 —— 而 Cas9 的 R-loop 正好整段都是那種骨架。這是你的系統一個誠實的、必須寫進 limitation 的空缺。
- RNA χOL3（Zgarbová 2011）：修 glycosidic χ 扭轉，主要是為了消掉 ff99 在長時間 RNA 模擬裡會長出來的 ladder-like 假結構。你的系統有 sgRNA，這條直接相關。
ff19SB + TIP3P 或 ff14SB + OPC 都不是被 validate 過的組合。力場與水模型是一整套 fit 出來的，蛋白參數在調的時候「看到的」就是某個特定的水。這跟你不會把某個泛函生出來的贗勢拿去配另一個泛函是同一類的顧慮 —— 不是一定會炸，但你失去了所有的驗證背書，而且錯的方向無法預測。
7. 水模型：便宜與正確之間的一個明確取捨。
- TIP3P（Jorgensen 1983，3-site）：最便宜、GPU 上最快、跟 ff14SB 配套。已知代價是動力學被系統性加速：自擴散係數是實驗值的兩倍以上，剪切黏度只有實驗值的三分之一上下（確切數值隨 cutoff、PME 設定與有限尺寸修正而異，要寫進 limitation 就自己引一次原始量測）。注意這些動力學缺陷不是原始論文的內容 —— 1983 那篇是 NPT Monte Carlo 研究，方法上根本不產生軌跡，只比了密度、位能與結構。TIP3P 的動力學問題是後來的 MD 文獻才揭露的，也就是說：它當初根本沒有對動力學 fit 過。
- OPC（Izadi 2014，3 電荷 / 4 點的剛性模型）：設計理念是方法論上的翻轉 —— 不先套用實驗的水分子幾何再調參數，而是除了對稱性之外不加任何幾何約束，直接最佳化點電荷分布去重現水的低階靜電多極矩。結果是一整組體相性質的重現度明顯優於常用剛性模型，小分子水合自由能的 RMSE 小於 1 kcal/mol。代價是多一個交互作用位點。搭配關係的方向要講對：OPC（Izadi 2014）比 ff19SB（Tian 2020）早，是後來的 ff19SB 被參數化成要搭 OPC（ff19SB 原文即以 OPC 為建議水模型），不是 OPC 為了配 ff19SB 而設計的；同理 ff14SB 當年是在 TIP3P 下參數化的。
決策準則： 如果你要的是接觸統計、界面埋藏面積、鹽橋數這類幾何量（也就是打分層用的聚合量），TIP3P 夠用，而且它讓你能多跑 replica。如果你要報時間尺度或速率，TIP3P 的時間軸不可信，那時候水模型就從「省錢選項」變成「結論的一部分」。這一條跟全站的立場是同一件事：先問穩不穩（replica 數），再問準不準（水模型）。
⚠ OPC 相對 TIP3P 的實際 wall-time 代價，這一頁沒有可靠數字可以給你（它跟原子數、GPU 型號、PME 設定都有關）。自己在你的系統上跑一次短 benchmark，不要照抄任何人的百分比。
8. 離子：中和加生理鹽度。
ionicStrength=0.15*molar, neutralize=True。Cas9 ＋ sgRNA ＋ dsDNA 是強負電系統，光是中和就要放進不少 Na⁺。0.15 M 是生理值，也是絕大多數文獻的預設，跟人家一致比較好比較。要留意但目前不用改的一件事： fixed-charge 一價離子在高濃度下有被報導過的 spurious clustering 傾向（Joung-Cheatham 參數改善了但沒有完全消除）。這在 0.15 M 到底有多嚴重，這一頁沒有查證過，若你之後要跑高鹽條件，先查一次文獻再投。
9. 盒型與 padding —— 這一節是本頁最值錢的地方。
先講清楚 padding 在物理上是什麼。它服務的是 minimum image convention：在週期性邊界下，任何一個原子都不應該「看見」自己的鏡像；否則你算的是一個人工的週期性晶體，不是一個溶液中的分子。這對高電荷的 protein–DNA 複合體特別敏感，因為長程靜電衰減得慢。定義要讀對：padding 是「solute 與它的週期鏡像之間至少相隔這麼多」，不是「solute 到盒壁有這麼厚」 —— 兩種讀法差一倍，而 OpenMM 做的是前者。
陷阱在於 OpenMM 怎麼把 padding 換算成盒子。 Modeller.addSolvent(padding=...) 的實作是：
其中 radius 是當下 Topology 裡所有原子的外接球半徑 —— 注意是所有原子，你保留下來的結晶水與已經放進去的離子也會一起把它撐大，不是只算蛋白與核酸 —— 而球心取的是包圍盒（bounding box）的中心，不是質心，原始碼裡就寫著 center = 0.5*(minRange+maxRange)、radius = max(norm(center-pos))。也就是說，它用的是分子最長的那條跨距，不是逐軸的尺寸。順帶注意盒寬只加一個 padding 不是每邊各一份：外接球直徑已經涵蓋分子在任何轉向下的最大跨距，padding 補的是球面到鏡像的那一段，只需要補一次。
對球狀蛋白這沒差。對 5F9R 這種又長又斜的完整 R-loop 複合體，這是災難。 實際數字：solute 最長軸只有約 13 nm，但立方盒建出來邊長是 19.87 nm、總共 762,027 個原子，是 Chai 版（273k 原子）的 2.8 倍。把公式倒推一次就懂了：padding 是 1.0 nm，所以 2 × radius = 18.87 nm —— 一個最長軸 13 nm 的分子，外接球直徑逼近 19 nm，因為它同時在另外兩軸也有可觀的厚度，外接球吃的是三個軸的平方和，不是最長那一軸。這就是「又長又斜」在成本上的代價。多出來的水全部塞在立方體的八個角落，那些位置的水分子在整條軌跡裡永遠不會碰到 solute，你只是付錢請它們待在那裡。
三個選項與它們的真實代價：
| 選項 | 原子數 | 省多少 | 代價 |
|---|---|---|---|
| A. 照跑立方盒 | 762,027 | 0% | 貴，而且貴在完全沒有用的角落水 |
| B. 菱形十二面體 | 535,265 | 30% | 零科學代價，改一行 boxShape='dodecahedron' |
| C. 對齊主軸的長方盒 | ~255,000 | ~67% | 必須假設分子在整條軌跡都不翻滾 |
這裡選 B。這個決定的邏輯值得記下來，因為它可以套用到 workflow 的每一個成本決策：
在「便宜」是硬約束的專案裡，最好的省錢決定，是那種不需要新增任何假設的省錢決定。
C 更省，但它把成本轉嫁成一個 reviewer 可以正面攻擊的假設 —— 87 ns 的自由 MD 沒有任何機制阻止分子轉向。B 是純幾何：最小影像距離一模一樣，翻滾一樣安全，只是不再付錢養角落的水。吞吐量因此快約 1.4 倍 —— 這個數字就是原子數的比（762,027 / 535,265 = 1.42），也就是「nonbonded 成本隨原子數走」的線性外推，不是獨立量測。要當成本估算用可以，要寫進 paper 就自己 benchmark 一次。
幾何上為什麼可以這樣：這對材料背景的人是熟面孔。 OpenMM 的 dodecahedron 盒向量是 (w, 0, 0)、(0, w, 0)、(w/2, w/2, w√2/2)，體積是 w³ × √2/2 ≈ 0.707 w³，而最小影像距離仍然是 w。這個形狀就是 FCC 晶格的 Wigner-Seitz cell（同一族裡 truncated octahedron 是 BCC 的，體積比約 0.77，所以菱形十二面體更省）。你只是把「怎麼用最小的體積填滿空間、又讓每個中心離彼此夠遠」這個你在畫 Brillouin zone 時就處理過的問題，從晶體搬到水盒。
10. 最小化與分階段平衡。
先破除一個 DFT 帶過來的直覺： MD 的 minimization 不是幾何最佳化。它的唯一目的是移除建模與加水造成的原子重疊，讓第一步積分不會炸。跑到「力收斂」沒有意義，因為你下一秒就要把系統加熱到 300 K 離開那個點。你的 build 腳本設 maxIterations=1500 就停，這是對的。（上機之後 production 腳本還會再做一次最小化，那次帶 k = 300 kcal/mol/Ų 的位置約束，目的是讓補出來的 loop 在 solute 被綁住的狀態下先消掉局部重疊 —— 一樣不是在找極小點。這兩次做的事不同，方法段要分開寫。）
分階段平衡才是真正在幹活的部分。你的專案用的是 Skeens 風格的階梯，把諧振位置約束加在 solute 重原子上（k 的單位是 kcal/mol/Ų）：
| 階段 | T (K) | k | 時間 | dt |
|---|---|---|---|---|
| heat_0_100 | 100 | 100 | 20 ps | 2 fs |
| npt_100_200 | 200 | 25 | 100 ps | 2 fs |
| npt_200_300 | 300 | 25 | 100 ps | 2 fs |
| release_300 | 300 | 5 | 200 ps | 2 fs |
| free_300 | 300 | 0 | 300 ps | 2 fs |
| density_300 | 300 | 0 | 5 ns | 4 fs |
為什麼要階梯而不是直接在 300 K 放開： 剛加完的水和剛補出來的 loop 帶著建模的殘餘應力，一次放開會用力推 solute。你會在前幾百 ps 看到一個很大的構形變化，然後很容易把它寫成「結構鬆弛到平衡構形」—— 那不是生物學，那是你的建模誤差在洩壓。約束階梯讓水先安頓下來、密度先穩定，再一級一級把 solute 放開。
HMR 與 timestep： 你用 LangevinMiddleIntegrator + constraints=HBonds + hydrogenMass=1.5*amu，最後一階跳到 4 fs。hydrogen mass repartitioning 的原理是把氫的質量從相鄰重原子那裡挪過來一些（總質量不變），壓低最高頻的 X–H 振動，讓積分器能用更長的步長。
第一，4 fs 的正當性主要不是靠 HMR，是靠積分器加 constraint。 OpenMM 使用手冊自己的說法是：在 HBonds constraint 之下，Verlet 類積分器約 2 fs、Langevin 類約 4 fs。你用的正是 Langevin 類，所以 4 fs 本身就落在手冊的範圍內，hydrogenMass=1.5 amu 是額外買的一點餘裕。它不是 Hopkins 那種把氫加重到幾個 amu 量級的正規 HMR —— 把「1.5 amu 對 4 fs 偏輕（社群慣例約 4 amu）」列為待補是對的，但理由是「HMR 的餘裕沒有真的買到」，不是「4 fs 本身違規」。方法段要寫的是積分器、constraint、hydrogen mass、步長四件事一起，不是「用了 HMR」四個字。Hopkins 原文採用的確切氫質量，這一頁沒有查證到，引用前自己核對一次原文。
第二，也是更少人講的：HMR 原始論文的驗證體系只有一個三殘基胜肽與 129 殘基的雞蛋白溶菌酶，完全沒有核酸。 你的系統是 535,265 原子的蛋白–RNA–DNA 三元複合體，遠在它的驗證範圍之外。這不代表 HMR 在這裡會錯，但代表你不能只引用那篇就當作背書。便宜的補救：拿一段短軌跡跑 2 fs 的對照，比較能量漂移與幾個關鍵幾何量的分布，證明 4 fs 這一組設定沒有改變你在乎的東西。
11. 序列化：讓叢集上零 build 風險。
最後一個實務動作。build 跑完之後把 system.xml（力場套用完的完整系統）與 state_min.xml（最小化後的座標）序列化存檔，叢集上的 production 只讀這兩個檔案，不再碰 PDBFixer 與力場。好處是：所有 replica 保證用的是完全同一個系統，差別只有初速度種子；而且叢集上不會因為套件版本差異在半夜掛掉。
06When什麼時候用
- 每一次你換了起始結構、力場、水模型、盒型、離子濃度 —— 那就是一個新系統，build 與平衡全部要重來，不能沿用。
- 要請人確認結構時 —— 先吐一個 solute-only 的 PDB（B-factor 標記 modelled / crystal），這比任何文字描述都有效。
- 要在叢集上跑多條 replica 之前 —— 在能除錯的機器上把 build 做完、序列化，再上傳。
- 要開始估成本時 —— build 完的
total atoms與 box vectors 就是你的成本單位。在投任何長跑之前先看這兩個數字。 - 要寫方法段時 —— 建系統的每一個決定都要在方法段出現一次。倒過來說：如果你寫不出來，代表你當初沒有做那個決定，是預設值幫你做的。
07When NOT什麼時候別用
-
不要在缺失殘基還沒被標記出來之前，就把系統送去跑生產。 為什麼會壞：補出來的 loop 沒有任何實驗座標約束，它在 300 K 下會自由甩動，而這些甩動會進入 RMSF、接觸數、界面統計、埋藏面積的每一個平均值。怎麼看出它壞了：把 RMSF 對殘基序號畫出來，看高峰是不是全部落在補出來的區段與鏈端 —— 你的 5′ 單股末端（res 1 到 11）是 6 到 8 Å，晶體本體只有 3.9 Å，這個對比一眼就看得出來。沒標記就沒辦法排除，沒排除就沒辦法報。
-
不要用預設的 protonation 猜測去跑跟活性中心有關的任何問題。 為什麼會壞：His 三態的重原子幾何幾乎相同，力場完全不會抱怨，但氫鍵網路、金屬配位、局部靜電全都不一樣。而且 pH 只在
addHydrogens那一瞬間存在，之後整條軌跡都無法改變 —— 這是單向門。 怎麼看出它壞了：檢查活性中心「該有的氫鍵在不在」，以及用 PROPKA/H++ 算出來的位移 pKa 跟你實際給的態一不一致。兩個工具給不同答案時，那個位點就是你方法段必須誠實交代的地方。 -
不要靠縮 padding 來省錢。 為什麼會壞：minimum image convention 一旦被破壞，solute 會透過 PME 感受到自己的週期鏡像，長程靜電被人工週期性汙染。高電荷的 protein–DNA 複合體是這件事最敏感的一類系統。怎麼看出它壞了：計算 solute 與其最近鏡像之間的最小距離隨時間的變化 —— 分子會翻滾，起始時滿足不代表全程滿足。如果這個距離曾經掉到接近你的 real-space cutoff（你用 1.0 nm）的量級，那段軌跡的靜電就有問題。正確的省錢方式是換盒型，不是縮 padding。
-
不要用對齊主軸的長方盒省錢，除非你能證明分子不翻滾。 為什麼會壞：長方盒的安全性建立在「solute 長軸持續對齊盒長軸」上，而自由 MD 裡沒有任何東西維持這件事。它一開始安全，然後在某個你不知道的時刻變得不安全。怎麼看出它壞了：追蹤 solute 慣性主軸與盒軸的夾角隨時間變化。這是本頁選菱形十二面體而不是長方盒的全部理由 —— 省 67% 卻要背一個會隨時間破掉的假設，比省 30% 但零假設差。
-
不要在起始結構本身可疑時，投入完整的多 replica 生產跑。 為什麼會壞：力場只會讓一個錯的起點在 300 K 下抖動，它沒有機制把結構拉回正確的 basin。怎麼看出它壞了：跨 replica 一致性會崩，而且從逐殘基層級先崩 —— Chai-1 起始是 0/8 共享、ρ ≈ −0.03，5F9R 晶體起始是 3/8 共享、Jaccard 0.48、Spearman +0.77。便宜的做法是先跑 2 條短的探路，但只用它來否證，不要用它來背書：像 Chai 那種 0/8、ρ ≈ 0 的災難級結果，2 條就足以叫停；反過來，2 條看起來一致不代表站得住 —— 你自己的數據就是反例，n = 2 時逐殘基是漂亮的 6/8，跑到 n = 4 掉到 3/8。小 n 會系統性高估一致性。
-
不要把「minimization 跑完了」當成系統健康的證明。 為什麼會壞：minimization 只保證沒有原子重疊，它不檢查 loop 建得對不對、protonation 選得對不對、離子有沒有卡在不該待的口袋、鏈的連接對不對。怎麼看出它壞了：看平衡階段的 box volume 與密度有沒有真的穩下來；看約束放開之後頭 300 ps 的 RMSD 有沒有跳一大階 —— 跳一大階通常代表你在洩建模應力，不是在觀察生物學。
-
不要宣稱只換亂數種子的 replica 涵蓋了建系統的不確定性。 為什麼會壞：你的 4 條 replica 共用同一個
system.xml、同一套補出來的 loop、同一組水與離子配置，差別只有初速度。它們是有效的統計獨立樣本（Lyapunov 不穩定性保證這件事），但它們對「loop 建錯了怎麼辦」這個問題完全沉默。怎麼看出它壞了：這個不會自己顯示出來，你只能主動說明。要涵蓋建系統的不確定性，得從不同的 build 重跑，那是另一組成本。
08Project Lens跟我的 project 多相關
為什麼是五星： 因為這是專案裡唯一一個「決定一次、影響全部」的環節，而且已經有兩個實測證據把它釘死。
證據一，起始結構決定了逐殘基分析能不能用。 Chai-1 起始：0/8 共享、ρ ≈ −0.03。5F9R 晶體起始：3/8 共享 [450, 695, 765]、pairwise Jaccard 0.48、Spearman +0.77。「元兇是起始結構，不是 metric」這句話是整個 reliability 分析最重要的結論，而它指向的是一個建系統決定，不是一個分析方法。
證據二，盒型決定了你能跑幾條 replica。 762,027 → 535,265 原子，省 30%，吞吐量約 1.4 倍（原子數比的線性外推）。在算力配額與磁碟空間都是硬約束的現實下，這 30% 直接兌換成 replica 數。而 replica 數決定了你的 CV 表可不可信 —— 你已經親眼看過小 n 陷阱：n = 2 時逐殘基一致性看起來是漂亮的 6/8，跑到 n = 4 掉到 3/8。 建系統省下來的算力，最後是在買統計力。
為什麼不是更低： 因為它便宜到接近免費。改盒型是一行程式。標記 modelled residue 是幾十行。請人確認結構只需要傳一個 PDB 檔。這些動作的成本，相對於它們能避免的浪費，比例是荒謬的。整張地圖上找不到第二個投報率這麼高的環節。
具體怎麼用：
- 把 build 當成一個要被 review 的產出。 你的 script 已經有正確的形狀：先吐 solute-only PDB（B-factor 標 modelled=1 / crystal=0）→ 等確認 → 才溶劑化 → 序列化成
system.xml+state_min.xml上叢集。把這個流程固定下來當標準。 - build 完必看三個數字：
total atoms、box vectors、以及 solute 的逐軸尺寸。前兩個 build report 會印，第三個要自己算（取座標的 min/max 相減，三行 numpy）。第一個決定成本，後兩個放在一起就是「有沒有在養角落水」的判準 —— 盒邊長明顯大於最長軸，就是外接球在替你決定成本。 - 維護一份「這個系統的建立決定」清單，跟軌跡放在一起。方法段直接從它長出來，而且下一個 session 接手時不用重新考古。
09Reviewer Thinkingreviewer 會問什麼
你補了多少殘基？補出來的區段有沒有進入你的統計？點開看參考答案
5F9R 補了 19 個 disordered residues：蛋白 N/C 端、sgRNA 3′ 端、non-target 股 5′ 的 11-nt。最尖銳的版本是「你的 non-target 股 30 nt 只有 19 個有實驗座標（63.3%），那你怎麼談 non-target 接觸？」能防守的回答需要三件東西：一張圖上把這些區段標成不同顏色；所有 domain 級平均值明確排除它們；以及一個可以直接給的 prepared_solute PDB（B-factor 欄位就是標記）。附帶兩個很有說服力的數字：補出來的區段 RMSF 是 6 到 8 Å 而晶體本體只有 3.9 Å；non-target 接觸的跨 replica CV 是 17.7%（marginal）而 target-DNA 接觸只有 9.4%（reproducible）。第二組對比要小心措辭 —— 你還沒做「排除 res 1–11 後重算」的對照，所以只能說它與「補出來的區段汙染了 non-target 統計」一致，不能說它證明了因果。主動端出這兩組對比、同時自己標明哪一組還只是假說，等於先替 reviewer 把問題問完了。
你的 protonation state 是怎麼決定的？H840 呢？點開看參考答案
最弱的答案是「用了 pH 7 的預設」。可以防守的答案是：一般殘基用標準 pKa 表，但活性中心與界面上的關鍵可解離基團另外用 PROPKA 或 H++ 檢查過，並且列出你實際採用的態。H840 是 HNH 的催化 His，如果你的結論跟 HNH 有關卻沒檢查過它，這一問會很難看。誠實的補充是：這個專案的 5F9R run 沒有加活性中心 Mg，focus 在 engagement 而非催化，所以 H840 的質子化態對現有結論的影響有限 —— 但要把這個 scope 講出來，不要讓它變成一個被抓到的疏漏。
你的盒子夠大嗎？minimum image 在整條軌跡都成立嗎？點開看參考答案
盒子邊長 19.87 nm 立方 / 換成菱形十二面體後最小影像距離不變、padding 1.0 nm。但「起始時成立」不等於「全程成立」，因為分子會翻滾。可防守的做法是實際算出 solute 與最近鏡像的最小距離對時間的曲線放進補充材料。這也順便把換盒型的決定變成優點而不是可疑點：同樣的最小影像距離、少 30% 的水。
為什麼選 TIP3P 而不是現在推薦的 OPC？點開看參考答案
誠實的回答有兩層。第一層是 scope：你報的是接觸數、埋藏面積、鹽橋這類幾何與統計量，不是溶劑動力學或速率，而 TIP3P 的主要缺陷（自擴散比實驗快兩倍以上、黏度只有實驗的三分之一上下）落在時間軸上 —— 要在方法段給數字就自己引一次原始量測，別引 Jorgensen 1983，那篇量不到動力學。第二層是成本：省下來的算力換成 replica 數，而 replica 數是這個 project 可信度的瓶頸。這是一個明說的取捨，不是一個沒想過的預設值 —— 這句話講不講得出來，決定了這一問是加分還是扣分。
你有沒有加活性中心的 Mg²⁺？沒有的話結論還成立嗎？點開看參考答案
沒有加，而且 5F9R 結構本身就沒有解出任何 Mg（非標準項目只有三個結晶用的 SO₄，加上 sgRNA 5′ 端的 GTP —— 後者是鏈的一部分，不是助劑）。這是刻意的 scope 決定：這一組 run 問的是 engagement 與界面接觸，不是催化。成立的部分是所有跟 protein–DNA 界面有關的量；不成立的部分是任何跟催化幾何、HNH 活性中心排列有關的推論。要碰後者就必須加 Mg，並且同時面對 fixed-charge 力場描述二價離子的已知問題。把界線劃出來，比含糊帶過安全得多。
你的 4 條 replica 是怎麼獨立的？點開看參考答案
同一個 system.xml 與 state_min.xml，只換亂數種子（初速度）。因為 MD 有 Lyapunov 不穩定性，幾十 ps 後軌跡就完全分開，所以這是有效的統計獨立取樣。但要主動說清楚它涵蓋什麼、不涵蓋什麼：它涵蓋「同一個系統的取樣不確定性」，不涵蓋建系統本身的不確定性（loop 怎麼建、水與離子怎麼放）。要涵蓋後者需要從不同 build 重跑。主動講出這個界線，比被 reviewer 問出來好得多。
10Common Mistakes最常犯的錯
- build 只跑到一半就去投 array job。 後果：4 條 job 全部在幾秒內崩掉（
FileNotFoundError: system.xml），而且因為死得太快，很容易誤以為是叢集問題。常見情境是 build script 設計成「先吐prepared_solute等確認 → 才溶劑化寫system.xml」，結果拿到確認檔就把 build 停了。正解：確認完把 build 一路跑完，看到system.xml存在再投。 - 在 login node 直接
python build.py。 後果：壓在登入節點上跑、只有 CPU、很慢、Ctrl-C 不一定停得下來，還可能被管理員砍。正解：任何要「算」的一律用排程系統送出（sbatch/srun），login node 只用來編輯、傳檔、看結果。 - build 完不看原子數就直接投長跑。 後果：762,027 個原子跑了才發現三分之一的水在角落發呆。正解：build report 一定印
total atoms與 box vectors，跟 solute 逐軸尺寸對照，差太多就檢查盒型。 - 把結晶助劑當成生物學配體留在系統裡。 後果：SO₄²⁻、PEG、甘油這些東西會佔住口袋、扭曲局部靜電，而它們只是結晶條件的產物。正解：
removeHeterogens，要留哪一個都必須有明確理由並寫進方法段。 - 反過來，把共價接在鏈上的非標準殘基當成助劑刪掉，還以為自己只是在清結晶條件。 後果：序列被悄悄改短，而方法段寫的是另一回事。5F9R 的 sgRNA 5′ GTP 就是這個例子（
removeHeterogens把 118 nt 變成 117 nt）。正解：build 前把HET與LINK記錄各印一遍，凡是有LINK接到聚合物鏈上的就不是助劑；刪了就照實寫。 - 沒處理 dangling 5′-phosphate。 後果：force field template match 直接失敗，build 掛掉 —— 這個算好的，因為它至少會報錯。正解：刪掉首殘基的
P/OP1/OP2變成 5′-OH，並在方法段註明。 - 混搭力場與水模型。 後果：
ff19SB + TIP3P、ff14SB + OPC都失去驗證背書，而且偏差方向無法預測。正解：整套用，換就整套換並重跑。 - 每條 replica 各自重跑 build，卻沒有記錄差異來源。 後果：取樣不確定性與建系統不確定性混在一起，CV 變得無法詮釋。正解：先固定一套 build 跑 replica 量取樣誤差；要量建系統誤差就另外設計一組實驗，兩者分開報。
11Further Reading讀哪幾篇
12Summary三句話
建系統是把整個位能面與邊界條件一次決定完的步驟，而它的錯誤幾乎都不會炸給你看，只會安靜地跑完並給出一個錯的結論。這個專案有兩個實測教訓：起始結構決定了逐殘基分析是隨機還是有訊號（Chai-1 起始 0/8、ρ ≈ −0.03，5F9R 晶體起始 3/8、Spearman +0.77），而盒型決定了成本（立方 762,027 原子 → 菱形十二面體 535,265 原子，省 30%、吞吐量約 1.4 倍，最小影像距離不變）。最值得記住的原則是：在便宜是硬約束的專案裡，最好的省錢決定是那種不需要新增任何假設的省錢決定。