Topics in Catalysis (2022) 65:118–140 https://doi.org/10.1007/s11244-021-01520-2 ORIGINAL PAPER Beyond Continuum Solvent Models in Computational Homogeneous Catalysis Gantulga Norjmaa1 · Gregori Ujaque1 · Agustí Lledós1 Accepted: 7 October 2021 / Published online: 16 November 2021 © The Author(s) 2021 Abstract In homogeneous catalysis solvent is an inherent part of the catalytic system. As such, it must be considered in the computa- tional modeling. The most common approach to include solvent effects in quantum mechanical calculations is by means of continuum solvent models. When they are properly used, average solvent effects are efficiently captured, mainly those related with solvent polarity. However, neglecting atomistic description of solvent molecules has its limitations, and continuum solvent models all alone cannot be applied to whatever situation. In many cases, inclusion of explicit solvent molecules in the quantum mechanical description of the system is mandatory. The purpose of this article is to highlight through selected examples what are the reasons that urge to go beyond the continuum models to the employment of micro-solvated (cluster- continuum) of fully explicit solvent models, in this way setting the limits of continuum solvent models in computational homogeneous catalysis. These examples showcase that inclusion of solvent molecules in the calculation not only can improve the description of already known mechanisms but can yield new mechanistic views of a reaction. With the aim of systema- tizing the use of explicit solvent models, after discussing the success and limitations of continuum solvent models, issues related with solvent coordination and solvent dynamics, solvent effects in reactions involving small, charged species, as well as reactions in protic solvents and the role of solvent as reagent itself are successively considered. Keywords Continuum solvent models · DFT calculations · Explicit solvent · Micro-solvation · Reaction mechanisms · Solvent effects · Solvent models 1 Introduction intermediates and transition states. The three main solvent properties accounting for its interactions with the reaction In homogenous catalysis the solvent is an inherent part of the system are polarity (related with the dielectric constant ε) catalytic system. Although for years it had been considered and Brönsted and Lewis acidity and basicity. These interac- just as a passive medium where reaction events take place, tions can have a deep impact on reaction rate and selectivity. its critical role on catalyzed reactions is nowadays well rec- With this in mind it is clear that solvent must be considered ognized [1]. Indeed, selection of the best solvent to carry when discussing and computing reaction mechanisms. out one reaction is crucial for the reaction performance and Three main approaches have been developed to deal with is a mandatory step in the optimization of reaction condi- solvent effects in computational studies of reaction mecha- tions [2]. At the macroscopic level, the solvent influences the nisms using quantum chemistry schemes: implicit (or con- solubility of the species present along the reaction pathway. tinuum), cluster/continuum (or hybrid implicit/explicit) and At the molecular level, solvent molecules interact directly fully explicit solvation models (Fig. 1). There is a vast lit- with substrates, products and catalyst molecules, as well as erature on this subject and here we will only briefly recall the main characteristics of each approach. Recent articles * Gregori Ujaque give a comparative overview of the three approaches [3, 4]. gregori.ujaque@uab.cat In continuum solvent models (CSM) no explicit solvent * Agustí Lledós molecules are introduced in the chemical system to be com- agusti.lledos@uab.es puted. In this implicit approach the solvent is a homogene- ous polarizable medium, described by its dielectric constant 1 Departament de Química, Universitat Autònoma de ε. The solute is placed in a cavity that is created in this Barcelona, 08193 Cerdanyola del Vallès, Catalonia, Spain 13 Vol:.(1234567890) Topics in Catalysis (2022) 65:118–140 119 Fig. 1  Schematic representation of solvent models medium. Electron density of the solute inside the cavity continuum solvation to account for the long-range elec- polarizes the dielectric, continuum, which in turns acts on trostatic effects, thus becoming a hybrid implicit–explicit the solute, modifying its electron distribution. In this way, solvation scheme. A critical discussion of several hybrid the electrostatic contribution to the solvation Gibbs energy discrete-continuum models has been recently published [9]. is evaluated. Evaluation of non-electrostatic contributions, Appealing as it is, this cost-effective approach has a main to be added to the electrostatic one, is required to obtain shortcoming: the difficulty in defining the number and loca- accurate solvation energies [5]. Differences in the theoretical tion of the explicit solvent molecules. Moreover, the con- foundations, the way of defining the cavity and the boundary formational space strongly increases when adding several conditions give rise to various implicit solvation models, solvent molecules, making it necessary extensive confor- as the original Polarizable Continuum Model (PCM), the mational sampling [10]. Recently automated methods to Conductor-like Screening Model (COSMO) and the Solva- determine how many explicit solvent molecules need to be tion Model based on Density (SMD) among others. A thor- added in hybrid cluster-continuum schemes to capture most ough review on the current state-of-the art of CSM have of the interaction between the solute and the environment been recently published [6]. These methods are computa- have been proposed [11, 12]. tionally efficient, keep the system size as in the gas phase A realistic description of the condensed phase environ- and avoid sampling over solvent degrees of freedom. With ment in which catalytic reactions occur is achieved with these advantages, continuum solvent models have become fully explicit solvent models: the solute is placed in a box the most employed approach in computational homogeneous solvated by as many solvent molecules as needed to repro- catalysis. Very efficiently they combine an accurate, quan- duce the density of the system. These computational models tum mechanical (QM) description of molecular systems with can simulate the dynamics of the solvent, but to sample the simple but physically rooted description of environment [7]. vast configurational space of solvent configurations statisti- However, they provide “averaged” solvent effects and do not cal methods are required, for instance within the ab initio capture specific interactions with solvent molecules, which molecular dynamics framework (AIMD) [13]. These cal- in some cases may be important. It is not an appropriate culations are generally tremendously demanding in com- model if solvent molecules can directly participate in the putational resources. Moreover, acceleration techniques as reaction. metadynamics or umbrella sampling must be employed to These limitations can be overcome, in principle, by access time scales at which reactive events (“rare events”) including a selected number of explicit solvent molecules happen. It should be pointed out that despite the tremendous into the atomistic part of the calculation. This is the main progress achieved in modeling chemical events in solution, feature of the cluster-continuum approach, that is also modeling of mixed solvents is still very challenging [14, 15]. referred as “micro-solvation” [8]. The explicit solvent mol- As commented above, nowadays practically all the com- ecules are treated in the same theoretical framework as the putational studies of reaction mechanisms in homogenous rest of the catalytic system. This “selected” explicit solva- catalysis performed with Density Functional Theory (DFT) tion is usually supplemented with the addition of implicit calculations consider the solvent, mostly by means of 13 120 Topics in Catalysis (2022) 65:118–140 continuum models. This methodology has been very use- reason to skip optimization in solution, regrettably a still ful to approach the computed model of the system to the common practice. Not doing that can lead to misconcep- chemical “real world” [16]. Thus, according to the success tions in reaction mechanisms. of plain continuum models, sometimes it is believed that Dehydrogenation of a gem-diol to form carboxylic acid is they can be applied to any situation. The main goal of this a proposed step in the ruthenium-catalyzed transformation chapter is to showcase, presenting selected examples, that of primary alcohols to carboxylic acid salts and ­H2 using there are situations in which proper modeling of the cata- water as the -OH group source. The reaction was carried out lytic system demands going beyond of implicit solvation; in in basic aqueous solution [25]. these cases, the inclusion of explicit solvent molecules in the From calculations in the gas-phase the dehydrogena- quantum mechanical framework is mandatory for a proper tion of the gem-diol proceeds over a concerted transition description of the system. The article is not aimed to collect state where the proton and hydride transfer simultaneously an exhaustive list of calculations including explicit solvent, (Fig. 2). Therefore, this would be the conclusion about the but to highlight with examples the reasons behind the choice reaction mechanism if the solvent is introduced by means of of micro-solvated or fully explicit solvent models, and the single-point calculations in a continuum medium at the gas- benefits provided by this choice. In this way we would like phase optimized geometries. However, optimizing in con- to contribute to set the limits of continuum solvent models tinuum water it is found that the same reaction proceeds via in reactivity studies. With the aid of examples, we will start two transition states, the first one describing hydride transfer showing the success and limitations of CSM. Then, we will and the second one proton transfer. An ion-pair intermedi- address first issues related with solvent coordination and ate connects both transition states [26]. When eight explicit solvent dynamics. This section will be followed by an over- water molecules are further added to the computations, the view of solvent effects in reactions involving small, charged first step (hydride transfer) proceeds similarly, but with a species. Next section will be devoted to reactions in protic barrier lowered about 5 kcal ­mol−1. Importantly, in the sec- solvents and finally we will focus on the role of solvent as ond step (proton transfer) the source of the proton is not the a reagent itself. C-H functionality, but the solvent [26]. Inclusion of solvent effects not only changes the reaction mechanism from con- certed to stepwise but changes the overall view of its mecha- 2 Success and Limitations of Continuum nism. This example demonstrates the importance of properly Solvent Models accounting for solvent effects in mechanistic studies. When correctly applied to suitable reactions, continuum The success of continuum models in modeling homogene- solvent models provide an accurate description of reactions ous catalysis processes is undeniable and has been already in solution. A recent study of the mechanism of peroxyfor- demonstrated on numerous cases [17–20]. Nowadays its mate decomposition (Scheme 1a) illustrates the potential of use has become mandatory when computing energy pro- CSM, when used correctly, to differentiate solvent effects in files, even for reactions in solvents with low dielectric a reaction strongly affected by solvent polarity [27]. constant. The current models have been carefully param- Experimental kinetic constants for this reaction are avail- eterized to reproduce known experimental solvation able for twenty solvents and vary by a factor of 140 [28], Gibbs energies. When they are correctly used [21], these making it an excellent case for benchmarking. The transition methods, coupled with accurate methods for computing state (Scheme 1b) shows an important degree of charge sepa- the potential energy (DFT methods including dispersion ration and its relative stability with respect to the reactant effects or highly-correlated post-Hartree–Fock approaches is expected to be deeply affected by the polarity of the sol- as CCSD(T) or DLPNO/CCSD(T)), deliver accurate elec- vent. DFT calculations, performed with the SMD continuum tronic energies in solution E solv, that include the Gibbs model in a set of 29 solvents (from heptane, ε = 1.91, to energies of solvation. Adding to this thermodynamically formamide, ε = 108.94) give Gibbs energy barriers in very “non-pure” magnitude, Esolv, thermal and entropic correc- good agreement with the experimental values (RMSD value tions for the solute, solution-phase Gibbs energies, suita- of 2.4 kcal ­mol−1 for the DLPNO/CCSD(T) energies with ble to build energy profiles, are obtained. It is important to implicit solvent). Moreover, optimization in solvent reveals stress that, as we [16, 22] and others [23, 24] have pointed subtle changes in the TS structure due to the solvent polar- out, proper use of CSM for modeling reaction mechanisms ity, that can be correlated with the barriers [27]. This work in solution implies optimization and characterization of all clearly shows that when the solvent influence is due to its the structures in solution. Only in the cases in which opti- polarity continuum models can capture solvent effects. mization in solution fails due to technical issues, gas-phase Computational studies aimed to compare implicit and optimization followed by single-point energy calculation explicit solvent models in reaction mechanisms are very in solution is justified. Saving computer time cannot be a relevant to delineate the limits of continuum models and to 13 Topics in Catalysis (2022) 65:118–140 121 Fig. 2  Ru-catalayzed dehy- drogenation of gem-diol into butyric acid computed with different solvent models [26] 13 122 Topics in Catalysis (2022) 65:118–140 Scheme 1  a Peroxyformate decomposition reaction. b Proposed transition state determine the origin of their shortfalls. Calculations with by implicit solvent models with a reduced computational explicit and implicit solvent models have been compared on effort [29]. non-catalyzed and Ag-catalyzed ring-closing reaction steps The Gibbs energy profile for the reductive elimination in a polar, aprotic solvent (N,N-dimethyl formamide, DMF, step of Suzuki–Miyaura couplings have been computed ε = 37.22) [29]. The reaction routes compared for the furan- using explicit and implicit solvation approaches in a diverse ring formation have different charge states (Scheme 2). set of solvents (benzene, toluene, DMF, ethanol and water) The main conclusion of this study is that for this type of [30]. Although a good agreement between both approaches reactions the effect of the DMF solvent can be adequately was obtained for aprotic solvents, the correlation was poor described by continuum solvent models: barriers computed for protic solvents. The difference was attributed to the for- without and with explicit solvent molecules are similar. mation of Pd···HO hydrogen bonds between Pd(PPh3)2 and Analysis of trajectories obtained with molecular dynamics protic solvent molecules. This work is a warning about the simulations with explicit solvent indicate that despite pair- use of implicit solvation for modeling reactions in protic wise interactions between the solute and the polar solvent solvents. The reductive elimination from a Au(III) complex molecules are significant, they do not induce any kind of to form a C–C bond in MeOH was shown to be also highly ordering in the solvation environment. Solvent molecules affected by micro-solvation [31]. Similar conclusions, for are collectively fluctuating around the solute molecule. This aqueous phase reactions, were inferred from the comparison mobile, fluctuating solvent shell can be efficiently substituted Scheme 2  Ring-closing reac- tions computed in DMF with implicit and explicit solvent models [29] 13 Topics in Catalysis (2022) 65:118–140 123 of different solvation models for the ­CO2 reduction by to zero are defined as “ambiguous coordinator”, with similar ­NaBH4 in aqueous solution [32]. possibilities to appear coordinated and uncoordinated. For instance, this is the case for acetonitrile and tetrahydrofuran (THF) [34]. Moreover, solvent coordination arises from an 3 Solvent Coordination and Solvent equilibrium between enthalpic and entropic terms, with the Dynamics latter being quite difficult to evaluate. We analyzed the speciation of Negishi coupling reagents, Solvent dynamics can affect the course of a reaction despite as ­ZnMe2, ZnMeCl, and ­ZnCl2 in THF solution, by experi- the solvent molecules may not be involved in the general mental (infrared (IR) and calorimetric experiments) and stoichiometric or catalytic process. One of the simplest DFT calculations [35]. IR observation suggests the metal ways to exhibit this behavior is by solvent coordination to a center undergoes the THF coordination equilibria depicted metal catalyst. Solvent molecules do move and, if they have in Fig. 3. Calculations sustain this proposal: the coordination basic character, they can coordinate to Lewis acid centers, of a second THF is clearly preferred over the coordination of as unsaturated transition metals [33]. Such a coordination, a single THF molecule to ­ZnMe2 (Scheme 3). that is a particular case of association-dissociation process The solvent was modeled using a cluster/continuum in a condensed phase, affects the speciation of reactive spe- approach employing (THF)n clusters. In this way for n ≥ 3, cies and can impact the reaction mechanism. From the point the total number of components (two) remains constant at of view of modeling this possibility raises the question of both sides of the equilibria, as shown in Scheme 3, eliminat- whether coordinated solvent molecules should be included ing a major contribution to entropy changes upon coordina- in the QM description of the system. For the case of good tion. The study was also performed to ZnMeCl and Z ­ nCl2; coordinating solvents (e.g., pyridine, dimethyl sulfoxide in all cases (THF)4 appears a good enough representation (DMSO), formamide, etc.), generally the answer is yes, but of the solvent. The conclusion of this study was that any the answer is not as easy for weaker coordinating solvents. computational modeling in THF should consider the spe- A scale that attempts to quantify the weakly coordinating cies ­ZnMe2(THF)2, ZnMeCl(THF)2, and ­ZnCl2(THF)2 as the character of a variety of solvents has been proposed and a ones participating in reactivity [35]. coordinating ability index has been calculated from struc- Solvent coordination can happen even with solvents with tural data base searches [34]. This index gives a first indica- weak coordinative properties, as toluene. We investigated tion of the feasibility of solvent coordination (positive values Pd(PPh3) speciation in toluene (Tol) by means of explicit means coordinating solvent and negative non-coordinating solvent simulations, using AIMD techniques [37]. Calcu- ones). Those solvents with a coordinating ability index close lations provided computational evidence for the formation Fig. 3  a THF coordination equilibria suggested by IR measurements. b Optimized structures [35]. 3D-structures were generated using CYLview [36] 13 124 Topics in Catalysis (2022) 65:118–140 Scheme 3  DFT calculated ΔG values (kcal ­mol−1) for reactions 1 and 2 in THF solvent. n is the number of solvent molecules in the solvent cluster of [Pd(PPh3)]-solvent adducts, as (Tol)Pd(PPh3)2, (Tol) in explicit solvent made possible solving this longstanding Pd(PPh3) and (Tol)2Pd(PPh3). An important conclusion from chemical puzzle. this study was that the bare Pd(PPh3) often proposed in cata- The transformation of Grignard reagent C ­ H3MgCl into lytic cycles, is actually not accessible in solution. Therefore, ­MgCl2 and Mg(CH3)2 in tetrahydrofuran solvent (Fig. 4a) solvent coordination should be considered when modeling was characterized at the molecular level by means of AIMD the reactivity of Pd(PR3) complexes in solution. simulations with enhanced-sampling metadynamics, using Solvent coordination issues are even more sensitive when the simulation box displayed in Fig. 4c containing 42 THF the reaction involves main group alkaline and alkaline-earth molecules as a solvent model [39]. The interchange of the metal cations, without clear coordinative preferences and methyl and chloride groups between two magnesium centers with a mainly electrostatic interaction with solvent mole- is associated with a change in the solvation number. Asym- cules. Extensive investigation of the most stable aggregates metric solvation on Mg ions is needed to promote the Mg-Cl formed in these systems is mandatory. For instance, in a and Mg-CH3 bond cleavage. Thus, solvent dynamics is driv- recent study of the Pummerer coupling between sulfoxides ing the ligand exchange in the Schlenk equilibrium. and turbo-organomagnesium amides (Scheme 4), several Using similar methodology than for the Schlenk equi- structural possibilities, including explicit THF molecules, librium, two competing mechanisms, the polar (or nucleo- have been checked. The investigation was able to deter- philic) implying heterolytic Mg-C bond breaking, and the mine the most stable species that Grignard reagent (isopro- radical, involving homolytic Mg-C bond breaking have been pylmagnesium chloride) and the magnesium amide (“R- computed for the Grignard reaction (Fig. 4b) [40]. Both MgX” + DIPAMgCl·LiCL) can form prior to the addition mechanisms are strongly dependent on the coordination of of the sulfoxide [38]. the Mg center. The solvent plays an essential role in the reac- Grignard reagents, with nominal formula RMgX, are tion: the more solvated the Mg species are, the more reac- seemingly simple organometallic species, however unrave- tive become. Thus, as for the Schlenk equilibrium, solvent ling its speciation in THF solution and obtaining a detailed dynamics is essential for representing the energy profile. understanding of the mechanism by which they operate in Inclusion of explicit solvent can impact the mechanistic the Grignard reaction have pushed solvent models to the view of a reaction, even not governing the reaction course. limit [39, 40]. This is because this system constitutes a pro- This has been the case for the allylic oxidation of cyclohex- totypical example of a reaction governed by solvent dynam- ene in acetonitrile (ε = 35.69) using a Cu(II) complex. Static ics. AIMD simulations and enhanced-sampling methods DFT calculations of the reaction mechanism with a con- tinuum description of acetonitrile resulted in an estimation Scheme 4  a Pummerer coupling between sulfoxides and Grignard reagent (isopropylmagnesium chloride) mediated by DIPAMgCl·LiCl. b New species are formed in solution upon mixing Grignard reagents and DIPAMgCl·LiCl [32] 13 Topics in Catalysis (2022) 65:118–140 125 Fig. 4  a Schlenk equilibrium of Grignard reagents. b Grignard reaction. c Simulation box used in the AIMD study (reproduced from [33]) of high barriers: 43.7 kcal ­mol−1 for the H-abstraction step water molecules yields a ­pKa of 4.4, closer to the experi- and 37.5 for the step in which cyclohexanone product is mental value. For this reason, seven explicit water molecules released, respectively [41]. Umbrella-sampling DFT-based were included in a computational study of the reduction of molecular dynamics (DFT-MD) simulations in a periodic ­CO2 to methanol by dihydropteridine [43]. Similar conclu- box containing 60 acetonitrile molecules opened up a new sions have been obtained for the calculation of redox poten- perspective on the mechanism for this reaction, identifying tials [44, 45]. new reaction branches and side processes. Explicit solvent In a computational study of the transmetalation step in environment of the simulation stabilizes reactive species that the Suzuki–Miyaura cross-coupling reaction we analyzed open two new reaction pathways involving dehydrogenation the bromide-by-hydroxide replacement leading to a pal- and re-hydrogenation of the -NH2 group bound to the Cu- ladium hydroxo complex [46]. Without including explicit site, with significantly lowered barriers relative to the prior solvent molecules, the transition state for this substitution is static DFT studies with implicit solvent [42]. located at 32.5 kcal ­mol−1 above the energy of the reactants (Fig. 5a). Such a high energy barrier does not agree with experimental results [47]. This poor result can be explained 4 Solvation of Small, Charged Reagents by the inaccurate representation of the base as a naked ­OH−. The reaction is performed in THF solvent, but the presence It is very common in homogenous catalysis processes to find of water in the experiments is highly likely: they are coming steps in which participate small and charged reagents, as for from hydration molecules of the base itself, or for instance, instance ­OH−, ­HCOO−, ­Cl− or ­H3O+. Charged species are from the solvent when using mixtures as THF/water [47]. also often engaged in ligand exchanges. As the following The inclusion of water molecules solvating the two exchang- examples highlight, inclusion of explicit solvent molecules ing anions significantly decreases the barrier. The influence is required for a realistic modeling of such processes. in the barrier of the number of water molecules considered A simple example illustrates the failure of continuum sol- (n) was evaluated. The series was considered converged at vent models to describe the stabilization of small anionic n = 3 because the change in the Gibbs energy profile does species by a polar solvent: calculation of p­ Ka of HCOOH in not significantly varies with respect to n = 4 (24.2 and water with such a model gives a value of 11.9, very far from 22.9 kcal ­mol−1 for n = 3 and 4, respectively) (Fig. 5b). the experimental value of 3.9. This failure is mainly due to Protodemetallation is another process highly influenced the absence of explicit H-bonding by solvating H ­ 2O mol- by the presence of solvent. We remark here protodeauration, ecules that stabilize the anionic conjugate base ­(HCOO−). A the formal exchange of Au(I) by a proton as a prototypi- cluster-continuum model with two explicit water molecules cal case. It is a very common final step in gold(I)-catalyzed H-bonding the anion only slightly improves the ­pKa value reactions of unsaturated C–C bonds. The barrier is very (9.8). However, enlarging anion micro-solvation to seven sensitive to the solvation of the electrophile. For example, 13 126 Topics in Catalysis (2022) 65:118–140 Fig. 5  Transition state for the bromide-by-hydroxide exchange in the Pd(PPh3)2 coordination sphere. a Without explicit water molecules. b With four explicit water molecules. Main distances in Å and Gibbs energy barriers in kcal ­mol−1. 3D-structures were generated using CYLview [36] the protodeauration of an Au(PMe3)(alkenyl) compound by pairs can be high energy species that could be wrongly a proton solvated by one, two, three, and four water mol- discarded in reaction pathways. For instance, complex ecules was investigated by DFT calculations (structure [Ir(pic)2(MeOC8H12)] (pic = picolinate ligand) undergoes 4 in Fig. 6). The authors found that the activation energy a fast isomerization between 2a and 2b structures (Fig. 8) increases when the number of the water molecules solvating [50]. A plausible pathway for this isomerization would entail the proton increases (Fig. 6) [48]. Solvation is modifying the cleavage of the C–OCH3 bond to form the ionic pair the electrophilicity of the protonating agent, overestimated [Ir(pic)2(cod)][OCH3]. Using a continuum solvent model using a purely continuum solvation model. When the num- such ion pair is found 28.4 kcal ­mol−1 above 2a, a value ber of water molecules solvating the proton increases, the that suggests discarding this mechanism. However, incor- electrophilicity of the proton decreases, as it is reflected by porating a small cluster of five MeOH molecules solvat- an increase in the LUMO energy of the protonating agent. It ing the methoxide anion, the ion pair was found to be only follows that the less electrophilic the protonating agent, the 11.0 kcal ­mol−1 above 2a. Using this cluster-continuum higher is the activation energy for protodeauration. model in DFT calculations the process in Fig. 8 takes place It must be stressed that explicit solvation not only modi- through a fairly accessible transition state of 20.9 kcal ­mol−1 fies acid–base properties of species as O ­ H− or ­H3O+, but (TSa, Fig. 8). also affects weak basic anions, as a chloride. In a combined The mechanistic investigation of the catalytic transfer experimental/computational investigation we disclosed hydrogenation of ketones to alcohols in water using a water- counteranion assistance in ruthenium-mediated alkyne to soluble ruthenium catalyst and aqueous HCOONa/HCOOH vinylidene isomerizations in methanol [49]. The theoretical as the hydrogen source at pH 4.4 illustrates the usefulness of study outlined the feasibility of hydrido − alkynyl species careful employing of cluster-continuum models to deal with as intermediates, in which the hydride has an acidic char- reaction mechanisms. acter. Due to the acidic character of this intermediate, basic In the experimental investigation, aqua, formato, and species as a chloride anion can act as a proton shuttle, tak- hydrido species were detected by 1H NMR experiments in ing the hydrogen from the metal to place it at the β-carbon, ­D2O [51]. Overall, the experimental evidence is consist- forming the vinylidene product. The basic character shown ent with the following sequence of reactions: (i) chloride by the chloride anion decreases when one solvent methanol by water ligand exchange, (ii) water by formate ligand molecule is located in its direct environment. Accordingly, exchange, (iii) decarboxylation of the formato complex the Gibbs energy barrier for migration increases from 19.5 to produce a hydrido species, and (iv) hydride and pro- to 23.2 kcal ­mol−1. Transition states in both models differ ton transfer to yield the reduced product (Fig. 9a). This considerably (Fig. 7). sequence involves several of the situations we have been The description of ion pairs can be also highly affected describing so far: a solvent molecule becomes a ligand, by the presence of explicit solvent molecules in the model. small and charged species are exchanged, and a charged Without inclusion of explicit solvent molecules, ion reagent (proton from the water medium) is added. DFT 13 Topics in Catalysis (2022) 65:118–140 127 Fig. 6  Energy profile calculated for protodeauration of the alkenyl complex 4 with an increasing number of water molecules solvating the ­H3O+ protonating agent [48] Fig. 7  Transition states for the Cl-assisted H-migration in Ru- mediated alkyne to vinylidene isomerizations in methanol, without and with a solvent molecule solvating the chloride anion. Main distances in Å and Gibbs energy barriers in kcal ­mol−1. 3D-structures were gen- erated using CYLview [36] 13 128 Topics in Catalysis (2022) 65:118–140 Fig. 8  a 2a to 2b isomerization. b Relevant structures in the isomerization process com- puted with a cluster-continuum solvent model including five MeOH molecules solvating methoxide anion. In blue, relative Gibbs energies, in kcal ­mol−1. 3D-structures were gen- erated using CYLview [36] Fig. 9  a Sequence of reactions for the transfer hydrogenation of cyclohexanone with a ruthe- nium water-soluble complex. b Transition state for the hydride transfer step with three and five water solvent molecules. 3D-structures were generated using CYLview [36] 13 Topics in Catalysis (2022) 65:118–140 129 Scheme 5  Ru-catalyzed hydrogen transfer via the metal–ligand bifunctional mechanism. Original proposed concerted mechanism calculations, using a cluster-continuum solvent model, pro- transfer hydrogenation in methanol solvent [52]. To model vided a detailed microscopic description of the full catalytic the solvated system a Ru complex and 40 methanol mol- cycle. Three explicit water molecules were included in the ecules were included in the simulation. From gas-phase QM description of the system in addition to the polariz- calculations, a concerted transfer of a proton and a hydride able continuum model (water, ε = 78.36.) All the intermedi- from a ruthenium complex with a N–H ligand had been pro- ates and transition states were located by full optimization posed (metal–ligand bifunctional mechanism, also called of the complex + substrate + water cluster species into the Noyori type mechanism, Scheme 5) [53]. AIMD calcula- solvent. Using this cluster-continuum model the transition tions showed that the reaction in solution follows a differ- states for the chloride by water, and water by formate ligand ent mechanism than in gas phase: the reaction is stepwise, exchanges, were located. The effect of the cluster size was hydride transfer happens first, and a methoxide intermediate assessed by increasing the number of water molecules to exists in solution for a short but finite time. five in selected steps (Fig. 9b). The barrier of the hydride A similar conclusion arises from a recent computational transfer with two more explicit solvent molecules increases study of hydrogenation of ketones with t-BuOK in tert- only 1.4 kcal ­mol−1, thus validating the employed model. butanol (Scheme 6) [54]. By analogy with Ru-catalyzed This reaction is an example of reactivity in protic solvent ketone hydrogenations (Scheme 5), experimentalists pro- media. In the next section we will go further on this issue. posed that ketone reduction proceeds through a six-mem- bered cyclic transition state involving the ketone, ­H2 and the alkoxide base [55]. A subsequent theoretical study of the 5 Reactions in Protic Solvents reaction in solution was performed through CSM calcula- tions on the optimized gas-phase structures [56], thus assum- According to the discussion presented in the precedent sec- ing the concerted transition state depicted in Scheme 6b. tions it is apparent that modeling catalytic processes in protic Static and dynamic DFT calculations (AIMD), including solvents, particularly alcohols and water, do demand inclu- explicit solvent molecules give a completely different pic- sion of explicit solvent molecules. Protic solvents provide a ture of the reaction, consisting of three steps: (i) cleavage strong hydrogen bonding environment that can influence the of the H–H bond to afford potassium hydride; (ii) addition relative stabilities of reactants, intermediates and transition of potassium hydride across the C=O bond of the ketone; states and products. In a protic solvent, solvent molecules (iii) H/K exchange to regenerate the catalyst and afford the can coordinate, can be directly involved as a proton shuttle product (Scheme 6c) [54]. in proton transfer steps, can generate small charged species Dehydrogenations in alcohol solvent are deeply influ- as hydroxide, proton or alkoxide and could compete with the enced by solvent molecules, similarly to the above com- catalyst as the proton source (i.e. metal complexes including mented hydrogenations. An AIMD study of the mecha- N–H ligands). These effects are not adequately accounted for nism of methanol dehydrogenation catalyzed by a highly with a continuum model. Even more important, this type of active PNP pincer ruthenium complex in methanol solvent solvent can drive the reaction towards new pathways, differ- unveiled this behavior [57]. In all the previously proposed ent from gas-phase or low polar solvents. mechanisms, protonation and deprotonation of the amido ligand were considered to be key steps of the catalytic 5.1 Alcohol Solvents cycle but relayed on a Noyori-type concerted mechanism (Scheme 5) [58, 59]. DFT-MD simulations with explicit In a pioneering work, Meijer and Handgraaf applied AIMD solvent evidenced the spontaneous protonation of the techniques to simulate the ruthenium-catalyzed intercon- ligand nitrogen atom of PNP-ruthenium complex by a sol- version of formaldehyde and methanol, as a model for vent molecule. Indeed, the amido moiety of the catalyst 13 130 Topics in Catalysis (2022) 65:118–140 Scheme 6  a Base-catalyzed hydrogenation of ketones. b Concerted transition state. c Revisited mechanism from static DFT and AIMD simulations, including solvent molecules remained protonated during the complete catalytic pro- Solvent can be crucial to stabilize reactive species. The cess, thus precluding a metal–ligand bifunctional mecha- Brönsted-acid catalyzed ring-opening of cyclopropanes nism. Protonation was accompanied by the coordination of takes place using triflic acid (HOTf) as the catalyst and hex- an additional ligand such as ­MeO− or MeOH to the metal afluoroisopropanol (HFIP) as a protic solvent [60]. The reac- center, in this way forming the resting state of the catalyst tion is initiated by protonation of cyclopropane by HOTf. (Scheme 7). Importantly, the nitrogen atom of the ligand However, without the inclusion of explicit solvation effects, was found to have a high proton affinity in polar, protic no protonated intermediate is found (no minimum in the solvents, indicating that ever under very basic conditions potential energy surface). Introducing three explicit HFIP the nitrogen atom should remain protonated [57]. molecules the proton of HTOf is transferred to the C=O Static DFT calculations including a small number of substituent of cyclopropane, enhancing the reactivity of explicit solvent molecules (cluster-continuum model) give cyclopropane ring (Scheme 8) [61]. a similar description of solvent effects in this reaction. The reliability of DFT calculations to elucidate reac- Using a more realistic solvent model which incorporates tion mechanism was challenged when Plata and Singleton explicit solvent molecules was critical to devise a new reported a detailed experimental study of the Morita–Bay- mechanistic view of this reaction. lis–Hillman (MBH) reaction and demonstrated that previ- ous computational studies using continuum solvent models 13 Topics in Catalysis (2022) 65:118–140 131 Scheme 7  Proposed mechanism for methanol oxidation by a PNP pincer ruthenium complex with explicit solvent [49] Scheme 8  Relative stabilities of neutral and protonated forms of cyclopropane-based reactant, without and with incorporation of HFIP solvent molecules failed to reproduce known mechanistic features of the reac- of the putative intermediates are either charged or zwitte- tion [62]. MBH reaction entails coupling of an activated rionic. Moreover, it takes place in protic solvents. Taking alkene derivative with an aldehyde to form a new C–C bond. into account the reaction’s scenario, it could be expected It is a multistep reaction in which while the reactants and that a very accurate description of solvation is required to products are neutral species with moderate polarity, many accurately reproduce by DFT calculations the main features 13 132 Topics in Catalysis (2022) 65:118–140 of the reaction mechanism. Indeed, three subsequent com- stabilizing the released C ­ l− anion. In the next step, H–H putational studies have validated this hypothesis. Harvey, bond breaking does not happen by the common oxidative Sunoj and coworkers employed high-level energy calcula- addition mechanism, but by an heterolytic splitting pathway tions combined with cluster-continuum calculations includ- in which the water solvent takes an active part, acting as ing methanol molecules, particularly important for species an external base by deprotonating the dihydrogen ligand. with partial charges, and molecular dynamics simulations Even in the presence of the more basic ­C5Me5 ligand the of the reacting system in a large box of explicit solvent [63]. heterolytic splitting is favored over the oxidative addition Basdogan and Keith used a stochastic filtering procedure [67]. The proton ends up either at a nitrogen atom of the with a global optimization protocol, based on static QM PTA ligand, yielding a monohydride with a protonated PTA calculations of micro-solvated clusters, to identify the low- ligand (for the complex with Cp’ = ­C5H5) or at the metal, energy geometries of the solute–solvent clusters. From the affording a dihydride (for the complex with Cp’ = ­C5Me5). globally optimized intermediates reaction pathways were Interestingly, when a small cluster of three water molecules calculated using the growing string method. They found was used in the cluster-continuum solvent model, only a con- that gradually adding more solvent molecules does not certed water-chain assisted deprotonation of the H ­ 2 ligand improve the agreement with experiment and that a micro- and proton transfer to one N atom of the PTA ligand was solvated model with five explicit methanol molecules was obtained (Fig. 10a) [66]. Introduction of an additional water the one in closer agreement with the experimental results molecule in the cluster allows the stabilization of a hydro- [64]. Sure et al. have devised a simple approach to address nium ion in the water medium than can protonate either of the main issues facing micro-solvation: when, where, and the basic centers (Fig. 10b) [67]. It seems that the dichotomy how many explicit solvent molecules have to be added to a between proton-shuttle and acid–base mechanisms [62] is solute, and they have successfully applied this protocol to highly dependent on the size of the solvent cluster in micro- the MBH reaction [65]. In this approach, that only employs solvated models. static DFT calculations, the proposed criterion to apply The C–H bond activation of methane in Shilov’s system micro-solvation is the thermodynamic preference for the over cis- and trans-PtCl2(H2O)(CH4) complexes was investi- association of a solvent (ΔGa < 0); the interaction sites are gated by ab initio molecular dynamics in explicit water [68]. found from screening charge densities. The three compu- Simulations within a realistic solvent representation revealed tational studies of the MBH reaction [63–65] have similar that the square pyramidal hydride complex ­PtHCl2(H2O) accuracies compared with the experimental results (absolute (CH3) formed in the oxidative addition step (Scheme 9b) is errors below 4 kcal ­mol−1 and MAD about 2 kcal ­mol−1) and not stable and spontaneously evolves into the square planar demonstrate that QM calculations can give a good descrip- methyl complex by the release of the proton to the bulk solu- tion of the MHB mechanism, provided that an accurate sol- tion (Scheme 9). Compared with previous modeling based vent modeling is employed. on minimal solvent clusters [69, 70], the explicit solvent modeling revealed that the platinum hydride is a transient 5.2 Water species with strong acid character, which spontaneously released a proton to the solution rather than coordinating Water medium can deeply affect the activation of strong a further water molecule to stabilize the metal center [68]. non-polar bonds, as H–H and C–H bonds. Computational C–H bond activation is a key step in the Ru-catalyzed studies with explicit solvent molecules have disclosed dehydrogenation of aqueous methanol. In this reaction, in active participation of solvent molecules in both processes, which three equivalents of ­H2 and ­CO2 are produced from resulting in reaction mechanisms different from those in ­CH3OH (Scheme 10a), C–H activation steps take place for non-polar solvents. Using a cluster-continuum modeling of methanol, methanediol and formic acid. Cluster-continuum the reaction media we unraveled the role of water in the static DFT and DFT-based molecular dynamics (DFT-MD) activation of dihydrogen by [Cp’Ru(PTA)2Cl] (Cp’ = ­C5H5, calculations have been performed to provide mechanistic ­C5Me5; PTA = 1,3,5-triaza-7-phosphaadamantane) water insights into this reaction [71]. The active form of the cata- soluble complexes [66, 67]. Water comes into play in all lyst (complex 2’) binds methanol and undergoes internal the steps of the reaction. The first step leads to the sub- proton transfer from the OH group of methanol to the amido stitution of an anionic C­ l− ligand by the weaker neutral moiety to form complex 2’A, that rearranges to form the ­H2 ligand. Initial dissociation of ­Cl− does not yield an zwitterionic agostic methoxide adduct 2’A’ (Scheme 10b). unsaturated [Cp’Ru(PTA)2]+ complex, but a monoaquo C-H bond breaking occurs from this intermediate. [Cp’Ru(PTA)2(H2O)]+ complex, detected by NMR. Subse- Constrained DFT-MD and static DFT micro-solvation quent displacement of the coordinated water by ­H2, results calculations using explicit solvent molecules yield very simi- in the complex {[Cp’Ru(PTA)2(η2-H2)]+. Thus, solvent con- lar results for this process: the solvent actively participates tributes to this step by coordinating to the metal as well as in the reaction, assisting the cleavage of the C–H bond by 13 Topics in Catalysis (2022) 65:118–140 133 Fig. 10  Heterolytic splitting of the H–H bond of tle mechanism) b Cluster of four water molecules: monohydride [CpRu(PTA)2(η2-H2)]+. a Cluster of three water molecules: tran- [Cp*Ru(PTA)2H] and solvated hydronium ion intermediate (acid– sition state of a concerted water-chain-assisted deprotonation of base mechanism). 3D-structures were generated using CYLview [36] the ­H2 ligand and proton transfer to the PTA ligand (proton shut- Scheme 9  Relevant species in the Shilov system from AIMD simulations in explicit water. a σ -complex; b Pt–H transient species; c square- planar methyl complex with a proton released to the bulk solvent means of hydrogen bonding interactions (Scheme 10c). In Conversely, the explicit solvent decreases the barrier for fact, an entirely uphill energy profile is obtained without the the formate dissociation because favors the formate to be in presence of explicit solvent molecules. Thus, the inclusion solution. Therefore, aqueous solvent provides a significant of explicit solvent is essential for the right modeling of the contribution to the reaction barriers. Moreover, in water the reaction. catalyst remains protonated, and a water molecule transfers The importance of hydrogen bond interactions between a proton instead of the catalyst to accomplish the hydrogena- the reacting system and the solvent has been also high- tion [73]. lighted in the AIMD study of Ru-catalyzed aqueous transfer A combined experimental and computational study of hydrogenation using formate as hydrogen donor. Explicit the ruthenium-catalyzed isomerization of allylic alcohols in aqueous solvent was included in the simulations (54 water aqueous medium (Fig. 11) disclosed another possible role molecules). Two key steps were analyzed: the hydride trans- that water solvent molecules can play: they not only can fer between formate and the catalyst and the dissociation coordinate but can act as cooperating ligand in metal–ligand of the ruthenium-formato complex [72]. Compared with bifunctional mechanisms [74]. a gas phase model, explicit aqueous solvent increases the During the mechanistic analysis cationic chloride free and barrier for the hydride transfer due to the stabilization of aqua species were detected in the ESI mass spectrum. Cal- the reactant complex by several hydrogen bonds between culations evidenced the accessibility of hydroxo-imidazole formate and solvent molecules. This stabilization is absent species in aqueous solution that can be the active species in the product because C ­ O2 does not form hydrogen bonds. in the isomerization process. The isomerization takes place 13 134 Topics in Catalysis (2022) 65:118–140 Scheme 10  a Elementary steps in the dehydrogenation of aqueous methanol to three equivalents of H ­ O2. b Relevant intermediates. c ­ 2 and C Micro-solvated transition state for the methanol C-H bond activation Fig. 11  a Ru-catalyzed isomerization of allylic alcohols in water. b Bis-allyl ruthenium(IV) catalysts. c Transition states for the hydrogen trans- fer steps. 3D-structures were generated using CYLview [36] through two sequential concerted hydrogen transfer steps, (acrolein hydrogenation: TS2, Fig. 11). Both steps disclose from the substrate to the catalyst (dehydrogenative step: an outer-sphere bifunctional catalyzed mechanism, with TS1, Fig. 11) and subsequent transfer back to the substrate a cooperative action of the metal and the hydroxo ligand 13 Topics in Catalysis (2022) 65:118–140 135 (Fig. 11c). Likewise to the previous dehydrogenation attack to the coordinated double bond (hydroxypalladation). described processes, this study also revealed that in water During these processes neighboring solvent molecules will the participation in the bifunctional catalysis of the N–H influence the reactivity of a water molecule, making manda- group of the azole ligand is not required, provided that a tory explicit solvent modeling. A proper description of the water ligand could play this role. Indeed, the related complex twofold role of water, and thus of the reaction mechanism, including a N-methylimidazole ligand, where both nitrogen can hardly be achieved by using static DFT calculations with atoms of the ring are substituted, shows a remarkable cata- a continuum or cluster-continuum solvent models. This has lytic activity in water [74]. been possible by using AIMD simulations with explicit solvent [75]. Two mechanisms had been postulated for the explanation of the experimental rate law: (a) anti nucleo- 6 Solvent as a Reagent philic attack by a water molecule from the bulk, and (b) syn nucleophilic attack by a coordinated hydroxide group In the previous section we have underlined the difficulties located cis to the alkene (Scheme 11). AIMD simulations, entailing modeling of reactions in protic solvents, and par- combined with rare event sampling techniques (as metady- ticularly water. This is even more challenging when water namics) were the approaches at choice to evaluate a long- is directly involved in the process as a reagent. The Wacker standing controversy regarding the Wacker’s mechanism process and O–O bond formation through water nucleophilic [75]. attack are both prototypical examples of reactions in which One of the questions revealed by the simulations was the water acts as a solvent and as a reagent. nature of the active species formed by a sequence of ligand The Wacker process is the aerobic oxidation of ethylene exchanges from the catalyst precursor ­[PdCl4]2 – (1). The to acetaldehyde catalyzed by palladium(II) chloride in aque- ligand exchanges take place in a stepwise fashion for the ous solution. The net reaction can be described as follows: formation of both cis and trans isomers, complexes 5 and [ ]2− 3, respectively (Scheme 12). Simulations show that the for- PdCl4 + C2 H4 + H2 O → CH3 CHO + Pd + 2 HCl + 2 Cl− mation of cis compound is uphill in Gibbs energy by about 12 kcal/mol, whereas the formation of trans intermediate is The Pd(II) catalyst is regenerated by oxidation of Pd(0) nearly thermoneutral on the Gibbs energy scale, by ~ 1 kcal/ by copper(II) reagent ­(CuCl2). In this reaction the water sol- mol [76, 77]. The Gibbs energy profiles obtained for the vent is initially involved in the ligand exchange that leads to formation of both isomers are kinetically feasible. However, the catalytically active species, and then in the nucleophilic the formation of the cis isomer is clearly thermodynamically unfavorable. Indeed, the energy difference between both iso- mers is equivalent to an equilibrium constant of ­10–9 for the trans to cis equilibrium; therefore, the concentration of cis isomer in solution should be negligible. According to the results obtained for the speciation of the ­[PdCl4]2− complex, the formation of the trans isomer is favored, thus supporting the outer sphere mechanism. Moreover, the outer sphere mechanism was also supported by the study of the nucleophilic addition step in two different reaction conditions, at low [78] and high ­Cl− concentrations Scheme 11  Postulated reaction pathways for the Wacker process: syn [76, 79]. Simulations performed under the first conditions, and anti addition of water to the Pd-coordinated ethene low ­Cl− concentration, were carried out on both cis- and Scheme 12  Ligand exchanges on the catalyst precursor [­ PdCl4]2− leading to the trans- and cis-[PdCl2(C2H4)(H2O)] species 13 136 Topics in Catalysis (2022) 65:118–140 solvent molecules always make calculations more com- plex. Moreover, both cluster-continuum and fully explicit solvent models have inherent limitations commented in the Introduction section. Thus, if a good description of the reac- tion medium is already achieved by means of a continuum solvent model alone, there is no need for explicit solvent. The purpose of this article is to give a warning about situa- tions which require a deeper analysis of solvent effects and that can demand going beyond the usual continuum solvent modeling. This is done highlighting selected examples from the literature, grouped by reaction scenarios, that help to Scheme 13  O–O bond formation pathway via water nucleophilic illustrate the reasons that urge to go beyond the continuum attack (WNA). M is a transition metal complex in high oxidation state models to the employment of micro-solvated of fully explicit solvent models. This section summarizes the examples analyzed. trans-[PdCl2(C2H4)(H2O)] complexes. For the trans isomer The success and limitations of CSM are clearly shown the obtained barrier for the hydroxypalladation step was by means of three different examples. One example per- compatible with experimental observations: the calculated forms a comparative analysis over non-catalyzed and Ag- barrier of 19 kcal ­mol−1 is slightly lower than the barrier catalyzed ring-closing reaction steps. The work shows that obtained experimentally, 22 kcal ­mol−1. Simulations start- pairwise interactions between the solute and the polar sol- ing from the cis-[PdCl2(C2H4)(H2O)] complex showed that vent molecules are significant, but they do not induce any isomerization to the trans complex and a subsequent outer kind of ordering in the solvation environment, therefore, the sphere nucleophilic addition is more feasible than the inner process can be safely described with a continuum model sphere mechanism. Moreover, metadynamics calculations that provides “averaged” solvent effects. A second exam- were also performed by allowing the inner and outer sphere ple analyses the decomposition of peroxyformate in many mechanisms take place simultaneously; simulations revealed different solvents. The main transition state, described as a that the outer sphere mechanism is again the preferred one. proton abstraction process, clearly shows that when the sol- Simulations at high ­Cl− concentration generated similar vent influence is due to its polarity, continuum models can conclusions regarding the inner/outer reaction mechanism capture solvent effects, showing the success of CSM. The comparison. third example, devoted to a metal catalyzed dehydrogenation Another appealing example of reaction in which water is process, shows that inclusion of solvent effects (by means of the reagent and the solvent is the O–O bond formation in the continuum and explicit solvent molecules) not only changes catalytic oxidation of water through the called water nucleo- the reaction mechanism from concerted to stepwise, but also philic attack mechanism (WNA). In the WNA mechanism modifies the overall view of its mechanism. This example the O–O bond is formed by a water attack on an oxo ligand demonstrates the importance of properly accounting for sol- (Scheme 13) [80]. vent effects in mechanistic studies. A series of recent papers, dealing with different catalysts, Another important aspect to pay attention when decid- point out the suitability of AIMD simulations with explicit ing the solvent model is the solvent capacity to coordinate solvent to tackle this challenging process [81–83]. reagents (especially when using metal catalysts); this is com- monly stated as solvent dynamics. Indeed, the association- dissociation process can affect the speciation of reactive 7 Summary and Conclusions species and they can modify the course of reaction. This is shown through several examples as the speciation of Negishi The use of continuum solvent models to study chemical pro- coupling reagents or the Schlenk equilibrium of Grignard cesses is a great success in the use of theoretical approaches reagents. In the former, the study concluded that species to model the chemistry performed in the lab. However, they ­ZnMe2(THF)2, ZnMeCl(THF)2, and Z ­ nCl2(THF)2 should cannot be applied to whatever situation. In many cases intro- be the ones considered in any computational modeling of duction of explicit solvent molecules in the QM description these zinc(II) species in THF solvent. The latter shows that of the system is mandatory. It is not in our aim, and we think the reaction is controlled by solvent dynamics: asymmetric it is not possible, to give a simple recipe to advise/propose solvation on Mg ions is required to promote the Mg–Cl and when explicit solvent molecules must be incorporated in the Mg–CH3 bond cleavage. The general advise is that when QM description of the system, along with the continuum good coordinating solvents are employed, they must be method or as a fully explicit model. Addition of explicit included in the model. For poorly coordinating solvents the 13 Topics in Catalysis (2022) 65:118–140 137 answer is not so clear, though the described example shows Intimately close to the previous examples are the reac- that “empty coordination sites” are not accessible even for tions where the solvent is a reagent in the process. The non-coordinating solvents. paradigmatic cases are those where the reaction takes place Another critical issue is related to the solvation of small in water solvent, and water molecules are also reactants in and charged species. As the described examples highlight, the process. The Wacker process and the water nucleophilic inclusion of explicit solvent molecules is required for a addition reaction were selected. In both cases the use of realistic modeling of processes involving such species. In AIMD procedures was employed (therefore using a fully one example, a protodeauration reaction step, computa- explicit solvent model) to gain a clearer view of the overall tions shows that solvation modifies the electrophilicity of chemical process. the protonating agent: it is overestimated by a pure con- Overall, the selected examples showed the successes tinuum model. The electrophilicity of the proton decreases and failures of the continuum model, and which cases the by increasing the number of solvent water molecules in the hybrid cluster-continuum model or even fully explicit mod- solvation sphere; this is reflected by the increase in the acti- els need to be considered. These models represent significant vation barrier for the protodeauration. Therefore, inclusion improvements in the computational modeling of catalytic of explicit solvent molecules is mandatory in these cases. systems but have inherent limitations. Certainly, proper In another example, related to the behavior of ion pairs, it description of solvation is a difficult and tricky issue that is shown that their description can be also highly affected need to be cautiously considered for the accurate computa- by the inclusion or not of explicit solvent molecules in the tional modeling of chemistry in solution. model. The example selected shows that the presence of explicit solvent molecules modifies the relative energy of Acknowledgements Financial support from Spanish Ministry of Sci- ence and Innovation (project PID-2020-116861GB-I00) is acknowl- the formed ion pair during the reaction pathway. Neglecting edged. Universitat Autònoma de Barcelona (UAB) is also acknowl- the incorporation of explicit solvent molecules could drive edged for a PIF grant to G.N. to a wrong conclusion about the mechanism. One of the paradigmatic cases where continuum models Funding Open Access Funding provided by Universitat Autonoma de need to be complemented are probably the reactions per- Barcelona. formed in protic solvents. Protic solvents provide a strong Open Access This article is licensed under a Creative Commons Attri- hydrogen bonding environment. In reaction profiles this can bution 4.0 International License, which permits use, sharing, adapta- largely influence the relative stabilities of reactants, inter- tion, distribution and reproduction in any medium or format, as long mediates, transition states and products. Moreover, protic as you give appropriate credit to the original author(s) and the source, solvents can also coordinate, or they can be directly involved provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are as proton shuttles in proton transfer steps. Note that protic included in the article's Creative Commons licence, unless indicated solvents can also generate small charged species as hydrox- otherwise in a credit line to the material. If material is not included in ide, proton or alkoxide, thus becoming small and charged the article's Creative Commons licence and your intended use is not species (earlier commented). Their behavior is generally not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a properly described by continuum models alone, therefore copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. becoming prominent candidates to employ cluster-contin- uum models. In addition, this type of solvents can drive the reaction towards new pathways, different from gas-phase or low polar solvents. References The transfer hydrogenation reaction in alcohol solvent was described as a paradigmatic case. A concerted transfer 1. Dyson PJ, Jessop PG (2016) Solvent effects in catalysis: rational of a proton and a hydride from the metal complex with a improvements of catalysts via manipulation of solvent interac- N–H ligand was initially proposed (metal–ligand bifunc- tions. Catal Sci Technol 6:3302–3316 2. Huxoll F, Jameel F, Bianga J, Seidensticker T, Stein M, Sadowski tional mechanism, also called Noyori type mechanism). G, Vogt D (2021) Solvent selection in homogeneous catalysis AIMD simulations disclosed that the reaction pathway in optimization of kinetics and reaction performance. ACS Catal solution was operating through a different mechanism than 11:590–594 in gas phase (even including CSM). The reaction is stepwise, 3. Varghese JJ, Mushrif SK (2019) Origins of complex solvent effects on chemical reactivity and computational tools to investi- hydride transfer happens first, and then the proton transfer; gate them: a review. React Chem Eng 4:165–206 moreover, the latter may not involve even the metal catalyst. 4. Basdogan Y, Maldonado AM, Keith JA (2020) Advances and chal- Several examples on how water solvent affect the activation lenges in modeling solvated reaction mechanisms for renewable of H–H or C–H are also commented. In all cases, the pres- fuels and chemicals. WIREs Comput Mol Sci 10:e1446 5. Luque FJ, Curutchet C, Muñoz-Muriedas J, Bidon-Chanal A, ence of explicit solvent molecules is mandatory to obtain Soteras I, Morreale A, Gelpí JL, Orozco M (2003) Continuum reliable results on the reaction mechanism. 13 138 Topics in Catalysis (2022) 65:118–140 solvation models: dissecting the free energy of solvation. Phys 28. Pincock RE (1964) Effects of nonpolar solvents on an ionic reac- Chem Chem Phys 5:3827–3836 tion. The ionic decomposition of tert-butylperoxy formate. J Am 6. Herbert JM (2021) Dielectric continuum methods for quantum Chem Soc 86:1820–1826 chemistry. WIREs Comput Mol Sci 11:e1519 29. Fehér PP, Stirling A (2019) Assessment of reactivities with 7. Mennucci B (2010) Continuum solvation models: what else can explicit and implicit sol-vent models: QM/MM and gas-phase we learn from them? J Phys Lett 1:1666–1674 evaluation of three different Ag-catalysed furan ring formation 8. Sunoj RB, Anand M (2012) Microsolvated transition state models routes. New J Chem 43:15706–15713 for improved insight into chemical properties and reaction mecha- 30. Zeifman AA, Novikov FN, Stroylov VS, Stroganov OV, Svitanko nisms. Phys Chem Chem Phys 14:12715–12736 IV, Chilov GG (2015) An explicit account of solvation is essential 9. Pliego JR Jr, Riveros JM (2020) Hybrid discrete-continuum solva- for modeling Suzuki-Miyaura coupling in protic solvents. Dalton tion methods. WIREs Comput Mol Sci 10:e1440 Trans 44:17795 10. Hadad C, Florez E, Acelas N, Merino G, Restreppo A (2019) 31. Norjmaa G, Maréchal J-D, Ujaque G (2019) Microsolvation and Microsolvation of small cations and anions. Int J Quant Chem encapsulation effects on supramolecular catalysis: C-C reduc- 119:e25766 tive elimination inside ­[Ga4L6]12– metallocage. J Am Chem Soc 11. Simm GN, Türtscher PL, Reiher M (2020) Systematic microsol- 141:13114–13123 vation approach with a cluster-continuum scheme and conforma- 32. Groenenboom MC, Keith JA (2016) Explicitly unraveling tional sampling. J Comput Chem 41:1144–1145 the roles of counterions, solvent molecules, and electron cor- 12. Basdogan Y, Groenenboom MC, Henderson E, De S, Rempe SB, relation in solution phase reaction pathways. J Phys Chem B Keith JA (2020) Machine learning-guided approach for studying 120:10797–10807 solvation environments. J Chem Theory Comput 16:633–642 33. Catalán J, Palomar J, Díaz C, de Paz JLG (1997) On sol- 13. Vidossich P, Lledós A, Ujaque G (2016) First-principles molecu- vent basicity: analysis of the SB scale. J Phys Chem A lar dynamics studies of organometallic complexes and homogene- 101:5183–5189 ous catalytic processes. Acc Chem Res 49:1271–1278 34. Díaz-Torres R, Alvarez S (2011) Coordinating ability of anions 14. Maldonado AM, Basdogan Y, Berryman JT, Rempe SB, Keith JA and solvents towards transition metals and lanthanides. Dalton (2020) First-principles modeling of chemistry in mixed solvents: Trans 40:10742–10750 where to go from here? J Chem Phys 152:130902 35. del Pozo J, Pérez-Iglesias M, Álvarez R, Lledós A, Casares JA, 15. Rossi K, Jurásková V, Wischert R, Garel L, Corminboeuf C, Espinet P (2017) Speciation of ZnMe2, ZnMeCl, and ZnCl2 in Ceriotti M (2020) Simulating solvation and acidity in complex tetrahydrofuran (THF), and its influence on mechanism calcula- mixtures with first-principles accuracy: the case of ­CH3SO3H and tions of catalytic processes. ACS Catal 7:3575–3583 ­H2O2 in phenol. J Chem Theory Comput 16:5139–5149 36. CYLview20 (2020) Legault CY. Université de Sherbrooke. http://​ 16. Eisenstein O, Ujaque G, Lledós A (2020) What makes a good www.​cylvi​ew.​org (computed) energy profile? Top Organomet Chem 67:1–38 37. Vidossich P, Ujaque G, Lledós A (2014) Palladium monophos- 17. Tomasi J, Mennucci B, Cammi R (2005) Quantum mechanical phine Pd(PPh3): is it really accessible in solution? Chem Commun continuum solvation models. Chem Rev 105:999–3093 50:661–663 18. Marenich AV, Cramer CJ, Truhlar DG (2009) Universal solvation 38. Planas F, Kohlhepp SV, Huang G, Mendoza A, Himo F (2021) model based on solute electron density and on a continuum model Computational and experimental study of turbo-organomagne- of the solvent defined by the bulk dielectric constant and atomic sium amide reagents: cubane aggregates as reactive intermediates surface tensions. J Phys Chem B 113:6378–6396 in Pummerer coupling. Chem Eur J 27:2767–2773 19. Klamt A (2011) The COSMO and COSMO-RS solvation models. 39. Peltzer RM, Eisenstein O, Nova A, Cascella M (2017) How sol- WIREs Comput Mol Sci 1:699–709 vent dynamics controls the Schlenk equilibrium of Grignard rea- 20. Lipparini F, Mennucci B (2016) Perspective: polarizable contin- gents: a computational study of ­CH3MgCl in tetrahydrofuran. J uum models for quantum-mechanical descriptions. J Chem Phys Phys Chem B 121:4226–4237 144:160901 40. Peltzer RM, Gauss J, Eisenstein O, Cascella M (2020) The Grig- 21. Ho J, Klamt A, Coote ML (2010) Comment on the correct use of nard reaction: unraveling a chemical puzzle. J Am Chem Soc continuum solvent models. J Phys Chem A 114:3442–13444 142:2984–2994 22. Lledós A (2021) Computational organometallic catalysis: where 41. Li Y, Lee T, Weerasiri K, Wang T, Buss EE, McKee ML, Gorden we are, where we are going. Eur J Inorg Chem 2021:2547–2555 AEV (2014) 2-Quinoxalinol diamine Cu(II) complex: facilitat- 23. Harvey JN, Himo F, Maseras F, Perrin L (2019) Scope and chal- ing catalytic oxidation through dual mechanisms. Dalton Trans lenge of computational methods for studying mechanism and 43:13578–13583 reactivity in homogeneous catalysis. ACS Catal 9:6803–6813 42. Le HM, Guagliardo M, Gorden AEV, Clark AE (2021) Ensemble 24. Maldonado AM, Hagiwara S, Choi TH, Eckert F, Schwarz K, effects on allylic oxidation within explicit solvation environments. Sundararaman R, Otani M, Keith JA (2021) Quantifying uncer- Dalton Trans 50:9259–9268 tainties in solvation procedures for modeling aqueous phase reac- 43. Lim C-H, Holder AM, Hynes JT, Musgrave CB (2017) Dihydrop- tion mechanisms. J Phys Chem A 125:154–164 teridine/pteridine as a 2H+/2e− redox mediator for the reduc- 25. Balaraman E, Khaskin E, Leitus G, Milstein D (2013) Catalytic tion of C ­ O2 to methanol: a computational study. J Phys Chem B transformation of alcohols to carboxylic acid salts and H ­ 2 using 121:4158–4167 water as the oxygen atom source. Nat Chem 5:122–125 44. Li J, Fisher CL, Chen JL, Bashford D, Noodleman L (1996) Cal- 26. Dub PA, Gordon JC (2017) Metal−ligand bifunctional cataly- culation of redox potentials and pKa values of hydrated transition sis: the “accepted” mechanism, the issue of concertedness, and metal cations by a combined density functional and continuum the function of the ligand in catalytic cycles involving hydrogen dielectric theory. Inorg Chem 35:4694–4702 atoms. ACS Catal 7:6635–6655 45. Uudsemaa M, Tamm T (2003) Density-functional theory calcula- 27. Garay-Ruíz F, Bo C (2021) Rationalizing the mechanism of per- tions of aqueous redox potentials of fourth-period transition met- oxyformate decomposition: computational insights to understand als. J Phys Chem A 107:9997–10003 solvent influence. Chem Eur J 27:11618–11626 46. Ortuño MA, Lledós A, Maseras F, Ujaque G (2014) The trans- metalation process in Suzuki-Miyaura reactions: calculations 13 Topics in Catalysis (2022) 65:118–140 139 indicate lower barrier via boronate intermediate. ChemCatChem 64. Basdogan Y, Keith JA (2018) A paramedic treatment for mod- 6:3132–3138 eling explicitly solvated chemical reaction mechanisms. Chem 47. Carrow BP, Hartwig JF (2011) Distinguishing between pathways Sci 9:5341–5346 for transmetalation in Suzuki−Miyaura reactions. J Am Chem Soc 65. Sure R, el Mahdali M, Plajer A, Deglmann P (2021) Towards 133:2116–2119 a converged strategy for including microsolvation in reaction 48. BabaAhmadi R, Ghanbari P, Rajabi NA, Hashmi ASK, Yates mechanism calculations. J Comput Aided Mol Des 35:473–492 BF, Ariafard A (2015) A theoretical study on the protodeauration 66. Rossin A, Gonsalvi L, Phillips AD, Maresca O, Lledós A, Peru- step of the gold(I)-catalyzed organic reactions. Organometallics zzini M (2007) Water-assisted H-H bond splitting mediated by 34:3186–3195 [CpRu(PTA)2Cl] (PTA = 1,3,5-triaza-7-phosphaadamantane). A 49. Jiménez-Tenorio M, Puerta MC, Valerga P, Ortuño MA, Ujaque DFT analysis. Organometallics 26:3289–3296 G, Lledós A (2013) Counteranion and solvent assistance in ruthe- 67. Kovács G, Rossin A, Gonsalvi L, Lledós A, Peruzzini M nium-mediated alkyne to vinylidene isomerizations. Inorg Chem (2010) Comparative DFT analysis of ligand and solvent effects 52:8919–8932 on the mechanism of ­H2 activation in water mediated by half- 50. Abril PM, del Río MP, López JA, Lledós A, Ciriano MA, Tejel sandwich Complexes [Cp’Ru(PTA)2Cl] (Cp’ =C5H5, ­C5Me5; C (2019) Inner-sphere oxygen activation promoting outer-sphere PTA = 1,3,5-triaza-7-phosphaadamantane). Organometallics nucleophilic attack on olefins. Chem Eur J 25:14546–14554 29:5121–5131 51. Carrión MC, Ruiz-Castañeda M, Espino G, Aliende C, Santos L, 68. Vidossich P, Ujaque G, Lledós A (2012) PtII as a proton shuttle Rodríguez AM, Manzano BR, Jalón FA, Lledós A (2014) Selec- during C-H bond activation in the Shilov process. Chem Commun tive catalytic deuterium labeling of alcohols during a transfer 48:1979 hydrogenation process of ketones using D ­ 2O as the only deuterium 69. Siegbahn PEM, Crabtree RH (1996) Modeling the solvent source. Theoretical and experimental demonstration of a Ru−H/ sphere: mechanism of the Shilov reaction. J Am Chem Soc D+ exchange as the key step. ACS Catal 4:1040–1053 118:4442–4450 52. Handgraaf J-W, Meijer EJ (2007) Realistic modeling of ruthe- 70. Zhu H, Ziegler T (2006) A theoretical study of the original Shilov nium-catalyzed transfer hydrogenation. J Am Chem Soc reaction involving methane activation by platinum tetrachlo- 129:3099–3103 ride ­(PtCl42-) in an acidic aqueous solution. J Organomet Chem 53. Yamakawa M, Ito H, Noyori R (2000) The metal−ligand bifunc- 691:4486–4497 tional catalysis: a theoretical study on the ruthenium(II)-catalyzed 71. Govindarajan N, Sinha V, Trincado M, Grützmacher H, Meijer EJ, hydrogen transfer between alcohols and carbonyl compounds. J de Bruin B (2020) An in-depth mechanistic study of Ru-catalyzed Am Chem Soc 122:1466–1478 aqueous methanol dehydrogenation and prospects for future cata- 54. Dub PA, Tkachenko NV (2021) Mechanism of potassium tert- lyst design. ChemCatChem 12:2610–2621 butoxide-catalyzed ketones hydrogenation in the solution phase. 72. Pavlova A, Rösler E, Meijer EJ (2016) Mechanistic aspects of J Phys Chem A 125:5726–5737 using formate as a hydrogen donor in aqueous transfer hydrogena- 55. Berkessel A, Schubert TJS, Müller TN (2002) Hydrogena- tion. ACS Catal 6:5350–5358 tion without a transition metal catalyst: on the mechanism of 73. Pavlova A, Meijer EJ (2012) Understanding the role of water in the base-catalyzed hydrogenation of ketones. J Am Chem Soc aqueous ruthenium-catalyzed transfer hydrogenation of ketones. 124:8693–8698 ChemPhysChem 13:3492–3496 56. Chan B, Radom L (2005) Base-catalyzed hydrogenation: rational- 74. Díez J, Gimeno J, Lledós A, Suárez FJ, Vicent C (2012) Imidazole izing the effects of catalyst and substrate structures and solvation. based ruthenium(IV) complexes as highly efficient bifunctional J Am Chem Soc 127:2443–2454 catalysts for the redox isomerization of allylic alcohols in aqueous 57. Sinha V, Govindarajan N, de Bruin B, Meijer EJ (2018) How sol- medium: water as cooperating ligand. ACS Catal 2:2087–2099 vent affect C-H activation and hydrogen production pathways in 75. Stirling A, Nair NN, Lledós A, Ujaque G (2014) Challenges in homogeneous Ru-catalyzed methanol dehydrogenation reactions. modelling homogeneous catalysis: new answers from ab initio ACS Catal 8:6908–6913 molecular dynamics to the controversy over the Wacker process. 58. Alberico E, Lennox AJJ, Vogt LK, Jiao H, Baumann W, Drexler Chem Soc Rev 43:4940–4952 H-J, Nielsen M, Spannenberg A, Checinski MP, Junge H, Beller 76. Nair NN (2011) Ligand exchanges and hydroxypalladation reac- M (2016) Unravelling the mechanism of basic aqueous methanol tions of the Wacker process in aqueous solution at high C ­ l− con- dehydrogenation catalyzed by Ru-PNP pincer complexes. J Am centration. J Phys Chem B 115:2312–2321 Chem Soc 138:14890–14904 77. Kovács G, Stirling A, Lledós A, Ujaque G (2012) The nature of 59. Yang X (2014) Mechanistic insights into ruthenium-catalyzed pro- [PdCl2 (C2H4)(H2O)] as an active species in the Wacker pro- duction of H ­ 2 and C­ O2 from methanol and water: a DFT study. cess: new insights from ab initio molecular dynamics simulations. ACS Catal 4:1129–1133 Chem Eur J 18:5612–5619 60. Richmond E, Vuković VD, Moran J (2018) Nucleophilic ring 78. Comas-Vives A, Stirling A, Lledós A, Ujaque G (2010) The opening of donor–acceptor cyclopropanes catalyzed by a Brønsted Wacker process: inner- or outer-sphere nucleophilic addition? acid in hexafluoroisopropanol. Org Lett 20:74–577 New insights from ab initio molecular dynamics. Chem Eur J 61. Zhou Y, Xue R-C, Feng Y, Zhang L (2020) How does HOTf/ 16:8738–8747 HFIP cooperative system catalyze the ring-opening reaction of 79. Imandi V, Kunnikuruvan S, Nair NN (2013) Hydroxypalladation Cyclopropanes? A DFT study. Asian J Org Chem 9:311–316 precedes the rate-determining step in the Wacker oxidation of 62. Plata RE, Singleton DA (2015) A case study of the mecha- ethene. Chem Eur J 19:4724–4731 nism of alcohol-mediated Morita Baylis-Hillman reactions. 80. Sala X, Maji S, Bofill R, García-Antón J, Escriche L, Llobet A The importance of experimental observations. J Am Chem Soc (2014) Molecular water oxidation mechanisms followed by transi- 137:3811–3826 tion metals: state of the art. Acc Chem Res 47:504–516 63. Liu Z, Patel C, Harvey JN, Sunoj RB (2017) Mechanism and 81. Hodel FH, Luber S (2016) What influences the water oxidation reactivity in the Morita–Baylis–Hillman reaction: the challenge of activity of a bioinspired molecular ­CoII4O4 cubane? An in-depth accurate computations. Phys Chem Chem Phys 19:30647–30657 exploration of catalytic pathways. ACS Catal 6:1505–1517 13 140 Topics in Catalysis (2022) 65:118–140 82. Bernasconi L, Kazaryan A, Belanzoni P, Baerends EJ (2017) Publisher's Note Springer Nature remains neutral with regard to Catalytic oxidation of water with high-spin iron(IV)−oxo spe- jurisdictional claims in published maps and institutional affiliations. cies: role of the water solvent. ACS Catal 7:4018–4025 83. Govindarajan N, Tiwari A, Ensing B, Meijer EJ (2018) Impact of the ligand flexibility and solvent on the O−O bond formation step in a highly active ruthenium water oxidation catalyst. Inorg Chem 57:13063–13066 13