Gaussian 过渡态、正反IRC与端点优化

本篇以9原子乙醇的羟基扭转为例，连接两个相反方向的旁式构象。它是构象互变，没有断键成键。全部数值来自鸢算科技集群实际 Gaussian 16 C.01 AVX CPU 计算，B3LYP/6-31G(d)、UltraFine 积分网格；厂商预编译程序按用户独立初始化，没有公共 Gaussian module。

下载TS、正反IRC、端点优化的完整输入与脚本。在新解压目录先重画实际数据：

bash run.sh

这一步无需 Gaussian，使用 ASE 3.29 Python 环境中的 NumPy、SciPy 和 Matplotlib。输出 analysis/、figures/，拒绝覆盖旧结果。原生重算需要个人 Gaussian 环境和已分配的 CPU 作业资源；参考输入每个1线程、1GB主机内存，不作为 Intel Xeon 8173M 性能数据。

1. 整个流程先做什么、后做什么#

1. 从0°附近的乙醇几何搜索无约束扭转 TS。
2. 做频率，确认它只有1个虚频，且振动方向确实是羟基扭转。
3. 从同一个真实 TS 检查点分别做 Forward 和 Reverse IRC。
4. 从每条路径保存的实际校正端点生成新的优化输入。
5. 分别优化两端并做频率，核对稳定构象和能量，再画路径图。

在已分配 CPU 的作业环境中：

bash compute.sh reference-ts
bash compute.sh ethanol-irc-forward
bash compute.sh ethanol-irc-reverse
bash compute.sh ethanol-endpoint-forward
bash compute.sh ethanol-endpoint-reverse
bash run.sh analyze

也可在新解压目录用 bash run.sh compute 按顺序自动完成。首先 yskit gaussian-init 初始化个人环境，普通提交方法见Gaussian输入与提交。沿用56核整节点模板时，输入线程/内存与实际分配一致；不要在客户登录终端直接运行正式任务。

2. 先搜索真正的过渡态#

完整 TS 输入：

%Chk=result.chk
%NProcShared=1
%Mem=1GB
#p B3LYP/6-31G(d) Opt=(TS,CalcFC,NoEigenTest,Tight,MaxCycles=60) Freq Int=UltraFine SCF=Tight NoSymm

Ethanol unconstrained torsional saddle search

0 1
C 1.20596950 -0.40274986 0.00000000
C -0.00474869 0.52305884 0.00000000
O -1.24895047 -0.17550944 -0.00000000
H -1.06147472 -1.12475847 -0.00000000
H 0.02590867 1.17573580 0.88421084
H 0.02590867 1.17573580 -0.88421084
H 2.13979760 0.17049228 -0.00000000
H 1.20674645 -1.04476455 0.88994247
H 1.20674645 -1.04476455 -0.88994247

Opt=TS 进行过渡态搜索；CalcFC 初始求 Hessian；Tight 收紧几何收敛。NoEigenTest 不是找到正确过渡态的证明，结束后仍必须核对完整振动结果和反应方向。

grep 'Optimization completed' work/reference-ts/output.log
grep 'Frequencies --' work/reference-ts/output.log
grep 'Normal termination' work/reference-ts/output.log

本例最终电子能 -155.0316048327 Hartree，唯一虚频 -353.5689 cm⁻¹，其余20个振动为正。模式位移使 C1–C2–O3–H4 二面角向两边变化，见IR/Raman与虚频判断。

3. 正反IRC分别运行#

Forward 的完整输入：

%OldChk=../reference-ts/result.chk
%Chk=result.chk
%NProcShared=1
%Mem=1GB
#p B3LYP/6-31G(d) IRC=(RCFC,Forward,MaxPoints=12,StepSize=10) Geom=AllCheck Guess=Read Int=UltraFine SCF=Tight NoSymm

Reverse 的完整输入：

%OldChk=../reference-ts/result.chk
%Chk=result.chk
%NProcShared=1
%Mem=1GB
#p B3LYP/6-31G(d) IRC=(RCFC,Reverse,MaxPoints=12,StepSize=10) Geom=AllCheck Guess=Read Int=UltraFine SCF=Tight NoSymm

两份 %OldChk 都读同一个 reference-ts/result.chk，各自 %Chk 写本目录的新结果。RCFC 从 TS 检查点读取 Cartesian 力常数；Geom=AllCheck 读几何、电荷、多重度等，不再重复填写标题和坐标。

MaxPoints=12 是每个方向最多12个路径点；输入中 StepSize=10，本次输出显示 supplied step 为0.100 Bohr，质量加权积分步约0.1073 √amu·Bohr。它不是10°，也不是10 Å。最后横坐标直接取原生保存的反应坐标，不能用点号乘一个未经核对的单位代替。

grep 'Point Number:' work/ethanol-irc-forward/output.log
grep 'NET REACTION COORDINATE' work/ethanol-irc-forward/output.log
grep 'Reaction path calculation complete' work/ethanol-irc-forward/output.log
grep 'Maximum number of steps reached' work/ethanol-irc-forward/output.log

Point Number: 0 是 TS 起点；随后1–12是沿路径推进的点。汇总表含起点共13行，不等于做了13个新的独立反应。达到本例点数上限时，路径尚未到达极小值，所以还需要下一步端点优化。

4. 自动从自己的IRC结果生成端点输入#

bash compute.sh ethanol-endpoint-forward
bash compute.sh ethanol-endpoint-reverse
cat work/ethanol-endpoint-forward/input.gjf

compute.sh 调用 generate_endpoint.py，从当前 work/ethanol-irc-*/result.fchk 的 IRC point 1 Geometries 数组最后一帧读取坐标，转成 Å，生成新的完整 GJF。它读取用户本次路径，不直接照搬参考构象，也不需要手工抄9个原子的坐标。

本次正向路径端点生成的完整优化输入：

%Chk=result.chk
%NProcShared=1
%Mem=1GB
#p B3LYP/6-31G(d) Opt=(Tight,CalcFC) Freq Int=UltraFine SCF=Tight NoSymm

Optimize actual IRC forward endpoint

0 1
C 1.25026857 -0.35575675 0.02420178
C 0.03622821 0.56575644 -0.04391596
O -1.19890644 -0.13877281 0.02917554
H -1.13369289 -0.90132021 -0.56547430
H 0.03063518 1.25681232 0.80580902
H 0.08458949 1.17834749 -0.95992917
H 2.18336942 0.21966956 0.01310664
H 1.22594719 -0.96142135 0.93664146
H 1.27355664 -1.03634335 -0.83787340

这里是普通无约束 Opt，没有 TS 或固定二面角行。分别检查两端的优化收敛、正常结束和21个正频率：

grep 'Optimization completed' work/ethanol-endpoint-forward/output.log
grep 'Frequencies --' work/ethanol-endpoint-forward/output.log
grep 'Frequencies --' work/ethanol-endpoint-reverse/output.log

5. 提取实际路径能量与构象结果#

head analysis/irc-forward.csv
head analysis/irc-reverse.csv
cat analysis/RESULTS.json

反向 FCHK 中保存的坐标为负，正向为正；脚本保留这个符号。日志里的 NET REACTION COORDINATE 可打印为路径长度的正值，不能只按该行强行把两边画到同一侧。

两端是相反的旁式分支，电子能在1×10⁻⁸ Hartree内一致；TS相对端点的电子势垒约 7.050 kJ/mol。图中路径只画本次实际12点，两端极小值单独展示，没有把未计算的IRC段补画出来。

6. 不要把最后一条SCF当作校正后的路径点#

本版本默认 HPC 积分器会在路径点做校正。Current cartesian coordinates / Total Energy 描述最后一次 SCF/梯度评价；它们与校正后的 IRC 点可以略有差异。本例最后端点两者的原子间距最大差约0.0029 Å。

配套脚本读取路径专用字段：

formchk work/ethanol-irc-forward/result.chk forward.fchk
grep '^IRC Reference Energy' forward.fchk
grep '^IRC point' forward.fchk

IRC point 1 Results for each geome 保存本次13个状态的能量偏移和带符号反应坐标，IRC point 1 Geometries 保存对应13帧几何。字段标题在 FCHK 中有固定宽度，显示可能被截短，包内按实际字段名提取。

实际路径能量的恢复方式为 IRC Reference Energy + Results第一列；该 reference 不一定等于 TS 总能，不能未经检查就直接把第一列称作相对TS能量。路径中保存的是原生积分器的校正结果，不把它与邻近梯度评价的SCF能量硬凑成完全相同。

公开包保留相关实际数值摘录，省略程序、二进制检查点和轨道系数；完整原生重算会生成自己的 result.chk 和 result.fchk。

7. 换体系时应检查什么#

先确认 TS 唯一虚频确实对应目标反应，再检查两个方向是否通向期望结构；需要时增加 IRC 点数或调整步长，并重新优化两端。Forward/Reverse 只是相对本次模式相位的方向名称，不自动等于反应物/产物。

本例选择 Opt=TS，没有把未执行的 QST2/QST3 当作已运行的路线。能垒是相同方法下的电子势能差，不是含ZPE/热修正的自由能势垒，也不直接给出反应速率。自由能、溶剂、受阻转子、隧穿和动力学修正要按目标问题另行建模。

继续看结构优化与构象、约束与扫描、热化学和温压。