Essays doi.org/10.1002/ejic.202100330 Computational Organometallic Catalysis: Where We Are, Where We Are Going Agustí Lledós*[a] Dedicated to Professor Joan Bertran, who opened the doors to the quantum world for me, on the occasion of his 90th birthday This essay gives my personal perspective of the current stage of closed-shell species. However, in most cases they are still used computational methods applied to modeling organometallic to “predict the past”, providing after-the-fact explanations and catalysis, as well as the new directions the field is taking. The missing out the full potential of contemporary simulation first part of the essay deals with what I consider the state-of- techniques. Simulations are mature enough to be incorporated the-art to build up energy profiles, regarding both chemical at the design stage and to guide the experimental exploration. and computational models. With a proper choice of the The new directions the field is taking, incorporating automated chemical model and computational methods, quantum me- exploration methods and combined with extensive data chanical calculations are nowadays able to provide accurate analysis and machine learning algorithms, approach the holy energy profiles of organometallic reactions in solution involving grail of catalyst discovering. 1. Introduction The tremendous increase in computing power, together with the development of user-friendly software and quantum Over almost forty years of theoretical research on organo- chemistry methods, specially DFT, makes it possible to metallic reactivity, I have witnessed the transformation of accurately compute energy profiles of complex molecular quantum mechanical (QM) calculations from something un- systems. However, there is not always a clear understanding of common, exclusively on specialist’s hands and suitable just for the sources of error and limitations inherent to these qualitative analysis, to becoming a technique fully integrated in techniques, nor about their quality control. Nowadays, calcu- the chemist’s toolbox and able to give quantitative measure- lations are usually employed in reactivity studies to measure a ments. In the last years, density functional (DFT) calculations physical quantity (the potential energy of an ensemble of nuclei have become the most common tool to study reaction and electrons at fixed nuclei positions). While for experimental mechanisms. A basic search in the Web of Science combining measurements such as elemental analysis, NMR chemical shifts the words “DFT” and “European Journal of Inorganic Chemistry” and X-ray diffraction analysis there is a general consensus about evidences that, in the last years, more than a ten percent of the the minimum standards required, the same does not apply to articles published in this journal include such calculations. A calculations. Even more demanding is the connection of these few recent reviews witness the importance that quantum measurements with the proposal of a reaction mechanism. chemical reaction modeling has gained in the field of homoge- Calculations must be performed using state-of-the-art computa- neous organometallic catalysis.[1–4] As an “old specialist” I am tional approaches. Of course, the standard to be chosen will pleased to cope with the “democratization” of these techniques, depend on the nature of the system to be computed and nowadays within the reach of every chemist. But this fact also evolves with time, very quickly in the case of computational raises some philosophical questions, such as what the role of techniques. The state-of-the-art can be defined as the best theoreticians is today or even if we are no longer needed. In possible calculation at any given time. The first part of this this essay I will give my personal view about the current essay will be devoted to what I consider the current state-of- situation of theoretical calculations applied to the reactivity of the-art methods. Using them, reliable energy profiles can be molecular transition metal systems, particularly related to constructed, but these do not answer all the mechanistic homogeneous catalysis, as well as the new developments I questions that chemists ask themselves: what is occurring in my foresee for the near future. flask and why, what will happen if…and even more, how can I improve one reaction or discover new reactions? These [a] Prof. Dr. A. Lledós questions will be the subject of the second part of this essay. Departament de Química The scope and challenges,[5–10] as well as pitfalls[11] of computa- Universitat Autònoma de Barcelona Campus UAB, 08193 Cerdanyola del Vallès, Catalonia, Spain tional methods for mechanistic studies in homogeneous E-mail: agusti@klingon.uab.es catalysis have also been addressed in some recent articles https://compbionanocat.wordpress.com/ collecting a number of examples, pointing out the current Part of the “RSEQ-GEQO Prize Winners” Special Collection. interest on the topic. © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Eur. J. Inorg. Chem. 2021, 2547 – 2555 2547 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 2. Building Up Energy Profiles: The Choice of the Model Until recently, computational tools were mainly used to explain chemical reactions after experimental results were obtained. The usual goal was to substantiate how the main product was formed and to understand the factors controlling the reaction. In order to do that, a sequence of intermediates and transition states (TS) along the pathway from reactants to products should be located and characterized. The connection of TSs with back and forward intermediates must be carefully checked. Ideally, this is done by two intrinsic reaction (IRC) calculations following both reaction path directions. In flat energy surfaces IRC often Scheme 1. Minimal and extended chemical models to perform Quantum- fails. In that case, displacement of the TS structure in both Mechanical based calculations on organometallic reactions. directions following the eigenvector corresponding to the imaginary frequency, and subsequent optimization of the resulting structures is recommended. QM calculations, com- bined with optimization techniques, provide absolute energies chosen to carry out the calculations. This is the computational and geometries of chemically relevant structures in the model (formerly called the level of calculation). It should never potential energy surface of the system (PES). The relative be overlooked that all QM methods applied to molecular energies of all these species, usually represented as an energy systems rely on approximations and it is safe to have in mind profile, describe the energy landscape for the reaction. Analysis some estimation about the error that the chosen approximation of their geometries discloses the bond-forming and bond- can introduce in the results. breaking events taking place along the reaction pathway and The choice of both models is the first and most critical allows the proposal of a mechanism for such reaction. In this decision to be taken when starting a computational project, approach the exploration of the PES to locate stationary and it will have a major impact on the quality of results. structures (minima and TSs) is biased by the previous knowl- Traditionally, particularly from the theoretical side, much more edge and the chemical intuition of the person performing the attention has been devoted to the computational model than calculations. to the chemical one, but a very accurate calculation in an One of the things that makes computational methods incomplete chemical model can yield completely wrong results different from experimental techniques is that they do not work regarding reaction mechanisms. with samples, but with models. Two types of models operate in a QM calculation: on one side, the chemical complexity found in a reaction flask (millions of molecules of different types) 3. The Chemical Model should be reduced and adapted to something that a computer can handle. This is known as the chemical model. As This article focuses on organometallic chemistry in solution. calculations provide an atomistic description of the reaction, That is, reactions performed in a solvent in a range of the most common way to do this is to include, in QM temperatures typically between 0–100 degrees Celsius. We calculations, one molecule of the organometallic reagent and want to simulate such a chemical reaction in a computer (i.e, in one molecule for each reactant. In many cases, this minimal silico). Beforehand, it is necessary to have a clear perception of model has to be extended in order to obtain reliable computa- the experimental conditions under which the reaction takes tional results (Scheme 1). On the other side, among the panoply place. As commented above, the most common chemical of the available computational methods (for DFT calculations model includes one molecule of the organometallic reagent, hundreds of different functionals and basis sets) one must be nowadays with the actual ligands without further simplification, and one molecule for each reagent. Even when this model is reliable, the conformational flexibility of most of the bulky Agustí Lledós was born in Barcelona in 1955. ligands employed in catalysis causes additional complications He obtained a Ph.D. in Chemistry from the and can be a source of error.[12] An example of the importance Universitat Autònoma de Barcelona (UAB) in 1984. After a postdoctoral stay at the Uni- of a wide exploration of the conformational space when versité de Paris-Sud, in 1987 he started computing energy profiles can be found in the recent computa- computational research on organometallic tional analysis of the Buchwald-Hartwig amination of aryl catalysis at UAB. He has been Full Professor of halides using Pd-Josiphos complex.[13] The Gibbs energy differ- Physical Chemistry at UAB since 1994. In 2018 ence (DFT calculations with M06 functional) between two he was awarded with the Gold Medal of the conformations of the same intermediate (15.3 kcal · mol 1) is Organometallic Chemistry Group of the Span- ish Royal Society of Chemistry (RSEQ). much larger than the barrier of the step in which this intermediate is involved (12.3 kcal · mol 1) (Figure 1). Recently Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2548 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 Figure 1. DFT-optimized geometries and relative Gibbs energies (kcal · mol 1, computed with M06 functional) in 1,4-dioxane solvent of two conformations of the Pd(Josiphos)(NH3)(Ph)Br intermediate in the Buchwald Hartwig amination of arylhalides with one of the arms of the chelate phosphine Figure 2. Extended chemical model for ketone hydrogenation catalyzed by a open.[13] 3D-structures were generated using CYLview.[15] iridium-hydride complex, in alcohol solvent and with added base (NaOMe).[19] 3D-structures were generated using CYLview.[15] developed methods combining low-level, semiempirical, DFT has been very prolific in showing the impact of counteranions calculations with metadynamics simulations, enable fast QM in both kinetics and selectivity of transformations.[22] Counter- exploration of the conformational space of large systems.[14] ions were generally considered as innocent partners and not Using this minimal model and without adding further included in the chemical model, but an increasing number of molecules in the computational description, the speciation of computational papers have started taking into account counter- reactive species in the reaction vessel arises as an important anions, particularly when dealing with proton transfer steps. issue. The active form of the catalysts is often different from the Calculations have largely contributed to the understanding of species initially added, and in most cases it is not known. In the role of anions in gold-catalyzed reactions.[23] However, the such cases, several possibilities should be carefully scrutinized. inclusion of countercations in models is still a less common Attention should also be paid to the speciation of simple practice.[19] It must be highlighted that explicit treatment of reagents, which can be very solvent dependent. Reagents such counterions poses an added challenge for calculations: its as sodium tert-butoxide and copper diacetate keep aggregated solvation. structures in non-polar solvent environments, and their descrip- Additives are commonly used in homogeneous catalysis to tion as monomeric species in QM calculations could be improve the outcome of synthetically relevant reactions. misleading.[16,17] The same issue arises with iodosylbenzene However, in most cases their role is not well understood, and (PhIO) a polymeric solid of low solubility.[18] they are not commonly included in calculations. The synergistic To build up a reliable energy profile, all the species that can use of calculations and experiments has allowed a deeper influence the energy of the system in each step of a catalytic understanding of additive effects in a number of organometallic cycle should be taken into account. The correct mechanism reactions.[24] cannot be represented using an incomplete model. Thus, when The reaction to be simulated in a computer proceeds in a no convincing mechanism arises from the PES exploration with solvent. Therefore, the solvent must be a consubstantial part of the minimal model, my first recommendation would be to the chemical model and it needs to be taken into account in extend the chemical model. The role in the reaction mechanism order to achieve a good computational description of the of these additional molecules that expand the minimal model reaction. There are three main approaches for introducing the should be scrutinized. Figure 2 displays an example of an solvent in a QM calculation: i) implicit solvent model, in which extended chemical model for a ketone hydrogenation catalyzed the solvent molecules are replaced by a continuum polarizable by an iridium-hydride complex, in alcohol solvent and in the medium characterized by its dielectric constant ɛ; ii) hybrid presence of NaOMe additive.[19] cluster-continuum model, in which a limited number of solvent A first extension of the initial model can be done by molecules are introduced in the QM description of the system, including additional molecules of one reagent, usually in excess and the rest is represented as a continuum, and iii) explicit in the experimental conditions. For instance, it is now well solvent model, which includes a huge number of solvent established in hydrofunctionalization reactions that a second molecules, usually suitable to reproduce the density of the molecule of the nucleophile can help in the proton migration chemical system. The most used and simplest solvent descrip- step.[20] When ionic species are present there are always tion is the continuum model. In most cases, particularly when counterions in the reaction medium that might influence the dealing with low-polar, non-coordinating solvents, this model is reaction.[21] In some cases, there are experimental evidences of enough to compute accurate energies in solution. However, it is such effects, particularly regarding counteranions. Gold catalysis striking, to my mind, that the most common way of using Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2549 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 continuum models in recent literature is still performing single- very far from the actual experimental conditions.[29] Moreover, point calculations at gas-phase optimized geometries. Today, the key chemical magnitudes that govern the direction and rate most of the time there are no technical reasons that can justify of a reaction (equilibrium and rate constants) are not related to this approach. Only in the cases in which optimization in differences in internal energies, but in Gibbs energies. Energy solution of relevant species fails to properly converge despite profiles of reactions carried out in a solvent have to be built repeated attempts, mainly due to numerical artifacts related upon G values in solution. Transforming energies from E to G with grids, gas-phase optimization, followed by single-point implies adding additional approximations to those inherent to energy calculation in solution, is justified. Assuming the same the QM calculation and raises more choices to the modeler. geometries in the gas-phase and in the solvent can be specially Currently, exploration of the PES of real systems is habitually misleading when there are charged species in the medium. performed using a medium-size basis set (BS1), usually valence Ionic mechanisms, implying steps in which ions are formed double-ξ + polarization. It is a common and highly recom- from neutral molecules, cannot be properly accounted for using mended practice to refine the energy values with a very this model. State-of-the-art calculations of reaction mechanisms extended basis set (BS2) in order to remove basis-set super- demand the optimization and characterization of the stationary position errors (BSSE). This is an approximation based on the points in solution. assumption, usually true, that the optimized geometries with A continuum model is not able to describe specific both basis set are very similar. In this way the Gibbs energies in interactions between solute and solvent. Usually, they can be solution are obtained from equation 1: important in protic solvents, such as alcohol solvents and water. Thus, the use of a purely continuum model is not recom- Gsol ¼ Eelec ðBS2Þ þ ½Gsol ðBS1Þ-Eelec ðBS1Þ� þ DG1atm!1M (1) mended for such solvents. In the same way, if there is the possibility for solvent coordination or, even more, solvent The first term (Eelec) is the QM energy in solution (usually molecules can participate in the reaction, the inclusion of DFT) and includes the Gibbs energy of solvation.[29] Its accuracy explicit solvent molecules in the chemical model is mandatory. depends on the level of theory chosen (functional + basis set). Hybrid cluster-continuum models can be a good alternative to The second term embraces the enthalpic and entropic contribu- tackle these issues. A major concern in these models is the tions of the solute. It demands the partition functions for the number and relative position of explicit solvent molecules to be solute, obtained from a frequency calculation performed at the included. The influence of this number and their conformations optimized geometries (BS1). The last term accounts for the on reaction barriers should be carefully checked until con- change of standard state from gas phase (1 atm) to solution vergence is reached.[25] Automated methods to determine how (1 M) and at 298 K amounts 1.89 kcal · mol 1 for each species. It many explicit solvent molecules need to be added in hybrid is only relevant when the number of moles changes in the cluster-continuum schemes to capture most of the interaction reaction (Δn ¼ 6 0). between the solute and the environment have been recently The accuracy in computing Eelec is known as the electronic proposed.[26] structure issue. The quality to price ratio have made DFT The elucidation of reaction mechanisms when solvent methods almost the only option for computational studies of molecules play a direct role in the mechanism and are involved organometallic reactivity. But which functional should be used? in multiple ligand exchange processes has proven to be This is a question to ask ourselves before starting calculations, challenging. This occurs, for instance, in two textbook organo- which does not have a general answer. Strictly speaking, it metallic reactions: the Wacker process and the Grignard depends on the particular reaction to be studied. Moreover, reaction. Calculations using explicit solvent models, placing the catalytic reactions are multistep and the selected functional reagents in a box with a large number of solvent molecules, should be able to provide accurate results for all the steps. To water and tetrahydrofuran (THF), respectively, and performing further complicate matters, kinetic experimental studies that ab initio Molecular Dynamics (AIMD) simulations with them, might serve as benchmarks are not very frequent. Despite that, have solved long standing controversies about the mechanism some benchmarking studies have been carried out for one of both reactions.[27,28] metal and one particular reaction[30] or several elementary steps of metal-catalyzed transformations.[31] In recent years local coupled-cluster methods, such as DLPNO-CCSD(T) have 4. The Computational Model emerged as accurate and effective methods for benchmarking in big-sized systems.[32] Once the chemical model has been decided, the next step is My personal impression is that the choice of the functional the choice of the computational model. While the first election is not as important as it is generally thought to be. It is worthy has consequences for the possible mechanisms to be explored to perform some calibration of the functional, but in general (for instance, mechanisms with participation of dimers cannot energy profiles with any non-local GGA or meta-GGA functional be explored if only one molecule of the catalyst is included in will differ by just a few kcal · mol 1 when computing closed-shell the chemical system), the second one will have impact on the species. The most relevant aspect to take into account is accuracy of the numbers provided by the calculations. QM probably the percentage of Hartree-Fock (HF) exchange in the calculations are devised to afford the internal energy (electronic functional, especially in reactions involving open-shell species. energy, Eelec) of a molecular assembly in gas phase and at 0 K, It is not optional, however, whether to take into account Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2550 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 dispersion interactions or not. The use of dispersion-corrected pointed out the validity of the IGRRHO and the unnecessariness functionals or the incorporation of empirical dispersion of further correction, provided that dispersion forces are taken parameters[33] to non-dispersion-corrected functionals is man- into account in the calculation. The best currently available way datory. As already mentioned when describing continuum to handle entropic effects for species in solution is the IGRRHO models, optimizations have to be performed taking into approach, combined with continuum solvent treatment of account dispersion. With the amount of evidence accumulated solvation Gibbs energy, and taking into account the change of on the importance of dispersion effects when describing real reference standard state (ΔG1atm!1M term in equation 1).[7–9] organometallic systems,[34] and the low cost that its inclusion entails, it is hard to understand the reason why in many recent articles dispersion interactions are introduced only as single- 5. Success… and Limitations point energy calculations at the geometries optimized without dispersion. In the previous sections I have provided the ground for what I Unfortunately, the similarity of energy profiles obtained consider an accurate computing of an energy profile in the using different functionals of similar quality for closed-shell current times.[9] Nowadays, taking the necessary precautions systems does not hold true when there are open-shell species regarding both the chemical and computational models, DFT (radicals) involved in a reaction mechanism. Comparison of calculations can afford in most cases a reliable quantitative DFT-computed relative energies of different spin-states of one description of the energy landscape for the formation of the molecule shows a great disparity, depending particularly on the major product of a transition metal promoted reaction in HF exchange in the functional. An accurate calculation of the solution. In this way the determining states of a catalytic cycle energy splitting between low-spin and high-spin states is can be found, and by using the energetic span model turnover difficult due to its strong functional dependence.[35] For frequencies (TOF) can be calculated from the energy profile.[37] instance, taking two functionals widely used in reactivity Commonly, the energetic information provided by the studies, such as wB97XD or B3LYP D3, both incorporating energy profile is used to infer and propose a reaction dispersion effects, with a very large basis set, the high-spin mechanism. It should never be overlooked that a reaction [triplet, formally copper(II)] and low-spin [singlet, formally mechanism cannot be proven, but only disproved.[38,39] As copper(III)] isomers of an intermediate formed in the Cu(I)- commented before, the mechanism search often suffers from promoted aerobic C H bond oxidation (Figure 3) differ by 5.4 the “wish bias”, defined as the tendency on the investigator to (wB97XD) or 15.4 (B3LYP D3) kcal · mol 1.[36] When leaving the reach a desired result.[40] When dealing computationally with field of “comfortable” closed-shell PES, careful analysis of the very well-known reactions, as it happens in most cases, no big functional dependence is needed. surprises can be expected, the computed profile is just one As commented previously, electronic energies must be more on this kind of reaction and the mechanistic proposal is transformed into Gibbs energies by adding enthalpic and sound….until proven otherwise. However, there is not a unique entropic contributions. This is usually done by means of pathway for the most interesting reactions. As not all the frequency calculations in which the different contributions to possible routes have been explored, there could be a better the molecular entropy are estimated from the partition mechanism, involving lower barriers than the initially proposed. functions. In order to do that, the usual approach is the ideal Another shortcoming is that often only the most productive gas/rigid rotor/harmonic oscillator approach (IGRRHO). It was part of the reaction is computationally explored. Rarely does thought that this approach overestimated entropic effects in one reaction give a 100 % yield. In the reaction vessel, the solution, and several alternatives have been devised, most of catalyst can be deactivated and side reactions can happen, them using empirical corrections, such as dividing by 1/2 or 2/3 both facts reducing the reaction yield. Improving the reaction the translational entropy. A number of recent works have yield is a main goal in synthetic chemistry. Computational analysis of the off-cycle part of reaction mechanisms can help improving catalysts, preventing for instance their deactivation, as it has been shown for cross-coupling reactions.[41] Performing calculations once the experimental part of a project has been finished, bringing after-the-fact explanations, is still the most common use of calculations in transition metal catalysis that can be found in the literature. This approach has proved to be very useful to understand experiments. However, proceeding solely in this way misses out the full potential of contemporary simulation techniques. Simulations are mature enough to be incorporated at the experimental design stage and to guide the experimental exploration throughout a project’s development. In this way, and in synergistic collabo- Figure 3. Functional dependence of the Gibbs energy difference between the low-spin and high-spin isomers of a copper-hydroxo complex.[36] Relative ration with experiments, calculations can evolve from being a Gibbs energies in THF solvent in kcal · mol 1. 3D-structures were generated technique to “predict the past” to becoming a powerful tool to using CYLview.[15] forecast the future. In chemistry that means to design new Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2551 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 catalysts and to discover new reactions. In the next section I will outline how the combination of accurate QM calculations, described before, with emerging computational tools can make this possible. 6. Designing Catalysts and Discovering New Reactions Beyond calculation of energy profiles, computational organo- metallic chemistry is taking new directions to overcome short- comings and to reach the goals pointed out in the previous section. On one side, calculations must be closer to the experiments, regarding both the outcome and the chemical models. On the other side, they have to explore the full PES of Figure 4. Microkinetic modeling of the isomerization of 2-phenyl-2-meth- yloxetane 2 a to allylic alcohol catalyzed by BF3 · Et2O.[45] II: zwitterionic the system, overcoming the human bias. Finally, they can be intermediate; 3 a: homoallylic alcohol product; 4 a: E/Z-allylic alcohol combined with extensive data analysis and machine learning products. algorithms. All together, these advances approach the holy grail of catalyst design and can lead to the longed-for discovery of new reactions. The future is now: recent literature collects an explosion of articles devoted to these topics. For instance, in the last two years three review journals have published special issues devoted to machine learning in chemistry, focusing on the design of catalysts and materials.[42] Here, I will mention just a few examples. Computational mechanistic studies are aimed to locate intermediates and transition states in a catalytic cycle and to determine their relative energies. In this way a Gibbs energy profile is obtained, and a mechanistic proposal is inferred from it. However, the direct outcome of experiments is the evolution of concentrations of reagents, intermediates, and products through time. These values depend on the rate constants of the reaction steps, related with Gibbs energy differences, but also on the concentrations of the reagents, not accounted for in the PES exploration. Microkinetics modeling is a simple and useful method for the introduction of concentration effects in the description provided by the Gibbs energy profiles, getting simulations closer to experimental conditions.[43,44] The energy profile defines the steps of a complex reaction Figure 5. Snapshot from an explicit water AIMD simulation of the hydration and their energy barriers. In the framework of the Transition of alkenes catalyzed by gold(I) complexes. State Theory, the rate constant of each reaction step can be computed from its Gibbs energy of activation. This allows the definition of a system of differential equations, which can be atically investigating how the potential energy varies as a solved numerically by using the initial concentrations of all the function of few nuclear coordinates, as done in typical PES compounds as starting conditions.[43,44] The outcome of the exploration on reduced model systems, is not viable because of microkinetic modeling is the evolution of the concentration of the large number of particles involved. First-principles molec- each species with time (Figure 4),[45] which is precisely what ular dynamics (also called ab initio Molecular Dynamics, AIMD), most experiments measure. The computational cost of this in which atomic forces are computed from electronic structure procedure, which can also be used as additional testing of a calculations (namely, density functional theory), is the techni- proposed reaction mechanism, is very low. que of choice to investigate chemical events in solution when Improving the solvent description is another steppingstone the solvent has a key role in the reaction and/or where solvent in approaching computer models to “real chemistry” in solution. reorganization is important.[46] This can be done employing extended explicit representations These simulations make use of intensive computational of the solvent, built to reproduce experimental densities (Fig- resources, but the inclusion of explicit solvent allows the ure 5). Dealing with large model systems, however, poses a characterization of molecular events that would otherwise not challenge: how to explore their configurational space? System- be accessible in reduced model systems, as it has already been Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2552 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 proved in Wacker[27] and Grignard[28] reactions. For some DFT and ML methods can be combined to enable the processes, such as transfer hydrogenation in alcohol solvent, prediction of reactivity within large chemical spaces containing the use of explicit solvent simulations may provide mechanistic thousands of complexes.[59] Around 2500 QM transition state insights that differ from those obtained in vacuo or in structures in the chemical space of IrL3X Vaska’s complex were continuum models.[47] employed to train and test ML models that predict the H2 As commented above, the usual way of proceeding in activation barrier. Different ML algorithms were used for differ- computational studies of organometallic reactivity has typically ent purposes, including neural networks for accuracy, Gaussian been to explore a limited number of pathways, inferred from processes for transferability, and Gradient boosting for previous knowledge and/or chemical intuition of the person explainability.[60] A reactivity descriptor database based on QM performing the calculations, that lead to the major product of calculations of 130,000 organic molecules was used to predict the reaction. That means to explore a very small fraction of the regio-selectivity for three general types of substitution PES of the reacting system. The emergence, in the last years, of reactions.[61] Trends surrounding the thermodynamics of the the automated reaction path methods, promises to revolution- hydroformylation reaction catalyzed by group 9 metals bearing ize the way of exploring the PES. With these methods, which phosphine ligands have been analyzed using a data-driven have been successfully applied to organometallic reactions,[48–50] inspired workflow (data-powered volcano plots). The total data an unbiased exploration of the full PES can be performed. They set used consisted of 1510 catalytic cycles derived from DFT do not only allow exploration of multiple reaction paths computations and 491 catalytic cycles derived from machine accounting for product selectivity in connection with micro- learned profiles.[62] Electroanalytical techniques have been kinetic models but also pave the way for the discovery of combined with parameterization tools, which include DFT unexpected reactions in the full potential energy surface of the calculations, to uncover reaction mechanism in redox system. catalysis.[63] Catalyst optimization and, furthermore, design and discov- ery, employing the usual approach of thoroughly PES explora- tion with accurate quantum mechanical calculations, becomes 7. Conclusions impractical when a large chemical space (thousands of catalysts and reactions) has to be explored. In the framework of the In this essay I wanted to give a personal insight on the current explosion of the use of data-driven, machine learning (ML) state of computational methods applied to homogenous algorithms in all fields of knowledge, they have also burst onto catalysis processes and to point out some of the new directions catalyst design. Although these techniques are still at an early they are taking. The basic assumption is that modern QM-based stage in this field, they have already achieved impressive methods are able to give accurate energy profiles of an successes, particularly in the area of enantioselective organometallic reaction in solution entailing closed-shell spe- catalysis,[51,52] proving its potential. A number of very recent cies, provided that a proper choice of the chemical and reviews prove the fact that they are called to change the way in computational methods has been made. This proper choice, which catalysts are discovered.[51–55] which describes what I consider the state-of-the-art of these The starting point of data-driven tools is establishing a techniques, has been analyzed in the first part of the article. In relationship between a quantitative description of reactants, my opinion that means, basically, that all the species that can catalysts and reaction conditions with a property (for instance, significantly affect the energy landscape of the reaction should activity: Quantitative Structure-Activity relationship, QSAR). The be included in the QM description of the chemical model to be key ingredients of these models are the quantities (descriptors) computed. In addition, the conformational space of big and that are correlated with the properties of interest. Multivariate flexible ligands, as well as the speciation issues of catalyst and linear regression models correlate data to molecular descriptors. reagents should be carefully analyzed. Regarding the computa- For non-linear relationships in complex chemical data, machine tional model, optimization and characterization of stationary learning offers a wide variety of algorithms that can be used to structures should be done in solution with functionals taking build accurate, transferable, and explainable predictive models. into account dispersion effects. Refinement of the energies with However, many of these algorithms require big data sets. very extended basis sets is highly recommended, as well as Computational results from QM calculations can be used to some functional calibration for the particular reaction to be generate these required data needed to feed the ML models computationally studied. Entropic contributions can be com- (quantum-based ML).[56] The calculation of over 200,000 organic puted using the IGRRHO approach, with no further corrections. radical species[57] and over 80,000 organometallic compounds[58] In this way, accurate relative Gibbs energies of the species are impressive examples of how DFT methods can be used to along the computed pathway for the formation of product(s) generate big data sets. These data sets can be exploited in can be obtained. Success of the computational protocol that I data-driven predictive ML models for the discovery of new have described cannot be assured when more than one spin- catalysts and reactions. However, the inclusion of many state can be involved in the reaction. Such energy profiles are complexes which are not realistic/synthetically feasible to still challenging and require careful functional calibration. generate training data could lead to bias in a ML model’s The translation from energy profile to “the” mechanism of prediction. the reaction is not as direct as it is usually thought. In the best cases, it can be considered a sound proposal…until a better Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2553 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 one comes out. In most cases calculation of an energy profile is Keywords: Catalyst discovery · Chemical and computational performed when the experimental research on the system has models · Density functional calculations · Energy profile · been concluded and is used only to explain the observed Reaction mechanisms reactivity, missing out the actual potential of QM calculations. On the other side, an energy profile is the outcome of a biased exploration of a very small region of the PES of the system. [1] T. Sperger, I. A. Sanhueza, I. Kalvet, F. Schoenebeck, Chem. Rev. 2015, Automated reaction path methods can explore the full 115, 9532–9586. [2] S. Santoro, M. Kalek, G. Huang, F. Himo, Acc. Chem. Res. 2016, 49, 1006– reactivity of the system, accounting for the formation of side 1018. products and catalyst deactivation pathways and even discover- [3] K. D. Vogiatzis, M. V. Polynski, J. K. Kirkland, J. Townsend, A. Hashemi, C. ing new reactions than can take place within the reaction Liu, E. A. Pidko, Chem. Rev. 2019, 119, 2453–2523. [4] S. Ahn, M. Hong, M. Sundararajan, D. H. Ess, M.-H. Baik, Chem. Rev. 2019, mixture. 119, 6509–6560. To fully benefit from the current potential of QM calcu- [5] L. Perrin, K. J. T. Carr, D. McKay, C. L. McMullin, S. A. Macgregor, O. lations in homogenous catalysis, these must be used to test and Eisenstein, Struct. Bond. 2016, 167, 1–38. [6] Z. Liu, C. Patel, J. N. Harvey, R. B. Sunoj, Phys. Chem. Chem. Phys. 2017, to guide the experimental research throughout all the project’s 19, 30647–30657. development, from its initial design to the final optimization. [7] J. N. Harvey, F. Himo, F. Maseras, L. Perrin, ACS Catal. 2019, 9, 6803– However, the exploration of large chemical spaces required to 6813. [8] J. N. Harvey in: Transition Metals in Coordination Environments. Chal- optimize and discover new catalysts with the usual techniques lenges and Advances in Computational Chemistry and Physics, Vol. 29 of PES exploration is impractical. Increasingly, this is done using (Eds: E. Broclawik, T. Borowski, M. Radoń), Springer, Cham, 2019, pp data-driven techniques and machine learning algorithms, able 289–313. [9] O. Eisenstein, G. Ujaque, A. Lledós, Top. Organomet. Chem. 2020, 67, 1– to explore very large chemical spaces, but these approaches 38. often require big data sets. In this field, QM calculations have [10] I. Funes-Ardoiz, F. Schoenebeck, Chem 2020, 6, 1904–1913. found their place to generate big data sets. [11] H. Ryu, J. Park, H. K. Kim, J. Y. Park, S.-T. Kim, M.-H. Baik, Organometallics 2018, 37, 3228–3239. In the first lines of this essay, I asked myself if theoreticians [12] M. Besora, A. A. C. Braga, G. Ujaque, F. Maseras, A. Lledós, Theor. Chem. are still needed. Throughout the text, I have sought to show Acc. 2011, 128, 639–646. that my answer is yes, but only for as long as we are able to do [13] P. Gómez-Orellana, A. Lledós, G. Ujaque, J. Org. Chem. 2021, 86, 4007– 4017. things that other chemists cannot. Routine, ex post calculation [14] S. Grimme, J. Chem. Theory Comput. 2019, 15, 2847–2862. of energy profiles is nowadays within the reach of any chemist. [15] CYLview20; C. Y. Legault, Université de Sherbrooke, 2020 (http:// However, exploring large chemical spaces and finding new, www.cylview.org). [16] P. Villuendas, S. Ruiz, P. Vidossich, A. Lledós, E. P. Urriolabeitia, Chem. synthetically relevant reactions using state-of-the-art methods Eur. J. 2018, 24, 13124–13135. and exploiting them with machine learning methods, without [17] S.-T. Kim, S. Kim, M.-H. Baik, Chem. Sci. 2020, 11, 1017–1025. losing the physical meaning of the quantities handled, can only [18] G. Barea, F. Maseras, A. Lledós, New J. Chem. 2003, 27, 811–817. [19] J. M. Hayes, E. Deydier, G. Ujaque, A. Lledós, R. Malacea-Kabbara, E. be done by researchers combining a deep knowledge of Manoury, S. Vincendeau, R. Poli, ACS Catal. 2015, 5, 4368–4376. computational methods with a molecular perspective. With this [20] A. Couce-Rios, A. Lledós, I. Fernández, G. Ujaque, ACS Catal. 2019, 9, in mind, computational chemists can play an ever more central 848–858. [21] A. Macchioni, Chem. Rev. 2005, 105, 2039–2074. role in homogenous catalysis. There is still plenty of room for [22] M. Jia, M. Bandini, ACS Catal. 2015, 5, 1638–1652. computational organometallic catalysis. [23] G. Ciancaleoni, L. Belpassi, D. Zuccaccia, F. Tarantelli, P. Belanzoni, ACS Catal. 2015, 5, 803–814. [24] T. Sperge, I. A. Sanhueza, F. Schoenebeck, Acc. Chem. Res. 2016, 49, 1311–1319. Acknowledgements [25] M. A. Ortuño, A. Lledós, F. Maseras, G. Ujaque, ChemCatChem 2014, 6, 3132–3138. [26] G. N. Simm, P. L. Türtscher, M. Reiher, J. Comput. Chem. 2020, 41, 1144– This essay summarizes forty years of human learning (HL), fed 1155. with discussions and talks with all the people who have been in [27] A. Stirling, N. N. Nair, A. Lledós, G. Ujaque, Chem. Soc. Rev. 2014, 43, my group and experimental and theoretical colleagues around the 4940–4952. [28] R. M. Peltzer, J. Gauss, O. Eisenstein, M. Cascella, J. Am. Chem. Soc. 2020, world. I am indebted to all of them. Special thanks to Gregori 142, 2984–2994. Ujaque (UAB), Giuseppe Sciortino (ICIQ), and David Balcells (Oslo) [29] Calculations including a continuum solvent afford a “non-pure” for a critical reading of the manuscript. Financial support from thermodynamical magnitude, Esolv. This is an electronic energy that also includes the Gibbs energy effects of the solvent, but not those of the Spanish Ministry of Science and Innovation (project PID-2020- solute. 116861GB-I00) is acknowledged. [30] K. Rohmann, M. Hölscher, W. Leitner, J. Am. Chem. Soc. 2016, 138, 433– 443. [31] K. H. Hopmann, Organometallics 2016, 35,3795–3807. [32] D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. L. Martin, F. Neese, J. Conflict of Interest Chem. Theory Comput. 2015, 11, 1525–1539. [33] S. Grimme, WIREs Comput. Mol. Sci. 2011, 1, 211–228. [34] See, for instance: E. Lyngvi, I. A. Sanhueza, F. Schoenebeck, Organo- The authors declare no conflict of interest. metallics 2015, 34, 805–812. [35] M. Swart, Top. Organomet. Chem. 2020, 67, 191–226. [36] M. Álvarez, F. Molina, M. R. Fructos, J. Urbano, E. Álvarez, M. Sodupe, A. Lledós, P. J. Pérez, Dalton Trans. 2020, 49, 14647–14655. [37] S. Kozuch, S. Shaik, Acc. Chem. Res. 2011, 44, 101–110. [38] O. Eisenstein, Faraday Discuss. 2019, 220, 489–495. Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2554 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH Essays doi.org/10.1002/ejic.202100330 [39] S. L. Scott, ACS Catal. 2019, 9, 4706–4708. [53] M. Foscato, V. R. Jensen, ACS Catal. 2020, 10, 2354–2377. [40] E. L. Wynder, I. T. Higgins, R. E. Harris, J. Clin. Epidemiol. 1990, 43, 619– [54] L. C. Gallegos, G. Luchini, P. C. St. John, S. Kim, R. S. Paton, Acc. Chem. 621. Res. 2021, 54, 827–836. [41] D. Balcells, A. Nova, ACS Catal. 2018, 8, 3499–3515. [55] D. J. Durand, N. Fey, Acc. Chem. Res. 2021, 54, 837–848. [42] a) D. Ess, L. Gagliardi, S. Hammes-Schiffer, Chem. Rev. 2019, 119, 6507– [56] O. A. von Lilienfeld, K.-R. Müller, A. Tkatchenko, Nat. Chem. Rev. 2020, 4, 6508 (Special Issue: Computational Design of Catalysts from Molecules 347–358. to Materials); b) H. J. Kulik, M. S. Sigman, Acc. Chem. Res. 2021, 54, [57] P. C. St. John, Y. G. Y. Kim, B. D. Etz, S. Kim, R. S. Paton, Sci. Data 2020, 2335 2336 (Special issue: Data Science Meets Chemistry); c) R. Gómez- 7 : 244. Bombarelli, A. Wiltschko, Trends Chem. 2021, 3, 75–76 (Special Issue: [58] D. Balcells, B. B. Skjelstad, J. Chem. Inf. Model. 2020, 60, 6135–6146. Machine Learning for Molecules and Materials). [59] G. d. P. Gomes, R. Pollice, A. Aspuru-Guzik, Trends Chem. 2021, 3, 96– [43] M. Besora, F. Maseras, WIREs Comput. Mol. Sci. 2018, 8, e1372. 110. [44] M. Jaraíz, Top. Organomet. Chem. 2020, 67, 81–106. [60] P. Friederich, G. d. P. Gomes, R. De Bin, A. Aspuru-Guzik, D. Balcells, [45] A. Cabré, S. Rafael, G. Sciortino, G. Ujaque, X. Verdaguer, A. Lledós, A. Chem. Sci. 2020, 11, 4584–4601. Riera, Angew. Chem. Int. Ed. 2020, 59, 7521–7527. [61] Y. Guan, C. W. Coley, H. Wu, D. Ranasinghe, E. Heid, T. J. Struble, L. [46] P. Vidossich, A. Lledós, G. Ujaque, Acc. Chem. Res. 2016, 49, 1271–1278. Pattanaik, W. H. Green, K. F. Jensen, Chem. Sci. 2021, 12, 2198–2208. [47] J. W. Handgraaf, E. J. Meijer, J. Am. Chem. Soc. 2007, 129, 3099–3103. [62] M. D. Wodrich, A. Fabrizio, B. Meyerac, C. Corminboeuf, Chem. Sci. 2020, [48] J. A. Varela, S. A. Vázquez, E. Martínez-Núñez, Chem. Sci. 2017, 8, 3843– 11, 12070–12080. 3851. [63] T. Tang, C. Sandford, S. D. Minteer, M. S. Sigman, Chem. Sci. 2021, 12, [49] A. L. Dewyer, A. J. Argüelles, M. Zimmerman, WIREs Comput. Mol. Sci. 4771–4778. 2018, 8, e1354. [50] M. Hatanaka, Y. Takayoshi, S. Maeda, Top. Organomet. Chem. 2020, 67, 57–80. [51] M. S. Sigman, K. C. Harper, E. N. Bess, A. Milo, Acc. Chem. Res. 2016, 49, 1292–1301. [52] A. F. Zahrt, S. V. Athavale, S. E. Denmark, Chem. Rev. 2020, 120, 1620– Manuscript received: April 23, 2021 1689. Revised manuscript received: May 24, 2021 Eur. J. Inorg. Chem. 2021, 2547 – 2555 www.eurjic.org 2555 © 2021 The Authors. European Journal of Inorganic Chemistry published by Wiley-VCH GmbH