ABACUS 轨迹、RDF与MSD实操

下载13个条件的完整输入、原生轨迹、续算状态与绘图脚本，解压进入附件根目录。先用已有数据生成图表：

bash analyze.sh

结果在reference/analysis、图片在reference/figures、ASE轨迹在reference/export。从头执行整套流程：

sbatch submit.slurm

脚本在新的work中依次运行13个条件，已有work时请新解压，避免覆盖已完成结果。想先跑一个条件，可在计算资源环境中加载下述模块、运行setup.py后用run_case.py指定目录名；普通用户不要在登录节点运行DFT。提交模板申请1张P100、1个MPI进程、1个线程，套餐用户按实际权限用yskit生成脚本。

本例ABACUS3.10.1模块为abacus/3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8环境，实际为P100双精度；后处理为ase/3.29.0。模型采用32原子fcc铝、3价电子LDA赝势Al.pz-vbc.UPF及PZ泛函，40Ry、Γ点、64条带，共96个价电子。不能用同名11价电子赝势替换后仍保留这些带数与精度设置。

这批轨迹只有40fs，用来指导输入、读取和续算，不用于扩散系数、熔点、平衡密度或收敛热力学性质。

1 从原生MD_dump到ASE轨迹#

附件analyze.py读取每帧晶胞、坐标、力和速度，生成reference/export/*.extxyz。原生LATTICE_VECTORS要乘同一帧的LATTICE_CONSTANT（Å）；位置已经是Å，速度已经是Å/fs，不能再次乘晶格常数。力是eV/Å。

先用ASE实际读取并导出首帧：

module load ase/3.29.0
python - <<'PY'
from ase.io import read,write
frames=read('reference/export/al-nve.extxyz',index=':')
print('frames',len(frames),'atoms',len(frames[0]))
print('volume/A3',frames[0].get_volume())
print('first velocity/A per fs',frames[0].arrays['velocity_A_per_fs'][0])
write('al-first.cif',frames[0])
PY

这里的velocity_A_per_fs是显式单位的数组，不能直接当作ASE内部时间单位的速度。需要ASE动量/动力学操作时应先正确换算。

2 用ASE邻居表提取RDF#

截断半径取3.6Å，严格小于本批每一帧最短晶胞面间距的一半。邻居表每个原子对按两个方向计数，排除自身，按球壳体积及N(N−1)/V归一化。有限N约定已明确，不能与采用N²归一化的工具不加说明地直接比数值。

module load ase/3.29.0
python - <<'PY'
import numpy as np
from ase.io import read
from ase.neighborlist import neighbor_list
a=read('reference/export/al-nve.extxyz',index=-1)
i,j,d=neighbor_list('ijd',a,3.6,self_interaction=False)
edges=np.linspace(0,3.6,181)
shell=4*np.pi/3*(edges[1:]**3-edges[:-1]**3)
g=np.histogram(d,bins=edges)[0]/(len(a)*(len(a)-1)*shell/a.get_volume())
print('coordination within3.3A',np.count_nonzero(d<3.3)/len(a))
print('RDF max',g.max())
PY

本例第一近邻配位数约12，符合短时间fcc结构。图中RDF使用全部实际帧的同一约定平均，不能从40fs晶体RDF推断液体平衡结构。

3 MSD必须处理真实跨边界#

原子起初就在周期边界附近，第一小步就可能从0附近回到8Å附近。直接减去包装坐标会得到巨大的假位移。对固定晶胞，用ASE的最小镜像把相邻帧位移连续累加：

module load ase/3.29.0
python - <<'PY'
import numpy as np
from ase.io import read
from ase.geometry import find_mic
frames=read('reference/export/al-nve.extxyz',index=':')
u=frames[0].positions.copy()
start=u.copy()
for previous,current in zip(frames[:-1],frames[1:]):
    step,_=find_mic(current.positions-previous.positions,current.cell,pbc=True)
    u+=step
msd=np.mean(np.sum((u-start)**2,axis=1))
naive=np.mean(np.sum((frames[-1].positions-start)**2,axis=1))
print('unwrapped MSD/A2',msd)
print('naive wrapped MSD/A2',naive)
PY

本批MSD保留原生的整体质心运动，随机浴产生的整体漂移没有被悄悄减掉。研究相对质心的自扩散时，应明确去漂移的定义并做长时统计；不能把整体漂移当成自扩散。

本例相邻步位移远小于半个晶胞，展开没有跨越多格的不确定性。存帧太稀疏时不能靠最小镜像找回丢失的轨迹。NPT还需要处理晶胞仿射变形；附件对其明确使用参考晶胞映射，不能照抄上面固定晶胞代码来解释扩散。

4 最后得到哪些文件#

conditions.csv汇总各条件，*-thermo.csv是逐帧温压/能量/体积，rdf.csv给出RDF，trajectories.npz保存实际数组，restart-comparison.csv给出续算误差，export/*.extxyz供ASE及可视化软件使用。图件为高清PNG。不要用本批短时MSD强行拟合扩散系数，也不要把未平衡的瞬时温压平均当成材料常数。

工具操作参阅ASE读写、ASE邻居表和ASE周期几何。