VASP液体自扩散系数
从解包裹坐标和多时间起点MSD，估算Na液体的自扩散系数。
阅读本篇目录
自扩散系数描述单个原子在长时间尺度上的随机迁移。不能只看某一步原子移动了多远；通常先计算均方位移MSD，再寻找近似线性时间区间。
1. 取得完整生产轨迹#
使用Na液体结构教程的三条2 ps NVE生产轨迹。体积和原子数相同，温度用生产阶段的实际平均值；输入文件名中的600、900、1200只是支路标识，不作为最终温度测量值。
短NVT阶段后，两个较低温支路的实际温度可能接近。若目标是严格比较指定温度，应延长热化并重新检查能量和结构分布，而不是把名义温度填进结果表。
2. 先修复周期跳跃，再去除整体漂移#
XDATCAR中的原子穿过晶胞边界后会回到另一侧。直接用包裹坐标作位移，会产生错误跳跃。先按相邻帧的最小分数坐标变化累积成连续坐标；再减去体系质心移动。
本例固定正交晶胞、相邻步位移远小于半盒长，适合这种解包裹方法。可变晶胞轨迹需要额外处理晶胞形变，不能照搬固定盒算法。
3. 多时间起点平均MSD#
对多个起始时刻和全部同种原子求平均：MSD(t)=⟨|rᵢ(t₀+t)−rᵢ(t₀)|²⟩。三维各向同性扩散的长时间关系为MSD≈6Dt+C。
本例用0.3–0.9 ps区间拟合斜率，D=斜率/6。Ų/ps转换为m²/s需要乘10⁻⁸。二维扩散的相应系数是4，不能仍除以6。
4. 本次短轨迹估计#
| 实测平均T / K | 全轨迹拟合D / m²·s⁻¹ | 线性拟合R² | 前/后半段斜率÷6 / m²·s⁻¹ |
|---|---|---|---|
| 793.7 | 不报告 | 0.8027 | 2.931e-09 / 7.819e-10 |
| 801.8 | 不报告 | 0.3604 | 4.409e-11 / -1.687e-09 |
| 1050.0 | 1.277e-08 | 0.9963 | 1.049e-08 / 1.680e-08 |
仅1050.1 K支路给出正且近似线性的扩散估计。793.7和801.8 K支路在所选区间的斜率分别为负，不能解释成负扩散系数；它们反映短轨迹中的振荡、回返或采样不足，本例不报告这两条的D，也不据此确认液相平衡。应延长热化和生产、检查结构与MSD是否进入扩散区间。
前后半段各只有1 ps，分别使用0.15–0.45 ps滞后区间；表中列的是拟合斜率除以6，负值同样不作为物理扩散系数，，用于显示对轨迹片段和拟合窗口的敏感性，不是严格统计误差棒。高R²仅说明所选区间接近直线，不能证明已经达到真正长时间扩散极限。
5. 扩散活化能需要怎样的数据#
常见做法是多个温度下得到可靠D，再拟合lnD对1/T；若服从Arrhenius关系，斜率为−Ea/kB。必须使用实测温度、保持可比较的热力学条件，并检查每个D的统计误差。
本次16原子、每条2 ps，两个低温支路温度接近且没有可报告的扩散估计，不据此给出可靠的材料扩散活化能。固定体积下的温度扫描也不等同于恒压实验。需要材料结论时，应先扩大体系、延长轨迹和热化，再补充间隔合适的温度点与独立重复。
输入、脚本与数据#
本例VASP 6.6.1 P100 Gamma版，PBE PAW64 Na（1个价电子），EDIFF=1E-6；构建使用NVHPC 25.3、CUDA 12.8、MKL 2023.2及配套MPI。只有Na16小体系、固定密度和短轨迹，不代表任意液体的通用收敛设置。下载不含POTCAR；按README顺序提交热化和生产任务，结束后再分析自己的完整轨迹。