过渡态优化
过渡态的搜索算法有很多种，sob老师对过渡态结构优化算法做了较为全面的综述博文（过渡态、反应路径的计算方法及相关问题），有对优化算法有兴趣的可以多多了解一下这方面的知识。在物理化学书中也经常介绍到，对于一个反应，存在着一个能垒，但是通过理论计算得到的能垒并不是反应所需的活化能！！！活化能是通过计算得到的自由能垒转换成动力学常数k并建立k与温度的关系得到的，也可以通过实验值的拟合，这一点务必要注意！（谈谈如何通过势垒判断反应是否容易发生）
不过不是所有反应都具有能垒，比如一些低温自发反应和自由基反应
对于计算所得的能垒，位置是势能面上的相邻两个极小点间的一阶鞍点位置，而过渡态优化任务就是为了找到这个一阶鞍点的位置。而在过渡态的位置做振动分析，会得到有且仅有一个虚频的结果，如结果不存在虚频，说明找的不是过渡态。
ORCA 5.0.2 Manual p.181
之后会通过呋喃和乙烯的DA加成反应作为算例进行过渡态优化的过程。
呋喃与乙烯的DA加成反应
呋喃与乙烯的DA加成反应，起始和终止位置的结构取自xTB在线手册的Reaction Path Methods小节本算例所用到的文件点击此处进行下载。
使用xTB进行过渡态搜索
首先使用xTB进行本例的计算，初始结构文件为start.xyz，终态（产物）结构文件为end.xyz，反应路径搜索的输入文件为path.inp。将三个文件放在同一目录下，使用以下命令调用xTB程序计算：
1 | xtb start.xyz --path end.xyz -I path.inp |
运算结果很快就出来了，大约10 s。在当前目录下xtbpath_ts.xyz为过渡态结构文件，xtbpath.xyz为反应路径轨迹文件，把xtbpath.xyz丢入VMD查看轨迹，首先载入xtbpath。xyz文件，随后在VMD的命令窗口输入bt，然后回车键入，这样就可以及时更新连接键的信息了，方便观看反应过程轨迹：
使用ORCA进行过渡态搜索
这里用两种算法来进行过渡态搜索，效果一样，但是耗费的时间相差甚巨！
optTS(EF算法)
ORCA默认的optTS关键词用的是EF算法，这是一种基于反应物或产物来搜索过渡态的算法，会沿着振动矢量方向前进寻找过渡态，即使是在稳定态也可以找寻过渡态，这对未知过渡态时较有用，但十分费时。本身我也不太喜欢直接从稳定态出发搜索，一是太慢，二是不一定能获得过渡态。所以本例中采用xTB优化得到的过渡态作为初猜进行过渡态搜索。
使用Multiwfn载入xtbpath_ts.xyz：
1 | Multiwfn xtbpath_ts.xyz |
之后依次键入如下内容：
oi // 转换为ORCA的输入文件
DA_reaction_optTS_ORCA.inp // 输入输入文件的文件名
0 // 进入计算任务选择界面
5 // 选择过渡态搜索+振动分析任务
4 // 选择计算使用的B3LYP-D3泛函，def2-TZVP基组（这一步无所谓，因为之后需要修改基组）
q // 退出Multiwfn
随后目录下就出现了DA_reaction_optTS_ORCA.inp文件，使用vim编辑器对其进行修改：
1 | vim DA_reaction_optTS_ORCA.inp |
按下键盘上的i进入编辑模式，将文件内容改为如下：
1 | ! B3LYP D3 TZVP def2/J RIJCOSX optTS freq tightSCF noautostart miniprint nopop |
修改完毕后按键盘上的Esc键退出编辑模式，随后在英文输入法下键入:wq保存退出文件。
这里主要是将基组从def2-TZVP更改为了TZVP（即def-TZVP），这里还多了一行%geom Calc_Hess true end，这是因为EF算法是根据振动矢量方向进行优化的，故在正式计算前需要获取Hessian矩阵的精确解。
随后进行orca的正常运算就行：
1 | orca DA_reaction_optTS_ORCA.inp | tee DA_reaction_optTS_ORCA.out |
运算约3 min结束，算上xTB的时间不过3.5 min。
1 | ...... |
观察一下振动分析的结果：
1 | ...... |
有且仅有一个虚频，说明过渡态找对了，为了验证是否真的是过渡态，之后还需要进行IRC的验证。
IRC验证
同样用Multiwfn载入上述optTS优化得到的过渡态结构文件DA_reaction_optTS_ORCA.xyz。
1 | Multiwfn DA_reaction_optTS_ORCA.xyz |
依次输入下述内容：
oi // 转换为ORCA的输入文件
DA_reaction_optTS_IRC_ORCA.inp // 输入输入文件的文件名
0 // 进入计算任务选择界面
4 // 选择振动分析任务
4 // 选择计算使用的B3LYP-D3泛函，def2-TZVP基组（这一步无所谓，因为之后需要修改基组）
q // 退出Multiwfn
随后用相同的泛函与基组设置IRC任务，同样用vim编辑器进行修改相应内容：
1 | vim DA_reaction_optTS_IRC_ORCA.inp |
按键盘上的i键进入编辑模式，修改后的DA_reaction_optTS_IRC_ORCA.inp输入文件内容如下：
1 | ! B3LYP D3 TZVP def2/J RIJCOSX freq IRC tightSCF noautostart miniprint nopop |
修改完毕后按键盘上的Esc键退出编辑模式，随后在英文输入法下键入:wq保存退出文件。
这里主要在第一行的关键词行添加了IRC关键词，并做了相应的IRC计算精细设置，因为默认的IRC走步为20步（即MaxIter 20），太短的步子可能不能使IRC走至收敛，为此稍微延长一下IRC的步子，这里设置的是50步（即MaxIter 50），对于这个反应路径已经错错有余了。
IRC的精细调控参照ORCA 5.0.2 Manual p.768-p.769
然后正常调用ORCA进行计算：
1 | orca DA_reaction_optTS_IRC_ORCA.inp | tee DA_reaction_optTS_IRC_ORCA.out |
运算约33 min完毕，验证时间还是比较长的。
1 | ...... |
我摘录了部分IRC计算输出文件的内容：
1 | ...... |
IRC的过程主要是对过渡态位置向前走至反应物或产物结构，直至收敛或超出IRC上限步长；随后向后走至反应物或产物结构，直至收敛或超出IRC上限步长，然后汇总成反应物到产物的反应路径。一般而言，挑选出Step和E两列数据做图分析即可，图像一定是过渡态能量高于两侧的反应物与产物侧的平稳能量：
随后再来看一下跑了IRC获得的反应路径两端是不是和start.xyz和end.xyz相近：
可以看到近乎是一致的，只是B3LYP-D3与def-TZVP算出来的结构和能量要更精准一些。
NEB-TS
本例还使用了一种过渡态优化的方法NEB，此法也是一种基于反应物与产物的过渡态搜索方法。说实话一般我是不推荐上来就做NEB的，因为此法耗时极高，尤其是NEB，因为这个算法在设计时就是为了更好的利用多核心的优势来做优化，ORCA默认设置的NEB在反应物与产物间的插点数为8，在计算过程中如果任务使用的核心多于8个，但不及8的整数倍，那么ORCA只会调用8个核心，分别使用1个核心算一个点，这个过程是同步进行的，所以对大一些的体系，若初始设定给每个核心的内存量不够的话，程序容易崩溃，此例仅仅只有15个原子，基组也只是中大型的def-TZVP，所以1000MB的内存还是错错有余的。
NEB-TS的精细调控及其他关键词参考ORCA 5.0.2 Manual p.780-p.785
本例的NEB-TS的输入文件DA_reaction_NEB-TS_ORCA.inp内容如下：
1 | ! B3LYP D3 TZVP def2/J RIJCOSX NEB-TS freq tightSCF noautostart miniprint nopop |
NEB-TS的输入文件是以反应物的输入文件为基础，在%neb的调控中设置产物Product的几何结构文件（.xyz）即可。
因为在使用NEB-TS计算之前，已经使用ORCA对初始的结构和终态的结构做了同水平的几何优化，这样在NEB-TS的输入文件中，可以不设置预优化参数PreOpt（默认是false，不写就完了），如果偷懒懒得分别优化反应物与产物的结构，可以在%neb的调控中加上这一句PreOpt true即可，这样在NEB-TS计算之前会对反应物与产物分别优化至极小点再进行过渡态的搜索。
也是一样正常调用ORCA进行计算：
1 | orca DA_reaction_NEB-TS_ORCA.inp | tee DA_reaction_NEB-TS_ORCA.out |
计算共花费了44.5 min，非常之耗时。
中间的过程就不看了，直接看最后得到的过渡态结构的振动分析：
1 | ...... |
与之前optTS采用的EF算法得到的过渡态结构的振动分析相比，两者的虚频相差不大，说明两种方法找到的过渡态差别不大，都可以用。这里也简单看一看NEB-TS优化得到的反应路径：
这里就不重新跑IRC了，感兴趣的可以参照之前的IRC例子跑一遍，NEB-TS获得的最终过渡态结构见文件DA_reaction_NEB-TS_ORCA.xyz或DA_reaction_NEB-TS_ORCA_NEB-TS_converged.xyz。
自由能垒计算
既然都跑出了过渡态的结构，这里顺便再分析一下呋喃和乙烯加成的自由能垒大小。关于能垒的说明，可以看看sob老师的这篇博文（谈谈如何通过势垒判断反应是否容易发生），这里严谨一些，讨论自由能垒。
自由能属于热力学统计量，因为此前在计算时都分析过反应物、产物、过渡态结构的振动分析，所以直接读取文件的中的自由能数值与自由能校正值。
这里用grep命令来直接获取吉布斯自由能的值：
1 | grep "Final Gibbs free energy" *.out |
返回：
1 | DA_reaction_end_ORCA.out:Final Gibbs free energy ... -308.46056087 Eh |
自由能校正值：
1 | grep "G-E(el)" *.out |
返回：
1 | DA_reaction_end_ORCA.out:G-E(el) ... 0.09999187 Eh 62.75 kcal/mol |
一般来说，振动分析所用的基组和泛函算得的电子能量较差，一般需要更进一步使用精度更高的单点能计算来获取，这里统一采用双杂化泛函PWPB95，以及高档的4-zeta基组def2-QZVPP。输入文件直接通过Multiwfn生成，以反应物的优化结构作为例子进行演示：
1 | Multiwfn DA_reaction_start_ORCA.xyz |
依次输入：
oi // 转换为ORCA的输入文件
DA_reaction_start_scf.inp // 输入输入文件的文件名
7 // 选择计算使用的PWPB95-D3双杂化泛函，def2-QZVPP基组
q // 退出Multiwfn
DA_reaction_start_scf.inp输入文件内容如下：
1 | ! PWPB95 D3 def2-QZVPP def2/J def2-QZVPP/C RIJCOSX tightSCF noautostart miniprint nopop |
然后直接调用ORCA进行计算：
1 | orca DA_reaction_start_scf.inp | tee DA_reaction_start_scf.out |
整个过程不过1.5 min，把剩下的DA_reaction_end_scf.inp与DA_reaction_ts_scf.inp分别算好后，同样使用grep命令来获取体系的电子能量：
1 | grep "FINAL SINGLE POINT ENERGY" *scf.out |
返回：
1 | DA_reaction_end_scf.out:FINAL SINGLE POINT ENERGY -308.556269309249 |
取此前的自由能校正值与高精度计算得到的电子能量进行相加即可得到体系在298.15 K，1 atm下的最终自由能大小：
如果想计算其他温度压力下的自由能，可以试试sob老师开发的Shermo程序，使用也是比较灵活，具体可以参考程序的手册，手册中也详细阐述了热力学量计算的方法，也是值得我们好好学习的材料。
当然也可以通过拟合阿伦尼乌斯公式获取活化能及反应速率常数，参考sob老师的这篇博文（谈谈如何通过势垒判断反应是否容易发生
），这样谈论一个反应是否容易发生就变得有理有据了，如果还需要深入分析，可以进一步利用Multiwfn程序对每个状态做感兴趣的波函数分析，得以做更为深入的机理讨论。