ONIOM计算(一):简介

2022-09-02 13:30:05 浏览数 (1)

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方式显示。

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

代码语言:javascript复制
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

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

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

代码语言:javascript复制
#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的必要,而是对体系重要区域使用较大基组,环境部分使用较小基组,这样的混合基组策略更为合适。

0 人点赞