Christoph
Plett
a,
Abylay
Katbashev
a,
Sebastian
Ehlert
b,
Stefan
Grimme
*a and
Markus
Bursch
*c
aMulliken Center for Theoretical Chemistry, Universität Bonn, Beringstr. 4, 53115 Bonn, Germany. E-mail: grimme@thch.uni-bonn.de
bMicrosoft Research AI4Science, Evert van de Beekstraat 254, 1118 CZ Schiphol, The Netherlands
cMax-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Mülheim an der Ruhr, Germany. E-mail: bursch@kofo.mpg.de
First published on 22nd June 2023
The computational treatment of large molecular structures is of increasing interest in fields of modern chemistry. Accordingly, efficient quantum chemical approaches are needed to perform sophisticated investigations on such systems. This engaged the development of the well-established “Our own N-layered integrated molecular orbital and molecular mechanics” (ONIOM) multi-layer scheme [L. W. Chung et al., Chem. Rev., 2015, 115, 5678–5796]. In this work, we present the specific implementation of the ONIOM scheme into the xtb semi-empirical extended tight-binding program package and its application to challenging transition-metal complexes. The efficient and broadly applicable GFNn-xTB and -FF methods are applied in the ONIOM framework to elucidate reaction energies, geometry optimizations, and explicit solvation effects for metal–organic systems with up to several hundreds of atoms. It is shown that an ONIOM-based combination of density functional theory, semi-empirical, and force-field methods can be used to drastically reduce the computational costs and thus enable the investigation of huge systems at almost no significant loss in accuracy.
This encouraged the continuous development of force field (FF) and semi-empirical quantum mechanical (SQM) methods that are routinely applied to systems of hundreds up to thousands of atoms.18–21 More recent developments in this field are the GFNn-xTB and force-field methods,22 GFN1-xTB,23 GFN2-xTB,24 and GFN-FF,25 which became widely used and well-established tools for the treatment of large systems.26,27 In contrast to most other semi-empirical approaches and FFs, the GFN methods are consistently parameterized for all elements up to Rn and are thus applicable to a large chemical space throughout the periodic table. These methods are already routinely used in black-box approaches for conformer sampling,28 energetic sorting of structure ensembles,29 or docking algorithms such as the recent aISS30 method with the latter being well-suited for molecule-binding to large structures.14 Nevertheless, complex and challenging electronic structures often require an even more accurate direct or post-processing description at a more sophisticated DFT or WFT level.31
To close this gap, various multi-layer schemes have emerged. They rely on the fact that treating the whole system highly accurately is often unnecessary, but rather a small reactive part or interaction site is of importance.32 For multi-layer schemes, computationally demanding but accurate methods are applied only to a small region of interest while non-negligible environmental effects are treated with computationally cheaper methods. Crucial environmental influences include, e.g., backbone strain in proteins,33 supramolecular stabilization of reactive species,34,35 reaction mechanisms including large structurally strained systems,36 or explicit solvation.37,38 Besides well-known QM/MM schemes,39 the so-called “Our own N-layered integrated molecular orbital and molecular mechanics” (ONIOM) method32,40 proved to be a valuable tool for multi-layer simulations. ONIOM was already successfully applied for modeling organic reactions such as Friedel–Crafts reactions,41 metal–organic catalysis,42 and zeolite reactivity to only name a few.43 In this context, the GFN methods represent a perfect match to simulate the outer region within the ONIOM scheme due to their broad parameterization and reasonable accuracy.19,44–46
Accordingly, we present here the implementation of the ONIOM scheme in the open source xtb program that allows the direct combination of the GFN methods with any desired DFT or WFT method via an interface to the prominent ORCA47,48 or TURBOMOLE49,50 program packages in an easily applicable and intuitive way. This enables the efficient and accurate treatment of not only organic structures but also of molecules containing, e.g., transition metals, heavy main group elements, or even lanthanoids.19,46 Besides the theoretical excursion into the ONIOM approach and its implementation in xtb, we demonstrate its efficient application by investigating molecular geometries and energies of challenging metal–organic structures.
EONIOM = Ehighinner + Elowouter | (1) |
The low-level description of the outer region is obtained by the subtraction of the energy of the inner region (Elowinner) from the energy of the whole system (Elowwhole) calculated with the corresponding low-level approach:
Elowouter = Elowwhole − Elowinner | (2) |
In this implementation, no electrostatic embedding is employed for the calculation of the inner region.
If requested, solvation effects can be accounted for by a simplified solvation embedding. Thereby, only Elowwhole is computed embedded in an implicit solvent model while the inner region is always computed in the gas phase.
The general idea behind the link atom saturation technique is the capping of the cleaved bonds to fill the incomplete valence shells of the interface atoms of the partitioned inner region. This is done via the introduction of a dummy atom, usually Hydrogen, on the bond vector between the boundary atoms. The exact position of a link atom is derived as:
LA = inner + k(outer − inner) | (3) |
Fig. 1 The fragmentation and subsequent saturation procedures within the ONIOM framework. Capping hydrogen atoms are highlighted in blue. |
(4) |
ONIOM = lowwhole − lowinnerJ(inner; whole) + highinnerJ(inner; whole) | (5) |
The formation of the Jacobian matrix requires the differentiation of eqn (3) with respect to the real coordinates (whole).
All DFT calculations employ Ahlrichs' def2 Gaussian type atomic orbital basis sets63 and some modified variants developed in the framework of the respective DFT-3c composite methods PBEh-3c64 (def2-mSVP), B97-3c65 (def-mTZVP), and r2SCAN-3c66 (def-mTZVPP). Matching def2 effective core potentials67 were used throughout. If not stated otherwise, the default calculation settings and convergence criteria of the respective codes were applied including the use of the resolution of the identity (RI) approximation with matching auxiliary basis sets.68 The D469–71 London dispersion correction was used throughout.
The systematic workflow for locating transition states was adopted from the work of Dohm et al.72 The double-ended Growing-String Method (GSM) by Zimmerman et al. was used to refine reaction paths73–75 at the GFN2-xTB potential energy surface (PES). The resulting saddle points were further optimized with corresponding density functionals to obtain final transition state geometries. All calculations were executed on an Intel® Xeon® CPU E5-2660 v4 @ 2.00 GHz machine.
As a demonstration, the following examples show the use of our implementation for geometry optimizations of challenging metal–organic systems. To allow for a fair assessment of the ONIOM results, the showcase systems were selected to still be treatable at a reasonable DFT level. Thus, a direct comparison of the multi-layer combinations and full DFT is possible. As the size of the investigated systems limits the choice of suitable methods, the TPSS76 functional in combination with the def2-SVP basis set was chosen as the DFT approach. Nevertheless, employing the ONIOM scheme would allow for the treatment of much larger structures, where a full DFT treatment even at this low level would be unfeasible.
As a first example, a 484 atoms large cutout of the prominent zirconium-based UiO-66 MOF, synthesized by Cavka et al.,77 was examined. Various combinations of GFN and DFT methods were tested within the ONIOM scheme (Fig. 2) for optimizing the geometry. For the ONIOM calculations, the Zr-containing nodes were treated with the high-level method. The resulting method combinations are denoted as “ONIOM(high-level method:low-level method)” in the following. The influence of the inner region size on the ONIOM performance was tested by including one, three, and six Zr nodes. In general, the included Zr nodes were treated as one single extended inner region, not as separated subregions. As a measure of structural quality, the heavy-atom root-mean-square deviation (hRMSD) from the crystal structure is used.
Despite having the largest deviation from the reference, the atomistic GFN-FF is highly efficient in terms of computational timings (5 s). It performs well for the purely organic linker and is still reasonable for the Zr-moiety. The semi-empirical GFN2-xTB shows an improved structure compared to GFN-FF by still remaining computationally efficient (3 min). Applying the ONIOM(GFN2-xTB:GFN-FF) scheme yields a better hRMSD compared to a full GFN-FF optimization already with only one Zr moiety in the inner region. By extending the GFN2-xTB region to the other Zr moieties, the geometry can be improved systematically while still being computationally very efficient. Applying GFN2-xTB to each of the six Zr nodes leads to a slightly lower hRMSD compared to the full GFN2-xTB optimization. The combination of GFN2-xTB for the metal-containing moieties and GFN-FF for the purely organic linker matches the strengths of both methods and yields a very reasonable geometry within seconds up to a few minutes.
The overall lowest hRMSDs are observed for the two DFT methods (TPSS and r2SCAN-3c), but with these methods, the geometry optimization takes several hours up to more than a week for r2SCAN-3c. The combination of DFT methods and GFN-FF does not yield improved hRMSDs over the ONIOM(GFN2-xTB:GFN-FF) approach for the inclusion of one or three Zr nodes. This may be due to asymmetrical distortions of the overall structure which is not the case when six nodes are included in the inner region. Then, the ONIOM(DFT:GFN-FF) combination yields excellent agreement with the reference and the full-DFT treatments. The computation time is reduced from 46 to only 11 hours at no loss in accuracy in case of employing TPSS. The reduction in computational time is even larger for the ONIOM(r2SCAN-3c:GFN-FF) combination (191 hours to 21 hours). Accordingly, the ONIOM approach can save days to weeks of computation time at no relevant loss in accuracy. Noteworthy here is the lower computational time for the inner region including six nodes with the ONIOM(r2SCAN-3c:GFN2-xTB) scheme compared to that including three nodes. This can be explained due to faster convergence of the geometries resulting from the consistent DFT treatment of the zirconium-centered moieties.
Another possible use case of our implementation is the modeling of explicit solvation. If an implicit solvation model becomes insufficient to describe the solvent effect, it is necessary to include explicit solvent molecules which can increase the computational costs tremendously. To still include a sufficient amount of explicit solvents to describe the system accurately, the ONIOM scheme can be used to employ a DFT or WFT method for the solute and a GFN method for the solvent. In this case, no covalent bonds are broken, which renders the GFN methods specifically suitable as they are designed to reproduce non-covalent interactions accurately.
As an example, the Bis(cyclopentadienyl)silylniobium complex [Cp2 Nb(H)2 (SiCliPr2)]78 (Fig. 3A) was solvated with 20 DMSO molecules using the quantum-cluster-growth (QCG) algorithm.30,79 This transition metal complex has previously proven to be a specifically difficult case for common semi-empirical methods in the TMG145 benchmark study.19 The solvated complex was optimized at three different levels, TPSS-D4(CPCM)/def2-SVP (Fig. 3C), GFN2-xTB(ALPB), and ONIOM(TPSS-D4/def2-SVP:GFN2-xTB(ALPB)). In the ONIOM scheme, the bis(cyclopentadienyl)silylniobium complex was chosen as the inner region. A comparison of the cutout solute geometry optimized within the ONIOM scheme and pure GFN2-xTB is shown in Fig. 3B. Optimizing the explicitly solvated complex with TPSS-D4/def2-SVP retains the major structure motif of the Bis(cyclopentadienyl)silylniobium complex with respect to the gas-phase TPSSh80-D3(BJ)-ATM reference structure.19 Using GFN2-xTB for optimizing the whole cluster yields a qualitatively wrong, distorted structure with an RMSD of 0.26 Å of the cutout solute in respect to the fully DFT optimized geometry. Thereby, a haptotropic shift of the cyclopentadienyl ligands occurs and the symmetric coordination of the terminal hydride ligands is distorted. Accordingly, GFN2-xTB is, similar to other semi-empirical and Force Field methods,19 not able to describe this complex system accurately. Applying the ONIOM(TPSS-D4/def2-SVP:GFN2-xTB) scheme instead yields the correct structure in excellent agreement with the full DFT optimization (RMSD of 0.01 Å) at a fraction of computational costs (59 min vs. 48 h 58 min). Despite the good performance of the ONIOM scheme demonstrated for this example, electrostatic embedding may enhance the results for cases with highly polar environments showing strong charge–dependent interactions. Nevertheless, the appropriate choice of the ONIOM boundary may reduce errors in this respect.
As a first example, we discuss the Rh-functionalized metal–organic cuboctahedron illustrated in Fig. 4. In the following, it is referred to as DALTES81 in correspondence with its reference code from the Cambridge Structural Database (CSD).82
Herein, we consider the cyanosilylation reaction facilitated by the dirhodium paddle-wheel nodes (Fig. 5A). The proposed mechanism depicted in Fig. 5B was adapted from Zhang et al.83 The first step is the coordination of an aldehyde substrate to the open Rh sites of DALTES, which is followed by the nucleophilic addition of the isomerized trimethylsilyl cyanide (iso-TMSCN) to the activated carbonyl compound. Afterwards, the trimethylsilyl group isomerizes into the product which is the rate-determining step of the catalytic reaction. Finally, the produced cyanohydrin derivative is cleaved from the metal cluster of the catalysts. To validate the ONIOM approach, a full quantum mechanical reference energy profile was calculated. Due to the size of the investigated system (404 atoms), the choice of the reference method was again limited to the (meta-)GGA/DZ level. Accordingly, the TPSS-D4/def2-SVP level with COSMO implicit solvation was chosen as a reference and QM component in the ONIOM approach. The energy profile of the catalytic cycle was then recomputed at the GFN2-xTB and ONIOM(TPSS(COSMO):GFN2-xTB(ddCOSMO)) level (Fig. 5). For the ONIOM scheme, only the respective dimetal node involved in the cyanosilylation and the corresponding reactants were included in the inner region.
The semi-empirical GFN2-xTB energy profile qualitatively agrees with the full-DFT results but the relative energies differ by up to 35 kcal mol−1. Applying the ONIOM approach instead yields very good agreement with the TPSS-D4/def2-SVP reference energies, only varying by 1–3 kcal mol−1. Accordingly, the ONIOM scheme can be used to effectively converge the results by combining GFN2-xTB with any suitable DFT method of choice. Using ONIOM reduces the average computational wall time of the energy evaluations from 30 minutes to only 40 seconds. This allows choosing a significantly more accurate method in the inner region that would be unfeasible in a full DFT approach. Similar results are obtained for other methods combinations such as ONIOM(r2SCAN-3c:GFN2-xTB) and ONIOM(ωB97X-3c84:GFN2-xTB) (see ESI†).
An even more challenging system is the recently synthesized porphyrin-based macromolecular spoked-wheel complex by Majewski et al.,85 illustrated in Fig. 6. Its size of 870 atoms renders most conventional DFT methods unfeasible, but the bond formations at the rim of the spoked-wheel complex can be computed employing DFT within the ONIOM approach. The rim of the spoked-wheel complex consists of 18 5,15-linked porphyrin units coordinated by Zn(II) or Ni(II) ions. By the addition of bis(trifluoroacteoxy)iodobenzene (PIFA), the meso–meso coupled Zn–porphyrin rings can be fused further forming a nanobelt involving three-fold single bond coupled Zn–porphyrin units. The corresponding reaction energy of this coupling reaction was computed including the Zinc-functionalized porphyrin rings in the inner region with 432 atoms. The obtained reaction energies are shown in Table 1. The reference method (ONIOM(TPSS:GFN2-xTB(GBSA))) yields an overall reaction energy of −44.9 kcal mol−1. GFN2-xTB and GFN1-xTB yield reasonable reaction energies of −37.6 and −34.9 kcal mol−1, respectively, and show the qualitatively right trend. Using only GFN-FF yields a much too high, qualitatively wrong interaction energy. Combing it with GFNn-xTB within the ONIOM scheme can reduce this error drastically and gives a qualitatively correct behavior. These results show that it is possible to get good and fast results with the ONIOM(GFNn-xTB:GFN-FF) scheme and that, if more accurate results are required for such large systems, ONIOM(DFT:GFN) can be a valuable tool as full DFT calculations are yet too costly.
Fig. 6 The molecular structure of the 5,15-linked porphyrin nanoring and the fusion reaction between its Zn-functionalized porphyrin units. |
Methods | ΔE/kcal mol−1 | SCF time |
---|---|---|
ONIOM(TPSS†:GFN2-xTB) | −44.9 | 25 m |
GFN2-xTB | −34.9 | 36 s |
ONIOM(GFN2-xTB:GFN-FF) | −24.4 | 4 s |
ONIOM(TPSS†:GFN1-xTB) | −47.0 | 24 m |
GFN1-xTB | −37.6 | 39 s |
ONIOM(GFN1-xTB:GFN-FF) | −26.6 | 4 s |
ONIOM(TPSS†:GFN-FF) | −35.3 | 15 m |
GFN-FF | 172.6 | 0.4 s |
It was shown that the ONIOM approach can be utilized to clearly improve on the already reasonable results of the GFN methods for challenging systems containing many transition metal atoms. By matching the strengths of the GFN methods and DFT, even critical cases that require a highly accurate quantum mechanical description can be treated. Accordingly, we demonstrated that by using the ONIOM framework in xtb, DFT quality reaction energies and molecular structures can be obtained at a fraction of the computational costs of conventional DFT. This allows also the use of computationally more expensive and more accurate DFT methods. The efficient, robust, and easy-to-use implementation of the ONIOM scheme into xtb opens up new possibilities with regard to the treatment of very large systems containing a variety of elements across the periodic table.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3cp02178e |
This journal is © the Owner Societies 2023 |