ORCA的几何优化以及振动分析
不同于高精度单点能的计算，几何优化及振动分析并不需要很大的基组，可以参考sob老师的博文（浅谈为什么优化和振动分析不需要用大基组），普遍来说量子化学计算的共识是利用相同的小基组配合适当的普通泛函做几何优化和振动分析，再采用大基组配合对应的普通的泛函乃至双杂化泛函进行高精度单点能的计算。这里说的小基组至少是2-zeta基组，大基组普遍是3-zeta以上的大基组，关于sob老师对基组的选择（谈谈量子化学中基组的选择）以及泛函的选择（简谈量子化学计算中DFT泛函的选择），具体内容可以参考对应的博文。所有计算用到的基组一定要加上极化！！！具体如何使用含极化表达的基组参考ORCA 5.0.2 Manual p.28-p.34，当然可以直接使用Multiwfn生成的ORCA输入文件，sob老师推荐的几档泛函和基组组合可以直接参考使用，当然更推荐大家都了解一些基组和泛函的基本知识，这对于做合规的基本的量子化学计算以及判断文献计算的方法是否可取是非常重要的。
知己知彼，百战不殆
对于简单的有机物计算，一般可以采用B3LYP(D3-BJ)杂化泛函配合def2-SVP或def-TZVP进行几何优化和振动分析，这对200个原子的有机物体系也错错有余。下边我贴出了科研室里可用的计算资源：
- 2 × 28核（Platinum 8173M）
- 6通道386GB内存（12 × 32GB）
- 4TB机械硬盘（存计算数据专用）
56核对于普通的DFT计算以及双杂化泛函计算都是能算得动的，而386GB的内存也保证了每个核心至少有6.8GB的内存资源可以使用，日常使用基本不用担心爆显存和算力不够，请放心大胆的使用！
因为几何优化与振动分析通常是一块进行的，所以将放在一块进行说明，当然单独计算也是可以的，只需要对应写上事宜的关键词即可。
本文算例中的文件点击此处下载获取。
甲醇的几何优化及振动分析
甲醇属于有机体系的分子，原子数目仅为6个，这里直接利用B3LYP泛函辅以DFT-D3色散校正，采取3-zeta中大型基组def-TZVP进行结构优化与振动分析。输入文件内容如下：
1 | ! B3LYP D3 TZVP def2/J RIJCOSX opt freq tightSCF noautostart miniprint nopop |
与之前在进行单点能计算时不一样，因为ORCA默认计算是单点能计算，此时的输入文件的关键词行多了opt与freq两个关键词，并且未做进一步精细控制的情况下，ORCA会首先对结构进行优化，优化至Energy（体系能量变化）、gradient（体系受力梯度）、step（与上次结构的位移变化）变化量都在限定的阈值内即几何结构优化收敛，此时结构应该落在了势能面上的某个极小点的位置了，这样能在输出文件里看到下边的提示：
1 | ...... |
几何结构优化结束后，因为关键词中还有freq关键词，表示还要进行解析振动分析（AnFreq），对普通的DFT泛函而言，ORCA目前都是支持梯度解析的，少数不支持的可以使用数值解析（NumFreq，也是一个关键词），而数值解析是很贵的计算，对超过20个原子以上的体系要慎重使用，一般情况下，ORCA5版本对常规振动分析计算任务也不需要进行数值解析，吃力不讨好。振动分析主要是对Hessian矩阵进行解析，之后可以得到振动频率、红外光谱、热力学量等数据：
1 | ...... |
可以从VIBRATIONAL FREQUENCIES的数据里直接看到，振动分析并未出现虚频，这说明甲醇的基态几何优化正确，结构与数据可以采用。
这里我还量了一下计算得到的甲醇结构与CRC Handbook of Chemistry and Physical里的物性数据进行了比对：
可以看到，偏差很低，除键角外基本在3%以内，可以说B3LYP结合DFT-D3在def-TZVP基组下对甲醇结构的描述很到位。
苯酚的几何优化与振动分析
在建模软件中，建好了苯酚的模型，并利用Multiwfn生成相应的输入文件phenol_opt_freq.inp。几何结构优化结果如下：
1 | ...... |
表示收敛都是YES，但是振动分析发现存在虚频：
1 | ----------------------- |
这时用VMD软件看一看优化的结构，VMD的初始载入配置设置可以参考sob老师的vmd.rc配置，用起VMD来也比较方便一些。
可以从图中看到苯酚基本没有怎么优化，用软件量了一下$\angle COH$，角度为109.46$^o$，C-O键长为1.390埃，O-H键长为0.963埃，与文献参考值差别也不大。
不过体系中存在虚频，这也说明该体系可能不稳定，即能量并不足够低。
Significant negative frequencies indicate saddle points of the energy hypersurface and prove that the optimization has not resulted in an energy minimum.
(ORCA 5.0.2 Mannual. p.198)
为此利用Grimme组开发的好用的xTB程序预优化了一下，使用的是xTB-GFN2方法，使用如下命令进行：
1 | xtb phenol_opt.xyz --opt --gfn 2 |
过程很快，不到2 s就完成了，此时在文件夹中得到的xtbopt.xyz即为优化好的苯酚结构，xtbopt.log是优化结构的轨迹文件，添加后缀.xyz即可用VMD观看优化结构的轨迹，这里直接取xtbopt.xyz，利用Multiwfn程序生成输入文件phenol_after_xTB_ORCA_opt_freq.inp供ORCA进行计算，摘取了部分结果：
1 | ----------------------- |
这次没有虚频出现了，最终的苯酚的几何结构是这样的：
再来分别看一下两次结构优化后的能量大小差距：
1 | phenol_opt_freq.out:FINAL SINGLE POINT ENERGY -307.410282918115 |
差了0.00587111582701 Eh的能量大小，这是非常大的差距，换算成kcal/mol就约有3.684 kcal/mol的差距，这比量子化学中常用的1 kcal/mol粗略能量差距判据足足大了3.6倍之多！可见最初所得的几何优化结构并不是很稳定，也不足以作为苯酚的另一种较稳定的构象，所以苯酚最终的稳定几何结构应该就是第二次通过xTB预优化再利用ORCA优化得到的结构，这与文献的结构相比也差距不大。
结构优化的其他说明
当然几何优化在不收敛的情况下计算任务也会正常结束，但这并不意味着结果可以用，一定要判断清楚是什么原因导致的几何优化不收敛，通常从两方面判断：
- 正在逐步收敛，但是优化次数达到了优化的步长上限值
- 化学直觉，对物质可能结构的判断
对于第一种情况，最有可能的是初猜结构离附近能量最低的结构过于遥远，可以尝试调整初猜结构，不过对这种情况最好的方式就是续算，即用上次计算结束的最后一帧优化结构进行接续计算，波函数的初猜也可以用最后一次得到的波函数文件。
对于第二种情况，一般直接根据可靠的文献参考以及化学直觉来重新设置初猜结构即可。
如果啥都不行，可以利用xTB做一下100 ps的分子动力学获取多帧的初猜结构，再进行构象搜索，这对单一物质的几何优化，尤其是柔性大分子而言是非常有效的方式（参考sob老师对瑞德西韦的构象搜索例子，之后我也会对之前计算的多元醇酯基础油做算例详解）