CP2K专题

CP2K 过渡态与虚频检查

NH3平面鞍点精修、解除约束的真实受力、一个虚频与Molden模式。

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

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 对称性帮助精修候选鞍点#

脚本选取收敛路径的最高图像,将四个原子的z设置到同一平面z=6 Å,再固定z、只优化x/y。这利用了本例NH₃翻转的平面对称性;一般反应不一定有这个平面,不能不加判断地把同样约束用于任意反应。

完整input.inp:saddle-refine#

&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 5.999872153900 5.999999999500 6.000000000000
   H 7.004749255500 5.999999998300 6.000000000000
   H 5.497432407200 6.870255648800 6.000000000000
   H 5.497432404900 5.129744351500 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
 &CONSTRAINT
  &FIXED_ATOMS
   LIST 1 2 3 4
   COMPONENTS_TO_FIX Z
  &END FIXED_ATOMS
 &END CONSTRAINT
&END MOTION

LIST 1 2 3 4指四个原子的1-based编号,COMPONENTS_TO_FIX Z仅用于平面精修。满足约束下的GEO_OPT停止条件后,还不能只凭它认定完整三维驻点。

3 去掉约束,再检查全方向受力#

脚本从最终XYZ生成独立 RUN_TYPE ENERGY_FORCE,其中不再包含MOTION约束;这样可以检查三个方向的真实受力,而不是被约束投影后写成零的分量。

完整input.inp:saddle-static#

&GLOBAL
 PROJECT inversion
 RUN_TYPE ENERGY_FORCE
 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.999869105800 5.999999999800 6.000000000000
   H 7.004736359900 6.000000001500 6.000000000000
   H 5.497438566900 6.870244944500 6.000000000000
   H 5.497438571400 5.129755052000 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
grep -A 7 'FORCES| Atomic forces' work/cases/NH3-inversion/saddle-static/output.log
cat work/analysis/RESULTS.json

本文最终最大原子力矢量模约 0.000320 eV/Å,三个方向一起参与检查。

4 全自由度振动#

振动阶段同样没有原子约束,使用DX=0.005 Bohr。四个原子做24个正负方向位移,分析非线性分子的3N−6=6个内部模式。

完整input.inp:saddle-vibrations#

&GLOBAL
 PROJECT inversion
 RUN_TYPE VIBRATIONAL_ANALYSIS
 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.999869105800 5.999999999800 6.000000000000
   H 7.004736359900 6.000000001500 6.000000000000
   H 5.497438566900 6.870244944500 6.000000000000
   H 5.497438571400 5.129755052000 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
&VIBRATIONAL_ANALYSIS
 DX 0.005
 NPROC_REP 1
 FULLY_PERIODIC .FALSE.
 &PRINT
  &MOLDEN_VIB ON
  &END MOLDEN_VIB
 &END PRINT
&END VIBRATIONAL_ANALYSIS
grep -E 'VIB\|Frequency' work/cases/NH3-inversion/saddle-vibrations/output.log
cat work/analysis/saddle-frequencies.csv
ls work/cases/NH3-inversion/saddle-vibrations/*.mol
内部模式 频率/cm⁻¹
1 -787.4702
2 1525.4604
3 1525.7576
4 3562.1867
5 3775.8021
6 3775.8308

真实Molden虚频模式的位移箭头和六个内部频率

本例恰有一个虚频,约 −787.47 cm⁻¹,另外五个内部模式为正。Molden特征向量显示氮与三个氢沿平面法向反向运动,三个氢的位移方向一致,对应伞形翻转。箭头仅按统一比例放大显示,不是又做了一次受扰动几何优化。

5 怎样判断自己的结果#

先检查解除约束后的驻点受力,再看内部虚频数量,最后用模式位移判断是否对应目标反应。只有“一个虚频”仍不够:如果虚频沿着别的键断裂或不希望的构型变化,应重新检查路径和候选结构。小虚频还需要结合SCF、优化精度和有限位移步长判断。

inversion-VIBRATIONS-1.mol保留原生模式,可用支持Molden振动格式的查看器播放。本例附件中的图直接读取该特征向量;没有用插值动画代替实际频率计算。

反应速率或自由能势垒还需要额外的端点振动、标准态、统计模型等信息;本例只给出电子路径、驻点受力与谐振动性质。稳定分子的热化学流程见气相RRHO热修正,不能把过渡态虚频按稳定分子正频率直接处理。

提取与重画#

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 VIBRATIONAL_ANALYSIS手册。

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

查看配图

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

100%查看原图 ↗