Efficient Computation of Geometries for Gold Complexes.

Computationally obtaining structural parameters along a reaction coordinate is commonly performed with Kohn-Sham density functional theory which generally provides a good balance between speed and accuracy. However, CPU times still range from inconvenient to prohibitive, depending on the size of the system under study. Herein, the tight binding GFN2-xTB method [C. Bannwarth, S. Ehlert, S. Grimme, J. Chem. Theory Comput. 2019, 15, 1652] is investigated as an alternative to produce reasonable geometries along a reaction path, that is, reactant, product and transition state structures for a series of transformations involving gold complexes. A small mean error (1 kcal/mol) was found, with respect to an efficient composite hybrid-GGA exchange-correlation functional (PBEh-3c) paired with a double-ζ basis set, which is 2-3 orders of magnitude slower. The outlined protocol may serve as a rapid tool to probe the viability of proposed mechanistic pathways in the field of gold catalysis.

Introduction

Gold catalysis is a rapidly expanding field (for representative review articles see Ref. [1]) in which computational modelling is often used to aid the interpretation of experimental results. Typically, one would like to obtain a set of chemically‐meaningful structures along a reaction path, to understand geometric and energetic changes. To this end, Kohn‐Sham density functional theory (KS‐DFT or simply DFT) is by far the most used method due to its favourable cost‐accuracy ratio and applicability to a wide‐range of organometallic systems. Particular exchange‐correlation (XC) functional and basis set (BS) combinations (e. g. B3LYP/6‐31G*) have found prolific use where scalar relativistic effects are treated through the use of a pseudo potential (ECP). More recently, the ‘3c’ composite methods from the Grimme group have been put forward as a replacement as they overcome certain drawbacks. These methods reparametrize successful DFT functionals and introduce (empirical) contributions that efficiently correct for inherent limitations of the functional and small BS, such as corrections for the dispersion interaction and the (intra‐ and inter‐ molecular) basis set superposition error. These methods greatly facilitate computation of large systems. However, some calculations still require days of CPU time. Recently, an extremely fast semi‐empirical tight binding method was proposed, coined GFN2‐xTB and referred to hereafter as xTB. It seeks to approximate the KS‐DFT energy in terms of fluctuations in the electron density and accounts for multipole electrostatics (up to 2nd order), (self‐consistent) dispersion and static electron‐electron correlation via Fermi smearing. xTB was parameterized for near‐equilibrium geometries but has been shown to reproduce other equilibrium properties, such as dipole moments, better than all other semi‐empirical approaches, indicating the effectiveness of its parameterization strategy and physically‐motivated design. If very large systems are to be treated the generally‐applicable and robust GFN‐FF force field is available. More recently, it was shown that xTB also produces good geometries for transition metal complexes. The article also alludes to applications of the xTB method for the study of reaction mechanisms which has previously been applied to, e. g. Zr complexes and lanthanoid complexes. While in the aforementioned studies, xTB is used to explore the reaction path, the transition states (TS) were still reoptimized with KS‐DFT, which has recently been further explored in a semi‐automated approach for organometallic complexes. With the field of gold chemistry and especially the area of catalysis in mind, a key question that arises from these reports is whether xTB can produce reasonable structures along a reaction path, including TS, for structures involving gold complexes. If so, it would be an effective tool to explore reactivity, especially for experimentalists who might only be interested in probing if a given scenario/path is viable or not, with the aim of supporting a mechanistic proposal rather than performing an extended computational study which would be much more time demanding and might be beyond the scope of a given study. In the present article we will demonstrate, for a series of transformations of gold complexes, that the xTB method provides meaningful structures that give insight into reaction mechanisms at very low computational cost. Some technical aspects are highlighted that make the use of the xTB method straightforward.

We begin by investigating the hydroamination of ethyne, shown in Figure 1, as a small model reaction following the procedure of Ciancaleoni et al. We compare xTB with a variety of the ‘3c’ composite methods (HF‐3c, PBEh‐3c and B97‐3c) and popular XC functional/basis set combinations combined with different dispersion correction schemes (B3LYP ‐D3(BJ)[ , ]/def2‐SVP, B3LYP ‐D3(BJ)[ , ]/def2‐TZVPP and B3LYP ‐D4 /def2‐TZVPP ) for this purpose. These structures, and the additional tri‐coordinated species (LAu(NH3)(C2H2)), with geometries optimized at each of the stated methods define data set [A]. This allowed for an appropriate choice of DFT functional to compare xTB against in the second data set, [B], which consists of a representative set of stoichiometric gold reactions (shown in Figure 2), which generally include more realistically‐sized systems. The geometries in data set [B] are optimized with PBEh‐3c and xTB and include 1) the intramolecular oxidative addition of a naphthyl iodide moiety; 2) a Friedel‐Crafts reaction with Au(I) & Au(III); 3) cyclopropanation involving a Au(I) carbene; 4) a 1,6‐enyne cycloisomerization with Au(I) and 5) the gold‐mediated Cope rearrangement of a 1,5‐diene.

Figure 1

The stationary points along the reaction path of the hydroamination of data set [A], L=Cl,PH3 +.

Figure 2

The small but diverse set of reactions selected for the main data set, [B], with larger Au complexes. [Au]=AuCl for reaction (2a), AuCl3 for (2b), AuPH3 + for (4) and AuPMe3 for (5).

We probe the similarity of the xTB and PBEh‐3c structures by evaluating their energies at the domain‐based pair natural orbital coupled cluster method, DLPNO‐CCSD(T), to provide a more accurate picture. We additionally use this opportunity to probe the accuracy of the double hybrid GGA XC functional B2‐PLYP‐D3(BJ) for the computation of accurate energetics of homogeneous gold systems, which has previously been shown to perform well for computational studies of gold chemistry.[ , , ] We also probe the Minnesota XC functional PW6B95‐D3(BJ), which also performs well across all systems studied. Overall, the combination of xTB with the relatively inexpensive PW6B95‐D3(BJ) functional (PW6B95‐D3(BJ)//xTB) provides an excellent balance between rapid evaluation of geometries and accurate energetics.

Computational Details

General Comments

All calculations were performed with ORCA 4.2.1. xTB refers to GFN2‐xTB v6.2.39 and was used interfaced with ORCA. For xTB calculations the default SCF convergence criterion of 10−6 au was used. All other calculations use VeryTightSCF settings. The nature of all the stationary points was confirmed via frequency analyses and all TS structures were found to connect reactants and products via intrinsic reaction coordinate (IRC) calculations.

Hydroamination, Data Set [A]

All calculations were performed in the gas phase. The reference geometries were taken from Ciancaleoni et al. which had been computed at the density‐fitting local coupled cluster (DF‐LCCSD(T)/TZVPP) level of theory. In the present study these geometries were reoptimized with the selected methods shown below (Table 1).

Table 1

Computational methods evaluated for data set [A].

Entry	Method	Basis Set
1	xTB [28]	–
2	HF‐3c [30]	MINIX [30]
3	PBEh‐3c [31]	def2‐mSVP [31]
4	B97‐3c [5a]	mTZVP [5a]
5	B3LYP [18] ‐D3(BJ)[ 7 , 19 ]	def2‐SVP [20]
6	B3LYP [18] ‐D3(BJ)[ 7 , 19 ]	def2‐TZVPP [20]

–

HF‐3c

MINIX

PBEh‐3c

def2‐mSVP

B97‐3c

mTZVP

B3LYP ‐D3(BJ)[ , ]

def2‐SVP

B3LYP ‐D3(BJ)[ , ]

def2‐TZVPP

For brevity, B3LYP‐D3/S and B3LYP‐D3/L are used throughout as labels to refer to B3LYP‐D3(BJ)/def2‐SVP and B3LYP‐D3(BJ)/def2‐TZVPP, respectively. Analytical second derivatives were used for KS‐DFT calculations and HF‐3c. For xTB calculations numerical second derivatives were employed. For all geometries, electronic energies were evaluated using the DLPNO‐CCSD(T) level of theory in combination with the def2‐TZVPP basis set and the def2‐TZVPP/C fitting basis set. The TightPNO setting was used throughout. Relative electronic energies were compared, as outlined in ref. [17]. For KS‐DFT, HF‐3c and DLPNO‐CCSD(T) calculations scalar relativistic effects were accounted for by replacing the 60 inner electrons of Au by the appropriate pseudo potential (ECP) for the specified basis set. DLPNO‐CCSD(T)/def2‐TZVPP single point energies obtained for the previously reported DF‐LCCSD(T)/TZVPP geometries from ref. [17] are used as reference values.

Larger Au complexes, Data Set [B]

The performance of xTB and DFT (PBEh‐3c) for structural parameters was further compared in data set [B]. For geometry optimizations with xTB the Generalized Born implicit solvation model with Surface Area contributions (GBSA, as implemented in the xTB code)[ , ] with the default (230 point Lebedev) grid, was added, unless stated otherwise (see discussion below). To maintain consistency, the conductor‐like polarizable continuum model (cPCM)[ , ] was used for the corresponding PBEh‐3c optimizations, using a modified (Gaussian van der Waals) cavity as implemented in ORCA. Water was specified in both cases. Numerical second derivatives were used for both DFT and xTB calculations. For reaction 1 shown in Figure 2 we also computed nudged elastic band (NEB) paths (as implemented in the ORCA code) with various methods. Initial geometries for the five reactions were taken from refs. [21-24].

Final Single Point Energies

For comparison between the xTB and DFT geometries we again computed electronic energies at the DLPNO‐CCSD(T)/def2‐TZVPP level of theory as outlined above. In addition, we used this opportunity to probe the performance of the PW6B95 and B2‐PLYP density functional approximations in combination with the D3(BJ) dispersion correction with Becke‐Johnson damping (PW6B95‐D3(BJ) and B2‐PLYP‐D3(BJ), respectively), employing the def2‐TZVPP basis set. For these calculations Grid6 and GridX6 were used. Calculations were accelerated by the use of the RIJCOSX approximation using Weigend's universal fitting basis set. For B2‐PLYP‐D3(BJ), the RI approximations was also employed for the MP2 step, using the def2‐TZVPP/C fitting basis set. Energies were referenced to the fully separated reactants for the bimolecular cyclopropanation (3).

Results and Discussion

Hydroamination

The reliability of the calculated structural parameters of the various methods was judged by evaluating the energy of each structure with local coupled cluster DLPNO‐CCSD(T)/def2‐TZVPP. The errors were calculated as deviations with respect to the previously reported reference geometries. The statistics of these results can be found in Figure 3. xTB represents a clear improvement over HF‐3c, despite being ∼2 orders of magnitude faster. The other tested DFT functionals outperform xTB to a similar extent, with the exception of B3LYP‐D3/L which produces the smallest errors (and has the largest basis set). The performance of the self‐consistent D4 model was also probed and found to give quite comparable results to B3LYP‐D3/L, indicating that the D3 dispersion correcting already performs very well. These results can be found in the Table S6. Despite this, PBEh‐3c was selected for the more extensive comparison against xTB (data set [B]) as it is significantly more efficient than B3LYP‐D3/L. This choice is further motivated by the common use of a hybrid GGA/double‐ζ BS combination in computational studies of Au chemistry (for representative examples see refs. [1b, 1e, 1 f, 23–24, 42]) and PBEh‐3c has previously found successful applications in this field.

Figure 3

Mean Error (ϵ), Standard Deviation (σ) and Maximum Error (M) of the DLPNO‐CCSD(T)/def2‐TZVPP energy evaluated at the geometry of the tested methods, with respect to the previously reported DF‐LCCSD(T) structures in data set [A].

Larger Au Complexes

Inclusion of Implicit Solvation

The oxidative addition (reaction 1 shown in Figure 2) was the first reaction of the larger set [B] to be investigated. Although xTB was able to locate the TS, its IRC led to the intermediate shown in Figure 4b. To investigate this further, a NEB path – which approximates the minimum energy path connecting two structures by iteratively optimizing the linear interpolation between them – was obtained from the optimized reactant (Figure 4a) and product geometries (Figure 4c). At the DFT level, a single TS is identified, leading to a downhill process, b→c, which – using xTB – exhibits an additional barrier. This makes the overall process (a→c) appear to be two‐step rather than concerted. We do note here, that intermediate b obtained from xTB may indeed be optimized to a minimum in the gas phase at the DFT level (see ESI). Nevertheless, neither the NEB nor the IRC result in this structure. The inclusion of the implicit solvation model (GBSA) for the xTB method remedies this qualitative difference and also leads to a NEB path that similarly displays a single TS. When the same DFT level of theory was used to evaluate the energy of these geometries (DFT/cPCM//xTB/GBSA), the reaction curves agree to a remarkable extent; the (electronic) barriers are within 3 kcal/mol and the TS lies at almost the same point along the path. Encouraged by this result, xTB/GBSA was used for all further calculations. An implicit solvation model (cPCM) was also used in the accompanying PBEh‐3c calculations, to allow for an appropriate comparison.

Figure 4

The NEB reaction profile of the oxidative addition (1), with xTB and DFT (PBEh‐3c), with and without their implicit solvation models (GBSA and cPCM, respectively) and single point DFT energies evaluated at the xTB geometries (DFT/cPCM//xTB/GBSA). *Minimum in the gas phase. See text for more details.

Overall Performance

Before comparing some key interatomic distances, the DLPNO‐CCSD(T) energy is used as a single parameter to judge the geometric performance of xTB vs. the selected double‐ζ DFT composite method, PBEh‐3c (Figure 5a). The xTB result is, in general, within a few kcal/mol and the ‘shape’ of the reaction is well reproduced. The only exception to this is the product ion pair complex (v2) of the cyclopropanation (reaction 3). This is because v2 DFT has L−Au bound to the newly formed cyclopropane moiety (Figure 5c), while v2 xTB has gold coordinating to the triflate anion (L−Au‐OTf, Figure 5b), as in the disassociated products. In other words, this ion pair is being overlooked by xTB – perhaps indicating an over‐stabilizing of the Au−O bonds. This was the only qualitative divergence found in the entire data set, and the interpretation of the end point of an IRC calculation as representative of the (often non‐isolatable) ion pair complex may also vary within standard DFT approaches. The unarguably crucial stationary points of this reaction exhibit only 0.4 and 1.1 kcal/mol deviations between xTB and DFT for the transition state and (separated) product geometries, respectively. The error statistics of this comparison can be found in Figure 6a and showcases xTB's good performance, with a low mean error (ϵ=−0.4 kcal/mol) and standard deviation (σ=2.0 kcal/mol). The Maximum Error (M=−8.1 kcal/mol) is due to the discussed difference in the cyclopropanated v2 structures.

Figure 5

a) The DLPNO‐CCSD(T) energies of the larger Au set [B] from Figure 2 evaluated at the xTB or DFT geometry, demonstrating their broad agreement. Transition state (TS), product (P), intermediate (int) and vdW complex (v) energies are all referenced to the reactant (R). The ion pairing product complex of the cyclopropanation, optimized with b) xTB and c) DFT, showing their qualitative differences.

Figure 6

Mean Error (ϵ), Standard Deviation (σ) and Maximum Error (M) of a) DLPNO‐CCSD(T)//xTB vs. DLPNO‐CCSD(T)//PBEh‐3c b) PW6B95‐D3(BJ)//xTB vs. DLPNO‐CCSD(T)//xTB c) B2‐PLYP‐D3(BJ)//xTB vs. DLPNO‐CCSD(T)//xTB & d) PW6B95‐D3(BJ)//xTB vs. DLPNO‐CCSD(T)//PBEh‐3c for data set [B].

Structural Performance

The performance of xTB (vs. PBEh‐3c) for obtaining structural parameters is examined explicitly in Table 2. The mean error of the transition states (ϵ=−2.3 pm) is marginally lower than that of the reactants and products (−2.9 & −3.0 pm, respectively), although there is a larger spread – seen in the greater standard deviation (σ=6.3 pm vs. 3.1 & 3.2 pm) and maximum error (M=−20.5 pm vs. −10.3 & −11.1 pm).

Table 2

Metal‐to‐ligand interatomic distances and error statistics (in pm) of the xTB optimized structures in data set [B],[a] DFT refers to the composite PBEh‐3c method.

Reaction	atom #1‐atom #2[b]		xTB				Deviation (w.r.t DFT)
			R	TS	P		R	TS	P
1	Au−P40 Au−I42 Au−I16 Au−C2 Au−Cl24 Au−C20 Au−C21 Au−Cl24 Au−Cl25 Au−Cl26 Au−C20 Au−C21 Au−C2 Au−C3 Au−P11 Au−P5 Au−C1 Au−O20 Au−P2 Au−O16		226.7	247.6	226.0		−3.6	8.5	−5.3
			262.6	268.6		−1.6	−3.5	−3.3
			279.6	263.9		−1.5	−20.5	−4.7
			209.0	205.8		−0.2	−5.1	−1.6
2a			230.4	236.5		−10.3	−0.9	1.5
			205.0	202.5		−1.4	−0.7	0.3
			257.0	293.1		−4.9	−2.6	−0.3
2b			230.7	227.2		−2.4	0.9	−3.5
			230.6	240.5		−1.4	0.2	0.3
			231.0	227.3		−2.0	1.0	−3.2
			221.3	198.8		0.3	6.3	−2.5
			237.2	288.4		−0.3	−7.7	−0.8
4			283.0	305.1		1.8	4.5	−2.0
			208.9	202.0		−5.1	−0.4	−0.8
			228.5	229.0		−1.7	−2.7	−5.3
5			TS−I	int		R	TS−I	int
			285.7	306.6		−5.7	−5.6	−3.8
			215.7	212.2		−9.1	−0.2	0.0
			227.8	228.6		−2.9	−3.1	−4.2
			TS−II	P			TS−II	P
			216.0	230.4			0.0	−11.1
			264.0	218.8			−14.2	−10.2
			231.6	230.5			−3.2	−1.9
Mean error (ϵ)							−2.9	−2.3	−3.0
Standard deviation (σ)							3.1	6.3	3.2
Maximum error (M)							−10.3	−20.5	−11.1

[a] The cyclopropanation (3) has been omitted due to the previously discussed qualitative differences, please see Table S8 for further information. [b] The subscripts in the atom labels correspond to the atom numbers in the xyz coordinate tables, which can be found in the Supporting Information.

Energetics

While DLPNO‐CCSD(T) energies are feasible, they are far from convenient due to long CPU times and high memory usage. In practice, one would often like to acquire accurate energies with a less costly approach. Here, a small and non‐exhaustive comparison is made between two such approaches: a DFT functional well‐parameterized for thermodynamics (PW6B95‐D3(BJ)) and a double hybrid (B2‐PLYP‐D3(BJ)) functional. Their deviations from the DLPNO‐CCSD(T) reference for the data set [B] are shown in Figures 6b and 6c, respectively. Both methods perform similarly well, with PW6B95‐D3(BJ)’s maximum error (M=−4.1 kcal/mol) a little lower than B2‐PLYP‐D3(BJ)’s (M=−5.4 kcal/mol). Notably, the maximum error corresponds in both cases to the transition state of reaction 1. Comparing PW6B95‐D3(BJ) and PW6B95‐D4 (i. e. replacing the D3 dispersion model with self‐consistent D4 model) slightly reduced the standard deviation (σ) at no additional computational cost. These results, along with those of B3LYP‐D3(BJ) can be found in the SI.

Since Figure 5 shows the good performance of xTB (vs. DFT) for structural parameters and Figure 6b shows the good performance of PW6B95‐D3(BJ) (vs. DLPNO‐CCSD(T)) for obtaining electronic energies, one may naturally ask if they might work together in a hybrid procedure of geometry optimization with xTB and single point energy evaluation with PW6B95‐D3(BJ) (PW6B95‐D3(BJ)//xTB). The error statistics of such an approach, with reference to the DLPNO‐CCSD(T) energies evaluated at the DFT geometries (DLPNO‐CCSD(T)//PBEh‐3c), can be found in Figure 6d and demonstrate that the relative energies for all the species in Figure 2 (and the TS connecting them) could have been obtained, on average, within 1.1 kcal/mol at a mere fraction of the computational effort. The large maximum error (M=−8.1 kcal/mol) is due to the previously discussed complications of the cyclopropanation (Figure 5), which might be judged to be of limited relevance for exploring reaction paths.

Conclusions

The extremely efficient tight‐binding method xTB was found to be capable of producing very reasonable geometries of five diverse reactions of Au complexes, in comparison to the well‐balanced KS‐DFT composite method PBEh‐3c. This holds for both minima (as xTB was designed to do) and transition state geometries (a far‐from equilibrium property that it was not originally intended for). PW6B95‐D3(BJ) and B2‐PLYP‐D3(BJ) reproduced the DLPNO‐CCSD(T) energies with similar accuracy. This is in line with the results of Dohm et al., and demonstrate that well‐parameterized (meta‐)hybrid and double‐hybrid density functionals can be effectively used to evaluate energies. The inclusion of an implicit solvation model for xTB optimizations was found to aid in reproducing DFT‐like potential energy surfaces. Some care should be taken when examining ion pair/loosely bound complexes of bimolecular reactions as the only tested structure of this kind was found to be qualitatively different to DFT but aside from these caveats, PW6B95‐D3(BJ)//xTB has been shown to be an effective tool at rapidly exploring the reactive space of Au complexes. Similarly, it has recently been shown that the conformational space of transition metal complexes may be efficiently probed with a similar computational procedure employing geometries from xTB, rendering this approach a comprehensive tool for the exploration of the chemical space of transition metal complexes. With the often more qualitative perspective of many experimental studies, asking simply if a proposed mechanistic pathway is feasible or not, we feel that xTB provides an excellent balance between speed and accuracy, especially when used in the proposed PW6B95‐D3(BJ)//xTB scheme.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supplementary

Click here for additional data file.

Supplementary

Click here for additional data file.