氫鍵分析Hydrogen bond analysis
方向性最強的非共價作用，也是最吃幾何判準的分析。
01一句話只有一句
用距離與角度門檻把一個連續的靜電作用切成有或沒有。
02Why should I care?我為什麼要讀這頁
因為這是整個分析章裡，你的 DFT 直覺最容易害你的一頁。
在 DFT 世界裡，氫鍵是一個有物理定義的東西：你可以在 QTAIM 裡找到鍵臨界點、算那裡的電子密度與 Laplacian；可以用 NBO 算 lone pair 到 σ* 的 donor–acceptor 穩定化能；可以做 counterpoise 修正過的二聚體掃描，得到一條連續的位能曲線。這些量都是連續的，而且都直接來自電子結構。
古典 MD 的 Hamiltonian 裡根本沒有「氫鍵」這一項。 力場只有 bonded 項、Lennard-Jones、和點電荷庫侖力。你在軌跡分析裡「數出來」的氫鍵，是你事後貼上去的標籤 —— 一個由你選的距離門檻與角度門檻定義出來的布林事件。力場不知道它的存在，也沒有針對它被參數化。
這件事有一個立即的後果：兩個軟體算同一條軌跡給你不同的氫鍵數，兩個都沒有錯，因為它們回答的是不同的問題。任何沒有寫出判準的氫鍵數字都是不可比、也不可重現的。
03What這是什麼
判準永遠是「一個距離 + 一個角度」，但是哪個距離、哪個角度，各家不同。三個原子是 donor（D，帶著氫的重原子）、hydrogen（H）、acceptor（A），常見的組合有：
- d(D···A) 加上 ∠(D–H···A)（頂點在氫上），要求接近 180°
- d(H···A) 加上 ∠(D–H···A)
- d(D···A) 加上 ∠(H–D···A)（頂點在 donor 上，也就是 D–H 向量與 D···A 向量的夾角），要求接近 0°
最後一種是 GROMACS 的寫法。它的 30° 與別人的 150° 方向上在要求同一件事（氫鍵要接近直線），但那是兩個不同頂點的角，數值不能直接相減互換 —— 三角形三個角加起來才是 180°，只有在 ∠(D···A–H) 趨近 0 的極限下兩者才互為補角。光是沒注意到這件事，就足以讓你把兩篇 paper 的數字比錯方向。
| 工具 | 距離判準（預設） | 角度判準（預設） | 額外行為 |
|---|---|---|---|
gmx hbond |
d(D···A) ≤ 0.35 nm | ∠(H–D···A) ≤ 30°（頂點在 donor） | 0.35 nm 來自 SPC 水 RDF 第一個極小值；donor／acceptor 預設元素只有 N 與 O（S 不算）；旗標在新版改名為 -hbr／-hba，舊版是 -r／-a |
| VMD HBonds | d(D···A) ≤ 3.0 Å | ∠(D–H···A) 偏離 180° ≤ 20°（即 ≥ 160°） | 預設偏嚴，常常一根氫鍵都畫不出來 |
MDAnalysis HydrogenBondAnalysis |
d(D···A) ≤ 3.0 Å | ∠(D–H···A) ≥ 150° | 需要拓樸裡有顯式氫 |
mdtraj baker_hubbard |
d(H···A) < 2.5 Å | ∠(D–H···A) > 120° | 預設 freq=0.1：佔據率低於 10% 的氫鍵直接不回報 |
注意 GROMACS 的 0.35 nm 與其他家的 3.0 Å 差了 17%。同樣是 donor–acceptor 距離，可及的球體積差了 (3.5/3.0)³ ≈ 1.6 倍。實際數量差沒有這麼誇張（介面上的原子對分布不是均勻的），但足以讓「氫鍵數 18 條 vs 24 條」這種比較完全失去意義。mdtraj 那個 freq=0.1 更陰險：它是靜默的，你不知道它幫你篩掉了什麼。
佔據率（occupancy） 是把時間維度收起來的方式：某一個特定的 donor–acceptor 對，在多少比例的幀裡滿足判準。它是介面分析裡真正該報的量，但它有一個在核酸介面上會爆炸的失效模式，見 When NOT 第 4 條。
04Why為什麼重要
物理上它是什麼。 一個具方向性、帶部分共價性質的靜電交互作用。氣相水二聚體的結合能是 5 kcal/mol 的量級；但在水裡，任何一個介面氫鍵的淨貢獻小得多，因為 donor 與 acceptor 原本都跟水成鍵，成鍵之前要先脫水。這就是為什麼「氫鍵越多結合越強」這種線性直覺在溶劑中幾乎總是錯的。
在古典 MD 裡它從哪裡來。 完全來自點電荷的空間排列加上 LJ 的排斥。方向性不是被寫進去的，是點電荷幾何的副產品。這意味著兩件事：氫鍵的品質完全繼承自力場的電荷模型（以及對核酸而言，骨架扭轉參數）；而力場調參的目標從來不是「氫鍵計數」，是密度、溶合自由能、螺旋參數這些巨觀量。你在量一個沒有被校準過的量。
在 Cas9 上它為什麼還是值得講。 因為 PAM 識別（PAM 是 Cas9 下刀前必須先認到的那三個鹼基 5′-NGG-3′）是整台機器上最像一道化學題的一段：R1333 與 R1335 各自對其中一個 guanine 做雙齒氫鍵。同一個結構還給了第二個同級的例子 —— phosphate lock loop（K1107 到 S1109）抓住 target 股上介於 PAM 與 +1 位之間的那個磷酸，把 +1 核苷酸轉出來與 guide RNA 配對，是局部解旋的物理起點。這一類「特定殘基的特定官能基、特定的氫鍵幾何、直接決定受質識別」的位置，有高解析度結構支持，是氫鍵分析最站得住的舞台。要在 paper 的 mechanism 段落把接觸數翻譯成化學語言，你需要這一層 —— 接觸數說「碰到了」，氫鍵說「以什麼方式碰」。
但要非常小心一個誘人的錯誤推論：鹼基配對本身就是氫鍵，所以 mismatch 就是氫鍵缺陷，所以 off-target 專一性可以用氫鍵計數來預測。 前兩步對，第三步錯得很徹底。高保真變體對 on/off-target 的結合親和力跟 WT 差不多，差別在 HNH 構形檢查點的動力學，不在配對能量。詳見 off-target 與 fidelity。
05How怎麼做
第 0 步：質子化狀態。這一步錯了，後面全部白算。
- His 有 HID / HIE / HIP 三種狀態。HID 與 HIE 的差別是那個質子在 δ 還是 ε 氮上 —— 同一個氮會因此從 donor 變成 acceptor。你數出來的那根氫鍵可能根本不存在。
- Asp / Glu 在低介電、被埋住或靠近其他負電荷的環境裡，pKa 可以顯著上移而實際帶質子。Cas9 的介面上到處是磷酸，這種環境不罕見。
- 工具：PROPKA、H++、PDB2PQR。不要把 PDB 直接丟進
tleap或pdb2gmx讓它按預設猜。 這是系統性誤差，跑再多 replica 也不會被平均掉。
第 1 步：接受一個事實 —— 氫的位置是建出來的，不是量出來的。
X-ray 在一般解析度下看不到氫（電子太少），只有超高解析度或中子繞射才行。所以你 pipeline 裡的每一個氫都是某個程式依幾何規則放上去、再被力場的鍵長鍵角約束住的。所有依賴氫位置的判準（d(H···A)、∠(D–H···A)）都繼承了那個建模步驟的不確定性。這就是為什麼有些人偏好只用重原子判準（d(D···A) 加上重原子幾何）—— 犧牲方向性，換取不依賴建模出來的座標。
第 2 步：算，並且分開報 direct 與 water-mediated。
gmx hbond -f traj.xtc -s topol.tpr -n index.ndx -num hbnum.xvg -r 0.35 -a 30，-r與-a就是判準，一定要在指令裡顯式寫出來，不要靠預設。⚠ 這兩個旗標在新版 GROMACS 已改名為-hbr／-hba，照抄舊指令會直接報錯 —— 先gmx hbond -h確認你那個版本叫什麼，並把版本號寫進 methods。對 535k 原子系統，先用 index 把 selection 縮到介面，否則會很慢。- MDAnalysis：
HydrogenBondAnalysis(u, donors_sel=..., hydrogens_sel=..., acceptors_sel=..., d_a_cutoff=3.0, d_h_a_angle_cutoff=150)，之後用.count_by_time()拿時間序列、.count_by_ids()拿逐對統計。 - mdtraj：
md.baker_hubbard(traj, freq=0.0)—— 把freq顯式設成 0，自己做佔據率統計，不要讓它替你篩。另有md.wernet_nilsson用的是錐狀判準，在水結構的文獻裡更常見。 - water bridge：要求同一個水分子同時對蛋白與 DNA 各成一根氫鍵。MDAnalysis 有
WaterBridgeAnalysis可以直接做。介面上這類作用不少，忽略它會讓你以為介面比實際上鬆。
第 3 步：敏感度掃描。 至少 3.0 / 3.5 Å × 130° / 150° 四組跑一次，證明結論不隨判準翻轉。這件事十分鐘，卻是這頁所有建議裡最能救你的一個。
第 4 步：報告。 每幀總數的時間序列（附跨 replica mean ± SD）＋ 殘基層級聚合佔據率表（附門檻）＋ 判準敏感度表。三者缺一，reviewer 都問得倒你。
06When什麼時候用
- 描述一個結構上已知的特定識別事件。 PAM 的 R1333 / R1335 對兩個 guanine、以及 phosphate lock loop（K1107 到 S1109）對 +1 磷酸，就是最典型的兩個例子。這種「有名有姓、有結構支持」的氫鍵是最站得住的用法。
- 解釋一個突變為什麼會壞掉。 把 donor 拿掉、把 acceptor 換位置，氫鍵分析講得出因果，接觸數講不出。Cas9 領域最乾淨的案例是 SpCas9-HF1：N497A／R661A／Q695A／Q926A 四個位點各自對 target 股磷酸骨架有一根直接氫鍵，突變等於逐根拆掉。但要分清楚 —— 「解釋這個變體為什麼變了」跟「預測哪個變體比較好」是兩件事，後者見 When NOT 第 2 與第 5 條。
- 當接觸數的機制註解。 打分交給接觸數，氫鍵負責把數字翻成化學。
- 蛋白內部與二級結構的穩定性檢查。 這是氫鍵最穩的舞台 —— backbone 的 α-helix / β-sheet 氫鍵被摺疊本身幾何約束住，佔據率高、判準敏感度低。
- 跟晶體結構的氫鍵清單做定性比對，當作「起始結構在 MD 裡有沒有鬆掉」的 sanity check。
07When NOT什麼時候別用
-
不要跨軟體或跨 paper 比較氫鍵數。 為什麼會壞：判準不同。GROMACS 用 d(D···A) ≤ 0.35 nm 加 ∠(H–D···A) ≤ 30°，VMD 用 3.0 Å 加 20°，mdtraj 的 Baker-Hubbard 用 d(H···A) < 2.5 Å 加 120° 而且預設只回報佔據率超過 10% 的。距離量的是不同的原子對、角度的頂點不同（不能直接換算，見 What）、還有靜默過濾。怎麼看出它壞了：把同一條軌跡餵給兩個工具各跑一次。這是五分鐘的實驗，而且幾乎沒人做；做過一次你就再也不會相信任何沒寫判準的氫鍵數字。
-
不要用氫鍵計數當變體排序分數。 為什麼會壞：它是布林事件的計數，門檻附近的 pair 會在 0 與 1 之間跳；而且完全沒有能量權重 —— 一根 2.8 Å 近乎直線的強氫鍵和一根 3.4 Å 彎曲的弱作用都算 1。怎麼看出它壞了：掃判準（3.0/3.5 Å × 130°/150°），看變體排序會不會翻。參考點：你自己的資料裡，定義寬鬆得多、完全不看角度的鹽橋聚合計數，跨 4 replica 的 CV 已經是 9.5%（勉強擠進 reproducible tier 的邊緣）。加上角度門檻只會更嚴、更容易跳。氫鍵這一項的 CV 目前尚未量測 —— 不要引用不存在的數字，但也不要假裝它會比 9.5% 好。
-
不要在質子化狀態沒定好前報氫鍵。 為什麼會壞：His 的兩個 tautomer 會把同一個氮從 donor 翻成 acceptor；被埋在介面裡、周圍全是磷酸的 Asp/Glu 若實際上是質子化的，你當它只能當 acceptor 就整組算反了。這是系統性偏差，不是統計雜訊，多跑 replica 沒有用。怎麼看出它壞了：對介面與催化中心附近的 His、Asp、Glu 跑 pKa 預測；然後換一組合理的替代質子化狀態重跑分析。若結論跟著變，那結論不是結論，是質子化假設的產物。
-
不要在蛋白–核酸介面逐對報氫鍵佔據率。 為什麼會壞：磷酸的兩個非橋接氧在力場裡幾乎等價，一個 Lys 或 Arg 側鏈可以在它們之間來回跳 —— 身分換了，物理上什麼都沒變；再加上介面暴露在溶劑裡，水隨時能插進來接手。結果是同一個殘基對的作用被切碎成一堆低佔據率的碎片，每一片看起來都「不重要」。怎麼看出它壞了：把同一個殘基對底下所有 donor–acceptor pair 的佔據率加總，跟其中最高的單一 pair 比。若殘基層級的總佔據率遠高於最高單一 pair，你看到的就是這個切碎效應。正解：報殘基層級聚合，或直接改用接觸數。
-
不要用氫鍵計數解釋 off-target 專一性。 為什麼會壞：這是這個 project 已經確立的結構性判斷 —— 高保真變體（SpCas9-HF1、eSpCas9(1.1)、HypaCas9）對 on-target 與 off-target 的結合親和力跟 WT 差不多，差別在 HNH 構形檢查點的動力學。氫鍵計數是平衡結構的描述，量不到跨越門檻有多難。怎麼看出它壞了：拿 SpCas9-HF1 當測試案例，然後看清楚你的數字在追什麼。它的四個突變 N497A／R661A／Q695A／Q926A 就是把四根對 target 股磷酸骨架的直接氫鍵拿掉的，所以你一定會算出氫鍵數變少 —— 那不是發現，那是你自己刪的。真正該被解釋的是它的表型：on-target 活性幾乎沒掉、off-target 掉到 GUIDE-seq 測不到，而這個 on／off 的鑑別力沒有寫在氫鍵數的差值裡。計數忠實地回報了你動的那隻手，卻回報不了那隻手造成的專一性。 更狠的旁證是文獻上的 ΔREC3 對照：刪掉 REC3 之後切割速率掉了約 1000 倍，對 on-target 的結合親和力卻接近 WT —— 結合面的幾何描述量與切割能力可以完全脫鉤（見 off-target 與 fidelity）。
-
不要忽略 water-mediated 氫鍵。 為什麼會壞：蛋白–DNA 介面有大量作用是透過一個橋接水完成的，直接 donor–acceptor 判準完全看不到它們。你會得到一個「介面很鬆」的錯誤印象，而且如果變體改變的恰好是水橋網路，你會完全漏掉效應。怎麼看出它壞了：分開統計 direct 與 water-mediated，看兩者比例；若 water-mediated 佔相當比例而你的敘事完全沒提，那個敘事是不完整的。
-
不要拿氫鍵佔據率當「介面穩定度」的單一指標。 為什麼會壞：佔據率低可能是判準太嚴、可能是水在競爭、可能是等價受體之間的切換、也可能是介面真的鬆 —— 四個原因給你同一個數字。怎麼看出它壞了：同時看接觸數。若接觸數穩（你的總 protein–DNA 接觸 CV 4.1%）而氫鍵佔據率碎，那介面是穩的，碎的是你的判準。
蛋白內部的氫鍵大多是 backbone 對 backbone，被二級結構的幾何鎖死，沒有水進得去，佔據率動輒 90% 以上、對判準幾乎不敏感。介面氫鍵剛好相反：side chain 對磷酸、暴露在溶劑中、受體有等價的多個氧、側鏈可自由旋轉、而且水隨時準備接手。同一個分析，在兩個地方的可靠度差一個等級。 用蛋白內部氫鍵的經驗去判斷介面氫鍵的可信度，是這一頁最容易犯、也最難自己發現的錯。
08Project Lens跟我的 project 多相關
為什麼不是更低（不是一星或兩星）： 因為它是把「聚合接觸數」翻譯成化學語言最直接的工具，而 methodology paper 的 mechanism 段落需要那個翻譯。PAM 識別的 R1333 / R1335 雙齒氫鍵、phosphate lock loop 對 +1 磷酸的鎖定、SpCas9-HF1 拆掉的那四根 target 股骨架氫鍵 —— 這三處是整台機器上最接近你熟悉世界的化學，都有結構支持，也都能接到你的 DFT 直覺。加上它的成本近乎為零（軌跡跑完順手就算），把它完全丟掉是浪費。
為什麼不是四星或五星： 三個獨立的理由，任何一個都足以把它擋在 Gate 3 打分軸之外。
- 判準依賴。 數字隨你選的距離與角度門檻移動，而門檻沒有第一原理依據。打分軸不能建立在一個你自己選的參數上，除非你能證明結論對它不敏感。
- 氫的位置是建模出來的。 依賴 ∠(D–H···A) 的判準，其精度上限受限於加氫程式與力場，而不是受限於你的取樣。
- 介面上被水與等價受體切碎。 見 When NOT 第 4 條。可靠性參考點是鹽橋（判準寬鬆得多、不看角度）的跨 replica CV 已經是 9.5%；氫鍵只會更嚴格、更容易跳。這一項的 CV 本身目前尚未量測，而這正是一個該補的 benchmark。
具體怎麼用：
- 機制敘事層與 QC 層，不進打分層。 打分交給聚合接觸數（CV 4.1%）。
- 若要把氫鍵寫進 paper，三件事缺一不可：殘基層級的聚合佔據率（不是逐對）、判準敏感度掃描、跨 replica 誤差。
- 一個具體可做的 benchmark 提案：拿現有的 4 條 replica，在兩三組判準下把介面氫鍵數的 CV 算出來，補進可靠性表。這是零額外 GPU 成本、可以直接寫進 methodology paper 的一格 —— 而且無論結果是綠是紅，都是有價值的資訊。目前這一格是空的。
09Reviewer Thinkingreviewer 會問什麼
你的氫鍵判準是什麼？換一組判準結論還一樣嗎？點開看參考答案
必須逐字寫出三件事：距離量的是哪兩個原子、角度量的是哪三個原子、以及門檻值。然後給敏感度掃描（例如 3.0/3.5 Å × 130°/150° 四組）。這一問跟接觸分析裡「你的 cutoff 怎麼選」是同一種責任，也跟你在 DFT 裡要交代 ENCUT 與 k-mesh 收斂測試是同一件事 —— 差別只在於，氫鍵判準沒有一個像能量收斂那樣的客觀終點，所以你只能靠證明「結論對它不敏感」。
你的質子化狀態是怎麼定的？介面上的 His 與 Asp 有沒有查過 pKa？點開看參考答案
如果答案是「用預設」，這一問會直接動搖整組結果，因為 tautomer 選錯會把 donor 與 acceptor 對調。可防守的回答：用 PROPKA 或 H++ 預測、逐一檢查介面與催化中心附近的可質子化殘基、並在 SI 列出最終指派。加分做法是對最不確定的那一兩個殘基做替代指派的敏感度測試。
晶體結構裡本來就沒有氫，你的氫是誰放的？點開看參考答案
這一問看起來刁鑽，其實是在測你知不知道自己在算什麼。誠實回答：X-ray 在一般解析度下看不到氫，所有氫位置都是建模程式依幾何規則加上去、再由力場約束的。因此依賴氫位置的角度判準有一層無法消除的建模不確定性；若要降低依賴，可以另外報一組只用重原子距離的判準做交叉檢查。
你有把 water-mediated 氫鍵算進去嗎？點開看參考答案
蛋白–DNA 介面的水橋不是少數例外。如果你的分析只看 direct donor–acceptor，就要明講這個限制，最好另外報一組 water bridge 統計。這一問答不好，reviewer 會合理懷疑你「介面比較鬆」的結論其實是分析看不到水。
為什麼你報氫鍵數而不是報結合自由能？點開看參考答案
不要辯稱氫鍵數等價於能量。正確回答：氫鍵數是幾何描述，不含去溶劑化代價、不含熵、也沒有能量權重，所以我們只用它做機制註解，量化排序交給可重現性更高的聚合接觸數。至於為什麼不做 MM-PBSA 或 FEP —— 那些排的是平衡 ΔG，而這個系統的專一性差別根本不在平衡 ΔG 上。
10Common Mistakes最常犯的錯
- 只寫「用 VMD 算氫鍵」，不寫判準。 後果：不可重現、不可比、reviewer 一問就沒得答。正解：明列距離量哪兩個原子、角度量哪三個原子、門檻各是多少。
- 用
md.baker_hubbard的預設卻不知道它內建freq=0.1。 後果：低佔據率的氫鍵被靜默丟掉，總數被系統性低估，而且介面上被切碎的那些恰好最容易低於 10%。正解：顯式設freq=0.0，自己做佔據率統計。 - 把 PDB 直接餵給
tleap/pdb2gmx，讓程式決定 His tautomer。 後果：系統性錯誤，跑再多 replica 都不會消失。正解：先跑 pKa 預測，逐一檢查介面與活性中心附近，並在 SI 列出指派。 - 逐對報介面氫鍵佔據率。 後果：被磷酸等價氧與側鏈翻轉切碎，每一根都「不重要」，讀者得到錯誤印象。正解：報殘基層級聚合佔據率。
- 把氫鍵數當結合強度的代理。 後果：忽略去溶劑化與熵，結論在水中幾乎必然不成立。正解：講清楚它是幾何描述；要能量就做 MM-PBSA / FEP，並且知道那些方法自己的限制。
- 拿別人 paper 的氫鍵數當對照組。 後果：判準不同、力場不同、水模型不同，數字不可比。正解：自己重跑一組對照，用同一套判準。
11Further Reading讀哪幾篇
12Summary三句話
氫鍵在古典 MD 裡不是 Hamiltonian 的一項，而是你事後用距離與角度門檻貼上去的布林標籤，所以不同軟體的預設判準會給出不同、且互相不可比的數字。它在蛋白內部最穩、在蛋白–核酸介面最不穩，因為磷酸的等價氧、可自由旋轉的側鏈、以及隨時準備接手的水，會把 per-pair 佔據率切成一堆碎片。把它放在機制敘事層、打分交給聚合接觸數，並且無論如何都要把判準與敏感度掃描寫出來。