MDAnalysis指定原子对与氢键存活

总体O—O RDF可以描述同类原子的平均分布，指定原子对的RDF则跟踪“这两个原子”怎样相遇。氢键存活又是第三个问题：初始存在的某条氢键在后续是否仍存在。本文分别演示，避免把不同统计量混在一起。

1 下载和运行#

下载水、泛素与珠簧聚合物的实际轨迹、分析脚本和完整数据，保存为mdanalysis-water-polymer-wave330-v2.zip：

module load pymatgen/2026.9.24
python -m zipfile -e mdanalysis-water-polymer-wave330-v2.zip .
cd mdanalysis-water-polymer
bash run.sh
source env.sh

run.sh只对已完成的真实轨迹做后处理，生成results/和六张plots/图片，不启动新的分子动力学。请在新的解压目录运行，避免覆盖自己的输出。MDAnalysis为2.9.0，环境通过现有qvasp 2.26的Python选择，没有调用qvasp或VASP；其他机器可通过ANALYSIS_PYTHON指定兼容环境运行run.sh。

水体系来自已有GROMACS水模拟，884个TIP3P水、2652原子、1001帧，时间0–100 ps。水分析每20帧取一帧，共51帧。轨迹为NPT，实际盒长约29.87494–30.29525 Å，不能当作固定体积轨迹交给只支持固定正交盒的LinearDensity。本例使用能明确处理所需空间网格或平均体积的工具；短轨迹统计仍要检查采样收敛。

2 InterRDF_s：保留每个指定配对#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.rdf import InterRDF_s
u = mda.Universe('inputs/water/water.gro', 'inputs/water/water.xtc')
o = u.select_atoms('name OW')
a = InterRDF_s(u, [[o[:2], o[2:4]]], nbins=50, range=(0,10), norm='rdf')
a.run(stop=1001, step=20)
print(a.results.rdf[0].shape)
print(a.get_cdf()[0].shape)
PYCODE

每组结果为(2,2,50)：OW1/OW2分别与OW3/OW4组成四对，保留各自50个径向区间。get_cdf()输出每个指定配对的累计平均出现次数；它不是对全部884个氧求和后的配位数。

norm='rdf'按每个配对的球壳体积和样本平均盒体积归一化。只有四个指定配对、51帧时曲线会很稀疏，不能用它替代充分统计的整体水RDF。整体O—O RDF流程见基础轨迹分析。

3 HydrogenBondAutoCorrel：先固定判据与取样#

本例选前五个水的HW1作为氢及对应OW作为供体，全部884个OW作为受体。这里只跟踪这个明确子集，不代表同时统计所有水的两颗氢。原生工具的距离阈值是氢—受体H—A < 3.0 Å，角度为D—H—A > 130°；不要把它误读成供体—受体距离。

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel import HydrogenBondAutoCorrel
u = mda.Universe('inputs/water/water.gro', 'inputs/water/water.xtc')
h = u.select_atoms('name HW1')[:5]
d = u.select_atoms('name OW')[:5]
a = u.select_atoms('name OW')
for kind in ('continuous', 'intermittent'):
    c = HydrogenBondAutoCorrel(u, hydrogens=h, donors=d, acceptors=a,
        bond_type=kind, dist_crit=3.0, angle_crit=130.0,
        sample_time=2.1, nruns=3, nsamples=10, pbc=True)
    c.run()
    print(kind)
    print(c.solution['time'])
    print(c.solution['results'])
PYCODE

三个时间起点分别建立自己的初始氢键集合，再平均归一化结果。受轨迹保存间隔和原生整数帧截断影响，实际横轴以solution['time']为准，不能只按sample_time想当然画成0–2.1 ps。

连续定义中，一条键在某次采样中断开后就不再算作存活；间歇定义允许同一初始配对后来重新满足判据。因此本例间歇曲线可以高于连续曲线，二者回答的问题不同。本例没有设置time_cut，也没有进行多指数寿命拟合。

本例实际存活概率的三个取样点如下：

0.1 ps保存间隔无法观测更短时间尺度的瞬间断裂和重连。只有五个HW1、三个起点和短窗口时，曲线适合学习操作，不适合直接作为收敛的水氢键寿命。正式比较不同材料或条件，应保持判据、原子子集、保存间隔及连续/间歇定义一致，并增加独立采样。

4 导出完整配对和曲线#

python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
print(z['site_RDF_count'].shape)
np.savetxt('my-hbond-survival.csv',
    np.column_stack([z['hbond_time_ps'],z['hbond_continuous'],z['hbond_intermittent']]),
    delimiter=',', header='lag_ps,continuous,intermittent', comments='')
PYCODE

计数、归一化RDF、CDF及两种存活曲线都保留完整数组。水的电荷、连接关系和密度见上一篇。