qNEP静态计算与电荷输出

qNEP是在NEP中引入电荷相关静电项的模型。本篇直接使用GPUMD官方已训练的PbTe参考模型，在鸢算科技P10016G和V10016G上实际执行静态计算。软件为GPUMD5.9.1，GCC13.3/CUDA12.9；分析用ASE3.29.0/Python3.11及NumPy/Matplotlib。

本次没有重新训练模型，也没有新做DFT或MD。选取的8帧来自该参考模型的训练数据，用于说明运行与数值检查，不作为独立模型精度测试。

1. 下载并查看参考结果#

下载模型、输入、实际原生结果和脚本，保存为qnep.zip：

unzip qnep.zip -d qnep
cd qnep
head -n 6 inputs/nep.txt
bash analysis.sh
ls analysis figures

bash analysis.sh只分析包内真实参考结果，不占GPU；输出在analysis/与figures/。模型文件开头：

nep4_charge2 2 Te Pb 
cutoff 8 4 73 8
n_max 6 6
basis_size 6 6
l_max 4 2 0
ANN 30 0

charge2表示本模型采用电荷模式2，即倒空间静电贡献路线。它不同于极化率模型的model_type 2；本例是能量型模型加charge_mode 2。关于两种电荷模式见官方charge_mode说明。

2. 建一个最小静态目录#

mkdir my-static
cp inputs/nep.txt my-static/
cp reference-data/p100-f24/baseline/model.xyz my-static/
cp reference-data/p100-f24/baseline/run.in my-static/
cp static.slurm my-static/
cd my-static
cat run.in
sbatch static.slurm

static.slurm是单卡P100提交模板。用V10016G时，同步选择V100分区、GRES类型和gpumd/5.9.1-v100模块；资源申请方法见环境指南。单个任务启动一个GPUMD进程，不使用mpirun复制同一作业。

本次实际交给GPUMD的完整run.in：

potential nep.txt
kspace ewald
velocity 1 seed 13579
ensemble nve
time_step 0
dump_xyz 1 dump.xyz force charge precision double
dump_thermo 1
run 1

time_step 0保持坐标，run 1触发一次求值，没有MD时间演化。velocity 1用于初始化这次静态调用，不表示在1K进行了材料动力学采样。force charge precision double将力和模型电荷写入扩展XYZ。

本例有250个原子，model.xyz开头如下，完整坐标在附件中：

250
pbc="T T T" Lattice="0.0 16.42598 16.42598 16.42598 0.0 16.42598 16.42598 16.42598 0.0" Properties=species:S:1:pos:R:3
Te 3.2716189999999998 3.0938159999999999 3.3547220000000002
Te 3.7261190000000002 6.5249110000000003 6.8534889999999997
Te 3.062335 9.8543559999999992 9.6550309999999993
...

Properties告诉程序各列含义。这里输入只有元素和坐标；DFT能量/力标签保留在原始参考帧中，不冒充位移后结构的标签。

3. 提取能量、力与电荷#

任务完成后，在my-static中：

head thermo.out
head -n 5 dump.xyz
module purge
module load ase/3.29.0
python - <<'PYTHON'
from ase.io import read
a = read('dump.xyz')
print('Atoms:', len(a))
print('Forces / eV Angstrom:', a.get_forces())
print('Model charges / e:', a.get_charges())
print('Total model charge / e:', a.get_charges().sum())
PYTHON

thermo.out第三列是整个构型的势能eV；力是eV/Å，电荷输出为模型电荷。本例P100第24帧的Ewald势能为 −930.79421934eV，最大模型电荷绝对值约 0.577067e，总电荷约 2.77×10⁻⁶e，数值上接近中性。

这些电荷是势函数内部变量，没有在本篇用目标原子电荷或BEC做独立物性验证，不能直接当作Bader电荷或据此指定氧化态。这份参考模型输出中Te偏正、Pb偏负；仅用能量/力拟合的内部电荷仍有符号等模型约定，应结合额外参考解释。当前本例运行路线会中性化原子电荷；任意净带电体系不能直接照搬。

4. 为什么nep和GPUMD的电荷文件不能直接相减#

GPUMD的dump.xyz写中性化后的电荷。当前NEP预测的charge_train.out写未平移的网络输出；静电计算内部另使用满足目标总电荷的数组。因此两个文件的电荷零点约定不同。

本例第24帧NEP原始输出之和约0.243112e。为比较目标总电荷0的同一数组定义，可在解压根目录新建派生数组：

python - <<'PYTHON'
import numpy as np
q_raw = np.loadtxt('reference-data/p100-f24/prediction/charge_train.out')
q_neutral = q_raw - q_raw.mean()
print(q_raw.sum(), q_neutral.sum())
PYTHON

不要覆盖原始文件。按同一约定比较后，16组结果的最大电荷分量差约 1.64×10⁻⁶e。这个对照说明输出约定与实现一致，不是对真实原子电荷准确度的证明。

5. 重新跑配套对照#

在解压根目录：

python3 prepare.py
python3 make_jobs.py
cd cases/p100-f24
sbatch job.slurm

每个case运行一个参考帧，包含Ewald静态计算、NEP预测、三个差分步长及三个PPPM网格。其他case在各自目录提交；全部16组结束后回到根目录bash analysis.sh --new，新结果写analysis-new/figures-new/。只跑单例时，直接看该case的原生输出即可，完整汇总需要16组完成。

实际项目不必每次都重复整套教学矩阵；先用自己的模型、坐标与最小run.in跑通，再按需要检查精度。检查点新增零电荷头的操作见GPUMDkit电荷模型；本篇使用的是已训练参考模型，不能只改文件头来代替训练。求解器和力检查继续下一篇。