GPUMDkit 空位迁移NEB实操

NEB用多个中间结构连接起点和终点，寻找势能面上的迁移路径。本例计算Si与PbTe空位附近的一次跳跃，使用NEP势、215原子、3个中间像。外围固定的小模型用于说明流程，结果不等同于尺寸、路径和模型均已收敛的DFT迁移势垒。

1. 下载和认识输入#

下载输入、最终结构、参考数据与绘图程序，解压进入资料包目录。

module load gpumdkit/1.5.8
python run.py si-calorine

这会在work/si-calorine计算，inputs/si-calorine保持原样。计算应在申请的计算资源中执行；该入口调用CPU上的Calorine/ASE，文件名中的NEP不表示一定占用GPU。

Si由3×3×3常规晶胞删去一个原子得到215原子；邻近原子跳入空位，保留相同原子编号。距迁移区域较远的196个原子固定，19个自由。PbTe同样215原子，Te子晶格的一次邻位跳跃，180固定、35自由。两者使用不同材料的模型，不能交换nep.txt。

2. 自己建体系时先优化端点#

起终态需要相同晶胞、相同原子数及顺序，仅把目标迁移原子移到终点位置。固定区域也应保持一致。不要分别按坐标排序两个结构，否则可能把不同原子错误地连接起来。

附件endpoints/si、endpoints/pbte保留原始端点和实际优化程序；在各自的新副本中运行：

python optimize_endpoints.py

关键部分为：

atoms.calc = CPUNEP("nep.txt")
atoms.set_constraint(FixAtoms(indices=fixed))
opt = BFGS(atoms, logfile="endpoint.log")
converged = opt.run(fmax=0.005, steps=300)
if not converged:
    raise RuntimeError("端点还未收敛")

这里0.005 eV/Å是本例端点标准；固定原子上的非零力不参与自由原子收敛判断。实际Si端点14步，最大自由原子力约0.00404 eV/Å；PbTe端点24步，约0.00377 eV/Å。初终态是对称等价的局域空位配置，能量差接近数值零。

3. 运行NEB#

在含四项输入的新目录执行：

gpumdkit.sh -calc neb init.xyz final.xyz 3 nep.txt

交互时先输入index，再粘贴fixed_indices.txt里的一整行。或直接使用附件答案：

gpumdkit.sh -calc neb init.xyz final.xyz 3 nep.txt < answers.in

3表示3个中间像，加起终态共5像。本入口用IDPP生成初始路径、FIRE优化、爬山NEB和0.01 eV/Å的NEB投影力阈值，弹簧常数0.5。端点优化标准和NEB投影力是两件事。

运行PbTe的完整案例：

python run.py pbte-calorine

4. 检查输出并读势垒#

tail neb.log
ls images

images/image_0.xyz到image_4.xyz是最终5个像；neb.xyz包含优化过程的全部历史像，不能把全部帧当作一条最终路径。nep_barrier.png是原生路径图。应同时检查收敛力、路径结构、固定原子和能量曲线，不能只看程序是否退出。

势垒是路径最高能量减初态能量，不是最高能量的绝对值。可从最终像独立提取：

from ase.io import read
from calorine.calculators import CPUNEP
import numpy as np
calc = CPUNEP("nep.txt")
energies = []
for i in range(5):
    a = read(f"images/image_{i}.xyz")
    a.calc = calc
    energies.append(a.get_potential_energy())
relative = np.array(energies) - energies[0]
np.savetxt("path_energy.dat", np.c_[range(5), relative],
           header="image relative_energy_eV")
print("Barrier / eV:", relative.max())

本例5像最高点是中间像。更复杂路径应增加像数、检查中间构型，并可能从多个初始路径寻找不同机制。

5. 可选NepTrainKit后端#

这个补充入口需要独立的NepTrainKit环境，不是共享GPUMDkit模块自动包含的依赖。本例版本3.1.5，CPU后端实际计算的PbTe路径与Calorine一致。可以在自己的目录创建环境：

module load gpumdkit/1.5.8
python -m venv --system-site-packages ~/venvs/neptrainkit-3.1.5
~/venvs/neptrainkit-3.1.5/bin/python -m pip install 'NepTrainKit==3.1.5'
export NEPTRAINKIT_PYTHON=~/venvs/neptrainkit-3.1.5/bin/python
python run.py pbte-neptrainkit

安装会下载Qt等依赖，并可能编译扩展；本例在G++13.3工具链下完成。无图形界面执行时设置QT_QPA_PLATFORM=offscreen。原生命令为：

"$NEPTRAINKIT_PYTHON" \
  "$GPUMDkit_path/Scripts/calculators/neb_calculation_neptrain.py" \
  init.xyz final.xyz 3 nep.txt < answers.in

本次适配器的CPU实现支持的angular n_max上限为8，不接受本例Si模型的10；PbTe模型为6，可以运行。不要只改模型首部数字绕过检查，那会改变模型含义。

6. 用于自己材料前要补什么#

先检查NEP是否包含缺陷和过渡区域的训练信息，再考察超胞、自由区域半径、像数、端点和NEB收敛。最终用于材料结论时，建议用相应DFT设置复核端点和鞍点附近的能量、力。NEB得到的是本模型的势能势垒；温度、振动自由能和跃迁速率还需要额外计算。

参考：ASE NEB、NepTrainKit、GPUMDkit。