分子动力学模拟氢键

我做过一个氢键网络分析项目，客户研究的是含羟基高分子在水溶液中的自组装，想知道氢键数量、几何分布和动态寿命。我这个项目让我把氢键的定义、判据和统计方法重新系统整理了一遍。

左侧图是氢键的几何判据。供体 D（donor）、氢 H、受体 A（acceptor）三点成键。我用的判据是：D-A 距离 r < 3.5 Å，D-H-A 夹角 θ > 120°。图中给了一个典型例子：r = 2.85 Å，θ ≈ 160°，这是一个非常强的近线性氢键。严格一点的判据会把距离收紧到 3.2 Å、角度收紧到 150°，但统计数量会少一些；宽松判据会包含更多弱相互作用，但可能把偶然靠近也当成氢键。

右侧图是氢键动力学。绿色实线是瞬时氢键数量随时间变化，平均值 mean = 24 条；蓝色虚线是氢键自相关函数 C(t)，描述一条氢键在 t 时刻仍然保持的概率。C(t) 从 1 开始衰减，前 5 ps 快速下降，然后趋于一个平台，说明氢键网络在不断断裂和重组。衰减时间 τ 可以从指数拟合得到，我这个体系 τ ≈ 2.5 ps，比纯水（约 1-2 ps）稍长，因为高分子链限制了水分子的运动。

分析流程：我用 GROMACS 的 gmx hbond 工具，输入轨迹和拓扑，自动识别供体和受体。输出包括氢键数量随时间变化、氢键寿命分布、氢键距离和角度分布。如果想更细，可以用 MDAnalysis 或 VMD 插件自定义几何判据。

踩坑经验：第一个是周期性边界：算 D-A 距离时必须用最小映像约定，否则跨盒子的氢键会被漏掉或算错。第二个是供体/受体列表：有些分子既有供体又有受体（如水、醇），要完整列出，否则数量会系统偏低。第三个是时间分辨率：如果轨迹输出间隔太大（比如 10 ps 一帧），会丢失快速断裂-重组事件，氢键寿命会被高估。我一般用 1-2 ps 输出一帧。

工程经验：氢键分析不能只看平均数量。寿命、角度分布、供体-受体配对类型，这些信息更能解释实验现象。比如客户后来发现，高分子链上的羟基与水的氢键寿命比水-水氢键长，说明高分子把水“锚定”了，这解释了体系黏度升高的实验现象。

更多 GROMACS 氢键分析与溶剂化结构的实战案例，我整理在 [分子动力学模拟氢键](https://www.keyanxueshu.com/category/md/) 栏目里。[分子动力学模拟氢键](https://www.keyanxueshu.com/)

氢键寿命的分析方法也值得关注。右侧自相关函数 C(t) 通常不是单指数衰减，而是多指数过程：快速衰减对应氢键断裂，慢速平台对应重新形成的氢键。我会用双指数拟合 C(t) = A1·exp(−t/τ1) + A2·exp(−t/τ2)，τ1 通常在 0.1-0.5 ps，τ2 在 1-10 ps。把这个拟合结果和客户讲清楚，比只说“平均寿命约 2.5 ps”更有信息量。另外，氢键角度分布也很重要：强氢键通常 θ>160°，弱氢键在 120°-150° 之间。如果体系中弱氢键比例高，说明氢键网络比较脆弱，可能对应较低的黏度或较快的质子转移。