VTST：instanton完整矩阵与零模检查

instanton路径搜索中的Lanczos只寻找低曲率方向，不能代替完整振动谱。本篇在同一个已达到0.02路径力标准的120K、8珠子NH₃模型上，实际计算完整矩阵，演示如何看零模和前因子的适用条件。

1. 下载并读取实际矩阵#

下载输入、原生矩阵、28次电子结构输出及绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd vtst-instanton-full-matrix
bash read.sh
bash plot.sh
head reference-results/01/eigvals.fyl

默认只读取已有结果。实际文件为：

本例三个H固定，N的x、y、z方向自由，共3个原子坐标自由度；完整8珠子对应 3×8=24 维矩阵。它是经过质量与虚时间尺度变换并包含弹簧的矩阵，不能直接把特征值当作普通分子振动频率的cm⁻¹。

2. 运行完整矩阵所需的输入#

VASP版本仍为 vasp/6.6.1-8173m-vtst，Intel oneAPI构建。与路径搜索的关键差别是：

Qdo_ins = .FALSE.
Qdo_pre = .TRUE.
QdR = 0.005
IMAGES = 4
Qtpz = 120
EDIFF = 1E-6

Qdo_ins 为假时直接围绕提供的路径做有限位移，不再搜索路径。因此要先准备合理的驻点路径。每个半路径图像计算1个基准点和3个自由坐标的正负位移，即7次求值；4个图像合计28次电子结构求值。本例约240秒，速度只供本小模型参考。

完整INCAR：

SYSTEM=NH3 fixed-H constrained N inversion ring-polymer teaching case
ENCUT=400
PREC=Accurate
EDIFF=1e-6
EDIFFG=-0.02
ISPIN=1
ISYM=0
ISMEAR=0
SIGMA=0.05
NELM=160
ALGO=Normal
LREAL=.FALSE.
LASPH=.TRUE.
IBRION=3
POTIM=0
NSW=80
ISIF=2
ICHAIN=4
IMAGES=4
IOPT=0
NCORE=1
LWAVE=.FALSE.
LCHARG=.FALSE.
Qltol=0.01
QdR=0.005
Qmaxmove=0.02
Qdt=0.05
Qnl=4
Qifcg=.TRUE.
Qtpz=120
Qdo_ins=.FALSE.
Qdo_pre=.TRUE.
Qftol=0.02

计算目录、POSCAR、KPOINTS和MODECAR均在 inputs 中。使用自己的许可势文件，H在前、N在后：

python3 prepare.py
cat /你的势目录/H/POTCAR /你的势目录/N/POTCAR > work-new/POTCAR
cd work-new
sbatch ../submit.sh
# 等作业结束后
cd ..
bash read.sh work-new
bash plot.sh

3. 用NumPy直接核对矩阵和特征值#

module load ase/3.29.0
python - <<'PYCODE'
import numpy as np
H = np.loadtxt('reference-results/01/premat.fyl')
e = np.loadtxt('reference-results/01/eigvals.fyl')
print('matrix shape:', H.shape)
print('max asymmetry:', np.max(abs(H-H.T)))
print('lowest five:', e[:5])
print('negative modes:', np.count_nonzero(e < 0))
print('matrix/eigenvalue difference:', np.max(abs(np.linalg.eigvalsh(H)-e)))
PYCODE

实际前五项为：

-0.49521142
-0.00610734
 0.20106512
 1.50478966
 1.52625564

脚本还利用全部正负位移的实际原子力重建矩阵；与原生矩阵最大差约 3.7×10⁻⁸，特征值最大差约 1.8×10⁻⁸。可下载完整谱CSV和各次位移求值CSV。

4. 为什么不能直接引用打印的前因子#

这组粗离散路径出现两个负特征值，其中第二个 −0.00610734 较接近零。它仍不等于已得到严格的零模或已满足定量速率要求。

当前VTST实现会对第一个特征值取相反数，并把第二个特征值替换为1后计算前因子。程序因此打印约 2.85×10⁷ s⁻¹，但本篇不把这个数作为可靠物理前因子，更不是反应速率。它缺少完整的路径、珠子数和零模收敛分析，也没有完整配分函数参照。

做客户自己的体系时，先确认结构与约束、路径驻点和完整谱，再增加珠子数并检查有限位移与电子精度。仅看到“正常退出”或有限的前因子数字，不足以判断速率可信。本例的实际用途，是展示完整矩阵怎么计算、怎么读，以及如何识别还需完善的近零模。