Beyond Continuum Solvent Models in Computational Homogeneous Catalysis

Abstract

In homogeneous catalysis solvent is an inherent part of the catalytic system. As such, it must be considered in the computational 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 systematizing 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.

1 Introduction

In homogenous catalysis the solvent is an inherent part of the catalytic system. Although for years it had been considered just as a passive medium where reaction events take place, its critical role on catalyzed reactions is nowadays well recognized [1]. Indeed, selection of the best solvent to carry out one reaction is crucial for the reaction performance and is a mandatory step in the optimization of reaction conditions [2]. At the macroscopic level, the solvent influences the solubility of the species present along the reaction pathway. At the molecular level, solvent molecules interact directly with substrates, products and catalyst molecules, as well as intermediates and transition states. The three main solvent properties accounting for its interactions with the reaction system are polarity (related with the dielectric constant ε) and Brönsted and Lewis acidity and basicity. These interactions can have a deep impact on reaction rate and selectivity. With this in mind it is clear that solvent must be considered when discussing and computing reaction mechanisms.

Three main approaches have been developed to deal with solvent effects in computational studies of reaction mechanisms using quantum chemistry schemes: implicit (or continuum), cluster/continuum (or hybrid implicit/explicit) and fully explicit solvation models (Fig. 1). There is a vast literature on this subject and here we will only briefly recall the main characteristics of each approach. Recent articles give a comparative overview of the three approaches [3, 4].

Fig. 1

Schematic representation of solvent models

In continuum solvent models (CSM) no explicit solvent molecules are introduced in the chemical system to be computed. In this implicit approach the solvent is a homogeneous polarizable medium, described by its dielectric constant ε. The solute is placed in a cavity that is created in this medium. Electron density of the solute inside the cavity polarizes the dielectric, continuum, which in turns acts on the solute, modifying its electron distribution. In this way, the electrostatic contribution to the solvation Gibbs energy is evaluated. Evaluation of non-electrostatic contributions, to be added to the electrostatic one, is required to obtain accurate solvation energies [5]. Differences in the theoretical foundations, the way of defining the cavity and the boundary conditions give rise to various implicit solvation models, as the original Polarizable Continuum Model (PCM), the Conductor-like Screening Model (COSMO) and the Solvation Model based on Density (SMD) among others. A thorough review on the current state-of-the art of CSM have been recently published [6]. These methods are computationally efficient, keep the system size as in the gas phase and avoid sampling over solvent degrees of freedom. With these advantages, continuum solvent models have become the most employed approach in computational homogeneous catalysis. Very efficiently they combine an accurate, quantum mechanical (QM) description of molecular systems with simple but physically rooted description of environment [7]. However, they provide “averaged” solvent effects and do not capture specific interactions with solvent molecules, which in some cases may be important. It is not an appropriate model if solvent molecules can directly participate in the reaction.

These limitations can be overcome, in principle, by including a selected number of explicit solvent molecules into the atomistic part of the calculation. This is the main feature of the cluster-continuum approach, that is also referred as “micro-solvation” [8]. The explicit solvent molecules are treated in the same theoretical framework as the rest of the catalytic system. This “selected” explicit solvation is usually supplemented with the addition of implicit continuum solvation to account for the long-range electrostatic effects, thus becoming a hybrid implicit–explicit solvation scheme. A critical discussion of several hybrid discrete-continuum models has been recently published [9]. Appealing as it is, this cost-effective approach has a main shortcoming: the difficulty in defining the number and location of the explicit solvent molecules. Moreover, the conformational space strongly increases when adding several solvent molecules, making it necessary extensive conformational sampling [10]. Recently automated methods to determine how many explicit solvent molecules need to be added in hybrid cluster-continuum schemes to capture most of the interaction between the solute and the environment have been proposed [11, 12].

A realistic description of the condensed phase environment in which catalytic reactions occur is achieved with fully explicit solvent models: the solute is placed in a box solvated by as many solvent molecules as needed to reproduce the density of the system. These computational models can simulate the dynamics of the solvent, but to sample the vast configurational space of solvent configurations statistical methods are required, for instance within the ab initio molecular dynamics framework (AIMD) [13]. These calculations are generally tremendously demanding in computational resources. Moreover, acceleration techniques as metadynamics or umbrella sampling must be employed to access time scales at which reactive events (“rare events”) happen. It should be pointed out that despite the tremendous progress achieved in modeling chemical events in solution, modeling of mixed solvents is still very challenging [14, 15].

As commented above, nowadays practically all the computational studies of reaction mechanisms in homogenous catalysis performed with Density Functional Theory (DFT) calculations consider the solvent, mostly by means of continuum models. This methodology has been very useful to approach the computed model of the system to the chemical “real world” [16]. Thus, according to the success of plain continuum models, sometimes it is believed that they can be applied to any situation. The main goal of this chapter is to showcase, presenting selected examples, that there are situations in which proper modeling of the catalytic system demands going beyond of implicit solvation; in these cases, the inclusion of explicit solvent molecules in the quantum mechanical framework is mandatory for a proper description of the system. The article is not aimed to collect an exhaustive list of calculations including explicit solvent, but to highlight with examples the reasons behind the choice of micro-solvated or fully explicit solvent models, and the benefits provided by this choice. In this way we would like to contribute to set the limits of continuum solvent models in reactivity studies. With the aid of examples, we will start showing the success and limitations of CSM. Then, we will address first issues related with solvent coordination and solvent dynamics. This section will be followed by an overview of solvent effects in reactions involving small, charged species. Next section will be devoted to reactions in protic solvents and finally we will focus on the role of solvent as a reagent itself.

2 Success and Limitations of Continuum Solvent Models

The success of continuum models in modeling homogeneous catalysis processes is undeniable and has been already demonstrated on numerous cases [17,18,19,20]. Nowadays its use has become mandatory when computing energy profiles, even for reactions in solvents with low dielectric constant. The current models have been carefully parameterized to reproduce known experimental solvation Gibbs energies. When they are correctly used [21], these methods, coupled with accurate methods for computing the potential energy (DFT methods including dispersion effects or highly-correlated post-Hartree–Fock approaches as CCSD(T) or DLPNO/CCSD(T)), deliver accurate electronic energies in solution E_solv, that include the Gibbs energies of solvation. Adding to this thermodynamically “non-pure” magnitude, E_solv, thermal and entropic corrections for the solute, solution-phase Gibbs energies, suitable to build energy profiles, are obtained. It is important to stress that, as we [16, 22] and others [23, 24] have pointed out, proper use of CSM for modeling reaction mechanisms in solution implies optimization and characterization of all the structures in solution. Only in the cases in which optimization in solution fails due to technical issues, gas-phase optimization followed by single-point energy calculation in solution is justified. Saving computer time cannot be a reason to skip optimization in solution, regrettably a still common practice. Not doing that can lead to misconceptions in reaction mechanisms.

Dehydrogenation of a gem-diol to form carboxylic acid is a proposed step in the ruthenium-catalyzed transformation of primary alcohols to carboxylic acid salts and H₂ using water as the -OH group source. The reaction was carried out in basic aqueous solution [25].

From calculations in the gas-phase the dehydrogenation of the gem-diol proceeds over a concerted transition state where the proton and hydride transfer simultaneously (Fig. 2). Therefore, this would be the conclusion about the reaction mechanism if the solvent is introduced by means of single-point calculations in a continuum medium at the gas-phase optimized geometries. However, optimizing in continuum water it is found that the same reaction proceeds via two transition states, the first one describing hydride transfer and the second one proton transfer. An ion-pair intermediate connects both transition states [26]. When eight explicit water molecules are further added to the computations, the first step (hydride transfer) proceeds similarly, but with a barrier lowered about 5 kcal mol⁻¹. Importantly, in the second step (proton transfer) the source of the proton is not the C-H functionality, but the solvent [26]. Inclusion of solvent effects not only changes the reaction mechanism from concerted to stepwise but changes the overall view of its mechanism. This example demonstrates the importance of properly accounting for solvent effects in mechanistic studies.

Fig. 2

Ru-catalayzed dehydrogenation of gem-diol into butyric acid computed with different solvent models [26]

When correctly applied to suitable reactions, continuum solvent models provide an accurate description of reactions in solution. A recent study of the mechanism of peroxyformate decomposition (Scheme 1a) illustrates the potential of CSM, when used correctly, to differentiate solvent effects in a reaction strongly affected by solvent polarity [27].

Scheme 1

a Peroxyformate decomposition reaction. b Proposed transition state

Experimental kinetic constants for this reaction are available for twenty solvents and vary by a factor of 140 [28], making it an excellent case for benchmarking. The transition state (Scheme 1b) shows an important degree of charge separation and its relative stability with respect to the reactant is expected to be deeply affected by the polarity of the solvent. DFT calculations, performed with the SMD continuum model in a set of 29 solvents (from heptane, ε = 1.91, to formamide, ε = 108.94) give Gibbs energy barriers in very good agreement with the experimental values (RMSD value of 2.4 kcal mol⁻¹ for the DLPNO/CCSD(T) energies with implicit solvent). Moreover, optimization in solvent reveals subtle changes in the TS structure due to the solvent polarity, that can be correlated with the barriers [27]. This work clearly shows that when the solvent influence is due to its polarity continuum models can capture solvent effects.

Computational studies aimed to compare implicit and explicit solvent models in reaction mechanisms are very relevant to delineate the limits of continuum models and to determine the origin of their shortfalls. Calculations with explicit and implicit solvent models have been compared on non-catalyzed and Ag-catalyzed ring-closing reaction steps in a polar, aprotic solvent (N,N-dimethyl formamide, DMF, ε = 37.22) [29]. The reaction routes compared for the furan-ring formation have different charge states (Scheme 2).

Scheme 2

Ring-closing reactions computed in DMF with implicit and explicit solvent models [29]

The main conclusion of this study is that for this type of reactions the effect of the DMF solvent can be adequately described by continuum solvent models: barriers computed without and with explicit solvent molecules are similar. Analysis of trajectories obtained with molecular dynamics simulations with explicit solvent indicate that despite pairwise interactions between the solute and the polar solvent molecules are significant, they do not induce any kind of ordering in the solvation environment. Solvent molecules are collectively fluctuating around the solute molecule. This mobile, fluctuating solvent shell can be efficiently substituted by implicit solvent models with a reduced computational effort [29].

The Gibbs energy profile for the reductive elimination step of Suzuki–Miyaura couplings have been computed using explicit and implicit solvation approaches in a diverse set of solvents (benzene, toluene, DMF, ethanol and water) [30]. Although a good agreement between both approaches was obtained for aprotic solvents, the correlation was poor for protic solvents. The difference was attributed to the formation of Pd···HO hydrogen bonds between Pd(PPh₃)₂ and protic solvent molecules. This work is a warning about the use of implicit solvation for modeling reactions in protic solvents. The reductive elimination from a Au(III) complex to form a C–C bond in MeOH was shown to be also highly affected by micro-solvation [31]. Similar conclusions, for aqueous phase reactions, were inferred from the comparison of different solvation models for the CO₂ reduction by NaBH₄ in aqueous solution [32].

3 Solvent Coordination and Solvent Dynamics

Solvent dynamics can affect the course of a reaction despite the solvent molecules may not be involved in the general stoichiometric or catalytic process. One of the simplest ways to exhibit this behavior is by solvent coordination to a metal catalyst. Solvent molecules do move and, if they have basic character, they can coordinate to Lewis acid centers, as unsaturated transition metals [33]. Such a coordination, that is a particular case of association-dissociation process in a condensed phase, affects the speciation of reactive species and can impact the reaction mechanism. From the point of view of modeling this possibility raises the question of whether coordinated solvent molecules should be included in the QM description of the system. For the case of good coordinating solvents (e.g., pyridine, dimethyl sulfoxide (DMSO), formamide, etc.), generally the answer is yes, but the answer is not as easy for weaker coordinating solvents. A scale that attempts to quantify the weakly coordinating character of a variety of solvents has been proposed and a coordinating ability index has been calculated from structural data base searches [34]. This index gives a first indication of the feasibility of solvent coordination (positive values means coordinating solvent and negative non-coordinating ones). Those solvents with a coordinating ability index close to zero are defined as “ambiguous coordinator”, with similar possibilities to appear coordinated and uncoordinated. For instance, this is the case for acetonitrile and tetrahydrofuran (THF) [34]. Moreover, solvent coordination arises from an equilibrium between enthalpic and entropic terms, with the latter being quite difficult to evaluate.

We analyzed the speciation of Negishi coupling reagents, as ZnMe₂, ZnMeCl, and ZnCl₂ in THF solution, by experimental (infrared (IR) and calorimetric experiments) and DFT calculations [35]. IR observation suggests the metal center undergoes the THF coordination equilibria depicted in Fig. 3. Calculations sustain this proposal: the coordination of a second THF is clearly preferred over the coordination of a single THF molecule to ZnMe₂ (Scheme 3).

Fig. 3

a THF coordination equilibria suggested by IR measurements. b Optimized structures [35]. 3D-structures were generated using CYLview [36]

Scheme 3

DFT calculated ΔG values (kcal mol⁻¹) for reactions 1 and 2 in THF solvent. n is the number of solvent molecules in the solvent cluster

The solvent was modeled using a cluster/continuum approach employing (THF)_n clusters. In this way for n ≥ 3, the total number of components (two) remains constant at both sides of the equilibria, as shown in Scheme 3, eliminating a major contribution to entropy changes upon coordination. The study was also performed to ZnMeCl and ZnCl₂; in all cases (THF)₄ appears a good enough representation of the solvent. The conclusion of this study was that any computational modeling in THF should consider the species ZnMe₂(THF)₂, ZnMeCl(THF)_2, and ZnCl₂(THF)₂ as the ones participating in reactivity [35].

Solvent coordination can happen even with solvents with weak coordinative properties, as toluene. We investigated Pd(PPh₃) speciation in toluene (Tol) by means of explicit solvent simulations, using AIMD techniques [37]. Calculations provided computational evidence for the formation of [Pd(PPh₃)]-solvent adducts, as (Tol)Pd(PPh₃)₂, (Tol)Pd(PPh₃) and (Tol)₂Pd(PPh₃). An important conclusion from this study was that the bare Pd(PPh₃) often proposed in catalytic cycles, is actually not accessible in solution. Therefore, solvent coordination should be considered when modeling the reactivity of Pd(PR₃) complexes in solution.

Solvent coordination issues are even more sensitive when the reaction involves main group alkaline and alkaline-earth metal cations, without clear coordinative preferences and with a mainly electrostatic interaction with solvent molecules. Extensive investigation of the most stable aggregates formed in these systems is mandatory. For instance, in a recent study of the Pummerer coupling between sulfoxides and turbo-organomagnesium amides (Scheme 4), several structural possibilities, including explicit THF molecules, have been checked. The investigation was able to determine the most stable species that Grignard reagent (isopropylmagnesium chloride) and the magnesium amide (“R-MgX” + DIPAMgCl·LiCL) can form prior to the addition of the sulfoxide [38].

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]

Grignard reagents, with nominal formula RMgX, are seemingly simple organometallic species, however unraveling its speciation in THF solution and obtaining a detailed understanding of the mechanism by which they operate in the Grignard reaction have pushed solvent models to the limit [39, 40]. This is because this system constitutes a prototypical example of a reaction governed by solvent dynamics. AIMD simulations and enhanced-sampling methods in explicit solvent made possible solving this longstanding chemical puzzle.

The transformation of Grignard reagent CH₃MgCl into MgCl₂ and Mg(CH₃)₂ in tetrahydrofuran solvent (Fig. 4a) was characterized at the molecular level by means of AIMD simulations with enhanced-sampling metadynamics, using the simulation box displayed in Fig. 4c containing 42 THF molecules as a solvent model [39]. The interchange of the methyl and chloride groups between two magnesium centers is associated with a change in the solvation number. Asymmetric solvation on Mg ions is needed to promote the Mg-Cl and Mg-CH₃ bond cleavage. Thus, solvent dynamics is driving the ligand exchange in the Schlenk equilibrium.

Fig. 4

a Schlenk equilibrium of Grignard reagents. b Grignard reaction. c Simulation box used in the AIMD study (reproduced from [33])

Using similar methodology than for the Schlenk equilibrium, two competing mechanisms, the polar (or nucleophilic) implying heterolytic Mg-C bond breaking, and the radical, involving homolytic Mg-C bond breaking have been computed for the Grignard reaction (Fig. 4b) [40]. Both mechanisms are strongly dependent on the coordination of the Mg center. The solvent plays an essential role in the reaction: the more solvated the Mg species are, the more reactive become. Thus, as for the Schlenk equilibrium, solvent dynamics is essential for representing the energy profile.

Inclusion of explicit solvent can impact the mechanistic view of a reaction, even not governing the reaction course. This has been the case for the allylic oxidation of cyclohexene in acetonitrile (ε = 35.69) using a Cu(II) complex. Static DFT calculations of the reaction mechanism with a continuum description of acetonitrile resulted in an estimation of high barriers: 43.7 kcal mol⁻¹ for the H-abstraction step and 37.5 for the step in which cyclohexanone product is released, respectively [41]. Umbrella-sampling DFT-based molecular dynamics (DFT-MD) simulations in a periodic box containing 60 acetonitrile molecules opened up a new perspective on the mechanism for this reaction, identifying new reaction branches and side processes. Explicit solvent environment of the simulation stabilizes reactive species that open two new reaction pathways involving dehydrogenation and re-hydrogenation of the -NH₂ group bound to the Cu-site, with significantly lowered barriers relative to the prior static DFT studies with implicit solvent [42].

4 Solvation of Small, Charged Reagents

It is very common in homogenous catalysis processes to find steps in which participate small and charged reagents, as for instance OH⁻, HCOO⁻, Cl⁻ or H₃O⁺. Charged species are also often engaged in ligand exchanges. As the following examples highlight, inclusion of explicit solvent molecules is required for a realistic modeling of such processes.

A simple example illustrates the failure of continuum solvent models to describe the stabilization of small anionic species by a polar solvent: calculation of pK_a of HCOOH in water with such a model gives a value of 11.9, very far from the experimental value of 3.9. This failure is mainly due to the absence of explicit H-bonding by solvating H₂O molecules that stabilize the anionic conjugate base (HCOO⁻). A cluster-continuum model with two explicit water molecules H-bonding the anion only slightly improves the pK_a value (9.8). However, enlarging anion micro-solvation to seven water molecules yields a pK_a of 4.4, closer to the experimental value. For this reason, seven explicit water molecules were included in a computational study of the reduction of CO₂ to methanol by dihydropteridine [43]. Similar conclusions have been obtained for the calculation of redox potentials [44, 45].

In a computational study of the transmetalation step in the Suzuki–Miyaura cross-coupling reaction we analyzed the bromide-by-hydroxide replacement leading to a palladium hydroxo complex [46]. Without including explicit solvent molecules, the transition state for this substitution is located at 32.5 kcal mol⁻¹ 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 by the inaccurate representation of the base as a naked OH⁻. The reaction is performed in THF solvent, but the presence of water in the experiments is highly likely: they are coming from hydration molecules of the base itself, or for instance, from the solvent when using mixtures as THF/water [47]. The inclusion of water molecules solvating the two exchanging anions significantly decreases the barrier. The influence in the barrier of the number of water molecules considered (n) was evaluated. The series was considered converged at n = 3 because the change in the Gibbs energy profile does not significantly varies with respect to n = 4 (24.2 and 22.9 kcal mol⁻¹ for n = 3 and 4, respectively) (Fig. 5b).

Fig. 5

Transition state for the bromide-by-hydroxide exchange in the Pd(PPh₃)₂ coordination sphere. a Without explicit water molecules. b With four explicit water molecules. Main distances in Å and Gibbs energy barriers in kcal mol⁻¹. 3D-structures were generated using CYLview [36]

Protodemetallation is another process highly influenced by the presence of solvent. We remark here protodeauration, the formal exchange of Au(I) by a proton as a prototypical case. It is a very common final step in gold(I)-catalyzed reactions of unsaturated C–C bonds. The barrier is very sensitive to the solvation of the electrophile. For example, the protodeauration of an Au(PMe₃)(alkenyl) compound by a proton solvated by one, two, three, and four water molecules was investigated by DFT calculations (structure 4 in Fig. 6). The authors found that the activation energy increases when the number of the water molecules solvating the proton increases (Fig. 6) [48]. Solvation is modifying the electrophilicity of the protonating agent, overestimated using a purely continuum solvation model. When the number of water molecules solvating the proton increases, the electrophilicity of the proton decreases, as it is reflected by an increase in the LUMO energy of the protonating agent. It follows that the less electrophilic the protonating agent, the higher is the activation energy for protodeauration.

Fig. 6

Energy profile calculated for protodeauration of the alkenyl complex 4 with an increasing number of water molecules solvating the H₃O⁺ protonating agent [48]

It must be stressed that explicit solvation not only modifies acid–base properties of species as OH⁻ or H₃O⁺, but also affects weak basic anions, as a chloride. In a combined experimental/computational investigation we disclosed counteranion assistance in ruthenium-mediated alkyne to vinylidene isomerizations in methanol [49]. The theoretical study outlined the feasibility of hydrido − alkynyl species as intermediates, in which the hydride has an acidic character. Due to the acidic character of this intermediate, basic species as a chloride anion can act as a proton shuttle, taking the hydrogen from the metal to place it at the β-carbon, forming the vinylidene product. The basic character shown by the chloride anion decreases when one solvent methanol molecule is located in its direct environment. Accordingly, the Gibbs energy barrier for migration increases from 19.5 to 23.2 kcal mol⁻¹. Transition states in both models differ considerably (Fig. 7).

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⁻¹. 3D-structures were generated using CYLview [36]

The description of ion pairs can be also highly affected by the presence of explicit solvent molecules in the model. Without inclusion of explicit solvent molecules, ion pairs can be high energy species that could be wrongly discarded in reaction pathways. For instance, complex [Ir(pic)₂(MeOC₈H₁₂)] (pic = picolinate ligand) undergoes a fast isomerization between 2a and 2b structures (Fig. 8) [50]. A plausible pathway for this isomerization would entail the cleavage of the C–OCH₃ bond to form the ionic pair [Ir(pic)₂(cod)][OCH₃]. Using a continuum solvent model such ion pair is found 28.4 kcal mol⁻¹ above 2a, a value that suggests discarding this mechanism. However, incorporating a small cluster of five MeOH molecules solvating the methoxide anion, the ion pair was found to be only 11.0 kcal mol⁻¹ above 2a. Using this cluster-continuum model in DFT calculations the process in Fig. 8 takes place through a fairly accessible transition state of 20.9 kcal mol⁻¹ (TSa, Fig. 8).

Fig. 8

a 2a to 2b isomerization. b Relevant structures in the isomerization process computed with a cluster-continuum solvent model including five MeOH molecules solvating methoxide anion. In blue, relative Gibbs energies, in kcal mol⁻¹. 3D-structures were generated using CYLview [36]

The mechanistic investigation of the catalytic transfer hydrogenation of ketones to alcohols in water using a water-soluble ruthenium catalyst and aqueous HCOONa/HCOOH as the hydrogen source at pH 4.4 illustrates the usefulness of careful employing of cluster-continuum models to deal with reaction mechanisms.

In the experimental investigation, aqua, formato, and hydrido species were detected by ¹H NMR experiments in D₂O [51]. Overall, the experimental evidence is consistent with the following sequence of reactions: (i) chloride by water ligand exchange, (ii) water by formate ligand exchange, (iii) decarboxylation of the formato complex to produce a hydrido species, and (iv) hydride and proton transfer to yield the reduced product (Fig. 9a). This sequence involves several of the situations we have been describing so far: a solvent molecule becomes a ligand, small and charged species are exchanged, and a charged reagent (proton from the water medium) is added. DFT calculations, using a cluster-continuum solvent model, provided a detailed microscopic description of the full catalytic cycle. Three explicit water molecules were included in the QM description of the system in addition to the polarizable continuum model (water, ε = 78.36.) All the intermediates and transition states were located by full optimization of the complex + substrate + water cluster species into the solvent. Using this cluster-continuum model the transition states for the chloride by water, and water by formate ligand exchanges, were located. The effect of the cluster size was assessed by increasing the number of water molecules to five in selected steps (Fig. 9b). The barrier of the hydride transfer with two more explicit solvent molecules increases only 1.4 kcal mol⁻¹, thus validating the employed model. This reaction is an example of reactivity in protic solvent media. In the next section we will go further on this issue.

Fig. 9

a Sequence of reactions for the transfer hydrogenation of cyclohexanone with a ruthenium 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]

5 Reactions in Protic Solvents

According to the discussion presented in the precedent sections it is apparent that modeling catalytic processes in protic solvents, particularly alcohols and water, do demand inclusion of explicit solvent molecules. Protic solvents provide a strong hydrogen bonding environment that can influence the relative stabilities of reactants, intermediates and transition states and products. In a protic solvent, solvent molecules can coordinate, can be directly involved as a proton shuttle in proton transfer steps, can generate small charged species as hydroxide, proton or alkoxide and could compete with the catalyst as the proton source (i.e. metal complexes including N–H ligands). These effects are not adequately accounted for with a continuum model. Even more important, this type of solvent can drive the reaction towards new pathways, different from gas-phase or low polar solvents.

5.1 Alcohol Solvents

In a pioneering work, Meijer and Handgraaf applied AIMD techniques to simulate the ruthenium-catalyzed interconversion of formaldehyde and methanol, as a model for transfer hydrogenation in methanol solvent [52]. To model the solvated system a Ru complex and 40 methanol molecules were included in the simulation. From gas-phase calculations, a concerted transfer of a proton and a hydride from a ruthenium complex with a N–H ligand had been proposed (metal–ligand bifunctional mechanism, also called Noyori type mechanism, Scheme 5) [53]. AIMD calculations showed that the reaction in solution follows a different mechanism than in gas phase: the reaction is stepwise, hydride transfer happens first, and a methoxide intermediate exists in solution for a short but finite time.

Scheme 5

Ru-catalyzed hydrogen transfer via the metal–ligand bifunctional mechanism. Original proposed concerted mechanism

A similar conclusion arises from a recent computational study of hydrogenation of ketones with t-BuOK in tert-butanol (Scheme 6) [54]. By analogy with Ru-catalyzed ketone hydrogenations (Scheme 5), experimentalists proposed that ketone reduction proceeds through a six-membered cyclic transition state involving the ketone, H₂ and the alkoxide base [55]. A subsequent theoretical study of the reaction in solution was performed through CSM calculations on the optimized gas-phase structures [56], thus assuming the concerted transition state depicted in Scheme 6b. Static and dynamic DFT calculations (AIMD), including explicit solvent molecules give a completely different picture of the reaction, consisting of three steps: (i) cleavage of the H–H bond to afford potassium hydride; (ii) addition of potassium hydride across the C=O bond of the ketone; (iii) H/K exchange to regenerate the catalyst and afford the product (Scheme 6c) [54].

Scheme 6

a Base-catalyzed hydrogenation of ketones. b Concerted transition state. c Revisited mechanism from static DFT and AIMD simulations, including solvent molecules

Dehydrogenations in alcohol solvent are deeply influenced by solvent molecules, similarly to the above commented hydrogenations. An AIMD study of the mechanism of methanol dehydrogenation catalyzed by a highly active PNP pincer ruthenium complex in methanol solvent unveiled this behavior [57]. In all the previously proposed mechanisms, protonation and deprotonation of the amido ligand were considered to be key steps of the catalytic cycle but relayed on a Noyori-type concerted mechanism (Scheme 5) [58, 59]. DFT-MD simulations with explicit solvent evidenced the spontaneous protonation of the ligand nitrogen atom of PNP-ruthenium complex by a solvent molecule. Indeed, the amido moiety of the catalyst remained protonated during the complete catalytic process, thus precluding a metal–ligand bifunctional mechanism. Protonation was accompanied by the coordination of an additional ligand such as MeO⁻ or MeOH to the metal center, in this way forming the resting state of the catalyst (Scheme 7). Importantly, the nitrogen atom of the ligand was found to have a high proton affinity in polar, protic solvents, indicating that ever under very basic conditions the nitrogen atom should remain protonated [57].

Scheme 7

Proposed mechanism for methanol oxidation by a PNP pincer ruthenium complex with explicit solvent [49]

Static DFT calculations including a small number of explicit solvent molecules (cluster-continuum model) give a similar description of solvent effects in this reaction. Using a more realistic solvent model which incorporates explicit solvent molecules was critical to devise a new mechanistic view of this reaction.

Solvent can be crucial to stabilize reactive species. The Brönsted-acid catalyzed ring-opening of cyclopropanes takes place using triflic acid (HOTf) as the catalyst and hexafluoroisopropanol (HFIP) as a protic solvent [60]. The reaction is initiated by protonation of cyclopropane by HOTf. However, without the inclusion of explicit solvation effects, no protonated intermediate is found (no minimum in the potential energy surface). Introducing three explicit HFIP molecules the proton of HTOf is transferred to the C=O substituent of cyclopropane, enhancing the reactivity of cyclopropane ring (Scheme 8) [61].

Scheme 8

Relative stabilities of neutral and protonated forms of cyclopropane-based reactant, without and with incorporation of HFIP solvent molecules

The reliability of DFT calculations to elucidate reaction mechanism was challenged when Plata and Singleton reported a detailed experimental study of the Morita–Baylis–Hillman (MBH) reaction and demonstrated that previous computational studies using continuum solvent models failed to reproduce known mechanistic features of the reaction [62]. MBH reaction entails coupling of an activated alkene derivative with an aldehyde to form a new C–C bond. It is a multistep reaction in which while the reactants and products are neutral species with moderate polarity, many of the putative intermediates are either charged or zwitterionic. Moreover, it takes place in protic solvents. Taking into account the reaction’s scenario, it could be expected that a very accurate description of solvation is required to accurately reproduce by DFT calculations the main features of the reaction mechanism. Indeed, three subsequent computational studies have validated this hypothesis. Harvey, Sunoj and coworkers employed high-level energy calculations combined with cluster-continuum calculations including methanol molecules, particularly important for species with partial charges, and molecular dynamics simulations of the reacting system in a large box of explicit solvent [63]. Basdogan and Keith used a stochastic filtering procedure with a global optimization protocol, based on static QM calculations of micro-solvated clusters, to identify the low-energy geometries of the solute–solvent clusters. From the globally optimized intermediates reaction pathways were calculated using the growing string method. They found that gradually adding more solvent molecules does not improve the agreement with experiment and that a micro-solvated model with five explicit methanol molecules was the one in closer agreement with the experimental results [64]. Sure et al. have devised a simple approach to address the main issues facing micro-solvation: when, where, and how many explicit solvent molecules have to be added to a solute, and they have successfully applied this protocol to the MBH reaction [65]. In this approach, that only employs static DFT calculations, the proposed criterion to apply micro-solvation is the thermodynamic preference for the association of a solvent (ΔG_a < 0); the interaction sites are found from screening charge densities. The three computational studies of the MBH reaction [63,64,65] have similar accuracies compared with the experimental results (absolute errors below 4 kcal mol⁻¹ and MAD about 2 kcal mol⁻¹) and demonstrate that QM calculations can give a good description of the MHB mechanism, provided that an accurate solvent modeling is employed.

5.2 Water

Water medium can deeply affect the activation of strong non-polar bonds, as H–H and C–H bonds. Computational studies with explicit solvent molecules have disclosed active participation of solvent molecules in both processes, resulting in reaction mechanisms different from those in non-polar solvents. Using a cluster-continuum modeling of the reaction media we unraveled the role of water in the activation of dihydrogen by [Cp’Ru(PTA)₂Cl] (Cp’ = C₅H₅, C₅Me₅; PTA = 1,3,5-triaza-7-phosphaadamantane) water soluble complexes [66, 67]. Water comes into play in all the steps of the reaction. The first step leads to the substitution of an anionic Cl⁻ ligand by the weaker neutral H₂ ligand. Initial dissociation of Cl⁻ does not yield an unsaturated [Cp’Ru(PTA)₂]⁺ complex, but a monoaquo [Cp’Ru(PTA)₂(H₂O)]⁺ complex, detected by NMR. Subsequent displacement of the coordinated water by H₂, results in the complex {[Cp’Ru(PTA)₂(η²-H₂)]⁺. Thus, solvent contributes to this step by coordinating to the metal as well as stabilizing the released Cl⁻ anion. In the next step, H–H bond breaking does not happen by the common oxidative addition mechanism, but by an heterolytic splitting pathway in which the water solvent takes an active part, acting as an external base by deprotonating the dihydrogen ligand. Even in the presence of the more basic C₅Me₅ ligand the heterolytic splitting is favored over the oxidative addition [67]. The proton ends up either at a nitrogen atom of the PTA ligand, yielding a monohydride with a protonated PTA ligand (for the complex with Cp’ = C₅H₅) or at the metal, affording a dihydride (for the complex with Cp’ = C₅Me₅). Interestingly, when a small cluster of three water molecules was used in the cluster-continuum solvent model, only a concerted water-chain assisted deprotonation of the H₂ ligand and proton transfer to one N atom of the PTA ligand was obtained (Fig. 10a) [66]. Introduction of an additional water molecule in the cluster allows the stabilization of a hydronium ion in the water medium than can protonate either of the basic centers (Fig. 10b) [67]. It seems that the dichotomy between proton-shuttle and acid–base mechanisms [62] is highly dependent on the size of the solvent cluster in micro-solvated models.

Fig. 10

Heterolytic splitting of the H–H bond of [CpRu(PTA)₂(η²-H₂)]⁺. a Cluster of three water molecules: transition state of a concerted water-chain-assisted deprotonation of the H₂ ligand and proton transfer to the PTA ligand (proton shuttle mechanism) b Cluster of four water molecules: monohydride [Cp*Ru(PTA)₂H] and solvated hydronium ion intermediate (acid–base mechanism). 3D-structures were generated using CYLview [36]

The C–H bond activation of methane in Shilov’s system over cis- and trans-PtCl₂(H₂O)(CH₄) complexes was investigated by ab initio molecular dynamics in explicit water [68]. Simulations within a realistic solvent representation revealed that the square pyramidal hydride complex PtHCl₂(H₂O)(CH₃) formed in the oxidative addition step (Scheme 9b) is not stable and spontaneously evolves into the square planar methyl complex by the release of the proton to the bulk solution (Scheme 9). Compared with previous modeling based on minimal solvent clusters [69, 70], the explicit solvent modeling revealed that the platinum hydride is a transient species with strong acid character, which spontaneously released a proton to the solution rather than coordinating a further water molecule to stabilize the metal center [68].

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

C–H bond activation is a key step in the Ru-catalyzed dehydrogenation of aqueous methanol. In this reaction, in which three equivalents of H₂ and CO₂ are produced from CH₃OH (Scheme 10a), C–H activation steps take place for methanol, methanediol and formic acid. Cluster-continuum static DFT and DFT-based molecular dynamics (DFT-MD) calculations have been performed to provide mechanistic insights into this reaction [71]. The active form of the catalyst (complex 2’) binds methanol and undergoes internal proton transfer from the OH group of methanol to the amido moiety to form complex 2’A, that rearranges to form the zwitterionic agostic methoxide adduct 2’A’ (Scheme 10b). C-H bond breaking occurs from this intermediate.

Scheme 10

a Elementary steps in the dehydrogenation of aqueous methanol to three equivalents of H₂ and CO_2. b Relevant intermediates. c Micro-solvated transition state for the methanol C-H bond activation

Constrained DFT-MD and static DFT micro-solvation calculations using explicit solvent molecules yield very similar results for this process: the solvent actively participates in the reaction, assisting the cleavage of the C–H bond by means of hydrogen bonding interactions (Scheme 10c). In fact, an entirely uphill energy profile is obtained without the presence of explicit solvent molecules. Thus, the inclusion of explicit solvent is essential for the right modeling of the reaction.

The importance of hydrogen bond interactions between the reacting system and the solvent has been also highlighted in the AIMD study of Ru-catalyzed aqueous transfer hydrogenation using formate as hydrogen donor. Explicit aqueous solvent was included in the simulations (54 water molecules). Two key steps were analyzed: the hydride transfer between formate and the catalyst and the dissociation of the ruthenium-formato complex [72]. Compared with a gas phase model, explicit aqueous solvent increases the barrier for the hydride transfer due to the stabilization of the reactant complex by several hydrogen bonds between formate and solvent molecules. This stabilization is absent in the product because CO₂ does not form hydrogen bonds. Conversely, the explicit solvent decreases the barrier for the formate dissociation because favors the formate to be in solution. Therefore, aqueous solvent provides a significant contribution to the reaction barriers. Moreover, in water the catalyst remains protonated, and a water molecule transfers a proton instead of the catalyst to accomplish the hydrogenation [73].

A combined experimental and computational study of the ruthenium-catalyzed isomerization of allylic alcohols in aqueous medium (Fig. 11) disclosed another possible role that water solvent molecules can play: they not only can coordinate but can act as cooperating ligand in metal–ligand bifunctional mechanisms [74].

Fig. 11

a Ru-catalyzed isomerization of allylic alcohols in water. b Bis-allyl ruthenium(IV) catalysts. c Transition states for the hydrogen transfer steps. 3D-structures were generated using CYLview [36]

During the mechanistic analysis cationic chloride free and aqua species were detected in the ESI mass spectrum. Calculations evidenced the accessibility of hydroxo-imidazole species in aqueous solution that can be the active species in the isomerization process. The isomerization takes place through two sequential concerted hydrogen transfer steps, from the substrate to the catalyst (dehydrogenative step: TS1, Fig. 11) and subsequent transfer back to the substrate (acrolein hydrogenation: TS2, Fig. 11). Both steps disclose an outer-sphere bifunctional catalyzed mechanism, with a cooperative action of the metal and the hydroxo ligand (Fig. 11c). Likewise to the previous dehydrogenation described processes, this study also revealed that in water the participation in the bifunctional catalysis of the N–H group of the azole ligand is not required, provided that a water ligand could play this role. Indeed, the related complex including a N-methylimidazole ligand, where both nitrogen atoms of the ring are substituted, shows a remarkable catalytic activity in water [74].

6 Solvent as a Reagent

In the previous section we have underlined the difficulties entailing modeling of reactions in protic solvents, and particularly water. This is even more challenging when water is directly involved in the process as a reagent. The Wacker process and O–O bond formation through water nucleophilic attack are both prototypical examples of reactions in which water acts as a solvent and as a reagent.

The Wacker process is the aerobic oxidation of ethylene to acetaldehyde catalyzed by palladium(II) chloride in aqueous solution. The net reaction can be described as follows:

$$\left[ {{\text{PdCl}}_{{4}} } \right]^{{{2 } - }} + {\text{C}}_{{2}} {\text{H}}_{{4}} + {\text{H}}_{{2}} {\text{O}} \to {\text{CH}}_{{3}} {\text{CHO}} + {\text{Pd}} + {\text{2 HCl}} + {\text{2 Cl}}^{ - }$$

The Pd(II) catalyst is regenerated by oxidation of Pd(0) by copper(II) reagent (CuCl₂). In this reaction the water solvent is initially involved in the ligand exchange that leads to the catalytically active species, and then in the nucleophilic attack to the coordinated double bond (hydroxypalladation). During these processes neighboring solvent molecules will influence the reactivity of a water molecule, making mandatory explicit solvent modeling. A proper description of the twofold role of water, and thus of the reaction mechanism, can hardly be achieved by using static DFT calculations with a continuum or cluster-continuum solvent models. This has 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 nucleophilic attack by a water molecule from the bulk, and (b) syn nucleophilic attack by a coordinated hydroxide group located cis to the alkene (Scheme 11). AIMD simulations, combined with rare event sampling techniques (as metadynamics) were the approaches at choice to evaluate a long-standing controversy regarding the Wacker’s mechanism [75].

Scheme 11

Postulated reaction pathways for the Wacker process: syn and anti addition of water to the Pd-coordinated ethene

One of the questions revealed by the simulations was the nature of the active species formed by a sequence of ligand exchanges from the catalyst precursor [PdCl₄]^{2 –} (1). The ligand exchanges take place in a stepwise fashion for the formation of both cis and trans isomers, complexes 5 and 3, respectively (Scheme 12). Simulations show that the formation of cis compound is uphill in Gibbs energy by about 12 kcal/mol, whereas the formation of trans intermediate is nearly thermoneutral on the Gibbs energy scale, by ~ 1 kcal/mol [76, 77]. The Gibbs energy profiles obtained for the formation of both isomers are kinetically feasible. However, the formation of the cis isomer is clearly thermodynamically unfavorable. Indeed, the energy difference between both isomers 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.

Scheme 12

Ligand exchanges on the catalyst precursor [PdCl₄]²⁻ leading to the trans- and cis-[PdCl₂(C₂H₄)(H₂O)] species

According to the results obtained for the speciation of the [PdCl₄]²⁻ 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 [76, 79]. Simulations performed under the first conditions, low Cl⁻ concentration, were carried out on both cis- and trans-[PdCl₂(C₂H₄)(H₂O)] complexes. For the trans isomer the obtained barrier for the hydroxypalladation step was compatible with experimental observations: the calculated barrier of 19 kcal mol⁻¹ is slightly lower than the barrier obtained experimentally, 22 kcal mol⁻¹. Simulations starting from the cis-[PdCl₂(C₂H₄)(H₂O)] complex showed that isomerization to the trans complex and a subsequent outer sphere nucleophilic addition is more feasible than the inner sphere mechanism. Moreover, metadynamics calculations were also performed by allowing the inner and outer sphere mechanisms take place simultaneously; simulations revealed that the outer sphere mechanism is again the preferred one. Simulations at high Cl⁻ concentration generated similar conclusions regarding the inner/outer reaction mechanism comparison.

Another appealing example of reaction in which water is the reagent and the solvent is the O–O bond formation in the catalytic oxidation of water through the called water nucleophilic attack mechanism (WNA). In the WNA mechanism the O–O bond is formed by a water attack on an oxo ligand (Scheme 13) [80].

Scheme 13

O–O bond formation pathway via water nucleophilic attack (WNA). M is a transition metal complex in high oxidation state

A series of recent papers, dealing with different catalysts, point out the suitability of AIMD simulations with explicit solvent to tackle this challenging process [81,82,83].

7 Summary and Conclusions

The use of continuum solvent models to study chemical processes is a great success in the use of theoretical approaches to model the chemistry performed in the lab. However, they cannot be applied to whatever situation. In many cases introduction of explicit solvent molecules in the QM description of the system is mandatory. It is not in our aim, and we think it is not possible, to give a simple recipe to advise/propose when explicit solvent molecules must be incorporated in the QM description of the system, along with the continuum method or as a fully explicit model. Addition of explicit solvent molecules always make calculations more complex. Moreover, both cluster-continuum and fully explicit solvent models have inherent limitations commented in the Introduction section. Thus, if a good description of the reaction 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 situations 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 illustrate the reasons that urge to go beyond the continuum models to the employment of micro-solvated of fully explicit solvent models. This section summarizes the examples analyzed.

The success and limitations of CSM are clearly shown by means of three different examples. One example performs a comparative analysis over non-catalyzed and Ag-catalyzed ring-closing reaction steps. The work shows that pairwise interactions between the solute and the polar solvent molecules are significant, but they do not induce any kind of ordering in the solvation environment, therefore, the process can be safely described with a continuum model that provides “averaged” solvent effects. A second example analyses the decomposition of peroxyformate in many different solvents. The main transition state, described as a proton abstraction process, clearly shows that when the solvent influence is due to its polarity, continuum models can capture solvent effects, showing the success of CSM. The third example, devoted to a metal catalyzed dehydrogenation process, shows that inclusion of solvent effects (by means of continuum and explicit solvent molecules) not only changes the reaction mechanism from concerted to stepwise, but also modifies the overall view of its mechanism. This example demonstrates the importance of properly accounting for solvent effects in mechanistic studies.

Another important aspect to pay attention when deciding the solvent model is the solvent capacity to coordinate reagents (especially when using metal catalysts); this is commonly stated as solvent dynamics. Indeed, the association-dissociation process can affect the speciation of reactive species and they can modify the course of reaction. This is shown through several examples as the speciation of Negishi coupling reagents or the Schlenk equilibrium of Grignard reagents. In the former, the study concluded that species ZnMe₂(THF)₂, ZnMeCl(THF)_2, and ZnCl₂(THF)₂ should be the ones considered in any computational modeling of these zinc(II) species in THF solvent. The latter shows that the reaction is controlled by solvent dynamics: asymmetric solvation on Mg ions is required to promote the Mg–Cl and Mg–CH₃ bond cleavage. The general advise is that when good coordinating solvents are employed, they must be included in the model. For poorly coordinating solvents the answer is not so clear, though the described example shows that “empty coordination sites” are not accessible even for non-coordinating solvents.

Another critical issue is related to the solvation of small and charged species. As the described examples highlight, inclusion of explicit solvent molecules is required for a realistic modeling of processes involving such species. In one example, a protodeauration reaction step, computations shows that solvation modifies the electrophilicity of the protonating agent: it is overestimated by a pure continuum model. The electrophilicity of the proton decreases by increasing the number of solvent water molecules in the solvation sphere; this is reflected by the increase in the activation barrier for the protodeauration. Therefore, inclusion of explicit solvent molecules is mandatory in these cases. In another example, related to the behavior of ion pairs, it is shown that their description can be also highly affected by the inclusion or not of explicit solvent molecules in the model. The example selected shows that the presence of explicit solvent molecules modifies the relative energy of the formed ion pair during the reaction pathway. Neglecting the incorporation of explicit solvent molecules could drive to a wrong conclusion about the mechanism.

One of the paradigmatic cases where continuum models need to be complemented are probably the reactions performed in protic solvents. Protic solvents provide a strong hydrogen bonding environment. In reaction profiles this can largely influence the relative stabilities of reactants, intermediates, transition states and products. Moreover, protic solvents can also coordinate, or they can be directly involved as proton shuttles in proton transfer steps. Note that protic solvents can also generate small charged species as hydroxide, proton or alkoxide, thus becoming small and charged species (earlier commented). Their behavior is generally not properly described by continuum models alone, therefore becoming prominent candidates to employ cluster-continuum models. In addition, this type of solvents can drive the reaction towards new pathways, different from gas-phase or low polar solvents.

The transfer hydrogenation reaction in alcohol solvent was described as a paradigmatic case. A concerted transfer of a proton and a hydride from the metal complex with a N–H ligand was initially proposed (metal–ligand bifunctional mechanism, also called Noyori type mechanism). AIMD simulations disclosed that the reaction pathway in solution was operating through a different mechanism than in gas phase (even including CSM). The reaction is stepwise, hydride transfer happens first, and then the proton transfer; moreover, the latter may not involve even the metal catalyst. Several examples on how water solvent affect the activation of H–H or C–H are also commented. In all cases, the presence of explicit solvent molecules is mandatory to obtain reliable results on the reaction mechanism.

Intimately close to the previous examples are the reactions where the solvent is a reagent in the process. The paradigmatic cases are those where the reaction takes place in water solvent, and water molecules are also reactants in the process. The Wacker process and the water nucleophilic addition reaction were selected. In both cases the use of AIMD procedures was employed (therefore using a fully explicit solvent model) to gain a clearer view of the overall chemical process.

Overall, the selected examples showed the successes and failures of the continuum model, and which cases the hybrid cluster-continuum model or even fully explicit models need to be considered. These models represent significant improvements in the computational modeling of catalytic systems but have inherent limitations. Certainly, proper description of solvation is a difficult and tricky issue that need to be cautiously considered for the accurate computational modeling of chemistry in solution.