CP2K专题

CP2K NEB反应路径与端点准备

NH3翻转的真实端点、七图像CI-NEB、完整输入与电子能垒。

更新于 2026-10-28约 9 分钟阅读按任务与环境选择适用配置
阅读本篇目录

本篇用气相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,后面不再额外平移。

实际NH3端点和最高路径图像

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优化器。

需要根据你的体系选择资源?联系讨论计算需求。实际可用软件与权限,以当前账户和运行环境为准。

查看配图

放大后可滚动查看细节,按 Esc 关闭。

100%查看原图 ↗