Marika Savarese

Éric Brémond

Ilaria Ciofini

Carlo Adamo
email: carlo.adamo@chimie-paristech.fr

Electron Spin Densities and Density Functional Approximations: Open-Shell Polycyclic Aromatic Hydrocarbons as case study

The way how different Density Functional Approximations (DFAs) are able to predict, in open shell systems, spin density that is the difference between the density of electrons with spin  and that of spin  is interesting to be investigated in order to assess their validity.

Here a large panel of functionals were tested on a set composed by seven -radicals expected to amplify DFA's errors in modeling electron delocalization and spin polarization effects due to their extended electronic conjugation coupled with their planar structures. Our results show that generally the DFAs performance follow a systematic improvement in going from semilocal to hybrid functionals. More problematic is, instead the case of double hybrid functionals, where the perturbative contribution to correlation damps the positive effect of the presence of a high percent of exact exchange. More interestingly, differences are observed in the spin delocalization and polarization patterns, thus restraining the possibility of applying some of current DFAs to study chemically relevant properties, like molecular magnetism or charge/electron transport.

1.Introduction

Along the years, the Kohn-Sham approach to Density-Functional Theory 1,2 (KS-DFT) has become a central player on the computational chemistry stage. Indeed, it is widely applied to provide deep insights on chemical systems, spanning from molecular to materials sciences.

Its renowned successes stem from decades of development and implementation improving concomitantly the theoretical background, numerical robustness, and system-size accessibility. [START_REF] Becke | Perspective: Fifty years of density-functional theory in chemical physics[END_REF] The greatest ongoing effort is probably focused on the development of new density functional approximations (DFAs) for the exchange-correlation contribution to the total energy, a key point for the improvement of DFT reliability and accuracy. [START_REF] Mardirossian | Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals[END_REF] In parallel with the proposition of new DFAs, the evaluation of their performances is becoming standardized with the definition and use of benchmark datasets, based on an increasingly large number of chemical systems and properties [START_REF] Goerigk | A Look at the Density Functional Theory Zoo with the Advanced GMTKN55 Database for General Main Group Thermochemistry, Kinetics and Noncovalent Interactions[END_REF][START_REF] Peverati | Quest for a Universal Density Functional: The Accuracy of Density Functionals Across a Broad Spectrum of Databases in Chemistry and Physics[END_REF][START_REF] Řezáč | A Well-balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures[END_REF][START_REF] Karton | W4-11: A high-confidence benchmark dataset for computational thermochemistry derived from first-principles W4 data[END_REF] . Such datasets are also used to train semi-empirical functionals, that is approaches tuned on selected properties and systems [START_REF] Zhao | The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-class Functionals and 12 Other Functionals[END_REF][START_REF] Santra | Minimally Empirical Double-Hybrid Functionals Trained against the GMTKN55 Database: revDSD-PBEP86-D4, revDOD-PBE-D4, and DOD-SCAN-D4[END_REF] . This process is not always linear, and it could happen that basic quantities are let out from DFA development or benchmarking. For instance, structural parameters (i.e. bond lengths), systematically considered in the '90 [START_REF] Johnson | The performance of a family of Density Functional methods[END_REF][START_REF] Adamo | Exchange functionals with improved long-range behavior and adiabatic connection methods without adjustable parameters: The mPW and mPW1PW models[END_REF] , were not necessarily taken into account in standard benchmarks and only recently their interest was renewed (see for instance reference 13 and 14). Electron density followed a similar trajectory and only recently was again raised at the center of the scientific debate concerning its place in the development of DFAs. It was argued that functionals developed until the early 2000s provide electron densities closer to the exact ones, while more recent DFAs provide larger errors for such fundamental quantity [START_REF] Medvedev | Density functional theory is straying from the path toward the exact functional[END_REF] . While some strong arguments were lifted against such considerations [START_REF] Kepp | Density functional theory is straying from the path toward the exact functional[END_REF][START_REF] Korth | Density Functional Theory: Not Quite the Right Answer for the Right Reason Yet Angew[END_REF] , it should be considered that 3D density evaluation on a numerical grid makes difficult and poorly accurate its comparison to a reference [START_REF] Datta | Accurate Prediction of Hyperfine Coupling Tensors for Main Group Elements Using a Unitary Group Based Rigorously Spin-Adapted Coupled-Cluster Theory[END_REF] , not only for small systems like atoms, but also for molecules [START_REF] Wang | Electron Densities of Several Samll Molecules Ac Calculated from Density Functional Theory[END_REF] .

It has also been argued that the energy functional itself provides a better measure of density errors [START_REF] Sim | Quantifying Density Errors in DFT[END_REF] .

A more challenging playground is represented, in our opinion, by spin densities, that is the difference, in open shell systems, between the density of electrons with spin  and that with spin  (  -  ). This property is particularly intriguing for several reasons. An accurate description of spin densities is required for the prediction of several spectroscopic parameters (e.g. EPR 21 ) as well as a number of properties like magnetic exchange, electron and charge transport and more generally properties of open shell compounds. It involves an accurate description of correlation effects in order to model both delocalization of the unpaired electron(s) and spin polarization, that is the possibility of regions characterized either by positive and negative spin densities [START_REF] Mcconnell | Theory of Isotropic Hyperfine Interactions in π Electron Radicals[END_REF] . Information about spin density delocalization and spin polarization is either obtained from the analysis of the spin density distributions or, in a more compact form, from the corresponding computed atomic spin population that can be derived from a Mulliken analysis for instance. In this case, information about the spin density is condensed by a single value representing the excess of positive () or negative () spin on an atomic site. This latter will be defined in the following as Mulliken spin population even though sometimes they are incorrectly referred in literature as Mulliken spin densities.

It should be recalled that, in the framework of DFT, open-shell molecules are treated in an unrestricted KS approach, even if, according to the Hohenberg-Kohn theorems, densities of  and  electrons are not needed to calculate any molecular properties. Furthermore, these densities (their sum or difference) are also used as variable in some functionals, such as LYP [START_REF] Lee | Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[END_REF] , just to mention one.

Despite such relevant role, few systematic investigations were carried out on spin densities.

In middle '90 [START_REF] Adamo | Structure and magnetic properties of benzyl, anilino, and phenoxyl radicals by density functional computations[END_REF][START_REF] Qin | Density-functional methods give accurate vibrational frequencies and spin densities for phenoxyl radical[END_REF][START_REF] Jensen | Density Functional and MP2 Calculations of Spin Densities of Oxidized 3-Methylindole: Models for Tryptophan Radicals[END_REF] , spin densities in organic radicals were investigated for a limited number of DFAs, with some lately works dating to first years of the 2000s [START_REF] Adamo | The mechanism of spin polarization in aromatic free radicals[END_REF][START_REF] Filatov | Calculation of spin-densities within the context of density functional theory. The crucial role of the correlation functional[END_REF] . More recently, the performances of different DFAs were tested for some transition metal complexes or organic radicals, in view of specific interests [START_REF] Riffet | Benchmarking DFT and TD-DFT Functionals for the Ground and Excited States of Hydrogen-Rich Peptide Radicals[END_REF][START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Jacob | Can DFT Accurately Predict Spin Densities? Analysis of Discrepancies in Iron Nitrosyl Complexes Katharina Boguslawski[END_REF][START_REF] Das | Polyradical Character of Triangular on-Kekul Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study[END_REF] . In general, it was argued [START_REF] Adamo | Validation of self-consistent hybrid density functionals for the study of structural and electronic characteristics of organic  radicals[END_REF] that the introduction of a fraction of exact-like exchange (EXX) into the functional leads to an improvement in spin densities, so that global hybrid functionals are, generally speaking, better than their parent semilocal approximations. Concerning, transition metal complexes, the situation could be more entangled, due an eventual multireference character (see reference 31 and discussion therein).

The observed deviations on properties of open shell systems are related to the approximate nature of the exchange-correlation functional used, which leads to the spurious selfinteraction and related delocalization errors [START_REF] Schattenberg | Lessons from the Spin-Polarization/Spin Contamination Dilemma of Transition-Metal Hyperfine Couplings for the Construction of Exchange-Correlation Functionals[END_REF] . These latter are in turn strictly connected to fractional electron/spin occupation [START_REF] Li | Localized Orbital Scaling Correction for Systematic Elimination of Delocalization Error in Density Functional Approximations[END_REF] .

In order to investigate how spin densities are reproduced by functionals spanning a large timeframe, we propose here to consider a set composed by seven open-shell polycyclic aromatic hydrocarbons (PAHs). Thanks to their planar structure and to the fully delocalized electron and spin density on their  system, PAHs are challenging systems for DFT, particularly when their size increases and induces delocalization errors.

A large number of functionals was tested in the present work, spanning over all the so-called Perdew's ladder. [START_REF] Perdew | Some Fundamental Issues in Ground-State Density Functional Theory: A Guide for the Perplexed[END_REF] The lowest rungs gather the semilocal (or non-hybrid) approximations, including the Generalized Gradient Approximations (GGAs). Highest rungs concern functionals having a nonlocal exchange and correlation contribution through the addition of EXX and second order perturbation theory (PT2) correlation, respectively. They belong to the global or range-separated hybrids (GH and RSH) and to double hybrid (DHs) classes of approximation. Out of introducing more and more nonlocal ingredients into their energy expressions, these successive approximations provide an (often systematic) increasing accuracy [START_REF] Brémond | Systematic Improvement of Density Functionals Through Parameter-free Hybridization Schemes[END_REF][START_REF] Brémond | Nonempirical Double-Hybrid Functionals: An Effective Tool for Chemists[END_REF][START_REF] Mehta | Semi-empirical or non-empirical double hybrid density functionals: which are more robust?[END_REF][START_REF] Martin | Empirical Double-Hybrid Density Functional Theory: A 'Third Way' in Between WFT and DFT Isr[END_REF] . All DFAs can be also classified as semiempirical, that is functionals containing parameters fitted on selected reference properties, or as non-empirical, that is approaches whose internal coefficients are determined using only theoretical arguments. It is worth to mention that electron density was, at the best of our knowledge, never considered in any training set for semi-empirical functionals, except in a recent work based on machine learning [START_REF] Schmidt | Machine Learning the Physical Nonlocal Exchange-Correlation Functional of Density-Functional Theory[END_REF] . As a consequence, it is absolutely not granted that even this class of functionals will enable for an accurate estimate of spin properties.

2.Computational Details

The neutral structures of the 7 PAHs reported in Figure 1 were fully optimized with the PBE0 GH functional [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF][START_REF] Ernzerhof | Assessment of the Perdew-Burke-Ernzerhof exchangecorrelation functional[END_REF] and the cc-pVTZ basis set, [START_REF] Dunning | Encyclopedia of Computational Chemistry[END_REF] a methodology proved as reliable to compute structural parameters. [START_REF] Brémond | Benchmarking Density Functionals on Structural Parameters of Small-/Medium-sized Organic Molecules[END_REF] Once optimized, the Mulliken atomic spin populations of the PAHs in their Doublet (D) or Triplet (T) spin states were computed for the neutral, cationic or anionic form for the molecules as reported in Figure 1. A large and representative panel of density functionals, reported in Table 1, was tested. Electron densities of the non-hybrid, globalhybrid, and range-separated hybrid classes of approximation are computed in a self-consistent fashion following the Unrestricted Kohn-Sham (UKS) procedure. Densities of Double-hybrid functionals were produced by taking into account the response of the molecular orbitals to the external perturbation via the Z-vector technique once the self-consistent KS density generated. [START_REF] Wiberg | Analysis of the effect of electron correlation on charge density distributions[END_REF][START_REF] Handy | On the evaluation of analytic energy derivatives for correlated wave functions[END_REF] For selected cases, the Restricted-Open shell Kohn-Sham (ROKS) method was applied for comparison.

All the calculations of spin densities were carried out using the EPR-II basis set, which was specifically developed to describe spin density at the nuclei, relevant for the description of EPR parameters, in conjunction with DFAs. [START_REF] Barone | Recent Advances in Density Functional Methods (Part I)[END_REF] Mulliken spin populations computed at the QCISD/EPR-II level were taken as reference.

Clearly, the quality of the performance analysis depends on the reliability of the reference values. While CASSCF (with large active space) and higher level approaches were used in literature for studying spin densities in small radicals (such as FeNO 2+ ) [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Radoń | Electronic Structure of Selected {FeNO} 7 Complexes in Heme and Non-Heme Architectures: A Density Functional and Multireference ab Initio Study[END_REF][START_REF] Boguslawski | Accurate ab Initio Spin Densities[END_REF] , some papers report multi-references calculations on PAHs, including phenalenyl and triangulene, two of the molecules considered in the present paper [START_REF] Das | Polyradical Character of Triangular on-Kekul Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study[END_REF][START_REF] Sandoval-Salinas | Triangular graphene nanofragments: open-shell character and doping[END_REF] . These works showed, using limited active spaces, that some polyradical character could be expected, more for phenalenyl than for triangulene. Therefore, it could be likely that a large spin contamination could affect the QCISD results. An easy way to assess the quality of QCISD results is to compare their unrestricted and restricted-open shell energies: if they are close, the UQCISD results are reliable, since the contamination has been removed. Missing a ROQCISD approach in Gaussian, we have considered the Coupled Cluster (CC) method, which has been showed to be quite independent on the choice of restricted or unrestricted reference [START_REF] Stanton | On the extent of spin contamination in openshell coupledcluster wave functions[END_REF] . Furthermore, the CCSD approach eliminates the contamination coming from the S+1 state [START_REF] Schlegel | Møller-Plesset perturbation theory with spin projection[END_REF] . The energy difference between ROCCSD and UCCSD calculation for phenalenyl is 0.7 kcal/mol and 0.3 kcal/mol for triangulene. Since differences up to 1.2 kcal/mol have been used as discriminating criterium for similar problems [START_REF] Menon | Consequences of Spin Contamination in Unrestricted Calculations on Open-Shell Species: Effect of Hartree-Fock and Møller-Plesset Contributions in Hybrid and Double-Hybrid Density Functional Theory Approaches[END_REF] , the UCCSD calculations can be considered reliable. Since the UCCSD and UQCISD spin populations are very close (see Tables 2 and3), these latter will then be considered as references in the following.

It is worth to notice also that QCISD approach is also commonly used as reference in large organic radicals [START_REF] Sui | Stabilities and Spin Distributions of Benzannulated Benzyl Radicals[END_REF][START_REF] Korcsok | ultiscale Modeling of Interfacial Oxidation Mechanism at Air/Organic Interface: Reactions of CH2═CH-Terminated Self-Assembled onolayer with OH•, O3, and HO2[END_REF] and even in systems dominated by large dynamic correlation effects, such as noncovalent dimers [START_REF] Sedlak | Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes[END_REF] and it could provide results very close to CC [START_REF] Chen | Toward a less costly but accurate calculation of the CCSD(T)/CBS noncovalent interaction energy?[END_REF] .

All the computations are performed with the release D.01 of the Gaussian'09 code. [START_REF] Frisch | Gaussian'09 Revision D.01[END_REF]

Results and Discussion

As mentioned in the introduction, it is expected that the selected PAHS represent an interesting test for the evaluation of spin densities and spin polarization effects. Indeed, the strong electronic conjugation coupled with the planarity of the system are likely to amplify electron delocalization and spin polarization effects which, in turn, will be largely affected by self-interaction and delocalization error of the DFAs. The erroneous behavior of the functionals should be, therefore, more easy to quantify.

phenalenyl and triangulene

Our survey starts with the phenalenyl and triangulene (Figure 1), two highly symmetric PAHs which are characterized by one and two unpaired electrons, respectively. Triangulene has a high spin-state -and open-shell radical character-in its ground state, accordingly to the socalled Ovchinnikov's rule [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Ovchinnikov | Multiplicity of the ground state of large alternant organic molecules with conjugated bonds[END_REF] . The spin-state was recently confirmed by high-level post-Hartree-Fock (post-HF) approaches and it is also in agreement with experimental findings [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Sandoval-Salinas | Triangular graphene nanofragments: open-shell character and doping[END_REF] . Of note being the molecule in its high-spin state, that is a triplet state, the system is expected to have very low spin-contamination.

To discriminate the density functionals listed in Table 1, we compare the Mulliken spin population of C atoms with those obtained at the QCISD (reference) level of theory, calculating the mean absolute deviation (MAD) over the total number of C atoms of each molecule. Figure 2 depicts these statistics computed with the 37 density functionals which are ranked according to the fraction of EXX. A first look at Figure 2 reveals that, irrespective of the method, the MADs calculated for the single radical PAH (i.e., phenalenyl) are very close to the one estimated for the diradical (i.e., triangulene). The difference between their MADs lies around ∼0.01 a.u., the single radical being more accurately modeled than the diradical.

Few exceptions can be noticed like BH&HLYP, LC-BLYP or TPSS-QIDH which are inversing this trend.

Regarding the different families of density functionals, the best candidates belong to the global, range-separated and double hybrids classes (Figure 2). In particular, ωB97X is found as the best performing for this benchmark subset, showing a MAD lower than 0.01 a.u. for both PAHs. It is closely followed by the M05 and M05-2X global hybrids and the TPSS0-DH double hybrid with an error of about 0.01 a.u.. In comparison, the semilocal approximations provide poor performance, the best one being MN12-L, which belongs to the latest generations of Minnesota density functionals (error of about 0.06 a.u. for both molecules).

Looking more in detail, the performance of the semilocal approximations increases, as it could be expected, in going from Local density to generalized gradient, to meta generalized gradient approximations (i.e., LDA, GGA and meta-GGA, respectively). This improvement is independent of the level of parameterization of the density functional, being the same for semiempirical and non-empirical functionals. In case of phenalenyl, for instance, the MADs range between 0.13 a.u. (SVWN5) and 0.06 (MN12-L), while the TPSS meta-GGA performs slightly better than its parent PBE GGA (0.09 versus 0.11 a.u.). It is interesting also to remark a systematic improvement for the successive generations of semilocal Minnesota density. The MADs of M06-L, M11-L and MN12-L are indeed calculated as 0.09, 0.07 and 0.06 a.u., respectively.

It has been argued that the introduction of a fraction of EXX leads to an improvement of spin densities and related properties [START_REF] Adamo | Validation of self-consistent hybrid density functionals for the study of structural and electronic characteristics of organic  radicals[END_REF] . While this trend is clear for functionals casting the same GGA part, this feature is not a systematic guarantee of the good behavior of the method (Figure 2). The parameterization of the exchange-correlation density-based terms plays, indeed, an important role which is well depicted here by the difference observed between M05 and M06 global hybrids. Both of them share the same fraction of EXX (i.e., 27%) but M05 overperforms M06 with MADs calculated as 0.01 and 0.07 a.u., respectively. With 56% and 54% of EXX, M05-2X and M06-2X follow the same behavior with errors of about 0.03 and 0.01 a.u., respectively. M06-HF which casts 100% of EXX, decreases again the MAD with an error of about 0.02 a.u., however the high percentage of EXX is not sufficient to reach the excellent performance of the M05 family of density functionals. On the other side, nonempirical global hybrids like PBE0 do not reach the excellent results provided by M05 but fairly lies around 0.04 a.u.. The strong importance of the parameterization of density functionals is also depicted by the behavior of range-separated hybrids. If LC-BLYP, LC-PBE and LC-ωPBE provide performance in line with global hybrids casting a low fraction of EXX (the MADs being larger than 0.06 a.u.), CAM-B3LYP, M11 and the ωB97 family provide an excellent estimate of the density with errors lower than 0.02 a.u.. However, highly parameterized density functionals like N12SX and MN12-SX are disappointing with errors reaching up to 0.10 a.u.

(for MN12-SX), underlining here again that a parameterization is not granting good results for spin density description.

In double hybrids functionals, instead, the beneficial effect of a large quantity of EXX is partially cancelled by the PT2 contribution. PAHs are weakly correlated systems which induce an artificial symmetry-breaking of the unrestricted Kohn-Sham orbitals used to compute the PT2 correlation energy term. [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF] It results from this failure a dramatic estimate of the density-based properties. PBE0-DH and TPSS0-DH are not much impacted by this error since they mix a low fraction of PT2 correlation (12.5%). B2-PLYP, mPW2PLYP, PBE-QIDH and TPSS-QIDH introduce however more than 25% of PT2 correlation, so they struggle to reach nice performance in the field. We notice that the observations made above concerning the performance ranking between GGA and meta-GGA are nevertheless respected with double hybrids, TPSS-QIDH (0.06 a.u.) being more accurate than PBE-QIDH (0.08 a.u.). We conclude also that a compromise between the fractions of EXX and PT2 exists since in going from B2-PLYP (EXX = 52%; PT2 = 27%) to PBE-QIDH (EXX ≈ 69%; PT2 ≈ 33%) the MAD is reduced of 0.05 a.u., the former providing one of the largest error on this test set (0.13 a.u.).

Starting from this general picture on the performances of the different DFAs, it is worth to better analyze the chemical effects underpinning the found average errors. Before we analyze the numerical results, it is worth to recall that phenalenyl allows for a large delocalization of the unpaired electron, as it clearly appears when writing the possible resonance structures reported in Figure 3. All the atoms involved in the delocalization of the unpaired electron are expected to bring a positive spin density. The spin-polarization then arises from the fact that the unpaired electron interacts in a different way with the two electrons of a spin-paired bond or inner shell, since the exchange is operative only for electrons with the same spins. This leads to alternance of positive () and negative () spin densities as those reported in Figure 4. Please note that -spin contributions are localized on the atoms contributing to the singly occupied molecular orbital (SOMO), as it can be found, for instance, using a simple Hückel approximation [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF] , except for the central carbon atom. This pattern should be respected by any UKS calculations, regardless of the used functionals. The spin density on the central atom then arises from a spin-polarization effect.

In Figure 5 are reported the spin density maps for the phenalenyl radical, computed with 4 functionals belonging to the GGA, GH and DH families, casting the PBE model (i.e. PBE, PBE0, PBE0-DH and PBE-QIDH). From these plots, it is clear that all of them give the expected spin pattern, sketched in Figure 4. However, PBE and PBE-QIDH have a small or null spin population on the central atom (C4, see Figure 4 for atom labelling), whereas all the other functionals have a significant -spin contribution (blue surfaces in Figure 5) in agreement with the qualitative sketch of Figure 3. The same two functionals show a significantly lower spin polarization, the contribution of -spin on C atoms being always quite low (red surfaces in Figure 5). While the PBE behavior is somehow expected and it could be related to the quality of the exchange and correlation contributions, more difficult is to explain the behavior of the PBE-QIDH functional. Indeed, its large quantity of EXX should lead to the accurate spin-polarization patterns. Figure 6 reports the spin density map for the PBE-QIDH and its hybrid counterpart obtained from the SCF KS orbitals, giving the so-called PBE-QIH approach. The spin density obtained with this last functional clearly shows the expected pattern, with larger contributions ( or ) at each carbon atom, including the central C4. The difference observed between PBE-QIDH and PBE-QIH can thus be directly ascribed to the electron density relaxation (through the Z-vector approach) which significantly damps the spin-polarization effects. The effect is lower for PBE0-DH than for PBE-QIDH due to the smaller respective PT2 contributions (12.5 vs. 33%).

It should be noticed that these results are in qualitative agreement with previously reported post-HF studies [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Das | Polyradical Character of Triangular on-Kekul Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study[END_REF] , both predicting a large localization of the spin density on the edge carbon atoms.

It is also interesting to directly analyze the spin Mulliken populations, reported in Table 2. As already discussed for the statistics reported in Figure 2, GGAs and DHs are the functionals providing large deviations (the MAD being between 28.0 and 13.2 %), with the only exception of PBE-DH (4%). The best results are obtained with B97X (0.8%) and M05(1.2%). More in details, it should be noticed that some functionals provide acceptable deviations but overestimate the spin polarization effects. This is particularly evident for BH&HLYP (4.6 %), which significantly miscalculates all the spin populations. Also, the already discussed behavior of PBE and PBE-QIDH clearly appears from the spin populations, the two functionals having a very similar behavior. However, an almost zero spin population on C4 is found for PBE, while PBE-QIDH indicates a small negative population. The MAD for this functional is comparable with that of LC-PBE (13.7 vs. 11.3 %), but this latter gives an overall better spin polarization pattern, the C4 bringing an -spin as it should. This result underlines how a small MAD may not necessarily imply the presence of a correct spin pattern. Furthermore, the comparison between PBE-QIH and PBE-QIDH further illustrates the effect of electron relaxation. Indeed, the spin population at the PBE-QIH level are significantly higher than those obtain at the PBE-QIDH level, showing the largest deviations among all the considered functionals (26.4 %). Indeed, the high EXX contribution leads to a significant spin contamination (<S 2 >= 1.22) as for the UHF case (<S 2 >= 1.95).

As discussed some years ago, the spin density can be decomposed into a direct contribution arising from only the singly occupied molecular orbital (SOMO) and an indirect contribution, which is due to the spin polarization [START_REF] Adamo | The mechanism of spin polarization in aromatic free radicals[END_REF] . The direct contribution can be easily computed at Restricted Open Shell (RO) level. The spin populations computed with the ROKS approach for selected functionals are also reported in Table 2. As it can be observed from these data, the results are the same independently from the considered approach (HF, GGA, GH or DHs) and follow the symmetry requirement of the Single Occupied Molecular Orbital (SOMO) [START_REF] Riffet | Benchmarking DFT and TD-DFT Functionals for the Ground and Excited States of Hydrogen-Rich Peptide Radicals[END_REF] . This symmetry requirement is stronger than in other aromatic molecules, where significant differences are observed between ROHF and ROKS results [START_REF] Suaud | Singly Occupied MOs in Mono-and Diradical Conjugated Hydrocarbons: Comparison between Variational Single-Reference, π-Fully Correlated and Hückel Descriptions[END_REF] .

It is interesting to notice that ROPBE and UPBE results are very close, thus pointing out the role the approximate functional in the poor representation of the spin polarization. The difference observed between ROPBE-QIDH and UPBE-QIDH results well evidences the already discussed relaxation effects in the latter which reduce the spin polarization magnitude.

The next molecule, triangulene has an even larger number of resonance structures for the two unpaired electrons, as can be seen from Figure 3. The spin pattern is, however, similar to that found for phenalenyl, since triangulene can be obtained from this latter by adding three aromatic rings at its base. Therefore, it is not surprising that the same trends can be evidenced for triangulene, even if they are moderated by the presence of two unpaired electrons which reduce the flexibility of the spin-polarization pattern. Indeed, all the considered unrestricted approaches give a significant -spin population at the C4 atom, albeit the magnitude depends on the functionals. As before, PBE and PBE-QIDH provide very close spin-density maps (see Figure 5) and this behavior can be ascribed to the density relaxation as clearly appears from the comparison of the maps obtained with and without PT2 contribution (Figure 6). Table 3 collects the MADs for the spin-population of triangulene, which are slightly lower than in case of phenalenyl, albeit the global trend of functional performances is the same. Indeed, B97x and M05 provide, once again, the lowest MADs (0.4 and 0.6%). Other functionals, provide relatively small errors but, due to a compensation between  and  spin contributions, wrong magnitude of the spin polarization. This is the case of BH&HLYP (2.5 %), LC-PBE (5.3 %) and PBE-QIDH (5.6%), for which the effect is more evident than in the case of phenalenyl.

In going from UKS to ROKS, the direct contribution to the spin density is evidenced and, as before, very close values for the spin populations are obtained independently from the chosen computational approach. As for phenalenyl, there is also a qualitative agreement with the data reported in literature, concerning the largest distribution of the spin densities on the edge carbons [START_REF] Trinquier | Spreading out spin density in polyphenalenyl radicals[END_REF][START_REF] Das | Polyradical Character of Triangular on-Kekul Structures, Zethrenes, p-Quinodimethane-Linked Bisphenalenyl, and the Clar Goblet in Comparison: An Extended Multireference Study[END_REF] .

Other PAH radicals

Phenalenyl and triangulene being difficult test cases for DFT, we continue this DFAs' survey with the radical cation and anion forms of phenanthrene (Figure 1). The skeleton of this PAH is composed by 3 fused benzene rings whose electron density is delocalized on its whole planar structure. Figure 4 reports the MADs of the C atom Mulliken spin populations of these two systems computed with the 37 investigated density functionals, taking QCISD populations as a reference. As a first global overview, we remark that the error range is here drastically decreased by a factor 2 or 3, depending on the method, in direct comparison with the previous systems. Moreover, we observe that the density of the anions form is systematically better modeled than that of the cations.

Despite of these few differences, all the trends and comments mentioned above can be applied on these two systems (Figure 7). More precisely, the meta-GGA class of approximation is a general improvement with respect to the GGAs, also when these approaches are involved in the hybrid and double-hybrid classes of approximation. The empirical parameterization of the density-based terms can induce large performance variations, sometimes more important than the fraction of EXX. However, global hybrids mixing more than 25% of EXX are often the guarantee of an excellent performance for this property. As example, with 25, 27 and 56% of EXX, PBE0, M05 and M05-2X provide MADs lower than 0.01 a.u.. Double hybrids mixing a large fraction of PT2 correlation tend to provide poor performance. PBE0-DH and TPSS0-DH are the only ones with a MAD lower than 0.01.

In order to validate the previous observations, the investigation is extended to larger molecular systems, i.e. anthracene, pyrene, perylene and coronene (Figure 1). Since QCISD densities are too computationally demanding for these systems, PBE0 densities, which are up to now recognized as a good compromise between accuracy and computational time, are taken as a reference for comparison. This is, of course, an arbitrary choice and other functionals, such as M05, could be also considered as reference in view of their comparable performances. The analysis also was restricted to a subset of 10 density functionals part of the GGA, meta-GGA, global-hybrid and range-separated hybrid classes of approximations (Figure 8).

Anthracene is a regioisomer of phenanthrene. Deviations with respect to the PBE0 spin density populations are here very low (lower than 0.02 a.u.). Only BLYP, LC-BLYP and M06-HF reach nevertheless this maximum, the other density functionals providing deviations close to 0.01 a.u.. Extending the size of the system confirms this trend and shows that semilocal approximations (e.g., BLYP and M06-L) or hybrid density functionals mixing a large fraction of EXX (e.g., LC-BLYP and M06-HF) are less accurate to model spin densities.

Conclusions

In this investigation, we survey a selection of 37 density functionals covering the semilocal, hybrid and double-hybrid classes of approximation to model the spin density populations of open-shell PAHs molecular systems. This property is chosen as a possible index allowing monitoring and quantification of the quality of the electron density produced by the density functionals. From this analysis, we show that the accuracy with respect to this property follows globally the approximation ranking provided by the Perdew's ladder.

Regarding semilocal approximations, the meta-GGA approximation is found as an improvement with respect to the GGA one which is itself an improvement with respect to LDA. The hybrid approximation including both global and range-separated hybrid density functionals further improves with respect to the previous ranks. However, double hybrids which are usually the best performing in energy benchmarking, provide here a more moderate performance. Only double hybrids mixing a low fraction of PT2 correlation (like PBE0-DH and TPSS0-DH) are found as excellent with respect to this property. Double hybrids with high fractions of PT2 correlation (more than 25%) provide poor performance due to the artificial symmetry-breaking occurring in these open-shell PAHs.

Beyond the analysis of the numerical performances, our results point out that, in some cases, significantly variations of the spin delocalization and polarization patterns, typically with GGA and DH functionals. This outcome may be considered as a warning for the applications of the concerned DFAs in chemical fields where such properties (spin delocalization and polarization) are relevant, such as molecular magnetism or charge/electron transport.

Figure captions Figure 1 .

 captions1 Figure captions Figure 1. Representation of the chemical structures, molecular formula and names of the selected PAHs. Corresponding computed spin state is reported: D (doublet); T (triplet).

Figure 2 .

 2 Figure 2. Mean absolute deviations (a.u.) of the C atom Mulliken spin populations computed for a panel of exchange-correlation density functionals on the (red) phenalenyl and (green) triangulene. The density functionals are ordered according to their classes of approximation and their fraction of exact-like exchange.

Figure 3 .

 3 Figure 3. Electronic resonance structures showing unpair electron delocalization for phenalenyl (up) and triangulene (down). Only non-symmetric resonance structures are reported for the sake of space.

Figure 4 .

 4 Figure 4. Atom numbering and sketch of the spin polarization mechanism for phenalenyl (left) and triangulene (right). Blue:  spin; Red:  spin.

Figure 5 .Figure 6 .Figure 7 .

 567 Figure 5. Spin density map (contour 0.005 a.u.) for the phenalenyl (up) and triangulene (down). Blue:  spin density; Red:  spin densityFigure 6. Spin density map (contour 0.005 a.u.) for the phenalenyl (left) and triangulene (right), obtained without (up) and with (down) MP2 contribution in the indicated double hybrid functionals. Blue: -spin density; red: -spin density Figure 7. Mean absolute deviations (a.u.) of the C atom Mulliken spin populations computed for a panel of exchange-correlation density functionals on the phenanthrene (red) cation and (green) anion PAHs.

Figure 8 .

 8 Figure 8. Mean absolute deviations (a.u.) of the C atom Mulliken spin populations computed for a selected panel of exchange-correlation density functionals on the (top left) anthracene, (top right) pyrene, (down left) perylene and (down right) coronene (red) cation and (green) PAHs.

Figure 1 .

 1 Figure 1.

Figure 7 32 Figure 8 Graphical

 328 Figure 7

Table 1 .

 1 List of the exchange-correlation functionals considered in the present work. EXX denotes the percent of exact-like exchange; PT2 denotes the percent of correlation contribution computed at the second-order perturbation theory level.

	non-hybrid	global hybrids		range-separated hybrids	double-hybrids
		ref.		EXX ref.		EXX	ref.		EXX PT2 ref.
	SVWN5	45,46 B3LYP	20	55-57 B97	0-100 66	B2-PLYP	53	27	74
	BLYP	47,23 B97-2	21	58	B97X	16-100 66	mPW2PLYP 55	25	75
	PBE	48	M05 a	27	59	B97XD	22-100 67	PBE-QIDH	69	33	76
	SOGGA11 49	M06 a	25	9	CAM-B3LYP 19-65 68	PBE0-DH	50	12.5 77,78
	N12	50	PBE0	25	42,43 LC-PBE	0-100 69,70 TPSS-QIDH a 69	33	76
	TPSS a	51	TPSSh a	10	60	LC-BLYP	0-100 71	TPSS0-DH a	50	35	77,78
	M06-L a	52	BH&HLYP	50	61	LC-PBE	0-100 71
	M11-L a	53	BMK a	48	62	M11 a	43-100 72
	MN12-L a 54	M05-2X a	56	63	MN12-SX a	25-0	73
			M06-2X a	54	9	N12-SX	25-0	73
			M06-HF a	100	64
			SOGGA11-X 40	65
	a)These functionals also depend on the kinetic energy density ()

Table 3 .

 3 Computed spin densities (a.u.) for triangulene (see Figure4for atom labelling). Mean Absolute Deviation (MAD, a.u and %) are computed with respect to the QCISD reference. All calculations were done at Unrestricted Kohn-Sham level, excepted those preceded by the prefix

	RO, which were carried out at Restricted Open shell Kohn-Sham level. See Figure 4 for atom
	numbering.
		C1	C2	C3	C4	C5	C6 MAD MAD %
	HF	-0.685 0.752 -0.600 0.638 0.768 -0.659 0.426	27.6
	ROHF	0.000 0.184 0.000 0.000 0.222 0.000 0.190	12.3
	PBE	-0.085 0.267 -0.053 0.087 0.294 -0.047 0.119	7.7
	ROPBE	0.000 0.185 0.000 0.000 0.217 0.000 0.191	12.3
	MN12L	-0.144 0.322 -0.105 0.124 0.355 -0.096 0.067	4.3
	PBE0	-0.173 0.338 -0.128 0.143 0.370 -0.119 0.046	3.0
	ROPBE0	0.000 0.185 0.000 0.000 0.215 0.000 0.191	12.3
	M05	-0.194 0.359 -0.176 0.177 0.412 -0.179 0.009	0.6
	M06-HF	-0.218 0.396 -0.214 0.197 0.422 -0.149 0.018	1.1
	BH&HLYP	-0.266 0.410 -0.214 0.225 0.445 -0.214 0.038	2.5
	B97	-0.222 0.376 -0.193 0.198 0.423 -0.211 0.013	0.8
	B97X	-0.208 0.366 -0.183 0.190 0.415 -0.191 0.006	0.4
	LC-PBE	-0.312 0.447 -0.256 0.269 0.483 -0.268 0.082	5.3
	ROLC-PBE	0.000 0.183 0.000 0.000 0.224 0.000 0.190	12.3
	PBE-QIDH	-0.143 0.316 -0.090 0.072 0.356 -0.053 0.086	5.6
	ROPBE-QIH 0.000 0.185 0.000 0.000 0.218	0.00 0.190	12.3
	PBE-QIH	-0.431 0.542 -0.360 0.379 0.572 -0.381 0.187	12.1
	PBE0-DH	-0.191 0.355 -0.141 0.144 0.388 -0.124 0.034	2.2
	CCSD	-0.207 0.372 -0.174 0.184 0.410 -0.168 0.005	0.3
	QCISD	-0.212 0.374 -0.178 0.191 0.414 -0.176

Acknowledgement

This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (CoG STRIGES; grant agreement No 648558) and by the ANR Grant E-Storic. E.B. gratefully acknowledges the GENCI-CINES for HPC resources (Projects AP010810360 and A0040810359). ANR (Agence Nationale de la Recherche) and CGI (Commissariat à l'Investissement d'Avenir) are gratefully acknowledged for their financial support of this work through Labex SEAM (Science and Engineering for Advanced Materials and devices) ANR-10-LABX 096, ANR-18-IDEX-0001.


Pierre-Adrien Payard

Alessandro Luca

Laurence Perego

Ilaria Grimaud

Ciofini

A DFT Protocol for the Prediction of 31 P NMR Chemical Shifts of Phosphine Ligands in Transition Metal Complex-es

A DFT Protocol for the Prediction of 31 P NMR Chemical Shifts of Phosphine Ligands in Transition Metal Complexes

Pierre-Adrien Payard,* [a] ‡ Luca Alessandro Perego,* [a][b] ‡ Laurence Grimaud [a] and Ilaria Ciofini* [b] [a] Laboratoire des Biomolécules, LBM, Département de Chimie, École Normale Supérieure, PSL University, Sorbonne Université, CNRS, 75005 Paris, France. [b] PSL Research University, Institute of Chemistry for Health and Life Sciences, ICLeHS, CNRS-Chimie ParisTech, 11 rue P. et M. Curie, F-75005 Paris 05 (France). DFT • 31 P NMR • phosphine • basis set • first row organometallic compounds • theoretical prediction ABSTRACT: While 31 P NMR is a major technique to characterize phosphine-ligated transition metal complexes -which are ubiquitous in catalysis -31 P NMR chemical shifts are difficult to predict using empirical rules or tabulated data. Aiming at filling this gap, we propose here guidelines enabling their prediction at a modest computational cost. Rooted in DFT (Density Functional Theory), our protocol features structural optimization and magnetic shielding tensor calculations performed at a global hybrid level using a tailored locally dense basis set. Validation on an experimental data series revealed that while a careful conformational analysis is required in the case of flexible phosphines, the use of the free ligand or another complex as a reference for chemical shifts often allows to solve this drawback. Applicability to various diamagnetic complexes of first-row transition metals is demonstrated, including large systems relevant to contemporary catalysis.

In organometallic and coordination chemistry, 31 P NMR is a key spectroscopy to disclose the structural features of metal complexes having phosphorus ligands. This tool is routinely used for the characterization of stable compounds as well as to evidence short-lived reaction intermediates. 1 Therefore, this technique is particularly useful to elucidate the mechanisms of catalytic reactions and to aid in the development of new transformations 2 promoted by phosphine-ligated metal catalysts. 1a Indeed, 31 P nuclei usually yield easily observable and sharp NMR signals over a wide chemical shift range, providing valuable chemical information due to their high sensitivity even to minor structural changes. 1a Unfortunately, while the chemical shifts of 1 H or 13 C in organic molecules has been traditionally estimated by the use of empirical correlation tables, 3 31 P NMR chemical shifts are notoriously difficult to predict based only on qualitative considerations and tabulated data. 1,[START_REF] Tong | Prediction of 31 P nuclear magnetic resonance chemical shifts for phosphines[END_REF] Indeed, whereas 1 H and 13 C NMR spectroscopies are mainly concerned with contributions due to σ-bonding or to well-defined and recurring πbonded structural elements, the influence of π-bonding in 31 P NMR is larger in magnitude, more frequent and less transferable from one compound to another, thus preventing the development of predictive empirical rules. 1b In the last decades, however, the estimation of 1 H and 13 C NMR shielding constants by using density functional theory (DFT) emerged and imposed as a very reliable and computationally affordable tool to support the structural elucidation of complex organic molecules. 5 By contrast, the use of DFT to predict 31 P NMR chemical shifts has not yet been validated extensively. 6 Van Wüllen first used DFT methods to calculate 31 P chemical shift using both IGLO and GIAO approaches on a series of small molecules. 6b Maryasin and co-workers described a protocol based on the use of a GIAO in conjunction with a rather large basis set (6-311++G(2d,2p)). 6e It seems however that a lower level of theory may be efficient as well: 6d, 6n Latypov and co-workers suggested "smaller" basis sets such as 6-31G(d) may be efficient in most of the cases. 6n Nonspecific and specific solvent effects, 6f such as hydrogen bonding, were also investigated. 6o Indeed, the use of DFT approaches allows to handle large systems such as polyoxometalates. 6m More recently, the case of solid state 31 P NMR has also been examined. 6p As we were writing these lines, Latypov and co-workers reported a study on the applicability of DFT to predict the 31 P NMR chemical shifts of model nickel complexes based on small/medium-sized organophosphorus ligands. [START_REF] Latypov | Quantum Chemical Calculations of 31 P NMR Chemical Shifts in Nickel Complexes: Scope and Limitations[END_REF] The development of an easily handled and reliable routine method would be of great help to assist structural assignments, particularly in the case of complexes featuring multiple phosphine, in the case of complex mixtures of products, or to complete the characterization of non-isolable reaction products and intermediates. 6j, 6l The aim of the present work is to establish and validate a DFT protocol based on the GIAO approach for the prediction of solution (isotropic) 31 P NMR chemical shifts of diamagnetic first-row transition metal complexes. Due to the comparative abundance of data available in the literature for nickel-phosphine complexes, Ni containing complexes were predominantly selected to build a benchmark set suitable to validate our computational approach. Moreover, nickel-based catalytic systems are currently under active development. 8 The full training set of complexes analyzed in this work is reported in Figure 1 (for the complete tables of chemical shifts of complexes and associated references, see the Supporting Information, Tables S1-S2). This series contains a variety of N i ( 0 First, the computational method was established by using a subset of six simple Ni(0) PMe3-ligated complexes (Figure 1, complexes 1-6). Since the chemical shift is critically dependent on the geometry of the complexes, we determined the level of theory needed both for structural optimizations and for magnetic shielding tensor calculations. The resulting methodology was then validated on a larger number of compounds using the full training set depicted in Figure 1. Finally, some concrete examples of application of our protocol to identify intermediates and off-cycles species in the case of "real-life" catalytic cycles are presented in the last part of this paper, including diamagnetic complexes of transition metals other than nickel (Fe(II) and Mn(I)). 0) complexes used to setup the computational protocol (1-6), in black the Ni complexes used to validate it (7-15), and in green some Ni(0) complexes that are discussed more in detail due to tri-aryl phosphine peculiar conformational issues (16-18), in violet Ni(II) complexes which deserve special care due to various degree of paramagnetic character (19-21), in blue complexes of metals other than Ni (22-25). Experimental 31 P NMR chemical shifts in solution are indicated in parentheses. See the Supporting Information, Table S2 for literature references and details of the conditions in which the spectra were recorded. When more than one inequivalent phosphorus nucleus is present, only the one highlighted in boldface was included in the analysis. o-Tol = o-Me-C6H4.

Results and Discussion

We recall that the chemical shift (# ) of a nucleus X is defined as: #

where # and # are the NMR isotropic magnetic shielding constants of a reference nucleus and of the nucleus X respectively and # is the chemical shift of a reference nucleus with respect to a primary reference (for the complete table of calculated isotropic magnetic constants of complexes and free phosphines see the Supporting Information, Tables S3-S4). In the case of 31 P NMR, the zero of the scale (

) is a solution of orthophosphoric acid (H3PO4, 85% in water). Considering the difficulties to model this system, PH3 6b, 6m or PPh3 6e were proposed as "secondary standards" to reference computed chemical shifts. 6e It may also be advantageous to choose different reference compounds for nuclei of the same kind residing in qualitatively different chemical environments. An analogous concept has already been applied to the prediction of 13 C NMR chemical shifts: the use of two different reference molecules, one for sp 3 -hybridized and one for sp and sp 2 hybridized carbons allows indeed to reduce the error. [START_REF] Sarotti | A Multi-standard Approach for GIAO 13 C NMR Calculations[END_REF] In our case, we are particularly interested in the prediction of the change in chemical shift associated to the formation of the metal-phosphorus bond. This variation can be defined as # where # and # are the chemical shifts of the phosphorus nucleus in the metal-bound and free ligand respectively and # and # are the corresponding isotropic shielding constants. Indeed, the chemical shift of the free phosphine ( #) is normally known and excess free ligand may be present in the analyzed solution and used as a secondary reference standard. Thus, the 31 P chemical shift of the phosphine in the metal complex can be rewritten as # .

Setup and Validation of the computational approach

Experimental studies on series of structurally related complexes, highlighted that the metal-phosphorus bond length [START_REF] Hunter | Correlations among 31 P NMR Coordination Chemical Shifts, Ru-P Bond Distances, and Enthalpies of Reaction in Cp′Ru(PR3)2Cl Complexes (Cp′ = η 5 -C5H5, η 5 -C5Me5[END_REF] and the cone angle of the phosphine [START_REF] Carturan | Influence of the phosphine ligands in the 1 H and 31 P NMR behaviour of [(η 3 -allyl)Py(phophine)Cl] complexes[END_REF] have a critical effect on the chemical shift. 12 As a consequence, obtaining reliable geometries is a pre-requisite to predict 31 P NMR chemical shift. The Ni-P bond length and the Ni-P-C angle were chosen as indicators to judge the convergence of the computational protocol for structural predictions focusing first on complex 1 (Figure 1).In order to setup a cost-effective method, we decided to fix the exchange correlation functional to the popular global hybrid B3LYP 13 and to screen different basis set combinations so as to obtain adequately accurate convergence of the geometric parameters at a modest computational cost. Solvent effects were accounted for with the PCM model, using the parameters of the solvent used for experimental measurements (for both optimization and NMR calculations; see the Supporting Information, Table S2 for literature references and details of the conditions (solvent) in which the spectra were recorded). First, we analyzed the effect of the basis set on the metal atom.

We considered the Los-Alamos family of basis and pseudo potential for the Ni atom and tested the LANL2DZ and LAN-L2TZ basis sets. [26] Starting with the triple-zeta Los-Alamos basis set (LANL2TZ) used with the associated pseudo-poten-

δ X δ X = σ ref -σ X + δ ref σ ref σ X δ ref δ ref = 0 δ LP -δ FP = σ FP -σ LP δ LP δ FP σ LP σ FP δ FP δ LP = σ FP -σ LP + δ FP
tial (LANL2), addition of diffuse (LANL2TZ+) or f-polarization functions (LANL2TZ(f)) lead to less than 1 mÅ change in the Ni-P bond length (see the Supporting Information, Table S5). By contrast, using the smaller LANL2DZ basis set lead to a variation of more than 20 mÅ. To estimate the impact of bond length on the isotropic shielding constant, a scan of σ along the Ni-P distance in complex 1 was performed. A modification of 1 mÅ induced a change of less than 0.1 ppm (see the Supporting Information, Figure S1), which validates the choice of LANL2TZ/LANL2 to describe Ni for geometry optimization. Similarly, the effect the basis set on P and O on the optimized Ni-P bond length and P-Ni-P angle was studied (see the Supporting Information, Table S6). In this case split-valence Pople basis sets (double or triple zeta) were evaluated. Starting with 6-31G, addition of diffuse functions (6-31+G) resulted in a noticeable change of the bond length (-61 mÅ). By contrast, the addition of d-polarization functions (6-31+G(d)) had only a negligible impact. Changing from double to triple zeta basis (6311+G(d)) does not significantly affect the results. Conse-quently, the B3LYP functional in conjunction with LAN2TZ/ LANL2 on Ni, 6-31+G (d) on P, O and other atoms, and 6-31G (d) on C and H have been retained as computational protocol for structural optimization. Hereafter this composite basis set will be referred to as BS1. All optimized structures are reported in the Supporting Information, §4 Cartesian coordinates and absolute energies).

Another caveat with the structural optimization of phosphine and phosphine-metal complexes is the possibility of populating several conformations at room temperature. Indeed, a conformational analysis and a subsequent Boltzmann average was necessary for the prediction of the chemical shift of flexible phosphites. 6e For instance, if we consider the very simple structure of complex 2 (Figure 1), there could be free rotation around the Ni-P bond and the cone angle of the phosphine may also vary. To estimate the influence of these effects on the isotropic shielding tensor, a scan of σ along the C-P-Ni-P dihedral was performed (Figure 2a). A 60° rotation of trimethylphosphine corresponds to a variation of σ of more than 6 ppm. 14 This process is predicted to be nearly barrierless and thus may have a significant impact on the chemical shift. Using a Boltzmann distribution and sampling conformations with a 5° rotation step, the chemical shift is indeed shifted downfield of 1.5 ppm compared to the chemical shift computed for minimal energy configuration.

The effect of the cone angle in complex 2 on the chemical shift was also investigated. As described in the literature [START_REF] Carturan | Influence of the phosphine ligands in the 1 H and 31 P NMR behaviour of [(η 3 -allyl)Py(phophine)Cl] complexes[END_REF] the effect of the cone angle on the isotropic shielding tensor is relevant (more than 60 ppm for a 15° distortion, see the Supporting Information Figure S2). If we consider only the structures accessible at room temperature (Figure 2b, corresponding to less than 10 kJ mol -1 for a variation of 5°) the variation of the isotropic shielding tensor in this area is more than 25 ppm. Using a Boltzmann-averaged sampling, the predicted chemical shift changed by nearly 1 ppm compared to the most stable staggered conformation. An analogous calculation on free PMe3 (see the Supporting Information, Figure S3) revealed a very similar behavior and the contribution to the chemical shift due to a variation of the cone angle is of the same order of magnitude (0.9 ppm). This may be explained by considering that in complex 2 the steric hindrance due to the interaction with the CO ligands is very weak and thus the possible deformations of the ligated phosphine are similar to those of the free ligand.

Based on these observations, we can conclude that it is advantageous to use the corresponding free phosphine as a reference for the calculation of the 31 P NMR chemical shifts of complexes and focus on the difference of chemical shift between the free and ligated phosphine. When the accessible conformational spaces are similar for the free and ligated phosphine, no conformational analysis would likely be needed due to transferability of dynamic conformational effects between the free and the metal-bound ligand. On the contrary, the rotation of the phosphine around the metal-P axis (a degree of freedom not relevant for the free phosphine) should be checked if thermally accessible.

The choice of the basis set for the calculation of the NMR shielding tensor was then examined. The use of unbalanced, "locally dense" basis sets for the cost-effective calculation of 13 C NMR chemical shifts in large organic molecules at the Hartree-Fock level was first proposed by Chesnut and Moore. 15 In this approach, high quality large basis set are employed only on the atom for which the NMR chemical shift is to be calculated, while computationally cheaper basis sets are used to describe all the other atoms. 16 Its applicability in conjunction with post-HF methods has been investigated later, [START_REF] Chesnut | The use of locally dense basis sets in correlated NMR chemical shielding calculations[END_REF] and this subject has been reexamined in a systematic fashion for the case of peptides recently. 18 In our case, the "locally dense" basis set concept would allow significant saving of computa- We investigated the suitability of this approach for metal complexes at DFT level. Thus, we kept the basis set employed for geometry optimization for all atoms but the P atom of interest and we investigated the effect of the basis on the latter. In this respect, we considered the family of the polarization consistent basis sets with segmented contraction pcsSeg-x (x = 0-3), which was designed by Jensen to describe the core properties of nuclei. [START_REF] Jensen | Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods[END_REF] This basis set family proved to be particularly efficient for NMR calculations. [START_REF] Jensen | Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods[END_REF] ! Figure 3. Error on the chemical shift calculated for complex 3 (with respect to experimental values), using polarization consistent basis sets with segmented contraction on the P atom of interest (pcSseg-x, x = 0, 1, 2 and 3) with the GIAO method, and the B3LYP functional. The following basis set was used for the other nuclei: LANL2TZ-LANL2 (Ni), 6-31G(d) (C and H), 6-31+G(d) (other P atoms).

Using the B3LYP functional in combination with the GIAO approach, the absolute error on the chemical shift of complex 3 (Figure 3) converged to less than 1 ppm using the pcsSeg-3 basis set on the P nucleus of interest (Figure 3, the 6-31+G (d) basis set was employed for the other three P centers). The B3-LYP functional was selected on the grounds of previous studies revealing its efficiency in calculation of NMR shielding tensors. 6d Differently from what demonstrated for related systems, [START_REF] Zhang | OPBE: A promising density functional for the calculation of nuclear shielding constants[END_REF] the PBE0 functional did not significantly outperform B3LYP (see the Supporting Information, Figure S4). The effect of the basis set on the other atoms was also investigated. Addition of diffuse functions to describe C and H significantly improves the results while addition of d-polarization functions is detrimental (see the Supporting Information, Table S8). Going from a double zeta basis set to a triple zeta did not lead to any significant improvement (see the Supporting Information, Tables S8) while sizably increasing the computational cost. Concerning P atoms other than that of interest, both d-polarization and diffuse functions are needed (see the Supporting Information, Table S9). The following composite basis set was selected and it will be referred as BS2 in the following: pcsSeg-3 on the P nucleus of interest, 6-31G (d) on C and H atoms, 6-31+G (d) on any other atoms (including P centers other than that of interest). Finally, we assessed the magnitude of solvent effects modelled using polarization continuum model (PCM, refer to Computational Details). The 31 P chemical shift of complex 3 varied of about 2 ppm passing from benzene to acetonitrile, which is in agreement with previous reports 6e, 6f (see the Supporting Information, Figure S5) and of comparable magnitude with the accuracy of our method (vide infra). We decided to keep PCM solvation for the rest of this study, as it does not significantly increase the computational cost.

The optimized protocol was then tested on complexes 1 to 6 (Table 1) and proved to be efficient in the prediction of chemical shifts. Indeed, based on geometric optimization at the B3-LYP/BS1 level and subsequent NMR calculation at the GIAO-B3LYP/BS2 level (vide supra the definition of BS1 and BS2 composite basis sets), it was possible to predict the chemical shift of these species with an average unsigned error of 1.6 ppm.

Scope and limitations

The efficiency of the method optimized on a small subset of Ni(0), PMe3-ligated complexes was then further tested on a larger set (Table 2,Figure Varying the nature of the phosphine did not drastically influence the accuracy of the prediction. However, in the case of highly flexible phosphines such of PEt3 (complex 9) or PF2CH3 (complex 10) the absolute error increased, thus suggesting the importance of conformational analysis and dynamics in these instances. ppm (Figure 5). Experimentally, the spectra of this complex features two signals at 20.7 and 50.7 ppm in 31 P NMR, matching with the second structure. This conclusion is in agreement with the computed thermodynamics, since the second geometry is predicted to be favored by 24.6 kJ mol -1 . The case of aryl phosphine was studied in detail (vide infra) due to their relevance as ligands in catalysis. The prediction of the 31 P chemical shift of complexes of aromatic phosphines is especially difficult due to their conformational freedom. 6e, 21 The case of P(o-Tol)3 has been analyzed in depth (complex 17).

Rotation of the o-Tol groups around the P-C bond can be considered as free at ambient temperature with an energy barrier estimated at only 7 kJ mol -1 (Figure S7). Two minimal energy structures are obtained, corresponding to the conformation with the methyl group outside ("out") or inside ("in") the cone formed by the aryl groups (Figure 4). The "in" conformer is thermodynamically favored by 16 kJ mol -1 , consequently only this conformer was considered. However, the metal-bound P(o-Tol)3 complex 17 has two nearly isoergonic conformations of the ligand, with the 'in' one slightly favored by 5 kJ mol -1 . In this case, both conformations have to be considered to obtain a reasonable estimation of the chemical shift at room temperature (Table 3). Both conformer are separated by more than 20 ppm due to the variation of the cone angle. The same reasoning was successfully applied to the case of (o-iPr-C6H4)3P in complex 18 (Table 3). Changing the oxidation state of the metal from Ni(0) to Ni(II) (complex 19-21) can affect the accuracy of the approach as evidenced by the values of Table 4. While the chemical shift of complex 19 is described correctly the error become significant for complexes 20 and 21. This failure can be directly related with the energy gap between singlet and triplet state (∆G(S!T), Table 4): the lower the triplet energy, the more the error. As already noticed by Latipov and co-workers, residual paramagnetic behavior hampers accurate chemical shift prediction by DFT calculations. [START_REF] Latypov | Quantum Chemical Calculations of 31 P NMR Chemical Shifts in Nickel Complexes: Scope and Limitations[END_REF] Finally, our protocol could also be applied to complexes based on first-row transition metal other than Ni (such as Cu(I) and Cr(0), complexes 22-25), but it substantially failed at predicting the 31 P NMR chemical shifts of complexes of second and third row metals (here exemplified by Mo(0), W(0), complexes 24-25), likely due to the lack of appropriate treatment of relativistic effects in the chosen model. 6f, 6m The data from Tables 1234were used to plot calculated vs experiment chemical shift (Figure 6). Overall our method performs well with a regression coefficient of R 2 = 0.99, the slope is very close to one (slope = 0.96). Consequently, no additional affine correction is needed, differently from previous reports. 6n The average unsigned error was 4.2 ppm, it compares favorably with previously reported methods. (in red, test set). Geometries optimized at the B3LYP/BS1 level of theory. 31 P NMR shielding constants were evaluated with the GIAO method at the B3LYP/BS2 level of theory.

Application to the structural elucidation of complex catalytic system

To demonstrate the interest of our approach for practical applications, we report here some examples showing how it could assist in the interpretation of experimental 31 P NMR spectra.

First, we attempted to discriminate between a mixture of structurally related Ni(II) complexes differentiated by their anionic ligands (Cl vs Br vs OH), organization (mono vs dinuclear form) and, finally, configuration (trans vs cis). The complexes selected for this test (Figure 7, Ar = 2-Me-4-F-C6H3) were described in a previous work on the mechanism of the nickelcatalyzed Suzuki-Miyaura cross-coupling, in which the role of these species in the catalytic cycle was elucidated. [START_REF] Payard | Taming Nickel-Catalyzed Suzuki-Miyaura Coupling: A Mechanistic Focus on Boron-to-Nickel Transmetalation[END_REF] More in details, the interaction of trans-[NiArCl(PPh3)2] (complex 27) with Br -and OH -has been described. The variation of the . This trend is in agreement with the experimental observations. [START_REF] Payard | Taming Nickel-Catalyzed Suzuki-Miyaura Coupling: A Mechanistic Focus on Boron-to-Nickel Transmetalation[END_REF] Lastly, the calculations allow to put forward an assignment of the signals of all isomers, which could not be deduced from experimental data (Figure 7). # Figure 7. ArNi(II) complexes potentially involved in the catalytic cycle of a Ni-catalyzed Suzuki coupling [START_REF] Passera | Mn(I) and Fe(II)/PN(H)P Catalysts for the Hydrogenation of Ketones: A Comparison by Experiment and Calculation[END_REF] Experimental [START_REF] Passera | Mn(I) and Fe(II)/PN(H)P Catalysts for the Hydrogenation of Ketones: A Comparison by Experiment and Calculation[END_REF] (blue) and calculated (red, B3LYP(GIAO)/BS2//B3LYP/BS1 level of theory) chemical shifts are reported. Complex 27 was used as a secondary reference (δ = σ28σ + δ28).

To highlight the applicability of our computational protocol to assist the interpretation of 31 P NMR of complexes of first-row transition metal relevant to catalysis, additional systems based on metals other than nickel were examined. The Mn(I) complexes 30a-b (Figure 8) featuring a P-stereogenic PN(H)P ligand are useful as catalyst for the asymmetric hydrogenation of prochiral ketones. [START_REF] Passera | Mn(I) and Fe(II)/PN(H)P Catalysts for the Hydrogenation of Ketones: A Comparison by Experiment and Calculation[END_REF] In the context of an investigation of the mechanism of this reaction, several potentials intermediates of the catalytic cycle were prepared and studied by 31 P NMR spectroscopy. More in details, treatment of 30 with KH gave a mixture of the amido complex 31, syn hydride 32 and anti hydride 33. Compound 32a reduced acetophenone to give the diastereoisomeric alcoxo complexes (S)-34a and (R)-34a. 31 P NMR chemical shifts of these complexes were calculated using our theoretical protocol (Figure 8) taking complexes 30a (for R = Cy) and 30b (for R = tBu) as references to anchor the computed chemical shifts, to avoid a burdensome conformational analysis of this flexible ligand. The mean unsigned error of the 11 calculated values was 3.4 ppm, which is in line with the accuracy observed in our benchmark set of Ni complexes. Remarkably higher accuracy was obtained series of complexes with R = tBu (31-32b) than for analogues with R = Cy (31-32a), most likely due to the reduced conformational freedom of the tBu residue compared to Cy. A last example was drawn from the field of iron catalysis, based on a recent report by Mezzetti and coworkers. [START_REF] De Luca | Asymmetric Transfer Hydrogenation with a Bifunctional Iron(II) Hydride: Experiment Meets Computation[END_REF] In the course of the study of the mechanism of the asymmetric transfer hydrogenation of ketones catalyzed by the iron(II) complex [Fe(CNCEt3)2(L)](BF4)2 (L is a chiral N2P2 macrocycle), the reaction of the hydride complex 35 [FeH(CNEt3)( 1)](BF4) with acetophenone was studied (Figure 8). [START_REF] De Luca | Asymmetric Transfer Hydrogenation with a Bifunctional Iron(II) Hydride: Experiment Meets Computation[END_REF] A complicated mixture of products containing the diastereoisomeric alkoxo complexes (S)-36 and (R)-36, and bromo complex 37 (from Br -present in the reaction medium) was obtained. Application of our procedure to the prediction of the 31 P NMR chemical shifts of these species relevant to catalysis (with 35 as the reference compound, 7 computed values, Figure 9) gave a mean unsigned error of 3.8 ppm, which is a performance similar to that observed for the systems discussed previously. This observation suggests that stereochemical analysis of the remote Et3CNC ligand, which has high conformational freedom, is not needed (i.e. this ligand is sufficiently far from the P, so that the detailed description of its conformations is not required to estimate the 31 P NMR chemical shifts). Unfortunately, the difference in the chemical shifts of diasteromers (R)-and (S)-36 (3.1 ppm in the best case) is so small that the error of the method would not allow to use the calculations to support assignment of the spectra. This example, however, highlights that our computational protocol can be applied to rather large systems (complexes (R)-and (S)-36 have 58 non-hydrogen atoms) due to its moderate computational cost, and that no significant loss of accuracy is observed passing from smaller model compounds to these larger systems encountered in modern catalysis. The data presented of Figures 789were added to the calculated versus experimental data plot (red, Figure 6), no difference in behavior was observed with the training set described in the first part of this paper (blue), the average unsigned error was 4.1 ppm, the slope closed to one (1.02) and the regression coefficient is only slightly lower (0.94).

# Figure 9. Fe(II) complexes of a chiral P2N2 macrocyclic ligand involved in the asymmetric transfer hydrogenation of prochiral ketones. [START_REF] De Luca | Asymmetric Transfer Hydrogenation with a Bifunctional Iron(II) Hydride: Experiment Meets Computation[END_REF] Experimental 24 (blue) and calculated (red, B3LYP(GIAO)/BS2//B3LYP/BS1 level of theory) chemical shifts are reported. Complex 35 was used as a secondary reference (δ = σ28σ + δ28).

Conclusions

Guidelines enabling the theoretical prediction of 31 P NMR chemical shifts of first row transition metal complexes at a modest computational cost have been established. First, adequately accurate geometries are to be obtained. At this end, our protocol employs a global hybrid functional (B3LYP) with a composite basis set (BS1). BS1 comprises a triple zeta basis set with an effective core potential (LAN2TZ/LANL2) on the metal and a double zeta basis set on other atoms (6-31G(d) for C and H, 6-31G+(d) for other atoms). Solvent effects are included using an implicit solvation model (PCM). These resulting optimized structures are submitted to NMR shielding tensor calculation using the GIAO approach with the B3LYP functional and a tailored composite basis set (BS2). BS2 is a composite "locally dense" basis set comprising a triple-zeta quality segmented polarization-consistent basis set (pcs-Seg3) on the P center of interest and double zeta basis set on other atoms (6-31G (d) for C and H, 6-31G+(d) for other atoms. Careful analysis of conformational and dynamic effects on the chemical shift is crucial to obtain reliable results. In the simplest cases, these effects are transferable from the free phosphine to the complex. Consequently, using the free phosphine as a reference to anchor the calculated chemical shifts allows to reduce errors due to systematic error cancelation. In selected cases, however, sampling of conformational space and Boltzmann averaging is recommended. Overall, this procedure allows estimating 31 P NMR chemical shifts of diamagnetic complexes featuring different phosphines and first-row metals with errors typically below 3 ppm. This optimized protocol has been applied to catalytically relevant systems and proved particularly robust for the estimation of chemical shift differences of structurally related complexes. In perspective, calculated 31 P NMR chemical shifts can be used as a tool to assist spectral assignment of stable complexes and transient catalytic intermediates.

Computational Details

All DFT calculations were performed with the Gaussian09 program (Rev. A.02). [START_REF] Frisch | Gaussian[END_REF] The structures of all minima were fully optimized using the B3LYP functional 13 without any symmetry constraint. Basis sets used for structural optimization (BS1: LAN2TZ/LANL2) 26 on Ni, 6-31G (d) on C and H and 6-31+G(d) on all other atoms) and NMR calculation (BS2: pcsSeg-3 [START_REF] Jensen | Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods[END_REF] on the P nucleus of interest, 6-31G(d) on C and H atoms, 6-31+G (d) on any other atoms (including P centers other than that of interest)) are detailed in the text (section "Setup and Validation of the computational approach). To simulate experimental conditions, bulk solvent effects were accounted for by using an implicit solvation model (IEF-PCM) as implemented in Gaussian, for both geometry optimization and shielding constant calculations. 27 The default spheres radii, static and optical dielectric constants for each solvent were used. All stationary points were characterized as minima by analytical frequency calculations.

ASSOCIATED CONTENT

Supporting Information

A listing of the contents of each file supplied as Supporting Information should be included. For instructions on what should be included in the Supporting Information as well as how to prepare this material for publication, refer to the journal's Instructions for Authors.

The Supporting Information is available free of charge on the

  ) c o m p l e x e s ( 1 -1 8 ) a n d N i ( I I ) c o m p l e x e s (complexes 19-21). It is complemented by a few complexes of other transition metals (Cr(0), Cu(I), Mo(0) and W(0), complexes 22-25).

#Figure 1 .

 1 Figure 1. Structures and experimental 31 P isotropic chemical shifts of the training set of complexes used in this study. In red are reported the Ni(0) complexes used to setup the computational protocol (1-6), in black the Ni complexes used to validate it(7-15), and in green some Ni(0) complexes that are discussed more in detail due to tri-aryl phosphine peculiar conformational issues(16-18), in violet Ni(II) complexes which deserve special care due to various degree of paramagnetic character(19-21), in blue complexes of metals other than Ni(22- 25). Experimental 31 P NMR chemical shifts in solution are indicated in parentheses. See the Supporting Information, TableS2for literature references and details of the conditions in which the spectra were recorded. When more than one inequivalent phosphorus nucleus is present, only the one highlighted in boldface was included in the analysis. o-Tol = o-Me-C6H4.

  a) Energy scan (ΔE) for rotation of PMe3 (θ = CNi-PC dihedral angle) around the Ni-P bond axes in complex 2 (right scale, black curve) and the associated evolution of 31 P NMR isotropic shielding tensor (left scale, blue curve, Δσ = σ(θ) -σ(θmin), with θmin the minimal energy angle). b) Energy scan (ΔE) for change in the PMe3 cone angle (ϕ) in complex 2 (right scale, black curve) and the associated evolution of 31 P NMR isotropic shielding tensor (left scale, blue curve, Δσ = σ(ϕ) -σ(ϕmin), with ϕmin the minimal energy angle). 31 P NMR shielding constants were evaluated at the B3LYP(GIAO)/BS2 level of theory.

  especially for large metal complexes featuring multiple phosphorus ligands.

  1) including various phosphine-type ligands (complexes 7-18, PF3, PCl3, PEt3, PF2CF3, PF(CF3)2, PCy3, PPh3, P(o-Tol), P(o-iPr-C6H4)3, bidentate phosphines), Ni(II) complexes (complexes 19-21) and finally complexes of other metals (complexes 22-25, Cu(I), Cr(0), Mo(0), W(0)).

  6e

# Figure 5 .

 5 Conformation of complex 15 [Ni(η 2 -CO2)(PCy3)2], calculated 31 P chemical shifts are indicated (Ni atom in yellow, P atoms in green, O atoms in red, C atoms in white, H atoms omitted for clarity). PCy3 was used as a secondary reference (δ = σPCy3 σ + δPCy3).Two different conformations can be optimized for [Ni(η 2 -CO2) (PCy3)2] (complex 15, Figure5), a complex relevant for activation of CO2. The first displays orthogonal PNiP and OCO planes while the latter shows coplanar PNiP and OCO. In the latter case, the two phosphorus centers are not equivalent. While this structure is predicted to possess two resonances at 11.1 and 62.5 ppm, the latter would have a single signal at 23

!

Figure 4 .

 4 Figure 4. Conformations of P(o-Tol)3 (P atoms in green, C atoms in white, H atoms omitted for clarity).

7

 7

# Figure 6 .

 6 Calculated versus experimental chemical shift of complexes 1-25 (in blue, training set, complexes 21, 24 and 25 were omitted for above discussed reasons, namely residual paramagneticity and relativistic effects) and complexes 26-37

  31 P NMR chemical shift upon substitution of the chlorido ligand by a bromido to give trans-[NiArBr(PPh3)2] (complex 26, experimental ∆δ = +0.8 ppm) is well reproduced by calculations (computed ∆δ = +0.6 ppm). Upon addition of OH -to 26, four new signals at higher field appeared. Substitution of Br -by OH -may potentially yield either the monomeric complex trans-[NiAr(OH)(PPh3)2] (28) or dimeric µ-OH bridged complexes [NiAr(µ-OH)(PPh3)]2 (29a-d). Calculations allow to clearly decide between these two hypotheses: a negative variation of -8 ppm is expected for the monomer versus a positive shift of +15 ppm for dimers. Finally, four isomers of [NiAr(µ-OH)(PPh3)]2 can be formed (see Figure 6): phosphine ligands may be placed trans (29a-b) or cis (29c-d) around the dinuclear core and the ortho methyl groups of Ar can be either anti (29b,d) or syn (29a,c) with respect to the Ni(µ-OH)2Ni plane. Complexes featuring the phosphine arranged in a trans fashion are predicted to resonate at lower field (37.4 and 38.0 ppm) than complexes with phosphine arranged in a cis conformation (34.2 and 36.2 ppm)

Figure 8 .

 8 Mn(I) complexes of a P-chiral PNP ligand involved in the asymmetric reduction of prochiral ketones.[START_REF] Passera | Mn(I) and Fe(II)/PN(H)P Catalysts for the Hydrogenation of Ketones: A Comparison by Experiment and Calculation[END_REF] Experimental 23 (blue) and calculated (red) chemical shifts are reported. Complex 30a and 30b were used respectively as a secondary reference for complexes 31-34a and 31-32b.

Table 1 .

 1 Experimental and calculated 31 P chemical shift (ppm) of complexes 1-6, PMe3 was used as a secondary reference (δ = σPMe3 -σ + δPMe3).

	Complex 1	2	3	4	5	6
	Exp.	-18.2 -19.9	-22.2	-26.5 -17.9	-22.2
	Calc.	-18.3 -17.4 a	-23.5	-25.2 -15.1	-20.7
	A b s . error	0.1	2.5	1.3	1.3	2.8	1.5
	a Before Boltzmann correction related to free rotation (vide
	supra, Figure 2a): -16.3 ppm. With correction corresponding to
	free rotation (no Boltzmann-weighed averaging): -18.0 ppm.
	With Boltzmann-averaged correction -17.4 ppm.

Table 2 .

 2 Experimental and calculated 31 P chemical shift of complexes 7-15. The associated free phosphine (L) was used as a secondary reference (δ = σLσ + δL).

	Complex	Phosphine (L)	Calculated δ (ppm)	Experimental δ (ppm)	Abs. error (ppm)
	7	PF3	136.5	136.0	0.5
	8	PCl3	176.3	182.5	6.2
	9	PEt3	-8.6	-3.0	5.6
	10	PF2CF3	170.5	184.8	14.3
	11	PF(CF3)2	155.1	156.7	1.6
	12	PEt3	19.7	19.4	0.3
	13	tBu2PCH2-CH2PtBu2	86.6	92.7	6.6
	14	tBu2PCH2-CH2PtBu2	89.2	95.6	6.4
	15	PCy3	11.1	20.7	9.6
			62.5	50.7	11.8

Table 3 .

 3 Experimental and calculated 31 P chemical shift of complexes 16-18 for different conformations and average values. The associated free phosphine (L) in its most stable conformation ("in") was used as a secondary reference (δ = σLσ + δL).

		Δ f G i n / out (kJ mol 1 )	er (ppm) Calc. δ -' o u t ' conform	er (ppm) Calc. δ -" i n ' conform	Average calc. δ	Exp. shift	A b s . error
	16	-	-	-	47.7	42.9	4.8
	17	5.5	30.65	18.36	24.5	29.4	4.9
	18	5.1	46.15	21.83	33.9	31.9	1.9

Table 4 .

 4 Experimental and calculated 31 P chemical shift of complexes 19-21. The associated free phosphine (L) was used as a secondary reference (δ = σLσ + δL).

	Complex Calculated δ (ppm)	Experimental δ (ppm)	Abs. error	ΔG(S!T )
				(ppm)	(kcal mol -
					1 )
	19	39.3	35.6	3.7	+17.7
	20	70.6	70.3	0.3	+6.7
		97.7	78.3	19.4
	21	42.8	85.5	42.7	+1.9

Table 4 .

 4 Experimental and calculated 31P chemical shift of complexes 22-25. PMe3 was used as a secondary reference ( δ = σPMe3σ + δPMe3).

	Comple	Metal center	Calculate	Experiment	Abs.
	x		d	δ	al δ (ppm)	error
			(ppm)			(ppm)
	22	Cu(I)	-43.4		-47.5	4.1
	23	Cr(0)	13.8		6.5	7.3
	24	Mo(0)	-2.4		-17.3	14.9

ACKNOWLEDGMENT

We thank CNRS, ENS, ENSCP, and the ANR program CD2I (project CuFeCCBond) for financial support. P.A. Payard is grateful to ENS Paris Saclay for a PhD grant.

AUTHOR INFORMATION

Corresponding Author *pierre-adrien.payard@ens.fr; luca.perego@ens.fr; ilaria.ciofini@chimie-paris.psl.eu Author Contributions ‡These authors contributed equally.