8.7 均方位移、扩散、周期展开与可信误差
本课输入均为未执行的教学示例，数值是待测试的起点；图为原创示意，不是实际计算结果。使用授权 VASP 和 PAW 数据，替换所有占位符，记录程序版本并自行检验收敛。
8.7.1 模型、单位与记录
能量用 eV，长度用 Å，力用 eV/Å；1 kbar = 0.1 GPa。明确每原子、每分子、每原胞或每计算胞的归一化。记录 PAW 名称、版本、ZVAL、ENMAX 及允许的哈希；不分发 POTCAR。SCF 收敛只验证所选电子问题，目标物性的收敛仍需单独检查。
原创示意图；曲线与空白数据区用于说明或填写真实结果，不代表已运行计算。
8.7.2 完整案例：步骤、解释与检验
从跨边界到长时间扩散
原子跨越周期边界时，实际可能只走了一小步，但包裹坐标会突然跳到盒子另一侧。只比较首末包裹坐标又可能把长距离运动抹掉。因此扩散需要连续展开的轨迹。“让每个分子保持完整”与“让长期运动路径连续”不是完全相同的操作。
均方位移为 \(\operatorname{MSD}(\tau)=\left\langle|\mathbf r_i(t_0+\tau)-\mathbf r_i(t_0)|^2\right\rangle_{i,t_0}\)，这里 r 必须采用定义明确的展开坐标。各向同性的长期扩散区间满足 \(D=\frac{1}{2d}\frac{d\,\operatorname{MSD}}{d\tau}\)。体相三维用 d=3，面内二维用 d=2；受限法向不能当成自由扩散方向。各向异性体系应给方向分量或扩散张量。
实施步骤与误差
先核对采样时钟。每 q 步保存一帧时，帧间隔是 q×POTIM；ML 输出和重启可能改变实际文件节奏。正交固定晶胞可在相邻保存帧位移足够小的前提下逐步累加最小像位移。若保存间隔内跨越超过半盒长，周期像历史可能已无法唯一恢复，需要更密输出，不能事后猜测。
一般斜晶胞或变晶胞应使用经过测试的方法。晶胞形变和粒子迁移不能混为一谈，要明确实验室坐标、参考晶胞或其他位移约定。若去除整体漂移，说明使用哪个质量加权群组；直接去掉移动离子自身的质心运动可能删除真实集体输运。
对多个时间起点平均 MSD，记录每个滞后时间可用的起点数。避开短时弹道区、中间笼困平台和长滞后噪声区。局部对数斜率可辅助判断，但不是自动批准扩散拟合的证书。如果没有可信的扩散区，就应写“当前采样不足以确定 D”。
误差来自独立重复或足够长的时间块，并通过改变拟合区间检验。相邻 MSD 点共享大量坐标，不能把普通直线回归输出的标准误直接当成全部统计误差。还要检验尺寸、恒温器、密度和电子设置；有限尺寸修正有模型适用条件，不能对所有体系机械加同一个修正。
上面的伪代码假设坐标已正确展开、时间单位为 ps。三维各向同性才除以 6；Ų/ps 转 m²/s 乘 10⁻⁸。分子扩散要使用定义可靠的质心；质子交换或化学反应体系不能盲目固定分子标签。自扩散也不等于离子电导率。
报告与练习
报告选择对象、展开和漂移处理、粒子数、温度与系综、晶胞、步长、输出间隔、生产时长、拟合窗、独立重复数及误差方法。不能拼接彼此无关的重复轨迹成为一条连续路径，也不能把“没有观察到跳跃”写成“扩散严格为零”。
练习：先制作跨周期边界的合成轨迹作为软件单元测试，明确它不是材料模拟结果；比较包裹与展开 MSD。然后处理真实轨迹，改变拟合窗口，给出带条件的扩散值，或如实说明采样仍不足。
8.7.3 未执行输入与分析框架
本课输入均为未执行的教学示例，数值是待测试的起点；图为原创示意，不是实际计算结果。使用授权 VASP 和 PAW 数据，替换所有占位符，记录程序版本并自行检验收敛。
8.7.3.1 输入片段 1
# UNEXECUTED conceptual analysis; arrays supplied by a validated trajectory reader
# positions: continuous unwrapped coordinates, shape (n_frames,n_particles,3), Angstrom
# dt_ps: actual time between saved frames, in ps
for lag in selected_lags:
dr = positions[lag:] - positions[:-lag]
msd[lag] = mean(sum(dr * dr, axis=-1))
time_ps[lag] = lag * dt_ps
# Fit only a justified window, and repeat the entire fit per independent replica/block.
# D_A2_per_ps = fitted_slope / 6.0 # ONLY for isotropic 3D diffusion
# D_m2_per_s = D_A2_per_ps * 1e-8