ORCA的简单使用之SCF计算

2022-04-18

在最开始，我想首先申明：
理论计算不是算命，理论计算也不能算出“我们所想要的值”！
理论计算是独立但又与实验相辅相成的存在，能提供我们一种独特的视角去观察一个现象。
理论计算没有巧合，遵循规则，大家都能重复出相近甚至一样的数值。

ORCA的简单介绍

ORCA目前是一款学术免费使用的量子化学计算程序，是继流行度最高的Gaussian的后起势力。这两款量子化学程序也是各有长处，ORCA得益于RI技术在稍大体系（150原子左右）的双精度杂化泛函计算中，能以忽略不计的误差实现较快的计算速度，这是Gaussian所不能及的，但Gaussian因为对泛函的解析梯度做得好，收敛算法也较强，在振动分析方面的稳健性还是胜于ORCA的。
作为一款量子化学计算程序，ORCA已经能实现大部分最为常用计算任务并拥有了不少好用的新特性，目前ORCA的版本已经来到了5.0.3版本，但据说5.0.3版本存在一些小bug，现使用得较多的还是5.0.2版本，本教程中使用的所有例子若不带说明均采用ORCA 5.0.2做计算。

常见的计算任务

量子化学计算中常见的计算有5种，分别是：单点能计算（SCF）、几何优化（OPT）、振动分析（FREQ）、过渡态搜索（TS）以及分子动力学（MD）。

单点能计算（SCF）

单点能又称电子能量，不同的计算程序计算得到的电子能量可能存在差异，一般不能直接进行比较，但相对的能量是可以进行比较的。在KS方程的DFT计算中，对电子能量的描述分为已知和未知两部分。
已知的部分可以由体系内原子的坐标配合所用的基组（Basis set）直接计算得到，由四种能量：

电子的动能
电子与原子核间的库仑作用力贡献的势能
电子与电子间的库仑作用力贡献的势能
原子核与原子核之间的库仑作用力贡献的势能
未知的部分只有一种能量：电子的交换相关势能，其值可以利用电子的密度来做交换相关函数的密度函数近似获取，这也是为什么称DFT计算的原因。不同的泛函都是为了更精确地处理这个值而设计的。

几何优化（OPT）

几何优化是为了获取结构在物理上具有意义的优化结构，这个优化的结构有两个特点：在势能面局部受力最小、能量最低。这里一定要注意是局部！因为对于稍大一些的体系，其势能面可能会存在非常多的凹面，这些凹面的最低点也就是所谓的优化结构所在的位置，这也是一个分子会具有如此多构象的原因，尤其是柔性大分子。一般而言，为了使得优化的结构更有物理意义，收敛限一般会设置得较单点能计算更严格一些。因为要计算最低的能量，所以几何优化的过程中每一步的优化过程都是包含单点能计算的，所以对于大分子体系，若初始的结构设置的不好，会让几何优化任务的耗时大大增加。这有个小技巧：可以通过半经验的方法如PM7或者xTB程序对初始结构做预优化，使得初始结构率先掉入势能面上离初始位置的凹面，来减少后续DFT计算几何优化的时间。

势能面说明

振动分析（FREQ）

振动分析是为了获取体系原子的间正振动的力常数，这是通过求解Hessian矩阵的本征值得到的。可以将其简单理解为获取热力学函数及分析原子振动的频率及振动模式的计算任务，在论文中常见的台阶图以及红外光谱等均可以通过该任务计算得到。对于几何优化得到的优化结构，其振动分析的结果中不应该出现虚频。

过渡态搜索（TS）

过渡态搜索是研究反应物初态末态之间的过渡态任务，对于研究反应机理过程有着很重要的作用。与几何优化得到的凹面点相对，过渡态搜索是为了获取势能面上相邻两个凹点间凸起的一阶鞍点。通过对过渡态结构进行振动分析，有且仅有一个虚频出现。通过与反应初态结构的热力学能量相减，可以获取反应所需的势能垒，结合公式可以用于分析反应过程中间的决速步。

分子动力学（MD）

一般而言，利用量子化学程序计算分子动力学的轨迹，又可以称为AIMD（ab initio Molecular Dynamics）。AIMD最大的用处就是在有限的时间内模拟化学反应过程、获取分子不同的构象初猜以及统计分析。AIMD比普通力场的MD而言精度要高好几个数量级，但是对计算资源和时间而言非常昂贵，一般模拟的时间在50 ps以内，而力场的MD利用GPU加速计算，能短时间来到ns级别的模拟。这部分比较有趣，有兴趣可以多参考一些资料进行学习。

使用ORCA进行简单的单点能计算

以下算例可以点击这里获取。以下计算均使用ORCA 5.0.2进行，输入文件均利用Multiwfn软件生成并稍作修改，使用与详细介绍可以看看sob老师的博文。

甲醇的单点能计算

甲醇单点能计算的输入文件内容如下：

! B3LYP D3 def2-TZVP def2/J RIJCOSX noautostart miniprint nopop
%maxcore  1000
%pal nprocs  6 end
* xyz   0   1
C     -2.04044573    2.68778863    0.02088746
H     -1.66508826    1.66542794   -0.00088209
H     -1.67995641    3.16820148    0.93760543
H     -3.13552932    2.65298710    0.04665737
O     -1.55170419    3.32975501   -1.14898890
H     -1.86780487    4.23863570   -1.15774844
 *

简单介绍一下输入文件的内容：

1	! B3LYP D3 def2-TZVP def2/J RIJCOSX noautostart miniprint nopop // 关键词行，可以有多列，以"!"为首。具体写法参考ORCA的官方手册

1
2
3

// 对计算参数的精细控制，以"%"为首，"end"结束
%maxcore  1000 // 每个核心可利用的内存大小（MB）
%pal nprocs  6 end // 本次计算采用的核心数，通常小分子建议6个核心以内，多了容易报错，计算大分子时，有多少物理核心就上多少物理核心，根据服务器资源进行合理调控

// 两个"*"围着的部分就是对体系原子的描述，包括坐标系、净电荷数目、电子自选多重度、原子的坐标
* xyz   0   1 // "xyz"为采用笛卡尔坐标系，"0"代表当前体系下净电荷数为0，"1"代表当前体系的电子自旋多重度为1，以下为对应的原子坐标描述
C     -2.04044573    2.68778863    0.02088746
H     -1.66508826    1.66542794   -0.00088209
H     -1.67995641    3.16820148    0.93760543
H     -3.13552932    2.65298710    0.04665737
O     -1.55170419    3.32975501   -1.14898890
H     -1.86780487    4.23863570   -1.15774844
 *

电子自选多重度的计算方式为：
$$自选多重度 = (\alpha电子数-\beta电子数)+1$$

将以上内容复制的到一个文本文件中，将其命名为methanol_scf_TZVP.inp，随后将其丢入服务器中，在该文件的当前目录执行下述代码完成甲醇的单点能计算。

1	orca methanol_scf_TZVP.inp \| tee methanol_scf_TZVP.out

计算结果很快就出来了，在服务器上使用vim编辑器查看一下运算的结果，输入如下命令：

1	vim methanol_scf_TZVP.out

......
 -------------------------------------------------------------------------------
                           DFT DISPERSION CORRECTION
                               DFTD3 V3.1  Rev 1
                           USING Becke-Johnson damping
 -------------------------------------------------------------------------------
 -------------------------   ----------------
 Dispersion correction           -0.003103441
 -------------------------   ----------------
 -------------------------   --------------------
 FINAL SINGLE POINT ENERGY      -115.712688760143
 -------------------------   --------------------
......
 Timings for individual modules:
 Sum of individual times         ...        2.532 sec (=   0.042 min)
 GTO integral calculation        ...        0.317 sec (=   0.005 min)  12.5 %
 SCF iterations                  ...        2.215 sec (=   0.037 min)  87.5 %
                              ****ORCA TERMINATED NORMALLY****
 TOTAL RUN TIME: 0 days 0 hours 0 minutes 2 seconds 800 msec

我摘取了单点能计算中关键的结果部分，一般在输出文件最后出现 ‘ORCA TERMINATED NORMALLY’ 即可说明计算任务正常完成了。
其中最关键的是 ‘FINAL SINGLE POINT ENERGY’ 的值，也就是甲醇的电子能量，以及运算完成后在当前目录下生成的methanol_scf_TZVP.gbw文件，这是记录了单点能计算后体系的二进制波函数文件，可以利用下述命令进行转化：

1 2	orca_2mkl methanol_scf_TZVP -molden // 一般使用文件名即可，不用加后缀，即orca_2mkl [文件名] -molden

运行后会提示

1 2	Reading the input file methanol_scf_TZVP.gbw ... done Writing the output file methanol_scf_TZVP.molden.input ... done

如此便得到了methanol_scf_TZVP.molden.input波函数文件，之后可以利用sob老师开发的Multiwfn软件进行各式各样的波函数分析，以获取更为丰富的数据进行深层次的可视化分析。

乙烯的单点能计算

为区分上述和甲醇例子的区别，这里用sob老师在Multiwfn中推荐的几档计算来计算乙烯的单点能量(从左往右，由低到高)，并让大家对高精度计算的耗时有个直观的感受。

项目	B3LYP(D3)/def2-SVP	B3LYP(D3)/def2-TZVP	PWPB95(D3)/def2-QZVPP	DLPNO-CCSD(T) with tightPNO	CCSD(T)/cc-pVTZ	CCSD(T)/CBS (cc-pVTZ->QZ extrapolation)
核数	6	6	12	12	12	12
时间$(s)$	1.935	2.800	7.245	20.345	9.523	147.619