GPUMDkit离子电导实操

本篇从LLZO固态电解质的真实GPUMD轨迹出发，先看Li是否扩散，再算自扩散系数和Nernst–Einstein电导，最后比较不同温度与随机初速度。它适合学习操作流程；192原子和100ps生产段尚不足以确定收敛的材料电导率。

下载8组实算输入、结果、轨迹和脚本。模型为作者公开的NEP25-ZBL，初始坐标来自该模型训练结构中的一个192原子周期晶胞。因此这不是独立模型精度测试，也不是上游12288原子长轨迹的复现。8组数据均由鸢算科技P100实际运行得到。

1. 先准备已有轨迹的后处理#

在附件根目录执行：

module load gpumdkit/1.5.8
python new-case.py ionic
cd work/ionic/seed91/1000K
head -2 model.xyz
head -2 msd.out

目录按seed91/950K、1000K、1100K、1200K排列；seed92保存另一组初速度。每条轨迹含56Li、24La、16Zr、96O。model.xyz里的group为Li=0、La=1、Zr=2、O=3。第一组必须是要计算电导的Li，因为当前电导脚本直接读取MSD的第1—3个数据列，不根据元素参数重新挑列。

2. 看清计算输入与时间#

1000K、seed91的run.in如下，其他目录只调整温度与seed：

potential nep.txt
velocity 1000 seed 91
time_step 1
ensemble nvt_nhc 1000 1000 100
run 10000
ensemble npt_mttk temp 1000 1000 aniso 0 0
run 20000
ensemble npt_mttk temp 1000 1000 aniso 0 0
compute_msd 10 1000 all_groups 0 save_every 20000
dump_thermo 1000
dump_xyz 1000 trajectory.xyz
run 100000

时间步1fs：先10ps NVT，再20ps NPT平衡，最后100ps NPT生产。MSD每10步采样，1000个相关点覆盖约10ps的相关滞后时间；它与100ps生产长度不同。save_every 20000每20ps保存一份累计统计，最后得到20、40、60、80、100ps五份MSD。

想重新计算时，回到附件根目录，检查Slurm脚本中的账号适用分区，然后提交：

cd ../../../..
cat compute.sbatch
sbatch compute.sbatch
squeue -u "$USER"

脚本采用gpumd/5.9.1-p100，8项数组、每项1张P100，结果写入新的compute-作业号/温度K-s种子目录。它只复制输入和NEP势，不把旧输出当新结果。普通用户应通过Slurm计算，不能直接在登录节点执行MD。若只想学习后处理，不必重新提交。

3. 画各元素MSD，再计算Li电导#

回到刚才的1000K工作目录：

cd work/ionic/seed91/1000K
gpumdkit.sh -plt msd_all msd.out Li La Zr O save
gpumdkit.sh -calc ionic-cond Li 1

Li 1表示载流离子Li和电荷数+1。图中Li的MSD随时间增长，其余三类在本次相关时间内接近平坦。总MSD三维扩散系数为斜率/6；分方向分别除2。当前工具拟合数据40%—80%段，本例约4—8ps。

Nernst–Einstein近似为σ=nq²D/(kBT)。它使用单粒子自扩散，没有计入不同离子运动之间的交叉关联；与完整集体电导之间可能存在明显差别。ionic-cond终端单位是mS/cm，后面的sigma绘图单位是S/cm，相差1000倍。

ionic-cond采用thermo后半段平均温度/体积；温度目录绘图脚本采用目录名温度和整个生产段平均体积。所以两者电导末位可能略有不同。不要为了对齐数字修改原轨迹。

4. 多温度绘图#

cd ..
gpumdkit.sh -plt sigma save
gpumdkit.sh -plt sigma_xyz save

输出Arrhenius_sigma.png和Arrhenius_sigma_xyz.png。要求每个数字K子目录均有msd.out、thermo.out、model.xyz；本例不使用replicate指令，离子总数就是56。当前sigma系列按Li/Na计数、假设单位价态，不应直接用于Mg²⁺等多价离子。

下面列出脚本采用的温度目录口径，D单位cm²/s，σ单位S/cm：

换到第二组目录：

cd ../seed92
gpumdkit.sh -plt sigma save
gpumdkit.sh -plt sigma_xyz save

电导Arrhenius拟合使用log10(σT)对1000/T，不能把扩散的log D拟合公式不加区别地搬过来。seed91得到0.480eV、R²=0.831；seed92为0.305eV、R²=0.904。两组差异以及个别温度点的非单调性说明短轨迹仍有统计噪声。正式研究应增加轨迹长度、独立初态和超胞尺寸，并检查势的适用范围；不要把这两个活化能平均后当作高精度结论。

工具会自动打印300K外推值，但这里只计算950—1200K，短采样也未收敛，本教程不把该外推解释为可靠室温电导。

5. 分段图只是绘图操作，不是相变判定#

gpumdkit.sh -plt sigma_PT 1050 save

输出Arrhenius_sigma_PT.png。1050K人为将本例分成950/1000K和1100/1200K两段，各段只有两个温度点；它不能证明发生相变，也不能给出稳健的两段活化能。真正存在相变证据时，再根据结构、热力学和更密集温度采样选择分段位置。

6. 最后保留什么#

保留model.xyz、nep.txt、run.in、gpumd.log、thermo.out、msd.out、msd_step*.out及trajectory.xyz，连同分析命令和拟合窗口一起归档。图片可以重画，原始轨迹和模型出处不可丢。

回到附件根目录，python plot_results.py可重画本篇结构/数据图。图和终端数字所用的数据保存在附件RESULTS.json，完整原始输出在llzo目录。

模型及原始公开教程：GPUMD官方LLZO示例、作者NEP25-ZBL资料。进一步阅读：扩散拟合、GPUMDkit全部绘图入口。