MDAnalysis水分子位移与停留概率

本篇沿用水分子取向教程的同一份20分子、201连续帧真实水轨迹。MSD描述水氧移动多远；停留概率描述选中的水是否一直留在指定区域，两者不能互相替代。

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 在连续坐标上计算氧原子的MSD#

MeanSquareDisplacement要求完整水分子按OW HW1 HW2顺序排列，然后取每个三原子组的氧位移。它不会自动替你修复周期边界：本例已在输入准备阶段完成整分子和NoJump处理。

from waterdynamics import MeanSquareDisplacement
msd = MeanSquareDisplacement(water(), 'all', t0=0, tf=201, dtmax=20)
msd.run()

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
np.savetxt('water-msd.dat', np.c_[f['lag_ps'], f['water_MSD_A2']],
           header='lag_ps MSD_A2')
print('2 ps MSD:', f['water_MSD_A2'][-1], 'angstrom^2')
PYCODE

本次2 ps氧MSD为6.920107 Ų。滞后范围0.1–2 ps，每个滞后按该接口自身的间隔选择时间原点。没有选择长时线性区并检验统计误差，因此本例不把MSD/(6*t)直接当作收敛扩散系数。

3 明确区域，再计算停留#

这里定义一个具体示范区域：连续坐标中，实验室固定z方向的10 < z < 20 Å区间，统计水氧。它不是把每帧重新折回周期盒后的同名区间，也不是某个表面的吸附层。

from waterdynamics import SurvivalProbability
selection = 'name OW and prop z > 10 and prop z < 20'
continuous = SurvivalProbability(water(), selection)
continuous.run(tau_max=20, start=0, stop=201, step=1, intermittency=0)
returned = SurvivalProbability(water(), selection)
returned.run(tau_max=20, start=0, stop=201, step=1, intermittency=2)

stop=201不包括第201帧，实际使用0–200帧；该接口的stop遵循右端不包含的规则。intermittency=0要求每一帧都在区域内；intermittency=2允许中间最多缺席2个保存帧后返回，即本例最多0.2 ps，而不是允许整个分析时长之外无限期返回。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
time = np.r_[0, f['lag_ps']]
np.savetxt('water-survival.dat', np.c_[time, f['survival_gap0'],
           f['survival_gap2']], header='lag_ps continuous gap2frames')
print('2 ps continuous:', f['survival_gap0'][-1])
print('2 ps gap2frames:', f['survival_gap2'][-1])
PYCODE

2 ps处，连续停留为0.813584，允许缺席2帧为0.854409。两者的差异来自短暂跨出区域后返回；概率按每个起点当时区域内的分子数归一化，再对有效起点取平均，不能改成“总幸存数/所有起点总人数”后仍称同一定义。本例对每一个滞后的所有起点都做了独立布尔集合核对。

应用到孔道或表面时，先确定区域是否随界面移动、是否需要周期折回、按原子还是按残基统计，再调整选择语句。短轨迹的示范曲线没有拟合停留寿命，也不表示稀有脱附事件已经采样充分。

参考：停留概率API及帧数约定。