Starting with Multiscale models - QM/XTB ONIOM#
从多尺度模型开始 - QM/XTB ONIOM #

Let's start our incursion on the world of multiscale calculations. Multiscale here means that we will use different levels of theory/methods for different parts of our system.
让我们开始探索多尺度计算的世界。这里的多尺度意味着我们将为系统的不同部分采用不同层次的理论/方法。

These multiscale models can be useful for a variety of problems: treating the active site of a protein on a higher level and the rest of the protein on a simpler level, modeling a chromophore on a higher level and the surrounding explicit solvent on a lower level, or even splitting a medium-sized system into smaller parts for improved efficiency.
这些多尺度模型可用于多种问题:在更高层次上处理蛋白质的活性位点,而在更简单的层次上处理蛋白质的其余部分;在更高层次上建模发色团,在更低层次上显式建模周围溶剂,甚至将中型系统分割成更小的部分以提高效率。

In ORCA, there are different approaches to run multiscale simulations, based on three main categories:
在 ORCA 中,运行多尺度模拟有不同的方法,主要基于以下三大类别:

  1. ONIOM methods mixing different "Quantum Mechanical" (QM) methods,
    ONIOM 方法混合了不同的“量子力学”(QM)方法,

  2. QM/MM methods mixing QM with "Molecular Mechanics" (MM) methods that rely on predefined force-fields, and
    QM/MM 方法将量子力学(QM)与依赖于预定义力场的分子力学(MM)方法相结合,

  3. Crystal-QM/MM methods which treat a small part of a periodic crystal (e.g. the unit cell) at a high QM level and embed this in a lower level model of the large super-cell model.
    晶体-QM/MM 方法,将周期性晶体的一小部分(如单胞)在高 QM 水平上处理,并将其嵌入大型超胞模型的较低级别模型中。

It is important to emphasize that our implementation is focused in all cases on a "QM-centric" perspective, meaning it is aimed at users familiar with QM approaches, and so will be these tutorials.
重要的是要强调,我们的实施在所有情况下都聚焦于“以量子力学为中心”的视角,这意味着它旨在面向熟悉量子力学方法的用户,因此这些教程也将如此。

Our first ONIOM single point#
我们的第一个 ONIOM 单点计算#

We will begin with a very simple approach: a QM/QM2 calculation. This means we will use two different levels of QM on different parts of the same system. In particular, let's start really simple and use QM/XTB, with the GFN2 method from Stefan Grimme (XTB2) [Grimme2019] being used together with a regular Density functional theory (DFT).
我们将从一个非常简单的方法开始:QM/QM2 计算。这意味着我们将在同一系统的不同部分使用两种不同级别的量子力学方法。具体来说,让我们从最简单的开始,使用 QM/XTB,其中结合了 Stefan Grimme 的 GFN2 方法(XTB2)[Grimme2019]与常规的密度泛函理论(DFT)。

For our test case, the Peptide...Peptide pair from the S66 benchmark set will be used [Hobza2011] (structure taken from here ). It is a complex of two simple amides, made to simulate peptide intermolecular interactions:
在我们的测试案例中,将使用 S66 基准集中的 Peptide...Peptide 对[Hobza2011](结构取自此处)。这是一个由两个简单酰胺组成的复合物,旨在模拟肽间相互作用:

../_images/pep2.png

Starting with the objective to compute only Single point energies, you will see that using QM/XTB in ORCA is really simple. All you have to do is prepare your regular QM input and just add the !QM/XTB flag, plus a list of atoms that should be treated at the higher QM level under %QMMM.
从仅计算单点能的目标出发,你会发现使用 ORCA 中的 QM/XTB 方法其实非常简单。只需准备常规的 QM 输入文件,并添加!QM/XTB 标志,以及在%QMMM 下指定应采用更高 QM 级别处理的那些原子列表即可。

In our first trial let's just separate the two molecules, treating one at revPBE level and the other at XTB2:
在我们的首次试验中,让我们仅分离这两种分子,将其中一个在 revPBE 水平上处理,另一个在 XTB2 水平上处理:

!QM/XTB revPBE D4 DEF2-TZVP
%QMMM QMATOMS {0:11} END END
* XYZ 0 1
C   -0.701502936  -0.290627698   2.406884396
H   -1.183295956   0.395647770   3.098874220
H    0.349561571  -0.030321572   2.307833035
H   -0.794056854  -1.291605451   2.824039291
C   -1.448546246  -0.244876636   1.091815299
O   -2.660450004  -0.428479088   1.034345768
N   -0.670056563   0.005916557   0.009776912
H    0.326675319   0.122563958   0.141592839
C   -1.227054574   0.089793737  -1.319967541
H   -2.292024256  -0.106501186  -1.240877562
H   -1.077801692   1.079940300  -1.748543540
H   -0.776628489  -0.647999190  -1.983372734
C    2.041774914  -2.351697968   0.686397615
H    2.599999718  -3.261701204   0.480489612
H    1.113083056  -2.358227418   0.122072198
H    1.782555985  -2.328251269   1.743338615
C    2.809410856  -1.097285929   0.350160876
O    2.264224213   0.004150876   0.293188485
N    4.136169069  -1.266099696   0.136412906
H    4.512490369  -2.193345392   0.213170232
C    5.023407253  -0.159633723  -0.152535629
H    4.409214873   0.731176053  -0.232359342
H    5.750821803  -0.020167991   0.644867678
H    5.548397546  -0.319615451  -1.091677965
*

In this case, the calculation will be such that the QM method (revPBE D4 DEF2-TZVP) will be applied to atoms 0 to 11. These atoms form the high level region. The other atoms which are not in that list (12-23) will be treated on XTB level only.
在这种情况下,计算将设定为对原子 0 至 11 应用 QM 方法(revPBE D4 DEF2-TZVP),这些原子构成高层次区域。而未列入该列表的其他原子(12-23)则仅在 XTB 层次上进行处理。

Important 重要

ORCA atom counting starts from 0!
ORCA 原子计数从 0 开始!

In a simple graphical representation, that is:
在简单的图形表示中,即:

../_images/pep2_qmmmregions.png

ONIOM output# ONIOM 输出#

Now, except for the regular ORCA output, there are the QM/QM2 specific parts:
现在,除了常规的 ORCA 输出外,还有 QM/QM2 特定的部分:

                     **************************
                     *    2-layered ONIOM     *
                     **************************
(...)
Multiscale model                       ... QM1/QM2
QM2 method                             ... XTB2
Coupling Scheme                        ... subtractive
Embedding Scheme                       ... electrostatic
(...)
Size of different subsystems (in atoms):
Size of QMMM System                    ... 24
Size of QM2 Subsystem                  ... 12
Size of QM1 Subsystem                  ... 12
(...)

There is quite some printing related to the ONIOM method, but what matters here is:
关于 ONIOM 方法的文献相当丰富,但在此重要的是:

  1. the Multiscale model is QM1/QM2, using XTB2 as QM2
    多尺度模型为 QM1/QM2,其中 QM2 采用 XTB2

  2. the default coupling scheme is subtractive. This is related to how the final energy is calculated.
    默认耦合方案为减法。这与最终能量的计算方式有关。

  3. the default embedding scheme is electrostatic. This is related to how both QM systems interact electrostatically.
    默认的嵌入方案是静电的。这与 QM 系统如何通过静电相互作用有关。

  4. the full size of the system is 24 atoms, and the size of each subsystem is 12.
    系统的总大小为 24 个原子,每个子系统的大小为 12。

Next come the calculations necessary for the subtractive scheme: full system using QM2, high level region (small system) using QM2 and high level region using QM1. First, both QM2 level calculations are printed:
接下来是减法方案所需的计算:全系统使用 QM2,高层次区域(小系统)使用 QM2,以及高层次区域使用 QM1。首先,打印出 QM2 层次的两种计算:

            ************************************
            *       QM/   XTB2 SP Energy       *
            ************************************

------------------------------------------------------------------------------
                          XTB2 - FULL SYSTEM
------------------------------------------------------------------------------
Running XTB calculation
---------------------------------   --------------------
FINAL SINGLE POINT ENERGY (L-QM2)      -33.983328801310
---------------------------------   --------------------
------------------------------------------------------------------------------
                          XTB2 - SMALL SYSTEM
------------------------------------------------------------------------------
Running XTB calculation
---------------------------------   --------------------
FINAL SINGLE POINT ENERGY (S-QM2)      -16.992757934120
---------------------------------   --------------------

and the high level QM1 calculation:
高级别 QM1 计算:

------------------------------------------------------------------------------
                                QM - SMALL SYSTEM
------------------------------------------------------------------------------

comes next, with the same printout as a regular single point calculation would have.
接下来是与常规单点计算相同的输出结果。

At the very end, we have the energy of the high level region using QM1 printed as:
最后,我们以 QM1 的形式打印出高能级区域的能量:

-------------------------   --------------------
FINAL SINGLE POINT ENERGY      -248.586744400306
-------------------------   --------------------

But this is still not our ONIOM energy! The final energy for our multiscale calculation is then given by:
但这还不是我们的 ONIOM 能量!最终的多尺度计算能量由以下公式给出:

--------------------------------------   --------------------
FINAL SINGLE POINT ENERGY (QM/QM2)          -265.577315267496
--------------------------------------   --------------------

Using ONIOM to compute binding energies#
使用 ONIOM 计算结合能

Now let's explore a bit more of the ONIOM method to compute a property. For example, let's compute a simple binding energy between these two amides. By "binding energy" here we mean simply the energy necessary to separate both molecules while keeping their geometry.
现在让我们进一步探讨 ONIOM 方法来计算一个性质。例如,我们来计算这两个酰胺之间的简单结合能。这里所说的“结合能”指的是在保持其几何结构不变的情况下,将两个分子分离所需的能量。

Our reference structure for the separate molecules will be:
我们参考的独立分子结构为:

24

C      -14.659345474      0.562404274      1.314809681
H      -15.141138494      1.248679742      2.006799505
H      -13.608280967      0.822710400      1.215758320
H      -14.751899392     -0.438573479      1.731964576
C      -15.406388784      0.608155336     -0.000259416
O      -16.618292542      0.424552884     -0.057728947
N      -14.627899101      0.858948529     -1.082297803
H      -13.631167219      0.975595930     -0.950481876
C      -15.184897112      0.942825709     -2.412042256
H      -16.249866794      0.746530786     -2.332952277
H      -15.035644230      1.932972272     -2.840618255
H      -14.734471027      0.205032782     -3.075447449
C       15.999617452     -3.204729940      1.778472330
H       16.557842256     -4.114733176      1.572564327
H       15.070925594     -3.211259390      1.214146913
H       15.740398523     -3.181283241      2.835413330
C       16.767253394     -1.950317901      1.442235591
O       16.222066751     -0.848881096      1.385263200
N       18.094011607     -2.119131668      1.228487621
H       18.470332907     -3.046377364      1.305244947
C       18.981249791     -1.012665695      0.939539086
H       18.367057411     -0.121855919      0.859715373
H       19.708664341     -0.873199963      1.736942393
H       19.506240084     -1.172647423      0.000396750

Now, when computing the energy difference EbindingONIOM=EfarONIOMEcloseONIOM we get a value of 10.25kcal/mol. This value is not too good when compared to either the pure DFT result (8.70kcal/mol) or the pure XTB value (8.03kcal/mol), which are closer to the reference value of  8.7kcal/mol [Hobza2011].
现在,在计算能量差 EbindingONIOM=EfarONIOMEcloseONIOM 时,我们得到一个值为 10.25kcal/mol 。与纯 DFT 结果( 8.70kcal/mol )或纯 XTB 值( 8.03kcal/mol )相比,这个值并不理想,因为它们更接近参考值  8.7kcal/mol [Hobza2011]。

Selection of the proper QM region#
选择合适的 QM 区域

This error actually is not intrinsic from the ONIOM process, but rather stems from a poor choice of the QM regions for this particular problem. Since we want to ultimately compute the energy related to the breaking of the hydrogen bond, which is our strongest binding point here, a more sensible choice of spaces would be to select the atoms directly involved in binding as our high level region:
此误差实际上并非源自 ONIOM 过程本身,而是源于对该特定问题 QM 区域的不当选择。由于我们最终目的是计算氢键断裂相关的能量,这是此处最强的结合点,因此更合理的做法是选择直接参与结合的原子作为高级别区域:

../_images/pep2_qmmmregions_2.png

Which in terms of the input simply translates to a slightly different input:
这意味着在输入方面仅需要稍作不同的输入:

! QM/XTB revPBE D4 DEF2-TZVP PAL16
%QMMM QMATOMS {4:7 16:19} END END
%MAXCORE 3000
* XYZ 0 1

And a resulting binding energy of EbindingONIOM=8.41kcal/mol, which deviates only by 0.3kcal/mol from the reference value.
并且得出的结合能为 EbindingONIOM=8.41kcal/mol ,仅与参考值偏离 0.3kcal/mol

Of course, this is a rather small example, but the same would apply for a much larger system. As long as the regions are chosen properly, the multiscale methods can greatly accelerate calculations with a still relatively small error!
当然,这是一个相当小的例子,但对于一个更大的系统,情况同样适用。只要区域选择得当,多尺度方法就能在保持相对较小误差的同时,极大地加速计算!

Comparison of different methods to compute the binding energy of the Peptide...Peptide cluster from the S66 set.
S66 数据集中肽簇结合能计算的不同方法比较。
#

Method 方法

Ebinding kcal/mol  Ebinding 千卡/摩尔

pure DFT 纯密度泛函理论

8.70

pure XTB 纯 XTB

8.03

naive ONIOM 朴素 ONIOM

10.25

proper ONIOM 适当的 ONIOM

8.41

Reference [Hobza2011] 参考文献 [Hobza2011]

8.72

Important 重要

Always choose your QM region such that it reflects the chemical phenomenon that you want to study!
始终选择你的 QM 区域,使其反映出你想要研究的化学现象!

Extra content# 额外内容 #

Here are some extra comments on the ONIOM or QM/MM schemes. A nice reference related to these can be found in the review by Lili Cao and Ulf Ryde [Ryde2018].
以下是对 ONIOM 或 QM/MM 方案的一些额外评论。关于这些方法的一个很好的参考资料可以在 Lili Cao 和 Ulf Ryde 的综述中找到[Ryde2018]。

What is "subtractive ONIOM"?#
什么是“减法 ONIOM”?

There are two major approaches to compute QM/MM or ONIOM energies: additive and subtractive. In the subtractive scheme, the energy is computed as:
计算 QM/MM 或 ONIOM 能量的方法主要有两种:加和法与减和法。在减和法方案中,能量计算公式为:

EQM1/QM2sub=E1QM1+E12QM2E1QM2

where the subscript refers to the system where the method is applied: 1 corresponds to high level region, 2 to low level region and 12 to the full system. The subtractive scheme can be applied for QM1/QM2 methods and QM/MM methods.
其中下标指代应用该方法的系统:1 对应高能级区域,2 对应低能级区域,12 对应整个系统。减法方案可应用于 QM1/QM2 方法和 QM/MM 方法。

In the additive scheme, which can only be used for QM/MM, the energy is computed as:
在加和方案中,仅适用于 QM/MM,能量计算如下:

EQM/MMadd=E1QM+E2MM+E1/2QM/MM

with the last term corresponding to the interaction energy term between systems 1 and 2.
最后一个项对应于系统 1 和系统 2 之间的相互作用能项。

In ORCA, the default for the ONIOM method is the subtractive scheme, while the default for QM/MM is the additive scheme.
在 ORCA 中,ONIOM 方法的默认设置为减法方案,而 QM/MM 方法的默认设置为加法方案。

What is "electrostatic embedding"?#
什么是“静电嵌入”?

When splitting a large system into two or more smaller ones we need to decide how to properly include the interaction between them? There are two different schemes to compute the electrostatic interaction between high and low level region: "mechanical embedding" and "electrostatic embedding".
将一个大系统分割成两个或多个较小的系统时,我们需要确定如何恰当地包含它们之间的相互作用?计算高低层次区域间静电相互作用有两种不同方案:“机械嵌入”和“静电嵌入”。

In the simpler mechanical embedding scheme, the electrostatic interaction of both systems is just considered using the lower QM2 or MM level method, without accounting for the polarization of the high level region by the low level region at the higher level QM method.
在较为简单的机械嵌入方案中,仅通过较低的 QM2 或 MM 级方法考虑两个系统的静电相互作用,而未在高级别 QM 方法中计入低级别区域对高级别区域的极化效应。

In the electrostatic embedding scheme (ORCA's default), the higher level calculation is done in the presence of the charges of the lower level region. This second approach is more accurate, for it lets the high level region polarize in the presence of the low level region on the higher level calculation.
在静电嵌入方案(ORCA 的默认设置)中,高层次计算是在低层次区域的电荷存在的情况下进行的。第二种方法更为精确,因为它允许高层次区域在高层次计算中因低层次区域的存在而极化。