ASE调用ABACUS：独立接口与结果单位

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

1. 准备环境和运行#

本例调用ABACUS3.10.1的CPU构建。ABACUS calculator来自ABACUS维护的独立ASE-ABACUS分支（版本标记3.26.0b1），共享标准ASE3.29.0本身没有这个calculator。请把分支放在自己的目录，并固定已测试提交：

git clone https://gitlab.com/1041176461/ase-abacus.git
cd ase-abacus
git checkout 65bdad3493fea002180cdc73664d688a2087a99b
export ASE_ABACUS_SOURCE="$PWD"
cd ..

回到资料包目录：

module load ase/3.29.0 abacus/3.10.1-8173m
bash run.sh abacus

run.sh只在这个进程通过PYTHONPATH选择分支，不覆盖共享ASE。独立分支的安装说明或API可能更新，复现时以固定提交为准。

Si原胞2原子，a=5.43 Å，第2原子沿x偏移0.06 Å，并取±0.005 Å。参数为PBE、平面波40 Ry、3×3×3 Γ中心网格、SCF阈值1E-8、CG求解、Gaussian展宽0.002 Ry，4个MPI进程。使用公开Si_ONCV_PBE-1.0.upf；这不是VASP的PAW赝势。

该分支的AbacusProfile参数名是command，用AbacusProfile(command="mpirun -np 4 abacus", ...)。本例PW计算无需数值原子轨道，但接口仍显式设置pseudo_dir和basis_dir，防止从缺省环境变量读取失败。

2. 本例完整计算程序#

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

from pathlib import Path
import json,os,numpy as np
from ase.io import write
from ase.build import bulk
from ase import Atoms
B=Path(__file__).resolve().parent
rows=[]
from ase.calculators.abacus import Abacus,AbacusProfile
profile=AbacusProfile(command='mpirun -np 4 abacus',pseudo_dir=str(B/'PP_ORB'),basis_dir=str(B/'PP_ORB'))
for label,delta in [('minus',-.005),('center',0),('plus',.005)]:
 d=B/label;d.mkdir();a=bulk('Si','diamond',a=5.43);a.positions[1,0]+=.06+delta
 a.calc=Abacus(profile=profile,directory=d,ntype=1,ecutwfc=40,scf_thr=1e-8,scf_nmax=150,smearing_method='gaussian',smearing_sigma=.002,basis_type='pw',ks_solver='cg',calculation='scf',pp={'Si':'Si_ONCV_PBE-1.0.upf'},kpts=(3,3,3),kmode='Gamma',xc='PBE',device='cpu',cal_force=1,cal_stress=1,symmetry=0)
 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.892815395 eV/Å
能量差分：-0.892862000 eV/Å
绝对差异：4.66e-05 eV/Å

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

ABACUS日志中的TOTAL-STRESS (KBAR)采用压力型符号；ASE张力为正且单位eV/ų。本例按σ_ASE = −σ_ABACUS × 0.0006241509独立换算，再按Voigt顺序xx, yy, zz, yz, xz, xy对照。TOTAL-FORCE已经是eV/Å，不应重复乘单位常数。

Si示例与VASP页的原子数和位移相同，但赝势、截断、展宽及计算实现不同；绝对总能量不应直接相减。要比较方法结果，应先统一物理条件并单独收敛。

参考：ASE-ABACUS源码。其他联用：VASP、ABACUS、CP2K、NEP。