MDAnalysis水分子取向与角分布

水分子取向关联回答“经过一段时间后，方向还记得多少”；角分布回答“水分子更倾向朝哪个方向”。两者不是同一个量。本例直接分析鸢算科技集群已有TIP3P水轨迹，选前20个完整水分子、连续201帧（0–20 ps，保存间隔0.1 ps），没有另外生成分子动力学。

1 下载、准备和运行#

下载水分子动力学与蛋白路径比较完整示例，保存为water-paths.zip。该包使用MDAnalysis 2.9.0，独立扩展waterdynamics 1.2.0和PathSimAnalysis 1.2.1；两项扩展只安装在示例目录的.deps中。

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

results/FULL-FIELDS.npz是完整数值，results/VALIDATION.json列出逐项核对，plots/保存4张图片。脚本发现已有同名输出目录会停止，第二次运行可以用bash run.sh --output results-new指定新目录；不要把旧结果当成新轨迹。

2 先让水分子完整，再处理跨周期边界#

每个水分子的原子顺序为OW HW1 HW2，必须选择完整的三个原子；不能只选氧来运行取向分析。准备数据时先根据真实的O–H连接做unwrap，再对连续帧做NoJump，最后才截取这20个分子。配套compact-inputs/water.npz保留了处理后的全部坐标、实际盒子和时间。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
w = np.load('compact-inputs/water.npz')
print('coordinates:', w['positions_A'].shape)
print('time:', w['time_ps'][0], w['time_ps'][-1])
print('first box:', w['dimensions_A_degree'][0])
PYCODE

输出坐标形状为(201, 60, 3)。本例是有体积波动的NPT轨迹，不能把盒长固定后再把结果当作同一轨迹。

3 用原生接口计算OH、HH和偶极方向#

analyze.py中的water()用包内坐标建立Universe，补入同一TIP3P模型的质量、电荷和O–H连接。关键调用如下：

from waterdynamics import WaterOrientationalRelaxation
analysis = WaterOrientationalRelaxation(
    water(), 'all', t0=0, tf=201, dtmax=20, order=2
)
analysis.run()
values = analysis.timeseries  # 20行，依次为OH、HH、偶极方向

order=1计算方向点积；order=2计算(3*cos(theta)**2-1)/2。dtmax=20是20帧，即本例2 ps，不是20 ps。该版本WOR按每个滞后帧数移动时间原点：例如滞后2帧比较0→2、2→4等；它并非所有相邻时间原点的普通滑动平均。配套文件也保存了每个滞后的实际原点数。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
print('lag (ps):', f['lag_ps'][-1])
print('P1 OH, HH, dipole:', f['WOR_P1'][-1])
print('P2 OH, HH, dipole:', f['WOR_P2'][-1])
PYCODE

2 ps处，P1三个方向约为0.384962、0.367707、0.354953，P2约为0.129059、0.102035、0.104082。这是20个分子的短时间示范；没有进行指数拟合，也不据此报告收敛的水分子转动弛豫时间。

4 角分布的横坐标是cos(theta)#

from waterdynamics import AngularDistribution
angular = AngularDistribution(water(), 'all', bins=20, axis='z')
angular.run()

axis='z'指定与z轴比较，也可改为x或y。这里统计每个方向4020个水分子—帧样本（20×201）。角度均匀和余弦均匀是两回事：各向同性方向的cos(theta)概率密度是0.5。

本版本angular.graph是供阅读的文本，边界值会被转成简短字符串。配套脚本读取原始数值直方图保存完整边界，用sum(density*bin_width)检查归一化，不从截短后的文本反推精确分箱。

在自己的体系里，先确认水模型、原子命名、完整分子选择和边界处理，再调整分析帧数。界面水建议以界面法向为轴，不能把本例体相水的小样本波动解释为界面定向。

参考：waterdynamics原生API。