5.3 RDF、氢键与取向分析

5.3 RDF、氢键与取向分析

本节摘要：局域结构统计三件套：径向分布函数 g(r) 给出"围绕参考原子，别人怎么堆"，氢键计数给出"谁和谁连着网络"，取向序参数给出"各向异性体系排得齐不齐"。本节逐项推导 g(r) 的归一化公式并用 Python 从零实现验证，讲氢键的几何判据与其解读边界，最后把序参数在脂膜双层上的用法说清。这三件套是 7.1 聚合物与 7.2 脂膜案例的全部观测量来源。

5.2 确认了整体稳定，本节把镜头推近到原子尺度的组织方式。承上：RMSD 看不出"水在聚合物周围怎么排"这类局域问题；启下：RDF 的峰位置与峰高是 7.1 判断聚合物-溶剂混溶性的硬指标，序参数是 7.2 脂膜的验收量。三个量共用同一套方法论：时间平均的空间统计——对每帧做空间统计、对帧取平均，正是 1.1 遍历性假设的兑现处。

RDF：归一化公式逐项拆

g(r) 的物理定义：距参考原子 r 处找到别的原子的概率密度，相对"完全均匀分布"的比值。写成可计算的形式：

g(r) = ( 1 / (N ρ) ) × ⟨ Σ_i Σ_{j≠i} δ(r − r_ij) / (4π r²) ⟩ / 归一化体积项

拆成操作步骤更好懂：对每帧、每个参考原子，把其他原子按距离落箱（bin）；每个箱内的计数除以"该箱在理想气体下的期望计数"（箱体积 4πr²Δr 乘数密度 ρ，再乘参考原子数 N）；最后对帧取平均。理想气体的 g(r) 处处为 1，液体在近距有峰（配位壳层）、远距趋于 1。Python 从零实现一个二维演示（二维便于画、便于逐行核对归一化）：

import numpy as np

rng = np.random.default_rng(7)
# 造一个 12x12 盒子的二维"液体"：90 个圆盘，含轻度短程排斥后的位形
box, n_part = 12.0, 90
pos = rng.uniform(0, box, size=(n_part, 2))
for _ in range(400):                          # 简单松弛，避免完全重叠
    d = pos[:, None, :] - pos[None, :, :]
    np.fill_diagonal(d[:, :, 0], np.inf)
    r = np.linalg.norm(d, axis=-1)
    close = np.argwhere(r < 0.9)
    if len(close) == 0:
        break
    i, j = close[0]
    push = (pos[i] - pos[j]) / (np.linalg.norm(pos[i] - pos[j]) + 1e-12) * 0.05
    pos[i] += push; pos[j] -= push

rho = n_part / (box ** 2)                     # 二维数密度
bins = np.linspace(0.05, 4.0, 80)
hist = np.zeros(len(bins) - 1)
d = pos[:, None, :] - pos[None, :, :]
dist = np.linalg.norm(d, axis=-1)
iu = np.triu_indices(n_part, k=1)             # 每对只数一次，最后乘 2 还原
hist, _ = np.histogram(dist[iu], bins=bins)
r_mid = 0.5 * (bins[:-1] + bins[1:])
dr = bins[1] - bins[0]
ideal = np.pi * ((r_mid + dr/2) ** 2 - (r_mid - dr/2) ** 2) * rho  # 环形面积 × 密度
g2d = 2 * hist / n_part / ideal
print("r      g(r)")
for r, g in list(zip(r_mid, g2d))[::6]:
    print(f"{r:5.2f}  {g:6.3f}")
​

核对口径：r 超过第一个峰后 g(r) 应在 1 附近摆动（找不到 1 说明归一化错了——十有八九是环面积公式或对数因子出错）。三维版本只是把环面积换成球壳面积 4πr²Δr，gmx rdf 的 -cn 选项还直接给出配位数（g(r) 从 0 积到第一谷），两者互验是分析报告的标准动作。

图 5-3 液态体系的 RDF 峰结构：从峰读组织

氢键：两要素判据与解读边界

GROMACS 的氢键判据有两个几何要素，缺一不可：给体-受体距离低于阈值（默认 0.35 nm）且氢-给体-受体夹角小于阈值（默认 30 度偏差）。只看距离会把大量"路过"的分子对误判成氢键——角度条件筛掉的正是不在成键取向上的。实操：

# 聚合物醚氧与水之间的氢键：两组各给一次
$ printf "PEO\nSOL\n" | gmx hbond -s md.tpr -f fix3.xtc -hbnum -o hbnum_peo_water.xvg
​

输出三类曲线要会读：总数目时间序列（涨落属正常，趋势性漂移要回查 4.3）、每个给体的平均键数（判断羟基是单供体还是双供体）、氢键寿命（由自相关推得，真正的动力学量）。解读边界同样要立住：平均氢键数不是强度的可靠指标——强氢键弱而少、弱氢键多而密，数目可以相同而网络强度天差地别；要谈强度得看 6.1 的自由能或光谱类观测量。脂膜案例（7.2）里还要防一个坑：头基间水桥氢键的计数对居中处理（5.1 配方三）极其敏感，居中做错，氢键统计直接翻车。

取向序参数：各向异性体系的体检表

RDF 与氢键对"排得齐不齐"不敏感，各向异性体系（脂膜、拉伸聚合物、界面水层）要用取向统计。脂膜的标准量是 C–H 键序参数：

S_CD = ⟨ (3 cos²θ − 1) / 2 ⟩

θ 是 C–H 键（或链轴）相对双层法向的夹角，尖括号对时间与等价原子取平均。取值 −0.5（与法向平行）到 1（严格垂直）到 0（各向同性）。流动相 DPPC 链中段的 S_CD 在 −0.2 上下，凝胶相绝对值明显更大——它直接给出链的"挺直程度"，是 7.2 判断相态、验证力场（与氘代实验的序参数直读对比）的核心量。gmx order 即为它而生；聚合物侧的对应物是沿拉伸方向或沿链轴的取向因子，七章案例里会用到简版。

演算：从 g(r) 积出配位数

配位数是 RDF 的第一个衍生量：把 g(r) 乘以数密度与球壳体积，从零积到第一谷，得到参考原子第一壳层里的平均邻居数。gmx rdf 的 -cn 给出它，但手算一遍才知道它对第一谷位置的敏感：

import numpy as np

# 沿用 5.3 的三维口径：给出一条模拟的 O-O RDF（示意数据）
r = np.linspace(0.20, 0.40, 200)
g = 1 + 2.4 * np.exp(-((r - 0.275) ** 2) / (2 * 0.018 ** 2))  # 第一峰
rho = 33.4                                                     # 水：分子每 nm3

r_max = 0.325                                                  # 第一谷位置（目测）
shell = 4 * np.pi * r ** 2 * (r[1] - r[0])
cn = np.sum(g * rho * shell[r <= r_max])
print(f"积分到 {r_max} nm 的配位数 ≈ {cn:.2f}")

# 敏感性：谷位置取 0.31 与 0.34，配位数差多少
for r_max in (0.31, 0.34):
    cn = np.sum(g * rho * shell[r <= r_max])
    print(f"谷位置 {r_max} nm -> 配位数 {cn:.2f}")
​

典型结果：水氧的第一壳层配位数约 4 到 5，谷位置从 0.31 挪到 0.34 会多吃进半个壳层的尾巴——报告配位数必须同时报告积分上限，只写一个数不给口径是这类分析最常见的打折。

三个高频疑问

RDF 该用哪两个组？ 按"化学问题"定参考-目标对，而不是默认全对全。问聚合物亲水性就用"醚氧-水氧"对；问离子缔合就用"阳离子-阴离子"对；全原子对全原子的 RDF 通常是一锅粥，峰互相叠看不出化学。一个实用习惯：每组 RDF 开跑前先在注释里写一句"这个峰若移动说明了什么"，分析时对着检验。

氢键寿命怎么算？ 用氢键存在的自相关函数：t 时刻成键、t+dt 仍成键的概率随 dt 衰减，拟合指数（或拉伸指数）得到寿命。gmx hbond 的 -life 选项直接给这个量。注意它对帧间隔敏感（5.3 前文说过的细帧距需求），寿命分析的轨迹帧距必须明显短于待测寿命。

序参数能用于聚合物吗？ 可以，且思路一样：选一个链轴或拉伸方向当"法向"，把 C-H 或链段矢量的 S_CD 公式换上去。熔体里各向同性、S 趋近 0；拉伸流动的体系 S 显著非零——它是"高分子取向度"的定量语言，与双折射实验同源可比。

速记：g(r) 是"相对均匀分布的堆积概率"，归一化核对口诀——远处必须趋于 1；峰位给距离、峰高给有序度、积到第一谷给配位数；氢键判据距离加角度两要素，数目不等于强度；序参数专治各向异性，脂膜链的挺直程度一测便知。三件套在手，第6章去够普通模拟够不着的量：自由能与罕见事件。