ASE调用NEP：势函数与力的核对

ASE用统一的Atoms对象保存结构，把能量和力的计算交给calculator。本篇使用真实程序计算三个邻近构型，再把返回的力与能量差分相互核对。下载本组输入、实际结果及绘图程序，解压进入资料包目录。

1. 准备环境和运行#

本例用GPUMDkit1.5.8环境中的Calorine CPUNEP，通过ASE直接计算NEP模型的能量、力和应力。它在CPU中执行，不会因模型名为NEP就自动启动GPUMD GPU程序。

module load gpumdkit/1.5.8
bash run.sh ml-nep

模型为官方Si NEP4示例，包内inputs/nep.txt保留原始参数。结构是a=5.43 Å的Si原胞2×2×2复制，共16原子，第2原子沿x位移0.06 Å，再取±0.005 Å。

最小调用形式：

from ase.build import bulk
from calorine.calculators import CPUNEP
atoms = bulk("Si", "diamond", a=5.43).repeat((2, 2, 2))
atoms.calc = CPUNEP("nep.txt")
energy = atoms.get_potential_energy()
forces = atoms.get_forces()
stress = atoms.get_stress()

运行目录应由自己拥有写权限。实际程序既检查模型支持的元素，也依赖原胞、周期边界和模型格式；不能把Si势用于其他元素，不能仅改文件头“扩展”模型。

2. 本例完整计算程序#

run.sh负责加载环境和复制输入，实际计算逻辑如下。ASE返回能量单位eV、力单位eV/Å，应力通常为eV/ų。

from pathlib import Path
import os,json,numpy as np
from ase.build import bulk
from ase.io import write
from calorine.calculators import CPUNEP
B=Path(__file__).resolve().parent;os.chdir(B);model=str(B/'nep.txt');rows=[]
for label,delta in [('minus',-.005),('center',0),('plus',.005)]:
 d=B/label;d.mkdir();a=bulk('Si','diamond',a=5.43).repeat((2,2,2));a.positions[1,0]+=.06+delta;a.calc=CPUNEP(model);e=a.get_potential_energy();f=a.get_forces();s=a.get_stress();write(d/'result.xyz',a);rows.append(dict(label=label,delta_A=delta,energy_eV=e,forces_eV_A=f.tolist(),stress_eV_A3=s.tolist()))
(B/'RESULTS.json').write_text(json.dumps(rows,indent=2)+'\n')

生成三个目录：minus、center、plus。程序保存RESULTS.json与每个构型的result.xyz；DFT案例还保留引擎原始输入输出。要继续结构优化时，可将已配置的calculator交给ASE优化器，但本例没有把三个单点误称为优化轨迹。

3. 实际结果和独立检查#

用中心差分检查第2原子的x方向力：

F ≈ −[E(+δ) − E(−δ)] / (2δ)，δ = 0.005 Å
程序返回：-0.779049384 eV/Å
能量差分：-0.779043807 eV/Å
绝对差异：5.58e-06 eV/Å

运行python check_results.py reference检查附件参考；复算后执行python check_results.py work检查自己的输出。程序分别核对原始输出和ASE结果；有限差分有截断误差，SCF和积分网格也有误差，不能要求两者逐位相同。相互吻合说明本例接口的主要单位与符号一致，不能证明材料精度已经收敛。

NEP的导数可以内部自洽，却仍可能偏离真实DFT或实验。能量差分核对与模型泛化能力是两种检查。本例只使用已有模型推断，没有重新训练，也没有验证缺陷、表面、高温或其他成分的可用性。NEP例16原子，DFT页Si例2原子，不用于直接比较绝对能量。

参考：Calorine calculators。其他联用：VASP、ABACUS、CP2K、NEP。