ABACUS分子键长扫描与能量和力
七个H2固定键长单点,展示结构修改、终端取数、势能曲线和力的方向。
阅读本篇目录
下载完整输入、实际结果、赝势与五张高清图。解压后在已分配P100资源中执行 bash run.sh,也可用 sbatch submit.slurm提交一张GPU。新结果在 work/cases,表格在 work/analysis,结构在 work/converted,图片在 work/figures。
1 键长扫描和结构优化有什么区别#
键长扫描由你事先指定一系列距离,每个距离各做一次固定结构单点。它能帮助理解势能曲线、压缩和拉伸产生的力;结构优化则让程序根据力寻找局部极小值。
这里使用H₂,键长0.55、0.65、0.75、0.85、1.10、1.50、2.00Å。所有条件使用同一个12Å立方盒、Γ点、PBE/ONCV、60Ry和相同电子设置。
2 一个点的完整输入#
module load abacus/3.10.1-p100
cat inputs/h2-molecule/INPUT
cat inputs/h2-molecule/STRU
cat inputs/h2-molecule/KPT
程序为ABACUS3.10.1 P100版,GNU13.3、HPC-X2.50、CUDA12.8,实际1MPI×1线程×1GPU。
INPUT_PARAMETERS
suffix tutorial
calculation scf
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
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
H₂放在z方向,两个H分别位于中心下方和上方。Direct坐标乘晶格矢量后才是Å;改变键长时不能把0.75Å直接填成分数坐标0.75。
例如12Å盒中r=1.10Å,两个z分数坐标分别为:
z1 = (6−1.10/2)/12
z2 = (6+1.10/2)/12
3 运行七个点#
完整包会按次序运行所有13个教学条件,包括本篇七个点、前面的空位与二维优化。每个计算独立使用同一张已分配GPU,不需要多卡才能复现。
ls inputs/h2-*
sbatch submit.slurm
若手动运行单点,保持赝势相对位置:
mkdir my-h2
cp inputs/h2-molecule/INPUT inputs/h2-molecule/STRU inputs/h2-molecule/KPT my-h2/
cd my-h2
mpirun -np 1 abacus > stdout.log
包根目录的 PP与 my-h2并列,因此 pseudo_dir ../PP能找到H赝势。
4 检查收敛和提取能量#
grep 'charge density convergence is achieved' work/cases/h2-molecule/OUT.tutorial/running_scf.log
grep 'FINAL_ETOT' work/cases/h2-molecule/OUT.tutorial/running_scf.log
grep -A6 'TOTAL-FORCE' work/cases/h2-molecule/OUT.tutorial/running_scf.log
| H–H距离/Å | 实际总能量/eV | 上方H的Fz/eV·Å⁻¹ |
|---|---|---|
| 0.55 | -30.626945795 | 13.963073 |
| 0.65 | -31.515588678 | 4.809122 |
| 0.75 | -31.730546158 | -0.007594 |
| 0.85 | -31.588901171 | -2.550863 |
| 1.10 | -30.615944414 | -4.596230 |
| 1.50 | -28.821375926 | -4.070552 |
| 2.00 | -27.159191597 | -2.601400 |

本次离散点中0.75Å能量最低,但它只是最低采样点,不能把采样间隔当作结构优化精度。较短距离时,上方H受到+z方向力、两原子排斥;较长距离时,上方H受到−z方向力、趋向靠近。
5 提取CSV并重新画图#
module load ase/3.29.0
python -B work/analyze.py
python -B work/plot_results.py
cat work/analysis/h2-scan.csv
本篇的扫描数据在JSON的 H2_curve部分。画图使用实际总能量减去七点中的最低能量,零点平移不会改变曲线形状或力。
6 不能把这条曲线直接叫解离能#
本例采用限制性闭壳层H₂电子态;键长很大时,电子态描述需要额外检查。解离能需要另算自旋极化自由H原子参照;零点能、温度和熵修正则决定它与相应热力学量的联系。
计算反应能时,必须分别获得各反应物与产物的正确自旋和已收敛能量,按配平系数相加;计算自由能还需要相应热力学修正。这里先把结构修改、单点、力和曲线的操作流程跑清楚。