MDAnalysis水桥网络与几何分析

氢键分析回答哪些原子形成氢键；水桥分析继续回答：两个选区是否经一个或多个水分子连接。本篇从准备文件、设置选区开始，到提取路径、统计数量和查看真实结构，演示原生WaterBridgeAnalysis的用法。

1 下载与运行#

下载水桥分析操作包，保存为mdanalysis-water-bridges.zip。包内提供已经完成计算的真实GRO/XTC水轨迹、完整分析脚本和参考数据，不需要重新做分子动力学模拟。

在集群解压后，普通批作业可提交sbatch job.slurm。以下直接运行命令用于已经分配的CPU资源内；大轨迹分析请提交作业处理。

python3 -m zipfile -e mdanalysis-water-bridges.zip .
cd mdanalysis-water-bridges
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

依次完成环境检查、轨迹分析、原样示例对照与绘图，结果在results/，两张图片在plots/。环境为MDAnalysis2.9.0、NumPy1.26.4、SciPy1.13.1、Python3.11；通过集群已有qvasp-2.26的Python环境调用MDAnalysis，本篇不运行qvasp或VASP。脚本设置1个线程。

2 读取真实轨迹，明确两个端点选区#

本例为884个TIP3P水分子、2652个原子，原轨迹0–100 ps共1001帧。本次每20帧取一帧，即每2 ps抽样一次，共51帧。长度使用Å、时间使用ps，盒尺寸取各帧实际值。

import MDAnalysis as mda
u = mda.Universe('input-data/water.gro', 'input-data/water.xtc')
print(len(u.atoms), len(u.residues), len(u.trajectory))

selection1 = 'resid 1:10'
selection2 = 'resid 11:20'
water = 'resname SOL and not resid 1:20'

端点是两组各10个固定编号的水分子，其余水可以充当中间连接。这样可以直接利用已有液态水轨迹理解算法；这里没有蛋白或配体，不把这两个选区解释为蛋白–配体结合位点。

实际研究蛋白–配体水桥时，改成对应残基与配体选区，并把中间水设成真实水残基；同时检查氢原子名称、供体/受体规则和拓扑。两个端点选区和中间水应按模型含义合理区分。

3 显式设置距离与角度#

from MDAnalysis.analysis.hydrogenbonds import WaterBridgeAnalysis

w = WaterBridgeAnalysis(
    u, selection1, selection2,
    water_selection=water,
    order=2,
    distance_type='heavy', distance=3.5, angle=150,
    pbc=True,
    forcefield='other', donors=['OW'], acceptors=['OW'],
    filter_first=True,
    update_selection=True, update_water_selection=True)
w.run(start=0, stop=1001, step=20)

这里明确采用供体–受体重原子距离≤3.5 Å，D–H–A角大于150°。本版本原生水桥角度实现采用严格大于；不要把输入距离与默认的氢–受体距离混为一谈。供体和受体重原子名称均为OW，氢按同一残基内的名称及与供体的距离识别。GRO不含完整力场拓扑，示例不凭空补出电荷。

order表示允许的最多中间水分子数量：

所以order=2的输出不能全部叫“二水桥”。一条路径含n条氢键时，中间水数量为n−1；附件逐条遍历原生网络，按这个定义分开统计。

4 实际逐帧计数#

同一对端点可能有多条连接路径；一条路径也可能在多个采样帧出现。下面统计的是每帧路径数再求和，不是唯一水分子数量，也不是不同事件的去重数。

只查看输出，不启动新模拟：

python3 -c "import json; v=json.load(open('results/VALIDATION.json')); print(v['summary']); print(v['results'][0])"

对自己的体系，常用的另一种指标是“某一对端点在多少帧存在至少一条水桥”。这需要先按端点和帧去重，再除以分析帧数；不能直接拿上表路径总数除以51当作存在概率。

5 查看网络里的真实水桥#

原生w.results.network按帧保存氢键连接网络。网络键中的原子编号从0开始；某一侧是供体氢时，还附带对应重原子的编号。读取路径时应保留供体氢、供体重原子、受体重原子三者，避免只按两个水残基编号去重而丢掉不同氢的连接。

本例第一帧的一条实际二水桥为：

SOL4 → SOL149 → SOL397 → SOL20

结构图只沿所选路径选择连续的周期镜像，没有优化或重造连接结构。蓝色和黄色为两个端点，红色为中间水氧，浅色为氢，虚线连接实际氢键。

6 输出文件与应用边界#

NATIVE-PATHS.json保存各帧完整路径和每条键的距离/角度；FIELDS.npz保存真实坐标、各帧盒尺寸、采样时间和分类计数；INDEPENDENT-PATHS.json保存按同一几何定义重新枚举得到的连接，便于理解输出含义。

source env.sh
python3 -c "import numpy as np; d=np.load('results/FIELDS.npz'); print('times/ps:', d['time_ps']); print('counts:', d['paths_by_order'])"
ls -lh plots/water-bridge-counts.png plots/water-bridge-geometry.png

这个例子用于演示选区、几何判据、周期边界和路径读取。2 ps抽样无法解析更短的断裂/重建过程，不能从这些离散计数直接得到氢键寿命、动力学速率或结合自由能。提高order会显著增加搜索量；先用低阶与少量帧确认选区，再处理正式轨迹。

参考：MDAnalysis2.9.0水桥分析与网络接口。