ASE 周期轨迹展开与均方位移

本例使用集群现有 ASE 3.29.0、Python 3.11、NumPy 2.3.5 与 Matplotlib 3.11.2，处理已有真实 GPUMD/VASP 输出，不需要新提交 GPU 或 DFT 任务。下载完整结构、原始结果和可运行脚本，在新解压目录执行：

bash run.sh

analysis/ 保存 CSV、结果汇总和可继续使用的数据，figures/ 保存三张高清图。脚本只用一个 CPU 线程；重复运行请另解压一份，以免混入旧输出。

1. 确认轨迹帧间隔与原子身份#

这里沿用Si/Ge热接触的实际轨迹：1664 原子，x 开放、y/z 周期，两端固定组为 0。生产段 20 ps，每 1000 个 1 fs 步存一帧，共 20 帧；首帧在生产 1 ps，末帧在 20 ps，彼此跨度为 19 ps。

head -n 5 data/contact.xyz
head -n 5 data/contact-initial.xyz

本篇时间横轴从第一张存储帧开始计零，因此是 0–19 ps。传给 ASE 的时间步应为存储帧间隔 1000 fs，不能误用积分步长 1 fs。原子顺序必须在所有帧中一致；若轨迹重排了原子，要先按稳定的原子编号恢复对应。

2. 先逐帧展开周期坐标#

原子在周期方向从接近盒长移动到接近 0，真实位移可能很小，原始包裹坐标却会相差整盒。直接对包裹坐标平方会把边界跳跃误算成几十 Ų 的 MSD。

from ase.io import read
from ase.geometry import find_mic
frames = read('data/contact.xyz', ':')
unwrapped = [frames[0].copy()]
for old, new in zip(frames[:-1], frames[1:]):
    step, distance = find_mic(new.positions - old.positions,
                              new.cell, new.pbc)
    f = new.copy()
    f.positions = unwrapped[-1].positions + step
    unwrapped.append(f)

find_mic 只在周期方向取最小镜像，不会把开放 x 方向的真实位移卷回去。本例相邻存储帧的最大原子位移约 0.4643 Å，远小于周期方向半盒长 11 Å。

如果存帧太稀，原子已跨过半盒，单靠包裹坐标就不能唯一还原位移。应提高输出频率，或在 GPUMD 中直接请求 unwrapped_position 等原生未包裹坐标。变胞轨迹还需进一步明确去除晶胞形变的方式，不能不加区分套用这份固定晶胞脚本。

3. 切出真正需要分析的原子#

initial = read('data/contact-initial.xyz')
mobile = initial.arrays['group'] != 0
trajectory = [f[mobile] for f in unwrapped]

剩下 1408 个移动原子，固定的两端不参与移动原子 MSD。分元素统计时，本例有 768 个移动 Ge、640 个移动 Si。

当前 ASE DiffusionCoefficient 的 atom_indices 分母仍使用完整结构中各元素的原子数。因此选择子集时，采用上面这种先切片构建新轨迹的方式，再调用分析器；不要只传子集索引却默认分母也跟着改了。

4. 计算MSD并查看曲线#

包内脚本保存每种元素的 〈|r(t)-r(0)|²〉：

head analysis/contact-msd.csv
ls analysis/contact-mobile-unwrapped.traj
ls figures

如需单独调用 ASE 的分段分析：

from ase import units
from ase.md.analysis import DiffusionCoefficient
analysis = DiffusionCoefficient(trajectory, timestep=1000 * units.fs)
analysis.calculate(number_of_segments=2)
slopes, spread = analysis.get_diffusion_coefficients()

这里 20 帧分为两个 10 帧段。ASE 内部逐方向用 MSD_component/2 对时间做含截距拟合，再平均方向与分段。程序返回的斜率单位是 Ų/ASE时间；乘 units.fs 才是 Ų/fs。spread 是这些方向/分段斜率的散布，不能自动解释为充分统计的材料误差条。

5. 为什么这里不报告扩散系数#

这份短轨迹是带温度梯度的固体热接触，不是平衡液体或长时间离子扩散采样。展开后的 MSD 大致在 0.02–0.03 Ų 附近波动，主要反映有限温度位移；把几段曲线强行拟合出斜率，不会自动变成材料自扩散系数。

要研究扩散，应选择合适的平衡采样，检查较长时间的线性区、多个时间起点、体系尺寸和独立采样；再考虑三维 MSD=6Dt、二维 MSD=4Dt 或各向异性分量的区别。此处重点是坐标展开、子集归一化和时间单位。

依据：ASE DiffusionCoefficient实际接口与单位。原始模型和生成命令见GPUMD热接触。