VTST：100K instanton两阶段实操

本篇接续instanton路径输入和完整矩阵读取，给出一个从路径到前因子输出的完整终端流程。模型仍为10Å晶胞内固定三个H、允许N三个方向移动的NH₃，4个计算半图像镜像成8珠子；本次温度为100K。

1. 下载并先读取实际结果#

下载全部输入、两阶段原生输出、矩阵与脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd vtst-instanton-two-stage
bash read.sh
bash plot.sh

默认读取附带结果，不重新运行VASP。CSV在 analysis-read，图片在 figures。

本例使用 vasp/6.6.1-8173m-vtst CPU构建，Intel oneAPI环境。提交模板按4MPI×14CPU申请本站56核整节点，实际每MPI使用1个线程，计费按整节点规则。8珠子是虚时间路径的离散结构数，不是8个真实MD时间点。

2. 第一阶段：只搜索路径#

先看本例完整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=100
Qdo_ins=.TRUE.
Qdo_pre=.FALSE.
Qftol=0.02

温度 Qtpz=100，电子收敛 EDIFF=1E-6，路径力分量标准 Qftol=0.02。这一阶段只打开 Qdo_ins，关闭 Qdo_pre。

python3 prepare.py
cat /你的势目录/H/POTCAR /你的势目录/N/POTCAR > work-new/path/POTCAR
cd work-new/path
sbatch ../../submit.sh
squeue -u "$USER"
# 等作业结束后
grep 'ut ' 01/insout.dat

资料包不提供POTCAR。势文件顺序须与POSCAR的 H N 对应。实际搜索约83秒，3次力求值，路径最大力分量为 0.0013855862 eV/Å，低于设置的0.02。

原子约束、目录、KPOINTS和MODECAR含义见路径输入篇。本例初始坐标来自已有真实能量和力数据构造的初猜；最终数字仍由本次VASP求得。

3. 第二阶段：单独计算矩阵和前因子输出#

不要把两个开关同时设为真。 当前实现会直接调用两个过程，不能把它们理解成“程序会自动先搜完路径再进入下一阶段”。

本例初始基准路径第一轮即满足力标准，准备脚本会确认这一条件，再复制未被Lanczos有限位移扰动的输入结构：

cd ../..
python3 prepare-prefactor.py work-new/path
cd work-new/prefactor
sbatch ../../submit.sh
squeue -u "$USER"

新INCAR中的关键差别是：

Qdo_ins = .FALSE.
Qdo_pre = .TRUE.

其余温度、珠子数、晶胞、势文件与有限位移保持对应。本例该阶段约245秒，每个图像有7次电子结构求值，形成24×24的完整矩阵。

不要把最后一次有限位移探测的CONTCAR直接当作驻点路径。如果用户改动后的路径经过多轮移动，资料包的小例准备脚本会停止并提示，需先提取对应的驻点基准结构再做矩阵。

4. 读取新算出的结果#

cd ../..
bash read.sh work-new
bash plot.sh
head analysis-read/VTST-instanton-full-eigenvalues.csv
cat analysis-read/DATA-VALIDATION.json

提交脚本同时记录运行退出码与耗时，读取器还检查实际OUTCAR和矩阵；仅有作业号或退出码为0，不能代替结果检查。

本例实际前五个无量纲特征值：

-0.29319458
 0.03237942
 0.68327793
 1.70680672
 1.91984103

100K例子有1个负方向，第二项为正的小值；120K例子有2个负值。两者来自实际不同温度计算，不能把其中一个温度的谱或前因子直接搬到另一个温度。

5. 程序打印的前因子怎样理解#

100K实际打印约 3.96927×10⁶ s⁻¹。读取脚本用实际矩阵特征值、质量、温度和程序常数独立复算，矩阵最大差约 3.1×10⁻⁸，特征值最大差约 1.6×10⁻⁸。

这个实现仍会把第二个模式按零模近似替换后计算前因子。因此这里给的是该离散模型和代码近似下的程序输出，没有给出反应速率，也没有声称珠子数、温度、路径精度或零模处理已经达到定量研究要求。固定H的模型也不能直接代替自由NH₃全原子速率。

6. 换成自己的体系#

先确定结构、约束和反应坐标，准备合理初猜；先运行路径搜索，检查路径力和结构；提取对应的基准路径；再单独运行矩阵阶段；检查完整谱、单位和零模，再考虑更细的珠子数及速率参照。

数据下载：实际基准路径、路径求值记录、完整谱、矩阵有限位移记录。