ABACUS 分子优化与反应能

本篇用 2H₂ + O₂ → 2H₂O 演示“分别优化反应物与产物，再按计量系数相减”。结果是电子反应能；零点振动能、温度修正和熵没有加入，不能直接称为反应自由能。

下载完整输入、真实结果与绘图脚本，解压后进入包含 run.sh 的目录。以下脚本使用 ABACUS 3.10.1 P100 版（GNU 13.3、HPC-X 2.50 / Open MPI 5.0.10rc2、CUDA 12.8、OpenBLAS 0.3.30），分析采用 ASE 3.29.0。每个条件使用一张 GPU、一个 MPI 进程和一个线程，按顺序计算，适合先学习小体系流程。

1 提交这一组计算#

sbatch submit.slurm reaction

这是核时 P100 模板。套餐用户先用 yskit 生成自己的分区脚本，保留脚本里的计算命令。run.sh 在分配的计算资源中调用；不要在登录节点直接执行科学计算。脚本中按组运行的命令很短：

bash run.sh reaction

脚本自动准备 work/cases，完成计算后生成 work/analysis、work/converted 和 work/figures。重新计算请新解压一份，避免覆盖已有结果；reference 是本篇的实际输入输出。只想重画参考结果时，在资料包目录执行：

module load ase/3.29.0
python analyze.py reference
python convert_structures.py reference
python plot_results.py reference

电子收敛参数 scf_thr 1e-8 指 ABACUS 的密度残差；几何优化的 force_thr_ev 0.02 单位为 eV/Å。优化完成后，脚本会用最终结构再做一次静态 SCF，重新核对能量与自由原子的力。

2 分别优化三个分子#

三个分子都采用PBE、同一套H/O赝势、60 Ry平面波截断、12 Å立方盒和Γ点。H₂与H₂O使用闭壳层；O₂采用共线自旋 nspin 2、nupdown 2，对应两个自旋通道电子数之差为2。不要把O₂错误地强制为闭壳层。该盒长是本例的有限周期近似，正式气相参考还需盒长收敛。

h2-relax：完整 INPUT#

INPUT_PARAMETERS
suffix tutorial
calculation relax
basis_type pw
ecutwfc 60
scf_thr 1e-8
scf_nmax 150
ks_solver cg
device gpu
precision double
cal_force 1
cal_stress 1
pseudo_dir ../PP
symmetry 0
smearing_method fixed
mixing_beta 0.3
nbands 4
relax_nmax 60
force_thr_ev 0.02
relax_method cg
out_stru 1

ATOMIC_SPECIES
H 1.00794 H_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.8897261246257702

LATTICE_VECTORS
12.000000000000 0.000000000000 0.000000000000
0.000000000000 12.000000000000 0.000000000000
0.000000000000 0.000000000000 12.000000000000

ATOMIC_POSITIONS
Direct
H
0.0
2
0.500000000000 0.500000000000 0.468750000000 1 1 1
0.500000000000 0.500000000000 0.531250000000 1 1 1

K_POINTS
0
Gamma
1 1 1 0 0 0

o2-triplet-relax：完整 INPUT#

INPUT_PARAMETERS
suffix tutorial
calculation relax
basis_type pw
ecutwfc 60
scf_thr 1e-8
scf_nmax 150
ks_solver cg
device gpu
precision double
cal_force 1
cal_stress 1
pseudo_dir ../PP
symmetry 0
smearing_method fixed
mixing_beta 0.3
nspin 2
nupdown 2
nbands 12
relax_nmax 60
force_thr_ev 0.02
relax_method cg
out_stru 1

ATOMIC_SPECIES
O 15.9994 O_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.8897261246257702

LATTICE_VECTORS
12.000000000000 0.000000000000 0.000000000000
0.000000000000 12.000000000000 0.000000000000
0.000000000000 0.000000000000 12.000000000000

ATOMIC_POSITIONS
Direct
O
1.0
2
0.500000000000 0.500000000000 0.448750000000 1 1 1
0.500000000000 0.500000000000 0.551250000000 1 1 1

K_POINTS
0
Gamma
1 1 1 0 0 0

water-relax：完整 INPUT#

INPUT_PARAMETERS
suffix tutorial
calculation relax
basis_type pw
ecutwfc 60
scf_thr 1e-8
scf_nmax 150
ks_solver cg
device gpu
precision double
cal_force 1
cal_stress 1
pseudo_dir ../PP
symmetry 0
smearing_method fixed
mixing_beta 0.3
nbands 12
relax_nmax 60
force_thr_ev 0.02
relax_method cg
out_stru 1

ATOMIC_SPECIES
O 15.9994 O_ONCV_PBE-1.0.upf
H 1.00794 H_ONCV_PBE-1.0.upf

LATTICE_CONSTANT
1.8897261246257702

LATTICE_VECTORS
12.000000000000 0.000000000000 0.000000000000
0.000000000000 12.000000000000 0.000000000000
0.000000000000 0.000000000000 12.000000000000

ATOMIC_POSITIONS
Direct
O
0.0
1
0.500000000000 0.500000000000 0.500000000000 1 1 1
H
0.0
2
0.500000000000 0.563603250000 0.450307583333 1 1 1
0.500000000000 0.436396750000 0.450307583333 1 1 1

K_POINTS
0
Gamma
1 1 1 0 0 0

3 从优化结构做静态能量#

程序完成每个 relax 后，读取最终 STRU_ION_D 或 STRU_NOW，在该例的 final-static 中以 calculation scf 重新计算。可以这样查看：

grep -E -i 'relaxation is converged|final etot' work/cases/*/OUT.tutorial/running_relax.log
grep -E 'final etot|total magnetism' work/cases/*/final-static/OUT.tutorial/running_scf.log
cat work/analysis/RESULTS.json

最终静态最大力分量：H₂约0.007594、O₂约0.013289、H₂O约0.008371 eV/Å，均满足本例0.02的几何目标。O₂实际总自旋人口差为2。自旋初猜和电子态是能量参考的一部分，不能只检查程序是否退出。

4 按配平系数计算差值#

本例采用：ΔE = 2E(H2O) − 2E(H2) − E(O2)。analyze.py 从三份最终静态日志提取能量并计算；生成的数据保存在 work/analysis。

负号表示这一组电子能参考下产物更低。它不包含气体标准态、零点能、热焓和熵；若需要自由能或自由能台阶图，应在一致标准态下补上这些项。

5 套用其他反应#

先配平所有元素和电荷，确定各分子的自旋，再按同一泛函、赝势与数值标准分别优化。比较的是按计量系数组合后的能量；不要比较原子数不同体系的裸总能，也不要把本例PBE水能量作为下篇PBE-D3吸附的水参考。

继续看吸附能、分子键长扫描及局域基组Counterpoise。