MDAnalysis的RDF、氢键、MSD与键角

下载本篇输入、完整脚本、数据与高清图。解压后在新目录运行，reference保存鸢算科技实际运行的参考结果。

1 运行方式#

先下载并按轨迹读取教程运行完整脚本。各分析输出独立CSV，配套plots.py会生成高清图。以下输入使用同一条884个水分子的100 ps轨迹。

2 O–O径向分布函数#

from MDAnalysis.analysis.rdf import InterRDF
o = u.select_atoms('name OW')
r = InterRDF(o, o, nbins=100, range=(0,10),
             exclusion_block=(1,1)).run(step=10)
print(r.results.bins, r.results.rdf)

每10帧取一帧，共101帧；距离0–10 Å，100个区间。exclusion_block=(1,1)排除原子与自身的配对。本次第一个峰的区间中心为 2.75 Å，峰高 2.64455。

截断应小于最短盒长的一半；峰位置还受bin宽、温度、力场与采样长度影响。此例演示提取，不代表已完成尺寸和采样收敛。

3 氢键数量#

因为GRO没有电荷和键，本例显式给出供体、氢和受体名称：

from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis
h = HydrogenBondAnalysis(u, donors_sel='name OW',
    hydrogens_sel='name HW1 HW2', acceptors_sel='name OW',
    d_h_cutoff=1.2, d_a_cutoff=3.5, d_h_a_angle_cutoff=150,
    update_selections=False).run(step=50)
print(h.count_by_time())

供体–氢配对距离≤1.2 Å，供体–受体距离≤3.5 Å，D–H–A角≥150°。21个采样帧平均有 1142.238条氢键，除以884得到每个水平均捐出 1.29212条。

这个数按捐出键计，每条氢键只计一次；若按每个分子“参与”的氢键计数，会采用另一种归一化。比较文献时先统一判据和计数方式。

4 展开轨迹，再计算MSD#

原轨迹盒长随时间变化。先按顺序遍历全部帧，应用支持变胞的NoJump，再抽取20个氧原子计算MSD：

from MDAnalysis.transformations.nojump import NoJump
u.trajectory.add_transformations(NoJump())
coords = []
oxygen = u.select_atoms('name OW')[:20]
for ts in u.trajectory:
    coords.append(oxygen.positions.copy())

不能先大步抽帧再展开，否则中间跨界信息可能丢失。NoJump要求相邻存储帧之间的运动可由半盒范围判断；轨迹输出过稀时仍可能产生歧义。脚本将展开结果装入内存Universe，再执行：

from MDAnalysis.analysis.msd import EinsteinMSD
m = EinsteinMSD(v, select='all', msd_type='xyz', fft=False).run()
msd = m.results.timeseries

这里v是脚本创建的20氧原子展开轨迹。0.1、1、10 ps滞后时间的MSD分别为 0.260028、3.508543、34.735633 Ų；与直接对所有时间起点求平均的结果一致。

长滞后时间可用时间起点更少，噪声通常更大。本例没有给出收敛扩散系数；正式研究应增加统计量，检查线性时间窗、误差、盒尺寸与恒温恒压设置。

5 周期边界下的键角#

from MDAnalysis.lib.distances import calc_angles
angle = calc_angles(H1_positions, O_positions, H2_positions,
                    box=u.dimensions)

函数返回弧度，乘180/π才是度。该水分子的101帧平均H–O–H角为 104.54323°，与刚性水几何相符。

参考：MDAnalysis分析类与results接口、NoJump与变胞轨迹。