ABACUS：D4与D4S能量、受力和应力

在ABACUS里加入D4，不只是给总能量加一项。本例同时读取原生色散能、色散力和色散应力，并比较D4与D4S两种模型。

模型是自建的6原子水双体，放在非正交的周期晶胞中；结构固定且未优化，原子受力较大。这里用于学习接口和单位，不把总能差解释成收敛的气相结合能，也不把该几何作为水双体平衡结构。

下载输入、脚本与参考数据。

unzip abacus-d4-wave690-*.zip
cd abacus-d4
bash results.sh
bash plot.sh

以上只读取已有参考数据；重新计算请使用新解压的目录，结果写入work/和analysis/。

1. 选择环境并准备输入#

先按个人环境教程编译一次，再执行：

bash use-toolchain.sh ../abacus-d4-env/toolchain
bash prepare.sh
cat inputs/d4/STRU

完整STRU：

ATOMIC_SPECIES
O 15.999 O_ONCV_PBE-1.0.upf
H 1.008 H_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.889726877774355

LATTICE_VECTORS
8.0 0 0
1.0 8.0 0
0.6 0.4 8.0

ATOMIC_POSITIONS
Cartesian
O
0
2
2.0 2.0 2.0 1 1 1
2.0 2.0 5.0 1 1 1
H
0
4
2.758602 2.0 2.504284 1 1 1
1.241398 2.0 2.504284 1 1 1
2.758602 2.0 5.504284 1 1 1
1.241398 2.0 5.504284 1 1 1

LATTICE_CONSTANT为每Å的Bohr数；晶格矢量和笛卡尔坐标按该长度因子换算。三个晶格矢量不互相垂直，因此可以同时观察应力的对角与非对角分量。

KPT：

K_POINTS
0
Gamma
1 1 1 0 0 0

本例为Γ点教学采样；实际材料仍要收敛电子k网格、截断能和晶胞尺寸。

2. D4输入怎么写#

INPUT_PARAMETERS
suffix d4
calculation scf
basis_type pw
ntype 2
nspin 1
nelec 16
nbands 12
ecutwfc 40
scf_thr 1e-6
scf_nmax 160
dft_functional pbe
ks_solver cg
symmetry -1
smearing_method fixed
mixing_type broyden
mixing_beta 0.3
cal_force 1
cal_stress 1
out_chg 0
vdw_method d4
vdw_d4_xc pbe
vdw_d4_model d4
vdw_cutoff_radius 95
vdw_radius_unit Bohr
vdw_cn_thr 40
vdw_cn_thr_unit Bohr
vdw_cutoff_width2 0
vdw_cutoff_width3 0

关键项：

• vdw_method d4调用原生D4接口。
• vdw_d4_xc pbe选择PBE对应的阻尼参数；要与实际泛函匹配。
• vdw_d4_model d4为标准D4，改成d4s选择D4S，其他物理设置保持相同。
• vdw_cutoff_radius 95、vdw_radius_unit Bohr设置两体距离截断；本例三体截断40 Bohr、配位数截断40 Bohr。
• vdw_cutoff_width2/3 0不使用平滑宽度。独立比较时必须使用同一组截断和宽度，不能拿另一个程序的默认值直接比。
• cal_force 1、cal_stress 1要求输出力和应力。电子收敛阈值为1e-6。

3. 提交五种条件#

sbatch job.slurm
squeue -u "$USER"

脚本在独立目录计算：

CPU脚本申请56核独占整节点，本例内部一个MPI进程、一个线程，仍按整节点计费。结束后会编译一个轻量C API读数程序，独立调用同一D4库，并提取数据和图片；不会再次运行这五个DFT条件。

4. 去哪里读能量、力和应力#

grep -E 'SCF IS CONVERGED|E_vdwD4|FINAL_ETOT' work/d4/OUT.d4/running_scf.log
grep -A10 'TOTAL-FORCE' work/d4/OUT.d4/running_scf.log
grep -A7 'TOTAL-STRESS' work/d4/OUT.d4/running_scf.log
cat analysis/D4-energy.csv
cat analysis/D4-force.csv
cat analysis/D4-stress.csv

E_vdwD4一行同时列Ry和eV；力为eV/Å，应力为kbar。D4库原生梯度是dE/dR，所以力需要取负号，不能把梯度当成力直接相加。

本例纯PBE总能−926.0198797247 eV，加入D4后总能−926.0497319027 eV，D4修正−0.0298521780 eV。比较的是同一个固定构型，不能把这一项当成双体相对于两个单分子的结合能。

5. 如何判断读到的是色散修正#

同一几何下：

ΔE = E(PBE+D4) − E(PBE)
ΔF = F(PBE+D4) − F(PBE)
Δσ = σ(PBE+D4) − σ(PBE)

这里的色散属于本例的加和修正；独立C API工具使用完全相同的元素、周期晶胞、电荷、泛函、截断和D4模型。它用于核对ABACUS接口的单位与符号，不是另一种电子结构理论。

实际最大差：力约9×10⁻¹¹ eV/Å，应力约7.6×10⁻⁸ kbar。正负位移的D4能量差分给出第一个O的x方向色散力−0.00021980 eV/Å，原生力差值−0.0002198402 eV/Å，两者相差约4.0×10⁻⁸ eV/Å。独立D4有限应变也与相应sigma分量一致。

6. 换成自己的体系#

先对目标体系选择合适的泛函、赝势、电子网格和结构优化条件，再加入匹配的D4参数。带电体系还需核对总电子数与D4电荷定义；本篇只演示中性体系。周期模型包含周期镜像色散，做分子、吸附或层间结合能时，所有参照必须使用一致的设置，并分别完成尺寸与数值收敛。