CP2K 过渡态与虚频检查
NH3平面鞍点精修、解除约束的真实受力、一个虚频与Molden模式。
阅读本篇目录
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 |

本例恰有一个虚频,约 −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手册。