ONIOM计算（一）：简介

ONIOM（our Own N-layer Integrated molecular Orbital molecular Mechanics）方法是由已故量子化学家Keiji Morokuma等人提出的一种用于大体系计算的杂化方法，它可以对体系的不同部分用不同的计算级别处理，例如对化学上比较关心的成键、断键区域用较高水平的计算方法，对周围环境用较低水平的方法。如果两者分别为量子力学QM和分子力学MM方法，这时候就是一种QM/MM方法。当然，在ONIOM方法中，可以是任意的方法的组合，两层可以都是QM方法。原理上来说，ONIOM方法可以将体系分成n层，而实际计算中一般将体系分成2层或3层。以下我们以两层ONIOM（记为ONIOM2）方法为例，简单说明ONIOM方法的原理，关于三层ONIOM以及更详细的方法原理可参看相关文献，如DOI: 10.1002/wcms.85

在ONIOM2中，整个体系被称为Real system，用高水平方法描述的部分（下图黄色部分）称为Model部分，而蓝色部分则用低水平方法描述：

下图是Morokuma在2000年ACS年会上报告的一张幻灯片，介绍了两层ONIOM方法的能量计算公式：

上图可以这样理解：我们的目标是得到高水平下整个体系的能量，与低水平下的Model部分对比，两者的差别在于两方面：一是体系的大小，二是计算的级别。因此，以低水平下Model部分的能量为起始，通过外推的方法得到近似的高水平下Real体系的能量。图中，E(low, real)−E(low, model)为体系大小的外推项，E(high, model)−E(low, model)为计算级别的外推项，整理之后即得到上图中右下部分的能量表达式。类似的思想也用在计算大基组下CCSD(T)能量，见《电子能量的基组外推以及ORCA中的自动实现》一文。

ONIOM2的能量公式也可以从另一个角度理解，将其写成：

E(ONIOM) = E(high, model) + E(low, real) − E(low, model)

其中第一项为高水平部分的能量，后两项是低水平下总体系的能量减去model部分的能量，这部分并不能简单地理解为低水平部分的能量，实际还近似地包含了高水平区域和低水平区域之间的相互作用。

Gaussian、Q-Chem、ORCA等程序均支持ONIOM方法。本文通过一个例子来展示如何用GaussView来设置输入文件，以及用Gaussian做ONIOM计算，后续我们还会推送更多关于ONIOM方法的应用实例。

该例是Exploring 3上的例2.9，计算了维生素E分子的甲基被氧化成亚甲基的两种异构体的能量差。

down构型

up构型

该分子含79个原子，在现在看来只能算是小分子范畴，实际的应用中完全没有必要用ONIOM方法，此处仅用之作为示例。我们将左边结构差异部分（结构式黑色和红色部分）作为高层，而右边的烷基链作为低层，因为这部分在两个异构体中应当是相似的结构。

ONIOM的输入文件可借助GaussView来完成。构建完分子后，选择Tools →Atom groups，打开下图所示的操作界面：

将Atom Group Class选为ONIOM Layer，则出现下方的表格，默认为三层ONIOM，如果是两层，则一层为High，另一层为Medium或Low都行。本例中，我们构建两层ONIOM。在默认状态下，所有的原子都被归到了High层。此时，我们选中低层部分的原子，如下图所示：

在Group Actions下拉菜单中选择Add Selected Atoms to “ONIOM Layer (Medium)”，则显示如下：

低层原子显示为蓝色高亮，点击OK，则得到下图：

这是GaussView默认的显示方法，高层用球棍模型，低层用Wireframe，如果有中层，则中层为Tube方式显示。

最后，保存坐标。此处我们展示坐标的前几行，如下：

C      0  -16.10401900   0.12571900   -0.41400100 H
C      0  -15.59713500  -1.12390400   -0.62002900 H
C      0  -14.27035200  -1.34352000   -0.49793200 H
C      0  -13.42160900  -0.34615800   -0.16252800 H
C      0  -13.92535100   0.88913800    0.03805900 H
C      0  -15.24543400   1.13243100   -0.08871700 H
O      0  -12.07871200  -0.62683300   -0.08528100 H
C      0  -11.15273600   0.41022100   -0.16937700 H
C      0  -11.59238300   1.50349500    0.79970400 H
C      0  -12.98561800   2.01182700    0.42650700 H
C      0   -9.79930900  -0.19816700    0.25643900 L H 8
C      0  -11.07995800   0.90071300   -1.62226700 H
C      0   -8.57973600   0.73902200    0.19849800 L
C      0   -7.29537600   0.03322800    0.66993500 L

在ONIOM的输入文件中，每行坐标的标准格式为

atom [freeze-code] coordinate-spec layer [link-atom [bonded-to [scale-fac1 [scale-fac2 [scale-fac3]]]]]

其中方括号表示可省。各项的含义如下：

• atom：元素符号
• freeze-code：0表示在做几何结构优化时不被冻结，-1表示冻结，若是其他负数，则标有相同负数的原子在结构优化时作为刚性片段。
• coordinate-spec：原子的XYZ坐标
• layer：原子所属的层，H、M、L分别表示高、中、低层
• link-atom：链接原子。当分层涉及到切断化学键时，需要给高（或中）层的悬键进行饱和，常用氢原子，避免高层出现非正常的电子结构。
• bonded-to：该原子在高层中所连接的原子的序号。本例中，低层的11号原子连接了高层的8号原子，则在计算高层时，11号原子会被替换成H原子。
• scale-fac：高层的原子与链接原子之间的键长会相对原始键长进行校正。本例中，8号与11号原子之间本为C-C键，实际在计算高层时，则为C-H键，C-H的键长要比C-C键短，因此需要乘上校正因子，使计算结果更合理。这些值一般不需要指定，程序内部会有内置的值。

本例中，如果我们对高层使用B3LYP/6-31G(d,p)级别，低层使用PM6半经验方法，则关键词写为ONIOM(B3LYP/6-31G(d,p):PM6)，也就是用冒号将不同层的计算级别隔开。

在电荷和自旋多重度部分，对两层ONIOM方法，其书写格式为

Charreal,low Spinreal,low Charmodel,high Spinmodel,high Charmodel,low Spinmodel,low

三套电荷和自旋多重度实际对应了计算公式中的三项。

综上，我们对该体系进行结构优化，并做频率计算，以比较两个异构体的自由能的差别，则输入文件如下：

#p opt freq oniom(b3lyp/6-31g(d,p):pm6)

[No Title]

0 1 0 1 0 1
C      0  -16.10401900   0.12571900   -0.41400100 H
C      0  -15.59713500  -1.12390400   -0.62002900 H
C      0  -14.27035200  -1.34352000   -0.49793200 H
C      0  -13.42160900  -0.34615800   -0.16252800 H
C      0  -13.92535100   0.88913800    0.03805900 H
C      0  -15.24543400   1.13243100   -0.08871700 H
O      0  -12.07871200  -0.62683300   -0.08528100 H
C      0  -11.15273600   0.41022100   -0.16937700 H
C      0  -11.59238300   1.50349500    0.79970400 H
C      0  -12.98561800   2.01182700    0.42650700 H
C      0   -9.79930900  -0.19816700    0.25643900 L H 8
C      0  -11.07995800   0.90071300   -1.62226700 H
C      0   -8.57973600   0.73902200    0.19849800 L
C      0   -7.29537600   0.03322800    0.66993500 L
[坐标略]

限于篇幅，此处我们不展示输出文件的内容。读者可练习此例，阅读输出文件，并与前面的公式进行对照，加深理解。最后，自由能也是从Sum of electronic and thermal Free Energies=后面读取出来，与普通的opt freq计算一样。本例由于体系不大，我们也用DFT直接做了计算，结果对比如下：

ONIOM(B3LYP/6-31G(d,p):PM6)

G(up): -655.001051 au

G(down): -654.996503 au

ΔG: 2.85 kcal/mol

B3LYP/6-31G(d,p)

G(up): -1283.889184 au

G(down): -1283.885332 au

ΔG: 2.42 kcal/mol

ONIOM方法分层注意事项

（1）高层区域不能太小，可以将高层区域定为化学活性中心（如成键/断键）区域向外延伸三个化学键以上。

（2）断键（此处指不同层之间的截断）的位置尽可能在单键位置，最好在非极性单键位置。尽量不要将双键、三键、配位键以及芳环等切断。

小 结

本文简单介绍了ONIOM方法的基本原理，并展示了一个简单的例子。后面我们还会推送更多关于ONIOM的应用实例，如在计算复杂体系化学反应中的应用。对于低层，本文的例子使用的是半经验量子化学方法，而关于使用分子力场的情况，我们以后再介绍。ONIOM方法有很高的普适性，已经得到了广泛的应用。它除了可用于分子基态的研究，也可用于激发态的计算。除了ONIOM方法外，线性标度算法是另一类可以用于大分子体系计算的方法，但是由于很多线性标度算法还不够“黑箱”，加之梯度等性质的计算相对不易，因此ONIOM方法依然有较大的用武之地。但是由于现在计算机的计算能力发展迅速，ORCA等程序中的RI等技术使得DFT方法已经能较轻松地处理含有三五百个原子的体系，因此，笔者认为对于此类中等大小（也可称之为大体系了）的体系，没有太多使用ONIOM的必要，而是对体系重要区域使用较大基组，环境部分使用较小基组，这样的混合基组策略更为合适。