The Nature of the Bonding in Symmetrical Pincer Palladacycles

The accuracy of DFT-optimised geometries of the symmetrical pincer palladacycles PdNCN and PdSCS, [ClPd{2,6-(Me 2 NCH 2 ) 2 C 6 H 3 }] and [ClPd{2,6-(MeSCH 2 ) 2 C 6 H 3 }] respectively, has been evaluated by investigating the performance of eight commonly used density functionals with four combinations of basis set, in reproducing their X-ray crystal structures. It was found that whilst the ωB97XD functional performed best over all, the PBE and TPSS functionals performed best when considering the palladium coordination geometry. The role of the donor atom in the stability and reactivity of the symmetric palladacycles, PdYCY, Y = N, S, or P, has been determined using Bader’s Atoms in Molecules method to elucidate the nature of the bonding, and using a model formation reaction, which involves the C-H activation of the pincer ligand YCY by PdCl 2 . The calculations reveal distinct differences in the bond strength and nature of the interaction of Pd with the donor atoms Y, which support differences in the thermodynamic stability.


Introduction
Palladacycles, where a Pd-C bond is intramolecularly stabilised by donor atoms typically from a sulfur, nitrogen or phosphorus donor, are an interesting class of compound. The discovery that they could display extremely high catalytic activity in Suzuki- Miyaura  greater activity than their symmetrical analogues. 9,10 The ability to fine-tune the properties of the ligand, and combine hard and soft Lewis bases with palladium results in hemilability. 11 Thus, determining the nature of the bonding in symmetrical pincer palladacycles is a vital first step to understanding the bonding and reactivity of unsymmetric pincer palladacycles.  12 and PdSCS (II) 13

pincer palladacycles
The structure and bonding of palladacycles can be investigated using density functional theory (DFT), however an appropriate choice of functional, basis set, and relativistic effective core potential (ECP) for Pd, is essential for obtaining accurate structures, and hence thermodynamic and kinetic data. 14 Several papers provide useful benchmarking studies of DFT geometries compared to experimental data including the recent study by Minenkov et al. 15 The authors compared DFT-optimised geometries with experimental crystal structures for 18 ruthenium complexes, and a further 10 transition metal complexes. The results showed that the use of functionals that account for dispersion decrease the statistical error between experiment and theory compared to those that do not include dispersion corrections, with the ωB97XD functional providing the best overall results. However for accuracy around the metal centre the PBE and TPSS functionals also performed very well, and for organic ligands B3LYP performed very well.
Waller et al. 16 investigated the ability of 15 functionals to describe the optimised geometries of 19 second-row transition metal complexes compared to experimental gas-phase data, showing that hybrid functionals such as the PBE hybrid, B3PW91 and B3P86 provide the most accurate geometries. Clearly the choice of functional is dependent on the chemical nature of the structures studied.
The purpose of this paper is two-fold. First, to determine an optimum DFT methodology for the study of pincer palladacycles that combines accuracy with computational speed. This will be achieved by evaluating the performance of a range of functionals in their ability to reproduce the structural features of two experimentally characterised symmetric pincer palladacycles ( Figure 1). Second, to determine the role of the donor atoms, Y = N, S or P, in the symmetric palladacycles, PdYCY. To compare stability and reactivity, the steps involved in a simple formation reaction, which involves the C-H activation of the pincer ligand YCY by PdCl 2 , will be calculated, and the strength and nature of the Pd-L bonding in PdYCY, where L = Y, C, Cl, evaluated.

Computational details
All calculations were performed with the Gaussian 09 package. 17 The neutral spin singlet complexes were studied. The geometry and electronic structure of PdNCN and PdSCS ( Figure 1) were calculated using the same eight density functionals considered by Minenkov et al. 15 The set of density functionals investigated were three generalized gradient approximation (GGA)  functionals: BP86, 18 PBE, 19,20 and B97D, 21 a hybrid-GGA functional: B3LYP, 18,22 two meta-GGA functionals: TPSS 23 and M06L, 24 and two hybrid meta-GGA functionals: ωB97XD 25 and   M06. 26 The basis sets tested were 6-31G(d) and 6-31+G(d,p) for all atoms except Pd for which the two relativistic ECPs, LanL2DZ 27,28 and SDD, 29 were considered. All optimised structures were compared with the X-ray crystal structures 12,13 (Figure 1) obtained from the Cambridge Structural Database (CSD number 720256 for PdNCN and 725124 for PdSCS). 30 The root mean square (RMS) errors between the calculated and experimental structures were calculated using the Quatfit program. 31 The optimum methodology for geometry optimisation in the present work was found to be  ECPs are not recommended for AIM analysis as bond paths cannot be traced, 39 hence the relativistic DGDZVP all electron basis set was used to treat palladium. 40

A. Method validation
To determine an optimum model chemistry the X-ray crystal structures were optimised using a range of functionals and basis sets, and differences between the calculated and observed structural data compared using the Quatfit program. 31 This provides quantitative root mean square (RMS) errors between the experimental and theoretical structures by analysing the difference between the Cartesian coordinates of each atom pair. Hydrogen atoms were excluded due to the known difficulty in accurately determining their position from X-ray crystallography. 15 Two calculations were performed: (i) with equal weighting given to every atom (except hydrogens which are excluded), and (ii) with zero weighting given to atoms that were not directly   When the basis set choice is studied, it is very clear that the SDD ECP outperforms LanL2DZ ( Figures 3 and 4). There is negligible difference between the 6-31G(d) and 6-31+G(d,p) basis sets, however the inclusion of the diffuse functions is desirable because this provides greater flexibility without a substantial increase in computational time. Therefore, the 6-31+G(d,p) [SDD] basis set was chosen as the optimum basis set for geometry optimisation. Minenkov et al. showed that functionals that do not account for dispersion interactions systematically overestimated internuclear distances, whereas when dispersion was included these overestimations were matched with underestimates, resulting in very small mean signed errors (MSE). 15 In the study by Waller et al. it was also found that standard DFT functionals overestimated bond distances. 16   Of key importance is the Pd-L bonding environment.  Taking into account the data from the MUE and MSE of all bond and interatomic distances, all the DFT functionals considered predict expanded structures, but the errors are reasonably small with ωB97XD having the smallest error. The non-hybrid functionals PBE, BP86 and TPSS provide the best results for the Pd-L environment, and provide advantageous calculation time compared to the other functionals which appear higher on Jacob's ladder. 41 The importance of the Pd-L environment, and the desire for an accurate but computationally efficient method for optimisation of large molecular structures led to the choice of the PBE functional, although the use of TPSS or ωB97XD would be equally appropriate. To confirm this, the Gibbs free energy of formation of the palladacycles via PdCl 2 using both the PBE-and ωB97XD-optimised geometries were compared and found to differ by less than 2 kJmol -1 (see Table 1). Therefore the methodology of choice for optimisation in the proceeding section is PBE/6-31+G(d,p)[SDD].
Single point energies using ωB97XD/6-311++G(2df,2p)[SDD] at the optimised geometry were performed to ensure accurate energetics as it has been shown that the inclusion of non-covalent interactions is essential for accurate thermodynamics and kinetics, for example, it has been shown that they constitute a significant proportion of the binding energy in various transition metal phosphine complexes. 32,33

B. Energy and mechanism of formation of PdYCY
The Gibbs free energy of formation was calculated for PdNCN and PdSCS palladacycles,  (Table 1), and show that for these systems there is negligible difference (< 2 kJmol -1 ) between the two datasets. The calculation time using PBE and ωB97XD, starting from the same PdNCN structure and with the same number of iterations, differed by almost a factor of 2.
Therefore as the results are comparable, the PBE functional was preferred for geometry optimisation due to its lower computation time.
Scheme 1. Model formation reaction of palladacycles from PdCl 2 The data in Table 1 show that the PdPCP complex is the most thermodynamically stable, and that the formation of all three complexes is spontaneous (∆G 0 < 0).
The simplest concerted formation reaction, 42 The AIM parameters at the key BCPs are provided in Table 2. The magnitude of the electron density ρ(r) at the BCP indicates the strength of the interaction. The Laplacian of the electron density ∇ 2 ρ(r) indicates the nature of the interaction; ∇ 2 ρ(r) < 0 indicates regions of local electronic charge concentration corresponding to a covalent-type interaction, whereas ∇ 2 ρ(r) > 0 indicates density depletion corresponding to an ionic or closed-shell interaction. 39 Additionally, the sum of the gradient kinetic electron density G(r) and the potential energy density V(r), provides the total electron energy density H(r), which further elucidates the nature of the bond.
The data in Table 2 indicate that, for all PdYCY complexes the Pd-C bond is the strongest, indicated by the largest ρ(r), and the Pd-Cl bond is the weakest. Furthermore, the Pd-P bonds are stronger than the analogous Pd-S and Pd-N bonds at the BCP. This is corroborated by the fact that PdPCP is the most thermodynamically stable (Table 1) In Suzuki-Miyaura catalysis, 47 and other types of catalysis where the active species are Pd(0) species, the active species must be generated from the palladacycle precatalysts. A necessary step in this activation is the reductive elimination of Pd from the Pd-C bond. 46 In an investigation into the reductive elimination of ethane from Pd(PR 3 ) 2 (Me) 2 by Sajith and Suresh it was found that the nature of the bond to be eliminated is more important than their strengths. 40 They found that a more ionic bond results in a higher activation barrier for the reductive elimination whereas a more covalent bond is easier to cleave. Therefore in the present work the nature of the bonds present in PdNCN, PdSCS and PdPCP, is determined, as this may prove important in fine-tuning the reactivity of pincer palladacycles.  Table 2). The bond degree parameter |H(r)/ρ(r)|, which can be interpreted as the total pressure per electron density unit, 51 also indicates that the Pd-P bond is the strongest (|H(r)/ρ(r)| = 0.356) and the Pd-N the weakest (|H(r)/ρ(r)| = 0.161).
The bond ellipticity, ε, can be used to determine the degree of σ and π character in the Pd-L bonds. It measures the extent to which the electron density is unequally distorted away from the bond axis. 56 The low ε of all the Pd-L bonds reflect their predominant σ character, with the Pd-Cl interactions exhibiting a greater π character contribution.
It is clear from this analysis that all of the Pd-L bonds are characterised by low ρ(r), positive ∇ 2 ρ(r) and negative H(r) indicating partial ionic and covalent character. However, it is also clear that there are subtle differences between the Pd-L interactions; the relative magnitude of the parameters indicates that the Pd-P bond is most covalent and the Pd-N bond most ionic. This varying nature of the bonds in each palladacycle could potentially have implications in reactivity. Table 2. AIM parameters (electron density ρ (r), Laplacian of the electron density ∇ 2 ρ (r), total energy H(r), ellipticity ε, and delocalisation index δ(Pd-L)) of the Pd-L Bond Critical Points in PdNCN, PdSCS ( Figure 1) and PdPCP (Figure 2). All values are in atomic units.

Conclusions
The importance of securing an appropriate DFT methodology for the particular system of interest is now well established. Therefore, the experimentally-characterised PdNCN and PdSCS complexes ( Figure 1)  with PdCl 2 , via C-H activation, was determined in order to give a measure of the relative stability and reactivity as a function of the donor atom in the absence of other effects (such as intra-or inter-molecular assisted formation and solvent effects). It was found that both thermodynamically and kinetically PdPCP was the most stable and PdNCN was the least stable. This was supported by a topological analysis of the electron density in the PdYCY complexes. The data suggested that Pd-P bonds were the strongest and the Pd-N bonds the weakest. Furthermore, it was found that the Pd-P bond had more covalent character than the Pd-S and that the Pd-N bond had the most ionic character.
The calculations clearly revealed distinct differences in the bond strength and nature of the interaction of Pd with the donor atoms N, S and P, which supported the differences in the stability, namely that the PdPCP was the most thermodynamically stable (in the absence of external effects). The varying nature of the bonds in each palladacycle could potentially have implications in their reactivity. Exploitation of these differences by considering unsymmetric YCY' complexes, both experimentally and theoretically, are underway.

Electronic Supplementary Information
Energies and Cartesian coordinates of the stationary points along the formation reaction pathway are provided.
the EPSRC UK National Service for Computational Chemistry Software (NSCCS) at Imperial College London in carrying out this work.