pymatgen周期轨迹与MSD

下载输入、实测结果和完整脚本。在解压目录运行 bash run.sh，结果写入 work/；reference/保留原始结构与CSV，脚本不会自动提交新的电子结构计算。

1 读取实际VASP轨迹#

module load pymatgen/2026.9.24
python - <<'READ'
from pymatgen.io.vasp.outputs import Xdatcar
frames = Xdatcar('inputs/na-md/XDATCAR').structures
print('帧数/原子数:', len(frames), len(frames[0]))
READ

本例16个Na原子、1000帧，实际INCAR中POTIM=2 fs，帧间隔2 fs，约2 ps生产轨迹。输入目录名中的1200是历史支路标识，不作为实测平均温度；严格按指定温度研究时应重新检查热化和温度统计。

2 先修复周期跳跃#

将下面内容保存为 msd.py：

import numpy as np
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core.trajectory import Trajectory
frames = Xdatcar('inputs/na-md/XDATCAR').structures
coords = np.array([s.frac_coords for s in frames])
t = Trajectory([str(x.specie) for x in frames[0]], coords,
               lattice=frames[0].lattice, time_step=2)
t.to_displacements()
u = np.cumsum(t.coords, axis=0) @ frames[0].lattice.matrix
u -= u.mean(axis=1, keepdims=True)
lags = np.arange(1,301)
msd = [np.mean(np.sum((u[k:]-u[:-k])**2, axis=2)) for k in lags]
np.savetxt('na-msd.csv', np.c_[lags*0.002,msd], delimiter=',',
           header='lag_ps,MSD_A2', comments='')

python msd.py

to_displacements()先按周期边界取得相邻帧位移，累积后得到连续坐标。原子质量相同，本例用平均原子位移去除整体平移；多组分体系要按实际质量处理质心漂移。

这里假设晶胞固定、相邻帧的分数位移不超过半盒长。轨迹输出过稀时不能唯一恢复跨越次数；可变晶胞也需要专门处理晶胞形变，不能照搬固定盒代码。

3 多时间起点平均#

每个滞后时间k，代码使用所有允许的帧对，而非只相对第0帧求位移。最大滞后0.6 ps时，本例MSD约5.0160 Ų；CSV保留全部300个滞后时间。

扩散系数需要长时间线性区斜率、维度因子和可靠统计误差。本篇只演示短轨迹数据处理，未据此新增一个“已收敛扩散系数”。延长轨迹、分块比较、检查温度与有限尺寸后再作物理解释。

4 简单的边界核对#

包内还包含一个明确构造的4帧示意：10 Å盒中分数坐标0.95→0.99→0.03→0.07。正确展开后相对首帧位移为0、0.4、0.8、1.2 Å，对应MSD为0、0.16、0.64、1.44 Ų。这是程序边界处理示意，不是额外MD模拟。

当前模块未安装可选的 pymatgen-analysis-diffusion 扩展，本篇使用已安装的核心Trajectory接口。需要更广泛的扩散、配位和RDF分析，也可参考MDAnalysis专题。