ABACUS 色散能量与原子力核对

下载完整输入、原生能量/力/能级及绘图脚本，解压进入附件根目录。先用已有结果生成表格和四张图：

bash analyze.sh

该步骤只读参考数据，用一个CPU核心处理，不启动DFT。完整重算：

sbatch submit.slurm

脚本申请一张P100，顺序运行20个条件，每个1个MPI进程、1个CPU线程、GPU双精度。参考环境为ABACUS3.10.1、GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8、OpenBLAS0.3.30。ABACUS这一构建自带的Libxc是6.2.2，不能套用其他软件扩展的Libxc版本。不要在登录节点执行DFT。本例固定结构、scf_thr=1e-6；截断能和k网格是教学设置，换体系后需自己检查。

1. 除了能量，还能提取什么#

继续采用双层石墨烯的完整输入。本例在3.35 Å附近，分别建立3.34和3.36 Å的几何，每个几何都计算PBE和PBE+D3(BJ)，中点也保留相同对照，共6个实际SCF结果用于这项检查。

所有输入都设cal_force 1，原子位置固定；小位移是为了数值求导，不是执行结构优化。

2. 两层怎样移动#

改变层间距h时，保持两层中心位置不变，下层移动−δh/2，上层移动+δh/2。每层有两个原子。比较的力应为整个上层两个原子的z向力之和，不能只取一个C原子的力。

在总平移力为零时，F_upper,z = −dE/dh。这里先从D3(BJ)总力减去同几何PBE总力，再与原生色散项的中心差分比较：

F_dispersion = F_upper(PBE+D3BJ) − F_upper(PBE)
F_difference = −[E_dispersion(h+δh) − E_dispersion(h−δh)] / (2δh)

δh约0.01 Å；脚本同时按ABACUS原生长度常数换算输入坐标，避免把Bohr与Å混用。

3. 提取原子力#

完整计算完成后，在附件根目录执行：

sed -n '/TOTAL-FORCE/,+6p' work/cases/graphene-d3_bj-3.35/OUT.tutorial/running_scf.log
sed -n '/TOTAL-FORCE/,+6p' work/cases/graphene-none-3.35/OUT.tutorial/running_scf.log

原生表单位是eV/Å，本例C3/C4属于上层。若重新排列原子顺序，必须重新识别层归属，不能机械相加第3和第4行。

python - <<'PY'
import json
with open('reference/analysis/RESULTS.json') as f:
    a=json.load(f)
for key,value in a['D3BJ_force_check'].items():
    print(key,value)
PY

负号表示色散部分把上层拉向下层。这个一致性检查针对色散能与色散附加力，不代表整体DFT力已精确到10⁻⁷ eV/Å，也不代表结构已经平衡；PBE部分仍可能有明显层间排斥。

4. 自己复用流程#

保留相同单胞、k点、泛函和数值设置，仅对目标距离做小位移；每个几何都使用自己的PBE/色散成对结果。位移太大时中心差分截断误差增大，太小时SCF/输出数值噪声可能占主导。需要高精度时再做不同δh和更严SCF的对照，不能把一次吻合推广到所有参数。