Reaction mechanism between Cu(II)-enolate complex and O2 as a test case for methodology used in DFT computational studies

The reaction mechanism of an intricate oxidation reaction of chlorodiketonate ligand of mononuclear Cu(II) complex was studied computationally employing five different models that differ in: a) basis set, b) the way that solvent corrections are included, and c) DFT functional. Qualitative and quantitative comparison of structures and enthalpy reaction profiles enabled us to assess how sensitive they are to the changes in computational methodology. Graphical abstract Comparison of enthalpy reaction profiles and molecular structures demonstrate how the qualitative picture on Cu(II)-catalyzed reaction changes upon variation of computational methodology Comparison of enthalpy reaction profiles and molecular structures demonstrate how the qualitative picture on Cu(II)-catalyzed reaction changes upon variation of computational methodology


Introduction
Using quantum chemistry methods to model (bio)inorganic reaction mechanisms for systems involving first row transition metals (TM) is a challenging task. First, such systems involve coordination bonds that are difficult to describe quantitatively with robust DFT methods. Second, electronic configurations with open 3d shells often result in several spin states that are close in energy and a so-called multi state reactivity may emerge [1]. Finally, the metal ion(s) may be redox active, and hence may promote radical mechanisms aside from the heterolytic processes. Even though recent years have witnessed major breakthroughs in the development of very promising post-HF methods, such as DMRG-CASPT2 or local CCSD(T) methods [2,3], their computational cost still remains relatively high when they are used with adequate one-electron basis sets. Therefore, it is foreseen that much less demanding DFT methods will continue to be widely used for geometry optimizations and energy calculations, whereas these more advanced post-HF methods will be employed to compute accurate electronic energies for (selected) structures optimized with DFT.
Embarking on a DFT computational study on a reaction mechanism for a TM system necessitates taking several decisions on the computational methodology, including which density functional and which basis set shall be used for geometry optimization and which for energies or weather a solvent model shall be applied during optimization or rather it suffices to use it only as an energy correction for structures optimized in vacuum.
To assess to what extent such particular decisions may affect conclusions one draws from computational results, we did a small yet systematic case study for a synthetic mononuclear copper(II)-chloroenolate complex, for which congruent experimental and computational results are already available [4,5]. This system has several particularly attractive features concerning testing variations of computational methodology. First, it is a mononuclear system of size that is not prohibitively large, so both double and triple-zeta one electron basis sets may be tested for geometry optimization and frequency calculations. Second, the reaction between this complex and dioxygen features several plausible yet different reaction types, i.e., reduction of dioxygen to peroxo species, nucleophilic substitution as a prelude accompanying oxidation of chlorine, aliphatic or aromatic carbons. Finally, experimental results are available and they can be used as a reference point when assessing various computational methodologies.
The model system is a cationic chelate complex of Cu(II) with 6-Ph 2 TPA and PhC(O)CClC(O)Ph chlorodiketonate as ligands (S), or its derivative with one chloride anion bound in the first coordination shell of copper (S cl ; see Scheme 1). Results of previous computational investigations on the mechanism of the reaction between these complexes and dioxygen have already been published, yet for the sake of completeness the key findings are summarized here. Thus, the X-ray crystal structure of a salt with Cu(II) complex S as a cation and ClO 4 − as an anion is know [4]. In acetonitrile, the reaction between S and O 2 proceeds with an initial slow lag phase that can be eliminated by addition of catalytic amounts of chloride salts.
Crystal structure of S Cl is not available, yet UV-vis and EPR spectroscopic data together with the DFT optimized structure suggest that chloride anion coordinates Cu in a complex with distorted square pyramidal geometry. The composition of the final reaction mixture is the same as the one obtained in an independent reaction between NaOCl, d i p h e n y l p r o p a n e t r i o n e ( t r i k e t o n e ) , a n d [ ( 6 -Ph 2 TPA)Cu(CH 3 CN)](ClO 4 ) 2 , which strongly indicates that triketone and hypochlorite are formed as reactive intermediates in the initial stage of the reaction [5].
From the DFT results obtained for the parent compound S, it follows that in the first step of the reaction one of the oxygen atoms of the chloroenolate detaches from Cu(II) making a place for dioxygen; O 2 binds first to copper and then to the central carbon (C2) of the enolate. Thus, formed peroxobridged intermediate P subsequently undergoes a remarkable and concerted, yet asynchronous transformation, whereby first, the oxygen atom proximal to copper (O1) makes an SN2 type nucleophilic attack on C2 expelling Cl − and second, this liberated chloride anion attacks a strained 3-membered CO 2 ring cleaving the O-O bond heterolitically and forming an adduct between hypochlorite and triketone (T). The two barriers connected with TS1 and TS2 are substantial and comparable in size. The reaction between S Cl and O 2 proceeds in a similar way, yet the presence of Cl ligand changes the mechanism of the second step, as in TS2 Cl the oxygen atom O1 is attacked by Cl2, which elicits O-O cleavage and dissociation of Cl1 as Cl − . Both barriers, connected with TS1 Cl and TS2 Cl , are substantially lower than for the parent compound lacking the chloride ligand, which successfully explains the catalytic effect exerted by the chloride ion.
Two alternative pathways were found for the decay of the peroxo intermediate P Cl (Scheme 2).
First, the O-O bond may be cleaved homolytically via TS2 Cl_homo . Second, instead of oxidizing Cl2, one of the phenyl appendages of the 6-Ph 2 TPA ligand may be oxidized to Scheme 1 Mechanisms for the initial stage of oxidation of S or S Cl formulated on the basis of a joined experimental and computational study [5] arene oxide (T cl_Ph ). In the original study, both of these processes involved barriers higher than that connected with TS2 Cl ; however, the latter reaction was also found for the parent compound P, and in this case TS2 Ph was lower in energy than TS2.
In this study we re-examined the above summarized mechanisms using several variations of the computational methodology originally applied to this system. From comparison of the results we were able to draw several conclusions that, in our opinion, will be useful when one considers a similar project.

Computational details
Three one-electron basis sets were used in this study. BS1 of double-ζ quality consisted of LANL2DZ basis and an effective core potential for Cu and 6-31G basis for other elements. BS2 of triple-ζ quality for the key atoms combining lacv3p + basis for Cu, 6-311G* for Cl, N, O, and three central carbon atoms of chlorodiketonate and 6-31G for other carbon and all hydrogen atoms. BS1 or BS2 were used for geometry optimization, frequency calculations, and computation of solvent effects. The final electronic energies were computed with BS3, which combined the lacv3p + basis for copper with cc-pVTZ(−f) basis for other elements.
Three hybrid density functionals were used for computations: B3LYP [6], TPSSh [7], and ωB97XD [8]. All three were combined with the D2 dispersion correction [9]; in the case of ωB97XD the D2 correction is included in the functional formula. For B3LYP and TPSSh, the final electronic energy, which was computed with the BS3 basis set, was combined with the D3 dispersion correction obtained with the DFT-D3 program [10]. Solvent (acetonitrile) effects were computed with the polarizable continuum model using the integral equation formalism (IEFPCM). Computations were done with the spinunrestricted formalism using Gaussian 09 [11].
For species with antiferromagnetic coupling between unpaired electrons, the electronic energy was corrected with the procedure proposed by Yamaguchi and co-workers [12].
All reported values combine: electronic energy computed with BS3, Yamaguchi correction, if applicable, D3 correction for the appropriate functional or D2 built-in ωB97XD, solvent effects computed either as single point corrections or used also for geometry optimization, and zero point and thermal corrections to enthalpy computed with the harmonic approximation at the same level of theory as used for geometry optimization.
The results taken from the previous study are used as a reference and are labeled as B3LYP-D3/BS3//B3LYP-D2/ BS1 [5]. In previous work, although the TPSSh functional was used, it was only to compute electronic energies for geometries of stationary points that were optimized at the B3LYP-D2/BS1 level; such energy profiles were very similar to those obtained with B3LYP only. In the present work, the TPSSh functional is used both to optimize geometries and compute final energies; hence, the results and conclusions may be different than from previous work.

Results and discussion
The representative enthalpy reaction profile, computed with the B3LYP-D3/BS3//B3LYP-D2/BS2 model, is presented in Fig. 1, and profiles computed with other models can be found in Supporting information. Relative enthalpy values obtained for all models are gathered in Table 1, whereas Table 2 presents key interatomic distances for stationary points along the reaction coordinates.
Concerning the qualitative aspects, all of the enthalpy profiles computed with different models reproduce the catalytic effect introduced by the chloride ligand, as the barriers connected with TS1 Cl and TS2 Cl are markedly lower than those Scheme 2 Alternative reactions of P Cl [5] corresponding to TS1 and TS2. Similarly, for all models the first elementary step, i.e., S → P, is more exothermic when Cl is bound in the first coordination shell of Cu(II), whereas T Cl_homo is the least stable of all products. Relative stability of other products is somewhat less preserved among the models. However, some trends may be noticed, as for almost all models T Cl_Ph is the lowest enthalpy product; the exception is the TPSSh-D3/BS3//TPSSh-D2/BS1 model, for which the most stable product is T Ph , with T Cl_Ph being the second most stable. Relative barrier heights for various decay channels of the peroxo species P indicate that for the parent system T Ph would be the most easily accessible product for all models except the ωB97XD/BS3//ωB97XD/BS1 model, for which T is predicted to be produced with the lowest barrier. Analogous analysis for the channels starting at P Cl shows that T Cl is the kinetically preferred product for all models except TPSSh-D3/BS3//TPSSh-D2/BS1, for which T Cl_homo is accessible with the lowest barrier. The latter prediction when confronted with the experimental observation of identical reaction outcomes from reactions starting with S/O 2 or triketone/ClO − may suggest that the TPSSh-D3/BS3// TPSSh-D2/BS1 model probably incorrectly favors the homolytic O-O bond cleavage over the heterolytic process leading to ClO − .
Comparing qualitative features of the enthalpy profiles computed with the three different hybrid functionals enables one to conclude that the choice of the functional affects the computed landscape more than using a bigger basis set or solvent model for geometry optimization. The most important differences can be noticed for the decay of the peroxo intermediate. One difference has already been mentioned above, i.e., the TPSSh model predicts the lowest barrier for homolytic O-O cleavage starting from P Cl . Another qualitative difference is that for the parent system the lowest barrier for decay    of the P species varies with the functional, i.e., in the case of B3LYP or TPSSh it is connected with TS2 Ph , whereas ωB97XD prefers TS2. Unfortunately, at present there are no experimental data to discriminate between these two alternatives, and both of them lead to generation of catalytic Cl − . More quantitative comparison between the different computational models may be based on the numeric data gathered in Tables 1 and 2. From the enthalpy values (Table 1), one can notice that certain reaction energies and activation barriers are relatively insensitive to the variations of the methodology (differences within 5 kcal mol -1 ), i.e., P → TS2, P → T, P → T Ph . On the other hand, several other barriers and reaction enthalpy values vary dramatically (differences above 10 kcal mol -1 ) with the model, i.e., S Cl → TS1 Cl , S Cl → P Cl , P Cl → T Cl , P Cl → TS2 Cl_Ph , P Cl → TS2 Cl_homo , and P Cl → T Cl_homo . Notably, the former set includes exclusively the reaction step for the parent system, whereas the latter set exclusively the derivative with Cu-bound Cl − , which might indicate that the presence of inorganic chloride makes the energy of the system more variable with the computational method used. Importantly, most of this variability comes from the change of the density functional. Indeed, if one compares the relative enthalpy values obtained with the three models based on the B3LYP functional, the differences are in almost all cases less than 5 kcal mol -1 , with the exceptions of P Cl → T Cl_Ph, P Cl → TS2 Cl_homo , and P Cl → T Cl_homo , where the spread is 5.8, 5.5, and 5.1 kcal mol -1 , respectively. Nevertheless, including the solvent model or using a larger basis set in the optimization step affects both energies and geometries quite noticeably; hence, whenever it is practically feasible, it is recommended to use these more complete methods.
From analysis of key inter-atomic distances gathered in Table 2 and view of superimposed optimized structures ( Fig.  S1 and S2), it can be inferred that bond lengths and overall conformation of complexes may vary quite substantially with the computational method used. A larger geometrical variance seems to be correlated with a larger spread of energies. For example, one can compare TS2 and TS2 Cl , i.e., two transition states for similar chemical transformations that proceed somewhat differently. For TS2 the variance of critical bond lengths, measured by standard deviation, is 0.08 Å for the O 1 -O 2 bond, 0.06 Å for the C 2 -Cl 1 bond, and 0.14 Å for the O 2 -Cl 1 bond. For TS2 Cl , the corresponding values are 0.14 Å for O 1 -O 2 , 0.19 Å for C 2 -Cl 1 , and 0.06 Å for O 1 -Cl 2 , i.e., larger for two out of three critical bonds. For the preceding species, i.e., P and P Cl , the variance of bond lengths are very comparable. The spread of computed values of enthalpy of activation are 4.6 and 7.3 kcal mol -1 for TS2 and TS2 Cl , respectively. Another example is T Cl_homo , which is characterized by large variability of bond lengths and for which the relative stability, measured by enthalpy of the P Cl → T Cl_homo step, spans a range of 22.0 kcal mol -1 .

Conclusions
Stationary points defining coordinates of several elementary reactions were optimized, and their relative enthalpy values were computed with the use of five computational models. From the comparison of the reaction enthalpy profiles, one can conclude that using a basis set of triple-zeta quality or a PCM solvent model (acetonitrile) for geometry optimization causes the relative values to change by less than 6 kcal mol -1 . On the other hand, changing the functional, even within the hybrid subfamily, brings about more pronounced changes in relative enthalpy values, such that the predicted kinetically preferred product changes. In addition, a correlation between variability of the key bond lengths and enthalpy variance was noticed. In the context of the mechanisms of the reaction between S and O 2 , the TPSSh-D3/BS3//TPSSh-D2/BS1 model is likely at odds with the available experimental data, since it predicts a kinetic preference for homolytic O-O bond cleavage, whereas the experimental evidence supports a process leading to triketone/hypochlorite, presumably via a heterolytic pathway.