CP2K NEB反应路径与端点准备
NH3翻转的真实端点、七图像CI-NEB、完整输入与电子能垒。
阅读本篇目录
本篇用气相NH₃的伞形翻转演示反应路径:氮原子从三个氢构成的平面一侧翻到另一侧。两端是同一分子的两个等价构型,中间经过近似平面结构;所有原子保持一致编号,中间图像允许自由调整。
下载完整输入、实际路径和日志、模式与绘图脚本,解压进入run.sh目录。本篇使用CP2K2026.2双精度P100构建:GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3;分析为ASE3.29.0。
1 提交完整流程#
sbatch submit.slurm
模板申请七张P100、七个MPI进程,每进程一个CPU核、一个OpenMP线程,NPROC_REP=1;不同图像/位移分配到不同GPU。在已分配这些资源后也可以执行 bash run.sh。这是核时P100模板,套餐用户用yskit生成所属分区脚本,保留资源分配后的运行命令;不要在登录节点计算。
参考结果中的鞍点精修与振动采用同版本CPU构建(Intel oneAPI2024.2、Intel MPI2021.13、MKL),端点与最终静态力采用P100构建;公开脚本已整套在P100复跑。这里比较科学流程与数值,不作CPU/GPU性能排名。
脚本顺序是“分别优化两端 → 放在共同参考位置 → 在此位置重新优化两端 → 插值七幅图像 → CI-NEB → 平面过渡态精修 → 解除约束的力检查 → 全自由度振动 → 提取与画图”。inputs可编辑,reference是本文实际结果,work为你的新计算;已有work时重新解压,避免覆盖。
2 先把端点准备好#
采用PBE、DZVP-MOLOPT-SR-GTH、N的GTH-PBE-q5和H的GTH-PBE-q1,12 Å非周期MT盒、400/50 Ry网格,EPS_SCF1E−8。中性分子有8个价电子,采用闭壳层基态。端点GEO_OPT的最大力目标1E−5 Hartree/Bohr,约0.000514 eV/Å。
完整input.inp:endpoint-up#
&GLOBAL
PROJECT inversion
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../../DATA/GTH_POTENTIALS
&MGRID
CUTOFF 400
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&END DFT
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&COORD
N 6.000000000000 6.000000000000 6.400000000000
H 6.950000000000 6.000000000000 6.000000000000
H 5.525000000000 6.822724133595 6.000000000000
H 5.525000000000 5.177275866405 6.000000000000
&END COORD
&KIND N
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q5
&END KIND
&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&MOTION
&GEO_OPT
OPTIMIZER BFGS
MAX_ITER 100
MAX_FORCE 1.0E-5
RMS_FORCE 5.0E-6
MAX_DR 1.0E-4
RMS_DR 5.0E-5
&END GEO_OPT
&END MOTION
另一端保留相同元素顺序与电子参数,把初始氮位置改到氢平面的另一侧;完整文件为inputs/endpoint-down/input.inp。不要在两个端点中交换等价氢的编号,否则插值可能制造无意义的原子交叉。
脚本把各自优化后的质心放到(6,6,6) Å,再在该位置分别重新优化。有限数值网格下,平移后的结构可能出现小残余力,不能直接把平移前的收敛结论搬过来。第二次优化后的端点坐标直接用于BAND,后面不再额外平移。

3 插值、弹簧与爬山图像#
脚本从两个实际端点按相同原子编号线性插值,生成image-00.xyz至image-06.xyz,共七幅图像,包含两端。NUMBER_OF_REPLICA 7不是七个独立反应,也不代表一个图像有七个MPI进程。
完整input.inp:band#
&GLOBAL
PROJECT inversion
RUN_TYPE BAND
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../../DATA/GTH_POTENTIALS
&MGRID
CUTOFF 400
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&END DFT
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&COORD
N 5.999831213100 6.000000001400 6.079774698000
H 6.942928647200 5.999999996200 5.688640762600
H 5.528373181400 6.816689639300 5.688408835700
H 5.528373177900 5.183310358700 5.688408850100
&END COORD
&KIND N
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q5
&END KIND
&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&MOTION
&BAND
BAND_TYPE CI-NEB
NUMBER_OF_REPLICA 7
NPROC_REP 1
K_SPRING 0.05
ALIGN_FRAMES .FALSE.
ROTATE_FRAMES .FALSE.
&CI_NEB
NSTEPS_IT 2
&END CI_NEB
&CONVERGENCE_CONTROL
MAX_FORCE 3.5E-4
RMS_FORCE 2.0E-4
MAX_DR 0.003
RMS_DR 0.0015
&END CONVERGENCE_CONTROL
&OPTIMIZE_BAND
OPT_TYPE DIIS
OPTIMIZE_END_POINTS .FALSE.
&DIIS
MAX_STEPS 80
NO_LS .TRUE.
N_DIIS 4
STEPSIZE 1.0
MAX_STEPSIZE 0.5
&END DIIS
&END OPTIMIZE_BAND
&CONVERGENCE_INFO ON
&END CONVERGENCE_INFO
&PROGRAM_RUN_INFO ON
&END PROGRAM_RUN_INFO
&REPLICA
COORD_FILE_NAME image-00.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-01.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-02.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-03.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-04.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-05.xyz
&END REPLICA
&REPLICA
COORD_FILE_NAME image-06.xyz
&END REPLICA
&END BAND
&END MOTION
BAND_TYPE CI-NEB启用爬山图像;K_SPRING 0.05维持路径图像分布。ALIGN_FRAMES FALSE与 ROTATE_FRAMES FALSE保持本例已经准备好的共同坐标框架;两端固定,内部图像调整。NPROC_REP 1每个副本环境用一个MPI,七MPI可以同时处理七幅图像。
DIIS的 STEPSIZE 1.0是优化器步长比例,并不是1 Å的原子位移。最大位移和力还分别通过CONVERGENCE_CONTROL判断;本例最大力分量阈值3.5E−4 Hartree/Bohr,约0.0180 eV/Å,同时要求RMS力及位移收敛。
4 查看路径是否结束,提取能垒#
grep -E 'STEP NUMBER|RMS FORCE|MAX FORCE|BAND TASK|Band task converged' work/cases/NH3-inversion/band/output.log
cat work/cases/NH3-inversion/PATH.json
cat work/analysis/path.csv
原生收敛路径XYZ的注释中Energy使用Hartree,必须转换后才能写eV;analyze.py已经处理。原生日志中的 Band task converged 与力、位移停止条件应一起查看。

本例收敛CI-NEB最高图像相对端点A的电子能垒约 0.221897 eV,随后独立精修平面鞍点得到约 0.221897 eV。参考路径采用四MPI/四P100分配七图像,完整公开脚本采用七MPI/七P100;它们使用相同电子模型和收敛条件。路径初猜会影响优化过程,最终仍要检查收敛和鞍点性质。
这里是有限基组和盒长下的气相PBE电子能垒,没有加入端点/过渡态零点能差。它不代表反应速率、隧穿分裂或溶液活化自由能。继续看过渡态与虚频确认最高点是哪一种运动。
提取与重画#
module load ase/3.29.0
python analyze.py work
python plot_results.py work
cat work/analysis/path.csv
cat work/analysis/saddle-frequencies.csv
ls work/figures
脚本读取原生收敛标记、最终路径XYZ、独立静态受力和Molden模式,生成CSV、JSON与三幅图。只查看本文参考,把work换成reference。七幅图像的路径距离是全部原子3N笛卡尔位移的累计长度;它不是时间,也不是单独一根键的键长。
输入定义见CP2K2026.2 BAND手册与DIIS优化器。